adding local optimization to qsat

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

src/qe/qe_arith.cpp

@@ -20,6 +20,7 @@ Revision History:
 --*/
 
 #include "qe_arith.h"
+#include "qe_mbp.h"
 #include "ast_util.h"
 #include "arith_decl_plugin.h"
 #include "ast_pp.h"
@@ -50,13 +51,184 @@ namespace qe {
         return is_divides(a, e1, e2, k, t) || is_divides(a, e2, e1, k, t);
     }
 
+    enum ineq_type {
+        t_eq,
+        t_lt,
+        t_le
+    };
+
+    struct tableau {
+        struct var {
+            unsigned m_id;
+            rational m_coeff;
+            var(unsigned id, rational const& c): m_id(id), m_coeff(c) {}
+        };
+        struct row {
+            vector<var> m_vars;         // variables with coefficients
+            rational    m_coeff;        // constant in inequality
+            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.
+        };
+        vector<row>             m_rows;
+        vector<unsigned_vector> m_var2rows;
+        vector<rational>        m_var2value;
+        row                     m_objective;
+        
+        void invariant() {
+            // variables in each row are sorted.
+        }
+        
+        mbp::bound_type maximize(rational& value) {
+            // tbd
+            return mbp::unbounded;
+        }
+        
+        rational get_coefficient(unsigned row_id, unsigned var_id) {
+            row const& r = m_rows[row_id];
+            unsigned lo = 0, hi = r.m_vars.size();
+            while (lo < hi) {
+                unsigned mid = lo + (hi - lo)/2;
+                SASSERT(mid < hi);
+                unsigned id = r.m_vars[mid].m_id;
+                if (id == var_id) {
+                    lo = mid;
+                    break;
+                }
+                if (id < var_id) {
+                    lo = mid + 1;
+                }
+                else {
+                    hi = mid - 1;
+                }
+            }
+            unsigned id = r.m_vars[lo].m_id;
+            if (id == var_id) {
+                return r.m_vars[lo].m_coeff;
+            }
+            else {
+                return rational::zero();
+            }
+        }
+
+        void resolve(unsigned row_id1, unsigned row_id2, unsigned x) {
+            // row1 is of the form a1*x + t1 <~ 0
+            // row2 is of the form a2*x + t2 <~ 0
+            // assume that a1, a2 have the same sign.
+            // if a1 is positive, then val(t1*a2/a1) <= val(t2*a1/a2)
+            //   replace row2 with the new inequality of the form:
+            //   t1 - a1*t2/a2 <~~ 0
+            //   where <~~ is strict if either <~1 or <~2 is strict.
+            // if a1 is negative, then ....
+            //   
+        }
+
+        void multiply(rational const& c, unsigned row_id) {
+            if (c.is_one()) {
+                return;
+            }
+            row& r = m_rows[row_id];
+            SASSERT(r.m_alive);
+            for (unsigned i = 0; i < r.m_vars.size(); ++i) {
+                r.m_vars[i].m_coeff *= c;
+            }
+            r.m_coeff *= c;
+            r.m_value *= c;
+        }
+
+        // subtract row2 from row1, store result in row2
+    
+        vector<var> m_new_vars;
+
+        void subtract(unsigned row_id1, unsigned row_id2) {
+            m_new_vars.reset();
+            row const& r1 = m_rows[row_id1];
+            row& r2 = m_rows[row_id2];
+            unsigned i = 0, j = 0;
+            for(; i < r1.m_vars.size() || j < r2.m_vars.size(); ) {
+                if (j == r2.m_vars.size()) {
+                    for (; i < r1.m_vars.size(); ++i) {
+                        m_new_vars.push_back(r1.m_vars[i]);
+                        m_var2rows[r1.m_vars[i].m_id].push_back(row_id2);
+                    }
+                }
+                else 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.neg();
+                    }
+                }
+                else {
+                    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 -= 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]);
+                        m_var2rows[r1.m_vars[i].m_id].push_back(row_id2);
+                        ++i;                        
+                    }
+                    else {
+                        m_new_vars.push_back(r2.m_vars[j]);
+                        m_new_vars.back().m_coeff.neg();
+                        ++j;                        
+                    }
+                }                
+            }
+            r2.m_coeff.neg();
+            r2.m_coeff += r1.m_coeff;
+            r2.m_vars.swap(m_new_vars);
+            r2.m_value.neg();
+            r2.m_value += r1.m_value;
+            if (r1.m_type == t_lt) {
+                r2.m_type = t_lt;
+            }
+        }
+
+        void display(std::ostream& out) const {
+            for (unsigned i = 0; i < m_rows.size(); ++i) {
+                display(out, m_rows[i]);
+            }
+        }        
+
+        void display(std::ostream& out, row const& r) const {
+            vector<var> const& vars = r.m_vars;
+            for (unsigned i = 0; i < vars.size(); ++i) {
+                if (i > 0 && vars[i].m_coeff.is_pos()) {
+                    out << "+ ";
+                }
+                out << vars[i].m_coeff << "* v" << vars[i].m_id << " ";                
+            }
+            out << r.m_coeff;
+            switch (r.m_type) {
+            case t_eq:
+                out << " = 0\n";
+                break;
+            case t_lt:
+                out << " < 0\n";
+                break;
+            case t_le:
+                out << " <= 0\n";
+                break;
+            }
+        }
+    };
+
+#if 0
+        obj_map<expr, unsigned> m_expr2var;
+        ptr_vector<expr>        m_var2expr;
+
+#endif
+    
     struct arith_project_plugin::imp {
 
-        enum ineq_type {
-            t_eq,
-            t_lt,
-            t_le
-        };
         ast_manager&      m;
         arith_util        a;
         th_rewriter       m_rw;
@@ -84,6 +256,62 @@ namespace qe {
             }            
         }
 
+        void insert_mul(expr* x, rational const& v, obj_map<expr, rational>& ts)
+        {
+            rational w;
+            if (ts.find(x, w)) {
+                ts.insert(x, w + v);
+            }
+            else {
+                ts.insert(x, v); 
+            }
+        }
+
+        void linearize(model& model, rational const& mul, expr* t, rational& c, obj_map<expr, rational>& ts) {
+            expr* t1, *t2, *t3;
+            rational mul1;
+            expr_ref val(m);
+            if (a.is_mul(t, t1, t2) && is_numeral(model, t1, mul1)) {
+                linearize(model, mul* mul1, t2, c, ts);
+            }
+            else if (a.is_mul(t, t1, t2) && is_numeral(model, t2, mul1)) {
+                linearize(model, mul* mul1, t1, c, ts);
+            }
+            else if (a.is_add(t)) {
+                app* ap = to_app(t);
+                for (unsigned i = 0; i < ap->get_num_args(); ++i) {
+                    linearize(model, mul, ap->get_arg(i), c, ts);
+                }
+            }
+            else if (a.is_sub(t, t1, t2)) {
+                linearize(model, mul, t1, c, ts);
+                linearize(model, -mul, t2, c, ts);
+            }
+            else if (a.is_uminus(t, t1)) {
+                linearize(model, -mul, t1, c, ts);
+            }
+            else if (a.is_numeral(t, mul1)) {
+                c += mul*mul1;
+            }
+            else if (extract_mod(model, t, val)) {
+                insert_mul(val, mul, ts);
+            }
+            else if (m.is_ite(t, t1, t2, t3)) {
+                VERIFY(model.eval(t1, val));
+                SASSERT(m.is_true(val) || m.is_false(val));
+                TRACE("qe", tout << mk_pp(t1, m) << " := " << val << "\n";);
+                if (m.is_true(val)) {
+                    linearize(model, mul, t2, c, ts);
+                }
+                else {
+                    linearize(model, mul, t3, c, ts);
+                }
+            }
+            else {
+                insert_mul(t, mul, ts);
+            }
+        }
+
         void is_linear(model& model, rational const& mul, expr* t, rational& c, expr_ref_vector& ts) {
             expr* t1, *t2, *t3;
             rational mul1;
@@ -853,6 +1081,19 @@ namespace qe {
             }
             return true;
         }
+
+        mbp::bound_type maximize(expr_ref_vector const& fmls, model& mdl, app* t, expr_ref& value, expr_ref& bound) {
+            obj_map<expr, rational> ts;
+            rational c(0), mul(1);
+            linearize(mdl, mul, t, c, ts);              
+                
+            
+            // pick variables one by one from ts.
+            // m_var = alloc(contains_app, m, v);
+            // perform upper or lower projection depending on sign of v.
+            // 
+            return mbp::unbounded;
+        }
     };
 
     arith_project_plugin::arith_project_plugin(ast_manager& m) {
@@ -875,6 +1116,10 @@ namespace qe {
         return m_imp->a.get_family_id();
     }
 
+    mbp::bound_type arith_project_plugin::maximize(expr_ref_vector const& fmls, model& mdl, app* t, expr_ref& value, expr_ref& bound) {
+        return m_imp->maximize(fmls, mdl, t, value, bound);
+    }
+
     bool arith_project(model& model, app* var, expr_ref_vector& lits) {
         ast_manager& m = lits.get_manager();
         arith_project_plugin ap(m);

src/qe/qe_arith.h

@@ -29,6 +29,7 @@ namespace qe {
         virtual bool operator()(model& model, app* var, app_ref_vector& vars, expr_ref_vector& lits);
         virtual bool solve(model& model, app_ref_vector& vars, expr_ref_vector& lits);
         virtual family_id get_family_id();
+        mbp::bound_type maximize(expr_ref_vector const& fmls, model& mdl, app* t, expr_ref& value, expr_ref& bound);
     };
 
     bool arith_project(model& model, app* var, expr_ref_vector& lits);

src/qe/qe_mbp.cpp

@@ -213,6 +213,11 @@ class mbp::impl {
 public:
 
 
+    mbp::bound_type maximize(expr_ref_vector const& fmls, model& mdl, app* t, expr_ref& value, expr_ref& bound) {
+        arith_project_plugin arith(m);
+        return arith.maximize(fmls, mdl, t, value, bound);
+    }
+
     void extract_literals(model& model, expr_ref_vector& fmls) {
         expr_ref val(m);
         for (unsigned i = 0; i < fmls.size(); ++i) {
@@ -415,3 +420,7 @@ void mbp::solve(model& model, app_ref_vector& vars, expr_ref_vector& fmls) {
 void mbp::extract_literals(model& model, expr_ref_vector& lits) {
     m_impl->extract_literals(model, lits);
 }
+
+mbp::bound_type mbp::maximize(expr_ref_vector const& fmls, model& mdl, app* t, expr_ref& value, expr_ref& bound) {
+    return m_impl->maximize(fmls, mdl, t, value, bound);
+}

src/qe/qe_mbp.h

@@ -70,6 +70,17 @@ namespace qe {
            Extract literals from formulas based on model.
          */
         void extract_literals(model& model, expr_ref_vector& lits);
+
+        /**
+           \brief 
+           Maximize objective t under current model for constraints in fmls.
+         */
+        enum bound_type {
+            unbounded,
+            strict,
+            non_strict
+        };
+        bound_type maximize(expr_ref_vector const& fmls, model& mdl, app* t, expr_ref& value, expr_ref& bound);
     };
 }
 

src/qe/qsat.cpp

@@ -537,6 +537,13 @@ namespace qe {
         }
     };
 
+    enum qsat_mode {
+        qsat_qe,
+        qsat_qe_rec,
+        qsat_sat,
+        qsat_maximize
+    };
+    
     class qsat : public tactic {
         
         struct stats {
@@ -559,8 +566,7 @@ namespace qe {
         vector<app_ref_vector>     m_vars;       // variables from alternating prefixes.
         unsigned                   m_level;
         model_ref                  m_model;
-        bool                       m_qelim;       // perform quantifier elimination
+        qsat_mode                  m_mode;
-        bool                       m_force_elim;  // force elimination of variables during projection.
         app_ref_vector             m_avars;       // variables to project
         app_ref_vector             m_free_vars;
 
@@ -584,12 +590,12 @@ namespace qe {
                     SASSERT(validate_model(asms));
                     TRACE("qe", k.display(tout); display(tout << "\n", *m_model.get()); display(tout, asms); );
                     push();
-                break;
+                    break;
                 case l_false:
                     switch (m_level) {
                     case 0: return l_false;
                     case 1: 
-                        if (!m_qelim) return l_true; 
+                        if (m_mode == qsat_sat) return l_true; 
                         if (m_model.get()) {
                             project_qe(asms);
                         }
@@ -672,7 +678,7 @@ namespace qe {
             m_pred_abs.get_free_vars(fml, vars);
             m_vars.push_back(vars);
             vars.reset();
-            if (m_qelim) {
+            if (m_mode != qsat_sat) {
                 is_forall = true;
                 hoist.pull_quantifier(is_forall, fml, vars);
                 m_vars.push_back(vars);
@@ -858,12 +864,18 @@ namespace qe {
             get_core(core, m_level);
             SASSERT(validate_core(core));
             get_vars(m_level);
-            m_mbp(m_force_elim, m_avars, mdl, core);
+            m_mbp(force_elim(), m_avars, mdl, core);
-            fml = negate_core(core);
+            if (m_mode == qsat_maximize) {
-            add_assumption(fml);
+                maximize(core, mdl);
-            m_answer.push_back(fml);
+                pop(1);
-            m_free_vars.append(m_avars);
+            }
-            pop(1);
+            else {
+                fml = negate_core(core);
+                add_assumption(fml);
+                m_answer.push_back(fml);
+                m_free_vars.append(m_avars);
+                pop(1);
+            }
         }
                 
         void project(expr_ref_vector& core) {
@@ -878,7 +890,7 @@ namespace qe {
             
             get_vars(m_level-1);
             SASSERT(validate_project(mdl, core));
-            m_mbp(m_force_elim, m_avars, mdl, core);
+            m_mbp(force_elim(), m_avars, mdl, core);
             m_free_vars.append(m_avars);
             fml = negate_core(core);
             unsigned num_scopes = 0;
@@ -889,7 +901,7 @@ namespace qe {
             if (level.max() == UINT_MAX) {
                 num_scopes = 2*(m_level/2);
             }
-            else if (m_qelim && !m_force_elim) {
+            else if (m_mode == qsat_qe_rec) {
                 num_scopes = 2;
             }
             else {
@@ -903,7 +915,7 @@ namespace qe {
             
             pop(num_scopes); 
             TRACE("qe", tout << "backtrack: " << num_scopes << " new level: " << m_level << "\nproject:\n" << core << "\n|->\n" << fml << "\n";);
-            if (m_level == 0 && m_qelim) {
+            if (m_level == 0 && m_mode != qsat_sat) {
                 add_assumption(fml);
             }
             else {
@@ -919,10 +931,14 @@ namespace qe {
             }
         } 
         
-        expr_ref negate_core(expr_ref_vector& core) {
+        expr_ref negate_core(expr_ref_vector const& core) {
             return ::push_not(::mk_and(core));
         }
 
+        bool force_elim() const {
+            return m_mode != qsat_qe_rec;
+        }
+        
         expr_ref elim_rec(expr* fml) {
             expr_ref tmp(m);
             expr_ref_vector     trail(m);
@@ -1135,7 +1151,7 @@ namespace qe {
 
     public:
         
-        qsat(ast_manager& m, params_ref const& p, bool qelim, bool force_elim):
+        qsat(ast_manager& m, params_ref const& p, qsat_mode mode):
             m(m),
             m_mbp(m),
             m_fa(m),
@@ -1144,10 +1160,10 @@ namespace qe {
             m_answer(m),
             m_asms(m),
             m_level(0),
-            m_qelim(qelim),
+            m_mode(mode),
-            m_force_elim(force_elim),
             m_avars(m),
-            m_free_vars(m)
+            m_free_vars(m),
+            m_value(m)
         {
             reset();
         }
@@ -1182,7 +1198,7 @@ namespace qe {
             // fail if cores.  (TBD)
             // fail if proofs. (TBD)
             
-            if (!m_force_elim) {
+            if (m_mode == qsat_qe_rec) {
                 fml = elim_rec(fml);
                 in->reset();
                 in->inc_depth();
@@ -1193,7 +1209,7 @@ namespace qe {
                 
             reset();
             TRACE("qe", tout << fml << "\n";);
-            if (m_qelim) {
+            if (m_mode != qsat_sat) {
                 fml = push_not(fml);
             }
             hoist(fml);
@@ -1211,11 +1227,12 @@ namespace qe {
             case l_false:
                 in->reset();
                 in->inc_depth();
-                if (m_qelim) {
+                if (m_mode == qsat_qe) {
                     fml = ::mk_and(m_answer);
                     in->assert_expr(fml);
                 }
                 else {
+                    SASSERT(m_mode == qsat_sat);
                     in->assert_expr(m.mk_false());
                 }
                 result.push_back(in.get());
@@ -1262,13 +1279,111 @@ namespace qe {
         }
         
         tactic * translate(ast_manager & m) {
-            return alloc(qsat, m, m_params, m_qelim, m_force_elim);
+            return alloc(qsat, m, m_params, m_mode);
         }        
+
+        app*             m_objective;
+        expr_ref         m_value;
+        mbp::bound_type  m_bound;
+        bool             m_was_sat;
+
+        lbool maximize(expr_ref_vector const& fmls, app* t, expr_ref& value, mbp::bound_type& bound) {
+            expr_ref_vector defs(m);
+            expr_ref fml = negate_core(fmls);
+            hoist(fml);
+            m_objective = t;
+            m_value = 0;
+            m_bound = mbp::unbounded;
+            m_was_sat = false;
+            m_pred_abs.abstract_atoms(fml, defs);
+            fml = m_pred_abs.mk_abstract(fml);
+            m_ex.assert_expr(mk_and(defs));
+            m_fa.assert_expr(mk_and(defs));
+            m_ex.assert_expr(fml);
+            m_fa.assert_expr(m.mk_not(fml));
+            lbool is_sat = check_sat();
+            switch (is_sat) {
+            case l_false:
+                if (!m_was_sat) {
+                    return l_false;
+                }
+                break;
+            case l_true:
+                UNREACHABLE();
+                break;
+            case l_undef:
+                std::string s = m_ex.k().last_failure_as_string();
+                if (s == "ok") {
+                    s = m_fa.k().last_failure_as_string();
+                }
+                throw tactic_exception(s.c_str()); 
+            }        
+            value = m_value;
+            bound = m_bound;
+            return l_true;
+        }
+
+        void maximize(expr_ref_vector const& core, model& mdl) {
+            TRACE("qe", tout << "maximize: " << core << "\n";);
+            m_was_sat |= !core.empty();
+            if (core.empty()) {
+                m_ex.assert_expr(m.mk_false());
+                m_fa.assert_expr(m.mk_false());
+                return;
+            }
+            expr_ref bound(m);
+            m_bound = m_mbp.maximize(core, mdl, m_objective, m_value, bound);
+            switch (m_bound) {
+            case mbp::unbounded:
+                m_ex.assert_expr(m.mk_false());
+                m_fa.assert_expr(m.mk_false());
+                break;
+            case mbp::strict:
+                m_ex.assert_expr(bound);
+                break;
+            case mbp::non_strict:
+                m_ex.assert_expr(bound);
+                break;
+            }
+        }
+
     };
 
+    lbool maximize(expr_ref_vector const& fmls, app* t, expr_ref& value, mbp::bound_type& bound, params_ref const& p) {
+        ast_manager& m = fmls.get_manager();
+        qsat qs(m, p, qsat_maximize);
+        return qs.maximize(fmls, t, value, bound);
+    }    
+};
+
+tactic * mk_qsat_tactic(ast_manager& m, params_ref const& p) {
+    return alloc(qe::qsat, m, p, qe::qsat_sat);
+}
+
+tactic * mk_qe2_tactic(ast_manager& m, params_ref const& p) {   
+    return alloc(qe::qsat, m, p, qe::qsat_qe);
+}
+
+tactic * mk_qe_rec_tactic(ast_manager& m, params_ref const& p) {   
+    return alloc(qe::qsat, m, p, qe::qsat_qe_rec);
+}
 
 
 
+
+#if 0
+
+    class min_max_opt {
+        struct imp;
+        imp* m_imp;
+    public:
+        min_max_opt(ast_manager& m);
+        ~min_max_opt();
+        void add(expr* e);
+        void add(expr_ref_vector const& fmls);
+        lbool check(svector<bool> const& is_max, app_ref_vector const& vars, app* t);
+    };
+
     struct min_max_opt::imp {
         ast_manager&    m;
         expr_ref_vector m_fmls;
@@ -1346,20 +1461,4 @@ namespace qe {
         return m_imp->check(is_max, vars, t);
     }
 
-
+#endif
-    
-};
-
-tactic * mk_qsat_tactic(ast_manager& m, params_ref const& p) {
-    return alloc(qe::qsat, m, p, false, true);
-}
-
-tactic * mk_qe2_tactic(ast_manager& m, params_ref const& p) {   
-    return alloc(qe::qsat, m, p, true, true);
-}
-
-tactic * mk_qe_rec_tactic(ast_manager& m, params_ref const& p) {   
-    return alloc(qe::qsat, m, p, true, false);
-}
-
-

src/qe/qsat.h

@@ -23,6 +23,7 @@ Revision History:
 
 #include "tactic.h"
 #include "filter_model_converter.h"
+#include "qe_mbp.h"
 
 namespace qe {
 
@@ -113,17 +114,7 @@ namespace qe {
         void collect_statistics(statistics& st) const;
     };
 
-    class min_max_opt {
+    lbool maximize(expr_ref_vector const& fmls, app* t, expr_ref& value, mbp::bound_type& bound, params_ref const& p);
-        struct imp;
-        imp* m_imp;
-    public:
-        min_max_opt(ast_manager& m);
-        ~min_max_opt();
-        void add(expr* e);
-        void add(expr_ref_vector const& fmls);
-        lbool check(svector<bool> const& is_max, app_ref_vector const& vars, app* t);
-    };
-
 
 }