3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2025-04-23 17:15:31 +00:00

cleanup in hnf.h

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

fixex in maximize_term

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

prepare for LIRA (#76)

* merge

Signed-off-by: Nikolaj Bjorner <nbjorner@microsoft.com>

* simple mixed integer example

Signed-off-by: Nikolaj Bjorner <nbjorner@microsoft.com>

Nikolaj's fixes in theory_lra.cpp

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

fixes in maximize_term, disable find_cube for mixed case

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

fix cube heuristic for the mixed case

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

fix the build

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

fix in find cube delta's calculation

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

remove a printout

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

remove a blank line

Signed-off-by: Lev Nachmanson <levnach@microsoft.com>

test credentials

Signed-off-by: Lev Nachmanson <levnach@microsoft.com>

avoid double checks to add terms in hnf_cutter

Signed-off-by: Lev Nachmanson <levnach@microsoft.com>

fixes in add terms to hnf_cutter

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

remove a variable used for debug only

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

remove a field

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

remove a warning and hide m_A_orig under debug

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

use var_register to deal with mapping between external and local variables

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

tighten the loop in explain_implied_bound

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

fixes in theory_lra and relaxing debug checks in pop(), for the case when push() has been done from not fully initialized state

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

suppress hnf cutter when the numbers become too large

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

remove some debug code

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

adjusting hnf

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>
This commit is contained in:
Lev Nachmanson 2018-06-08 15:05:08 -07:00
parent 9ba4026bc6
commit eeaca949e0
13 changed files with 424 additions and 260 deletions

View file

@ -126,33 +126,35 @@ bool prepare_pivot_for_lower_triangle(M &m, unsigned r) {
}
template <typename M>
void pivot_column_non_fractional(M &m, unsigned r) {
void pivot_column_non_fractional(M &m, unsigned r, bool & overflow, const mpq & big_number) {
lp_assert(!is_zero(m[r][r]));
for (unsigned j = r + 1; j < m.column_count(); j++) {
for (unsigned i = r + 1; i < m.row_count(); i++) {
for (unsigned i = r + 1; i < m.row_count(); i++) {
m[i][j] =
(r > 0) ?
(m[r][r]*m[i][j] - m[i][r]*m[r][j]) / m[r-1][r-1] :
(m[r][r]*m[i][j] - m[i][r]*m[r][j]);
if (m[i][j] >= big_number) {
overflow = true;
return;
}
lp_assert(is_int(m[i][j]));
}
}
// debug
for (unsigned k = r + 1; k < m.row_count(); k++) {
m[k][r] = zero_of_type<mpq>();
}
}
// returns the rank of the matrix
template <typename M>
unsigned to_lower_triangle_non_fractional(M &m) {
unsigned to_lower_triangle_non_fractional(M &m, bool & overflow, const mpq& big_number) {
unsigned i = 0;
for (; i < m.row_count(); i++) {
if (!prepare_pivot_for_lower_triangle(m, i)) {
return i;
}
pivot_column_non_fractional(m, i);
pivot_column_non_fractional(m, i, overflow, big_number);
if (overflow)
return 0;
}
lp_assert(i == m.row_count());
return i;
@ -184,9 +186,12 @@ mpq gcd_of_row_starting_from_diagonal(const M& m, unsigned i) {
// The gcd of these minors is the return value.
template <typename M>
mpq determinant_of_rectangular_matrix(const M& m, svector<unsigned> & basis_rows) {
mpq determinant_of_rectangular_matrix(const M& m, svector<unsigned> & basis_rows, const mpq& big_number) {
auto m_copy = m;
unsigned rank = to_lower_triangle_non_fractional(m_copy);
bool overflow = false;
unsigned rank = to_lower_triangle_non_fractional(m_copy, overflow, big_number);
if (overflow)
return big_number;
if (rank == 0)
return one_of_type<mpq>();
@ -216,9 +221,8 @@ class hnf {
M m_H;
M m_U;
M m_U_reverse;
#endif
M m_A_orig;
#endif
M m_W;
vector<mpq> m_buffer;
unsigned m_m;
@ -529,14 +533,14 @@ private:
// with some changes related to that we create a low triangle matrix
// with non-positive elements under the diagonal
void process_column_in_row_modulo() {
mpq& aii = m_W[m_i][m_i];
const mpq& aii = m_W[m_i][m_i];
const mpq& aij = m_W[m_i][m_j];
mpq d, p,q;
hnf_calc::extended_gcd_minimal_uv(aii, aij, d, p, q);
if (is_zero(d))
return;
mpq aii_over_d = aii / d;
mpq aij_over_d = aij / d;
mpq aii_over_d = mod_R(aii / d);
mpq aij_over_d = mod_R(aij / d);
buffer_p_col_i_plus_q_col_j_W_modulo(p, q);
pivot_column_i_to_column_j_W_modulo(- aij_over_d, aii_over_d);
copy_buffer_to_col_i_W_modulo();

View file

@ -31,15 +31,18 @@ class hnf_cutter {
vector<const lar_term*> m_terms;
svector<constraint_index> m_constraints_for_explanation;
vector<mpq> m_right_sides;
unsigned m_row_count;
unsigned m_column_count;
lp_settings & m_settings;
mpq m_abs_max;
bool m_overflow;
public:
hnf_cutter(lp_settings & settings) : m_row_count(0), m_column_count(0), m_settings(settings) {}
unsigned row_count() const {
return m_row_count;
const mpq & abs_max() const { return m_abs_max; }
hnf_cutter(lp_settings & settings) : m_settings(settings),
m_abs_max(zero_of_type<mpq>()) {}
unsigned terms_count() const {
return m_terms.size();
}
const vector<const lar_term*>& terms() const { return m_terms; }
@ -52,8 +55,9 @@ public:
m_var_register.clear();
m_terms.clear();
m_constraints_for_explanation.clear();
m_right_sides.clear();
m_row_count = m_column_count = 0;
m_right_sides.clear();
m_abs_max = zero_of_type<mpq>();
m_overflow = false;
}
void add_term(const lar_term* t, const mpq &rs, constraint_index ci) {
m_terms.push_back(t);
@ -61,12 +65,12 @@ public:
m_constraints_for_explanation.push_back(ci);
for (const auto &p : *t) {
m_var_register.add_var(p.var());
}
if (m_terms.size() <= m_var_register.size()) { // capture the maximal acceptable sizes
m_row_count = m_terms.size();
m_column_count = m_var_register.size();
mpq t = abs(ceil(p.coeff()));
if (t > m_abs_max)
m_abs_max = t;
}
}
void print(std::ostream & out) {
out << "terms = " << m_terms.size() << ", var = " << m_var_register.size() << std::endl;
}
@ -76,10 +80,8 @@ public:
}
void init_matrix_A() {
m_A = general_matrix(m_row_count, m_column_count);
// use the last suitable counts to make the number
// of rows less than or equal to the number of columns
for (unsigned i = 0; i < m_row_count; i++)
m_A = general_matrix(terms_count(), vars().size());
for (unsigned i = 0; i < terms_count(); i++)
initialize_row(i);
}
@ -123,7 +125,6 @@ public:
return ret;
}
// fills e_i*H_minus_1
void get_ei_H_minus_1(unsigned i, const general_matrix& H, vector<mpq> & row) {
// we solve x = ei * H_min_1
@ -141,29 +142,28 @@ public:
row[k] = -t / H[k][k];
}
// test region
vector<mpq> ei(H.row_count(), zero_of_type<mpq>());
ei[i] = one_of_type<mpq>();
vector<mpq> pr = row * H;
pr.shrink(ei.size());
lp_assert(ei == pr);
// end test region
// // test region
// vector<mpq> ei(H.row_count(), zero_of_type<mpq>());
// ei[i] = one_of_type<mpq>();
// vector<mpq> pr = row * H;
// pr.shrink(ei.size());
// lp_assert(ei == pr);
// // end test region
}
void fill_term(const vector<mpq> & 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_var_to_user_var(j));
t.add_monomial(row[j], m_var_register.local_to_external(j));
}
}
#ifdef Z3DEBUG
vector<mpq> transform_to_local_columns(const vector<impq> & x) const {
vector<mpq> ret;
lp_assert(m_column_count <= m_var_register.size());
for (unsigned j = 0; j < m_column_count;j++) {
lp_assert(is_zero(x[m_var_register.local_var_to_user_var(j)].y));
ret.push_back(x[m_var_register.local_var_to_user_var(j)].x);
for (unsigned j = 0; j < vars().size(); j++) {
lp_assert(is_zero(x[m_var_register.local_to_external(j)].y));
ret.push_back(x[m_var_register.local_to_external(j)].x);
}
return ret;
}
@ -176,6 +176,8 @@ public:
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
#ifdef Z3DEBUG
,
@ -185,7 +187,16 @@ public:
// we suppose that x0 has at least one non integer element
init_matrix_A();
svector<unsigned> basis_rows;
mpq d = hnf_calc::determinant_of_rectangular_matrix(m_A, basis_rows);
mpq big_number = m_abs_max.expt(3);
mpq d = hnf_calc::determinant_of_rectangular_matrix(m_A, basis_rows, big_number);
// std::cout << "max = " << m_abs_max << ", d = " << d << ", d/max = " << ceil (d /m_abs_max) << std::endl;
//std::cout << "max cube " << m_abs_max * m_abs_max * m_abs_max << std::endl;
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()) {

View file

@ -416,6 +416,8 @@ impq int_solver::get_cube_delta_for_term(const lar_term& t) const {
bool seen_minus = false;
bool seen_plus = false;
for(const auto & p : t) {
if (!is_int(p.var()))
goto usual_delta;
const mpq & c = p.coeff();
if (c == one_of_type<mpq>()) {
seen_plus = true;
@ -431,9 +433,10 @@ impq int_solver::get_cube_delta_for_term(const lar_term& t) const {
}
usual_delta:
mpq delta = zero_of_type<mpq>();
for (const auto & p : t) {
delta += abs(p.coeff());
}
for (const auto & p : t)
if (is_int(p.var()))
delta += abs(p.coeff());
delta *= mpq(1, 2);
return impq(delta);
}
@ -443,7 +446,6 @@ bool int_solver::tighten_term_for_cube(unsigned i) {
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))
@ -533,20 +535,28 @@ lia_move int_solver::gomory_cut() {
void int_solver::try_add_term_to_A_for_hnf(unsigned i) {
mpq rs;
const lar_term* t = m_lar_solver->terms()[i];
#if Z3DEBUG
for (const auto & p : *t) {
if (!is_int(p.var())) {
lp_assert(false);
}
}
#endif
constraint_index ci;
if (m_lar_solver->get_equality_and_right_side_for_term_on_current_x(i, rs, ci)) {
if (!hnf_cutter_is_full() && m_lar_solver->get_equality_and_right_side_for_term_on_current_x(i, rs, ci)) {
m_hnf_cutter.add_term(t, rs, ci);
}
}
bool int_solver::hnf_has_non_integral_var() const {
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;
}
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() == false)
return true;
@ -555,10 +565,10 @@ bool int_solver::hnf_has_non_integral_var() const {
bool int_solver::init_terms_for_hnf_cut() {
m_hnf_cutter.clear();
for (unsigned i = 0; i < m_lar_solver->terms().size() && m_hnf_cutter.row_count() < settings().limit_on_rows_for_hnf_cutter; i++) {
for (unsigned i = 0; i < m_lar_solver->terms().size() && !hnf_cutter_is_full(); i++) {
try_add_term_to_A_for_hnf(i);
}
return hnf_has_non_integral_var();
return hnf_has_var_with_non_integral_value();
}
lia_move int_solver::make_hnf_cut() {
@ -583,7 +593,7 @@ lia_move int_solver::make_hnf_cut() {
for (unsigned i : m_hnf_cutter.constraints_for_explanation()) {
m_ex->push_justification(i);
}
}
}
return r;
}
@ -609,7 +619,6 @@ lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex, bool & upper) {
r = patch_nbasic_columns();
if (r != lia_move::undef) return r;
++m_branch_cut_counter;
r = find_cube();
if (r != lia_move::undef) return r;
@ -697,14 +706,12 @@ void int_solver::set_value_for_nbasic_column(unsigned j, const impq & new_val) {
m_lar_solver->change_basic_columns_dependend_on_a_given_nb_column(j, delta);
}
void int_solver::patch_nbasic_column(unsigned j) {
void int_solver::patch_nbasic_column(unsigned j, bool patch_only_int_vals) {
auto & lcs = m_lar_solver->m_mpq_lar_core_solver;
impq & val = lcs.m_r_x[j];
bool val_is_int = val.is_int();
if (settings().m_int_patch_only_integer_values) {
if (!val_is_int)
if (patch_only_int_vals && !val_is_int)
return;
}
bool inf_l, inf_u;
impq l, u;
@ -754,7 +761,7 @@ lia_move int_solver::patch_nbasic_columns() {
settings().st().m_patches++;
lp_assert(is_feasible());
for (unsigned j : m_lar_solver->m_mpq_lar_core_solver.m_r_nbasis) {
patch_nbasic_column(j);
patch_nbasic_column(j, settings().m_int_patch_only_integer_values);
}
lp_assert(is_feasible());
if (!has_inf_int()) {
@ -1128,12 +1135,6 @@ bool int_solver::at_upper(unsigned j) const {
}
}
lp_settings& int_solver::settings() {
return m_lar_solver->settings();
}
void int_solver::display_row_info(std::ostream & out, unsigned row_index) const {
auto & rslv = m_lar_solver->m_mpq_lar_core_solver.m_r_solver;
for (auto &c: rslv.m_A.m_rows[row_index]) {

View file

@ -50,6 +50,7 @@ public:
// main function to check that the solution provided by lar_solver is valid for integral values,
// or provide a way of how it can be adjusted.
lia_move check(lar_term& t, mpq& k, explanation& ex, bool & upper);
lia_move check_(lar_term& t, mpq& k, explanation& ex, bool & upper);
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;
@ -100,6 +101,7 @@ private:
int find_inf_int_boxed_base_column_with_smallest_range(unsigned&);
int get_kth_inf_int(unsigned) const;
lp_settings& settings();
const lp_settings& settings() const;
bool move_non_basic_columns_to_bounds();
void branch_infeasible_int_var(unsigned);
lia_move mk_gomory_cut(unsigned inf_col, const row_strip<mpq>& row);
@ -157,7 +159,8 @@ public:
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_non_integral_var() const;
void patch_nbasic_column(unsigned j);
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);
};
}

View file

@ -31,8 +31,7 @@ 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_int_solver(nullptr),
m_has_int_var(false)
m_int_solver(nullptr)
{}
void lar_solver::set_track_pivoted_rows(bool v) {
@ -199,7 +198,7 @@ void lar_solver::calculate_implied_bounds_for_row(unsigned i, bound_propagator &
}
unsigned lar_solver::adjust_column_index_to_term_index(unsigned j) const {
unsigned ext_var_or_term = m_columns_to_ext_vars_or_term_indices[j];
unsigned ext_var_or_term = m_var_register.local_to_external(j);
return ext_var_or_term < m_terms_start_index ? j : ext_var_or_term;
}
@ -212,21 +211,14 @@ void lar_solver::explain_implied_bound(implied_bound & ib, bound_propagator & bp
unsigned i = ib.m_row_or_term_index;
int bound_sign = ib.m_is_lower_bound? 1: -1;
int j_sign = (ib.m_coeff_before_j_is_pos ? 1 :-1) * bound_sign;
unsigned m_j = ib.m_j;
if (is_term(m_j)) {
auto it = m_ext_vars_to_columns.find(m_j);
lp_assert(it != m_ext_vars_to_columns.end());
m_j = it->second.internal_j();
unsigned bound_j = ib.m_j;
if (is_term(bound_j)) {
bound_j = m_var_register.external_to_local(bound_j);
}
for (auto const& r : A_r().m_rows[i]) {
unsigned j = r.m_j;
mpq const& a = r.get_val();
if (j == m_j) continue;
if (is_term(j)) {
auto it = m_ext_vars_to_columns.find(j);
lp_assert(it != m_ext_vars_to_columns.end());
j = it->second.internal_j();
}
mpq const& a = r.get_val();
int a_sign = is_pos(a)? 1: -1;
int sign = j_sign * a_sign;
const ul_pair & ul = m_columns_to_ul_pairs[j];
@ -240,7 +232,7 @@ void lar_solver::explain_implied_bound(implied_bound & ib, bound_propagator & bp
bool lar_solver::term_is_used_as_row(unsigned term) const {
lp_assert(is_term(term));
return contains(m_ext_vars_to_columns, term);
return m_var_register.external_is_used(term);
}
void lar_solver::propagate_bounds_on_terms(bound_propagator & bp) {
@ -273,10 +265,6 @@ lp_status lar_solver::get_status() const { return m_status;}
void lar_solver::set_status(lp_status s) {m_status = s;}
bool lar_solver::has_int_var() const {
return m_has_int_var;
}
lp_status lar_solver::find_feasible_solution() {
m_settings.st().m_make_feasible++;
if (A_r().column_count() > m_settings.st().m_max_cols)
@ -352,12 +340,9 @@ void lar_solver::pop(unsigned k) {
TRACE("arith_int", tout << "pop" << std::endl;);
TRACE("lar_solver", tout << "k = " << k << std::endl;);
int n_was = static_cast<int>(m_ext_vars_to_columns.size());
m_infeasible_column_index.pop(k);
unsigned n = m_columns_to_ul_pairs.peek_size(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);
m_var_register.shrink(n);
TRACE("arith_int", tout << "pop" << std::endl;);
if (m_settings.use_tableau()) {
pop_tableau();
@ -375,7 +360,6 @@ void lar_solver::pop(unsigned k) {
lp_assert(ax_is_correct());
lp_assert(m_mpq_lar_core_solver.m_r_solver.inf_set_is_correct());
m_constraint_count.pop(k);
for (unsigned i = m_constraint_count; i < m_constraints.size(); i++)
delete m_constraints[i];
@ -404,18 +388,17 @@ vector<constraint_index> lar_solver::get_all_constraint_indices() const {
return ret;
}
bool lar_solver::maximize_term_on_tableau(const vector<std::pair<mpq, var_index>> & term,
bool lar_solver::maximize_term_on_tableau(const lar_term & term,
impq &term_max) {
if (settings().simplex_strategy() == simplex_strategy_enum::undecided)
decide_on_strategy_and_adjust_initial_state();
m_mpq_lar_core_solver.m_r_solver.set_status(lp_status::FEASIBLE);
m_mpq_lar_core_solver.solve();
if (m_mpq_lar_core_solver.m_r_solver.get_status() == lp_status::UNBOUNDED)
return false;
term_max = 0;
for (auto & p : term)
term_max += p.first * m_mpq_lar_core_solver.m_r_x[p.second];
term_max = term.apply(m_mpq_lar_core_solver.m_r_x);
return true;
}
@ -433,13 +416,13 @@ bool lar_solver::reduced_costs_are_zeroes_for_r_solver() const {
return true;
}
void lar_solver::set_costs_to_zero(const vector<std::pair<mpq, var_index>> & term) {
void lar_solver::set_costs_to_zero(const lar_term& term) {
auto & rslv = m_mpq_lar_core_solver.m_r_solver;
auto & jset = m_mpq_lar_core_solver.m_r_solver.m_inf_set; // hijack this set that should be empty right now
lp_assert(jset.m_index.size()==0);
for (auto & p : term) {
unsigned j = p.second;
for (const auto & p : term) {
unsigned j = p.var();
rslv.m_costs[j] = zero_of_type<mpq>();
int i = rslv.m_basis_heading[j];
if (i < 0)
@ -459,25 +442,25 @@ void lar_solver::set_costs_to_zero(const vector<std::pair<mpq, var_index>> & ter
lp_assert(costs_are_zeros_for_r_solver());
}
void lar_solver::prepare_costs_for_r_solver(const vector<std::pair<mpq, var_index>> & term) {
void lar_solver::prepare_costs_for_r_solver(const lar_term & term) {
auto & rslv = m_mpq_lar_core_solver.m_r_solver;
rslv.m_using_infeas_costs = false;
lp_assert(costs_are_zeros_for_r_solver());
lp_assert(reduced_costs_are_zeroes_for_r_solver());
rslv.m_costs.resize(A_r().column_count(), zero_of_type<mpq>());
for (auto & p : term) {
unsigned j = p.second;
rslv.m_costs[j] = p.first;
for (const auto & p : term) {
unsigned j = p.var();
rslv.m_costs[j] = p.coeff();
if (rslv.m_basis_heading[j] < 0)
rslv.m_d[j] += p.first;
rslv.m_d[j] += p.coeff();
else
rslv.update_reduced_cost_for_basic_column_cost_change(- p.first, j);
rslv.update_reduced_cost_for_basic_column_cost_change(- p.coeff(), j);
}
lp_assert(rslv.reduced_costs_are_correct_tableau());
}
bool lar_solver::maximize_term_on_corrected_r_solver(const vector<std::pair<mpq, var_index>> & term,
bool lar_solver::maximize_term_on_corrected_r_solver(lar_term & term,
impq &term_max) {
settings().backup_costs = false;
switch (settings().simplex_strategy()) {
@ -514,16 +497,35 @@ bool lar_solver::remove_from_basis(unsigned j) {
return m_mpq_lar_core_solver.m_r_solver.remove_from_basis(j);
}
lar_term lar_solver::get_term_to_maximize(unsigned ext_j) const {
unsigned local_j;
if (m_var_register.external_is_used(ext_j, local_j)) {
lar_term r;
r. add_monomial(one_of_type<mpq>(), local_j);
return r;
}
return get_term(ext_j);
}
// starting from a given feasible state look for the maximum of the term
// return true if found and false if unbounded
lp_status lar_solver::maximize_term(const vector<std::pair<mpq, var_index>> & term,
lp_status lar_solver::maximize_term(unsigned ext_j,
impq &term_max) {
lp_assert(m_mpq_lar_core_solver.m_r_solver.current_x_is_feasible());
bool was_feasible = m_mpq_lar_core_solver.m_r_solver.calc_current_x_is_feasible_include_non_basis();
impq prev_value;
lar_term term = get_term_to_maximize(ext_j);
auto backup = m_mpq_lar_core_solver.m_r_x;
if (was_feasible) {
prev_value = term.apply(m_mpq_lar_core_solver.m_r_x);
}
m_mpq_lar_core_solver.m_r_solver.m_look_for_feasible_solution_only = false;
if (!maximize_term_on_corrected_r_solver(term, term_max))
if (!maximize_term_on_corrected_r_solver(term, term_max)) {
m_mpq_lar_core_solver.m_r_x = backup;
return lp_status::UNBOUNDED;
}
impq opt_val = term_max;
bool change = false;
for (unsigned j = 0; j < m_mpq_lar_core_solver.m_r_x.size(); j++) {
if (!column_is_int(j))
@ -534,18 +536,19 @@ lp_status lar_solver::maximize_term(const vector<std::pair<mpq, var_index>> & te
if (!remove_from_basis(j)) // consider a special version of remove_from_basis that would not remove inf_int columns
return lp_status::FEASIBLE; // it should not happen
}
m_int_solver->patch_nbasic_column(j);
m_int_solver->patch_nbasic_column(j, false);
if (!column_value_is_integer(j))
return lp_status::FEASIBLE;
change = true;
}
if (change) {
term_max = zero_of_type<impq>();
for (const auto& p : term) {
term_max += p.first * m_mpq_lar_core_solver.m_r_x[p.second];
}
term_max = term.apply(m_mpq_lar_core_solver.m_r_x);
}
return lp_status::OPTIMAL;
if (was_feasible && term_max < prev_value) {
term_max = prev_value;
m_mpq_lar_core_solver.m_r_x = backup;
}
return term_max == opt_val? lp_status::OPTIMAL :lp_status::FEASIBLE;
}
@ -907,10 +910,10 @@ column_type lar_solver::get_column_type(const column_info<mpq> & ci) {
std::string lar_solver::get_column_name(unsigned j) const {
if (j >= m_terms_start_index)
return std::string("_t") + T_to_string(j);
if (j >= m_columns_to_ext_vars_or_term_indices.size())
if (j >= m_var_register.size())
return std::string("_s") + T_to_string(j);
return std::string("v") + T_to_string(m_columns_to_ext_vars_or_term_indices[j]);
return std::string("v") + T_to_string(m_var_register.local_to_external(j));
}
bool lar_solver::all_constrained_variables_are_registered(const vector<std::pair<mpq, var_index>>& left_side) {
@ -1392,7 +1395,8 @@ void lar_solver::pop_tableau() {
lp_assert(m_mpq_lar_core_solver.m_r_solver.basis_heading_is_correct());
// We remove last variables starting from m_column_names.size() to m_vec_of_canonic_left_sides.size().
// At this moment m_column_names is already popped
while (A_r().column_count() > m_columns_to_ext_vars_or_term_indices.size())
unsigned size = m_var_register.size();
while (A_r().column_count() > size)
remove_last_column_from_tableau();
lp_assert(m_mpq_lar_core_solver.m_r_solver.m_costs.size() == A_r().column_count());
lp_assert(m_mpq_lar_core_solver.m_r_solver.m_basis.size() == A_r().row_count());
@ -1469,9 +1473,7 @@ bool lar_solver::var_is_int(var_index v) const {
}
bool lar_solver::column_is_int(unsigned j) const {
unsigned ext_var = m_columns_to_ext_vars_or_term_indices[j];
lp_assert(contains(m_ext_vars_to_columns, ext_var));
return m_ext_vars_to_columns.find(ext_var)->second.is_integer();
return m_var_register.local_is_int(j);
}
bool lar_solver::column_is_fixed(unsigned j) const {
@ -1480,9 +1482,7 @@ bool lar_solver::column_is_fixed(unsigned j) const {
bool lar_solver::ext_var_is_int(var_index ext_var) const {
auto it = m_ext_vars_to_columns.find(ext_var);
lp_assert(it != m_ext_vars_to_columns.end());
return it->second.is_integer();
return m_var_register.external_is_int(ext_var);
}
// below is the initialization functionality of lar_solver
@ -1492,30 +1492,22 @@ bool lar_solver::strategy_is_undecided() const {
}
var_index lar_solver::add_var(unsigned ext_j, bool is_int) {
if (is_int)
m_has_int_var = true;
TRACE("add_var", tout << "adding var " << ext_j << (is_int? " int" : " nonint") << std::endl;);
var_index i;
var_index local_j;
lp_assert(ext_j < m_terms_start_index);
auto it = m_ext_vars_to_columns.find(ext_j);
if (it != m_ext_vars_to_columns.end()) {
return it->second.internal_j();
}
if (m_var_register.external_is_used(ext_j, local_j))
return local_j;
lp_assert(m_columns_to_ul_pairs.size() == A_r().column_count());
i = A_r().column_count();
local_j = A_r().column_count();
m_columns_to_ul_pairs.push_back(ul_pair(static_cast<unsigned>(-1)));
add_non_basic_var_to_core_fields(ext_j, is_int);
lp_assert(sizes_are_correct());
return i;
return local_j;
}
void lar_solver::register_new_ext_var_index(unsigned ext_v, bool is_int) {
lp_assert(!contains(m_ext_vars_to_columns, ext_v));
unsigned j = static_cast<unsigned>(m_ext_vars_to_columns.size());
m_ext_vars_to_columns.insert(std::make_pair(ext_v, ext_var_info(j, is_int)));
lp_assert(m_columns_to_ext_vars_or_term_indices.size() == j);
m_columns_to_ext_vars_or_term_indices.push_back(ext_v);
lp_assert(!m_var_register.external_is_used(ext_v));
m_var_register.add_var(ext_v, is_int);
}
void lar_solver::add_non_basic_var_to_core_fields(unsigned ext_j, bool is_int) {
@ -1631,7 +1623,7 @@ var_index lar_solver::add_term(const vector<std::pair<mpq, var_index>> & coeffs,
if (m_settings.bound_propagation())
m_rows_with_changed_bounds.insert(A_r().row_count() - 1);
}
lp_assert(m_ext_vars_to_columns.size() == A_r().column_count());
lp_assert(m_var_register.size() == A_r().column_count());
return ret;
}
@ -1716,9 +1708,8 @@ void lar_solver::add_var_bound_on_constraint_for_term(var_index j, lconstraint_k
lp_assert(is_term(j));
unsigned adjusted_term_index = adjust_term_index(j);
// lp_assert(!term_is_int(m_terms[adjusted_term_index]) || right_side.is_int());
auto it = m_ext_vars_to_columns.find(j);
if (it != m_ext_vars_to_columns.end()) {
unsigned term_j = it->second.internal_j();
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);
@ -1804,10 +1795,10 @@ void lar_solver::adjust_initial_state_for_lu() {
}
void lar_solver::adjust_initial_state_for_tableau_rows() {
for (unsigned j = 0; j < m_terms.size(); j++) {
if (contains(m_ext_vars_to_columns, j + m_terms_start_index))
for (unsigned i = 0; i < m_terms.size(); i++) {
if (m_var_register.external_is_used(i + m_terms_start_index))
continue;
add_row_from_term_no_constraint(m_terms[j], j + m_terms_start_index);
add_row_from_term_no_constraint(m_terms[i], i + m_terms_start_index);
}
}
@ -2122,21 +2113,18 @@ void lar_solver::update_fixed_column_type_and_bound(var_index j, lconstraint_kin
}
bool lar_solver::column_corresponds_to_term(unsigned j) const {
return m_columns_to_ext_vars_or_term_indices[j] >= m_terms_start_index;
return m_var_register.local_to_external(j) >= m_terms_start_index;
}
var_index lar_solver:: to_var_index(unsigned ext_j) const {
auto it = m_ext_vars_to_columns.find(ext_j);
lp_assert(it != m_ext_vars_to_columns.end());
return it->second.internal_j();
var_index lar_solver::to_column(unsigned ext_j) const {
return m_var_register.external_to_local(ext_j);
}
bool lar_solver::tighten_term_bounds_by_delta(unsigned term_index, const impq& delta) {
unsigned tj = term_index + m_terms_start_index;
auto it = m_ext_vars_to_columns.find(tj);
if (it == m_ext_vars_to_columns.end())
return true;
unsigned j = it->second.internal_j();
unsigned j;
if (m_var_register.external_is_used(tj, j) == false)
return true; // the term is not a column so it has no bounds
auto & slv = m_mpq_lar_core_solver.m_r_solver;
TRACE("cube", tout << "delta = " << delta << std::endl;
m_int_solver->display_column(tout, j); );
@ -2166,7 +2154,7 @@ void lar_solver::update_delta_for_terms(const impq & delta, unsigned j, const ve
for (unsigned i : terms_of_var) {
lar_term & t = *m_terms[i];
auto it = t.m_coeffs.find(j);
unsigned tj = to_var_index(i + m_terms_start_index);
unsigned tj = to_column(i + m_terms_start_index);
TRACE("cube",
tout << "t.apply = " << t.apply(m_mpq_lar_core_solver.m_r_x) << ", m_mpq_lar_core_solver.m_r_x[tj]= " << m_mpq_lar_core_solver.m_r_x[tj];);
TRACE("cube", print_term_as_indices(t, tout);
@ -2199,6 +2187,7 @@ void lar_solver::round_to_integer_solution() {
fill_vars_to_terms(vars_to_terms);
for (unsigned j = 0; j < column_count(); j++) {
if (column_is_int(j)) continue;
if (column_corresponds_to_term(j)) continue;
TRACE("cube", m_int_solver->display_column(tout, j););
impq& v = m_mpq_lar_core_solver.m_r_x[j];
@ -2219,16 +2208,17 @@ void lar_solver::round_to_integer_solution() {
bool lar_solver::get_equality_and_right_side_for_term_on_current_x(unsigned term_index, mpq & rs, constraint_index& ci) const {
unsigned tj = term_index + m_terms_start_index;
auto it = m_ext_vars_to_columns.find(tj);
if (it == m_ext_vars_to_columns.end())
unsigned j;
bool is_int;
if (m_var_register.external_is_used(tj, j, is_int) == false)
return false; // the term does not have bound because it does not correspond to a column
if (!is_int) // todo - allow for the next version of hnf
return false;
unsigned j = it->second.internal_j();
impq term_val;
bool term_val_ready = false;
mpq b;
bool is_strict;
if (has_upper_bound(j, ci, b, is_strict)) {
lp_assert(!is_strict);
if (has_upper_bound(j, ci, b, is_strict) && !is_strict) {
lp_assert(b.is_int());
term_val = terms()[term_index]->apply(m_mpq_lar_core_solver.m_r_x);
term_val_ready = true;
@ -2237,8 +2227,7 @@ bool lar_solver::get_equality_and_right_side_for_term_on_current_x(unsigned term
return true;
}
}
if (has_lower_bound(j, ci, b, is_strict)) {
lp_assert(!is_strict);
if (has_lower_bound(j, ci, b, is_strict) && !is_strict) {
if (!term_val_ready)
term_val = terms()[term_index]->apply(m_mpq_lar_core_solver.m_r_x);
lp_assert(b.is_int());

View file

@ -47,15 +47,6 @@ namespace lp {
class lar_solver : public column_namer {
class ext_var_info {
unsigned m_internal_j; // the internal index
bool m_is_integer;
public:
ext_var_info(unsigned j, var_index internal_j): ext_var_info(j, false) {}
ext_var_info(unsigned j , bool is_int) : m_internal_j(j), m_is_integer(is_int) {}
unsigned internal_j() const { return m_internal_j;}
bool is_integer() const {return m_is_integer;}
};
#if Z3DEBUG_CHECK_UNIQUE_TERMS
struct term_hasher {
std::size_t operator()(const lar_term *t) const
@ -100,8 +91,7 @@ class lar_solver : public column_namer {
lp_settings m_settings;
lp_status m_status;
stacked_value<simplex_strategy_enum> m_simplex_strategy;
std::unordered_map<unsigned, ext_var_info> m_ext_vars_to_columns;
vector<unsigned> m_columns_to_ext_vars_or_term_indices;
var_register m_var_register;
stacked_vector<ul_pair> m_columns_to_ul_pairs;
vector<lar_base_constraint*> m_constraints;
private:
@ -119,8 +109,6 @@ public:
lar_core_solver m_mpq_lar_core_solver;
private:
int_solver * m_int_solver;
bool m_has_int_var;
public :
unsigned terms_start_index() const { return m_terms_start_index; }
@ -311,21 +299,21 @@ public:
vector<constraint_index> get_all_constraint_indices() const;
bool maximize_term_on_tableau(const vector<std::pair<mpq, var_index>> & term,
bool maximize_term_on_tableau(const lar_term & term,
impq &term_max);
bool costs_are_zeros_for_r_solver() const;
bool reduced_costs_are_zeroes_for_r_solver() const;
void set_costs_to_zero(const vector<std::pair<mpq, var_index>> & term);
void set_costs_to_zero(const lar_term & term);
void prepare_costs_for_r_solver(const vector<std::pair<mpq, var_index>> & term);
void prepare_costs_for_r_solver(const lar_term & term);
bool maximize_term_on_corrected_r_solver(const vector<std::pair<mpq, var_index>> & term,
bool maximize_term_on_corrected_r_solver(lar_term & term,
impq &term_max);
// starting from a given feasible state look for the maximum of the term
// return true if found and false if unbounded
lp_status maximize_term(const vector<std::pair<mpq, var_index>> & term,
lp_status maximize_term(unsigned ext_j ,
impq &term_max);
@ -578,7 +566,7 @@ public:
lar_core_solver & get_core_solver() { return m_mpq_lar_core_solver; }
bool column_corresponds_to_term(unsigned) const;
void catch_up_in_updating_int_solver();
var_index to_var_index(unsigned ext_j) const;
var_index to_column(unsigned ext_j) const;
bool tighten_term_bounds_by_delta(unsigned, const impq&);
void round_to_integer_solution();
void update_delta_for_terms(const impq & delta, unsigned j, const vector<unsigned>&);
@ -588,5 +576,6 @@ public:
const vector<unsigned> & r_nbasis() const { return m_mpq_lar_core_solver.r_nbasis(); }
bool get_equality_and_right_side_for_term_on_current_x(unsigned i, mpq &rs, constraint_index& ci) const;
bool remove_from_basis(unsigned);
lar_term get_term_to_maximize(unsigned ext_j) const;
};
}

View file

@ -973,6 +973,8 @@ template <typename T, typename X> bool lp_core_solver_base<T, X>::remove_from_ba
indexed_vector<T> w(m_basis.size()); // the buffer
unsigned i = m_basis_heading[basic_j];
for (auto &c : m_A.m_rows[i]) {
if (c.var() == basic_j)
continue;
if (pivot_column_general(c.var(), basic_j, w))
return true;
}

View file

@ -504,8 +504,7 @@ public:
}
X theta = (this->m_x[leaving] - new_val_for_leaving) / a_ent;
advance_on_entering_and_leaving_tableau_rows(entering, leaving, theta );
X xleaving = this->m_x[leaving];
lp_assert(xleaving == new_val_for_leaving);
lp_assert(this->m_x[leaving] == new_val_for_leaving);
if (this->current_x_is_feasible())
this->set_status(lp_status::OPTIMAL);
}

View file

@ -200,6 +200,7 @@ public:
bool m_int_pivot_fixed_vars_from_basis;
bool m_int_patch_only_integer_values;
unsigned limit_on_rows_for_hnf_cutter;
unsigned limit_on_columns_for_hnf_cutter;
unsigned random_next() { return m_rand(); }
void set_random_seed(unsigned s) { m_rand.set_seed(s); }
@ -219,10 +220,10 @@ public:
reps_in_scaler(20),
pivot_epsilon(0.00000001),
positive_price_epsilon(1e-7),
entering_diag_epsilon ( 1e-8),
c_partial_pivoting ( 10), // this is the constant c from page 410
depth_of_rook_search ( 4),
using_partial_pivoting ( true),
entering_diag_epsilon (1e-8),
c_partial_pivoting (10), // this is the constant c from page 410
depth_of_rook_search (4),
using_partial_pivoting (true),
// dissertation of Achim Koberstein
// if Bx - b is different at any component more that refactor_epsilon then we refactor
refactor_tolerance ( 1e-4),
@ -264,7 +265,8 @@ public:
m_chase_cut_solver_cycle_on_var(10),
m_int_pivot_fixed_vars_from_basis(false),
m_int_patch_only_integer_values(true),
limit_on_rows_for_hnf_cutter(75)
limit_on_rows_for_hnf_cutter(75),
limit_on_columns_for_hnf_cutter(150)
{}
void set_resource_limit(lp_resource_limit& lim) { m_resource_limit = &lim; }

View file

@ -285,7 +285,9 @@ class numeric_traits<lp::numeric_pair<T>> {
return numeric_traits<T>::is_neg(p.x) ||
(numeric_traits<T>::is_zero(p.x) && numeric_traits<T>::is_neg(p.y));
}
static bool is_int(const numeric_pair<T> & p) {
return numeric_traits<T>::is_int(p.x) && numeric_traits<T>::is_zero(p.y);
}
};
template <>

View file

@ -18,33 +18,94 @@ Revision History:
--*/
#pragma once
namespace lp {
class ext_var_info {
unsigned m_internal_j; // the internal index
bool m_is_integer;
public:
ext_var_info() {}
ext_var_info(unsigned j): ext_var_info(j, true) {}
ext_var_info(unsigned j , bool is_int) : m_internal_j(j), m_is_integer(is_int) {}
unsigned internal_j() const { return m_internal_j;}
bool is_integer() const {return m_is_integer;}
};
class var_register {
svector<unsigned> m_local_vars_to_external;
std::unordered_map<unsigned, unsigned> m_external_vars_to_local;
svector<unsigned> m_local_to_external;
std::unordered_map<unsigned, ext_var_info> m_external_to_local;
public:
unsigned add_var(unsigned user_var) {
auto t = m_external_vars_to_local.find(user_var);
if (t != m_external_vars_to_local.end()) {
return t->second;
return add_var(user_var, true);
}
unsigned add_var(unsigned user_var, bool is_int) {
auto t = m_external_to_local.find(user_var);
if (t != m_external_to_local.end()) {
return t->second.internal_j();
}
unsigned ret = size();
m_external_vars_to_local[user_var] = ret;
m_local_vars_to_external.push_back(user_var);
return ret;
unsigned j = size();
m_external_to_local[user_var] = ext_var_info(j, is_int);
m_local_to_external.push_back(user_var);
return j;
}
const svector<unsigned> & vars() const { return m_local_vars_to_external; }
const svector<unsigned> & vars() const { return m_local_to_external; }
unsigned local_var_to_user_var(unsigned local_var) const {
return m_local_vars_to_external[local_var];
unsigned local_to_external(unsigned local_var) const {
return m_local_to_external[local_var];
}
unsigned size() const {
return m_local_vars_to_external.size();
return m_local_to_external.size();
}
void clear() {
m_local_vars_to_external.clear();
m_external_vars_to_local.clear();
m_local_to_external.clear();
m_external_to_local.clear();
}
unsigned external_to_local(unsigned j) const {
auto it = m_external_to_local.find(j);
lp_assert(it != m_external_to_local.end());
return it->second.internal_j();
}
bool external_is_int(unsigned j) const {
auto it = m_external_to_local.find(j);
lp_assert(it != m_external_to_local.end());
return it->second.is_integer();
}
bool external_is_used(unsigned ext_j) const {
auto it = m_external_to_local.find(ext_j);
return it != m_external_to_local.end();
}
bool external_is_used(unsigned ext_j, unsigned & local_j ) const {
auto it = m_external_to_local.find(ext_j);
if ( it == m_external_to_local.end())
return false;
local_j = it->second.internal_j();
return true;
}
bool external_is_used(unsigned ext_j, unsigned & local_j, bool & is_int ) const {
auto it = m_external_to_local.find(ext_j);
if ( it == m_external_to_local.end())
return false;
local_j = it->second.internal_j();
is_int = it->second.is_integer();
return true;
}
bool local_is_int(unsigned j) const {
return external_is_int(m_local_to_external[j]);
}
void shrink(unsigned shrunk_size) {
for (unsigned j = size(); j-- > shrunk_size;) {
m_external_to_local.erase(m_local_to_external[j]);
}
m_local_to_external.resize(shrunk_size);
}
};
}