From be9f172cc062f6892deba4e40d7d5c8d218edab7 Mon Sep 17 00:00:00 2001 From: Nikolaj Bjorner Date: Wed, 4 Aug 2021 14:02:32 -0700 Subject: [PATCH] adding deps Signed-off-by: Nikolaj Bjorner --- src/math/bigfix/u256.h | 12 +- src/math/polysat/fixplex.h | 55 +++++--- src/math/polysat/fixplex_def.h | 218 +++++++++++++++++------------ src/math/polysat/linear_solver.cpp | 96 ++++++++----- src/math/polysat/linear_solver.h | 34 +++-- src/math/polysat/solver.cpp | 17 ++- src/test/fixplex.cpp | 30 ++-- 7 files changed, 268 insertions(+), 194 deletions(-) diff --git a/src/math/bigfix/u256.h b/src/math/bigfix/u256.h index d0addf591..41a146d48 100644 --- a/src/math/bigfix/u256.h +++ b/src/math/bigfix/u256.h @@ -70,12 +70,8 @@ inline std::ostream& operator<<(std::ostream& out, u256 const& u) { inline bool operator<(uint64_t n, u256 const& y) { return y > n; } inline bool operator<=(uint64_t n, u256 const& y) { return y >= n; } inline bool operator>(uint64_t n, u256 const& y) { return y < n; } -inline unsigned trailing_zeros(u256 const& n) { - NOT_IMPLEMENTED_YET(); - return 0; -} -inline u256 operator-(uint64_t n, u256 const& y) { - u256 x(n); - return x - y; -} +inline unsigned trailing_zeros(u256 const& n) { return n.trailing_zeros(); } + +inline u256 operator-(uint64_t n, u256 const& y) { return u256(n) - y; } inline bool operator>=(uint64_t n, u256 const& y) { return y <= n; } +inline rational to_rational(u256 const& x) { return x.to_rational(); } diff --git a/src/math/polysat/fixplex.h b/src/math/polysat/fixplex.h index 45ce69ba6..d0f94ff91 100644 --- a/src/math/polysat/fixplex.h +++ b/src/math/polysat/fixplex.h @@ -27,6 +27,8 @@ Author: #include "util/lbool.h" #include "util/uint_set.h" +inline rational to_rational(uint64_t n) { return rational(n, rational::ui64()); } + namespace polysat { typedef unsigned var_t; @@ -38,12 +40,14 @@ namespace polysat { virtual void del_row(var_t base_var) = 0; virtual std::ostream& display(std::ostream& out) const = 0; virtual void collect_statistics(::statistics & st) const = 0; - virtual void set_bounds(var_t v, rational const& lo, rational const& hi) = 0; - virtual void set_value(var_t v, rational const& val) = 0; + virtual void set_bounds(var_t v, rational const& lo, rational const& hi, unsigned dep) = 0; + virtual void set_value(var_t v, rational const& val, unsigned dep) = 0; + virtual rational get_value(var_t v) = 0; virtual void restore_bound() = 0; - virtual void add_le(var_t v, var_t w) = 0; - virtual void add_lt(var_t v, var_t w) = 0; + virtual void add_le(var_t v, var_t w, unsigned dep) = 0; + virtual void add_lt(var_t v, var_t w, unsigned dep) = 0; virtual void restore_ineq() = 0; + virtual void get_infeasible_deps(unsigned_vector& deps) = 0; }; @@ -92,7 +96,9 @@ namespace polysat { struct var_info : public mod_interval { unsigned m_base2row:29; unsigned m_is_base:1; - numeral m_value { 0 }; + numeral m_value = 0; + unsigned m_lo_dep = UINT_MAX; + unsigned m_hi_dep = UINT_MAX; var_info(): m_base2row(0), m_is_base(false) @@ -115,10 +121,10 @@ namespace polysat { numeral m_base_coeff; }; - struct stashed_bound : mod_interval { + struct stashed_bound : var_info { var_t m_var; - stashed_bound(var_t v, mod_interval const& i): - mod_interval(i), + stashed_bound(var_t v, var_info const& i): + var_info(i), m_var(v) {} }; @@ -131,12 +137,13 @@ namespace polysat { }; struct ineq { - var_t v { null_var }; - var_t w { null_var }; - bool strict { false }; - bool is_active { true }; - ineq(var_t v, var_t w, bool s): - v(v), w(w), strict(s) {} + var_t v = null_var; + var_t w = null_var; + bool strict = false; + bool is_active = true; + unsigned dep = UINT_MAX; + ineq(var_t v, var_t w, unsigned dep, bool s): + v(v), w(w), strict(s), dep(dep) {} }; static const var_t null_var = UINT_MAX; @@ -179,14 +186,15 @@ namespace polysat { std::ostream& display(std::ostream& out) const override; void collect_statistics(::statistics & st) const override; void del_row(var_t base_var) override; - void set_bounds(var_t v, rational const& lo, rational const& hi) override; - void set_value(var_t v, rational const& val) override; + void set_bounds(var_t v, rational const& lo, rational const& hi, unsigned dep) override; + void set_value(var_t v, rational const& val, unsigned dep) override; + rational get_value(var_t v) override; void restore_bound() override; - void add_le(var_t v, var_t w) override { add_ineq(v, w, false); } - void add_lt(var_t v, var_t w) override { add_ineq(v, w, true); } + void add_le(var_t v, var_t w, unsigned dep) override { add_ineq(v, w, dep, false); } + void add_lt(var_t v, var_t w, unsigned dep) override { add_ineq(v, w, dep, true); } virtual void restore_ineq() override; - void set_bounds(var_t v, numeral const& lo, numeral const& hi); + void set_bounds(var_t v, numeral const& lo, numeral const& hi, unsigned dep); void unset_bounds(var_t v) { m_vars[v].set_free(); } @@ -202,10 +210,13 @@ namespace polysat { void reset_eqs() { m_var_eqs.reset(); } void add_row(var_t base, unsigned num_vars, var_t const* vars, numeral const* coeffs); - row get_infeasible_row(); + + void get_infeasible_deps(unsigned_vector& deps) override; private: + row get_infeasible_row(); + std::ostream& display_row(std::ostream& out, row const& r, bool values = true); var_t get_base_var(row const& r) const { return m_rows[r.id()].m_base; } @@ -258,7 +269,7 @@ namespace polysat { var_t select_error_var(bool least); // facilities for handling inequalities - void add_ineq(var_t v, var_t w, bool strict); + void add_ineq(var_t v, var_t w, unsigned dep, bool strict); void touch_var(var_t x); lbool check_ineqs(); @@ -284,7 +295,6 @@ namespace polysat { numeral const& old_value_y); }; - template struct generic_uint_ext { typedef uint_type numeral; @@ -302,6 +312,7 @@ namespace polysat { } }; numeral from_rational(rational const& n) { return static_cast(n.get_uint64()); } + rational to_rational(numeral const& n) const { return ::to_rational(n); } void reset() {} void reset(numeral& n) { n = 0; } void del(numeral const& n) {} diff --git a/src/math/polysat/fixplex_def.h b/src/math/polysat/fixplex_def.h index f26f303b6..2ce807b58 100644 --- a/src/math/polysat/fixplex_def.h +++ b/src/math/polysat/fixplex_def.h @@ -70,16 +70,16 @@ namespace polysat { template fixplex::~fixplex() { reset(); - } - + } + template void fixplex::ensure_var(var_t v) { while (v >= m_vars.size()) { M.ensure_var(m_vars.size()); - m_vars.push_back(var_info()); + m_vars.push_back(var_info()); } - if (m_to_patch.get_bounds() <= v) - m_to_patch.set_bounds(2*v+1); + if (m_to_patch.get_bounds() <= v) + m_to_patch.set_bounds(2 * v + 1); } template @@ -105,7 +105,7 @@ namespace polysat { SASSERT(well_formed()); while ((v = select_var_to_fix()) != null_var) { TRACE("polysat", display(tout << "v" << v << "\n");); - if (!m_limit.inc() || num_iterations > m_max_iterations) + if (!m_limit.inc() || num_iterations > m_max_iterations) return l_undef; check_blands_rule(v, num_repeated); switch (make_var_feasible(v)) { @@ -138,20 +138,20 @@ namespace polysat { template void fixplex::add_row(var_t base_var, unsigned num_vars, var_t const* vars, numeral const* coeffs) { - for (unsigned i = 0; i < num_vars; ++i) + for (unsigned i = 0; i < num_vars; ++i) ensure_var(vars[i]); m_base_vars.reset(); row r = M.mk_row(); - for (unsigned i = 0; i < num_vars; ++i) - if (coeffs[i] != 0) + for (unsigned i = 0; i < num_vars; ++i) + if (coeffs[i] != 0) M.add_var(r, coeffs[i], vars[i]); numeral base_coeff = 0; numeral value = 0; for (auto const& e : M.row_entries(r)) { var_t v = e.var(); - if (v == base_var) + if (v == base_var) base_coeff = e.coeff(); else { if (is_base(v)) @@ -161,7 +161,7 @@ namespace polysat { } SASSERT(base_coeff != 0); SASSERT(!is_base(base_var)); - while (m_rows.size() <= r.id()) + while (m_rows.size() <= r.id()) m_rows.push_back(row_info()); m_rows[r.id()].m_base = base_var; m_rows[r.id()].m_base_coeff = base_coeff; @@ -189,9 +189,9 @@ namespace polysat { /** * Eliminate base variable v from all rows except where v is basic. * Return false if elimination required to multiply a non-basic row with an even number. - * This happens when the parity in the non-basic row is smaller than the parity of v in - * the basic row. It is expected to be a corner case and we don't try to solve this - * inside of fixplex. Instead to deal with the corner case we assume the layer around + * This happens when the parity in the non-basic row is smaller than the parity of v in + * the basic row. It is expected to be a corner case and we don't try to solve this + * inside of fixplex. Instead to deal with the corner case we assume the layer around * fixplex uses a solution from fixplex as a starting point for a complete search (in polysat). */ template @@ -248,14 +248,14 @@ namespace polysat { var_t old_base = row2base(r); numeral new_value; var_info& vi = m_vars[old_base]; - if (!vi.contains(value(old_base))) + if (!vi.contains(value(old_base))) new_value = vi.lo; - else + else new_value = value(old_base); // need to move var such that old_base comes in bound. - pivot(old_base, var, coeff, new_value); + pivot(old_base, var, coeff, new_value); SASSERT(is_base(var)); - SASSERT(base2row(var).id() == r.id()); + SASSERT(base2row(var).id() == r.id()); SASSERT(vi.contains(value(old_base))); } del_row(r); @@ -288,17 +288,17 @@ namespace polysat { ri.m_value += delta * c.get_row_entry().coeff(); set_base_value(s); add_patch(s); - } - } + } + } /** * Attempt to improve assigment to make x feasible. - * + * * return l_false if x is base variable of infeasible row * return l_true if it is possible to find an assignment that improves * return l_undef if the row could not be used for an improvement - * + * * delta - improvement over previous gap to feasible bound. * new_value - the new value to assign to x within its bounds. */ @@ -320,37 +320,37 @@ namespace polysat { return l_undef; } - + pivot(x, y, b, new_value); return l_true; } template - var_t fixplex::select_pivot(var_t x, numeral const& new_value, numeral & out_b) { + var_t fixplex::select_pivot(var_t x, numeral const& new_value, numeral& out_b) { if (m_bland) return select_pivot_blands(x, new_value, out_b); return select_pivot_core(x, new_value, out_b); } /** - \brief Select a variable y in the row r defining the base var x, + \brief Select a variable y in the row r defining the base var x, s.t. y can be used to patch the error in x_i. Return null_var if there is no y. Otherwise, return y and store its coefficient in out_b. The routine gives up if the coefficients of all free variables do not have the minimal - number of trailing zeros. + number of trailing zeros. */ template - var_t fixplex::select_pivot_core(var_t x, numeral const& new_value, numeral & out_b) { + var_t fixplex::select_pivot_core(var_t x, numeral const& new_value, numeral& out_b) { SASSERT(is_base(x)); - var_t max = get_num_vars(); + var_t max = get_num_vars(); var_t result = max; row r = base2row(x); int n = 0; unsigned best_col_sz = UINT_MAX; - int best_so_far = INT_MAX; + int best_so_far = INT_MAX; numeral a = row2base_coeff(r); numeral row_value = row2value(r) + a * new_value; numeral delta_y = 0; @@ -358,13 +358,13 @@ namespace polysat { bool best_in_bounds = false; for (auto const& r : M.row_entries(r)) { - var_t y = r.var(); - numeral const & b = r.coeff(); - if (x == y) + var_t y = r.var(); + numeral const& b = r.coeff(); + if (x == y) continue; if (!has_minimal_trailing_zeros(y, b)) continue; - numeral new_y_value = solve_for(row_value - b*value(y), b); + numeral new_y_value = solve_for(row_value - b * value(y), b); bool _in_bounds = in_bounds(y, new_y_value); if (!_in_bounds) { if (lo(y) - new_y_value < new_y_value - hi(y)) @@ -372,7 +372,7 @@ namespace polysat { else delta_y = new_y_value - hi(y) - 1; } - int num = get_num_non_free_dep_vars(y, best_so_far); + int num = get_num_non_free_dep_vars(y, best_so_far); unsigned col_sz = M.column_size(y); bool is_improvement = false, is_plateau = false; @@ -391,23 +391,23 @@ namespace polysat { is_plateau = true; else if (best_in_bounds && _in_bounds && best_so_far == num && col_sz == best_col_sz) is_plateau = true; - + if (is_improvement) { - result = y; - out_b = b; - best_so_far = num; - best_col_sz = col_sz; + result = y; + out_b = b; + best_so_far = num; + best_col_sz = col_sz; best_in_bounds = _in_bounds; - delta_best = delta_y; - n = 1; - } + delta_best = delta_y; + n = 1; + } else if (is_plateau) { n++; if (m_random() % n == 0) { - result = y; - out_b = b; + result = y; + out_b = b; } - } + } } if (result == max) return null_var; @@ -417,18 +417,18 @@ namespace polysat { } template - var_t fixplex::select_pivot_blands(var_t x, numeral const& new_value, numeral & out_b) { + var_t fixplex::select_pivot_blands(var_t x, numeral const& new_value, numeral& out_b) { SASSERT(is_base(x)); unsigned max = get_num_vars(); var_t result = max; row r = base2row(x); for (auto const& c : M.col_entries(r)) { - var_t y = c.var(); + var_t y = c.var(); if (x == y || y >= result) continue; - numeral const & b = c.coeff(); + numeral const& b = c.coeff(); if (can_improve(y, b)) { - out_b = b; + out_b = b; result = y; } } @@ -436,8 +436,8 @@ namespace polysat { } /** - * determine whether setting x := new_value - * allows to change the value of y in a direction + * determine whether setting x := new_value + * allows to change the value of y in a direction * that reduces or maintains the overall error. */ template @@ -452,30 +452,30 @@ namespace polysat { /** * Compute delta to add to the value, such that value + delta is either lo(v), or hi(v) - 1 - * A pre-condition is that value is not in the interval [lo(v),hi(v)[, + * A pre-condition is that value is not in the interval [lo(v),hi(v)[, * and therefore as a consequence lo(v) != hi(v). */ template - typename fixplex::numeral - fixplex::value2delta(var_t v, numeral const& value) const { + typename fixplex::numeral + fixplex::value2delta(var_t v, numeral const& value) const { SASSERT(!in_bounds(v)); SASSERT(lo(v) != hi(v)); if (lo(v) - value < value - hi(v)) return lo(v) - value; else - return hi(v) - value - 1; + return hi(v) - value - 1; } 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)); if (lo(v) - value < value - hi(v)) return lo(v) - value; else - return value - hi(v) - 1; + return value - hi(v) - 1; } /** @@ -486,46 +486,60 @@ namespace polysat { * - the variable v is queued to patch if v is basic. */ template - void fixplex::set_bounds(var_t v, numeral const& l, numeral const& h) { - m_vars[v] = mod_interval(l, h); + void fixplex::set_bounds(var_t v, numeral const& l, numeral const& h, unsigned dep) { + auto lo0 = m_vars[v].lo; + auto hi0 = m_vars[v].hi; + m_vars[v] &= mod_interval(l, h); + if (lo0 != m_vars[v].lo) + m_vars[v].m_lo_dep = dep; + if (hi0 != m_vars[v].hi) + m_vars[v].m_hi_dep = dep; if (in_bounds(v)) return; - if (is_base(v)) + if (is_base(v)) add_patch(v); else update_value(v, value2delta(v, value(v))); } template - void fixplex::set_bounds(var_t v, rational const& _lo, rational const& _hi) { + void fixplex::set_bounds(var_t v, rational const& _lo, rational const& _hi, unsigned dep) { numeral lo = m.from_rational(_lo); numeral hi = m.from_rational(_hi); m_stashed_bounds.push_back(stashed_bound(v, m_vars[v])); - set_bounds(v, lo, hi); + set_bounds(v, lo, hi, dep); } template - void fixplex::set_value(var_t v, rational const& _val) { + void fixplex::set_value(var_t v, rational const& _val, unsigned dep) { numeral val = m.from_rational(_val); m_stashed_bounds.push_back(stashed_bound(v, m_vars[v])); - set_bounds(v, val, val + 1); + set_bounds(v, val, val + 1, dep); + } + + template + rational fixplex::get_value(var_t v) { + return m.to_rational(m_vars[v].m_value); } template void fixplex::restore_bound() { auto const& b = m_stashed_bounds.back(); - set_bounds(b.m_var, b.lo, b.hi); + m_vars[b.m_var].lo = b.lo; + m_vars[b.m_var].hi = b.hi; + m_vars[b.m_var].m_lo_dep = b.m_lo_dep; + m_vars[b.m_var].m_hi_dep = b.m_hi_dep; m_stashed_bounds.pop_back(); } template - void fixplex::add_ineq(var_t v, var_t w, bool strict) { + void fixplex::add_ineq(var_t v, var_t w, unsigned dep, bool strict) { unsigned idx = m_ineqs.size(); m_var2ineqs.reserve(std::max(v, w) + 1); m_var2ineqs[v].push_back(idx); m_var2ineqs[w].push_back(idx); m_ineqs_to_check.push_back(idx); - m_ineqs.push_back(ineq(v, w, strict)); + m_ineqs.push_back(ineq(v, w, dep, strict)); } template @@ -577,7 +591,7 @@ namespace polysat { m_ineqs_to_check.reset(); m_vars_to_untouch.reset(); return l_true; - } + } /** @@ -603,15 +617,15 @@ namespace polysat { * A row is linear infeasible if it can be established * that none of the available assignments within current * bounds let the row add up to 0. - * + * * Assume the row is of the form: * ax + by + cz = 0 * with bounds * x : [lo_x, hi_x[, y : [lo_y, hi_y[, z : [lo_z : hi_z[ - * + * * Let range = [lo_x, hi_x[ + [lo_y, hi_y[ + [lo_z : hi_z[ * Claim. If range does not contain 0, then the row is infeasible. - * + * */ template bool fixplex::is_infeasible_row(var_t x) { @@ -624,7 +638,7 @@ namespace polysat { range += m_vars[v] * c; if (range.is_free()) return false; - } + } return !range.contains(0); } @@ -650,14 +664,14 @@ namespace polysat { var_t v = e.var(); auto c = e.coeff(); if (is_fixed(v)) - fixed += value(v)*c; - else + fixed += value(v) * c; + else parity = std::min(m.trailing_zeros(c), parity); } if (m.trailing_zeros(fixed) < parity) return true; - + return false; } @@ -676,23 +690,23 @@ namespace polysat { Effect: base(r_x) := y - value(x) := new_value + value(x) := new_value value(r_x) := value(r_x) - b*value(y) + a*new_value value(y) := -value(r_x) / b base_coeff(r_x) := b - + Let r be a row where y has coefficient c != 0. Assume: trailing_zeros(c) >= trailing_zeros(b) - + z = base(r) d = base_coeff(r) b1 = (b >> tz(b)) - c1 = (c >> (tz(c) - tz(b))) + c1 = (c >> (tz(c) - tz(b))) r <- b1 * r - c1 * r_x value(r) := b1 * value(r) - b1 * old_value(y) - c1 * value(r_x) value(z) := - value(r) / d base_coeff(r) := b1 * base_coeff(r) - */ + */ template void fixplex::pivot(var_t x, var_t y, numeral const& b, numeral const& new_value) { ++m_stats.m_num_pivots; @@ -705,7 +719,7 @@ namespace polysat { numeral const& a = row_x.m_base_coeff; numeral old_value_y = yI.m_value; row_x.m_base = y; - row_x.m_value = row_x.m_value - b*old_value_y + a*new_value; + row_x.m_value = row_x.m_value - b * old_value_y + a * new_value; row_x.m_base_coeff = b; yI.m_base2row = rx; yI.m_is_base = true; @@ -718,7 +732,7 @@ namespace polysat { SASSERT(well_formed_row(r_x)); unsigned tz_b = m.trailing_zeros(b); - + for (auto col : M.col_entries(y)) { row r_z = col.get_row(); unsigned rz = r_z.id(); @@ -740,12 +754,12 @@ namespace polysat { * * returns true if elimination preserves equivalence (is lossless). */ - template + template bool fixplex::eliminate_var( - row const& r_y, - row const& r_z, - numeral const& c, - unsigned tz_b, + row const& r_y, + row const& r_z, + numeral const& c, + unsigned tz_b, numeral const& old_value_y) { numeral b = row2base_coeff(r_y); @@ -770,7 +784,7 @@ namespace polysat { return tz_b <= tz_c; } - template + template bool fixplex::is_feasible() const { for (unsigned i = m_vars.size(); i-- > 0; ) if (!in_bounds(i)) @@ -779,12 +793,34 @@ namespace polysat { } template - typename fixplex::row - fixplex::get_infeasible_row() { + typename fixplex::row + fixplex::get_infeasible_row() { SASSERT(is_base(m_infeasible_var)); return base2row(m_infeasible_var); } + /*** + * Extract dependencies for infeasible row. + * A safe approximation is to extract dependencies for all bounds. + * + * Different modes of infeasibility may not be based on a row: + * - inequalities + * - parity constraints between two rows. + */ + template + void fixplex::get_infeasible_deps(unsigned_vector& deps) { + auto row = get_infeasible_row(); + for (auto const& e : M.row_entries(row)) { + var_t v = e.var(); + auto lo = m_vars[v].m_lo_dep; + auto hi = m_vars[v].m_hi_dep; + if (lo != UINT_MAX) + deps.push_back(lo); + if (hi != UINT_MAX) + deps.push_back(hi); + } + } + /** \brief Return the number of base variables that are non free and are v dependent. The function adds 1 to the result if v is non free. diff --git a/src/math/polysat/linear_solver.cpp b/src/math/polysat/linear_solver.cpp index 420c7d0ae..2549b5586 100644 --- a/src/math/polysat/linear_solver.cpp +++ b/src/math/polysat/linear_solver.cpp @@ -60,11 +60,16 @@ namespace polysat { m_rows.pop_back(); break; } + case trail_i::set_active_i: + m_active.pop_back(); + break; + default: UNREACHABLE(); } m_trail.pop_back(); } + m_unsat_f = nullptr; } fixplex_base& linear_solver::sz2fixplex(unsigned sz) { @@ -78,10 +83,10 @@ namespace polysat { b = alloc(fixplex, s.m_lim); break; case 128: - NOT_IMPLEMENTED_YET(); + NOT_IMPLEMENTED_YET(); break; case 256: - b = alloc(fixplex>, s.m_lim); + b = alloc(fixplex>, s.m_lim); break; default: NOT_IMPLEMENTED_YET(); @@ -117,8 +122,16 @@ namespace polysat { auto pr = std::make_pair(v, v); m_bool_var2row.setx(c.bvar(), pr, pr); } + + void linear_solver::new_le(ule_constraint& c) { + var_t v = internalize_pdd(c.lhs()); + var_t w = internalize_pdd(c.rhs()); + auto pr = std::make_pair(v, w); + m_bool_var2row.setx(c.bvar(), pr, pr); + } - void linear_solver::assert_eq(eq_constraint& c) { + void linear_solver::assert_eq(bool is_positive, eq_constraint& c) { + unsigned c_dep = constraint_idx2dep(m_active.size() - 1); var_t v = m_bool_var2row[c.bvar()].first; unsigned sz = c.p().power_of_2(); auto& fp = sz2fixplex(sz); @@ -126,17 +139,10 @@ namespace polysat { m_rows.push_back(std::make_pair(v, sz)); rational z(0), o(1); SASSERT(!c.is_undef()); - if (c.is_positive()) - fp.set_bounds(v, z, z); + if (is_positive) + fp.set_bounds(v, z, z, c_dep); else - fp.set_bounds(v, o, z); -} - - void linear_solver::new_le(ule_constraint& c) { - var_t v = internalize_pdd(c.lhs()); - var_t w = internalize_pdd(c.rhs()); - auto pr = std::make_pair(v, w); - m_bool_var2row.setx(c.bvar(), pr, pr); + fp.set_bounds(v, o, z, c_dep); } // @@ -152,44 +158,45 @@ namespace polysat { // inequality graph (with offsets) // - void linear_solver::assert_le(ule_constraint& c) { + void linear_solver::assert_le(bool is_positive, ule_constraint& c) { auto [v, w] = m_bool_var2row[c.bvar()]; unsigned sz = c.lhs().power_of_2(); auto& fp = sz2fixplex(sz); + unsigned c_dep = constraint_idx2dep(m_active.size() - 1); rational z(0); if (c.rhs().is_val()) { bool is_max_value = false; - if (c.is_positive()) + if (is_positive) // v <= rhs - fp.set_bounds(v, z, c.rhs().val()); + fp.set_bounds(v, z, c.rhs().val(), c_dep); else if (is_max_value) throw default_exception("conflict not implemented"); else // rhs < v - fp.set_bounds(v, c.rhs().val() + 1, z); + fp.set_bounds(v, c.rhs().val() + 1, z, c_dep); m_trail.push_back(trail_i::set_bound_i); m_rows.push_back(std::make_pair(v, sz)); return; } if (c.lhs().is_val()) { - if (c.is_positive()) + if (is_positive) // w >= lhs - fp.set_bounds(w, c.lhs().val(), z); + fp.set_bounds(w, c.lhs().val(), z, c_dep); else if (c.lhs().val() == 0) throw default_exception("conflict not implemented"); else // w < lhs - fp.set_bounds(w, z, c.lhs().val() - 1); + fp.set_bounds(w, z, c.lhs().val() - 1, c_dep); m_trail.push_back(trail_i::set_bound_i); m_rows.push_back(std::make_pair(w, sz)); return; } - if (c.is_positive()) - fp.add_le(v, w); + if (is_positive) + fp.add_le(v, w, c_dep); else - fp.add_lt(w, v); + fp.add_lt(w, v, c_dep); m_trail.push_back(trail_i::add_ineq_i); m_rows.push_back(std::make_pair(v, sz)); } @@ -209,14 +216,16 @@ namespace polysat { } } - void linear_solver::activate_constraint(constraint& c) { + void linear_solver::activate_constraint(bool is_positive, constraint& c) { + m_active.push_back(&c); + m_trail.push_back(trail_i::set_active_i); SASSERT(!c.is_undef()); switch (c.kind()) { case ckind_t::eq_t: - assert_eq(c.to_eq()); + assert_eq(is_positive, c.to_eq()); break; case ckind_t::ule_t: - assert_le(c.to_ule()); + assert_le(is_positive, c.to_ule()); break; default: UNREACHABLE(); @@ -261,34 +270,37 @@ namespace polysat { return m_sz2num_vars[sz]++; } - void linear_solver::set_value(pvar v, rational const& value) { + void linear_solver::set_value(pvar v, rational const& value, unsigned dep) { unsigned sz = s.size(v); auto& fp = sz2fixplex(sz); var_t w = pvar2var(sz, v); m_trail.push_back(trail_i::set_bound_i); m_rows.push_back(std::make_pair(w, sz)); - fp.set_value(w, value); + fp.set_value(w, value, external_dep2dep(dep)); } - void linear_solver::set_bound(pvar v, rational const& lo, rational const& hi) { + void linear_solver::set_bound(pvar v, rational const& lo, rational const& hi, unsigned dep) { unsigned sz = s.size(v); auto& fp = sz2fixplex(sz); var_t w = pvar2var(sz, v); m_trail.push_back(trail_i::set_bound_i); m_rows.push_back(std::make_pair(w, sz)); - fp.set_bounds(w, lo, hi); + fp.set_bounds(w, lo, hi, external_dep2dep(dep)); } - // check integer modular feasibility under current bounds. - // and inequalities + /** + * check integer modular feasibility under current bounds. + * and inequalities + */ 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(); if (r == l_false) { - // m_unsat_f = f; + m_unsat_f = f; return r; } if (r == l_undef) @@ -297,13 +309,25 @@ namespace polysat { return res; } - void linear_solver::unsat_core(unsigned_vector& deps) { - NOT_IMPLEMENTED_YET(); + void linear_solver::unsat_core(ptr_vector& cs, unsigned_vector& deps) { + SASSERT(m_unsat_f); + deps.reset(); + cs.reset(); + m_unsat_f->get_infeasible_deps(deps); + unsigned j = 0; + for (unsigned dep : deps) { + if (is_constraint_dep(dep)) + cs.push_back(m_active[dep2constraint_idx(dep)]); + else + deps[j++] = dep2external_dep(dep); + } + deps.shrink(j); } // current value assigned to (linear) variable according to tableau. rational linear_solver::value(pvar v) { - return rational(0); + unsigned sz = s.size(v); + return sz2fixplex(sz).get_value(pvar2var(sz, v)); } }; diff --git a/src/math/polysat/linear_solver.h b/src/math/polysat/linear_solver.h index 8de262f91..51bc05fc0 100644 --- a/src/math/polysat/linear_solver.h +++ b/src/math/polysat/linear_solver.h @@ -43,7 +43,8 @@ namespace polysat { add_mono_i, set_bound_i, add_ineq_i, - add_row_i + add_row_i, + set_active_i }; struct mono_info { @@ -81,6 +82,7 @@ namespace polysat { svector> m_rows; unsigned_vector m_var2ext; unsigned_vector m_ext2var; + ptr_vector m_active; svector m_vars; vector m_coeffs; @@ -90,29 +92,31 @@ namespace polysat { unsigned_vector m_sz2num_vars; small_object_allocator m_alloc; svector m_monomials; + fixplex_base* m_unsat_f = nullptr; fixplex_base& sz2fixplex(unsigned sz); void linearize(pdd const& p); var_t fresh_var(unsigned sz); - var_t internalize_pdd(pdd const& p); void new_eq(eq_constraint& eq); void new_le(ule_constraint& le); - void assert_eq(eq_constraint& eq); - void assert_le(ule_constraint& le); + void assert_eq(bool is_positive, eq_constraint& eq); + void assert_le(bool is_positive, ule_constraint& le); // bind monomial to variable. var_t mono2var(unsigned sz, unsigned_vector const& m); var_t pvar2var(unsigned sz, pvar v); - unsigned_vector var2mono(unsigned sz, var_t v) { throw default_exception("nyi"); } - // - // TBD trail object for - // removing variables - // undoing variable bounds bounds - // removing rows from fixplex - // + // unsigned_vector var2mono(unsigned sz, var_t v) { throw default_exception("nyi"); } + + // distinguish constraint and justification dependencies + unsigned external_dep2dep(unsigned dep) const { return UINT_MAX - dep; } + unsigned constraint_idx2dep(unsigned idx) const { return idx; } + bool is_constraint_dep(unsigned dep) const { return dep < UINT_MAX / 2; } + unsigned dep2constraint_idx(unsigned dep) const { return dep; } + unsigned dep2external_dep(unsigned dep) const { return UINT_MAX - dep; } + public: linear_solver(solver& s): s(s), @@ -122,13 +126,13 @@ namespace polysat { void push(); void pop(unsigned n); void new_constraint(constraint& c); - void set_value(pvar v, rational const& value); - void set_bound(pvar v, rational const& lo, rational const& hi); - void activate_constraint(constraint& c); + void set_value(pvar v, rational const& value, unsigned dep); + void set_bound(pvar v, rational const& lo, rational const& hi, unsigned dep); + void activate_constraint(bool is_positive, constraint& c); // check integer modular feasibility under current bounds. lbool check(); - void unsat_core(unsigned_vector& deps); + void unsat_core(ptr_vector& constraints, unsigned_vector& deps); // current value assigned to (linear) variable according to tableau. rational value(pvar v); diff --git a/src/math/polysat/solver.cpp b/src/math/polysat/solver.cpp index 0633ecd9d..e5c6a31eb 100644 --- a/src/math/polysat/solver.cpp +++ b/src/math/polysat/solver.cpp @@ -383,7 +383,8 @@ namespace polysat { m_trail.push_back(trail_instr_t::assign_i); m_justification[v] = j; #if ENABLE_LINEAR_SOLVER - m_linear_solver.set_value(v, val); + // TODO: convert justification into a format that can be tracked in a depdendency core. + m_linear_solver.set_value(v, val, UINT_MAX); #endif } @@ -942,7 +943,12 @@ namespace polysat { m_bvars.assign(lit, m_level, reason, lemma); } - /// Activate constraint immediately + /** + * Activate constraint immediately + * Activation and de-activation of constraints follows the scope controlled by push/pop. + * constraints activated within the linear solver are de-activated when the linear + * solver is popped. + */ void solver::activate_constraint(constraint& c, bool is_true) { LOG("Activating constraint: " << c << " ; is_true = " << is_true); SASSERT(m_bvars.value(c.bvar()) == to_lbool(is_true)); @@ -950,18 +956,15 @@ namespace polysat { add_watch(c); c.narrow(*this); #if ENABLE_LINEAR_SOLVER - m_linear_solver.activate_constraint(c); + m_linear_solver.activate_constraint(c.is_positive(), c); #endif } - /// Deactivate constraint immediately + /// Deactivate constraint void solver::deactivate_constraint(constraint& c) { LOG("Deactivating constraint: " << c); erase_watch(c); c.unassign(); -#if ENABLE_LINEAR_SOLVER - m_linear_solver.deactivate_constraint(c); // TODO add this method to linear solver? -#endif } void solver::backjump(unsigned new_level) { diff --git a/src/test/fixplex.cpp b/src/test/fixplex.cpp index 4fefbfcbd..d1d35670c 100644 --- a/src/test/fixplex.cpp +++ b/src/test/fixplex.cpp @@ -28,7 +28,7 @@ namespace polysat { var_t ys[3] = { x, y, z }; numeral coeffs[3] = { 2, 1, 4 }; fp.add_row(x, 3, ys, coeffs); - fp.set_bounds(x, 1, 2); + fp.set_bounds(x, 1, 2, 0); fp.run(); } @@ -38,7 +38,7 @@ namespace polysat { var_t ys[3] = { x, y, z }; numeral coeffs[3] = { 1, 2, 4 }; fp.add_row(x, 3, ys, coeffs); - fp.set_bounds(x, 3, 4); + fp.set_bounds(x, 3, 4, 1); fp.run(); } @@ -48,9 +48,9 @@ namespace polysat { var_t ys[3] = { x, y, z }; numeral coeffs[3] = { 1, 1, 1 }; fp.add_row(x, 3, ys, coeffs); - fp.set_bounds(x, 3, 4); - fp.set_bounds(y, 3, 4); - fp.set_bounds(z, 3, 4); + fp.set_bounds(x, 3, 4, 1); + fp.set_bounds(y, 3, 4, 2); + fp.set_bounds(z, 3, 4, 3); fp.run(); } @@ -60,22 +60,22 @@ namespace polysat { var_t ys[3] = { x, y, z }; numeral coeffs[3] = { 1, 1, 1 }; fp.add_row(x, 3, ys, coeffs); - fp.set_bounds(x, 3, 4); - fp.set_bounds(y, 3, 6); + fp.set_bounds(x, 3, 4, 1); + fp.set_bounds(y, 3, 6, 2); fp.run(); fp.propagate_bounds(); fp.reset(); coeffs[2] = 0ull - 1; fp.add_row(x, 3, ys, coeffs); - fp.set_bounds(x, 3, 4); - fp.set_bounds(y, 3, 6); + fp.set_bounds(x, 3, 4, 1); + fp.set_bounds(y, 3, 6, 2); fp.run(); fp.propagate_bounds(); fp.reset(); fp.add_row(x, 3, ys, coeffs); - fp.set_bounds(x, 3, 4); - fp.set_bounds(y, 3, 6); - fp.set_bounds(z, 1, 8); + fp.set_bounds(x, 3, 4, 1); + fp.set_bounds(y, 3, 6, 2); + fp.set_bounds(z, 1, 8, 3); fp.run(); fp.propagate_bounds(); fp.reset(); @@ -88,8 +88,8 @@ namespace polysat { var_t ys[3] = { x, y, z }; numeral coeffs[3] = { 1, 1, 1 }; fp.add_row(x, 3, ys, coeffs); - fp.set_bounds(x, 3, 4); - fp.set_bounds(y, 3, 6); + fp.set_bounds(x, 3, 4, 1); + fp.set_bounds(y, 3, 6, 2); var_t ys2[3] = { u, y, z }; fp.add_row(u, 3, ys2, coeffs); fp.run(); @@ -107,7 +107,7 @@ namespace polysat { numeral coeffs2[3] = { 1, m1, 1 }; fp.add_row(x, 3, ys1, coeffs1); fp.add_row(z, 3, ys2, coeffs2); - fp.set_bounds(u, 1, 2); + fp.set_bounds(u, 1, 2, 1); fp.run(); fp.propagate_eqs(); for (auto e : fp.var_eqs())