experiment with cuts in dio

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

src/math/lp/bound_analyzer_on_row.h

@@ -148,8 +148,7 @@ namespace lp {
         mpq m_total, m_bound;
         void limit_all_monoids_from_above() {
             int strict = 0;
-            m_total.reset();
+            m_total = m_rs.x;
-            lp_assert(is_zero(m_total));
             for (const auto& p : m_row) {
                 bool str;
                 m_total -= monoid_min(p.coeff(), p.var(), str);
@@ -174,8 +173,7 @@ namespace lp {
 
         void limit_all_monoids_from_below() {
             int strict = 0;
-            m_total.reset();
+            m_total = m_rs.x;
-            lp_assert(is_zero(m_total));
             for (const auto &p : m_row) {
                 bool str;
                 m_total -= monoid_max(p.coeff(), p.var(), str);
@@ -205,7 +203,7 @@ namespace lp {
             // every other monoid is impossible to limit from below
             mpq u_coeff;
             unsigned j;
-            m_bound = -m_rs.x;
+            m_bound = m_rs.x;
             bool strict = false;
             for (const auto& p : m_row) {
                 j = p.var();
@@ -234,7 +232,7 @@ namespace lp {
             // every other monoid is impossible to limit from above
             mpq l_coeff;
             unsigned j;
-            m_bound = -m_rs.x;
+            m_bound = m_rs.x;
             bool strict = false;
             for (const auto &p : m_row) {
                 j = p.var();
@@ -275,31 +273,25 @@ namespace lp {
         //     */
         // }
 
-        void limit_j(unsigned bound_j, const mpq& u, bool coeff_before_j_is_pos, bool is_lower_bound, bool strict)
+        void limit_j(unsigned bound_j, const mpq& u, bool coeff_before_j_is_pos, bool is_lower_bound, bool strict) {
-            {
+            auto* lar = &m_bp.lp();
-                auto* lar = &m_bp.lp();
+            int bound_sign = (is_lower_bound ? 1 : -1);
-                const auto& row = this->m_row;
+            int j_sign = (coeff_before_j_is_pos ? 1 : -1) * bound_sign;
-                auto explain = [row, bound_j, coeff_before_j_is_pos, is_lower_bound, strict, lar]() {
+                
-                    (void) strict;
+            u_dependency* dep = nullptr;
-                    TRACE("bound_analyzer", tout << "explain_bound_on_var_on_coeff, bound_j = " << bound_j << ", coeff_before_j_is_pos = " << coeff_before_j_is_pos << ", is_lower_bound = " << is_lower_bound << ", strict = " << strict << "\n";);
+            for (auto const& r : m_row) {
-                    int bound_sign = (is_lower_bound ? 1 : -1);
+                unsigned j = r.var();
-                    int j_sign = (coeff_before_j_is_pos ? 1 : -1) * bound_sign;
+                if (j == bound_j)
-
+                    continue;
-                    u_dependency* ret = nullptr;
+                mpq const& a = r.coeff();
-                    for (auto const& r : row) {
+                int a_sign = is_pos(a) ? 1 : -1;
-                        unsigned j = r.var();
+                int sign = j_sign * a_sign;
-                        if (j == bound_j)
+                u_dependency* witness = sign > 0 ? lar->get_column_upper_bound_witness(j) : lar->get_column_lower_bound_witness(j);
-                            continue;
+                dep = lar->join_deps(dep, witness);
-                        mpq const& a = r.coeff();
-                        int a_sign = is_pos(a) ? 1 : -1;
-                        int sign = j_sign * a_sign;
-                        u_dependency* witness = sign > 0 ? lar->get_column_upper_bound_witness(j) : lar->get_column_lower_bound_witness(j);
-                        ret = lar->join_deps(ret, witness);
-                    }
-                    return ret;
-                };
-                m_bp.add_bound(u, bound_j, is_lower_bound, strict, explain);
             }
+                
+            m_bp.add_bound(u, bound_j, is_lower_bound, strict, dep);
+        }
 
         void advance_u(unsigned j) {
             m_column_of_u = (m_column_of_u == -1) ? j : -2;

src/math/lp/dioph_eq.cpp

@@ -245,27 +245,16 @@ namespace lp {
         std::ostream& print_S(std::ostream& out) {
             out << "S:\n";
             for (const auto& p : m_k2s.m_map) {
-                print_entry(p.second, out, false, false);
+                print_entry(p.second, out, false, false, true);
             }
             return out;
         }
-
-        std::ostream& print_column_bounds(std::ostream& out, unsigned j) {
-            out << "j" << j << ":= " << lra.get_column_value(j).x << ": ";
-            if (lra.column_has_lower_bound(j))
-                out << lra.column_lower_bound(j).x << " <= ";
-            out << "x" << j;
-            if (lra.column_has_upper_bound(j))
-                out << " <= " << lra.column_upper_bound(j).x;
-            out << "\n";
-            return out;
-        }
             
         std::ostream& print_bounds(std::ostream& out) {
             out << "bounds:\n";
             for (unsigned v = 0; v < m_var_register.size(); ++v) {
                 unsigned j = m_var_register.local_to_external(v);
-                print_column_bounds(out, j);
+                lra.print_column_info(j, out);
             }
             return out;
         }
@@ -1446,12 +1435,9 @@ namespace lp {
             m_c = mpq(0);
             m_lspace.clear();
             SASSERT(get_extended_term_value(lar_t).is_zero());
-            TRACE("dio", tout << "t:";print_lar_term_L(lar_t, tout) << std::endl; tout << "value:" << get_extended_term_value(lar_t) << std::endl;);
             for (const auto& p : lar_t) {
-                SASSERT(p.coeff().is_int());
                 if (is_fixed(p.j())) {
                     m_c += p.coeff() * lia.lower_bound(p.j()).x;
-                    SASSERT(lia.lower_bound(p.j()).x == lra.get_column_value(p.j()).x); 
                 }
                 else {
                     unsigned lj = lar_solver_to_local(p.j());
@@ -1460,9 +1446,6 @@ namespace lp {
                         q.push(lar_solver_to_local(p.j()));
                 }
             }
-            TRACE("dio", tout << "m_espace:"; print_term_o(create_term_from_subst_space(), tout) << std::endl;
-                  tout << "m_c:" << m_c << std::endl;
-                  tout << "get_value_of_subs_space:" << get_value_of_espace() << std::endl;);
         }
         unsigned lar_solver_to_local(unsigned j) const {
             return m_var_register.external_to_local(j);
@@ -1637,7 +1620,136 @@ namespace lp {
             }
         }
 
+        struct prop_bound {
+            mpq m_bound;
+            unsigned m_j;
+            bool m_is_low;
+            bool m_strict;
+            u_dependency* m_dep;
+        };
+
+        std_vector<prop_bound> m_prop_bounds;
+
+
+        void add_bound(mpq const& bound, unsigned j, bool is_low, bool strict, u_dependency* dep) {
+            m_prop_bounds.push_back({bound, j, is_low, strict, dep});
+        }
+
+        lconstraint_kind get_bound_kind(const prop_bound& pb) const {
+            if (!pb.m_is_low) {
+                return pb.m_strict ? lp::LT : lp::LE; 
+            }
+            else {
+                return pb.m_strict ? lp::GT : lp::GE;
+            }
+        }
+
+        bool cut(const prop_bound& pb) const {            
+            const auto & v = lra.get_column_value(pb.m_j);
+            switch (get_bound_kind(pb)) {
+            case LT:
+                return v >= pb.m_bound;
+            case LE:
+                return v > pb.m_bound;
+            case GT:
+                return v <= pb.m_bound;
+            case GE:
+                return v < pb.m_bound;
+            default:
+                UNREACHABLE();
+            }
+            return false;
+        }
+
+        std::ostream& print_int_inf_vars(std::ostream& out) const {
+            out << "Integer infeasible variables:\n";
+            for (unsigned j = 0; j < lra.column_count(); j++) {
+                if (!lia.column_is_int_inf(j)) continue;
+                out << "x" << j << " = " << lra.get_column_value(j).x << " bounds:\n";
+                if (lra.column_has_lower_bound(j))
+                    out << lra.get_lower_bound(j).x << " <= ";
+                else 
+                    out << "-oo" << " <= ";
+
+                out << "x" << j;
+                if (lra.column_has_upper_bound(j)) 
+                    out << " <= " << lra.get_upper_bound(j).x;
+                else 
+                    out << " <= " << "oo";    
+                out << "\n";
+                
+            }
+            return out;
+        }
+
+        // return true iff lar_solvre is feasible or cancelled
+        bool create_cut(const prop_bound& pb) {
+            lra.update_column_type_and_bound(pb.m_j, get_bound_kind(pb), pb.m_bound, pb.m_dep);
+            auto st = lra.find_feasible_solution();
+            if ((int)st >= (int)(lp_status::FEASIBLE) || !lia.settings().get_cancel_flag())
+                return true;
+            if (st == lp_status::INFEASIBLE) {
+                lra.get_infeasibility_explanation(m_infeas_explanation);
+                return false;
+            }
+            return true;
+        }
+
+        bool row_has_int_inf(unsigned ei) {
+            for (const auto& p: m_e_matrix.m_rows[ei]) {
+                unsigned j = local_to_lar_solver(p.var());
+                if (lia.column_is_int_inf(j))
+                    return true;
+            }
+            return false;
+        }
+        int prop_calls = 0;
+        lia_move tighten_S_row(unsigned ei) {
+            if (!row_has_int_inf(ei)) return lia_move::undef;
+            // SASSERT(entry_invariant(ei));
+            m_espace.clear();
+            for (const auto& p: m_e_matrix.m_rows[ei]) {
+                m_espace.add(p.coeff(), p.var());
+            }
+            remove_fresh_from_espace();
+            std_vector<iv> row_in_lar_solver_indices;
+            for (const auto & p: m_espace.m_data) {
+                row_in_lar_solver_indices.push_back({p.coeff(), local_to_lar_solver(p.var())});
+            }
+            m_prop_bounds.clear();
+            bound_analyzer_on_row<std_vector<iv>, dioph_eq::imp>::analyze_row(row_in_lar_solver_indices, impq(- m_sum_of_fixed[ei]), *this);
+            for (auto& pb: m_prop_bounds) {
+                if (lra.settings().get_cancel_flag())
+                    return lia_move::undef;
+
+                if (cut(pb) && lia.column_is_int_inf(pb.m_j)) {
+                    pb.m_dep = lra.join_deps(pb.m_dep, explain_fixed_in_meta_term(m_l_matrix.m_rows[ei])); // todo; not efficient
+                    TRACE("dio", tout << "cut on variable x" << pb.m_j << ", with " << (pb.m_is_low? "lower": "upper") << " bound:" << pb.m_bound << "\n"; tout << "row " << ei << " in lar indices:\n";
+                          print_lar_term_L(row_in_lar_solver_indices, tout) << " + " << m_sum_of_fixed[ei] << "\n";
+                          tout << "bounds:\n";
+                          for (const auto& p: row_in_lar_solver_indices)
+                              lra.print_column_info(p.var(), tout););
+                    if (!create_cut(pb))
+                        return lia_move::conflict;                    
+                }
+            }
+            return lia_move::undef;
+        }
+        
         lia_move propagate_bounds_on_tightened_columns() {
+            TRACE("dio", print_int_inf_vars(tout); print_S(tout););
+            std::unordered_set<unsigned> rows;
+            for (unsigned tcol : m_tightened_columns) {
+                for (const auto& p: m_e_matrix.m_columns[lar_solver_to_local(tcol)] ) {
+                    if (lra.settings().get_cancel_flag())
+                        return lia_move::undef;
+                    if (contains(rows, p.var())) continue;
+                    rows.insert(p.var());
+                    lia_move r = tighten_S_row(p.var());
+                    if (r == lia_move::conflict)
+                        return r;                        
+                }
+            }
             return lia_move::undef;
         }
         // m_espace contains the coefficients of the term
@@ -1751,6 +1863,9 @@ namespace lp {
                 lra.stats().m_dio_rewrite_conflicts++;
                 return lia_move::conflict;
             }
+            TRACE("dio", print_int_inf_vars(tout) << "\n"; 
+                 print_S(tout));
+
             return lia_move::undef;
         }
 
@@ -1763,7 +1878,14 @@ namespace lp {
             if (ret == lia_move::conflict) {
                 lra.stats().m_dio_tighten_conflicts++;
                 return lia_move::conflict;
+            }           
+            if (ret == lia_move::undef) {
+                TRACE("dio", print_int_inf_vars(tout); print_S(tout););
+                ret = propagate_bounds_on_tightened_columns();
+                if (ret == lia_move::conflict)
+                    lra.stats().m_dio_bound_propagation_conflicts++;
             }
+
             return ret;
         }
 
@@ -1920,7 +2042,7 @@ namespace lp {
                     lra_pop();
                     continue;
                 }
-                auto st = lra.find_feasible_solution();
+                auto st = lra.find_feasible_solution();                
                 TRACE("dio_br", tout << "st:" << lp_status_to_string(st) << std::endl;);
                 if ((int)st >= (int)(lp_status::FEASIBLE)) {
                     // have a feasible solution
@@ -2285,17 +2407,13 @@ namespace lp {
         }
 
         bool entry_invariant(unsigned ei) const {
-            if (belongs_to_s(ei)) {
-                auto it = m_k2s.m_rev_map.find(ei);
-                SASSERT(it != m_k2s.m_rev_map.end());
-                unsigned j = it->second;
-                get_sign_in_e_row(ei, j);
-            }
-
             for (const auto& p : m_e_matrix.m_rows[ei]) {
                 if (!p.coeff().is_int()) {
                     return false;
                 }
+                unsigned j = local_to_lar_solver(p.var());
+                if (is_fixed(j))
+                    return false;
             }
 
             bool ret =
@@ -2432,12 +2550,13 @@ namespace lp {
         }
 
         std::ostream& print_entry(unsigned i, std::ostream& out,
-                                  bool print_dep = false, bool print_in_S = true) {
+                                  bool print_dep = false, bool print_in_S = true, bool print_column_info = false) {
             unsigned j = m_k2s.has_val(i) ? m_k2s.get_key(i) : UINT_MAX;
             out << "entry " << i << ": ";
             if (j != UINT_MAX) 
                 out << "x" << j << " ";               
             out << "{\n";
+            
             print_term_o(get_term_from_entry(i), out << "\t") << ",\n";
             // out << "\tstatus:" << (int)e.m_entry_status;
             if (print_dep) {
@@ -2453,13 +2572,15 @@ namespace lp {
                 out << "in F\n";
             }
             else {
-                if (local_to_lar_solver(j) == UINT_MAX) {
+                if (print_in_S) {
-                    out << "FRESH\n";
-                }
-                else if (print_in_S) {
                     out << "in S\n";
                 }
             }
+            if (print_column_info) {
+                for (const auto& p : get_term_from_entry(i)) {
+                    lra.print_column_info(local_to_lar_solver(p.var()), out) << "\n";
+                }
+            }
             out << "}\n";
             return out;
         }

src/math/lp/implied_bound.h

@@ -33,11 +33,11 @@ class implied_bound {
     bool m_is_lower_bound;
     bool m_strict;
     private:
-    std::function<u_dependency*()> m_explain_bound = nullptr;
+    u_dependency* m_dep = nullptr;
     public:
     // s is expected to be the pointer to lp_bound_propagator.
-    u_dependency* explain_implied() const { return m_explain_bound(); }
+    u_dependency* dep() const { return m_dep; }
-    void set_explain(std::function<u_dependency*()> f) { m_explain_bound = f; }
+    void set_dep(u_dependency* dep) { m_dep = dep; }
     lconstraint_kind kind() const {
         lconstraint_kind k = m_is_lower_bound? GE : LE;
         if (m_strict)
@@ -49,12 +49,12 @@ class implied_bound {
                   unsigned j,
                   bool is_lower_bound,
                   bool is_strict,
-                  std::function<u_dependency*()> get_dep):
+                  u_dependency* dep):
         m_bound(a),
         m_j(j),
         m_is_lower_bound(is_lower_bound),
         m_strict(is_strict),
-        m_explain_bound(get_dep) {
+        m_dep(dep) {
     }
 };
 }

src/math/lp/int_solver.cpp

@@ -227,8 +227,8 @@ namespace lp {
 
         unsigned m_bounds_refine_depth = 0;
        
-        void add_bound(mpq const& bound, unsigned j, bool is_low, bool strict, const std::function<u_dependency* ()>& explain_bound) {
+        void add_bound(mpq const& bound, unsigned j, bool is_low, bool strict, u_dependency* dep) {
-            m_prop_bounds.push_back({bound, j, is_low, strict, explain_bound()});
+            m_prop_bounds.push_back({bound, j, is_low, strict, dep});
         }        
 
         lconstraint_kind get_bound_kind(bool upper, bool is_strict) {
@@ -280,7 +280,7 @@ namespace lp {
             return lia_move::undef;
         }
         bool should_tighten_bounds() {
-            return m_number_of_calls % 4 == 0;
+            return false && m_number_of_calls % 4 == 0;
         }
      // needed for the template bound_analyzer_on_row.h
         const lar_solver& lp() const { return lra; }

src/math/lp/lar_solver.h

@@ -312,7 +312,7 @@ public:
     std::function<void(lpvar)> m_fixed_var_eh;
     template <typename T>
     void explain_implied_bound(const implied_bound& ib, lp_bound_propagator<T>& bp) {
-        u_dependency* dep = ib.explain_implied();
+        u_dependency* dep = ib.dep();
         for (auto ci : flatten(dep))
             bp.consume(mpq(1), ci); // TODO: flatten should provide the coefficients
         /*

src/math/lp/lp_bound_propagator.h

@@ -125,7 +125,7 @@ public:
         return (*m_column_types)[j] == column_type::fixed && get_lower_bound(j).y.is_zero();
     }
 
-    void add_bound(mpq const& v, unsigned j, bool is_low, bool strict, std::function<u_dependency* ()> explain_bound) {
+    void add_bound(mpq const& v, unsigned j, bool is_low, bool strict, u_dependency* dep) {
         lconstraint_kind kind = is_low ? GE : LE;
         if (strict)
             kind = static_cast<lconstraint_kind>(kind / 2);
@@ -140,12 +140,12 @@ public:
 
                     found_bound.m_bound = v;
                     found_bound.m_strict = strict;
-                    found_bound.set_explain(explain_bound);
+                    found_bound.set_dep(dep);
                     TRACE("add_bound", lp().print_implied_bound(found_bound, tout););
                 }
             } else {
                 m_improved_lower_bounds.insert(j, static_cast<unsigned>(m_ibounds.size()));
-                m_ibounds.push_back(implied_bound(v, j, is_low, strict, explain_bound));
+                m_ibounds.push_back(implied_bound(v, j, is_low, strict, dep));
                 TRACE("add_bound", lp().print_implied_bound(m_ibounds.back(), tout););
             }
         } else {  // the upper bound case
@@ -156,12 +156,12 @@ public:
 
                     found_bound.m_bound = v;
                     found_bound.m_strict = strict;
-                    found_bound.set_explain(explain_bound);
+                    found_bound.set_dep(dep);
                     TRACE("add_bound", lp().print_implied_bound(found_bound, tout););
                 }
             } else {
                 m_improved_upper_bounds.insert(j, static_cast<unsigned>(m_ibounds.size()));
-                m_ibounds.push_back(implied_bound(v, j, is_low, strict, explain_bound));
+                m_ibounds.push_back(implied_bound(v, j, is_low, strict, dep));
                 TRACE("add_bound", lp().print_implied_bound(m_ibounds.back(), tout););
             }
         }

src/math/lp/lp_settings.h

@@ -138,6 +138,7 @@ struct statistics {
     unsigned m_dio_rewrite_conflicts = 0;
     unsigned m_dio_branching_sats = 0;
     unsigned m_dio_branching_conflicts = 0;
+    unsigned m_dio_bound_propagation_conflicts = 0;
     unsigned m_bounds_tightening_conflicts = 0;
     unsigned m_bounds_tightenings = 0;
     ::statistics m_st = {};
@@ -184,6 +185,7 @@ struct statistics {
         st.update("arith-dio-branching-depth", m_dio_branching_depth);
         st.update("arith-dio-branching-conflicts", m_dio_branching_conflicts);
         st.update("arith-bounds-tightening-conflicts", m_bounds_tightening_conflicts);
+        st.update("arith-dio-propagation-conflicts", m_dio_bound_propagation_conflicts);
         st.update("arith-bounds-tightenings", m_bounds_tightenings);
         st.copy(m_st);
     }