mirrors/z3
3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2025-04-13 12:28:44 +00:00
Code Activity

solve send-more-money_lev.smt2

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

handle integer vars in random_update

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

call the assert in gomory_cut and branching to a correct place

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

fixes in goromy cut

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

disable x values tracking in random_update

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

more fixes in gomory cut

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

change in mk_bound by Nikolaj

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

fixes in gomory cut and setup

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

fixes in int_solver

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

change a printout

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

fix by Nikolaj in treating terms returned by int_solver

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

fix syntax

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

fix a free coefficient bug in bound propagaion and simplify gomory cut

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

avoid tracking pivoted rows during int_solver::check()
This commit is contained in:
Lev Nachmanson 2017-07-27 10:49:00 -07:00 committed by Lev Nachmanson
parent aba7dcab3e
commit db8f01894f
31 changed files with 894 additions and 767 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

4
src/ast/arith_decl_plugin.h
View file

@ -259,9 +259,9 @@ public:
bool is_uminus(expr const * n) const { return is_app_of(n, m_afid, OP_UMINUS); }
bool is_mul(expr const * n) const { return is_app_of(n, m_afid, OP_MUL); }
bool is_div(expr const * n) const { return is_app_of(n, m_afid, OP_DIV); }
//bool is_div0(expr const * n) const { return is_app_of(n, m_afid, OP_DIV_0); }
bool is_div0(expr const * n) const { return is_app_of(n, m_afid, OP_DIV_0); }
bool is_idiv(expr const * n) const { return is_app_of(n, m_afid, OP_IDIV); }
//bool is_idiv0(expr const * n) const { return is_app_of(n, m_afid, OP_IDIV_0); }
bool is_idiv0(expr const * n) const { return is_app_of(n, m_afid, OP_IDIV_0); }
bool is_mod(expr const * n) const { return is_app_of(n, m_afid, OP_MOD); }
bool is_rem(expr const * n) const { return is_app_of(n, m_afid, OP_REM); }
bool is_to_real(expr const * n) const { return is_app_of(n, m_afid, OP_TO_REAL); }

21
src/smt/smt_setup.cpp
View file

@ -730,8 +730,12 @@ namespace smt {
}
void setup::setup_i_arith() {
// m_context.register_plugin(alloc(smt::theory_lra, m_manager, m_params));
m_context.register_plugin(alloc(smt::theory_i_arith, m_manager, m_params));
if (AS_LRA == m_params.m_arith_mode) {
setup_r_arith();
}
else {
m_context.register_plugin(alloc(smt::theory_i_arith, m_manager, m_params));
}
}
void setup::setup_r_arith() {
@ -739,14 +743,21 @@ namespace smt {
}
void setup::setup_mi_arith() {
if (m_params.m_arith_mode == AS_OPTINF) {
switch (m_params.m_arith_mode) {
case AS_OPTINF:
m_context.register_plugin(alloc(smt::theory_inf_arith, m_manager, m_params));
}
else {
break;
case AS_LRA:
setup_r_arith();
break;
default:
m_context.register_plugin(alloc(smt::theory_mi_arith, m_manager, m_params));
break;
}
}
void setup::setup_arith() {
static_features st(m_manager);
IF_VERBOSE(100, verbose_stream() << "(smt.collecting-features)\n";);

3
src/smt/theory_arith_int.h
View file

@ -1385,6 +1385,7 @@ namespace smt {
m_branch_cut_counter++;
// TODO: add giveup code
TRACE("gomory_cut", tout << m_branch_cut_counter << ", " << m_params.m_arith_branch_cut_ratio << std::endl;);
if (m_branch_cut_counter % m_params.m_arith_branch_cut_ratio == 0) {
TRACE("opt_verbose", display(tout););
move_non_base_vars_to_bounds();
@ -1399,7 +1400,7 @@ namespace smt {
SASSERT(is_base(int_var));
row const & r = m_rows[get_var_row(int_var)];
if (!mk_gomory_cut(r)) {
// silent failure
TRACE("gomory_cut", tout << "silent failure\n";);
}
return FC_CONTINUE;
}

15
src/smt/theory_arith_pp.h
View file

@ -389,8 +389,19 @@ namespace smt {
void theory_arith<Ext>::display_vars(std::ostream & out) const {
out << "vars:\n";
int n = get_num_vars();
for (theory_var v = 0; v < n; v++)
display_var(out, v);
int inf_vars = 0;
int int_inf_vars = 0;
for (theory_var v = 0; v < n; v++) {
if ((lower(v) && lower(v)->get_value() > get_value(v))
|| (upper(v) && upper(v)->get_value() < get_value(v)))
inf_vars++;
if (is_int(v) && !get_value(v).is_int())
int_inf_vars++;
}
out << "infeasibles = " << inf_vars << " int_inf = " << int_inf_vars << std::endl;
for (theory_var v = 0; v < n; v++) {
display_var(out, v);
}
}
template<typename Ext>

220
src/smt/theory_lra.cpp
View file

@ -77,6 +77,7 @@ namespace lp_api {
if (is_true) return inf_rational(m_value); // v >= value or v <= value
if (m_is_int) {
SASSERT(m_value.is_int());
if (m_bound_kind == lower_t) return inf_rational(m_value - rational::one()); // v <= value - 1
return inf_rational(m_value + rational::one()); // v >= value + 1
}
@ -86,7 +87,7 @@ namespace lp_api {
}
}
virtual std::ostream& display(std::ostream& out) const {
return out << "v" << get_var() << " " << get_bound_kind() << " " << m_value;
return out << m_value << " " << get_bound_kind() << " v" << get_var();
}
};
@ -303,7 +304,10 @@ namespace smt {
m_solver->settings().simplex_strategy() = static_cast<lp::simplex_strategy_enum>(lp.simplex_strategy());
reset_variable_values();
m_solver->settings().bound_propagation() = BP_NONE != propagation_mode();
m_solver->set_propagate_bounds_on_pivoted_rows_mode(lp.bprop_on_pivoted_rows());
m_solver->set_track_pivoted_rows(lp.bprop_on_pivoted_rows());
m_solver->settings().m_int_branch_cut_gomory_threshold = ctx().get_fparams().m_arith_branch_cut_ratio;
m_solver->settings().m_run_gcd_test = ctx().get_fparams().m_arith_gcd_test;
m_solver->settings().set_random_seed(ctx().get_fparams().m_random_seed);
//m_solver->settings().set_ostream(0);
m_lia = alloc(lp::int_solver, m_solver.get());
}
@ -318,6 +322,9 @@ namespace smt {
}
void found_not_handled(expr* n) {
if (a.is_div0(n)) {
return;
}
m_not_handled = n;
if (is_app(n) && is_underspecified(to_app(n))) {
TRACE("arith", tout << "Unhandled: " << mk_pp(n, m) << "\n";);
@ -608,8 +615,8 @@ namespace smt {
}
bool all_zeros(vector<rational> const& v) const {
for (unsigned i = 0; i < v.size(); ++i) {
if (!v[i].is_zero()) {
for (rational const& r : v) {
if (!r.is_zero()) {
return false;
}
}
@ -708,7 +715,7 @@ namespace smt {
m_var_index2theory_var.setx(vi, v, UINT_MAX);
}
m_var_trail.push_back(v);
TRACE("arith", tout << "v" << v << " := " << mk_pp(term, m) << " slack: " << vi << " scopes: " << m_scopes.size() << "\n";
TRACE("arith_verbose", tout << "v" << v << " := " << mk_pp(term, m) << " slack: " << vi << " scopes: " << m_scopes.size() << "\n";
m_solver->print_term(m_solver->get_term(vi), tout); tout << "\n";);
}
rational val;
@ -777,7 +784,7 @@ namespace smt {
updt_unassigned_bounds(v, +1);
m_bounds_trail.push_back(v);
m_bool_var2bound.insert(bv, b);
TRACE("arith", tout << "Internalized " << mk_pp(atom, m) << "\n";);
TRACE("arith_verbose", tout << "Internalized " << mk_pp(atom, m) << "\n";);
mk_bound_axioms(*b);
//add_use_lists(b);
return true;
@ -801,7 +808,7 @@ namespace smt {
if (n1->get_th_var(get_id()) != null_theory_var &&
n2->get_th_var(get_id()) != null_theory_var &&
n1 != n2) {
TRACE("arith", tout << mk_pp(atom, m) << "\n";);
TRACE("arith_verbose", tout << mk_pp(atom, m) << "\n";);
m_arith_eq_adapter.mk_axioms(n1, n2);
}
}
@ -1037,38 +1044,58 @@ namespace smt {
return m_solver->var_is_registered(m_theory_var2var_index[v]);
}
mutable vector<std::pair<lp::var_index, rational>> m_todo_terms;
lp::impq get_ivalue(theory_var v) const {
lp_assert(can_get_ivalue(v));
lp::var_index vi = m_theory_var2var_index[v];
if (!m_solver->is_term(vi))
return m_solver->get_value(vi);
const lp::lar_term& term = m_solver->get_term(vi);
lp::impq result(term.m_v);
for (const auto & i: term.m_coeffs) {
result += m_solver->get_value(i.first) * i.second;
return m_solver->get_column_value(vi);
m_todo_terms.push_back(std::make_pair(vi, rational::one()));
lp::impq result(0);
while (!m_todo_terms.empty()) {
vi = m_todo_terms.back().first;
rational coeff = m_todo_terms.back().second;
m_todo_terms.pop_back();
if (m_solver->is_term(vi)) {
const lp::lar_term& term = m_solver->get_term(vi);
result += term.m_v * coeff;
for (const auto & i: term.m_coeffs) {
m_todo_terms.push_back(std::make_pair(i.first, coeff * i.second));
}
}
else {
result += m_solver->get_column_value(vi) * coeff;
}
}
return result;
}
rational get_value(theory_var v) const {
if (!can_get_value(v)) return rational::zero();
lp::var_index vi = m_theory_var2var_index[v];
if (m_variable_values.count(vi) > 0) {
return m_variable_values[vi];
}
if (m_solver->is_term(vi)) {
const lp::lar_term& term = m_solver->get_term(vi);
rational result = term.m_v;
for (auto i = term.m_coeffs.begin(); i != term.m_coeffs.end(); ++i) {
result += m_variable_values[i->first] * i->second;
m_todo_terms.push_back(std::make_pair(vi, rational::one()));
rational result(0);
while (!m_todo_terms.empty()) {
lp::var_index wi = m_todo_terms.back().first;
rational coeff = m_todo_terms.back().second;
m_todo_terms.pop_back();
if (m_solver->is_term(wi)) {
const lp::lar_term& term = m_solver->get_term(wi);
result += term.m_v * coeff;
for (auto const& i : term.m_coeffs) {
m_todo_terms.push_back(std::make_pair(i.first, i.second * coeff));
}
}
else {
result += m_variable_values[wi] * coeff;
}
m_variable_values[vi] = result;
return result;
}
UNREACHABLE();
return m_variable_values[vi];
m_variable_values[vi] = result;
return result;
}
void init_variable_values() {
@ -1230,9 +1257,12 @@ namespace smt {
// create a bound atom representing term <= k
app_ref mk_bound(lp::lar_term const& term, rational const& k) {
SASSERT(k.is_int());
app_ref t = mk_term(term, true);
app_ref atom(a.mk_le(t, a.mk_numeral(k, true)), m);
app_ref t = mk_term(term, k.is_int());
app_ref atom(a.mk_le(t, a.mk_numeral(k, k.is_int())), m);
expr_ref atom1(m);
proof_ref atomp(m);
ctx().get_simplifier()(atom, atom1, atomp);
atom = to_app(atom1);
TRACE("arith", tout << atom << "\n";
m_solver->print_term(term, tout << "bound atom: "); tout << " <= " << k << "\n";
display(tout);
@ -1243,7 +1273,10 @@ namespace smt {
}
lbool check_lia() {
if (m.canceled()) return l_undef;
if (m.canceled()) {
TRACE("arith", tout << "canceled\n";);
return l_undef;
}
lp::lar_term term;
lp::mpq k;
lp::explanation ex; // TBD, this should be streamlined accross different explanations
@ -1288,7 +1321,10 @@ namespace smt {
lbool check_nra() {
m_use_nra_model = false;
if (m.canceled()) return l_undef;
if (m.canceled()) {
TRACE("arith", tout << "canceled\n";);
return l_undef;
}
if (!m_nra) return l_true;
if (!m_nra->need_check()) return l_true;
m_a1 = 0; m_a2 = 0;
@ -1476,6 +1512,7 @@ namespace smt {
void consume(rational const& v, unsigned j) override {
m_imp.set_evidence(j);
m_imp.m_explanation.push_back(std::make_pair(v, j));
}
};
@ -1512,6 +1549,7 @@ namespace smt {
if (lit == null_literal) {
continue;
}
TRACE("arith", tout << lit << " bound: " << *b << " first: " << first << "\n";);
m_solver->settings().st().m_num_of_implied_bounds ++;
if (first) {
@ -1519,6 +1557,7 @@ namespace smt {
m_core.reset();
m_eqs.reset();
m_params.reset();
m_explanation.clear();
local_bound_propagator bp(*this);
m_solver->explain_implied_bound(be, bp);
}
@ -1529,7 +1568,7 @@ namespace smt {
tout << "\n --> ";
ctx().display_literal_verbose(tout, lit);
tout << "\n";
// display_evidence(tout, m_explanation);
display_evidence(tout, m_explanation);
m_solver->print_implied_bound(be, tout);
);
DEBUG_CODE(
@ -1547,8 +1586,8 @@ namespace smt {
SASSERT(validate_assign(lit));
if (m_core.size() < small_lemma_size() && m_eqs.empty()) {
m_core2.reset();
for (unsigned i = 0; i < m_core.size(); ++i) {
m_core2.push_back(~m_core[i]);
for (auto const& c : m_core) {
m_core2.push_back(~c);
}
m_core2.push_back(lit);
justification * js = nullptr;
@ -1569,27 +1608,27 @@ namespace smt {
literal is_bound_implied(lp::lconstraint_kind k, rational const& value, lp_api::bound const& b) const {
if ((k == lp::LE || k == lp::LT) && b.get_bound_kind() == lp_api::upper_t && value <= b.get_value()) {
// v <= value <= b.get_value() => v <= b.get_value()
TRACE("arith", tout << "v <= value <= b.get_value() => v <= b.get_value() \n";);
return literal(b.get_bv(), false);
}
if ((k == lp::GE || k == lp::GT) && b.get_bound_kind() == lp_api::lower_t && b.get_value() <= value) {
// b.get_value() <= value <= v => b.get_value() <= v
TRACE("arith", tout << "b.get_value() <= value <= v => b.get_value() <= v \n";);
return literal(b.get_bv(), false);
}
if (k == lp::LE && b.get_bound_kind() == lp_api::lower_t && value < b.get_value()) {
// v <= value < b.get_value() => v < b.get_value()
TRACE("arith", tout << "v <= value < b.get_value() => v < b.get_value()\n";);
return literal(b.get_bv(), true);
}
if (k == lp::LT && b.get_bound_kind() == lp_api::lower_t && value <= b.get_value()) {
// v < value <= b.get_value() => v < b.get_value()
TRACE("arith", tout << "v < value <= b.get_value() => v < b.get_value()\n";);
return literal(b.get_bv(), true);
}
if (k == lp::GE && b.get_bound_kind() == lp_api::upper_t && b.get_value() < value) {
// b.get_value() < value <= v => b.get_value() < v
TRACE("arith", tout << "b.get_value() < value <= v => b.get_value() < v\n";);
return literal(b.get_bv(), true);
}
if (k == lp::GT && b.get_bound_kind() == lp_api::upper_t && b.get_value() <= value) {
// b.get_value() <= value < v => b.get_value() < v
TRACE("arith", tout << "b.get_value() <= value < v => b.get_value() < v\n";);
return literal(b.get_bv(), true);
}
@ -1923,16 +1962,29 @@ namespace smt {
++m_stats.m_bounds_propagations;
}
svector<lp::var_index> m_todo_vars;
void add_use_lists(lp_api::bound* b) {
theory_var v = b->get_var();
lp::var_index vi = get_var_index(v);
if (m_solver->is_term(vi)) {
if (!m_solver->is_term(vi)) {
return;
}
m_todo_vars.push_back(vi);
while (!m_todo_vars.empty()) {
vi = m_todo_vars.back();
m_todo_vars.pop_back();
lp::lar_term const& term = m_solver->get_term(vi);
for (auto i = term.m_coeffs.begin(); i != term.m_coeffs.end(); ++i) {
lp::var_index wi = i->first;
unsigned w = m_var_index2theory_var[wi];
m_use_list.reserve(w + 1, ptr_vector<lp_api::bound>());
m_use_list[w].push_back(b);
if (m_solver->is_term(wi)) {
m_todo_vars.push_back(wi);
}
else {
unsigned w = m_var_index2theory_var[wi];
m_use_list.reserve(w + 1, ptr_vector<lp_api::bound>());
m_use_list[w].push_back(b);
}
}
}
}
@ -1940,13 +1992,24 @@ namespace smt {
void del_use_lists(lp_api::bound* b) {
theory_var v = b->get_var();
lp::var_index vi = m_theory_var2var_index[v];
if (m_solver->is_term(vi)) {
if (!m_solver->is_term(vi)) {
return;
}
m_todo_vars.push_back(vi);
while (!m_todo_vars.empty()) {
vi = m_todo_vars.back();
m_todo_vars.pop_back();
lp::lar_term const& term = m_solver->get_term(vi);
for (auto i = term.m_coeffs.begin(); i != term.m_coeffs.end(); ++i) {
lp::var_index wi = i->first;
unsigned w = m_var_index2theory_var[wi];
SASSERT(m_use_list[w].back() == b);
m_use_list[w].pop_back();
if (m_solver->is_term(wi)) {
m_todo_vars.push_back(wi);
}
else {
unsigned w = m_var_index2theory_var[wi];
SASSERT(m_use_list[w].back() == b);
m_use_list[w].pop_back();
}
}
}
}
@ -2033,6 +2096,9 @@ namespace smt {
lp::constraint_index ci;
rational value;
bool is_strict;
if (m_solver->is_term(wi)) {
return false;
}
if (coeff.second.is_neg() == is_lub) {
// -3*x ... <= lub based on lower bound for x.
if (!m_solver->has_lower_bound(wi, ci, value, is_strict)) {
@ -2378,15 +2444,31 @@ namespace smt {
SASSERT(m_use_nra_model);
lp::var_index vi = m_theory_var2var_index[v];
if (m_solver->is_term(vi)) {
lp::lar_term const& term = m_solver->get_term(vi);
scoped_anum r1(m_nra->am());
m_nra->am().set(r, term.m_v.to_mpq());
for (auto const coeff : term.m_coeffs) {
lp::var_index wi = coeff.first;
m_nra->am().set(r1, coeff.second.to_mpq());
m_nra->am().mul(m_nra->value(wi), r1, r1);
m_nra->am().add(r1, r, r);
m_todo_terms.push_back(std::make_pair(vi, rational::one()));
m_nra->am().set(r, 0);
while (!m_todo_terms.empty()) {
rational wcoeff = m_todo_terms.back().second;
lp::var_index wi = m_todo_terms.back().first; // todo : got a warning "wi is not used"
m_todo_terms.pop_back();
lp::lar_term const& term = m_solver->get_term(vi);
scoped_anum r1(m_nra->am());
rational c1 = term.m_v * wcoeff;
m_nra->am().set(r1, c1.to_mpq());
m_nra->am().add(r, r1, r);
for (auto const coeff : term.m_coeffs) {
lp::var_index wi = coeff.first;
c1 = coeff.second * wcoeff;
if (m_solver->is_term(wi)) {
m_todo_terms.push_back(std::make_pair(wi, c1));
}
else {
m_nra->am().set(r1, c1.to_mpq());
m_nra->am().mul(m_nra->value(wi), r1, r1);
m_nra->am().add(r1, r, r);
}
}
}
return r;
}
@ -2431,11 +2513,22 @@ namespace smt {
// Auxiliary verification utilities.
struct scoped_arith_mode {
smt_params& p;
scoped_arith_mode(smt_params& p) : p(p) {
p.m_arith_mode = AS_ARITH;
}
~scoped_arith_mode() {
p.m_arith_mode = AS_LRA;
}
};
bool validate_conflict() {
if (dump_lemmas()) {
ctx().display_lemma_as_smt_problem(m_core.size(), m_core.c_ptr(), m_eqs.size(), m_eqs.c_ptr(), false_literal);
}
if (m_arith_params.m_arith_mode == AS_LRA) return true;
if (m_arith_params.m_arith_mode != AS_LRA) return true;
scoped_arith_mode _sa(ctx().get_fparams());
context nctx(m, ctx().get_fparams(), ctx().get_params());
add_background(nctx);
bool result = l_true != nctx.check();
@ -2448,7 +2541,8 @@ namespace smt {
if (dump_lemmas()) {
ctx().display_lemma_as_smt_problem(m_core.size(), m_core.c_ptr(), m_eqs.size(), m_eqs.c_ptr(), lit);
}
if (m_arith_params.m_arith_mode == AS_LRA) return true;
if (m_arith_params.m_arith_mode != AS_LRA) return true;
scoped_arith_mode _sa(ctx().get_fparams());
context nctx(m, ctx().get_fparams(), ctx().get_params());
m_core.push_back(~lit);
add_background(nctx);
@ -2486,10 +2580,14 @@ namespace smt {
theory_lra::inf_eps maximize(theory_var v, expr_ref& blocker, bool& has_shared) {
lp::var_index vi = m_theory_var2var_index.get(v, UINT_MAX);
vector<std::pair<rational, lp::var_index> > coeffs;
rational coeff;
rational coeff(0);
//
// TBD: change API for maximize_term to take a proper term as input.
//
if (m_solver->is_term(vi)) {
const lp::lar_term& term = m_solver->get_term(vi);
for (auto & ti : term.m_coeffs) {
SASSERT(!m_solver->is_term(ti.first));
coeffs.push_back(std::make_pair(ti.second, ti.first));
}
coeff = term.m_v;
@ -2547,7 +2645,13 @@ namespace smt {
app_ref mk_term(lp::lar_term const& term, bool is_int) {
expr_ref_vector args(m);
for (auto & ti : term.m_coeffs) {
theory_var w = m_var_index2theory_var[ti.first];
theory_var w;
if (m_solver->is_term(ti.first)) {
w = m_term_index2theory_var[m_solver->adjust_term_index(ti.first)];
}
else {
w = m_var_index2theory_var[ti.first];
}
expr* o = get_enode(w)->get_owner();
if (ti.second.is_one()) {
args.push_back(o);

4
src/smt/theory_str.cpp
View file

@ -11074,8 +11074,8 @@ namespace smt {
ctx.force_phase(lit);
}
zstring aStr = gen_val_string(len, options[i - l]);
expr * strAst;
zstring aStr = gen_val_string(len, options[static_cast<int>(i) - static_cast<int>(l)]);
expr * strAst;
if (m_params.m_UseFastValueTesterCache) {
if (!valueTesterCache.find(aStr, strAst)) {
strAst = mk_string(aStr);

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

@ -1124,7 +1124,7 @@ void update_settings(argument_parser & args_parser, lp_settings& settings) {
settings.harris_feasibility_tolerance = d;
}
if (get_int_from_args_parser("--random_seed", args_parser, n)) {
settings.random_seed(n);
settings.set_random_seed(n);
}
if (get_int_from_args_parser("--simplex_strategy", args_parser, n)) {
settings.simplex_strategy() = static_cast<simplex_strategy_enum>(n);

19
src/util/lp/bound_propagator.cpp
View file

@ -15,8 +15,13 @@ const impq & bound_propagator::get_low_bound(unsigned j) const {
const impq & bound_propagator::get_upper_bound(unsigned j) const {
return m_lar_solver.m_mpq_lar_core_solver.m_r_upper_bounds()[j];
}
void bound_propagator::try_add_bound(const mpq & v, unsigned j, bool is_low, bool coeff_before_j_is_pos, unsigned row_or_term_index, bool strict) {
j = m_lar_solver.adjust_column_index_to_term_index(j);
void bound_propagator::try_add_bound(mpq v, unsigned j, bool is_low, bool coeff_before_j_is_pos, unsigned row_or_term_index, bool strict) {
j = m_lar_solver.adjust_column_index_to_term_index(j);
if (m_lar_solver.is_term(j)) {
// lp treats terms as not having a free coefficient, restoring it below for the outside consumption
v += m_lar_solver.get_term(j).m_v;
}
lconstraint_kind kind = is_low? GE : LE;
if (strict)
kind = static_cast<lconstraint_kind>(kind / 2);
@ -27,20 +32,26 @@ void bound_propagator::try_add_bound(const mpq & v, unsigned j, bool is_low, boo
if (is_low) {
if (try_get_val(m_improved_low_bounds, j, k)) {
auto & found_bound = m_ibounds[k];
if (v > found_bound.m_bound || (v == found_bound.m_bound && found_bound.m_strict == false && strict))
if (v > found_bound.m_bound || (v == found_bound.m_bound && found_bound.m_strict == false && strict)) {
found_bound = implied_bound(v, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict);
TRACE("try_add_bound", m_lar_solver.print_implied_bound(found_bound, tout););
}
} else {
m_improved_low_bounds[j] = m_ibounds.size();
m_ibounds.push_back(implied_bound(v, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict));
TRACE("try_add_bound", m_lar_solver.print_implied_bound(m_ibounds.back(), tout););
}
} else { // the upper bound case
if (try_get_val(m_improved_upper_bounds, j, k)) {
auto & found_bound = m_ibounds[k];
if (v < found_bound.m_bound || (v == found_bound.m_bound && found_bound.m_strict == false && strict))
if (v < found_bound.m_bound || (v == found_bound.m_bound && found_bound.m_strict == false && strict)) {
found_bound = implied_bound(v, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict);
TRACE("try_add_bound", m_lar_solver.print_implied_bound(found_bound, tout););
}
} else {
m_improved_upper_bounds[j] = m_ibounds.size();
m_ibounds.push_back(implied_bound(v, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict));
TRACE("try_add_bound", m_lar_solver.print_implied_bound(m_ibounds.back(), tout););
}
}
}

4
src/util/lp/column_namer.h
View file

@ -40,7 +40,7 @@ public:
T a;
unsigned i;
while (it->next(a, i)) {
coeff.emplace_back(a, i);
coeff.push_back(std::make_pair(a, i));
}
print_linear_combination_of_column_indices_only(coeff, out);
}
@ -65,7 +65,7 @@ public:
else if (val != numeric_traits<T>::one())
out << T_to_string(val);
out << "_" << it.second;
out << "v" << it.second;
}
}

4
src/util/lp/implied_bound.h
View file

@ -28,7 +28,6 @@ struct implied_bound {
bool m_coeff_before_j_is_pos;
unsigned m_row_or_term_index;
bool m_strict;
lconstraint_kind kind() const {
lconstraint_kind k = m_is_low_bound? GE : LE;
if (m_strict)
@ -52,6 +51,7 @@ struct implied_bound {
m_is_low_bound(low_bound),
m_coeff_before_j_is_pos(coeff_before_j_is_pos),
m_row_or_term_index(row_or_term_index),
m_strict(strict) {}
m_strict(strict) {
}
};
}

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

@ -8,7 +8,6 @@
namespace lp {
void int_solver::fix_non_base_columns() {
lp_assert(is_feasible() && inf_int_set_is_correct());
auto & lcs = m_lar_solver->m_mpq_lar_core_solver;
bool change = false;
for (unsigned j : lcs.m_r_nbasis) {
@ -66,6 +65,10 @@ 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;
@ -81,11 +84,12 @@ int int_solver::find_inf_int_boxed_base_column_with_smallest_range() {
mpq range;
mpq new_range;
mpq small_range_thresold(1024);
unsigned n = 0;
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;
@ -104,8 +108,8 @@ int int_solver::find_inf_int_boxed_base_column_with_smallest_range() {
continue;
}
if (new_range == range) {
n++;
if (settings().random_next() % n == 0) {
lp_assert(n >= 1);
if (settings().random_next() % (++n) == 0) {
result = j;
continue;
}
@ -115,32 +119,30 @@ int int_solver::find_inf_int_boxed_base_column_with_smallest_range() {
}
bool int_solver::is_gomory_cut_target() {
m_iter_on_gomory_row->reset();
bool int_solver::is_gomory_cut_target(linear_combination_iterator<mpq> &iter) {
unsigned j;
TRACE("gomory_cut", m_lar_solver->print_linear_iterator(m_iter_on_gomory_row, tout);
m_iter_on_gomory_row->reset();
);
while (m_iter_on_gomory_row->next(j)) {
// All non base variables must be at their bounds and assigned to rationals (that is, infinitesimals are not allowed).
if (j != m_gomory_cut_inf_column && (!at_bound(j) || !is_zero(get_value(j).y))) {
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 << "at_bound: " << at_bound(j) << "\n";
tout << "infinitesimal: " << !is_zero(get_value(j).y) << "\n";);
iter.reset();
return false;
}
}
m_iter_on_gomory_row->reset();
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) {
mpq f_0 = fractional_part(get_value(m_gomory_cut_inf_column));
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_lower(x_j)) {
if (at_low(x_j)) {
if (a.is_pos()) {
new_a = a / (1 - f_0);
}
@ -148,8 +150,8 @@ void int_solver::real_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, l
new_a = a / f_0;
new_a.neg();
}
k += lower_bound(x_j).x * k; // k.addmul(new_a, lower_bound(x_j).x); // is it a faster operation
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 {
@ -161,11 +163,11 @@ void int_solver::real_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, l
else {
new_a = a / (mpq(1) - f_0);
}
k += upper_bound(x_j).x * k; // k.addmul(new_a, upper_bound(x_j).get_rational());
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", tout << a << "*v" << x_j << " k: " << k << "\n";);
pol.add_monoid(new_a, x_j);
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 {
@ -177,152 +179,203 @@ constraint_index int_solver::column_low_bound_constraint(unsigned j) const {
}
void int_solver::int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term & pol, explanation& expl, mpq & lcm_den) {
mpq f_0 = fractional_part(get_value(m_gomory_cut_inf_column));
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<mpq>() && f_0 < one_of_type<mpq>());
mpq f_j = fractional_part(a);
TRACE("gomory_cut_detail",
tout << a << "*v" << x_j << "\n";
tout << "fractional_part: " << fractional_part(a) << "\n";
tout << a << " x_j" << x_j << " k = " << k << "\n";
tout << "f_j: " << f_j << "\n";
tout << "f_0: " << f_0 << "\n";
tout << "one_minus_f_0: " << 1 - f_0 << "\n";);
if (!f_j.is_zero()) {
mpq new_a;
if (at_lower(x_j)) {
auto one_minus_f_0 = 1 - f_0;
if (f_j <= one_minus_f_0) {
new_a = f_j / one_minus_f_0;
}
else {
new_a = (1 - f_j) / f_0;
}
k.addmul(new_a, lower_bound(x_j).x);
expl.push_justification(column_low_bound_constraint(x_j), new_a);
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 {
SASSERT(at_upper(x_j));
if (f_j <= f_0) {
new_a = f_j / f_0;
}
else {
new_a = (mpq(1) - f_j) / 1 - f_0;
}
new_a.neg(); // the upper terms are inverted
k.addmul(new_a, upper_bound(x_j).x);
expl.push_justification(column_upper_bound_constraint(x_j), new_a);
new_a = (1 - f_j) / f_0;
}
TRACE("gomory_cut_detail", tout << "new_a: " << new_a << " k: " << k << "\n";);
pol.add_monoid(new_a, x_j);
lcm_den = lcm(lcm_den, denominator(new_a));
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",
display_row_info(tout,
m_lar_solver->m_mpq_lar_core_solver.m_r_heading[m_gomory_cut_inf_column]););
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_for_size_gt_1(
vector<std::pair<mpq, unsigned>> & pol,
lar_term & t,
mpq &k,
unsigned num_ints,
mpq & lcm_den) {
if (num_ints > 0) {
lcm_den = lcm(lcm_den, denominator(k));
TRACE("gomory_cut_detail", tout << "k: " << k << " lcm_den: " << lcm_den << "\n";
linear_combination_iterator_on_vector<mpq> pi(pol);
m_lar_solver->print_linear_iterator(&pi, tout);
tout << "\nk: " << k << "\n";);
lp_assert(lcm_den.is_pos());
if (!lcm_den.is_one()) {
// normalize coefficients of integer parameters to be integers.
for (auto & pi: pol) {
pi.first *= lcm_den;
SASSERT(!is_int(pi.second) || pi.first.is_int());
}
k *= lcm_den;
}
TRACE("gomory_cut_detail", tout << "after *lcm_den\n";
for (unsigned i = 0; i < pol.size(); i++) {
tout << pol[i].first << " * v" << pol[i].second << "\n";
}
tout << "k: " << k << "\n";);
void int_solver::gomory_cut_adjust_t_and_k(vector<std::pair<mpq, unsigned>> & 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;
}
t.clear();
// negate everything to return -pol <= -k
for (const auto & pi: pol)
t.add_monoid(-pi.first, pi.second);
t.add_monomial(-pi.first, pi.second);
k.neg();
}
void int_solver::gomory_cut_adjust_t_and_k_for_size_1(const vector<std::pair<mpq, unsigned>> & pol, lar_term& t, mpq &k) {
lp_assert(pol.size() == 1);
unsigned j = pol[0].second;
k /= pol[0].first;
bool is_lower = pol[0].first.is_pos();
if (is_int(j) && !k.is_int()) {
k = is_lower?ceil(k):floor(k);
}
if (is_lower) { // returning -t <= -k which is equivalent to t >= k
k.neg();
t.negate();
}
}
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;
}
lia_move int_solver::report_gomory_cut(lar_term& t, mpq& k, mpq &lcm_den, unsigned num_ints) {
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();
if (pol.size() == 1)
gomory_cut_adjust_t_and_k_for_size_1(pol, t, k);
else
gomory_cut_adjust_t_and_k_for_size_gt_1(pol, t, k, num_ints, lcm_den);
m_lar_solver->subs_terms_for_debugging(t); // todo: remove later
return lia_move::cut;
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) {
lia_move int_solver::mk_gomory_cut(lar_term& t, mpq& k, explanation & expl, unsigned inf_col, linear_combination_iterator<mpq>& iter) {
lp_assert(column_is_int_inf(m_gomory_cut_inf_column));
lp_assert(column_is_int_inf(inf_col));
TRACE("gomory_cut", tout << "applying cut at:\n"; m_lar_solver->print_linear_iterator(m_iter_on_gomory_row, tout); tout << std::endl; m_iter_on_gomory_row->reset(););
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 num_ints = 0;
unsigned x_j;
mpq a;
while (m_iter_on_gomory_row->next(a, x_j)) {
if (x_j == m_gomory_cut_inf_column)
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);
real_case_in_gomory_cut(a, x_j, k, t, expl, inf_col);
else {
num_ints++;
int_case_in_gomory_cut(a, x_j, k, t, expl, lcm_den);
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);
return report_gomory_cut(t, k, lcm_den, num_ints);
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() {
@ -331,61 +384,53 @@ void int_solver::init_check_data() {
m_old_values_data.resize(n);
}
int int_solver::find_next_free_var_in_gomory_row() {
lp_assert(m_iter_on_gomory_row != nullptr);
int int_solver::find_free_var_in_gomory_row(linear_combination_iterator<mpq>& iter) {
unsigned j;
while(m_iter_on_gomory_row->next(j)) {
if (j != m_gomory_cut_inf_column && is_free(j))
while(iter.next(j)) {
if (!is_base(j) && is_free(j))
return static_cast<int>(j);
}
iter.reset();
return -1;
}
lia_move int_solver::proceed_with_gomory_cut(lar_term& t, mpq& k, explanation& ex) {
int j = find_next_free_var_in_gomory_row();
if (j != -1) {
m_found_free_var_in_gomory_row = true;
lia_move int_solver::proceed_with_gomory_cut(lar_term& t, mpq& k, explanation& ex, unsigned j, linear_combination_iterator<mpq>& iter) {
int free_j = find_free_var_in_gomory_row(iter);
if (free_j != -1) {
lp_assert(t.is_empty());
t.add_monoid(mpq(1), j);
t.add_monomial(mpq(1), m_lar_solver->adjust_column_index_to_term_index(free_j));
k = zero_of_type<mpq>();
return lia_move::branch; // branch on a free column
}
if (m_found_free_var_in_gomory_row || !is_gomory_cut_target()) {
m_found_free_var_in_gomory_row = false;
delete m_iter_on_gomory_row;
m_iter_on_gomory_row = nullptr;
return lia_move::continue_with_check;
}
lia_move ret = mk_gomory_cut(t, k, ex);
delete m_iter_on_gomory_row;
m_iter_on_gomory_row = nullptr;
return ret;
if (!is_gomory_cut_target(iter))
return create_branch_on_column(j, t, k);
return mk_gomory_cut(t, k, ex, j, iter);
}
lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex) {
if (m_iter_on_gomory_row != nullptr) {
auto ret = proceed_with_gomory_cut(t, k, ex);
TRACE("gomory_cut", tout << "term t = "; m_lar_solver->print_term_as_indices(t, tout););
if (ret != lia_move::continue_with_check)
return ret;
}
unsigned int_solver::row_of_basic_column(unsigned j) const {
return m_lar_solver->m_mpq_lar_core_solver.m_r_heading[j];
}
lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex) {
init_check_data();
lp_assert(inf_int_set_is_correct());
// currently it is a reimplementation of
// it is mostly a reimplementation of
// final_check_status theory_arith<Ext>::check_int_feasibility()
// from theory_arith_int.h
if (m_lar_solver->model_is_int_feasible())
if (!has_inf_int())
return lia_move::ok;
if (!gcd_test(ex))
return lia_move::conflict;
if (settings().m_run_gcd_test)
if (!gcd_test(ex))
return lia_move::conflict;
/*
if (m_params.m_arith_euclidean_solver)
apply_euclidean_solver();
*/
bool track_pivoted_rows = m_lar_solver->get_track_pivoted_rows();
m_lar_solver->set_track_pivoted_rows(false);
m_lar_solver->pivot_fixed_vars_from_basis();
lean_assert(m_lar_solver->m_mpq_lar_core_solver.r_basis_is_OK());
patch_int_infeasible_columns();
@ -395,93 +440,90 @@ lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex) {
lean_assert(is_feasible());
TRACE("arith_int_rows", trace_inf_rows(););
if (find_inf_int_base_column() == -1)
if (!has_inf_int()) {
m_lar_solver->set_track_pivoted_rows(track_pivoted_rows);
return lia_move::ok;
if ((++m_branch_cut_counter) % settings().m_int_branch_cut_threshold == 0) {
move_non_base_vars_to_bounds(); // todo track changed variables
lp_status st = m_lar_solver->find_feasible_solution();
if (st != lp_status::FEASIBLE && st != lp_status::OPTIMAL) {
return lia_move::give_up;
}
lp_assert(inf_int_set_is_correct());
// init_inf_int_set(); // todo - can we avoid this call?
int j = find_inf_int_base_column();
if (j != -1) {
TRACE("arith_int", tout << "j = " << j << " does not have an integer assignment: " << get_value(j) << "\n";);
unsigned row_index = m_lar_solver->m_mpq_lar_core_solver.m_r_heading[j];
if (!mk_gomory_cut(row_index, ex)) {
}
TRACE("gomory_cut", tout << m_branch_cut_counter+1 << ", " << settings().m_int_branch_cut_gomory_threshold << std::endl;);
if ((++m_branch_cut_counter) % settings().m_int_branch_cut_gomory_threshold == 0) {
if (move_non_base_vars_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";);
m_lar_solver->set_track_pivoted_rows(track_pivoted_rows);
return lia_move::give_up;
// silent failure
}
return lia_move::cut;
}
}
else {
int j = find_inf_int_base_column();
if (j != -1) {
TRACE("arith_int", tout << "j" << j << " does not have an integer assignment: " << get_value(j) << "\n";);
lp_assert(t.is_empty());
t.add_monoid(1, j);
k = floor(get_value(j));
TRACE("arith_int", tout << "branching v" << j << " = " << get_value(j) << "\n";
display_column(tout, j);
tout << "k = " << k << std::endl;
);
return lia_move::branch;
}
lp_assert(j != -1);
TRACE("arith_int", tout << "j = " << j << " does not have an integer assignment: " << get_value(j) << "\n";);
auto iter_on_gomory_row = m_lar_solver->get_iterator_on_row(row_of_basic_column(j));
lia_move ret = proceed_with_gomory_cut(t, k, ex, j, *iter_on_gomory_row);
delete iter_on_gomory_row;
m_lar_solver->set_track_pivoted_rows(track_pivoted_rows);
return ret;
}
lp_assert(m_lar_solver->m_mpq_lar_core_solver.r_basis_is_OK());
// return true;
return lia_move::give_up;
m_lar_solver->set_track_pivoted_rows(track_pivoted_rows);
return create_branch_on_column(find_inf_int_base_column(), t, k);
}
void int_solver::move_non_base_vars_to_bounds() {
bool int_solver::move_non_base_vars_to_bounds() {
auto & lcs = m_lar_solver->m_mpq_lar_core_solver;
bool change = false;
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])
set_value(j, lcs.m_r_low_bounds()[j]);
if (val != lcs.m_r_low_bounds()[j] && val != lcs.m_r_upper_bounds()[j]) {
set_value_for_nbasic_column(j, lcs.m_r_low_bounds()[j]);
change = true;
}
break;
case column_type::low_bound:
if (val != lcs.m_r_low_bounds()[j])
set_value(j, lcs.m_r_low_bounds()[j]);
if (val != lcs.m_r_low_bounds()[j]) {
set_value_for_nbasic_column(j, lcs.m_r_low_bounds()[j]);
change = true;
}
break;
case column_type::upper_bound:
if (val != lcs.m_r_upper_bounds()[j])
set_value(j, lcs.m_r_upper_bounds()[j]);
if (val != lcs.m_r_upper_bounds()[j]) {
set_value_for_nbasic_column(j, lcs.m_r_upper_bounds()[j]);
change = true;
}
break;
default:
if (is_int(j) && !val.is_int()) {
set_value(j, impq(floor(val)));
if (is_int(j) && !val.is_int()) {
set_value_for_nbasic_column(j, impq(floor(val)));
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_old_values_set.contains(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;
m_lar_solver->change_basic_x_by_delta_on_column(j, delta);
auto * it = get_column_iterator(j);
update_column_in_int_inf_set(j);
unsigned i;
while (it->next(i))
update_column_in_int_inf_set(m_lar_solver->m_mpq_lar_core_solver.m_r_basis[i]);
delete it;
m_lar_solver->change_basic_columns_dependend_on_a_given_nb_column(j, delta);
}
void int_solver::patch_int_infeasible_columns() {
@ -660,7 +702,7 @@ bool int_solver::ext_gcd_test(iterator_on_row<mpq> & it,
else {
// l += ncoeff * upper_bound(j).get_rational();
l.addmul(ncoeff, m_lar_solver->column_upper_bound(j).x);
// u += ncoeff * lower_bound(j).get_rational();
// 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);
@ -702,10 +744,11 @@ int_solver::int_solver(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<impq>());
m_old_values_data.resize(lar_slv->A_r().column_count(), zero_of_type<impq>());
m_lar_solver->set_int_solver(this);
}
bool int_solver::lower(unsigned j) const {
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:
@ -716,7 +759,7 @@ bool int_solver::lower(unsigned j) const {
}
}
bool int_solver::upper(unsigned j) const {
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:
@ -727,14 +770,6 @@ bool int_solver::upper(unsigned j) const {
}
}
const impq& int_solver::lower_bound(unsigned j) const {
return m_lar_solver->m_mpq_lar_core_solver.m_r_low_bounds()[j];
}
const impq& int_solver::upper_bound(unsigned j) const {
return m_lar_solver->m_mpq_lar_core_solver.m_r_upper_bounds()[j];
}
void set_lower(impq & l,
bool & inf_l,
@ -754,73 +789,62 @@ void set_upper(impq & u,
}
}
bool int_solver::get_freedom_interval_for_column(unsigned x_j, bool & inf_l, impq & l, bool & inf_u, impq & u, mpq & m) {
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[x_j] >= 0) // the basic var
if (lcs.m_r_heading[j] >= 0) // the basic var
return false;
impq const & x_j_val = lcs.m_r_x[x_j];
linear_combination_iterator<mpq> *it = get_column_iterator(x_j);
impq const & xj = get_value(j);
linear_combination_iterator<mpq> *it = get_column_iterator(j);
inf_l = true;
inf_u = true;
l = u = zero_of_type<impq>();
m = mpq(1);
if (lower(x_j)) {
set_lower(l, inf_l, lower_bound(x_j));
if (has_low(j)) {
set_lower(l, inf_l, low_bound(j));
}
if (upper(x_j)) {
set_upper(u, inf_u, upper_bound(x_j));
if (has_upper(j)) {
set_upper(u, inf_u, upper_bound(j));
}
mpq a_ij; unsigned i;
while (it->next(a_ij, i)) {
unsigned x_i = lcs.m_r_basis[i];
impq const & x_i_val = lcs.m_r_x[x_i];
if (is_int(x_i) && is_int(x_j) && !a_ij.is_int())
m = lcm(m, denominator(a_ij));
bool x_i_lower = lower(x_i);
bool x_i_upper = upper(x_i);
if (a_ij.is_neg()) {
if (x_i_lower) {
impq new_l = x_j_val + ((x_i_val - lcs.m_r_low_bounds()[x_i]) / a_ij);
set_lower(l, inf_l, new_l);
if (!inf_l && !inf_u && l == u) break;;
}
if (x_i_upper) {
impq new_u = x_j_val + ((x_i_val - lcs.m_r_upper_bounds()[x_i]) / a_ij);
set_upper(u, inf_u, new_u);
if (!inf_l && !inf_u && l == u) break;;
}
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 (x_i_upper) {
impq new_l = x_j_val + ((x_i_val - lcs.m_r_upper_bounds()[x_i]) / a_ij);
set_lower(l, inf_l, new_l);
if (!inf_l && !inf_u && l == u) break;;
}
if (x_i_lower) {
impq new_u = x_j_val + ((x_i_val - lcs.m_r_low_bounds()[x_i]) / a_ij);
set_upper(u, inf_u, new_u);
if (!inf_l && !inf_u && l == u) break;;
}
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(x_j);
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(x_j) << "\n";
tout << "val = " << get_value(j) << "\n";
);
lp_assert(inf_l || l <= get_value(x_j));
lp_assert(inf_u || u >= get_value(x_j));
lp_assert(inf_l || l <= get_value(j));
lp_assert(inf_u || u >= get_value(j));
return true;
}
@ -834,7 +858,7 @@ bool int_solver::is_real(unsigned j) const {
}
bool int_solver::value_is_int(unsigned j) const {
return m_lar_solver->m_mpq_lar_core_solver.m_r_x[j].is_int();
return m_lar_solver->column_value_is_int(j);
}
@ -856,8 +880,10 @@ void int_solver::display_column(std::ostream & out, unsigned j) const {
bool int_solver::inf_int_set_is_correct() const {
for (unsigned j = 0; j < m_lar_solver->A_r().column_count(); j++) {
if (inf_int_set().contains(j) != (is_int(j) && (!value_is_int(j))))
if (inf_int_set().contains(j) != (is_int(j) && (!value_is_int(j)))) {
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;
}
@ -881,6 +907,10 @@ 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;
}
@ -902,7 +932,7 @@ bool int_solver::at_bound(unsigned j) const {
}
}
bool int_solver::at_lower(unsigned j) const {
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:
@ -950,4 +980,115 @@ void int_solver::display_row_info(std::ostream & out, unsigned row_index) const
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<unsigned>(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) const {
lp_assert(t.is_empty());
lp_assert(j != -1);
t.add_monomial(mpq(1), j);
k = floor(get_value(j));
TRACE("arith_int", tout << "branching v" << j << " = " << get_value(j) << "\n";
display_column(tout, j);
tout << "k = " << k << std::endl;
);
lp_assert(current_solution_is_inf_on_cut(t, k));
m_lar_solver->subs_term_columns(t);
return lia_move::branch;
}
const impq& int_solver::upper_bound(unsigned j) const {
return m_lar_solver->column_upper_bound(j);
}
}

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

@ -34,9 +34,6 @@ public:
int_set m_old_values_set;
vector<impq> m_old_values_data;
unsigned m_branch_cut_counter;
linear_combination_iterator<mpq>* m_iter_on_gomory_row;
unsigned m_gomory_cut_inf_column;
bool m_found_free_var_in_gomory_row;
// methods
int_solver(lar_solver* lp);
@ -80,16 +77,17 @@ private:
void patch_int_infeasible_columns();
bool get_freedom_interval_for_column(unsigned j, bool & inf_l, impq & l, bool & inf_u, impq & u, mpq & m);
linear_combination_iterator<mpq> * get_column_iterator(unsigned j);
bool lower(unsigned j) const;
bool upper(unsigned j) const;
const impq & lower_bound(unsigned j) const;
const impq & low_bound(unsigned j) const;
const impq & upper_bound(unsigned j) const;
bool is_int(unsigned j) const;
bool is_real(unsigned j) const;
bool is_base(unsigned j) const;
bool is_boxed(unsigned j) const;
bool is_fixed(unsigned j) const;
bool is_free(unsigned j) const;
bool value_is_int(unsigned j) const;
void set_value(unsigned j, const impq & new_val);
void set_value_for_nbasic_column(unsigned j, const impq & new_val);
void set_value_for_nbasic_column_ignore_old_values(unsigned j, const impq & new_val);
void fix_non_base_columns();
void failed();
bool is_feasible() const;
@ -102,35 +100,44 @@ private:
int find_inf_int_base_column();
int find_inf_int_boxed_base_column_with_smallest_range();
lp_settings& settings();
void move_non_base_vars_to_bounds();
bool move_non_base_vars_to_bounds();
void branch_infeasible_int_var(unsigned);
lia_move mk_gomory_cut(lar_term& t, mpq& k,explanation & ex);
lia_move mk_gomory_cut(lar_term& t, mpq& k,explanation & ex, unsigned inf_col, linear_combination_iterator<mpq>& iter);
lia_move report_conflict_from_gomory_cut(mpq & k);
lia_move report_gomory_cut(lar_term& t, mpq& k, mpq& lcm_den, unsigned num_ints);
void adjust_term_and_k_for_some_ints_case_gomory(lar_term& t, mpq& k, mpq& lcm_den);
void init_check_data();
bool constrain_free_vars(linear_combination_iterator<mpq> * r);
lia_move proceed_with_gomory_cut(lar_term& t, mpq& k, explanation& ex);
int find_next_free_var_in_gomory_row();
bool is_gomory_cut_target();
lia_move proceed_with_gomory_cut(lar_term& t, mpq& k, explanation& ex, unsigned j, linear_combination_iterator<mpq>& iter);
int find_free_var_in_gomory_row(linear_combination_iterator<mpq>& iter);
bool is_gomory_cut_target(linear_combination_iterator<mpq> &iter);
bool at_bound(unsigned j) const;
bool at_lower(unsigned j) const;
bool at_low(unsigned j) const;
bool at_upper(unsigned j) const;
bool has_low(unsigned j) const;
bool has_upper(unsigned j) const;
unsigned row_of_basic_column(unsigned j) const;
inline static bool is_rational(const impq & n) {
return is_zero(n.y);
}
inline static
mpq fractional_part(const impq & n) {
lp_assert(is_rational);
lp_assert(is_rational(n));
return n.x - floor(n.x);
}
void real_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& t, explanation & ex);
void int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& t, explanation& ex, mpq & lcm_den);
void real_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& t, explanation & ex, unsigned inf_column);
void int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& t, explanation& ex, mpq & lcm_den, unsigned inf_column);
constraint_index column_upper_bound_constraint(unsigned j) const;
constraint_index column_low_bound_constraint(unsigned j) const;
void display_row_info(std::ostream & out, unsigned row_index) const;
void gomory_cut_adjust_t_and_k_for_size_1(const vector<std::pair<mpq, unsigned>> & pol, lar_term & t, mpq &k);
void gomory_cut_adjust_t_and_k_for_size_gt_1(vector<std::pair<mpq, unsigned>> & pol, lar_term & t, mpq &k, unsigned num_ints, mpq &lcm_den);
void gomory_cut_adjust_t_and_k(vector<std::pair<mpq, unsigned>> & pol, lar_term & t, mpq &k, bool num_ints, mpq &lcm_den);
bool current_solution_is_inf_on_cut(const lar_term& t, const mpq& k) const;
public:
bool shift_var(unsigned j, unsigned range);
private:
unsigned random();
bool non_basic_columns_are_at_bounds() const;
bool has_inf_int() const;
lia_move create_branch_on_column(int j, lar_term& t, mpq& k) const;
};
}

4
src/util/lp/iterator_on_column.h
View file

@ -57,7 +57,9 @@ struct iterator_on_column:linear_combination_iterator<T> {
m_i = -1;
}
linear_combination_iterator<mpq> * clone() override {
bool is_reset() const { return m_i == -1;}
linear_combination_iterator<mpq> * clone() {
iterator_on_column * r = new iterator_on_column(m_column, m_A);
return r;
}

5
src/util/lp/iterator_on_indexed_vector.h
View file

@ -46,7 +46,10 @@ struct iterator_on_indexed_vector:linear_combination_iterator<T> {
void reset() override {
m_offset = 0;
}
linear_combination_iterator<T>* clone() override {
bool is_reset() const { return m_offset == 0;}
linear_combination_iterator<T>* clone() {
return new iterator_on_indexed_vector(m_v);
}
};

5
src/util/lp/iterator_on_pivot_row.h
View file

@ -51,7 +51,10 @@ struct iterator_on_pivot_row:linear_combination_iterator<T> {
m_basis_returned = false;
m_it.reset();
}
linear_combination_iterator<T> * clone() override {
bool is_reset() const { return m_basis_returned == false && m_it.is_reset();}
linear_combination_iterator<T> * clone() {
iterator_on_pivot_row * r = new iterator_on_pivot_row(m_v, m_basis_j);
return r;
}

5
src/util/lp/iterator_on_row.h
View file

@ -45,7 +45,10 @@ struct iterator_on_row:linear_combination_iterator<T> {
void reset() override {
m_i = 0;
}
linear_combination_iterator<T>* clone() override {
bool is_reset() const { return m_i == 0;}
linear_combination_iterator<T>* clone() {
return new iterator_on_row(m_row);
}
};

5
src/util/lp/iterator_on_term_with_basis_var.h
View file

@ -60,7 +60,10 @@ struct iterator_on_term_with_basis_var:linear_combination_iterator<mpq> {
m_i++;
return true;
}
void reset() override {
bool is_reset() const { return m_term_j_returned == false && m_i == m_term.m_coeffs.begin();}
void reset() {
m_term_j_returned = false;
m_i = m_term.m_coeffs.begin();
}

5
src/util/lp/lar_core_solver.h
View file

@ -155,7 +155,7 @@ public:
void fill_evidence(unsigned row);
unsigned get_number_of_non_ints() const;
void solve();
@ -344,8 +344,7 @@ public:
for (const auto & cc : m_r_solver.m_A.m_columns[j]){
unsigned i = cc.m_i;
unsigned jb = m_r_solver.m_basis[i];
m_r_solver.m_x[jb] -= delta * m_r_solver.m_A.get_val(cc);
m_r_solver.update_column_in_inf_set(jb);
m_r_solver.update_x_with_delta_and_track_feasibility(jb, - delta * m_r_solver.m_A.get_val(cc));
}
lp_assert(m_r_solver.A_mult_x_is_off() == false);
}

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

@ -255,14 +255,23 @@ void lar_core_solver::fill_not_improvable_zero_sum() {
}
}
unsigned lar_core_solver::get_number_of_non_ints() const {
unsigned n = 0;
for (auto & x : m_r_solver.m_x) {
if (x.is_int() == false)
n++;
}
return n;
}
void lar_core_solver::solve() {
lp_assert(m_r_solver.non_basic_columns_are_set_correctly());
lp_assert(m_r_solver.inf_set_is_correct());
if (m_r_solver.current_x_is_feasible() && m_r_solver.m_look_for_feasible_solution_only) {
m_r_solver.set_status(lp_status::OPTIMAL);
return;
}
TRACE("find_feas_stats", tout << "infeasibles = " << m_r_solver.m_inf_set.size() << ", int_infs = " << get_number_of_non_ints() << std::endl;);
if (m_r_solver.current_x_is_feasible() && m_r_solver.m_look_for_feasible_solution_only) {
m_r_solver.set_status(lp_status::OPTIMAL);
return;
}
++settings().st().m_need_to_solve_inf;
lp_assert(!m_r_solver.A_mult_x_is_off());
lp_assert((!settings().use_tableau()) || r_basis_is_OK());
@ -278,6 +287,7 @@ void lar_core_solver::solve() {
solve_on_signature_tableau(solution_signature, changes_of_basis);
else
solve_on_signature(solution_signature, changes_of_basis);
lp_assert(!settings().use_tableau() || r_basis_is_OK());
} else {
if (!settings().use_tableau()) {

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

@ -31,23 +31,22 @@ 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, const impq & x){
if (!var_is_int(j)) {
lp_assert(m_inf_int_set.contains(j) == false);
return;
}
if (m_mpq_lar_core_solver.m_r_x[j].is_int())
m_inf_int_set.erase(j);
else
m_inf_int_set.insert(j);
})
{
}
m_tracker_of_x_change([&](unsigned j) {
call_assignment_tracker(j);
}
),
m_int_solver(nullptr)
{}
void lar_solver::set_propagate_bounds_on_pivoted_rows_mode(bool v) {
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;
@ -55,8 +54,6 @@ lar_solver::~lar_solver(){
delete t;
}
numeric_pair<mpq> const& lar_solver::get_value(var_index vi) const { return m_mpq_lar_core_solver.m_r_x[vi]; }
bool lar_solver::is_term(var_index j) const {
return j >= m_terms_start_index && j - m_terms_start_index < m_terms.size();
}
@ -89,12 +86,6 @@ void lar_solver::print_implied_bound(const implied_bound& be, std::ostream & out
out << get_column_name(v);
}
out << " " << lconstraint_kind_string(be.kind()) << " " << be.m_bound << std::endl;
// for (auto & p : be.m_explanation) {
// out << p.first << " : ";
// print_constraint(p.second, out);
// }
// m_mpq_lar_core_solver.m_r_solver.print_column_info(be.m_j< m_terms_start_index? be.m_j : adjust_term_index(be.m_j), out);
out << "end of implied bound" << std::endl;
}
@ -294,6 +285,8 @@ bool lar_solver::has_int_var() const {
}
lp_status lar_solver::find_feasible_solution() {
TRACE("arith_int", );
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();
@ -307,10 +300,13 @@ lp_status lar_solver::find_feasible_solution() {
}
m_mpq_lar_core_solver.m_r_solver.m_look_for_feasible_solution_only = true;
return solve();
auto ret = solve();
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;
}
@ -321,6 +317,7 @@ lp_status lar_solver::solve() {
}
m_columns_with_changed_bound.clear();
lp_assert(inf_int_set_is_correct());
return m_status;
}
@ -368,6 +365,8 @@ void lar_solver::shrink_inf_set_after_pop(unsigned n, int_set & set) {
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<int>(m_ext_vars_to_columns.size());
@ -376,15 +375,23 @@ void lar_solver::pop(unsigned 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();
}
TRACE("arith_int", tout << "pop" << std::endl;);
lp_assert(inf_int_set_is_correct());
TRACE("arith_int", tout << "pop" << std::endl;);
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();
m_status = m_mpq_lar_core_solver.m_r_solver.current_x_is_feasible()? lp_status::OPTIMAL: lp_status::UNKNOWN;
lp_assert(m_settings.simplex_strategy() == simplex_strategy_enum::undecided ||
@ -655,15 +662,15 @@ bool lar_solver::costs_are_used() const {
return m_settings.simplex_strategy() != simplex_strategy_enum::tableau_rows;
}
void lar_solver::change_basic_x_by_delta_on_column(unsigned j, const numeric_pair<mpq> & delta) {
void lar_solver::change_basic_columns_dependend_on_a_given_nb_column(unsigned j, const numeric_pair<mpq> & 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];
m_mpq_lar_core_solver.m_r_x[bj] -= A_r().get_val(c) * delta;
if (tableau_with_costs()) {
m_basic_columns_with_changed_cost.insert(bj);
}
m_mpq_lar_core_solver.m_r_solver.update_column_in_inf_set(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;);
@ -675,26 +682,26 @@ void lar_solver::change_basic_x_by_delta_on_column(unsigned j, const numeric_pai
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_x[bj] -= m_column_buffer[i] * delta;
m_mpq_lar_core_solver.m_r_solver.update_column_in_inf_set(bj);
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.update_column_in_inf_set(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.update_column_in_inf_set(j);
m_mpq_lar_core_solver.m_r_solver.track_column_feasibility(j);
}
} else {
numeric_pair<mpq> delta;
if (m_mpq_lar_core_solver.m_r_solver.make_column_feasible(j, delta))
change_basic_x_by_delta_on_column(j, delta);
change_basic_columns_dependend_on_a_given_nb_column(j, delta);
}
}
@ -722,7 +729,6 @@ void lar_solver::update_x_and_inf_costs_for_columns_with_changed_bounds() {
}
void lar_solver::update_x_and_inf_costs_for_columns_with_changed_bounds_tableau() {
lp_assert(ax_is_correct());
for (auto j : m_columns_with_changed_bound.m_index)
update_x_and_inf_costs_for_column_with_changed_bounds(j);
@ -751,6 +757,7 @@ void lar_solver::solve_with_core_solver() {
update_x_and_inf_costs_for_columns_with_changed_bounds();
m_mpq_lar_core_solver.solve();
set_status(m_mpq_lar_core_solver.m_r_solver.get_status());
lp_assert(inf_int_set_is_correct());
lp_assert(m_status != lp_status::OPTIMAL || all_constraints_hold());
}
@ -928,6 +935,7 @@ bool lar_solver::all_constraints_hold() const {
}
bool lar_solver::constraint_holds(const lar_base_constraint & constr, std::unordered_map<var_index, mpq> & 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;
@ -1246,8 +1254,9 @@ void lar_solver::fill_var_set_for_random_update(unsigned sz, var_index const * v
void lar_solver::random_update(unsigned sz, var_index const * vars) {
vector<unsigned> column_list;
fill_var_set_for_random_update(sz, vars, column_list);
random_updater ru(m_mpq_lar_core_solver, column_list);
random_updater ru(*this, column_list);
ru.update();
lp_assert(inf_int_set_is_correct());
}
@ -1378,6 +1387,7 @@ void lar_solver::pop_tableau() {
}
void lar_solver::clean_inf_set_of_r_solver_after_pop() {
lp_assert(inf_int_set_is_correct());
vector<unsigned> became_feas;
clean_popped_elements(A_r().column_count(), m_mpq_lar_core_solver.m_r_solver.m_inf_set);
std::unordered_set<unsigned> basic_columns_with_changed_cost;
@ -1388,8 +1398,10 @@ void lar_solver::clean_inf_set_of_r_solver_after_pop() {
}
// some basic columns might become non-basic - these columns need to be made feasible
numeric_pair<mpq> delta;
if (m_mpq_lar_core_solver.m_r_solver.make_column_feasible(j, delta))
change_basic_x_by_delta_on_column(j, 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);
}
@ -1492,6 +1504,7 @@ var_index lar_solver::add_var(unsigned ext_j, bool is_int) {
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;
}
@ -1507,6 +1520,7 @@ 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);
@ -1574,7 +1588,7 @@ var_index lar_solver::add_term(const vector<std::pair<mpq, var_index>> & coeffs,
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_for_term(m_terms.back(), ret);
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);
}
@ -1582,12 +1596,6 @@ var_index lar_solver::add_term(const vector<std::pair<mpq, var_index>> & coeffs,
return ret;
}
void lar_solver::add_row_for_term(const lar_term * term, unsigned term_ext_index) {
lp_assert(sizes_are_correct());
add_row_from_term_no_constraint(term, term_ext_index);
lp_assert(sizes_are_correct());
}
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
@ -1604,9 +1612,10 @@ void lar_solver::add_row_from_term_no_constraint(const lar_term * term, unsigned
else {
fill_last_row_of_A_r(A_r(), term);
}
m_mpq_lar_core_solver.m_r_x[j] = get_basic_var_value_from_row_directly(A_r().row_count() - 1);
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() {
@ -1615,6 +1624,7 @@ void lar_solver::add_basic_var_to_core_fields() {
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);
@ -1639,6 +1649,7 @@ constraint_index lar_solver::add_var_bound(var_index j, lconstraint_kind kind, c
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;
}

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

@ -81,9 +81,18 @@ class lar_solver : public column_namer {
indexed_vector<mpq> m_column_buffer;
public:
lar_core_solver m_mpq_lar_core_solver;
private:
std::function<void (unsigned)> m_tracker_of_x_change;
int_solver * m_int_solver;
public:
void set_int_solver(int_solver * int_slv) {
m_int_solver = int_slv;
}
int_solver * get_int_solver() {
return m_int_solver;
}
unsigned constraint_count() const;
const lar_base_constraint& get_constraint(unsigned ci) const;
std::function<void (unsigned, const impq&)> m_tracker_of_x_change;
int_set m_inf_int_set;
////////////////// methods ////////////////////////////////
static_matrix<mpq, numeric_pair<mpq>> & A_r() { return m_mpq_lar_core_solver.m_r_A;}
@ -94,6 +103,14 @@ public:
static bool valid_index(unsigned j){ return static_cast<int>(j) >= 0;}
bool column_is_int(unsigned j) const;
bool column_value_is_int(unsigned j) const {
return m_mpq_lar_core_solver.m_r_x[j].is_int();
}
const impq& get_column_value(unsigned j) const {
return m_mpq_lar_core_solver.m_r_x[j];
}
bool is_term(var_index j) const;
bool column_is_fixed(unsigned j) const;
public:
@ -134,87 +151,13 @@ public:
void clear();
lar_solver();
void set_propagate_bounds_on_pivoted_rows_mode(bool v);
void set_track_pivoted_rows(bool v);
bool get_track_pivoted_rows() const;
virtual ~lar_solver();
void 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_orig_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;
// for (auto & p : be.m_explanation) {
// out << p.first << " : ";
// print_constraint(p.second, out);
// }
// m_mpq_lar_core_solver.m_r_solver.print_column_info(be.m_j< m_terms_start_index? be.m_j : adjust_term_index(be.m_j), out);
out << "end of implied bound" << std::endl;
}
bool implied_bound_is_correctly_explained(implied_bound const & be, const vector<std::pair<mpq, unsigned>> & explanation) const {
std::unordered_map<unsigned, mpq> coeff_map;
auto rs_of_evidence = zero_of_type<mpq>();
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;
}
SASSERT(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<lconstraint_kind>((static_cast<int>(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<mpq>()) {
kind = static_cast<lconstraint_kind>(-kind);
}
rs_of_evidence /= ratio;
} else {
const lar_term * t = m_orig_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<mpq>()) {
kind = static_cast<lconstraint_kind>(-kind);
}
rs_of_evidence /= ratio;
rs_of_evidence += t->m_v * ratio;
}
return kind == be.kind() && rs_of_evidence == be.m_bound;
}
unsigned adjust_term_index(unsigned j) const;
void analyze_new_bounds_on_row(
unsigned row_index,
@ -707,32 +650,9 @@ public:
bool tableau_with_costs() const;
bool costs_are_used() const {
return m_settings.simplex_strategy() != simplex_strategy_enum::tableau_rows;
}
void change_basic_x_by_delta_on_column(unsigned j, const numeric_pair<mpq> & delta) {
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];
m_mpq_lar_core_solver.m_r_x[bj] -= A_r().get_val(c) * delta;
if (tableau_with_costs()) {
m_basic_columns_with_changed_cost.insert(bj);
}
m_mpq_lar_core_solver.m_r_solver.update_column_in_inf_set(bj);
}
} 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_x[bj] -= m_column_buffer[i] * delta;
m_mpq_lar_core_solver.m_r_solver.update_column_in_inf_set(bj);
}
}
}
bool costs_are_used() const;
void change_basic_columns_dependend_on_a_given_nb_column(unsigned j, const numeric_pair<mpq> & delta);
void update_x_and_inf_costs_for_column_with_changed_bounds(unsigned j);
void detect_rows_with_changed_bounds_for_column(unsigned j) {
@ -1261,8 +1181,8 @@ public:
void clean_inf_set_of_r_solver_after_pop();
void shrink_explanation_to_minimum(vector<std::pair<mpq, constraint_index>> & explanation) const;
bool column_represents_row_in_tableau(unsigned j) {
return m_vars_to_ul_pairs()[j].m_i != static_cast<row_index>(-1);
bool column_value_is_integer(unsigned j) const {
return get_column_value(j).is_int();
}
bool column_is_real(unsigned j) const {
@ -1442,6 +1362,43 @@ bool model_is_int_feasible() const;
t = lar_term(pol_after_subs, v);
}
bool inf_int_set_is_correct_for_column(unsigned j) const {
if (m_inf_int_set.contains(j) != (column_is_int(j) && (!column_value_is_integer(j)))) {
TRACE("arith_int",
tout << "j= " << j <<
" inf_int_set().contains(j) = " << m_inf_int_set.contains(j) <<
", column_is_int(j) = " << column_is_int(j) <<
"\n column_value_is_integer(j) = " << column_value_is_integer(j) <<
", val = " << get_column_value(j) << std::endl;);
return false;
}
return true;
}
bool inf_int_set_is_correct() const {
if (!has_int_var())
return true;
for (unsigned j = 0; j < A_r().column_count(); j++) {
if (inf_int_set_is_correct_for_column(j) == false)
return false;
}
return true;
}
bool has_int_var() const;
void call_assignment_tracker(unsigned j) {
if (!var_is_int(j)) {
lp_assert(m_inf_int_set.contains(j) == false);
return;
}
if (m_mpq_lar_core_solver.m_r_x[j].is_int())
m_inf_int_set.erase(j);
else
m_inf_int_set.insert(j);
}
lar_core_solver & get_core_solver() { return m_mpq_lar_core_solver; }
};
}

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

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

3
src/util/lp/linear_combination_iterator.h
View file

@ -27,6 +27,7 @@ struct linear_combination_iterator {
virtual linear_combination_iterator * clone() = 0;
virtual ~linear_combination_iterator(){}
virtual unsigned size() const = 0;
virtual bool is_reset() const = 0;
};
template <typename T>
struct linear_combination_iterator_on_vector : linear_combination_iterator<T> {
@ -50,7 +51,7 @@ struct linear_combination_iterator_on_vector : linear_combination_iterator<T> {
m_offset++;
return true;
}
bool is_reset() const { return m_offset == 0;}
void reset() {m_offset = 0;}
linear_combination_iterator<T> * clone() {
return new linear_combination_iterator_on_vector(m_vector);

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

@ -32,7 +32,7 @@ public:
column_type get_column_type(unsigned) const;
const impq & get_low_bound(unsigned) const;
const impq & get_upper_bound(unsigned) const;
void try_add_bound(const mpq & v, unsigned j, bool is_low, bool coeff_before_j_is_pos, unsigned row_or_term_index, bool strict);
void try_add_bound(mpq v, unsigned j, bool is_low, bool coeff_before_j_is_pos, unsigned row_or_term_index, bool strict);
virtual bool bound_is_interesting(unsigned vi,
lp::lconstraint_kind kind,
const rational & bval) {return true;}

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

@ -45,6 +45,9 @@ public:
TRACE("feas",
if (m_inf_set.size()) {
tout << "column " << m_inf_set.m_index[0] << " is infeasible" << std::endl;
print_column_info(m_inf_set.m_index[0], tout);
} else {
tout << "x is feasible\n";
}
);
return m_inf_set.size() == 0;
@ -84,9 +87,9 @@ public:
bool m_tracing_basis_changes;
int_set* m_pivoted_rows;
bool m_look_for_feasible_solution_only;
std::function<void (unsigned, const X &)> * m_tracker_of_x_change;
std::function<void (unsigned)> * m_tracker_of_x_change;
void set_tracker_of_x(std::function<void (unsigned, const X&)>* tracker) {
void set_tracker_of_x(std::function<void (unsigned)>* tracker) {
m_tracker_of_x_change = tracker;
}
@ -401,6 +404,7 @@ public:
}
bool make_column_feasible(unsigned j, numeric_pair<mpq> & delta) {
bool ret = false;
lp_assert(m_basis_heading[j] < 0);
auto & x = m_x[j];
switch (m_column_types[j]) {
@ -408,43 +412,39 @@ public:
lp_assert(m_low_bounds[j] == m_upper_bounds[j]);
if (x != m_low_bounds[j]) {
delta = m_low_bounds[j] - x;
x = m_low_bounds[j];
return true;
ret = true;;
}
break;
case column_type::boxed:
if (x < m_low_bounds[j]) {
delta = m_low_bounds[j] - x;
x = m_low_bounds[j];
return true;
ret = true;;
}
if (x > m_upper_bounds[j]) {
delta = m_upper_bounds[j] - x;
x = m_upper_bounds[j];
return true;
ret = true;
}
break;
case column_type::low_bound:
if (x < m_low_bounds[j]) {
delta = m_low_bounds[j] - x;
x = m_low_bounds[j];
return true;
ret = true;
}
break;
case column_type::upper_bound:
if (x > m_upper_bounds[j]) {
delta = m_upper_bounds[j] - x;
x = m_upper_bounds[j];
return true;
ret = true;
}
break;
case column_type::free_column:
break;
default:
lp_assert(false);
break;
}
return false;
if (ret)
add_delta_to_x_and_call_tracker(j, delta);
return ret;
}
@ -584,25 +584,26 @@ public:
}
void print_column_info(unsigned j, std::ostream & out) const {
out << "column_index = " << j << ", name = "<< column_name(j) << " type = " << column_type_to_string(m_column_types[j]) << std::endl;
out << "j = " << j << ", name = "<< column_name(j);
switch (m_column_types[j]) {
case column_type::fixed:
case column_type::boxed:
out << "(" << m_low_bounds[j] << ", " << m_upper_bounds[j] << ")" << std::endl;
out << " [" << m_low_bounds[j] << ", " << m_upper_bounds[j] << "]";
break;
case column_type::low_bound:
out << m_low_bounds[j] << std::endl;
out << " [" << m_low_bounds[j] << "," << "oo" << "]";
break;
case column_type::upper_bound:
out << m_upper_bounds[j] << std::endl;
out << " [-oo, " << m_upper_bounds[j] << ']';
break;
case column_type::free_column:
out << " [-oo, oo]";
break;
default:
lp_assert(false);
}
out << "basis heading = " << m_basis_heading[j] << std::endl;
out << "x = " << m_x[j] << std::endl;
// out << "basis heading = " << m_basis_heading[j] << std::endl;
out << " x = " << m_x[j] << std::endl;
/*
std::cout << "cost = " << m_costs[j] << std::endl;
std:: cout << "m_d = " << m_d[j] << std::endl;*/
@ -689,22 +690,43 @@ public:
return m_inf_set.contains(j);
}
void update_column_in_inf_set(unsigned j) {
if (column_is_feasible(j)) {
void update_x_with_feasibility_tracking(unsigned j, const X & v) {
m_x[j] = v;
track_column_feasibility(j);
}
void update_x_with_delta_and_track_feasibility(unsigned j, const X & del) {
m_x[j] += del;
track_column_feasibility(j);
}
void update_x_and_call_tracker(unsigned j, const X & v) {
m_x[j] = v;
if (m_tracker_of_x_change != nullptr)
(*m_tracker_of_x_change)(j);
}
void add_delta_to_x_and_call_tracker(unsigned j, const X & delta) {
m_x[j] += delta;
if (m_tracker_of_x_change != nullptr)
(*m_tracker_of_x_change)(j);
}
void track_column_feasibility(unsigned j) {
if (column_is_feasible(j))
remove_column_from_inf_set(j);
} else {
else
insert_column_into_inf_set(j);
}
}
void insert_column_into_inf_set(unsigned j) {
if (m_tracker_of_x_change != nullptr)
(*m_tracker_of_x_change)(j, m_x[j]);
(*m_tracker_of_x_change)(j);
m_inf_set.insert(j);
lp_assert(!column_is_feasible(j));
}
void remove_column_from_inf_set(unsigned j) {
if (m_tracker_of_x_change != nullptr)
(*m_tracker_of_x_change)(j, m_x[j]);
(*m_tracker_of_x_change)(j);
m_inf_set.erase(j);
lp_assert(column_is_feasible(j));
}

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

@ -445,7 +445,7 @@ public:
void advance_on_entering_and_leaving_tableau_rows(int entering, int leaving, const X &theta ) {
this->update_basis_and_x_tableau(entering, leaving, theta);
this->update_column_in_inf_set(entering);
this->track_column_feasibility(entering);
}

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

@ -864,7 +864,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::print_column
template <typename T, typename X> unsigned lp_primal_core_solver<T, X>::solve() {
if (numeric_traits<T>::precise() && this->m_settings.use_tableau())
return solve_with_tableau();
init_run();
if (this->current_x_is_feasible() && this->m_look_for_feasible_solution_only) {
this->set_status(lp_status::FEASIBLE);
@ -880,6 +880,7 @@ template <typename T, typename X> unsigned lp_primal_core_solver<T, X>::solve()
return this->total_iterations();
}
one_iteration();
lp_assert(!this->m_using_infeas_costs || this->costs_on_nbasis_are_zeros());
switch (this->get_status()) {
case lp_status::OPTIMAL: // double check that we are at optimum

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

@ -112,8 +112,9 @@ unsigned lp_primal_core_solver<T, X>::solve_with_tableau() {
if (this->print_statistics_with_iterations_and_nonzeroes_and_cost_and_check_that_the_time_is_over((this->m_using_infeas_costs? "inf t" : "feas t"), * this->m_settings.get_message_ostream())) {
return this->total_iterations();
}
if (this->m_settings.use_tableau_rows())
if (this->m_settings.use_tableau_rows()) {
one_iteration_tableau_rows();
}
else
one_iteration_tableau();
switch (this->get_status()) {
@ -352,22 +353,20 @@ update_basis_and_x_tableau(int entering, int leaving, X const & tt) {
}
template <typename T, typename X> void lp_primal_core_solver<T, X>::
update_x_tableau(unsigned entering, const X& delta) {
this->add_delta_to_x_and_call_tracker(entering, delta);
if (!this->m_using_infeas_costs) {
this->m_x[entering] += delta;
for (const auto & c : this->m_A.m_columns[entering]) {
unsigned i = c.m_i;
this->m_x[this->m_basis[i]] -= delta * this->m_A.get_val(c);
this->update_column_in_inf_set(this->m_basis[i]);
this->update_x_with_delta_and_track_feasibility(this->m_basis[i], - delta * this->m_A.get_val(c));
}
} else { // m_using_infeas_costs == true
this->m_x[entering] += delta;
lp_assert(this->column_is_feasible(entering));
lp_assert(this->m_costs[entering] == zero_of_type<T>());
// m_d[entering] can change because of the cost change for basic columns.
for (const auto & c : this->m_A.m_columns[entering]) {
unsigned i = c.m_i;
unsigned j = this->m_basis[i];
this->m_x[j] -= delta * this->m_A.get_val(c);
this->add_delta_to_x_and_call_tracker(j, -delta * this->m_A.get_val(c));
update_inf_cost_for_column_tableau(j);
if (is_zero(this->m_costs[j]))
this->remove_column_from_inf_set(j);

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

@ -221,7 +221,8 @@ public:
max_row_length_for_bound_propagation(300),
backup_costs(true),
column_number_threshold_for_using_lu_in_lar_solver(4000),
m_int_branch_cut_threshold(10000000)
m_int_branch_cut_gomory_threshold(4),
m_run_gcd_test(true)
{}
void set_resource_limit(lp_resource_limit& lim) { m_resource_limit = &lim; }
@ -324,11 +325,12 @@ public:
#endif
bool use_breakpoints_in_feasibility_search;
unsigned random_next() { return m_rand(); }
void random_seed(unsigned s) { m_rand.set_seed(s); }
void set_random_seed(unsigned s) { m_rand.set_seed(s); }
unsigned max_row_length_for_bound_propagation;
bool backup_costs;
unsigned column_number_threshold_for_using_lu_in_lar_solver;
unsigned m_int_branch_cut_threshold;
unsigned m_int_branch_cut_gomory_threshold;
bool m_run_gcd_test;
}; // end of lp_settings class
@ -351,7 +353,7 @@ inline std::string T_to_string(const numeric_pair<mpq> & t) {
inline std::string T_to_string(const mpq & t) {
std::ostringstream strs;
strs << t.get_double();
strs << t;
return strs.str();
}

64
src/util/lp/random_updater.h
View file

@ -30,67 +30,19 @@ Revision History:
// The class searches for a feasible solution with as many different values of variables as it can find
namespace lp {
template <typename T> struct numeric_pair; // forward definition
class lar_core_solver; // forward definition
class lar_solver; // forward definition
class random_updater {
struct interval {
bool upper_bound_is_set;
numeric_pair<mpq> upper_bound;
bool low_bound_is_set;
numeric_pair<mpq> low_bound;
interval() : upper_bound_is_set(false),
low_bound_is_set(false) {}
void set_low_bound(const numeric_pair<mpq> & v) {
if (low_bound_is_set) {
low_bound = std::max(v, low_bound);
} else {
low_bound = v;
low_bound_is_set = true;
}
}
void set_upper_bound(const numeric_pair<mpq> & v) {
if (upper_bound_is_set) {
upper_bound = std::min(v, upper_bound);
} else {
upper_bound = v;
upper_bound_is_set = true;
}
}
bool is_empty() const {return
upper_bound_is_set && low_bound_is_set && low_bound >= upper_bound;
}
bool low_bound_holds(const numeric_pair<mpq> & a) const {
return low_bound_is_set == false || a >= low_bound;
}
bool upper_bound_holds(const numeric_pair<mpq> & a) const {
return upper_bound_is_set == false || a <= upper_bound;
}
bool contains(const numeric_pair<mpq> & a) const {
return low_bound_holds(a) && upper_bound_holds(a);
}
std::string lbs() { return low_bound_is_set ? T_to_string(low_bound):std::string("inf");}
std::string rbs() { return upper_bound_is_set? T_to_string(upper_bound):std::string("inf");}
std::string to_str() { return std::string("[")+ lbs() + ", " + rbs() + "]";}
};
std::set<var_index> m_var_set;
lar_core_solver & m_core_solver;
unsigned range;
linear_combination_iterator<mpq>* m_column_j; // the actual column
interval find_shift_interval(unsigned j);
interval get_interval_of_non_basic_var(unsigned j);
lar_solver & m_lar_solver;
unsigned m_range;
void add_column_to_sets(unsigned j);
void random_shift_var(unsigned j);
bool random_shift_var(unsigned j);
std::unordered_map<numeric_pair<mpq>, unsigned> m_values; // it maps a value to the number of time it occurs
void diminish_interval_to_leave_basic_vars_feasible(numeric_pair<mpq> &nb_x, interval & inter);
void shift_var(unsigned j, interval & r);
void diminish_interval_for_basic_var(numeric_pair<mpq> &nb_x, unsigned j, mpq & a, interval & r);
numeric_pair<mpq> get_random_from_interval(interval & r);
void add_value(numeric_pair<mpq>& v);
void remove_value(numeric_pair<mpq> & v);
bool shift_var(unsigned j);
void add_value(const numeric_pair<mpq>& v);
void remove_value(const numeric_pair<mpq> & v);
public:
random_updater(lar_core_solver & core_solver, const vector<unsigned> & column_list);
random_updater(lar_solver & solver, const vector<unsigned> & column_list);
void update();
};
}

173
src/util/lp/random_updater.hpp
View file

@ -26,156 +26,25 @@ namespace lp {
random_updater::random_updater(
lar_core_solver & lar_core_solver,
lar_solver & lar_solver,
const vector<unsigned> & column_indices) :
m_core_solver(lar_core_solver),
range(100000) {
m_lar_solver(lar_solver),
m_range(100000) {
for (unsigned j : column_indices)
add_column_to_sets(j);
}
random_updater::interval random_updater::get_interval_of_non_basic_var(unsigned j) {
interval ret;
switch (m_core_solver.get_column_type(j)) {
case column_type::free_column:
break;
case column_type::low_bound:
ret.set_low_bound(m_core_solver.m_r_low_bounds[j]);
break;
case column_type::upper_bound:
ret.set_upper_bound(m_core_solver.m_r_upper_bounds[j]);
break;
case column_type::boxed:
case column_type::fixed:
ret.set_low_bound(m_core_solver.m_r_low_bounds[j]);
ret.set_upper_bound(m_core_solver.m_r_upper_bounds[j]);
break;
default:
lp_assert(false);
}
return ret;
bool random_updater::shift_var(unsigned v) {
return m_lar_solver.get_int_solver()->shift_var(v, m_range);
}
void random_updater::diminish_interval_for_basic_var(numeric_pair<mpq>& nb_x, unsigned j,
mpq & a,
interval & r) {
lp_assert(m_core_solver.m_r_heading[j] >= 0);
numeric_pair<mpq> delta;
lp_assert(a != zero_of_type<mpq>());
switch (m_core_solver.get_column_type(j)) {
case column_type::free_column:
break;
case column_type::low_bound:
delta = m_core_solver.m_r_x[j] - m_core_solver.m_r_low_bounds[j];
lp_assert(delta >= zero_of_type<numeric_pair<mpq>>());
if (a > 0) {
r.set_upper_bound(nb_x + delta / a);
} else {
r.set_low_bound(nb_x + delta / a);
}
break;
case column_type::upper_bound:
delta = m_core_solver.m_r_upper_bounds()[j] - m_core_solver.m_r_x[j];
lp_assert(delta >= zero_of_type<numeric_pair<mpq>>());
if (a > 0) {
r.set_low_bound(nb_x - delta / a);
} else {
r.set_upper_bound(nb_x - delta / a);
}
break;
case column_type::boxed:
if (a > 0) {
delta = m_core_solver.m_r_x[j] - m_core_solver.m_r_low_bounds[j];
lp_assert(delta >= zero_of_type<numeric_pair<mpq>>());
r.set_upper_bound(nb_x + delta / a);
delta = m_core_solver.m_r_upper_bounds()[j] - m_core_solver.m_r_x[j];
lp_assert(delta >= zero_of_type<numeric_pair<mpq>>());
r.set_low_bound(nb_x - delta / a);
} else { // a < 0
delta = m_core_solver.m_r_upper_bounds()[j] - m_core_solver.m_r_x[j];
lp_assert(delta >= zero_of_type<numeric_pair<mpq>>());
r.set_upper_bound(nb_x - delta / a);
delta = m_core_solver.m_r_x[j] - m_core_solver.m_r_low_bounds[j];
lp_assert(delta >= zero_of_type<numeric_pair<mpq>>());
r.set_low_bound(nb_x + delta / a);
}
break;
case column_type::fixed:
r.set_low_bound(nb_x);
r.set_upper_bound(nb_x);
break;
default:
lp_assert(false);
}
}
void random_updater::diminish_interval_to_leave_basic_vars_feasible(numeric_pair<mpq> &nb_x, interval & r) {
m_column_j->reset();
unsigned i;
mpq a;
while (m_column_j->next(a, i)) {
diminish_interval_for_basic_var(nb_x, m_core_solver.m_r_basis[i], a, r);
if (r.is_empty())
break;
}
}
random_updater::interval random_updater::find_shift_interval(unsigned j) {
interval ret = get_interval_of_non_basic_var(j);
diminish_interval_to_leave_basic_vars_feasible(m_core_solver.m_r_x[j], ret);
return ret;
}
void random_updater::shift_var(unsigned j, interval & r) {
lp_assert(r.contains(m_core_solver.m_r_x[j]));
lp_assert(m_core_solver.m_r_solver.column_is_feasible(j));
auto old_x = m_core_solver.m_r_x[j];
remove_value(old_x);
auto new_val = m_core_solver.m_r_x[j] = get_random_from_interval(r);
add_value(new_val);
lp_assert(r.contains(m_core_solver.m_r_x[j]));
lp_assert(m_core_solver.m_r_solver.column_is_feasible(j));
auto delta = m_core_solver.m_r_x[j] - old_x;
unsigned i;
m_column_j->reset();
mpq a;
while(m_column_j->next(a, i)) {
unsigned bj = m_core_solver.m_r_basis[i];
m_core_solver.m_r_x[bj] -= a * delta;
lp_assert(m_core_solver.m_r_solver.column_is_feasible(bj));
}
lp_assert(m_core_solver.m_r_solver.A_mult_x_is_off() == false);
}
numeric_pair<mpq> random_updater::get_random_from_interval(interval & r) {
unsigned rand = m_core_solver.settings().random_next();
if ((!r.low_bound_is_set) && (!r.upper_bound_is_set))
return numeric_pair<mpq>(rand % range, 0);
if (r.low_bound_is_set && (!r.upper_bound_is_set))
return r.low_bound + numeric_pair<mpq>(rand % range, 0);
if ((!r.low_bound_is_set) && r.upper_bound_is_set)
return r.upper_bound - numeric_pair<mpq>(rand % range, 0);
lp_assert(r.low_bound_is_set && r.upper_bound_is_set);
return r.low_bound + (rand % range) * (r.upper_bound - r.low_bound)/ range;
}
void random_updater::random_shift_var(unsigned j) {
m_column_j = m_core_solver.get_column_iterator(j);
if (m_column_j->size() >= 50) {
delete m_column_j;
return;
}
interval interv = find_shift_interval(j);
if (interv.is_empty()) {
delete m_column_j;
return;
bool random_updater::random_shift_var(unsigned j) {
if (m_lar_solver.A_r().m_columns.size() >= 50) {
return false;
}
shift_var(j, interv);
delete m_column_j;
return shift_var(j);
}
void random_updater::update() {
@ -183,11 +52,15 @@ void random_updater::update() {
if (m_var_set.size() <= m_values.size()) {
break; // we are done
}
random_shift_var(j);
auto old_x = m_lar_solver.get_column_value(j);
if (random_shift_var(j)) {
remove_value(old_x);
add_value(m_lar_solver.get_column_value(j));
}
}
}
void random_updater::add_value(numeric_pair<mpq>& v) {
void random_updater::add_value(const numeric_pair<mpq>& v) {
auto it = m_values.find(v);
if (it == m_values.end()) {
m_values[v] = 1;
@ -196,7 +69,7 @@ void random_updater::add_value(numeric_pair<mpq>& v) {
}
}
void random_updater::remove_value(numeric_pair<mpq>& v) {
void random_updater::remove_value(const numeric_pair<mpq>& v) {
std::unordered_map<numeric_pair<mpq>, unsigned>::iterator it = m_values.find(v);
lp_assert(it != m_values.end());
it->second--;
@ -205,16 +78,16 @@ void random_updater::remove_value(numeric_pair<mpq>& v) {
}
void random_updater::add_column_to_sets(unsigned j) {
if (m_core_solver.m_r_heading[j] < 0) {
if (m_lar_solver.get_core_solver().m_r_heading[j] < 0) {
m_var_set.insert(j);
add_value(m_core_solver.m_r_x[j]);
add_value(m_lar_solver.get_core_solver().m_r_x[j]);
} else {
unsigned row = m_core_solver.m_r_heading[j];
for (auto row_c : m_core_solver.m_r_A.m_rows[row]) {
unsigned row = m_lar_solver.get_core_solver().m_r_heading[j];
for (auto row_c : m_lar_solver.get_core_solver().m_r_A.m_rows[row]) {
unsigned cj = row_c.m_j;
if (m_core_solver.m_r_heading[cj] < 0) {
if (m_lar_solver.get_core_solver().m_r_heading[cj] < 0) {
m_var_set.insert(cj);
add_value(m_core_solver.m_r_x[cj]);
add_value(m_lar_solver.get_core_solver().m_r_x[cj]);
}
}
}