Dev (#51)

* initial skeletons for nra solver

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

* initial skeletons for nra solver

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

* adding more nlsat

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

* nlsat integration

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

* adding constraints

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

* adding nra solver

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

* add missing initialization

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

* adding nra solver

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

contrib/cmake/src/CMakeLists.txt

@@ -35,10 +35,10 @@ endforeach()
 # raised if you try to declare a component is dependent on another component
 # that has not yet been declared.
 add_subdirectory(util)
-add_subdirectory(util/lp)
 add_subdirectory(math/polynomial)
 add_subdirectory(sat)
 add_subdirectory(nlsat)
+add_subdirectory(util/lp)
 add_subdirectory(math/hilbert)
 add_subdirectory(math/simplex)
 add_subdirectory(math/automata)

contrib/cmake/src/smt/CMakeLists.txt

@@ -70,6 +70,7 @@ z3_add_component(smt
     euclid
     fpa
     grobner
+    nlsat
     lp
     macros
     normal_forms

contrib/cmake/src/test/lp/CMakeLists.txt

@@ -1,4 +1,4 @@
-add_executable(lp_tst lp_main.cpp lp.cpp $<TARGET_OBJECTS:util>   $<TARGET_OBJECTS:lp>)
+add_executable(lp_tst lp_main.cpp lp.cpp $<TARGET_OBJECTS:util> $<TARGET_OBJECTS:polynomial> $<TARGET_OBJECTS:nlsat>  $<TARGET_OBJECTS:lp> )
 target_compile_definitions(lp_tst PRIVATE ${Z3_COMPONENT_CXX_DEFINES})
 target_compile_options(lp_tst PRIVATE ${Z3_COMPONENT_CXX_FLAGS})
 target_include_directories(lp_tst PRIVATE ${Z3_COMPONENT_EXTRA_INCLUDE_DIRS})

contrib/cmake/src/util/lp/CMakeLists.txt

@@ -30,6 +30,8 @@ z3_add_component(lp
     random_updater_instances.cpp      
   COMPONENT_DEPENDENCIES
     util
+    polynomial
+    nlsat
   PYG_FILES
     lp_params.pyg
 )

src/nlsat/nlsat_solver.h

@@ -21,10 +21,10 @@ Revision History:
 #ifndef NLSAT_SOLVER_H_
 #define NLSAT_SOLVER_H_
 
-#include"nlsat_types.h"
-#include"params.h"
-#include"statistics.h"
-#include"rlimit.h"
+#include"nlsat/nlsat_types.h"
+#include"util/params.h"
+#include"util/statistics.h"
+#include"util/rlimit.h"
 
 namespace nlsat {
 

src/nlsat/nlsat_types.h

@@ -19,10 +19,10 @@ Revision History:
 #ifndef NLSAT_TYPES_H_
 #define NLSAT_TYPES_H_
 
-#include"polynomial.h"
-#include"buffer.h"
-#include"sat_types.h"
-#include"z3_exception.h"
+#include"math/polynomial/polynomial.h"
+#include"util/buffer.h"
+#include"sat/sat_types.h"
+#include"util/z3_exception.h"
 
 namespace algebraic_numbers {
     class anum;

src/smt/theory_lra.cpp

@@ -415,6 +415,31 @@ namespace smt {
             }
         }
 
+        void internalize_mul(app* t) {
+            SASSERT(a.is_mul(t));
+            mk_enode(t);
+            theory_var v = mk_var(t);
+            svector<lean::var_index> vars;
+            ptr_vector<expr> todo;
+            todo.push_back(t);
+            while (!todo.empty()) {
+                expr* n = todo.back();
+                todo.pop_back();
+                expr* n1, *n2;
+                if (a.is_mul(n, n1, n2)) {
+                    todo.push_back(n1);
+                    todo.push_back(n2);
+                }
+                else {
+                    if (!ctx().e_internalized(n)) {
+                        ctx().internalize(n, false);
+                    }
+                    vars.push_back(get_var_index(mk_var(n)));
+                }
+            }
+            // m_solver->add_monomial(get_var_index(v), vars);
+        }
+
         enode * mk_enode(app * n) {
             if (ctx().e_internalized(n)) {
                 return get_enode(n);
@@ -1166,18 +1191,41 @@ namespace smt {
             else if (m_solver->get_status() != lean::lp_status::OPTIMAL) {
                 is_sat = make_feasible();
             }
+            final_check_status st = FC_DONE;
             switch (is_sat) {
             case l_true:
+                
                 if (delayed_assume_eqs()) {
                     return FC_CONTINUE;
                 }
                 if (assume_eqs()) {
                     return FC_CONTINUE;
                 }
-                if (m_not_handled != 0) {                    
-                    return FC_GIVEUP;
+
+                switch (check_lia()) {
+                case l_true:
+                    break;
+                case l_false:
+                    return FC_CONTINUE;
+                case l_undef:
+                    st = FC_GIVEUP;
+                    break;
                 }
-                return FC_DONE;
+                
+                switch (check_nra()) {
+                case l_true:
+                    break;
+                case l_false:
+                    return FC_CONTINUE;
+                case l_undef:
+                    st = FC_GIVEUP;
+                    break;
+                }
+                if (m_not_handled != 0) {                    
+                    st = FC_GIVEUP;
+                }
+
+                return st;
             case l_false:
                 set_conflict();
                 return FC_CONTINUE;
@@ -1190,6 +1238,28 @@ namespace smt {
             return FC_GIVEUP;
         }
 
+        lbool check_lia() {
+            if (m.canceled()) return l_undef;
+            return l_true;
+        }
+
+        lbool check_nra() {
+            if (m.canceled()) return l_undef;
+            return l_true;
+            // TBD:
+            switch (m_solver->check_nra(m_variable_values, m_explanation)) {
+            case lean::final_check_status::DONE:
+                return l_true;
+            case lean::final_check_status::CONTINUE:
+                return l_true; // ?? why have a continue at this level ??
+            case lean::final_check_status::UNSAT: 
+                set_conflict1();
+                return l_false;
+            case lean::final_check_status::GIVEUP:
+                return l_undef;
+            }
+            return l_true;
+        }
 
         /**
            \brief We must redefine this method, because theory of arithmetic contains
@@ -2197,11 +2267,15 @@ namespace smt {
         }
 
         void set_conflict() {
+            m_solver->get_infeasibility_explanation(m_explanation);
+            set_conflict1();
+        }
+
+        void set_conflict1() {
             m_eqs.reset();
             m_core.reset();
             m_params.reset();
             m_explanation.clear();
-            m_solver->get_infeasibility_explanation(m_explanation);
             // m_solver->shrink_explanation_to_minimum(m_explanation); // todo, enable when perf is fixed
             /*
             static unsigned cn = 0;

src/util/lp/lar_solver.cpp

@@ -31,7 +31,9 @@ lar_solver::lar_solver() : m_status(OPTIMAL),
                            m_infeasible_column_index(-1),
                            m_terms_start_index(1000000),
                            m_mpq_lar_core_solver(m_settings, *this)
-{}
+{
+    m_nra = alloc(nra::solver, *this);
+}
     
 void lar_solver::set_propagate_bounds_on_pivoted_rows_mode(bool v) {
     m_mpq_lar_core_solver.m_r_solver.m_pivoted_rows = v? (& m_rows_with_changed_bounds) : nullptr;
@@ -330,6 +332,7 @@ void lar_solver::push() {
     m_term_count.push();
     m_constraint_count = m_constraints.size();
     m_constraint_count.push();
+    m_nra->push();
 }
 
 void lar_solver::clean_large_elements_after_pop(unsigned n, int_set& set) {
@@ -385,6 +388,7 @@ void lar_solver::pop(unsigned k) {
     m_settings.simplex_strategy() = m_simplex_strategy;
     lean_assert(sizes_are_correct());
     lean_assert((!m_settings.use_tableau()) || m_mpq_lar_core_solver.m_r_solver.reduced_costs_are_correct_tableau());
+    m_nra->pop(k);
 }
     
 vector<constraint_index> lar_solver::get_all_constraint_indices() const {
@@ -1084,6 +1088,15 @@ void lar_solver::get_infeasibility_explanation(vector<std::pair<mpq, constraint_
     lean_assert(explanation_is_correct(explanation));
 }
 
+final_check_status lar_solver::check_nra(nra_model_t& model, explanation_t& explanation) {
+    return m_nra->check(model, explanation);
+}
+
+void lar_solver::add_monomial(var_index v, svector<var_index> const& vars) {
+    m_nra->add_monomial(v, vars.size(), vars.c_ptr());
+}
+
+
 void lar_solver::get_infeasibility_explanation_for_inf_sign(
                                                             vector<std::pair<mpq, constraint_index>> & explanation,
                                                             const vector<std::pair<mpq, unsigned>> & inf_row,

src/util/lp/lar_solver.h

@@ -31,9 +31,11 @@
 #include "util/lp/quick_xplain.h"
 #include "util/lp/conversion_helper.h"
 #include "util/lp/int_solver.h"
+#include "util/lp/nra_solver.h"
 
 namespace lean {
 
+
 class lar_solver : public column_namer {
     class ext_var_info {
         unsigned m_ext_j; // the external index
@@ -61,7 +63,8 @@ class lar_solver : public column_namer {
     stacked_value<unsigned> m_term_count;
     vector<lar_term*> m_terms;
     const var_index m_terms_start_index;
-    indexed_vector<mpq> m_column_buffer;    
+    indexed_vector<mpq> m_column_buffer;
+    scoped_ptr<nra::solver> m_nra;
 public:
     lar_core_solver m_mpq_lar_core_solver;
     unsigned constraint_count() const;
@@ -200,10 +203,14 @@ public:
     void set_status(lp_status s);
 
     lp_status find_feasible_solution();
+   
+    final_check_status check_nra(nra_model_t& model, explanation_t& explanation);    
+
+    void add_monomial(var_index v, svector<var_index> const& vars);
     
     lp_status solve();
 
-    void fill_explanation_from_infeasible_column(vector<std::pair<mpq, constraint_index>> & evidence) const;
+    void fill_explanation_from_infeasible_column(explanation_t & evidence) const;
 
     
     unsigned get_total_iterations() const;

src/util/lp/lp_settings.h

@@ -18,11 +18,17 @@ typedef unsigned constraint_index;
 typedef unsigned row_index;
 
 enum class final_check_status {
-	DONE,
-	CONTINUE,
-	GIVEUP
+    DONE,
+    CONTINUE,
+    UNSAT,
+    GIVEUP
 };
 
+typedef vector<std::pair<mpq, constraint_index>> explanation_t;
+
+typedef std::unordered_map<lean::var_index, rational> nra_model_t;
+
+
 enum class column_type  {
     free_column = 0,
         low_bound = 1,

src/util/lp/nra_solver.cpp

@@ -4,64 +4,254 @@
 */
 
 #pragma once
+#include "util/lp/lar_solver.h"
 #include "util/lp/nra_solver.h"
+#include "nlsat/nlsat_solver.h"
+#include "math/polynomial/polynomial.h"
+#include "math/polynomial/algebraic_numbers.h"
+#include "util/map.h"
 
-namespace lp {
 
-    struct nra_solver::imp {
+namespace nra {
+
+    struct solver::imp {
         lean::lar_solver& s;
+        reslimit      m_limit;  // TBD: extract from lar_solver
+        params_ref    m_params; // TBD: pass from outside
+        u_map<polynomial::var> m_lp2nl;  // map from lar_solver variables to nlsat::solver variables
 
-        svector<lean::var_index>         m_vars;
-        vector<svector<lean::var_index>> m_monomials;
+        struct mon_eq {
+            mon_eq(lean::var_index v, svector<lean::var_index> const& vs):
+                m_v(v), m_vs(vs) {}
+            lean::var_index          m_v;
+            svector<lean::var_index> m_vs;
+        };
+
+        vector<mon_eq>                   m_monomials;
         unsigned_vector                  m_lim;
+        mutable std::unordered_map<lean::var_index, rational> m_variable_values; // current model
 
-        imp(lean::lar_solver& s): s(s) {
-            m_lim.push_back(0);
+        imp(lean::lar_solver& s): 
+            s(s) {
         }
 
-        lean::final_check_status check_feasible() {
-            return lean::final_check_status::GIVEUP;
+        lean::final_check_status check_feasible(lean::nra_model_t& m, lean::explanation_t& ex) {
+            if (m_monomials.empty()) {
+                return lean::final_check_status::DONE;
+            }
+            if (check_assignments()) {
+                return lean::final_check_status::DONE;
+            }
+            switch (check_nlsat(m, ex)) {
+            case l_undef: return lean::final_check_status::GIVEUP;
+            case l_true: lean::final_check_status::DONE;
+            case l_false: lean::final_check_status::UNSAT;
+            }
+            return lean::final_check_status::DONE;
         }
 
         void add(lean::var_index v, unsigned sz, lean::var_index const* vs) {
-            m_vars.push_back(v);
-            m_monomials.push_back(svector<lean::var_index>(sz, vs));
+            m_monomials.push_back(mon_eq(v, svector<lean::var_index>(sz, vs)));
         }
 
         void push() {
-            m_lim.push_back(m_vars.size());
+            m_lim.push_back(m_monomials.size());
         }
 
         void pop(unsigned n) {
+            if (n == 0) return;
             SASSERT(n < m_lim.size());
+            m_monomials.shrink(m_lim[m_lim.size() - n]);
             m_lim.shrink(m_lim.size() - n);       
-            m_vars.shrink(m_lim.back());
-            m_monomials.shrink(m_lim.back());
+        }
+
+        /*
+          \brief Check if polynomials are well defined.
+          multiply values for vs and check if they are equal to value for v.
+          epsilon has been computed.
+        */
+        bool check_assignment(mon_eq const& m) const {
+            rational r1 = m_variable_values[m.m_v];
+            rational r2(1);
+            for (auto w : m.m_vs) {
+                r2 *= m_variable_values[w];
+            }
+            return r1 == r2;
+        }
+
+        bool check_assignments() const {
+            s.get_model(m_variable_values);
+            for (auto const& m : m_monomials) {
+                if (!check_assignment(m)) return false;
+            }
+            return true;
+        }
+
+
+        /**
+           \brief one-shot nlsat check.
+           A one shot checker is the least functionality that can 
+           enable non-linear reasoning. 
+           In addition to checking satisfiability we would also need
+           to identify equalities in the model that should be assumed
+           with the remaining solver.
+           
+           TBD: use partial model from lra_solver to prime the state of nlsat_solver.
+        */
+        lbool check_nlsat(lean::nra_model_t& model, lean::explanation_t& ex) {
+            nlsat::solver solver(m_limit, m_params);
+            m_lp2nl.reset();
+            // add linear inequalities from lra_solver
+            for (unsigned i = 0; i < s.constraint_count(); ++i) {
+                add_constraint(solver, i);
+            }
+
+            // add polynomial definitions.
+            for (auto const& m : m_monomials) {
+                add_monomial_eq(solver, m);
+            }
+            // TBD: add variable bounds?
+
+            lbool r = solver.check(); 
+            switch (r) {
+            case l_true: {
+                nlsat::anum_manager& am = solver.am();
+                model.clear();
+                for (auto kv : m_lp2nl) {
+                    kv.m_key;
+                    nlsat::anum const& v = solver.value(kv.m_value);
+                    if (is_int(kv.m_key) && !am.is_int(v)) {
+                        // the nlsat solver should already have returned unknown.
+                        TRACE("lp", tout << "Value is not integer " << kv.m_key << "\n";);
+                        return l_undef;
+                    }
+                    if (!am.is_rational(v)) {
+                        // TBD extract and convert model.
+                        TRACE("lp", tout << "Cannot handle algebraic numbers\n";);
+                        return l_undef;
+                    }
+                    rational r;
+                    am.to_rational(v, r);
+                    model[kv.m_key] = r;
+                }
+                break;
+            }
+            case l_false: {
+                ex.reset();
+                vector<nlsat::assumption, false> core;
+                solver.get_core(core);
+                for (auto c : core) {
+                    unsigned idx = static_cast<unsigned>(static_cast<imp*>(c) - this);
+                    ex.push_back(std::pair<rational, unsigned>(rational(1), idx));
+                }
+                break;
+            }
+            case l_undef:
+                break;
+            }            
+            return r;
+        }                
+
+        void add_monomial_eq(nlsat::solver& solver, mon_eq const& m) {
+            polynomial::manager& pm = solver.pm();
+            svector<polynomial::var> vars;
+            for (auto v : m.m_vs) {
+                vars.push_back(lp2nl(solver, v));
+            }
+            polynomial::monomial_ref m1(pm.mk_monomial(vars.size(), vars.c_ptr()), pm);
+            polynomial::monomial_ref m2(pm.mk_monomial(lp2nl(solver, m.m_v), 1), pm);
+            polynomial::monomial* mls[2] = { m1, m2 };
+            polynomial::scoped_numeral_vector coeffs(pm.m());
+            coeffs.push_back(mpz(1));
+            coeffs.push_back(mpz(-1));
+            polynomial::polynomial_ref p(pm.mk_polynomial(2, coeffs.c_ptr(), mls),  pm);
+            polynomial::polynomial* ps[1] = { p };
+            bool even[1] = { false };
+            nlsat::literal lit = solver.mk_ineq_literal(nlsat::atom::kind::EQ, 1, ps, even);
+            solver.mk_clause(1, &lit, 0);
+        }
+
+        void add_constraint(nlsat::solver& solver, unsigned idx) {
+            lean::lar_base_constraint const& c = s.get_constraint(idx);
+            polynomial::manager& pm = solver.pm();
+            auto k = c.m_kind;
+            auto rhs = c.m_right_side;
+            auto lhs = c.get_left_side_coefficients();
+            unsigned sz = lhs.size();
+            svector<polynomial::var> vars;
+            rational den = denominator(rhs);
+            for (auto kv : lhs) {
+                vars.push_back(lp2nl(solver, kv.second));
+                den = lcm(den, denominator(kv.first));
+            }
+            vector<rational> coeffs;
+            for (auto kv : lhs) {
+                coeffs.push_back(den * kv.first);
+            }
+            rhs *= den;
+            polynomial::polynomial_ref p(pm.mk_linear(sz, coeffs.c_ptr(), vars.c_ptr(), -rhs), pm);
+            polynomial::polynomial* ps[1] = { p };
+            bool is_even[1] = { false };
+            nlsat::literal lit;
+            switch (k) {
+            case lean::lconstraint_kind::LE:
+                lit = ~solver.mk_ineq_literal(nlsat::atom::kind::GT, 1, ps, is_even);
+                break;
+            case lean::lconstraint_kind::GE:
+                lit = ~solver.mk_ineq_literal(nlsat::atom::kind::LT, 1, ps, is_even);
+                break;
+            case lean::lconstraint_kind::LT:
+                lit = solver.mk_ineq_literal(nlsat::atom::kind::LT, 1, ps, is_even);
+                break;
+            case lean::lconstraint_kind::GT:
+                lit = solver.mk_ineq_literal(nlsat::atom::kind::GT, 1, ps, is_even);
+                break;
+            case lean::lconstraint_kind::EQ:
+                lit = solver.mk_ineq_literal(nlsat::atom::kind::EQ, 1, ps, is_even);
+                break;
+            }
+
+            nlsat::assumption a = this + idx;
+            solver.mk_clause(1, &lit, a);
+        }               
+
+        bool is_int(lean::var_index v) {
+            // TBD: is it s.column_is_integer(v), if then the function should take a var_index and not unsigned; s.is_int(v);
+            return false;
+        }
+
+        polynomial::var lp2nl(nlsat::solver& solver, lean::var_index v) {
+            polynomial::var r;
+            if (!m_lp2nl.find(v, r)) {
+                r = solver.mk_var(is_int(v));
+                m_lp2nl.insert(v, r);
+            }
+            return r;
         }
 
     };
 
-    nra_solver::nra_solver(lean::lar_solver& s) {
+    solver::solver(lean::lar_solver& s) {
         m_imp = alloc(imp, s);
     }
 
-    nra_solver::~nra_solver() {
+    solver::~solver() {
         dealloc(m_imp);
     }
 
-    void nra_solver::add_monomial(lean::var_index v, unsigned sz, lean::var_index const* vs) {
+    void solver::add_monomial(lean::var_index v, unsigned sz, lean::var_index const* vs) {
         m_imp->add(v, sz, vs);
     }
 
-    lean::final_check_status nra_solver::check_feasible() {
-        return m_imp->check_feasible();
+    lean::final_check_status solver::check(lean::nra_model_t& m, lean::explanation_t& ex) {
+        return m_imp->check_feasible(m, ex);
     }
 
-    void nra_solver::push() {
+    void solver::push() {
         m_imp->push();
     }
 
-    void nra_solver::pop(unsigned n) {
+    void solver::pop(unsigned n) {
         m_imp->pop(n);
     }
 }

src/util/lp/nra_solver.h

@@ -6,23 +6,22 @@
 #pragma once
 #include "util/vector.h"
 #include "util/lp/lp_settings.h"
-#include "util/lp/lar_solver.h"
 
 namespace lean {
     class lar_solver;
 }
 
-namespace lp {
+namespace nra {
 
-    class nra_solver {
+    class solver {
         struct imp;
         imp* m_imp;
 
     public:
 
-        nra_solver(lean::lar_solver& s);
+        solver(lean::lar_solver& s);
         
-        ~nra_solver();
+        ~solver();
 
         /*
           \brief Add a definition v = vs[0]*vs[1]*...*vs[sz-1]
@@ -34,7 +33,7 @@ namespace lp {
           \brief Check feasiblity of linear constraints augmented by polynomial definitions
           that are added.
          */
-        lean::final_check_status check_feasible();
+        lean::final_check_status check(lean::nra_model_t& m, lean::explanation_t& ex);
 
         /*
           \brief push and pop scope.