factor out coi, use polynomial elaboration for nlsat solver

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

src/ast/sls/sls_arith_base.cpp

@@ -2620,8 +2620,10 @@ namespace sls {
                 display(out, ad) << "\n";
             }
         };
-        for (var_t v = 0; v < m_vars.size(); ++v) {
+        for (var_t v = 0; v < m_vars.size(); ++v) {            
             if (!eval_is_correct(v)) {
+//                if (m.rlimit().is_canceled())
+//                    return;
                 report_error(verbose_stream(), v);
                 TRACE(arith, report_error(tout, v));
                 UNREACHABLE();

src/math/lp/CMakeLists.txt

@@ -24,6 +24,7 @@ z3_add_component(lp
     monomial_bounds.cpp
     nex_creator.cpp
     nla_basics_lemmas.cpp
+    nla_coi.cpp
     nla_common.cpp
     nla_core.cpp
     nla_divisions.cpp

src/math/lp/nla_coi.cpp

@@ -0,0 +1,80 @@
+
+#include "math/lp/nla_core.h"
+#include "math/lp/nla_coi.h"
+
+namespace nla {
+
+  void coi::init() {
+      indexed_uint_set visited;
+      unsigned_vector todo;
+      vector<occurs> var2occurs;
+      m_term_set.reset();
+      m_mon_set.reset();
+      m_constraint_set.reset();
+      m_var_set.reset();
+      auto& lra = c.lra_solver();
+
+      for (auto ci : lra.constraints().indices()) {
+          auto const& c = lra.constraints()[ci];
+          if (c.is_auxiliary())
+              continue;
+          for (auto const& [coeff, v] : c.coeffs()) {
+              var2occurs.reserve(v + 1);
+              var2occurs[v].constraints.push_back(ci);
+          }
+      }
+
+      for (auto const& m : c.emons()) {
+          for (auto v : m.vars()) {
+              var2occurs.reserve(v + 1);
+              var2occurs[v].monics.push_back(m.var());
+          }
+      }
+
+      for (const auto *t :  lra.terms() ) {
+          for (auto const iv : *t) {
+              auto v = iv.j();
+              var2occurs.reserve(v + 1);
+              var2occurs[v].terms.push_back(t->j());
+          }
+      }
+
+      for (auto const& m : c.to_refine())
+          todo.push_back(m);
+
+      for (unsigned i = 0; i < todo.size(); ++i) {
+          auto v = todo[i];
+          if (visited.contains(v))
+              continue;
+          visited.insert(v);
+          m_var_set.insert(v);
+          var2occurs.reserve(v + 1);
+          for (auto ci : var2occurs[v].constraints) {
+              m_constraint_set.insert(ci);
+              auto const& c = lra.constraints()[ci];
+              for (auto const& [coeff, w] : c.coeffs())
+                  todo.push_back(w);
+          }
+          for (auto w : var2occurs[v].monics)
+              todo.push_back(w);
+
+          for (auto ti : var2occurs[v].terms) {
+              for (auto iv : lra.get_term(ti))
+                  todo.push_back(iv.j());
+              todo.push_back(ti);
+          }
+
+          if (lra.column_has_term(v)) {
+              m_term_set.insert(v);
+              for (auto kv : lra.get_term(v))
+                  todo.push_back(kv.j());
+          }            
+
+          if (c.is_monic_var(v)) {
+              m_mon_set.insert(v);
+              for (auto w : c.emons()[v])
+                  todo.push_back(w);
+          }
+      }
+  }
+}

src/math/lp/nla_coi.h

@@ -0,0 +1,30 @@
+
+#pragma once
+
+namespace nla {
+
+    class core;
+
+    class coi {
+        core& c;
+        indexed_uint_set m_mon_set, m_constraint_set;
+        indexed_uint_set m_term_set, m_var_set;
+
+        struct occurs {
+            unsigned_vector constraints;
+            unsigned_vector monics;
+            unsigned_vector terms;
+        };
+
+    public:
+        coi(core& c) : c(c) {}
+
+        void init();        
+
+        indexed_uint_set const& mons() const { return m_mon_set; }
+        indexed_uint_set const& constraints() const { return m_constraint_set; }
+        indexed_uint_set& terms() { return m_term_set; }
+        indexed_uint_set const &vars() { return m_var_set; }
+
+    };
+}

src/math/lp/nla_core.h

@@ -450,6 +450,12 @@ public:
     nla_throttle& throttle() { return m_throttle; }
     const nla_throttle& throttle() const { return m_throttle; }
 
+    lp::lar_solver& lra_solver() { return lra; }
+
+    indexed_uint_set const& to_refine() const {
+        return m_to_refine;
+    }
+
 };  // end of core
 
 struct pp_mon {

src/math/lp/nra_solver.cpp

@@ -9,6 +9,7 @@
 #include <fstream>
 #include "math/lp/lar_solver.h"
 #include "math/lp/nra_solver.h"
+#include "math/lp/nla_coi.h"
 #include "nlsat/nlsat_solver.h"
 #include "math/polynomial/polynomial.h"
 #include "math/polynomial/algebraic_numbers.h"
@@ -25,112 +26,150 @@ typedef nla::mon_eq mon_eq;
 typedef nla::variable_map_type variable_map_type;
 
 struct solver::imp {
+    struct model_bound {
+        lp::lpvar v;
+        rational r;
+        bool is_lower;
+    };
+
     lp::lar_solver&           lra;
     reslimit&                 m_limit;  
     params_ref                m_params; 
     u_map<polynomial::var>    m_lp2nl;  // map from lar_solver variables to nlsat::solver variables        
-    indexed_uint_set          m_term_set;
     scoped_ptr<nlsat::solver> m_nlsat;
     scoped_ptr<scoped_anum_vector>   m_values; // values provided by LRA solver
     scoped_ptr<scoped_anum> m_tmp1, m_tmp2;
+    nla::coi                  m_coi;
     nla::core&                m_nla_core;
+    unsigned                  m_max_constraint_index = 0;
+    vector<model_bound>       m_model_bounds;
     
     imp(lp::lar_solver& s, reslimit& lim, params_ref const& p, nla::core& nla_core): 
         lra(s), 
         m_limit(lim),
         m_params(p),
+        m_coi(nla_core),
         m_nla_core(nla_core) {}
 
     bool need_check() {
         return m_nla_core.m_to_refine.size() != 0;
     }
 
-    indexed_uint_set m_mon_set, m_constraint_set;
-
-    struct occurs {
-        unsigned_vector constraints;
-        unsigned_vector monics;
-        unsigned_vector terms;
-    };
-
-    void init_cone_of_influence() {
-        indexed_uint_set visited;
-        unsigned_vector todo;
-        vector<occurs> var2occurs;
-        m_term_set.reset();
-        m_mon_set.reset();
-        m_constraint_set.reset();
-
-        for (auto ci : lra.constraints().indices()) {
-            auto const& c = lra.constraints()[ci];
-            if (c.is_auxiliary())
-                continue;
-            for (auto const& [coeff, v] : c.coeffs()) {
-                var2occurs.reserve(v + 1);
-                var2occurs[v].constraints.push_back(ci);
-            }
-        }
-
-        for (auto const& m : m_nla_core.emons()) {
-            for (auto v : m.vars()) {
-                var2occurs.reserve(v + 1);
-                var2occurs[v].monics.push_back(m.var());
-            }
-        }
-
-        for (const auto *t :  lra.terms() ) {
-            for (auto const iv : *t) {
-                auto v = iv.j();
-                var2occurs.reserve(v + 1);
-                var2occurs[v].terms.push_back(t->j());
-            }
-        }
-
-        for (auto const& m : m_nla_core.m_to_refine)
-            todo.push_back(m);
-
-        for (unsigned i = 0; i < todo.size(); ++i) {
-            auto v = todo[i];
-            if (visited.contains(v))
-                continue;
-            visited.insert(v);
-            var2occurs.reserve(v + 1);
-            for (auto ci : var2occurs[v].constraints) {
-                m_constraint_set.insert(ci);
-                auto const& c = lra.constraints()[ci];
-                for (auto const& [coeff, w] : c.coeffs())
-                    todo.push_back(w);
-            }
-            for (auto w : var2occurs[v].monics)
-                todo.push_back(w);
-
-            for (auto ti : var2occurs[v].terms) {
-                for (auto iv : lra.get_term(ti))
-                    todo.push_back(iv.j());
-                todo.push_back(ti);
-            }
-
-            if (lra.column_has_term(v)) {
-                m_term_set.insert(v);
-                for (auto kv : lra.get_term(v))
-                    todo.push_back(kv.j());
-            }            
-
-            if (m_nla_core.is_monic_var(v)) {
-                m_mon_set.insert(v);
-                for (auto w : m_nla_core.emons()[v])
-                    todo.push_back(w);
-            }
-        }
-    }
-
     void reset() {
         m_values = nullptr;
         m_tmp1 = nullptr; m_tmp2 = nullptr;
         m_nlsat = alloc(nlsat::solver, m_limit, m_params, false);
         m_values = alloc(scoped_anum_vector, am());
-        m_term_set.reset();
         m_lp2nl.reset();
+        m_model_bounds.reset();
+        m_max_constraint_index = 0;
+    }
+
+    // Create polynomial definition for variable v used in setup_assignment_solver.
+    // Side-effects: updates m_vars2mon when v is a monic variable.
+    void mk_definition(unsigned v, polynomial_ref_vector &definitions, vector<rational>& denominators) {
+        auto &pm = m_nlsat->pm();
+        polynomial::polynomial_ref p(pm);
+        rational den(1);
+        if (m_nla_core.emons().is_monic_var(v)) {
+            auto const &m = m_nla_core.emons()[v];
+            for (auto v2 : m.vars()) {
+                polynomial_ref pw(definitions.get(v2), m_nlsat->pm());
+                if (!p)
+                    p = pw;
+                else
+                    p = p * pw;
+            }
+        }
+        else if (lra.column_has_term(v)) {
+            for (auto const &[w, coeff] : lra.get_term(v)) {
+                den = lcm(denominator(coeff), den);
+            }
+            for (auto const &[w, coeff] : lra.get_term(v)) {
+                auto coeff1 = den * coeff;
+                polynomial_ref pw(definitions.get(w), m_nlsat->pm());
+                if (!p)
+                    p = constant(coeff1) * pw;
+                else
+                    p = p + (constant(coeff1) * pw);
+            }
+        }
+        else {
+            p = pm.mk_polynomial(lp2nl(v));  // nlsat var index equals v (verified above when created)
+        }
+        definitions.push_back(p);
+        denominators.push_back(den);
+    }
+
+    void setup_solver_poly() {
+        m_coi.init();
+        auto &pm = m_nlsat->pm();
+        polynomial_ref_vector definitions(pm);
+        vector<rational> denominators;
+        for (unsigned v = 0; v < lra.number_of_vars(); ++v) {
+            if (m_coi.vars().contains(v)) {
+                auto j = m_nlsat->mk_var(lra.var_is_int(v));
+                m_lp2nl.insert(v, j);  // we don't really need this. It is going to be the identify map.
+                mk_definition(v, definitions, denominators);
+            }
+            else {
+                definitions.push_back(nullptr);
+                denominators.push_back(rational(0));
+            }
+        }
+
+        // we rely on that all information encoded into the tableau is present as a constraint.
+        for (auto ci : m_coi.constraints()) {
+            auto &c = lra.constraints()[ci];
+            auto &pm = m_nlsat->pm();
+            auto k = c.kind();
+            auto rhs = c.rhs();
+            auto lhs = c.coeffs();
+            rational den = denominator(rhs);
+            for (auto [coeff, v] : lhs)
+                den = lcm(lcm(den, denominator(coeff)), denominators[v]);
+            polynomial::polynomial_ref p(pm);
+            p = pm.mk_const(-den * rhs);
+
+            for (auto [coeff, v] : lhs) {
+                polynomial_ref poly(pm);
+                poly = definitions.get(v);
+                poly = poly * constant(den * coeff / denominators[v]);
+                p = p + poly;
+            }
+            auto lit = add_constraint(p, ci, k);
+        }
+        definitions.reset();
+    }
+
+    void setup_solver_terms() {
+        m_coi.init();
+        // add linear inequalities from lra_solver
+        for (auto ci : m_coi.constraints())
+            add_constraint(ci);
+
+        // add polynomial definitions.
+        for (auto const &m : m_coi.mons())
+            add_monic_eq(m_nla_core.emons()[m]);
+
+        // add term definitions.
+        for (unsigned i : m_coi.terms())
+            add_term(i);
+    }
+
+
+
+    polynomial::polynomial_ref sub(polynomial::polynomial *a, polynomial::polynomial *b) {
+        return polynomial_ref(m_nlsat->pm().sub(a, b), m_nlsat->pm());
+    }
+    polynomial::polynomial_ref mul(polynomial::polynomial *a, polynomial::polynomial *b) {
+        return polynomial_ref(m_nlsat->pm().mul(a, b), m_nlsat->pm());
+    }
+    polynomial::polynomial_ref var(lp::lpvar v) {
+        return polynomial_ref(m_nlsat->pm().mk_polynomial(lp2nl(v)), m_nlsat->pm());
+    }
+    polynomial::polynomial_ref constant(rational const& r) {
+        return polynomial_ref(m_nlsat->pm().mk_const(r), m_nlsat->pm());
     }
 
     /**
@@ -147,24 +186,14 @@ struct solver::imp {
     lbool check() {
         SASSERT(need_check());
         reset();
-        vector<nlsat::assumption, false> core;
-
-        init_cone_of_influence();
-        // add linear inequalities from lra_solver
-        for (auto ci : m_constraint_set)
-            add_constraint(ci);
+        vector<nlsat::assumption, false> core;        
         
-        // add polynomial definitions.
-        for (auto const& m : m_mon_set)
-            add_monic_eq(m_nla_core.emons()[m]);
+        smt_params_helper p(m_params);
 
-        // add term definitions.
-        for (unsigned i : m_term_set)
-            add_term(i);
+	setup_solver_poly();
 
         TRACE(nra, m_nlsat->display(tout));
 
-        smt_params_helper p(m_params);
         if (p.arith_nl_log()) {
             static unsigned id = 0;
             std::stringstream strm;
@@ -196,7 +225,8 @@ struct solver::imp {
             }
         }
         m_nlsat->collect_statistics(st);
-        TRACE(nra,
+        TRACE(nra, tout << "nra result " << r << "\n");
+        CTRACE(nra, false,
               m_nlsat->display(tout << r << "\n");
               display(tout);
               for (auto [j, x] : m_lp2nl) tout << "j" << j << " := x" << x << "\n";);
@@ -223,14 +253,26 @@ struct solver::imp {
         case l_false: {
             lp::explanation ex;
             m_nlsat->get_core(core);
+            nla::lemma_builder lemma(m_nla_core, __FUNCTION__);
             for (auto c : core) {
                 unsigned idx = static_cast<unsigned>(static_cast<imp*>(c) - this);
-                ex.push_back(idx);
-                TRACE(nra, lra.display_constraint(tout << "ex: " << idx << ": ", idx) << "\n";);
+                TRACE(nra, tout << "core index " << idx << " " << m_max_constraint_index << "\n");
+                if (idx <= m_max_constraint_index)
+                    ex.push_back(idx);
+                else {
+                    idx -= m_max_constraint_index;
+                    auto const& [v, bound, is_lower] = m_model_bounds[idx];
+                    TRACE(nra, tout << "bound violated for v" << v << (is_lower ? " >= " : " <= ") << bound << "\n");
+                    if (is_lower)
+                        lemma |= nla::ineq(v, lp::lconstraint_kind::LE, bound - 1);
+                    else
+                        lemma |= nla::ineq(v, lp::lconstraint_kind::GE, bound + 1);
+                }
             }
-            nla::lemma_builder lemma(m_nla_core, __FUNCTION__);
+
             lemma &= ex;
             m_nla_core.set_use_nra_model(true);
+            TRACE(nra, tout << lemma << "\n");
             break;
         }
         case l_undef:
@@ -272,13 +314,27 @@ struct solver::imp {
         coeffs.push_back(mpz(1));
         coeffs.push_back(mpz(-1));
         polynomial::polynomial_ref p(pm.mk_polynomial(2, coeffs.data(), mls),  pm);
-        polynomial::polynomial* ps[1] = { p };
-        bool even[1] = { false };
-        nlsat::literal lit = m_nlsat->mk_ineq_literal(nlsat::atom::kind::EQ, 1, ps, even);
+        auto lit = mk_literal(p.get(), lp::lconstraint_kind::EQ);
+        nlsat::assumption a = nullptr;
         m_nlsat->mk_clause(1, &lit, nullptr);
     }
 
+    nlsat::literal mk_literal(polynomial::polynomial *p, lp::lconstraint_kind k) {
+        polynomial::polynomial *ps[1] = { p };
+        bool is_even[1] = { false };
+        switch (k) {
+        case lp::lconstraint_kind::LE: return ~m_nlsat->mk_ineq_literal(nlsat::atom::kind::GT, 1, ps, is_even);
+        case lp::lconstraint_kind::GE: return ~m_nlsat->mk_ineq_literal(nlsat::atom::kind::LT, 1, ps, is_even);
+        case lp::lconstraint_kind::LT: return m_nlsat->mk_ineq_literal(nlsat::atom::kind::LT, 1, ps, is_even);
+        case lp::lconstraint_kind::GT: return m_nlsat->mk_ineq_literal(nlsat::atom::kind::GT, 1, ps, is_even);
+        case lp::lconstraint_kind::EQ: return m_nlsat->mk_ineq_literal(nlsat::atom::kind::EQ, 1, ps, is_even);
+        default: UNREACHABLE();  // unreachable
+        }
+        throw default_exception("uexpected operator");
+    }
+
     void add_constraint(unsigned idx) {
+        m_max_constraint_index = std::max(m_max_constraint_index, idx);
         auto& c = lra.constraints()[idx];
         auto& pm = m_nlsat->pm();
         auto k = c.kind();
@@ -297,30 +353,27 @@ struct solver::imp {
         }
         rhs *= den;
         polynomial::polynomial_ref p(pm.mk_linear(sz, coeffs.data(), vars.data(), -rhs), pm);
-        polynomial::polynomial* ps[1] = { p };
-        bool is_even[1] = { false };
+        nlsat::literal lit = mk_literal(p.get(), k);
+        nlsat::assumption a = this + idx;
+        m_nlsat->mk_clause(1, &lit, a);
+    }
+
+    nlsat::literal add_constraint(polynomial::polynomial *p, unsigned idx, lp::lconstraint_kind k) {
+        m_max_constraint_index = std::max(m_max_constraint_index, idx);
+        polynomial::polynomial *ps[1] = {p};
+        bool is_even[1] = {false};
         nlsat::literal lit;
         nlsat::assumption a = this + idx;
         switch (k) {
-        case lp::lconstraint_kind::LE:
-            lit = ~m_nlsat->mk_ineq_literal(nlsat::atom::kind::GT, 1, ps, is_even);
-            break;
-        case lp::lconstraint_kind::GE:
-            lit = ~m_nlsat->mk_ineq_literal(nlsat::atom::kind::LT, 1, ps, is_even);
-            break;
-        case lp::lconstraint_kind::LT:
-            lit = m_nlsat->mk_ineq_literal(nlsat::atom::kind::LT, 1, ps, is_even);
-            break;
-        case lp::lconstraint_kind::GT:
-            lit = m_nlsat->mk_ineq_literal(nlsat::atom::kind::GT, 1, ps, is_even);
-            break;
-        case lp::lconstraint_kind::EQ:
-            lit = m_nlsat->mk_ineq_literal(nlsat::atom::kind::EQ, 1, ps, is_even);                
-            break;
-        default:
-            UNREACHABLE(); // unreachable
+        case lp::lconstraint_kind::LE: lit = ~m_nlsat->mk_ineq_literal(nlsat::atom::kind::GT, 1, ps, is_even); break;
+        case lp::lconstraint_kind::GE: lit = ~m_nlsat->mk_ineq_literal(nlsat::atom::kind::LT, 1, ps, is_even); break;
+        case lp::lconstraint_kind::LT: lit = m_nlsat->mk_ineq_literal(nlsat::atom::kind::LT, 1, ps, is_even); break;
+        case lp::lconstraint_kind::GT: lit = m_nlsat->mk_ineq_literal(nlsat::atom::kind::GT, 1, ps, is_even); break;
+        case lp::lconstraint_kind::EQ: lit = m_nlsat->mk_ineq_literal(nlsat::atom::kind::EQ, 1, ps, is_even); break;
+        default: UNREACHABLE();  // unreachable
         }
         m_nlsat->mk_clause(1, &lit, a);
+        return lit;
     }
 
     bool check_monic(mon_eq const& m) {
@@ -370,7 +423,7 @@ struct solver::imp {
         for (auto const& m : m_nla_core.emons())
             if (any_of(m.vars(), [&](lp::lpvar v) { return m_lp2nl.contains(v); }))
                 add_monic_eq_bound(m);
-        for (unsigned i : m_term_set)
+        for (unsigned i : m_coi.terms())
             add_term(i);
         for (auto const& [v, w] : m_lp2nl) {
             if (lra.column_has_lower_bound(v))
@@ -418,6 +471,7 @@ struct solver::imp {
                 ex.push_back(ci);
             nla::lemma_builder lemma(m_nla_core, __FUNCTION__);
             lemma &= ex;
+            TRACE(nra, tout << lemma << "\n");
             break;
         }
         case l_undef:
@@ -554,8 +608,8 @@ struct solver::imp {
         if (!m_lp2nl.find(v, r)) {
             r = m_nlsat->mk_var(is_int(v));
             m_lp2nl.insert(v, r);
-            if (!m_term_set.contains(v) && lra.column_has_term(v)) {
-                m_term_set.insert(v);
+            if (!m_coi.terms().contains(v) && lra.column_has_term(v)) {
+                m_coi.terms().insert(v);
             }
         }
         return r;
@@ -586,20 +640,55 @@ struct solver::imp {
         m_nlsat->mk_clause(1, &lit, nullptr);
     }
 
-    nlsat::anum const& value(lp::lpvar v)  {
-        polynomial::var pv;
-        if (m_lp2nl.find(v, pv))
-            return m_nlsat->value(pv);
-        else {
-            for (unsigned w = m_values->size(); w <= v; ++w) {
-                scoped_anum a(am());
-                am().set(a, m_nla_core.val(w).to_mpq());
+    nlsat::anum const &value(lp::lpvar v) {
+        init_values(v + 1);
+        return (*m_values)[v];
+    }
+
+    void init_values(unsigned sz) {
+        if (m_values->size() >= sz)
+            return;
+        unsigned w;
+        scoped_anum a(am());
+        for (unsigned v = m_values->size(); v < sz; ++v) {
+            if (m_nla_core.emons().is_monic_var(v)) {
+                am().set(a, 1);
+                auto &m = m_nla_core.emon(v);
+                for (auto x : m.vars()) 
+                    am().mul(a, (*m_values)[x], a);
                 m_values->push_back(a);
-            }           
-            return (*m_values)[v];
+            }      
+            else if (lra.column_has_term(v)) {
+                scoped_anum b(am());
+                am().set(a, 0);
+                for (auto const &[w, coeff] : lra.get_term(v)) {
+                    am().set(b, coeff.to_mpq());
+                    am().mul(b, (*m_values)[w], b);
+                    am().add(a, b, a);
+                }
+                m_values->push_back(a);
+            }
+            else if (m_lp2nl.find(v, w)) {
+                m_values->push_back(m_nlsat->value(w));
+            }
+            else {
+                am().set(a, m_nla_core.val(v).to_mpq());
+                m_values->push_back(a);
+            }
         }
     }
 
+
+    void set_value(lp::lpvar v, rational const& value) {
+        if (!m_values)
+            m_values = alloc(scoped_anum_vector, am());
+        scoped_anum a(am());
+        am().set(a, value.to_mpq());
+        while (m_values->size() <= v)
+            m_values->push_back(a);
+        am().set((*m_values)[v], a);
+    }
+
     nlsat::anum_manager& am() {
         return m_nlsat->am();
     }
@@ -680,4 +769,8 @@ void solver::updt_params(params_ref& p) {
     m_imp->updt_params(p);
 }
 
+void solver::set_value(lp::lpvar v, rational const& value) {
+    m_imp->set_value(v, value);
+}
+
 }

src/math/lp/nra_solver.h

@@ -59,6 +59,8 @@ namespace nra {
 
         nlsat::anum_manager& am();        
 
+        void set_value(lp::lpvar v, rational const &value);
+
         scoped_anum& tmp1();
 
         scoped_anum& tmp2();