From 77c22ce201d6cce5047a2737110dd321fed016b7 Mon Sep 17 00:00:00 2001 From: akif Date: Wed, 29 Jul 2026 15:21:03 +0000 Subject: [PATCH 01/15] Add general mod-2 zero-half separation Signed-off-by: akif --- cpp/src/cuts/cuts.cpp | 494 ++++++++++++++++++++++++++++++++++++- cpp/src/cuts/cuts.hpp | 27 +- cpp/tests/mip/cuts_test.cu | 58 ++++- 3 files changed, 567 insertions(+), 12 deletions(-) diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index be45ffeecd..e171015519 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -1126,6 +1126,112 @@ std::vector> find_violated_odd_cycles_for_test( 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 +std::vector> find_mod2_row_combinations( + const std::vector>& parity_rows, + const std::vector& rhs_parity, + i_t max_combination_size, + i_t max_combinations, + f_t* work_estimate, + f_t max_work_estimate) +{ + cuopt_assert(parity_rows.size() == rhs_parity.size(), + "GF(2) parity row and rhs sizes must match"); + if (max_combination_size <= 0 || max_combinations <= 0) { return {}; } + + struct basis_row_t { + std::vector parity; + std::vector combination; + bool rhs{false}; + }; + + std::vector permutation(parity_rows.size()); + std::iota(permutation.begin(), permutation.end(), 0); + std::stable_sort(permutation.begin(), permutation.end(), [&](i_t a, i_t b) { + if (parity_rows[a].size() != parity_rows[b].size()) { + return parity_rows[a].size() < parity_rows[b].size(); + } + // Prefer an even-rhs representative for a pivot. If an otherwise + // identical odd-rhs row arrives later, their dependency is immediately a + // valid zero-half aggregation. + return rhs_parity[a] < rhs_parity[b]; + }); + + i_t max_index = -1; + for (const auto& row : parity_rows) { + cuopt_assert(std::is_sorted(row.begin(), row.end()), "GF(2) parity rows must be sorted"); + cuopt_assert(std::adjacent_find(row.begin(), row.end()) == row.end(), + "GF(2) parity rows must not contain duplicates"); + if (!row.empty()) { + cuopt_assert(row.front() >= 0, "GF(2) parity index must be nonnegative"); + max_index = std::max(max_index, row.back()); + } + } + + std::vector pivot_to_basis(static_cast(max_index + 1), -1); + std::vector basis; + basis.reserve(std::min(parity_rows.size(), static_cast(max_index + 1))); + std::vector> combinations; + combinations.reserve(std::min(static_cast(max_combinations), parity_rows.size())); + + std::vector parity_tmp; + std::vector combination_tmp; + for (const i_t candidate : permutation) { + basis_row_t current; + current.parity = parity_rows[candidate]; + current.combination = {candidate}; + current.rhs = rhs_parity[candidate] != 0; + + bool abandoned = false; + while (!current.parity.empty()) { + const i_t pivot = current.parity.front(); + const i_t basis_index = pivot_to_basis[pivot]; + if (basis_index < 0) { break; } + + const auto& pivot_row = basis[basis_index]; + symmetric_difference_sorted(current.parity, pivot_row.parity, parity_tmp); + symmetric_difference_sorted(current.combination, pivot_row.combination, combination_tmp); + if (add_work_estimate( + static_cast(current.parity.size() + pivot_row.parity.size() + + current.combination.size() + pivot_row.combination.size()), + work_estimate, + max_work_estimate) || + combination_tmp.size() > static_cast(max_combination_size)) { + abandoned = true; + break; + } + current.parity.swap(parity_tmp); + current.combination.swap(combination_tmp); + current.rhs = current.rhs != pivot_row.rhs; + } + if (abandoned) { continue; } + + if (current.parity.empty()) { + if (current.rhs && !current.combination.empty()) { + combinations.push_back(std::move(current.combination)); + if (combinations.size() >= static_cast(max_combinations)) { break; } + } + continue; + } + + const i_t pivot = current.parity.front(); + pivot_to_basis[pivot] = static_cast(basis.size()); + basis.push_back(std::move(current)); + if (work_estimate != nullptr && *work_estimate > max_work_estimate) { break; } + } + return combinations; +} + // 64-bit integer mixer (SplitMix64). Used as the building block for the // cousin filter's per-slot independent hash family. inline uint64_t splitmix64_mix(uint64_t x) @@ -1143,6 +1249,21 @@ inline uint64_t hash64_with_seed(uint64_t value, uint64_t seed) } // 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) +{ + double work_estimate = 0.0; + return find_mod2_row_combinations(parity_rows, + rhs_parity, + max_combination_size, + max_combinations, + &work_estimate, + std::numeric_limits::infinity()); +} + template void cut_pool_t::add_cut(cut_type_t cut_type, const inequality_t& cut) { @@ -3604,7 +3725,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) { @@ -3855,9 +3977,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 +4007,219 @@ bool cut_generation_t::generate_zero_half_cuts( static_cast(sub_cg_.ready), sub_cg_.vertices.size()); + // General zero-half separation. The previous implementation only searched + // for odd cycles in the conflict graph. A zero-half separator must also find + // small GF(2) dependencies among tight rows: + // + // sum_i (a_i x >= b_i), with even aggregate integer coefficients and an + // odd aggregate rhs. + // + // Dividing that aggregation by two and applying c-MIR yields a valid + // zero-half cut. The sparse elimination below bounds both combination size + // and work so this path remains predictable on large models. + 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); + + struct mod2_candidate_t { + inequality_t transformed_inequality; + std::vector parity; + bool rhs_parity{false}; + bool reversible{false}; + }; + + constexpr i_t max_integral_scale = 1000; + constexpr i_t max_combination_size = 64; + constexpr i_t max_row_combinations = 1000; + const i_t max_integer_row_length = 1000 + lp.num_cols / 10; + const f_t row_tight_tol = std::max(settings.primal_tol, static_cast(1e-7)); + const f_t coefficient_integral_tol = static_cast(1e-6); + const f_t min_violation = std::max(settings.primal_tol, static_cast(1e-6)); + f_t mod2_work_estimate = 0.0; + const f_t max_mod2_work_estimate = static_cast(1e8); + std::vector mod2_candidates; + mod2_candidates.reserve(lp.num_rows); + + auto integral_scale = [&](const inequality_t& inequality) { + for (i_t scale = 1; scale <= max_integral_scale; ++scale) { + bool integral = true; + const f_t scaled_rhs = static_cast(scale) * inequality.rhs; + if (std::abs(scaled_rhs - std::round(scaled_rhs)) > + coefficient_integral_tol * std::max(static_cast(1.0), std::abs(scaled_rhs))) { + integral = false; + } + for (i_t k = 0; integral && k < static_cast(inequality.size()); ++k) { + 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 = static_cast(scale) * inequality.coeff(k); + if (std::abs(scaled_coefficient - std::round(scaled_coefficient)) > + coefficient_integral_tol * + std::max(static_cast(1.0), std::abs(scaled_coefficient))) { + integral = false; + } + } + if (integral) { return scale; } + } + return i_t{0}; + }; + + for (i_t row = 0; row < lp.num_rows; ++row) { + if (toc(start_time) >= settings.time_limit || mod2_work_estimate > max_mod2_work_estimate) { + break; + } + const i_t slack = complemented_mir.slack_cols(row); + if (slack < 0 || transformed_xstar[slack] > row_tight_tol) { continue; } + + inequality_t inequality(Arow, row, lp.rhs[row]); + if (inequality.size() > static_cast(max_integer_row_length)) { continue; } + complemented_mir.transform_inequality(variable_bounds, var_types, inequality); + inequality.sort(); + mod2_work_estimate += static_cast(4 * inequality.size()); + + // Every LP row is an equality after slack insertion. Choose the direction + // in which the zero-valued transformed slack has a negative coefficient; + // it can then be removed while preserving a valid >= inequality. + i_t slack_position = -1; + for (i_t k = 0; k < static_cast(inequality.size()); ++k) { + if (inequality.index(k) == slack) { + slack_position = k; + break; + } + } + if (slack_position < 0 || inequality.coeff(slack_position) == 0.0) { continue; } + 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); + + // As in general mod-2 separators, only rows whose continuous variables + // are at their selected bounds participate in the parity system. + bool continuous_at_bounds = true; + for (i_t k = 0; k < static_cast(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 = integral_scale(inequality); + if (scale == 0) { continue; } + if (scale != 1) { inequality.scale(static_cast(scale)); } + + mod2_candidate_t candidate; + candidate.transformed_inequality = std::move(inequality); + candidate.rhs_parity = + (std::llabs(static_cast(std::llround(candidate.transformed_inequality.rhs))) % + 2) != 0; + candidate.reversible = std::abs(lp.upper[slack] - lp.lower[slack]) <= row_tight_tol; + for (i_t k = 0; k < static_cast(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 = + static_cast(std::llround(candidate.transformed_inequality.coeff(k))); + if ((std::llabs(coefficient) % 2) != 0) { candidate.parity.push_back(j); } + } + if (candidate.parity.size() > static_cast(max_integer_row_length)) { continue; } + mod2_candidates.push_back(std::move(candidate)); + } + + std::vector> parity_rows; + std::vector rhs_parity; + parity_rows.reserve(mod2_candidates.size()); + rhs_parity.reserve(mod2_candidates.size()); + for (const auto& candidate : mod2_candidates) { + parity_rows.push_back(candidate.parity); + rhs_parity.push_back(candidate.rhs_parity); + } + + auto row_combinations = find_mod2_row_combinations(parity_rows, + rhs_parity, + max_combination_size, + max_row_combinations, + &mod2_work_estimate, + max_mod2_work_estimate); + scratch_pad_t aggregate_pad(lp.num_cols); + i_t mod2_cuts_added = 0; + for (const auto& combination : row_combinations) { + if (toc(start_time) >= settings.time_limit || mod2_work_estimate > max_mod2_work_estimate) { + break; + } + + inequality_t aggregate(lp.num_cols); + bool reversible = true; + for (const i_t candidate_index : combination) { + const auto& candidate = mod2_candidates[candidate_index]; + aggregate.rhs += candidate.transformed_inequality.rhs; + reversible = reversible && candidate.reversible; + for (i_t k = 0; k < static_cast(candidate.transformed_inequality.size()); ++k) { + aggregate_pad.add_to_pad(candidate.transformed_inequality.index(k), + candidate.transformed_inequality.coeff(k)); + } + mod2_work_estimate += static_cast(candidate.transformed_inequality.size()); + } + aggregate_pad.get_pad(aggregate.vector.i, aggregate.vector.x); + aggregate_pad.clear_pad(); + aggregate.sort(); + aggregate.scale(static_cast(0.5)); + + auto generate_from_aggregate = [&](const inequality_t& oriented_aggregate) { + auto add_transformed_cut = [&](inequality_t transformed_cut) { + 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); + complemented_mir.remove_small_coefficients(lp.lower, lp.upper, transformed_cut); + const f_t violation = complemented_mir.compute_violation(transformed_cut, xstar); + mod2_work_estimate += static_cast(10 * transformed_cut.size()); + if (violation > min_violation) { + cut_pool_.add_cut(cut_type_t::ZERO_HALF, transformed_cut); + mod2_cuts_added++; + } + }; + + inequality_t mir_cut(lp.num_cols); + if (complemented_mir.generate_cut_nonnegative_maintain_indicies( + oriented_aggregate, var_types, mir_cut)) { + add_transformed_cut(std::move(mir_cut)); + } + + inequality_t lifted_cover_cut(lp.num_cols); + if (complemented_mir.generate_lifted_mixed_binary_cover(oriented_aggregate, + var_types, + transformed_xstar, + lifted_cover_cut, + mod2_work_estimate)) { + add_transformed_cut(std::move(lifted_cover_cut)); + } + }; + + generate_from_aggregate(aggregate); + if (reversible) { + aggregate.negate(); + generate_from_aggregate(aggregate); + } + } + ZERO_HALF_DEBUG("general mod2 candidates=%lld combinations=%lld cuts=%lld work=%g", + static_cast(mod2_candidates.size()), + static_cast(row_combinations.size()), + static_cast(mod2_cuts_added), + static_cast(mod2_work_estimate)); + // 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,7 +4237,6 @@ 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; // shortest path of length >= 0.5 - min_violation cannot yield a violated cut const f_t cutoff = static_cast(0.5) - min_violation; @@ -5304,7 +5635,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 +5700,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, @@ -5682,6 +6018,146 @@ bool complemented_mixed_integer_rounding_cut_t:: 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) +{ + constexpr f_t tolerance = static_cast(1e-6); + + // Work in <= form. All variables have already been shifted or complemented + // to nonnegative variables by bound_substitution()/transform_inequality(). + 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 < static_cast(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]; + // Positive continuous coefficients can be relaxed from a <= row. + if (aj > 0.0) { base.vector.x[k] = 0.0; } + continue; + } + + const f_t upper = new_upper(j); + // This lifting function is for binary variables. General integer rows + // continue to use the c-MIR zero-half cut. + if (upper == inf || std::abs(upper - static_cast(1.0)) > tolerance) { return false; } + is_integral[k] = 1; + if (aj < 0.0) { + // z = 1 - w makes the knapsack coefficient positive. + 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 < static_cast(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(), [&](i_t a, i_t b) { + 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); + }); + + 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(static_cast(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; } + + auto lifting_function = [&](f_t coefficient) { + for (size_t h = 0; h < p; ++h) { + if (coefficient <= prefix[h] - lambda + tolerance) { return static_cast(h) * lambda; } + if (coefficient <= prefix[h] + tolerance) { + return static_cast(h + 1) * lambda + coefficient - prefix[h]; + } + } + return static_cast(p) * lambda + coefficient - prefix[p - 1]; + }; + + transformed_cut = base; + transformed_cut.rhs = -lambda; + for (i_t k = 0; k < static_cast(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] = lifting_function(base.coeff(k)); + } + } + + // Undo the extra sign-complementations used to obtain a positive knapsack + // row, then return to cuOpt's >= cut convention. + for (i_t k = 0; k < static_cast(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(); + + work_estimate += static_cast(12 * base.size()) + + static_cast(cover.size()) * + std::log2(static_cast(cover.size()) + static_cast(1.0)); + return true; +} + template f_t complemented_mixed_integer_rounding_cut_t::compute_violation( const inequality_t& cut, const std::vector& xstar) diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index 78091c85f6..0f13261a61 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -286,6 +286,16 @@ 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); + template class cut_pool_t { public: @@ -709,12 +719,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 +969,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 +1013,14 @@ 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 compute_violation(const inequality_t& cut, const std::vector& xstar); f_t new_upper(i_t j) const { return transformed_upper_[j]; } diff --git a/cpp/tests/mip/cuts_test.cu b/cpp/tests/mip/cuts_test.cu index b4fc3e8cc7..08b978216d 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,40 @@ 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_separator_no_cycle_for_4_cycle) { // Even cycle: 0-1-2-3-0 @@ -1615,6 +1649,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 From cdd0e19b903791ecb47c0865809cda70c3a86b91 Mon Sep 17 00:00:00 2001 From: akif Date: Tue, 4 Aug 2026 08:34:25 +0000 Subject: [PATCH 02/15] Bound root cut work and report work units --- cpp/src/branch_and_bound/branch_and_bound.cpp | 133 +++++++- cpp/src/branch_and_bound/branch_and_bound.hpp | 1 + cpp/src/cuts/cuts.cpp | 313 ++++++++++++++---- cpp/src/cuts/cuts.hpp | 48 ++- cpp/src/dual_simplex/phase2.cpp | 17 + cpp/src/utilities/work_limit_context.hpp | 5 +- cpp/tests/mip/cuts_test.cu | 36 ++ 7 files changed, 482 insertions(+), 71 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 4dc6bc67a8..2981a759e7 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -69,6 +69,8 @@ using simplex::variable_type_t; namespace { +bool cut_ab_logging_enabled() { return std::getenv("CUOPT_CUT_AB_MODE") != nullptr; } + template bool is_fractional(f_t x, variable_type_t var_type, f_t integer_tol) { @@ -277,6 +279,11 @@ branch_and_bound_t::branch_and_bound_t( solver_status_(mip_status_t::UNSET) { exploration_stats_.start_time = start_time; + work_unit_context_.deterministic = settings_.deterministic; + // Root LP and cut work contributes to the public work limit, but it must not + // enter the outer heuristic/B&B horizon barrier. The heuristic may finish + // before root processing, leaving no matching barrier participant. + work_unit_context_.sync_on_horizon = false; #ifdef PRINT_CONSTRAINT_MATRIX settings_.log.printf("A"); original_problem_.A.print_matrix(); @@ -2794,7 +2801,7 @@ lp_status_t branch_and_bound_t::solve_root_relaxation( nonbasic_list, root_vstatus_, edge_norms_, - nullptr); + &work_unit_context_); } // Wait for the root relaxation solution to be sent by the diversity manager or dual simplex @@ -2957,12 +2964,32 @@ auto branch_and_bound_t::do_cut_pass( nonbasic_list, variable_bounds, exploration_stats_.start_time); + const auto& separator_work = cut_generation.last_work_stats(); + work_unit_context_.record_work_sync_on_horizon(separator_work.work_units()); + settings_.log.debug( + "Cut pass %d separator work %.4f units: gomory %.4f, knapsack %.4f, flow %.4f, " + "MIR %.4f, implied %.4f, graph %.4f, clique %.4f, zero-half %.4f\n", + cut_pass, + separator_work.work_units(), + separator_work.gomory / 1e8, + separator_work.knapsack / 1e8, + separator_work.flow_cover / 1e8, + separator_work.mir / 1e8, + separator_work.implied_bound / 1e8, + separator_work.conflict_graph / 1e8, + separator_work.clique / 1e8, + separator_work.zero_half / 1e8); if (!problem_feasible) { if (settings_.heuristic_preemption_callback != nullptr) { settings_.heuristic_preemption_callback(); } return {cut_pass_action_t::RETURN, mip_status_t::INFEASIBLE}; } + if (work_unit_context_.global_work_units_elapsed >= settings_.work_limit) { + solver_status_ = mip_status_t::WORK_LIMIT; + set_final_solution(solution, root_objective_); + return {cut_pass_action_t::RETURN, solver_status_}; + } if (toc(exploration_stats_.start_time) >= settings_.time_limit) { solver_status_ = mip_status_t::TIME_LIMIT; set_final_solution(solution, root_objective_); @@ -2974,15 +3001,48 @@ auto branch_and_bound_t::do_cut_pass( } // Score the cuts f_t score_start_time = tic(); - cut_pool.score_cuts(root_relax_soln_.x); + constexpr f_t max_cut_scoring_work = static_cast(1e8); + const f_t cut_scoring_work = cut_pool.score_cuts(root_relax_soln_.x, max_cut_scoring_work); + work_unit_context_.record_work_sync_on_horizon(cut_scoring_work / static_cast(1e8)); f_t score_time = toc(score_start_time); if (score_time > 1.0) { settings_.log.debug("Cut scoring time %.2f seconds\n", score_time); } + settings_.log.debug( + "Cut pass %d scoring work %.4f units\n", cut_pass, cut_scoring_work / static_cast(1e8)); + settings_.log.printf( + "Cut pass %d work units: total %.4f, separation %.4f " + "(Gomory %.4f, knapsack %.4f, flow %.4f, MIR %.4f, implied %.4f, graph %.4f, " + "clique %.4f, zero-half %.4f), scoring %.4f\n", + cut_pass, + separator_work.work_units() + cut_scoring_work / static_cast(1e8), + separator_work.work_units(), + separator_work.gomory / 1e8, + separator_work.knapsack / 1e8, + separator_work.flow_cover / 1e8, + separator_work.mir / 1e8, + separator_work.implied_bound / 1e8, + separator_work.conflict_graph / 1e8, + separator_work.clique / 1e8, + separator_work.zero_half / 1e8, + cut_scoring_work / static_cast(1e8)); + if (work_unit_context_.global_work_units_elapsed >= settings_.work_limit) { + solver_status_ = mip_status_t::WORK_LIMIT; + set_final_solution(solution, root_objective_); + return {cut_pass_action_t::RETURN, solver_status_}; + } // Get the best cuts from the cut pool csr_matrix_t cuts_to_add(0, original_lp_.num_cols, 0); std::vector cut_rhs; std::vector cut_types; i_t num_cuts = cut_pool.get_best_cuts(cuts_to_add, cut_rhs, cut_types); if (num_cuts == 0) { return {cut_pass_action_t::BREAK, mip_status_t::UNSET}; } + const f_t cut_assembly_work = + static_cast(5 * cuts_to_add.row_start[cuts_to_add.m] + 4 * num_cuts); + work_unit_context_.record_work_sync_on_horizon(cut_assembly_work / static_cast(1e8)); + if (work_unit_context_.global_work_units_elapsed >= settings_.work_limit) { + solver_status_ = mip_status_t::WORK_LIMIT; + set_final_solution(solution, root_objective_); + return {cut_pass_action_t::RETURN, solver_status_}; + } cut_info.record_cut_types(cut_types); #ifdef PRINT_CUT_POOL_TYPES cut_pool.print_cutpool_types(); @@ -3064,6 +3124,8 @@ auto branch_and_bound_t::do_cut_pass( std::vector new_upper = original_lp_.upper; bool feasible = node_presolve.bounds_strengthening(settings_, bounds_changed, new_lower, new_upper); + work_unit_context_.record_work_sync_on_horizon(node_presolve.last_nnz_processed / + static_cast(1e8)); mutex_original_lp_.lock(); original_lp_.lower = new_lower; original_lp_.upper = new_upper; @@ -3080,6 +3142,12 @@ auto branch_and_bound_t::do_cut_pass( return {cut_pass_action_t::RETURN, mip_status_t::INFEASIBLE}; } + if (work_unit_context_.global_work_units_elapsed >= settings_.work_limit) { + solver_status_ = mip_status_t::WORK_LIMIT; + set_final_solution(solution, root_objective_); + return {cut_pass_action_t::RETURN, solver_status_}; + } + if (toc(exploration_stats_.start_time) >= settings_.time_limit) { solver_status_ = mip_status_t::TIME_LIMIT; set_final_solution(solution, root_objective_); @@ -3102,7 +3170,8 @@ auto branch_and_bound_t::do_cut_pass( nonbasic_list, root_relax_soln_, iter, - edge_norms_); + edge_norms_, + &work_unit_context_); exploration_stats_.total_simplex_iters += iter; f_t dual_phase2_time = toc(dual_phase2_start_time); if (dual_phase2_time > 1.0) { @@ -3113,6 +3182,11 @@ auto branch_and_bound_t::do_cut_pass( set_final_solution(solution, root_objective_); return {cut_pass_action_t::RETURN, solver_status_}; } + if (cut_status == dual_status_t::WORK_LIMIT) { + solver_status_ = mip_status_t::WORK_LIMIT; + set_final_solution(solution, root_objective_); + return {cut_pass_action_t::RETURN, solver_status_}; + } if (cut_status != dual_status_t::OPTIMAL) { settings_.log.printf("Numerical issue at root node. Resolving from scratch\n"); @@ -3125,12 +3199,17 @@ auto branch_and_bound_t::do_cut_pass( basic_list, nonbasic_list, root_vstatus_, - edge_norms_); + edge_norms_, + &work_unit_context_); if (scratch_status == lp_status_t::OPTIMAL) { // We recovered cut_status = convert_lp_status_to_dual_status(scratch_status); exploration_stats_.total_simplex_iters += root_relax_soln_.iterations; root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); + } else if (scratch_status == lp_status_t::WORK_LIMIT) { + solver_status_ = mip_status_t::WORK_LIMIT; + set_final_solution(solution, root_objective_); + return {cut_pass_action_t::RETURN, solver_status_}; } else { settings_.log.printf("Cut status %s\n", simplex::dual_status_to_string(cut_status).c_str()); #ifdef WRITE_CUT_INFEASIBLE_MPS @@ -3301,7 +3380,8 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut basic_list, nonbasic_list, root_vstatus_, - edge_norms_); + edge_norms_, + &work_unit_context_); root_relax_solved_by = DualSimplex; exploration_stats_.total_simplex_iters = root_relax_soln_.iterations; @@ -3367,6 +3447,15 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut return solver_status_; } + if (work_unit_context_.global_work_units_elapsed >= settings_.work_limit) { + settings_.log.printf("\n"); + solver_status_ = mip_status_t::WORK_LIMIT; + set_final_solution(solution, -inf); + signal_extend_cliques_.store(true, std::memory_order_release); +#pragma omp taskwait depend(in : *clique_signal) + return solver_status_; + } + assert(root_status == lp_status_t::OPTIMAL); settings_.log.printf("\n"); settings_.log.print_format("Root relaxation solution found in {} iterations and {:.2f}s by {}\n", @@ -3379,6 +3468,10 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut set_uninitialized_steepest_edge_norms(original_lp_, basic_list, edge_norms_); root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); + if (cut_ab_logging_enabled()) { + settings_.log.printf("CUT_AB_ROOT_NO_CUTS %.17g\n", + compute_user_objective(original_lp_, root_objective_)); + } if (settings_.set_simplex_solution_callback != nullptr) { std::vector original_x; @@ -3401,6 +3494,10 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut cut_info_t cut_info; if (num_fractional == 0) { + if (cut_ab_logging_enabled()) { + settings_.log.printf("CUT_AB_ROOT_WITH_CUTS %.17g\n", + compute_user_objective(original_lp_, root_objective_)); + } if (settings_.benchmark_info_ptr != nullptr) { const double v = static_cast(compute_user_objective(original_lp_, root_objective_)); settings_.benchmark_info_ptr->root_lp_no_cuts = v; @@ -3475,6 +3572,10 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut settings_.benchmark_info_ptr->root_lp_with_cuts = compute_user_objective(original_lp_, root_objective_); } + if (cut_ab_logging_enabled()) { + settings_.log.printf("CUT_AB_ROOT_WITH_CUTS %.17g\n", + compute_user_objective(original_lp_, root_objective_)); + } set_solution_at_root(solution, cut_info); if (settings_.benchmark_info_ptr != nullptr) { settings_.benchmark_info_ptr->cut_generation_time_sec = toc(cut_generation_start_time); @@ -3508,6 +3609,10 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut root_fj_cpu_worker.stop(); if (cut_pass_result.action == cut_pass_action_t::RETURN) { + if (cut_ab_logging_enabled()) { + settings_.log.printf("CUT_AB_ROOT_WITH_CUTS %.17g\n", + compute_user_objective(original_lp_, root_objective_)); + } if (settings_.benchmark_info_ptr != nullptr) { settings_.benchmark_info_ptr->cut_generation_time_sec = toc(cut_generation_start_time); } @@ -3533,6 +3638,10 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut settings_.benchmark_info_ptr->root_lp_with_cuts = compute_user_objective(original_lp_, root_objective_); } + if (cut_ab_logging_enabled()) { + settings_.log.printf("CUT_AB_ROOT_WITH_CUTS %.17g\n", + compute_user_objective(original_lp_, root_objective_)); + } print_cut_info(settings_, cut_info); f_t cut_generation_time = toc(cut_generation_start_time); @@ -3563,6 +3672,16 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut settings_.log.printf("\n"); } + // Benchmark-only root stop: CUT_AB measurements compare separator families + // at the completed root LP, before strong branching or tree exploration. + if (cut_ab_logging_enabled()) { + solver_status_ = mip_status_t::NODE_LIMIT; + set_final_solution(solution, root_objective_); + signal_extend_cliques_.store(true, std::memory_order_release); +#pragma omp taskwait depend(in : *clique_signal) + return solver_status_; + } + if (enable_root_cut_cpufj && cut_info.has_cuts()) { f_t root_cut_cpufj_build_start_time = tic(); // In deterministic mode this CPUFJ is built on the B&B task while the LS deterministic @@ -3702,6 +3821,8 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut } print_table_header(); + deterministic_root_work_offset_ = work_unit_context_.global_work_units_elapsed; + #pragma omp taskgroup { if (settings_.deterministic) { @@ -4189,7 +4310,7 @@ void branch_and_bound_t::deterministic_sync_callback() } // Stop early if next horizon exceeds work limit - if (deterministic_current_horizon_ > settings_.work_limit) { + if (deterministic_root_work_offset_ + deterministic_current_horizon_ > settings_.work_limit) { deterministic_global_termination_status_ = mip_status_t::WORK_LIMIT; } diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index 96b8a6d8fe..a33c6c84b8 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -487,6 +487,7 @@ class branch_and_bound_t { mip_status_t deterministic_global_termination_status_{mip_status_t::UNSET}; double deterministic_horizon_step_{5.0}; // Work unit step per horizon (tunable) double deterministic_current_horizon_{0.0}; // Current horizon target + double deterministic_root_work_offset_{0.0}; bool deterministic_mode_enabled_{false}; int deterministic_horizon_number_{0}; // Current horizon number (for debugging) diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index e171015519..49d99389e1 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -1155,8 +1155,25 @@ std::vector> find_mod2_row_combinations( bool rhs{false}; }; + i_t max_index = -1; + for (const auto& row : parity_rows) { + if (add_work_estimate(static_cast(row.size() + 1), work_estimate, max_work_estimate)) { + return {}; + } + cuopt_assert(std::is_sorted(row.begin(), row.end()), "GF(2) parity rows must be sorted"); + cuopt_assert(std::adjacent_find(row.begin(), row.end()) == row.end(), + "GF(2) parity rows must not contain duplicates"); + if (!row.empty()) { + cuopt_assert(row.front() >= 0, "GF(2) parity index must be nonnegative"); + max_index = std::max(max_index, row.back()); + } + } + std::vector permutation(parity_rows.size()); std::iota(permutation.begin(), permutation.end(), 0); + const f_t sort_work = static_cast(permutation.size()) * + std::log2(static_cast(permutation.size()) + static_cast(1.0)); + if (add_work_estimate(sort_work, work_estimate, max_work_estimate)) { return {}; } std::stable_sort(permutation.begin(), permutation.end(), [&](i_t a, i_t b) { if (parity_rows[a].size() != parity_rows[b].size()) { return parity_rows[a].size() < parity_rows[b].size(); @@ -1167,17 +1184,9 @@ std::vector> find_mod2_row_combinations( return rhs_parity[a] < rhs_parity[b]; }); - i_t max_index = -1; - for (const auto& row : parity_rows) { - cuopt_assert(std::is_sorted(row.begin(), row.end()), "GF(2) parity rows must be sorted"); - cuopt_assert(std::adjacent_find(row.begin(), row.end()) == row.end(), - "GF(2) parity rows must not contain duplicates"); - if (!row.empty()) { - cuopt_assert(row.front() >= 0, "GF(2) parity index must be nonnegative"); - max_index = std::max(max_index, row.back()); - } + if (add_work_estimate(static_cast(max_index + 1), work_estimate, max_work_estimate)) { + return {}; } - std::vector pivot_to_basis(static_cast(max_index + 1), -1); std::vector basis; basis.reserve(std::min(parity_rows.size(), static_cast(max_index + 1))); @@ -1187,6 +1196,10 @@ std::vector> find_mod2_row_combinations( std::vector parity_tmp; std::vector combination_tmp; for (const i_t candidate : permutation) { + if (add_work_estimate( + static_cast(parity_rows[candidate].size() + 2), work_estimate, max_work_estimate)) { + break; + } basis_row_t current; current.parity = parity_rows[candidate]; current.combination = {candidate}; @@ -1199,14 +1212,16 @@ std::vector> find_mod2_row_combinations( if (basis_index < 0) { break; } const auto& pivot_row = basis[basis_index]; - symmetric_difference_sorted(current.parity, pivot_row.parity, parity_tmp); - symmetric_difference_sorted(current.combination, pivot_row.combination, combination_tmp); if (add_work_estimate( static_cast(current.parity.size() + pivot_row.parity.size() + current.combination.size() + pivot_row.combination.size()), work_estimate, - max_work_estimate) || - combination_tmp.size() > static_cast(max_combination_size)) { + max_work_estimate)) { + return combinations; + } + symmetric_difference_sorted(current.parity, pivot_row.parity, parity_tmp); + symmetric_difference_sorted(current.combination, pivot_row.combination, combination_tmp); + if (combination_tmp.size() > static_cast(max_combination_size)) { abandoned = true; break; } @@ -1254,14 +1269,32 @@ std::vector> find_mod2_row_combinations_for_test( 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) { double work_estimate = 0.0; - return find_mod2_row_combinations(parity_rows, - rhs_parity, - max_combination_size, - max_combinations, - &work_estimate, - std::numeric_limits::infinity()); + auto combinations = find_mod2_row_combinations(parity_rows, + rhs_parity, + max_combination_size, + max_combinations, + &work_estimate, + max_work_estimate); + if (work_estimate_out != nullptr) { *work_estimate_out = work_estimate; } + return combinations; } template @@ -1347,10 +1380,21 @@ f_t cut_pool_t::cut_orthogonality(i_t i, i_t j) template void cut_pool_t::check_for_duplicate_cuts() +{ + f_t work_estimate = 0.0; + check_for_duplicate_cuts(work_estimate, std::numeric_limits::infinity()); +} + +template +bool cut_pool_t::check_for_duplicate_cuts(f_t& work_estimate, f_t max_work_estimate) { // 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 + const f_t setup_work = static_cast(5 * cut_storage_.m + 2 * cut_storage_.n) + + static_cast(4 * cut_storage_.row_start[cut_storage_.m]); + if (add_work_estimate(setup_work, &work_estimate, max_work_estimate)) { return false; } + std::vector divisors(cut_storage_.m, 0.0); std::vector sets(cut_storage_.m, 0); @@ -1381,6 +1425,10 @@ void cut_pool_t::check_for_duplicate_cuts() new_rows++; } else if (sets[r] < new_set_0) { // Look over indices a_ij with i > r + if (add_work_estimate( + static_cast(6 * (col_end - (p + 1))), &work_estimate, max_work_estimate)) { + return false; + } 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]; @@ -1421,6 +1469,10 @@ void cut_pool_t::check_for_duplicate_cuts() const i_t set_r = sets[r]; if (set_r > 0 && set_r < sentinel && cuts_to_remove[r] == 0) { // This cut has a duplicate + if (add_work_estimate( + static_cast(5 * (m - (r + 1))), &work_estimate, max_work_estimate)) { + return false; + } for (i_t i = r + 1; i < m; i++) { if (sets[i] == set_r) { const f_t f_r = divisors[r]; @@ -1475,12 +1527,24 @@ void cut_pool_t::check_for_duplicate_cuts() cut_type_.resize(write); cut_age_.resize(write); } + return true; } template -void cut_pool_t::score_cuts(std::vector& x_relax) +f_t cut_pool_t::score_cuts(std::vector& x_relax, f_t max_work_estimate) { - check_for_duplicate_cuts(); + f_t work_estimate = 0.0; + best_cuts_.clear(); + scored_cuts_ = 0; + + const f_t duplicate_work_limit = std::isfinite(max_work_estimate) + ? max_work_estimate / static_cast(2.0) + : max_work_estimate; + check_for_duplicate_cuts(work_estimate, duplicate_work_limit); + + const f_t distance_work = static_cast(5 * cut_storage_.row_start[cut_storage_.m]) + + static_cast(3 * cut_storage_.m); + if (add_work_estimate(distance_work, &work_estimate, max_work_estimate)) { return work_estimate; } cut_distances_.resize(cut_storage_.m, 0.0); cut_norms_.resize(cut_storage_.m, 0.0); @@ -1500,13 +1564,14 @@ void cut_pool_t::score_cuts(std::vector& x_relax) } std::vector sorted_indices; + const f_t sort_work = static_cast(cut_storage_.m) * + std::log2(static_cast(cut_storage_.m) + static_cast(1.0)); + if (add_work_estimate(sort_work, &work_estimate, max_work_estimate)) { return work_estimate; } best_score_last_permutation(cut_distances_, sorted_indices); const i_t max_cuts = 2000; const f_t min_orthogonality = settings_.cut_min_orthogonality; best_cuts_.reserve(std::min(max_cuts, cut_storage_.m)); - best_cuts_.clear(); - scored_cuts_ = 0; if (!sorted_indices.empty()) { const i_t i = sorted_indices.back(); @@ -1515,7 +1580,8 @@ void cut_pool_t::score_cuts(std::vector& x_relax) scored_cuts_++; } - while (scored_cuts_ < max_cuts && !sorted_indices.empty()) { + bool work_limit_reached = false; + while (scored_cuts_ < max_cuts && !sorted_indices.empty() && !work_limit_reached) { const i_t i = sorted_indices.back(); sorted_indices.pop_back(); @@ -1525,13 +1591,21 @@ void cut_pool_t::score_cuts(std::vector& x_relax) const i_t best_cuts_size = best_cuts_.size(); for (i_t k = 0; k < best_cuts_size; k++) { const i_t j = best_cuts_[k]; + const i_t i_nz = cut_storage_.row_start[i + 1] - cut_storage_.row_start[i]; + const i_t j_nz = cut_storage_.row_start[j + 1] - cut_storage_.row_start[j]; + if (add_work_estimate( + static_cast(4 * (i_nz + j_nz)), &work_estimate, max_work_estimate)) { + work_limit_reached = true; + break; + } cut_ortho = std::min(cut_ortho, cut_orthogonality(i, j)); } - if (cut_ortho >= min_orthogonality) { + if (!work_limit_reached && cut_ortho >= min_orthogonality) { best_cuts_.push_back(i); scored_cuts_++; } } + return work_estimate; } template @@ -3216,6 +3290,12 @@ void cut_generation_t::generate_implied_bound_cuts( { if (probing_implied_bound_.zero_offsets.empty()) { return; } + last_work_stats_.implied_bound += + static_cast( + 4 * std::min(lp.num_cols, static_cast(probing_implied_bound_.zero_offsets.size()) - 1)) + + static_cast(20 * (probing_implied_bound_.zero_variables.size() + + probing_implied_bound_.one_variables.size())); + 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; @@ -3386,8 +3466,8 @@ void cut_generation_t::prepare_fractional_sub_conflict_graph( } const f_t bound_tol = settings.primal_tol; - f_t work_estimate = 0.0; - const f_t max_work_estimate = 1e7; + f_t& work_estimate = last_work_stats_.conflict_graph; + const f_t max_work_estimate = work_estimate + static_cast(1e7); sub_cg_.num_vars = num_vars; sub_cg_.vertices.reserve(static_cast(num_vars) * 2); @@ -3631,6 +3711,8 @@ bool cut_generation_t::generate_cuts(const lp_problem_t& lp, variable_bounds_t& variable_bounds, f_t start_time) { + last_work_stats_ = {}; + // Generate Gomory and CG Cuts if (settings.mixed_integer_gomory_cuts != 0 || settings.strong_chvatal_gomory_cuts != 0) { if (toc(start_time) >= settings.time_limit) { return true; } @@ -3757,6 +3839,9 @@ void cut_generation_t::generate_knapsack_cuts( if (knapsack_generation_.num_knapsack_constraints() > 0) { for (i_t knapsack_row : knapsack_generation_.get_knapsack_constraints()) { if (toc(start_time) >= settings.time_limit) { return; } + const f_t row_size = static_cast(Arow.row_length(knapsack_row)); + last_work_stats_.knapsack += + static_cast(20.0) * row_size + row_size * std::log2(row_size + static_cast(1.0)); 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); @@ -3776,8 +3861,12 @@ void cut_generation_t::generate_flow_cover_cuts( f_t start_time) { if (flow_cover_generation_.num_constraints() > 0) { + last_work_stats_.flow_cover += static_cast(4 * lp.num_cols); for (const auto& flow_cover_row : flow_cover_generation_.get_constraints()) { if (toc(start_time) >= settings.time_limit) { return; } + const f_t row_size = static_cast(Arow.row_length(flow_cover_row.row)); + last_work_stats_.flow_cover += + static_cast(30.0) * row_size + row_size * std::log2(row_size + static_cast(1.0)); inequality_t cut(lp.num_cols); i_t status = flow_cover_generation_.generate_cut( lp, settings, Arow, variable_bounds, var_types, xstar, flow_cover_row, cut); @@ -3826,8 +3915,8 @@ bool cut_generation_t::generate_clique_cuts( const f_t min_weight = 1.0 + min_violation; // TODO this can be problem dependent const i_t max_calls = 100000; - f_t work_estimate = 0.0; - const f_t max_work_estimate = 1e8; + f_t& work_estimate = last_work_stats_.clique; + const f_t max_work_estimate = work_estimate + static_cast(1e8); const std::vector& vertices = sub_cg_.vertices; const std::vector& weights = sub_cg_.weights; @@ -4017,6 +4106,28 @@ bool cut_generation_t::generate_zero_half_cuts( // Dividing that aggregation by two and applying c-MIR yields a valid // zero-half cut. The sparse elimination below bounds both combination size // and work so this path remains predictable on large models. + constexpr i_t max_integral_scale = 1000; + constexpr i_t max_combination_size = 64; + constexpr i_t max_row_combinations = 1000; + const i_t max_integer_row_length = 1000 + lp.num_cols / 10; + const f_t row_tight_tol = std::max(settings.primal_tol, static_cast(1e-7)); + const f_t coefficient_integral_tol = static_cast(1e-6); + const f_t min_violation = std::max(settings.primal_tol, static_cast(1e-6)); + f_t& mod2_work_estimate = last_work_stats_.zero_half; + const f_t max_mod2_work_estimate = mod2_work_estimate + static_cast(1e8); + bool mod2_work_limit_reached = false; + + auto charge_mod2_work = [&](f_t work) { + const bool limit = add_work_estimate(work, &mod2_work_estimate, max_mod2_work_estimate); + mod2_work_limit_reached = mod2_work_limit_reached || limit; + return limit; + }; + + if (charge_mod2_work(static_cast(3 * lp.num_cols) + + static_cast(variable_bounds.upper_variables.size() + + variable_bounds.lower_variables.size()))) { + return true; + } complemented_mixed_integer_rounding_cut_t complemented_mir(lp, settings, new_slacks); std::vector transformed_xstar; complemented_mir.bound_substitution( @@ -4029,20 +4140,15 @@ bool cut_generation_t::generate_zero_half_cuts( bool reversible{false}; }; - constexpr i_t max_integral_scale = 1000; - constexpr i_t max_combination_size = 64; - constexpr i_t max_row_combinations = 1000; - const i_t max_integer_row_length = 1000 + lp.num_cols / 10; - const f_t row_tight_tol = std::max(settings.primal_tol, static_cast(1e-7)); - const f_t coefficient_integral_tol = static_cast(1e-6); - const f_t min_violation = std::max(settings.primal_tol, static_cast(1e-6)); - f_t mod2_work_estimate = 0.0; - const f_t max_mod2_work_estimate = static_cast(1e8); std::vector mod2_candidates; mod2_candidates.reserve(lp.num_rows); auto integral_scale = [&](const inequality_t& inequality) { for (i_t scale = 1; scale <= max_integral_scale; ++scale) { + if (toc(start_time) >= settings.time_limit || + charge_mod2_work(static_cast(inequality.size() + 1))) { + return i_t{0}; + } bool integral = true; const f_t scaled_rhs = static_cast(scale) * inequality.rhs; if (std::abs(scaled_rhs - std::round(scaled_rhs)) > @@ -4067,17 +4173,22 @@ bool cut_generation_t::generate_zero_half_cuts( }; for (i_t row = 0; row < lp.num_rows; ++row) { - if (toc(start_time) >= settings.time_limit || mod2_work_estimate > max_mod2_work_estimate) { + if (toc(start_time) >= settings.time_limit || mod2_work_limit_reached || + mod2_work_estimate > max_mod2_work_estimate) { 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 = static_cast(8 * row_length + 5) + + static_cast(row_length) * + std::log2(static_cast(row_length) + static_cast(1.0)); + if (charge_mod2_work(row_work)) { break; } inequality_t inequality(Arow, row, lp.rhs[row]); - if (inequality.size() > static_cast(max_integer_row_length)) { continue; } complemented_mir.transform_inequality(variable_bounds, var_types, inequality); inequality.sort(); - mod2_work_estimate += static_cast(4 * inequality.size()); // Every LP row is an equality after slack insertion. Choose the direction // in which the zero-valued transformed slack has a negative coefficient; @@ -4138,23 +4249,37 @@ bool cut_generation_t::generate_zero_half_cuts( parity_rows.reserve(mod2_candidates.size()); rhs_parity.reserve(mod2_candidates.size()); for (const auto& candidate : mod2_candidates) { + if (charge_mod2_work(static_cast(candidate.parity.size() + 1))) { break; } parity_rows.push_back(candidate.parity); rhs_parity.push_back(candidate.rhs_parity); } - auto row_combinations = find_mod2_row_combinations(parity_rows, + auto row_combinations = find_mod2_row_combinations(parity_rows, rhs_parity, max_combination_size, max_row_combinations, &mod2_work_estimate, max_mod2_work_estimate); + mod2_work_limit_reached = mod2_work_limit_reached || mod2_work_estimate > max_mod2_work_estimate; scratch_pad_t aggregate_pad(lp.num_cols); i_t mod2_cuts_added = 0; + for (const auto& combination : row_combinations) { - if (toc(start_time) >= settings.time_limit || mod2_work_estimate > max_mod2_work_estimate) { + if (toc(start_time) >= settings.time_limit || mod2_work_limit_reached || + mod2_work_estimate > max_mod2_work_estimate) { break; } + size_t aggregate_input_nz = 0; + for (const i_t candidate_index : combination) { + aggregate_input_nz += mod2_candidates[candidate_index].transformed_inequality.size(); + } + const f_t aggregate_work = + static_cast(4 * aggregate_input_nz + 1) + + static_cast(aggregate_input_nz) * + std::log2(static_cast(aggregate_input_nz) + static_cast(1.0)); + if (charge_mod2_work(aggregate_work)) { break; } + inequality_t aggregate(lp.num_cols); bool reversible = true; for (const i_t candidate_index : combination) { @@ -4165,7 +4290,6 @@ bool cut_generation_t::generate_zero_half_cuts( aggregate_pad.add_to_pad(candidate.transformed_inequality.index(k), candidate.transformed_inequality.coeff(k)); } - mod2_work_estimate += static_cast(candidate.transformed_inequality.size()); } aggregate_pad.get_pad(aggregate.vector.i, aggregate.vector.x); aggregate_pad.clear_pad(); @@ -4174,18 +4298,25 @@ bool cut_generation_t::generate_zero_half_cuts( auto generate_from_aggregate = [&](const inequality_t& oriented_aggregate) { auto add_transformed_cut = [&](inequality_t transformed_cut) { + if (toc(start_time) >= settings.time_limit || + charge_mod2_work(static_cast(10 * transformed_cut.size() + 1))) { + return; + } 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); complemented_mir.remove_small_coefficients(lp.lower, lp.upper, transformed_cut); const f_t violation = complemented_mir.compute_violation(transformed_cut, xstar); - mod2_work_estimate += static_cast(10 * transformed_cut.size()); if (violation > min_violation) { cut_pool_.add_cut(cut_type_t::ZERO_HALF, transformed_cut); mod2_cuts_added++; } }; + if (toc(start_time) >= settings.time_limit || + charge_mod2_work(static_cast(3 * oriented_aggregate.size() + 1))) { + return; + } inequality_t mir_cut(lp.num_cols); if (complemented_mir.generate_cut_nonnegative_maintain_indicies( oriented_aggregate, var_types, mir_cut)) { @@ -4193,13 +4324,17 @@ bool cut_generation_t::generate_zero_half_cuts( } inequality_t lifted_cover_cut(lp.num_cols); - if (complemented_mir.generate_lifted_mixed_binary_cover(oriented_aggregate, + if (toc(start_time) < settings.time_limit && + complemented_mir.generate_lifted_mixed_binary_cover(oriented_aggregate, var_types, transformed_xstar, lifted_cover_cut, - mod2_work_estimate)) { + mod2_work_estimate, + max_mod2_work_estimate)) { add_transformed_cut(std::move(lifted_cover_cut)); } + mod2_work_limit_reached = + mod2_work_limit_reached || mod2_work_estimate > max_mod2_work_estimate; }; generate_from_aggregate(aggregate); @@ -4237,11 +4372,11 @@ 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 bound_tol = settings.primal_tol; + 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; - const f_t max_work_estimate = 1e8; + f_t& work_estimate = last_work_stats_.zero_half; + const f_t max_work_estimate = work_estimate + static_cast(1e8); const std::vector& vertices = sub_cg_.vertices; const std::vector& weights = sub_cg_.weights; @@ -4407,7 +4542,10 @@ void cut_generation_t::generate_mir_cuts( 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; + f_t& work_estimate = last_work_stats_.mir; + work_estimate += static_cast(4 * lp.num_cols + 3 * lp.num_rows) + + static_cast(4 * Arow.row_start[Arow.m]); + const f_t max_work_estimate = work_estimate + static_cast(2e9); i_t num_cuts = 0; while (num_cuts < max_cuts && !score_queue.empty()) { if (toc(start_time) >= settings.time_limit) { break; } @@ -4426,7 +4564,7 @@ void cut_generation_t::generate_mir_cuts( const f_t slack_value = xstar[slack]; if (max_score <= 0.0) { break; } - if (work_estimate > 2e9) { break; } + if (work_estimate > max_work_estimate) { break; } inequality_t inequality(Arow, i, lp.rhs[i]); work_estimate += inequality.size(); @@ -4658,6 +4796,14 @@ void cut_generation_t::generate_gomory_cuts( const std::vector& nonbasic_list, f_t start_time) { + f_t& work_estimate = last_work_stats_.gomory; + const f_t max_work_estimate = work_estimate + static_cast(1e8); + const f_t predicted_setup_work = + static_cast(2.0) * static_cast(lp.num_rows) * static_cast(lp.num_rows) + + static_cast(12 * lp.num_cols) + static_cast(5 * Arow.row_start[Arow.m]); + if (add_work_estimate(predicted_setup_work, &work_estimate, max_work_estimate)) { return; } + + const f_t basis_setup_start_work = basis_update.work_estimate(); tableau_equality_t tableau(lp, basis_update, nonbasic_list); mixed_integer_gomory_cut_t gomory_cut; complemented_mixed_integer_rounding_cut_t complemented_mir(lp, settings, new_slacks); @@ -4668,17 +4814,38 @@ void cut_generation_t::generate_gomory_cuts( std::vector transformed_xstar; complemented_mir.bound_substitution(lp, variable_bounds, var_types, xstar, transformed_xstar); + const f_t actual_setup_work = basis_update.work_estimate() - basis_setup_start_work; + if (actual_setup_work > predicted_setup_work) { + work_estimate += actual_setup_work - predicted_setup_work; + } + work_estimate += + static_cast(4 * lp.num_rows) + static_cast(variable_bounds.upper_variables.size() + + variable_bounds.lower_variables.size()); + if (work_estimate > max_work_estimate) { return; } + for (i_t i = 0; i < lp.num_rows; i++) { - if (toc(start_time) >= settings.time_limit) { break; } + if (toc(start_time) >= settings.time_limit || work_estimate > max_work_estimate) { break; } inequality_t inequality(lp.num_cols); const i_t j = basic_list[i]; if (var_types[j] != variable_type_t::INTEGER) { continue; } const f_t x_j = xstar[j]; if (fractional_part(x_j) < 0.05 || fractional_part(x_j) > 0.95) { continue; } - i_t tableau_status = tableau.generate_base_equality( - lp, settings, Arow, var_types, basis_update, xstar, basic_list, nonbasic_list, i, inequality); + i_t tableau_status = tableau.generate_base_equality(lp, + settings, + Arow, + var_types, + basis_update, + xstar, + basic_list, + nonbasic_list, + i, + inequality, + work_estimate, + max_work_estimate); if (tableau_status == 0) { + const f_t cut_work = static_cast(120 * inequality.size() + 10); + if (add_work_estimate(cut_work, &work_estimate, max_work_estimate)) { break; } // Generate a CG cut const bool generate_cg_cut = settings.strong_chvatal_gomory_cuts != 0; if (generate_cg_cut) { @@ -4788,7 +4955,9 @@ i_t tableau_equality_t::generate_base_equality( const std::vector& basic_list, const std::vector& nonbasic_list, i_t i, - inequality_t& inequality) + inequality_t& inequality, + f_t& work_estimate, + f_t max_work_estimate) { // Let's look for Gomory cuts const i_t j = basic_list[i]; @@ -4804,7 +4973,16 @@ i_t tableau_equality_t::generate_base_equality( e_i.i[0] = i; e_i.x[0] = 1.0; sparse_vector_t u_bar(lp.num_rows, 0); + const f_t predicted_basis_work = static_cast(3 * lp.num_rows + 4); + if (add_work_estimate(predicted_basis_work, &work_estimate, max_work_estimate)) { return -2; } + const f_t basis_work_start = basis_update.work_estimate(); basis_update.b_transpose_solve(e_i, u_bar); + const f_t actual_basis_work = basis_update.work_estimate() - basis_work_start; + if (actual_basis_work > predicted_basis_work && + add_work_estimate( + actual_basis_work - predicted_basis_work, &work_estimate, max_work_estimate)) { + return -2; + } #ifdef CHECK_B_TRANSPOSE_SOLVE std::vector u_bar_dense(lp.num_rows); @@ -4832,6 +5010,11 @@ i_t tableau_equality_t::generate_base_equality( // Compute a_bar = N^T u_bar // TODO: This is similar to a function in phase2 of dual simplex. See if it can be reused. const i_t nz_ubar = u_bar.i.size(); + f_t tableau_multiply_work = static_cast(3 * nz_ubar + 2); + for (const i_t row : u_bar.i) { + tableau_multiply_work += static_cast(6 * (Arow.row_start[row + 1] - Arow.row_start[row])); + } + if (add_work_estimate(tableau_multiply_work, &work_estimate, max_work_estimate)) { return -2; } std::vector abar_indices; abar_indices.reserve(nz_ubar); for (i_t k = 0; k < nz_ubar; k++) { @@ -4892,6 +5075,10 @@ i_t tableau_equality_t::generate_base_equality( // Check that the tableau equality is satisfied const f_t tableau_tol = 1e-6; + if (add_work_estimate( + static_cast(4 * a_bar.i.size() + 2), &work_estimate, max_work_estimate)) { + return -2; + } f_t a_bar_dot_xstar = a_bar.dot(xstar); if (std::abs(a_bar_dot_xstar - b_bar_[i]) > tableau_tol) { settings.log.debug("bad tableau equality. error %e\n", std::abs(a_bar_dot_xstar - b_bar_[i])); @@ -6024,10 +6211,17 @@ bool complemented_mixed_integer_rounding_cut_t::generate_lifted_mixed_ const std::vector& var_types, const std::vector& transformed_xstar, inequality_t& transformed_cut, - f_t& work_estimate) + f_t& work_estimate, + f_t max_work_estimate) { constexpr f_t tolerance = static_cast(1e-6); + const f_t estimated_work = + static_cast(12 * transformed_inequality.size()) + + static_cast(transformed_inequality.size()) * + std::log2(static_cast(transformed_inequality.size()) + static_cast(1.0)); + if (add_work_estimate(estimated_work, &work_estimate, max_work_estimate)) { return false; } + // Work in <= form. All variables have already been shifted or complemented // to nonnegative variables by bound_substitution()/transform_inequality(). inequality_t base = transformed_inequality; @@ -6152,9 +6346,6 @@ bool complemented_mixed_integer_rounding_cut_t::generate_lifted_mixed_ transformed_cut = std::move(squeezed_cut); transformed_cut.negate(); - work_estimate += static_cast(12 * base.size()) + - static_cast(cover.size()) * - std::log2(static_cast(cover.size()) + static_cast(1.0)); return true; } diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index 0f13261a61..e023618664 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -17,6 +17,7 @@ #include #include #include +#include #include #include #include @@ -54,6 +55,29 @@ struct cut_gap_closure_t { f_t gap_closed_ratio{0.0}; }; +// Deterministic operation estimates for one separator pass. Raw estimates are +// kept per phase so expensive cut families remain attributable; callers can +// convert the total to solver work units with total() / 1e8. +template +struct cut_work_stats_t { + f_t gomory{0.0}; + f_t knapsack{0.0}; + f_t flow_cover{0.0}; + f_t mir{0.0}; + f_t implied_bound{0.0}; + f_t conflict_graph{0.0}; + f_t clique{0.0}; + f_t zero_half{0.0}; + + f_t total() const + { + return gomory + knapsack + flow_cover + mir + implied_bound + conflict_graph + clique + + zero_half; + } + + f_t work_units() const { return total() / static_cast(1e8); } +}; + template cut_gap_closure_t compute_cut_gap_closure(f_t objective_reference, f_t objective_before_cuts, @@ -296,6 +320,14 @@ std::vector> find_mod2_row_combinations_for_test( 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: @@ -314,7 +346,10 @@ class cut_pool_t { // We expect that the cut is violated by the current relaxation xstar. void add_cut(cut_type_t cut_type, const inequality_t& cut); - void score_cuts(std::vector& x_relax); + // Returns a deterministic operation estimate. Orthogonality selection is + // stopped at max_work_estimate, preserving the best cuts found so far. + f_t score_cuts(std::vector& x_relax, + f_t max_work_estimate = std::numeric_limits::infinity()); // We return the cuts in the form best_cuts*x <= best_rhs i_t get_best_cuts(csr_matrix_t& best_cuts, @@ -330,6 +365,7 @@ class cut_pool_t { void print_cutpool_types() { print_cut_types("In cut pool", cut_type_, settings_); } void check_for_duplicate_cuts(); + bool check_for_duplicate_cuts(f_t& work_estimate, f_t max_work_estimate); private: f_t cut_distance(i_t row, const std::vector& x, f_t& cut_violation, f_t& cut_norm); @@ -669,6 +705,8 @@ class cut_generation_t { variable_bounds_t& variable_bounds, f_t start_time); + const cut_work_stats_t& last_work_stats() const { return last_work_stats_; } + private: // Generate all mixed integer gomory cuts void generate_gomory_cuts(const simplex::lp_problem_t& lp, @@ -751,6 +789,7 @@ class cut_generation_t { std::shared_ptr> clique_table_; omp_atomic_t* signal_extend_{nullptr}; fractional_conflict_subgraph_t sub_cg_; + cut_work_stats_t last_work_stats_; }; template @@ -839,7 +878,9 @@ class tableau_equality_t { const std::vector& basic_list, const std::vector& nonbasic_list, i_t i, - inequality_t& inequality); + inequality_t& inequality, + f_t& work_estimate, + f_t max_work_estimate); private: std::vector b_bar_; @@ -1019,7 +1060,8 @@ class complemented_mixed_integer_rounding_cut_t { const std::vector& var_types, const std::vector& transformed_xstar, inequality_t& transformed_cut, - f_t& work_estimate); + f_t& work_estimate, + f_t max_work_estimate); f_t compute_violation(const inequality_t& cut, const std::vector& xstar); diff --git a/cpp/src/dual_simplex/phase2.cpp b/cpp/src/dual_simplex/phase2.cpp index a5f10c3229..8c11bec81b 100644 --- a/cpp/src/dual_simplex/phase2.cpp +++ b/cpp/src/dual_simplex/phase2.cpp @@ -2559,6 +2559,23 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, f_t phase2_work_estimate = 0.0; ft.clear_work_estimate(); + struct final_work_flush_t { + work_limit_context_t* context; + f_t& phase_work; + basis_update_mpf_t& basis_update; + + ~final_work_flush_t() + { + if (context == nullptr) { return; } + phase_work += basis_update.work_estimate(); + basis_update.clear_work_estimate(); + if (phase_work > static_cast(0.0)) { + context->record_work_sync_on_horizon(phase_work / static_cast(1e8)); + phase_work = 0.0; + } + } + } final_work_flush{work_unit_context, phase2_work_estimate, ft}; + std::vector& x = sol.x; std::vector& y = sol.y; std::vector& z = sol.z; diff --git a/cpp/src/utilities/work_limit_context.hpp b/cpp/src/utilities/work_limit_context.hpp index 3463847123..f7e2b537f1 100644 --- a/cpp/src/utilities/work_limit_context.hpp +++ b/cpp/src/utilities/work_limit_context.hpp @@ -18,6 +18,7 @@ struct work_limit_context_t { double global_work_units_elapsed{0.0}; double total_sync_time{0.0}; // Total time spent waiting at sync barriers (seconds) bool deterministic{false}; + bool sync_on_horizon{true}; work_unit_scheduler_t* scheduler{nullptr}; std::string name; @@ -27,7 +28,9 @@ struct work_limit_context_t { { if (!deterministic) return; global_work_units_elapsed += work; - if (scheduler) { scheduler->on_work_recorded(*this, global_work_units_elapsed); } + if (scheduler && sync_on_horizon) { + scheduler->on_work_recorded(*this, global_work_units_elapsed); + } } }; diff --git a/cpp/tests/mip/cuts_test.cu b/cpp/tests/mip/cuts_test.cu index 08b978216d..08987dd687 100644 --- a/cpp/tests/mip/cuts_test.cu +++ b/cpp/tests/mip/cuts_test.cu @@ -981,6 +981,22 @@ TEST(cuts, test_duplicate_cuts_detection) cut_pool.add_cut(mip::cut_type_t::MIXED_INTEGER_GOMORY, cut8); cut_pool.check_for_duplicate_cuts(); + + std::vector xstar(4, 0.0); + const double scoring_work = cut_pool.score_cuts(xstar, 1.0); + EXPECT_GT(scoring_work, 1.0); + EXPECT_TRUE(std::isfinite(scoring_work)); +} + +TEST(cuts, cut_work_stats_total) +{ + mip::cut_work_stats_t stats; + stats.gomory = 1.0; + stats.mir = 2.0; + stats.zero_half = 3.0; + stats.conflict_graph = 4.0; + EXPECT_DOUBLE_EQ(stats.total(), 10.0); + EXPECT_DOUBLE_EQ(stats.work_units(), 1e-7); } TEST(cuts, clique_phase1_smoke_conflict_graph_edges) @@ -1516,6 +1532,26 @@ TEST(cuts, zero_half_unit_mod2_row_finder_single_pair_and_four_row_dependencies) } } +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); + + // The budget admits preprocessing and the first basis row, then expires at + // the next symmetric difference. The search must return immediately rather + // than repeating an over-budget reduction for every remaining input row. + 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 From 8b6907a6b1fd2c3ea0e530dbab686c6fa4235115 Mon Sep 17 00:00:00 2001 From: akif Date: Wed, 5 Aug 2026 14:08:30 +0200 Subject: [PATCH 03/15] Coarsen zero-half work accounting Track work at meaningful phase boundaries and document the convention so cut generation avoids inaccurate or redundant limit checks. Signed-off-by: Akif Corduk --- cpp/src/cuts/cuts.cpp | 382 +++++++++++------- skills/cuopt-developer/SKILL.md | 2 +- .../cuopt-developer/references/conventions.md | 17 +- 3 files changed, 250 insertions(+), 151 deletions(-) diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index 49d99389e1..6b3012f753 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -1155,11 +1155,10 @@ std::vector> find_mod2_row_combinations( bool rhs{false}; }; - i_t max_index = -1; + i_t max_index = -1; + f_t input_scan_work = 0.0; for (const auto& row : parity_rows) { - if (add_work_estimate(static_cast(row.size() + 1), work_estimate, max_work_estimate)) { - return {}; - } + input_scan_work += static_cast(row.size() + 1); cuopt_assert(std::is_sorted(row.begin(), row.end()), "GF(2) parity rows must be sorted"); cuopt_assert(std::adjacent_find(row.begin(), row.end()) == row.end(), "GF(2) parity rows must not contain duplicates"); @@ -1168,6 +1167,7 @@ std::vector> find_mod2_row_combinations( max_index = std::max(max_index, row.back()); } } + if (add_work_estimate(input_scan_work, work_estimate, max_work_estimate)) { return {}; } std::vector permutation(parity_rows.size()); std::iota(permutation.begin(), permutation.end(), 0); @@ -1196,10 +1196,7 @@ std::vector> find_mod2_row_combinations( std::vector parity_tmp; std::vector combination_tmp; for (const i_t candidate : permutation) { - if (add_work_estimate( - static_cast(parity_rows[candidate].size() + 2), work_estimate, max_work_estimate)) { - break; - } + f_t candidate_work = static_cast(parity_rows[candidate].size() + 2); basis_row_t current; current.parity = parity_rows[candidate]; current.combination = {candidate}; @@ -1212,13 +1209,8 @@ std::vector> find_mod2_row_combinations( if (basis_index < 0) { break; } const auto& pivot_row = basis[basis_index]; - if (add_work_estimate( - static_cast(current.parity.size() + pivot_row.parity.size() + - current.combination.size() + pivot_row.combination.size()), - work_estimate, - max_work_estimate)) { - return combinations; - } + candidate_work += static_cast(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() > static_cast(max_combination_size)) { @@ -1229,6 +1221,7 @@ std::vector> find_mod2_row_combinations( 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; } if (current.parity.empty()) { @@ -1242,11 +1235,148 @@ std::vector> find_mod2_row_combinations( const i_t pivot = current.parity.front(); pivot_to_basis[pivot] = static_cast(basis.size()); basis.push_back(std::move(current)); - if (work_estimate != nullptr && *work_estimate > max_work_estimate) { break; } } return combinations; } +template +struct mod2_candidate_t { + inequality_t transformed_inequality; + std::vector parity; + bool rhs_parity{false}; + bool reversible{false}; +}; + +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) +{ + if (toc(start_time) >= time_limit) { return i_t{0}; } + f_t scale_work = 0.0; + for (i_t scale = 1; scale <= max_integral_scale; ++scale) { + 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 (integral) { + work_estimate += scale_work; + return scale; + } + } + work_estimate += scale_work; + return i_t{0}; +} + +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)(10 * 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); + complemented_mir.remove_small_coefficients(lp.lower, lp.upper, transformed_cut); + if (complemented_mir.compute_violation(transformed_cut, xstar) > min_violation) { + cut_pool.add_cut(cut_type_t::ZERO_HALF, transformed_cut); + ++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) +{ + work_estimate += (f_t)(3 * oriented_aggregate.size() + 1); + inequality_t mir_cut(lp.num_cols); + if (complemented_mir.generate_cut_nonnegative_maintain_indicies( + oriented_aggregate, var_types, mir_cut)) { + 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); + if (toc(start_time) < settings.time_limit && + complemented_mir.generate_lifted_mixed_binary_cover(oriented_aggregate, + var_types, + transformed_xstar, + lifted_cover_cut, + work_estimate, + max_work_estimate)) { + 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; } +} + // 64-bit integer mixer (SplitMix64). Used as the building block for the // cousin filter's per-slot independent hash family. inline uint64_t splitmix64_mix(uint64_t x) @@ -1590,7 +1720,7 @@ f_t cut_pool_t::score_cuts(std::vector& x_relax, f_t max_work_est f_t cut_ortho = 1.0; const i_t best_cuts_size = best_cuts_.size(); for (i_t k = 0; k < best_cuts_size; k++) { - const i_t j = best_cuts_[k]; + const i_t j = best_cuts_[k]; const i_t i_nz = cut_storage_.row_start[i + 1] - cut_storage_.row_start[i]; const i_t j_nz = cut_storage_.row_start[j + 1] - cut_storage_.row_start[j]; if (add_work_estimate( @@ -1598,7 +1728,7 @@ f_t cut_pool_t::score_cuts(std::vector& x_relax, f_t max_work_est work_limit_reached = true; break; } - cut_ortho = std::min(cut_ortho, cut_orthogonality(i, j)); + cut_ortho = std::min(cut_ortho, cut_orthogonality(i, j)); } if (!work_limit_reached && cut_ortho >= min_orthogonality) { best_cuts_.push_back(i); @@ -4106,26 +4236,22 @@ bool cut_generation_t::generate_zero_half_cuts( // Dividing that aggregation by two and applying c-MIR yields a valid // zero-half cut. The sparse elimination below bounds both combination size // and work so this path remains predictable on large models. - constexpr i_t max_integral_scale = 1000; - constexpr i_t max_combination_size = 64; - constexpr i_t max_row_combinations = 1000; - const i_t max_integer_row_length = 1000 + lp.num_cols / 10; - const f_t row_tight_tol = std::max(settings.primal_tol, static_cast(1e-7)); - const f_t coefficient_integral_tol = static_cast(1e-6); - const f_t min_violation = std::max(settings.primal_tol, static_cast(1e-6)); - f_t& mod2_work_estimate = last_work_stats_.zero_half; - const f_t max_mod2_work_estimate = mod2_work_estimate + static_cast(1e8); - bool mod2_work_limit_reached = false; - - auto charge_mod2_work = [&](f_t work) { - const bool limit = add_work_estimate(work, &mod2_work_estimate, max_mod2_work_estimate); - mod2_work_limit_reached = mod2_work_limit_reached || limit; - return limit; - }; - - if (charge_mod2_work(static_cast(3 * lp.num_cols) + - static_cast(variable_bounds.upper_variables.size() + - variable_bounds.lower_variables.size()))) { + constexpr i_t max_integral_scale = 1000; + constexpr i_t max_combination_size = 64; + constexpr i_t max_row_combinations = 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; + constexpr f_t min_violation = (f_t)1e-6; + f_t& mod2_work_estimate = last_work_stats_.zero_half; + const f_t max_mod2_work_estimate = mod2_work_estimate + (f_t)1e8; + bool mod2_work_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()), + &mod2_work_estimate, + max_mod2_work_estimate, + &mod2_work_limit_reached)) { return true; } complemented_mixed_integer_rounding_cut_t complemented_mir(lp, settings, new_slacks); @@ -4133,59 +4259,22 @@ bool cut_generation_t::generate_zero_half_cuts( complemented_mir.bound_substitution( lp, variable_bounds, var_types, xstar, transformed_xstar, true); - struct mod2_candidate_t { - inequality_t transformed_inequality; - std::vector parity; - bool rhs_parity{false}; - bool reversible{false}; - }; - - std::vector mod2_candidates; + std::vector> mod2_candidates; mod2_candidates.reserve(lp.num_rows); - auto integral_scale = [&](const inequality_t& inequality) { - for (i_t scale = 1; scale <= max_integral_scale; ++scale) { - if (toc(start_time) >= settings.time_limit || - charge_mod2_work(static_cast(inequality.size() + 1))) { - return i_t{0}; - } - bool integral = true; - const f_t scaled_rhs = static_cast(scale) * inequality.rhs; - if (std::abs(scaled_rhs - std::round(scaled_rhs)) > - coefficient_integral_tol * std::max(static_cast(1.0), std::abs(scaled_rhs))) { - integral = false; - } - for (i_t k = 0; integral && k < static_cast(inequality.size()); ++k) { - 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 = static_cast(scale) * inequality.coeff(k); - if (std::abs(scaled_coefficient - std::round(scaled_coefficient)) > - coefficient_integral_tol * - std::max(static_cast(1.0), std::abs(scaled_coefficient))) { - integral = false; - } - } - if (integral) { return scale; } - } - return i_t{0}; - }; - for (i_t row = 0; row < lp.num_rows; ++row) { - if (toc(start_time) >= settings.time_limit || mod2_work_limit_reached || - mod2_work_estimate > max_mod2_work_estimate) { - break; - } + if (toc(start_time) >= settings.time_limit || mod2_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 = static_cast(8 * row_length + 5) + - static_cast(row_length) * - std::log2(static_cast(row_length) + static_cast(1.0)); - if (charge_mod2_work(row_work)) { break; } + 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, &mod2_work_estimate, max_mod2_work_estimate, &mod2_work_limit_reached)) { + break; + } inequality_t inequality(Arow, row, lp.rhs[row]); complemented_mir.transform_inequality(variable_bounds, var_types, inequality); inequality.sort(); @@ -4221,11 +4310,23 @@ bool cut_generation_t::generate_zero_half_cuts( } if (!continuous_at_bounds) { continue; } - const i_t scale = integral_scale(inequality); + const i_t scale = mod2_integral_scale(inequality, + var_types, + transformed_xstar, + max_integral_scale, + row_tight_tol, + coefficient_integral_tol, + start_time, + settings.time_limit, + mod2_work_estimate); + if (mod2_work_estimate > max_mod2_work_estimate) { + mod2_work_limit_reached = true; + break; + } if (scale == 0) { continue; } - if (scale != 1) { inequality.scale(static_cast(scale)); } + if (scale != 1) { inequality.scale((f_t)scale); } - mod2_candidate_t candidate; + mod2_candidate_t candidate; candidate.transformed_inequality = std::move(inequality); candidate.rhs_parity = (std::llabs(static_cast(std::llround(candidate.transformed_inequality.rhs))) % @@ -4249,36 +4350,34 @@ bool cut_generation_t::generate_zero_half_cuts( parity_rows.reserve(mod2_candidates.size()); rhs_parity.reserve(mod2_candidates.size()); for (const auto& candidate : mod2_candidates) { - if (charge_mod2_work(static_cast(candidate.parity.size() + 1))) { break; } parity_rows.push_back(candidate.parity); rhs_parity.push_back(candidate.rhs_parity); } - auto row_combinations = find_mod2_row_combinations(parity_rows, + auto row_combinations = find_mod2_row_combinations(parity_rows, rhs_parity, max_combination_size, max_row_combinations, &mod2_work_estimate, max_mod2_work_estimate); - mod2_work_limit_reached = mod2_work_limit_reached || mod2_work_estimate > max_mod2_work_estimate; + if (mod2_work_estimate > max_mod2_work_estimate) { mod2_work_limit_reached = true; } scratch_pad_t aggregate_pad(lp.num_cols); i_t mod2_cuts_added = 0; for (const auto& combination : row_combinations) { - if (toc(start_time) >= settings.time_limit || mod2_work_limit_reached || - mod2_work_estimate > max_mod2_work_estimate) { - break; - } + if (toc(start_time) >= settings.time_limit || mod2_work_limit_reached) { break; } size_t aggregate_input_nz = 0; for (const i_t candidate_index : combination) { aggregate_input_nz += mod2_candidates[candidate_index].transformed_inequality.size(); } const f_t aggregate_work = - static_cast(4 * aggregate_input_nz + 1) + - static_cast(aggregate_input_nz) * - std::log2(static_cast(aggregate_input_nz) + static_cast(1.0)); - if (charge_mod2_work(aggregate_work)) { break; } + (f_t)(4 * aggregate_input_nz + 1) + + (f_t)aggregate_input_nz * std::log2((f_t)aggregate_input_nz + (f_t)1.0); + if (add_work_estimate( + aggregate_work, &mod2_work_estimate, max_mod2_work_estimate, &mod2_work_limit_reached)) { + break; + } inequality_t aggregate(lp.num_cols); bool reversible = true; @@ -4286,7 +4385,7 @@ bool cut_generation_t::generate_zero_half_cuts( const auto& candidate = mod2_candidates[candidate_index]; aggregate.rhs += candidate.transformed_inequality.rhs; reversible = reversible && candidate.reversible; - for (i_t k = 0; k < static_cast(candidate.transformed_inequality.size()); ++k) { + 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)); } @@ -4294,53 +4393,42 @@ bool cut_generation_t::generate_zero_half_cuts( aggregate_pad.get_pad(aggregate.vector.i, aggregate.vector.x); aggregate_pad.clear_pad(); aggregate.sort(); - aggregate.scale(static_cast(0.5)); - - auto generate_from_aggregate = [&](const inequality_t& oriented_aggregate) { - auto add_transformed_cut = [&](inequality_t transformed_cut) { - if (toc(start_time) >= settings.time_limit || - charge_mod2_work(static_cast(10 * transformed_cut.size() + 1))) { - return; - } - 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); - 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) { - cut_pool_.add_cut(cut_type_t::ZERO_HALF, transformed_cut); - mod2_cuts_added++; - } - }; - - if (toc(start_time) >= settings.time_limit || - charge_mod2_work(static_cast(3 * oriented_aggregate.size() + 1))) { - return; - } - inequality_t mir_cut(lp.num_cols); - if (complemented_mir.generate_cut_nonnegative_maintain_indicies( - oriented_aggregate, var_types, mir_cut)) { - add_transformed_cut(std::move(mir_cut)); - } - - inequality_t lifted_cover_cut(lp.num_cols); - if (toc(start_time) < settings.time_limit && - complemented_mir.generate_lifted_mixed_binary_cover(oriented_aggregate, - var_types, - transformed_xstar, - lifted_cover_cut, - mod2_work_estimate, - max_mod2_work_estimate)) { - add_transformed_cut(std::move(lifted_cover_cut)); - } - mod2_work_limit_reached = - mod2_work_limit_reached || mod2_work_estimate > max_mod2_work_estimate; - }; - - generate_from_aggregate(aggregate); - if (reversible) { + aggregate.scale((f_t)0.5); + + mod2_generate_cuts_from_aggregate(complemented_mir, + cut_pool_, + lp, + settings, + Arow, + variable_bounds, + var_types, + xstar, + transformed_xstar, + aggregate, + min_violation, + start_time, + mod2_work_estimate, + max_mod2_work_estimate, + mod2_work_limit_reached, + mod2_cuts_added); + if (reversible && toc(start_time) < settings.time_limit && !mod2_work_limit_reached) { aggregate.negate(); - generate_from_aggregate(aggregate); + mod2_generate_cuts_from_aggregate(complemented_mir, + cut_pool_, + lp, + settings, + Arow, + variable_bounds, + var_types, + xstar, + transformed_xstar, + aggregate, + min_violation, + start_time, + mod2_work_estimate, + max_mod2_work_estimate, + mod2_work_limit_reached, + mod2_cuts_added); } } ZERO_HALF_DEBUG("general mod2 candidates=%lld combinations=%lld cuts=%lld work=%g", @@ -4546,7 +4634,7 @@ void cut_generation_t::generate_mir_cuts( work_estimate += static_cast(4 * lp.num_cols + 3 * lp.num_rows) + static_cast(4 * Arow.row_start[Arow.m]); const f_t max_work_estimate = work_estimate + static_cast(2e9); - i_t num_cuts = 0; + i_t num_cuts = 0; while (num_cuts < max_cuts && !score_queue.empty()) { if (toc(start_time) >= settings.time_limit) { break; } // Get the row with the highest score from the queue @@ -5009,7 +5097,7 @@ i_t tableau_equality_t::generate_base_equality( // Compute a_bar = N^T u_bar // TODO: This is similar to a function in phase2 of dual simplex. See if it can be reused. - const i_t nz_ubar = u_bar.i.size(); + const i_t nz_ubar = u_bar.i.size(); f_t tableau_multiply_work = static_cast(3 * nz_ubar + 2); for (const i_t row : u_bar.i) { tableau_multiply_work += static_cast(6 * (Arow.row_start[row + 1] - Arow.row_start[row])); @@ -5079,7 +5167,7 @@ i_t tableau_equality_t::generate_base_equality( static_cast(4 * a_bar.i.size() + 2), &work_estimate, max_work_estimate)) { return -2; } - f_t a_bar_dot_xstar = a_bar.dot(xstar); + f_t a_bar_dot_xstar = a_bar.dot(xstar); if (std::abs(a_bar_dot_xstar - b_bar_[i]) > tableau_tol) { settings.log.debug("bad tableau equality. error %e\n", std::abs(a_bar_dot_xstar - b_bar_[i])); return -1; diff --git a/skills/cuopt-developer/SKILL.md b/skills/cuopt-developer/SKILL.md index 5188663ed9..eb8de49a63 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.) From f1b72ce35acf33678af1749f2cf485d86ed22d5b Mon Sep 17 00:00:00 2001 From: akif Date: Wed, 5 Aug 2026 14:59:24 +0200 Subject: [PATCH 04/15] Isolate mod-2 zero-half cut generation Move the separator into a dedicated compilation unit and retain candidate parity data without duplicate storage. Signed-off-by: akif --- cpp/src/cuts/CMakeLists.txt | 1 + cpp/src/cuts/cuts.cpp | 649 +---------------------------- cpp/src/cuts/cuts.hpp | 12 + cpp/src/cuts/zero_half_mod2.cpp | 715 ++++++++++++++++++++++++++++++++ 4 files changed, 740 insertions(+), 637 deletions(-) create mode 100644 cpp/src/cuts/zero_half_mod2.cpp 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 6b3012f753..982764c4a8 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -1126,257 +1126,6 @@ std::vector> find_violated_odd_cycles_for_test( 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 -std::vector> find_mod2_row_combinations( - const std::vector>& parity_rows, - const std::vector& rhs_parity, - i_t max_combination_size, - i_t max_combinations, - f_t* work_estimate, - f_t max_work_estimate) -{ - cuopt_assert(parity_rows.size() == rhs_parity.size(), - "GF(2) parity row and rhs sizes must match"); - if (max_combination_size <= 0 || max_combinations <= 0) { return {}; } - - struct basis_row_t { - std::vector parity; - std::vector combination; - bool rhs{false}; - }; - - i_t max_index = -1; - f_t input_scan_work = 0.0; - for (const auto& row : parity_rows) { - input_scan_work += static_cast(row.size() + 1); - cuopt_assert(std::is_sorted(row.begin(), row.end()), "GF(2) parity rows must be sorted"); - cuopt_assert(std::adjacent_find(row.begin(), row.end()) == row.end(), - "GF(2) parity rows must not contain duplicates"); - if (!row.empty()) { - cuopt_assert(row.front() >= 0, "GF(2) parity index must be nonnegative"); - max_index = std::max(max_index, row.back()); - } - } - if (add_work_estimate(input_scan_work, work_estimate, max_work_estimate)) { return {}; } - - std::vector permutation(parity_rows.size()); - std::iota(permutation.begin(), permutation.end(), 0); - const f_t sort_work = static_cast(permutation.size()) * - std::log2(static_cast(permutation.size()) + static_cast(1.0)); - if (add_work_estimate(sort_work, work_estimate, max_work_estimate)) { return {}; } - std::stable_sort(permutation.begin(), permutation.end(), [&](i_t a, i_t b) { - if (parity_rows[a].size() != parity_rows[b].size()) { - return parity_rows[a].size() < parity_rows[b].size(); - } - // Prefer an even-rhs representative for a pivot. If an otherwise - // identical odd-rhs row arrives later, their dependency is immediately a - // valid zero-half aggregation. - return rhs_parity[a] < rhs_parity[b]; - }); - - if (add_work_estimate(static_cast(max_index + 1), work_estimate, max_work_estimate)) { - return {}; - } - std::vector pivot_to_basis(static_cast(max_index + 1), -1); - std::vector basis; - basis.reserve(std::min(parity_rows.size(), static_cast(max_index + 1))); - std::vector> combinations; - combinations.reserve(std::min(static_cast(max_combinations), parity_rows.size())); - - std::vector parity_tmp; - std::vector combination_tmp; - for (const i_t candidate : permutation) { - f_t candidate_work = static_cast(parity_rows[candidate].size() + 2); - basis_row_t current; - current.parity = parity_rows[candidate]; - current.combination = {candidate}; - current.rhs = rhs_parity[candidate] != 0; - - bool abandoned = false; - while (!current.parity.empty()) { - const i_t pivot = current.parity.front(); - const i_t basis_index = pivot_to_basis[pivot]; - if (basis_index < 0) { break; } - - const auto& pivot_row = basis[basis_index]; - candidate_work += static_cast(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() > static_cast(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; } - - if (current.parity.empty()) { - if (current.rhs && !current.combination.empty()) { - combinations.push_back(std::move(current.combination)); - if (combinations.size() >= static_cast(max_combinations)) { break; } - } - continue; - } - - const i_t pivot = current.parity.front(); - pivot_to_basis[pivot] = static_cast(basis.size()); - basis.push_back(std::move(current)); - } - return combinations; -} - -template -struct mod2_candidate_t { - inequality_t transformed_inequality; - std::vector parity; - bool rhs_parity{false}; - bool reversible{false}; -}; - -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) -{ - if (toc(start_time) >= time_limit) { return i_t{0}; } - f_t scale_work = 0.0; - for (i_t scale = 1; scale <= max_integral_scale; ++scale) { - 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 (integral) { - work_estimate += scale_work; - return scale; - } - } - work_estimate += scale_work; - return i_t{0}; -} - -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)(10 * 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); - complemented_mir.remove_small_coefficients(lp.lower, lp.upper, transformed_cut); - if (complemented_mir.compute_violation(transformed_cut, xstar) > min_violation) { - cut_pool.add_cut(cut_type_t::ZERO_HALF, transformed_cut); - ++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) -{ - work_estimate += (f_t)(3 * oriented_aggregate.size() + 1); - inequality_t mir_cut(lp.num_cols); - if (complemented_mir.generate_cut_nonnegative_maintain_indicies( - oriented_aggregate, var_types, mir_cut)) { - 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); - if (toc(start_time) < settings.time_limit && - complemented_mir.generate_lifted_mixed_binary_cover(oriented_aggregate, - var_types, - transformed_xstar, - lifted_cover_cut, - work_estimate, - max_work_estimate)) { - 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; } -} - // 64-bit integer mixer (SplitMix64). Used as the building block for the // cousin filter's per-slot independent hash family. inline uint64_t splitmix64_mix(uint64_t x) @@ -1394,39 +1143,6 @@ inline uint64_t hash64_with_seed(uint64_t value, uint64_t seed) } // 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) -{ - double work_estimate = 0.0; - auto combinations = find_mod2_row_combinations(parity_rows, - rhs_parity, - max_combination_size, - max_combinations, - &work_estimate, - max_work_estimate); - if (work_estimate_out != nullptr) { *work_estimate_out = work_estimate; } - return combinations; -} - template void cut_pool_t::add_cut(cut_type_t cut_type, const inequality_t& cut) { @@ -4226,216 +3942,18 @@ bool cut_generation_t::generate_zero_half_cuts( static_cast(sub_cg_.ready), sub_cg_.vertices.size()); - // General zero-half separation. The previous implementation only searched - // for odd cycles in the conflict graph. A zero-half separator must also find - // small GF(2) dependencies among tight rows: - // - // sum_i (a_i x >= b_i), with even aggregate integer coefficients and an - // odd aggregate rhs. - // - // Dividing that aggregation by two and applying c-MIR yields a valid - // zero-half cut. The sparse elimination below bounds both combination size - // and work so this path remains predictable on large models. - constexpr i_t max_integral_scale = 1000; - constexpr i_t max_combination_size = 64; - constexpr i_t max_row_combinations = 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; - constexpr f_t min_violation = (f_t)1e-6; - f_t& mod2_work_estimate = last_work_stats_.zero_half; - const f_t max_mod2_work_estimate = mod2_work_estimate + (f_t)1e8; - bool mod2_work_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()), - &mod2_work_estimate, - max_mod2_work_estimate, - &mod2_work_limit_reached)) { + if (!generate_mod2_zero_half_cuts(cut_pool_, + lp, + settings, + Arow, + new_slacks, + var_types, + xstar, + variable_bounds, + start_time, + last_work_stats_.zero_half)) { return true; } - 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); - - std::vector> mod2_candidates; - mod2_candidates.reserve(lp.num_rows); - - for (i_t row = 0; row < lp.num_rows; ++row) { - if (toc(start_time) >= settings.time_limit || mod2_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, &mod2_work_estimate, max_mod2_work_estimate, &mod2_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. Choose the direction - // in which the zero-valued transformed slack has a negative coefficient; - // it can then be removed while preserving a valid >= inequality. - i_t slack_position = -1; - for (i_t k = 0; k < static_cast(inequality.size()); ++k) { - if (inequality.index(k) == slack) { - slack_position = k; - break; - } - } - if (slack_position < 0 || inequality.coeff(slack_position) == 0.0) { continue; } - 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); - - // As in general mod-2 separators, only rows whose continuous variables - // are at their selected bounds participate in the parity system. - bool continuous_at_bounds = true; - for (i_t k = 0; k < static_cast(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, - settings.time_limit, - mod2_work_estimate); - if (mod2_work_estimate > max_mod2_work_estimate) { - mod2_work_limit_reached = true; - break; - } - 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::llabs(static_cast(std::llround(candidate.transformed_inequality.rhs))) % - 2) != 0; - candidate.reversible = std::abs(lp.upper[slack] - lp.lower[slack]) <= row_tight_tol; - for (i_t k = 0; k < static_cast(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 = - static_cast(std::llround(candidate.transformed_inequality.coeff(k))); - if ((std::llabs(coefficient) % 2) != 0) { candidate.parity.push_back(j); } - } - if (candidate.parity.size() > static_cast(max_integer_row_length)) { continue; } - mod2_candidates.push_back(std::move(candidate)); - } - - std::vector> parity_rows; - std::vector rhs_parity; - parity_rows.reserve(mod2_candidates.size()); - rhs_parity.reserve(mod2_candidates.size()); - for (const auto& candidate : mod2_candidates) { - parity_rows.push_back(candidate.parity); - rhs_parity.push_back(candidate.rhs_parity); - } - - auto row_combinations = find_mod2_row_combinations(parity_rows, - rhs_parity, - max_combination_size, - max_row_combinations, - &mod2_work_estimate, - max_mod2_work_estimate); - if (mod2_work_estimate > max_mod2_work_estimate) { mod2_work_limit_reached = true; } - scratch_pad_t aggregate_pad(lp.num_cols); - i_t mod2_cuts_added = 0; - - for (const auto& combination : row_combinations) { - if (toc(start_time) >= settings.time_limit || mod2_work_limit_reached) { break; } - - size_t aggregate_input_nz = 0; - for (const i_t candidate_index : combination) { - aggregate_input_nz += mod2_candidates[candidate_index].transformed_inequality.size(); - } - const f_t aggregate_work = - (f_t)(4 * aggregate_input_nz + 1) + - (f_t)aggregate_input_nz * std::log2((f_t)aggregate_input_nz + (f_t)1.0); - if (add_work_estimate( - aggregate_work, &mod2_work_estimate, max_mod2_work_estimate, &mod2_work_limit_reached)) { - break; - } - - inequality_t aggregate(lp.num_cols); - bool reversible = true; - for (const i_t candidate_index : combination) { - const auto& candidate = mod2_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(); - aggregate.sort(); - aggregate.scale((f_t)0.5); - - mod2_generate_cuts_from_aggregate(complemented_mir, - cut_pool_, - lp, - settings, - Arow, - variable_bounds, - var_types, - xstar, - transformed_xstar, - aggregate, - min_violation, - start_time, - mod2_work_estimate, - max_mod2_work_estimate, - mod2_work_limit_reached, - mod2_cuts_added); - if (reversible && toc(start_time) < settings.time_limit && !mod2_work_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, - mod2_work_estimate, - max_mod2_work_estimate, - mod2_work_limit_reached, - mod2_cuts_added); - } - } - ZERO_HALF_DEBUG("general mod2 candidates=%lld combinations=%lld cuts=%lld work=%g", - static_cast(mod2_candidates.size()), - static_cast(row_combinations.size()), - static_cast(mod2_cuts_added), - static_cast(mod2_work_estimate)); // The fractional conflict-graph subgraph is built once per cut pass in // prepare_fractional_sub_conflict_graph() and remains a complementary @@ -4460,7 +3978,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 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 = last_work_stats_.zero_half; @@ -6293,150 +5812,6 @@ bool complemented_mixed_integer_rounding_cut_t:: 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 = static_cast(1e-6); - - const f_t estimated_work = - static_cast(12 * transformed_inequality.size()) + - static_cast(transformed_inequality.size()) * - std::log2(static_cast(transformed_inequality.size()) + static_cast(1.0)); - if (add_work_estimate(estimated_work, &work_estimate, max_work_estimate)) { return false; } - - // Work in <= form. All variables have already been shifted or complemented - // to nonnegative variables by bound_substitution()/transform_inequality(). - 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 < static_cast(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]; - // Positive continuous coefficients can be relaxed from a <= row. - if (aj > 0.0) { base.vector.x[k] = 0.0; } - continue; - } - - const f_t upper = new_upper(j); - // This lifting function is for binary variables. General integer rows - // continue to use the c-MIR zero-half cut. - if (upper == inf || std::abs(upper - static_cast(1.0)) > tolerance) { return false; } - is_integral[k] = 1; - if (aj < 0.0) { - // z = 1 - w makes the knapsack coefficient positive. - 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 < static_cast(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(), [&](i_t a, i_t b) { - 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); - }); - - 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(static_cast(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; } - - auto lifting_function = [&](f_t coefficient) { - for (size_t h = 0; h < p; ++h) { - if (coefficient <= prefix[h] - lambda + tolerance) { return static_cast(h) * lambda; } - if (coefficient <= prefix[h] + tolerance) { - return static_cast(h + 1) * lambda + coefficient - prefix[h]; - } - } - return static_cast(p) * lambda + coefficient - prefix[p - 1]; - }; - - transformed_cut = base; - transformed_cut.rhs = -lambda; - for (i_t k = 0; k < static_cast(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] = lifting_function(base.coeff(k)); - } - } - - // Undo the extra sign-complementations used to obtain a positive knapsack - // row, then return to cuOpt's >= cut convention. - for (i_t k = 0; k < static_cast(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; -} - template f_t complemented_mixed_integer_rounding_cut_t::compute_violation( const inequality_t& cut, const std::vector& xstar) diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index e023618664..b9b772d839 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -392,6 +392,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; diff --git a/cpp/src/cuts/zero_half_mod2.cpp b/cpp/src/cuts/zero_half_mod2.cpp new file mode 100644 index 0000000000..4c8d5f1180 --- /dev/null +++ b/cpp/src/cuts/zero_half_mod2.cpp @@ -0,0 +1,715 @@ +/* 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_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) +{ + if (max_combination_size <= 0 || max_combinations <= 0) { return {}; } + + struct basis_row_t { + std::vector parity; + std::vector combination; + bool rhs{false}; + }; + + 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 {}; } + 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) { + f_t candidate_work = (f_t)(rows[candidate].parity.size() + 2); + 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]; + 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; } + + 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) +{ + if (toc(start_time) >= time_limit) { return i_t{0}; } + f_t scale_work = 0.0; + for (i_t scale = 1; scale <= max_integral_scale; ++scale) { + 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 (integral) { + work_estimate += scale_work; + return scale; + } + } + work_estimate += scale_work; + 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); + if (work_estimate > max_work_estimate) { + work_limit_reached = true; + break; + } + // 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)(10 * 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); + complemented_mir.remove_small_coefficients(lp.lower, lp.upper, transformed_cut); + if (complemented_mir.compute_violation(transformed_cut, xstar) > min_violation) { + cut_pool.add_cut(cut_type_t::ZERO_HALF, transformed_cut); + ++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) +{ + work_estimate += (f_t)(3 * oriented_aggregate.size() + 1); + inequality_t mir_cut(lp.num_cols); + if (complemented_mir.generate_cut_nonnegative_maintain_indicies( + oriented_aggregate, var_types, mir_cut)) { + 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); + if (toc(start_time) < settings.time_limit && + complemented_mir.generate_lifted_mixed_binary_cover(oriented_aggregate, + var_types, + transformed_xstar, + lifted_cover_cut, + work_estimate, + max_work_estimate)) { + 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; + const f_t max_work_estimate = work_estimate + (f_t)1e8; + bool work_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()), + &work_estimate, + max_work_estimate, + &work_limit_reached)) { + 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, + work_estimate, + max_work_estimate, + work_limit_reached); + + auto row_combinations = find_mod2_row_combinations( + candidates, max_combination_size, max_row_combinations, &work_estimate, max_work_estimate); + if (work_estimate > max_work_estimate) { work_limit_reached = true; } + scratch_pad_t aggregate_pad(lp.num_cols); + + for (const auto& combination : row_combinations) { + if (toc(start_time) >= settings.time_limit || work_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(); + } + const f_t aggregate_work = + (f_t)(4 * aggregate_input_nz + 1) + + (f_t)aggregate_input_nz * std::log2((f_t)aggregate_input_nz + (f_t)1.0); + if (add_work_estimate(aggregate_work, &work_estimate, max_work_estimate, &work_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(); + 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, + work_estimate, + max_work_estimate, + work_limit_reached, + cuts_added); + if (reversible && toc(start_time) < settings.time_limit && !work_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, + work_estimate, + max_work_estimate, + work_limit_reached, + cuts_added); + } + } + 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; } + + 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 From 6c6830476a1561b5bdfb12a5248b549e07d2619d Mon Sep 17 00:00:00 2001 From: akif Date: Thu, 6 Aug 2026 11:54:41 +0200 Subject: [PATCH 05/15] few more cleanup --- cpp/src/cuts/zero_half_mod2.cpp | 37 +++++++++++++++++++++------------ 1 file changed, 24 insertions(+), 13 deletions(-) diff --git a/cpp/src/cuts/zero_half_mod2.cpp b/cpp/src/cuts/zero_half_mod2.cpp index 4c8d5f1180..49d8381c74 100644 --- a/cpp/src/cuts/zero_half_mod2.cpp +++ b/cpp/src/cuts/zero_half_mod2.cpp @@ -49,6 +49,13 @@ struct mod2_candidate_t : mod2_parity_row_t { 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; @@ -69,13 +76,8 @@ std::vector> find_mod2_row_combinations(const std::vector parity; - std::vector combination; - bool rhs{false}; - }; + 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; @@ -96,11 +98,12 @@ std::vector> find_mod2_row_combinations(const std::vector{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; + 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())); @@ -109,7 +112,7 @@ std::vector> find_mod2_row_combinations(const std::vector combination_tmp; for (const i_t candidate : permutation) { f_t candidate_work = (f_t)(rows[candidate].parity.size() + 2); - basis_row_t current; + mod2_basis_row_t current; current.parity = rows[candidate].parity; current.combination = {candidate}; current.rhs = rows[candidate].rhs_parity; @@ -118,6 +121,7 @@ std::vector> find_mod2_row_combinations(const std::vector> find_mod2_row_combinations(const std::vector mir_cut(lp.num_cols); - if (complemented_mir.generate_cut_nonnegative_maintain_indicies( - oriented_aggregate, var_types, mir_cut)) { + 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, @@ -362,13 +368,17 @@ void mod2_generate_cuts_from_aggregate( return; } inequality_t lifted_cover_cut(lp.num_cols); - if (toc(start_time) < settings.time_limit && + 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)) { + max_work_estimate); + } + if (lifted_cover_cut_generated) { mod2_add_transformed_zero_half_cut(complemented_mir, cut_pool, lp, @@ -548,6 +558,7 @@ bool generate_mod2_zero_half_cuts(cut_pool_t& cut_pool, max_work_estimate, work_limit_reached, cuts_added); + // if the final inequality is reversable, try the reversed version as well if (reversible && toc(start_time) < settings.time_limit && !work_limit_reached) { aggregate.negate(); mod2_generate_cuts_from_aggregate(complemented_mir, From 1227a489bd24d79f99c2e018e5a125e53b739cbc Mon Sep 17 00:00:00 2001 From: akif Date: Thu, 6 Aug 2026 14:06:16 +0200 Subject: [PATCH 06/15] more clean up --- cpp/src/cuts/zero_half_mod2.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/cpp/src/cuts/zero_half_mod2.cpp b/cpp/src/cuts/zero_half_mod2.cpp index 49d8381c74..ac55131563 100644 --- a/cpp/src/cuts/zero_half_mod2.cpp +++ b/cpp/src/cuts/zero_half_mod2.cpp @@ -320,7 +320,8 @@ void mod2_add_transformed_zero_half_cut( complemented_mir.remove_small_coefficients(lp.lower, lp.upper, transformed_cut); complemented_mir.substitute_slacks(lp, Arow, transformed_cut); complemented_mir.remove_small_coefficients(lp.lower, lp.upper, transformed_cut); - if (complemented_mir.compute_violation(transformed_cut, xstar) > min_violation) { + const f_t violation = complemented_mir.compute_violation(transformed_cut, xstar); + if (violation > min_violation) { cut_pool.add_cut(cut_type_t::ZERO_HALF, transformed_cut); ++cuts_added; } From 1bfaf1800be6fee415bd35b2da873fc22faa0a7e Mon Sep 17 00:00:00 2001 From: akif Date: Fri, 7 Aug 2026 15:15:40 +0200 Subject: [PATCH 07/15] Pin papilo to the change_coefficient capacity-guard fix Clique merging could insert a coefficient into a row or column with no spare space left, tripping the changeRow size assertion during sub-MIP presolve. Repin to the fork revision that rejects an insertion only when the range is actually full, instead of when one slot remains. Signed-off-by: akif --- cpp/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index a5c8acc655..d301f0bf0c 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 From 0d39a0f44bc69ba776669f2ced0b890d9c9f77d5 Mon Sep 17 00:00:00 2001 From: akif Date: Mon, 10 Aug 2026 18:06:08 +0200 Subject: [PATCH 08/15] remove work unit related changes --- cpp/src/branch_and_bound/branch_and_bound.cpp | 133 +-------------- cpp/src/branch_and_bound/branch_and_bound.hpp | 1 - cpp/src/cuts/cuts.cpp | 157 +++--------------- cpp/src/cuts/cuts.hpp | 38 +---- cpp/src/dual_simplex/phase2.cpp | 17 -- cpp/src/utilities/work_limit_context.hpp | 5 +- cpp/tests/mip/cuts_test.cu | 19 --- 7 files changed, 34 insertions(+), 336 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 2981a759e7..4dc6bc67a8 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -69,8 +69,6 @@ using simplex::variable_type_t; namespace { -bool cut_ab_logging_enabled() { return std::getenv("CUOPT_CUT_AB_MODE") != nullptr; } - template bool is_fractional(f_t x, variable_type_t var_type, f_t integer_tol) { @@ -279,11 +277,6 @@ branch_and_bound_t::branch_and_bound_t( solver_status_(mip_status_t::UNSET) { exploration_stats_.start_time = start_time; - work_unit_context_.deterministic = settings_.deterministic; - // Root LP and cut work contributes to the public work limit, but it must not - // enter the outer heuristic/B&B horizon barrier. The heuristic may finish - // before root processing, leaving no matching barrier participant. - work_unit_context_.sync_on_horizon = false; #ifdef PRINT_CONSTRAINT_MATRIX settings_.log.printf("A"); original_problem_.A.print_matrix(); @@ -2801,7 +2794,7 @@ lp_status_t branch_and_bound_t::solve_root_relaxation( nonbasic_list, root_vstatus_, edge_norms_, - &work_unit_context_); + nullptr); } // Wait for the root relaxation solution to be sent by the diversity manager or dual simplex @@ -2964,32 +2957,12 @@ auto branch_and_bound_t::do_cut_pass( nonbasic_list, variable_bounds, exploration_stats_.start_time); - const auto& separator_work = cut_generation.last_work_stats(); - work_unit_context_.record_work_sync_on_horizon(separator_work.work_units()); - settings_.log.debug( - "Cut pass %d separator work %.4f units: gomory %.4f, knapsack %.4f, flow %.4f, " - "MIR %.4f, implied %.4f, graph %.4f, clique %.4f, zero-half %.4f\n", - cut_pass, - separator_work.work_units(), - separator_work.gomory / 1e8, - separator_work.knapsack / 1e8, - separator_work.flow_cover / 1e8, - separator_work.mir / 1e8, - separator_work.implied_bound / 1e8, - separator_work.conflict_graph / 1e8, - separator_work.clique / 1e8, - separator_work.zero_half / 1e8); if (!problem_feasible) { if (settings_.heuristic_preemption_callback != nullptr) { settings_.heuristic_preemption_callback(); } return {cut_pass_action_t::RETURN, mip_status_t::INFEASIBLE}; } - if (work_unit_context_.global_work_units_elapsed >= settings_.work_limit) { - solver_status_ = mip_status_t::WORK_LIMIT; - set_final_solution(solution, root_objective_); - return {cut_pass_action_t::RETURN, solver_status_}; - } if (toc(exploration_stats_.start_time) >= settings_.time_limit) { solver_status_ = mip_status_t::TIME_LIMIT; set_final_solution(solution, root_objective_); @@ -3001,48 +2974,15 @@ auto branch_and_bound_t::do_cut_pass( } // Score the cuts f_t score_start_time = tic(); - constexpr f_t max_cut_scoring_work = static_cast(1e8); - const f_t cut_scoring_work = cut_pool.score_cuts(root_relax_soln_.x, max_cut_scoring_work); - work_unit_context_.record_work_sync_on_horizon(cut_scoring_work / static_cast(1e8)); + cut_pool.score_cuts(root_relax_soln_.x); f_t score_time = toc(score_start_time); if (score_time > 1.0) { settings_.log.debug("Cut scoring time %.2f seconds\n", score_time); } - settings_.log.debug( - "Cut pass %d scoring work %.4f units\n", cut_pass, cut_scoring_work / static_cast(1e8)); - settings_.log.printf( - "Cut pass %d work units: total %.4f, separation %.4f " - "(Gomory %.4f, knapsack %.4f, flow %.4f, MIR %.4f, implied %.4f, graph %.4f, " - "clique %.4f, zero-half %.4f), scoring %.4f\n", - cut_pass, - separator_work.work_units() + cut_scoring_work / static_cast(1e8), - separator_work.work_units(), - separator_work.gomory / 1e8, - separator_work.knapsack / 1e8, - separator_work.flow_cover / 1e8, - separator_work.mir / 1e8, - separator_work.implied_bound / 1e8, - separator_work.conflict_graph / 1e8, - separator_work.clique / 1e8, - separator_work.zero_half / 1e8, - cut_scoring_work / static_cast(1e8)); - if (work_unit_context_.global_work_units_elapsed >= settings_.work_limit) { - solver_status_ = mip_status_t::WORK_LIMIT; - set_final_solution(solution, root_objective_); - return {cut_pass_action_t::RETURN, solver_status_}; - } // Get the best cuts from the cut pool csr_matrix_t cuts_to_add(0, original_lp_.num_cols, 0); std::vector cut_rhs; std::vector cut_types; i_t num_cuts = cut_pool.get_best_cuts(cuts_to_add, cut_rhs, cut_types); if (num_cuts == 0) { return {cut_pass_action_t::BREAK, mip_status_t::UNSET}; } - const f_t cut_assembly_work = - static_cast(5 * cuts_to_add.row_start[cuts_to_add.m] + 4 * num_cuts); - work_unit_context_.record_work_sync_on_horizon(cut_assembly_work / static_cast(1e8)); - if (work_unit_context_.global_work_units_elapsed >= settings_.work_limit) { - solver_status_ = mip_status_t::WORK_LIMIT; - set_final_solution(solution, root_objective_); - return {cut_pass_action_t::RETURN, solver_status_}; - } cut_info.record_cut_types(cut_types); #ifdef PRINT_CUT_POOL_TYPES cut_pool.print_cutpool_types(); @@ -3124,8 +3064,6 @@ auto branch_and_bound_t::do_cut_pass( std::vector new_upper = original_lp_.upper; bool feasible = node_presolve.bounds_strengthening(settings_, bounds_changed, new_lower, new_upper); - work_unit_context_.record_work_sync_on_horizon(node_presolve.last_nnz_processed / - static_cast(1e8)); mutex_original_lp_.lock(); original_lp_.lower = new_lower; original_lp_.upper = new_upper; @@ -3142,12 +3080,6 @@ auto branch_and_bound_t::do_cut_pass( return {cut_pass_action_t::RETURN, mip_status_t::INFEASIBLE}; } - if (work_unit_context_.global_work_units_elapsed >= settings_.work_limit) { - solver_status_ = mip_status_t::WORK_LIMIT; - set_final_solution(solution, root_objective_); - return {cut_pass_action_t::RETURN, solver_status_}; - } - if (toc(exploration_stats_.start_time) >= settings_.time_limit) { solver_status_ = mip_status_t::TIME_LIMIT; set_final_solution(solution, root_objective_); @@ -3170,8 +3102,7 @@ auto branch_and_bound_t::do_cut_pass( nonbasic_list, root_relax_soln_, iter, - edge_norms_, - &work_unit_context_); + edge_norms_); exploration_stats_.total_simplex_iters += iter; f_t dual_phase2_time = toc(dual_phase2_start_time); if (dual_phase2_time > 1.0) { @@ -3182,11 +3113,6 @@ auto branch_and_bound_t::do_cut_pass( set_final_solution(solution, root_objective_); return {cut_pass_action_t::RETURN, solver_status_}; } - if (cut_status == dual_status_t::WORK_LIMIT) { - solver_status_ = mip_status_t::WORK_LIMIT; - set_final_solution(solution, root_objective_); - return {cut_pass_action_t::RETURN, solver_status_}; - } if (cut_status != dual_status_t::OPTIMAL) { settings_.log.printf("Numerical issue at root node. Resolving from scratch\n"); @@ -3199,17 +3125,12 @@ auto branch_and_bound_t::do_cut_pass( basic_list, nonbasic_list, root_vstatus_, - edge_norms_, - &work_unit_context_); + edge_norms_); if (scratch_status == lp_status_t::OPTIMAL) { // We recovered cut_status = convert_lp_status_to_dual_status(scratch_status); exploration_stats_.total_simplex_iters += root_relax_soln_.iterations; root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); - } else if (scratch_status == lp_status_t::WORK_LIMIT) { - solver_status_ = mip_status_t::WORK_LIMIT; - set_final_solution(solution, root_objective_); - return {cut_pass_action_t::RETURN, solver_status_}; } else { settings_.log.printf("Cut status %s\n", simplex::dual_status_to_string(cut_status).c_str()); #ifdef WRITE_CUT_INFEASIBLE_MPS @@ -3380,8 +3301,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut basic_list, nonbasic_list, root_vstatus_, - edge_norms_, - &work_unit_context_); + edge_norms_); root_relax_solved_by = DualSimplex; exploration_stats_.total_simplex_iters = root_relax_soln_.iterations; @@ -3447,15 +3367,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut return solver_status_; } - if (work_unit_context_.global_work_units_elapsed >= settings_.work_limit) { - settings_.log.printf("\n"); - solver_status_ = mip_status_t::WORK_LIMIT; - set_final_solution(solution, -inf); - signal_extend_cliques_.store(true, std::memory_order_release); -#pragma omp taskwait depend(in : *clique_signal) - return solver_status_; - } - assert(root_status == lp_status_t::OPTIMAL); settings_.log.printf("\n"); settings_.log.print_format("Root relaxation solution found in {} iterations and {:.2f}s by {}\n", @@ -3468,10 +3379,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut set_uninitialized_steepest_edge_norms(original_lp_, basic_list, edge_norms_); root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); - if (cut_ab_logging_enabled()) { - settings_.log.printf("CUT_AB_ROOT_NO_CUTS %.17g\n", - compute_user_objective(original_lp_, root_objective_)); - } if (settings_.set_simplex_solution_callback != nullptr) { std::vector original_x; @@ -3494,10 +3401,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut cut_info_t cut_info; if (num_fractional == 0) { - if (cut_ab_logging_enabled()) { - settings_.log.printf("CUT_AB_ROOT_WITH_CUTS %.17g\n", - compute_user_objective(original_lp_, root_objective_)); - } if (settings_.benchmark_info_ptr != nullptr) { const double v = static_cast(compute_user_objective(original_lp_, root_objective_)); settings_.benchmark_info_ptr->root_lp_no_cuts = v; @@ -3572,10 +3475,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut settings_.benchmark_info_ptr->root_lp_with_cuts = compute_user_objective(original_lp_, root_objective_); } - if (cut_ab_logging_enabled()) { - settings_.log.printf("CUT_AB_ROOT_WITH_CUTS %.17g\n", - compute_user_objective(original_lp_, root_objective_)); - } set_solution_at_root(solution, cut_info); if (settings_.benchmark_info_ptr != nullptr) { settings_.benchmark_info_ptr->cut_generation_time_sec = toc(cut_generation_start_time); @@ -3609,10 +3508,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut root_fj_cpu_worker.stop(); if (cut_pass_result.action == cut_pass_action_t::RETURN) { - if (cut_ab_logging_enabled()) { - settings_.log.printf("CUT_AB_ROOT_WITH_CUTS %.17g\n", - compute_user_objective(original_lp_, root_objective_)); - } if (settings_.benchmark_info_ptr != nullptr) { settings_.benchmark_info_ptr->cut_generation_time_sec = toc(cut_generation_start_time); } @@ -3638,10 +3533,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut settings_.benchmark_info_ptr->root_lp_with_cuts = compute_user_objective(original_lp_, root_objective_); } - if (cut_ab_logging_enabled()) { - settings_.log.printf("CUT_AB_ROOT_WITH_CUTS %.17g\n", - compute_user_objective(original_lp_, root_objective_)); - } print_cut_info(settings_, cut_info); f_t cut_generation_time = toc(cut_generation_start_time); @@ -3672,16 +3563,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut settings_.log.printf("\n"); } - // Benchmark-only root stop: CUT_AB measurements compare separator families - // at the completed root LP, before strong branching or tree exploration. - if (cut_ab_logging_enabled()) { - solver_status_ = mip_status_t::NODE_LIMIT; - set_final_solution(solution, root_objective_); - signal_extend_cliques_.store(true, std::memory_order_release); -#pragma omp taskwait depend(in : *clique_signal) - return solver_status_; - } - if (enable_root_cut_cpufj && cut_info.has_cuts()) { f_t root_cut_cpufj_build_start_time = tic(); // In deterministic mode this CPUFJ is built on the B&B task while the LS deterministic @@ -3821,8 +3702,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut } print_table_header(); - deterministic_root_work_offset_ = work_unit_context_.global_work_units_elapsed; - #pragma omp taskgroup { if (settings_.deterministic) { @@ -4310,7 +4189,7 @@ void branch_and_bound_t::deterministic_sync_callback() } // Stop early if next horizon exceeds work limit - if (deterministic_root_work_offset_ + deterministic_current_horizon_ > settings_.work_limit) { + if (deterministic_current_horizon_ > settings_.work_limit) { deterministic_global_termination_status_ = mip_status_t::WORK_LIMIT; } diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index a33c6c84b8..96b8a6d8fe 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -487,7 +487,6 @@ class branch_and_bound_t { mip_status_t deterministic_global_termination_status_{mip_status_t::UNSET}; double deterministic_horizon_step_{5.0}; // Work unit step per horizon (tunable) double deterministic_current_horizon_{0.0}; // Current horizon target - double deterministic_root_work_offset_{0.0}; bool deterministic_mode_enabled_{false}; int deterministic_horizon_number_{0}; // Current horizon number (for debugging) diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index 982764c4a8..0a518ed2da 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -1226,21 +1226,10 @@ f_t cut_pool_t::cut_orthogonality(i_t i, i_t j) template void cut_pool_t::check_for_duplicate_cuts() -{ - f_t work_estimate = 0.0; - check_for_duplicate_cuts(work_estimate, std::numeric_limits::infinity()); -} - -template -bool cut_pool_t::check_for_duplicate_cuts(f_t& work_estimate, f_t max_work_estimate) { // 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 - const f_t setup_work = static_cast(5 * cut_storage_.m + 2 * cut_storage_.n) + - static_cast(4 * cut_storage_.row_start[cut_storage_.m]); - if (add_work_estimate(setup_work, &work_estimate, max_work_estimate)) { return false; } - std::vector divisors(cut_storage_.m, 0.0); std::vector sets(cut_storage_.m, 0); @@ -1271,10 +1260,6 @@ bool cut_pool_t::check_for_duplicate_cuts(f_t& work_estimate, f_t max_ new_rows++; } else if (sets[r] < new_set_0) { // Look over indices a_ij with i > r - if (add_work_estimate( - static_cast(6 * (col_end - (p + 1))), &work_estimate, max_work_estimate)) { - return false; - } 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]; @@ -1315,10 +1300,6 @@ bool cut_pool_t::check_for_duplicate_cuts(f_t& work_estimate, f_t max_ const i_t set_r = sets[r]; if (set_r > 0 && set_r < sentinel && cuts_to_remove[r] == 0) { // This cut has a duplicate - if (add_work_estimate( - static_cast(5 * (m - (r + 1))), &work_estimate, max_work_estimate)) { - return false; - } for (i_t i = r + 1; i < m; i++) { if (sets[i] == set_r) { const f_t f_r = divisors[r]; @@ -1373,24 +1354,12 @@ bool cut_pool_t::check_for_duplicate_cuts(f_t& work_estimate, f_t max_ cut_type_.resize(write); cut_age_.resize(write); } - return true; } template -f_t cut_pool_t::score_cuts(std::vector& x_relax, f_t max_work_estimate) +void cut_pool_t::score_cuts(std::vector& x_relax) { - f_t work_estimate = 0.0; - best_cuts_.clear(); - scored_cuts_ = 0; - - const f_t duplicate_work_limit = std::isfinite(max_work_estimate) - ? max_work_estimate / static_cast(2.0) - : max_work_estimate; - check_for_duplicate_cuts(work_estimate, duplicate_work_limit); - - const f_t distance_work = static_cast(5 * cut_storage_.row_start[cut_storage_.m]) + - static_cast(3 * cut_storage_.m); - if (add_work_estimate(distance_work, &work_estimate, max_work_estimate)) { return work_estimate; } + check_for_duplicate_cuts(); cut_distances_.resize(cut_storage_.m, 0.0); cut_norms_.resize(cut_storage_.m, 0.0); @@ -1410,14 +1379,13 @@ f_t cut_pool_t::score_cuts(std::vector& x_relax, f_t max_work_est } std::vector sorted_indices; - const f_t sort_work = static_cast(cut_storage_.m) * - std::log2(static_cast(cut_storage_.m) + static_cast(1.0)); - if (add_work_estimate(sort_work, &work_estimate, max_work_estimate)) { return work_estimate; } best_score_last_permutation(cut_distances_, sorted_indices); const i_t max_cuts = 2000; const f_t min_orthogonality = settings_.cut_min_orthogonality; best_cuts_.reserve(std::min(max_cuts, cut_storage_.m)); + best_cuts_.clear(); + scored_cuts_ = 0; if (!sorted_indices.empty()) { const i_t i = sorted_indices.back(); @@ -1426,8 +1394,7 @@ f_t cut_pool_t::score_cuts(std::vector& x_relax, f_t max_work_est scored_cuts_++; } - bool work_limit_reached = false; - while (scored_cuts_ < max_cuts && !sorted_indices.empty() && !work_limit_reached) { + while (scored_cuts_ < max_cuts && !sorted_indices.empty()) { const i_t i = sorted_indices.back(); sorted_indices.pop_back(); @@ -1436,22 +1403,14 @@ f_t cut_pool_t::score_cuts(std::vector& x_relax, f_t max_work_est f_t cut_ortho = 1.0; const i_t best_cuts_size = best_cuts_.size(); for (i_t k = 0; k < best_cuts_size; k++) { - const i_t j = best_cuts_[k]; - const i_t i_nz = cut_storage_.row_start[i + 1] - cut_storage_.row_start[i]; - const i_t j_nz = cut_storage_.row_start[j + 1] - cut_storage_.row_start[j]; - if (add_work_estimate( - static_cast(4 * (i_nz + j_nz)), &work_estimate, max_work_estimate)) { - work_limit_reached = true; - break; - } - cut_ortho = std::min(cut_ortho, cut_orthogonality(i, j)); + const i_t j = best_cuts_[k]; + cut_ortho = std::min(cut_ortho, cut_orthogonality(i, j)); } - if (!work_limit_reached && cut_ortho >= min_orthogonality) { + if (cut_ortho >= min_orthogonality) { best_cuts_.push_back(i); scored_cuts_++; } } - return work_estimate; } template @@ -3136,12 +3095,6 @@ void cut_generation_t::generate_implied_bound_cuts( { if (probing_implied_bound_.zero_offsets.empty()) { return; } - last_work_stats_.implied_bound += - static_cast( - 4 * std::min(lp.num_cols, static_cast(probing_implied_bound_.zero_offsets.size()) - 1)) + - static_cast(20 * (probing_implied_bound_.zero_variables.size() + - probing_implied_bound_.one_variables.size())); - 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; @@ -3312,8 +3265,8 @@ void cut_generation_t::prepare_fractional_sub_conflict_graph( } const f_t bound_tol = settings.primal_tol; - f_t& work_estimate = last_work_stats_.conflict_graph; - const f_t max_work_estimate = work_estimate + static_cast(1e7); + f_t work_estimate = 0.0; + const f_t max_work_estimate = 1e7; sub_cg_.num_vars = num_vars; sub_cg_.vertices.reserve(static_cast(num_vars) * 2); @@ -3557,8 +3510,6 @@ bool cut_generation_t::generate_cuts(const lp_problem_t& lp, variable_bounds_t& variable_bounds, f_t start_time) { - last_work_stats_ = {}; - // Generate Gomory and CG Cuts if (settings.mixed_integer_gomory_cuts != 0 || settings.strong_chvatal_gomory_cuts != 0) { if (toc(start_time) >= settings.time_limit) { return true; } @@ -3685,9 +3636,6 @@ void cut_generation_t::generate_knapsack_cuts( if (knapsack_generation_.num_knapsack_constraints() > 0) { for (i_t knapsack_row : knapsack_generation_.get_knapsack_constraints()) { if (toc(start_time) >= settings.time_limit) { return; } - const f_t row_size = static_cast(Arow.row_length(knapsack_row)); - last_work_stats_.knapsack += - static_cast(20.0) * row_size + row_size * std::log2(row_size + static_cast(1.0)); 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); @@ -3707,12 +3655,8 @@ void cut_generation_t::generate_flow_cover_cuts( f_t start_time) { if (flow_cover_generation_.num_constraints() > 0) { - last_work_stats_.flow_cover += static_cast(4 * lp.num_cols); for (const auto& flow_cover_row : flow_cover_generation_.get_constraints()) { if (toc(start_time) >= settings.time_limit) { return; } - const f_t row_size = static_cast(Arow.row_length(flow_cover_row.row)); - last_work_stats_.flow_cover += - static_cast(30.0) * row_size + row_size * std::log2(row_size + static_cast(1.0)); inequality_t cut(lp.num_cols); i_t status = flow_cover_generation_.generate_cut( lp, settings, Arow, variable_bounds, var_types, xstar, flow_cover_row, cut); @@ -3761,8 +3705,8 @@ bool cut_generation_t::generate_clique_cuts( const f_t min_weight = 1.0 + min_violation; // TODO this can be problem dependent const i_t max_calls = 100000; - f_t& work_estimate = last_work_stats_.clique; - const f_t max_work_estimate = work_estimate + static_cast(1e8); + f_t work_estimate = 0.0; + const f_t max_work_estimate = 1e8; const std::vector& vertices = sub_cg_.vertices; const std::vector& weights = sub_cg_.weights; @@ -3942,6 +3886,7 @@ bool cut_generation_t::generate_zero_half_cuts( static_cast(sub_cg_.ready), sub_cg_.vertices.size()); + f_t mod2_work_estimate = 0.0; if (!generate_mod2_zero_half_cuts(cut_pool_, lp, settings, @@ -3951,7 +3896,7 @@ bool cut_generation_t::generate_zero_half_cuts( xstar, variable_bounds, start_time, - last_work_stats_.zero_half)) { + mod2_work_estimate)) { return true; } @@ -3982,8 +3927,8 @@ bool cut_generation_t::generate_zero_half_cuts( 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 = last_work_stats_.zero_half; - const f_t max_work_estimate = work_estimate + static_cast(1e8); + f_t work_estimate = 0.0; + const f_t max_work_estimate = 1e8; const std::vector& vertices = sub_cg_.vertices; const std::vector& weights = sub_cg_.weights; @@ -4149,11 +4094,8 @@ void cut_generation_t::generate_mir_cuts( 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 = last_work_stats_.mir; - work_estimate += static_cast(4 * lp.num_cols + 3 * lp.num_rows) + - static_cast(4 * Arow.row_start[Arow.m]); - const f_t max_work_estimate = work_estimate + static_cast(2e9); - i_t num_cuts = 0; + f_t work_estimate = 0.0; + i_t num_cuts = 0; while (num_cuts < max_cuts && !score_queue.empty()) { if (toc(start_time) >= settings.time_limit) { break; } // Get the row with the highest score from the queue @@ -4171,7 +4113,7 @@ void cut_generation_t::generate_mir_cuts( const f_t slack_value = xstar[slack]; if (max_score <= 0.0) { break; } - if (work_estimate > max_work_estimate) { break; } + if (work_estimate > 2e9) { break; } inequality_t inequality(Arow, i, lp.rhs[i]); work_estimate += inequality.size(); @@ -4403,14 +4345,6 @@ void cut_generation_t::generate_gomory_cuts( const std::vector& nonbasic_list, f_t start_time) { - f_t& work_estimate = last_work_stats_.gomory; - const f_t max_work_estimate = work_estimate + static_cast(1e8); - const f_t predicted_setup_work = - static_cast(2.0) * static_cast(lp.num_rows) * static_cast(lp.num_rows) + - static_cast(12 * lp.num_cols) + static_cast(5 * Arow.row_start[Arow.m]); - if (add_work_estimate(predicted_setup_work, &work_estimate, max_work_estimate)) { return; } - - const f_t basis_setup_start_work = basis_update.work_estimate(); tableau_equality_t tableau(lp, basis_update, nonbasic_list); mixed_integer_gomory_cut_t gomory_cut; complemented_mixed_integer_rounding_cut_t complemented_mir(lp, settings, new_slacks); @@ -4421,38 +4355,17 @@ void cut_generation_t::generate_gomory_cuts( std::vector transformed_xstar; complemented_mir.bound_substitution(lp, variable_bounds, var_types, xstar, transformed_xstar); - const f_t actual_setup_work = basis_update.work_estimate() - basis_setup_start_work; - if (actual_setup_work > predicted_setup_work) { - work_estimate += actual_setup_work - predicted_setup_work; - } - work_estimate += - static_cast(4 * lp.num_rows) + static_cast(variable_bounds.upper_variables.size() + - variable_bounds.lower_variables.size()); - if (work_estimate > max_work_estimate) { return; } - for (i_t i = 0; i < lp.num_rows; i++) { - if (toc(start_time) >= settings.time_limit || work_estimate > max_work_estimate) { break; } + if (toc(start_time) >= settings.time_limit) { break; } inequality_t inequality(lp.num_cols); const i_t j = basic_list[i]; if (var_types[j] != variable_type_t::INTEGER) { continue; } const f_t x_j = xstar[j]; if (fractional_part(x_j) < 0.05 || fractional_part(x_j) > 0.95) { continue; } - i_t tableau_status = tableau.generate_base_equality(lp, - settings, - Arow, - var_types, - basis_update, - xstar, - basic_list, - nonbasic_list, - i, - inequality, - work_estimate, - max_work_estimate); + i_t tableau_status = tableau.generate_base_equality( + lp, settings, Arow, var_types, basis_update, xstar, basic_list, nonbasic_list, i, inequality); if (tableau_status == 0) { - const f_t cut_work = static_cast(120 * inequality.size() + 10); - if (add_work_estimate(cut_work, &work_estimate, max_work_estimate)) { break; } // Generate a CG cut const bool generate_cg_cut = settings.strong_chvatal_gomory_cuts != 0; if (generate_cg_cut) { @@ -4562,9 +4475,7 @@ i_t tableau_equality_t::generate_base_equality( const std::vector& basic_list, const std::vector& nonbasic_list, i_t i, - inequality_t& inequality, - f_t& work_estimate, - f_t max_work_estimate) + inequality_t& inequality) { // Let's look for Gomory cuts const i_t j = basic_list[i]; @@ -4580,16 +4491,7 @@ i_t tableau_equality_t::generate_base_equality( e_i.i[0] = i; e_i.x[0] = 1.0; sparse_vector_t u_bar(lp.num_rows, 0); - const f_t predicted_basis_work = static_cast(3 * lp.num_rows + 4); - if (add_work_estimate(predicted_basis_work, &work_estimate, max_work_estimate)) { return -2; } - const f_t basis_work_start = basis_update.work_estimate(); basis_update.b_transpose_solve(e_i, u_bar); - const f_t actual_basis_work = basis_update.work_estimate() - basis_work_start; - if (actual_basis_work > predicted_basis_work && - add_work_estimate( - actual_basis_work - predicted_basis_work, &work_estimate, max_work_estimate)) { - return -2; - } #ifdef CHECK_B_TRANSPOSE_SOLVE std::vector u_bar_dense(lp.num_rows); @@ -4616,12 +4518,7 @@ i_t tableau_equality_t::generate_base_equality( // Compute a_bar = N^T u_bar // TODO: This is similar to a function in phase2 of dual simplex. See if it can be reused. - const i_t nz_ubar = u_bar.i.size(); - f_t tableau_multiply_work = static_cast(3 * nz_ubar + 2); - for (const i_t row : u_bar.i) { - tableau_multiply_work += static_cast(6 * (Arow.row_start[row + 1] - Arow.row_start[row])); - } - if (add_work_estimate(tableau_multiply_work, &work_estimate, max_work_estimate)) { return -2; } + const i_t nz_ubar = u_bar.i.size(); std::vector abar_indices; abar_indices.reserve(nz_ubar); for (i_t k = 0; k < nz_ubar; k++) { @@ -4682,11 +4579,7 @@ i_t tableau_equality_t::generate_base_equality( // Check that the tableau equality is satisfied const f_t tableau_tol = 1e-6; - if (add_work_estimate( - static_cast(4 * a_bar.i.size() + 2), &work_estimate, max_work_estimate)) { - return -2; - } - f_t a_bar_dot_xstar = a_bar.dot(xstar); + f_t a_bar_dot_xstar = a_bar.dot(xstar); if (std::abs(a_bar_dot_xstar - b_bar_[i]) > tableau_tol) { settings.log.debug("bad tableau equality. error %e\n", std::abs(a_bar_dot_xstar - b_bar_[i])); return -1; diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index b9b772d839..c052935806 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -17,7 +17,6 @@ #include #include #include -#include #include #include #include @@ -55,29 +54,6 @@ struct cut_gap_closure_t { f_t gap_closed_ratio{0.0}; }; -// Deterministic operation estimates for one separator pass. Raw estimates are -// kept per phase so expensive cut families remain attributable; callers can -// convert the total to solver work units with total() / 1e8. -template -struct cut_work_stats_t { - f_t gomory{0.0}; - f_t knapsack{0.0}; - f_t flow_cover{0.0}; - f_t mir{0.0}; - f_t implied_bound{0.0}; - f_t conflict_graph{0.0}; - f_t clique{0.0}; - f_t zero_half{0.0}; - - f_t total() const - { - return gomory + knapsack + flow_cover + mir + implied_bound + conflict_graph + clique + - zero_half; - } - - f_t work_units() const { return total() / static_cast(1e8); } -}; - template cut_gap_closure_t compute_cut_gap_closure(f_t objective_reference, f_t objective_before_cuts, @@ -319,7 +295,6 @@ std::vector> find_mod2_row_combinations_for_test( 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, @@ -346,10 +321,7 @@ class cut_pool_t { // We expect that the cut is violated by the current relaxation xstar. void add_cut(cut_type_t cut_type, const inequality_t& cut); - // Returns a deterministic operation estimate. Orthogonality selection is - // stopped at max_work_estimate, preserving the best cuts found so far. - f_t score_cuts(std::vector& x_relax, - f_t max_work_estimate = std::numeric_limits::infinity()); + void score_cuts(std::vector& x_relax); // We return the cuts in the form best_cuts*x <= best_rhs i_t get_best_cuts(csr_matrix_t& best_cuts, @@ -365,7 +337,6 @@ class cut_pool_t { void print_cutpool_types() { print_cut_types("In cut pool", cut_type_, settings_); } void check_for_duplicate_cuts(); - bool check_for_duplicate_cuts(f_t& work_estimate, f_t max_work_estimate); private: f_t cut_distance(i_t row, const std::vector& x, f_t& cut_violation, f_t& cut_norm); @@ -717,8 +688,6 @@ class cut_generation_t { variable_bounds_t& variable_bounds, f_t start_time); - const cut_work_stats_t& last_work_stats() const { return last_work_stats_; } - private: // Generate all mixed integer gomory cuts void generate_gomory_cuts(const simplex::lp_problem_t& lp, @@ -801,7 +770,6 @@ class cut_generation_t { std::shared_ptr> clique_table_; omp_atomic_t* signal_extend_{nullptr}; fractional_conflict_subgraph_t sub_cg_; - cut_work_stats_t last_work_stats_; }; template @@ -890,9 +858,7 @@ class tableau_equality_t { const std::vector& basic_list, const std::vector& nonbasic_list, i_t i, - inequality_t& inequality, - f_t& work_estimate, - f_t max_work_estimate); + inequality_t& inequality); private: std::vector b_bar_; diff --git a/cpp/src/dual_simplex/phase2.cpp b/cpp/src/dual_simplex/phase2.cpp index 8c11bec81b..a5f10c3229 100644 --- a/cpp/src/dual_simplex/phase2.cpp +++ b/cpp/src/dual_simplex/phase2.cpp @@ -2559,23 +2559,6 @@ dual_status_t dual_phase2_with_advanced_basis(i_t phase, f_t phase2_work_estimate = 0.0; ft.clear_work_estimate(); - struct final_work_flush_t { - work_limit_context_t* context; - f_t& phase_work; - basis_update_mpf_t& basis_update; - - ~final_work_flush_t() - { - if (context == nullptr) { return; } - phase_work += basis_update.work_estimate(); - basis_update.clear_work_estimate(); - if (phase_work > static_cast(0.0)) { - context->record_work_sync_on_horizon(phase_work / static_cast(1e8)); - phase_work = 0.0; - } - } - } final_work_flush{work_unit_context, phase2_work_estimate, ft}; - std::vector& x = sol.x; std::vector& y = sol.y; std::vector& z = sol.z; diff --git a/cpp/src/utilities/work_limit_context.hpp b/cpp/src/utilities/work_limit_context.hpp index f7e2b537f1..3463847123 100644 --- a/cpp/src/utilities/work_limit_context.hpp +++ b/cpp/src/utilities/work_limit_context.hpp @@ -18,7 +18,6 @@ struct work_limit_context_t { double global_work_units_elapsed{0.0}; double total_sync_time{0.0}; // Total time spent waiting at sync barriers (seconds) bool deterministic{false}; - bool sync_on_horizon{true}; work_unit_scheduler_t* scheduler{nullptr}; std::string name; @@ -28,9 +27,7 @@ struct work_limit_context_t { { if (!deterministic) return; global_work_units_elapsed += work; - if (scheduler && sync_on_horizon) { - scheduler->on_work_recorded(*this, global_work_units_elapsed); - } + if (scheduler) { scheduler->on_work_recorded(*this, global_work_units_elapsed); } } }; diff --git a/cpp/tests/mip/cuts_test.cu b/cpp/tests/mip/cuts_test.cu index 08987dd687..5af0754e3d 100644 --- a/cpp/tests/mip/cuts_test.cu +++ b/cpp/tests/mip/cuts_test.cu @@ -981,22 +981,6 @@ TEST(cuts, test_duplicate_cuts_detection) cut_pool.add_cut(mip::cut_type_t::MIXED_INTEGER_GOMORY, cut8); cut_pool.check_for_duplicate_cuts(); - - std::vector xstar(4, 0.0); - const double scoring_work = cut_pool.score_cuts(xstar, 1.0); - EXPECT_GT(scoring_work, 1.0); - EXPECT_TRUE(std::isfinite(scoring_work)); -} - -TEST(cuts, cut_work_stats_total) -{ - mip::cut_work_stats_t stats; - stats.gomory = 1.0; - stats.mir = 2.0; - stats.zero_half = 3.0; - stats.conflict_graph = 4.0; - EXPECT_DOUBLE_EQ(stats.total(), 10.0); - EXPECT_DOUBLE_EQ(stats.work_units(), 1e-7); } TEST(cuts, clique_phase1_smoke_conflict_graph_edges) @@ -1539,9 +1523,6 @@ TEST(cuts, zero_half_unit_mod2_row_finder_stops_at_work_limit) const std::vector> parity_rows(256, support); const std::vector rhs_parity(256, 0); - // The budget admits preprocessing and the first basis row, then expires at - // the next symmetric difference. The search must return immediately rather - // than repeating an over-budget reduction for every remaining input row. constexpr double max_work = 18900.0; double work = 0.0; const auto combinations = From 25f49238d7caf80e0a0e246c775eb817afa00c9b Mon Sep 17 00:00:00 2001 From: akif Date: Tue, 11 Aug 2026 15:32:27 +0200 Subject: [PATCH 09/15] bound cut generation and B&B timeout work Limit cut pool growth and mod-2 work so expensive cut phases remain bounded, while broadcasting B&B timeouts to active node solves. --- cpp/src/branch_and_bound/branch_and_bound.cpp | 36 +++-- cpp/src/branch_and_bound/branch_and_bound.hpp | 2 + cpp/src/cuts/cuts.cpp | 128 +++++++++++++----- cpp/src/cuts/cuts.hpp | 16 ++- cpp/src/cuts/zero_half_mod2.cpp | 98 ++++++++++---- cpp/tests/mip/cuts_test.cu | 39 ++++++ 6 files changed, 244 insertions(+), 75 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 4dc6bc67a8..462ba43b0f 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; } @@ -2930,6 +2936,8 @@ auto branch_and_bound_t::do_cut_pass( f_t& last_objective, f_t root_relax_objective, i_t& cut_pool_size, + f_t& cut_scoring_time, + i_t& max_scoring_pool_size, [[maybe_unused]] const std::vector& saved_solution) -> cut_pass_result_t { #ifdef PRINT_FRACTIONAL_INFO @@ -2973,10 +2981,10 @@ auto branch_and_bound_t::do_cut_pass( settings_.log.debug("Cut generation time %.2f seconds\n", cut_generation_time); } // Score the cuts - f_t score_start_time = tic(); + max_scoring_pool_size = std::max(max_scoring_pool_size, cut_pool.pool_size()); + f_t score_start_time = tic(); cut_pool.score_cuts(root_relax_soln_.x); - f_t score_time = toc(score_start_time); - if (score_time > 1.0) { settings_.log.debug("Cut scoring time %.2f seconds\n", score_time); } + cut_scoring_time += toc(score_start_time); // Get the best cuts from the cut pool csr_matrix_t cuts_to_add(0, original_lp_.num_cols, 0); std::vector cut_rhs; @@ -3456,6 +3464,8 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut f_t cut_generation_start_time = tic(); i_t cut_pool_size = 0; + f_t cut_scoring_time = 0.0; + i_t max_scoring_pool_size = 0; for (i_t cut_pass = 0; cut_pass < settings_.max_cut_passes; cut_pass++) { if (toc(exploration_stats_.start_time) >= settings_.time_limit) { solver_status_ = mip_status_t::TIME_LIMIT; @@ -3504,6 +3514,8 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut last_objective, root_relax_objective, cut_pool_size, + cut_scoring_time, + max_scoring_pool_size, saved_solution); root_fj_cpu_worker.stop(); @@ -3554,7 +3566,9 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut mutex_upper_.unlock(); settings_.log.printf("Cut generation time: %.2f seconds\n", cut_generation_time); - settings_.log.printf("Cut pool size : %d\n", cut_pool_size); + settings_.log.printf("Cut scoring time : %.2f seconds\n", cut_scoring_time); + settings_.log.printf("Cut scoring max pool: %d\n", max_scoring_pool_size); + settings_.log.printf("Cut pool size : %d\n", cut_pool_size); settings_.log.printf("Size with cuts : %d constraints, %d variables, %d nonzeros\n", original_lp_.num_rows, original_lp_.num_cols, diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index 96b8a6d8fe..6612090952 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -323,6 +323,8 @@ class branch_and_bound_t { f_t& last_objective, f_t root_relax_objective, i_t& cut_pool_size, + f_t& cut_scoring_time, + i_t& max_scoring_pool_size, const std::vector& saved_solution); // Set the solution when found at the root node diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index 0a518ed2da..97c5922c47 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -1146,6 +1146,8 @@ inline uint64_t hash64_with_seed(uint64_t value, uint64_t seed) template void cut_pool_t::add_cut(cut_type_t cut_type, const inequality_t& cut) { + if (generation_limit_reached(cut_type)) { return; } + // TODO: Add fast duplicate check and only add if the cut is not already in the pool for (i_t p = 0; p < cut.size(); p++) { @@ -1168,6 +1170,7 @@ void cut_pool_t::add_cut(cut_type_t cut_type, const inequality_t @@ -3095,18 +3098,30 @@ 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; for (i_t j = 0; j < n_cols; j++) { + if (toc(start_time) >= settings.time_limit || + cut_pool_.generation_limit_reached(IMPLIED_BOUND)) { + 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]; + if (add_work_estimate( + (f_t)(zero_end - zero_begin + one_end - one_begin), &work_estimate, max_work_estimate)) { + return; + } for (i_t p = zero_begin; p < zero_end; p++) { const i_t i = probing_implied_bound_.zero_variables[p]; if (i == j) { continue; } @@ -3151,8 +3166,6 @@ void cut_generation_t::generate_implied_bound_cuts( } // 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++) { const i_t i = probing_implied_bound_.one_variables[p]; if (i == j) { continue; } @@ -3514,16 +3527,22 @@ bool cut_generation_t::generate_cuts(const lp_problem_t& lp, if (settings.mixed_integer_gomory_cuts != 0 || settings.strong_chvatal_gomory_cuts != 0) { if (toc(start_time) >= settings.time_limit) { return true; } f_t cut_start_time = tic(); - generate_gomory_cuts(lp, - settings, - Arow, - new_slacks, - var_types, - basis_update, - xstar, - basic_list, - nonbasic_list, - start_time); + if (!cut_pool_.pool_limit_reached() && + !((settings.strong_chvatal_gomory_cuts == 0 || + cut_pool_.generation_limit_reached(CHVATAL_GOMORY)) && + (settings.mixed_integer_gomory_cuts == 0 || + cut_pool_.generation_limit_reached(MIXED_INTEGER_GOMORY)))) { + generate_gomory_cuts(lp, + settings, + Arow, + new_slacks, + var_types, + basis_update, + xstar, + basic_list, + nonbasic_list, + start_time); + } f_t cut_generation_time = toc(cut_start_time); if (cut_generation_time > 1.0) { settings.log.debug("Gomory and CG cut generation time %.2f seconds\n", cut_generation_time); @@ -3534,7 +3553,9 @@ bool cut_generation_t::generate_cuts(const lp_problem_t& lp, if (settings.knapsack_cuts != 0) { if (toc(start_time) >= settings.time_limit) { return true; } f_t cut_start_time = tic(); - generate_knapsack_cuts(lp, settings, Arow, new_slacks, var_types, xstar, start_time); + if (!cut_pool_.generation_limit_reached(KNAPSACK)) { + generate_knapsack_cuts(lp, settings, Arow, new_slacks, var_types, xstar, start_time); + } f_t cut_generation_time = toc(cut_start_time); if (cut_generation_time > 1.0) { settings.log.debug("Knapsack cut generation time %.2f seconds\n", cut_generation_time); @@ -3545,7 +3566,9 @@ bool cut_generation_t::generate_cuts(const lp_problem_t& lp, if (settings.flow_cover_cuts != 0) { if (toc(start_time) >= settings.time_limit) { return true; } f_t cut_start_time = tic(); - generate_flow_cover_cuts(lp, settings, Arow, var_types, xstar, variable_bounds, start_time); + if (!cut_pool_.generation_limit_reached(FLOW_COVER)) { + generate_flow_cover_cuts(lp, settings, Arow, var_types, xstar, variable_bounds, start_time); + } f_t cut_generation_time = toc(cut_start_time); if (cut_generation_time > 1.0) { settings.log.debug("Flow cover cut generation time %.2f seconds\n", cut_generation_time); @@ -3556,8 +3579,14 @@ bool cut_generation_t::generate_cuts(const lp_problem_t& lp, if (settings.mir_cuts != 0 || settings.strong_chvatal_gomory_cuts != 0) { if (toc(start_time) >= settings.time_limit) { return true; } f_t cut_start_time = tic(); - generate_mir_cuts( - lp, settings, Arow, new_slacks, var_types, xstar, ystar, variable_bounds, start_time); + const bool mir_families_limited = + (settings.strong_chvatal_gomory_cuts == 0 || + cut_pool_.generation_limit_reached(CHVATAL_GOMORY)) && + (settings.mir_cuts == 0 || cut_pool_.generation_limit_reached(MIXED_INTEGER_ROUNDING)); + if (!cut_pool_.pool_limit_reached() && !mir_families_limited) { + generate_mir_cuts( + lp, settings, Arow, new_slacks, var_types, xstar, ystar, variable_bounds, start_time); + } f_t cut_generation_time = toc(cut_start_time); if (cut_generation_time > 1.0) { settings.log.debug("MIR and CG cut generation time %.2f seconds\n", cut_generation_time); @@ -3568,7 +3597,9 @@ bool cut_generation_t::generate_cuts(const lp_problem_t& lp, if (settings.implied_bound_cuts != 0) { if (toc(start_time) >= settings.time_limit) { return true; } f_t cut_start_time = tic(); - generate_implied_bound_cuts(lp, settings, var_types, xstar, start_time); + if (!cut_pool_.generation_limit_reached(IMPLIED_BOUND)) { + generate_implied_bound_cuts(lp, settings, var_types, xstar, start_time); + } f_t cut_generation_time = toc(cut_start_time); if (cut_generation_time > 1.0) { settings.log.debug("Implied bounds cut generation time %.2f seconds\n", cut_generation_time); @@ -3588,7 +3619,10 @@ bool cut_generation_t::generate_cuts(const lp_problem_t& lp, if (settings.clique_cuts != 0) { if (toc(start_time) >= settings.time_limit) { return true; } f_t cut_start_time = tic(); - bool feasible = generate_clique_cuts(lp, settings, var_types, xstar, zstar, start_time); + bool feasible = true; + if (!cut_pool_.generation_limit_reached(CLIQUE)) { + feasible = generate_clique_cuts(lp, settings, var_types, xstar, zstar, start_time); + } if (!feasible) { settings.log.printf("Clique cuts proved infeasible\n"); return false; @@ -3604,8 +3638,11 @@ 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, Arow, new_slacks, var_types, xstar, zstar, variable_bounds, start_time); + bool feasible = true; + if (!cut_pool_.generation_limit_reached(ZERO_HALF)) { + 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) { @@ -3635,7 +3672,9 @@ void cut_generation_t::generate_knapsack_cuts( { if (knapsack_generation_.num_knapsack_constraints() > 0) { for (i_t knapsack_row : knapsack_generation_.get_knapsack_constraints()) { - if (toc(start_time) >= settings.time_limit) { return; } + if (toc(start_time) >= settings.time_limit || cut_pool_.generation_limit_reached(KNAPSACK)) { + 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); @@ -3656,7 +3695,10 @@ void cut_generation_t::generate_flow_cover_cuts( { if (flow_cover_generation_.num_constraints() > 0) { for (const auto& flow_cover_row : flow_cover_generation_.get_constraints()) { - if (toc(start_time) >= settings.time_limit) { return; } + if (toc(start_time) >= settings.time_limit || + cut_pool_.generation_limit_reached(FLOW_COVER)) { + return; + } inequality_t cut(lp.num_cols); i_t status = flow_cover_generation_.generate_cut( lp, settings, Arow, variable_bounds, var_types, xstar, flow_cover_row, cut); @@ -3771,7 +3813,9 @@ bool cut_generation_t::generate_clique_cuts( size_t extension_gain = 0; #endif for (std::vector& clique_local : ctx.cliques) { - if (toc(start_time) >= settings.time_limit) { return true; } + if (toc(start_time) >= settings.time_limit || cut_pool_.generation_limit_reached(CLIQUE)) { + return true; + } #if DEBUG_CLIQUE_CUTS candidate_cliques++; #endif @@ -3952,6 +3996,7 @@ bool cut_generation_t::generate_zero_half_cuts( for (i_t s = 0; s < num_local; ++s) { if (toc(start_time) >= settings.time_limit) { break; } if (work_estimate > max_work_estimate) { break; } + if (cut_pool_.generation_limit_reached(ZERO_HALF)) { break; } if (already_used[s]) { continue; } ZERO_HALF_DEBUG("separation loop s=%lld / %lld", static_cast(s), @@ -4093,10 +4138,13 @@ void cut_generation_t::generate_mir_cuts( 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()) { + f_t work_estimate = 0.0; + while (!score_queue.empty()) { + const bool mir_families_limited = + (settings.strong_chvatal_gomory_cuts == 0 || + cut_pool_.generation_limit_reached(CHVATAL_GOMORY)) && + (settings.mir_cuts == 0 || cut_pool_.generation_limit_reached(MIXED_INTEGER_ROUNDING)); + if (cut_pool_.pool_limit_reached() || mir_families_limited) { break; } 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(); @@ -4356,7 +4404,13 @@ void cut_generation_t::generate_gomory_cuts( complemented_mir.bound_substitution(lp, variable_bounds, var_types, xstar, transformed_xstar); for (i_t i = 0; i < lp.num_rows; i++) { - if (toc(start_time) >= settings.time_limit) { break; } + if (toc(start_time) >= settings.time_limit || cut_pool_.pool_limit_reached() || + ((settings.strong_chvatal_gomory_cuts == 0 || + cut_pool_.generation_limit_reached(CHVATAL_GOMORY)) && + (settings.mixed_integer_gomory_cuts == 0 || + cut_pool_.generation_limit_reached(MIXED_INTEGER_GOMORY)))) { + break; + } inequality_t inequality(lp.num_cols); const i_t j = basic_list[i]; if (var_types[j] != variable_type_t::INTEGER) { continue; } @@ -5716,7 +5770,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 @@ -5724,6 +5781,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); @@ -5766,6 +5824,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) { @@ -5787,6 +5846,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 c052935806..f5f8c73532 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -306,6 +306,9 @@ std::vector> find_mod2_row_combinations_for_test( template class cut_pool_t { public: + static constexpr i_t max_pool_size = 150000; + static constexpr i_t max_cut_family_size = 20000; + cut_pool_t(i_t original_vars, const simplex::simplex_solver_settings_t& settings) : original_vars_(original_vars), settings_(settings), @@ -321,6 +324,15 @@ class cut_pool_t { // We expect that the cut is violated by the current relaxation xstar. void add_cut(cut_type_t cut_type, const inequality_t& cut); + bool generation_limit_reached(cut_type_t cut_type) const + { + return cut_storage_.m >= max_pool_size || cut_type_counts_[cut_type] >= max_cut_family_size; + } + + bool pool_limit_reached() const { return cut_storage_.m >= max_pool_size; } + + i_t cut_family_size(cut_type_t cut_type) const { return cut_type_counts_[cut_type]; } + void score_cuts(std::vector& x_relax); // We return the cuts in the form best_cuts*x <= best_rhs @@ -350,6 +362,7 @@ class cut_pool_t { std::vector rhs_storage_; std::vector cut_age_; std::vector cut_type_; + std::array cut_type_counts_{}; i_t scored_cuts_; std::vector cut_distances_; @@ -1054,7 +1067,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 index ac55131563..74f34ea540 100644 --- a/cpp/src/cuts/zero_half_mod2.cpp +++ b/cpp/src/cuts/zero_half_mod2.cpp @@ -315,15 +315,16 @@ void mod2_add_transformed_zero_half_cut( f_t& work_estimate, i_t& cuts_added) { - work_estimate += (f_t)(10 * transformed_cut.size() + 1); + 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); + 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); - ++cuts_added; + if (cut_pool.pool_size() > pool_size) { ++cuts_added; } } } @@ -477,17 +478,24 @@ bool generate_mod2_zero_half_cuts(cut_pool_t& cut_pool, 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; - const f_t max_work_estimate = work_estimate + (f_t)1e8; - bool work_limit_reached = false; + 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()), - &work_estimate, - max_work_estimate, - &work_limit_reached)) { + &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); @@ -503,26 +511,38 @@ bool generate_mod2_zero_half_cuts(cut_pool_t& cut_pool, transformed_xstar, start_time, settings.time_limit, - work_estimate, - max_work_estimate, - work_limit_reached); - - auto row_combinations = find_mod2_row_combinations( - candidates, max_combination_size, max_row_combinations, &work_estimate, max_work_estimate); - if (work_estimate > max_work_estimate) { work_limit_reached = true; } + candidate_work, + candidate_work_limit, + candidate_limit_reached); + + auto row_combinations = find_mod2_row_combinations(candidates, + max_combination_size, + max_row_combinations, + &combination_work, + combination_work_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 || work_limit_reached) { break; } + if (toc(start_time) >= settings.time_limit || generation_limit_reached || + cut_pool.generation_limit_reached(ZERO_HALF)) { + break; + } size_t aggregate_input_nz = 0; for (const i_t candidate_index : combination) { aggregate_input_nz += candidates[candidate_index].transformed_inequality.size(); } - const f_t aggregate_work = - (f_t)(4 * aggregate_input_nz + 1) + - (f_t)aggregate_input_nz * std::log2((f_t)aggregate_input_nz + (f_t)1.0); - if (add_work_estimate(aggregate_work, &work_estimate, max_work_estimate, &work_limit_reached)) { + if (add_work_estimate((f_t)(2 * aggregate_input_nz + 1), + &generation_work, + generation_work_limit, + &generation_limit_reached)) { break; } @@ -539,6 +559,15 @@ bool generate_mod2_zero_half_cuts(cut_pool_t& cut_pool, } 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); @@ -555,12 +584,12 @@ bool generate_mod2_zero_half_cuts(cut_pool_t& cut_pool, aggregate, min_violation, start_time, - work_estimate, - max_work_estimate, - work_limit_reached, + 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 && !work_limit_reached) { + if (reversible && toc(start_time) < settings.time_limit && !generation_limit_reached) { aggregate.negate(); mod2_generate_cuts_from_aggregate(complemented_mir, cut_pool, @@ -574,12 +603,13 @@ bool generate_mod2_zero_half_cuts(cut_pool_t& cut_pool, aggregate, min_violation, start_time, - work_estimate, - max_work_estimate, - work_limit_reached, + generation_work, + generation_work_limit, + generation_limit_reached, cuts_added); } } + work_estimate = candidate_work + combination_work + generation_work; return true; } @@ -671,6 +701,14 @@ bool complemented_mixed_integer_rounding_cut_t::generate_lifted_mixed_ } 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) { diff --git a/cpp/tests/mip/cuts_test.cu b/cpp/tests/mip/cuts_test.cu index 5af0754e3d..9c89c148f6 100644 --- a/cpp/tests/mip/cuts_test.cu +++ b/cpp/tests/mip/cuts_test.cu @@ -983,6 +983,45 @@ TEST(cuts, test_duplicate_cuts_detection) cut_pool.check_for_duplicate_cuts(); } +TEST(cuts, cut_pool_enforces_cut_family_limit) +{ + simplex::simplex_solver_settings_t settings; + mip::cut_pool_t cut_pool(1, settings); + mip::inequality_t cut; + cut.push_back(0, 1.0); + cut.rhs = 1.0; + + for (int i = 0; i <= cut_pool.max_cut_family_size; ++i) { + cut_pool.add_cut(mip::cut_type_t::MIXED_INTEGER_GOMORY, cut); + } + + EXPECT_EQ(cut_pool.pool_size(), cut_pool.max_cut_family_size); + EXPECT_EQ(cut_pool.cut_family_size(mip::cut_type_t::MIXED_INTEGER_GOMORY), + cut_pool.max_cut_family_size); + EXPECT_TRUE(cut_pool.generation_limit_reached(mip::cut_type_t::MIXED_INTEGER_GOMORY)); + + cut_pool.add_cut(mip::cut_type_t::MIXED_INTEGER_ROUNDING, cut); + EXPECT_EQ(cut_pool.pool_size(), cut_pool.max_cut_family_size + 1); +} + +TEST(cuts, cut_pool_enforces_total_limit) +{ + simplex::simplex_solver_settings_t settings; + mip::cut_pool_t cut_pool(1, settings); + mip::inequality_t cut; + cut.push_back(0, 1.0); + cut.rhs = 1.0; + + for (int i = 0; i < cut_pool.max_pool_size; ++i) { + const auto cut_type = (mip::cut_type_t)(i / cut_pool.max_cut_family_size); + cut_pool.add_cut(cut_type, cut); + } + cut_pool.add_cut(mip::cut_type_t::FLOW_COVER, cut); + + EXPECT_EQ(cut_pool.pool_size(), cut_pool.max_pool_size); + EXPECT_TRUE(cut_pool.pool_limit_reached()); +} + TEST(cuts, clique_phase1_smoke_conflict_graph_edges) { const raft::handle_t handle{}; From db89ce1cbbbb1914fe8e9e50fcb112605f35f411 Mon Sep 17 00:00:00 2001 From: akif Date: Wed, 12 Aug 2026 13:49:04 +0200 Subject: [PATCH 10/15] Bound long-running cut generation work Add cooperative time and work gates to knapsack lifting and mod-2 separation so root cuts cannot substantially overrun the solve deadline. Signed-off-by: akif --- cpp/src/cuts/cuts.cpp | 21 ++++++++++------ cpp/src/cuts/cuts.hpp | 9 ++++--- cpp/src/cuts/zero_half_mod2.cpp | 44 +++++++++++++++++++++------------ 3 files changed, 48 insertions(+), 26 deletions(-) diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index 97c5922c47..1ee37ffd4b 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -2327,7 +2327,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 @@ -2502,7 +2503,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: @@ -2682,7 +2684,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 @@ -2798,14 +2801,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; @@ -3024,8 +3028,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 @@ -3053,6 +3059,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); @@ -3677,7 +3684,7 @@ void cut_generation_t::generate_knapsack_cuts( } 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); } } } diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index f5f8c73532..095600672d 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -585,7 +585,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_; } @@ -610,7 +611,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, @@ -621,7 +623,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_; diff --git a/cpp/src/cuts/zero_half_mod2.cpp b/cpp/src/cuts/zero_half_mod2.cpp index 74f34ea540..c9927634cc 100644 --- a/cpp/src/cuts/zero_half_mod2.cpp +++ b/cpp/src/cuts/zero_half_mod2.cpp @@ -74,7 +74,9 @@ std::vector> find_mod2_row_combinations(const std::vector 0, "Maximum GF(2) combination size must be positive"); cuopt_assert(max_combinations > 0, "Maximum number of GF(2) combinations must be positive"); @@ -111,6 +113,7 @@ std::vector> find_mod2_row_combinations(const 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; @@ -165,12 +168,13 @@ i_t mod2_integral_scale(const inequality_t& inequality, f_t coefficient_integral_tol, f_t start_time, f_t time_limit, - f_t& work_estimate) + f_t& work_estimate, + f_t max_work_estimate, + bool& work_limit_reached) { - if (toc(start_time) >= time_limit) { return i_t{0}; } - f_t scale_work = 0.0; for (i_t scale = 1; scale <= max_integral_scale; ++scale) { - scale_work += 1.0; + 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)) > @@ -189,12 +193,11 @@ i_t mod2_integral_scale(const inequality_t& inequality, integral = false; } } - if (integral) { - work_estimate += scale_work; - return scale; + if (add_work_estimate(scale_work, &work_estimate, max_work_estimate, &work_limit_reached)) { + return i_t{0}; } + if (integral) { return scale; } } - work_estimate += scale_work; return i_t{0}; } @@ -273,11 +276,9 @@ std::vector> mod2_collect_candidates( coefficient_integral_tol, start_time, time_limit, - work_estimate); - if (work_estimate > max_work_estimate) { - work_limit_reached = true; - break; - } + 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); } @@ -347,7 +348,12 @@ void mod2_generate_cuts_from_aggregate( bool& work_limit_reached, i_t& cuts_added) { - work_estimate += (f_t)(3 * oriented_aggregate.size() + 1); + 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); @@ -515,11 +521,17 @@ bool generate_mod2_zero_half_cuts(cut_pool_t& cut_pool, 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); + combination_work_limit, + start_time, + settings.time_limit); if (add_work_estimate((f_t)(2 * lp.num_cols), &generation_work, generation_work_limit, From 7c2ed268bd407174db621d07c3ab1bb891cf99e5 Mon Sep 17 00:00:00 2001 From: akif Date: Wed, 12 Aug 2026 15:39:30 +0200 Subject: [PATCH 11/15] Remove cut pool size limits Rely on generator work and wall-time budgets instead of truncating accepted cuts by family or total pool size. Signed-off-by: akif --- cpp/src/cuts/cuts.cpp | 95 ++++++++------------------------- cpp/src/cuts/cuts.hpp | 13 ----- cpp/src/cuts/zero_half_mod2.cpp | 5 +- cpp/tests/mip/cuts_test.cu | 39 -------------- 4 files changed, 24 insertions(+), 128 deletions(-) diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index 1ee37ffd4b..ba20c00de9 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -1146,8 +1146,6 @@ inline uint64_t hash64_with_seed(uint64_t value, uint64_t seed) template void cut_pool_t::add_cut(cut_type_t cut_type, const inequality_t& cut) { - if (generation_limit_reached(cut_type)) { return; } - // TODO: Add fast duplicate check and only add if the cut is not already in the pool for (i_t p = 0; p < cut.size(); p++) { @@ -1170,7 +1168,6 @@ void cut_pool_t::add_cut(cut_type_t cut_type, const inequality_t @@ -3113,10 +3110,7 @@ void cut_generation_t::generate_implied_bound_cuts( const f_t max_work_estimate = 1e8; for (i_t j = 0; j < n_cols; j++) { - if (toc(start_time) >= settings.time_limit || - cut_pool_.generation_limit_reached(IMPLIED_BOUND)) { - return; - } + if (toc(start_time) >= settings.time_limit) { return; } if (var_types[j] == variable_type_t::CONTINUOUS) { continue; } const f_t xstar_j = xstar[j]; @@ -3534,22 +3528,16 @@ bool cut_generation_t::generate_cuts(const lp_problem_t& lp, if (settings.mixed_integer_gomory_cuts != 0 || settings.strong_chvatal_gomory_cuts != 0) { if (toc(start_time) >= settings.time_limit) { return true; } f_t cut_start_time = tic(); - if (!cut_pool_.pool_limit_reached() && - !((settings.strong_chvatal_gomory_cuts == 0 || - cut_pool_.generation_limit_reached(CHVATAL_GOMORY)) && - (settings.mixed_integer_gomory_cuts == 0 || - cut_pool_.generation_limit_reached(MIXED_INTEGER_GOMORY)))) { - generate_gomory_cuts(lp, - settings, - Arow, - new_slacks, - var_types, - basis_update, - xstar, - basic_list, - nonbasic_list, - start_time); - } + generate_gomory_cuts(lp, + settings, + Arow, + new_slacks, + var_types, + basis_update, + xstar, + basic_list, + nonbasic_list, + start_time); f_t cut_generation_time = toc(cut_start_time); if (cut_generation_time > 1.0) { settings.log.debug("Gomory and CG cut generation time %.2f seconds\n", cut_generation_time); @@ -3560,9 +3548,7 @@ bool cut_generation_t::generate_cuts(const lp_problem_t& lp, if (settings.knapsack_cuts != 0) { if (toc(start_time) >= settings.time_limit) { return true; } f_t cut_start_time = tic(); - if (!cut_pool_.generation_limit_reached(KNAPSACK)) { - generate_knapsack_cuts(lp, settings, Arow, new_slacks, var_types, xstar, start_time); - } + generate_knapsack_cuts(lp, settings, Arow, new_slacks, var_types, xstar, start_time); f_t cut_generation_time = toc(cut_start_time); if (cut_generation_time > 1.0) { settings.log.debug("Knapsack cut generation time %.2f seconds\n", cut_generation_time); @@ -3573,9 +3559,7 @@ bool cut_generation_t::generate_cuts(const lp_problem_t& lp, if (settings.flow_cover_cuts != 0) { if (toc(start_time) >= settings.time_limit) { return true; } f_t cut_start_time = tic(); - if (!cut_pool_.generation_limit_reached(FLOW_COVER)) { - generate_flow_cover_cuts(lp, settings, Arow, var_types, xstar, variable_bounds, start_time); - } + generate_flow_cover_cuts(lp, settings, Arow, var_types, xstar, variable_bounds, start_time); f_t cut_generation_time = toc(cut_start_time); if (cut_generation_time > 1.0) { settings.log.debug("Flow cover cut generation time %.2f seconds\n", cut_generation_time); @@ -3586,14 +3570,8 @@ bool cut_generation_t::generate_cuts(const lp_problem_t& lp, if (settings.mir_cuts != 0 || settings.strong_chvatal_gomory_cuts != 0) { if (toc(start_time) >= settings.time_limit) { return true; } f_t cut_start_time = tic(); - const bool mir_families_limited = - (settings.strong_chvatal_gomory_cuts == 0 || - cut_pool_.generation_limit_reached(CHVATAL_GOMORY)) && - (settings.mir_cuts == 0 || cut_pool_.generation_limit_reached(MIXED_INTEGER_ROUNDING)); - if (!cut_pool_.pool_limit_reached() && !mir_families_limited) { - generate_mir_cuts( - lp, settings, Arow, new_slacks, var_types, xstar, ystar, variable_bounds, start_time); - } + generate_mir_cuts( + lp, settings, Arow, new_slacks, var_types, xstar, ystar, variable_bounds, start_time); f_t cut_generation_time = toc(cut_start_time); if (cut_generation_time > 1.0) { settings.log.debug("MIR and CG cut generation time %.2f seconds\n", cut_generation_time); @@ -3604,9 +3582,7 @@ bool cut_generation_t::generate_cuts(const lp_problem_t& lp, if (settings.implied_bound_cuts != 0) { if (toc(start_time) >= settings.time_limit) { return true; } f_t cut_start_time = tic(); - if (!cut_pool_.generation_limit_reached(IMPLIED_BOUND)) { - generate_implied_bound_cuts(lp, settings, var_types, xstar, start_time); - } + generate_implied_bound_cuts(lp, settings, var_types, xstar, start_time); f_t cut_generation_time = toc(cut_start_time); if (cut_generation_time > 1.0) { settings.log.debug("Implied bounds cut generation time %.2f seconds\n", cut_generation_time); @@ -3626,10 +3602,7 @@ bool cut_generation_t::generate_cuts(const lp_problem_t& lp, if (settings.clique_cuts != 0) { if (toc(start_time) >= settings.time_limit) { return true; } f_t cut_start_time = tic(); - bool feasible = true; - if (!cut_pool_.generation_limit_reached(CLIQUE)) { - feasible = generate_clique_cuts(lp, settings, var_types, xstar, zstar, start_time); - } + bool feasible = generate_clique_cuts(lp, settings, var_types, xstar, zstar, start_time); if (!feasible) { settings.log.printf("Clique cuts proved infeasible\n"); return false; @@ -3645,11 +3618,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 = true; - if (!cut_pool_.generation_limit_reached(ZERO_HALF)) { - feasible = generate_zero_half_cuts( - lp, settings, Arow, new_slacks, var_types, xstar, zstar, variable_bounds, 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) { @@ -3679,9 +3649,7 @@ void cut_generation_t::generate_knapsack_cuts( { if (knapsack_generation_.num_knapsack_constraints() > 0) { for (i_t knapsack_row : knapsack_generation_.get_knapsack_constraints()) { - if (toc(start_time) >= settings.time_limit || cut_pool_.generation_limit_reached(KNAPSACK)) { - return; - } + 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, start_time); @@ -3702,10 +3670,7 @@ void cut_generation_t::generate_flow_cover_cuts( { if (flow_cover_generation_.num_constraints() > 0) { for (const auto& flow_cover_row : flow_cover_generation_.get_constraints()) { - if (toc(start_time) >= settings.time_limit || - cut_pool_.generation_limit_reached(FLOW_COVER)) { - return; - } + if (toc(start_time) >= settings.time_limit) { return; } inequality_t cut(lp.num_cols); i_t status = flow_cover_generation_.generate_cut( lp, settings, Arow, variable_bounds, var_types, xstar, flow_cover_row, cut); @@ -3820,9 +3785,7 @@ bool cut_generation_t::generate_clique_cuts( size_t extension_gain = 0; #endif for (std::vector& clique_local : ctx.cliques) { - if (toc(start_time) >= settings.time_limit || cut_pool_.generation_limit_reached(CLIQUE)) { - return true; - } + if (toc(start_time) >= settings.time_limit) { return true; } #if DEBUG_CLIQUE_CUTS candidate_cliques++; #endif @@ -4003,7 +3966,6 @@ bool cut_generation_t::generate_zero_half_cuts( for (i_t s = 0; s < num_local; ++s) { if (toc(start_time) >= settings.time_limit) { break; } if (work_estimate > max_work_estimate) { break; } - if (cut_pool_.generation_limit_reached(ZERO_HALF)) { break; } if (already_used[s]) { continue; } ZERO_HALF_DEBUG("separation loop s=%lld / %lld", static_cast(s), @@ -4147,11 +4109,6 @@ void cut_generation_t::generate_mir_cuts( f_t work_estimate = 0.0; while (!score_queue.empty()) { - const bool mir_families_limited = - (settings.strong_chvatal_gomory_cuts == 0 || - cut_pool_.generation_limit_reached(CHVATAL_GOMORY)) && - (settings.mir_cuts == 0 || cut_pool_.generation_limit_reached(MIXED_INTEGER_ROUNDING)); - if (cut_pool_.pool_limit_reached() || mir_families_limited) { break; } 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(); @@ -4411,13 +4368,7 @@ void cut_generation_t::generate_gomory_cuts( complemented_mir.bound_substitution(lp, variable_bounds, var_types, xstar, transformed_xstar); for (i_t i = 0; i < lp.num_rows; i++) { - if (toc(start_time) >= settings.time_limit || cut_pool_.pool_limit_reached() || - ((settings.strong_chvatal_gomory_cuts == 0 || - cut_pool_.generation_limit_reached(CHVATAL_GOMORY)) && - (settings.mixed_integer_gomory_cuts == 0 || - cut_pool_.generation_limit_reached(MIXED_INTEGER_GOMORY)))) { - break; - } + if (toc(start_time) >= settings.time_limit) { break; } inequality_t inequality(lp.num_cols); const i_t j = basic_list[i]; if (var_types[j] != variable_type_t::INTEGER) { continue; } diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index 095600672d..fcb6080178 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -306,9 +306,6 @@ std::vector> find_mod2_row_combinations_for_test( template class cut_pool_t { public: - static constexpr i_t max_pool_size = 150000; - static constexpr i_t max_cut_family_size = 20000; - cut_pool_t(i_t original_vars, const simplex::simplex_solver_settings_t& settings) : original_vars_(original_vars), settings_(settings), @@ -324,15 +321,6 @@ class cut_pool_t { // We expect that the cut is violated by the current relaxation xstar. void add_cut(cut_type_t cut_type, const inequality_t& cut); - bool generation_limit_reached(cut_type_t cut_type) const - { - return cut_storage_.m >= max_pool_size || cut_type_counts_[cut_type] >= max_cut_family_size; - } - - bool pool_limit_reached() const { return cut_storage_.m >= max_pool_size; } - - i_t cut_family_size(cut_type_t cut_type) const { return cut_type_counts_[cut_type]; } - void score_cuts(std::vector& x_relax); // We return the cuts in the form best_cuts*x <= best_rhs @@ -362,7 +350,6 @@ class cut_pool_t { std::vector rhs_storage_; std::vector cut_age_; std::vector cut_type_; - std::array cut_type_counts_{}; i_t scored_cuts_; std::vector cut_distances_; diff --git a/cpp/src/cuts/zero_half_mod2.cpp b/cpp/src/cuts/zero_half_mod2.cpp index c9927634cc..13e2413fa1 100644 --- a/cpp/src/cuts/zero_half_mod2.cpp +++ b/cpp/src/cuts/zero_half_mod2.cpp @@ -542,10 +542,7 @@ bool generate_mod2_zero_half_cuts(cut_pool_t& cut_pool, scratch_pad_t aggregate_pad(lp.num_cols); for (const auto& combination : row_combinations) { - if (toc(start_time) >= settings.time_limit || generation_limit_reached || - cut_pool.generation_limit_reached(ZERO_HALF)) { - break; - } + if (toc(start_time) >= settings.time_limit || generation_limit_reached) { break; } size_t aggregate_input_nz = 0; for (const i_t candidate_index : combination) { diff --git a/cpp/tests/mip/cuts_test.cu b/cpp/tests/mip/cuts_test.cu index 9c89c148f6..5af0754e3d 100644 --- a/cpp/tests/mip/cuts_test.cu +++ b/cpp/tests/mip/cuts_test.cu @@ -983,45 +983,6 @@ TEST(cuts, test_duplicate_cuts_detection) cut_pool.check_for_duplicate_cuts(); } -TEST(cuts, cut_pool_enforces_cut_family_limit) -{ - simplex::simplex_solver_settings_t settings; - mip::cut_pool_t cut_pool(1, settings); - mip::inequality_t cut; - cut.push_back(0, 1.0); - cut.rhs = 1.0; - - for (int i = 0; i <= cut_pool.max_cut_family_size; ++i) { - cut_pool.add_cut(mip::cut_type_t::MIXED_INTEGER_GOMORY, cut); - } - - EXPECT_EQ(cut_pool.pool_size(), cut_pool.max_cut_family_size); - EXPECT_EQ(cut_pool.cut_family_size(mip::cut_type_t::MIXED_INTEGER_GOMORY), - cut_pool.max_cut_family_size); - EXPECT_TRUE(cut_pool.generation_limit_reached(mip::cut_type_t::MIXED_INTEGER_GOMORY)); - - cut_pool.add_cut(mip::cut_type_t::MIXED_INTEGER_ROUNDING, cut); - EXPECT_EQ(cut_pool.pool_size(), cut_pool.max_cut_family_size + 1); -} - -TEST(cuts, cut_pool_enforces_total_limit) -{ - simplex::simplex_solver_settings_t settings; - mip::cut_pool_t cut_pool(1, settings); - mip::inequality_t cut; - cut.push_back(0, 1.0); - cut.rhs = 1.0; - - for (int i = 0; i < cut_pool.max_pool_size; ++i) { - const auto cut_type = (mip::cut_type_t)(i / cut_pool.max_cut_family_size); - cut_pool.add_cut(cut_type, cut); - } - cut_pool.add_cut(mip::cut_type_t::FLOW_COVER, cut); - - EXPECT_EQ(cut_pool.pool_size(), cut_pool.max_pool_size); - EXPECT_TRUE(cut_pool.pool_limit_reached()); -} - TEST(cuts, clique_phase1_smoke_conflict_graph_edges) { const raft::handle_t handle{}; From ee7ce6b963dc5399dc0331b5b58f16e56b246fa5 Mon Sep 17 00:00:00 2001 From: akif Date: Wed, 12 Aug 2026 15:41:53 +0200 Subject: [PATCH 12/15] Reduce mod-2 zero-half work budget Limit candidate, dependency, and generation phases to a combined 10 million work units to reduce root cut cost. Signed-off-by: akif --- cpp/src/cuts/zero_half_mod2.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cpp/src/cuts/zero_half_mod2.cpp b/cpp/src/cuts/zero_half_mod2.cpp index 13e2413fa1..8a3b4dea23 100644 --- a/cpp/src/cuts/zero_half_mod2.cpp +++ b/cpp/src/cuts/zero_half_mod2.cpp @@ -487,9 +487,9 @@ bool generate_mod2_zero_half_cuts(cut_pool_t& cut_pool, 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; + constexpr f_t candidate_work_limit = (f_t)3e6; + constexpr f_t combination_work_limit = (f_t)3e6; + constexpr f_t generation_work_limit = (f_t)4e6; f_t candidate_work = 0.0; f_t combination_work = 0.0; f_t generation_work = 0.0; From 9772b7531e5ca9fbda15a58543f6f2798ffbaa70 Mon Sep 17 00:00:00 2001 From: akif Date: Thu, 13 Aug 2026 15:28:50 +0200 Subject: [PATCH 13/15] Revert \"Reduce mod-2 zero-half work budget\" Restore the previous zero-half work budget behavior. --- cpp/src/cuts/zero_half_mod2.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cpp/src/cuts/zero_half_mod2.cpp b/cpp/src/cuts/zero_half_mod2.cpp index 8a3b4dea23..13e2413fa1 100644 --- a/cpp/src/cuts/zero_half_mod2.cpp +++ b/cpp/src/cuts/zero_half_mod2.cpp @@ -487,9 +487,9 @@ bool generate_mod2_zero_half_cuts(cut_pool_t& cut_pool, 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)3e6; - constexpr f_t combination_work_limit = (f_t)3e6; - constexpr f_t generation_work_limit = (f_t)4e6; + 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; From 4e373e265a565e612b7787023667d510ef42ffb4 Mon Sep 17 00:00:00 2001 From: akif Date: Fri, 14 Aug 2026 16:32:23 +0200 Subject: [PATCH 14/15] cleanups and revert some changes --- cpp/src/branch_and_bound/branch_and_bound.cpp | 7 +-- cpp/src/branch_and_bound/branch_and_bound.hpp | 1 - cpp/src/cuts/cuts.cpp | 60 +++++++++++-------- 3 files changed, 35 insertions(+), 33 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 462ba43b0f..31bbe256bb 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -2937,7 +2937,6 @@ auto branch_and_bound_t::do_cut_pass( f_t root_relax_objective, i_t& cut_pool_size, f_t& cut_scoring_time, - i_t& max_scoring_pool_size, [[maybe_unused]] const std::vector& saved_solution) -> cut_pass_result_t { #ifdef PRINT_FRACTIONAL_INFO @@ -2981,8 +2980,7 @@ auto branch_and_bound_t::do_cut_pass( settings_.log.debug("Cut generation time %.2f seconds\n", cut_generation_time); } // Score the cuts - max_scoring_pool_size = std::max(max_scoring_pool_size, cut_pool.pool_size()); - f_t score_start_time = tic(); + f_t score_start_time = tic(); cut_pool.score_cuts(root_relax_soln_.x); cut_scoring_time += toc(score_start_time); // Get the best cuts from the cut pool @@ -3465,7 +3463,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut f_t cut_generation_start_time = tic(); i_t cut_pool_size = 0; f_t cut_scoring_time = 0.0; - i_t max_scoring_pool_size = 0; for (i_t cut_pass = 0; cut_pass < settings_.max_cut_passes; cut_pass++) { if (toc(exploration_stats_.start_time) >= settings_.time_limit) { solver_status_ = mip_status_t::TIME_LIMIT; @@ -3515,7 +3512,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut root_relax_objective, cut_pool_size, cut_scoring_time, - max_scoring_pool_size, saved_solution); root_fj_cpu_worker.stop(); @@ -3567,7 +3563,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut settings_.log.printf("Cut generation time: %.2f seconds\n", cut_generation_time); settings_.log.printf("Cut scoring time : %.2f seconds\n", cut_scoring_time); - settings_.log.printf("Cut scoring max pool: %d\n", max_scoring_pool_size); settings_.log.printf("Cut pool size : %d\n", cut_pool_size); settings_.log.printf("Size with cuts : %d constraints, %d variables, %d nonzeros\n", original_lp_.num_rows, diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index 6612090952..a56d382db2 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -324,7 +324,6 @@ class branch_and_bound_t { f_t root_relax_objective, i_t& cut_pool_size, f_t& cut_scoring_time, - i_t& max_scoring_pool_size, const std::vector& saved_solution); // Set the solution when found at the root node diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index ba20c00de9..3fa79f88c7 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -3102,15 +3102,17 @@ 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); - f_t work_estimate = 0.0; - const f_t max_work_estimate = 1e8; + 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 (toc(start_time) >= settings.time_limit) { return; } + 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]; @@ -3119,11 +3121,8 @@ void cut_generation_t::generate_implied_bound_cuts( 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]; - if (add_work_estimate( - (f_t)(zero_end - zero_begin + one_end - one_begin), &work_estimate, max_work_estimate)) { - return; - } 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]; @@ -3143,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++; } } @@ -3161,13 +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 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]; @@ -3187,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++; } } @@ -3204,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) { @@ -3900,19 +3906,18 @@ bool cut_generation_t::generate_zero_half_cuts( static_cast(sub_cg_.ready), sub_cg_.vertices.size()); - f_t mod2_work_estimate = 0.0; - if (!generate_mod2_zero_half_cuts(cut_pool_, - lp, - settings, - Arow, - new_slacks, - var_types, - xstar, - variable_bounds, - start_time, - mod2_work_estimate)) { - return true; - } + 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() and remains a complementary @@ -4102,19 +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); - f_t work_estimate = 0.0; - while (!score_queue.empty()) { + f_t work_estimate = 0.0; + 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; From cb33afff1ac41cea388e966159a27bf4ecda8da8 Mon Sep 17 00:00:00 2001 From: akif Date: Fri, 14 Aug 2026 17:56:46 +0200 Subject: [PATCH 15/15] restore cut scoring timer --- cpp/src/branch_and_bound/branch_and_bound.cpp | 9 +++------ cpp/src/branch_and_bound/branch_and_bound.hpp | 1 - 2 files changed, 3 insertions(+), 7 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 31bbe256bb..6fdf374c2e 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -2936,7 +2936,6 @@ auto branch_and_bound_t::do_cut_pass( f_t& last_objective, f_t root_relax_objective, i_t& cut_pool_size, - f_t& cut_scoring_time, [[maybe_unused]] const std::vector& saved_solution) -> cut_pass_result_t { #ifdef PRINT_FRACTIONAL_INFO @@ -2982,7 +2981,8 @@ auto branch_and_bound_t::do_cut_pass( // Score the cuts f_t score_start_time = tic(); cut_pool.score_cuts(root_relax_soln_.x); - cut_scoring_time += toc(score_start_time); + f_t score_time = toc(score_start_time); + if (score_time > 1.0) { settings_.log.debug("Cut scoring time %.2f seconds\n", score_time); } // Get the best cuts from the cut pool csr_matrix_t cuts_to_add(0, original_lp_.num_cols, 0); std::vector cut_rhs; @@ -3462,7 +3462,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut f_t cut_generation_start_time = tic(); i_t cut_pool_size = 0; - f_t cut_scoring_time = 0.0; for (i_t cut_pass = 0; cut_pass < settings_.max_cut_passes; cut_pass++) { if (toc(exploration_stats_.start_time) >= settings_.time_limit) { solver_status_ = mip_status_t::TIME_LIMIT; @@ -3511,7 +3510,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut last_objective, root_relax_objective, cut_pool_size, - cut_scoring_time, saved_solution); root_fj_cpu_worker.stop(); @@ -3562,8 +3560,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut mutex_upper_.unlock(); settings_.log.printf("Cut generation time: %.2f seconds\n", cut_generation_time); - settings_.log.printf("Cut scoring time : %.2f seconds\n", cut_scoring_time); - settings_.log.printf("Cut pool size : %d\n", cut_pool_size); + settings_.log.printf("Cut pool size : %d\n", cut_pool_size); settings_.log.printf("Size with cuts : %d constraints, %d variables, %d nonzeros\n", original_lp_.num_rows, original_lp_.num_cols, diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index a56d382db2..96b8a6d8fe 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -323,7 +323,6 @@ class branch_and_bound_t { f_t& last_objective, f_t root_relax_objective, i_t& cut_pool_size, - f_t& cut_scoring_time, const std::vector& saved_solution); // Set the solution when found at the root node