a simple version of choosing a column for gomory cut

src/math/lp/gomory.cpp

@@ -410,69 +410,11 @@ lia_move gomory::cut(lar_term & t, mpq & k, explanation* ex, unsigned basic_inf_
     return cc.cut();
 }
     
-bool gomory::is_gomory_cut_target(const row_strip<mpq>& row) {
+lia_move gomory::get_cut(lpvar j) {
-    // All non base variables must be at their bounds and assigned to rationals (that is, infinitesimals are not allowed).
-    unsigned j;
-    for (const auto & p : row) {
-        j = p.var();
-        if (!lia.is_base(j) && (!lia.at_bound(j) || !is_zero(lia.get_value(j).y))) {
-            TRACE("gomory_cut", tout << "row is not gomory cut target:\n";
-                  lia.display_column(tout, j);
-                  tout << "infinitesimal: " << !is_zero(lia.get_value(j).y) << "\n";);
-            return false;
-        }
-    }
-    return true;
-}
-
-int gomory::find_basic_var() {
-    unsigned n = 0;
-    int result = -1;
-    unsigned min_row_size = UINT_MAX;
-
-#if 1
-    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)));
-        // Prefer smaller row size
-        if (min_row_size == UINT_MAX || 
-            2*row.size() < min_row_size || 
-            (4*row.size() < 5*min_row_size && lia.random() % (++n) == 0)) {
-            result = j;
-            n = 1;
-            min_row_size = std::min(min_row_size, row.size());
-        }
-    }
-    return result;
-}
-    
-lia_move gomory::operator()() {
-    int j = find_basic_var();
-    if (j == -1)
-        return lia_move::undef;
     unsigned r = lia.row_of_basic_column(j);
     const row_strip<mpq>& row = lra.get_row(r);
     SASSERT(lra.row_is_correct(r));
-    SASSERT(is_gomory_cut_target(row));
+    SASSERT(lia.is_gomory_cut_target(j));
     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/gomory.h

@@ -27,11 +27,9 @@ namespace lp {
     class gomory {
         class int_solver& lia;
         class lar_solver& lra;
-        int find_basic_var();
-        bool is_gomory_cut_target(const row_strip<mpq>& row);
         lia_move cut(lar_term & t, mpq & k, explanation* ex, unsigned basic_inf_int_j, const row_strip<mpq>& row);
     public:
         gomory(int_solver& lia);
-        lia_move operator()();
+        lia_move get_cut(lpvar j);
     };
 }

src/math/lp/int_branch.cpp

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

src/math/lp/int_solver.cpp

@@ -198,15 +198,9 @@ namespace lp {
         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)(); };
+        std::function<lia_move(lpvar)> gomory_fn = [&](lpvar j) { return gomory(*this).get_cut(j); };
-        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_cut(2, gomory_fn);
 
-#endif
-        m_cut_vars.reset();
         if (r == lia_move::undef) r = int_branch(*this)();
         if (settings().get_cancel_flag()) r = lia_move::undef;        
         return r;
@@ -635,7 +629,7 @@ namespace lp {
     }
 
 
-    int int_solver::select_int_infeasible_var(bool check_bounded) {
+    int int_solver::select_int_infeasible_var() {
         int r_small_box = -1;
         int r_small_value = -1;
         int r_any_value = -1;
@@ -648,18 +642,6 @@ namespace lp {
         lar_core_solver & lcs = lra.m_mpq_lar_core_solver;
         unsigned prev_usage = 0;
 
-        auto check_bounded_fn = [&](unsigned j) {
-            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;
@@ -670,9 +652,7 @@ namespace lp {
         for (unsigned j : lra.r_basis()) {
             if (!column_is_int_inf(j))
                 continue;
-            if (!check_bounded_fn(j))
+             if (m_cut_vars.contains(j))
-                continue;
-            if (m_cut_vars.contains(j))
                 continue;
 
             SASSERT(!is_fixed(j));
@@ -849,19 +829,71 @@ namespace lp {
 
 #endif
     }
+    // return the minimal distance from the column value to an integer
+    mpq get_gomory_score(const int_solver& lia, lpvar j) {
+        const mpq& val = lia.get_value(j).x;
+        auto l = val - floor(val);
+        if (l <= mpq(1, 2))
+            return l;
+        return mpq(1) - l;
+    }
 
+    unsigned_vector int_solver::gomory_select_int_infeasible_vars(unsigned num_cuts) {
+        SASSERT(m_cut_vars.size() == 0&& num_cuts >= 0);
         
-    lia_move int_solver::local_cut(unsigned num_cuts, std::function<lia_move(void)>& cut_fn) {
+        std::list<lpvar> sorted_vars;
+        std::unordered_map<lpvar, mpq> score;
+        for (lpvar j : lra.r_basis()) {
+            if (!column_is_int_inf(j) || !is_gomory_cut_target(j))
+                continue;
+            SASSERT(!is_fixed(j));            
+            sorted_vars.push_back(j);
+            score[j] = get_gomory_score(*this, j);
+        }
+        // prefer the columns with the values close to integers
+        sorted_vars.sort([&](lpvar j, lpvar k) {
+            auto diff = score[j] - score[k];
+            if (diff.is_neg())
+                return true;
+            if (diff.is_pos())
+                return false;
+            return lra.usage_in_terms(j) > lra.usage_in_terms(k);
+        });
+        unsigned_vector ret;
+        unsigned n = static_cast<unsigned>(sorted_vars.size());
+
+        while (num_cuts-- && n > 0) {
+            unsigned k = random() % n;
+           
+            double k_ratio = k / (double) n;
+            k_ratio *= k_ratio*k_ratio;  // square k_ratio to make it smaller
+            k = static_cast<unsigned>(std::floor(k_ratio * n));
+            // these operations move k to the beginning of the indices range
+            SASSERT(0 <= k && k < n);
+            auto it = sorted_vars.begin();
+            while(k--) it++;
+
+            ret.push_back(*it);
+            sorted_vars.erase(it);
+            n--;            
+        }
+        return ret;
+    }
+    
+    lia_move int_solver::local_cut(unsigned num_cuts, std::function<lia_move(lpvar)>& cut_fn) {
         
         struct ex { explanation m_ex; lar_term m_term; mpq m_k; bool m_is_upper; };
+        unsigned_vector columns_for_cuts = gomory_select_int_infeasible_vars(num_cuts);
+
         vector<ex> cuts;
-        for (unsigned i = 0; i < num_cuts && has_inf_int(); ++i) {
+        
+        for (unsigned j : columns_for_cuts) {
             m_ex->clear();
             m_t.clear();
             m_k.reset();
-            auto r = cut_fn();
+            auto r = cut_fn(j);
             if (r != lia_move::cut)
-                break;
+                continue;
             cuts.push_back({ *m_ex, m_t, m_k, is_upper() });
             if (settings().get_cancel_flag())
                 return lia_move::undef;
@@ -938,6 +970,22 @@ namespace lp {
         return lia_move::undef;
     }
 
+    bool int_solver::is_gomory_cut_target(lpvar k) {
+        SASSERT(is_base(k));
+        // All non base variables must be at their bounds and assigned to rationals (that is, infinitesimals are not allowed).
+        const row_strip<mpq>& row = lra.get_row(row_of_basic_column(k));
+        unsigned j;
+        for (const auto & p : row) {
+            j = p.var();
+            if ( k != j && (!at_bound(j) || !is_zero(get_value(j).y))) {
+                TRACE("gomory_cut", tout << "row is not gomory cut target:\n";
+                      display_column(tout, j);
+                      tout << "infinitesimal: " << !is_zero(get_value(j).y) << "\n";);
+                return false;
+            }
+        }
+        return true;
+    }
 
 
 }

src/math/lp/int_solver.h

@@ -111,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_cut(unsigned num_cuts, std::function<lia_move(void)>& cut_fn);
+    lia_move local_cut(unsigned num_cuts, std::function<lia_move(lpvar)>& cut_fn);
         
 public:
     std::ostream& display_column(std::ostream & out, unsigned j) const;
@@ -132,7 +132,8 @@ public:
     bool all_columns_are_bounded() const;
     lia_move hnf_cut();
 
-    int select_int_infeasible_var(bool check_bounded);
+    int select_int_infeasible_var();
-
+    unsigned_vector gomory_select_int_infeasible_vars(unsigned num_cuts);
+    bool is_gomory_cut_target(lpvar);
 };
 }

src/math/lp/lar_solver.h

@@ -532,6 +532,11 @@ public:
     inline const impq& get_lower_bound(column_index j) const {
         return m_mpq_lar_core_solver.m_r_solver.m_lower_bounds[j];
     }
+
+    inline mpq bound_span_x(column_index j) const {
+        return m_mpq_lar_core_solver.m_r_solver.m_upper_bounds[j].x - m_mpq_lar_core_solver.m_r_solver.m_lower_bounds[j].x;
+    }
+    
     bool has_lower_bound(var_index var, u_dependency*& ci, mpq& value, bool& is_strict) const;
     bool has_upper_bound(var_index var, u_dependency*& ci, mpq& value, bool& is_strict) const;
     bool has_value(var_index var, mpq& value) const;