diff --git a/src/math/lp/lp_core_solver_base.cpp b/src/math/lp/lp_core_solver_base.cpp index b76976d371c..e87dc5dd2f6 100644 --- a/src/math/lp/lp_core_solver_base.cpp +++ b/src/math/lp/lp_core_solver_base.cpp @@ -83,9 +83,6 @@ template std::string lp::lp_core_solver_base::column_name(unsi template void lp::lp_core_solver_base::pretty_print(std::ostream & out); template std::string lp::lp_core_solver_base >::column_name(unsigned int) const; template void lp::lp_core_solver_base >::pretty_print(std::ostream & out); -template int lp::lp_core_solver_base::pivots_in_column_and_row_are_different(int, int) const; -template int lp::lp_core_solver_base >::pivots_in_column_and_row_are_different(int, int) const; -template int lp::lp_core_solver_base::pivots_in_column_and_row_are_different(int, int) const; template bool lp::lp_core_solver_base::calc_current_x_is_feasible_include_non_basis(void)const; template bool lp::lp_core_solver_base::calc_current_x_is_feasible_include_non_basis(void)const; template bool lp::lp_core_solver_base >::calc_current_x_is_feasible_include_non_basis() const; diff --git a/src/math/lp/lp_core_solver_base.h b/src/math/lp/lp_core_solver_base.h index 5f723c1185f..fff50090591 100644 --- a/src/math/lp/lp_core_solver_base.h +++ b/src/math/lp/lp_core_solver_base.h @@ -416,7 +416,6 @@ class lp_core_solver_base { non_basic_column_value_position get_non_basic_column_value_position(unsigned j) const; - int pivots_in_column_and_row_are_different(int entering, int leaving) const; void pivot_fixed_vars_from_basis(); bool remove_from_basis(unsigned j); bool remove_from_basis(unsigned j, const impq&); diff --git a/src/math/lp/lp_core_solver_base_def.h b/src/math/lp/lp_core_solver_base_def.h index 8b87728e090..c5910f4274c 100644 --- a/src/math/lp/lp_core_solver_base_def.h +++ b/src/math/lp/lp_core_solver_base_def.h @@ -549,23 +549,6 @@ get_non_basic_column_value_position(unsigned j) const { return at_lower_bound; } -template int lp_core_solver_base::pivots_in_column_and_row_are_different(int entering, int leaving) const { - const T & column_p = this->m_ed[this->m_basis_heading[leaving]]; - const T & row_p = this->m_pivot_row[entering]; - if (is_zero(column_p) || is_zero(row_p)) return true; // pivots cannot be zero - // the pivots have to have the same sign - if (column_p < 0) { - if (row_p > 0) - return 2; - } else { // column_p > 0 - if (row_p < 0) - return 2; - } - T diff_normalized = abs((column_p - row_p) / (numeric_traits::one() + abs(row_p))); - if ( !this->m_settings.abs_val_is_smaller_than_harris_tolerance(diff_normalized / T(10))) - return 1; - return 0; -} template void lp_core_solver_base::transpose_rows_tableau(unsigned i, unsigned j) { transpose_basis(i, j); m_A.transpose_rows(i, j); diff --git a/src/math/lp/lp_primal_core_solver.h b/src/math/lp/lp_primal_core_solver.h index 2c7360712d6..568eefa9829 100644 --- a/src/math/lp/lp_primal_core_solver.h +++ b/src/math/lp/lp_primal_core_solver.h @@ -43,20 +43,16 @@ class lp_primal_core_solver:public lp_core_solver_base { public: // m_sign_of_entering is set to 1 if the entering variable needs // to grow and is set to -1 otherwise - unsigned m_column_norm_update_counter; T m_enter_price_eps; int m_sign_of_entering_delta; - indexed_vector m_beta; // see Swietanowski working vector beta for column norms T m_epsilon_of_reduced_cost; vector m_costs_backup; - T m_converted_harris_eps; 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; - // T m_converted_harris_eps = convert_struct::convert(this->m_settings.harris_feasibility_tolerance); std::list m_non_basis_list; void sort_non_basis(); void sort_non_basis_rational(); @@ -277,14 +273,9 @@ class lp_primal_core_solver:public lp_core_solver_base { return convert_struct::convert(std::numeric_limits::max()); } - bool get_harris_theta(X & theta); - - void restore_harris_eps() { m_converted_harris_eps = convert_struct::convert(this->m_settings.harris_feasibility_tolerance); } - void zero_harris_eps() { m_converted_harris_eps = zero_of_type(); } - int find_leaving_on_harris_theta(X const & harris_theta, X & t); + 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(unsigned entering, X & t); int find_leaving_and_t_precise(unsigned entering, X & t); int find_leaving_and_t_tableau(unsigned entering, X & t); @@ -318,10 +309,7 @@ class lp_primal_core_solver:public lp_core_solver_base { 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); }; - - X harris_eps_for_bound(const X & bound) const { return ( convert_struct::convert(1) + abs(bound)/10) * m_converted_harris_eps/3; - } - + void get_bound_on_variable_and_update_leaving_precisely(unsigned j, vector & leavings, T m, X & t, T & abs_of_d_of_leaving); vector m_lower_bounds_dummy; // needed for the base class only @@ -514,10 +502,7 @@ class lp_primal_core_solver:public lp_core_solver_base { // } 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); - const X& eps = harris_eps_for_bound(this->m_lower_bounds[j]); - limit_theta((this->m_lower_bounds[j] - this->m_x[j] - eps) / m, theta, unlimited); - if (theta < zero_of_type()) theta = zero_of_type(); + lp_assert(false); } bool limit_inf_on_bound_m_neg(const T & m, const X & x, const X & bound, X & theta, bool & unlimited) { @@ -531,16 +516,7 @@ class lp_primal_core_solver:public lp_core_solver_base { theta = zero_of_type(); unlimited = false; } - } else { - const X& eps = harris_eps_for_bound(bound); - if (this->below_bound(x, bound)) return false; - if (this->above_bound(x, bound)) { - limit_theta((bound - x - eps) / m, theta, unlimited); - } else { - theta = zero_of_type(); - unlimited = false; - } - } + } return true; } @@ -555,16 +531,7 @@ class lp_primal_core_solver:public lp_core_solver_base { theta = zero_of_type(); unlimited = false; } - } else { - const X& eps = harris_eps_for_bound(bound); - if (this->above_bound(x, bound)) return false; - if (this->below_bound(x, bound)) { - limit_theta((bound - x + eps) / m, theta, unlimited); - } else { - theta = zero_of_type(); - unlimited = false; - } - } + } return true; } @@ -576,82 +543,32 @@ class lp_primal_core_solver:public lp_core_solver_base { limit_theta((bound - x) / m, theta, unlimited); } } - else { - // x gets larger - lp_assert(m > 0); - const X& eps = harris_eps_for_bound(bound); - if (this->below_bound(x, bound)) { - limit_theta((bound - x + eps) / 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); - const X& eps = harris_eps_for_bound(bound); - if (this->above_bound(x, bound)) { - limit_theta((bound - x - eps) / m, theta, unlimited); - } + lp_assert(false); + } void limit_theta_on_basis_column_for_inf_case_m_pos_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 & lbound = this->m_lower_bounds[j]; - - if (this->below_bound(x, lbound)) { - const X& eps = harris_eps_for_bound(this->m_upper_bounds[j]); - limit_theta((lbound - x + eps) / m, theta, unlimited); - } else { - const X & ubound = this->m_upper_bounds[j]; - if (this->below_bound(x, ubound)){ - const X& eps = harris_eps_for_bound(ubound); - limit_theta((ubound - x + eps) / m, theta, unlimited); - } else if (!this->above_bound(x, ubound)) { - theta = zero_of_type(); - unlimited = false; - } - } + lp_assert(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)) { - const X& eps = harris_eps_for_bound(ubound); - limit_theta((ubound - x - eps) / m, theta, unlimited); - } else { - const X & lbound = this->m_lower_bounds[j]; - if (this->above_bound(x, lbound)){ - const X& eps = harris_eps_for_bound(lbound); - limit_theta((lbound - x - eps) / m, theta, unlimited); - } else if (!this->below_bound(x, lbound)) { - theta = zero_of_type(); - unlimited = false; - } - } + lp_assert(false); + } void limit_theta_on_basis_column_for_feas_case_m_pos(unsigned j, const T & m, X & theta, bool & unlimited) { - lp_assert(m > 0); - const T& eps = harris_eps_for_bound(this->m_upper_bounds[j]); - if (this->below_bound(this->m_x[j], this->m_upper_bounds[j])) { - limit_theta((this->m_upper_bounds[j] - this->m_x[j] + eps) / m, theta, unlimited); - if (theta < zero_of_type()) { - theta = zero_of_type(); - unlimited = false; - } - } + lp_assert(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); - const X& eps = harris_eps_for_bound(this->m_upper_bounds[j]); - limit_theta( (this->m_upper_bounds[j] - this->m_x[j] + eps) / m, theta, unlimited); - if (theta < zero_of_type()) { - theta = zero_of_type(); - } + lp_assert(false); + } // j is a basic column or the entering, in any case x[j] has to stay feasible. @@ -722,14 +639,12 @@ class lp_primal_core_solver:public lp_core_solver_base { 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); bool done(); 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); // j is the basic column, x is the value at x[j] // d is the coefficient before m_entering in the row with j as the basis column @@ -743,13 +658,6 @@ 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->using_infeas_costs()) { - init_infeasibility_costs_for_changed_basis_only(); - this->m_costs[leaving] = zero_of_type(); - this->remove_column_from_inf_set(leaving); - } - } void init_inf_set() { this->clear_inf_set(); @@ -885,15 +793,10 @@ class lp_primal_core_solver:public lp_core_solver_base { column_type_array, lower_bound_values, upper_bound_values), - m_beta(A.row_count()), m_epsilon_of_reduced_cost(T(1)/T(10000000)), m_bland_mode_threshold(1000) { - if (!(numeric_traits::precise())) { - m_converted_harris_eps = convert_struct::convert(this->m_settings.harris_feasibility_tolerance); - } else { - m_converted_harris_eps = zero_of_type(); - } + this->set_status(lp_status::UNKNOWN); } @@ -919,9 +822,7 @@ class lp_primal_core_solver:public lp_core_solver_base { column_names, column_type_array, m_lower_bounds_dummy, - upper_bound_values), - m_beta(A.row_count()), - m_converted_harris_eps(convert_struct::convert(this->m_settings.harris_feasibility_tolerance)) { + upper_bound_values) { lp_assert(initial_x_is_correct()); m_lower_bounds_dummy.resize(A.column_count(), zero_of_type()); m_enter_price_eps = numeric_traits::precise() ? numeric_traits::zero() : T(1e-5); diff --git a/src/math/lp/lp_primal_core_solver_def.h b/src/math/lp/lp_primal_core_solver_def.h index 7a027f8d38c..22bd4db52e1 100644 --- a/src/math/lp/lp_primal_core_solver_def.h +++ b/src/math/lp/lp_primal_core_solver_def.h @@ -228,54 +228,8 @@ int lp_primal_core_solver::choose_entering_column(unsigned number_of_benef } -template bool lp_primal_core_solver::get_harris_theta(X & theta) { - lp_assert(this->m_ed.is_OK()); - bool unlimited = true; - for (unsigned i : this->m_ed.m_index) { - if (this->m_settings.abs_val_is_smaller_than_pivot_tolerance(this->m_ed[i])) continue; - limit_theta_on_basis_column(this->m_basis[i], - this->m_ed[i] * m_sign_of_entering_delta, theta, unlimited); - if (!unlimited && is_zero(theta)) break; - } - return unlimited; -} -template int lp_primal_core_solver:: -find_leaving_on_harris_theta(X const & harris_theta, X & t) { - int leaving = -1; - T pivot_abs_max = zero_of_type(); - // we know already that there is no bound flip on entering - // we also know that harris_theta is limited, so we will find a leaving - zero_harris_eps(); - unsigned steps = this->m_ed.m_index.size(); - unsigned k = this->m_settings.random_next() % steps; - unsigned initial_k = k; - do { - unsigned i = this->m_ed.m_index[k]; - const T & ed = this->m_ed[i]; - if (this->m_settings.abs_val_is_smaller_than_pivot_tolerance(ed)) { - if (++k == steps) - k = 0; - continue; - } - X ratio; - unsigned j = this->m_basis[i]; - bool unlimited = true; - limit_theta_on_basis_column(j, - ed * m_sign_of_entering_delta, ratio, unlimited); - if ((!unlimited) && ratio <= harris_theta) { - if (leaving == -1 || abs(ed) > pivot_abs_max) { - t = ratio; - leaving = j; - pivot_abs_max = abs(ed); - } - } - if (++k == steps) k = 0; - } while (k != initial_k); - if (!this->precise()) - restore_harris_eps(); - return leaving; -} - template bool lp_primal_core_solver::try_jump_to_another_bound_on_entering(unsigned entering, const X & theta, @@ -398,17 +352,6 @@ template int lp_primal_core_solver::find_leaving_ } -template int lp_primal_core_solver::find_leaving_and_t(unsigned entering, X & t) { - X theta = zero_of_type(); - bool unlimited = get_harris_theta(theta); - lp_assert(unlimited || theta >= zero_of_type()); - if (try_jump_to_another_bound_on_entering(entering, theta, t, unlimited)) return entering; - if (unlimited) - return -1; - return find_leaving_on_harris_theta(theta, t); -} - - // m is the multiplier. updating t in a way that holds the following // x[j] + t * m >= m_lower_bounds[j] ( if m < 0 ) @@ -687,47 +630,10 @@ template void lp_primal_core_solver::init_column_ } } -// debug only -template T lp_primal_core_solver::calculate_column_norm_exactly(unsigned j) { - lp_assert(false); -} -template void lp_primal_core_solver::update_or_init_column_norms(unsigned entering, unsigned leaving) { - lp_assert(numeric_traits::precise() == false); - lp_assert(m_column_norm_update_counter <= this->m_settings.column_norms_update_frequency); - if (m_column_norm_update_counter == this->m_settings.column_norms_update_frequency) { - m_column_norm_update_counter = 0; - init_column_norms(); - } else { - m_column_norm_update_counter++; - update_column_norms(entering, leaving); - } -} -// following Swietanowski - A new steepest ... -template void lp_primal_core_solver::update_column_norms(unsigned entering, unsigned leaving) { - lp_assert(numeric_traits::precise() == false); - T pivot = this->m_pivot_row[entering]; - T g_ent = calculate_norm_of_entering_exactly() / pivot / pivot; - if (!numeric_traits::precise()) { - if (g_ent < T(0.000001)) - g_ent = T(0.000001); - } - this->m_column_norms[leaving] = g_ent; - for (unsigned j : this->m_pivot_row.m_index) { - if (j == leaving) - continue; - const T & t = this->m_pivot_row[j]; - T s = this->m_A.dot_product_with_column(m_beta.m_data, j); - T k = -2 / pivot; - T tp = t/pivot; - if (this->m_column_types[j] != column_type::fixed) { // a fixed columns do not enter the basis, we don't use the norm of a fixed column - this->m_column_norms[j] = std::max(this->m_column_norms[j] + t * (t * g_ent + k * s), // see Istvan Maros, page 196 - 1 + tp * tp); - } - } -} +// following Swietanowski - A new steepest ... template T lp_primal_core_solver::calculate_norm_of_entering_exactly() { T r = numeric_traits::one(); @@ -771,13 +677,6 @@ template bool lp_primal_core_solver::done() { return false; } -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() { 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 b63c54bfd10..a63a6be17ad 100644 --- a/src/math/lp/lp_primal_core_solver_tableau_def.h +++ b/src/math/lp/lp_primal_core_solver_tableau_def.h @@ -295,10 +295,7 @@ template void lp_primal_core_solver::init_run_tab 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); - if (!numeric_traits::precise()) { - this->m_column_norm_update_counter = 0; - init_column_norms(); - } + if (this->m_settings.simplex_strategy() == simplex_strategy_enum::tableau_rows) init_tableau_rows(); lp_assert(this->reduced_costs_are_correct_tableau());