diff --git a/src/math/lp/core_solver_pretty_printer_def.h b/src/math/lp/core_solver_pretty_printer_def.h index 931333ee725..27aa6c75d72 100644 --- a/src/math/lp/core_solver_pretty_printer_def.h +++ b/src/math/lp/core_solver_pretty_printer_def.h @@ -139,7 +139,7 @@ template void core_solver_pretty_printer::adjust_ case column_type::free_column: break; default: - lp_assert(false); + UNREACHABLE(); break; } } diff --git a/src/math/lp/general_matrix.h b/src/math/lp/general_matrix.h index 749fa30dc43..a4f6693a211 100644 --- a/src/math/lp/general_matrix.h +++ b/src/math/lp/general_matrix.h @@ -114,9 +114,6 @@ class general_matrix { } } - void copy_column_to_indexed_vector(unsigned entering, indexed_vector &w ) const { - lp_assert(false); // not implemented - } general_matrix operator*(const general_matrix & m) const { lp_assert(m.row_count() == column_count()); general_matrix ret(row_count(), m.column_count()); @@ -172,24 +169,7 @@ class general_matrix { return r; } - // bool create_upper_triangle(general_matrix& m, vector& x) { - // for (unsigned i = 1; i < m.row_count(); i++) { - // lp_assert(false); // to be continued - // } - // } - - // bool solve_A_x_equal_b(const general_matrix& m, vector& x, const vector& b) const { - // auto m_copy = m; - // // for square matrices - // lp_assert(row_count() == b.size()); - // lp_assert(x.size() == column_count()); - // lp_assert(row_count() == column_count()); - // x = b; - // create_upper_triangle(copy_of_m, x); - // solve_on_triangle(copy_of_m, x); - // } - // - + void transpose_rows(unsigned i, unsigned l) { lp_assert(i != l); m_row_permutation.transpose_from_right(i, l); diff --git a/src/math/lp/lar_constraints.h b/src/math/lp/lar_constraints.h index 8e6311683b5..f8cffbe5793 100644 --- a/src/math/lp/lar_constraints.h +++ b/src/math/lp/lar_constraints.h @@ -44,7 +44,7 @@ inline std::string lconstraint_kind_string(lconstraint_kind t) { case EQ: return std::string("="); case NE: return std::string("!="); } - lp_unreachable(); + UNREACHABLE(); return std::string(); // it is unreachable } diff --git a/src/math/lp/lar_core_solver.h b/src/math/lp/lar_core_solver.h index 4aa62f0df09..06ef4d50ba3 100644 --- a/src/math/lp/lar_core_solver.h +++ b/src/math/lp/lar_core_solver.h @@ -159,7 +159,7 @@ class lar_core_solver { case column_type::fixed: return true; default: - lp_assert(false); + UNREACHABLE(); } return false; } @@ -174,7 +174,7 @@ class lar_core_solver { case column_type::fixed: return true; default: - lp_assert(false); + UNREACHABLE(); } return false; } diff --git a/src/math/lp/lar_solver.cpp b/src/math/lp/lar_solver.cpp index 1c9f5461beb..9fab5e4371d 100644 --- a/src/math/lp/lar_solver.cpp +++ b/src/math/lp/lar_solver.cpp @@ -466,7 +466,7 @@ namespace lp { default: - lp_unreachable(); // wrong mode + UNREACHABLE(); // wrong mode } return false; } @@ -802,7 +802,7 @@ namespace lp { case GT: return left_side_val > constr.rhs(); case EQ: return left_side_val == constr.rhs(); default: - lp_unreachable(); + UNREACHABLE(); } return false; // it is unreachable } @@ -857,7 +857,7 @@ namespace lp { case EQ: lp_assert(rs != zero_of_type()); break; default: - lp_assert(false); + UNREACHABLE(); return false; } #endif @@ -1816,7 +1816,7 @@ namespace lp { adjust_initial_state_for_tableau_rows(); break; case simplex_strategy_enum::tableau_costs: - lp_assert(false); // not implemented + UNREACHABLE(); // not implemented case simplex_strategy_enum::undecided: adjust_initial_state_for_tableau_rows(); break; @@ -1905,7 +1905,7 @@ namespace lp { } default: - lp_unreachable(); + UNREACHABLE(); } if (m_mpq_lar_core_solver.m_r_upper_bounds[j] == m_mpq_lar_core_solver.m_r_lower_bounds[j]) { m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed; @@ -1959,7 +1959,7 @@ namespace lp { } default: - lp_unreachable(); + UNREACHABLE(); } } @@ -2009,7 +2009,7 @@ namespace lp { } default: - lp_unreachable(); + UNREACHABLE(); } } void lar_solver::update_bound_with_no_ub_no_lb(var_index j, lconstraint_kind kind, const mpq& right_side, constraint_index ci) { @@ -2050,7 +2050,7 @@ namespace lp { } default: - lp_unreachable(); + UNREACHABLE(); } } diff --git a/src/math/lp/lar_term.h b/src/math/lp/lar_term.h index 3ef424e243e..fc73f949f1f 100644 --- a/src/math/lp/lar_term.h +++ b/src/math/lp/lar_term.h @@ -179,7 +179,7 @@ class lar_term { return p.coeff().is_one(); } } - lp_unreachable(); + UNREACHABLE(); return false; } diff --git a/src/math/lp/lia_move.h b/src/math/lp/lia_move.h index 65da5826e7b..ca61d7b7aba 100644 --- a/src/math/lp/lia_move.h +++ b/src/math/lp/lia_move.h @@ -45,7 +45,7 @@ inline std::string lia_move_to_string(lia_move m) { case lia_move::unsat: return "unsat"; default: - lp_assert(false); + UNREACHABLE(); }; return "strange"; } diff --git a/src/math/lp/lp_core_solver_base.cpp b/src/math/lp/lp_core_solver_base.cpp index 61eae344f6c..9a4c375f638 100644 --- a/src/math/lp/lp_core_solver_base.cpp +++ b/src/math/lp/lp_core_solver_base.cpp @@ -27,7 +27,6 @@ template bool lp::lp_core_solver_base >::prin template bool lp::lp_core_solver_base::basis_heading_is_correct() const ; template bool lp::lp_core_solver_base::column_is_dual_feasible(unsigned int) const; template bool lp::lp_core_solver_base::print_statistics_with_iterations_and_nonzeroes_and_cost_and_check_that_the_time_is_over(char const*, std::ostream &); -template void lp::lp_core_solver_base::set_non_basic_x_to_correct_bounds(); template void lp::lp_core_solver_base::add_delta_to_entering(unsigned int, const lp::mpq&); template void lp::lp_core_solver_base >::init(); template void lp::lp_core_solver_base >::init_basis_heading_and_non_basic_columns_vector(); @@ -61,7 +60,6 @@ template bool lp::lp_core_solver_base >::calc template bool lp::lp_core_solver_base::column_is_feasible(unsigned int) const; // template void lp::lp_core_solver_base >::print_linear_combination_of_column_indices(vector, std::allocator > > const&, std::ostream&) const; template bool lp::lp_core_solver_base >::column_is_feasible(unsigned int) const; -template bool lp::lp_core_solver_base >::snap_non_basic_x_to_bound(); template bool lp::lp_core_solver_base>::pivot_column_tableau(unsigned int, unsigned int); template bool lp::lp_core_solver_base::pivot_column_tableau(unsigned int, unsigned int); template void lp::lp_core_solver_base >::transpose_rows_tableau(unsigned int, unsigned int); @@ -70,6 +68,5 @@ template bool lp::lp_core_solver_base::inf_set_is_correct() co template bool lp::lp_core_solver_base >::infeasibility_costs_are_correct() const; template bool lp::lp_core_solver_base::infeasibility_costs_are_correct() const; template bool lp::lp_core_solver_base >::remove_from_basis(unsigned int); -template bool lp::lp_core_solver_base >::remove_from_basis(unsigned int, lp::numeric_pair const&); diff --git a/src/math/lp/lp_core_solver_base.h b/src/math/lp/lp_core_solver_base.h index 2e751330b69..fc167324234 100644 --- a/src/math/lp/lp_core_solver_base.h +++ b/src/math/lp/lp_core_solver_base.h @@ -173,8 +173,6 @@ class lp_core_solver_base { void set_total_iterations(unsigned s) { m_total_iterations = s; } - void set_non_basic_x_to_correct_bounds(); - bool at_bound(const X &x, const X & bound) const { return !below_bound(x, bound) && !above_bound(x, bound); } @@ -300,45 +298,6 @@ class lp_core_solver_base { std::string column_name(unsigned column) const; - bool snap_non_basic_x_to_bound() { - bool ret = false; - for (unsigned j : non_basis()) - ret = snap_column_to_bound(j) || ret; - return ret; - } - - bool snap_column_to_bound(unsigned j) { - switch (m_column_types[j]) { - case column_type::fixed: - if (x_is_at_bound(j)) - break; - m_x[j] = m_lower_bounds[j]; - return true; - case column_type::boxed: - if (x_is_at_bound(j)) - break; // we should preserve x if possible - // snap randomly - if (m_settings.random_next() % 2 == 1) - m_x[j] = m_lower_bounds[j]; - else - m_x[j] = m_upper_bounds[j]; - return true; - case column_type::lower_bound: - if (x_is_at_lower_bound(j)) - break; - m_x[j] = m_lower_bounds[j]; - return true; - case column_type::upper_bound: - if (x_is_at_upper_bound(j)) - break; - m_x[j] = m_upper_bounds[j]; - return true; - default: - break; - } - return false; - } - bool make_column_feasible(unsigned j, numeric_pair & delta) { bool ret = false; lp_assert(m_basis_heading[j] < 0); @@ -384,7 +343,6 @@ class lp_core_solver_base { } bool remove_from_basis(unsigned j); - bool remove_from_basis(unsigned j, const impq&); bool pivot_column_general(unsigned j, unsigned j_basic, indexed_vector & w); void init_basic_part_of_basis_heading() { unsigned m = m_basis.size(); @@ -456,31 +414,6 @@ class lp_core_solver_base { change_basis_unconditionally(leaving, entering); } - bool non_basic_column_is_set_correctly(unsigned j) const { - if (j >= this->m_n()) - return false; - switch (this->m_column_types[j]) { - case column_type::fixed: - case column_type::boxed: - if (!this->x_is_at_bound(j)) - return false; - break; - case column_type::lower_bound: - if (!this->x_is_at_lower_bound(j)) - return false; - break; - case column_type::upper_bound: - if (!this->x_is_at_upper_bound(j)) - return false; - break; - case column_type::free_column: - break; - default: - lp_assert(false); - break; - } - return true; - } bool non_basic_columns_are_set_correctly() const { for (unsigned j : this->m_nbasis) if (!column_is_feasible(j)) { @@ -540,13 +473,11 @@ class lp_core_solver_base { out << "[-oo, oo]"; break; default: - lp_assert(false); + UNREACHABLE(); } return out << "\n"; } - bool column_is_free(unsigned j) const { return this->m_column_types[j] == column_type::free_column; } - bool column_is_fixed(unsigned j) const { return this->m_column_types[j] == column_type::fixed; } @@ -579,16 +510,6 @@ class lp_core_solver_base { } } - // only check for basic columns - bool calc_current_x_is_feasible() const { - unsigned i = this->m_m(); - while (i--) { - if (!column_is_feasible(m_basis[i])) - return false; - } - return true; - } - void transpose_rows_tableau(unsigned i, unsigned ii); void pivot_to_reduced_costs_tableau(unsigned i, unsigned j); diff --git a/src/math/lp/lp_core_solver_base_def.h b/src/math/lp/lp_core_solver_base_def.h index a1255ade3e2..cb0084a168a 100644 --- a/src/math/lp/lp_core_solver_base_def.h +++ b/src/math/lp/lp_core_solver_base_def.h @@ -166,26 +166,6 @@ print_statistics_with_cost_and_check_that_the_time_is_over(X cost, std::ostream return time_is_over(); } -template void lp_core_solver_base:: -set_non_basic_x_to_correct_bounds() { - for (unsigned j : non_basis()) { - switch (m_column_types[j]) { - case column_type::boxed: - m_x[j] = m_d[j] < 0? m_upper_bounds[j]: m_lower_bounds[j]; - break; - case column_type::lower_bound: - m_x[j] = m_lower_bounds[j]; - lp_assert(column_is_dual_feasible(j)); - break; - case column_type::upper_bound: - m_x[j] = m_upper_bounds[j]; - lp_assert(column_is_dual_feasible(j)); - break; - default: - break; - } - } -} template bool lp_core_solver_base:: column_is_dual_feasible(unsigned j) const { switch (m_column_types[j]) { @@ -201,9 +181,9 @@ column_is_dual_feasible(unsigned j) const { case column_type::free_column: return numeric_traits::is_zero(m_d[j]); default: - lp_unreachable(); + UNREACHABLE(); } - lp_unreachable(); + UNREACHABLE(); return false; } template bool lp_core_solver_base:: @@ -257,7 +237,7 @@ template bool lp_core_solver_base::column_is_feas return true; break; default: - lp_unreachable(); + UNREACHABLE(); } return false; // it is unreachable } @@ -453,18 +433,6 @@ template bool lp_core_solver_base::remove_from_ba return false; } -template bool lp_core_solver_base::remove_from_basis(unsigned basic_j, const impq& val) { - indexed_vector 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; - } - return false; -} - template bool lp_core_solver_base::infeasibility_costs_are_correct() const { @@ -513,7 +481,7 @@ lp_core_solver_base::infeasibility_cost_is_correct_for_column(unsigned j) case column_type::free_column: return is_zero(this->m_costs[j]); default: - lp_assert(false); + UNREACHABLE(); return true; } } diff --git a/src/math/lp/lp_primal_core_solver.h b/src/math/lp/lp_primal_core_solver.h index 628777f3328..641f6be2160 100644 --- a/src/math/lp/lp_primal_core_solver.h +++ b/src/math/lp/lp_primal_core_solver.h @@ -19,704 +19,680 @@ Revision History: --*/ #pragma once -#include -#include -#include -#include -#include -#include "util/vector.h" -#include -#include -#include -#include -#include "math/lp/static_matrix.h" #include "math/lp/core_solver_pretty_printer.h" #include "math/lp/lp_core_solver_base.h" +#include "math/lp/static_matrix.h" #include "math/lp/u_set.h" +#include "util/vector.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include namespace lp { -// This core solver solves (Ax=b, lower_bound_values \leq x \leq upper_bound_values, maximize costs*x ) -// The right side b is given implicitly by x and the basis +// This core solver solves (Ax=b, lower_bound_values \leq x \leq +// upper_bound_values, maximize costs*x ) The right side b is given implicitly +// by x and the basis template -class lp_primal_core_solver:public lp_core_solver_base { +class lp_primal_core_solver : public lp_core_solver_base { public: - int m_sign_of_entering_delta; - vector m_costs_backup; - unsigned m_inf_row_index_for_tableau; - bool m_bland_mode_tableau; - u_set m_left_basis_tableau; - unsigned m_bland_mode_threshold; - unsigned m_left_basis_repeated; - vector m_leaving_candidates; - - std::list m_non_basis_list; - void sort_non_basis(); - int choose_entering_column(unsigned number_of_benefitial_columns_to_go_over); - int choose_entering_column_tableau(); - int choose_entering_column_presize(unsigned number_of_benefitial_columns_to_go_over); - - - bool needs_to_grow(unsigned bj) const { - lp_assert(!this->column_is_feasible(bj)); - switch(this->m_column_types[bj]) { - case column_type::free_column: - return false; - case column_type::fixed: - case column_type::lower_bound: - case column_type::boxed: - return this-> x_below_low_bound(bj); - default: - return false; - } - lp_assert(false); // unreachable - return false; - } - - int inf_sign_of_column(unsigned bj) const { - lp_assert(!this->column_is_feasible(bj)); - switch(this->m_column_types[bj]) { - case column_type::free_column: - return 0; - case column_type::lower_bound: - return 1; - case column_type::fixed: - case column_type::boxed: - return this->x_above_upper_bound(bj)? -1: 1; - default: - return -1; - } - lp_assert(false); // unreachable - return 0; - - } - - - bool monoid_can_decrease(const row_cell & rc) const { - unsigned j = rc.var(); - lp_assert(this->column_is_feasible(j)); - switch (this->m_column_types[j]) { - case column_type::free_column: - return true; - case column_type::fixed: - return false; - case column_type::lower_bound: - if (is_pos(rc.coeff())) { - return this->x_above_lower_bound(j); - } - - return true; - case column_type::upper_bound: - if (is_pos(rc.coeff())) { - return true; - } - - return this->x_below_upper_bound(j); - case column_type::boxed: - if (is_pos(rc.coeff())) { - return this->x_above_lower_bound(j); - } - - return this->x_below_upper_bound(j); - default: - return false; - } - lp_assert(false); // unreachable - return false; - } - - bool monoid_can_increase(const row_cell & rc) const { - unsigned j = rc.var(); - lp_assert(this->column_is_feasible(j)); - switch (this->m_column_types[j]) { - case column_type::free_column: - return true; - case column_type::fixed: - return false; - case column_type::lower_bound: - if (is_neg(rc.coeff())) { - return this->x_above_lower_bound(j); - } - - return true; - case column_type::upper_bound: - if (is_neg(rc.coeff())) { - return true; - } - - return this->x_below_upper_bound(j); - case column_type::boxed: - if (is_neg(rc.coeff())) { - return this->x_above_lower_bound(j); - } - - return this->x_below_upper_bound(j); - default: - return false; - } - lp_assert(false); // unreachable - return false; - } + int m_sign_of_entering_delta; + vector m_costs_backup; + unsigned m_inf_row_index_for_tableau; + bool m_bland_mode_tableau; + u_set m_left_basis_tableau; + unsigned m_bland_mode_threshold; + unsigned m_left_basis_repeated; + vector m_leaving_candidates; + + std::list m_non_basis_list; + void sort_non_basis(); + int choose_entering_column_tableau(); + + bool needs_to_grow(unsigned bj) const { + lp_assert(!this->column_is_feasible(bj)); + switch (this->m_column_types[bj]) { + case column_type::free_column: + return false; + case column_type::fixed: + case column_type::lower_bound: + case column_type::boxed: + return this->x_below_low_bound(bj); + default: + return false; + } + UNREACHABLE(); // unreachable + return false; + } + + int inf_sign_of_column(unsigned bj) const { + lp_assert(!this->column_is_feasible(bj)); + switch (this->m_column_types[bj]) { + case column_type::free_column: + return 0; + case column_type::lower_bound: + return 1; + case column_type::fixed: + case column_type::boxed: + return this->x_above_upper_bound(bj) ? -1 : 1; + default: + return -1; + } + UNREACHABLE(); // unreachable + return 0; + } + + bool monoid_can_decrease(const row_cell &rc) const { + unsigned j = rc.var(); + lp_assert(this->column_is_feasible(j)); + switch (this->m_column_types[j]) { + case column_type::free_column: + return true; + case column_type::fixed: + return false; + case column_type::lower_bound: + if (is_pos(rc.coeff())) { + return this->x_above_lower_bound(j); + } - unsigned get_number_of_basic_vars_that_might_become_inf(unsigned j) const { // consider looking at the signs here: todo - unsigned r = 0; - for (const auto & cc : this->m_A.m_columns[j]) { - unsigned k = this->m_basis[cc.var()]; - if (this->m_column_types[k] != column_type::free_column) - r++; - } - return r; - } + return true; + case column_type::upper_bound: + if (is_pos(rc.coeff())) { + return true; + } - - int find_beneficial_column_in_row_tableau_rows_bland_mode(int i, T & a_ent) { - int j = -1; - unsigned bj = this->m_basis[i]; - bool bj_needs_to_grow = needs_to_grow(bj); - for (const row_cell& rc : this->m_A.m_rows[i]) { - if (rc.var() == bj) - continue; - if (bj_needs_to_grow) { - if (!monoid_can_decrease(rc)) - continue; - } else { - if (!monoid_can_increase(rc)) - continue; - } - if (rc.var() < static_cast(j) ) { - j = rc.var(); - a_ent = rc.coeff(); - } - } - if (j == -1) { - m_inf_row_index_for_tableau = i; - } - - return j; - } - - int find_beneficial_column_in_row_tableau_rows(int i, T & a_ent) { - if (m_bland_mode_tableau) - return find_beneficial_column_in_row_tableau_rows_bland_mode(i, a_ent); - // a short row produces short infeasibility explanation and benefits at least one pivot operation - int choice = -1; - int nchoices = 0; - unsigned num_of_non_free_basics = 1000000; - unsigned len = 100000000; - unsigned bj = this->m_basis[i]; - bool bj_needs_to_grow = needs_to_grow(bj); - for (unsigned k = 0; k < this->m_A.m_rows[i].size(); k++) { - const row_cell& rc = this->m_A.m_rows[i][k]; - unsigned j = rc.var(); - if (j == bj) - continue; - if (bj_needs_to_grow) { - if (!monoid_can_decrease(rc)) - continue; - } else { - if (!monoid_can_increase(rc)) - continue; - } - unsigned damage = get_number_of_basic_vars_that_might_become_inf(j); - if (damage < num_of_non_free_basics) { - num_of_non_free_basics = damage; - len = this->m_A.m_columns[j].size(); - choice = k; - nchoices = 1; - } else if (damage == num_of_non_free_basics && - this->m_A.m_columns[j].size() <= len && (this->m_settings.random_next() % (++nchoices))) { - choice = k; - len = this->m_A.m_columns[j].size(); - } - } - + return this->x_below_upper_bound(j); + case column_type::boxed: + if (is_pos(rc.coeff())) { + return this->x_above_lower_bound(j); + } - if (choice == -1) { - m_inf_row_index_for_tableau = i; - return -1; - } - const row_cell& rc = this->m_A.m_rows[i][choice]; - a_ent = rc.coeff(); - return rc.var(); - } - - bool try_jump_to_another_bound_on_entering(unsigned entering, const X & theta, X & t, bool & unlimited); - bool try_jump_to_another_bound_on_entering_unlimited(unsigned entering, X & t); - int find_leaving_and_t_tableau(unsigned entering, X & t); - - void limit_theta(const X & lim, X & theta, bool & unlimited) { - if (unlimited) { - theta = lim; - unlimited = false; - } else { - theta = std::min(lim, theta); - } - } + return this->x_below_upper_bound(j); + default: + return false; + } + UNREACHABLE(); // unreachable + return false; + } + + bool monoid_can_increase(const row_cell &rc) const { + unsigned j = rc.var(); + lp_assert(this->column_is_feasible(j)); + switch (this->m_column_types[j]) { + case column_type::free_column: + return true; + case column_type::fixed: + return false; + case column_type::lower_bound: + if (is_neg(rc.coeff())) { + return this->x_above_lower_bound(j); + } - void limit_theta_on_basis_column_for_inf_case_m_neg_upper_bound(unsigned j, const T & m, X & theta, bool & unlimited) { - lp_assert(m < 0 && this->m_column_types[j] == column_type::upper_bound); - limit_inf_on_upper_bound_m_neg(m, this->m_x[j], this->m_upper_bounds[j], theta, unlimited); - } + return true; + case column_type::upper_bound: + if (is_neg(rc.coeff())) { + return true; + } + return this->x_below_upper_bound(j); + case column_type::boxed: + if (is_neg(rc.coeff())) { + return this->x_above_lower_bound(j); + } - void limit_theta_on_basis_column_for_inf_case_m_neg_lower_bound(unsigned j, const T & m, X & theta, bool & unlimited) { - lp_assert(m < 0 && this->m_column_types[j] == column_type::lower_bound); - limit_inf_on_bound_m_neg(m, this->m_x[j], this->m_lower_bounds[j], theta, unlimited); + return this->x_below_upper_bound(j); + default: + return false; + } + UNREACHABLE(); // unreachable + return false; + } + + unsigned get_number_of_basic_vars_that_might_become_inf( + unsigned j) const { // consider looking at the signs here: todo + unsigned r = 0; + for (const auto &cc : this->m_A.m_columns[j]) { + unsigned k = this->m_basis[cc.var()]; + if (this->m_column_types[k] != column_type::free_column) + r++; + } + return r; + } + + int find_beneficial_column_in_row_tableau_rows_bland_mode(int i, T &a_ent) { + int j = -1; + unsigned bj = this->m_basis[i]; + bool bj_needs_to_grow = needs_to_grow(bj); + for (const row_cell &rc : this->m_A.m_rows[i]) { + if (rc.var() == bj) + continue; + if (bj_needs_to_grow) { + if (!monoid_can_decrease(rc)) + continue; + } else { + if (!monoid_can_increase(rc)) + continue; + } + if (rc.var() < static_cast(j)) { + j = rc.var(); + a_ent = rc.coeff(); + } } - - - void limit_theta_on_basis_column_for_inf_case_m_pos_lower_bound(unsigned j, const T & m, X & theta, bool & unlimited) { - lp_assert(m > 0 && this->m_column_types[j] == column_type::lower_bound); - limit_inf_on_lower_bound_m_pos(m, this->m_x[j], this->m_lower_bounds[j], theta, unlimited); + if (j == -1) { + m_inf_row_index_for_tableau = i; + } + + return j; + } + + int find_beneficial_column_in_row_tableau_rows(int i, T &a_ent) { + if (m_bland_mode_tableau) + return find_beneficial_column_in_row_tableau_rows_bland_mode(i, a_ent); + // a short row produces short infeasibility explanation and benefits at + // least one pivot operation + int choice = -1; + int nchoices = 0; + unsigned num_of_non_free_basics = 1000000; + unsigned len = 100000000; + unsigned bj = this->m_basis[i]; + bool bj_needs_to_grow = needs_to_grow(bj); + for (unsigned k = 0; k < this->m_A.m_rows[i].size(); k++) { + const row_cell &rc = this->m_A.m_rows[i][k]; + unsigned j = rc.var(); + if (j == bj) + continue; + if (bj_needs_to_grow) { + if (!monoid_can_decrease(rc)) + continue; + } else { + if (!monoid_can_increase(rc)) + continue; + } + unsigned damage = get_number_of_basic_vars_that_might_become_inf(j); + if (damage < num_of_non_free_basics) { + num_of_non_free_basics = damage; + len = this->m_A.m_columns[j].size(); + choice = k; + nchoices = 1; + } else if (damage == num_of_non_free_basics && + this->m_A.m_columns[j].size() <= len && + (this->m_settings.random_next() % (++nchoices))) { + choice = k; + len = this->m_A.m_columns[j].size(); + } } - void limit_theta_on_basis_column_for_inf_case_m_pos_upper_bound(unsigned j, const T & m, X & theta, bool & unlimited) { - lp_assert(m > 0 && this->m_column_types[j] == column_type::upper_bound); - limit_inf_on_bound_m_pos(m, this->m_x[j], this->m_upper_bounds[j], theta, unlimited); - }; - - void get_bound_on_variable_and_update_leaving_precisely(unsigned j, vector & leavings, T m, X & t, T & abs_of_d_of_leaving); - - X get_max_bound(vector & b); + if (choice == -1) { + m_inf_row_index_for_tableau = i; + return -1; + } + const row_cell &rc = this->m_A.m_rows[i][choice]; + a_ent = rc.coeff(); + return rc.var(); + } + + bool try_jump_to_another_bound_on_entering(unsigned entering, const X &theta, + X &t, bool &unlimited); + bool try_jump_to_another_bound_on_entering_unlimited(unsigned entering, X &t); + int find_leaving_and_t_tableau(unsigned entering, X &t); + + void limit_theta(const X &lim, X &theta, bool &unlimited) { + if (unlimited) { + theta = lim; + unlimited = false; + } else { + theta = std::min(lim, theta); + } + } + + void limit_theta_on_basis_column_for_inf_case_m_neg_upper_bound( + unsigned j, const T &m, X &theta, bool &unlimited) { + lp_assert(m < 0 && this->m_column_types[j] == column_type::upper_bound); + limit_inf_on_upper_bound_m_neg(m, this->m_x[j], this->m_upper_bounds[j], + theta, unlimited); + } + + void limit_theta_on_basis_column_for_inf_case_m_neg_lower_bound( + unsigned j, const T &m, X &theta, bool &unlimited) { + lp_assert(m < 0 && this->m_column_types[j] == column_type::lower_bound); + limit_inf_on_bound_m_neg(m, this->m_x[j], this->m_lower_bounds[j], theta, + unlimited); + } + + void limit_theta_on_basis_column_for_inf_case_m_pos_lower_bound( + unsigned j, const T &m, X &theta, bool &unlimited) { + lp_assert(m > 0 && this->m_column_types[j] == column_type::lower_bound); + limit_inf_on_lower_bound_m_pos(m, this->m_x[j], this->m_lower_bounds[j], + theta, unlimited); + } + + void limit_theta_on_basis_column_for_inf_case_m_pos_upper_bound( + unsigned j, const T &m, X &theta, bool &unlimited) { + lp_assert(m > 0 && this->m_column_types[j] == column_type::upper_bound); + limit_inf_on_bound_m_pos(m, this->m_x[j], this->m_upper_bounds[j], theta, + unlimited); + }; + + void get_bound_on_variable_and_update_leaving_precisely( + unsigned j, vector &leavings, T m, X &t, + T &abs_of_d_of_leaving); + + X get_max_bound(vector &b); #ifdef Z3DEBUG - void check_Ax_equal_b(); - void check_the_bounds(); - void check_bound(unsigned i); - void check_correctness(); + void check_Ax_equal_b(); + void check_the_bounds(); + void check_bound(unsigned i); + void check_correctness(); #endif - // from page 183 of Istvan Maros's book - // the basis structures have not changed yet - void update_reduced_costs_from_pivot_row(unsigned entering, unsigned leaving); + // from page 183 of Istvan Maros's book + // the basis structures have not changed yet + void update_reduced_costs_from_pivot_row(unsigned entering, unsigned leaving); - // return 0 if the reduced cost at entering is close enough to the refreshed - // 1 if it is way off, and 2 if it is unprofitable - int refresh_reduced_cost_at_entering_and_check_that_it_is_off(unsigned entering); + // return 0 if the reduced cost at entering is close enough to the refreshed + // 1 if it is way off, and 2 if it is unprofitable + int refresh_reduced_cost_at_entering_and_check_that_it_is_off( + unsigned entering); - void backup_and_normalize_costs(); + void backup_and_normalize_costs(); - void advance_on_entering_and_leaving(int entering, int leaving, X & t); - void advance_on_entering_and_leaving_tableau(int entering, int leaving, X & t); - void advance_on_entering_equal_leaving(int entering, X & t); - void advance_on_entering_equal_leaving_tableau(int entering, X & t); + void advance_on_entering_and_leaving_tableau(int entering, int leaving, X &t); + void advance_on_entering_equal_leaving_tableau(int entering, X &t); - bool need_to_switch_costs() const { - if (this->m_settings.simplex_strategy() == simplex_strategy_enum::tableau_rows) - return false; - // lp_assert(calc_current_x_is_feasible() == current_x_is_feasible()); - return this->current_x_is_feasible() == this->using_infeas_costs(); - } + bool need_to_switch_costs() const { + if (this->m_settings.simplex_strategy() == + simplex_strategy_enum::tableau_rows) + return false; + // lp_assert(calc_current_x_is_feasible() == + // current_x_is_feasible()); + return this->current_x_is_feasible() == this->using_infeas_costs(); + } + void advance_on_entering_tableau(int entering); - void advance_on_entering(int entering); - void advance_on_entering_tableau(int entering); - void advance_on_entering_precise(int entering); - void push_forward_offset_in_non_basis(unsigned & offset_in_nb); - - unsigned get_number_of_non_basic_column_to_try_for_enter(); - - - // returns the number of iterations - unsigned solve(); - - - - void find_feasible_solution(); - - // bool is_tiny() const {return this->m_m < 10 && this->m_n < 20;} - - void one_iteration_tableau(); - - // this version assumes that the leaving already has the right value, and does not update it - void update_x_tableau_rows(unsigned entering, unsigned leaving, const X& delta) { - this->add_delta_to_x(entering, delta); - if (!this->using_infeas_costs()) { - for (const auto & c : this->m_A.m_columns[entering]) { - if (leaving != this->m_basis[c.var()]) { - this->add_delta_to_x_and_track_feasibility(this->m_basis[c.var()], - delta * this->m_A.get_val(c)); - } - } - } else { // using_infeas_costs() == true - lp_assert(this->column_is_feasible(entering)); - lp_assert(this->m_costs[entering] == zero_of_type()); - // m_d[entering] can change because of the cost change for basic columns. - for (const auto & c : this->m_A.m_columns[entering]) { - unsigned j = this->m_basis[c.var()]; - if (j != leaving) - this->add_delta_to_x(j, -delta * this->m_A.get_val(c)); - update_inf_cost_for_column_tableau(j); - if (is_zero(this->m_costs[j])) - this->remove_column_from_inf_set(j); - else - this->insert_column_into_inf_set(j); - } - } - } - - void update_basis_and_x_tableau_rows(int entering, int leaving, X const & tt) { - lp_assert(entering != leaving); - update_x_tableau_rows(entering, leaving, tt); - this->pivot_column_tableau(entering, this->m_basis_heading[leaving]); - this->change_basis(entering, leaving); - } + void push_forward_offset_in_non_basis(unsigned &offset_in_nb); - - void advance_on_entering_and_leaving_tableau_rows(int entering, int leaving, const X &theta ) { - update_basis_and_x_tableau_rows(entering, leaving, theta); - this->track_column_feasibility(entering); - } + unsigned get_number_of_non_basic_column_to_try_for_enter(); - int find_smallest_inf_column() { - int j = -1; - for (unsigned k : this->inf_set()) { - if (k < static_cast(j)) { - j = k; - } - } - return j; - } + // returns the number of iterations + unsigned solve(); - const X& get_val_for_leaving(unsigned j) const { - lp_assert(!this->column_is_feasible(j)); - switch (this->m_column_types[j]) { - case column_type::fixed: - case column_type::upper_bound: - return this->m_upper_bounds[j]; - case column_type::lower_bound: - return this->m_lower_bounds[j]; - break; - case column_type::boxed: - if (this->x_above_upper_bound(j)) - return this->m_upper_bounds[j]; - else - return this->m_lower_bounds[j]; - break; - default: - UNREACHABLE(); - return this->m_lower_bounds[j]; - } - } - - - void one_iteration_tableau_rows() { - int leaving = find_smallest_inf_column(); - if (leaving == -1) { - this->set_status(lp_status::OPTIMAL); - return; - } - - SASSERT(this->column_is_base(leaving)); - - if (!m_bland_mode_tableau) { - if (m_left_basis_tableau.contains(leaving)) { - if (++m_left_basis_repeated > m_bland_mode_threshold) { - m_bland_mode_tableau = true; - } - } else { - m_left_basis_tableau.insert(leaving); - } - } - T a_ent; - int entering = find_beneficial_column_in_row_tableau_rows(this->m_basis_heading[leaving], a_ent); - if (entering == -1) { - this->set_status(lp_status::INFEASIBLE); - return; - } - const X& new_val_for_leaving = get_val_for_leaving(leaving); - X theta = (this->m_x[leaving] - new_val_for_leaving) / a_ent; - this->m_x[leaving] = new_val_for_leaving; - this->remove_column_from_inf_set(leaving); - advance_on_entering_and_leaving_tableau_rows(entering, leaving, theta ); - if (this->current_x_is_feasible()) - this->set_status(lp_status::OPTIMAL); - } - - void decide_on_status_when_cannot_find_entering() { - lp_assert(!need_to_switch_costs()); - this->set_status(this->current_x_is_feasible()? lp_status::OPTIMAL: lp_status::INFEASIBLE); - } + void find_feasible_solution(); - void limit_theta_on_basis_column_for_feas_case_m_neg_no_check(unsigned j, const T & m, X & theta, bool & unlimited) { - lp_assert(m < 0); - limit_theta((this->m_lower_bounds[j] - this->m_x[j]) / m, theta, unlimited); - if (theta < zero_of_type()) theta = zero_of_type(); - } + // bool is_tiny() const {return this->m_m < 10 && this->m_n < 20;} - bool limit_inf_on_bound_m_neg(const T & m, const X & x, const X & bound, X & theta, bool & unlimited) { - // x gets smaller - lp_assert(m < 0); - if (this->below_bound(x, bound)) return false; - if (this->above_bound(x, bound)) { - limit_theta((bound - x) / m, theta, unlimited); - } else { - theta = zero_of_type(); - unlimited = false; - } - return true; - } + void one_iteration_tableau(); - bool limit_inf_on_bound_m_pos(const T & m, const X & x, const X & bound, X & theta, bool & unlimited) { - // x gets larger - lp_assert(m > 0); - if (this->above_bound(x, bound)) return false; - if (this->below_bound(x, bound)) { - limit_theta((bound - x) / m, theta, unlimited); - } else { - theta = zero_of_type(); - unlimited = false; + // this version assumes that the leaving already has the right value, and does + // not update it + void update_x_tableau_rows(unsigned entering, unsigned leaving, + const X &delta) { + this->add_delta_to_x(entering, delta); + if (!this->using_infeas_costs()) { + for (const auto &c : this->m_A.m_columns[entering]) { + if (leaving != this->m_basis[c.var()]) { + this->add_delta_to_x_and_track_feasibility( + this->m_basis[c.var()], -delta * this->m_A.get_val(c)); } - - return true; + } + } else { // using_infeas_costs() == true + lp_assert(this->column_is_feasible(entering)); + lp_assert(this->m_costs[entering] == zero_of_type()); + // m_d[entering] can change because of the cost change for basic columns. + for (const auto &c : this->m_A.m_columns[entering]) { + unsigned j = this->m_basis[c.var()]; + if (j != leaving) + this->add_delta_to_x(j, -delta * this->m_A.get_val(c)); + update_inf_cost_for_column_tableau(j); + if (is_zero(this->m_costs[j])) + this->remove_column_from_inf_set(j); + else + this->insert_column_into_inf_set(j); + } } - - void limit_inf_on_lower_bound_m_pos(const T & m, const X & x, const X & bound, X & theta, bool & unlimited) { - // x gets larger - lp_assert(m > 0); - if (this->below_bound(x, bound)) { - limit_theta((bound - x) / m, theta, unlimited); + } + + void update_basis_and_x_tableau_rows(int entering, int leaving, X const &tt) { + lp_assert(entering != leaving); + update_x_tableau_rows(entering, leaving, tt); + this->pivot_column_tableau(entering, this->m_basis_heading[leaving]); + this->change_basis(entering, leaving); + } + + void advance_on_entering_and_leaving_tableau_rows(int entering, int leaving, + const X &theta) { + update_basis_and_x_tableau_rows(entering, leaving, theta); + this->track_column_feasibility(entering); + } + + int find_smallest_inf_column() { + int j = -1; + for (unsigned k : this->inf_set()) { + if (k < static_cast(j)) { + j = k; } - } - - void limit_inf_on_upper_bound_m_neg(const T & m, const X & x, const X & bound, X & theta, bool & unlimited) { - // x gets smaller - lp_assert(m < 0); - if (this->above_bound(x, bound)) { - limit_theta((bound - x) / m, theta, unlimited); + return j; + } + + const X &get_val_for_leaving(unsigned j) const { + lp_assert(!this->column_is_feasible(j)); + switch (this->m_column_types[j]) { + case column_type::fixed: + case column_type::upper_bound: + return this->m_upper_bounds[j]; + case column_type::lower_bound: + return this->m_lower_bounds[j]; + break; + case column_type::boxed: + if (this->x_above_upper_bound(j)) + return this->m_upper_bounds[j]; + else + return this->m_lower_bounds[j]; + break; + default: + UNREACHABLE(); + return this->m_lower_bounds[j]; + } + } + + void one_iteration_tableau_rows() { + int leaving = find_smallest_inf_column(); + if (leaving == -1) { + this->set_status(lp_status::OPTIMAL); + return; + } + + SASSERT(this->column_is_base(leaving)); + + if (!m_bland_mode_tableau) { + if (m_left_basis_tableau.contains(leaving)) { + if (++m_left_basis_repeated > m_bland_mode_threshold) { + m_bland_mode_tableau = true; } + } else { + m_left_basis_tableau.insert(leaving); + } } - - void limit_theta_on_basis_column_for_inf_case_m_pos_boxed(unsigned j, const T & m, X & theta, bool & unlimited) { - const X & x = this->m_x[j]; - const X & lbound = this->m_lower_bounds[j]; - - if (this->below_bound(x, lbound)) { - limit_theta((lbound - x) / m, theta, unlimited); + T a_ent; + int entering = find_beneficial_column_in_row_tableau_rows( + this->m_basis_heading[leaving], a_ent); + if (entering == -1) { + this->set_status(lp_status::INFEASIBLE); + return; + } + const X &new_val_for_leaving = get_val_for_leaving(leaving); + X theta = (this->m_x[leaving] - new_val_for_leaving) / a_ent; + this->m_x[leaving] = new_val_for_leaving; + this->remove_column_from_inf_set(leaving); + advance_on_entering_and_leaving_tableau_rows(entering, leaving, theta); + if (this->current_x_is_feasible()) + this->set_status(lp_status::OPTIMAL); + } + + void decide_on_status_when_cannot_find_entering() { + lp_assert(!need_to_switch_costs()); + this->set_status(this->current_x_is_feasible() ? lp_status::OPTIMAL + : lp_status::INFEASIBLE); + } + + void limit_theta_on_basis_column_for_feas_case_m_neg_no_check( + unsigned j, const T &m, X &theta, bool &unlimited) { + lp_assert(m < 0); + limit_theta((this->m_lower_bounds[j] - this->m_x[j]) / m, theta, unlimited); + if (theta < zero_of_type()) + theta = zero_of_type(); + } + + bool limit_inf_on_bound_m_neg(const T &m, const X &x, const X &bound, + X &theta, bool &unlimited) { + // x gets smaller + lp_assert(m < 0); + if (this->below_bound(x, bound)) + return false; + if (this->above_bound(x, bound)) { + limit_theta((bound - x) / m, theta, unlimited); + } else { + theta = zero_of_type(); + unlimited = false; + } + return true; + } + + bool limit_inf_on_bound_m_pos(const T &m, const X &x, const X &bound, + X &theta, bool &unlimited) { + // x gets larger + lp_assert(m > 0); + if (this->above_bound(x, bound)) + return false; + if (this->below_bound(x, bound)) { + limit_theta((bound - x) / m, theta, unlimited); + } else { + theta = zero_of_type(); + unlimited = false; + } + + return true; + } + + void limit_inf_on_lower_bound_m_pos(const T &m, const X &x, const X &bound, + X &theta, bool &unlimited) { + // x gets larger + lp_assert(m > 0); + if (this->below_bound(x, bound)) { + limit_theta((bound - x) / m, theta, unlimited); + } + } + + void limit_inf_on_upper_bound_m_neg(const T &m, const X &x, const X &bound, + X &theta, bool &unlimited) { + // x gets smaller + lp_assert(m < 0); + if (this->above_bound(x, bound)) { + limit_theta((bound - x) / m, theta, unlimited); + } + } + + void limit_theta_on_basis_column_for_inf_case_m_pos_boxed(unsigned j, + const T &m, + X &theta, + bool &unlimited) { + const X &x = this->m_x[j]; + const X &lbound = this->m_lower_bounds[j]; + + if (this->below_bound(x, lbound)) { + limit_theta((lbound - x) / m, theta, unlimited); + } else { + const X &ubound = this->m_upper_bounds[j]; + if (this->below_bound(x, ubound)) { + limit_theta((ubound - x) / m, theta, unlimited); + } else if (!this->above_bound(x, ubound)) { + theta = zero_of_type(); + unlimited = false; + } + } + } + + void limit_theta_on_basis_column_for_inf_case_m_neg_boxed(unsigned j, + const T &m, + X &theta, + bool &unlimited) { + // lp_assert(m < 0 && this->m_column_type[j] == column_type::boxed); + const X &x = this->m_x[j]; + const X &ubound = this->m_upper_bounds[j]; + if (this->above_bound(x, ubound)) { + limit_theta((ubound - x) / m, theta, unlimited); + } else { + const X &lbound = this->m_lower_bounds[j]; + if (this->above_bound(x, lbound)) { + limit_theta((lbound - x) / m, theta, unlimited); + } else if (!this->below_bound(x, lbound)) { + theta = zero_of_type(); + unlimited = false; + } + } + } + + void limit_theta_on_basis_column_for_feas_case_m_pos_no_check( + unsigned j, const T &m, X &theta, bool &unlimited) { + lp_assert(m > 0); + limit_theta((this->m_upper_bounds[j] - this->m_x[j]) / m, theta, unlimited); + if (theta < zero_of_type()) { + theta = zero_of_type(); + } + } + + // j is a basic column or the entering, in any case x[j] has to stay feasible. + // m is the multiplier. updating t in a way that holds the following + // x[j] + t * m >= this->m_lower_bounds[j]( if m < 0 ) + // or + // x[j] + t * m <= this->m_upper_bounds[j] ( if m > 0) + void limit_theta_on_basis_column(unsigned j, T m, X &theta, bool &unlimited) { + switch (this->m_column_types[j]) { + case column_type::free_column: + break; + case column_type::upper_bound: + if (this->current_x_is_feasible()) { + if (m > 0) + limit_theta_on_basis_column_for_feas_case_m_pos_no_check(j, m, theta, + unlimited); + } else { // inside of feasibility_loop + if (m > 0) + limit_theta_on_basis_column_for_inf_case_m_pos_upper_bound( + j, m, theta, unlimited); + else + limit_theta_on_basis_column_for_inf_case_m_neg_upper_bound( + j, m, theta, unlimited); + } + break; + case column_type::lower_bound: + if (this->current_x_is_feasible()) { + if (m < 0) + limit_theta_on_basis_column_for_feas_case_m_neg_no_check(j, m, theta, + unlimited); + } else { + if (m < 0) + limit_theta_on_basis_column_for_inf_case_m_neg_lower_bound( + j, m, theta, unlimited); + else + limit_theta_on_basis_column_for_inf_case_m_pos_lower_bound( + j, m, theta, unlimited); + } + break; + // case fixed: + // if (get_this->current_x_is_feasible()) { + // theta = zero_of_type(); + // break; + // } + // if (m < 0) + // limit_theta_on_basis_column_for_inf_case_m_neg_fixed(j, m, + // theta); + // else + // limit_theta_on_basis_column_for_inf_case_m_pos_fixed(j, m, + // theta); + // break; + case column_type::fixed: + case column_type::boxed: + if (this->current_x_is_feasible()) { + if (m > 0) { + limit_theta_on_basis_column_for_feas_case_m_pos_no_check(j, m, theta, + unlimited); } else { - const X & ubound = this->m_upper_bounds[j]; - if (this->below_bound(x, ubound)){ - limit_theta((ubound - x) / m, theta, unlimited); - } else if (!this->above_bound(x, ubound)) { - theta = zero_of_type(); - unlimited = false; - } + limit_theta_on_basis_column_for_feas_case_m_neg_no_check(j, m, theta, + unlimited); } - } - - void limit_theta_on_basis_column_for_inf_case_m_neg_boxed(unsigned j, const T & m, X & theta, bool & unlimited) { - // lp_assert(m < 0 && this->m_column_type[j] == column_type::boxed); - const X & x = this->m_x[j]; - const X & ubound = this->m_upper_bounds[j]; - if (this->above_bound(x, ubound)) { - limit_theta((ubound - x) / m, theta, unlimited); + } else { + if (m > 0) { + limit_theta_on_basis_column_for_inf_case_m_pos_boxed(j, m, theta, + unlimited); } else { - const X & lbound = this->m_lower_bounds[j]; - if (this->above_bound(x, lbound)){ - limit_theta((lbound - x) / m, theta, unlimited); - } else if (!this->below_bound(x, lbound)) { - theta = zero_of_type(); - unlimited = false; - } + limit_theta_on_basis_column_for_inf_case_m_neg_boxed(j, m, theta, + unlimited); } + } + + break; + default: + UNREACHABLE(); } - void limit_theta_on_basis_column_for_feas_case_m_pos(unsigned j, const T & m, X & theta, bool & unlimited) { - lp_assert(m > 0); - if (this->below_bound(this->m_x[j], this->m_upper_bounds[j])) { - limit_theta((this->m_upper_bounds[j] - this->m_x[j]) / m, theta, unlimited); - if (theta < zero_of_type()) { - theta = zero_of_type(); - unlimited = false; - } - } + if (!unlimited && theta < zero_of_type()) { + theta = zero_of_type(); } + } - void limit_theta_on_basis_column_for_feas_case_m_pos_no_check(unsigned j, const T & m, X & theta, bool & unlimited ) { - lp_assert(m > 0); - limit_theta( (this->m_upper_bounds[j] - this->m_x[j]) / m, theta, unlimited); - if (theta < zero_of_type()) { - theta = zero_of_type(); - } - } + bool column_is_benefitial_for_entering_basis(unsigned j) const; + void init_infeasibility_costs(); - // j is a basic column or the entering, in any case x[j] has to stay feasible. - // m is the multiplier. updating t in a way that holds the following - // x[j] + t * m >= this->m_lower_bounds[j]( if m < 0 ) - // or - // x[j] + t * m <= this->m_upper_bounds[j] ( if m > 0) - void limit_theta_on_basis_column(unsigned j, T m, X & theta, bool & unlimited) { - switch (this->m_column_types[j]) { - case column_type::free_column: break; - case column_type::upper_bound: - if (this->current_x_is_feasible()) { - if (m > 0) - limit_theta_on_basis_column_for_feas_case_m_pos_no_check(j, m, theta, unlimited); - } else { // inside of feasibility_loop - if (m > 0) - limit_theta_on_basis_column_for_inf_case_m_pos_upper_bound(j, m, theta, unlimited); - else - limit_theta_on_basis_column_for_inf_case_m_neg_upper_bound(j, m, theta, unlimited); - } - break; - case column_type::lower_bound: - if (this->current_x_is_feasible()) { - if (m < 0) - limit_theta_on_basis_column_for_feas_case_m_neg_no_check(j, m, theta, unlimited); - } else { - if (m < 0) - limit_theta_on_basis_column_for_inf_case_m_neg_lower_bound(j, m, theta, unlimited); - else - limit_theta_on_basis_column_for_inf_case_m_pos_lower_bound(j, m, theta, unlimited); - } - break; - // case fixed: - // if (get_this->current_x_is_feasible()) { - // theta = zero_of_type(); - // break; - // } - // if (m < 0) - // limit_theta_on_basis_column_for_inf_case_m_neg_fixed(j, m, theta); - // else - // limit_theta_on_basis_column_for_inf_case_m_pos_fixed(j, m, theta); - // break; - case column_type::fixed: - case column_type::boxed: - if (this->current_x_is_feasible()) { - if (m > 0) { - limit_theta_on_basis_column_for_feas_case_m_pos_no_check(j, m, theta, unlimited); - } else { - limit_theta_on_basis_column_for_feas_case_m_neg_no_check(j, m, theta, unlimited); - } - } else { - if (m > 0) { - limit_theta_on_basis_column_for_inf_case_m_pos_boxed(j, m, theta, unlimited); - } else { - limit_theta_on_basis_column_for_inf_case_m_neg_boxed(j, m, theta, unlimited); - } - } - - break; - default: - lp_unreachable(); - } - if (!unlimited && theta < zero_of_type()) { - theta = zero_of_type(); - } - } + void init_infeasibility_cost_for_column(unsigned j); + T get_infeasibility_cost_for_column(unsigned j) const; + void init_infeasibility_costs_for_changed_basis_only(); - - bool column_is_benefitial_for_entering_basis(unsigned j) const; - bool column_is_benefitial_for_entering_basis_precise(unsigned j) const; - bool can_enter_basis(unsigned j); - void init_infeasibility_costs(); - - void init_infeasibility_cost_for_column(unsigned j); - T get_infeasibility_cost_for_column(unsigned j) const; - void init_infeasibility_costs_for_changed_basis_only(); - - void print_column(unsigned j, std::ostream & out); - - void print_bound_info_and_x(unsigned j, std::ostream & out); - - bool basis_column_is_set_correctly(unsigned j) const { - return this->m_A.m_columns[j].size() == 1; - - } - - bool basis_columns_are_set_correctly() const { - for (unsigned j : this->m_basis) - if(!basis_column_is_set_correctly(j)) - return false; - - return this->m_basis_heading.size() == this->m_A.column_count() && this->m_basis.size() == this->m_A.row_count(); - } + void print_column(unsigned j, std::ostream &out); - void init_run_tableau(); - void update_x_tableau(unsigned entering, const X & delta); - void update_inf_cost_for_column_tableau(unsigned j); - -// the delta is between the old and the new cost (old - new) - void update_reduced_cost_for_basic_column_cost_change(const T & delta, unsigned j) { - lp_assert(this->m_basis_heading[j] >= 0); - unsigned i = static_cast(this->m_basis_heading[j]); - for (const row_cell & rc : this->m_A.m_rows[i]) { - unsigned k = rc.var(); - if (k == j) - continue; - this->m_d[k] += delta * rc.coeff(); - } - } - - bool update_basis_and_x_tableau(int entering, int leaving, X const & tt); - void init_reduced_costs_tableau(); - void init_tableau_rows() { - m_bland_mode_tableau = false; - m_left_basis_tableau.clear(); - m_left_basis_tableau.resize(this->m_A.column_count()); - m_left_basis_repeated = 0; - } -// stage1 constructor - lp_primal_core_solver(static_matrix & A, - vector & b, // the right side vector - vector & x, // the number of elements in x needs to be at least as large as the number of columns in A - vector & basis, - vector & nbasis, - vector & heading, - vector & costs, - const vector & column_type_array, - const vector & lower_bound_values, - const vector & upper_bound_values, - lp_settings & settings, - const column_namer& column_names): - lp_core_solver_base(A, // b, - basis, - nbasis, - heading, - x, - costs, - settings, - column_names, - column_type_array, - lower_bound_values, - upper_bound_values), - m_bland_mode_threshold(1000) { - this->set_status(lp_status::UNKNOWN); - } + void print_bound_info_and_x(unsigned j, std::ostream &out); - + bool basis_column_is_set_correctly(unsigned j) const { + return this->m_A.m_columns[j].size() == 1; + } - bool initial_x_is_correct() { - std::set basis_set; - for (unsigned i = 0; i < this->m_A.row_count(); i++) { - basis_set.insert(this->m_basis[i]); - } - for (unsigned j = 0; j < this->m_n(); j++) { - if (this->column_has_lower_bound(j) && this->m_x[j] < numeric_traits::zero()) { - LP_OUT(this->m_settings, "low bound for variable " << j << " does not hold: this->m_x[" << j << "] = " << this->m_x[j] << " is negative " << std::endl); - return false; - } - - if (this->column_has_upper_bound(j) && this->m_x[j] > this->m_upper_bounds[j]) { - LP_OUT(this->m_settings, "upper bound for " << j << " does not hold: " << this->m_upper_bounds[j] << ">" << this->m_x[j] << std::endl); - return false; - } - - if (basis_set.find(j) != basis_set.end()) continue; - if (this->m_column_types[j] == column_type::lower_bound) { - if (numeric_traits::zero() != this->m_x[j]) { - LP_OUT(this->m_settings, "only low bound is set for " << j << " but low bound value " << numeric_traits::zero() << " is not equal to " << this->m_x[j] << std::endl); - return false; - } - } - if (this->m_column_types[j] == column_type::boxed) { - if (this->m_upper_bounds[j] != this->m_x[j] && !numeric_traits::is_zero(this->m_x[j])) { - return false; - } - } - } - return true; - } + bool basis_columns_are_set_correctly() const { + for (unsigned j : this->m_basis) + if (!basis_column_is_set_correctly(j)) + return false; + return this->m_basis_heading.size() == this->m_A.column_count() && + this->m_basis.size() == this->m_A.row_count(); + } + + void init_run_tableau(); + void update_x_tableau(unsigned entering, const X &delta); + void update_inf_cost_for_column_tableau(unsigned j); + + // the delta is between the old and the new cost (old - new) + void update_reduced_cost_for_basic_column_cost_change(const T &delta, + unsigned j) { + lp_assert(this->m_basis_heading[j] >= 0); + unsigned i = static_cast(this->m_basis_heading[j]); + for (const row_cell &rc : this->m_A.m_rows[i]) { + unsigned k = rc.var(); + if (k == j) + continue; + this->m_d[k] += delta * rc.coeff(); + } + } + + bool update_basis_and_x_tableau(int entering, int leaving, X const &tt); + void init_reduced_costs_tableau(); + void init_tableau_rows() { + m_bland_mode_tableau = false; + m_left_basis_tableau.clear(); + m_left_basis_tableau.resize(this->m_A.column_count()); + m_left_basis_repeated = 0; + } + // stage1 constructor + lp_primal_core_solver( + static_matrix &A, + vector &b, // the right side vector + vector &x, // the number of elements in x needs to be at least as large + // as the number of columns in A + vector &basis, vector &nbasis, vector &heading, + vector &costs, const vector &column_type_array, + const vector &lower_bound_values, const vector &upper_bound_values, + lp_settings &settings, const column_namer &column_names) + : lp_core_solver_base(A, // b, + basis, nbasis, heading, x, costs, settings, + column_names, column_type_array, + lower_bound_values, upper_bound_values), + m_bland_mode_threshold(1000) { + this->set_status(lp_status::UNKNOWN); + } - friend core_solver_pretty_printer; + friend core_solver_pretty_printer; }; -} +} // namespace lp diff --git a/src/math/lp/lp_primal_core_solver_def.h b/src/math/lp/lp_primal_core_solver_def.h index cc8ad88b392..4ee8305fa1a 100644 --- a/src/math/lp/lp_primal_core_solver_def.h +++ b/src/math/lp/lp_primal_core_solver_def.h @@ -53,10 +53,6 @@ void lp_primal_core_solver::sort_non_basis() { template bool lp_primal_core_solver::column_is_benefitial_for_entering_basis(unsigned j) const { - return column_is_benefitial_for_entering_basis_precise(j); -} -template -bool lp_primal_core_solver::column_is_benefitial_for_entering_basis_precise(unsigned j) const { const T& dj = this->m_d[j]; TRACE("lar_solver", tout << "dj=" << dj << "\n";); switch (this->m_column_types[j]) { @@ -88,56 +84,12 @@ bool lp_primal_core_solver::column_is_benefitial_for_entering_basis_precis } break; default: - lp_unreachable(); + UNREACHABLE(); break; } return false; } -template -int lp_primal_core_solver::choose_entering_column_presize(unsigned number_of_benefitial_columns_to_go_over) { // at this moment m_y = cB * B(-1) - if (number_of_benefitial_columns_to_go_over == 0) - return -1; - if (this->m_basis_sort_counter == 0) { - sort_non_basis(); - this->m_basis_sort_counter = 20; - } - else { - this->m_basis_sort_counter--; - } - unsigned j_nz = this->m_m() + 1; // this number is greater than the max column size - std::list::iterator entering_iter = m_non_basis_list.end(); - for (auto non_basis_iter = m_non_basis_list.begin(); number_of_benefitial_columns_to_go_over && non_basis_iter != m_non_basis_list.end(); ++non_basis_iter) { - unsigned j = *non_basis_iter; - if (!column_is_benefitial_for_entering_basis(j)) - continue; - - // if we are here then j is a candidate to enter the basis - unsigned t = this->m_columns_nz[j]; - if (t < j_nz) { - j_nz = t; - entering_iter = non_basis_iter; - if (number_of_benefitial_columns_to_go_over) - number_of_benefitial_columns_to_go_over--; - } else if (t == j_nz && this->m_settings.random_next() % 2 == 0) { - entering_iter = non_basis_iter; - } - }// while (number_of_benefitial_columns_to_go_over && initial_offset_in_non_basis != offset_in_nb); - if (entering_iter == m_non_basis_list.end()) - return -1; - unsigned entering = *entering_iter; - m_sign_of_entering_delta = this->m_d[entering] > 0 ? 1 : -1; - m_non_basis_list.erase(entering_iter); - m_non_basis_list.push_back(entering); - return entering; -} - - -template -int lp_primal_core_solver::choose_entering_column(unsigned number_of_benefitial_columns_to_go_over) { // at this moment m_y = cB * B(-1) - return choose_entering_column_presize(number_of_benefitial_columns_to_go_over); -} - template bool lp_primal_core_solver::try_jump_to_another_bound_on_entering(unsigned entering, const X & theta, X & t, @@ -278,24 +230,6 @@ template void lp_primal_core_solver::backup_an -template -void lp_primal_core_solver::advance_on_entering_equal_leaving(int entering, X & t) { - -} - -template void lp_primal_core_solver::advance_on_entering_and_leaving(int entering, int leaving, X & t) { - -} - - -template void lp_primal_core_solver::advance_on_entering_precise(int entering) { - lp_assert(false); -} - -template void lp_primal_core_solver::advance_on_entering(int entering) { - lp_assert(false); -} - template void lp_primal_core_solver::push_forward_offset_in_non_basis(unsigned & offset_in_nb) { if (++offset_in_nb == this->m_nbasis.size()) offset_in_nb = 0; @@ -377,7 +311,7 @@ lp_primal_core_solver::get_infeasibility_cost_for_column(unsigned j) const ret = numeric_traits::zero(); break; default: - lp_assert(false); + UNREACHABLE(); ret = numeric_traits::zero(); // does not matter break; } @@ -427,7 +361,7 @@ lp_primal_core_solver::init_infeasibility_cost_for_column(unsigned j) { this->m_costs[j] = numeric_traits::zero(); break; default: - lp_assert(false); + UNREACHABLE(); break; } @@ -458,7 +392,7 @@ template void lp_primal_core_solver::print_column out << "( _" << this->m_x[j] << "_)" << std::endl; break; default: - lp_unreachable(); + UNREACHABLE(); } } @@ -480,7 +414,7 @@ template void lp_primal_core_solver::print_bound_ out << "inf, inf" << std::endl; break; default: - lp_assert(false); + UNREACHABLE(); break; } } diff --git a/src/math/lp/lp_primal_core_solver_tableau_def.h b/src/math/lp/lp_primal_core_solver_tableau_def.h index c7b604b9553..241164c023f 100644 --- a/src/math/lp/lp_primal_core_solver_tableau_def.h +++ b/src/math/lp/lp_primal_core_solver_tableau_def.h @@ -122,7 +122,7 @@ unsigned lp_primal_core_solver::solve() { } break; case lp_status::TENTATIVE_UNBOUNDED: - lp_assert(false); + UNREACHABLE(); break; case lp_status::UNBOUNDED: if (this->current_x_is_infeasible()) { @@ -132,7 +132,7 @@ unsigned lp_primal_core_solver::solve() { break; case lp_status::UNSTABLE: - lp_assert(false); + UNREACHABLE(); break; default: diff --git a/src/math/lp/lp_settings_def.h b/src/math/lp/lp_settings_def.h index a439466d18c..a19558949c3 100644 --- a/src/math/lp/lp_settings_def.h +++ b/src/math/lp/lp_settings_def.h @@ -31,7 +31,7 @@ std::string column_type_to_string(column_type t) { case column_type::lower_bound: return "lower_bound"; case column_type::upper_bound: return "upper_bound"; case column_type::free_column: return "free_column"; - default: lp_unreachable(); + default: UNREACHABLE(); } return "unknown"; // it is unreachable } @@ -50,7 +50,7 @@ const char* lp_status_to_string(lp_status status) { case lp_status::UNSTABLE: return "UNSTABLE"; case lp_status::CANCELLED: return "CANCELLED"; default: - lp_unreachable(); + UNREACHABLE(); } return "UNKNOWN"; // it is unreachable } @@ -63,7 +63,7 @@ lp_status lp_status_from_string(std::string status) { if (status == "FEASIBLE") return lp_status::FEASIBLE; if (status == "TIME_EXHAUSTED") return lp_status::TIME_EXHAUSTED; if (status == "EMPTY") return lp_status::EMPTY; - lp_unreachable(); + UNREACHABLE(); return lp_status::UNKNOWN; // it is unreachable } diff --git a/src/math/lp/lp_utils.h b/src/math/lp/lp_utils.h index ad5ba380d91..3c1383cb39e 100644 --- a/src/math/lp/lp_utils.h +++ b/src/math/lp/lp_utils.h @@ -141,7 +141,6 @@ inline void throw_exception(std::string && str) { typedef z3_exception exception; #define lp_assert(_x_) { SASSERT(_x_); } -inline void lp_unreachable() { lp_assert(false); } template inline X zero_of_type() { return numeric_traits::zero(); } template inline X one_of_type() { return numeric_traits::one(); } template inline bool is_zero(const X & v) { return numeric_traits::is_zero(v); } diff --git a/src/math/lp/nra_solver.cpp b/src/math/lp/nra_solver.cpp index 56a84d1f0d7..1f4e0b76abb 100644 --- a/src/math/lp/nra_solver.cpp +++ b/src/math/lp/nra_solver.cpp @@ -171,7 +171,7 @@ struct solver::imp { lit = m_nlsat->mk_ineq_literal(nlsat::atom::kind::EQ, 1, ps, is_even); break; default: - lp_assert(false); // unreachable + UNREACHABLE(); // unreachable } m_nlsat->mk_clause(1, &lit, a); } diff --git a/src/math/lp/numeric_pair.h b/src/math/lp/numeric_pair.h index a64825cd8f1..25127400627 100644 --- a/src/math/lp/numeric_pair.h +++ b/src/math/lp/numeric_pair.h @@ -107,8 +107,8 @@ class numeric_traits { template struct convert_struct { static X convert(const Y & y){ return X(y);} - static bool below_bound_numeric(const X &, const X &, const Y &) { /*lp_unreachable();*/ return false;} - static bool above_bound_numeric(const X &, const X &, const Y &) { /*lp_unreachable();*/ return false; } + static bool below_bound_numeric(const X &, const X &, const Y &) { /*UNREACHABLE();*/ return false;} + static bool above_bound_numeric(const X &, const X &, const Y &) { /*UNREACHABLE();*/ return false; } }; @@ -190,7 +190,7 @@ struct numeric_pair { } numeric_pair operator/(const numeric_pair &) const { - // lp_unreachable(); + // UNREACHABLE(); } @@ -199,7 +199,7 @@ struct numeric_pair { } numeric_pair operator*(const numeric_pair & /*a*/) const { - // lp_unreachable(); + // UNREACHABLE(); } numeric_pair& operator+=(const numeric_pair & a) { @@ -275,14 +275,14 @@ numeric_pair operator/(const numeric_pair & r, const X & a) { return numeric_pair(r.x / a, r.y / a); } -template double get_double(const lp::numeric_pair & ) { /* lp_unreachable(); */ return 0;} +template double get_double(const lp::numeric_pair & ) { /* UNREACHABLE(); */ return 0;} template class numeric_traits> { public: static lp::numeric_pair zero() { return lp::numeric_pair(numeric_traits::zero(), numeric_traits::zero()); } static bool is_zero(const lp::numeric_pair & v) { return numeric_traits::is_zero(v.x) && numeric_traits::is_zero(v.y); } static double get_double(const lp::numeric_pair & v){ return numeric_traits::get_double(v.x); } // just return the double of the first coordinate - static double one() { /*lp_unreachable();*/ return 0;} + static double one() { /*UNREACHABLE();*/ return 0;} static bool is_pos(const numeric_pair &p) { return numeric_traits::is_pos(p.x) || (numeric_traits::is_zero(p.x) && numeric_traits::is_pos(p.y)); diff --git a/src/math/lp/static_matrix.cpp b/src/math/lp/static_matrix.cpp index 28a23b0c394..efb6e07cf75 100644 --- a/src/math/lp/static_matrix.cpp +++ b/src/math/lp/static_matrix.cpp @@ -31,7 +31,6 @@ template std::set> lp::static_matrix::add_column_to_vector(mpq const&, unsigned int, mpq*) const; template void static_matrix::add_columns_at_the_end(unsigned int); template bool static_matrix::is_correct() const; -template void static_matrix::copy_column_to_indexed_vector(unsigned int, indexed_vector&) const; template mpq static_matrix::get_balance() const; template mpq static_matrix::get_elem(unsigned int, unsigned int) const; @@ -47,7 +46,6 @@ template static_matrix::static_matrix(unsigned int, unsigned int); #ifdef Z3DEBUG template bool static_matrix >::is_correct() const; #endif -template void static_matrix >::copy_column_to_indexed_vector(unsigned int, indexed_vector&) const; template mpq static_matrix >::get_elem(unsigned int, unsigned int) const; template void static_matrix >::init_empty_matrix(unsigned int, unsigned int); template void static_matrix >::set(unsigned int, unsigned int, mpq const&); diff --git a/src/math/lp/static_matrix.h b/src/math/lp/static_matrix.h index d7e4370a344..f79ff36ac38 100644 --- a/src/math/lp/static_matrix.h +++ b/src/math/lp/static_matrix.h @@ -168,8 +168,6 @@ class static_matrix std::set> get_domain(); - void copy_column_to_indexed_vector(unsigned j, indexed_vector & v) const; - T get_max_abs_in_row(unsigned row) const; void add_column_to_vector (const T & a, unsigned j, T * v) const { for (const auto & it : m_columns[j]) { @@ -222,7 +220,7 @@ class static_matrix virtual void set_number_of_columns(unsigned /*n*/) { } #endif - T get_max_val_in_row(unsigned /* i */) const { lp_unreachable(); } + T get_max_val_in_row(unsigned /* i */) const { UNREACHABLE(); } T get_balance() const; diff --git a/src/math/lp/static_matrix_def.h b/src/math/lp/static_matrix_def.h index af2eac36017..76c1dec546c 100644 --- a/src/math/lp/static_matrix_def.h +++ b/src/math/lp/static_matrix_def.h @@ -174,14 +174,6 @@ std::set> static_matrix::get_domain() { return ret; } -template void static_matrix::copy_column_to_indexed_vector (unsigned j, indexed_vector & v) const { - lp_assert(j < m_columns.size()); - for (auto & it : m_columns[j]) { - const T& val = get_val(it); - if (!is_zero(val)) - v.set_value(val, it.var()); - } -} template T static_matrix::get_max_abs_in_row(unsigned row) const { T ret = numeric_traits::zero(); diff --git a/src/test/lp/gomory_test.h b/src/test/lp/gomory_test.h index 890ff90e352..c64c0103653 100644 --- a/src/test/lp/gomory_test.h +++ b/src/test/lp/gomory_test.h @@ -130,7 +130,7 @@ struct gomory_test { void report_conflict_from_gomory_cut(mpq &k) { - lp_assert(false); + UNREACHABLE(); } void adjust_term_and_k_for_some_ints_case_gomory(lar_term& t, mpq& k, mpq &lcm_den) { diff --git a/src/test/lp/lp.cpp b/src/test/lp/lp.cpp index 8ad73f0c003..9120d64cfdf 100644 --- a/src/test/lp/lp.cpp +++ b/src/test/lp/lp.cpp @@ -1365,7 +1365,7 @@ void test_gomory_cut_0() { if (j == 2) return zero_of_type(); if (j == 3) return mpq(3); - lp_assert(false); + UNREACHABLE(); return zero_of_type(); }, [](unsigned j) { // at_low_p @@ -1375,7 +1375,7 @@ void test_gomory_cut_0() { return true; if (j == 3) return true; - lp_assert(false); + UNREACHABLE(); return false; }, [](unsigned j) { // at_upper @@ -1385,31 +1385,31 @@ void test_gomory_cut_0() { return true; if (j == 3) return false; - lp_assert(false); + UNREACHABLE(); return false; }, [](unsigned j) { // lower_bound if (j == 1) { - lp_assert(false); //unlimited from below + UNREACHABLE(); //unlimited from below return impq(0); } if (j == 2) return impq(0); if (j == 3) return impq(3); - lp_assert(false); + UNREACHABLE(); return impq(0); }, [](unsigned j) { // upper if (j == 1) { - lp_assert(false); //unlimited from above + UNREACHABLE(); //unlimited from above return impq(0); } if (j == 2) return impq(0); if (j == 3) return impq(10); - lp_assert(false); + UNREACHABLE(); return impq(0); }, [] (unsigned) { return 0; }, @@ -1437,7 +1437,7 @@ void test_gomory_cut_1() { return mpq(4363334, 2730001); if (j == 3) return mpq(1); - lp_assert(false); + UNREACHABLE(); return zero_of_type(); }, [](unsigned j) { // at_low_p @@ -1447,7 +1447,7 @@ void test_gomory_cut_1() { return false; if (j == 3) return true; - lp_assert(false); + UNREACHABLE(); return false; }, [](unsigned j) { // at_upper @@ -1457,19 +1457,19 @@ void test_gomory_cut_1() { return false; if (j == 3) return true; - lp_assert(false); + UNREACHABLE(); return false; }, [](unsigned j) { // lower_bound if (j == 1) { - lp_assert(false); //unlimited from below + UNREACHABLE(); //unlimited from below return impq(0); } if (j == 2) return impq(1); if (j == 3) return impq(1); - lp_assert(false); + UNREACHABLE(); return impq(0); }, [](unsigned j) { // upper @@ -1480,7 +1480,7 @@ void test_gomory_cut_1() { return impq(3333); if (j == 3) return impq(10000); - lp_assert(false); + UNREACHABLE(); return impq(0); }, [] (unsigned) { return 0; }, diff --git a/src/test/lp/smt_reader.h b/src/test/lp/smt_reader.h index 75edb23b987..7843d5714d2 100644 --- a/src/test/lp/smt_reader.h +++ b/src/test/lp/smt_reader.h @@ -272,7 +272,7 @@ namespace lp { } else if (el.m_head == "+") { add_sum(c, el.m_elems); } else { - lp_assert(false); // unexpected input + UNREACHABLE(); // unexpected input } }