diff --git a/src/math/lp/nla_stellensatz.cpp b/src/math/lp/nla_stellensatz.cpp index 3556fd693..82d7067ea 100644 --- a/src/math/lp/nla_stellensatz.cpp +++ b/src/math/lp/nla_stellensatz.cpp @@ -95,20 +95,32 @@ namespace nla { lbool stellensatz::saturate() { init_solver(); TRACE(arith, display(tout << "stellensatz before saturation\n")); + start_saturate: lbool r; #if 1 r = conflict_saturation(); - if (r == l_false) - return r; #else r = model_repair(); - if (r == l_false) - return r; #endif - r = m_solver.solve(); - // IF_VERBOSE(0, verbose_stream() << "stellensatz " << r << "\n"); + if (r != l_false) + r = m_solver.solve(m_core); TRACE(arith, display(tout << "stellensatz after saturation " << r << "\n")); + if (r == l_false) { + while (backtrack(m_core)) { + lbool rb = m_solver.solve(m_core); + if (rb == l_false) + continue; + if (rb == l_undef) // TODO: if there is a core that doesn't depend on new monomials we could use this for conflicts + return rb; + + m_solver.update_values(m_values); + goto start_saturate; + } + conflict(m_core); + } + if (r == l_true) + r = l_undef; return r; } @@ -282,6 +294,8 @@ namespace nla { return true; if (p_value == 0) return true; + if (m_tabu[ci].contains(x)) + return true; unsigned num_x = m_occurs[x].size(); for (unsigned i = 0; i < f.degree; ++i) f.p *= to_pdd(x); @@ -289,6 +303,8 @@ namespace nla { for (unsigned cx = 0; cx < num_x; ++cx) { auto other_ci = m_occurs[x][cx]; + if (m_tabu[other_ci].contains(x)) + continue; if (!resolve_variable(x, ci, other_ci, p_value, f, _m1, _f_p)) return false; } @@ -336,10 +352,12 @@ namespace nla { f_p = f_p.mul(m1); g_p = g_p.mul(m2); + #if 0 if (!has_term_offset(f_p)) return true; if (!has_term_offset(g_p)) return true; + #endif TRACE(arith, tout << "m1 " << m1 << " m2 " << m2 << " m1m2: " << m1m2 << "\n"); @@ -378,39 +396,44 @@ namespace nla { return true; if (m_constraints[new_ci].p.degree() <= 3) init_occurs(new_ci); - TRACE(arith, tout << "eliminate j" << x << ":\n"; display_constraint(tout << "ci: ", ci) << "\n"; + TRACE(arith, tout << "eliminate j" << x << ":\n"; + display_constraint(tout << "ci: ", ci) << "\n"; display_constraint(tout << "other_ci: ", other_ci) << "\n"; - display_constraint(tout << "ci_a: ", ci_a) << "\n"; display_constraint(tout << "ci_b: ", ci_b) << "\n"; + display_constraint(tout << "ci_a: ", ci_a) << "\n"; + display_constraint(tout << "ci_b: ", ci_b) << "\n"; display_constraint(tout << "new_ci: ", new_ci) << "\n"); if (conflict(new_ci)) return false; - auto const &[new_p, new_k, new_j] = m_constraints[new_ci]; - if (new_p.degree() <= 3 && !new_p.free_vars().contains(x)) { - uint_set new_tabu(m_tabu[ci]); - new_tabu.insert(x); - add_active(new_ci, new_tabu); - } - return true; + uint_set new_tabu(m_tabu[ci]); + new_tabu.insert(x); + add_active(new_ci, new_tabu); + return factor(new_ci) != l_false; } bool stellensatz::conflict(lp::constraint_index ci) { auto const &[new_p, new_k, new_j] = m_constraints[ci]; if (new_p.is_val() && ((new_p.val() < 0 && new_k == lp::lconstraint_kind::GE) || (new_p.val() <= 0 && new_k == lp::lconstraint_kind::GT))) { - lp::explanation ex; - lemma_builder new_lemma(c(), "stellensatz"); - m_constraints_in_conflict.reset(); - explain_constraint(new_lemma, ci, ex); - new_lemma &= ex; - IF_VERBOSE(2, verbose_stream() << "stellensatz unsat \n" << new_lemma << "\n"); - TRACE(arith, tout << "unsat\n" << new_lemma << "\n"); - c().lra_solver().settings().stats().m_nla_stellensatz++; + m_core.reset(); + m_core.push_back(ci); return true; } return false; } + void stellensatz::conflict(svector const& core) { + lp::explanation ex; + lemma_builder new_lemma(c(), "stellensatz"); + m_constraints_in_conflict.reset(); + for (auto ci : core) + explain_constraint(new_lemma, ci, ex); + new_lemma &= ex; + IF_VERBOSE(2, verbose_stream() << "stellensatz unsat \n" << new_lemma << "\n"); + TRACE(arith, tout << "unsat\n" << new_lemma << "\n"); + c().lra_solver().settings().stats().m_nla_stellensatz++; + } + lbool stellensatz::model_repair() { for (lp::constraint_index ci = 0; ci < m_constraints.size(); ++ci) m_active.insert(ci); @@ -467,7 +490,8 @@ namespace nla { if (!resolve_variable(v, inf, s, p_value, f, m, f_p)) return false; for (auto i : infs) - assume_ge(v, i, inf); + if (!assume_ge(v, i, inf)) + return false; } else { auto f = factor(v, sup); @@ -478,20 +502,23 @@ namespace nla { for (auto i : infs) if (!resolve_variable(v, sup, i, p_value, f, m, f_p)) return false; - for (auto s : sups) - assume_ge(v, sup, s); + for (auto s : sups) + if (!assume_ge(v, sup, s)) + return false; } for (auto ci : infs) - m_active.remove(ci); + if (m_active.contains(ci)) + m_active.remove(ci); for (auto ci : sups) - m_active.remove(ci); + if (m_active.contains(ci)) + m_active.remove(ci); return true; } // lo <= hi - void stellensatz::assume_ge(lpvar v, lp::constraint_index lo, lp::constraint_index hi) { + bool stellensatz::assume_ge(lpvar v, lp::constraint_index lo, lp::constraint_index hi) { if (lo == hi) - return; + return true; auto const &[plo, klo, jlo] = m_constraints[lo]; auto const &[phi, khi, jhi] = m_constraints[hi]; auto f = factor(v, lo); @@ -520,6 +547,7 @@ namespace nla { SASSERT(value(r) >= 0); auto new_ci = assume(r, join(klo, khi)); m_active.insert(new_ci); + return factor(new_ci) != l_false; } std::pair stellensatz::find_bounds(lpvar v) { @@ -560,8 +588,6 @@ namespace nla { } bool stellensatz::vanishing(lpvar x, factorization const &f, lp::constraint_index ci) { - if (f.q.is_zero()) - return false; if (f.p.is_zero()) return false; auto p_value = value(f.p); @@ -591,7 +617,6 @@ namespace nla { auto [vars, q] = p.var_factors(); // p = vars * q auto add_new = [&](lp::constraint_index new_ci) { - SASSERT(!constraint_is_true(new_ci)); TRACE(arith, display_constraint(tout << "factor ", new_ci) << "\n"); if (conflict(new_ci)) return l_false; @@ -601,33 +626,12 @@ namespace nla { return l_true; }; - auto subst = [&](lpvar v, dd::pdd p) { - auto q = pddm.mk_var(v) - p; - auto new_ci = substitute(ci, assume(q, lp::lconstraint_kind::EQ), v, p); - TRACE(arith, tout << "j" << v << " " << p << "\n"; - display_constraint(tout, ci) << "\n"; - display_constraint(tout, new_ci) << "\n"); - return add_new(new_ci); - }; - - TRACE(arith, tout << "factor (" << ci << ") " << p << " q: " << q << " vars: " << vars << "\n"); - if (false && vars.empty()) { - for (auto v : p.free_vars()) { - if (value(v) == 0) - return subst(v, pddm.mk_val(0)); - if (value(v) == 1) - return subst(v, pddm.mk_val(1)); - if (value(v) == -1) - return subst(v, pddm.mk_val(rational(-1))); - } - return l_undef; - } + TRACE(arith, tout << "factor (" << ci << ") " << p << " q: " << q << " vars: " << vars << "\n"); if (vars.empty()) return l_undef; for (auto v : vars) { if (value(v) == 0) return l_undef; -// return subst(v, pddm.mk_val(0)); } vector muls; muls.push_back(q); @@ -640,6 +644,8 @@ namespace nla { auto k = value(vars[i]) > 0 ? lp::lconstraint_kind::GT : lp::lconstraint_kind::LT; new_ci = divide(new_ci, assume(pv, k), muls[i]); } + if (m_active.contains(ci)) + m_active.remove(ci); return add_new(new_ci); } @@ -669,54 +675,132 @@ namespace nla { } + // + // check if core depends on an assumption + // identify the maximal assumption + // undo m_constraints down to max_ci. + // negate max_ci + // propagate it using remaining external and assumptions. + // find a new satisfying assignment (if any) before continuing. + // + bool stellensatz::backtrack(svector const &core) { + m_constraints_in_conflict.reset(); + svector external, assumptions; + for (auto ci : core) + explain_constraint(ci, external, assumptions); + if (assumptions.empty()) + return false; + lp::constraint_index max_ci = 0; + for (auto ci : assumptions) + max_ci = std::max(ci, max_ci); + TRACE(arith, tout << "max " << max_ci << " external " << external << " assumptions " << assumptions << "\n"; + display_constraint(tout, max_ci);); + // TODO, we can in reality replay all constraints that don't depend on max_ci + vector replay; + unsigned i = 0; + for (auto ci : external) { + if (ci > max_ci) + replay.push_back(m_constraints[ci]); + else + external[i++] = ci; + } + external.shrink(i); + auto [p, k, j] = m_constraints[max_ci]; + while (m_constraints.size() > max_ci) { + auto const& [_p, _k, _j] = m_constraints.back(); + m_constraint_index.erase({_p.index(), _k}); + m_constraints.pop_back(); + auto ci = m_constraints.size(); + if (!m_occurs_trail.empty() && m_occurs_trail.back() == ci) { + remove_occurs(ci); + m_occurs_trail.pop_back(); + } + } + for (auto const &[_p, _k, _j] : replay) { + auto ci = add_constraint(_p, _k, _j); + init_occurs(ci); + external.push_back(ci); + } + assumptions.erase(max_ci); + external.append(assumptions); + propagation_justification prop{external}; + auto new_ci = add_constraint(p, negate(k), prop); + TRACE(arith, display_constraint(tout << "backtrack ", new_ci) << "\n"); + init_occurs(new_ci); + return true; + } + + void stellensatz::explain_constraint(lp::constraint_index ci, svector &external, + svector &assumptions) { + if (m_constraints_in_conflict.contains(ci)) + return; + m_constraints_in_conflict.insert(ci); + auto &[p, k, b] = m_constraints[ci]; + if (std::holds_alternative(b)) { + external.push_back(ci); + } + else if (std::holds_alternative(b)) { + auto &m = std::get(b); + explain_constraint(m.left, external, assumptions); + explain_constraint(m.right, external, assumptions); + } + else if (std::holds_alternative(b)) { + auto &m = std::get(b); + explain_constraint(m.left, external, assumptions); + explain_constraint(m.right, external, assumptions); + } + else if (std::holds_alternative(b)) { + auto &m = std::get(b); + explain_constraint(m.ci, external, assumptions); + explain_constraint(m.divisor, external, assumptions); + } + else if (std::holds_alternative(b)) { + auto &m = std::get(b); + explain_constraint(m.ci, external, assumptions); + explain_constraint(m.ci_eq, external, assumptions); + } + else if (std::holds_alternative(b)) { + auto &m = std::get(b); + for (auto c : m.cs) + explain_constraint(c, external, assumptions); + } + else if (std::holds_alternative(b)) { + auto &m = std::get(b); + explain_constraint(m.left, external, assumptions); + explain_constraint(m.right, external, assumptions); + } + else if (std::holds_alternative(b)) { + auto &m = std::get(b); + explain_constraint(m.ci, external, assumptions); + } + else if (std::holds_alternative(b)) { + assumptions.push_back(ci); + } + else if (std::holds_alternative(b)) { + auto &m = std::get(b); + explain_constraint(m.ci, external, assumptions); + } + else + UNREACHABLE(); + } // // a constraint can be explained by a set of bounds used as justifications // and by an original constraint. // void stellensatz::explain_constraint(lemma_builder &new_lemma, lp::constraint_index ci, lp::explanation &ex) { - if (m_constraints_in_conflict.contains(ci)) - return; - m_constraints_in_conflict.insert(ci); - auto &[p, k, b] = m_constraints[ci]; - if (std::holds_alternative(b)) { + svector external, assumptions; + explain_constraint(ci, external, assumptions); + for (auto ci : external) { + auto &[p, k, b] = m_constraints[ci]; auto dep = std::get(b); c().lra_solver().push_explanation(dep.dep, ex); - } - else if (std::holds_alternative(b)) { - auto& m = std::get(b); - explain_constraint(new_lemma, m.left, ex); - explain_constraint(new_lemma, m.right, ex); - } - else if (std::holds_alternative(b)) { - auto &m = std::get(b); - explain_constraint(new_lemma, m.ci, ex); - explain_constraint(new_lemma, m.divisor, ex); - } - else if (std::holds_alternative(b)) { - auto &m = std::get(b); - explain_constraint(new_lemma, m.ci, ex); - explain_constraint(new_lemma, m.ci_eq, ex); - } - else if (std::holds_alternative(b)) { - auto& m = std::get(b); - explain_constraint(new_lemma, m.left, ex); - explain_constraint(new_lemma, m.right, ex); - } - else if (std::holds_alternative(b)) { - auto& m = std::get(b); - explain_constraint(new_lemma, m.ci, ex); - } - else if (std::holds_alternative(b)) { - auto [t, coeff] = to_term_offset(p); + } + for (auto ci : assumptions) { + auto &[p, k, b] = m_constraints[ci]; + auto [t, coeff] = to_term_offset(p); new_lemma |= ineq(t, negate(k), -coeff); } - else if (std::holds_alternative(b)) { - auto& m = std::get(b); - explain_constraint(new_lemma, m.ci, ex); - } - else - UNREACHABLE(); } rational stellensatz::value(dd::pdd const& p) const { @@ -765,6 +849,37 @@ namespace nla { } lp::constraint_index stellensatz::assume(dd::pdd const& p, lp::lconstraint_kind k) { + if (k == lp::lconstraint_kind::EQ) { + auto left = assume(p, lp::lconstraint_kind::GE); + auto right = assume(-p, lp::lconstraint_kind::GE); + return add_constraint(p, k, eq_justification{left, right}); + } + u_dependency *d = nullptr; + auto has_bound = [&](rational a, lpvar x, rational b) { + rational bound; + bool is_strict = false; + if (a == 1 && k == lp::lconstraint_kind::GE && c().lra_solver().has_lower_bound(x, d, bound, is_strict) && + bound >= -b) { + return true; + } + if (a == 1 && k == lp::lconstraint_kind::GT && c().lra_solver().has_lower_bound(x, d, bound, is_strict) && + (bound > -b || (is_strict && bound >= -b))) { + return true; + } + if (a == -1 && k == lp::lconstraint_kind::GE && c().lra_solver().has_upper_bound(x, d, bound, is_strict) && + bound <= -b) { + return true; + } + if (a == -1 && k == lp::lconstraint_kind::GT && c().lra_solver().has_upper_bound(x, d, bound, is_strict) && + (bound < -b || (is_strict && bound <= -b))) { + return true; + } + return false; + }; + + if (p.is_unilinear() && has_bound(p.hi().val(), p.var(), p.lo().val())) + // ax + b ~k~ 0 + return add_constraint(p, k, external_justification(d)); return add_constraint(p, k, assumption_justification()); } @@ -820,12 +935,19 @@ namespace nla { void stellensatz::init_occurs(lp::constraint_index ci) { if (ci == lp::null_ci) return; + m_occurs_trail.push_back(ci); auto const &con = m_constraints[ci]; for (auto v : con.p.free_vars()) m_occurs[v].push_back(ci); } + void stellensatz::remove_occurs(lp::constraint_index ci) { + auto const &con = m_constraints[ci]; + for (auto v : con.p.free_vars()) + m_occurs[v].pop_back(); + } + bool stellensatz::is_int(svector const& vars) const { return all_of(vars, [&](lpvar v) { return c().lra_solver().var_is_int(v); }); } @@ -922,10 +1044,20 @@ namespace nla { auto m = std::get(j); out << "(" << m.left << ") + (" << m.right << ")"; } + else if (std::holds_alternative(j)) { + auto &m = std::get(j); + out << "(" << m.left << ") & (" << m.right << ")"; + } else if (std::holds_alternative(j)) { auto m = std::get(j); out << "(" << m.ci << ") (" << m.ci_eq << ") by j" << m.v << " := " << m.p; } + else if (std::holds_alternative(j)) { + auto &m = std::get(j); + out << "propagate "; + for (auto c : m.cs) + out << "(" << c << ") "; + } else if (std::holds_alternative(j)) { auto m = std::get(j); out << "(" << m.left << ") * (" << m.right << ")"; @@ -970,6 +1102,15 @@ namespace nla { todo.push_back(m.left); todo.push_back(m.right); } + if (std::holds_alternative(j)) { + auto m = std::get(j); + todo.push_back(m.left); + todo.push_back(m.right); + } + if (std::holds_alternative(j)) { + auto m = std::get(j); + todo.append(m.cs); + } else if (std::holds_alternative(j)) { auto m = std::get(j); todo.push_back(m.ci); @@ -1052,40 +1193,19 @@ namespace nla { return to; } - - // - // convert a conflict from m_solver.lra()/lia() into - // a conflict for c().lra_solver() - // - void stellensatz::solver::add_lemma() { - TRACE(arith, s.display(tout); s.display_lemma(tout, m_ex)); - auto &src = s.c().lra_solver(); - lp::explanation ex2; - lemma_builder new_lemma(s.c(), "stellensatz"); - s.m_constraints_in_conflict.reset(); - for (auto p : m_ex) - s.explain_constraint(new_lemma, m_internal2external_constraints[p.ci()], ex2); - new_lemma &= ex2; - IF_VERBOSE(2, verbose_stream() << "stellensatz unsat \n" << new_lemma << "\n"); - TRACE(arith, tout << "unsat\n" << new_lemma << "\n"); - s.c().lra_solver().settings().stats().m_nla_stellensatz++; - } - - lbool stellensatz::solver::solve() { + lbool stellensatz::solver::solve(svector& core) { init(); lbool r = solve_lra(); - // IF_VERBOSE(0, verbose_stream() << "solve lra " << r << "\n"); - if (r != l_true) - return r; - r = solve_lia(); - // IF_VERBOSE(0, verbose_stream() << "solve lia " << r << "\n"); - if (r != l_true) - return r; - return l_undef; - if (update_values()) - return l_true; - - return l_undef; + if (r == l_true) + r = solve_lia(); + + if (r == l_false) { + core.reset(); + for (auto p : m_ex) + core.push_back(m_internal2external_constraints[p.ci()]); + return l_false; + } + return r; } lbool stellensatz::solver::solve_lra() { @@ -1094,7 +1214,6 @@ namespace nla { return l_true; if (status == lp::lp_status::INFEASIBLE) { lra_solver->get_infeasibility_explanation(m_ex); - add_lemma(); return l_false; } return l_undef; @@ -1105,7 +1224,6 @@ namespace nla { case lp::lia_move::sat: return l_true; case lp::lia_move::conflict: - add_lemma(); return l_false; default: // TODO: an option is to perform (bounded) search here to get an LIA verdict. return l_undef; @@ -1113,37 +1231,12 @@ namespace nla { return l_undef; } - // update m_values - // return true if the new assignment satisfies the products. - // return false if value constraints are not satisfied on monomials and there is a false constaint. - bool stellensatz::solver::update_values() { - return false; - #if 0 - std::unordered_map values; - lra_solver->get_model(values); - unsigned sz = lra_solver->number_of_vars(); - for (unsigned i = 0; i < sz; ++i) { - auto const &value = values[i]; - bool is_int = lra_solver->var_is_int(i); - SASSERT(!is_int || value.is_int()); - if (!s.is_mon_var(i)) - continue; - rational val(1); - for (auto w : s.c().emons()[i]) - val *= values[w]; - values[i] = val; - } - auto indices = lra_solver->constraints().indices(); - bool satisfies_products = all_of(indices, [&](auto ci) { - NOT_IMPLEMENTED_YET(); - // todo: wrong, needs to be at level of lra constraint evaluation - return s.constraint_is_true(ci);}); - if (satisfies_products) { - TRACE(arith, tout << "found satisfying model\n"); - for (auto v : s.m_coi.vars()) - s.c().lra_solver().set_column_value(v, lp::impq(values[v], rational(0))); - } - return satisfies_products; - #endif + // update values using the model + void stellensatz::solver::update_values(vector& values) { + std::unordered_map _values; + lra_solver->get_model(_values); + unsigned sz = values.size(); + for (unsigned i = 0; i < sz; ++i) + values[i] = _values[i]; } } diff --git a/src/math/lp/nla_stellensatz.h b/src/math/lp/nla_stellensatz.h index 4448bb05e..6ecfb950d 100644 --- a/src/math/lp/nla_stellensatz.h +++ b/src/math/lp/nla_stellensatz.h @@ -41,9 +41,16 @@ namespace nla { lpvar v; dd::pdd p; }; + struct eq_justification { + lp::constraint_index left; + lp::constraint_index right; + }; struct gcd_justification { lp::constraint_index ci; }; + struct propagation_justification { + svector cs; + }; using justification = std::variant< external_justification, @@ -52,6 +59,8 @@ namespace nla { substitute_justification, addition_justification, division_justification, + eq_justification, + propagation_justification, multiplication_poly_justification, multiplication_justification>; @@ -97,13 +106,13 @@ namespace nla { monomial_factory m_monomial_factory; lbool solve_lra(); lbool solve_lia(); - bool update_values(); + void init(); - void add_lemma(); term_offset to_term_offset(dd::pdd const &p); public: solver(stellensatz &s) : s(s) {}; - lbool solve(); + lbool solve(svector& core); + void update_values(vector& values); }; solver m_solver; @@ -123,6 +132,7 @@ namespace nla { indexed_uint_set m_active; vector m_tabu; vector m_values; + svector m_core, m_occurs_trail; struct constraint_key { unsigned pdd; @@ -157,10 +167,12 @@ namespace nla { void init_vars(); void init_occurs(); void init_occurs(lp::constraint_index ci); + void remove_occurs(lp::constraint_index ci); lbool conflict_saturation(); lbool factor(lp::constraint_index ci); bool conflict(lp::constraint_index ci); + void conflict(svector const& core); bool vanishing(lpvar x, factorization const& f, lp::constraint_index ci); lp::lpvar select_variable_to_eliminate(lp::constraint_index ci); unsigned degree_of_var_in_constraint(lpvar v, lp::constraint_index ci) const; @@ -173,11 +185,11 @@ namespace nla { bool model_repair(lp::lpvar v); struct bound_info { rational value; - lp::constraint_index bound; + lp::constraint_index bound = lp::null_ci; svector bounds; }; std::pair find_bounds(lpvar v); - void assume_ge(lpvar v, lp::constraint_index lo, lp::constraint_index hi); + bool assume_ge(lpvar v, lp::constraint_index lo, lp::constraint_index hi); bool constraint_is_true(lp::constraint_index ci) const; bool is_new_constraint(lp::constraint_index ci) const; @@ -200,6 +212,9 @@ namespace nla { // lemmas indexed_uint_set m_constraints_in_conflict; void explain_constraint(lemma_builder& new_lemma, lp::constraint_index ci, lp::explanation &ex); + void explain_constraint(lp::constraint_index ci, svector &external, + svector &assumptions); + bool backtrack(svector const& core); struct pp_j { stellensatz const &s; diff --git a/src/math/lp/nra_solver.cpp b/src/math/lp/nra_solver.cpp index 1658863c0..6e233514d 100644 --- a/src/math/lp/nra_solver.cpp +++ b/src/math/lp/nra_solver.cpp @@ -26,6 +26,12 @@ typedef nla::mon_eq mon_eq; typedef nla::variable_map_type variable_map_type; struct solver::imp { + struct model_bound { + lp::lpvar v; + rational r; + bool is_lower; + }; + lp::lar_solver& lra; reslimit& m_limit; params_ref m_params; @@ -35,6 +41,8 @@ struct solver::imp { scoped_ptr m_tmp1, m_tmp2; nla::coi m_coi; nla::core& m_nla_core; + unsigned m_max_constraint_index = 0; + vector m_model_bounds; imp(lp::lar_solver& s, reslimit& lim, params_ref const& p, nla::core& nla_core): lra(s), @@ -53,6 +61,8 @@ struct solver::imp { m_nlsat = alloc(nlsat::solver, m_limit, m_params, false); m_values = alloc(scoped_anum_vector, am()); m_lp2nl.reset(); + m_model_bounds.reset(); + m_max_constraint_index = 0; } // Create polynomial definition for variable v used in setup_assignment_solver. @@ -136,6 +146,20 @@ struct solver::imp { add_term(i); } + + void setup_model_bounds() { + for (unsigned v = 0; v < lra.number_of_vars(); ++v) { + if (!lra.column_is_int(v)) + continue; + auto value = m_nla_core.val(v); + SASSERT(value.is_int()); + unsigned sz = m_model_bounds.size(); + m_model_bounds.push_back({v, value, true}); + m_model_bounds.push_back({v, value, false}); + m_nlsat->track_model_value(v, value, this + sz, this + sz + 1); + } + } + /** \brief one-shot nlsat check. A one shot checker is the least functionality that can @@ -151,8 +175,6 @@ struct solver::imp { SASSERT(need_check()); reset(); vector core; - - smt_params_helper p(m_params); @@ -161,6 +183,9 @@ struct solver::imp { else setup_solver_terms(); + if (p.arith_nl_nra_model_bounds()) + setup_model_bounds(); + TRACE(nra, m_nlsat->display(tout)); if (p.arith_nl_log()) { @@ -222,11 +247,21 @@ struct solver::imp { case l_false: { lp::explanation ex; m_nlsat->get_core(core); + nla::lemma_builder lemma(m_nla_core, __FUNCTION__); for (auto c : core) { unsigned idx = static_cast(static_cast(c) - this); - ex.push_back(idx); + if (idx <= m_max_constraint_index) + ex.push_back(idx); + else { + idx -= m_max_constraint_index; + auto const& [v, bound, is_lower] = m_model_bounds[idx]; + if (is_lower) + lemma |= nla::ineq(v, lp::lconstraint_kind::LE, bound - 1); + else + lemma |= nla::ineq(v, lp::lconstraint_kind::GE, bound + 1); + } } - nla::lemma_builder lemma(m_nla_core, __FUNCTION__); + lemma &= ex; m_nla_core.set_use_nra_model(true); TRACE(nra, tout << lemma << "\n"); @@ -278,6 +313,7 @@ struct solver::imp { } void add_constraint(unsigned idx) { + m_max_constraint_index = std::max(m_max_constraint_index, idx); auto& c = lra.constraints()[idx]; auto& pm = m_nlsat->pm(); auto k = c.kind(); diff --git a/src/nlsat/nlsat_params.pyg b/src/nlsat/nlsat_params.pyg index b035f4189..ebe18a1a5 100644 --- a/src/nlsat/nlsat_params.pyg +++ b/src/nlsat/nlsat_params.pyg @@ -8,6 +8,7 @@ def_module_params('nlsat', ('cell_sample', BOOL, True, "cell sample projection"), ('lazy', UINT, 0, "how lazy the solver is."), ('reorder', BOOL, True, "reorder variables."), + ('model_bounds', BOOL, False, "constraint integer search using model bounds."), ('log_lemmas', BOOL, False, "display lemmas as self-contained SMT formulas"), ('dump_mathematica', BOOL, False, "display lemmas as matematica"), ('check_lemmas', BOOL, False, "check lemmas on the fly using an independent nlsat solver"), diff --git a/src/nlsat/nlsat_solver.cpp b/src/nlsat/nlsat_solver.cpp index bad981011..2016d19bc 100644 --- a/src/nlsat/nlsat_solver.cpp +++ b/src/nlsat/nlsat_solver.cpp @@ -1918,6 +1918,36 @@ namespace nlsat { << " :propagations " << m_stats.m_propagations << " :clauses " << m_clauses.size() << " :learned " << m_learned.size() << ")\n"); + + + #if 0 + // todo, review, test + // cut on the first model value + if (!m_model_values.empty()) { + bool found = false; + for (auto const &[x, bound] : bounds) { + for (auto const &[mx, mvalue, lo, hi] : m_model_values) { + if (mx == x) { + polynomial_ref p(m_pm); + rational one(1); + bool is_even = false; + p = m_pm.mk_linear(1, &one, &x, mvalue); + auto* p1 = p.get(); + if (mvalue < bound) { + mk_external_clause(1, &mk_ineq_literal(atom::GT, 1, &p1, &is_even), assumption(lo)); + } + else { + mk_external_clause(1, &mk_ineq_literal(atom::LT, 1, &p1, &is_even), assumption(hi)); + } + found = true; + break; + } + } + if (found) + break; + } + } + #endif for (auto const& b : bounds) { var x = b.first; rational lo = b.second; @@ -1944,6 +1974,19 @@ namespace nlsat { return r; } + struct model_value { + var x; + rational v; + _assumption_set lo, hi; + }; + vector m_model_values; + void track_model_value(var x, rational const &v, _assumption_set lo, _assumption_set hi) { + m_model_values.push_back(model_value{x, v, lo, hi}); + } + void reset_model_values() { + m_model_values.reset(); + } + bool m_reordered = false; bool simple_check() { literal_vector learned_unit; @@ -4388,6 +4431,14 @@ namespace nlsat { assumption solver::join(assumption a, assumption b) { return (m_imp->m_asm.mk_join(static_cast(a), static_cast(b))); } + + void solver::track_model_value(var x, rational const& v, assumption lo, assumption hi) { + m_imp->track_model_value(x, v, static_cast(lo), static_cast(hi)); + } + + void solver::reset_model_values() { + m_imp->reset_model_values(); + } }; diff --git a/src/nlsat/nlsat_solver.h b/src/nlsat/nlsat_solver.h index ae403fc83..7f8c37a8a 100644 --- a/src/nlsat/nlsat_solver.h +++ b/src/nlsat/nlsat_solver.h @@ -133,6 +133,10 @@ namespace nlsat { */ void mk_clause(unsigned num_lits, literal * lits, assumption a = nullptr); + void track_model_value(var x, rational const &v, assumption lo, assumption hi); + + void reset_model_values(); + // ----------------------- // // Basic diff --git a/src/params/smt_params_helper.pyg b/src/params/smt_params_helper.pyg index f343082d0..0aaa81f99 100644 --- a/src/params/smt_params_helper.pyg +++ b/src/params/smt_params_helper.pyg @@ -92,6 +92,7 @@ def_module_params(module_name='smt', ('arith.nl.stellensatz', BOOL, False, 'enable stellensatz saturation'), ('arith.nl.linearize', BOOL, True, 'enable NL linearization'), ('arith.nl.nra.poly', BOOL, False, 'use polynomial translation'), + ('arith.nl.nra.model_bounds', BOOL, False, 'use model-based integer bounds for nra solver'), ('arith.nl.log', BOOL, False, 'Log lemmas sent to nra solver'), ('arith.propagate_eqs', BOOL, True, 'propagate (cheap) equalities'), ('arith.epsilon', DOUBLE, 1.0, 'initial value of epsilon used for model generation of infinitesimals'),