From 3de5af3cb0ad8b2481598a61065dbdd51525daba Mon Sep 17 00:00:00 2001 From: Nikolaj Bjorner Date: Fri, 10 Nov 2023 16:38:55 +0100 Subject: [PATCH] fix bug in bound simplification in Gomory for mixed integer linear cuts, enable fixed variable redution after bugfix, add notes that rewriting bounds does not work Signed-off-by: Nikolaj Bjorner --- src/math/lp/gomory.cpp | 209 +++++++++++++++---------------------- src/math/lp/hnf_cutter.cpp | 40 +++---- src/math/lp/int_solver.cpp | 116 ++------------------ src/math/lp/int_solver.h | 6 +- src/smt/theory_lra.cpp | 11 +- 5 files changed, 121 insertions(+), 261 deletions(-) diff --git a/src/math/lp/gomory.cpp b/src/math/lp/gomory.cpp index cc752c93194..2c096320a69 100644 --- a/src/math/lp/gomory.cpp +++ b/src/math/lp/gomory.cpp @@ -32,7 +32,6 @@ class create_cut { unsigned m_inf_col; // a basis column which has to be an integer but has a non integral value const row_strip& m_row; int_solver& lia; - mpq m_lcm_den = { mpq(1) }; mpq m_f; mpq m_one_minus_f; mpq m_fj; @@ -131,120 +130,7 @@ class create_cut { // conflict 0 >= k where k is positive return lia_move::conflict; } - - void divd(mpq& r, mpq const& d) { - r /= d; - if (!r.is_int()) - r = ceil(r); - } - - bool can_divide_by(vector> const& p, mpq const& d) { - mpq lhs(0), rhs(m_k); - mpq max_c(abs(m_k)); - for (auto const& [c, v] : p) { - auto c1 = c; - max_c = std::max(max_c, abs(c1)); - divd(c1, d); - if (c1 == 0) - return false; - VERIFY(lia.get_value(v).y == 0); - lhs += c1 * lia.get_value(v).x; - } - if (max_c == 1) - return false; - divd(rhs, d); - return lhs < rhs; - } - void adjust_term_and_k_for_some_ints_case_gomory() { - lp_assert(!m_t.is_empty()); - // k = 1 + sum of m_t at bounds - lar_term t = lia.lra.unfold_nested_subterms(m_t); - auto pol = t.coeffs_as_vector(); - m_t.clear(); - if (pol.size() == 1) { - TRACE("gomory_cut_detail", tout << "pol.size() is 1" << std::endl;); - auto const& [a, v] = pol[0]; - lp_assert(is_int(v)); - if (a.is_pos()) { // we have av >= k - divd(m_k, a); - m_t.add_monomial(mpq(1), v); - } - else { - // av >= k - // a/-a*v >= k / - a - // -v >= k / - a - // -v >= ceil(k / -a) - divd(m_k, -a); - m_t.add_monomial(-mpq(1), v); - } - } - else { - m_lcm_den = denominator(m_k); - for (auto const& [c, v] : pol) - m_lcm_den = lcm(m_lcm_den, denominator(c)); - lp_assert(m_lcm_den.is_pos()); - TRACE("gomory_cut_detail", tout << "pol.size() > 1 den: " << m_lcm_den << std::endl;); - if (!m_lcm_den.is_one()) { - // normalize coefficients of integer parameters to be integers. - for (auto & [c,v]: pol) { - c *= m_lcm_den; - SASSERT(!is_int(v) || c.is_int()); - } - m_k *= m_lcm_den; - } -#if 0 - unsigned j = 0, i = 0; - for (auto & [c, v] : pol) { - if (lia.is_fixed(v)) { - push_explanation(column_lower_bound_constraint(v)); - push_explanation(column_upper_bound_constraint(v)); - m_k -= c; - IF_VERBOSE(0, verbose_stream() << "got fixed " << v << "\n"); - } - else - pol[j++] = pol[i]; - ++i; - } - pol.shrink(j); -#endif - - // gcd reduction is loss-less: - mpq g(1); - for (const auto & [c, v] : pol) - g = gcd(g, c); - - if (g != 1) { - for (auto & [c, v] : pol) - c /= g; - divd(m_k, g); - } - -#if 0 - // TODO: create self-contained rounding mode to weaken cuts - // whose cofficients are considered too large - // (larger than bounds from the input) - mpq min_c = abs(m_k); - for (const auto & [c, v] : pol) - min_c = std::min(min_c, abs(c)); - - if (min_c > 1 && can_divide_by(pol, min_c)) { - for (auto& [c, v] : pol) - divd(c, min_c); - divd(m_k, min_c); - } -#endif - - for (const auto & [c, v]: pol) - m_t.add_monomial(c, v); - VERIFY(m_t.size() > 0); - } - - TRACE("gomory_cut_detail", tout << "k = " << m_k << std::endl;); - lp_assert(m_k.is_int()); - - - } std::string var_name(unsigned j) const { return std::string("x") + std::to_string(j); @@ -359,12 +245,12 @@ class create_cut { // gomory will be t >= k and the current solution has a property t < k m_k = 1; m_t.clear(); - bool some_int_columns = false; mpq m_f = fractional_part(get_value(m_inf_col)); TRACE("gomory_cut_detail", tout << "m_f: " << m_f << ", "; 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()); + bool some_int_columns = false; #if SMALL_CUTS m_abs_max = 0; for (const auto & p : m_row) { @@ -408,13 +294,7 @@ class create_cut { if (m_t.is_empty()) return report_conflict_from_gomory_cut(); if (some_int_columns) - adjust_term_and_k_for_some_ints_case_gomory(); - if (!lia.current_solution_is_inf_on_cut()) { - m_ex->clear(); - m_t.clear(); - m_k = 1; - return lia_move::undef; - } + simplify_inequality(); lp_assert(lia.current_solution_is_inf_on_cut()); // checks that indices are columns TRACE("gomory_cut", print_linear_combination_of_column_indices_only(m_t.coeffs_as_vector(), tout << "gomory cut: "); tout << " >= " << m_k << std::endl;); TRACE("gomory_cut_detail", dump_cut_and_constraints_as_smt_lemma(tout); @@ -423,6 +303,91 @@ class create_cut { return lia_move::cut; } + // TODO: use this also for HNF cuts? + mpq m_lcm_den = { mpq(1) }; + + void simplify_inequality() { + + auto divd = [](mpq& r, mpq const& d) { + r /= d; + if (!r.is_int()) + r = ceil(r); + }; + SASSERT(!lia.m_upper); + lp_assert(!m_t.is_empty()); + // k = 1 + sum of m_t at bounds + lar_term t = lia.lra.unfold_nested_subterms(m_t); + auto pol = t.coeffs_as_vector(); + m_t.clear(); + if (pol.size() == 1) { + TRACE("gomory_cut_detail", tout << "pol.size() is 1" << std::endl;); + auto const& [a, v] = pol[0]; + lp_assert(is_int(v)); + if (a.is_pos()) { // we have av >= k + divd(m_k, a); + m_t.add_monomial(mpq(1), v); + } + else { + // av >= k + // a/-a*v >= k / - a + // -v >= k / - a + // -v >= ceil(k / -a) + divd(m_k, -a); + m_t.add_monomial(-mpq(1), v); + } + } + else { + m_lcm_den = denominator(m_k); + for (auto const& [c, v] : pol) + m_lcm_den = lcm(m_lcm_den, denominator(c)); + lp_assert(m_lcm_den.is_pos()); + bool int_row = true; + TRACE("gomory_cut_detail", tout << "pol.size() > 1 den: " << m_lcm_den << std::endl;); + if (!m_lcm_den.is_one()) { + // normalize coefficients of integer parameters to be integers. + for (auto & [c,v]: pol) { + c *= m_lcm_den; + SASSERT(!is_int(v) || c.is_int()); + int_row &= is_int(v); + } + m_k *= m_lcm_den; + } + unsigned j = 0, i = 0; + for (auto & [c, v] : pol) { + if (lia.is_fixed(v)) { + push_explanation(column_lower_bound_constraint(v)); + push_explanation(column_upper_bound_constraint(v)); + m_k -= c * lower_bound(v).x; + } + else + pol[j++] = pol[i]; + ++i; + } + pol.shrink(j); + + // gcd reduction is loss-less: + mpq g(1); + for (const auto & [c, v] : pol) + g = gcd(g, c); + if (!int_row) + g = gcd(g, m_k); + + if (g != 1) { + for (auto & [c, v] : pol) + c /= g; + divd(m_k, g); + } + + for (const auto & [c, v]: pol) + m_t.add_monomial(c, v); + VERIFY(m_t.size() > 0); + } + + TRACE("gomory_cut_detail", tout << "k = " << m_k << std::endl;); + lp_assert(m_k.is_int()); + } + + create_cut(lar_term & t, mpq & k, explanation* ex, unsigned basic_inf_int_j, const row_strip& row, int_solver& lia) : m_t(t), m_k(k), diff --git a/src/math/lp/hnf_cutter.cpp b/src/math/lp/hnf_cutter.cpp index d688eeb63a6..75095ca515e 100644 --- a/src/math/lp/hnf_cutter.cpp +++ b/src/math/lp/hnf_cutter.cpp @@ -83,14 +83,13 @@ namespace lp { // consider return from here if b[i] is not an integer and return i } } - + vector hnf_cutter::create_b(const svector & basis_rows) { if (basis_rows.size() == m_right_sides.size()) return m_right_sides; vector b; - for (unsigned i : basis_rows) { - b.push_back(m_right_sides[i]); - } + for (unsigned i : basis_rows) + b.push_back(m_right_sides[i]); return b; } @@ -98,16 +97,15 @@ namespace lp { int ret = -1; int n = 0; for (int i = 0; i < static_cast(b.size()); i++) { - if (is_integer(b[i])) continue; - if (n == 0 ) { + if (is_integer(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_settings.random_next() % (++n) == 0) + ret = i; } return ret; } @@ -240,10 +238,7 @@ branch y_i >= ceil(y0_i) is impossible. } bool hnf_cutter::hnf_has_var_with_non_integral_value() const { - for (unsigned j : vars()) - if (!lia.get_value(j).is_int()) - return true; - return false; + return any_of(vars(), [&](unsigned j) { return !lia.get_value(j).is_int(); }); } bool hnf_cutter::init_terms_for_hnf_cut() { @@ -252,17 +247,15 @@ branch y_i >= ceil(y0_i) is impossible. try_add_term_to_A_for_hnf(tv::term(i)); return hnf_has_var_with_non_integral_value(); } - + lia_move hnf_cutter::make_hnf_cut() { - if (!init_terms_for_hnf_cut()) { + if (!init_terms_for_hnf_cut()) return lia_move::undef; - } lia.settings().stats().m_hnf_cutter_calls++; TRACE("hnf_cut", tout << "settings().stats().m_hnf_cutter_calls = " << lia.settings().stats().m_hnf_cutter_calls << "\n"; - for (u_dependency* d : constraints_for_explanation()) { - for (auto ci : lra.flatten(d)) - lra.constraints().display(tout, ci); - } + for (u_dependency* d : constraints_for_explanation()) + for (auto ci : lra.flatten(d)) + lra.constraints().display(tout, ci); tout << lra.constraints(); ); #ifdef Z3DEBUG @@ -273,14 +266,14 @@ branch y_i >= ceil(y0_i) is impossible. , x0 #endif ); - + if (r == lia_move::cut) { TRACE("hnf_cut", lra.print_term(lia.m_t, tout << "cut:"); tout << " <= " << lia.m_k << std::endl; for (auto* dep : constraints_for_explanation()) for (auto ci : lra.flatten(dep)) - lra.constraints().display(tout, ci); + lra.constraints().display(tout, ci); ); lp_assert(lia.current_solution_is_inf_on_cut()); lia.settings().stats().m_hnf_cuts++; @@ -291,5 +284,4 @@ branch y_i >= ceil(y0_i) is impossible. } return r; } - } diff --git a/src/math/lp/int_solver.cpp b/src/math/lp/int_solver.cpp index 602bb0a39f4..c06c891b0ac 100644 --- a/src/math/lp/int_solver.cpp +++ b/src/math/lp/int_solver.cpp @@ -197,11 +197,14 @@ namespace lp { if (r == lia_move::undef && should_find_cube()) r = int_cube(*this)(); if (r == lia_move::undef) lra.move_non_basic_columns_to_bounds(); if (r == lia_move::undef && should_hnf_cut()) r = hnf_cut(); + + std::function gomory_fn = [&]() { return gomory(*this)(); }; m_cut_vars.reset(); #if 0 if (r == lia_move::undef && should_gomory_cut()) r = gomory(*this)(); #else - if (r == lia_move::undef && should_gomory_cut()) r = local_gomory(2); + if (r == lia_move::undef && should_gomory_cut()) r = local_cut(2, gomory_fn); + #endif m_cut_vars.reset(); if (r == lia_move::undef) r = int_branch(*this)(); @@ -711,6 +714,8 @@ namespace lp { return; +#if 0 + // in-processing simplification can go here, such as bounds improvements. if (!lra.is_feasible()) { @@ -728,113 +733,9 @@ namespace lp { if (has_inf_int()) local_gomory(5); -#if 0 stopwatch sw; explanation exp1, exp2; - // - // tighten integer bounds - // It is a weak method as it ownly strengthens bounds if - // variables are already at one of the to-be discovered bounds. - // - sw.start(); - unsigned changes = 0; - auto const& constraints = lra.constraints(); - auto print_var = [this](lpvar j) { - if (lra.column_corresponds_to_term(j)) { - std::stringstream strm; - lra.print_column_info(j, strm); - return strm.str(); - } - else - return std::string("j") + std::to_string(j); - }; - unsigned start = random(); - unsigned num_checks = 0; - for (lpvar j0 = 0; j0 < lra.column_count(); ++j0) { - lpvar j = (j0 + start) % lra.column_count(); - - if (num_checks > 1000) - break; - if (is_fixed(j)) - continue; - if (!lra.column_is_int(j)) - continue; - rational value = get_value(j).x; - bool tight_lower = false, tight_upper = false; - u_dependency* dep; - - if (!value.is_int()) - continue; - - bool at_up = at_upper(j); - - if (!at_lower(j)) { - ++num_checks; - lra.push(); - auto k = lp::lconstraint_kind::LE; - lra.update_column_type_and_bound(j, k, (value - 1).to_mpq(), nullptr); - lra.find_feasible_solution(); - if (!lra.is_feasible()) { - tight_upper = true; - ++changes; - lra.get_infeasibility_explanation(exp1); -#if 0 - display_column(std::cout, j); - std::cout << print_var(j) << " >= " << value << "\n"; - unsigned i = 0; - for (auto p : exp1) { - std::cout << "(" << p.ci() << ")"; - constraints.display(std::cout, print_var, p.ci()); - if (++i < exp1.size()) - std::cout << " "; - } -#endif - } - lra.pop(1); - if (tight_upper) { - dep = nullptr; - for (auto& cc : exp1) - dep = lra.join_deps(dep, constraints[cc.ci()].dep()); - lra.update_column_type_and_bound(j, lp::lconstraint_kind::GE, value.to_mpq(), dep); - } - } - - if (!at_up) { - ++num_checks; - lra.push(); - auto k = lp::lconstraint_kind::GE; - lra.update_column_type_and_bound(j, k, (value + 1).to_mpq(), nullptr); - lra.find_feasible_solution(); - if (!lra.is_feasible()) { - tight_lower = true; - ++changes; - lra.get_infeasibility_explanation(exp1); -#if 0 - display_column(std::cout, j); - std::cout << print_var(j) << " <= " << value << "\n"; - unsigned i = 0; - for (auto p : exp1) { - std::cout << "(" << p.ci() << ")"; - constraints.display(std::cout, print_var, p.ci()); - if (++i < exp1.size()) - std::cout << " "; - } -#endif - } - lra.pop(1); - if (tight_lower) { - dep = nullptr; - for (auto& cc : exp1) - dep = lra.join_deps(dep, constraints[cc.ci()].dep()); - lra.update_column_type_and_bound(j, lp::lconstraint_kind::LE, value.to_mpq(), dep); - } - } - } - sw.stop(); - std::cout << "changes " << changes << " columns " << lra.column_count() << " time: " << sw.get_seconds() << "\n"; - std::cout.flush(); - // // identify equalities // @@ -948,7 +849,8 @@ namespace lp { #endif } - lia_move int_solver::local_gomory(unsigned num_cuts) { + + lia_move int_solver::local_cut(unsigned num_cuts, std::function& cut_fn) { struct ex { explanation m_ex; lar_term m_term; mpq m_k; bool m_is_upper; }; vector cuts; @@ -956,7 +858,7 @@ namespace lp { m_ex->clear(); m_t.clear(); m_k.reset(); - auto r = gomory(*this)(); + auto r = cut_fn(); if (r != lia_move::cut) break; cuts.push_back({ *m_ex, m_t, m_k, is_upper() }); diff --git a/src/math/lp/int_solver.h b/src/math/lp/int_solver.h index 5a7d253d882..0832bd45e4b 100644 --- a/src/math/lp/int_solver.h +++ b/src/math/lp/int_solver.h @@ -68,7 +68,7 @@ class int_solver { hnf_cutter m_hnf_cutter; unsigned m_hnf_cut_period; unsigned_vector m_cut_vars; // variables that should not be selected for cuts - + vector m_equalities; public: int_solver(lar_solver& lp); @@ -111,8 +111,8 @@ class int_solver { bool has_upper(unsigned j) const; unsigned row_of_basic_column(unsigned j) const; bool cut_indices_are_columns() const; - lia_move local_gomory(unsigned num_cuts); - + lia_move local_cut(unsigned num_cuts, std::function& cut_fn); + public: std::ostream& display_column(std::ostream & out, unsigned j) const; u_dependency* column_upper_bound_constraint(unsigned j) const; diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp index 10ec08247d3..377db4fc2bc 100644 --- a/src/smt/theory_lra.cpp +++ b/src/smt/theory_lra.cpp @@ -1784,14 +1784,15 @@ class theory_lra::imp { expr_ref atom(m); t = coeffs2app(coeffs, rational::zero(), is_int); - if (lower_bound) { + if (lower_bound) atom = a.mk_ge(t, a.mk_numeral(offset, is_int)); - } - else { - atom = a.mk_le(t, a.mk_numeral(offset, is_int)); - } + else + atom = a.mk_le(t, a.mk_numeral(offset, is_int)); // ctx().get_rewriter()(atom); + // Note: it is not safe to rewrite atom because the rewriter can + // destroy structure, such as (div x 24) >= 0 becomes x >= 0 and the internal variable + // corresponding to (div x 24) is not constrained. TRACE("arith", tout << t << ": " << atom << "\n"; lp().print_term(term, tout << "bound atom: ") << (lower_bound?" >= ":" <= ") << k << "\n";); ctx().internalize(atom, true);