diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index 4ce11b830b..9ee51e92d9 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -277,7 +277,7 @@ FetchContent_Declare( # This is the reason we are using the development branch # from Oct 12, 2025. Once these changes are merged into the main branch, #we can switch to the main branch. - GIT_TAG "32b3a87dbf4955d5a2803be74145c389ea31434d" + GIT_TAG "55d5edece584885061639ecf3a6eb8a4629be9f2" GIT_PROGRESS TRUE EXCLUDE_FROM_ALL SYSTEM diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 4dc6bc67a8..6fdf374c2e 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -1565,8 +1565,9 @@ dual_status_t branch_and_bound_t::solve_node_lp( } else { lp_settings.cut_off = cutoff + settings_.dual_tol; } - lp_settings.inside_mip = 2; - lp_settings.time_limit = settings_.time_limit - toc(exploration_stats_.start_time); + lp_settings.inside_mip = 2; + lp_settings.time_limit = settings_.time_limit - toc(exploration_stats_.start_time); + if (lp_settings.time_limit <= 0.0) { return dual_status_t::TIME_LIMIT; } lp_settings.scale_columns = false; lp_settings.iteration_limit = iter_limit; @@ -1727,7 +1728,8 @@ void branch_and_bound_t::plunge_with(bfs_worker_t* worker, } if (now > settings_.time_limit) { - solver_status_ = mip_status_t::TIME_LIMIT; + node_concurrent_halt_ = 1; + solver_status_ = mip_status_t::TIME_LIMIT; stack.push_front(node_ptr); --exploration_stats_.nodes_being_solved; break; @@ -1760,7 +1762,8 @@ void branch_and_bound_t::plunge_with(bfs_worker_t* worker, --exploration_stats_.nodes_being_solved; if (lp_status == dual_status_t::TIME_LIMIT) { - solver_status_ = mip_status_t::TIME_LIMIT; + node_concurrent_halt_ = 1; + solver_status_ = mip_status_t::TIME_LIMIT; stack.push_front(node_ptr); break; } @@ -1948,7 +1951,8 @@ void branch_and_bound_t::best_first_search_with(bfs_worker_t } if (toc(exploration_stats_.start_time) > settings_.time_limit) { - solver_status_ = mip_status_t::TIME_LIMIT; + node_concurrent_halt_ = 1; + solver_status_ = mip_status_t::TIME_LIMIT; break; } @@ -2044,7 +2048,8 @@ void branch_and_bound_t::dive_with(diving_worker_t* worker, } if (toc(exploration_stats_.start_time) > settings_.time_limit) { - solver_status_ = mip_status_t::TIME_LIMIT; + node_concurrent_halt_ = 1; + solver_status_ = mip_status_t::TIME_LIMIT; break; } if (dive_stats.nodes_explored >= diving_node_limit) { break; } @@ -2063,7 +2068,8 @@ void branch_and_bound_t::dive_with(diving_worker_t* worker, ++dive_stats.nodes_explored; if (lp_status == dual_status_t::TIME_LIMIT) { - solver_status_ = mip_status_t::TIME_LIMIT; + node_concurrent_halt_ = 1; + solver_status_ = mip_status_t::TIME_LIMIT; break; } if (lp_status == dual_status_t::CONCURRENT_LIMIT) { break; } diff --git a/cpp/src/cuts/CMakeLists.txt b/cpp/src/cuts/CMakeLists.txt index 813ac88a59..2d4412d00a 100644 --- a/cpp/src/cuts/CMakeLists.txt +++ b/cpp/src/cuts/CMakeLists.txt @@ -6,6 +6,7 @@ set(CUTS_SRC_FILES ${CMAKE_CURRENT_SOURCE_DIR}/cuts.cpp ${CMAKE_CURRENT_SOURCE_DIR}/objective_step.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/zero_half_mod2.cpp ) set(CUOPT_SRC_FILES ${CUOPT_SRC_FILES} diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index be45ffeecd..3fa79f88c7 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -2324,7 +2324,8 @@ i_t knapsack_generation_t::generate_knapsack_cut( const std::vector& var_types, const std::vector& xstar, i_t knapsack_row, - inequality_t& cut) + inequality_t& cut, + f_t start_time) { const bool verbose = false; // Get the row associated with the knapsack constraint @@ -2499,7 +2500,8 @@ i_t knapsack_generation_t::generate_knapsack_cut( // Lift the cut inequality_t lifted_cut(lp.num_cols); - lift_knapsack_cut(knapsack_inequality, minimal_cover_cut, c1_partition, c2_partition, lifted_cut); + lift_knapsack_cut( + knapsack_inequality, minimal_cover_cut, c1_partition, c2_partition, lifted_cut, start_time); lifted_cut.negate(); // The cut is now in the form: @@ -2679,7 +2681,8 @@ void knapsack_generation_t::lift_knapsack_cut( const inequality_t& base_cut, const std::vector& c1_partition, const std::vector& c2_partition, - inequality_t& lifted_cut) + inequality_t& lifted_cut, + f_t start_time) { // The base cut is in the form: sum_{j in cover} x_j <= |cover| - 1 @@ -2795,14 +2798,15 @@ void knapsack_generation_t::lift_knapsack_cut( best_score_last_permutation(remaining_coefficients, permutation); while (permutation.size() > 0) { + if (toc(start_time) >= settings_.time_limit) { break; } const i_t h = permutation.back(); const i_t k = remaining_variables[h]; const f_t a_k = remaining_coefficients[h]; f_t capacity = knapsack_inequality.rhs - a_k; - f_t objective = - exact_knapsack_problem_integer_values_fraction_values(values, weights, capacity, solution); + f_t objective = exact_knapsack_problem_integer_values_fraction_values( + values, weights, capacity, solution, start_time); if (std::isnan(objective)) { settings_.log.debug("lifting knapsack problem failed\n"); break; @@ -3021,8 +3025,10 @@ f_t knapsack_generation_t::exact_knapsack_problem_integer_values_fract const std::vector& values, const std::vector& weights, f_t rhs, - std::vector& solution) + std::vector& solution, + f_t start_time) { + if (toc(start_time) >= settings_.time_limit) { return std::numeric_limits::quiet_NaN(); } // Solve the knapsack problem // maximize sum_{j=0}^n values[j] * solution[j] // subject to sum_{j=0}^n weights[j] * solution[j] <= rhs @@ -3050,6 +3056,7 @@ f_t knapsack_generation_t::exact_knapsack_problem_integer_values_fract // 4. Dynamic programming for (i_t j = 1; j <= n; ++j) { + if (toc(start_time) >= settings_.time_limit) { return std::numeric_limits::quiet_NaN(); } for (i_t v = 0; v <= sum_value; ++v) { // Do not take item i-1 dp(j, v) = dp(j - 1, v); @@ -3095,19 +3102,27 @@ void cut_generation_t::generate_implied_bound_cuts( { if (probing_implied_bound_.zero_offsets.empty()) { return; } - const f_t tol = 1e-4; - i_t num_cuts = 0; - const i_t pib_cols = static_cast(probing_implied_bound_.zero_offsets.size()) - 1; - const i_t n_cols = std::min(lp.num_cols, pib_cols); + const f_t tol = 1e-4; + i_t num_cuts = 0; + const i_t pib_cols = static_cast(probing_implied_bound_.zero_offsets.size()) - 1; + const i_t n_cols = std::min(lp.num_cols, pib_cols); + f_t work_estimate = 0.0; + const f_t max_work_estimate = 1e8; + constexpr f_t implication_work = 16.0; + constexpr f_t generated_cut_work = 16.0; for (i_t j = 0; j < n_cols; j++) { + if (work_estimate > max_work_estimate || toc(start_time) >= settings.time_limit) { return; } if (var_types[j] == variable_type_t::CONTINUOUS) { continue; } const f_t xstar_j = xstar[j]; // x_j = 0 implications const i_t zero_begin = probing_implied_bound_.zero_offsets[j]; const i_t zero_end = probing_implied_bound_.zero_offsets[j + 1]; + const i_t one_begin = probing_implied_bound_.one_offsets[j]; + const i_t one_end = probing_implied_bound_.one_offsets[j + 1]; for (i_t p = zero_begin; p < zero_end; p++) { + work_estimate += implication_work; const i_t i = probing_implied_bound_.zero_variables[p]; if (i == j) { continue; } const f_t l_i = lp.lower[i]; @@ -3127,6 +3142,7 @@ void cut_generation_t::generate_implied_bound_cuts( cut.push_back(j, coeff_j); cut.rhs = -b_ub; cut_pool_.add_cut(cut_type_t::IMPLIED_BOUND, cut); + work_estimate += generated_cut_work; num_cuts++; } } @@ -3145,15 +3161,16 @@ void cut_generation_t::generate_implied_bound_cuts( cut.push_back(j, coeff_j); cut.rhs = b_lb; cut_pool_.add_cut(cut_type_t::IMPLIED_BOUND, cut); + work_estimate += generated_cut_work; num_cuts++; } } } + if (work_estimate > max_work_estimate || toc(start_time) >= settings.time_limit) { return; } // x_j = 1 implications - const i_t one_begin = probing_implied_bound_.one_offsets[j]; - const i_t one_end = probing_implied_bound_.one_offsets[j + 1]; for (i_t p = one_begin; p < one_end; p++) { + work_estimate += implication_work; const i_t i = probing_implied_bound_.one_variables[p]; if (i == j) { continue; } const f_t l_i = lp.lower[i]; @@ -3173,6 +3190,7 @@ void cut_generation_t::generate_implied_bound_cuts( cut.push_back(j, coeff_j); cut.rhs = -u_i; cut_pool_.add_cut(cut_type_t::IMPLIED_BOUND, cut); + work_estimate += generated_cut_work; num_cuts++; } } @@ -3190,10 +3208,12 @@ void cut_generation_t::generate_implied_bound_cuts( cut.push_back(j, coeff_j); cut.rhs = rhs_val; cut_pool_.add_cut(cut_type_t::IMPLIED_BOUND, cut); + work_estimate += generated_cut_work; num_cuts++; } } } + if (work_estimate > max_work_estimate || toc(start_time) >= settings.time_limit) { return; } } if (num_cuts > 0) { @@ -3604,7 +3624,8 @@ bool cut_generation_t::generate_cuts(const lp_problem_t& lp, if (toc(start_time) >= settings.time_limit) { return true; } ZERO_HALF_DEBUG("generate_cuts: about to call generate_zero_half_cuts"); f_t cut_start_time = tic(); - bool feasible = generate_zero_half_cuts(lp, settings, var_types, xstar, zstar, start_time); + bool feasible = generate_zero_half_cuts( + lp, settings, Arow, new_slacks, var_types, xstar, zstar, variable_bounds, start_time); ZERO_HALF_DEBUG("generate_cuts: returned from generate_zero_half_cuts feasible=%d", static_cast(feasible)); if (!feasible) { @@ -3637,7 +3658,7 @@ void cut_generation_t::generate_knapsack_cuts( if (toc(start_time) >= settings.time_limit) { return; } inequality_t cut(lp.num_cols); i_t knapsack_status = knapsack_generation_.generate_knapsack_cut( - lp, settings, Arow, new_slacks, var_types, xstar, knapsack_row, cut); + lp, settings, Arow, new_slacks, var_types, xstar, knapsack_row, cut, start_time); if (knapsack_status == 0) { cut_pool_.add_cut(cut_type_t::KNAPSACK, cut); } } } @@ -3855,9 +3876,12 @@ template bool cut_generation_t::generate_zero_half_cuts( const lp_problem_t& lp, const simplex_solver_settings_t& settings, + csr_matrix_t& Arow, + const std::vector& new_slacks, const std::vector& var_types, const std::vector& xstar, const std::vector& reduced_costs, + variable_bounds_t& variable_bounds, f_t start_time) { if (settings.zero_half_cuts == 0) { return true; } @@ -3882,12 +3906,25 @@ bool cut_generation_t::generate_zero_half_cuts( static_cast(sub_cg_.ready), sub_cg_.vertices.size()); + f_t mod2_work_estimate = 0.0; + const bool mod2_completed = generate_mod2_zero_half_cuts(cut_pool_, + lp, + settings, + Arow, + new_slacks, + var_types, + xstar, + variable_bounds, + start_time, + mod2_work_estimate); + if (!mod2_completed) { return true; } + // The fractional conflict-graph subgraph is built once per cut pass in - // prepare_fractional_sub_conflict_graph() (called from generate_cuts) and shared with - // the clique-cut separator. Skip if the build was unable to produce a - // useable sub-CG (clique table missing/empty, work/time budget hit, etc.). + // prepare_fractional_sub_conflict_graph() and remains a complementary + // odd-cycle / odd-wheel separator. If no conflict graph is available, the + // general row-parity cuts above are still retained. if (!sub_cg_.ready) { - ZERO_HALF_DEBUG("sub_cg_ not ready, skipping"); + ZERO_HALF_DEBUG("sub_cg_ not ready, skipping odd-cycle path"); return true; } if (sub_cg_.empty_subgraph()) { @@ -3905,8 +3942,8 @@ bool cut_generation_t::generate_zero_half_cuts( cuopt_assert(user_problem_.var_types.size() == static_cast(num_vars), "Zero-half user problem var_types size mismatch"); - const f_t min_violation = std::max(settings.primal_tol, static_cast(1e-6)); - const f_t bound_tol = settings.primal_tol; + constexpr f_t min_violation = (f_t)1e-6; + const f_t bound_tol = settings.primal_tol; // shortest path of length >= 0.5 - min_violation cannot yield a violated cut const f_t cutoff = static_cast(0.5) - min_violation; f_t work_estimate = 0.0; @@ -4070,21 +4107,22 @@ void cut_generation_t::generate_mir_cuts( // at the beginning of each iteration of the for loop below std::vector aggregated_rows; std::vector aggregated_mark(lp.num_rows, 0); + const i_t max_cuts = std::min(lp.num_rows, 100000); // Transform the relaxation solution std::vector transformed_xstar; complemented_mir.bound_substitution(lp, variable_bounds, var_types, xstar, transformed_xstar); - const i_t max_cuts = std::min(lp.num_rows, 100000); f_t work_estimate = 0.0; - i_t num_cuts = 0; - while (num_cuts < max_cuts && !score_queue.empty()) { + i_t cuts_processed = 0; + while (cuts_processed < max_cuts && !score_queue.empty()) { if (toc(start_time) >= settings.time_limit) { break; } // Get the row with the highest score from the queue auto [max_score, i] = score_queue.top(); score_queue.pop(); // skip stale score entries if (max_score != scores[i]) { continue; } + ++cuts_processed; // Add the current row to the aggregated set aggregated_mark[i] = 1; @@ -5304,7 +5342,8 @@ void complemented_mixed_integer_rounding_cut_t::bound_substitution( const variable_bounds_t& variable_bounds, const std::vector& var_types, const std::vector& xstar, - std::vector& transformed_xstar) + std::vector& transformed_xstar, + bool prefer_variable_bound_on_tie) { transformed_xstar.resize(lp.num_cols); // Perform bound substitution for continuous variables @@ -5368,8 +5407,12 @@ void complemented_mixed_integer_rounding_cut_t::bound_substitution( bound_changed_[j] = 0; continue; } - if (has_finite_lower_bound && - (!has_finite_upper_bound || (xstar_j - lb_star_[j] <= ub_star_[j] - xstar_j))) { + const f_t lower_distance = xstar_j - lb_star_[j]; + const f_t upper_distance = ub_star_[j] - xstar_j; + const bool prefer_upper_variable_bound = + prefer_variable_bound_on_tie && ub_variable_[j] >= 0 && lower_distance == upper_distance; + if (has_finite_lower_bound && (!has_finite_upper_bound || (lower_distance <= upper_distance && + !prefer_upper_variable_bound))) { // Use the lower bound // lb_star_j <= x_j <= ub_star_j // v_j = x_j - lb_star_j, @@ -5693,7 +5736,10 @@ f_t complemented_mixed_integer_rounding_cut_t::compute_violation( template void complemented_mixed_integer_rounding_cut_t::substitute_slacks( - const lp_problem_t& lp, csr_matrix_t& Arow, inequality_t& cut) + const lp_problem_t& lp, + csr_matrix_t& Arow, + inequality_t& cut, + f_t* work_estimate) { // Remove slacks from the cut // So that the cut is only over the original variables @@ -5701,6 +5747,7 @@ void complemented_mixed_integer_rounding_cut_t::substitute_slacks( i_t cut_nz = 0; std::vector cut_indices; cut_indices.reserve(cut.size()); + if (work_estimate != nullptr) { *work_estimate += cut.size(); } for (i_t k = 0; k < cut.size(); k++) { const i_t j = cut.index(k); @@ -5743,6 +5790,7 @@ void complemented_mixed_integer_rounding_cut_t::substitute_slacks( cut.rhs -= cj * lp.rhs[i] / alpha; const i_t row_start = Arow.row_start[i]; const i_t row_end = Arow.row_start[i + 1]; + if (work_estimate != nullptr) { *work_estimate += row_end - row_start; } for (i_t q = row_start; q < row_end; q++) { const i_t h = Arow.j[q]; if (h != j) { @@ -5764,6 +5812,9 @@ void complemented_mixed_integer_rounding_cut_t::substitute_slacks( if (found_slack) { scratch_pad_.get_pad(cut.vector.i, cut.vector.x); + if (work_estimate != nullptr) { + *work_estimate += 2 * cut.size() + cut.size() * std::log2((f_t)cut.size() + (f_t)1.0); + } // Sort the cut cut.sort(); } diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index 78091c85f6..fcb6080178 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -286,6 +286,23 @@ std::vector> find_violated_odd_cycles_for_test( double min_violation, double time_limit); +// Test-only helper to run the production sparse GF(2) row-dependency finder used +// by general zero-half cuts. Each parity row contains the integer-variable +// indices with an odd coefficient. A returned combination has even aggregate +// parity and an odd aggregate rhs. +std::vector> find_mod2_row_combinations_for_test( + const std::vector>& parity_rows, + const std::vector& rhs_parity, + int max_combination_size, + int max_combinations); +std::vector> find_mod2_row_combinations_for_test( + const std::vector>& parity_rows, + const std::vector& rhs_parity, + int max_combination_size, + int max_combinations, + double max_work_estimate, + double* work_estimate); + template class cut_pool_t { public: @@ -346,6 +363,18 @@ class cut_pool_t { template class variable_bounds_t; +template +bool generate_mod2_zero_half_cuts(cut_pool_t& cut_pool, + const simplex::lp_problem_t& lp, + const simplex::simplex_solver_settings_t& settings, + csr_matrix_t& Arow, + const std::vector& new_slacks, + const std::vector& var_types, + const std::vector& xstar, + variable_bounds_t& variable_bounds, + f_t start_time, + f_t& work_estimate); + template struct flow_cover_row_t { i_t row; @@ -543,7 +572,8 @@ class knapsack_generation_t { const std::vector& var_types, const std::vector& xstar, i_t knapsack_row, - inequality_t& cut); + inequality_t& cut, + f_t start_time); i_t num_knapsack_constraints() const { return knapsack_constraints_.size(); } const std::vector& get_knapsack_constraints() const { return knapsack_constraints_; } @@ -568,7 +598,8 @@ class knapsack_generation_t { const inequality_t& base_cut, const std::vector& c1_partition, const std::vector& c2_partition, - inequality_t& lifted_cut); + inequality_t& lifted_cut, + f_t start_time); // Solve a 0-1 knapsack problem using dynamic programming f_t solve_knapsack_problem(const std::vector& values, @@ -579,7 +610,8 @@ class knapsack_generation_t { f_t exact_knapsack_problem_integer_values_fraction_values(const std::vector& values, const std::vector& weights, f_t rhs, - std::vector& solution); + std::vector& solution, + f_t start_time); std::vector is_slack_; std::vector knapsack_constraints_; @@ -709,12 +741,16 @@ class cut_generation_t { const std::vector& reduced_costs, f_t start_time); - // Generate zero-half (odd-cycle / odd-wheel) cuts from the conflict graph + // Generate general row-parity zero-half cuts and conflict-graph + // odd-cycle / odd-wheel cuts. bool generate_zero_half_cuts(const simplex::lp_problem_t& lp, const simplex::simplex_solver_settings_t& settings, + csr_matrix_t& Arow, + const std::vector& new_slacks, const std::vector& var_types, const std::vector& xstar, const std::vector& reduced_costs, + variable_bounds_t& variable_bounds, f_t start_time); // Generate implied bounds cuts from probing implications @@ -955,7 +991,8 @@ class complemented_mixed_integer_rounding_cut_t { const variable_bounds_t& variable_bounds, const std::vector& var_types, const std::vector& xstar, - std::vector& transformed_xstar); + std::vector& transformed_xstar, + bool prefer_variable_bound_on_tie = false); // Converts an inequality of the form: sum_j a_j x_j >= beta // with l_j <= x_j <= u_j into the form: @@ -998,6 +1035,15 @@ class complemented_mixed_integer_rounding_cut_t { const std::vector& var_types, inequality_t& cut); + // Generate a lifted mixed-binary cover inequality from a transformed + // nonnegative >= row. Returns the cut in the same >= convention. + bool generate_lifted_mixed_binary_cover(const inequality_t& transformed_inequality, + const std::vector& var_types, + const std::vector& transformed_xstar, + inequality_t& transformed_cut, + f_t& work_estimate, + f_t max_work_estimate); + f_t compute_violation(const inequality_t& cut, const std::vector& xstar); f_t new_upper(i_t j) const { return transformed_upper_[j]; } @@ -1011,7 +1057,8 @@ class complemented_mixed_integer_rounding_cut_t { void substitute_slacks(const simplex::lp_problem_t& lp, csr_matrix_t& Arow, - inequality_t& cut); + inequality_t& cut, + f_t* work_estimate = nullptr); // Combine the pivot row with the inequality to eliminate the variable j // The new inequality is returned in inequality and inequality_rhs diff --git a/cpp/src/cuts/zero_half_mod2.cpp b/cpp/src/cuts/zero_half_mod2.cpp new file mode 100644 index 0000000000..13e2413fa1 --- /dev/null +++ b/cpp/src/cuts/zero_half_mod2.cpp @@ -0,0 +1,774 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace cuopt::mathematical_optimization::mip { + +using simplex::lp_problem_t; +using simplex::simplex_solver_settings_t; +using simplex::variable_type_t; + +namespace { + +template +void symmetric_difference_sorted(const std::vector& a, + const std::vector& b, + std::vector& result) +{ + result.clear(); + result.reserve(a.size() + b.size()); + std::set_symmetric_difference(a.begin(), a.end(), b.begin(), b.end(), std::back_inserter(result)); +} + +template +struct mod2_parity_row_t { + std::vector parity; + bool rhs_parity{false}; +}; + +template +struct mod2_candidate_t : mod2_parity_row_t { + inequality_t transformed_inequality; + bool reversible{false}; +}; + +template +struct mod2_basis_row_t { + std::vector parity; + std::vector combination; + bool rhs{false}; +}; + +template +struct mod2_row_order_t { + const std::vector& rows; + + bool operator()(i_t a, i_t b) const + { + if (rows[a].parity.size() != rows[b].parity.size()) { + return rows[a].parity.size() < rows[b].parity.size(); + } + return rows[a].rhs_parity < rows[b].rhs_parity; + } +}; + +template +std::vector> find_mod2_row_combinations(const std::vector& rows, + i_t max_combination_size, + i_t max_combinations, + f_t* work_estimate, + f_t max_work_estimate, + f_t start_time = 0.0, + f_t time_limit = inf) +{ + cuopt_assert(max_combination_size > 0, "Maximum GF(2) combination size must be positive"); + cuopt_assert(max_combinations > 0, "Maximum number of GF(2) combinations must be positive"); + + i_t max_index = -1; + f_t input_scan_work = 0.0; + for (const auto& row : rows) { + input_scan_work += (f_t)(row.parity.size() + 1); + cuopt_assert(std::is_sorted(row.parity.begin(), row.parity.end()), + "GF(2) parity rows must be sorted"); + cuopt_assert(std::adjacent_find(row.parity.begin(), row.parity.end()) == row.parity.end(), + "GF(2) parity rows must not contain duplicates"); + if (!row.parity.empty()) { + cuopt_assert(row.parity.front() >= 0, "GF(2) parity index must be nonnegative"); + max_index = std::max(max_index, row.parity.back()); + } + } + if (add_work_estimate(input_scan_work, work_estimate, max_work_estimate)) { return {}; } + + std::vector permutation(rows.size()); + std::iota(permutation.begin(), permutation.end(), 0); + const f_t sort_work = (f_t)permutation.size() * std::log2((f_t)permutation.size() + (f_t)1.0); + if (add_work_estimate(sort_work, work_estimate, max_work_estimate)) { return {}; } + // this is to process small/sparse rows first, for faster perf and smaller combinations + std::stable_sort(permutation.begin(), permutation.end(), mod2_row_order_t{rows}); + + if (add_work_estimate((f_t)(max_index + 1), work_estimate, max_work_estimate)) { return {}; } + std::vector pivot_to_basis((size_t)(max_index + 1), -1); + std::vector> basis; + basis.reserve(std::min(rows.size(), (size_t)(max_index + 1))); + std::vector> combinations; + combinations.reserve(std::min((size_t)max_combinations, rows.size())); + + std::vector parity_tmp; + std::vector combination_tmp; + for (const i_t candidate : permutation) { + if (toc(start_time) >= time_limit) { break; } + f_t candidate_work = (f_t)(rows[candidate].parity.size() + 2); + mod2_basis_row_t current; + current.parity = rows[candidate].parity; + current.combination = {candidate}; + current.rhs = rows[candidate].rhs_parity; + + bool abandoned = false; + while (!current.parity.empty()) { + const i_t pivot = current.parity.front(); + const i_t basis_index = pivot_to_basis[pivot]; + // pivot has not been seen before + if (basis_index < 0) { break; } + + const auto& pivot_row = basis[basis_index]; + candidate_work += (f_t)(current.parity.size() + pivot_row.parity.size() + + current.combination.size() + pivot_row.combination.size()); + symmetric_difference_sorted(current.parity, pivot_row.parity, parity_tmp); + symmetric_difference_sorted(current.combination, pivot_row.combination, combination_tmp); + if (combination_tmp.size() > (size_t)max_combination_size) { + abandoned = true; + break; + } + current.parity.swap(parity_tmp); + current.combination.swap(combination_tmp); + current.rhs = current.rhs != pivot_row.rhs; + } + if (add_work_estimate(candidate_work, work_estimate, max_work_estimate)) { break; } + if (abandoned) { continue; } + + // when reduced, add to combinations and continue, don't add to basis + if (current.parity.empty()) { + if (current.rhs && !current.combination.empty()) { + combinations.push_back(std::move(current.combination)); + if (combinations.size() >= (size_t)max_combinations) { break; } + } + continue; + } + + const i_t pivot = current.parity.front(); + pivot_to_basis[pivot] = (i_t)basis.size(); + basis.push_back(std::move(current)); + } + return combinations; +} + +template +i_t mod2_integral_scale(const inequality_t& inequality, + const std::vector& var_types, + const std::vector& transformed_xstar, + i_t max_integral_scale, + f_t row_tight_tol, + f_t coefficient_integral_tol, + f_t start_time, + f_t time_limit, + f_t& work_estimate, + f_t max_work_estimate, + bool& work_limit_reached) +{ + for (i_t scale = 1; scale <= max_integral_scale; ++scale) { + if (toc(start_time) >= time_limit || work_limit_reached) { return i_t{0}; } + f_t scale_work = 1.0; + bool integral = true; + const f_t scaled_rhs = (f_t)scale * inequality.rhs; + if (std::abs(scaled_rhs - std::round(scaled_rhs)) > + coefficient_integral_tol * std::max((f_t)1.0, std::abs(scaled_rhs))) { + integral = false; + } + for (i_t k = 0; integral && k < (i_t)inequality.size(); ++k) { + scale_work += 1.0; + const i_t j = inequality.index(k); + if (var_types[j] == variable_type_t::CONTINUOUS || transformed_xstar[j] <= row_tight_tol) { + continue; + } + const f_t scaled_coefficient = (f_t)scale * inequality.coeff(k); + if (std::abs(scaled_coefficient - std::round(scaled_coefficient)) > + coefficient_integral_tol * std::max((f_t)1.0, std::abs(scaled_coefficient))) { + integral = false; + } + } + if (add_work_estimate(scale_work, &work_estimate, max_work_estimate, &work_limit_reached)) { + return i_t{0}; + } + if (integral) { return scale; } + } + return i_t{0}; +} + +template +std::vector> mod2_collect_candidates( + complemented_mixed_integer_rounding_cut_t& complemented_mir, + const lp_problem_t& lp, + csr_matrix_t& Arow, + const variable_bounds_t& variable_bounds, + const std::vector& var_types, + const std::vector& transformed_xstar, + f_t start_time, + f_t time_limit, + f_t& work_estimate, + f_t max_work_estimate, + bool& work_limit_reached) +{ + constexpr i_t max_integral_scale = 1000; + const i_t max_integer_row_length = 1000 + lp.num_cols / 10; + constexpr f_t row_tight_tol = (f_t)1e-6; + constexpr f_t coefficient_integral_tol = (f_t)1e-6; + + std::vector> candidates; + candidates.reserve(lp.num_rows); + for (i_t row = 0; row < lp.num_rows; ++row) { + if (toc(start_time) >= time_limit || work_limit_reached) { break; } + const i_t slack = complemented_mir.slack_cols(row); + if (slack < 0 || transformed_xstar[slack] > row_tight_tol) { continue; } + + const i_t row_length = Arow.row_start[row + 1] - Arow.row_start[row]; + if (row_length > max_integer_row_length) { continue; } + const f_t row_work = + (f_t)(8 * row_length + 5) + (f_t)row_length * std::log2((f_t)row_length + (f_t)1.0); + if (add_work_estimate(row_work, &work_estimate, max_work_estimate, &work_limit_reached)) { + break; + } + inequality_t inequality(Arow, row, lp.rhs[row]); + complemented_mir.transform_inequality(variable_bounds, var_types, inequality); + inequality.sort(); + + // Every LP row is an equality after slack insertion. Remove a zero-valued transformed slack + // in the direction that preserves a valid >= inequality. + i_t slack_position = -1; + for (i_t k = 0; k < (i_t)inequality.size(); ++k) { + if (inequality.index(k) == slack) { + slack_position = k; + break; + } + } + if (slack_position < 0 || inequality.coeff(slack_position) == 0.0) { continue; } + // we want a row that is a.x >= b + if (inequality.coeff(slack_position) > 0.0) { inequality.negate(); } + inequality.vector.x[slack_position] = 0.0; + inequality_t squeezed_inequality(lp.num_cols); + inequality.squeeze(squeezed_inequality); + inequality = std::move(squeezed_inequality); + + // Continuous variables must be at their selected bounds to participate in the parity system. + bool continuous_at_bounds = true; + for (i_t k = 0; k < (i_t)inequality.size(); ++k) { + const i_t j = inequality.index(k); + if (var_types[j] == variable_type_t::CONTINUOUS && + std::abs(inequality.coeff(k)) > coefficient_integral_tol && + transformed_xstar[j] > row_tight_tol) { + continuous_at_bounds = false; + break; + } + } + if (!continuous_at_bounds) { continue; } + + const i_t scale = mod2_integral_scale(inequality, + var_types, + transformed_xstar, + max_integral_scale, + row_tight_tol, + coefficient_integral_tol, + start_time, + time_limit, + work_estimate, + max_work_estimate, + work_limit_reached); + // no integral scale found or time limit reached + if (scale == 0) { continue; } + if (scale != 1) { inequality.scale((f_t)scale); } + + mod2_candidate_t candidate; + candidate.transformed_inequality = std::move(inequality); + candidate.rhs_parity = (std::abs(std::llround(candidate.transformed_inequality.rhs)) % 2) != 0; + // checks if this could be safely reversed + candidate.reversible = std::abs(lp.upper[slack] - lp.lower[slack]) <= row_tight_tol; + for (i_t k = 0; k < (i_t)candidate.transformed_inequality.size(); ++k) { + const i_t j = candidate.transformed_inequality.index(k); + if (var_types[j] == variable_type_t::CONTINUOUS || transformed_xstar[j] <= row_tight_tol) { + continue; + } + const auto coefficient = std::llround(candidate.transformed_inequality.coeff(k)); + if ((std::abs(coefficient) % 2) != 0) { candidate.parity.push_back(j); } + } + if (candidate.parity.size() > (size_t)max_integer_row_length) { continue; } + candidates.push_back(std::move(candidate)); + } + return candidates; +} + +template +void mod2_add_transformed_zero_half_cut( + complemented_mixed_integer_rounding_cut_t& complemented_mir, + cut_pool_t& cut_pool, + const lp_problem_t& lp, + csr_matrix_t& Arow, + const variable_bounds_t& variable_bounds, + const std::vector& var_types, + const std::vector& xstar, + inequality_t transformed_cut, + f_t min_violation, + f_t& work_estimate, + i_t& cuts_added) +{ + work_estimate += (f_t)(4 * transformed_cut.size() + 1); + complemented_mir.untransform_inequality(variable_bounds, var_types, transformed_cut); + complemented_mir.remove_small_coefficients(lp.lower, lp.upper, transformed_cut); + complemented_mir.substitute_slacks(lp, Arow, transformed_cut, &work_estimate); + complemented_mir.remove_small_coefficients(lp.lower, lp.upper, transformed_cut); + const f_t violation = complemented_mir.compute_violation(transformed_cut, xstar); + if (violation > min_violation) { + const i_t pool_size = cut_pool.pool_size(); + cut_pool.add_cut(cut_type_t::ZERO_HALF, transformed_cut); + if (cut_pool.pool_size() > pool_size) { ++cuts_added; } + } +} + +template +void mod2_generate_cuts_from_aggregate( + complemented_mixed_integer_rounding_cut_t& complemented_mir, + cut_pool_t& cut_pool, + const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + csr_matrix_t& Arow, + const variable_bounds_t& variable_bounds, + const std::vector& var_types, + const std::vector& xstar, + const std::vector& transformed_xstar, + const inequality_t& oriented_aggregate, + f_t min_violation, + f_t start_time, + f_t& work_estimate, + f_t max_work_estimate, + bool& work_limit_reached, + i_t& cuts_added) +{ + if (add_work_estimate((f_t)(3 * oriented_aggregate.size() + 1), + &work_estimate, + max_work_estimate, + &work_limit_reached)) { + return; + } + inequality_t mir_cut(lp.num_cols); + const bool mir_cut_generated = complemented_mir.generate_cut_nonnegative_maintain_indicies( + oriented_aggregate, var_types, mir_cut); + if (mir_cut_generated) { + mod2_add_transformed_zero_half_cut(complemented_mir, + cut_pool, + lp, + Arow, + variable_bounds, + var_types, + xstar, + std::move(mir_cut), + min_violation, + work_estimate, + cuts_added); + } + + if (work_estimate > max_work_estimate) { + work_limit_reached = true; + return; + } + inequality_t lifted_cover_cut(lp.num_cols); + bool lifted_cover_cut_generated = false; + if (toc(start_time) < settings.time_limit) { + lifted_cover_cut_generated = + complemented_mir.generate_lifted_mixed_binary_cover(oriented_aggregate, + var_types, + transformed_xstar, + lifted_cover_cut, + work_estimate, + max_work_estimate); + } + if (lifted_cover_cut_generated) { + mod2_add_transformed_zero_half_cut(complemented_mir, + cut_pool, + lp, + Arow, + variable_bounds, + var_types, + xstar, + std::move(lifted_cover_cut), + min_violation, + work_estimate, + cuts_added); + } + if (work_estimate > max_work_estimate) { work_limit_reached = true; } +} + +template +struct lifted_cover_order_t { + const std::vector& solution_value; + const inequality_t& base; + f_t tolerance; + + bool operator()(int a, int b) const + { + const bool a_at_upper = solution_value[a] >= 1.0 - tolerance; + const bool b_at_upper = solution_value[b] >= 1.0 - tolerance; + if (a_at_upper != b_at_upper) { return a_at_upper; } + const f_t contribution_a = solution_value[a] * base.coeff(a); + const f_t contribution_b = solution_value[b] * base.coeff(b); + if (contribution_a != contribution_b) { return contribution_a > contribution_b; } + return base.coeff(a) > base.coeff(b); + } +}; + +template +f_t lifted_cover_coefficient( + f_t coefficient, const std::vector& prefix, size_t p, f_t lambda, f_t tolerance) +{ + for (size_t h = 0; h < p; ++h) { + if (coefficient <= prefix[h] - lambda + tolerance) { return (f_t)h * lambda; } + if (coefficient <= prefix[h] + tolerance) { + return (f_t)(h + 1) * lambda + coefficient - prefix[h]; + } + } + return (f_t)p * lambda + coefficient - prefix[p - 1]; +} + +} // namespace + +std::vector> find_mod2_row_combinations_for_test( + const std::vector>& parity_rows, + const std::vector& rhs_parity, + int max_combination_size, + int max_combinations) +{ + return find_mod2_row_combinations_for_test(parity_rows, + rhs_parity, + max_combination_size, + max_combinations, + std::numeric_limits::infinity(), + nullptr); +} + +std::vector> find_mod2_row_combinations_for_test( + const std::vector>& parity_rows, + const std::vector& rhs_parity, + int max_combination_size, + int max_combinations, + double max_work_estimate, + double* work_estimate_out) +{ + cuopt_assert(parity_rows.size() == rhs_parity.size(), + "GF(2) parity row and rhs sizes must match"); + std::vector> rows; + rows.reserve(parity_rows.size()); + for (size_t i = 0; i < parity_rows.size(); ++i) { + rows.push_back({parity_rows[i], rhs_parity[i] != 0}); + } + + double work_estimate = 0.0; + auto combinations = find_mod2_row_combinations( + rows, max_combination_size, max_combinations, &work_estimate, max_work_estimate); + if (work_estimate_out != nullptr) { *work_estimate_out = work_estimate; } + return combinations; +} + +template +bool generate_mod2_zero_half_cuts(cut_pool_t& cut_pool, + const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + csr_matrix_t& Arow, + const std::vector& new_slacks, + const std::vector& var_types, + const std::vector& xstar, + variable_bounds_t& variable_bounds, + f_t start_time, + f_t& work_estimate) +{ + constexpr i_t max_combination_size = 64; + constexpr i_t max_row_combinations = 1000; + constexpr f_t min_violation = (f_t)1e-6; + constexpr f_t candidate_work_limit = (f_t)3e7; + constexpr f_t combination_work_limit = (f_t)3e7; + constexpr f_t generation_work_limit = (f_t)4e7; + f_t candidate_work = 0.0; + f_t combination_work = 0.0; + f_t generation_work = 0.0; + bool candidate_limit_reached = false; + bool generation_limit_reached = false; + + if (add_work_estimate((f_t)(3 * lp.num_cols) + (f_t)(variable_bounds.upper_variables.size() + + variable_bounds.lower_variables.size()), + &candidate_work, + candidate_work_limit, + &candidate_limit_reached)) { + work_estimate = candidate_work; + return false; + } + complemented_mixed_integer_rounding_cut_t complemented_mir(lp, settings, new_slacks); + std::vector transformed_xstar; + complemented_mir.bound_substitution( + lp, variable_bounds, var_types, xstar, transformed_xstar, true); + + auto candidates = mod2_collect_candidates(complemented_mir, + lp, + Arow, + variable_bounds, + var_types, + transformed_xstar, + start_time, + settings.time_limit, + candidate_work, + candidate_work_limit, + candidate_limit_reached); + + if (toc(start_time) >= settings.time_limit) { + work_estimate = candidate_work; + return true; + } + auto row_combinations = find_mod2_row_combinations(candidates, + max_combination_size, + max_row_combinations, + &combination_work, + combination_work_limit, + start_time, + settings.time_limit); + if (add_work_estimate((f_t)(2 * lp.num_cols), + &generation_work, + generation_work_limit, + &generation_limit_reached)) { + work_estimate = candidate_work + combination_work + generation_work; + return true; + } + scratch_pad_t aggregate_pad(lp.num_cols); + + for (const auto& combination : row_combinations) { + if (toc(start_time) >= settings.time_limit || generation_limit_reached) { break; } + + size_t aggregate_input_nz = 0; + for (const i_t candidate_index : combination) { + aggregate_input_nz += candidates[candidate_index].transformed_inequality.size(); + } + if (add_work_estimate((f_t)(2 * aggregate_input_nz + 1), + &generation_work, + generation_work_limit, + &generation_limit_reached)) { + break; + } + + inequality_t aggregate(lp.num_cols); + bool reversible = true; + for (const i_t candidate_index : combination) { + const auto& candidate = candidates[candidate_index]; + aggregate.rhs += candidate.transformed_inequality.rhs; + reversible = reversible && candidate.reversible; + for (i_t k = 0; k < (i_t)candidate.transformed_inequality.size(); ++k) { + aggregate_pad.add_to_pad(candidate.transformed_inequality.index(k), + candidate.transformed_inequality.coeff(k)); + } + } + aggregate_pad.get_pad(aggregate.vector.i, aggregate.vector.x); + aggregate_pad.clear_pad(); + const f_t aggregate_output_work = + (f_t)(3 * aggregate.size()) + + (f_t)aggregate.size() * std::log2((f_t)aggregate.size() + (f_t)1.0); + if (add_work_estimate(aggregate_output_work, + &generation_work, + generation_work_limit, + &generation_limit_reached)) { + break; + } + aggregate.sort(); + aggregate.scale((f_t)0.5); + + i_t cuts_added = 0; + mod2_generate_cuts_from_aggregate(complemented_mir, + cut_pool, + lp, + settings, + Arow, + variable_bounds, + var_types, + xstar, + transformed_xstar, + aggregate, + min_violation, + start_time, + generation_work, + generation_work_limit, + generation_limit_reached, + cuts_added); + // if the final inequality is reversable, try the reversed version as well + if (reversible && toc(start_time) < settings.time_limit && !generation_limit_reached) { + aggregate.negate(); + mod2_generate_cuts_from_aggregate(complemented_mir, + cut_pool, + lp, + settings, + Arow, + variable_bounds, + var_types, + xstar, + transformed_xstar, + aggregate, + min_violation, + start_time, + generation_work, + generation_work_limit, + generation_limit_reached, + cuts_added); + } + } + work_estimate = candidate_work + combination_work + generation_work; + return true; +} + +template +bool complemented_mixed_integer_rounding_cut_t::generate_lifted_mixed_binary_cover( + const inequality_t& transformed_inequality, + const std::vector& var_types, + const std::vector& transformed_xstar, + inequality_t& transformed_cut, + f_t& work_estimate, + f_t max_work_estimate) +{ + constexpr f_t tolerance = (f_t)1e-6; + + const f_t estimated_work = + (f_t)(12 * transformed_inequality.size()) + + (f_t)transformed_inequality.size() * std::log2((f_t)transformed_inequality.size() + (f_t)1.0); + if (add_work_estimate(estimated_work, &work_estimate, max_work_estimate)) { return false; } + + inequality_t base = transformed_inequality; + base.negate(); + + std::vector locally_complemented(base.size(), 0); + std::vector solution_value(base.size(), 0.0); + std::vector is_integral(base.size(), 0); + for (i_t k = 0; k < (i_t)base.size(); ++k) { + const i_t j = base.index(k); + f_t aj = base.coeff(k); + if (var_types[j] == variable_type_t::CONTINUOUS) { + solution_value[k] = transformed_xstar[j]; + if (aj > 0.0) { base.vector.x[k] = 0.0; } + continue; + } + + const f_t upper = new_upper(j); + if (upper == inf || std::abs(upper - (f_t)1.0) > tolerance) { return false; } + is_integral[k] = 1; + if (aj < 0.0) { + base.rhs -= aj * upper; + base.vector.x[k] = -aj; + solution_value[k] = upper - transformed_xstar[j]; + locally_complemented[k] = 1; + } else { + solution_value[k] = transformed_xstar[j]; + } + } + + std::vector cover; + cover.reserve(base.size()); + for (i_t k = 0; k < (i_t)base.size(); ++k) { + if (is_integral[k] && base.coeff(k) > tolerance && solution_value[k] > tolerance) { + cover.push_back(k); + } + } + if (cover.empty()) { return false; } + + std::stable_sort( + cover.begin(), cover.end(), lifted_cover_order_t{solution_value, base, tolerance}); + + f_t cover_weight = 0.0; + size_t cover_size = 0; + for (; cover_size < cover.size(); ++cover_size) { + cover_weight += base.coeff(cover[cover_size]); + if (cover_weight - base.rhs > tolerance * std::max((f_t)1.0, std::abs(base.rhs))) { + ++cover_size; + break; + } + } + if (cover_size == 0 || cover_size > cover.size()) { return false; } + cover.resize(cover_size); + + const f_t lambda = cover_weight - base.rhs; + if (lambda <= tolerance) { return false; } + std::sort( + cover.begin(), cover.end(), [&](i_t a, i_t b) { return base.coeff(a) > base.coeff(b); }); + + std::vector prefix(cover.size(), 0.0); + std::vector in_cover(base.size(), 0); + f_t prefix_sum = 0.0; + size_t p = cover.size(); + for (size_t h = 0; h < cover.size(); ++h) { + const i_t k = cover[h]; + in_cover[k] = 1; + if (base.coeff(k) - lambda <= tolerance && p == cover.size()) { p = h; } + if (h < p) { + prefix_sum += base.coeff(k); + prefix[h] = prefix_sum; + } + } + if (p == 0) { return false; } + + size_t non_cover_count = 0; + for (i_t k = 0; k < (i_t)base.size(); ++k) { + if (is_integral[k] && !in_cover[k]) { non_cover_count++; } + } + if (add_work_estimate((f_t)(non_cover_count * p), &work_estimate, max_work_estimate)) { + return false; + } + + transformed_cut = base; + transformed_cut.rhs = -lambda; + for (i_t k = 0; k < (i_t)base.size(); ++k) { + if (!is_integral[k]) { + if (base.coeff(k) >= 0.0) { transformed_cut.vector.x[k] = 0.0; } + continue; + } + if (in_cover[k]) { + transformed_cut.vector.x[k] = std::min(base.coeff(k), lambda); + transformed_cut.rhs += transformed_cut.coeff(k); + } else { + transformed_cut.vector.x[k] = + lifted_cover_coefficient(base.coeff(k), prefix, p, lambda, tolerance); + } + } + + for (i_t k = 0; k < (i_t)transformed_cut.size(); ++k) { + if (!locally_complemented[k]) { continue; } + const i_t j = transformed_cut.index(k); + const f_t coefficient = transformed_cut.coeff(k); + transformed_cut.rhs -= coefficient * new_upper(j); + transformed_cut.vector.x[k] = -coefficient; + } + inequality_t squeezed_cut(transformed_cut.vector.n); + transformed_cut.squeeze(squeezed_cut); + transformed_cut = std::move(squeezed_cut); + transformed_cut.negate(); + return true; +} + +#ifdef DUAL_SIMPLEX_INSTANTIATE_DOUBLE +template bool generate_mod2_zero_half_cuts( + cut_pool_t& cut_pool, + const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + csr_matrix_t& Arow, + const std::vector& new_slacks, + const std::vector& var_types, + const std::vector& xstar, + variable_bounds_t& variable_bounds, + double start_time, + double& work_estimate); + +template bool +complemented_mixed_integer_rounding_cut_t::generate_lifted_mixed_binary_cover( + const inequality_t& transformed_inequality, + const std::vector& var_types, + const std::vector& transformed_xstar, + inequality_t& transformed_cut, + double& work_estimate, + double max_work_estimate); +#endif + +} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/tests/mip/cuts_test.cu b/cpp/tests/mip/cuts_test.cu index b4fc3e8cc7..5af0754e3d 100644 --- a/cpp/tests/mip/cuts_test.cu +++ b/cpp/tests/mip/cuts_test.cu @@ -422,7 +422,7 @@ void disable_non_clique_cuts(mip_solver_settings_t& settings) void disable_non_zero_half_cuts(mip_solver_settings_t& settings) { - settings.clique_cuts = 1; + settings.clique_cuts = 0; settings.zero_half_cuts = 1; settings.max_cut_passes = 10; settings.mixed_integer_gomory_cuts = 0; @@ -1482,6 +1482,57 @@ TEST(cuts, zero_half_unit_separator_simple_pentagon) EXPECT_TRUE(found); } +TEST(cuts, zero_half_unit_mod2_row_finder_single_pair_and_four_row_dependencies) +{ + // Empty parity with odd rhs is a one-row zero-half aggregation. + { + const std::vector> parity_rows = {{}}; + const std::vector rhs_parity = {1}; + const auto combinations = + mip::find_mod2_row_combinations_for_test(parity_rows, rhs_parity, 8, 8); + ASSERT_EQ(combinations.size(), 1); + EXPECT_EQ(combinations.front(), std::vector{0}); + } + + // Equal parity and opposite rhs form a two-row dependency. + { + const std::vector> parity_rows = {{0, 2}, {0, 2}}; + const std::vector rhs_parity = {0, 1}; + const auto combinations = + mip::find_mod2_row_combinations_for_test(parity_rows, rhs_parity, 8, 8); + ASSERT_EQ(combinations.size(), 1); + EXPECT_EQ(combinations.front(), (std::vector{0, 1})); + } + + // Four edges of an even cycle cancel in GF(2); the odd aggregate rhs makes + // the dependency eligible for a zero-half cut. + { + const std::vector> parity_rows = {{0, 1}, {1, 2}, {2, 3}, {0, 3}}; + const std::vector rhs_parity = {1, 0, 0, 0}; + const auto combinations = + mip::find_mod2_row_combinations_for_test(parity_rows, rhs_parity, 8, 8); + ASSERT_EQ(combinations.size(), 1); + EXPECT_EQ(combinations.front(), (std::vector{0, 1, 2, 3})); + } +} + +TEST(cuts, zero_half_unit_mod2_row_finder_stops_at_work_limit) +{ + std::vector support(64); + std::iota(support.begin(), support.end(), 0); + const std::vector> parity_rows(256, support); + const std::vector rhs_parity(256, 0); + + constexpr double max_work = 18900.0; + double work = 0.0; + const auto combinations = + mip::find_mod2_row_combinations_for_test(parity_rows, rhs_parity, 64, 1000, max_work, &work); + + EXPECT_TRUE(combinations.empty()); + EXPECT_GT(work, max_work); + EXPECT_LT(work, 19100.0); +} + TEST(cuts, zero_half_unit_separator_no_cycle_for_4_cycle) { // Even cycle: 0-1-2-3-0 @@ -1615,6 +1666,28 @@ TEST(cuts, zero_half_end_to_end_pentagon_tightens_lp_relaxation) EXPECT_NEAR(mip_solution.get_objective_value(), -2.0, kCliqueTestTol); } +TEST(cuts, zero_half_end_to_end_general_row_parity_closes_triangle_root_gap) +{ + const raft::handle_t handle{}; + auto mip_problem = create_pairwise_triangle_set_packing_problem(); + + mip_solver_settings_t settings; + settings.time_limit = 10.0; + settings.presolver = presolver_t::None; + settings.node_limit = 0; + disable_non_zero_half_cuts(settings); + + benchmark_info_t benchmark_info; + settings.benchmark_info_ptr = &benchmark_info; + auto mip_solution = solve_mip(&handle, mip_problem, settings); + + EXPECT_NE(mip_solution.get_termination_status(), mip_termination_status_t::Infeasible); + ASSERT_FALSE(std::isnan(benchmark_info.root_lp_no_cuts)); + ASSERT_FALSE(std::isnan(benchmark_info.root_lp_with_cuts)); + EXPECT_NEAR(benchmark_info.root_lp_no_cuts, -1.5, kCliqueTestTol); + EXPECT_NEAR(benchmark_info.root_lp_with_cuts, -1.0, kCliqueTestTol); +} + TEST(cuts, zero_half_unit_separator_seven_cycle_violated_below_half) { // 7-cycle: 0-1-2-3-4-5-6-0, all weights 0.4. Each edge weight = (1-0.4-0.4)/2 = 0.1 diff --git a/skills/cuopt-developer/SKILL.md b/skills/cuopt-developer/SKILL.md index aa488e064b..4aae65c8be 100644 --- a/skills/cuopt-developer/SKILL.md +++ b/skills/cuopt-developer/SKILL.md @@ -224,7 +224,7 @@ For pre-commit setup, DCO sign-off (`git commit -s`), the fork-based PR workflow ## Coding Conventions -For C++ naming (`snake_case`, `d_`/`h_` prefixes, `_t` suffix), file extensions (`.hpp`/`.cpp`/`.cu`/`.cuh` and which compiler each uses), include order, Python style, error handling (`CUOPT_EXPECTS`, `RAFT_CUDA_TRY`), memory management (RMM patterns, no raw `new`/`delete`), test-impact rules, and volatile-comment rules (hardware names and self-referential issue/PR numbers in comments or skip messages go stale; issue links to a separate tracking issue are fine), see [references/conventions.md](references/conventions.md). +For C++ naming (`snake_case`, `d_`/`h_` prefixes, `_t` suffix), file extensions (`.hpp`/`.cpp`/`.cu`/`.cuh` and which compiler each uses), include order, Python style, error handling (`CUOPT_EXPECTS`, `RAFT_CUDA_TRY`), memory management (RMM patterns, no raw `new`/`delete`), test-impact rules, volatile-comment rules (hardware names and self-referential issue/PR numbers in comments or skip messages go stale; issue links to a separate tracking issue are fine), **no large local lambdas** (extract named helpers instead), and **coarse work-estimate / time-limit gating** (phase/outer-loop only; no fine inner-loop or double checks), see [references/conventions.md](references/conventions.md). ## OpenMP task/runtime compatibility diff --git a/skills/cuopt-developer/references/conventions.md b/skills/cuopt-developer/references/conventions.md index e9963d4824..8ab14eb205 100644 --- a/skills/cuopt-developer/references/conventions.md +++ b/skills/cuopt-developer/references/conventions.md @@ -55,9 +55,20 @@ This applies to all comment types: inline comments, block comments, suppression ## C++ Implementation Style -- Prefer direct loops or named helpers for performance-critical traversal logic. Reserve lambdas for - short predicates and callbacks; large local lambdas obscure control flow and can lead to repeated - scans. +- **Never write large lambdas inside a function body.** If a lambda is more than a short + predicate/comparator (roughly more than ~3–5 lines, or it has nested lambdas, local state, + or non-trivial control flow), extract it as a named free function, file-local helper in an + anonymous namespace, or private method. Large local lambdas obscure control flow, hide reuse, + and make work/time accounting harder to reason about. +- Prefer direct loops or named helpers for performance-critical traversal logic. Reserve + in-function lambdas for short predicates and callbacks only. +- Keep work-estimate and time-limit checks at phase or outer-loop boundaries. Accumulate the work + performed by cheap inner loops and charge it once when the phase completes; do not gate every + inner iteration. Do not separately test a sticky limit or raw estimate immediately before or + after `add_work_estimate` when that call already provides the gate for the same work. +- Charge the operation that is actually performed. For vector copies and sparse traversals, base + work on the number of visited or copied entries rather than only the number of containers. + Avoid charging the same traversal in both its caller and callee. ### Suppression comments (`// NOSONAR`, `// NOLINT`, etc.)