diff --git a/src/math/interval/dep_intervals.h b/src/math/interval/dep_intervals.h index f5851c1f0..c928f6426 100644 --- a/src/math/interval/dep_intervals.h +++ b/src/math/interval/dep_intervals.h @@ -210,6 +210,19 @@ public: } } + template + void neg(const interval &a, interval &b) { + if (with_deps == wd) { + interval_deps_combine_rule combine_rule; + m_imanager.neg(a, b, combine_rule); + combine_deps(a, combine_rule, b); + } + else { + m_imanager.neg(a, b); + } + TRACE(dep_intervals, tout << "neg of "; display(tout, a) << " = "; display(tout, b) << "\n";); + } + template void power(const interval& a, unsigned n, interval& b) { if (n == 1 && &a == &b) diff --git a/src/math/lp/nla_stellensatz2.cpp b/src/math/lp/nla_stellensatz2.cpp index 24c165b4f..6babcef99 100644 --- a/src/math/lp/nla_stellensatz2.cpp +++ b/src/math/lp/nla_stellensatz2.cpp @@ -4,23 +4,6 @@ vanishing(v, p) = (q, r) with q = 0, r >= 0 if p := q*v + r, M(q) = 0 - TODOs: - - Assumptions are added to impose sign constraints on non-linear multipliers - or to identify vanishing polynomials. - The sign of the assumption is based on the current candidate interpretation. - It could be that the negation of this assumption is asserted (and is - false under the current candiate interpretation). - The procedure could loop around refining a false assumption. - To handle this we could also track bounds of polynomials. - Then we can check if a bound on a polynomial is already asserted, - either as a case split or as an input assertion. - To retract consequences of temporary assumed bounds we would associate - levels with constraints and also undo and replay constraints during backjumping/backtracking. - A bound assumption corresponds to a decision. A conflict that depends on a bound assumption - reverts the assumptiom, and justifies it by the conflict dependencies. - In this way we would eliminate the double nested loop: saturate looping over the search loop. - --*/ #include @@ -306,6 +289,7 @@ namespace nla { m_justifications.push_back(j); m_constraint_index.insert({p.index(), k}, ci); push_bound(ci); + push_constraint(ci); ++c().lp_settings().stats().m_stellensatz.m_num_constraints; m_has_occurs.reserve(ci + 1); return ci; @@ -885,20 +869,25 @@ namespace nla { bool stellensatz2::constraint_is_bound_conflict(lp::constraint_index ci) { auto const &c = m_constraints[ci]; - auto const& [p, k] = c; + auto const &[p, k] = c; scoped_dep_interval iv(m_di); - interval(p, iv); - if (iv->m_upper_inf) + calculate_interval(iv, p); + return constraint_is_bound_conflict(ci, iv); + } + + bool stellensatz2::constraint_is_bound_conflict(lp::constraint_index ci, dep_interval const &iv) { + auto const &[p, k] = m_constraints[ci]; + if (iv.m_upper_inf) return false; - if (iv->m_upper > 0) + if (iv.m_upper > 0) return false; - bool is_negative = iv->m_upper < 0 || iv->m_upper_open; - SASSERT(is_negative || iv->m_upper == 0); + bool is_negative = iv.m_upper < 0 || iv.m_upper_open; + SASSERT(is_negative || iv.m_upper == 0); bool is_conflict = k == lp::lconstraint_kind::GT || (k == lp::lconstraint_kind::GE && is_negative); if (!is_conflict) return false; m_conflict_dep.reset(); - m_dm.linearize(iv->m_lower_dep, m_conflict_dep); + m_dm.linearize(iv.m_lower_dep, m_conflict_dep); TRACE(arith, tout << "constraint is bound conflict: "; display_constraint(tout, ci) << "\n";); return true; } @@ -923,20 +912,27 @@ namespace nla { // // Based on current bounds for variables, compute bounds of polynomial // - void stellensatz2::interval(dd::pdd p, scoped_dep_interval& iv) { - auto &m = iv.m(); + void stellensatz2::calculate_interval(scoped_dep_interval& out, dd::pdd p) { + auto &m = out.m(); if (p.is_val()) { - m.set_value(iv, p.val()); + m.set_value(out, p.val()); return; } - scoped_dep_interval x(m), lo(m), hi(m); - auto v = p.var(); + scoped_dep_interval lo(m), hi(m), x(m); + calculate_interval(x, p.var()); + calculate_interval(hi, p.hi()); + calculate_interval(lo, p.lo()); + calculate_interval(out, x, lo, hi); + } + + void stellensatz2::calculate_interval(scoped_dep_interval& x, lpvar v) { + auto &m = x.m(); auto lb = get_lower(v); if (lb == UINT_MAX) { m.set_lower_is_inf(x, true); } else { - auto const& lo = m_bounds[lb]; + auto const &lo = m_bounds[lb]; m.set_lower_is_inf(x, false); m.set_lower_is_open(x, lo.m_kind == lp::lconstraint_kind::GT); m.set_lower(x, lo.m_value); @@ -947,16 +943,19 @@ namespace nla { m.set_upper_is_inf(x, true); } else { - auto const& hi = m_bounds[ub]; + auto const &hi = m_bounds[ub]; m.set_upper_is_inf(x, false); m.set_upper_is_open(x, hi.m_kind == lp::lconstraint_kind::LT); m.set_upper(x, hi.m_value); m.set_upper_dep(x, hi.d); } - interval(p.hi(), hi); - interval(p.lo(), lo); - m.mul(x, hi, hi); - m.add(lo, hi, iv); + } + + void stellensatz2::calculate_interval(scoped_dep_interval& out, dep_interval const& x, dep_interval const& lo, dep_interval const& hi) { + auto &m = out.m(); + scoped_dep_interval tmp(m); + m.mul(x, hi, tmp); + m.add(lo, tmp, out); } // @@ -1148,14 +1147,17 @@ namespace nla { m_justifications.shrink(head); m_levels.shrink(head); // re-insert bounds - for (; m_bounds.size() >= ci;) + for (; m_bounds.size() >= ci;) { pop_bound(); - for (; m_bounds.size() < head;) + pop_constraint(m_bounds.size()); + } + for (; m_bounds.size() < head;) { push_bound(m_bounds.size()); + push_constraint(m_bounds.size() - 1); + } SASSERT(well_formed_last_bound()); - reset_conflict(); - m_prop_qhead = ci; + reset_conflict(); } stellensatz2::justification stellensatz2::translate_j(std::function const& f, @@ -1168,6 +1170,12 @@ namespace nla { cs.push_back(f(ci)); return bound_propagation_justification{f(arg.ci), cs}; } + else if constexpr (std::is_same_v) { + svector cs; + for (auto ci : arg.cs) + cs.push_back(f(ci)); + return propagation_justification{cs}; + } else if constexpr (std::is_same_v) { return assumption_justification{}; } @@ -1218,6 +1226,10 @@ namespace nla { auto const &bpj = std::get(j); return index == 0 ? bpj.ci : bpj.cs[index - 1]; } + if (std::holds_alternative(j)) { + auto const &bpj = std::get(j); + return bpj.cs[index]; + } if (std::holds_alternative(j)) { UNREACHABLE(); } @@ -1264,6 +1276,8 @@ namespace nla { unsigned stellensatz2::num_constraints(justification const &j) { if (std::holds_alternative(j)) return 1 + std::get(j).cs.size(); + if (std::holds_alternative(j)) + return std::get(j).cs.size(); if (std::holds_alternative(j)) return 0; if (std::holds_alternative(j)) @@ -1286,57 +1300,25 @@ namespace nla { return 0; } - void stellensatz2::propagate() { - if (m_prop_qhead == m_constraints.size()) + void stellensatz2::propagate_constraint(lpvar x, lp::lconstraint_kind k, rational const& value, + lp::constraint_index ci, svector& cs) { + TRACE(arith, display_constraint(tout << "constraint is propagating ", ci) << "\n"; + tout << "v" << x << " " << k << " " << value << "\n";); + + // block repeated bounds propagation + if (propagation_cycle(x, value, k, ci, cs)) return; - - for (; m_prop_qhead < m_constraints.size(); ++m_prop_qhead) { - if (!c().reslim().inc()) - return; - auto q = m_bounds[m_prop_qhead].q; - if (is_bound_conflict(q)) - return; - - if (!q.is_var()) - continue; - auto v = q.var(); - for (unsigned i = 0; i < m_occurs[v].size(); ++i) { - lp::constraint_index ci = m_occurs[v][i]; + ci = add_constraint(pddm.mk_var(x) - value, k, bound_propagation_justification{ci, cs}); - TRACE(arith, tout << "propagate v" << v << " (" << ci << ")\n"); + if (constraint_is_bound_conflict(ci)) + return; - // detect conflict or propagate dep_intervals + set_in_bounds(x); - if (constraint_is_bound_conflict(ci)) - return; - - lpvar w; - rational value; - lp::lconstraint_kind k; - svector cs; - if (!constraint_is_propagating(ci, cs, w, k, value)) - continue; - - TRACE(arith, display_constraint(tout << "constraint is propagating ", ci) << "\n"; - tout << "v" << w << " " << k << " " << value << "\n";); - - // block repeated bounds propagation - if (propagation_cycle(w, value, k, ci, cs)) - continue; - - ci = add_constraint(pddm.mk_var(w) - value, k, bound_propagation_justification{ci, cs}); - - if (constraint_is_bound_conflict(ci)) - return; - - set_in_bounds(w); - - CTRACE(arith, !well_formed_last_bound(), display(tout)); - SASSERT(well_formed_last_bound()); - SASSERT(well_formed_var(w)); - } - } + CTRACE(arith, !well_formed_last_bound(), display(tout)); + SASSERT(well_formed_last_bound()); + SASSERT(well_formed_var(x)); } bool stellensatz2::propagation_cycle(lpvar v, rational const &value, lp::lconstraint_kind k, lp::constraint_index ci, svector& cs) const { @@ -1427,8 +1409,7 @@ namespace nla { unsigned stellensatz2::max_level(constraint const &c) const { unsigned level = 0; - auto const &vars = c.p.free_vars(); - for (auto v : vars) + for (auto v : c.p.free_vars()) level = std::max(level, m_var2level[v]); return level; } @@ -1449,99 +1430,108 @@ namespace nla { lp::lconstraint_kind &k, rational &value) { auto [p, ck] = m_constraints[ci]; cs.reset(); + k = ck; for (auto x : p.free_vars()) { auto f = factor(x, ci); if (f.degree > 1) continue; scoped_dep_interval ivp(m_di), ivq(m_di); - interval(f.p, ivp); - interval(-f.q, ivq); + calculate_interval(ivp, f.p); + calculate_interval(ivq, -f.q); + TRACE(arith_verbose, tout << "variable v" << x << " in " << p << "\n"; m_di.display(tout << "interval: " << f.p << ": ", ivp) << "\n"; m_di.display(tout << "interval: " << -f.q << ": ", ivq) << "\n"; display_constraint(tout, ci) << "\n"); - - // hi_p > 0 - // 0 <= lo_q <= -q - // => x >= lo_q / hi_p - // - if (!ivp->m_upper_inf && ivp->m_upper > 0 && !ivq->m_lower_inf && ivq->m_lower >= 0) { - // v >= -q / p - auto quot = ivq->m_lower / ivp->m_upper; - if (var_is_int(x)) - quot = ceil(quot); - bool is_strict = !var_is_int(x) && (ck == lp::lconstraint_kind::GT || ivq->m_lower_open || ivp->m_lower_open); - if (!has_lo(x) || quot > lo_val(x) || (!lo_is_strict(x) && quot == lo_val(x) && is_strict)) { - TRACE(arith, tout << "new lower bound v" << x << " " << quot << "\n"); - auto d = m_dm.mk_join(ivq->m_lower_dep, ivp->m_upper_dep); - m_dm.linearize(d, cs); - k = is_strict ? lp::lconstraint_kind::GT : lp::lconstraint_kind::GE; - v = x; - value = quot; - return true; - } - } - // hi_p >= p >= lo_p > 0 - // 0 > hi_q >= -q >= lo_q - // => x >= lo_q / lo_p - if (!ivp->m_lower_inf && ivp->m_lower > 0 && !ivq->m_lower_inf && ivq->m_lower <= 0) { - auto quot = ivq->m_lower / ivp->m_lower; - if (var_is_int(x)) - quot = ceil(quot); - bool is_strict = - !var_is_int(x) && (ck == lp::lconstraint_kind::GT || ivq->m_lower_open || ivp->m_lower_open); - if (!has_lo(x) || quot > lo_val(x) || (!lo_is_strict(x) && quot == lo_val(x) && is_strict)) { - TRACE(arith, tout << "new lower bound v" << x << " " << quot << "\n"); - auto d = m_dm.mk_join(ivp->m_lower_dep, ivq->m_lower_dep); - m_dm.linearize(d, cs); - k = is_strict ? lp::lconstraint_kind::GT : lp::lconstraint_kind::GE; - v = x; - value = quot; - return true; - } + v = x; + k = ck; + if (constraint_is_propagating(ivp, ivq, v, cs, k, value)) + return true; + } + + return false; + } + + bool stellensatz2::constraint_is_propagating(dep_interval const &ivp, dep_interval const &ivq, lpvar x, + svector &cs, lp::lconstraint_kind &k, + rational &value) { + // hi_p > 0 + // 0 <= lo_q <= -q + // => x >= lo_q / hi_p + // + if (!ivp.m_upper_inf && ivp.m_upper > 0 && !ivq.m_lower_inf && ivq.m_lower >= 0) { + // v >= -q / p + auto quot = ivq.m_lower / ivp.m_upper; + if (var_is_int(x)) + quot = ceil(quot); + bool is_strict = + !var_is_int(x) && (k == lp::lconstraint_kind::GT || ivq.m_lower_open || ivp.m_lower_open); + if (!has_lo(x) || quot > lo_val(x) || (!lo_is_strict(x) && quot == lo_val(x) && is_strict)) { + TRACE(arith, tout << "new lower bound v" << x << " " << quot << "\n"); + auto d = m_dm.mk_join(ivq.m_lower_dep, ivp.m_upper_dep); + m_dm.linearize(d, cs); + k = is_strict ? lp::lconstraint_kind::GT : lp::lconstraint_kind::GE; + value = quot; + return true; } - - // p <= hi_p < 0 - // lo_q <= q <= hi_q < 0 - // => x <= lo_q / hi_p - if (!ivp->m_upper_inf && ivp->m_upper < 0 && !ivq->m_upper_inf && !ivq->m_lower_inf && ivq->m_upper <= 0) { - // v <= -q / p - auto quot = ivq->m_lower / ivp->m_upper; - if (var_is_int(x)) - quot = floor(quot); - bool is_strict = - !var_is_int(x) && (ck == lp::lconstraint_kind::GT || ivq->m_lower_open || ivp->m_upper_open); - if (!has_hi(x) || quot < hi_val(x) || (!hi_is_strict(x) && quot == hi_val(x) && is_strict)) { - TRACE(arith, tout << "new upper bound v" << x << " " << quot << "\n"); - auto d = m_dm.mk_join(ivq->m_upper_dep, m_dm.mk_join(ivq->m_lower_dep, ivp->m_upper_dep)); - m_dm.linearize(d, cs); - k = is_strict ? lp::lconstraint_kind::LT : lp::lconstraint_kind::LE; - v = x; - value = quot; - return true; - } + } + + // hi_p >= p >= lo_p > 0 + // 0 > hi_q >= -q >= lo_q + // => x >= lo_q / lo_p + if (!ivp.m_lower_inf && ivp.m_lower > 0 && !ivq.m_lower_inf && ivq.m_lower <= 0) { + auto quot = ivq.m_lower / ivp.m_lower; + if (var_is_int(x)) + quot = ceil(quot); + bool is_strict = + !var_is_int(x) && (k == lp::lconstraint_kind::GT || ivq.m_lower_open || ivp.m_lower_open); + if (!has_lo(x) || quot > lo_val(x) || (!lo_is_strict(x) && quot == lo_val(x) && is_strict)) { + TRACE(arith, tout << "new lower bound v" << x << " " << quot << "\n"); + auto d = m_dm.mk_join(ivp.m_lower_dep, ivq.m_lower_dep); + m_dm.linearize(d, cs); + k = is_strict ? lp::lconstraint_kind::GT : lp::lconstraint_kind::GE; + value = quot; + return true; } - // lo_p <= p < 0 - // 0 <= lo_q <= -q - // => x <= lo_q / lo_p - // - if (!ivp->m_upper_inf && ivp->m_upper < 0 && !ivp->m_lower_inf && - !ivq->m_lower_inf && ivq->m_lower >= 0) { - auto quot = ivq->m_lower / ivp->m_lower; - if (var_is_int(x)) - quot = floor(quot); - bool is_strict = - !var_is_int(x) && (ck == lp::lconstraint_kind::GT || ivq->m_lower_open || ivp->m_lower_open); - if (!has_hi(x) || quot < hi_val(x) || (!hi_is_strict(x) && quot == hi_val(x) && is_strict)) { - TRACE(arith, tout << "new upper bound v" << x << " " << quot << "\n"); - auto d = m_dm.mk_join(ivp->m_upper_dep, m_dm.mk_join(ivp->m_lower_dep, ivq->m_lower_dep)); - m_dm.linearize(d, cs); - k = is_strict ? lp::lconstraint_kind::LT : lp::lconstraint_kind::LE; - v = x; - value = quot; - return true; - } + } + + // p <= hi_p < 0 + // lo_q <= q <= hi_q < 0 + // => x <= lo_q / hi_p + if (!ivp.m_upper_inf && ivp.m_upper < 0 && !ivq.m_upper_inf && !ivq.m_lower_inf && ivq.m_upper <= 0) { + // v <= -q / p + auto quot = ivq.m_lower / ivp.m_upper; + if (var_is_int(x)) + quot = floor(quot); + bool is_strict = + !var_is_int(x) && (k == lp::lconstraint_kind::GT || ivq.m_lower_open || ivp.m_upper_open); + if (!has_hi(x) || quot < hi_val(x) || (!hi_is_strict(x) && quot == hi_val(x) && is_strict)) { + TRACE(arith, tout << "new upper bound v" << x << " " << quot << "\n"); + auto d = m_dm.mk_join(ivq.m_upper_dep, m_dm.mk_join(ivq.m_lower_dep, ivp.m_upper_dep)); + m_dm.linearize(d, cs); + k = is_strict ? lp::lconstraint_kind::LT : lp::lconstraint_kind::LE; + value = quot; + return true; + } + } + // lo_p <= p < 0 + // 0 <= lo_q <= -q + // => x <= lo_q / lo_p + // + if (!ivp.m_upper_inf && ivp.m_upper < 0 && !ivp.m_lower_inf && !ivq.m_lower_inf && ivq.m_lower >= 0) { + auto quot = ivq.m_lower / ivp.m_lower; + if (var_is_int(x)) + quot = floor(quot); + bool is_strict = + !var_is_int(x) && (k == lp::lconstraint_kind::GT || ivq.m_lower_open || ivp.m_lower_open); + if (!has_hi(x) || quot < hi_val(x) || (!hi_is_strict(x) && quot == hi_val(x) && is_strict)) { + TRACE(arith, tout << "new upper bound v" << x << " " << quot << "\n"); + auto d = m_dm.mk_join(ivp.m_upper_dep, m_dm.mk_join(ivp.m_lower_dep, ivq.m_lower_dep)); + m_dm.linearize(d, cs); + k = is_strict ? lp::lconstraint_kind::LT : lp::lconstraint_kind::LE; + value = quot; + return true; } } return false; @@ -1913,8 +1903,7 @@ namespace nla { m_var2level[w] = l; m_level2var[l0] = v; m_level2var[l] = w; - } - + } std::ostream &stellensatz2::display(std::ostream &out) const { for (unsigned v = 0; v < num_vars(); ++v) @@ -1924,6 +1913,58 @@ namespace nla { display_constraint(out, ci) << "\n"; display(out << "\t<- ", m_justifications[ci]) << "\n"; } + + // Display propagation data structures + out << "\n=== Propagation State ===\n"; + + // Display polynomial queue + out << "Polynomial queue (qhead=" << m_prop_qhead << "):\n"; + for (unsigned i = 0; i < m_polynomial_queue.size(); ++i) { + out << " [" << i << "]" << (i < m_prop_qhead ? " (processed)" : "") << " "; + m_polynomial_queue[i].display(out) << "\n"; + } + + // Display intervals + out << "\nIntervals:\n"; + for (unsigned idx = 0; idx < m_intervals.size(); ++idx) { + auto const& ivs = m_intervals[idx]; + if (ivs.empty()) + continue; + auto const& iv = *ivs.back(); + out << " [" << idx << "] "; + m_di.display(out, iv) << "\n"; + } + + // Display parents + out << "\nParents:\n"; + for (unsigned idx = 0; idx < m_parents.size(); ++idx) { + auto const& parents = m_parents[idx]; + if (parents.empty()) + continue; + out << " [" << idx << "] -> { "; + bool first = true; + for (auto const& p : parents) { + if (!first) out << ", "; + first = false; + p.display(out); + } + out << " }\n"; + } + + // Display factors + out << "\nFactors:\n"; + for (unsigned idx = 0; idx < m_factors.size(); ++idx) { + auto const& factors = m_factors[idx]; + if (factors.empty()) + continue; + out << " [" << idx << "]:\n"; + for (auto const& [x, f, ci] : factors) { + out << " var=j" << x << " deg=" << f.degree << " p="; + f.p.display(out) << " q="; + f.q.display(out) << " ci=" << ci << "\n"; + } + } + return out; } @@ -2138,4 +2179,205 @@ namespace nla { } } } -} + + // incremental bounds propagation + // When adding a new bound x >= r, then + // set up to walk all parents of x in constraints + // recompute bounds for every of these sub-expressions + // propagate recomputed bounds up if they are stronger + // For sub-expressions that correspond to factors, perform + // bounds propagations + // For sub-expressions that correspond to constraints check that bounds have + // non-empty intersection with the bound imposed by the constraint + + // 1. Life-time management of sub-expressions and insertion into propagation queue + void stellensatz2::push_constraint(lp::constraint_index ci) { + auto [p, k] = m_constraints[ci]; + m_scopes.push_back( + {m_parent_trail.size(), + m_factor_trail.size(), + m_polynomial_queue.size(), + m_interval_trail.size(), + m_prop_qhead}); + + insert_parents(p); + insert_parents(p, ci); + for (auto x : p.free_vars()) { + auto f = factor(x, ci); + if (f.degree > 1) + continue; + insert_parents(f.p); + insert_parents(f.q); + insert_factor(f.p, x, f, ci); + insert_factor(f.q, x, f, ci); + } + + // 2 new bound is added, then update intervals and queue up for propagation + auto const &b = m_bounds[ci]; + auto &m = m_di.num_manager(); + scoped_dep_interval new_iv(m_di); + bool is_strict = k == lp::lconstraint_kind::GT; + + if (b.q.is_var()) { + auto const &iv = get_interval(b.q); + m_di.set_lower_is_inf(new_iv, false); + m_di.set_lower(new_iv, b.m_value); + m_di.set_lower_is_open(new_iv, is_strict); + m_di.set_lower_dep(new_iv, b.d); + m_di.copy_upper_bound(iv, new_iv); + if (update_interval(new_iv, b.q)) + m_polynomial_queue.push_back(b.q); + } + if (b.mq.is_var()) { + auto const &iv = get_interval(b.mq); + m_di.set_upper_is_inf(new_iv, false); + m_di.set_upper_is_open(new_iv, is_strict); + m_di.set_upper_dep(new_iv, b.d); + m_di.set_upper(new_iv, -b.m_value); + m_di.copy_lower_bound(iv, new_iv); + if (update_interval(new_iv, b.mq)) + m_polynomial_queue.push_back(b.mq); + } + } + + + void stellensatz2::insert_parents(dd::pdd const &p, lp::constraint_index ci) { + m_parent_constraints.reserve(p.index() + 1); + m_parent_constraints[p.index()].push_back(ci); + } + + void stellensatz2::insert_parents(dd::pdd const &p) { + if (m_is_parent.get(p.index(), false) || p.is_val()) + return; + m_is_parent.setx(p.index(), true, false); + m_parent_trail.push_back(p); + m_intervals.reserve(p.index() + 1); + m_intervals[p.index()].push_back(alloc(dep_interval)); + insert_child(p.lo(), p); + insert_child(p.hi(), p); + } + + void stellensatz2::insert_child(dd::pdd const &child, dd::pdd const &parent) { + m_parents.reserve(child.index() + 1); + m_parents[child.index()].push_back(parent); + insert_parents(child); + } + + void stellensatz2::insert_factor(dd::pdd const &p, lpvar x, factorization const &f, lp::constraint_index ci) { + m_factors[p.index()].push_back({x, f, ci}); + m_factor_trail.push_back(p.index()); + } + + void stellensatz2::pop_constraint(lp::constraint_index ci) { + auto const &s = m_scopes[ci]; + m_parent_constraints[ci].pop_back(); + while (m_factor_trail.size() >= s.factors_lim) { + auto p_index = m_factor_trail.back(); + m_factors[p_index].pop_back(); + m_factor_trail.pop_back(); + } + while (m_parent_trail.size() >= s.parents_lim) { + auto p = m_parent_trail.back(); + m_is_parent[p.index()] = false; + m_parents[p.lo().index()].pop_back(); + m_parents[p.hi().index()].pop_back(); + m_parent_trail.pop_back(); + } + m_polynomial_queue.shrink(s.polynomial_lim); + + while (m_interval_trail.size() > s.interval_lim) { + auto p_index = m_interval_trail.back(); + m_intervals[p_index].pop_back(); + m_interval_trail.pop_back(); + } + m_prop_qhead = s.qhead; + m_scopes.pop_back(); + } + + + void stellensatz2::propagate() { + while (m_prop_qhead < m_polynomial_queue.size()) { + auto p = m_polynomial_queue[m_prop_qhead++]; + propagate_intervals(p); + } + } + + bool stellensatz2::is_better(dep_interval const &new_iv, dep_interval const &old_iv) { + auto &m = m_di.num_manager(); + if (old_iv.m_lower_inf && !new_iv.m_lower_inf) + return true; + if (old_iv.m_upper_inf && !new_iv.m_upper_inf) + return true; + if (!new_iv.m_lower_inf && !old_iv.m_lower_inf && + (m.gt(new_iv.m_lower, old_iv.m_lower) || + (m.eq(new_iv.m_lower, old_iv.m_lower) && new_iv.m_lower_open && !old_iv.m_lower_open))) + return true; + + if (!new_iv.m_upper_inf && !old_iv.m_upper_inf && + (m.gt(old_iv.m_upper, new_iv.m_upper) || + (m.eq(new_iv.m_upper, old_iv.m_upper) && new_iv.m_upper_open && !old_iv.m_upper_open))) + return true; + + return false; + } + + dep_interval const &stellensatz2::get_interval(dd::pdd const &p) { + auto &ivs = m_intervals[p.index()]; + if (ivs.empty() && p.is_val()) { + ivs.push_back(alloc(dep_interval)); + m_di.set_value(*ivs.back(), p.val()); + return *ivs.back(); + } + SASSERT(!ivs.empty()); + return *ivs.back(); + } + + bool stellensatz2::update_interval(dep_interval const &new_iv, dd::pdd const &p) { + auto &old_iv = get_interval(p); + + if (!is_better(new_iv, old_iv)) + return false; + + m_intervals[p.index()].push_back(alloc(dep_interval)); + m_di.set(*m_intervals[p.index()].back(), new_iv); + m_interval_trail.push_back(p.index()); + return true; + } + + void stellensatz2::propagate_intervals(dd::pdd const &p) { + // for each parent of p, propagate updated intervals + for (auto parent : m_parents[p.index()]) { + SASSERT(!parent.is_val()); + scoped_dep_interval new_iv(m_di); + auto &x = get_interval(pddm.mk_var(parent.var())); + auto &lo = get_interval(parent.lo()); + auto &hi = get_interval(parent.hi()); + calculate_interval(new_iv, x, lo, hi); + if (update_interval(new_iv, parent)) + m_polynomial_queue.push_back(parent); + } + + // check constraints + for (auto ci : m_parent_constraints[p.index()]) { + auto const &[poly, k] = m_constraints[ci]; + auto &iv_poly = get_interval(poly); + if (constraint_is_bound_conflict(ci, iv_poly)) + return; // conflict detected + } + + // propagate to factors + for (auto const &[x, f, ci] : m_factors[p.index()]) { + scoped_dep_interval iv_mq(m_di); + auto &iv_p = get_interval(f.p); + auto &iv_q = get_interval(f.q); + m_di.neg(iv_q, iv_mq); + auto [p, k] = m_constraints[ci]; // p == f.p * x + f.q + SASSERT(p == (f.p * pddm.mk_var(x)) + f.q); + rational value; + svector cs; + if (constraint_is_propagating(iv_p, iv_mq, x, cs, k, value)) + propagate_constraint(x, k, value, ci, cs); + } + + } +} \ No newline at end of file diff --git a/src/math/lp/nla_stellensatz2.h b/src/math/lp/nla_stellensatz2.h index 3b535f000..f095a76a0 100644 --- a/src/math/lp/nla_stellensatz2.h +++ b/src/math/lp/nla_stellensatz2.h @@ -222,7 +222,11 @@ namespace nla { bool should_propagate() const { return m_prop_qhead < m_constraints.size(); } // assuming variables have bounds determine if polynomial has lower/upper bounds - void interval(dd::pdd p, scoped_dep_interval &iv); + void calculate_interval(scoped_dep_interval &out, dd::pdd p); + void calculate_interval(scoped_dep_interval &out, lpvar x); + void calculate_interval(scoped_dep_interval &out, dep_interval const& x, dep_interval const&lo, dep_interval const&hi); + void retrieve_interval(scoped_dep_interval &out, dd::pdd const &p); + void retrieve_interval(scoped_dep_interval &out, lpvar v); void set_conflict(lp::constraint_index ci) { m_conflict = ci; @@ -240,6 +244,27 @@ namespace nla { indexed_uint_set m_tabu; unsigned_vector m_var2level, m_level2var; + // propagation + struct scope { + unsigned parents_lim, factors_lim, polynomial_lim, interval_lim, qhead; + }; + + struct factor_prop { + lpvar x; + factorization f; + lp::constraint_index ci; + }; + vector m_scopes; + vector m_parent_trail; + vector> m_parents; + vector> m_factors; + vector m_polynomial_queue; + unsigned_vector m_interval_trail; + unsigned_vector m_factor_trail; + vector> m_parent_constraints; + vector> m_intervals; + bool_vector m_is_parent; + void push_bound(lp::constraint_index ci); unsigned get_lower(lpvar v) const { return get_lower(pddm.mk_var(v)); } @@ -320,6 +345,21 @@ namespace nla { lp::constraint_index resolve_variable(lpvar x, lp::constraint_index ci, lp::constraint_index other_ci); lp::constraint_index resolve(lp::constraint_index c1, lp::constraint_index c2); + // propagation + void push_constraint(lp::constraint_index ci); + void insert_parents(dd::pdd const &p, lp::constraint_index ci); + void insert_parents(dd::pdd const &p); + void insert_child(dd::pdd const &child, dd::pdd const &parent); + void insert_factor(dd::pdd const &p, lpvar x, factorization const &f, lp::constraint_index ci); + void pop_constraint(lp::constraint_index ci); + bool is_better(dep_interval const &new_iv, dep_interval const &old_iv); + bool update_interval(dep_interval const &new_iv, dd::pdd const &p); + dep_interval const &get_interval(dd::pdd const &p); + void propagate_intervals(dd::pdd const &p); + void propagate_constraint(lpvar x, lp::lconstraint_kind k, rational const &value, lp::constraint_index ci, svector &cs); + + // constraints + bool propagation_cycle(lpvar v, rational const& value, lp::lconstraint_kind k, lp::constraint_index ci, svector& cs) const; bool constraint_is_true(lp::constraint_index ci) const; bool constraint_is_true(constraint const &c) const; @@ -330,11 +370,15 @@ namespace nla { bool constraint_is_trivial(lp::constraint_index ci) const; bool constraint_is_trivial(constraint const& c) const; bool constraint_is_bound_conflict(lp::constraint_index ci); + bool constraint_is_bound_conflict(lp::constraint_index ci, dep_interval const &iv); bool var_is_bound_conflict(lpvar v); bool is_bound_conflict(dd::pdd const &p); bool constraint_is_propagating(lp::constraint_index ci, svector &cs, lpvar &v, lp::lconstraint_kind &k, rational &value); + bool constraint_is_propagating(dep_interval const &ivp, dep_interval const& ivq, lpvar x, + svector &cs, lp::lconstraint_kind &k, rational &value); + lbool sign(dd::pdd const &p);