diff --git a/src/math/polysat/CMakeLists.txt b/src/math/polysat/CMakeLists.txt index f1a28ef58..ddc5e7edc 100644 --- a/src/math/polysat/CMakeLists.txt +++ b/src/math/polysat/CMakeLists.txt @@ -15,7 +15,6 @@ z3_add_component(polysat solver.cpp ule_constraint.cpp viable.cpp - viable2.cpp variable_elimination.cpp COMPONENT_DEPENDENCIES util diff --git a/src/math/polysat/conflict.cpp b/src/math/polysat/conflict.cpp index 267bfb0c6..081c5767d 100644 --- a/src/math/polysat/conflict.cpp +++ b/src/math/polysat/conflict.cpp @@ -100,12 +100,6 @@ namespace polysat { LOG("Conflict: v" << v); SASSERT(empty()); m_conflict_var = v; -#if !NEW_VIABLE - for (auto c : s.m_cjust[v]) { - c->set_var_dependent(); // ?? - insert(c); - } -#endif SASSERT(!empty()); } @@ -300,23 +294,13 @@ namespace polysat { return false; } -#if NEW_VIABLE if (conflict_var() == v && s.m_viable.resolve(v, *this)) return true; -#else - if (conflict_var() == v && s.m_forbidden_intervals.perform(v, cjust_v, *this)) - return true; -#endif m_vars.remove(v); -#if NEW_VIABLE for (auto const& c : s.m_viable.get_constraints(v)) insert(c); -#else - for (auto c : cjust_v) - insert(c); -#endif for (auto* engine : ex_engines) if (engine->try_explain(v, *this)) diff --git a/src/math/polysat/forbidden_intervals.cpp b/src/math/polysat/forbidden_intervals.cpp index 70945e030..a06089b07 100644 --- a/src/math/polysat/forbidden_intervals.cpp +++ b/src/math/polysat/forbidden_intervals.cpp @@ -22,147 +22,6 @@ Author: namespace polysat { - /** - * Find a sequence of intervals that covers all of Z_modulus. - * - * \returns true iff such a covering exists - * \param longest_i: the longest interval (as index into 'records') - * \param out_seq: will contain the covering (as list of indices into 'records') - */ - static bool find_covering_sequence(vector const& records, unsigned longest_i, rational modulus, unsigned_vector& out_seq) { - rational baseline = records[longest_i].interval.hi_val(); - while (!records[longest_i].interval.currently_contains(baseline)) { - rational best_extent = rational::zero(); - unsigned furthest_i = UINT_MAX; - for (unsigned i = records.size(); i-- > 0; ) { - auto const& interval = records[i].interval; - if (interval.currently_contains(baseline)) { - rational extent = mod(interval.hi_val() - baseline, modulus); - if (extent > best_extent) { - best_extent = extent; - furthest_i = i; - } - } - } - if (furthest_i == UINT_MAX) { - // There's a hole we can't cover. - // This can happen if a constraint didn't produce an interval - // (but not necessarily, values may be covered by multiple constraints) - return false; - } - SASSERT(best_extent > 0); - out_seq.push_back(furthest_i); - baseline = records[furthest_i].interval.hi_val(); - } - SASSERT(out_seq.size() > 0); - if (!records[out_seq[0]].interval.currently_contains(baseline)) - out_seq.push_back(longest_i); - return true; - } - - /** - * A single constraint implies an empty interval. - * We assume that neg_cond is a consequence of src that - * does not mention the variable v to be eliminated. - */ - void forbidden_intervals::full_interval_conflict(signed_constraint src, vector const& side_cond, conflict& core) { - core.reset(); - for (auto c : side_cond) - core.insert(c); - core.insert(src); - core.set_bailout(); - } - - bool forbidden_intervals::perform(pvar v, vector const& just, conflict& core) { - - // Extract forbidden intervals from conflicting constraints - vector records; - rational longest_len; - unsigned longest_i = UINT_MAX; - for (signed_constraint c : just) { - LOG_H3("Computing forbidden interval for: " << c); - eval_interval interval = eval_interval::full(); - vector side_cond; - if (get_interval(c, v, interval, side_cond)) { - LOG("interval: " << interval); - LOG("neg_cond: " << side_cond); - if (interval.is_currently_empty()) - continue; - if (interval.is_full()) { - // We have a single interval covering the whole domain - // => the side conditions of that interval are enough to produce a conflict - full_interval_conflict(c, side_cond, core); - revert_core(core); - return true; - } - else { - auto const len = interval.current_len(); - if (len > longest_len) { - longest_len = len; - longest_i = records.size(); - } - } - records.push_back({ std::move(interval), std::move(side_cond), c }); - } - } - - if (records.empty()) { - LOG("aborted (no intervals)"); - return false; - } - - - SASSERT(longest_i != UINT_MAX); - LOG("longest: i=" << longest_i << "; " << records[longest_i].interval); - - rational const modulus = rational::power_of_two(s.size(v)); - - // Select a sequence of covering intervals - unsigned_vector seq; - if (!find_covering_sequence(records, longest_i, modulus, seq)) { - LOG("aborted (intervals do not cover domain)"); - return false; - } - LOG("seq: " << seq); - SASSERT(seq.size() >= 2); // otherwise has_full should have been true - - // Update the conflict state - // Idea: - // - If the src constraints hold, and - // - if the side conditions hold, and - // - the upper bound of each interval is contained in the next interval, - // then the forbidden intervals cover the whole domain and we have a conflict. - // - - core.reset(); - // Add side conditions and interval constraints - for (unsigned seq_i = seq.size(); seq_i-- > 0; ) { - unsigned const i = seq[seq_i]; - unsigned const next_i = seq[(seq_i + 1) % seq.size()]; - // Build constraint: upper bound of each interval is not contained in the next interval, - // using the equivalence: t \in [l;h[ <=> t-l < h-l - auto const& hi = records[i].interval.hi(); - auto const& next_lo = records[next_i].interval.lo(); - auto const& next_hi = records[next_i].interval.hi(); - auto lhs = hi - next_lo; - auto rhs = next_hi - next_lo; - signed_constraint c = s.m_constraints.ult(lhs, rhs); - LOG("constraint: " << c); - core.insert(c); - // Side conditions - // TODO: check whether the condition is subsumed by c? maybe at the end do a "lemma reduction" step, to try and reduce branching? - for (auto sc : records[i].side_cond) - core.insert(sc); - - core.insert(records[i].src); - } - core.set_bailout(); - revert_core(core); - return true; - } - - - /** Precondition: all variables other than v are assigned. * * \param[out] out_interval The forbidden interval for this constraint @@ -464,15 +323,4 @@ namespace polysat { } return false; } - - void forbidden_intervals::revert_core(conflict& core) { - for (auto c : core) { - if (c.bvalue(s) == l_false) { - core.reset(); - core.set(~c); - return; - } - } - } - } diff --git a/src/math/polysat/forbidden_intervals.h b/src/math/polysat/forbidden_intervals.h index d8e067c35..31ba21bb4 100644 --- a/src/math/polysat/forbidden_intervals.h +++ b/src/math/polysat/forbidden_intervals.h @@ -22,8 +22,6 @@ namespace polysat { class forbidden_intervals { solver& s; - void revert_core(conflict& core); - void full_interval_conflict(signed_constraint c, vector const & side_cond, conflict& core); void push_eq(bool is_trivial, pdd const& p, vector& side_cond); eval_interval to_interval(signed_constraint const& c, bool is_trivial, rational const& coeff, @@ -60,7 +58,6 @@ namespace polysat { bool coefficient_is_01(dd::pdd_manager& m, rational const& r) { return r.is_zero() || r.is_one() || r == m.max_value(); }; public: forbidden_intervals(solver& s) :s(s) {} - bool perform(pvar v, vector const& just, conflict& core); bool get_interval(signed_constraint const& c, pvar v, eval_interval& out_interval, vector& side_cond); }; } diff --git a/src/math/polysat/saturation.cpp b/src/math/polysat/saturation.cpp index f0d11719f..23b4c1f8e 100644 --- a/src/math/polysat/saturation.cpp +++ b/src/math/polysat/saturation.cpp @@ -180,17 +180,10 @@ namespace polysat { if (x_max * y_max > pddm.max_value()) push_omega_bisect(x, x_max, y, y_max); else { -#if NEW_VIABLE for (auto const& c : s.m_viable.get_constraints(y.var())) m_new_constraints.insert(c); for (auto const& c : s.m_viable.get_constraints(x.var())) m_new_constraints.insert(c); -#else - for (auto c : s.m_cjust[y.var()]) - m_new_constraints.insert(c); - for (auto c : s.m_cjust[x.var()]) - m_new_constraints.insert(c); -#endif } } diff --git a/src/math/polysat/solver.cpp b/src/math/polysat/solver.cpp index 3c6a5e1d4..1d09a600c 100644 --- a/src/math/polysat/solver.cpp +++ b/src/math/polysat/solver.cpp @@ -98,9 +98,6 @@ namespace polysat { m_value.push_back(rational::zero()); m_justification.push_back(justification::unassigned()); m_viable.push(sz); -#if !NEW_VIABLE - m_cjust.push_back({}); -#endif m_pwatch.push_back({}); m_activity.push_back(0); m_vars.push_back(sz2pdd(sz).mk_var(v)); @@ -118,9 +115,6 @@ namespace polysat { // TODO also remove v from all learned constraints. pvar v = m_value.size() - 1; m_viable.pop(); -#if !NEW_VIABLE - m_cjust.pop_back(); -#endif m_value.pop_back(); m_justification.pop_back(); m_pwatch.pop_back(); @@ -345,15 +339,6 @@ namespace polysat { m_search.pop(); break; } - case trail_instr_t::just_i: { -#if !NEW_VIABLE - auto v = m_cjust_trail.back(); - LOG_V("Undo just_i"); - m_cjust[v].pop_back(); - m_cjust_trail.pop_back(); -#endif - break; - } default: UNREACHABLE(); } @@ -480,12 +465,8 @@ namespace polysat { LOG_H2("Resolve conflict"); LOG("\n" << *this); LOG("search state: " << m_search); -#if NEW_VIABLE - m_viable.log(); -#else - for (pvar v = 0; v < m_cjust.size(); ++v) - LOG("cjust[v" << v << "]: " << m_cjust[v]); -#endif + for (pvar v = 0; v < m_justification.size(); ++v) + LOG("v" << v << " " << viable::var_pp(m_viable, v)); ++m_stats.m_num_conflicts; SASSERT(is_conflict()); @@ -829,10 +810,6 @@ namespace polysat { pvar v = item.var(); auto const& j = m_justification[v]; out << "\t" << assignment_pp(*this, v, get_value(v)) << " @" << j.level(); -#if !NEW_VIABLE - if (j.is_propagation()) - out << " " << m_cjust[v]; -#endif out << "\n"; } else { diff --git a/src/math/polysat/solver.h b/src/math/polysat/solver.h index 3917bfd84..b0afdf01f 100644 --- a/src/math/polysat/solver.h +++ b/src/math/polysat/solver.h @@ -32,11 +32,9 @@ Author: #include "math/polysat/forbidden_intervals.h" #include "math/polysat/trail.h" #include "math/polysat/viable.h" -#include "math/polysat/viable2.h" #include "math/polysat/log.h" -#define NEW_VIABLE 1 namespace polysat { @@ -77,11 +75,7 @@ namespace polysat { params_ref m_params; scoped_ptr_vector m_pdd; -#if NEW_VIABLE - viable2 m_viable; -#else viable m_viable; // viable sets per variable -#endif linear_solver m_linear_solver; conflict m_conflict; forbidden_intervals m_forbidden_intervals; @@ -101,9 +95,6 @@ namespace polysat { // Per variable information vector m_value; // assigned value vector m_justification; // justification for variable assignment -#if !NEW_VIABLE - vector m_cjust; // constraints justifying variable range. -#endif vector m_pwatch; // watch list datastructure into constraints. unsigned_vector m_activity; @@ -132,17 +123,6 @@ namespace polysat { m_qhead_trail.pop_back(); } -#if !NEW_VIABLE - void push_cjust(pvar v, signed_constraint c) { - if (m_cjust[v].contains(c)) // TODO: better check (flag on constraint?) - return; - LOG_V("cjust[v" << v << "] += " << c); - SASSERT(c); - m_cjust[v].push_back(c); - m_trail.push_back(trail_instr_t::just_i); - m_cjust_trail.push_back(v); - } -#endif unsigned size(pvar v) const { return m_size[v]; } diff --git a/src/math/polysat/ule_constraint.cpp b/src/math/polysat/ule_constraint.cpp index 3c72657bc..08ed8f96d 100644 --- a/src/math/polysat/ule_constraint.cpp +++ b/src/math/polysat/ule_constraint.cpp @@ -129,9 +129,6 @@ namespace polysat { return; } -#if NEW_VIABLE - // setup with viable2: - // we no longer need cjust pvar v = null_var; if (p.is_unilinear()) v = p.var(); @@ -152,64 +149,6 @@ namespace polysat { break; } } -#else - - // p <= 0, e.g., p == 0 - if (q.is_zero() && p.is_unilinear()) { - // a*x + b == 0 - pvar v = p.var(); - s.push_cjust(v, { this, is_positive }); - - rational a = p.hi().val(); - rational b = p.lo().val(); - s.m_viable.intersect_eq(a, v, b, is_positive); - - rational val; - if (s.m_viable.find_viable(v, val) == dd::find_t::singleton) - s.propagate(v, val, { this, is_positive }); - return; - } - - pvar v = null_var; - rational a, b, c, d; - if (p.is_unilinear() && q.is_unilinear() && p.var() == q.var()) { - // a*x + b <=u c*x + d - v = p.var(); - a = p.hi().val(); - b = p.lo().val(); - c = q.hi().val(); - d = q.lo().val(); - } - else if (p.is_unilinear() && q.is_val()) { - // a*x + b <=u d - v = p.var(); - a = p.hi().val(); - b = p.lo().val(); - c = rational::zero(); - d = q.val(); - } - else if (p.is_val() && q.is_unilinear()) { - // b <=u c*x + d - v = q.var(); - a = rational::zero(); - b = p.val(); - c = q.hi().val(); - d = q.lo().val(); - } - if (v != null_var) { - s.push_cjust(v, {this, is_positive}); - - s.m_viable.intersect_ule(v, a, b, c, d, is_positive); - - rational val; - if (s.m_viable.find_viable(v, val) == dd::find_t::singleton) - s.propagate(v, val, {this, is_positive}); - - return; - } - - // TODO: other cheap constraints possible? -#endif } bool ule_constraint::is_always_false(bool is_positive, pdd const& lhs, pdd const& rhs) const { diff --git a/src/math/polysat/viable.cpp b/src/math/polysat/viable.cpp index 835c1b7a8..5ccb7dc05 100644 --- a/src/math/polysat/viable.cpp +++ b/src/math/polysat/viable.cpp @@ -10,133 +10,349 @@ Author: Nikolaj Bjorner (nbjorner) 2021-03-19 Jakob Rath 2021-04-6 -Notes: - --*/ +#include "util/debug.h" #include "math/polysat/viable.h" #include "math/polysat/solver.h" + namespace polysat { - viable::viable(solver& s): - s(s), - m_bdd(1000) - {} + viable::viable(solver& s) : s(s) {} viable::~viable() { + for (entry* e : m_alloc) + dealloc(e); } - void viable::push_viable(pvar v) { - s.m_trail.push_back(trail_instr_t::viable_add_i); - m_viable_trail.push_back(std::make_pair(v, m_viable[v])); - + viable::entry* viable::alloc_entry() { + if (m_alloc.empty()) + return alloc(entry); + auto* e = m_alloc.back(); + e->side_cond.reset(); + m_alloc.pop_back(); + return e; } void viable::pop_viable() { - auto const & p = m_viable_trail.back(); - m_viable.set(p.first, p.second); - m_viable_trail.pop_back(); + auto& [v, e] = m_trail.back(); + e->remove_from(m_viable[v], e); + m_alloc.push_back(e); + m_trail.pop_back(); } + void viable::push_viable() { + auto& [v, e] = m_trail.back(); + SASSERT(e->prev() != e || !m_viable[v]); + SASSERT(e->prev() != e || e->next() == e); + if (e->prev() != e) { + e->prev()->insert_after(e); + if (e->interval.lo_val() < e->next()->interval.lo_val()) + m_viable[v] = e; + } + else + m_viable[v] = e; + m_trail.pop_back(); + } - // a*v + b == 0 or a*v + b != 0 - void viable::intersect_eq(rational const& a, pvar v, rational const& b, bool is_positive) { - - bddv const& x = var2bits(v).var(); - if (b == 0 && a.is_odd()) { - // hacky test optimizing special case. - // general case is compute inverse(a)*-b for equality 2^k*a*x + b == 0 - // then constrain x. - // - intersect_viable(v, is_positive ? x.all0() : !x.all0()); - } - else if (a.is_odd()) { - rational a_inv; - VERIFY(a.mult_inverse(x.size(), a_inv)); - bdd eq = x == mod(a_inv * -b, rational::power_of_two(x.size())); - intersect_viable(v, is_positive ? eq : !eq); - } + void viable::intersect(pvar v, signed_constraint const& c) { + auto& fi = s.m_forbidden_intervals; + entry* ne = alloc_entry(); + if (!fi.get_interval(c, v, ne->interval, ne->side_cond) || ne->interval.is_currently_empty()) + m_alloc.push_back(ne); else { - IF_VERBOSE(10, verbose_stream() << a << "*x + " << b << "\n"); - - bddv lhs = a * x + b; - bdd xs = is_positive ? lhs.all0() : !lhs.all0(); - intersect_viable(v, xs); + ne->src = c; + intersect(v, ne); } } - void viable::intersect_ule(pvar v, rational const& a, rational const& b, rational const& c, rational const& d, bool is_positive) { - bddv const& x = var2bits(v).var(); - // hacky special case - if (a == 1 && b == 0 && c == 0 && d == 0) - // x <= 0 - intersect_viable(v, is_positive ? x.all0() : !x.all0()); + void viable::intersect(pvar v, entry* ne) { + entry* e = m_viable[v]; + if (e && e->interval.is_full()) + return; + + if (ne->interval.is_currently_empty()) { + m_alloc.push_back(ne); + return; + } + + auto create_entry = [&]() { + m_trail.push_back({ v, ne }); + s.m_trail.push_back(trail_instr_t::viable_add_i); + ne->init(ne); + return ne; + }; + + auto remove_entry = [&](entry* e) { + m_trail.push_back({ v, e }); + s.m_trail.push_back(trail_instr_t::viable_rem_i); + e->remove_from(m_viable[v], e); + }; + + //LOG("intersect " << ne->interval); + + if (!e) + m_viable[v] = create_entry(); else { - IF_VERBOSE(10, verbose_stream() << a << "*x + " << b << (is_positive ? " <= " : " > ") << c << "*x + " << d << "\n"); - bddv l = a * x + b; - bddv r = c * x + d; - bdd xs = is_positive ? (l <= r) : (l > r); - intersect_viable(v, xs); + entry* first = e; + do { + if (e->interval.contains(ne->interval)) { + m_alloc.push_back(ne); + return; + } + while (ne->interval.contains(e->interval)) { + entry* n = e->next(); + remove_entry(e); + if (!m_viable[v]) { + m_viable[v] = create_entry(); + return; + } + if (e == first) + first = n; + e = n; + } + SASSERT(e->interval.lo_val() != ne->interval.lo_val()); + if (e->interval.lo_val() > ne->interval.lo_val()) { + if (first->prev()->interval.contains(ne->interval)) { + m_alloc.push_back(ne); + return; + } + e->insert_before(create_entry()); + if (e == first) + m_viable[v] = e->prev(); + SASSERT(well_formed(m_viable[v])); + return; + } + e = e->next(); + } + while (e != first); + // otherwise, append to end of list + first->insert_before(create_entry()); } + SASSERT(well_formed(m_viable[v])); } - bool viable::has_viable(pvar v) { - return !m_viable[v].is_false(); + bool viable::has_viable(pvar v) { + auto* e = m_viable[v]; + if (!e) + return true; + entry* first = e; + auto const& max_value = s.var2pdd(v).max_value(); + do { + if (e->interval.is_full()) + return false; + entry* n = e->next(); + if (n == e) + return true; + if (e->interval.hi_val() < n->interval.lo_val()) + return true; + if (n == first) + return e->interval.lo_val() <= e->interval.hi_val(); + e = n; + } + while (e != first); + return false; } - bool viable::is_viable(pvar v, rational const& val) { - return var2bits(v).contains(m_viable[v], val); + bool viable::is_viable(pvar v, rational const& val) { + auto* e = m_viable[v]; + if (!e) + return true; + entry* first = e; + entry* last = first->prev(); + if (last->interval.currently_contains(val)) + return false; + for (; e != last; e = e->next()) { + if (e->interval.currently_contains(val)) + return false; + if (val < e->interval.lo_val()) + return true; + } + return true; } - void viable::intersect_viable(pvar v, bdd vals) { - push_viable(v); - m_viable[v] &= vals; - if (m_viable[v].is_false()) - s.set_conflict(v); + rational viable::min_viable(pvar v) { + rational lo(0); + auto* e = m_viable[v]; + if (!e) + return lo; + entry* first = e; + entry* last = first->prev(); + if (last->interval.currently_contains(lo)) + lo = last->interval.hi_val(); + do { + if (!e->interval.currently_contains(lo)) + break; + lo = e->interval.hi_val(); + e = e->next(); + } + while (e != first); + SASSERT(is_viable(v, lo)); + return lo; } - dd::find_t viable::find_viable(pvar v, rational & val) { - return var2bits(v).find_hint(m_viable[v], s.m_value[v], val); + rational viable::max_viable(pvar v) { + rational hi = s.var2pdd(v).max_value(); + auto* e = m_viable[v]; + if (!e) + return hi; + entry* last = e->prev(); + e = last; + do { + if (!e->interval.currently_contains(hi)) + break; + hi = e->interval.lo_val() - 1; + e = e->prev(); + } + while (e != last); + SASSERT(is_viable(v, hi)); + return hi; } - rational viable::min_viable(pvar v) { - return var2bits(v).min(m_viable[v]); + dd::find_t viable::find_viable(pvar v, rational& lo) { + lo = 0; + auto* e = m_viable[v]; + if (!e) + return dd::find_t::multiple; + if (e->interval.is_full()) + return dd::find_t::empty; + + entry* first = e; + entry* last = first->prev(); + if (last->interval.currently_contains(lo)) + lo = last->interval.hi_val(); + do { + if (!e->interval.currently_contains(lo)) + break; + lo = e->interval.hi_val(); + e = e->next(); + } + while (e != first); + + if (e->interval.currently_contains(lo)) + return dd::find_t::empty; + + rational hi = s.var2pdd(v).max_value(); + e = last; + do { + if (!e->interval.currently_contains(hi)) + break; + hi = e->interval.lo_val() - 1; + e = e->prev(); + } + while (e != last); + if (lo == hi) + return dd::find_t::singleton; + else + return dd::find_t::multiple; } - rational viable::max_viable(pvar v) { - return var2bits(v).max(m_viable[v]); - } + bool viable::resolve(pvar v, conflict& core) { + if (has_viable(v)) + return false; + auto* e = m_viable[v]; + entry* first = e; + SASSERT(e); + core.reset(); + do { + // Build constraint: upper bound of each interval is not contained in the next interval, + // using the equivalence: t \in [l;h[ <=> t-l < h-l + entry* n = e->next(); + if (!e->interval.is_full()) { + auto const& hi = e->interval.hi(); + auto const& next_lo = n->interval.lo(); + auto const& next_hi = n->interval.hi(); + auto lhs = hi - next_lo; + auto rhs = next_hi - next_lo; + signed_constraint c = s.m_constraints.ult(lhs, rhs); + core.insert(c); + } + for (auto sc : e->side_cond) + core.insert(sc); + e->src->set_var_dependent(); // ? + core.insert(e->src); + e = n; + } + while (e != first); - dd::fdd const& viable::sz2bits(unsigned sz) { - m_bits.reserve(sz + 1); - auto* bits = m_bits[sz]; - if (!bits) { - m_bits.set(sz, alloc(dd::fdd, m_bdd, sz)); - bits = m_bits[sz]; + core.set_bailout(); + for (auto c : core) { + if (c.bvalue(s) == l_false) { + core.reset(); + core.set(~c); + break; + } } - return *bits; - } - - void viable::log() { - // only for small problems - for (pvar v = 0; v < std::min(10u, m_viable.size()); ++v) - log(v); + return true; } void viable::log(pvar v) { - if (s.size(v) <= 5) { - vector xs; - for (rational x = rational::zero(); x < rational::power_of_two(s.size(v)); x += 1) - if (is_viable(v, x)) - xs.push_back(x); - - LOG("Viable for v" << v << ": " << xs); - } + if (!well_formed(m_viable[v])) + LOG("v" << v << " not well formed"); + auto* e = m_viable[v]; + if (!e) + return; + entry* first = e; + do { + LOG("v" << v << ": " << e->interval << " " << e->side_cond << " " << e->src); + e = e->next(); + } + while (e != first); } - dd::fdd const& viable::var2bits(pvar v) { return sz2bits(s.size(v)); } + void viable::log() { + for (pvar v = 0; v < std::min(10u, m_viable.size()); ++v) + log(v); + } + + std::ostream& viable::display(std::ostream& out, pvar v) const { + auto* e = m_viable[v]; + if (!e) + return out; + entry* first = e; + do { + out << "v" << v << ": " << e->interval << " " << e->side_cond << " " << e->src << "\n"; + e = e->next(); + } + while (e != first); + return out; + } + + std::ostream& viable::display(std::ostream& out) const { + for (pvar v = 0; v < m_viable.size(); ++v) + display(out, v); + return out; + } + + /* + * Lower bounds are strictly ascending. + * intervals don't contain each-other (since lower bounds are ascending, + * it suffices to check containment in one direction). + */ + bool viable::well_formed(entry* e) { + if (!e) + return true; + entry* first = e; + while (true) { + if (e->interval.is_full()) + return e->next() == e; + if (e->interval.is_currently_empty()) + return false; + + auto* n = e->next(); + if (n != e && e->interval.contains(n->interval)) + return false; + + if (n == first) + break; + if (e->interval.lo_val() >= n->interval.lo_val()) + return false; + e = n; + } + return true; + } } diff --git a/src/math/polysat/viable.h b/src/math/polysat/viable.h index e2792ec65..20023da4f 100644 --- a/src/math/polysat/viable.h +++ b/src/math/polysat/viable.h @@ -4,77 +4,67 @@ Copyright (c) 2021 Microsoft Corporation Module Name: maintain viable domains - + It uses the interval extraction functions from forbidden intervals. + An empty viable set corresponds directly to a conflict that does not rely on + the non-viable variable. Author: Nikolaj Bjorner (nbjorner) 2021-03-19 Jakob Rath 2021-04-6 -Notes: - - NEW_VIABLE uses cheaper book-keeping, but is partial. - - --*/ #pragma once -#include -#include "math/dd/dd_bdd.h" -#include "math/polysat/types.h" +#include +#include "util/dlist.h" +#include "util/small_object_allocator.h" +#include "math/polysat/types.h" +#include "math/polysat/conflict.h" +#include "math/polysat/constraint.h" namespace polysat { class solver; - class viable { - typedef dd::bdd bdd; - typedef dd::fdd fdd; - solver& s; - dd::bdd_manager m_bdd; - scoped_ptr_vector m_bits; - vector m_viable; // set of viable values. - vector> m_viable_trail; + solver& s; + + struct entry : public dll_base, public fi_record { + public: + entry() : fi_record({ eval_interval::full(), {}, {} }) {} + }; + + ptr_vector m_alloc; + ptr_vector m_viable; // set of viable values. + svector> m_trail; // undo stack + bool well_formed(entry* e); - /** - * Register all values that are not contained in vals as non-viable. - */ - void intersect_viable(pvar v, bdd vals); + entry* alloc_entry(); - dd::bdd_manager& get_bdd() { return m_bdd; } - dd::fdd const& sz2bits(unsigned sz); - dd::fdd const& var2bits(pvar v); - - - void push_viable(pvar v); + void intersect(pvar v, entry* e); public: viable(solver& s); ~viable(); - void push(unsigned num_bits) { - m_viable.push_back(m_bdd.mk_true()); - } + // declare and remove var + void push(unsigned) { m_viable.push_back(nullptr); } - void pop() { - m_viable.pop_back(); - } + void pop() { m_viable.pop_back(); } void pop_viable(); - void push_viable() {} + void push_viable(); + /** * update state of viable for pvar v * based on affine constraints */ - - void intersect_eq(rational const& a, pvar v, rational const& b, bool is_positive); - - void intersect_ule(pvar v, rational const& a, rational const& b, rational const& c, rational const& d, bool is_positive); + void intersect(pvar v, signed_constraint const& c); /** * Check whether variable v has any viable values left according to m_viable. @@ -98,6 +88,12 @@ namespace polysat { */ dd::find_t find_viable(pvar v, rational & val); + /** + * Retrieve the unsat core for v. + * \pre there are no viable values for v + */ + bool resolve(pvar v, conflict& core); + /** Log all viable values for the given variable. * (Inefficient, but useful for debugging small instances.) */ @@ -106,7 +102,64 @@ namespace polysat { /** Like log(v) but for all variables */ void log(); + std::ostream& display(std::ostream& out) const; + + class iterator { + entry* curr = nullptr; + bool visited = false; + public: + iterator(entry* curr, bool visited) : + curr(curr), visited(visited || !curr) {} + + iterator& operator++() { + visited = true; + curr = curr->next(); + return *this; + } + + signed_constraint& operator*() { + return curr->src; + } + + bool operator==(iterator const& other) const { + return visited == other.visited && curr == other.curr; + } + + bool operator!=(iterator const& other) const { + return !(*this == other); + } + }; + + class constraints { + viable& v; + pvar var; + public: + constraints(viable& v, pvar var) : v(v), var(var) {} + iterator begin() const { return iterator(v.m_viable[var], false); } + iterator end() const { return iterator(v.m_viable[var], true); } + }; + + constraints get_constraints(pvar v) { + return constraints(*this, v); + } + + std::ostream& display(std::ostream& out, pvar v) const; + + struct var_pp { + viable& v; + pvar var; + var_pp(viable& v, pvar var) : v(v), var(var) {} + }; + }; + + inline std::ostream& operator<<(std::ostream& out, viable const& v) { + return v.display(out); + } + + inline std::ostream& operator<<(std::ostream& out, viable::var_pp const& v) { + return v.v.display(out, v.var); + } } diff --git a/src/math/polysat/viable2.cpp b/src/math/polysat/viable2.cpp deleted file mode 100644 index b969f274b..000000000 --- a/src/math/polysat/viable2.cpp +++ /dev/null @@ -1,359 +0,0 @@ -/*++ -Copyright (c) 2021 Microsoft Corporation - -Module Name: - - maintain viable domains - -Author: - - Nikolaj Bjorner (nbjorner) 2021-03-19 - Jakob Rath 2021-04-6 - ---*/ - - -#include "util/debug.h" -#include "math/polysat/viable2.h" -#include "math/polysat/solver.h" - - -namespace polysat { - - viable2::viable2(solver& s) : s(s) {} - - viable2::~viable2() { - for (entry* e : m_alloc) - dealloc(e); - } - - viable2::entry* viable2::alloc_entry() { - if (m_alloc.empty()) - return alloc(entry); - auto* e = m_alloc.back(); - e->side_cond.reset(); - m_alloc.pop_back(); - return e; - } - - void viable2::pop_viable() { - auto& [v, e] = m_trail.back(); - e->remove_from(m_viable[v], e); - m_alloc.push_back(e); - m_trail.pop_back(); - } - - void viable2::push_viable() { - auto& [v, e] = m_trail.back(); - SASSERT(e->prev() != e || !m_viable[v]); - SASSERT(e->prev() != e || e->next() == e); - if (e->prev() != e) { - e->prev()->insert_after(e); - if (e->interval.lo_val() < e->next()->interval.lo_val()) - m_viable[v] = e; - } - else - m_viable[v] = e; - m_trail.pop_back(); - } - - void viable2::intersect(pvar v, signed_constraint const& c) { - auto& fi = s.m_forbidden_intervals; - entry* ne = alloc_entry(); - if (!fi.get_interval(c, v, ne->interval, ne->side_cond) || ne->interval.is_currently_empty()) - m_alloc.push_back(ne); - else { - ne->src = c; - intersect(v, ne); - } - } - - void viable2::intersect(pvar v, entry* ne) { - entry* e = m_viable[v]; - if (e && e->interval.is_full()) - return; - - if (ne->interval.is_currently_empty()) { - m_alloc.push_back(ne); - return; - } - - auto create_entry = [&]() { - m_trail.push_back({ v, ne }); - s.m_trail.push_back(trail_instr_t::viable_add_i); - ne->init(ne); - return ne; - }; - - auto remove_entry = [&](entry* e) { - m_trail.push_back({ v, e }); - s.m_trail.push_back(trail_instr_t::viable_rem_i); - e->remove_from(m_viable[v], e); - }; - - //LOG("intersect " << ne->interval); - - if (!e) - m_viable[v] = create_entry(); - else { - entry* first = e; - do { - if (e->interval.contains(ne->interval)) { - m_alloc.push_back(ne); - return; - } - while (ne->interval.contains(e->interval)) { - entry* n = e->next(); - remove_entry(e); - if (!m_viable[v]) { - m_viable[v] = create_entry(); - return; - } - if (e == first) - first = n; - e = n; - } - SASSERT(e->interval.lo_val() != ne->interval.lo_val()); - if (e->interval.lo_val() > ne->interval.lo_val()) { - if (first->prev()->interval.contains(ne->interval)) { - m_alloc.push_back(ne); - return; - } - e->insert_before(create_entry()); - if (e == first) - m_viable[v] = e->prev(); - SASSERT(well_formed(m_viable[v])); - return; - } - e = e->next(); - } - while (e != first); - // otherwise, append to end of list - first->insert_before(create_entry()); - } - SASSERT(well_formed(m_viable[v])); - } - - bool viable2::has_viable(pvar v) { - auto* e = m_viable[v]; - if (!e) - return true; - entry* first = e; - auto const& max_value = s.var2pdd(v).max_value(); - do { - if (e->interval.is_full()) - return false; - entry* n = e->next(); - if (n == e) - return true; - if (e->interval.hi_val() < n->interval.lo_val()) - return true; - if (n == first) - return e->interval.lo_val() <= e->interval.hi_val(); - e = n; - } - while (e != first); - return false; - } - - bool viable2::is_viable(pvar v, rational const& val) { - auto* e = m_viable[v]; - if (!e) - return true; - entry* first = e; - entry* last = first->prev(); - if (last->interval.currently_contains(val)) - return false; - for (; e != last; e = e->next()) { - if (e->interval.currently_contains(val)) - return false; - if (val < e->interval.lo_val()) - return true; - } - return true; - } - - rational viable2::min_viable(pvar v) { - rational lo(0); - auto* e = m_viable[v]; - if (!e) - return lo; - entry* first = e; - entry* last = first->prev(); - if (last->interval.currently_contains(lo)) - lo = last->interval.hi_val(); - do { - if (!e->interval.currently_contains(lo)) - break; - lo = e->interval.hi_val(); - e = e->next(); - } - while (e != first); - SASSERT(is_viable(v, lo)); - return lo; - } - - rational viable2::max_viable(pvar v) { - rational hi = s.var2pdd(v).max_value(); - auto* e = m_viable[v]; - if (!e) - return hi; - entry* last = e->prev(); - e = last; - do { - if (!e->interval.currently_contains(hi)) - break; - hi = e->interval.lo_val() - 1; - e = e->prev(); - } - while (e != last); - SASSERT(is_viable(v, hi)); - return hi; - } - - dd::find_t viable2::find_viable(pvar v, rational& lo) { - lo = 0; - auto* e = m_viable[v]; - if (!e) - return dd::find_t::multiple; - if (e->interval.is_full()) - return dd::find_t::empty; - - entry* first = e; - entry* last = first->prev(); - if (last->interval.currently_contains(lo)) - lo = last->interval.hi_val(); - do { - if (!e->interval.currently_contains(lo)) - break; - lo = e->interval.hi_val(); - e = e->next(); - } - while (e != first); - - if (e->interval.currently_contains(lo)) - return dd::find_t::empty; - - rational hi = s.var2pdd(v).max_value(); - e = last; - do { - if (!e->interval.currently_contains(hi)) - break; - hi = e->interval.lo_val() - 1; - e = e->prev(); - } - while (e != last); - if (lo == hi) - return dd::find_t::singleton; - else - return dd::find_t::multiple; - } - - bool viable2::resolve(pvar v, conflict& core) { - if (has_viable(v)) - return false; - auto* e = m_viable[v]; - entry* first = e; - SASSERT(e); - core.reset(); - do { - // Build constraint: upper bound of each interval is not contained in the next interval, - // using the equivalence: t \in [l;h[ <=> t-l < h-l - entry* n = e->next(); - if (!e->interval.is_full()) { - auto const& hi = e->interval.hi(); - auto const& next_lo = n->interval.lo(); - auto const& next_hi = n->interval.hi(); - auto lhs = hi - next_lo; - auto rhs = next_hi - next_lo; - signed_constraint c = s.m_constraints.ult(lhs, rhs); - core.insert(c); - } - for (auto sc : e->side_cond) - core.insert(sc); - e->src->set_var_dependent(); // ? - core.insert(e->src); - e = n; - } - while (e != first); - - core.set_bailout(); - for (auto c : core) { - if (c.bvalue(s) == l_false) { - core.reset(); - core.set(~c); - break; - } - } - return true; - } - - void viable2::log(pvar v) { - if (!well_formed(m_viable[v])) - LOG("v" << v << " not well formed"); - auto* e = m_viable[v]; - if (!e) - return; - entry* first = e; - do { - LOG("v" << v << ": " << e->interval << " " << e->side_cond << " " << e->src); - e = e->next(); - } - while (e != first); - } - - void viable2::log() { - for (pvar v = 0; v < std::min(10u, m_viable.size()); ++v) - log(v); - } - - std::ostream& viable2::display(std::ostream& out, pvar v) const { - auto* e = m_viable[v]; - if (!e) - return out; - entry* first = e; - do { - out << "v" << v << ": " << e->interval << " " << e->side_cond << " " << e->src << "\n"; - e = e->next(); - } - while (e != first); - return out; - } - - std::ostream& viable2::display(std::ostream& out) const { - for (pvar v = 0; v < m_viable.size(); ++v) - display(out, v); - return out; - } - - /* - * Lower bounds are strictly ascending. - * intervals don't contain each-other (since lower bounds are ascending, - * it suffices to check containment in one direction). - */ - bool viable2::well_formed(entry* e) { - if (!e) - return true; - entry* first = e; - while (true) { - if (e->interval.is_full()) - return e->next() == e; - if (e->interval.is_currently_empty()) - return false; - - auto* n = e->next(); - if (n != e && e->interval.contains(n->interval)) - return false; - - if (n == first) - break; - if (e->interval.lo_val() >= n->interval.lo_val()) - return false; - e = n; - } - return true; - } - - -} - diff --git a/src/math/polysat/viable2.h b/src/math/polysat/viable2.h deleted file mode 100644 index 759d00b3a..000000000 --- a/src/math/polysat/viable2.h +++ /dev/null @@ -1,155 +0,0 @@ -/*++ -Copyright (c) 2021 Microsoft Corporation - -Module Name: - - maintain viable domains - It uses the interval extraction functions from forbidden intervals. - An empty viable set corresponds directly to a conflict that does not rely on - the non-viable variable. - -Author: - - Nikolaj Bjorner (nbjorner) 2021-03-19 - Jakob Rath 2021-04-6 - ---*/ -#pragma once - - -#include -#include "util/dlist.h" -#include "util/small_object_allocator.h" -#include "math/polysat/types.h" -#include "math/polysat/conflict.h" -#include "math/polysat/constraint.h" - -namespace polysat { - - class solver; - - class viable2 { - solver& s; - - struct entry : public dll_base, public fi_record { - public: - entry() : fi_record({ eval_interval::full(), {}, {} }) {} - }; - - ptr_vector m_alloc; - ptr_vector m_viable; // set of viable values. - svector> m_trail; // undo stack - - bool well_formed(entry* e); - - entry* alloc_entry(); - - void intersect(pvar v, entry* e); - - std::ostream& display(std::ostream& out, pvar v) const; - - public: - viable2(solver& s); - - ~viable2(); - - // declare and remove var - void push(unsigned) { m_viable.push_back(nullptr); } - - void pop() { m_viable.pop_back(); } - - void pop_viable(); - - void push_viable(); - - /** - * update state of viable for pvar v - * based on affine constraints - */ - void intersect(pvar v, signed_constraint const& c); - - /** - * Check whether variable v has any viable values left according to m_viable. - */ - bool has_viable(pvar v); - - /** - * check if value is viable according to m_viable. - */ - bool is_viable(pvar v, rational const& val); - - /* - * Extract min and max viable values for v - */ - rational min_viable(pvar v); - - rational max_viable(pvar v); - - /** - * Find a next viable value for variable. - */ - dd::find_t find_viable(pvar v, rational & val); - - /** - * Retrieve the unsat core for v. - * \pre there are no viable values for v - */ - bool resolve(pvar v, conflict& core); - - /** Log all viable values for the given variable. - * (Inefficient, but useful for debugging small instances.) - */ - void log(pvar v); - - /** Like log(v) but for all variables */ - void log(); - - std::ostream& display(std::ostream& out) const; - - class iterator { - entry* curr = nullptr; - bool visited = false; - public: - iterator(entry* curr, bool visited) : - curr(curr), visited(visited || !curr) {} - - iterator& operator++() { - visited = true; - curr = curr->next(); - return *this; - } - - signed_constraint& operator*() { - return curr->src; - } - - bool operator==(iterator const& other) const { - return visited == other.visited && curr == other.curr; - } - - bool operator!=(iterator const& other) const { - return !(*this == other); - } - }; - - class constraints { - viable2& viable; - pvar var; - public: - constraints(viable2& viable, pvar var) : viable(viable), var(var) {} - iterator begin() const { return iterator(viable.m_viable[var], false); } - iterator end() const { return iterator(viable.m_viable[var], true); } - }; - - constraints get_constraints(pvar v) { - return constraints(*this, v); - } - - }; - - inline std::ostream& operator<<(std::ostream& out, viable2 const& v) { - return v.display(out); - } -} - - diff --git a/src/math/polysat/viable_set.h b/src/math/polysat/viable_set.h deleted file mode 100644 index 69f513ff3..000000000 --- a/src/math/polysat/viable_set.h +++ /dev/null @@ -1,56 +0,0 @@ -/*++ -Copyright (c) 2021 Microsoft Corporation - -Module Name: - - set of viable values as wrap-around interval - -Author: - - Nikolaj Bjorner (nbjorner) 2021-03-19 - Jakob Rath 2021-04-6 - -Notes: - - - replace BDDs by viable sets that emulate affine relations. - viable_set has an interval of feasible values. - it also can use ternary bit-vectors. - or we could also just use a vector of lbool instead of ternary bit-vectors - updating them at individual positions is relatively cheap instead of copying the - vectors every time a range is narrowed. - - ---*/ -#pragma once - - -#include - -#include "math/dd/dd_bdd.h" -#include "math/polysat/types.h" -#include "math/interval/mod_interval.h" - - -namespace polysat { - - - class viable_set : public mod_interval { - unsigned m_num_bits; - rational p2() const { return rational::power_of_two(m_num_bits); } - bool is_max(rational const& a) const override; - void intersect_eq(rational const& a, bool is_positive); - bool narrow(std::function& eval); - public: - viable_set(unsigned num_bits): m_num_bits(num_bits) {} - ~viable_set() override {} - dd::find_t find_hint(rational const& c, rational& val) const; - bool intersect_eq(rational const& a, rational const& b, bool is_positive); - bool intersect_le(rational const& a, rational const& b, rational const& c, rational const& d, bool is_positive); - rational prev(rational const& p) const; - rational next(rational const& p) const; - }; - -} - - diff --git a/src/math/polysat/viable_set_def.h b/src/math/polysat/viable_set_def.h deleted file mode 100644 index 6695997ab..000000000 --- a/src/math/polysat/viable_set_def.h +++ /dev/null @@ -1,121 +0,0 @@ -/*++ -Copyright (c) 2021 Microsoft Corporation - -Module Name: - - set of viable values as wrap-around interval - -Author: - - Nikolaj Bjorner (nbjorner) 2021-03-19 - Jakob Rath 2021-04-6 - - ---*/ -#pragma once - - -#include "math/polysat/viable_set.h" -#include "math/interval/mod_interval_def.h" - -namespace polysat { - - dd::find_t viable_set::find_hint(rational const& d, rational& val) const { - if (is_empty()) - return dd::find_t::empty; - if (contains(d)) - val = d; - else - val = lo; - if (lo + 1 == hi || hi == 0 && is_max(lo)) - return dd::find_t::singleton; - return dd::find_t::multiple; - } - - bool viable_set::is_max(rational const& a) const { - return a + 1 == rational::power_of_two(m_num_bits); - } - - void viable_set::intersect_eq(rational const& a, bool is_positive) { - if (is_positive) - intersect_fixed(a); - else - intersect_diff(a); - } - - bool viable_set::intersect_eq(rational const& a, rational const& b, bool is_positive) { - if (!a.is_odd()) { - std::function eval = [&](rational const& x) { - return is_positive == (mod(a * x + b, p2()) == 0); - }; - return narrow(eval); - } - if (b == 0) - intersect_eq(b, is_positive); - else { - rational a_inv; - VERIFY(a.mult_inverse(m_num_bits, a_inv)); - intersect_eq(mod(a_inv * -b, p2()), is_positive); - } - return true; - } - - bool viable_set::intersect_le(rational const& a, rational const& b, rational const& c, rational const& d, bool is_positive) { - // x <= 0 - if (a.is_odd() && b == 0 && c == 0 && d == 0) - intersect_eq(b, is_positive); - else if (a == 1 && b == 0 && c == 0) { - // x <= d or x > d - if (is_positive) - intersect_ule(d); - else - intersect_ugt(d); - } - else if (a == 0 && c == 1 && d == 0) { - // x >= b or x < b - if (is_positive) - intersect_uge(b); - else - intersect_ult(b); - } - // TBD: can also handle wrap-around semantics (for signed comparison) - else { - std::function eval = [&](rational const& x) { - return is_positive == mod(a * x + b, p2()) <= mod(c * x + d, p2()); - }; - return narrow(eval); - } - return true; - } - - rational viable_set::prev(rational const& p) const { - if (p > 0) - return p - 1; - else - return rational::power_of_two(m_num_bits) - 1; - } - - rational viable_set::next(rational const& p) const { - if (is_max(p)) - return rational(0); - else - return p + 1; - } - - bool viable_set::narrow(std::function& eval) { - unsigned budget = 10; - while (budget > 0 && !is_empty() && !eval(lo)) { - --budget; - intersect_diff(lo); - } - while (budget > 0 && !is_empty() && !eval(prev(hi))) { - --budget; - intersect_diff(prev(hi)); - } - return 0 < budget; - } - - -} - - diff --git a/src/test/viable.cpp b/src/test/viable.cpp index 352f0c76f..dceda0be6 100644 --- a/src/test/viable.cpp +++ b/src/test/viable.cpp @@ -1,6 +1,6 @@ #include "math/polysat/log.h" #include "math/polysat/solver.h" -#include "math/polysat/viable2.h" +#include "math/polysat/viable.h" namespace polysat { @@ -10,7 +10,7 @@ namespace polysat { class scoped_solverv : public solver_scopev, public solver { public: - viable2 v; + viable v; scoped_solverv(): solver(lim), v(*this) {} ~scoped_solverv() { for (unsigned i = m_trail.size(); i-- > 0;) {