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,63 @@ 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 group
+  // entries by their current set in compressed storage. 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;
+  compressed_set_groups_t<i_t> set_groups;
   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];
+
+    set_groups.build_with_entry_indices(col_start, col_end, new_set_0, cut_storage_csc.i, sets);
+
     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;  // To enable us to find this new set later.
         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 old_set = sets[r];
+        bool matched      = false;
+        // Loop over all indices a_ij with i > r where i is in the same set as r.
+        for (const i_t q : set_groups.entries_after(old_set, p)) {
           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)) {
-              sets[r] = new_set;
-              sets[i] = new_set;
-            }
+          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,45 +1302,46 @@ 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
+  // 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.
+  set_groups.build(0, m, new_set, sets);
+
   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; }
+    // This cut has a duplicate. The set members are in row order, preserving the legacy
+    // strongest-cut selection order without scanning unrelated rows.
+    for (const i_t i : set_groups.entries_after(set_r, 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.0 && f_i > 0.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.
+        }
+      } else if (f_r < 0.0 && f_i < 0.0) {
+        // Dividing by a negative divisor reverses the inequality:
+        // 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.
         }
       }
     }

cpp/src/cuts/cuts.hpp

@@ -19,12 +19,14 @@
 #include <future>
 #include <memory>
 #include <numeric>
+#include <span>
 #include <string>
 #include <unordered_map>
 #include <utility>
 #include <vector>
 
 #include <cmath>
+#include <cstddef>
 #include <cstdint>
 
 namespace cuopt::mathematical_optimization::mip {
@@ -303,6 +305,105 @@ std::vector<std::vector<int>> find_mod2_row_combinations_for_test(
   double max_work_estimate,
   double* work_estimate);
 
+// Groups integer entries by sparse set id using prefix offsets and one flat member array.
+// Repeated builds reuse the allocated storage. Entry order within each set is preserved.
+template <typename i_t>
+class compressed_set_groups_t {
+ public:
+  // Groups entries in [first, last) whose set ids are in [1, set_id_limit).
+  void build(i_t first, i_t last, i_t set_id_limit, const std::vector<i_t>& set_ids)
+  {
+    build_impl<false>(first, last, set_id_limit, set_ids, nullptr);
+  }
+
+  void build_with_entry_indices(i_t first,
+                                i_t last,
+                                i_t set_id_limit,
+                                const std::vector<i_t>& entry_indices,
+                                const std::vector<i_t>& set_ids)
+  {
+    build_impl<true>(first, last, set_id_limit, set_ids, &entry_indices);
+  }
+
+  // Entry order within a set matches input order, so position + 1 is its first later entry.
+  std::span<const i_t> entries_after(i_t set_id, i_t entry) const
+  {
+    if (set_id <= 0 || set_id >= static_cast<i_t>(bucket_by_set_.size())) { return {}; }
+    const i_t bucket = bucket_by_set_[set_id];
+    if (bucket < 0) { return {}; }
+    const i_t position = position_by_entry_[entry - first_entry_];
+    const i_t end      = set_starts_[bucket + 1];
+    return {entries_.data() + position + 1, static_cast<std::size_t>(end - position - 1)};
+  }
+
+ private:
+  template <bool use_entry_indices>
+  void build_impl(i_t first,
+                  i_t last,
+                  i_t set_id_limit,
+                  const std::vector<i_t>& set_ids,
+                  const std::vector<i_t>* entry_indices)
+  {
+    first_entry_ = first;
+    for (const i_t set_id : active_set_ids_) {
+      bucket_by_set_[set_id] = -1;
+    }
+    active_set_ids_.clear();
+    set_counts_.clear();
+    if (bucket_by_set_.size() < static_cast<std::size_t>(set_id_limit)) {
+      bucket_by_set_.resize(set_id_limit, -1);
+    }
+
+    for (i_t entry = first; entry < last; entry++) {
+      i_t set_id;
+      if constexpr (use_entry_indices) {
+        set_id = set_ids[(*entry_indices)[entry]];
+      } else {
+        set_id = set_ids[entry];
+      }
+      if (set_id > 0 && set_id < set_id_limit) {
+        i_t& bucket = bucket_by_set_[set_id];
+        if (bucket < 0) {
+          bucket = static_cast<i_t>(active_set_ids_.size());
+          active_set_ids_.push_back(set_id);
+          set_counts_.push_back(0);
+        }
+        set_counts_[bucket]++;
+      }
+    }
+
+    set_starts_.resize(set_counts_.size() + 1);
+    set_starts_[0] = 0;
+    std::inclusive_scan(set_counts_.begin(), set_counts_.end(), set_starts_.begin() + 1);
+    entries_.resize(set_starts_.back());
+    position_by_entry_.resize(last - first);
+    next_entry_in_set_ = set_starts_;
+    for (i_t entry = first; entry < last; entry++) {
+      i_t set_id;
+      if constexpr (use_entry_indices) {
+        set_id = set_ids[(*entry_indices)[entry]];
+      } else {
+        set_id = set_ids[entry];
+      }
+      if (set_id > 0 && set_id < set_id_limit) {
+        const i_t bucket                  = bucket_by_set_[set_id];
+        const i_t position                = next_entry_in_set_[bucket]++;
+        entries_[position]                = entry;
+        position_by_entry_[entry - first] = position;
+      }
+    }
+  }
+
+  i_t first_entry_{0};
+  std::vector<i_t> bucket_by_set_;
+  std::vector<i_t> active_set_ids_;
+  std::vector<i_t> set_counts_;
+  std::vector<i_t> set_starts_;
+  std::vector<i_t> next_entry_in_set_;
+  std::vector<i_t> entries_;
+  std::vector<i_t> position_by_entry_;
+};
+
 template <typename i_t, typename f_t>
 class cut_pool_t {
  public:

cpp/tests/mip/cuts_test.cu

@@ -981,6 +981,7 @@ TEST(cuts, test_duplicate_cuts_detection)
   cut_pool.add_cut(mip::cut_type_t::MIXED_INTEGER_GOMORY, cut8);
 
   cut_pool.check_for_duplicate_cuts();
+  EXPECT_EQ(cut_pool.pool_size(), 5);
 }
 
 TEST(cuts, clique_phase1_smoke_conflict_graph_edges)