diff --git a/src/math/polysat/fixplex.h b/src/math/polysat/fixplex.h index f6af50ccd..8cbf49956 100644 --- a/src/math/polysat/fixplex.h +++ b/src/math/polysat/fixplex.h @@ -252,12 +252,12 @@ namespace polysat { void fixed_var_eh(row const& r, var_t x); void eq_eh(var_t x, var_t y, row const& r1, row const& r2); lbool propagate_bounds(row const& r); - lbool propagate_bounds(ineq const& i); - lbool new_bound(row const& r, var_t x, mod_interval const& range); - lbool new_bound(ineq const& i, var_t x, numeral const& lo, numeral const& hi, u_dependency* a = nullptr, u_dependency* b = nullptr); - lbool conflict(ineq const& i, u_dependency* a = nullptr, u_dependency* b = nullptr); - lbool conflict(u_dependency* a); - lbool conflict(u_dependency* a, u_dependency* b) { return conflict(m_deps.mk_join(a, b)); } + bool propagate_bounds(ineq const& i); + bool new_bound(row const& r, var_t x, mod_interval const& range); + bool new_bound(ineq const& i, var_t x, numeral const& lo, numeral const& hi, u_dependency* a = nullptr, u_dependency* b = nullptr, u_dependency* c = nullptr, u_dependency* d = nullptr); + void conflict(ineq const& i, u_dependency* a = nullptr, u_dependency* b = nullptr, u_dependency* c = nullptr, u_dependency* d = nullptr); + void conflict(u_dependency* a); + void conflict(u_dependency* a, u_dependency* b, u_dependency* c = nullptr, u_dependency* d = nullptr) { conflict(m_deps.mk_join(m_deps.mk_join(a, b), m_deps.mk_join(c, d))); } u_dependency* row2dep(row const& r); void pivot(var_t x_i, var_t x_j, numeral const& b, numeral const& value); numeral value2delta(var_t v, numeral const& new_value) const; @@ -297,7 +297,9 @@ namespace polysat { // facilities for handling inequalities void add_ineq(var_t v, var_t w, unsigned dep, bool strict); void touch_var(var_t x); - lbool check_ineqs(); + bool ineqs_are_violated(); + bool ineqs_are_satisfied(); + void reset_ineqs_to_check(); bool is_solved(row const& r) const; bool is_solved(var_t v) const { SASSERT(is_base(v)); return is_solved(base2row(v)); } diff --git a/src/math/polysat/fixplex_def.h b/src/math/polysat/fixplex_def.h index 722621189..3c61494ec 100644 --- a/src/math/polysat/fixplex_def.h +++ b/src/math/polysat/fixplex_def.h @@ -149,13 +149,17 @@ namespace polysat { return l_false; case l_undef: m_to_patch.insert(v); - if (l_false == check_ineqs()) + if (ineqs_are_violated()) return l_false; return l_undef; } } SASSERT(well_formed()); - return check_ineqs(); + if (ineqs_are_violated()) + return l_false; + if (ineqs_are_satisfied()) + return l_true; + return l_undef; } template @@ -352,7 +356,6 @@ namespace polysat { return l_undef; } - pivot(x, y, b, new_value); return l_true; @@ -499,8 +502,7 @@ namespace polysat { } template - typename fixplex::numeral - fixplex::value2error(var_t v, numeral const& value) const { + typename fixplex::numeral fixplex::value2error(var_t v, numeral const& value) const { if (in_bounds(v)) return 0; SASSERT(lo(v) != hi(v)); @@ -519,6 +521,7 @@ namespace polysat { */ template void fixplex::set_bounds(var_t v, numeral const& l, numeral const& h, unsigned dep) { + ensure_var(v); update_bounds(v, l, h, mk_leaf(dep)); if (in_bounds(v)) return; @@ -571,6 +574,8 @@ namespace polysat { template void fixplex::add_ineq(var_t v, var_t w, unsigned dep, bool strict) { + ensure_var(v); + ensure_var(w); unsigned idx = m_ineqs.size(); m_var2ineqs.reserve(std::max(v, w) + 1); m_var2ineqs[v].push_back(idx); @@ -604,31 +609,51 @@ namespace polysat { } template - lbool fixplex::check_ineqs() { - m_vars_to_untouch.reset(); + void fixplex::reset_ineqs_to_check() { + for (auto idx : m_ineqs_to_check) { + if (idx >= m_ineqs.size()) + continue; + auto& ineq = m_ineqs[idx]; + m_var_is_touched.setx(ineq.v, false, false); + m_var_is_touched.setx(ineq.w, false, false); + ineq.is_active = false; + } + m_ineqs_to_check.reset(); + } + + /** + * Check if the current assignment satisfies the inequalities + */ + template + bool fixplex::ineqs_are_satisfied() { for (auto idx : m_ineqs_to_check) { if (idx >= m_ineqs.size()) continue; auto& ineq = m_ineqs[idx]; var_t v = ineq.v; var_t w = ineq.w; - if (ineq.strict && value(v) >= value(w)) { - IF_VERBOSE(0, verbose_stream() << "violated v" << v << " < v" << w << "\n"); - return l_false; - } - if (!ineq.strict && value(v) > value(w)) { - IF_VERBOSE(0, verbose_stream() << "violated v" << v << " <= v" << w << "\n"); - return l_false; - } - ineq.is_active = false; - m_vars_to_untouch.push_back(v); - m_vars_to_untouch.push_back(w); + if (ineq.strict && value(v) >= value(w)) + return false; + if (!ineq.strict && value(v) > value(w)) + return false; } - for (auto v : m_vars_to_untouch) - m_var_is_touched.set(v, false); - m_ineqs_to_check.reset(); - m_vars_to_untouch.reset(); - return l_true; + reset_ineqs_to_check(); + return true; + } + + /** + * Propagate bounds and check if the current inequalities are satisfied + */ + template + bool fixplex::ineqs_are_violated() { + for (unsigned i = 0; i < m_ineqs_to_check.size(); ++i) { + unsigned idx = m_ineqs_to_check[i]; + if (idx >= m_ineqs.size()) + continue; + if (!propagate_bounds(m_ineqs[idx])) + return true; + } + return false; } @@ -1087,16 +1112,16 @@ namespace polysat { m_var_eqs.push_back(var_eq(x, y, r1, r2)); } - template lbool fixplex::propagate_bounds() { lbool r = l_true; for (unsigned i = 0; r == l_true && i < m_rows.size(); ++i) r = propagate_bounds(row(i)); + if (r != l_true) + return r; for (auto ineq : m_ineqs) { - if (r != l_true) - break; - r = propagate_bounds(ineq); + if (!propagate_bounds(ineq)) + return l_false; } return r; } @@ -1132,7 +1157,7 @@ namespace polysat { if (free_v != null_var) { range = (-range) * free_c; - lbool res = new_bound(r, free_v, range); + lbool res = new_bound(r, free_v, range) ? l_true : l_false; SASSERT(in_bounds(free_v)); return res; } @@ -1140,7 +1165,7 @@ namespace polysat { var_t v = e.var(); SASSERT(!is_free(v)); auto range1 = range - m_vars[v] * e.coeff(); - lbool res = new_bound(r, v, range1); + lbool res = new_bound(r, v, range1) ? l_true : l_false; if (res != l_true) return res; // SASSERT(in_bounds(v)); @@ -1149,56 +1174,61 @@ namespace polysat { } template - lbool fixplex::propagate_bounds(ineq const& i) { + bool fixplex::propagate_bounds(ineq const& i) { // v < w & lo(v) + 1 = 0 -> conflict - // v < w & hi(w) = 0 -> conflict - // v < w & lo(v) >= hi(w) -> conflict - // v <= w & lo(v) > hi(w) -> conflict - // v < w & hi(v) >= hi(w) -> hi(v) := hi(w) - 1 - // v < w & lo(v) >= lo(w) -> lo(w) := lo(v) + 1 + // v < w & lo(w) = 0 & hi(w) = 1 -> conflict + // v < w & hi(w) != 0 & lo(w) <= hi(w) & hi(w) - 1 <= lo(v) -> conflict + // v <= w & hi(w) != 0 & lo(w) <= hi(w) & hi(w) <= lo(v) -> conflict + // v < w & hi(w) != 0 & lo(w) <= hi(w) <= hi(v) -> hi(v) := hi(w) - 1 + // v < w & lo(w) <= lo(v) -> lo(w) := lo(v) + 1 // v <= w & hi(v) > hi(w) -> hi(v) := hi(w) // v <= w & lo(v) > lo(w) -> lo(w) := lo(w) var_t v = i.v, w = i.w; bool s = i.strict; - if (s && lo(v) + 1 == 0) - return conflict(i, m_vars[v].m_lo_dep); - else if (s && hi(w) == 0) - return conflict(i, m_vars[v].m_hi_dep); - else if (s && lo(v) >= hi(w)) - return conflict(i, m_vars[v].m_lo_dep, m_vars[w].m_hi_dep); - else if (!s && lo(v) > hi(w)) - return conflict(i, m_vars[v].m_lo_dep, m_vars[w].m_hi_dep); - else if (s && hi(v) >= hi(w)) { - SASSERT(lo(v) < hi(w)); - SASSERT(hi(w) != 0); - return new_bound(i, v, lo(v), hi(w) - 1, m_vars[v].m_hi_dep, m_vars[w].m_hi_dep); - } - else if (s && lo(v) >= lo(w)) { - SASSERT(lo(v) + 1 != 0); - SASSERT(hi(w) > lo(v)); - return new_bound(i, w, lo(v) + 1, hi(w), m_vars[v].m_lo_dep, m_vars[w].m_lo_dep); - } - else if (!s && hi(v) > hi(w)) { - SASSERT(lo(v) <= hi(w)); - return new_bound(i, v, lo(v), hi(w), m_vars[v].m_hi_dep, m_vars[w].m_hi_dep); - } - else if (!s && lo(v) > lo(w)) { - SASSERT(lo(v) <= hi(w)); - return new_bound(i, w, lo(v), hi(w), m_vars[v].m_lo_dep, m_vars[w].m_lo_dep); - } - return l_true; + auto* vlo = m_vars[v].m_lo_dep, *vhi = m_vars[v].m_hi_dep; + auto* wlo = m_vars[w].m_lo_dep, *whi = m_vars[w].m_hi_dep; + if (s && lo(v) + 1 == 0 && is_fixed(v)) + return conflict(i, vlo, vhi), false; + if (s && lo(w) == 0 && is_fixed(w)) + return conflict(i, wlo, whi), false; + if (s && hi(w) != 0 && lo(w) <= hi(w) && lo(v) <= hi(v) && hi(w) - 1 <= lo(v)) + return conflict(i, vlo, wlo, whi), false; + if (s && hi(v) == 0 && lo(w) < hi(w) && hi(w) - 1 <= lo(v)) + return conflict(i, vlo, vhi, wlo, whi), false; + if (!s && hi(w) != 0 && lo(w) <= hi(w) && hi(w) <= lo(v) && lo(v) <= hi(v)) + return conflict(i, vlo, vhi, wlo, whi), false; + if (!s && hi(w) != 0 && lo(w) <= hi(w) && hi(w) <= lo(v) && hi(v) == 0) + return conflict(i, vlo, vhi, wlo, whi), false; + if (s && hi(w) != 0 && lo(w) <= hi(w) && hi(w) <= hi(v) && !new_bound(i, v, lo(v), hi(w) - 1, wlo, vhi, whi)) + return false; + if (s && lo(w) <= lo(v) && !new_bound(i, w, lo(v) + 1, hi(w), vlo, wlo)) + return false; + if (s && hi(w) != 0 && hi(w) - 1 <= lo(v) && lo(v) <= hi(v) && hi(w) < lo(w) && !new_bound(i, w, lo(w), 0, wlo, whi, vlo, vhi)) + return false; + if (s && hi(w) == 1 && !is_fixed(w) && !new_bound(i, w, lo(w), 0, wlo, whi)) + return false; + if (!s && hi(v) > hi(w) && !new_bound(i, v, lo(v), hi(w), vhi, whi)) + return false; + if (!s && lo(v) > lo(w) && !new_bound(i, w, lo(v), hi(w), vlo, wlo)) + return false; + if (!s && hi(w) != 0 && hi(w) < lo(w) && hi(w) <= lo(v) && lo(v) <= hi(v) && !new_bound(i, w, lo(w), 0, vlo, vhi, wlo, whi)) + return false; + if (hi(w) != 0 && lo(w) <= hi(w) && hi(w) <= lo(v) && !new_bound(i, v, 0, hi(v), wlo, vlo, whi)) + return false; + + return true; } template - lbool fixplex::conflict(ineq const& i, u_dependency* a, u_dependency* b) { - return conflict(a, m_deps.mk_join(mk_leaf(i.dep), b)); + void fixplex::conflict(ineq const& i, u_dependency* a, u_dependency* b, u_dependency* c, u_dependency* d) { + conflict(a, m_deps.mk_join(mk_leaf(i.dep), m_deps.mk_join(b, m_deps.mk_join(c, d)))); } template - lbool fixplex::conflict(u_dependency* a) { + void fixplex::conflict(u_dependency* a) { m_unsat_core.reset(); m_deps.linearize(a, m_unsat_core); - return l_false; + TRACE("polysat", tout << "core: " << m_unsat_core << "\n";); } template @@ -1213,31 +1243,31 @@ namespace polysat { } template - lbool fixplex::new_bound(ineq const& i, var_t x, numeral const& l, numeral const& h, u_dependency* a, u_dependency* b) { + bool fixplex::new_bound(ineq const& i, var_t x, numeral const& l, numeral const& h, u_dependency* a, u_dependency* b, u_dependency* c, u_dependency* d) { bool was_fixed = lo(x) + 1 == hi(x); - u_dependency* dep = m_deps.mk_join(mk_leaf(i.dep), m_deps.mk_join(a, b)); + u_dependency* dep = m_deps.mk_join(mk_leaf(i.dep), m_deps.mk_join(a, m_deps.mk_join(b, m_deps.mk_join(c, d)))); update_bounds(x, l, h, dep); if (m_vars[x].is_empty()) - return conflict(m_vars[x].m_lo_dep, m_vars[x].m_hi_dep); + return conflict(m_vars[x].m_lo_dep, m_vars[x].m_hi_dep), false; else if (!was_fixed && lo(x) + 1 == hi(x)) { // TBD: track based on inequality not row // fixed_var_eh(r, x); } - return l_true; + return true; } template - lbool fixplex::new_bound(row const& r, var_t x, mod_interval const& range) { + bool fixplex::new_bound(row const& r, var_t x, mod_interval const& range) { if (range.is_free()) return l_true; bool was_fixed = lo(x) + 1 == hi(x); update_bounds(x, range.lo, range.hi, row2dep(r)); IF_VERBOSE(0, verbose_stream() << "new-bound v" << x << " " << m_vars[x] << "\n"); - if (m_vars[x].is_empty()) - return conflict(m_vars[x].m_lo_dep, m_vars[x].m_hi_dep); + if (m_vars[x].is_empty()) + return conflict(m_vars[x].m_lo_dep, m_vars[x].m_hi_dep), false; else if (!was_fixed && lo(x) + 1 == hi(x)) fixed_var_eh(r, x); - return l_true; + return true; } template @@ -1249,6 +1279,12 @@ namespace polysat { if (vi.m_is_base) out << "b:" << vi.m_base2row << " " << pp(m_rows[vi.m_base2row].m_value) << " "; out << "\n"; } + for (auto const& i : m_ineqs) { + if (i.strict) + out << "v" << i.v << " < v" << i.w << "\n"; + else + out << "v" << i.v << " <= v" << i.w << "\n"; + } return out; } diff --git a/src/math/polysat/linear_solver.cpp b/src/math/polysat/linear_solver.cpp index 535d9c670..e2e127442 100644 --- a/src/math/polysat/linear_solver.cpp +++ b/src/math/polysat/linear_solver.cpp @@ -286,12 +286,10 @@ namespace polysat { lbool linear_solver::check() { lbool res = l_true; m_unsat_f = nullptr; - for (auto* f : m_fix) { - if (!f) - continue; - lbool r = f->make_feasible(); + for (auto* fp : m_fix_ptr) { + lbool r = fp->make_feasible(); if (r == l_false) { - m_unsat_f = f; + m_unsat_f = fp; return r; } if (r == l_undef) @@ -313,12 +311,13 @@ namespace polysat { } // current value assigned to (linear) variable according to tableau. - rational linear_solver::value(pvar v) { + bool linear_solver::value(pvar v, rational& val) { unsigned sz = s.size(v); auto* fp = sz2fixplex(sz); if (!fp) - return rational::zero(); - return fp->get_value(pvar2var(sz, v)); + return false; + val = fp->get_value(pvar2var(sz, v)); + return true; } }; diff --git a/src/math/polysat/linear_solver.h b/src/math/polysat/linear_solver.h index 88602a4e3..4ec9d855c 100644 --- a/src/math/polysat/linear_solver.h +++ b/src/math/polysat/linear_solver.h @@ -41,8 +41,6 @@ namespace polysat { inc_level_i, add_var_i, add_mono_i, - add_ineq_i, - add_row_i, set_active_i }; @@ -135,7 +133,7 @@ namespace polysat { void unsat_core(ptr_vector& constraints, unsigned_vector& deps); // current value assigned to (linear) variable according to tableau. - rational value(pvar v); + bool value(pvar v, rational& val); std::ostream& display(std::ostream& out) const { return out; } void collect_statistics(::statistics & st) const {} diff --git a/src/test/fixplex.cpp b/src/test/fixplex.cpp index d1d35670c..c24e1851d 100644 --- a/src/test/fixplex.cpp +++ b/src/test/fixplex.cpp @@ -1,5 +1,10 @@ #include "math/polysat/fixplex.h" #include "math/polysat/fixplex_def.h" +#include "ast/bv_decl_plugin.h" +#include "ast/reg_decl_plugins.h" +#include "smt/smt_kernel.h" +#include "smt/params/smt_params.h" + namespace polysat { @@ -118,7 +123,6 @@ namespace polysat { static void test_gcd() { std::cout << "gcd\n"; uint64_ext::manager e; - uint64_t x = 0, y = 0, z = 0, a = 0, b = 0; uint64_t g = e.gcd(15, 27); std::cout << g << "\n"; std::cout << 15 << " " << e.mul_inverse(15) << " " << 15*e.mul_inverse(15) << "\n"; @@ -126,9 +130,172 @@ namespace polysat { std::cout << 60 << " " << e.mul_inverse(60) << " " << 60*e.mul_inverse(60) << "\n"; std::cout << 29 << " " << e.mul_inverse(29) << " " << 29*e.mul_inverse(29) << "\n"; } + + static void test_ineq1() { + std::cout << "ineq1\n"; + scoped_fp fp; + var_t x = 0, y = 1; + fp.add_lt(x, y, 1); + fp.add_le(y, x, 2); + fp.run(); + } + + static void test_ineqs() { + var_t x = 0, y = 1; + unsigned nb = 6; + uint64_t bounds[6] = { 0, 1, 2, 10 , (uint64_t)-2, (uint64_t)-1 }; + ast_manager m; + reg_decl_plugins(m); + bv_util bv(m); + expr_ref X(m.mk_const("x", bv.mk_sort(64)), m); + expr_ref Y(m.mk_const("y", bv.mk_sort(64)), m); + smt_params pa; + smt::kernel solver(m, pa); + + auto mk_ult = [&](expr* a, expr* b) { + return m.mk_not(bv.mk_ule(b, a)); + }; + + auto add_bound = [&](expr* x, uint64_t i, uint64_t j) { + expr_ref I(bv.mk_numeral(i, 64), m); + expr_ref J(bv.mk_numeral(j, 64), m); + if (i < j) { + solver.assert_expr(bv.mk_ule(I, x)); + solver.assert_expr(mk_ult(x, J)); + } + else if (i > j && j != 0) + solver.assert_expr(m.mk_or(bv.mk_ule(I, x), mk_ult(x, J))); + else if (i > j && j == 0) + solver.assert_expr(bv.mk_ule(I, x)); + }; + + auto test_le = [&](bool test_le, uint64_t i, uint64_t j, uint64_t k, uint64_t l) { + if (i == j && i != 0) + return; + if (k == l && k != 0) + return; + scoped_fp fp; + fp.set_bounds(x, i, j, 1); + fp.set_bounds(y, k, l, 2); + if (test_le) + fp.add_le(x, y, 3); + else + fp.add_lt(x, y, 3); + + lbool feas = fp.make_feasible(); + + // validate result + + + solver.push(); + if (test_le) + solver.assert_expr(bv.mk_ule(X, Y)); + else + solver.assert_expr(mk_ult(X, Y)); + add_bound(X, i, j); + add_bound(Y, k, l); + + lbool feas2 = solver.check(); + + + if (feas == feas2) { + solver.pop(1); + return; + } + + + if (feas2 == l_false && feas == l_true) { + std::cout << "BUG!\n"; + solver.pop(1); + return; + } + + bool bad = false; + + switch (feas) { + case l_false: + for (unsigned c : fp.get_unsat_core()) + std::cout << c << "\n"; + std::cout << "\n"; + break; + case l_true: + case l_undef: + + // Check for missed bounds: + solver.push(); + solver.assert_expr(m.mk_eq(X, bv.mk_numeral(fp.lo(x), 64))); + if (l_true != solver.check()) { + std::cout << "missed lower bound on x\n"; + bad = true; + } + solver.pop(1); + + solver.push(); + solver.assert_expr(m.mk_eq(Y, bv.mk_numeral(fp.lo(y), 64))); + if (l_true != solver.check()) { + std::cout << "missed lower bound on y\n"; + bad = true; + } + solver.pop(1); + + if (fp.lo(x) != fp.hi(x) && fp.hi(x) != 0) { + solver.push(); + solver.assert_expr(m.mk_eq(X, bv.mk_numeral(fp.hi(x)-1, 64))); + if (l_true != solver.check()) { + std::cout << "missed upper bound on x\n"; + bad = true; + } + solver.pop(1); + } + + if (fp.lo(y) != fp.hi(y) && fp.hi(y) != 0) { + solver.push(); + solver.assert_expr(m.mk_eq(Y, bv.mk_numeral(fp.hi(y) - 1, 64))); + if (l_true != solver.check()) { + std::cout << "missed upper bound on y\n"; + bad = true; + } + solver.pop(1); + } + + if (bad) { + std::cout << fp << "\n"; + std::cout << feas << " " << feas2 << "\n"; + } + + break; + } + + // check assignment is valid when returning satisfiable. + if (false && feas == l_true) { + rational vx = fp.get_value(x); + rational vy = fp.get_value(y); + SASSERT(vx <= vy); + SASSERT(i >= j || (rational(i, rational::ui64()) <= vx && vx < rational(j, rational::ui64()))); + SASSERT(k >= l || (rational(k, rational::ui64()) <= vy && vy < rational(l, rational::ui64()))); + + } + + solver.pop(1); + + return; + }; + for (unsigned i = 0; i < nb; ++i) + for (unsigned j = 0; j < nb; ++j) + for (unsigned k = 0; k < nb; ++k) + for (unsigned l = 0; l < nb; ++l) { + test_le(true, bounds[i], bounds[j], bounds[k], bounds[l]); + test_le(false, bounds[i], bounds[j], bounds[k], bounds[l]); + } + } } void tst_fixplex() { + + polysat::test_ineq1(); + polysat::test_ineqs(); + return; + polysat::test1(); polysat::test2(); polysat::test3();