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/smt/theory_lra.cpp

@@ -1166,18 +1166,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 +1213,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 +2242,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

@@ -330,6 +330,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 +386,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 +1086,10 @@ 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::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,12 @@ public:
     void set_status(lp_status s);
 
     lp_status find_feasible_solution();
+   
+    final_check_status check_nra(nra_model_t& model, explanation_t& explanation);    
     
     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,
+    DONE,
-	CONTINUE,
+    CONTINUE,
-	GIVEUP
+    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,14 +4,17 @@
 */
 
 #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
@@ -27,13 +30,25 @@ namespace lp {
         vector<mon_eq>                   m_monomials;
         unsigned_vector                  m_lim;
         mutable std::unordered_map<lean::var_index, rational> m_variable_values; // current model
+        vector<std::pair<rational, unsigned>> m_core;
 
         imp(lean::lar_solver& s): 
             s(s) {
         }
 
-        lean::final_check_status check_feasible() {
+        lean::final_check_status check_feasible(lean::nra_model_t& m, lean::explanation_t& ex) {
-            return lean::final_check_status::GIVEUP;
+            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) {
@@ -84,21 +99,57 @@ namespace lp {
            
            TBD: use partial model from lra_solver to prime the state of nlsat_solver.
         */
-        lbool check_nlsat() {
+        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, s.get_constraint(i));
+                add_constraint(solver, i);
             }
+
             // add polynomial definitions.
             for (auto const& m : m_monomials) {
                 add_monomial_eq(solver, m);
             }
-            lbool r = solver.check(); // TBD: get assumptions from literals that are asserted above level 0.
+            // TBD: add variable bounds?
-            if (r == l_true) {
+
-                // TBD extract model.
+            lbool r = solver.check(); 
-                // check interface equalities
+            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;
         }                
 
@@ -121,7 +172,8 @@ namespace lp {
             solver.mk_clause(1, &lit, 0);
         }
 
-        void add_constraint(nlsat::solver& solver, lean::lar_base_constraint const& c) {
+        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;
@@ -139,58 +191,68 @@ namespace lp {
             }
             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, );
+                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;
             }
 
-            // solver.mk_clause();
+            nlsat::assumption a = this + idx;
+            solver.mk_clause(1, &lit, a);
+        }               
 
-            // c.get_free_coeff_of_left_side();
+        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;
         }
 
-        
-        // translate var_index into polynomial::var that are declared on nlsat::solver.
-        
-
         polynomial::var lp2nl(nlsat::solver& solver, lean::var_index v) {
             polynomial::var r;
             if (!m_lp2nl.find(v, r)) {
-                r = solver.mk_var(false); // TBD: is it s.column_is_integer(v), if then the function should take a var_index and not unsigned; s.is_int(v);
+                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() {
+    lean::final_check_status solver::check(lean::nra_model_t& m, lean::explanation_t& ex) {
-        return m_imp->check_feasible();
+        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. 
@@ -43,5 +42,6 @@ namespace lp {
         void push();
         
         void pop(unsigned n);
+        
     };
 }