From bd8e5eee4b78587bf45ea9c12691914869fb4378 Mon Sep 17 00:00:00 2001 From: Nikolaj Bjorner Date: Sun, 29 Oct 2023 10:21:24 -0700 Subject: [PATCH] add simplification experiment (disabled) for tracking, some reshuffling of equation/fixed_equation structs Signed-off-by: Nikolaj Bjorner --- src/math/lp/explanation.h | 14 ++ src/math/lp/int_solver.cpp | 231 ++++++++++++++++++++++++++++++ src/math/lp/int_solver.h | 6 +- src/math/lp/lp_types.h | 1 + src/math/lp/monomial_bounds.cpp | 4 +- src/math/lp/nla_core.cpp | 13 +- src/math/lp/nla_core.h | 9 +- src/math/lp/nla_solver.cpp | 8 +- src/math/lp/nla_solver.h | 5 +- src/math/lp/nla_types.h | 12 -- src/sat/smt/arith_diagnostics.cpp | 1 - src/smt/theory_lra.cpp | 21 ++- 12 files changed, 286 insertions(+), 39 deletions(-) diff --git a/src/math/lp/explanation.h b/src/math/lp/explanation.h index b7f721a30..960c5fb4a 100644 --- a/src/math/lp/explanation.h +++ b/src/math/lp/explanation.h @@ -121,4 +121,18 @@ public: } }; + + struct equality { + lp::lpvar i, j; + lp::explanation e; + equality(lp::lpvar i, lp::lpvar j, lp::explanation const& e):i(i),j(j),e(e) {} + }; + + struct fixed_equality { + lp::lpvar v; + rational k; + lp::explanation e; + fixed_equality(lp::lpvar v, rational const& k, lp::explanation const& e):v(v),k(k),e(e) {} + }; + } diff --git a/src/math/lp/int_solver.cpp b/src/math/lp/int_solver.cpp index e248a0081..35a082da9 100644 --- a/src/math/lp/int_solver.cpp +++ b/src/math/lp/int_solver.cpp @@ -683,5 +683,236 @@ namespace lp { return result; } + void int_solver::simplify(std::function& is_root) { + +#if 0 + + // in-processing simplification can go here, such as bounds improvements. + + if (!lra.is_feasible()) { + lra.find_feasible_solution(); + if (!lra.is_feasible()) + return; + } + + stopwatch sw; + explanation exp1, exp2; + + // + // tighten integer bounds + // It is a weak method as it ownly strengthens bounds if + // variables are already at one of the to-be discovered bounds. + // + sw.start(); + unsigned changes = 0; + auto const& constraints = lra.constraints(); + auto print_var = [this](lpvar j) { + if (lra.column_corresponds_to_term(j)) { + std::stringstream strm; + lra.print_column_info(j, strm); + return strm.str(); + } + else + return std::string("j") + std::to_string(j); + }; + unsigned start = random(); + unsigned num_checks = 0; + for (lpvar j0 = 0; j0 < lra.column_count(); ++j0) { + lpvar j = (j0 + start) % lra.column_count(); + + if (num_checks > 1000) + break; + if (is_fixed(j)) + continue; + if (!lra.column_is_int(j)) + continue; + rational value = get_value(j).x; + bool tight_lower = false, tight_upper = false; + u_dependency* dep; + + if (!value.is_int()) + continue; + + bool at_up = at_upper(j); + + if (!at_lower(j)) { + ++num_checks; + lra.push(); + auto k = lp::lconstraint_kind::LE; + lra.update_column_type_and_bound(j, k, (value - 1).to_mpq(), nullptr); + lra.find_feasible_solution(); + if (!lra.is_feasible()) { + tight_upper = true; + ++changes; + lra.get_infeasibility_explanation(exp1); +#if 0 + display_column(std::cout, j); + std::cout << print_var(j) << " >= " << value << "\n"; + unsigned i = 0; + for (auto p : exp1) { + std::cout << "(" << p.ci() << ")"; + constraints.display(std::cout, print_var, p.ci()); + if (++i < exp1.size()) + std::cout << " "; + } +#endif + } + lra.pop(1); + if (tight_upper) { + dep = nullptr; + for (auto& cc : exp1) + dep = lra.join_deps(dep, constraints[cc.ci()].dep()); + lra.update_column_type_and_bound(j, lp::lconstraint_kind::GE, value.to_mpq(), dep); + } + } + + if (!at_up) { + ++num_checks; + lra.push(); + auto k = lp::lconstraint_kind::GE; + lra.update_column_type_and_bound(j, k, (value + 1).to_mpq(), nullptr); + lra.find_feasible_solution(); + if (!lra.is_feasible()) { + tight_lower = true; + ++changes; + lra.get_infeasibility_explanation(exp1); +#if 0 + display_column(std::cout, j); + std::cout << print_var(j) << " <= " << value << "\n"; + unsigned i = 0; + for (auto p : exp1) { + std::cout << "(" << p.ci() << ")"; + constraints.display(std::cout, print_var, p.ci()); + if (++i < exp1.size()) + std::cout << " "; + } +#endif + } + lra.pop(1); + if (tight_lower) { + dep = nullptr; + for (auto& cc : exp1) + dep = lra.join_deps(dep, constraints[cc.ci()].dep()); + lra.update_column_type_and_bound(j, lp::lconstraint_kind::LE, value.to_mpq(), dep); + } + } + } + sw.stop(); + std::cout << "changes " << changes << " columns " << lra.column_count() << " time: " << sw.get_seconds() << "\n"; + std::cout.flush(); + + // + // identify equalities + // + + m_equalities.reset(); + map value2roots; + + vector> coeffs; + coeffs.push_back({-rational::one(), 0}); + coeffs.push_back({rational::one(), 0}); + + num_checks = 0; + + // make sure values are sampled with respect to the same state of the Simplex. + vector values; + for (lpvar j = 0; j < lra.column_count(); ++j) + values.push_back(get_value(j).x); + + sw.reset(); + sw.start(); + start = random(); + for (lpvar j0 = 0; j0 < lra.column_count(); ++j0) { + lpvar j = (j0 + start) % lra.column_count(); + if (is_fixed(j)) + continue; + if (!lra.column_is_int(j)) + continue; + if (!is_root(j)) + continue; + rational value = values[j]; + if (!value2roots.contains(value)) { + unsigned_vector vec; + vec.push_back(j); + value2roots.insert(value, vec); + continue; + } + auto& roots = value2roots.find(value); + bool has_eq = false; + // + // Super inefficient check. There are better ways. + // 1. call into equality finder: + // the cheap equality finder can also be used. + // 2. value sweeping: + // update partitions of values based on feasible tableaus + // instead of having just the values vector use the values + // collected when the find_feasible_solution succeeds with + // a new assignment. + // 3. a more expensive equality finder: + // use the tableau to extract equalities from tight rows. + // If x = y is implied, there is a set of rows that link x and y + // and such that the variables are at their bounds. + // 4. retain information between calls: + // If simplification is invoked at the same backtracking level (or above) + // form the previous call and it is established that x <= y (but not x == y), then no need to + // recheck the inequality x <= y. + for (auto k : roots) { + bool le = false, ge = false; + u_dependency* dep = nullptr; + lra.push(); + coeffs[0].second = j; + coeffs[1].second = k; + lp::lpvar term_index = lra.add_term(coeffs, UINT_MAX); + term_index = lra.map_term_index_to_column_index(term_index); + lra.push(); + lra.update_column_type_and_bound(term_index, lp::lconstraint_kind::GE, mpq(1), nullptr); + lra.find_feasible_solution(); + if (!lra.is_feasible()) { + lra.get_infeasibility_explanation(exp1); + le = true; + } + lra.pop(1); + ++num_checks; + if (le) { + lra.push(); + lra.update_column_type_and_bound(term_index, lp::lconstraint_kind::LE, mpq(-1), nullptr); + lra.find_feasible_solution(); + if (!lra.is_feasible()) { + lra.get_infeasibility_explanation(exp2); + exp1.add_expl(exp2); + ge = true; + } + lra.pop(1); + ++num_checks; + } + lra.pop(1); + if (le && ge) { + has_eq = true; + m_equalities.push_back({j, k, exp1}); + break; + } + // artificial throttle. + if (num_checks > 10000) + break; + } + if (!has_eq) + roots.push_back(j); + + // artificial throttle. + if (num_checks > 10000) + break; + } + + sw.stop(); + std::cout << "equalities " << m_equalities.size() << " num checks " << num_checks << " time: " << sw.get_seconds() << "\n"; + std::cout.flush(); + + // + // Cuts? Eg. for 0-1 variables or bounded integers? + // + +#endif + } + } diff --git a/src/math/lp/int_solver.h b/src/math/lp/int_solver.h index 1bd3f1c3d..f1faedf35 100644 --- a/src/math/lp/int_solver.h +++ b/src/math/lp/int_solver.h @@ -67,6 +67,8 @@ class int_solver { bool m_upper; // we have a cut m_t*x <= k if m_upper is true nad m_t*x >= k otherwise hnf_cutter m_hnf_cutter; unsigned m_hnf_cut_period; + + vector m_equalities; public: int_solver(lar_solver& lp); @@ -84,7 +86,9 @@ public: const impq & get_value(unsigned j) const; bool at_lower(unsigned j) const; bool at_upper(unsigned j) const; - + void simplify(std::function& is_root); + vector const& equalities() const { return m_equalities; } + private: // lia_move patch_nbasic_columns(); bool get_freedom_interval_for_column(unsigned j, bool & inf_l, impq & l, bool & inf_u, impq & u, mpq & m); diff --git a/src/math/lp/lp_types.h b/src/math/lp/lp_types.h index c69dbe10b..f770fe956 100644 --- a/src/math/lp/lp_types.h +++ b/src/math/lp/lp_types.h @@ -92,6 +92,7 @@ public: }; + } inline std::ostream& operator<<(std::ostream& out, lp::tv const& t) { diff --git a/src/math/lp/monomial_bounds.cpp b/src/math/lp/monomial_bounds.cpp index e54f82682..9c99fa8c3 100644 --- a/src/math/lp/monomial_bounds.cpp +++ b/src/math/lp/monomial_bounds.cpp @@ -391,8 +391,8 @@ namespace nla { void monomial_bounds::propagate_nonfixed(monic const& m, rational const& k, lpvar w) { vector> coeffs; - coeffs.push_back(std::make_pair(-k, w)); - coeffs.push_back(std::make_pair(rational::one(), m.var())); + coeffs.push_back({-k, w}); + coeffs.push_back({rational::one(), m.var()}); lp::lpvar term_index = c().lra.add_term(coeffs, UINT_MAX); auto* dep = explain_fixed(m, k); term_index = c().lra.map_term_index_to_column_index(term_index); diff --git a/src/math/lp/nla_core.cpp b/src/math/lp/nla_core.cpp index 2ec4a3a10..d9fc43057 100644 --- a/src/math/lp/nla_core.cpp +++ b/src/math/lp/nla_core.cpp @@ -17,7 +17,7 @@ Author: #include "math/grobner/pdd_solver.h" #include "math/dd/pdd_interval.h" #include "math/dd/pdd_eval.h" -namespace nla { +using namespace nla; typedef lp::lar_term term; @@ -1785,20 +1785,17 @@ void core::set_use_nra_model(bool m) { } } -void core::collect_statistics(::statistics & st) { -} - void core::propagate() { clear(); m_monomial_bounds.unit_propagate(); m_monics_with_changed_bounds.reset(); } - void core::simplify() { - // in-processing simplifiation can go here, such as bounds improvements. - } +void core::simplify() { + // in-processing simplifiation can go here, such as bounds improvements. + +} -} // end of nla diff --git a/src/math/lp/nla_core.h b/src/math/lp/nla_core.h index 25156e8b6..29effd4e4 100644 --- a/src/math/lp/nla_core.h +++ b/src/math/lp/nla_core.h @@ -74,8 +74,8 @@ class core { std::function m_relevant; vector m_lemmas; vector m_literals; - vector m_equalities; - vector m_fixed_equalities; + vector m_equalities; + vector m_fixed_equalities; indexed_uint_set m_to_refine; indexed_uint_set m_monics_with_changed_bounds; tangents m_tangents; @@ -420,11 +420,10 @@ public: bool has_real(const monic& m) const; void set_use_nra_model(bool m); bool use_nra_model() const { return m_use_nra_model; } - void collect_statistics(::statistics&); vector const& lemmas() const { return m_lemmas; } vector const& literals() const { return m_literals; } - vector const& equalities() const { return m_equalities; } - vector const& fixed_equalities() const { return m_fixed_equalities; } + vector const& equalities() const { return m_equalities; } + vector const& fixed_equalities() const { return m_fixed_equalities; } bool check_feasible() const { return m_check_feasible; } void add_fixed_equality(lp::lpvar v, rational const& k, lp::explanation const& e) { m_fixed_equalities.push_back({v, k, e}); } diff --git a/src/math/lp/nla_solver.cpp b/src/math/lp/nla_solver.cpp index f4d09810e..4b501a39e 100644 --- a/src/math/lp/nla_solver.cpp +++ b/src/math/lp/nla_solver.cpp @@ -88,10 +88,6 @@ namespace nla { return m_core->m_nra.value(v); } - void solver::collect_statistics(::statistics & st) { - m_core->collect_statistics(st); - } - // ensure r = x^y, add abstraction/refinement lemmas lbool solver::check_power(lpvar r, lpvar x, lpvar y) { return m_core->check_power(r, x, y); @@ -109,11 +105,11 @@ namespace nla { return m_core->literals(); } - vector const& solver::equalities() const { + vector const& solver::equalities() const { return m_core->equalities(); } - vector const& solver::fixed_equalities() const { + vector const& solver::fixed_equalities() const { return m_core->fixed_equalities(); } diff --git a/src/math/lp/nla_solver.h b/src/math/lp/nla_solver.h index 9e650c127..7da05c3e4 100644 --- a/src/math/lp/nla_solver.h +++ b/src/math/lp/nla_solver.h @@ -47,11 +47,10 @@ namespace nla { core& get_core(); nlsat::anum_manager& am(); nlsat::anum const& am_value(lp::var_index v) const; - void collect_statistics(::statistics & st); vector const& lemmas() const; vector const& literals() const; - vector const& fixed_equalities() const; - vector const& equalities() const; + vector const& fixed_equalities() const; + vector const& equalities() const; bool check_feasible() const { return m_core->check_feasible(); } }; } diff --git a/src/math/lp/nla_types.h b/src/math/lp/nla_types.h index 3930a62a9..cb62c5a30 100644 --- a/src/math/lp/nla_types.h +++ b/src/math/lp/nla_types.h @@ -25,18 +25,6 @@ namespace nla { typedef lp::var_index lpvar; const lpvar null_lpvar = UINT_MAX; - struct equality { - lp::lpvar i, j; - lp::explanation e; - equality(lp::lpvar i, lp::lpvar j, lp::explanation const& e):i(i),j(j),e(e) {} - }; - - struct fixed_equality { - lp::lpvar v; - rational k; - lp::explanation e; - fixed_equality(lp::lpvar v, rational const& k, lp::explanation const& e):v(v),k(k),e(e) {} - }; inline int rat_sign(const rational& r) { return r.is_pos()? 1 : ( r.is_neg()? -1 : 0); } diff --git a/src/sat/smt/arith_diagnostics.cpp b/src/sat/smt/arith_diagnostics.cpp index e9d989545..47e3ee551 100644 --- a/src/sat/smt/arith_diagnostics.cpp +++ b/src/sat/smt/arith_diagnostics.cpp @@ -90,7 +90,6 @@ namespace arith { void solver::collect_statistics(statistics& st) const { m_stats.collect_statistics(st); lp().settings().stats().collect_statistics(st); - if (m_nla) m_nla->collect_statistics(st); } void solver::explain_assumptions(lp::explanation const& e) { diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp index 4de1e732b..fffa85d44 100644 --- a/src/smt/theory_lra.cpp +++ b/src/smt/theory_lra.cpp @@ -1091,6 +1091,26 @@ public: void restart_eh() { m_arith_eq_adapter.restart_eh(); +#if 0 + // experiment + if (m_lia) { + std::function is_root = [&](unsigned j) { + theory_var v = lp().local_to_external(j); + if (v < 0) + return false; + auto* n = get_enode(v); + if (!th.is_relevant_and_shared(n)) + return false; + if (n->is_root()) + return true; + theory_var w = n->get_root()->get_th_var(get_id()); + return w == v; + }; + m_lia->simplify(is_root); + for (auto const& [i, j, e] : m_lia->equalities()) + add_eq(i, j, e, false); + } +#endif if (m_nla) m_nla->simplify(); } @@ -3843,7 +3863,6 @@ public: m_arith_eq_adapter.collect_statistics(st); m_stats.collect_statistics(st); lp().settings().stats().collect_statistics(st); - if (m_nla) m_nla->collect_statistics(st); } /*