diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index 7acd7dee0a..ec32afb824 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,25 @@ 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 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; @@ -1249,33 +1255,47 @@ 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]; + + std::unordered_map> 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::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> 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 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 new_cut_storage(0, 0, 0);