fix bug in bound simplification in Gomory for mixed integer linear cuts, enable fixed variable redution after bugfix, add notes that rewriting bounds does not work

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

src/math/lp/gomory.cpp

@@ -32,7 +32,6 @@ class create_cut {
     unsigned              m_inf_col; // a basis column which has to be an integer but has a non integral value
     const row_strip<mpq>& m_row;
     int_solver&           lia;
-    mpq                   m_lcm_den = { mpq(1) };
     mpq                   m_f;
     mpq                   m_one_minus_f;
     mpq                   m_fj;
@@ -131,120 +130,7 @@ class create_cut {
         // conflict 0 >= k where k is positive
         return lia_move::conflict;
     }
-
-    void divd(mpq& r, mpq const& d) {
-        r /= d;
-        if (!r.is_int())
-            r = ceil(r);
-    }
-
-    bool can_divide_by(vector<std::pair<mpq, lpvar>> const& p, mpq const& d) {
-        mpq lhs(0), rhs(m_k);
-        mpq max_c(abs(m_k));
-        for (auto const& [c, v] : p) {            
-            auto c1 = c;
-            max_c = std::max(max_c, abs(c1));
-            divd(c1, d);
-            if (c1 == 0)
-                return false;
-            VERIFY(lia.get_value(v).y == 0);
-            lhs += c1 * lia.get_value(v).x;
-        }
-        if (max_c == 1)
-            return false;
-        divd(rhs, d);
-        return lhs < rhs;        
-    }
     
-    void adjust_term_and_k_for_some_ints_case_gomory() {
-        lp_assert(!m_t.is_empty());
-        // k = 1 + sum of m_t at bounds
-        lar_term t = lia.lra.unfold_nested_subterms(m_t);
-        auto pol = t.coeffs_as_vector();
-        m_t.clear();
-        if (pol.size() == 1) {
-            TRACE("gomory_cut_detail", tout << "pol.size() is 1" << std::endl;);
-            auto const& [a, v] = pol[0];
-            lp_assert(is_int(v));
-            if (a.is_pos()) { // we have av >= k
-                divd(m_k, a);
-                m_t.add_monomial(mpq(1), v);
-            }
-            else {
-                // av >= k
-                // a/-a*v >= k / - a
-                // -v >= k / - a
-                // -v >= ceil(k / -a)
-                divd(m_k, -a);
-                m_t.add_monomial(-mpq(1), v);
-            }
-        }
-        else {
-            m_lcm_den = denominator(m_k);
-            for (auto const& [c, v] : pol)
-                m_lcm_den = lcm(m_lcm_den, denominator(c));
-            lp_assert(m_lcm_den.is_pos());
-            TRACE("gomory_cut_detail", tout << "pol.size() > 1 den: " << m_lcm_den << std::endl;);
-            if (!m_lcm_den.is_one()) {
-                // normalize coefficients of integer parameters to be integers.
-                for (auto & [c,v]: pol) {
-                    c *= m_lcm_den;
-                    SASSERT(!is_int(v) || c.is_int());
-                }
-                m_k *= m_lcm_den;
-            }
-#if 0
-            unsigned j = 0, i = 0;
-            for (auto & [c, v] : pol) {
-                if (lia.is_fixed(v)) {
-                    push_explanation(column_lower_bound_constraint(v));
-                    push_explanation(column_upper_bound_constraint(v));
-                    m_k -= c;
-                    IF_VERBOSE(0, verbose_stream() << "got fixed " << v << "\n");
-                }
-                else
-                    pol[j++] = pol[i];
-                ++i;
-            }
-            pol.shrink(j);
-#endif
-
-            // gcd reduction is loss-less:
-            mpq g(1);
-            for (const auto & [c, v] : pol)
-                g = gcd(g, c);
-
-            if (g != 1) {
-                for (auto & [c, v] : pol)
-                    c /= g;
-                divd(m_k, g);
-            }
-                
-#if 0
-            // TODO: create self-contained rounding mode to weaken cuts
-            // whose cofficients are considered too large
-            // (larger than bounds from the input)
-            mpq min_c = abs(m_k);
-            for (const auto & [c, v] : pol) 
-                min_c = std::min(min_c, abs(c));
-                
-            if (min_c > 1 && can_divide_by(pol, min_c)) {
-                for (auto& [c, v] : pol)
-                    divd(c, min_c);
-                divd(m_k, min_c);
-            }
-#endif
-
-            for (const auto & [c, v]: pol) 
-                m_t.add_monomial(c, v);
-            VERIFY(m_t.size() > 0);
-        }
-
-        TRACE("gomory_cut_detail", tout << "k = " << m_k << std::endl;);
-        lp_assert(m_k.is_int());
-
-        
-    }
 
     std::string var_name(unsigned j) const {
         return std::string("x") + std::to_string(j);
@@ -359,12 +245,12 @@ public:
         // gomory will be   t >= k and the current solution has a property t < k
         m_k = 1;
         m_t.clear();
-        bool some_int_columns = false;
         mpq m_f  = fractional_part(get_value(m_inf_col));
         TRACE("gomory_cut_detail", tout << "m_f: " << m_f << ", ";
               tout << "1 - m_f: " << 1 - m_f << ", get_value(m_inf_col).x - m_f = " << get_value(m_inf_col).x - m_f << "\n";);
         lp_assert(m_f.is_pos() && (get_value(m_inf_col).x - m_f).is_int());  
 
+        bool some_int_columns = false;
 #if SMALL_CUTS
         m_abs_max = 0;
         for (const auto & p : m_row) {
@@ -408,13 +294,7 @@ public:
         if (m_t.is_empty())
             return report_conflict_from_gomory_cut();
         if (some_int_columns)
-            adjust_term_and_k_for_some_ints_case_gomory();
+            simplify_inequality();
-        if (!lia.current_solution_is_inf_on_cut()) {
-            m_ex->clear();
-            m_t.clear();
-            m_k = 1;
-            return lia_move::undef;
-        }
         lp_assert(lia.current_solution_is_inf_on_cut());  // checks that indices are columns
         TRACE("gomory_cut", print_linear_combination_of_column_indices_only(m_t.coeffs_as_vector(), tout << "gomory cut: "); tout << " >= " << m_k << std::endl;);
         TRACE("gomory_cut_detail", dump_cut_and_constraints_as_smt_lemma(tout);
@@ -423,6 +303,91 @@ public:
         return lia_move::cut;
     }
 
+    // TODO: use this also for HNF cuts?
+    mpq                 m_lcm_den = { mpq(1) };
+    
+    void simplify_inequality() {
+
+        auto divd = [](mpq& r, mpq const& d) {
+            r /= d;
+            if (!r.is_int())
+                r = ceil(r);
+        };
+        SASSERT(!lia.m_upper);
+        lp_assert(!m_t.is_empty());
+        // k = 1 + sum of m_t at bounds
+        lar_term t = lia.lra.unfold_nested_subterms(m_t);
+        auto pol = t.coeffs_as_vector();
+        m_t.clear();
+        if (pol.size() == 1) {
+            TRACE("gomory_cut_detail", tout << "pol.size() is 1" << std::endl;);
+            auto const& [a, v] = pol[0];
+            lp_assert(is_int(v));
+            if (a.is_pos()) { // we have av >= k
+                divd(m_k, a);
+                m_t.add_monomial(mpq(1), v);
+            }
+            else {
+                // av >= k
+                // a/-a*v >= k / - a
+                // -v >= k / - a
+                // -v >= ceil(k / -a)
+                divd(m_k, -a);
+                m_t.add_monomial(-mpq(1), v);
+            }
+        }
+        else {
+            m_lcm_den = denominator(m_k);
+            for (auto const& [c, v] : pol)
+                m_lcm_den = lcm(m_lcm_den, denominator(c));
+            lp_assert(m_lcm_den.is_pos());
+            bool int_row = true;
+            TRACE("gomory_cut_detail", tout << "pol.size() > 1 den: " << m_lcm_den << std::endl;);
+            if (!m_lcm_den.is_one()) {
+                // normalize coefficients of integer parameters to be integers.
+                for (auto & [c,v]: pol) {
+                    c *= m_lcm_den;
+                    SASSERT(!is_int(v) || c.is_int());
+                    int_row &= is_int(v);
+                }
+                m_k *= m_lcm_den;
+            }
+            unsigned j = 0, i = 0;
+            for (auto & [c, v] : pol) {
+                if (lia.is_fixed(v)) {
+                    push_explanation(column_lower_bound_constraint(v));
+                    push_explanation(column_upper_bound_constraint(v));
+                    m_k -= c * lower_bound(v).x;
+                }
+                else
+                    pol[j++] = pol[i];
+                ++i;
+            }
+            pol.shrink(j);
+
+            // gcd reduction is loss-less:
+            mpq g(1);
+            for (const auto & [c, v] : pol)
+                g = gcd(g, c);
+            if (!int_row)
+                g = gcd(g, m_k);
+
+            if (g != 1) {
+                for (auto & [c, v] : pol)
+                    c /= g;
+                divd(m_k, g);
+            }
+                
+            for (const auto & [c, v]: pol) 
+                m_t.add_monomial(c, v);
+            VERIFY(m_t.size() > 0);
+        }
+
+        TRACE("gomory_cut_detail", tout << "k = " << m_k << std::endl;);
+        lp_assert(m_k.is_int());
+    }
+
+
     create_cut(lar_term & t, mpq & k, explanation* ex, unsigned basic_inf_int_j, const row_strip<mpq>& row, int_solver& lia) :
         m_t(t),
         m_k(k),

src/math/lp/hnf_cutter.cpp

@@ -83,14 +83,13 @@ namespace lp  {
             // consider return from here if b[i] is not an integer and return i
         }
     }
-
+    
     vector<mpq> hnf_cutter::create_b(const svector<unsigned> & basis_rows) {
         if (basis_rows.size() == m_right_sides.size())
             return m_right_sides;
         vector<mpq> b;
-        for (unsigned i : basis_rows) {
+        for (unsigned i : basis_rows) 
-            b.push_back(m_right_sides[i]);
+            b.push_back(m_right_sides[i]);        
-        }
         return b;
     }
 
@@ -98,16 +97,15 @@ namespace lp  {
         int ret = -1;
         int n = 0;
         for (int i = 0; i < static_cast<int>(b.size()); i++) {
-            if (is_integer(b[i])) continue;
+            if (is_integer(b[i]))
-            if (n == 0 ) {
+                continue;
+            if (n == 0) {
                 lp_assert(ret == -1);
                 n = 1;
                 ret = i;
-            } else {
-                if (m_settings.random_next() % (++n) == 0) {
-                    ret = i;
-                }
             }
+            else if (m_settings.random_next() % (++n) == 0) 
+                ret = i;            
         }
         return ret;
     }
@@ -240,10 +238,7 @@ branch y_i >= ceil(y0_i) is impossible.
     }
 
     bool hnf_cutter::hnf_has_var_with_non_integral_value() const {
-        for (unsigned j : vars())
+        return any_of(vars(), [&](unsigned j) { return !lia.get_value(j).is_int(); });
-            if (!lia.get_value(j).is_int())
-                return true;
-        return false;
     }
 
     bool hnf_cutter::init_terms_for_hnf_cut() {
@@ -252,17 +247,15 @@ branch y_i >= ceil(y0_i) is impossible.
             try_add_term_to_A_for_hnf(tv::term(i));
         return hnf_has_var_with_non_integral_value();
     }
-
+    
     lia_move hnf_cutter::make_hnf_cut() {
-        if (!init_terms_for_hnf_cut()) {
+        if (!init_terms_for_hnf_cut()) 
             return lia_move::undef;
-        }
         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 (u_dependency* d : constraints_for_explanation()) {
+              for (u_dependency* d : constraints_for_explanation()) 
-                for (auto ci : lra.flatten(d))
+                  for (auto ci : lra.flatten(d))
-                  lra.constraints().display(tout, ci);
+                      lra.constraints().display(tout, ci);
-              }              
               tout << lra.constraints();
               );
 #ifdef Z3DEBUG
@@ -273,14 +266,14 @@ branch y_i >= ceil(y0_i) is impossible.
                                  , x0
 #endif
                                  );
-
+        
         if (r == lia_move::cut) {      
             TRACE("hnf_cut",
                   lra.print_term(lia.m_t, tout << "cut:"); 
                   tout << " <= " << lia.m_k << std::endl;
                   for (auto* dep : constraints_for_explanation()) 
                       for (auto ci : lra.flatten(dep))
-                        lra.constraints().display(tout, ci);                  
+                          lra.constraints().display(tout, ci);                  
                   );
             lp_assert(lia.current_solution_is_inf_on_cut());
             lia.settings().stats().m_hnf_cuts++;
@@ -291,5 +284,4 @@ branch y_i >= ceil(y0_i) is impossible.
         } 
         return r;
     }
-
 }

src/math/lp/int_solver.cpp

@@ -197,11 +197,14 @@ namespace lp {
         if (r == lia_move::undef && should_find_cube()) r = int_cube(*this)();
         if (r == lia_move::undef) lra.move_non_basic_columns_to_bounds();
         if (r == lia_move::undef && should_hnf_cut()) r = hnf_cut();
+
+        std::function<lia_move(void)> gomory_fn = [&]() { return gomory(*this)(); };
         m_cut_vars.reset();
 #if 0
         if (r == lia_move::undef && should_gomory_cut()) r = gomory(*this)();
 #else
-        if (r == lia_move::undef && should_gomory_cut()) r = local_gomory(2);
+        if (r == lia_move::undef && should_gomory_cut()) r = local_cut(2, gomory_fn);
+        
 #endif
         m_cut_vars.reset();
         if (r == lia_move::undef) r = int_branch(*this)();
@@ -711,6 +714,8 @@ namespace lp {
 
         return;
 
+#if 0
+
         // in-processing simplification can go here, such as bounds improvements.
 
         if (!lra.is_feasible()) {
@@ -728,113 +733,9 @@ namespace lp {
         if (has_inf_int()) 
             local_gomory(5);            
 
-#if 0
         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
         //
@@ -948,7 +849,8 @@ namespace lp {
 #endif
     }
 
-    lia_move int_solver::local_gomory(unsigned num_cuts) {
+
+    lia_move int_solver::local_cut(unsigned num_cuts, std::function<lia_move(void)>& cut_fn) {
 
         struct ex { explanation m_ex; lar_term m_term; mpq m_k; bool m_is_upper; };
         vector<ex> cuts;
@@ -956,7 +858,7 @@ namespace lp {
             m_ex->clear();
             m_t.clear();
             m_k.reset();
-            auto r = gomory(*this)();
+            auto r = cut_fn();
             if (r != lia_move::cut)
                 break;
             cuts.push_back({ *m_ex, m_t, m_k, is_upper() });

src/math/lp/int_solver.h

@@ -68,7 +68,7 @@ class int_solver {
     hnf_cutter          m_hnf_cutter;
     unsigned            m_hnf_cut_period;
     unsigned_vector     m_cut_vars;        // variables that should not be selected for cuts
-
+    
     vector<equality>       m_equalities;
 public:
     int_solver(lar_solver& lp);
@@ -111,8 +111,8 @@ private:
     bool has_upper(unsigned j) const;
     unsigned row_of_basic_column(unsigned j) const;
     bool cut_indices_are_columns() const;
-    lia_move local_gomory(unsigned num_cuts);
+    lia_move local_cut(unsigned num_cuts, std::function<lia_move(void)>& cut_fn);
-    
+        
 public:
     std::ostream& display_column(std::ostream & out, unsigned j) const;
     u_dependency* column_upper_bound_constraint(unsigned j) const;

src/smt/theory_lra.cpp

@@ -1784,14 +1784,15 @@ public:
 
         expr_ref atom(m);
         t = coeffs2app(coeffs, rational::zero(), is_int);
-        if (lower_bound) {
+        if (lower_bound) 
             atom = a.mk_ge(t, a.mk_numeral(offset, is_int));
-        }
+        else 
-        else {
+            atom = a.mk_le(t, a.mk_numeral(offset, is_int));        
-            atom = a.mk_le(t, a.mk_numeral(offset, is_int));
-        }
 
         // ctx().get_rewriter()(atom);
+        // Note: it is not safe to rewrite atom because the rewriter can
+        // destroy structure, such as (div x 24) >= 0 becomes x >= 0 and the internal variable
+        // corresponding to (div x 24) is not constrained.
         TRACE("arith", tout << t << ": " << atom << "\n";
               lp().print_term(term, tout << "bound atom: ") << (lower_bound?" >= ":" <= ") << k << "\n";);
         ctx().internalize(atom, true);