diff --git a/src/shell/lp_frontend.cpp b/src/shell/lp_frontend.cpp index 3e68f237b..9145f8adf 100644 --- a/src/shell/lp_frontend.cpp +++ b/src/shell/lp_frontend.cpp @@ -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(); } diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp index 6460af4b9..894844008 100644 --- a/src/smt/theory_lra.cpp +++ b/src/smt/theory_lra.cpp @@ -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(); diff --git a/src/test/lp/lp.cpp b/src/test/lp/lp.cpp index 7fccebbbc..adcf9b969 100644 --- a/src/test/lp/lp.cpp +++ b/src/test/lp/lp.cpp @@ -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 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> 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 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 model; - lp_assert(solver.get_status() == INFEASIBLE); + lp_assert(solver.get_status() == lp_status::INFEASIBLE); } void test_bound_propagation_one_small_sample1() { /* diff --git a/src/util/lp/int_solver.cpp b/src/util/lp/int_solver.cpp index b152d8bcf..3345a6389 100644 --- a/src/util/lp/int_solver.cpp +++ b/src/util/lp/int_solver.cpp @@ -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 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::const_iterator it = r.begin_entries(); - typename vector::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> & 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 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> & 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(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(); + 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::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; +} } diff --git a/src/util/lp/int_solver.h b/src/util/lp/int_solver.h index af331eacd..9820e484c 100644 --- a/src/util/lp/int_solver.h +++ b/src/util/lp/int_solver.h @@ -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 @@ -23,6 +22,9 @@ enum class lia_move { struct explanation { vector> 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 m_old_values_data; - int_set m_inf_int_set; unsigned m_branch_cut_counter; + linear_combination_iterator* 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 * 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> & pol, lar_term & t, mpq &k); + void gomory_cut_adjust_t_and_k_for_size_gt_1(vector> & pol, lar_term & t, mpq &k, unsigned num_ints, mpq &lcm_den); }; } diff --git a/src/util/lp/lar_core_solver.hpp b/src/util/lp/lar_core_solver.hpp index 95a31c072..4f1a8cfc2 100644 --- a/src/util/lp/lar_core_solver.hpp +++ b/src/util/lp/lar_core_solver.hpp @@ -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 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()); diff --git a/src/util/lp/lar_solver.cpp b/src/util/lp/lar_solver.cpp index 0b5508ec9..712f8a001 100644 --- a/src/util/lp/lar_solver.cpp +++ b/src/util/lp/lar_solver.cpp @@ -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> & 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::one(), ul.upper_bound_witness())); evidence.push_back(std::make_pair(-numeric_traits::one(), ul.low_bound_witness())); } @@ -325,8 +343,7 @@ vector 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 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(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 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 & coeffs, const mpq & a, unsigned j) { +void lar_solver::register_monoid_in_map(std::unordered_map & 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 & void lar_solver::substitute_terms_in_linear_expression(const vector>& left_side_with_terms, - vector> &left_side, mpq & right_side) const { + vector> &left_side, mpq & free_coeff) const { std::unordered_map 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= 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(-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(-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 & 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::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_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(-1); + return m_columns_to_ul_pairs()[j].m_i != static_cast(-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 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 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(-1))); + m_columns_to_ul_pairs.push_back(ul_pair(static_cast(-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(right_side, zero_of_type()); 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(right_side, zero_of_type()); 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(right_side, zero_of_type()); 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(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); } diff --git a/src/util/lp/lar_solver.h b/src/util/lp/lar_solver.h index 7f700bd75..41794fef6 100644 --- a/src/util/lp/lar_solver.h +++ b/src/util/lp/lar_solver.h @@ -61,11 +61,11 @@ class lar_solver : public column_namer { }; //////////////////// fields ////////////////////////// lp_settings m_settings; - stacked_value m_status; + lp_status m_status; stacked_value m_simplex_strategy; std::unordered_map m_ext_vars_to_columns; vector m_columns_to_ext_vars_or_term_indices; - stacked_vector m_vars_to_ul_pairs; + stacked_vector m_columns_to_ul_pairs; vector m_constraints; stacked_value 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 m_tracker_of_x_change; + int_set m_inf_int_set; ////////////////// methods //////////////////////////////// static_matrix> & A_r() { return m_mpq_lar_core_solver.m_r_A;} static_matrix> 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 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>& left_side_with_terms, - vector> &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(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>& left_side_with_terms, + vector> &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> & evidence, lconstraint_kind & the_kind_of_sum) const; static void register_in_map(std::unordered_map & coeffs, const lar_base_constraint & cn, const mpq & a); - static void register_one_coeff_in_map(std::unordered_map & coeffs, const mpq & a, unsigned j); + static void register_monoid_in_map(std::unordered_map & coeffs, const mpq & a, unsigned j); bool the_left_sides_sum_to_zero(const vector> & 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_map) const; void fill_var_set_for_random_update(unsigned sz, var_index const * vars, vector& column_list) { @@ -1283,22 +1265,12 @@ public: return m_vars_to_ul_pairs()[j].m_i != static_cast(-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(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 * 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> 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> 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; }; } diff --git a/src/util/lp/lar_term.h b/src/util/lp/lar_term.h index 659f56f17..1c4a38801 100644 --- a/src/util/lp/lar_term.h +++ b/src/util/lp/lar_term.h @@ -25,7 +25,7 @@ struct lar_term { std::unordered_map 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>& 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(); + } + }; } diff --git a/src/util/lp/lp_core_solver_base.h b/src/util/lp/lp_core_solver_base.h index 5accb4172..6bc3f2daf 100644 --- a/src/util/lp/lp_core_solver_base.h +++ b/src/util/lp/lp_core_solver_base.h @@ -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 * m_tracker_of_x_change; + + void set_tracker_of_x(std::function* 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]; + } }; } diff --git a/src/util/lp/lp_core_solver_base.hpp b/src/util/lp/lp_core_solver_base.hpp index a7fdaa2df..51f86a8b5 100644 --- a/src/util/lp/lp_core_solver_base.hpp +++ b/src/util/lp/lp_core_solver_base.hpp @@ -39,7 +39,7 @@ lp_core_solver_base(static_matrix & A, const vector & 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 & 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 void lp_core_solver_base::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); diff --git a/src/util/lp/lp_dual_core_solver.hpp b/src/util/lp/lp_dual_core_solver.hpp index 5c73df04a..765b80db8 100644 --- a/src/util/lp/lp_dual_core_solver.hpp +++ b/src/util/lp/lp_dual_core_solver.hpp @@ -97,11 +97,11 @@ template void lp_dual_core_solver::start_with_ini } template bool lp_dual_core_solver::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 void lp_dual_core_solver::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 void lp_dual_core_solver::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 void lp_dual_core_solver::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 bool lp_dual_core_solver::delta_keeps_th } template void lp_dual_core_solver::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 bool lp_dual_core_solver::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 void lp_dual_core_solver::update_xb_afte template void lp_dual_core_solver::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 void lp_dual_core_solver::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); } diff --git a/src/util/lp/lp_dual_simplex.hpp b/src/util/lp/lp_dual_simplex.hpp index e334ded90..96fb65399 100644 --- a/src/util/lp/lp_dual_simplex.hpp +++ b/src/util/lp/lp_dual_simplex.hpp @@ -22,23 +22,23 @@ namespace lp{ template void lp_dual_simplex::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 void lp_dual_simplex::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 void lp_dual_simplex::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 void lp_dual_simplex::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 void lp_dual_simplex::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(); } } diff --git a/src/util/lp/lp_primal_core_solver.h b/src/util/lp/lp_primal_core_solver.h index c5273cacf..8df141551 100644 --- a/src/util/lp/lp_primal_core_solver.h +++ b/src/util/lp/lp_primal_core_solver.h @@ -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(); } - this->set_status(UNKNOWN); + this->set_status(lp_status::UNKNOWN); } // constructor diff --git a/src/util/lp/lp_primal_core_solver.hpp b/src/util/lp/lp_primal_core_solver.hpp index ceec1edc9..000d274e1 100644 --- a/src/util/lp/lp_primal_core_solver.hpp +++ b/src/util/lp/lp_primal_core_solver.hpp @@ -713,14 +713,14 @@ template void lp_primal_core_solver::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 void lp_primal_core_solver::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 void lp_primal_core_solver::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 void lp_primal_core_solver::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 void lp_primal_core_solver::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 void lp_primal_core_solver::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 void lp_primal_core_solver::push_forw template unsigned lp_primal_core_solver::get_number_of_non_basic_column_to_try_for_enter() { unsigned ret = static_cast(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 unsigned lp_primal_core_solver::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::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 unsigned lp_primal_core_solver::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::precise()) { @@ -892,7 +892,7 @@ template unsigned lp_primal_core_solver::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 unsigned lp_primal_core_solver::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 unsigned lp_primal_core_solver::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::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 unsigned lp_primal_core_solver::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 unsigned lp_primal_core_solver::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 T lp_primal_core_solver::calculate_no template void lp_primal_core_solver::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 void lp_primal_core_solver::fill_breakpo template bool lp_primal_core_solver::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::init_infeasibility_cost_for_column(unsigned j) { } if (numeric_traits::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]; diff --git a/src/util/lp/lp_primal_core_solver_tableau.h b/src/util/lp/lp_primal_core_solver_tableau.h index 51ace29a6..b6566f7c3 100644 --- a/src/util/lp/lp_primal_core_solver_tableau.h +++ b/src/util/lp/lp_primal_core_solver_tableau.h @@ -35,7 +35,7 @@ template void lp_primal_core_solver::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 unsigned lp_primal_core_solver::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::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::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::precise()) { @@ -127,7 +127,7 @@ unsigned lp_primal_core_solver::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::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::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::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::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); diff --git a/src/util/lp/lp_primal_simplex.hpp b/src/util/lp/lp_primal_simplex.hpp index 63ccd469d..c179a9fd1 100644 --- a/src/util/lp/lp_primal_simplex.hpp +++ b/src/util/lp/lp_primal_simplex.hpp @@ -231,7 +231,7 @@ template void lp_primal_simplex::fill_A_x_and_bas template void lp_primal_simplex::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(); diff --git a/src/util/lp/lp_settings.h b/src/util/lp/lp_settings.h index a3801f81a..063001d53 100644 --- a/src/util/lp/lp_settings.h +++ b/src/util/lp/lp_settings.h @@ -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, diff --git a/src/util/lp/lp_settings.hpp b/src/util/lp/lp_settings.hpp index 3a25e64df..68b08490b 100644 --- a/src/util/lp/lp_settings.hpp +++ b/src/util/lp/lp_settings.hpp @@ -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(); } diff --git a/src/util/lp/lp_solver.hpp b/src/util/lp/lp_solver.hpp index bfa69c5fc..071e4232f 100644 --- a/src/util/lp/lp_solver.hpp +++ b/src/util/lp/lp_solver.hpp @@ -238,7 +238,7 @@ template bool lp_solver::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 bool lp_solver::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 bool lp_solver::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 bool lp_solver::row_ge_is_obsolete(std: T rs = m_constraints[row_index].m_rs; if (row_is_zero(row)) { if (rs > zero_of_type()) - m_status = INFEASIBLE; + m_status = lp_status::INFEASIBLE; return true; } @@ -288,7 +288,7 @@ template bool lp_solver::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 bool lp_solver::row_le_is_obsolete(std:: T rs = m_constraints[row_index].m_rs; if (row_is_zero(row)) { if (rs < zero_of_type()) - m_status = INFEASIBLE; + m_status = lp_status::INFEASIBLE; return true; } diff --git a/src/util/lp/quick_xplain.cpp b/src/util/lp/quick_xplain.cpp index 256619221..c4a89ae40 100644 --- a/src/util/lp/quick_xplain.cpp +++ b/src/util/lp/quick_xplain.cpp @@ -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 & 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());