From d00fcc87c9a822d3a995f527e5bca69c400f60b0 Mon Sep 17 00:00:00 2001 From: Lev Nachmanson Date: Mon, 6 Mar 2023 10:56:44 -0800 Subject: [PATCH] Revert "rm dealing with doubles" This reverts commit 547254abe786c80231ca78bcd245e6ddb5a15c47. --- src/math/lp/lp_core_solver_base.cpp | 3 + src/math/lp/lp_core_solver_base.h | 1 + src/math/lp/lp_core_solver_base_def.h | 17 +++ src/math/lp/lp_primal_core_solver.h | 137 +++++++++++++++--- src/math/lp/lp_primal_core_solver_def.h | 103 ++++++++++++- .../lp/lp_primal_core_solver_tableau_def.h | 5 +- 6 files changed, 245 insertions(+), 21 deletions(-) diff --git a/src/math/lp/lp_core_solver_base.cpp b/src/math/lp/lp_core_solver_base.cpp index e87dc5dd2..b76976d37 100644 --- a/src/math/lp/lp_core_solver_base.cpp +++ b/src/math/lp/lp_core_solver_base.cpp @@ -83,6 +83,9 @@ 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 fff500905..5f723c118 100644 --- a/src/math/lp/lp_core_solver_base.h +++ b/src/math/lp/lp_core_solver_base.h @@ -416,6 +416,7 @@ public: 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 c5910f427..8b87728e0 100644 --- a/src/math/lp/lp_core_solver_base_def.h +++ b/src/math/lp/lp_core_solver_base_def.h @@ -549,6 +549,23 @@ 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 568eefa98..2c7360712 100644 --- a/src/math/lp/lp_primal_core_solver.h +++ b/src/math/lp/lp_primal_core_solver.h @@ -43,16 +43,20 @@ 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(); @@ -273,9 +277,14 @@ public: 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); @@ -309,7 +318,10 @@ public: 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 @@ -502,7 +514,10 @@ public: // } void limit_theta_on_basis_column_for_feas_case_m_neg_no_check(unsigned j, const T & m, X & theta, bool & unlimited) { - lp_assert(false); + 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(); } bool limit_inf_on_bound_m_neg(const T & m, const X & x, const X & bound, X & theta, bool & unlimited) { @@ -516,7 +531,16 @@ public: 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; } @@ -531,7 +555,16 @@ public: 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; } @@ -543,32 +576,82 @@ public: 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(false); - + 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); + } } void limit_theta_on_basis_column_for_inf_case_m_pos_boxed(unsigned j, const T & m, X & theta, bool & unlimited) { - lp_assert(false); - + // 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; + } + } } void limit_theta_on_basis_column_for_inf_case_m_neg_boxed(unsigned j, const T & m, X & theta, bool & unlimited) { - lp_assert(false); - + // 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; + } + } } void limit_theta_on_basis_column_for_feas_case_m_pos(unsigned j, const T & m, X & theta, bool & unlimited) { - lp_assert(false); - + 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; + } + } } void limit_theta_on_basis_column_for_feas_case_m_pos_no_check(unsigned j, const T & m, X & theta, bool & unlimited ) { - lp_assert(false); - + 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(); + } } // j is a basic column or the entering, in any case x[j] has to stay feasible. @@ -639,12 +722,14 @@ public: 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 @@ -658,6 +743,13 @@ public: 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(); @@ -793,10 +885,15 @@ public: 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); } @@ -822,7 +919,9 @@ public: column_names, column_type_array, m_lower_bounds_dummy, - upper_bound_values) { + upper_bound_values), + m_beta(A.row_count()), + m_converted_harris_eps(convert_struct::convert(this->m_settings.harris_feasibility_tolerance)) { 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 22bd4db52..7a027f8d3 100644 --- a/src/math/lp/lp_primal_core_solver_def.h +++ b/src/math/lp/lp_primal_core_solver_def.h @@ -228,8 +228,54 @@ 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, @@ -352,6 +398,17 @@ 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 ) @@ -630,10 +687,47 @@ 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); + } + } +} template T lp_primal_core_solver::calculate_norm_of_entering_exactly() { T r = numeric_traits::one(); @@ -677,6 +771,13 @@ 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 a63a6be17..b63c54bfd 100644 --- a/src/math/lp/lp_primal_core_solver_tableau_def.h +++ b/src/math/lp/lp_primal_core_solver_tableau_def.h @@ -295,7 +295,10 @@ 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());