diff --git a/src/smt/params/theory_arith_params.h b/src/smt/params/theory_arith_params.h index eb1459058..9d09b4c38 100644 --- a/src/smt/params/theory_arith_params.h +++ b/src/smt/params/theory_arith_params.h @@ -25,11 +25,11 @@ Revision History: enum arith_solver_id { AS_NO_ARITH, // 0 AS_DIFF_LOGIC, // 1 - AS_ARITH, // 2 + AS_OLD_ARITH, // 2 AS_DENSE_DIFF_LOGIC, // 3 AS_UTVPI, // 4 AS_OPTINF, // 5 - AS_LRA // 6 + AS_NEW_ARITH // 6 }; enum bound_prop_mode { @@ -113,7 +113,7 @@ struct theory_arith_params { theory_arith_params(params_ref const & p = params_ref()): m_arith_eq2ineq(false), m_arith_process_all_eqs(false), - m_arith_mode(AS_ARITH), + m_arith_mode(AS_NEW_ARITH), m_arith_auto_config_simplex(false), m_arith_blands_rule_threshold(1000), m_arith_propagate_eqs(true), diff --git a/src/smt/smt_setup.cpp b/src/smt/smt_setup.cpp index d55a50d4a..9390d91ba 100644 --- a/src/smt/smt_setup.cpp +++ b/src/smt/smt_setup.cpp @@ -730,11 +730,11 @@ namespace smt { } void setup::setup_i_arith() { - if (AS_LRA == m_params.m_arith_mode) { - setup_r_arith(); + if (AS_OLD_ARITH == m_params.m_arith_mode) { + m_context.register_plugin(alloc(smt::theory_i_arith, m_manager, m_params)); } else { - m_context.register_plugin(alloc(smt::theory_i_arith, m_manager, m_params)); + setup_r_arith(); } } @@ -747,7 +747,7 @@ namespace smt { case AS_OPTINF: m_context.register_plugin(alloc(smt::theory_inf_arith, m_manager, m_params)); break; - case AS_LRA: + case AS_NEW_ARITH: setup_r_arith(); break; default: @@ -811,12 +811,15 @@ namespace smt { case AS_OPTINF: m_context.register_plugin(alloc(smt::theory_inf_arith, m_manager, m_params)); break; - case AS_LRA: - setup_r_arith(); + case AS_OLD_ARITH: + if (m_params.m_arith_int_only && int_only) + m_context.register_plugin(alloc(smt::theory_i_arith, m_manager, m_params)); + else + m_context.register_plugin(alloc(smt::theory_mi_arith, m_manager, m_params)); break; default: if (m_params.m_arith_int_only && int_only) - m_context.register_plugin(alloc(smt::theory_i_arith, m_manager, m_params)); + setup_i_arith(); else m_context.register_plugin(alloc(smt::theory_mi_arith, m_manager, m_params)); break; diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp index de6606964..07a8aecae 100644 --- a/src/smt/theory_lra.cpp +++ b/src/smt/theory_lra.cpp @@ -307,7 +307,7 @@ class theory_lra::imp { m_solver->settings().bound_propagation() = BP_NONE != propagation_mode(); m_solver->set_track_pivoted_rows(lp.bprop_on_pivoted_rows()); m_solver->settings().m_int_gomory_cut_period = ctx().get_fparams().m_arith_branch_cut_ratio; - m_solver->settings().m_int_cuts_etc_period = ctx().get_fparams().m_arith_branch_cut_ratio; + m_solver->settings().m_hnf_cut_period = ctx().get_fparams().m_arith_branch_cut_ratio; m_solver->settings().m_int_chase_cut_solver_period = std::max(8u, ctx().get_fparams().m_arith_branch_cut_ratio); m_solver->settings().m_int_run_gcd_test = ctx().get_fparams().m_arith_gcd_test; @@ -1302,7 +1302,7 @@ public: case lp::lia_move::cut: { ++m_stats.m_gomory_cuts; // m_explanation implies term <= k - app_ref b = mk_bound(term, k, false); + app_ref b = mk_bound(term, k, !upper); m_eqs.reset(); m_core.reset(); m_params.reset(); @@ -2533,10 +2533,10 @@ public: struct scoped_arith_mode { smt_params& p; scoped_arith_mode(smt_params& p) : p(p) { - p.m_arith_mode = AS_ARITH; + p.m_arith_mode = AS_OLD_ARITH; } ~scoped_arith_mode() { - p.m_arith_mode = AS_LRA; + p.m_arith_mode = AS_NEW_ARITH; } }; @@ -2544,7 +2544,7 @@ public: if (dump_lemmas()) { ctx().display_lemma_as_smt_problem(m_core.size(), m_core.c_ptr(), m_eqs.size(), m_eqs.c_ptr(), false_literal); } - if (m_arith_params.m_arith_mode != AS_LRA) return true; + if (m_arith_params.m_arith_mode != AS_NEW_ARITH) return true; scoped_arith_mode _sa(ctx().get_fparams()); context nctx(m, ctx().get_fparams(), ctx().get_params()); add_background(nctx); @@ -2559,7 +2559,7 @@ public: if (dump_lemmas()) { ctx().display_lemma_as_smt_problem(m_core.size(), m_core.c_ptr(), m_eqs.size(), m_eqs.c_ptr(), lit); } - if (m_arith_params.m_arith_mode != AS_LRA) return true; + if (m_arith_params.m_arith_mode != AS_NEW_ARITH) return true; scoped_arith_mode _sa(ctx().get_fparams()); context nctx(m, ctx().get_fparams(), ctx().get_params()); m_core.push_back(~lit); @@ -2574,7 +2574,7 @@ public: } bool validate_eq(enode* x, enode* y) { - if (m_arith_params.m_arith_mode == AS_LRA) return true; + if (m_arith_params.m_arith_mode == AS_NEW_ARITH) return true; context nctx(m, ctx().get_fparams(), ctx().get_params()); add_background(nctx); nctx.assert_expr(m.mk_not(m.mk_eq(x->get_owner(), y->get_owner()))); diff --git a/src/test/lp/lp.cpp b/src/test/lp/lp.cpp index b40bd545a..140aebf3d 100644 --- a/src/test/lp/lp.cpp +++ b/src/test/lp/lp.cpp @@ -47,7 +47,6 @@ #include "util/lp/stacked_unordered_set.h" #include "util/lp/int_set.h" #include "util/stopwatch.h" -#include "util/lp/integer_domain.h" #include "util/lp/stacked_map.h" #include #include "test/lp/gomory_test.h" @@ -3121,152 +3120,6 @@ void get_random_interval(bool& neg_inf, bool& pos_inf, int& x, int &y) { lp_assert((neg_inf || (0 <= x && x <= 100)) && (pos_inf || (0 <= y && y <= 100))); } -void test_integer_domain_intersection(integer_domain & d) { - // int x, y; bool neg_inf, pos_inf; - // get_random_interval(neg_inf, pos_inf, x, y); - // if (neg_inf) { - // if (!pos_inf) { - // d.intersect_with_upper_bound(y); - // } - // } - // else if (pos_inf) - // d.intersect_with_lower_bound(x); - // else - // d.intersect_with_interval(x, y); -} - -void test_integer_domain_union(integer_domain & d) { - int x, y; bool neg_inf, pos_inf; - get_random_interval(neg_inf, pos_inf, x, y); - if (neg_inf) { - if (!pos_inf) { - d.unite_with_interval_neg_inf_x(y); - } - else - d.init_to_contain_all(); - } - else if (pos_inf) - d.unite_with_interval_x_pos_inf(x); - else - d.unite_with_interval(x, y); - - lp_assert(d.is_correct()); -} - - -void test_integer_domain_randomly(integer_domain & d) { - int i = my_random() % 10; - if (i == 0) - test_integer_domain_intersection(d); - else - test_integer_domain_union(d); -} - -void test_integer_domain() { -#ifdef Z3DEBUG - - std::cout << "test_integer_domain\n"; - unsigned e0 = 0; - unsigned e1 = 1; - unsigned e2 = 2; - unsigned e3 = 3; // these are explanations - unsigned e4 = 4; - unsigned e5 = 5; - integer_domain d; - unsigned l0 = 0, l1 = 1, l2 = 3; - unsigned u0 = 10, u1 = 9, u2 = 8; - d.push(); - d.intersect_with_lower_bound(l0, e0); - unsigned b; - unsigned e; - bool r = d.get_lower_bound_with_expl(b, e); - lp_assert(r && b == l0 && e == e0); - d.push(); - d.intersect_with_upper_bound(u0, e1); - r = d.get_upper_bound_with_expl(b, e); - lp_assert(r && b == u0 && e == e1); - r = d.get_lower_bound_with_expl(b, e); - lp_assert(r && b == l0 && e == e0); - d.pop(); - r = d.get_upper_bound_with_expl(b, e); - lp_assert(!r); - d.intersect_with_upper_bound(u0, e1); - d.push(); - d.intersect_with_lower_bound(l1, e2); - d.intersect_with_upper_bound(u1, e3); - d.push(); - d.intersect_with_lower_bound(l2, e4); - d.intersect_with_upper_bound(u2, e5); - lp_assert(d.is_empty() == false); - d.print(std::cout); - d.pop(); - r = d.get_lower_bound_with_expl(b, e); - lp_assert(r && b == l1 && e == e2); - d.print(std::cout); - d.pop(2); - d.print(std::cout); - lp_assert(d.has_neg_inf() && d.has_pos_inf()); -#endif - // integer_domain d; - // std::vector> stack; - // for (int i = 0; i < 10000; i++) { - // test_integer_domain_randomly(d); - // stack.push_back(d); - // d.push(); - // if (i > 0 && i%100 == 0) { - // if (stack.size() == 0) continue; - // unsigned k = my_random() % stack.size(); - // if (k == 0) - // k = 1; - // d.pop(k); - // d.restore_domain(); - // for (unsigned j = 0; j + 1 < k; j++) { - // stack.pop_back(); - // } - // std::cout<<"comparing i = " << i << std::endl; - // lp_assert(d == *stack.rbegin()); - // stack.pop_back(); - // } - // //d.print(std::cout); - // } -} - - - -void test_resolve_with_tight_constraint(chase_cut_solver& cs, - lp::chase_cut_solver::polynomial&i , - unsigned int j, - chase_cut_solver::polynomial& ti) { - - // std::cout << "resolve constraint "; - // cs.print_polynomial(std::cout, i); - // std::cout << " for " << cs.get_column_name(j) << " by using poly "; - // cs.print_polynomial(std::cout, ti); - // std::cout << std::endl; - // bool j_coeff_is_one = ti.coeff(j) == 1; - // chase_cut_solver::polynomial result; - // cs.resolve(i, j, j_coeff_is_one, ti); - // std::cout << "resolve result is "; - // cs.print_polynomial(std::cout, i); - // std::cout << std::endl; -} - -typedef chase_cut_solver::monomial mono; - -void test_resolve(chase_cut_solver& cs, unsigned constraint_index, unsigned i0) { - var_index x = 0; - var_index y = 1; - var_index z = 2; - std::cout << "test_resolve\n"; - - chase_cut_solver::polynomial i; i += mono(2, x);i += mono(-3,y); - i+= mono(4, z); - i.m_a = 5; - chase_cut_solver::polynomial ti; ti += mono(1, x); ti+= mono(1,y);ti.m_a = 3; - test_resolve_with_tight_constraint(cs, i, x, ti); - test_resolve_with_tight_constraint(cs, i, y ,ti); -} - void test_gomory_cut_0() { gomory_test g( @@ -3688,15 +3541,29 @@ void test_larger_generated_hnf() { std::cout << "test_larger_generated_rank_hnf" << std::endl; general_matrix A; vector v; - v.push_back(zero_of_type()); - v.push_back(zero_of_type()); - v.push_back(zero_of_type()); - A.push_row(v); + v.clear(); + v.push_back(mpq(5)); + v.push_back(mpq(6)); + v.push_back(mpq(3)); + v.push_back(mpq(1)); A.push_row(v); v.clear(); v.push_back(mpq(5)); - v.push_back(mpq(26)); + v.push_back(mpq(2)); v.push_back(mpq(3)); + v.push_back(mpq(7)); + A.push_row(v); + v.clear(); + v.push_back(mpq(5)); + v.push_back(mpq(6)); + v.push_back(mpq(3)); + v.push_back(mpq(1)); + A.push_row(v); + v.clear(); + v.push_back(mpq(5)); + v.push_back(mpq(2)); + v.push_back(mpq(3)); + v.push_back(mpq(7)); A.push_row(v); call_hnf(A); std::cout << "test_larger_generated_rank_hnf passed" << std::endl; @@ -3750,10 +3617,6 @@ void test_lp_local(int argn, char**argv) { test_gomory_cut(); return finalize(0); } - if (args_parser.option_is_used("-intd")) { - test_integer_domain(); - return finalize(0); - } if (args_parser.option_is_used("--test_mpq")) { test_rationals(); diff --git a/src/util/lp/CMakeLists.txt b/src/util/lp/CMakeLists.txt index a946328f6..c67c55759 100644 --- a/src/util/lp/CMakeLists.txt +++ b/src/util/lp/CMakeLists.txt @@ -4,7 +4,6 @@ z3_add_component(lp binary_heap_priority_queue.cpp binary_heap_upair_queue.cpp bound_propagator.cpp - chase_cut_solver.cpp core_solver_pretty_printer.cpp dense_matrix.cpp eta_matrix.cpp diff --git a/src/util/lp/chase_cut_solver.cpp b/src/util/lp/chase_cut_solver.cpp deleted file mode 100644 index 1369e0710..000000000 --- a/src/util/lp/chase_cut_solver.cpp +++ /dev/null @@ -1,23 +0,0 @@ -/* - Copyright (c) 2017 Microsoft Corporation - Author: Nikolaj Bjorner, Lev Nachmanson -*/ -#include "util/lp/chase_cut_solver.h" -namespace lp { - mpq polynomial::m_local_zero = zero_of_type(); - - size_t constraint_hash::operator() (const constraint* c) const { return c->id(); } - - bool constraint_equal::operator() (const constraint* a, const constraint * b) const { return a->id() == b->id(); } - - std::ostream& operator<<(std::ostream& out, pp_poly const& p) { - p.s.print_polynomial(out, p.p); - return out; - } - - std::ostream& operator<<(std::ostream& out, pp_constraint const& c) { - c.s.print_constraint(out, c.c); - return out; - } -} - diff --git a/src/util/lp/chase_cut_solver.h b/src/util/lp/chase_cut_solver.h deleted file mode 100644 index 4b41468fe..000000000 --- a/src/util/lp/chase_cut_solver.h +++ /dev/null @@ -1,2534 +0,0 @@ -/*++ - Copyright (c) 2017 Microsoft Corporation - - Module Name: - - - - Abstract: - - - - Author: - Nikolaj Bjorner (nbjorner) - Lev Nachmanson (levnach) - - Revision History: - - - --*/ - -#pragma once -#include "util/vector.h" -#include "util/trace.h" -#include "util/lp/lp_settings.h" -#include "util/lp/column_namer.h" -#include "util/lp/integer_domain.h" -#include "util/lp/lp_utils.h" -#include -#include "util/lp/int_set.h" -#include "util/lp/stacked_vector.h" -#include "util/lp/monomial.h" -#include "util/lp/polynomial.h" -#include "util/lp/constraint.h" -#include "util/lp/active_set.h" -#include "util/lp/indexer_of_constraints.h" -#include "util/lp/lar_term.h" -namespace lp { -class chase_cut_solver; //forward definition - -struct pp_poly { - chase_cut_solver const& s; - polynomial const& p; - pp_poly(chase_cut_solver const& s, polynomial const& p): s(s), p(p) {} -}; - -struct pp_constraint { - chase_cut_solver const& s; - constraint const& c; - pp_constraint(chase_cut_solver const& s, constraint const& c): s(s), c(c) {} -}; - -std::ostream& operator<<(std::ostream& out, pp_poly const& p); -std::ostream& operator<<(std::ostream& out, pp_constraint const& p); - -class chase_cut_solver : public column_namer { -public: - enum class lbool { l_false, l_true, l_undef }; - inline std::string lbool_to_string(lbool l) { - switch(l) { - case lbool::l_true: return std::string("true"); - case lbool::l_false: return std::string("false"); - case lbool::l_undef: return std::string("undef"); - default: - return std::string("what is it?"); - } - } - - typedef lp::polynomial polynomial; - typedef lp::monomial monomial; - - vector> to_pairs(const vector& ms) const { - vector> ret; - for (const auto p : ms) - ret.push_back(p.to_pair()); - return ret; - } - - - typedef const constraint const_constr; - - class literal { - int m_decision_context_index; // points to the trail element if a decision has been made - unsigned m_var; - bool m_is_lower; - mpq m_bound; - constraint* m_constraint; // nullptr if it is a decided literal - constraint* m_tight_constr; // nullptr if it is not calculated - unsigned m_prev_var_level; - public: - private: - literal( // creates an implied bound - unsigned var_index, - bool is_lower, - const mpq & bound, - constraint * constr, - unsigned prev_var_level) : - m_decision_context_index(-1), - m_var(var_index), - m_is_lower(is_lower), - m_bound(bound), - m_constraint(constr), - m_tight_constr(nullptr), - m_prev_var_level(prev_var_level) { - } - literal( // creates a decided bound - int trail_index, - unsigned var_index, - bool is_lower, - const mpq & bound, - unsigned prev_var_level): - m_decision_context_index(trail_index), - m_var(var_index), - m_is_lower(is_lower), - m_bound(bound), - m_constraint(nullptr), - m_tight_constr(nullptr), - m_prev_var_level(prev_var_level) { - } - public: - const constraint * tight_constr() const { return m_tight_constr; } - constraint*& tight_constr () { return m_tight_constr; } - const mpq & bound() const { return m_bound; } - bool is_lower() const { return m_is_lower; } - bool is_upper() const { return !m_is_lower; } - int decision_context_index() const { return m_decision_context_index; } - const_constr * cnstr() const { return m_constraint; } - constraint * cnstr() { return m_constraint; } - literal() {} - - bool tight_constraint_is_calculated() const { - return m_tight_constr != nullptr; - } - - bool is_decided() const { return m_decision_context_index != -1; } - - bool is_implied() const { return !is_decided();} - - // TBD: would be nice with a designated type for variables? - - unsigned var() const { return m_var; } - - static literal make_implied_literal(unsigned var_index, bool is_lower, const mpq & bound, constraint * c, unsigned var_level) { - return literal(var_index, is_lower, bound, c, var_level); - } - static literal make_decided_literal(unsigned var_index, bool is_lower, - const mpq & bound, int decision_context_index, unsigned var_level) { - return literal(decision_context_index, var_index, is_lower, bound, var_level); - } - - unsigned prev_var_level() const { return m_prev_var_level; } - }; - - enum class bound_type { - LOWER, UPPER, UNDEF - }; - - struct bound_result { - mpq m_bound; - bound_type m_type; - - bound_result(const mpq & b, bound_type bt): m_bound(b), m_type(bt) {} - bound_result() : m_type(bound_type::UNDEF) { - } - void print( std::ostream & out) const { - if (m_type == bound_type::LOWER) { - out << "lower bound = "; - } - else if (m_type == bound_type::UPPER) { - out << "upper bound = "; - } - else { - out << "undef"; - return; - } - out << m_bound; - } - const mpq & bound() const { return m_bound; } - }; - - class var_info { - unsigned m_internal_j; // it is just the index into m_var_infos of this var_info, if the var_info is not active then this value is set to (unsigned)-1 - integer_domain m_domain; - // the map of constraints using the var: bound_type = UNDEF stands for a div constraint - std::unordered_map m_dependent_constraints; // pair (c, true) means that c has a monomial (a, m_internal_j) for a > 0, and (c, false) means that a < 0 - unsigned m_external_stack_level; - unsigned m_number_of_bound_propagations; - unsigned m_number_of_asserts; - - public: - unsigned number_of_asserts() const { return m_number_of_asserts; } - - var_info(unsigned user_var_index) : - m_internal_j(user_var_index), - m_external_stack_level(static_cast(-1)), - m_number_of_bound_propagations(0), - m_number_of_asserts(0) - {} - var_info() : m_internal_j(static_cast(-1)), - m_external_stack_level(static_cast(-1)), - m_number_of_bound_propagations(0), - m_number_of_asserts(0) {} - - bool is_active() const { return m_internal_j != static_cast(-1); } - - const integer_domain & domain() const { return m_domain; } - unsigned internal_j() const { - return m_internal_j; - } - void activate(unsigned internal_j) { - m_internal_j = internal_j; - } - void add_dependent_constraint(constraint* i, bool coeff_is_pos) { - lp_assert(m_dependent_constraints.find(i) == m_dependent_constraints.end()); - if (i->is_assert()) - m_number_of_asserts++; - m_dependent_constraints[i] = coeff_is_pos; - } - void remove_depended_constraint(constraint* i) { - lp_assert(m_dependent_constraints.find(i) != m_dependent_constraints.end()); - if (i->is_assert()) - m_number_of_asserts--; - m_dependent_constraints.erase(i); - } - - void conditional_push(unsigned external_level) { - if (external_level != m_external_stack_level) { - m_domain.push(); - m_external_stack_level = external_level; - } - } - - bool intersect_with_lower_bound(const mpq & b, unsigned explanation, unsigned stack_level) { - conditional_push(stack_level); - bool ret = m_domain.intersect_with_bound(b, true, explanation); - lp_assert(!m_domain.is_empty()); - return ret; - } - - bool intersect_with_upper_bound(const mpq & b, unsigned explanation, unsigned external_level) { - conditional_push(external_level); - bool ret = m_domain.intersect_with_bound(b, false, explanation); - lp_assert(!m_domain.is_empty()); - return ret; - } - - bool is_fixed() const { return m_domain.is_fixed();} - bool get_upper_bound(mpq & b) const { return m_domain.get_upper_bound(b); } - bool get_lower_bound(mpq & b) const { return m_domain.get_lower_bound(b); } - bool get_upper_bound_with_expl(mpq & b, unsigned& expl) const { return m_domain.get_upper_bound_with_expl(b, expl); } - bool get_lower_bound_with_expl(mpq & b, unsigned& expl) const { return m_domain.get_lower_bound_with_expl(b, expl); } - void print_var_domain(std::ostream &out) const { m_domain.print(out); } - std::unordered_map & dependent_constraints() { return m_dependent_constraints; } - const std::unordered_map & dependent_constraints() const { return m_dependent_constraints; } - int get_lower_bound_expl() const { return m_domain.get_lower_bound_expl();} - int get_upper_bound_expl() const { return m_domain.get_upper_bound_expl();} - public: - void pop(unsigned k, unsigned external_level) { - m_domain.pop(k); - m_external_stack_level = external_level; - } - unsigned external_stack_level() const { return m_external_stack_level; } - unsigned &external_stack_level() { return m_external_stack_level; } - - unsigned number_of_bound_propagations() const { return m_number_of_bound_propagations; } - unsigned & number_of_bound_propagations() { return m_number_of_bound_propagations; } - - }; // end of var_info - - struct lemmas_container { - std::unordered_set m_lemmas; - unsigned size() const { return m_lemmas.size(); } - void submit_assert_origin_for_delete(constraint_index ci) { m_deleted_assert_origins.insert(ci);} - vector remove_lemmas_depending_on_submitted_origins() { - vector r; - for (auto * l : m_lemmas) { - for (unsigned o : l->assert_origins()) - if (m_deleted_assert_origins.find(o) != m_deleted_assert_origins.end()) { - r.push_back(l); - break; // from the inner loop - } - } - for (auto * l : r) - m_lemmas.erase(l); - - m_deleted_assert_origins.clear(); - return r; - } - std::unordered_set m_deleted_assert_origins; - void add_lemma(constraint* l) { m_lemmas.insert(l); } - }; - - bool lhs_is_int(const vector & lhs) const { - for (auto & p : lhs) { - if (numeric_traits::is_int(p.coeff()) == false) return false; - } - return true; - } - -public: - - bool all_fixed_except_j(const polynomial & p, var_index j) const { - for (auto &m : p.coeffs()) - if (m.var() != j && m_var_infos[m.var()].is_fixed() == false) - return false; - return true; - } - - bool lower_bound_exists(const var_info & v) const { - return v.get_lower_bound_expl() != -1; - } - - bool upper_bound_exists(const var_info & v) const { - return v.get_upper_bound_expl() != -1; - } - - bool lower_bound_exists(const integer_domain & v) const { - return v.get_lower_bound_expl() != -1; - } - - bool upper_bound_exists(const integer_domain & v) const { - return v.get_upper_bound_expl() != -1; - } - - - bool is_cc(var_index j, const constraint*&lower, const constraint*&upper) const { - const var_info & vj = m_var_infos[j]; - if (lower_bound_exists(vj) && upper_bound_exists(vj)) - return false; - if (vj.domain().is_empty()) - return false; - - unsigned upper_bounds = 0; - unsigned lower_bounds = 0; - for (auto p : vj.dependent_constraints()) { - constraint* c = p.first; - const mpq& coeff = c->poly().coeff(j); - if (coeff == one_of_type() || coeff == - one_of_type()) - continue; - if (!all_fixed_except_j(c->poly(), j)) continue; - if (!p.second) { - upper_bounds++; - upper = c; - if (lower_bounds) return true; - } else { - lower_bounds++; - lower = c; - if (upper_bounds) - return true; - } - } - return false; - } - - std::string get_column_name(unsigned j) const { - return m_var_name_function(m_var_infos[j].internal_j()); - } - - - ~chase_cut_solver() { - for (constraint * c : m_asserts) - delete c; - for (constraint * c : m_lemmas_container.m_lemmas) - delete c; - } - - struct scope { - unsigned m_trail_size; - scope() {} - scope( unsigned trail_size) : - m_trail_size(trail_size) {} - }; - - // fields - vector m_var_infos; - vector m_asserts; - lemmas_container m_lemmas_container; - vector m_v; // the values of the variables - std::function m_var_name_function; - std::function m_print_constraint_function; - std::function m_number_of_variables_function; - std::function m_var_value_function; - active_set m_active_set; - vector m_trail; - lp_settings & m_settings; - std::set m_U; // the set of conflicting cores - unsigned m_bounded_search_calls; - unsigned m_number_of_conflicts; - vector m_scopes; - std::unordered_map m_user_vars_to_chase_cut_solver_vars; - std::unordered_set m_explanation; // if this collection is not empty we have a conflict - // the number of decisions in the current trail - unsigned m_decision_level; - bool m_stuck_state; - bool m_cancelled; - // debug fields - unsigned m_number_of_constraints_tried_for_propagaton; - unsigned m_pushes_to_trail; - indexer_of_constraints m_indexer_of_constraints; - - - bool is_lower_bound(literal & l) const { - return l.is_lower(); - } - - // bool lower_for_var(unsigned j, mpq & lower) const { - // bool ret = false; - // for (unsigned i : m_var_infos[j].m_literals) - // if (is_lower_bound(m_trail[i])) { - // if (ret == false) { - // ret = true; - // lower = get_bound(m_trail[i]); - // } else { - // lower = std::max(lower, get_bound(m_trail[i])); - // } - // } - // return ret; - // } - - bool is_upper_bound(literal & l) const { - return !l.is_lower(); - } - - bool at_base_lvl() const { return m_decision_level == 0; } - - void simplify_problem() { - // no-op - } - - void gc() { - // no-op - } - - void restart() { - // no-op for now - } - - bool check_inconsistent() { - if (at_base_lvl()) { - for (constraint *c : m_lemmas_container.m_lemmas ) { - if (lower_is_pos(c)) { - TRACE("check_inconsistent_int", tout << pp_poly(*this, c->poly()) << "\n";); - lp_assert(false); // not implemented - return true; - } - } - for (unsigned j = m_asserts.size(); j--; ) { - if (lower_is_pos(m_asserts[j])) { - TRACE("check_inconsistent_int", tout << pp_poly(*this, m_asserts[j]->poly()) << "\n";); - lp_assert(false); // not implemented - return true; - } - } - } - return false; - } - - void cleanup() { } - - bool all_constraints_hold() const { - for (auto c : m_asserts) { - if (is_pos(value(c->poly()))) { - TRACE("all_constraints_hold_int", tout << "constraint does not hold\n"; - tout << pp_constraint(*this, *c) << "value = " << value(c->poly()) << std::endl;); - - return false; - } - } - for (auto c : m_lemmas_container.m_lemmas) { - if (is_pos(value(c->poly()))) { - TRACE("all_constraints_hold_int", tout << "constraint does not hold\n"; - print_constraint(tout, *c);); - return false; - } - } - return true; - } - - lbool final_check() { - if (!all_vars_are_fixed()) { - TRACE("final_check_int", tout << "return undef" << std::endl;); - m_stuck_state = true; - return lbool::l_undef; // we are stuck - } - lp_assert(all_constraints_hold()); - lp_assert(at_base_lvl()); - // there are no more case splits, and all clauses are satisfied. - // prepare the model for external consumption. - return lbool::l_true; - } - - void find_new_conflicting_cores_under_constraint(var_index j, const_constr* c) { - // lp_assert(false); - } - - void find_new_conflicting_cores(var_index j) { - for (auto p: m_var_infos[j].dependent_constraints()) - find_new_conflicting_cores_under_constraint(j, p.first); - } - - void set_value_for_fixed_var_and_check_for_conf_cores(unsigned j) { - TRACE("set_value_for_fixed_var_and_check_for_conf_cores", tout << "j = " << j << std::endl;); - if (m_v.size() <= j) - m_v.resize(j + 1); - lp_assert(m_var_infos[j].is_fixed()); - lp_assert(m_var_infos[j].is_active()); - m_var_infos[j].domain().get_upper_bound(m_v[j]); // sets the value of the variable - find_new_conflicting_cores(j); - } - - void restrict_var_domain_with_bound_result(var_index j, const bound_result & br, unsigned trail_index) { - TRACE("restrict_var_domain_with_bound_result", tout << "j = " << j << std::endl; - br.print(tout);); - auto & vi = m_var_infos[j]; - lp_assert(!vi.is_fixed()); - lp_assert(m_trail.back().var() == j); - if (br.m_type == bound_type::UPPER) { - vi.intersect_with_upper_bound(br.m_bound, trail_index, m_scopes.size()); - } else { - vi.intersect_with_lower_bound(br.m_bound, trail_index, m_scopes.size()); - } - if (vi.is_fixed()) { - TRACE("d.is_fixed", tout << "j = " << j << std::endl;); - set_value_for_fixed_var_and_check_for_conf_cores(j); - } - } - - // 'j' is a variable that changed - void add_changed_var(unsigned j, bool lower_bound_got_tighter, int priority) { - TRACE("add_changed_var_int", tout << "j = " << j << "\n";); - for (auto & p: m_var_infos[j].dependent_constraints()) { - TRACE("add_changed_var_int", tout << pp_constraint(*this, *p.first) << "\n";); - if (p.second == lower_bound_got_tighter) - m_active_set.add_constraint(p.first, 2); // 2 is a priority - } - } - - unsigned number_of_vars() const { return static_cast( m_var_infos.size()); } - - const mpq & var_value(unsigned j) const { - return m_v[j]; - } - - bool var_is_active(unsigned j) const { - return m_var_infos[j].is_active(); - } - - struct test_bound_struct { - mpq m_lower_bound; - mpq m_upper_bound; - int m_expl_lower; - int m_expl_upper; - test_bound_struct() :m_expl_lower(-1), m_expl_upper(-1) {} - }; - - bool state_is_correct() const { - return var_infos_are_correct() && constraint_indices_are_correct(); - } - - bool constraint_indices_are_correct() const { - std::unordered_set r; - unsigned n = m_indexer_of_constraints.max(); - for (auto c: m_asserts) { - if (c->id() >= n) { - std::cout << "c->id >= n\n"; - return false; - } - if (r.find(c->id()) != r.end()){ - std::cout << "id found\n"; - return false; - } - r.insert(c->id()); - } - - for (auto c: m_lemmas_container.m_lemmas) { - if (c->id() >= n) { - std::cout << "c->id >= n\n"; - return false; - } - if (r.find(c->id()) != r.end()){ - std::cout << "id found\n"; - return false; - } - r.insert(c->id()); - } - return true; - } - - bool var_infos_are_correct() const { - vector bounds = get_bounds_from_trail(); - for (unsigned j = 0; j < m_var_infos.size(); j++) - if (!var_info_is_correct(j, bounds[j])) { - TRACE("var_infos_are_correct", print_var_info(tout, j); tout << " var_info is incorrect j = " << j;); - return false; - } - return true; - } - - bound_result bound_test_on_poly(const polynomial &p, const mpq& a, unsigned j, vector& bs) const { - lp_assert(!is_zero(a)); - if (numeric_traits::is_pos(a)) { - TRACE("bound_test_on_poly", tout << var_name(j) << " a is pos, p = "; print_polynomial(tout, p);); - bound_result r = lower_without_test(p, j, bs); - if (r.m_type == bound_type::UNDEF) - return r; - lp_assert(is_int(r.m_bound)); - r.m_bound = - ceil_ratio(r.m_bound , a); - r.m_type = bound_type::UPPER; - return r; - } - else { - TRACE("bound_test_on_poly", tout << var_name(j) << " a is neg, p = "; print_polynomial(tout, p);); - bound_result r = lower_without_test(p, j, bs); - if (r.m_type == bound_type::UNDEF) - return r; - r.m_bound = -floor_ratio(r.m_bound, a); - r.m_type = bound_type::LOWER; - return r; - } - - } - - bound_result bound_test(unsigned j, const polynomial &p, vector& bs) const { - const mpq &a = p.coeff(j); - bound_result br = bound_test_on_poly(p, a, j, bs); - TRACE("bound_test", tout << var_name(j) << ", p = "; print_polynomial(tout, p); - tout << ", br = "; br.print(tout); ); - return br; - } - - void bound_test_literal(const literal & l, vector & bs) const { - bound_result br = l.cnstr() != nullptr? - bound_test(l.var(), l.cnstr()->poly(), bs) : - bound_result(l.bound(), l.is_lower() ? bound_type::LOWER: bound_type::UPPER); - lp_assert(br.m_type != bound_type::UNDEF); - if (l.is_lower()) { - lp_assert(br.m_type == bound_type::LOWER); - lp_assert(br.m_bound >= l.bound()); - } else { - lp_assert(br.m_type == bound_type::UPPER); - lp_assert(br.m_bound <= l.bound()); - } - - if (l.tight_constr() != nullptr) { - br = bound_test(l.var(), l.tight_constr()->poly(), bs); - lp_assert(br.m_type != bound_type::UNDEF); - if (l.is_lower()) { - lp_assert(br.m_type == bound_type::LOWER); - lp_assert(br.m_bound >= l.bound()); - } else { - lp_assert(br.m_type == bound_type::UPPER); - lp_assert(br.m_bound <= l.bound()); - } - } - - } - - void get_bounds_from_trail_elem(unsigned j, vector& bs) const { - const literal & l = m_trail[j]; - auto & t = bs[l.var()]; - TRACE("get_bounds_from_trail_elem", tout << "j = " << j << ", literal = "; print_literal(tout, l);); - bound_test_literal(l, bs); - if (l.is_lower()) { - if (t.m_expl_lower == -1) { - t.m_expl_lower = j; - t.m_lower_bound = l.bound(); - TRACE("get_bounds_from_trail_elem", tout << var_name(l.var()) << " m_lower_bound = " << t.m_lower_bound << std::endl;); - } else { - lp_assert(t.m_lower_bound < l.bound()); - t.m_lower_bound = l.bound(); - TRACE("get_bounds_from_trail_elem", tout << var_name(l.var()) << " m_lower_bound = " << t.m_lower_bound << std::endl;); - t.m_expl_lower = j; - } - } else { - if (t.m_expl_upper == -1) { - t.m_expl_upper = j; - t.m_upper_bound = l.bound(); - TRACE("get_bounds_from_trail_elem", tout << var_name(l.var()) << " m_u = " << t.m_upper_bound << std::endl;); - } else { - lp_assert(t.m_upper_bound > l.bound()); - t.m_upper_bound = l.bound(); - t.m_expl_upper = j; - TRACE("get_bounds_from_trail_elem", tout << var_name(l.var()) << " m_u = " << t.m_upper_bound << std::endl;); - } - } - } - - vector get_bounds_from_trail() const { - vector ret(m_var_infos.size()); - for (unsigned j = 0; j < m_trail.size(); j++) { - get_bounds_from_trail_elem(j, ret); - } - return ret; - } - - bool var_info_is_correct(unsigned j, const test_bound_struct& t) const { - const var_info & v = m_var_infos[j]; - if (v.external_stack_level() != static_cast(-1) && v.external_stack_level() > m_scopes.size()) { - - TRACE("var_info_is_correct", tout << "incorrect: the level is too high:"; - print_var_info(tout, j); - tout << "external_level = " << v.external_stack_level(); - tout << "\nm_scopes.size() = " << m_scopes.size(); ); - return false; - } - std::unordered_set deps; - for (const auto c: m_asserts) { - if (!is_zero(c->coeff(j))) - deps.insert(c); - } - for (const auto c: m_lemmas_container.m_lemmas) { - if (!is_zero(c->coeff(j))) - deps.insert(c); - } - - for (auto p : v.dependent_constraints()) { - if (deps.find(p.first) == deps.end()) { - TRACE("var_info_is_correct", tout << "deps.find(p.first) == deps.end()";); - return false; - } - } - if (v.is_fixed() && m_v.size() <= j) { - TRACE("var_info_is_correct", tout << "m_v is too short, j="<< j <= m_trail.size()) { - TRACE("var_bound_is_correct_by_trail", tout<< "expl =" << expl << " of " << var_name(j) << ", j = "<< j << " points out of the trail, trail.size() = " << m_trail.size(); tout << "return false";); - return false; - } - const literal & l = m_trail[expl]; - - if (! (b == l.bound() && !l.is_lower() && j == l.var())) { - TRACE("var_bound_is_correct_by_trail", tout<< "expl=" << expl << std::endl; print_var_info(tout, j); tout << "b=" << b<<", expl = " << expl << std::endl; print_literal(tout, l); tout << "return false";); - return false; - } - return (t.m_expl_upper == static_cast(expl) && t.m_upper_bound == b); - } - CTRACE("var_bound_is_correct_by_trail", t.m_expl_upper != -1 , - tout << "literal index = " << t.m_expl_upper << "\n"; - print_literal(tout, m_trail[t.m_expl_upper]); - tout << "will return false"; - ); - - return t.m_expl_upper == -1; - } - - // bool run_over_trail_with_upper_bound(unsigned j, const mpq & b) const { - // bool limited = false; - // mpq found_b; - // run_over_trail_with_upper_bound_on_literal(j, b, - // } - - bool var_lower_bound_is_correct_by_trail(unsigned j, const test_bound_struct& t) const { - const var_info & v = m_var_infos[j]; - mpq b; - unsigned expl; - if (v.get_lower_bound_with_expl(b, expl)) { - if (expl >= m_trail.size()) { - TRACE("var_bound_is_correct_by_trail", tout<< "expl =" << expl << " of " << var_name(j) << ", j = "<< j << " points out of the trail, trail.size() = " << m_trail.size(); - tout << ", "; print_var_info(tout, j);); - return false; - } - const literal & l = m_trail[expl]; - - - if (! (t.m_expl_lower == static_cast(expl) && t.m_lower_bound == b)) { - TRACE("var_bound_is_correct_by_trail", print_var_info(tout, j); tout << "b=" << b<<", expl = " << expl << std::endl; print_literal(tout, l); tout << "return " << (b == l.bound() && l.is_lower());); - return false; - } - return b == l.bound() && l.is_lower(); - } - - CTRACE("var_bound_is_correct_by_trail", !(t.m_expl_lower == -1), tout << "return " << (t.m_expl_lower == -1) << " t.m_expl_lower = " << t.m_expl_lower << ", "; print_var_info(tout, v); print_trail(tout);); - return t.m_expl_lower == -1; - } - - void push_literal_to_trail(literal & l) { - m_pushes_to_trail ++; - m_trail.push_back(l); - add_changed_var(l.var(), l.is_lower(), 0); - } - - unsigned find_large_enough_j(unsigned i) { - unsigned r = 0; - for (const auto & p : m_asserts[i]->poly().m_coeffs) { - r = std::max(r, p.var() + 1); - } - return r; - } - - std::string var_name(unsigned j) const { - return get_column_name(j); - } - - void trace_print_domain_change(std::ostream& out, unsigned j, const mpq& v, const monomial & p, const_constr* c) const { - out << "trail.size() = " << m_trail.size() << "\n"; - out << "change in domain of " << var_name(j) << ", v = " << v << ", domain becomes "; - print_var_domain(out, j); - mpq lb; - bool r = lower(c, lb); - if (r) - out << "lower_of_constraint = " << lb << "\n"; - else - out << "no lower bound for constraint\n"; - } - - - - bool new_lower_bound_is_relevant(unsigned j, const mpq & v) const { - mpq lb; - const var_info & vi = m_var_infos[j]; - auto & d = vi.domain(); - bool has_bound = d.get_lower_bound(lb); - if (!has_bound) - return true; - if (v <= lb) - return false; - if (upper_bound_exists(d)) - return true; - - if (too_many_propagations_for_var(vi)) - return false; - - // int delta = 2; - return v >= lb + 2 * abs(lb); - - } - - bool too_many_propagations_for_var(const var_info& vi) const { - return vi.number_of_bound_propagations() > m_settings.m_chase_cut_solver_cycle_on_var * vi.number_of_asserts(); - } - - bool new_upper_bound_is_relevant(unsigned j, const mpq & v) const { - auto & vi = m_var_infos[j]; - auto & d = vi.domain(); - mpq b; - bool has_bound = d.get_upper_bound(b); - if (!has_bound) - return true; - if (v >= b) - return false; - if (lower_bound_exists(d)) - return true; - - if (too_many_propagations_for_var(vi)) - return false; - - // delta = 2 - return v <= b - 2 * abs(b); // returns false if the improvement is small - } - - void intersect_var_info_with_upper_bound(unsigned j, const mpq & v) { - lp_assert(m_trail.back().var() == j && m_trail.back().is_upper()); - m_var_infos[j].intersect_with_upper_bound(v, m_trail.size() - 1, m_scopes.size()); // m_trail.size() - 1 is the explanation - if (m_var_infos[j].is_fixed()) - set_value_for_fixed_var_and_check_for_conf_cores(j); - } - - void intersect_var_info_with_lower_bound(unsigned j, const mpq& v) { - lp_assert(m_trail.back().var() == j && m_trail.back().is_lower()); - m_var_infos[j].intersect_with_lower_bound(v, m_trail.size() - 1, m_scopes.size()); // m_trail.size() - 1 is the explanation - if (m_var_infos[j].is_fixed()) - set_value_for_fixed_var_and_check_for_conf_cores(j); - } - - void propagate_monomial_on_lower(const monomial & p, - const mpq& lower_val, - constraint* c) { - unsigned j = p.var(); - if (m_var_infos[j].is_fixed()) - return; - if (is_pos(p.coeff())) { - mpq m; - get_var_lower_bound(p.var(), m); - mpq v = floor(- lower_val / p.coeff()) + m; - if (new_upper_bound_is_relevant(j, v)) { - add_bound(v, j, false, c); - intersect_var_info_with_upper_bound(j, v); - } - } else { - mpq m; - get_var_upper_bound(p.var(), m); - mpq v = ceil( - lower_val / p.coeff()) + m; - if (new_lower_bound_is_relevant(j, v)) { - add_bound(v, j, true, c); - intersect_var_info_with_lower_bound(j, v); - } - } - } - - void propagate_monomial_on_right_side(const monomial & p, - const mpq& rs, - constraint *c) { - unsigned j = p.var(); - - if (is_pos(p.coeff())) { - mpq v = floor(rs / p.coeff()); - if (new_upper_bound_is_relevant(j, v)) { - add_bound(v, j, false, c); - intersect_var_info_with_upper_bound(j, v); - TRACE("ba_int_change", trace_print_domain_change(tout, j, v, p, c);); - } - } else { - mpq v = ceil(rs / p.coeff()); - if (new_lower_bound_is_relevant(j, v)) { - TRACE("ba_int_change", print_var_info(tout, j);); - add_bound(v, j, true , c); - intersect_var_info_with_lower_bound(j, v); - TRACE("ba_int_change", trace_print_domain_change(tout, j, v, p, c);); - } - } - } - - void print_var_info(std::ostream & out, const var_info & vi) const { - out << m_var_name_function(vi.internal_j()) << " "; - if (vi.internal_j() < m_v.size()) { - out << "m_v[" << vi.internal_j() << "] = " << m_v[vi.internal_j()] << std::endl; - } - if (vi.is_active()) out << ", active var "; - print_var_domain(out, vi); - out << ", propagations = " << vi.number_of_bound_propagations() << ", deps = " << vi.dependent_constraints().size(); - out << ", asserts = " << vi.number_of_asserts() << std::endl; - // out << "external levels: "; - // for (auto j : vi.external_stack_level()) - // out << j << " "; - } - - void print_var_info(std::ostream & out, unsigned j) const { - print_var_info(out, m_var_infos[j]); - } - - void print_var_domain(std::ostream & out, unsigned j) const { - m_var_infos[j].print_var_domain(out); - } - - void print_var_domain(std::ostream & out, const var_info & vi) const { - vi.print_var_domain(out); - } - - // b is the value of lower - void propagate_constraint_on_lower(constraint* c, const mpq & b) { - for (const auto & p: c->coeffs()) - propagate_monomial_on_lower(p, b, c); - } - - - void explain_literal(unsigned trail_index, - std::unordered_set & visited_literals, - std::unordered_set & visited_constraints) { - if (visited_literals.find(trail_index) != visited_literals.end()) - return; - visited_literals.insert(trail_index); - const literal & l = m_trail[trail_index]; - const_constr* c = l.tight_constr(); - if (c == nullptr) - c = l.cnstr(); - lp_assert(c != nullptr); - add_premises(c, visited_constraints); - explain_bound(c, trail_index, visited_literals, visited_constraints); - } - - void explain_bound(const_constr * c, - unsigned trail_lim, - std::unordered_set & visited_literals, - std::unordered_set visited_constraints) { - add_premises(c, visited_constraints); - TRACE("fill_conflict_explanation", tout << "visited_constraints\n"; - for (auto cc: visited_constraints) - trace_print_constraint(tout, *cc);); - for (const monomial & m : c->poly().coeffs()) { - int trail_index = find_literal_index_after(m.var(), is_pos(m.coeff()), trail_lim); - if (trail_index == -1) - continue; - explain_literal(trail_index, visited_literals, visited_constraints); - } - - } - - void dumb_exlplanations() { - m_explanation.clear(); - for (auto c: m_asserts) - for (auto j : c->assert_origins()) - m_explanation.insert(j); - } - - void add_premises(const_constr *c, std::unordered_set & visited_constraints) { - if (visited_constraints.find(c) != visited_constraints.end()) - return; - - visited_constraints.insert(c); - - for (auto j : c->assert_origins()) - m_explanation.insert(j); - } - - void fill_conflict_explanation(const constraint *confl) { - //dumb_exlplanations(); - TRACE("fill_conflict_explanation", - trace_print_constraint(tout, *confl); - print_trail(tout); - ); - m_explanation.clear(); - std::unordered_set visited_literals; - std::unordered_set visited_constraints; - explain_bound(confl, m_trail.size(), visited_literals, visited_constraints); - TRACE("fill_conflict_explanation", tout << "m_explanation = "; - for (unsigned j : m_explanation) { - tout << j << " "; - }); - } - - void trace_print_constraint(std::ostream& out, const_constr& i) const { - trace_print_constraint(out, &i); - } - - void trace_print_constraint(std::ostream& out, const_constr* i) const { - print_constraint(out, *i); - out << "id = " << i->id() << ", "; - for (auto p : i->poly().m_coeffs){ - out << "domain of " << var_name(p.var()) << " = "; - print_var_domain(out, p.var()); - } - if (i->assert_origins().size()) { - out << (i->assert_origins().size() > 1?"origins: ":"origin: ") ; - for (auto o : i->assert_origins()) - out << o << ", "; - out << "\n"; - } - } - - - void propagate_constraint_only_one_unlim(constraint* i, unsigned the_only_unlim) { - mpq rs = - i->poly().m_a; - for (unsigned j = 0; j < i->poly().m_coeffs.size(); j++) { - if (j == the_only_unlim) continue; - mpq m; - lower_monomial(i->poly().m_coeffs[j], m); - rs -= m; - } - - // we cannot get a conflict here because the monomial i.poly().m_coeffs[the_only_unlim] - // is unlimited from below and we are adding an upper bound for it - propagate_monomial_on_right_side(i->poly().m_coeffs[the_only_unlim], rs, i); - } - - bool conflict() const { return m_explanation.size() > 0; } - - bool get_var_lower_bound(var_index i, mpq & bound) const { - return m_var_infos[i].get_lower_bound(bound); - } - - bool get_var_lower_bound_test(var_index i, mpq & bound, const vector& bs) const { - const auto & t = bs[i]; - if (t.m_expl_lower == -1) - return false; - bound = t.m_lower_bound; - return true; - } - bool get_var_upper_bound_test(var_index i, mpq & bound, const vector& bs) const { - const auto & t = bs[i]; - if (t.m_expl_upper == -1) - return false; - bound = t.m_upper_bound; - return true; - } - - - - bool get_var_upper_bound(var_index i, mpq & bound) const { - const var_info & v = m_var_infos[i]; - return v.get_upper_bound(bound); - } - - // returns the reason for the lower bound, which is the index of the literal, - // or -1 if the bound does not exist - int get_lower_for_monomial_expl(const monomial & p) const { - lp_assert(p.coeff() != 0); - - const auto & v = m_var_infos[p.var()]; - return p.coeff() > 0? v.get_lower_bound_expl():v.get_upper_bound_expl(); - } - - bool lower_for_monomial_exists(const monomial &p ) const { - return get_lower_for_monomial_expl(p) != -1; - } - - bool lower_for_monomial_exists_test(const monomial &p, const vector& bs ) const { - return is_pos(p.coeff())? bs[p.var()].m_expl_lower != -1: bs[p.var()].m_expl_upper != -1; - } - - bool upper_for_monomial_exists(const monomial &p ) const { - return get_upper_for_monomial_expl(p) != -1; - } - bool upper_for_monomial_exists_test(const monomial &p, const vector & bs ) const { - const auto &t = bs[p.var()]; - TRACE("upper_for_monomial_exists_test", tout << p.coeff() << ", " << var_name(p.var()); - tout << "t.m_expl_lower = " << t.m_expl_lower << ", t.m_expl_upper = " << t.m_expl_upper;); - return is_neg(p.coeff())? t.m_expl_lower != -1: t.m_expl_upper != -1; - } - - int get_upper_for_monomial_expl(const monomial & p) const { - lp_assert(p.coeff() != 0); - const auto &v = m_var_infos[p.var()]; - return p.coeff() > 0? v.get_upper_bound_expl() : v.get_lower_bound_expl(); - } - - bool lower_monomial_test(const monomial & p, mpq & lb, vector& bs) const { - lp_assert(p.coeff() != 0); - mpq var_bound; - if (p.coeff() > 0) { - if (!get_var_lower_bound_test(p.var(), var_bound, bs)) - return false; - lb = p.coeff() * var_bound; - } - else { - if (!get_var_upper_bound_test(p.var(), var_bound, bs)) - return false; - lb = p.coeff() * var_bound; - } - return true; - } - // finds the lower bound of the monomial, - // otherwise returns false - bool lower_monomial(const monomial & p, mpq & lb) const { - lp_assert(p.coeff() != 0); - mpq var_bound; - if (p.coeff() > 0) { - if (!get_var_lower_bound(p.var(), var_bound)) - return false; - lb = p.coeff() * var_bound; - } - else { - if (!get_var_upper_bound(p.var(), var_bound)) - return false; - lb = p.coeff() * var_bound; - } - return true; - } - - bool upper_monomial(const monomial & p, mpq & lb) const { - lp_assert(p.coeff() != 0); - mpq var_bound; - if (p.coeff() > 0) { - if (!get_var_upper_bound(p.var(), var_bound)) - return false; - } - else { - if (!get_var_lower_bound(p.var(), var_bound)) - return false; - } - lb = p.coeff() * var_bound; - return true; - } - - bool upper_monomial_test(const monomial & p, mpq & lb, const vector & bs) const { - lp_assert(p.coeff() != 0); - mpq var_bound; - if (p.coeff() > 0) { - if (!get_var_upper_bound_test(p.var(), var_bound, bs)) - return false; - } - else { - if (!get_var_lower_bound_test(p.var(), var_bound, bs)) - return false; - } - lb = p.coeff() * var_bound; - return true; - } - - - - // returns false if not limited from below - // otherwise the answer is put into lb - bool lower(const_constr* c, mpq & lb) const { - return lower(c->poly(), lb); - } - - // returns false if not limited from below - // otherwise the answer is put into lb - mpq lower_val(const_constr * c) const { - return lower_val(*c); - } - - mpq lower_val(const_constr & i) const { - mpq lb; -#if Z3DEBUG - bool r = -#endif - lower(i.poly(), lb); - lp_assert(r); - return lb; - } - - - // returns false if not limited from below - // otherwise the answer is put into lb - bool lower(const polynomial & f, mpq & lb) const { - lb = f.m_a; - mpq lm; - for (const auto & p : f.m_coeffs) { - if (lower_monomial(p, lm)) { - lb += lm; - } else { - return false; - } - } - return true; - } - - - // returns the number of lower unlimited monomials - 0, 1, >=2 - // if there is only one lower unlimited then the index is put into the_only_unlimited - int lower_analize(const_constr * f, unsigned & the_only_unlimited) const { - int ret = 0; - for (unsigned j = 0; j < f->poly().m_coeffs.size(); j++) { - int k; - if ((k = get_lower_for_monomial_expl(f->poly().m_coeffs[j])) == -1) { - if (ret == 1) - return 2; - ret ++; - the_only_unlimited = j; - } - } - return ret; - } - - - - bound_result lower_without(const polynomial & p, var_index j) const { - for (const auto & t: p.m_coeffs) { - if (t.var() == j) - continue; - if (!lower_for_monomial_exists(t)) { - return bound_result(); - } - } - // if we are here then there is a lower bound for p - mpq bound = p.m_a; - for (const auto & t: p.m_coeffs) { - if (t.var() == j) - continue; - - mpq l; - lower_monomial(t, l); - bound += l; - } - return bound_result(bound,bound_type::LOWER); - } - - bound_result lower_without_test(const polynomial & p, var_index j, vector& bs ) const { - for (const auto & t: p.m_coeffs) { - if (t.var() == j) - continue; - if (!lower_for_monomial_exists_test(t, bs)) { - return bound_result(); - } - } - // if we are here then there is a lower bound for p - mpq bound = p.m_a; - for (const auto & t: p.m_coeffs) { - if (t.var() == j) - continue; - - mpq l; - lower_monomial_test(t, l, bs); - bound += l; - } - return bound_result(bound,bound_type::LOWER); - } - - // a is the coefficient before j - bound_result bound_on_polynomial(const polynomial & p, const mpq& a, var_index j) const { - lp_assert(!is_zero(a)); - bound_result r = lower_without(p, j); - if (r.m_type == bound_type::UNDEF) - return r; - if (numeric_traits::is_pos(a)) { - r.m_bound = - ceil_ratio(r.m_bound , a); - r.m_type = bound_type::UPPER; - return r; - } - else { - r.m_bound = -floor_ratio(r.m_bound, a); - r.m_type = bound_type::LOWER; - return r; - } - } - - - - bound_result bound(const_constr * q, var_index j) const { - const mpq& a = q->poly().coeff(j); - return bound_on_polynomial(q->poly(), a, j); - } - - bound_result bound(const polynomial& q, var_index j) const { - const mpq& a = q.coeff(j); - return bound_on_polynomial(q, a, j); - } - - - bound_result bound(unsigned constraint_index, var_index j) const { - return bound(m_asserts[constraint_index], j); - } - - - void print_constraint(std::ostream & out, const_constr & i ) const { - out << (i.is_lemma()? "lemma ": "assert "); - out << "id = " << i.id() << ", "; - if (!i.is_ineq()) { - out << i.divider() << " | "; - } - out << pp_poly(*this, i.poly()); - if (i.is_ineq()) { - out << " <= 0"; - } - out << " active = " << m_active_set.contains(&i) << " "; - mpq b; - if (lower(&i, b)) { - out << ", lower = " << b; - if (is_pos(b)) - out << " INF"; - } else { - out << ", no lower"; - } - - bool all_vars_are_fixed = true; - for (const auto & p : i.poly().coeffs()) { - if (!m_var_infos[p.var()].is_fixed()) { - all_vars_are_fixed = false; - break; - } - } - - if (all_vars_are_fixed) { - auto v = value(i.poly()); - out << ", value = " << v; - if (is_pos(v)) - out << " INF"; - } - out << std::endl; - } - - void print_literal(std::ostream & out, const literal & t) const { - out << (t.is_decided()? "decided ": "implied "); - out << var_name(t.var()) << " "; - if (t.is_lower()) - out << ">= "; - else - out << "<= "; - out << t.bound(); - - if (t.is_decided() == false) { - out << " by constraint " << pp_constraint(*this, *(t.cnstr())); - } else { - out << " decided on trail element " << t.decision_context_index(); - if (m_trail[t.decision_context_index()].tight_constr() != nullptr) { - out << " with tight ineq " << pp_constraint(*this, *m_trail[t.decision_context_index()].tight_constr()); - } - out << "\n"; - } - if (t.tight_constr() != nullptr) { - out << "tight_constr() = " << pp_constraint(*this, *t.tight_constr()) << "\n"; - } - } - - void print_polynomial(std::ostream & out, const polynomial & p) const { - vector> pairs = to_pairs(p.m_coeffs); - this->print_linear_combination_of_column_indices_std(pairs, out); - if (!is_zero(p.m_a)) { - if (p.m_a < 0) { - out << " - " << -p.m_a; - } else { - out << " + " << p.m_a; - } - } - } - - void trace_print_polynomial(std::ostream & out, const polynomial & p) const { - vector> pairs = to_pairs(p.m_coeffs); - this->print_linear_combination_of_column_indices_std(pairs, out); - if (!is_zero(p.m_a)) { - if (p.m_a < 0) { - out << " - " << -p.m_a; - } else { - out << " + " << p.m_a; - } - } - out << "\n"; - for (auto m : p.coeffs()) { - unsigned j = m.var(); - out << "domain of " << var_name(j) << " = "; - print_var_domain(out, j); - } - mpq b; - if (lower(p, b)) { - out << ", lower = " << b; - if (is_pos(b)) - out << " INF"; - } else { - out << ", no lower"; - } - - bool all_vars_are_fixed = true; - for (const auto & pp : p.coeffs()) { - if (!m_var_infos[pp.var()].is_fixed()) { - all_vars_are_fixed = false; - break; - } - } - - if (all_vars_are_fixed) { - auto v = value(p); - out << ", value = " << v; - if (is_pos(v)) - out << " INF"; - } - } - - // trying to improve constraint "ie" by eliminating var j by using the tight inequality - // for j. the left side of the inequality is passed as a parameter. - bool resolve(polynomial & ie, unsigned j, bool sign_j_in_ti_is_pos, const polynomial & ti) const { - TRACE("resolve_int", tout << "ie = " << pp_poly(*this, ie); - tout << ", j = " << j << "(" << var_name(j) << ")" << ", sign_j_in_ti_is_pos = " << sign_j_in_ti_is_pos << ", ti = " << pp_poly(*this, ti) << "\n";); - lp_assert(ti.var_coeff_is_unit(j)); - lp_assert(is_pos(ti.coeff(j)) == sign_j_in_ti_is_pos); - auto &coeffs = ie.m_coeffs; - // todo: implement a more efficient version - bool found = false; - mpq a; - for (const auto & c : coeffs) { - if (c.var() == j) { - a = c.coeff(); - found = true; - break; - } - } - - if (!found) { - TRACE("resolve_int", tout << " no change\n";); - return false; - } - - if (is_neg(a)) { - if (!sign_j_in_ti_is_pos) - return false; - a = -a; - } else { - if (sign_j_in_ti_is_pos) - return false; - } - - for (auto & c : ti.m_coeffs) { - ie += monomial(a * c.coeff(), c.var()); - } - - ie.m_a += a * ti.m_a; - TRACE("resolve_int", tout << "ie = " << pp_poly(*this, ie) << "\n";); - return true; - } - - // returns true iff p imposes a better bound on j - bool improves(var_index j, const_constr * p) const { - auto a = p->coeff(j); - if (is_zero(a)) - return false; - const auto& dom = m_var_infos[j].domain(); - if (dom.is_empty()) - return false; - if (is_pos(a)) { - bound_result new_upper = bound(p, j); - if (new_upper.m_type == bound_type::UNDEF) - return false; - return dom.improves_with_upper_bound(new_upper.bound()); - } - - lp_assert(is_neg(a)); - bound_result new_lower = bound(p, j); - if (new_lower.m_type == bound_type::UNDEF) - return false; - return dom.improves_with_lower_bound(new_lower.bound()); - } - - - // returns true iff br imposes a better bound on j - bool improves(var_index j, const bound_result & br) const { - if (br.m_type == bound_type::UNDEF) - return false; - const auto& dom = m_var_infos[j].domain(); - if (dom.is_empty()) - return false; - if (br.m_type == bound_type::UPPER) { - mpq b; - bool j_has_upper_bound = get_var_upper_bound(j, b); - return (!j_has_upper_bound || br.bound() < b) && - !dom.intersection_with_upper_bound_is_empty(br.bound()); - } - - if (br.m_type == bound_type::UNDEF) - return false; - mpq b; - bool lower_bound_exists = get_var_lower_bound(j, b); - return (!lower_bound_exists || br.bound() > b) && - !dom.intersection_with_lower_bound_is_empty(br.bound()); - } - - - void add_bound(mpq v, unsigned j, bool is_lower, constraint * c) { - literal l = literal::make_implied_literal(j, is_lower, v, c, m_var_infos[j].external_stack_level()); - push_literal_to_trail(l); - m_var_infos[j].number_of_bound_propagations()++; - } - - mpq value(const polynomial & p) const { - mpq ret= p.m_a; - for (const auto & t:p.m_coeffs) - ret += t.coeff() * m_v[t.var()]; - return ret; - } - - bool flip_coin() { - return m_settings.random_next() % 2 == 0; - } - - void pop_to_external_level() { - while (m_decision_level > 0) { - pop(); - } - } - - void print_state_stats(std::ostream & out) const { - if (m_bounded_search_calls % 10 ) - return; - out << "search_calls = " << m_bounded_search_calls << ", "; - out << "vars = " << m_var_infos.size() << ","; - out << "asserts = "<< m_asserts.size() << ", "; - out << "lemmas = " << m_lemmas_container.size() << ", "; - out << "trail = " << m_trail.size() << ", props = " << m_number_of_constraints_tried_for_propagaton << ", "; - out << "cnfls = " << m_number_of_conflicts << ", "; - } - - - lbool check() { - init_search(); - TRACE("check_int", tout << "starting check" << "\n";); - while (!m_stuck_state && !cancel()) { - TRACE("cs_ch", tout << "inside loop\n";); - lbool r = bounded_search(); - if (cancel()) { - break; - } - if (r != lbool::l_undef) { - TRACE("check_int", tout << "return " << (r == lbool::l_true ? "true" : "false") << "\n"; ); - pop_to_external_level(); - return r; - } - restart(); - simplify_problem(); - if (check_inconsistent()) { - TRACE("check_int", tout << "return " << (r == lbool::l_true ? "true" : "false") << "\n"; ); - pop_to_external_level(); - return lbool::l_false; - } - gc(); - TRACE("check_int", tout << "continue\n"; print_state_stats(tout); ); - - } - TRACE("check_int", tout << "return undef\n";); - pop_to_external_level(); - return lbool::l_undef; - } - - unsigned find_unused_index() const { - for (unsigned j = m_var_infos.size(); ; j++) - if (m_user_vars_to_chase_cut_solver_vars.find(j) == m_user_vars_to_chase_cut_solver_vars.end()) - return j; - - } - - - void init_search() { - lp_assert(m_explanation.size() == 0); - m_number_of_conflicts = 0; - m_bounded_search_calls = 0; - m_stuck_state = false; - m_cancelled = false; - for (auto & vi : m_var_infos) - vi.number_of_bound_propagations() = 0; - m_number_of_constraints_tried_for_propagaton = 0; - m_pushes_to_trail = 0; - } - - constraint* propagate_constraint(constraint* c) { - m_number_of_constraints_tried_for_propagaton ++; - lp_assert(c->is_ineq()); - TRACE("ba_int", trace_print_constraint(tout, c);); - // consider a special case for a constraint with just two variables - unsigned the_only_unlim; - int r = lower_analize(c, the_only_unlim); - if (r == 0) { - mpq b; - lower(c->poly(), b); - if (is_pos(b)) { - TRACE("cs_inconsistent", tout << "incostistent constraint "; - trace_print_constraint(tout, c); - tout << "\nlevel = " << m_decision_level << std::endl;); - return c; - } - propagate_constraint_on_lower(c, b); - } else if (r == 1 && - !too_many_propagations_for_var(c->poly().coeffs()[the_only_unlim].var())) { - // otherwise the new bound, even if it exists, will be rejected! - lp_assert(!lower_is_pos(c->poly())); - propagate_constraint_only_one_unlim(c, the_only_unlim); - } - lp_assert(!lower_is_pos(c->poly())); - return nullptr; - } - - void print_trail(std::ostream & out) const { print_trail(out, m_trail.size()); } - - void print_trail(std::ostream & out, unsigned last_index) const { - out << "trail\n"; - for (unsigned j = 0; j <= last_index && j < m_trail.size(); j++ ) { - out << "offset = " << j << ", "; - print_literal(out, m_trail[j]); - } - out << "end of trail\n"; - } - void print_state(std::ostream & out) const { - out << "asserts total " << m_asserts.size() << "\n"; - for (const auto i: m_asserts) { - print_constraint(out, *i); - } - out << "end of constraints\n"; - out << "lemmas total " << m_lemmas_container.size() << "\n"; - for (const auto i: m_lemmas_container.m_lemmas) { - print_constraint(out, *i); - } - out << "end of constraints\n"; - print_trail(out); - out << "var_infos\n"; - for (const auto & v: m_var_infos) { - if (v.is_active()) - print_var_info(out, v); - } - out << "end of var_infos\n"; - out << "level = " << m_decision_level << "\n"; - out << "end of state dump, bounded_search_calls = " << m_bounded_search_calls << ", number_of_conflicts = " << m_number_of_conflicts << std::endl; - } - lbool bounded_search() { - m_bounded_search_calls++; - lbool is_sat = propagate_and_backjump_step(); - if (is_sat != lbool::l_undef) { - TRACE("decide_int", tout << "returning " << (int)is_sat << "\n";); - return is_sat; - } - gc(); - if (cancel()) - return lbool::l_undef; - if (!decide()) { - TRACE("decide_int", tout << "going to final_check()\n";); - lbool is_sat = final_check(); - if (is_sat != lbool::l_undef) { - return is_sat; - } - } - TRACE("decide_int", tout << "returning undef\n";); - return lbool::l_undef; - } - - bool pick_bound_kind(unsigned j) { - var_info & vi = m_var_infos[j]; - if (lower_bound_exists(vi) && upper_bound_exists(vi)) - return flip_coin(); - if (lower_bound_exists(vi)) - return true; - lp_assert(upper_bound_exists(vi)); - return false; - } - - bool decide() { - int j = find_var_for_deciding(); - if (j < 0) - return false; - TRACE("decide_int", tout << "going to decide " << var_name(j) << " var index = " << j << "\n"; - tout << "domain = "; print_var_domain(tout, j); tout << ", m_decision_level="< bound || m_number_of_conflicts > bound) { - m_cancelled = true; - return true; - } - return false; - } - - bool cancel_when_propagation_speed_is_too_slow() { - if (m_number_of_constraints_tried_for_propagaton >= 10000) { - m_cancelled = m_pushes_to_trail <= 50; - m_pushes_to_trail = 0; - m_number_of_constraints_tried_for_propagaton = 0; - return m_cancelled; - } - return false; - } - - - lbool propagate_and_backjump_step() { - do { - if (cancel_when_propagation_speed_is_too_slow()) - return lbool::l_undef; - constraint* c = propagate(); - if (cancel()) - return lbool::l_undef; - TRACE("cs_dec", tout << "trail = \n"; print_trail(tout); tout << "end of trail\n";); - - if (c != nullptr) { - m_number_of_conflicts++; - TRACE("propagate_and_backjump_step_int", tout << "incostistent_constraint "; trace_print_constraint(tout, c); tout << "m_number_of_conflicts = " << m_number_of_conflicts << std::endl; tout << "chase_cut_solver_calls " << m_settings.st().m_chase_cut_solver_calls << std::endl;); - if (at_base_lvl()) { - fill_conflict_explanation(c); - return lbool::l_false; - } - handle_conflicting_cores(); // testing only - resolve_conflict(c); - } - } - while (!m_active_set.is_empty()); - - return !all_vars_are_fixed()? lbool::l_undef :lbool::l_true; - } - - bool decision_is_redundant_for_constraint(const polynomial& i, literal & l) const { - const mpq & coeff = i.coeff(l.var()); - if (is_zero(coeff)) - return true; - return is_pos(coeff)? ! l.is_lower(): l.is_lower(); - } - - bool is_divizible(const mpq & a, const mpq & b) const { - lp_assert(!is_zero(b)); - return is_zero(a % b); - } - - void create_div_ndiv_parts_for_tightening(const polynomial & p, const mpq & coeff, polynomial & div_part, polynomial & ndiv_part) { - for (const auto &m : p.m_coeffs) { - if (is_divizible(m.coeff(), coeff)){ - div_part.m_coeffs.push_back(m); - } else { - ndiv_part.m_coeffs.push_back(m); - } - } - - TRACE("tight", - tout << "div_part = " << pp_poly(*this, div_part) << "\n"; - tout << "ndiv_part = " << pp_poly(*this, ndiv_part) << "\n";); - } - - - void add_tight_constr_of_literal(polynomial &ndiv_part, - const mpq & c, - literal &l) { - lp_assert(is_pos(c)); - ndiv_part.add(c, m_trail[l.decision_context_index()].tight_constr()->poly()); - lp_assert(is_zero(ndiv_part.coeff(l.var()))); - TRACE("tight", tout << "ndiv_part = " << pp_poly(*this, ndiv_part) << "\n";); - } - - void decided_lower(const mpq & a, - const mpq & c, - polynomial &div_part, - polynomial &ndiv_part, - literal &l) { - mpq k = is_pos(a)?ceil( c / a): floor(c / a); - ndiv_part += monomial(-c, l.var()); // it will be moved to div_part - mpq a_k = a * k; - mpq m = a_k - c; - TRACE("tight", tout << "c = " << c << ", a = " << a << - ", c / a = " << c/a << ", k = " << - k << ", a * k = " << a * k << ", m = " << m << "\n"; ); - lp_assert(!is_neg(m)); - create_tight_constr_under_literal(l.decision_context_index()); - const literal & lex = m_trail[l.decision_context_index()]; - lp_assert(lex.var() == l.var()); - for (const monomial & n : lex.tight_constr()->poly().coeffs()) { - if (n.var() == l.var()) { - lp_assert(n.coeff() == one_of_type()); - div_part += monomial(a_k, l.var()); - } else { - ndiv_part += monomial(m * n.coeff(), n.var()); - } - } - ndiv_part += m * lex.tight_constr()->poly().m_a; - TRACE("tight", tout << "Decided-Lower "; - tout << "div_part = " << pp_poly(*this, div_part) << "\n"; - tout << "ndiv_part = " << pp_poly(*this, ndiv_part) << "\n";); - } - - void decided_upper(const mpq & a, - const mpq & c, - polynomial &div_part, - polynomial &r, - literal &l) { - // we would like to have c - ak > 0, or ak < c - mpq k = is_pos(a)? floor( c / a): ceil(c / a); - r += monomial(-c, l.var()); // it will be moved to div_part - - mpq a_k = a * k; - mpq m = c - a_k; - TRACE("tight", tout << "c = " << c << ", a = " << a << - ", c / a = " << c/a << ", k = " << - k << ", a * k = " << a * k << ", m = " << m << "\n"; ); - lp_assert(!is_neg(m)); - create_tight_constr_under_literal(l.decision_context_index()); - const literal & lex = m_trail[l.decision_context_index()]; - lp_assert(lex.var() == l.var()); - for (const monomial & n : lex.tight_constr()->poly().coeffs()) { - if (n.var() == l.var()) { - lp_assert(n.coeff() == -one_of_type()); - div_part += monomial(a_k, l.var()); - } else { - r += monomial(m * n.coeff(), n.var()); - } - } - r += m * lex.tight_constr()->poly().m_a; - TRACE("tight", tout << "Decided-Lower "; - tout << "div_part = " << pp_poly(*this, div_part) << "\n"; - tout << "r = " << pp_poly(*this, r) << "\n";); - } - - // returns true iff there was a change - bool tighten_on_literal(const mpq & a, - polynomial & div_part, - polynomial & ndiv_part, - int index_of_literal) { - bool change = false; - literal & l = m_trail[index_of_literal]; - if (l.tight_constr() == nullptr) { - create_tight_constr_under_literal(index_of_literal); - } - TRACE("tight", - tout << "index_of_literal = " << index_of_literal << ", "; - print_literal(tout, m_trail[index_of_literal]);); - if (l.is_implied()) { // Resolve-Implied - change = resolve(ndiv_part, l.var(), !l.is_lower(), l.tight_constr()->poly()); - TRACE("tight", tout << "after resolve ndiv_part = " << pp_poly(*this, ndiv_part) << "\n";); - } else { - create_tight_constr_under_literal(l.decision_context_index()); - TRACE("tight", - tout << "index_of_literal = " << index_of_literal << ", "; - print_literal(tout, m_trail[index_of_literal]); tout << "\n"; - tout << "div_part = " << pp_poly(*this, div_part) << "\n"; - tout << "ndiv_part = " << pp_poly(*this, ndiv_part) << "\n"; - tout << "a = " << a << "\n"; - ); - mpq c = ndiv_part.coeff(l.var()); - if (is_zero(c)) - return false; - if (l.is_lower()) { - if (is_neg(c)) - add_tight_constr_of_literal(ndiv_part, -c, l); - else - decided_lower(a, c, div_part, ndiv_part, l); - } else { - if (is_pos(c)) // Decided-Upper-Pos - add_tight_constr_of_literal(ndiv_part, c, l); - else - decided_upper(a, c, div_part, ndiv_part, l); - } - change = true; - } - return change; - } - - void eliminate_ndiv_part_monomials(const mpq& a, polynomial & div_part, polynomial& ndiv_part, unsigned index_of_literal) { - if (ndiv_part.number_of_monomials() == 0) - return; - int k = index_of_literal - 1; - lp_assert(k >= 0); - literal& l = m_trail[index_of_literal]; - while (ndiv_part.number_of_monomials() > 0) { - if (tighten_on_literal(a, div_part, ndiv_part, k)) { - literal & lk = m_trail[k]; - l.tight_constr()->add_predecessor( - lk.is_implied()? - lk.tight_constr() : - m_trail[lk.decision_context_index()].tight_constr() - ); - } - k--; - } - } - - // see page 88 of Cutting to Chase - void tighten(constraint * c, - unsigned j_of_var, - const mpq& a, - unsigned index_of_literal) { - polynomial div_part, ndiv_part; - ndiv_part.m_a = c->poly().m_a; - create_div_ndiv_parts_for_tightening(c->poly(), a, div_part, ndiv_part); - eliminate_ndiv_part_monomials(a, div_part, ndiv_part, index_of_literal); - mpq abs_a = abs(a); - polynomial & p = c->poly(); - p.clear(); - for (const auto & m : div_part.m_coeffs) { - p.m_coeffs.push_back(monomial(m.coeff() / abs_a, m.var())); - } - p.m_a = ceil(ndiv_part.m_a / abs_a); - lp_assert(p.m_a >= ndiv_part.m_a / abs_a); - TRACE("tight", tout << "index_of_literal = " << index_of_literal << ", got tight p = " << pp_poly(*this, p) << "\n";); - } - - void create_tight_constr_under_literal(unsigned index_of_literal) { - literal & l = m_trail[index_of_literal]; - if (l.tight_constr() != nullptr) - return; - if (l.is_decided()) { - create_tight_constr_under_literal(l.decision_context_index()); - return; - } - TRACE("tight", tout << "index_of_literal = " << index_of_literal << "\n";); - const_constr* c = l.cnstr(); - lp_assert(c->is_ineq()); - unsigned j = l.var(); - const mpq& a = c->poly().coeff(j); - lp_assert(!is_zero(a)); - if (a == one_of_type() || a == - one_of_type()) { - l.tight_constr() = l.cnstr(); - return; - } - l.tight_constr() = new constraint( get_new_constraint_id(), - c->assert_origins(), - c->poly(), - true); - l.tight_constr()->add_predecessor(c); - tighten(l.tight_constr(), j, a, index_of_literal); - add_lemma_as_not_active(l.tight_constr()); - lp_assert(constraint_indices_are_correct()); - } - - void backjump(unsigned index_of_literal, - bool lemma_has_been_modified, - constraint* lemma, - constraint * orig_conflict) { - polynomial &p = lemma->poly(); - lp_assert(m_trail[index_of_literal].is_decided()); - unsigned j = m_trail[index_of_literal].var(); - while(m_trail.size() > index_of_literal) { pop(); } - lp_assert(m_trail.size() == index_of_literal); - bound_result br = bound_on_polynomial(p, p.coeff(j), j); - - TRACE("int_backjump", br.print(tout); - print_var_info(tout, j); - tout << "p = " << pp_poly(*this, p) << "\n";); - if (!improves(j, br)) { - TRACE("int_backjump", br.print(tout); tout << "\nimproves is false\n"; - tout << "var info after pop = "; print_var_info(tout, j);); - if (lemma_has_been_modified) { // flip_coin() gives - // priority 0 or 1, so it is be dequeued fast - add_lemma(lemma, flip_coin()); - } else { - m_active_set.add_constraint(orig_conflict, flip_coin()); - } - } else { - constraint *c; - if (lemma_has_been_modified) { - c = lemma; - add_lemma_as_not_active(lemma); - } else { - c = orig_conflict; - } - add_bound(br.bound(), j, br.m_type == bound_type::LOWER, c); - restrict_var_domain_with_bound_result(j, br, m_trail.size() - 1); - lp_assert(!m_var_infos[j].domain().is_empty()); - TRACE("int_backjump", tout << "done resolving:\nvar info after restricton = "; - print_var_info(tout, j); - tout << "new literal = "; print_literal(tout, m_trail.back());); - lp_assert(!lower_is_pos(c->poly())); - } - } - - void print_resolvent(std::ostream& out, const polynomial& p, literal &l) const { - out << "p = "; trace_print_polynomial(out, p); - mpq rr; - bool bb = lower(p, rr); - if (!bb) { - out << "\nlower(p) is not defined\n"; - } else { - out << "\nlower(p) = " << rr << "\n"; - } - - out << "tight_constr = " << pp_constraint(*this, *l.tight_constr()) << "\n"; - out << "constraint = " << pp_constraint(*this, *l.cnstr()) << "\n"; - out << "var domains" << "\n"; - for (auto & m : l.tight_constr()->poly().coeffs()) { - out << "var = " << m.var() << " " << var_name(m.var()) << " "; - print_var_domain(out, m.var()); - out << " "; - } - out << "\n"; - } - - void trace_resolve_print(std::ostream& out, const polynomial & p, literal & l, unsigned index_of_literal) { - out << "index_of_literal = " << index_of_literal <<", p = " << pp_poly(*this, p) << "\n"; - out << "l = "; print_literal(out, l); - out << "lower(p) = " << lower_no_check(p) << "\n"; - for (auto & m : p.coeffs()) { - out << var_name(m.var()) << " "; - print_var_domain(out, m.var()); - out << " "; - } - out << "\nm_decision_level = " << m_decision_level << "\n"; - } - - - bool resolve_decided_literal(unsigned index_of_literal, literal& l, bool lemma_has_been_modified, constraint* lemma, constraint *orig_conflict) { - if (decision_is_redundant_for_constraint(lemma->poly(), l)) { - do { pop();} while(m_trail.size() > index_of_literal); - lp_assert(m_trail.size() == index_of_literal); - TRACE("resolve_decided_literal", tout << "n "; print_literal(tout, l); if (m_decision_level == 0) tout << ", done resolving";); - lp_assert(lower_is_pos(lemma->poly())); - return m_decision_level == 0; - } - handle_conflicting_cores(); - backjump(index_of_literal, lemma_has_been_modified, lemma, orig_conflict); - return true; // done - } - - // applying Resolve rule - void resolve_implied_literal(literal & l, - bool &lemma_has_been_modified, - constraint* lemma) { - polynomial & p = lemma->poly(); - lp_assert(lower_is_pos(p)); - if (!resolve(p, l.var(), !l.is_lower(), l.tight_constr()->poly())) - return; - lemma_has_been_modified = true; - lemma->add_predecessor(l.tight_constr()); - TRACE("resolve_implied_literal", tout << "resolved p: "; print_resolvent(tout, p, l);); - lp_assert(lower_is_pos(p)); - } - - // returns true iff resolved - bool resolve_conflict_for_inequality_on_trail_element(unsigned index_of_literal, bool & lemma_has_been_modified, constraint* lemma, constraint * orig_conflict) { - lp_assert(lower_is_pos(lemma->poly())); - literal & l = m_trail[index_of_literal]; - TRACE("resolve_conflict_for_inequality_on_trail_element", - tout << "index_of_literal = " << index_of_literal <<", p = " << pp_poly(*this, lemma->poly()) << "\n"; - tout << "\nm_decision_level = " << m_decision_level << "\n"; - print_literal(tout, l);); - if (l.is_decided()) { - return resolve_decided_literal(index_of_literal, l, lemma_has_been_modified, lemma, orig_conflict); - } else { - create_tight_constr_under_literal(index_of_literal); - resolve_implied_literal(l, lemma_has_been_modified, lemma); - return false; - } - } - - bool lower_is_pos(const_constr* c) const { return lower_is_pos(c->poly()); } - - bool lower_is_pos(const polynomial & p) const { - mpq b; - bool r = lower(p, b); - return r && is_pos(b); - } - - mpq lower_no_check(const polynomial & p) const { - mpq b; -#if Z3DEBUG - bool r = -#endif - lower(p, b); -#if Z3DEBUG - lp_assert(r); -#endif - return b; - } - - unsigned get_new_constraint_id() { return m_indexer_of_constraints.get_new_index(); } - - void resolve_conflict_for_inequality(constraint * c) { - // Create a new constraint, almost copy of "c", that becomes a lemma and could be modified later - constraint *lemma = new constraint(get_new_constraint_id(), c->poly(), true); - lemma->add_predecessor(c); - - TRACE("resolve_conflict_for_inequality", print_constraint(tout, *lemma);); - lp_assert(lower_is_pos(lemma->poly())); - bool done = false; - unsigned j = m_trail.size() - 1; - bool lemma_has_been_modified = false; - unsigned number_of_lemmas = m_lemmas_container.size(); - while (!done) { - if (cancel()) { - return; - } - done = resolve_conflict_for_inequality_on_trail_element(j--, lemma_has_been_modified, lemma, c); - if (j >= m_trail.size()) { - TRACE("resolve_conflict_for_inequality", tout << "adjust j";); - lp_assert(m_trail.size()); - j = m_trail.size() - 1; - } - } - - if (number_of_lemmas == m_lemmas_container.size()) { - if (lemma_has_been_modified && lemma->poly().number_of_monomials() != 0) { - add_lemma_as_not_active(lemma); - } - else { - delete_constraint(lemma); - } - } - } - - void delete_constraint(constraint * c) { - m_indexer_of_constraints.release_index(c->id()); - delete c; - } - - void resolve_conflict(constraint * i) { - lp_assert(!at_base_lvl()); - TRACE("int_resolve_confl", tout << "inconstistent_constraint = "; - print_constraint(tout, *i); print_state(tout);); - if (i->is_ineq()) { - resolve_conflict_for_inequality(i); - } else { - lp_assert(false); // not implemented - } - } - - - void print_scope(std::ostream& out) const { - for (const scope & s : m_scopes) { - out << ", trail_size = " << s.m_trail_size << "\n"; - } - } - -public: - unsigned number_of_asserts() const { return m_asserts.size(); } - - void push() { - m_scopes.push_back(scope(m_trail.size())); - TRACE("pp_cs", tout << "level = " << m_scopes.size() << ", trail size = " << m_trail.size();); - } - - void pop_last_assert() { - constraint * i = m_asserts.back();; - for (auto & p: i->poly().m_coeffs) { - m_var_infos[p.var()].remove_depended_constraint(i); - } - lp_assert(i->assert_origins().size() == 1); - for (constraint_index ci : i->assert_origins()) - m_lemmas_container.submit_assert_origin_for_delete(ci); - m_active_set.remove_constraint(i); - delete_constraint(i); - m_asserts.pop_back(); - } - - void pop_constraints() { - vector gone_lemmas = m_lemmas_container.remove_lemmas_depending_on_submitted_origins(); - - for (constraint * i : gone_lemmas) { - for (auto & p: i->poly().m_coeffs) { - m_var_infos[p.var()].remove_depended_constraint(i); - } - m_active_set.remove_constraint(i); - delete_constraint(i); - } - } - -public: - void pop_trail(unsigned k) { - unsigned new_scope_size = m_scopes.size() - k; - scope s = m_scopes[new_scope_size]; - m_scopes.shrink(new_scope_size); - for (unsigned j = m_trail.size(); j-- > s.m_trail_size; ) { - literal& lit = m_trail[j]; - if (lit.is_decided()) - m_decision_level--; - else - m_active_set.add_constraint(lit.cnstr(), 2); // 2 is a priority; this constraint will not be processed urgentlty - var_info & vi = m_var_infos[lit.var()]; - if (vi.external_stack_level() != lit.prev_var_level()) - vi.pop(1, lit.prev_var_level()); - m_trail.pop_back(); - } - } - void pop(unsigned k) { - pop_trail(k); - pop_constraints(); - lp_assert(var_infos_are_correct()); - } -public: - void pop() { pop(1); } - - chase_cut_solver(std::function var_name_function, - std::function print_constraint_function, - std::function number_of_variables_function, - std::function var_value_function, - lp_settings & settings - ) : m_var_name_function(var_name_function), - m_print_constraint_function(print_constraint_function), - m_number_of_variables_function(number_of_variables_function), - m_var_value_function(var_value_function), - m_settings(settings), - m_decision_level(0) - {} - - - int find_conflicting_core(const constraint* &lower, const constraint* & upper) const { - for (unsigned j = 0; j < m_var_infos.size(); j++) { - if (is_cc(j, lower, upper)) - return j; - } - return -1; - } - - void list_confl_cores() { - const constraint* lower; const constraint* upper; - for (unsigned j = 0; j < m_var_infos.size(); j++) { - if (is_cc(j, lower, upper)) { - std::cout << "confl core = "; print_var_info(std::cout, j); - std::cout << "lower = "; print_constraint(std::cout, *lower); - std::cout << "upper = "; print_constraint(std::cout, *upper); - } - } - } - - void handle_conflicting_cores() { - return; - const constraint* lower; - const constraint* upper; - int j = find_conflicting_core(lower, upper); - - if (j >=0) { - std::cout << "confl core = "; print_var_info(std::cout, j); - std::cout << "lower = "; print_constraint(std::cout, *lower); - std::cout << "upper = "; print_constraint(std::cout, *upper); - lp_assert(false); // not implemented - } - } - - constraint* find_constraint_to_propagate() { - handle_conflicting_cores(); - return m_active_set.remove_constraint(); - } - - // returns nullptr if there is no conflict, or a conflict constraint otherwise - constraint* propagate_constraints_on_active_set() { - constraint *c; - while ((c = find_constraint_to_propagate()) != nullptr) { - c = propagate_constraint(c); - if (cancel()) { - return nullptr; - } - if (c != nullptr) { - return c; - } - } - return nullptr; - } - - - // returns -1 if there is no conflict and the index of the conflict constraint otherwise - constraint * propagate() { - constraint* cnf = propagate_constraints_on_active_set();; - if (cnf != nullptr){ - lp_assert(lower_is_pos(cnf)); - return cnf; - } - handle_conflicting_cores(); - return nullptr; - } - - - // walk the trail backward and find the last implied bound on j (of the right kind) - int find_literal_index_after(unsigned j, bool is_lower, unsigned trail_lim) const { - for (unsigned k = trail_lim; k-- > 0;) { - const auto & l = m_trail[k]; - lp_assert(!l.is_decided()); - if (l.var() == j && l.is_lower() == is_lower) - return k; - } - TRACE("find_literal", tout << "cannot find a literal for " << var_name(j)<< " j = " << j << " is_lower = " << is_lower << ", trail_lim = " << trail_lim;); - return -1; - } - - void decide_var_on_bound(unsigned j, bool decide_on_lower) { - mpq b; - vector lhs; - unsigned decision_context_index; - var_info & vi = m_var_infos[j]; - unsigned var_level = vi.external_stack_level(); - push(); - if (decide_on_lower) { - vi.domain().get_lower_bound_with_expl(b, decision_context_index); - vi.intersect_with_upper_bound(b, m_trail.size(), m_scopes.size()); - } - else { - vi.domain().get_upper_bound_with_expl(b, decision_context_index); - vi.intersect_with_lower_bound(b, m_trail.size(), m_scopes.size()); - } - if (j >= m_v.size()) - m_v.resize(j + 1); - m_v[j] = b; - TRACE("decide_var_on_bound", tout<< "j="<< j<<" ";print_var_info(tout, j);); - add_changed_var(j, !decide_on_lower, 1); - m_decision_level++; - literal nl = literal::make_decided_literal(j, !decide_on_lower, b, decision_context_index, var_level); - push_literal_to_trail(nl); - } - - bool consistent(const_constr * i) const { - // an option could be to check that upper(i.poly()) <= 0 - bool ret = value(i->poly()) <= zero_of_type(); - CTRACE("chase_cut_solver_state_inconsistent", !ret, - tout << "inconsistent constraint " << pp_constraint(*this, *i) << "\n"; - tout << "value = " << value(i->poly()) << '\n';); - return ret; - } - - int find_var_for_deciding() const { - unsigned j = m_settings.random_next() % m_var_infos.size(); - - for (unsigned k = 0; k < m_var_infos.size(); k++, j++) { - if (j == m_var_infos.size()) - j = 0; - const auto & d = m_var_infos[j].domain(); - lp_assert(!d.is_empty()); - if (!d.is_fixed() && (lower_bound_exists(d) || upper_bound_exists(d))) - return j; - } - - // start using the rational solution for bounds and branches - - return -1; - } - - bool there_is_var_with_empty_domain() const { - for (unsigned j = 0; j < m_var_infos.size(); j++) { - const auto & d = m_var_infos[j].domain(); - if (d.is_empty()) - return true; - } - return false; - } - - bool all_vars_are_fixed() const { - for (unsigned j = 0; j < m_var_infos.size(); j++) { - if (m_var_infos[j].is_active() && ! m_var_infos[j].is_fixed()) - return false; - } - return true; - } - - bool consistent() const { - if (!all_vars_are_fixed()) { - // this check could be removed if we use upper bound to check if an constraint holds - return false; // ignore the variables values and only return true if every variable is fixed - } - - for (const_constr* i : m_asserts) { - if (!consistent(i)) - return false; - } - return true; - } - - - void simplify_ineq(polynomial & p) const { - TRACE("simplify_ineq_int", tout << "p = " << pp_poly(*this, p) << "\n";); - auto & ms = p.m_coeffs; - if (ms.size() == 0) - return; - mpq g; - if (ms.size() == 1) { - g = abs(ms[0].coeff()); - } else { - g = gcd(ms[0].coeff(), ms[1].coeff()); - for (unsigned j = 2; j < ms.size(); j++) { - g = gcd(g, ms[j].coeff()); - } - lp_assert(is_pos(g)); - } - if (g != one_of_type()) { - for (auto & m : ms) - m.coeff() /= g; - p.m_a = ceil(p.m_a /g); - } - TRACE("simplify_ineq_int", tout << "p = " << pp_poly(*this, p) << "\n";); - } - - void add_lemma_common(constraint* lemma) { - lp_assert(lemma->poly().number_of_monomials() > 0); - m_lemmas_container.add_lemma(lemma); - polynomial & p = lemma->poly(); - simplify_ineq(p); - for (const auto & m : p.coeffs()) { - m_var_infos[m.var()].add_dependent_constraint(lemma, is_pos(m.coeff())); - } - } - - void add_lemma_as_not_active(constraint * lemma) { - add_lemma_common(lemma); - TRACE("add_lemma_int", trace_print_constraint(tout, lemma);); - } - - void add_lemma(constraint * lemma, int priority) { - add_lemma_common(lemma); - m_active_set.add_constraint(lemma, priority); - lp_assert(constraint_indices_are_correct()); - TRACE("add_lemma_int", trace_print_constraint(tout, lemma);); - } - - - unsigned add_ineq(const vector & lhs, - const mpq& free_coeff, - constraint_index origin) { - lp_assert(lhs_is_int(lhs)); - lp_assert(is_int(free_coeff)); - for (auto & p : lhs) { - if (p.var() >= m_var_infos.size()) { - m_var_infos.resize(m_number_of_variables_function()); - } - - var_info & vi = m_var_infos[p.var()]; - - if (!vi.is_active()) { - vi.activate(p.var()); - } - } - - constraint * c = new constraint(get_new_constraint_id(), origin, polynomial(lhs, free_coeff), true); - - lp_assert(constraint_indices_are_correct()); - - m_asserts.push_back(c); - add_constraint_to_dependend_for_its_vars(c); - m_active_set.add_constraint(c, 2); // 2 is a priority - - TRACE("add_ineq_int",tout << "explanation :"; m_print_constraint_function(origin, tout); tout << "\n"; - tout << "m_asserts[" << m_asserts.size() - 1 << "] = "; - print_constraint(tout, *m_asserts.back()); tout << "\n";); - - return m_asserts.size() - 1; - } - - - void add_constraint_to_dependend_for_its_vars(constraint * c) { - for (auto & p : c->poly().coeffs()) { - m_var_infos[p.var()].add_dependent_constraint(c, is_pos(p.coeff())); - } - } - - bool var_has_no_bounds(const var_info& vi) const { - return !lower_bound_exists(vi) && !upper_bound_exists(vi); - } - - unsigned number_of_constraints() const { return m_asserts.size() + m_lemmas_container.size(); } - - void copy_poly_coeffs_to_term(polynomial& poly, lar_term & t) { - for (auto & p :poly.m_coeffs) - t.add_monomial(p.coeff(), p.var()); - } - - bool try_getting_cut(lar_term& t, mpq &k, vector& x) { - // todo : create an efficient version based on var_info bounds - - vector short_lemmas; - unsigned size = static_cast(-1); - for (constraint *c : m_lemmas_container.m_lemmas ) { - if (!c->is_ineq()) continue; - const auto & p = c->poly(); - if (p.number_of_monomials() > size) - continue; - if (is_pos(c->poly().value(x))) { - if (p.number_of_monomials() < size) { - size = p.number_of_monomials(); - // even is size == 1 it makes sense to look for a random cut - short_lemmas.clear(); - } - short_lemmas.push_back(c); - } - } - unsigned n = short_lemmas.size(); - if (n == 0) return false; - constraint *c = short_lemmas[m_settings.random_next() % n]; - k = - c->poly().m_a; - copy_poly_coeffs_to_term(c->poly(), t); - TRACE("chase_cut_solver_cuts", trace_print_constraint(tout, *c);); - return true; - } -}; - -inline polynomial operator*(const mpq & a, polynomial & p) { - polynomial ret; - ret.m_a = p.m_a * a; - - for (const auto & t: p.m_coeffs) - ret.m_coeffs.push_back(monomial(a * t.coeff(), t.var())); - - return ret; -} -} diff --git a/src/util/lp/general_matrix.h b/src/util/lp/general_matrix.h index 73521424c..715f2cb08 100644 --- a/src/util/lp/general_matrix.h +++ b/src/util/lp/general_matrix.h @@ -18,6 +18,7 @@ Revision History: --*/ #pragma once +#include namespace lp { class general_matrix { // fields @@ -72,12 +73,18 @@ public: #ifdef Z3DEBUG void print(std::ostream & out, unsigned blanks = 0) const { - print_matrix(m_data, out, blanks); + unsigned m = row_count(); + unsigned n = column_count(); + general_matrix g(m, n); + for (unsigned i = 0; i < m; i++) + for (unsigned j = 0; j < n; j++) + g[i][j] = (*this)[i][j]; + print_matrix(g.m_data, out, blanks); } void print(std::ostream & out, const char * ss) const { std::string s(ss); out << s; - print(out, s.size()); + print(out, static_cast(s.size())); } void print_submatrix(std::ostream & out, unsigned k, unsigned blanks = 0) const { diff --git a/src/util/lp/hnf.h b/src/util/lp/hnf.h index aa31d5d01..7a1b448c7 100644 --- a/src/util/lp/hnf.h +++ b/src/util/lp/hnf.h @@ -20,7 +20,6 @@ $1$ at $i$-th position. Then we need to find the row vector $e_iU^{-1}=t$. Notic We find $e_iH^{-1} = f$ by solving $e_i = fH$ and then $fA$ gives us $t$. Author: - Nikolaj Bjorner (nbjorner) Lev Nachmanson (levnach) Revision History: @@ -104,15 +103,13 @@ void extended_gcd_minimal_uv(const mpq & a, const mpq & b, mpq & d, mpq & u, mpq -template bool prepare_pivot_for_lower_triangle(M &m, unsigned r, svector & basis_rows) { - lp_assert(m.row_count() <= m.column_count()); +template bool prepare_pivot_for_lower_triangle(M &m, unsigned r) { for (unsigned i = r; i < m.row_count(); i++) { for (unsigned j = r; j < m.column_count(); j++) { if (!is_zero(m[i][j])) { if (i != r) { m.transpose_rows(i, r); } - basis_rows.push_back(i); if (j != r) { m.transpose_columns(j, r); } @@ -123,7 +120,8 @@ template bool prepare_pivot_for_lower_triangle(M &m, unsigned r, sv return false; } -template void pivot_column_non_fractional(M &m, unsigned & r) { +template void pivot_column_non_fractional(M &m, unsigned r) { + lp_assert(!is_zero(m[r][r])); lp_assert(m.row_count() <= m.column_count()); for (unsigned j = r + 1; j < m.column_count(); j++) { for (unsigned i = r + 1; i < m.row_count(); i++) { @@ -134,37 +132,36 @@ template void pivot_column_non_fractional(M &m, unsigned & r) { lp_assert(is_int(m[i][j])); } } - + // debug + for (unsigned k = r + 1; k < m.row_count(); k++) { + m[k][r] = zero_of_type(); + } + } // returns the rank of the matrix -template void to_lower_triangle_non_fractional(M &m, svector & basis_rows ) { +template unsigned to_lower_triangle_non_fractional(M &m) { lp_assert(m.row_count() <= m.column_count()); unsigned i = 0; - for (; i < m.row_count() - 1; i++) { - if (!prepare_pivot_for_lower_triangle(m, i, basis_rows)) { - return; + for (; i < m.row_count(); i++) { + if (!prepare_pivot_for_lower_triangle(m, i)) { + return i; } pivot_column_non_fractional(m, i); } - lp_assert(i == m.row_count() - 1); - // go over the last row and try to find a non-zero in the row to the right of diagonal - for (unsigned j = i; j < m.column_count(); j++) { - if (!is_zero(m[i][j])) { - basis_rows.push_back(i); - break; - } - } + lp_assert(i == m.row_count()); + return i; } template mpq gcd_of_row_starting_from_diagonal(const M& m, unsigned i) { mpq g = zero_of_type(); unsigned j = i; - for (; j < m.column_count() && is_zero(j); j++) { + for (; j < m.column_count() && is_zero(g); j++) { const auto & t = m[i][j]; if (!is_zero(t)) - g = t; + g = abs(t); } + lp_assert(!is_zero(g)); for (; j < m.column_count(); j++) { const auto & t = m[i][j]; if (!is_zero(t)) @@ -183,10 +180,15 @@ template mpq determinant_of_rectangular_matrix(const M& m, svector< // m[r-1][r-1], m[r-1][r], ..., m[r-1]m[m.column_count() - 1] give the determinants of all minors of rank r. // The gcd of these minors is the return value auto mc = m; - to_lower_triangle_non_fractional(mc, basis_rows); - if (basis_rows.size() == 0) + unsigned rank = to_lower_triangle_non_fractional(mc); + if (rank == 0) return one_of_type(); - return gcd_of_row_starting_from_diagonal(mc, basis_rows.size() - 1); + + for (unsigned i = 0; i < rank; i++) { + basis_rows.push_back(mc.adjust_row(i)); + } + TRACE("hnf_calc", tout << "basis_rows = "; print_vector(basis_rows, tout); mc.print(tout, "mc = ");); + return gcd_of_row_starting_from_diagonal(mc, rank - 1); } template mpq determinant(const M& m) { @@ -499,9 +501,7 @@ private: m_W[k][j] -= u * m_W[k][m_i]; // m_W[k][j] = mod_R_balanced(m_W[k][j]); } - } - - + } bool is_unit_matrix(const M& u) const { unsigned m = u.row_count(); @@ -597,8 +597,8 @@ public: hnf(M & A, const mpq & d) : #ifdef Z3DEBUG m_H(A), -#endif m_A_orig(A), +#endif m_W(A), m_buffer(std::max(A.row_count(), A.column_count())), m_m(A.row_count()), diff --git a/src/util/lp/hnf_cutter.h b/src/util/lp/hnf_cutter.h index 8677356f0..c2a8672f7 100644 --- a/src/util/lp/hnf_cutter.h +++ b/src/util/lp/hnf_cutter.h @@ -11,7 +11,6 @@ Abstract: Author: Lev Nachmanson (levnach) - Nikolaj Bjorner (nbjorner) Revision History: @@ -33,13 +32,21 @@ class hnf_cutter { vector m_right_sides; unsigned m_row_count; unsigned m_column_count; - std::function m_random_next; + lp_settings & m_settings; public: - hnf_cutter(std::function random) : m_random_next(random) {} + hnf_cutter(lp_settings & settings) : m_row_count(0), m_column_count(0), m_settings(settings) {} + unsigned row_count() const { + return m_row_count; + } + + const vector& terms() const { return m_terms; } + const vector & right_sides() const { return m_right_sides; } void clear() { + // m_A will be filled from scratch in init_matrix_A m_var_register.clear(); m_terms.clear(); + m_right_sides.clear(); m_row_count = m_column_count = 0; } void add_term(const lar_term* t, const mpq &rs) { @@ -95,15 +102,14 @@ public: int ret = -1; int n = 0; for (int i = 0; i < static_cast(b.size()); i++) { - if (!is_int(b[i])) { - if (n == 0 ) { - lp_assert(ret == -1); - n = 1; + if (is_int(b[i])) continue; + if (n == 0 ) { + lp_assert(ret == -1); + n = 1; + ret = i; + } else { + if (m_settings.random_next() % (++n) == 0) { ret = i; - } else { - if (m_random_next() % (++n) == 0) { - ret = i; - } } } } @@ -144,11 +150,28 @@ public: t.add_monomial(row[j], m_var_register.local_var_to_user_var(j)); } } - - lia_move create_cut(lar_term& t, mpq& k, explanation& ex, bool & upper) { +#ifdef Z3DEBUG + vector transform_to_local_columns(const vector & x) const { + vector ret; + lp_assert(m_column_count <= m_var_register.size()); + for (unsigned j = 0; j < m_column_count;j++) { + lp_assert(is_zero(x[m_var_register.local_var_to_user_var(j)].y)); + ret.push_back(x[m_var_register.local_var_to_user_var(j)].x); + } + return ret; + } +#endif + lia_move create_cut(lar_term& t, mpq& k, explanation& ex, bool & upper + #ifdef Z3DEBUG + , + const vector & x0 + #endif + ) { init_matrix_A(); svector basis_rows; mpq d = hnf_calc::determinant_of_rectangular_matrix(m_A, basis_rows); + if (m_settings.get_cancel_flag()) + return lia_move::undef; if (basis_rows.size() < m_A.row_count()) m_A.shrink_to_rank(basis_rows); @@ -156,10 +179,10 @@ public: // general_matrix A_orig = m_A; vector b = create_b(basis_rows); - vector bcopy = b; + lp_assert(m_A * x0 == b); + // vector bcopy = b; find_h_minus_1_b(h.W(), b); - - lp_assert(bcopy == h.W().take_first_n_columns(b.size()) * b); + // lp_assert(bcopy == h.W().take_first_n_columns(b.size()) * b); int cut_row = find_cut_row_index(b); if (cut_row == -1) { return lia_move::undef; @@ -181,7 +204,6 @@ public: vector row(m_A.column_count()); get_ei_H_minus_1(cut_row, h.W(), row); vector f = row * m_A; - fill_term(f, t); k = floor(b[cut_row]); upper = true; diff --git a/src/util/lp/int_solver.cpp b/src/util/lp/int_solver.cpp index 94121d27f..73fdacbaf 100644 --- a/src/util/lp/int_solver.cpp +++ b/src/util/lp/int_solver.cpp @@ -5,9 +5,9 @@ #include "util/lp/int_solver.h" #include "util/lp/lar_solver.h" -#include "util/lp/chase_cut_solver.h" #include "util/lp/lp_utils.h" #include +#include "util/lp/monomial.h" namespace lp { @@ -32,22 +32,15 @@ void int_solver::trace_inf_rows() const { ); } -bool int_solver::all_columns_are_bounded() const { - for (unsigned j = 0; j < m_lar_solver->column_count(); j++) - if (m_lar_solver->column_is_bounded(j) == false) - return false; - return true; -} - bool int_solver::has_inf_int() const { return m_lar_solver->has_inf_int(); } int int_solver::find_inf_int_base_column() { - unsigned inf_int_count; + unsigned inf_int_count = 0; int j = find_inf_int_boxed_base_column_with_smallest_range(inf_int_count); - if (j != -1) - return j; + if (j != -1) + return j; if (inf_int_count == 0) return -1; unsigned k = random() % inf_int_count; @@ -55,22 +48,17 @@ int int_solver::find_inf_int_base_column() { } int int_solver::get_kth_inf_int(unsigned k) const { - unsigned inf_int_count = 0; - for (unsigned j : m_lar_solver->r_basis()) { - if (! column_is_int_inf(j) ) - continue; - if (inf_int_count++ == k) + for (unsigned j : m_lar_solver->r_basis()) + if (column_is_int_inf(j) && k-- == 0) return j; - } lp_assert(false); return -1; } int int_solver::find_inf_int_nbasis_column() const { for (unsigned j : m_lar_solver->r_nbasis()) - if (! column_is_int_inf(j) ) - return j; - + if (!column_is_int_inf(j)) + return j; return -1; } @@ -80,7 +68,7 @@ int int_solver::find_inf_int_boxed_base_column_with_smallest_range(unsigned & in mpq range; mpq new_range; mpq small_range_thresold(1024); - unsigned n; + unsigned n = 0; lar_core_solver & lcs = m_lar_solver->m_mpq_lar_core_solver; for (unsigned j : m_lar_solver->r_basis()) { @@ -93,24 +81,14 @@ int int_solver::find_inf_int_boxed_base_column_with_smallest_range(unsigned & in new_range = lcs.m_r_upper_bounds()[j].x - lcs.m_r_lower_bounds()[j].x; if (new_range > small_range_thresold) continue; - if (result == -1) { + if (result == -1 || new_range < range) { result = j; range = new_range; n = 1; - continue; } - if (new_range < range) { - n = 1; + else if (new_range == range && settings().random_next() % (++n) == 0) { + lp_assert(n > 1); result = j; - range = new_range; - continue; - } - if (new_range == range) { - lp_assert(n >= 1); - if (settings().random_next() % (++n) == 0) { - result = j; - continue; - } } } return result; @@ -272,18 +250,10 @@ void int_solver::gomory_cut_adjust_t_and_k(vector> & po bool int_solver::current_solution_is_inf_on_cut() const { const auto & x = m_lar_solver->m_mpq_lar_core_solver.m_r_x; impq v = m_t->apply(x); - mpq sign = !(*m_upper) ? one_of_type() : -one_of_type(); - TRACE("current_solution_is_inf_on_cut", - if (is_pos(sign)) { - tout << "v = " << v << " k = " << (*m_k) << std::endl; - if (v <=(*m_k)) { - tout << "v <= k - it should not happen!\n"; - } - } else { - if (v >= (*m_k)) { - tout << "v > k - it should not happen!\n"; - } - } + mpq sign = *m_upper ? one_of_type() : -one_of_type(); + CTRACE("current_solution_is_inf_on_cut", v * sign <= (*m_k) * sign, + tout << "m_upper = " << *m_upper << std::endl; + tout << "v = " << v << ", k = " << (*m_k) << std::endl; ); return v * sign > (*m_k) * sign; } @@ -362,8 +332,7 @@ lia_move int_solver::mk_gomory_cut( unsigned inf_col, const row_strip & row a.neg(); if (is_real(x_j)) real_case_in_gomory_cut(a, x_j, f_0, one_min_f_0); - else { - if (a.is_int()) continue; // f_j will be zero and no monomial will be added + else if (!a.is_int()) { // f_j will be zero and no monomial will be added some_int_columns = true; int_case_in_gomory_cut(a, x_j, lcm_den, f_0, one_min_f_0); } @@ -392,14 +361,14 @@ int int_solver::find_free_var_in_gomory_row(const row_strip& row) { lia_move int_solver::proceed_with_gomory_cut(unsigned j) { const row_strip& row = m_lar_solver->get_row(row_of_basic_column(j)); - int free_j = find_free_var_in_gomory_row(row); - if (free_j != -1) - return lia_move::undef; - if (!is_gomory_cut_target(row)) { - return lia_move::undef; - } - *m_upper = false; + if (-1 != find_free_var_in_gomory_row(row)) + return lia_move::undef; + + if (!is_gomory_cut_target(row)) + return lia_move::undef; + + *m_upper = true; return mk_gomory_cut(j, row); } @@ -420,29 +389,7 @@ unsigned int_solver::row_of_basic_column(unsigned j) const { // } -typedef chase_cut_solver::monomial mono; - -// it produces an inequality coeff*x <= rs -template -void int_solver::get_int_coeffs_from_constraint(const lar_base_constraint* c, - vector& coeffs, T & rs) { - lp_assert(c->m_kind != EQ); // it is not implemented, we need to create two inequalities in this case - int sign = ((int)c->m_kind > 0) ? -1 : 1; - vector> lhs = c->get_left_side_coefficients(); - T den = denominator(c->m_right_side); - for (auto & kv : lhs) { - lp_assert(!is_term(kv.second)); - lp_assert(is_int(kv.second)); // not implemented for real vars! - den = lcm(den, denominator(kv.first)); - } - lp_assert(den > 0); - for (auto& kv : lhs) { - coeffs.push_back(mono(den * kv.first * sign, kv.second)); - } - rs = den * c->m_right_side * sign; - if (kind_is_strict(c->m_kind)) - rs--; -} +typedef monomial mono; // this will allow to enable and disable tracking of the pivot rows @@ -462,36 +409,6 @@ struct pivoted_rows_tracking_control { } }; -void int_solver::copy_explanations_from_chase_cut_solver() { - TRACE("propagate_and_backjump_step_int", - for (unsigned j: m_chase_cut_solver.m_explanation) - m_lar_solver->print_constraint(m_lar_solver->constraints()[j], tout);); - - for (unsigned j : m_chase_cut_solver.m_explanation) { - m_ex->push_justification(j); - } - m_chase_cut_solver.m_explanation.clear(); -} - -void int_solver::copy_values_from_chase_cut_solver() { - for (unsigned j = 0; j < m_lar_solver->A_r().column_count() && j < m_chase_cut_solver.number_of_vars(); j++) { - if (!m_chase_cut_solver.var_is_active(j)) - continue; - if (!is_int(j)) { - continue; - } - m_lar_solver->m_mpq_lar_core_solver.m_r_x[j] = m_chase_cut_solver.var_value(j); - lp_assert(m_lar_solver->column_value_is_int(j)); - } -} - -void int_solver::catch_up_in_adding_constraints_to_chase_cut_solver() { - lp_assert(m_chase_cut_solver.number_of_asserts() <= m_lar_solver->constraints().size()); - for (unsigned j = m_chase_cut_solver.number_of_asserts(); j < m_lar_solver->constraints().size(); j++) { - add_constraint_to_chase_cut_solver(j, m_lar_solver->constraints()[j]); - } -} - impq int_solver::get_cube_delta_for_term(const lar_term& t) const { if (t.size() == 2) { bool seen_minus = false; @@ -541,9 +458,9 @@ bool int_solver::tighten_terms_for_cube() { return true; } -bool int_solver::find_cube() { +lia_move int_solver::find_cube() { if (m_branch_cut_counter % settings().m_int_find_cube_period != 0) - return false; + return lia_move::undef; settings().st().m_cube_calls++; TRACE("cube", @@ -551,26 +468,25 @@ bool int_solver::find_cube() { display_column(tout, j); m_lar_solver->print_terms(tout); ); - m_lar_solver->push(); - if(!tighten_terms_for_cube()) { - m_lar_solver->pop(); - return false; + + lar_solver::scoped_push _sp(*m_lar_solver); + if (!tighten_terms_for_cube()) { + return lia_move::undef; } lp_status st = m_lar_solver->find_feasible_solution(); if (st != lp_status::FEASIBLE && st != lp_status::OPTIMAL) { TRACE("cube", tout << "cannot find a feasiblie solution";); - m_lar_solver->pop(); + _sp.pop(); move_non_basic_columns_to_bounds(); find_feasible_solution(); - lp_assert(m_chase_cut_solver.cancel() || is_feasible()); // it can happen that we found an integer solution here - return !m_lar_solver->r_basis_has_inf_int(); + return !m_lar_solver->r_basis_has_inf_int()? lia_move::sat: lia_move::undef; } - m_lar_solver->pop(); + _sp.pop(); m_lar_solver->round_to_integer_solution(); - lp_assert(m_chase_cut_solver.cancel() || is_feasible()); - return true; + settings().st().m_cube_success++; + return lia_move::sat; } void int_solver::find_feasible_solution() { @@ -589,55 +505,26 @@ lia_move int_solver::run_gcd_test() { return lia_move::undef; } -lia_move int_solver::call_chase_cut_solver() { - if ((m_branch_cut_counter) % settings().m_int_chase_cut_solver_period != 0 || !all_columns_are_bounded()) - return lia_move::undef; - TRACE("check_main_int", tout<<"chase_cut_solver";); - catch_up_in_adding_constraints_to_chase_cut_solver(); - auto check_res = m_chase_cut_solver.check(); - settings().st().m_chase_cut_solver_calls++; - switch (check_res) { - case chase_cut_solver::lbool::l_false: - copy_explanations_from_chase_cut_solver(); - settings().st().m_chase_cut_solver_false++; - return lia_move::conflict; - case chase_cut_solver::lbool::l_true: - settings().st().m_chase_cut_solver_true++; - copy_values_from_chase_cut_solver(); - lp_assert(m_lar_solver->all_constraints_hold()); - return lia_move::sat; - case chase_cut_solver::lbool::l_undef: - settings().st().m_chase_cut_solver_undef++; - if (m_chase_cut_solver.try_getting_cut(*m_t, *m_k, m_lar_solver->m_mpq_lar_core_solver.m_r_x)) { - m_lar_solver->subs_term_columns(*m_t); - TRACE("chase_cut_solver_cuts", - tout<<"precut from chase_cut_solver:"; m_lar_solver->print_term(*m_t, tout); tout << " <= " << *m_k << std::endl;); - - return lia_move::cut; - } - default: - return lia_move::undef; - } -} - lia_move int_solver::gomory_cut() { - TRACE("check_main_int", tout << "gomory";); + if ((m_branch_cut_counter) % settings().m_int_gomory_cut_period != 0) + return lia_move::undef; + if (move_non_basic_columns_to_bounds()) { - lp_status st = m_lar_solver->find_feasible_solution(); - if (st != lp_status::FEASIBLE && st != lp_status::OPTIMAL) { - TRACE("arith_int", tout << "give_up\n";); - return lia_move::undef; - } +#if Z3DEBUG + lp_status st = +#endif + m_lar_solver->find_feasible_solution(); +#if Z3DEBUG + lp_assert(st == lp_status::FEASIBLE || st == lp_status::OPTIMAL); +#endif } + int j = find_inf_int_base_column(); if (j == -1) { j = find_inf_int_nbasis_column(); return j == -1? lia_move::sat : create_branch_on_column(j); } - lia_move r = proceed_with_gomory_cut(j); - if (r != lia_move::undef) - return r; - return create_branch_on_column(j); + return proceed_with_gomory_cut(j); } @@ -645,39 +532,49 @@ bool int_solver::try_add_term_to_A_for_hnf(unsigned i) { mpq rs; const lar_term* t = m_lar_solver->terms()[i]; for (const auto & p : *t) { - if (!is_int(p.var())) + if (!is_int(p.var())) { + lp_assert(false); return false; // todo : the mix case! + } } - if (m_lar_solver->get_equality_for_term_on_corrent_x(i, rs)) { + bool has_bounds; + if (m_lar_solver->get_equality_and_right_side_for_term_on_corrent_x(i, rs, has_bounds)) { m_hnf_cutter.add_term(t, rs); return true; - } else { - return false; } + return !has_bounds; } bool int_solver::hnf_matrix_is_empty() const { return true; } -bool int_solver::prepare_matrix_A_for_hnf_cut() { +bool int_solver::init_terms_for_hnf_cut() { m_hnf_cutter.clear(); for (unsigned i = 0; i < m_lar_solver->terms().size(); i++) { - bool r = try_add_term_to_A_for_hnf(i); - if (!r && settings().hnf_cutter_exit_if_x_is_not_on_bound_or_mixed ) - return false; + try_add_term_to_A_for_hnf(i); } - return true; + return m_hnf_cutter.row_count() < settings().limit_on_rows_for_hnf_cutter; } lia_move int_solver::make_hnf_cut() { - if (!prepare_matrix_A_for_hnf_cut()) { + if (!init_terms_for_hnf_cut()) { return lia_move::undef; } settings().st().m_hnf_cutter_calls++; - lia_move r = m_hnf_cutter.create_cut(*m_t, *m_k, *m_ex, *m_upper); + TRACE("hnf_cut", tout << "settings().st().m_hnf_cutter_calls = " << settings().st().m_hnf_cutter_calls;); +#ifdef Z3DEBUG + vector x0 = m_hnf_cutter.transform_to_local_columns(m_lar_solver->m_mpq_lar_core_solver.m_r_x); +#endif + lia_move r = m_hnf_cutter.create_cut(*m_t, *m_k, *m_ex, *m_upper +#ifdef Z3DEBUG + , x0 +#endif + ); CTRACE("hnf_cut", r == lia_move::cut, tout<< "cut:"; m_lar_solver->print_term(*m_t, tout); tout << " <= " << *m_k << std::endl;); - if (r == lia_move::cut) + if (r == lia_move::cut) { + lp_assert(current_solution_is_inf_on_cut()); settings().st().m_hnf_cuts++; + } return r; } @@ -689,43 +586,42 @@ lia_move int_solver::hnf_cut() { } lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex, bool & upper) { - if (!has_inf_int()) - return lia_move::sat; + if (!has_inf_int()) return lia_move::sat; + m_t = &t; m_k = &k; m_ex = &ex; m_upper = &upper; - if (run_gcd_test() == lia_move::conflict) - return lia_move::conflict; - + lia_move r = run_gcd_test(); + if (r != lia_move::undef) return r; + pivoted_rows_tracking_control pc(m_lar_solver); + if(settings().m_int_pivot_fixed_vars_from_basis) m_lar_solver->pivot_fixed_vars_from_basis(); - if (patch_nbasic_columns() == lia_move::sat) - return lia_move::sat; + r = patch_nbasic_columns(); + if (r != lia_move::undef) return r; ++m_branch_cut_counter; - if (find_cube()){ - settings().st().m_cube_success++; - return lia_move::sat; - } - lia_move r = call_chase_cut_solver(); - if (r != lia_move::undef) - return r; - - r = hnf_cut(); - if (r != lia_move::undef) - return r; + r = find_cube(); + if (r != lia_move::undef) return r; - if ((m_branch_cut_counter) % settings().m_int_gomory_cut_period == 0) { - return gomory_cut(); - } - int j = find_inf_int_base_column(); - if (j == -1) { - j = find_inf_int_nbasis_column(); - if (j == -1) - return lia_move::sat; - } - return create_branch_on_column(j); + r = hnf_cut(); + if (r != lia_move::undef) return r; + + r = gomory_cut(); + return (r == lia_move::undef)? branch_or_sat() : r; +} + +lia_move int_solver::branch_or_sat() { + int j = find_any_inf_int_column_basis_first(); + return j == -1? lia_move::sat : create_branch_on_column(j); +} + +int int_solver::find_any_inf_int_column_basis_first() { + int j = find_inf_int_base_column(); + if (j != -1) + return j; + return find_inf_int_nbasis_column(); } bool int_solver::move_non_basic_column_to_bounds(unsigned j) { @@ -751,7 +647,7 @@ bool int_solver::move_non_basic_column_to_bounds(unsigned j) { if (val != lcs.m_r_upper_bounds()[j]) { set_value_for_nbasic_column(j, lcs.m_r_upper_bounds()[j]); return true; - } + } break; default: if (is_int(j) && !val.is_int()) { @@ -845,6 +741,7 @@ void int_solver::patch_nbasic_column(unsigned j) { tout << "patching with 0\n";); } } + lia_move int_solver::patch_nbasic_columns() { settings().st().m_patches++; lp_assert(is_feasible()); @@ -931,6 +828,7 @@ bool int_solver::gcd_test_for_row(static_matrix> & A, uns return true; } + void int_solver::add_to_explanation_from_fixed_or_boxed_column(unsigned j) { constraint_index lc, uc; m_lar_solver->get_bound_constraint_witnesses_for_column(j, lc, uc); @@ -948,10 +846,8 @@ void int_solver::fill_explanation_from_fixed_columns(const row_strip & row) bool int_solver::gcd_test() { auto & A = m_lar_solver->A_r(); // getting the matrix for (unsigned i = 0; i < A.row_count(); i++) - if (!gcd_test_for_row(A, i)) { - return false; - } - + if (!gcd_test_for_row(A, i)) + return false; return true; } @@ -1022,19 +918,15 @@ linear_combination_iterator * int_solver::get_column_iterator(unsigned j) { } */ + int_solver::int_solver(lar_solver* lar_slv) : m_lar_solver(lar_slv), m_branch_cut_counter(0), - m_chase_cut_solver([this](unsigned j) {return m_lar_solver->get_column_name(j);}, - [this](unsigned j, std::ostream &o) {m_lar_solver->print_constraint(j, o);}, - [this]() {return m_lar_solver->A_r().column_count();}, - [this](unsigned j) {return get_value(j);}, - settings()), - m_hnf_cutter([this](){ return settings().random_next(); }) -{ + m_hnf_cutter(settings()) { m_lar_solver->set_int_solver(this); } + bool int_solver::has_low(unsigned j) const { switch (m_lar_solver->m_mpq_lar_core_solver.m_column_types()[j]) { case column_type::fixed: @@ -1067,6 +959,7 @@ void set_lower(impq & l, } } + void set_upper(impq & u, bool & inf_u, impq const & v) { @@ -1386,22 +1279,6 @@ bool int_solver::is_term(unsigned j) const { return m_lar_solver->column_corresponds_to_term(j); } -void int_solver::add_constraint_to_chase_cut_solver(unsigned ci, const lar_base_constraint * c) { - vector coeffs; - mpq rs; - get_int_coeffs_from_constraint(c, coeffs, rs); - m_chase_cut_solver.add_ineq(coeffs, -rs, ci); -} - -void int_solver::pop(unsigned k) { - m_chase_cut_solver.pop_trail(k); - while (m_chase_cut_solver.number_of_asserts() > m_lar_solver->constraints().size()) - m_chase_cut_solver.pop_last_assert(); - m_chase_cut_solver.pop_constraints(); -} - -void int_solver::push() { m_chase_cut_solver.push(); } - unsigned int_solver::column_count() const { return m_lar_solver->column_count(); } } diff --git a/src/util/lp/int_solver.h b/src/util/lp/int_solver.h index 210f350e1..c2aa32360 100644 --- a/src/util/lp/int_solver.h +++ b/src/util/lp/int_solver.h @@ -22,7 +22,6 @@ Revision History: #include "util/lp/static_matrix.h" #include "util/lp/int_set.h" #include "util/lp/lar_term.h" -#include "util/lp/chase_cut_solver.h" #include "util/lp/lar_constraints.h" #include "util/lp/hnf_cutter.h" #include "util/lp/lia_move.h" @@ -38,10 +37,9 @@ struct lp_constraint; class int_solver { public: // fields - lar_solver * m_lar_solver; + lar_solver *m_lar_solver; unsigned m_branch_cut_counter; - chase_cut_solver m_chase_cut_solver; - lar_term* m_t; // the term to return in the cut + lar_term *m_t; // the term to return in the cut mpq *m_k; // the right side of the cut explanation *m_ex; // the conflict explanation bool *m_upper; // we have a cut m_t*x <= k if m_upper is true nad m_t*x >= k otherwise @@ -96,6 +94,8 @@ private: const impq & get_value(unsigned j) const; bool column_is_int_inf(unsigned j) const; void trace_inf_rows() const; + lia_move branch_or_sat(); + int find_any_inf_int_column_basis_first(); int find_inf_int_base_column(); int find_inf_int_boxed_base_column_with_smallest_range(unsigned&); int get_kth_inf_int(unsigned) const; @@ -139,20 +139,10 @@ private: unsigned random(); bool has_inf_int() const; lia_move create_branch_on_column(int j); - void catch_up_in_adding_constraints_to_chase_cut_solver(); public: - template - void fill_chase_cut_solver_vars(); - template - void get_int_coeffs_from_constraint(const lar_base_constraint* c, vector& coeff, T & rs); bool is_term(unsigned j) const; - void add_constraint_to_chase_cut_solver(unsigned,const lar_base_constraint*); - void copy_explanations_from_chase_cut_solver(); - void pop(unsigned); - void push(); - void copy_values_from_chase_cut_solver(); bool left_branch_is_more_narrow_than_right(unsigned); - bool find_cube(); + lia_move find_cube(); bool tighten_terms_for_cube(); bool tighten_term_for_cube(unsigned); unsigned column_count() const; @@ -161,11 +151,10 @@ public: void find_feasible_solution(); int find_inf_int_nbasis_column() const; lia_move run_gcd_test(); - lia_move call_chase_cut_solver(); lia_move gomory_cut(); lia_move hnf_cut(); lia_move make_hnf_cut(); - bool prepare_matrix_A_for_hnf_cut(); + bool init_terms_for_hnf_cut(); bool hnf_matrix_is_empty() const; bool try_add_term_to_A_for_hnf(unsigned term_index); }; diff --git a/src/util/lp/integer_domain.h b/src/util/lp/integer_domain.h deleted file mode 100644 index 44f3c7f9e..000000000 --- a/src/util/lp/integer_domain.h +++ /dev/null @@ -1,1050 +0,0 @@ -/* - Copyright (c) 2017 Microsoft Corporation - Author: Lev Nachmanson -*/ -#pragma once -#include -#include "util/trace.h" -#include "util/lp/stacked_value.h" -#include "util/lp/stacked_map.h" -namespace lp { -enum class endpoint_kind { START, END, STEND }; -// represents the set of disjoint intervals of integer number -template -class integer_domain { -#ifdef Z3DEBUG - // std::set m_domain; -#endif - struct endpoint { int m_start_expl; int m_end_expl; - endpoint() : m_start_expl(-1), m_end_expl(-1) {} - endpoint(int s, int e) : m_start_expl(s), m_end_expl(e) {} - endpoint_kind kind() const { - lp_assert(m_start_expl != -1 || m_end_expl != -1); - if (m_end_expl == -1) { - return endpoint_kind::START; - } - if (m_start_expl == -1) { - return endpoint_kind::END; - } - - return endpoint_kind::STEND; - } - bool operator==(const endpoint& e) const { - return m_start_expl == e.m_start_expl && m_end_expl == e.m_end_expl; - } - bool operator!=(const endpoint & e) const { return !(*this == e); } - - void print(std::ostream & out) const { - if (m_start_expl != -1 && m_end_expl != -1) - out << "(" << m_start_expl << ", " << m_end_expl << ")"; - else { - if (m_start_expl != -1) { - out << "(" << m_start_expl << ")"; - } - else if (m_end_expl != -1) { - out << "(" << m_end_expl << ")"; - } - } - } - }; - stacked_map m_endpoints; - stacked_value m_empty; - typedef typename std::map::iterator iter; - typedef typename std::map::const_iterator const_iter; - typedef typename std::map::reverse_iterator riter; - -public: - // the default constructor creates a set containing all integer numbers - integer_domain() : integer_domain(false) {} - - - // if is_empty = false then the constructor creates a set containing all integer numbers, - // otherwise it creates an empty set - integer_domain(bool is_empty) : m_empty(is_empty) { -#if Z3DEBUG - // if (!is_empty) { - // for (int i = 0; i <= 100; i++) - // m_domain.insert(i); - // } -#endif - } - - void init_to_contain_all() { - m_empty = false; - m_endpoints.clear(); -#if Z3DEBUG - // for (int i = 0; i <= 100; i++) - // m_domain.insert(i); -#endif - } - - // copy constructor - integer_domain(const integer_domain & t) : -#if Z3DEBUG - // m_domain(t.m_domain), -#endif - m_endpoints(t.m_endpoints), - m_empty(t.m_empty) - { - } - - // needed for debug only - void restore_domain() { -#if Z3DEBUG - // for (int i = 0; i <= 100; i++) - // if (contains(i)) - // m_domain.insert(i); - // else - // m_domain.erase(i); -#endif - } - - bool operator==(const integer_domain & t) { - return m_empty == t.m_empty && m_endpoints() == t.m_endpoints(); - } - bool contains_all() const { - return !m_empty && m_endpoints.empty(); - } - - bool contains(const T & x) const { - if (contains_all()) - return true; - bool neg_inf; - const_iter l; - bool found_left_point = get_left_point(x, neg_inf, l); - if (!found_left_point) - return has_neg_inf(); - if (neg_inf) - return true; - if (pos(l) == x) - return true; - return is_proper_start(l); - } - - void handle_right_point_in_union(const_iter &r, const T &y) { - if (pos(r) == y) { - if (is_proper_start(r)) - erase(r); - else - set_end(y); - } - else if (pos(r) == y + 1) { - if (is_proper_start(r)) { - erase(r); - } - else { - set_end(r->first); - } - } - else if (!is_proper_end(r)) - set_end(y); - - lp_assert(is_correct()); - } - void handle_left_point_in_union(const_iter& l, const T &x, const T & y) { - if (pos(l) == x || pos(l) + 1 == x) { - if (is_proper_end(l)) { - l++; - erase(std::prev(l)); - } - else { - if (is_one_point_interval(l)) { - set_start(pos(l)); - } - } - } - else { - if (!is_proper_start(l)) { - set_start(x); - } - } - - while (l!= m_endpoints.end() && pos(l) <= x) - l++; - remove_from_the_left(y, l); - } - - void unite_with_interval(const T& x, const T& y) { - lp_assert(false); - // TRACE("disj_intervals", tout << "unite_with_interval(" << x << ", " << y << ")\n";); - // #if Z3DEBUG - // // for (int i = std::max(x, 0); i <= std::min(100, y); i++) - // // m_domain.insert(i); - // #endif - - // lp_assert(x <= y); - // if (x == y) { - // unite_with_one_point_interval(x); - // return; - // } - - // const_iter l, r; - // bool neg_inf, pos_inf; - // bool found_left_point = get_left_point(x, neg_inf, l); - // bool found_right_point = get_right_point(y, pos_inf, r); - // m_empty = false; - - // if (!found_left_point) { - // if (!found_right_point) { - // m_endpoints.clear(); - // set_start(x); - // set_end(y); - // return; - // } - // // found_right_point is true - // if (pos_inf) { - // m_endpoints.clear(); - // set_start(x); - // return; - // } - // remove_from_the_left(y); - - // if (pos(m_endpoints.begin()) == y || pos(m_endpoints.begin()) == y + 1) { - // if (is_proper_start(m_endpoints.begin())) - // m_endpoints.erase(m_endpoints.begin()); - // else - // set_end(pos(m_endpoints.begin())); - // set_start(x); - // } - // else { - // lp_assert(pos(m_endpoints.begin()) > y + 1); - // if (is_start(m_endpoints.begin())) - // set_end(y); - // set_start(x); - // } - // return; - // } - - // lp_assert(found_left_point); - // if (!found_right_point) { - // bool neg_inf = has_neg_inf(); - // remove_from_the_right(x); - // if (m_endpoints.empty()) { - // if (!neg_inf) - // set_start_end(x, y); - // else - // set_end(y); - // return; - // } - // if (pos(m_endpoints.rbegin()) == x || pos(m_endpoints.rbegin()) == x - 1) { - // if (is_proper_end(m_endpoints.rbegin())) { - // m_endpoints.erase(m_endpoints.rbegin()); - // } - // else if (is_one_point_interval(m_endpoints.rbegin())) { - // set_start(pos(m_endpoints.rbegin())); - // } - // set_end(y); - // } - // else { - // if (is_end(m_endpoints.rbegin())) { - // set_start(x); - // set_end(y); - // } - // else { - // set_end(y); - // } - // } - // return; - // } - - // // found_right_point and found_left_point - // if (!neg_inf) - // handle_left_point_in_union(l, x, y); - // else { - // remove_from_the_left(y); - // } - // if (!pos_inf) - // handle_right_point_in_union(r, y); - // else - // remove_from_the_right(x); - } - - bool has_pos_inf() const { - if (m_empty) - return false; - - if (m_endpoints.empty()) - return true; - - lp_assert(m_endpoints.rbegin() != m_endpoints.rend()); - return m_endpoints.rbegin()->second.kind() == endpoint_kind::START; - } - - bool has_neg_inf() const { - if (m_empty) - return false; - - if (m_endpoints.empty()) - return true; - auto it = m_endpoints.begin(); - return is_proper_end(it->second.kind());//m_endpoints.begin()); - } - - bool is_correct() const { - if (m_empty) { - if (m_endpoints.size() > 0) { - TRACE("disj_intervals", tout << "is empty is true but m_endpoints.size() = " << m_endpoints.size() << std::endl;); - return false; - } - return true; - } - bool expect_end; - bool prev = false; - T prev_x; - for (auto t : m_endpoints()) { - if (prev && ((expect_end && !is_end(t.second)) || (!expect_end && !is_start(t.second)))) { - TRACE("disj_intervals", tout << "x = " << t.first << "\n";); - if (expect_end) { - TRACE("disj_intervals", tout << "expecting an interval end\n";); - } else { - TRACE("disj_intervals", tout << "expecting an interval start\n";); - } - return false; - } - - if (prev) { - if (t.first - prev_x <= 1 && !expect_end) { - TRACE("disj_intervals", tout << "the sequence is not increasing or the gap is too small: " << prev_x << ", " << t.first << std::endl;); - return false; - } - } - if (t.second.kind() == endpoint_kind::STEND) { - expect_end = false; // swallow a point interval - } - else { - if (prev) - expect_end = !expect_end; - else - expect_end = is_start(t.second); - } - prev = true; - prev_x = t.first; - } -#if Z3DEBUG - // for (int i = 0; i <= 100; i++ ) { - // if ( (m_domain.find(i) != m_domain.end()) != contains(i)) { - // TRACE("disj_intervals", tout << "incorrect value of contains(" << i << ") is = " << contains(i) << std::endl;); - // return false; - // } - // } -#endif - return true; - } -public: - void print(std::ostream & out) const { - if (m_empty) { - out << "empty\n"; - return; - } - if (m_endpoints.empty()){ - out << "[-oo,oo]\n"; - return; - } - bool first = true; - for (auto t : m_endpoints()) { - if (first) { - if (t.second.kind() == endpoint_kind::END) { - out << "[-oo," << t.first; t.second.print(out); out << "]"; - } - else if (t.second.kind() == endpoint_kind::START) { - out << "[" << t.first; t.second.print(out); out << ","; - } else if (t.second.kind() == endpoint_kind::STEND) { - out << "[" << t.first; t.second.print(out); out << "]"; - } - first = false; - } else { - if (t.second.kind() == endpoint_kind::START) { - out << "[" << t.first; t.second.print(out); out << ","; - } - else if (t.second.kind() == endpoint_kind::END) { - out << t.first; t.second.print(out); out << "]"; - } - else if (t.second.kind() == endpoint_kind::STEND) { - out << "[" << t.first; t.second.print(out); out << "]";; - } - } - } - if (has_pos_inf()) - out << "oo]"; - - out << "\n"; - } - - void push() { m_endpoints.push(); m_empty.push(); } - void pop() { m_endpoints.pop(); m_empty.pop(); } - void pop(unsigned k) { while(k--) pop(); } - - bool intersect_with_bound(const T & x, bool is_lower, unsigned explanation) { - return is_lower? intersect_with_lower_bound(x, explanation) : intersect_with_upper_bound(x, explanation); - } - // we intersect the existing set with the half open to the right interval - // returns true if the domain changes - bool intersect_with_lower_bound(const T& x, unsigned explanation) { -#ifdef Z3DEBUG - // for (int i = 0; i < x; i++) - // m_domain.erase(i); -#endif - TRACE("disj_intervals", tout << "intersect_with_lower_bound(" << x << ")\n";); - - if (m_empty) - return false; - if (m_endpoints.empty()) { - set_start(x, explanation); - return true; - } - bool pos_inf = has_pos_inf(); - auto it = m_endpoints.begin(); - while (it != m_endpoints.end() && pos(it) < x) { - m_endpoints.erase(it); - it = m_endpoints.begin(); - } - if (m_endpoints.empty()) { - if (!pos_inf) { - m_empty = true; - return true; - } - set_start(x, explanation); - return true; - } - lp_assert(pos(it) >= x); - if (pos(it) == x) { - if (is_proper_end(it)) { - set_start(x, explanation); - return true; - } - } - else { // x(it) > x - if (is_proper_end(it)) { - set_start(x, explanation); - return true; - } - } - return false; - } -public: - bool intersection_with_upper_bound_is_empty(const T& x) const { - if (has_neg_inf()) - return false; - if (m_empty) - return true; - T b; - lp_assert(get_lower_bound(b)); - get_lower_bound(b); - return x < b; - } - - bool intersection_with_lower_bound_is_empty(const T& x) const { - if (has_pos_inf()) - return false; - if (m_empty) - return true; - T b; - lp_assert(get_upper_bound(b)); - get_upper_bound(b); - return x > b; - } - - // we intersect the existing set with the half open interval - // returns true if there is a change - bool intersect_with_upper_bound(const T& x, unsigned explanation) { -#ifdef Z3DEBUG - // for (int i = 100; i > x; i--) - // m_domain.erase(i); -#endif - TRACE("disj_intervals", tout << "intersect_with_upper_bound(" << x << ")\n";); - if (m_empty) - return false; - if (m_endpoints.empty()) { - set_end(x, explanation); - return true; - } - bool neg_inf = has_neg_inf(); - auto it = m_endpoints.rbegin(); - - bool change = false; - while (!m_endpoints.empty() && pos(it) > x) { - m_endpoints.erase(std::prev(m_endpoints.end())); - it = m_endpoints.rbegin(); - change = true; - } - if (m_endpoints.empty()) { - if (!neg_inf) { - m_empty = true; - return true; - } - set_end(x, explanation); - change = true; - } - lp_assert(pos(it) <= x); - if (pos(it) == x) { - if (is_one_point_interval(it)) {} - else if (is_proper_end(it)) {} - else {// is_proper_start(it->second) - set_end(x, explanation); - change = true; - } - } - else { // pos(it) < x} - if (is_proper_start(it)) { - set_end(x, explanation); - change = true; - } - } - lp_assert(is_correct()); - return change; - } -public: - void intersect_with_interval(const T& x, const T & y) { -#ifdef Z3DEBUG - // for (int i = 0; i <= 100; i++) - // if (i < x || i > y) - // m_domain.erase(i); -#endif - - TRACE("disj_intervals", tout << "intersect_with_interval(" << x << ", " << y <<")\n";); - if (m_empty) - return; - lp_assert(x <= y); - intersect_with_lower_bound(x); - intersect_with_upper_bound(y); - } - - // add an intervar [x, inf] - void unite_with_interval_x_pos_inf(const T& x) { - lp_assert(false); - // if (contains_all()) - // return; - // #if Z3DEBUG - // // for (int i = x; i <= 100; i++) - // // m_domain.insert(i); - // #endif - // TRACE("disj_intervals", tout << "unite_with_interval_x_pos_inf(" << x << ")\n";); - // if (m_empty) { - // set_start(x); - // m_empty = false; - // return; - // } - // bool neg_inf = has_neg_inf(); - // remove_from_the_right(x); - // if (m_endpoints.empty()) { - // if (!neg_inf) - // set_start(x); - // return; - // } - // auto it = m_endpoints.rbegin(); - // lp_assert(pos(it) <= x); - // if (pos(it) == x) { - // if (is_proper_end(it)) { - // m_endpoints.erase(x); - // } else { - // set_start(x); - // } - // } else if (pos(it) == x - 1 && is_end(it)) { - // if (is_proper_start(it)) { - // // do nothing - // } - // else if (is_proper_end(it)) { - // m_endpoints.erase(it); - // } - // else { - // lp_assert(is_one_point_interval(it)); - // set_start(it); - // } - // } else { - // if (!has_pos_inf()) - // set_start(x); - // } - } - - // add an interval [-inf, x] - void unite_with_interval_neg_inf_x(const T& x) { - lp_assert(false); // not implemented - // #if Z3DEBUG - // // for (int i = 0; i <= x; i++) - // // m_domain.insert(i); - // #endif - // TRACE("disj_intervals", tout << "unite_with_interval_neg_inf_x(" << x << ")\n";); - // if (m_empty) { - // set_end(x); - // m_empty = false; - // return; - // } - // bool pos_inf; - // const_iter r; - // bool found_right_point = get_right_point(x, pos_inf, r); - // if (!found_right_point) { - // m_endpoints.clear(); - // set_end(x); - // return; - // } - // if (pos_inf) { - // m_endpoints.clear(); - // return; - // } - // lp_assert(pos(r) >= x); - // if (pos(r) == x || pos(r) == x + 1) { - // if (is_proper_start(r)) - // erase(r); - // else if (is_one_point_interval(r)) { - // set_end(pos(r)); - // } // do nothing for the proper end - // } else { - // if (!is_proper_end(r)) - // set_end(x); - // } - - // while (!m_endpoints.empty() && m_endpoints.begin()->first < x) { - // m_endpoints.erase(m_endpoints.begin()); - // } - // lp_assert(is_correct()); - } - -private: - bool is_start(endpoint_kind x) const { return x == endpoint_kind::START || x == endpoint_kind::STEND; } - bool is_start(const iter & it) const { return is_start(it->second); } - bool is_start(const const_iter & it) const { return is_start(it->second); } - bool is_start(const riter & it) const { return is_start(it->second); } - bool is_start(const endpoint & e) const { return is_start(e.kind()); } - - bool is_end(endpoint_kind x) const { return x == endpoint_kind::END || x == endpoint_kind::STEND; } - bool is_end(const iter & it) const { return is_end(it->second); } - bool is_end(const const_iter & it) const { return is_end(it->second); } - bool is_end(const riter & it) const { return is_end(it->second); } - bool is_end(const endpoint& e) const { return is_end(e.kind()); } - - - T pos(const iter & it) const { - return it->first; - } - T pos(const const_iter & it) const { - return it->first; - } - T pos(const riter & it) const { - return it->first; - } - - T bound_kind(iter & it) const { - return it->second; - } - - T bound_kind(riter & it) const { - return it->second; - } - - bool is_proper_start(endpoint_kind x) const { return x == endpoint_kind::START; } - bool is_proper_start(const riter &x) const { return is_proper_start(x->second);} - bool is_proper_start(const iter &x) const { return is_proper_start(x->second);} - bool is_proper_start(const const_iter &x) const { return is_proper_start(x->second);} - bool is_proper_start(const endpoint &x) const { return is_proper_start(x.kind());} - - bool is_proper_end(endpoint_kind x) const { return x == endpoint_kind::END; } - bool is_proper_end(const iter & it) const { return is_proper_end(it->second.kind()); } - bool is_proper_end(const const_iter & it) const { return is_proper_end(it->second); } - bool is_proper_end(const riter & it) const { return is_proper_end(it->second); } - bool is_proper_end(const endpoint & x) const { return is_proper_end(x.kind()); } - - bool is_one_point_interval(const endpoint & x) const { return is_one_point_interval(x.kind()); } - bool is_one_point_interval(endpoint_kind x) const { return x == endpoint_kind::STEND; } - bool is_one_point_interval(const iter & it) const { return is_one_point_interval(it->second); } - bool is_one_point_interval(const const_iter & it) const { - return is_one_point_interval(it->second); - } - bool is_one_point_interval(const riter & it) const { - return is_one_point_interval(it->second); - } - - void erase(const iter& it) { - m_endpoints.erase(it); - } - void erase(const const_iter& it) { - m_endpoints.erase(it); - } - - void erase(const riter& it) { - m_endpoints.erase(it); - } - - void erase(const T& x) { - m_endpoints.erase(x); - } - - /* void set_one_point_interval(const T& x, unsigned explanation) { - auto it = m_endpoints().find(x); - set_one_point_interval(it, explanation); - }*/ - - void set_start(const_iter &t, unsigned explanation) { - lp_assert(t != m_endpoints.end()); - endpoint e = t->second; - e.m_start_expl = explanation; - m_endpoints[t->first] = e; - } - - void set_start(const T& x, unsigned explanation) { - endpoint e = get_endpoint(x); - e.m_start_expl = explanation; - m_endpoints[x] = e; - } - - endpoint get_endpoint(const T& x) const { - auto it = m_endpoints().find(x); - if (it == m_endpoints().end()) - return endpoint(); - return it->second; - } - - void set_end(const T& x, unsigned expl) { - endpoint e = get_endpoint(x); - e.m_end_expl = expl; - m_endpoints[x] = e; - } - - void set_end(const_iter& t, unsigned expl) { - endpoint e = t->second; - e.m_end_expl = expl; - m_endpoints[t->first] = e; - } - - void set_end(riter& t, unsigned explanation) { - endpoint e = t->second; - e.m_end_expl = expl; - m_endpoints[t->first] = e; - } - -private: - /* void set_start_end(const T& x, const T & y, unsigned expl) { - set_start(x, expl); - set_end(y, expl); - }*/ - - void unite_with_one_point_interval(const T &x) { - TRACE("disj_intervals", tout << "unite_with_one_point_interval(" << x << ")\n";); - if (m_empty) { - m_empty = false; - set_one_point_interval(x); - return; - } - bool has_left_neig, has_right_neig; - bool neg_inf, pos_inf; - const_iter l, r; - has_left_neig = get_left_point(x, neg_inf, l); - has_right_neig = get_right_point(x, pos_inf, r); - if (!has_left_neig) { - if (!has_right_neig) { - set_one_point_interval(x); - return; - } - lp_assert(!neg_inf); - // if (pos(r) == x ) nothing happens - if (pos(r) == x + 1) { - if (is_one_point_interval(r)) { - set_start(x); - set_end(x + 1); - } else { - lp_assert(is_proper_start(r)); - erase(r); - set_start(x); - } - } else { - lp_assert(pos(r) > x); - set_one_point_interval(x); - } - return; - } - - // has_left_neig - if (!has_right_neig) { - if (pos_inf) - return; - // if (pos(l) == x nothing happens - if (pos(l) == x - 1) { - if (is_proper_end(l)) { - erase(l); - set_end(x); - } else { - lp_assert(is_one_point_interval(l)); - set_start(l->first); - set_end(x); - } - } else { - lp_assert(pos(l) < x - 1); - set_one_point_interval(x); - } - return; - } - // has both neighbors - if (neg_inf || pos_inf) - return; - if (pos(l) == x || pos(r) == x) - return; - - // now the cases pos(l) == pos(r) or pos(l) + 1 == pos(r) are impossible - if (is_proper_start(l)) { - lp_assert(is_proper_end(r)); - return; - } - - if (pos(l) + 2 < pos(r)) { // we can glue only to one neighbor - if (pos(l) + 1 == x) { - if (is_proper_end(l)) { - erase(l); - set_end(x); - } - else { - lp_assert(!is_proper_start(l)); - set_start(pos(l)); - set_end(x); - } - } - else if (x + 1 == pos(r)) { - if (is_proper_start(r)) { - erase(r); - } - else if (is_one_point_interval(r)) { - set_end(pos(r)); - } - set_start(x); - } - else { - set_one_point_interval(x); - } - } else { - lp_assert(pos(l) + 2 == pos(r)); - lp_assert(pos(l) + 1 == x); // x is just in between l and r - if (is_proper_end(l)) { - erase(l); - } else { - lp_assert(is_one_point_interval(l)); - set_start(l->first); - } - if (is_proper_start(r)) { - erase(r); - } else { - lp_assert(is_one_point_interval(r)); - set_end(r->first); - } - } - - } - // return false if there are no points y in m_endpoints with pos(y) <= x - // and negative infiniti is not true - bool get_left_point(const T& x, bool &neg_inf, const_iter & l) const { - if (m_empty) { - neg_inf = false; - return false; - } - - if (m_endpoints.empty()) { - neg_inf = true; - return true; - } - - l = m_endpoints.lower_bound(x); - if (l == m_endpoints.end()) { - l--; - neg_inf = false; - return true; - } - if (pos(l) == x) { - neg_inf = false; - return true; - } - if (l == m_endpoints.begin()) - return neg_inf = has_neg_inf(); - l--; - neg_inf = false; - return true; - } - - // return false iff there are no points y in m_endpoints with pos(y) >= x, and positive infinity is not true - bool get_right_point(const T& x, bool &pos_inf, const_iter & r) const { - if (m_empty) { - pos_inf = false; - return false; - } - - if (m_endpoints.empty()) { - pos_inf = true; - return true; - } - - r = m_endpoints.lower_bound(x); - if (r == m_endpoints.end()) { - return pos_inf = has_pos_inf();; - } - - pos_inf = false; - return true; - } - - - - // inserts x, if needed and return a pointer to the leftmost point that is a candidate to removal - /* - iter insert_start_for_unite_with_interval(const T& x) { - auto lbx = m_endpoints.lower_bound(x); - if (lbx == m_endpoints.end()) { - if (!has_pos_inf()) { - T lastx = pos(m_endpoints.rbegin()); - if (lastx + 1 < x) - set_start_end(x, y); - else { - m_endpoints.erase(lastx); - set_end(y); - } - } - return; - } - if (pos(lbx) == x) { - if(is_proper_end(lbx)) { - m_endpoints.erase(x); - } else { set_start(x);} - } - if (pos(lbx) > x) { - if (pos(lbx) > y + 1) { - if (is_end(lbx)) - return; - set_start_end(x, y); - return; - } - if (pos(lpx) == y + 1) { - if (is_end(lbx)) - } - - } - - lp_assert(false); // not implemented - } - */ - - void remove_from_the_left(const T & y ) { - while (!m_endpoints.empty() && pos(m_endpoints.begin()) < y) { - m_endpoints.erase(m_endpoints.begin()); - } - } - - void remove_from_the_left(const T & y, const_iter& l ) { - while (l!= m_endpoints.end() && pos(l) < y) { - l++; - erase(std::prev(l)); - } - } - - - void remove_from_the_right(const T & x) { - while (!m_endpoints.empty() && pos(m_endpoints.rbegin()) > x) { - m_endpoints.erase(m_endpoints.rbegin()); - } - } - - void remove_from_the_right(const T & x, riter & r) { - while (!m_endpoints.empty() && pos(r) > x) { - r--; - m_endpoints.erase(std::next(r)); - } - } -public: - bool get_lower_bound_with_expl(T& b, unsigned & expl) const { - if (m_empty) - return false; - if (has_neg_inf()) - return false; - expl = m_endpoints.begin()->second.m_start_expl; - if (expl == static_cast(-1)) - return false; - b = pos(m_endpoints.begin()); - return true; - } - - bool get_lower_bound(T& b) const { - if (m_empty) - return false; - if (has_neg_inf()) - return false; - b = pos(m_endpoints.begin()); - return true; - } - - int get_lower_bound_expl() const { - if (m_empty) - return -1; - if (has_neg_inf()) - return -1; - return m_endpoints.begin()->second.m_start_expl; - } - - int get_upper_bound_expl() const { - if (m_empty) - return -1; - if (has_pos_inf()) - return -1; - return m_endpoints.rbegin()->second.m_end_expl; - } - - bool get_upper_bound_with_expl(T& b, unsigned & expl) const { - if (m_empty) - return false; - if (has_pos_inf()) - return false; - expl = m_endpoints.rbegin()->second.m_end_expl; - if (expl == static_cast(-1)) - return false; - b = m_endpoints.rbegin()->first; - return true; - } - - bool get_upper_bound_and_kind_with_expl(T& b, endpoint_kind & kind, unsigned & expl) const { - if (m_empty) - return false; - if (has_pos_inf()) - return false; - b = m_endpoints.rbegin()->first; - kind = m_endpoints.rbegin()->second.kind(); - expl = m_endpoints.rbegin()->second.m_explanation; - return true; - } - - bool get_upper_bound(T& b) const { - if (m_empty) - return false; - if (has_pos_inf()) - return false; - b = m_endpoints.rbegin()->first; - return true; - } - - bool is_empty() const { return m_empty; } - - bool is_fixed() const { - if (has_pos_inf() || has_neg_inf()) - return false; - T l; - get_lower_bound(l); - T u; - get_upper_bound(u); - return l==u; - } - - - bool improves_with_lower_bound(const T & v) const { - T b; - bool lower_bound_exists = get_lower_bound(b); - return (!lower_bound_exists || v > b) && - !intersection_with_lower_bound_is_empty(v); - } - - bool improves_with_upper_bound(const T & v) const { - T b; - bool upper_bound_exists = get_upper_bound(b); - return (!upper_bound_exists || v < b) && - !intersection_with_upper_bound_is_empty(v); - } - - // returns true if adding the bound b narrows the domain, but does not make it empty - bool improves(const T & v, bool is_lower_bound) const { - if (is_lower_bound) - return improves_with_lower_bound(v); - return improves_with_upper_bound(v); - } -}; -} diff --git a/src/util/lp/lar_solver.cpp b/src/util/lp/lar_solver.cpp index c8de1807c..d6d556303 100644 --- a/src/util/lp/lar_solver.cpp +++ b/src/util/lp/lar_solver.cpp @@ -331,7 +331,6 @@ void lar_solver::push() { m_term_count.push(); m_constraint_count = m_constraints.size(); m_constraint_count.push(); - m_int_solver->push(); } void lar_solver::clean_popped_elements(unsigned n, int_set& set) { @@ -396,7 +395,6 @@ void lar_solver::pop(unsigned k) { lp_assert(sizes_are_correct()); lp_assert((!m_settings.use_tableau()) || m_mpq_lar_core_solver.m_r_solver.reduced_costs_are_correct_tableau()); m_status = m_mpq_lar_core_solver.m_r_solver.current_x_is_feasible()? lp_status::OPTIMAL: lp_status::UNKNOWN; - m_int_solver->pop(k); } vector lar_solver::get_all_constraint_indices() const { @@ -1476,17 +1474,9 @@ bool lar_solver::strategy_is_undecided() const { return m_settings.simplex_strategy() == simplex_strategy_enum::undecided; } -void lar_solver::catch_up_in_updating_int_solver() { - for (unsigned i = 0; i < constraints().size(); i++) { - m_int_solver->add_constraint_to_chase_cut_solver(i, constraints()[i]); - } -} - var_index lar_solver::add_var(unsigned ext_j, bool is_int) { if (is_int) m_has_int_var = true; - if (is_int && !has_int_var()) - catch_up_in_updating_int_solver(); TRACE("add_var", tout << "adding var " << ext_j << (is_int? " int" : " nonint") << std::endl;); var_index i; @@ -2217,9 +2207,10 @@ void lar_solver::round_to_integer_solution() { } } -bool lar_solver::get_equality_for_term_on_corrent_x(unsigned term_index, mpq & rs) const { +bool lar_solver::get_equality_for_term_on_corrent_x(unsigned term_index, mpq & rs, bool & has_bounds) const { unsigned tj = term_index + m_terms_start_index; auto it = m_ext_vars_to_columns.find(tj); + has_bounds = false; if (it == m_ext_vars_to_columns.end()) return false; unsigned j = it->second.internal_j(); @@ -2227,6 +2218,7 @@ bool lar_solver::get_equality_for_term_on_corrent_x(unsigned term_index, mpq & r impq term_val; bool term_val_ready = false; if (slv.column_has_upper_bound(j)) { + has_bounds = true; const impq & b = slv.m_upper_bounds[j]; lp_assert(is_zero(b.y) && is_int(b.x)); term_val = terms()[term_index]->apply(m_mpq_lar_core_solver.m_r_x); @@ -2237,6 +2229,7 @@ bool lar_solver::get_equality_for_term_on_corrent_x(unsigned term_index, mpq & r } } if (slv.column_has_lower_bound(j)) { + has_bounds = true; if (!term_val_ready) term_val = terms()[term_index]->apply(m_mpq_lar_core_solver.m_r_x); const impq & b = slv.m_lower_bounds[j]; @@ -2247,7 +2240,7 @@ bool lar_solver::get_equality_for_term_on_corrent_x(unsigned term_index, mpq & r return true; } } - return false; + return false; } diff --git a/src/util/lp/lar_solver.h b/src/util/lp/lar_solver.h index 304483817..c7b62ee1d 100644 --- a/src/util/lp/lar_solver.h +++ b/src/util/lp/lar_solver.h @@ -298,9 +298,17 @@ public: static void clean_popped_elements(unsigned n, int_set& set); static void shrink_inf_set_after_pop(unsigned n, int_set & set); - void pop(unsigned k); + + class scoped_push { + lar_solver& m_solver; + bool m_pop; + public: + scoped_push(lar_solver& s):m_solver(s), m_pop(true) { s.push(); } + ~scoped_push() { if (m_pop) m_solver.pop(); } + void pop() { SASSERT(m_pop); m_solver.pop(); m_pop = false; } + }; vector get_all_constraint_indices() const; @@ -579,6 +587,6 @@ public: unsigned column_count() const { return A_r().column_count(); } const vector & r_basis() const { return m_mpq_lar_core_solver.r_basis(); } const vector & r_nbasis() const { return m_mpq_lar_core_solver.r_nbasis(); } - bool get_equality_for_term_on_corrent_x(unsigned i, mpq &rs) const; + bool get_equality_and_right_side_for_term_on_corrent_x(unsigned i, mpq &rs, bool & has_bounds) const; }; } diff --git a/src/util/lp/lp_settings.h b/src/util/lp/lp_settings.h index eb460d1ad..d8c877e8e 100644 --- a/src/util/lp/lp_settings.h +++ b/src/util/lp/lp_settings.h @@ -133,54 +133,77 @@ private: }; default_lp_resource_limit m_default_resource_limit; - lp_resource_limit* m_resource_limit; + lp_resource_limit* m_resource_limit; // used for debug output - std::ostream* m_debug_out; + std::ostream* m_debug_out; // used for messages, for example, the computation progress messages - std::ostream* m_message_out; + std::ostream* m_message_out; - stats m_stats; - random_gen m_rand; + stats m_stats; + random_gen m_rand; public: - unsigned reps_in_scaler; + unsigned reps_in_scaler; // when the absolute value of an element is less than pivot_epsilon // in pivoting, we treat it as a zero - double pivot_epsilon; + double pivot_epsilon; // see Chatal, page 115 - double positive_price_epsilon; + double positive_price_epsilon; // a quatation "if some choice of the entering vairable leads to an eta matrix // whose diagonal element in the eta column is less than e2 (entering_diag_epsilon) in magnitude, the this choice is rejected ... - double entering_diag_epsilon; - int c_partial_pivoting; // this is the constant c from page 410 - unsigned depth_of_rook_search; - bool using_partial_pivoting; + double entering_diag_epsilon; + int c_partial_pivoting; // this is the constant c from page 410 + unsigned depth_of_rook_search; + bool using_partial_pivoting; // dissertation of Achim Koberstein // if Bx - b is different at any component more that refactor_epsilon then we refactor - double refactor_tolerance; - double pivot_tolerance; - double zero_tolerance; - double drop_tolerance; - double tolerance_for_artificials; - double can_be_taken_to_basis_tolerance; + double refactor_tolerance; + double pivot_tolerance; + double zero_tolerance; + double drop_tolerance; + double tolerance_for_artificials; + double can_be_taken_to_basis_tolerance; - unsigned percent_of_entering_to_check; // we try to find a profitable column in a percentage of the columns - bool use_scaling; - double scaling_maximum; - double scaling_minimum; - double harris_feasibility_tolerance; // page 179 of Istvan Maros - double ignore_epsilon_of_harris; - unsigned max_number_of_iterations_with_no_improvements; - unsigned max_total_number_of_iterations; - double time_limit; // the maximum time limit of the total run time in seconds + unsigned percent_of_entering_to_check; // we try to find a profitable column in a percentage of the columns + bool use_scaling; + double scaling_maximum; + double scaling_minimum; + double harris_feasibility_tolerance; // page 179 of Istvan Maros + double ignore_epsilon_of_harris; + unsigned max_number_of_iterations_with_no_improvements; + unsigned max_total_number_of_iterations; + double time_limit; // the maximum time limit of the total run time in seconds // dual section - double dual_feasibility_tolerance; // // page 71 of the PhD thesis of Achim Koberstein - double primal_feasibility_tolerance; // page 71 of the PhD thesis of Achim Koberstein - double relative_primal_feasibility_tolerance; // page 71 of the PhD thesis of Achim Koberstein - bool hnf_cutter_exit_if_x_is_not_on_bound_or_mixed = true; - - bool m_bound_propagation; + double dual_feasibility_tolerance; // // page 71 of the PhD thesis of Achim Koberstein + double primal_feasibility_tolerance; // page 71 of the PhD thesis of Achim Koberstein + double relative_primal_feasibility_tolerance; // page 71 of the PhD thesis of Achim Koberstein + // end of dual section + bool m_bound_propagation; + bool presolve_with_double_solver_for_lar; + simplex_strategy_enum m_simplex_strategy; + int report_frequency; + bool print_statistics; + unsigned column_norms_update_frequency; + bool scale_with_ratio; + double density_threshold; + bool use_breakpoints_in_feasibility_search; + unsigned max_row_length_for_bound_propagation; + bool backup_costs; + unsigned column_number_threshold_for_using_lu_in_lar_solver; + unsigned m_int_gomory_cut_period; + unsigned m_int_chase_cut_solver_period; + unsigned m_int_find_cube_period; + unsigned m_hnf_cut_period; + bool m_int_run_gcd_test; + unsigned m_chase_cut_solver_cycle_on_var; + bool m_int_pivot_fixed_vars_from_basis; + bool m_int_patch_only_integer_values; + unsigned limit_on_rows_for_hnf_cutter; + + unsigned random_next() { return m_rand(); } + void set_random_seed(unsigned s) { m_rand.set_seed(s); } + bool bound_progation() const { return m_bound_propagation; } @@ -237,12 +260,12 @@ public: m_int_gomory_cut_period(4), m_int_chase_cut_solver_period(8), m_int_find_cube_period(4), - m_int_cuts_etc_period(4), m_hnf_cut_period(4), m_int_run_gcd_test(true), m_chase_cut_solver_cycle_on_var(10), m_int_pivot_fixed_vars_from_basis(false), - m_int_patch_only_integer_values(true) + m_int_patch_only_integer_values(true), + limit_on_rows_for_hnf_cutter(100) {} void set_resource_limit(lp_resource_limit& lim) { m_resource_limit = &lim; } @@ -312,8 +335,6 @@ public: return is_eps_small_general(t, tolerance_for_artificials); } // the method of lar solver to use - bool presolve_with_double_solver_for_lar; - simplex_strategy_enum m_simplex_strategy; simplex_strategy_enum simplex_strategy() const { return m_simplex_strategy; } @@ -335,29 +356,9 @@ public: return m_simplex_strategy == simplex_strategy_enum::tableau_rows; } - int report_frequency; - bool print_statistics; - unsigned column_norms_update_frequency; - bool scale_with_ratio; - double density_threshold; // need to tune it up, todo #ifdef Z3DEBUG static unsigned ddd; // used for debugging #endif - bool use_breakpoints_in_feasibility_search; - unsigned random_next() { return m_rand(); } - void set_random_seed(unsigned s) { m_rand.set_seed(s); } - unsigned max_row_length_for_bound_propagation; - bool backup_costs; - unsigned column_number_threshold_for_using_lu_in_lar_solver; - unsigned m_int_gomory_cut_period; - unsigned m_int_chase_cut_solver_period; - unsigned m_int_find_cube_period; - unsigned m_int_cuts_etc_period; - unsigned m_hnf_cut_period; - bool m_int_run_gcd_test; - unsigned m_chase_cut_solver_cycle_on_var; - bool m_int_pivot_fixed_vars_from_basis; - bool m_int_patch_only_integer_values; }; // end of lp_settings class diff --git a/src/util/lp/stacked_map.h b/src/util/lp/stacked_map.h index caa80636c..e1a8963fb 100644 --- a/src/util/lp/stacked_map.h +++ b/src/util/lp/stacked_map.h @@ -24,6 +24,7 @@ Revision History: #include #include #include +#include namespace lp { template