diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index 7acd7dee0a..e9f51666dc 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -14,6 +14,7 @@ #include #include +#include #include #include #include @@ -21,6 +22,7 @@ #include #include #include +#include #include #include @@ -1227,21 +1229,26 @@ f_t cut_pool_t::cut_orthogonality(i_t i, i_t j) template void cut_pool_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 divisors(cut_storage_.m, 0.0); - std::vector sets(cut_storage_.m, 0); + const i_t m = cut_storage_.m; + constexpr f_t duplicate_tolerance = 1e-10; + std::vector divisors(m, 0.0); + std::vector sets(m, 0); csc_matrix_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::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 set_groups; for (i_t j = 0; j < n; j++) { i_t r0 = -1; i_t new_rows = 0; @@ -1249,33 +1256,36 @@ void cut_pool_t::check_for_duplicate_cuts() 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::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 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. } } } diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index fcb6080178..0f6d4d634b 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -19,12 +19,14 @@ #include #include #include +#include #include #include #include #include #include +#include #include namespace cuopt::mathematical_optimization::mip { @@ -303,6 +305,105 @@ std::vector> 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 +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& set_ids) + { + build_impl(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& entry_indices, + const std::vector& set_ids) + { + build_impl(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 entries_after(i_t set_id, i_t entry) const + { + if (set_id <= 0 || set_id >= static_cast(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(end - position - 1)}; + } + + private: + template + void build_impl(i_t first, + i_t last, + i_t set_id_limit, + const std::vector& set_ids, + const std::vector* 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(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(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 bucket_by_set_; + std::vector active_set_ids_; + std::vector set_counts_; + std::vector set_starts_; + std::vector next_entry_in_set_; + std::vector entries_; + std::vector position_by_entry_; +}; + template class cut_pool_t { public: diff --git a/cpp/tests/mip/cuts_test.cu b/cpp/tests/mip/cuts_test.cu index 5af0754e3d..4932fbcf41 100644 --- a/cpp/tests/mip/cuts_test.cu +++ b/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)