diff --git a/cpp/src/dual_simplex/bounds_strengthening.cpp b/cpp/src/dual_simplex/bounds_strengthening.cpp index 6506ce4873..e137343eee 100644 --- a/cpp/src/dual_simplex/bounds_strengthening.cpp +++ b/cpp/src/dual_simplex/bounds_strengthening.cpp @@ -211,9 +211,14 @@ bool bounds_strengthening_t::bounds_strengthening( bool lb_updated = std::abs(new_lb - old_lb) > 1e3 * settings.primal_tol; bool ub_updated = std::abs(new_ub - old_ub) > 1e3 * settings.primal_tol; - - new_lb = std::max(new_lb, lower_bounds[k]); - new_ub = std::min(new_ub, upper_bounds[k]); + if (lb_updated) + new_lb = std::max(new_lb, lower_bounds[k]); + else + new_lb = old_lb; + if (ub_updated) + new_ub = std::min(new_ub, upper_bounds[k]); + else + new_ub = old_ub; if (new_lb > new_ub + settings.primal_tol) { settings.log.debug( diff --git a/cpp/src/mip_heuristics/CMakeLists.txt b/cpp/src/mip_heuristics/CMakeLists.txt index a200d4265b..22816acbd4 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 0000000000..5c950a2af0 --- /dev/null +++ b/cpp/src/mip_heuristics/presolve/single_lock_dual_aggregation.cpp @@ -0,0 +1,483 @@ +/* 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. + +namespace { + +enum lock_dir { UP = 0, DOWN = 1 }; +enum bound_side { LOWER = 0, UPPER = 1 }; + +struct candidate_t { + int col; + int locking_row; + lock_dir dir; +}; + +template +bool is_binary_or_implied(int col, + const papilo::Flags* col_flags, + const f_t* lower_bounds, + const f_t* upper_bounds) +{ + 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; +} + +template +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}; + } + } +}; + +// ========================================================================= +// 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. +// ========================================================================= + +bool time_exceeded(const papilo::Timer& timer, double tlim) +{ + return tlim != std::numeric_limits::max() && timer.getTime() >= tlim; +} + +template +void compute_single_locks(const papilo::Problem& problem, + const papilo::Timer& timer, + double tlim, + std::vector locks[2], + std::vector lock_row[2]) +{ + const auto& constraint_matrix = problem.getConstraintMatrix(); + const auto& row_flags = constraint_matrix.getRowFlags(); + const int nrows = constraint_matrix.getNRows(); + const int ncols = problem.getNCols(); + + locks[UP].assign(ncols, 0); + locks[DOWN].assign(ncols, 0); + lock_row[UP].assign(ncols, -1); + lock_row[DOWN].assign(ncols, -1); + + for (int row = 0; row < nrows; ++row) { + if (time_exceeded(timer, tlim)) return; + 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) ? 'R' : (!lhs_inf && !rhs_inf) ? 'E' : lhs_inf ? 'L' : 'G'; + 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. +// ========================================================================= + +template +std::vector collect_candidates(const papilo::Problem& problem, + const std::vector locks[2], + const std::vector lock_row[2]) +{ + const auto& constraint_matrix = problem.getConstraintMatrix(); + 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 ncols = problem.getNCols(); + const int nrows = constraint_matrix.getNRows(); + + 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}); + } + return candidates; +} + +// ========================================================================= +// Step 3: Mini-Probing +// +// 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. +// +// Candidates are sorted by lock_row so all K candidates sharing a row are +// processed together in a single O(L) scan. +// ========================================================================= + +// can't use "check_if_substitution_generates_huge_or_small_coefficients" directly since it expects +// equality rows only +template +bool substitution_numerically_stable(const papilo::ConstraintMatrix& constraint_matrix, + int cand_col) +{ + auto col_coeff = constraint_matrix.getColumnCoefficients(cand_col); + const int* rows = col_coeff.getIndices(); + const f_t* col_vals = col_coeff.getValues(); + const int col_length = col_coeff.getLength(); + + for (int k = 0; k < col_length; ++k) { + int r = rows[k]; + f_t abs_cand = std::abs(col_vals[k]); + auto row_coeff = constraint_matrix.getRowCoefficients(r); + const f_t* rvals = row_coeff.getValues(); + const int rlength = row_coeff.getLength(); + + f_t row_max = 0; + for (int p = 0; p < rlength; ++p) + row_max = std::max(row_max, std::abs(rvals[p])); + + if (abs_cand > 1e6 * row_max || abs_cand * 1e6 < row_max) return false; + } + return true; +} + +template +int try_substitutions_for_row(const papilo::Problem& problem, + const papilo::Num& num, + papilo::Reductions& reductions, + typename std::vector::iterator cand_begin, + typename std::vector::iterator cand_end, + int row, + std::vector& dense_row_coefs, + std::vector& substituted) +{ + 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; + + 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(); + + bool has_lhs = !row_flags[row].test(papilo::RowFlag::kLhsInf); + bool has_rhs = !row_flags[row].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. + // anti=false: direct (y_unfav=0 for upward, y_unfav=1 for downward) + // anti=true: complement (y_unfav=1 for upward, y_unfav=0 for downward) + auto evaluate = [&](f_t cand_coeff, bool is_upward, int y_col, f_t y_coef, bool anti) -> bool { + if (y_col < 0) return false; + f_t cand_test = is_upward ? cand_coeff : f_t{0}; + f_t y_test = (is_upward == anti) ? y_coef : f_t{0}; + 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 (num.isFeasGT(probed_min, rhs_values[row])) 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 (num.isFeasLT(probed_max, lhs_values[row])) 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}}; + }; + + int n_substitutions = 0; + + for (auto ci = cand_begin; ci != cand_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; + bool is_anti = false; + + // For LEQ upward direct: y=0 zeroes out y's contribution, so the best + // master is the one with the most negative coefficient (maximizes + // probed_min). For anti (complement): y=1, so pick most positive instead. + auto try_prove = [&](bool check, const top2_t& direct_trk, const top2_t& anti_trk) { + if (!check || proven) return; + auto [yd, ycd] = pick_master(direct_trk, cand); + if (evaluate(cand_coeff, is_upward, yd, ycd, false)) { + proven = true; + master_col = yd; + is_anti = false; + return; + } + auto [ya, yca] = pick_master(anti_trk, cand); + if (evaluate(cand_coeff, is_upward, ya, yca, true)) { + proven = true; + master_col = ya; + is_anti = true; + return; + } + }; + try_prove(use_leq_check, is_upward ? neg_y : pos_y, is_upward ? pos_y : neg_y); + try_prove(use_geq_check, is_upward ? pos_y : neg_y, is_upward ? neg_y : pos_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. + // For direct: favorable y=1 (upward) or y=0 (downward). + // For anti (complement): favorable state flips. + f_t y_coef_val = dense_row_coefs[master_col]; + f_t fav_y_contrib = (is_upward != is_anti) ? y_coef_val : f_t{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 ? num.isFeasGT(fav, bound) : num.isFeasLT(fav, bound)) proven = false; + }; + check_side( + has_rhs, can_reach_pos_inf, A_max, std::max(f_t{0}, y_coef_val), rhs_values[row], UPPER); + check_side( + has_lhs, can_reach_neg_inf, A_min, std::min(f_t{0}, y_coef_val), lhs_values[row], LOWER); + if (!proven) continue; + + if (!substitution_numerically_stable(constraint_matrix, cand)) continue; + + substituted[cand] = true; + if (is_anti) + reductions.replaceCol(cand, master_col, f_t{-1}, f_t{1}); // x = 1 - y + else + reductions.replaceCol(cand, master_col, f_t{1}, f_t{0}); // x = y + ++n_substitutions; + } + + for (int j = 0; j < length; ++j) + dense_row_coefs[cols[j]] = 0; + + return n_substitutions; +} + +} // namespace + +// ========================================================================= +// Top-level entry point +// ========================================================================= + +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 int ncols = problem.getNCols(); + const double tlim = problemUpdate.getPresolveOptions().tlim; + + std::vector locks[2], lock_row[2]; + compute_single_locks(problem, timer, tlim, locks, lock_row); + + auto candidates = collect_candidates(problem, locks, lock_row); + + if (this->is_time_exceeded(timer, tlim) || candidates.empty()) + return papilo::PresolveStatus::kUnchanged; + + // Well, technically O(K log K). But could be O(K) if it made a difference + std::sort(candidates.begin(), candidates.end(), [](const candidate_t& a, const candidate_t& b) { + return a.locking_row < b.locking_row; + }); + + int n_substitutions = 0; + std::vector dense_row_coefs(ncols, f_t{0}); + std::vector substituted(ncols, 0); + + 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; }); + + n_substitutions += try_substitutions_for_row( + problem, num, reductions, cand_it, row_end, r, dense_row_coefs, substituted); + + 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 || PDLP_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 0000000000..1d047922ee --- /dev/null +++ b/cpp/src/mip_heuristics/presolve/single_lock_dual_aggregation.hpp @@ -0,0 +1,42 @@ +/* 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; +}; + +} // 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 fae143de29..2264316a3f 100644 --- a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp +++ b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp @@ -30,6 +30,7 @@ #endif #include #include +#include #include #include #include @@ -541,6 +542,8 @@ void set_presolve_methods(papilo::Presolve& presolver, presolver.addPresolveMethod(uptr(new papilo::SimpleSubstitution())); presolver.addPresolveMethod(uptr(new papilo::Sparsify())); presolver.addPresolveMethod(uptr(new papilo::Substitution())); + presolver.addPresolveMethod( + uptr(new cuopt::linear_programming::detail::SingleLockDualAggregation())); } else { CUOPT_LOG_INFO("Disabling the presolver methods that do not support dual postsolve"); } diff --git a/cpp/src/pdlp/termination_strategy/infeasibility_information.cu b/cpp/src/pdlp/termination_strategy/infeasibility_information.cu index dbb35b732d..0e001b802f 100644 --- a/cpp/src/pdlp/termination_strategy/infeasibility_information.cu +++ b/cpp/src/pdlp/termination_strategy/infeasibility_information.cu @@ -24,6 +24,8 @@ #include #include +#include + namespace cuopt::linear_programming::detail { template infeasibility_information_t::infeasibility_information_t( diff --git a/cpp/tests/mip/CMakeLists.txt b/cpp/tests/mip/CMakeLists.txt index 2f2139890f..c992719bb4 100644 --- a/cpp/tests/mip/CMakeLists.txt +++ b/cpp/tests/mip/CMakeLists.txt @@ -37,7 +37,7 @@ ConfigureTest(EMPTY_FIXED_PROBLEMS_TEST ${CMAKE_CURRENT_SOURCE_DIR}/empty_fixed_problems_test.cu ) ConfigureTest(PRESOLVE_TEST - ${CMAKE_CURRENT_SOURCE_DIR}/presolve_test.cu + ${CMAKE_CURRENT_SOURCE_DIR}/presolve_test.cpp ) # Disable for now # ConfigureTest(FEASIBILITY_JUMP_TEST diff --git a/cpp/tests/mip/presolve_test.cpp b/cpp/tests/mip/presolve_test.cpp new file mode 100644 index 0000000000..02e4959be4 --- /dev/null +++ b/cpp/tests/mip/presolve_test.cpp @@ -0,0 +1,217 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#include + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include + +namespace cuopt::linear_programming::detail::test { + +struct papilo_harness_t { + papilo::Num num; + papilo::Reductions reductions; + papilo::Statistics stats; + papilo::PresolveOptions options; + papilo::Message msg; + double timer_acc{0}; + + papilo_harness_t() { options.tlim = std::numeric_limits::max(); } + + papilo::PresolveStatus run(papilo::Problem& problem, + SingleLockDualAggregation& presolver) + { + papilo::PostsolveStorage postsolve(problem.getNRows(), problem.getNCols()); + papilo::ProblemUpdate update(problem, postsolve, stats, options, num, msg); + papilo::Timer timer(timer_acc); + int reason = 0; + return presolver.execute(problem, update, num, reductions, timer, reason); + } + + bool has_replace(int col1, int col2, double factor, double offset) const + { + const auto& txns = reductions.getTransactions(); + const auto& reds = reductions.getReductions(); + for (const auto& t : txns) { + if (t.end - t.start != 2) continue; + const auto& r0 = reds[t.start]; + const auto& r1 = reds[t.start + 1]; + if (r0.row == papilo::ColReduction::REPLACE && r0.col == col1 && + std::abs(r0.newval - factor) < 1e-9 && r1.col == col2 && + std::abs(r1.newval - offset) < 1e-9) + return true; + } + return false; + } +}; + +papilo::Problem build_problem(int nrows, + int ncols, + const std::vector>& entries, + const std::vector& obj, + const std::vector& lb, + const std::vector& ub, + const std::vector& integer, + const std::vector& row_lhs, + const std::vector& row_rhs, + const std::vector& lhs_inf, + const std::vector& rhs_inf) +{ + papilo::ProblemBuilder builder; + builder.setNumRows(nrows); + builder.setNumCols(ncols); + for (auto& [r, c, v] : entries) + builder.addEntry(r, c, v); + for (int c = 0; c < ncols; ++c) { + builder.setObj(c, obj[c]); + builder.setColLb(c, lb[c]); + builder.setColUb(c, ub[c]); + builder.setColIntegral(c, integer[c]); + } + for (int r = 0; r < nrows; ++r) { + builder.setRowLhs(r, row_lhs[r]); + builder.setRowRhs(r, row_rhs[r]); + builder.setRowLhsInf(r, lhs_inf[r]); + builder.setRowRhsInf(r, rhs_inf[r]); + } + return builder.build(); +} + +// x has one up-lock in a LEQ row. Probe proves y=0 => x=0 via activity. +// Favorable-state check passes. Result: x = y (direct substitution). +// +// min -x +// s.t. 3x - 4y <= 1 (the locking row) +// x + y >= 0 (GEQ slack: positive coeffs add down-locks only) +// x, y in {0,1} +// +// On the locking row: +// A_min=-4, A_max=3. Probe(x=1,y=0): probed_min = -4-0-(-4)+3 = 3 > 1 => proven. +// Favorable(y=1): A_max - max(0,-4) + (-4) = 3-0-4 = -1 <= 1 => safe. +TEST(SingleLockDualAggregation, DirectSubstitution) +{ + auto problem = build_problem(2, + 2, + {{0, 0, 3.0}, {0, 1, -4.0}, {1, 0, 1.0}, {1, 1, 1.0}}, + {-1.0, 0.0}, + {0.0, 0.0}, + {1.0, 1.0}, + {true, true}, + {0.0, 0.0}, + {1.0, 0.0}, + {true, false}, + {false, true}); + + SingleLockDualAggregation presolver; + papilo_harness_t h; + auto status = h.run(problem, presolver); + + EXPECT_EQ(status, papilo::PresolveStatus::kReduced); + EXPECT_TRUE(h.has_replace(0, 1, 1.0, 0.0)); // x = y +} + +// Direct master has no negative binary coeff, so direct probe fails. +// Anti probe (y=1 unfavorable for upward) succeeds. Result: x = 1 - y. +// +// min -x +// s.t. 3x + 4y <= 5 (the locking row) +// x + y >= 0 (GEQ slack: positive coeffs add down-locks only) +// x, y in {0,1} +// +// On the locking row: +// A_min=0, A_max=7. Direct: no neg binary master. Anti master: y (coeff +4). +// Probe(x=1,y=1): probed_min = 0-0-0+(3+4) = 7 > 5 => proven. +// Favorable(y=0 for anti-upward): A_max - max(0,4) + 0 = 7-4 = 3 <= 5 => safe. +TEST(SingleLockDualAggregation, AntiSubstitution) +{ + auto problem = build_problem(2, + 2, + {{0, 0, 3.0}, {0, 1, 4.0}, {1, 0, 1.0}, {1, 1, 1.0}}, + {-1.0, 0.0}, + {0.0, 0.0}, + {1.0, 1.0}, + {true, true}, + {0.0, 0.0}, + {5.0, 0.0}, + {true, false}, + {false, true}); + + SingleLockDualAggregation presolver; + papilo_harness_t h; + auto status = h.run(problem, presolver); + + EXPECT_EQ(status, papilo::PresolveStatus::kReduced); + EXPECT_TRUE(h.has_replace(0, 1, -1.0, 1.0)); // x = 1 - y +} + +// A free row (both sides infinite) produces no locks and no candidates. +TEST(SingleLockDualAggregation, FreeRowNonCase) +{ + auto problem = build_problem(1, + 2, + {{0, 0, 2.0}, {0, 1, 3.0}}, + {-1.0, 0.0}, + {0.0, 0.0}, + {1.0, 1.0}, + {true, true}, + {0.0}, + {0.0}, + {true}, + {true}); + + SingleLockDualAggregation presolver; + papilo_harness_t h; + auto status = h.run(problem, presolver); + + EXPECT_EQ(status, papilo::PresolveStatus::kUnchanged); + EXPECT_EQ(h.reductions.size(), 0u); +} + +// Probe proves the implication but the favorable-state check fails, +// so the substitution is correctly rejected. +// +// min -x +// s.t. 3x - 2y <= 0 +// x, y in {0,1} +// +// A_min=-2, A_max=3. Probe(x=1,y=0): probed_min = -2-0-(-2)+3 = 3 > 0 => proven. +// Favorable(y=1): A_max - max(0,-2) + (-2) = 3-0-2 = 1 > 0 => FAILS. +TEST(SingleLockDualAggregation, FavorableStateRejects) +{ + auto problem = build_problem(1, + 2, + {{0, 0, 3.0}, {0, 1, -2.0}}, + {-1.0, 0.0}, + {0.0, 0.0}, + {1.0, 1.0}, + {true, true}, + {0.0}, + {0.0}, + {true}, + {false}); + + SingleLockDualAggregation presolver; + papilo_harness_t h; + auto status = h.run(problem, presolver); + + EXPECT_EQ(status, papilo::PresolveStatus::kUnchanged); + EXPECT_EQ(h.reductions.size(), 0u); +} + +} // namespace cuopt::linear_programming::detail::test diff --git a/cpp/tests/mip/presolve_test.cu b/cpp/tests/mip/presolve_test.cu deleted file mode 100644 index cf8abd1b69..0000000000 --- a/cpp/tests/mip/presolve_test.cu +++ /dev/null @@ -1,65 +0,0 @@ -/* clang-format off */ -/* - * SPDX-FileCopyrightText: Copyright (c) 2024-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. - * SPDX-License-Identifier: Apache-2.0 - */ -/* clang-format on */ - -#include "../linear_programming/utilities/pdlp_test_utilities.cuh" - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include - -#include - -#include -#include -#include -#include - -namespace cuopt::linear_programming::test { - -TEST(problem, find_implied_integers) -{ - const raft::handle_t handle_{}; - - auto path = make_path_absolute("mip/fiball.mps"); - auto mps_data_model = cuopt::mps_parser::parse_mps(path, false); - auto op_problem = mps_data_model_to_optimization_problem(&handle_, mps_data_model); - auto presolver = std::make_unique>(); - auto result = presolver->apply(op_problem, - cuopt::linear_programming::problem_category_t::MIP, - cuopt::linear_programming::presolver_t::Papilo, - false, - 1e-6, - 1e-12, - 20, - 1); - ASSERT_TRUE(result.has_value()); - - auto problem = detail::problem_t(result->reduced_problem); - problem.set_implied_integers(result->implied_integer_indices); - ASSERT_TRUE(result->implied_integer_indices.size() > 0); - auto var_types = host_copy(problem.variable_types, handle_.get_stream()); - // Find the index of the one continuous variable - auto it = std::find_if(var_types.begin(), var_types.end(), [](var_t var_type) { - return var_type == var_t::CONTINUOUS; - }); - ASSERT_NE(it, var_types.end()); - ASSERT_EQ(problem.presolve_data.var_flags.size(), var_types.size()); - // Ensure it is an implied integer - EXPECT_EQ(problem.presolve_data.var_flags.element(it - var_types.begin(), handle_.get_stream()), - ((int)detail::problem_t::var_flags_t::VAR_IMPLIED_INTEGER)); -} - -} // namespace cuopt::linear_programming::test diff --git a/cpp/tests/mip/problem_test.cu b/cpp/tests/mip/problem_test.cu index 92fa6d41d1..f925bf8558 100644 --- a/cpp/tests/mip/problem_test.cu +++ b/cpp/tests/mip/problem_test.cu @@ -8,6 +8,7 @@ #include "../linear_programming/utilities/pdlp_test_utilities.cuh" #include +#include #include #include #include @@ -484,4 +485,37 @@ TEST(optimization_problem_t_DeathTest, test_check_problem_validity) } #endif +TEST(problem, find_implied_integers) +{ + const raft::handle_t handle_{}; + + auto path = make_path_absolute("mip/fiball.mps"); + auto mps_data_model = cuopt::mps_parser::parse_mps(path, false); + auto op_problem = mps_data_model_to_optimization_problem(&handle_, mps_data_model); + auto presolver = std::make_unique>(); + auto result = presolver->apply(op_problem, + cuopt::linear_programming::problem_category_t::MIP, + cuopt::linear_programming::presolver_t::Papilo, + false, + 1e-6, + 1e-12, + 20, + 1); + ASSERT_TRUE(result.has_value()); + + auto problem = detail::problem_t(result->reduced_problem); + problem.set_implied_integers(result->implied_integer_indices); + ASSERT_TRUE(result->implied_integer_indices.size() > 0); + auto var_types = host_copy(problem.variable_types, handle_.get_stream()); + // Find the index of the one continuous variable + auto it = std::find_if(var_types.begin(), var_types.end(), [](var_t var_type) { + return var_type == var_t::CONTINUOUS; + }); + ASSERT_NE(it, var_types.end()); + ASSERT_EQ(problem.presolve_data.var_flags.size(), var_types.size()); + // Ensure it is an implied integer + EXPECT_EQ(problem.presolve_data.var_flags.element(it - var_types.begin(), handle_.get_stream()), + ((int)detail::problem_t::var_flags_t::VAR_IMPLIED_INTEGER)); +} + } // namespace cuopt::linear_programming::test diff --git a/cpp/tests/mip/termination_test.cu b/cpp/tests/mip/termination_test.cu index 5f21d294f9..b0bf32e344 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 3040f0f543..483d06814c 100755 --- a/datasets/mip/download_miplib_test_dataset.sh +++ b/datasets/mip/download_miplib_test_dataset.sh @@ -24,6 +24,7 @@ INSTANCES=( "swath1" "enlight_hard" "enlight11" + "neos-787933" "supportcase22" )