diff --git a/.gitignore b/.gitignore index 3fe3a3110..ffc50c1ba 100644 --- a/.gitignore +++ b/.gitignore @@ -91,4 +91,5 @@ examples/**/obj CMakeSettings.json # Editor temp files *.swp -.DS_Store \ No newline at end of file +.DS_Store +dbg/** diff --git a/src/api/api_ast.cpp b/src/api/api_ast.cpp index 647cc8631..df3059d4b 100644 --- a/src/api/api_ast.cpp +++ b/src/api/api_ast.cpp @@ -1178,7 +1178,7 @@ extern "C" { case OP_BSMOD: return Z3_OP_BSMOD; case OP_BSDIV0: return Z3_OP_BSDIV0; case OP_BUDIV0: return Z3_OP_BUDIV0; - case OP_BSREM0: return Z3_OP_BUREM0; + case OP_BSREM0: return Z3_OP_BSREM0; case OP_BUREM0: return Z3_OP_BUREM0; case OP_BSMOD0: return Z3_OP_BSMOD0; case OP_ULEQ: return Z3_OP_ULEQ; diff --git a/src/ast/array_decl_plugin.cpp b/src/ast/array_decl_plugin.cpp index b336dd2e5..975eaa07a 100644 --- a/src/ast/array_decl_plugin.cpp +++ b/src/ast/array_decl_plugin.cpp @@ -326,7 +326,9 @@ func_decl * array_decl_plugin::mk_array_ext(unsigned arity, sort * const * domai } sort * r = to_sort(s->get_parameter(i).get_ast()); parameter param(i); - return m_manager->mk_func_decl(m_array_ext_sym, arity, domain, r, func_decl_info(m_family_id, OP_ARRAY_EXT, 1, ¶m)); + func_decl_info info(func_decl_info(m_family_id, OP_ARRAY_EXT, 1, ¶m)); + info.set_commutative(true); + return m_manager->mk_func_decl(m_array_ext_sym, arity, domain, r, info); } diff --git a/src/ast/decl_collector.cpp b/src/ast/decl_collector.cpp index ecad96521..5b634abbd 100644 --- a/src/ast/decl_collector.cpp +++ b/src/ast/decl_collector.cpp @@ -48,6 +48,7 @@ bool decl_collector::is_bool(sort * s) { } void decl_collector::visit_func(func_decl * n) { + func_decl* g; if (!m_visited.is_marked(n)) { family_id fid = n->get_family_id(); if (fid == null_family_id) @@ -57,6 +58,8 @@ void decl_collector::visit_func(func_decl * n) { recfun::util u(m()); m_todo.push_back(u.get_def(n).get_rhs()); } + else if (m_ar_util.is_as_array(n, g)) + m_todo.push_back(g); m_visited.mark(n, true); m_trail.push_back(n); } @@ -65,7 +68,8 @@ void decl_collector::visit_func(func_decl * n) { decl_collector::decl_collector(ast_manager & m): m_manager(m), m_trail(m), - m_dt_util(m) { + m_dt_util(m), + m_ar_util(m) { m_basic_fid = m_manager.get_basic_family_id(); m_dt_fid = m_dt_util.get_family_id(); recfun::util rec_util(m); @@ -156,8 +160,15 @@ void decl_collector::collect_deps(sort* s, sort_set& set) { for (unsigned i = 0; i < num_cnstr; i++) { func_decl * cnstr = m_dt_util.get_datatype_constructors(s)->get(i); set.insert(cnstr->get_range()); - for (unsigned j = 0; j < cnstr->get_arity(); ++j) - set.insert(cnstr->get_domain(j)); + for (unsigned j = 0; j < cnstr->get_arity(); ++j) { + sort* n = cnstr->get_domain(j); + set.insert(n); + for (unsigned i = n->get_num_parameters(); i-- > 0; ) { + parameter const& p = n->get_parameter(i); + if (p.is_ast() && is_sort(p.get_ast())) + set.insert(to_sort(p.get_ast())); + } + } } } diff --git a/src/ast/decl_collector.h b/src/ast/decl_collector.h index 876b97188..31cabe796 100644 --- a/src/ast/decl_collector.h +++ b/src/ast/decl_collector.h @@ -23,6 +23,7 @@ Revision History: #include "util/lim_vector.h" #include "ast/ast.h" #include "ast/datatype_decl_plugin.h" +#include "ast/array_decl_plugin.h" class decl_collector { ast_manager & m_manager; @@ -35,6 +36,7 @@ class decl_collector { family_id m_basic_fid; family_id m_dt_fid; datatype_util m_dt_util; + array_util m_ar_util; family_id m_rec_fid; ptr_vector m_todo; diff --git a/src/math/simplex/model_based_opt.cpp b/src/math/simplex/model_based_opt.cpp index 0223e8bff..55f2616b8 100644 --- a/src/math/simplex/model_based_opt.cpp +++ b/src/math/simplex/model_based_opt.cpp @@ -1,1300 +1,1651 @@ -/*++ -Copyright (c) 2016 Microsoft Corporation - -Module Name: - - model_based_opt.cpp - -Abstract: - - Model-based optimization and projection for linear real, integer arithmetic. - -Author: - - Nikolaj Bjorner (nbjorner) 2016-27-4 - -Revision History: - - ---*/ - -#include "math/simplex/model_based_opt.h" -#include "util/uint_set.h" -#include "util/z3_exception.h" - -std::ostream& operator<<(std::ostream& out, opt::ineq_type ie) { - switch (ie) { - case opt::t_eq: return out << " = "; - case opt::t_lt: return out << " < "; - case opt::t_le: return out << " <= "; - case opt::t_mod: return out << " mod "; - } - return out; -} - - -namespace opt { - - /** - * Convert a row ax + coeffs + coeff = value into a definition for x - * x = (value - coeffs - coeff)/a - * as backdrop we have existing assignments to x and other variables that - * satisfy the equality with value, and such that value satisfies - * the row constraint ( = , <= , < , mod) - */ - model_based_opt::def::def(row const& r, unsigned x) { - for (var const & v : r.m_vars) { - if (v.m_id != x) { - m_vars.push_back(v); - } - else { - m_div = -v.m_coeff; - } - } - m_coeff = r.m_coeff; - switch (r.m_type) { - case opt::t_lt: - m_coeff += m_div; - break; - case opt::t_le: - // for: ax >= t, then x := (t + a - 1) div a - if (m_div.is_pos()) { - m_coeff += m_div; - m_coeff -= rational::one(); - } - break; - default: - break; - } - normalize(); - SASSERT(m_div.is_pos()); - } - - model_based_opt::def model_based_opt::def::operator+(def const& other) const { - def result; - vector const& vs1 = m_vars; - vector const& vs2 = other.m_vars; - vector & vs = result.m_vars; - rational c1(1), c2(1); - if (m_div != other.m_div) { - c1 = other.m_div; - c2 = m_div; - } - unsigned i = 0, j = 0; - while (i < vs1.size() || j < vs2.size()) { - unsigned v1 = UINT_MAX, v2 = UINT_MAX; - if (i < vs1.size()) v1 = vs1[i].m_id; - if (j < vs2.size()) v2 = vs2[j].m_id; - if (v1 == v2) { - vs.push_back(vs1[i]); - vs.back().m_coeff *= c1; - vs.back().m_coeff += c2 * vs2[j].m_coeff; - ++i; ++j; - if (vs.back().m_coeff.is_zero()) { - vs.pop_back(); - } - } - else if (v1 < v2) { - vs.push_back(vs1[i]); - vs.back().m_coeff *= c1; - } - else { - vs.push_back(vs2[j]); - vs.back().m_coeff *= c2; - } - } - result.m_div = c1*m_div; - result.m_coeff = (m_coeff*c1) + (other.m_coeff*c2); - result.normalize(); - return result; - } - - model_based_opt::def model_based_opt::def::operator/(rational const& r) const { - def result(*this); - result.m_div *= r; - result.normalize(); - return result; - } - - model_based_opt::def model_based_opt::def::operator*(rational const& n) const { - def result(*this); - for (var& v : result.m_vars) { - v.m_coeff *= n; - } - result.m_coeff *= n; - result.normalize(); - return result; - } - - model_based_opt::def model_based_opt::def::operator+(rational const& n) const { - def result(*this); - result.m_coeff += n * result.m_div; - result.normalize(); - return result; - } - - void model_based_opt::def::normalize() { - if (!m_div.is_int()) { - rational den = denominator(m_div); - SASSERT(den > 1); - for (var& v : m_vars) - v.m_coeff *= den; - m_coeff *= den; - m_div *= den; - - } - if (m_div.is_neg()) { - for (var& v : m_vars) - v.m_coeff.neg(); - m_coeff.neg(); - m_div.neg(); - } - if (m_div.is_one()) - return; - rational g(m_div); - if (!m_coeff.is_int()) - return; - g = gcd(g, m_coeff); - for (var const& v : m_vars) { - if (!v.m_coeff.is_int()) - return; - g = gcd(g, abs(v.m_coeff)); - if (g.is_one()) - break; - } - if (!g.is_one()) { - for (var& v : m_vars) - v.m_coeff /= g; - m_coeff /= g; - m_div /= g; - } - } - - model_based_opt::model_based_opt() { - m_rows.push_back(row()); - } - - bool model_based_opt::invariant() { - for (unsigned i = 0; i < m_rows.size(); ++i) { - if (!invariant(i, m_rows[i])) { - return false; - } - } - return true; - } - -#define PASSERT(_e_) { CTRACE("qe", !(_e_), display(tout, r); display(tout);); SASSERT(_e_); } - - bool model_based_opt::invariant(unsigned index, row const& r) { - vector const& vars = r.m_vars; - for (unsigned i = 0; i < vars.size(); ++i) { - // variables in each row are sorted and have non-zero coefficients - PASSERT(i + 1 == vars.size() || vars[i].m_id < vars[i+1].m_id); - PASSERT(!vars[i].m_coeff.is_zero()); - PASSERT(index == 0 || m_var2row_ids[vars[i].m_id].contains(index)); - } - - PASSERT(r.m_value == eval(r)); - PASSERT(r.m_type != t_eq || r.m_value.is_zero()); - // values satisfy constraints - PASSERT(index == 0 || r.m_type != t_lt || r.m_value.is_neg()); - PASSERT(index == 0 || r.m_type != t_le || !r.m_value.is_pos()); - PASSERT(index == 0 || r.m_type != t_mod || (mod(r.m_value, r.m_mod).is_zero())); - return true; - } - - // a1*x + obj - // a2*x + t2 <= 0 - // a3*x + t3 <= 0 - // a4*x + t4 <= 0 - // a1 > 0, a2 > 0, a3 > 0, a4 < 0 - // x <= -t2/a2 - // x <= -t2/a3 - // determine lub among these. - // then resolve lub with others - // e.g., -t2/a2 <= -t3/a3, then - // replace inequality a3*x + t3 <= 0 by -t2/a2 + t3/a3 <= 0 - // mark a4 as invalid. - // - - // a1 < 0, a2 < 0, a3 < 0, a4 > 0 - // x >= t2/a2 - // x >= t3/a3 - // determine glb among these - // the resolve glb with others. - // e.g. t2/a2 >= t3/a3 - // then replace a3*x + t3 by t3/a3 - t2/a2 <= 0 - // - inf_eps model_based_opt::maximize() { - SASSERT(invariant()); - unsigned_vector bound_trail, bound_vars; - TRACE("opt", display(tout << "tableau\n");); - while (!objective().m_vars.empty()) { - var v = objective().m_vars.back(); - unsigned x = v.m_id; - rational const& coeff = v.m_coeff; - unsigned bound_row_index; - rational bound_coeff; - if (find_bound(x, bound_row_index, bound_coeff, coeff.is_pos())) { - SASSERT(!bound_coeff.is_zero()); - TRACE("opt", display(tout << "update: " << v << " ", objective()); - for (unsigned above : m_above) { - display(tout << "resolve: ", m_rows[above]); - }); - for (unsigned above : m_above) { - resolve(bound_row_index, bound_coeff, above, x); - } - for (unsigned below : m_below) { - resolve(bound_row_index, bound_coeff, below, x); - } - // coeff*x + objective <= ub - // a2*x + t2 <= 0 - // => coeff*x <= -t2*coeff/a2 - // objective + t2*coeff/a2 <= ub - - mul_add(false, m_objective_id, - coeff/bound_coeff, bound_row_index); - retire_row(bound_row_index); - bound_trail.push_back(bound_row_index); - bound_vars.push_back(x); - } - else { - TRACE("opt", display(tout << "unbound: " << v << " ", objective());); - update_values(bound_vars, bound_trail); - return inf_eps::infinity(); - } - } - - // - // update the evaluation of variables to satisfy the bound. - // - - update_values(bound_vars, bound_trail); - - rational value = objective().m_value; - if (objective().m_type == t_lt) { - return inf_eps(inf_rational(value, rational(-1))); - } - else { - return inf_eps(inf_rational(value)); - } - } - - - void model_based_opt::update_value(unsigned x, rational const& val) { - rational old_val = m_var2value[x]; - m_var2value[x] = val; - SASSERT(val.is_int() || !is_int(x)); - unsigned_vector const& row_ids = m_var2row_ids[x]; - for (unsigned row_id : row_ids) { - rational coeff = get_coefficient(row_id, x); - if (coeff.is_zero()) { - continue; - } - row & r = m_rows[row_id]; - rational delta = coeff * (val - old_val); - r.m_value += delta; - SASSERT(invariant(row_id, r)); - } - } - - - void model_based_opt::update_values(unsigned_vector const& bound_vars, unsigned_vector const& bound_trail) { - for (unsigned i = bound_trail.size(); i-- > 0; ) { - unsigned x = bound_vars[i]; - row& r = m_rows[bound_trail[i]]; - rational val = r.m_coeff; - rational old_x_val = m_var2value[x]; - rational new_x_val; - rational x_coeff, eps(0); - vector const& vars = r.m_vars; - for (var const& v : vars) { - if (x == v.m_id) { - x_coeff = v.m_coeff; - } - else { - val += m_var2value[v.m_id]*v.m_coeff; - } - } - SASSERT(!x_coeff.is_zero()); - new_x_val = -val/x_coeff; - - if (r.m_type == t_lt) { - eps = abs(old_x_val - new_x_val)/rational(2); - eps = std::min(rational::one(), eps); - SASSERT(!eps.is_zero()); - - // - // ax + t < 0 - // <=> x < -t/a - // <=> x := -t/a - epsilon - // - if (x_coeff.is_pos()) { - new_x_val -= eps; - } - // - // -ax + t < 0 - // <=> -ax < -t - // <=> -x < -t/a - // <=> x > t/a - // <=> x := t/a + epsilon - // - else { - new_x_val += eps; - } - } - TRACE("opt", display(tout << "v" << x - << " coeff_x: " << x_coeff - << " old_x_val: " << old_x_val - << " new_x_val: " << new_x_val - << " eps: " << eps << " ", r); ); - m_var2value[x] = new_x_val; - - r.m_value = eval(r); - SASSERT(invariant(bound_trail[i], r)); - } - - // update and check bounds for all other affected rows. - for (unsigned i = bound_trail.size(); i-- > 0; ) { - unsigned x = bound_vars[i]; - unsigned_vector const& row_ids = m_var2row_ids[x]; - for (unsigned row_id : row_ids) { - row & r = m_rows[row_id]; - r.m_value = eval(r); - SASSERT(invariant(row_id, r)); - } - } - SASSERT(invariant()); - } - - bool model_based_opt::find_bound(unsigned x, unsigned& bound_row_index, rational& bound_coeff, bool is_pos) { - bound_row_index = UINT_MAX; - rational lub_val; - rational const& x_val = m_var2value[x]; - unsigned_vector const& row_ids = m_var2row_ids[x]; - uint_set visited; - m_above.reset(); - m_below.reset(); - for (unsigned row_id : row_ids) { - SASSERT(row_id != m_objective_id); - if (visited.contains(row_id)) { - continue; - } - visited.insert(row_id); - row& r = m_rows[row_id]; - if (r.m_alive) { - rational a = get_coefficient(row_id, x); - if (a.is_zero()) { - // skip - } - else if (a.is_pos() == is_pos || r.m_type == t_eq) { - rational value = x_val - (r.m_value/a); - if (bound_row_index == UINT_MAX) { - lub_val = value; - bound_row_index = row_id; - bound_coeff = a; - } - else if ((value == lub_val && r.m_type == opt::t_lt) || - (is_pos && value < lub_val) || - - (!is_pos && value > lub_val)) { - m_above.push_back(bound_row_index); - lub_val = value; - bound_row_index = row_id; - bound_coeff = a; - } - else { - m_above.push_back(row_id); - } - } - else { - m_below.push_back(row_id); - } - } - } - return bound_row_index != UINT_MAX; - } - - void model_based_opt::retire_row(unsigned row_id) { - m_rows[row_id].m_alive = false; - m_retired_rows.push_back(row_id); - } - - rational model_based_opt::eval(unsigned x) const { - return m_var2value[x]; - } - - rational model_based_opt::eval(def const& d) const { - vector const& vars = d.m_vars; - rational val = d.m_coeff; - for (var const& v : vars) { - val += v.m_coeff * eval(v.m_id); - } - val /= d.m_div; - return val; - } - - rational model_based_opt::eval(row const& r) const { - vector const& vars = r.m_vars; - rational val = r.m_coeff; - for (var const& v : vars) { - val += v.m_coeff * eval(v.m_id); - } - return val; - } - - rational model_based_opt::eval(vector const& coeffs) const { - rational val(0); - for (var const& v : coeffs) - val += v.m_coeff * eval(v.m_id); - return val; - } - - rational model_based_opt::get_coefficient(unsigned row_id, unsigned var_id) const { - return m_rows[row_id].get_coefficient(var_id); - } - - rational model_based_opt::row::get_coefficient(unsigned var_id) const { - if (m_vars.empty()) { - return rational::zero(); - } - unsigned lo = 0, hi = m_vars.size(); - while (lo < hi) { - unsigned mid = lo + (hi - lo)/2; - SASSERT(mid < hi); - unsigned id = m_vars[mid].m_id; - if (id == var_id) { - lo = mid; - break; - } - if (id < var_id) { - lo = mid + 1; - } - else { - hi = mid; - } - } - if (lo == m_vars.size()) { - return rational::zero(); - } - unsigned id = m_vars[lo].m_id; - if (id == var_id) { - return m_vars[lo].m_coeff; - } - else { - return rational::zero(); - } - } - - model_based_opt::row& model_based_opt::row::normalize() { -#if 0 - if (m_type == t_mod) - return *this; - rational D(denominator(abs(m_coeff))); - if (D == 0) - D = 1; - for (auto const& [id, coeff] : m_vars) - if (coeff != 0) - D = lcm(D, denominator(abs(coeff))); - if (D == 1) - return *this; - SASSERT(D > 0); - for (auto & [id, coeff] : m_vars) - coeff *= D; - m_coeff *= D; -#endif - return *this; - } - - // - // Let - // row1: t1 + a1*x <= 0 - // row2: t2 + a2*x <= 0 - // - // assume a1, a2 have the same signs: - // (t2 + a2*x) <= (t1 + a1*x)*a2/a1 - // <=> t2*a1/a2 - t1 <= 0 - // <=> t2 - t1*a2/a1 <= 0 - // - // assume a1 > 0, -a2 < 0: - // t1 + a1*x <= 0, t2 - a2*x <= 0 - // t2/a2 <= -t1/a1 - // t2 + t1*a2/a1 <= 0 - // assume -a1 < 0, a2 > 0: - // t1 - a1*x <= 0, t2 + a2*x <= 0 - // t1/a1 <= -t2/a2 - // t2 + t1*a2/a1 <= 0 - // - // the resolvent is the same in all cases (simpler proof should exist) - // - - void model_based_opt::resolve(unsigned row_src, rational const& a1, unsigned row_dst, unsigned x) { - - SASSERT(a1 == get_coefficient(row_src, x)); - SASSERT(!a1.is_zero()); - SASSERT(row_src != row_dst); - - if (m_rows[row_dst].m_alive) { - rational a2 = get_coefficient(row_dst, x); - if (is_int(x)) { - TRACE("opt", - tout << a1 << " " << a2 << ": "; - display(tout, m_rows[row_dst]); - display(tout, m_rows[row_src]);); - if (a1.is_pos() != a2.is_pos() || m_rows[row_src].m_type == opt::t_eq) { - mul_add(x, a1, row_src, a2, row_dst); - } - else { - mul(row_dst, abs(a1)); - mul_add(false, row_dst, -abs(a2), row_src); - } - TRACE("opt", display(tout, m_rows[row_dst]);); - normalize(row_dst); - } - else { - mul_add(row_dst != m_objective_id && a1.is_pos() == a2.is_pos(), row_dst, -a2/a1, row_src); - } - } - } - - /** - * a1 > 0 - * a1*x + r1 = value - * a2*x + r2 <= 0 - * ------------------ - * a1*r2 - a2*r1 <= value - */ - void model_based_opt::solve(unsigned row_src, rational const& a1, unsigned row_dst, unsigned x) { - SASSERT(a1 == get_coefficient(row_src, x)); - SASSERT(a1.is_pos()); - SASSERT(row_src != row_dst); - if (!m_rows[row_dst].m_alive) return; - rational a2 = get_coefficient(row_dst, x); - mul(row_dst, a1); - mul_add(false, row_dst, -a2, row_src); - SASSERT(get_coefficient(row_dst, x).is_zero()); - } - - // resolution for integer rows. - void model_based_opt::mul_add( - unsigned x, rational const& src_c, unsigned row_src, rational const& dst_c, unsigned row_dst) { - row& dst = m_rows[row_dst]; - row const& src = m_rows[row_src]; - SASSERT(is_int(x)); - SASSERT(t_le == dst.m_type && t_le == src.m_type); - SASSERT(src_c.is_int()); - SASSERT(dst_c.is_int()); - SASSERT(m_var2value[x].is_int()); - - rational abs_src_c = abs(src_c); - rational abs_dst_c = abs(dst_c); - rational x_val = m_var2value[x]; - rational slack = (abs_src_c - rational::one()) * (abs_dst_c - rational::one()); - rational dst_val = dst.m_value - x_val*dst_c; - rational src_val = src.m_value - x_val*src_c; - rational distance = abs_src_c * dst_val + abs_dst_c * src_val + slack; - bool use_case1 = distance.is_nonpos() || abs_src_c.is_one() || abs_dst_c.is_one(); - -#if 0 - if (distance.is_nonpos() && !abs_src_c.is_one() && !abs_dst_c.is_one()) { - unsigned r = copy_row(row_src); - mul_add(false, r, rational::one(), row_dst); - del_var(r, x); - add(r, slack); - TRACE("qe", tout << m_rows[r];); - SASSERT(!m_rows[r].m_value.is_pos()); - } -#endif - - if (use_case1) { - TRACE("opt", tout << "slack: " << slack << " " << src_c << " " << dst_val << " " << dst_c << " " << src_val << "\n";); - // dst <- abs_src_c*dst + abs_dst_c*src + slack - mul(row_dst, abs_src_c); - add(row_dst, slack); - mul_add(false, row_dst, abs_dst_c, row_src); - return; - } - - // - // create finite disjunction for |b|. - // exists x, z in [0 .. |b|-2] . b*x + s + z = 0 && ax + t <= 0 && bx + s <= 0 - // <=> - // exists x, z in [0 .. |b|-2] . b*x = -z - s && ax + t <= 0 && bx + s <= 0 - // <=> - // exists x, z in [0 .. |b|-2] . b*x = -z - s && a|b|x + |b|t <= 0 && bx + s <= 0 - // <=> - // exists x, z in [0 .. |b|-2] . b*x = -z - s && a|b|x + |b|t <= 0 && -z - s + s <= 0 - // <=> - // exists x, z in [0 .. |b|-2] . b*x = -z - s && a|b|x + |b|t <= 0 && -z <= 0 - // <=> - // exists x, z in [0 .. |b|-2] . b*x = -z - s && a|b|x + |b|t <= 0 - // <=> - // exists x, z in [0 .. |b|-2] . b*x = -z - s && a*n_sign(b)(s + z) + |b|t <= 0 - // <=> - // exists z in [0 .. |b|-2] . |b| | (z + s) && a*n_sign(b)(s + z) + |b|t <= 0 - // - - TRACE("qe", tout << "finite disjunction " << distance << " " << src_c << " " << dst_c << "\n";); - vector coeffs; - if (abs_dst_c <= abs_src_c) { - rational z = mod(dst_val, abs_dst_c); - if (!z.is_zero()) z = abs_dst_c - z; - mk_coeffs_without(coeffs, dst.m_vars, x); - add_divides(coeffs, dst.m_coeff + z, abs_dst_c); - add(row_dst, z); - mul(row_dst, src_c * n_sign(dst_c)); - mul_add(false, row_dst, abs_dst_c, row_src); - } - else { - // z := b - (s + bx) mod b - // := b - s mod b - // b | s + z <=> b | s + b - s mod b <=> b | s - s mod b - rational z = mod(src_val, abs_src_c); - if (!z.is_zero()) z = abs_src_c - z; - mk_coeffs_without(coeffs, src.m_vars, x); - add_divides(coeffs, src.m_coeff + z, abs_src_c); - mul(row_dst, abs_src_c); - add(row_dst, z * dst_c * n_sign(src_c)); - mul_add(false, row_dst, dst_c * n_sign(src_c), row_src); - } - } - - void model_based_opt::mk_coeffs_without(vector& dst, vector const& src, unsigned x) { - for (var const & v : src) { - if (v.m_id != x) dst.push_back(v); - } - } - - rational model_based_opt::n_sign(rational const& b) const { - return rational(b.is_pos()?-1:1); - } - - void model_based_opt::mul(unsigned dst, rational const& c) { - if (c.is_one()) return; - row& r = m_rows[dst]; - for (auto & v : r.m_vars) { - v.m_coeff *= c; - } - r.m_coeff *= c; - r.m_value *= c; - } - - void model_based_opt::add(unsigned dst, rational const& c) { - row& r = m_rows[dst]; - r.m_coeff += c; - r.m_value += c; - } - - void model_based_opt::sub(unsigned dst, rational const& c) { - row& r = m_rows[dst]; - r.m_coeff -= c; - r.m_value -= c; - } - - void model_based_opt::del_var(unsigned dst, unsigned x) { - row& r = m_rows[dst]; - unsigned j = 0; - for (var & v : r.m_vars) { - if (v.m_id == x) { - r.m_value -= eval(x)*r.m_coeff; - } - else { - r.m_vars[j++] = v; - } - } - r.m_vars.shrink(j); - } - - - void model_based_opt::normalize(unsigned row_id) { - row& r = m_rows[row_id]; - if (r.m_vars.empty()) { - retire_row(row_id); - return; - } - if (r.m_type == t_mod) return; - rational g(abs(r.m_vars[0].m_coeff)); - bool all_int = g.is_int(); - for (unsigned i = 1; all_int && !g.is_one() && i < r.m_vars.size(); ++i) { - rational const& coeff = r.m_vars[i].m_coeff; - if (coeff.is_int()) { - g = gcd(g, abs(coeff)); - } - else { - all_int = false; - } - } - if (all_int && !r.m_coeff.is_zero()) { - if (r.m_coeff.is_int()) { - g = gcd(g, abs(r.m_coeff)); - } - else { - all_int = false; - } - } - if (all_int && !g.is_one()) { - SASSERT(!g.is_zero()); - mul(row_id, rational::one()/g); - } - } - - // - // set row1 <- row1 + c*row2 - // - void model_based_opt::mul_add(bool same_sign, unsigned row_id1, rational const& c, unsigned row_id2) { - if (c.is_zero()) { - return; - } - - m_new_vars.reset(); - row& r1 = m_rows[row_id1]; - row const& r2 = m_rows[row_id2]; - unsigned i = 0, j = 0; - while (i < r1.m_vars.size() || j < r2.m_vars.size()) { - if (j == r2.m_vars.size()) { - m_new_vars.append(r1.m_vars.size() - i, r1.m_vars.data() + i); - break; - } - if (i == r1.m_vars.size()) { - for (; j < r2.m_vars.size(); ++j) { - m_new_vars.push_back(r2.m_vars[j]); - m_new_vars.back().m_coeff *= c; - if (row_id1 != m_objective_id) { - m_var2row_ids[r2.m_vars[j].m_id].push_back(row_id1); - } - } - break; - } - - unsigned v1 = r1.m_vars[i].m_id; - unsigned v2 = r2.m_vars[j].m_id; - if (v1 == v2) { - m_new_vars.push_back(r1.m_vars[i]); - m_new_vars.back().m_coeff += c*r2.m_vars[j].m_coeff; - ++i; - ++j; - if (m_new_vars.back().m_coeff.is_zero()) { - m_new_vars.pop_back(); - } - } - else if (v1 < v2) { - m_new_vars.push_back(r1.m_vars[i]); - ++i; - } - else { - m_new_vars.push_back(r2.m_vars[j]); - m_new_vars.back().m_coeff *= c; - if (row_id1 != m_objective_id) { - m_var2row_ids[r2.m_vars[j].m_id].push_back(row_id1); - } - ++j; - } - } - r1.m_coeff += c*r2.m_coeff; - r1.m_vars.swap(m_new_vars); - r1.m_value += c*r2.m_value; - - if (!same_sign && r2.m_type == t_lt) { - r1.m_type = t_lt; - } - else if (same_sign && r1.m_type == t_lt && r2.m_type == t_lt) { - r1.m_type = t_le; - } - SASSERT(invariant(row_id1, r1)); - } - - void model_based_opt::display(std::ostream& out) const { - for (auto const& r : m_rows) { - display(out, r); - } - for (unsigned i = 0; i < m_var2row_ids.size(); ++i) { - unsigned_vector const& rows = m_var2row_ids[i]; - out << i << ": "; - for (auto const& r : rows) { - out << r << " "; - } - out << "\n"; - } - } - - void model_based_opt::display(std::ostream& out, vector const& vars, rational const& coeff) { - unsigned i = 0; - for (var const& v : vars) { - if (i > 0 && v.m_coeff.is_pos()) { - out << "+ "; - } - ++i; - if (v.m_coeff.is_one()) { - out << "v" << v.m_id << " "; - } - else { - out << v.m_coeff << "*v" << v.m_id << " "; - } - } - if (coeff.is_pos()) { - out << " + " << coeff << " "; - } - else if (coeff.is_neg()) { - out << coeff << " "; - } - } - - std::ostream& model_based_opt::display(std::ostream& out, row const& r) { - out << (r.m_alive?"a":"d") << " "; - display(out, r.m_vars, r.m_coeff); - if (r.m_type == opt::t_mod) { - out << r.m_type << " " << r.m_mod << " = 0; value: " << r.m_value << "\n"; - } - else { - out << r.m_type << " 0; value: " << r.m_value << "\n"; - } - return out; - } - - std::ostream& model_based_opt::display(std::ostream& out, def const& r) { - display(out, r.m_vars, r.m_coeff); - if (!r.m_div.is_one()) { - out << " / " << r.m_div; - } - return out; - } - - unsigned model_based_opt::add_var(rational const& value, bool is_int) { - unsigned v = m_var2value.size(); - m_var2value.push_back(value); - m_var2is_int.push_back(is_int); - SASSERT(value.is_int() || !is_int); - m_var2row_ids.push_back(unsigned_vector()); - return v; - } - - rational model_based_opt::get_value(unsigned var) { - return m_var2value[var]; - } - - void model_based_opt::set_row(unsigned row_id, vector const& coeffs, rational const& c, rational const& m, ineq_type rel) { - row& r = m_rows[row_id]; - rational val(c); - SASSERT(r.m_vars.empty()); - r.m_vars.append(coeffs.size(), coeffs.data()); - bool is_int_row = !coeffs.empty(); - std::sort(r.m_vars.begin(), r.m_vars.end(), var::compare()); - for (auto const& c : coeffs) { - val += m_var2value[c.m_id] * c.m_coeff; - SASSERT(!is_int(c.m_id) || c.m_coeff.is_int()); - is_int_row &= is_int(c.m_id); - } - r.m_alive = true; - r.m_coeff = c; - r.m_value = val; - r.m_type = rel; - r.m_mod = m; - if (is_int_row && rel == t_lt) { - r.m_type = t_le; - r.m_coeff += rational::one(); - r.m_value += rational::one(); - } - } - - unsigned model_based_opt::new_row() { - unsigned row_id = 0; - if (m_retired_rows.empty()) { - row_id = m_rows.size(); - m_rows.push_back(row()); - } - else { - row_id = m_retired_rows.back(); - m_retired_rows.pop_back(); - m_rows[row_id].reset(); - m_rows[row_id].m_alive = true; - } - return row_id; - } - - unsigned model_based_opt::copy_row(unsigned src) { - unsigned dst = new_row(); - row const& r = m_rows[src]; - set_row(dst, r.m_vars, r.m_coeff, r.m_mod, r.m_type); - for (auto const& v : r.m_vars) { - m_var2row_ids[v.m_id].push_back(dst); - } - SASSERT(invariant(dst, m_rows[dst])); - return dst; - } - - void model_based_opt::add_constraint(vector const& coeffs, rational const& c, ineq_type rel) { - add_constraint(coeffs, c, rational::zero(), rel); - } - - void model_based_opt::add_divides(vector const& coeffs, rational const& c, rational const& m) { - add_constraint(coeffs, c, m, t_mod); - } - - void model_based_opt::add_constraint(vector const& coeffs, rational const& c, rational const& m, ineq_type rel) { - unsigned row_id = new_row(); - set_row(row_id, coeffs, c, m, rel); - for (var const& coeff : coeffs) { - m_var2row_ids[coeff.m_id].push_back(row_id); - } - SASSERT(invariant(row_id, m_rows[row_id])); - } - - void model_based_opt::set_objective(vector const& coeffs, rational const& c) { - set_row(m_objective_id, coeffs, c, rational::zero(), t_le); - } - - void model_based_opt::get_live_rows(vector& rows) { - for (row & r : m_rows) { - if (r.m_alive) { - rows.push_back(r.normalize()); - } - } - } - - // - // pick glb and lub representative. - // The representative is picked such that it - // represents the fewest inequalities. - // The constraints that enforce a glb or lub are not forced. - // The constraints that separate the glb from ub or the lub from lb - // are not forced. - // In other words, suppose there are - // . N inequalities of the form t <= x - // . M inequalities of the form s >= x - // . t0 is glb among N under valuation. - // . s0 is lub among M under valuation. - // If N < M - // create the inequalities: - // t <= t0 for each t other than t0 (N-1 inequalities). - // t0 <= s for each s (M inequalities). - // If N >= M the construction is symmetric. - // - model_based_opt::def model_based_opt::project(unsigned x, bool compute_def) { - unsigned_vector& lub_rows = m_lub; - unsigned_vector& glb_rows = m_glb; - unsigned_vector& mod_rows = m_mod; - unsigned lub_index = UINT_MAX, glb_index = UINT_MAX; - bool lub_strict = false, glb_strict = false; - rational lub_val, glb_val; - rational const& x_val = m_var2value[x]; - unsigned_vector const& row_ids = m_var2row_ids[x]; - uint_set visited; - lub_rows.reset(); - glb_rows.reset(); - mod_rows.reset(); - bool lub_is_unit = true, glb_is_unit = true; - unsigned eq_row = UINT_MAX; - // select the lub and glb. - for (unsigned row_id : row_ids) { - if (visited.contains(row_id)) { - continue; - } - visited.insert(row_id); - row& r = m_rows[row_id]; - if (!r.m_alive) { - continue; - } - rational a = get_coefficient(row_id, x); - if (a.is_zero()) { - continue; - } - if (r.m_type == t_eq) { - eq_row = row_id; - continue; - } - if (r.m_type == t_mod) { - mod_rows.push_back(row_id); - } - else if (a.is_pos()) { - rational lub_value = x_val - (r.m_value/a); - if (lub_rows.empty() || - lub_value < lub_val || - (lub_value == lub_val && r.m_type == t_lt && !lub_strict)) { - lub_val = lub_value; - lub_index = row_id; - lub_strict = r.m_type == t_lt; - } - lub_rows.push_back(row_id); - lub_is_unit &= a.is_one(); - } - else { - SASSERT(a.is_neg()); - rational glb_value = x_val - (r.m_value/a); - if (glb_rows.empty() || - glb_value > glb_val || - (glb_value == glb_val && r.m_type == t_lt && !glb_strict)) { - glb_val = glb_value; - glb_index = row_id; - glb_strict = r.m_type == t_lt; - } - glb_rows.push_back(row_id); - glb_is_unit &= a.is_minus_one(); - } - } - - if (!mod_rows.empty()) { - return solve_mod(x, mod_rows, compute_def); - } - - if (eq_row != UINT_MAX) { - return solve_for(eq_row, x, compute_def); - } - - def result; - unsigned lub_size = lub_rows.size(); - unsigned glb_size = glb_rows.size(); - unsigned row_index = (lub_size <= glb_size) ? lub_index : glb_index; - - // There are only upper or only lower bounds. - if (row_index == UINT_MAX) { - if (compute_def) { - if (lub_index != UINT_MAX) { - result = solve_for(lub_index, x, true); - } - else if (glb_index != UINT_MAX) { - result = solve_for(glb_index, x, true); - } - else { - result = def() + m_var2value[x]; - } - SASSERT(eval(result) == eval(x)); - } - else { - for (unsigned row_id : lub_rows) retire_row(row_id); - for (unsigned row_id : glb_rows) retire_row(row_id); - } - return result; - } - - SASSERT(lub_index != UINT_MAX); - SASSERT(glb_index != UINT_MAX); - if (compute_def) { - if (lub_size <= glb_size) { - result = def(m_rows[lub_index], x); - } - else { - result = def(m_rows[glb_index], x); - } - } - - // The number of matching lower and upper bounds is small. - if ((lub_size <= 2 || glb_size <= 2) && - (lub_size <= 3 && glb_size <= 3) && - (!is_int(x) || lub_is_unit || glb_is_unit)) { - for (unsigned i = 0; i < lub_size; ++i) { - unsigned row_id1 = lub_rows[i]; - bool last = i + 1 == lub_size; - rational coeff = get_coefficient(row_id1, x); - for (unsigned row_id2 : glb_rows) { - if (last) { - resolve(row_id1, coeff, row_id2, x); - } - else { - unsigned row_id3 = copy_row(row_id2); - resolve(row_id1, coeff, row_id3, x); - } - } - } - for (unsigned row_id : lub_rows) retire_row(row_id); - - return result; - } - - // General case. - rational coeff = get_coefficient(row_index, x); - for (unsigned row_id : lub_rows) { - if (row_id != row_index) { - resolve(row_index, coeff, row_id, x); - } - } - for (unsigned row_id : glb_rows) { - if (row_id != row_index) { - resolve(row_index, coeff, row_id, x); - } - } - retire_row(row_index); - return result; - } - - // - // compute D and u. - // - // D = lcm(d1, d2) - // u = eval(x) mod D - // - // d1 | (a1x + t1) & d2 | (a2x + t2) - // = - // d1 | (a1(D*x' + u) + t1) & d2 | (a2(D*x' + u) + t2) - // = - // d1 | (a1*u + t1) & d2 | (a2*u + t2) - // - // x := D*x' + u - // - - model_based_opt::def model_based_opt::solve_mod(unsigned x, unsigned_vector const& mod_rows, bool compute_def) { - SASSERT(!mod_rows.empty()); - rational D(1); - for (unsigned idx : mod_rows) { - D = lcm(D, m_rows[idx].m_mod); - } - if (D.is_zero()) { - throw default_exception("modulo 0 is not defined"); - } - if (D.is_neg()) D = abs(D); - TRACE("opt1", display(tout << "lcm: " << D << " x: v" << x << " tableau\n");); - rational val_x = m_var2value[x]; - rational u = mod(val_x, D); - SASSERT(u.is_nonneg() && u < D); - for (unsigned idx : mod_rows) { - replace_var(idx, x, u); - SASSERT(invariant(idx, m_rows[idx])); - normalize(idx); - } - TRACE("opt1", display(tout << "tableau after replace x under mod\n");); - // - // update inequalities such that u is added to t and - // D is multiplied to coefficient of x. - // the interpretation of the new version of x is (x-u)/D - // - // a*x + t <= 0 - // a*(D*x' + u) + t <= 0 - // a*D*x' + a*u + t <= 0 - // - rational new_val = (val_x - u) / D; - SASSERT(new_val.is_int()); - unsigned y = add_var(new_val, true); - unsigned_vector const& row_ids = m_var2row_ids[x]; - uint_set visited; - for (unsigned row_id : row_ids) { - if (!visited.contains(row_id)) { - // x |-> D*y + u - replace_var(row_id, x, D, y, u); - visited.insert(row_id); - normalize(row_id); - } - } - TRACE("opt1", display(tout << "tableau after replace x by y := v" << y << "\n");); - def result = project(y, compute_def); - if (compute_def) { - result = (result * D) + u; - m_var2value[x] = eval(result); - } - TRACE("opt1", display(tout << "tableau after project y" << y << "\n");); - - return result; - } - - // update row with: x |-> C - void model_based_opt::replace_var(unsigned row_id, unsigned x, rational const& C) { - row& r = m_rows[row_id]; - SASSERT(!get_coefficient(row_id, x).is_zero()); - unsigned sz = r.m_vars.size(); - unsigned i = 0, j = 0; - rational coeff(0); - for (; i < sz; ++i) { - if (r.m_vars[i].m_id == x) { - coeff = r.m_vars[i].m_coeff; - } - else { - if (i != j) { - r.m_vars[j] = r.m_vars[i]; - } - ++j; - } - } - if (j != sz) { - r.m_vars.shrink(j); - } - r.m_coeff += coeff*C; - r.m_value += coeff*(C - m_var2value[x]); - } - - // update row with: x |-> A*y + B - void model_based_opt::replace_var(unsigned row_id, unsigned x, rational const& A, unsigned y, rational const& B) { - row& r = m_rows[row_id]; - rational coeff = get_coefficient(row_id, x); - if (coeff.is_zero()) return; - if (!r.m_alive) return; - replace_var(row_id, x, B); - r.m_vars.push_back(var(y, coeff*A)); - r.m_value += coeff*A*m_var2value[y]; - if (!r.m_vars.empty() && r.m_vars.back().m_id > y) { - std::sort(r.m_vars.begin(), r.m_vars.end(), var::compare()); - } - m_var2row_ids[y].push_back(row_id); - SASSERT(invariant(row_id, r)); - } - - // 3x + t = 0 & 7 | (c*x + s) & ax <= u - // 3 | -t & 21 | (-ct + 3s) & a-t <= 3u - - model_based_opt::def model_based_opt::solve_for(unsigned row_id1, unsigned x, bool compute_def) { - TRACE("opt", tout << "v" << x << " := " << eval(x) << "\n" << m_rows[row_id1] << "\n";); - rational a = get_coefficient(row_id1, x), b; - row& r1 = m_rows[row_id1]; - ineq_type ty = r1.m_type; - SASSERT(!a.is_zero()); - SASSERT(r1.m_alive); - if (a.is_neg()) { - a.neg(); - r1.neg(); - } - SASSERT(a.is_pos()); - if (ty == t_lt) { - SASSERT(compute_def); - r1.m_coeff -= r1.m_value; - r1.m_type = t_le; - r1.m_value = 0; - } - - if (m_var2is_int[x] && !a.is_one()) { - r1.m_coeff -= r1.m_value; - r1.m_value = 0; - vector coeffs; - mk_coeffs_without(coeffs, r1.m_vars, x); - rational c = mod(-eval(coeffs), a); - add_divides(coeffs, c, a); - } - unsigned_vector const& row_ids = m_var2row_ids[x]; - uint_set visited; - visited.insert(row_id1); - for (unsigned row_id2 : row_ids) { - if (!visited.contains(row_id2)) { - visited.insert(row_id2); - b = get_coefficient(row_id2, x); - if (b.is_zero()) - continue; - row& dst = m_rows[row_id2]; - switch (dst.m_type) { - case t_eq: - case t_lt: - case t_le: - solve(row_id1, a, row_id2, x); - break; - case t_mod: - // mod reduction already done. - UNREACHABLE(); - break; - } - } - } - def result; - if (compute_def) { - result = def(m_rows[row_id1], x); - m_var2value[x] = eval(result); - TRACE("opt1", tout << "updated eval " << x << " := " << eval(x) << "\n";); - } - retire_row(row_id1); - return result; - } - - vector model_based_opt::project(unsigned num_vars, unsigned const* vars, bool compute_def) { - vector result; - for (unsigned i = 0; i < num_vars; ++i) { - result.push_back(project(vars[i], compute_def)); - TRACE("opt", display(tout << "After projecting: v" << vars[i] << "\n");); - } - return result; - } - -} - +/*++ +Copyright (c) 2016 Microsoft Corporation + +Module Name: + + model_based_opt.cpp + +Abstract: + + Model-based optimization and projection for linear real, integer arithmetic. + +Author: + + Nikolaj Bjorner (nbjorner) 2016-27-4 + +Revision History: + + +--*/ + +#include "math/simplex/model_based_opt.h" +#include "util/uint_set.h" +#include "util/z3_exception.h" + +std::ostream& operator<<(std::ostream& out, opt::ineq_type ie) { + switch (ie) { + case opt::t_eq: return out << " = "; + case opt::t_lt: return out << " < "; + case opt::t_le: return out << " <= "; + case opt::t_divides: return out << " divides "; + case opt::t_mod: return out << " mod "; + case opt::t_div: return out << " div "; + } + return out; +} + + +namespace opt { + + /** + * Convert a row ax + coeffs + coeff = value into a definition for x + * x = (value - coeffs - coeff)/a + * as backdrop we have existing assignments to x and other variables that + * satisfy the equality with value, and such that value satisfies + * the row constraint ( = , <= , < , mod) + */ + model_based_opt::def::def(row const& r, unsigned x) { + for (var const & v : r.m_vars) { + if (v.m_id != x) { + m_vars.push_back(v); + } + else { + m_div = -v.m_coeff; + } + } + m_coeff = r.m_coeff; + switch (r.m_type) { + case opt::t_lt: + m_coeff += m_div; + break; + case opt::t_le: + // for: ax >= t, then x := (t + a - 1) div a + if (m_div.is_pos()) { + m_coeff += m_div; + m_coeff -= rational::one(); + } + break; + default: + break; + } + normalize(); + SASSERT(m_div.is_pos()); + } + + model_based_opt::def model_based_opt::def::operator+(def const& other) const { + def result; + vector const& vs1 = m_vars; + vector const& vs2 = other.m_vars; + vector & vs = result.m_vars; + rational c1(1), c2(1); + if (m_div != other.m_div) { + c1 = other.m_div; + c2 = m_div; + } + unsigned i = 0, j = 0; + while (i < vs1.size() || j < vs2.size()) { + unsigned v1 = UINT_MAX, v2 = UINT_MAX; + if (i < vs1.size()) v1 = vs1[i].m_id; + if (j < vs2.size()) v2 = vs2[j].m_id; + if (v1 == v2) { + vs.push_back(vs1[i]); + vs.back().m_coeff *= c1; + vs.back().m_coeff += c2 * vs2[j].m_coeff; + ++i; ++j; + if (vs.back().m_coeff.is_zero()) { + vs.pop_back(); + } + } + else if (v1 < v2) { + vs.push_back(vs1[i]); + vs.back().m_coeff *= c1; + } + else { + vs.push_back(vs2[j]); + vs.back().m_coeff *= c2; + } + } + result.m_div = c1*m_div; + result.m_coeff = (m_coeff*c1) + (other.m_coeff*c2); + result.normalize(); + return result; + } + + model_based_opt::def model_based_opt::def::operator/(rational const& r) const { + def result(*this); + result.m_div *= r; + result.normalize(); + return result; + } + + model_based_opt::def model_based_opt::def::operator*(rational const& n) const { + def result(*this); + for (var& v : result.m_vars) { + v.m_coeff *= n; + } + result.m_coeff *= n; + result.normalize(); + return result; + } + + model_based_opt::def model_based_opt::def::operator+(rational const& n) const { + def result(*this); + result.m_coeff += n * result.m_div; + result.normalize(); + return result; + } + + void model_based_opt::def::normalize() { + if (!m_div.is_int()) { + rational den = denominator(m_div); + SASSERT(den > 1); + for (var& v : m_vars) + v.m_coeff *= den; + m_coeff *= den; + m_div *= den; + + } + if (m_div.is_neg()) { + for (var& v : m_vars) + v.m_coeff.neg(); + m_coeff.neg(); + m_div.neg(); + } + if (m_div.is_one()) + return; + rational g(m_div); + if (!m_coeff.is_int()) + return; + g = gcd(g, m_coeff); + for (var const& v : m_vars) { + if (!v.m_coeff.is_int()) + return; + g = gcd(g, abs(v.m_coeff)); + if (g.is_one()) + break; + } + if (!g.is_one()) { + for (var& v : m_vars) + v.m_coeff /= g; + m_coeff /= g; + m_div /= g; + } + } + + model_based_opt::model_based_opt() { + m_rows.push_back(row()); + } + + bool model_based_opt::invariant() { + for (unsigned i = 0; i < m_rows.size(); ++i) { + if (!invariant(i, m_rows[i])) { + return false; + } + } + return true; + } + +#define PASSERT(_e_) { CTRACE("qe", !(_e_), display(tout, r); display(tout);); SASSERT(_e_); } + + bool model_based_opt::invariant(unsigned index, row const& r) { + vector const& vars = r.m_vars; + for (unsigned i = 0; i < vars.size(); ++i) { + // variables in each row are sorted and have non-zero coefficients + PASSERT(i + 1 == vars.size() || vars[i].m_id < vars[i+1].m_id); + PASSERT(!vars[i].m_coeff.is_zero()); + PASSERT(index == 0 || m_var2row_ids[vars[i].m_id].contains(index)); + } + + PASSERT(r.m_value == eval(r)); + PASSERT(r.m_type != t_eq || r.m_value.is_zero()); + // values satisfy constraints + PASSERT(index == 0 || r.m_type != t_lt || r.m_value.is_neg()); + PASSERT(index == 0 || r.m_type != t_le || !r.m_value.is_pos()); + PASSERT(index == 0 || r.m_type != t_divides || (mod(r.m_value, r.m_mod).is_zero())); + PASSERT(index == 0 || r.m_type != t_mod || r.m_id < m_var2value.size()); + PASSERT(index == 0 || r.m_type != t_div || r.m_id < m_var2value.size()); + return true; + } + + // a1*x + obj + // a2*x + t2 <= 0 + // a3*x + t3 <= 0 + // a4*x + t4 <= 0 + // a1 > 0, a2 > 0, a3 > 0, a4 < 0 + // x <= -t2/a2 + // x <= -t2/a3 + // determine lub among these. + // then resolve lub with others + // e.g., -t2/a2 <= -t3/a3, then + // replace inequality a3*x + t3 <= 0 by -t2/a2 + t3/a3 <= 0 + // mark a4 as invalid. + // + + // a1 < 0, a2 < 0, a3 < 0, a4 > 0 + // x >= t2/a2 + // x >= t3/a3 + // determine glb among these + // the resolve glb with others. + // e.g. t2/a2 >= t3/a3 + // then replace a3*x + t3 by t3/a3 - t2/a2 <= 0 + // + inf_eps model_based_opt::maximize() { + SASSERT(invariant()); + unsigned_vector bound_trail, bound_vars; + TRACE("opt", display(tout << "tableau\n");); + while (!objective().m_vars.empty()) { + var v = objective().m_vars.back(); + unsigned x = v.m_id; + rational const& coeff = v.m_coeff; + unsigned bound_row_index; + rational bound_coeff; + if (find_bound(x, bound_row_index, bound_coeff, coeff.is_pos())) { + SASSERT(!bound_coeff.is_zero()); + TRACE("opt", display(tout << "update: " << v << " ", objective()); + for (unsigned above : m_above) { + display(tout << "resolve: ", m_rows[above]); + }); + for (unsigned above : m_above) { + resolve(bound_row_index, bound_coeff, above, x); + } + for (unsigned below : m_below) { + resolve(bound_row_index, bound_coeff, below, x); + } + // coeff*x + objective <= ub + // a2*x + t2 <= 0 + // => coeff*x <= -t2*coeff/a2 + // objective + t2*coeff/a2 <= ub + + mul_add(false, m_objective_id, - coeff/bound_coeff, bound_row_index); + retire_row(bound_row_index); + bound_trail.push_back(bound_row_index); + bound_vars.push_back(x); + } + else { + TRACE("opt", display(tout << "unbound: " << v << " ", objective());); + update_values(bound_vars, bound_trail); + return inf_eps::infinity(); + } + } + + // + // update the evaluation of variables to satisfy the bound. + // + + update_values(bound_vars, bound_trail); + + rational value = objective().m_value; + if (objective().m_type == t_lt) { + return inf_eps(inf_rational(value, rational(-1))); + } + else { + return inf_eps(inf_rational(value)); + } + } + + + void model_based_opt::update_value(unsigned x, rational const& val) { + rational old_val = m_var2value[x]; + m_var2value[x] = val; + SASSERT(val.is_int() || !is_int(x)); + unsigned_vector const& row_ids = m_var2row_ids[x]; + for (unsigned row_id : row_ids) { + rational coeff = get_coefficient(row_id, x); + if (coeff.is_zero()) { + continue; + } + row & r = m_rows[row_id]; + rational delta = coeff * (val - old_val); + r.m_value += delta; + SASSERT(invariant(row_id, r)); + } + } + + + void model_based_opt::update_values(unsigned_vector const& bound_vars, unsigned_vector const& bound_trail) { + for (unsigned i = bound_trail.size(); i-- > 0; ) { + unsigned x = bound_vars[i]; + row& r = m_rows[bound_trail[i]]; + rational val = r.m_coeff; + rational old_x_val = m_var2value[x]; + rational new_x_val; + rational x_coeff, eps(0); + vector const& vars = r.m_vars; + for (var const& v : vars) { + if (x == v.m_id) { + x_coeff = v.m_coeff; + } + else { + val += m_var2value[v.m_id]*v.m_coeff; + } + } + SASSERT(!x_coeff.is_zero()); + new_x_val = -val/x_coeff; + + if (r.m_type == t_lt) { + eps = abs(old_x_val - new_x_val)/rational(2); + eps = std::min(rational::one(), eps); + SASSERT(!eps.is_zero()); + + // + // ax + t < 0 + // <=> x < -t/a + // <=> x := -t/a - epsilon + // + if (x_coeff.is_pos()) { + new_x_val -= eps; + } + // + // -ax + t < 0 + // <=> -ax < -t + // <=> -x < -t/a + // <=> x > t/a + // <=> x := t/a + epsilon + // + else { + new_x_val += eps; + } + } + TRACE("opt", display(tout << "v" << x + << " coeff_x: " << x_coeff + << " old_x_val: " << old_x_val + << " new_x_val: " << new_x_val + << " eps: " << eps << " ", r); ); + m_var2value[x] = new_x_val; + + r.m_value = eval(r); + SASSERT(invariant(bound_trail[i], r)); + } + + // update and check bounds for all other affected rows. + for (unsigned i = bound_trail.size(); i-- > 0; ) { + unsigned x = bound_vars[i]; + unsigned_vector const& row_ids = m_var2row_ids[x]; + for (unsigned row_id : row_ids) { + row & r = m_rows[row_id]; + r.m_value = eval(r); + SASSERT(invariant(row_id, r)); + } + } + SASSERT(invariant()); + } + + bool model_based_opt::find_bound(unsigned x, unsigned& bound_row_index, rational& bound_coeff, bool is_pos) { + bound_row_index = UINT_MAX; + rational lub_val; + rational const& x_val = m_var2value[x]; + unsigned_vector const& row_ids = m_var2row_ids[x]; + uint_set visited; + m_above.reset(); + m_below.reset(); + for (unsigned row_id : row_ids) { + SASSERT(row_id != m_objective_id); + if (visited.contains(row_id)) { + continue; + } + visited.insert(row_id); + row& r = m_rows[row_id]; + if (r.m_alive) { + rational a = get_coefficient(row_id, x); + if (a.is_zero()) { + // skip + } + else if (a.is_pos() == is_pos || r.m_type == t_eq) { + rational value = x_val - (r.m_value/a); + if (bound_row_index == UINT_MAX) { + lub_val = value; + bound_row_index = row_id; + bound_coeff = a; + } + else if ((value == lub_val && r.m_type == opt::t_lt) || + (is_pos && value < lub_val) || + + (!is_pos && value > lub_val)) { + m_above.push_back(bound_row_index); + lub_val = value; + bound_row_index = row_id; + bound_coeff = a; + } + else { + m_above.push_back(row_id); + } + } + else { + m_below.push_back(row_id); + } + } + } + return bound_row_index != UINT_MAX; + } + + void model_based_opt::retire_row(unsigned row_id) { + m_rows[row_id].m_alive = false; + m_retired_rows.push_back(row_id); + } + + rational model_based_opt::eval(unsigned x) const { + return m_var2value[x]; + } + + rational model_based_opt::eval(def const& d) const { + vector const& vars = d.m_vars; + rational val = d.m_coeff; + for (var const& v : vars) { + val += v.m_coeff * eval(v.m_id); + } + val /= d.m_div; + return val; + } + + rational model_based_opt::eval(row const& r) const { + vector const& vars = r.m_vars; + rational val = r.m_coeff; + for (var const& v : vars) { + val += v.m_coeff * eval(v.m_id); + } + return val; + } + + rational model_based_opt::eval(vector const& coeffs) const { + rational val(0); + for (var const& v : coeffs) + val += v.m_coeff * eval(v.m_id); + return val; + } + + rational model_based_opt::get_coefficient(unsigned row_id, unsigned var_id) const { + return m_rows[row_id].get_coefficient(var_id); + } + + rational model_based_opt::row::get_coefficient(unsigned var_id) const { + if (m_vars.empty()) + return rational::zero(); + unsigned lo = 0, hi = m_vars.size(); + while (lo < hi) { + unsigned mid = lo + (hi - lo)/2; + SASSERT(mid < hi); + unsigned id = m_vars[mid].m_id; + if (id == var_id) { + lo = mid; + break; + } + if (id < var_id) + lo = mid + 1; + else + hi = mid; + } + if (lo == m_vars.size()) + return rational::zero(); + unsigned id = m_vars[lo].m_id; + if (id == var_id) + return m_vars[lo].m_coeff; + else + return rational::zero(); + } + + model_based_opt::row& model_based_opt::row::normalize() { +#if 0 + if (m_type == t_divides || m_type == t_mod || m_type == t_div) + return *this; + rational D(denominator(abs(m_coeff))); + if (D == 0) + D = 1; + for (auto const& [id, coeff] : m_vars) + if (coeff != 0) + D = lcm(D, denominator(abs(coeff))); + if (D == 1) + return *this; + SASSERT(D > 0); + for (auto & [id, coeff] : m_vars) + coeff *= D; + m_coeff *= D; +#endif + return *this; + } + + // + // Let + // row1: t1 + a1*x <= 0 + // row2: t2 + a2*x <= 0 + // + // assume a1, a2 have the same signs: + // (t2 + a2*x) <= (t1 + a1*x)*a2/a1 + // <=> t2*a1/a2 - t1 <= 0 + // <=> t2 - t1*a2/a1 <= 0 + // + // assume a1 > 0, -a2 < 0: + // t1 + a1*x <= 0, t2 - a2*x <= 0 + // t2/a2 <= -t1/a1 + // t2 + t1*a2/a1 <= 0 + // assume -a1 < 0, a2 > 0: + // t1 - a1*x <= 0, t2 + a2*x <= 0 + // t1/a1 <= -t2/a2 + // t2 + t1*a2/a1 <= 0 + // + // the resolvent is the same in all cases (simpler proof should exist) + // + // assume a1 < 0, -a1 = a2: + // t1 <= a2*div(t2, a2) + // + + void model_based_opt::resolve(unsigned row_src, rational const& a1, unsigned row_dst, unsigned x) { + + SASSERT(a1 == get_coefficient(row_src, x)); + SASSERT(!a1.is_zero()); + SASSERT(row_src != row_dst); + + if (m_rows[row_dst].m_alive) { + rational a2 = get_coefficient(row_dst, x); + if (is_int(x)) { + TRACE("opt", + tout << x << ": " << a1 << " " << a2 << ": "; + display(tout, m_rows[row_dst]); + display(tout, m_rows[row_src]);); + if (a1.is_pos() != a2.is_pos() || m_rows[row_src].m_type == opt::t_eq) { + mul_add(x, a1, row_src, a2, row_dst); + } + else { + mul(row_dst, abs(a1)); + mul_add(false, row_dst, -abs(a2), row_src); + } + TRACE("opt", display(tout << "result ", m_rows[row_dst]);); + normalize(row_dst); + } + else { + mul_add(row_dst != m_objective_id && a1.is_pos() == a2.is_pos(), row_dst, -a2/a1, row_src); + } + } + } + + /** + * a1 > 0 + * a1*x + r1 = value + * a2*x + r2 <= 0 + * ------------------ + * a1*r2 - a2*r1 <= value + */ + void model_based_opt::solve(unsigned row_src, rational const& a1, unsigned row_dst, unsigned x) { + SASSERT(a1 == get_coefficient(row_src, x)); + SASSERT(a1.is_pos()); + SASSERT(row_src != row_dst); + if (!m_rows[row_dst].m_alive) return; + rational a2 = get_coefficient(row_dst, x); + mul(row_dst, a1); + mul_add(false, row_dst, -a2, row_src); + normalize(row_dst); + SASSERT(get_coefficient(row_dst, x).is_zero()); + } + + // resolution for integer rows. + void model_based_opt::mul_add( + unsigned x, rational src_c, unsigned row_src, rational dst_c, unsigned row_dst) { + row& dst = m_rows[row_dst]; + row const& src = m_rows[row_src]; + SASSERT(is_int(x)); + SASSERT(t_le == dst.m_type && t_le == src.m_type); + SASSERT(src_c.is_int()); + SASSERT(dst_c.is_int()); + SASSERT(m_var2value[x].is_int()); + + rational abs_src_c = abs(src_c); + rational abs_dst_c = abs(dst_c); + rational x_val = m_var2value[x]; + rational slack = (abs_src_c - rational::one()) * (abs_dst_c - rational::one()); + rational dst_val = dst.m_value - x_val*dst_c; + rational src_val = src.m_value - x_val*src_c; + rational distance = abs_src_c * dst_val + abs_dst_c * src_val + slack; + bool use_case1 = distance.is_nonpos() || abs_src_c.is_one() || abs_dst_c.is_one(); + bool use_case2 = false && abs_src_c == abs_dst_c && src_c.is_pos() != dst_c.is_pos() && !abs_src_c.is_one() && t_le == dst.m_type && t_le == src.m_type; + bool use_case3 = false && src_c.is_pos() != dst_c.is_pos() && t_le == dst.m_type && t_le == src.m_type; + + + if (use_case1) { + TRACE("opt", tout << "slack: " << slack << " " << src_c << " " << dst_val << " " << dst_c << " " << src_val << "\n";); + // dst <- abs_src_c*dst + abs_dst_c*src + slack + mul(row_dst, abs_src_c); + add(row_dst, slack); + mul_add(false, row_dst, abs_dst_c, row_src); + return; + } + + if (use_case2 || use_case3) { + // case2: + // x*src_c + s <= 0 + // -x*src_c + t <= 0 + // + // -src_c*div(-s, src_c) + t <= 0 + // + // Example: + // t <= 100*x <= s + // Then t <= 100*div(s, 100) + // + // case3: + // x*src_c + s <= 0 + // -x*dst_c + t <= 0 + // t <= x*dst_c, x*src_c <= -s -> + // t <= dst_c*div(-s, src_c) -> + // -dst_c*div(-s,src_c) + t <= 0 + // + + bool swapped = false; + if (src_c < 0) { + std::swap(row_src, row_dst); + std::swap(src_c, dst_c); + std::swap(abs_src_c, abs_dst_c); + swapped = true; + } + vector src_coeffs, dst_coeffs; + rational src_coeff = m_rows[row_src].m_coeff; + rational dst_coeff = m_rows[row_dst].m_coeff; + for (auto const& v : m_rows[row_src].m_vars) + if (v.m_id != x) + src_coeffs.push_back(var(v.m_id, -v.m_coeff)); + for (auto const& v : m_rows[row_dst].m_vars) + if (v.m_id != x) + dst_coeffs.push_back(v); + unsigned v = UINT_MAX; + if (src_coeffs.empty()) + dst_coeff -= abs_dst_c*div(-src_coeff, abs_src_c); + else + v = add_div(src_coeffs, -src_coeff, abs_src_c); + if (v != UINT_MAX) dst_coeffs.push_back(var(v, -abs_dst_c)); + if (swapped) + std::swap(row_src, row_dst); + retire_row(row_dst); + add_constraint(dst_coeffs, dst_coeff, t_le); + return; + } + + // + // create finite disjunction for |b|. + // exists x, z in [0 .. |b|-2] . b*x + s + z = 0 && ax + t <= 0 && bx + s <= 0 + // <=> + // exists x, z in [0 .. |b|-2] . b*x = -z - s && ax + t <= 0 && bx + s <= 0 + // <=> + // exists x, z in [0 .. |b|-2] . b*x = -z - s && a|b|x + |b|t <= 0 && bx + s <= 0 + // <=> + // exists x, z in [0 .. |b|-2] . b*x = -z - s && a|b|x + |b|t <= 0 && -z - s + s <= 0 + // <=> + // exists x, z in [0 .. |b|-2] . b*x = -z - s && a|b|x + |b|t <= 0 && -z <= 0 + // <=> + // exists x, z in [0 .. |b|-2] . b*x = -z - s && a|b|x + |b|t <= 0 + // <=> + // exists x, z in [0 .. |b|-2] . b*x = -z - s && a*n_sign(b)(s + z) + |b|t <= 0 + // <=> + // exists z in [0 .. |b|-2] . |b| | (z + s) && a*n_sign(b)(s + z) + |b|t <= 0 + // + + TRACE("qe", tout << "finite disjunction " << distance << " " << src_c << " " << dst_c << "\n";); + vector coeffs; + if (abs_dst_c <= abs_src_c) { + rational z = mod(dst_val, abs_dst_c); + if (!z.is_zero()) z = abs_dst_c - z; + mk_coeffs_without(coeffs, dst.m_vars, x); + add_divides(coeffs, dst.m_coeff + z, abs_dst_c); + add(row_dst, z); + mul(row_dst, src_c * n_sign(dst_c)); + mul_add(false, row_dst, abs_dst_c, row_src); + } + else { + // z := b - (s + bx) mod b + // := b - s mod b + // b | s + z <=> b | s + b - s mod b <=> b | s - s mod b + rational z = mod(src_val, abs_src_c); + if (!z.is_zero()) z = abs_src_c - z; + mk_coeffs_without(coeffs, src.m_vars, x); + add_divides(coeffs, src.m_coeff + z, abs_src_c); + mul(row_dst, abs_src_c); + add(row_dst, z * dst_c * n_sign(src_c)); + mul_add(false, row_dst, dst_c * n_sign(src_c), row_src); + } + } + + void model_based_opt::mk_coeffs_without(vector& dst, vector const& src, unsigned x) { + for (var const & v : src) { + if (v.m_id != x) dst.push_back(v); + } + } + + rational model_based_opt::n_sign(rational const& b) const { + return rational(b.is_pos()?-1:1); + } + + void model_based_opt::mul(unsigned dst, rational const& c) { + if (c.is_one()) + return; + row& r = m_rows[dst]; + for (auto & v : r.m_vars) + v.m_coeff *= c; + r.m_mod *= c; + r.m_coeff *= c; + if (r.m_type != t_div && r.m_type != t_mod) + r.m_value *= c; + } + + void model_based_opt::add(unsigned dst, rational const& c) { + row& r = m_rows[dst]; + r.m_coeff += c; + r.m_value += c; + } + + void model_based_opt::sub(unsigned dst, rational const& c) { + row& r = m_rows[dst]; + r.m_coeff -= c; + r.m_value -= c; + } + + void model_based_opt::normalize(unsigned row_id) { + row& r = m_rows[row_id]; + if (r.m_vars.empty()) { + retire_row(row_id); + return; + } + if (r.m_type == t_divides) + return; + if (r.m_type == t_mod) + return; + if (r.m_type == t_div) + return; + rational g(abs(r.m_vars[0].m_coeff)); + bool all_int = g.is_int(); + for (unsigned i = 1; all_int && !g.is_one() && i < r.m_vars.size(); ++i) { + rational const& coeff = r.m_vars[i].m_coeff; + if (coeff.is_int()) { + g = gcd(g, abs(coeff)); + } + else { + all_int = false; + } + } + if (all_int && !r.m_coeff.is_zero()) { + if (r.m_coeff.is_int()) { + g = gcd(g, abs(r.m_coeff)); + } + else { + all_int = false; + } + } + if (all_int && !g.is_one()) { + SASSERT(!g.is_zero()); + mul(row_id, rational::one()/g); + } + } + + // + // set row1 <- row1 + c*row2 + // + void model_based_opt::mul_add(bool same_sign, unsigned row_id1, rational const& c, unsigned row_id2) { + if (c.is_zero()) + return; + + + m_new_vars.reset(); + row& r1 = m_rows[row_id1]; + row const& r2 = m_rows[row_id2]; + unsigned i = 0, j = 0; + while (i < r1.m_vars.size() || j < r2.m_vars.size()) { + if (j == r2.m_vars.size()) { + m_new_vars.append(r1.m_vars.size() - i, r1.m_vars.data() + i); + break; + } + if (i == r1.m_vars.size()) { + for (; j < r2.m_vars.size(); ++j) { + m_new_vars.push_back(r2.m_vars[j]); + m_new_vars.back().m_coeff *= c; + if (row_id1 != m_objective_id) + m_var2row_ids[r2.m_vars[j].m_id].push_back(row_id1); + } + break; + } + + unsigned v1 = r1.m_vars[i].m_id; + unsigned v2 = r2.m_vars[j].m_id; + if (v1 == v2) { + m_new_vars.push_back(r1.m_vars[i]); + m_new_vars.back().m_coeff += c*r2.m_vars[j].m_coeff; + ++i; + ++j; + if (m_new_vars.back().m_coeff.is_zero()) + m_new_vars.pop_back(); + } + else if (v1 < v2) { + m_new_vars.push_back(r1.m_vars[i]); + ++i; + } + else { + m_new_vars.push_back(r2.m_vars[j]); + m_new_vars.back().m_coeff *= c; + if (row_id1 != m_objective_id) + m_var2row_ids[r2.m_vars[j].m_id].push_back(row_id1); + ++j; + } + } + r1.m_coeff += c*r2.m_coeff; + r1.m_vars.swap(m_new_vars); + r1.m_value += c*r2.m_value; + + if (!same_sign && r2.m_type == t_lt) + r1.m_type = t_lt; + else if (same_sign && r1.m_type == t_lt && r2.m_type == t_lt) + r1.m_type = t_le; + SASSERT(invariant(row_id1, r1)); + } + + void model_based_opt::display(std::ostream& out) const { + for (auto const& r : m_rows) + display(out, r); + for (unsigned i = 0; i < m_var2row_ids.size(); ++i) { + unsigned_vector const& rows = m_var2row_ids[i]; + out << i << ": "; + for (auto const& r : rows) + out << r << " "; + out << "\n"; + } + } + + void model_based_opt::display(std::ostream& out, vector const& vars, rational const& coeff) { + unsigned i = 0; + for (var const& v : vars) { + if (i > 0 && v.m_coeff.is_pos()) + out << "+ "; + ++i; + if (v.m_coeff.is_one()) + out << "v" << v.m_id << " "; + else + out << v.m_coeff << "*v" << v.m_id << " "; + } + if (coeff.is_pos()) + out << " + " << coeff << " "; + else if (coeff.is_neg()) + out << coeff << " "; + } + + std::ostream& model_based_opt::display(std::ostream& out, row const& r) { + out << (r.m_alive?"a":"d") << " "; + display(out, r.m_vars, r.m_coeff); + switch (r.m_type) { + case opt::t_divides: + out << r.m_type << " " << r.m_mod << " = 0; value: " << r.m_value << "\n"; + break; + case opt::t_mod: + out << r.m_type << " " << r.m_mod << " = v" << r.m_id << " ; mod: " << mod(r.m_value, r.m_mod) << "\n"; + break; + case opt::t_div: + out << r.m_type << " " << r.m_mod << " = v" << r.m_id << " ; div: " << div(r.m_value, r.m_mod) << "\n"; + break; + default: + out << r.m_type << " 0; value: " << r.m_value << "\n"; + break; + } + return out; + } + + std::ostream& model_based_opt::display(std::ostream& out, def const& r) { + display(out, r.m_vars, r.m_coeff); + if (!r.m_div.is_one()) { + out << " / " << r.m_div; + } + return out; + } + + unsigned model_based_opt::add_var(rational const& value, bool is_int) { + unsigned v = m_var2value.size(); + m_var2value.push_back(value); + m_var2is_int.push_back(is_int); + SASSERT(value.is_int() || !is_int); + m_var2row_ids.push_back(unsigned_vector()); + return v; + } + + rational model_based_opt::get_value(unsigned var) { + return m_var2value[var]; + } + + void model_based_opt::set_row(unsigned row_id, vector const& coeffs, rational const& c, rational const& m, ineq_type rel) { + row& r = m_rows[row_id]; + rational val(c); + SASSERT(r.m_vars.empty()); + r.m_vars.append(coeffs.size(), coeffs.data()); + bool is_int_row = !coeffs.empty(); + std::sort(r.m_vars.begin(), r.m_vars.end(), var::compare()); + for (auto const& c : coeffs) { + val += m_var2value[c.m_id] * c.m_coeff; + SASSERT(!is_int(c.m_id) || c.m_coeff.is_int()); + is_int_row &= is_int(c.m_id); + } + r.m_alive = true; + r.m_coeff = c; + r.m_value = val; + r.m_type = rel; + r.m_mod = m; + if (is_int_row && rel == t_lt) { + r.m_type = t_le; + r.m_coeff += rational::one(); + r.m_value += rational::one(); + } + } + + unsigned model_based_opt::new_row() { + unsigned row_id = 0; + if (m_retired_rows.empty()) { + row_id = m_rows.size(); + m_rows.push_back(row()); + } + else { + row_id = m_retired_rows.back(); + m_retired_rows.pop_back(); + m_rows[row_id].reset(); + m_rows[row_id].m_alive = true; + } + return row_id; + } + + unsigned model_based_opt::copy_row(unsigned src, unsigned excl) { + unsigned dst = new_row(); + row const& r = m_rows[src]; + set_row(dst, r.m_vars, r.m_coeff, r.m_mod, r.m_type); + for (auto const& v : r.m_vars) { + if (v.m_id != excl) + m_var2row_ids[v.m_id].push_back(dst); + } + SASSERT(invariant(dst, m_rows[dst])); + return dst; + } + + void model_based_opt::add_lower_bound(unsigned x, rational const& lo) { + vector coeffs; + coeffs.push_back(var(x, rational::minus_one())); + add_constraint(coeffs, lo, t_le); + } + + void model_based_opt::add_upper_bound(unsigned x, rational const& hi) { + vector coeffs; + coeffs.push_back(var(x, rational::one())); + add_constraint(coeffs, -hi, t_le); + } + + void model_based_opt::add_constraint(vector const& coeffs, rational const& c, ineq_type rel) { + add_constraint(coeffs, c, rational::zero(), rel, 0); + } + + void model_based_opt::add_divides(vector const& coeffs, rational const& c, rational const& m) { + rational g(c); + for (auto const& [v, coeff] : coeffs) + g = gcd(coeff, g); + if ((g/m).is_int()) + return; + add_constraint(coeffs, c, m, t_divides, 0); + } + + unsigned model_based_opt::add_mod(vector const& coeffs, rational const& c, rational const& m) { + rational value = c; + for (auto const& var : coeffs) + value += var.m_coeff * m_var2value[var.m_id]; + unsigned v = add_var(mod(value, m), true); + add_constraint(coeffs, c, m, t_mod, v); + return v; + } + + unsigned model_based_opt::add_div(vector const& coeffs, rational const& c, rational const& m) { + rational value = c; + for (auto const& var : coeffs) + value += var.m_coeff * m_var2value[var.m_id]; + unsigned v = add_var(div(value, m), true); + add_constraint(coeffs, c, m, t_div, v); + return v; + } + + void model_based_opt::add_constraint(vector const& coeffs, rational const& c, rational const& m, ineq_type rel, unsigned id) { + auto const& r = m_rows.back(); + if (r.m_vars == coeffs && r.m_coeff == c && r.m_mod == m && r.m_type == rel && r.m_id == id && r.m_alive) + return; + unsigned row_id = new_row(); + set_row(row_id, coeffs, c, m, rel); + m_rows[row_id].m_id = id; + for (var const& coeff : coeffs) + m_var2row_ids[coeff.m_id].push_back(row_id); + SASSERT(invariant(row_id, m_rows[row_id])); + normalize(row_id); + } + + void model_based_opt::set_objective(vector const& coeffs, rational const& c) { + set_row(m_objective_id, coeffs, c, rational::zero(), t_le); + } + + void model_based_opt::get_live_rows(vector& rows) { + for (row & r : m_rows) + if (r.m_alive) + rows.push_back(r.normalize()); + } + + // + // pick glb and lub representative. + // The representative is picked such that it + // represents the fewest inequalities. + // The constraints that enforce a glb or lub are not forced. + // The constraints that separate the glb from ub or the lub from lb + // are not forced. + // In other words, suppose there are + // . N inequalities of the form t <= x + // . M inequalities of the form s >= x + // . t0 is glb among N under valuation. + // . s0 is lub among M under valuation. + // If N < M + // create the inequalities: + // t <= t0 for each t other than t0 (N-1 inequalities). + // t0 <= s for each s (M inequalities). + // If N >= M the construction is symmetric. + // + model_based_opt::def model_based_opt::project(unsigned x, bool compute_def) { + unsigned_vector& lub_rows = m_lub; + unsigned_vector& glb_rows = m_glb; + unsigned_vector& divide_rows = m_divides; + unsigned_vector& mod_rows = m_mod; + unsigned_vector& div_rows = m_div; + unsigned lub_index = UINT_MAX, glb_index = UINT_MAX; + bool lub_strict = false, glb_strict = false; + rational lub_val, glb_val; + rational const& x_val = m_var2value[x]; + unsigned_vector const& row_ids = m_var2row_ids[x]; + uint_set visited; + lub_rows.reset(); + glb_rows.reset(); + divide_rows.reset(); + mod_rows.reset(); + div_rows.reset(); + bool lub_is_unit = true, glb_is_unit = true; + unsigned eq_row = UINT_MAX; + // select the lub and glb. + for (unsigned row_id : row_ids) { + if (visited.contains(row_id)) { + continue; + } + visited.insert(row_id); + row& r = m_rows[row_id]; + if (!r.m_alive) { + continue; + } + rational a = get_coefficient(row_id, x); + if (a.is_zero()) { + continue; + } + if (r.m_type == t_eq) { + eq_row = row_id; + continue; + } + if (r.m_type == t_mod) + mod_rows.push_back(row_id); + else if (r.m_type == t_div) + div_rows.push_back(row_id); + else if (r.m_type == t_divides) + divide_rows.push_back(row_id); + else if (a.is_pos()) { + rational lub_value = x_val - (r.m_value/a); + if (lub_rows.empty() || + lub_value < lub_val || + (lub_value == lub_val && r.m_type == t_lt && !lub_strict)) { + lub_val = lub_value; + lub_index = row_id; + lub_strict = r.m_type == t_lt; + } + lub_rows.push_back(row_id); + lub_is_unit &= a.is_one(); + } + else { + SASSERT(a.is_neg()); + rational glb_value = x_val - (r.m_value/a); + if (glb_rows.empty() || + glb_value > glb_val || + (glb_value == glb_val && r.m_type == t_lt && !glb_strict)) { + glb_val = glb_value; + glb_index = row_id; + glb_strict = r.m_type == t_lt; + } + glb_rows.push_back(row_id); + glb_is_unit &= a.is_minus_one(); + } + } + + if (!mod_rows.empty()) + return solve_mod(x, mod_rows, compute_def); + + if (!div_rows.empty()) + return solve_div(x, div_rows, compute_def); + + if (!divide_rows.empty()) + return solve_divides(x, divide_rows, compute_def); + + if (eq_row != UINT_MAX) + return solve_for(eq_row, x, compute_def); + + def result; + unsigned lub_size = lub_rows.size(); + unsigned glb_size = glb_rows.size(); + unsigned row_index = (lub_size <= glb_size) ? lub_index : glb_index; + + // There are only upper or only lower bounds. + if (row_index == UINT_MAX) { + if (compute_def) { + if (lub_index != UINT_MAX) + result = solve_for(lub_index, x, true); + else if (glb_index != UINT_MAX) + result = solve_for(glb_index, x, true); + else + result = def() + m_var2value[x]; + SASSERT(eval(result) == eval(x)); + } + else { + for (unsigned row_id : lub_rows) retire_row(row_id); + for (unsigned row_id : glb_rows) retire_row(row_id); + } + return result; + } + + SASSERT(lub_index != UINT_MAX); + SASSERT(glb_index != UINT_MAX); + if (compute_def) { + if (lub_size <= glb_size) + result = def(m_rows[lub_index], x); + else + result = def(m_rows[glb_index], x); + } + + // The number of matching lower and upper bounds is small. + if ((lub_size <= 2 || glb_size <= 2) && + (lub_size <= 3 && glb_size <= 3) && + (!is_int(x) || lub_is_unit || glb_is_unit)) { + for (unsigned i = 0; i < lub_size; ++i) { + unsigned row_id1 = lub_rows[i]; + bool last = i + 1 == lub_size; + rational coeff = get_coefficient(row_id1, x); + for (unsigned row_id2 : glb_rows) { + if (last) { + resolve(row_id1, coeff, row_id2, x); + } + else { + unsigned row_id3 = copy_row(row_id2); + resolve(row_id1, coeff, row_id3, x); + } + } + } + for (unsigned row_id : lub_rows) + retire_row(row_id); + + return result; + } + + // General case. + rational coeff = get_coefficient(row_index, x); + + for (unsigned row_id : lub_rows) + if (row_id != row_index) + resolve(row_index, coeff, row_id, x); + + for (unsigned row_id : glb_rows) + if (row_id != row_index) + resolve(row_index, coeff, row_id, x); + retire_row(row_index); + return result; + } + + + // + // Given v = a*x + b mod K + // + // - remove v = a*x + b mod K + // + // case a = 1: + // - add w = b mod K + // - x |-> K*y + z, 0 <= z < K + // - if z.value + w.value < K: + // add z + w - v = 0 + // - if z.value + w.value >= K: + // add z + w - v - K = 0 + // + // case a != 1, gcd(a, K) = 1 + // - x |-> x*y + a^-1*z, 0 <= z < K + // - add w = b mod K + // if z.value + w.value < K + // add z + w - v = 0 + // if z.value + w.value >= K + // add z + w - v - K = 0 + // + // case a != 1, gcd(a,K) = g != 1 + // - x |-> x*y + a^-1*z, 0 <= z < K + // a*x + b mod K = v is now + // g*z + b mod K = v + // - add w = b mod K + // - 0 <= g*z.value + w.value < K*(g+1) + // - add g*z + w - v - k*K = 0 for suitable k from 0 .. g based on model + // + + model_based_opt::def model_based_opt::solve_mod(unsigned x, unsigned_vector const& _mod_rows, bool compute_def) { + def result; + unsigned_vector mod_rows(_mod_rows); + rational K(1); + for (unsigned ri : mod_rows) + K = lcm(K, m_rows[ri].m_mod); + + rational x_value = m_var2value[x]; + rational y_value = div(x_value, K); + rational z_value = mod(x_value, K); + SASSERT(x_value == K * y_value + z_value); + SASSERT(0 <= z_value && z_value < K); + // add new variables + unsigned z = add_var(z_value, true); + unsigned y = add_var(y_value, true); + + uint_set visited; + for (unsigned ri : mod_rows) { + m_rows[ri].m_alive = false; + visited.insert(ri); + } + + // replace x by K*y + z in other rows. + for (unsigned ri : m_var2row_ids[x]) { + if (visited.contains(ri)) + continue; + replace_var(ri, x, K, y, rational::one(), z); + visited.insert(ri); + normalize(ri); + } + + // add bounds for z + add_lower_bound(z, rational::zero()); + add_upper_bound(z, K - 1); + + for (unsigned ri : mod_rows) { + rational a = get_coefficient(ri, x); + // add w = b mod K + vector coeffs = m_rows[ri].m_vars; + rational coeff = m_rows[ri].m_coeff; + unsigned v = m_rows[ri].m_id; + + unsigned w = UINT_MAX; + rational offset(0); + if (coeffs.empty()) + offset = mod(coeff, K); + else + w = add_mod(coeffs, coeff, K); + + rational w_value = w == UINT_MAX ? offset : m_var2value[w]; + + // add v = a*z + w - k*K = 0, for k = (a*z_value + w_value) div K + // claim: (= (mod x K) (- x (* K (div x K)))))) is a theorem for every x, K != 0 + rational kK = div(a * z_value + w_value, K) * K; + vector mod_coeffs; + mod_coeffs.push_back(var(v, rational::minus_one())); + mod_coeffs.push_back(var(z, a)); + if (w != UINT_MAX) mod_coeffs.push_back(var(w, rational::one())); + add_constraint(mod_coeffs, kK + offset, t_eq); + add_lower_bound(v, rational::zero()); + add_upper_bound(v, K - 1); + + // allow to recycle row. + m_retired_rows.push_back(ri); + + project(v, false); + } + + def y_def = project(y, compute_def); + def z_def = project(z, compute_def); + + if (compute_def) { + result = (y_def * K) + z_def; + m_var2value[x] = eval(result); + } + TRACE("opt", display(tout << "solve_mod\n")); + return result; + } + + // + // Given v = a*x + b div K + // Replace x |-> K*y + z + // - w = b div K + // - v = ((a*K*y + a*z) + b) div K + // = a*y + (a*z + b) div K + // = a*y + b div K + (b mod K + a*z) div K + // = a*y + b div K + k + // where k := (b.value mod K + a*z.value) div K + // k is between 0 and a + // + // - k*K <= b mod K + a*z < (k+1)*K + // + // A better version using a^-1 + // - v = (a*K*y + a^-1*a*z + b) div K + // = a*y + ((K*A + g)*z + b) div K where we write a*a^-1 = K*A + g + // = a*y + A + (g*z + b) div K + // - k*K <= b Kod m + gz < (k+1)*K + // where k is between 0 and g + // when gcd(a, K) = 1, then there are only two cases. + // + model_based_opt::def model_based_opt::solve_div(unsigned x, unsigned_vector const& _div_rows, bool compute_def) { + def result; + unsigned_vector div_rows(_div_rows); + SASSERT(!div_rows.empty()); + + rational K(1); + for (unsigned ri : div_rows) + K = lcm(K, m_rows[ri].m_mod); + + rational x_value = m_var2value[x]; + rational z_value = mod(x_value, K); + rational y_value = div(x_value, K); + SASSERT(x_value == K * y_value + z_value); + SASSERT(0 <= z_value && z_value < K); + // add new variables + unsigned z = add_var(z_value, true); + unsigned y = add_var(y_value, true); + + uint_set visited; + for (unsigned ri : div_rows) { + row& r = m_rows[ri]; + mul(ri, K / r.m_mod); + r.m_alive = false; + visited.insert(ri); + } + + // replace x by K*y + z in other rows. + for (unsigned ri : m_var2row_ids[x]) { + if (visited.contains(ri)) + continue; + replace_var(ri, x, K, y, rational::one(), z); + visited.insert(ri); + normalize(ri); + } + + // add bounds for z + add_lower_bound(z, rational::zero()); + add_upper_bound(z, K - 1); + + TRACE("opt", display(tout)); + + // solve for x_value = K*y_value + z_value, 0 <= z_value < K. + + for (unsigned ri : div_rows) { + + rational a = get_coefficient(ri, x); + replace_var(ri, x, rational::zero()); + + // add w = b div m + vector coeffs = m_rows[ri].m_vars; + rational coeff = m_rows[ri].m_coeff; + unsigned w = UINT_MAX; + rational offset(0); + if (coeffs.empty()) + offset = div(coeff, K); + else + w = add_div(coeffs, coeff, K); + + // + // w = b div K + // v = a*y + w + k + // k = (a*z_value + (b_value mod K)) div K + // k*K <= a*z + b mod K < (k+1)*K + // + /** + * It is based on the following claim (tested for select values of a, K) + * (define-const K Int 13) + * (declare-const b Int) + * (define-const a Int -11) + * (declare-const y Int) + * (declare-const z Int) + * (define-const w Int (div b K)) + * (define-const k1 Int (+ (* a z) (mod b K))) + * (define-const k Int (div k1 K)) + * (define-const x Int (+ (* K y) z)) + * (define-const u Int (+ (* a x) b)) + * (define-const v Int (+ (* a y) w k)) + * (assert (<= 0 z)) + * (assert (< z K)) + * (assert (<= (* K k) k1)) + * (assert (< k1 (* K (+ k 1)))) + * (assert (not (= (div u K) v))) + * (check-sat) + */ + unsigned v = m_rows[ri].m_id; + rational b_value = eval(coeffs) + coeff; + rational k = div(a * z_value + mod(b_value, K), K); + vector div_coeffs; + div_coeffs.push_back(var(v, rational::minus_one())); + div_coeffs.push_back(var(y, a)); + if (w != UINT_MAX) div_coeffs.push_back(var(w, rational::one())); + add_constraint(div_coeffs, k + offset, t_eq); + + unsigned u = UINT_MAX; + offset = 0; + if (coeffs.empty()) + offset = mod(coeff, K); + else + u = add_mod(coeffs, coeff, K); + + // add a*z + (b mod K) < (k + 1)*K + vector bound_coeffs; + bound_coeffs.push_back(var(z, a)); + if (u != UINT_MAX) bound_coeffs.push_back(var(u, rational::one())); + add_constraint(bound_coeffs, 1 - K * (k + 1) + offset, t_le); + + // add k*K <= az + (b mod K) + for (auto& c : bound_coeffs) + c.m_coeff.neg(); + add_constraint(bound_coeffs, k * K - offset, t_le); + // allow to recycle row. + m_retired_rows.push_back(ri); + project(v, false); + } + + TRACE("opt", display(tout << "solve_div reduced " << y << " " << z << "\n")); + // project internal variables. + + def y_def = project(y, compute_def); + def z_def = project(z, compute_def); + + if (compute_def) { + result = (y_def * K) + z_def; + m_var2value[x] = eval(result); + } + TRACE("opt", display(tout << "solve_div done v" << x << "\n")); + return result; + } + + // + // compute D and u. + // + // D = lcm(d1, d2) + // u = eval(x) mod D + // + // d1 | (a1x + t1) & d2 | (a2x + t2) + // = + // d1 | (a1(D*x' + u) + t1) & d2 | (a2(D*x' + u) + t2) + // = + // d1 | (a1*u + t1) & d2 | (a2*u + t2) + // + // x := D*x' + u + // + + model_based_opt::def model_based_opt::solve_divides(unsigned x, unsigned_vector const& divide_rows, bool compute_def) { + SASSERT(!divide_rows.empty()); + rational D(1); + for (unsigned idx : divide_rows) { + D = lcm(D, m_rows[idx].m_mod); + } + if (D.is_zero()) { + throw default_exception("modulo 0 is not defined"); + } + if (D.is_neg()) D = abs(D); + TRACE("opt1", display(tout << "lcm: " << D << " x: v" << x << " tableau\n");); + rational val_x = m_var2value[x]; + rational u = mod(val_x, D); + SASSERT(u.is_nonneg() && u < D); + for (unsigned idx : divide_rows) { + replace_var(idx, x, u); + SASSERT(invariant(idx, m_rows[idx])); + normalize(idx); + } + TRACE("opt1", display(tout << "tableau after replace x under mod\n");); + // + // update inequalities such that u is added to t and + // D is multiplied to coefficient of x. + // the interpretation of the new version of x is (x-u)/D + // + // a*x + t <= 0 + // a*(D*x' + u) + t <= 0 + // a*D*x' + a*u + t <= 0 + // + rational new_val = (val_x - u) / D; + SASSERT(new_val.is_int()); + unsigned y = add_var(new_val, true); + unsigned_vector const& row_ids = m_var2row_ids[x]; + uint_set visited; + for (unsigned row_id : row_ids) { + if (!visited.contains(row_id)) { + // x |-> D*y + u + replace_var(row_id, x, D, y, u); + visited.insert(row_id); + normalize(row_id); + } + } + TRACE("opt1", display(tout << "tableau after replace x by y := v" << y << "\n");); + def result = project(y, compute_def); + if (compute_def) { + result = (result * D) + u; + m_var2value[x] = eval(result); + } + TRACE("opt1", display(tout << "tableau after project y" << y << "\n");); + + return result; + } + + // update row with: x |-> C + void model_based_opt::replace_var(unsigned row_id, unsigned x, rational const& C) { + row& r = m_rows[row_id]; + SASSERT(!get_coefficient(row_id, x).is_zero()); + unsigned sz = r.m_vars.size(); + unsigned i = 0, j = 0; + rational coeff(0); + for (; i < sz; ++i) { + if (r.m_vars[i].m_id == x) { + coeff = r.m_vars[i].m_coeff; + } + else { + if (i != j) { + r.m_vars[j] = r.m_vars[i]; + } + ++j; + } + } + if (j != sz) { + r.m_vars.shrink(j); + } + r.m_coeff += coeff*C; + r.m_value += coeff*(C - m_var2value[x]); + } + + // update row with: x |-> A*y + B + void model_based_opt::replace_var(unsigned row_id, unsigned x, rational const& A, unsigned y, rational const& B) { + row& r = m_rows[row_id]; + rational coeff = get_coefficient(row_id, x); + if (coeff.is_zero()) return; + if (!r.m_alive) return; + replace_var(row_id, x, B); + r.m_vars.push_back(var(y, coeff*A)); + r.m_value += coeff*A*m_var2value[y]; + if (!r.m_vars.empty() && r.m_vars.back().m_id > y) + std::sort(r.m_vars.begin(), r.m_vars.end(), var::compare()); + m_var2row_ids[y].push_back(row_id); + SASSERT(invariant(row_id, r)); + } + + // update row with: x |-> A*y + B*z + void model_based_opt::replace_var(unsigned row_id, unsigned x, rational const& A, unsigned y, rational const& B, unsigned z) { + row& r = m_rows[row_id]; + rational coeff = get_coefficient(row_id, x); + if (coeff.is_zero() || !r.m_alive) + return; + replace_var(row_id, x, rational::zero()); + if (A != 0) r.m_vars.push_back(var(y, coeff*A)); + if (B != 0) r.m_vars.push_back(var(z, coeff*B)); + r.m_value += coeff*A*m_var2value[y]; + r.m_value += coeff*B*m_var2value[z]; + std::sort(r.m_vars.begin(), r.m_vars.end(), var::compare()); + if (A != 0) m_var2row_ids[y].push_back(row_id); + if (B != 0) m_var2row_ids[z].push_back(row_id); + SASSERT(invariant(row_id, r)); + } + + // 3x + t = 0 & 7 | (c*x + s) & ax <= u + // 3 | -t & 21 | (-ct + 3s) & a-t <= 3u + + model_based_opt::def model_based_opt::solve_for(unsigned row_id1, unsigned x, bool compute_def) { + TRACE("opt", tout << "v" << x << " := " << eval(x) << "\n" << m_rows[row_id1] << "\n"; + display(tout)); + rational a = get_coefficient(row_id1, x), b; + row& r1 = m_rows[row_id1]; + ineq_type ty = r1.m_type; + SASSERT(!a.is_zero()); + SASSERT(r1.m_alive); + if (a.is_neg()) { + a.neg(); + r1.neg(); + } + SASSERT(a.is_pos()); + if (ty == t_lt) { + SASSERT(compute_def); + r1.m_coeff -= r1.m_value; + r1.m_type = t_le; + r1.m_value = 0; + } + + if (m_var2is_int[x] && !a.is_one()) { + r1.m_coeff -= r1.m_value; + r1.m_value = 0; + vector coeffs; + mk_coeffs_without(coeffs, r1.m_vars, x); + rational c = mod(-eval(coeffs), a); + add_divides(coeffs, c, a); + } + unsigned_vector const& row_ids = m_var2row_ids[x]; + uint_set visited; + visited.insert(row_id1); + for (unsigned row_id2 : row_ids) { + if (!visited.contains(row_id2)) { + visited.insert(row_id2); + b = get_coefficient(row_id2, x); + if (b.is_zero()) + continue; + row& dst = m_rows[row_id2]; + switch (dst.m_type) { + case t_eq: + case t_lt: + case t_le: + solve(row_id1, a, row_id2, x); + break; + case t_divides: + // mod reduction already done. + UNREACHABLE(); + break; + } + } + } + def result; + if (compute_def) { + result = def(m_rows[row_id1], x); + m_var2value[x] = eval(result); + TRACE("opt1", tout << "updated eval " << x << " := " << eval(x) << "\n";); + } + retire_row(row_id1); + TRACE("opt", display(tout << "solved v" << x << "\n")); + return result; + } + + vector model_based_opt::project(unsigned num_vars, unsigned const* vars, bool compute_def) { + vector result; + for (unsigned i = 0; i < num_vars; ++i) { + result.push_back(project(vars[i], compute_def)); + TRACE("opt", display(tout << "After projecting: v" << vars[i] << "\n");); + } + return result; + } + +} + diff --git a/src/math/simplex/model_based_opt.h b/src/math/simplex/model_based_opt.h index cb8c18542..a62c265a9 100644 --- a/src/math/simplex/model_based_opt.h +++ b/src/math/simplex/model_based_opt.h @@ -30,7 +30,9 @@ namespace opt { t_eq, t_lt, t_le, - t_mod + t_divides, + t_mod, + t_div }; @@ -48,6 +50,14 @@ namespace opt { return x.m_id < y.m_id; } }; + + bool operator==(var const& other) const { + return m_id == other.m_id && m_coeff == other.m_coeff; + } + + bool operator!=(var const& other) const { + return !(*this == other); + } }; struct row { row(): m_type(t_le), m_value(0), m_alive(false) {} @@ -57,6 +67,7 @@ namespace opt { ineq_type m_type; // inequality type rational m_value; // value of m_vars + m_coeff under interpretation of m_var2value. bool m_alive; // rows can be marked dead if they have been processed. + unsigned m_id; // variable defined by row (used for mod_t and div_t) void reset() { m_vars.reset(); m_coeff.reset(); m_value.reset(); } row& normalize(); @@ -85,9 +96,9 @@ namespace opt { static const unsigned m_objective_id = 0; vector m_var2row_ids; vector m_var2value; - bool_vector m_var2is_int; + bool_vector m_var2is_int; vector m_new_vars; - unsigned_vector m_lub, m_glb, m_mod; + unsigned_vector m_lub, m_glb, m_divides, m_mod, m_div; unsigned_vector m_above, m_below; unsigned_vector m_retired_rows; @@ -114,7 +125,7 @@ namespace opt { void mul_add(bool same_sign, unsigned row_id1, rational const& c, unsigned row_id2); - void mul_add(unsigned x, rational const& a1, unsigned row_src, rational const& a2, unsigned row_dst); + void mul_add(unsigned x, rational a1, unsigned row_src, rational a2, unsigned row_dst); void mul(unsigned dst, rational const& c); @@ -122,14 +133,18 @@ namespace opt { void sub(unsigned dst, rational const& c); - void del_var(unsigned dst, unsigned x); + void set_row(unsigned row_id, vector const& coeffs, rational const& c, rational const& m, ineq_type rel); - void set_row(unsigned row_id, vector const& coeffs, rational const& c, rational const& m, ineq_type rel); + void add_lower_bound(unsigned x, rational const& lo); - void add_constraint(vector const& coeffs, rational const& c, rational const& m, ineq_type r); + void add_upper_bound(unsigned x, rational const& hi); + + void add_constraint(vector const& coeffs, rational const& c, rational const& m, ineq_type r, unsigned id); void replace_var(unsigned row_id, unsigned x, rational const& A, unsigned y, rational const& B); + void replace_var(unsigned row_id, unsigned x, rational const& A, unsigned y, rational const& B, unsigned z); + void replace_var(unsigned row_id, unsigned x, rational const& C); void normalize(unsigned row_id); @@ -138,7 +153,7 @@ namespace opt { unsigned new_row(); - unsigned copy_row(unsigned row_id); + unsigned copy_row(unsigned row_id, unsigned excl = UINT_MAX); rational n_sign(rational const& b) const; @@ -150,8 +165,12 @@ namespace opt { def solve_for(unsigned row_id, unsigned x, bool compute_def); - def solve_mod(unsigned x, unsigned_vector const& mod_rows, bool compute_def); - + def solve_divides(unsigned x, unsigned_vector const& divide_rows, bool compute_def); + + def solve_mod(unsigned x, unsigned_vector const& divide_rows, bool compute_def); + + def solve_div(unsigned x, unsigned_vector const& divide_rows, bool compute_def); + bool is_int(unsigned x) const { return m_var2is_int[x]; } void retire_row(unsigned row_id); @@ -173,6 +192,17 @@ namespace opt { // add a divisibility constraint. The row should divide m. void add_divides(vector const& coeffs, rational const& c, rational const& m); + + // add sub-expression for modulus + // v = add_mod(coeffs, m) means + // v = coeffs mod m + unsigned add_mod(vector const& coeffs, rational const& c, rational const& m); + + // add sub-expression for div + // v = add_div(coeffs, m) means + // v = coeffs div m + unsigned add_div(vector const& coeffs, rational const& c, rational const& m); + // Set the objective function (linear). void set_objective(vector const& coeffs, rational const& c); diff --git a/src/qe/mbp/mbp_arith.cpp b/src/qe/mbp/mbp_arith.cpp index 377b62ac0..713a63c14 100644 --- a/src/qe/mbp/mbp_arith.cpp +++ b/src/qe/mbp/mbp_arith.cpp @@ -23,7 +23,6 @@ Revision History: #include "ast/ast_util.h" #include "ast/arith_decl_plugin.h" #include "ast/ast_pp.h" -#include "ast/rewriter/th_rewriter.h" #include "ast/expr_functors.h" #include "ast/rewriter/expr_safe_replace.h" #include "math/simplex/model_based_opt.h" @@ -32,13 +31,13 @@ Revision History: #include "model/model_v2_pp.h" namespace mbp { - + struct arith_project_plugin::imp { - ast_manager& m; + ast_manager& m; arith_util a; - bool m_check_purified { true }; // check that variables are properly pure - bool m_apply_projection { false }; + bool m_check_purified = true; // check that variables are properly pure + bool m_apply_projection = false; imp(ast_manager& m) : @@ -48,10 +47,10 @@ namespace mbp { void insert_mul(expr* x, rational const& v, obj_map& ts) { rational w; - if (ts.find(x, w)) - ts.insert(x, w + v); - else - ts.insert(x, v); + if (ts.find(x, w)) + ts.insert(x, w + v); + else + ts.insert(x, v); } @@ -65,15 +64,15 @@ namespace mbp { rational c(0), mul(1); expr_ref t(m); opt::ineq_type ty = opt::t_le; - expr* e1, *e2; - DEBUG_CODE(expr_ref val(m); - eval(lit, val); - CTRACE("qe", !m.is_true(val), tout << mk_pp(lit, m) << " := " << val << "\n";); - SASSERT(m.limit().is_canceled() || !m.is_false(val));); + expr* e1, * e2; + DEBUG_CODE(expr_ref val(m); + eval(lit, val); + CTRACE("qe", !m.is_true(val), tout << mk_pp(lit, m) << " := " << val << "\n";); + SASSERT(m.limit().is_canceled() || !m.is_false(val));); - if (!m.inc()) + if (!m.inc()) return false; - + TRACE("opt", tout << mk_pp(lit, m) << " " << a.is_lt(lit) << " " << a.is_gt(lit) << "\n";); bool is_not = m.is_not(lit, lit); if (is_not) { @@ -86,37 +85,35 @@ namespace mbp { ty = is_not ? opt::t_lt : opt::t_le; } else if ((a.is_lt(lit, e1, e2) || a.is_gt(lit, e2, e1))) { - linearize(mbo, eval, mul, e1, c, fmls, ts, tids); + linearize(mbo, eval, mul, e1, c, fmls, ts, tids); linearize(mbo, eval, -mul, e2, c, fmls, ts, tids); - ty = is_not ? opt::t_le: opt::t_lt; + ty = is_not ? opt::t_le : opt::t_lt; } else if (m.is_eq(lit, e1, e2) && !is_not && is_arith(e1)) { - linearize(mbo, eval, mul, e1, c, fmls, ts, tids); + linearize(mbo, eval, mul, e1, c, fmls, ts, tids); linearize(mbo, eval, -mul, e2, c, fmls, ts, tids); ty = opt::t_eq; - } + } else if (m.is_eq(lit, e1, e2) && is_not && is_arith(e1)) { - + rational r1, r2; - expr_ref val1 = eval(e1); + expr_ref val1 = eval(e1); expr_ref val2 = eval(e2); - //TRACE("qe", tout << mk_pp(e1, m) << " " << val1 << "\n";); - //TRACE("qe", tout << mk_pp(e2, m) << " " << val2 << "\n";); if (!a.is_numeral(val1, r1)) return false; if (!a.is_numeral(val2, r2)) return false; SASSERT(r1 != r2); if (r1 < r2) { std::swap(e1, e2); - } + } ty = opt::t_lt; - linearize(mbo, eval, mul, e1, c, fmls, ts, tids); - linearize(mbo, eval, -mul, e2, c, fmls, ts, tids); - } + linearize(mbo, eval, mul, e1, c, fmls, ts, tids); + linearize(mbo, eval, -mul, e2, c, fmls, ts, tids); + } else if (m.is_distinct(lit) && !is_not && is_arith(to_app(lit)->get_arg(0))) { expr_ref val(m); rational r; app* alit = to_app(lit); - vector > nums; + vector > nums; for (expr* arg : *alit) { val = eval(arg); TRACE("qe", tout << mk_pp(arg, m) << " " << val << "\n";); @@ -125,8 +122,8 @@ namespace mbp { } std::sort(nums.begin(), nums.end(), compare_second()); for (unsigned i = 0; i + 1 < nums.size(); ++i) { - SASSERT(nums[i].second < nums[i+1].second); - expr_ref fml(a.mk_lt(nums[i].first, nums[i+1].first), m); + SASSERT(nums[i].second < nums[i + 1].second); + expr_ref fml(a.mk_lt(nums[i].first, nums[i + 1].first), m); if (!linearize(mbo, eval, fml, fmls, tids)) { return false; } @@ -139,14 +136,14 @@ namespace mbp { map values; bool found_eq = false; for (unsigned i = 0; !found_eq && i < to_app(lit)->get_num_args(); ++i) { - expr* arg1 = to_app(lit)->get_arg(i), *arg2 = nullptr; + expr* arg1 = to_app(lit)->get_arg(i), * arg2 = nullptr; rational r; expr_ref val = eval(arg1); TRACE("qe", tout << mk_pp(arg1, m) << " " << val << "\n";); if (!a.is_numeral(val, r)) return false; if (values.find(r, arg2)) { ty = opt::t_eq; - linearize(mbo, eval, mul, arg1, c, fmls, ts, tids); + linearize(mbo, eval, mul, arg1, c, fmls, ts, tids); linearize(mbo, eval, -mul, arg2, c, fmls, ts, tids); found_eq = true; } @@ -169,27 +166,37 @@ namespace mbp { // // convert linear arithmetic term into an inequality for mbo. // - void linearize(opt::model_based_opt& mbo, model_evaluator& eval, rational const& mul, expr* t, rational& c, - expr_ref_vector& fmls, obj_map& ts, obj_map& tids) { - expr* t1, *t2, *t3; + void linearize(opt::model_based_opt& mbo, model_evaluator& eval, rational const& mul, expr* t, rational& c, + expr_ref_vector& fmls, obj_map& ts, obj_map& tids) { + expr* t1, * t2, * t3; rational mul1; expr_ref val(m); - if (a.is_mul(t, t1, t2) && is_numeral(t1, mul1)) - linearize(mbo, eval, mul* mul1, t2, c, fmls, ts, tids); - else if (a.is_mul(t, t1, t2) && is_numeral(t2, mul1)) - linearize(mbo, eval, mul* mul1, t1, c, fmls, ts, tids); + + auto add_def = [&](expr* t1, rational const& m, vars& coeffs) { + obj_map ts0; + rational mul0(1), c0(0); + linearize(mbo, eval, mul0, t1, c0, fmls, ts0, tids); + extract_coefficients(mbo, eval, ts0, tids, coeffs); + insert_mul(t, mul, ts); + return c0; + }; + + if (a.is_mul(t, t1, t2) && is_numeral(t1, mul1)) + linearize(mbo, eval, mul * mul1, t2, c, fmls, ts, tids); + else if (a.is_mul(t, t1, t2) && is_numeral(t2, mul1)) + linearize(mbo, eval, mul * mul1, t1, c, fmls, ts, tids); else if (a.is_uminus(t, t1)) linearize(mbo, eval, -mul, t1, c, fmls, ts, tids); - else if (a.is_numeral(t, mul1)) - c += mul * mul1; + else if (a.is_numeral(t, mul1)) + c += mul * mul1; else if (a.is_add(t)) { - for (expr* arg : *to_app(t)) - linearize(mbo, eval, mul, arg, c, fmls, ts, tids); + for (expr* arg : *to_app(t)) + linearize(mbo, eval, mul, arg, c, fmls, ts, tids); } else if (a.is_sub(t, t1, t2)) { - linearize(mbo, eval, mul, t1, c, fmls, ts, tids); + linearize(mbo, eval, mul, t1, c, fmls, ts, tids); linearize(mbo, eval, -mul, t2, c, fmls, ts, tids); - } + } else if (m.is_ite(t, t1, t2, t3)) { val = eval(t1); @@ -210,21 +217,31 @@ namespace mbp { else if (a.is_mod(t, t1, t2) && is_numeral(t2, mul1) && !mul1.is_zero()) { rational r; val = eval(t); - if (!a.is_numeral(val, r)) { + if (!a.is_numeral(val, r)) throw default_exception("mbp evaluation didn't produce an integer"); - } - c += mul*r; - // t1 mod mul1 == r - rational c0(-r), mul0(1); - obj_map ts0; - linearize(mbo, eval, mul0, t1, c0, fmls, ts0, tids); + c += mul * r; + + rational c0(-r), mul0(1); + obj_map ts0; + linearize(mbo, eval, mul0, t1, c0, fmls, ts0, tids); + vars coeffs; + extract_coefficients(mbo, eval, ts0, tids, coeffs); + mbo.add_divides(coeffs, c0, mul1); + } + else if (false && a.is_mod(t, t1, t2) && is_numeral(t2, mul1) && mul1 > 0) { + // v = t1 mod mul1 vars coeffs; - extract_coefficients(mbo, eval, ts0, tids, coeffs); - mbo.add_divides(coeffs, c0, mul1); + rational c0 = add_def(t1, mul1, coeffs); + tids.insert(t, mbo.add_mod(coeffs, c0, mul1)); } - else { + else if (false && a.is_idiv(t, t1, t2) && is_numeral(t2, mul1) && mul1 > 0) { + // v = t1 div mul1 + vars coeffs; + rational c0 = add_def(t1, mul1, coeffs); + tids.insert(t, mbo.add_div(coeffs, c0, mul1)); + } + else insert_mul(t, mul, ts); - } } bool is_numeral(expr* t, rational& r) { @@ -232,8 +249,8 @@ namespace mbp { } struct compare_second { - bool operator()(std::pair const& a, - std::pair const& b) const { + bool operator()(std::pair const& a, + std::pair const& b) const { return a.second < b.second; } }; @@ -243,14 +260,14 @@ namespace mbp { } rational n_sign(rational const& b) { - return rational(b.is_pos()?-1:1); + return rational(b.is_pos() ? -1 : 1); } bool operator()(model& model, app* v, app_ref_vector& vars, expr_ref_vector& lits) { app_ref_vector vs(m); vs.push_back(v); vector defs; - return project(model, vs, lits, defs, false) && vs.empty(); + return project(model, vs, lits, defs, false) && vs.empty(); } typedef opt::model_based_opt::var var; @@ -267,10 +284,10 @@ namespace mbp { bool project(model& model, app_ref_vector& vars, expr_ref_vector& fmls, vector& result, bool compute_def) { bool has_arith = false; - for (expr* v : vars) - has_arith |= is_arith(v); - if (!has_arith) - return true; + for (expr* v : vars) + has_arith |= is_arith(v); + if (!has_arith) + return true; model_evaluator eval(model); TRACE("qe", tout << model;); eval.set_model_completion(true); @@ -281,7 +298,7 @@ namespace mbp { expr_ref_vector pinned(m); unsigned j = 0; TRACE("qe", tout << "vars: " << vars << "\n"; - for (expr* f : fmls) tout << mk_pp(f, m) << " := " << model(f) << "\n";); + for (expr* f : fmls) tout << mk_pp(f, m) << " := " << model(f) << "\n";); for (unsigned i = 0; i < fmls.size(); ++i) { expr* fml = fmls.get(i); if (!linearize(mbo, eval, fml, fmls, tids)) { @@ -294,8 +311,8 @@ namespace mbp { } fmls.shrink(j); TRACE("qe", tout << "formulas\n" << fmls << "\n"; - for (auto const& [e, id] : tids) - tout << mk_pp(e, m) << " -> " << id << "\n";); + for (auto const& [e, id] : tids) + tout << mk_pp(e, m) << " -> " << id << "\n";); // fmls holds residue, // mbo holds linear inequalities that are in scope @@ -306,7 +323,7 @@ namespace mbp { // return those to fmls. expr_mark var_mark, fmls_mark; - for (app * v : vars) { + for (app* v : vars) { var_mark.mark(v); if (is_arith(v) && !tids.contains(v)) { rational r; @@ -321,60 +338,70 @@ namespace mbp { } // bail on variables in non-linear sub-terms + auto is_pure = [&](expr* e) { + expr* x, * y; + rational r; + if (a.is_mod(e, x, y) && a.is_numeral(y)) + return true; + if (a.is_idiv(e, x, y) && a.is_numeral(y, r) && r > 0) + return true; + return false; + }; + for (auto& kv : tids) { expr* e = kv.m_key; - if (is_arith(e) && !var_mark.is_marked(e)) - mark_rec(fmls_mark, e); + if (is_arith(e) && !is_pure(e) && !var_mark.is_marked(e)) + mark_rec(fmls_mark, e); } if (m_check_purified) { - for (expr* fml : fmls) - mark_rec(fmls_mark, fml); + for (expr* fml : fmls) + mark_rec(fmls_mark, fml); for (auto& kv : tids) { expr* e = kv.m_key; - if (!var_mark.is_marked(e)) - mark_rec(fmls_mark, e); + if (!var_mark.is_marked(e) && !is_pure(e)) + mark_rec(fmls_mark, e); } } ptr_vector index2expr; - for (auto& kv : tids) - index2expr.setx(kv.m_value, kv.m_key, nullptr); + for (auto& kv : tids) + index2expr.setx(kv.m_value, kv.m_key, nullptr); j = 0; unsigned_vector real_vars; for (app* v : vars) { - if (is_arith(v) && !fmls_mark.is_marked(v)) - real_vars.push_back(tids.find(v)); - else - vars[j++] = v; + if (is_arith(v) && !fmls_mark.is_marked(v)) + real_vars.push_back(tids.find(v)); + else + vars[j++] = v; } vars.shrink(j); - - TRACE("qe", tout << "remaining vars: " << vars << "\n"; - for (unsigned v : real_vars) tout << "v" << v << " " << mk_pp(index2expr[v], m) << "\n"; - mbo.display(tout);); + + TRACE("qe", tout << "remaining vars: " << vars << "\n"; + for (unsigned v : real_vars) tout << "v" << v << " " << mk_pp(index2expr[v], m) << "\n"; + mbo.display(tout);); vector defs = mbo.project(real_vars.size(), real_vars.data(), compute_def); vector rows; mbo.get_live_rows(rows); rows2fmls(rows, index2expr, fmls); TRACE("qe", mbo.display(tout << "mbo result\n"); - for (auto const& d : defs) tout << "def: " << d << "\n"; - tout << fmls << "\n";); - - if (compute_def) - optdefs2mbpdef(defs, index2expr, real_vars, result); + for (auto const& d : defs) tout << "def: " << d << "\n"; + tout << fmls << "\n";); + + if (compute_def) + optdefs2mbpdef(defs, index2expr, real_vars, result); if (m_apply_projection && !apply_projection(eval, result, fmls)) return false; TRACE("qe", for (auto const& [v, t] : result) tout << v << " := " << t << "\n"; - for (auto* f : fmls) - tout << mk_pp(f, m) << " := " << eval(f) << "\n"; - tout << "fmls:" << fmls << "\n";); + for (auto* f : fmls) + tout << mk_pp(f, m) << " := " << eval(f) << "\n"; + tout << "fmls:" << fmls << "\n";); return true; - } + } void optdefs2mbpdef(vector const& defs, ptr_vector const& index2expr, unsigned_vector const& real_vars, vector& result) { SASSERT(defs.size() == real_vars.size()); @@ -384,31 +411,92 @@ namespace mbp { bool is_int = a.is_int(x); expr_ref_vector ts(m); expr_ref t(m); - for (var const& v : d.m_vars) - ts.push_back(var2expr(index2expr, v)); + for (var const& v : d.m_vars) + ts.push_back(var2expr(index2expr, v)); if (!d.m_coeff.is_zero()) ts.push_back(a.mk_numeral(d.m_coeff, is_int)); if (ts.empty()) ts.push_back(a.mk_numeral(rational(0), is_int)); t = mk_add(ts); - if (!d.m_div.is_one() && is_int) - t = a.mk_idiv(t, a.mk_numeral(d.m_div, is_int)); - else if (!d.m_div.is_one() && !is_int) - t = a.mk_div(t, a.mk_numeral(d.m_div, is_int)); + if (!d.m_div.is_one() && is_int) + t = a.mk_idiv(t, a.mk_numeral(d.m_div, is_int)); + else if (!d.m_div.is_one() && !is_int) + t = a.mk_div(t, a.mk_numeral(d.m_div, is_int)); result.push_back(def(expr_ref(x, m), t)); } } - void rows2fmls(vector const& rows, ptr_vector const& index2expr, expr_ref_vector& fmls) { - for (row const& r : rows) { - expr_ref_vector ts(m); - expr_ref t(m), s(m), val(m); - if (r.m_vars.empty()) { + expr_ref id2expr(u_map const& def_vars, ptr_vector const& index2expr, unsigned id) { + row r; + if (def_vars.find(id, r)) + return row2expr(def_vars, index2expr, r); + return expr_ref(index2expr[id], m); + } + + expr_ref row2expr(u_map const& def_vars, ptr_vector const& index2expr, row const& r) { + expr_ref_vector ts(m); + expr_ref t(m); + rational n; + for (var const& v : r.m_vars) { + t = id2expr(def_vars, index2expr, v.m_id); + if (a.is_numeral(t, n) && n == 0) continue; + else if (a.is_numeral(t, n)) + t = a.mk_numeral(v.m_coeff * n, a.is_int(t)); + else if (!v.m_coeff.is_one()) + t = a.mk_mul(a.mk_numeral(v.m_coeff, a.is_int(t)), t); + ts.push_back(t); + } + switch (r.m_type) { + case opt::t_mod: + if (ts.empty()) { + t = a.mk_int(mod(r.m_coeff, r.m_mod)); + return t; } - if (r.m_vars.size() == 1 && r.m_vars[0].m_coeff.is_neg() && r.m_type != opt::t_mod) { + ts.push_back(a.mk_int(r.m_coeff)); + t = mk_add(ts); + t = a.mk_mod(t, a.mk_int(r.m_mod)); + return t; + case opt::t_div: + if (ts.empty()) { + t = a.mk_int(div(r.m_coeff, r.m_mod)); + return t; + } + ts.push_back(a.mk_int(r.m_coeff)); + t = mk_add(ts); + t = a.mk_idiv(t, a.mk_int(r.m_mod)); + return t; + case opt::t_divides: + ts.push_back(a.mk_int(r.m_coeff)); + return mk_add(ts); + default: + return mk_add(ts); + } + } + + void rows2fmls(vector const& rows, ptr_vector const& index2expr, expr_ref_vector& fmls) { + + u_map def_vars; + for (row const& r : rows) { + if (r.m_type == opt::t_mod) + def_vars.insert(r.m_id, r); + else if (r.m_type == opt::t_div) + def_vars.insert(r.m_id, r); + } + + for (row const& r : rows) { + expr_ref t(m), s(m), val(m); + + if (r.m_vars.empty()) + continue; + if (r.m_type == opt::t_mod) + continue; + if (r.m_type == opt::t_div) + continue; + + if (r.m_vars.size() == 1 && r.m_vars[0].m_coeff.is_neg() && r.m_type != opt::t_divides) { var const& v = r.m_vars[0]; - t = index2expr[v.m_id]; + t = id2expr(def_vars, index2expr, v.m_id); if (!v.m_coeff.is_minus_one()) { t = a.mk_mul(a.mk_numeral(-v.m_coeff, a.is_int(t)), t); } @@ -422,24 +510,14 @@ namespace mbp { fmls.push_back(t); continue; } - for (var const& v : r.m_vars) { - t = index2expr[v.m_id]; - if (!v.m_coeff.is_one()) { - t = a.mk_mul(a.mk_numeral(v.m_coeff, a.is_int(t)), t); - } - ts.push_back(t); - } - t = mk_add(ts); + t = row2expr(def_vars, index2expr, r); s = a.mk_numeral(-r.m_coeff, r.m_coeff.is_int() && a.is_int(t)); switch (r.m_type) { case opt::t_lt: t = a.mk_lt(t, s); break; case opt::t_le: t = a.mk_le(t, s); break; case opt::t_eq: t = a.mk_eq(t, s); break; - case opt::t_mod: - if (!r.m_coeff.is_zero()) { - t = a.mk_sub(t, s); - } - t = a.mk_eq(a.mk_mod(t, a.mk_numeral(r.m_mod, true)), a.mk_int(0)); + case opt::t_divides: + t = a.mk_eq(a.mk_mod(t, a.mk_int(r.m_mod)), a.mk_int(0)); break; } fmls.push_back(t); @@ -468,11 +546,11 @@ namespace mbp { SASSERT(validate_model(eval, fmls0)); // extract linear constraints - - for (expr * fml : fmls) { + + for (expr* fml : fmls) { linearize(mbo, eval, fml, fmls, tids); } - + // find optimal value value = mbo.maximize(); @@ -543,7 +621,7 @@ namespace mbp { } CTRACE("qe", kv.m_value.is_zero(), tout << mk_pp(v, m) << " has coefficeint 0\n";); if (!kv.m_value.is_zero()) { - coeffs.push_back(var(id, kv.m_value)); + coeffs.push_back(var(id, kv.m_value)); } } } @@ -571,7 +649,7 @@ namespace mbp { }; - arith_project_plugin::arith_project_plugin(ast_manager& m):project_plugin(m) { + arith_project_plugin::arith_project_plugin(ast_manager& m) :project_plugin(m) { m_imp = alloc(imp, m); } diff --git a/src/qe/qsat.cpp b/src/qe/qsat.cpp index 767b96e5d..6023273c0 100644 --- a/src/qe/qsat.cpp +++ b/src/qe/qsat.cpp @@ -241,9 +241,8 @@ namespace qe { while (sz0 != todo.size()) { app* a = to_app(todo.back()); todo.pop_back(); - if (mark.is_marked(a)) { + if (mark.is_marked(a)) continue; - } mark.mark(a); if (m_lit2pred.find(a, p)) { @@ -284,9 +283,8 @@ namespace qe { m_elevel.insert(r, l); eq = m.mk_eq(r, a); defs.push_back(eq); - if (!is_predicate(a, l.max())) { + if (!is_predicate(a, l.max())) insert(r, l); - } level.merge(l); } } @@ -637,57 +635,55 @@ namespace qe { check_cancel(); expr_ref_vector asms(m_asms); m_pred_abs.get_assumptions(m_model.get(), asms); - if (m_model.get()) { + if (m_model.get()) validate_assumptions(*m_model.get(), asms); - } TRACE("qe", tout << asms << "\n";); solver& s = get_kernel(m_level).s(); lbool res = s.check_sat(asms); switch (res) { case l_true: s.get_model(m_model); + CTRACE("qe", !m_model, tout << "no model\n"); if (!m_model) return l_undef; SASSERT(validate_defs("check_sat")); SASSERT(!m_model.get() || validate_assumptions(*m_model.get(), asms)); SASSERT(validate_model(asms)); TRACE("qe", s.display(tout); display(tout << "\n", *m_model.get()); display(tout, asms); ); - if (m_level == 0) { + if (m_level == 0) m_model_save = m_model; - } push(); - if (m_level == 1 && m_mode == qsat_maximize) { + if (m_level == 1 && m_mode == qsat_maximize) maximize_model(); - } break; case l_false: switch (m_level) { case 0: return l_false; case 1: - if (m_mode == qsat_sat) { + if (m_mode == qsat_sat) return l_true; - } if (m_model.get()) { SASSERT(validate_assumptions(*m_model.get(), asms)); - if (!project_qe(asms)) return l_undef; + if (!project_qe(asms)) + return l_undef; } - else { + else pop(1); - } break; default: if (m_model.get()) { - if (!project(asms)) return l_undef; + if (!project(asms)) + return l_undef; } - else { + else pop(1); - } break; } break; case l_undef: + TRACE("qe", tout << "check-sat is undef\n"); return res; } } @@ -833,11 +829,10 @@ namespace qe { } } - bool get_core(expr_ref_vector& core, unsigned level) { + void get_core(expr_ref_vector& core, unsigned level) { SASSERT(validate_defs("get_core")); get_kernel(level).get_core(core); m_pred_abs.pred2lit(core); - return true; } bool minimize_core(expr_ref_vector& core, unsigned level) { @@ -905,9 +900,7 @@ namespace qe { SASSERT(m_level == 1); expr_ref fml(m); model& mdl = *m_model.get(); - if (!get_core(core, m_level)) { - return false; - } + get_core(core, m_level); SASSERT(validate_core(mdl, core)); get_vars(m_level); SASSERT(validate_assumptions(mdl, core)); @@ -927,7 +920,7 @@ namespace qe { } bool project(expr_ref_vector& core) { - if (!get_core(core, m_level)) return false; + get_core(core, m_level); TRACE("qe", display(tout); display(tout << "core\n", core);); SASSERT(m_level >= 2); expr_ref fml(m); @@ -950,14 +943,17 @@ namespace qe { if (level.max() == UINT_MAX) { num_scopes = 2*(m_level/2); } + else if (level.max() + 2 > m_level) { + // fishy - this can happen. + TRACE("qe", tout << "max-level: " << level.max() << " level: " << m_level << "\n"); + return false; + } else { - if (level.max() + 2 > m_level) return false; SASSERT(level.max() + 2 <= m_level); num_scopes = m_level - level.max(); SASSERT(num_scopes >= 2); - if ((num_scopes % 2) != 0) { + if ((num_scopes % 2) != 0) --num_scopes; - } } pop(num_scopes); diff --git a/src/sat/smt/xor_solver.cpp b/src/sat/smt/xor_solver.cpp new file mode 100644 index 000000000..27354eac2 --- /dev/null +++ b/src/sat/smt/xor_solver.cpp @@ -0,0 +1,79 @@ +/*++ +Copyright (c) 2014 Microsoft Corporation + +Module Name: + + xor_solver.h + +Abstract: + + XOR solver. + Interface outline. + +--*/ + + +#include "sat/smt/xor_solver.h" +#include "sat/sat_simplifier_params.hpp" +#include "sat/sat_xor_finder.h" + +namespace xr { + + solver::solver(euf::solver& ctx): + th_solver(ctx.get_manager(), symbol("xor-solver"), ctx.get_manager().get_family_id("xor-solver")) + {} + + euf::th_solver* solver::clone(euf::solver& ctx) { + // and relevant copy internal state + return alloc(solver, ctx); + } + + void solver::asserted(sat::literal l) { + + } + + bool solver::unit_propagate() { + return false; + } + + void solver::get_antecedents(sat::literal l, sat::ext_justification_idx idx, + sat::literal_vector & r, bool probing) { + + } + + sat::check_result solver::check() { + return sat::check_result::CR_DONE; + } + + void solver::push() { + } + + void solver::pop(unsigned n) { + } + + // inprocessing + // pre_simplify: decompile all xor constraints to allow other in-processing. + // simplify: recompile clauses to xor constraints + // literals that get added to xor constraints are tagged with the theory. + void solver::pre_simplify() { + + } + + void solver::simplify() { + + } + + std::ostream& solver::display(std::ostream& out) const { + return out; + } + + std::ostream& solver::display_justification(std::ostream& out, sat::ext_justification_idx idx) const { + return out; + } + + std::ostream& solver::display_constraint(std::ostream& out, sat::ext_constraint_idx idx) const { + return out; + } + +} + diff --git a/src/sat/smt/xor_solver.h b/src/sat/smt/xor_solver.h new file mode 100644 index 000000000..3da30c580 --- /dev/null +++ b/src/sat/smt/xor_solver.h @@ -0,0 +1,48 @@ +/*++ +Copyright (c) 2014 Microsoft Corporation + +Module Name: + + xor_solver.h + +Abstract: + + XOR solver. + Interface outline. + +--*/ + +#pragma once + +#include "sat/smt/euf_solver.h" + +namespace xr { + class solver : public euf::th_solver { + public: + solver(euf::solver& ctx); + + th_solver* clone(euf::solver& ctx) override; + + sat::literal internalize(expr* e, bool sign, bool root, bool redundant) override { UNREACHABLE(); return sat::null_literal; } + + void internalize(expr* e, bool redundant) override { UNREACHABLE(); } + + + void asserted(sat::literal l) override; + bool unit_propagate() override; + void get_antecedents(sat::literal l, sat::ext_justification_idx idx, sat::literal_vector & r, bool probing) override; + + void pre_simplify() override; + void simplify() override; + + sat::check_result check() override; + void push() override; + void pop(unsigned n) override; + + std::ostream& display(std::ostream& out) const override; + std::ostream& display_justification(std::ostream& out, sat::ext_justification_idx idx) const override; + std::ostream& display_constraint(std::ostream& out, sat::ext_constraint_idx idx) const override; + + }; + +} diff --git a/src/smt/qi_queue.cpp b/src/smt/qi_queue.cpp index de52cdb01..f8fbe739e 100644 --- a/src/smt/qi_queue.cpp +++ b/src/smt/qi_queue.cpp @@ -23,6 +23,7 @@ Revision History: #include "ast/rewriter/var_subst.h" #include "smt/smt_context.h" #include "smt/qi_queue.h" +#include namespace smt { @@ -228,15 +229,12 @@ namespace smt { TRACE("qi_queue", tout << "new instance:\n" << mk_pp(instance, m) << "\n";); - TRACE("qi_queue_instance", tout << "new instance:\n" << mk_pp(instance, m) << "\n";); expr_ref s_instance(m); proof_ref pr(m); m_context.get_rewriter()(instance, s_instance, pr); TRACE("qi_queue_bug", tout << "new instance after simplification:\n" << s_instance << "\n";); if (m.is_true(s_instance)) { - TRACE("checker", tout << "reduced to true, before:\n" << mk_ll_pp(instance, m);); - STRACE("instance", tout << "Instance reduced to true\n";); stat -> inc_num_instances_simplify_true(); if (m.has_trace_stream()) { @@ -250,9 +248,11 @@ namespace smt { std::cout << "instantiate\n"; enode_vector _bindings(num_bindings, bindings); for (auto * b : _bindings) - std::cout << enode_pp(b, m_context) << " "; + std::cout << mk_pp(b->get_expr(), m) << " "; std::cout << "\n"; std::cout << mk_pp(q, m) << "\n"; + std::cout << "instance\n"; + std::cout << instance << "\n"; #endif TRACE("qi_queue", tout << "simplified instance:\n" << s_instance << "\n";); diff --git a/src/smt/smt_model_checker.cpp b/src/smt/smt_model_checker.cpp index 039416f73..0ffe8108d 100644 --- a/src/smt/smt_model_checker.cpp +++ b/src/smt/smt_model_checker.cpp @@ -218,7 +218,7 @@ namespace smt { TRACE("model_checker", tout << "Got some value " << sk_value << "\n";); if (use_inv) { - unsigned sk_term_gen = 0; + unsigned sk_term_gen = 0; expr * sk_term = m_model_finder.get_inv(q, i, sk_value, sk_term_gen); if (sk_term != nullptr) { TRACE("model_checker", tout << "Found inverse " << mk_pp(sk_term, m) << "\n";); @@ -243,6 +243,8 @@ namespace smt { func_decl * f = nullptr; if (autil.is_as_array(sk_value, f) && cex->get_func_interp(f) && cex->get_func_interp(f)->get_interp()) { expr_ref body(cex->get_func_interp(f)->get_interp(), m); + if (contains_model_value(body)) + return false; ptr_vector sorts(f->get_arity(), f->get_domain()); svector names; for (unsigned i = 0; i < f->get_arity(); ++i) { diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp index b8001174e..7a46e1a1e 100644 --- a/src/smt/theory_lra.cpp +++ b/src/smt/theory_lra.cpp @@ -3507,16 +3507,15 @@ public: case lp::lp_status::OPTIMAL: { init_variable_values(); TRACE("arith", display(tout << st << " v" << v << " vi: " << vi << "\n");); - inf_rational val = get_value(v); - // inf_rational val(term_max.x, term_max.y); + auto val = value(v); blocker = mk_gt(v); - return inf_eps(rational::zero(), val); + return val; } case lp::lp_status::FEASIBLE: { - inf_rational val = get_value(v); + auto val = value(v); TRACE("arith", display(tout << st << " v" << v << " vi: " << vi << "\n");); blocker = mk_gt(v); - return inf_eps(rational::zero(), val); + return val; } default: SASSERT(st == lp::lp_status::UNBOUNDED); @@ -3533,23 +3532,19 @@ public: rational r = val.x; expr_ref e(m); if (a.is_int(obj->get_sort())) { - if (r.is_int()) { + if (r.is_int()) r += rational::one(); - } - else { + else r = ceil(r); - } e = a.mk_numeral(r, obj->get_sort()); e = a.mk_ge(obj, e); } else { e = a.mk_numeral(r, obj->get_sort()); - if (val.y.is_neg()) { + if (val.y.is_neg()) e = a.mk_ge(obj, e); - } - else { + else e = a.mk_gt(obj, e); - } } TRACE("opt", tout << "v" << v << " " << val << " " << r << " " << e << "\n";); return e;