From e705e5a3094a0ffde405d32afaae9d26492fe6fe Mon Sep 17 00:00:00 2001 From: Lev Date: Thu, 13 Sep 2018 11:38:22 -0700 Subject: [PATCH 01/11] branch on inf basic in gomory Signed-off-by: Lev --- src/util/lp/gomory.h | 46 ++++++++++++++++++++++++++++++++++++++ src/util/lp/int_solver.cpp | 17 +++----------- src/util/lp/int_solver.h | 1 - 3 files changed, 49 insertions(+), 15 deletions(-) create mode 100644 src/util/lp/gomory.h diff --git a/src/util/lp/gomory.h b/src/util/lp/gomory.h new file mode 100644 index 000000000..87879f9f1 --- /dev/null +++ b/src/util/lp/gomory.h @@ -0,0 +1,46 @@ +/*++ +Copyright (c) 2017 Microsoft Corporation + +Module Name: + + + +Abstract: + + + +Author: + Nikolaj Bjorner (nbjorner) + Lev Nachmanson (levnach) + +Revision History: + + +--*/ +#pragma once +#include "util/lp/lar_term.h" +#include "util/lp/lia_move.h" +#include "util/lp/explanation.h" + +namespace lp { +class gomory { + lar_term & m_t; // the term to return in the cut + mpq & m_k; // the right side of the cut + explanation& m_ex; // the conflict explanation + bool & m_upper; // we have a cut m_t*x <= k if m_upper is true nad m_t*x >= k otherwise + unsigned m_basic_inf_int_j; // a basis column which has to be an integer but has a not integral value + const row_strip& m_row +public : + gomory(lar_term & m_t, + mpq & m_k, + explanation& m_ex, + bool & m_upper, + unsigned m_basic_inf_int_j ) : + m_t(t), + m_k(k), + m_ex(ex), + m_upper(upper), + m_basic_inf_int_j(j) { + } +}; +} diff --git a/src/util/lp/int_solver.cpp b/src/util/lp/int_solver.cpp index a77c202a0..c416c63f3 100644 --- a/src/util/lp/int_solver.cpp +++ b/src/util/lp/int_solver.cpp @@ -102,6 +102,8 @@ bool int_solver::is_gomory_cut_target(const row_strip& row) { for (const auto & p : row) { j = p.var(); if (is_base(j)) continue; + if (is_free(j)) + return false; if (!at_bound(j)) return false; if (!is_zero(get_value(j).y)) { @@ -350,24 +352,11 @@ lia_move int_solver::mk_gomory_cut( unsigned inf_col, const row_strip & row return lia_move::cut; } -int int_solver::find_free_var_in_gomory_row(const row_strip& row) { - unsigned j; - for (const auto & p : row) { - j = p.var(); - if (!is_base(j) && is_free(j)) - return static_cast(j); - } - return -1; -} - lia_move int_solver::proceed_with_gomory_cut(unsigned j) { const row_strip& row = m_lar_solver->get_row(row_of_basic_column(j)); - if (-1 != find_free_var_in_gomory_row(row)) - return lia_move::undef; - if (!is_gomory_cut_target(row)) - return lia_move::undef; + return create_branch_on_column(j); *m_upper = true; return mk_gomory_cut(j, row); diff --git a/src/util/lp/int_solver.h b/src/util/lp/int_solver.h index ec708918d..82fcb6eb4 100644 --- a/src/util/lp/int_solver.h +++ b/src/util/lp/int_solver.h @@ -107,7 +107,6 @@ private: lia_move report_conflict_from_gomory_cut(); void adjust_term_and_k_for_some_ints_case_gomory(mpq& lcm_den); lia_move proceed_with_gomory_cut(unsigned j); - int find_free_var_in_gomory_row(const row_strip& ); bool is_gomory_cut_target(const row_strip&); bool at_bound(unsigned j) const; bool at_low(unsigned j) const; From 5dee39721a73d60c67395a78451aecaf481be2f0 Mon Sep 17 00:00:00 2001 From: Lev Date: Fri, 14 Sep 2018 11:52:14 -0700 Subject: [PATCH 02/11] rebase Signed-off-by: Lev --- src/util/lp/bound_propagator.cpp | 4 ---- src/util/lp/lar_constraints.h | 2 +- src/util/lp/lar_solver.cpp | 25 ++++++++++--------------- src/util/lp/lar_solver.h | 8 +++----- src/util/lp/lar_term.h | 9 +++------ 5 files changed, 17 insertions(+), 31 deletions(-) diff --git a/src/util/lp/bound_propagator.cpp b/src/util/lp/bound_propagator.cpp index c4fa2aefa..a5c7c976a 100644 --- a/src/util/lp/bound_propagator.cpp +++ b/src/util/lp/bound_propagator.cpp @@ -17,10 +17,6 @@ const impq & bound_propagator::get_upper_bound(unsigned j) const { } 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) diff --git a/src/util/lp/lar_constraints.h b/src/util/lp/lar_constraints.h index ac15028bb..e06c6d1c7 100644 --- a/src/util/lp/lar_constraints.h +++ b/src/util/lp/lar_constraints.h @@ -75,7 +75,7 @@ struct lar_term_constraint: public lar_base_constraint { } unsigned size() const override { return m_term->size();} lar_term_constraint(const lar_term *t, lconstraint_kind kind, const mpq& right_side) : lar_base_constraint(kind, right_side), m_term(t) { } - mpq get_free_coeff_of_left_side() const override { return m_term->m_v;} + mpq get_free_coeff_of_left_side() const override { return zero_of_type();} }; diff --git a/src/util/lp/lar_solver.cpp b/src/util/lp/lar_solver.cpp index 1340d1826..bc4c17d65 100644 --- a/src/util/lp/lar_solver.cpp +++ b/src/util/lp/lar_solver.cpp @@ -137,7 +137,6 @@ bool lar_solver::implied_bound_is_correctly_explained(implied_bound const & be, 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; @@ -613,7 +612,6 @@ void lar_solver::substitute_terms_in_linear_expression(const vector(); for (auto const& cv : t) { impq const& r = get_column_value(cv.var()); if (!numeric_traits::is_zero(r.y)) return false; @@ -1497,7 +1495,7 @@ bool lar_solver::term_is_int(const lar_term * t) const { for (auto const & p : t->m_coeffs) if (! (column_is_int(p.first) && p.second.is_int())) return false; - return t->m_v.is_int(); + return true; } bool lar_solver::var_is_int(var_index v) const { @@ -1598,9 +1596,8 @@ void lar_solver::add_new_var_to_core_fields_for_mpq(bool register_in_basis) { } -var_index lar_solver::add_term_undecided(const vector> & coeffs, - const mpq &m_v) { - push_and_register_term(new lar_term(coeffs, m_v)); +var_index lar_solver::add_term_undecided(const vector> & coeffs) { + push_and_register_term(new lar_term(coeffs)); return m_terms_start_index + m_terms.size() - 1; } @@ -1643,12 +1640,11 @@ void lar_solver::push_and_register_term(lar_term* t) { } // terms -var_index lar_solver::add_term(const vector> & coeffs, - const mpq &m_v) { +var_index lar_solver::add_term(const vector> & coeffs) { if (strategy_is_undecided()) - return add_term_undecided(coeffs, m_v); + return add_term_undecided(coeffs); - push_and_register_term(new lar_term(coeffs, m_v)); + push_and_register_term(new lar_term(coeffs)); unsigned adjusted_term_index = m_terms.size() - 1; var_index ret = m_terms_start_index + adjusted_term_index; if (use_tableau() && !coeffs.empty()) { @@ -1743,9 +1739,8 @@ void lar_solver::add_var_bound_on_constraint_for_term(var_index j, lconstraint_k // lp_assert(!term_is_int(m_terms[adjusted_term_index]) || right_side.is_int()); unsigned term_j; if (m_var_register.external_is_used(j, term_j)) { - mpq rs = right_side - m_terms[adjusted_term_index]->m_v; m_constraints.push_back(new lar_term_constraint(m_terms[adjusted_term_index], kind, right_side)); - update_column_type_and_bound(term_j, kind, rs, ci); + update_column_type_and_bound(term_j, kind, right_side, ci); } else { add_constraint_from_term_and_create_new_column_row(j, m_terms[adjusted_term_index], kind, right_side); @@ -1756,7 +1751,7 @@ constraint_index lar_solver::add_constraint(const vector> left_side; mpq rs = -right_side_parm; substitute_terms_in_linear_expression(left_side_with_terms, left_side, rs); - unsigned term_index = add_term(left_side, zero_of_type()); + unsigned term_index = add_term(left_side); constraint_index ci = m_constraints.size(); add_var_bound_on_constraint_for_term(term_index, kind_par, -rs, ci); return ci; @@ -1767,7 +1762,7 @@ void lar_solver::add_constraint_from_term_and_create_new_column_row(unsigned ter add_row_from_term_no_constraint(term, term_j); unsigned j = A_r().column_count() - 1; - update_column_type_and_bound(j, kind, right_side - term->m_v, m_constraints.size()); + update_column_type_and_bound(j, kind, right_side, m_constraints.size()); m_constraints.push_back(new lar_term_constraint(term, kind, right_side)); lp_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_costs.size()); } diff --git a/src/util/lp/lar_solver.h b/src/util/lp/lar_solver.h index 9afb70c72..2e85513ff 100644 --- a/src/util/lp/lar_solver.h +++ b/src/util/lp/lar_solver.h @@ -164,13 +164,11 @@ public: // terms - var_index add_term(const vector> & coeffs, - const mpq &m_v); + var_index add_term(const vector> & coeffs); - var_index add_term_undecided(const vector> & coeffs, - const mpq &m_v); + var_index add_term_undecided(const vector> & coeffs); - bool term_coeffs_are_ok(const vector> & coeffs, const mpq& v); + bool term_coeffs_are_ok(const vector> & coeffs); void push_and_register_term(lar_term* t); void add_row_for_term(const lar_term * term, unsigned term_ext_index); diff --git a/src/util/lp/lar_term.h b/src/util/lp/lar_term.h index 519847848..47861baf8 100644 --- a/src/util/lp/lar_term.h +++ b/src/util/lp/lar_term.h @@ -23,7 +23,6 @@ namespace lp { struct lar_term { // the term evaluates to sum of m_coeffs + m_v std::unordered_map m_coeffs; - mpq m_v; lar_term() {} void add_monomial(const mpq& c, unsigned j) { auto it = m_coeffs.find(j); @@ -37,7 +36,7 @@ struct lar_term { } bool is_empty() const { - return m_coeffs.size() == 0 && is_zero(m_v); + return m_coeffs.size() == 0; } unsigned size() const { return static_cast(m_coeffs.size()); } @@ -46,8 +45,7 @@ struct lar_term { return m_coeffs; } - lar_term(const vector>& coeffs, - const mpq & v) : m_v(v) { + lar_term(const vector>& coeffs) { for (const auto & p : coeffs) { add_monomial(p.first, p.second); } @@ -87,7 +85,7 @@ struct lar_term { template T apply(const vector& x) const { - T ret = T(m_v); + T ret = zero_of_type(); for (const auto & t : m_coeffs) { ret += t.second * x[t.first]; } @@ -96,7 +94,6 @@ struct lar_term { void clear() { m_coeffs.clear(); - m_v = zero_of_type(); } struct ival { From 57357b7ecebec0d6065658210f1bf4c5a573c6e2 Mon Sep 17 00:00:00 2001 From: Lev Nachmanson Date: Thu, 13 Sep 2018 17:39:06 -0700 Subject: [PATCH 03/11] does not compile Signed-off-by: Lev Nachmanson --- src/smt/theory_lra.cpp | 29 +++++++++++++++++++++-------- 1 file changed, 21 insertions(+), 8 deletions(-) diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp index f8e2f9fe1..a5190459f 100644 --- a/src/smt/theory_lra.cpp +++ b/src/smt/theory_lra.cpp @@ -54,6 +54,17 @@ std::ostream& operator<<(std::ostream& out, bound_kind const& k) { return out; } +struct term_info { + rational m_offset; + lp::var_index m_var; + + rational const& offset() const { return m_offset; } + lp::var_index var() const { return m_var; } + + term_info(): m_var(-1) {} + term_info(lp:var_index vi, rational const& offset): m_offset(offset), m_var(vi) {} +}; + class bound { smt::bool_var m_bv; smt::theory_var m_var; @@ -777,8 +788,8 @@ class theory_lra::imp { lp::var_index vi = m_theory_var2var_index.get(v, UINT_MAX); TRACE("arith", tout << mk_pp(term, m) << " " << v << " " << vi << "\n";); if (vi == UINT_MAX) { - vi = m_solver->add_term(m_left_side, st.coeff()); - m_theory_var2var_index.setx(v, vi, UINT_MAX); + vi = m_solver->add_term(m_left_side); + m_theory_var2var_index.setx(v, term_info(vi, st.coeff()), term_info(0, rational::zero())); if (m_solver->is_term(vi)) { m_term_index2theory_var.setx(m_solver->adjust_term_index(vi), v, UINT_MAX); } @@ -1187,17 +1198,19 @@ public: lp::impq get_ivalue(theory_var v) const { SASSERT(can_get_ivalue(v)); lp::var_index vi = m_theory_var2var_index[v]; - if (!m_solver->is_term(vi)) - return m_solver->get_column_value(vi); - m_todo_terms.push_back(std::make_pair(vi, rational::one())); - lp::impq result(0); + term_info ti = m_theory_var2var_index[v]; + if (!m_solver->is_term(ti.var())) + return m_solver->get_column_value(ti.var()) + ti.offset(); + m_todo_terms.push_back(std::make_pair(ti, rational::one())); + lp::impq result(ti); while (!m_todo_terms.empty()) { - vi = m_todo_terms.back().first; + ti = m_todo_terms.back().first; + vi = ti.var(); rational coeff = m_todo_terms.back().second; m_todo_terms.pop_back(); + result += ti.offset() * coeff; 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_todo_terms.push_back(std::make_pair(i.var(), coeff * i.coeff())); } From 22213a9e73567983fa058c2c0a24b45ba74c932b Mon Sep 17 00:00:00 2001 From: Lev Date: Fri, 14 Sep 2018 11:53:54 -0700 Subject: [PATCH 04/11] rebase Signed-off-by: Lev --- src/smt/theory_lra.cpp | 29 ++++++++--------------------- src/util/lp/bound_propagator.cpp | 4 ++++ src/util/lp/lar_constraints.h | 2 +- src/util/lp/lar_solver.cpp | 25 +++++++++++++++---------- src/util/lp/lar_solver.h | 8 +++++--- src/util/lp/lar_term.h | 9 ++++++--- 6 files changed, 39 insertions(+), 38 deletions(-) diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp index a5190459f..f8e2f9fe1 100644 --- a/src/smt/theory_lra.cpp +++ b/src/smt/theory_lra.cpp @@ -54,17 +54,6 @@ std::ostream& operator<<(std::ostream& out, bound_kind const& k) { return out; } -struct term_info { - rational m_offset; - lp::var_index m_var; - - rational const& offset() const { return m_offset; } - lp::var_index var() const { return m_var; } - - term_info(): m_var(-1) {} - term_info(lp:var_index vi, rational const& offset): m_offset(offset), m_var(vi) {} -}; - class bound { smt::bool_var m_bv; smt::theory_var m_var; @@ -788,8 +777,8 @@ class theory_lra::imp { lp::var_index vi = m_theory_var2var_index.get(v, UINT_MAX); TRACE("arith", tout << mk_pp(term, m) << " " << v << " " << vi << "\n";); if (vi == UINT_MAX) { - vi = m_solver->add_term(m_left_side); - m_theory_var2var_index.setx(v, term_info(vi, st.coeff()), term_info(0, rational::zero())); + vi = m_solver->add_term(m_left_side, st.coeff()); + m_theory_var2var_index.setx(v, vi, UINT_MAX); if (m_solver->is_term(vi)) { m_term_index2theory_var.setx(m_solver->adjust_term_index(vi), v, UINT_MAX); } @@ -1198,19 +1187,17 @@ public: lp::impq get_ivalue(theory_var v) const { SASSERT(can_get_ivalue(v)); lp::var_index vi = m_theory_var2var_index[v]; - term_info ti = m_theory_var2var_index[v]; - if (!m_solver->is_term(ti.var())) - return m_solver->get_column_value(ti.var()) + ti.offset(); - m_todo_terms.push_back(std::make_pair(ti, rational::one())); - lp::impq result(ti); + if (!m_solver->is_term(vi)) + 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()) { - ti = m_todo_terms.back().first; - vi = ti.var(); + vi = m_todo_terms.back().first; rational coeff = m_todo_terms.back().second; m_todo_terms.pop_back(); - result += ti.offset() * coeff; 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_todo_terms.push_back(std::make_pair(i.var(), coeff * i.coeff())); } diff --git a/src/util/lp/bound_propagator.cpp b/src/util/lp/bound_propagator.cpp index a5c7c976a..c4fa2aefa 100644 --- a/src/util/lp/bound_propagator.cpp +++ b/src/util/lp/bound_propagator.cpp @@ -17,6 +17,10 @@ const impq & bound_propagator::get_upper_bound(unsigned j) const { } 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) diff --git a/src/util/lp/lar_constraints.h b/src/util/lp/lar_constraints.h index e06c6d1c7..ac15028bb 100644 --- a/src/util/lp/lar_constraints.h +++ b/src/util/lp/lar_constraints.h @@ -75,7 +75,7 @@ struct lar_term_constraint: public lar_base_constraint { } unsigned size() const override { return m_term->size();} lar_term_constraint(const lar_term *t, lconstraint_kind kind, const mpq& right_side) : lar_base_constraint(kind, right_side), m_term(t) { } - mpq get_free_coeff_of_left_side() const override { return zero_of_type();} + mpq get_free_coeff_of_left_side() const override { return m_term->m_v;} }; diff --git a/src/util/lp/lar_solver.cpp b/src/util/lp/lar_solver.cpp index bc4c17d65..1340d1826 100644 --- a/src/util/lp/lar_solver.cpp +++ b/src/util/lp/lar_solver.cpp @@ -137,6 +137,7 @@ bool lar_solver::implied_bound_is_correctly_explained(implied_bound const & be, 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; @@ -612,6 +613,7 @@ void lar_solver::substitute_terms_in_linear_expression(const vector(); + value = t.m_v; for (auto const& cv : t) { impq const& r = get_column_value(cv.var()); if (!numeric_traits::is_zero(r.y)) return false; @@ -1495,7 +1497,7 @@ bool lar_solver::term_is_int(const lar_term * t) const { for (auto const & p : t->m_coeffs) if (! (column_is_int(p.first) && p.second.is_int())) return false; - return true; + return t->m_v.is_int(); } bool lar_solver::var_is_int(var_index v) const { @@ -1596,8 +1598,9 @@ void lar_solver::add_new_var_to_core_fields_for_mpq(bool register_in_basis) { } -var_index lar_solver::add_term_undecided(const vector> & coeffs) { - push_and_register_term(new lar_term(coeffs)); +var_index lar_solver::add_term_undecided(const vector> & coeffs, + const mpq &m_v) { + push_and_register_term(new lar_term(coeffs, m_v)); return m_terms_start_index + m_terms.size() - 1; } @@ -1640,11 +1643,12 @@ void lar_solver::push_and_register_term(lar_term* t) { } // terms -var_index lar_solver::add_term(const vector> & coeffs) { +var_index lar_solver::add_term(const vector> & coeffs, + const mpq &m_v) { if (strategy_is_undecided()) - return add_term_undecided(coeffs); + return add_term_undecided(coeffs, m_v); - push_and_register_term(new lar_term(coeffs)); + push_and_register_term(new lar_term(coeffs, m_v)); unsigned adjusted_term_index = m_terms.size() - 1; var_index ret = m_terms_start_index + adjusted_term_index; if (use_tableau() && !coeffs.empty()) { @@ -1739,8 +1743,9 @@ void lar_solver::add_var_bound_on_constraint_for_term(var_index j, lconstraint_k // lp_assert(!term_is_int(m_terms[adjusted_term_index]) || right_side.is_int()); unsigned term_j; if (m_var_register.external_is_used(j, term_j)) { + mpq rs = right_side - m_terms[adjusted_term_index]->m_v; m_constraints.push_back(new lar_term_constraint(m_terms[adjusted_term_index], kind, right_side)); - update_column_type_and_bound(term_j, kind, right_side, ci); + update_column_type_and_bound(term_j, kind, rs, ci); } else { add_constraint_from_term_and_create_new_column_row(j, m_terms[adjusted_term_index], kind, right_side); @@ -1751,7 +1756,7 @@ constraint_index lar_solver::add_constraint(const vector> left_side; mpq rs = -right_side_parm; substitute_terms_in_linear_expression(left_side_with_terms, left_side, rs); - unsigned term_index = add_term(left_side); + unsigned term_index = add_term(left_side, zero_of_type()); constraint_index ci = m_constraints.size(); add_var_bound_on_constraint_for_term(term_index, kind_par, -rs, ci); return ci; @@ -1762,7 +1767,7 @@ void lar_solver::add_constraint_from_term_and_create_new_column_row(unsigned ter add_row_from_term_no_constraint(term, term_j); unsigned j = A_r().column_count() - 1; - update_column_type_and_bound(j, kind, right_side, m_constraints.size()); + update_column_type_and_bound(j, kind, right_side - term->m_v, m_constraints.size()); m_constraints.push_back(new lar_term_constraint(term, kind, right_side)); lp_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_costs.size()); } diff --git a/src/util/lp/lar_solver.h b/src/util/lp/lar_solver.h index 2e85513ff..9afb70c72 100644 --- a/src/util/lp/lar_solver.h +++ b/src/util/lp/lar_solver.h @@ -164,11 +164,13 @@ public: // terms - var_index add_term(const vector> & coeffs); + var_index add_term(const vector> & coeffs, + const mpq &m_v); - var_index add_term_undecided(const vector> & coeffs); + var_index add_term_undecided(const vector> & coeffs, + const mpq &m_v); - bool term_coeffs_are_ok(const vector> & coeffs); + bool term_coeffs_are_ok(const vector> & coeffs, const mpq& v); void push_and_register_term(lar_term* t); void add_row_for_term(const lar_term * term, unsigned term_ext_index); diff --git a/src/util/lp/lar_term.h b/src/util/lp/lar_term.h index 47861baf8..519847848 100644 --- a/src/util/lp/lar_term.h +++ b/src/util/lp/lar_term.h @@ -23,6 +23,7 @@ namespace lp { struct lar_term { // the term evaluates to sum of m_coeffs + m_v std::unordered_map m_coeffs; + mpq m_v; lar_term() {} void add_monomial(const mpq& c, unsigned j) { auto it = m_coeffs.find(j); @@ -36,7 +37,7 @@ struct lar_term { } bool is_empty() const { - return m_coeffs.size() == 0; + return m_coeffs.size() == 0 && is_zero(m_v); } unsigned size() const { return static_cast(m_coeffs.size()); } @@ -45,7 +46,8 @@ struct lar_term { return m_coeffs; } - lar_term(const vector>& coeffs) { + lar_term(const vector>& coeffs, + const mpq & v) : m_v(v) { for (const auto & p : coeffs) { add_monomial(p.first, p.second); } @@ -85,7 +87,7 @@ struct lar_term { template T apply(const vector& x) const { - T ret = zero_of_type(); + T ret = T(m_v); for (const auto & t : m_coeffs) { ret += t.second * x[t.first]; } @@ -94,6 +96,7 @@ struct lar_term { void clear() { m_coeffs.clear(); + m_v = zero_of_type(); } struct ival { From 257ba6218fcf41b531b55f7f7d45ae092e189bf2 Mon Sep 17 00:00:00 2001 From: Lev Date: Fri, 14 Sep 2018 11:48:17 -0700 Subject: [PATCH 05/11] remove gomory.h Signed-off-by: Lev --- src/util/lp/gomory.h | 46 -------------------------------------------- 1 file changed, 46 deletions(-) delete mode 100644 src/util/lp/gomory.h diff --git a/src/util/lp/gomory.h b/src/util/lp/gomory.h deleted file mode 100644 index 87879f9f1..000000000 --- a/src/util/lp/gomory.h +++ /dev/null @@ -1,46 +0,0 @@ -/*++ -Copyright (c) 2017 Microsoft Corporation - -Module Name: - - - -Abstract: - - - -Author: - Nikolaj Bjorner (nbjorner) - Lev Nachmanson (levnach) - -Revision History: - - ---*/ -#pragma once -#include "util/lp/lar_term.h" -#include "util/lp/lia_move.h" -#include "util/lp/explanation.h" - -namespace lp { -class gomory { - lar_term & m_t; // the term to return in the cut - mpq & m_k; // the right side of the cut - explanation& m_ex; // the conflict explanation - bool & m_upper; // we have a cut m_t*x <= k if m_upper is true nad m_t*x >= k otherwise - unsigned m_basic_inf_int_j; // a basis column which has to be an integer but has a not integral value - const row_strip& m_row -public : - gomory(lar_term & m_t, - mpq & m_k, - explanation& m_ex, - bool & m_upper, - unsigned m_basic_inf_int_j ) : - m_t(t), - m_k(k), - m_ex(ex), - m_upper(upper), - m_basic_inf_int_j(j) { - } -}; -} From 26764b076f2eab2656a5093946877376999a5530 Mon Sep 17 00:00:00 2001 From: Lev Date: Fri, 14 Sep 2018 12:39:46 -0700 Subject: [PATCH 06/11] adjust cuts and branch (m_t and m_k) for terms Signed-off-by: Lev --- src/util/lp/int_solver.cpp | 18 ++++++++++++------ src/util/lp/int_solver.h | 2 +- src/util/lp/lar_solver.cpp | 10 ++++++++++ src/util/lp/lar_solver.h | 1 + 4 files changed, 24 insertions(+), 7 deletions(-) diff --git a/src/util/lp/int_solver.cpp b/src/util/lp/int_solver.cpp index c416c63f3..f12e93103 100644 --- a/src/util/lp/int_solver.cpp +++ b/src/util/lp/int_solver.cpp @@ -383,19 +383,25 @@ typedef monomial mono; // this will allow to enable and disable tracking of the pivot rows -struct pivoted_rows_tracking_control { - lar_solver * m_lar_solver; - bool m_track_pivoted_rows; - pivoted_rows_tracking_control(lar_solver* ls) : +struct check_return_helper { + lar_solver * m_lar_solver; + const lia_move & m_r; + bool m_track_pivoted_rows; + check_return_helper(lar_solver* ls, const lia_move& r) : m_lar_solver(ls), + m_r(r), m_track_pivoted_rows(ls->get_track_pivoted_rows()) { TRACE("pivoted_rows", tout << "pivoted rows = " << ls->m_mpq_lar_core_solver.m_r_solver.m_pivoted_rows->size() << std::endl;); m_lar_solver->set_track_pivoted_rows(false); } - ~pivoted_rows_tracking_control() { + ~check_return_helper() { TRACE("pivoted_rows", tout << "pivoted rows = " << m_lar_solver->m_mpq_lar_core_solver.m_r_solver.m_pivoted_rows->size() << std::endl;); m_lar_solver->set_track_pivoted_rows(m_track_pivoted_rows); + if (m_r == lia_move::cut || m_r == lia_move::branch) { + int_solver * s = m_lar_solver->get_int_solver(); + m_lar_solver->adjust_cut_for_terms(*(s->m_t), *(s->m_k)); + } } }; @@ -615,7 +621,7 @@ lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex, bool & upper) { lia_move r = run_gcd_test(); if (r != lia_move::undef) return r; - pivoted_rows_tracking_control pc(m_lar_solver); + check_return_helper pc(m_lar_solver, r); if(settings().m_int_pivot_fixed_vars_from_basis) m_lar_solver->pivot_fixed_vars_from_basis(); diff --git a/src/util/lp/int_solver.h b/src/util/lp/int_solver.h index 82fcb6eb4..414ca3006 100644 --- a/src/util/lp/int_solver.h +++ b/src/util/lp/int_solver.h @@ -160,5 +160,5 @@ public: bool hnf_has_var_with_non_integral_value() const; bool hnf_cutter_is_full() const; void patch_nbasic_column(unsigned j, bool patch_only_int_vals); -}; + }; } diff --git a/src/util/lp/lar_solver.cpp b/src/util/lp/lar_solver.cpp index 1340d1826..89044c7d6 100644 --- a/src/util/lp/lar_solver.cpp +++ b/src/util/lp/lar_solver.cpp @@ -2265,6 +2265,16 @@ void lar_solver::set_cut_strategy(unsigned cut_frequency) { } } +void lar_solver::adjust_cut_for_terms(const lar_term& t, mpq & rs) { + for (const auto& p : t) { + if (!is_term(p.var())) continue; + const lar_term & p_term = get_term(p.var()); + if (p_term.m_v.is_zero()) continue; + rs -= p.coeff() * p_term.m_v; + } +} + + } // namespace lp diff --git a/src/util/lp/lar_solver.h b/src/util/lp/lar_solver.h index 9afb70c72..6ef0ea596 100644 --- a/src/util/lp/lar_solver.h +++ b/src/util/lp/lar_solver.h @@ -584,5 +584,6 @@ public: lar_term get_term_to_maximize(unsigned ext_j) const; void set_cut_strategy(unsigned cut_frequency); bool sum_first_coords(const lar_term& t, mpq & val) const; + void adjust_cut_for_terms(const lar_term& t, mpq & rs); }; } From 324396e4039d417057f466341fba45108c57dd24 Mon Sep 17 00:00:00 2001 From: Lev Date: Fri, 14 Sep 2018 17:12:49 -0700 Subject: [PATCH 07/11] separate the gomory cut functionality in a separate file Signed-off-by: Lev --- src/util/lp/CMakeLists.txt | 1 + src/util/lp/gomory.cpp | 227 +++++++++++++++++++++++++++++++++++++ src/util/lp/gomory.h | 36 ++++++ src/util/lp/int_solver.cpp | 221 +----------------------------------- src/util/lp/int_solver.h | 25 ++-- src/util/lp/lar_solver.cpp | 1 + src/util/lp/lar_solver.h | 5 +- 7 files changed, 284 insertions(+), 232 deletions(-) create mode 100644 src/util/lp/gomory.cpp create mode 100644 src/util/lp/gomory.h diff --git a/src/util/lp/CMakeLists.txt b/src/util/lp/CMakeLists.txt index bde6ed93a..539c68712 100644 --- a/src/util/lp/CMakeLists.txt +++ b/src/util/lp/CMakeLists.txt @@ -6,6 +6,7 @@ z3_add_component(lp core_solver_pretty_printer.cpp dense_matrix.cpp eta_matrix.cpp + gomory.cpp indexed_vector.cpp int_solver.cpp lar_solver.cpp diff --git a/src/util/lp/gomory.cpp b/src/util/lp/gomory.cpp new file mode 100644 index 000000000..dd3d7bbed --- /dev/null +++ b/src/util/lp/gomory.cpp @@ -0,0 +1,227 @@ +/*++ + Copyright (c) 2017 Microsoft Corporation + + Module Name: + + + + Abstract: + + + + Author: + Nikolaj Bjorner (nbjorner) + Lev Nachmanson (levnach) + + Revision History: + + + --*/ +#include "util/lp/gomory.h" +#include "util/lp/int_solver.h" +#include "util/lp/lar_solver.h" +namespace lp { + +class gomory::imp { + lar_term & m_t; // the term to return in the cut + mpq & m_k; // the right side of the cut + explanation& m_ex; // the conflict explanation + unsigned m_inf_col; // a basis column which has to be an integer but has a not integral value + const row_strip& m_row; + const int_solver& m_int_solver; + + + const impq & get_value(unsigned j) const { return m_int_solver.get_value(j); } + bool is_real(unsigned j) const { return m_int_solver.is_real(j); } + bool at_lower(unsigned j) const { return m_int_solver.at_lower(j); } + bool at_upper(unsigned j) const { return m_int_solver.at_upper(j); } + const impq & lower_bound(unsigned j) const { return m_int_solver.lower_bound(j); } + const impq & upper_bound(unsigned j) const { return m_int_solver.upper_bound(j); } + constraint_index column_lower_bound_constraint(unsigned j) const { return m_int_solver.column_lower_bound_constraint(j); } + constraint_index column_upper_bound_constraint(unsigned j) const { return m_int_solver.column_upper_bound_constraint(j); } + + void int_case_in_gomory_cut(const mpq & a, unsigned x_j, + mpq & lcm_den, const mpq& f_0, const mpq& one_minus_f_0) { + lp_assert(is_int(x_j)); + lp_assert(!a.is_int()); + mpq f_j = int_solver::fractional_part(a); + TRACE("gomory_cut_detail", + tout << a << " x_j" << x_j << " k = " << m_k << "\n"; + tout << "f_j: " << f_j << "\n"; + tout << "f_0: " << f_0 << "\n"; + tout << "1 - f_0: " << 1 - f_0 << "\n"; + tout << "at_lower(" << x_j << ") = " << at_lower(x_j) << std::endl; + ); + lp_assert (!f_j.is_zero()); + mpq new_a; + if (at_lower(x_j)) { + if (f_j <= one_minus_f_0) { + new_a = f_j / one_minus_f_0; + } + else { + new_a = (1 - f_j) / f_0; + } + m_k.addmul(new_a, lower_bound(x_j).x); + m_ex.push_justification(column_lower_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) / one_minus_f_0; + } + new_a.neg(); // the upper terms are inverted + m_k.addmul(new_a, upper_bound(x_j).x); + m_ex.push_justification(column_upper_bound_constraint(x_j), new_a); + } + TRACE("gomory_cut_detail", tout << "new_a: " << new_a << " k: " << m_k << "\n";); + m_t.add_monomial(new_a, x_j); + lcm_den = lcm(lcm_den, denominator(new_a)); + } + + void real_case_in_gomory_cut(const mpq & a, unsigned x_j, const mpq& f_0, const mpq& one_minus_f_0) { + TRACE("gomory_cut_detail_real", tout << "real\n";); + mpq new_a; + if (at_lower(x_j)) { + if (a.is_pos()) { + new_a = a / one_minus_f_0; + } + else { + new_a = a / f_0; + new_a.neg(); + } + m_k.addmul(new_a, lower_bound(x_j).x); // is it a faster operation than + // k += lower_bound(x_j).x * new_a; + m_ex.push_justification(column_lower_bound_constraint(x_j), new_a); + } + else { + lp_assert(at_upper(x_j)); + if (a.is_pos()) { + new_a = a / f_0; + new_a.neg(); // the upper terms are inverted. + } + else { + new_a = a / one_minus_f_0; + } + m_k.addmul(new_a, upper_bound(x_j).x); // k += upper_bound(x_j).x * new_a; + m_ex.push_justification(column_upper_bound_constraint(x_j), new_a); + } + TRACE("gomory_cut_detail_real", tout << a << "*v" << x_j << " k: " << m_k << "\n";); + m_t.add_monomial(new_a, x_j); + } + + lia_move report_conflict_from_gomory_cut() { + lp_assert(m_k.is_pos()); + // conflict 0 >= k where k is positive + m_k.neg(); // returning 0 <= -k + return lia_move::conflict; + } + + void adjust_term_and_k_for_some_ints_case_gomory(mpq &lcm_den) { + lp_assert(!m_t.is_empty()); + auto pol = m_t.coeffs_as_vector(); + m_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; + m_k /= a; + if (a.is_pos()) { // we have av >= k + if (!m_k.is_int()) + m_k = ceil(m_k); + // switch size + m_t.add_monomial(- mpq(1), v); + m_k.neg(); + } else { + if (!m_k.is_int()) + m_k = floor(m_k); + m_t.add_monomial(mpq(1), v); + } + } else { + TRACE("gomory_cut_detail", tout << "pol.size() > 1" << std::endl;); + lcm_den = lcm(lcm_den, denominator(m_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()); + } + m_k *= lcm_den; + } + // negate everything to return -pol <= -m_k + for (const auto & pi: pol) + m_t.add_monomial(-pi.first, pi.second); + m_k.neg(); + } + TRACE("gomory_cut_detail", tout << "k = " << m_k << std::endl;); + lp_assert(m_k.is_int()); + } +public: + lia_move create_cut() { + TRACE("gomory_cut", + tout << "applying cut at:\n"; m_int_solver.m_lar_solver->print_row(m_row, tout); tout << std::endl; + for (auto & p : m_row) { + m_int_solver.m_lar_solver->m_mpq_lar_core_solver.m_r_solver.print_column_info(p.var(), tout); + } + tout << "inf_col = " << m_inf_col << std::endl; + ); + + // gomory will be t <= k and the current solution has a property t > k + m_k = 1; + mpq lcm_den(1); + unsigned x_j; + mpq a; + bool some_int_columns = false; + mpq f_0 = int_solver::fractional_part(get_value(m_inf_col)); + mpq one_min_f_0 = 1 - f_0; + for (const auto & p : m_row) { + x_j = p.var(); + if (x_j == m_inf_col) + continue; + // make the format compatible with the format used in: Integrating Simplex with DPLL(T) + a = p.coeff(); + a.neg(); + if (is_real(x_j)) + real_case_in_gomory_cut(a, x_j, f_0, one_min_f_0); + else if (!a.is_int()) { // f_j will be zero and no monomial will be added + some_int_columns = true; + int_case_in_gomory_cut(a, x_j, lcm_den, f_0, one_min_f_0); + } + } + + if (m_t.is_empty()) + return report_conflict_from_gomory_cut(); + if (some_int_columns) + adjust_term_and_k_for_some_ints_case_gomory(lcm_den); + lp_assert(m_int_solver.current_solution_is_inf_on_cut()); + m_int_solver.m_lar_solver->subs_term_columns(m_t, m_k); + TRACE("gomory_cut", tout<<"gomory cut:"; m_int_solver.m_lar_solver->print_term(m_t, tout); tout << " <= " << m_k << std::endl;); + return lia_move::cut; + } + imp(lar_term & t, mpq & k, explanation& ex, unsigned basic_inf_int_j, const row_strip& row, const int_solver& int_slv ) : + m_t(t), + m_k(k), + m_ex(ex), + m_inf_col(basic_inf_int_j), + m_row(row), + m_int_solver(int_slv) + { + } + +}; + +lia_move gomory::create_cut() { + return m_imp->create_cut(); +} + +gomory::gomory(lar_term & t, mpq & k, explanation& ex, unsigned basic_inf_int_j, const row_strip& row, const int_solver& s) { + m_imp = alloc(imp, t, k, ex, basic_inf_int_j, row, s); +} + +gomory::~gomory() { dealloc(m_imp); } + +} diff --git a/src/util/lp/gomory.h b/src/util/lp/gomory.h new file mode 100644 index 000000000..b7946d6b0 --- /dev/null +++ b/src/util/lp/gomory.h @@ -0,0 +1,36 @@ +/*++ +Copyright (c) 2017 Microsoft Corporation + +Module Name: + + + +Abstract: + + + +Author: + Nikolaj Bjorner (nbjorner) + Lev Nachmanson (levnach) + +Revision History: + + +--*/ +#pragma once +#include "util/lp/lar_term.h" +#include "util/lp/lia_move.h" +#include "util/lp/explanation.h" +#include "util/lp/static_matrix.h" + +namespace lp { +class int_solver; +class gomory { + class imp; + imp *m_imp; +public : + gomory(lar_term & t, mpq & k, explanation& ex, unsigned basic_inf_int_j, const row_strip& row, const int_solver& s); + lia_move create_cut(); + ~gomory(); +}; +} diff --git a/src/util/lp/int_solver.cpp b/src/util/lp/int_solver.cpp index f12e93103..cc01a0038 100644 --- a/src/util/lp/int_solver.cpp +++ b/src/util/lp/int_solver.cpp @@ -8,6 +8,7 @@ #include "util/lp/lp_utils.h" #include #include "util/lp/monomial.h" +#include "util/lp/gomory.h" namespace lp { @@ -101,12 +102,7 @@ bool int_solver::is_gomory_cut_target(const row_strip& row) { unsigned j; for (const auto & p : row) { j = p.var(); - if (is_base(j)) continue; - if (is_free(j)) - return false; - if (!at_bound(j)) - return false; - if (!is_zero(get_value(j).y)) { + if (!is_base(j) && (!at_bound(j) || !is_zero(get_value(j).y))) { TRACE("gomory_cut", tout << "row is not gomory cut target:\n"; display_column(tout, j); tout << "infinitesimal: " << !is_zero(get_value(j).y) << "\n";); @@ -117,36 +113,6 @@ bool int_solver::is_gomory_cut_target(const row_strip& row) { } -void int_solver::real_case_in_gomory_cut(const mpq & a, unsigned x_j, const mpq& f_0, const mpq& one_minus_f_0) { - TRACE("gomory_cut_detail_real", tout << "real\n";); - mpq new_a; - if (at_low(x_j)) { - if (a.is_pos()) { - new_a = a / one_minus_f_0; - } - else { - new_a = a / f_0; - new_a.neg(); - } - m_k->addmul(new_a, lower_bound(x_j).x); // is it a faster operation than - // k += lower_bound(x_j).x * new_a; - m_ex->push_justification(column_lower_bound_constraint(x_j), new_a); - } - else { - lp_assert(at_upper(x_j)); - if (a.is_pos()) { - new_a = a / f_0; - new_a.neg(); // the upper terms are inverted. - } - else { - new_a = a / one_minus_f_0; - } - m_k->addmul(new_a, upper_bound(x_j).x); // k += upper_bound(x_j).x * new_a; - m_ex->push_justification(column_upper_bound_constraint(x_j), new_a); - } - TRACE("gomory_cut_detail_real", tout << a << "*v" << x_j << " k: " << *m_k << "\n";); - m_t->add_monomial(new_a, x_j); -} constraint_index int_solver::column_upper_bound_constraint(unsigned j) const { return m_lar_solver->get_column_upper_bound_witness(j); @@ -157,99 +123,6 @@ constraint_index int_solver::column_lower_bound_constraint(unsigned j) const { } -void int_solver::int_case_in_gomory_cut(const mpq & a, unsigned x_j, - mpq & lcm_den, const mpq& f_0, const mpq& one_minus_f_0) { - lp_assert(is_int(x_j)); - lp_assert(!a.is_int()); - mpq f_j = fractional_part(a); - TRACE("gomory_cut_detail", - tout << a << " x_j" << x_j << " k = " << *m_k << "\n"; - tout << "f_j: " << f_j << "\n"; - tout << "f_0: " << f_0 << "\n"; - 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)) { - if (f_j <= one_minus_f_0) { - new_a = f_j / one_minus_f_0; - } - else { - new_a = (1 - f_j) / f_0; - } - m_k->addmul(new_a, lower_bound(x_j).x); - m_ex->push_justification(column_lower_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) / one_minus_f_0; - } - new_a.neg(); // the upper terms are inverted - m_k->addmul(new_a, upper_bound(x_j).x); - m_ex->push_justification(column_upper_bound_constraint(x_j), new_a); - } - TRACE("gomory_cut_detail", tout << "new_a: " << new_a << " k: " << *m_k << "\n";); - m_t->add_monomial(new_a, x_j); - lcm_den = lcm(lcm_den, denominator(new_a)); -} - -lia_move int_solver::report_conflict_from_gomory_cut() { - TRACE("empty_pol",); - lp_assert(m_k->is_pos()); - // conflict 0 >= k where k is positive - m_k->neg(); // returning 0 <= -k - return lia_move::conflict; -} - -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; - } - // negate everything to return -pol <= -k - for (const auto & pi: pol) - t.add_monomial(-pi.first, pi.second); - k.neg(); - } -} - bool int_solver::current_solution_is_inf_on_cut() const { const auto & x = m_lar_solver->m_mpq_lar_core_solver.m_r_x; impq v = m_t->apply(x); @@ -261,95 +134,11 @@ bool int_solver::current_solution_is_inf_on_cut() const { return v * sign > (*m_k) * sign; } -void int_solver::adjust_term_and_k_for_some_ints_case_gomory(mpq &lcm_den) { - lp_assert(!m_t->is_empty()); - auto pol = m_t->coeffs_as_vector(); - m_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; - (*m_k) /= a; - if (a.is_pos()) { // we have av >= k - if (!(*m_k).is_int()) - (*m_k) = ceil((*m_k)); - // switch size - m_t->add_monomial(- mpq(1), v); - (*m_k).neg(); - } else { - if (!(*m_k).is_int()) - (*m_k) = floor((*m_k)); - m_t->add_monomial(mpq(1), v); - } - } else { - TRACE("gomory_cut_detail", tout << "pol.size() > 1" << std::endl;); - lcm_den = lcm(lcm_den, denominator((*m_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()); - } - (*m_k) *= lcm_den; - } - // negate everything to return -pol <= -(*m_k) - for (const auto & pi: pol) - m_t->add_monomial(-pi.first, pi.second); - (*m_k).neg(); - } - TRACE("gomory_cut_detail", tout << "k = " << (*m_k) << std::endl;); - lp_assert((*m_k).is_int()); -} - - - - lia_move int_solver::mk_gomory_cut( unsigned inf_col, const row_strip & row) { - lp_assert(column_is_int_inf(inf_col)); - TRACE("gomory_cut", - tout << "applying cut at:\n"; m_lar_solver->print_row(row, tout); tout << std::endl; - for (auto & p : row) { - m_lar_solver->m_mpq_lar_core_solver.m_r_solver.print_column_info(p.var(), tout); - } - tout << "inf_col = " << inf_col << std::endl; - ); - - // gomory will be t <= k and the current solution has a property t > k - *m_k = 1; - mpq lcm_den(1); - unsigned x_j; - mpq a; - bool some_int_columns = false; - mpq f_0 = int_solver::fractional_part(get_value(inf_col)); - mpq one_min_f_0 = 1 - f_0; - for (const auto & p : row) { - x_j = p.var(); - if (x_j == inf_col) - continue; - // make the format compatible with the format used in: Integrating Simplex with DPLL(T) - a = p.coeff(); - a.neg(); - if (is_real(x_j)) - real_case_in_gomory_cut(a, x_j, f_0, one_min_f_0); - else if (!a.is_int()) { // f_j will be zero and no monomial will be added - some_int_columns = true; - int_case_in_gomory_cut(a, x_j, lcm_den, f_0, one_min_f_0); - } - } - - if (m_t->is_empty()) - return report_conflict_from_gomory_cut(); - if (some_int_columns) - adjust_term_and_k_for_some_ints_case_gomory(lcm_den); - - lp_assert(current_solution_is_inf_on_cut()); - m_lar_solver->subs_term_columns(*m_t); - TRACE("gomory_cut", tout<<"precut:"; m_lar_solver->print_term(*m_t, tout); tout << " <= " << *m_k << std::endl;); - return lia_move::cut; + gomory gc(*m_t, *m_k, *m_ex, inf_col, row, *this); + return gc.create_cut(); } lia_move int_solver::proceed_with_gomory_cut(unsigned j) { @@ -1121,7 +910,7 @@ bool int_solver::at_bound(unsigned j) const { } } -bool int_solver::at_low(unsigned j) const { +bool int_solver::at_lower(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: diff --git a/src/util/lp/int_solver.h b/src/util/lp/int_solver.h index 414ca3006..dfe51711c 100644 --- a/src/util/lp/int_solver.h +++ b/src/util/lp/int_solver.h @@ -53,6 +53,13 @@ public: bool move_non_basic_column_to_bounds(unsigned j); lia_move check_wrapper(lar_term& t, mpq& k, explanation& ex); bool is_base(unsigned j) const; + bool is_real(unsigned j) const; + const impq & lower_bound(unsigned j) const; + const impq & upper_bound(unsigned j) const; + bool is_int(unsigned j) const; + const impq & get_value(unsigned j) const; + bool at_lower(unsigned j) const; + bool at_upper(unsigned j) const; private: @@ -79,10 +86,7 @@ private: void add_to_explanation_from_fixed_or_boxed_column(unsigned j); lia_move patch_nbasic_columns(); bool get_freedom_interval_for_column(unsigned j, bool & inf_l, impq & l, bool & inf_u, impq & u, mpq & m); - const impq & lower_bound(unsigned j) const; - const impq & upper_bound(unsigned j) const; - bool is_int(unsigned j) const; - bool is_real(unsigned j) const; +private: bool is_boxed(unsigned j) const; bool is_fixed(unsigned j) const; bool is_free(unsigned j) const; @@ -91,7 +95,6 @@ private: void set_value_for_nbasic_column_ignore_old_values(unsigned j, const impq & new_val); bool non_basic_columns_are_at_bounds() const; bool is_feasible() const; - const impq & get_value(unsigned j) const; bool column_is_int_inf(unsigned j) const; void trace_inf_rows() const; lia_move branch_or_sat(); @@ -104,13 +107,9 @@ private: bool move_non_basic_columns_to_bounds(); void branch_infeasible_int_var(unsigned); lia_move mk_gomory_cut(unsigned inf_col, const row_strip& row); - lia_move report_conflict_from_gomory_cut(); - void adjust_term_and_k_for_some_ints_case_gomory(mpq& lcm_den); lia_move proceed_with_gomory_cut(unsigned j); bool is_gomory_cut_target(const row_strip&); bool at_bound(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; @@ -125,17 +124,13 @@ public: lp_assert(is_rational(n)); return n.x - floor(n.x); } -private: - void real_case_in_gomory_cut(const mpq & a, unsigned x_j, const mpq& f_0, const mpq& one_minus_f_0); - void int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & lcm_den, const mpq& f_0, const mpq& one_minus_f_0); constraint_index column_upper_bound_constraint(unsigned j) const; constraint_index column_lower_bound_constraint(unsigned j) const; - void display_row_info(std::ostream & out, unsigned row_index) const; - 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; -public: + bool shift_var(unsigned j, unsigned range); private: + void display_row_info(std::ostream & out, unsigned row_index) const; unsigned random(); bool has_inf_int() const; lia_move create_branch_on_column(int j); diff --git a/src/util/lp/lar_solver.cpp b/src/util/lp/lar_solver.cpp index 89044c7d6..3a53b6068 100644 --- a/src/util/lp/lar_solver.cpp +++ b/src/util/lp/lar_solver.cpp @@ -1645,6 +1645,7 @@ void lar_solver::push_and_register_term(lar_term* t) { // terms var_index lar_solver::add_term(const vector> & coeffs, const mpq &m_v) { + TRACE("add_term_lar_solver", print_linear_combination_of_column_indices(coeffs, tout);); if (strategy_is_undecided()) return add_term_undecided(coeffs, m_v); diff --git a/src/util/lp/lar_solver.h b/src/util/lp/lar_solver.h index 6ef0ea596..3c0ed4fbf 100644 --- a/src/util/lp/lar_solver.h +++ b/src/util/lp/lar_solver.h @@ -535,7 +535,7 @@ public: return m_columns_to_ul_pairs()[j].lower_bound_witness(); } - void subs_term_columns(lar_term& t) { + void subs_term_columns(lar_term& t, mpq & rs) { vector> columns_to_subs; for (const auto & m : t.m_coeffs) { unsigned tj = adjust_column_index_to_term_index(m.first); @@ -545,9 +545,12 @@ public: for (const auto & p : columns_to_subs) { auto it = t.m_coeffs.find(p.first); lp_assert(it != t.m_coeffs.end()); + const lar_term& lt = get_term(p.second); mpq v = it->second; t.m_coeffs.erase(it); t.m_coeffs[p.second] = v; + if (lt.m_v.is_zero()) continue; + rs -= v * lt.m_v; } } From 03d55426bb6c64c248d5024f1b423b2aa26964e5 Mon Sep 17 00:00:00 2001 From: Lev Date: Sat, 15 Sep 2018 17:15:46 -0700 Subject: [PATCH 08/11] fixes in gomory cut Signed-off-by: Lev --- src/smt/theory_lra.cpp | 3 +- src/util/lp/gomory.cpp | 93 ++++++++++++++++++++------------------ src/util/lp/lar_solver.cpp | 2 +- 3 files changed, 53 insertions(+), 45 deletions(-) diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp index f8e2f9fe1..d976df56a 100644 --- a/src/smt/theory_lra.cpp +++ b/src/smt/theory_lra.cpp @@ -1704,10 +1704,11 @@ public: TRACE("arith", tout << "canceled\n";); return l_undef; } + /* if (!check_idiv_bounds()) { TRACE("arith", tout << "idiv bounds check\n";); return l_false; - } + }*/ lp::lar_term term; lp::mpq k; lp::explanation ex; // TBD, this should be streamlined accross different explanations diff --git a/src/util/lp/gomory.cpp b/src/util/lp/gomory.cpp index dd3d7bbed..ec6877536 100644 --- a/src/util/lp/gomory.cpp +++ b/src/util/lp/gomory.cpp @@ -39,57 +39,59 @@ class gomory::imp { const impq & upper_bound(unsigned j) const { return m_int_solver.upper_bound(j); } constraint_index column_lower_bound_constraint(unsigned j) const { return m_int_solver.column_lower_bound_constraint(j); } constraint_index column_upper_bound_constraint(unsigned j) const { return m_int_solver.column_upper_bound_constraint(j); } - - void int_case_in_gomory_cut(const mpq & a, unsigned x_j, - mpq & lcm_den, const mpq& f_0, const mpq& one_minus_f_0) { - lp_assert(is_int(x_j)); - lp_assert(!a.is_int()); - mpq f_j = int_solver::fractional_part(a); + bool column_is_fixed(unsigned j) const { return m_int_solver.m_lar_solver->column_is_fixed(j); } + void int_case_in_gomory_cut(const mpq & a, unsigned j, + mpq & lcm_den, const mpq& f0, const mpq& one_minus_f0) { + lp_assert(is_int(j) && !a.is_int()); + mpq fj = int_solver::fractional_part(a); + lp_assert(fj.is_pos()); TRACE("gomory_cut_detail", - tout << a << " x_j" << x_j << " k = " << m_k << "\n"; - tout << "f_j: " << f_j << "\n"; - tout << "f_0: " << f_0 << "\n"; - tout << "1 - f_0: " << 1 - f_0 << "\n"; - tout << "at_lower(" << x_j << ") = " << at_lower(x_j) << std::endl; + tout << a << " j=" << j << " k = " << m_k; + tout << ", fj: " << fj << ", "; + tout << "f0: " << f0 << ", "; + tout << "1 - f0: " << 1 - f0 << ", "; + tout << (at_lower(j)?"at_lower":"at_upper")<< std::endl; ); - lp_assert (!f_j.is_zero()); mpq new_a; - if (at_lower(x_j)) { - if (f_j <= one_minus_f_0) { - new_a = f_j / one_minus_f_0; + mpq one_minus_fj = 1 - fj; + if (at_lower(j)) { + bool go_for_pos_a = fj / one_minus_f0 < one_minus_fj / f0; + if (go_for_pos_a) { + new_a = fj / one_minus_f0; } else { - new_a = (1 - f_j) / f_0; + new_a = one_minus_fj / f0; } - m_k.addmul(new_a, lower_bound(x_j).x); - m_ex.push_justification(column_lower_bound_constraint(x_j), new_a); + m_k.addmul(new_a, lower_bound(j).x); + m_ex.push_justification(column_lower_bound_constraint(j), new_a); } else { - lp_assert(at_upper(x_j)); - if (f_j <= f_0) { - new_a = f_j / f_0; + bool go_for_pos_a = fj / f0 < one_minus_fj / one_minus_f0; + lp_assert(at_upper(j)); + // the upper terms are inverted + if (go_for_pos_a) { + new_a = - fj / f0; } else { - new_a = (mpq(1) - f_j) / one_minus_f_0; + new_a = - one_minus_fj / one_minus_f0; } - new_a.neg(); // the upper terms are inverted - m_k.addmul(new_a, upper_bound(x_j).x); - m_ex.push_justification(column_upper_bound_constraint(x_j), new_a); + m_k.addmul(new_a, upper_bound(j).x); + m_ex.push_justification(column_upper_bound_constraint(j), new_a); } TRACE("gomory_cut_detail", tout << "new_a: " << new_a << " k: " << m_k << "\n";); - m_t.add_monomial(new_a, x_j); + m_t.add_monomial(new_a, j); lcm_den = lcm(lcm_den, denominator(new_a)); } - void real_case_in_gomory_cut(const mpq & a, unsigned x_j, const mpq& f_0, const mpq& one_minus_f_0) { + void real_case_in_gomory_cut(const mpq & a, unsigned x_j, const mpq& f0, const mpq& one_minus_f0) { TRACE("gomory_cut_detail_real", tout << "real\n";); mpq new_a; if (at_lower(x_j)) { if (a.is_pos()) { - new_a = a / one_minus_f_0; + new_a = a / one_minus_f0; } else { - new_a = a / f_0; + new_a = a / f0; new_a.neg(); } m_k.addmul(new_a, lower_bound(x_j).x); // is it a faster operation than @@ -99,11 +101,11 @@ class gomory::imp { else { lp_assert(at_upper(x_j)); if (a.is_pos()) { - new_a = a / f_0; + new_a = a / f0; new_a.neg(); // the upper terms are inverted. } else { - new_a = a / one_minus_f_0; + new_a = a / one_minus_f0; } m_k.addmul(new_a, upper_bound(x_j).x); // k += upper_bound(x_j).x * new_a; m_ex.push_justification(column_upper_bound_constraint(x_j), new_a); @@ -173,23 +175,28 @@ public: // gomory will be t <= k and the current solution has a property t > k m_k = 1; mpq lcm_den(1); - unsigned x_j; - mpq a; bool some_int_columns = false; - mpq f_0 = int_solver::fractional_part(get_value(m_inf_col)); - mpq one_min_f_0 = 1 - f_0; + mpq f0 = int_solver::fractional_part(get_value(m_inf_col)); + mpq one_min_f0 = 1 - f0; for (const auto & p : m_row) { - x_j = p.var(); - if (x_j == m_inf_col) + unsigned j = p.var(); + if (column_is_fixed(j)) { + m_ex.push_justification(column_lower_bound_constraint(j)); + m_ex.push_justification(column_upper_bound_constraint(j)); continue; + } + if (j == m_inf_col) { + lp_assert(p.coeff() == one_of_type()); + TRACE("gomory_cut_detail", tout << "seeing basic var";); + continue; + } // make the format compatible with the format used in: Integrating Simplex with DPLL(T) - a = p.coeff(); - a.neg(); - if (is_real(x_j)) - real_case_in_gomory_cut(a, x_j, f_0, one_min_f_0); - else if (!a.is_int()) { // f_j will be zero and no monomial will be added + mpq a = - p.coeff(); + if (is_real(j)) + real_case_in_gomory_cut(a, j, f0, one_min_f0); + else if (!a.is_int()) { // fj will be zero and no monomial will be added some_int_columns = true; - int_case_in_gomory_cut(a, x_j, lcm_den, f_0, one_min_f_0); + int_case_in_gomory_cut(a, j, lcm_den, f0, one_min_f0); } } diff --git a/src/util/lp/lar_solver.cpp b/src/util/lp/lar_solver.cpp index 3a53b6068..09c3bd5b9 100644 --- a/src/util/lp/lar_solver.cpp +++ b/src/util/lp/lar_solver.cpp @@ -1645,7 +1645,6 @@ void lar_solver::push_and_register_term(lar_term* t) { // terms var_index lar_solver::add_term(const vector> & coeffs, const mpq &m_v) { - TRACE("add_term_lar_solver", print_linear_combination_of_column_indices(coeffs, tout);); if (strategy_is_undecided()) return add_term_undecided(coeffs, m_v); @@ -1657,6 +1656,7 @@ var_index lar_solver::add_term(const vector> & coeffs, if (m_settings.bound_propagation()) m_rows_with_changed_bounds.insert(A_r().row_count() - 1); } + CTRACE("add_term_lar_solver", !m_v.is_zero(), print_term(*m_terms.back(), tout);); lp_assert(m_var_register.size() == A_r().column_count()); return ret; } From 8c122ba9bd6840d08453b33a3ec5f9c388c3883f Mon Sep 17 00:00:00 2001 From: Lev Date: Sat, 15 Sep 2018 17:33:35 -0700 Subject: [PATCH 09/11] fixes in gomory cut Signed-off-by: Lev --- src/smt/theory_lra.cpp | 3 +-- src/util/lp/gomory.cpp | 22 +++++++--------------- 2 files changed, 8 insertions(+), 17 deletions(-) diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp index d976df56a..f8e2f9fe1 100644 --- a/src/smt/theory_lra.cpp +++ b/src/smt/theory_lra.cpp @@ -1704,11 +1704,10 @@ public: TRACE("arith", tout << "canceled\n";); return l_undef; } - /* if (!check_idiv_bounds()) { TRACE("arith", tout << "idiv bounds check\n";); return l_false; - }*/ + } lp::lar_term term; lp::mpq k; lp::explanation ex; // TBD, this should be streamlined accross different explanations diff --git a/src/util/lp/gomory.cpp b/src/util/lp/gomory.cpp index ec6877536..244244da0 100644 --- a/src/util/lp/gomory.cpp +++ b/src/util/lp/gomory.cpp @@ -55,26 +55,18 @@ class gomory::imp { mpq new_a; mpq one_minus_fj = 1 - fj; if (at_lower(j)) { - bool go_for_pos_a = fj / one_minus_f0 < one_minus_fj / f0; - if (go_for_pos_a) { - new_a = fj / one_minus_f0; - } - else { - new_a = one_minus_fj / f0; - } + mpq fj_over_one_min_f0 = fj / one_minus_f0; + mpq one_minus_fj_over_f0 = one_minus_fj / f0; + new_a = fj_over_one_min_f0 < one_minus_fj_over_f0? fj_over_one_min_f0 : one_minus_fj_over_f0; m_k.addmul(new_a, lower_bound(j).x); m_ex.push_justification(column_lower_bound_constraint(j), new_a); } else { - bool go_for_pos_a = fj / f0 < one_minus_fj / one_minus_f0; + mpq fj_over_f0 = fj / f0; + mpq one_minus_fj_over_one_minus_f0 = one_minus_fj / one_minus_f0; lp_assert(at_upper(j)); - // the upper terms are inverted - if (go_for_pos_a) { - new_a = - fj / f0; - } - else { - new_a = - one_minus_fj / one_minus_f0; - } + // the upper terms are inverted - therefore we have - + new_a = - (fj_over_f0 < one_minus_fj_over_one_minus_f0? fj_over_f0 : one_minus_fj_over_one_minus_f0); m_k.addmul(new_a, upper_bound(j).x); m_ex.push_justification(column_upper_bound_constraint(j), new_a); } From 34bdea750c77639c1d64e74731da8a3292319c29 Mon Sep 17 00:00:00 2001 From: Lev Date: Sat, 15 Sep 2018 17:46:16 -0700 Subject: [PATCH 10/11] fixes in gomory cut Signed-off-by: Lev --- src/util/lp/gomory.cpp | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/src/util/lp/gomory.cpp b/src/util/lp/gomory.cpp index 244244da0..29dc98bd6 100644 --- a/src/util/lp/gomory.cpp +++ b/src/util/lp/gomory.cpp @@ -55,18 +55,14 @@ class gomory::imp { mpq new_a; mpq one_minus_fj = 1 - fj; if (at_lower(j)) { - mpq fj_over_one_min_f0 = fj / one_minus_f0; - mpq one_minus_fj_over_f0 = one_minus_fj / f0; - new_a = fj_over_one_min_f0 < one_minus_fj_over_f0? fj_over_one_min_f0 : one_minus_fj_over_f0; + new_a = fj < one_minus_f0? fj / one_minus_f0 : one_minus_fj / f0; m_k.addmul(new_a, lower_bound(j).x); m_ex.push_justification(column_lower_bound_constraint(j), new_a); } else { - mpq fj_over_f0 = fj / f0; - mpq one_minus_fj_over_one_minus_f0 = one_minus_fj / one_minus_f0; lp_assert(at_upper(j)); // the upper terms are inverted - therefore we have - - new_a = - (fj_over_f0 < one_minus_fj_over_one_minus_f0? fj_over_f0 : one_minus_fj_over_one_minus_f0); + new_a = - (fj < f0? fj / f0 : one_minus_fj / one_minus_f0); m_k.addmul(new_a, upper_bound(j).x); m_ex.push_justification(column_upper_bound_constraint(j), new_a); } From 106b677201fd385a60aa510b6abe827bd0f2ee91 Mon Sep 17 00:00:00 2001 From: Lev Date: Sat, 15 Sep 2018 17:47:54 -0700 Subject: [PATCH 11/11] fixes in gomory cut Signed-off-by: Lev --- src/util/lp/gomory.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/util/lp/gomory.cpp b/src/util/lp/gomory.cpp index 29dc98bd6..55098dccf 100644 --- a/src/util/lp/gomory.cpp +++ b/src/util/lp/gomory.cpp @@ -61,7 +61,7 @@ class gomory::imp { } else { lp_assert(at_upper(j)); - // the upper terms are inverted - therefore we have - + // the upper terms are inverted: therefore we have the minus new_a = - (fj < f0? fj / f0 : one_minus_fj / one_minus_f0); m_k.addmul(new_a, upper_bound(j).x); m_ex.push_justification(column_upper_bound_constraint(j), new_a);