diff --git a/src/math/lp/lar_core_solver_def.h b/src/math/lp/lar_core_solver_def.h index e22d77c1b30..550b6fe3674 100644 --- a/src/math/lp/lar_core_solver_def.h +++ b/src/math/lp/lar_core_solver_def.h @@ -40,7 +40,6 @@ void lar_core_solver::prefix_r() { if (m_r_solver.m_settings.simplex_strategy() != simplex_strategy_enum::tableau_rows) { m_r_solver.m_costs.resize(m_r_solver.m_n()); m_r_solver.m_d.resize(m_r_solver.m_n()); - m_r_solver.set_using_infeas_costs(true); } } diff --git a/src/math/lp/lar_solver.cpp b/src/math/lp/lar_solver.cpp index 9fab5e4371d..21eec31d839 100644 --- a/src/math/lp/lar_solver.cpp +++ b/src/math/lp/lar_solver.cpp @@ -191,9 +191,11 @@ namespace lp { stats().m_max_rows = A_r().row_count(); if (strategy_is_undecided()) decide_on_strategy_and_adjust_initial_state(); - + auto strategy_was = settings().simplex_strategy(); + settings().set_simplex_strategy(simplex_strategy_enum::tableau_rows); m_mpq_lar_core_solver.m_r_solver.m_look_for_feasible_solution_only = true; auto ret = solve(); + settings().set_simplex_strategy(strategy_was); return ret; } @@ -355,7 +357,6 @@ namespace lp { m_basic_columns_with_changed_cost.resize(m_mpq_lar_core_solver.m_r_x.size()); move_non_basic_columns_to_bounds(false); auto& rslv = m_mpq_lar_core_solver.m_r_solver; - rslv.set_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()); @@ -490,11 +491,11 @@ namespace lp { lp_status lar_solver::maximize_term(unsigned j_or_term, impq& term_max) { TRACE("lar_solver", print_values(tout);); + lar_term term = get_term_to_maximize(j_or_term); if (term.is_empty()) { return lp_status::UNBOUNDED; } - impq prev_value; auto backup = m_mpq_lar_core_solver.m_r_x; if (m_mpq_lar_core_solver.m_r_solver.calc_current_x_is_feasible_include_non_basis()) { @@ -710,14 +711,6 @@ namespace lp { void lar_solver::update_x_and_inf_costs_for_columns_with_changed_bounds_tableau() { for (auto j : m_columns_with_changed_bounds) update_x_and_inf_costs_for_column_with_changed_bounds(j); - - if (tableau_with_costs()) { - if (m_mpq_lar_core_solver.m_r_solver.using_infeas_costs()) { - for (unsigned j : m_basic_columns_with_changed_cost) - m_mpq_lar_core_solver.m_r_solver.update_inf_cost_for_column_tableau(j); - lp_assert(m_mpq_lar_core_solver.m_r_solver.reduced_costs_are_correct_tableau()); - } - } } @@ -1344,14 +1337,6 @@ namespace lp { for (unsigned j : became_feas) m_mpq_lar_core_solver.m_r_solver.remove_column_from_inf_set(j); - - if (use_tableau_costs()) { - for (unsigned j : became_feas) - m_mpq_lar_core_solver.m_r_solver.update_inf_cost_for_column_tableau(j); - for (unsigned j : basic_columns_with_changed_cost) - m_mpq_lar_core_solver.m_r_solver.update_inf_cost_for_column_tableau(j); - lp_assert(m_mpq_lar_core_solver.m_r_solver.reduced_costs_are_correct_tableau()); - } } bool lar_solver::model_is_int_feasible() const { diff --git a/src/math/lp/lp_core_solver_base.cpp b/src/math/lp/lp_core_solver_base.cpp index 9a4c375f638..f1ae95ea0b9 100644 --- a/src/math/lp/lp_core_solver_base.cpp +++ b/src/math/lp/lp_core_solver_base.cpp @@ -65,8 +65,6 @@ template bool lp::lp_core_solver_base::pivot_column_tableau(un template void lp::lp_core_solver_base >::transpose_rows_tableau(unsigned int, unsigned int); template bool lp::lp_core_solver_base >::inf_set_is_correct() const; template bool lp::lp_core_solver_base::inf_set_is_correct() const; -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); diff --git a/src/math/lp/lp_core_solver_base.h b/src/math/lp/lp_core_solver_base.h index fc167324234..fb0c2850795 100644 --- a/src/math/lp/lp_core_solver_base.h +++ b/src/math/lp/lp_core_solver_base.h @@ -63,15 +63,12 @@ class lp_core_solver_base { bool current_x_is_infeasible() const { return m_inf_set.size() != 0; } private: u_set m_inf_set; - bool m_using_infeas_costs; public: const u_set& inf_set() const { return m_inf_set; } u_set& inf_set() { return m_inf_set; } void inf_set_increase_size_by_one() { m_inf_set.increase_size_by_one(); } bool inf_set_contains(unsigned j) const { return m_inf_set.contains(j); } - unsigned inf_set_size() const { return m_inf_set.size(); } - bool using_infeas_costs() const { return m_using_infeas_costs; } - void set_using_infeas_costs(bool val) { m_using_infeas_costs = val; } + unsigned inf_set_size() const { return m_inf_set.size(); } indexed_vector m_pivot_row; // this is the real pivot row of the simplex tableu static_matrix & m_A; // the matrix A // vector const & m_b; // the right side @@ -198,11 +195,7 @@ class lp_core_solver_base { if (m_settings.simplex_strategy() == simplex_strategy_enum::tableau_rows) return true; CASSERT("check_static_matrix", m_A.is_correct()); - if (m_using_infeas_costs) { - if (infeasibility_costs_are_correct() == false) { - return false; - } - } + unsigned n = m_A.column_count(); for (unsigned j = 0; j < n; j++) { @@ -236,8 +229,6 @@ class lp_core_solver_base { return below_bound(m_x[p], m_lower_bounds[p]); } - bool infeasibility_costs_are_correct() const; - bool infeasibility_cost_is_correct_for_column(unsigned j) const; bool x_above_lower_bound(unsigned p) const { return above_bound(m_x[p], m_lower_bounds[p]); diff --git a/src/math/lp/lp_core_solver_base_def.h b/src/math/lp/lp_core_solver_base_def.h index cb0084a168a..8619c926e8d 100644 --- a/src/math/lp/lp_core_solver_base_def.h +++ b/src/math/lp/lp_core_solver_base_def.h @@ -43,7 +43,6 @@ lp_core_solver_base(static_matrix & A, m_iters_with_no_cost_growing(0), m_status(lp_status::FEASIBLE), m_inf_set(A.column_count()), - m_using_infeas_costs(false), m_pivot_row(A.column_count()), m_A(A), m_basis(basis), @@ -94,8 +93,6 @@ pivot_to_reduced_costs_tableau(unsigned i, unsigned j) { - - // template void lp_core_solver_base:: // update_index_of_ed() { // m_index_of_ed.clear(); @@ -434,57 +431,4 @@ template bool lp_core_solver_base::remove_from_ba } -template bool -lp_core_solver_base::infeasibility_costs_are_correct() const { - if (! this->m_using_infeas_costs) - return true; - lp_assert(costs_on_nbasis_are_zeros()); - for (unsigned j :this->m_basis) { - if (!infeasibility_cost_is_correct_for_column(j)) { - TRACE("lar_solver", tout << "incorrect cost for column " << j << std::endl;); - return false; - } - if (!is_zero(m_d[j])) { - TRACE("lar_solver", tout << "non zero inf cost for basis j = " << j << std::endl;); - return false; - } - } - return true; -} - -template bool -lp_core_solver_base::infeasibility_cost_is_correct_for_column(unsigned j) const { - T r = -one_of_type(); - - switch (this->m_column_types[j]) { - case column_type::fixed: - case column_type::boxed: - if (this->x_above_upper_bound(j)) { - return (this->m_costs[j] == r); - } - if (this->x_below_low_bound(j)) { - return (this->m_costs[j] == -r); - } - return is_zero(this->m_costs[j]); - - case column_type::lower_bound: - if (this->x_below_low_bound(j)) { - return this->m_costs[j] == -r; - } - return is_zero(this->m_costs[j]); - - case column_type::upper_bound: - if (this->x_above_upper_bound(j)) { - return this->m_costs[j] == r; - } - return is_zero(this->m_costs[j]); - case column_type::free_column: - return is_zero(this->m_costs[j]); - default: - UNREACHABLE(); - return true; - } -} - - } diff --git a/src/math/lp/lp_primal_core_solver.cpp b/src/math/lp/lp_primal_core_solver.cpp index 22042668da3..efbfd27e1cc 100644 --- a/src/math/lp/lp_primal_core_solver.cpp +++ b/src/math/lp/lp_primal_core_solver.cpp @@ -33,6 +33,5 @@ template unsigned lp_primal_core_solver::solve(); template unsigned lp_primal_core_solver >::solve(); template bool lp::lp_primal_core_solver::update_basis_and_x_tableau(int, int, lp::mpq const&); template bool lp::lp_primal_core_solver >::update_basis_and_x_tableau(int, int, lp::numeric_pair const&); -template void lp::lp_primal_core_solver >::update_inf_cost_for_column_tableau(unsigned); } diff --git a/src/math/lp/lp_primal_core_solver.h b/src/math/lp/lp_primal_core_solver.h index 641f6be2160..207428985bf 100644 --- a/src/math/lp/lp_primal_core_solver.h +++ b/src/math/lp/lp_primal_core_solver.h @@ -334,28 +334,12 @@ class lp_primal_core_solver : public lp_core_solver_base { 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); - } - } + 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)); + } + } } void update_basis_and_x_tableau_rows(int entering, int leaving, X const &tt) { @@ -437,7 +421,6 @@ class lp_primal_core_solver : public lp_core_solver_base { } 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); } @@ -628,11 +611,6 @@ class lp_primal_core_solver : public lp_core_solver_base { bool column_is_benefitial_for_entering_basis(unsigned j) const; 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); @@ -652,8 +630,6 @@ class lp_primal_core_solver : public lp_core_solver_base { 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) { diff --git a/src/math/lp/lp_primal_core_solver_def.h b/src/math/lp/lp_primal_core_solver_def.h index 4ee8305fa1a..c3c545fdd8a 100644 --- a/src/math/lp/lp_primal_core_solver_def.h +++ b/src/math/lp/lp_primal_core_solver_def.h @@ -258,122 +258,10 @@ template void lp_primal_core_solver::find_feas solve(); } -template -void lp_primal_core_solver::init_infeasibility_costs_for_changed_basis_only() { - for (unsigned i : this->m_ed.m_index) - init_infeasibility_cost_for_column(this->m_basis[i]); - this->set_using_infeas_costs(true); -} - -template -void lp_primal_core_solver::init_infeasibility_costs() { - lp_assert(this->m_x.size() >= this->m_n()); - lp_assert(this->m_column_types.size() >= this->m_n()); - for (unsigned j = this->m_n(); j--;) - init_infeasibility_cost_for_column(j); - this->set_using_infeas_costs(true); -} -template T -lp_primal_core_solver::get_infeasibility_cost_for_column(unsigned j) const { - if (this->m_basis_heading[j] < 0) { - return zero_of_type(); - } - T ret; - // j is a basis column - switch (this->m_column_types[j]) { - case column_type::fixed: - case column_type::boxed: - if (this->x_above_upper_bound(j)) { - ret = 1; - } else if (this->x_below_low_bound(j)) { - ret = -1; - } else { - ret = numeric_traits::zero(); - } - break; - case column_type::lower_bound: - if (this->x_below_low_bound(j)) { - ret = -1; - } else { - ret = numeric_traits::zero(); - } - break; - case column_type::upper_bound: - if (this->x_above_upper_bound(j)) { - ret = 1; - } else { - ret = numeric_traits::zero(); - } - break; - case column_type::free_column: - ret = numeric_traits::zero(); - break; - default: - UNREACHABLE(); - ret = numeric_traits::zero(); // does not matter - break; - } - - ret = - ret; - - return ret; -} -// changed m_inf_set too! -template void -lp_primal_core_solver::init_infeasibility_cost_for_column(unsigned j) { - - if (this->m_basis_heading[j] < 0) { - this->m_costs[j] = numeric_traits::zero(); - this->remove_column_from_inf_set(j); - return; - } - // j is a basis column - switch (this->m_column_types[j]) { - case column_type::fixed: - case column_type::boxed: - if (this->x_above_upper_bound(j)) { - this->m_costs[j] = 1; - } else if (this->x_below_low_bound(j)) { - this->m_costs[j] = -1; - } else { - this->m_costs[j] = numeric_traits::zero(); - } - break; - case column_type::lower_bound: - if (this->x_below_low_bound(j)) { - this->m_costs[j] = -1; - } else { - this->m_costs[j] = numeric_traits::zero(); - } - break; - case column_type::upper_bound: - if (this->x_above_upper_bound(j)) { - this->m_costs[j] = 1; - } else { - this->m_costs[j] = numeric_traits::zero(); - } - break; - case column_type::free_column: - this->m_costs[j] = numeric_traits::zero(); - break; - default: - UNREACHABLE(); - break; - } - - if (numeric_traits::is_zero(this->m_costs[j])) { - this->remove_column_from_inf_set(j); - } else { - this->insert_column_into_inf_set(j); - } - this->m_costs[j] = - this->m_costs[j]; - -} - template void lp_primal_core_solver::print_column(unsigned j, std::ostream & out) { out << this->column_name(j) << " " << j << " " << column_type_to_string(this->m_column_type[j]) << " x = " << this->m_x[j] << " " << "c = " << this->m_costs[j];; 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 241164c023f..898abb1525b 100644 --- a/src/math/lp/lp_primal_core_solver_tableau_def.h +++ b/src/math/lp/lp_primal_core_solver_tableau_def.h @@ -96,7 +96,7 @@ unsigned lp_primal_core_solver::solve() { } do { - if (this->print_statistics_with_iterations_and_nonzeroes_and_cost_and_check_that_the_time_is_over((this->using_infeas_costs()? "inf t" : "feas t"), * this->m_settings.get_message_ostream())) { + if (this->print_statistics_with_iterations_and_nonzeroes_and_cost_and_check_that_the_time_is_over( "feas t", * this->m_settings.get_message_ostream())) { return this->total_iterations(); } if (this->m_settings.use_tableau_rows()) { @@ -107,28 +107,12 @@ unsigned lp_primal_core_solver::solve() { TRACE("lar_solver", tout << "one iteration tableau " << this->get_status() << "\n";); switch (this->get_status()) { case lp_status::OPTIMAL: // check again that we are at optimum - case lp_status::INFEASIBLE: - if (this->m_look_for_feasible_solution_only && this->current_x_is_feasible()) - break; - { // precise case - if ((!this->infeasibility_costs_are_correct())) { - init_reduced_costs_tableau(); // forcing recalc - if (choose_entering_column_tableau() == -1) { - decide_on_status_when_cannot_find_entering(); - break; - } - this->set_status(lp_status::UNKNOWN); - } - } break; case lp_status::TENTATIVE_UNBOUNDED: UNREACHABLE(); break; case lp_status::UNBOUNDED: - if (this->current_x_is_infeasible()) { - init_reduced_costs_tableau(); - this->set_status(lp_status::UNKNOWN); - } + lp_assert (this->current_x_is_feasible()); break; case lp_status::UNSTABLE: @@ -169,7 +153,7 @@ template void lp_primal_core_solver::advance_on_en lp_assert((this->m_settings.simplex_strategy() == simplex_strategy_enum::tableau_rows) || m_non_basis_list.back() == static_cast(entering)); - lp_assert(this->using_infeas_costs() || !is_neg(t)); + lp_assert(!is_neg(t)); lp_assert(entering != leaving || !is_zero(t)); // otherwise nothing changes if (entering == leaving) { advance_on_entering_equal_leaving_tableau(entering, t); @@ -191,11 +175,6 @@ template void lp_primal_core_solver::advance_on_en return; if (this->m_settings.simplex_strategy() != simplex_strategy_enum::tableau_rows) { - if (need_to_switch_costs()) { - this->init_reduced_costs_tableau(); - } - - lp_assert(!need_to_switch_costs()); std::list::iterator it = m_non_basis_list.end(); it--; * it = static_cast(leaving); @@ -208,9 +187,7 @@ void lp_primal_core_solver::advance_on_entering_equal_leaving_tableau(int if (this->m_look_for_feasible_solution_only && this->current_x_is_feasible()) return; - if (need_to_switch_costs()) { - init_reduced_costs_tableau(); - } + this->iters_with_no_cost_growing() = 0; } template int lp_primal_core_solver::find_leaving_and_t_tableau(unsigned entering, X & t) { @@ -299,54 +276,19 @@ update_basis_and_x_tableau(int entering, int leaving, X const & tt) { this->change_basis(entering, leaving); return true; } + template void lp_primal_core_solver:: update_x_tableau(unsigned entering, 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]) { - unsigned i = c.var(); - this->add_delta_to_x_and_track_feasibility(this->m_basis[i], - 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 i = c.var(); - unsigned j = this->m_basis[i]; - 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); - } + for (const auto & c : this->m_A.m_columns[entering]) { + unsigned i = c.var(); + this->add_delta_to_x_and_track_feasibility(this->m_basis[i], - delta * this->m_A.get_val(c)); } } -template void lp_primal_core_solver:: -update_inf_cost_for_column_tableau(unsigned j) { - lp_assert(this->m_settings.simplex_strategy() != simplex_strategy_enum::tableau_rows); - - lp_assert(this->using_infeas_costs()); - - T new_cost = get_infeasibility_cost_for_column(j); - T delta = this->m_costs[j] - new_cost; - if (is_zero(delta)) - return; - this->m_costs[j] = new_cost; - update_reduced_cost_for_basic_column_cost_change(delta, j); -} template void lp_primal_core_solver::init_reduced_costs_tableau() { - if (this->current_x_is_infeasible() && !this->using_infeas_costs()) { - init_infeasibility_costs(); - } else if (this->current_x_is_feasible() && this->using_infeas_costs()) { - if (this->m_look_for_feasible_solution_only) - return; - this->m_costs = m_costs_backup; - this->set_using_infeas_costs(false); - } + unsigned size = this->m_basis_heading.size(); for (unsigned j = 0; j < size; j++) { if (this->m_basis_heading[j] >= 0)