From fc606907428bf7b1b8b335f367d1766b18872d4a Mon Sep 17 00:00:00 2001 From: Nikolaj Bjorner Date: Tue, 20 Apr 2021 12:21:50 -0700 Subject: [PATCH] tidy, initial polysat Signed-off-by: Nikolaj Bjorner --- src/math/dd/dd_bdd.cpp | 48 +++--- src/math/dd/dd_bdd.h | 5 +- src/math/polysat/fixplex.h | 195 ++++++++++++++++++++++++ src/math/polysat/fixplex_def.h | 264 +++++++++++++++++++++++++++++++++ src/math/polysat/solver.cpp | 1 + src/math/polysat/solver.h | 1 + src/test/bdd.cpp | 4 +- src/test/polysat.cpp | 10 +- 8 files changed, 496 insertions(+), 32 deletions(-) create mode 100644 src/math/polysat/fixplex.h create mode 100644 src/math/polysat/fixplex_def.h diff --git a/src/math/dd/dd_bdd.cpp b/src/math/dd/dd_bdd.cpp index 41706266d..d36a82e63 100644 --- a/src/math/dd/dd_bdd.cpp +++ b/src/math/dd/dd_bdd.cpp @@ -976,17 +976,10 @@ namespace dd { bddv a_shifted = a; bddv result = mk_zero(a.size()); for (unsigned i = 0; i < b.size(); ++i) { -#if 1 bddv s = a_shifted; for (unsigned j = i; j < b.size(); ++j) s[j] &= b[i]; result = mk_add(result, s); -#else - // From BuDDy's bvec_mul. It seems to compute more intermediate BDDs than the version above? - bddv added = mk_add(result, a_shifted); - for (unsigned j = 0; j < result.size(); ++j) - result[j] = mk_ite(b[i], added[j], result[j]); -#endif a_shifted.shl(); } return result; @@ -1014,31 +1007,40 @@ namespace dd { return mk_mul(a, [b](unsigned i) { return b[i]; }); } + bddv bdd_manager::mk_concat(bddv const& a, bddv const& b) { + bddv result = a; + result.m_bits.append(b.m_bits); + return result; + } + + + /** + * Quotient remainder + * + * rem, div have size 2*|a| = worksize. + * Initialization: + * rem := a ++ false + * div := false ++ b + */ void bdd_manager::mk_quot_rem(bddv const& a, bddv const& b, bddv& quot, bddv& rem) { SASSERT(a.size() == b.size()); quot = mk_zero(a.size()); - // We work with double-size vectors unsigned worksize = a.size() + b.size(); - // Extend dividend to worksize - rem = a; - for (unsigned i = b.size(); i-- > 0; ) - rem.push_back(mk_false()); - // Shift divisor to the left - bddv div(this); - for (unsigned i = a.size(); i-- > 0; ) - div.push_back(mk_false()); - div.m_bits.append(b.m_bits); + rem = a.append(mk_zero(b.size())); + bddv div = mk_zero(a.size()).append(b); + // // Keep shifting divisor to the right and subtract whenever it is // smaller than the remaining value - for (int i = 0; i <= b.size(); ++i) { + // + for (unsigned i = 0; i <= b.size(); ++i) { bdd divLteRem = div <= rem; bddv remSubDiv = rem - div; - for (int j = 0; j < worksize; ++j) + for (unsigned j = 0; j < worksize; ++j) rem[j] = mk_ite(divLteRem, remSubDiv[j], rem[j]); if (i > 0) - quot[b.size()-i] = divLteRem; + quot[b.size() - i] = divLteRem; div.shr(); } @@ -1105,14 +1107,14 @@ namespace dd { void bddv::shl() { for (unsigned j = size(); j-- > 1;) - m_bits[j] = m_bits[j-1]; + m_bits[j] = m_bits[j - 1]; m_bits[0] = m->mk_false(); } void bddv::shr() { for (unsigned j = 1; j < size(); ++j) - m_bits[j-1] = m_bits[j]; - m_bits[size()-1] = m->mk_false(); + m_bits[j - 1] = m_bits[j]; + m_bits[size() - 1] = m->mk_false(); } } diff --git a/src/math/dd/dd_bdd.h b/src/math/dd/dd_bdd.h index b86b2ebab..1d618c3e3 100644 --- a/src/math/dd/dd_bdd.h +++ b/src/math/dd/dd_bdd.h @@ -251,6 +251,7 @@ namespace dd { bddv mk_mul(bddv const& a, bddv const& b); bddv mk_mul(bddv const& a, bool_vector const& b); bddv mk_mul(bddv const& a, rational const& val); + bddv mk_concat(bddv const& a, bddv const& b); void mk_quot_rem(bddv const& a, bddv const& b, bddv& quot, bddv& rem); bool is_constv(bddv const& a); rational to_val(bddv const& a); @@ -303,7 +304,7 @@ namespace dd { vector m_bits; bdd_manager* m; - bddv(vector bits, bdd_manager* m): m_bits(bits), m(m) { SASSERT(m); } + bddv(vector const& bits, bdd_manager* m): m_bits(bits), m(m) { SASSERT(m); } bddv(vector&& bits, bdd_manager* m): m_bits(std::move(bits)), m(m) { SASSERT(m); } bddv(bdd_manager* m): m_bits(), m(m) { SASSERT(m); } @@ -347,7 +348,7 @@ namespace dd { bddv operator*(bddv const& other) const { return m->mk_mul(*this, other); } bddv operator*(rational const& other) const { return m->mk_mul(*this, other); } bddv operator*(bool_vector const& other) const { return m->mk_mul(*this, other); } - + bddv append(bddv const& other) const { return m->mk_concat(*this, other); } void quot_rem(bddv const& divisor, bddv& quot, bddv& rem) const { m->mk_quot_rem(*this, divisor, quot, rem); } bool is_const() const { return m->is_constv(*this); } diff --git a/src/math/polysat/fixplex.h b/src/math/polysat/fixplex.h new file mode 100644 index 000000000..82691bc75 --- /dev/null +++ b/src/math/polysat/fixplex.h @@ -0,0 +1,195 @@ +/*++ +Copyright (c) 2014 Microsoft Corporation + +Module Name: + + fixplex.h + +Abstract: + + Fixed-precision unsigned integer simplex tableau. + +Author: + + Nikolaj Bjorner (nbjorner) 2021-03-19 + Jakob Rath 2021-04-6 + +--*/ + +#pragma once + +#include "math/simplex/sparse_matrix.h" +#include "util/heap.h" +#include "util/lbool.h" +#include "util/uint_set.h" + +namespace polysat { + + template + class fixplex { + + typedef unsigned var_t; + typedef typename Ext::numeral numeral; + typedef typename Ext::scoped_numeral scoped_numeral; + typedef typename Ext::manager manager; + typedef simplex::sparse_matrix matrix; + struct var_lt { + bool operator()(var_t v1, var_t v2) const { return v1 < v2; } + }; + typedef heap var_heap; + + struct stats { + unsigned m_num_pivots; + unsigned m_num_infeasible; + unsigned m_num_checks; + stats() { reset(); } + void reset() { + memset(this, 0, sizeof(*this)); + } + }; + + enum pivot_strategy_t { + S_BLAND, + S_GREATEST_ERROR, + S_LEAST_ERROR, + S_DEFAULT + }; + + struct var_info { + unsigned m_base2row:29; + unsigned m_is_base:1; + numeral m_lo; + numeral m_hi; + numeral m_value; + var_info(): + m_base2row(0), + m_is_base(false) + {} + }; + + static const var_t null_var; + reslimit& m_limit; + mutable manager m; + mutable matrix M; + unsigned m_max_iterations; + var_heap m_to_patch; + vector m_vars; + svector m_row2base; + bool m_bland; + unsigned m_blands_rule_threshold; + random_gen m_random; + uint_set m_left_basis; + unsigned m_infeasible_var; + unsigned_vector m_base_vars; + unsigned_vector m_to_fix_base; + stats m_stats; + + public: + fixplex(reslimit& lim): + m_limit(lim), + M(m), + m_max_iterations(UINT_MAX), + m_to_patch(1024), + m_bland(false), + m_blands_rule_threshold(1000) {} + + ~fixplex(); + + typedef typename matrix::row row; + typedef typename matrix::row_iterator row_iterator; + typedef typename matrix::col_iterator col_iterator; + + var_t get_base_var(row const& r) const { return m_row2base[r.id()]; } + numeral const& get_lo(var_t var) const { return m_vars[var].m_lo; } + numeral const& get_hi(var_t var) const { return m_vars[var].m_hi; } + void set_max_iterations(unsigned n) { m_max_iterations = n; } + row_iterator row_begin(row const& r) { return M.row_begin(r); } + row_iterator row_end(row const& r) { return M.row_end(r); } + unsigned get_num_vars() const { return m_vars.size(); } + void ensure_var(var_t v); + void reset(); + lbool make_feasible(); + row add_row(var_t base, unsigned num_vars, var_t const* vars, numeral const* coeffs); + +#if 0 + row get_infeasible_row(); + void del_row(var_t base_var); + void set_lo(var_t var, numeral const& b); + void set_hi(var_t var, numeral const& b); + bool in_range(var_t var, numeral const& b) const; + void unset_lo(var_t var); + void unset_hi(var_t var); + void set_value(var_t var, numeral const& b); + numeral const& get_value(var_t v); + void display(std::ostream& out) const; + void display_row(std::ostream& out, row const& r, bool values = true); + + + + void collect_statistics(::statistics & st) const; + +#endif + + private: + + void gauss_jordan(); + bool gauss_jordan(row const& r); + +#if 0 + void del_row(row const& r); + var_t select_var_to_fix(); + pivot_strategy_t pivot_strategy(); + var_t select_smallest_var() { return m_to_patch.empty()?null_var:m_to_patch.erase_min(); } + var_t select_error_var(bool least); + void check_blands_rule(var_t v, unsigned& num_repeated); + bool make_var_feasible(var_t x_i); + void update_and_pivot(var_t x_i, var_t x_j, numeral const& a_ij, numeral const& new_value); + void update_value(var_t v, numeral const& delta); + void update_value_core(var_t v, numeral const& delta); + void pivot(var_t x_i, var_t x_j, numeral const& a_ij); + void move_to_bound(var_t x, bool to_lower); + var_t select_pivot(var_t x_i, bool is_below, scoped_numeral& out_a_ij); + var_t select_pivot_blands(var_t x_i, bool is_below, scoped_numeral& out_a_ij); + var_t select_pivot_core(var_t x_i, bool is_below, scoped_numeral& out_a_ij); + int get_num_non_free_dep_vars(var_t x_j, int best_so_far); + + var_t pick_var_to_leave(var_t x_j, bool is_pos, + scoped_numeral& gain, scoped_numeral& new_a_ij, bool& inc); + + + void select_pivot_primal(var_t v, var_t& x_i, var_t& x_j, scoped_numeral& a_ij, bool& inc_x_i, bool& inc_x_j); + + + bool at_lower(var_t v) const; + bool at_upper(var_t v) const; + bool above_lower(var_t v) const; + bool below_upper(var_t v) const; + bool outside_bounds(var_t v) const { return below_lower(v) || above_upper(v); } + bool is_free(var_t v) const { return m_vars[v].m_lo == m_vars[v].m_hi; } + bool is_non_free(var_t v) const { return !is_free(v); } + bool is_base(var_t x) const { return m_vars[x].m_is_base; } + void add_patch(var_t v); + + bool well_formed() const; + bool well_formed_row(row const& r) const; + bool is_feasible() const; + +#endif + }; + + struct uint64_ext { + typedef uint64_t numeral; + struct manager { + void reset() {} + void reset(numeral& n) {} + bool is_zero(numeral const& n) const { return n == 0; } + }; + struct scoped_numeral { + scoped_numeral(manager& m) { n = 0; } + numeral n; + numeral& operator()() { return n; } + }; + }; + +}; + diff --git a/src/math/polysat/fixplex_def.h b/src/math/polysat/fixplex_def.h new file mode 100644 index 000000000..5593bf26b --- /dev/null +++ b/src/math/polysat/fixplex_def.h @@ -0,0 +1,264 @@ +/*++ +Copyright (c) 2021 Microsoft Corporation + +Module Name: + + fixplex_def.h + +Abstract: + + Fixed-precision unsigned integer simplex tableau. + +Author: + + Nikolaj Bjorner (nbjorner) 2021-03-19 + Jakob Rath 2021-04-6 + +--*/ + +#pragma once + +#include "math/polysat/fixplex.h" +#include "math/simplex/sparse_matrix_def.h" + +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()); + } + if (m_to_patch.get_bounds() <= v) + m_to_patch.set_bounds(2*v+1); + } + + + template + void fixplex::reset() { + M.reset(); + m_to_patch.reset(); + m_vars.reset(); + m_row2base.reset(); + m_left_basis.reset(); + m_base_vars.reset(); + } + + template + lbool fixplex::make_feasible() { + ++m_stats.m_num_checks; + m_left_basis.reset(); + m_infeasible_var = null_var; + unsigned num_iterations = 0; + unsigned num_repeated = 0; + var_t v = null_var; + m_bland = false; + SASSERT(well_formed()); + while ((v = select_var_to_fix()) != null_var) { + TRACE("simplex", display(tout << "v" << v << "\n");); + if (!m_limit.inc() || num_iterations > m_max_iterations) { + return l_undef; + } + check_blands_rule(v, num_repeated); + if (!make_var_feasible(v)) { + m_to_patch.insert(v); + m_infeasible_var = v; + ++m_stats.m_num_infeasible; + return l_false; + } + ++num_iterations; + } + SASSERT(well_formed()); + return l_true; + } + + template + typename fixplex::row + fixplex::add_row(var_t base_var, unsigned num_vars, var_t const* vars, numeral const* coeffs) { + row r = M.mk_row(); + 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; + bool has_base = false; + for (auto const& e : M.row_entries(r)) { + var_t v = e.m_var; + if (v == base_var) + base_coeff = e.m_coeff; + else { + has_base |= is_base(v); + value += e.m_coeff * m_vars[v].m_value; + } + } + SASSERT(base_coeff != 0); + SASSERT(!is_base(base_var)); + while (m_row2base.size() <= r.id()) + m_row2base.push_back(null_var); + m_row2base[r.id()] = base_var; + m_vars[base_var].m_base2row = r.id(); + m_vars[base_var].m_is_base = true; + m_vars[base_var].m_base_coeff = base_coeff; + m_vars[base_vars].m_value = value / base_coeff; + // TBD: record when base_coeff does not divide value + add_patch(base_var); + if (has_base) { + m_to_fix_base.push_back(r.id()); + gauss_jordan(); + } + SASSERT(well_formed_row(r)); + SASSERT(well_formed()); + return r; + } + + template + void fixplex::gauss_jordan() { + while (!m_to_fix_base.empty()) { + auto rid = m_to_fix_base.back(); + if (gauss_jordan(m_rows[rid])) { + m_to_fix_base.pop_back(); + } + } + } + + template + bool fixplex::gauss_jordan(row const& r) { + auto base_var = m_row2base[r.id()]; + unsigned other_base = null_var; + numeral c1; + for (auto const& e : M.row_entries(r)) { + var_t v = e.m_var; + if (is_base(v) && v != base_var) { + other_base = v; + c1 = e.m_coeff; + break; + } + } + if (null_var == other_base) + return true; + + auto c2 = m_vars[other_base].m_base_coeff; + auto r2 = m_vars[other_base].m_base2row; + unsigned exp1 = trailing_zeros(c1); // exponent of two for v in r + unsigned exp2 = trailing_zeros(c2); // exponent of two for v in r2 + if (exp1 >= exp2) { + // eliminate v from r + } + else { + // eliminate v from r2, + // make v the new base for r + // perform gauss-jordan for both r and r2 (add to queue) + } + + NOT_IMPLEMENTED_YET(); + return false; + } + +#if 0 + /** + \brief Make x_j the new base variable for row of x_i. + x_i is assumed base variable of row r_i. + x_j is assumed to have coefficient a_ij in r_i. + + a_ii*x_i + a_ij*x_j + r_i = 0 + + current value of x_i is v_i + new value of x_i is new_value + a_ii*(x_i + new_value - x_i) + a_ij*((x_i - new_value)*a_ii/a_ij + x_j) + r_i = 0 + + Let r_k be a row where x_j has coefficient x_kj != 0. + r_k <- r_k * a_ij - r_i * a_kj + */ + template + void fixplex::update_and_pivot(var_t x_i, var_t x_j, numeral const& a_ij, numeral const& new_value) { + SASSERT(is_base(x_i)); + SASSERT(!is_base(x_j)); + var_info& x_iI = m_vars[x_i]; + numeral const& a_ii = x_iI.m_base_coeff; + auto theta = x_iI.m_value - new_value; + theta *= a_ii; + theta /= a_ij; + update_value(x_j, theta); + SASSERT(new_value == x_iI.m_value); + pivot(x_i, x_j, a_ij); + } + + template + void fixplex::pivot(var_t x_i, var_t x_j, numeral const& a_ij) { + ++m_stats.m_num_pivots; + var_info& x_iI = m_vars[x_i]; + var_info& x_jI = m_vars[x_j]; + unsigned r_i = x_iI.m_base2row; + m_row2base[r_i] = x_j; + x_jI.m_base2row = r_i; + x_jI.m_base_coeff = a_ij; + x_jI.m_is_base = true; + x_iI.m_is_base = false; + add_patch(x_j); + SASSERT(well_formed_row(row(r_i))); + + col_iterator it = M.col_begin(x_j), end = M.col_end(x_j); + scoped_numeral a_kj(m), g(m); + for (; it != end; ++it) { + row r_k = it.get_row(); + if (r_k.id() != r_i) { + a_kj = it.get_row_entry().m_coeff; + a_kj.neg(); + M.mul(r_k, a_ij); + M.add(r_k, a_kj, row(r_i)); + var_t s = m_row2base[r_k.id()]; + numeral& coeff = m_vars[s].m_base_coeff; + m.mul(coeff, a_ij, coeff); + M.gcd_normalize(r_k, g); + if (!m.is_one(g)) { + m.div(coeff, g, coeff); + } + SASSERT(well_formed_row(row(r_k))); + } + } + SASSERT(well_formed()); + } + + template + void fixplex::update_value(var_t v, eps_numeral const& delta) { + if (em.is_zero(delta)) { + return; + } + update_value_core(v, delta); + col_iterator it = M.col_begin(v), end = M.col_end(v); + + // v <- v + delta + // s*s_coeff + v*v_coeff + R = 0 + // -> + // (v + delta)*v_coeff + (s - delta*v_coeff/s_coeff)*v + R = 0 + for (; it != end; ++it) { + row r = it.get_row(); + var_t s = m_row2base[r.id()]; + var_info& si = m_vars[s]; + scoped_eps_numeral delta2(em); + numeral const& coeff = it.get_row_entry().m_coeff; + em.mul(delta, coeff, delta2); + em.div(delta2, si.m_base_coeff, delta2); + delta2.neg(); + update_value_core(s, delta2); + } + } + + template + void fixplex::update_value_core(var_t v, eps_numeral const& delta) { + eps_numeral& val = m_vars[v].m_value; + em.add(val, delta, val); + if (is_base(v)) { + add_patch(v); + } + } +#endif + +} + diff --git a/src/math/polysat/solver.cpp b/src/math/polysat/solver.cpp index 3a2ffa8de..65bde89c4 100644 --- a/src/math/polysat/solver.cpp +++ b/src/math/polysat/solver.cpp @@ -66,6 +66,7 @@ namespace polysat { solver::solver(reslimit& lim): m_lim(lim), m_bdd(1000), + m_fixplex(m_lim), m_dm(m_value_manager, m_alloc), m_free_vars(m_activity) { diff --git a/src/math/polysat/solver.h b/src/math/polysat/solver.h index 7136ee390..83b87d94f 100644 --- a/src/math/polysat/solver.h +++ b/src/math/polysat/solver.h @@ -48,6 +48,7 @@ namespace polysat { scoped_ptr_vector m_pdd; scoped_ptr_vector m_bits; dd::bdd_manager m_bdd; + fixplex m_fixplex; dep_value_manager m_value_manager; small_object_allocator m_alloc; poly_dep_manager m_dm; diff --git a/src/test/bdd.cpp b/src/test/bdd.cpp index 5fcacb0e8..17f92d887 100644 --- a/src/test/bdd.cpp +++ b/src/test/bdd.cpp @@ -339,7 +339,7 @@ public: bddv const& x = x_dom.var(); vector num; - for (unsigned n = 0; n < (1 << x_dom.num_bits()); ++n) { + for (unsigned n = 0; n < (1ul << x_dom.num_bits()); ++n) { num.push_back(x == rational(n)); SASSERT(x_dom.contains(num[n], rational(n))); rational r; @@ -359,7 +359,7 @@ public: SASSERT(old_levels != m.m_var2level); // ensure that reorder actually did something. // Should still give the correct answer after reordering - for (unsigned n = 0; n < (1 << x_dom.num_bits()); ++n) { + for (unsigned n = 0; n < (1ul << x_dom.num_bits()); ++n) { SASSERT(x_dom.contains(num[n], rational(n))); rational r; SASSERT_EQ(x_dom.find(num[n], r), find_t::singleton); diff --git a/src/test/polysat.cpp b/src/test/polysat.cpp index 0cf3d57a4..63fc99e88 100644 --- a/src/test/polysat.cpp +++ b/src/test/polysat.cpp @@ -148,7 +148,7 @@ namespace polysat { */ static void test_monot1() { scoped_solver s; - auto bw = 5; + auto bw = 5; auto tb1 = s.var(s.add_var(bw)); auto tb2 = s.var(s.add_var(bw)); @@ -179,13 +179,13 @@ namespace polysat { s.add_eq(elastic1 + a - elastic2); - // tb2 = ((v * base2) / elastic2); + // tb2 = ((v * base2) / elastic2); s.add_eq((tb2 * elastic2) + rem3 - (v * base2)); - // quot4 = v / (elastic1 + a); - s.add_eq((quot4 * (elastic1 + a)) + rem4 - v); + // quot4 = v / (elastic1 + a); + s.add_eq((quot4 * (elastic1 + a)) + rem4 - v); - s.add_eq(quot4 + 1 - err); + s.add_eq(quot4 + 1 - err); s.add_ult(tb1, tb2); s.check();