From 4ecf186580fdcdfd381b95bc32fbadceb15a4c54 Mon Sep 17 00:00:00 2001 From: Nikolaj Bjorner Date: Wed, 7 Feb 2018 15:43:33 -0800 Subject: [PATCH] remove files Signed-off-by: Nikolaj Bjorner --- src/util/lp/int_solver.cpp | 1155 -------------------- src/util/lp/lar_solver.cpp | 2089 ------------------------------------ 2 files changed, 3244 deletions(-) delete mode 100644 src/util/lp/int_solver.cpp delete mode 100644 src/util/lp/lar_solver.cpp diff --git a/src/util/lp/int_solver.cpp b/src/util/lp/int_solver.cpp deleted file mode 100644 index 87fafee77..000000000 --- a/src/util/lp/int_solver.cpp +++ /dev/null @@ -1,1155 +0,0 @@ -/* - Copyright (c) 2017 Microsoft Corporation - Author: Lev Nachmanson -*/ - -#include "util/lp/int_solver.h" -#include "util/lp/lar_solver.h" -#include "util/lp/cut_solver.h" -#include -namespace lp { - -void int_solver::failed() { - auto & lcs = m_lar_solver->m_mpq_lar_core_solver; - - for (unsigned j : m_old_values_set.m_index) { - lcs.m_r_x[j] = m_old_values_data[j]; - lp_assert(lcs.m_r_solver.column_is_feasible(j)); - lcs.m_r_solver.remove_column_from_inf_set(j); - } - lp_assert(lcs.m_r_solver.calc_current_x_is_feasible_include_non_basis()); - lp_assert(lcs.m_r_solver.current_x_is_feasible()); - m_old_values_set.clear(); -} - -void int_solver::trace_inf_rows() const { - TRACE("arith_int_rows", - unsigned num = m_lar_solver->A_r().column_count(); - for (unsigned v = 0; v < num; v++) { - if (is_int(v) && !get_value(v).is_int()) { - display_column(tout, v); - } - } - - num = 0; - for (unsigned i = 0; i < m_lar_solver->A_r().row_count(); i++) { - unsigned j = m_lar_solver->m_mpq_lar_core_solver.m_r_basis[i]; - if (column_is_int_inf(j)) { - num++; - iterator_on_row it(m_lar_solver->A_r().m_rows[i]); - m_lar_solver->print_linear_iterator(&it, tout); - tout << "\n"; - } - } - 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; -} - -bool int_solver::has_inf_int() const { - return !inf_int_set().is_empty(); -} - -int int_solver::find_inf_int_base_column() { - 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() % 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() { - int result = -1; - mpq range; - mpq new_range; - mpq small_range_thresold(1024); - unsigned n; - lar_core_solver & lcs = m_lar_solver->m_mpq_lar_core_solver; - - for (int j : inf_int_set().m_index) { - lp_assert(is_base(j) && column_is_int_inf(j)); - lp_assert(!is_fixed(j)); - if (!is_boxed(j)) - continue; - new_range = lcs.m_r_upper_bounds()[j].x - lcs.m_r_low_bounds()[j].x; - if (new_range > small_range_thresold) - continue; - if (result == -1) { - result = j; - range = new_range; - n = 1; - continue; - } - if (new_range < range) { - n = 1; - result = j; - range = new_range; - continue; - } - if (new_range == range) { - lp_assert(n >= 1); - if (settings().random_next() % (++n) == 0) { - result = j; - continue; - } - } - } - return result; - -} - -bool int_solver::is_gomory_cut_target(linear_combination_iterator &iter) { - unsigned j; - lp_assert(iter.is_reset()); - // All non base variables must be at their bounds and assigned to rationals (that is, infinitesimals are not allowed). - while (iter.next(j)) { - if (is_base(j)) continue; - if (!is_zero(get_value(j).y)) { - TRACE("gomory_cut", tout << "row is not gomory cut target:\n"; - display_column(tout, j); - tout << "infinitesimal: " << !is_zero(get_value(j).y) << "\n";); - iter.reset(); - return false; - } - } - iter.reset(); - return true; -} - - -void int_solver::real_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& pol, explanation & expl, unsigned gomory_cut_inf_column) { - TRACE("gomory_cut_detail_real", tout << "real\n";); - mpq f_0 = fractional_part(get_value(gomory_cut_inf_column)); - mpq new_a; - if (at_low(x_j)) { - if (a.is_pos()) { - new_a = a / (1 - f_0); - } - else { - new_a = a / f_0; - new_a.neg(); - } - k.addmul(new_a, low_bound(x_j).x); // is it a faster operation than - // k += lower_bound(x_j).x * new_a; - 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.addmul(new_a, upper_bound(x_j).x); // k += upper_bound(x_j).x * new_a; - expl.push_justification(column_upper_bound_constraint(x_j), new_a); - } - TRACE("gomory_cut_detail_real", tout << a << "*v" << x_j << " k: " << k << "\n";); - pol.add_monomial(new_a, x_j); -} - -constraint_index int_solver::column_upper_bound_constraint(unsigned j) const { - return m_lar_solver->get_column_upper_bound_witness(j); -} - -constraint_index int_solver::column_low_bound_constraint(unsigned j) const { - return m_lar_solver->get_column_low_bound_witness(j); -} - - -void int_solver::int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term & t, explanation& expl, mpq & lcm_den, unsigned inf_column) { - lp_assert(is_int(x_j)); - lp_assert(!a.is_int()); - mpq f_0 = fractional_part(get_value(inf_column)); - lp_assert(f_0 > zero_of_type() && f_0 < one_of_type()); - mpq f_j = fractional_part(a); - TRACE("gomory_cut_detail", - tout << a << " x_j" << x_j << " k = " << k << "\n"; - tout << "f_j: " << f_j << "\n"; - tout << "f_0: " << f_0 << "\n"; - tout << "1 - f_0: " << 1 - f_0 << "\n"; - tout << "at_low(" << x_j << ") = " << at_low(x_j) << std::endl; - ); - lp_assert (!f_j.is_zero()); - mpq new_a; - if (at_low(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 { - new_a = (1 - f_j) / f_0; - } - k.addmul(new_a, low_bound(x_j).x); - expl.push_justification(column_low_bound_constraint(x_j), new_a); - } - else { - lp_assert(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); - } - TRACE("gomory_cut_detail", tout << "new_a: " << new_a << " k: " << k << "\n";); - t.add_monomial(new_a, x_j); - lcm_den = lcm(lcm_den, denominator(new_a)); -} - -lia_move int_solver::report_conflict_from_gomory_cut(mpq & k) { - TRACE("empty_pol",); - 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(vector> & pol, - lar_term & t, - mpq &k, - bool some_ints, - mpq & lcm_den) { - if (!some_ints) - return; - - t.clear(); - if (pol.size() == 1) { - unsigned v = pol[0].second; - lp_assert(is_int(v)); - bool k_is_int = k.is_int(); - const mpq& a = pol[0].first; - k /= a; - if (a.is_pos()) { // we have av >= k - if (!k_is_int) - k = ceil(k); - // switch size - t.add_monomial(- mpq(1), v); - k.neg(); - } else { - if (!k_is_int) - k = floor(k); - t.add_monomial(mpq(1), v); - } - } else if (some_ints) { - lcm_den = lcm(lcm_den, denominator(k)); - lp_assert(lcm_den.is_pos()); - if (!lcm_den.is_one()) { - // normalize coefficients of integer parameters to be integers. - for (auto & pi: pol) { - pi.first *= lcm_den; - SASSERT(!is_int(pi.second) || pi.first.is_int()); - } - k *= lcm_den; - } - // negate everything to return -pol <= -k - for (const auto & pi: pol) - t.add_monomial(-pi.first, pi.second); - k.neg(); - } -} - -bool int_solver::current_solution_is_inf_on_cut(const lar_term& t, const mpq& k) const { - const auto & x = m_lar_solver->m_mpq_lar_core_solver.m_r_x; - impq v = t.apply(x); - TRACE( - "current_solution_is_inf_on_cut", tout << "v = " << v << " k = " << k << std::endl; - if (v <=k) { - tout << "v <= k - it should not happen!\n"; - } - ); - - return v > k; -} - -void int_solver::adjust_term_and_k_for_some_ints_case_gomory(lar_term& t, mpq& k, mpq &lcm_den) { - lp_assert(!t.is_empty()); - auto pol = t.coeffs_as_vector(); - t.clear(); - if (pol.size() == 1) { - TRACE("gomory_cut_detail", tout << "pol.size() is 1" << std::endl;); - unsigned v = pol[0].second; - lp_assert(is_int(v)); - const mpq& a = pol[0].first; - k /= a; - if (a.is_pos()) { // we have av >= k - if (!k.is_int()) - k = ceil(k); - // switch size - t.add_monomial(- mpq(1), v); - k.neg(); - } else { - if (!k.is_int()) - k = floor(k); - t.add_monomial(mpq(1), v); - } - } else { - TRACE("gomory_cut_detail", tout << "pol.size() > 1" << std::endl;); - lcm_den = lcm(lcm_den, denominator(k)); - lp_assert(lcm_den.is_pos()); - if (!lcm_den.is_one()) { - // normalize coefficients of integer parameters to be integers. - for (auto & pi: pol) { - pi.first *= lcm_den; - SASSERT(!is_int(pi.second) || pi.first.is_int()); - } - k *= lcm_den; - } - // negate everything to return -pol <= -k - for (const auto & pi: pol) - t.add_monomial(-pi.first, pi.second); - k.neg(); - } - TRACE("gomory_cut_detail", tout << "k = " << k << std::endl;); - lp_assert(k.is_int()); -} - - - - -lia_move int_solver::mk_gomory_cut(lar_term& t, mpq& k, explanation & expl, unsigned inf_col, linear_combination_iterator& iter) { - - lp_assert(column_is_int_inf(inf_col)); - - TRACE("gomory_cut", - tout << "applying cut at:\n"; m_lar_solver->print_linear_iterator_indices_only(&iter, tout); tout << std::endl; - iter.reset(); - unsigned j; - while(iter.next(j)) { - m_lar_solver->m_mpq_lar_core_solver.m_r_solver.print_column_info(j, tout); - } - iter.reset(); - tout << "inf_col = " << inf_col << std::endl; - ); - - // gomory will be t >= k - k = 1; - mpq lcm_den(1); - unsigned x_j; - mpq a; - bool some_int_columns = false; - lp_assert(iter.is_reset()); - while (iter.next(a, x_j)) { - if (x_j == inf_col) - 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, inf_col); - else { - if (a.is_int()) continue; // f_j will be zero and no monomial will be added - some_int_columns = true; - int_case_in_gomory_cut(a, x_j, k, t, expl, lcm_den, inf_col); - } - } - - if (t.is_empty()) - return report_conflict_from_gomory_cut(k); - if (some_int_columns) - adjust_term_and_k_for_some_ints_case_gomory(t, k, lcm_den); - - lp_assert(current_solution_is_inf_on_cut(t, k)); - m_lar_solver->subs_term_columns(t); - return lia_move::cut; - -} - -void int_solver::init_check_data() { - 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_free_var_in_gomory_row(linear_combination_iterator& iter) { - unsigned j; - while(iter.next(j)) { - if (!is_base(j) && is_free(j)) - return static_cast(j); - } - iter.reset(); - return -1; -} - -lia_move int_solver::proceed_with_gomory_cut(lar_term& t, mpq& k, explanation& ex, unsigned j) { - lia_move ret; - linear_combination_iterator* iter = m_lar_solver->get_iterator_on_row(row_of_basic_column(j)); - int free_j = find_free_var_in_gomory_row(*iter); - if (free_j != -1) { - ret = create_branch_on_column(j, t, k, true); - } else if (!is_gomory_cut_target(*iter)) { - ret = create_branch_on_column(j, t, k, false); - } else { - ret = mk_gomory_cut(t, k, ex, j, *iter); - } - delete iter; - return ret; -} - - -unsigned int_solver::row_of_basic_column(unsigned j) const { - return m_lar_solver->m_mpq_lar_core_solver.m_r_heading[j]; -} - -template -void int_solver::fill_cut_solver(cut_solver & cs) { - for (lar_base_constraint * c : m_lar_solver->constraints()) - fill_cut_solver_for_constraint(c, cs); -} - -template -void int_solver::fill_cut_solver_for_constraint(const lar_base_constraint* c, cut_solver & cs) { - vector> coeffs; - T rs; - get_int_coeffs_from_constraint(c, coeffs, rs); - cs.add_ineq(coeffs, rs); -} - - -// it produces an inequality coeff*x <= rs -template -void int_solver::get_int_coeffs_from_constraint(const lar_base_constraint* c, vector>& coeffs, T & rs) { - lp_assert(c->m_kind != EQ); // it is not implemented, we need to create two inequalities in this case - int sign = ((int)c->m_kind > 0) ? -1 : 1; - vector> lhs = c->get_left_side_coefficients(); - T den = denominator(c->m_right_side); - for (auto & kv : lhs) { - den = lcm(den, denominator(kv.first)); - } - lp_assert(den > 0); - for (auto& kv : lhs) { - coeffs.push_back(std::make_pair(den * kv.first * sign, kv.second)); - } - rs = den * c->m_right_side * sign; - if (kind_is_strict(c->m_kind)) - rs--; -} - - -// this will allow to enable and disable tracking of the pivot rows -struct pivoted_rows_tracking_control { - lar_solver * m_lar_solver; - bool m_track_pivoted_rows; - pivoted_rows_tracking_control(lar_solver* ls) : - m_lar_solver(ls), - m_track_pivoted_rows(ls->get_track_pivoted_rows()) - { - TRACE("pivoted_rows", tout << "pivoted rows = " << ls->m_mpq_lar_core_solver.m_r_solver.m_pivoted_rows->size() << std::endl;); - m_lar_solver->set_track_pivoted_rows(false); - } - ~pivoted_rows_tracking_control() { - TRACE("pivoted_rows", tout << "pivoted rows = " << m_lar_solver->m_mpq_lar_core_solver.m_r_solver.m_pivoted_rows->size() << std::endl;); - m_lar_solver->set_track_pivoted_rows(m_track_pivoted_rows); - } -}; - -lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex) { - init_check_data(); - lp_assert(inf_int_set_is_correct()); - // it is mostly a reimplementation of - // final_check_status theory_arith::check_int_feasibility() - // from theory_arith_int.h - if (!has_inf_int()) - return lia_move::ok; - if (settings().m_run_gcd_test) - if (!gcd_test(ex)) - return lia_move::conflict; - pivoted_rows_tracking_control pc(m_lar_solver); - /* if (m_params.m_arith_euclidean_solver) apply_euclidean_solver(); */ - //m_lar_solver->pivot_fixed_vars_from_basis(); - patch_int_infeasible_nbasic_columns(); - if (!has_inf_int()) - return lia_move::ok; - - // lp_assert(non_basic_columns_are_at_bounds()); - TRACE("gomory_cut", tout << m_branch_cut_counter+1 << ", " << settings().m_int_branch_cut_gomory_threshold << std::endl;); - if (++m_branch_cut_counter > 0) { // testing cut_solver - cut_solver cs; - fill_cut_solver(cs); - } else - if ((++m_branch_cut_counter) % settings().m_int_branch_cut_gomory_threshold == 0) { - if (move_non_basic_columns_to_bounds()) { - lp_status st = m_lar_solver->find_feasible_solution(); - lp_assert(non_basic_columns_are_at_bounds()); - if (st != lp_status::FEASIBLE && st != lp_status::OPTIMAL) { - TRACE("arith_int", tout << "give_up\n";); - return lia_move::give_up; - } - } - int j = find_inf_int_base_column(); - if (j == -1) return lia_move::ok; - TRACE("arith_int", tout << "j = " << j << " does not have an integer assignment: " << get_value(j) << "\n";); - return proceed_with_gomory_cut(t, k, ex, j); - } - return create_branch_on_column(find_inf_int_base_column(), t, k, false); -} - -bool int_solver::move_non_basic_column_to_bounds(unsigned j) { - auto & lcs = m_lar_solver->m_mpq_lar_core_solver; - auto & val = lcs.m_r_x[j]; - switch (lcs.m_column_types()[j]) { - case column_type::boxed: - if (val != lcs.m_r_low_bounds()[j] && val != lcs.m_r_upper_bounds()[j]) { - if (random() % 2 == 0) - set_value_for_nbasic_column(j, lcs.m_r_low_bounds()[j]); - else - set_value_for_nbasic_column(j, lcs.m_r_upper_bounds()[j]); - return true; - } - break; - case column_type::low_bound: - if (val != lcs.m_r_low_bounds()[j]) { - set_value_for_nbasic_column(j, lcs.m_r_low_bounds()[j]); - return true; - } - break; - case column_type::upper_bound: - if (val != lcs.m_r_upper_bounds()[j]) { - set_value_for_nbasic_column(j, lcs.m_r_upper_bounds()[j]); - return true; - } - break; - default: - if (is_int(j) && !val.is_int()) { - set_value_for_nbasic_column(j, impq(floor(val))); - return true; - } - break; - } - return false; -} - -bool int_solver::move_non_basic_columns_to_bounds() { - auto & lcs = m_lar_solver->m_mpq_lar_core_solver; - bool change = false; - for (unsigned j : lcs.m_r_nbasis) { - if (move_non_basic_column_to_bounds(j)) - change = true; - } - return change; -} - -void int_solver::set_value_for_nbasic_column_ignore_old_values(unsigned j, const impq & new_val) { - lp_assert(!is_base(j)); - auto & x = m_lar_solver->m_mpq_lar_core_solver.m_r_x[j]; - auto delta = new_val - x; - x = new_val; - update_column_in_int_inf_set(j); - m_lar_solver->change_basic_columns_dependend_on_a_given_nb_column(j, delta); -} - - -void int_solver::set_value_for_nbasic_column(unsigned j, const impq & new_val) { - lp_assert(!is_base(j)); - auto & x = m_lar_solver->m_mpq_lar_core_solver.m_r_x[j]; - if (m_lar_solver->has_int_var() && !m_old_values_set.contains(j)) { - m_old_values_set.insert(j); - m_old_values_data[j] = x; - } - auto delta = new_val - x; - x = new_val; - update_column_in_int_inf_set(j); - m_lar_solver->change_basic_columns_dependend_on_a_given_nb_column(j, delta); -} - -void int_solver::patch_int_infeasible_non_basic_column(unsigned j) { - if (!is_int(j)) return; - bool inf_l, inf_u; - impq l, u; - mpq m; - if (!get_value(j).is_int() || !get_freedom_interval_for_column(j, inf_l, l, inf_u, u, m)) { - move_non_basic_column_to_bounds(j); - return; - } - auto & lcs = m_lar_solver->m_mpq_lar_core_solver; - impq & val = lcs.m_r_x[j]; - bool val_is_int = val.is_int(); - bool m_is_one = m.is_one(); - if (m.is_one() && val_is_int) - return; - // check whether value of j is already a multiple of m. - if (val_is_int && (val.x / m).is_int()) - return; - TRACE("patch_int", - tout << "TARGET j" << j << " -> ["; - if (inf_l) tout << "-oo"; else tout << l; - tout << ", "; - if (inf_u) tout << "oo"; else tout << u; - tout << "]"; - tout << ", m: " << m << ", val: " << val << ", is_int: " << m_lar_solver->column_is_int(j) << "\n";); - if (!inf_l) { - l = m_is_one ? ceil(l) : m * ceil(l / m); - if (inf_u || l <= u) { - TRACE("patch_int", - tout << "patching with l: " << l << '\n';); - - set_value_for_nbasic_column(j, l); - } - else { - TRACE("patch_int", - tout << "not patching " << l << "\n";); - } - } - else if (!inf_u) { - u = m_is_one ? floor(u) : m * floor(u / m); - set_value_for_nbasic_column(j, u); - TRACE("patch_int", - tout << "patching with u: " << u << '\n';); - } - else { - set_value_for_nbasic_column(j, impq(0)); - TRACE("patch_int", - tout << "patching with 0\n";); - } -} -void int_solver::patch_int_infeasible_nbasic_columns() { - lp_assert(is_feasible()); - for (unsigned j : m_lar_solver->m_mpq_lar_core_solver.m_r_nbasis) { - patch_int_infeasible_non_basic_column(j); - if (!is_feasible()) - break; - } - if (!is_feasible()) { - move_non_basic_columns_to_bounds(); - m_lar_solver->find_feasible_solution(); - } - lp_assert(is_feasible() && inf_int_set_is_correct()); -} - -mpq get_denominators_lcm(iterator_on_row &it) { - mpq r(1); - mpq a; - unsigned j; - while (it.next(a, j)) { - r = lcm(r, denominator(a)); - } - return r; -} - -bool int_solver::gcd_test_for_row(static_matrix> & A, unsigned i, explanation & ex) { - iterator_on_row it(A.m_rows[i]); - mpq lcm_den = get_denominators_lcm(it); - mpq consts(0); - mpq gcds(0); - mpq least_coeff(0); - bool least_coeff_is_bounded = false; - mpq a; - unsigned j; - while (it.next(a, j)) { - if (m_lar_solver->column_is_fixed(j)) { - mpq aux = lcm_den * a; - consts += aux * m_lar_solver->column_low_bound(j).x; - } - else if (m_lar_solver->column_is_real(j)) { - return true; - } - else if (gcds.is_zero()) { - gcds = abs(lcm_den * a); - least_coeff = gcds; - least_coeff_is_bounded = m_lar_solver->column_is_bounded(j); - } - else { - mpq aux = abs(lcm_den * a); - gcds = gcd(gcds, aux); - if (aux < least_coeff) { - least_coeff = aux; - least_coeff_is_bounded = m_lar_solver->column_is_bounded(j); - } - else if (least_coeff_is_bounded && aux == least_coeff) { - least_coeff_is_bounded = m_lar_solver->column_is_bounded(j); - } - } - SASSERT(gcds.is_int()); - SASSERT(least_coeff.is_int()); - TRACE("gcd_test_bug", tout << "coeff: " << a << ", gcds: " << gcds - << " least_coeff: " << least_coeff << " consts: " << consts << "\n";); - - } - - if (gcds.is_zero()) { - // All variables are fixed. - // This theory guarantees that the assignment satisfies each row, and - // fixed integer variables are assigned to integer values. - return true; - } - - if (!(consts / gcds).is_int()) - fill_explanation_from_fixed_columns(it, ex); - - if (least_coeff.is_one() && !least_coeff_is_bounded) { - SASSERT(gcds.is_one()); - return true; - } - - if (least_coeff_is_bounded) { - return ext_gcd_test(it, least_coeff, lcm_den, consts, ex); - } - return true; -} - -void int_solver::add_to_explanation_from_fixed_or_boxed_column(unsigned j, explanation & ex) { - constraint_index lc, uc; - m_lar_solver->get_bound_constraint_witnesses_for_column(j, lc, uc); - ex.m_explanation.push_back(std::make_pair(mpq(1), lc)); - ex.m_explanation.push_back(std::make_pair(mpq(1), uc)); -} -void int_solver::fill_explanation_from_fixed_columns(iterator_on_row & it, explanation & ex) { - it.reset(); - unsigned j; - while (it.next(j)) { - if (!m_lar_solver->column_is_fixed(j)) - continue; - add_to_explanation_from_fixed_or_boxed_column(j, ex); - } -} - -bool int_solver::gcd_test(explanation & ex) { - auto & A = m_lar_solver->A_r(); // getting the matrix - for (unsigned i = 0; i < A.row_count(); i++) - if (!gcd_test_for_row(A, i, ex)) { - std::cout << "false from gcd_test\n" ; - return false; - } - - return true; -} - -bool int_solver::ext_gcd_test(iterator_on_row & it, - mpq const & least_coeff, - mpq const & lcm_den, - mpq const & consts, explanation& ex) { - mpq gcds(0); - mpq l(consts); - mpq u(consts); - - it.reset(); - mpq a; - unsigned j; - while (it.next(a, j)) { - if (m_lar_solver->column_is_fixed(j)) - continue; - SASSERT(!m_lar_solver->column_is_real(j)); - mpq ncoeff = lcm_den * a; - SASSERT(ncoeff.is_int()); - mpq abs_ncoeff = abs(ncoeff); - if (abs_ncoeff == least_coeff) { - SASSERT(m_lar_solver->column_is_bounded(j)); - if (ncoeff.is_pos()) { - // l += ncoeff * m_lar_solver->column_low_bound(j).x; - l.addmul(ncoeff, m_lar_solver->column_low_bound(j).x); - // u += ncoeff * m_lar_solver->column_upper_bound(j).x; - u.addmul(ncoeff, m_lar_solver->column_upper_bound(j).x); - } - else { - // l += ncoeff * upper_bound(j).get_rational(); - l.addmul(ncoeff, m_lar_solver->column_upper_bound(j).x); - // u += ncoeff * low_bound(j).get_rational(); - u.addmul(ncoeff, m_lar_solver->column_low_bound(j).x); - } - add_to_explanation_from_fixed_or_boxed_column(j, ex); - } - else if (gcds.is_zero()) { - gcds = abs_ncoeff; - } - else { - gcds = gcd(gcds, abs_ncoeff); - } - SASSERT(gcds.is_int()); - } - - if (gcds.is_zero()) { - return true; - } - - mpq l1 = ceil(l/gcds); - mpq u1 = floor(u/gcds); - - if (u1 < l1) { - fill_explanation_from_fixed_columns(it, ex); - return false; - } - - return true; - -} - -linear_combination_iterator * int_solver::get_column_iterator(unsigned j) { - if (m_lar_solver->use_tableau()) - return new iterator_on_column(m_lar_solver->A_r().m_columns[j], m_lar_solver->A_r()); - return new iterator_on_indexed_vector(m_lar_solver->get_column_in_lu_mode(j)); -} - - -int_solver::int_solver(lar_solver* lar_slv) : - m_lar_solver(lar_slv), - m_branch_cut_counter(0) { - lp_assert(m_old_values_set.size() == 0); - m_old_values_set.resize(lar_slv->A_r().column_count()); - m_old_values_data.resize(lar_slv->A_r().column_count(), zero_of_type()); - m_lar_solver->set_int_solver(this); -} - -bool int_solver::has_low(unsigned j) const { - switch (m_lar_solver->m_mpq_lar_core_solver.m_column_types()[j]) { - case column_type::fixed: - case column_type::boxed: - case column_type::low_bound: - return true; - default: - return false; - } -} - -bool int_solver::has_upper(unsigned j) const { - switch (m_lar_solver->m_mpq_lar_core_solver.m_column_types()[j]) { - case column_type::fixed: - case column_type::boxed: - case column_type::upper_bound: - return true; - default: - return false; - } -} - - -void set_lower(impq & l, - bool & inf_l, - impq const & v ) { - if (inf_l || v > l) { - l = v; - inf_l = false; - } -} - -void set_upper(impq & u, - bool & inf_u, - impq const & v) { - if (inf_u || v < u) { - u = v; - inf_u = false; - } -} - -bool int_solver::get_freedom_interval_for_column(unsigned j, bool & inf_l, impq & l, bool & inf_u, impq & u, mpq & m) { - auto & lcs = m_lar_solver->m_mpq_lar_core_solver; - if (lcs.m_r_heading[j] >= 0) // the basic var - return false; - - impq const & xj = get_value(j); - linear_combination_iterator *it = get_column_iterator(j); - - inf_l = true; - inf_u = true; - l = u = zero_of_type(); - m = mpq(1); - - if (has_low(j)) { - set_lower(l, inf_l, low_bound(j)); - } - if (has_upper(j)) { - set_upper(u, inf_u, upper_bound(j)); - } - - mpq a; // the coefficient in the column - unsigned row_index; - while (it->next(a, row_index)) { - unsigned i = lcs.m_r_basis[row_index]; - impq const & xi = get_value(i); - if (is_int(i) && is_int(j) && !a.is_int()) - m = lcm(m, denominator(a)); - if (a.is_neg()) { - if (has_low(i)) - set_lower(l, inf_l, xj + (xi - lcs.m_r_low_bounds()[i]) / a); - - if (has_upper(i)) - set_upper(u, inf_u, xj + (xi - lcs.m_r_upper_bounds()[i]) / a); - } - else { - if (has_upper(i)) - set_lower(l, inf_l, xj + (xi - lcs.m_r_upper_bounds()[i]) / a); - if (has_low(i)) - set_upper(u, inf_u, xj + (xi - lcs.m_r_low_bounds()[i]) / a); - } - if (!inf_l && !inf_u && l == u) break;; - } - - delete it; - TRACE("freedom_interval", - tout << "freedom variable for:\n"; - tout << m_lar_solver->get_column_name(j); - tout << "["; - if (inf_l) tout << "-oo"; else tout << l; - tout << "; "; - if (inf_u) tout << "oo"; else tout << u; - tout << "]\n"; - tout << "val = " << get_value(j) << "\n"; - ); - lp_assert(inf_l || l <= get_value(j)); - lp_assert(inf_u || u >= get_value(j)); - return true; - -} - -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->column_value_is_int(j); -} - - - -bool int_solver::is_feasible() const { - auto & lcs = m_lar_solver->m_mpq_lar_core_solver; - lp_assert( - lcs.m_r_solver.calc_current_x_is_feasible_include_non_basis() == - lcs.m_r_solver.current_x_is_feasible()); - return lcs.m_r_solver.current_x_is_feasible(); -} -const impq & int_solver::get_value(unsigned j) const { - return m_lar_solver->m_mpq_lar_core_solver.m_r_x[j]; -} - -void int_solver::display_column(std::ostream & out, unsigned j) const { - m_lar_solver->m_mpq_lar_core_solver.m_r_solver.print_column_info(j, out); -} - -bool int_solver::inf_int_set_is_correct() const { - for (unsigned j = 0; j < m_lar_solver->A_r().column_count(); j++) { - if (inf_int_set().contains(j) != (is_int(j) && (!value_is_int(j)))) { - TRACE("arith_int", tout << "j= " << j << " inf_int_set().contains(j) = " << inf_int_set().contains(j) << ", is_int(j) = " << is_int(j) << "\nvalue_is_int(j) = " << value_is_int(j) << ", val = " << get_value(j) << std::endl;); - return false; - } - } - return true; -} - -bool int_solver::column_is_int_inf(unsigned j) const { - return is_int(j) && (!value_is_int(j)); -} - -void int_solver::update_column_in_int_inf_set(unsigned j) { - if (is_int(j) && (!value_is_int(j))) - inf_int_set().insert(j); - else - inf_int_set().erase(j); -} - -bool int_solver::is_base(unsigned j) const { - return m_lar_solver->m_mpq_lar_core_solver.m_r_heading[j] >= 0; -} - -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_fixed(unsigned j) const { - return m_lar_solver->m_mpq_lar_core_solver.m_column_types[j] == column_type::fixed; -} - -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_low(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; -} - -unsigned int_solver::random() { - return m_lar_solver->get_core_solver().settings().random_next(); -} - -bool int_solver::shift_var(unsigned j, unsigned range) { - if (is_fixed(j) || is_base(j)) - return false; - - bool inf_l, inf_u; - impq l, u; - mpq m; - get_freedom_interval_for_column(j, inf_l, l, inf_u, u, m); - if (inf_l && inf_u) { - impq new_val = impq(random() % (range + 1)); - set_value_for_nbasic_column_ignore_old_values(j, new_val); - return true; - } - if (is_int(j)) { - if (!inf_l) { - l = ceil(l); - if (!m.is_one()) - l = m*ceil(l/m); - } - if (!inf_u) { - u = floor(u); - if (!m.is_one()) - u = m*floor(u/m); - } - } - if (!inf_l && !inf_u && l >= u) - return false; - if (inf_u) { - SASSERT(!inf_l); - impq delta = impq(random() % (range + 1)); - impq new_val = l + m*delta; - set_value_for_nbasic_column_ignore_old_values(j, new_val); - return true; - } - if (inf_l) { - SASSERT(!inf_u); - impq delta = impq(random() % (range + 1)); - impq new_val = u - m*delta; - set_value_for_nbasic_column_ignore_old_values(j, new_val); - return true; - } - if (!is_int(j)) { - SASSERT(!inf_l && !inf_u); - mpq delta = mpq(random() % (range + 1)); - impq new_val = l + ((delta * (u - l)) / mpq(range)); - set_value_for_nbasic_column_ignore_old_values(j, new_val); - return true; - } - else { - mpq r = (u.x - l.x) / m; - if (r < mpq(range)) - range = static_cast(r.get_uint64()); - impq new_val = l + m * (impq(random() % (range + 1))); - set_value_for_nbasic_column_ignore_old_values(j, new_val); - return true; - } -} - -bool int_solver::non_basic_columns_are_at_bounds() const { - auto & lcs = m_lar_solver->m_mpq_lar_core_solver; - for (unsigned j :lcs.m_r_nbasis) { - auto & val = lcs.m_r_x[j]; - switch (lcs.m_column_types()[j]) { - case column_type::boxed: - if (val != lcs.m_r_low_bounds()[j] && val != lcs.m_r_upper_bounds()[j]) - return false; - break; - case column_type::low_bound: - if (val != lcs.m_r_low_bounds()[j]) - return false; - break; - case column_type::upper_bound: - if (val != lcs.m_r_upper_bounds()[j]) - return false; - break; - default: - if (is_int(j) && !val.is_int()) { - return false; - } - } - } - return true; -} - -const impq& int_solver::low_bound(unsigned j) const { - return m_lar_solver->column_low_bound(j); -} - -lia_move int_solver::create_branch_on_column(int j, lar_term& t, mpq& k, bool free_column) const { - lp_assert(t.is_empty()); - lp_assert(j != -1); - t.add_monomial(mpq(1), m_lar_solver->adjust_column_index_to_term_index(j)); - k = free_column? mpq(0) : floor(get_value(j)); - TRACE("arith_int", tout << "branching v" << j << " = " << get_value(j) << "\n"; - display_column(tout, j); - tout << "k = " << k << std::endl; - ); - return lia_move::branch; - -} - -const impq& int_solver::upper_bound(unsigned j) const { - return m_lar_solver->column_upper_bound(j); -} -void int_solver::display_inf_or_int_inf_columns(std::ostream & out) const { - out << "int inf\n"; - for (unsigned j : m_lar_solver->m_inf_int_set.m_index) { - display_column(out, j); - } - out << "regular inf\n"; - for (unsigned j : m_lar_solver->m_mpq_lar_core_solver.m_r_solver.m_inf_set.m_index) { - display_column(out, j); - } -} -} diff --git a/src/util/lp/lar_solver.cpp b/src/util/lp/lar_solver.cpp deleted file mode 100644 index 6846717af..000000000 --- a/src/util/lp/lar_solver.cpp +++ /dev/null @@ -1,2089 +0,0 @@ -#include "util/lp/lar_solver.h" -/* - Copyright (c) 2017 Microsoft Corporation - Author: Nikolaj Bjorner, Lev Nachmanson -*/ - -namespace lp { - -unsigned lar_solver::constraint_count() const { - return m_constraints.size(); -} -const lar_base_constraint& lar_solver::get_constraint(unsigned ci) const { - return *(m_constraints[ci]); -} - -////////////////// methods //////////////////////////////// -static_matrix> & lar_solver::A_r() { return m_mpq_lar_core_solver.m_r_A;} -static_matrix> const & lar_solver::A_r() const { return m_mpq_lar_core_solver.m_r_A;} -static_matrix & lar_solver::A_d() { return m_mpq_lar_core_solver.m_d_A;} -static_matrix const & lar_solver::A_d() const { return m_mpq_lar_core_solver.m_d_A;} - -lp_settings & lar_solver::settings() { return m_settings;} - -lp_settings const & lar_solver::settings() const { return m_settings;} - -void clear() {lp_assert(false); // not implemented -} - - -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_tracker_of_x_change([&](unsigned j) { - call_assignment_tracker(j); - } - ), - m_int_solver(nullptr) -{} - -void lar_solver::set_track_pivoted_rows(bool v) { - m_mpq_lar_core_solver.m_r_solver.m_pivoted_rows = v? (& m_rows_with_changed_bounds) : nullptr; -} - -bool lar_solver::get_track_pivoted_rows() const { - return m_mpq_lar_core_solver.m_r_solver.m_pivoted_rows != nullptr; -} - - -lar_solver::~lar_solver(){ - for (auto c : m_constraints) - delete c; - for (auto t : m_terms) - delete t; -} - -bool lar_solver::is_term(var_index j) const { - return j >= m_terms_start_index && j - m_terms_start_index < m_terms.size(); -} - -unsigned lar_solver::adjust_term_index(unsigned j) const { - lp_assert(is_term(j)); - return j - m_terms_start_index; -} - - -bool lar_solver::use_lu() const { return m_settings.simplex_strategy() == simplex_strategy_enum::lu; } - -bool lar_solver::sizes_are_correct() const { - lp_assert(strategy_is_undecided() || !m_mpq_lar_core_solver.need_to_presolve_with_double_solver() || A_r().column_count() == A_d().column_count()); - lp_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_column_types.size()); - lp_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_costs.size()); - lp_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_x.size()); - return true; -} - - -void lar_solver::print_implied_bound(const implied_bound& be, std::ostream & out) const { - out << "implied bound\n"; - unsigned v = be.m_j; - if (is_term(v)) { - out << "it is a term number " << be.m_j << std::endl; - print_term(*m_terms[be.m_j - m_terms_start_index], out); - } - else { - out << get_column_name(v); - } - out << " " << lconstraint_kind_string(be.kind()) << " " << be.m_bound << std::endl; - out << "end of implied bound" << std::endl; -} - -bool lar_solver::implied_bound_is_correctly_explained(implied_bound const & be, const vector> & explanation) const { - std::unordered_map coeff_map; - auto rs_of_evidence = zero_of_type(); - unsigned n_of_G = 0, n_of_L = 0; - bool strict = false; - for (auto & it : explanation) { - mpq coeff = it.first; - constraint_index con_ind = it.second; - const auto & constr = *m_constraints[con_ind]; - lconstraint_kind kind = coeff.is_pos() ? constr.m_kind : flip_kind(constr.m_kind); - register_in_map(coeff_map, constr, coeff); - if (kind == GT || kind == LT) - strict = true; - if (kind == GE || kind == GT) n_of_G++; - else if (kind == LE || kind == LT) n_of_L++; - rs_of_evidence += coeff*constr.m_right_side; - } - lp_assert(n_of_G == 0 || n_of_L == 0); - lconstraint_kind kind = n_of_G ? GE : (n_of_L ? LE : EQ); - if (strict) - kind = static_cast((static_cast(kind) / 2)); - - if (!is_term(be.m_j)) { - if (coeff_map.size() != 1) - return false; - auto it = coeff_map.find(be.m_j); - if (it == coeff_map.end()) return false; - mpq ratio = it->second; - if (ratio < zero_of_type()) { - kind = static_cast(-kind); - } - rs_of_evidence /= ratio; - } else { - const lar_term * t = m_terms[adjust_term_index(be.m_j)]; - const auto first_coeff = *t->m_coeffs.begin(); - unsigned j = first_coeff.first; - auto it = coeff_map.find(j); - if (it == coeff_map.end()) - return false; - mpq ratio = it->second / first_coeff.second; - for (auto & p : t->m_coeffs) { - it = coeff_map.find(p.first); - if (it == coeff_map.end()) - return false; - if (p.second * ratio != it->second) - return false; - } - if (ratio < zero_of_type()) { - kind = static_cast(-kind); - } - rs_of_evidence /= ratio; - rs_of_evidence += t->m_v * ratio; - } - - return kind == be.kind() && rs_of_evidence == be.m_bound; -} - - -void lar_solver::analyze_new_bounds_on_row( - unsigned row_index, - bound_propagator & bp) { - lp_assert(!use_tableau()); - iterator_on_pivot_row it(m_mpq_lar_core_solver.get_pivot_row(), m_mpq_lar_core_solver.m_r_basis[row_index]); - - bound_analyzer_on_row ra_pos(it, - zero_of_type>(), - row_index, - bp - ); - ra_pos.analyze(); -} - -void lar_solver::analyze_new_bounds_on_row_tableau( - unsigned row_index, - bound_propagator & bp - ) { - - if (A_r().m_rows[row_index].size() > settings().max_row_length_for_bound_propagation) - return; - iterator_on_row it(A_r().m_rows[row_index]); - lp_assert(use_tableau()); - bound_analyzer_on_row::analyze_row(it, - zero_of_type>(), - row_index, - bp - ); -} - - -void lar_solver::substitute_basis_var_in_terms_for_row(unsigned i) { - // todo : create a map from term basic vars to the rows where they are used - unsigned basis_j = m_mpq_lar_core_solver.m_r_solver.m_basis[i]; - for (unsigned k = 0; k < m_terms.size(); k++) { - if (term_is_used_as_row(k)) - continue; - if (!m_terms[k]->contains(basis_j)) - continue; - m_terms[k]->subst(basis_j, m_mpq_lar_core_solver.m_r_solver.m_pivot_row); - } -} - -void lar_solver::calculate_implied_bounds_for_row(unsigned i, bound_propagator & bp) { - if(use_tableau()) { - analyze_new_bounds_on_row_tableau(i, bp); - } else { - m_mpq_lar_core_solver.calculate_pivot_row(i); - substitute_basis_var_in_terms_for_row(i); - analyze_new_bounds_on_row(i, bp); - } -} - - -linear_combination_iterator * lar_solver::create_new_iter_from_term(unsigned term_index) const { - lp_assert(false); // not implemented - return nullptr; - // new linear_combination_iterator_on_vector(m_terms[adjust_term_index(term_index)]->coeffs_as_vector()); -} - -unsigned lar_solver::adjust_column_index_to_term_index(unsigned j) const { - unsigned ext_var_or_term = m_columns_to_ext_vars_or_term_indices[j]; - return ext_var_or_term < m_terms_start_index ? j : ext_var_or_term; -} - -void lar_solver::propagate_bounds_on_a_term(const lar_term& t, bound_propagator & bp, unsigned term_offset) { - lp_assert(false); // not implemented -} - - -void lar_solver::explain_implied_bound(implied_bound & ib, bound_propagator & bp) { - unsigned i = ib.m_row_or_term_index; - int bound_sign = ib.m_is_low_bound? 1: -1; - int j_sign = (ib.m_coeff_before_j_is_pos ? 1 :-1) * bound_sign; - unsigned m_j = ib.m_j; - if (is_term(m_j)) { - auto it = m_ext_vars_to_columns.find(m_j); - lp_assert(it != m_ext_vars_to_columns.end()); - m_j = it->second.ext_j(); - } - for (auto const& r : A_r().m_rows[i]) { - unsigned j = r.m_j; - mpq const& a = r.get_val(); - if (j == m_j) continue; - if (is_term(j)) { - auto it = m_ext_vars_to_columns.find(j); - lp_assert(it != m_ext_vars_to_columns.end()); - j = it->second.ext_j(); - } - int a_sign = is_pos(a)? 1: -1; - int sign = j_sign * a_sign; - 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); - } - // lp_assert(implied_bound_is_correctly_explained(ib, explanation)); -} - - -bool lar_solver::term_is_used_as_row(unsigned term) const { - lp_assert(is_term(term)); - return contains(m_ext_vars_to_columns, term); -} - -void lar_solver::propagate_bounds_on_terms(bound_propagator & bp) { - for (unsigned i = 0; i < m_terms.size(); i++) { - if (term_is_used_as_row(i + m_terms_start_index)) - continue; // this term is used a left side of a constraint, - // it was processed as a touched row if needed - propagate_bounds_on_a_term(*m_terms[i], bp, i); - } -} - - -// goes over touched rows and tries to induce bounds -void lar_solver::propagate_bounds_for_touched_rows(bound_propagator & bp) { - if (!use_tableau()) - return; // todo: consider to remove the restriction - - for (unsigned i : m_rows_with_changed_bounds.m_index) { - calculate_implied_bounds_for_row(i, bp); - } - m_rows_with_changed_bounds.clear(); - if (!use_tableau()) { - propagate_bounds_on_terms(bp); - } -} - -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() { - lp_assert(inf_int_set_is_correct()); - m_settings.st().m_make_feasible++; - if (A_r().column_count() > m_settings.st().m_max_cols) - m_settings.st().m_max_cols = A_r().column_count(); - if (A_r().row_count() > m_settings.st().m_max_rows) - m_settings.st().m_max_rows = A_r().row_count(); - 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; - auto ret = solve(); - TRACE("intinf", m_int_solver->display_inf_or_int_inf_columns(tout);); - lp_assert(inf_int_set_is_correct()); - return ret; -} - -lp_status lar_solver::solve() { - lp_assert(inf_int_set_is_correct()); - if (m_status == lp_status::INFEASIBLE) { - return m_status; - } - solve_with_core_solver(); - if (m_status != lp_status::INFEASIBLE) { - if (m_settings.bound_propagation()) - detect_rows_with_changed_bounds(); - } - - m_columns_with_changed_bound.clear(); - lp_assert(inf_int_set_is_correct()); - return m_status; -} - -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_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())); -} - - -unsigned lar_solver::get_total_iterations() const { return m_mpq_lar_core_solver.m_r_solver.total_iterations(); } -vector lar_solver::get_list_of_all_var_indices() const { - vector ret; - for (unsigned j = 0; j < m_mpq_lar_core_solver.m_r_heading.size(); j++) - ret.push_back(j); - return ret; -} -void lar_solver::push() { - m_simplex_strategy = m_settings.simplex_strategy(); - m_simplex_strategy.push(); - m_columns_to_ul_pairs.push(); - m_infeasible_column_index.push(); - m_mpq_lar_core_solver.push(); - m_term_count = m_terms.size(); - m_term_count.push(); - m_constraint_count = m_constraints.size(); - m_constraint_count.push(); -} - -void lar_solver::clean_popped_elements(unsigned n, int_set& set) { - vector to_remove; - for (unsigned j: set.m_index) - if (j >= n) - to_remove.push_back(j); - for (unsigned j : to_remove) - set.erase(j); -} - -void lar_solver::shrink_inf_set_after_pop(unsigned n, int_set & set) { - clean_popped_elements(n, set); - set.resize(n); -} - - -void lar_solver::pop(unsigned k) { - TRACE("arith_int", tout << "pop" << std::endl;); - lp_assert(inf_int_set_is_correct()); - TRACE("lar_solver", tout << "k = " << k << std::endl;); - - int n_was = static_cast(m_ext_vars_to_columns.size()); - m_infeasible_column_index.pop(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); - TRACE("arith_int", tout << "pop" << std::endl;); - if (m_settings.use_tableau()) { - pop_tableau(); - } - lp_assert(A_r().column_count() == n); - m_columns_to_ul_pairs.pop(k); - - m_mpq_lar_core_solver.pop(k); - clean_popped_elements(n, m_columns_with_changed_bound); - clean_popped_elements(n, m_inf_int_set); - unsigned m = A_r().row_count(); - lp_assert(inf_int_set_is_correct()); - clean_popped_elements(m, m_rows_with_changed_bounds); - lp_assert(inf_int_set_is_correct()); - clean_inf_set_of_r_solver_after_pop(); - 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()); - - - lp_assert(ax_is_correct()); - lp_assert(m_mpq_lar_core_solver.m_r_solver.inf_set_is_correct()); - m_constraint_count.pop(k); - for (unsigned i = m_constraint_count; i < m_constraints.size(); i++) - delete m_constraints[i]; - - m_constraints.resize(m_constraint_count); - m_term_count.pop(k); - for (unsigned i = m_term_count; i < m_terms.size(); i++) { - delete m_terms[i]; - } - m_terms.resize(m_term_count); - m_simplex_strategy.pop(k); - m_settings.simplex_strategy() = m_simplex_strategy; - lp_assert(sizes_are_correct()); - lp_assert((!m_settings.use_tableau()) || m_mpq_lar_core_solver.m_r_solver.reduced_costs_are_correct_tableau()); - m_status = m_mpq_lar_core_solver.m_r_solver.current_x_is_feasible()? lp_status::OPTIMAL: lp_status::UNKNOWN; - -} - -vector lar_solver::get_all_constraint_indices() const { - vector ret; - constraint_index i = 0; - while ( i < m_constraints.size()) - ret.push_back(i++); - return ret; -} - -bool lar_solver::maximize_term_on_tableau(const vector> & term, - impq &term_max) { - if (settings().simplex_strategy() == simplex_strategy_enum::undecided) - decide_on_strategy_and_adjust_initial_state(); - - m_mpq_lar_core_solver.solve(); - if (m_mpq_lar_core_solver.m_r_solver.get_status() == lp_status::UNBOUNDED) - return false; - - term_max = 0; - for (auto & p : term) - term_max += p.first * m_mpq_lar_core_solver.m_r_x[p.second]; - - return true; -} - -bool lar_solver::costs_are_zeros_for_r_solver() const { - for (unsigned j = 0; j < m_mpq_lar_core_solver.m_r_solver.m_costs.size(); j++) { - lp_assert(is_zero(m_mpq_lar_core_solver.m_r_solver.m_costs[j])); - } - return true; -} -bool lar_solver::reduced_costs_are_zeroes_for_r_solver() const { - for (unsigned j = 0; j < m_mpq_lar_core_solver.m_r_solver.m_d.size(); j++) { - lp_assert(is_zero(m_mpq_lar_core_solver.m_r_solver.m_d[j])); - } - return true; -} - -void lar_solver::set_costs_to_zero(const vector> & term) { - auto & rslv = m_mpq_lar_core_solver.m_r_solver; - auto & jset = m_mpq_lar_core_solver.m_r_solver.m_inf_set; // hijack this set that should be empty right now - lp_assert(jset.m_index.size()==0); - - for (auto & p : term) { - unsigned j = p.second; - rslv.m_costs[j] = zero_of_type(); - int i = rslv.m_basis_heading[j]; - if (i < 0) - jset.insert(j); - else { - for (auto & rc : A_r().m_rows[i]) - jset.insert(rc.m_j); - } - } - - for (unsigned j : jset.m_index) - rslv.m_d[j] = zero_of_type(); - - jset.clear(); - - lp_assert(reduced_costs_are_zeroes_for_r_solver()); - lp_assert(costs_are_zeros_for_r_solver()); -} - -void lar_solver::prepare_costs_for_r_solver(const vector> & term) { - - auto & rslv = m_mpq_lar_core_solver.m_r_solver; - rslv.m_using_infeas_costs = false; - lp_assert(costs_are_zeros_for_r_solver()); - lp_assert(reduced_costs_are_zeroes_for_r_solver()); - rslv.m_costs.resize(A_r().column_count(), zero_of_type()); - for (auto & p : term) { - unsigned j = p.second; - rslv.m_costs[j] = p.first; - if (rslv.m_basis_heading[j] < 0) - rslv.m_d[j] += p.first; - else - rslv.update_reduced_cost_for_basic_column_cost_change(- p.first, j); - } - lp_assert(rslv.reduced_costs_are_correct_tableau()); -} - -bool lar_solver::maximize_term_on_corrected_r_solver(const vector> & term, - impq &term_max) { - settings().backup_costs = false; - switch (settings().simplex_strategy()) { - case simplex_strategy_enum::tableau_rows: - prepare_costs_for_r_solver(term); - settings().simplex_strategy() = simplex_strategy_enum::tableau_costs; - { - 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(lp_status::OPTIMAL); - return ret; - } - case simplex_strategy_enum::tableau_costs: - prepare_costs_for_r_solver(term); - { - bool ret = maximize_term_on_tableau(term, term_max); - set_costs_to_zero(term); - m_mpq_lar_core_solver.m_r_solver.set_status(lp_status::OPTIMAL); - return ret; - } - - case simplex_strategy_enum::lu: - lp_assert(false); // not implemented - return false; - default: - lp_unreachable(); // wrong mode - } - return false; -} -// starting from a given feasible state look for the maximum of the term -// return true if found and false if unbounded -bool lar_solver::maximize_term(const vector> & term, - impq &term_max) { - lp_assert(m_mpq_lar_core_solver.m_r_solver.current_x_is_feasible()); - m_mpq_lar_core_solver.m_r_solver.m_look_for_feasible_solution_only = false; - return maximize_term_on_corrected_r_solver(term, term_max); -} - - - -const lar_term & lar_solver::get_term(unsigned j) const { - lp_assert(j >= m_terms_start_index); - return *m_terms[j - m_terms_start_index]; -} - -void lar_solver::pop_core_solver_params() { - pop_core_solver_params(1); -} - -void lar_solver::pop_core_solver_params(unsigned k) { - A_r().pop(k); - A_d().pop(k); -} - - -void lar_solver::set_upper_bound_witness(var_index j, constraint_index ci) { - ul_pair ul = m_columns_to_ul_pairs[j]; - ul.upper_bound_witness() = ci; - m_columns_to_ul_pairs[j] = ul; -} - -void lar_solver::set_low_bound_witness(var_index j, constraint_index ci) { - ul_pair ul = m_columns_to_ul_pairs[j]; - ul.low_bound_witness() = ci; - m_columns_to_ul_pairs[j] = ul; -} - -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; - } else { - it->second += a; - } -} - - -void lar_solver::substitute_terms_in_linear_expression(const vector>& left_side_with_terms, - 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_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_monoid_in_map(coeffs, t.first * p.second , p.first); - } - free_coeff += t.first * term.m_v; - } - } - - for (auto & p : coeffs) - left_side.push_back(std::make_pair(p.second, p.first)); -} - - -void lar_solver::detect_rows_of_bound_change_column_for_nbasic_column(unsigned j) { - if (A_r().row_count() != m_column_buffer.data_size()) - m_column_buffer.resize(A_r().row_count()); - else - m_column_buffer.clear(); - lp_assert(m_column_buffer.size() == 0 && m_column_buffer.is_OK()); - - m_mpq_lar_core_solver.m_r_solver.solve_Bd(j, m_column_buffer); - for (unsigned i : m_column_buffer.m_index) - m_rows_with_changed_bounds.insert(i); -} - - - -void lar_solver::detect_rows_of_bound_change_column_for_nbasic_column_tableau(unsigned j) { - for (auto & rc : m_mpq_lar_core_solver.m_r_A.m_columns[j]) - m_rows_with_changed_bounds.insert(rc.m_i); -} - -bool lar_solver::use_tableau() const { return m_settings.use_tableau(); } - -bool lar_solver::use_tableau_costs() const { - return m_settings.simplex_strategy() == simplex_strategy_enum::tableau_costs; -} - -void lar_solver::detect_rows_of_column_with_bound_change(unsigned j) { - if (m_mpq_lar_core_solver.m_r_heading[j] >= 0) { // it is a basic column - // just mark the row at touched and exit - m_rows_with_changed_bounds.insert(m_mpq_lar_core_solver.m_r_heading[j]); - return; - } - - if (use_tableau()) - detect_rows_of_bound_change_column_for_nbasic_column_tableau(j); - else - detect_rows_of_bound_change_column_for_nbasic_column(j); -} - -void lar_solver::adjust_x_of_column(unsigned j) { - lp_assert(false); -} - -bool lar_solver::row_is_correct(unsigned i) const { - numeric_pair r = zero_of_type>(); - for (const auto & c : A_r().m_rows[i]) - r += c.m_value * m_mpq_lar_core_solver.m_r_x[c.m_j]; - return is_zero(r); -} - -bool lar_solver::ax_is_correct() const { - for (unsigned i = 0; i < A_r().row_count(); i++) { - if (!row_is_correct(i)) - return false; - } - return true; -} - -bool lar_solver::tableau_with_costs() const { - return m_settings.simplex_strategy() == simplex_strategy_enum::tableau_costs; -} - -bool lar_solver::costs_are_used() const { - return m_settings.simplex_strategy() != simplex_strategy_enum::tableau_rows; -} - -void lar_solver::change_basic_columns_dependend_on_a_given_nb_column(unsigned j, const numeric_pair & delta) { - lp_assert(inf_int_set_is_correct()); - if (use_tableau()) { - for (const auto & c : A_r().m_columns[j]) { - unsigned bj = m_mpq_lar_core_solver.m_r_basis[c.m_i]; - if (tableau_with_costs()) { - m_basic_columns_with_changed_cost.insert(bj); - } - m_mpq_lar_core_solver.m_r_solver.update_x_with_delta_and_track_feasibility(bj, - A_r().get_val(c) * delta); - TRACE("change_x_del", - tout << "changed basis column " << bj << ", it is " << - ( m_mpq_lar_core_solver.m_r_solver.column_is_feasible(bj)? "feas":"inf") << std::endl;); - - } - } else { - m_column_buffer.clear(); - m_column_buffer.resize(A_r().row_count()); - m_mpq_lar_core_solver.m_r_solver.solve_Bd(j, m_column_buffer); - for (unsigned i : m_column_buffer.m_index) { - unsigned bj = m_mpq_lar_core_solver.m_r_basis[i]; - m_mpq_lar_core_solver.m_r_solver.update_x_with_delta_and_track_feasibility(bj, -m_column_buffer[i] * delta); - } - } - lp_assert(inf_int_set_is_correct()); -} - -void lar_solver::update_x_and_inf_costs_for_column_with_changed_bounds(unsigned j) { - if (m_mpq_lar_core_solver.m_r_heading[j] >= 0) { - if (costs_are_used()) { - bool was_infeas = m_mpq_lar_core_solver.m_r_solver.m_inf_set.contains(j); - m_mpq_lar_core_solver.m_r_solver.track_column_feasibility(j); - if (was_infeas != m_mpq_lar_core_solver.m_r_solver.m_inf_set.contains(j)) - m_basic_columns_with_changed_cost.insert(j); - } else { - m_mpq_lar_core_solver.m_r_solver.track_column_feasibility(j); - } - } else { - numeric_pair delta; - if (m_mpq_lar_core_solver.m_r_solver.make_column_feasible(j, delta)) - change_basic_columns_dependend_on_a_given_nb_column(j, delta); - } -} - - -void lar_solver::detect_rows_with_changed_bounds_for_column(unsigned j) { - if (m_mpq_lar_core_solver.m_r_heading[j] >= 0) { - m_rows_with_changed_bounds.insert(m_mpq_lar_core_solver.m_r_heading[j]); - return; - } - - if (use_tableau()) - detect_rows_of_bound_change_column_for_nbasic_column_tableau(j); - else - detect_rows_of_bound_change_column_for_nbasic_column(j); -} - -void lar_solver::detect_rows_with_changed_bounds() { - for (auto j : m_columns_with_changed_bound.m_index) - detect_rows_with_changed_bounds_for_column(j); -} - -void lar_solver::update_x_and_inf_costs_for_columns_with_changed_bounds() { - for (auto j : m_columns_with_changed_bound.m_index) - update_x_and_inf_costs_for_column_with_changed_bounds(j); -} - -void lar_solver::update_x_and_inf_costs_for_columns_with_changed_bounds_tableau() { - for (auto j : m_columns_with_changed_bound.m_index) - update_x_and_inf_costs_for_column_with_changed_bounds(j); - - if (tableau_with_costs()) { - for (unsigned j : m_basic_columns_with_changed_cost.m_index) - m_mpq_lar_core_solver.m_r_solver.update_inf_cost_for_column_tableau(j); - lp_assert(m_mpq_lar_core_solver.m_r_solver.reduced_costs_are_correct_tableau()); - } -} - - -void lar_solver::solve_with_core_solver() { - if (!use_tableau()) - add_last_rows_to_lu(m_mpq_lar_core_solver.m_r_solver); - if (m_mpq_lar_core_solver.need_to_presolve_with_double_solver()) { - add_last_rows_to_lu(m_mpq_lar_core_solver.m_d_solver); - } - m_mpq_lar_core_solver.prefix_r(); - if (costs_are_used()) { - m_basic_columns_with_changed_cost.clear(); - m_basic_columns_with_changed_cost.resize(m_mpq_lar_core_solver.m_r_x.size()); - } - if (use_tableau()) - update_x_and_inf_costs_for_columns_with_changed_bounds_tableau(); - else - update_x_and_inf_costs_for_columns_with_changed_bounds(); - TRACE("intinf", m_int_solver->display_inf_or_int_inf_columns(tout);); - m_mpq_lar_core_solver.solve(); - set_status(m_mpq_lar_core_solver.m_r_solver.get_status()); - lp_assert(inf_int_set_is_correct()); - lp_assert(m_status != lp_status::OPTIMAL || all_constraints_hold()); -} - - -numeric_pair lar_solver::get_basic_var_value_from_row_directly(unsigned i) { - numeric_pair r = zero_of_type>(); - - unsigned bj = m_mpq_lar_core_solver.m_r_solver.m_basis[i]; - for (const auto & c: A_r().m_rows[i]) { - if (c.m_j == bj) continue; - const auto & x = m_mpq_lar_core_solver.m_r_x[c.m_j]; - if (!is_zero(x)) - r -= c.m_value * x; - } - return r; -} - -numeric_pair lar_solver::get_basic_var_value_from_row(unsigned i) { - if (settings().use_tableau()) { - return get_basic_var_value_from_row_directly(i); - } - - numeric_pair r = zero_of_type>(); - m_mpq_lar_core_solver.calculate_pivot_row(i); - for (unsigned j : m_mpq_lar_core_solver.m_r_solver.m_pivot_row.m_index) { - lp_assert(m_mpq_lar_core_solver.m_r_solver.m_basis_heading[j] < 0); - r -= m_mpq_lar_core_solver.m_r_solver.m_pivot_row.m_data[j] * m_mpq_lar_core_solver.m_r_x[j]; - } - return r; -} - -template -void lar_solver::add_last_rows_to_lu(lp_primal_core_solver & s) { - auto & f = s.m_factorization; - if (f != nullptr) { - auto columns_to_replace = f->get_set_of_columns_to_replace_for_add_last_rows(s.m_basis_heading); - if (f->m_refactor_counter + columns_to_replace.size() >= 200 || f->has_dense_submatrix()) { - delete f; - f = nullptr; - } else { - f->add_last_rows_to_B(s.m_basis_heading, columns_to_replace); - } - } - if (f == nullptr) { - init_factorization(f, s.m_A, s.m_basis, m_settings); - if (f->get_status() != LU_status::OK) { - delete f; - f = nullptr; - } - } - -} - -bool lar_solver::x_is_correct() const { - if (m_mpq_lar_core_solver.m_r_x.size() != A_r().column_count()) { - // std::cout << "the size is off " << m_r_solver.m_x.size() << ", " << A().column_count() << std::endl; - return false; - } - for (unsigned i = 0; i < A_r().row_count(); i++) { - numeric_pair delta = A_r().dot_product_with_row(i, m_mpq_lar_core_solver.m_r_x); - if (!delta.is_zero()) { - // std::cout << "x is off ("; - // std::cout << "m_b[" << i << "] = " << m_b[i] << " "; - // std::cout << "left side = " << A().dot_product_with_row(i, m_r_solver.m_x) << ' '; - // std::cout << "delta = " << delta << ' '; - // std::cout << "iters = " << total_iterations() << ")" << std::endl; - // std::cout << "row " << i << " is off" << std::endl; - return false; - } - } - return true;; - -} - -bool lar_solver::var_is_registered(var_index vj) const { - if (vj >= m_terms_start_index) { - if (vj - m_terms_start_index >= m_terms.size()) - return false; - } else if ( vj >= A_r().column_count()) { - return false; - } - return true; -} - -unsigned lar_solver::constraint_stack_size() const { - return m_constraint_count.stack_size(); -} - -void lar_solver::fill_last_row_of_A_r(static_matrix> & A, const lar_term * ls) { - lp_assert(A.row_count() > 0); - lp_assert(A.column_count() > 0); - unsigned last_row = A.row_count() - 1; - lp_assert(A.m_rows[last_row].size() == 0); - for (auto & t : ls->m_coeffs) { - lp_assert(!is_zero(t.second)); - var_index j = t.first; - A.set(last_row, j, - t.second); - } - unsigned basis_j = A.column_count() - 1; - A.set(last_row, basis_j, mpq(1)); -} - -template -void lar_solver::create_matrix_A(static_matrix & matr) { - lp_assert(false); // not implemented - /* - unsigned m = number_or_nontrivial_left_sides(); - unsigned n = m_vec_of_canonic_left_sides.size(); - if (matr.row_count() == m && matr.column_count() == n) - return; - matr.init_empty_matrix(m, n); - copy_from_mpq_matrix(matr); - */ -} - -template -void lar_solver::copy_from_mpq_matrix(static_matrix & matr) { - matr.m_rows.resize(A_r().row_count()); - matr.m_columns.resize(A_r().column_count()); - for (unsigned i = 0; i < matr.row_count(); i++) { - for (auto & it : A_r().m_rows[i]) { - matr.set(i, it.m_j, convert_struct::convert(it.get_val())); - } - } -} - - -bool lar_solver::try_to_set_fixed(column_info & ci) { - if (ci.upper_bound_is_set() && ci.low_bound_is_set() && ci.get_upper_bound() == ci.get_low_bound() && !ci.is_fixed()) { - ci.set_fixed_value(ci.get_upper_bound()); - return true; - } - return false; -} - -column_type lar_solver::get_column_type(const column_info & ci) { - auto ret = ci.get_column_type_no_flipping(); - if (ret == column_type::boxed) { // changing boxed to fixed because of the no span - if (ci.get_low_bound() == ci.get_upper_bound()) - ret = column_type::fixed; - } - return ret; -} - -std::string lar_solver::get_column_name(unsigned j) const { - if (j >= m_terms_start_index) - return std::string("_t") + T_to_string(j); - if (j >= m_columns_to_ext_vars_or_term_indices.size()) - return std::string("_s") + T_to_string(j); - - return std::string("v") + T_to_string(m_columns_to_ext_vars_or_term_indices[j]); -} - -bool lar_solver::all_constrained_variables_are_registered(const vector>& left_side) { - for (auto it : left_side) { - if (! var_is_registered(it.second)) - return false; - } - return true; -} - -bool lar_solver::all_constraints_hold() const { - if (m_settings.get_cancel_flag()) - return true; - std::unordered_map var_map; - get_model_do_not_care_about_diff_vars(var_map); - - for (unsigned i = 0; i < m_constraints.size(); i++) { - if (!constraint_holds(*m_constraints[i], var_map)) { - print_constraint(i, std::cout); - return false; - } - } - return true; -} - -bool lar_solver::constraint_holds(const lar_base_constraint & constr, std::unordered_map & var_map) const { - return true; - mpq left_side_val = get_left_side_val(constr, var_map); - switch (constr.m_kind) { - case LE: return left_side_val <= constr.m_right_side; - case LT: return left_side_val < constr.m_right_side; - case GE: return left_side_val >= constr.m_right_side; - case GT: return left_side_val > constr.m_right_side; - case EQ: return left_side_val == constr.m_right_side; - default: - lp_unreachable(); - } - return false; // it is unreachable -} - -bool lar_solver::the_relations_are_of_same_type(const vector> & evidence, lconstraint_kind & the_kind_of_sum) const { - unsigned n_of_G = 0, n_of_L = 0; - bool strict = false; - for (auto & it : evidence) { - mpq coeff = it.first; - constraint_index con_ind = it.second; - lconstraint_kind kind = coeff.is_pos() ? - m_constraints[con_ind]->m_kind : - flip_kind(m_constraints[con_ind]->m_kind); - if (kind == GT || kind == LT) - strict = true; - if (kind == GE || kind == GT) n_of_G++; - else if (kind == LE || kind == LT) n_of_L++; - } - the_kind_of_sum = n_of_G ? GE : (n_of_L ? LE : EQ); - if (strict) - the_kind_of_sum = static_cast((static_cast(the_kind_of_sum) / 2)); - - return n_of_G == 0 || n_of_L == 0; -} - -void lar_solver::register_in_map(std::unordered_map & coeffs, const lar_base_constraint & cn, const mpq & a) { - for (auto & it : cn.get_left_side_coefficients()) { - unsigned j = it.second; - auto p = coeffs.find(j); - if (p == coeffs.end()) - coeffs[j] = it.first * a; - else { - p->second += it.first * a; - if (p->second.is_zero()) - coeffs.erase(p); - } - } -} - -bool lar_solver::the_left_sides_sum_to_zero(const vector> & evidence) const { - std::unordered_map coeff_map; - for (auto & it : evidence) { - mpq coeff = it.first; - constraint_index con_ind = it.second; - lp_assert(con_ind < m_constraints.size()); - register_in_map(coeff_map, *m_constraints[con_ind], coeff); - } - - if (!coeff_map.empty()) { - std::cout << "left side = "; - vector> t; - for (auto & it : coeff_map) { - t.push_back(std::make_pair(it.second, it.first)); - } - print_linear_combination_of_column_indices(t, std::cout); - std::cout << std::endl; - return false; - } - - return true; -} - -bool lar_solver::the_right_sides_do_not_sum_to_zero(const vector> & evidence) { - mpq ret = numeric_traits::zero(); - for (auto & it : evidence) { - mpq coeff = it.first; - constraint_index con_ind = it.second; - lp_assert(con_ind < m_constraints.size()); - const lar_constraint & constr = *m_constraints[con_ind]; - ret += constr.m_right_side * coeff; - } - return !numeric_traits::is_zero(ret); -} - -bool lar_solver::explanation_is_correct(const vector>& explanation) const { -#ifdef LEAN_DEBUG - lconstraint_kind kind; - lp_assert(the_relations_are_of_same_type(explanation, kind)); - lp_assert(the_left_sides_sum_to_zero(explanation)); - mpq rs = sum_of_right_sides_of_explanation(explanation); - switch (kind) { - case LE: lp_assert(rs < zero_of_type()); - break; - case LT: lp_assert(rs <= zero_of_type()); - break; - case GE: lp_assert(rs > zero_of_type()); - break; - case GT: lp_assert(rs >= zero_of_type()); - break; - case EQ: lp_assert(rs != zero_of_type()); - break; - default: - lp_assert(false); - return false; - } -#endif - return true; -} - -bool lar_solver::inf_explanation_is_correct() const { -#ifdef LEAN_DEBUG - vector> explanation; - get_infeasibility_explanation(explanation); - return explanation_is_correct(explanation); -#endif - return true; -} - -mpq lar_solver::sum_of_right_sides_of_explanation(const vector> & explanation) const { - mpq ret = numeric_traits::zero(); - for (auto & it : explanation) { - mpq coeff = it.first; - constraint_index con_ind = it.second; - lp_assert(con_ind < m_constraints.size()); - ret += (m_constraints[con_ind]->m_right_side - m_constraints[con_ind]->get_free_coeff_of_left_side()) * coeff; - } - return ret; -} - -bool lar_solver::has_lower_bound(var_index var, constraint_index& ci, mpq& value, bool& is_strict) { - - 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_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]; - value = p.x; - is_strict = p.y.is_pos(); - return true; - } - else { - return false; - } -} - -bool lar_solver::has_upper_bound(var_index var, constraint_index& ci, mpq& value, bool& is_strict) { - - 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_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]; - value = p.x; - is_strict = p.y.is_neg(); - return true; - } - else { - return false; - } -} - -void lar_solver::get_infeasibility_explanation(vector> & explanation) const { - explanation.clear(); - if (m_infeasible_column_index != -1) { - fill_explanation_from_infeasible_column(explanation); - return; - } - if (m_mpq_lar_core_solver.get_infeasible_sum_sign() == 0) { - return; - } - // the infeasibility sign - int inf_sign; - auto inf_row = m_mpq_lar_core_solver.get_infeasibility_info(inf_sign); - get_infeasibility_explanation_for_inf_sign(explanation, inf_row, inf_sign); - lp_assert(explanation_is_correct(explanation)); -} - - - -void lar_solver::get_infeasibility_explanation_for_inf_sign( - vector> & explanation, - const vector> & inf_row, - int inf_sign) const { - - for (auto & it : inf_row) { - mpq coeff = it.first; - unsigned j = it.second; - - int adj_sign = coeff.is_pos() ? inf_sign : -inf_sign; - 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()); - explanation.push_back(std::make_pair(coeff, bound_constr_i)); - } -} - -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 == lp_status::OPTIMAL); - unsigned i; - do { - // different pairs have to produce different singleton values - std::unordered_set set_of_different_pairs; - std::unordered_set set_of_different_singles; - delta = m_mpq_lar_core_solver.find_delta_for_strict_bounds(delta); - for (i = 0; i < m_mpq_lar_core_solver.m_r_x.size(); i++ ) { - const numeric_pair & rp = m_mpq_lar_core_solver.m_r_x[i]; - set_of_different_pairs.insert(rp); - mpq x = rp.x + delta * rp.y; - set_of_different_singles.insert(x); - if (set_of_different_pairs.size() - != set_of_different_singles.size()) { - delta /= mpq(2); - break; - } - - variable_values[i] = x; - } - } while (i != m_mpq_lar_core_solver.m_r_x.size()); -} - -void lar_solver::get_model_do_not_care_about_diff_vars(std::unordered_map & variable_values) const { - mpq delta = mpq(1); - delta = m_mpq_lar_core_solver.find_delta_for_strict_bounds(delta); - for (unsigned i = 0; i < m_mpq_lar_core_solver.m_r_x.size(); i++ ) { - const impq & rp = m_mpq_lar_core_solver.m_r_x[i]; - variable_values[i] = rp.x + delta * rp.y; - } -} - - -std::string lar_solver::get_variable_name(var_index vi) const { - return get_column_name(vi); -} - -// ********** print region start -void lar_solver::print_constraint(constraint_index ci, std::ostream & out) const { - if (ci >= m_constraints.size()) { - out << "constraint " << T_to_string(ci) << " is not found"; - out << std::endl; - return; - } - - print_constraint(m_constraints[ci], out); -} - -void lar_solver::print_constraints(std::ostream& out) const { - for (auto c : m_constraints) { - print_constraint(c, out); - } -} - -void lar_solver::print_terms(std::ostream& out) const { - for (auto it : m_terms) { - print_term(*it, out); - out << "\n"; - } -} - -void lar_solver::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 lar_solver::print_term(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(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()) { - var_index j = it.second; - auto vi = var_map.find(j); - lp_assert(vi != var_map.end()); - ret += it.first * vi->second; - } - return ret; -} - -void lar_solver::print_constraint(const lar_base_constraint * c, std::ostream & out) const { - print_left_side_of_constraint(c, out); - out << " " << lconstraint_kind_string(c->m_kind) << " " << c->m_right_side << std::endl; -} - -void lar_solver::fill_var_set_for_random_update(unsigned sz, var_index const * vars, vector& column_list) { - for (unsigned i = 0; i < sz; i++) { - var_index var = vars[i]; - if (var >= m_terms_start_index) { // handle the term - for (auto & it : m_terms[var - m_terms_start_index]->m_coeffs) { - column_list.push_back(it.first); - } - } else { - column_list.push_back(var); - } - } -} - -void lar_solver::random_update(unsigned sz, var_index const * vars) { - vector column_list; - fill_var_set_for_random_update(sz, vars, column_list); - random_updater ru(*this, column_list); - ru.update(); - lp_assert(inf_int_set_is_correct()); -} - - -void lar_solver::pivot_fixed_vars_from_basis() { - m_mpq_lar_core_solver.m_r_solver.pivot_fixed_vars_from_basis(); -} - -void lar_solver::pop() { - pop(1); -} - -bool lar_solver::column_represents_row_in_tableau(unsigned j) { - 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) { - // i, j - is the indices of the bottom-right element of the tableau - lp_assert(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; - } - - lp_assert(non_zero_column_cell_index != -1); - lp_assert(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 lar_solver::remove_last_row_and_column_from_tableau(unsigned j) { - lp_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_costs.size()); - auto & slv = m_mpq_lar_core_solver.m_r_solver; - unsigned i = A_r().row_count() - 1; //last row index - make_sure_that_the_bottom_right_elem_not_zero_in_tableau(i, j); - if (slv.m_basis_heading[j] < 0) { - slv.pivot_column_tableau(j, i); - } - - auto & last_row = A_r().m_rows[i]; - mpq &cost_j = m_mpq_lar_core_solver.m_r_solver.m_costs[j]; - bool cost_is_nz = !is_zero(cost_j); - for (unsigned k = last_row.size(); k-- > 0;) { - auto &rc = last_row[k]; - if (cost_is_nz) { - m_mpq_lar_core_solver.m_r_solver.m_d[rc.m_j] += cost_j*rc.get_val(); - } - - A_r().remove_element(last_row, rc); - } - lp_assert(last_row.size() == 0); - lp_assert(A_r().m_columns[j].size() == 0); - A_r().m_rows.pop_back(); - A_r().m_columns.pop_back(); - slv.m_b.pop_back(); -} - -void lar_solver::remove_last_column_from_A() { - // the last column has to be empty - lp_assert(A_r().m_columns.back().size() == 0); - A_r().m_columns.pop_back(); -} - -void lar_solver::remove_last_column_from_basis_tableau(unsigned j) { - auto& rslv = m_mpq_lar_core_solver.m_r_solver; - int i = rslv.m_basis_heading[j]; - if (i >= 0) { // j is a basic var - int last_pos = static_cast(rslv.m_basis.size()) - 1; - lp_assert(last_pos >= 0); - if (i != last_pos) { - unsigned j_at_last_pos = rslv.m_basis[last_pos]; - rslv.m_basis[i] = j_at_last_pos; - rslv.m_basis_heading[j_at_last_pos] = i; - } - rslv.m_basis.pop_back(); // remove j from the basis - } else { - int last_pos = static_cast(rslv.m_nbasis.size()) - 1; - lp_assert(last_pos >= 0); - i = - 1 - i; - if (i != last_pos) { - unsigned j_at_last_pos = rslv.m_nbasis[last_pos]; - rslv.m_nbasis[i] = j_at_last_pos; - rslv.m_basis_heading[j_at_last_pos] = - i - 1; - } - rslv.m_nbasis.pop_back(); // remove j from the basis - } - rslv.m_basis_heading.pop_back(); - lp_assert(rslv.m_basis.size() == A_r().row_count()); - lp_assert(rslv.basis_heading_is_correct()); -} - -void lar_solver::remove_last_column_from_tableau() { - auto& rslv = m_mpq_lar_core_solver.m_r_solver; - unsigned j = A_r().column_count() - 1; - lp_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_costs.size()); - if (column_represents_row_in_tableau(j)) { - remove_last_row_and_column_from_tableau(j); - if (rslv.m_basis_heading[j] < 0) - rslv.change_basis_unconditionally(j, rslv.m_basis[A_r().row_count()]); // A_r().row_count() is the index of the last row in the basis still - } - else { - remove_last_column_from_A(); - } - rslv.m_x.pop_back(); - rslv.m_d.pop_back(); - rslv.m_costs.pop_back(); - - remove_last_column_from_basis_tableau(j); - lp_assert(m_mpq_lar_core_solver.r_basis_is_OK()); - lp_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_costs.size()); -} - -void lar_solver::pop_tableau() { - lp_assert(m_mpq_lar_core_solver.m_r_solver.m_costs.size() == A_r().column_count()); - - lp_assert(m_mpq_lar_core_solver.m_r_solver.m_basis.size() == A_r().row_count()); - lp_assert(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 - while (A_r().column_count() > m_columns_to_ext_vars_or_term_indices.size()) - remove_last_column_from_tableau(); - lp_assert(m_mpq_lar_core_solver.m_r_solver.m_costs.size() == A_r().column_count()); - lp_assert(m_mpq_lar_core_solver.m_r_solver.m_basis.size() == A_r().row_count()); - lp_assert(m_mpq_lar_core_solver.m_r_solver.basis_heading_is_correct()); -} - -void lar_solver::clean_inf_set_of_r_solver_after_pop() { - lp_assert(inf_int_set_is_correct()); - vector became_feas; - 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) { - if (m_mpq_lar_core_solver.m_r_heading[j] >= 0) { - continue; - } - // some basic columns might become non-basic - these columns need to be made feasible - numeric_pair delta; - if (m_mpq_lar_core_solver.m_r_solver.make_column_feasible(j, delta)) { - - change_basic_columns_dependend_on_a_given_nb_column(j, delta); - } - became_feas.push_back(j); - } - - for (unsigned j : became_feas) { - lp_assert(m_mpq_lar_core_solver.m_r_solver.m_basis_heading[j] < 0); - m_mpq_lar_core_solver.m_r_solver.m_d[j] -= m_mpq_lar_core_solver.m_r_solver.m_costs[j]; - m_mpq_lar_core_solver.m_r_solver.m_costs[j] = zero_of_type(); - m_mpq_lar_core_solver.m_r_solver.m_inf_set.erase(j); - } - became_feas.clear(); - for (unsigned j : m_mpq_lar_core_solver.m_r_solver.m_inf_set.m_index) { - lp_assert(m_mpq_lar_core_solver.m_r_heading[j] >= 0); - if (m_mpq_lar_core_solver.m_r_solver.column_is_feasible(j)) - became_feas.push_back(j); - } - for (unsigned j : became_feas) - m_mpq_lar_core_solver.m_r_solver.m_inf_set.erase(j); - - - if (use_tableau_costs()) { - for (unsigned j : became_feas) - m_mpq_lar_core_solver.m_r_solver.update_inf_cost_for_column_tableau(j); - for (unsigned j : basic_columns_with_changed_cost) - m_mpq_lar_core_solver.m_r_solver.update_inf_cost_for_column_tableau(j); - lp_assert(m_mpq_lar_core_solver.m_r_solver.reduced_costs_are_correct_tableau()); - } -} - -void lar_solver::shrink_explanation_to_minimum(vector> & explanation) const { - // implementing quickXplain - quick_xplain::run(explanation, *this); - lp_assert(this->explanation_is_correct(explanation)); -} - -bool lar_solver::model_is_int_feasible() const { - unsigned n = A_r().column_count(); - for (unsigned j = 0; j < n; j++) { - if (column_is_int(j) && !column_value_is_integer(j)) - return false; - } - return true; -} - -bool lar_solver::term_is_int(const lar_term * t) const { - for (auto const & p : t->m_coeffs) - if (! (column_is_int(p.first) && p.second.is_int())) - return false; - return t->m_v.is_int(); -} - -bool lar_solver::var_is_int(var_index v) const { - if (is_term(v)) { - lar_term const& t = get_term(v); - return term_is_int(&t); - } - else { - return column_is_int(v); - } -} - -bool lar_solver::column_is_int(unsigned j) const { - unsigned ext_var = m_columns_to_ext_vars_or_term_indices[j]; - lp_assert(contains(m_ext_vars_to_columns, ext_var)); - return m_ext_vars_to_columns.find(ext_var)->second.is_integer(); -} - -bool lar_solver::column_is_fixed(unsigned j) const { - return m_mpq_lar_core_solver.column_is_fixed(j); -} - - -bool lar_solver::ext_var_is_int(var_index ext_var) const { - auto it = m_ext_vars_to_columns.find(ext_var); - lp_assert(it != m_ext_vars_to_columns.end()); - return it == m_ext_vars_to_columns.end() || it->second.is_integer(); -} - -// below is the initialization functionality of lar_solver - -bool lar_solver::strategy_is_undecided() const { - return m_settings.simplex_strategy() == simplex_strategy_enum::undecided; -} - -var_index lar_solver::add_var(unsigned ext_j, bool is_int) { - TRACE("add_var", tout << "adding var " << ext_j << (is_int? " int" : " nonint") << std::endl;); - var_index i; - lp_assert(ext_j < m_terms_start_index); - - if (ext_j >= m_terms_start_index) - throw 0; // todo : what is the right way to exit? - auto it = m_ext_vars_to_columns.find(ext_j); - if (it != m_ext_vars_to_columns.end()) { - return it->second.ext_j(); - } - lp_assert(m_columns_to_ul_pairs.size() == A_r().column_count()); - i = A_r().column_count(); - 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); - } - lp_assert(inf_int_set_is_correct()); - return i; -} - -void lar_solver::register_new_ext_var_index(unsigned ext_v, bool is_int) { - lp_assert(!contains(m_ext_vars_to_columns, ext_v)); - unsigned j = static_cast(m_ext_vars_to_columns.size()); - m_ext_vars_to_columns.insert(std::make_pair(ext_v, ext_var_info(j, is_int))); - lp_assert(m_columns_to_ext_vars_or_term_indices.size() == j); - m_columns_to_ext_vars_or_term_indices.push_back(ext_v); -} - -void lar_solver::add_non_basic_var_to_core_fields(unsigned ext_j, bool is_int) { - register_new_ext_var_index(ext_j, is_int); - m_mpq_lar_core_solver.m_column_types.push_back(column_type::free_column); - m_columns_with_changed_bound.increase_size_by_one(); - m_inf_int_set.increase_size_by_one(); - add_new_var_to_core_fields_for_mpq(false); - if (use_lu()) - add_new_var_to_core_fields_for_doubles(false); -} - -void lar_solver::add_new_var_to_core_fields_for_doubles(bool register_in_basis) { - unsigned j = A_d().column_count(); - A_d().add_column(); - lp_assert(m_mpq_lar_core_solver.m_d_x.size() == j); - // lp_assert(m_mpq_lar_core_solver.m_d_low_bounds.size() == j && m_mpq_lar_core_solver.m_d_upper_bounds.size() == j); // restore later - m_mpq_lar_core_solver.m_d_x.resize(j + 1); - m_mpq_lar_core_solver.m_d_low_bounds.resize(j + 1); - m_mpq_lar_core_solver.m_d_upper_bounds.resize(j + 1); - lp_assert(m_mpq_lar_core_solver.m_d_heading.size() == j); // as A().column_count() on the entry to the method - if (register_in_basis) { - A_d().add_row(); - m_mpq_lar_core_solver.m_d_heading.push_back(m_mpq_lar_core_solver.m_d_basis.size()); - m_mpq_lar_core_solver.m_d_basis.push_back(j); - } - else { - m_mpq_lar_core_solver.m_d_heading.push_back(-static_cast(m_mpq_lar_core_solver.m_d_nbasis.size()) - 1); - m_mpq_lar_core_solver.m_d_nbasis.push_back(j); - } -} - -void lar_solver::add_new_var_to_core_fields_for_mpq(bool register_in_basis) { - unsigned j = A_r().column_count(); - A_r().add_column(); - lp_assert(m_mpq_lar_core_solver.m_r_x.size() == j); - // lp_assert(m_mpq_lar_core_solver.m_r_low_bounds.size() == j && m_mpq_lar_core_solver.m_r_upper_bounds.size() == j); // restore later - m_mpq_lar_core_solver.m_r_x.resize(j + 1); - m_mpq_lar_core_solver.m_r_low_bounds.increase_size_by_one(); - m_mpq_lar_core_solver.m_r_upper_bounds.increase_size_by_one(); - m_mpq_lar_core_solver.m_r_solver.m_inf_set.increase_size_by_one(); - m_mpq_lar_core_solver.m_r_solver.m_costs.resize(j + 1); - m_mpq_lar_core_solver.m_r_solver.m_d.resize(j + 1); - lp_assert(m_mpq_lar_core_solver.m_r_heading.size() == j); // as A().column_count() on the entry to the method - if (register_in_basis) { - A_r().add_row(); - m_mpq_lar_core_solver.m_r_heading.push_back(m_mpq_lar_core_solver.m_r_basis.size()); - m_mpq_lar_core_solver.m_r_basis.push_back(j); - if (m_settings.bound_propagation()) - m_rows_with_changed_bounds.insert(A_r().row_count() - 1); - } - else { - m_mpq_lar_core_solver.m_r_heading.push_back(-static_cast(m_mpq_lar_core_solver.m_r_nbasis.size()) - 1); - m_mpq_lar_core_solver.m_r_nbasis.push_back(j); - } -} - - -var_index lar_solver::add_term_undecided(const vector> & coeffs, - const mpq &m_v) { - m_terms.push_back(new lar_term(coeffs, m_v)); - return m_terms_start_index + m_terms.size() - 1; -} - -// terms -var_index lar_solver::add_term(const vector> & coeffs, - const mpq &m_v) { - if (strategy_is_undecided()) - return add_term_undecided(coeffs, m_v); - - m_terms.push_back(new lar_term(coeffs, m_v)); - unsigned adjusted_term_index = m_terms.size() - 1; - var_index ret = m_terms_start_index + adjusted_term_index; - if (use_tableau() && !coeffs.empty()) { - add_row_from_term_no_constraint(m_terms.back(), ret); - if (m_settings.bound_propagation()) - m_rows_with_changed_bounds.insert(A_r().row_count() - 1); - } - lp_assert(m_ext_vars_to_columns.size() == A_r().column_count()); - return ret; -} - -void lar_solver::add_row_from_term_no_constraint(const lar_term * term, unsigned term_ext_index) { - register_new_ext_var_index(term_ext_index, term_is_int(term)); - // j will be a new variable - unsigned j = A_r().column_count(); - ul_pair ul(j); - 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); - A_r().fill_last_row_with_pivoting(it, - m_mpq_lar_core_solver.m_r_solver.m_basis_heading); - m_mpq_lar_core_solver.m_r_solver.m_b.resize(A_r().column_count(), zero_of_type()); - } - else { - fill_last_row_of_A_r(A_r(), term); - } - m_mpq_lar_core_solver.m_r_solver.update_x_and_call_tracker(j, get_basic_var_value_from_row_directly(A_r().row_count() - 1)); - if (use_lu()) - fill_last_row_of_A_d(A_d(), term); - lp_assert(inf_int_set_is_correct()); -} - -void lar_solver::add_basic_var_to_core_fields() { - bool use_lu = m_mpq_lar_core_solver.need_to_presolve_with_double_solver(); - lp_assert(!use_lu || A_r().column_count() == A_d().column_count()); - m_mpq_lar_core_solver.m_column_types.push_back(column_type::free_column); - m_columns_with_changed_bound.increase_size_by_one(); - m_rows_with_changed_bounds.increase_size_by_one(); - m_inf_int_set.increase_size_by_one(); - add_new_var_to_core_fields_for_mpq(true); - if (use_lu) - add_new_var_to_core_fields_for_doubles(true); -} - -bool lar_solver::bound_is_integer_if_needed(unsigned j, const mpq & right_side) const { - if (!column_is_int(j)) - return true; - return right_side.is_int(); -} - -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)); - auto vc = new lar_var_constraint(j, kind, right_side); - m_constraints.push_back(vc); - update_column_type_and_bound(j, kind, right_side, ci); - } - else { - add_var_bound_on_constraint_for_term(j, kind, right_side, ci); - } - lp_assert(sizes_are_correct()); - lp_assert(inf_int_set_is_correct()); - return ci; -} - -void lar_solver::update_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_index) { - switch (m_mpq_lar_core_solver.m_column_types[j]) { - case column_type::free_column: - update_free_column_type_and_bound(j, kind, right_side, constr_index); - break; - case column_type::boxed: - update_boxed_column_type_and_bound(j, kind, right_side, constr_index); - break; - case column_type::low_bound: - update_low_bound_column_type_and_bound(j, kind, right_side, constr_index); - break; - case column_type::upper_bound: - update_upper_bound_column_type_and_bound(j, kind, right_side, constr_index); - break; - case column_type::fixed: - update_fixed_column_type_and_bound(j, kind, right_side, constr_index); - break; - default: - lp_assert(false); // cannot be here - } -} - -void lar_solver::add_var_bound_on_constraint_for_term(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) { - lp_assert(is_term(j)); - unsigned adjusted_term_index = adjust_term_index(j); - lp_assert(!term_is_int(m_terms[adjusted_term_index]) || right_side.is_int()); - auto it = m_ext_vars_to_columns.find(j); - if (it != m_ext_vars_to_columns.end()) { - unsigned term_j = it->second.ext_j(); - mpq rs = right_side - m_terms[adjusted_term_index]->m_v; - m_constraints.push_back(new lar_term_constraint(m_terms[adjusted_term_index], kind, right_side)); - update_column_type_and_bound(term_j, kind, rs, ci); - } - else { - add_constraint_from_term_and_create_new_column_row(j, m_terms[adjusted_term_index], kind, right_side); - } -} - -constraint_index lar_solver::add_constraint(const vector>& left_side_with_terms, lconstraint_kind kind_par, const mpq& right_side_parm) { - vector> left_side; - mpq rs = -right_side_parm; - substitute_terms_in_linear_expression(left_side_with_terms, left_side, rs); - unsigned term_index = add_term(left_side, zero_of_type()); - constraint_index ci = m_constraints.size(); - add_var_bound_on_constraint_for_term(term_index, kind_par, -rs, ci); - return ci; -} - -void lar_solver::add_constraint_from_term_and_create_new_column_row(unsigned term_j, const lar_term* term, - lconstraint_kind kind, const mpq & right_side) { - - add_row_from_term_no_constraint(term, term_j); - unsigned j = A_r().column_count() - 1; - update_column_type_and_bound(j, kind, right_side - term->m_v, m_constraints.size()); - m_constraints.push_back(new lar_term_constraint(term, kind, right_side)); - lp_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_costs.size()); -} - -void lar_solver::decide_on_strategy_and_adjust_initial_state() { - lp_assert(strategy_is_undecided()); - 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 { - m_settings.simplex_strategy() = simplex_strategy_enum::tableau_rows; // todo: when to switch to tableau_costs? - } - adjust_initial_state(); -} - -void lar_solver::adjust_initial_state() { - switch (m_settings.simplex_strategy()) { - case simplex_strategy_enum::lu: - adjust_initial_state_for_lu(); - break; - case simplex_strategy_enum::tableau_rows: - adjust_initial_state_for_tableau_rows(); - break; - case simplex_strategy_enum::tableau_costs: - lp_assert(false); // not implemented - case simplex_strategy_enum::undecided: - adjust_initial_state_for_tableau_rows(); - break; - } -} - -void lar_solver::adjust_initial_state_for_lu() { - copy_from_mpq_matrix(A_d()); - unsigned n = A_d().column_count(); - m_mpq_lar_core_solver.m_d_x.resize(n); - m_mpq_lar_core_solver.m_d_low_bounds.resize(n); - m_mpq_lar_core_solver.m_d_upper_bounds.resize(n); - m_mpq_lar_core_solver.m_d_heading = m_mpq_lar_core_solver.m_r_heading; - m_mpq_lar_core_solver.m_d_basis = m_mpq_lar_core_solver.m_r_basis; - - /* - unsigned j = A_d().column_count(); - A_d().add_column(); - lp_assert(m_mpq_lar_core_solver.m_d_x.size() == j); - // lp_assert(m_mpq_lar_core_solver.m_d_low_bounds.size() == j && m_mpq_lar_core_solver.m_d_upper_bounds.size() == j); // restore later - m_mpq_lar_core_solver.m_d_x.resize(j + 1 ); - m_mpq_lar_core_solver.m_d_low_bounds.resize(j + 1); - m_mpq_lar_core_solver.m_d_upper_bounds.resize(j + 1); - lp_assert(m_mpq_lar_core_solver.m_d_heading.size() == j); // as A().column_count() on the entry to the method - if (register_in_basis) { - A_d().add_row(); - m_mpq_lar_core_solver.m_d_heading.push_back(m_mpq_lar_core_solver.m_d_basis.size()); - m_mpq_lar_core_solver.m_d_basis.push_back(j); - }else { - m_mpq_lar_core_solver.m_d_heading.push_back(- static_cast(m_mpq_lar_core_solver.m_d_nbasis.size()) - 1); - m_mpq_lar_core_solver.m_d_nbasis.push_back(j); - }*/ -} - -void lar_solver::adjust_initial_state_for_tableau_rows() { - for (unsigned j = 0; j < m_terms.size(); j++) { - if (contains(m_ext_vars_to_columns, j + m_terms_start_index)) - continue; - add_row_from_term_no_constraint(m_terms[j], j + m_terms_start_index); - } -} - -// this fills the last row of A_d and sets the basis column: -1 in the last column of the row -void lar_solver::fill_last_row_of_A_d(static_matrix & A, const lar_term* ls) { - lp_assert(A.row_count() > 0); - lp_assert(A.column_count() > 0); - unsigned last_row = A.row_count() - 1; - lp_assert(A.m_rows[last_row].empty()); - - for (auto & t : ls->m_coeffs) { - lp_assert(!is_zero(t.second)); - var_index j = t.first; - A.set(last_row, j, -t.second.get_double()); - } - - unsigned basis_j = A.column_count() - 1; - A.set(last_row, basis_j, -1); -} - -void lar_solver::update_free_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_ind) { - mpq y_of_bound(0); - switch (kind) { - case LT: - y_of_bound = -1; - case LE: - m_mpq_lar_core_solver.m_column_types[j] = column_type::upper_bound; - lp_assert(m_mpq_lar_core_solver.m_column_types()[j] == column_type::upper_bound); - lp_assert(m_mpq_lar_core_solver.m_r_upper_bounds.size() > j); - { - auto up = numeric_pair(right_side, y_of_bound); - m_mpq_lar_core_solver.m_r_upper_bounds[j] = up; - } - set_upper_bound_witness(j, constr_ind); - break; - case GT: - y_of_bound = 1; - case GE: - m_mpq_lar_core_solver.m_column_types[j] = column_type::low_bound; - lp_assert(m_mpq_lar_core_solver.m_r_upper_bounds.size() > j); - { - auto low = numeric_pair(right_side, y_of_bound); - m_mpq_lar_core_solver.m_r_low_bounds[j] = low; - } - set_low_bound_witness(j, constr_ind); - break; - case EQ: - 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] = numeric_pair(right_side, zero_of_type()); - set_upper_bound_witness(j, constr_ind); - set_low_bound_witness(j, constr_ind); - break; - - default: - lp_unreachable(); - - } - m_columns_with_changed_bound.insert(j); -} - -void lar_solver::update_upper_bound_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) { - lp_assert(m_mpq_lar_core_solver.m_column_types()[j] == column_type::upper_bound); - mpq y_of_bound(0); - switch (kind) { - case LT: - y_of_bound = -1; - case LE: - { - auto up = numeric_pair(right_side, y_of_bound); - if (up < m_mpq_lar_core_solver.m_r_upper_bounds()[j]) { - m_mpq_lar_core_solver.m_r_upper_bounds[j] = up; - set_upper_bound_witness(j, ci); - m_columns_with_changed_bound.insert(j); - } - } - break; - case GT: - y_of_bound = 1; - case GE: - m_mpq_lar_core_solver.m_column_types[j] = column_type::boxed; - { - auto low = numeric_pair(right_side, y_of_bound); - m_mpq_lar_core_solver.m_r_low_bounds[j] = low; - 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 = lp_status::INFEASIBLE; - m_infeasible_column_index = j; - } - else { - m_mpq_lar_core_solver.m_column_types[j] = m_mpq_lar_core_solver.m_r_low_bounds()[j] < m_mpq_lar_core_solver.m_r_upper_bounds()[j] ? column_type::boxed : column_type::fixed; - } - } - break; - case EQ: - { - auto v = numeric_pair(right_side, zero_of_type()); - if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) { - m_status = lp_status::INFEASIBLE; - set_low_bound_witness(j, ci); - m_infeasible_column_index = j; - } - else { - m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = v; - m_columns_with_changed_bound.insert(j); - set_low_bound_witness(j, ci); - set_upper_bound_witness(j, ci); - m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed; - } - break; - } - break; - - default: - lp_unreachable(); - - } -} - -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 == 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: - y_of_bound = -1; - case LE: - { - auto up = numeric_pair(right_side, y_of_bound); - if (up < m_mpq_lar_core_solver.m_r_upper_bounds[j]) { - m_mpq_lar_core_solver.m_r_upper_bounds[j] = up; - set_upper_bound_witness(j, ci); - m_columns_with_changed_bound.insert(j); - } - - if (up < m_mpq_lar_core_solver.m_r_low_bounds[j]) { - m_status = lp_status::INFEASIBLE; - lp_assert(false); - m_infeasible_column_index = j; - } - else { - if (m_mpq_lar_core_solver.m_r_low_bounds()[j] == m_mpq_lar_core_solver.m_r_upper_bounds()[j]) - m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed; - } - } - break; - case GT: - y_of_bound = 1; - case GE: - { - auto low = numeric_pair(right_side, y_of_bound); - if (low > m_mpq_lar_core_solver.m_r_low_bounds[j]) { - m_mpq_lar_core_solver.m_r_low_bounds[j] = low; - m_columns_with_changed_bound.insert(j); - set_low_bound_witness(j, ci); - } - if (low > m_mpq_lar_core_solver.m_r_upper_bounds[j]) { - m_status = lp_status::INFEASIBLE; - m_infeasible_column_index = j; - } - else if (low == m_mpq_lar_core_solver.m_r_upper_bounds[j]) { - m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed; - } - } - break; - case EQ: - { - auto v = numeric_pair(right_side, zero_of_type()); - if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) { - 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 = lp_status::INFEASIBLE; - m_infeasible_column_index = j; - set_low_bound_witness(j, ci); - } - else { - m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = v; - set_low_bound_witness(j, ci); - set_upper_bound_witness(j, ci); - m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed; - m_columns_with_changed_bound.insert(j); - } - - break; - } - - default: - lp_unreachable(); - - } -} -void lar_solver::update_low_bound_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) { - lp_assert(m_mpq_lar_core_solver.m_column_types()[j] == column_type::low_bound); - mpq y_of_bound(0); - switch (kind) { - case LT: - y_of_bound = -1; - case LE: - { - auto up = numeric_pair(right_side, y_of_bound); - m_mpq_lar_core_solver.m_r_upper_bounds[j] = up; - set_upper_bound_witness(j, ci); - m_columns_with_changed_bound.insert(j); - - if (up < m_mpq_lar_core_solver.m_r_low_bounds[j]) { - m_status = lp_status::INFEASIBLE; - m_infeasible_column_index = j; - } - else { - m_mpq_lar_core_solver.m_column_types[j] = m_mpq_lar_core_solver.m_r_low_bounds()[j] < m_mpq_lar_core_solver.m_r_upper_bounds()[j] ? column_type::boxed : column_type::fixed; - } - } - break; - case GT: - y_of_bound = 1; - case GE: - { - auto low = numeric_pair(right_side, y_of_bound); - if (low > m_mpq_lar_core_solver.m_r_low_bounds[j]) { - m_mpq_lar_core_solver.m_r_low_bounds[j] = low; - m_columns_with_changed_bound.insert(j); - set_low_bound_witness(j, ci); - } - } - break; - case EQ: - { - auto v = numeric_pair(right_side, zero_of_type()); - if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) { - m_status = lp_status::INFEASIBLE; - m_infeasible_column_index = j; - set_upper_bound_witness(j, ci); - } - else { - m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = v; - set_low_bound_witness(j, ci); - set_upper_bound_witness(j, ci); - m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed; - } - m_columns_with_changed_bound.insert(j); - break; - } - - default: - lp_unreachable(); - - } -} - -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 == 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 = lp_status::INFEASIBLE; - m_infeasible_column_index = j; - set_upper_bound_witness(j, ci); - } - break; - case LE: - { - if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) { - m_status = lp_status::INFEASIBLE; - m_infeasible_column_index = j; - set_upper_bound_witness(j, ci); - } - } - break; - case GT: - { - if (v >= m_mpq_lar_core_solver.m_r_upper_bounds[j]) { - m_status = lp_status::INFEASIBLE; - m_infeasible_column_index = j; - set_low_bound_witness(j, ci); - } - } - break; - case GE: - { - if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) { - m_status = lp_status::INFEASIBLE; - m_infeasible_column_index = j; - set_low_bound_witness(j, ci); - } - } - break; - case EQ: - { - if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) { - 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 = lp_status::INFEASIBLE; - m_infeasible_column_index = j; - set_low_bound_witness(j, ci); - } - break; - } - - default: - lp_unreachable(); - - } -} - - -} // namespace lp - -