Skip to content

Commit

Permalink
rp precise
Browse files Browse the repository at this point in the history
  • Loading branch information
@levnach
levnach committed Mar 8, 2023
1 parent 569f5be commit 9ec8263
Show file tree
Hide file tree
Showing 15 changed files with 48 additions and 350 deletions.
4 changes: 1 addition & 3 deletions src/math/lp/core_solver_pretty_printer_def.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -88,7 +88,6 @@ template <typename T, typename X> T core_solver_pretty_printer<T, X>::current_co
}

template <typename T, typename X> void core_solver_pretty_printer<T, X>::init_m_A_and_signs() {
if (numeric_traits<T>::precise() ) {
for (unsigned column = 0; column < ncols(); column++) {
vector<T> t(nrows(), zero_of_type<T>());
for (const auto & c : m_core_solver.m_A.m_columns[column]){
Expand All @@ -115,8 +114,7 @@ template <typename T, typename X> void core_solver_pretty_printer<T, X>::init_m_
name);
m_rs[row] += t[row] * m_core_solver.m_x[column];
}
}
}
}
}

template <typename T, typename X> void core_solver_pretty_printer<T, X>::init_column_widths() {
Expand Down
45 changes: 12 additions & 33 deletions src/math/lp/lar_solver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -281,9 +281,9 @@ namespace lp {
unsigned m = A_r().row_count();
clean_popped_elements(m, m_rows_with_changed_bounds);
clean_inf_set_of_r_solver_after_pop();
lp_assert(m_settings.simplex_strategy() == simplex_strategy_enum::undecided ||
( m_mpq_lar_core_solver.m_r_solver.reduced_costs_are_correct_tableau()));

lp_assert(
m_settings.simplex_strategy() == simplex_strategy_enum::undecided ||
m_mpq_lar_core_solver.m_r_solver.reduced_costs_are_correct_tableau());

m_constraints.pop(k);
m_term_count.pop(k);
Expand Down Expand Up @@ -627,10 +627,6 @@ namespace lp {
m_rows_with_changed_bounds.insert(rid);
}

void lar_solver::detect_rows_of_bound_change_column_for_nbasic_column(unsigned j) {
lp_assert(false);
}



void lar_solver::detect_rows_of_bound_change_column_for_nbasic_column_tableau(unsigned j) {
Expand Down Expand Up @@ -676,20 +672,16 @@ namespace lp {
}

void lar_solver::change_basic_columns_dependend_on_a_given_nb_column(unsigned j, const numeric_pair<mpq>& delta) {
{
for (const auto& c : A_r().m_columns[j]) {
unsigned bj = m_mpq_lar_core_solver.m_r_basis[c.var()];
if (tableau_with_costs()) {
m_basic_columns_with_changed_cost.insert(bj);
}
m_mpq_lar_core_solver.m_r_solver.add_delta_to_x_and_track_feasibility(bj, -A_r().get_val(c) * delta);
TRACE("change_x_del",
tout << "changed basis column " << bj << ", it is " <<
(m_mpq_lar_core_solver.m_r_solver.column_is_feasible(bj) ? "feas" : "inf") << std::endl;);

}
for (const auto& c : A_r().m_columns[j]) {
unsigned bj = m_mpq_lar_core_solver.m_r_basis[c.var()];
if (tableau_with_costs()) {
m_basic_columns_with_changed_cost.insert(bj);
}
m_mpq_lar_core_solver.m_r_solver.add_delta_to_x_and_track_feasibility(bj, -A_r().get_val(c) * delta);
TRACE("change_x_del",
tout << "changed basis column " << bj << ", it is " <<
(m_mpq_lar_core_solver.m_r_solver.column_is_feasible(bj) ? "feas" : "inf") << std::endl;);
}

}

void lar_solver::update_x_and_inf_costs_for_column_with_changed_bounds(unsigned j) {
Expand Down Expand Up @@ -819,19 +811,6 @@ namespace lp {
A.set(last_row, basis_j, mpq(1));
}

template <typename U, typename V>
void lar_solver::create_matrix_A(static_matrix<U, V>& matr) {
lp_assert(false); // not implemented
/*
unsigned m = number_or_nontrivial_left_sides();
unsigned n = m_vec_of_canonic_left_sides.size();
if (matr.row_count() == m && matr.column_count() == n)
return;
matr.init_empty_matrix(m, n);
copy_from_mpq_matrix(matr);
*/
}

template <typename U, typename V>
void lar_solver::copy_from_mpq_matrix(static_matrix<U, V>& matr) {
matr.m_rows.resize(A_r().row_count());
Expand Down
3 changes: 1 addition & 2 deletions src/math/lp/lar_solver.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -205,10 +205,9 @@ class lar_solver : public column_namer {
void set_lower_bound_witness(var_index j, constraint_index ci);
void substitute_terms_in_linear_expression( const vector<std::pair<mpq, var_index>>& left_side_with_terms,
vector<std::pair<mpq, var_index>> &left_side) const;
void detect_rows_of_bound_change_column_for_nbasic_column(unsigned j);

void detect_rows_of_bound_change_column_for_nbasic_column_tableau(unsigned j);
bool use_tableau_costs() const;
void detect_rows_of_column_with_bound_change(unsigned j);
bool tableau_with_costs() const;
bool costs_are_used() const;
void change_basic_columns_dependend_on_a_given_nb_column(unsigned j, const numeric_pair<mpq> & delta);
Expand Down
3 changes: 0 additions & 3 deletions src/math/lp/lp_core_solver_base.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,15 +43,13 @@ template bool lp::lp_core_solver_base<double, double>::print_statistics_with_ite
template bool lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::print_statistics_with_iterations_and_nonzeroes_and_cost_and_check_that_the_time_is_over(char const*, std::ostream &);
template void lp::lp_core_solver_base<double, double>::restore_x(unsigned int, double const&);
template void lp::lp_core_solver_base<double, double>::set_non_basic_x_to_correct_bounds();
template void lp::lp_core_solver_base<double, double>::solve_Ax_eq_b();
template void lp::lp_core_solver_base<double, double>::add_delta_to_entering(unsigned int, const double&);
template bool lp::lp_core_solver_base<lp::mpq, lp::mpq>::basis_heading_is_correct() const ;
template bool lp::lp_core_solver_base<lp::mpq, lp::mpq>::column_is_dual_feasible(unsigned int) const;
template void lp::lp_core_solver_base<lp::mpq, lp::mpq>::fill_reduced_costs_from_m_y_by_rows();
template bool lp::lp_core_solver_base<lp::mpq, lp::mpq>::print_statistics_with_iterations_and_nonzeroes_and_cost_and_check_that_the_time_is_over(char const*, std::ostream &);
template void lp::lp_core_solver_base<lp::mpq, lp::mpq>::restore_x(unsigned int, lp::mpq const&);
template void lp::lp_core_solver_base<lp::mpq, lp::mpq>::set_non_basic_x_to_correct_bounds();
template void lp::lp_core_solver_base<lp::mpq, lp::mpq>::solve_Ax_eq_b();
template void lp::lp_core_solver_base<lp::mpq, lp::mpq>::add_delta_to_entering(unsigned int, const lp::mpq&);
template void lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::init();
template void lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::init_basis_heading_and_non_basic_columns_vector();
Expand All @@ -61,7 +59,6 @@ template lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::lp_core_s
const vector<lp::numeric_pair<lp::mpq> >&,
const vector<lp::numeric_pair<lp::mpq> >&);
template bool lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::print_statistics_with_cost_and_check_that_the_time_is_over(lp::numeric_pair<lp::mpq>, std::ostream&);
template void lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::solve_Ax_eq_b();

template void lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::add_delta_to_entering(unsigned int, const lp::numeric_pair<lp::mpq>&);
template lp::lp_core_solver_base<lp::mpq, lp::mpq>::lp_core_solver_base(
Expand Down
9 changes: 2 additions & 7 deletions src/math/lp/lp_core_solver_base.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -238,11 +238,11 @@ class lp_core_solver_base {
}

bool below_bound(const X & x, const X & bound) const {
return precise()? x < bound : below_bound_numeric<X>(x, bound, m_settings.primal_feasibility_tolerance);
return x < bound ;
}

bool above_bound(const X & x, const X & bound) const {
return precise()? x > bound : above_bound_numeric<X>(x, bound, m_settings.primal_feasibility_tolerance);
return x > bound ;
}

bool x_below_low_bound(unsigned p) const {
Expand Down Expand Up @@ -323,9 +323,6 @@ class lp_core_solver_base {

void find_error_in_BxB(vector<X>& rs);

// recalculates the projection of x to B, such that Ax = b, whereab is the right side
void solve_Ax_eq_b();

bool snap_non_basic_x_to_bound() {
bool ret = false;
for (unsigned j : non_basis())
Expand Down Expand Up @@ -628,8 +625,6 @@ class lp_core_solver_base {
bool pivot_column_tableau(unsigned j, unsigned row_index);
bool divide_row_by_pivot(unsigned pivot_row, unsigned pivot_col);

bool precise() const { return numeric_traits<T>::precise(); }

simplex_strategy_enum simplex_strategy() const { return
m_settings.simplex_strategy();
}
Expand Down
18 changes: 2 additions & 16 deletions src/math/lp/lp_core_solver_base_def.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -233,18 +233,12 @@ column_is_dual_feasible(unsigned j) const {
}
template <typename T, typename X> bool lp_core_solver_base<T, X>::
d_is_not_negative(unsigned j) const {
if (numeric_traits<T>::precise()) {
return m_d[j] >= numeric_traits<T>::zero();
}
return m_d[j] > -T(0.00001);
return m_d[j] >= numeric_traits<T>::zero();
}

template <typename T, typename X> bool lp_core_solver_base<T, X>::
d_is_not_positive(unsigned j) const {
if (numeric_traits<T>::precise()) {
return m_d[j] <= numeric_traits<T>::zero();
}
return m_d[j] < T(0.00001);
return m_d[j] <= numeric_traits<T>::zero();
}


Expand Down Expand Up @@ -319,7 +313,6 @@ template <typename T, typename X> bool lp_core_solver_base<T, X>::inf_set_is_cor

template <typename T, typename X> bool lp_core_solver_base<T, X>::
divide_row_by_pivot(unsigned pivot_row, unsigned pivot_col) {
lp_assert(numeric_traits<T>::precise());
int pivot_index = -1;
auto & row = m_A.m_rows[pivot_row];
unsigned size = row.size();
Expand Down Expand Up @@ -517,13 +510,6 @@ find_error_in_BxB(vector<X>& rs){
}
}

// recalculates the projection of x to B, such that Ax = b
template <typename T, typename X> void lp_core_solver_base<T, X>::
solve_Ax_eq_b() {
lp_assert(false);
}


template <typename T, typename X> non_basic_column_value_position lp_core_solver_base<T, X>::
get_non_basic_column_value_position(unsigned j) const {
switch (m_column_types[j]) {
Expand Down
93 changes: 15 additions & 78 deletions src/math/lp/lp_primal_core_solver.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -511,67 +511,37 @@ class lp_primal_core_solver:public lp_core_solver_base<T, X> {
bool limit_inf_on_bound_m_neg(const T & m, const X & x, const X & bound, X & theta, bool & unlimited) {
// x gets smaller
lp_assert(m < 0);
if (numeric_traits<T>::precise()) {
if (this->below_bound(x, bound)) return false;
if (this->above_bound(x, bound)) {
limit_theta((bound - x) / m, theta, unlimited);
} else {
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;
}

bool limit_inf_on_bound_m_pos(const T & m, const X & x, const X & bound, X & theta, bool & unlimited) {
// x gets larger
lp_assert(m > 0);
if (numeric_traits<T>::precise()) {
if (this->above_bound(x, bound)) return false;
if (this->below_bound(x, bound)) {
limit_theta((bound - x) / m, theta, unlimited);
} else {
theta = zero_of_type<X>();
unlimited = false;
}
if (this->above_bound(x, bound)) return false;
if (this->below_bound(x, bound)) {
limit_theta((bound - x) / m, theta, unlimited);
} 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;
}
theta = zero_of_type<X>();
unlimited = false;
}

return true;
}

void limit_inf_on_lower_bound_m_pos(const T & m, const X & x, const X & bound, X & theta, bool & unlimited) {
if (numeric_traits<T>::precise()) {
// x gets larger
lp_assert(m > 0);
if (this->below_bound(x, bound)) {
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);
}
}
lp_assert(m > 0);
if (this->below_bound(x, bound)) {
limit_theta((bound - x) / m, theta, unlimited);
}

}

void limit_inf_on_upper_bound_m_neg(const T & m, const X & x, const X & bound, X & theta, bool & unlimited) {
Expand Down Expand Up @@ -877,46 +847,13 @@ class lp_primal_core_solver:public lp_core_solver_base<T, X> {
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>();
}

m_converted_harris_eps = zero_of_type<T>();

this->set_status(lp_status::UNKNOWN);
}

// constructor
lp_primal_core_solver(static_matrix<T, X> & A,
vector<X> & b, // the right side vector
vector<X> & x, // the number of elements in x needs to be at least as large as the number of columns in A
vector<unsigned> & basis,
vector<unsigned> & nbasis,
vector<int> & heading,
vector<T> & costs,
const vector<column_type> & column_type_array,
const vector<X> & upper_bound_values,
lp_settings & settings,
const column_namer& column_names):
lp_core_solver_base<T, X>(A, // b,
basis,
nbasis,
heading,
x,
costs,
settings,
column_names,
column_type_array,
m_lower_bounds_dummy,
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);
#ifdef Z3DEBUG
lp_assert(false);
#endif
}


bool initial_x_is_correct() {
std::set<unsigned> basis_set;
Expand Down

0 comments on commit 9ec8263

Please sign in to comment.