add simplification experiment (disabled) for tracking, some reshuffling of equation/fixed_equation structs

Signed-off-by: Nikolaj Bjorner <nbjorner@microsoft.com>

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) {}
+    };
+
 }

src/math/lp/int_solver.cpp

@@ -683,5 +683,236 @@ namespace lp {
         return result;
     }
 
+    void int_solver::simplify(std::function<bool(unsigned)>& 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<rational, unsigned_vector, rational::hash_proc, rational::eq_proc> value2roots;
+
+        vector<std::pair<lp::mpq, unsigned>> 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<rational> 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
+    }
+
 
 }

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<equality>       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<bool(unsigned)>& is_root);
+    vector<equality> 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);

src/math/lp/lp_types.h

@@ -92,6 +92,7 @@ public:
 
 };
 
+
 }
 
 inline std::ostream& operator<<(std::ostream& out, lp::tv const& t) {

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<std::pair<lp::mpq, unsigned>> 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);

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
 

src/math/lp/nla_core.h

@@ -74,8 +74,8 @@ class core {
     std::function<bool(lpvar)> m_relevant;
     vector<lemma>            m_lemmas;
     vector<ineq>             m_literals;
-    vector<equality>         m_equalities;
-    vector<fixed_equality>   m_fixed_equalities;
+    vector<lp::equality>       m_equalities;
+    vector<lp::fixed_equality> 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<nla::lemma> const& lemmas() const { return m_lemmas; }
     vector<nla::ineq> const& literals() const { return m_literals; }
-    vector<equality> const& equalities() const { return m_equalities; }
-    vector<fixed_equality> const& fixed_equalities() const { return m_fixed_equalities; }
+    vector<lp::equality> const& equalities() const { return m_equalities; }
+    vector<lp::fixed_equality> 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}); }

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<nla::equality> const& solver::equalities() const {
+    vector<lp::equality> const& solver::equalities() const {
         return m_core->equalities();
     }
 
-    vector<nla::fixed_equality> const& solver::fixed_equalities() const {
+    vector<lp::fixed_equality> const& solver::fixed_equalities() const {
         return m_core->fixed_equalities();
     }
 

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<nla::lemma> const& lemmas() const;
         vector<nla::ineq> const& literals() const;
-        vector<nla::fixed_equality> const& fixed_equalities() const;
-        vector<nla::equality> const& equalities() const;
+        vector<lp::fixed_equality> const& fixed_equalities() const;
+        vector<lp::equality> const& equalities() const;
         bool check_feasible() const { return m_core->check_feasible(); }
     };
 }

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); }

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) {

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<bool(unsigned)> 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);
     }        
 
     /*