outline for adding monomials

src/math/lp/explanation.h

@@ -23,10 +23,10 @@ Revision History:
 #include "util/hashtable.h"
 namespace lp {
 class explanation {
-    typedef vector<std::pair<unsigned, mpq>> pair_vec;
+    typedef vector<std::pair<constraint_index, mpq>> pair_vec;
     typedef    hashtable<unsigned, u_hash, u_eq>  ci_set;
     // Only one of the fields below is used. The first call adding an entry decides which one it is.
-    vector<std::pair<constraint_index, mpq>> m_vector;
+    pair_vec m_vector;
     ci_set   m_set;   
 public:
     explanation() = default;

src/math/lp/lp_types.h

@@ -19,8 +19,8 @@ Revision History:
 --*/
 #pragma once
 
-#include <sstream>
 #include <limits.h>
+#include <sstream>
 #include "util/debug.h"
 #include "util/dependency.h"
 
@@ -30,12 +30,17 @@ namespace nla {
 
 namespace lp {
 
-typedef unsigned constraint_index;
+    typedef unsigned constraint_index;
-typedef unsigned row_index;
+    typedef unsigned row_index;
-enum lconstraint_kind { LE = -2, LT = -1 , GE = 2, GT = 1, EQ = 0, NE = 3 };
+    enum lconstraint_kind {
-typedef unsigned lpvar;
+        LE = -2,
-const lpvar null_lpvar = UINT_MAX;
+        LT = -1,
-const constraint_index null_ci = UINT_MAX;
+        GE = 2,
-}
+        GT = 1,
-
+        EQ = 0,
-
+        NE = 3
+    };
+    typedef unsigned lpvar;
+    const lpvar null_lpvar = UINT_MAX;
+    const constraint_index null_ci = UINT_MAX;
+}  // namespace lp

src/math/lp/mon_eq.cpp

@@ -9,28 +9,20 @@ namespace nla {
 template <typename T>
 bool check_assignment(T const& m, variable_map_type & vars) {
     rational r1 = vars[m.var()];
-    if (r1.is_zero()) {
+    if (r1.is_zero()) 
-        for (auto w : m.vars()) {
+        return any_of(m.vars(), [&](auto w) { return vars[w].is_zero(); });    
-            if (vars[w].is_zero())
-                return true;
-        }
-        return false;
-    }
     rational r2(1);
-    for (auto w : m.vars()) {
+    for (auto w : m.vars()) 
         r2 *= vars[w];    
-    }
     return r1 == r2;
 }
+
 template <typename K>
 bool check_assignments(const K & monomials,
                        const lp::lar_solver& s,
                        variable_map_type & vars) {
     s.get_model(vars);
-    for (auto const& m : monomials) {
+    return all_of(monomials, [&](auto const& m) { return check_assignment(m, vars); });
-        if (!check_assignment(m, vars)) return false;
-    }
-    return true;
 }
 
 template bool check_assignments<vector<mon_eq>>(const vector<mon_eq>&,

src/math/lp/nla_basics_lemmas.cpp

@@ -245,11 +245,7 @@ bool basics::basic_lemma(bool derived) {
     if (derived) 
         return false;
     const auto& mon_inds_to_ref = c().m_to_refine;
-    TRACE(nla_solver, tout << "mon_inds_to_ref = "; print_vector(mon_inds_to_ref, tout) << "\n";);
+    for (auto v : mon_inds_to_ref) {
-    unsigned start = c().random();
-    unsigned sz = mon_inds_to_ref.size();
-    for (unsigned j = 0; j < sz; ++j) {
-        lpvar v = mon_inds_to_ref[(j + start) % mon_inds_to_ref.size()];
         const monic& r = c().emons()[v];
         SASSERT (!c().check_monic(c().emons()[v]));
         basic_lemma_for_mon(r, derived);

src/math/lp/nla_core.cpp

@@ -35,6 +35,7 @@ core::core(lp::lar_solver& s, params_ref const& p, reslimit & lim) :
     m_divisions(*this),
     m_intervals(this, lim),
     m_monomial_bounds(this),
+    m_mul_saturate(this),
     m_horner(this),
     m_grobner(this),
     m_emons(m_evars),
@@ -1332,7 +1333,6 @@ lbool core::check() {
             return l_false;
     }
     
-    
     if (no_effect() && should_run_bounded_nlsat()) 
         ret = bounded_nlsat();
                 
@@ -1348,6 +1348,9 @@ lbool core::check() {
     if (no_effect())
         m_order.order_lemma();
 
+    if (false && no_effect())
+        ret = m_mul_saturate.saturate();
+
     if (no_effect()) {
         unsigned num_calls = lp_settings().stats().m_nla_calls;
         if (!conflict_found() && params().arith_nl_nra() && num_calls % 50 == 0 && num_calls > 500)

src/math/lp/nla_core.h

@@ -21,6 +21,7 @@
 #include "math/lp/nla_grobner.h"
 #include "math/lp/nla_powers.h"
 #include "math/lp/nla_divisions.h"
+#include "math/lp/nla_mul_saturate.h"
 #include "math/lp/emonics.h"
 #include "math/lp/nex.h"
 #include "math/lp/horner.h"
@@ -90,6 +91,7 @@ class core {
     divisions                m_divisions;
     intervals                m_intervals; 
     monomial_bounds          m_monomial_bounds;
+    mul_saturate             m_mul_saturate;
     unsigned                 m_conflicts;
     bool                     m_check_feasible = false;
     horner                   m_horner;
@@ -237,11 +239,16 @@ public:
     std::ostream & display_row(std::ostream& out, lp::row_strip<lp::mpq> const& row) const;
     std::ostream & display(std::ostream& out);
     std::ostream& display_smt(std::ostream& out);
+    std::ostream& display_coeff(std::ostream& out, bool first, lp::mpq const& p) const;
+    std::ostream& display_constraint(std::ostream& out, lp::constraint_index ci) const;
+    std::ostream& display_constraint(std::ostream& out, lp::lar_term const& lhs, lp::lconstraint_kind k, lp::mpq const& rhs) const;
+    std::ostream& display_constraint(std::ostream& out, vector<std::pair<rational, lpvar>> const& lhs, lp::lconstraint_kind k, lp::mpq const& rhs) const;
     std::ostream & print_ineq(const ineq & in, std::ostream & out) const;
     std::ostream & print_var(lpvar j, std::ostream & out) const;
     std::ostream & print_monics(std::ostream & out) const;    
     std::ostream & print_ineqs(const lemma& l, std::ostream & out) const;    
     std::ostream & print_factorization(const factorization& f, std::ostream& out) const;
+    
     template <typename T>
     std::ostream& print_product(const T & m, std::ostream& out) const;    
     template <typename T>

src/math/lp/nla_mul_saturate.cpp

@@ -10,7 +10,11 @@
   Check if the system with new constraints is LP feasible.
   If it is not, then produce a lemma that explains the infeasibility.
 
-  The lemma is in terms of the original constraints and bounds.
+  Strategy 1: The lemma is in terms of the original constraints and bounds.
+  Strategy 2: Attempt to eliminate new monomials from the lemma by relying on Farkas multipliers.
+              If it succeeds to eliminate new monomials we have a lemma that is a linear 
+              combination of existing variables. 
+  Strategy 3: The lemma uses the new constraints.
 
   --*/
 
@@ -21,78 +25,144 @@
 namespace nla {
 
     mul_saturate::mul_saturate(core* core) : 
-        common(core), 
+        common(core) {}
-        lra(m_core.lra) {}
 
     lbool mul_saturate::saturate() {
-        lra.push();
+        lp::explanation ex;
-        for (auto j : c().m_to_refine) {
+        init_solver();
-            auto& m = c().emons()[j];
+        add_multiply_constraints();
-            for (auto& con : lra.constraints().active()) {                
+        lbool r = solve(ex);
-                for (auto v : m.vars()) {
+        if (r == l_false)
-                    for (auto [coeff, u] : con.coeffs()) {
+            add_lemma(ex);
-                        if (u == v) 
+        return r;
-                            // multiply by remaining vars
-                            multiply_constraint(con, m, v);
-                            // add new constraint
     }
-                }
-            }
-        }       
-        // record new monomials that are created and recursively down-saturate with respect to these.
 
-        auto st = lra.solve();
+    void mul_saturate::init_solver() {
+        local_solver = alloc(lp::lar_solver);
+    }
+
+    void mul_saturate::add_lemma(lp::explanation const& ex1) {
+        lp::explanation ex2;
+        for (auto p : ex1) {
+            lp::constraint_index src_ci;
+            if (m_new_mul_constraints.find(p.ci(), src_ci))
+                ex2.add_pair(src_ci, mpq(1));
+            else
+                ex2.add_pair(p.ci(), p.coeff());
+        }
+        lemma_builder new_lemma(c(), "stellensatz");
+        new_lemma &= ex2;
+        for (auto [v, sign] : m_var_signs) {
+            if (sign)
+                new_lemma.explain_existing_upper_bound(v);
+            else
+                new_lemma.explain_existing_lower_bound(v);
+        }
+        IF_VERBOSE(1, verbose_stream() << "unsat \n" << new_lemma << "\n");
+    }
+
+    lbool mul_saturate::solve(lp::explanation& ex) {
+        for (auto const& [new_ci, old_ci] : m_new_mul_constraints)
+            local_solver->activate(new_ci);
+        auto st = local_solver->solve();
         lbool r = l_undef;
         if (st == lp::lp_status::INFEASIBLE) {
-            // now we need to filter new constraints into bounds and old constraints.
+            local_solver->get_infeasibility_explanation(ex);
-
+            IF_VERBOSE(0, c().print_explanation(ex, verbose_stream()) << "\n";);
             r = l_false;
         }
         if (st == lp::lp_status::OPTIMAL || st == lp::lp_status::FEASIBLE) {
             // TODO: check model just in case it got lucky.
+            IF_VERBOSE(1, verbose_stream() << "saturation is LP feasible, maybe it is a model of the NLA problem\n");
         }
-        lra.pop(1);
+        IF_VERBOSE(0, local_solver->display(verbose_stream()); c().display(verbose_stream()));
         return r;
     }
 
+    // record new monomials that are created and recursively down-saturate with respect to these.
+    void mul_saturate::add_multiply_constraints() {
+        m_new_mul_constraints.reset();
+        m_seen_vars.reset();
+        m_var_signs.reset();
+        for (auto j : c().m_to_refine) {
+            for (auto con_id : local_solver->constraints().indices()) {
+                unsigned num_vars = c().emon(j).vars().size();
+                for (unsigned i = 0; i < num_vars; ++i) {
+                    auto v = c().emon(j).vars()[i];
+                    for (auto [coeff, u] : local_solver->constraints()[con_id].coeffs())
+                        if (u == v)
+                            add_multiply_constraint(con_id, j, v);
+                }
+            }
+        }
+    }
+
     // multiply by remaining vars
-    void mul_saturate::multiply_constraint(lp::lar_base_constraint const& con, monic const& m, lpvar x) {
+    void mul_saturate::add_multiply_constraint(lp::constraint_index old_ci, lp::lpvar mi, lpvar x) {
+        lp::lar_base_constraint const& con = local_solver->constraints()[old_ci];
         auto const& lhs = con.coeffs();
         auto const& rhs = con.rhs();
         auto k = con.kind();        
+        if (k == lp::lconstraint_kind::NE || k == lp::lconstraint_kind::EQ)
+            return;  // not supported
         auto sign = false;
         svector<lpvar> vars;
         bool first = true;
-        for (auto v : m.vars()) {
+        for (auto v : c().emon(mi).vars()) {
             if (v != x || !first) 
                 vars.push_back(v);            
             else
                 first = false;
         }
-        vector<std::pair<rational, lpvar>> new_lhs;
+        // compute sign of vars
+        for (auto v : vars) {
+            if (m_seen_vars.contains(v))
+                continue;
+            // retrieve bounds of v
+            // if v has non-negative lower bound add as positive
+            // if v has non-positive upper bound add as negative
+            // otherwise, fail
+            if (local_solver->column_has_lower_bound(v) && !local_solver->get_lower_bound(v).is_neg()) {
+                m_var_signs.push_back({v, false});
+                m_seen_vars.insert(v);
+            } 
+            else if (local_solver->column_has_upper_bound(v) && !local_solver->get_upper_bound(v).is_pos()) {
+                m_var_signs.push_back({v, true});
+                m_seen_vars.insert(v);
+                sign = !sign;
+            }
+            else 
+                return;
+        }
+        lp::lar_term new_lhs;        
         rational new_rhs(rhs);
         for (auto [coeff, v] : lhs) {
+            #if 0
             vars.push_back(v);
-
+            lpvar new_monic_var = c().m_add_monomial(vars);
-            auto new_m = c().emons().find_canonical(vars);
+            auto const& new_m = c().emons()[new_monic_var];            
-            if (!new_m) {
+            verbose_stream() << vars << " v " << new_m.var() << " coeff " << coeff << "\n";
-                bool is_int = lra.var_is_int(x);  // assume all vars in monic have the same type, can be changed for MIP
+            new_lhs.add_monomial(coeff, new_m.var());
-                lpvar new_monic_var = 0; // lra.add_var(is_int);
-                c().emons().add(new_monic_var, vars);
-                new_m = c().emons().find_canonical(vars);
-                SASSERT(new_m);
-            }
-            new_lhs.push_back({coeff, new_m->var()});
             vars.pop_back();
+            #endif
         }
         if (rhs != 0) {
-            new_lhs.push_back({-rhs, m.var()});
+            if (vars.size() == 1) {
+                new_lhs.add_monomial(-rhs, vars[0]);
+                verbose_stream() << "rhs mul " << -rhs << " * j" << vars[0] << "\n";
+            } 
+            else {
+#if 0
+                lpvar new_monic_var = c().m_add_monomial(vars);
+                auto const& new_m = c().emons()[new_monic_var];
+                verbose_stream() << vars << " v " << new_m.var() << " coeff " << coeff << "\n";
+                new_lhs.add_monomial(-rhs, new_m.var());
+                verbose_stream() << "rhs mul " << -rhs << " * j" << new_m.var() << "\n";
+#endif
+            }
             new_rhs = 0;
         }
-        // compute sign of vars
+
-        for (auto v : vars) 
-            if (c().val(v).is_neg())
-                sign = !sign;
         if (sign) {
             switch (k) {
             case lp::lconstraint_kind::LE: k = lp::lconstraint_kind::GE; break;
@@ -102,6 +172,10 @@ namespace nla {
             default: break;
             }
         }   
-        // instead of adding a constraint here, add row to tableau based on the new_lhs, new_rhs, k.
+        c().display_constraint(verbose_stream(), old_ci) << " -> ";
+        c().display_constraint(verbose_stream(), new_lhs, k, new_rhs) << "\n";
+        // TODO: 
+        // auto new_ci = lra.m_add_constraint(new_lhs, k, new_rhs);
+        // m_new_mul_constraints.insert(new_ci, old_ci);        
     }
 }

src/math/lp/nla_mul_saturate.h

@@ -9,8 +9,16 @@ namespace nla {
     class core;
     class lar_solver;
     class mul_saturate : common {
-        lp::lar_solver& lra;
+        // source of multiplication constraint
-        void multiply_constraint(lp::lar_base_constraint const& c, monic const& m, lpvar x);
+        u_map<lp::constraint_index> m_new_mul_constraints;
+        svector<std::pair<lpvar, bool>> m_var_signs;
+        tracked_uint_set m_seen_vars;
+        scoped_ptr<lp::lar_solver> local_solver;
+        void init_solver();
+        void add_multiply_constraints();
+        void add_multiply_constraint(lp::constraint_index con_id, lp::lpvar mi, lpvar x);
+        lbool solve(lp::explanation& ex);
+        void add_lemma(lp::explanation const& ex1);
     public:
         mul_saturate(core* core);
 

src/math/lp/nla_pp.cpp

@@ -50,7 +50,8 @@ std::ostream& core::print_factor(const factor& f, std::ostream& out) const {
         out << "- ";
     if (f.is_var()) {
         out << "VAR,  " << pp(f.var());
-    } else {
+    } 
+    else {
         out << "MON, v" << m_emons[f.var()] << " = ";
         print_product(m_emons[f.var()].rvars(), out);
     }
@@ -61,7 +62,8 @@ std::ostream& core::print_factor(const factor& f, std::ostream& out) const {
 std::ostream& core::print_factor_with_vars(const factor& f, std::ostream& out) const {
     if (f.is_var()) {
         out << pp(f.var());
-    } else {
+    } 
+    else {
         out << " MON = " << pp_mon_with_vars(*this, m_emons[f.var()]);
     }
     return out;
@@ -133,9 +135,8 @@ std::ostream& core::print_var(lpvar j, std::ostream& out) const {
     lra.print_column_info(j, out);
     signed_var jr = m_evars.find(j);
     out << "root=";
-    if (jr.sign()) {
+    if (jr.sign()) 
         out << "-";    
-    }
 
     out << lra.get_variable_name(jr.var()) << "\n";
     return out;
@@ -245,10 +246,9 @@ std::string core::var_str(lpvar j) const {
     return result;
 }
 
-std::ostream& core::display_row(std::ostream& out, lp::row_strip<lp::mpq> const& row) const {
+std::ostream& core::display_coeff(std::ostream& out, bool first, lp::mpq const& p) const {
-    auto display_coeff = [&](bool first, lp::mpq const& p) {
     if (first && p == 1)
-            return;
+        return out;
     if (first && p > 0)
         out << p;
     else if (p == 1)
@@ -261,7 +261,10 @@ std::ostream& core::display_row(std::ostream& out, lp::row_strip<lp::mpq> const&
         out << p << " * ";
     else
         out << " - " << -p << " * ";
-    };
+    return out;
+}
+
+std::ostream& core::display_row(std::ostream& out, lp::row_strip<lp::mpq> const& row) const {
     auto display_var = [&](bool first, lp::mpq p, lp::lpvar v) {
         if (is_monic_var(v)) {
             for (auto w : m_emons[v].vars())
@@ -270,7 +273,7 @@ std::ostream& core::display_row(std::ostream& out, lp::row_strip<lp::mpq> const&
         else 
             p *= m_evars.find(v).rsign();
         
-        display_coeff(first, p);
+        display_coeff(out, first, p);
         if (is_monic_var(v)) {
             bool first = true;
             for (auto w : m_emons[v].vars())
@@ -356,6 +359,28 @@ std::ostream& core::display_declarations_smt(std::ostream& out) const {
     return out;
 }
 
+std::ostream& core::display_constraint(std::ostream& out, lp::constraint_index ci) const {
+    auto const& c = lra.constraints()[ci];
+    return display_constraint(out, c.coeffs(), c.kind(), c.rhs());
+}
+
+std::ostream& core::display_constraint(std::ostream& out, lp::lar_term const& lhs, lp::lconstraint_kind k, lp::mpq const& rhs) const {
+    return display_constraint(out, lhs.coeffs_as_vector(), k, rhs);
+}
+
+std::ostream& core::display_constraint(std::ostream& out, vector<std::pair<rational, lpvar>> const & lhs, lp::lconstraint_kind k, lp::mpq const& rhs) const {
+    bool first = true;
+    for (auto [coeff, v] : lhs) {
+        display_coeff(out, first, coeff);
+        first = false;
+        if (is_monic_var(v))
+            print_product(m_emons[v], out);      
+        else
+            out << "j" << v;
+    }
+    return out << " " << k << " " << rhs;    
+}
+
 std::ostream& core::display_constraint_smt(std::ostream& out, unsigned id, lp::lar_base_constraint const& c) const {
     auto k = c.kind();
     auto rhs = c.rhs();

src/smt/smt_context.h

@@ -70,8 +70,8 @@ namespace smt {
 
     struct enode_pp {
         context const& ctx;
-        enode*   n;
+        enode const*   n;
-        enode_pp(enode* n, context const& ctx): ctx(ctx), n(n) {}
+        enode_pp(enode const* n, context const& ctx): ctx(ctx), n(n) {}
     };
 
     struct replay_unit {
@@ -1414,7 +1414,7 @@ namespace smt {
 
         void display_asserted_formulas(std::ostream & out) const;
 
-        enode_pp pp(enode* n) { return enode_pp(n, *this); }
+        enode_pp pp(enode const* n) { return enode_pp(n, *this); }
 
         std::ostream& display_literal(std::ostream & out, literal l) const;
 

src/smt/smt_context_pp.cpp

@@ -692,7 +692,7 @@ namespace smt {
 
     std::ostream& operator<<(std::ostream& out, enode_pp const& p) {
         ast_manager& m = p.ctx.get_manager();
-        enode* n = p.n;
+        enode const* n = p.n;
         return out << n->get_owner_id() << ": " << mk_bounded_pp(n->get_expr(), m);
     }
 

src/smt/smt_theory.cpp

@@ -59,6 +59,11 @@ namespace smt {
         }
     }
 
+    bool theory::is_attached_to_var(enode const *n) const {
+        theory_var v = n->get_th_var(get_id());
+        return v != null_theory_var && get_enode(v) == n;
+    }
+
     void theory::display_var2enode(std::ostream & out) const {
         unsigned sz = m_var2enode.size();
         for (unsigned v = 0; v < sz; v++) {

src/smt/smt_theory.h

@@ -97,10 +97,7 @@ namespace smt {
            but it may be inherited from another enode n' during an
            equivalence class merge. That is, get_enode(v) != n.
         */
-        bool is_attached_to_var(enode const * n) const {
+        bool is_attached_to_var(enode const *n) const;
-            theory_var v = n->get_th_var(get_id());
-            return v != null_theory_var && get_enode(v) == n;
-        }
 
         struct scoped_trace_stream {
             ast_manager& m;

src/smt/theory_lra.cpp

@@ -890,7 +890,7 @@ public:
             mk_is_int_axiom(n);        
     }
 
-    bool internalize_atom(app * atom, bool gate_ctx) {
+    void internalize_atom(app * atom) {
         TRACE(arith_internalize, tout << bpp(atom) << "\n";);
         SASSERT(!ctx().b_internalized(atom));
         expr* n1, *n2;
@@ -918,12 +918,12 @@ public:
         }
         else if (a.is_is_int(atom)) {
             internalize_is_int(atom);
-            return true;
+            return;
         }
         else {
             TRACE(arith, tout << "Could not internalize " << mk_pp(atom, m) << "\n";);
             found_unsupported(atom);
-            return true;
+            return;
         }
 
         if (is_int(v) && !r.is_int()) 
@@ -936,7 +936,6 @@ public:
         m_bool_var2bound.insert(bv, b);
         mk_bound_axioms(*b);
         TRACE(arith_internalize, tout << "Internalized " << bv << ": " << bpp(atom) << "\n";);
-        return true;
     }
         
     bool internalize_term(app * term) {
@@ -1730,15 +1729,10 @@ public:
             return atom;
         }
     }
-    // create a bound atom representing term >= k is lower_bound is true, and term <= k if it is false
-    expr_ref mk_bound(lp::lar_term const& term, rational const& k, bool lower_bound) {
-        rational offset;
-        expr_ref t(m);
-        return mk_bound(term, k, lower_bound, offset, t);
-    }
 
-    expr_ref mk_bound(lp::lar_term const& term, rational const& k, bool lower_bound, rational& offset, expr_ref& t) {
+    expr_ref mk_bound(lp::lar_term const& term, rational const& k, bool lower_bound) {
-        offset = k;
+        rational offset = k;
+        expr_ref t(m);
         u_map<rational> coeffs;
         term2coeffs(term, coeffs);
         bool is_int = true;
@@ -1925,9 +1919,7 @@ public:
             TRACE(arith, tout << "branch\n";);
             bool u = m_lia->is_upper();
             auto const & k = m_lia->offset();
-            rational offset;
+            expr_ref b = mk_bound(m_lia->get_term(), k, !u);
-            expr_ref t(m);
-            expr_ref b = mk_bound(m_lia->get_term(), k, !u, offset, t);
             if (m.has_trace_stream()) {
                 app_ref body(m);
                 body = m.mk_or(b, m.mk_not(b));
@@ -4192,16 +4184,13 @@ public:
     unsigned init_range() const { return 5; }
     unsigned max_range() const { return 20; }
     
-
     void setup() {
         m_bounded_range_lit = null_literal;
         m_bound_terms.reset();
         m_bound_predicate = nullptr;
     }
 
-
     void validate_model(proto_model& mdl) {
-
         rational r1, r2;
         expr_ref res(m);
         if (!m_model_is_initialized)
@@ -4240,7 +4229,8 @@ void theory_lra::init() {
     m_imp->init();
 }    
 bool theory_lra::internalize_atom(app * atom, bool gate_ctx) {
-    return m_imp->internalize_atom(atom, gate_ctx);
+    m_imp->internalize_atom(atom);
+    return true;
 }
 bool theory_lra::internalize_term(app * term) {
     return m_imp->internalize_term(term);