From 8d293171d5b5928f8931dcc10b5c61ca1c5e454e Mon Sep 17 00:00:00 2001 From: Nikolaj Bjorner Date: Sat, 8 Feb 2020 15:27:11 -0800 Subject: [PATCH] separate int-cube functionalty Signed-off-by: Nikolaj Bjorner --- src/math/lp/CMakeLists.txt | 1 + src/math/lp/int_cube.cpp | 117 ++++++++++++++++++++++++++++ src/math/lp/int_cube.h | 41 ++++++++++ src/math/lp/int_solver.cpp | 155 ++++++++----------------------------- src/math/lp/int_solver.h | 1 + 5 files changed, 191 insertions(+), 124 deletions(-) create mode 100644 src/math/lp/int_cube.cpp create mode 100644 src/math/lp/int_cube.h diff --git a/src/math/lp/CMakeLists.txt b/src/math/lp/CMakeLists.txt index 549c16b21..05665d606 100644 --- a/src/math/lp/CMakeLists.txt +++ b/src/math/lp/CMakeLists.txt @@ -13,6 +13,7 @@ z3_add_component(lp horner.cpp indexed_vector.cpp int_solver.cpp + int_cube.cpp lar_solver.cpp lar_core_solver.cpp lp_core_solver_base.cpp diff --git a/src/math/lp/int_cube.cpp b/src/math/lp/int_cube.cpp new file mode 100644 index 000000000..b3aba8cd5 --- /dev/null +++ b/src/math/lp/int_cube.cpp @@ -0,0 +1,117 @@ +/*++ +Copyright (c) 2017 Microsoft Corporation + +Module Name: + + int_cube.cpp + +Abstract: + + Cube finder + +Author: + Nikolaj Bjorner (nbjorner) + Lev Nachmanson (levnach) + +Revision History: +--*/ +#pragma once + + +#include "math/lp/int_solver.h" +#include "math/lp/lar_solver.h" +#include "math/lp/int_cube.h" + +namespace lp { + + int_cube::int_cube(int_solver& lia):lia(lia), lra(*lia.m_lar_solver) {} + + lia_move int_cube::operator()() { + lia.settings().stats().m_cube_calls++; + TRACE("cube", + for (unsigned j = 0; j < lra.A_r().column_count(); j++) + display_column(tout, j); + tout << lra.constraints(); + ); + + lra.push(); + if (!tighten_terms_for_cube()) { + lra.pop(); + return lia_move::undef; + } + + lp_status st = lra.find_feasible_solution(); + if (st != lp_status::FEASIBLE && st != lp_status::OPTIMAL) { + TRACE("cube", tout << "cannot find a feasiblie solution";); + lra.pop(); + lra.move_non_basic_columns_to_bounds(); + find_feasible_solution(); + // it can happen that we found an integer solution here + return !lra.r_basis_has_inf_int()? lia_move::sat: lia_move::undef; + } + lra.pop(); + lra.round_to_integer_solution(); + lra.set_status(lp_status::FEASIBLE); + lp_assert(lia.settings().get_cancel_flag() || is_feasible()); + TRACE("cube", tout << "success";); + lia.settings().stats().m_cube_success++; + return lia_move::sat; + } + + bool int_cube::tighten_term_for_cube(unsigned i) { + unsigned ti = i + lra.terms_start_index(); + if (!lra.term_is_used_as_row(ti)) + return true; + const lar_term* t = lra.terms()[i]; + impq delta = get_cube_delta_for_term(*t); + TRACE("cube", lra.print_term_as_indices(*t, tout); tout << ", delta = " << delta;); + if (is_zero(delta)) + return true; + return lra.tighten_term_bounds_by_delta(i, delta); + } + + bool int_cube::tighten_terms_for_cube() { + for (unsigned i = 0; i < lra.terms().size(); i++) + if (!tighten_term_for_cube(i)) { + TRACE("cube", tout << "cannot tighten";); + return false; + } + return true; + } + + void int_cube::find_feasible_solution() { + lra.find_feasible_solution(); + lp_assert(lp_status::OPTIMAL == lra.get_status() || lp_status::FEASIBLE == lra.get_status()); + } + + impq int_cube::get_cube_delta_for_term(const lar_term& t) const { + if (t.size() == 2) { + bool seen_minus = false; + bool seen_plus = false; + for(const auto & p : t) { + if (!lia.column_is_int(p.var())) + goto usual_delta; + const mpq & c = p.coeff(); + if (c == one_of_type()) { + seen_plus = true; + } else if (c == -one_of_type()) { + seen_minus = true; + } else { + goto usual_delta; + } + } + if (seen_minus && seen_plus) + return zero_of_type(); + return impq(0, 1); + } + usual_delta: + mpq delta = zero_of_type(); + for (const auto & p : t) + if (lia.column_is_int(p.var())) + delta += abs(p.coeff()); + + delta *= mpq(1, 2); + return impq(delta); + } + +} diff --git a/src/math/lp/int_cube.h b/src/math/lp/int_cube.h new file mode 100644 index 000000000..18af0cbc1 --- /dev/null +++ b/src/math/lp/int_cube.h @@ -0,0 +1,41 @@ +/*++ +Copyright (c) 2020 Microsoft Corporation + +Module Name: + + int_cube.h + +Abstract: + + Cube finder + + This routine attempts to find a feasible integer solution + by tightnening bounds and running an LRA solver on the + tighter system. + +Author: + Nikolaj Bjorner (nbjorner) + Lev Nachmanson (levnach) + +Revision History: +--*/ +#pragma once + +#include "math/lp/lia_move.h" + +namespace lp { + class int_solver; + class lar_solver; + class int_cube { + class int_solver& lia; + class lar_solver& lra; + bool tighten_term_for_cube(unsigned i); + bool tighten_terms_for_cube(); + void find_feasible_solution(); + impq get_cube_delta_for_term(const lar_term& t) const; + public: + int_cube(int_solver& lia); + ~int_cube() {} + lia_move operator()(); + }; +} diff --git a/src/math/lp/int_solver.cpp b/src/math/lp/int_solver.cpp index cc7d5f240..44a45b86b 100644 --- a/src/math/lp/int_solver.cpp +++ b/src/math/lp/int_solver.cpp @@ -9,6 +9,7 @@ #include #include "math/lp/monic.h" #include "math/lp/gomory.h" +#include "math/lp/int_cube.h" namespace lp { @@ -97,12 +98,6 @@ int int_solver::find_inf_int_boxed_base_column_with_smallest_range(unsigned & in return result; } - - -constraint_index int_solver::column_upper_bound_constraint(unsigned j) const { - return m_lar_solver->get_column_upper_bound_witness(j); -} - 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); @@ -114,15 +109,38 @@ bool int_solver::current_solution_is_inf_on_cut() const { return v * sign > impq(m_k) * sign; } +constraint_index int_solver::column_upper_bound_constraint(unsigned j) const { + return m_lar_solver->get_column_upper_bound_witness(j); +} + constraint_index int_solver::column_lower_bound_constraint(unsigned j) const { return m_lar_solver->get_column_lower_bound_witness(j); } - unsigned int_solver::row_of_basic_column(unsigned j) const { return m_lar_solver->row_of_basic_column(j); } +lp_settings& int_solver::settings() { + return m_lar_solver->settings(); +} + +const lp_settings& int_solver::settings() const { + return m_lar_solver->settings(); +} + +bool int_solver::column_is_int(unsigned j) const { + return m_lar_solver->column_is_int(j); +} + +bool int_solver::is_real(unsigned j) const { + return !column_is_int(j); +} + +bool int_solver::value_is_int(unsigned j) const { + return m_lar_solver->column_value_is_int(j); +} + // this will allow to enable and disable tracking of the pivot rows struct check_return_helper { @@ -141,102 +159,20 @@ struct check_return_helper { } }; - - -impq int_solver::get_cube_delta_for_term(const lar_term& t) const { - if (t.size() == 2) { - bool seen_minus = false; - bool seen_plus = false; - for(const auto & p : t) { - if (!column_is_int(p.var())) - goto usual_delta; - const mpq & c = p.coeff(); - if (c == one_of_type()) { - seen_plus = true; - } else if (c == -one_of_type()) { - seen_minus = true; - } else { - goto usual_delta; - } - } - if (seen_minus && seen_plus) - return zero_of_type(); - return impq(0, 1); - } - usual_delta: - mpq delta = zero_of_type(); - for (const auto & p : t) - if (column_is_int(p.var())) - delta += abs(p.coeff()); - - delta *= mpq(1, 2); - return impq(delta); -} - -bool int_solver::tighten_term_for_cube(unsigned i) { - unsigned ti = i + m_lar_solver->terms_start_index(); - if (!m_lar_solver->term_is_used_as_row(ti)) - return true; - const lar_term* t = m_lar_solver->terms()[i]; - impq delta = get_cube_delta_for_term(*t); - TRACE("cube", m_lar_solver->print_term_as_indices(*t, tout); tout << ", delta = " << delta;); - if (is_zero(delta)) - return true; - return m_lar_solver->tighten_term_bounds_by_delta(i, delta); -} - -bool int_solver::tighten_terms_for_cube() { - for (unsigned i = 0; i < m_lar_solver->terms().size(); i++) - if (!tighten_term_for_cube(i)) { - TRACE("cube", tout << "cannot tighten";); - return false; - } - return true; -} - bool int_solver::should_find_cube() { return m_number_of_calls % settings().m_int_find_cube_period == 0; } lia_move int_solver::find_cube() { - if (!should_find_cube()) - return lia_move::undef; - - settings().stats().m_cube_calls++; - TRACE("cube", - for (unsigned j = 0; j < m_lar_solver->A_r().column_count(); j++) - display_column(tout, j); - tout << m_lar_solver->constraints(); - ); - - m_lar_solver->push(); - if (!tighten_terms_for_cube()) { - m_lar_solver->pop(); + int_cube ic(*this); + if (should_find_cube()) { + return ic(); + } + else { return lia_move::undef; } - - lp_status st = m_lar_solver->find_feasible_solution(); - if (st != lp_status::FEASIBLE && st != lp_status::OPTIMAL) { - TRACE("cube", tout << "cannot find a feasiblie solution";); - m_lar_solver->pop(); - m_lar_solver->move_non_basic_columns_to_bounds(); - find_feasible_solution(); - // it can happen that we found an integer solution here - return !m_lar_solver->r_basis_has_inf_int()? lia_move::sat: lia_move::undef; - } - m_lar_solver->pop(); - m_lar_solver->round_to_integer_solution(); - m_lar_solver->set_status(lp_status::FEASIBLE); - lp_assert(settings().get_cancel_flag() || is_feasible()); - TRACE("cube", tout << "success";); - settings().stats().m_cube_success++; - return lia_move::sat; } -void int_solver::find_feasible_solution() { - m_lar_solver->find_feasible_solution(); - lp_assert(lp_status::OPTIMAL == m_lar_solver->get_status() || lp_status::FEASIBLE == m_lar_solver->get_status()); -} bool int_solver::should_run_gcd_test() { return settings().m_int_run_gcd_test; @@ -287,14 +223,6 @@ bool int_solver::hnf_cutter_is_full() const { m_hnf_cutter.vars().size() >= settings().limit_on_columns_for_hnf_cutter; } -lp_settings& int_solver::settings() { - return m_lar_solver->settings(); -} - -const lp_settings& int_solver::settings() const { - return m_lar_solver->settings(); -} - bool int_solver::hnf_has_var_with_non_integral_value() const { for (unsigned j : m_hnf_cutter.vars()) if (!get_value(j).is_int()) @@ -634,14 +562,6 @@ bool int_solver::ext_gcd_test(const row_strip & row, return true; } -/* -linear_combination_iterator * int_solver::get_column_iterator(unsigned j) { - if (m_lar_solver->use_tableau()) - return new iterator_on_column(m_lar_solver->A_r().m_columns[j], m_lar_solver->A_r()); - return new iterator_on_indexed_vector(m_lar_solver->get_column_in_lu_mode(j)); -} -*/ - int_solver::int_solver(lar_solver* lar_slv) : m_lar_solver(lar_slv), @@ -674,7 +594,6 @@ bool int_solver::has_upper(unsigned j) const { } } - static void set_lower(impq & l, bool & inf_l, impq const & v ) { if (inf_l || v > l) { l = v; @@ -682,7 +601,6 @@ static void set_lower(impq & l, bool & inf_l, impq const & v ) { } } - static void set_upper(impq & u, bool & inf_u, impq const & v) { if (inf_u || v < u) { u = v; @@ -778,17 +696,6 @@ bool int_solver::get_freedom_interval_for_column(unsigned j, bool & inf_l, impq return (inf_l || inf_u || l <= u); } -bool int_solver::column_is_int(unsigned j) const { - return m_lar_solver->column_is_int(j); -} - -bool int_solver::is_real(unsigned j) const { - return !column_is_int(j); -} - -bool int_solver::value_is_int(unsigned j) const { - return m_lar_solver->column_value_is_int(j); -} bool int_solver::is_feasible() const { auto & lcs = m_lar_solver->m_mpq_lar_core_solver; @@ -802,7 +709,7 @@ const impq & int_solver::get_value(unsigned j) const { } std::ostream& int_solver::display_column(std::ostream & out, unsigned j) const { - m_lar_solver->m_mpq_lar_core_solver.m_r_solver.print_column_info(j, out); + return m_lar_solver->m_mpq_lar_core_solver.m_r_solver.print_column_info(j, out); return out; } diff --git a/src/math/lp/int_solver.h b/src/math/lp/int_solver.h index 97aa441ca..a437d4c19 100644 --- a/src/math/lp/int_solver.h +++ b/src/math/lp/int_solver.h @@ -36,6 +36,7 @@ struct lp_constraint; class int_solver { friend class gomory; + friend class int_cube; public: // fields lar_solver *m_lar_solver;