first round for combined mbi

Signed-off-by: Nikolaj Bjorner <nbjorner@microsoft.com>

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
+     * Convert a row ax + coeffs + coeff = value into a definition for x
-     *    x  = -(coeffs + coeff)/a 
+     *    x  = (value - coeffs - coeff)/a 
-     * For rows ax + coeffs + coeff < 0 convert into
+     * as backdrop we have existing assignments to x and other variables that 
-     *    x  = -(coeffs + coeff - a)/a 
+     * 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;
+        m_coeff = r.m_coeff - r.m_value;
-        if (is_pos) m_coeff.neg();
-        if (r.m_type == opt::t_lt) m_coeff += abs(c);
-        m_div = abs(c);
         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<var> 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<var> 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 m_rows[row_id].get_coefficient(var_id);
-        return r.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<var> const& vars, rational const& coeff) {
-        for (unsigned i = 0; i < vars.size(); ++i) {
+        unsigned i = 0;
-            if (i > 0 && vars[i].m_coeff.is_pos()) {
+        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;
-        def result(m_rows[row_id1], x);
+        if (compute_def) {
+            result = def(m_rows[row_id1], x);
+            SASSERT(eval(result) == eval(x));
+        }
         retire_row(row_id1);
         return result;
     }

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);
 

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<expr, rational>& 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) {
+            if (m_check_purified) {
-                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);
+                    }
+                }
             }
+
             ptr_vector<expr> index2expr;
             for (auto& kv : tids) {
-                expr* e = kv.m_key;
+                index2expr.setx(kv.m_value, kv.m_key, 0);
-                if (!var_mark.is_marked(e)) {
-                    mark_rec(fmls_mark, e);
-                }
-                index2expr.setx(kv.m_value, e, 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<row> 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);
     }

src/qe/qe_arith.h

@@ -33,6 +33,12 @@ namespace qe {
         vector<def> 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);

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<func_decl> 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 {
             }
         }
     }
+
 };

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
 };

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();
 }