diff --git a/src/math/simplex/model_based_opt.cpp b/src/math/simplex/model_based_opt.cpp index 0223e8bff..be2b4fe2b 100644 --- a/src/math/simplex/model_based_opt.cpp +++ b/src/math/simplex/model_based_opt.cpp @@ -1,1300 +1,1564 @@ -/*++ -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())); + 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) + 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 << 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 = 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_case2 = distance.is_nonpos() || abs_src_c.is_one() || abs_dst_c.is_one(); + + if (use_case1 && false) { + // + // 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) + // + + if (src_c < 0) + std::swap(row_src, row_dst); + + 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 = add_div(src_coeffs, -src_coeff, abs_src_c); + dst_coeffs.push_back(var(v, -abs_src_c)); + add_constraint(dst_coeffs, dst_coeff, t_le); + retire_row(row_dst); + return; + } + + if (use_case2) { + 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::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; + 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); + } + + void model_based_opt::add_divides(vector const& coeffs, rational const& c, rational const& m) { + add_constraint(coeffs, c, m, t_divides); + } + + unsigned model_based_opt::add_mod(vector const& coeffs, rational const& c, rational const& m) { + add_constraint(coeffs, c, m, t_mod); + unsigned row_id = m_rows.size() - 1; + auto& r = m_rows[row_id]; + unsigned v = add_var(mod(r.m_value, m), true); + r.m_id = v; + m_var2row_ids[v].push_back(row_id); + return v; + } + + unsigned model_based_opt::add_div(vector const& coeffs, rational const& c, rational const& m) { + add_constraint(coeffs, c, m, t_div); + unsigned row_id = m_rows.size() - 1; + auto& r = m_rows[row_id]; + unsigned v = add_var(div(r.m_value, m), true); + r.m_id = v; + m_var2row_ids[v].push_back(row_id); + return v; + } + + 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& 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 m + // + // - remove v = a*x + b mod m + // + // case a = 1: + // - add w = b mod m + // - x |-> m*y + z, 0 <= z < m + // - if z.value + w.value < m: + // add z + w - v = 0 + // - if z.value + w.value >= m: + // add z + w - v - m = 0 + // + // case a != 1, gcd(a, m) = 1 + // - x |-> x*y + a^-1*z, 0 <= z < m + // - add w = b mod m + // if z.value + w.value < m + // add z + w - v = 0 + // if z.value + w.value >= m + // add z + w - v - m = 0 + // + // case a != 1, gcd(a,m) = g != 1 + // - x |-> x*y + a^-1*z, 0 <= z < m + // a*x + b mod m = v is now + // g*z + b mod m = v + // - add w = b mod m + // - 0 <= g*z.value + w.value < m*(g+1) + // - add g*z + w - v - k*m = 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; + SASSERT(!mod_rows.empty()); + unsigned row_index = mod_rows.back(); + rational a = get_coefficient(row_index, x); + rational m = m_rows[row_index].m_mod; + unsigned v = m_rows[row_index].m_id; + rational x_value = m_var2value[x]; + // disable row + m_rows[row_index].m_alive = false; + replace_var(row_index, x, rational::zero()); + + + // compute a_inv + rational a_inv, m_inv; + rational g = gcd(a, m, a_inv, m_inv); + if (a_inv.is_neg()) + a_inv = mod(a_inv, m); + SASSERT(mod(a_inv * a, m) == g); + + // solve for x_value = m*y_value + a^-1*z_value, 0 <= z_value < m. + rational z_value = mod(x_value, m); + rational y_value = div(x_value, m) - div(a_inv*z_value, m); + SASSERT(x_value == m*y_value + a_inv*z_value); + SASSERT(0 <= z_value && z_value < m); + + // add new variables + unsigned y = add_var(y_value, true); + unsigned z = add_var(z_value, true); + // TODO: we could recycle x by either y or z. + + // replace x by m*y + a^-1*z in other rows. + unsigned_vector const& row_ids = m_var2row_ids[x]; + uint_set visited; + visited.insert(row_index); + for (unsigned row_id : row_ids) { + if (visited.contains(row_id)) + continue; + replace_var(row_id, x, m, y, a_inv, z); + visited.insert(row_id); + normalize(row_id); + } + + // add bounds for z + add_lower_bound(z, rational::zero()); + add_upper_bound(z, m - 1); + + + // add w = b mod m + vector coeffs = m_rows[row_index].m_vars; + rational coeff = m_rows[row_index].m_coeff; + + unsigned w = add_mod(coeffs, coeff, m); + + rational w_value = m_var2value[w]; + + // add g*z + w - v - k*m = 0, for k = (g*z_value + w_value) div m + rational km = div(g*z_value + w_value, m)*m; + vector mod_coeffs; + mod_coeffs.push_back(var(z, g)); + mod_coeffs.push_back(var(w, rational::one())); + mod_coeffs.push_back(var(v, rational::minus_one())); + add_constraint(mod_coeffs, km, t_eq); + add_lower_bound(v, rational::zero()); + add_upper_bound(v, m - 1); + + // allow to recycle row. + m_retired_rows.push_back(row_index); + + project(v, false); + def y_def = project(y, compute_def); + def z_def = project(z, compute_def); + + if (compute_def) { + result = (y_def * m) + z_def; + m_var2value[x] = eval(result); + } + + return result; + } + + // + // Given v = a*x + b div m + // Replace x |-> m*y + a_inv*z + // - w = b div m + // - v = ((m*y + g*z) + b) div m + // = a*y + (a_inv*z + b) div m + // = a*y + b div m + (b mod m + g*z) div m + // = a*y + b div m + k + // where k := (b.value mod m + g*z.value) div m + // + model_based_opt::def model_based_opt::solve_div(unsigned x, unsigned_vector const& div_rows, bool compute_def) { + def result; + SASSERT(!div_rows.empty()); + unsigned row_index = div_rows.back(); + rational a = get_coefficient(row_index, x); + rational m = m_rows[row_index].m_mod; + unsigned v = m_rows[row_index].m_id; + rational x_value = m_var2value[x]; + // disable row + m_rows[row_index].m_alive = false; + replace_var(row_index, x, rational::zero()); + rational b_value = m_rows[row_index].m_value; + + // compute a_inv + rational a_inv, m_inv; + rational g = gcd(a, m, a_inv, m_inv); + if (a_inv.is_neg()) + a_inv = mod(a_inv, m); + SASSERT(mod(a_inv * a, m) == g); + + // solve for x_value = m*y_value + a^-1*z_value, 0 <= z_value < m. + rational z_value = mod(x_value, m); + rational y_value = div(x_value, m) - div(a_inv*z_value, m); + SASSERT(x_value == m*y_value + a_inv*z_value); + SASSERT(0 <= z_value && z_value < m); + + // add new variables + unsigned y = add_var(y_value, true); + unsigned z = add_var(z_value, true); + // TODO: we could recycle x by either y or z. + + // replace x by m*y + a^-1*z in other rows. + unsigned_vector const& row_ids = m_var2row_ids[x]; + uint_set visited; + visited.insert(row_index); + for (unsigned row_id : row_ids) { + if (visited.contains(row_id)) + continue; + replace_var(row_id, x, m, y, a_inv, z); + visited.insert(row_id); + normalize(row_id); + } + + // add bounds for z + add_lower_bound(z, rational::zero()); + add_upper_bound(z, m - 1); + + // add w = b div m + vector coeffs = m_rows[row_index].m_vars; + rational coeff = m_rows[row_index].m_coeff; + + unsigned w = add_div(coeffs, coeff, m); + rational k = div(g*z_value + mod(b_value, m), m); + vector div_coeffs; + div_coeffs.push_back(var(v, rational::minus_one())); + div_coeffs.push_back(var(y, a)); + div_coeffs.push_back(var(w, rational::one())); + add_constraint(div_coeffs, k, t_eq); + + // allow to recycle row. + m_retired_rows.push_back(row_index); + project(v, false); + def y_def = project(y, compute_def); + def z_def = project(z, compute_def); + + if (compute_def) { + result = (y_def * m) + z_def; + m_var2value[x] = eval(result); + } + + 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()); + r.m_vars.push_back(var(y, coeff*A)); + 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()); + m_var2row_ids[y].push_back(row_id); + 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";); + 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); + 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..9e1a67721 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 }; @@ -57,6 +59,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 +88,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; @@ -122,14 +125,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_upper_bound(unsigned x, rational const& hi); void add_constraint(vector const& coeffs, rational const& c, rational const& m, ineq_type r); 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 +145,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 +157,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 +184,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..7674a1577 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,30 @@ 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; + c += mul * r; // t1 mod mul1 == 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); + rational c0 = add_def(t1, mul1, coeffs); + mbo.add_divides(coeffs, c0 - r, mul1); } - else { + else if (false && a.is_mod(t, t1, t2) && is_numeral(t2, mul1) && !mul1.is_zero()) { + // v = t1 mod mul1 + vars coeffs; + rational c0 = add_def(t1, mul1, coeffs); + unsigned v = mbo.add_mod(coeffs, c0, mul1); + tids.insert(t, v); + } + 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); + unsigned v = mbo.add_div(coeffs, c0, mul1); + tids.insert(t, v); + } + else insert_mul(t, mul, ts); - } } bool is_numeral(expr* t, rational& r) { @@ -232,8 +248,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 +259,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 +283,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 +297,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 +310,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 +322,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 +337,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 +410,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 +509,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 +545,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 +620,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 +648,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); }