mirrors/z3
3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2025-04-22 16:45:31 +00:00
Code Activity

fix a bug in the lar_solver::m_status update during push/pop

Browse source 
Signed-off-by: Lev Nachmanson <levnach@microsoft.com>

progress in gomory cut

Signed-off-by: Lev Nachmanson <levnach@microsoft.com>

the first version of Gomory cut, probably broken

Signed-off-by: Lev Nachmanson <levnach@microsoft.com>

rename a function

Signed-off-by: Lev Nachmanson <levnach@microsoft.com>

gomory cut worked on a toy example

Signed-off-by: Lev Nachmanson <levnach@microsoft.com>

track the set of integer variables that are not set to integer values

Signed-off-by: Lev Nachmanson <levnach@microsoft.com>
This commit is contained in:
Lev Nachmanson 2017-07-10 16:34:23 -07:00 • committed by Lev Nachmanson
parent 0164ea9abb
commit aba7dcab3e
21 changed files with 682 additions and 455 deletions
Download patch file Download diff file Expand all files Collapse all files

2
src/shell/lp_frontend.cpp
View file

@ -84,7 +84,7 @@ void run_solver(lp_params & params, char const * mps_file_name) {
solver->find_maximal_solution();
*(solver->settings().get_message_ostream()) << "status is " << lp_status_to_string(solver->get_status()) << std::endl;
if (solver->get_status() == lp::OPTIMAL) {
if (solver->get_status() == lp::lp_status::OPTIMAL) {
if (params.min()) {
solver->flip_costs();
}

3
src/smt/theory_lra.cpp
View file

@ -107,6 +107,7 @@ namespace lp_api {
unsigned m_bound_propagations1;
unsigned m_bound_propagations2;
unsigned m_assert_diseq;
unsigned m_gomory_cuts;
stats() { reset(); }
void reset() {
memset(this, 0, sizeof(*this));
@ -1174,6 +1175,7 @@ namespace smt {
}
final_check_status final_check_eh() {
TRACE("lar_solver", tout << "ddd=" <<++lp::lp_settings::ddd << std::endl;);
m_use_nra_model = false;
lbool is_sat = l_true;
if (m_solver->get_status() != lp::lp_status::OPTIMAL) {
@ -1256,6 +1258,7 @@ namespace smt {
return l_false;
}
case lp::lia_move::cut: {
++m_stats.m_gomory_cuts;
// m_explanation implies term <= k
app_ref b = mk_bound(term, k);
m_eqs.reset();

24
src/test/lp/lp.cpp
View file

@ -1672,7 +1672,7 @@ void solve_some_mps(argument_parser & args_parser) {
}
if (!solve_for_rational) {
solve_mps(file_names[6], false, 0, time_limit, false, dual, compare_with_primal, args_parser);
solve_mps_with_known_solution(file_names[3], nullptr, INFEASIBLE, dual); // chvatal: 135(d)
solve_mps_with_known_solution(file_names[3], nullptr, lp_status::INFEASIBLE, dual); // chvatal: 135(d)
std::unordered_map<std::string, double> sol;
sol["X1"] = 0;
sol["X2"] = 6;
@ -1682,8 +1682,8 @@ void solve_some_mps(argument_parser & args_parser) {
sol["X6"] = 1;
sol["X7"] = 1;
sol["X8"] = 0;
solve_mps_with_known_solution(file_names[9], &sol, OPTIMAL, dual);
solve_mps_with_known_solution(file_names[0], &sol, OPTIMAL, dual);
solve_mps_with_known_solution(file_names[9], &sol, lp_status::OPTIMAL, dual);
solve_mps_with_known_solution(file_names[0], &sol, lp_status::OPTIMAL, dual);
sol.clear();
sol["X1"] = 25.0/14.0;
// sol["X2"] = 0;
@ -1692,10 +1692,10 @@ void solve_some_mps(argument_parser & args_parser) {
// sol["X5"] = 0;
// sol["X6"] = 0;
// sol["X7"] = 9.0/14.0;
solve_mps_with_known_solution(file_names[5], &sol, OPTIMAL, dual); // chvatal: 135(e)
solve_mps_with_known_solution(file_names[4], &sol, OPTIMAL, dual); // chvatal: 135(e)
solve_mps_with_known_solution(file_names[2], nullptr, UNBOUNDED, dual); // chvatal: 135(c)
solve_mps_with_known_solution(file_names[1], nullptr, UNBOUNDED, dual); // chvatal: 135(b)
solve_mps_with_known_solution(file_names[5], &sol, lp_status::OPTIMAL, dual); // chvatal: 135(e)
solve_mps_with_known_solution(file_names[4], &sol, lp_status::OPTIMAL, dual); // chvatal: 135(e)
solve_mps_with_known_solution(file_names[2], nullptr, lp_status::UNBOUNDED, dual); // chvatal: 135(c)
solve_mps_with_known_solution(file_names[1], nullptr, lp_status::UNBOUNDED, dual); // chvatal: 135(b)
solve_mps(file_names[8], false, 0, time_limit, false, dual, compare_with_primal, args_parser);
// return;
for (auto& s : file_names) {
@ -1761,7 +1761,7 @@ void solve_rational() {
expected_sol["x7"] = lp::mpq(1);
expected_sol["x8"] = lp::mpq(0);
solver.find_maximal_solution();
lp_assert(solver.get_status() == OPTIMAL);
lp_assert(solver.get_status() == lp_status::OPTIMAL);
for (auto it : expected_sol) {
lp_assert(it.second == solver.get_column_value_by_name(it.first));
}
@ -2449,12 +2449,12 @@ void run_lar_solver(argument_parser & args_parser, lar_solver * solver, mps_read
sw.start();
lp_status status = solver->solve();
std::cout << "status is " << lp_status_to_string(status) << ", processed for " << sw.get_current_seconds() <<" seconds, and " << solver->get_total_iterations() << " iterations" << std::endl;
if (solver->get_status() == INFEASIBLE) {
if (solver->get_status() == lp_status::INFEASIBLE) {
vector<std::pair<lp::mpq, constraint_index>> evidence;
solver->get_infeasibility_explanation(evidence);
}
if (args_parser.option_is_used("--randomize_lar")) {
if (solver->get_status() != OPTIMAL) {
if (solver->get_status() != lp_status::OPTIMAL) {
std::cout << "cannot check randomize on an infeazible problem" << std::endl;
return;
}
@ -2733,7 +2733,7 @@ void test_term() {
auto status = solver.solve();
std::cout << lp_status_to_string(status) << std::endl;
std::unordered_map<var_index, mpq> model;
if (status != OPTIMAL)
if (status != lp_status::OPTIMAL)
return;
solver.get_model(model);
@ -2761,7 +2761,7 @@ void test_evidence_for_total_inf_simple(argument_parser & args_parser) {
auto status = solver.solve();
std::cout << lp_status_to_string(status) << std::endl;
std::unordered_map<var_index, mpq> model;
lp_assert(solver.get_status() == INFEASIBLE);
lp_assert(solver.get_status() == lp_status::INFEASIBLE);
}
void test_bound_propagation_one_small_sample1() {
/*

489
src/util/lp/int_solver.cpp
View file

@ -19,9 +19,8 @@ void int_solver::fix_non_base_columns() {
}
if (!change)
return;
if (m_lar_solver->find_feasible_solution() == INFEASIBLE)
if (m_lar_solver->find_feasible_solution() == lp_status::INFEASIBLE)
failed();
init_inf_int_set();
lp_assert(is_feasible() && inf_int_set_is_correct());
}
@ -59,14 +58,22 @@ void int_solver::trace_inf_rows() const {
tout << "num of int infeasible: " << num << "\n";
}
int_set& int_solver::inf_int_set() {
return m_lar_solver->m_inf_int_set;
}
const int_set& int_solver::inf_int_set() const {
return m_lar_solver->m_inf_int_set;
}
int int_solver::find_inf_int_base_column() {
if (m_inf_int_set.is_empty())
if (inf_int_set().is_empty())
return -1;
int j = find_inf_int_boxed_base_column_with_smallest_range();
if (j != -1)
return j;
unsigned k = settings().random_next() % m_inf_int_set.m_index.size();
return m_inf_int_set.m_index[k];
unsigned k = settings().random_next() % inf_int_set().m_index.size();
return inf_int_set().m_index[k];
}
int int_solver::find_inf_int_boxed_base_column_with_smallest_range() {
@ -77,7 +84,7 @@ int int_solver::find_inf_int_boxed_base_column_with_smallest_range() {
unsigned n = 0;
lar_core_solver & lcs = m_lar_solver->m_mpq_lar_core_solver;
for (int j : m_inf_int_set.m_index) {
for (int j : inf_int_set().m_index) {
lp_assert(is_base(j) && column_is_int_inf(j));
if (!is_boxed(j))
continue;
@ -108,197 +115,263 @@ int int_solver::find_inf_int_boxed_base_column_with_smallest_range() {
}
bool int_solver::mk_gomory_cut(unsigned row_index, explanation & ex) {
lp_assert(false);
return true;
/*
const auto & row = m_lar_solver->A_r().m_rows[row_index];
// The following assertion is wrong. It may be violated in mixed-integer problems.
// SASSERT(!all_coeff_int(r));
theory_var x_i = r.get_base_var();
SASSERT(is_int(x_i));
// The following assertion is wrong. It may be violated in mixed-real-interger problems.
// The check is_gomory_cut_target will discard rows where any variable contains infinitesimals.
// SASSERT(m_value[x_i].is_rational()); // infinitesimals are not used for integer variables
SASSERT(!m_value[x_i].is_int()); // the base variable is not assigned to an integer value.
bool int_solver::is_gomory_cut_target() {
m_iter_on_gomory_row->reset();
unsigned j;
TRACE("gomory_cut", m_lar_solver->print_linear_iterator(m_iter_on_gomory_row, tout);
m_iter_on_gomory_row->reset();
);
if (constrain_free_vars(r) || !is_gomory_cut_target(r)) {
TRACE("gomory_cut", tout << "failed to apply gomory cut:\n";
tout << "constrain_free_vars(r): " << constrain_free_vars(r) << "\n";);
return false;
while (m_iter_on_gomory_row->next(j)) {
// All non base variables must be at their bounds and assigned to rationals (that is, infinitesimals are not allowed).
if (j != m_gomory_cut_inf_column && (!at_bound(j) || !is_zero(get_value(j).y))) {
TRACE("gomory_cut", tout << "row is not gomory cut target:\n";
display_column(tout, j);
tout << "at_bound: " << at_bound(j) << "\n";
tout << "infinitesimal: " << !is_zero(get_value(j).y) << "\n";);
return false;
}
}
m_iter_on_gomory_row->reset();
return true;
}
TRACE("gomory_cut", tout << "applying cut at:\n"; display_row_info(tout, r););
void int_solver::real_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& pol, explanation & expl) {
mpq f_0 = fractional_part(get_value(m_gomory_cut_inf_column));
mpq new_a;
if (at_lower(x_j)) {
if (a.is_pos()) {
new_a = a / (1 - f_0);
}
else {
new_a = a / f_0;
new_a.neg();
}
k += lower_bound(x_j).x * k; // k.addmul(new_a, lower_bound(x_j).x); // is it a faster operation
antecedents ante(*this);
expl.push_justification(column_low_bound_constraint(x_j), new_a);
}
else {
lp_assert(at_upper(x_j));
if (a.is_pos()) {
new_a = a / f_0;
new_a.neg(); // the upper terms are inverted.
}
else {
new_a = a / (mpq(1) - f_0);
}
k += upper_bound(x_j).x * k; // k.addmul(new_a, upper_bound(x_j).get_rational());
expl.push_justification(column_upper_bound_constraint(x_j), new_a);
}
TRACE("gomory_cut_detail", tout << a << "*v" << x_j << " k: " << k << "\n";);
pol.add_monoid(new_a, x_j);
}
m_stats.m_gomory_cuts++;
constraint_index int_solver::column_upper_bound_constraint(unsigned j) const {
return m_lar_solver->get_column_upper_bound_witness(j);
}
// gomory will be pol >= k
numeral k(1);
buffer<row_entry> pol;
numeral f_0 = Ext::fractional_part(m_value[x_i]);
numeral one_minus_f_0 = numeral(1) - f_0;
SASSERT(!f_0.is_zero());
SASSERT(!one_minus_f_0.is_zero());
numeral lcm_den(1);
unsigned num_ints = 0;
constraint_index int_solver::column_low_bound_constraint(unsigned j) const {
return m_lar_solver->get_column_low_bound_witness(j);
}
typename vector<row_entry>::const_iterator it = r.begin_entries();
typename vector<row_entry>::const_iterator end = r.end_entries();
for (; it != end; ++it) {
if (!it->is_dead() && it->m_var != x_i) {
theory_var x_j = it->m_var;
numeral a_ij = it->m_coeff;
a_ij.neg(); // make the used format compatible with the format used in: Integrating Simplex with DPLL(T)
if (is_real(x_j)) {
numeral new_a_ij;
if (at_lower(x_j)) {
if (a_ij.is_pos()) {
new_a_ij = a_ij / one_minus_f_0;
}
else {
new_a_ij = a_ij / f_0;
new_a_ij.neg();
}
k.addmul(new_a_ij, lower_bound(x_j).get_rational());
lower(x_j)->push_justification(ante, new_a_ij, coeffs_enabled());
}
else {
SASSERT(at_upper(x_j));
if (a_ij.is_pos()) {
new_a_ij = a_ij / f_0;
new_a_ij.neg(); // the upper terms are inverted.
}
else {
new_a_ij = a_ij / one_minus_f_0;
}
k.addmul(new_a_ij, upper_bound(x_j).get_rational());
upper(x_j)->push_justification(ante, new_a_ij, coeffs_enabled());
}
TRACE("gomory_cut_detail", tout << a_ij << "*v" << x_j << " k: " << k << "\n";);
pol.push_back(row_entry(new_a_ij, x_j));
void int_solver::int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term & pol, explanation& expl, mpq & lcm_den) {
mpq f_0 = fractional_part(get_value(m_gomory_cut_inf_column));
lp_assert(is_int(x_j));
mpq f_j = fractional_part(a);
TRACE("gomory_cut_detail",
tout << a << "*v" << x_j << "\n";
tout << "fractional_part: " << fractional_part(a) << "\n";
tout << "f_j: " << f_j << "\n";
tout << "f_0: " << f_0 << "\n";
tout << "one_minus_f_0: " << 1 - f_0 << "\n";);
if (!f_j.is_zero()) {
mpq new_a;
if (at_lower(x_j)) {
auto one_minus_f_0 = 1 - f_0;
if (f_j <= one_minus_f_0) {
new_a = f_j / one_minus_f_0;
}
else {
++num_ints;
SASSERT(is_int(x_j));
numeral f_j = Ext::fractional_part(a_ij);
TRACE("gomory_cut_detail",
tout << a_ij << "*v" << x_j << "\n";
tout << "fractional_part: " << Ext::fractional_part(a_ij) << "\n";
tout << "f_j: " << f_j << "\n";
tout << "f_0: " << f_0 << "\n";
tout << "one_minus_f_0: " << one_minus_f_0 << "\n";);
if (!f_j.is_zero()) {
numeral new_a_ij;
if (at_lower(x_j)) {
if (f_j <= one_minus_f_0) {
new_a_ij = f_j / one_minus_f_0;
}
else {
new_a_ij = (numeral(1) - f_j) / f_0;
}
k.addmul(new_a_ij, lower_bound(x_j).get_rational());
lower(x_j)->push_justification(ante, new_a_ij, coeffs_enabled());
}
else {
SASSERT(at_upper(x_j));
if (f_j <= f_0) {
new_a_ij = f_j / f_0;
}
else {
new_a_ij = (numeral(1) - f_j) / one_minus_f_0;
}
new_a_ij.neg(); // the upper terms are inverted
k.addmul(new_a_ij, upper_bound(x_j).get_rational());
upper(x_j)->push_justification(ante, new_a_ij, coeffs_enabled());
}
TRACE("gomory_cut_detail", tout << "new_a_ij: " << new_a_ij << " k: " << k << "\n";);
pol.push_back(row_entry(new_a_ij, x_j));
lcm_den = lcm(lcm_den, denominator(new_a_ij));
}
new_a = (1 - f_j) / f_0;
}
k.addmul(new_a, lower_bound(x_j).x);
expl.push_justification(column_low_bound_constraint(x_j), new_a);
}
}
CTRACE("empty_pol", pol.empty(), display_row_info(tout, r););
expr_ref bound(get_manager());
if (pol.empty()) {
SASSERT(k.is_pos());
// conflict 0 >= k where k is positive
set_conflict(ante, ante, "gomory-cut");
return true;
}
else if (pol.size() == 1) {
theory_var v = pol[0].m_var;
k /= pol[0].m_coeff;
bool is_lower = pol[0].m_coeff.is_pos();
if (is_int(v) && !k.is_int()) {
k = is_lower?ceil(k):floor(k);
else {
SASSERT(at_upper(x_j));
if (f_j <= f_0) {
new_a = f_j / f_0;
}
else {
new_a = (mpq(1) - f_j) / 1 - f_0;
}
new_a.neg(); // the upper terms are inverted
k.addmul(new_a, upper_bound(x_j).x);
expl.push_justification(column_upper_bound_constraint(x_j), new_a);
}
rational _k = k.to_rational();
if (is_lower)
bound = m_util.mk_ge(get_enode(v)->get_owner(), m_util.mk_numeral(_k, is_int(v)));
else
bound = m_util.mk_le(get_enode(v)->get_owner(), m_util.mk_numeral(_k, is_int(v)));
TRACE("gomory_cut_detail", tout << "new_a: " << new_a << " k: " << k << "\n";);
pol.add_monoid(new_a, x_j);
lcm_den = lcm(lcm_den, denominator(new_a));
}
else {
}
lia_move int_solver::report_conflict_from_gomory_cut(mpq & k) {
TRACE("empty_pol",
display_row_info(tout,
m_lar_solver->m_mpq_lar_core_solver.m_r_heading[m_gomory_cut_inf_column]););
lp_assert(k.is_pos());
// conflict 0 >= k where k is positive
k.neg(); // returning 0 <= -k
return lia_move::conflict;
}
void int_solver::gomory_cut_adjust_t_and_k_for_size_gt_1(
vector<std::pair<mpq, unsigned>> & pol,
lar_term & t,
mpq &k,
unsigned num_ints,
mpq & lcm_den) {
if (num_ints > 0) {
lcm_den = lcm(lcm_den, denominator(k));
TRACE("gomory_cut_detail", tout << "k: " << k << " lcm_den: " << lcm_den << "\n";
for (unsigned i = 0; i < pol.size(); i++) {
tout << pol[i].m_coeff << " " << pol[i].m_var << "\n";
}
tout << "k: " << k << "\n";);
SASSERT(lcm_den.is_pos());
linear_combination_iterator_on_vector<mpq> pi(pol);
m_lar_solver->print_linear_iterator(&pi, tout);
tout << "\nk: " << k << "\n";);
lp_assert(lcm_den.is_pos());
if (!lcm_den.is_one()) {
// normalize coefficients of integer parameters to be integers.
unsigned n = pol.size();
for (unsigned i = 0; i < n; i++) {
pol[i].m_coeff *= lcm_den;
SASSERT(!is_int(pol[i].m_var) || pol[i].m_coeff.is_int());
for (auto & pi: pol) {
pi.first *= lcm_den;
SASSERT(!is_int(pi.second) || pi.first.is_int());
}
k *= lcm_den;
}
TRACE("gomory_cut_detail", tout << "after *lcm\n";
TRACE("gomory_cut_detail", tout << "after *lcm_den\n";
for (unsigned i = 0; i < pol.size(); i++) {
tout << pol[i].m_coeff << " * v" << pol[i].m_var << "\n";
tout << pol[i].first << " * v" << pol[i].second << "\n";
}
tout << "k: " << k << "\n";);
}
mk_polynomial_ge(pol.size(), pol.c_ptr(), k.to_rational(), bound);
}
TRACE("gomory_cut", tout << "new cut:\n" << bound << "\n"; ante.display(tout););
literal l = null_literal;
context & ctx = get_context();
ctx.internalize(bound, true);
l = ctx.get_literal(bound);
ctx.mark_as_relevant(l);
dump_lemmas(l, ante);
ctx.assign(l, ctx.mk_justification(
gomory_cut_justification(
get_id(), ctx.get_region(),
ante.lits().size(), ante.lits().c_ptr(),
ante.eqs().size(), ante.eqs().c_ptr(), ante, l)));
return true;
*/
t.clear();
// negate everything to return -pol <= -k
for (const auto & pi: pol)
t.add_monoid(-pi.first, pi.second);
k.neg();
}
void int_solver::gomory_cut_adjust_t_and_k_for_size_1(const vector<std::pair<mpq, unsigned>> & pol, lar_term& t, mpq &k) {
lp_assert(pol.size() == 1);
unsigned j = pol[0].second;
k /= pol[0].first;
bool is_lower = pol[0].first.is_pos();
if (is_int(j) && !k.is_int()) {
k = is_lower?ceil(k):floor(k);
}
if (is_lower) { // returning -t <= -k which is equivalent to t >= k
k.neg();
t.negate();
}
}
lia_move int_solver::report_gomory_cut(lar_term& t, mpq& k, mpq &lcm_den, unsigned num_ints) {
lp_assert(!t.is_empty());
auto pol = t.coeffs_as_vector();
if (pol.size() == 1)
gomory_cut_adjust_t_and_k_for_size_1(pol, t, k);
else
gomory_cut_adjust_t_and_k_for_size_gt_1(pol, t, k, num_ints, lcm_den);
m_lar_solver->subs_terms_for_debugging(t); // todo: remove later
return lia_move::cut;
}
lia_move int_solver::mk_gomory_cut(lar_term& t, mpq& k, explanation & expl) {
lp_assert(column_is_int_inf(m_gomory_cut_inf_column));
TRACE("gomory_cut", tout << "applying cut at:\n"; m_lar_solver->print_linear_iterator(m_iter_on_gomory_row, tout); tout << std::endl; m_iter_on_gomory_row->reset(););
// gomory will be t >= k
k = 1;
mpq lcm_den(1);
unsigned num_ints = 0;
unsigned x_j;
mpq a;
while (m_iter_on_gomory_row->next(a, x_j)) {
if (x_j == m_gomory_cut_inf_column)
continue;
// make the format compatible with the format used in: Integrating Simplex with DPLL(T)
a.neg();
if (is_real(x_j))
real_case_in_gomory_cut(a, x_j, k, t, expl);
else {
num_ints++;
int_case_in_gomory_cut(a, x_j, k, t, expl, lcm_den);
}
}
if (t.is_empty())
return report_conflict_from_gomory_cut(k);
return report_gomory_cut(t, k, lcm_den, num_ints);
}
void int_solver::init_check_data() {
init_inf_int_set();
unsigned n = m_lar_solver->A_r().column_count();
m_old_values_set.resize(n);
m_old_values_data.resize(n);
}
int int_solver::find_next_free_var_in_gomory_row() {
lp_assert(m_iter_on_gomory_row != nullptr);
unsigned j;
while(m_iter_on_gomory_row->next(j)) {
if (j != m_gomory_cut_inf_column && is_free(j))
return static_cast<int>(j);
}
return -1;
}
lia_move int_solver::proceed_with_gomory_cut(lar_term& t, mpq& k, explanation& ex) {
int j = find_next_free_var_in_gomory_row();
if (j != -1) {
m_found_free_var_in_gomory_row = true;
lp_assert(t.is_empty());
t.add_monoid(mpq(1), j);
k = zero_of_type<mpq>();
return lia_move::branch; // branch on a free column
}
if (m_found_free_var_in_gomory_row || !is_gomory_cut_target()) {
m_found_free_var_in_gomory_row = false;
delete m_iter_on_gomory_row;
m_iter_on_gomory_row = nullptr;
return lia_move::continue_with_check;
}
lia_move ret = mk_gomory_cut(t, k, ex);
delete m_iter_on_gomory_row;
m_iter_on_gomory_row = nullptr;
return ret;
}
lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex) {
lp_assert(m_lar_solver->m_mpq_lar_core_solver.r_basis_is_OK());
lp_assert(is_feasible());
if (m_iter_on_gomory_row != nullptr) {
auto ret = proceed_with_gomory_cut(t, k, ex);
TRACE("gomory_cut", tout << "term t = "; m_lar_solver->print_term_as_indices(t, tout););
if (ret != lia_move::continue_with_check)
return ret;
}
init_check_data();
lp_assert(inf_int_set_is_correct());
// currently it is a reimplementation of
@ -327,11 +400,13 @@ lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex) {
if ((++m_branch_cut_counter) % settings().m_int_branch_cut_threshold == 0) {
move_non_base_vars_to_bounds();
move_non_base_vars_to_bounds(); // todo track changed variables
lp_status st = m_lar_solver->find_feasible_solution();
if (st != lp_status::FEASIBLE && st != lp_status::OPTIMAL) {
return lia_move::give_up;
}
lp_assert(inf_int_set_is_correct());
// init_inf_int_set(); // todo - can we avoid this call?
int j = find_inf_int_base_column();
if (j != -1) {
TRACE("arith_int", tout << "j = " << j << " does not have an integer assignment: " << get_value(j) << "\n";);
@ -349,7 +424,7 @@ lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex) {
TRACE("arith_int", tout << "j" << j << " does not have an integer assignment: " << get_value(j) << "\n";);
lp_assert(t.is_empty());
t.add_to_map(j, mpq(1));
t.add_monoid(1, j);
k = floor(get_value(j));
TRACE("arith_int", tout << "branching v" << j << " = " << get_value(j) << "\n";
display_column(tout, j);
@ -754,6 +829,10 @@ bool int_solver::is_int(unsigned j) const {
return m_lar_solver->column_is_int(j);
}
bool int_solver::is_real(unsigned j) const {
return !is_int(j);
}
bool int_solver::value_is_int(unsigned j) const {
return m_lar_solver->m_mpq_lar_core_solver.m_r_x[j].is_int();
}
@ -777,7 +856,7 @@ void int_solver::display_column(std::ostream & out, unsigned j) const {
bool int_solver::inf_int_set_is_correct() const {
for (unsigned j = 0; j < m_lar_solver->A_r().column_count(); j++) {
if (m_inf_int_set.contains(j) != is_int(j) && (!value_is_int(j)))
if (inf_int_set().contains(j) != (is_int(j) && (!value_is_int(j))))
return false;
}
return true;
@ -787,20 +866,11 @@ bool int_solver::column_is_int_inf(unsigned j) const {
return is_int(j) && (!value_is_int(j));
}
void int_solver::init_inf_int_set() {
m_inf_int_set.clear();
m_inf_int_set.resize(m_lar_solver->A_r().column_count());
for (unsigned j : m_lar_solver->m_mpq_lar_core_solver.m_r_basis) {
if (column_is_int_inf(j))
m_inf_int_set.insert(j);
}
}
void int_solver::update_column_in_inf_set_set(unsigned j) {
void int_solver::update_column_in_int_inf_set(unsigned j) {
if (is_int(j) && (!value_is_int(j)))
m_inf_int_set.insert(j);
inf_int_set().insert(j);
else
m_inf_int_set.erase(j);
inf_int_set().erase(j);
}
bool int_solver::is_base(unsigned j) const {
@ -811,8 +881,73 @@ bool int_solver::is_boxed(unsigned j) const {
return m_lar_solver->m_mpq_lar_core_solver.m_column_types[j] == column_type::boxed;
}
bool int_solver::is_free(unsigned j) const {
return m_lar_solver->m_mpq_lar_core_solver.m_column_types[j] == column_type::free_column;
}
bool int_solver::at_bound(unsigned j) const {
auto & mpq_solver = m_lar_solver->m_mpq_lar_core_solver.m_r_solver;
switch (mpq_solver.m_column_types[j] ) {
case column_type::fixed:
case column_type::boxed:
return
mpq_solver.m_low_bounds[j] == get_value(j) ||
mpq_solver.m_upper_bounds[j] == get_value(j);
case column_type::low_bound:
return mpq_solver.m_low_bounds[j] == get_value(j);
case column_type::upper_bound:
return mpq_solver.m_upper_bounds[j] == get_value(j);
default:
return false;
}
}
bool int_solver::at_lower(unsigned j) const {
auto & mpq_solver = m_lar_solver->m_mpq_lar_core_solver.m_r_solver;
switch (mpq_solver.m_column_types[j] ) {
case column_type::fixed:
case column_type::boxed:
case column_type::low_bound:
return mpq_solver.m_low_bounds[j] == get_value(j);
default:
return false;
}
}
bool int_solver::at_upper(unsigned j) const {
auto & mpq_solver = m_lar_solver->m_mpq_lar_core_solver.m_r_solver;
switch (mpq_solver.m_column_types[j] ) {
case column_type::fixed:
case column_type::boxed:
case column_type::upper_bound:
return mpq_solver.m_upper_bounds[j] == get_value(j);
default:
return false;
}
}
lp_settings& int_solver::settings() {
return m_lar_solver->settings();
}
void int_solver::display_row_info(std::ostream & out, unsigned row_index) const {
auto & rslv = m_lar_solver->m_mpq_lar_core_solver.m_r_solver;
auto it = m_lar_solver->get_iterator_on_row(row_index);
mpq a;
unsigned j;
while (it->next(a, j)) {
if (numeric_traits<mpq>::is_pos(a))
out << "+";
out << a << rslv.column_name(j) << " ";
}
it->reset();
while(it->next(j)) {
rslv.print_column_bound_info(j, out);
}
rslv.print_column_bound_info(rslv.m_basis[row_index], out);
delete it;
}
}

43
src/util/lp/int_solver.h
View file

@ -8,7 +8,6 @@
#include "util/lp/iterator_on_row.h"
#include "util/lp/int_set.h"
#include "util/lp/lar_term.h"
namespace lp {
class lar_solver;
template <typename T, typename X>
@ -23,6 +22,9 @@ enum class lia_move {
struct explanation {
vector<std::pair<mpq, constraint_index>> m_explanation;
void push_justification(constraint_index j, const mpq& v) {
m_explanation.push_back(std::make_pair(v, j));
}
};
class int_solver {
@ -31,10 +33,15 @@ public:
lar_solver *m_lar_solver;
int_set m_old_values_set;
vector<impq> m_old_values_data;
int_set m_inf_int_set;
unsigned m_branch_cut_counter;
linear_combination_iterator<mpq>* m_iter_on_gomory_row;
unsigned m_gomory_cut_inf_column;
bool m_found_free_var_in_gomory_row;
// methods
int_solver(lar_solver* lp);
int_set& inf_int_set();
const int_set& inf_int_set() const;
// main function to check that solution provided by lar_solver is valid for integral values,
// or provide a way of how it can be adjusted.
lia_move check(lar_term& t, mpq& k, explanation& ex);
@ -78,6 +85,7 @@ private:
const impq & lower_bound(unsigned j) const;
const impq & upper_bound(unsigned j) const;
bool is_int(unsigned j) const;
bool is_real(unsigned j) const;
bool is_base(unsigned j) const;
bool is_boxed(unsigned j) const;
bool value_is_int(unsigned j) const;
@ -88,8 +96,7 @@ private:
const impq & get_value(unsigned j) const;
void display_column(std::ostream & out, unsigned j) const;
bool inf_int_set_is_correct() const;
void init_inf_int_set();
void update_column_in_inf_set_set(unsigned j);
void update_column_in_int_inf_set(unsigned j);
bool column_is_int_inf(unsigned j) const;
void trace_inf_rows() const;
int find_inf_int_base_column();
@ -97,7 +104,33 @@ private:
lp_settings& settings();
void move_non_base_vars_to_bounds();
void branch_infeasible_int_var(unsigned);
bool mk_gomory_cut(unsigned row_index, explanation & ex);
lia_move mk_gomory_cut(lar_term& t, mpq& k,explanation & ex);
lia_move report_conflict_from_gomory_cut(mpq & k);
lia_move report_gomory_cut(lar_term& t, mpq& k, mpq& lcm_den, unsigned num_ints);
void init_check_data();
bool constrain_free_vars(linear_combination_iterator<mpq> * r);
lia_move proceed_with_gomory_cut(lar_term& t, mpq& k, explanation& ex);
int find_next_free_var_in_gomory_row();
bool is_gomory_cut_target();
bool at_bound(unsigned j) const;
bool at_lower(unsigned j) const;
bool at_upper(unsigned j) const;
inline static bool is_rational(const impq & n) {
return is_zero(n.y);
}
inline static
mpq fractional_part(const impq & n) {
lp_assert(is_rational);
return n.x - floor(n.x);
}
void real_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& t, explanation & ex);
void int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& t, explanation& ex, mpq & lcm_den);
constraint_index column_upper_bound_constraint(unsigned j) const;
constraint_index column_low_bound_constraint(unsigned j) const;
void display_row_info(std::ostream & out, unsigned row_index) const;
void gomory_cut_adjust_t_and_k_for_size_1(const vector<std::pair<mpq, unsigned>> & pol, lar_term & t, mpq &k);
void gomory_cut_adjust_t_and_k_for_size_gt_1(vector<std::pair<mpq, unsigned>> & pol, lar_term & t, mpq &k, unsigned num_ints, mpq &lcm_den);
};
}

12
src/util/lp/lar_core_solver.hpp
View file

@ -260,7 +260,7 @@ void lar_core_solver::solve() {
lp_assert(m_r_solver.non_basic_columns_are_set_correctly());
lp_assert(m_r_solver.inf_set_is_correct());
if (m_r_solver.current_x_is_feasible() && m_r_solver.m_look_for_feasible_solution_only) {
m_r_solver.set_status(OPTIMAL);
m_r_solver.set_status(lp_status::OPTIMAL);
return;
}
++settings().st().m_need_to_solve_inf;
@ -270,8 +270,8 @@ void lar_core_solver::solve() {
prefix_d();
lar_solution_signature solution_signature;
vector<unsigned> changes_of_basis = find_solution_signature_with_doubles(solution_signature);
if (m_d_solver.get_status() == TIME_EXHAUSTED) {
m_r_solver.set_status(TIME_EXHAUSTED);
if (m_d_solver.get_status() == lp_status::TIME_EXHAUSTED) {
m_r_solver.set_status(lp_status::TIME_EXHAUSTED);
return;
}
if (settings().use_tableau())
@ -292,10 +292,10 @@ void lar_core_solver::solve() {
m_r_solver.solve();
lp_assert(!settings().use_tableau() || r_basis_is_OK());
}
if (m_r_solver.get_status() == INFEASIBLE) {
if (m_r_solver.get_status() == lp_status::INFEASIBLE) {
fill_not_improvable_zero_sum();
} else if (m_r_solver.get_status() != UNBOUNDED) {
m_r_solver.set_status(OPTIMAL);
} else if (m_r_solver.get_status() != lp_status::UNBOUNDED) {
m_r_solver.set_status(lp_status::OPTIMAL);
}
lp_assert(r_basis_is_OK());
lp_assert(m_r_solver.non_basic_columns_are_set_correctly());

144
src/util/lp/lar_solver.cpp
View file

@ -27,10 +27,20 @@ void clear() {lp_assert(false); // not implemented
}
lar_solver::lar_solver() : m_status(OPTIMAL),
lar_solver::lar_solver() : m_status(lp_status::OPTIMAL),
m_infeasible_column_index(-1),
m_terms_start_index(1000000),
m_mpq_lar_core_solver(m_settings, *this)
m_mpq_lar_core_solver(m_settings, *this),
m_tracker_of_x_change([&](unsigned j, const impq & x){
if (!var_is_int(j)) {
lp_assert(m_inf_int_set.contains(j) == false);
return;
}
if (m_mpq_lar_core_solver.m_r_x[j].is_int())
m_inf_int_set.erase(j);
else
m_inf_int_set.insert(j);
})
{
}
@ -237,7 +247,7 @@ void lar_solver::explain_implied_bound(implied_bound & ib, bound_propagator & bp
}
int a_sign = is_pos(a)? 1: -1;
int sign = j_sign * a_sign;
const ul_pair & ul = m_vars_to_ul_pairs[j];
const ul_pair & ul = m_columns_to_ul_pairs[j];
auto witness = sign > 0? ul.upper_bound_witness(): ul.low_bound_witness();
lp_assert(is_valid(witness));
bp.consume(a, witness);
@ -279,6 +289,10 @@ lp_status lar_solver::get_status() const { return m_status;}
void lar_solver::set_status(lp_status s) {m_status = s;}
bool lar_solver::has_int_var() const {
return m_mpq_lar_core_solver.m_r_solver.m_tracker_of_x_change != nullptr;
}
lp_status lar_solver::find_feasible_solution() {
m_settings.st().m_make_feasible++;
if (A_r().column_count() > m_settings.st().m_max_cols)
@ -288,16 +302,20 @@ lp_status lar_solver::find_feasible_solution() {
if (strategy_is_undecided())
decide_on_strategy_and_adjust_initial_state();
if (has_int_var()) {
m_inf_int_set.resize(A_r().column_count());
}
m_mpq_lar_core_solver.m_r_solver.m_look_for_feasible_solution_only = true;
return solve();
}
lp_status lar_solver::solve() {
if (m_status == INFEASIBLE) {
if (m_status == lp_status::INFEASIBLE) {
return m_status;
}
solve_with_core_solver();
if (m_status != INFEASIBLE) {
if (m_status != lp_status::INFEASIBLE) {
if (m_settings.bound_propagation())
detect_rows_with_changed_bounds();
}
@ -309,7 +327,7 @@ lp_status lar_solver::solve() {
void lar_solver::fill_explanation_from_infeasible_column(vector<std::pair<mpq, constraint_index>> & evidence) const{
// this is the case when the lower bound is in conflict with the upper one
const ul_pair & ul = m_vars_to_ul_pairs[m_infeasible_column_index];
const ul_pair & ul = m_columns_to_ul_pairs[m_infeasible_column_index];
evidence.push_back(std::make_pair(numeric_traits<mpq>::one(), ul.upper_bound_witness()));
evidence.push_back(std::make_pair(-numeric_traits<mpq>::one(), ul.low_bound_witness()));
}
@ -325,8 +343,7 @@ vector<unsigned> lar_solver::get_list_of_all_var_indices() const {
void lar_solver::push() {
m_simplex_strategy = m_settings.simplex_strategy();
m_simplex_strategy.push();
m_status.push();
m_vars_to_ul_pairs.push();
m_columns_to_ul_pairs.push();
m_infeasible_column_index.push();
m_mpq_lar_core_solver.push();
m_term_count = m_terms.size();
@ -335,7 +352,7 @@ void lar_solver::push() {
m_constraint_count.push();
}
void lar_solver::clean_large_elements_after_pop(unsigned n, int_set& set) {
void lar_solver::clean_popped_elements(unsigned n, int_set& set) {
vector<int> to_remove;
for (unsigned j: set.m_index)
if (j >= n)
@ -345,29 +362,31 @@ void lar_solver::clean_large_elements_after_pop(unsigned n, int_set& set) {
}
void lar_solver::shrink_inf_set_after_pop(unsigned n, int_set & set) {
clean_large_elements_after_pop(n, set);
clean_popped_elements(n, set);
set.resize(n);
}
void lar_solver::pop(unsigned k) {
TRACE("lar_solver", tout << "k = " << k << std::endl;);
int n_was = static_cast<int>(m_ext_vars_to_columns.size());
m_status.pop(k);
m_infeasible_column_index.pop(k);
unsigned n = m_vars_to_ul_pairs.peek_size(k);
unsigned n = m_columns_to_ul_pairs.peek_size(k);
for (unsigned j = n_was; j-- > n;)
m_ext_vars_to_columns.erase(m_columns_to_ext_vars_or_term_indices[j]);
m_columns_to_ext_vars_or_term_indices.resize(n);
if (m_settings.use_tableau()) {
pop_tableau();
}
m_vars_to_ul_pairs.pop(k);
m_columns_to_ul_pairs.pop(k);
m_mpq_lar_core_solver.pop(k);
clean_large_elements_after_pop(n, m_columns_with_changed_bound);
clean_popped_elements(n, m_columns_with_changed_bound);
unsigned m = A_r().row_count();
clean_large_elements_after_pop(m, m_rows_with_changed_bounds);
clean_popped_elements(m, m_rows_with_changed_bounds);
clean_inf_set_of_r_solver_after_pop();
m_status = m_mpq_lar_core_solver.m_r_solver.current_x_is_feasible()? lp_status::OPTIMAL: lp_status::UNKNOWN;
lp_assert(m_settings.simplex_strategy() == simplex_strategy_enum::undecided ||
(!use_tableau()) || m_mpq_lar_core_solver.m_r_solver.reduced_costs_are_correct_tableau());
@ -404,7 +423,7 @@ bool lar_solver::maximize_term_on_tableau(const vector<std::pair<mpq, var_index>
decide_on_strategy_and_adjust_initial_state();
m_mpq_lar_core_solver.solve();
if (m_mpq_lar_core_solver.m_r_solver.get_status() == UNBOUNDED)
if (m_mpq_lar_core_solver.m_r_solver.get_status() == lp_status::UNBOUNDED)
return false;
term_max = 0;
@ -482,7 +501,7 @@ bool lar_solver::maximize_term_on_corrected_r_solver(const vector<std::pair<mpq,
bool ret = maximize_term_on_tableau(term, term_max);
settings().simplex_strategy() = simplex_strategy_enum::tableau_rows;
set_costs_to_zero(term);
m_mpq_lar_core_solver.m_r_solver.set_status(OPTIMAL);
m_mpq_lar_core_solver.m_r_solver.set_status(lp_status::OPTIMAL);
return ret;
}
case simplex_strategy_enum::tableau_costs:
@ -490,7 +509,7 @@ bool lar_solver::maximize_term_on_corrected_r_solver(const vector<std::pair<mpq,
{
bool ret = maximize_term_on_tableau(term, term_max);
set_costs_to_zero(term);
m_mpq_lar_core_solver.m_r_solver.set_status(OPTIMAL);
m_mpq_lar_core_solver.m_r_solver.set_status(lp_status::OPTIMAL);
return ret;
}
@ -529,18 +548,18 @@ void lar_solver::pop_core_solver_params(unsigned k) {
void lar_solver::set_upper_bound_witness(var_index j, constraint_index ci) {
ul_pair ul = m_vars_to_ul_pairs[j];
ul_pair ul = m_columns_to_ul_pairs[j];
ul.upper_bound_witness() = ci;
m_vars_to_ul_pairs[j] = ul;
m_columns_to_ul_pairs[j] = ul;
}
void lar_solver::set_low_bound_witness(var_index j, constraint_index ci) {
ul_pair ul = m_vars_to_ul_pairs[j];
ul_pair ul = m_columns_to_ul_pairs[j];
ul.low_bound_witness() = ci;
m_vars_to_ul_pairs[j] = ul;
m_columns_to_ul_pairs[j] = ul;
}
void lar_solver::register_one_coeff_in_map(std::unordered_map<var_index, mpq> & coeffs, const mpq & a, unsigned j) {
void lar_solver::register_monoid_in_map(std::unordered_map<var_index, mpq> & coeffs, const mpq & a, unsigned j) {
auto it = coeffs.find(j);
if (it == coeffs.end()) {
coeffs[j] = a;
@ -551,18 +570,18 @@ void lar_solver::register_one_coeff_in_map(std::unordered_map<var_index, mpq> &
void lar_solver::substitute_terms_in_linear_expression(const vector<std::pair<mpq, var_index>>& left_side_with_terms,
vector<std::pair<mpq, var_index>> &left_side, mpq & right_side) const {
vector<std::pair<mpq, var_index>> &left_side, mpq & free_coeff) const {
std::unordered_map<var_index, mpq> coeffs;
for (auto & t : left_side_with_terms) {
unsigned j = t.second;
if (!is_term(j)) {
register_one_coeff_in_map(coeffs, t.first, j);
register_monoid_in_map(coeffs, t.first, j);
} else {
const lar_term & term = * m_terms[adjust_term_index(t.second)];
for (auto & p : term.coeffs()){
register_one_coeff_in_map(coeffs, t.first * p.second , p.first);
register_monoid_in_map(coeffs, t.first * p.second , p.first);
}
right_side += t.first * term.m_v;
free_coeff += t.first * term.m_v;
}
}
@ -732,7 +751,7 @@ void lar_solver::solve_with_core_solver() {
update_x_and_inf_costs_for_columns_with_changed_bounds();
m_mpq_lar_core_solver.solve();
set_status(m_mpq_lar_core_solver.m_r_solver.get_status());
lp_assert(m_status != OPTIMAL || all_constraints_hold());
lp_assert(m_status != lp_status::OPTIMAL || all_constraints_hold());
}
@ -1039,11 +1058,11 @@ mpq lar_solver::sum_of_right_sides_of_explanation(const vector<std::pair<mpq, un
bool lar_solver::has_lower_bound(var_index var, constraint_index& ci, mpq& value, bool& is_strict) {
if (var >= m_vars_to_ul_pairs.size()) {
if (var >= m_columns_to_ul_pairs.size()) {
// TBD: bounds on terms could also be used, caller may have to track these.
return false;
}
const ul_pair & ul = m_vars_to_ul_pairs[var];
const ul_pair & ul = m_columns_to_ul_pairs[var];
ci = ul.low_bound_witness();
if (ci != static_cast<constraint_index>(-1)) {
auto& p = m_mpq_lar_core_solver.m_r_low_bounds()[var];
@ -1058,11 +1077,11 @@ bool lar_solver::has_lower_bound(var_index var, constraint_index& ci, mpq& value
bool lar_solver::has_upper_bound(var_index var, constraint_index& ci, mpq& value, bool& is_strict) {
if (var >= m_vars_to_ul_pairs.size()) {
if (var >= m_columns_to_ul_pairs.size()) {
// TBD: bounds on terms could also be used, caller may have to track these.
return false;
}
const ul_pair & ul = m_vars_to_ul_pairs[var];
const ul_pair & ul = m_columns_to_ul_pairs[var];
ci = ul.upper_bound_witness();
if (ci != static_cast<constraint_index>(-1)) {
auto& p = m_mpq_lar_core_solver.m_r_upper_bounds()[var];
@ -1103,7 +1122,7 @@ void lar_solver::get_infeasibility_explanation_for_inf_sign(
unsigned j = it.second;
int adj_sign = coeff.is_pos() ? inf_sign : -inf_sign;
const ul_pair & ul = m_vars_to_ul_pairs[j];
const ul_pair & ul = m_columns_to_ul_pairs[j];
constraint_index bound_constr_i = adj_sign < 0 ? ul.upper_bound_witness() : ul.low_bound_witness();
lp_assert(bound_constr_i < m_constraints.size());
@ -1113,7 +1132,7 @@ void lar_solver::get_infeasibility_explanation_for_inf_sign(
void lar_solver::get_model(std::unordered_map<var_index, mpq> & variable_values) const {
mpq delta = mpq(1, 2); // start from 0.5 to have less clashes
lp_assert(m_status == OPTIMAL);
lp_assert(m_status == lp_status::OPTIMAL);
unsigned i;
do {
@ -1188,6 +1207,13 @@ void lar_solver::print_term(lar_term const& term, std::ostream & out) const {
print_linear_combination_of_column_indices(term.coeffs_as_vector(), out);
}
void lar_solver::print_term_as_indices(lar_term const& term, std::ostream & out) const {
if (!numeric_traits<mpq>::is_zero(term.m_v)) {
out << term.m_v << " + ";
}
print_linear_combination_of_column_indices_only(term.coeffs_as_vector(), out);
}
mpq lar_solver::get_left_side_val(const lar_base_constraint & cns, const std::unordered_map<var_index, mpq> & var_map) const {
mpq ret = cns.get_free_coeff_of_left_side();
for (auto & it : cns.get_left_side_coefficients()) {
@ -1234,7 +1260,7 @@ void lar_solver::pop() {
}
bool lar_solver::column_represents_row_in_tableau(unsigned j) {
return m_vars_to_ul_pairs()[j].m_i != static_cast<row_index>(-1);
return m_columns_to_ul_pairs()[j].m_i != static_cast<row_index>(-1);
}
void lar_solver::make_sure_that_the_bottom_right_elem_not_zero_in_tableau(unsigned i, unsigned j) {
@ -1353,7 +1379,7 @@ void lar_solver::pop_tableau() {
void lar_solver::clean_inf_set_of_r_solver_after_pop() {
vector<unsigned> became_feas;
clean_large_elements_after_pop(A_r().column_count(), m_mpq_lar_core_solver.m_r_solver.m_inf_set);
clean_popped_elements(A_r().column_count(), m_mpq_lar_core_solver.m_r_solver.m_inf_set);
std::unordered_set<unsigned> basic_columns_with_changed_cost;
auto inf_index_copy = m_mpq_lar_core_solver.m_r_solver.m_inf_set.m_index;
for (auto j: inf_index_copy) {
@ -1458,11 +1484,14 @@ var_index lar_solver::add_var(unsigned ext_j, bool is_int) {
if (it != m_ext_vars_to_columns.end()) {
return it->second.ext_j();
}
lp_assert(m_vars_to_ul_pairs.size() == A_r().column_count());
lp_assert(m_columns_to_ul_pairs.size() == A_r().column_count());
i = A_r().column_count();
m_vars_to_ul_pairs.push_back(ul_pair(static_cast<unsigned>(-1)));
m_columns_to_ul_pairs.push_back(ul_pair(static_cast<unsigned>(-1)));
add_non_basic_var_to_core_fields(ext_j, is_int);
lp_assert(sizes_are_correct());
if (is_int) {
m_mpq_lar_core_solver.m_r_solver.set_tracker_of_x(& m_tracker_of_x_change);
}
return i;
}
@ -1564,7 +1593,7 @@ void lar_solver::add_row_from_term_no_constraint(const lar_term * term, unsigned
// j will be a new variable
unsigned j = A_r().column_count();
ul_pair ul(j);
m_vars_to_ul_pairs.push_back(ul);
m_columns_to_ul_pairs.push_back(ul);
add_basic_var_to_core_fields();
if (use_tableau()) {
auto it = iterator_on_term_with_basis_var(*term, j);
@ -1598,6 +1627,7 @@ bool lar_solver::bound_is_integer_if_needed(unsigned j, const mpq & right_side)
}
constraint_index lar_solver::add_var_bound(var_index j, lconstraint_kind kind, const mpq & right_side) {
TRACE("lar_solver", tout << "j = " << j << std::endl;);
constraint_index ci = m_constraints.size();
if (!is_term(j)) { // j is a var
lp_assert(bound_is_integer_if_needed(j, right_side));
@ -1672,7 +1702,7 @@ void lar_solver::add_constraint_from_term_and_create_new_column_row(unsigned ter
void lar_solver::decide_on_strategy_and_adjust_initial_state() {
lp_assert(strategy_is_undecided());
if (m_vars_to_ul_pairs.size() > m_settings.column_number_threshold_for_using_lu_in_lar_solver) {
if (m_columns_to_ul_pairs.size() > m_settings.column_number_threshold_for_using_lu_in_lar_solver) {
m_settings.simplex_strategy() = simplex_strategy_enum::lu;
}
else {
@ -1816,7 +1846,7 @@ void lar_solver::update_upper_bound_column_type_and_bound(var_index j, lconstrai
set_low_bound_witness(j, ci);
m_columns_with_changed_bound.insert(j);
if (low > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
m_status = INFEASIBLE;
m_status = lp_status::INFEASIBLE;
m_infeasible_column_index = j;
}
else {
@ -1828,7 +1858,7 @@ void lar_solver::update_upper_bound_column_type_and_bound(var_index j, lconstrai
{
auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
m_status = INFEASIBLE;
m_status = lp_status::INFEASIBLE;
set_low_bound_witness(j, ci);
m_infeasible_column_index = j;
}
@ -1850,7 +1880,7 @@ void lar_solver::update_upper_bound_column_type_and_bound(var_index j, lconstrai
}
void lar_solver::update_boxed_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
lp_assert(m_status == INFEASIBLE || (m_mpq_lar_core_solver.m_column_types()[j] == column_type::boxed && m_mpq_lar_core_solver.m_r_low_bounds()[j] < m_mpq_lar_core_solver.m_r_upper_bounds()[j]));
lp_assert(m_status == lp_status::INFEASIBLE || (m_mpq_lar_core_solver.m_column_types()[j] == column_type::boxed && m_mpq_lar_core_solver.m_r_low_bounds()[j] < m_mpq_lar_core_solver.m_r_upper_bounds()[j]));
mpq y_of_bound(0);
switch (kind) {
case LT:
@ -1865,7 +1895,7 @@ void lar_solver::update_boxed_column_type_and_bound(var_index j, lconstraint_kin
}
if (up < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
m_status = INFEASIBLE;
m_status = lp_status::INFEASIBLE;
lp_assert(false);
m_infeasible_column_index = j;
}
@ -1886,7 +1916,7 @@ void lar_solver::update_boxed_column_type_and_bound(var_index j, lconstraint_kin
set_low_bound_witness(j, ci);
}
if (low > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
m_status = INFEASIBLE;
m_status = lp_status::INFEASIBLE;
m_infeasible_column_index = j;
}
else if (low == m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
@ -1898,12 +1928,12 @@ void lar_solver::update_boxed_column_type_and_bound(var_index j, lconstraint_kin
{
auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
m_status = INFEASIBLE;
m_status = lp_status::INFEASIBLE;
m_infeasible_column_index = j;
set_upper_bound_witness(j, ci);
}
else if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
m_status = INFEASIBLE;
m_status = lp_status::INFEASIBLE;
m_infeasible_column_index = j;
set_low_bound_witness(j, ci);
}
@ -1937,7 +1967,7 @@ void lar_solver::update_low_bound_column_type_and_bound(var_index j, lconstraint
m_columns_with_changed_bound.insert(j);
if (up < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
m_status = INFEASIBLE;
m_status = lp_status::INFEASIBLE;
m_infeasible_column_index = j;
}
else {
@ -1961,7 +1991,7 @@ void lar_solver::update_low_bound_column_type_and_bound(var_index j, lconstraint
{
auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
m_status = INFEASIBLE;
m_status = lp_status::INFEASIBLE;
m_infeasible_column_index = j;
set_upper_bound_witness(j, ci);
}
@ -1982,15 +2012,15 @@ void lar_solver::update_low_bound_column_type_and_bound(var_index j, lconstraint
}
void lar_solver::update_fixed_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
lp_assert(m_status == INFEASIBLE || (m_mpq_lar_core_solver.m_column_types()[j] == column_type::fixed && m_mpq_lar_core_solver.m_r_low_bounds()[j] == m_mpq_lar_core_solver.m_r_upper_bounds()[j]));
lp_assert(m_status == INFEASIBLE || (m_mpq_lar_core_solver.m_r_low_bounds()[j].y.is_zero() && m_mpq_lar_core_solver.m_r_upper_bounds()[j].y.is_zero()));
lp_assert(m_status == lp_status::INFEASIBLE || (m_mpq_lar_core_solver.m_column_types()[j] == column_type::fixed && m_mpq_lar_core_solver.m_r_low_bounds()[j] == m_mpq_lar_core_solver.m_r_upper_bounds()[j]));
lp_assert(m_status == lp_status::INFEASIBLE || (m_mpq_lar_core_solver.m_r_low_bounds()[j].y.is_zero() && m_mpq_lar_core_solver.m_r_upper_bounds()[j].y.is_zero()));
auto v = numeric_pair<mpq>(right_side, mpq(0));
mpq y_of_bound(0);
switch (kind) {
case LT:
if (v <= m_mpq_lar_core_solver.m_r_low_bounds[j]) {
m_status = INFEASIBLE;
m_status = lp_status::INFEASIBLE;
m_infeasible_column_index = j;
set_upper_bound_witness(j, ci);
}
@ -1998,7 +2028,7 @@ void lar_solver::update_fixed_column_type_and_bound(var_index j, lconstraint_kin
case LE:
{
if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
m_status = INFEASIBLE;
m_status = lp_status::INFEASIBLE;
m_infeasible_column_index = j;
set_upper_bound_witness(j, ci);
}
@ -2007,7 +2037,7 @@ void lar_solver::update_fixed_column_type_and_bound(var_index j, lconstraint_kin
case GT:
{
if (v >= m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
m_status = INFEASIBLE;
m_status = lp_status::INFEASIBLE;
m_infeasible_column_index = j;
set_low_bound_witness(j, ci);
}
@ -2016,7 +2046,7 @@ void lar_solver::update_fixed_column_type_and_bound(var_index j, lconstraint_kin
case GE:
{
if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
m_status = INFEASIBLE;
m_status = lp_status::INFEASIBLE;
m_infeasible_column_index = j;
set_low_bound_witness(j, ci);
}
@ -2025,12 +2055,12 @@ void lar_solver::update_fixed_column_type_and_bound(var_index j, lconstraint_kin
case EQ:
{
if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
m_status = INFEASIBLE;
m_status = lp_status::INFEASIBLE;
m_infeasible_column_index = j;
set_upper_bound_witness(j, ci);
}
else if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
m_status = INFEASIBLE;
m_status = lp_status::INFEASIBLE;
m_infeasible_column_index = j;
set_low_bound_witness(j, ci);
}

117
src/util/lp/lar_solver.h
View file

@ -61,11 +61,11 @@ class lar_solver : public column_namer {
};
//////////////////// fields //////////////////////////
lp_settings m_settings;
stacked_value<lp_status> m_status;
lp_status m_status;
stacked_value<simplex_strategy_enum> m_simplex_strategy;
std::unordered_map<unsigned, ext_var_info> m_ext_vars_to_columns;
vector<unsigned> m_columns_to_ext_vars_or_term_indices;
stacked_vector<ul_pair> m_vars_to_ul_pairs;
stacked_vector<ul_pair> m_columns_to_ul_pairs;
vector<lar_base_constraint*> m_constraints;
stacked_value<unsigned> m_constraint_count;
// the set of column indices j such that bounds have changed for j
@ -83,7 +83,8 @@ public:
lar_core_solver m_mpq_lar_core_solver;
unsigned constraint_count() const;
const lar_base_constraint& get_constraint(unsigned ci) const;
std::function<void (unsigned, const impq&)> m_tracker_of_x_change;
int_set m_inf_int_set;
////////////////// methods ////////////////////////////////
static_matrix<mpq, numeric_pair<mpq>> & A_r() { return m_mpq_lar_core_solver.m_r_A;}
static_matrix<mpq, numeric_pair<mpq>> const & A_r() const { return m_mpq_lar_core_solver.m_r_A;}
@ -131,14 +132,11 @@ public:
bool use_lu() const { return m_settings.simplex_strategy() == simplex_strategy_enum::lu; }
bool sizes_are_correct() const {
SASSERT(strategy_is_undecided() || !m_mpq_lar_core_solver.need_to_presolve_with_double_solver() || A_r().column_count() == A_d().column_count());
SASSERT(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_column_types.size());
SASSERT(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_costs.size());
SASSERT(A_r().column_count() == m_mpq_lar_core_solver.m_r_x.size());
return true;
}
void clear();
lar_solver();
void set_propagate_bounds_on_pivoted_rows_mode(bool v);
virtual ~lar_solver();
void print_implied_bound(const implied_bound& be, std::ostream & out) const {
out << "implied bound\n";
@ -462,7 +460,7 @@ public:
vector<unsigned> get_list_of_all_var_indices() const;
void push();
static void clean_large_elements_after_pop(unsigned n, int_set& set);
static void clean_popped_elements(unsigned n, int_set& set);
static void shrink_inf_set_after_pop(unsigned n, int_set & set);
@ -646,20 +644,8 @@ public:
void set_low_bound_witness(var_index j, constraint_index ci);
void substitute_terms(const mpq & mult,
const vector<std::pair<mpq, var_index>>& left_side_with_terms,
vector<std::pair<mpq, var_index>> &left_side, mpq & right_side) const {
for (auto & t : left_side_with_terms) {
if (t.second < m_terms_start_index) {
SASSERT(t.second < A_r().column_count());
left_side.push_back(std::pair<mpq, var_index>(mult * t.first, t.second));
} else {
const lar_term & term = * m_terms[adjust_term_index(t.second)];
substitute_terms(mult * t.first, left_side_with_terms, left_side, right_side);
right_side -= mult * term.m_v;
}
}
}
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, mpq & free_coeff) const;
void detect_rows_of_bound_change_column_for_nbasic_column(unsigned j) {
@ -1020,7 +1006,7 @@ public:
bool the_relations_are_of_same_type(const vector<std::pair<mpq, unsigned>> & evidence, lconstraint_kind & the_kind_of_sum) const;
static void register_in_map(std::unordered_map<var_index, mpq> & coeffs, const lar_base_constraint & cn, const mpq & a);
static void register_one_coeff_in_map(std::unordered_map<var_index, mpq> & coeffs, const mpq & a, unsigned j);
static void register_monoid_in_map(std::unordered_map<var_index, mpq> & coeffs, const mpq & a, unsigned j);
bool the_left_sides_sum_to_zero(const vector<std::pair<mpq, unsigned>> & evidence) const;
@ -1227,13 +1213,7 @@ public:
void print_constraints(std::ostream& out) const ;
void print_left_side_of_constraint(const lar_base_constraint * c, std::ostream & out) const {
print_linear_combination_of_column_indices(c->get_left_side_coefficients(), out);
mpq free_coeff = c->get_free_coeff_of_left_side();
if (!is_zero(free_coeff))
out << " + " << free_coeff;
}
void print_terms(std::ostream& out) const;
void print_left_side_of_constraint(const lar_base_constraint * c, std::ostream & out) const;
@ -1248,6 +1228,8 @@ public:
return ret;
}
void print_term_as_indices(lar_term const& term, std::ostream & out) const;
mpq get_left_side_val(const lar_base_constraint & cns, const std::unordered_map<var_index, mpq> & var_map) const;
void fill_var_set_for_random_update(unsigned sz, var_index const * vars, vector<unsigned>& column_list) {
@ -1283,22 +1265,12 @@ public:
return m_vars_to_ul_pairs()[j].m_i != static_cast<row_index>(-1);
}
void make_sure_that_the_bottom_right_elem_not_zero_in_tableau(unsigned i, unsigned j) {
// i, j - is the indices of the bottom-right element of the tableau
SASSERT(A_r().row_count() == i + 1 && A_r().column_count() == j + 1);
auto & last_column = A_r().m_columns[j];
int non_zero_column_cell_index = -1;
for (unsigned k = last_column.size(); k-- > 0;){
auto & cc = last_column[k];
if (cc.m_i == i)
return;
non_zero_column_cell_index = k;
}
bool column_is_real(unsigned j) const {
return !column_is_int(j);
}
bool model_is_int_feasible() const;
SASSERT(non_zero_column_cell_index != -1);
SASSERT(static_cast<unsigned>(non_zero_column_cell_index) != i);
m_mpq_lar_core_solver.m_r_solver.transpose_rows_tableau(last_column[non_zero_column_cell_index].m_i, i);
}
void remove_last_row_and_column_from_tableau(unsigned j) {
SASSERT(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_costs.size());
@ -1383,19 +1355,10 @@ public:
SASSERT(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_costs.size());
}
void pop_tableau() {
SASSERT(m_mpq_lar_core_solver.m_r_solver.m_costs.size() == A_r().column_count());
SASSERT(m_mpq_lar_core_solver.m_r_solver.m_basis.size() == A_r().row_count());
SASSERT(m_mpq_lar_core_solver.m_r_solver.basis_heading_is_correct());
// We remove last variables starting from m_column_names.size() to m_vec_of_canonic_left_sides.size().
// At this moment m_column_names is already popped
for (unsigned j = A_r().column_count(); j-- > m_columns_to_ext_vars_or_term_indices.size();)
remove_column_from_tableau(j);
SASSERT(m_mpq_lar_core_solver.m_r_solver.m_costs.size() == A_r().column_count());
SASSERT(m_mpq_lar_core_solver.m_r_solver.m_basis.size() == A_r().row_count());
SASSERT(m_mpq_lar_core_solver.m_r_solver.basis_heading_is_correct());
void get_bound_constraint_witnesses_for_column(unsigned j, constraint_index & lc, constraint_index & uc) const {
const ul_pair & ul = m_columns_to_ul_pairs[j];
lc = ul.low_bound_witness();
uc = ul.upper_bound_witness();
}
@ -1448,5 +1411,37 @@ public:
quick_xplain::run(explanation, *this);
SASSERT(this->explanation_is_correct(explanation));
}
bool bound_is_integer_if_needed(unsigned j, const mpq & right_side) const;
linear_combination_iterator<mpq> * get_iterator_on_row(unsigned i) {
return m_mpq_lar_core_solver.m_r_solver.get_iterator_on_row(i);
}
unsigned get_base_column_in_row(unsigned row_index) const {
return m_mpq_lar_core_solver.m_r_solver.get_base_column_in_row(row_index);
}
constraint_index get_column_upper_bound_witness(unsigned j) const {
return m_columns_to_ul_pairs()[j].upper_bound_witness();
}
constraint_index get_column_low_bound_witness(unsigned j) const {
return m_columns_to_ul_pairs()[j].low_bound_witness();
}
void subs_terms_for_debugging(lar_term& t) {
vector<std::pair<mpq, unsigned>> pol;
for (const auto & m : t.m_coeffs) {
pol.push_back(std::make_pair(m.second, adjust_column_index_to_term_index(m.first)));
}
mpq v = t.m_v;
vector<std::pair<mpq, unsigned>> pol_after_subs;
substitute_terms_in_linear_expression(pol, pol_after_subs, v);
t.clear();
t = lar_term(pol_after_subs, v);
}
bool has_int_var() const;
};
}

17
src/util/lp/lar_term.h
View file

@ -25,7 +25,7 @@ struct lar_term {
std::unordered_map<unsigned, mpq> m_coeffs;
mpq m_v;
lar_term() {}
void add_to_map(unsigned j, const mpq& c) {
void add_monoid(const mpq& c, unsigned j) {
auto it = m_coeffs.find(j);
if (it == m_coeffs.end()) {
m_coeffs.emplace(j, c);
@ -49,7 +49,7 @@ struct lar_term {
lar_term(const vector<std::pair<mpq, unsigned>>& coeffs,
const mpq & v) : m_v(v) {
for (const auto & p : coeffs) {
add_to_map(p.second, p.first);
add_monoid(p.first, p.second);
}
}
bool operator==(const lar_term & a) const { return false; } // take care not to create identical terms
@ -71,7 +71,7 @@ struct lar_term {
if (it == m_coeffs.end()) return;
const mpq & b = it->second;
for (unsigned it_j :li.m_index) {
add_to_map(it_j, - b * li.m_data[it_j]);
add_monoid(- b * li.m_data[it_j], it_j);
}
m_coeffs.erase(it);
}
@ -79,5 +79,16 @@ struct lar_term {
bool contains(unsigned j) const {
return m_coeffs.find(j) != m_coeffs.end();
}
void negate() {
for (auto & t : m_coeffs)
t.second.neg();
}
void clear() {
m_coeffs.clear();
m_v = zero_of_type<mpq>();
}
};
}

26
src/util/lp/lp_core_solver_base.h
View file

@ -41,7 +41,14 @@ class lp_core_solver_base {
private:
lp_status m_status;
public:
bool current_x_is_feasible() const { return m_inf_set.size() == 0; }
bool current_x_is_feasible() const {
TRACE("feas",
if (m_inf_set.size()) {
tout << "column " << m_inf_set.m_index[0] << " is infeasible" << std::endl;
}
);
return m_inf_set.size() == 0;
}
bool current_x_is_infeasible() const { return m_inf_set.size() != 0; }
int_set m_inf_set;
bool m_using_infeas_costs;
@ -77,6 +84,12 @@ public:
bool m_tracing_basis_changes;
int_set* m_pivoted_rows;
bool m_look_for_feasible_solution_only;
std::function<void (unsigned, const X &)> * m_tracker_of_x_change;
void set_tracker_of_x(std::function<void (unsigned, const X&)>* tracker) {
m_tracker_of_x_change = tracker;
}
void start_tracing_basis_changes() {
m_trace_of_basis_change_vector.resize(0);
m_tracing_basis_changes = true;
@ -678,16 +691,20 @@ public:
void update_column_in_inf_set(unsigned j) {
if (column_is_feasible(j)) {
m_inf_set.erase(j);
remove_column_from_inf_set(j);
} else {
m_inf_set.insert(j);
insert_column_into_inf_set(j);
}
}
void insert_column_into_inf_set(unsigned j) {
if (m_tracker_of_x_change != nullptr)
(*m_tracker_of_x_change)(j, m_x[j]);
m_inf_set.insert(j);
lp_assert(!column_is_feasible(j));
}
void remove_column_from_inf_set(unsigned j) {
if (m_tracker_of_x_change != nullptr)
(*m_tracker_of_x_change)(j, m_x[j]);
m_inf_set.erase(j);
lp_assert(column_is_feasible(j));
}
@ -715,5 +732,8 @@ public:
}
void calculate_pivot_row(unsigned i);
unsigned get_base_column_in_row(unsigned row_index) const {
return m_basis[row_index];
}
};
}

16
src/util/lp/lp_core_solver_base.hpp
View file

@ -39,7 +39,7 @@ lp_core_solver_base(static_matrix<T, X> & A,
const vector<X> & upper_bound_values):
m_total_iterations(0),
m_iters_with_no_cost_growing(0),
m_status(FEASIBLE),
m_status(lp_status::FEASIBLE),
m_inf_set(A.column_count()),
m_using_infeas_costs(false),
m_pivot_row_of_B_1(A.row_count()),
@ -67,7 +67,8 @@ lp_core_solver_base(static_matrix<T, X> & A,
m_steepest_edge_coefficients(A.column_count()),
m_tracing_basis_changes(false),
m_pivoted_rows(nullptr),
m_look_for_feasible_solution_only(false) {
m_look_for_feasible_solution_only(false),
m_tracker_of_x_change(nullptr) {
lp_assert(bounds_for_boxed_are_set_correctly());
init();
init_basis_heading_and_non_basic_columns_vector();
@ -549,7 +550,7 @@ update_basis_and_x(int entering, int leaving, X const & tt) {
if (!find_x_by_solving()) {
restore_x(entering, tt);
if(A_mult_x_is_off()) {
m_status = FLOATING_POINT_ERROR;
m_status = lp_status::FLOATING_POINT_ERROR;
m_iters_with_no_cost_growing++;
return false;
}
@ -559,7 +560,7 @@ update_basis_and_x(int entering, int leaving, X const & tt) {
if (m_factorization->get_status() != LU_status::OK) {
std::stringstream s;
// s << "failing refactor on off_result for entering = " << entering << ", leaving = " << leaving << " total_iterations = " << total_iterations();
m_status = FLOATING_POINT_ERROR;
m_status = lp_status::FLOATING_POINT_ERROR;
return false;
}
return false;
@ -581,19 +582,19 @@ update_basis_and_x(int entering, int leaving, X const & tt) {
init_lu();
if (m_factorization->get_status() != LU_status::OK) {
if (m_look_for_feasible_solution_only && !precise()) {
m_status = UNSTABLE;
m_status = lp_status::UNSTABLE;
delete m_factorization;
m_factorization = nullptr;
return false;
}
// LP_OUT(m_settings, "failing refactor for entering = " << entering << ", leaving = " << leaving << " total_iterations = " << total_iterations() << std::endl);
restore_x_and_refactor(entering, leaving, tt);
if (m_status == FLOATING_POINT_ERROR)
if (m_status == lp_status::FLOATING_POINT_ERROR)
return false;
lp_assert(!A_mult_x_is_off());
m_iters_with_no_cost_growing++;
// LP_OUT(m_settings, "rolled back after failing of init_factorization()" << std::endl);
m_status = UNSTABLE;
m_status = lp_status::UNSTABLE;
return false;
}
return true;
@ -975,7 +976,6 @@ template <typename T, typename X> void lp_core_solver_base<T, X>::pivot_fixed_v
unsigned basic_j = m_basis[i];
if (get_column_type(basic_j) != column_type::fixed) continue;
//todo run over the row here!!!!! call get_iterator_on_row();
T a;
unsigned j;
auto * it = get_iterator_on_row(i);

30
src/util/lp/lp_dual_core_solver.hpp
View file

@ -97,11 +97,11 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::start_with_ini
}
template <typename T, typename X> bool lp_dual_core_solver<T, X>::done() {
if (this->get_status() == OPTIMAL) {
if (this->get_status() == lp_status::OPTIMAL) {
return true;
}
if (this->total_iterations() > this->m_settings.max_total_number_of_iterations) { // debug !!!!
this->set_status(ITERATIONS_EXHAUSTED);
this->set_status(lp_status::ITERATIONS_EXHAUSTED);
return true;
}
return false; // todo, need to be more cases
@ -185,8 +185,8 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::pricing_loop(u
}
} while (i != initial_offset_in_rows && rows_left);
if (m_r == -1) {
if (this->get_status() != UNSTABLE) {
this->set_status(OPTIMAL);
if (this->get_status() != lp_status::UNSTABLE) {
this->set_status(lp_status::OPTIMAL);
}
} else {
m_p = this->m_basis[m_r];
@ -196,10 +196,10 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::pricing_loop(u
return;
}
// failure in advance_on_known_p
if (this->get_status() == FLOATING_POINT_ERROR) {
if (this->get_status() == lp_status::FLOATING_POINT_ERROR) {
return;
}
this->set_status(UNSTABLE);
this->set_status(lp_status::UNSTABLE);
m_forbidden_rows.insert(m_r);
}
}
@ -481,12 +481,12 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::revert_to_prev
this->change_basis_unconditionally(m_p, m_q);
init_factorization(this->m_factorization, this->m_A, this->m_basis, this->m_settings);
if (this->m_factorization->get_status() != LU_status::OK) {
this->set_status(FLOATING_POINT_ERROR); // complete failure
this->set_status(lp_status::FLOATING_POINT_ERROR); // complete failure
return;
}
recover_leaving();
if (!this->find_x_by_solving()) {
this->set_status(FLOATING_POINT_ERROR);
this->set_status(lp_status::FLOATING_POINT_ERROR);
return;
}
recalculate_xB_and_d();
@ -566,10 +566,10 @@ template <typename T, typename X> bool lp_dual_core_solver<T, X>::delta_keeps_th
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::set_status_to_tentative_dual_unbounded_or_dual_unbounded() {
if (this->get_status() == TENTATIVE_DUAL_UNBOUNDED) {
this->set_status(DUAL_UNBOUNDED);
if (this->get_status() == lp_status::TENTATIVE_DUAL_UNBOUNDED) {
this->set_status(lp_status::DUAL_UNBOUNDED);
} else {
this->set_status(TENTATIVE_DUAL_UNBOUNDED);
this->set_status(lp_status::TENTATIVE_DUAL_UNBOUNDED);
}
}
@ -675,7 +675,7 @@ template <typename T, typename X> bool lp_dual_core_solver<T, X>::ratio_test() {
set_status_to_tentative_dual_unbounded_or_dual_unbounded();
return false;
}
this->set_status(FEASIBLE);
this->set_status(lp_status::FEASIBLE);
find_q_and_tight_set();
if (!tight_breakpoinst_are_all_boxed()) break;
T del = m_delta - delta_lost_on_flips_of_tight_breakpoints() * initial_delta_sign;
@ -731,10 +731,10 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::update_xb_afte
template <typename T, typename X> void lp_dual_core_solver<T, X>::one_iteration() {
unsigned number_of_rows_to_try = get_number_of_rows_to_try_for_leaving();
unsigned offset_in_rows = this->m_settings.random_next() % this->m_m();
if (this->get_status() == TENTATIVE_DUAL_UNBOUNDED) {
if (this->get_status() == lp_status::TENTATIVE_DUAL_UNBOUNDED) {
number_of_rows_to_try = this->m_m();
} else {
this->set_status(FEASIBLE);
this->set_status(lp_status::FEASIBLE);
}
pricing_loop(number_of_rows_to_try, offset_in_rows);
lp_assert(problem_is_dual_feasible());
@ -751,7 +751,7 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::solve() { // s
return;
}
one_iteration();
} while (this->get_status() != FLOATING_POINT_ERROR && this->get_status() != DUAL_UNBOUNDED && this->get_status() != OPTIMAL &&
} while (this->get_status() != lp_status::FLOATING_POINT_ERROR && this->get_status() != lp_status::DUAL_UNBOUNDED && this->get_status() != lp_status::OPTIMAL &&
this->iters_with_no_cost_growing() <= this->m_settings.max_number_of_iterations_with_no_improvements
&& this->total_iterations() <= this->m_settings.max_total_number_of_iterations);
}

44
src/util/lp/lp_dual_simplex.hpp
View file

@ -22,23 +22,23 @@ namespace lp{
template <typename T, typename X> void lp_dual_simplex<T, X>::decide_on_status_after_stage1() {
switch (m_core_solver->get_status()) {
case OPTIMAL:
case lp_status::OPTIMAL:
if (this->m_settings.abs_val_is_smaller_than_artificial_tolerance(m_core_solver->get_cost())) {
this->m_status = FEASIBLE;
this->m_status = lp_status::FEASIBLE;
} else {
this->m_status = UNBOUNDED;
this->m_status = lp_status::UNBOUNDED;
}
break;
case DUAL_UNBOUNDED:
case lp_status::DUAL_UNBOUNDED:
lp_unreachable();
case ITERATIONS_EXHAUSTED:
this->m_status = ITERATIONS_EXHAUSTED;
case lp_status::ITERATIONS_EXHAUSTED:
this->m_status = lp_status::ITERATIONS_EXHAUSTED;
break;
case TIME_EXHAUSTED:
this->m_status = TIME_EXHAUSTED;
case lp_status::TIME_EXHAUSTED:
this->m_status = lp_status::TIME_EXHAUSTED;
break;
case FLOATING_POINT_ERROR:
this->m_status = FLOATING_POINT_ERROR;
case lp_status::FLOATING_POINT_ERROR:
this->m_status = lp_status::FLOATING_POINT_ERROR;
break;
default:
lp_unreachable();
@ -114,20 +114,20 @@ template <typename T, typename X> void lp_dual_simplex<T, X>::solve_for_stage2()
m_core_solver->solve_yB(m_core_solver->m_y);
m_core_solver->fill_reduced_costs_from_m_y_by_rows();
m_core_solver->start_with_initial_basis_and_make_it_dual_feasible();
m_core_solver->set_status(FEASIBLE);
m_core_solver->set_status(lp_status::FEASIBLE);
m_core_solver->solve();
switch (m_core_solver->get_status()) {
case OPTIMAL:
this->m_status = OPTIMAL;
case lp_status::OPTIMAL:
this->m_status = lp_status::OPTIMAL;
break;
case DUAL_UNBOUNDED:
this->m_status = INFEASIBLE;
case lp_status::DUAL_UNBOUNDED:
this->m_status = lp_status::INFEASIBLE;
break;
case TIME_EXHAUSTED:
this->m_status = TIME_EXHAUSTED;
case lp_status::TIME_EXHAUSTED:
this->m_status = lp_status::TIME_EXHAUSTED;
break;
case FLOATING_POINT_ERROR:
this->m_status = FLOATING_POINT_ERROR;
case lp_status::FLOATING_POINT_ERROR:
this->m_status = lp_status::FLOATING_POINT_ERROR;
break;
default:
lp_unreachable();
@ -166,7 +166,7 @@ template <typename T, typename X> void lp_dual_simplex<T, X>::stage1() {
m_core_solver->start_with_initial_basis_and_make_it_dual_feasible();
if (this->m_settings.abs_val_is_smaller_than_artificial_tolerance(m_core_solver->get_cost())) {
// skipping stage 1
m_core_solver->set_status(OPTIMAL);
m_core_solver->set_status(lp_status::OPTIMAL);
m_core_solver->set_total_iterations(0);
} else {
m_core_solver->solve();
@ -351,7 +351,7 @@ template <typename T, typename X> void lp_dual_simplex<T, X>::find_maximal_solut
this->flip_costs(); // do it for now, todo ( remove the flipping)
this->cleanup();
if (this->m_status == INFEASIBLE) {
if (this->m_status == lp_status::INFEASIBLE) {
return;
}
this->fill_matrix_A_and_init_right_side();
@ -361,7 +361,7 @@ template <typename T, typename X> void lp_dual_simplex<T, X>::find_maximal_solut
fill_first_stage_solver_fields();
copy_m_b_aside_and_set_it_to_zeros();
stage1();
if (this->m_status == FEASIBLE) {
if (this->m_status == lp_status::FEASIBLE) {
stage2();
}
}

12
src/util/lp/lp_primal_core_solver.h
View file

@ -484,7 +484,7 @@ public:
X new_val_for_leaving;
int leaving = find_leaving_tableau_rows(new_val_for_leaving);
if (leaving == -1) {
this->set_status(OPTIMAL);
this->set_status(lp_status::OPTIMAL);
return;
}
@ -500,14 +500,14 @@ public:
T a_ent;
int entering = find_beneficial_column_in_row_tableau_rows(this->m_basis_heading[leaving], a_ent);
if (entering == -1) {
this->set_status(INFEASIBLE);
this->set_status(lp_status::INFEASIBLE);
return;
}
X theta = (this->m_x[leaving] - new_val_for_leaving) / a_ent;
advance_on_entering_and_leaving_tableau_rows(entering, leaving, theta );
lp_assert(this->m_x[leaving] == new_val_for_leaving);
if (this->current_x_is_feasible())
this->set_status(OPTIMAL);
this->set_status(lp_status::OPTIMAL);
}
void fill_breakpoints_array(unsigned entering);
@ -523,7 +523,7 @@ public:
void decide_on_status_when_cannot_find_entering() {
lp_assert(!need_to_switch_costs());
this->set_status(this->current_x_is_feasible()? OPTIMAL: INFEASIBLE);
this->set_status(this->current_x_is_feasible()? lp_status::OPTIMAL: lp_status::INFEASIBLE);
}
// void limit_theta_on_basis_column_for_feas_case_m_neg(unsigned j, const T & m, X & theta) {
@ -794,7 +794,7 @@ public:
if (this->m_basis_heading[j] < 0)
continue;
if (!this->column_is_feasible(j))
this->m_inf_set.insert(j);
this->insert_column_into_inf_set(j);
}
}
@ -930,7 +930,7 @@ public:
} else {
m_converted_harris_eps = zero_of_type<T>();
}
this->set_status(UNKNOWN);
this->set_status(lp_status::UNKNOWN);
}
// constructor

78
src/util/lp/lp_primal_core_solver.hpp
View file

@ -713,14 +713,14 @@ template <typename T, typename X>void lp_primal_core_solver<T, X>::advance_on_en
int pivot_compare_result = this->pivots_in_column_and_row_are_different(entering, leaving);
if (!pivot_compare_result){;}
else if (pivot_compare_result == 2) { // the sign is changed, cannot continue
this->set_status(UNSTABLE);
this->set_status(lp_status::UNSTABLE);
this->iters_with_no_cost_growing()++;
return;
} else {
lp_assert(pivot_compare_result == 1);
this->init_lu();
if (this->m_factorization == nullptr || this->m_factorization->get_status() != LU_status::OK) {
this->set_status(UNSTABLE);
this->set_status(lp_status::UNSTABLE);
this->iters_with_no_cost_growing()++;
return;
}
@ -732,10 +732,10 @@ template <typename T, typename X>void lp_primal_core_solver<T, X>::advance_on_en
t = -t;
}
if (!this->update_basis_and_x(entering, leaving, t)) {
if (this->get_status() == FLOATING_POINT_ERROR)
if (this->get_status() == lp_status::FLOATING_POINT_ERROR)
return;
if (this->m_look_for_feasible_solution_only) {
this->set_status(FLOATING_POINT_ERROR);
this->set_status(lp_status::FLOATING_POINT_ERROR);
return;
}
init_reduced_costs();
@ -748,7 +748,7 @@ template <typename T, typename X>void lp_primal_core_solver<T, X>::advance_on_en
}
if (this->current_x_is_feasible()) {
this->set_status(FEASIBLE);
this->set_status(lp_status::FEASIBLE);
if (this->m_look_for_feasible_solution_only)
return;
}
@ -775,7 +775,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::advance_on_e
X t;
int leaving = find_leaving_and_t_precise(entering, t);
if (leaving == -1) {
this->set_status(UNBOUNDED);
this->set_status(lp_status::UNBOUNDED);
return;
}
advance_on_entering_and_leaving(entering, leaving, t);
@ -791,7 +791,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::advance_on_e
int refresh_result = refresh_reduced_cost_at_entering_and_check_that_it_is_off(entering);
if (refresh_result) {
if (this->m_look_for_feasible_solution_only) {
this->set_status(FLOATING_POINT_ERROR);
this->set_status(lp_status::FLOATING_POINT_ERROR);
return;
}
@ -814,19 +814,19 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::advance_on_e
// }
if (this->get_status() == UNSTABLE) {
this->set_status(FLOATING_POINT_ERROR);
if (this->get_status() == lp_status::UNSTABLE) {
this->set_status(lp_status::FLOATING_POINT_ERROR);
return;
}
init_infeasibility_costs();
this->set_status(UNSTABLE);
this->set_status(lp_status::UNSTABLE);
return;
}
if (this->get_status() == TENTATIVE_UNBOUNDED) {
this->set_status(UNBOUNDED);
if (this->get_status() == lp_status::TENTATIVE_UNBOUNDED) {
this->set_status(lp_status::UNBOUNDED);
} else {
this->set_status(TENTATIVE_UNBOUNDED);
this->set_status(lp_status::TENTATIVE_UNBOUNDED);
}
return;
}
@ -840,7 +840,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::push_forw
template <typename T, typename X> unsigned lp_primal_core_solver<T, X>::get_number_of_non_basic_column_to_try_for_enter() {
unsigned ret = static_cast<unsigned>(this->m_nbasis.size());
if (this->get_status() == TENTATIVE_UNBOUNDED)
if (this->get_status() == lp_status::TENTATIVE_UNBOUNDED)
return ret; // we really need to find entering with a large reduced cost
if (ret > 300) {
ret = (unsigned)(ret * this->m_settings.percent_of_entering_to_check / 100);
@ -867,12 +867,12 @@ template <typename T, typename X> unsigned lp_primal_core_solver<T, X>::solve()
init_run();
if (this->current_x_is_feasible() && this->m_look_for_feasible_solution_only) {
this->set_status(FEASIBLE);
this->set_status(lp_status::FEASIBLE);
return 0;
}
if ((!numeric_traits<T>::precise()) && this->A_mult_x_is_off()) {
this->set_status(FLOATING_POINT_ERROR);
this->set_status(lp_status::FLOATING_POINT_ERROR);
return 0;
}
do {
@ -882,8 +882,8 @@ template <typename T, typename X> unsigned lp_primal_core_solver<T, X>::solve()
one_iteration();
lp_assert(!this->m_using_infeas_costs || this->costs_on_nbasis_are_zeros());
switch (this->get_status()) {
case OPTIMAL: // double check that we are at optimum
case INFEASIBLE:
case lp_status::OPTIMAL: // double check that we are at optimum
case lp_status::INFEASIBLE:
if (this->m_look_for_feasible_solution_only && this->current_x_is_feasible())
break;
if (!numeric_traits<T>::precise()) {
@ -892,7 +892,7 @@ template <typename T, typename X> unsigned lp_primal_core_solver<T, X>::solve()
this->init_lu();
if (this->m_factorization->get_status() != LU_status::OK) {
this->set_status (FLOATING_POINT_ERROR);
this->set_status (lp_status::FLOATING_POINT_ERROR);
break;
}
init_reduced_costs();
@ -900,7 +900,7 @@ template <typename T, typename X> unsigned lp_primal_core_solver<T, X>::solve()
decide_on_status_when_cannot_find_entering();
break;
}
this->set_status(UNKNOWN);
this->set_status(lp_status::UNKNOWN);
} else { // precise case
if (this->m_look_for_feasible_solution_only) { // todo: keep the reduced costs correct all the time!
init_reduced_costs();
@ -908,31 +908,31 @@ template <typename T, typename X> unsigned lp_primal_core_solver<T, X>::solve()
decide_on_status_when_cannot_find_entering();
break;
}
this->set_status(UNKNOWN);
this->set_status(lp_status::UNKNOWN);
}
}
break;
case TENTATIVE_UNBOUNDED:
case lp_status::TENTATIVE_UNBOUNDED:
this->init_lu();
if (this->m_factorization->get_status() != LU_status::OK) {
this->set_status(FLOATING_POINT_ERROR);
this->set_status(lp_status::FLOATING_POINT_ERROR);
break;
}
init_reduced_costs();
break;
case UNBOUNDED:
case lp_status::UNBOUNDED:
if (this->current_x_is_infeasible()) {
init_reduced_costs();
this->set_status(UNKNOWN);
this->set_status(lp_status::UNKNOWN);
}
break;
case UNSTABLE:
case lp_status::UNSTABLE:
lp_assert(! (numeric_traits<T>::precise()));
this->init_lu();
if (this->m_factorization->get_status() != LU_status::OK) {
this->set_status(FLOATING_POINT_ERROR);
this->set_status(lp_status::FLOATING_POINT_ERROR);
break;
}
init_reduced_costs();
@ -941,13 +941,13 @@ template <typename T, typename X> unsigned lp_primal_core_solver<T, X>::solve()
default:
break; // do nothing
}
} while (this->get_status() != FLOATING_POINT_ERROR
} while (this->get_status() != lp_status::FLOATING_POINT_ERROR
&&
this->get_status() != UNBOUNDED
this->get_status() != lp_status::UNBOUNDED
&&
this->get_status() != OPTIMAL
this->get_status() != lp_status::OPTIMAL
&&
this->get_status() != INFEASIBLE
this->get_status() != lp_status::INFEASIBLE
&&
this->iters_with_no_cost_growing() <= this->m_settings.max_number_of_iterations_with_no_improvements
&&
@ -955,7 +955,7 @@ template <typename T, typename X> unsigned lp_primal_core_solver<T, X>::solve()
&&
!(this->current_x_is_feasible() && this->m_look_for_feasible_solution_only));
lp_assert(this->get_status() == FLOATING_POINT_ERROR
lp_assert(this->get_status() == lp_status::FLOATING_POINT_ERROR
||
this->current_x_is_feasible() == false
||
@ -1043,7 +1043,7 @@ template <typename T, typename X> T lp_primal_core_solver<T, X>::calculate_no
template <typename T, typename X> void lp_primal_core_solver<T, X>::find_feasible_solution() {
this->m_look_for_feasible_solution_only = true;
lp_assert(this->non_basic_columns_are_set_correctly());
this->set_status(UNKNOWN);
this->set_status(lp_status::UNKNOWN);
solve();
}
@ -1087,15 +1087,15 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::fill_breakpo
template <typename T, typename X> bool lp_primal_core_solver<T, X>::done() {
if (this->get_status() == OPTIMAL || this->get_status() == FLOATING_POINT_ERROR) return true;
if (this->get_status() == INFEASIBLE) {
if (this->get_status() == lp_status::OPTIMAL || this->get_status() == lp_status::FLOATING_POINT_ERROR) return true;
if (this->get_status() == lp_status::INFEASIBLE) {
return true;
}
if (this->m_iters_with_no_cost_growing >= this->m_settings.max_number_of_iterations_with_no_improvements) {
this->get_status() = ITERATIONS_EXHAUSTED; return true;
this->get_status() = lp_status::ITERATIONS_EXHAUSTED; return true;
}
if (this->total_iterations() >= this->m_settings.max_total_number_of_iterations) {
this->get_status() = ITERATIONS_EXHAUSTED; return true;
this->get_status() = lp_status::ITERATIONS_EXHAUSTED; return true;
}
return false;
}
@ -1212,9 +1212,9 @@ lp_primal_core_solver<T, X>::init_infeasibility_cost_for_column(unsigned j) {
}
if (numeric_traits<T>::is_zero(this->m_costs[j])) {
this->m_inf_set.erase(j);
this->remove_column_from_inf_set(j);
} else {
this->m_inf_set.insert(j);
this->insert_column_into_inf_set(j);
}
if (!this->m_settings.use_breakpoints_in_feasibility_search) {
this->m_costs[j] = - this->m_costs[j];

34
src/util/lp/lp_primal_core_solver_tableau.h
View file

@ -35,7 +35,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::advance_on_e
X t;
int leaving = find_leaving_and_t_tableau(entering, t);
if (leaving == -1) {
this->set_status(UNBOUNDED);
this->set_status(lp_status::UNBOUNDED);
return;
}
advance_on_entering_and_leaving_tableau(entering, leaving, t);
@ -100,12 +100,12 @@ template <typename T, typename X>
unsigned lp_primal_core_solver<T, X>::solve_with_tableau() {
init_run_tableau();
if (this->current_x_is_feasible() && this->m_look_for_feasible_solution_only) {
this->set_status(FEASIBLE);
this->set_status(lp_status::FEASIBLE);
return 0;
}
if ((!numeric_traits<T>::precise()) && this->A_mult_x_is_off()) {
this->set_status(FLOATING_POINT_ERROR);
this->set_status(lp_status::FLOATING_POINT_ERROR);
return 0;
}
do {
@ -117,8 +117,8 @@ unsigned lp_primal_core_solver<T, X>::solve_with_tableau() {
else
one_iteration_tableau();
switch (this->get_status()) {
case OPTIMAL: // double check that we are at optimum
case INFEASIBLE:
case lp_status::OPTIMAL: // double check that we are at optimum
case lp_status::INFEASIBLE:
if (this->m_look_for_feasible_solution_only && this->current_x_is_feasible())
break;
if (!numeric_traits<T>::precise()) {
@ -127,7 +127,7 @@ unsigned lp_primal_core_solver<T, X>::solve_with_tableau() {
this->init_lu();
if (this->m_factorization->get_status() != LU_status::OK) {
this->set_status(FLOATING_POINT_ERROR);
this->set_status(lp_status::FLOATING_POINT_ERROR);
break;
}
init_reduced_costs();
@ -135,7 +135,7 @@ unsigned lp_primal_core_solver<T, X>::solve_with_tableau() {
decide_on_status_when_cannot_find_entering();
break;
}
this->set_status(UNKNOWN);
this->set_status(lp_status::UNKNOWN);
} else { // precise case
if ((!this->infeasibility_costs_are_correct())) {
init_reduced_costs_tableau(); // forcing recalc
@ -143,31 +143,31 @@ unsigned lp_primal_core_solver<T, X>::solve_with_tableau() {
decide_on_status_when_cannot_find_entering();
break;
}
this->set_status(UNKNOWN);
this->set_status(lp_status::UNKNOWN);
}
}
break;
case TENTATIVE_UNBOUNDED:
case lp_status::TENTATIVE_UNBOUNDED:
this->init_lu();
if (this->m_factorization->get_status() != LU_status::OK) {
this->set_status(FLOATING_POINT_ERROR);
this->set_status(lp_status::FLOATING_POINT_ERROR);
break;
}
init_reduced_costs();
break;
case UNBOUNDED:
case lp_status::UNBOUNDED:
if (this->current_x_is_infeasible()) {
init_reduced_costs();
this->set_status(UNKNOWN);
this->set_status(lp_status::UNKNOWN);
}
break;
case UNSTABLE:
case lp_status::UNSTABLE:
lp_assert(! (numeric_traits<T>::precise()));
this->init_lu();
if (this->m_factorization->get_status() != LU_status::OK) {
this->set_status(FLOATING_POINT_ERROR);
this->set_status(lp_status::FLOATING_POINT_ERROR);
break;
}
init_reduced_costs();
@ -196,7 +196,7 @@ unsigned lp_primal_core_solver<T, X>::solve_with_tableau() {
this->set_status(CANCELLED);
}
lp_assert(this->get_status() == FLOATING_POINT_ERROR
lp_assert(this->get_status() == lp_status::FLOATING_POINT_ERROR
||
this->current_x_is_feasible() == false
||
@ -370,9 +370,9 @@ update_x_tableau(unsigned entering, const X& delta) {
this->m_x[j] -= delta * this->m_A.get_val(c);
update_inf_cost_for_column_tableau(j);
if (is_zero(this->m_costs[j]))
this->m_inf_set.erase(j);
this->remove_column_from_inf_set(j);
else
this->m_inf_set.insert(j);
this->insert_column_into_inf_set(j);
}
}
lp_assert(this->A_mult_x_is_off() == false);

2
src/util/lp/lp_primal_simplex.hpp
View file

@ -231,7 +231,7 @@ template <typename T, typename X> void lp_primal_simplex<T, X>::fill_A_x_and_bas
template <typename T, typename X> void lp_primal_simplex<T, X>::solve_with_total_inf() {
int total_vars = this->m_A->column_count() + this->row_count();
if (total_vars == 0) {
this->m_status = OPTIMAL;
this->m_status = lp_status::OPTIMAL;
return;
}
m_low_bounds.clear();

2
src/util/lp/lp_settings.h
View file

@ -51,7 +51,7 @@ enum class simplex_strategy_enum {
std::string column_type_to_string(column_type t);
enum lp_status {
enum class lp_status {
UNKNOWN,
INFEASIBLE,
TENTATIVE_UNBOUNDED,

24
src/util/lp/lp_settings.hpp
View file

@ -36,18 +36,18 @@ std::string column_type_to_string(column_type t) {
const char* lp_status_to_string(lp_status status) {
switch (status) {
case UNKNOWN: return "UNKNOWN";
case INFEASIBLE: return "INFEASIBLE";
case UNBOUNDED: return "UNBOUNDED";
case TENTATIVE_DUAL_UNBOUNDED: return "TENTATIVE_DUAL_UNBOUNDED";
case DUAL_UNBOUNDED: return "DUAL_UNBOUNDED";
case OPTIMAL: return "OPTIMAL";
case FEASIBLE: return "FEASIBLE";
case FLOATING_POINT_ERROR: return "FLOATING_POINT_ERROR";
case TIME_EXHAUSTED: return "TIME_EXHAUSTED";
case ITERATIONS_EXHAUSTED: return "ITERATIONS_EXHAUSTED";
case EMPTY: return "EMPTY";
case UNSTABLE: return "UNSTABLE";
case lp_status::UNKNOWN: return "UNKNOWN";
case lp_status::INFEASIBLE: return "INFEASIBLE";
case lp_status::UNBOUNDED: return "UNBOUNDED";
case lp_status::TENTATIVE_DUAL_UNBOUNDED: return "TENTATIVE_DUAL_UNBOUNDED";
case lp_status::DUAL_UNBOUNDED: return "DUAL_UNBOUNDED";
case lp_status::OPTIMAL: return "OPTIMAL";
case lp_status::FEASIBLE: return "FEASIBLE";
case lp_status::FLOATING_POINT_ERROR: return "FLOATING_POINT_ERROR";
case lp_status::TIME_EXHAUSTED: return "TIME_EXHAUSTED";
case lp_status::ITERATIONS_EXHAUSTED: return "ITERATIONS_EXHAUSTED";
case lp_status::EMPTY: return "EMPTY";
case lp_status::UNSTABLE: return "UNSTABLE";
default:
lp_unreachable();
}

12
src/util/lp/lp_solver.hpp
View file

@ -238,7 +238,7 @@ template <typename T, typename X> bool lp_solver<T, X>::row_e_is_obsolete(std
T rs = m_constraints[row_index].m_rs;
if (row_is_zero(row)) {
if (!is_zero(rs))
m_status = INFEASIBLE;
m_status = lp_status::INFEASIBLE;
return true;
}
@ -248,7 +248,7 @@ template <typename T, typename X> bool lp_solver<T, X>::row_e_is_obsolete(std
T diff = low_bound - rs;
if (!val_is_smaller_than_eps(diff, m_settings.refactor_tolerance)){
// low_bound > rs + m_settings.refactor_epsilon
m_status = INFEASIBLE;
m_status = lp_status::INFEASIBLE;
return true;
}
if (val_is_smaller_than_eps(-diff, m_settings.refactor_tolerance)){
@ -263,7 +263,7 @@ template <typename T, typename X> bool lp_solver<T, X>::row_e_is_obsolete(std
T diff = rs - upper_bound;
if (!val_is_smaller_than_eps(diff, m_settings.refactor_tolerance)) {
// upper_bound < rs - m_settings.refactor_tolerance
m_status = INFEASIBLE;
m_status = lp_status::INFEASIBLE;
return true;
}
if (val_is_smaller_than_eps(-diff, m_settings.refactor_tolerance)){
@ -279,7 +279,7 @@ template <typename T, typename X> bool lp_solver<T, X>::row_ge_is_obsolete(std:
T rs = m_constraints[row_index].m_rs;
if (row_is_zero(row)) {
if (rs > zero_of_type<X>())
m_status = INFEASIBLE;
m_status = lp_status::INFEASIBLE;
return true;
}
@ -288,7 +288,7 @@ template <typename T, typename X> bool lp_solver<T, X>::row_ge_is_obsolete(std:
T diff = rs - upper_bound;
if (!val_is_smaller_than_eps(diff, m_settings.refactor_tolerance)) {
// upper_bound < rs - m_settings.refactor_tolerance
m_status = INFEASIBLE;
m_status = lp_status::INFEASIBLE;
return true;
}
if (val_is_smaller_than_eps(-diff, m_settings.refactor_tolerance)){
@ -305,7 +305,7 @@ template <typename T, typename X> bool lp_solver<T, X>::row_le_is_obsolete(std::
T rs = m_constraints[row_index].m_rs;
if (row_is_zero(row)) {
if (rs < zero_of_type<X>())
m_status = INFEASIBLE;
m_status = lp_status::INFEASIBLE;
return true;
}

6
src/util/lp/quick_xplain.cpp
View file

@ -44,7 +44,7 @@ void quick_xplain::copy_constraint_and_add_constraint_vars(const lar_constraint&
bool quick_xplain::infeasible() {
m_qsol.solve();
return m_qsol.get_status() == INFEASIBLE;
return m_qsol.get_status() == lp_status::INFEASIBLE;
}
// u - unexplored constraints
@ -115,7 +115,7 @@ bool quick_xplain::is_feasible(const vector<unsigned> & x, unsigned k) const {
l.add_constraint(ls, c.m_kind, c.m_right_side);
}
l.solve();
return l.get_status() != INFEASIBLE;
return l.get_status() != lp_status::INFEASIBLE;
}
bool quick_xplain::x_is_minimal() const {
@ -142,7 +142,7 @@ void quick_xplain::solve() {
for (unsigned i : m_x)
add_constraint_to_qsol(i);
m_qsol.solve();
lp_assert(m_qsol.get_status() == INFEASIBLE);
lp_assert(m_qsol.get_status() == lp_status::INFEASIBLE);
m_qsol.get_infeasibility_explanation(m_explanation);
lp_assert(m_qsol.explanation_is_correct(m_explanation));
lp_assert(x_is_minimal());