From 19c19ac09a485a5f9cf6af6531b401bbf079b626 Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Tue, 22 Sep 2026 17:38:17 +0200 Subject: [PATCH 01/12] activated capacity constraints presolve reduction Signed-off-by: Nicolas L. Guidotti --- .../mathematical_optimization/constants.h | 3 + .../mip/solver_settings.hpp | 8 + cpp/src/math_optimization/solver_settings.cu | 1 + cpp/src/mip_heuristics/CMakeLists.txt | 1 + .../presolve/activated_capacity.cpp | 225 ++++++++++++++++++ .../presolve/activated_capacity.hpp | 51 ++++ .../presolve/third_party_presolve.cpp | 17 +- .../presolve/third_party_presolve.hpp | 3 + cpp/src/mip_heuristics/solve.cu | 3 +- 9 files changed, 307 insertions(+), 5 deletions(-) create mode 100644 cpp/src/mip_heuristics/presolve/activated_capacity.cpp create mode 100644 cpp/src/mip_heuristics/presolve/activated_capacity.hpp diff --git a/cpp/include/cuopt/mathematical_optimization/constants.h b/cpp/include/cuopt/mathematical_optimization/constants.h index 3656791a98..431043f48e 100644 --- a/cpp/include/cuopt/mathematical_optimization/constants.h +++ b/cpp/include/cuopt/mathematical_optimization/constants.h @@ -155,6 +155,9 @@ /* @brief Block bounded-variable-elimination step of cuOpt's internal MIP presolve */ #define CUOPT_MIP_HYPER_BLOCK_BVE "mip_hyper_block_bve" +/* @brief Tie a group capacity row to the group's activation variable during presolve */ +#define CUOPT_MIP_HYPER_ACTIVATED_CAPACITY "mip_hyper_activated_capacity" + /* @brief QCQP (barrier) scaling hyper-parameters */ #define CUOPT_QCQP_HYPER_RUIZ_EQUILIBRATION "qcqp_hyper_ruiz_equilibration" diff --git a/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp b/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp index 7ed45f1f9b..e81a096bac 100644 --- a/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp +++ b/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp @@ -175,6 +175,14 @@ class mip_solver_settings_t { * no-op when no certified reduction exists. */ bool block_bve{true}; + /** + * @brief Strengthen a group capacity row against the group's activation variable. + * + * Where a row caps how many members of a group may be selected and every member is linked to the + * same activation binary by a variable-upper-bound row, rewrites the cap as a multiple of that + * activation, so a fractionally open group is not handed the full capacity allowance. + */ + bool activated_capacity{false}; /** * @brief Determinism mode for MIP solver. * diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index 5a4ab72c32..cd6c04170a 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -260,6 +260,7 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings // Recursive sub-MIP (RINS) hyper-parameters (hidden from default --help: name contains "hyper_") {CUOPT_MIP_HYPER_SUBMIP_ENABLE_CPUFJ, &mip_settings.submip_params.enable_cpufj, true, "run CPU FJ over the sub-MIP"}, {CUOPT_MIP_HYPER_BLOCK_BVE, &mip_settings.block_bve, true, "eliminate blocks of binaries in cuOpt's MIP presolve (needs " CUOPT_MIP_PROBING ")"}, + {CUOPT_MIP_HYPER_ACTIVATED_CAPACITY, &mip_settings.activated_capacity, false, "tie a group capacity row to the group's activation variable during presolve"}, }; // String parameters string_parameters = { diff --git a/cpp/src/mip_heuristics/CMakeLists.txt b/cpp/src/mip_heuristics/CMakeLists.txt index 187017fb14..4196697c04 100644 --- a/cpp/src/mip_heuristics/CMakeLists.txt +++ b/cpp/src/mip_heuristics/CMakeLists.txt @@ -15,6 +15,7 @@ set(MIP_LP_NECESSARY_FILES ${CMAKE_CURRENT_SOURCE_DIR}/presolve/single_lock_dual_aggregation.cpp ${CMAKE_CURRENT_SOURCE_DIR}/presolve/gf2_presolve.cpp ${CMAKE_CURRENT_SOURCE_DIR}/presolve/bhw_coeff_reduce.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/presolve/activated_capacity.cpp ${CMAKE_CURRENT_SOURCE_DIR}/solution/solution.cu ${CMAKE_CURRENT_SOURCE_DIR}/presolve/conflict_graph/clique_table.cu ) diff --git a/cpp/src/mip_heuristics/presolve/activated_capacity.cpp b/cpp/src/mip_heuristics/presolve/activated_capacity.cpp new file mode 100644 index 0000000000..9862330394 --- /dev/null +++ b/cpp/src/mip_heuristics/presolve/activated_capacity.cpp @@ -0,0 +1,225 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#include "activated_capacity.hpp" + +#include +#include +#include + +#include +#include +#include +#include + +// A group capacity row caps how many members of a group may be selected, but says nothing about +// whether the group is open. Where every member is gated by the same activation z, the cap is only +// available once z is paid for: +// +// sum_{i in S} x_i - s <= K, x_i <= z for all i in S, s >= 0 +// => sum_{i in S} x_i - s <= K z. +// +// At z = 0 the gates force every x_i to 0 and the relaxing terms are non-positive, so the +// strengthened row reads 0 - s <= 0; at z = 1 it is the original row. Valid for integral z only, +// which is why this presolver is registered on the MIP path alone. + +namespace cuopt::mathematical_optimization::mip { + +template +papilo::PresolveStatus ActivatedCapacity::execute( + const papilo::Problem& problem, + const papilo::ProblemUpdate& problemUpdate, + const papilo::Num& num, + papilo::Reductions& reductions, + const papilo::Timer& timer, + int& reason_of_infeasibility) +{ + const auto& constraint_matrix = problem.getConstraintMatrix(); + const auto& lhs_values = constraint_matrix.getLeftHandSides(); + const auto& rhs_values = constraint_matrix.getRightHandSides(); + const auto& row_flags = constraint_matrix.getRowFlags(); + const auto& domains = problem.getVariableDomains(); + const auto& col_flags = domains.flags; + const auto& lower_bounds = domains.lower_bounds; + const auto& upper_bounds = domains.upper_bounds; + const auto& presolve_options = problemUpdate.getPresolveOptions(); + + const int num_rows = constraint_matrix.getNRows(); + + auto is_free_binary = [&](int col) { + const auto& flags = col_flags[col]; + return flags.test(papilo::ColFlag::kIntegral) && !flags.test(papilo::ColFlag::kLbInf) && + !flags.test(papilo::ColFlag::kUbInf) && !flags.test(papilo::ColFlag::kFixed) && + num.isZero(lower_bounds[col]) && num.isEq(upper_bounds[col], f_t{1}); + }; + + // Orientation of a one-sided row, or 0 when the row is an equation, a range or free. + // +1 means the stored row reads a.x <= side, -1 means a.x >= side; multiplying by the direction + // puts it in <= form either way. + auto orientation = [&](int row) { + const auto& row_flag = row_flags[row]; + if (row_flag.test(papilo::RowFlag::kRedundant)) return 0; + const bool lhs_infinite = row_flag.test(papilo::RowFlag::kLhsInf); + const bool rhs_infinite = row_flag.test(papilo::RowFlag::kRhsInf); + if (lhs_infinite == rhs_infinite) return 0; + return lhs_infinite ? 1 : -1; + }; + + // Pass one: every two-term row x - z <= 0 over free binaries, as sorted (x, z) pairs. + std::vector> gates; + for (int row = 0; row < num_rows; ++row) { + const int direction = orientation(row); + if (direction == 0) continue; + auto row_coefficients = constraint_matrix.getRowCoefficients(row); + if (row_coefficients.getLength() != 2) continue; + const f_t side = direction == 1 ? rhs_values[row] : lhs_values[row]; + if (!num.isZero(side)) continue; + + const int* indices = row_coefficients.getIndices(); + const f_t* values = row_coefficients.getValues(); + int gated = -1, activation = -1; + for (int j = 0; j < 2; ++j) { + if (!is_free_binary(indices[j])) break; + const f_t v = direction * values[j]; + if (num.isEq(v, f_t{1})) + gated = indices[j]; + else if (num.isEq(v, f_t{-1})) + activation = indices[j]; + } + if (gated >= 0 && activation >= 0) gates.emplace_back(gated, activation); + } + std::sort(gates.begin(), gates.end()); + gates.erase(std::unique(gates.begin(), gates.end()), gates.end()); + + auto gates_of = [&](int col) { + const auto lo = + std::lower_bound(gates.begin(), gates.end(), col, [](const std::pair& g, int c) { + return g.first < c; + }); + const auto hi = std::upper_bound( + lo, gates.end(), col, [](int c, const std::pair& g) { return c < g.first; }); + return std::make_pair(lo, hi); + }; + + papilo::PresolveStatus status = papilo::PresolveStatus::kUnchanged; + int rows_strengthened = 0; + std::vector selections; + std::vector common; + std::vector candidates; + std::vector intersection; + std::vector support; + + // Pass two: the capacity rows themselves. + for (int row = 0; row < num_rows && !gates.empty(); ++row) { + if (reductions.size() >= presolve_options.max_reduction_seq) break; + if (papilo::PresolveMethod::is_interrupted( + timer, presolve_options.tlim, presolve_options.early_exit_callback)) + break; + + const int direction = orientation(row); + if (direction == 0) continue; + auto row_coefficients = constraint_matrix.getRowCoefficients(row); + const int len = row_coefficients.getLength(); + if (len < 2 || len > ACTIVATED_CAPACITY_MAX_LEN) continue; + + const f_t side = direction == 1 ? rhs_values[row] : lhs_values[row]; + const f_t capacity = direction * side; + if (!num.isGT(capacity, f_t{0})) continue; + + const int* indices = row_coefficients.getIndices(); + const f_t* values = row_coefficients.getValues(); + selections.clear(); + support.assign(indices, indices + len); + bool usable = true; + for (int j = 0; j < len && usable; ++j) { + const int col = indices[j]; + const f_t v = direction * values[j]; + if (col_flags[col].test(papilo::ColFlag::kIntegral)) { + // A negative coefficient on an integral column is what this presolver itself writes, so + // rejecting it here is also what keeps a rewritten row from matching a second time. + usable = is_free_binary(col) && num.isEq(v, f_t{1}); + if (usable) selections.push_back(col); + } else { + // A relaxing term: non-positive over the whole box, so it cannot violate the row at z = 0. + usable = num.isLT(v, f_t{0}) && !col_flags[col].test(papilo::ColFlag::kLbInf) && + !num.isLT(lower_bounds[col], f_t{0}); + } + } + if (!usable || selections.size() < 2) continue; + // At or above its own support size the cap is implied by the gates already, and so is K z. + const f_t n_selections = selections.size(); + if (!num.isLT(capacity, n_selections)) continue; + + auto [lo, hi] = gates_of(selections[0]); + common.clear(); + for (auto it = lo; it != hi; ++it) + common.push_back(it->second); + std::sort(common.begin(), common.end()); + for (size_t k = 1; k < selections.size() && !common.empty(); ++k) { + auto [klo, khi] = gates_of(selections[k]); + candidates.clear(); + for (auto it = klo; it != khi; ++it) + candidates.push_back(it->second); + std::sort(candidates.begin(), candidates.end()); + intersection.clear(); + std::set_intersection(common.begin(), + common.end(), + candidates.begin(), + candidates.end(), + std::back_inserter(intersection)); + common.swap(intersection); + } + if (common.empty()) continue; + + std::sort(support.begin(), support.end()); + int activation = -1; + for (int z : common) { + if (std::binary_search(support.begin(), support.end(), z)) continue; + activation = z; + break; + } + if (activation < 0) continue; + + cuopt_assert(is_free_binary(activation), "the activation of a gate row is a free binary"); + + papilo::TransactionGuard guard{reductions}; + reductions.lockRow(row); + reductions.changeMatrixEntry(row, activation, direction * -capacity); + if (direction == 1) + reductions.changeRowRHS(row, f_t{0}); + else + reductions.changeRowLHS(row, f_t{0}); + ++rows_strengthened; + status = papilo::PresolveStatus::kReduced; + } + + // Proposed, not applied: the activation is a nonzero the row does not have yet, and + // ConstraintMatrix::change_coefficient refuses the insert when the row or the column has no slack + // space left in its range, which rejects the whole transaction. The presolved nonzero count is + // the figure to check this against. + if (rows_strengthened > 0) { + CUOPT_LOG_INFO("Activated capacity: proposed %d strengthened rows against %zu gates", + rows_strengthened, + gates.size()); + } + + return status; +} + +#define INSTANTIATE(F_TYPE) template class ActivatedCapacity; + +#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT +INSTANTIATE(float) +#endif + +#if MIP_INSTANTIATE_DOUBLE +INSTANTIATE(double) +#endif + +#undef INSTANTIATE + +} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/presolve/activated_capacity.hpp b/cpp/src/mip_heuristics/presolve/activated_capacity.hpp new file mode 100644 index 0000000000..6d1f3e8e35 --- /dev/null +++ b/cpp/src/mip_heuristics/presolve/activated_capacity.hpp @@ -0,0 +1,51 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#pragma once + +#if !defined(__clang__) +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wstringop-overflow" // ignore boost error for pip wheel build +#pragma GCC diagnostic ignored "-Wnarrowing" +#endif +#include +#include +#include +#include +#if !defined(__clang__) +#pragma GCC diagnostic pop +#endif + +#include + +namespace cuopt::mathematical_optimization::mip { + +// A capacity row wider than this is not worth the gate-set intersection. +static constexpr int ACTIVATED_CAPACITY_MAX_LEN = 4096; + +template +class ActivatedCapacity : public papilo::PresolveMethod { + public: + ActivatedCapacity() : papilo::PresolveMethod() + { + this->setName("activatedcapacity"); + this->setType(papilo::PresolverType::kIntegralCols); + this->setTiming(papilo::PresolverTiming::kMedium); + // The vacuous copies of these capacity rows are dropped by the cheaper presolvers first, which + // keeps the gate-set intersection off rows that would gain nothing. + this->setDelayed(true); + } + + papilo::PresolveStatus execute(const papilo::Problem& problem, + const papilo::ProblemUpdate& problemUpdate, + const papilo::Num& num, + papilo::Reductions& reductions, + const papilo::Timer& timer, + int& reason_of_infeasibility) override; +}; + +} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp index 7bf7bd76d6..a9e27721e0 100644 --- a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp +++ b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp @@ -40,6 +40,7 @@ #include #include #include +#include #include #include #include @@ -731,7 +732,8 @@ void set_presolve_methods( papilo::Presolve& presolver, problem_category_t category, bool dual_postsolve, - std::optional> const& method_allowlist = std::nullopt) + std::optional> const& method_allowlist = std::nullopt, + bool activated_capacity = false) { using uptr = std::unique_ptr>; @@ -747,6 +749,9 @@ void set_presolve_methods( // cuOpt custom GF2 presolver maybe_add(uptr(new cuopt::mathematical_optimization::mip::GF2Presolve())); maybe_add(uptr(new cuopt::mathematical_optimization::mip::BHWCoeffReduce())); + if (activated_capacity) { + maybe_add(uptr(new cuopt::mathematical_optimization::mip::ActivatedCapacity())); + } } // fast presolvers maybe_add(uptr(new papilo::SingletonCols())); @@ -935,7 +940,8 @@ third_party_presolve_status_t third_party_presolve_t::apply_papilo( CUOPT_LOG_INFO("\nRunning Papilo presolve (git hash %s)", PAPILO_GITHASH); if (category == problem_category_t::MIP) { dual_postsolve = false; } papilo::Presolve papilo_presolver; - set_presolve_methods(papilo_presolver, category, dual_postsolve, reduction_allowlist_); + set_presolve_methods( + papilo_presolver, category, dual_postsolve, reduction_allowlist_, activated_capacity_); set_presolve_options(papilo_presolver, category, absolute_tolerance, @@ -1221,8 +1227,11 @@ third_party_presolve_status_t third_party_presolve_t::apply_to_subprob papilo_problem.getConstraintMatrix().getNnz()); papilo::Presolve papilo_presolver; - set_presolve_methods( - papilo_presolver, problem_category_t::MIP, dual_postsolve, reduction_allowlist_); + set_presolve_methods(papilo_presolver, + problem_category_t::MIP, + dual_postsolve, + reduction_allowlist_, + activated_capacity_); set_presolve_options(papilo_presolver, problem_category_t::MIP, settings.primal_tol, diff --git a/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp b/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp index f1426749c3..c81b6d55a2 100644 --- a/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp +++ b/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp @@ -115,6 +115,8 @@ class third_party_presolve_t { reduction_allowlist_ = std::move(allowlist); } + void set_activated_capacity(bool enabled) { activated_capacity_ = enabled; } + // Apply the presolve on an simplex::user_problem in-place. Used in sub MIP and (in the future) // restarts. third_party_presolve_status_t apply_to_subproblem( @@ -225,6 +227,7 @@ class third_party_presolve_t { f_t original_objective_scaling_factor_{1}; std::optional> reduction_allowlist_{}; + bool activated_capacity_{false}; }; // Just for testing the conversion: user_problem -> Papilo problem -> user_problem. diff --git a/cpp/src/mip_heuristics/solve.cu b/cpp/src/mip_heuristics/solve.cu index 6a6c6d795c..f6bbdb5691 100644 --- a/cpp/src/mip_heuristics/solve.cu +++ b/cpp/src/mip_heuristics/solve.cu @@ -621,7 +621,8 @@ mip_solution_t solve_mip_helper( ? std::numeric_limits::infinity() : timer.remaining_time(); - presolver = std::make_unique>(); + presolver = std::make_unique>(); + presolver->set_activated_capacity(settings.activated_capacity); auto result = presolver->apply_presolve_from_op_problem( op_problem, cuopt::mathematical_optimization::problem_category_t::MIP, From c57c8a1b6d1e25320228187f08666ab7d772b6d8 Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Tue, 22 Sep 2026 17:36:48 +0200 Subject: [PATCH 02/12] group cover cuts Signed-off-by: Nicolas L. Guidotti --- .../mathematical_optimization/constants.h | 1 + .../mip/solver_settings.hpp | 3 + cpp/src/branch_and_bound/branch_and_bound.cpp | 1 + cpp/src/cuts/cuts.cpp | 196 ++++++++++++++++++ cpp/src/cuts/cuts.hpp | 21 +- .../dual_simplex/simplex_solver_settings.hpp | 2 + cpp/src/math_optimization/solver_settings.cu | 1 + .../diversity/recombiners/sub_mip.cuh | 1 + cpp/src/mip_heuristics/solver.cu | 1 + 9 files changed, 225 insertions(+), 2 deletions(-) diff --git a/cpp/include/cuopt/mathematical_optimization/constants.h b/cpp/include/cuopt/mathematical_optimization/constants.h index 431043f48e..cce4e79abb 100644 --- a/cpp/include/cuopt/mathematical_optimization/constants.h +++ b/cpp/include/cuopt/mathematical_optimization/constants.h @@ -80,6 +80,7 @@ #define CUOPT_MIP_CLIQUE_CUTS "mip_clique_cuts" #define CUOPT_MIP_ZERO_HALF_CUTS "mip_zero_half_cuts" #define CUOPT_MIP_STRONG_CHVATAL_GOMORY_CUTS "mip_strong_chvatal_gomory_cuts" +#define CUOPT_MIP_GROUP_COVER_CUTS "mip_group_cover_cuts" #define CUOPT_MIP_REDUCED_COST_STRENGTHENING "mip_reduced_cost_strengthening" #define CUOPT_MIP_RINS "mip_rins" #define CUOPT_MIP_RENS "mip_rens" diff --git a/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp b/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp index e81a096bac..52f6e1a0bf 100644 --- a/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp +++ b/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp @@ -135,6 +135,9 @@ class mip_solver_settings_t { i_t clique_cuts = -1; i_t zero_half_cuts = -1; i_t implied_bound_cuts = -1; + // Aggregate an enabler row through its variable-upper-bound gates, counting each group once. + // 0 = disable, >0 = enable. Off by default while the separator is a prototype. + i_t group_cover_cuts = 0; i_t strong_chvatal_gomory_cuts = -1; i_t reduced_cost_strengthening = -1; i_t objective_step = 1; // 0 = disable objective step tightening, 1 = enable diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index a320cc0602..f2ae9c6761 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -2321,6 +2321,7 @@ void branch_and_bound_t::solve_submip(diving_worker_t* worke submip_settings.reliability_branching = 0; submip_settings.clique_cuts = 0; submip_settings.zero_half_cuts = 0; + submip_settings.group_cover_cuts = 0; submip_settings.inside_submip = 1; submip_settings.strong_branching_simplex_iteration_limit = 50; submip_settings.inside_root_node = 0; diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index e9f51666dc..47c4e9a9dc 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -3232,6 +3232,191 @@ void cut_generation_t::generate_implied_bound_cuts( } } +// A group cover cut aggregates an enabler row through the gates behind it. Where the model carries +// +// y <= sum_{j in S} x_j (enabler) and x_j <= z_{g(j)} for every j in S, +// +// a binary y that is one forces some x_j to one, which forces its own group activation to one, so +// +// y <= sum_{g in D} z_g, D = the distinct groups covering S. +// +// Counting each group once is where the strength is: the enabler alone lets y reach one against a +// whole group held at 1/|S|, while the cut holds y down to that group's own activation. Valid for +// integral y and x only. +template +void cut_generation_t::build_group_cover_candidates( + const simplex_solver_settings_t& settings) +{ + group_cover_built_ = true; + + const i_t num_rows = user_problem_.num_rows; + const i_t num_cols = user_problem_.num_cols; + if (num_rows <= 0 || num_cols <= 0) { return; } + if (user_problem_.var_types.size() != (size_t)num_cols || + user_problem_.row_sense.size() != (size_t)num_rows || + user_problem_.rhs.size() != (size_t)num_rows) { + return; + } + + csr_matrix_t Arow(num_rows, num_cols, user_problem_.A.col_start[num_cols]); + user_problem_.A.to_compressed_row(Arow); + + auto is_binary = [&](i_t col) { + return user_problem_.var_types[col] != variable_type_t::CONTINUOUS && + user_problem_.lower[col] == 0.0 && user_problem_.upper[col] == 1.0; + }; + // Row sense in <= orientation: +1 when the row reads a.x <= rhs, -1 when a.x >= rhs, 0 otherwise. + auto direction_of = [&](i_t row) { + if (user_problem_.row_sense[row] == 'L') { return 1; } + if (user_problem_.row_sense[row] == 'G') { return -1; } + return 0; + }; + const bool has_ranges = user_problem_.num_range_rows > 0; + std::vector is_range(has_ranges ? num_rows : 0, 0); + for (i_t k = 0; k < user_problem_.num_range_rows; ++k) { + is_range[user_problem_.range_rows[k]] = 1; + } + + // Pass one: the gates, as (selection, activation) pairs. + std::vector> gates; + for (i_t row = 0; row < num_rows; ++row) { + if (has_ranges && is_range[row]) { continue; } + const i_t direction = direction_of(row); + if (direction == 0) { continue; } + if (Arow.row_length(row) != 2) { continue; } + if (direction * user_problem_.rhs[row] != 0.0) { continue; } + + i_t gated = -1, activation = -1; + for (i_t p = Arow.row_start[row]; p < Arow.row_start[row + 1]; ++p) { + const i_t col = Arow.j[p]; + if (!is_binary(col)) { + gated = activation = -1; + break; + } + const f_t v = direction * Arow.x[p]; + if (v == 1.0) { + gated = col; + } else if (v == -1.0) { + activation = col; + } + } + if (gated >= 0 && activation >= 0) { gates.emplace_back(gated, activation); } + } + if (gates.empty()) { return; } + std::sort(gates.begin(), gates.end()); + gates.erase(std::unique(gates.begin(), gates.end()), gates.end()); + + auto gates_of = [&](i_t col) { + const auto lo = + std::lower_bound(gates.begin(), gates.end(), col, [](const std::pair& g, i_t c) { + return g.first < c; + }); + const auto hi = std::upper_bound( + lo, gates.end(), col, [](i_t c, const std::pair& g) { return c < g.first; }); + return std::make_pair(lo, hi); + }; + + // Pass two: the enabler rows, one +1 head against a tail of -1 selections. + std::vector groups; + std::vector>> candidates; + for (i_t row = 0; row < num_rows; ++row) { + if (has_ranges && is_range[row]) { continue; } + const i_t direction = direction_of(row); + if (direction == 0) { continue; } + const i_t len = Arow.row_length(row); + if (len < 3) { continue; } + if (direction * user_problem_.rhs[row] != 0.0) { continue; } + + i_t head = -1; + bool usable = true; + groups.clear(); + for (i_t p = Arow.row_start[row]; p < Arow.row_start[row + 1] && usable; ++p) { + const i_t col = Arow.j[p]; + const f_t v = direction * Arow.x[p]; + if (!is_binary(col)) { + usable = false; + } else if (v == 1.0) { + usable = head < 0; + head = col; + } else if (v == -1.0) { + auto [lo, hi] = gates_of(col); + // A tail member with no gate leaves nothing to aggregate through; keep the candidate by + // standing in the selection itself, which is what the enabler already bounds the head by. + if (lo == hi) { + groups.push_back(col); + } else { + for (auto it = lo; it != hi; ++it) { + groups.push_back(it->second); + } + } + } else { + usable = false; + } + } + if (!usable || head < 0 || groups.empty()) { continue; } + + std::sort(groups.begin(), groups.end()); + groups.erase(std::unique(groups.begin(), groups.end()), groups.end()); + // Nothing was merged, so the cut is the sum of the gates the LP already has. + if (groups.size() >= (size_t)(len - 1)) { continue; } + if (std::binary_search(groups.begin(), groups.end(), head)) { continue; } + + candidates.emplace_back(head, groups); + } + if (candidates.empty()) { return; } + + // Two enabler rows over the same group set give the same cut; keep one. + std::sort(candidates.begin(), candidates.end()); + candidates.erase(std::unique(candidates.begin(), candidates.end()), candidates.end()); + + group_cover_heads_.reserve(candidates.size()); + group_cover_offsets_.reserve(candidates.size() + 1); + group_cover_offsets_.push_back(0); + for (const auto& [head, group_set] : candidates) { + group_cover_heads_.push_back(head); + group_cover_groups_.insert(group_cover_groups_.end(), group_set.begin(), group_set.end()); + group_cover_offsets_.push_back(group_cover_groups_.size()); + } + + settings.log.print_format( + "Group cover: {} candidate cuts over {} gates\n", group_cover_heads_.size(), gates.size()); +} + +template +void cut_generation_t::generate_group_cover_cuts( + const simplex_solver_settings_t& settings, + const std::vector& xstar, + f_t start_time) +{ + if (!group_cover_built_) { build_group_cover_candidates(settings); } + if (group_cover_heads_.empty()) { return; } + + const f_t tol = 1e-4; + i_t num_cuts = 0; + const i_t n = group_cover_heads_.size(); + for (i_t k = 0; k < n; ++k) { + if ((k & 0xFF) == 0 && toc(start_time) >= settings.time_limit) { return; } + const i_t head = group_cover_heads_[k]; + f_t activity = xstar[head]; + for (i_t p = group_cover_offsets_[k]; p < group_cover_offsets_[k + 1]; ++p) { + activity -= xstar[group_cover_groups_[p]]; + } + if (activity <= tol) { continue; } + + // add_cut expects cut'x >= rhs, so the cut head - sum_g z_g <= 0 is emitted negated. + inequality_t cut; + cut.push_back(head, -1.0); + for (i_t p = group_cover_offsets_[k]; p < group_cover_offsets_[k + 1]; ++p) { + cut.push_back(group_cover_groups_[p], 1.0); + } + cut.rhs = 0.0; + cut_pool_.add_cut(cut_type_t::GROUP_COVER, cut); + num_cuts++; + } + + if (num_cuts > 0) { settings.log.debug("Generated %d group cover cuts\n", num_cuts); } +} + namespace { // Total probing-edge budget from the byte cap and the remaining work headroom @@ -3606,6 +3791,17 @@ bool cut_generation_t::generate_cuts(const lp_problem_t& lp, } } + // Generate group cover cuts + if (settings.group_cover_cuts != 0) { + if (toc(start_time) >= settings.time_limit) { return true; } + f_t cut_start_time = tic(); + generate_group_cover_cuts(settings, xstar, start_time); + f_t cut_generation_time = toc(cut_start_time); + if (cut_generation_time > 1.0) { + settings.log.debug("Group cover cut generation time %.2f seconds\n", cut_generation_time); + } + } + // Build the fractional conflict-graph subgraph once (resolving the async // clique-table future on the way) so both clique-cut and zero-half cut // separators consume the same vertex/weight/adjacency tables instead of diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index ca87e26c39..a4cf3d592d 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -45,7 +45,8 @@ enum cut_type_t : int8_t { IMPLIED_BOUND = 5, ZERO_HALF = 6, FLOW_COVER = 7, - MAX_CUT_TYPE = 8 + GROUP_COVER = 8, + MAX_CUT_TYPE = 9 }; template @@ -186,7 +187,8 @@ struct cut_info_t { "Clique ", "Implied Bounds", "Zero-Half ", - "Flow Cover "}; + "Flow Cover ", + "Group Cover "}; std::array num_cuts = {0}; }; @@ -861,6 +863,14 @@ class cut_generation_t { const std::vector& xstar, f_t start_time); + // Generate group cover cuts from the variable-upper-bound gates behind each enabler row + void generate_group_cover_cuts(const simplex::simplex_solver_settings_t& settings, + const std::vector& xstar, + f_t start_time); + + // Scan the user problem for gate and enabler rows. Called once, on the first cut pass. + void build_group_cover_candidates(const simplex::simplex_solver_settings_t& settings); + void prepare_fractional_sub_conflict_graph( const simplex::simplex_solver_settings_t& settings, const std::vector& xstar, @@ -876,6 +886,13 @@ class cut_generation_t { std::shared_ptr>& clique_table_; omp_atomic_t* signal_extend_{nullptr}; fractional_conflict_subgraph_t sub_cg_; + // One candidate per enabler row whose selections are all gated: the head, and the distinct group + // activations behind its tail. Built once from user_problem_, then separated against xstar each + // pass. group_cover_groups_ is the flat store the offsets index into. + std::vector group_cover_heads_; + std::vector group_cover_offsets_; + std::vector group_cover_groups_; + bool group_cover_built_{false}; }; template diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index 8b3eba56d3..64cb7f5726 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -98,6 +98,7 @@ struct simplex_solver_settings_t { knapsack_cuts(-1), flow_cover_cuts(-1), implied_bound_cuts(-1), + group_cover_cuts(-1), clique_cuts(-1), zero_half_cuts(-1), strong_chvatal_gomory_cuts(-1), @@ -210,6 +211,7 @@ struct simplex_solver_settings_t { i_t knapsack_cuts; // -1 automatic, 0 to disable, >0 to enable knapsack cuts i_t flow_cover_cuts; // -1 automatic, 0 to disable, >0 to enable flow cover cuts i_t implied_bound_cuts; // -1 automatic, 0 to disable, >0 to enable implied bound cuts + i_t group_cover_cuts; // 0 to disable, >0 to enable group cover cuts i_t clique_cuts; // -1 automatic, 0 to disable, >0 to enable clique cuts i_t zero_half_cuts; // -1 automatic, 0 to disable, >0 to enable zero-half cuts i_t strong_chvatal_gomory_cuts; // -1 automatic, 0 to disable, >0 to enable strong Chvatal Gomory diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index cd6c04170a..75e2166ac5 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -189,6 +189,7 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_MIP_CLIQUE_CUTS, &mip_settings.clique_cuts, -1, 1, -1}, {CUOPT_MIP_ZERO_HALF_CUTS, &mip_settings.zero_half_cuts, -1, 1, -1}, {CUOPT_MIP_IMPLIED_BOUND_CUTS, &mip_settings.implied_bound_cuts, -1, 1, -1}, + {CUOPT_MIP_GROUP_COVER_CUTS, &mip_settings.group_cover_cuts, -1, 1, -1}, {CUOPT_MIP_STRONG_CHVATAL_GOMORY_CUTS, &mip_settings.strong_chvatal_gomory_cuts, -1, 1, -1}, {CUOPT_MIP_REDUCED_COST_STRENGTHENING, &mip_settings.reduced_cost_strengthening, -1, std::numeric_limits::max(), -1}, {CUOPT_MIP_RINS, &mip_settings.submip_params.rins, -1, 1, -1}, diff --git a/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh b/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh index 7990465b81..643bf694b4 100644 --- a/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh +++ b/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh @@ -112,6 +112,7 @@ class sub_mip_recombiner_t : public recombiner_t { branch_and_bound_settings.max_cut_passes = 0; branch_and_bound_settings.clique_cuts = 0; branch_and_bound_settings.zero_half_cuts = 0; + branch_and_bound_settings.group_cover_cuts = 0; branch_and_bound_settings.inside_submip = 1; branch_and_bound_settings.submip_settings.rins = 0; branch_and_bound_settings.submip_settings.rens = 0; diff --git a/cpp/src/mip_heuristics/solver.cu b/cpp/src/mip_heuristics/solver.cu index 166afd832b..535163d53c 100644 --- a/cpp/src/mip_heuristics/solver.cu +++ b/cpp/src/mip_heuristics/solver.cu @@ -369,6 +369,7 @@ solution_t mip_solver_t::run_solver() branch_and_bound_settings.knapsack_cuts = context.settings.knapsack_cuts; branch_and_bound_settings.flow_cover_cuts = context.settings.flow_cover_cuts; branch_and_bound_settings.implied_bound_cuts = context.settings.implied_bound_cuts; + branch_and_bound_settings.group_cover_cuts = context.settings.group_cover_cuts; branch_and_bound_settings.clique_cuts = context.settings.clique_cuts; branch_and_bound_settings.zero_half_cuts = context.settings.zero_half_cuts; branch_and_bound_settings.strong_chvatal_gomory_cuts = From 9bb432ffe80a9e0ccbf8e073ca1db3cd49fe5a54 Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Tue, 22 Sep 2026 18:05:44 +0200 Subject: [PATCH 03/12] enable the group cover cuts and activated capacity presolve by default. removed unnecessary parameters. Signed-off-by: Nicolas L. Guidotti --- .../mathematical_optimization/constants.h | 3 --- .../mip/solver_settings.hpp | 12 +---------- cpp/src/math_optimization/solver_settings.cu | 1 - .../presolve/third_party_presolve.cpp | 21 +++++++------------ .../presolve/third_party_presolve.hpp | 3 --- cpp/src/mip_heuristics/solve.cu | 3 +-- 6 files changed, 9 insertions(+), 34 deletions(-) diff --git a/cpp/include/cuopt/mathematical_optimization/constants.h b/cpp/include/cuopt/mathematical_optimization/constants.h index cce4e79abb..ad365fc0d1 100644 --- a/cpp/include/cuopt/mathematical_optimization/constants.h +++ b/cpp/include/cuopt/mathematical_optimization/constants.h @@ -156,9 +156,6 @@ /* @brief Block bounded-variable-elimination step of cuOpt's internal MIP presolve */ #define CUOPT_MIP_HYPER_BLOCK_BVE "mip_hyper_block_bve" -/* @brief Tie a group capacity row to the group's activation variable during presolve */ -#define CUOPT_MIP_HYPER_ACTIVATED_CAPACITY "mip_hyper_activated_capacity" - /* @brief QCQP (barrier) scaling hyper-parameters */ #define CUOPT_QCQP_HYPER_RUIZ_EQUILIBRATION "qcqp_hyper_ruiz_equilibration" diff --git a/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp b/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp index 52f6e1a0bf..4a60ac991c 100644 --- a/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp +++ b/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp @@ -135,9 +135,7 @@ class mip_solver_settings_t { i_t clique_cuts = -1; i_t zero_half_cuts = -1; i_t implied_bound_cuts = -1; - // Aggregate an enabler row through its variable-upper-bound gates, counting each group once. - // 0 = disable, >0 = enable. Off by default while the separator is a prototype. - i_t group_cover_cuts = 0; + i_t group_cover_cuts = -1; i_t strong_chvatal_gomory_cuts = -1; i_t reduced_cost_strengthening = -1; i_t objective_step = 1; // 0 = disable objective step tightening, 1 = enable @@ -178,14 +176,6 @@ class mip_solver_settings_t { * no-op when no certified reduction exists. */ bool block_bve{true}; - /** - * @brief Strengthen a group capacity row against the group's activation variable. - * - * Where a row caps how many members of a group may be selected and every member is linked to the - * same activation binary by a variable-upper-bound row, rewrites the cap as a multiple of that - * activation, so a fractionally open group is not handed the full capacity allowance. - */ - bool activated_capacity{false}; /** * @brief Determinism mode for MIP solver. * diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index 75e2166ac5..1a8908d562 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -261,7 +261,6 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings // Recursive sub-MIP (RINS) hyper-parameters (hidden from default --help: name contains "hyper_") {CUOPT_MIP_HYPER_SUBMIP_ENABLE_CPUFJ, &mip_settings.submip_params.enable_cpufj, true, "run CPU FJ over the sub-MIP"}, {CUOPT_MIP_HYPER_BLOCK_BVE, &mip_settings.block_bve, true, "eliminate blocks of binaries in cuOpt's MIP presolve (needs " CUOPT_MIP_PROBING ")"}, - {CUOPT_MIP_HYPER_ACTIVATED_CAPACITY, &mip_settings.activated_capacity, false, "tie a group capacity row to the group's activation variable during presolve"}, }; // String parameters string_parameters = { diff --git a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp index a9e27721e0..b3f35e0b40 100644 --- a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp +++ b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp @@ -732,8 +732,7 @@ void set_presolve_methods( papilo::Presolve& presolver, problem_category_t category, bool dual_postsolve, - std::optional> const& method_allowlist = std::nullopt, - bool activated_capacity = false) + std::optional> const& method_allowlist = std::nullopt) { using uptr = std::unique_ptr>; @@ -747,11 +746,9 @@ void set_presolve_methods( if (category == problem_category_t::MIP) { // cuOpt custom GF2 presolver - maybe_add(uptr(new cuopt::mathematical_optimization::mip::GF2Presolve())); - maybe_add(uptr(new cuopt::mathematical_optimization::mip::BHWCoeffReduce())); - if (activated_capacity) { - maybe_add(uptr(new cuopt::mathematical_optimization::mip::ActivatedCapacity())); - } + maybe_add(uptr(new GF2Presolve())); + maybe_add(uptr(new BHWCoeffReduce())); + maybe_add(uptr(new ActivatedCapacity())); } // fast presolvers maybe_add(uptr(new papilo::SingletonCols())); @@ -940,8 +937,7 @@ third_party_presolve_status_t third_party_presolve_t::apply_papilo( CUOPT_LOG_INFO("\nRunning Papilo presolve (git hash %s)", PAPILO_GITHASH); if (category == problem_category_t::MIP) { dual_postsolve = false; } papilo::Presolve papilo_presolver; - set_presolve_methods( - papilo_presolver, category, dual_postsolve, reduction_allowlist_, activated_capacity_); + set_presolve_methods(papilo_presolver, category, dual_postsolve, reduction_allowlist_); set_presolve_options(papilo_presolver, category, absolute_tolerance, @@ -1227,11 +1223,8 @@ third_party_presolve_status_t third_party_presolve_t::apply_to_subprob papilo_problem.getConstraintMatrix().getNnz()); papilo::Presolve papilo_presolver; - set_presolve_methods(papilo_presolver, - problem_category_t::MIP, - dual_postsolve, - reduction_allowlist_, - activated_capacity_); + set_presolve_methods( + papilo_presolver, problem_category_t::MIP, dual_postsolve, reduction_allowlist_); set_presolve_options(papilo_presolver, problem_category_t::MIP, settings.primal_tol, diff --git a/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp b/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp index c81b6d55a2..f1426749c3 100644 --- a/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp +++ b/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp @@ -115,8 +115,6 @@ class third_party_presolve_t { reduction_allowlist_ = std::move(allowlist); } - void set_activated_capacity(bool enabled) { activated_capacity_ = enabled; } - // Apply the presolve on an simplex::user_problem in-place. Used in sub MIP and (in the future) // restarts. third_party_presolve_status_t apply_to_subproblem( @@ -227,7 +225,6 @@ class third_party_presolve_t { f_t original_objective_scaling_factor_{1}; std::optional> reduction_allowlist_{}; - bool activated_capacity_{false}; }; // Just for testing the conversion: user_problem -> Papilo problem -> user_problem. diff --git a/cpp/src/mip_heuristics/solve.cu b/cpp/src/mip_heuristics/solve.cu index f6bbdb5691..6a6c6d795c 100644 --- a/cpp/src/mip_heuristics/solve.cu +++ b/cpp/src/mip_heuristics/solve.cu @@ -621,8 +621,7 @@ mip_solution_t solve_mip_helper( ? std::numeric_limits::infinity() : timer.remaining_time(); - presolver = std::make_unique>(); - presolver->set_activated_capacity(settings.activated_capacity); + presolver = std::make_unique>(); auto result = presolver->apply_presolve_from_op_problem( op_problem, cuopt::mathematical_optimization::problem_category_t::MIP, From ed6e7ba391348a45b47e392cb04f90ac9b20ef17 Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Wed, 23 Sep 2026 16:00:06 +0200 Subject: [PATCH 04/12] converted to activated capacity to a cut Signed-off-by: Nicolas L. Guidotti --- .../mathematical_optimization/constants.h | 1 + .../mip/solver_settings.hpp | 1 + cpp/src/branch_and_bound/branch_and_bound.cpp | 1 + cpp/src/cuts/cuts.cpp | 317 ++++++++++++++---- cpp/src/cuts/cuts.hpp | 33 +- .../dual_simplex/simplex_solver_settings.hpp | 2 + cpp/src/math_optimization/solver_settings.cu | 1 + cpp/src/mip_heuristics/CMakeLists.txt | 1 - .../diversity/recombiners/sub_mip.cuh | 1 + .../presolve/activated_capacity.cpp | 225 ------------- .../presolve/activated_capacity.hpp | 51 --- .../presolve/third_party_presolve.cpp | 2 - cpp/src/mip_heuristics/solver.cu | 13 +- 13 files changed, 304 insertions(+), 345 deletions(-) delete mode 100644 cpp/src/mip_heuristics/presolve/activated_capacity.cpp delete mode 100644 cpp/src/mip_heuristics/presolve/activated_capacity.hpp diff --git a/cpp/include/cuopt/mathematical_optimization/constants.h b/cpp/include/cuopt/mathematical_optimization/constants.h index ad365fc0d1..888eed9a15 100644 --- a/cpp/include/cuopt/mathematical_optimization/constants.h +++ b/cpp/include/cuopt/mathematical_optimization/constants.h @@ -81,6 +81,7 @@ #define CUOPT_MIP_ZERO_HALF_CUTS "mip_zero_half_cuts" #define CUOPT_MIP_STRONG_CHVATAL_GOMORY_CUTS "mip_strong_chvatal_gomory_cuts" #define CUOPT_MIP_GROUP_COVER_CUTS "mip_group_cover_cuts" +#define CUOPT_MIP_ACTIVATED_CAPACITY_CUTS "mip_activated_capacity_cuts" #define CUOPT_MIP_REDUCED_COST_STRENGTHENING "mip_reduced_cost_strengthening" #define CUOPT_MIP_RINS "mip_rins" #define CUOPT_MIP_RENS "mip_rens" diff --git a/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp b/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp index 4a60ac991c..e177623549 100644 --- a/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp +++ b/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp @@ -136,6 +136,7 @@ class mip_solver_settings_t { i_t zero_half_cuts = -1; i_t implied_bound_cuts = -1; i_t group_cover_cuts = -1; + i_t activated_capacity_cuts = -1; i_t strong_chvatal_gomory_cuts = -1; i_t reduced_cost_strengthening = -1; i_t objective_step = 1; // 0 = disable objective step tightening, 1 = enable diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index f2ae9c6761..0e778d2eea 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -2322,6 +2322,7 @@ void branch_and_bound_t::solve_submip(diving_worker_t* worke submip_settings.clique_cuts = 0; submip_settings.zero_half_cuts = 0; submip_settings.group_cover_cuts = 0; + submip_settings.activated_capacity_cuts = 0; submip_settings.inside_submip = 1; submip_settings.strong_branching_simplex_iteration_limit = 50; submip_settings.inside_root_node = 0; diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index 47c4e9a9dc..26c27a2af4 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -3232,6 +3232,81 @@ void cut_generation_t::generate_implied_bound_cuts( } } +template +void cut_generation_t::build_gate_table(const csr_matrix_t& Arow) +{ + if (gates_built_) { return; } + gates_built_ = true; + + const i_t num_rows = user_problem_.num_rows; + const i_t num_cols = user_problem_.num_cols; + + const bool has_ranges = user_problem_.num_range_rows > 0; + std::vector is_range(has_ranges ? num_rows : 0, 0); + for (i_t k = 0; k < user_problem_.num_range_rows; ++k) { + is_range[user_problem_.range_rows[k]] = 1; + } + + std::vector> gates; + for (i_t row = 0; row < num_rows; ++row) { + if (has_ranges && is_range[row]) { continue; } + if (Arow.row_length(row) != 2) { continue; } + const char sense = user_problem_.row_sense[row]; + if (sense != 'L' && sense != 'G') { continue; } + // Row sense in <= orientation: +1 when the row reads a.x <= rhs, -1 when a.x >= rhs. + const i_t direction = sense == 'L' ? 1 : -1; + if (direction * user_problem_.rhs[row] != 0.0) { continue; } + + i_t gated = -1, activation = -1; + for (i_t p = Arow.row_start[row]; p < Arow.row_start[row + 1]; ++p) { + const i_t col = Arow.j[p]; + if (user_problem_.var_types[col] == variable_type_t::CONTINUOUS || + user_problem_.lower[col] != 0.0 || user_problem_.upper[col] != 1.0) { + gated = activation = -1; + break; + } + const f_t v = direction * Arow.x[p]; + if (v == 1.0) { + gated = col; + } else if (v == -1.0) { + activation = col; + } + } + if (gated >= 0 && activation >= 0) { gates.emplace_back(gated, activation); } + } + if (gates.empty()) { return; } + + gate_offsets_.assign(num_cols + 1, 0); + for (const auto& [gated, activation] : gates) { + ++gate_offsets_[gated + 1]; + } + for (i_t col = 0; col < num_cols; ++col) { + gate_offsets_[col + 1] += gate_offsets_[col]; + } + gate_activations_.resize(gates.size()); + std::vector cursor(gate_offsets_.begin(), gate_offsets_.end() - 1); + for (const auto& [gated, activation] : gates) { + gate_activations_[cursor[gated]++] = activation; + } + + // Sort each column's span and drop duplicate gate rows, compacting in place. The write cursor + // trails the read cursor because it only advances on a kept entry. + i_t out = 0; + for (i_t col = 0; col < num_cols; ++col) { + const i_t begin = gate_offsets_[col]; + const i_t end = gate_offsets_[col + 1]; + gate_offsets_[col] = out; + std::sort(gate_activations_.begin() + begin, gate_activations_.begin() + end); + for (i_t p = begin; p < end; ++p) { + if (p == begin || gate_activations_[p] != gate_activations_[p - 1]) { + gate_activations_[out++] = gate_activations_[p]; + } + } + } + gate_offsets_[num_cols] = out; + gate_activations_.resize(out); +} + // A group cover cut aggregates an enabler row through the gates behind it. Where the model carries // // y <= sum_{j in S} x_j (enabler) and x_j <= z_{g(j)} for every j in S, @@ -3261,69 +3336,25 @@ void cut_generation_t::build_group_cover_candidates( csr_matrix_t Arow(num_rows, num_cols, user_problem_.A.col_start[num_cols]); user_problem_.A.to_compressed_row(Arow); - auto is_binary = [&](i_t col) { - return user_problem_.var_types[col] != variable_type_t::CONTINUOUS && - user_problem_.lower[col] == 0.0 && user_problem_.upper[col] == 1.0; - }; - // Row sense in <= orientation: +1 when the row reads a.x <= rhs, -1 when a.x >= rhs, 0 otherwise. - auto direction_of = [&](i_t row) { - if (user_problem_.row_sense[row] == 'L') { return 1; } - if (user_problem_.row_sense[row] == 'G') { return -1; } - return 0; - }; const bool has_ranges = user_problem_.num_range_rows > 0; std::vector is_range(has_ranges ? num_rows : 0, 0); for (i_t k = 0; k < user_problem_.num_range_rows; ++k) { is_range[user_problem_.range_rows[k]] = 1; } - // Pass one: the gates, as (selection, activation) pairs. - std::vector> gates; - for (i_t row = 0; row < num_rows; ++row) { - if (has_ranges && is_range[row]) { continue; } - const i_t direction = direction_of(row); - if (direction == 0) { continue; } - if (Arow.row_length(row) != 2) { continue; } - if (direction * user_problem_.rhs[row] != 0.0) { continue; } - - i_t gated = -1, activation = -1; - for (i_t p = Arow.row_start[row]; p < Arow.row_start[row + 1]; ++p) { - const i_t col = Arow.j[p]; - if (!is_binary(col)) { - gated = activation = -1; - break; - } - const f_t v = direction * Arow.x[p]; - if (v == 1.0) { - gated = col; - } else if (v == -1.0) { - activation = col; - } - } - if (gated >= 0 && activation >= 0) { gates.emplace_back(gated, activation); } - } - if (gates.empty()) { return; } - std::sort(gates.begin(), gates.end()); - gates.erase(std::unique(gates.begin(), gates.end()), gates.end()); - - auto gates_of = [&](i_t col) { - const auto lo = - std::lower_bound(gates.begin(), gates.end(), col, [](const std::pair& g, i_t c) { - return g.first < c; - }); - const auto hi = std::upper_bound( - lo, gates.end(), col, [](i_t c, const std::pair& g) { return c < g.first; }); - return std::make_pair(lo, hi); - }; + build_gate_table(Arow); + if (gate_activations_.empty()) { return; } // Pass two: the enabler rows, one +1 head against a tail of -1 selections. std::vector groups; std::vector>> candidates; for (i_t row = 0; row < num_rows; ++row) { if (has_ranges && is_range[row]) { continue; } - const i_t direction = direction_of(row); - if (direction == 0) { continue; } - const i_t len = Arow.row_length(row); + const char sense = user_problem_.row_sense[row]; + if (sense != 'L' && sense != 'G') { continue; } + // Row sense in <= orientation: +1 when the row reads a.x <= rhs, -1 when a.x >= rhs. + const i_t direction = sense == 'L' ? 1 : -1; + const i_t len = Arow.row_length(row); if (len < 3) { continue; } if (direction * user_problem_.rhs[row] != 0.0) { continue; } @@ -3333,20 +3364,20 @@ void cut_generation_t::build_group_cover_candidates( for (i_t p = Arow.row_start[row]; p < Arow.row_start[row + 1] && usable; ++p) { const i_t col = Arow.j[p]; const f_t v = direction * Arow.x[p]; - if (!is_binary(col)) { + if (user_problem_.var_types[col] == variable_type_t::CONTINUOUS || + user_problem_.lower[col] != 0.0 || user_problem_.upper[col] != 1.0) { usable = false; } else if (v == 1.0) { usable = head < 0; head = col; } else if (v == -1.0) { - auto [lo, hi] = gates_of(col); // A tail member with no gate leaves nothing to aggregate through; keep the candidate by // standing in the selection itself, which is what the enabler already bounds the head by. - if (lo == hi) { + if (gate_offsets_[col] == gate_offsets_[col + 1]) { groups.push_back(col); } else { - for (auto it = lo; it != hi; ++it) { - groups.push_back(it->second); + for (i_t p = gate_offsets_[col]; p < gate_offsets_[col + 1]; ++p) { + groups.push_back(gate_activations_[p]); } } } else { @@ -3378,8 +3409,9 @@ void cut_generation_t::build_group_cover_candidates( group_cover_offsets_.push_back(group_cover_groups_.size()); } - settings.log.print_format( - "Group cover: {} candidate cuts over {} gates\n", group_cover_heads_.size(), gates.size()); + settings.log.print_format("Group cover: {} candidate cuts over {} gates\n", + group_cover_heads_.size(), + gate_activations_.size()); } template @@ -3417,6 +3449,163 @@ void cut_generation_t::generate_group_cover_cuts( if (num_cuts > 0) { settings.log.debug("Generated %d group cover cuts\n", num_cuts); } } +// An activated capacity cut ties a group capacity row to the group's activation. Where the model +// carries +// +// sum_{i in S} x_i - s <= K (capacity) and x_i <= z for every i in S, +// +// z = 0 forces every x_i to zero while the relaxing term stays non-positive, so the row holds at a +// right-hand side of zero, and z = 1 leaves it unchanged: +// +// sum_{i in S} x_i - s <= K z. +// +// Valid for integral z only. Unlike the enabler rows behind a group cover cut, the capacity row +// stays in the model; the cut is the strengthened copy of it. +template +void cut_generation_t::build_activated_capacity_candidates( + const simplex_solver_settings_t& settings) +{ + activated_capacity_built_ = true; + + const i_t num_rows = user_problem_.num_rows; + const i_t num_cols = user_problem_.num_cols; + if (num_rows <= 0 || num_cols <= 0) { return; } + if (user_problem_.var_types.size() != (size_t)num_cols || + user_problem_.row_sense.size() != (size_t)num_rows || + user_problem_.rhs.size() != (size_t)num_rows) { + return; + } + + csr_matrix_t Arow(num_rows, num_cols, user_problem_.A.col_start[num_cols]); + user_problem_.A.to_compressed_row(Arow); + + const bool has_ranges = user_problem_.num_range_rows > 0; + std::vector is_range(has_ranges ? num_rows : 0, 0); + for (i_t k = 0; k < user_problem_.num_range_rows; ++k) { + is_range[user_problem_.range_rows[k]] = 1; + } + + build_gate_table(Arow); + if (gate_activations_.empty()) { return; } + + std::vector selections; + std::vector support; + std::vector common; + std::vector intersection; + + activated_capacity_offsets_.push_back(0); + for (i_t row = 0; row < num_rows; ++row) { + if (has_ranges && is_range[row]) { continue; } + const char sense = user_problem_.row_sense[row]; + if (sense != 'L' && sense != 'G') { continue; } + // Row sense in <= orientation: +1 when the row reads a.x <= rhs, -1 when a.x >= rhs. + const i_t direction = sense == 'L' ? 1 : -1; + if (Arow.row_length(row) < 2) { continue; } + + const f_t capacity = direction * user_problem_.rhs[row]; + if (capacity <= 0.0) { continue; } + + selections.clear(); + support.clear(); + bool usable = true; + for (i_t p = Arow.row_start[row]; p < Arow.row_start[row + 1] && usable; ++p) { + const i_t col = Arow.j[p]; + const f_t v = direction * Arow.x[p]; + support.push_back(col); + if (user_problem_.var_types[col] != variable_type_t::CONTINUOUS) { + usable = user_problem_.lower[col] == 0.0 && user_problem_.upper[col] == 1.0 && v == 1.0; + if (usable) { selections.push_back(col); } + } else { + // A relaxing term: non-positive over the whole box, so it cannot violate the row at z = 0. + usable = v < 0.0 && user_problem_.lower[col] >= 0.0; + } + } + if (!usable || selections.size() < 2) { continue; } + // At or above its own support size the cap is implied by the gates already, and so is K z. + const f_t n_selections = selections.size(); + if (capacity >= n_selections) { continue; } + + common.assign(gate_activations_.begin() + gate_offsets_[selections[0]], + gate_activations_.begin() + gate_offsets_[selections[0] + 1]); + for (size_t k = 1; k < selections.size() && !common.empty(); ++k) { + const i_t col = selections[k]; + intersection.clear(); + std::set_intersection(common.begin(), + common.end(), + gate_activations_.begin() + gate_offsets_[col], + gate_activations_.begin() + gate_offsets_[col + 1], + std::back_inserter(intersection)); + common.swap(intersection); + } + if (common.empty()) { continue; } + + std::sort(support.begin(), support.end()); + i_t activation = -1; + for (i_t z : common) { + if (std::binary_search(support.begin(), support.end(), z)) { continue; } + activation = z; + break; + } + if (activation < 0) { continue; } + + activated_capacity_activations_.push_back(activation); + activated_capacity_caps_.push_back(capacity); + for (i_t p = Arow.row_start[row]; p < Arow.row_start[row + 1]; ++p) { + activated_capacity_cols_.push_back(Arow.j[p]); + activated_capacity_coeffs_.push_back(direction * Arow.x[p]); + } + activated_capacity_offsets_.push_back(activated_capacity_cols_.size()); + } + if (activated_capacity_activations_.empty()) { return; } + + settings.log.print_format("Activated capacity: {} candidate cuts over {} gates\n", + activated_capacity_activations_.size(), + gate_activations_.size()); +} + +template +void cut_generation_t::generate_activated_capacity_cuts( + const simplex_solver_settings_t& settings, + const std::vector& xstar, + f_t start_time) +{ + if (!activated_capacity_built_) { build_activated_capacity_candidates(settings); } + if (activated_capacity_activations_.empty()) { return; } + + // The strengthened row carries a coefficient of K against +-1 everywhere else, so it costs the + // node LP far more per row than a unit-coefficient cut and is held to a higher bar than + // cut_pool_t's global min_cut_distance_ of 1e-4. Measured over the first root passes: the cuts + // that carry tsmc-setcover-3 have distance 0.23 upwards, while on tsmc-setcover-2, where the + // family buys no bound the group cover cuts do not already have, three quarters sit below 0.14. + const f_t min_distance = 0.2; + i_t num_cuts = 0; + const i_t n = activated_capacity_activations_.size(); + for (i_t k = 0; k < n; ++k) { + if ((k & 0xFF) == 0 && toc(start_time) >= settings.time_limit) { return; } + const i_t activation = activated_capacity_activations_[k]; + const f_t capacity = activated_capacity_caps_[k]; + f_t activity = -capacity * xstar[activation]; + f_t norm = capacity * capacity; + for (i_t p = activated_capacity_offsets_[k]; p < activated_capacity_offsets_[k + 1]; ++p) { + activity += activated_capacity_coeffs_[p] * xstar[activated_capacity_cols_[p]]; + norm += activated_capacity_coeffs_[p] * activated_capacity_coeffs_[p]; + } + if (activity <= min_distance * std::sqrt(norm)) { continue; } + + // add_cut expects cut'x >= rhs, so the row sum_p a_p x_p - K z <= 0 is emitted negated. + inequality_t cut; + for (i_t p = activated_capacity_offsets_[k]; p < activated_capacity_offsets_[k + 1]; ++p) { + cut.push_back(activated_capacity_cols_[p], -activated_capacity_coeffs_[p]); + } + cut.push_back(activation, capacity); + cut.rhs = 0.0; + cut_pool_.add_cut(cut_type_t::ACTIVATED_CAPACITY, cut); + num_cuts++; + } + + if (num_cuts > 0) { settings.log.debug("Generated %d activated capacity cuts\n", num_cuts); } +} + namespace { // Total probing-edge budget from the byte cap and the remaining work headroom @@ -3802,6 +3991,18 @@ bool cut_generation_t::generate_cuts(const lp_problem_t& lp, } } + // Generate activated capacity cuts + if (settings.activated_capacity_cuts != 0) { + if (toc(start_time) >= settings.time_limit) { return true; } + f_t cut_start_time = tic(); + generate_activated_capacity_cuts(settings, xstar, start_time); + f_t cut_generation_time = toc(cut_start_time); + if (cut_generation_time > 1.0) { + settings.log.debug("Activated capacity cut generation time %.2f seconds\n", + cut_generation_time); + } + } + // Build the fractional conflict-graph subgraph once (resolving the async // clique-table future on the way) so both clique-cut and zero-half cut // separators consume the same vertex/weight/adjacency tables instead of diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index a4cf3d592d..b3e382761f 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -46,7 +46,8 @@ enum cut_type_t : int8_t { ZERO_HALF = 6, FLOW_COVER = 7, GROUP_COVER = 8, - MAX_CUT_TYPE = 9 + ACTIVATED_CAPACITY = 9, + MAX_CUT_TYPE = 10 }; template @@ -188,7 +189,8 @@ struct cut_info_t { "Implied Bounds", "Zero-Half ", "Flow Cover ", - "Group Cover "}; + "Group Cover ", + "Activated Cap "}; std::array num_cuts = {0}; }; @@ -871,6 +873,20 @@ class cut_generation_t { // Scan the user problem for gate and enabler rows. Called once, on the first cut pass. void build_group_cover_candidates(const simplex::simplex_solver_settings_t& settings); + // Generate activated capacity cuts, tying a group capacity row to the group's activation + void generate_activated_capacity_cuts( + const simplex::simplex_solver_settings_t& settings, + const std::vector& xstar, + f_t start_time); + + // Scan the user problem for gate and capacity rows. Called once, on the first cut pass. + void build_activated_capacity_candidates( + const simplex::simplex_solver_settings_t& settings); + + // The variable-upper-bound gates x_j <= z, in CSR form over the columns. Shared by both + // structural separators and built on whichever of them runs first. + void build_gate_table(const csr_matrix_t& Arow); + void prepare_fractional_sub_conflict_graph( const simplex::simplex_solver_settings_t& settings, const std::vector& xstar, @@ -893,6 +909,19 @@ class cut_generation_t { std::vector group_cover_offsets_; std::vector group_cover_groups_; bool group_cover_built_{false}; + // The activations gating column j are gate_activations_[gate_offsets_[j] .. gate_offsets_[j+1]), + // sorted and deduplicated, so both separators can intersect the spans directly. + std::vector gate_offsets_; + std::vector gate_activations_; + bool gates_built_{false}; + // One candidate per capacity row whose members share an activation: that activation, the cap K, + // and the row itself in <= orientation, so separation needs no second pass over the matrix. + std::vector activated_capacity_activations_; + std::vector activated_capacity_caps_; + std::vector activated_capacity_offsets_; + std::vector activated_capacity_cols_; + std::vector activated_capacity_coeffs_; + bool activated_capacity_built_{false}; }; template diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index 64cb7f5726..227a9f1a36 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -99,6 +99,7 @@ struct simplex_solver_settings_t { flow_cover_cuts(-1), implied_bound_cuts(-1), group_cover_cuts(-1), + activated_capacity_cuts(-1), clique_cuts(-1), zero_half_cuts(-1), strong_chvatal_gomory_cuts(-1), @@ -212,6 +213,7 @@ struct simplex_solver_settings_t { i_t flow_cover_cuts; // -1 automatic, 0 to disable, >0 to enable flow cover cuts i_t implied_bound_cuts; // -1 automatic, 0 to disable, >0 to enable implied bound cuts i_t group_cover_cuts; // 0 to disable, >0 to enable group cover cuts + i_t activated_capacity_cuts; // 0 to disable, >0 to enable activated capacity cuts i_t clique_cuts; // -1 automatic, 0 to disable, >0 to enable clique cuts i_t zero_half_cuts; // -1 automatic, 0 to disable, >0 to enable zero-half cuts i_t strong_chvatal_gomory_cuts; // -1 automatic, 0 to disable, >0 to enable strong Chvatal Gomory diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index 1a8908d562..0a556ef0da 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -190,6 +190,7 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_MIP_ZERO_HALF_CUTS, &mip_settings.zero_half_cuts, -1, 1, -1}, {CUOPT_MIP_IMPLIED_BOUND_CUTS, &mip_settings.implied_bound_cuts, -1, 1, -1}, {CUOPT_MIP_GROUP_COVER_CUTS, &mip_settings.group_cover_cuts, -1, 1, -1}, + {CUOPT_MIP_ACTIVATED_CAPACITY_CUTS, &mip_settings.activated_capacity_cuts, -1, 1, -1}, {CUOPT_MIP_STRONG_CHVATAL_GOMORY_CUTS, &mip_settings.strong_chvatal_gomory_cuts, -1, 1, -1}, {CUOPT_MIP_REDUCED_COST_STRENGTHENING, &mip_settings.reduced_cost_strengthening, -1, std::numeric_limits::max(), -1}, {CUOPT_MIP_RINS, &mip_settings.submip_params.rins, -1, 1, -1}, diff --git a/cpp/src/mip_heuristics/CMakeLists.txt b/cpp/src/mip_heuristics/CMakeLists.txt index 4196697c04..187017fb14 100644 --- a/cpp/src/mip_heuristics/CMakeLists.txt +++ b/cpp/src/mip_heuristics/CMakeLists.txt @@ -15,7 +15,6 @@ set(MIP_LP_NECESSARY_FILES ${CMAKE_CURRENT_SOURCE_DIR}/presolve/single_lock_dual_aggregation.cpp ${CMAKE_CURRENT_SOURCE_DIR}/presolve/gf2_presolve.cpp ${CMAKE_CURRENT_SOURCE_DIR}/presolve/bhw_coeff_reduce.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/presolve/activated_capacity.cpp ${CMAKE_CURRENT_SOURCE_DIR}/solution/solution.cu ${CMAKE_CURRENT_SOURCE_DIR}/presolve/conflict_graph/clique_table.cu ) diff --git a/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh b/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh index 643bf694b4..c6f8453c9e 100644 --- a/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh +++ b/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh @@ -113,6 +113,7 @@ class sub_mip_recombiner_t : public recombiner_t { branch_and_bound_settings.clique_cuts = 0; branch_and_bound_settings.zero_half_cuts = 0; branch_and_bound_settings.group_cover_cuts = 0; + branch_and_bound_settings.activated_capacity_cuts = 0; branch_and_bound_settings.inside_submip = 1; branch_and_bound_settings.submip_settings.rins = 0; branch_and_bound_settings.submip_settings.rens = 0; diff --git a/cpp/src/mip_heuristics/presolve/activated_capacity.cpp b/cpp/src/mip_heuristics/presolve/activated_capacity.cpp deleted file mode 100644 index 9862330394..0000000000 --- a/cpp/src/mip_heuristics/presolve/activated_capacity.cpp +++ /dev/null @@ -1,225 +0,0 @@ -/* clang-format off */ -/* - * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. - * SPDX-License-Identifier: Apache-2.0 - */ -/* clang-format on */ - -#include "activated_capacity.hpp" - -#include -#include -#include - -#include -#include -#include -#include - -// A group capacity row caps how many members of a group may be selected, but says nothing about -// whether the group is open. Where every member is gated by the same activation z, the cap is only -// available once z is paid for: -// -// sum_{i in S} x_i - s <= K, x_i <= z for all i in S, s >= 0 -// => sum_{i in S} x_i - s <= K z. -// -// At z = 0 the gates force every x_i to 0 and the relaxing terms are non-positive, so the -// strengthened row reads 0 - s <= 0; at z = 1 it is the original row. Valid for integral z only, -// which is why this presolver is registered on the MIP path alone. - -namespace cuopt::mathematical_optimization::mip { - -template -papilo::PresolveStatus ActivatedCapacity::execute( - const papilo::Problem& problem, - const papilo::ProblemUpdate& problemUpdate, - const papilo::Num& num, - papilo::Reductions& reductions, - const papilo::Timer& timer, - int& reason_of_infeasibility) -{ - const auto& constraint_matrix = problem.getConstraintMatrix(); - const auto& lhs_values = constraint_matrix.getLeftHandSides(); - const auto& rhs_values = constraint_matrix.getRightHandSides(); - const auto& row_flags = constraint_matrix.getRowFlags(); - const auto& domains = problem.getVariableDomains(); - const auto& col_flags = domains.flags; - const auto& lower_bounds = domains.lower_bounds; - const auto& upper_bounds = domains.upper_bounds; - const auto& presolve_options = problemUpdate.getPresolveOptions(); - - const int num_rows = constraint_matrix.getNRows(); - - auto is_free_binary = [&](int col) { - const auto& flags = col_flags[col]; - return flags.test(papilo::ColFlag::kIntegral) && !flags.test(papilo::ColFlag::kLbInf) && - !flags.test(papilo::ColFlag::kUbInf) && !flags.test(papilo::ColFlag::kFixed) && - num.isZero(lower_bounds[col]) && num.isEq(upper_bounds[col], f_t{1}); - }; - - // Orientation of a one-sided row, or 0 when the row is an equation, a range or free. - // +1 means the stored row reads a.x <= side, -1 means a.x >= side; multiplying by the direction - // puts it in <= form either way. - auto orientation = [&](int row) { - const auto& row_flag = row_flags[row]; - if (row_flag.test(papilo::RowFlag::kRedundant)) return 0; - const bool lhs_infinite = row_flag.test(papilo::RowFlag::kLhsInf); - const bool rhs_infinite = row_flag.test(papilo::RowFlag::kRhsInf); - if (lhs_infinite == rhs_infinite) return 0; - return lhs_infinite ? 1 : -1; - }; - - // Pass one: every two-term row x - z <= 0 over free binaries, as sorted (x, z) pairs. - std::vector> gates; - for (int row = 0; row < num_rows; ++row) { - const int direction = orientation(row); - if (direction == 0) continue; - auto row_coefficients = constraint_matrix.getRowCoefficients(row); - if (row_coefficients.getLength() != 2) continue; - const f_t side = direction == 1 ? rhs_values[row] : lhs_values[row]; - if (!num.isZero(side)) continue; - - const int* indices = row_coefficients.getIndices(); - const f_t* values = row_coefficients.getValues(); - int gated = -1, activation = -1; - for (int j = 0; j < 2; ++j) { - if (!is_free_binary(indices[j])) break; - const f_t v = direction * values[j]; - if (num.isEq(v, f_t{1})) - gated = indices[j]; - else if (num.isEq(v, f_t{-1})) - activation = indices[j]; - } - if (gated >= 0 && activation >= 0) gates.emplace_back(gated, activation); - } - std::sort(gates.begin(), gates.end()); - gates.erase(std::unique(gates.begin(), gates.end()), gates.end()); - - auto gates_of = [&](int col) { - const auto lo = - std::lower_bound(gates.begin(), gates.end(), col, [](const std::pair& g, int c) { - return g.first < c; - }); - const auto hi = std::upper_bound( - lo, gates.end(), col, [](int c, const std::pair& g) { return c < g.first; }); - return std::make_pair(lo, hi); - }; - - papilo::PresolveStatus status = papilo::PresolveStatus::kUnchanged; - int rows_strengthened = 0; - std::vector selections; - std::vector common; - std::vector candidates; - std::vector intersection; - std::vector support; - - // Pass two: the capacity rows themselves. - for (int row = 0; row < num_rows && !gates.empty(); ++row) { - if (reductions.size() >= presolve_options.max_reduction_seq) break; - if (papilo::PresolveMethod::is_interrupted( - timer, presolve_options.tlim, presolve_options.early_exit_callback)) - break; - - const int direction = orientation(row); - if (direction == 0) continue; - auto row_coefficients = constraint_matrix.getRowCoefficients(row); - const int len = row_coefficients.getLength(); - if (len < 2 || len > ACTIVATED_CAPACITY_MAX_LEN) continue; - - const f_t side = direction == 1 ? rhs_values[row] : lhs_values[row]; - const f_t capacity = direction * side; - if (!num.isGT(capacity, f_t{0})) continue; - - const int* indices = row_coefficients.getIndices(); - const f_t* values = row_coefficients.getValues(); - selections.clear(); - support.assign(indices, indices + len); - bool usable = true; - for (int j = 0; j < len && usable; ++j) { - const int col = indices[j]; - const f_t v = direction * values[j]; - if (col_flags[col].test(papilo::ColFlag::kIntegral)) { - // A negative coefficient on an integral column is what this presolver itself writes, so - // rejecting it here is also what keeps a rewritten row from matching a second time. - usable = is_free_binary(col) && num.isEq(v, f_t{1}); - if (usable) selections.push_back(col); - } else { - // A relaxing term: non-positive over the whole box, so it cannot violate the row at z = 0. - usable = num.isLT(v, f_t{0}) && !col_flags[col].test(papilo::ColFlag::kLbInf) && - !num.isLT(lower_bounds[col], f_t{0}); - } - } - if (!usable || selections.size() < 2) continue; - // At or above its own support size the cap is implied by the gates already, and so is K z. - const f_t n_selections = selections.size(); - if (!num.isLT(capacity, n_selections)) continue; - - auto [lo, hi] = gates_of(selections[0]); - common.clear(); - for (auto it = lo; it != hi; ++it) - common.push_back(it->second); - std::sort(common.begin(), common.end()); - for (size_t k = 1; k < selections.size() && !common.empty(); ++k) { - auto [klo, khi] = gates_of(selections[k]); - candidates.clear(); - for (auto it = klo; it != khi; ++it) - candidates.push_back(it->second); - std::sort(candidates.begin(), candidates.end()); - intersection.clear(); - std::set_intersection(common.begin(), - common.end(), - candidates.begin(), - candidates.end(), - std::back_inserter(intersection)); - common.swap(intersection); - } - if (common.empty()) continue; - - std::sort(support.begin(), support.end()); - int activation = -1; - for (int z : common) { - if (std::binary_search(support.begin(), support.end(), z)) continue; - activation = z; - break; - } - if (activation < 0) continue; - - cuopt_assert(is_free_binary(activation), "the activation of a gate row is a free binary"); - - papilo::TransactionGuard guard{reductions}; - reductions.lockRow(row); - reductions.changeMatrixEntry(row, activation, direction * -capacity); - if (direction == 1) - reductions.changeRowRHS(row, f_t{0}); - else - reductions.changeRowLHS(row, f_t{0}); - ++rows_strengthened; - status = papilo::PresolveStatus::kReduced; - } - - // Proposed, not applied: the activation is a nonzero the row does not have yet, and - // ConstraintMatrix::change_coefficient refuses the insert when the row or the column has no slack - // space left in its range, which rejects the whole transaction. The presolved nonzero count is - // the figure to check this against. - if (rows_strengthened > 0) { - CUOPT_LOG_INFO("Activated capacity: proposed %d strengthened rows against %zu gates", - rows_strengthened, - gates.size()); - } - - return status; -} - -#define INSTANTIATE(F_TYPE) template class ActivatedCapacity; - -#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT -INSTANTIATE(float) -#endif - -#if MIP_INSTANTIATE_DOUBLE -INSTANTIATE(double) -#endif - -#undef INSTANTIATE - -} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/presolve/activated_capacity.hpp b/cpp/src/mip_heuristics/presolve/activated_capacity.hpp deleted file mode 100644 index 6d1f3e8e35..0000000000 --- a/cpp/src/mip_heuristics/presolve/activated_capacity.hpp +++ /dev/null @@ -1,51 +0,0 @@ -/* clang-format off */ -/* - * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. - * SPDX-License-Identifier: Apache-2.0 - */ -/* clang-format on */ - -#pragma once - -#if !defined(__clang__) -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wstringop-overflow" // ignore boost error for pip wheel build -#pragma GCC diagnostic ignored "-Wnarrowing" -#endif -#include -#include -#include -#include -#if !defined(__clang__) -#pragma GCC diagnostic pop -#endif - -#include - -namespace cuopt::mathematical_optimization::mip { - -// A capacity row wider than this is not worth the gate-set intersection. -static constexpr int ACTIVATED_CAPACITY_MAX_LEN = 4096; - -template -class ActivatedCapacity : public papilo::PresolveMethod { - public: - ActivatedCapacity() : papilo::PresolveMethod() - { - this->setName("activatedcapacity"); - this->setType(papilo::PresolverType::kIntegralCols); - this->setTiming(papilo::PresolverTiming::kMedium); - // The vacuous copies of these capacity rows are dropped by the cheaper presolvers first, which - // keeps the gate-set intersection off rows that would gain nothing. - this->setDelayed(true); - } - - papilo::PresolveStatus execute(const papilo::Problem& problem, - const papilo::ProblemUpdate& problemUpdate, - const papilo::Num& num, - papilo::Reductions& reductions, - const papilo::Timer& timer, - int& reason_of_infeasibility) override; -}; - -} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp index b3f35e0b40..baf311889c 100644 --- a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp +++ b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp @@ -40,7 +40,6 @@ #include #include #include -#include #include #include #include @@ -748,7 +747,6 @@ void set_presolve_methods( // cuOpt custom GF2 presolver maybe_add(uptr(new GF2Presolve())); maybe_add(uptr(new BHWCoeffReduce())); - maybe_add(uptr(new ActivatedCapacity())); } // fast presolvers maybe_add(uptr(new papilo::SingletonCols())); diff --git a/cpp/src/mip_heuristics/solver.cu b/cpp/src/mip_heuristics/solver.cu index 535163d53c..c5675db075 100644 --- a/cpp/src/mip_heuristics/solver.cu +++ b/cpp/src/mip_heuristics/solver.cu @@ -366,12 +366,13 @@ solution_t mip_solver_t::run_solver() } branch_and_bound_settings.mixed_integer_gomory_cuts = context.settings.mixed_integer_gomory_cuts; - branch_and_bound_settings.knapsack_cuts = context.settings.knapsack_cuts; - branch_and_bound_settings.flow_cover_cuts = context.settings.flow_cover_cuts; - branch_and_bound_settings.implied_bound_cuts = context.settings.implied_bound_cuts; - branch_and_bound_settings.group_cover_cuts = context.settings.group_cover_cuts; - branch_and_bound_settings.clique_cuts = context.settings.clique_cuts; - branch_and_bound_settings.zero_half_cuts = context.settings.zero_half_cuts; + branch_and_bound_settings.knapsack_cuts = context.settings.knapsack_cuts; + branch_and_bound_settings.flow_cover_cuts = context.settings.flow_cover_cuts; + branch_and_bound_settings.implied_bound_cuts = context.settings.implied_bound_cuts; + branch_and_bound_settings.group_cover_cuts = context.settings.group_cover_cuts; + branch_and_bound_settings.activated_capacity_cuts = context.settings.activated_capacity_cuts; + branch_and_bound_settings.clique_cuts = context.settings.clique_cuts; + branch_and_bound_settings.zero_half_cuts = context.settings.zero_half_cuts; branch_and_bound_settings.strong_chvatal_gomory_cuts = context.settings.strong_chvatal_gomory_cuts; branch_and_bound_settings.cut_change_threshold = context.settings.cut_change_threshold; From d3a7a2f4aeb64571d569fd5c618d0db741d22f73 Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Wed, 23 Sep 2026 16:42:54 +0200 Subject: [PATCH 05/12] renamed variables/cuts/etc. to be closer to the textbook nomenclature. Signed-off-by: Nicolas L. Guidotti --- .../mathematical_optimization/constants.h | 4 +- .../mip/solver_settings.hpp | 4 +- cpp/src/branch_and_bound/branch_and_bound.cpp | 4 +- cpp/src/cuts/cuts.cpp | 308 +++++++++--------- cpp/src/cuts/cuts.hpp | 82 ++--- .../dual_simplex/simplex_solver_settings.hpp | 8 +- cpp/src/math_optimization/solver_settings.cu | 4 +- .../diversity/recombiners/sub_mip.cuh | 4 +- cpp/src/mip_heuristics/solver.cu | 14 +- 9 files changed, 219 insertions(+), 213 deletions(-) diff --git a/cpp/include/cuopt/mathematical_optimization/constants.h b/cpp/include/cuopt/mathematical_optimization/constants.h index 888eed9a15..7ddc5eaaf8 100644 --- a/cpp/include/cuopt/mathematical_optimization/constants.h +++ b/cpp/include/cuopt/mathematical_optimization/constants.h @@ -80,8 +80,8 @@ #define CUOPT_MIP_CLIQUE_CUTS "mip_clique_cuts" #define CUOPT_MIP_ZERO_HALF_CUTS "mip_zero_half_cuts" #define CUOPT_MIP_STRONG_CHVATAL_GOMORY_CUTS "mip_strong_chvatal_gomory_cuts" -#define CUOPT_MIP_GROUP_COVER_CUTS "mip_group_cover_cuts" -#define CUOPT_MIP_ACTIVATED_CAPACITY_CUTS "mip_activated_capacity_cuts" +#define CUOPT_MIP_IMPLIED_INDICATOR_CUTS "mip_implied_indicator_cuts" +#define CUOPT_MIP_CAPACITY_LIFTING_CUTS "mip_capacity_lifting_cuts" #define CUOPT_MIP_REDUCED_COST_STRENGTHENING "mip_reduced_cost_strengthening" #define CUOPT_MIP_RINS "mip_rins" #define CUOPT_MIP_RENS "mip_rens" diff --git a/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp b/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp index e177623549..0d28258164 100644 --- a/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp +++ b/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp @@ -135,8 +135,8 @@ class mip_solver_settings_t { i_t clique_cuts = -1; i_t zero_half_cuts = -1; i_t implied_bound_cuts = -1; - i_t group_cover_cuts = -1; - i_t activated_capacity_cuts = -1; + i_t implied_indicator_cuts = -1; + i_t capacity_lifting_cuts = -1; i_t strong_chvatal_gomory_cuts = -1; i_t reduced_cost_strengthening = -1; i_t objective_step = 1; // 0 = disable objective step tightening, 1 = enable diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 0e778d2eea..6a170154ae 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -2321,8 +2321,8 @@ void branch_and_bound_t::solve_submip(diving_worker_t* worke submip_settings.reliability_branching = 0; submip_settings.clique_cuts = 0; submip_settings.zero_half_cuts = 0; - submip_settings.group_cover_cuts = 0; - submip_settings.activated_capacity_cuts = 0; + submip_settings.implied_indicator_cuts = 0; + submip_settings.capacity_lifting_cuts = 0; submip_settings.inside_submip = 1; submip_settings.strong_branching_simplex_iteration_limit = 50; submip_settings.inside_root_node = 0; diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index 26c27a2af4..986ad4510d 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -3232,11 +3232,16 @@ void cut_generation_t::generate_implied_bound_cuts( } } +// The variable upper bounds x_j <= u_j z linking a member to the indicator that must be on before +// it can be chosen. Only the binary unit form x_j <= z is recognised here -- both columns binary, +// coefficients +-1, right-hand side zero -- which is narrower than either separator needs: the +// implied indicator cut holds for any u_j > 0, and the capacity lifting cut does not need the +// member to be integral at all. Widening the scan would find more structure than it does today. template -void cut_generation_t::build_gate_table(const csr_matrix_t& Arow) +void cut_generation_t::build_vub_table(const csr_matrix_t& Arow) { - if (gates_built_) { return; } - gates_built_ = true; + if (vub_built_) { return; } + vub_built_ = true; const i_t num_rows = user_problem_.num_rows; const i_t num_cols = user_problem_.num_cols; @@ -3247,7 +3252,7 @@ void cut_generation_t::build_gate_table(const csr_matrix_t& is_range[user_problem_.range_rows[k]] = 1; } - std::vector> gates; + std::vector> vubs; for (i_t row = 0; row < num_rows; ++row) { if (has_ranges && is_range[row]) { continue; } if (Arow.row_length(row) != 2) { continue; } @@ -3257,72 +3262,73 @@ void cut_generation_t::build_gate_table(const csr_matrix_t& const i_t direction = sense == 'L' ? 1 : -1; if (direction * user_problem_.rhs[row] != 0.0) { continue; } - i_t gated = -1, activation = -1; + i_t member = -1, indicator = -1; for (i_t p = Arow.row_start[row]; p < Arow.row_start[row + 1]; ++p) { const i_t col = Arow.j[p]; if (user_problem_.var_types[col] == variable_type_t::CONTINUOUS || user_problem_.lower[col] != 0.0 || user_problem_.upper[col] != 1.0) { - gated = activation = -1; + member = indicator = -1; break; } const f_t v = direction * Arow.x[p]; if (v == 1.0) { - gated = col; + member = col; } else if (v == -1.0) { - activation = col; + indicator = col; } } - if (gated >= 0 && activation >= 0) { gates.emplace_back(gated, activation); } + if (member >= 0 && indicator >= 0) { vubs.emplace_back(member, indicator); } } - if (gates.empty()) { return; } + if (vubs.empty()) { return; } - gate_offsets_.assign(num_cols + 1, 0); - for (const auto& [gated, activation] : gates) { - ++gate_offsets_[gated + 1]; + vub_offsets_.assign(num_cols + 1, 0); + for (const auto& [member, indicator] : vubs) { + ++vub_offsets_[member + 1]; } for (i_t col = 0; col < num_cols; ++col) { - gate_offsets_[col + 1] += gate_offsets_[col]; + vub_offsets_[col + 1] += vub_offsets_[col]; } - gate_activations_.resize(gates.size()); - std::vector cursor(gate_offsets_.begin(), gate_offsets_.end() - 1); - for (const auto& [gated, activation] : gates) { - gate_activations_[cursor[gated]++] = activation; + vub_indicators_.resize(vubs.size()); + std::vector cursor(vub_offsets_.begin(), vub_offsets_.end() - 1); + for (const auto& [member, indicator] : vubs) { + vub_indicators_[cursor[member]++] = indicator; } - // Sort each column's span and drop duplicate gate rows, compacting in place. The write cursor - // trails the read cursor because it only advances on a kept entry. + // Sort each column's span and drop rows stating the same bound twice, compacting in place. The + // write cursor trails the read cursor because it only advances on a kept entry. i_t out = 0; for (i_t col = 0; col < num_cols; ++col) { - const i_t begin = gate_offsets_[col]; - const i_t end = gate_offsets_[col + 1]; - gate_offsets_[col] = out; - std::sort(gate_activations_.begin() + begin, gate_activations_.begin() + end); + const i_t begin = vub_offsets_[col]; + const i_t end = vub_offsets_[col + 1]; + vub_offsets_[col] = out; + std::sort(vub_indicators_.begin() + begin, vub_indicators_.begin() + end); for (i_t p = begin; p < end; ++p) { - if (p == begin || gate_activations_[p] != gate_activations_[p - 1]) { - gate_activations_[out++] = gate_activations_[p]; + if (p == begin || vub_indicators_[p] != vub_indicators_[p - 1]) { + vub_indicators_[out++] = vub_indicators_[p]; } } } - gate_offsets_[num_cols] = out; - gate_activations_.resize(out); + vub_offsets_[num_cols] = out; + vub_indicators_.resize(out); } -// A group cover cut aggregates an enabler row through the gates behind it. Where the model carries +// An implied indicator cut aggregates an implication row over the indicators of its members. Where +// the model carries // -// y <= sum_{j in S} x_j (enabler) and x_j <= z_{g(j)} for every j in S, +// y <= sum_{j in S} x_j (implication) and x_j <= z_{g(j)} for every j in S, // -// a binary y that is one forces some x_j to one, which forces its own group activation to one, so +// a binary y that is one forces some single x_j to one, which forces that member's indicator to +// one, so // -// y <= sum_{g in D} z_g, D = the distinct groups covering S. +// y <= sum_{g in D} z_g, D = the distinct indicators over S. // -// Counting each group once is where the strength is: the enabler alone lets y reach one against a -// whole group held at 1/|S|, while the cut holds y down to that group's own activation. Valid for -// integral y and x only. +// Counting each indicator once is where the strength is: chaining the members through their own +// variable upper bounds gives y <= sum_j z_{g(j)}, which counts one indicator once per member. template -void cut_generation_t::build_group_cover_candidates( +void cut_generation_t::build_implied_indicator_candidates( const simplex_solver_settings_t& settings) { - group_cover_built_ = true; + implied_indicator_built_ = true; const i_t num_rows = user_problem_.num_rows; const i_t num_cols = user_problem_.num_cols; @@ -3342,12 +3348,12 @@ void cut_generation_t::build_group_cover_candidates( is_range[user_problem_.range_rows[k]] = 1; } - build_gate_table(Arow); - if (gate_activations_.empty()) { return; } + build_vub_table(Arow); + if (vub_indicators_.empty()) { return; } - // Pass two: the enabler rows, one +1 head against a tail of -1 selections. - std::vector groups; - std::vector>> candidates; + // The implication rows, one +1 head against a tail of -1 members. + std::vector indicators; + implied_indicator_offsets_.push_back(0); for (i_t row = 0; row < num_rows; ++row) { if (has_ranges && is_range[row]) { continue; } const char sense = user_problem_.row_sense[row]; @@ -3360,10 +3366,10 @@ void cut_generation_t::build_group_cover_candidates( i_t head = -1; bool usable = true; - groups.clear(); - for (i_t p = Arow.row_start[row]; p < Arow.row_start[row + 1] && usable; ++p) { - const i_t col = Arow.j[p]; - const f_t v = direction * Arow.x[p]; + indicators.clear(); + for (i_t k = Arow.row_start[row]; k < Arow.row_start[row + 1] && usable; ++k) { + const i_t col = Arow.j[k]; + const f_t v = direction * Arow.x[k]; if (user_problem_.var_types[col] == variable_type_t::CONTINUOUS || user_problem_.lower[col] != 0.0 || user_problem_.upper[col] != 1.0) { usable = false; @@ -3371,101 +3377,96 @@ void cut_generation_t::build_group_cover_candidates( usable = head < 0; head = col; } else if (v == -1.0) { - // A tail member with no gate leaves nothing to aggregate through; keep the candidate by - // standing in the selection itself, which is what the enabler already bounds the head by. - if (gate_offsets_[col] == gate_offsets_[col + 1]) { - groups.push_back(col); + // A member with no variable upper bound leaves nothing to aggregate through; keep the + // candidate by standing in the member itself, which the implication row already bounds + // the head by. + if (vub_offsets_[col] == vub_offsets_[col + 1]) { + indicators.push_back(col); } else { - for (i_t p = gate_offsets_[col]; p < gate_offsets_[col + 1]; ++p) { - groups.push_back(gate_activations_[p]); + for (i_t p = vub_offsets_[col]; p < vub_offsets_[col + 1]; ++p) { + indicators.push_back(vub_indicators_[p]); } } } else { usable = false; } } - if (!usable || head < 0 || groups.empty()) { continue; } - - std::sort(groups.begin(), groups.end()); - groups.erase(std::unique(groups.begin(), groups.end()), groups.end()); - // Nothing was merged, so the cut is the sum of the gates the LP already has. - if (groups.size() >= (size_t)(len - 1)) { continue; } - if (std::binary_search(groups.begin(), groups.end(), head)) { continue; } + if (!usable || head < 0 || indicators.empty()) { continue; } - candidates.emplace_back(head, groups); - } - if (candidates.empty()) { return; } + std::sort(indicators.begin(), indicators.end()); + indicators.erase(std::unique(indicators.begin(), indicators.end()), indicators.end()); + // Nothing was merged, so the cut is the sum of the variable upper bounds the LP already has. + if (indicators.size() >= (size_t)(len - 1)) { continue; } + if (std::binary_search(indicators.begin(), indicators.end(), head)) { continue; } - // Two enabler rows over the same group set give the same cut; keep one. - std::sort(candidates.begin(), candidates.end()); - candidates.erase(std::unique(candidates.begin(), candidates.end()), candidates.end()); - - group_cover_heads_.reserve(candidates.size()); - group_cover_offsets_.reserve(candidates.size() + 1); - group_cover_offsets_.push_back(0); - for (const auto& [head, group_set] : candidates) { - group_cover_heads_.push_back(head); - group_cover_groups_.insert(group_cover_groups_.end(), group_set.begin(), group_set.end()); - group_cover_offsets_.push_back(group_cover_groups_.size()); + implied_indicator_heads_.push_back(head); + implied_indicator_indicators_.insert( + implied_indicator_indicators_.end(), indicators.begin(), indicators.end()); + implied_indicator_offsets_.push_back(implied_indicator_indicators_.size()); } + if (implied_indicator_heads_.empty()) { return; } - settings.log.print_format("Group cover: {} candidate cuts over {} gates\n", - group_cover_heads_.size(), - gate_activations_.size()); + settings.log.print_format("Implied indicator: {} candidate cuts over {} variable upper bounds\n", + implied_indicator_heads_.size(), + vub_indicators_.size()); } template -void cut_generation_t::generate_group_cover_cuts( +void cut_generation_t::generate_implied_indicator_cuts( const simplex_solver_settings_t& settings, const std::vector& xstar, f_t start_time) { - if (!group_cover_built_) { build_group_cover_candidates(settings); } - if (group_cover_heads_.empty()) { return; } + if (!implied_indicator_built_) { build_implied_indicator_candidates(settings); } + if (implied_indicator_heads_.empty()) { return; } const f_t tol = 1e-4; i_t num_cuts = 0; - const i_t n = group_cover_heads_.size(); + const i_t n = implied_indicator_heads_.size(); for (i_t k = 0; k < n; ++k) { if ((k & 0xFF) == 0 && toc(start_time) >= settings.time_limit) { return; } - const i_t head = group_cover_heads_[k]; + const i_t head = implied_indicator_heads_[k]; f_t activity = xstar[head]; - for (i_t p = group_cover_offsets_[k]; p < group_cover_offsets_[k + 1]; ++p) { - activity -= xstar[group_cover_groups_[p]]; + for (i_t p = implied_indicator_offsets_[k]; p < implied_indicator_offsets_[k + 1]; ++p) { + activity -= xstar[implied_indicator_indicators_[p]]; } if (activity <= tol) { continue; } // add_cut expects cut'x >= rhs, so the cut head - sum_g z_g <= 0 is emitted negated. inequality_t cut; cut.push_back(head, -1.0); - for (i_t p = group_cover_offsets_[k]; p < group_cover_offsets_[k + 1]; ++p) { - cut.push_back(group_cover_groups_[p], 1.0); + for (i_t p = implied_indicator_offsets_[k]; p < implied_indicator_offsets_[k + 1]; ++p) { + cut.push_back(implied_indicator_indicators_[p], 1.0); } cut.rhs = 0.0; - cut_pool_.add_cut(cut_type_t::GROUP_COVER, cut); + cut_pool_.add_cut(cut_type_t::IMPLIED_INDICATOR, cut); num_cuts++; } - if (num_cuts > 0) { settings.log.debug("Generated %d group cover cuts\n", num_cuts); } + if (num_cuts > 0) { settings.log.debug("Generated %d implied indicator cuts\n", num_cuts); } } -// An activated capacity cut ties a group capacity row to the group's activation. Where the model +// A capacity lifting cut lifts the complemented indicator into a capacity row. Where the model // carries // // sum_{i in S} x_i - s <= K (capacity) and x_i <= z for every i in S, // -// z = 0 forces every x_i to zero while the relaxing term stays non-positive, so the row holds at a -// right-hand side of zero, and z = 1 leaves it unchanged: +// sequential lifting of zbar = 1 - z into the capacity row asks for the largest coefficient a +// keeping sum_i x_i - s + a zbar <= K valid. At zbar = 1 the variable upper bounds force every +// x_i to zero and -s is non-positive, so the left side maxes out at 0 and a may rise to K: // // sum_{i in S} x_i - s <= K z. // -// Valid for integral z only. Unlike the enabler rows behind a group cover cut, the capacity row -// stays in the model; the cut is the strengthened copy of it. +// K is therefore the maximal lifting coefficient, not a choice. Valid for integral z only. Unlike +// the implication rows behind an implied indicator cut, the capacity row stays in the model; the +// cut is the lifted copy of it. The capacity row also bounds the violation, since it already gives +// sum_i x*_i - s* <= K and hence a violation of at most K(1 - z*): the cut can only bite where the +// row is near-tight and the indicator is well below one. template -void cut_generation_t::build_activated_capacity_candidates( +void cut_generation_t::build_capacity_lifting_candidates( const simplex_solver_settings_t& settings) { - activated_capacity_built_ = true; + capacity_lifting_built_ = true; const i_t num_rows = user_problem_.num_rows; const i_t num_cols = user_problem_.num_cols; @@ -3485,15 +3486,15 @@ void cut_generation_t::build_activated_capacity_candidates( is_range[user_problem_.range_rows[k]] = 1; } - build_gate_table(Arow); - if (gate_activations_.empty()) { return; } + build_vub_table(Arow); + if (vub_indicators_.empty()) { return; } - std::vector selections; + std::vector members; std::vector support; std::vector common; std::vector intersection; - activated_capacity_offsets_.push_back(0); + capacity_lifting_offsets_.push_back(0); for (i_t row = 0; row < num_rows; ++row) { if (has_ranges && is_range[row]) { continue; } const char sense = user_problem_.row_sense[row]; @@ -3505,7 +3506,7 @@ void cut_generation_t::build_activated_capacity_candidates( const f_t capacity = direction * user_problem_.rhs[row]; if (capacity <= 0.0) { continue; } - selections.clear(); + members.clear(); support.clear(); bool usable = true; for (i_t p = Arow.row_start[row]; p < Arow.row_start[row + 1] && usable; ++p) { @@ -3514,96 +3515,98 @@ void cut_generation_t::build_activated_capacity_candidates( support.push_back(col); if (user_problem_.var_types[col] != variable_type_t::CONTINUOUS) { usable = user_problem_.lower[col] == 0.0 && user_problem_.upper[col] == 1.0 && v == 1.0; - if (usable) { selections.push_back(col); } + if (usable) { members.push_back(col); } } else { // A relaxing term: non-positive over the whole box, so it cannot violate the row at z = 0. usable = v < 0.0 && user_problem_.lower[col] >= 0.0; } } - if (!usable || selections.size() < 2) { continue; } - // At or above its own support size the cap is implied by the gates already, and so is K z. - const f_t n_selections = selections.size(); - if (capacity >= n_selections) { continue; } + if (!usable || members.size() < 2) { continue; } + // At or above its own support size the cap is implied by the variable upper bounds, and so is + // K z, leaving nothing to lift. + const f_t n_members = members.size(); + if (capacity >= n_members) { continue; } - common.assign(gate_activations_.begin() + gate_offsets_[selections[0]], - gate_activations_.begin() + gate_offsets_[selections[0] + 1]); - for (size_t k = 1; k < selections.size() && !common.empty(); ++k) { - const i_t col = selections[k]; + common.assign(vub_indicators_.begin() + vub_offsets_[members[0]], + vub_indicators_.begin() + vub_offsets_[members[0] + 1]); + for (size_t k = 1; k < members.size() && !common.empty(); ++k) { + const i_t col = members[k]; intersection.clear(); std::set_intersection(common.begin(), common.end(), - gate_activations_.begin() + gate_offsets_[col], - gate_activations_.begin() + gate_offsets_[col + 1], + vub_indicators_.begin() + vub_offsets_[col], + vub_indicators_.begin() + vub_offsets_[col + 1], std::back_inserter(intersection)); common.swap(intersection); } if (common.empty()) { continue; } std::sort(support.begin(), support.end()); - i_t activation = -1; + i_t indicator = -1; for (i_t z : common) { if (std::binary_search(support.begin(), support.end(), z)) { continue; } - activation = z; + indicator = z; break; } - if (activation < 0) { continue; } + if (indicator < 0) { continue; } - activated_capacity_activations_.push_back(activation); - activated_capacity_caps_.push_back(capacity); + capacity_lifting_indicators_.push_back(indicator); + capacity_lifting_caps_.push_back(capacity); for (i_t p = Arow.row_start[row]; p < Arow.row_start[row + 1]; ++p) { - activated_capacity_cols_.push_back(Arow.j[p]); - activated_capacity_coeffs_.push_back(direction * Arow.x[p]); + capacity_lifting_cols_.push_back(Arow.j[p]); + capacity_lifting_coeffs_.push_back(direction * Arow.x[p]); } - activated_capacity_offsets_.push_back(activated_capacity_cols_.size()); + capacity_lifting_offsets_.push_back(capacity_lifting_cols_.size()); } - if (activated_capacity_activations_.empty()) { return; } + if (capacity_lifting_indicators_.empty()) { return; } - settings.log.print_format("Activated capacity: {} candidate cuts over {} gates\n", - activated_capacity_activations_.size(), - gate_activations_.size()); + settings.log.print_format("Capacity lifting: {} candidate cuts over {} variable upper bounds\n", + capacity_lifting_indicators_.size(), + vub_indicators_.size()); } template -void cut_generation_t::generate_activated_capacity_cuts( +void cut_generation_t::generate_capacity_lifting_cuts( const simplex_solver_settings_t& settings, const std::vector& xstar, f_t start_time) { - if (!activated_capacity_built_) { build_activated_capacity_candidates(settings); } - if (activated_capacity_activations_.empty()) { return; } - - // The strengthened row carries a coefficient of K against +-1 everywhere else, so it costs the - // node LP far more per row than a unit-coefficient cut and is held to a higher bar than - // cut_pool_t's global min_cut_distance_ of 1e-4. Measured over the first root passes: the cuts - // that carry tsmc-setcover-3 have distance 0.23 upwards, while on tsmc-setcover-2, where the - // family buys no bound the group cover cuts do not already have, three quarters sit below 0.14. + if (!capacity_lifting_built_) { build_capacity_lifting_candidates(settings); } + if (capacity_lifting_indicators_.empty()) { return; } + + // The lifted row carries a coefficient of K against +-1 everywhere else, so it costs the node LP + // far more per row than a unit-coefficient cut and is held to a higher bar than cut_pool_t's + // global min_cut_distance_ of 1e-4. Measured over the first root passes: the cuts that carry + // tsmc-setcover-3, where the cap is 3 of 25 members and the rows bind constantly, have distance + // 0.23 upwards; on tsmc-setcover-2 the cap is 12 of 25, the rows rarely bind, and three quarters + // sit below 0.14 while buying no bound the implied indicator cuts do not already have. const f_t min_distance = 0.2; i_t num_cuts = 0; - const i_t n = activated_capacity_activations_.size(); + const i_t n = capacity_lifting_indicators_.size(); for (i_t k = 0; k < n; ++k) { if ((k & 0xFF) == 0 && toc(start_time) >= settings.time_limit) { return; } - const i_t activation = activated_capacity_activations_[k]; - const f_t capacity = activated_capacity_caps_[k]; - f_t activity = -capacity * xstar[activation]; - f_t norm = capacity * capacity; - for (i_t p = activated_capacity_offsets_[k]; p < activated_capacity_offsets_[k + 1]; ++p) { - activity += activated_capacity_coeffs_[p] * xstar[activated_capacity_cols_[p]]; - norm += activated_capacity_coeffs_[p] * activated_capacity_coeffs_[p]; + const i_t indicator = capacity_lifting_indicators_[k]; + const f_t capacity = capacity_lifting_caps_[k]; + f_t activity = -capacity * xstar[indicator]; + f_t norm = capacity * capacity; + for (i_t p = capacity_lifting_offsets_[k]; p < capacity_lifting_offsets_[k + 1]; ++p) { + activity += capacity_lifting_coeffs_[p] * xstar[capacity_lifting_cols_[p]]; + norm += capacity_lifting_coeffs_[p] * capacity_lifting_coeffs_[p]; } if (activity <= min_distance * std::sqrt(norm)) { continue; } // add_cut expects cut'x >= rhs, so the row sum_p a_p x_p - K z <= 0 is emitted negated. inequality_t cut; - for (i_t p = activated_capacity_offsets_[k]; p < activated_capacity_offsets_[k + 1]; ++p) { - cut.push_back(activated_capacity_cols_[p], -activated_capacity_coeffs_[p]); + for (i_t p = capacity_lifting_offsets_[k]; p < capacity_lifting_offsets_[k + 1]; ++p) { + cut.push_back(capacity_lifting_cols_[p], -capacity_lifting_coeffs_[p]); } - cut.push_back(activation, capacity); + cut.push_back(indicator, capacity); cut.rhs = 0.0; - cut_pool_.add_cut(cut_type_t::ACTIVATED_CAPACITY, cut); + cut_pool_.add_cut(cut_type_t::CAPACITY_LIFTING, cut); num_cuts++; } - if (num_cuts > 0) { settings.log.debug("Generated %d activated capacity cuts\n", num_cuts); } + if (num_cuts > 0) { settings.log.debug("Generated %d capacity lifting cuts\n", num_cuts); } } namespace { @@ -3980,25 +3983,26 @@ bool cut_generation_t::generate_cuts(const lp_problem_t& lp, } } - // Generate group cover cuts - if (settings.group_cover_cuts != 0) { + // Generate implied indicator cuts + if (settings.implied_indicator_cuts != 0) { if (toc(start_time) >= settings.time_limit) { return true; } f_t cut_start_time = tic(); - generate_group_cover_cuts(settings, xstar, start_time); + generate_implied_indicator_cuts(settings, xstar, start_time); f_t cut_generation_time = toc(cut_start_time); if (cut_generation_time > 1.0) { - settings.log.debug("Group cover cut generation time %.2f seconds\n", cut_generation_time); + settings.log.debug("Implied indicator cut generation time %.2f seconds\n", + cut_generation_time); } } - // Generate activated capacity cuts - if (settings.activated_capacity_cuts != 0) { + // Generate capacity lifting cuts + if (settings.capacity_lifting_cuts != 0) { if (toc(start_time) >= settings.time_limit) { return true; } f_t cut_start_time = tic(); - generate_activated_capacity_cuts(settings, xstar, start_time); + generate_capacity_lifting_cuts(settings, xstar, start_time); f_t cut_generation_time = toc(cut_start_time); if (cut_generation_time > 1.0) { - settings.log.debug("Activated capacity cut generation time %.2f seconds\n", + settings.log.debug("Capacity lifting cut generation time %.2f seconds\n", cut_generation_time); } } diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index b3e382761f..463c85e63c 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -45,8 +45,8 @@ enum cut_type_t : int8_t { IMPLIED_BOUND = 5, ZERO_HALF = 6, FLOW_COVER = 7, - GROUP_COVER = 8, - ACTIVATED_CAPACITY = 9, + IMPLIED_INDICATOR = 8, + CAPACITY_LIFTING = 9, MAX_CUT_TYPE = 10 }; @@ -189,8 +189,8 @@ struct cut_info_t { "Implied Bounds", "Zero-Half ", "Flow Cover ", - "Group Cover ", - "Activated Cap "}; + "Implied Ind. ", + "Capacity Lift "}; std::array num_cuts = {0}; }; @@ -865,27 +865,27 @@ class cut_generation_t { const std::vector& xstar, f_t start_time); - // Generate group cover cuts from the variable-upper-bound gates behind each enabler row - void generate_group_cover_cuts(const simplex::simplex_solver_settings_t& settings, - const std::vector& xstar, - f_t start_time); + // Generate implied indicator cuts by aggregating an implication row over its members' indicators + void generate_implied_indicator_cuts(const simplex::simplex_solver_settings_t& settings, + const std::vector& xstar, + f_t start_time); - // Scan the user problem for gate and enabler rows. Called once, on the first cut pass. - void build_group_cover_candidates(const simplex::simplex_solver_settings_t& settings); + // Scan the user problem for implication rows. Called once, on the first cut pass. + void build_implied_indicator_candidates( + const simplex::simplex_solver_settings_t& settings); - // Generate activated capacity cuts, tying a group capacity row to the group's activation - void generate_activated_capacity_cuts( - const simplex::simplex_solver_settings_t& settings, - const std::vector& xstar, - f_t start_time); + // Generate capacity lifting cuts by lifting the complemented indicator into a capacity row + void generate_capacity_lifting_cuts(const simplex::simplex_solver_settings_t& settings, + const std::vector& xstar, + f_t start_time); - // Scan the user problem for gate and capacity rows. Called once, on the first cut pass. - void build_activated_capacity_candidates( + // Scan the user problem for capacity rows. Called once, on the first cut pass. + void build_capacity_lifting_candidates( const simplex::simplex_solver_settings_t& settings); - // The variable-upper-bound gates x_j <= z, in CSR form over the columns. Shared by both - // structural separators and built on whichever of them runs first. - void build_gate_table(const csr_matrix_t& Arow); + // Build the variable-upper-bound relation. Shared by both structural separators and built on + // whichever of them runs first. + void build_vub_table(const csr_matrix_t& Arow); void prepare_fractional_sub_conflict_graph( const simplex::simplex_solver_settings_t& settings, @@ -902,26 +902,28 @@ class cut_generation_t { std::shared_ptr>& clique_table_; omp_atomic_t* signal_extend_{nullptr}; fractional_conflict_subgraph_t sub_cg_; - // One candidate per enabler row whose selections are all gated: the head, and the distinct group - // activations behind its tail. Built once from user_problem_, then separated against xstar each - // pass. group_cover_groups_ is the flat store the offsets index into. - std::vector group_cover_heads_; - std::vector group_cover_offsets_; - std::vector group_cover_groups_; - bool group_cover_built_{false}; - // The activations gating column j are gate_activations_[gate_offsets_[j] .. gate_offsets_[j+1]), - // sorted and deduplicated, so both separators can intersect the spans directly. - std::vector gate_offsets_; - std::vector gate_activations_; - bool gates_built_{false}; - // One candidate per capacity row whose members share an activation: that activation, the cap K, - // and the row itself in <= orientation, so separation needs no second pass over the matrix. - std::vector activated_capacity_activations_; - std::vector activated_capacity_caps_; - std::vector activated_capacity_offsets_; - std::vector activated_capacity_cols_; - std::vector activated_capacity_coeffs_; - bool activated_capacity_built_{false}; + // One candidate per implication row whose members are all bounded by an indicator: the head, and + // the distinct indicators behind its tail. Built once from user_problem_, then separated against + // xstar each pass. implied_indicator_indicators_ is the flat store the offsets index into. + std::vector implied_indicator_heads_; + std::vector implied_indicator_offsets_; + std::vector implied_indicator_indicators_; + bool implied_indicator_built_{false}; + // The variable upper bounds x_j <= z, in CSR form over the columns: the indicators bounding + // column j are vub_indicators_[vub_offsets_[j] .. vub_offsets_[j+1]), sorted and deduplicated so + // both separators can intersect the spans directly. Only the binary form x_j <= z is recognised, + // which is narrower than either cut needs -- see build_vub_table. + std::vector vub_offsets_; + std::vector vub_indicators_; + bool vub_built_{false}; + // One candidate per capacity row whose members share an indicator: that indicator, the cap K, and + // the row itself in <= orientation, so separation needs no second pass over the matrix. + std::vector capacity_lifting_indicators_; + std::vector capacity_lifting_caps_; + std::vector capacity_lifting_offsets_; + std::vector capacity_lifting_cols_; + std::vector capacity_lifting_coeffs_; + bool capacity_lifting_built_{false}; }; template diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index 227a9f1a36..b7b92f309a 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -98,8 +98,8 @@ struct simplex_solver_settings_t { knapsack_cuts(-1), flow_cover_cuts(-1), implied_bound_cuts(-1), - group_cover_cuts(-1), - activated_capacity_cuts(-1), + implied_indicator_cuts(-1), + capacity_lifting_cuts(-1), clique_cuts(-1), zero_half_cuts(-1), strong_chvatal_gomory_cuts(-1), @@ -212,8 +212,8 @@ struct simplex_solver_settings_t { i_t knapsack_cuts; // -1 automatic, 0 to disable, >0 to enable knapsack cuts i_t flow_cover_cuts; // -1 automatic, 0 to disable, >0 to enable flow cover cuts i_t implied_bound_cuts; // -1 automatic, 0 to disable, >0 to enable implied bound cuts - i_t group_cover_cuts; // 0 to disable, >0 to enable group cover cuts - i_t activated_capacity_cuts; // 0 to disable, >0 to enable activated capacity cuts + i_t implied_indicator_cuts; // 0 to disable, >0 to enable implied indicator cuts + i_t capacity_lifting_cuts; // 0 to disable, >0 to enable capacity lifting cuts i_t clique_cuts; // -1 automatic, 0 to disable, >0 to enable clique cuts i_t zero_half_cuts; // -1 automatic, 0 to disable, >0 to enable zero-half cuts i_t strong_chvatal_gomory_cuts; // -1 automatic, 0 to disable, >0 to enable strong Chvatal Gomory diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index 0a556ef0da..6cb9d503fb 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -189,8 +189,8 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_MIP_CLIQUE_CUTS, &mip_settings.clique_cuts, -1, 1, -1}, {CUOPT_MIP_ZERO_HALF_CUTS, &mip_settings.zero_half_cuts, -1, 1, -1}, {CUOPT_MIP_IMPLIED_BOUND_CUTS, &mip_settings.implied_bound_cuts, -1, 1, -1}, - {CUOPT_MIP_GROUP_COVER_CUTS, &mip_settings.group_cover_cuts, -1, 1, -1}, - {CUOPT_MIP_ACTIVATED_CAPACITY_CUTS, &mip_settings.activated_capacity_cuts, -1, 1, -1}, + {CUOPT_MIP_IMPLIED_INDICATOR_CUTS, &mip_settings.implied_indicator_cuts, -1, 1, -1}, + {CUOPT_MIP_CAPACITY_LIFTING_CUTS, &mip_settings.capacity_lifting_cuts, -1, 1, -1}, {CUOPT_MIP_STRONG_CHVATAL_GOMORY_CUTS, &mip_settings.strong_chvatal_gomory_cuts, -1, 1, -1}, {CUOPT_MIP_REDUCED_COST_STRENGTHENING, &mip_settings.reduced_cost_strengthening, -1, std::numeric_limits::max(), -1}, {CUOPT_MIP_RINS, &mip_settings.submip_params.rins, -1, 1, -1}, diff --git a/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh b/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh index c6f8453c9e..76ab9dc07d 100644 --- a/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh +++ b/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh @@ -112,8 +112,8 @@ class sub_mip_recombiner_t : public recombiner_t { branch_and_bound_settings.max_cut_passes = 0; branch_and_bound_settings.clique_cuts = 0; branch_and_bound_settings.zero_half_cuts = 0; - branch_and_bound_settings.group_cover_cuts = 0; - branch_and_bound_settings.activated_capacity_cuts = 0; + branch_and_bound_settings.implied_indicator_cuts = 0; + branch_and_bound_settings.capacity_lifting_cuts = 0; branch_and_bound_settings.inside_submip = 1; branch_and_bound_settings.submip_settings.rins = 0; branch_and_bound_settings.submip_settings.rens = 0; diff --git a/cpp/src/mip_heuristics/solver.cu b/cpp/src/mip_heuristics/solver.cu index c5675db075..2477412d8c 100644 --- a/cpp/src/mip_heuristics/solver.cu +++ b/cpp/src/mip_heuristics/solver.cu @@ -366,13 +366,13 @@ solution_t mip_solver_t::run_solver() } branch_and_bound_settings.mixed_integer_gomory_cuts = context.settings.mixed_integer_gomory_cuts; - branch_and_bound_settings.knapsack_cuts = context.settings.knapsack_cuts; - branch_and_bound_settings.flow_cover_cuts = context.settings.flow_cover_cuts; - branch_and_bound_settings.implied_bound_cuts = context.settings.implied_bound_cuts; - branch_and_bound_settings.group_cover_cuts = context.settings.group_cover_cuts; - branch_and_bound_settings.activated_capacity_cuts = context.settings.activated_capacity_cuts; - branch_and_bound_settings.clique_cuts = context.settings.clique_cuts; - branch_and_bound_settings.zero_half_cuts = context.settings.zero_half_cuts; + branch_and_bound_settings.knapsack_cuts = context.settings.knapsack_cuts; + branch_and_bound_settings.flow_cover_cuts = context.settings.flow_cover_cuts; + branch_and_bound_settings.implied_bound_cuts = context.settings.implied_bound_cuts; + branch_and_bound_settings.implied_indicator_cuts = context.settings.implied_indicator_cuts; + branch_and_bound_settings.capacity_lifting_cuts = context.settings.capacity_lifting_cuts; + branch_and_bound_settings.clique_cuts = context.settings.clique_cuts; + branch_and_bound_settings.zero_half_cuts = context.settings.zero_half_cuts; branch_and_bound_settings.strong_chvatal_gomory_cuts = context.settings.strong_chvatal_gomory_cuts; branch_and_bound_settings.cut_change_threshold = context.settings.cut_change_threshold; From 01584e1f74bcf5c05dd64bf059bddf793ae3874c Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Mon, 28 Sep 2026 12:05:20 +0200 Subject: [PATCH 06/12] transform the cuts into presolve reductions Signed-off-by: Nicolas L. Guidotti --- .../mathematical_optimization/constants.h | 2 - .../mip/solver_settings.hpp | 2 - cpp/src/branch_and_bound/branch_and_bound.cpp | 2 - cpp/src/cuts/cuts.cpp | 401 ------------------ cpp/src/cuts/cuts.hpp | 52 +-- .../dual_simplex/simplex_solver_settings.hpp | 4 - cpp/src/math_optimization/solver_settings.cu | 2 - cpp/src/mip_heuristics/CMakeLists.txt | 1 + .../diversity/recombiners/sub_mip.cuh | 2 - .../presolve/indicator_strengthening.cpp | 266 ++++++++++++ .../presolve/indicator_strengthening.hpp | 30 ++ .../presolve/third_party_presolve.cpp | 7 + cpp/src/mip_heuristics/solver.cu | 12 +- 13 files changed, 311 insertions(+), 472 deletions(-) create mode 100644 cpp/src/mip_heuristics/presolve/indicator_strengthening.cpp create mode 100644 cpp/src/mip_heuristics/presolve/indicator_strengthening.hpp diff --git a/cpp/include/cuopt/mathematical_optimization/constants.h b/cpp/include/cuopt/mathematical_optimization/constants.h index 7ddc5eaaf8..3656791a98 100644 --- a/cpp/include/cuopt/mathematical_optimization/constants.h +++ b/cpp/include/cuopt/mathematical_optimization/constants.h @@ -80,8 +80,6 @@ #define CUOPT_MIP_CLIQUE_CUTS "mip_clique_cuts" #define CUOPT_MIP_ZERO_HALF_CUTS "mip_zero_half_cuts" #define CUOPT_MIP_STRONG_CHVATAL_GOMORY_CUTS "mip_strong_chvatal_gomory_cuts" -#define CUOPT_MIP_IMPLIED_INDICATOR_CUTS "mip_implied_indicator_cuts" -#define CUOPT_MIP_CAPACITY_LIFTING_CUTS "mip_capacity_lifting_cuts" #define CUOPT_MIP_REDUCED_COST_STRENGTHENING "mip_reduced_cost_strengthening" #define CUOPT_MIP_RINS "mip_rins" #define CUOPT_MIP_RENS "mip_rens" diff --git a/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp b/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp index 0d28258164..7ed45f1f9b 100644 --- a/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp +++ b/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp @@ -135,8 +135,6 @@ class mip_solver_settings_t { i_t clique_cuts = -1; i_t zero_half_cuts = -1; i_t implied_bound_cuts = -1; - i_t implied_indicator_cuts = -1; - i_t capacity_lifting_cuts = -1; i_t strong_chvatal_gomory_cuts = -1; i_t reduced_cost_strengthening = -1; i_t objective_step = 1; // 0 = disable objective step tightening, 1 = enable diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 6a170154ae..a320cc0602 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -2321,8 +2321,6 @@ void branch_and_bound_t::solve_submip(diving_worker_t* worke submip_settings.reliability_branching = 0; submip_settings.clique_cuts = 0; submip_settings.zero_half_cuts = 0; - submip_settings.implied_indicator_cuts = 0; - submip_settings.capacity_lifting_cuts = 0; submip_settings.inside_submip = 1; submip_settings.strong_branching_simplex_iteration_limit = 50; submip_settings.inside_root_node = 0; diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index 986ad4510d..e9f51666dc 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -3232,383 +3232,6 @@ void cut_generation_t::generate_implied_bound_cuts( } } -// The variable upper bounds x_j <= u_j z linking a member to the indicator that must be on before -// it can be chosen. Only the binary unit form x_j <= z is recognised here -- both columns binary, -// coefficients +-1, right-hand side zero -- which is narrower than either separator needs: the -// implied indicator cut holds for any u_j > 0, and the capacity lifting cut does not need the -// member to be integral at all. Widening the scan would find more structure than it does today. -template -void cut_generation_t::build_vub_table(const csr_matrix_t& Arow) -{ - if (vub_built_) { return; } - vub_built_ = true; - - const i_t num_rows = user_problem_.num_rows; - const i_t num_cols = user_problem_.num_cols; - - const bool has_ranges = user_problem_.num_range_rows > 0; - std::vector is_range(has_ranges ? num_rows : 0, 0); - for (i_t k = 0; k < user_problem_.num_range_rows; ++k) { - is_range[user_problem_.range_rows[k]] = 1; - } - - std::vector> vubs; - for (i_t row = 0; row < num_rows; ++row) { - if (has_ranges && is_range[row]) { continue; } - if (Arow.row_length(row) != 2) { continue; } - const char sense = user_problem_.row_sense[row]; - if (sense != 'L' && sense != 'G') { continue; } - // Row sense in <= orientation: +1 when the row reads a.x <= rhs, -1 when a.x >= rhs. - const i_t direction = sense == 'L' ? 1 : -1; - if (direction * user_problem_.rhs[row] != 0.0) { continue; } - - i_t member = -1, indicator = -1; - for (i_t p = Arow.row_start[row]; p < Arow.row_start[row + 1]; ++p) { - const i_t col = Arow.j[p]; - if (user_problem_.var_types[col] == variable_type_t::CONTINUOUS || - user_problem_.lower[col] != 0.0 || user_problem_.upper[col] != 1.0) { - member = indicator = -1; - break; - } - const f_t v = direction * Arow.x[p]; - if (v == 1.0) { - member = col; - } else if (v == -1.0) { - indicator = col; - } - } - if (member >= 0 && indicator >= 0) { vubs.emplace_back(member, indicator); } - } - if (vubs.empty()) { return; } - - vub_offsets_.assign(num_cols + 1, 0); - for (const auto& [member, indicator] : vubs) { - ++vub_offsets_[member + 1]; - } - for (i_t col = 0; col < num_cols; ++col) { - vub_offsets_[col + 1] += vub_offsets_[col]; - } - vub_indicators_.resize(vubs.size()); - std::vector cursor(vub_offsets_.begin(), vub_offsets_.end() - 1); - for (const auto& [member, indicator] : vubs) { - vub_indicators_[cursor[member]++] = indicator; - } - - // Sort each column's span and drop rows stating the same bound twice, compacting in place. The - // write cursor trails the read cursor because it only advances on a kept entry. - i_t out = 0; - for (i_t col = 0; col < num_cols; ++col) { - const i_t begin = vub_offsets_[col]; - const i_t end = vub_offsets_[col + 1]; - vub_offsets_[col] = out; - std::sort(vub_indicators_.begin() + begin, vub_indicators_.begin() + end); - for (i_t p = begin; p < end; ++p) { - if (p == begin || vub_indicators_[p] != vub_indicators_[p - 1]) { - vub_indicators_[out++] = vub_indicators_[p]; - } - } - } - vub_offsets_[num_cols] = out; - vub_indicators_.resize(out); -} - -// An implied indicator cut aggregates an implication row over the indicators of its members. Where -// the model carries -// -// y <= sum_{j in S} x_j (implication) and x_j <= z_{g(j)} for every j in S, -// -// a binary y that is one forces some single x_j to one, which forces that member's indicator to -// one, so -// -// y <= sum_{g in D} z_g, D = the distinct indicators over S. -// -// Counting each indicator once is where the strength is: chaining the members through their own -// variable upper bounds gives y <= sum_j z_{g(j)}, which counts one indicator once per member. -template -void cut_generation_t::build_implied_indicator_candidates( - const simplex_solver_settings_t& settings) -{ - implied_indicator_built_ = true; - - const i_t num_rows = user_problem_.num_rows; - const i_t num_cols = user_problem_.num_cols; - if (num_rows <= 0 || num_cols <= 0) { return; } - if (user_problem_.var_types.size() != (size_t)num_cols || - user_problem_.row_sense.size() != (size_t)num_rows || - user_problem_.rhs.size() != (size_t)num_rows) { - return; - } - - csr_matrix_t Arow(num_rows, num_cols, user_problem_.A.col_start[num_cols]); - user_problem_.A.to_compressed_row(Arow); - - const bool has_ranges = user_problem_.num_range_rows > 0; - std::vector is_range(has_ranges ? num_rows : 0, 0); - for (i_t k = 0; k < user_problem_.num_range_rows; ++k) { - is_range[user_problem_.range_rows[k]] = 1; - } - - build_vub_table(Arow); - if (vub_indicators_.empty()) { return; } - - // The implication rows, one +1 head against a tail of -1 members. - std::vector indicators; - implied_indicator_offsets_.push_back(0); - for (i_t row = 0; row < num_rows; ++row) { - if (has_ranges && is_range[row]) { continue; } - const char sense = user_problem_.row_sense[row]; - if (sense != 'L' && sense != 'G') { continue; } - // Row sense in <= orientation: +1 when the row reads a.x <= rhs, -1 when a.x >= rhs. - const i_t direction = sense == 'L' ? 1 : -1; - const i_t len = Arow.row_length(row); - if (len < 3) { continue; } - if (direction * user_problem_.rhs[row] != 0.0) { continue; } - - i_t head = -1; - bool usable = true; - indicators.clear(); - for (i_t k = Arow.row_start[row]; k < Arow.row_start[row + 1] && usable; ++k) { - const i_t col = Arow.j[k]; - const f_t v = direction * Arow.x[k]; - if (user_problem_.var_types[col] == variable_type_t::CONTINUOUS || - user_problem_.lower[col] != 0.0 || user_problem_.upper[col] != 1.0) { - usable = false; - } else if (v == 1.0) { - usable = head < 0; - head = col; - } else if (v == -1.0) { - // A member with no variable upper bound leaves nothing to aggregate through; keep the - // candidate by standing in the member itself, which the implication row already bounds - // the head by. - if (vub_offsets_[col] == vub_offsets_[col + 1]) { - indicators.push_back(col); - } else { - for (i_t p = vub_offsets_[col]; p < vub_offsets_[col + 1]; ++p) { - indicators.push_back(vub_indicators_[p]); - } - } - } else { - usable = false; - } - } - if (!usable || head < 0 || indicators.empty()) { continue; } - - std::sort(indicators.begin(), indicators.end()); - indicators.erase(std::unique(indicators.begin(), indicators.end()), indicators.end()); - // Nothing was merged, so the cut is the sum of the variable upper bounds the LP already has. - if (indicators.size() >= (size_t)(len - 1)) { continue; } - if (std::binary_search(indicators.begin(), indicators.end(), head)) { continue; } - - implied_indicator_heads_.push_back(head); - implied_indicator_indicators_.insert( - implied_indicator_indicators_.end(), indicators.begin(), indicators.end()); - implied_indicator_offsets_.push_back(implied_indicator_indicators_.size()); - } - if (implied_indicator_heads_.empty()) { return; } - - settings.log.print_format("Implied indicator: {} candidate cuts over {} variable upper bounds\n", - implied_indicator_heads_.size(), - vub_indicators_.size()); -} - -template -void cut_generation_t::generate_implied_indicator_cuts( - const simplex_solver_settings_t& settings, - const std::vector& xstar, - f_t start_time) -{ - if (!implied_indicator_built_) { build_implied_indicator_candidates(settings); } - if (implied_indicator_heads_.empty()) { return; } - - const f_t tol = 1e-4; - i_t num_cuts = 0; - const i_t n = implied_indicator_heads_.size(); - for (i_t k = 0; k < n; ++k) { - if ((k & 0xFF) == 0 && toc(start_time) >= settings.time_limit) { return; } - const i_t head = implied_indicator_heads_[k]; - f_t activity = xstar[head]; - for (i_t p = implied_indicator_offsets_[k]; p < implied_indicator_offsets_[k + 1]; ++p) { - activity -= xstar[implied_indicator_indicators_[p]]; - } - if (activity <= tol) { continue; } - - // add_cut expects cut'x >= rhs, so the cut head - sum_g z_g <= 0 is emitted negated. - inequality_t cut; - cut.push_back(head, -1.0); - for (i_t p = implied_indicator_offsets_[k]; p < implied_indicator_offsets_[k + 1]; ++p) { - cut.push_back(implied_indicator_indicators_[p], 1.0); - } - cut.rhs = 0.0; - cut_pool_.add_cut(cut_type_t::IMPLIED_INDICATOR, cut); - num_cuts++; - } - - if (num_cuts > 0) { settings.log.debug("Generated %d implied indicator cuts\n", num_cuts); } -} - -// A capacity lifting cut lifts the complemented indicator into a capacity row. Where the model -// carries -// -// sum_{i in S} x_i - s <= K (capacity) and x_i <= z for every i in S, -// -// sequential lifting of zbar = 1 - z into the capacity row asks for the largest coefficient a -// keeping sum_i x_i - s + a zbar <= K valid. At zbar = 1 the variable upper bounds force every -// x_i to zero and -s is non-positive, so the left side maxes out at 0 and a may rise to K: -// -// sum_{i in S} x_i - s <= K z. -// -// K is therefore the maximal lifting coefficient, not a choice. Valid for integral z only. Unlike -// the implication rows behind an implied indicator cut, the capacity row stays in the model; the -// cut is the lifted copy of it. The capacity row also bounds the violation, since it already gives -// sum_i x*_i - s* <= K and hence a violation of at most K(1 - z*): the cut can only bite where the -// row is near-tight and the indicator is well below one. -template -void cut_generation_t::build_capacity_lifting_candidates( - const simplex_solver_settings_t& settings) -{ - capacity_lifting_built_ = true; - - const i_t num_rows = user_problem_.num_rows; - const i_t num_cols = user_problem_.num_cols; - if (num_rows <= 0 || num_cols <= 0) { return; } - if (user_problem_.var_types.size() != (size_t)num_cols || - user_problem_.row_sense.size() != (size_t)num_rows || - user_problem_.rhs.size() != (size_t)num_rows) { - return; - } - - csr_matrix_t Arow(num_rows, num_cols, user_problem_.A.col_start[num_cols]); - user_problem_.A.to_compressed_row(Arow); - - const bool has_ranges = user_problem_.num_range_rows > 0; - std::vector is_range(has_ranges ? num_rows : 0, 0); - for (i_t k = 0; k < user_problem_.num_range_rows; ++k) { - is_range[user_problem_.range_rows[k]] = 1; - } - - build_vub_table(Arow); - if (vub_indicators_.empty()) { return; } - - std::vector members; - std::vector support; - std::vector common; - std::vector intersection; - - capacity_lifting_offsets_.push_back(0); - for (i_t row = 0; row < num_rows; ++row) { - if (has_ranges && is_range[row]) { continue; } - const char sense = user_problem_.row_sense[row]; - if (sense != 'L' && sense != 'G') { continue; } - // Row sense in <= orientation: +1 when the row reads a.x <= rhs, -1 when a.x >= rhs. - const i_t direction = sense == 'L' ? 1 : -1; - if (Arow.row_length(row) < 2) { continue; } - - const f_t capacity = direction * user_problem_.rhs[row]; - if (capacity <= 0.0) { continue; } - - members.clear(); - support.clear(); - bool usable = true; - for (i_t p = Arow.row_start[row]; p < Arow.row_start[row + 1] && usable; ++p) { - const i_t col = Arow.j[p]; - const f_t v = direction * Arow.x[p]; - support.push_back(col); - if (user_problem_.var_types[col] != variable_type_t::CONTINUOUS) { - usable = user_problem_.lower[col] == 0.0 && user_problem_.upper[col] == 1.0 && v == 1.0; - if (usable) { members.push_back(col); } - } else { - // A relaxing term: non-positive over the whole box, so it cannot violate the row at z = 0. - usable = v < 0.0 && user_problem_.lower[col] >= 0.0; - } - } - if (!usable || members.size() < 2) { continue; } - // At or above its own support size the cap is implied by the variable upper bounds, and so is - // K z, leaving nothing to lift. - const f_t n_members = members.size(); - if (capacity >= n_members) { continue; } - - common.assign(vub_indicators_.begin() + vub_offsets_[members[0]], - vub_indicators_.begin() + vub_offsets_[members[0] + 1]); - for (size_t k = 1; k < members.size() && !common.empty(); ++k) { - const i_t col = members[k]; - intersection.clear(); - std::set_intersection(common.begin(), - common.end(), - vub_indicators_.begin() + vub_offsets_[col], - vub_indicators_.begin() + vub_offsets_[col + 1], - std::back_inserter(intersection)); - common.swap(intersection); - } - if (common.empty()) { continue; } - - std::sort(support.begin(), support.end()); - i_t indicator = -1; - for (i_t z : common) { - if (std::binary_search(support.begin(), support.end(), z)) { continue; } - indicator = z; - break; - } - if (indicator < 0) { continue; } - - capacity_lifting_indicators_.push_back(indicator); - capacity_lifting_caps_.push_back(capacity); - for (i_t p = Arow.row_start[row]; p < Arow.row_start[row + 1]; ++p) { - capacity_lifting_cols_.push_back(Arow.j[p]); - capacity_lifting_coeffs_.push_back(direction * Arow.x[p]); - } - capacity_lifting_offsets_.push_back(capacity_lifting_cols_.size()); - } - if (capacity_lifting_indicators_.empty()) { return; } - - settings.log.print_format("Capacity lifting: {} candidate cuts over {} variable upper bounds\n", - capacity_lifting_indicators_.size(), - vub_indicators_.size()); -} - -template -void cut_generation_t::generate_capacity_lifting_cuts( - const simplex_solver_settings_t& settings, - const std::vector& xstar, - f_t start_time) -{ - if (!capacity_lifting_built_) { build_capacity_lifting_candidates(settings); } - if (capacity_lifting_indicators_.empty()) { return; } - - // The lifted row carries a coefficient of K against +-1 everywhere else, so it costs the node LP - // far more per row than a unit-coefficient cut and is held to a higher bar than cut_pool_t's - // global min_cut_distance_ of 1e-4. Measured over the first root passes: the cuts that carry - // tsmc-setcover-3, where the cap is 3 of 25 members and the rows bind constantly, have distance - // 0.23 upwards; on tsmc-setcover-2 the cap is 12 of 25, the rows rarely bind, and three quarters - // sit below 0.14 while buying no bound the implied indicator cuts do not already have. - const f_t min_distance = 0.2; - i_t num_cuts = 0; - const i_t n = capacity_lifting_indicators_.size(); - for (i_t k = 0; k < n; ++k) { - if ((k & 0xFF) == 0 && toc(start_time) >= settings.time_limit) { return; } - const i_t indicator = capacity_lifting_indicators_[k]; - const f_t capacity = capacity_lifting_caps_[k]; - f_t activity = -capacity * xstar[indicator]; - f_t norm = capacity * capacity; - for (i_t p = capacity_lifting_offsets_[k]; p < capacity_lifting_offsets_[k + 1]; ++p) { - activity += capacity_lifting_coeffs_[p] * xstar[capacity_lifting_cols_[p]]; - norm += capacity_lifting_coeffs_[p] * capacity_lifting_coeffs_[p]; - } - if (activity <= min_distance * std::sqrt(norm)) { continue; } - - // add_cut expects cut'x >= rhs, so the row sum_p a_p x_p - K z <= 0 is emitted negated. - inequality_t cut; - for (i_t p = capacity_lifting_offsets_[k]; p < capacity_lifting_offsets_[k + 1]; ++p) { - cut.push_back(capacity_lifting_cols_[p], -capacity_lifting_coeffs_[p]); - } - cut.push_back(indicator, capacity); - cut.rhs = 0.0; - cut_pool_.add_cut(cut_type_t::CAPACITY_LIFTING, cut); - num_cuts++; - } - - if (num_cuts > 0) { settings.log.debug("Generated %d capacity lifting cuts\n", num_cuts); } -} - namespace { // Total probing-edge budget from the byte cap and the remaining work headroom @@ -3983,30 +3606,6 @@ bool cut_generation_t::generate_cuts(const lp_problem_t& lp, } } - // Generate implied indicator cuts - if (settings.implied_indicator_cuts != 0) { - if (toc(start_time) >= settings.time_limit) { return true; } - f_t cut_start_time = tic(); - generate_implied_indicator_cuts(settings, xstar, start_time); - f_t cut_generation_time = toc(cut_start_time); - if (cut_generation_time > 1.0) { - settings.log.debug("Implied indicator cut generation time %.2f seconds\n", - cut_generation_time); - } - } - - // Generate capacity lifting cuts - if (settings.capacity_lifting_cuts != 0) { - if (toc(start_time) >= settings.time_limit) { return true; } - f_t cut_start_time = tic(); - generate_capacity_lifting_cuts(settings, xstar, start_time); - f_t cut_generation_time = toc(cut_start_time); - if (cut_generation_time > 1.0) { - settings.log.debug("Capacity lifting cut generation time %.2f seconds\n", - cut_generation_time); - } - } - // Build the fractional conflict-graph subgraph once (resolving the async // clique-table future on the way) so both clique-cut and zero-half cut // separators consume the same vertex/weight/adjacency tables instead of diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index 463c85e63c..ca87e26c39 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -45,9 +45,7 @@ enum cut_type_t : int8_t { IMPLIED_BOUND = 5, ZERO_HALF = 6, FLOW_COVER = 7, - IMPLIED_INDICATOR = 8, - CAPACITY_LIFTING = 9, - MAX_CUT_TYPE = 10 + MAX_CUT_TYPE = 8 }; template @@ -188,9 +186,7 @@ struct cut_info_t { "Clique ", "Implied Bounds", "Zero-Half ", - "Flow Cover ", - "Implied Ind. ", - "Capacity Lift "}; + "Flow Cover "}; std::array num_cuts = {0}; }; @@ -865,28 +861,6 @@ class cut_generation_t { const std::vector& xstar, f_t start_time); - // Generate implied indicator cuts by aggregating an implication row over its members' indicators - void generate_implied_indicator_cuts(const simplex::simplex_solver_settings_t& settings, - const std::vector& xstar, - f_t start_time); - - // Scan the user problem for implication rows. Called once, on the first cut pass. - void build_implied_indicator_candidates( - const simplex::simplex_solver_settings_t& settings); - - // Generate capacity lifting cuts by lifting the complemented indicator into a capacity row - void generate_capacity_lifting_cuts(const simplex::simplex_solver_settings_t& settings, - const std::vector& xstar, - f_t start_time); - - // Scan the user problem for capacity rows. Called once, on the first cut pass. - void build_capacity_lifting_candidates( - const simplex::simplex_solver_settings_t& settings); - - // Build the variable-upper-bound relation. Shared by both structural separators and built on - // whichever of them runs first. - void build_vub_table(const csr_matrix_t& Arow); - void prepare_fractional_sub_conflict_graph( const simplex::simplex_solver_settings_t& settings, const std::vector& xstar, @@ -902,28 +876,6 @@ class cut_generation_t { std::shared_ptr>& clique_table_; omp_atomic_t* signal_extend_{nullptr}; fractional_conflict_subgraph_t sub_cg_; - // One candidate per implication row whose members are all bounded by an indicator: the head, and - // the distinct indicators behind its tail. Built once from user_problem_, then separated against - // xstar each pass. implied_indicator_indicators_ is the flat store the offsets index into. - std::vector implied_indicator_heads_; - std::vector implied_indicator_offsets_; - std::vector implied_indicator_indicators_; - bool implied_indicator_built_{false}; - // The variable upper bounds x_j <= z, in CSR form over the columns: the indicators bounding - // column j are vub_indicators_[vub_offsets_[j] .. vub_offsets_[j+1]), sorted and deduplicated so - // both separators can intersect the spans directly. Only the binary form x_j <= z is recognised, - // which is narrower than either cut needs -- see build_vub_table. - std::vector vub_offsets_; - std::vector vub_indicators_; - bool vub_built_{false}; - // One candidate per capacity row whose members share an indicator: that indicator, the cap K, and - // the row itself in <= orientation, so separation needs no second pass over the matrix. - std::vector capacity_lifting_indicators_; - std::vector capacity_lifting_caps_; - std::vector capacity_lifting_offsets_; - std::vector capacity_lifting_cols_; - std::vector capacity_lifting_coeffs_; - bool capacity_lifting_built_{false}; }; template diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index b7b92f309a..8b3eba56d3 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -98,8 +98,6 @@ struct simplex_solver_settings_t { knapsack_cuts(-1), flow_cover_cuts(-1), implied_bound_cuts(-1), - implied_indicator_cuts(-1), - capacity_lifting_cuts(-1), clique_cuts(-1), zero_half_cuts(-1), strong_chvatal_gomory_cuts(-1), @@ -212,8 +210,6 @@ struct simplex_solver_settings_t { i_t knapsack_cuts; // -1 automatic, 0 to disable, >0 to enable knapsack cuts i_t flow_cover_cuts; // -1 automatic, 0 to disable, >0 to enable flow cover cuts i_t implied_bound_cuts; // -1 automatic, 0 to disable, >0 to enable implied bound cuts - i_t implied_indicator_cuts; // 0 to disable, >0 to enable implied indicator cuts - i_t capacity_lifting_cuts; // 0 to disable, >0 to enable capacity lifting cuts i_t clique_cuts; // -1 automatic, 0 to disable, >0 to enable clique cuts i_t zero_half_cuts; // -1 automatic, 0 to disable, >0 to enable zero-half cuts i_t strong_chvatal_gomory_cuts; // -1 automatic, 0 to disable, >0 to enable strong Chvatal Gomory diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index 6cb9d503fb..5a4ab72c32 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -189,8 +189,6 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_MIP_CLIQUE_CUTS, &mip_settings.clique_cuts, -1, 1, -1}, {CUOPT_MIP_ZERO_HALF_CUTS, &mip_settings.zero_half_cuts, -1, 1, -1}, {CUOPT_MIP_IMPLIED_BOUND_CUTS, &mip_settings.implied_bound_cuts, -1, 1, -1}, - {CUOPT_MIP_IMPLIED_INDICATOR_CUTS, &mip_settings.implied_indicator_cuts, -1, 1, -1}, - {CUOPT_MIP_CAPACITY_LIFTING_CUTS, &mip_settings.capacity_lifting_cuts, -1, 1, -1}, {CUOPT_MIP_STRONG_CHVATAL_GOMORY_CUTS, &mip_settings.strong_chvatal_gomory_cuts, -1, 1, -1}, {CUOPT_MIP_REDUCED_COST_STRENGTHENING, &mip_settings.reduced_cost_strengthening, -1, std::numeric_limits::max(), -1}, {CUOPT_MIP_RINS, &mip_settings.submip_params.rins, -1, 1, -1}, diff --git a/cpp/src/mip_heuristics/CMakeLists.txt b/cpp/src/mip_heuristics/CMakeLists.txt index 187017fb14..2c999dcab1 100644 --- a/cpp/src/mip_heuristics/CMakeLists.txt +++ b/cpp/src/mip_heuristics/CMakeLists.txt @@ -15,6 +15,7 @@ set(MIP_LP_NECESSARY_FILES ${CMAKE_CURRENT_SOURCE_DIR}/presolve/single_lock_dual_aggregation.cpp ${CMAKE_CURRENT_SOURCE_DIR}/presolve/gf2_presolve.cpp ${CMAKE_CURRENT_SOURCE_DIR}/presolve/bhw_coeff_reduce.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/presolve/indicator_strengthening.cpp ${CMAKE_CURRENT_SOURCE_DIR}/solution/solution.cu ${CMAKE_CURRENT_SOURCE_DIR}/presolve/conflict_graph/clique_table.cu ) diff --git a/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh b/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh index 76ab9dc07d..7990465b81 100644 --- a/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh +++ b/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh @@ -112,8 +112,6 @@ class sub_mip_recombiner_t : public recombiner_t { branch_and_bound_settings.max_cut_passes = 0; branch_and_bound_settings.clique_cuts = 0; branch_and_bound_settings.zero_half_cuts = 0; - branch_and_bound_settings.implied_indicator_cuts = 0; - branch_and_bound_settings.capacity_lifting_cuts = 0; branch_and_bound_settings.inside_submip = 1; branch_and_bound_settings.submip_settings.rins = 0; branch_and_bound_settings.submip_settings.rens = 0; diff --git a/cpp/src/mip_heuristics/presolve/indicator_strengthening.cpp b/cpp/src/mip_heuristics/presolve/indicator_strengthening.cpp new file mode 100644 index 0000000000..5939dcdad0 --- /dev/null +++ b/cpp/src/mip_heuristics/presolve/indicator_strengthening.cpp @@ -0,0 +1,266 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#include "indicator_strengthening.hpp" + +#include +#include + +#include +#include +#include +#include +#include + +namespace cuopt::mathematical_optimization::mip { + +template +void strengthen_indicators(papilo::Problem& problem) +{ + const auto& constraint_matrix = problem.getConstraintMatrix(); + const auto& lhs_values = constraint_matrix.getLeftHandSides(); + const auto& rhs_values = constraint_matrix.getRightHandSides(); + const auto& row_flags = constraint_matrix.getRowFlags(); + const auto& domains = problem.getVariableDomains(); + const auto& col_flags = domains.flags; + const auto& lower_bounds = domains.lower_bounds; + const auto& upper_bounds = domains.upper_bounds; + + const int num_rows = constraint_matrix.getNRows(); + const int num_cols = problem.getNCols(); + if (num_rows <= 0 || num_cols <= 0) { return; } + + auto is_free_binary = [&](int col) { + const auto& flags = col_flags[col]; + return flags.test(papilo::ColFlag::kIntegral) && !flags.test(papilo::ColFlag::kLbInf) && + !flags.test(papilo::ColFlag::kUbInf) && !flags.test(papilo::ColFlag::kFixed) && + lower_bounds[col] == 0.0 && upper_bounds[col] == 1.0; + }; + + // +1 when the stored row reads a.x <= rhs, -1 when it reads a.x >= lhs, 0 for equations, ranges + // and free rows. + auto orientation = [&](int row) { + const bool lhs_infinite = row_flags[row].test(papilo::RowFlag::kLhsInf); + const bool rhs_infinite = row_flags[row].test(papilo::RowFlag::kRhsInf); + if (lhs_infinite == rhs_infinite) { return 0; } + return lhs_infinite ? 1 : -1; + }; + + // The variable upper bounds x <= z over binaries, in CSR form over the members: the indicators + // bounding column j are vub_indicators[vub_offsets[j] .. vub_offsets[j+1]), sorted and unique. + std::vector> vubs; + for (int row = 0; row < num_rows; ++row) { + const int direction = orientation(row); + if (direction == 0) { continue; } + auto row_coefficients = constraint_matrix.getRowCoefficients(row); + if (row_coefficients.getLength() != 2) { continue; } + const f_t side = direction == 1 ? rhs_values[row] : lhs_values[row]; + if (side != 0.0) { continue; } + + const int* indices = row_coefficients.getIndices(); + const f_t* values = row_coefficients.getValues(); + int member = -1, indicator = -1; + for (int p = 0; p < 2; ++p) { + if (!is_free_binary(indices[p])) { + member = indicator = -1; + break; + } + const f_t v = direction * values[p]; + if (v == 1.0) { + member = indices[p]; + } else if (v == -1.0) { + indicator = indices[p]; + } + } + if (member >= 0 && indicator >= 0) { vubs.emplace_back(member, indicator); } + } + if (vubs.empty()) { return; } + + std::sort(vubs.begin(), vubs.end()); + vubs.erase(std::unique(vubs.begin(), vubs.end()), vubs.end()); + std::vector vub_offsets(num_cols + 1, 0); + std::vector vub_indicators(vubs.size()); + for (size_t k = 0; k < vubs.size(); ++k) { + ++vub_offsets[vubs[k].first + 1]; + vub_indicators[k] = vubs[k].second; + } + for (int col = 0; col < num_cols; ++col) { + vub_offsets[col + 1] += vub_offsets[col]; + } + + papilo::Vec> entries; + entries.reserve(constraint_matrix.getNnz() + num_rows); + papilo::Vec lhs(lhs_values.begin(), lhs_values.end()); + papilo::Vec rhs(rhs_values.begin(), rhs_values.end()); + papilo::Vec flags(row_flags.begin(), row_flags.end()); + + std::vector implied_heads; + std::vector implied_offsets{0}; + std::vector implied_indicators; + int num_lifted = 0; + + std::vector indicators; + std::vector members; + std::vector support; + std::vector common; + std::vector intersection; + + for (int row = 0; row < num_rows; ++row) { + auto row_coefficients = constraint_matrix.getRowCoefficients(row); + const int len = row_coefficients.getLength(); + const int* indices = row_coefficients.getIndices(); + const f_t* values = row_coefficients.getValues(); + for (int p = 0; p < len; ++p) { + entries.emplace_back(row, indices[p], values[p]); + } + + const int direction = orientation(row); + if (direction == 0 || len < 2) { continue; } + const f_t capacity = direction * (direction == 1 ? rhs_values[row] : lhs_values[row]); + + // Implication row y - sum_{j in S} x_j <= 0: one +1 head against a tail of -1 members. + if (capacity == 0.0) { + if (len < 3) { continue; } + int head = -1; + bool usable = true; + indicators.clear(); + for (int p = 0; p < len && usable; ++p) { + const int col = indices[p]; + const f_t v = direction * values[p]; + if (!is_free_binary(col)) { + usable = false; + } else if (v == 1.0) { + usable = head < 0; + head = col; + } else if (v == -1.0) { + if (vub_offsets[col] == vub_offsets[col + 1]) { + indicators.push_back(col); + } else { + indicators.insert(indicators.end(), + vub_indicators.begin() + vub_offsets[col], + vub_indicators.begin() + vub_offsets[col + 1]); + } + } else { + usable = false; + } + } + if (!usable || head < 0 || indicators.empty()) { continue; } + + std::sort(indicators.begin(), indicators.end()); + indicators.erase(std::unique(indicators.begin(), indicators.end()), indicators.end()); + const size_t num_members = len - 1; + if (indicators.size() >= num_members) { continue; } + if (std::binary_search(indicators.begin(), indicators.end(), head)) { continue; } + + implied_heads.push_back(head); + implied_indicators.insert(implied_indicators.end(), indicators.begin(), indicators.end()); + implied_offsets.push_back(implied_indicators.size()); + continue; + } + if (capacity < 0.0) { continue; } + + // Capacity row sum_{i in S} x_i - s <= K with every x_i bounded by a common indicator z. + members.clear(); + support.assign(indices, indices + len); + bool usable = true; + for (int p = 0; p < len && usable; ++p) { + const int col = indices[p]; + const f_t v = direction * values[p]; + if (col_flags[col].test(papilo::ColFlag::kIntegral)) { + usable = is_free_binary(col) && v == 1.0; + if (usable) { members.push_back(col); } + } else { + usable = + v < 0.0 && !col_flags[col].test(papilo::ColFlag::kLbInf) && lower_bounds[col] >= 0.0; + } + } + if (!usable || members.size() < 2) { continue; } + const f_t n_members = members.size(); + if (capacity >= n_members) { continue; } + + common.assign(vub_indicators.begin() + vub_offsets[members[0]], + vub_indicators.begin() + vub_offsets[members[0] + 1]); + for (size_t k = 1; k < members.size() && !common.empty(); ++k) { + const int col = members[k]; + intersection.clear(); + std::set_intersection(common.begin(), + common.end(), + vub_indicators.begin() + vub_offsets[col], + vub_indicators.begin() + vub_offsets[col + 1], + std::back_inserter(intersection)); + common.swap(intersection); + } + if (common.empty()) { continue; } + + std::sort(support.begin(), support.end()); + int indicator = -1; + for (int z : common) { + if (std::binary_search(support.begin(), support.end(), z)) { continue; } + indicator = z; + break; + } + if (indicator < 0) { continue; } + + entries.emplace_back(row, indicator, direction * -capacity); + if (direction == 1) { + rhs[row] = 0.0; + } else { + lhs[row] = 0.0; + } + ++num_lifted; + } + + const int num_implied = implied_heads.size(); + if (num_implied == 0 && num_lifted == 0) { return; } + + for (int k = 0; k < num_implied; ++k) { + const int row = num_rows + k; + entries.emplace_back(row, implied_heads[k], f_t{1}); + for (int p = implied_offsets[k]; p < implied_offsets[k + 1]; ++p) { + entries.emplace_back(row, implied_indicators[p], f_t{-1}); + } + lhs.push_back(0.0); + rhs.push_back(0.0); + papilo::RowFlags row_flag; + row_flag.set(papilo::RowFlag::kLhsInf); + flags.push_back(row_flag); + } + + if (!problem.getConstraintNames().empty()) { + papilo::Vec names = problem.getConstraintNames(); + for (int k = 0; k < num_implied; ++k) { + names.push_back("implied_indicator_" + std::to_string(k)); + } + problem.setConstraintNames(std::move(names)); + } + + const int num_vubs = vub_indicators.size(); + papilo::SparseStorage storage( + std::move(entries), num_rows + num_implied, num_cols, false, 4.0, 30); + problem.setConstraintMatrix(std::move(storage), std::move(lhs), std::move(rhs), std::move(flags)); + + CUOPT_LOG_INFO( + "Indicator strengthening: %d implied indicator rows added, %d capacity rows lifted over %d " + "variable upper bounds", + num_implied, + num_lifted, + num_vubs); +} + +#define INSTANTIATE(F_TYPE) template void strengthen_indicators(papilo::Problem&); + +#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT +INSTANTIATE(float) +#endif + +#if MIP_INSTANTIATE_DOUBLE +INSTANTIATE(double) +#endif + +#undef INSTANTIATE + +} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/presolve/indicator_strengthening.hpp b/cpp/src/mip_heuristics/presolve/indicator_strengthening.hpp new file mode 100644 index 0000000000..daa574196c --- /dev/null +++ b/cpp/src/mip_heuristics/presolve/indicator_strengthening.hpp @@ -0,0 +1,30 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#pragma once + +#if !defined(__clang__) +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wstringop-overflow" // ignore boost error for pip wheel build +#pragma GCC diagnostic ignored "-Wnarrowing" +#endif +#include +#include +#if !defined(__clang__) +#pragma GCC diagnostic pop +#endif + +namespace cuopt::mathematical_optimization::mip { + +// Adds an implied indicator row y <= sum_{g in D} z_g for every implication row +// y <= sum_{j in S} x_j whose members are bounded by indicators x_j <= z_g, and lifts every +// capacity row sum_{i in S} x_i - s <= K whose members share an indicator z into sum_{i in S} x_i - +// s <= K z. Both need an integral indicator, so this is for MIPs only. +template +void strengthen_indicators(papilo::Problem& problem); + +} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp index baf311889c..424e33ceaf 100644 --- a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp +++ b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp @@ -42,6 +42,7 @@ #include #include #include +#include #include #include #include @@ -922,6 +923,12 @@ third_party_presolve_status_t third_party_presolve_t::apply_papilo( { raft::common::nvtx::range fun_scope("Apply Papilo presolve on host"); + if (category == problem_category_t::MIP && + (!reduction_allowlist_.has_value() || + reduction_allowlist_->count("indicatorstrengthening") > 0)) { + strengthen_indicators(papilo_problem); + } + // Capture original dimensions before papilo.apply() mutates papilo_problem // in place into its reduced form. const i_t original_n_vars = static_cast(papilo_problem.getNCols()); diff --git a/cpp/src/mip_heuristics/solver.cu b/cpp/src/mip_heuristics/solver.cu index 2477412d8c..166afd832b 100644 --- a/cpp/src/mip_heuristics/solver.cu +++ b/cpp/src/mip_heuristics/solver.cu @@ -366,13 +366,11 @@ solution_t mip_solver_t::run_solver() } branch_and_bound_settings.mixed_integer_gomory_cuts = context.settings.mixed_integer_gomory_cuts; - branch_and_bound_settings.knapsack_cuts = context.settings.knapsack_cuts; - branch_and_bound_settings.flow_cover_cuts = context.settings.flow_cover_cuts; - branch_and_bound_settings.implied_bound_cuts = context.settings.implied_bound_cuts; - branch_and_bound_settings.implied_indicator_cuts = context.settings.implied_indicator_cuts; - branch_and_bound_settings.capacity_lifting_cuts = context.settings.capacity_lifting_cuts; - branch_and_bound_settings.clique_cuts = context.settings.clique_cuts; - branch_and_bound_settings.zero_half_cuts = context.settings.zero_half_cuts; + branch_and_bound_settings.knapsack_cuts = context.settings.knapsack_cuts; + branch_and_bound_settings.flow_cover_cuts = context.settings.flow_cover_cuts; + branch_and_bound_settings.implied_bound_cuts = context.settings.implied_bound_cuts; + branch_and_bound_settings.clique_cuts = context.settings.clique_cuts; + branch_and_bound_settings.zero_half_cuts = context.settings.zero_half_cuts; branch_and_bound_settings.strong_chvatal_gomory_cuts = context.settings.strong_chvatal_gomory_cuts; branch_and_bound_settings.cut_change_threshold = context.settings.cut_change_threshold; From 74e2f1f211e77a3acd9c4881c3eb5d7fe8a1927d Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Mon, 28 Sep 2026 14:53:57 +0200 Subject: [PATCH 07/12] simplified the code Signed-off-by: Nicolas L. Guidotti --- .../presolve/indicator_strengthening.cpp | 177 +++++++----------- 1 file changed, 69 insertions(+), 108 deletions(-) diff --git a/cpp/src/mip_heuristics/presolve/indicator_strengthening.cpp b/cpp/src/mip_heuristics/presolve/indicator_strengthening.cpp index 5939dcdad0..fa43db3a6d 100644 --- a/cpp/src/mip_heuristics/presolve/indicator_strengthening.cpp +++ b/cpp/src/mip_heuristics/presolve/indicator_strengthening.cpp @@ -11,7 +11,6 @@ #include #include -#include #include #include #include @@ -32,65 +31,52 @@ void strengthen_indicators(papilo::Problem& problem) const int num_rows = constraint_matrix.getNRows(); const int num_cols = problem.getNCols(); - if (num_rows <= 0 || num_cols <= 0) { return; } - auto is_free_binary = [&](int col) { - const auto& flags = col_flags[col]; - return flags.test(papilo::ColFlag::kIntegral) && !flags.test(papilo::ColFlag::kLbInf) && - !flags.test(papilo::ColFlag::kUbInf) && !flags.test(papilo::ColFlag::kFixed) && - lower_bounds[col] == 0.0 && upper_bounds[col] == 1.0; - }; + std::vector is_binary(num_cols); + for (int col = 0; col < num_cols; ++col) { + is_binary[col] = col_flags[col].test(papilo::ColFlag::kIntegral) && + !col_flags[col].test(papilo::ColFlag::kLbInf, papilo::ColFlag::kUbInf) && + lower_bounds[col] == 0.0 && upper_bounds[col] == 1.0; + } // +1 when the stored row reads a.x <= rhs, -1 when it reads a.x >= lhs, 0 for equations, ranges // and free rows. - auto orientation = [&](int row) { + std::vector orientation(num_rows, 0); + for (int row = 0; row < num_rows; ++row) { const bool lhs_infinite = row_flags[row].test(papilo::RowFlag::kLhsInf); const bool rhs_infinite = row_flags[row].test(papilo::RowFlag::kRhsInf); - if (lhs_infinite == rhs_infinite) { return 0; } - return lhs_infinite ? 1 : -1; - }; - - // The variable upper bounds x <= z over binaries, in CSR form over the members: the indicators - // bounding column j are vub_indicators[vub_offsets[j] .. vub_offsets[j+1]), sorted and unique. - std::vector> vubs; - for (int row = 0; row < num_rows; ++row) { - const int direction = orientation(row); - if (direction == 0) { continue; } - auto row_coefficients = constraint_matrix.getRowCoefficients(row); - if (row_coefficients.getLength() != 2) { continue; } - const f_t side = direction == 1 ? rhs_values[row] : lhs_values[row]; - if (side != 0.0) { continue; } - - const int* indices = row_coefficients.getIndices(); - const f_t* values = row_coefficients.getValues(); - int member = -1, indicator = -1; - for (int p = 0; p < 2; ++p) { - if (!is_free_binary(indices[p])) { - member = indicator = -1; - break; - } - const f_t v = direction * values[p]; - if (v == 1.0) { - member = indices[p]; - } else if (v == -1.0) { - indicator = indices[p]; - } - } - if (member >= 0 && indicator >= 0) { vubs.emplace_back(member, indicator); } + if (lhs_infinite != rhs_infinite) { orientation[row] = lhs_infinite ? 1 : -1; } } - if (vubs.empty()) { return; } - std::sort(vubs.begin(), vubs.end()); - vubs.erase(std::unique(vubs.begin(), vubs.end()), vubs.end()); + // The variable upper bounds x <= z over binaries, in CSR form over the members: the indicators + // bounding column j are vub_indicators[vub_offsets[j] .. vub_offsets[j+1]). std::vector vub_offsets(num_cols + 1, 0); - std::vector vub_indicators(vubs.size()); - for (size_t k = 0; k < vubs.size(); ++k) { - ++vub_offsets[vubs[k].first + 1]; - vub_indicators[k] = vubs[k].second; - } + std::vector vub_indicators; for (int col = 0; col < num_cols; ++col) { - vub_offsets[col + 1] += vub_offsets[col]; + if (is_binary[col]) { + auto col_coefficients = constraint_matrix.getColumnCoefficients(col); + const int* rows = col_coefficients.getIndices(); + const f_t* col_values = col_coefficients.getValues(); + for (int p = 0; p < col_coefficients.getLength(); ++p) { + const int row = rows[p]; + const int direction = orientation[row]; + if (direction == 0 || direction * col_values[p] != 1.0) { continue; } + auto row_coefficients = constraint_matrix.getRowCoefficients(row); + if (row_coefficients.getLength() != 2) { continue; } + const f_t side = direction == 1 ? rhs_values[row] : lhs_values[row]; + if (side != 0.0) { continue; } + + const int* indices = row_coefficients.getIndices(); + const f_t* values = row_coefficients.getValues(); + const int other = indices[0] == col ? 1 : 0; + if (is_binary[indices[other]] && direction * values[other] == -1.0) { + vub_indicators.push_back(indices[other]); + } + } + } + vub_offsets[col + 1] = vub_indicators.size(); } + if (vub_indicators.empty()) { return; } papilo::Vec> entries; entries.reserve(constraint_matrix.getNnz() + num_rows); @@ -98,16 +84,14 @@ void strengthen_indicators(papilo::Problem& problem) papilo::Vec rhs(rhs_values.begin(), rhs_values.end()); papilo::Vec flags(row_flags.begin(), row_flags.end()); - std::vector implied_heads; - std::vector implied_offsets{0}; - std::vector implied_indicators; - int num_lifted = 0; + papilo::Vec> implied_entries; + papilo::RowFlags implied_flags; + implied_flags.set(papilo::RowFlag::kLhsInf); + int num_implied = 0; + int num_lifted = 0; std::vector indicators; std::vector members; - std::vector support; - std::vector common; - std::vector intersection; for (int row = 0; row < num_rows; ++row) { auto row_coefficients = constraint_matrix.getRowCoefficients(row); @@ -118,8 +102,8 @@ void strengthen_indicators(papilo::Problem& problem) entries.emplace_back(row, indices[p], values[p]); } - const int direction = orientation(row); - if (direction == 0 || len < 2) { continue; } + const int direction = orientation[row]; + if (direction == 0) { continue; } const f_t capacity = direction * (direction == 1 ? rhs_values[row] : lhs_values[row]); // Implication row y - sum_{j in S} x_j <= 0: one +1 head against a tail of -1 members. @@ -131,7 +115,7 @@ void strengthen_indicators(papilo::Problem& problem) for (int p = 0; p < len && usable; ++p) { const int col = indices[p]; const f_t v = direction * values[p]; - if (!is_free_binary(col)) { + if (!is_binary[col]) { usable = false; } else if (v == 1.0) { usable = head < 0; @@ -148,7 +132,7 @@ void strengthen_indicators(papilo::Problem& problem) usable = false; } } - if (!usable || head < 0 || indicators.empty()) { continue; } + if (!usable || head < 0) { continue; } std::sort(indicators.begin(), indicators.end()); indicators.erase(std::unique(indicators.begin(), indicators.end()), indicators.end()); @@ -156,22 +140,27 @@ void strengthen_indicators(papilo::Problem& problem) if (indicators.size() >= num_members) { continue; } if (std::binary_search(indicators.begin(), indicators.end(), head)) { continue; } - implied_heads.push_back(head); - implied_indicators.insert(implied_indicators.end(), indicators.begin(), indicators.end()); - implied_offsets.push_back(implied_indicators.size()); + const int implied_row = num_rows + num_implied; + implied_entries.emplace_back(implied_row, head, f_t{1}); + for (int z : indicators) { + implied_entries.emplace_back(implied_row, z, f_t{-1}); + } + lhs.push_back(0.0); + rhs.push_back(0.0); + flags.push_back(implied_flags); + ++num_implied; continue; } if (capacity < 0.0) { continue; } // Capacity row sum_{i in S} x_i - s <= K with every x_i bounded by a common indicator z. members.clear(); - support.assign(indices, indices + len); bool usable = true; for (int p = 0; p < len && usable; ++p) { const int col = indices[p]; const f_t v = direction * values[p]; if (col_flags[col].test(papilo::ColFlag::kIntegral)) { - usable = is_free_binary(col) && v == 1.0; + usable = is_binary[col] && v == 1.0; if (usable) { members.push_back(col); } } else { usable = @@ -182,26 +171,16 @@ void strengthen_indicators(papilo::Problem& problem) const f_t n_members = members.size(); if (capacity >= n_members) { continue; } - common.assign(vub_indicators.begin() + vub_offsets[members[0]], - vub_indicators.begin() + vub_offsets[members[0] + 1]); - for (size_t k = 1; k < members.size() && !common.empty(); ++k) { - const int col = members[k]; - intersection.clear(); - std::set_intersection(common.begin(), - common.end(), - vub_indicators.begin() + vub_offsets[col], - vub_indicators.begin() + vub_offsets[col + 1], - std::back_inserter(intersection)); - common.swap(intersection); - } - if (common.empty()) { continue; } - - std::sort(support.begin(), support.end()); int indicator = -1; - for (int z : common) { - if (std::binary_search(support.begin(), support.end(), z)) { continue; } - indicator = z; - break; + for (int p = vub_offsets[members[0]]; p < vub_offsets[members[0] + 1] && indicator < 0; ++p) { + const int z = vub_indicators[p]; + bool shared = true; + for (size_t k = 1; k < members.size() && shared; ++k) { + const auto span_begin = vub_indicators.begin() + vub_offsets[members[k]]; + const auto span_end = vub_indicators.begin() + vub_offsets[members[k] + 1]; + shared = std::find(span_begin, span_end, z) != span_end; + } + if (shared) { indicator = z; } } if (indicator < 0) { continue; } @@ -214,21 +193,8 @@ void strengthen_indicators(papilo::Problem& problem) ++num_lifted; } - const int num_implied = implied_heads.size(); if (num_implied == 0 && num_lifted == 0) { return; } - - for (int k = 0; k < num_implied; ++k) { - const int row = num_rows + k; - entries.emplace_back(row, implied_heads[k], f_t{1}); - for (int p = implied_offsets[k]; p < implied_offsets[k + 1]; ++p) { - entries.emplace_back(row, implied_indicators[p], f_t{-1}); - } - lhs.push_back(0.0); - rhs.push_back(0.0); - papilo::RowFlags row_flag; - row_flag.set(papilo::RowFlag::kLhsInf); - flags.push_back(row_flag); - } + entries.insert(entries.end(), implied_entries.begin(), implied_entries.end()); if (!problem.getConstraintNames().empty()) { papilo::Vec names = problem.getConstraintNames(); @@ -238,29 +204,24 @@ void strengthen_indicators(papilo::Problem& problem) problem.setConstraintNames(std::move(names)); } - const int num_vubs = vub_indicators.size(); papilo::SparseStorage storage( std::move(entries), num_rows + num_implied, num_cols, false, 4.0, 30); problem.setConstraintMatrix(std::move(storage), std::move(lhs), std::move(rhs), std::move(flags)); - CUOPT_LOG_INFO( - "Indicator strengthening: %d implied indicator rows added, %d capacity rows lifted over %d " + CUOPT_LOG_DEBUG( + "Indicator strengthening: %d implied indicator rows added, %d capacity rows lifted over %zu " "variable upper bounds", num_implied, num_lifted, - num_vubs); + vub_indicators.size()); } -#define INSTANTIATE(F_TYPE) template void strengthen_indicators(papilo::Problem&); - #if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT -INSTANTIATE(float) +template void strengthen_indicators(papilo::Problem&); #endif #if MIP_INSTANTIATE_DOUBLE -INSTANTIATE(double) +template void strengthen_indicators(papilo::Problem&); #endif -#undef INSTANTIATE - } // namespace cuopt::mathematical_optimization::mip From e69fcda695bd1388fc43a11a1e68cfd593437ba3 Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Tue, 29 Sep 2026 16:57:50 +0200 Subject: [PATCH 08/12] refactor presolve to use an object-oriented style. separated the two reductions. Signed-off-by: Nicolas L. Guidotti --- .../presolve/indicator_strengthening.cpp | 405 ++++++++++++------ .../presolve/indicator_strengthening.hpp | 2 +- .../presolve/third_party_presolve.cpp | 2 +- 3 files changed, 281 insertions(+), 128 deletions(-) diff --git a/cpp/src/mip_heuristics/presolve/indicator_strengthening.cpp b/cpp/src/mip_heuristics/presolve/indicator_strengthening.cpp index fa43db3a6d..c5f172c348 100644 --- a/cpp/src/mip_heuristics/presolve/indicator_strengthening.cpp +++ b/cpp/src/mip_heuristics/presolve/indicator_strengthening.cpp @@ -12,156 +12,227 @@ #include #include +#include #include #include namespace cuopt::mathematical_optimization::mip { -template -void strengthen_indicators(papilo::Problem& problem) +namespace { + +template +class indicator_strengthening_t { + public: + indicator_strengthening_t(const papilo::Problem& problem); + + // Implication row y - sum_{j in S} x_j <= 0: one +1 head against a tail of -1 members. + i_t add_implied_indicator_rows(papilo::Vec>& implied_entries) const; + + // Capacity row sum_{i in S} x_i - s <= K with every x_i bounded by a common indicator z. + i_t lift_capacity_rows(papilo::Vec>& lifted_entries) const; + + i_t num_vubs() const { return vub_indicators.size(); } + i_t row_orientation(i_t row) const { return orientation[row]; } + + private: + const papilo::Problem& problem; + + std::vector is_binary; + + // +1 when the stored row reads a.x <= rhs, -1 when it reads a.x >= lhs, 0 for equations, ranges + // and free rows. + std::vector orientation; + // The variable upper bounds x <= z over binaries, in CSR form over the members: the indicators + // bounding column j are vub_indicators[vub_offsets[j] .. vub_offsets[j+1]). + std::vector vub_offsets; + std::vector vub_indicators; + + i_t max_row_size = 0; +}; + +template +indicator_strengthening_t::indicator_strengthening_t(const papilo::Problem& problem) + : problem(problem) { const auto& constraint_matrix = problem.getConstraintMatrix(); const auto& lhs_values = constraint_matrix.getLeftHandSides(); const auto& rhs_values = constraint_matrix.getRightHandSides(); const auto& row_flags = constraint_matrix.getRowFlags(); + const auto& row_sizes = constraint_matrix.getRowSizes(); const auto& domains = problem.getVariableDomains(); const auto& col_flags = domains.flags; const auto& lower_bounds = domains.lower_bounds; const auto& upper_bounds = domains.upper_bounds; - const int num_rows = constraint_matrix.getNRows(); - const int num_cols = problem.getNCols(); + const i_t num_rows = constraint_matrix.getNRows(); + const i_t num_cols = problem.getNCols(); - std::vector is_binary(num_cols); - for (int col = 0; col < num_cols; ++col) { + is_binary.resize(num_cols); + for (i_t col = 0; col < num_cols; ++col) { is_binary[col] = col_flags[col].test(papilo::ColFlag::kIntegral) && !col_flags[col].test(papilo::ColFlag::kLbInf, papilo::ColFlag::kUbInf) && lower_bounds[col] == 0.0 && upper_bounds[col] == 1.0; } - // +1 when the stored row reads a.x <= rhs, -1 when it reads a.x >= lhs, 0 for equations, ranges - // and free rows. - std::vector orientation(num_rows, 0); - for (int row = 0; row < num_rows; ++row) { + orientation.assign(num_rows, 0); + std::vector> vubs; + vubs.reserve(std::count(row_sizes.begin(), row_sizes.end(), 2)); + for (i_t row = 0; row < num_rows; ++row) { + max_row_size = std::max(max_row_size, row_sizes[row]); const bool lhs_infinite = row_flags[row].test(papilo::RowFlag::kLhsInf); const bool rhs_infinite = row_flags[row].test(papilo::RowFlag::kRhsInf); - if (lhs_infinite != rhs_infinite) { orientation[row] = lhs_infinite ? 1 : -1; } - } + if (lhs_infinite == rhs_infinite) { continue; } + const i_t direction = lhs_infinite ? 1 : -1; + orientation[row] = direction; + if (row_sizes[row] != 2) { continue; } + const f_t side = direction == 1 ? rhs_values[row] : lhs_values[row]; + if (side != 0.0) { continue; } - // The variable upper bounds x <= z over binaries, in CSR form over the members: the indicators - // bounding column j are vub_indicators[vub_offsets[j] .. vub_offsets[j+1]). - std::vector vub_offsets(num_cols + 1, 0); - std::vector vub_indicators; - for (int col = 0; col < num_cols; ++col) { - if (is_binary[col]) { - auto col_coefficients = constraint_matrix.getColumnCoefficients(col); - const int* rows = col_coefficients.getIndices(); - const f_t* col_values = col_coefficients.getValues(); - for (int p = 0; p < col_coefficients.getLength(); ++p) { - const int row = rows[p]; - const int direction = orientation[row]; - if (direction == 0 || direction * col_values[p] != 1.0) { continue; } - auto row_coefficients = constraint_matrix.getRowCoefficients(row); - if (row_coefficients.getLength() != 2) { continue; } - const f_t side = direction == 1 ? rhs_values[row] : lhs_values[row]; - if (side != 0.0) { continue; } - - const int* indices = row_coefficients.getIndices(); - const f_t* values = row_coefficients.getValues(); - const int other = indices[0] == col ? 1 : 0; - if (is_binary[indices[other]] && direction * values[other] == -1.0) { - vub_indicators.push_back(indices[other]); - } - } + auto row_coefficients = constraint_matrix.getRowCoefficients(row); + const i_t* indices = row_coefficients.getIndices(); + const f_t* values = row_coefficients.getValues(); + if (!is_binary[indices[0]] || !is_binary[indices[1]]) { continue; } + const f_t v0 = direction * values[0]; + const f_t v1 = direction * values[1]; + if (v0 == 1.0 && v1 == -1.0) { + vubs.emplace_back(indices[0], indices[1]); + } else if (v0 == -1.0 && v1 == 1.0) { + vubs.emplace_back(indices[1], indices[0]); } - vub_offsets[col + 1] = vub_indicators.size(); } - if (vub_indicators.empty()) { return; } - papilo::Vec> entries; - entries.reserve(constraint_matrix.getNnz() + num_rows); - papilo::Vec lhs(lhs_values.begin(), lhs_values.end()); - papilo::Vec rhs(rhs_values.begin(), rhs_values.end()); - papilo::Vec flags(row_flags.begin(), row_flags.end()); + vub_offsets.assign(num_cols + 1, 0); + for (const auto& vub : vubs) { + ++vub_offsets[vub.first + 1]; + } + for (i_t col = 0; col < num_cols; ++col) { + vub_offsets[col + 1] += vub_offsets[col]; + } + vub_indicators.resize(vubs.size()); + std::vector next(vub_offsets.begin(), vub_offsets.end() - 1); + for (const auto& [member, indicator] : vubs) { + vub_indicators[next[member]++] = indicator; + } +} - papilo::Vec> implied_entries; - papilo::RowFlags implied_flags; - implied_flags.set(papilo::RowFlag::kLhsInf); - int num_implied = 0; - int num_lifted = 0; +template +i_t indicator_strengthening_t::add_implied_indicator_rows( + papilo::Vec>& implied_entries) const +{ + const auto& constraint_matrix = problem.getConstraintMatrix(); + const auto& lhs_values = constraint_matrix.getLeftHandSides(); + const auto& rhs_values = constraint_matrix.getRightHandSides(); + const auto& row_sizes = constraint_matrix.getRowSizes(); + + const i_t num_rows = constraint_matrix.getNRows(); + const i_t num_cols = problem.getNCols(); + + i_t num_implied = 0; + + std::vector indicators; + indicators.reserve(max_row_size); + std::vector mark(num_cols, -1); - std::vector indicators; - std::vector members; + for (i_t row = 0; row < num_rows; ++row) { + const i_t direction = orientation[row]; + if (direction == 0 || row_sizes[row] < 3) { continue; } + const f_t side = direction == 1 ? rhs_values[row] : lhs_values[row]; + if (side != 0.0) { continue; } - for (int row = 0; row < num_rows; ++row) { auto row_coefficients = constraint_matrix.getRowCoefficients(row); - const int len = row_coefficients.getLength(); - const int* indices = row_coefficients.getIndices(); + const i_t len = row_coefficients.getLength(); + const i_t* indices = row_coefficients.getIndices(); const f_t* values = row_coefficients.getValues(); - for (int p = 0; p < len; ++p) { - entries.emplace_back(row, indices[p], values[p]); - } - const int direction = orientation[row]; - if (direction == 0) { continue; } - const f_t capacity = direction * (direction == 1 ? rhs_values[row] : lhs_values[row]); + const size_t num_members = len - 1; + i_t head = -1; + bool usable = true; + indicators.clear(); + for (i_t p = 0; p < len; ++p) { + const i_t col = indices[p]; + const f_t v = direction * values[p]; + if (!is_binary[col]) { + usable = false; + break; + } - // Implication row y - sum_{j in S} x_j <= 0: one +1 head against a tail of -1 members. - if (capacity == 0.0) { - if (len < 3) { continue; } - int head = -1; - bool usable = true; - indicators.clear(); - for (int p = 0; p < len && usable; ++p) { - const int col = indices[p]; - const f_t v = direction * values[p]; - if (!is_binary[col]) { - usable = false; - } else if (v == 1.0) { - usable = head < 0; - head = col; - } else if (v == -1.0) { - if (vub_offsets[col] == vub_offsets[col + 1]) { - indicators.push_back(col); - } else { - indicators.insert(indicators.end(), - vub_indicators.begin() + vub_offsets[col], - vub_indicators.begin() + vub_offsets[col + 1]); - } - } else { - usable = false; + if (v == 1.0) { + usable = head < 0; + head = col; + } else if (v == -1.0) { + const bool bounded = vub_offsets[col] < vub_offsets[col + 1]; + const i_t* first = bounded ? vub_indicators.data() + vub_offsets[col] : &col; + const i_t* last = bounded ? vub_indicators.data() + vub_offsets[col + 1] : &col + 1; + for (const i_t* z = first; z != last && usable; ++z) { + if (mark[*z] == row) { continue; } + mark[*z] = row; + indicators.push_back(*z); + usable = indicators.size() < num_members; } + } else { + usable = false; + break; } - if (!usable || head < 0) { continue; } - - std::sort(indicators.begin(), indicators.end()); - indicators.erase(std::unique(indicators.begin(), indicators.end()), indicators.end()); - const size_t num_members = len - 1; - if (indicators.size() >= num_members) { continue; } - if (std::binary_search(indicators.begin(), indicators.end(), head)) { continue; } - - const int implied_row = num_rows + num_implied; - implied_entries.emplace_back(implied_row, head, f_t{1}); - for (int z : indicators) { - implied_entries.emplace_back(implied_row, z, f_t{-1}); - } - lhs.push_back(0.0); - rhs.push_back(0.0); - flags.push_back(implied_flags); - ++num_implied; - continue; } - if (capacity < 0.0) { continue; } + if (!usable || head < 0 || mark[head] == row) { continue; } + + std::sort(indicators.begin(), indicators.end()); + const i_t implied_row = num_rows + num_implied; + const auto split = std::lower_bound(indicators.begin(), indicators.end(), head); + for (auto it = indicators.begin(); it != split; ++it) { + implied_entries.emplace_back(implied_row, *it, f_t{-1}); + } + implied_entries.emplace_back(implied_row, head, f_t{1}); + for (auto it = split; it != indicators.end(); ++it) { + implied_entries.emplace_back(implied_row, *it, f_t{-1}); + } + ++num_implied; + } + + return num_implied; +} + +template +i_t indicator_strengthening_t::lift_capacity_rows( + papilo::Vec>& lifted_entries) const +{ + const auto& constraint_matrix = problem.getConstraintMatrix(); + const auto& lhs_values = constraint_matrix.getLeftHandSides(); + const auto& rhs_values = constraint_matrix.getRightHandSides(); + const auto& domains = problem.getVariableDomains(); + const auto& col_flags = domains.flags; + const auto& lower_bounds = domains.lower_bounds; + + const i_t num_rows = constraint_matrix.getNRows(); + + i_t num_lifted = 0; + + std::vector members; + members.reserve(max_row_size); + + for (i_t row = 0; row < num_rows; ++row) { + const i_t direction = orientation[row]; + if (direction == 0) { continue; } + const f_t capacity = direction * (direction == 1 ? rhs_values[row] : lhs_values[row]); + if (capacity <= 0.0) { continue; } + + auto row_coefficients = constraint_matrix.getRowCoefficients(row); + const i_t len = row_coefficients.getLength(); + const i_t* indices = row_coefficients.getIndices(); + const f_t* values = row_coefficients.getValues(); - // Capacity row sum_{i in S} x_i - s <= K with every x_i bounded by a common indicator z. members.clear(); + i_t pivot = -1; bool usable = true; - for (int p = 0; p < len && usable; ++p) { - const int col = indices[p]; + for (i_t p = 0; p < len && usable; ++p) { + const i_t col = indices[p]; const f_t v = direction * values[p]; if (col_flags[col].test(papilo::ColFlag::kIntegral)) { - usable = is_binary[col] && v == 1.0; - if (usable) { members.push_back(col); } + const i_t num_vubs = vub_offsets[col + 1] - vub_offsets[col]; + usable = is_binary[col] && v == 1.0 && num_vubs > 0; + if (!usable) { continue; } + members.push_back(col); + if (pivot < 0 || num_vubs < vub_offsets[pivot + 1] - vub_offsets[pivot]) { pivot = col; } } else { usable = v < 0.0 && !col_flags[col].test(papilo::ColFlag::kLbInf) && lower_bounds[col] >= 0.0; @@ -171,11 +242,12 @@ void strengthen_indicators(papilo::Problem& problem) const f_t n_members = members.size(); if (capacity >= n_members) { continue; } - int indicator = -1; - for (int p = vub_offsets[members[0]]; p < vub_offsets[members[0] + 1] && indicator < 0; ++p) { - const int z = vub_indicators[p]; + i_t indicator = -1; + for (i_t p = vub_offsets[pivot]; p < vub_offsets[pivot + 1] && indicator < 0; ++p) { + const i_t z = vub_indicators[p]; bool shared = true; - for (size_t k = 1; k < members.size() && shared; ++k) { + for (size_t k = 0; k < members.size() && shared; ++k) { + if (members[k] == pivot) { continue; } const auto span_begin = vub_indicators.begin() + vub_offsets[members[k]]; const auto span_end = vub_indicators.begin() + vub_offsets[members[k] + 1]; shared = std::find(span_begin, span_end, z) != span_end; @@ -184,44 +256,125 @@ void strengthen_indicators(papilo::Problem& problem) } if (indicator < 0) { continue; } - entries.emplace_back(row, indicator, direction * -capacity); - if (direction == 1) { - rhs[row] = 0.0; - } else { - lhs[row] = 0.0; - } + lifted_entries.emplace_back(row, indicator, direction * -capacity); ++num_lifted; } + return num_lifted; +} + +} // namespace + +template +void strengthen_indicators(papilo::Problem& problem) +{ + const indicator_strengthening_t strengthening(problem); + if (strengthening.num_vubs() == 0) { return; } + + const auto& constraint_matrix = problem.getConstraintMatrix(); + const auto& lhs_values = constraint_matrix.getLeftHandSides(); + const auto& rhs_values = constraint_matrix.getRightHandSides(); + const auto& row_flags = constraint_matrix.getRowFlags(); + const auto& row_sizes = constraint_matrix.getRowSizes(); + + const i_t num_rows = constraint_matrix.getNRows(); + const i_t num_cols = problem.getNCols(); + + i_t max_implied_entries = 0; + i_t max_lifted_rows = 0; + for (i_t row = 0; row < num_rows; ++row) { + const i_t direction = strengthening.row_orientation(row); + if (direction == 0) { continue; } + const f_t capacity = direction * (direction == 1 ? rhs_values[row] : lhs_values[row]); + if (capacity == 0.0 && row_sizes[row] >= 3) { + max_implied_entries += row_sizes[row] - 1; + } else if (capacity > 0.0) { + ++max_lifted_rows; + } + } + + papilo::Vec> implied_entries; + papilo::Vec> lifted_entries; + implied_entries.reserve(max_implied_entries); + lifted_entries.reserve(max_lifted_rows); + const i_t num_implied = strengthening.add_implied_indicator_rows(implied_entries); + const i_t num_lifted = strengthening.lift_capacity_rows(lifted_entries); if (num_implied == 0 && num_lifted == 0) { return; } + + papilo::Vec> entries; + entries.reserve(constraint_matrix.getNnz() + lifted_entries.size() + implied_entries.size()); + auto lifted = lifted_entries.begin(); + for (i_t row = 0; row < num_rows; ++row) { + auto row_coefficients = constraint_matrix.getRowCoefficients(row); + const i_t len = row_coefficients.getLength(); + const i_t* indices = row_coefficients.getIndices(); + const f_t* values = row_coefficients.getValues(); + i_t p = 0; + if (lifted != lifted_entries.end() && std::get<0>(*lifted) == row) { + for (; p < len && indices[p] < std::get<1>(*lifted); ++p) { + entries.emplace_back(row, indices[p], values[p]); + } + entries.push_back(*lifted); + ++lifted; + } + for (; p < len; ++p) { + entries.emplace_back(row, indices[p], values[p]); + } + } entries.insert(entries.end(), implied_entries.begin(), implied_entries.end()); - if (!problem.getConstraintNames().empty()) { - papilo::Vec names = problem.getConstraintNames(); - for (int k = 0; k < num_implied; ++k) { + papilo::Vec lhs; + papilo::Vec rhs; + papilo::Vec flags; + lhs.reserve(num_rows + num_implied); + rhs.reserve(num_rows + num_implied); + flags.reserve(num_rows + num_implied); + lhs.assign(lhs_values.begin(), lhs_values.end()); + rhs.assign(rhs_values.begin(), rhs_values.end()); + flags.assign(row_flags.begin(), row_flags.end()); + for (const auto& entry : lifted_entries) { + const i_t row = std::get<0>(entry); + if (strengthening.row_orientation(row) == 1) { + rhs[row] = 0.0; + } else { + lhs[row] = 0.0; + } + } + papilo::RowFlags implied_flags; + implied_flags.set(papilo::RowFlag::kLhsInf); + lhs.resize(num_rows + num_implied, 0.0); + rhs.resize(num_rows + num_implied, 0.0); + flags.resize(num_rows + num_implied, implied_flags); + + const auto& constraint_names = problem.getConstraintNames(); + if (!constraint_names.empty()) { + papilo::Vec names; + names.reserve(constraint_names.size() + num_implied); + names.assign(constraint_names.begin(), constraint_names.end()); + for (i_t k = 0; k < num_implied; ++k) { names.push_back("implied_indicator_" + std::to_string(k)); } problem.setConstraintNames(std::move(names)); } papilo::SparseStorage storage( - std::move(entries), num_rows + num_implied, num_cols, false, 4.0, 30); + std::move(entries), num_rows + num_implied, num_cols, true, 4.0, 30); problem.setConstraintMatrix(std::move(storage), std::move(lhs), std::move(rhs), std::move(flags)); CUOPT_LOG_DEBUG( - "Indicator strengthening: %d implied indicator rows added, %d capacity rows lifted over %zu " + "Indicator strengthening: %d implied indicator rows added, %d capacity rows lifted over %d " "variable upper bounds", num_implied, num_lifted, - vub_indicators.size()); + strengthening.num_vubs()); } #if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT -template void strengthen_indicators(papilo::Problem&); +template void strengthen_indicators(papilo::Problem&); #endif #if MIP_INSTANTIATE_DOUBLE -template void strengthen_indicators(papilo::Problem&); +template void strengthen_indicators(papilo::Problem&); #endif } // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/presolve/indicator_strengthening.hpp b/cpp/src/mip_heuristics/presolve/indicator_strengthening.hpp index daa574196c..203447f2df 100644 --- a/cpp/src/mip_heuristics/presolve/indicator_strengthening.hpp +++ b/cpp/src/mip_heuristics/presolve/indicator_strengthening.hpp @@ -24,7 +24,7 @@ namespace cuopt::mathematical_optimization::mip { // y <= sum_{j in S} x_j whose members are bounded by indicators x_j <= z_g, and lifts every // capacity row sum_{i in S} x_i - s <= K whose members share an indicator z into sum_{i in S} x_i - // s <= K z. Both need an integral indicator, so this is for MIPs only. -template +template void strengthen_indicators(papilo::Problem& problem); } // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp index 424e33ceaf..a4389c6337 100644 --- a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp +++ b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp @@ -926,7 +926,7 @@ third_party_presolve_status_t third_party_presolve_t::apply_papilo( if (category == problem_category_t::MIP && (!reduction_allowlist_.has_value() || reduction_allowlist_->count("indicatorstrengthening") > 0)) { - strengthen_indicators(papilo_problem); + strengthen_indicators(papilo_problem); } // Capture original dimensions before papilo.apply() mutates papilo_problem From 3e2acdad4fd511b3f03a2616c4b80bd0e83cebcb Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Tue, 29 Sep 2026 18:07:10 +0200 Subject: [PATCH 09/12] addressed code rabbit Signed-off-by: Nicolas L. Guidotti --- .../mip_heuristics/presolve/indicator_strengthening.cpp | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/cpp/src/mip_heuristics/presolve/indicator_strengthening.cpp b/cpp/src/mip_heuristics/presolve/indicator_strengthening.cpp index c5f172c348..b55eff2f09 100644 --- a/cpp/src/mip_heuristics/presolve/indicator_strengthening.cpp +++ b/cpp/src/mip_heuristics/presolve/indicator_strengthening.cpp @@ -148,15 +148,12 @@ i_t indicator_strengthening_t::add_implied_indicator_rows( i_t head = -1; bool usable = true; indicators.clear(); - for (i_t p = 0; p < len; ++p) { + for (i_t p = 0; p < len && usable; ++p) { const i_t col = indices[p]; const f_t v = direction * values[p]; if (!is_binary[col]) { usable = false; - break; - } - - if (v == 1.0) { + } else if (v == 1.0) { usable = head < 0; head = col; } else if (v == -1.0) { @@ -171,7 +168,6 @@ i_t indicator_strengthening_t::add_implied_indicator_rows( } } else { usable = false; - break; } } if (!usable || head < 0 || mark[head] == row) { continue; } From 3a618d15da206b300cd8fb9f7a28f6cefabe8f11 Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Thu, 1 Oct 2026 15:20:33 +0200 Subject: [PATCH 10/12] added hyper parameter and unit test. added additional comments. Signed-off-by: Nicolas L. Guidotti --- .../mathematical_optimization/constants.h | 4 + .../mip/solver_settings.hpp | 7 + cpp/src/math_optimization/solver_settings.cu | 1 + .../presolve/indicator_strengthening.cpp | 95 +++-- .../presolve/third_party_presolve.cpp | 19 +- .../presolve/third_party_presolve.hpp | 3 + cpp/src/mip_heuristics/solve.cu | 3 +- cpp/tests/internal/CMakeLists.txt | 1 + .../mip/indicator_strengthening_test.cpp | 326 ++++++++++++++++++ 9 files changed, 413 insertions(+), 46 deletions(-) create mode 100644 cpp/tests/mip/indicator_strengthening_test.cpp diff --git a/cpp/include/cuopt/mathematical_optimization/constants.h b/cpp/include/cuopt/mathematical_optimization/constants.h index 01afb9f98f..8fc0784899 100644 --- a/cpp/include/cuopt/mathematical_optimization/constants.h +++ b/cpp/include/cuopt/mathematical_optimization/constants.h @@ -158,6 +158,10 @@ /* @brief Block bounded-variable-elimination step of cuOpt's internal MIP presolve */ #define CUOPT_MIP_HYPER_BLOCK_BVE "mip_hyper_block_bve" +/* @brief Indicator-strengthening step that runs before Papilo presolve on MIPs */ +#define CUOPT_MIP_HYPER_PRESOLVE_INDICATOR_STRENGTHENING \ + "mip_hyper_presolve_indicator_strengthening" + /* @brief QCQP (barrier) scaling hyper-parameters */ #define CUOPT_QCQP_HYPER_RUIZ_EQUILIBRATION "qcqp_hyper_ruiz_equilibration" diff --git a/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp b/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp index f3dc5fe340..07204846e6 100644 --- a/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp +++ b/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp @@ -176,6 +176,13 @@ class mip_solver_settings_t { * no-op when no certified reduction exists. */ bool block_bve{true}; + /** + * @brief Enable the indicator-strengthening step of presolve (MIP only). + * + * Runs before Papilo and only when the higher-level presolve is enabled. It appends implied + * indicator rows and lifts capacity rows by the indicator that bounds all of their members. + */ + bool indicator_strengthening{true}; /** * @brief Determinism mode for MIP solver. * diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index 49a3e4fd2b..dff8f4c6cb 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -264,6 +264,7 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings // Recursive sub-MIP (RINS) hyper-parameters (hidden from default --help: name contains "hyper_") {CUOPT_MIP_HYPER_SUBMIP_ENABLE_CPUFJ, &mip_settings.submip_params.enable_cpufj, true, "run CPU FJ over the sub-MIP"}, {CUOPT_MIP_HYPER_BLOCK_BVE, &mip_settings.block_bve, true, "eliminate blocks of binaries in cuOpt's MIP presolve (needs " CUOPT_MIP_PROBING ")"}, + {CUOPT_MIP_HYPER_PRESOLVE_INDICATOR_STRENGTHENING, &mip_settings.indicator_strengthening, true, "append implied indicator rows and lift capacity rows before Papilo presolve"}, // PDLP scaling hyper-parameter (hidden from default --help: name contains "hyper_") {CUOPT_PDLP_HYPER_ENABLE_CURTIS_REID_SCALING, &pdlp_settings.hyper_params.do_curtis_reid_scaling, true, "Curtis-Reid prescaling, run before Ruiz/Pock-Chambolle scaling"}, }; diff --git a/cpp/src/mip_heuristics/presolve/indicator_strengthening.cpp b/cpp/src/mip_heuristics/presolve/indicator_strengthening.cpp index b55eff2f09..6ef2c31e86 100644 --- a/cpp/src/mip_heuristics/presolve/indicator_strengthening.cpp +++ b/cpp/src/mip_heuristics/presolve/indicator_strengthening.cpp @@ -20,6 +20,12 @@ namespace cuopt::mathematical_optimization::mip { namespace { +// Indicator variables $z$ are binaries that must be paid for before anything they own may be used, +// while member variables $x_j$ are the variables owned by a given indicator. +// They are linked via the following constraints: +// Link -> x_j - z <= 0 +// Disjunction -> y - \sum_{j \in \mathcal{S}} x_j <= 0 +// Capacity -> \sum_{i \in \mathcal{S}} x_i - s <= K template class indicator_strengthening_t { public: @@ -31,7 +37,7 @@ class indicator_strengthening_t { // Capacity row sum_{i in S} x_i - s <= K with every x_i bounded by a common indicator z. i_t lift_capacity_rows(papilo::Vec>& lifted_entries) const; - i_t num_vubs() const { return vub_indicators.size(); } + i_t num_variable_upper_bounds() const { return variable_upper_bound_indicators.size(); } i_t row_orientation(i_t row) const { return orientation[row]; } private: @@ -39,13 +45,17 @@ class indicator_strengthening_t { std::vector is_binary; - // +1 when the stored row reads a.x <= rhs, -1 when it reads a.x >= lhs, 0 for equations, ranges - // and free rows. + // +1 when the stored row reads A[i, :]^T x <= rhs, -1 when it reads A[i, :]^T x >= lhs, 0 for + // equations, ranges and free rows. std::vector orientation; - // The variable upper bounds x <= z over binaries, in CSR form over the members: the indicators - // bounding column j are vub_indicators[vub_offsets[j] .. vub_offsets[j+1]). - std::vector vub_offsets; - std::vector vub_indicators; + + // CSR storing the implication graph for binary variables of the form x <= z, where x is a + // member variable and z, an indicator. Each row of the CSR corresponds to a member variable + // x and each column the indicator z that owns x. Here the coefficients of the CSR are irrelevant, + // so we just store the offsets (`variable_upper_bound_offsets`) and the nonzero columns + // (`variable_upper_bound_indicators`). + std::vector variable_upper_bound_offsets; + std::vector variable_upper_bound_indicators; i_t max_row_size = 0; }; @@ -75,8 +85,8 @@ indicator_strengthening_t::indicator_strengthening_t(const papilo::Pro } orientation.assign(num_rows, 0); - std::vector> vubs; - vubs.reserve(std::count(row_sizes.begin(), row_sizes.end(), 2)); + std::vector> variable_upper_bounds; + variable_upper_bounds.reserve(std::count(row_sizes.begin(), row_sizes.end(), 2)); for (i_t row = 0; row < num_rows; ++row) { max_row_size = std::max(max_row_size, row_sizes[row]); const bool lhs_infinite = row_flags[row].test(papilo::RowFlag::kLhsInf); @@ -95,23 +105,24 @@ indicator_strengthening_t::indicator_strengthening_t(const papilo::Pro const f_t v0 = direction * values[0]; const f_t v1 = direction * values[1]; if (v0 == 1.0 && v1 == -1.0) { - vubs.emplace_back(indices[0], indices[1]); + variable_upper_bounds.emplace_back(indices[0], indices[1]); } else if (v0 == -1.0 && v1 == 1.0) { - vubs.emplace_back(indices[1], indices[0]); + variable_upper_bounds.emplace_back(indices[1], indices[0]); } } - vub_offsets.assign(num_cols + 1, 0); - for (const auto& vub : vubs) { - ++vub_offsets[vub.first + 1]; + variable_upper_bound_offsets.assign(num_cols + 1, 0); + for (const auto& variable_upper_bound : variable_upper_bounds) { + ++variable_upper_bound_offsets[variable_upper_bound.first + 1]; } for (i_t col = 0; col < num_cols; ++col) { - vub_offsets[col + 1] += vub_offsets[col]; + variable_upper_bound_offsets[col + 1] += variable_upper_bound_offsets[col]; } - vub_indicators.resize(vubs.size()); - std::vector next(vub_offsets.begin(), vub_offsets.end() - 1); - for (const auto& [member, indicator] : vubs) { - vub_indicators[next[member]++] = indicator; + variable_upper_bound_indicators.resize(variable_upper_bounds.size()); + std::vector next(variable_upper_bound_offsets.begin(), + variable_upper_bound_offsets.end() - 1); + for (const auto& [member, indicator] : variable_upper_bounds) { + variable_upper_bound_indicators[next[member]++] = indicator; } } @@ -157,13 +168,14 @@ i_t indicator_strengthening_t::add_implied_indicator_rows( usable = head < 0; head = col; } else if (v == -1.0) { - const bool bounded = vub_offsets[col] < vub_offsets[col + 1]; - const i_t* first = bounded ? vub_indicators.data() + vub_offsets[col] : &col; - const i_t* last = bounded ? vub_indicators.data() + vub_offsets[col + 1] : &col + 1; - for (const i_t* z = first; z != last && usable; ++z) { - if (mark[*z] == row) { continue; } - mark[*z] = row; - indicators.push_back(*z); + const i_t start = variable_upper_bound_offsets[col]; + const i_t end = variable_upper_bound_offsets[col + 1]; + const i_t count = std::max(end - start, 1); + for (i_t k = 0; k < count && usable; ++k) { + const i_t z = start < end ? variable_upper_bound_indicators[start + k] : col; + if (mark[z] == row) { continue; } + mark[z] = row; + indicators.push_back(z); usable = indicators.size() < num_members; } } else { @@ -172,6 +184,8 @@ i_t indicator_strengthening_t::add_implied_indicator_rows( } if (!usable || head < 0 || mark[head] == row) { continue; } + // In the previous loop we insert the indicators in a random order, however, Papilo + // expects the triplets (row, col, val) to be sorted by row and then by column. std::sort(indicators.begin(), indicators.end()); const i_t implied_row = num_rows + num_implied; const auto split = std::lower_bound(indicators.begin(), indicators.end(), head); @@ -224,11 +238,16 @@ i_t indicator_strengthening_t::lift_capacity_rows( const i_t col = indices[p]; const f_t v = direction * values[p]; if (col_flags[col].test(papilo::ColFlag::kIntegral)) { - const i_t num_vubs = vub_offsets[col + 1] - vub_offsets[col]; - usable = is_binary[col] && v == 1.0 && num_vubs > 0; + const i_t num_indicators = + variable_upper_bound_offsets[col + 1] - variable_upper_bound_offsets[col]; + usable = is_binary[col] && v == 1.0 && num_indicators > 0; if (!usable) { continue; } members.push_back(col); - if (pivot < 0 || num_vubs < vub_offsets[pivot + 1] - vub_offsets[pivot]) { pivot = col; } + + const i_t num_indicators_pivot = + pivot >= 0 ? variable_upper_bound_offsets[pivot + 1] - variable_upper_bound_offsets[pivot] + : 0; + if (pivot < 0 || num_indicators < num_indicators_pivot) { pivot = col; } } else { usable = v < 0.0 && !col_flags[col].test(papilo::ColFlag::kLbInf) && lower_bounds[col] >= 0.0; @@ -239,14 +258,18 @@ i_t indicator_strengthening_t::lift_capacity_rows( if (capacity >= n_members) { continue; } i_t indicator = -1; - for (i_t p = vub_offsets[pivot]; p < vub_offsets[pivot + 1] && indicator < 0; ++p) { - const i_t z = vub_indicators[p]; + for (i_t p = variable_upper_bound_offsets[pivot]; + p < variable_upper_bound_offsets[pivot + 1] && indicator < 0; + ++p) { + const i_t z = variable_upper_bound_indicators[p]; bool shared = true; for (size_t k = 0; k < members.size() && shared; ++k) { if (members[k] == pivot) { continue; } - const auto span_begin = vub_indicators.begin() + vub_offsets[members[k]]; - const auto span_end = vub_indicators.begin() + vub_offsets[members[k] + 1]; - shared = std::find(span_begin, span_end, z) != span_end; + const auto span_begin = + variable_upper_bound_indicators.begin() + variable_upper_bound_offsets[members[k]]; + const auto span_end = + variable_upper_bound_indicators.begin() + variable_upper_bound_offsets[members[k] + 1]; + shared = std::find(span_begin, span_end, z) != span_end; } if (shared) { indicator = z; } } @@ -265,7 +288,7 @@ template void strengthen_indicators(papilo::Problem& problem) { const indicator_strengthening_t strengthening(problem); - if (strengthening.num_vubs() == 0) { return; } + if (strengthening.num_variable_upper_bounds() == 0) { return; } const auto& constraint_matrix = problem.getConstraintMatrix(); const auto& lhs_values = constraint_matrix.getLeftHandSides(); @@ -362,7 +385,7 @@ void strengthen_indicators(papilo::Problem& problem) "variable upper bounds", num_implied, num_lifted, - strengthening.num_vubs()); + strengthening.num_variable_upper_bounds()); } #if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT diff --git a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp index 40df3f6e2e..81189c8e09 100644 --- a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp +++ b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp @@ -926,22 +926,23 @@ third_party_presolve_status_t third_party_presolve_t::apply_papilo( { raft::common::nvtx::range fun_scope("Apply Papilo presolve on host"); - if (category == problem_category_t::MIP && - (!reduction_allowlist_.has_value() || - reduction_allowlist_->count("indicatorstrengthening") > 0)) { - strengthen_indicators(papilo_problem); - } - // Capture original dimensions before papilo.apply() mutates papilo_problem // in place into its reduced form. - const i_t original_n_vars = static_cast(papilo_problem.getNCols()); - const i_t original_n_cons = static_cast(papilo_problem.getNRows()); - const i_t original_nnz = static_cast(papilo_problem.getConstraintMatrix().getNnz()); + const i_t original_n_vars = papilo_problem.getNCols(); + const i_t original_n_cons = papilo_problem.getNRows(); + const i_t original_nnz = papilo_problem.getConstraintMatrix().getNnz(); CUOPT_LOG_DEBUG("Original problem: %d constraints, %d variables, %d nonzeros", original_n_cons, original_n_vars, original_nnz); + + if (category == problem_category_t::MIP && indicator_strengthening_ && + (!reduction_allowlist_.has_value() || + reduction_allowlist_->count("indicatorstrengthening") > 0)) { + strengthen_indicators(papilo_problem); + } + CUOPT_LOG_INFO("\nRunning Papilo presolve (git hash %s)", PAPILO_GITHASH); if (category == problem_category_t::MIP) { dual_postsolve = false; } papilo::Presolve papilo_presolver; diff --git a/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp b/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp index f1426749c3..8dfb69fe3f 100644 --- a/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp +++ b/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp @@ -115,6 +115,8 @@ class third_party_presolve_t { reduction_allowlist_ = std::move(allowlist); } + void set_indicator_strengthening(bool enabled) { indicator_strengthening_ = enabled; } + // Apply the presolve on an simplex::user_problem in-place. Used in sub MIP and (in the future) // restarts. third_party_presolve_status_t apply_to_subproblem( @@ -225,6 +227,7 @@ class third_party_presolve_t { f_t original_objective_scaling_factor_{1}; std::optional> reduction_allowlist_{}; + bool indicator_strengthening_{true}; }; // Just for testing the conversion: user_problem -> Papilo problem -> user_problem. diff --git a/cpp/src/mip_heuristics/solve.cu b/cpp/src/mip_heuristics/solve.cu index 263122bd03..e6b6243d96 100644 --- a/cpp/src/mip_heuristics/solve.cu +++ b/cpp/src/mip_heuristics/solve.cu @@ -637,7 +637,8 @@ mip_solution_t solve_mip_helper( ? std::numeric_limits::infinity() : timer.remaining_time(); - presolver = std::make_unique>(); + presolver = std::make_unique>(); + presolver->set_indicator_strengthening(settings.indicator_strengthening); auto result = presolver->apply_presolve_from_op_problem( op_problem, cuopt::mathematical_optimization::problem_category_t::MIP, diff --git a/cpp/tests/internal/CMakeLists.txt b/cpp/tests/internal/CMakeLists.txt index 3cd1235717..3a5293a6d2 100644 --- a/cpp/tests/internal/CMakeLists.txt +++ b/cpp/tests/internal/CMakeLists.txt @@ -32,6 +32,7 @@ ConfigureTest(NUMOPT_INTERNAL_TEST ${CUOPT_TEST_DIR}/mip/block_bve_test.cu ${CUOPT_TEST_DIR}/mip/bhw_coeff_reduce_test.cpp ${CUOPT_TEST_DIR}/mip/gf2_presolve_test.cpp + ${CUOPT_TEST_DIR}/mip/indicator_strengthening_test.cpp ${CUOPT_TEST_DIR}/mip/arc_flow_test.cu ${CUOPT_TEST_DIR}/mip/single_lock_dual_aggregation_test.cpp ${CUOPT_TEST_DIR}/mip/termination_test.cu diff --git a/cpp/tests/mip/indicator_strengthening_test.cpp b/cpp/tests/mip/indicator_strengthening_test.cpp new file mode 100644 index 0000000000..9a90df418e --- /dev/null +++ b/cpp/tests/mip/indicator_strengthening_test.cpp @@ -0,0 +1,326 @@ +/* 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 + +namespace cuopt::mathematical_optimization::mip::test { + +namespace { + +constexpr double kInf = std::numeric_limits::infinity(); +constexpr double kTol = 1e-9; + +using entries_t = std::vector>; + +struct row_t { + entries_t entries; + double lhs; + double rhs; +}; + +// Every column has a lower bound of 0. A row side of +-kInf is left unbounded. +papilo::Problem build_problem(const std::vector& upper_bounds, + const std::vector& integral, + const std::vector& rows) +{ + const int num_cols = upper_bounds.size(); + const int num_rows = rows.size(); + + papilo::ProblemBuilder builder; + builder.setNumCols(num_cols); + builder.setNumRows(num_rows); + for (int col = 0; col < num_cols; ++col) { + builder.setObj(col, 0.0); + builder.setColLb(col, 0.0); + builder.setColUb(col, upper_bounds[col]); + builder.setColIntegral(col, integral[col]); + } + for (int row = 0; row < num_rows; ++row) { + for (const auto& [col, value] : rows[row].entries) { + builder.addEntry(row, col, value); + } + builder.setRowLhsInf(row, rows[row].lhs == -kInf); + builder.setRowRhsInf(row, rows[row].rhs == kInf); + if (rows[row].lhs != -kInf) { builder.setRowLhs(row, rows[row].lhs); } + if (rows[row].rhs != kInf) { builder.setRowRhs(row, rows[row].rhs); } + } + return builder.build(); +} + +// The stored (column, value) pairs of a row, in storage order. +entries_t row_entries(const papilo::Problem& problem, int row) +{ + const auto coefficients = problem.getConstraintMatrix().getRowCoefficients(row); + const int len = coefficients.getLength(); + entries_t entries; + entries.reserve(len); + for (int p = 0; p < len; ++p) { + entries.emplace_back(coefficients.getIndices()[p], coefficients.getValues()[p]); + } + return entries; +} + +double row_activity(const papilo::Problem& problem, + int row, + const std::vector& point) +{ + double activity = 0.0; + for (const auto& [col, value] : row_entries(problem, row)) { + activity += value * point[col]; + } + return activity; +} + +bool is_feasible(const papilo::Problem& problem, const std::vector& point) +{ + const auto& matrix = problem.getConstraintMatrix(); + const auto& flags = matrix.getRowFlags(); + for (int row = 0; row < matrix.getNRows(); ++row) { + const double activity = row_activity(problem, row, point); + if (!flags[row].test(papilo::RowFlag::kLhsInf) && + activity < matrix.getLeftHandSides()[row] - kTol) { + return false; + } + if (!flags[row].test(papilo::RowFlag::kRhsInf) && + activity > matrix.getRightHandSides()[row] + kTol) { + return false; + } + } + return true; +} + +// Every integer point in the box [0, upper_bounds]. The continuous columns in these tests only meet +// integral row sides, so the integer grid reaches every vertex of their feasible intervals. +std::vector> integer_points(const std::vector& upper_bounds) +{ + const int num_cols = upper_bounds.size(); + size_t num_points = 1; + for (int upper_bound : upper_bounds) { + num_points *= upper_bound + 1; + } + std::vector> points(num_points, std::vector(num_cols)); + for (size_t k = 0; k < num_points; ++k) { + size_t rest = k; + for (int col = 0; col < num_cols; ++col) { + points[k][col] = rest % (upper_bounds[col] + 1); + rest /= upper_bounds[col] + 1; + } + } + return points; +} + +} // namespace + +// Reduction 1 alone. Two groups, z0 owning {x0, x1} and z1 owning {x2, x3}. +// +// y0 <= x0 + x1 + x2 members span {z0, z1} => y0 <= z0 + z1 +// x2 + x3 - y1 >= 0 members span {z1} => y1 <= z1 (stored as >=) +// y2 <= x0 + x2 one indicator per member, nothing to aggregate +TEST(IndicatorStrengthening, ImpliedIndicatorRows) +{ + constexpr int z0 = 0; + constexpr int z1 = 1; + constexpr int x0 = 2; + constexpr int x1 = 3; + constexpr int x2 = 4; + constexpr int x3 = 5; + constexpr int y0 = 6; + constexpr int y1 = 7; + constexpr int y2 = 8; + + const std::vector upper_bounds(9, 1); + const std::vector integral(9, true); + const std::vector rows{ + {{{x0, 1.0}, {z0, -1.0}}, -kInf, 0.0}, + {{{x1, 1.0}, {z0, -1.0}}, -kInf, 0.0}, + {{{x2, 1.0}, {z1, -1.0}}, -kInf, 0.0}, + {{{x3, 1.0}, {z1, -1.0}}, -kInf, 0.0}, + {{{y0, 1.0}, {x0, -1.0}, {x1, -1.0}, {x2, -1.0}}, -kInf, 0.0}, + {{{x2, 1.0}, {x3, 1.0}, {y1, -1.0}}, 0.0, kInf}, + {{{y2, 1.0}, {x0, -1.0}, {x2, -1.0}}, -kInf, 0.0}, + }; + const int num_rows = rows.size(); + + const auto original = build_problem(upper_bounds, integral, rows); + auto strengthened = build_problem(upper_bounds, integral, rows); + strengthen_indicators(strengthened); + + const auto& matrix = strengthened.getConstraintMatrix(); + ASSERT_EQ(matrix.getNRows(), num_rows + 2); + for (int row = 0; row < num_rows; ++row) { + EXPECT_EQ(row_entries(strengthened, row), row_entries(original, row)) << "row " << row; + } + + const int implied_y0 = num_rows; + const int implied_y1 = num_rows + 1; + EXPECT_EQ(row_entries(strengthened, implied_y0), (entries_t{{z0, -1.0}, {z1, -1.0}, {y0, 1.0}})); + EXPECT_EQ(row_entries(strengthened, implied_y1), (entries_t{{z1, -1.0}, {y1, 1.0}})); + for (int row : {implied_y0, implied_y1}) { + EXPECT_TRUE(matrix.getRowFlags()[row].test(papilo::RowFlag::kLhsInf)) << "row " << row; + EXPECT_FALSE(matrix.getRowFlags()[row].test(papilo::RowFlag::kRhsInf)) << "row " << row; + EXPECT_EQ(matrix.getRightHandSides()[row], 0.0) << "row " << row; + } + + for (const auto& point : integer_points(upper_bounds)) { + EXPECT_EQ(is_feasible(strengthened, point), is_feasible(original, point)); + } + + // LP point: y0 = 1 is covered by x0 + x1 + x2 = 1, but z0 + z1 = 2/3. + std::vector fractional(9, 0.0); + fractional[z0] = fractional[z1] = 1.0 / 3.0; + fractional[x0] = fractional[x1] = fractional[x2] = 1.0 / 3.0; + fractional[y0] = 1.0; + EXPECT_TRUE(is_feasible(original, fractional)); + EXPECT_GT(row_activity(strengthened, implied_y0, fractional), kTol); +} + +// Reduction 2 alone. z0 owns {x0, x1, x2}, z1 owns {x3}, s is a continuous slack in [0, 2]. +// +// x0 + x1 + x2 - s <= 1 => x0 + x1 + x2 - s <= z0 +// -x0 - x1 - x2 >= -2 => -x0 - x1 - x2 + 2 z0 >= 0 +// x0 + x3 <= 1 members share no indicator, unchanged +// x0 + x1 <= 2 capacity is not binding, unchanged +TEST(IndicatorStrengthening, LiftedCapacityRows) +{ + constexpr int z0 = 0; + constexpr int z1 = 1; + constexpr int x0 = 2; + constexpr int x1 = 3; + constexpr int x2 = 4; + constexpr int x3 = 5; + constexpr int s = 6; + + const std::vector upper_bounds{1, 1, 1, 1, 1, 1, 2}; + const std::vector integral{true, true, true, true, true, true, false}; + const std::vector rows{ + {{{x0, 1.0}, {z0, -1.0}}, -kInf, 0.0}, + {{{x1, 1.0}, {z0, -1.0}}, -kInf, 0.0}, + {{{x2, 1.0}, {z0, -1.0}}, -kInf, 0.0}, + {{{x3, 1.0}, {z1, -1.0}}, -kInf, 0.0}, + {{{x0, 1.0}, {x1, 1.0}, {x2, 1.0}, {s, -1.0}}, -kInf, 1.0}, + {{{x0, -1.0}, {x1, -1.0}, {x2, -1.0}}, -2.0, kInf}, + {{{x0, 1.0}, {x3, 1.0}}, -kInf, 1.0}, + {{{x0, 1.0}, {x1, 1.0}}, -kInf, 2.0}, + }; + const int num_rows = rows.size(); + const int lifted_le = 4; + const int lifted_ge = 5; + + const auto original = build_problem(upper_bounds, integral, rows); + auto strengthened = build_problem(upper_bounds, integral, rows); + strengthen_indicators(strengthened); + + const auto& matrix = strengthened.getConstraintMatrix(); + ASSERT_EQ(matrix.getNRows(), num_rows); + for (int row = 0; row < num_rows; ++row) { + if (row == lifted_le || row == lifted_ge) { continue; } + EXPECT_EQ(row_entries(strengthened, row), row_entries(original, row)) << "row " << row; + EXPECT_EQ(matrix.getRightHandSides()[row], + original.getConstraintMatrix().getRightHandSides()[row]) + << "row " << row; + } + + EXPECT_EQ(row_entries(strengthened, lifted_le), + (entries_t{{z0, -1.0}, {x0, 1.0}, {x1, 1.0}, {x2, 1.0}, {s, -1.0}})); + EXPECT_TRUE(matrix.getRowFlags()[lifted_le].test(papilo::RowFlag::kLhsInf)); + EXPECT_EQ(matrix.getRightHandSides()[lifted_le], 0.0); + + EXPECT_EQ(row_entries(strengthened, lifted_ge), + (entries_t{{z0, 2.0}, {x0, -1.0}, {x1, -1.0}, {x2, -1.0}})); + EXPECT_TRUE(matrix.getRowFlags()[lifted_ge].test(papilo::RowFlag::kRhsInf)); + EXPECT_EQ(matrix.getLeftHandSides()[lifted_ge], 0.0); + + for (const auto& point : integer_points(upper_bounds)) { + EXPECT_EQ(is_feasible(strengthened, point), is_feasible(original, point)); + } + + // LP point: x0 + x1 + x2 = 1 fits the capacity, but z0 = 1/3 does not pay for it. + std::vector fractional(7, 0.0); + fractional[z0] = 1.0 / 3.0; + fractional[x0] = fractional[x1] = fractional[x2] = 1.0 / 3.0; + EXPECT_TRUE(is_feasible(original, fractional)); + EXPECT_GT(row_activity(strengthened, lifted_le, fractional), kTol); + EXPECT_LT(row_activity(strengthened, lifted_ge, fractional), -kTol); +} + +// Both reductions on one set-cover-like model: z0 owns {x0, x1}, z1 owns {x2, x3}, y must be +// covered by one of the members, and each group has a capacity row. +// +// y <= x0 + x1 + x2 + x3 => appended y <= z0 + z1 +// x0 + x1 - s <= 1 => x0 + x1 - s <= z0 +// x2 + x3 <= 1 => x2 + x3 <= z1 +TEST(IndicatorStrengthening, ImpliedAndLiftedRows) +{ + constexpr int z0 = 0; + constexpr int z1 = 1; + constexpr int x0 = 2; + constexpr int x1 = 3; + constexpr int x2 = 4; + constexpr int x3 = 5; + constexpr int y = 6; + constexpr int s = 7; + + const std::vector upper_bounds{1, 1, 1, 1, 1, 1, 1, 2}; + const std::vector integral{true, true, true, true, true, true, true, false}; + const std::vector rows{ + {{{x0, 1.0}, {z0, -1.0}}, -kInf, 0.0}, + {{{x1, 1.0}, {z0, -1.0}}, -kInf, 0.0}, + {{{x2, 1.0}, {z1, -1.0}}, -kInf, 0.0}, + {{{x3, 1.0}, {z1, -1.0}}, -kInf, 0.0}, + {{{y, 1.0}, {x0, -1.0}, {x1, -1.0}, {x2, -1.0}, {x3, -1.0}}, -kInf, 0.0}, + {{{x0, 1.0}, {x1, 1.0}, {s, -1.0}}, -kInf, 1.0}, + {{{x2, 1.0}, {x3, 1.0}}, -kInf, 1.0}, + }; + const int num_rows = rows.size(); + const int lifted_z0 = 5; + const int lifted_z1 = 6; + const int implied_y = num_rows; + + const auto original = build_problem(upper_bounds, integral, rows); + auto strengthened = build_problem(upper_bounds, integral, rows); + strengthen_indicators(strengthened); + + const auto& matrix = strengthened.getConstraintMatrix(); + ASSERT_EQ(matrix.getNRows(), num_rows + 1); + for (int row = 0; row < lifted_z0; ++row) { + EXPECT_EQ(row_entries(strengthened, row), row_entries(original, row)) << "row " << row; + } + + EXPECT_EQ(row_entries(strengthened, lifted_z0), + (entries_t{{z0, -1.0}, {x0, 1.0}, {x1, 1.0}, {s, -1.0}})); + EXPECT_EQ(row_entries(strengthened, lifted_z1), (entries_t{{z1, -1.0}, {x2, 1.0}, {x3, 1.0}})); + EXPECT_EQ(row_entries(strengthened, implied_y), (entries_t{{z0, -1.0}, {z1, -1.0}, {y, 1.0}})); + for (int row : {lifted_z0, lifted_z1, implied_y}) { + EXPECT_TRUE(matrix.getRowFlags()[row].test(papilo::RowFlag::kLhsInf)) << "row " << row; + EXPECT_FALSE(matrix.getRowFlags()[row].test(papilo::RowFlag::kRhsInf)) << "row " << row; + EXPECT_EQ(matrix.getRightHandSides()[row], 0.0) << "row " << row; + } + + for (const auto& point : integer_points(upper_bounds)) { + EXPECT_EQ(is_feasible(strengthened, point), is_feasible(original, point)); + } + + // LP point: every original row holds, and each of the three new rows is violated by 1/3. + std::vector fractional(8, 0.0); + fractional[z0] = fractional[z1] = 1.0 / 3.0; + fractional[x0] = fractional[x1] = fractional[x2] = fractional[x3] = 1.0 / 3.0; + fractional[y] = 1.0; + EXPECT_TRUE(is_feasible(original, fractional)); + for (int row : {lifted_z0, lifted_z1, implied_y}) { + EXPECT_GT(row_activity(strengthened, row, fractional), kTol) << "row " << row; + } +} + +} // namespace cuopt::mathematical_optimization::mip::test From c170fddadf1320eb5df6402cc34f0a5d88d4054f Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Thu, 1 Oct 2026 18:11:53 +0200 Subject: [PATCH 11/12] add more comments. tightened the formulation for implied indicator. Signed-off-by: Nicolas L. Guidotti --- .../presolve/indicator_strengthening.cpp | 67 +++++++++++++------ 1 file changed, 46 insertions(+), 21 deletions(-) diff --git a/cpp/src/mip_heuristics/presolve/indicator_strengthening.cpp b/cpp/src/mip_heuristics/presolve/indicator_strengthening.cpp index 6ef2c31e86..b78cd3fa8a 100644 --- a/cpp/src/mip_heuristics/presolve/indicator_strengthening.cpp +++ b/cpp/src/mip_heuristics/presolve/indicator_strengthening.cpp @@ -30,13 +30,8 @@ template class indicator_strengthening_t { public: indicator_strengthening_t(const papilo::Problem& problem); - - // Implication row y - sum_{j in S} x_j <= 0: one +1 head against a tail of -1 members. i_t add_implied_indicator_rows(papilo::Vec>& implied_entries) const; - - // Capacity row sum_{i in S} x_i - s <= K with every x_i bounded by a common indicator z. - i_t lift_capacity_rows(papilo::Vec>& lifted_entries) const; - + i_t find_lift_capacity_rows(papilo::Vec>& lifted_entries) const; i_t num_variable_upper_bounds() const { return variable_upper_bound_indicators.size(); } i_t row_orientation(i_t row) const { return orientation[row]; } @@ -77,6 +72,7 @@ indicator_strengthening_t::indicator_strengthening_t(const papilo::Pro const i_t num_rows = constraint_matrix.getNRows(); const i_t num_cols = problem.getNCols(); + // A binary is an integral column with finite bounds [0, 1]. is_binary.resize(num_cols); for (i_t col = 0; col < num_cols; ++col) { is_binary[col] = col_flags[col].test(papilo::ColFlag::kIntegral) && @@ -84,6 +80,8 @@ indicator_strengthening_t::indicator_strengthening_t(const papilo::Pro lower_bounds[col] == 0.0 && upper_bounds[col] == 1.0; } + // Record the orientation of every one-sided row, and collect each two-variable row x - z <= 0 + // between binaries as a variable upper bound (member x, indicator z). orientation.assign(num_rows, 0); std::vector> variable_upper_bounds; variable_upper_bounds.reserve(std::count(row_sizes.begin(), row_sizes.end(), 2)); @@ -111,6 +109,8 @@ indicator_strengthening_t::indicator_strengthening_t(const papilo::Pro } } + // Build the CSR from the (member, indicator) pairs: count the indicators of each member, + // prefix-sum the counts into offsets, then scatter the indicators into place. variable_upper_bound_offsets.assign(num_cols + 1, 0); for (const auto& variable_upper_bound : variable_upper_bounds) { ++variable_upper_bound_offsets[variable_upper_bound.first + 1]; @@ -144,6 +144,11 @@ i_t indicator_strengthening_t::add_implied_indicator_rows( indicators.reserve(max_row_size); std::vector mark(num_cols, -1); + // Loop over the constraints and identify the disjunctions (y - \sum_{j \in \mathcal{S}} x_j <= + // 0), + // $|\mathcal{S}| >= 2$ then construct a new constraint (y - \sum_{i \in \mathcal{D}} z_i <= 0), + // where + // $\mathcal{D}$ contains the first indicator associated with a member variable x_j. for (i_t row = 0; row < num_rows; ++row) { const i_t direction = orientation[row]; if (direction == 0 || row_sizes[row] < 3) { continue; } @@ -155,6 +160,9 @@ i_t indicator_strengthening_t::add_implied_indicator_rows( const i_t* indices = row_coefficients.getIndices(); const f_t* values = row_coefficients.getValues(); + // Every entry must be binary: the single +1 is the head y, and each -1 member is replaced by + // its indicator, deduplicated with `mark`. The row is rejected if no two members share an + // indicator, since the new row would then be implied by the LP relaxation. const size_t num_members = len - 1; i_t head = -1; bool usable = true; @@ -170,10 +178,10 @@ i_t indicator_strengthening_t::add_implied_indicator_rows( } else if (v == -1.0) { const i_t start = variable_upper_bound_offsets[col]; const i_t end = variable_upper_bound_offsets[col + 1]; - const i_t count = std::max(end - start, 1); - for (i_t k = 0; k < count && usable; ++k) { - const i_t z = start < end ? variable_upper_bound_indicators[start + k] : col; - if (mark[z] == row) { continue; } + + // Take the first indicator that owns this variable, or the member itself if it has none. + const i_t z = start < end ? variable_upper_bound_indicators[start] : col; + if (mark[z] != row) { mark[z] = row; indicators.push_back(z); usable = indicators.size() < num_members; @@ -182,6 +190,7 @@ i_t indicator_strengthening_t::add_implied_indicator_rows( usable = false; } } + // Skip rows without a head, or whose head is also one of the indicators (a trivial row). if (!usable || head < 0 || mark[head] == row) { continue; } // In the previous loop we insert the indicators in a random order, however, Papilo @@ -203,7 +212,7 @@ i_t indicator_strengthening_t::add_implied_indicator_rows( } template -i_t indicator_strengthening_t::lift_capacity_rows( +i_t indicator_strengthening_t::find_lift_capacity_rows( papilo::Vec>& lifted_entries) const { const auto& constraint_matrix = problem.getConstraintMatrix(); @@ -220,7 +229,11 @@ i_t indicator_strengthening_t::lift_capacity_rows( std::vector members; members.reserve(max_row_size); + // Loop over the constraints and identify the capacity row \sum_{i \in \mathcal{S}} x_i - s <= K, + // and then tighten it as \sum_{i \in \mathcal{S}} x_i - s <= Kz, for a indicator z with + // x_i <= z for (i_t row = 0; row < num_rows; ++row) { + // A capacity row is a one-sided row with a positive right-hand side K. const i_t direction = orientation[row]; if (direction == 0) { continue; } const f_t capacity = direction * (direction == 1 ? rhs_values[row] : lhs_values[row]); @@ -231,6 +244,8 @@ i_t indicator_strengthening_t::lift_capacity_rows( const i_t* indices = row_coefficients.getIndices(); const f_t* values = row_coefficients.getValues(); + // Every integral entry must be a +1 binary member with at least one indicator, and every + // continuous entry a -s slack with s >= 0. The pivot is the member with the fewest indicators. members.clear(); i_t pivot = -1; bool usable = true; @@ -253,28 +268,32 @@ i_t indicator_strengthening_t::lift_capacity_rows( v < 0.0 && !col_flags[col].test(papilo::ColFlag::kLbInf) && lower_bounds[col] >= 0.0; } } + // Lifting needs at least two members, and only tightens the row when K < |S|. if (!usable || members.size() < 2) { continue; } const f_t n_members = members.size(); if (capacity >= n_members) { continue; } - i_t indicator = -1; - for (i_t p = variable_upper_bound_offsets[pivot]; - p < variable_upper_bound_offsets[pivot + 1] && indicator < 0; - ++p) { + // A shared indicator must be in the pivot's list, so only its candidates are tried; each is + // searched for in the lists of the other members, and the first one found in all of them wins. + i_t indicator = -1; + i_t pivot_start = variable_upper_bound_offsets[pivot]; + i_t pivot_end = variable_upper_bound_offsets[pivot + 1]; + + for (i_t p = pivot_start; p < pivot_end && indicator < 0; ++p) { const i_t z = variable_upper_bound_indicators[p]; bool shared = true; for (size_t k = 0; k < members.size() && shared; ++k) { if (members[k] == pivot) { continue; } - const auto span_begin = - variable_upper_bound_indicators.begin() + variable_upper_bound_offsets[members[k]]; - const auto span_end = - variable_upper_bound_indicators.begin() + variable_upper_bound_offsets[members[k] + 1]; - shared = std::find(span_begin, span_end, z) != span_end; + i_t start = variable_upper_bound_offsets[members[k]]; + i_t end = variable_upper_bound_offsets[members[k] + 1]; + std::span indicator_span(variable_upper_bound_indicators.begin() + start, end - start); + shared = std::find(indicator_span.begin(), indicator_span.end(), z) != indicator_span.end(); } if (shared) { indicator = z; } } if (indicator < 0) { continue; } + // Add the -K z entry to the row; `strengthen_indicators` then sets its right-hand side to 0. lifted_entries.emplace_back(row, indicator, direction * -capacity); ++num_lifted; } @@ -299,6 +318,7 @@ void strengthen_indicators(papilo::Problem& problem) const i_t num_rows = constraint_matrix.getNRows(); const i_t num_cols = problem.getNCols(); + // Bound the number of implied entries and lifted rows from above to reserve the buffers. i_t max_implied_entries = 0; i_t max_lifted_rows = 0; for (i_t row = 0; row < num_rows; ++row) { @@ -317,9 +337,11 @@ void strengthen_indicators(papilo::Problem& problem) implied_entries.reserve(max_implied_entries); lifted_entries.reserve(max_lifted_rows); const i_t num_implied = strengthening.add_implied_indicator_rows(implied_entries); - const i_t num_lifted = strengthening.lift_capacity_rows(lifted_entries); + const i_t num_lifted = strengthening.find_lift_capacity_rows(lifted_entries); if (num_implied == 0 && num_lifted == 0) { return; } + // Rebuild the matrix as row-sorted triplets, inserting each lifted entry at its column position + // within its row, then append the implied rows after the original ones. papilo::Vec> entries; entries.reserve(constraint_matrix.getNnz() + lifted_entries.size() + implied_entries.size()); auto lifted = lifted_entries.begin(); @@ -342,6 +364,7 @@ void strengthen_indicators(papilo::Problem& problem) } entries.insert(entries.end(), implied_entries.begin(), implied_entries.end()); + // Lifted rows get a zero right-hand side, and the implied rows read y - \sum_{z} z <= 0. papilo::Vec lhs; papilo::Vec rhs; papilo::Vec flags; @@ -365,6 +388,7 @@ void strengthen_indicators(papilo::Problem& problem) rhs.resize(num_rows + num_implied, 0.0); flags.resize(num_rows + num_implied, implied_flags); + // Name the implied rows if the problem carries constraint names. const auto& constraint_names = problem.getConstraintNames(); if (!constraint_names.empty()) { papilo::Vec names; @@ -376,6 +400,7 @@ void strengthen_indicators(papilo::Problem& problem) problem.setConstraintNames(std::move(names)); } + // Replace the constraint matrix with the strengthened one. papilo::SparseStorage storage( std::move(entries), num_rows + num_implied, num_cols, true, 4.0, 30); problem.setConstraintMatrix(std::move(storage), std::move(lhs), std::move(rhs), std::move(flags)); From 5cfad1d4fb27cfc2406911f20b1e79f258fcd52b Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Fri, 2 Oct 2026 10:10:53 +0200 Subject: [PATCH 12/12] add missing include Signed-off-by: Nicolas L. Guidotti --- cpp/src/mip_heuristics/presolve/indicator_strengthening.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/cpp/src/mip_heuristics/presolve/indicator_strengthening.cpp b/cpp/src/mip_heuristics/presolve/indicator_strengthening.cpp index b78cd3fa8a..d3fa383ec4 100644 --- a/cpp/src/mip_heuristics/presolve/indicator_strengthening.cpp +++ b/cpp/src/mip_heuristics/presolve/indicator_strengthening.cpp @@ -11,6 +11,7 @@ #include #include +#include #include #include #include