Skip to content

Commit

Permalink
Revert "rm dealing with doubles"
Browse files Browse the repository at this point in the history 
This reverts commit 547254a.
  • Loading branch information
@levnach
levnach committed Mar 8, 2023
1 parent a418918 commit d00fcc8
Show file tree
Hide file tree
Showing 6 changed files with 245 additions and 21 deletions.
3 changes: 3 additions & 0 deletions src/math/lp/lp_core_solver_base.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,6 +83,9 @@ template std::string lp::lp_core_solver_base<lp::mpq, lp::mpq>::column_name(unsi
template void lp::lp_core_solver_base<lp::mpq, lp::mpq>::pretty_print(std::ostream & out);
template std::string lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::column_name(unsigned int) const;
template void lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::pretty_print(std::ostream & out);
template int lp::lp_core_solver_base<double, double>::pivots_in_column_and_row_are_different(int, int) const;
template int lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::pivots_in_column_and_row_are_different(int, int) const;
template int lp::lp_core_solver_base<lp::mpq, lp::mpq>::pivots_in_column_and_row_are_different(int, int) const;
template bool lp::lp_core_solver_base<double, double>::calc_current_x_is_feasible_include_non_basis(void)const;
template bool lp::lp_core_solver_base<lp::mpq, lp::mpq>::calc_current_x_is_feasible_include_non_basis(void)const;
template bool lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::calc_current_x_is_feasible_include_non_basis() const;
Expand Down
1 change: 1 addition & 0 deletions src/math/lp/lp_core_solver_base.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -416,6 +416,7 @@ 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&);
Expand Down
17 changes: 17 additions & 0 deletions src/math/lp/lp_core_solver_base_def.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -549,6 +549,23 @@ get_non_basic_column_value_position(unsigned j) const {
return at_lower_bound;
}

template <typename T, typename X> int lp_core_solver_base<T, X>::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<T>::one() + abs(row_p)));
if ( !this->m_settings.abs_val_is_smaller_than_harris_tolerance(diff_normalized / T(10)))
return 1;
return 0;
}
template <typename T, typename X> void lp_core_solver_base<T, X>::transpose_rows_tableau(unsigned i, unsigned j) {
transpose_basis(i, j);
m_A.transpose_rows(i, j);
Expand Down
137 changes: 118 additions & 19 deletions src/math/lp/lp_primal_core_solver.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,16 +43,20 @@ class lp_primal_core_solver:public lp_core_solver_base<T, X> {
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<T> m_beta; // see Swietanowski working vector beta for column norms
T m_epsilon_of_reduced_cost;
vector<T> 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<unsigned> m_leaving_candidates;
// T m_converted_harris_eps = convert_struct<T, double>::convert(this->m_settings.harris_feasibility_tolerance);
std::list<unsigned> m_non_basis_list;
void sort_non_basis();
void sort_non_basis_rational();
Expand Down Expand Up @@ -273,9 +277,14 @@ class lp_primal_core_solver:public lp_core_solver_base<T, X> {
return convert_struct<X, unsigned>::convert(std::numeric_limits<unsigned>::max());
}


bool get_harris_theta(X & theta);

void restore_harris_eps() { m_converted_harris_eps = convert_struct<T, double>::convert(this->m_settings.harris_feasibility_tolerance); }
void zero_harris_eps() { m_converted_harris_eps = zero_of_type<T>(); }
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);

Expand Down Expand Up @@ -309,7 +318,10 @@ class lp_primal_core_solver:public lp_core_solver_base<T, X> {
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<X, int>::convert(1) + abs(bound)/10) * m_converted_harris_eps/3;
}

void get_bound_on_variable_and_update_leaving_precisely(unsigned j, vector<unsigned> & leavings, T m, X & t, T & abs_of_d_of_leaving);

vector<T> m_lower_bounds_dummy; // needed for the base class only
Expand Down Expand Up @@ -502,7 +514,10 @@ class lp_primal_core_solver:public lp_core_solver_base<T, X> {
// }

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<X>()) theta = zero_of_type<X>();
}

bool limit_inf_on_bound_m_neg(const T & m, const X & x, const X & bound, X & theta, bool & unlimited) {
Expand All @@ -516,7 +531,16 @@ class lp_primal_core_solver:public lp_core_solver_base<T, X> {
theta = zero_of_type<X>();
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<X>();
unlimited = false;
}
}
return true;
}

Expand All @@ -531,7 +555,16 @@ class lp_primal_core_solver:public lp_core_solver_base<T, X> {
theta = zero_of_type<X>();
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<X>();
unlimited = false;
}
}
return true;
}

Expand All @@ -543,32 +576,82 @@ class lp_primal_core_solver:public lp_core_solver_base<T, X> {
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<X>();
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<X>();
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<X>()) {
theta = zero_of_type<X>();
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<X>()) {
theta = zero_of_type<X>();
}
}

// j is a basic column or the entering, in any case x[j] has to stay feasible.
Expand Down Expand Up @@ -639,12 +722,14 @@ class lp_primal_core_solver:public lp_core_solver_base<T, X> {

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
Expand All @@ -658,6 +743,13 @@ class lp_primal_core_solver:public lp_core_solver_base<T, X> {

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<T>();
this->remove_column_from_inf_set(leaving);
}
}

void init_inf_set() {
this->clear_inf_set();
Expand Down Expand Up @@ -793,10 +885,15 @@ class lp_primal_core_solver:public lp_core_solver_base<T, X> {
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<T>::precise())) {
m_converted_harris_eps = convert_struct<T, double>::convert(this->m_settings.harris_feasibility_tolerance);
} else {
m_converted_harris_eps = zero_of_type<T>();
}
this->set_status(lp_status::UNKNOWN);
}

Expand All @@ -822,7 +919,9 @@ class lp_primal_core_solver:public lp_core_solver_base<T, X> {
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<T, double>::convert(this->m_settings.harris_feasibility_tolerance)) {
lp_assert(initial_x_is_correct());
m_lower_bounds_dummy.resize(A.column_count(), zero_of_type<T>());
m_enter_price_eps = numeric_traits<T>::precise() ? numeric_traits<T>::zero() : T(1e-5);
Expand Down

0 comments on commit d00fcc8

Please sign in to comment.