diff --git a/src/math/lp/lar_core_solver_def.h b/src/math/lp/lar_core_solver_def.h index 425812db9c5..4ff4197d2db 100644 --- a/src/math/lp/lar_core_solver_def.h +++ b/src/math/lp/lar_core_solver_def.h @@ -199,7 +199,7 @@ void lar_core_solver::calculate_pivot_row(unsigned i) { m_r_solver.m_breakpoint_indices_queue.resize(m_r_solver.m_n()); 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.m_using_infeas_costs = true; + 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 bce1e9ef333..a6383a4001b 100644 --- a/src/math/lp/lar_solver.cpp +++ b/src/math/lp/lar_solver.cpp @@ -448,7 +448,7 @@ void lar_solver::prepare_costs_for_r_solver(const lar_term & term) { if (move_non_basic_columns_to_bounds()) find_feasible_solution(); auto & rslv = m_mpq_lar_core_solver.m_r_solver; - rslv.m_using_infeas_costs = false; + 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()); @@ -574,16 +574,19 @@ lar_term lar_solver::get_term_to_maximize(unsigned j_or_term) const { lp_status lar_solver::maximize_term(unsigned j_or_term, impq &term_max) { TRACE("lar_solver", print_values(tout);); - 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(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 (was_feasible) { + if (m_mpq_lar_core_solver.m_r_solver.calc_current_x_is_feasible_include_non_basis()) { prev_value = term.apply(m_mpq_lar_core_solver.m_r_x); + } else { + m_mpq_lar_core_solver.m_r_solver.m_look_for_feasible_solution_only = false; + if (solve() != lp_status::OPTIMAL) + return lp_status::UNBOUNDED; } m_mpq_lar_core_solver.m_r_solver.m_look_for_feasible_solution_only = false; @@ -618,7 +621,7 @@ lp_status lar_solver::maximize_term(unsigned j_or_term, if (change) { term_max = term.apply(m_mpq_lar_core_solver.m_r_x); } - if (was_feasible && term_max < prev_value) { + if (term_max < prev_value) { term_max = prev_value; m_mpq_lar_core_solver.m_r_x = backup; } @@ -822,7 +825,7 @@ void lar_solver::update_x_and_inf_costs_for_columns_with_changed_bounds_tableau( update_x_and_inf_costs_for_column_with_changed_bounds(j); if (tableau_with_costs()) { - if (m_mpq_lar_core_solver.m_r_solver.m_using_infeas_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()); diff --git a/src/math/lp/lp_core_solver_base.h b/src/math/lp/lp_core_solver_base.h index 666a4c1a795..02330fd5b9c 100644 --- a/src/math/lp/lp_core_solver_base.h +++ b/src/math/lp/lp_core_solver_base.h @@ -52,9 +52,11 @@ class lp_core_solver_base { } bool current_x_is_infeasible() const { return m_inf_set.size() != 0; } u_set m_inf_set; +private: bool m_using_infeas_costs; - - +public: + bool using_infeas_costs() const { return m_using_infeas_costs; } + void set_using_infeas_costs(bool val) { m_using_infeas_costs = val; } vector m_columns_nz; // m_columns_nz[i] keeps an approximate value of non zeroes the i-th column vector m_rows_nz; // m_rows_nz[i] keeps an approximate value of non zeroes in the i-th row indexed_vector m_pivot_row_of_B_1; // the pivot row of the reverse of B diff --git a/src/math/lp/lp_primal_core_solver.h b/src/math/lp/lp_primal_core_solver.h index 85326cbaabe..a8f4fcb7ae1 100644 --- a/src/math/lp/lp_primal_core_solver.h +++ b/src/math/lp/lp_primal_core_solver.h @@ -403,7 +403,7 @@ class lp_primal_core_solver:public lp_core_solver_base { 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->m_using_infeas_costs; + return this->current_x_is_feasible() == this->using_infeas_costs(); } @@ -445,13 +445,13 @@ class lp_primal_core_solver:public lp_core_solver_base { // 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->m_using_infeas_costs) { + 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 { // m_using_infeas_costs == 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. @@ -818,7 +818,7 @@ class lp_primal_core_solver:public lp_core_solver_base { void print_bound_info_and_x(unsigned j, std::ostream & out); void init_infeasibility_after_update_x_if_inf(unsigned leaving) { - if (this->m_using_infeas_costs) { + if (this->using_infeas_costs()) { init_infeasibility_costs_for_changed_basis_only(); this->m_costs[leaving] = zero_of_type(); this->m_inf_set.erase(leaving); diff --git a/src/math/lp/lp_primal_core_solver_def.h b/src/math/lp/lp_primal_core_solver_def.h index 015ac3694aa..85c1c4fe8bd 100644 --- a/src/math/lp/lp_primal_core_solver_def.h +++ b/src/math/lp/lp_primal_core_solver_def.h @@ -116,7 +116,7 @@ template bool lp_primal_core_solver::column_is_benefitial_for_entering_basis(unsigned j) const { if (numeric_traits::precise()) return column_is_benefitial_for_entering_basis_precise(j); - if (this->m_using_infeas_costs && this->m_settings.use_breakpoints_in_feasibility_search) + if (this->using_infeas_costs() && this->m_settings.use_breakpoints_in_feasibility_search) return column_is_benefitial_for_entering_on_breakpoints(j); const T& dj = this->m_d[j]; switch (this->m_column_types[j]) { @@ -150,7 +150,7 @@ bool lp_primal_core_solver::column_is_benefitial_for_entering_basis(unsign template bool lp_primal_core_solver::column_is_benefitial_for_entering_basis_precise(unsigned j) const { lp_assert (numeric_traits::precise()); - if (this->m_using_infeas_costs && this->m_settings.use_breakpoints_in_feasibility_search) + if (this->using_infeas_costs() && this->m_settings.use_breakpoints_in_feasibility_search) return column_is_benefitial_for_entering_on_breakpoints(j); const T& dj = this->m_d[j]; TRACE("lar_solver", tout << "dj=" << dj << "\n";); @@ -223,7 +223,7 @@ int lp_primal_core_solver::choose_entering_column_presize(unsigned number_ return -1; unsigned entering = *entering_iter; m_sign_of_entering_delta = this->m_d[entering] > 0 ? 1 : -1; - if (this->m_using_infeas_costs && this->m_settings.use_breakpoints_in_feasibility_search) + if (this->using_infeas_costs() && this->m_settings.use_breakpoints_in_feasibility_search) m_sign_of_entering_delta = -m_sign_of_entering_delta; m_non_basis_list.erase(entering_iter); m_non_basis_list.push_back(entering); @@ -263,7 +263,7 @@ int lp_primal_core_solver::choose_entering_column(unsigned number_of_benef if (entering_iter != m_non_basis_list.end()) { unsigned entering = *entering_iter; m_sign_of_entering_delta = this->m_d[entering] > 0? 1 : -1; - if (this->m_using_infeas_costs && this->m_settings.use_breakpoints_in_feasibility_search) + if (this->using_infeas_costs() && this->m_settings.use_breakpoints_in_feasibility_search) m_sign_of_entering_delta = - m_sign_of_entering_delta; m_non_basis_list.erase(entering_iter); m_non_basis_list.push_back(entering); @@ -648,7 +648,7 @@ template void lp_primal_core_solver::init_run( init_inf_set(); if (this->current_x_is_feasible() && this->m_look_for_feasible_solution_only) return; - this->m_using_infeas_costs = false; + this->set_using_infeas_costs(false); if (this->m_settings.backup_costs) backup_and_normalize_costs(); m_epsilon_of_reduced_cost = numeric_traits::precise()? zero_of_type(): T(1)/T(10000000); @@ -685,7 +685,7 @@ void lp_primal_core_solver::advance_on_entering_equal_leaving(int entering return; } } - if (this->m_using_infeas_costs) { + if (this->using_infeas_costs()) { lp_assert(is_zero(this->m_costs[entering])); init_infeasibility_costs_for_changed_basis_only(); } @@ -700,7 +700,7 @@ void lp_primal_core_solver::advance_on_entering_equal_leaving(int entering template void lp_primal_core_solver::advance_on_entering_and_leaving(int entering, int leaving, X & t) { lp_assert(entering >= 0 && m_non_basis_list.back() == static_cast(entering)); - lp_assert(this->m_using_infeas_costs || t >= zero_of_type()); + lp_assert(this->using_infeas_costs() || t >= zero_of_type()); lp_assert(leaving >= 0 && entering >= 0); lp_assert(entering != leaving || !is_zero(t)); // otherwise nothing changes if (entering == leaving) { @@ -879,13 +879,13 @@ template unsigned lp_primal_core_solver::solve() return 0; } do { - if (this->print_statistics_with_iterations_and_nonzeroes_and_cost_and_check_that_the_time_is_over((this->m_using_infeas_costs? "inf" : "feas"), * this->m_settings.get_message_ostream())) { + if (this->print_statistics_with_iterations_and_nonzeroes_and_cost_and_check_that_the_time_is_over((this->using_infeas_costs()? "inf" : "feas"), * this->m_settings.get_message_ostream())) { return this->total_iterations(); } one_iteration(); TRACE("lar_solver", tout << "one iteration: " << this->get_status() << "\n";); - lp_assert(!this->m_using_infeas_costs || this->costs_on_nbasis_are_zeros()); + lp_assert(!this->using_infeas_costs() || this->costs_on_nbasis_are_zeros()); switch (this->get_status()) { case lp_status::OPTIMAL: // double check that we are at optimum case lp_status::INFEASIBLE: @@ -1109,7 +1109,7 @@ 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->m_using_infeas_costs = true; + this->set_using_infeas_costs(true); } @@ -1119,7 +1119,7 @@ void lp_primal_core_solver::init_infeasibility_costs() { lp_assert(this->m_column_types.size() >= this->m_n()); for (unsigned j = this->m_n(); j--;) init_infeasibility_cost_for_column(j); - this->m_using_infeas_costs = true; + this->set_using_infeas_costs(true); } template T @@ -1297,13 +1297,13 @@ template void lp_primal_core_solver::print_breakp template void lp_primal_core_solver::init_reduced_costs() { lp_assert(!this->use_tableau()); - if (this->current_x_is_infeasible() && !this->m_using_infeas_costs) { + if (this->current_x_is_infeasible() && !this->using_infeas_costs()) { init_infeasibility_costs(); - } else if (this->current_x_is_feasible() && this->m_using_infeas_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->m_using_infeas_costs = false; + this->set_using_infeas_costs(false); } this->init_reduced_costs_for_one_iteration(); 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 59d23a74583..054c94cc04a 100644 --- a/src/math/lp/lp_primal_core_solver_tableau_def.h +++ b/src/math/lp/lp_primal_core_solver_tableau_def.h @@ -86,7 +86,7 @@ template int lp_primal_core_solver::choose_enteri return -1; unsigned entering = *entering_iter; m_sign_of_entering_delta = this->m_d[entering] > 0 ? 1 : -1; - if (this->m_using_infeas_costs && this->m_settings.use_breakpoints_in_feasibility_search) + if (this->using_infeas_costs() && this->m_settings.use_breakpoints_in_feasibility_search) m_sign_of_entering_delta = -m_sign_of_entering_delta; m_non_basis_list.erase(entering_iter); m_non_basis_list.push_back(entering); @@ -110,7 +110,7 @@ unsigned lp_primal_core_solver::solve_with_tableau() { return 0; } do { - if (this->print_statistics_with_iterations_and_nonzeroes_and_cost_and_check_that_the_time_is_over((this->m_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((this->using_infeas_costs()? "inf t" : "feas t"), * this->m_settings.get_message_ostream())) { return this->total_iterations(); } if (this->m_settings.use_tableau_rows()) { @@ -216,7 +216,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->m_using_infeas_costs || !is_neg(t)); + lp_assert(this->using_infeas_costs() || !is_neg(t)); lp_assert(entering != leaving || !is_zero(t)); // otherwise nothing changes if (entering == leaving) { advance_on_entering_equal_leaving_tableau(entering, t); @@ -359,12 +359,12 @@ update_basis_and_x_tableau(int entering, int leaving, X const & tt) { template void lp_primal_core_solver:: update_x_tableau(unsigned entering, const X& delta) { this->add_delta_to_x(entering, delta); - if (!this->m_using_infeas_costs) { + 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 { // m_using_infeas_costs == 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. @@ -386,7 +386,7 @@ 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->m_using_infeas_costs); + lp_assert(this->using_infeas_costs()); T new_cost = get_infeasibility_cost_for_column(j); T delta = this->m_costs[j] - new_cost; @@ -397,13 +397,13 @@ update_inf_cost_for_column_tableau(unsigned j) { } template void lp_primal_core_solver::init_reduced_costs_tableau() { - if (this->current_x_is_infeasible() && !this->m_using_infeas_costs) { + if (this->current_x_is_infeasible() && !this->using_infeas_costs()) { init_infeasibility_costs(); - } else if (this->current_x_is_feasible() && this->m_using_infeas_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->m_using_infeas_costs = false; + this->set_using_infeas_costs(false); } unsigned size = this->m_basis_heading.size(); for (unsigned j = 0; j < size; j++) {