From db8f01894f14fdbfdbad3cdcdb9398ffaa526e60 Mon Sep 17 00:00:00 2001 From: Lev Nachmanson Date: Thu, 27 Jul 2017 10:49:00 -0700 Subject: [PATCH] solve send-more-money_lev.smt2 Signed-off-by: Lev Nachmanson handle integer vars in random_update Signed-off-by: Lev Nachmanson call the assert in gomory_cut and branching to a correct place Signed-off-by: Lev Nachmanson fixes in goromy cut Signed-off-by: Lev Nachmanson disable x values tracking in random_update Signed-off-by: Lev Nachmanson more fixes in gomory cut Signed-off-by: Lev Nachmanson change in mk_bound by Nikolaj Signed-off-by: Lev Nachmanson fixes in gomory cut and setup Signed-off-by: Lev Nachmanson fixes in int_solver Signed-off-by: Lev Nachmanson change a printout Signed-off-by: Lev Nachmanson fix by Nikolaj in treating terms returned by int_solver Signed-off-by: Lev Nachmanson fix syntax Signed-off-by: Lev Nachmanson fix a free coefficient bug in bound propagaion and simplify gomory cut Signed-off-by: Lev Nachmanson avoid tracking pivoted rows during int_solver::check() --- src/ast/arith_decl_plugin.h | 4 +- src/smt/smt_setup.cpp | 21 +- src/smt/theory_arith_int.h | 3 +- src/smt/theory_arith_pp.h | 15 +- src/smt/theory_lra.cpp | 220 ++++-- src/smt/theory_str.cpp | 4 +- src/test/lp/lp.cpp | 2 +- src/util/lp/bound_propagator.cpp | 19 +- src/util/lp/column_namer.h | 4 +- src/util/lp/implied_bound.h | 4 +- src/util/lp/int_solver.cpp | 653 +++++++++++------- src/util/lp/int_solver.h | 47 +- src/util/lp/iterator_on_column.h | 4 +- src/util/lp/iterator_on_indexed_vector.h | 5 +- src/util/lp/iterator_on_pivot_row.h | 5 +- src/util/lp/iterator_on_row.h | 5 +- src/util/lp/iterator_on_term_with_basis_var.h | 5 +- src/util/lp/lar_core_solver.h | 5 +- src/util/lp/lar_core_solver.hpp | 18 +- src/util/lp/lar_solver.cpp | 95 +-- src/util/lp/lar_solver.h | 171 ++--- src/util/lp/lar_term.h | 6 +- src/util/lp/linear_combination_iterator.h | 3 +- src/util/lp/lp_bound_propagator.h | 2 +- src/util/lp/lp_core_solver_base.h | 78 ++- src/util/lp/lp_primal_core_solver.h | 2 +- src/util/lp/lp_primal_core_solver.hpp | 3 +- src/util/lp/lp_primal_core_solver_tableau.h | 11 +- src/util/lp/lp_settings.h | 10 +- src/util/lp/random_updater.h | 64 +- src/util/lp/random_updater.hpp | 173 +---- 31 files changed, 894 insertions(+), 767 deletions(-) diff --git a/src/ast/arith_decl_plugin.h b/src/ast/arith_decl_plugin.h index d7340297b..341b3fb12 100644 --- a/src/ast/arith_decl_plugin.h +++ b/src/ast/arith_decl_plugin.h @@ -259,9 +259,9 @@ public: bool is_uminus(expr const * n) const { return is_app_of(n, m_afid, OP_UMINUS); } bool is_mul(expr const * n) const { return is_app_of(n, m_afid, OP_MUL); } bool is_div(expr const * n) const { return is_app_of(n, m_afid, OP_DIV); } - //bool is_div0(expr const * n) const { return is_app_of(n, m_afid, OP_DIV_0); } + bool is_div0(expr const * n) const { return is_app_of(n, m_afid, OP_DIV_0); } bool is_idiv(expr const * n) const { return is_app_of(n, m_afid, OP_IDIV); } - //bool is_idiv0(expr const * n) const { return is_app_of(n, m_afid, OP_IDIV_0); } + bool is_idiv0(expr const * n) const { return is_app_of(n, m_afid, OP_IDIV_0); } bool is_mod(expr const * n) const { return is_app_of(n, m_afid, OP_MOD); } bool is_rem(expr const * n) const { return is_app_of(n, m_afid, OP_REM); } bool is_to_real(expr const * n) const { return is_app_of(n, m_afid, OP_TO_REAL); } diff --git a/src/smt/smt_setup.cpp b/src/smt/smt_setup.cpp index b90f08cd3..d2571de91 100644 --- a/src/smt/smt_setup.cpp +++ b/src/smt/smt_setup.cpp @@ -730,8 +730,12 @@ namespace smt { } void setup::setup_i_arith() { - // m_context.register_plugin(alloc(smt::theory_lra, m_manager, m_params)); - m_context.register_plugin(alloc(smt::theory_i_arith, m_manager, m_params)); + if (AS_LRA == m_params.m_arith_mode) { + setup_r_arith(); + } + else { + m_context.register_plugin(alloc(smt::theory_i_arith, m_manager, m_params)); + } } void setup::setup_r_arith() { @@ -739,14 +743,21 @@ namespace smt { } void setup::setup_mi_arith() { - if (m_params.m_arith_mode == AS_OPTINF) { + switch (m_params.m_arith_mode) { + case AS_OPTINF: m_context.register_plugin(alloc(smt::theory_inf_arith, m_manager, m_params)); - } - else { + break; + case AS_LRA: + setup_r_arith(); + break; + default: m_context.register_plugin(alloc(smt::theory_mi_arith, m_manager, m_params)); + break; } } + + void setup::setup_arith() { static_features st(m_manager); IF_VERBOSE(100, verbose_stream() << "(smt.collecting-features)\n";); diff --git a/src/smt/theory_arith_int.h b/src/smt/theory_arith_int.h index b5281f218..63b8e3cc1 100644 --- a/src/smt/theory_arith_int.h +++ b/src/smt/theory_arith_int.h @@ -1385,6 +1385,7 @@ namespace smt { m_branch_cut_counter++; // TODO: add giveup code + TRACE("gomory_cut", tout << m_branch_cut_counter << ", " << m_params.m_arith_branch_cut_ratio << std::endl;); if (m_branch_cut_counter % m_params.m_arith_branch_cut_ratio == 0) { TRACE("opt_verbose", display(tout);); move_non_base_vars_to_bounds(); @@ -1399,7 +1400,7 @@ namespace smt { SASSERT(is_base(int_var)); row const & r = m_rows[get_var_row(int_var)]; if (!mk_gomory_cut(r)) { - // silent failure + TRACE("gomory_cut", tout << "silent failure\n";); } return FC_CONTINUE; } diff --git a/src/smt/theory_arith_pp.h b/src/smt/theory_arith_pp.h index a218445f5..6381da9f1 100644 --- a/src/smt/theory_arith_pp.h +++ b/src/smt/theory_arith_pp.h @@ -389,8 +389,19 @@ namespace smt { void theory_arith::display_vars(std::ostream & out) const { out << "vars:\n"; int n = get_num_vars(); - for (theory_var v = 0; v < n; v++) - display_var(out, v); + int inf_vars = 0; + int int_inf_vars = 0; + for (theory_var v = 0; v < n; v++) { + if ((lower(v) && lower(v)->get_value() > get_value(v)) + || (upper(v) && upper(v)->get_value() < get_value(v))) + inf_vars++; + if (is_int(v) && !get_value(v).is_int()) + int_inf_vars++; + } + out << "infeasibles = " << inf_vars << " int_inf = " << int_inf_vars << std::endl; + for (theory_var v = 0; v < n; v++) { + display_var(out, v); + } } template diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp index 894844008..e8ab2b88e 100644 --- a/src/smt/theory_lra.cpp +++ b/src/smt/theory_lra.cpp @@ -77,6 +77,7 @@ namespace lp_api { if (is_true) return inf_rational(m_value); // v >= value or v <= value if (m_is_int) { + SASSERT(m_value.is_int()); if (m_bound_kind == lower_t) return inf_rational(m_value - rational::one()); // v <= value - 1 return inf_rational(m_value + rational::one()); // v >= value + 1 } @@ -86,7 +87,7 @@ namespace lp_api { } } virtual std::ostream& display(std::ostream& out) const { - return out << "v" << get_var() << " " << get_bound_kind() << " " << m_value; + return out << m_value << " " << get_bound_kind() << " v" << get_var(); } }; @@ -303,7 +304,10 @@ namespace smt { m_solver->settings().simplex_strategy() = static_cast(lp.simplex_strategy()); reset_variable_values(); m_solver->settings().bound_propagation() = BP_NONE != propagation_mode(); - m_solver->set_propagate_bounds_on_pivoted_rows_mode(lp.bprop_on_pivoted_rows()); + m_solver->set_track_pivoted_rows(lp.bprop_on_pivoted_rows()); + m_solver->settings().m_int_branch_cut_gomory_threshold = ctx().get_fparams().m_arith_branch_cut_ratio; + m_solver->settings().m_run_gcd_test = ctx().get_fparams().m_arith_gcd_test; + m_solver->settings().set_random_seed(ctx().get_fparams().m_random_seed); //m_solver->settings().set_ostream(0); m_lia = alloc(lp::int_solver, m_solver.get()); } @@ -318,6 +322,9 @@ namespace smt { } void found_not_handled(expr* n) { + if (a.is_div0(n)) { + return; + } m_not_handled = n; if (is_app(n) && is_underspecified(to_app(n))) { TRACE("arith", tout << "Unhandled: " << mk_pp(n, m) << "\n";); @@ -608,8 +615,8 @@ namespace smt { } bool all_zeros(vector const& v) const { - for (unsigned i = 0; i < v.size(); ++i) { - if (!v[i].is_zero()) { + for (rational const& r : v) { + if (!r.is_zero()) { return false; } } @@ -708,7 +715,7 @@ namespace smt { m_var_index2theory_var.setx(vi, v, UINT_MAX); } m_var_trail.push_back(v); - TRACE("arith", tout << "v" << v << " := " << mk_pp(term, m) << " slack: " << vi << " scopes: " << m_scopes.size() << "\n"; + TRACE("arith_verbose", tout << "v" << v << " := " << mk_pp(term, m) << " slack: " << vi << " scopes: " << m_scopes.size() << "\n"; m_solver->print_term(m_solver->get_term(vi), tout); tout << "\n";); } rational val; @@ -777,7 +784,7 @@ namespace smt { updt_unassigned_bounds(v, +1); m_bounds_trail.push_back(v); m_bool_var2bound.insert(bv, b); - TRACE("arith", tout << "Internalized " << mk_pp(atom, m) << "\n";); + TRACE("arith_verbose", tout << "Internalized " << mk_pp(atom, m) << "\n";); mk_bound_axioms(*b); //add_use_lists(b); return true; @@ -801,7 +808,7 @@ namespace smt { if (n1->get_th_var(get_id()) != null_theory_var && n2->get_th_var(get_id()) != null_theory_var && n1 != n2) { - TRACE("arith", tout << mk_pp(atom, m) << "\n";); + TRACE("arith_verbose", tout << mk_pp(atom, m) << "\n";); m_arith_eq_adapter.mk_axioms(n1, n2); } } @@ -1037,38 +1044,58 @@ namespace smt { return m_solver->var_is_registered(m_theory_var2var_index[v]); } + mutable vector> m_todo_terms; + lp::impq get_ivalue(theory_var v) const { lp_assert(can_get_ivalue(v)); lp::var_index vi = m_theory_var2var_index[v]; if (!m_solver->is_term(vi)) - return m_solver->get_value(vi); - - const lp::lar_term& term = m_solver->get_term(vi); - lp::impq result(term.m_v); - for (const auto & i: term.m_coeffs) { - result += m_solver->get_value(i.first) * i.second; + return m_solver->get_column_value(vi); + m_todo_terms.push_back(std::make_pair(vi, rational::one())); + lp::impq result(0); + while (!m_todo_terms.empty()) { + vi = m_todo_terms.back().first; + rational coeff = m_todo_terms.back().second; + m_todo_terms.pop_back(); + if (m_solver->is_term(vi)) { + const lp::lar_term& term = m_solver->get_term(vi); + result += term.m_v * coeff; + for (const auto & i: term.m_coeffs) { + m_todo_terms.push_back(std::make_pair(i.first, coeff * i.second)); + } + } + else { + result += m_solver->get_column_value(vi) * coeff; + } } return result; } - rational get_value(theory_var v) const { if (!can_get_value(v)) return rational::zero(); lp::var_index vi = m_theory_var2var_index[v]; if (m_variable_values.count(vi) > 0) { return m_variable_values[vi]; } - if (m_solver->is_term(vi)) { - const lp::lar_term& term = m_solver->get_term(vi); - rational result = term.m_v; - for (auto i = term.m_coeffs.begin(); i != term.m_coeffs.end(); ++i) { - result += m_variable_values[i->first] * i->second; + m_todo_terms.push_back(std::make_pair(vi, rational::one())); + rational result(0); + while (!m_todo_terms.empty()) { + lp::var_index wi = m_todo_terms.back().first; + rational coeff = m_todo_terms.back().second; + m_todo_terms.pop_back(); + if (m_solver->is_term(wi)) { + const lp::lar_term& term = m_solver->get_term(wi); + result += term.m_v * coeff; + for (auto const& i : term.m_coeffs) { + m_todo_terms.push_back(std::make_pair(i.first, i.second * coeff)); + } + } + else { + result += m_variable_values[wi] * coeff; } - m_variable_values[vi] = result; - return result; } - UNREACHABLE(); - return m_variable_values[vi]; + m_variable_values[vi] = result; + return result; } void init_variable_values() { @@ -1230,9 +1257,12 @@ namespace smt { // create a bound atom representing term <= k app_ref mk_bound(lp::lar_term const& term, rational const& k) { - SASSERT(k.is_int()); - app_ref t = mk_term(term, true); - app_ref atom(a.mk_le(t, a.mk_numeral(k, true)), m); + app_ref t = mk_term(term, k.is_int()); + app_ref atom(a.mk_le(t, a.mk_numeral(k, k.is_int())), m); + expr_ref atom1(m); + proof_ref atomp(m); + ctx().get_simplifier()(atom, atom1, atomp); + atom = to_app(atom1); TRACE("arith", tout << atom << "\n"; m_solver->print_term(term, tout << "bound atom: "); tout << " <= " << k << "\n"; display(tout); @@ -1243,7 +1273,10 @@ namespace smt { } lbool check_lia() { - if (m.canceled()) return l_undef; + if (m.canceled()) { + TRACE("arith", tout << "canceled\n";); + return l_undef; + } lp::lar_term term; lp::mpq k; lp::explanation ex; // TBD, this should be streamlined accross different explanations @@ -1288,7 +1321,10 @@ namespace smt { lbool check_nra() { m_use_nra_model = false; - if (m.canceled()) return l_undef; + if (m.canceled()) { + TRACE("arith", tout << "canceled\n";); + return l_undef; + } if (!m_nra) return l_true; if (!m_nra->need_check()) return l_true; m_a1 = 0; m_a2 = 0; @@ -1476,6 +1512,7 @@ namespace smt { void consume(rational const& v, unsigned j) override { m_imp.set_evidence(j); + m_imp.m_explanation.push_back(std::make_pair(v, j)); } }; @@ -1512,6 +1549,7 @@ namespace smt { if (lit == null_literal) { continue; } + TRACE("arith", tout << lit << " bound: " << *b << " first: " << first << "\n";); m_solver->settings().st().m_num_of_implied_bounds ++; if (first) { @@ -1519,6 +1557,7 @@ namespace smt { m_core.reset(); m_eqs.reset(); m_params.reset(); + m_explanation.clear(); local_bound_propagator bp(*this); m_solver->explain_implied_bound(be, bp); } @@ -1529,7 +1568,7 @@ namespace smt { tout << "\n --> "; ctx().display_literal_verbose(tout, lit); tout << "\n"; - // display_evidence(tout, m_explanation); + display_evidence(tout, m_explanation); m_solver->print_implied_bound(be, tout); ); DEBUG_CODE( @@ -1547,8 +1586,8 @@ namespace smt { SASSERT(validate_assign(lit)); if (m_core.size() < small_lemma_size() && m_eqs.empty()) { m_core2.reset(); - for (unsigned i = 0; i < m_core.size(); ++i) { - m_core2.push_back(~m_core[i]); + for (auto const& c : m_core) { + m_core2.push_back(~c); } m_core2.push_back(lit); justification * js = nullptr; @@ -1569,27 +1608,27 @@ namespace smt { literal is_bound_implied(lp::lconstraint_kind k, rational const& value, lp_api::bound const& b) const { if ((k == lp::LE || k == lp::LT) && b.get_bound_kind() == lp_api::upper_t && value <= b.get_value()) { - // v <= value <= b.get_value() => v <= b.get_value() + TRACE("arith", tout << "v <= value <= b.get_value() => v <= b.get_value() \n";); return literal(b.get_bv(), false); } if ((k == lp::GE || k == lp::GT) && b.get_bound_kind() == lp_api::lower_t && b.get_value() <= value) { - // b.get_value() <= value <= v => b.get_value() <= v + TRACE("arith", tout << "b.get_value() <= value <= v => b.get_value() <= v \n";); return literal(b.get_bv(), false); } if (k == lp::LE && b.get_bound_kind() == lp_api::lower_t && value < b.get_value()) { - // v <= value < b.get_value() => v < b.get_value() + TRACE("arith", tout << "v <= value < b.get_value() => v < b.get_value()\n";); return literal(b.get_bv(), true); } if (k == lp::LT && b.get_bound_kind() == lp_api::lower_t && value <= b.get_value()) { - // v < value <= b.get_value() => v < b.get_value() + TRACE("arith", tout << "v < value <= b.get_value() => v < b.get_value()\n";); return literal(b.get_bv(), true); } if (k == lp::GE && b.get_bound_kind() == lp_api::upper_t && b.get_value() < value) { - // b.get_value() < value <= v => b.get_value() < v + TRACE("arith", tout << "b.get_value() < value <= v => b.get_value() < v\n";); return literal(b.get_bv(), true); } if (k == lp::GT && b.get_bound_kind() == lp_api::upper_t && b.get_value() <= value) { - // b.get_value() <= value < v => b.get_value() < v + TRACE("arith", tout << "b.get_value() <= value < v => b.get_value() < v\n";); return literal(b.get_bv(), true); } @@ -1923,16 +1962,29 @@ namespace smt { ++m_stats.m_bounds_propagations; } + svector m_todo_vars; + void add_use_lists(lp_api::bound* b) { theory_var v = b->get_var(); lp::var_index vi = get_var_index(v); - if (m_solver->is_term(vi)) { + if (!m_solver->is_term(vi)) { + return; + } + m_todo_vars.push_back(vi); + while (!m_todo_vars.empty()) { + vi = m_todo_vars.back(); + m_todo_vars.pop_back(); lp::lar_term const& term = m_solver->get_term(vi); for (auto i = term.m_coeffs.begin(); i != term.m_coeffs.end(); ++i) { lp::var_index wi = i->first; - unsigned w = m_var_index2theory_var[wi]; - m_use_list.reserve(w + 1, ptr_vector()); - m_use_list[w].push_back(b); + if (m_solver->is_term(wi)) { + m_todo_vars.push_back(wi); + } + else { + unsigned w = m_var_index2theory_var[wi]; + m_use_list.reserve(w + 1, ptr_vector()); + m_use_list[w].push_back(b); + } } } } @@ -1940,13 +1992,24 @@ namespace smt { void del_use_lists(lp_api::bound* b) { theory_var v = b->get_var(); lp::var_index vi = m_theory_var2var_index[v]; - if (m_solver->is_term(vi)) { + if (!m_solver->is_term(vi)) { + return; + } + m_todo_vars.push_back(vi); + while (!m_todo_vars.empty()) { + vi = m_todo_vars.back(); + m_todo_vars.pop_back(); lp::lar_term const& term = m_solver->get_term(vi); for (auto i = term.m_coeffs.begin(); i != term.m_coeffs.end(); ++i) { lp::var_index wi = i->first; - unsigned w = m_var_index2theory_var[wi]; - SASSERT(m_use_list[w].back() == b); - m_use_list[w].pop_back(); + if (m_solver->is_term(wi)) { + m_todo_vars.push_back(wi); + } + else { + unsigned w = m_var_index2theory_var[wi]; + SASSERT(m_use_list[w].back() == b); + m_use_list[w].pop_back(); + } } } } @@ -2033,6 +2096,9 @@ namespace smt { lp::constraint_index ci; rational value; bool is_strict; + if (m_solver->is_term(wi)) { + return false; + } if (coeff.second.is_neg() == is_lub) { // -3*x ... <= lub based on lower bound for x. if (!m_solver->has_lower_bound(wi, ci, value, is_strict)) { @@ -2378,15 +2444,31 @@ namespace smt { SASSERT(m_use_nra_model); lp::var_index vi = m_theory_var2var_index[v]; if (m_solver->is_term(vi)) { - lp::lar_term const& term = m_solver->get_term(vi); - scoped_anum r1(m_nra->am()); - m_nra->am().set(r, term.m_v.to_mpq()); - - for (auto const coeff : term.m_coeffs) { - lp::var_index wi = coeff.first; - m_nra->am().set(r1, coeff.second.to_mpq()); - m_nra->am().mul(m_nra->value(wi), r1, r1); - m_nra->am().add(r1, r, r); + + m_todo_terms.push_back(std::make_pair(vi, rational::one())); + + m_nra->am().set(r, 0); + while (!m_todo_terms.empty()) { + rational wcoeff = m_todo_terms.back().second; + lp::var_index wi = m_todo_terms.back().first; // todo : got a warning "wi is not used" + m_todo_terms.pop_back(); + lp::lar_term const& term = m_solver->get_term(vi); + scoped_anum r1(m_nra->am()); + rational c1 = term.m_v * wcoeff; + m_nra->am().set(r1, c1.to_mpq()); + m_nra->am().add(r, r1, r); + for (auto const coeff : term.m_coeffs) { + lp::var_index wi = coeff.first; + c1 = coeff.second * wcoeff; + if (m_solver->is_term(wi)) { + m_todo_terms.push_back(std::make_pair(wi, c1)); + } + else { + m_nra->am().set(r1, c1.to_mpq()); + m_nra->am().mul(m_nra->value(wi), r1, r1); + m_nra->am().add(r1, r, r); + } + } } return r; } @@ -2431,11 +2513,22 @@ namespace smt { // Auxiliary verification utilities. + struct scoped_arith_mode { + smt_params& p; + scoped_arith_mode(smt_params& p) : p(p) { + p.m_arith_mode = AS_ARITH; + } + ~scoped_arith_mode() { + p.m_arith_mode = AS_LRA; + } + }; + bool validate_conflict() { 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_LRA) return true; + scoped_arith_mode _sa(ctx().get_fparams()); context nctx(m, ctx().get_fparams(), ctx().get_params()); add_background(nctx); bool result = l_true != nctx.check(); @@ -2448,7 +2541,8 @@ namespace smt { 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_LRA) return true; + scoped_arith_mode _sa(ctx().get_fparams()); context nctx(m, ctx().get_fparams(), ctx().get_params()); m_core.push_back(~lit); add_background(nctx); @@ -2486,10 +2580,14 @@ namespace smt { theory_lra::inf_eps maximize(theory_var v, expr_ref& blocker, bool& has_shared) { lp::var_index vi = m_theory_var2var_index.get(v, UINT_MAX); vector > coeffs; - rational coeff; + rational coeff(0); + // + // TBD: change API for maximize_term to take a proper term as input. + // if (m_solver->is_term(vi)) { const lp::lar_term& term = m_solver->get_term(vi); for (auto & ti : term.m_coeffs) { + SASSERT(!m_solver->is_term(ti.first)); coeffs.push_back(std::make_pair(ti.second, ti.first)); } coeff = term.m_v; @@ -2547,7 +2645,13 @@ namespace smt { app_ref mk_term(lp::lar_term const& term, bool is_int) { expr_ref_vector args(m); for (auto & ti : term.m_coeffs) { - theory_var w = m_var_index2theory_var[ti.first]; + theory_var w; + if (m_solver->is_term(ti.first)) { + w = m_term_index2theory_var[m_solver->adjust_term_index(ti.first)]; + } + else { + w = m_var_index2theory_var[ti.first]; + } expr* o = get_enode(w)->get_owner(); if (ti.second.is_one()) { args.push_back(o); diff --git a/src/smt/theory_str.cpp b/src/smt/theory_str.cpp index f2432b4fb..4e9b85367 100644 --- a/src/smt/theory_str.cpp +++ b/src/smt/theory_str.cpp @@ -11074,8 +11074,8 @@ namespace smt { ctx.force_phase(lit); } - zstring aStr = gen_val_string(len, options[i - l]); - expr * strAst; + zstring aStr = gen_val_string(len, options[static_cast(i) - static_cast(l)]); + expr * strAst; if (m_params.m_UseFastValueTesterCache) { if (!valueTesterCache.find(aStr, strAst)) { strAst = mk_string(aStr); diff --git a/src/test/lp/lp.cpp b/src/test/lp/lp.cpp index adcf9b969..911c127a1 100644 --- a/src/test/lp/lp.cpp +++ b/src/test/lp/lp.cpp @@ -1124,7 +1124,7 @@ void update_settings(argument_parser & args_parser, lp_settings& settings) { settings.harris_feasibility_tolerance = d; } if (get_int_from_args_parser("--random_seed", args_parser, n)) { - settings.random_seed(n); + settings.set_random_seed(n); } if (get_int_from_args_parser("--simplex_strategy", args_parser, n)) { settings.simplex_strategy() = static_cast(n); diff --git a/src/util/lp/bound_propagator.cpp b/src/util/lp/bound_propagator.cpp index cfb3aa68f..aac57af48 100644 --- a/src/util/lp/bound_propagator.cpp +++ b/src/util/lp/bound_propagator.cpp @@ -15,8 +15,13 @@ const impq & bound_propagator::get_low_bound(unsigned j) const { const impq & bound_propagator::get_upper_bound(unsigned j) const { return m_lar_solver.m_mpq_lar_core_solver.m_r_upper_bounds()[j]; } -void bound_propagator::try_add_bound(const mpq & v, unsigned j, bool is_low, bool coeff_before_j_is_pos, unsigned row_or_term_index, bool strict) { - j = m_lar_solver.adjust_column_index_to_term_index(j); +void bound_propagator::try_add_bound(mpq v, unsigned j, bool is_low, bool coeff_before_j_is_pos, unsigned row_or_term_index, bool strict) { + j = m_lar_solver.adjust_column_index_to_term_index(j); + if (m_lar_solver.is_term(j)) { + // lp treats terms as not having a free coefficient, restoring it below for the outside consumption + v += m_lar_solver.get_term(j).m_v; + } + lconstraint_kind kind = is_low? GE : LE; if (strict) kind = static_cast(kind / 2); @@ -27,20 +32,26 @@ void bound_propagator::try_add_bound(const mpq & v, unsigned j, bool is_low, boo if (is_low) { if (try_get_val(m_improved_low_bounds, j, k)) { auto & found_bound = m_ibounds[k]; - if (v > found_bound.m_bound || (v == found_bound.m_bound && found_bound.m_strict == false && strict)) + if (v > found_bound.m_bound || (v == found_bound.m_bound && found_bound.m_strict == false && strict)) { found_bound = implied_bound(v, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict); + TRACE("try_add_bound", m_lar_solver.print_implied_bound(found_bound, tout);); + } } else { m_improved_low_bounds[j] = m_ibounds.size(); m_ibounds.push_back(implied_bound(v, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict)); + TRACE("try_add_bound", m_lar_solver.print_implied_bound(m_ibounds.back(), tout);); } } else { // the upper bound case if (try_get_val(m_improved_upper_bounds, j, k)) { auto & found_bound = m_ibounds[k]; - if (v < found_bound.m_bound || (v == found_bound.m_bound && found_bound.m_strict == false && strict)) + if (v < found_bound.m_bound || (v == found_bound.m_bound && found_bound.m_strict == false && strict)) { found_bound = implied_bound(v, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict); + TRACE("try_add_bound", m_lar_solver.print_implied_bound(found_bound, tout);); + } } else { m_improved_upper_bounds[j] = m_ibounds.size(); m_ibounds.push_back(implied_bound(v, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict)); + TRACE("try_add_bound", m_lar_solver.print_implied_bound(m_ibounds.back(), tout);); } } } diff --git a/src/util/lp/column_namer.h b/src/util/lp/column_namer.h index 97d371f48..deb393796 100644 --- a/src/util/lp/column_namer.h +++ b/src/util/lp/column_namer.h @@ -40,7 +40,7 @@ public: T a; unsigned i; while (it->next(a, i)) { - coeff.emplace_back(a, i); + coeff.push_back(std::make_pair(a, i)); } print_linear_combination_of_column_indices_only(coeff, out); } @@ -65,7 +65,7 @@ public: else if (val != numeric_traits::one()) out << T_to_string(val); - out << "_" << it.second; + out << "v" << it.second; } } diff --git a/src/util/lp/implied_bound.h b/src/util/lp/implied_bound.h index f1c711ffa..e28e9cf74 100644 --- a/src/util/lp/implied_bound.h +++ b/src/util/lp/implied_bound.h @@ -28,7 +28,6 @@ struct implied_bound { bool m_coeff_before_j_is_pos; unsigned m_row_or_term_index; bool m_strict; - lconstraint_kind kind() const { lconstraint_kind k = m_is_low_bound? GE : LE; if (m_strict) @@ -52,6 +51,7 @@ struct implied_bound { m_is_low_bound(low_bound), m_coeff_before_j_is_pos(coeff_before_j_is_pos), m_row_or_term_index(row_or_term_index), - m_strict(strict) {} + m_strict(strict) { + } }; } diff --git a/src/util/lp/int_solver.cpp b/src/util/lp/int_solver.cpp index 3345a6389..d717c4047 100644 --- a/src/util/lp/int_solver.cpp +++ b/src/util/lp/int_solver.cpp @@ -8,7 +8,6 @@ namespace lp { void int_solver::fix_non_base_columns() { - lp_assert(is_feasible() && inf_int_set_is_correct()); auto & lcs = m_lar_solver->m_mpq_lar_core_solver; bool change = false; for (unsigned j : lcs.m_r_nbasis) { @@ -66,6 +65,10 @@ 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(); +} + int int_solver::find_inf_int_base_column() { if (inf_int_set().is_empty()) return -1; @@ -81,11 +84,12 @@ int int_solver::find_inf_int_boxed_base_column_with_smallest_range() { mpq range; mpq new_range; mpq small_range_thresold(1024); - unsigned n = 0; + 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)); if (!is_boxed(j)) continue; new_range = lcs.m_r_upper_bounds()[j].x - lcs.m_r_low_bounds()[j].x; @@ -104,8 +108,8 @@ int int_solver::find_inf_int_boxed_base_column_with_smallest_range() { continue; } if (new_range == range) { - n++; - if (settings().random_next() % n == 0) { + lp_assert(n >= 1); + if (settings().random_next() % (++n) == 0) { result = j; continue; } @@ -115,32 +119,30 @@ int int_solver::find_inf_int_boxed_base_column_with_smallest_range() { } -bool int_solver::is_gomory_cut_target() { - m_iter_on_gomory_row->reset(); +bool int_solver::is_gomory_cut_target(linear_combination_iterator &iter) { unsigned j; - TRACE("gomory_cut", m_lar_solver->print_linear_iterator(m_iter_on_gomory_row, tout); - m_iter_on_gomory_row->reset(); - ); - - while (m_iter_on_gomory_row->next(j)) { - // All non base variables must be at their bounds and assigned to rationals (that is, infinitesimals are not allowed). - if (j != m_gomory_cut_inf_column && (!at_bound(j) || !is_zero(get_value(j).y))) { + lp_assert(iter.is_reset()); + // All non base variables must be at their bounds and assigned to rationals (that is, infinitesimals are not allowed). + while (iter.next(j)) { + if (is_base(j)) continue; + if (!is_zero(get_value(j).y)) { TRACE("gomory_cut", tout << "row is not gomory cut target:\n"; display_column(tout, j); - tout << "at_bound: " << at_bound(j) << "\n"; tout << "infinitesimal: " << !is_zero(get_value(j).y) << "\n";); + iter.reset(); return false; } } - m_iter_on_gomory_row->reset(); + iter.reset(); return true; } -void int_solver::real_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& pol, explanation & expl) { - mpq f_0 = fractional_part(get_value(m_gomory_cut_inf_column)); +void int_solver::real_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& pol, explanation & expl, unsigned gomory_cut_inf_column) { + TRACE("gomory_cut_detail_real", tout << "real\n";); + mpq f_0 = fractional_part(get_value(gomory_cut_inf_column)); mpq new_a; - if (at_lower(x_j)) { + if (at_low(x_j)) { if (a.is_pos()) { new_a = a / (1 - f_0); } @@ -148,8 +150,8 @@ void int_solver::real_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, l new_a = a / f_0; new_a.neg(); } - k += lower_bound(x_j).x * k; // k.addmul(new_a, lower_bound(x_j).x); // is it a faster operation - + k.addmul(new_a, low_bound(x_j).x); // is it a faster operation than + // k += lower_bound(x_j).x * new_a; expl.push_justification(column_low_bound_constraint(x_j), new_a); } else { @@ -161,11 +163,11 @@ void int_solver::real_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, l else { new_a = a / (mpq(1) - f_0); } - k += upper_bound(x_j).x * k; // k.addmul(new_a, upper_bound(x_j).get_rational()); + k.addmul(new_a, upper_bound(x_j).x); // k += upper_bound(x_j).x * new_a; expl.push_justification(column_upper_bound_constraint(x_j), new_a); } - TRACE("gomory_cut_detail", tout << a << "*v" << x_j << " k: " << k << "\n";); - pol.add_monoid(new_a, x_j); + TRACE("gomory_cut_detail_real", tout << a << "*v" << x_j << " k: " << k << "\n";); + pol.add_monomial(new_a, x_j); } constraint_index int_solver::column_upper_bound_constraint(unsigned j) const { @@ -177,152 +179,203 @@ constraint_index int_solver::column_low_bound_constraint(unsigned j) const { } -void int_solver::int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term & pol, explanation& expl, mpq & lcm_den) { - mpq f_0 = fractional_part(get_value(m_gomory_cut_inf_column)); +void int_solver::int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term & t, explanation& expl, mpq & lcm_den, unsigned inf_column) { lp_assert(is_int(x_j)); + lp_assert(!a.is_int()); + mpq f_0 = fractional_part(get_value(inf_column)); + lp_assert(f_0 > zero_of_type() && f_0 < one_of_type()); mpq f_j = fractional_part(a); TRACE("gomory_cut_detail", - tout << a << "*v" << x_j << "\n"; - tout << "fractional_part: " << fractional_part(a) << "\n"; + tout << a << " x_j" << x_j << " k = " << k << "\n"; tout << "f_j: " << f_j << "\n"; tout << "f_0: " << f_0 << "\n"; - tout << "one_minus_f_0: " << 1 - f_0 << "\n";); - if (!f_j.is_zero()) { - mpq new_a; - if (at_lower(x_j)) { - auto one_minus_f_0 = 1 - f_0; - if (f_j <= one_minus_f_0) { - new_a = f_j / one_minus_f_0; - } - else { - new_a = (1 - f_j) / f_0; - } - k.addmul(new_a, lower_bound(x_j).x); - expl.push_justification(column_low_bound_constraint(x_j), new_a); + tout << "1 - f_0: " << 1 - f_0 << "\n"; + tout << "at_low(" << x_j << ") = " << at_low(x_j) << std::endl; + ); + lp_assert (!f_j.is_zero()); + mpq new_a; + if (at_low(x_j)) { + auto one_minus_f_0 = 1 - f_0; + if (f_j <= one_minus_f_0) { + new_a = f_j / one_minus_f_0; } else { - SASSERT(at_upper(x_j)); - if (f_j <= f_0) { - new_a = f_j / f_0; - } - else { - new_a = (mpq(1) - f_j) / 1 - f_0; - } - new_a.neg(); // the upper terms are inverted - k.addmul(new_a, upper_bound(x_j).x); - expl.push_justification(column_upper_bound_constraint(x_j), new_a); + new_a = (1 - f_j) / f_0; } - TRACE("gomory_cut_detail", tout << "new_a: " << new_a << " k: " << k << "\n";); - pol.add_monoid(new_a, x_j); - lcm_den = lcm(lcm_den, denominator(new_a)); + k.addmul(new_a, low_bound(x_j).x); + expl.push_justification(column_low_bound_constraint(x_j), new_a); } + else { + lp_assert(at_upper(x_j)); + if (f_j <= f_0) { + new_a = f_j / f_0; + } + else { + new_a = (mpq(1) - f_j) / (1 - f_0); + } + new_a.neg(); // the upper terms are inverted + k.addmul(new_a, upper_bound(x_j).x); + expl.push_justification(column_upper_bound_constraint(x_j), new_a); + } + TRACE("gomory_cut_detail", tout << "new_a: " << new_a << " k: " << k << "\n";); + t.add_monomial(new_a, x_j); + lcm_den = lcm(lcm_den, denominator(new_a)); } lia_move int_solver::report_conflict_from_gomory_cut(mpq & k) { - TRACE("empty_pol", - display_row_info(tout, - m_lar_solver->m_mpq_lar_core_solver.m_r_heading[m_gomory_cut_inf_column]);); + TRACE("empty_pol",); lp_assert(k.is_pos()); // conflict 0 >= k where k is positive k.neg(); // returning 0 <= -k return lia_move::conflict; } -void int_solver::gomory_cut_adjust_t_and_k_for_size_gt_1( - vector> & pol, - lar_term & t, - mpq &k, - unsigned num_ints, - mpq & lcm_den) { - if (num_ints > 0) { - lcm_den = lcm(lcm_den, denominator(k)); - TRACE("gomory_cut_detail", tout << "k: " << k << " lcm_den: " << lcm_den << "\n"; - linear_combination_iterator_on_vector pi(pol); - m_lar_solver->print_linear_iterator(&pi, tout); - tout << "\nk: " << k << "\n";); - lp_assert(lcm_den.is_pos()); - if (!lcm_den.is_one()) { - // normalize coefficients of integer parameters to be integers. - for (auto & pi: pol) { - pi.first *= lcm_den; - SASSERT(!is_int(pi.second) || pi.first.is_int()); - } - k *= lcm_den; - } - TRACE("gomory_cut_detail", tout << "after *lcm_den\n"; - for (unsigned i = 0; i < pol.size(); i++) { - tout << pol[i].first << " * v" << pol[i].second << "\n"; - } - tout << "k: " << k << "\n";); +void int_solver::gomory_cut_adjust_t_and_k(vector> & pol, + lar_term & t, + mpq &k, + bool some_ints, + mpq & lcm_den) { + if (!some_ints) + return; + + t.clear(); + if (pol.size() == 1) { + unsigned v = pol[0].second; + lp_assert(is_int(v)); + bool k_is_int = k.is_int(); + const mpq& a = pol[0].first; + k /= a; + if (a.is_pos()) { // we have av >= k + if (!k_is_int) + k = ceil(k); + // switch size + t.add_monomial(- mpq(1), v); + k.neg(); + } else { + if (!k_is_int) + k = floor(k); + t.add_monomial(mpq(1), v); + } + } else if (some_ints) { + lcm_den = lcm(lcm_den, denominator(k)); + lp_assert(lcm_den.is_pos()); + if (!lcm_den.is_one()) { + // normalize coefficients of integer parameters to be integers. + for (auto & pi: pol) { + pi.first *= lcm_den; + SASSERT(!is_int(pi.second) || pi.first.is_int()); + } + k *= lcm_den; } - t.clear(); // negate everything to return -pol <= -k for (const auto & pi: pol) - t.add_monoid(-pi.first, pi.second); + t.add_monomial(-pi.first, pi.second); k.neg(); -} - - -void int_solver::gomory_cut_adjust_t_and_k_for_size_1(const vector> & pol, lar_term& t, mpq &k) { - lp_assert(pol.size() == 1); - unsigned j = pol[0].second; - k /= pol[0].first; - bool is_lower = pol[0].first.is_pos(); - if (is_int(j) && !k.is_int()) { - k = is_lower?ceil(k):floor(k); - } - if (is_lower) { // returning -t <= -k which is equivalent to t >= k - k.neg(); - t.negate(); } } +bool int_solver::current_solution_is_inf_on_cut(const lar_term& t, const mpq& k) const { + const auto & x = m_lar_solver->m_mpq_lar_core_solver.m_r_x; + impq v = t.apply(x); + TRACE( + "current_solution_is_inf_on_cut", tout << "v = " << v << " k = " << k << std::endl; + if (v <=k) { + tout << "v <= k - it should not happen!\n"; + } + ); + + return v > k; +} - -lia_move int_solver::report_gomory_cut(lar_term& t, mpq& k, mpq &lcm_den, unsigned num_ints) { +void int_solver::adjust_term_and_k_for_some_ints_case_gomory(lar_term& t, mpq& k, mpq &lcm_den) { lp_assert(!t.is_empty()); auto pol = t.coeffs_as_vector(); - if (pol.size() == 1) - gomory_cut_adjust_t_and_k_for_size_1(pol, t, k); - else - gomory_cut_adjust_t_and_k_for_size_gt_1(pol, t, k, num_ints, lcm_den); - m_lar_solver->subs_terms_for_debugging(t); // todo: remove later - return lia_move::cut; + t.clear(); + if (pol.size() == 1) { + TRACE("gomory_cut_detail", tout << "pol.size() is 1" << std::endl;); + unsigned v = pol[0].second; + lp_assert(is_int(v)); + const mpq& a = pol[0].first; + k /= a; + if (a.is_pos()) { // we have av >= k + if (!k.is_int()) + k = ceil(k); + // switch size + t.add_monomial(- mpq(1), v); + k.neg(); + } else { + if (!k.is_int()) + k = floor(k); + t.add_monomial(mpq(1), v); + } + } else { + TRACE("gomory_cut_detail", tout << "pol.size() > 1" << std::endl;); + lcm_den = lcm(lcm_den, denominator(k)); + lp_assert(lcm_den.is_pos()); + if (!lcm_den.is_one()) { + // normalize coefficients of integer parameters to be integers. + for (auto & pi: pol) { + pi.first *= lcm_den; + SASSERT(!is_int(pi.second) || pi.first.is_int()); + } + k *= lcm_den; + } + // negate everything to return -pol <= -k + for (const auto & pi: pol) + t.add_monomial(-pi.first, pi.second); + k.neg(); + } + TRACE("gomory_cut_detail", tout << "k = " << k << std::endl;); + lp_assert(k.is_int()); } -lia_move int_solver::mk_gomory_cut(lar_term& t, mpq& k, explanation & expl) { +lia_move int_solver::mk_gomory_cut(lar_term& t, mpq& k, explanation & expl, unsigned inf_col, linear_combination_iterator& iter) { - lp_assert(column_is_int_inf(m_gomory_cut_inf_column)); + lp_assert(column_is_int_inf(inf_col)); - TRACE("gomory_cut", tout << "applying cut at:\n"; m_lar_solver->print_linear_iterator(m_iter_on_gomory_row, tout); tout << std::endl; m_iter_on_gomory_row->reset();); + TRACE("gomory_cut", + tout << "applying cut at:\n"; m_lar_solver->print_linear_iterator_indices_only(&iter, tout); tout << std::endl; + iter.reset(); + unsigned j; + while(iter.next(j)) { + m_lar_solver->m_mpq_lar_core_solver.m_r_solver.print_column_info(j, tout); + } + iter.reset(); + tout << "inf_col = " << inf_col << std::endl; + ); // gomory will be t >= k k = 1; mpq lcm_den(1); - unsigned num_ints = 0; unsigned x_j; mpq a; - while (m_iter_on_gomory_row->next(a, x_j)) { - if (x_j == m_gomory_cut_inf_column) + bool some_int_columns = false; + lp_assert(iter.is_reset()); + while (iter.next(a, x_j)) { + if (x_j == inf_col) continue; // make the format compatible with the format used in: Integrating Simplex with DPLL(T) a.neg(); if (is_real(x_j)) - real_case_in_gomory_cut(a, x_j, k, t, expl); + real_case_in_gomory_cut(a, x_j, k, t, expl, inf_col); else { - num_ints++; - int_case_in_gomory_cut(a, x_j, k, t, expl, lcm_den); + if (a.is_int()) continue; // f_j will be zero and no monomial will be added + some_int_columns = true; + int_case_in_gomory_cut(a, x_j, k, t, expl, lcm_den, inf_col); } } if (t.is_empty()) return report_conflict_from_gomory_cut(k); + if (some_int_columns) + adjust_term_and_k_for_some_ints_case_gomory(t, k, lcm_den); - return report_gomory_cut(t, k, lcm_den, num_ints); - + lp_assert(current_solution_is_inf_on_cut(t, k)); + m_lar_solver->subs_term_columns(t); + return lia_move::cut; } void int_solver::init_check_data() { @@ -331,61 +384,53 @@ void int_solver::init_check_data() { m_old_values_data.resize(n); } -int int_solver::find_next_free_var_in_gomory_row() { - lp_assert(m_iter_on_gomory_row != nullptr); +int int_solver::find_free_var_in_gomory_row(linear_combination_iterator& iter) { unsigned j; - while(m_iter_on_gomory_row->next(j)) { - if (j != m_gomory_cut_inf_column && is_free(j)) + while(iter.next(j)) { + if (!is_base(j) && is_free(j)) return static_cast(j); } + iter.reset(); return -1; } -lia_move int_solver::proceed_with_gomory_cut(lar_term& t, mpq& k, explanation& ex) { - int j = find_next_free_var_in_gomory_row(); - if (j != -1) { - m_found_free_var_in_gomory_row = true; +lia_move int_solver::proceed_with_gomory_cut(lar_term& t, mpq& k, explanation& ex, unsigned j, linear_combination_iterator& iter) { + int free_j = find_free_var_in_gomory_row(iter); + if (free_j != -1) { lp_assert(t.is_empty()); - t.add_monoid(mpq(1), j); + t.add_monomial(mpq(1), m_lar_solver->adjust_column_index_to_term_index(free_j)); k = zero_of_type(); return lia_move::branch; // branch on a free column } - if (m_found_free_var_in_gomory_row || !is_gomory_cut_target()) { - m_found_free_var_in_gomory_row = false; - delete m_iter_on_gomory_row; - m_iter_on_gomory_row = nullptr; - return lia_move::continue_with_check; - } - lia_move ret = mk_gomory_cut(t, k, ex); - delete m_iter_on_gomory_row; - m_iter_on_gomory_row = nullptr; - return ret; + if (!is_gomory_cut_target(iter)) + return create_branch_on_column(j, t, k); + + return mk_gomory_cut(t, k, ex, j, iter); } -lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex) { - if (m_iter_on_gomory_row != nullptr) { - auto ret = proceed_with_gomory_cut(t, k, ex); - TRACE("gomory_cut", tout << "term t = "; m_lar_solver->print_term_as_indices(t, tout);); - if (ret != lia_move::continue_with_check) - return ret; - } +unsigned int_solver::row_of_basic_column(unsigned j) const { + return m_lar_solver->m_mpq_lar_core_solver.m_r_heading[j]; +} +lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex) { init_check_data(); lp_assert(inf_int_set_is_correct()); - // currently it is a reimplementation of + // it is mostly a reimplementation of // final_check_status theory_arith::check_int_feasibility() // from theory_arith_int.h - if (m_lar_solver->model_is_int_feasible()) + if (!has_inf_int()) return lia_move::ok; - if (!gcd_test(ex)) - return lia_move::conflict; + if (settings().m_run_gcd_test) + if (!gcd_test(ex)) + return lia_move::conflict; /* if (m_params.m_arith_euclidean_solver) apply_euclidean_solver(); - */ + bool track_pivoted_rows = m_lar_solver->get_track_pivoted_rows(); + m_lar_solver->set_track_pivoted_rows(false); m_lar_solver->pivot_fixed_vars_from_basis(); lean_assert(m_lar_solver->m_mpq_lar_core_solver.r_basis_is_OK()); patch_int_infeasible_columns(); @@ -395,93 +440,90 @@ lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex) { lean_assert(is_feasible()); TRACE("arith_int_rows", trace_inf_rows();); - if (find_inf_int_base_column() == -1) + if (!has_inf_int()) { + m_lar_solver->set_track_pivoted_rows(track_pivoted_rows); return lia_move::ok; - - - if ((++m_branch_cut_counter) % settings().m_int_branch_cut_threshold == 0) { - move_non_base_vars_to_bounds(); // todo track changed variables - lp_status st = m_lar_solver->find_feasible_solution(); - if (st != lp_status::FEASIBLE && st != lp_status::OPTIMAL) { - return lia_move::give_up; - } - lp_assert(inf_int_set_is_correct()); - // init_inf_int_set(); // todo - can we avoid this call? - int j = find_inf_int_base_column(); - if (j != -1) { - TRACE("arith_int", tout << "j = " << j << " does not have an integer assignment: " << get_value(j) << "\n";); - unsigned row_index = m_lar_solver->m_mpq_lar_core_solver.m_r_heading[j]; - if (!mk_gomory_cut(row_index, ex)) { + } + TRACE("gomory_cut", tout << m_branch_cut_counter+1 << ", " << settings().m_int_branch_cut_gomory_threshold << std::endl;); + if ((++m_branch_cut_counter) % settings().m_int_branch_cut_gomory_threshold == 0) { + if (move_non_base_vars_to_bounds()) { + lp_status st = m_lar_solver->find_feasible_solution(); + lp_assert(non_basic_columns_are_at_bounds()); + if (st != lp_status::FEASIBLE && st != lp_status::OPTIMAL) { + TRACE("arith_int", tout << "give_up\n";); + m_lar_solver->set_track_pivoted_rows(track_pivoted_rows); return lia_move::give_up; - // silent failure } - return lia_move::cut; } - } - else { int j = find_inf_int_base_column(); - if (j != -1) { - TRACE("arith_int", tout << "j" << j << " does not have an integer assignment: " << get_value(j) << "\n";); - - lp_assert(t.is_empty()); - t.add_monoid(1, j); - k = floor(get_value(j)); - TRACE("arith_int", tout << "branching v" << j << " = " << get_value(j) << "\n"; - display_column(tout, j); - tout << "k = " << k << std::endl; - ); - return lia_move::branch; - } + lp_assert(j != -1); + TRACE("arith_int", tout << "j = " << j << " does not have an integer assignment: " << get_value(j) << "\n";); + auto iter_on_gomory_row = m_lar_solver->get_iterator_on_row(row_of_basic_column(j)); + lia_move ret = proceed_with_gomory_cut(t, k, ex, j, *iter_on_gomory_row); + delete iter_on_gomory_row; + m_lar_solver->set_track_pivoted_rows(track_pivoted_rows); + return ret; } - lp_assert(m_lar_solver->m_mpq_lar_core_solver.r_basis_is_OK()); - // return true; - return lia_move::give_up; + + m_lar_solver->set_track_pivoted_rows(track_pivoted_rows); + return create_branch_on_column(find_inf_int_base_column(), t, k); } -void int_solver::move_non_base_vars_to_bounds() { +bool int_solver::move_non_base_vars_to_bounds() { auto & lcs = m_lar_solver->m_mpq_lar_core_solver; + bool change = false; for (unsigned j : lcs.m_r_nbasis) { auto & val = lcs.m_r_x[j]; switch (lcs.m_column_types()[j]) { case column_type::boxed: - if (val != lcs.m_r_low_bounds()[j] && val != lcs.m_r_upper_bounds()[j]) - set_value(j, lcs.m_r_low_bounds()[j]); + if (val != lcs.m_r_low_bounds()[j] && val != lcs.m_r_upper_bounds()[j]) { + set_value_for_nbasic_column(j, lcs.m_r_low_bounds()[j]); + change = true; + } break; case column_type::low_bound: - if (val != lcs.m_r_low_bounds()[j]) - set_value(j, lcs.m_r_low_bounds()[j]); + if (val != lcs.m_r_low_bounds()[j]) { + set_value_for_nbasic_column(j, lcs.m_r_low_bounds()[j]); + change = true; + } break; case column_type::upper_bound: - if (val != lcs.m_r_upper_bounds()[j]) - set_value(j, lcs.m_r_upper_bounds()[j]); + if (val != lcs.m_r_upper_bounds()[j]) { + set_value_for_nbasic_column(j, lcs.m_r_upper_bounds()[j]); + change = true; + } break; default: - if (is_int(j) && !val.is_int()) { - set_value(j, impq(floor(val))); + if (is_int(j) && !val.is_int()) { + set_value_for_nbasic_column(j, impq(floor(val))); + change = true; } } } + return change; } +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]; + 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); +} void int_solver::set_value_for_nbasic_column(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]; - if (!m_old_values_set.contains(j)) { + if (m_lar_solver->has_int_var() && !m_old_values_set.contains(j)) { m_old_values_set.insert(j); m_old_values_data[j] = x; } auto delta = new_val - x; x = new_val; - m_lar_solver->change_basic_x_by_delta_on_column(j, delta); - - auto * it = get_column_iterator(j); update_column_in_int_inf_set(j); - unsigned i; - while (it->next(i)) - update_column_in_int_inf_set(m_lar_solver->m_mpq_lar_core_solver.m_r_basis[i]); - delete it; + m_lar_solver->change_basic_columns_dependend_on_a_given_nb_column(j, delta); } void int_solver::patch_int_infeasible_columns() { @@ -660,7 +702,7 @@ bool int_solver::ext_gcd_test(iterator_on_row & it, else { // l += ncoeff * upper_bound(j).get_rational(); l.addmul(ncoeff, m_lar_solver->column_upper_bound(j).x); - // u += ncoeff * lower_bound(j).get_rational(); + // u += ncoeff * low_bound(j).get_rational(); u.addmul(ncoeff, m_lar_solver->column_low_bound(j).x); } add_to_explanation_from_fixed_or_boxed_column(j, ex); @@ -702,10 +744,11 @@ int_solver::int_solver(lar_solver* lar_slv) : m_branch_cut_counter(0) { lp_assert(m_old_values_set.size() == 0); m_old_values_set.resize(lar_slv->A_r().column_count()); - m_old_values_data.resize(lar_slv->A_r().column_count(), zero_of_type()); + m_old_values_data.resize(lar_slv->A_r().column_count(), zero_of_type()); + m_lar_solver->set_int_solver(this); } -bool int_solver::lower(unsigned j) const { +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: case column_type::boxed: @@ -716,7 +759,7 @@ bool int_solver::lower(unsigned j) const { } } -bool int_solver::upper(unsigned j) const { +bool int_solver::has_upper(unsigned j) const { switch (m_lar_solver->m_mpq_lar_core_solver.m_column_types()[j]) { case column_type::fixed: case column_type::boxed: @@ -727,14 +770,6 @@ bool int_solver::upper(unsigned j) const { } } -const impq& int_solver::lower_bound(unsigned j) const { - return m_lar_solver->m_mpq_lar_core_solver.m_r_low_bounds()[j]; -} - -const impq& int_solver::upper_bound(unsigned j) const { - return m_lar_solver->m_mpq_lar_core_solver.m_r_upper_bounds()[j]; -} - void set_lower(impq & l, bool & inf_l, @@ -754,73 +789,62 @@ void set_upper(impq & u, } } -bool int_solver::get_freedom_interval_for_column(unsigned x_j, bool & inf_l, impq & l, bool & inf_u, impq & u, mpq & m) { +bool int_solver::get_freedom_interval_for_column(unsigned j, 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[x_j] >= 0) // the basic var + if (lcs.m_r_heading[j] >= 0) // the basic var return false; - impq const & x_j_val = lcs.m_r_x[x_j]; - linear_combination_iterator *it = get_column_iterator(x_j); + impq const & xj = get_value(j); + linear_combination_iterator *it = get_column_iterator(j); inf_l = true; inf_u = true; l = u = zero_of_type(); m = mpq(1); - if (lower(x_j)) { - set_lower(l, inf_l, lower_bound(x_j)); + if (has_low(j)) { + set_lower(l, inf_l, low_bound(j)); } - if (upper(x_j)) { - set_upper(u, inf_u, upper_bound(x_j)); + if (has_upper(j)) { + set_upper(u, inf_u, upper_bound(j)); } - mpq a_ij; unsigned i; - while (it->next(a_ij, i)) { - unsigned x_i = lcs.m_r_basis[i]; - impq const & x_i_val = lcs.m_r_x[x_i]; - if (is_int(x_i) && is_int(x_j) && !a_ij.is_int()) - m = lcm(m, denominator(a_ij)); - bool x_i_lower = lower(x_i); - bool x_i_upper = upper(x_i); - if (a_ij.is_neg()) { - if (x_i_lower) { - impq new_l = x_j_val + ((x_i_val - lcs.m_r_low_bounds()[x_i]) / a_ij); - set_lower(l, inf_l, new_l); - if (!inf_l && !inf_u && l == u) break;; - } - if (x_i_upper) { - impq new_u = x_j_val + ((x_i_val - lcs.m_r_upper_bounds()[x_i]) / a_ij); - set_upper(u, inf_u, new_u); - if (!inf_l && !inf_u && l == u) break;; - } + mpq a; // the coefficient in the column + unsigned row_index; + while (it->next(a, row_index)) { + unsigned i = lcs.m_r_basis[row_index]; + impq const & xi = get_value(i); + if (is_int(i) && is_int(j) && !a.is_int()) + m = lcm(m, denominator(a)); + if (a.is_neg()) { + if (has_low(i)) + set_lower(l, inf_l, xj + (xi - lcs.m_r_low_bounds()[i]) / a); + + if (has_upper(i)) + set_upper(u, inf_u, xj + (xi - lcs.m_r_upper_bounds()[i]) / a); } else { - if (x_i_upper) { - impq new_l = x_j_val + ((x_i_val - lcs.m_r_upper_bounds()[x_i]) / a_ij); - set_lower(l, inf_l, new_l); - if (!inf_l && !inf_u && l == u) break;; - } - if (x_i_lower) { - impq new_u = x_j_val + ((x_i_val - lcs.m_r_low_bounds()[x_i]) / a_ij); - set_upper(u, inf_u, new_u); - if (!inf_l && !inf_u && l == u) break;; - } + if (has_upper(i)) + set_lower(l, inf_l, xj + (xi - lcs.m_r_upper_bounds()[i]) / a); + if (has_low(i)) + set_upper(u, inf_u, xj + (xi - lcs.m_r_low_bounds()[i]) / a); } + if (!inf_l && !inf_u && l == u) break;; } delete it; TRACE("freedom_interval", tout << "freedom variable for:\n"; - tout << m_lar_solver->get_column_name(x_j); + tout << m_lar_solver->get_column_name(j); tout << "["; if (inf_l) tout << "-oo"; else tout << l; tout << "; "; if (inf_u) tout << "oo"; else tout << u; tout << "]\n"; - tout << "val = " << get_value(x_j) << "\n"; + tout << "val = " << get_value(j) << "\n"; ); - lp_assert(inf_l || l <= get_value(x_j)); - lp_assert(inf_u || u >= get_value(x_j)); + lp_assert(inf_l || l <= get_value(j)); + lp_assert(inf_u || u >= get_value(j)); return true; } @@ -834,7 +858,7 @@ bool int_solver::is_real(unsigned j) const { } bool int_solver::value_is_int(unsigned j) const { - return m_lar_solver->m_mpq_lar_core_solver.m_r_x[j].is_int(); + return m_lar_solver->column_value_is_int(j); } @@ -856,8 +880,10 @@ void int_solver::display_column(std::ostream & out, unsigned j) const { 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)))) + 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; } @@ -881,6 +907,10 @@ bool int_solver::is_boxed(unsigned j) const { return m_lar_solver->m_mpq_lar_core_solver.m_column_types[j] == column_type::boxed; } +bool int_solver::is_fixed(unsigned j) const { + return m_lar_solver->m_mpq_lar_core_solver.m_column_types[j] == column_type::fixed; +} + bool int_solver::is_free(unsigned j) const { return m_lar_solver->m_mpq_lar_core_solver.m_column_types[j] == column_type::free_column; } @@ -902,7 +932,7 @@ bool int_solver::at_bound(unsigned j) const { } } -bool int_solver::at_lower(unsigned j) const { +bool int_solver::at_low(unsigned j) const { auto & mpq_solver = m_lar_solver->m_mpq_lar_core_solver.m_r_solver; switch (mpq_solver.m_column_types[j] ) { case column_type::fixed: @@ -950,4 +980,115 @@ void int_solver::display_row_info(std::ostream & out, unsigned row_index) const rslv.print_column_bound_info(rslv.m_basis[row_index], out); delete it; } + +unsigned int_solver::random() { + return m_lar_solver->get_core_solver().settings().random_next(); +} + +bool int_solver::shift_var(unsigned j, unsigned range) { + if (is_fixed(j) || is_base(j)) + return false; + + bool inf_l, inf_u; + impq l, u; + mpq m; + get_freedom_interval_for_column(j, inf_l, l, inf_u, u, m); + if (inf_l && inf_u) { + impq new_val = impq(random() % (range + 1)); + set_value_for_nbasic_column_ignore_old_values(j, new_val); + return true; + } + if (is_int(j)) { + if (!inf_l) { + l = ceil(l); + if (!m.is_one()) + l = m*ceil(l/m); + } + if (!inf_u) { + u = floor(u); + if (!m.is_one()) + u = m*floor(u/m); + } + } + if (!inf_l && !inf_u && l >= u) + return false; + if (inf_u) { + SASSERT(!inf_l); + impq delta = impq(random() % (range + 1)); + impq new_val = l + m*delta; + set_value_for_nbasic_column_ignore_old_values(j, new_val); + return true; + } + if (inf_l) { + SASSERT(!inf_u); + impq delta = impq(random() % (range + 1)); + impq new_val = u - m*delta; + set_value_for_nbasic_column_ignore_old_values(j, new_val); + return true; + } + if (!is_int(j)) { + SASSERT(!inf_l && !inf_u); + mpq delta = mpq(random() % (range + 1)); + impq new_val = l + ((delta * (u - l)) / mpq(range)); + set_value_for_nbasic_column_ignore_old_values(j, new_val); + return true; + } + else { + mpq r = (u.x - l.x) / m; + if (r < mpq(range)) + range = static_cast(r.get_uint64()); + impq new_val = l + m * (impq(random() % (range + 1))); + set_value_for_nbasic_column_ignore_old_values(j, new_val); + return true; + } +} + +bool int_solver::non_basic_columns_are_at_bounds() const { + auto & lcs = m_lar_solver->m_mpq_lar_core_solver; + for (unsigned j :lcs.m_r_nbasis) { + auto & val = lcs.m_r_x[j]; + switch (lcs.m_column_types()[j]) { + case column_type::boxed: + if (val != lcs.m_r_low_bounds()[j] && val != lcs.m_r_upper_bounds()[j]) + return false; + break; + case column_type::low_bound: + if (val != lcs.m_r_low_bounds()[j]) + return false; + break; + case column_type::upper_bound: + if (val != lcs.m_r_upper_bounds()[j]) + return false; + break; + default: + if (is_int(j) && !val.is_int()) { + return false; + } + } + } + return true; +} + +const impq& int_solver::low_bound(unsigned j) const { + return m_lar_solver->column_low_bound(j); +} + +lia_move int_solver::create_branch_on_column(int j, lar_term& t, mpq& k) const { + lp_assert(t.is_empty()); + lp_assert(j != -1); + t.add_monomial(mpq(1), j); + k = floor(get_value(j)); + TRACE("arith_int", tout << "branching v" << j << " = " << get_value(j) << "\n"; + display_column(tout, j); + tout << "k = " << k << std::endl; + ); + lp_assert(current_solution_is_inf_on_cut(t, k)); + m_lar_solver->subs_term_columns(t); + return lia_move::branch; + +} + +const impq& int_solver::upper_bound(unsigned j) const { + return m_lar_solver->column_upper_bound(j); +} } diff --git a/src/util/lp/int_solver.h b/src/util/lp/int_solver.h index 9820e484c..12723effc 100644 --- a/src/util/lp/int_solver.h +++ b/src/util/lp/int_solver.h @@ -34,9 +34,6 @@ public: int_set m_old_values_set; vector m_old_values_data; unsigned m_branch_cut_counter; - linear_combination_iterator* m_iter_on_gomory_row; - unsigned m_gomory_cut_inf_column; - bool m_found_free_var_in_gomory_row; // methods int_solver(lar_solver* lp); @@ -80,16 +77,17 @@ private: void patch_int_infeasible_columns(); bool get_freedom_interval_for_column(unsigned j, bool & inf_l, impq & l, bool & inf_u, impq & u, mpq & m); linear_combination_iterator * get_column_iterator(unsigned j); - bool lower(unsigned j) const; - bool upper(unsigned j) const; - const impq & lower_bound(unsigned j) const; + const impq & low_bound(unsigned j) const; const impq & upper_bound(unsigned j) const; bool is_int(unsigned j) const; bool is_real(unsigned j) const; bool is_base(unsigned j) const; bool is_boxed(unsigned j) const; + bool is_fixed(unsigned j) const; + bool is_free(unsigned j) const; bool value_is_int(unsigned j) const; - void set_value(unsigned j, const impq & new_val); + void set_value_for_nbasic_column(unsigned j, const impq & new_val); + void set_value_for_nbasic_column_ignore_old_values(unsigned j, const impq & new_val); void fix_non_base_columns(); void failed(); bool is_feasible() const; @@ -102,35 +100,44 @@ private: int find_inf_int_base_column(); int find_inf_int_boxed_base_column_with_smallest_range(); lp_settings& settings(); - void move_non_base_vars_to_bounds(); + bool move_non_base_vars_to_bounds(); void branch_infeasible_int_var(unsigned); - lia_move mk_gomory_cut(lar_term& t, mpq& k,explanation & ex); + lia_move mk_gomory_cut(lar_term& t, mpq& k,explanation & ex, unsigned inf_col, linear_combination_iterator& iter); lia_move report_conflict_from_gomory_cut(mpq & k); - lia_move report_gomory_cut(lar_term& t, mpq& k, mpq& lcm_den, unsigned num_ints); + void adjust_term_and_k_for_some_ints_case_gomory(lar_term& t, mpq& k, mpq& lcm_den); void init_check_data(); bool constrain_free_vars(linear_combination_iterator * r); - lia_move proceed_with_gomory_cut(lar_term& t, mpq& k, explanation& ex); - int find_next_free_var_in_gomory_row(); - bool is_gomory_cut_target(); + lia_move proceed_with_gomory_cut(lar_term& t, mpq& k, explanation& ex, unsigned j, linear_combination_iterator& iter); + int find_free_var_in_gomory_row(linear_combination_iterator& iter); + bool is_gomory_cut_target(linear_combination_iterator &iter); bool at_bound(unsigned j) const; - bool at_lower(unsigned j) const; + bool at_low(unsigned j) const; bool at_upper(unsigned j) const; - + bool has_low(unsigned j) const; + bool has_upper(unsigned j) const; + unsigned row_of_basic_column(unsigned j) const; inline static bool is_rational(const impq & n) { return is_zero(n.y); } inline static mpq fractional_part(const impq & n) { - lp_assert(is_rational); + lp_assert(is_rational(n)); return n.x - floor(n.x); } - void real_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& t, explanation & ex); - void int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& t, explanation& ex, mpq & lcm_den); + void real_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& t, explanation & ex, unsigned inf_column); + void int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& t, explanation& ex, mpq & lcm_den, unsigned inf_column); constraint_index column_upper_bound_constraint(unsigned j) const; constraint_index column_low_bound_constraint(unsigned j) const; void display_row_info(std::ostream & out, unsigned row_index) const; - void gomory_cut_adjust_t_and_k_for_size_1(const vector> & pol, lar_term & t, mpq &k); - void gomory_cut_adjust_t_and_k_for_size_gt_1(vector> & pol, lar_term & t, mpq &k, unsigned num_ints, mpq &lcm_den); + void gomory_cut_adjust_t_and_k(vector> & pol, lar_term & t, mpq &k, bool num_ints, mpq &lcm_den); + bool current_solution_is_inf_on_cut(const lar_term& t, const mpq& k) const; +public: + bool shift_var(unsigned j, unsigned range); +private: + unsigned random(); + bool non_basic_columns_are_at_bounds() const; + bool has_inf_int() const; + lia_move create_branch_on_column(int j, lar_term& t, mpq& k) const; }; } diff --git a/src/util/lp/iterator_on_column.h b/src/util/lp/iterator_on_column.h index c9112a064..0d3b99ef5 100644 --- a/src/util/lp/iterator_on_column.h +++ b/src/util/lp/iterator_on_column.h @@ -57,7 +57,9 @@ struct iterator_on_column:linear_combination_iterator { m_i = -1; } - linear_combination_iterator * clone() override { + bool is_reset() const { return m_i == -1;} + + linear_combination_iterator * clone() { iterator_on_column * r = new iterator_on_column(m_column, m_A); return r; } diff --git a/src/util/lp/iterator_on_indexed_vector.h b/src/util/lp/iterator_on_indexed_vector.h index 6cb98b8f2..b0ef87506 100644 --- a/src/util/lp/iterator_on_indexed_vector.h +++ b/src/util/lp/iterator_on_indexed_vector.h @@ -46,7 +46,10 @@ struct iterator_on_indexed_vector:linear_combination_iterator { void reset() override { m_offset = 0; } - linear_combination_iterator* clone() override { + + bool is_reset() const { return m_offset == 0;} + + linear_combination_iterator* clone() { return new iterator_on_indexed_vector(m_v); } }; diff --git a/src/util/lp/iterator_on_pivot_row.h b/src/util/lp/iterator_on_pivot_row.h index 721502bc2..e20533c55 100644 --- a/src/util/lp/iterator_on_pivot_row.h +++ b/src/util/lp/iterator_on_pivot_row.h @@ -51,7 +51,10 @@ struct iterator_on_pivot_row:linear_combination_iterator { m_basis_returned = false; m_it.reset(); } - linear_combination_iterator * clone() override { + + bool is_reset() const { return m_basis_returned == false && m_it.is_reset();} + + linear_combination_iterator * clone() { iterator_on_pivot_row * r = new iterator_on_pivot_row(m_v, m_basis_j); return r; } diff --git a/src/util/lp/iterator_on_row.h b/src/util/lp/iterator_on_row.h index 55fbda907..36b1d22ad 100644 --- a/src/util/lp/iterator_on_row.h +++ b/src/util/lp/iterator_on_row.h @@ -45,7 +45,10 @@ struct iterator_on_row:linear_combination_iterator { void reset() override { m_i = 0; } - linear_combination_iterator* clone() override { + + bool is_reset() const { return m_i == 0;} + + linear_combination_iterator* clone() { return new iterator_on_row(m_row); } }; diff --git a/src/util/lp/iterator_on_term_with_basis_var.h b/src/util/lp/iterator_on_term_with_basis_var.h index 85d11cf36..cf36046cb 100644 --- a/src/util/lp/iterator_on_term_with_basis_var.h +++ b/src/util/lp/iterator_on_term_with_basis_var.h @@ -60,7 +60,10 @@ struct iterator_on_term_with_basis_var:linear_combination_iterator { m_i++; return true; } - void reset() override { + + bool is_reset() const { return m_term_j_returned == false && m_i == m_term.m_coeffs.begin();} + + void reset() { m_term_j_returned = false; m_i = m_term.m_coeffs.begin(); } diff --git a/src/util/lp/lar_core_solver.h b/src/util/lp/lar_core_solver.h index 4e90e785a..3daf2a3c1 100644 --- a/src/util/lp/lar_core_solver.h +++ b/src/util/lp/lar_core_solver.h @@ -155,7 +155,7 @@ public: void fill_evidence(unsigned row); - + unsigned get_number_of_non_ints() const; void solve(); @@ -344,8 +344,7 @@ public: for (const auto & cc : m_r_solver.m_A.m_columns[j]){ unsigned i = cc.m_i; unsigned jb = m_r_solver.m_basis[i]; - m_r_solver.m_x[jb] -= delta * m_r_solver.m_A.get_val(cc); - m_r_solver.update_column_in_inf_set(jb); + m_r_solver.update_x_with_delta_and_track_feasibility(jb, - delta * m_r_solver.m_A.get_val(cc)); } lp_assert(m_r_solver.A_mult_x_is_off() == false); } diff --git a/src/util/lp/lar_core_solver.hpp b/src/util/lp/lar_core_solver.hpp index 4f1a8cfc2..19833f5eb 100644 --- a/src/util/lp/lar_core_solver.hpp +++ b/src/util/lp/lar_core_solver.hpp @@ -255,14 +255,23 @@ void lar_core_solver::fill_not_improvable_zero_sum() { } } +unsigned lar_core_solver::get_number_of_non_ints() const { + unsigned n = 0; + for (auto & x : m_r_solver.m_x) { + if (x.is_int() == false) + n++; + } + return n; +} void lar_core_solver::solve() { lp_assert(m_r_solver.non_basic_columns_are_set_correctly()); lp_assert(m_r_solver.inf_set_is_correct()); - if (m_r_solver.current_x_is_feasible() && m_r_solver.m_look_for_feasible_solution_only) { - m_r_solver.set_status(lp_status::OPTIMAL); - return; - } + TRACE("find_feas_stats", tout << "infeasibles = " << m_r_solver.m_inf_set.size() << ", int_infs = " << get_number_of_non_ints() << std::endl;); + if (m_r_solver.current_x_is_feasible() && m_r_solver.m_look_for_feasible_solution_only) { + m_r_solver.set_status(lp_status::OPTIMAL); + return; + } ++settings().st().m_need_to_solve_inf; lp_assert(!m_r_solver.A_mult_x_is_off()); lp_assert((!settings().use_tableau()) || r_basis_is_OK()); @@ -278,6 +287,7 @@ void lar_core_solver::solve() { solve_on_signature_tableau(solution_signature, changes_of_basis); else solve_on_signature(solution_signature, changes_of_basis); + lp_assert(!settings().use_tableau() || r_basis_is_OK()); } else { if (!settings().use_tableau()) { diff --git a/src/util/lp/lar_solver.cpp b/src/util/lp/lar_solver.cpp index 712f8a001..191d1ab6f 100644 --- a/src/util/lp/lar_solver.cpp +++ b/src/util/lp/lar_solver.cpp @@ -31,23 +31,22 @@ 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, const impq & x){ - if (!var_is_int(j)) { - lp_assert(m_inf_int_set.contains(j) == false); - return; - } - 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); - }) -{ -} + m_tracker_of_x_change([&](unsigned j) { + call_assignment_tracker(j); + } + ), + m_int_solver(nullptr) +{} -void lar_solver::set_propagate_bounds_on_pivoted_rows_mode(bool v) { +void lar_solver::set_track_pivoted_rows(bool v) { m_mpq_lar_core_solver.m_r_solver.m_pivoted_rows = v? (& m_rows_with_changed_bounds) : nullptr; } +bool lar_solver::get_track_pivoted_rows() const { + return m_mpq_lar_core_solver.m_r_solver.m_pivoted_rows != nullptr; +} + + lar_solver::~lar_solver(){ for (auto c : m_constraints) delete c; @@ -55,8 +54,6 @@ lar_solver::~lar_solver(){ delete t; } -numeric_pair const& lar_solver::get_value(var_index vi) const { return m_mpq_lar_core_solver.m_r_x[vi]; } - bool lar_solver::is_term(var_index j) const { return j >= m_terms_start_index && j - m_terms_start_index < m_terms.size(); } @@ -89,12 +86,6 @@ void lar_solver::print_implied_bound(const implied_bound& be, std::ostream & out out << get_column_name(v); } out << " " << lconstraint_kind_string(be.kind()) << " " << be.m_bound << std::endl; - // for (auto & p : be.m_explanation) { - // out << p.first << " : "; - // print_constraint(p.second, out); - // } - - // m_mpq_lar_core_solver.m_r_solver.print_column_info(be.m_j< m_terms_start_index? be.m_j : adjust_term_index(be.m_j), out); out << "end of implied bound" << std::endl; } @@ -294,6 +285,8 @@ bool lar_solver::has_int_var() const { } lp_status lar_solver::find_feasible_solution() { + TRACE("arith_int", ); + 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(); @@ -307,10 +300,13 @@ lp_status lar_solver::find_feasible_solution() { } m_mpq_lar_core_solver.m_r_solver.m_look_for_feasible_solution_only = true; - return solve(); + auto ret = solve(); + 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; } @@ -321,6 +317,7 @@ lp_status lar_solver::solve() { } m_columns_with_changed_bound.clear(); + lp_assert(inf_int_set_is_correct()); return m_status; } @@ -368,6 +365,8 @@ 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()); @@ -376,15 +375,23 @@ void lar_solver::pop(unsigned k) { for (unsigned j = n_was; j-- > n;) m_ext_vars_to_columns.erase(m_columns_to_ext_vars_or_term_indices[j]); m_columns_to_ext_vars_or_term_indices.resize(n); + TRACE("arith_int", tout << "pop" << std::endl;); if (m_settings.use_tableau()) { pop_tableau(); } + TRACE("arith_int", tout << "pop" << std::endl;); + lp_assert(inf_int_set_is_correct()); + TRACE("arith_int", tout << "pop" << std::endl;); + m_columns_to_ul_pairs.pop(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 || @@ -655,15 +662,15 @@ bool lar_solver::costs_are_used() const { return m_settings.simplex_strategy() != simplex_strategy_enum::tableau_rows; } -void lar_solver::change_basic_x_by_delta_on_column(unsigned j, const numeric_pair & delta) { +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]; - m_mpq_lar_core_solver.m_r_x[bj] -= A_r().get_val(c) * delta; if (tableau_with_costs()) { m_basic_columns_with_changed_cost.insert(bj); } - m_mpq_lar_core_solver.m_r_solver.update_column_in_inf_set(bj); + m_mpq_lar_core_solver.m_r_solver.update_x_with_delta_and_track_feasibility(bj, - A_r().get_val(c) * delta); TRACE("change_x_del", tout << "changed basis column " << bj << ", it is " << ( m_mpq_lar_core_solver.m_r_solver.column_is_feasible(bj)? "feas":"inf") << std::endl;); @@ -675,26 +682,26 @@ void lar_solver::change_basic_x_by_delta_on_column(unsigned j, const numeric_pai m_mpq_lar_core_solver.m_r_solver.solve_Bd(j, m_column_buffer); for (unsigned i : m_column_buffer.m_index) { unsigned bj = m_mpq_lar_core_solver.m_r_basis[i]; - m_mpq_lar_core_solver.m_r_x[bj] -= m_column_buffer[i] * delta; - m_mpq_lar_core_solver.m_r_solver.update_column_in_inf_set(bj); + 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) { if (m_mpq_lar_core_solver.m_r_heading[j] >= 0) { if (costs_are_used()) { bool was_infeas = m_mpq_lar_core_solver.m_r_solver.m_inf_set.contains(j); - m_mpq_lar_core_solver.m_r_solver.update_column_in_inf_set(j); + m_mpq_lar_core_solver.m_r_solver.track_column_feasibility(j); if (was_infeas != m_mpq_lar_core_solver.m_r_solver.m_inf_set.contains(j)) m_basic_columns_with_changed_cost.insert(j); } else { - m_mpq_lar_core_solver.m_r_solver.update_column_in_inf_set(j); + m_mpq_lar_core_solver.m_r_solver.track_column_feasibility(j); } } else { numeric_pair delta; if (m_mpq_lar_core_solver.m_r_solver.make_column_feasible(j, delta)) - change_basic_x_by_delta_on_column(j, delta); + change_basic_columns_dependend_on_a_given_nb_column(j, delta); } } @@ -722,7 +729,6 @@ void lar_solver::update_x_and_inf_costs_for_columns_with_changed_bounds() { } void lar_solver::update_x_and_inf_costs_for_columns_with_changed_bounds_tableau() { - lp_assert(ax_is_correct()); for (auto j : m_columns_with_changed_bound.m_index) update_x_and_inf_costs_for_column_with_changed_bounds(j); @@ -751,6 +757,7 @@ void lar_solver::solve_with_core_solver() { update_x_and_inf_costs_for_columns_with_changed_bounds(); 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()); } @@ -928,6 +935,7 @@ bool lar_solver::all_constraints_hold() const { } bool lar_solver::constraint_holds(const lar_base_constraint & constr, std::unordered_map & var_map) const { + return true; mpq left_side_val = get_left_side_val(constr, var_map); switch (constr.m_kind) { case LE: return left_side_val <= constr.m_right_side; @@ -1246,8 +1254,9 @@ void lar_solver::fill_var_set_for_random_update(unsigned sz, var_index const * v void lar_solver::random_update(unsigned sz, var_index const * vars) { vector column_list; fill_var_set_for_random_update(sz, vars, column_list); - random_updater ru(m_mpq_lar_core_solver, column_list); + random_updater ru(*this, column_list); ru.update(); + lp_assert(inf_int_set_is_correct()); } @@ -1378,6 +1387,7 @@ 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; @@ -1388,8 +1398,10 @@ void lar_solver::clean_inf_set_of_r_solver_after_pop() { } // some basic columns might become non-basic - these columns need to be made feasible numeric_pair delta; - if (m_mpq_lar_core_solver.m_r_solver.make_column_feasible(j, delta)) - change_basic_x_by_delta_on_column(j, delta); + if (m_mpq_lar_core_solver.m_r_solver.make_column_feasible(j, delta)) { + + change_basic_columns_dependend_on_a_given_nb_column(j, delta); + } became_feas.push_back(j); } @@ -1492,6 +1504,7 @@ var_index lar_solver::add_var(unsigned ext_j, bool is_int) { 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; } @@ -1507,6 +1520,7 @@ 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); @@ -1574,7 +1588,7 @@ var_index lar_solver::add_term(const vector> & coeffs, unsigned adjusted_term_index = m_terms.size() - 1; var_index ret = m_terms_start_index + adjusted_term_index; if (use_tableau() && !coeffs.empty()) { - add_row_for_term(m_terms.back(), ret); + add_row_from_term_no_constraint(m_terms.back(), ret); if (m_settings.bound_propagation()) m_rows_with_changed_bounds.insert(A_r().row_count() - 1); } @@ -1582,12 +1596,6 @@ var_index lar_solver::add_term(const vector> & coeffs, return ret; } -void lar_solver::add_row_for_term(const lar_term * term, unsigned term_ext_index) { - lp_assert(sizes_are_correct()); - add_row_from_term_no_constraint(term, term_ext_index); - lp_assert(sizes_are_correct()); -} - void lar_solver::add_row_from_term_no_constraint(const lar_term * term, unsigned term_ext_index) { register_new_ext_var_index(term_ext_index, term_is_int(term)); // j will be a new variable @@ -1604,9 +1612,10 @@ void lar_solver::add_row_from_term_no_constraint(const lar_term * term, unsigned else { fill_last_row_of_A_r(A_r(), term); } - m_mpq_lar_core_solver.m_r_x[j] = get_basic_var_value_from_row_directly(A_r().row_count() - 1); + 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() { @@ -1615,6 +1624,7 @@ 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); @@ -1639,6 +1649,7 @@ 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; } diff --git a/src/util/lp/lar_solver.h b/src/util/lp/lar_solver.h index 41794fef6..47a1936e1 100644 --- a/src/util/lp/lar_solver.h +++ b/src/util/lp/lar_solver.h @@ -81,9 +81,18 @@ class lar_solver : public column_namer { indexed_vector m_column_buffer; public: 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; + } + int_solver * get_int_solver() { + return m_int_solver; + } unsigned constraint_count() const; const lar_base_constraint& get_constraint(unsigned ci) const; - std::function m_tracker_of_x_change; int_set m_inf_int_set; ////////////////// methods //////////////////////////////// static_matrix> & A_r() { return m_mpq_lar_core_solver.m_r_A;} @@ -94,6 +103,14 @@ public: static bool valid_index(unsigned j){ return static_cast(j) >= 0;} bool column_is_int(unsigned j) const; + bool column_value_is_int(unsigned j) const { + return m_mpq_lar_core_solver.m_r_x[j].is_int(); + } + + const impq& get_column_value(unsigned j) const { + return m_mpq_lar_core_solver.m_r_x[j]; + } + bool is_term(var_index j) const; bool column_is_fixed(unsigned j) const; public: @@ -134,87 +151,13 @@ public: void clear(); lar_solver(); - void set_propagate_bounds_on_pivoted_rows_mode(bool v); + void set_track_pivoted_rows(bool v); + bool get_track_pivoted_rows() const; + virtual ~lar_solver(); - void print_implied_bound(const implied_bound& be, std::ostream & out) const { - out << "implied bound\n"; - unsigned v = be.m_j; - if (is_term(v)) { - out << "it is a term number " << be.m_j << std::endl; - print_term(*m_orig_terms[be.m_j - m_terms_start_index], out); - } - else { - out << get_column_name(v); - } - out << " " << lconstraint_kind_string(be.kind()) << " " << be.m_bound << std::endl; - // for (auto & p : be.m_explanation) { - // out << p.first << " : "; - // print_constraint(p.second, out); - // } - - // m_mpq_lar_core_solver.m_r_solver.print_column_info(be.m_j< m_terms_start_index? be.m_j : adjust_term_index(be.m_j), out); - out << "end of implied bound" << std::endl; - } - - bool implied_bound_is_correctly_explained(implied_bound const & be, const vector> & explanation) const { - std::unordered_map coeff_map; - auto rs_of_evidence = zero_of_type(); - unsigned n_of_G = 0, n_of_L = 0; - bool strict = false; - for (auto & it : explanation) { - mpq coeff = it.first; - constraint_index con_ind = it.second; - const auto & constr = *m_constraints[con_ind]; - lconstraint_kind kind = coeff.is_pos() ? constr.m_kind : flip_kind(constr.m_kind); - register_in_map(coeff_map, constr, coeff); - if (kind == GT || kind == LT) - strict = true; - if (kind == GE || kind == GT) n_of_G++; - else if (kind == LE || kind == LT) n_of_L++; - rs_of_evidence += coeff*constr.m_right_side; - } - SASSERT(n_of_G == 0 || n_of_L == 0); - lconstraint_kind kind = n_of_G ? GE : (n_of_L ? LE : EQ); - if (strict) - kind = static_cast((static_cast(kind) / 2)); - - if (!is_term(be.m_j)) { - if (coeff_map.size() != 1) - return false; - auto it = coeff_map.find(be.m_j); - if (it == coeff_map.end()) return false; - mpq ratio = it->second; - if (ratio < zero_of_type()) { - kind = static_cast(-kind); - } - rs_of_evidence /= ratio; - } else { - const lar_term * t = m_orig_terms[adjust_term_index(be.m_j)]; - const auto first_coeff = *t->m_coeffs.begin(); - unsigned j = first_coeff.first; - auto it = coeff_map.find(j); - if (it == coeff_map.end()) - return false; - mpq ratio = it->second / first_coeff.second; - for (auto & p : t->m_coeffs) { - it = coeff_map.find(p.first); - if (it == coeff_map.end()) - return false; - if (p.second * ratio != it->second) - return false; - } - if (ratio < zero_of_type()) { - kind = static_cast(-kind); - } - rs_of_evidence /= ratio; - rs_of_evidence += t->m_v * ratio; - } - - return kind == be.kind() && rs_of_evidence == be.m_bound; - } - + unsigned adjust_term_index(unsigned j) const; void analyze_new_bounds_on_row( unsigned row_index, @@ -707,32 +650,9 @@ public: bool tableau_with_costs() const; - bool costs_are_used() const { - return m_settings.simplex_strategy() != simplex_strategy_enum::tableau_rows; - } - - void change_basic_x_by_delta_on_column(unsigned j, const numeric_pair & delta) { - 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]; - m_mpq_lar_core_solver.m_r_x[bj] -= A_r().get_val(c) * delta; - if (tableau_with_costs()) { - m_basic_columns_with_changed_cost.insert(bj); - } - m_mpq_lar_core_solver.m_r_solver.update_column_in_inf_set(bj); - } - } else { - m_column_buffer.clear(); - m_column_buffer.resize(A_r().row_count()); - m_mpq_lar_core_solver.m_r_solver.solve_Bd(j, m_column_buffer); - for (unsigned i : m_column_buffer.m_index) { - unsigned bj = m_mpq_lar_core_solver.m_r_basis[i]; - m_mpq_lar_core_solver.m_r_x[bj] -= m_column_buffer[i] * delta; - m_mpq_lar_core_solver.m_r_solver.update_column_in_inf_set(bj); - } - } - } - + bool costs_are_used() const; + + void change_basic_columns_dependend_on_a_given_nb_column(unsigned j, const numeric_pair & delta); void update_x_and_inf_costs_for_column_with_changed_bounds(unsigned j); void detect_rows_with_changed_bounds_for_column(unsigned j) { @@ -1261,8 +1181,8 @@ public: void clean_inf_set_of_r_solver_after_pop(); void shrink_explanation_to_minimum(vector> & explanation) const; - bool column_represents_row_in_tableau(unsigned j) { - return m_vars_to_ul_pairs()[j].m_i != static_cast(-1); + bool column_value_is_integer(unsigned j) const { + return get_column_value(j).is_int(); } bool column_is_real(unsigned j) const { @@ -1442,6 +1362,43 @@ bool model_is_int_feasible() const; t = lar_term(pol_after_subs, v); } + 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; + } + 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); + } + + lar_core_solver & get_core_solver() { return m_mpq_lar_core_solver; } }; } diff --git a/src/util/lp/lar_term.h b/src/util/lp/lar_term.h index 1c4a38801..3f83a3fba 100644 --- a/src/util/lp/lar_term.h +++ b/src/util/lp/lar_term.h @@ -25,7 +25,7 @@ struct lar_term { std::unordered_map m_coeffs; mpq m_v; lar_term() {} - void add_monoid(const mpq& c, unsigned j) { + void add_monomial(const mpq& c, unsigned j) { auto it = m_coeffs.find(j); if (it == m_coeffs.end()) { m_coeffs.emplace(j, c); @@ -49,7 +49,7 @@ struct lar_term { lar_term(const vector>& coeffs, const mpq & v) : m_v(v) { for (const auto & p : coeffs) { - add_monoid(p.first, p.second); + add_monomial(p.first, p.second); } } bool operator==(const lar_term & a) const { return false; } // take care not to create identical terms @@ -71,7 +71,7 @@ struct lar_term { if (it == m_coeffs.end()) return; const mpq & b = it->second; for (unsigned it_j :li.m_index) { - add_monoid(- b * li.m_data[it_j], it_j); + add_monomial(- b * li.m_data[it_j], it_j); } m_coeffs.erase(it); } diff --git a/src/util/lp/linear_combination_iterator.h b/src/util/lp/linear_combination_iterator.h index 417bdcf82..739cefdc0 100644 --- a/src/util/lp/linear_combination_iterator.h +++ b/src/util/lp/linear_combination_iterator.h @@ -27,6 +27,7 @@ struct linear_combination_iterator { virtual linear_combination_iterator * clone() = 0; virtual ~linear_combination_iterator(){} virtual unsigned size() const = 0; + virtual bool is_reset() const = 0; }; template struct linear_combination_iterator_on_vector : linear_combination_iterator { @@ -50,7 +51,7 @@ struct linear_combination_iterator_on_vector : linear_combination_iterator { m_offset++; return true; } - + bool is_reset() const { return m_offset == 0;} void reset() {m_offset = 0;} linear_combination_iterator * clone() { return new linear_combination_iterator_on_vector(m_vector); diff --git a/src/util/lp/lp_bound_propagator.h b/src/util/lp/lp_bound_propagator.h index 76870f457..357d9c381 100644 --- a/src/util/lp/lp_bound_propagator.h +++ b/src/util/lp/lp_bound_propagator.h @@ -32,7 +32,7 @@ public: column_type get_column_type(unsigned) const; const impq & get_low_bound(unsigned) const; const impq & get_upper_bound(unsigned) const; - void try_add_bound(const mpq & v, unsigned j, bool is_low, bool coeff_before_j_is_pos, unsigned row_or_term_index, bool strict); + void try_add_bound(mpq v, unsigned j, bool is_low, bool coeff_before_j_is_pos, unsigned row_or_term_index, bool strict); virtual bool bound_is_interesting(unsigned vi, lp::lconstraint_kind kind, const rational & bval) {return true;} diff --git a/src/util/lp/lp_core_solver_base.h b/src/util/lp/lp_core_solver_base.h index 6bc3f2daf..4b239e5de 100644 --- a/src/util/lp/lp_core_solver_base.h +++ b/src/util/lp/lp_core_solver_base.h @@ -45,6 +45,9 @@ public: TRACE("feas", if (m_inf_set.size()) { tout << "column " << m_inf_set.m_index[0] << " is infeasible" << std::endl; + print_column_info(m_inf_set.m_index[0], tout); + } else { + tout << "x is feasible\n"; } ); return m_inf_set.size() == 0; @@ -84,9 +87,9 @@ 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; + std::function * m_tracker_of_x_change; - void set_tracker_of_x(std::function* tracker) { + void set_tracker_of_x(std::function* tracker) { m_tracker_of_x_change = tracker; } @@ -401,6 +404,7 @@ public: } bool make_column_feasible(unsigned j, numeric_pair & delta) { + bool ret = false; lp_assert(m_basis_heading[j] < 0); auto & x = m_x[j]; switch (m_column_types[j]) { @@ -408,43 +412,39 @@ public: lp_assert(m_low_bounds[j] == m_upper_bounds[j]); if (x != m_low_bounds[j]) { delta = m_low_bounds[j] - x; - x = m_low_bounds[j]; - return true; + ret = true;; } break; case column_type::boxed: if (x < m_low_bounds[j]) { delta = m_low_bounds[j] - x; - x = m_low_bounds[j]; - return true; + ret = true;; } if (x > m_upper_bounds[j]) { delta = m_upper_bounds[j] - x; - x = m_upper_bounds[j]; - return true; + ret = true; } break; case column_type::low_bound: if (x < m_low_bounds[j]) { delta = m_low_bounds[j] - x; - x = m_low_bounds[j]; - return true; + ret = true; } break; case column_type::upper_bound: if (x > m_upper_bounds[j]) { delta = m_upper_bounds[j] - x; - x = m_upper_bounds[j]; - return true; + ret = true; } break; - case column_type::free_column: - break; default: - lp_assert(false); break; } - return false; + if (ret) + add_delta_to_x_and_call_tracker(j, delta); + + return ret; + } @@ -584,25 +584,26 @@ public: } void print_column_info(unsigned j, std::ostream & out) const { - out << "column_index = " << j << ", name = "<< column_name(j) << " type = " << column_type_to_string(m_column_types[j]) << std::endl; + out << "j = " << j << ", name = "<< column_name(j); switch (m_column_types[j]) { case column_type::fixed: case column_type::boxed: - out << "(" << m_low_bounds[j] << ", " << m_upper_bounds[j] << ")" << std::endl; + out << " [" << m_low_bounds[j] << ", " << m_upper_bounds[j] << "]"; break; case column_type::low_bound: - out << m_low_bounds[j] << std::endl; + out << " [" << m_low_bounds[j] << "," << "oo" << "]"; break; case column_type::upper_bound: - out << m_upper_bounds[j] << std::endl; + out << " [-oo, " << m_upper_bounds[j] << ']'; break; case column_type::free_column: + out << " [-oo, oo]"; break; default: lp_assert(false); } - out << "basis heading = " << m_basis_heading[j] << std::endl; - out << "x = " << m_x[j] << std::endl; + // out << "basis heading = " << m_basis_heading[j] << std::endl; + out << " x = " << m_x[j] << std::endl; /* std::cout << "cost = " << m_costs[j] << std::endl; std:: cout << "m_d = " << m_d[j] << std::endl;*/ @@ -689,22 +690,43 @@ public: return m_inf_set.contains(j); } - void update_column_in_inf_set(unsigned j) { - if (column_is_feasible(j)) { + void update_x_with_feasibility_tracking(unsigned j, const X & v) { + m_x[j] = v; + track_column_feasibility(j); + } + + void update_x_with_delta_and_track_feasibility(unsigned j, const X & del) { + m_x[j] += del; + track_column_feasibility(j); + } + + 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) { + if (column_is_feasible(j)) remove_column_from_inf_set(j); - } else { + else 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_x[j]); + (*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_x[j]); + (*m_tracker_of_x_change)(j); m_inf_set.erase(j); lp_assert(column_is_feasible(j)); } diff --git a/src/util/lp/lp_primal_core_solver.h b/src/util/lp/lp_primal_core_solver.h index 8df141551..cee5b1ef2 100644 --- a/src/util/lp/lp_primal_core_solver.h +++ b/src/util/lp/lp_primal_core_solver.h @@ -445,7 +445,7 @@ public: void advance_on_entering_and_leaving_tableau_rows(int entering, int leaving, const X &theta ) { this->update_basis_and_x_tableau(entering, leaving, theta); - this->update_column_in_inf_set(entering); + this->track_column_feasibility(entering); } diff --git a/src/util/lp/lp_primal_core_solver.hpp b/src/util/lp/lp_primal_core_solver.hpp index 000d274e1..eda85175b 100644 --- a/src/util/lp/lp_primal_core_solver.hpp +++ b/src/util/lp/lp_primal_core_solver.hpp @@ -864,7 +864,7 @@ template void lp_primal_core_solver::print_column template unsigned lp_primal_core_solver::solve() { if (numeric_traits::precise() && this->m_settings.use_tableau()) return solve_with_tableau(); - + init_run(); if (this->current_x_is_feasible() && this->m_look_for_feasible_solution_only) { this->set_status(lp_status::FEASIBLE); @@ -880,6 +880,7 @@ template unsigned lp_primal_core_solver::solve() return this->total_iterations(); } one_iteration(); + lp_assert(!this->m_using_infeas_costs || this->costs_on_nbasis_are_zeros()); switch (this->get_status()) { case lp_status::OPTIMAL: // double check that we are at optimum diff --git a/src/util/lp/lp_primal_core_solver_tableau.h b/src/util/lp/lp_primal_core_solver_tableau.h index b6566f7c3..77cebc821 100644 --- a/src/util/lp/lp_primal_core_solver_tableau.h +++ b/src/util/lp/lp_primal_core_solver_tableau.h @@ -112,8 +112,9 @@ unsigned lp_primal_core_solver::solve_with_tableau() { if (this->print_statistics_with_iterations_and_nonzeroes_and_cost_and_check_that_the_time_is_over((this->m_using_infeas_costs? "inf t" : "feas t"), * this->m_settings.get_message_ostream())) { return this->total_iterations(); } - if (this->m_settings.use_tableau_rows()) + if (this->m_settings.use_tableau_rows()) { one_iteration_tableau_rows(); + } else one_iteration_tableau(); switch (this->get_status()) { @@ -352,22 +353,20 @@ update_basis_and_x_tableau(int entering, int leaving, X const & tt) { } template void lp_primal_core_solver:: update_x_tableau(unsigned entering, const X& delta) { + this->add_delta_to_x_and_call_tracker(entering, delta); if (!this->m_using_infeas_costs) { - this->m_x[entering] += delta; for (const auto & c : this->m_A.m_columns[entering]) { unsigned i = c.m_i; - this->m_x[this->m_basis[i]] -= delta * this->m_A.get_val(c); - this->update_column_in_inf_set(this->m_basis[i]); + this->update_x_with_delta_and_track_feasibility(this->m_basis[i], - delta * this->m_A.get_val(c)); } } else { // m_using_infeas_costs == true - this->m_x[entering] += delta; lp_assert(this->column_is_feasible(entering)); lp_assert(this->m_costs[entering] == zero_of_type()); // m_d[entering] can change because of the cost change for basic columns. for (const auto & c : this->m_A.m_columns[entering]) { unsigned i = c.m_i; unsigned j = this->m_basis[i]; - this->m_x[j] -= delta * this->m_A.get_val(c); + this->add_delta_to_x_and_call_tracker(j, -delta * this->m_A.get_val(c)); update_inf_cost_for_column_tableau(j); if (is_zero(this->m_costs[j])) this->remove_column_from_inf_set(j); diff --git a/src/util/lp/lp_settings.h b/src/util/lp/lp_settings.h index 063001d53..08dbb86af 100644 --- a/src/util/lp/lp_settings.h +++ b/src/util/lp/lp_settings.h @@ -221,7 +221,8 @@ public: max_row_length_for_bound_propagation(300), backup_costs(true), column_number_threshold_for_using_lu_in_lar_solver(4000), - m_int_branch_cut_threshold(10000000) + m_int_branch_cut_gomory_threshold(4), + m_run_gcd_test(true) {} void set_resource_limit(lp_resource_limit& lim) { m_resource_limit = &lim; } @@ -324,11 +325,12 @@ public: #endif bool use_breakpoints_in_feasibility_search; unsigned random_next() { return m_rand(); } - void random_seed(unsigned s) { m_rand.set_seed(s); } + 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_threshold; + unsigned m_int_branch_cut_gomory_threshold; + bool m_run_gcd_test; }; // end of lp_settings class @@ -351,7 +353,7 @@ inline std::string T_to_string(const numeric_pair & t) { inline std::string T_to_string(const mpq & t) { std::ostringstream strs; - strs << t.get_double(); + strs << t; return strs.str(); } diff --git a/src/util/lp/random_updater.h b/src/util/lp/random_updater.h index 6b03ad941..441038c00 100644 --- a/src/util/lp/random_updater.h +++ b/src/util/lp/random_updater.h @@ -30,67 +30,19 @@ Revision History: // The class searches for a feasible solution with as many different values of variables as it can find namespace lp { template struct numeric_pair; // forward definition -class lar_core_solver; // forward definition +class lar_solver; // forward definition class random_updater { - struct interval { - bool upper_bound_is_set; - numeric_pair upper_bound; - bool low_bound_is_set; - numeric_pair low_bound; - interval() : upper_bound_is_set(false), - low_bound_is_set(false) {} - - void set_low_bound(const numeric_pair & v) { - if (low_bound_is_set) { - low_bound = std::max(v, low_bound); - } else { - low_bound = v; - low_bound_is_set = true; - } - } - void set_upper_bound(const numeric_pair & v) { - if (upper_bound_is_set) { - upper_bound = std::min(v, upper_bound); - } else { - upper_bound = v; - upper_bound_is_set = true; - } - } - bool is_empty() const {return - upper_bound_is_set && low_bound_is_set && low_bound >= upper_bound; - } - - bool low_bound_holds(const numeric_pair & a) const { - return low_bound_is_set == false || a >= low_bound; - } - bool upper_bound_holds(const numeric_pair & a) const { - return upper_bound_is_set == false || a <= upper_bound; - } - - bool contains(const numeric_pair & a) const { - return low_bound_holds(a) && upper_bound_holds(a); - } - std::string lbs() { return low_bound_is_set ? T_to_string(low_bound):std::string("inf");} - std::string rbs() { return upper_bound_is_set? T_to_string(upper_bound):std::string("inf");} - std::string to_str() { return std::string("[")+ lbs() + ", " + rbs() + "]";} - }; std::set m_var_set; - lar_core_solver & m_core_solver; - unsigned range; - linear_combination_iterator* m_column_j; // the actual column - interval find_shift_interval(unsigned j); - interval get_interval_of_non_basic_var(unsigned j); + lar_solver & m_lar_solver; + unsigned m_range; void add_column_to_sets(unsigned j); - void random_shift_var(unsigned j); + bool random_shift_var(unsigned j); std::unordered_map, unsigned> m_values; // it maps a value to the number of time it occurs - void diminish_interval_to_leave_basic_vars_feasible(numeric_pair &nb_x, interval & inter); - void shift_var(unsigned j, interval & r); - void diminish_interval_for_basic_var(numeric_pair &nb_x, unsigned j, mpq & a, interval & r); - numeric_pair get_random_from_interval(interval & r); - void add_value(numeric_pair& v); - void remove_value(numeric_pair & v); + bool shift_var(unsigned j); + void add_value(const numeric_pair& v); + void remove_value(const numeric_pair & v); public: - random_updater(lar_core_solver & core_solver, const vector & column_list); + random_updater(lar_solver & solver, const vector & column_list); void update(); }; } diff --git a/src/util/lp/random_updater.hpp b/src/util/lp/random_updater.hpp index 9305569f8..664ad4eca 100644 --- a/src/util/lp/random_updater.hpp +++ b/src/util/lp/random_updater.hpp @@ -26,156 +26,25 @@ namespace lp { random_updater::random_updater( - lar_core_solver & lar_core_solver, + lar_solver & lar_solver, const vector & column_indices) : - m_core_solver(lar_core_solver), - range(100000) { + m_lar_solver(lar_solver), + m_range(100000) { for (unsigned j : column_indices) add_column_to_sets(j); } -random_updater::interval random_updater::get_interval_of_non_basic_var(unsigned j) { - interval ret; - switch (m_core_solver.get_column_type(j)) { - case column_type::free_column: - break; - case column_type::low_bound: - ret.set_low_bound(m_core_solver.m_r_low_bounds[j]); - break; - case column_type::upper_bound: - ret.set_upper_bound(m_core_solver.m_r_upper_bounds[j]); - break; - case column_type::boxed: - case column_type::fixed: - ret.set_low_bound(m_core_solver.m_r_low_bounds[j]); - ret.set_upper_bound(m_core_solver.m_r_upper_bounds[j]); - break; - default: - lp_assert(false); - } - return ret; + +bool random_updater::shift_var(unsigned v) { + return m_lar_solver.get_int_solver()->shift_var(v, m_range); } -void random_updater::diminish_interval_for_basic_var(numeric_pair& nb_x, unsigned j, - mpq & a, - interval & r) { - lp_assert(m_core_solver.m_r_heading[j] >= 0); - numeric_pair delta; - lp_assert(a != zero_of_type()); - switch (m_core_solver.get_column_type(j)) { - case column_type::free_column: - break; - case column_type::low_bound: - delta = m_core_solver.m_r_x[j] - m_core_solver.m_r_low_bounds[j]; - lp_assert(delta >= zero_of_type>()); - if (a > 0) { - r.set_upper_bound(nb_x + delta / a); - } else { - r.set_low_bound(nb_x + delta / a); - } - break; - case column_type::upper_bound: - delta = m_core_solver.m_r_upper_bounds()[j] - m_core_solver.m_r_x[j]; - lp_assert(delta >= zero_of_type>()); - if (a > 0) { - r.set_low_bound(nb_x - delta / a); - } else { - r.set_upper_bound(nb_x - delta / a); - } - break; - case column_type::boxed: - if (a > 0) { - delta = m_core_solver.m_r_x[j] - m_core_solver.m_r_low_bounds[j]; - lp_assert(delta >= zero_of_type>()); - r.set_upper_bound(nb_x + delta / a); - delta = m_core_solver.m_r_upper_bounds()[j] - m_core_solver.m_r_x[j]; - lp_assert(delta >= zero_of_type>()); - r.set_low_bound(nb_x - delta / a); - } else { // a < 0 - delta = m_core_solver.m_r_upper_bounds()[j] - m_core_solver.m_r_x[j]; - lp_assert(delta >= zero_of_type>()); - r.set_upper_bound(nb_x - delta / a); - delta = m_core_solver.m_r_x[j] - m_core_solver.m_r_low_bounds[j]; - lp_assert(delta >= zero_of_type>()); - r.set_low_bound(nb_x + delta / a); - } - break; - case column_type::fixed: - r.set_low_bound(nb_x); - r.set_upper_bound(nb_x); - break; - default: - lp_assert(false); - } -} - - -void random_updater::diminish_interval_to_leave_basic_vars_feasible(numeric_pair &nb_x, interval & r) { - m_column_j->reset(); - unsigned i; - mpq a; - while (m_column_j->next(a, i)) { - diminish_interval_for_basic_var(nb_x, m_core_solver.m_r_basis[i], a, r); - if (r.is_empty()) - break; - } -} - -random_updater::interval random_updater::find_shift_interval(unsigned j) { - interval ret = get_interval_of_non_basic_var(j); - diminish_interval_to_leave_basic_vars_feasible(m_core_solver.m_r_x[j], ret); - return ret; -} - -void random_updater::shift_var(unsigned j, interval & r) { - lp_assert(r.contains(m_core_solver.m_r_x[j])); - lp_assert(m_core_solver.m_r_solver.column_is_feasible(j)); - auto old_x = m_core_solver.m_r_x[j]; - remove_value(old_x); - auto new_val = m_core_solver.m_r_x[j] = get_random_from_interval(r); - add_value(new_val); - - lp_assert(r.contains(m_core_solver.m_r_x[j])); - lp_assert(m_core_solver.m_r_solver.column_is_feasible(j)); - auto delta = m_core_solver.m_r_x[j] - old_x; - - unsigned i; - m_column_j->reset(); - mpq a; - while(m_column_j->next(a, i)) { - unsigned bj = m_core_solver.m_r_basis[i]; - m_core_solver.m_r_x[bj] -= a * delta; - lp_assert(m_core_solver.m_r_solver.column_is_feasible(bj)); - } - lp_assert(m_core_solver.m_r_solver.A_mult_x_is_off() == false); -} - -numeric_pair random_updater::get_random_from_interval(interval & r) { - unsigned rand = m_core_solver.settings().random_next(); - if ((!r.low_bound_is_set) && (!r.upper_bound_is_set)) - return numeric_pair(rand % range, 0); - if (r.low_bound_is_set && (!r.upper_bound_is_set)) - return r.low_bound + numeric_pair(rand % range, 0); - if ((!r.low_bound_is_set) && r.upper_bound_is_set) - return r.upper_bound - numeric_pair(rand % range, 0); - lp_assert(r.low_bound_is_set && r.upper_bound_is_set); - return r.low_bound + (rand % range) * (r.upper_bound - r.low_bound)/ range; -} - -void random_updater::random_shift_var(unsigned j) { - m_column_j = m_core_solver.get_column_iterator(j); - if (m_column_j->size() >= 50) { - delete m_column_j; - return; - } - interval interv = find_shift_interval(j); - if (interv.is_empty()) { - delete m_column_j; - return; +bool random_updater::random_shift_var(unsigned j) { + if (m_lar_solver.A_r().m_columns.size() >= 50) { + return false; } - shift_var(j, interv); - delete m_column_j; + return shift_var(j); } void random_updater::update() { @@ -183,11 +52,15 @@ void random_updater::update() { if (m_var_set.size() <= m_values.size()) { break; // we are done } - random_shift_var(j); + auto old_x = m_lar_solver.get_column_value(j); + if (random_shift_var(j)) { + remove_value(old_x); + add_value(m_lar_solver.get_column_value(j)); + } } } -void random_updater::add_value(numeric_pair& v) { +void random_updater::add_value(const numeric_pair& v) { auto it = m_values.find(v); if (it == m_values.end()) { m_values[v] = 1; @@ -196,7 +69,7 @@ void random_updater::add_value(numeric_pair& v) { } } -void random_updater::remove_value(numeric_pair& v) { +void random_updater::remove_value(const numeric_pair& v) { std::unordered_map, unsigned>::iterator it = m_values.find(v); lp_assert(it != m_values.end()); it->second--; @@ -205,16 +78,16 @@ void random_updater::remove_value(numeric_pair& v) { } void random_updater::add_column_to_sets(unsigned j) { - if (m_core_solver.m_r_heading[j] < 0) { + if (m_lar_solver.get_core_solver().m_r_heading[j] < 0) { m_var_set.insert(j); - add_value(m_core_solver.m_r_x[j]); + add_value(m_lar_solver.get_core_solver().m_r_x[j]); } else { - unsigned row = m_core_solver.m_r_heading[j]; - for (auto row_c : m_core_solver.m_r_A.m_rows[row]) { + unsigned row = m_lar_solver.get_core_solver().m_r_heading[j]; + for (auto row_c : m_lar_solver.get_core_solver().m_r_A.m_rows[row]) { unsigned cj = row_c.m_j; - if (m_core_solver.m_r_heading[cj] < 0) { + if (m_lar_solver.get_core_solver().m_r_heading[cj] < 0) { m_var_set.insert(cj); - add_value(m_core_solver.m_r_x[cj]); + add_value(m_lar_solver.get_core_solver().m_r_x[cj]); } } }