merge

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

src/util/lp/init_lar_solver.cpp

@@ -2,40 +2,44 @@
   Copyright (c) 2017 Microsoft Corporation
   Author: Lev Nachmanson
 */
+// this file represents the intiialization functionality of lar_solver
 
-// here we are inside lean::lar_solver class
 
-bool strategy_is_undecided() const {
+#include "util/lp/lar_solver.h"
+namespace lean {
+
+bool lar_solver::strategy_is_undecided() const {
     return m_settings.simplex_strategy() == simplex_strategy_enum::undecided;
 }
 
-var_index add_var(unsigned ext_j) {
+var_index lar_solver::add_var(unsigned ext_j, bool is_integer) {
     var_index i;
     lean_assert (ext_j < m_terms_start_index); 
 
     if (ext_j >= m_terms_start_index)
         throw 0; // todo : what is the right way to exit?
-            
-    if (try_get_val(m_ext_vars_to_columns, ext_j, i)) {
-        return i;
+    auto it = m_ext_vars_to_columns.find(ext_j);
+    if (it != m_ext_vars_to_columns.end()) {
+        return it->second.ext_j();
     }
     lean_assert(m_vars_to_ul_pairs.size() == A_r().column_count());
     i = A_r().column_count();
     m_vars_to_ul_pairs.push_back (ul_pair(static_cast<unsigned>(-1)));
     add_non_basic_var_to_core_fields(ext_j);
     lean_assert(sizes_are_correct());
+    lean_assert(!column_is_integer(i));
     return i;
 }
 
-void register_new_ext_var_index(unsigned ext_v) {
+void lar_solver::register_new_ext_var_index(unsigned ext_v) {
     lean_assert(!contains(m_ext_vars_to_columns, ext_v));
     unsigned j = static_cast<unsigned>(m_ext_vars_to_columns.size());
-    m_ext_vars_to_columns[ext_v] = j;
+    m_ext_vars_to_columns.insert(std::make_pair(ext_v, ext_var_info(j)));
     lean_assert(m_columns_to_ext_vars_or_term_indices.size() == j);
     m_columns_to_ext_vars_or_term_indices.push_back(ext_v);
 }
 
-void add_non_basic_var_to_core_fields(unsigned ext_j) {
+void lar_solver::add_non_basic_var_to_core_fields(unsigned ext_j) {
     register_new_ext_var_index(ext_j);
     m_mpq_lar_core_solver.m_column_types.push_back(column_type::free_column);
     m_columns_with_changed_bound.increase_size_by_one();
@@ -44,7 +48,7 @@ void add_non_basic_var_to_core_fields(unsigned ext_j) {
         add_new_var_to_core_fields_for_doubles(false);
 }
 
-void add_new_var_to_core_fields_for_doubles(bool register_in_basis) {
+void lar_solver::add_new_var_to_core_fields_for_doubles(bool register_in_basis) {
     unsigned j = A_d().column_count();
     A_d().add_column();
     lean_assert(m_mpq_lar_core_solver.m_d_x.size() == j);
@@ -63,7 +67,7 @@ void add_new_var_to_core_fields_for_doubles(bool register_in_basis) {
     }
 }
 
-void add_new_var_to_core_fields_for_mpq(bool register_in_basis) {
+void lar_solver::add_new_var_to_core_fields_for_mpq(bool register_in_basis) {
     unsigned j = A_r().column_count();
     A_r().add_column();
     lean_assert(m_mpq_lar_core_solver.m_r_x.size() == j);
@@ -88,25 +92,23 @@ void add_new_var_to_core_fields_for_mpq(bool register_in_basis) {
 }
 
 
-var_index add_term_undecided(const vector<std::pair<mpq, var_index>> & coeffs,
+var_index lar_solver::add_term_undecided(const vector<std::pair<mpq, var_index>> & coeffs,
                              const mpq &m_v) {
     m_terms.push_back(new lar_term(coeffs, m_v));
-    m_orig_terms.push_back(new lar_term(coeffs, m_v));
     return m_terms_start_index +  m_terms.size() - 1;
 }
 
 // terms
-var_index add_term(const vector<std::pair<mpq, var_index>> & coeffs,
+var_index lar_solver::add_term(const vector<std::pair<mpq, var_index>> & coeffs,
                    const mpq &m_v) {
     if (strategy_is_undecided())
         return add_term_undecided(coeffs, m_v);
 
     m_terms.push_back(new lar_term(coeffs, m_v));
-    m_orig_terms.push_back(new lar_term(coeffs, m_v));
     unsigned adjusted_term_index = m_terms.size() - 1;
     var_index ret = m_terms_start_index + adjusted_term_index;
     if (use_tableau() && !coeffs.empty()) {
-        add_row_for_term(m_orig_terms.back(), ret);
+        add_row_for_term(m_terms.back(), ret);
         if (m_settings.bound_propagation())
             m_rows_with_changed_bounds.insert(A_r().row_count() - 1);
     } 
@@ -114,13 +116,13 @@ var_index add_term(const vector<std::pair<mpq, var_index>> & coeffs,
     return ret;
 }
 
-void add_row_for_term(const lar_term * term, unsigned term_ext_index) {
+void lar_solver::add_row_for_term(const lar_term * term, unsigned term_ext_index) {
     lean_assert(sizes_are_correct());
     add_row_from_term_no_constraint(term, term_ext_index);
     lean_assert(sizes_are_correct());
 }
 
-void add_row_from_term_no_constraint(const lar_term * term, unsigned term_ext_index) {
+void lar_solver::add_row_from_term_no_constraint(const lar_term * term, unsigned term_ext_index) {
     register_new_ext_var_index(term_ext_index);
     // j will be a new variable
 	unsigned j = A_r().column_count();
@@ -140,7 +142,7 @@ void add_row_from_term_no_constraint(const lar_term * term, unsigned term_ext_in
         fill_last_row_of_A_d(A_d(), term);
 }
 
-void add_basic_var_to_core_fields() {
+void lar_solver::add_basic_var_to_core_fields() {
     bool use_lu = m_mpq_lar_core_solver.need_to_presolve_with_double_solver();
     lean_assert(!use_lu || A_r().column_count() == A_d().column_count());
     m_mpq_lar_core_solver.m_column_types.push_back(column_type::free_column);
@@ -151,7 +153,7 @@ void add_basic_var_to_core_fields() {
         add_new_var_to_core_fields_for_doubles(true);
 }
 
-constraint_index add_var_bound(var_index j, lconstraint_kind kind, const mpq & right_side)  {
+constraint_index lar_solver::add_var_bound(var_index j, lconstraint_kind kind, const mpq & right_side)  {
 	constraint_index ci = m_constraints.size();
     if (!is_term(j)) { // j is a var
         auto vc = new lar_var_constraint(j, kind, right_side);
@@ -164,7 +166,7 @@ constraint_index add_var_bound(var_index j, lconstraint_kind kind, const mpq & r
     return ci;
 }
 
-void update_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_index) {
+void lar_solver::update_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_index) {
     switch(m_mpq_lar_core_solver.m_column_types[j]) {
     case column_type::free_column:
         update_free_column_type_and_bound(j, kind, right_side, constr_index);
@@ -186,22 +188,32 @@ void update_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq
     }
 }
 
-void add_var_bound_on_constraint_for_term(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
+void lar_solver::add_var_bound_on_constraint_for_term(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
     lean_assert(is_term(j));
     unsigned adjusted_term_index = adjust_term_index(j);
-    unsigned term_j;
-    if (try_get_val(m_ext_vars_to_columns, j, term_j)) {
-        mpq rs = right_side - m_orig_terms[adjusted_term_index]->m_v;
-        m_constraints.push_back(new lar_term_constraint(m_orig_terms[adjusted_term_index], kind, right_side));
+    auto it = m_ext_vars_to_columns.find(j);
+    if (it != m_ext_vars_to_columns.end()) {
+        unsigned term_j = it->second.ext_j();
+        mpq rs = right_side - m_terms[adjusted_term_index]->m_v;
+        m_constraints.push_back(new lar_term_constraint(m_terms[adjusted_term_index], kind, right_side));
         update_column_type_and_bound(term_j, kind, rs, ci);
     }
     else {
-        add_constraint_from_term_and_create_new_column_row(j, m_orig_terms[adjusted_term_index], kind, right_side);
+        add_constraint_from_term_and_create_new_column_row(j, m_terms[adjusted_term_index], kind, right_side);
     }
 }
 
+constraint_index lar_solver::add_constraint(const vector<std::pair<mpq, var_index>>& left_side_with_terms, lconstraint_kind kind_par, const mpq& right_side_parm) {
+    vector<std::pair<mpq, var_index>> left_side;
+    mpq rs = - right_side_parm;
+    substitute_terms_in_linear_expression(left_side_with_terms, left_side, rs);
+    unsigned term_index = add_term(left_side, zero_of_type<mpq>());
+    constraint_index ci = m_constraints.size();
+    add_var_bound_on_constraint_for_term(term_index, kind_par, -rs, ci);
+    return ci;
+}
 
-void add_constraint_from_term_and_create_new_column_row(unsigned term_j, const lar_term* term,
+void lar_solver::add_constraint_from_term_and_create_new_column_row(unsigned term_j, const lar_term* term,
                                                         lconstraint_kind kind, const mpq & right_side) {
 
     add_row_from_term_no_constraint(term, term_j);
@@ -211,7 +223,7 @@ void add_constraint_from_term_and_create_new_column_row(unsigned term_j, const l
     lean_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_costs.size());
 }
 
-void decide_on_strategy_and_adjust_initial_state() {
+void lar_solver::decide_on_strategy_and_adjust_initial_state() {
 	lean_assert(strategy_is_undecided());
 	if (m_vars_to_ul_pairs.size() > m_settings.column_number_threshold_for_using_lu_in_lar_solver) {
         m_settings.simplex_strategy() = simplex_strategy_enum::lu;
@@ -221,7 +233,7 @@ void decide_on_strategy_and_adjust_initial_state() {
     adjust_initial_state();
 }
 
-void adjust_initial_state() {
+void lar_solver::adjust_initial_state() {
     switch (m_settings.simplex_strategy()) {
     case simplex_strategy_enum::lu:
         adjust_initial_state_for_lu();
@@ -237,7 +249,7 @@ void adjust_initial_state() {
     }
 }
 
-void adjust_initial_state_for_lu() {
+void lar_solver::adjust_initial_state_for_lu() {
     copy_from_mpq_matrix(A_d());
 	unsigned n = A_d().column_count();
 	m_mpq_lar_core_solver.m_d_x.resize(n);
@@ -265,7 +277,7 @@ void adjust_initial_state_for_lu() {
         }*/
 }
 
-void adjust_initial_state_for_tableau_rows() {
+void lar_solver::adjust_initial_state_for_tableau_rows() {
     for (unsigned j = 0; j < m_terms.size(); j++) {
         if (contains(m_ext_vars_to_columns, j + m_terms_start_index))
             continue;
@@ -274,7 +286,7 @@ void adjust_initial_state_for_tableau_rows() {
 }
 
 // this fills the last row of A_d and sets the basis column: -1 in the last column of the row
-void fill_last_row_of_A_d(static_matrix<double, double> & A, const lar_term* ls) {
+void lar_solver::fill_last_row_of_A_d(static_matrix<double, double> & A, const lar_term* ls) {
     lean_assert(A.row_count() > 0);
     lean_assert(A.column_count() > 0);
     unsigned last_row = A.row_count() - 1;
@@ -290,7 +302,7 @@ void fill_last_row_of_A_d(static_matrix<double, double> & A, const lar_term* ls)
     A.set(last_row, basis_j, - 1 );
 }
 
-void update_free_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_ind) {
+void lar_solver::update_free_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_ind) {
     mpq y_of_bound(0);
     switch (kind) {
     case LT:
@@ -330,7 +342,7 @@ void update_free_column_type_and_bound(var_index j, lconstraint_kind kind, const
     m_columns_with_changed_bound.insert(j);
 }
 
-void update_upper_bound_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
+void lar_solver::update_upper_bound_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
     lean_assert(m_mpq_lar_core_solver.m_column_types()[j] == column_type::upper_bound);
     mpq y_of_bound(0);
     switch (kind) {
@@ -387,7 +399,7 @@ void update_upper_bound_column_type_and_bound(var_index j, lconstraint_kind kind
     }
 }
     
-void update_boxed_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
+void lar_solver::update_boxed_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
     lean_assert(m_status == INFEASIBLE || (m_mpq_lar_core_solver.m_column_types()[j] == column_type::boxed && m_mpq_lar_core_solver.m_r_low_bounds()[j] < m_mpq_lar_core_solver.m_r_upper_bounds()[j]));
     mpq y_of_bound(0);
     switch (kind) {
@@ -457,7 +469,7 @@ void update_boxed_column_type_and_bound(var_index j, lconstraint_kind kind, cons
                 
     }
 }
-void update_low_bound_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
+void lar_solver::update_low_bound_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
     lean_assert(m_mpq_lar_core_solver.m_column_types()[j] == column_type::low_bound);
     mpq y_of_bound(0);
     switch (kind) {
@@ -513,7 +525,7 @@ void update_low_bound_column_type_and_bound(var_index j, lconstraint_kind kind,
     }
 }
 
-void update_fixed_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
+void lar_solver::update_fixed_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
     lean_assert(m_status == INFEASIBLE || (m_mpq_lar_core_solver.m_column_types()[j] == column_type::fixed && m_mpq_lar_core_solver.m_r_low_bounds()[j] == m_mpq_lar_core_solver.m_r_upper_bounds()[j]));
     lean_assert(m_status == INFEASIBLE || (m_mpq_lar_core_solver.m_r_low_bounds()[j].y.is_zero() && m_mpq_lar_core_solver.m_r_upper_bounds()[j].y.is_zero()));
     auto v = numeric_pair<mpq>(right_side, mpq(0));
@@ -574,3 +586,4 @@ void update_fixed_column_type_and_bound(var_index j, lconstraint_kind kind, cons
     }
 }
     
+}

src/util/lp/int_solver.cpp

@@ -0,0 +1,6 @@
+/*
+  Copyright (c) 2017 Microsoft Corporation
+  Author: Lev Nachmanson
+*/
+
+#include "util/lp/int_solver.h"

src/util/lp/int_solver.h

@@ -0,0 +1,44 @@
+/*
+  Copyright (c) 2017 Microsoft Corporation
+  Author: Lev Nachmanson
+*/
+#pragma once
+#include "util/lp/lp_settings.h"
+class lemma; // forward definition
+namespace lean {
+class lar_solver;
+template <typename T, typename X>
+struct lp_constraint;
+
+class int_solver {
+public:
+    lar_solver *m_solver;
+    
+    int_solver(lar_solver* lp);
+    bool check();// main function to check that solution provided by lar_solver is valid for integral values or can be adjusted.
+private:
+
+    // how to tighten bounds for integer variables.
+    
+    
+    // gcd test
+    // 5*x + 3*y + 6*z = 5
+    // suppose x is fixed at 2.
+    // so we have 10 + 3(y + 2z) = 5
+    //             5 = -3(y + 2z)
+    // this is unsolvable because 5/3 is not an integer.
+    // so we create a lemma that rules out this condition.
+    // 
+    bool gcd_test(lemma& lemma); // returns false in case of failure. Creates a theory lemma in case of failure.
+
+    // create goromy cuts
+    // either creates a conflict or a bound.
+
+    // branch and bound: 
+    // decide what to branch and bound on
+    // creates a fresh inequality.
+
+    bool branch(const lp_constraint<mpq, mpq> & new_inequality);
+    
+};
+}

src/util/lp/lar_solver_instances.cpp

@@ -0,0 +1,17 @@
+/*
+  Copyright (c) 2017 Microsoft Corporation
+  Author: Lev Nachmanson
+*/
+
+#include "util/lp/lar_solver.cpp"
+#include "util/lp/init_lar_solver.cpp"
+
+
+
+
+template void  lean::lar_solver::copy_from_mpq_matrix<double,double>(class lean::static_matrix<double,double> &);
+
+
+
+
+

src/util/lp/lp_settings.h

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

src/util/lp/mps_reader.h

@@ -810,7 +810,7 @@ public:
             auto kind = get_lar_relation_from_row(row->m_type);
             vector<std::pair<mpq, var_index>> ls;
             for (auto s : row->m_row_columns) {
-                var_index i = solver->add_var(get_var_index(s.first));
+                var_index i = solver->add_var(get_var_index(s.first), false);
                 ls.push_back(std::make_pair(s.second, i));
             }
             solver->add_constraint(ls, kind, row->m_right_side);
@@ -828,20 +828,20 @@ public:
 
     void create_low_constraint_for_var(column* col, bound * b, lar_solver *solver) {
         vector<std::pair<mpq, var_index>> ls;
-        var_index i = solver->add_var(col->m_index);
+        var_index i = solver->add_var(col->m_index, false);
         ls.push_back(std::make_pair(numeric_traits<T>::one(), i));
         solver->add_constraint(ls, GE, b->m_low);
     }
 
     void create_upper_constraint_for_var(column* col, bound * b, lar_solver *solver) {
-        var_index i = solver->add_var(col->m_index);
+        var_index i = solver->add_var(col->m_index, false);
         vector<std::pair<mpq, var_index>> ls;
         ls.push_back(std::make_pair(numeric_traits<T>::one(), i));
         solver->add_constraint(ls, LE, b->m_upper);
     }
 
     void create_equality_contraint_for_var(column* col, bound * b, lar_solver *solver) {
-        var_index i = solver->add_var(col->m_index);
+        var_index i = solver->add_var(col->m_index, false);
         vector<std::pair<mpq, var_index>> ls;
         ls.push_back(std::make_pair(numeric_traits<T>::one(), i));
         solver->add_constraint(ls, EQ, b->m_fixed_value);
@@ -850,7 +850,7 @@ public:
     void fill_lar_solver_on_columns(lar_solver * solver) {
         for (auto s : m_columns) {
             mps_reader::column * col = s.second;
-            solver->add_var(col->m_index);
+            solver->add_var(col->m_index, false);
             auto b = col->m_bound;
             if (b == nullptr) return;
 

src/util/lp/nra_solver.cpp

@@ -0,0 +1,262 @@
+/*
+  Copyright (c) 2017 Microsoft Corporation
+  Author: Lev Nachmanson
+*/
+
+#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 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
+        nlsat::solver m_nlsat;
+
+        struct mon_eq {
+            mon_eq(lean::var_index v, unsigned sz, lean::var_index const* vs):
+                m_v(v), m_vs(sz, 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, reslimit& lim, params_ref const& p): 
+            s(s), 
+            m_limit(lim),
+            m_params(p),
+            m_nlsat(m_limit, m_params) {
+        }
+
+        bool need_check() {
+            return !m_monomials.empty() && !check_assignments();
+        }
+
+        void add(lean::var_index v, unsigned sz, lean::var_index const* vs) {
+            m_monomials.push_back(mon_eq(v, sz, vs));
+        }
+
+        void push() {
+            m_lim.push_back(m_monomials.size());
+        }
+
+        void pop(unsigned n) {
+            if (n == 0) return;
+            m_monomials.shrink(m_lim[m_lim.size() - n]);
+            m_lim.shrink(m_lim.size() - n);       
+        }
+
+        /*
+          \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(lean::explanation_t& ex) {
+            SASSERT(need_check());
+            m_nlsat.reset();
+            m_lp2nl.reset();
+            vector<nlsat::assumption, false> core;
+
+            // add linear inequalities from lra_solver
+            for (unsigned i = 0; i < s.constraint_count(); ++i) {
+                add_constraint(i);
+            }
+
+            // add polynomial definitions.
+            for (auto const& m : m_monomials) {
+                add_monomial_eq(m);
+            }
+            // TBD: add variable bounds?
+
+            lbool r = m_nlsat.check(); 
+            TRACE("arith", m_nlsat.display(tout << r << "\n"););
+            switch (r) {
+            case l_true: 
+                break;
+            case l_false: 
+                ex.reset();
+                m_nlsat.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));
+                    TRACE("arith", tout << "ex: " << idx << "\n";);
+                }
+                break;
+
+            case l_undef:
+                break;
+            }            
+            return r;
+        }                
+
+        void add_monomial_eq(mon_eq const& m) {
+            polynomial::manager& pm = m_nlsat.pm();
+            svector<polynomial::var> vars;
+            for (auto v : m.m_vs) {
+                vars.push_back(lp2nl(v));
+            }
+            polynomial::monomial_ref m1(pm.mk_monomial(vars.size(), vars.c_ptr()), pm);
+            polynomial::monomial_ref m2(pm.mk_monomial(lp2nl(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 = m_nlsat.mk_ineq_literal(nlsat::atom::kind::EQ, 1, ps, even);
+            m_nlsat.mk_clause(1, &lit, 0);
+        }
+
+        void add_constraint(unsigned idx) {
+            auto& c = s.get_constraint(idx);
+            auto& pm = m_nlsat.pm();
+            auto k = c.m_kind;
+            auto rhs = c.m_right_side;
+            auto lhs = c.get_left_side_coefficients();
+            auto sz = lhs.size();
+            svector<polynomial::var> vars;
+            rational den = denominator(rhs);
+            for (auto kv : lhs) {
+                vars.push_back(lp2nl(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 = ~m_nlsat.mk_ineq_literal(nlsat::atom::kind::GT, 1, ps, is_even);
+                break;
+            case lean::lconstraint_kind::GE:
+                lit = ~m_nlsat.mk_ineq_literal(nlsat::atom::kind::LT, 1, ps, is_even);
+                break;
+            case lean::lconstraint_kind::LT:
+                lit = m_nlsat.mk_ineq_literal(nlsat::atom::kind::LT, 1, ps, is_even);
+                break;
+            case lean::lconstraint_kind::GT:
+                lit = m_nlsat.mk_ineq_literal(nlsat::atom::kind::GT, 1, ps, is_even);
+                break;
+            case lean::lconstraint_kind::EQ:
+                lit = m_nlsat.mk_ineq_literal(nlsat::atom::kind::EQ, 1, ps, is_even);
+                break;
+            }
+
+            nlsat::assumption a = this + idx;
+            m_nlsat.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(lean::var_index v) {
+            polynomial::var r;
+            if (!m_lp2nl.find(v, r)) {
+                r = m_nlsat.mk_var(is_int(v));
+                m_lp2nl.insert(v, r);
+            }
+            return r;
+        }
+
+        nlsat::anum const& value(lean::var_index v) const {
+           return m_nlsat.value(m_lp2nl.find(v));
+        }
+
+        std::ostream& display(std::ostream& out) const {
+            for (auto m : m_monomials) {
+                out << "v" << m.m_v << " = ";
+                for (auto v : m.m_vs) {
+                    out << "v" << v << " ";
+                }
+                out << "\n";
+            }
+            return out;
+        }
+
+    };
+
+    solver::solver(lean::lar_solver& s, reslimit& lim, params_ref const& p) {
+        m_imp = alloc(imp, s, lim, p);
+    }
+
+    solver::~solver() {
+        dealloc(m_imp);
+    }
+
+    void solver::add_monomial(lean::var_index v, unsigned sz, lean::var_index const* vs) {
+        m_imp->add(v, sz, vs);
+    }
+
+    lbool solver::check(lean::explanation_t& ex) {
+        return m_imp->check(ex);
+    }
+
+    bool solver::need_check() {
+        return m_imp->need_check();
+    }
+
+    void solver::push() {
+        m_imp->push();
+    }
+
+    void solver::pop(unsigned n) {
+        m_imp->pop(n);
+    }
+
+    std::ostream& solver::display(std::ostream& out) const {
+        return m_imp->display(out);
+    }
+
+    nlsat::anum const& solver::value(lean::var_index v) const {
+        return m_imp->value(v);
+    }
+
+}

src/util/lp/nra_solver.h

@@ -0,0 +1,68 @@
+/*
+  Copyright (c) 2017 Microsoft Corporation
+  Author: Lev Nachmanson
+*/
+
+#pragma once
+#include "util/vector.h"
+#include "util/lp/lp_settings.h"
+#include "util/rlimit.h"
+#include "util/params.h"
+#include "nlsat/nlsat_solver.h"
+
+namespace lean {
+    class lar_solver;
+}
+
+
+namespace nra {
+
+
+
+    class solver {
+        struct imp;
+        imp* m_imp;
+
+    public:
+
+        solver(lean::lar_solver& s, reslimit& lim, params_ref const& p = params_ref());
+        
+        ~solver();
+
+        /*
+          \brief Add a definition v = vs[0]*vs[1]*...*vs[sz-1]
+          The variable v is equal to the product of variables vs.
+        */
+        void add_monomial(lean::var_index v, unsigned sz, lean::var_index const* vs);
+
+        /*
+          \brief Check feasiblity of linear constraints augmented by polynomial definitions
+          that are added.
+         */
+        lbool check(lean::explanation_t& ex);
+
+        /*
+          \brief determine whether nra check is needed.
+        */
+        bool need_check();
+
+        /*
+          \brief Access model.
+        */
+        nlsat::anum const& value(lean::var_index v) const;
+
+        /*
+          \brief push and pop scope. 
+          Monomial definitions are retraced when popping scope.
+        */
+        void push();
+        
+        void pop(unsigned n);
+
+        /*
+          \brief display state
+         */
+        std::ostream& display(std::ostream& out) const;
+
+    };
+}

src/util/lp/quick_xplain.cpp

@@ -20,7 +20,7 @@ void quick_xplain::copy_constraint_and_add_constraint_vars(const lar_constraint&
     vector < std::pair<mpq, unsigned>> ls;
     for (auto & p : lar_c.get_left_side_coefficients()) {
         unsigned j = p.second;
-        unsigned lj = m_qsol.add_var(j);
+        unsigned lj = m_qsol.add_var(j, false);
         ls.push_back(std::make_pair(p.first, lj));
     }
     m_constraints_in_local_vars.push_back(lar_constraint(ls, lar_c.m_kind, lar_c.m_right_side));
@@ -94,7 +94,7 @@ bool quick_xplain::is_feasible(const vector<unsigned> & x, unsigned k) const {
         vector < std::pair<mpq, unsigned>> ls;
         const lar_constraint & c = m_constraints_in_local_vars[i];
         for (auto & p : c.get_left_side_coefficients()) {
-            unsigned lj = l.add_var(p.second);
+            unsigned lj = l.add_var(p.second, false);
             ls.push_back(std::make_pair(p.first, lj));
         }
         l.add_constraint(ls, c.m_kind, c.m_right_side);