From 03385bf78d79a6eb884ab2d29a5ad685b89037f2 Mon Sep 17 00:00:00 2001 From: Nikolaj Bjorner Date: Fri, 12 Aug 2022 10:20:43 -0400 Subject: [PATCH] improve quantifier elimination for arithmetic This update changes the handling of mod and adds support for nested div terms. Simple use cases that are handled using small results are given below. ``` (declare-const x Int) (declare-const y Int) (declare-const z Int) (assert (exists ((x Int)) (and (<= y (* 10 x)) (<= (* 10 x) z)))) (apply qe2) (reset) (declare-const y Int) (assert (exists ((x Int)) (and (> x 0) (= (div x 41) y)))) (apply qe2) (reset) (declare-const y Int) (assert (exists ((x Int)) (= (mod x 41) y))) (apply qe2) (reset) ``` The main idea is to introduce definition rows for mod/div terms. Elimination of variables under mod/div is defined by rewriting the variable to multiples of the mod/divisior and remainder. The functionality is disabled in this push. --- src/math/simplex/model_based_opt.cpp | 2864 ++++++++++++++------------ src/math/simplex/model_based_opt.h | 38 +- src/qe/mbp/mbp_arith.cpp | 321 +-- 3 files changed, 1793 insertions(+), 1430 deletions(-) diff --git a/src/math/simplex/model_based_opt.cpp b/src/math/simplex/model_based_opt.cpp index 0223e8bff83..be2b4fe2b02 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 cb8c1854233..9e1a67721b9 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 377b62ac050..7674a157703 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 { - insert_mul(t, mul, ts); + 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)); } } + 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; + } + 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_vector ts(m); expr_ref t(m), s(m), val(m); - if (r.m_vars.empty()) { + + if (r.m_vars.empty()) continue; - } - if (r.m_vars.size() == 1 && r.m_vars[0].m_coeff.is_neg() && r.m_type != opt::t_mod) { + 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); }