diff --git a/src/math/lp/CMakeLists.txt b/src/math/lp/CMakeLists.txt index 3bd1104c1..7edfe69ca 100644 --- a/src/math/lp/CMakeLists.txt +++ b/src/math/lp/CMakeLists.txt @@ -10,6 +10,7 @@ z3_add_component(lp factorization.cpp factorization_factory_imp.cpp gomory.cpp + hnf_cutter.cpp horner.cpp indexed_vector.cpp int_branch.cpp diff --git a/src/math/lp/hnf_cutter.cpp b/src/math/lp/hnf_cutter.cpp new file mode 100644 index 000000000..73689057f --- /dev/null +++ b/src/math/lp/hnf_cutter.cpp @@ -0,0 +1,275 @@ +/*++ +Copyright (c) 2017 Microsoft Corporation + +Module Name: + + hnf_cutter.cpp + +Author: + Lev Nachmanson (levnach) + +--*/ +#include "math/lp/int_solver.h" +#include "math/lp/lar_solver.h" +#include "math/lp/hnf_cutter.h" + +namespace lp { + + hnf_cutter::hnf_cutter(int_solver& lia): + lia(lia), + lra(lia.lra), + m_settings(lia.settings()), + m_abs_max(zero_of_type()), + m_var_register(0) {} + + bool hnf_cutter::is_full() const { + return + terms_count() >= lia.settings().limit_on_rows_for_hnf_cutter || + vars().size() >= lia.settings().limit_on_columns_for_hnf_cutter; + } + + void hnf_cutter::clear() { + // m_A will be filled from scratch in init_matrix_A + m_var_register.clear(); + m_terms.clear(); + m_terms_upper.clear(); + m_constraints_for_explanation.clear(); + m_right_sides.clear(); + m_abs_max = zero_of_type(); + m_overflow = false; + } + + void hnf_cutter::add_term(const lar_term* t, const mpq &rs, constraint_index ci, bool upper_bound) { + m_terms.push_back(t); + m_terms_upper.push_back(upper_bound); + if (upper_bound) + m_right_sides.push_back(rs); + else + m_right_sides.push_back(-rs); + + m_constraints_for_explanation.push_back(ci); + + for (const auto &p : *t) { + m_var_register.add_var(p.var()); + mpq t = abs(ceil(p.coeff())); + if (t > m_abs_max) + m_abs_max = t; + } + } + + void hnf_cutter::print(std::ostream & out) { + out << "terms = " << m_terms.size() << ", var = " << m_var_register.size() << std::endl; + } + + void hnf_cutter::initialize_row(unsigned i) { + mpq sign = m_terms_upper[i]? one_of_type(): - one_of_type(); + m_A.init_row_from_container(i, * m_terms[i], [this](unsigned j) { return m_var_register.add_var(j);}, sign); + } + + void hnf_cutter::init_matrix_A() { + m_A = general_matrix(terms_count(), vars().size()); + for (unsigned i = 0; i < terms_count(); i++) + initialize_row(i); + } + + // todo: as we need only one row i with non integral b[i] need to optimize later + void hnf_cutter::find_h_minus_1_b(const general_matrix& H, vector & b) { + // the solution will be put into b + for (unsigned i = 0; i < H.row_count() ;i++) { + for (unsigned j = 0; j < i; j++) { + b[i] -= H[i][j]*b[j]; + } + b[i] /= H[i][i]; + // consider return from here if b[i] is not an integer and return i + } + } + + vector hnf_cutter::create_b(const svector & basis_rows) { + if (basis_rows.size() == m_right_sides.size()) + return m_right_sides; + vector b; + for (unsigned i : basis_rows) { + b.push_back(m_right_sides[i]); + } + return b; + } + + int hnf_cutter::find_cut_row_index(const vector & b) { + int ret = -1; + int n = 0; + for (int i = 0; i < static_cast(b.size()); i++) { + if (is_integer(b[i])) continue; + if (n == 0 ) { + lp_assert(ret == -1); + n = 1; + ret = i; + } else { + if (m_settings.random_next() % (++n) == 0) { + ret = i; + } + } + } + return ret; + } + + // fills e_i*H_minus_1 + void hnf_cutter::get_ei_H_minus_1(unsigned i, const general_matrix& H, vector & row) { + // we solve x = ei * H_min_1 + // or x * H = ei + unsigned m = H.row_count(); + for (unsigned k = i + 1; k < m; k++) { + row[k] = zero_of_type(); + } + row[i] = one_of_type() / H[i][i]; + for(int k = i - 1; k >= 0; k--) { + mpq t = zero_of_type(); + for (unsigned l = k + 1; l <= i; l++) { + t += H[l][k]*row[l]; + } + row[k] = -t / H[k][k]; + } + + // // test region + // vector ei(H.row_count(), zero_of_type()); + // ei[i] = one_of_type(); + // vector pr = row * H; + // pr.shrink(ei.size()); + // lp_assert(ei == pr); + // // end test region + + } + + void hnf_cutter::fill_term(const vector & row, lar_term& t) { + for (unsigned j = 0; j < row.size(); j++) { + if (!is_zero(row[j])) + t.add_monomial(row[j], m_var_register.local_to_external(j)); + } + } +#ifdef Z3DEBUG + vector hnf_cutter::transform_to_local_columns(const vector & x) const { + vector ret; + for (unsigned j = 0; j < vars().size(); j++) { + ret.push_back(x[m_var_register.local_to_external(j)].x); + } + return ret; + } +#endif + void hnf_cutter::shrink_explanation(const svector& basis_rows) { + svector new_expl; + for (unsigned i : basis_rows) { + new_expl.push_back(m_constraints_for_explanation[i]); + } + m_constraints_for_explanation = new_expl; + } + + bool hnf_cutter::overflow() const { return m_overflow; } + + lia_move hnf_cutter::create_cut(lar_term& t, mpq& k, explanation* ex, bool & upper, const vector & x0) { + // we suppose that x0 has at least one non integer element + (void)x0; + + init_matrix_A(); + svector basis_rows; + mpq big_number = m_abs_max.expt(3); + mpq d = hnf_calc::determinant_of_rectangular_matrix(m_A, basis_rows, big_number); + + if (d >= big_number) { + return lia_move::undef; + } + + if (m_settings.get_cancel_flag()) { + return lia_move::undef; + } + + if (basis_rows.size() < m_A.row_count()) { + m_A.shrink_to_rank(basis_rows); + shrink_explanation(basis_rows); + } + + hnf h(m_A, d); + vector b = create_b(basis_rows); + lp_assert(m_A * x0 == b); + find_h_minus_1_b(h.W(), b); + int cut_row = find_cut_row_index(b); + + if (cut_row == -1) { + return lia_move::undef; + } + + // the matrix is not square - we can get + // all integers in b's projection + + vector row(m_A.column_count()); + get_ei_H_minus_1(cut_row, h.W(), row); + vector f = row * m_A; + fill_term(f, t); + k = floor(b[cut_row]); + upper = true; + return lia_move::cut; + } + + svector hnf_cutter::vars() const { return m_var_register.vars(); } + + void hnf_cutter::try_add_term_to_A_for_hnf(unsigned i) { + mpq rs; + const lar_term* t = lra.terms()[i]; + constraint_index ci; + bool upper_bound; + if (!is_full() && lra.get_equality_and_right_side_for_term_on_current_x(i, rs, ci, upper_bound)) { + add_term(t, rs, ci, upper_bound); + } + } + + + bool hnf_cutter::hnf_has_var_with_non_integral_value() const { + for (unsigned j : vars()) + if (!lia.get_value(j).is_int()) + return true; + return false; + } + + bool hnf_cutter::init_terms_for_hnf_cut() { + clear(); + for (unsigned i = 0; i < lra.terms().size() && !is_full(); i++) { + try_add_term_to_A_for_hnf(i); + } + return hnf_has_var_with_non_integral_value(); + } + + lia_move hnf_cutter::make_hnf_cut() { + if (!init_terms_for_hnf_cut()) { + return lia_move::undef; + } + lia.settings().stats().m_hnf_cutter_calls++; + TRACE("hnf_cut", tout << "settings().stats().m_hnf_cutter_calls = " << lia.settings().stats().m_hnf_cutter_calls << "\n"; + for (unsigned i : constraints_for_explanation()) { + lra.constraints().display(tout, i); + } + tout << lra.constraints(); + ); +#ifdef Z3DEBUG + vector x0 = transform_to_local_columns(lra.m_mpq_lar_core_solver.m_r_x); +#else + vector x0; +#endif + lia_move r = create_cut(lia.m_t, lia.m_k, lia.m_ex, lia.m_upper, x0); + + if (r == lia_move::cut) { + TRACE("hnf_cut", + lra.print_term(lia.m_t, tout << "cut:"); + tout << " <= " << lia.m_k << std::endl; + for (unsigned i : constraints_for_explanation()) { + lra.constraints().display(tout, i); + } + ); + lp_assert(lia.current_solution_is_inf_on_cut()); + lia.settings().stats().m_hnf_cuts++; + lia.m_ex->clear(); + for (unsigned i : constraints_for_explanation()) { + lia.m_ex->push_justification(i); + } + } + return r; + } + +} diff --git a/src/math/lp/hnf_cutter.h b/src/math/lp/hnf_cutter.h index c155f3336..374995d64 100644 --- a/src/math/lp/hnf_cutter.h +++ b/src/math/lp/hnf_cutter.h @@ -3,19 +3,17 @@ Copyright (c) 2017 Microsoft Corporation Module Name: - + hnf_cutter.h Abstract: - + Cuts (branches) from Hermite matrices Author: Lev Nachmanson (levnach) -Revision History: - - --*/ + #pragma once #include "math/lp/lar_term.h" #include "math/lp/hnf.h" @@ -25,212 +23,64 @@ Revision History: #include "math/lp/explanation.h" namespace lp { +class int_solver; +class lar_solver; + class hnf_cutter { + int_solver& lia; + lar_solver& lra; + lp_settings & m_settings; general_matrix m_A; vector m_terms; vector m_terms_upper; svector m_constraints_for_explanation; vector m_right_sides; - lp_settings & m_settings; mpq m_abs_max; bool m_overflow; var_register m_var_register; + public: + + hnf_cutter(int_solver& lia); + + lia_move hnf_cutter::make_hnf_cut(); + + bool hnf_cutter::init_terms_for_hnf_cut(); + bool hnf_cutter::hnf_has_var_with_non_integral_value() const; + void hnf_cutter::try_add_term_to_A_for_hnf(unsigned i); + + unsigned terms_count() const { return m_terms.size(); } const mpq & abs_max() const { return m_abs_max; } - - hnf_cutter(lp_settings & settings) : m_settings(settings), - m_abs_max(zero_of_type()), - m_var_register(0) {} - - unsigned terms_count() const { - return m_terms.size(); - } - const vector& terms() const { return m_terms; } - const svector& constraints_for_explanation() const { - return m_constraints_for_explanation; - } + const svector& constraints_for_explanation() const { return m_constraints_for_explanation; } const vector & right_sides() const { return m_right_sides; } - void clear() { - // m_A will be filled from scratch in init_matrix_A - m_var_register.clear(); - m_terms.clear(); - m_terms_upper.clear(); - m_constraints_for_explanation.clear(); - m_right_sides.clear(); - m_abs_max = zero_of_type(); - m_overflow = false; - } - void add_term(const lar_term* t, const mpq &rs, constraint_index ci, bool upper_bound) { - m_terms.push_back(t); - m_terms_upper.push_back(upper_bound); - if (upper_bound) - m_right_sides.push_back(rs); - else - m_right_sides.push_back(-rs); - - m_constraints_for_explanation.push_back(ci); - - for (const auto &p : *t) { - m_var_register.add_var(p.var()); - mpq t = abs(ceil(p.coeff())); - if (t > m_abs_max) - m_abs_max = t; - } - } + + bool is_full() const; + + void clear(); + void add_term(const lar_term* t, const mpq &rs, constraint_index ci, bool upper_bound); - void print(std::ostream & out) { - out << "terms = " << m_terms.size() << ", var = " << m_var_register.size() << std::endl; - } + void print(std::ostream & out); - void initialize_row(unsigned i) { - mpq sign = m_terms_upper[i]? one_of_type(): - one_of_type(); - m_A.init_row_from_container(i, * m_terms[i], [this](unsigned j) { return m_var_register.add_var(j);}, sign); - } + void initialize_row(unsigned i); + void init_matrix_A(); + void find_h_minus_1_b(const general_matrix& H, vector & b); - void init_matrix_A() { - m_A = general_matrix(terms_count(), vars().size()); - for (unsigned i = 0; i < terms_count(); i++) - initialize_row(i); - } + vector create_b(const svector & basis_rows); - // todo: as we need only one row i with non integral b[i] need to optimize later - void find_h_minus_1_b(const general_matrix& H, vector & b) { - // the solution will be put into b - for (unsigned i = 0; i < H.row_count() ;i++) { - for (unsigned j = 0; j < i; j++) { - b[i] -= H[i][j]*b[j]; - } - b[i] /= H[i][i]; - // consider return from here if b[i] is not an integer and return i - } - } - - vector create_b(const svector & basis_rows) { - if (basis_rows.size() == m_right_sides.size()) - return m_right_sides; - vector b; - for (unsigned i : basis_rows) { - b.push_back(m_right_sides[i]); - } - return b; - } - - int find_cut_row_index(const vector & b) { - int ret = -1; - int n = 0; - for (int i = 0; i < static_cast(b.size()); i++) { - if (is_integer(b[i])) continue; - if (n == 0 ) { - lp_assert(ret == -1); - n = 1; - ret = i; - } else { - if (m_settings.random_next() % (++n) == 0) { - ret = i; - } - } - } - return ret; - } + int find_cut_row_index(const vector & b); // fills e_i*H_minus_1 - void get_ei_H_minus_1(unsigned i, const general_matrix& H, vector & row) { - // we solve x = ei * H_min_1 - // or x * H = ei - unsigned m = H.row_count(); - for (unsigned k = i + 1; k < m; k++) { - row[k] = zero_of_type(); - } - row[i] = one_of_type() / H[i][i]; - for(int k = i - 1; k >= 0; k--) { - mpq t = zero_of_type(); - for (unsigned l = k + 1; l <= i; l++) { - t += H[l][k]*row[l]; - } - row[k] = -t / H[k][k]; - } + void get_ei_H_minus_1(unsigned i, const general_matrix& H, vector & row); - // // test region - // vector ei(H.row_count(), zero_of_type()); - // ei[i] = one_of_type(); - // vector pr = row * H; - // pr.shrink(ei.size()); - // lp_assert(ei == pr); - // // end test region - - } - - void fill_term(const vector & row, lar_term& t) { - for (unsigned j = 0; j < row.size(); j++) { - if (!is_zero(row[j])) - t.add_monomial(row[j], m_var_register.local_to_external(j)); - } - } + void fill_term(const vector & row, lar_term& t); #ifdef Z3DEBUG - vector transform_to_local_columns(const vector & x) const { - vector ret; - for (unsigned j = 0; j < vars().size(); j++) { - ret.push_back(x[m_var_register.local_to_external(j)].x); - } - return ret; - } + vector transform_to_local_columns(const vector & x) const; #endif - void shrink_explanation(const svector& basis_rows) { - svector new_expl; - for (unsigned i : basis_rows) { - new_expl.push_back(m_constraints_for_explanation[i]); - } - m_constraints_for_explanation = new_expl; - } - - bool overflow() const { return m_overflow; } - - lia_move create_cut(lar_term& t, mpq& k, explanation* ex, bool & upper, const vector & x0) { - // we suppose that x0 has at least one non integer element - (void)x0; - - init_matrix_A(); - svector basis_rows; - mpq big_number = m_abs_max.expt(3); - mpq d = hnf_calc::determinant_of_rectangular_matrix(m_A, basis_rows, big_number); - - if (d >= big_number) { - return lia_move::undef; - } - - if (m_settings.get_cancel_flag()) { - return lia_move::undef; - } - - if (basis_rows.size() < m_A.row_count()) { - m_A.shrink_to_rank(basis_rows); - shrink_explanation(basis_rows); - } - - hnf h(m_A, d); - vector b = create_b(basis_rows); - lp_assert(m_A * x0 == b); - find_h_minus_1_b(h.W(), b); - int cut_row = find_cut_row_index(b); - - if (cut_row == -1) { - return lia_move::undef; - } - - // the matrix is not square - we can get - // all integers in b's projection - - vector row(m_A.column_count()); - get_ei_H_minus_1(cut_row, h.W(), row); - vector f = row * m_A; - fill_term(f, t); - k = floor(b[cut_row]); - upper = true; - return lia_move::cut; - } - - svector vars() const { return m_var_register.vars(); } + void shrink_explanation(const svector& basis_rows); + bool overflow() const; + lia_move create_cut(lar_term& t, mpq& k, explanation* ex, bool & upper, const vector & x0); + svector vars() const; }; } diff --git a/src/math/lp/int_solver.cpp b/src/math/lp/int_solver.cpp index d41c8d024..3a718a3f4 100644 --- a/src/math/lp/int_solver.cpp +++ b/src/math/lp/int_solver.cpp @@ -18,7 +18,7 @@ namespace lp { int_solver::int_solver(lar_solver& lar_slv) : lra(lar_slv), m_number_of_calls(0), - m_hnf_cutter(settings()), + m_hnf_cutter(*this), m_hnf_cut_period(settings().hnf_cut_period()) { lra.set_int_solver(this); } @@ -172,80 +172,12 @@ bool int_solver::should_gomory_cut() { return m_number_of_calls % settings().m_int_gomory_cut_period == 0; } -void int_solver::try_add_term_to_A_for_hnf(unsigned i) { - mpq rs; - const lar_term* t = lra.terms()[i]; - constraint_index ci; - bool upper_bound; - if (!hnf_cutter_is_full() && lra.get_equality_and_right_side_for_term_on_current_x(i, rs, ci, upper_bound)) { - m_hnf_cutter.add_term(t, rs, ci, upper_bound); - } -} - -bool int_solver::hnf_cutter_is_full() const { - return - m_hnf_cutter.terms_count() >= settings().limit_on_rows_for_hnf_cutter - || - m_hnf_cutter.vars().size() >= settings().limit_on_columns_for_hnf_cutter; -} - -bool int_solver::hnf_has_var_with_non_integral_value() const { - for (unsigned j : m_hnf_cutter.vars()) - if (!get_value(j).is_int()) - return true; - return false; -} - -bool int_solver::init_terms_for_hnf_cut() { - m_hnf_cutter.clear(); - for (unsigned i = 0; i < lra.terms().size() && !hnf_cutter_is_full(); i++) { - try_add_term_to_A_for_hnf(i); - } - return hnf_has_var_with_non_integral_value(); -} - -lia_move int_solver::make_hnf_cut() { - if (!init_terms_for_hnf_cut()) { - return lia_move::undef; - } - settings().stats().m_hnf_cutter_calls++; - TRACE("hnf_cut", tout << "settings().stats().m_hnf_cutter_calls = " << settings().stats().m_hnf_cutter_calls << "\n"; - for (unsigned i : m_hnf_cutter.constraints_for_explanation()) { - lra.constraints().display(tout, i); - } - tout << lra.constraints(); - ); -#ifdef Z3DEBUG - vector x0 = m_hnf_cutter.transform_to_local_columns(lra.m_mpq_lar_core_solver.m_r_x); -#else - vector x0; -#endif - lia_move r = m_hnf_cutter.create_cut(m_t, m_k, m_ex, m_upper, x0); - - if (r == lia_move::cut) { - TRACE("hnf_cut", - lra.print_term(m_t, tout << "cut:"); - tout << " <= " << m_k << std::endl; - for (unsigned i : m_hnf_cutter.constraints_for_explanation()) { - lra.constraints().display(tout, i); - } - ); - lp_assert(current_solution_is_inf_on_cut()); - settings().stats().m_hnf_cuts++; - m_ex->clear(); - for (unsigned i : m_hnf_cutter.constraints_for_explanation()) { - m_ex->push_justification(i); - } - } - return r; -} - bool int_solver::should_hnf_cut() { return settings().m_enable_hnf && m_number_of_calls % m_hnf_cut_period == 0; } lia_move int_solver::hnf_cut() { - lia_move r = make_hnf_cut(); + lia_move r = m_hnf_cutter.make_hnf_cut(); if (r == lia_move::undef) { m_hnf_cut_period *= 2; } diff --git a/src/math/lp/int_solver.h b/src/math/lp/int_solver.h index 7ebc9ce5d..154f59ab4 100644 --- a/src/math/lp/int_solver.h +++ b/src/math/lp/int_solver.h @@ -36,6 +36,7 @@ class int_solver { friend class int_cube; friend class int_branch; friend class int_gcd_test; + friend class hnf_cutter; lar_solver& lra; unsigned m_number_of_calls; @@ -55,7 +56,6 @@ public: lar_term const& get_term() const { return m_t; } mpq const& get_offset() const { return m_k; } bool is_upper() const { return m_upper; } - //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; @@ -106,12 +106,6 @@ public: bool all_columns_are_bounded() const; void find_feasible_solution(); lia_move hnf_cut(); - lia_move make_hnf_cut(); - bool init_terms_for_hnf_cut(); - bool hnf_matrix_is_empty() const; - void try_add_term_to_A_for_hnf(unsigned term_index); - bool hnf_has_var_with_non_integral_value() const; - bool hnf_cutter_is_full() const; void patch_nbasic_column(unsigned j); }; }