lp: add best-of-N efficacy-based Gomory row selection

Add lp.gomory_efficacy_select (default false) and lp.gomory_candidate_rows
(default 3). When enabled, Gomory cut rows are no longer chosen by how close
the basic variable is to an integer. Instead, gomory_candidate_rows
integer-infeasible GMI-clean rows are sampled uniformly at random, a cut is
built for each, and the most efficacious cuts whose efficacy
(k - t(x*))/||t||_2 is at least gomory_cut_efficacy_threshold are added (up to
num_cuts of them, best first).

Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>

src/math/lp/gomory.cpp

@@ -21,6 +21,8 @@
 #include "math/lp/int_solver.h"
 #include "math/lp/lar_solver.h"
 #include "math/lp/lp_utils.h"
+#include <algorithm>
+#include <cmath>
 
 namespace lp {
     
@@ -479,11 +481,10 @@ public:
         return ret;
     }
 
-    // Returns true if the cut t >= k has efficacy at or above the configured threshold,
-    // i.e. the cut is worth keeping. efficacy = (k - t(x*)) / ||t||_2, where x* is the
-    // current LP solution. We compare efficacy^2 against threshold^2 to avoid a sqrt and
-    // keep the arithmetic exact until the final, deterministic double comparison.
-    bool gomory::cut_has_enough_efficacy(const lar_term& t, const mpq& k) {
+    // Efficacy of the cut t >= k: the distance from the current LP solution x* to the cut
+    // hyperplane, normalized by the coefficient norm, efficacy = (k - t(x*)) / ||t||_2.
+    // Returns 0 when the cut has zero norm or is not violated by x* (so it never qualifies).
+    double gomory::cut_efficacy(const lar_term& t, const mpq& k) {
         mpq val(0);     // t(x*)
         mpq norm2(0);   // ||t||_2^2
         for (lar_term::ival p : t) {
@@ -492,19 +493,46 @@ public:
             norm2 += a * a;
         }
         if (norm2.is_zero())
-            return false;
+            return 0.0;
         mpq violation = k - val; // positive, since x* violates the cut t >= k
         if (!violation.is_pos())
-            return false;
-        double thr = lia.settings().gomory_cut_efficacy_threshold();
-        // efficacy >= thr  <=>  violation^2 >= thr^2 * norm2
-        return (violation * violation).get_double() >= thr * thr * norm2.get_double();
+            return 0.0;
+        return violation.get_double() / std::sqrt(norm2.get_double());
+    }
+
+    // Returns true if the cut t >= k has efficacy at or above the configured threshold.
+    bool gomory::cut_has_enough_efficacy(const lar_term& t, const mpq& k) {
+        return cut_efficacy(t, k) >= lia.settings().gomory_cut_efficacy_threshold();
+    }
+
+    // Uniformly sample up to num_rows integer-infeasible rows that are valid Gomory cut
+    // targets, without ranking them by how close the basic variable is to an integer.
+    // This is the selection used by gomory_efficacy_select: row quality is judged later by
+    // the efficacy of the cut each row produces, not by the variable's fractionality.
+    unsigned_vector gomory::gomory_select_random_rows(unsigned num_rows) {
+        unsigned_vector eligible;
+        for (lpvar j : lra.r_basis())
+            if (lia.column_is_int_inf(j) && is_gomory_cut_target(j))
+                eligible.push_back(j);
+        unsigned_vector ret;
+        unsigned n = static_cast<unsigned>(eligible.size());
+        while (num_rows-- && n > 0) {
+            unsigned k = lia.settings().random_next() % n;
+            ret.push_back(eligible[k]);
+            eligible[k] = eligible[--n]; // swap-remove to sample without replacement
+        }
+        return ret;
     }
 
     lia_move gomory::get_gomory_cuts(unsigned num_cuts) {
         struct cut_result {lar_term t; mpq k; u_dependency *dep;};
+        struct scored_cut {lar_term t; mpq k; u_dependency *dep; double eff;};
         vector<cut_result> big_cuts;
-        unsigned_vector columns_for_cuts = gomory_select_int_infeasible_vars(num_cuts);
+        vector<scored_cut> scored_cuts; // candidates for best-of-N efficacy selection
+        const bool eff_select = lia.settings().gomory_efficacy_select();
+        unsigned_vector columns_for_cuts = eff_select
+            ? gomory_select_random_rows(lia.settings().gomory_candidate_rows())
+            : gomory_select_int_infeasible_vars(num_cuts);
         bool has_small_cut = false;
 
         // define inline helper functions
@@ -542,12 +570,21 @@ public:
             else if (cc.m_polarity == row_polarity::MIN)
                 lra.update_column_type_and_bound(j, lp::lconstraint_kind::GE, ceil(lra.get_column_value(j).x), add_deps(cc.m_dep, row, j));
 
+            // Best-of-N selection: collect every candidate cut that clears the efficacy
+            // threshold, then (after the loop) add the most efficacious ones. The row was
+            // picked at random rather than by the basic variable's fractionality, so the
+            // cut's efficacy is what decides whether and how strongly it is preferred.
+            if (eff_select) {
+                double e = cut_efficacy(cc.m_t, cc.m_k);
+                if (e >= lia.settings().gomory_cut_efficacy_threshold())
+                    scored_cuts.push_back({cc.m_t, cc.m_k, cc.m_dep, e});
+                continue;
+            }
+
             // Option A: optionally discard cuts whose efficacy (the distance from the LP
             // solution to the cut hyperplane, normalized by the coefficient norm) is too small.
             // The cut is m_t >= m_k and the current LP solution violates it, so the
             // violation m_k - m_t(x*) is positive. efficacy = violation / ||m_t||_2.
-            // To stay rational and deterministic we compare efficacy^2 against threshold^2,
-            // i.e. violation^2 >= threshold^2 * ||m_t||_2^2.
             if (lia.settings().gomory_cut_efficacy_filter() && !cut_has_enough_efficacy(cc.m_t, cc.m_k))
                 continue;
 
@@ -561,6 +598,26 @@ public:
                 return lia_move::cancelled;
         }
 
+        // best-of-N: add up to num_cuts of the most efficacious collected candidates
+        if (eff_select) {
+            std::stable_sort(scored_cuts.begin(), scored_cuts.end(),
+                             [](scored_cut const& a, scored_cut const& b) { return a.eff > b.eff; });
+            unsigned added = 0;
+            for (auto const& c : scored_cuts) {
+                if (added >= num_cuts)
+                    break;
+                ++added;
+                if (!is_small_cut(c.t)) {
+                    big_cuts.push_back({c.t, c.k, c.dep});
+                    continue;
+                }
+                has_small_cut = true;
+                add_cut(c.t, c.k, c.dep);
+                if (lia.settings().get_cancel_flag())
+                    return lia_move::cancelled;
+            }
+        }
+
         if (big_cuts.size()) {
             lra.push();        
             for (auto const& cut : big_cuts) 

src/math/lp/gomory.h

@@ -28,8 +28,10 @@ namespace lp {
         class int_solver& lia;
         class lar_solver& lra;
         unsigned_vector gomory_select_int_infeasible_vars(unsigned num_cuts);
+        unsigned_vector gomory_select_random_rows(unsigned num_rows);
         bool is_gomory_cut_target(lpvar j); 
         bool cut_has_enough_efficacy(const lar_term& t, const mpq& k);
+        double cut_efficacy(const lar_term& t, const mpq& k);
         u_dependency* add_deps(u_dependency*, const row_strip<mpq>&, lpvar);
     public:
         lia_move get_gomory_cuts(unsigned num_cuts);

src/math/lp/lp_params_helper.pyg

@@ -16,6 +16,8 @@ def_module_params(module_name='lp',
                           ('int_hammer_period', UINT, 4, 'period (in final_check calls) for the integer cut/cube heuristics (find_cube, hnf, gomory); a smaller value calls them more often'),
                           ('random_hammers', BOOL, True, 'draw the periodic integer heuristic gates (find_cube, lcube, hnf, gomory, dio) at random with the same 1/period rate instead of a deterministic every-k-th-call modulus'),
                           ('gomory_cut_efficacy_filter', BOOL, False, 'discard a generated Gomory cut whose efficacy (distance from the LP solution to the cut hyperplane, normalized by the cut coefficient norm) is below gomory_cut_efficacy_threshold'),
-                          ('gomory_cut_efficacy_threshold', DOUBLE, 0.01, 'minimal efficacy (normalized violation) required to keep a Gomory cut when gomory_cut_efficacy_filter is enabled'),
+                          ('gomory_cut_efficacy_threshold', DOUBLE, 0.01, 'minimal efficacy (normalized violation) required to keep a Gomory cut when gomory_cut_efficacy_filter or gomory_efficacy_select is enabled'),
+                          ('gomory_efficacy_select', BOOL, False, 'select Gomory cut rows by cut efficacy instead of by how close the basic variable is to an integer: pick gomory_candidate_rows integer-infeasible rows at random, build their cuts, and add the most efficacious ones whose efficacy is at least gomory_cut_efficacy_threshold'),
+                          ('gomory_candidate_rows', UINT, 3, 'number of integer-infeasible rows sampled at random to build candidate Gomory cuts from when gomory_efficacy_select is enabled'),
                          ))
                          

src/math/lp/lp_settings.cpp

@@ -50,6 +50,8 @@ void lp::lp_settings::updt_params(params_ref const& _p) {
     m_lcube_flips = lp_p.lcube_flips();
     m_gomory_cut_efficacy_filter = lp_p.gomory_cut_efficacy_filter();
     m_gomory_cut_efficacy_threshold = lp_p.gomory_cut_efficacy_threshold();
+    m_gomory_efficacy_select = lp_p.gomory_efficacy_select();
+    m_gomory_candidate_rows = lp_p.gomory_candidate_rows();
     unsigned hammer_period = lp_p.int_hammer_period();
     SASSERT(hammer_period != 0);
     m_int_find_cube_period = hammer_period;

src/math/lp/lp_settings.h

@@ -272,11 +272,15 @@ private:
     unsigned         m_lcube_flips = 16;
     bool             m_gomory_cut_efficacy_filter = false;
     double           m_gomory_cut_efficacy_threshold = 0.01;
+    bool             m_gomory_efficacy_select = false;
+    unsigned         m_gomory_candidate_rows = 3;
 public:
     bool lcube() const { return m_lcube; }
     unsigned lcube_flips() const { return m_lcube_flips; }
     bool gomory_cut_efficacy_filter() const { return m_gomory_cut_efficacy_filter; }
     double gomory_cut_efficacy_threshold() const { return m_gomory_cut_efficacy_threshold; }
+    bool gomory_efficacy_select() const { return m_gomory_efficacy_select; }
+    unsigned gomory_candidate_rows() const { return m_gomory_candidate_rows; }
     unsigned dio_calls_period() const { return m_dio_calls_period; }
     unsigned & dio_calls_period() { return m_dio_calls_period; }
     unsigned dio_calls_period_decrease() const { return m_dio_calls_period_decrease; }