diff --git a/cpp/src/mip_heuristics/CMakeLists.txt b/cpp/src/mip_heuristics/CMakeLists.txt index 538e3c49a..35a4f122e 100644 --- a/cpp/src/mip_heuristics/CMakeLists.txt +++ b/cpp/src/mip_heuristics/CMakeLists.txt @@ -12,6 +12,7 @@ set(MIP_LP_NECESSARY_FILES ${CMAKE_CURRENT_SOURCE_DIR}/solver_solution.cu ${CMAKE_CURRENT_SOURCE_DIR}/local_search/rounding/simple_rounding.cu ${CMAKE_CURRENT_SOURCE_DIR}/presolve/third_party_presolve.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/presolve/single_lock_dual_aggregation.cpp ${CMAKE_CURRENT_SOURCE_DIR}/presolve/gf2_presolve.cpp ${CMAKE_CURRENT_SOURCE_DIR}/solution/solution.cu ) diff --git a/cpp/src/mip_heuristics/presolve/single_lock_dual_aggregation.cpp b/cpp/src/mip_heuristics/presolve/single_lock_dual_aggregation.cpp new file mode 100644 index 000000000..76f7b98d2 --- /dev/null +++ b/cpp/src/mip_heuristics/presolve/single_lock_dual_aggregation.cpp @@ -0,0 +1,359 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#include "single_lock_dual_aggregation.hpp" + +#include +#include + +#include +#include + +namespace cuopt::linear_programming::detail { + +// Single-Lock Dual Aggregation +// +// For a binary variable x with exactly one "up-lock" (one constraint preventing +// it from increasing), we try to prove an implication y=0 => x=0 via activity +// bounds on the locking row. If additionally the row is non-binding when y=1 +// (no capacity competition), we can substitute x = y, eliminating a variable. +// +// Symmetric logic applies for "down-lock" candidates (one constraint preventing +// decrease), proving y=1 => x=1. + +template +papilo::PresolveStatus SingleLockDualAggregation::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& objective = problem.getObjective().coefficients; + + const int nrows = constraint_matrix.getNRows(); + const int ncols = problem.getNCols(); + const double tlim = problemUpdate.getPresolveOptions().tlim; + const f_t tol = num.getFeasTol(); + + // ========================================================================= + // Step 1: Lock Counting — O(nnz) + // + // An "up-lock" on column j means a constraint prevents j from increasing: + // - a_j > 0 in a <= row, or a_j < 0 in a >= row. + // "Down-lock" is the reverse. Equality rows lock both directions. + // We record the row index of the first lock; a second lock invalidates it. + // ========================================================================= + + enum lock_dir { UP = 0, DOWN = 1 }; + enum bound_side { LOWER = 0, UPPER = 1 }; + std::vector locks[2] = {std::vector(ncols, 0), std::vector(ncols, 0)}; + std::vector lock_row[2] = {std::vector(ncols, -1), std::vector(ncols, -1)}; + + for (int row = 0; row < nrows; ++row) { + if (this->is_time_exceeded(timer, tlim)) return papilo::PresolveStatus::kUnchanged; + if (row_flags[row].test(papilo::RowFlag::kRedundant)) continue; + + // Row direction: 'E' = equality, 'L' = <=, 'G' = >=, 'R' = ranged/free (skip) + bool lhs_inf = row_flags[row].test(papilo::RowFlag::kLhsInf); + bool rhs_inf = row_flags[row].test(papilo::RowFlag::kRhsInf); + char row_dir = (!lhs_inf && !rhs_inf) ? 'E' : (lhs_inf ? 'L' : (rhs_inf ? 'G' : 'R')); + if (row_dir == 'R') continue; + + auto row_coeff = constraint_matrix.getRowCoefficients(row); + const int* cols = row_coeff.getIndices(); + const f_t* coefs = row_coeff.getValues(); + const int length = row_coeff.getLength(); + + // Record the index of the locking row. + // If more than one lock exists, mark the col as excluded from the search. + auto record_lock = [&](lock_dir dir, int col) { + if (locks[dir][col]++ == 0) + lock_row[dir][col] = row; + else + lock_row[dir][col] = -1; + }; + + if (row_dir == 'E') { + // Equality: locks both directions + for (int j = 0; j < length; ++j) { + record_lock(UP, cols[j]); + record_lock(DOWN, cols[j]); + } + } else { + // One-sided: directions swap between L (<=) and G (>=) + lock_dir pos_dir = (row_dir == 'L') ? UP : DOWN; + lock_dir neg_dir = (row_dir == 'L') ? DOWN : UP; + for (int j = 0; j < length; ++j) { + if (coefs[j] > 0) + record_lock(pos_dir, cols[j]); + else if (coefs[j] < 0) + record_lock(neg_dir, cols[j]); + } + } + } + + // ========================================================================= + // Step 2: Candidate Identification — O(ncols) + // + // Upward candidates: binary, single up-lock, c <= 0 (objective doesn't + // penalize increase — needed so x pushes against the lock or is indifferent). + // Downward: symmetric with single down-lock, c >= 0. + // ========================================================================= + + struct candidate_t { + int col; + int locking_row; + lock_dir dir; + }; + std::vector candidates; + candidates.reserve(std::min(ncols, nrows)); + + for (int col = 0; col < ncols; ++col) { + if (col_flags[col].test(papilo::ColFlag::kFixed, papilo::ColFlag::kSubstituted)) continue; + if (!is_binary_or_implied(col, col_flags.data(), lower_bounds.data(), upper_bounds.data())) + continue; + // Skip singletons: PaPILO's stuffing presolver handles these. + if (constraint_matrix.getColumnCoefficients(col).getLength() <= 1) continue; + + // can be turned into strict checks if we need to guarantee + // that we never cut off any optimal solution + if (locks[UP][col] == 1 && objective[col] <= 0) + candidates.push_back({col, lock_row[UP][col], UP}); + else if (locks[DOWN][col] == 1 && objective[col] >= 0) + candidates.push_back({col, lock_row[DOWN][col], DOWN}); + } + + if (this->is_time_exceeded(timer, tlim) || candidates.empty()) + return papilo::PresolveStatus::kUnchanged; + + // ========================================================================= + // Step 3: Mini-Probing — O(L + K) per row + // + // For each locking row (L nonzeros, K candidates), we prove implications by + // fixing two variables and checking if the row's activity bounds are violated: + // - Fix candidate x to its "bad" bound (ub for upward, lb for downward) + // - Fix master y to its "unfavorable" bound (0 for upward, 1 for downward) + // - If the resulting minimum (LEQ) or maximum (GEQ) activity exceeds the + // row's bound, the combination is infeasible, proving y_unfav => x_safe. + // + // The master y is the binary variable in the row whose coefficient best + // amplifies the violation. We track the top-2 most extreme coefficients + // (neg_y for most negative, pos_y for most positive) so that if the + // candidate itself is the top-1 extremum, we can fall back to top-2. + // This keeps master selection O(1) per candidate instead of O(L). + // + // Candidates are sorted by lock_row so all K candidates sharing a row are + // processed together in a single O(L) scan, yielding O(L+K) per row group. + // + // dense_row_coefs[] is an ncols-sized scratch array giving O(1) coefficient + // lookup by column index; populated and cleaned per row in O(L). + // ========================================================================= + + std::sort(candidates.begin(), candidates.end(), [](const candidate_t& a, const candidate_t& b) { + return a.locking_row < b.locking_row; + }); + + struct top2_t { + std::pair top1{-1, 0}, top2{-1, 0}; + + void update(int idx, f_t val, bound_side side) + { + auto better = [side](f_t a, f_t b) { return side == LOWER ? a < b : a > b; }; + if (top1.first == -1 || better(val, top1.second)) { + top2 = top1; + top1 = {idx, val}; + } else if (top2.first == -1 || better(val, top2.second)) { + top2 = {idx, val}; + } + } + }; + + int n_substitutions = 0; + std::vector dense_row_coefs(ncols, f_t{0}); + std::vector substituted(ncols, false); + + auto cand_it = candidates.begin(); + while (cand_it != candidates.end()) { + if (this->is_time_exceeded(timer, tlim)) break; + + int r = cand_it->locking_row; + if (r < 0) { + ++cand_it; + continue; + } + + // advance row_end to the first candidate with a different locking_row + auto row_end = std::find_if( + cand_it, candidates.end(), [r](const candidate_t& c) { return c.locking_row != r; }); + + auto row_coeff = constraint_matrix.getRowCoefficients(r); + const int* cols = row_coeff.getIndices(); + const f_t* coefs = row_coeff.getValues(); + const int length = row_coeff.getLength(); + + bool has_lhs = !row_flags[r].test(papilo::RowFlag::kLhsInf); + bool has_rhs = !row_flags[r].test(papilo::RowFlag::kRhsInf); + + // A_min / A_max: tightest possible activity of the row over all variable bounds + f_t A_min = 0, A_max = 0; + bool can_reach_neg_inf = false, can_reach_pos_inf = false; + top2_t neg_y, pos_y; + + for (int j = 0; j < length; ++j) { + int col = cols[j]; + f_t coef = coefs[j]; + bool lb_inf = col_flags[col].test(papilo::ColFlag::kLbInf); + bool ub_inf = col_flags[col].test(papilo::ColFlag::kUbInf); + + dense_row_coefs[col] = coef; + + // coef > 0: min activity uses lb, max uses ub; coef < 0: swapped + bool min_inf = (coef > 0) ? lb_inf : ub_inf; + bool max_inf = (coef > 0) ? ub_inf : lb_inf; + f_t min_bound = (coef > 0) ? lower_bounds[col] : upper_bounds[col]; + f_t max_bound = (coef > 0) ? upper_bounds[col] : lower_bounds[col]; + + if (min_inf) + can_reach_neg_inf = true; + else + A_min += coef * min_bound; + if (max_inf) + can_reach_pos_inf = true; + else + A_max += coef * max_bound; + + if (col_flags[col].test(papilo::ColFlag::kFixed, papilo::ColFlag::kSubstituted)) continue; + if (!is_binary_or_implied(col, col_flags.data(), lower_bounds.data(), upper_bounds.data())) + continue; + if (lower_bounds[col] == upper_bounds[col]) continue; + + neg_y.update(col, coef, LOWER); + pos_y.update(col, coef, UPPER); + } + + // LEQ probe needs finite A_min; GEQ probe needs finite A_max + bool use_leq_check = has_rhs && !can_reach_neg_inf; + bool use_geq_check = has_lhs && !can_reach_pos_inf; + + // Probe: replace cand and y's min/max contributions with their fixed test + // values, then check if the resulting activity violates the row bound. + // Both candidate and master are binary [0,1], so min/max contributions simplify + auto evaluate = [&](f_t cand_coeff, bool is_upward, int y_col, f_t y_coef) -> bool { + if (y_col < 0) return false; + f_t cand_test = is_upward ? cand_coeff : f_t{0}; + f_t y_test = is_upward ? f_t{0} : y_coef; + f_t test = cand_test + y_test; + + if (use_leq_check) { + f_t probed_min = A_min - std::min(f_t{0}, cand_coeff) - std::min(f_t{0}, y_coef) + test; + if (probed_min > rhs_values[r] + tol) return true; + } + if (use_geq_check) { + f_t probed_max = A_max - std::max(f_t{0}, cand_coeff) - std::max(f_t{0}, y_coef) + test; + if (probed_max < lhs_values[r] - tol) return true; + } + return false; + }; + + // Return the best master from the top-2 tracker, skipping excluded columns. + auto pick_master = [&substituted](const top2_t& t, int exclude) -> std::pair { + if (t.top1.first >= 0 && t.top1.first != exclude && !substituted[t.top1.first]) return t.top1; + if (t.top2.first >= 0 && t.top2.first != exclude && !substituted[t.top2.first]) return t.top2; + return {-1, f_t{0}}; + }; + + for (auto ci = cand_it; ci != row_end; ++ci) { + auto [cand, locking_row, dir] = *ci; + if (substituted[cand]) continue; + + bool is_upward = (dir == UP); + f_t cand_coeff = dense_row_coefs[cand]; + + bool proven = false; + int master_col = -1; + + // For LEQ upward: y=0 zeroes out y's contribution, so the best master + // is the one with the most negative coefficient (maximizes probed_min). + // For LEQ downward: y=1 adds y's coefficient, so pick the most positive. + auto try_prove = [&](bool check, const top2_t& tracker) { + if (!check || proven) return; + auto [y, yc] = pick_master(tracker, cand); + if (evaluate(cand_coeff, is_upward, y, yc)) { + proven = true; + master_col = y; + } + }; + try_prove(use_leq_check, is_upward ? neg_y : pos_y); + try_prove(use_geq_check, is_upward ? pos_y : neg_y); + if (!proven) continue; + + // The probe proves a one-directional implication (e.g. y=0 => x=0). + // The substitution x=y also asserts the reverse (y=1 => x=1), which is + // only safe if forcing x to its bound doesn't starve other variables of + // capacity in the locking row. Verify the row becomes globally non-binding + // when y is in its favorable state. + f_t y_coef_val = dense_row_coefs[master_col]; + f_t fav_y_contrib = is_upward ? y_coef_val : 0.0; + + auto check_side = + [&](bool active, bool unbounded, f_t activity, f_t orig_y, f_t bound, bound_side side) { + if (!active || !proven) return; + if (unbounded) { + proven = false; + return; + } + f_t fav = activity - orig_y + fav_y_contrib; + if (side == UPPER ? fav > bound + tol : fav < bound - tol) proven = false; + }; + check_side( + has_rhs, can_reach_pos_inf, A_max, std::max(0.0, y_coef_val), rhs_values[r], UPPER); + check_side( + has_lhs, can_reach_neg_inf, A_min, std::min(0.0, y_coef_val), lhs_values[r], LOWER); + if (!proven) continue; + + substituted[cand] = true; + reductions.replaceCol(cand, master_col, f_t{1}, f_t{0}); + ++n_substitutions; + } + + for (int j = 0; j < length; ++j) + dense_row_coefs[cols[j]] = 0; + cand_it = row_end; + } + + if (n_substitutions == 0) return papilo::PresolveStatus::kUnchanged; + + CUOPT_LOG_DEBUG("Single-lock dual aggregation: %d candidates, %d substitutions", + (int)candidates.size(), + n_substitutions); + + return papilo::PresolveStatus::kReduced; +} + +#define INSTANTIATE(F_TYPE) template class SingleLockDualAggregation; + +#if MIP_INSTANTIATE_FLOAT +INSTANTIATE(float) +#endif + +#if MIP_INSTANTIATE_DOUBLE +INSTANTIATE(double) +#endif + +#undef INSTANTIATE + +} // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip_heuristics/presolve/single_lock_dual_aggregation.hpp b/cpp/src/mip_heuristics/presolve/single_lock_dual_aggregation.hpp new file mode 100644 index 000000000..dcf3b0923 --- /dev/null +++ b/cpp/src/mip_heuristics/presolve/single_lock_dual_aggregation.hpp @@ -0,0 +1,56 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2025-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 +#endif +#include +#include +#include +#include +#if !defined(__clang__) +#pragma GCC diagnostic pop +#endif + +namespace cuopt::linear_programming::detail { + +template +class SingleLockDualAggregation : public papilo::PresolveMethod { + public: + SingleLockDualAggregation() : papilo::PresolveMethod() + { + this->setName("single_lock_dual_aggregation"); + this->setType(papilo::PresolverType::kIntegralCols); + this->setTiming(papilo::PresolverTiming::kMedium); + } + + 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; + + private: + bool is_binary_or_implied(int col, + const papilo::Flags* col_flags, + const f_t* lower_bounds, + const f_t* upper_bounds) const + { + if (!col_flags[col].test(papilo::ColFlag::kIntegral) && + !col_flags[col].test(papilo::ColFlag::kImplInt)) + return false; + if (col_flags[col].test(papilo::ColFlag::kLbInf)) return false; + if (col_flags[col].test(papilo::ColFlag::kUbInf)) return false; + return lower_bounds[col] == 0.0 && upper_bounds[col] == 1.0; + } +}; + +} // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp index 5a89393a6..ca731f33f 100644 --- a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp +++ b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp @@ -21,6 +21,7 @@ #endif #include #include +#include #include #include #include @@ -496,8 +497,10 @@ void set_presolve_methods(papilo::Presolve& presolver, using uptr = std::unique_ptr>; if (category == problem_category_t::MIP) { - // cuOpt custom GF2 presolver + // cuOpt custom presolvers presolver.addPresolveMethod(uptr(new cuopt::linear_programming::detail::GF2Presolve())); + presolver.addPresolveMethod( + uptr(new cuopt::linear_programming::detail::SingleLockDualAggregation())); } // fast presolvers presolver.addPresolveMethod(uptr(new papilo::SingletonCols())); @@ -636,7 +639,6 @@ std::optional> third_party_presolve_t( papilo_presolver, category, op_problem.get_n_constraints(), op_problem.get_n_variables()); - // Disable papilo logs papilo_presolver.setVerbosityLevel(papilo::VerbosityLevel::kQuiet); auto result = papilo_presolver.apply(papilo_problem); diff --git a/cpp/tests/mip/termination_test.cu b/cpp/tests/mip/termination_test.cu index 5f21d294f..b0bf32e34 100644 --- a/cpp/tests/mip/termination_test.cu +++ b/cpp/tests/mip/termination_test.cu @@ -118,6 +118,13 @@ TEST(termination_status, gf2_presolve_infeasible) EXPECT_EQ(termination_status, mip_termination_status_t::Infeasible); } +TEST(termination_status, slda_presolve_optimal) +{ + auto [termination_status, obj_val, lb] = test_mps_file("mip/neos-787933.mps", 30, false); + EXPECT_EQ(termination_status, mip_termination_status_t::Optimal); + EXPECT_NEAR(obj_val, 30.0, 1e-6); +} + TEST(termination_status, bb_infeasible_test) { // First, check that presolve doesn't reduce the problem to infeasibility diff --git a/datasets/mip/download_miplib_test_dataset.sh b/datasets/mip/download_miplib_test_dataset.sh index dc2dd7966..45a90ac01 100755 --- a/datasets/mip/download_miplib_test_dataset.sh +++ b/datasets/mip/download_miplib_test_dataset.sh @@ -23,6 +23,7 @@ INSTANCES=( "swath1" "enlight_hard" "enlight11" + "neos-787933" "supportcase22" )