mirrors/z3
3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2025-04-08 18:31:49 +00:00
Code Activity

rm dealing with doubles

Browse source
This commit is contained in:
Lev Nachmanson 2023-03-06 10:23:52 -08:00
parent 5379fb820c
commit 547254abe7
6 changed files with 21 additions and 245 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

3
src/math/lp/lp_core_solver_base.cpp
View file

@ -83,9 +83,6 @@ 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;

1
src/math/lp/lp_core_solver_base.h
View file

@ -416,7 +416,6 @@ 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&);

17
src/math/lp/lp_core_solver_base_def.h
View file

@ -549,23 +549,6 @@ 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);

137
src/math/lp/lp_primal_core_solver.h
View file

@ -43,20 +43,16 @@ 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();
@ -277,14 +273,9 @@ public:
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);
@ -318,10 +309,7 @@ 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<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
@ -514,10 +502,7 @@ 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(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>();
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 @@ public:
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;
}
@ -555,16 +531,7 @@ public:
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;
}
@ -576,82 +543,32 @@ 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(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<X>();
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<X>();
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<X>()) {
theta = zero_of_type<X>();
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<X>()) {
theta = zero_of_type<X>();
}
lp_assert(false);
}
// j is a basic column or the entering, in any case x[j] has to stay feasible.
@ -722,14 +639,12 @@ 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
@ -743,13 +658,6 @@ 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<T>();
this->remove_column_from_inf_set(leaving);
}
}
void init_inf_set() {
this->clear_inf_set();
@ -885,15 +793,10 @@ 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<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);
}
@ -919,9 +822,7 @@ public:
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)) {
upper_bound_values) {
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);

103
src/math/lp/lp_primal_core_solver_def.h
View file

@ -228,54 +228,8 @@ int lp_primal_core_solver<T, X>::choose_entering_column(unsigned number_of_benef
}
template <typename T, typename X> bool lp_primal_core_solver<T, X>::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<X>(theta)) break;
}
return unlimited;
}
template <typename T, typename X> int lp_primal_core_solver<T, X>::
find_leaving_on_harris_theta(X const & harris_theta, X & t) {
int leaving = -1;
T pivot_abs_max = zero_of_type<T>();
// 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 <typename T, typename X> bool lp_primal_core_solver<T, X>::try_jump_to_another_bound_on_entering(unsigned entering,
const X & theta,
@ -398,17 +352,6 @@ template <typename T, typename X> int lp_primal_core_solver<T, X>::find_leaving_
}
template <typename T, typename X> int lp_primal_core_solver<T, X>::find_leaving_and_t(unsigned entering, X & t) {
X theta = zero_of_type<X>();
bool unlimited = get_harris_theta(theta);
lp_assert(unlimited || theta >= zero_of_type<X>());
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 <typename T, typename X> void lp_primal_core_solver<T, X>::init_column_
}
}
// debug only
template <typename T, typename X> T lp_primal_core_solver<T, X>::calculate_column_norm_exactly(unsigned j) {
lp_assert(false);
}
template <typename T, typename X> void lp_primal_core_solver<T, X>::update_or_init_column_norms(unsigned entering, unsigned leaving) {
lp_assert(numeric_traits<T>::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 <typename T, typename X> void lp_primal_core_solver<T, X>::update_column_norms(unsigned entering, unsigned leaving) {
lp_assert(numeric_traits<T>::precise() == false);
T pivot = this->m_pivot_row[entering];
T g_ent = calculate_norm_of_entering_exactly() / pivot / pivot;
if (!numeric_traits<T>::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 <typename T, typename X> T lp_primal_core_solver<T, X>::calculate_norm_of_entering_exactly() {
T r = numeric_traits<T>::one();
@ -771,13 +677,6 @@ template <typename T, typename X> bool lp_primal_core_solver<T, X>::done() {
return false;
}
template <typename T, typename X>
void lp_primal_core_solver<T, X>::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 <typename T, typename X>
void lp_primal_core_solver<T, X>::init_infeasibility_costs() {

5
src/math/lp/lp_primal_core_solver_tableau_def.h
View file

@ -295,10 +295,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::init_run_tab
if (this->m_settings.backup_costs)
backup_and_normalize_costs();
m_epsilon_of_reduced_cost = numeric_traits<X>::precise() ? zero_of_type<T>() : T(1) / T(10000000);
if (!numeric_traits<X>::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());