diff --git a/src/math/lp/lar_core_solver.h b/src/math/lp/lar_core_solver.h index 4a5cb20b1..d58cbb9f7 100644 --- a/src/math/lp/lar_core_solver.h +++ b/src/math/lp/lar_core_solver.h @@ -123,6 +123,7 @@ public: m_stacked_simplex_strategy.pop(k); m_r_solver.m_settings.set_simplex_strategy(m_stacked_simplex_strategy); + m_infeasible_linear_combination.reset(); lp_assert(m_r_solver.basis_heading_is_correct()); } diff --git a/src/math/lp/lar_core_solver_def.h b/src/math/lp/lar_core_solver_def.h index b6743873f..e0ffed16e 100644 --- a/src/math/lp/lar_core_solver_def.h +++ b/src/math/lp/lar_core_solver_def.h @@ -49,7 +49,7 @@ void lar_core_solver::fill_not_improvable_zero_sum_from_inf_row() { m_infeasible_sum_sign = m_r_solver.inf_sign_of_column(bj); m_infeasible_linear_combination.clear(); for (auto & rc : m_r_solver.m_A.m_rows[m_r_solver.m_inf_row_index_for_tableau]) - m_infeasible_linear_combination.push_back(std::make_pair(rc.coeff(), rc.var())); + m_infeasible_linear_combination.push_back({rc.coeff(), rc.var()}); } void lar_core_solver::fill_not_improvable_zero_sum() { diff --git a/src/math/lp/lar_solver.cpp b/src/math/lp/lar_solver.cpp index 069709126..4c26aaff6 100644 --- a/src/math/lp/lar_solver.cpp +++ b/src/math/lp/lar_solver.cpp @@ -179,7 +179,7 @@ namespace lp { lp_status lar_solver::get_status() const { return m_status; } - void lar_solver::set_status(lp_status s) { + void lar_solver::set_status(lp_status s) { TRACE("lar_solver", tout << "setting status to " << s << "\n";); m_status = s; } @@ -224,9 +224,9 @@ namespace lp { } void lar_solver::push() { + m_trail.push_scope(); m_simplex_strategy = m_settings.simplex_strategy(); m_simplex_strategy.push(); - m_columns_to_ul_pairs.push(); m_crossed_bounds_column = null_lpvar; m_crossed_bounds_deps = nullptr; m_mpq_lar_core_solver.push(); @@ -260,16 +260,16 @@ namespace lp { TRACE("lar_solver", tout << "k = " << k << std::endl;); m_crossed_bounds_column = null_lpvar; m_crossed_bounds_deps = nullptr; - unsigned n = m_columns_to_ul_pairs.peek_size(k); + m_trail.pop_scope(k); + unsigned n = m_columns_to_ul_pairs.size(); m_var_register.shrink(n); - pop_tableau(n); + + lp_assert(m_mpq_lar_core_solver.m_r_solver.m_costs.size() == A_r().column_count()); + lp_assert(m_mpq_lar_core_solver.m_r_solver.m_basis.size() == A_r().row_count()); + lp_assert(m_mpq_lar_core_solver.m_r_solver.basis_heading_is_correct()); + lp_assert(A_r().column_count() == n); - TRACE("lar_solver_details", - for (unsigned j = 0; j < n; j++) { - print_column_info(j, tout) << "\n"; - } - ); - m_columns_to_ul_pairs.pop(k); + TRACE("lar_solver_details", for (unsigned j = 0; j < n; j++) print_column_info(j, tout) << "\n";); m_mpq_lar_core_solver.pop(k); remove_non_fixed_from_fixed_var_table(); @@ -612,26 +612,23 @@ namespace lp { } + void lar_solver::set_upper_bound_witness(var_index j, u_dependency* dep) { - ul_pair ul = m_columns_to_ul_pairs[j]; - ul.upper_bound_witness() = dep; - m_columns_to_ul_pairs[j] = ul; + m_trail.push(vector_value_trail(m_columns_to_ul_pairs, j)); + m_columns_to_ul_pairs[j].upper_bound_witness() = dep; } void lar_solver::set_lower_bound_witness(var_index j, u_dependency* dep) { - ul_pair ul = m_columns_to_ul_pairs[j]; - ul.lower_bound_witness() = dep; - m_columns_to_ul_pairs[j] = ul; + m_trail.push(vector_value_trail(m_columns_to_ul_pairs, j)); + m_columns_to_ul_pairs[j].lower_bound_witness() = dep; } void lar_solver::register_monoid_in_map(std::unordered_map& coeffs, const mpq& a, unsigned j) { auto it = coeffs.find(j); - if (it == coeffs.end()) { + if (it == coeffs.end()) coeffs[j] = a; - } - else { + else it->second += a; - } } @@ -1269,7 +1266,7 @@ namespace lp { } bool lar_solver::column_represents_row_in_tableau(unsigned j) { - return m_columns_to_ul_pairs()[j].associated_with_row(); + return m_columns_to_ul_pairs[j].associated_with_row(); } void lar_solver::make_sure_that_the_bottom_right_elem_not_zero_in_tableau(unsigned i, unsigned j) { @@ -1371,20 +1368,6 @@ namespace lp { lp_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_costs.size()); } - void lar_solver::pop_tableau(unsigned old_size) { - lp_assert(m_mpq_lar_core_solver.m_r_solver.m_costs.size() == A_r().column_count()); - - lp_assert(m_mpq_lar_core_solver.m_r_solver.m_basis.size() == A_r().row_count()); - lp_assert(m_mpq_lar_core_solver.m_r_solver.basis_heading_is_correct()); - // We remove last variables starting from m_column_names.size() to m_vec_of_canonic_left_sides.size(). - // At this moment m_column_names is already popped - - while (A_r().column_count() > old_size) - remove_last_column_from_tableau(); - lp_assert(m_mpq_lar_core_solver.m_r_solver.m_costs.size() == A_r().column_count()); - lp_assert(m_mpq_lar_core_solver.m_r_solver.m_basis.size() == A_r().row_count()); - lp_assert(m_mpq_lar_core_solver.m_r_solver.basis_heading_is_correct()); - } void lar_solver::clean_inf_heap_of_r_solver_after_pop() { vector became_feas; @@ -1490,6 +1473,15 @@ namespace lp { return j; } + struct lar_solver::add_column : public trail { + lar_solver& s; + add_column(lar_solver& s) : s(s) {} + virtual void undo() { + s.remove_last_column_from_tableau(); + s.m_columns_to_ul_pairs.pop_back(); + } + }; + var_index lar_solver::add_var(unsigned ext_j, bool is_int) { TRACE("add_var", tout << "adding var " << ext_j << (is_int ? " int" : " nonint") << std::endl;); var_index local_j; @@ -1500,9 +1492,9 @@ namespace lp { lp_assert(m_columns_to_ul_pairs.size() == A_r().column_count()); local_j = A_r().column_count(); m_columns_to_ul_pairs.push_back(ul_pair(false)); // not associated with a row - while (m_usage_in_terms.size() <= ext_j) { + m_trail.push(add_column(*this)); + while (m_usage_in_terms.size() <= ext_j) m_usage_in_terms.push_back(0); - } add_non_basic_var_to_core_fields(ext_j, is_int); lp_assert(sizes_are_correct()); return local_j; @@ -1653,6 +1645,7 @@ namespace lp { unsigned j = A_r().column_count(); ul_pair ul(true); // to mark this column as associated_with_row m_columns_to_ul_pairs.push_back(ul); + m_trail.push(add_column(*this)); add_basic_var_to_core_fields(); A_r().fill_last_row_with_pivoting(*term, @@ -2367,7 +2360,7 @@ namespace lp { SASSERT(m_crossed_bounds_deps == nullptr); set_status(lp_status::INFEASIBLE); m_crossed_bounds_column = j; - const auto& ul = this->m_columns_to_ul_pairs()[j]; + const auto& ul = this->m_columns_to_ul_pairs[j]; u_dependency* bdep = lower_bound? ul.lower_bound_witness() : ul.upper_bound_witness(); SASSERT(bdep != nullptr); m_crossed_bounds_deps = m_dependencies.mk_join(bdep, dep); diff --git a/src/math/lp/lar_solver.h b/src/math/lp/lar_solver.h index 949cb0910..6c7b30664 100644 --- a/src/math/lp/lar_solver.h +++ b/src/math/lp/lar_solver.h @@ -42,6 +42,7 @@ #include "util/debug.h" #include "util/stacked_value.h" #include "util/vector.h" +#include "util/trail.h" namespace lp { @@ -74,6 +75,7 @@ class lar_solver : public column_namer { }; //////////////////// fields ////////////////////////// + trail_stack m_trail; lp_settings m_settings; lp_status m_status = lp_status::UNKNOWN; stacked_value m_simplex_strategy; @@ -85,7 +87,8 @@ class lar_solver : public column_namer { bool m_need_register_terms = false; var_register m_var_register; var_register m_term_register; - stacked_vector m_columns_to_ul_pairs; + struct add_column; + svector m_columns_to_ul_pairs; constraint_set m_constraints; // the set of column indices j such that bounds have changed for j indexed_uint_set m_columns_with_changed_bounds; @@ -240,7 +243,6 @@ class lar_solver : public column_namer { void remove_last_column_from_basis_tableau(unsigned j); void remove_last_column_from_tableau(); - void pop_tableau(unsigned old_size); void clean_inf_heap_of_r_solver_after_pop(); inline bool column_value_is_integer(unsigned j) const { return get_column_value(j).is_int(); } bool model_is_int_feasible() const; @@ -393,6 +395,7 @@ class lar_solver : public column_namer { inline column_index to_column_index(unsigned v) const { return column_index(external_to_column_index(v)); } bool external_is_used(unsigned) const; void pop(unsigned k); + unsigned num_scopes() const { return m_term_count.stack_size(); } bool compare_values(var_index j, lconstraint_kind kind, const mpq& right_side); var_index add_term(const vector>& coeffs, unsigned ext_i); void register_existing_terms(); @@ -503,7 +506,7 @@ class lar_solver : public column_namer { if (tv::is_term(j)) { j = m_var_register.external_to_local(j); } - return m_columns_to_ul_pairs()[j].upper_bound_witness(); + return m_columns_to_ul_pairs[j].upper_bound_witness(); } inline const impq& get_upper_bound(column_index j) const { @@ -591,7 +594,7 @@ class lar_solver : public column_namer { if (tv::is_term(j)) { j = m_var_register.external_to_local(j); } - return m_columns_to_ul_pairs()[j].lower_bound_witness(); + return m_columns_to_ul_pairs[j].lower_bound_witness(); } inline tv column2tv(column_index const& c) const { diff --git a/src/math/lp/monomial_bounds.cpp b/src/math/lp/monomial_bounds.cpp index 0d6fb47e7..08ae6f249 100644 --- a/src/math/lp/monomial_bounds.cpp +++ b/src/math/lp/monomial_bounds.cpp @@ -12,8 +12,6 @@ #include "math/lp/nla_intervals.h" #include "math/lp/numeric_pair.h" -#define UNIT_PROPAGATE_BOUNDS 0 - namespace nla { monomial_bounds::monomial_bounds(core* c): @@ -21,16 +19,17 @@ namespace nla { dep(c->m_intervals.get_dep_intervals()) {} void monomial_bounds::propagate() { - for (lpvar v : c().m_to_refine) + for (lpvar v : c().m_to_refine) { propagate(c().emon(v)); + if (add_lemma()) + break; + } } - bool monomial_bounds::is_too_big(mpq const& q) const { return rational(q).bitsize() > 256; } - /** * Accumulate product of variables in monomial starting at position 'start' */ @@ -53,8 +52,9 @@ namespace nla { * a bounds axiom. */ bool monomial_bounds::propagate_value(dep_interval& range, lpvar v) { - auto val = c().val(v); - if (dep.is_below(range, val)) { + // auto val = c().val(v); + bool propagated = false; + if (should_propagate_upper(range, v, 1)) { auto const& upper = dep.upper(range); auto cmp = dep.upper_is_open(range) ? llc::LT : llc::LE; ++c().lra.settings().stats().m_nla_propagate_bounds; @@ -72,9 +72,9 @@ namespace nla { lemma |= ineq(v, cmp, upper); TRACE("nla_solver", dep.display(tout << c().val(v) << " > ", range) << "\n" << lemma << "\n";); } - return true; + propagated = true; } - else if (dep.is_above(range, val)) { + if (should_propagate_lower(range, v, 1)) { auto const& lower = dep.lower(range); auto cmp = dep.lower_is_open(range) ? llc::GT : llc::GE; ++c().lra.settings().stats().m_nla_propagate_bounds; @@ -92,11 +92,45 @@ namespace nla { lemma |= ineq(v, cmp, lower); TRACE("nla_solver", dep.display(tout << c().val(v) << " < ", range) << "\n" << lemma << "\n";); } - return true; + propagated = true; } - else { + return propagated; + } + + bool monomial_bounds::should_propagate_lower(dep_interval const& range, lpvar v, unsigned p) { + if (dep.lower_is_inf(range)) return false; - } + u_dependency* d = nullptr; + rational bound; + bool is_strict; + if (!c().has_lower_bound(v, d, bound, is_strict)) + return true; + auto const& lower = dep.lower(range); + if (p > 1) + bound = power(bound, p); + if (bound < lower) + return true; + if (bound > lower) + return false; + return !is_strict && dep.lower_is_open(range); + } + + bool monomial_bounds::should_propagate_upper(dep_interval const& range, lpvar v, unsigned p) { + if (dep.upper_is_inf(range)) + return false; + u_dependency* d = nullptr; + rational bound; + bool is_strict; + if (!c().has_upper_bound(v, d, bound, is_strict)) + return true; + auto const& upper = dep.upper(range); + if (p > 1) + bound = power(bound, p); + if (bound > upper) + return true; + if (bound < upper) + return false; + return !is_strict && dep.upper_is_open(range); } /** @@ -135,68 +169,86 @@ namespace nla { SASSERT(p > 0); if (p == 1) return propagate_value(range, v); - auto val_v = c().val(v); - auto val = power(val_v, p); rational r; - if (dep.is_below(range, val)) { + if (should_propagate_upper(range, v, p)) { // v.upper^p > range.upper + SASSERT(c().has_upper_bound(v)); lp::explanation ex; dep.get_upper_dep(range, ex); + // p even, range.upper < 0, v^p >= 0 -> infeasible if (p % 2 == 0 && rational(dep.upper(range)).is_neg()) { ++c().lra.settings().stats().m_nla_propagate_bounds; new_lemma lemma(c(), "range requires a non-negative upper bound"); lemma &= ex; return true; } - else if (rational(dep.upper(range)).root(p, r)) { - // v = -2, [-4,-3]^3 < v^3 -> add bound v <= -3 - // v = -2, [-1,+1]^2 < v^2 -> add bound v >= -1 - if ((p % 2 == 1) || val_v.is_pos()) { - ++c().lra.settings().stats().m_nla_propagate_bounds; - auto le = dep.upper_is_open(range) ? llc::LT : llc::LE; - if (c().params().arith_nl_internal_bounds()) { - auto* d = dep.get_upper_dep(range); - propagate_bound(v, le, r, d); - } - else { - new_lemma lemma(c(), "propagate value - root case - upper bound of range is below value"); - lemma &= ex; - lemma |= ineq(v, le, r); - } - return true; - } - if (p % 2 == 0 && val_v.is_neg()) { - ++c().lra.settings().stats().m_nla_propagate_bounds; - SASSERT(!r.is_neg()); - auto ge = dep.upper_is_open(range) ? llc::GT : llc::GE; - if (c().params().arith_nl_internal_bounds()) { - auto* d = dep.get_upper_dep(range); - propagate_bound(v, ge, -r, d); - } - else { - new_lemma lemma(c(), "propagate value - root case - upper bound of range is below negative value"); - lemma &= ex; - lemma |= ineq(v, ge, -r); - } - return true; - } - } - // TBD: add bounds as long as difference to val is above some epsilon. - } - else if (dep.is_above(range, val)) { - if (rational(dep.lower(range)).root(p, r)) { + // v.upper < 0, but v^p > range.upper -> infeasible. + if (p % 2 == 0 && c().get_upper_bound(v) < 0) { ++c().lra.settings().stats().m_nla_propagate_bounds; - lp::explanation ex; - dep.get_lower_dep(range, ex); - auto ge = dep.lower_is_open(range) ? llc::GT : llc::GE; - auto le = dep.lower_is_open(range) ? llc::LT : llc::LE; - new_lemma lemma(c(), "propagate value - root case - lower bound of range is above value"); + new_lemma lemma(c(), "range requires a non-negative upper bound"); lemma &= ex; - lemma |= ineq(v, ge, r); - if (p % 2 == 0) - lemma |= ineq(v, le, -r); + lemma.explain_existing_upper_bound(v); return true; } - // TBD: add bounds as long as difference to val is above some epsilon. + SASSERT(p % 2 == 1 || c().get_upper_bound(v) >= 0); + + if (rational(dep.upper(range)).root(p, r)) { + // v = -2, [-4,-3]^3 < v^3 -> add bound v <= -3 + // v = -2, [-1,+1]^2 < v^2 -> add bound v >= -1 + ++c().lra.settings().stats().m_nla_propagate_bounds; + auto le = dep.upper_is_open(range) ? llc::LT : llc::LE; + if (c().params().arith_nl_internal_bounds()) { + auto* d = dep.get_upper_dep(range); + if (p % 2 == 0) + d = dep.mk_join(d, c().lra.get_column_upper_bound_witness(v)); + propagate_bound(v, le, r, d); + } + else { + new_lemma lemma(c(), "propagate value - root case - upper bound of range is below value"); + lemma &= ex; + if (p % 2 == 0) + lemma.explain_existing_upper_bound(v); + lemma |= ineq(v, le, r); + } + return true; + } + } + + if (should_propagate_lower(range, v, p)) { // v.lower^p < range.lower + // + // range.lower < 0 -> v.lower >= root(p, range.lower) + // range.lower >= 0, p odd -> v.lower >= root(p, range.lower) + // range.lower >= 0, p even, v.lower >= 0 -> v.lower >= root(p, range.lower) + // default: + // v.lower >= root(p, range.lower) || (p even & v.upper <= -root(p, range.lower)) + // + // pre-condition: p even -> range.lower >= 0 + // + if (rational(dep.lower(range)).root(p, r)) { + ++c().lra.settings().stats().m_nla_propagate_bounds; + auto ge = dep.lower_is_open(range) ? llc::GT : llc::GE; + auto le = dep.lower_is_open(range) ? llc::LT : llc::LE; + if (c().params().arith_nl_internal_bounds()) { + if (rational(dep.lower(range)).is_neg() || p % 2 == 1) { + auto* d = dep.get_lower_dep(range); + propagate_bound(v, ge, r, d); + return true; + } + if (c().get_lower_bound(v) >= 0) { + auto* d = dep.get_lower_dep(range); + d = dep.mk_join(d, c().lra.get_column_lower_bound_witness(v)); + propagate_bound(v, ge, r, d); + return true; + } + } + lp::explanation ex; + dep.get_lower_dep(range, ex); + new_lemma lemma(c(), "propagate value - root case - lower bound of range is above value"); + lemma &= ex; + lemma |= ineq(v, ge, r); + if (p % 2 == 0) + lemma |= ineq(v, le, -r); + return true; + } } return false; } @@ -316,28 +368,30 @@ namespace nla { continue; monic& m = c().emon(v); unit_propagate(m); - if (c().lra.get_status() == lp::lp_status::INFEASIBLE) { - lp::explanation exp; - c().lra.get_infeasibility_explanation(exp); - new_lemma lemma(c(), "propagate fixed - infeasible lra"); - lemma &= exp; + if (add_lemma()) break; - } if (c().m_conflicts > 0) break; } } + bool monomial_bounds::add_lemma() { + if (c().lra.get_status() != lp::lp_status::INFEASIBLE) + return false; + lp::explanation exp; + c().lra.get_infeasibility_explanation(exp); + new_lemma lemma(c(), "propagate fixed - infeasible lra"); + lemma &= exp; + return true; + } + void monomial_bounds::unit_propagate(monic & m) { if (m.is_propagated()) return; lpvar w, fixed_to_zero; - if (!is_linear(m, w, fixed_to_zero)) { - if (c().params().arith_nl_internal_bounds()) - propagate(m); + if (!is_linear(m, w, fixed_to_zero)) return; - } c().emons().set_propagated(m); diff --git a/src/math/lp/monomial_bounds.h b/src/math/lp/monomial_bounds.h index 554fb6109..19043e072 100644 --- a/src/math/lp/monomial_bounds.h +++ b/src/math/lp/monomial_bounds.h @@ -17,6 +17,8 @@ namespace nla { class monomial_bounds : common { dep_intervals& dep; + bool should_propagate_lower(dep_interval const& range, lpvar v, unsigned p); + bool should_propagate_upper(dep_interval const& range, lpvar v, unsigned p); void propagate_bound(lpvar v, lp::lconstraint_kind cmp, rational const& q, u_dependency* d); void var2interval(lpvar v, scoped_dep_interval& i); bool is_too_big(mpq const& q) const; @@ -34,6 +36,7 @@ namespace nla { void analyze_monomial(monic const& m, unsigned& num_free, lpvar& free_v, unsigned& power) const; bool is_free(lpvar v) const; bool is_zero(lpvar v) const; + bool add_lemma(); // monomial propagation void unit_propagate(monic & m); diff --git a/src/math/lp/ul_pair.h b/src/math/lp/ul_pair.h index 37fd6e9ed..354070281 100644 --- a/src/math/lp/ul_pair.h +++ b/src/math/lp/ul_pair.h @@ -44,7 +44,7 @@ inline bool compare(const std::pair & a, const std::paircheck(); - switch (r) { case l_false: add_lemmas(); @@ -3293,9 +3292,8 @@ public: for (auto ev : m_explanation) set_evidence(ev.ci(), m_core, m_eqs); - SASSERT(!m_core.empty() || !m_eqs.empty()); - // SASSERT(validate_conflict(m_core, m_eqs)); + // VERIFY(validate_conflict(m_core, m_eqs)); if (is_conflict) { ctx().set_conflict( ctx().mk_justification( @@ -3509,6 +3507,8 @@ public: bool validate_conflict(literal_vector const& core, svector const& eqs) { if (params().m_arith_mode != arith_solver_id::AS_NEW_ARITH) return true; + + VERIFY(!m_core.empty() || !m_eqs.empty()); scoped_arith_mode _sa(ctx().get_fparams()); context nctx(m, ctx().get_fparams(), ctx().get_params()); add_background(nctx);