diff --git a/src/math/lp/gomory.cpp b/src/math/lp/gomory.cpp index bc32698ed..5a789a610 100644 --- a/src/math/lp/gomory.cpp +++ b/src/math/lp/gomory.cpp @@ -21,6 +21,8 @@ #include "math/lp/int_solver.h" #include "math/lp/lar_solver.h" #include "math/lp/lp_utils.h" + +#define SMALL_CUTS 0 namespace lp { class gomory::imp { @@ -35,7 +37,11 @@ class gomory::imp { mpq m_one_minus_f; mpq m_fj; mpq m_one_minus_fj; - +#if SMALL_CUTS + mpq m_abs_max, m_big_number; +#endif + struct found_big {}; + const impq & get_value(unsigned j) const { return m_int_solver.get_value(j); } bool is_real(unsigned j) const { return m_int_solver.is_real(j); } bool at_lower(unsigned j) const { return m_int_solver.at_lower(j); } @@ -71,6 +77,9 @@ class gomory::imp { m_t.add_monomial(new_a, j); m_lcm_den = lcm(m_lcm_den, denominator(new_a)); TRACE("gomory_cut_detail", tout << "new_a = " << new_a << ", k = " << m_k << ", lcm_den = " << m_lcm_den << "\n";); +#if SMALL_CUTS + if (numerator(new_a) > m_big_number) throw found_big(); +#endif } void real_case_in_gomory_cut(const mpq & a, unsigned j) { @@ -100,6 +109,9 @@ class gomory::imp { } TRACE("gomory_cut_detail_real", tout << a << "*v" << j << " k: " << m_k << "\n";); m_t.add_monomial(new_a, j); +#if SMALL_CUTS + if (numerator(new_a) > m_big_number) throw found_big(); +#endif } lia_move report_conflict_from_gomory_cut() { @@ -280,6 +292,14 @@ public: tout << "1 - m_f: " << 1 - m_f << ", get_value(m_inf_col).x - m_f = " << get_value(m_inf_col).x - m_f << "\n";); lp_assert(m_f.is_pos() && (get_value(m_inf_col).x - m_f).is_int()); +#if SMALL_CUTS + m_abs_max = 0; + for (const auto & p : m_row) { + mpq t = abs(ceil(p.coeff())); + if (t > m_abs_max) m_abs_max = t; + } + m_big_number = m_abs_max; // .expt(2); +#endif mpq one_min_m_f = 1 - m_f; for (const auto & p : m_row) { unsigned j = p.var(); @@ -290,20 +310,24 @@ public: } // use -p.coeff() to make the format compatible with the format used in: Integrating Simplex with DPLL(T) - if (is_real(j)) { - real_case_in_gomory_cut(- p.coeff(), j); - } else { - if (p.coeff().is_int()) { - // m_fj will be zero and no monomial will be added - continue; + try { + if (is_real(j)) { + real_case_in_gomory_cut(- p.coeff(), j); + } + else if (!p.coeff().is_int()) { + some_int_columns = true; + m_fj = fractional_part(-p.coeff()); + m_one_minus_fj = 1 - m_fj; + int_case_in_gomory_cut(j); } - some_int_columns = true; - m_fj = fractional_part(-p.coeff()); - m_one_minus_fj = 1 - m_fj; - int_case_in_gomory_cut(j); + } + catch (found_big) { + m_ex->clear(); + m_t.clear(); + m_k = 1; + return lia_move::undef; } } - if (m_t.is_empty()) return report_conflict_from_gomory_cut(); if (some_int_columns) diff --git a/src/math/lp/int_solver.cpp b/src/math/lp/int_solver.cpp index 853b24f55..3e78af77e 100644 --- a/src/math/lp/int_solver.cpp +++ b/src/math/lp/int_solver.cpp @@ -12,26 +12,25 @@ namespace lp { -void int_solver::trace_inf_rows() const { - TRACE("arith_int_rows", - unsigned num = m_lar_solver->A_r().column_count(); - for (unsigned v = 0; v < num; v++) { - if (column_is_int(v) && !get_value(v).is_int()) { - display_column(tout, v); - } - } +std::ostream& int_solver::display_inf_rows(std::ostream& out) const { + unsigned num = m_lar_solver->A_r().column_count(); + for (unsigned v = 0; v < num; v++) { + if (column_is_int(v) && !get_value(v).is_int()) { + display_column(out, v); + } + } - num = 0; - for (unsigned i = 0; i < m_lar_solver->A_r().row_count(); i++) { - unsigned j = m_lar_solver->m_mpq_lar_core_solver.m_r_basis[i]; - if (column_is_int_inf(j)) { - num++; - m_lar_solver->print_row(m_lar_solver->A_r().m_rows[i], tout); - tout << "\n"; - } - } - tout << "num of int infeasible: " << num << "\n"; - ); + num = 0; + for (unsigned i = 0; i < m_lar_solver->A_r().row_count(); i++) { + unsigned j = m_lar_solver->m_mpq_lar_core_solver.m_r_basis[i]; + if (column_is_int_inf(j)) { + num++; + m_lar_solver->print_row(m_lar_solver->A_r().m_rows[i], out); + out << "\n"; + } + } + out << "num of int infeasible: " << num << "\n"; + return out; } bool int_solver::has_inf_int() const { @@ -134,7 +133,6 @@ constraint_index int_solver::column_lower_bound_constraint(unsigned j) const { } lia_move int_solver::mk_gomory_cut( unsigned inf_col, const row_strip & row) { - lp_assert(column_is_int_inf(inf_col)); gomory gc(m_t, m_k, m_ex, inf_col, row, *this); return gc.create_cut(); @@ -157,18 +155,6 @@ unsigned int_solver::row_of_basic_column(unsigned j) const { return m_lar_solver->row_of_basic_column(j); } -// template -// void int_solver::fill_cut_solver_for_constraint(constraint_index ci, cut_solver & cs) { -// const lar_base_constraint* c = m_lar_solver->constraints()[ci]; -// vector> coeffs; -// T rs; -// get_int_coeffs_from_constraint(c, coeffs, rs); -// vector explanation; -// explanation.push_back(ci); -// cs.add_ineq(coeffs, -rs, explanation); -// } - - // this will allow to enable and disable tracking of the pivot rows struct check_return_helper { @@ -240,8 +226,12 @@ bool int_solver::tighten_terms_for_cube() { return true; } +bool int_solver::should_find_cube() { + return m_number_of_calls % settings().m_int_find_cube_period == 0; +} + lia_move int_solver::find_cube() { - if (m_number_of_calls % settings().m_int_find_cube_period != 0) + if (!should_find_cube()) return lia_move::undef; settings().stats().m_cube_calls++; @@ -280,8 +270,12 @@ void int_solver::find_feasible_solution() { lp_assert(lp_status::OPTIMAL == m_lar_solver->get_status() || lp_status::FEASIBLE == m_lar_solver->get_status()); } +bool int_solver::should_run_gcd_test() { + return settings().m_int_run_gcd_test; +} + lia_move int_solver::run_gcd_test() { - if (settings().m_int_run_gcd_test) { + if (should_run_gcd_test()) { settings().stats().m_gcd_calls++; TRACE("int_solver", tout << "gcd-test " << settings().stats().m_gcd_calls << "\n";); if (!gcd_test()) { @@ -293,19 +287,19 @@ lia_move int_solver::run_gcd_test() { return lia_move::undef; } +bool int_solver::should_gomory_cut() { + return m_number_of_calls % settings().m_int_gomory_cut_period == 0; +} + lia_move int_solver::gomory_cut() { TRACE("int_solver", tout << "gomory " << m_number_of_calls << "\n";); - if ((m_number_of_calls) % settings().m_int_gomory_cut_period != 0) + if (!should_gomory_cut()) return lia_move::undef; if (m_lar_solver->move_non_basic_columns_to_bounds()) { -#if Z3DEBUG - lp_status st = -#endif - m_lar_solver->find_feasible_solution(); -#if Z3DEBUG + lp_status st = m_lar_solver->find_feasible_solution(); + (void)st; lp_assert(st == lp_status::FEASIBLE || st == lp_status::OPTIMAL); -#endif } int j = find_inf_int_base_column(); @@ -316,7 +310,6 @@ lia_move int_solver::gomory_cut() { return proceed_with_gomory_cut(j); } - void int_solver::try_add_term_to_A_for_hnf(unsigned i) { mpq rs; const lar_term* t = m_lar_solver->terms()[i]; @@ -344,7 +337,7 @@ const lp_settings& int_solver::settings() const { bool int_solver::hnf_has_var_with_non_integral_value() const { for (unsigned j : m_hnf_cutter.vars()) - if (get_value(j).is_int() == false) + if (!get_value(j).is_int()) return true; return false; } @@ -393,14 +386,22 @@ lia_move int_solver::make_hnf_cut() { return r; } +bool int_solver::should_hnf_cut() { + return settings().m_enable_hnf && m_number_of_calls % m_hnf_cut_period == 0; +} + lia_move int_solver::hnf_cut() { - if (!settings().m_enable_hnf) { - return lia_move::undef; + lia_move r = lia_move::undef; + if (should_hnf_cut()) { + r = make_hnf_cut(); + if (r == lia_move::undef) { + m_hnf_cut_period *= 2; + } + else { + m_hnf_cut_period = settings().hnf_cut_period(); + } } - if ((m_number_of_calls) % settings().hnf_cut_period() == 0) { - return make_hnf_cut(); - } - return lia_move::undef; + return r; } lia_move int_solver::check(lp::explanation * e) { @@ -441,12 +442,9 @@ lia_move int_solver::branch_or_sat() { 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(); + return j != -1 ? j : find_inf_int_nbasis_column(); } - void int_solver::set_value_for_nbasic_column_ignore_old_values(unsigned j, const impq & new_val) { lp_assert(!is_base(j)); auto & x = m_lar_solver->m_mpq_lar_core_solver.m_r_x[j]; @@ -455,8 +453,6 @@ void int_solver::set_value_for_nbasic_column_ignore_old_values(unsigned j, const m_lar_solver->change_basic_columns_dependend_on_a_given_nb_column(j, delta); } - - void int_solver::patch_nbasic_column(unsigned j, bool patch_only_int_vals) { auto & lcs = m_lar_solver->m_mpq_lar_core_solver; impq & val = lcs.m_r_x[j]; @@ -694,12 +690,13 @@ linear_combination_iterator * int_solver::get_column_iterator(unsigned j) { int_solver::int_solver(lar_solver* lar_slv) : m_lar_solver(lar_slv), m_number_of_calls(0), - m_hnf_cutter(settings()) { + m_hnf_cutter(settings()), + m_hnf_cut_period(settings().hnf_cut_period()) { m_lar_solver->set_int_solver(this); } -bool int_solver::has_low(unsigned j) const { +bool int_solver::has_lower(unsigned j) const { switch (m_lar_solver->m_mpq_lar_core_solver.m_column_types()[j]) { case column_type::fixed: case column_type::boxed: @@ -722,9 +719,7 @@ bool int_solver::has_upper(unsigned j) const { } -void set_lower(impq & l, - bool & inf_l, - impq const & v ) { +static void set_lower(impq & l, bool & inf_l, impq const & v ) { if (inf_l || v > l) { l = v; inf_l = false; @@ -732,16 +727,14 @@ void set_lower(impq & l, } -void set_upper(impq & u, - bool & inf_u, - impq const & v) { +static void set_upper(impq & u, bool & inf_u, impq const & v) { if (inf_u || v < u) { u = v; inf_u = false; } } - bool int_solver::get_freedom_interval_for_column(unsigned j, bool val_is_int, bool & inf_l, impq & l, bool & inf_u, impq & u, mpq & m) { +bool int_solver::get_freedom_interval_for_column(unsigned j, bool val_is_int, bool & inf_l, impq & l, bool & inf_u, impq & u, mpq & m) { auto & lcs = m_lar_solver->m_mpq_lar_core_solver; if (lcs.m_r_heading[j] >= 0) // the basic var return false; @@ -753,7 +746,7 @@ void set_upper(impq & u, l = u = zero_of_type(); m = mpq(1); - if (has_low(j)) { + if (has_lower(j)) { set_lower(l, inf_l, lower_bound(j) - xj); } if (has_upper(j)) { @@ -794,7 +787,7 @@ void set_upper(impq & u, if (a.is_neg()) { - if (has_low(i)) { + if (has_lower(i)) { SET_BOUND(set_lower, l, inf_l, a, xi, lcs.m_r_lower_bounds()[i]); } if (has_upper(i)) { @@ -805,7 +798,7 @@ void set_upper(impq & u, if (has_upper(i)) { SET_BOUND(set_lower, l, inf_l, a, xi, lcs.m_r_upper_bounds()[i]); } - if (has_low(i)) { + if (has_lower(i)) { SET_BOUND(set_upper, u, inf_u, a, xi, lcs.m_r_lower_bounds()[i]); } } @@ -839,9 +832,7 @@ bool int_solver::is_real(unsigned j) const { bool int_solver::value_is_int(unsigned j) const { return m_lar_solver->column_value_is_int(j); -} - - +} bool int_solver::is_feasible() const { auto & lcs = m_lar_solver->m_mpq_lar_core_solver; @@ -1031,18 +1022,16 @@ lia_move int_solver::create_branch_on_column(int j) { lp_assert(j != -1); m_t.add_monomial(mpq(1), m_lar_solver->adjust_column_index_to_term_index(j)); if (is_free(j)) { - m_upper = true; + m_upper = random() % 2; m_k = mpq(0); } else { m_upper = left_branch_is_more_narrow_than_right(j); m_k = m_upper? floor(get_value(j)) : ceil(get_value(j)); } - TRACE("int_solver", tout << "branching v" << j << " = " << get_value(j) << "\n"; - display_column(tout, j); - tout << "k = " << m_k << std::endl; - - ); + TRACE("int_solver", + display_column(tout << "branching v" << j << " = " << get_value(j) << "\n", j); + tout << "k = " << m_k << std::endl;); return lia_move::branch; } diff --git a/src/math/lp/int_solver.h b/src/math/lp/int_solver.h index 042491272..df60afee3 100644 --- a/src/math/lp/int_solver.h +++ b/src/math/lp/int_solver.h @@ -44,6 +44,7 @@ public: 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 hnf_cutter m_hnf_cutter; + unsigned m_hnf_cut_period; // methods int_solver(lar_solver* lp); @@ -97,7 +98,12 @@ private: bool non_basic_columns_are_at_bounds() const; bool is_feasible() const; bool column_is_int_inf(unsigned j) const; - void trace_inf_rows() const; + std::ostream& display_inf_rows(std::ostream&) const; + bool should_find_cube(); + bool should_run_gcd_test(); + bool should_gomory_cut(); + bool should_hnf_cut(); + lia_move branch_or_sat(); int find_any_inf_int_column_basis_first(); int find_inf_int_base_column(); @@ -110,7 +116,7 @@ private: lia_move proceed_with_gomory_cut(unsigned j); bool is_gomory_cut_target(const row_strip&); bool at_bound(unsigned j) const; - bool has_low(unsigned j) const; + bool has_lower(unsigned j) const; bool has_upper(unsigned j) const; unsigned row_of_basic_column(unsigned j) const;