Improve duplicate cut detection

cpp/src/cuts/cuts.cpp

@@ -14,13 +14,15 @@
 #include <utilities/logger.hpp>
 #include <utilities/macros.cuh>
 
+#include <algorithm>
 #include <array>
 #include <cstdint>
 #include <cstdio>
 #include <cstdlib>
 #include <limits>
 #include <stdexcept>
 #include <tuple>
+#include <unordered_map>
 #include <unordered_set>
 
 #include <linear_algebra/dense_matrix.hpp>
@@ -1227,55 +1229,73 @@ f_t cut_pool_t<i_t, f_t>::cut_orthogonality(i_t i, i_t j)
 template <typename i_t, typename f_t>
 void cut_pool_t<i_t, f_t>::check_for_duplicate_cuts()
 {
-  // Algorithm from Finding Duplicate Rows in a Linear Programming Model
-  // by J. A. Tomlin and J.S. Welch
-  // Operations Research Letters Volume 5, Number 1, June 1986
-  std::vector<f_t> divisors(cut_storage_.m, 0.0);
-  std::vector<i_t> sets(cut_storage_.m, 0);
+  const i_t m = cut_storage_.m;
 
+  constexpr f_t duplicate_tolerance = 1e-10;
+  std::vector<f_t> divisors(m, 0.0);
+  std::vector<i_t> sets(m, 0);
   csc_matrix_t<i_t, f_t> cut_storage_csc(0, 0, 1);
   cut_storage_.to_compressed_col(cut_storage_csc);
-  i_t n = cut_storage_csc.n;
-  i_t m = cut_storage_csc.m;
-
+  const i_t n        = cut_storage_csc.n;
   const i_t sentinel = std::numeric_limits<i_t>::max();
 
+  // Algorithm from Finding Duplicate Rows in a Linear Programming Model
+  // by J. A. Tomlin and J.S. Welch
+  // Operations Research Letters Volume 5, Number 1, June 1986.
+  //
+  // Preserve the legacy partition refinement and row-ordered deletion semantics, but index
+  // entries by their current set. This avoids scanning unrelated later entries without changing
+  // the first matching partner or the resulting removal mask.
   i_t new_set                        = 1;
-  i_t remaining_potential_duplicates = cut_storage_.m;
+  i_t remaining_potential_duplicates = m;
   for (i_t j = 0; j < n; j++) {
     i_t r0        = -1;
     i_t new_rows  = 0;
     i_t new_set_0 = new_set;
     new_set++;
     const i_t col_start = cut_storage_csc.col_start[j];
     const i_t col_end   = cut_storage_csc.col_start[j + 1];
+
+    std::unordered_map<i_t, std::vector<i_t>> positions_by_set;
+    positions_by_set.reserve(col_end - col_start);
+    for (i_t p = col_start; p < col_end; p++) {
+      const i_t set = sets[cut_storage_csc.i[p]];
+      if (set > 0 && set < new_set_0) { positions_by_set[set].push_back(p); }
+    }
+
     for (i_t p = col_start; p < col_end; p++) {
       const i_t r    = cut_storage_csc.i[p];
       const f_t a_rj = cut_storage_csc.x[p];
       const f_t f_r  = divisors[r];
       if (sets[r] == 0) {
-        r0          = r;  // To enable use to find this new set later
+        r0          = r;
         sets[r]     = new_set_0;
         divisors[r] = a_rj;
         new_rows++;
       } else if (sets[r] < new_set_0) {
-        // Look over indices a_ij with i > r
-        for (i_t q = p + 1; q < col_end; q++) {
-          const i_t i    = cut_storage_csc.i[q];
-          const f_t a_ij = cut_storage_csc.x[q];
-          if (sets[i] == sets[r]) {
-            // These two rows are currently in the same set
-            // Check to see if the coefficients still match
-            const f_t f_i     = divisors[i];
-            const f_t val     = (a_rj / f_r) * (f_i / a_ij);
-            const f_t epsilon = 1e-10;
-            if ((val >= 1.0 - epsilon && val <= 1.0 + epsilon)) {
+        const i_t old_set = sets[r];
+        bool matched      = false;
+        const auto set_positions = positions_by_set.find(old_set);
+        if (set_positions != positions_by_set.end()) {
+          auto q_position = std::upper_bound(
+            set_positions->second.begin(), set_positions->second.end(), p);
+          for (; q_position != set_positions->second.end(); ++q_position) {
+            const i_t q    = *q_position;
+            const i_t i    = cut_storage_csc.i[q];
+            const f_t a_ij = cut_storage_csc.x[q];
+            if (sets[i] != old_set) { continue; }
+            const f_t f_i = divisors[i];
+            const f_t val = (a_rj / f_r) * (f_i / a_ij);
+            if (val >= 1.0 - duplicate_tolerance &&
+                val <= 1.0 + duplicate_tolerance) {
               sets[r] = new_set;
               sets[i] = new_set;
+              matched = true;
+              break;
             }
           }
         }
-        if (sets[r] >= new_set_0) {  // This is only true if a match was found inside the above loop
+        if (matched) {
           new_set++;
         } else {
           sets[r] = sentinel;
@@ -1292,50 +1312,42 @@ void cut_pool_t<i_t, f_t>::check_for_duplicate_cuts()
     }
   }
 
-  // The cuts are stored in the form: sum_j d_ij x_j >= rhs_i
-  // We now look for cuts that are duplicates of each other and remove them
+  std::unordered_map<i_t, std::vector<i_t>> rows_by_set;
+  rows_by_set.reserve(m);
+  for (i_t r = 0; r < m; r++) {
+    if (sets[r] > 0 && sets[r] < sentinel) { rows_by_set[sets[r]].push_back(r); }
+  }
+
   std::vector<i_t> cuts_to_remove(m, 0);
-  i_t num_cuts_to_remove = 0;
   for (i_t r = 0; r < m; r++) {
     const i_t set_r = sets[r];
-    if (set_r > 0 && set_r < sentinel && cuts_to_remove[r] == 0) {
-      // This cut has a duplicate
-      for (i_t i = r + 1; i < m; i++) {
-        if (sets[i] == set_r) {
-          const f_t f_r     = divisors[r];
-          const f_t f_i     = divisors[i];
-          const f_t theta_r = rhs_storage_[r] / f_r;
-          const f_t theta_i = rhs_storage_[i] / f_i;
-          if (f_r > 0 && f_i > 0) {
-            // We have sum_j d_rj / f_r x_j >= rhs_r / f_r = theta_r
-            //    and  sum_j d_ij / f_i x_j >= rhs_i / f_i = theta_i
-            if (theta_r <= theta_i) {
-              // Cut i is either the same or stronger than cut r
-              if (cuts_to_remove[r] == 0) { num_cuts_to_remove++; }
-              cuts_to_remove[r] = 1;  // Remove row r
-            } else {
-              // theta_r > theta_i, so cut r is stricly stronger than cut i
-              if (cuts_to_remove[i] == 0) { num_cuts_to_remove++; }
-              cuts_to_remove[i] = 1;  // Remove row i
-            }
-          } else if (f_r < 0 && f_i < 0) {
-            // We have sum_j d_rj / f_r x_j <= rhs_r / f_r = theta_r
-            //    and  sum_j d_ij / f_i x_j <= rhs_i / f_i = theta_i
-            if (theta_r >= theta_i) {
-              // Cut i is either the same or stronger than cut r
-              if (cuts_to_remove[r] == 0) { num_cuts_to_remove++; }
-              cuts_to_remove[r] = 1;  // Remove row r
-            } else {
-              // theta_r < theta_i, so cut r is strictly stronger than cut i
-              if (cuts_to_remove[i] == 0) { num_cuts_to_remove++; }
-              cuts_to_remove[i] = 1;  // Remove row i
-            }
-          }
+    if (set_r <= 0 || set_r >= sentinel || cuts_to_remove[r] != 0) { continue; }
+    const auto& members = rows_by_set.at(set_r);
+    auto member         = std::upper_bound(members.begin(), members.end(), r);
+    for (; member != members.end(); ++member) {
+      const i_t i       = *member;
+      const f_t f_r     = divisors[r];
+      const f_t f_i     = divisors[i];
+      const f_t theta_r = rhs_storage_[r] / f_r;
+      const f_t theta_i = rhs_storage_[i] / f_i;
+      if (f_r > 0.0 && f_i > 0.0) {
+        if (theta_r <= theta_i) {
+          cuts_to_remove[r] = 1;
+        } else {
+          cuts_to_remove[i] = 1;
+        }
+      } else if (f_r < 0.0 && f_i < 0.0) {
+        if (theta_r >= theta_i) {
+          cuts_to_remove[r] = 1;
+        } else {
+          cuts_to_remove[i] = 1;
         }
       }
     }
   }
 
+  const i_t num_cuts_to_remove =
+    std::accumulate(cuts_to_remove.begin(), cuts_to_remove.end(), i_t{0});
   if (num_cuts_to_remove > 0) {
     settings_.log.debug("Removing %d duplicate cuts\n", num_cuts_to_remove);
     csr_matrix_t<i_t, f_t> new_cut_storage(0, 0, 0);