Gomory cut / branch and bound improvements

Improve fairness of cut generation by switching to find_infeasible_int_var with cascading priorities, allow stronger cuts by inlining terms.

src/math/lp/gomory.cpp

@@ -31,7 +31,7 @@ class create_cut {
     explanation*          m_ex; // the conflict explanation
     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;
-    const int_solver&     lia;
+    int_solver&           lia;
     mpq                   m_lcm_den;
     mpq                   m_f;
     mpq                   m_one_minus_f;
@@ -133,26 +133,50 @@ class create_cut {
         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
-        auto pol = m_t.coeffs_as_vector();
+        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;);
-            unsigned v = pol[0].second;
+            auto const& [a, v] = pol[0];
             lp_assert(is_int(v));
-            const mpq& a = pol[0].first;
             if (a.is_pos()) { // we have av >= k
-                m_k /= a;
+                divd(m_k, a);
-                if (!m_k.is_int())
-                    m_k = ceil(m_k);
                 m_t.add_monomial(mpq(1), v);
             }
             else {
-                m_k /= -a;
+                // av >= k
-                if (!m_k.is_int())
+                // a/-a*v >= k / - a
-                    m_k = ceil(m_k);
+                // -v >= k / - a
+                // -v >= ceil(k / -a)
+                divd(m_k, -a);
                 m_t.add_monomial(-mpq(1), v);
             }
         }
@@ -162,17 +186,48 @@ class create_cut {
             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 & pi: pol) {
+                for (auto & [c,v]: pol) {
-                    pi.first *= m_lcm_den;
+                    c *= m_lcm_den;
-                    SASSERT(!is_int(pi.second) || pi.first.is_int());
+                    SASSERT(!is_int(v) || c.is_int());
                 }
                 m_k *= m_lcm_den;
             }
-            for (const auto & pi: pol)
+
-                m_t.add_monomial(pi.first, pi.second);
+            // 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 {
@@ -191,10 +246,9 @@ class create_cut {
     
     template <typename T>
     void dump_coeff(std::ostream & out, const T& c) const {
-        out << "( * ";
+        out << "(* ";
         dump_coeff_val(out, c.coeff());
-        auto t = lia.lra.column2tv(c.column());
+        out << " " << var_name(c.column().index()) << ")";
-        out << " " << var_name(t.id()) << ")";
     }
     
     std::ostream& dump_row_coefficients(std::ostream & out) const {
@@ -208,7 +262,7 @@ class create_cut {
     
     void dump_the_row(std::ostream& out) const {
         out << "; the row, excluding fixed vars\n";
-        out << "(assert ( = ( +";
+        out << "(assert (= (+";
         dump_row_coefficients(out) << ") 0))\n";
     }
 
@@ -259,12 +313,12 @@ class create_cut {
         return dump_term_coefficients(out << "(+ ") << ")";
     }
     
-    std::ostream& dump_term_le_k(std::ostream & out) const {
+    std::ostream& dump_term_ge_k(std::ostream & out) const {
-        return dump_term_sum(out << "(<= ") << " " << m_k << ")";
+        return dump_term_sum(out << "(>= ") << " " << m_k << ")";
     }
 
     void dump_the_cut_assert(std::ostream & out) const {
-        dump_term_le_k(out << "(assert (not ") << "))\n";
+        dump_term_ge_k(out << "(assert (not ") << "))\n";
     }
 
     void dump_cut_and_constraints_as_smt_lemma(std::ostream& out) const {
@@ -310,7 +364,6 @@ public:
             unsigned j = p.var();
             if (j == m_inf_col) {
                 lp_assert(p.coeff() == one_of_type<mpq>());
-                TRACE("gomory_cut_detail", tout << "seeing basic var\n";);
                 continue;
             }
 
@@ -341,13 +394,21 @@ public:
             return report_conflict_from_gomory_cut();
         if (some_int_columns)
             adjust_term_and_k_for_some_ints_case_gomory();
-        TRACE("gomory_cut_detail", dump_cut_and_constraints_as_smt_lemma(tout););
+        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", 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);
+              lia.lra.display(tout));
+        lia.settings().stats().m_gomory_cuts++;
         return lia_move::cut;
     }
 
-    create_cut(lar_term & t, mpq & k, explanation* ex, unsigned basic_inf_int_j, const row_strip<mpq>& row, const int_solver& lia) :
+    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),
         m_ex(ex),
@@ -385,22 +446,26 @@ int gomory::find_basic_var() {
     int result = -1;
     unsigned min_row_size = UINT_MAX;
 
-#if 0
+#if 1
-    result = lia.select_int_infeasible_var();
+    result = lia.select_int_infeasible_var(true);
 
     if (result == -1)
         return result;
 
+    TRACE("gomory_cut", tout << "row: " << result << "\n");
     const row_strip<mpq>& row = lra.get_row(lia.row_of_basic_column(result));
     if (is_gomory_cut_target(row)) 
         return result;
     result = -1;
+
+    UNREACHABLE();
 #endif
 
     for (unsigned j : lra.r_basis()) {
         if (!lia.column_is_int_inf(j))
             continue;
         const row_strip<mpq>& row = lra.get_row(lia.row_of_basic_column(j));
+        TRACE("gomory_cut", tout << "try j" << j << "\n");
         if (!is_gomory_cut_target(row)) 
             continue;
         IF_VERBOSE(20, lia.display_row_info(verbose_stream(), lia.row_of_basic_column(j)));
@@ -417,7 +482,6 @@ int gomory::find_basic_var() {
 }
     
 lia_move gomory::operator()() {
-    lra.move_non_basic_columns_to_bounds();
     int j = find_basic_var();
     if (j == -1)
         return lia_move::undef;
@@ -426,6 +490,7 @@ lia_move gomory::operator()() {
     SASSERT(lra.row_is_correct(r));
     SASSERT(is_gomory_cut_target(row));
     lia.m_upper = false;
+    lia.m_cut_vars.push_back(j);
     return cut(lia.m_t, lia.m_k, lia.m_ex, j, row);
 }
 

src/math/lp/int_branch.cpp

@@ -52,8 +52,8 @@ lia_move int_branch::create_branch_on_column(int j) {
 
 int int_branch::find_inf_int_base_column() {
 
-#if 0
+#if 1
-    return lia.select_int_infeasible_var();
+    return lia.select_int_infeasible_var(false);
 #endif
 
     int result = -1;

src/math/lp/int_solver.cpp

@@ -180,6 +180,8 @@ namespace lp {
         m_ex = e;
         m_ex->clear();
         m_upper = false;
+        m_cut_vars.reset();
+        
         lia_move r = lia_move::undef;
 
         if (m_gcd.should_apply())
@@ -193,12 +195,15 @@ namespace lp {
         ++m_number_of_calls;
         if (r == lia_move::undef && m_patcher.should_apply()) r = m_patcher();
         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();
-#if 1
+        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();
+        if (r == lia_move::undef && should_gomory_cut()) r = local_gomory(2);
 #endif
+        m_cut_vars.reset();
         if (r == lia_move::undef) r = int_branch(*this)();
         return r;
     }
@@ -626,71 +631,85 @@ namespace lp {
     }
 
 
-    int int_solver::select_int_infeasible_var() {
+    int int_solver::select_int_infeasible_var(bool check_bounded) {
-        int result = -1;
+        int r_small_box = -1;
+        int r_small_value = -1;
+        int r_any_value = -1;
+        unsigned n_small_box = 1;
+        unsigned n_small_value = 1;
+        unsigned n_any_value = 1;
         mpq range;
         mpq new_range;
         mpq small_value(1024);
-        unsigned n = 0;
         lar_core_solver & lcs = lra.m_mpq_lar_core_solver;
-        unsigned prev_usage = 0; // to quiet down the compile
+        unsigned prev_usage = 0;
 
-        enum state { small_box, is_small_value,  any_value, not_found };
+        auto check_bounded_fn = [&](unsigned j) {
-        state st = not_found;
+            if (!check_bounded)
+                return true;
+            auto const& row = lra.get_row(row_of_basic_column(j));
+            for (const auto & p : row) {
+                unsigned j = p.var();
+                if (!is_base(j) && (!at_bound(j) || !is_zero(get_value(j).y)))
+                    return false;
+            }
+            return true;
+        };
 
+        auto add_column = [&](bool improved, int& result, unsigned& n, unsigned j) {
+            if (result == -1)
+                result = j;
+            else if (improved && ((random() % (++n)) == 0))
+                result = j;            
+        };
+        
         for (unsigned j : lra.r_basis()) {
             if (!column_is_int_inf(j))
                 continue;
+            if (!check_bounded_fn(j))
+                continue;
+            if (m_cut_vars.contains(j))
+                continue;
+
+            SASSERT(!is_fixed(j));
+
             unsigned usage = lra.usage_in_terms(j);
             if (is_boxed(j) && (new_range = lcs.m_r_upper_bounds()[j].x - lcs.m_r_lower_bounds()[j].x - rational(2*usage)) <= small_value) {
-                SASSERT(!is_fixed(j));
+
-                if (st != small_box) {
+                bool improved = new_range <= range || r_small_box == -1;
-                    n = 0;
+                if (improved)
-                    st = small_box;
-                }
-                if (n == 0 || new_range < range) {
-                    result = j;
                     range = new_range;
-                    n = 1;
+                add_column(improved, r_small_box, n_small_box, j);
-                }
-                else if (new_range == range && (random() % (++n) == 0)) {
-                    result = j;
-                }
                 continue;
             }
-            if (st == small_box)
-                continue;
             impq const& value = get_value(j);
             if (abs(value.x) < small_value ||
                 (has_upper(j) && small_value > upper_bound(j).x - value.x) ||
                 (has_lower(j) && small_value > value.x - lower_bound(j).x)) {
-                if (st != is_small_value) {
+                TRACE("gomory_cut", tout << "small j" << j << "\n");
-                    n = 0;
+                add_column(true, r_small_value, n_small_value, j);
-                    st = is_small_value;
-                }
-                if (random() % (++n) == 0)
-                    result = j;
-            }
-            if (st == is_small_value)
                 continue;
-            SASSERT(st == not_found || st == any_value);
+            }
-            st = any_value;
+            TRACE("gomory_cut", tout << "any j" << j << "\n");
-            if (n == 0 || usage > prev_usage) {
+            add_column(usage >= prev_usage, r_any_value, n_any_value, j);
-                result = j;
+            if (usage > prev_usage) 
                 prev_usage = usage;
-                n = 1;
-            } 
-            else if (usage > 0 && usage == prev_usage && (random() % (++n) == 0)) 
-                result = j;
         }
-    
+
-        return result;
+        if (r_small_box != -1 && (random() % 3 != 0))
+            return r_small_box;
+        if (r_small_value != -1 && (random() % 3) != 0)
+            return r_small_value;
+        if (r_any_value != -1)
+            return r_any_value;
+        if (r_small_box != -1)
+            return r_small_box;
+        return r_small_value;
     }
 
     void int_solver::simplify(std::function<bool(unsigned)>& is_root) {
 
         return;
-#if 1
 
         // in-processing simplification can go here, such as bounds improvements.
 
@@ -701,17 +720,13 @@ namespace lp {
         }
 
         
-#endif
-
-#if 1
         lp::explanation exp;
         m_ex = &exp;
         m_t.clear();
         m_k.reset();
 
         if (has_inf_int()) 
-            local_gomory();            
+            local_gomory(5);            
-#endif
 
 #if 0
         stopwatch sw;
@@ -933,35 +948,85 @@ namespace lp {
 #endif
     }
 
-    lia_move int_solver::local_gomory() {
+    lia_move int_solver::local_gomory(unsigned num_cuts) {
-        for (unsigned i = 0; i < 2 && has_inf_int() && !settings().get_cancel_flag(); ++i) {
+
+        struct ex { explanation m_ex; lar_term m_term; mpq m_k; bool m_is_upper; };
+        vector<ex> cuts;
+        for (unsigned i = 0; i < num_cuts && has_inf_int() && !settings().get_cancel_flag(); ++i) {
             m_ex->clear();
             m_t.clear();
             m_k.reset();
             auto r = gomory(*this)();
-            IF_VERBOSE(3, verbose_stream() << i << " " << r << "\n");
+            if (r != lia_move::cut)
-            if (r != lia_move::cut) 
+                break;
-                return r;
+            cuts.push_back({ *m_ex, m_t, m_k, is_upper() });
-            u_dependency* dep = nullptr;
-            for (auto c : *m_ex) 
-                dep = lra.join_deps(lra.dep_manager().mk_leaf(c.ci()), dep);
-            lp::lpvar term_index = lra.add_term(get_term().coeffs_as_vector(), UINT_MAX);
-            term_index = lra.map_term_index_to_column_index(term_index);
-            lra.update_column_type_and_bound(term_index, is_upper() ? lp::lconstraint_kind::LE : lp::lconstraint_kind::GE, get_offset(), dep);
-            lra.find_feasible_solution();
-            if (!lra.is_feasible()) {
-                lra.get_infeasibility_explanation(*m_ex);
-                return lia_move::conflict;
-            }
-            //r = m_patcher();
-            //if (r != lia_move::undef)
-            //    return r;
         }
+        m_cut_vars.reset();
+
+        auto is_small_cut = [&](ex const& cut) {
+            return all_of(cut.m_term, [&](auto ci) { return ci.coeff().is_small(); });
+        };
+
+        auto add_cut = [&](ex const& cut) {
+            u_dependency* dep = nullptr;
+            for (auto c : cut.m_ex) 
+                dep = lra.join_deps(lra.dep_manager().mk_leaf(c.ci()), dep);
+            lp::lpvar term_index = lra.add_term(cut.m_term.coeffs_as_vector(), UINT_MAX);
+            term_index = lra.map_term_index_to_column_index(term_index);
+            lra.update_column_type_and_bound(term_index,
+                                             cut.m_is_upper ? lp::lconstraint_kind::LE : lp::lconstraint_kind::GE,
+                                             cut.m_k, dep);
+        };
+
+        auto _check_feasible = [&](void) {
+            auto st = lra.find_feasible_solution();
+            if (!lra.is_feasible()) {                
+                lra.get_infeasibility_explanation(*m_ex);
+                return false;
+            }
+            return true;
+        };
+
+        bool has_small = false, has_large = false;
+
+        for (auto const& cut : cuts) {
+            if (!is_small_cut(cut)) {
+                has_large = true;
+                continue;
+            }
+            has_small = true;
+            add_cut(cut);
+        }
+
+        if (has_large) {
+            lra.push();
+        
+            for (auto const& cut : cuts) 
+                if (!is_small_cut(cut))
+                    add_cut(cut);
+
+            bool feas = _check_feasible();
+            lra.pop(1);
+
+            if (!feas)
+                return lia_move::conflict;
+            
+        }
+
+        if (!_check_feasible())
+            return lia_move::conflict;
+        
+                
         m_ex->clear();
         m_t.clear();
         m_k.reset();
         if (!has_inf_int())
             return lia_move::sat;
+
+        if (has_small || has_large)
+            return lia_move::continue_with_check;
+        
+        lra.move_non_basic_columns_to_bounds();
         return lia_move::undef;
     }
 

src/math/lp/int_solver.h

@@ -63,10 +63,11 @@ class int_solver {
     unsigned            m_number_of_calls;
     lar_term            m_t;               // the term to return in the cut
     mpq                 m_k;               // the right side of the cut
+    bool                m_upper;           // cut is an upper bound
     explanation         *m_ex;             // the conflict explanation
-    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;
+    unsigned_vector     m_cut_vars;        // variables that should not be selected for cuts
 
     vector<equality>       m_equalities;
 public:
@@ -110,7 +111,7 @@ 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();
+    lia_move local_gomory(unsigned num_cuts);
     
 public:
     std::ostream& display_column(std::ostream & out, unsigned j) const;
@@ -131,7 +132,7 @@ public:
     bool all_columns_are_bounded() const;
     lia_move hnf_cut();
 
-    int select_int_infeasible_var();
+    int select_int_infeasible_var(bool check_bounded);
 
 };
 }

src/smt/theory_lra.cpp

@@ -1678,7 +1678,6 @@ public:
                 return FC_CONTINUE;
             }
             
-
             for (expr* e : m_not_handled) {
                 if (!ctx().is_relevant(e))
                     continue;