From 5e3df9ee77c51ffba07e727f8ee680b9382f66e9 Mon Sep 17 00:00:00 2001 From: Nikolaj Bjorner Date: Sat, 19 Aug 2023 17:44:09 -0700 Subject: [PATCH] Arith min max (#6864) * prepare for dependencies Signed-off-by: Nikolaj Bjorner * snapshot Signed-off-by: Nikolaj Bjorner * more refactoring Signed-off-by: Nikolaj Bjorner * more refactoring Signed-off-by: Nikolaj Bjorner * build Signed-off-by: Nikolaj Bjorner * pass in u_dependency_manager Signed-off-by: Nikolaj Bjorner * address NYIs Signed-off-by: Nikolaj Bjorner * more refactoring names Signed-off-by: Nikolaj Bjorner * eq_explanation update Signed-off-by: Nikolaj Bjorner * add outline of bounds improvement functionality Signed-off-by: Nikolaj Bjorner * fix unit tests Signed-off-by: Nikolaj Bjorner * remove unused structs Signed-off-by: Nikolaj Bjorner * more bounds Signed-off-by: Nikolaj Bjorner * more bounds Signed-off-by: Nikolaj Bjorner * convert more internals to use u_dependency instead of constraint_index Signed-off-by: Nikolaj Bjorner * convert more internals to use u_dependency instead of constraint_index Signed-off-by: Nikolaj Bjorner * remember to push/pop scopes Signed-off-by: Nikolaj Bjorner * use the main function for updating bounds Signed-off-by: Nikolaj Bjorner * na Signed-off-by: Nikolaj Bjorner * na Signed-off-by: Nikolaj Bjorner * remove reset of shared dep manager Signed-off-by: Nikolaj Bjorner * disable improve-bounds, add statistics Signed-off-by: Nikolaj Bjorner --------- Signed-off-by: Nikolaj Bjorner --- src/ast/simplifiers/bound_simplifier.h | 3 +- src/math/grobner/pdd_solver.cpp | 5 +- src/math/grobner/pdd_solver.h | 5 +- src/math/interval/dep_intervals.h | 17 ++- src/math/lp/emonics.cpp | 7 +- src/math/lp/gomory.cpp | 20 +-- src/math/lp/hnf_cutter.cpp | 27 ++-- src/math/lp/hnf_cutter.h | 6 +- src/math/lp/horner.cpp | 2 +- src/math/lp/int_gcd_test.cpp | 7 +- src/math/lp/int_solver.cpp | 4 +- src/math/lp/int_solver.h | 4 +- src/math/lp/lar_constraints.h | 42 +++--- src/math/lp/lar_core_solver.h | 23 ++-- src/math/lp/lar_core_solver_def.h | 15 +-- src/math/lp/lar_solver.cpp | 163 ++++++++++++++-------- src/math/lp/lar_solver.h | 63 ++++++--- src/math/lp/lp_bound_propagator.h | 17 ++- src/math/lp/lp_primal_core_solver.h | 171 +++++++++++------------ src/math/lp/lp_types.h | 3 + src/math/lp/monomial_bounds.cpp | 10 +- src/math/lp/nla_common.cpp | 12 +- src/math/lp/nla_common.h | 15 --- src/math/lp/nla_core.cpp | 180 ++++++++++++++----------- src/math/lp/nla_core.h | 24 ++-- src/math/lp/nla_defs.h | 9 +- src/math/lp/nla_grobner.cpp | 41 +++--- src/math/lp/nla_grobner.h | 2 +- src/math/lp/nla_intervals.cpp | 38 +++--- src/math/lp/nla_intervals.h | 7 +- src/math/lp/nra_solver.cpp | 4 +- src/math/lp/ul_pair.h | 37 +++-- src/math/lp/var_eqs.h | 33 +++-- src/sat/sat_anf_simplifier.cpp | 3 +- src/sat/smt/arith_axioms.cpp | 6 +- src/sat/smt/arith_solver.cpp | 46 ++++--- src/sat/smt/arith_solver.h | 8 +- src/smt/theory_lra.cpp | 58 ++++---- src/test/pdd_solver.cpp | 7 +- src/util/dependency.h | 15 +-- 40 files changed, 630 insertions(+), 529 deletions(-) diff --git a/src/ast/simplifiers/bound_simplifier.h b/src/ast/simplifiers/bound_simplifier.h index 0e3fff239..9bd4b1908 100644 --- a/src/ast/simplifiers/bound_simplifier.h +++ b/src/ast/simplifiers/bound_simplifier.h @@ -37,6 +37,7 @@ class bound_simplifier : public dependent_expr_simplifier { unsynch_mpq_manager nm; small_object_allocator m_alloc; bound_propagator bp; + u_dependency_manager m_dep_manager; dep_intervals m_interval; ptr_vector m_var2expr; unsigned_vector m_expr2var; @@ -105,7 +106,7 @@ public: a(m), m_rewriter(m), bp(nm, m_alloc, p), - m_interval(m.limit()), + m_interval(m_dep_manager, m.limit()), m_trail(m), m_num_buffer(nm) { updt_params(p); diff --git a/src/math/grobner/pdd_solver.cpp b/src/math/grobner/pdd_solver.cpp index 63c5ad835..e689c78a9 100644 --- a/src/math/grobner/pdd_solver.cpp +++ b/src/math/grobner/pdd_solver.cpp @@ -59,9 +59,10 @@ namespace dd { */ - solver::solver(reslimit& lim, pdd_manager& m) : + solver::solver(reslimit& lim, u_dependency_manager& dm, pdd_manager& m) : m(m), - m_limit(lim) + m_limit(lim), + m_dep_manager(dm) {} solver::~solver() { diff --git a/src/math/grobner/pdd_solver.h b/src/math/grobner/pdd_solver.h index 40f8fdce2..bc20c21b4 100644 --- a/src/math/grobner/pdd_solver.h +++ b/src/math/grobner/pdd_solver.h @@ -112,6 +112,7 @@ private: pdd_manager& m; reslimit& m_limit; + u_dependency_manager& m_dep_manager; stats m_stats; config m_config; print_dep_t m_print_dep; @@ -119,12 +120,11 @@ private: equation_vector m_processed; equation_vector m_to_simplify; vector> m_subst; - mutable u_dependency_manager m_dep_manager; equation_vector m_all_eqs; equation* m_conflict = nullptr; bool m_too_complex; public: - solver(reslimit& lim, pdd_manager& m); + solver(reslimit& lim, u_dependency_manager& dm, pdd_manager& m); ~solver(); pdd_manager& get_manager() { return m; } @@ -144,7 +144,6 @@ public: void saturate(); equation_vector const& equations(); - u_dependency_manager& dep() const { return m_dep_manager; } void collect_statistics(statistics & st) const; std::ostream& display(std::ostream& out, const equation& eq) const; diff --git a/src/math/interval/dep_intervals.h b/src/math/interval/dep_intervals.h index d641a294d..f9768b6f4 100644 --- a/src/math/interval/dep_intervals.h +++ b/src/math/interval/dep_intervals.h @@ -27,6 +27,7 @@ #include "math/interval/interval.h" class dep_intervals { + public: enum with_deps_t { with_deps, without_deps }; @@ -142,8 +143,9 @@ private: public: typedef interval_manager::interval interval; + u_dependency_manager& m_dep_manager; mutable unsynch_mpq_manager m_num_manager; - mutable u_dependency_manager m_dep_manager; + im_config m_config; mutable interval_manager m_imanager; @@ -158,9 +160,10 @@ public: public: u_dependency_manager& dep_manager() { return m_dep_manager; } - dep_intervals(reslimit& lim) : - m_config(m_num_manager, m_dep_manager), - m_imanager(lim, im_config(m_num_manager, m_dep_manager)) + dep_intervals(u_dependency_manager& dm, reslimit& lim) : + m_dep_manager(dm), + m_config(m_num_manager, dm), + m_imanager(lim, im_config(m_num_manager, dm)) {} std::ostream& display(std::ostream& out, const interval& i) const; @@ -335,15 +338,17 @@ public: bool is_empty(interval const& a) const; void set_interval_for_scalar(interval&, const rational&); + template void linearize(u_dependency* dep, T& expl) const { vector v; m_dep_manager.linearize(dep, v); - for (unsigned ci: v) + for (auto ci: v) expl.push_back(ci); } - void reset() { m_dep_manager.reset(); } + + void reset() { } void del(interval& i) { m_imanager.del(i); } diff --git a/src/math/lp/emonics.cpp b/src/math/lp/emonics.cpp index bcdb81dd8..a8ac63689 100644 --- a/src/math/lp/emonics.cpp +++ b/src/math/lp/emonics.cpp @@ -517,11 +517,10 @@ bool emonics::invariant() const { TRACE("nla_solver_mons", display(tout);); // the variable index contains exactly the active monomials unsigned mons = 0; - for (lpvar v = 0; v < m_var2index.size(); v++) { - if (is_monic_var(v)) { + for (lpvar v = 0; v < m_var2index.size(); v++) + if (is_monic_var(v)) mons++; - } - } + if (m_monics.size() != mons) { TRACE("nla_solver_mons", tout << "missmatch of monic vars\n";); return false; diff --git a/src/math/lp/gomory.cpp b/src/math/lp/gomory.cpp index 2ecbc49ac..775025018 100644 --- a/src/math/lp/gomory.cpp +++ b/src/math/lp/gomory.cpp @@ -50,9 +50,13 @@ class create_cut { bool at_upper(unsigned j) const { return lia.at_upper(j); } const impq & lower_bound(unsigned j) const { return lia.lower_bound(j); } const impq & upper_bound(unsigned j) const { return lia.upper_bound(j); } - constraint_index column_lower_bound_constraint(unsigned j) const { return lia.column_lower_bound_constraint(j); } - constraint_index column_upper_bound_constraint(unsigned j) const { return lia.column_upper_bound_constraint(j); } + u_dependency* column_lower_bound_constraint(unsigned j) const { return lia.column_lower_bound_constraint(j); } + u_dependency* column_upper_bound_constraint(unsigned j) const { return lia.column_upper_bound_constraint(j); } bool column_is_fixed(unsigned j) const { return lia.lra.column_is_fixed(j); } + void push_explanation(u_dependency* d) { + for (auto ci : lia.lra.flatten(d)) + m_ex->push_back(ci); + } void int_case_in_gomory_cut(unsigned j) { lp_assert(is_int(j) && m_fj.is_pos()); @@ -67,7 +71,7 @@ class create_cut { new_a = m_fj <= m_one_minus_f ? m_fj / m_one_minus_f : ((1 - m_fj) / m_f); lp_assert(new_a.is_pos()); m_k.addmul(new_a, lower_bound(j).x); - m_ex->push_back(column_lower_bound_constraint(j)); + push_explanation(column_lower_bound_constraint(j)); } else { lp_assert(at_upper(j)); @@ -75,7 +79,7 @@ class create_cut { new_a = - (m_fj <= m_f ? m_fj / m_f : ((1 - m_fj) / m_one_minus_f)); lp_assert(new_a.is_neg()); m_k.addmul(new_a, upper_bound(j).x); - m_ex->push_back(column_upper_bound_constraint(j)); + push_explanation(column_upper_bound_constraint(j)); } m_t.add_monomial(new_a, j); m_lcm_den = lcm(m_lcm_den, denominator(new_a)); @@ -99,7 +103,7 @@ class create_cut { new_a = - a / m_f; m_k.addmul(new_a, lower_bound(j).x); // is it a faster operation than // k += lower_bound(j).x * new_a; - m_ex->push_back(column_lower_bound_constraint(j)); + push_explanation(column_lower_bound_constraint(j)); } else { lp_assert(at_upper(j)); @@ -110,7 +114,7 @@ class create_cut { // the delta is positive works again m_one_minus_f new_a = a / m_one_minus_f; m_k.addmul(new_a, upper_bound(j).x); // k += upper_bound(j).x * new_a; - m_ex->push_back(column_upper_bound_constraint(j)); + push_explanation(column_upper_bound_constraint(j)); } m_t.add_monomial(new_a, j); TRACE("gomory_cut_detail_real", tout << "add " << new_a << "*v" << j << ", k: " << m_k << "\n"; @@ -313,8 +317,8 @@ public: // use -p.coeff() to make the format compatible with the format used in: Integrating Simplex with DPLL(T) try { if (lia.is_fixed(j)) { - m_ex->push_back(column_lower_bound_constraint(j)); - m_ex->push_back(column_upper_bound_constraint(j)); + push_explanation(column_lower_bound_constraint(j)); + push_explanation(column_upper_bound_constraint(j)); continue; } if (is_real(j)) diff --git a/src/math/lp/hnf_cutter.cpp b/src/math/lp/hnf_cutter.cpp index 3c4ea10ab..d688eeb63 100644 --- a/src/math/lp/hnf_cutter.cpp +++ b/src/math/lp/hnf_cutter.cpp @@ -39,7 +39,7 @@ namespace lp { m_overflow = false; } - void hnf_cutter::add_term(const lar_term* t, const mpq &rs, constraint_index ci, bool upper_bound) { + void hnf_cutter::add_term(const lar_term* t, const mpq &rs, u_dependency* ci, bool upper_bound) { m_terms.push_back(t); m_terms_upper.push_back(upper_bound); if (upper_bound) @@ -146,7 +146,7 @@ namespace lp { } #endif void hnf_cutter::shrink_explanation(const svector& basis_rows) { - svector new_expl; + ptr_vector new_expl; for (unsigned i : basis_rows) { new_expl.push_back(m_constraints_for_explanation[i]); } @@ -232,10 +232,10 @@ branch y_i >= ceil(y0_i) is impossible. void hnf_cutter::try_add_term_to_A_for_hnf(tv const &i) { mpq rs; const lar_term& t = lra.get_term(i); - constraint_index ci; + u_dependency* dep; bool upper_bound; - if (!is_full() && lra.get_equality_and_right_side_for_term_on_current_x(i, rs, ci, upper_bound)) { - add_term(&t, rs, ci, upper_bound); + if (!is_full() && lra.get_equality_and_right_side_for_term_on_current_x(i, rs, dep, upper_bound)) { + add_term(&t, rs, dep, upper_bound); } } @@ -259,8 +259,9 @@ branch y_i >= ceil(y0_i) is impossible. } lia.settings().stats().m_hnf_cutter_calls++; TRACE("hnf_cut", tout << "settings().stats().m_hnf_cutter_calls = " << lia.settings().stats().m_hnf_cutter_calls << "\n"; - for (unsigned i : constraints_for_explanation()) { - lra.constraints().display(tout, i); + for (u_dependency* d : constraints_for_explanation()) { + for (auto ci : lra.flatten(d)) + lra.constraints().display(tout, ci); } tout << lra.constraints(); ); @@ -277,16 +278,16 @@ branch y_i >= ceil(y0_i) is impossible. TRACE("hnf_cut", lra.print_term(lia.m_t, tout << "cut:"); tout << " <= " << lia.m_k << std::endl; - for (unsigned i : constraints_for_explanation()) { - lra.constraints().display(tout, i); - } + for (auto* dep : constraints_for_explanation()) + for (auto ci : lra.flatten(dep)) + lra.constraints().display(tout, ci); ); lp_assert(lia.current_solution_is_inf_on_cut()); lia.settings().stats().m_hnf_cuts++; lia.m_ex->clear(); - for (unsigned i : constraints_for_explanation()) { - lia.m_ex->push_back(i); - } + for (u_dependency* dep : constraints_for_explanation()) + for (auto ci : lia.lra.flatten(dep)) + lia.m_ex->push_back(ci); } return r; } diff --git a/src/math/lp/hnf_cutter.h b/src/math/lp/hnf_cutter.h index b3530ea29..74fb52327 100644 --- a/src/math/lp/hnf_cutter.h +++ b/src/math/lp/hnf_cutter.h @@ -34,7 +34,7 @@ class hnf_cutter { general_matrix m_A; vector m_terms; vector m_terms_upper; - svector m_constraints_for_explanation; + ptr_vector m_constraints_for_explanation; vector m_right_sides; mpq m_abs_max; bool m_overflow; @@ -55,13 +55,13 @@ private: unsigned terms_count() const { return m_terms.size(); } const mpq & abs_max() const { return m_abs_max; } const vector& terms() const { return m_terms; } - const svector& constraints_for_explanation() const { return m_constraints_for_explanation; } + const ptr_vector& constraints_for_explanation() const { return m_constraints_for_explanation; } const vector & right_sides() const { return m_right_sides; } bool is_full() const; void clear(); - void add_term(const lar_term* t, const mpq &rs, constraint_index ci, bool upper_bound); + void add_term(const lar_term* t, const mpq &rs, u_dependency* ci, bool upper_bound); void print(std::ostream & out); diff --git a/src/math/lp/horner.cpp b/src/math/lp/horner.cpp index 40f16d709..0cd62ecaf 100644 --- a/src/math/lp/horner.cpp +++ b/src/math/lp/horner.cpp @@ -103,7 +103,7 @@ bool horner::horner_lemmas() { return false; } c().lp_settings().stats().m_horner_calls++; - const auto& matrix = c().m_lar_solver.A_r(); + const auto& matrix = c().lra.A_r(); // choose only rows that depend on m_to_refine variables std::set rows_to_check; for (lpvar j : c().m_to_refine) { diff --git a/src/math/lp/int_gcd_test.cpp b/src/math/lp/int_gcd_test.cpp index 4801cc436..f3b1b6389 100644 --- a/src/math/lp/int_gcd_test.cpp +++ b/src/math/lp/int_gcd_test.cpp @@ -250,10 +250,9 @@ namespace lp { } void int_gcd_test::add_to_explanation_from_fixed_or_boxed_column(unsigned j) { - constraint_index lc, uc; - lra.get_bound_constraint_witnesses_for_column(j, lc, uc); - lia.m_ex->push_back(lc); - lia.m_ex->push_back(uc); + auto* deps = lra.get_bound_constraint_witnesses_for_column(j); + for (auto d : lra.flatten(deps)) + lia.m_ex->push_back(d); } bool int_gcd_test::accumulate_parity(const row_strip & row, unsigned least_idx) { diff --git a/src/math/lp/int_solver.cpp b/src/math/lp/int_solver.cpp index dbb5dbae2..0ab7d1e40 100644 --- a/src/math/lp/int_solver.cpp +++ b/src/math/lp/int_solver.cpp @@ -242,11 +242,11 @@ namespace lp { return lra.has_inf_int(); } - constraint_index int_solver::column_upper_bound_constraint(unsigned j) const { + u_dependency* int_solver::column_upper_bound_constraint(unsigned j) const { return lra.get_column_upper_bound_witness(j); } - constraint_index int_solver::column_lower_bound_constraint(unsigned j) const { + u_dependency* int_solver::column_lower_bound_constraint(unsigned j) const { return lra.get_column_lower_bound_witness(j); } diff --git a/src/math/lp/int_solver.h b/src/math/lp/int_solver.h index c14fd0067..ab29ab7f4 100644 --- a/src/math/lp/int_solver.h +++ b/src/math/lp/int_solver.h @@ -113,8 +113,8 @@ private: public: std::ostream& display_column(std::ostream & out, unsigned j) const; - constraint_index column_upper_bound_constraint(unsigned j) const; - constraint_index column_lower_bound_constraint(unsigned j) const; + u_dependency* column_upper_bound_constraint(unsigned j) const; + u_dependency* column_lower_bound_constraint(unsigned j) const; bool current_solution_is_inf_on_cut() const; bool shift_var(unsigned j, unsigned range); diff --git a/src/math/lp/lar_constraints.h b/src/math/lp/lar_constraints.h index f8cffbe57..0a16353c1 100644 --- a/src/math/lp/lar_constraints.h +++ b/src/math/lp/lar_constraints.h @@ -1,21 +1,10 @@ /*++ Copyright (c) 2017 Microsoft Corporation -Module Name: - - - -Abstract: - - - Author: Lev Nachmanson (levnach) -Revision History: - - --*/ #pragma once @@ -53,15 +42,19 @@ class lar_base_constraint { mpq m_right_side; bool m_active; unsigned m_j; -public: + u_dependency* m_dep; + + public: virtual vector> coeffs() const = 0; - lar_base_constraint(unsigned j, lconstraint_kind kind, const mpq& right_side) :m_kind(kind), m_right_side(right_side), m_active(false), m_j(j) {} + lar_base_constraint(unsigned j, lconstraint_kind kind, u_dependency* dep, const mpq& right_side) : + m_kind(kind), m_right_side(right_side), m_active(false), m_j(j), m_dep(dep) {} virtual ~lar_base_constraint() = default; lconstraint_kind kind() const { return m_kind; } mpq const& rhs() const { return m_right_side; } unsigned column() const { return m_j; } + u_dependency* dep() const { return m_dep; } void activate() { m_active = true; } void deactivate() { m_active = false; } @@ -73,8 +66,8 @@ public: class lar_var_constraint: public lar_base_constraint { public: - lar_var_constraint(unsigned j, lconstraint_kind kind, const mpq& right_side) : - lar_base_constraint(j, kind, right_side) {} + lar_var_constraint(unsigned j, lconstraint_kind kind, u_dependency* dep, const mpq& right_side) : + lar_base_constraint(j, kind, dep, right_side) {} vector> coeffs() const override { vector> ret; @@ -88,8 +81,8 @@ public: class lar_term_constraint: public lar_base_constraint { const lar_term * m_term; public: - lar_term_constraint(unsigned j, const lar_term *t, lconstraint_kind kind, const mpq& right_side) : - lar_base_constraint(j, kind, right_side), m_term(t) {} + lar_term_constraint(unsigned j, const lar_term* t, lconstraint_kind kind, u_dependency* dep, const mpq& right_side) : + lar_base_constraint(j, kind, dep, right_side), m_term(t) {} vector> coeffs() const override { return m_term->coeffs_as_vector(); } unsigned size() const override { return m_term->size();} @@ -98,10 +91,11 @@ public: class constraint_set { region m_region; column_namer& m_namer; + u_dependency_manager& m_dep_manager; vector m_constraints; stacked_value m_constraint_count; unsigned_vector m_active; - stacked_value m_active_lim; + stacked_value m_active_lim; constraint_index add(lar_base_constraint* c) { constraint_index ci = m_constraints.size(); @@ -137,8 +131,13 @@ class constraint_set { return out << "constraint " << T_to_string(ci) << " is not found" << std::endl; } + u_dependency* mk_dep() { + return m_dep_manager.mk_leaf(m_constraints.size()); + } + public: - constraint_set(column_namer& cn): + constraint_set(u_dependency_manager& d, column_namer& cn): + m_dep_manager(d), m_namer(cn) {} ~constraint_set() { @@ -169,11 +168,12 @@ public: } constraint_index add_var_constraint(var_index j, lconstraint_kind k, mpq const& rhs) { - return add(new (m_region) lar_var_constraint(j, k, rhs)); + return add(new (m_region) lar_var_constraint(j, k, mk_dep(), rhs)); } constraint_index add_term_constraint(unsigned j, const lar_term* t, lconstraint_kind k, mpq const& rhs) { - return add(new (m_region) lar_term_constraint(j, t, k, rhs)); + auto* dep = mk_dep(); + return add(new (m_region) lar_term_constraint(j, t, k, dep, rhs)); } // future behavior uses activation bit. diff --git a/src/math/lp/lar_core_solver.h b/src/math/lp/lar_core_solver.h index 0780d21af..4a5cb20b1 100644 --- a/src/math/lp/lar_core_solver.h +++ b/src/math/lp/lar_core_solver.h @@ -21,9 +21,10 @@ class lar_core_solver { int m_infeasible_sum_sign; // todo: get rid of this field vector> m_right_sides_dummy; vector m_costs_dummy; - -public: stacked_value m_stacked_simplex_strategy; + +public: + stacked_vector m_column_types; // r - solver fields, for rational numbers vector> m_r_x; // the solution @@ -43,8 +44,6 @@ public: const column_namer & column_names ); - lp_settings & settings() { return m_r_solver.m_settings;} - const lp_settings & settings() const { return m_r_solver.m_settings;} int get_infeasible_sum_sign() const { return m_infeasible_sum_sign; } @@ -57,8 +56,7 @@ public: void fill_not_improvable_zero_sum_from_inf_row(); column_type get_column_type(unsigned j) { return m_column_types[j];} - - + void print_pivot_row(std::ostream & out, unsigned row_index) const { for (unsigned j : m_r_solver.m_pivot_row.m_index) { if (numeric_traits::is_pos(m_r_solver.m_pivot_row.m_data[j])) @@ -68,9 +66,9 @@ public: out << " +" << m_r_solver.column_name(m_r_solver.m_basis[row_index]) << std::endl; - for (unsigned j : m_r_solver.m_pivot_row.m_index) { + for (unsigned j : m_r_solver.m_pivot_row.m_index) m_r_solver.print_column_bound_info(j, out); - } + m_r_solver.print_column_bound_info(m_r_solver.m_basis[row_index], out); } @@ -110,10 +108,7 @@ public: m_column_types.push(); // rational m_r_lower_bounds.push(); - m_r_upper_bounds.push(); - - - + m_r_upper_bounds.push(); } void pop(unsigned k) { @@ -127,11 +122,9 @@ public: m_r_solver.m_d.resize(m_r_A.column_count()); m_stacked_simplex_strategy.pop(k); - settings().set_simplex_strategy(m_stacked_simplex_strategy); + m_r_solver.m_settings.set_simplex_strategy(m_stacked_simplex_strategy); lp_assert(m_r_solver.basis_heading_is_correct()); } - - bool r_basis_is_OK() const { #ifdef Z3DEBUG diff --git a/src/math/lp/lar_core_solver_def.h b/src/math/lp/lar_core_solver_def.h index 942c87b3b..b6743873f 100644 --- a/src/math/lp/lar_core_solver_def.h +++ b/src/math/lp/lar_core_solver_def.h @@ -88,19 +88,18 @@ void lar_core_solver::solve() { lp_assert(m_r_solver.inf_heap_is_correct()); TRACE("find_feas_stats", tout << "infeasibles = " << m_r_solver.inf_heap_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); - TRACE("lar_solver", tout << m_r_solver.get_status() << "\n";); - return; + m_r_solver.set_status(lp_status::OPTIMAL); + TRACE("lar_solver", tout << m_r_solver.get_status() << "\n";); + return; } - ++settings().stats().m_need_to_solve_inf; + ++m_r_solver.m_settings.stats().m_need_to_solve_inf; lp_assert( r_basis_is_OK()); - - + if (m_r_solver.m_look_for_feasible_solution_only) //todo : should it be set? m_r_solver.find_feasible_solution(); - else { + else m_r_solver.solve(); - } + lp_assert(r_basis_is_OK()); switch (m_r_solver.get_status()) diff --git a/src/math/lp/lar_solver.cpp b/src/math/lp/lar_solver.cpp index ddfdb1437..6159c436e 100644 --- a/src/math/lp/lar_solver.cpp +++ b/src/math/lp/lar_solver.cpp @@ -26,7 +26,7 @@ namespace lp { m_mpq_lar_core_solver(m_settings, *this), m_var_register(false), m_term_register(true), - m_constraints(*this) {} + m_constraints(m_dependencies, *this) {} // start or ends tracking the rows that were changed by solve() void lar_solver::track_touched_rows(bool v) { @@ -218,8 +218,14 @@ namespace lp { // this is the case when the lower bound is in conflict with the upper one const ul_pair& ul = m_columns_to_ul_pairs[m_crossed_bounds_column]; - evidence.add_pair(ul.upper_bound_witness(), numeric_traits::one()); - evidence.add_pair(ul.lower_bound_witness(), -numeric_traits::one()); + svector deps; + m_dependencies.linearize(ul.upper_bound_witness(), deps); + for (auto d : deps) + evidence.add_pair(d, numeric_traits::one()); + deps.reset(); + m_dependencies.linearize(ul.lower_bound_witness(), deps); + for (auto d : deps) + evidence.add_pair(d, -numeric_traits::one()); } void lar_solver::push() { @@ -232,6 +238,7 @@ namespace lp { m_term_count.push(); m_constraints.push(); m_usage_in_terms.push(); + m_dependencies.push_scope(); } void lar_solver::clean_popped_elements(unsigned n, indexed_uint_set& set) { @@ -297,6 +304,7 @@ namespace lp { lp_assert(sizes_are_correct()); lp_assert(m_mpq_lar_core_solver.m_r_solver.reduced_costs_are_correct_tableau()); m_usage_in_terms.pop(k); + m_dependencies.pop_scope(k); set_status(lp_status::UNKNOWN); } @@ -320,6 +328,39 @@ namespace lp { } } + bool lar_solver::improve_bound(lpvar j, bool improve_lower_bound) { + lar_term term = get_term_to_maximize(j); + if (improve_lower_bound) + term.negate(); + impq bound; + if (!maximize_term_on_tableau(term, bound)) + return false; + + return false; + // TODO + if (improve_lower_bound) { + if (column_has_lower_bound(j) && bound.x == column_lower_bound(j).x) + return false; + SASSERT(!column_has_lower_bound(j) || column_lower_bound(j).x < bound.x); + + // TODO - explain new lower bound. + // Seems the relevant information is in the "costs" that are used when + // setting the optimization objecive. The costs are cleared after a call so + // maybe have some way of extracting bound dependencies from the costs. + u_dependency* dep = nullptr; + update_column_type_and_bound(j, bound.y > 0 ? lconstraint_kind::GT : lconstraint_kind::GE, bound.x, dep); + } + else { + if (column_has_upper_bound(j) && bound.x == column_upper_bound(j).x) + return false; + SASSERT(!column_has_upper_bound(j) || column_upper_bound(j).x > bound.x); + // similar for upper bounds + u_dependency* dep = nullptr; + update_column_type_and_bound(j, bound.y < 0 ? lconstraint_kind::LT : lconstraint_kind::LE, bound.x, dep); + } + return true; + } + bool lar_solver::costs_are_zeros_for_r_solver() const { for (unsigned j = 0; j < m_mpq_lar_core_solver.m_r_solver.m_costs.size(); j++) { lp_assert(is_zero(m_mpq_lar_core_solver.m_r_solver.m_costs[j])); @@ -349,7 +390,7 @@ namespace lp { jset.insert(rc.var()); } } - + for (unsigned j : jset) rslv.m_d[j] = zero_of_type(); @@ -577,15 +618,15 @@ namespace lp { } - void lar_solver::set_upper_bound_witness(var_index j, constraint_index ci) { + 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() = ci; + ul.upper_bound_witness() = dep; m_columns_to_ul_pairs[j] = ul; } - void lar_solver::set_lower_bound_witness(var_index j, constraint_index ci) { + 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() = ci; + ul.lower_bound_witness() = dep; m_columns_to_ul_pairs[j] = ul; } @@ -920,7 +961,7 @@ namespace lp { return ret; } - bool lar_solver::has_lower_bound(var_index var, constraint_index& ci, mpq& value, bool& is_strict) const { + bool lar_solver::has_lower_bound(var_index var, u_dependency*& ci, mpq& value, bool& is_strict) const { if (var >= m_columns_to_ul_pairs.size()) { // TBD: bounds on terms could also be used, caller may have to track these. @@ -928,7 +969,7 @@ namespace lp { } const ul_pair& ul = m_columns_to_ul_pairs[var]; ci = ul.lower_bound_witness(); - if (ci != null_ci) { + if (ci != nullptr) { auto& p = m_mpq_lar_core_solver.m_r_lower_bounds()[var]; value = p.x; is_strict = p.y.is_pos(); @@ -939,7 +980,7 @@ namespace lp { } } - bool lar_solver::has_upper_bound(var_index var, constraint_index& ci, mpq& value, bool& is_strict) const { + bool lar_solver::has_upper_bound(var_index var, u_dependency*& ci, mpq& value, bool& is_strict) const { if (var >= m_columns_to_ul_pairs.size()) { // TBD: bounds on terms could also be used, caller may have to track these. @@ -947,7 +988,7 @@ namespace lp { } const ul_pair& ul = m_columns_to_ul_pairs[var]; ci = ul.upper_bound_witness(); - if (ci != null_ci) { + if (ci != nullptr) { auto& p = m_mpq_lar_core_solver.m_r_upper_bounds()[var]; value = p.x; is_strict = p.y.is_neg(); @@ -1005,9 +1046,13 @@ namespace lp { int adj_sign = coeff.is_pos() ? inf_sign : -inf_sign; const ul_pair& ul = m_columns_to_ul_pairs[j]; - constraint_index bound_constr_i = adj_sign < 0 ? ul.upper_bound_witness() : ul.lower_bound_witness(); - lp_assert(m_constraints.valid_index(bound_constr_i)); - exp.add_pair(bound_constr_i, coeff); + u_dependency* bound_constr_i = adj_sign < 0 ? ul.upper_bound_witness() : ul.lower_bound_witness(); + svector deps; + m_dependencies.linearize(bound_constr_i, deps); + for (auto d : deps) { + lp_assert(m_constraints.valid_index(d)); + exp.add_pair(d, coeff); + } } } @@ -1698,12 +1743,12 @@ namespace lp { void lar_solver::activate_check_on_equal(constraint_index ci, unsigned& equal_column) { auto const& c = m_constraints[ci]; - update_column_type_and_bound_check_on_equal(c.column(), c.kind(), c.rhs(), ci, equal_column); + update_column_type_and_bound_check_on_equal(c.column(), c.rhs(), ci, equal_column); } void lar_solver::activate(constraint_index ci) { auto const& c = m_constraints[ci]; - update_column_type_and_bound(c.column(), c.kind(), c.rhs(), ci); + update_column_type_and_bound(c.column(), c.rhs(), ci); } mpq lar_solver::adjust_bound_for_int(lpvar j, lconstraint_kind& k, const mpq& bound) { @@ -1765,31 +1810,37 @@ namespace lp { } void lar_solver::update_column_type_and_bound(unsigned j, - lconstraint_kind kind, const mpq& right_side, constraint_index constr_index) { TRACE("lar_solver_feas", tout << "j = " << j << " was " << (this->column_is_feasible(j)?"feas":"non-feas") << std::endl;); m_constraints.activate(constr_index); - if (column_has_upper_bound(j)) - update_column_type_and_bound_with_ub(j, kind, right_side, constr_index); - else - update_column_type_and_bound_with_no_ub(j, kind, right_side, constr_index); - - if (is_base(j) && column_is_fixed(j)) - m_fixed_base_var_set.insert(j); - - TRACE("lar_solver_feas", tout << "j = " << j << " became " << (this->column_is_feasible(j)?"feas":"non-feas") << ", and " << (this->column_is_bounded(j)? "bounded":"non-bounded") << std::endl;); + lconstraint_kind kind = m_constraints[constr_index].kind(); + u_dependency* dep = m_constraints[constr_index].dep(); + update_column_type_and_bound(j, kind, right_side, dep); } + + void lar_solver::update_column_type_and_bound(unsigned j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep) { + if (column_has_upper_bound(j)) + update_column_type_and_bound_with_ub(j, kind, right_side, dep); + else + update_column_type_and_bound_with_no_ub(j, kind, right_side, dep); + + if (is_base(j) && column_is_fixed(j)) + m_fixed_base_var_set.insert(j); + + TRACE("lar_solver_feas", tout << "j = " << j << " became " << (this->column_is_feasible(j) ? "feas" : "non-feas") << ", and " << (this->column_is_bounded(j) ? "bounded" : "non-bounded") << std::endl;); + } + void lar_solver::insert_to_columns_with_changed_bounds(unsigned j) { m_columns_with_changed_bounds.insert(j); TRACE("lar_solver", tout << "column " << j << (column_is_feasible(j) ? " feas" : " non-feas") << "\n";); } + void lar_solver::update_column_type_and_bound_check_on_equal(unsigned j, - lconstraint_kind kind, const mpq& right_side, constraint_index constr_index, unsigned& equal_to_j) { - update_column_type_and_bound(j, kind, right_side, constr_index); + update_column_type_and_bound(j, right_side, constr_index); equal_to_j = null_lpvar; if (column_is_fixed(j)) { register_in_fixed_var_table(j, equal_to_j); @@ -1852,28 +1903,28 @@ namespace lp { } - void lar_solver::update_column_type_and_bound_with_ub(unsigned j, lp::lconstraint_kind kind, const mpq& right_side, unsigned constraint_index) { + void lar_solver::update_column_type_and_bound_with_ub(unsigned j, lp::lconstraint_kind kind, const mpq& right_side, u_dependency* dep) { SASSERT(column_has_upper_bound(j)); if (column_has_lower_bound(j)) { - update_bound_with_ub_lb(j, kind, right_side, constraint_index); + update_bound_with_ub_lb(j, kind, right_side, dep); } else { - update_bound_with_ub_no_lb(j, kind, right_side, constraint_index); + update_bound_with_ub_no_lb(j, kind, right_side, dep); } } - void lar_solver::update_column_type_and_bound_with_no_ub(unsigned j, lp::lconstraint_kind kind, const mpq& right_side, unsigned constraint_index) { + void lar_solver::update_column_type_and_bound_with_no_ub(unsigned j, lp::lconstraint_kind kind, const mpq& right_side, u_dependency* dep) { SASSERT(!column_has_upper_bound(j)); if (column_has_lower_bound(j)) { - update_bound_with_no_ub_lb(j, kind, right_side, constraint_index); + update_bound_with_no_ub_lb(j, kind, right_side, dep); } else { - update_bound_with_no_ub_no_lb(j, kind, right_side, constraint_index); + update_bound_with_no_ub_no_lb(j, kind, right_side, dep); } } // clang-format on - void lar_solver::update_bound_with_ub_lb(var_index j, lconstraint_kind kind, const mpq& right_side, constraint_index ci) { + void lar_solver::update_bound_with_ub_lb(var_index j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep) { lp_assert(column_has_lower_bound(j) && column_has_upper_bound(j)); lp_assert(m_mpq_lar_core_solver.m_column_types[j] == column_type::boxed || m_mpq_lar_core_solver.m_column_types[j] == column_type::fixed); @@ -1889,7 +1940,7 @@ namespace lp { } if (up >= m_mpq_lar_core_solver.m_r_upper_bounds[j]) return; m_mpq_lar_core_solver.m_r_upper_bounds[j] = up; - set_upper_bound_witness(j, ci); + set_upper_bound_witness(j, dep); insert_to_columns_with_changed_bounds(j); break; } @@ -1904,7 +1955,7 @@ namespace lp { return; } m_mpq_lar_core_solver.m_r_lower_bounds[j] = low; - set_lower_bound_witness(j, ci); + set_lower_bound_witness(j, dep); m_mpq_lar_core_solver.m_column_types[j] = (low == m_mpq_lar_core_solver.m_r_upper_bounds[j] ? column_type::fixed : column_type::boxed); insert_to_columns_with_changed_bounds(j); break; @@ -1914,8 +1965,8 @@ namespace lp { if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j] || v < m_mpq_lar_core_solver.m_r_lower_bounds[j]) { set_infeasible_column(j); } - set_upper_bound_witness(j, ci); - set_lower_bound_witness(j, ci); + set_upper_bound_witness(j, dep); + set_lower_bound_witness(j, dep); m_mpq_lar_core_solver.m_r_upper_bounds[j] = m_mpq_lar_core_solver.m_r_lower_bounds[j] = v; break; } @@ -1928,7 +1979,7 @@ namespace lp { } } // clang-format off - void lar_solver::update_bound_with_no_ub_lb(var_index j, lconstraint_kind kind, const mpq& right_side, constraint_index ci) { + void lar_solver::update_bound_with_no_ub_lb(var_index j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep) { lp_assert(column_has_lower_bound(j) && !column_has_upper_bound(j)); lp_assert(m_mpq_lar_core_solver.m_column_types[j] == column_type::lower_bound); @@ -1942,7 +1993,7 @@ namespace lp { set_infeasible_column(j); } m_mpq_lar_core_solver.m_r_upper_bounds[j] = up; - set_upper_bound_witness(j, ci); + set_upper_bound_witness(j, dep); m_mpq_lar_core_solver.m_column_types[j] = (up == m_mpq_lar_core_solver.m_r_lower_bounds[j] ? column_type::fixed : column_type::boxed); insert_to_columns_with_changed_bounds(j); break; @@ -1955,7 +2006,7 @@ namespace lp { return; } m_mpq_lar_core_solver.m_r_lower_bounds[j] = low; - set_lower_bound_witness(j, ci); + set_lower_bound_witness(j, dep); insert_to_columns_with_changed_bounds(j); break; } @@ -1965,8 +2016,8 @@ namespace lp { set_infeasible_column(j); } - set_upper_bound_witness(j, ci); - set_lower_bound_witness(j, ci); + set_upper_bound_witness(j, dep); + set_lower_bound_witness(j, dep); m_mpq_lar_core_solver.m_r_upper_bounds[j] = m_mpq_lar_core_solver.m_r_lower_bounds[j] = v; m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed; break; @@ -1977,7 +2028,7 @@ namespace lp { } } // clang-format off - void lar_solver::update_bound_with_ub_no_lb(var_index j, lconstraint_kind kind, const mpq& right_side, constraint_index ci) { + void lar_solver::update_bound_with_ub_no_lb(var_index j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep) { lp_assert(!column_has_lower_bound(j) && column_has_upper_bound(j)); lp_assert(m_mpq_lar_core_solver.m_column_types[j] == column_type::upper_bound); mpq y_of_bound(0); @@ -1989,7 +2040,7 @@ namespace lp { auto up = numeric_pair(right_side, y_of_bound); if (up >= m_mpq_lar_core_solver.m_r_upper_bounds[j]) return; m_mpq_lar_core_solver.m_r_upper_bounds[j] = up; - set_upper_bound_witness(j, ci); + set_upper_bound_witness(j, dep); insert_to_columns_with_changed_bounds(j); } break; @@ -2002,7 +2053,7 @@ namespace lp { set_infeasible_column(j); } m_mpq_lar_core_solver.m_r_lower_bounds[j] = low; - set_lower_bound_witness(j, ci); + set_lower_bound_witness(j, dep); m_mpq_lar_core_solver.m_column_types[j] = (low == m_mpq_lar_core_solver.m_r_upper_bounds[j] ? column_type::fixed : column_type::boxed); insert_to_columns_with_changed_bounds(j); @@ -2015,8 +2066,8 @@ namespace lp { set_infeasible_column(j); } - set_upper_bound_witness(j, ci); - set_lower_bound_witness(j, ci); + set_upper_bound_witness(j, dep); + set_lower_bound_witness(j, dep); m_mpq_lar_core_solver.m_r_upper_bounds[j] = m_mpq_lar_core_solver.m_r_lower_bounds[j] = v; m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed; break; @@ -2027,7 +2078,7 @@ namespace lp { } } // clang-format on - void lar_solver::update_bound_with_no_ub_no_lb(var_index j, lconstraint_kind kind, const mpq& right_side, constraint_index ci) { + void lar_solver::update_bound_with_no_ub_no_lb(var_index j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep) { lp_assert(!column_has_lower_bound(j) && !column_has_upper_bound(j)); mpq y_of_bound(0); @@ -2037,7 +2088,7 @@ namespace lp { case LE: { auto up = numeric_pair(right_side, y_of_bound); m_mpq_lar_core_solver.m_r_upper_bounds[j] = up; - set_upper_bound_witness(j, ci); + set_upper_bound_witness(j, dep); m_mpq_lar_core_solver.m_column_types[j] = column_type::upper_bound; } break; case GT: @@ -2045,14 +2096,14 @@ namespace lp { case GE: { auto low = numeric_pair(right_side, y_of_bound); m_mpq_lar_core_solver.m_r_lower_bounds[j] = low; - set_lower_bound_witness(j, ci); + set_lower_bound_witness(j, dep); m_mpq_lar_core_solver.m_column_types[j] = column_type::lower_bound; } break; case EQ: { auto v = numeric_pair(right_side, zero_of_type()); - set_upper_bound_witness(j, ci); - set_lower_bound_witness(j, ci); + set_upper_bound_witness(j, dep); + set_lower_bound_witness(j, dep); m_mpq_lar_core_solver.m_r_upper_bounds[j] = m_mpq_lar_core_solver.m_r_lower_bounds[j] = v; m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed; break; @@ -2165,7 +2216,7 @@ namespace lp { return true; } - bool lar_solver::get_equality_and_right_side_for_term_on_current_x(tv const& t, mpq& rs, constraint_index& ci, bool& upper_bound) const { + bool lar_solver::get_equality_and_right_side_for_term_on_current_x(tv const& t, mpq& rs, u_dependency*& ci, bool& upper_bound) const { lp_assert(t.is_term()); unsigned j; bool is_int; diff --git a/src/math/lp/lar_solver.h b/src/math/lp/lar_solver.h index 7c1f9575b..4b60d7d62 100644 --- a/src/math/lp/lar_solver.h +++ b/src/math/lp/lar_solver.h @@ -47,6 +47,8 @@ namespace lp { class int_branch; class int_solver; + + class lar_solver : public column_namer { struct term_hasher { std::size_t operator()(const lar_term& t) const { @@ -88,6 +90,8 @@ class lar_solver : public column_namer { indexed_uint_set m_columns_with_changed_bounds; indexed_uint_set m_touched_rows; unsigned_vector m_row_bounds_to_replay; + u_dependency_manager m_dependencies; + svector m_tmp_dependencies; indexed_uint_set m_basic_columns_with_changed_cost; // these are basic columns with the value changed, so the corresponding row in the tableau @@ -136,14 +140,15 @@ class lar_solver : public column_namer { inline void clear_columns_with_changed_bounds() { m_columns_with_changed_bounds.reset(); } void insert_to_columns_with_changed_bounds(unsigned j); - void update_column_type_and_bound_check_on_equal(unsigned j, lconstraint_kind kind, const mpq& right_side, constraint_index constr_index, unsigned&); - void update_column_type_and_bound(unsigned j, lconstraint_kind kind, const mpq& right_side, constraint_index constr_index); - void update_column_type_and_bound_with_ub(var_index j, lconstraint_kind kind, const mpq& right_side, constraint_index constr_index); - void update_column_type_and_bound_with_no_ub(var_index j, lconstraint_kind kind, const mpq& right_side, constraint_index constr_index); - void update_bound_with_ub_lb(var_index j, lconstraint_kind kind, const mpq& right_side, constraint_index constr_index); - void update_bound_with_no_ub_lb(var_index j, lconstraint_kind kind, const mpq& right_side, constraint_index constr_index); - void update_bound_with_ub_no_lb(var_index j, lconstraint_kind kind, const mpq& right_side, constraint_index constr_index); - void update_bound_with_no_ub_no_lb(var_index j, lconstraint_kind kind, const mpq& right_side, constraint_index constr_index); + void update_column_type_and_bound_check_on_equal(unsigned j, const mpq& right_side, constraint_index ci, unsigned&); + void update_column_type_and_bound(unsigned j, const mpq& right_side, constraint_index ci); + void update_column_type_and_bound(unsigned j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep); + void update_column_type_and_bound_with_ub(var_index j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep); + void update_column_type_and_bound_with_no_ub(var_index j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep); + void update_bound_with_ub_lb(var_index j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep); + void update_bound_with_no_ub_lb(var_index j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep); + void update_bound_with_ub_no_lb(var_index j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep); + void update_bound_with_no_ub_no_lb(var_index j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep); void register_in_fixed_var_table(unsigned, unsigned&); void remove_non_fixed_from_fixed_var_table(); constraint_index add_var_bound_on_constraint_for_term(var_index j, lconstraint_kind kind, const mpq& right_side); @@ -186,8 +191,8 @@ class lar_solver : public column_namer { bool maximize_term_on_corrected_r_solver(lar_term& term, impq& term_max); void pop_core_solver_params(); void pop_core_solver_params(unsigned k); - void set_upper_bound_witness(var_index j, constraint_index ci); - void set_lower_bound_witness(var_index j, constraint_index ci); + void set_upper_bound_witness(var_index j, u_dependency* ci); + void set_lower_bound_witness(var_index j, u_dependency* ci); void substitute_terms_in_linear_expression(const vector>& left_side_with_terms, vector>& left_side) const; @@ -283,6 +288,8 @@ class lar_solver : public column_namer { lp_status maximize_term(unsigned j_or_term, impq& term_max); + bool improve_bound(lpvar j, bool is_lower); + inline core_solver_pretty_printer pp(std::ostream& out) const { return core_solver_pretty_printer(m_mpq_lar_core_solver.m_r_solver, out); } @@ -310,9 +317,10 @@ class lar_solver : public column_namer { int a_sign = is_pos(a) ? 1 : -1; int sign = j_sign * a_sign; const ul_pair& ul = m_columns_to_ul_pairs[j]; - auto witness = sign > 0 ? ul.upper_bound_witness() : ul.lower_bound_witness(); - lp_assert(is_valid(witness)); - bp.consume(a, witness); + auto* witness = sign > 0 ? ul.upper_bound_witness() : ul.lower_bound_witness(); + lp_assert(witness); + for (auto ci : flatten(witness)) + bp.consume(a, ci); } } @@ -460,7 +468,20 @@ class lar_solver : public column_namer { return m_mpq_lar_core_solver.m_r_solver.column_has_lower_bound(j); } - inline constraint_index get_column_upper_bound_witness(unsigned j) const { + svector const& flatten(u_dependency* d) { + m_tmp_dependencies.reset(); + m_dependencies.linearize(d, m_tmp_dependencies); + return m_tmp_dependencies; + } + + void push_explanation(u_dependency* d, explanation& ex) { + for (auto ci : flatten(d)) + ex.push_back(ci); + } + + u_dependency_manager& dep_manager() { return m_dependencies; } + + inline u_dependency* get_column_upper_bound_witness(unsigned j) const { if (tv::is_term(j)) { j = m_var_register.external_to_local(j); } @@ -474,8 +495,8 @@ class lar_solver : public column_namer { inline const impq& get_lower_bound(column_index j) const { return m_mpq_lar_core_solver.m_r_solver.m_lower_bounds[j]; } - bool has_lower_bound(var_index var, constraint_index& ci, mpq& value, bool& is_strict) const; - bool has_upper_bound(var_index var, constraint_index& ci, mpq& value, bool& is_strict) const; + bool has_lower_bound(var_index var, u_dependency*& ci, mpq& value, bool& is_strict) const; + bool has_upper_bound(var_index var, u_dependency*& ci, mpq& value, bool& is_strict) const; bool has_value(var_index var, mpq& value) const; bool fetch_normalized_term_column(const lar_term& t, std::pair&) const; unsigned map_term_index_to_column_index(unsigned j) const; @@ -530,15 +551,15 @@ class lar_solver : public column_namer { std::pair add_equality(lpvar j, lpvar k); - inline void get_bound_constraint_witnesses_for_column(unsigned j, constraint_index& lc, constraint_index& uc) const { + u_dependency* get_bound_constraint_witnesses_for_column(unsigned j) { const ul_pair& ul = m_columns_to_ul_pairs[j]; - lc = ul.lower_bound_witness(); - uc = ul.upper_bound_witness(); + return m_dependencies.mk_join(ul.lower_bound_witness(), ul.upper_bound_witness()); } inline constraint_set const& constraints() const { return m_constraints; } void push(); void pop(); - inline constraint_index get_column_lower_bound_witness(unsigned j) const { + + inline u_dependency* get_column_lower_bound_witness(unsigned j) const { if (tv::is_term(j)) { j = m_var_register.external_to_local(j); } @@ -601,7 +622,7 @@ class lar_solver : public column_namer { inline const column_strip& get_column(unsigned i) const { return A_r().m_columns[i]; } bool row_is_correct(unsigned i) const; bool ax_is_correct() const; - bool get_equality_and_right_side_for_term_on_current_x(tv const& t, mpq& rs, constraint_index& ci, bool& upper_bound) const; + bool get_equality_and_right_side_for_term_on_current_x(tv const& t, mpq& rs, u_dependency*& ci, bool& upper_bound) const; bool var_is_int(var_index v) const; inline const vector& r_heading() const { return m_mpq_lar_core_solver.m_r_heading; } inline const vector& r_basis() const { return m_mpq_lar_core_solver.r_basis(); } diff --git a/src/math/lp/lp_bound_propagator.h b/src/math/lp/lp_bound_propagator.h index 46670ad91..f676b0e1d 100644 --- a/src/math/lp/lp_bound_propagator.h +++ b/src/math/lp/lp_bound_propagator.h @@ -228,33 +228,32 @@ class lp_bound_propagator { return lp().column_is_int(j); } - void explain_fixed_in_row(unsigned row, explanation& ex) const { + void explain_fixed_in_row(unsigned row, explanation& ex) { TRACE("eq", tout << lp().get_row(row) << std::endl); for (const auto& c : lp().get_row(row)) if (lp().is_fixed(c.var())) explain_fixed_column(c.var(), ex); } - unsigned explain_fixed_in_row_and_get_base(unsigned row, explanation& ex) const { + unsigned explain_fixed_in_row_and_get_base(unsigned row, explanation& ex) { unsigned base = UINT_MAX; TRACE("eq", tout << lp().get_row(row) << std::endl); for (const auto& c : lp().get_row(row)) { if (lp().is_fixed(c.var())) { explain_fixed_column(c.var(), ex); - } else if (lp().is_base(c.var())) { + } + else if (lp().is_base(c.var())) { base = c.var(); } } return base; } - void explain_fixed_column(unsigned j, explanation& ex) const { + void explain_fixed_column(unsigned j, explanation& ex) { SASSERT(column_is_fixed(j)); - constraint_index lc, uc; - lp().get_bound_constraint_witnesses_for_column(j, lc, uc); - ex.push_back(lc); - if (lc != uc) - ex.push_back(uc); + auto* deps = lp().get_bound_constraint_witnesses_for_column(j); + for (auto ci : lp().flatten(deps)) + ex.push_back(ci); } #ifdef Z3DEBUG bool all_fixed_in_row(unsigned row) const { diff --git a/src/math/lp/lp_primal_core_solver.h b/src/math/lp/lp_primal_core_solver.h index c074da16a..9871ec691 100644 --- a/src/math/lp/lp_primal_core_solver.h +++ b/src/math/lp/lp_primal_core_solver.h @@ -138,95 +138,98 @@ namespace lp { return false; } -/** - * Return the number of base non-free variables depending on the column j, - * different from bj, - * but take the min with the (bound+1). - * This function is used to select the pivot variable. -*/ - unsigned get_num_of_not_free_basic_dependent_vars( - unsigned j, unsigned bound, unsigned bj) const { // consider looking at the signs here: todo - unsigned r = 0; - for (const auto &cc : this->m_A.m_columns[j]) { - unsigned basic_for_row = this->m_basis[cc.var()]; - if (basic_for_row == bj) - continue; - - // std::cout << this->m_A.m_rows[cc.var()] << std::endl; - if (this->m_column_types[basic_for_row] != column_type::free_column) - if (r++ > bound) return r; - } - return r; - } - int find_beneficial_entering_in_row_tableau_rows_bland_mode(int i, T &a_ent) { - int j = -1; - unsigned bj = this->m_basis[i]; - bool bj_needs_to_grow = needs_to_grow(bj); - for (const row_cell &rc : this->m_A.m_rows[i]) { - if (rc.var() == bj) - continue; - if (bj_needs_to_grow) { - if (!monoid_can_decrease(rc)) - continue; - } else { - if (!monoid_can_increase(rc)) + /** + * Return the number of base non-free variables depending on the column j, + * different from bj, + * but take the min with the (bound+1). + * This function is used to select the pivot variable. + */ + unsigned get_num_of_not_free_basic_dependent_vars(unsigned j, unsigned bound, unsigned bj) const { + // consider looking at the signs here: todo + unsigned r = 0; + for (const auto &cc : this->m_A.m_columns[j]) { + unsigned basic_for_row = this->m_basis[cc.var()]; + if (basic_for_row == bj) continue; + + // std::cout << this->m_A.m_rows[cc.var()] << std::endl; + if (this->m_column_types[basic_for_row] != column_type::free_column) + if (r++ > bound) return r; } - if (rc.var() < static_cast(j)) { - j = rc.var(); - a_ent = rc.coeff(); - } - } - if (j == -1) - m_inf_row_index_for_tableau = i; - return j; - } - int find_beneficial_entering_tableau_rows(int i, T &a_ent) { - if (m_bland_mode_tableau) - return find_beneficial_entering_in_row_tableau_rows_bland_mode(i, a_ent); - // a short row produces short infeasibility explanation and benefits at - // least one pivot operation - int choice = -1; - int nchoices = 0; - unsigned min_non_free_so_far = -1; - unsigned best_col_sz = -1; - unsigned bj = this->m_basis[i]; - bool bj_needs_to_grow = needs_to_grow(bj); - for (unsigned k = 0; k < this->m_A.m_rows[i].size(); k++) { - const row_cell &rc = this->m_A.m_rows[i][k]; - unsigned j = rc.var(); - if (j == bj) - continue; - if (bj_needs_to_grow) { - if (!monoid_can_decrease(rc)) - continue; - } else { - if (!monoid_can_increase(rc)) - continue; - } - unsigned not_free = get_num_of_not_free_basic_dependent_vars(j, min_non_free_so_far, bj); - unsigned col_sz = this->m_A.m_columns[j].size(); - if (not_free < min_non_free_so_far || (not_free == min_non_free_so_far && col_sz < best_col_sz)) { - min_non_free_so_far = not_free; - best_col_sz = this->m_A.m_columns[j].size(); - choice = k; - nchoices = 1; - } else if (not_free == min_non_free_so_far && - col_sz == best_col_sz) { - if (this->m_settings.random_next(++nchoices) == 0) { - choice = k; - } - } + return r; } - if (choice == -1) { - m_inf_row_index_for_tableau = i; - return -1; + int find_beneficial_entering_in_row_tableau_rows_bland_mode(int i, T &a_ent) { + int j = -1; + unsigned bj = this->m_basis[i]; + bool bj_needs_to_grow = needs_to_grow(bj); + for (const row_cell &rc : this->m_A.m_rows[i]) { + if (rc.var() == bj) + continue; + if (bj_needs_to_grow) { + if (!monoid_can_decrease(rc)) + continue; + } + else { + if (!monoid_can_increase(rc)) + continue; + } + if (rc.var() < static_cast(j)) { + j = rc.var(); + a_ent = rc.coeff(); + } + } + if (j == -1) + m_inf_row_index_for_tableau = i; + return j; + } + + int find_beneficial_entering_tableau_rows(int i, T &a_ent) { + if (m_bland_mode_tableau) + return find_beneficial_entering_in_row_tableau_rows_bland_mode(i, a_ent); + // a short row produces short infeasibility explanation and benefits at + // least one pivot operation + int choice = -1; + int nchoices = 0; + unsigned min_non_free_so_far = -1; + unsigned best_col_sz = -1; + unsigned bj = this->m_basis[i]; + bool bj_needs_to_grow = needs_to_grow(bj); + for (unsigned k = 0; k < this->m_A.m_rows[i].size(); k++) { + const row_cell &rc = this->m_A.m_rows[i][k]; + unsigned j = rc.var(); + if (j == bj) + continue; + if (bj_needs_to_grow) { + if (!monoid_can_decrease(rc)) + continue; + } else { + if (!monoid_can_increase(rc)) + continue; + } + unsigned not_free = get_num_of_not_free_basic_dependent_vars(j, min_non_free_so_far, bj); + unsigned col_sz = this->m_A.m_columns[j].size(); + if (not_free < min_non_free_so_far || (not_free == min_non_free_so_far && col_sz < best_col_sz)) { + min_non_free_so_far = not_free; + best_col_sz = this->m_A.m_columns[j].size(); + choice = k; + nchoices = 1; + } + else if (not_free == min_non_free_so_far && + col_sz == best_col_sz) { + if (this->m_settings.random_next(++nchoices) == 0) + choice = k; + } + } + + if (choice == -1) { + m_inf_row_index_for_tableau = i; + return -1; + } + const row_cell &rc = this->m_A.m_rows[i][choice]; + a_ent = rc.coeff(); + return rc.var(); } - const row_cell &rc = this->m_A.m_rows[i][choice]; - a_ent = rc.coeff(); - return rc.var(); - } bool try_jump_to_another_bound_on_entering(unsigned entering, const X &theta, X &t, bool &unlimited); diff --git a/src/math/lp/lp_types.h b/src/math/lp/lp_types.h index ac046e1a1..c69dbe10b 100644 --- a/src/math/lp/lp_types.h +++ b/src/math/lp/lp_types.h @@ -22,6 +22,8 @@ Revision History: #include #include #include "util/debug.h" +#include "util/dependency.h" + namespace nla { class core; } @@ -36,6 +38,7 @@ typedef unsigned lpvar; const lpvar null_lpvar = UINT_MAX; const constraint_index null_ci = UINT_MAX; + class column_index { unsigned m_index; friend class lar_solver; diff --git a/src/math/lp/monomial_bounds.cpp b/src/math/lp/monomial_bounds.cpp index 1ed0956dc..a89d27a9b 100644 --- a/src/math/lp/monomial_bounds.cpp +++ b/src/math/lp/monomial_bounds.cpp @@ -149,22 +149,22 @@ namespace nla { } void monomial_bounds::var2interval(lpvar v, scoped_dep_interval& i) { - lp::constraint_index ci; + u_dependency* d = nullptr; rational bound; bool is_strict; - if (c().has_lower_bound(v, ci, bound, is_strict)) { + if (c().has_lower_bound(v, d, bound, is_strict)) { dep.set_lower_is_open(i, is_strict); dep.set_lower(i, bound); - dep.set_lower_dep(i, dep.mk_leaf(ci)); + dep.set_lower_dep(i, d); dep.set_lower_is_inf(i, false); } else { dep.set_lower_is_inf(i, true); } - if (c().has_upper_bound(v, ci, bound, is_strict)) { + if (c().has_upper_bound(v, d, bound, is_strict)) { dep.set_upper_is_open(i, is_strict); dep.set_upper(i, bound); - dep.set_upper_dep(i, dep.mk_leaf(ci)); + dep.set_upper_dep(i, d); dep.set_upper_is_inf(i, false); } else { diff --git a/src/math/lp/nla_common.cpp b/src/math/lp/nla_common.cpp index 00f2ff4ee..eceacbcb4 100644 --- a/src/math/lp/nla_common.cpp +++ b/src/math/lp/nla_common.cpp @@ -60,11 +60,9 @@ unsigned common::random() { } void common::add_deps_of_fixed(lpvar j, u_dependency*& dep) { - unsigned lc, uc; - auto& dep_manager = c().m_intervals.get_dep_intervals().dep_manager(); - c().m_lar_solver.get_bound_constraint_witnesses_for_column(j, lc, uc); - dep = dep_manager.mk_join(dep, dep_manager.mk_leaf(lc)); - dep = dep_manager.mk_join(dep, dep_manager.mk_leaf(uc)); + auto& dm = c().lra.dep_manager(); + auto* deps = c().lra.get_bound_constraint_witnesses_for_column(j); + dep = dm.mk_join(dep, deps); } @@ -73,7 +71,7 @@ nex * common::nexvar(const rational & coeff, lpvar j, nex_creator& cn, u_depende SASSERT(!coeff.is_zero()); if (c().params().arith_nl_horner_subs_fixed() == 1 && c().var_is_fixed(j)) { add_deps_of_fixed(j, dep); - return cn.mk_scalar(coeff * c().m_lar_solver.column_lower_bound(j).x); + return cn.mk_scalar(coeff * c().lra.column_lower_bound(j).x); } if (c().params().arith_nl_horner_subs_fixed() == 2 && c().var_is_fixed_to_zero(j)) { add_deps_of_fixed(j, dep); @@ -91,7 +89,7 @@ nex * common::nexvar(const rational & coeff, lpvar j, nex_creator& cn, u_depende for (lpvar k : m.vars()) { if (c().params().arith_nl_horner_subs_fixed() == 1 && c().var_is_fixed(k)) { add_deps_of_fixed(k, dep); - mf *= c().m_lar_solver.column_lower_bound(k).x; + mf *= c().lra.column_lower_bound(k).x; } else if (c().params().arith_nl_horner_subs_fixed() == 2 && c().var_is_fixed_to_zero(k)) { dep = initial_dep; diff --git a/src/math/lp/nla_common.h b/src/math/lp/nla_common.h index 1302c3909..731b4766e 100644 --- a/src/math/lp/nla_common.h +++ b/src/math/lp/nla_common.h @@ -80,21 +80,6 @@ struct common { bool check_monic(const monic&) const; unsigned random(); void add_deps_of_fixed(lpvar j, u_dependency*& dep); - class ci_value_manager { - public: - void inc_ref(lp::constraint_index const & v) { - } - - void dec_ref(lp::constraint_index const & v) { - } - }; - - struct u_dependency_config { - typedef ci_value_manager value_manager; - typedef region allocator; - static const bool ref_count = false; - typedef lp::constraint_index value; - }; nex* nexvar(const rational& coeff, lpvar j, nex_creator&, u_dependency*&); template diff --git a/src/math/lp/nla_core.cpp b/src/math/lp/nla_core.cpp index 54f955760..49691d6b8 100644 --- a/src/math/lp/nla_core.cpp +++ b/src/math/lp/nla_core.cpp @@ -23,7 +23,7 @@ typedef lp::lar_term term; core::core(lp::lar_solver& s, params_ref const& p, reslimit & lim) : m_evars(), - m_lar_solver(s), + lra(s), m_reslim(lim), m_params(p), m_tangents(this), @@ -69,7 +69,7 @@ lp::lar_term core::subs_terms_to_columns(const lp::lar_term& t) const { for (lp::lar_term::ival p : t) { lpvar j = p.column(); if (lp::tv::is_term(j)) - j = m_lar_solver.map_term_index_to_column_index(j); + j = lra.map_term_index_to_column_index(j); r.add_monomial(p.coeff(), j); } return r; @@ -133,7 +133,7 @@ void core::add_monic(lpvar v, unsigned sz, lpvar const* vs) { for (unsigned i = 0; i < sz; i++) { lpvar j = vs[i]; if (lp::tv::is_term(j)) - j = m_lar_solver.map_term_index_to_column_index(j); + j = lra.map_term_index_to_column_index(j); m_add_buffer[i] = j; } m_emons.add(v, m_add_buffer); @@ -154,7 +154,7 @@ void core::pop(unsigned n) { rational core::product_value(const monic& m) const { rational r(1); for (auto j : m.vars()) { - r *= m_lar_solver.get_column_value(j).x; + r *= lra.get_column_value(j).x; } return r; } @@ -167,10 +167,10 @@ bool core::check_monic(const monic& m) const { if (!is_relevant(m.var())) return true; #endif - if (m_lar_solver.column_is_int(m.var()) && !m_lar_solver.get_column_value(m.var()).is_int()) + if (lra.column_is_int(m.var()) && !lra.get_column_value(m.var()).is_int()) return true; - bool ret = product_value(m) == m_lar_solver.get_column_value(m.var()).x; + bool ret = product_value(m) == lra.get_column_value(m.var()).x; CTRACE("nla_solver_check_monic", !ret, print_monic(m, tout) << '\n';); return ret; } @@ -182,7 +182,7 @@ std::ostream& core::print_product(const T & m, std::ostream& out) const { for (lpvar v : m) { if (!first) out << "*"; else first = false; if (lp_settings().print_external_var_name()) - out << "(" << m_lar_solver.get_variable_name(v) << "=" << val(v) << ")"; + out << "(" << lra.get_variable_name(v) << "=" << val(v) << ")"; else out << "(j" << v << " = " << val(v) << ")"; @@ -228,7 +228,7 @@ std::ostream & core::print_factor_with_vars(const factor& f, std::ostream& out) std::ostream& core::print_monic(const monic& m, std::ostream& out) const { if (lp_settings().print_external_var_name()) - out << "([" << m.var() << "] = " << m_lar_solver.get_variable_name(m.var()) << " = " << val(m.var()) << " = "; + out << "([" << m.var() << "] = " << lra.get_variable_name(m.var()) << " = " << val(m.var()) << " = "; else out << "(j" << m.var() << " = " << val(m.var()) << " = "; print_product(m.vars(), out) << ")\n"; @@ -271,7 +271,7 @@ std::ostream& core::print_explanation(const lp::explanation& exp, std::ostream& unsigned i = 0; for (auto p : exp) { out << "(" << p.ci() << ")"; - m_lar_solver.constraints().display(out, [this](lpvar j) { return var_str(j);}, p.ci()); + lra.constraints().display(out, [this](lpvar j) { return var_str(j);}, p.ci()); if (++i < exp.size()) out << " "; } @@ -316,21 +316,20 @@ bool core::explain_lower_bound(const lp::lar_term& t, const rational& rs, lp::ex bool core::explain_coeff_lower_bound(const lp::lar_term::ival& p, rational& bound, lp::explanation& e) const { const rational& a = p.coeff(); SASSERT(!a.is_zero()); - unsigned c; // the index for the lower or the upper bound if (a.is_pos()) { - unsigned c = m_lar_solver.get_column_lower_bound_witness(p.column()); - if (c + 1 == 0) + auto* dep = lra.get_column_lower_bound_witness(p.column()); + if (!dep) return false; - bound = a * m_lar_solver.get_lower_bound(p.column()).x; - e.push_back(c); + bound = a * lra.get_lower_bound(p.column()).x; + lra.push_explanation(dep, e); return true; } // a.is_neg() - c = m_lar_solver.get_column_upper_bound_witness(p.column()); - if (c + 1 == 0) + auto* dep = lra.get_column_upper_bound_witness(p.column()); + if (!dep) return false; - bound = a * m_lar_solver.get_upper_bound(p.column()).x; - e.push_back(c); + bound = a * lra.get_upper_bound(p.column()).x; + lra.push_explanation(dep, e); return true; } @@ -338,21 +337,20 @@ bool core::explain_coeff_upper_bound(const lp::lar_term::ival& p, rational& boun const rational& a = p.coeff(); lpvar j = p.column(); SASSERT(!a.is_zero()); - unsigned c; // the index for the lower or the upper bound if (a.is_neg()) { - unsigned c = m_lar_solver.get_column_lower_bound_witness(j); - if (c + 1 == 0) + auto *dep = lra.get_column_lower_bound_witness(j); + if (!dep) return false; - bound = a * m_lar_solver.get_lower_bound(j).x; - e.push_back(c); + bound = a * lra.get_lower_bound(j).x; + lra.push_explanation(dep, e); return true; } // a.is_pos() - c = m_lar_solver.get_column_upper_bound_witness(j); - if (c + 1 == 0) + auto* dep = lra.get_column_upper_bound_witness(j); + if (!dep) return false; - bound = a * m_lar_solver.get_upper_bound(j).x; - e.push_back(c); + bound = a * lra.get_upper_bound(j).x; + lra.push_explanation(dep, e); return true; } @@ -413,12 +411,12 @@ bool core::explain_by_equiv(const lp::lar_term& t, lp::explanation& e) const { return false; m_evars.explain(signed_var(i, false), signed_var(j, sign), e); - TRACE("nla_solver", tout << "explained :"; m_lar_solver.print_term_as_indices(t, tout);); + TRACE("nla_solver", tout << "explained :"; lra.print_term_as_indices(t, tout);); return true; } void core::mk_ineq_no_expl_check(new_lemma& lemma, lp::lar_term& t, llc cmp, const rational& rs) { - TRACE("nla_solver_details", m_lar_solver.print_term_as_indices(t, tout << "t = ");); + TRACE("nla_solver_details", lra.print_term_as_indices(t, tout << "t = ");); lemma |= ineq(cmp, t, rs); CTRACE("nla_solver", ineq_holds(ineq(cmp, t, rs)), print_ineq(ineq(cmp, t, rs), tout) << "\n";); SASSERT(!ineq_holds(ineq(cmp, t, rs))); @@ -485,18 +483,18 @@ int core::vars_sign(const svector& v) { } bool core::has_upper_bound(lpvar j) const { - return m_lar_solver.column_has_upper_bound(j); + return lra.column_has_upper_bound(j); } bool core::has_lower_bound(lpvar j) const { - return m_lar_solver.column_has_lower_bound(j); + return lra.column_has_lower_bound(j); } const rational& core::get_upper_bound(unsigned j) const { - return m_lar_solver.get_upper_bound(j).x; + return lra.get_upper_bound(j).x; } const rational& core::get_lower_bound(unsigned j) const { - return m_lar_solver.get_lower_bound(j).x; + return lra.get_lower_bound(j).x; } bool core::zero_is_an_inner_point_of_bounds(lpvar j) const { @@ -540,26 +538,26 @@ bool core::sign_contradiction(const monic& m) const { bool core::var_is_fixed_to_zero(lpvar j) const { return - m_lar_solver.column_is_fixed(j) && - m_lar_solver.get_lower_bound(j) == lp::zero_of_type(); + lra.column_is_fixed(j) && + lra.get_lower_bound(j) == lp::zero_of_type(); } bool core::var_is_fixed_to_val(lpvar j, const rational& v) const { return - m_lar_solver.column_is_fixed(j) && - m_lar_solver.get_lower_bound(j) == lp::impq(v); + lra.column_is_fixed(j) && + lra.get_lower_bound(j) == lp::impq(v); } bool core::var_is_fixed(lpvar j) const { - return m_lar_solver.column_is_fixed(j); + return lra.column_is_fixed(j); } bool core::var_is_free(lpvar j) const { - return m_lar_solver.column_is_free(j); + return lra.column_is_free(j); } std::ostream & core::print_ineq(const ineq & in, std::ostream & out) const { - m_lar_solver.print_term_as_indices(in.term(), out); + lra.print_term_as_indices(in.term(), out); out << " " << lconstraint_kind_string(in.cmp()) << " " << in.rs(); return out; } @@ -569,14 +567,14 @@ std::ostream & core::print_var(lpvar j, std::ostream & out) const { print_monic(m_emons[j], out); } - m_lar_solver.print_column_info(j, out); + lra.print_column_info(j, out); signed_var jr = m_evars.find(j); out << "root="; if (jr.sign()) { out << "-"; } - out << m_lar_solver.get_variable_name(jr.var()) << "\n"; + out << lra.get_variable_name(jr.var()) << "\n"; return out; } @@ -643,11 +641,11 @@ void core::trace_print_monic_and_factorization(const monic& rm, const factorizat bool core::var_has_positive_lower_bound(lpvar j) const { - return m_lar_solver.column_has_lower_bound(j) && m_lar_solver.get_lower_bound(j) > lp::zero_of_type(); + return lra.column_has_lower_bound(j) && lra.get_lower_bound(j) > lp::zero_of_type(); } bool core::var_has_negative_upper_bound(lpvar j) const { - return m_lar_solver.column_has_upper_bound(j) && m_lar_solver.get_upper_bound(j) < lp::zero_of_type(); + return lra.column_has_upper_bound(j) && lra.get_upper_bound(j) < lp::zero_of_type(); } bool core::var_is_separated_from_zero(lpvar j) const { @@ -684,10 +682,10 @@ template bool core::mon_has_zero(const unsigned_vector& product lp::lp_settings& core::lp_settings() { - return m_lar_solver.settings(); + return lra.settings(); } const lp::lp_settings& core::lp_settings() const { - return m_lar_solver.settings(); + return lra.settings(); } unsigned core::random() { return lp_settings().random_next(); } @@ -695,7 +693,7 @@ unsigned core::random() { return lp_settings().random_next(); } // we look for octagon constraints here, with a left part +-x +- y void core::collect_equivs() { - const lp::lar_solver& s = m_lar_solver; + const lp::lar_solver& s = lra; for (unsigned i = 0; i < s.terms().size(); i++) { if (!s.term_is_used_as_row(i)) @@ -737,7 +735,7 @@ bool core::is_octagon_term(const lp::lar_term& t, bool & sign, lpvar& i, lpvar & return true; } -void core::add_equivalence_maybe(const lp::lar_term *t, lpci c0, lpci c1) { +void core::add_equivalence_maybe(const lp::lar_term* t, u_dependency* c0, u_dependency* c1) { bool sign; lpvar i, j; if (!is_octagon_term(*t, sign, i, j)) @@ -866,7 +864,7 @@ std::unordered_set core::collect_vars(const lemma& l) const { } } for (auto p : l.expl()) { - const auto& c = m_lar_solver.constraints()[p.ci()]; + const auto& c = lra.constraints()[p.ci()]; for (const auto& r : c.coeffs()) { insert_j(r.second); } @@ -1125,8 +1123,8 @@ new_lemma& new_lemma::explain_equiv(lpvar a, lpvar b) { new_lemma& new_lemma::explain_var_separated_from_zero(lpvar j) { SASSERT(c.var_is_separated_from_zero(j)); - if (c.m_lar_solver.column_has_upper_bound(j) && - (c.m_lar_solver.get_upper_bound(j)< lp::zero_of_type())) + if (c.lra.column_has_upper_bound(j) && + (c.lra.get_upper_bound(j)< lp::zero_of_type())) explain_existing_upper_bound(j); else explain_existing_lower_bound(j); @@ -1136,7 +1134,7 @@ new_lemma& new_lemma::explain_var_separated_from_zero(lpvar j) { new_lemma& new_lemma::explain_existing_lower_bound(lpvar j) { SASSERT(c.has_lower_bound(j)); lp::explanation ex; - ex.push_back(c.m_lar_solver.get_column_lower_bound_witness(j)); + c.lra.push_explanation(c.lra.get_column_lower_bound_witness(j), ex); *this &= ex; TRACE("nla_solver", tout << j << ": " << *this << "\n";); return *this; @@ -1145,7 +1143,7 @@ new_lemma& new_lemma::explain_existing_lower_bound(lpvar j) { new_lemma& new_lemma::explain_existing_upper_bound(lpvar j) { SASSERT(c.has_upper_bound(j)); lp::explanation ex; - ex.push_back(c.m_lar_solver.get_column_upper_bound_witness(j)); + c.lra.push_explanation(c.lra.get_column_upper_bound_witness(j), ex); *this &= ex; return *this; } @@ -1155,7 +1153,7 @@ std::ostream& new_lemma::display(std::ostream & out) const { for (auto p : lemma.expl()) { out << "(" << p.ci() << ") "; - c.m_lar_solver.constraints().display(out, [this](lpvar j) { return c.var_str(j);}, p.ci()); + c.lra.constraints().display(out, [this](lpvar j) { return c.var_str(j);}, p.ci()); } out << " ==> "; if (lemma.ineqs().empty()) { @@ -1303,7 +1301,7 @@ bool core::has_real(const monic& m) const { bool core::is_patch_blocked(lpvar u, const lp::impq& ival) const { TRACE("nla_solver", tout << "u = " << u << '\n';); if (m_cautious_patching && - (!m_lar_solver.inside_bounds(u, ival) || (var_is_int(u) && ival.is_int() == false))) { + (!lra.inside_bounds(u, ival) || (var_is_int(u) && ival.is_int() == false))) { TRACE("nla_solver", tout << "u = " << u << " blocked, for feas or integr\n";); return true; // block } @@ -1334,7 +1332,7 @@ bool core::is_patch_blocked(lpvar u, const lp::impq& ival) const { bool core::try_to_patch(const rational& v) { auto is_blocked = [this](lpvar u, const lp::impq& iv) { return is_patch_blocked(u, iv); }; auto change_report = [this](lpvar u) { update_to_refine_of_var(u); }; - return m_lar_solver.try_to_patch(m_patched_var, v, is_blocked, change_report); + return lra.try_to_patch(m_patched_var, v, is_blocked, change_report); } bool in_power(const svector& vs, unsigned l) { @@ -1343,7 +1341,7 @@ bool in_power(const svector& vs, unsigned l) { } bool core::to_refine_is_correct() const { - for (unsigned j = 0; j < m_lar_solver.number_of_vars(); j++) { + for (unsigned j = 0; j < lra.number_of_vars(); j++) { if (!is_monic_var(j)) continue; bool valid = check_monic(emons()[j]); if (valid == m_to_refine.contains(j)) { @@ -1435,17 +1433,17 @@ void core::patch_monomials() { NOT_IMPLEMENTED_YET(); m_cautious_patching = false; patch_monomials_on_to_refine(); - m_lar_solver.push(); + lra.push(); save_tableau(); constrain_nl_in_tableau(); if (solve_tableau() && integrality_holds()) { - m_lar_solver.pop(1); + lra.pop(1); } else { - m_lar_solver.pop(); + lra.pop(); restore_tableau(); - m_lar_solver.clear_inf_heap(); + lra.clear_inf_heap(); } - SASSERT(m_lar_solver.ax_is_correct()); + SASSERT(lra.ax_is_correct()); } void core::constrain_nl_in_tableau() { @@ -1507,11 +1505,11 @@ void core::check_bounded_divisions(vector& l_vec) { lbool core::check(vector& l_vec) { lp_settings().stats().m_nla_calls++; TRACE("nla_solver", tout << "calls = " << lp_settings().stats().m_nla_calls << "\n";); - m_lar_solver.get_rid_of_inf_eps(); + lra.get_rid_of_inf_eps(); m_lemma_vec = &l_vec; - if (!(m_lar_solver.get_status() == lp::lp_status::OPTIMAL || - m_lar_solver.get_status() == lp::lp_status::FEASIBLE)) { - TRACE("nla_solver", tout << "unknown because of the m_lar_solver.m_status = " << m_lar_solver.get_status() << "\n";); + if (!(lra.get_status() == lp::lp_status::OPTIMAL || + lra.get_status() == lp::lp_status::FEASIBLE)) { + TRACE("nla_solver", tout << "unknown because of the lra.m_status = " << lra.get_status() << "\n";); return l_undef; } @@ -1528,6 +1526,9 @@ lbool core::check(vector& l_vec) { if (l_vec.empty() && !done()) m_monomial_bounds(); + + if (l_vec.empty() && !done() && improve_bounds()) + return l_false; if (l_vec.empty() && !done() && run_horner) m_horner.horner_lemmas(); @@ -1635,21 +1636,21 @@ bool core::no_lemmas_hold() const { } lbool core::test_check(vector& l) { - m_lar_solver.set_status(lp::lp_status::OPTIMAL); + lra.set_status(lp::lp_status::OPTIMAL); return check(l); } std::ostream& core::print_terms(std::ostream& out) const { - for (unsigned i = 0; i< m_lar_solver.terms().size(); i++) { + for (unsigned i = 0; i< lra.terms().size(); i++) { unsigned ext = lp::tv::mask_term(i); - if (!m_lar_solver.var_is_registered(ext)) { + if (!lra.var_is_registered(ext)) { out << "term is not registered\n"; continue; } - const lp::lar_term & t = *m_lar_solver.terms()[i]; + const lp::lar_term & t = *lra.terms()[i]; out << "term:"; print_term(t, out) << std::endl; - lpvar j = m_lar_solver.external_to_local(ext); + lpvar j = lra.external_to_local(ext); print_var(j, out); } return out; @@ -1677,7 +1678,7 @@ std::ostream& core::print_term( const lp::lar_term& t, std::ostream& out) const std::unordered_set core::get_vars_of_expr_with_opening_terms(const nex *e ) { auto ret = get_vars_of_expr(e); - auto & ls = m_lar_solver; + auto & ls = lra; svector added; for (auto j : ret) { added.push_back(j); @@ -1685,7 +1686,7 @@ std::unordered_set core::get_vars_of_expr_with_opening_terms(const nex *e for (unsigned i = 0; i < added.size(); ++i) { lpvar j = added[i]; if (ls.column_corresponds_to_term(j)) { - const auto& t = m_lar_solver.get_term(lp::tv::raw(ls.local_to_external(j))); + const auto& t = lra.get_term(lp::tv::raw(ls.local_to_external(j))); for (auto p : t) { if (ret.find(p.column()) == ret.end()) { added.push_back(p.column()); @@ -1705,7 +1706,7 @@ bool core::is_nl_var(lpvar j) const { unsigned core::get_var_weight(lpvar j) const { unsigned k; - switch (m_lar_solver.get_column_type(j)) { + switch (lra.get_column_type(j)) { case lp::column_type::fixed: k = 0; @@ -1734,7 +1735,7 @@ unsigned core::get_var_weight(lpvar j) const { void core::set_active_vars_weights(nex_creator& nc) { - nc.set_number_of_vars(m_lar_solver.column_count()); + nc.set_number_of_vars(lra.column_count()); for (lpvar j : active_var_set()) nc.set_var_weight(j, get_var_weight(j)); } @@ -1744,8 +1745,8 @@ bool core::influences_nl_var(lpvar j) const { j = lp::tv::unmask_term(j); if (is_nl_var(j)) return true; - for (const auto & c : m_lar_solver.A_r().m_columns[j]) { - lpvar basic_in_row = m_lar_solver.r_basis()[c.var()]; + for (const auto & c : lra.A_r().m_columns[j]) { + lpvar basic_in_row = lra.r_basis()[c.var()]; if (is_nl_var(basic_in_row)) return true; } @@ -1762,9 +1763,32 @@ void core::set_use_nra_model(bool m) { void core::collect_statistics(::statistics & st) { st.update("arith-nla-explanations", m_stats.m_nla_explanations); st.update("arith-nla-lemmas", m_stats.m_nla_lemmas); - st.update("arith-nra-calls", m_stats.m_nra_calls); + st.update("arith-nra-calls", m_stats.m_nra_calls); + st.update("arith-bounds-improvements", m_stats.m_bounds_improvements); } +bool core::improve_bounds() { + return false; + + uint_set seen; + bool bounds_improved = false; + auto insert = [&](lpvar v) { + if (seen.contains(v)) + return; + seen.insert(v); + if (lra.improve_bound(v, false)) + bounds_improved = true, m_stats.m_bounds_improvements++; + if (lra.improve_bound(v, true)) + bounds_improved = true, m_stats.m_bounds_improvements++; + }; + for (auto & m : m_emons) { + insert(m.var()); + for (auto v : m.vars()) + insert(v); + } + return bounds_improved; +} + } // end of nla diff --git a/src/math/lp/nla_core.h b/src/math/lp/nla_core.h index 78f46bb41..0bd4a5814 100644 --- a/src/math/lp/nla_core.h +++ b/src/math/lp/nla_core.h @@ -64,6 +64,7 @@ class core { unsigned m_nla_explanations; unsigned m_nla_lemmas; unsigned m_nra_calls; + unsigned m_bounds_improvements; stats() { reset(); } void reset() { memset(this, 0, sizeof(*this)); @@ -79,7 +80,7 @@ class core { var_eqs m_evars; - lp::lar_solver& m_lar_solver; + lp::lar_solver& lra; reslimit& m_reslim; smt_params_helper m_params; std::function m_relevant; @@ -110,6 +111,9 @@ class core { void check_weighted(unsigned sz, std::pair>* checks); + // try to improve bounds for variables in monomials. + bool improve_bounds(); + public: // constructor core(lp::lar_solver& s, params_ref const& p, reslimit&); @@ -142,13 +146,13 @@ public: bool ineq_holds(const ineq& n) const; bool lemma_holds(const lemma& l) const; bool is_monic_var(lpvar j) const { return m_emons.is_monic_var(j); } - const rational& val(lpvar j) const { return m_lar_solver.get_column_value(j).x; } + const rational& val(lpvar j) const { return lra.get_column_value(j).x; } - const rational& var_val(const monic& m) const { return m_lar_solver.get_column_value(m.var()).x; } + const rational& var_val(const monic& m) const { return lra.get_column_value(m.var()).x; } rational mul_val(const monic& m) const { rational r(1); - for (lpvar v : m.vars()) r *= m_lar_solver.get_column_value(v).x; + for (lpvar v : m.vars()) r *= lra.get_column_value(v).x; return r; } @@ -287,16 +291,16 @@ public: } const rational& get_upper_bound(unsigned j) const; const rational& get_lower_bound(unsigned j) const; - bool has_lower_bound(lp::var_index var, lp::constraint_index& ci, lp::mpq& value, bool& is_strict) const { - return m_lar_solver.has_lower_bound(var, ci, value, is_strict); + bool has_lower_bound(lp::var_index var, u_dependency*& ci, lp::mpq& value, bool& is_strict) const { + return lra.has_lower_bound(var, ci, value, is_strict); } - bool has_upper_bound(lp::var_index var, lp::constraint_index& ci, lp::mpq& value, bool& is_strict) const { - return m_lar_solver.has_upper_bound(var, ci, value, is_strict); + bool has_upper_bound(lp::var_index var, u_dependency*& ci, lp::mpq& value, bool& is_strict) const { + return lra.has_upper_bound(var, ci, value, is_strict); } bool zero_is_an_inner_point_of_bounds(lpvar j) const; - bool var_is_int(lpvar j) const { return m_lar_solver.column_is_int(j); } + bool var_is_int(lpvar j) const { return lra.column_is_int(j); } int rat_sign(const monic& m) const; inline int rat_sign(lpvar j) const { return nla::rat_sign(val(j)); } @@ -338,7 +342,7 @@ public: bool is_octagon_term(const lp::lar_term& t, bool & sign, lpvar& i, lpvar &j) const; - void add_equivalence_maybe(const lp::lar_term *t, lpci c0, lpci c1); + void add_equivalence_maybe(const lp::lar_term* t, u_dependency* c0, u_dependency* c1); void init_vars_equivalence(); diff --git a/src/math/lp/nla_defs.h b/src/math/lp/nla_defs.h index df9158b42..e9807eaeb 100644 --- a/src/math/lp/nla_defs.h +++ b/src/math/lp/nla_defs.h @@ -14,11 +14,10 @@ #include "math/lp/explanation.h" namespace nla { -typedef lp::constraint_index lpci; -typedef lp::lconstraint_kind llc; -typedef lp::constraint_index lpci; -typedef lp::explanation expl_set; -typedef lp::var_index lpvar; + typedef lp::constraint_index lpci; + typedef lp::lconstraint_kind llc; + typedef lp::explanation expl_set; + typedef lp::var_index lpvar; struct from_index_dummy{}; class signed_var { diff --git a/src/math/lp/nla_grobner.cpp b/src/math/lp/nla_grobner.cpp index eefb07fac..f2fe8d9b2 100644 --- a/src/math/lp/nla_grobner.cpp +++ b/src/math/lp/nla_grobner.cpp @@ -22,9 +22,9 @@ namespace nla { grobner::grobner(core* c): common(c), - m_pdd_manager(m_core.m_lar_solver.number_of_vars()), - m_solver(m_core.m_reslim, m_pdd_manager), - m_lar_solver(m_core.m_lar_solver), + m_pdd_manager(m_core.lra.number_of_vars()), + m_solver(m_core.m_reslim, m_core.lra.dep_manager(), m_pdd_manager), + lra(m_core.lra), m_quota(m_core.params().arith_nl_gr_q()) {} @@ -211,12 +211,12 @@ namespace nla { TRACE("grobner", tout << "base vars: "; for (lpvar j : c().active_var_set()) - if (m_lar_solver.is_base(j)) + if (lra.is_base(j)) tout << "j" << j << " "; tout << "\n"); for (lpvar j : c().active_var_set()) { - if (m_lar_solver.is_base(j)) - add_row(m_lar_solver.basic2row(j)); + if (lra.is_base(j)) + add_row(lra.basic2row(j)); if (c().is_monic_var(j) && c().var_is_fixed(j)) add_fixed_monic(j); @@ -267,12 +267,12 @@ namespace nla { } } - for (unsigned j = 0; j < m_lar_solver.number_of_vars(); ++j) { - if (m_lar_solver.column_has_lower_bound(j) || m_lar_solver.column_has_upper_bound(j)) { + for (unsigned j = 0; j < lra.number_of_vars(); ++j) { + if (lra.column_has_lower_bound(j) || lra.column_has_upper_bound(j)) { out << j << ": ["; - if (m_lar_solver.column_has_lower_bound(j)) out << m_lar_solver.get_lower_bound(j); + if (lra.column_has_lower_bound(j)) out << lra.get_lower_bound(j); out << ".."; - if (m_lar_solver.column_has_upper_bound(j)) out << m_lar_solver.get_upper_bound(j); + if (lra.column_has_upper_bound(j)) out << lra.get_upper_bound(j); out << "]\n"; } } @@ -343,14 +343,14 @@ namespace nla { if (c().var_is_fixed(j)) return; - const auto& matrix = m_lar_solver.A_r(); + const auto& matrix = lra.A_r(); for (auto & s : matrix.m_columns[j]) { unsigned row = s.var(); if (m_rows.contains(row)) continue; m_rows.insert(row); - unsigned k = m_lar_solver.get_base_column_in_row(row); - if (m_lar_solver.column_is_free(k) && k != j) + unsigned k = lra.get_base_column_in_row(row); + if (lra.column_is_free(k) && k != j) continue; CTRACE("grobner", matrix.m_rows[row].size() > c().params().arith_nl_grobner_row_length_limit(), tout << "ignore the row " << row << " with the size " << matrix.m_rows[row].size() << "\n";); @@ -362,11 +362,10 @@ namespace nla { } const rational& grobner::val_of_fixed_var_with_deps(lpvar j, u_dependency*& dep) { - unsigned lc, uc; - m_lar_solver.get_bound_constraint_witnesses_for_column(j, lc, uc); - dep = c().m_intervals.mk_join(dep, c().m_intervals.mk_leaf(lc)); - dep = c().m_intervals.mk_join(dep, c().m_intervals.mk_leaf(uc)); - return m_lar_solver.column_lower_bound(j).x; + auto* d = lra.get_bound_constraint_witnesses_for_column(j); + if (d) + dep = c().m_intervals.mk_join(dep, d); + return lra.column_lower_bound(j).x; } dd::pdd grobner::pdd_expr(const rational& coeff, lpvar j, u_dependency*& dep) { @@ -415,7 +414,7 @@ namespace nla { SASSERT(r.hi().is_val()); v = r.var(); rational val = r.hi().val(); - switch (m_lar_solver.get_column_type(v)) { + switch (lra.get_column_type(v)) { case lp::column_type::lower_bound: if (val > 0) num_lo++, lo = v, lc = val; else num_hi++, hi = v, hc = val; break; @@ -508,7 +507,7 @@ namespace nla { void grobner::display_matrix_of_m_rows(std::ostream & out) const { - const auto& matrix = m_lar_solver.A_r(); + const auto& matrix = lra.A_r(); out << m_rows.size() << " rows" << "\n"; out << "the matrix\n"; for (const auto & r : matrix.m_rows) @@ -516,7 +515,7 @@ namespace nla { } void grobner::set_level2var() { - unsigned n = m_lar_solver.column_count(); + unsigned n = lra.column_count(); unsigned_vector sorted_vars(n), weighted_vars(n); for (unsigned j = 0; j < n; j++) { sorted_vars[j] = j; diff --git a/src/math/lp/nla_grobner.h b/src/math/lp/nla_grobner.h index f8e21d1a7..c7d41413c 100644 --- a/src/math/lp/nla_grobner.h +++ b/src/math/lp/nla_grobner.h @@ -21,7 +21,7 @@ namespace nla { class grobner : common { dd::pdd_manager m_pdd_manager; dd::solver m_solver; - lp::lar_solver& m_lar_solver; + lp::lar_solver& lra; indexed_uint_set m_rows; unsigned m_quota = 0; diff --git a/src/math/lp/nla_intervals.cpp b/src/math/lp/nla_intervals.cpp index 4ffbcb7e3..791251fee 100644 --- a/src/math/lp/nla_intervals.cpp +++ b/src/math/lp/nla_intervals.cpp @@ -4,6 +4,11 @@ #include "util/mpq.h" namespace nla { + +intervals::intervals(core* c, reslimit& lim): + m_dep_intervals(c->lra.dep_manager(), lim), + m_core(c) {} + typedef enum dep_intervals::with_deps_t e_with_deps; const nex* intervals::get_inf_interval_child(const nex_sum& e) const { @@ -173,7 +178,7 @@ lp::lar_term intervals::expression_to_normalized_term(const nex_sum* e, rational // where m_terms[k] corresponds to the returned lpvar lpvar intervals::find_term_column(const lp::lar_term & norm_t, rational& a) const { std::pair a_j; - if (m_core->m_lar_solver.fetch_normalized_term_column(norm_t, a_j)) { + if (m_core->lra.fetch_normalized_term_column(norm_t, a_j)) { a /= a_j.first; return a_j.second; } @@ -206,19 +211,10 @@ void intervals::set_zero_interval_deps_for_mult(interval& a) { a.m_upper_dep = a.m_lower_dep; } -u_dependency *intervals::mk_dep(lp::constraint_index ci) { - return m_dep_intervals.mk_leaf(ci); -} - -u_dependency *intervals::mk_dep(const lp::explanation& expl) { +u_dependency* intervals::mk_dep(const lp::explanation& expl) { u_dependency * r = nullptr; - for (auto p : expl) { - if (r == nullptr) { - r = m_dep_intervals.mk_leaf(p.ci()); - } else { - r = m_dep_intervals.mk_join(r, m_dep_intervals.mk_leaf(p.ci())); - } - } + for (auto p : expl) + r = m_dep_intervals.mk_join(r, m_dep_intervals.mk_leaf(p.ci())); return r; } @@ -249,25 +245,25 @@ std::ostream& intervals::display(std::ostream& out, const interval& i) const { template void intervals::set_var_interval(lpvar v, interval& b) { TRACE("nla_intervals_details", m_core->print_var(v, tout) << "\n";); - lp::constraint_index ci; + u_dependency* dep = nullptr; rational val; bool is_strict; - if (ls().has_lower_bound(v, ci, val, is_strict)) { + if (ls().has_lower_bound(v, dep, val, is_strict)) { m_dep_intervals.set_lower(b, val); m_dep_intervals.set_lower_is_open(b, is_strict); m_dep_intervals.set_lower_is_inf(b, false); - if (wd == e_with_deps::with_deps) b.m_lower_dep = mk_dep(ci); + if (wd == e_with_deps::with_deps) b.m_lower_dep = dep; } else { m_dep_intervals.set_lower_is_open(b, true); m_dep_intervals.set_lower_is_inf(b, true); if (wd == e_with_deps::with_deps) b.m_lower_dep = nullptr; } - if (ls().has_upper_bound(v, ci, val, is_strict)) { + if (ls().has_upper_bound(v, dep, val, is_strict)) { m_dep_intervals.set_upper(b, val); m_dep_intervals.set_upper_is_open(b, is_strict); m_dep_intervals.set_upper_is_inf(b, false); - if (wd == e_with_deps::with_deps) b.m_upper_dep = mk_dep(ci); + if (wd == e_with_deps::with_deps) b.m_upper_dep = dep; } else { m_dep_intervals.set_upper_is_open(b, true); @@ -303,7 +299,7 @@ bool intervals::interval_from_term(const nex& e, scoped_dep_interval& i) { m_dep_intervals.set(i, bi); TRACE("nla_intervals", - m_core->m_lar_solver.print_column_info(j, tout) << "\n"; + m_core->lra.print_column_info(j, tout) << "\n"; tout << "a=" << a << ", b=" << b << "\n"; tout << e << ", interval = "; display(tout, i);); return true; @@ -476,9 +472,9 @@ bool intervals::interval_of_expr(const nex* e, unsigned p, scoped_dep_interval& } -lp::lar_solver& intervals::ls() { return m_core->m_lar_solver; } +lp::lar_solver& intervals::ls() { return m_core->lra; } -const lp::lar_solver& intervals::ls() const { return m_core->m_lar_solver; } +const lp::lar_solver& intervals::ls() const { return m_core->lra; } } // end of nla namespace diff --git a/src/math/lp/nla_intervals.h b/src/math/lp/nla_intervals.h index 0545e9933..514049aca 100644 --- a/src/math/lp/nla_intervals.h +++ b/src/math/lp/nla_intervals.h @@ -25,16 +25,13 @@ class intervals { public: typedef dep_intervals::interval interval; private: - u_dependency* mk_dep(lp::constraint_index ci); u_dependency* mk_dep(lp::explanation const&); lp::lar_solver& ls(); const lp::lar_solver& ls() const; public: - intervals(core* c, reslimit& lim) : - m_dep_intervals(lim), - m_core(c) - {} + intervals(core* c, reslimit& lim); + dep_intervals& get_dep_intervals() { return m_dep_intervals; } u_dependency* mk_join(u_dependency* a, u_dependency* b) { return m_dep_intervals.mk_join(a, b); } u_dependency* mk_leaf(lp::constraint_index ci) { return m_dep_intervals.mk_leaf(ci); } diff --git a/src/math/lp/nra_solver.cpp b/src/math/lp/nra_solver.cpp index 377b318dd..dc4681559 100644 --- a/src/math/lp/nra_solver.cpp +++ b/src/math/lp/nra_solver.cpp @@ -186,7 +186,7 @@ struct solver::imp { for (auto const& eq : eqs) add_eq(eq); for (auto const& [v, w] : m_lp2nl) { - auto& ls = m_nla_core.m_lar_solver; + auto& ls = m_nla_core.lra; if (ls.column_has_lower_bound(v)) add_lb(ls.get_lower_bound(v), w); if (ls.column_has_upper_bound(v)) @@ -209,7 +209,7 @@ struct solver::imp { IF_VERBOSE(0, verbose_stream() << "check-nra " << r << "\n"; m_nlsat->display(verbose_stream()); for (auto const& [v, w] : m_lp2nl) { - auto& ls = m_nla_core.m_lar_solver; + auto& ls = m_nla_core.lra; if (ls.column_has_lower_bound(v)) verbose_stream() << w << " >= " << ls.get_lower_bound(v) << "\n"; if (ls.column_has_upper_bound(v)) diff --git a/src/math/lp/ul_pair.h b/src/math/lp/ul_pair.h index abfb4483b..37fd6e9ed 100644 --- a/src/math/lp/ul_pair.h +++ b/src/math/lp/ul_pair.h @@ -1,25 +1,19 @@ /*++ Copyright (c) 2017 Microsoft Corporation -Module Name: - - - Abstract: - + justifications for upper or lower bounds Author: Lev Nachmanson (levnach) -Revision History: - - --*/ #pragma once #include "util/vector.h" +#include "util/dependency.h" #include #include #include @@ -48,14 +42,20 @@ inline bool compare(const std::pair & a, const std::pair cs) { + eq_justification(std::initializer_list cs) { int i = 0; - for (lpci c: cs) { + for (auto c: cs) { m_cs[i++] = c; } for (; i < 4; i++) { - m_cs[i] = -1; + m_cs[i] = nullptr; } } - void explain(lp::explanation& e) const { - for (lpci c : m_cs) - if (c + 1 != 0) // c != -1 - e.push_back(c); + u_dependency* const* begin() const { return m_cs; } + u_dependency* const* end() const { + unsigned i = 0; + for (; i < 4 && m_cs[i]; ++i) + ; + return m_cs + i; } }; @@ -202,7 +204,7 @@ public: } for (eq_justification const& j : m_justtrail) { - j.explain(e); + explain_eq(j, e); } m_stats.m_num_explains += m_justtrail.size(); m_stats.m_num_explain_calls++; @@ -216,6 +218,17 @@ public: // IF_VERBOSE(2, verbose_stream() << (double)m_stats.m_num_explains / m_stats.m_num_explain_calls << "\n"); } + void explain_eq(eq_justification const& eq, lp::explanation& e) const { + u_dependency_manager dm; + unsigned_vector deps; + for (auto* dep : eq) { + deps.reset(); + dm.linearize(dep, deps); + for (auto ci : deps) + e.push_back(ci); + } + } + void explain_bfs(signed_var v1, signed_var v2, lp::explanation& e) const { SASSERT(find(v1) == find(v2)); if (v1 == v2) { @@ -249,7 +262,7 @@ public: } while (head != 0) { - m_justtrail[head].explain(e); + explain_eq(m_justtrail[head], e); head = m_todo[head].m_index; ++m_stats.m_num_explains; } diff --git a/src/sat/sat_anf_simplifier.cpp b/src/sat/sat_anf_simplifier.cpp index 6d28d9f98..5ed45171c 100644 --- a/src/sat/sat_anf_simplifier.cpp +++ b/src/sat/sat_anf_simplifier.cpp @@ -42,7 +42,8 @@ namespace sat { void anf_simplifier::operator()() { dd::pdd_manager m(20, dd::pdd_manager::semantics::mod2_e); - pdd_solver solver(s.rlimit(), m); + u_dependency_manager dm; + pdd_solver solver(s.rlimit(), dm, m); report _report(*this); configure_solver(solver); clauses2anf(solver); diff --git a/src/sat/smt/arith_axioms.cpp b/src/sat/smt/arith_axioms.cpp index c2a96177d..d17154361 100644 --- a/src/sat/smt/arith_axioms.cpp +++ b/src/sat/smt/arith_axioms.cpp @@ -524,7 +524,7 @@ namespace arith { return all_divs_valid; } - void solver::fixed_var_eh(theory_var v, lp::constraint_index ci1, lp::constraint_index ci2, rational const& bound) { + void solver::fixed_var_eh(theory_var v, u_dependency* dep, rational const& bound) { theory_var w = euf::null_theory_var; enode* x = var2enode(v); if (bound.is_zero()) @@ -540,8 +540,8 @@ namespace arith { return; reset_evidence(); m_explanation.clear(); - consume(rational::one(), ci1); - consume(rational::one(), ci2); + for (auto ci : lp().flatten(dep)) + consume(rational::one(), ci); ++m_stats.m_fixed_eqs; auto* hint = explain_implied_eq(m_explanation, x, y); auto* jst = euf::th_explain::propagate(*this, m_core, m_eqs, x, y, hint); diff --git a/src/sat/smt/arith_solver.cpp b/src/sat/smt/arith_solver.cpp index 0e97c3503..bda3bebcd 100644 --- a/src/sat/smt/arith_solver.cpp +++ b/src/sat/smt/arith_solver.cpp @@ -402,12 +402,13 @@ namespace arith { } void solver::propagate_eqs(lp::tv t, lp::constraint_index ci1, lp::lconstraint_kind k, api_bound& b, rational const& value) { - lp::constraint_index ci2; - if (k == lp::GE && set_lower_bound(t, ci1, value) && has_upper_bound(t.index(), ci2, value)) { - fixed_var_eh(b.get_var(), ci1, ci2, value); + u_dependency* dep; + auto& dm = lp().dep_manager(); + if (k == lp::GE && set_lower_bound(t, ci1, value) && has_upper_bound(t.index(), dep, value)) { + fixed_var_eh(b.get_var(), dm.mk_join(dm.mk_leaf(ci1), dep), value); } - else if (k == lp::LE && set_upper_bound(t, ci1, value) && has_lower_bound(t.index(), ci2, value)) { - fixed_var_eh(b.get_var(), ci1, ci2, value); + else if (k == lp::LE && set_upper_bound(t, ci1, value) && has_lower_bound(t.index(), dep, value)) { + fixed_var_eh(b.get_var(), dm.mk_join(dm.mk_leaf(ci1), dep), value); } } @@ -433,11 +434,12 @@ namespace arith { // m_solver already tracks bounds on proper variables, but not on terms. bool is_strict = false; rational b; + u_dependency* dep = nullptr; if (is_lower) { - return lp().has_lower_bound(tv.id(), ci, b, is_strict) && !is_strict && b == v; + return lp().has_lower_bound(tv.id(), dep, b, is_strict) && !is_strict && b == v; } else { - return lp().has_upper_bound(tv.id(), ci, b, is_strict) && !is_strict && b == v; + return lp().has_upper_bound(tv.id(), dep, b, is_strict) && !is_strict && b == v; } } } @@ -692,7 +694,7 @@ namespace arith { void solver::report_equality_of_fixed_vars(unsigned vi1, unsigned vi2) { rational bound; - lp::constraint_index ci1, ci2, ci3, ci4; + u_dependency* ci1 = nullptr, *ci2 = nullptr, *ci3 = nullptr, *ci4 = nullptr; theory_var v1 = lp().local_to_external(vi1); theory_var v2 = lp().local_to_external(vi2); TRACE("arith", tout << "fixed: " << mk_pp(var2expr(v1), m) << " " << mk_pp(var2expr(v2), m) << "\n";); @@ -714,10 +716,10 @@ namespace arith { ++m_stats.m_fixed_eqs; reset_evidence(); m_explanation.clear(); - consume(rational::one(), ci1); - consume(rational::one(), ci2); - consume(rational::one(), ci3); - consume(rational::one(), ci4); + auto& dm = lp().dep_manager(); + auto* d = dm.mk_join(dm.mk_join(ci1, ci2), dm.mk_join(ci3, ci4)); + for (auto ci : lp().flatten(d)) + consume(rational::one(), ci); enode* x = var2enode(v1); enode* y = var2enode(v2); auto* ex = explain_implied_eq(m_explanation, x, y); @@ -729,26 +731,28 @@ namespace arith { return x == y || var2enode(x)->get_root() == var2enode(y)->get_root(); } - bool solver::has_upper_bound(lpvar vi, lp::constraint_index& ci, rational const& bound) { return has_bound(vi, ci, bound, false); } + bool solver::has_upper_bound(lpvar vi, u_dependency*& ci, rational const& bound) { return has_bound(vi, ci, bound, false); } - bool solver::has_lower_bound(lpvar vi, lp::constraint_index& ci, rational const& bound) { return has_bound(vi, ci, bound, true); } + bool solver::has_lower_bound(lpvar vi, u_dependency*& ci, rational const& bound) { return has_bound(vi, ci, bound, true); } - bool solver::has_bound(lpvar vi, lp::constraint_index& ci, rational const& bound, bool is_lower) { + bool solver::has_bound(lpvar vi, u_dependency*& dep, rational const& bound, bool is_lower) { if (lp::tv::is_term(vi)) { theory_var v = lp().local_to_external(vi); rational val; TRACE("arith", tout << lp().get_variable_name(vi) << " " << v << "\n";); if (v != euf::null_theory_var && a.is_numeral(var2expr(v), val) && bound == val) { - ci = UINT_MAX; + dep = nullptr; return bound == val; } auto& vec = is_lower ? m_lower_terms : m_upper_terms; lpvar ti = lp::tv::unmask_term(vi); if (vec.size() > ti) { - constraint_bound& b = vec[ti]; - ci = b.first; - return ci != UINT_MAX && bound == b.second; + auto& [ci, coeff] = vec[ti]; + if (ci == UINT_MAX) + return false; + dep = lp().dep_manager().mk_leaf(ci); + return bound == coeff; } else { return false; @@ -758,10 +762,10 @@ namespace arith { bool is_strict = false; rational b; if (is_lower) { - return lp().has_lower_bound(vi, ci, b, is_strict) && b == bound && !is_strict; + return lp().has_lower_bound(vi, dep, b, is_strict) && b == bound && !is_strict; } else { - return lp().has_upper_bound(vi, ci, b, is_strict) && b == bound && !is_strict; + return lp().has_upper_bound(vi, dep, b, is_strict) && b == bound && !is_strict; } } } diff --git a/src/sat/smt/arith_solver.h b/src/sat/smt/arith_solver.h index 20252ede9..1b6f58782 100644 --- a/src/sat/smt/arith_solver.h +++ b/src/sat/smt/arith_solver.h @@ -359,7 +359,7 @@ namespace arith { void assert_idiv_mod_axioms(theory_var u, theory_var v, theory_var w, rational const& r); api_bound* mk_var_bound(sat::literal lit, theory_var v, lp_api::bound_kind bk, rational const& bound); lp::lconstraint_kind bound2constraint_kind(bool is_int, lp_api::bound_kind bk, bool is_true); - void fixed_var_eh(theory_var v1, lp::constraint_index ci1, lp::constraint_index ci2, rational const& bound); + void fixed_var_eh(theory_var v1, u_dependency* dep, rational const& bound); bool set_upper_bound(lp::tv t, lp::constraint_index ci, rational const& v) { return set_bound(t, ci, v, false); } bool set_lower_bound(lp::tv t, lp::constraint_index ci, rational const& v) { return set_bound(t, ci, v, true); } bool set_bound(lp::tv tv, lp::constraint_index ci, rational const& v, bool is_lower); @@ -412,9 +412,9 @@ namespace arith { nlsat::anum const& nl_value(theory_var v, scoped_anum& r) const; - bool has_bound(lpvar vi, lp::constraint_index& ci, rational const& bound, bool is_lower); - bool has_lower_bound(lpvar vi, lp::constraint_index& ci, rational const& bound); - bool has_upper_bound(lpvar vi, lp::constraint_index& ci, rational const& bound); + bool has_bound(lpvar vi, u_dependency*& ci, rational const& bound, bool is_lower); + bool has_lower_bound(lpvar vi, u_dependency*& ci, rational const& bound); + bool has_upper_bound(lpvar vi, u_dependency*& ci, rational const& bound); /* * Facility to put a small box around integer variables used in branch and bounds. diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp index 6dc14fb13..08fa39864 100644 --- a/src/smt/theory_lra.cpp +++ b/src/smt/theory_lra.cpp @@ -2863,7 +2863,7 @@ public: lp::lar_term const& term = lp().get_term(ti); for (auto const mono : term) { auto wi = lp().column2tv(mono.column()); - lp::constraint_index ci; + u_dependency* ci = nullptr; rational value; bool is_strict; if (wi.is_term()) { @@ -2966,12 +2966,13 @@ public: vector m_upper_terms; void propagate_eqs(lp::tv t, lp::constraint_index ci1, lp::lconstraint_kind k, api_bound& b, rational const& value) { - lp::constraint_index ci2; + u_dependency* ci2 = nullptr; + auto pair = [&]() { return lp().dep_manager().mk_join(lp().dep_manager().mk_leaf(ci1), ci2); }; if (k == lp::GE && set_lower_bound(t, ci1, value) && has_upper_bound(t.index(), ci2, value)) { - fixed_var_eh(b.get_var(), t, ci1, ci2, value); + fixed_var_eh(b.get_var(), t, pair(), value); } else if (k == lp::LE && set_upper_bound(t, ci1, value) && has_lower_bound(t.index(), ci2, value)) { - fixed_var_eh(b.get_var(), t, ci1, ci2, value); + fixed_var_eh(b.get_var(), t, pair(), value); } } @@ -3012,11 +3013,12 @@ public: // m_solver already tracks bounds on proper variables, but not on terms. bool is_strict = false; rational b; + u_dependency* dep = nullptr; if (is_lower) { - return lp().has_lower_bound(tv.id(), ci, b, is_strict) && !is_strict && b == v; + return lp().has_lower_bound(tv.id(), dep, b, is_strict) && !is_strict && b == v; } else { - return lp().has_upper_bound(tv.id(), ci, b, is_strict) && !is_strict && b == v; + return lp().has_upper_bound(tv.id(), dep, b, is_strict) && !is_strict && b == v; } } } @@ -3024,35 +3026,37 @@ public: bool var_has_bound(lpvar vi, bool is_lower) { bool is_strict = false; rational b; - lp::constraint_index ci; + u_dependency* dep; if (is_lower) { - return lp().has_lower_bound(vi, ci, b, is_strict); + return lp().has_lower_bound(vi, dep, b, is_strict); } else { - return lp().has_upper_bound(vi, ci, b, is_strict); + return lp().has_upper_bound(vi, dep, b, is_strict); } } - bool has_upper_bound(lpvar vi, lp::constraint_index& ci, rational const& bound) { return has_bound(vi, ci, bound, false); } + bool has_upper_bound(lpvar vi, u_dependency*& ci, rational const& bound) { return has_bound(vi, ci, bound, false); } - bool has_lower_bound(lpvar vi, lp::constraint_index& ci, rational const& bound) { return has_bound(vi, ci, bound, true); } + bool has_lower_bound(lpvar vi, u_dependency*& ci, rational const& bound) { return has_bound(vi, ci, bound, true); } - bool has_bound(lpvar vi, lp::constraint_index& ci, rational const& bound, bool is_lower) { + bool has_bound(lpvar vi, u_dependency*& dep, rational const& bound, bool is_lower) { if (lp::tv::is_term(vi)) { theory_var v = lp().local_to_external(vi); rational val; TRACE("arith", tout << lp().get_variable_name(vi) << " " << v << "\n";); if (v != null_theory_var && a.is_numeral(get_owner(v), val) && bound == val) { - ci = UINT_MAX; + dep = nullptr; return bound == val; } auto& vec = is_lower ? m_lower_terms : m_upper_terms; lpvar ti = lp::tv::unmask_term(vi); if (vec.size() > ti) { - constraint_bound& b = vec[ti]; - ci = b.first; - return ci != UINT_MAX && bound == b.second; + auto const& [ci, coeff] = vec[ti]; + if (ci == UINT_MAX) + return false; + dep = lp().dep_manager().mk_leaf(ci); + return bound == coeff; } else { return false; @@ -3062,10 +3066,10 @@ public: bool is_strict = false; rational b; if (is_lower) { - return lp().has_lower_bound(vi, ci, b, is_strict) && b == bound && !is_strict; + return lp().has_lower_bound(vi, dep, b, is_strict) && b == bound && !is_strict; } else { - return lp().has_upper_bound(vi, ci, b, is_strict) && b == bound && !is_strict; + return lp().has_upper_bound(vi, dep, b, is_strict) && b == bound && !is_strict; } } } @@ -3078,7 +3082,7 @@ public: void report_equality_of_fixed_vars(unsigned vi1, unsigned vi2) { rational bound(0); - lp::constraint_index ci1, ci2, ci3, ci4; + u_dependency* ci1 = nullptr, *ci2 = nullptr, *ci3 = nullptr, *ci4 = nullptr; theory_var v1 = lp().local_to_external(vi1); theory_var v2 = lp().local_to_external(vi2); TRACE("arith", tout << "fixed: " << mk_pp(get_owner(v1), m) << " " << mk_pp(get_owner(v2), m) << "\n";); @@ -3130,7 +3134,7 @@ public: ctx().assign_eq(x, y, eq_justification(js)); } - void fixed_var_eh(theory_var v, lp::tv t, lp::constraint_index ci1, lp::constraint_index ci2, rational const& bound) { + void fixed_var_eh(theory_var v, lp::tv t, u_dependency* dep, rational const& bound) { theory_var w = null_theory_var; enode* x = get_enode(v); if (m_value2var.find(bound, w)) @@ -3147,8 +3151,7 @@ public: if (x->get_root() == y->get_root()) return; reset_evidence(); - set_evidence(ci1, m_core, m_eqs); - set_evidence(ci2, m_core, m_eqs); + set_evidence(dep, m_core, m_eqs); ++m_stats.m_fixed_eqs; assign_eq(v, w); } @@ -3182,6 +3185,11 @@ public: // lp::constraint_index const null_constraint_index = UINT_MAX; // not sure what a correct fix is + void set_evidence(u_dependency* dep, literal_vector& core, svector& eqs) { + for (auto ci : lp().flatten(dep)) + set_evidence(ci, core, eqs); + } + void set_evidence(lp::constraint_index idx, literal_vector& core, svector& eqs) { if (idx == UINT_MAX) { return; @@ -3402,7 +3410,7 @@ public: if (!is_registered_var(v)) return false; lpvar vi = get_lpvar(v); - lp::constraint_index ci; + u_dependency* ci; return lp().has_lower_bound(vi, ci, val, is_strict); } @@ -3421,8 +3429,8 @@ public: if (!is_registered_var(v)) return false; lpvar vi = get_lpvar(v); - lp::constraint_index ci; - return lp().has_upper_bound(vi, ci, val, is_strict); + u_dependency* dep = nullptr; + return lp().has_upper_bound(vi, dep, val, is_strict); } diff --git a/src/test/pdd_solver.cpp b/src/test/pdd_solver.cpp index a86504266..33b0abf77 100644 --- a/src/test/pdd_solver.cpp +++ b/src/test/pdd_solver.cpp @@ -20,13 +20,15 @@ namespace dd { } void test1() { pdd_manager m(4); + u_dependency_manager dm; reslimit lim; pdd v0 = m.mk_var(0); pdd v1 = m.mk_var(1); pdd v2 = m.mk_var(2); pdd v3 = m.mk_var(3); - solver gb(lim, m); + + solver gb(lim, dm, m); gb.add(v1*v2 + v1*v3); gb.add(v1 - 1); gb.display(std::cout); @@ -198,10 +200,11 @@ namespace dd { void test_simplify(expr_ref_vector& fmls, bool use_mod2) { ast_manager& m = fmls.get_manager(); unsigned_vector id2var; + u_dependency_manager dm; collect_id2var(id2var, fmls); pdd_manager p(id2var.size(), use_mod2 ? pdd_manager::mod2_e : pdd_manager::zero_one_vars_e); - solver g(m.limit(), p); + solver g(m.limit(), dm, p); for (expr* e : subterms::ground(fmls)) { add_def(id2var, to_app(e), m, p, g); diff --git a/src/util/dependency.h b/src/util/dependency.h index 7ccba716a..57057460c 100644 --- a/src/util/dependency.h +++ b/src/util/dependency.h @@ -69,7 +69,7 @@ private: value_manager & m_vmanager; allocator & m_allocator; - ptr_vector m_todo; + mutable ptr_vector m_todo; void inc_ref(value const & v) { if (C::ref_count) @@ -106,12 +106,9 @@ private: } } - void unmark_todo() { - typename ptr_vector::iterator it = m_todo.begin(); - typename ptr_vector::iterator end = m_todo.end(); - for (; it != end; ++it) { - (*it)->unmark(); - } + void unmark_todo() const { + for (auto* d : m_todo) + d->unmark(); m_todo.reset(); } @@ -193,7 +190,7 @@ public: return false; } - void linearize(dependency * d, vector & vs) { + void linearize(dependency * d, vector & vs) const { if (d) { m_todo.reset(); d->mark(); @@ -300,7 +297,7 @@ public: return m_dep_manager.contains(d, v); } - void linearize(dependency * d, vector & vs) { + void linearize(dependency * d, vector & vs) const { return m_dep_manager.linearize(d, vs); }