diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index 7acd7dee0a..0bb4f3733c 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -1485,6 +1485,18 @@ bool flow_cover_is_zero_one_integer_variable(const flow_cover_context_t +static bool flow_cover_valid_endpoint_terms(f_t endpoint_term, + f_t binary_coefficient, + bool in_n2, + f_t bound_tol) +{ + return in_n2 ? endpoint_term <= bound_tol && + endpoint_term + binary_coefficient <= bound_tol + : endpoint_term >= -bound_tol && + endpoint_term + binary_coefficient >= -bound_tol; +} + // Per-arc feasibility tolerances shared by the arc-acceptance gate and the assertion in // build_single_node_flow_relaxation, so the two sites cannot drift apart. template @@ -1685,6 +1697,78 @@ flow_cover_generation_t::flow_cover_generation_t( } } +// Repeated (variable, side, coefficient) queries share the same a=0 candidate because xstar and +// implied bounds are fixed for the pass. Cache that candidate and index nonzero-a candidates by +// controller so each row examines only bounds affected by its binary coefficients, to avoid an +// O(N^2) occurrence-by-bound scan when both grow with the model. +template +void flow_cover_generation_t::preprocess_cut_pass( + const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + const variable_bounds_t& variable_bounds, + const std::vector& var_types, + const std::vector& xstar) +{ + cut_pass_preprocessed = false; + cuopt_assert(var_types.size() >= (size_t)lp.num_cols, ""); + cuopt_assert(xstar.size() >= (size_t)lp.num_cols, ""); + for (flow_cover_bound_side_t side : + {flow_cover_bound_side_t::UPPER, flow_cover_bound_side_t::LOWER}) { + const bool use_upper_bound = side == flow_cover_bound_side_t::UPPER; + const auto& offsets = use_upper_bound ? variable_bounds.upper_offsets + : variable_bounds.lower_offsets; + const auto& variables = use_upper_bound ? variable_bounds.upper_variables + : variable_bounds.lower_variables; + const auto& weights = use_upper_bound ? variable_bounds.upper_weights + : variable_bounds.lower_weights; + const auto& biases = use_upper_bound ? variable_bounds.upper_biases + : variable_bounds.lower_biases; + auto& preprocessed = use_upper_bound ? upper_implied_bounds : lower_implied_bounds; + cuopt_assert(offsets.size() >= (size_t)(lp.num_cols + 1), ""); + cuopt_assert(variables.size() == weights.size(), ""); + cuopt_assert(variables.size() == biases.size(), ""); + + preprocessed.bounds.resize(variables.size()); + preprocessed.by_controller.clear(); + preprocessed.by_controller.resize(lp.num_cols); + preprocessed.zero_candidate_cache.clear(); + preprocessed.zero_candidate_cache.resize(lp.num_cols); + + const f_t bound_tol = settings.primal_tol; + for (i_t j = 0; j < lp.num_cols; j++) { + cuopt_assert(offsets[j] <= offsets[j + 1], ""); + cuopt_assert(offsets[j] >= 0, ""); + cuopt_assert(offsets[j + 1] <= (i_t)variables.size(), ""); + for (i_t p = offsets[j]; p < offsets[j + 1]; p++) { + const i_t x_col = variables[p]; + cuopt_assert(x_col >= 0 && x_col < lp.num_cols, ""); + const bool zero_one_eligible = + var_types[x_col] == variable_type_t::INTEGER && + std::abs(lp.lower[x_col]) <= bound_tol && + std::abs(lp.upper[x_col] - 1.0) <= bound_tol; + const bool finite_bound = std::isfinite(weights[p]) && std::isfinite(biases[p]); + + auto& bound = preprocessed.bounds[p]; + bound.zero_one_eligible = zero_one_eligible && finite_bound; + bound.active_bound = finite_bound ? weights[p] * xstar[x_col] + biases[p] : inf; + bound.distance = finite_bound ? std::abs(bound.active_bound - xstar[j]) : inf; + if (zero_one_eligible && finite_bound) { + auto& group = preprocessed.by_controller[j][x_col]; + if (group.bounds.empty()) { + group.minimum_alpha = biases[p]; + group.maximum_alpha = biases[p]; + } else { + group.minimum_alpha = std::min(group.minimum_alpha, biases[p]); + group.maximum_alpha = std::max(group.maximum_alpha, biases[p]); + } + group.bounds.push_back(p); + } + } + } + } + cut_pass_preprocessed = true; +} + template bool flow_cover_generation_t::normalize_row_side( const flow_cover_context_t& context, @@ -1772,6 +1856,56 @@ bool flow_cover_generation_t::normalize_row_side( return true; } +template +bool flow_cover_generation_t::try_add_implied_bound_candidate( + const flow_cover_context_t& context, + i_t variable, + f_t coefficient, + bool use_upper_bound, + i_t bound, + f_t binary_coefficient) +{ + const auto& preprocessed = + use_upper_bound ? upper_implied_bounds : lower_implied_bounds; + const auto& bound_variables = use_upper_bound ? context.variable_bounds.upper_variables + : context.variable_bounds.lower_variables; + const auto& bound_weights = use_upper_bound ? context.variable_bounds.upper_weights + : context.variable_bounds.lower_weights; + const auto& bound_biases = use_upper_bound ? context.variable_bounds.upper_biases + : context.variable_bounds.lower_biases; + cuopt_assert(preprocessed.bounds[bound].zero_one_eligible, ""); + + const f_t endpoint = + use_upper_bound ? context.lp.lower[variable] : context.lp.upper[variable]; + const bool in_n2 = use_upper_bound ? coefficient < 0.0 : coefficient > 0.0; + const i_t x_col = bound_variables[bound]; + const f_t gamma = bound_weights[bound]; + const f_t alpha = bound_biases[bound]; + const f_t signed_capacity = coefficient * gamma + binary_coefficient; + const f_t endpoint_term = coefficient * (endpoint - alpha); + const f_t bound_tol = context.settings.primal_tol; + const bool valid_endpoint = + flow_cover_valid_endpoint_terms(endpoint_term, binary_coefficient, in_n2, bound_tol) && + (in_n2 ? signed_capacity <= bound_tol : signed_capacity >= -bound_tol); + if (!valid_endpoint) { return false; } + + flow_cover_arc_spec_t spec; + spec.u = in_n2 ? -signed_capacity : signed_capacity; + spec.in_n2 = in_n2; + spec.x_col = x_col; + spec.fixed_x = 0.0; + spec.y_const = in_n2 ? coefficient * alpha : -coefficient * alpha; + spec.y_col = variable; + spec.y_coeff = in_n2 ? -coefficient : coefficient; + spec.y_x_coeff = in_n2 ? -binary_coefficient : binary_coefficient; + spec.b_shift = coefficient * alpha; + spec.active_bound = preprocessed.bounds[bound].active_bound; + spec.absorbs_binary_coeff = std::abs(binary_coefficient) > static_cast(1e-6); + const size_t size = candidates.size(); + flow_cover_try_add_candidate(context, spec, context.xstar[variable], candidates); + return candidates.size() > size; +} + template bool flow_cover_generation_t::build_single_node_flow_relaxation( const flow_cover_context_t& context, f_t b, f_t& single_node_flow_b) @@ -1779,7 +1913,6 @@ bool flow_cover_generation_t::build_single_node_flow_relaxation( auto& scratch = *this; const f_t coefficient_tol = static_cast(1e-6); const f_t feasibility_tol = context.settings.primal_tol; - const f_t bound_tol = context.settings.primal_tol; f_t b_shift = 0.0; scratch.arcs.reserve(scratch.continuous_terms.size() + scratch.binary_columns.size()); @@ -1791,52 +1924,83 @@ bool flow_cover_generation_t::build_single_node_flow_relaxation( if (use_upper_bound && lower_j <= -inf) { return; } if (!use_upper_bound && upper_j >= inf) { return; } - const i_t start = use_upper_bound ? context.variable_bounds.upper_offsets[j] - : context.variable_bounds.lower_offsets[j]; - const i_t end = use_upper_bound ? context.variable_bounds.upper_offsets[j + 1] - : context.variable_bounds.lower_offsets[j + 1]; - const f_t endpoint = use_upper_bound ? lower_j : upper_j; - - for (i_t p = start; p < end; p++) { - const i_t x_col = use_upper_bound ? context.variable_bounds.upper_variables[p] - : context.variable_bounds.lower_variables[p]; - if (!flow_cover_is_zero_one_integer_variable(context, x_col)) { continue; } - const f_t gamma = use_upper_bound ? context.variable_bounds.upper_weights[p] - : context.variable_bounds.lower_weights[p]; - const f_t alpha = use_upper_bound ? context.variable_bounds.upper_biases[p] - : context.variable_bounds.lower_biases[p]; - if (!std::isfinite(gamma) || !std::isfinite(alpha)) { continue; } - - const f_t direct_coeff = scratch.binary_coefficients_touched[x_col] - ? scratch.binary_coefficients[x_col] - : static_cast(0.0); - const std::array a_values = {direct_coeff, 0.0}; - const i_t num_a_values = std::abs(direct_coeff) > coefficient_tol ? 2 : 1; - for (i_t h = 0; h < num_a_values; h++) { - const f_t a = a_values[h]; - const bool in_n2 = use_upper_bound ? c < 0.0 : c > 0.0; - const f_t signed_capacity = c * gamma + a; - const f_t endpoint_term = c * (endpoint - alpha); - const bool valid_endpoint = - in_n2 ? (endpoint_term <= bound_tol && endpoint_term + a <= bound_tol && - signed_capacity <= bound_tol) - : (endpoint_term >= -bound_tol && endpoint_term + a >= -bound_tol && - signed_capacity >= -bound_tol); - if (!valid_endpoint) { continue; } + auto& preprocessed = use_upper_bound ? upper_implied_bounds : lower_implied_bounds; + const auto& source_offsets = use_upper_bound ? context.variable_bounds.upper_offsets + : context.variable_bounds.lower_offsets; + const auto& bound_variables = use_upper_bound ? context.variable_bounds.upper_variables + : context.variable_bounds.lower_variables; + const i_t source_start = source_offsets[j]; + const i_t source_end = source_offsets[j + 1]; + + // Small nonzero coefficients exclude a=0 for their controller, making the cache row-dependent. + bool has_small_direct_coeff = false; + for (i_t x_col : scratch.binary_columns) { + const f_t direct_coeff = scratch.binary_coefficients[x_col]; + if (direct_coeff != 0.0 && std::abs(direct_coeff) <= coefficient_tol) { + has_small_direct_coeff = true; + break; + } + } - flow_cover_arc_spec_t spec; - spec.u = in_n2 ? -signed_capacity : signed_capacity; - spec.in_n2 = in_n2; - spec.x_col = x_col; - spec.fixed_x = 0.0; - spec.y_const = in_n2 ? c * alpha : -c * alpha; - spec.y_col = j; - spec.y_coeff = in_n2 ? -c : c; - spec.y_x_coeff = in_n2 ? -a : a; - spec.b_shift = c * alpha; - spec.active_bound = gamma * context.xstar[x_col] + alpha; - spec.absorbs_binary_coeff = std::abs(a) > coefficient_tol; - flow_cover_try_add_candidate(context, spec, context.xstar[j], scratch.candidates); + auto& zero_candidate_cache = preprocessed.zero_candidate_cache[j]; + auto cached_zero = zero_candidate_cache.end(); + if (!has_small_direct_coeff) { cached_zero = zero_candidate_cache.find(c); } + if (has_small_direct_coeff || cached_zero == zero_candidate_cache.end()) { + // Compute the best a=0 candidate when it is uncached or row-specific. + i_t best_zero = -1; + f_t best_distance = inf; + for (i_t p = source_start; p < source_end; p++) { + if (!preprocessed.bounds[p].zero_one_eligible) { continue; } + const i_t x_col = bound_variables[p]; + const f_t direct_coeff = scratch.binary_coefficients_touched[x_col] + ? scratch.binary_coefficients[x_col] + : 0.0; + if (direct_coeff != 0.0 && std::abs(direct_coeff) <= coefficient_tol) { continue; } + // The normal acceptance gate keeps cached and uncached candidate validity identical. + if (!try_add_implied_bound_candidate(context, j, c, use_upper_bound, p, 0.0)) { + continue; + } + scratch.candidates.pop_back(); + if (best_zero < 0 || preprocessed.bounds[p].distance < best_distance) { + best_zero = p; + best_distance = preprocessed.bounds[p].distance; + } + } + if (!has_small_direct_coeff) { zero_candidate_cache.emplace(c, best_zero); } + if (best_zero >= 0) { + const bool added = + try_add_implied_bound_candidate(context, j, c, use_upper_bound, best_zero, 0.0); + cuopt_assert(added, ""); + } + } else { + // Reuse the pass-wide a=0 candidate for this variable, side, and coefficient. + if (cached_zero->second >= 0) { + const i_t p = cached_zero->second; + const bool added = + try_add_implied_bound_candidate(context, j, c, use_upper_bound, p, 0.0); + cuopt_assert(added, ""); + } + } + + for (i_t x_col : scratch.binary_columns) { + const f_t direct_coeff = scratch.binary_coefficients[x_col]; + if (direct_coeff == 0.0) { continue; } + const auto controller = preprocessed.by_controller[j].find(x_col); + if (controller == preprocessed.by_controller[j].end()) { continue; } + const auto& group = controller->second; + cuopt_assert(!group.bounds.empty(), ""); + cuopt_assert(group.minimum_alpha <= group.maximum_alpha, ""); + // Endpoint feasibility is monotone in alpha, so the side's extreme bounds the whole group. + const f_t endpoint = use_upper_bound ? lower_j : upper_j; + const f_t alpha = use_upper_bound ? group.minimum_alpha : group.maximum_alpha; + const f_t endpoint_term = c * (endpoint - alpha); + const bool in_n2 = use_upper_bound ? c < 0.0 : c > 0.0; + if (!flow_cover_valid_endpoint_terms( + endpoint_term, direct_coeff, in_n2, context.settings.primal_tol)) { + continue; + } + for (i_t p : group.bounds) { + try_add_implied_bound_candidate(context, j, c, use_upper_bound, p, direct_coeff); } } }; @@ -2280,6 +2444,7 @@ i_t flow_cover_generation_t::generate_cut( const flow_cover_row_t& flow_cover_row, inequality_t& cut) { + cuopt_assert(cut_pass_preprocessed, ""); flow_cover_context_t context{lp, settings, Arow, variable_bounds, var_types, xstar}; clear_cut_state(lp.num_cols); @@ -3675,6 +3840,10 @@ void cut_generation_t::generate_flow_cover_cuts( f_t start_time) { if (flow_cover_generation_.num_constraints() > 0) { + if (toc(start_time) >= settings.time_limit) { return; } + flow_cover_generation_.preprocess_cut_pass( + lp, settings, variable_bounds, var_types, xstar); + if (toc(start_time) >= settings.time_limit) { return; } for (const auto& flow_cover_row : flow_cover_generation_.get_constraints()) { if (toc(start_time) >= settings.time_limit) { return; } inequality_t cut(lp.num_cols); diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index 17938c964f..416df9f59f 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -441,6 +441,12 @@ class flow_cover_generation_t { csr_matrix_t& Arow, const std::vector& new_slacks); + void preprocess_cut_pass(const simplex::lp_problem_t& lp, + const simplex::simplex_solver_settings_t& settings, + const variable_bounds_t& variable_bounds, + const std::vector& var_types, + const std::vector& xstar); + i_t generate_cut(const simplex::lp_problem_t& lp, const simplex::simplex_solver_settings_t& settings, csr_matrix_t& Arow, @@ -486,6 +492,33 @@ class flow_cover_generation_t { const flow_cover_evaluation_t& simple_generalized_inequality, inequality_t& cut); + struct implied_bound_t { + f_t active_bound; + f_t distance; + uint8_t zero_one_eligible; + }; + + struct implied_bound_group_t { + std::vector bounds; + f_t minimum_alpha; + f_t maximum_alpha; + }; + + struct implied_bound_index_t { + std::vector bounds; + // Nonzero-a candidates only need bounds controlled by binaries in the current row. + std::vector> by_controller; + // The best a=0 bound is invariant for a fixed row coefficient during a cut pass. + std::vector> zero_candidate_cache; + }; + + bool try_add_implied_bound_candidate(const flow_cover_context_t& context, + i_t variable, + f_t coefficient, + bool use_upper_bound, + i_t bound, + f_t binary_coefficient); + void clear_cut_state(i_t num_cols) { continuous_terms.clear(); @@ -531,6 +564,9 @@ class flow_cover_generation_t { std::vector is_slack_; std::vector> flow_cover_constraints_; + implied_bound_index_t upper_implied_bounds; + implied_bound_index_t lower_implied_bounds; + bool cut_pass_preprocessed{false}; std::vector> continuous_terms; std::vector binary_columns; std::vector binary_coefficients; diff --git a/cpp/tests/mip/cuts_test.cu b/cpp/tests/mip/cuts_test.cu index f440eddfc8..db15adaa8a 100644 --- a/cpp/tests/mip/cuts_test.cu +++ b/cpp/tests/mip/cuts_test.cu @@ -1912,6 +1912,8 @@ TEST(cuts, flow_cover_generates_valid_single_node_flow_cut) test_problem.var_types, test_problem.Arow, test_problem.new_slacks); + generator.preprocess_cut_pass( + test_problem.lp, test_problem.settings, variable_bounds, test_problem.var_types, xstar); ASSERT_GT(generator.num_constraints(), 0); int generated_cuts = 0;