From 4726d32e2f5cfb7c2f6cd057f680b049f0a19df0 Mon Sep 17 00:00:00 2001 From: Nikolaj Bjorner Date: Wed, 24 May 2017 16:32:14 -0700 Subject: [PATCH] Dev (#51) * initial skeletons for nra solver Signed-off-by: Nikolaj Bjorner * initial skeletons for nra solver Signed-off-by: Nikolaj Bjorner * adding more nlsat Signed-off-by: Nikolaj Bjorner * nlsat integration Signed-off-by: Nikolaj Bjorner * adding constraints Signed-off-by: Nikolaj Bjorner * adding nra solver Signed-off-by: Nikolaj Bjorner * add missing initialization Signed-off-by: Nikolaj Bjorner * adding nra solver Signed-off-by: Nikolaj Bjorner --- contrib/cmake/src/CMakeLists.txt | 2 +- contrib/cmake/src/smt/CMakeLists.txt | 1 + contrib/cmake/src/test/lp/CMakeLists.txt | 2 +- contrib/cmake/src/util/lp/CMakeLists.txt | 2 + src/nlsat/nlsat_solver.h | 8 +- src/nlsat/nlsat_types.h | 8 +- src/smt/theory_lra.cpp | 82 +++++++- src/util/lp/lar_solver.cpp | 15 +- src/util/lp/lar_solver.h | 11 +- src/util/lp/lp_settings.h | 12 +- src/util/lp/nra_solver.cpp | 230 +++++++++++++++++++++-- src/util/lp/nra_solver.h | 11 +- 12 files changed, 338 insertions(+), 46 deletions(-) diff --git a/contrib/cmake/src/CMakeLists.txt b/contrib/cmake/src/CMakeLists.txt index dd440b34d..a98f92ca3 100644 --- a/contrib/cmake/src/CMakeLists.txt +++ b/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) diff --git a/contrib/cmake/src/smt/CMakeLists.txt b/contrib/cmake/src/smt/CMakeLists.txt index 41890dd05..3db66eb3e 100644 --- a/contrib/cmake/src/smt/CMakeLists.txt +++ b/contrib/cmake/src/smt/CMakeLists.txt @@ -70,6 +70,7 @@ z3_add_component(smt euclid fpa grobner + nlsat lp macros normal_forms diff --git a/contrib/cmake/src/test/lp/CMakeLists.txt b/contrib/cmake/src/test/lp/CMakeLists.txt index 99f8747b1..6683a1758 100644 --- a/contrib/cmake/src/test/lp/CMakeLists.txt +++ b/contrib/cmake/src/test/lp/CMakeLists.txt @@ -1,4 +1,4 @@ -add_executable(lp_tst lp_main.cpp lp.cpp $ $) +add_executable(lp_tst lp_main.cpp lp.cpp $ $ $ $ ) 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}) diff --git a/contrib/cmake/src/util/lp/CMakeLists.txt b/contrib/cmake/src/util/lp/CMakeLists.txt index 75051d68c..e33764f69 100644 --- a/contrib/cmake/src/util/lp/CMakeLists.txt +++ b/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 ) diff --git a/src/nlsat/nlsat_solver.h b/src/nlsat/nlsat_solver.h index 3668629cd..08902e709 100644 --- a/src/nlsat/nlsat_solver.h +++ b/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 { diff --git a/src/nlsat/nlsat_types.h b/src/nlsat/nlsat_types.h index 11e063a17..70da98e32 100644 --- a/src/nlsat/nlsat_types.h +++ b/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; diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp index 0072aea7d..cd663f354 100644 --- a/src/smt/theory_lra.cpp +++ b/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 vars; + ptr_vector 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; diff --git a/src/util/lp/lar_solver.cpp b/src/util/lp/lar_solver.cpp index 14d3535af..b3da4e8ba 100644 --- a/src/util/lp/lar_solver.cpp +++ b/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 lar_solver::get_all_constraint_indices() const { @@ -1084,6 +1088,15 @@ void lar_solver::get_infeasibility_explanation(vectorcheck(model, explanation); +} + +void lar_solver::add_monomial(var_index v, svector const& vars) { + m_nra->add_monomial(v, vars.size(), vars.c_ptr()); +} + + void lar_solver::get_infeasibility_explanation_for_inf_sign( vector> & explanation, const vector> & inf_row, diff --git a/src/util/lp/lar_solver.h b/src/util/lp/lar_solver.h index a47240fdd..6937455ae 100644 --- a/src/util/lp/lar_solver.h +++ b/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 m_term_count; vector m_terms; const var_index m_terms_start_index; - indexed_vector m_column_buffer; + indexed_vector m_column_buffer; + scoped_ptr 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 const& vars); lp_status solve(); - void fill_explanation_from_infeasible_column(vector> & evidence) const; + void fill_explanation_from_infeasible_column(explanation_t & evidence) const; unsigned get_total_iterations() const; diff --git a/src/util/lp/lp_settings.h b/src/util/lp/lp_settings.h index c081a84ff..f07fbeef7 100644 --- a/src/util/lp/lp_settings.h +++ b/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> explanation_t; + +typedef std::unordered_map nra_model_t; + + enum class column_type { free_column = 0, low_bound = 1, diff --git a/src/util/lp/nra_solver.cpp b/src/util/lp/nra_solver.cpp index 08b17030c..2308b9c45 100644 --- a/src/util/lp/nra_solver.cpp +++ b/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 m_lp2nl; // map from lar_solver variables to nlsat::solver variables - svector m_vars; - vector> m_monomials; + struct mon_eq { + mon_eq(lean::var_index v, svector const& vs): + m_v(v), m_vs(vs) {} + lean::var_index m_v; + svector m_vs; + }; + + vector m_monomials; unsigned_vector m_lim; + mutable std::unordered_map 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(sz, vs)); + m_monomials.push_back(mon_eq(v, svector(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 core; + solver.get_core(core); + for (auto c : core) { + unsigned idx = static_cast(static_cast(c) - this); + ex.push_back(std::pair(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 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 vars; + rational den = denominator(rhs); + for (auto kv : lhs) { + vars.push_back(lp2nl(solver, kv.second)); + den = lcm(den, denominator(kv.first)); + } + vector 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); } } diff --git a/src/util/lp/nra_solver.h b/src/util/lp/nra_solver.h index bae76b0d1..87820728a 100644 --- a/src/util/lp/nra_solver.h +++ b/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.