From b62d73f20958fe295d2ec728dd10e2d19218d0c9 Mon Sep 17 00:00:00 2001 From: Nikolaj Bjorner Date: Wed, 13 Jun 2018 07:55:45 -0700 Subject: [PATCH] first round for combined mbi Signed-off-by: Nikolaj Bjorner --- src/math/simplex/model_based_opt.cpp | 95 ++++++++++++------ src/math/simplex/model_based_opt.h | 6 +- src/qe/qe_arith.cpp | 53 ++++++---- src/qe/qe_arith.h | 6 ++ src/qe/qe_mbi.cpp | 139 ++++++++++++++++++++++++++- src/qe/qe_mbi.h | 30 +++--- src/test/model_based_opt.cpp | 32 ++++++ 7 files changed, 297 insertions(+), 64 deletions(-) diff --git a/src/math/simplex/model_based_opt.cpp b/src/math/simplex/model_based_opt.cpp index 2b175506f..56d51648d 100644 --- a/src/math/simplex/model_based_opt.cpp +++ b/src/math/simplex/model_based_opt.cpp @@ -36,24 +36,22 @@ std::ostream& operator<<(std::ostream& out, opt::ineq_type ie) { namespace opt { /** - * Convert a row ax + coeffs + coeff = 0 into a definition for x - * x = -(coeffs + coeff)/a - * For rows ax + coeffs + coeff < 0 convert into - * x = -(coeffs + coeff - a)/a + * 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) { - rational c = r.get_coefficient(x); - bool is_pos = c.is_pos(); for (var const & v : r.m_vars) { if (v.m_id != x) { m_vars.push_back(v); - if (is_pos) m_vars.back().m_coeff.neg(); + } + else { + m_div = -v.m_coeff; } } - m_coeff = r.m_coeff; - if (is_pos) m_coeff.neg(); - if (r.m_type == opt::t_lt) m_coeff += abs(c); - m_div = abs(c); + m_coeff = r.m_coeff - r.m_value; normalize(); SASSERT(m_div.is_pos()); } @@ -129,6 +127,9 @@ namespace opt { g = gcd(g, abs(v.m_coeff)); if (g.is_one()) break; } + if (m_div.is_neg()) { + g.neg(); + } if (!g.is_one()) { for (var& v : m_vars) { v.m_coeff /= g; @@ -138,7 +139,6 @@ namespace opt { } } - model_based_opt::model_based_opt() { m_rows.push_back(row()); } @@ -163,7 +163,7 @@ namespace opt { PASSERT(index == 0 || m_var2row_ids[vars[i].m_id].contains(index)); } - PASSERT(r.m_value == get_row_value(r)); + 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()); @@ -317,7 +317,7 @@ namespace opt { << " eps: " << eps << " ", r); ); m_var2value[x] = new_x_val; - r.m_value = get_row_value(r); + r.m_value = eval(r); SASSERT(invariant(bound_trail[i], r)); } @@ -327,7 +327,7 @@ namespace opt { unsigned_vector const& row_ids = m_var2row_ids[x]; for (unsigned row_id : row_ids) { row & r = m_rows[row_id]; - r.m_value = get_row_value(r); + r.m_value = eval(r); SASSERT(invariant(row_id, r)); } } @@ -385,19 +385,32 @@ namespace opt { 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::get_row_value(row const& r) const { + 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 * m_var2value[v.m_id]; + val += v.m_coeff * eval(v.m_id); } return val; } rational model_based_opt::get_coefficient(unsigned row_id, unsigned var_id) const { - row const& r = m_rows[row_id]; - return r.get_coefficient(var_id); + return m_rows[row_id].get_coefficient(var_id); } rational model_based_opt::row::get_coefficient(unsigned var_id) const { @@ -639,7 +652,7 @@ namespace opt { row& r1 = m_rows[row_id1]; row const& r2 = m_rows[row_id2]; unsigned i = 0, j = 0; - for(; i < r1.m_vars.size() || j < r2.m_vars.size(); ) { + 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.c_ptr() + i); break; @@ -707,11 +720,18 @@ namespace opt { } void model_based_opt::display(std::ostream& out, vector const& vars, rational const& coeff) { - for (unsigned i = 0; i < vars.size(); ++i) { - if (i > 0 && vars[i].m_coeff.is_pos()) { + unsigned i = 0; + for (var const& v : vars) { + if (i > 0 && v.m_coeff.is_pos()) { out << "+ "; } - out << vars[i].m_coeff << "* v" << vars[i].m_id << " "; + ++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 << " "; @@ -921,7 +941,7 @@ namespace opt { return solve_for(eq_row, x, compute_def); } - def result; + 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; @@ -935,6 +955,11 @@ namespace opt { else if (glb_index != UINT_MAX) { result = solve_for(glb_index, x, true); } + else { + result = def(); + result.m_coeff = eval(x); + } + SASSERT(eval(result) == eval(x)); } else { for (unsigned row_id : lub_rows) retire_row(row_id); @@ -946,9 +971,19 @@ namespace opt { SASSERT(lub_index != UINT_MAX); SASSERT(glb_index != UINT_MAX); if (compute_def) { +#if 0 def d1(m_rows[lub_index], x); - def d2(m_rows[lub_index], x); + def d2(m_rows[glb_index], x); result = (d1 + d2)/2; +#else + if (lub_size <= glb_size) { + result = def(m_rows[lub_index], x); + } + else { + result = def(m_rows[glb_index], x); + } +#endif + SASSERT(eval(result) == eval(x)); } // The number of matching lower and upper bounds is small. @@ -1045,8 +1080,9 @@ namespace opt { } def result = project(y, compute_def); if (compute_def) { - result = (result * D) + u; + result = (result * D) + u; } + SASSERT(!compute_def || eval(result) == eval(x)); return result; } @@ -1139,8 +1175,11 @@ namespace opt { } } } - // TBD: -t div a - def result(m_rows[row_id1], x); + def result; + if (compute_def) { + result = def(m_rows[row_id1], x); + SASSERT(eval(result) == eval(x)); + } retire_row(row_id1); return result; } diff --git a/src/math/simplex/model_based_opt.h b/src/math/simplex/model_based_opt.h index 268d3d81d..a1137d2fa 100644 --- a/src/math/simplex/model_based_opt.h +++ b/src/math/simplex/model_based_opt.h @@ -100,7 +100,11 @@ namespace opt { rational get_coefficient(unsigned row_id, unsigned var_id) const; - rational get_row_value(row const& r) const; + rational eval(row const& r) const; + + rational eval(unsigned x) const; + + rational eval(def const& d) const; void resolve(unsigned row_src, rational const& a1, unsigned row_dst, unsigned x); diff --git a/src/qe/qe_arith.cpp b/src/qe/qe_arith.cpp index 0b5ad72ca..1709711d7 100644 --- a/src/qe/qe_arith.cpp +++ b/src/qe/qe_arith.cpp @@ -37,6 +37,7 @@ namespace qe { ast_manager& m; arith_util a; + bool m_check_purified; // check that variables are properly pure void insert_mul(expr* x, rational const& v, obj_map& ts) { TRACE("qe", tout << "Adding variable " << mk_pp(x, m) << " " << v << "\n";); @@ -264,7 +265,7 @@ namespace qe { } imp(ast_manager& m): - m(m), a(m) {} + m(m), a(m), m_check_purified(true) {} ~imp() {} @@ -347,17 +348,23 @@ namespace qe { tids.insert(v, mbo.add_var(r, a.is_int(v))); } } - for (expr* fml : fmls) { - mark_rec(fmls_mark, fml); + if (m_check_purified) { + 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); + } + } } + ptr_vector index2expr; for (auto& kv : tids) { - expr* e = kv.m_key; - if (!var_mark.is_marked(e)) { - mark_rec(fmls_mark, e); - } - index2expr.setx(kv.m_value, e, 0); + index2expr.setx(kv.m_value, kv.m_key, 0); } + j = 0; unsigned_vector real_vars; for (app* v : vars) { @@ -369,6 +376,7 @@ namespace qe { } } vars.shrink(j); + TRACE("qe", tout << "remaining vars: " << vars << "\n"; for (unsigned v : real_vars) { tout << "v" << v << " " << mk_pp(index2expr[v], m) << "\n"; @@ -379,10 +387,9 @@ namespace qe { vector rows; mbo.get_live_rows(rows); - for (unsigned i = 0; i < rows.size(); ++i) { + for (row const& r : rows) { expr_ref_vector ts(m); expr_ref t(m), s(m), val(m); - row const& r = rows[i]; if (r.m_vars.size() == 0) { continue; } @@ -411,25 +418,19 @@ namespace qe { } ts.push_back(t); } + t = mk_add(ts); s = a.mk_numeral(-r.m_coeff, a.is_int(t)); - if (ts.size() == 1) { - t = ts[0].get(); - } - else { - t = a.mk_add(ts.size(), ts.c_ptr()); - } 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: { + 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)); break; } - } fmls.push_back(t); val = eval(t); CTRACE("qe", !m.is_true(val), tout << "Evaluated " << t << " to " << val << "\n";); @@ -447,19 +448,29 @@ namespace qe { ts.push_back(var2expr(index2expr, v)); } ts.push_back(a.mk_numeral(d.m_coeff, is_int)); - t = a.mk_add(ts.size(), ts.c_ptr()); + 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)); } + SASSERT(eval(t) == eval(x)); result.push_back(def(expr_ref(x, m), t)); } } return result; } + expr_ref mk_add(expr_ref_vector const& ts) { + if (ts.size() == 1) { + return expr_ref(ts.get(0), m); + } + else { + return expr_ref(a.mk_add(ts.size(), ts.c_ptr()), m); + } + } + opt::inf_eps maximize(expr_ref_vector const& fmls0, model& mdl, app* t, expr_ref& ge, expr_ref& gt) { SASSERT(a.is_real(t)); expr_ref_vector fmls(fmls0); @@ -577,6 +588,10 @@ namespace qe { return m_imp->project(model, vars, lits); } + void arith_project_plugin::set_check_purified(bool check_purified) { + m_imp->m_check_purified = check_purified; + } + bool arith_project_plugin::solve(model& model, app_ref_vector& vars, expr_ref_vector& lits) { return m_imp->solve(model, vars, lits); } diff --git a/src/qe/qe_arith.h b/src/qe/qe_arith.h index c01d7bbb6..b55e63fcf 100644 --- a/src/qe/qe_arith.h +++ b/src/qe/qe_arith.h @@ -33,6 +33,12 @@ namespace qe { vector project(model& model, app_ref_vector& vars, expr_ref_vector& lits) override; opt::inf_eps maximize(expr_ref_vector const& fmls, model& mdl, app* t, expr_ref& ge, expr_ref& gt); + + /** + * \brief check if formulas are purified, or leave it to caller to ensure that + * arithmetic variables nested under foreign functions are handled properly. + */ + void set_check_purified(bool check_purified); }; bool arith_project(model& model, app* var, expr_ref_vector& lits); diff --git a/src/qe/qe_mbi.cpp b/src/qe/qe_mbi.cpp index ae989ba02..20eaf9675 100644 --- a/src/qe/qe_mbi.cpp +++ b/src/qe/qe_mbi.cpp @@ -77,10 +77,12 @@ Notes: #include "ast/ast_util.h" #include "ast/for_each_expr.h" #include "ast/rewriter/bool_rewriter.h" +#include "ast/arith_decl_plugin.h" #include "model/model_evaluator.h" #include "solver/solver.h" #include "qe/qe_mbi.h" #include "qe/qe_term_graph.h" +#include "qe/qe_arith.h" namespace qe { @@ -141,7 +143,7 @@ namespace qe { } // ------------------------------- - // euf_mbi, TBD + // euf_mbi struct euf_mbi_plugin::is_atom_proc { ast_manager& m; @@ -235,6 +237,140 @@ namespace qe { } + // ------------------------------- + // euf_arith_mbi + + struct euf_arith_mbi_plugin::is_atom_proc { + ast_manager& m; + expr_ref_vector& m_atoms; + is_atom_proc(expr_ref_vector& atoms): m(atoms.m()), m_atoms(atoms) {} + void operator()(app* a) { + if (m.is_eq(a)) { + m_atoms.push_back(a); + } + else if (m.is_bool(a) && a->get_family_id() != m.get_basic_family_id()) { + m_atoms.push_back(a); + } + } + void operator()(expr*) {} + }; + + struct euf_arith_mbi_plugin::is_arith_var_proc { + ast_manager& m; + app_ref_vector& m_vars; + arith_util arith; + obj_hashtable m_exclude; + is_arith_var_proc(app_ref_vector& vars, func_decl_ref_vector const& excl): + m(vars.m()), m_vars(vars), arith(m) { + for (func_decl* f : excl) m_exclude.insert(f); + } + void operator()(app* a) { + if (arith.is_int_real(a) && a->get_family_id() != a->get_family_id() && !m_exclude.contains(a->get_decl())) { + m_vars.push_back(a); + } + } + void operator()(expr*) {} + + }; + + euf_arith_mbi_plugin::euf_arith_mbi_plugin(solver* s, solver* sNot): + m(s->get_manager()), + m_atoms(m), + m_solver(s), + m_dual_solver(sNot) { + params_ref p; + p.set_bool("core.minimize", true); + m_solver->updt_params(p); + m_dual_solver->updt_params(p); + expr_ref_vector fmls(m); + m_solver->get_assertions(fmls); + expr_fast_mark1 marks; + is_atom_proc proc(m_atoms); + for (expr* e : fmls) { + quick_for_each_expr(proc, marks, e); + } + } + + mbi_result euf_arith_mbi_plugin::operator()(func_decl_ref_vector const& vars, expr_ref_vector& lits, model_ref& mdl) { + lbool r = m_solver->check_sat(lits); + + // populate set of arithmetic variables to be projected. + app_ref_vector avars(m); + is_arith_var_proc _proc(avars, vars); + for (expr* l : lits) quick_for_each_expr(_proc, l); + switch (r) { + case l_false: + lits.reset(); + m_solver->get_unsat_core(lits); + // optionally minimize core using superposition. + return mbi_unsat; + case l_true: { + m_solver->get_model(mdl); + model_evaluator mev(*mdl.get()); + lits.reset(); + for (expr* e : m_atoms) { + if (mev.is_true(e)) { + lits.push_back(e); + } + else if (mev.is_false(e)) { + lits.push_back(m.mk_not(e)); + } + } + TRACE("qe", tout << "atoms from model: " << lits << "\n";); + r = m_dual_solver->check_sat(lits); + expr_ref_vector core(m); + term_graph tg(m); + switch (r) { + case l_false: { + // use the dual solver to find a 'small' implicant + m_dual_solver->get_unsat_core(core); + TRACE("qe", tout << "core: " << core << "\n";); + lits.reset(); + lits.append(core); + arith_util a(m); + + // 1. project arithmetic variables using mdl that satisfies core. + // ground any remaining arithmetic variables using model. + arith_project_plugin ap(m); + ap.set_check_purified(false); + + auto defs = ap.project(*mdl.get(), avars, lits); + // 2. Add the projected definitions to the remaining (euf) literals + for (auto const& def : defs) { + lits.push_back(m.mk_eq(def.var, def.term)); + } + + // 3. Project the remaining literals with respect to EUF. + tg.set_vars(vars, false); + tg.add_lits(lits); + lits.reset(); + lits.append(tg.project(*mdl)); + TRACE("qe", tout << "project: " << lits << "\n";); + return mbi_sat; + } + case l_undef: + return mbi_undef; + case l_true: + UNREACHABLE(); + return mbi_undef; + } + return mbi_sat; + } + default: + // TBD: if not running solver to completion, then: + // 1. extract unit literals from m_solver. + // 2. run a cc over the units + // 3. extract equalities or other assignments over the congruence classes + // 4. ensure that at least some progress is made over original lits. + return mbi_undef; + } + } + + void euf_arith_mbi_plugin::block(expr_ref_vector const& lits) { + m_solver->assert_expr(mk_not(mk_and(lits))); + } + + /** -------------------------------------------------------------- * ping-pong interpolation of Gurfinkel & Vizel * compute a binary interpolant. @@ -318,4 +454,5 @@ namespace qe { } } } + }; diff --git a/src/qe/qe_mbi.h b/src/qe/qe_mbi.h index 671099015..6d0d8dcf6 100644 --- a/src/qe/qe_mbi.h +++ b/src/qe/qe_mbi.h @@ -86,29 +86,29 @@ namespace qe { void block(expr_ref_vector const& lits) override; }; + class euf_arith_mbi_plugin : mbi_plugin { + ast_manager& m; + expr_ref_vector m_atoms; + solver_ref m_solver; + solver_ref m_dual_solver; + struct is_atom_proc; + struct is_arith_var_proc; + public: + euf_arith_mbi_plugin(solver* s, solver* sNot); + ~euf_arith_mbi_plugin() override {} + mbi_result operator()(func_decl_ref_vector const& vars, expr_ref_vector& lits, model_ref& mdl) override; + void block(expr_ref_vector const& lits) override; + }; + /** * use cases for interpolation. */ class interpolator { - ast_manager& m; + ast_manager& m; public: interpolator(ast_manager& m):m(m) {} lbool pingpong(mbi_plugin& a, mbi_plugin& b, func_decl_ref_vector const& vars, expr_ref& itp); lbool pogo(mbi_plugin& a, mbi_plugin& b, func_decl_ref_vector const& vars, expr_ref& itp); }; -#if 0 - TBD some other scenario, e.g., Farkas, Cute/Beautiful/maybe even Selfless - - class solver_mbi_plugin : public mbi_plugin { - solver_ref m_solver; - public: - solver_mbi_plugin(solver* s); - ~solver_mbi_plugin() override {} - mbi_result operator()(func_decl_ref_vector const& vars, expr_ref_vector& lits, model_ref& mdl) override; - void block(expr_ref_vector const& lits) override; - }; - - -#endif }; diff --git a/src/test/model_based_opt.cpp b/src/test/model_based_opt.cpp index f2f6004ff..1bae7435f 100644 --- a/src/test/model_based_opt.cpp +++ b/src/test/model_based_opt.cpp @@ -360,6 +360,37 @@ static void test10() { mbo.display(std::cout); } +static void test11() { + { + opt::model_based_opt mbo; + unsigned s = mbo.add_var(rational(2), true); + unsigned a = mbo.add_var(rational(1), true); + unsigned t = mbo.add_var(rational(3), true); + + add_ineq(mbo, s, 1, a, -2, 0, opt::t_le); // s - 2a <= 0 + add_ineq(mbo, a, 2, t, -1, 0, opt::t_le); // 2a - t <= 0 + + mbo.display(std::cout); + display_project(mbo.project(1, &a, true)); + mbo.display(std::cout); + } + + { + opt::model_based_opt mbo; + unsigned s = mbo.add_var(rational(3), true); + unsigned a = mbo.add_var(rational(2), true); + unsigned t = mbo.add_var(rational(4), true); + + add_ineq(mbo, s, 1, a, -2, 0, opt::t_le); // s - 2a <= 0 + add_ineq(mbo, a, 2, t, -1, 0, opt::t_le); // 2a - t <= 0 + + mbo.display(std::cout); + display_project(mbo.project(1, &a, true)); + mbo.display(std::cout); + } + +} + // test with mix of upper and lower bounds void tst_model_based_opt() { @@ -374,4 +405,5 @@ void tst_model_based_opt() { test7(); test8(); test9(); + test11(); }