From 7e82ab595ec12beccc11660d80333f0c07e56983 Mon Sep 17 00:00:00 2001 From: Lev Nachmanson Date: Thu, 22 Mar 2018 20:35:13 -0700 Subject: [PATCH] follow the smalles branch Signed-off-by: Lev Nachmanson correction in the sign of gomory_cut Signed-off-by: Lev Nachmanson fix in the gomory cut sign Signed-off-by: Lev Nachmanson try using lemmas of cut_solver as cuts Signed-off-by: Lev Nachmanson add find_cube() proposed by Nikolaj Signed-off-by: Lev Nachmanson restore m_int_branch_cut_solver to 8 Signed-off-by: Lev Nachmanson accept empty lar_terms in theory_lra and also do not create empty lar_terms/lemmas Signed-off-by: Lev Nachmanson qflia_tactic Signed-off-by: Lev Nachmanson call find_feasible solution to recover for a failure in find_cube Signed-off-by: Lev Nachmanson do not tighten unused terms Signed-off-by: Lev Nachmanson get rid of inf_int_set Signed-off-by: Lev Nachmanson fix a bug with an accidental solution in cube Signed-off-by: Lev Nachmanson get rid of inf_int_set Signed-off-by: Lev Nachmanson bug fix with has_int_var() for lar_solver Signed-off-by: Lev Nachmanson fix in find_inf_int_base_column Signed-off-by: Lev Nachmanson --- src/smt/theory_lra.cpp | 30 +++-- src/util/lp/CMakeLists.txt | 1 - src/util/lp/cut_solver.h | 40 ++++++- src/util/lp/int_solver.cpp | 165 +++++++++++++++++--------- src/util/lp/int_solver.h | 15 +-- src/util/lp/lar_core_solver.h | 1 + src/util/lp/lar_solver.cpp | 131 +++++++++++++------- src/util/lp/lar_solver.h | 63 ++++------ src/util/lp/lp_core_solver_base.h | 13 -- src/util/lp/lp_core_solver_base_def.h | 3 +- src/util/lp/lp_settings.h | 17 +-- src/util/lp/polynomial.h | 9 ++ src/util/lp/stacked_vector.h | 8 ++ 13 files changed, 314 insertions(+), 182 deletions(-) diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp index c240d243a..54bf366b8 100644 --- a/src/smt/theory_lra.cpp +++ b/src/smt/theory_lra.cpp @@ -1095,7 +1095,7 @@ public: } void init_variable_values() { - if (m_solver.get() && th.get_num_vars() > 0) { + if (!m.canceled() && m_solver.get() && th.get_num_vars() > 0) { m_solver->get_model(m_variable_values); } } @@ -1255,7 +1255,7 @@ public: return FC_GIVEUP; } - // create a bound atom representing term <= k + // create a bound atom representing term >= k is lower_bound is true, and term <= k if it is false app_ref mk_bound(lp::lar_term const& term, rational const& k, bool lower_bound) { app_ref t = mk_term(term, k.is_int()); app_ref atom(m); @@ -1289,7 +1289,7 @@ public: case lp::lia_move::ok: return l_true; case lp::lia_move::branch: { - app_ref b = mk_bound(term, k, upper); + app_ref b = mk_bound(term, k, !upper); // branch on term >= k + 1 // branch on term <= k // TBD: ctx().force_phase(ctx().get_literal(b)); @@ -1300,7 +1300,7 @@ public: case lp::lia_move::cut: { ++m_stats.m_gomory_cuts; // m_explanation implies term <= k - app_ref b = mk_bound(term, k, upper); + app_ref b = mk_bound(term, k, false); m_eqs.reset(); m_core.reset(); m_params.reset(); @@ -2683,10 +2683,14 @@ public: if (!term.m_v.is_zero()) { args.push_back(a.mk_numeral(term.m_v, is_int)); } - if (args.size() == 1) { + switch (args.size()) { + case 0: + return app_ref(a.mk_numeral(rational::zero(), is_int), m); + case 1: return app_ref(to_app(args[0].get()), m); + default: + return app_ref(a.mk_add(args.size(), args.c_ptr()), m); } - return app_ref(a.mk_add(args.size(), args.c_ptr()), m); } app_ref mk_obj(theory_var v) { @@ -2806,12 +2810,14 @@ public: st.update("arith-make-feasible", m_solver->settings().st().m_make_feasible); st.update("arith-max-columns", m_solver->settings().st().m_max_cols); st.update("arith-max-rows", m_solver->settings().st().m_max_rows); - st.update("cut_solver-calls", m_solver->settings().st().m_cut_solver_calls); - st.update("cut_solver-true", m_solver->settings().st().m_cut_solver_true); - st.update("cut_solver-false", m_solver->settings().st().m_cut_solver_false); - st.update("cut_solver-undef", m_solver->settings().st().m_cut_solver_undef); - st.update("gcd_calls", m_solver->settings().st().m_gcd_calls); - st.update("gcd_conflict", m_solver->settings().st().m_gcd_conflicts); + st.update("cut-solver-calls", m_solver->settings().st().m_cut_solver_calls); + st.update("cut-solver-true", m_solver->settings().st().m_cut_solver_true); + st.update("cut-solver-false", m_solver->settings().st().m_cut_solver_false); + st.update("cut-solver-undef", m_solver->settings().st().m_cut_solver_undef); + st.update("gcd-calls", m_solver->settings().st().m_gcd_calls); + st.update("gcd-conflict", m_solver->settings().st().m_gcd_conflicts); + st.update("cube-calls", m_solver->settings().st().m_cube_calls); + st.update("cube-success", m_solver->settings().st().m_cube_success); } }; diff --git a/src/util/lp/CMakeLists.txt b/src/util/lp/CMakeLists.txt index 0a821f9b7..2acf08ca8 100644 --- a/src/util/lp/CMakeLists.txt +++ b/src/util/lp/CMakeLists.txt @@ -23,7 +23,6 @@ z3_add_component(lp matrix.cpp nra_solver.cpp permutation_matrix.cpp - quick_xplain.cpp row_eta_matrix.cpp scaler.cpp sparse_matrix.cpp diff --git a/src/util/lp/cut_solver.h b/src/util/lp/cut_solver.h index e39244006..d7c0fdcad 100644 --- a/src/util/lp/cut_solver.h +++ b/src/util/lp/cut_solver.h @@ -33,7 +33,7 @@ #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 cut_solver; //forward definition @@ -397,6 +397,7 @@ public: unsigned m_pushes_to_trail; indexer_of_constraints m_indexer_of_constraints; + bool is_lower_bound(literal & l) const { return l.is_lower(); } @@ -2132,11 +2133,11 @@ public: } if (number_of_lemmas == m_lemmas_container.size()) { - if (lemma_has_been_modified) { + if (lemma_has_been_modified && lemma->poly().number_of_monomials() != 0) { add_lemma_as_not_active(lemma); } else { - delete_constraint(lemma);; + delete_constraint(lemma); } } } @@ -2420,6 +2421,7 @@ public: } 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); @@ -2485,7 +2487,39 @@ public: } 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("cut_solver_cuts", trace_print_constraint(tout, *c);); + return true; + } }; inline polynomial operator*(const mpq & a, polynomial & p) { diff --git a/src/util/lp/int_solver.cpp b/src/util/lp/int_solver.cpp index 1b948e2d6..03ab1424c 100644 --- a/src/util/lp/int_solver.cpp +++ b/src/util/lp/int_solver.cpp @@ -44,29 +44,35 @@ void int_solver::trace_inf_rows() const { ); } -int_set& int_solver::inf_int_set() { - return m_lar_solver->m_inf_int_set; -} - -const int_set& int_solver::inf_int_set() const { - return m_lar_solver->m_inf_int_set; -} - bool int_solver::has_inf_int() const { - return !inf_int_set().is_empty(); + return m_lar_solver->has_inf_int(); } int int_solver::find_inf_int_base_column() { - if (inf_int_set().is_empty()) + unsigned inf_int_count; + int j = find_inf_int_boxed_base_column_with_smallest_range(inf_int_count); + if (j != -1) + return j; + if (inf_int_count == 0) return -1; - int j = find_inf_int_boxed_base_column_with_smallest_range(); - if (j != -1) - return j; - unsigned k = settings().random_next() % inf_int_set().m_index.size(); - return inf_int_set().m_index[k]; + unsigned k = random() % inf_int_count; + return get_kth_inf_int(k); } -int int_solver::find_inf_int_boxed_base_column_with_smallest_range() { +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) + return j; + } + lp_assert(false); + return -1; +} + +int int_solver::find_inf_int_boxed_base_column_with_smallest_range(unsigned & inf_int_count) { + inf_int_count = 0; int result = -1; mpq range; mpq new_range; @@ -74,11 +80,13 @@ int int_solver::find_inf_int_boxed_base_column_with_smallest_range() { unsigned n; lar_core_solver & lcs = m_lar_solver->m_mpq_lar_core_solver; - for (int j : inf_int_set().m_index) { - lp_assert(is_base(j) && column_is_int_inf(j)); - lp_assert(!is_fixed(j)); + for (unsigned j : m_lar_solver->r_basis()) { + if (!column_is_int_inf(j)) + continue; + inf_int_count++; if (!is_boxed(j)) continue; + lp_assert(!is_fixed(j)); new_range = lcs.m_r_upper_bounds()[j].x - lcs.m_r_lower_bounds()[j].x; if (new_range > small_range_thresold) continue; @@ -325,7 +333,7 @@ lia_move int_solver::mk_gomory_cut(lar_term& t, mpq& k, explanation & expl, unsi tout << "inf_col = " << inf_col << std::endl; ); - // gomory will be t >= k + // gomory will be t <= k and the current solution has a property t > k k = 1; mpq lcm_den(1); unsigned x_j; @@ -356,7 +364,7 @@ lia_move int_solver::mk_gomory_cut(lar_term& t, mpq& k, explanation & expl, unsi lp_assert(current_solution_is_inf_on_cut(t, k)); m_lar_solver->subs_term_columns(t); - TRACE("gomory_cut", tout<<"precut:"; m_lar_solver->print_term(t, tout); tout << ">= " << k << std::endl;); + TRACE("gomory_cut", tout<<"precut:"; m_lar_solver->print_term(t, tout); tout << " <= " << k << std::endl;); return lia_move::cut; } @@ -476,15 +484,69 @@ void int_solver::copy_values_from_cut_solver() { } void int_solver::catch_up_in_adding_constraints_to_cut_solver() { - lp_assert(m_cut_solver.number_of_asserts() <= m_lar_solver->constraints().size()); + lp_assert(m_cut_solver.number_of_asserts() <= m_lar_solver->constraints().size()); for (unsigned j = m_cut_solver.number_of_asserts(); j < m_lar_solver->constraints().size(); j++) { add_constraint_to_cut_solver(j, m_lar_solver->constraints()[j]); } } +bool int_solver::tighten_term_for_cube(unsigned i) { + unsigned ti = i + m_lar_solver->terms_start_index(); + if (!m_lar_solver->term_is_used_as_row(ti)) + return true; + const lar_term* t = m_lar_solver->terms()[i]; + mpq delta = zero_of_type(); + for (const auto & p : *t) { + delta += abs(p.coeff()); + } + delta *= mpq(1, 2); + TRACE("cube", m_lar_solver->print_term_as_indices(*t, tout); tout << ", delta = " << delta;); + + return m_lar_solver->tighten_term_bounds_by_delta(i, delta); +} + +bool int_solver::tighten_terms_for_cube() { + for (unsigned i = 0; i < m_lar_solver->terms().size(); i++) + if (!tighten_term_for_cube(i)) { + TRACE("cube", tout << "cannot tighten";); + return false; + } + return true; +} + +bool int_solver::find_cube() { + if (m_branch_cut_counter % settings().m_int_branch_find_cube != 0) + return false; + + settings().st().m_cube_calls++; + TRACE("cube", + for (unsigned j = 0; j < m_lar_solver->A_r().column_count(); j++) + 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; + } + + 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(); + move_non_basic_columns_to_bounds(); + m_lar_solver->find_feasible_solution(); + lp_assert(m_lar_solver->get_status() == lp_status::OPTIMAL); + // it can happen that we found an integer solution here + return !m_lar_solver->r_basis_has_inf_int(); + } + m_lar_solver->round_to_integer_solution(); + m_lar_solver->pop(); + return true; +} + lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex, bool & upper) { init_check_data(); - lp_assert(inf_int_set_is_correct()); // it is a reimplementation of // final_check_status theory_arith::check_int_feasibility() // from theory_arith_int.h with the addition of cut_solver @@ -507,7 +569,14 @@ lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex, bool & upper) { if (!has_inf_int()) return lia_move::ok; - if ((++m_branch_cut_counter) % settings().m_int_branch_cut_solver == 0) { + ++m_branch_cut_counter; + if (find_cube()){ + settings().st().m_cube_success++; + return lia_move::ok; + } + TRACE("cube", tout << "cube did not succeed";); + + if ((m_branch_cut_counter) % settings().m_int_branch_cut_solver == 0) { TRACE("check_main_int", tout<<"cut_solver";); catch_up_in_adding_constraints_to_cut_solver(); auto check_res = m_cut_solver.check(); @@ -520,9 +589,18 @@ lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex, bool & upper) { case cut_solver::lbool::l_true: settings().st().m_cut_solver_true++; copy_values_from_cut_solver(); + lp_assert(m_lar_solver->all_constraints_hold()); return lia_move::ok; case cut_solver::lbool::l_undef: settings().st().m_cut_solver_undef++; + if (m_cut_solver.try_getting_cut(t, k, m_lar_solver->m_mpq_lar_core_solver.m_r_x)) { + m_lar_solver->subs_term_columns(t); + TRACE("cut_solver_cuts", + tout<<"precut from cut_solver:"; m_lar_solver->print_term(t, tout); tout << " <= " << k << std::endl;); + + + return lia_move::cut; + } break; default: return lia_move::give_up; @@ -598,7 +676,6 @@ void int_solver::set_value_for_nbasic_column_ignore_old_values(unsigned j, const auto & x = m_lar_solver->m_mpq_lar_core_solver.m_r_x[j]; auto delta = new_val - x; x = new_val; - update_column_in_int_inf_set(j); m_lar_solver->change_basic_columns_dependend_on_a_given_nb_column(j, delta); } @@ -612,7 +689,6 @@ void int_solver::set_value_for_nbasic_column(unsigned j, const impq & new_val) { } auto delta = new_val - x; x = new_val; - update_column_in_int_inf_set(j); m_lar_solver->change_basic_columns_dependend_on_a_given_nb_column(j, delta); } @@ -677,7 +753,6 @@ void int_solver::patch_int_infeasible_nbasic_columns() { move_non_basic_columns_to_bounds(); m_lar_solver->find_feasible_solution(); } - lp_assert(is_feasible() && inf_int_set_is_correct()); } mpq get_denominators_lcm(const row_strip & row) { @@ -990,26 +1065,9 @@ void int_solver::display_column(std::ostream & out, unsigned j) const { m_lar_solver->m_mpq_lar_core_solver.m_r_solver.print_column_info(j, out); } -bool int_solver::inf_int_set_is_correct() const { - for (unsigned j = 0; j < m_lar_solver->A_r().column_count(); j++) { - if (inf_int_set().contains(j) != (is_int(j) && (!value_is_int(j)))) { - TRACE("arith_int", tout << "j= " << j << " inf_int_set().contains(j) = " << inf_int_set().contains(j) << ", is_int(j) = " << is_int(j) << "\nvalue_is_int(j) = " << value_is_int(j) << ", val = " << get_value(j) << std::endl;); - return false; - } - } - return true; -} - bool int_solver::column_is_int_inf(unsigned j) const { return is_int(j) && (!value_is_int(j)); } - -void int_solver::update_column_in_int_inf_set(unsigned j) { - if (is_int(j) && (!value_is_int(j))) - inf_int_set().insert(j); - else - inf_int_set().erase(j); -} bool int_solver::is_base(unsigned j) const { return m_lar_solver->m_mpq_lar_core_solver.m_r_heading[j] >= 0; @@ -1181,6 +1239,7 @@ const impq& int_solver::lower_bound(unsigned j) const { } lia_move int_solver::create_branch_on_column(int j, lar_term& t, mpq& k, bool free_column, bool & upper) { + TRACE("check_main_int", tout << "branching" << std::endl;); lp_assert(t.is_empty()); lp_assert(j != -1); t.add_monomial(mpq(1), m_lar_solver->adjust_column_index_to_term_index(j)); @@ -1201,7 +1260,6 @@ lia_move int_solver::create_branch_on_column(int j, lar_term& t, mpq& k, bool fr } bool int_solver::left_branch_is_more_narrow_than_right(unsigned j) { - return settings().random_next() % 2; switch (m_lar_solver->m_mpq_lar_core_solver.m_r_solver.m_column_types[j] ) { case column_type::fixed: return false; @@ -1222,16 +1280,7 @@ bool int_solver::left_branch_is_more_narrow_than_right(unsigned j) { const impq& int_solver::upper_bound(unsigned j) const { return m_lar_solver->column_upper_bound(j); } -void int_solver::display_inf_or_int_inf_columns(std::ostream & out) const { - out << "int inf\n"; - for (unsigned j : m_lar_solver->m_inf_int_set.m_index) { - display_column(out, j); - } - out << "regular inf\n"; - for (unsigned j : m_lar_solver->m_mpq_lar_core_solver.m_r_solver.m_inf_set.m_index) { - display_column(out, j); - } -} + bool int_solver::is_term(unsigned j) const { return m_lar_solver->column_corresponds_to_term(j); } @@ -1250,8 +1299,8 @@ void int_solver::pop(unsigned k) { m_cut_solver.pop_constraints(); } -void int_solver::push() { - m_cut_solver.push(); -} +void int_solver::push() { m_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 badb55beb..f422273f1 100644 --- a/src/util/lp/int_solver.h +++ b/src/util/lp/int_solver.h @@ -59,8 +59,7 @@ public: cut_solver m_cut_solver; // methods int_solver(lar_solver* lp); - int_set& inf_int_set(); - const int_set& inf_int_set() const; + // main function to check that the solution provided by lar_solver is valid for integral values, // or provide a way of how it can be adjusted. lia_move check(lar_term& t, mpq& k, explanation& ex, bool & upper); @@ -115,13 +114,11 @@ private: void failed(); bool is_feasible() const; const impq & get_value(unsigned j) const; - void display_column(std::ostream & out, unsigned j) const; - bool inf_int_set_is_correct() const; - void update_column_in_int_inf_set(unsigned j); bool column_is_int_inf(unsigned j) const; void trace_inf_rows() const; int find_inf_int_base_column(); - int find_inf_int_boxed_base_column_with_smallest_range(); + int find_inf_int_boxed_base_column_with_smallest_range(unsigned&); + int get_kth_inf_int(unsigned) const; lp_settings& settings(); bool move_non_basic_columns_to_bounds(); void branch_infeasible_int_var(unsigned); @@ -143,6 +140,7 @@ private: } public: + void display_column(std::ostream & out, unsigned j) const; inline static mpq fractional_part(const impq & n) { lp_assert(is_rational(n)); @@ -164,7 +162,6 @@ private: lia_move create_branch_on_column(int j, lar_term& t, mpq& k, bool free_column, bool & upper); void catch_up_in_adding_constraints_to_cut_solver(); public: - void display_inf_or_int_inf_columns(std::ostream & out) const; template void fill_cut_solver_vars(); template @@ -176,5 +173,9 @@ public: void push(); void copy_values_from_cut_solver(); bool left_branch_is_more_narrow_than_right(unsigned); + bool find_cube(); + bool tighten_terms_for_cube(); + bool tighten_term_for_cube(unsigned); + unsigned column_count() const; }; } diff --git a/src/util/lp/lar_core_solver.h b/src/util/lp/lar_core_solver.h index 70277c848..0ddf03113 100644 --- a/src/util/lp/lar_core_solver.h +++ b/src/util/lp/lar_core_solver.h @@ -830,5 +830,6 @@ public: } } + const vector& r_basis() const { return m_r_basis; } }; } diff --git a/src/util/lp/lar_solver.cpp b/src/util/lp/lar_solver.cpp index f0c9cf09a..0b85f85b0 100644 --- a/src/util/lp/lar_solver.cpp +++ b/src/util/lp/lar_solver.cpp @@ -31,11 +31,8 @@ lar_solver::lar_solver() : m_status(lp_status::OPTIMAL), m_infeasible_column_index(-1), m_terms_start_index(1000000), m_mpq_lar_core_solver(m_settings, *this), - m_tracker_of_x_change([&](unsigned j) { - call_assignment_tracker(j); - } - ), - m_int_solver(nullptr) + m_int_solver(nullptr), + m_has_int_var(false) {} void lar_solver::set_track_pivoted_rows(bool v) { @@ -277,11 +274,10 @@ lp_status lar_solver::get_status() const { return m_status;} void lar_solver::set_status(lp_status s) {m_status = s;} bool lar_solver::has_int_var() const { - return m_mpq_lar_core_solver.m_r_solver.m_tracker_of_x_change != nullptr; + return m_has_int_var; } lp_status lar_solver::find_feasible_solution() { - lp_assert(inf_int_set_is_correct()); m_settings.st().m_make_feasible++; if (A_r().column_count() > m_settings.st().m_max_cols) m_settings.st().m_max_cols = A_r().column_count(); @@ -290,19 +286,12 @@ lp_status lar_solver::find_feasible_solution() { if (strategy_is_undecided()) decide_on_strategy_and_adjust_initial_state(); - if (has_int_var()) { - m_inf_int_set.resize(A_r().column_count()); - } - m_mpq_lar_core_solver.m_r_solver.m_look_for_feasible_solution_only = true; auto ret = solve(); - TRACE("intinf", m_int_solver->display_inf_or_int_inf_columns(tout);); - lp_assert(inf_int_set_is_correct()); return ret; } lp_status lar_solver::solve() { - lp_assert(inf_int_set_is_correct()); if (m_status == lp_status::INFEASIBLE) { return m_status; } @@ -313,7 +302,6 @@ lp_status lar_solver::solve() { } m_columns_with_changed_bound.clear(); - lp_assert(inf_int_set_is_correct()); return m_status; } @@ -363,7 +351,6 @@ void lar_solver::shrink_inf_set_after_pop(unsigned n, int_set & set) { void lar_solver::pop(unsigned k) { TRACE("arith_int", tout << "pop" << std::endl;); - lp_assert(inf_int_set_is_correct()); TRACE("lar_solver", tout << "k = " << k << std::endl;); int n_was = static_cast(m_ext_vars_to_columns.size()); @@ -381,11 +368,8 @@ void lar_solver::pop(unsigned k) { m_mpq_lar_core_solver.pop(k); clean_popped_elements(n, m_columns_with_changed_bound); - clean_popped_elements(n, m_inf_int_set); unsigned m = A_r().row_count(); - lp_assert(inf_int_set_is_correct()); clean_popped_elements(m, m_rows_with_changed_bounds); - lp_assert(inf_int_set_is_correct()); clean_inf_set_of_r_solver_after_pop(); m_status = m_mpq_lar_core_solver.m_r_solver.current_x_is_feasible()? lp_status::OPTIMAL: lp_status::UNKNOWN; lp_assert(m_settings.simplex_strategy() == simplex_strategy_enum::undecided || @@ -663,7 +647,6 @@ bool lar_solver::costs_are_used() const { } void lar_solver::change_basic_columns_dependend_on_a_given_nb_column(unsigned j, const numeric_pair & delta) { - lp_assert(inf_int_set_is_correct()); if (use_tableau()) { for (const auto & c : A_r().m_columns[j]) { unsigned bj = m_mpq_lar_core_solver.m_r_basis[c.m_i]; @@ -685,7 +668,6 @@ void lar_solver::change_basic_columns_dependend_on_a_given_nb_column(unsigned j, m_mpq_lar_core_solver.m_r_solver.update_x_with_delta_and_track_feasibility(bj, -m_column_buffer[i] * delta); } } - lp_assert(inf_int_set_is_correct()); } void lar_solver::update_x_and_inf_costs_for_column_with_changed_bounds(unsigned j) { @@ -755,10 +737,8 @@ void lar_solver::solve_with_core_solver() { update_x_and_inf_costs_for_columns_with_changed_bounds_tableau(); else update_x_and_inf_costs_for_columns_with_changed_bounds(); - TRACE("intinf", m_int_solver->display_inf_or_int_inf_columns(tout);); m_mpq_lar_core_solver.solve(); set_status(m_mpq_lar_core_solver.m_r_solver.get_status()); - lp_assert(inf_int_set_is_correct()); lp_assert(m_status != lp_status::OPTIMAL || all_constraints_hold()); } @@ -1140,8 +1120,8 @@ void lar_solver::get_infeasibility_explanation_for_inf_sign( } void lar_solver::get_model(std::unordered_map & variable_values) const { - mpq delta = mpq(1, 2); // start from 0.5 to have less clashes lp_assert(m_status == lp_status::OPTIMAL); + mpq delta = mpq(1, 2); // start from 0.5 to have less clashes unsigned i; do { // different pairs have to produce different singleton values @@ -1276,7 +1256,6 @@ void lar_solver::random_update(unsigned sz, var_index const * vars) { fill_var_set_for_random_update(sz, vars, column_list); random_updater ru(*this, column_list); ru.update(); - lp_assert(inf_int_set_is_correct()); } @@ -1406,7 +1385,6 @@ void lar_solver::pop_tableau() { } void lar_solver::clean_inf_set_of_r_solver_after_pop() { - lp_assert(inf_int_set_is_correct()); vector became_feas; clean_popped_elements(A_r().column_count(), m_mpq_lar_core_solver.m_r_solver.m_inf_set); std::unordered_set basic_columns_with_changed_cost; @@ -1449,12 +1427,6 @@ void lar_solver::clean_inf_set_of_r_solver_after_pop() { } } -void lar_solver::shrink_explanation_to_minimum(vector> & explanation) const { - // implementing quickXplain - quick_xplain::run(explanation, *this); - lp_assert(this->explanation_is_correct(explanation)); -} - bool lar_solver::model_is_int_feasible() const { unsigned n = A_r().column_count(); for (unsigned j = 0; j < n; j++) { @@ -1495,7 +1467,7 @@ bool lar_solver::column_is_fixed(unsigned j) const { bool lar_solver::ext_var_is_int(var_index ext_var) const { auto it = m_ext_vars_to_columns.find(ext_var); lp_assert(it != m_ext_vars_to_columns.end()); - return it == m_ext_vars_to_columns.end() || it->second.is_integer(); + return it->second.is_integer(); } // below is the initialization functionality of lar_solver @@ -1511,6 +1483,8 @@ void lar_solver::catch_up_in_updating_int_solver() { } 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(); @@ -1529,10 +1503,6 @@ var_index lar_solver::add_var(unsigned ext_j, bool is_int) { m_columns_to_ul_pairs.push_back(ul_pair(static_cast(-1))); add_non_basic_var_to_core_fields(ext_j, is_int); lp_assert(sizes_are_correct()); - if (is_int) { - m_mpq_lar_core_solver.m_r_solver.set_tracker_of_x(& m_tracker_of_x_change); - } - lp_assert(inf_int_set_is_correct()); return i; } @@ -1548,7 +1518,6 @@ void lar_solver::add_non_basic_var_to_core_fields(unsigned ext_j, bool is_int) { register_new_ext_var_index(ext_j, is_int); m_mpq_lar_core_solver.m_column_types.push_back(column_type::free_column); m_columns_with_changed_bound.increase_size_by_one(); - m_inf_int_set.increase_size_by_one(); add_new_var_to_core_fields_for_mpq(false); if (use_lu()) add_new_var_to_core_fields_for_doubles(false); @@ -1686,7 +1655,6 @@ void lar_solver::add_row_from_term_no_constraint(const lar_term * term, unsigned m_mpq_lar_core_solver.m_r_solver.update_x_and_call_tracker(j, get_basic_var_value_from_row_directly(A_r().row_count() - 1)); if (use_lu()) fill_last_row_of_A_d(A_d(), term); - lp_assert(inf_int_set_is_correct()); } void lar_solver::add_basic_var_to_core_fields() { @@ -1695,7 +1663,6 @@ void lar_solver::add_basic_var_to_core_fields() { m_mpq_lar_core_solver.m_column_types.push_back(column_type::free_column); m_columns_with_changed_bound.increase_size_by_one(); m_rows_with_changed_bounds.increase_size_by_one(); - m_inf_int_set.increase_size_by_one(); add_new_var_to_core_fields_for_mpq(true); if (use_lu) add_new_var_to_core_fields_for_doubles(true); @@ -1720,7 +1687,6 @@ constraint_index lar_solver::add_var_bound(var_index j, lconstraint_kind kind, c add_var_bound_on_constraint_for_term(j, kind, right_side, ci); } lp_assert(sizes_are_correct()); - lp_assert(inf_int_set_is_correct()); return ci; } @@ -1749,7 +1715,7 @@ void lar_solver::update_column_type_and_bound(var_index j, lconstraint_kind kind void lar_solver::add_var_bound_on_constraint_for_term(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) { lp_assert(is_term(j)); unsigned adjusted_term_index = adjust_term_index(j); - lp_assert(!term_is_int(m_terms[adjusted_term_index]) || right_side.is_int()); + // lp_assert(!term_is_int(m_terms[adjusted_term_index]) || right_side.is_int()); auto it = m_ext_vars_to_columns.find(j); if (it != m_ext_vars_to_columns.end()) { unsigned term_j = it->second.internal_j(); @@ -2164,6 +2130,87 @@ var_index lar_solver:: to_var_index(unsigned ext_j) const { lp_assert(it != m_ext_vars_to_columns.end()); return it->second.internal_j(); } + +bool lar_solver::tighten_term_bounds_by_delta(unsigned term_index, const mpq& delta) { + unsigned tj = term_index + m_terms_start_index; + auto it = m_ext_vars_to_columns.find(tj); + if (it == m_ext_vars_to_columns.end()) + return true; + unsigned j = it->second.internal_j(); + auto & slv = m_mpq_lar_core_solver.m_r_solver; + TRACE("cube", tout << "delta = " << delta << std::endl; + m_int_solver->display_column(tout, j); ); + if (slv.column_has_upper_bound(j) && slv.column_has_lower_bound(j)) { + if (slv.m_upper_bounds[j].x - delta < slv.m_lower_bounds[j].x + delta) { + TRACE("cube", tout << "cannot tighten, delta = " << delta;); + return false; + } + } + TRACE("cube", tout << "can tighten";); + if (slv.column_has_upper_bound(j)) { + add_var_bound(tj, lconstraint_kind::LE, slv.m_upper_bounds[j].x - delta); + } + if (slv.column_has_lower_bound(j)) { + add_var_bound(tj, lconstraint_kind::GE, slv.m_lower_bounds[j].x + delta); + } + return true; +} + +void lar_solver::update_delta_for_terms(const impq & delta, unsigned j, const vector& terms_of_var) { + for (unsigned i : terms_of_var) { + lar_term & t = *m_terms[i]; + auto it = t.m_coeffs.find(j); + unsigned tj = to_var_index(i + m_terms_start_index); + TRACE("cube", + tout << "t.apply = " << t.apply(m_mpq_lar_core_solver.m_r_x) << ", m_mpq_lar_core_solver.m_r_x[tj]= " << m_mpq_lar_core_solver.m_r_x[tj];); + TRACE("cube", print_term_as_indices(t, tout); + tout << ", it->second = " << it->second; + tout << ", tj = " << tj << ", "; + m_int_solver->display_column(tout, tj); + ); + + m_mpq_lar_core_solver.m_r_x[tj] += it->second * delta; + lp_assert(t.apply(m_mpq_lar_core_solver.m_r_x) == m_mpq_lar_core_solver.m_r_x[tj]); + TRACE("cube", m_int_solver->display_column(tout, tj); ); + } +} + + +void lar_solver::fill_vars_to_terms(vector> & vars_to_terms) { + for (unsigned j = 0; j < m_terms.size(); j++) { + if (!term_is_used_as_row(j + m_terms_start_index)) + continue; + for (const auto & p : *m_terms[j]) { + if (p.var() >= vars_to_terms.size()) + vars_to_terms.resize(p.var() + 1); + vars_to_terms[p.var()].push_back(j); + } + } +} + +void lar_solver::round_to_integer_solution() { + vector> vars_to_terms; + fill_vars_to_terms(vars_to_terms); + + for (unsigned j = 0; j < column_count(); j++) { + if (column_corresponds_to_term(j)) continue; + TRACE("cube", m_int_solver->display_column(tout, j);); + impq& v = m_mpq_lar_core_solver.m_r_x[j]; + if (v.is_int()) + continue; + impq flv = floor(v); + auto del = flv - v; // del is negative + if (del < - mpq(1, 2)) { + del = impq(one_of_type()) + del; + v = ceil(v); + } else { + v = flv; + } + TRACE("cube", m_int_solver->display_column(tout, j); tout << "v = " << v << " ,del = " << del;); + update_delta_for_terms(del, j, vars_to_terms[j]); + } +} + } // namespace lp diff --git a/src/util/lp/lar_solver.h b/src/util/lp/lar_solver.h index 662d0bb3d..d7850277b 100644 --- a/src/util/lp/lar_solver.h +++ b/src/util/lp/lar_solver.h @@ -39,7 +39,6 @@ Revision History: #include "util/lp/stacked_unordered_set.h" #include "util/lp/implied_bound.h" #include "util/lp/bound_analyzer_on_row.h" -#include "util/lp/quick_xplain.h" #include "util/lp/conversion_helper.h" #include "util/lp/int_solver.h" #include "util/lp/nra_solver.h" @@ -117,17 +116,19 @@ private: vector m_terms; const var_index m_terms_start_index; indexed_vector m_column_buffer; +public: + lar_core_solver m_mpq_lar_core_solver; +private: + int_solver * m_int_solver; + bool m_has_int_var; public : + unsigned terms_start_index() const { return m_terms_start_index; } + const vector terms() const { return m_terms; } const vector& constraints() const { return m_constraints; } - lar_core_solver m_mpq_lar_core_solver; -private: - std::function m_tracker_of_x_change; - int_solver * m_int_solver; -public: void set_int_solver(int_solver * int_slv) { m_int_solver = int_slv; } @@ -136,7 +137,6 @@ public: } unsigned constraint_count() const; const lar_base_constraint& get_constraint(unsigned ci) const; - int_set m_inf_int_set; ////////////////// methods //////////////////////////////// static_matrix> & A_r(); static_matrix> const & A_r() const; @@ -551,43 +551,32 @@ public: } } - bool inf_int_set_is_correct_for_column(unsigned j) const { - if (m_inf_int_set.contains(j) != (column_is_int(j) && (!column_value_is_integer(j)))) { - TRACE("arith_int", - tout << "j= " << j << - " inf_int_set().contains(j) = " << m_inf_int_set.contains(j) << - ", column_is_int(j) = " << column_is_int(j) << - "\n column_value_is_integer(j) = " << column_value_is_integer(j) << - ", val = " << get_column_value(j) << std::endl;); - return false; - } - return true; - } - - bool inf_int_set_is_correct() const { - if (!has_int_var()) - return true; - for (unsigned j = 0; j < A_r().column_count(); j++) { - if (inf_int_set_is_correct_for_column(j) == false) - return false; - } - return true; - } bool has_int_var() const; - void call_assignment_tracker(unsigned j) { - if (!var_is_int(j)) { - lp_assert(m_inf_int_set.contains(j) == false); - return; + bool has_inf_int() const { + for (unsigned j = 0; j < column_count(); j++) { + if (column_is_int(j) && ! column_value_is_int(j)) + return true; } - if (m_mpq_lar_core_solver.m_r_x[j].is_int()) - m_inf_int_set.erase(j); - else - m_inf_int_set.insert(j); + return false; } + bool r_basis_has_inf_int() const { + for (unsigned j : r_basis()) { + if (column_is_int(j) && ! column_value_is_int(j)) + return true; + } + return false; + } + lar_core_solver & get_core_solver() { return m_mpq_lar_core_solver; } bool column_corresponds_to_term(unsigned) const; void catch_up_in_updating_int_solver(); var_index to_var_index(unsigned ext_j) const; + bool tighten_term_bounds_by_delta(unsigned, const mpq&); + void round_to_integer_solution(); + void update_delta_for_terms(const impq & delta, unsigned j, const vector&); + void fill_vars_to_terms(vector> & vars_to_terms); + unsigned column_count() const { return A_r().column_count(); } + const vector & r_basis() const { return m_mpq_lar_core_solver.r_basis(); } }; } diff --git a/src/util/lp/lp_core_solver_base.h b/src/util/lp/lp_core_solver_base.h index f1a02b372..8fd8dc938 100644 --- a/src/util/lp/lp_core_solver_base.h +++ b/src/util/lp/lp_core_solver_base.h @@ -85,12 +85,7 @@ public: bool m_tracing_basis_changes; int_set* m_pivoted_rows; bool m_look_for_feasible_solution_only; - std::function * m_tracker_of_x_change; - void set_tracker_of_x(std::function* tracker) { - m_tracker_of_x_change = tracker; - } - void start_tracing_basis_changes() { m_trace_of_basis_change_vector.resize(0); m_tracing_basis_changes = true; @@ -705,14 +700,10 @@ public: void update_x_and_call_tracker(unsigned j, const X & v) { m_x[j] = v; - if (m_tracker_of_x_change != nullptr) - (*m_tracker_of_x_change)(j); } void add_delta_to_x_and_call_tracker(unsigned j, const X & delta) { m_x[j] += delta; - if (m_tracker_of_x_change != nullptr) - (*m_tracker_of_x_change)(j); } void track_column_feasibility(unsigned j) { @@ -722,14 +713,10 @@ public: insert_column_into_inf_set(j); } void insert_column_into_inf_set(unsigned j) { - if (m_tracker_of_x_change != nullptr) - (*m_tracker_of_x_change)(j); m_inf_set.insert(j); lp_assert(!column_is_feasible(j)); } void remove_column_from_inf_set(unsigned j) { - if (m_tracker_of_x_change != nullptr) - (*m_tracker_of_x_change)(j); m_inf_set.erase(j); lp_assert(column_is_feasible(j)); } diff --git a/src/util/lp/lp_core_solver_base_def.h b/src/util/lp/lp_core_solver_base_def.h index abadabf32..4577b26ec 100644 --- a/src/util/lp/lp_core_solver_base_def.h +++ b/src/util/lp/lp_core_solver_base_def.h @@ -67,8 +67,7 @@ lp_core_solver_base(static_matrix & A, m_steepest_edge_coefficients(A.column_count()), m_tracing_basis_changes(false), m_pivoted_rows(nullptr), - m_look_for_feasible_solution_only(false), - m_tracker_of_x_change(nullptr) { + m_look_for_feasible_solution_only(false) { lp_assert(bounds_for_boxed_are_set_correctly()); init(); init_basis_heading_and_non_basic_columns_vector(); diff --git a/src/util/lp/lp_settings.h b/src/util/lp/lp_settings.h index 460097c58..68f1b31bc 100644 --- a/src/util/lp/lp_settings.h +++ b/src/util/lp/lp_settings.h @@ -108,6 +108,8 @@ struct stats { unsigned m_cut_solver_undef; unsigned m_gcd_calls; unsigned m_gcd_conflicts; + unsigned m_cube_calls; + unsigned m_cube_success; stats() { reset(); } void reset() { memset(this, 0, sizeof(*this)); } }; @@ -334,13 +336,14 @@ public: 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_branch_cut_gomory_threshold; - unsigned m_int_branch_cut_solver; - bool m_run_gcd_test; - unsigned m_cut_solver_cycle_on_var; + unsigned max_row_length_for_bound_propagation; + bool backup_costs; + unsigned column_number_threshold_for_using_lu_in_lar_solver; + unsigned m_int_branch_cut_gomory_threshold; + unsigned m_int_branch_cut_solver; + unsigned m_int_branch_find_cube; + bool m_run_gcd_test; + unsigned m_cut_solver_cycle_on_var; }; // end of lp_settings class diff --git a/src/util/lp/polynomial.h b/src/util/lp/polynomial.h index e88cc4a9c..f4f9615ea 100644 --- a/src/util/lp/polynomial.h +++ b/src/util/lp/polynomial.h @@ -101,6 +101,15 @@ struct polynomial { const mpq & a = coeff(j); return a == 1 || a == -1; } + + template // c plays a role of a map from indices to impq + mpq value(const c& v) const { + mpq r = m_a; + for (auto & p : m_coeffs) + r += v[p.var()].x * p.coeff(); + return r; + } + const vector & coeffs() const { return m_coeffs; } }; } diff --git a/src/util/lp/stacked_vector.h b/src/util/lp/stacked_vector.h index 511c7cf18..539e932b2 100644 --- a/src/util/lp/stacked_vector.h +++ b/src/util/lp/stacked_vector.h @@ -52,6 +52,14 @@ public: bool operator==(B const& other) const { return m_vec.m_vector[m_i] == other; } + B& operator+=(B const &delta) { + // not tracking the change here! + return m_vec.m_vector[m_i] += delta; + } + B& operator-=(B const &delta) { + // not tracking the change here! + return m_vec.m_vector[m_i] -= delta; + } }; class ref_const {