From 61b1e6a92199499d54b76870edd3f435e6f79c3f Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Tue, 4 Aug 2026 06:30:30 -0700 Subject: [PATCH 1/3] fix gf2 underdetemrined handling, support rectangular cases, more tests --- .../mip_heuristics/presolve/gf2_presolve.cpp | 183 ++++++++++++------ .../presolve/third_party_presolve.cpp | 91 +++++---- .../presolve/third_party_presolve.hpp | 10 + cpp/tests/mip/presolve_test.cu | 154 +++++++++++++++ 4 files changed, 348 insertions(+), 90 deletions(-) diff --git a/cpp/src/mip_heuristics/presolve/gf2_presolve.cpp b/cpp/src/mip_heuristics/presolve/gf2_presolve.cpp index c984ef306c..0da0b97af1 100644 --- a/cpp/src/mip_heuristics/presolve/gf2_presolve.cpp +++ b/cpp/src/mip_heuristics/presolve/gf2_presolve.cpp @@ -8,8 +8,10 @@ #include "gf2_presolve.hpp" #include +#include #include +#include #include #if GF2_PRESOLVE_DEBUG @@ -33,55 +35,106 @@ static inline i_t positive_modulo(i_t i, i_t n) return (i % n + n) % n; } +enum class gf2_status_t { Feasible, Infeasible }; + +static constexpr int GF2_WORD_BITS = 64; + +static inline int gf2_nwords(int N) { return (N + GF2_WORD_BITS - 1) / GF2_WORD_BITS; } + +static inline bool gf2_test_bit(const std::vector& row, int col) +{ + return (row[col / GF2_WORD_BITS] >> (col % GF2_WORD_BITS)) & uint64_t{1}; +} + +static inline void gf2_set_bit(std::vector& row, int col) +{ + row[col / GF2_WORD_BITS] |= (uint64_t{1} << (col % GF2_WORD_BITS)); +} + // this is kind-of a stopgap implementation (as in practice MIPLIB2017 only contains a couple of GF2 // problems and they're small) but cuDSS could be used for this since A is likely to be sparse and // low-bandwidth (i think?) unlikely to occur in real-world problems however. doubt it'd be worth -// the effort trashes A and b, return true if solved -static bool gf2_solve(std::vector>& A, std::vector& b, std::vector& x) +// the effort +// trashes A and b. A is m x n (m rows, each packed over n columns). +static gf2_status_t gf2_solve(std::vector>& A, + int n_cols, + std::vector& b, + std::vector& x, + std::vector& determined) { - int i, j, k; - const int N = A.size(); - for (i = 0; i < N; i++) { - // Find pivot + const int m = (int)A.size(); + const int n = n_cols; + const int nwords = gf2_nwords(n); + cuopt_assert(m > 0, ""); + cuopt_assert(n >= 0, ""); + cuopt_assert((int)b.size() == m, ""); + cuopt_assert((int)A[0].size() == nwords, ""); + + // pivot_row_of_col[c] = row holding the pivot for column c, or -1 if free + std::vector pivot_row_of_col(n, -1); + int next_pivot_row = 0; + + for (int col = 0; col < n; col++) { int pivot = -1; - for (j = i; j < N; j++) { - if (A[j][i]) { - pivot = j; + for (int r = next_pivot_row; r < m; r++) { + if (gf2_test_bit(A[r], col)) { + pivot = r; break; } } - if (pivot == -1) return false; // No solution - - // Swap current row with pivot row if needed - if (pivot != i) { - for (k = 0; k < N; k++) { - int temp = A[i][k]; - A[i][k] = A[pivot][k]; - A[pivot][k] = temp; - } - int temp = b[i]; - b[i] = b[pivot]; - b[pivot] = temp; + if (pivot == -1) continue; // free column + + if (pivot != next_pivot_row) { + std::swap(A[next_pivot_row], A[pivot]); + std::swap(b[next_pivot_row], b[pivot]); } - // Eliminate downwards - for (j = i + 1; j < N; j++) { - if (A[j][i]) { - for (k = i; k < N; k++) - A[j][k] ^= A[i][k]; - b[j] ^= b[i]; + // Eliminate column from all other rows (RREF) + for (int r = 0; r < m; r++) { + if (r != next_pivot_row && gf2_test_bit(A[r], col)) { + for (int w = 0; w < nwords; w++) + A[r][w] ^= A[next_pivot_row][w]; + b[r] ^= b[next_pivot_row]; } } + + pivot_row_of_col[col] = next_pivot_row; + next_pivot_row++; + } + + const int rank = next_pivot_row; + for (int r = rank; r < m; r++) { + for (int w = 0; w < nwords; w++) { + cuopt_assert(A[r][w] == 0, "RREF unused row must be zero"); + } + if (b[r]) return gf2_status_t::Infeasible; + } + + std::vector free_mask(nwords, 0); + for (int c = 0; c < n; c++) { + if (pivot_row_of_col[c] == -1) gf2_set_bit(free_mask, c); } - // Back-substitution - for (i = N - 1; i >= 0; i--) { - x[i] = b[i]; - for (j = i + 1; j < N; j++) - x[i] ^= (A[i][j] & x[j]); - if (!A[i][i] && x[i]) return false; // No solution + determined.assign(n, 0); + x.assign(n, 0); + + for (int col = 0; col < n; col++) { + int row = pivot_row_of_col[col]; + if (row == -1) continue; // free: x=0, determined=false + + bool has_free_support = false; + for (int w = 0; w < nwords; w++) { + if (A[row][w] & free_mask[w]) { + has_free_support = true; + break; + } + } + // Particular solution with free vars = 0: x[pivot] = b[row] + x[col] = b[row]; + determined[col] = !has_free_support; } - return true; // Success + + return gf2_status_t::Feasible; } template @@ -161,16 +214,20 @@ papilo::PresolveStatus GF2Presolve::execute(const papilo::Problem& pro if (key_var_idx != -1) { NOT_GF2("multiple key variables", var_idx); } key_var_idx = var_idx; key_var_coeff = coeff; - gf2_key_vars.insert({var_idx, gf2_key_vars.size()}); } else { // Binary variable constraint_bin_vars.push_back({var_idx, coeff}); - gf2_bin_vars.insert({var_idx, gf2_bin_vars.size()}); } } if (key_var_idx == -1) NOT_GF2("missing key variable"); + // Commit to global maps only after the row is fully accepted + gf2_key_vars.insert({(size_t)key_var_idx, gf2_key_vars.size()}); + for (auto [bin_var, _] : constraint_bin_vars) { + gf2_bin_vars.insert({bin_var, gf2_bin_vars.size()}); + } + gf2_constraints.emplace_back((size_t)cstr_idx, std::move(constraint_bin_vars), std::pair{key_var_idx, key_var_coeff}, @@ -183,12 +240,11 @@ papilo::PresolveStatus GF2Presolve::execute(const papilo::Problem& pro // If no GF2 constraints found, return unchanged if (gf2_constraints.empty()) { return papilo::PresolveStatus::kUnchanged; } - // Skip if that would cause computational explosion (O(n^3) with simple gaussian elimination) - if (gf2_constraints.size() > 1000) { return papilo::PresolveStatus::kUnchanged; } + // one unique key per GF2 row. #bins may differ from #rows. + if (gf2_key_vars.size() != gf2_constraints.size()) { return papilo::PresolveStatus::kUnchanged; } - // Validate structure - if (gf2_key_vars.size() != gf2_constraints.size() || - gf2_bin_vars.size() != gf2_constraints.size()) { + // Skip if that would cause computational explosion (dense GE ~ O(m * n * min(m,n))) + if (gf2_constraints.size() > 1000 || gf2_bin_vars.size() > 1000) { return papilo::PresolveStatus::kUnchanged; } @@ -198,40 +254,57 @@ papilo::PresolveStatus GF2Presolve::execute(const papilo::Problem& pro gf2_bin_vars_invmap.insert({gf2_idx, var_idx}); } - // Build binary matrix - // Could be a flat vector but. oh well. in practice N is small - std::vector> A(gf2_constraints.size(), - std::vector(gf2_constraints.size(), 0)); - std::vector b(gf2_constraints.size()); - for (size_t gf2_cstr_idx = 0; gf2_cstr_idx < gf2_constraints.size(); ++gf2_cstr_idx) { + // Build binary matrix as packed uint64_t words + const int m = (int)gf2_constraints.size(); + const int n = (int)gf2_bin_vars.size(); + const int nwords = gf2_nwords(n); + std::vector> A(m, std::vector(nwords, 0)); + std::vector b(m); + for (int gf2_cstr_idx = 0; gf2_cstr_idx < m; ++gf2_cstr_idx) { const auto& cons = gf2_constraints[gf2_cstr_idx]; for (auto [bin_var, _] : cons.bin_vars) { - A[gf2_cstr_idx][gf2_bin_vars[bin_var]] = 1; + gf2_set_bit(A[gf2_cstr_idx], (int)gf2_bin_vars[bin_var]); } b[gf2_cstr_idx] = cons.rhs; } - std::vector solution(gf2_constraints.size()); - bool feasible = gf2_solve(A, b, solution); - if (!feasible) { return papilo::PresolveStatus::kInfeasible; } + std::vector solution(n); + std::vector determined(n); + gf2_status_t gf2_status = gf2_solve(A, n, b, solution, determined); + if (gf2_status == gf2_status_t::Infeasible) { return papilo::PresolveStatus::kInfeasible; } std::unordered_map fixings; - // Fix binary variables - for (size_t sol_idx = 0; sol_idx < gf2_constraints.size(); ++sol_idx) { - fixings[gf2_bin_vars_invmap[sol_idx]] = solution[sol_idx]; + + // Fix only uniquely determined binaries + for (int sol_idx = 0; sol_idx < n; ++sol_idx) { + if (determined[sol_idx]) { fixings[gf2_bin_vars_invmap[sol_idx]] = solution[sol_idx]; } } - // Compute fixings for key variables by solving for the constraint + // Fix key only when every binary in that constraint is uniquely determined for (const auto& cons : gf2_constraints) { + bool all_bins_determined = true; + for (auto [bin_var, _] : cons.bin_vars) { + cuopt_assert(gf2_bin_vars.count(bin_var), ""); + if (!determined[gf2_bin_vars[bin_var]]) { + all_bins_determined = false; + break; + } + } + if (!all_bins_determined) continue; + auto [key_var_idx, key_var_coeff] = cons.key_var; f_t constraint_rhs = lhs_values[cons.cstr_idx]; // equality constraint f_t lhs = -constraint_rhs; for (auto [bin_var, coeff] : cons.bin_vars) { + cuopt_assert(fixings.count(bin_var), ""); lhs += fixings[bin_var] * coeff; } fixings[key_var_idx] = std::round(-lhs / key_var_coeff); } + // necessary because Papilo asserts on empty TransactionGuard + if (fixings.empty()) { return papilo::PresolveStatus::kUnchanged; } + papilo::PresolveStatus status = papilo::PresolveStatus::kUnchanged; papilo::TransactionGuard rg{reductions}; for (const auto& [var_idx, fixing] : fixings) { diff --git a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp index dc8410e7a1..6a5bd341fb 100644 --- a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp +++ b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp @@ -657,46 +657,55 @@ void check_postsolve_status(const papilo::PostsolveStatus& status) } template -void set_presolve_methods(papilo::Presolve& presolver, - problem_category_t category, - bool dual_postsolve) +void set_presolve_methods( + papilo::Presolve& presolver, + problem_category_t category, + bool dual_postsolve, + std::optional> const& method_allowlist = std::nullopt) { using uptr = std::unique_ptr>; + auto maybe_add = [&](uptr method) { + if (method_allowlist.has_value()) { + const std::string& name = method->getName(); + if (!method_allowlist->count(name)) { return; } + } + presolver.addPresolveMethod(std::move(method)); + }; + if (category == problem_category_t::MIP) { // cuOpt custom GF2 presolver - presolver.addPresolveMethod( - uptr(new cuopt::mathematical_optimization::mip::GF2Presolve())); + maybe_add(uptr(new cuopt::mathematical_optimization::mip::GF2Presolve())); } // fast presolvers - presolver.addPresolveMethod(uptr(new papilo::SingletonCols())); - presolver.addPresolveMethod(uptr(new papilo::CoefficientStrengthening())); - presolver.addPresolveMethod(uptr(new papilo::ConstraintPropagation())); + maybe_add(uptr(new papilo::SingletonCols())); + maybe_add(uptr(new papilo::CoefficientStrengthening())); + maybe_add(uptr(new papilo::ConstraintPropagation())); // medium presolvers - presolver.addPresolveMethod(uptr(new papilo::FixContinuous())); - presolver.addPresolveMethod(uptr(new papilo::SimpleProbing())); - presolver.addPresolveMethod(uptr(new papilo::ParallelRowDetection())); - presolver.addPresolveMethod(uptr(new papilo::ParallelColDetection())); - presolver.addPresolveMethod(uptr(new papilo::DualFix())); - presolver.addPresolveMethod(uptr(new papilo::SimplifyInequalities())); - presolver.addPresolveMethod(uptr(new papilo::CliqueMerging())); + maybe_add(uptr(new papilo::FixContinuous())); + maybe_add(uptr(new papilo::SimpleProbing())); + maybe_add(uptr(new papilo::ParallelRowDetection())); + maybe_add(uptr(new papilo::ParallelColDetection())); + maybe_add(uptr(new papilo::DualFix())); + maybe_add(uptr(new papilo::SimplifyInequalities())); + maybe_add(uptr(new papilo::CliqueMerging())); // exhaustive presolvers - presolver.addPresolveMethod(uptr(new papilo::ImplIntDetection())); - presolver.addPresolveMethod(uptr(new papilo::DominatedCols())); - presolver.addPresolveMethod(uptr(new papilo::Probing())); + maybe_add(uptr(new papilo::ImplIntDetection())); + maybe_add(uptr(new papilo::DominatedCols())); + maybe_add(uptr(new papilo::Probing())); if (!dual_postsolve) { // SingletonStuffing causes dual crushing failures on: // tr12-30, ns1208400, gmu-35-50, dws008-01, neos-1445765, // neos-5107597-kakapo, rocI-4-11, traininstance2, traininstance6, // radiationm18-12-05, rococoB10-011000, b1c1s1 - presolver.addPresolveMethod(uptr(new papilo::SingletonStuffing())); - presolver.addPresolveMethod(uptr(new papilo::DualInfer())); - presolver.addPresolveMethod(uptr(new papilo::SimpleSubstitution())); - presolver.addPresolveMethod(uptr(new papilo::Sparsify())); - presolver.addPresolveMethod(uptr(new papilo::Substitution())); + maybe_add(uptr(new papilo::SingletonStuffing())); + maybe_add(uptr(new papilo::DualInfer())); + maybe_add(uptr(new papilo::SimpleSubstitution())); + maybe_add(uptr(new papilo::Sparsify())); + maybe_add(uptr(new papilo::Substitution())); } else { CUOPT_LOG_INFO("Disabling the presolver methods that do not support dual postsolve"); } @@ -721,22 +730,31 @@ void set_presolve_options(papilo::Presolve& presolver, } template -void set_presolve_parameters(papilo::Presolve& presolver, - problem_category_t category, - int nrows, - int ncols) +void set_presolve_parameters( + papilo::Presolve& presolver, + problem_category_t category, + int nrows, + int ncols, + std::optional> const& method_allowlist = std::nullopt) { // It looks like a copy. But this copy has the pointers to relevant variables in papilo auto params = presolver.getParameters(); if (category == problem_category_t::MIP) { + auto reduction_allowed = [&](char const* name) { + return !method_allowlist.has_value() || method_allowlist->count(name) > 0; + }; // Papilo has work unit measurements for probing. Because of this when the first batch fails to // produce any reductions, the algorithm stops. To avoid stopping the algorithm, we set a // minimum badge size to a huge value. The time limit makes sure that we exit if it takes too // long - int min_badgesize = std::max(ncols / 2, 32); - params.setParameter("probing.minbadgesize", min_badgesize); - params.setParameter("cliquemerging.enabled", true); - params.setParameter("cliquemerging.maxcalls", 50); + if (reduction_allowed("probing")) { + int min_badgesize = std::max(ncols / 2, 32); + params.setParameter("probing.minbadgesize", min_badgesize); + } + if (reduction_allowed("cliquemerging")) { + params.setParameter("cliquemerging.enabled", true); + params.setParameter("cliquemerging.maxcalls", 50); + } } } @@ -835,7 +853,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); + set_presolve_methods(papilo_presolver, category, dual_postsolve, reduction_allowlist_); set_presolve_options(papilo_presolver, category, absolute_tolerance, @@ -843,7 +861,8 @@ third_party_presolve_status_t third_party_presolve_t::apply_papilo( time_limit, dual_postsolve, num_cpu_threads); - set_presolve_parameters(papilo_presolver, category, original_n_cons, original_n_vars); + set_presolve_parameters( + papilo_presolver, category, original_n_cons, original_n_vars, reduction_allowlist_); papilo_presolver.setVerbosityLevel(papilo::VerbosityLevel::kQuiet); auto result = papilo_presolver.apply(papilo_problem); @@ -1067,7 +1086,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); + 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, @@ -1075,7 +1095,8 @@ third_party_presolve_status_t third_party_presolve_t::apply_to_subprob time_limit, dual_postsolve, num_threads); - set_presolve_parameters(papilo_presolver, problem_category_t::MIP, orig_rows, orig_cols); + set_presolve_parameters( + papilo_presolver, problem_category_t::MIP, orig_rows, orig_cols, reduction_allowlist_); // Disable papilo logs papilo_presolver.setVerbosityLevel(papilo::VerbosityLevel::kQuiet); diff --git a/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp b/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp index 7ed62ef07c..5eec96221d 100644 --- a/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp +++ b/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp @@ -9,6 +9,8 @@ #include #include +#include +#include #include #include @@ -96,6 +98,12 @@ class third_party_presolve_t { double time_limit, i_t num_cpu_threads = 0); + // If set, only Papilo methods whose getName() is listed are registered + void set_reduction_allowlist(std::optional> allowlist) + { + reduction_allowlist_ = std::move(allowlist); + } + // 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( @@ -196,6 +204,8 @@ class third_party_presolve_t { std::vector original_objective_coefficients_{}; f_t original_objective_offset_{0}; f_t original_objective_scaling_factor_{1}; + + std::optional> reduction_allowlist_{}; }; // Just for testing the conversion: user_problem -> Papilo problem -> user_problem. diff --git a/cpp/tests/mip/presolve_test.cu b/cpp/tests/mip/presolve_test.cu index 7067479093..6a2d760446 100644 --- a/cpp/tests/mip/presolve_test.cu +++ b/cpp/tests/mip/presolve_test.cu @@ -16,6 +16,7 @@ #include #include #include +#include #include #include @@ -151,4 +152,157 @@ TEST(gf2_presolve, uses_compact_constraint_indices) EXPECT_EQ(result.status, mip::third_party_presolve_status_t::REDUCED); } +static mip::third_party_presolve_device_result_t run_gf2_presolve( + std::string_view lp_text) +{ + const raft::handle_t handle{}; + auto mps_data_model = cuopt::test::parse_inline_lp(lp_text); + auto op_problem = mps_data_model_to_optimization_problem(&handle, mps_data_model); + auto presolver = std::make_unique>(); + presolver->set_reduction_allowlist(std::unordered_set{"gf2presolve"}); + return presolver->apply_presolve_from_op_problem( + op_problem, problem_category_t::MIP, presolver_t::Papilo, false, 1e-6, 1e-12, 20, 1); +} + +// Consistent singular: both rows x⊕y = 1. Check we do not return infeasible. +TEST(gf2_presolve, consistent_singular_unchanged) +{ + auto result = run_gf2_presolve(R"LP( +Minimize + obj: x0 + x1 + k0 + k1 +Subject To + c0: x0 + x1 + 2 k0 = 1 + c1: x0 + x1 + 2 k1 = 1 +Binaries + x0 + x1 + k0 + k1 +End +)LP"); + EXPECT_EQ(result.status, mip::third_party_presolve_status_t::UNCHANGED); +} + +// Inconsistent singular: x⊕y = 1 and x⊕y = 0. +TEST(gf2_presolve, inconsistent_singular_infeasible) +{ + auto result = run_gf2_presolve(R"LP( +Minimize + obj: x0 + x1 + k0 + k1 +Subject To + c0: x0 + x1 + 2 k0 = 1 + c1: x0 + x1 + 2 k1 = 0 +Binaries + x0 + x1 + k0 + k1 +End +)LP"); + EXPECT_EQ(result.status, mip::third_party_presolve_status_t::INFEASIBLE); +} + +// Partially determined: y is unique over GF(2); x,z free with x⊕z = 1. +// Rows: [1,0,1], [0,1,0], [1,1,1] with rhs [1,1,0] (row2 = row0⊕row1). +TEST(gf2_presolve, partial_determination_reduces) +{ + auto result = run_gf2_presolve(R"LP( +Minimize + obj: x0 + x1 + x2 + k0 + k1 + k2 +Subject To + c0: x0 - x2 + 2 k0 = 1 + c1: x1 + 2 k1 = 1 + c2: x0 + x1 - x2 + 2 k2 = 0 +Binaries + x0 + x1 + x2 + k0 + k1 + k2 +End +)LP"); + EXPECT_EQ(result.status, mip::third_party_presolve_status_t::REDUCED); +} + +// Fat (n > m): x0 fixed; x1⊕x2 free. +TEST(gf2_presolve, more_bins_than_rows_reduces) +{ + auto result = run_gf2_presolve(R"LP( +Minimize + obj: x0 + x1 + x2 + k0 + k1 +Subject To + c0: x0 + 2 k0 = 1 + c1: x1 + x2 + 2 k1 = 1 +Binaries + x0 + x1 + x2 + k0 + k1 +End +)LP"); + EXPECT_EQ(result.status, mip::third_party_presolve_status_t::REDUCED); +} + +// Tall consistent (m > n): x0=1, x1=0 uniquely (third row redundant). +TEST(gf2_presolve, more_rows_than_bins_reduces) +{ + auto result = run_gf2_presolve(R"LP( +Minimize + obj: x0 + x1 + k0 + k1 + k2 +Subject To + c0: x0 + 2 k0 = 1 + c1: x0 + x1 + 2 k1 = 1 + c2: x1 + 2 k2 = 0 +Binaries + x0 + x1 + k0 + k1 + k2 +End +)LP"); + EXPECT_EQ(result.status, mip::third_party_presolve_status_t::REDUCED); +} + +// Tall inconsistent (m > n): x0 = 1 and x0 = 0. +TEST(gf2_presolve, more_rows_than_bins_infeasible) +{ + auto result = run_gf2_presolve(R"LP( +Minimize + obj: x0 + k0 + k1 +Subject To + c0: x0 + 2 k0 = 1 + c1: x0 + 2 k1 = 0 +Binaries + x0 + k0 + k1 +End +)LP"); + EXPECT_EQ(result.status, mip::third_party_presolve_status_t::INFEASIBLE); +} + +// Near-miss row must not leak key/bin vars into the maps and suppress a valid GF2 reduction. +TEST(gf2_presolve, near_miss_row_does_not_suppress_reduction) +{ + auto result = run_gf2_presolve(R"LP( +Minimize + obj: x0 + x1 + k0 + k_bad + w +Subject To + c0: x0 + 2 k0 = 1 + c_bad: x0 + x1 + 2 k_bad + 3 w = 1 +Binaries + x0 + x1 + k0 + k_bad +Generals + w +End +)LP"); + EXPECT_EQ(result.status, mip::third_party_presolve_status_t::REDUCED); +} + } // namespace cuopt::mathematical_optimization::test From 41d7566d429dc0bedc2058a2bf95bf25a9e3d3d5 Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Tue, 4 Aug 2026 06:52:42 -0700 Subject: [PATCH 2/3] more gf2 tests --- .../mip_heuristics/presolve/gf2_presolve.cpp | 13 +- .../mip_heuristics/presolve/gf2_presolve.hpp | 14 + cpp/tests/internal/CMakeLists.txt | 1 + cpp/tests/mip/gf2_presolve_test.cpp | 561 ++++++++++++++++++ cpp/tests/mip/presolve_test.cu | 240 -------- 5 files changed, 581 insertions(+), 248 deletions(-) create mode 100644 cpp/tests/mip/gf2_presolve_test.cpp diff --git a/cpp/src/mip_heuristics/presolve/gf2_presolve.cpp b/cpp/src/mip_heuristics/presolve/gf2_presolve.cpp index 0da0b97af1..8ee3ed4877 100644 --- a/cpp/src/mip_heuristics/presolve/gf2_presolve.cpp +++ b/cpp/src/mip_heuristics/presolve/gf2_presolve.cpp @@ -35,8 +35,6 @@ static inline i_t positive_modulo(i_t i, i_t n) return (i % n + n) % n; } -enum class gf2_status_t { Feasible, Infeasible }; - static constexpr int GF2_WORD_BITS = 64; static inline int gf2_nwords(int N) { return (N + GF2_WORD_BITS - 1) / GF2_WORD_BITS; } @@ -55,12 +53,11 @@ static inline void gf2_set_bit(std::vector& row, int col) // problems and they're small) but cuDSS could be used for this since A is likely to be sparse and // low-bandwidth (i think?) unlikely to occur in real-world problems however. doubt it'd be worth // the effort -// trashes A and b. A is m x n (m rows, each packed over n columns). -static gf2_status_t gf2_solve(std::vector>& A, - int n_cols, - std::vector& b, - std::vector& x, - std::vector& determined) +gf2_status_t gf2_solve(std::vector>& A, + int n_cols, + std::vector& b, + std::vector& x, + std::vector& determined) { const int m = (int)A.size(); const int n = n_cols; diff --git a/cpp/src/mip_heuristics/presolve/gf2_presolve.hpp b/cpp/src/mip_heuristics/presolve/gf2_presolve.hpp index 1ae4bd5f06..6bac9c5f08 100644 --- a/cpp/src/mip_heuristics/presolve/gf2_presolve.hpp +++ b/cpp/src/mip_heuristics/presolve/gf2_presolve.hpp @@ -20,8 +20,22 @@ #pragma GCC diagnostic pop #endif +#include +#include + namespace cuopt::mathematical_optimization::mip { +enum class gf2_status_t { Feasible, Infeasible }; + +// Solves A x = b over GF(2). A is m x n, each row packed into ceil(n/64) words (column c lives at +// word c/64, bit c%64). Trashes A and b. On Feasible, x is the solution obtained by setting the +// free variables to 0, and determined[c] is set iff x[c] is the same in every solution. +gf2_status_t gf2_solve(std::vector>& A, + int n_cols, + std::vector& b, + std::vector& x, + std::vector& determined); + template class GF2Presolve : public papilo::PresolveMethod { public: diff --git a/cpp/tests/internal/CMakeLists.txt b/cpp/tests/internal/CMakeLists.txt index c580e0117a..2856c56f38 100644 --- a/cpp/tests/internal/CMakeLists.txt +++ b/cpp/tests/internal/CMakeLists.txt @@ -27,6 +27,7 @@ ConfigureTest(NUMOPT_INTERNAL_TEST ${CUOPT_TEST_DIR}/mip/integer_with_real_bounds.cu ${CUOPT_TEST_DIR}/mip/empty_fixed_problems_test.cu ${CUOPT_TEST_DIR}/mip/presolve_test.cu + ${CUOPT_TEST_DIR}/mip/gf2_presolve_test.cpp ${CUOPT_TEST_DIR}/mip/termination_test.cu ${CUOPT_TEST_DIR}/mip/determinism_test.cu # socp diff --git a/cpp/tests/mip/gf2_presolve_test.cpp b/cpp/tests/mip/gf2_presolve_test.cpp new file mode 100644 index 0000000000..2957ee6267 --- /dev/null +++ b/cpp/tests/mip/gf2_presolve_test.cpp @@ -0,0 +1,561 @@ +/* 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 +#include +#include +#include +#include + +namespace cuopt::mathematical_optimization::test { + +using mip::gf2_solve; +using mip::gf2_status_t; + +namespace { + +// A GF(2) system in a form that stays readable as test data: one '0'/'1' string per row. +struct gf2_system_t { + std::vector rows; + std::vector rhs; + + int n_cols() const { return rows.empty() ? 0 : (int)rows[0].size(); } + int n_rows() const { return (int)rows.size(); } +}; + +// Packed without going through gf2_set_bit, so a layout mistake in the solver shows up as a wrong +// result instead of cancelling out against a shared helper. +std::vector> pack(const gf2_system_t& system) +{ + const int n = system.n_cols(); + const int nwords = (n + 63) / 64; + std::vector> A(system.n_rows(), std::vector(nwords, 0)); + for (int r = 0; r < system.n_rows(); ++r) { + for (int c = 0; c < n; ++c) { + if (system.rows[r][c] == '1') { A[r][c >> 6] |= uint64_t{1} << (c & 63); } + } + } + return A; +} + +std::string describe(const gf2_system_t& system) +{ + std::string out = "system:"; + for (int r = 0; r < system.n_rows(); ++r) { + out += "\n [" + system.rows[r] + "] = " + std::to_string(system.rhs[r]); + } + return out; +} + +bool satisfies(const gf2_system_t& system, const std::vector& x) +{ + for (int r = 0; r < system.n_rows(); ++r) { + int parity = 0; + for (int c = 0; c < system.n_cols(); ++c) { + if (system.rows[r][c] == '1') { parity ^= x[c] & 1; } + } + if (parity != (system.rhs[r] & 1)) { return false; } + } + return true; +} + +// Every x in {0,1}^n satisfying the system, evaluated on the unpacked form. +std::vector> all_solutions(const gf2_system_t& system) +{ + const int n = system.n_cols(); + std::vector> solutions; + for (uint64_t assignment = 0; assignment < (uint64_t{1} << n); ++assignment) { + std::vector x(n); + for (int c = 0; c < n; ++c) { + x[c] = (assignment >> c) & 1; + } + if (satisfies(system, x)) { solutions.push_back(std::move(x)); } + } + return solutions; +} + +// Checks gf2_solve against exhaustive enumeration and returns the number of solutions. +// +// Both directions of the determinedness contract matter: reporting a varying column as determined +// would fix a variable the presolver has no right to fix, and reporting a constant column as free +// silently loses a reduction. +size_t check_against_enumeration(const gf2_system_t& system) +{ + const int n = system.n_cols(); + auto A = pack(system); + auto b = system.rhs; + std::vector x; + std::vector determined; + const gf2_status_t status = gf2_solve(A, n, b, x, determined); + + const auto solutions = all_solutions(system); + if (solutions.empty()) { + EXPECT_EQ(status, gf2_status_t::Infeasible) << describe(system); + return 0; + } + + EXPECT_EQ(status, gf2_status_t::Feasible) << describe(system); + EXPECT_EQ((int)x.size(), n) << describe(system); + EXPECT_EQ((int)determined.size(), n) << describe(system); + if (status != gf2_status_t::Feasible || (int)x.size() != n || (int)determined.size() != n) { + return solutions.size(); + } + + EXPECT_TRUE(satisfies(system, x)) << describe(system) << "\nreported x is not a solution"; + + for (int c = 0; c < n; ++c) { + const bool constant = + std::all_of(solutions.begin(), solutions.end(), [&](const std::vector& solution) { + return solution[c] == solutions[0][c]; + }); + EXPECT_EQ(determined[c] != 0, constant) << describe(system) << "\ncolumn " << c; + if (constant) { EXPECT_EQ(x[c], solutions[0][c]) << describe(system) << "\ncolumn " << c; } + } + return solutions.size(); +} + +// Builds a random system. Rows are often XORs of earlier rows: a uniformly random GF(2) matrix is +// almost always full rank, which is the case rank-deficiency handling is least concerned with. +gf2_system_t random_system(std::mt19937& rng) +{ + const int m = 1 + (int)(rng() % 6); + const int n = (int)(rng() % 11); + const int density = 15 + (int)(rng() % 71); + + gf2_system_t system; + system.rows.assign(m, std::string(n, '0')); + system.rhs.assign(m, 0); + + for (int r = 0; r < m; ++r) { + if (r > 0 && (rng() % 2) == 0) { + system.rows[r] = system.rows[rng() % r]; + if ((rng() % 2) == 0) { + const std::string& other = system.rows[rng() % r]; + for (int c = 0; c < n; ++c) { + system.rows[r][c] = (char)('0' + ((system.rows[r][c] - '0') ^ (other[c] - '0'))); + } + } + } else { + for (int c = 0; c < n; ++c) { + system.rows[r][c] = ((int)(rng() % 100) < density) ? '1' : '0'; + } + } + } + + // Half the systems get a planted solution so the feasible path stays well covered; the rest get + // a random rhs, which is usually inconsistent once the rows are dependent. + if ((rng() % 2) == 0) { + std::vector planted(n); + for (int c = 0; c < n; ++c) { + planted[c] = (int)(rng() % 2); + } + for (int r = 0; r < m; ++r) { + int parity = 0; + for (int c = 0; c < n; ++c) { + if (system.rows[r][c] == '1') { parity ^= planted[c]; } + } + system.rhs[r] = parity; + } + } else { + for (int r = 0; r < m; ++r) { + system.rhs[r] = (int)(rng() % 2); + } + } + return system; +} + +struct gf2_golden_case_t { + const char* name; + gf2_system_t system; + gf2_status_t status; + // One character per column: '0'/'1' where the column is uniquely determined, '.' where it is + // not. x is unconstrained on '.' columns, so the test must not pin it there. + const char* expected; +}; + +} // namespace + +TEST(gf2_solve, golden_cases) +{ + const std::vector cases = { + {"identity", {{"100", "010", "001"}, {1, 0, 1}}, gf2_status_t::Feasible, "101"}, + {"full_rank_mixed", {{"110", "011", "001"}, {1, 1, 0}}, gf2_status_t::Feasible, "010"}, + // Duplicate rows: rank 1 of 2, and column 2 appears in no row. + {"consistent_singular", {{"110", "110"}, {1, 1}}, gf2_status_t::Feasible, "..."}, + {"inconsistent_singular", {{"110", "110"}, {1, 0}}, gf2_status_t::Infeasible, ""}, + // More rows than columns, third row redundant after elimination. + {"tall_consistent", {{"10", "11", "01"}, {1, 1, 0}}, gf2_status_t::Feasible, "10"}, + {"tall_inconsistent", {{"1", "1"}, {1, 0}}, gf2_status_t::Infeasible, ""}, + // More columns than rows. + {"fat", {{"100", "011"}, {1, 1}}, gf2_status_t::Feasible, "1.."}, + // Rows carrying only a key variable reach gf2_solve with no columns at all. + {"no_columns_consistent", {{""}, {0}}, gf2_status_t::Feasible, ""}, + {"no_columns_inconsistent", {{""}, {1}}, gf2_status_t::Infeasible, ""}, + // Column 2 is all zero, as happens for a binary left in the map by a rejected row. + {"zero_column", {{"100", "010"}, {1, 1}}, gf2_status_t::Feasible, "11."}, + // The free column (1) sits below the determined pivot column (2). + {"free_col_before_pivot", {{"110", "001"}, {1, 1}}, gf2_status_t::Feasible, "..1"}, + // J - I at even dimension is nonsingular over GF(2); the shape the enlight instances hit. + {"j_minus_i_4", + {{"0111", "1011", "1101", "1110"}, {1, 1, 1, 1}}, + gf2_status_t::Feasible, + "1111"}, + }; + + for (const auto& test_case : cases) { + SCOPED_TRACE(test_case.name); + + auto A = pack(test_case.system); + auto b = test_case.system.rhs; + std::vector x; + std::vector determined; + const gf2_status_t status = gf2_solve(A, test_case.system.n_cols(), b, x, determined); + + EXPECT_EQ(status, test_case.status) << describe(test_case.system); + + if (test_case.status == gf2_status_t::Feasible && status == gf2_status_t::Feasible) { + const std::string expected{test_case.expected}; + ASSERT_EQ(x.size(), expected.size()) << describe(test_case.system); + ASSERT_EQ(determined.size(), expected.size()) << describe(test_case.system); + for (size_t c = 0; c < expected.size(); ++c) { + if (expected[c] == '.') { + EXPECT_EQ(determined[c], 0) << describe(test_case.system) << "\ncolumn " << c; + } else { + EXPECT_NE(determined[c], 0) << describe(test_case.system) << "\ncolumn " << c; + EXPECT_EQ(x[c], expected[c] - '0') << describe(test_case.system) << "\ncolumn " << c; + } + } + } + + // Double entry: the hand-written expectations above must also agree with enumeration. + check_against_enumeration(test_case.system); + } +} + +TEST(gf2_solve, matches_enumeration_on_random_systems) +{ + std::mt19937 rng{20260804u}; + int infeasible = 0; + int unique = 0; + int underdetermined = 0; + + for (int iteration = 0; iteration < 1500; ++iteration) { + SCOPED_TRACE(iteration); + const size_t n_solutions = check_against_enumeration(random_system(rng)); + if (n_solutions == 0) { + infeasible++; + } else if (n_solutions == 1) { + unique++; + } else { + underdetermined++; + } + } + + // A generator drifting into one bucket would gut the test without failing it. + EXPECT_GT(infeasible, 100); + EXPECT_GT(unique, 50); + EXPECT_GT(underdetermined, 200); +} + +// Guards against mixing up the row and column spaces, which is the failure mode the m x n +// generalization introduces. determined is a property of the solution set, so it must permute with +// the columns. x must not be compared on undetermined columns: it holds b[pivot_row] there, which +// legitimately depends on which column won the pivot. +TEST(gf2_solve, determinedness_permutes_with_the_columns) +{ + std::mt19937 rng{20260805u}; + + for (int iteration = 0; iteration < 300; ++iteration) { + SCOPED_TRACE(iteration); + const gf2_system_t system = random_system(rng); + const int n = system.n_cols(); + const int m = system.n_rows(); + + auto A = pack(system); + auto b = system.rhs; + std::vector x; + std::vector determined; + const gf2_status_t status = gf2_solve(A, n, b, x, determined); + + std::vector row_perm(m); + std::vector col_perm(n); + for (int r = 0; r < m; ++r) { + row_perm[r] = r; + } + for (int c = 0; c < n; ++c) { + col_perm[c] = c; + } + std::shuffle(row_perm.begin(), row_perm.end(), rng); + std::shuffle(col_perm.begin(), col_perm.end(), rng); + + gf2_system_t permuted; + permuted.rows.assign(m, std::string(n, '0')); + permuted.rhs.assign(m, 0); + for (int r = 0; r < m; ++r) { + for (int c = 0; c < n; ++c) { + permuted.rows[r][c] = system.rows[row_perm[r]][col_perm[c]]; + } + permuted.rhs[r] = system.rhs[row_perm[r]]; + } + + auto permuted_A = pack(permuted); + auto permuted_b = permuted.rhs; + std::vector permuted_x; + std::vector permuted_determined; + const gf2_status_t permuted_status = + gf2_solve(permuted_A, n, permuted_b, permuted_x, permuted_determined); + + ASSERT_EQ(status, permuted_status) << describe(system); + if (status != gf2_status_t::Feasible) { continue; } + + EXPECT_TRUE(satisfies(permuted, permuted_x)) << describe(permuted); + for (int c = 0; c < n; ++c) { + EXPECT_EQ(permuted_determined[c] != 0, determined[col_perm[c]] != 0) + << describe(system) << "\ncolumn " << c; + if (determined[col_perm[c]]) { + EXPECT_EQ(permuted_x[c], x[col_perm[c]]) << describe(system) << "\ncolumn " << c; + } + } + } +} + +namespace { + +mip::third_party_presolve_device_result_t run_gf2_presolve(std::string_view lp_text) +{ + const raft::handle_t handle{}; + auto mps_data_model = cuopt::test::parse_inline_lp(lp_text); + auto op_problem = mps_data_model_to_optimization_problem(&handle, mps_data_model); + auto presolver = std::make_unique>(); + presolver->set_reduction_allowlist(std::unordered_set{"gf2presolve"}); + return presolver->apply_presolve_from_op_problem( + op_problem, problem_category_t::MIP, presolver_t::Papilo, false, 1e-6, 1e-12, 20, 1); +} + +} // namespace + +TEST(gf2_presolve, uses_compact_constraint_indices) +{ + constexpr int num_packing_vars = 128; + constexpr int num_gf2_vars = 6; + constexpr int num_key_vars = 6; + constexpr int num_packing_rows = 128; + constexpr int num_key_rows = 2 * num_key_vars; + constexpr int num_gf2_rows = 6; + constexpr int num_vars = num_packing_vars + num_gf2_vars + num_key_vars; + constexpr int num_rows = num_packing_rows + num_key_rows + num_gf2_rows; + constexpr int x_offset = num_packing_vars; + constexpr int y_offset = x_offset + num_gf2_vars; + + std::vector values; + std::vector indices; + std::vector offsets{0}; + std::vector constraint_lb(num_rows, 1.0); + std::vector constraint_ub(num_rows, 2.0); + + auto add_entry = [&](int column, double value) { + indices.push_back(column); + values.push_back(value); + }; + auto finish_row = [&] { offsets.push_back(static_cast(values.size())); }; + + // A normal binary MIP block keeps the GF2 rows at high raw row indices. + for (int row = 0; row < num_packing_rows; ++row) { + std::array columns{row, (row + 1) % num_packing_vars, (row + 2) % num_packing_vars}; + std::sort(columns.begin(), columns.end()); + for (int column : columns) { + add_entry(column, 1.0); + } + finish_row(); + } + + // Keep every GF2 key column non-singleton without forcing it. + for (int key = 0; key < num_key_vars; ++key) { + add_entry(3 * key, 1.0); + add_entry(3 * key + 1, 1.0); + add_entry(y_offset + key, 1.0); + finish_row(); + + add_entry(3 * key + 1, 1.0); + add_entry(3 * key + 2, 1.0); + add_entry(y_offset + key, 1.0); + finish_row(); + } + + // Over GF(2), this is J-I for even dimension 6, hence nonsingular. Three positive and two + // negative coefficients per row prevent ordinary bound propagation from fixing the key. + for (int row = 0; row < num_gf2_rows; ++row) { + int term = 0; + for (int col = 0; col < num_gf2_vars; ++col) { + if (col == row) { continue; } + add_entry(x_offset + col, term < 3 ? 1.0 : -1.0); + ++term; + } + add_entry(y_offset + row, 2.0); + finish_row(); + constraint_lb[num_packing_rows + num_key_rows + row] = 1.0; + constraint_ub[num_packing_rows + num_key_rows + row] = 1.0; + } + + const raft::handle_t handle_{}; + optimization_problem_t problem(&handle_); + std::vector objective(num_vars, 1.0); + std::vector variable_lb(num_vars, 0.0); + std::vector variable_ub(num_vars, 1.0); + std::vector variable_types(num_vars, var_t::INTEGER); + problem.set_csr_constraint_matrix( + values.data(), values.size(), indices.data(), indices.size(), offsets.data(), offsets.size()); + problem.set_objective_coefficients(objective.data(), objective.size()); + problem.set_variable_lower_bounds(variable_lb.data(), variable_lb.size()); + problem.set_variable_upper_bounds(variable_ub.data(), variable_ub.size()); + problem.set_variable_types(variable_types.data(), variable_types.size()); + problem.set_constraint_lower_bounds(constraint_lb.data(), constraint_lb.size()); + problem.set_constraint_upper_bounds(constraint_ub.data(), constraint_ub.size()); + + auto presolver = std::make_unique>(); + auto result = presolver->apply_presolve_from_op_problem( + problem, problem_category_t::MIP, presolver_t::Papilo, false, 1e-6, 1e-12, 20, 1); + + EXPECT_EQ(result.status, mip::third_party_presolve_status_t::REDUCED); +} + +// Consistent singular: both rows x0 xor x1 = 1. Check we do not return infeasible. +TEST(gf2_presolve, consistent_singular_unchanged) +{ + auto result = run_gf2_presolve(R"LP( +Minimize + obj: x0 + x1 + k0 + k1 +Subject To + c0: x0 + x1 + 2 k0 = 1 + c1: x0 + x1 + 2 k1 = 1 +Binaries + x0 + x1 + k0 + k1 +End +)LP"); + EXPECT_EQ(result.status, mip::third_party_presolve_status_t::UNCHANGED); +} + +// Inconsistent singular: x0 xor x1 = 1 and x0 xor x1 = 0. +TEST(gf2_presolve, inconsistent_singular_infeasible) +{ + auto result = run_gf2_presolve(R"LP( +Minimize + obj: x0 + x1 + k0 + k1 +Subject To + c0: x0 + x1 + 2 k0 = 1 + c1: x0 + x1 + 2 k1 = 0 +Binaries + x0 + x1 + k0 + k1 +End +)LP"); + EXPECT_EQ(result.status, mip::third_party_presolve_status_t::INFEASIBLE); +} + +// Partially determined: x1 is unique over GF(2); x0, x2 free with x0 xor x2 = 1. +TEST(gf2_presolve, partial_determination_reduces) +{ + auto result = run_gf2_presolve(R"LP( +Minimize + obj: x0 + x1 + x2 + k0 + k1 + k2 +Subject To + c0: x0 - x2 + 2 k0 = 1 + c1: x1 + 2 k1 = 1 + c2: x0 + x1 - x2 + 2 k2 = 0 +Binaries + x0 + x1 + x2 + k0 + k1 + k2 +End +)LP"); + EXPECT_EQ(result.status, mip::third_party_presolve_status_t::REDUCED); +} + +// Fat (n > m): x0 fixed; x1 xor x2 free. +TEST(gf2_presolve, more_bins_than_rows_reduces) +{ + auto result = run_gf2_presolve(R"LP( +Minimize + obj: x0 + x1 + x2 + k0 + k1 +Subject To + c0: x0 + 2 k0 = 1 + c1: x1 + x2 + 2 k1 = 1 +Binaries + x0 + x1 + x2 + k0 + k1 +End +)LP"); + EXPECT_EQ(result.status, mip::third_party_presolve_status_t::REDUCED); +} + +// Tall consistent (m > n): x0 = 1, x1 = 0 uniquely, third row redundant. +TEST(gf2_presolve, more_rows_than_bins_reduces) +{ + auto result = run_gf2_presolve(R"LP( +Minimize + obj: x0 + x1 + k0 + k1 + k2 +Subject To + c0: x0 + 2 k0 = 1 + c1: x0 + x1 + 2 k1 = 1 + c2: x1 + 2 k2 = 0 +Binaries + x0 + x1 + k0 + k1 + k2 +End +)LP"); + EXPECT_EQ(result.status, mip::third_party_presolve_status_t::OPTIMAL); +} + +// Tall inconsistent (m > n): x0 = 1 and x0 = 0. +TEST(gf2_presolve, more_rows_than_bins_infeasible) +{ + auto result = run_gf2_presolve(R"LP( +Minimize + obj: x0 + k0 + k1 +Subject To + c0: x0 + 2 k0 = 1 + c1: x0 + 2 k1 = 0 +Binaries + x0 + k0 + k1 +End +)LP"); + EXPECT_EQ(result.status, mip::third_party_presolve_status_t::INFEASIBLE); +} + +} // namespace cuopt::mathematical_optimization::test diff --git a/cpp/tests/mip/presolve_test.cu b/cpp/tests/mip/presolve_test.cu index 6a2d760446..3ac4b1d9f3 100644 --- a/cpp/tests/mip/presolve_test.cu +++ b/cpp/tests/mip/presolve_test.cu @@ -16,7 +16,6 @@ #include #include #include -#include #include #include @@ -24,7 +23,6 @@ #include #include -#include #include #include #include @@ -67,242 +65,4 @@ TEST(problem, find_implied_integers) ((int)mip::problem_t::var_flags_t::VAR_IMPLIED_INTEGER)); } -TEST(gf2_presolve, uses_compact_constraint_indices) -{ - constexpr int num_packing_vars = 128; - constexpr int num_gf2_vars = 6; - constexpr int num_key_vars = 6; - constexpr int num_packing_rows = 128; - constexpr int num_key_rows = 2 * num_key_vars; - constexpr int num_gf2_rows = 6; - constexpr int num_vars = num_packing_vars + num_gf2_vars + num_key_vars; - constexpr int num_rows = num_packing_rows + num_key_rows + num_gf2_rows; - constexpr int x_offset = num_packing_vars; - constexpr int y_offset = x_offset + num_gf2_vars; - - std::vector values; - std::vector indices; - std::vector offsets{0}; - std::vector constraint_lb(num_rows, 1.0); - std::vector constraint_ub(num_rows, 2.0); - - auto add_entry = [&](int column, double value) { - indices.push_back(column); - values.push_back(value); - }; - auto finish_row = [&] { offsets.push_back(static_cast(values.size())); }; - - // A normal binary MIP block keeps the GF2 rows at high raw row indices. - for (int row = 0; row < num_packing_rows; ++row) { - std::array columns{row, (row + 1) % num_packing_vars, (row + 2) % num_packing_vars}; - std::sort(columns.begin(), columns.end()); - for (int column : columns) { - add_entry(column, 1.0); - } - finish_row(); - } - - // Keep every GF2 key column non-singleton without forcing it. - for (int key = 0; key < num_key_vars; ++key) { - add_entry(3 * key, 1.0); - add_entry(3 * key + 1, 1.0); - add_entry(y_offset + key, 1.0); - finish_row(); - - add_entry(3 * key + 1, 1.0); - add_entry(3 * key + 2, 1.0); - add_entry(y_offset + key, 1.0); - finish_row(); - } - - // Over GF(2), this is J-I for even dimension 6, hence nonsingular. Three positive and two - // negative coefficients per row prevent ordinary bound propagation from fixing the key. - for (int row = 0; row < num_gf2_rows; ++row) { - int term = 0; - for (int col = 0; col < num_gf2_vars; ++col) { - if (col == row) { continue; } - add_entry(x_offset + col, term < 3 ? 1.0 : -1.0); - ++term; - } - add_entry(y_offset + row, 2.0); - finish_row(); - constraint_lb[num_packing_rows + num_key_rows + row] = 1.0; - constraint_ub[num_packing_rows + num_key_rows + row] = 1.0; - } - - const raft::handle_t handle_{}; - optimization_problem_t problem(&handle_); - std::vector objective(num_vars, 1.0); - std::vector variable_lb(num_vars, 0.0); - std::vector variable_ub(num_vars, 1.0); - std::vector variable_types(num_vars, var_t::INTEGER); - problem.set_csr_constraint_matrix( - values.data(), values.size(), indices.data(), indices.size(), offsets.data(), offsets.size()); - problem.set_objective_coefficients(objective.data(), objective.size()); - problem.set_variable_lower_bounds(variable_lb.data(), variable_lb.size()); - problem.set_variable_upper_bounds(variable_ub.data(), variable_ub.size()); - problem.set_variable_types(variable_types.data(), variable_types.size()); - problem.set_constraint_lower_bounds(constraint_lb.data(), constraint_lb.size()); - problem.set_constraint_upper_bounds(constraint_ub.data(), constraint_ub.size()); - - auto presolver = std::make_unique>(); - auto result = presolver->apply_presolve_from_op_problem( - problem, problem_category_t::MIP, presolver_t::Papilo, false, 1e-6, 1e-12, 20, 1); - - EXPECT_EQ(result.status, mip::third_party_presolve_status_t::REDUCED); -} - -static mip::third_party_presolve_device_result_t run_gf2_presolve( - std::string_view lp_text) -{ - const raft::handle_t handle{}; - auto mps_data_model = cuopt::test::parse_inline_lp(lp_text); - auto op_problem = mps_data_model_to_optimization_problem(&handle, mps_data_model); - auto presolver = std::make_unique>(); - presolver->set_reduction_allowlist(std::unordered_set{"gf2presolve"}); - return presolver->apply_presolve_from_op_problem( - op_problem, problem_category_t::MIP, presolver_t::Papilo, false, 1e-6, 1e-12, 20, 1); -} - -// Consistent singular: both rows x⊕y = 1. Check we do not return infeasible. -TEST(gf2_presolve, consistent_singular_unchanged) -{ - auto result = run_gf2_presolve(R"LP( -Minimize - obj: x0 + x1 + k0 + k1 -Subject To - c0: x0 + x1 + 2 k0 = 1 - c1: x0 + x1 + 2 k1 = 1 -Binaries - x0 - x1 - k0 - k1 -End -)LP"); - EXPECT_EQ(result.status, mip::third_party_presolve_status_t::UNCHANGED); -} - -// Inconsistent singular: x⊕y = 1 and x⊕y = 0. -TEST(gf2_presolve, inconsistent_singular_infeasible) -{ - auto result = run_gf2_presolve(R"LP( -Minimize - obj: x0 + x1 + k0 + k1 -Subject To - c0: x0 + x1 + 2 k0 = 1 - c1: x0 + x1 + 2 k1 = 0 -Binaries - x0 - x1 - k0 - k1 -End -)LP"); - EXPECT_EQ(result.status, mip::third_party_presolve_status_t::INFEASIBLE); -} - -// Partially determined: y is unique over GF(2); x,z free with x⊕z = 1. -// Rows: [1,0,1], [0,1,0], [1,1,1] with rhs [1,1,0] (row2 = row0⊕row1). -TEST(gf2_presolve, partial_determination_reduces) -{ - auto result = run_gf2_presolve(R"LP( -Minimize - obj: x0 + x1 + x2 + k0 + k1 + k2 -Subject To - c0: x0 - x2 + 2 k0 = 1 - c1: x1 + 2 k1 = 1 - c2: x0 + x1 - x2 + 2 k2 = 0 -Binaries - x0 - x1 - x2 - k0 - k1 - k2 -End -)LP"); - EXPECT_EQ(result.status, mip::third_party_presolve_status_t::REDUCED); -} - -// Fat (n > m): x0 fixed; x1⊕x2 free. -TEST(gf2_presolve, more_bins_than_rows_reduces) -{ - auto result = run_gf2_presolve(R"LP( -Minimize - obj: x0 + x1 + x2 + k0 + k1 -Subject To - c0: x0 + 2 k0 = 1 - c1: x1 + x2 + 2 k1 = 1 -Binaries - x0 - x1 - x2 - k0 - k1 -End -)LP"); - EXPECT_EQ(result.status, mip::third_party_presolve_status_t::REDUCED); -} - -// Tall consistent (m > n): x0=1, x1=0 uniquely (third row redundant). -TEST(gf2_presolve, more_rows_than_bins_reduces) -{ - auto result = run_gf2_presolve(R"LP( -Minimize - obj: x0 + x1 + k0 + k1 + k2 -Subject To - c0: x0 + 2 k0 = 1 - c1: x0 + x1 + 2 k1 = 1 - c2: x1 + 2 k2 = 0 -Binaries - x0 - x1 - k0 - k1 - k2 -End -)LP"); - EXPECT_EQ(result.status, mip::third_party_presolve_status_t::REDUCED); -} - -// Tall inconsistent (m > n): x0 = 1 and x0 = 0. -TEST(gf2_presolve, more_rows_than_bins_infeasible) -{ - auto result = run_gf2_presolve(R"LP( -Minimize - obj: x0 + k0 + k1 -Subject To - c0: x0 + 2 k0 = 1 - c1: x0 + 2 k1 = 0 -Binaries - x0 - k0 - k1 -End -)LP"); - EXPECT_EQ(result.status, mip::third_party_presolve_status_t::INFEASIBLE); -} - -// Near-miss row must not leak key/bin vars into the maps and suppress a valid GF2 reduction. -TEST(gf2_presolve, near_miss_row_does_not_suppress_reduction) -{ - auto result = run_gf2_presolve(R"LP( -Minimize - obj: x0 + x1 + k0 + k_bad + w -Subject To - c0: x0 + 2 k0 = 1 - c_bad: x0 + x1 + 2 k_bad + 3 w = 1 -Binaries - x0 - x1 - k0 - k_bad -Generals - w -End -)LP"); - EXPECT_EQ(result.status, mip::third_party_presolve_status_t::REDUCED); -} - } // namespace cuopt::mathematical_optimization::test From 1e6415de20716006065489d737d32efc65c4e02c Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Tue, 4 Aug 2026 10:33:00 -0700 Subject: [PATCH 3/3] check the key value for adequacy --- .../mip_heuristics/presolve/gf2_presolve.cpp | 23 +++++++++- cpp/tests/mip/gf2_presolve_test.cpp | 42 ++++++++++++++++++- 2 files changed, 62 insertions(+), 3 deletions(-) diff --git a/cpp/src/mip_heuristics/presolve/gf2_presolve.cpp b/cpp/src/mip_heuristics/presolve/gf2_presolve.cpp index 8ee3ed4877..8060f534ea 100644 --- a/cpp/src/mip_heuristics/presolve/gf2_presolve.cpp +++ b/cpp/src/mip_heuristics/presolve/gf2_presolve.cpp @@ -290,13 +290,32 @@ papilo::PresolveStatus GF2Presolve::execute(const papilo::Problem& pro if (!all_bins_determined) continue; auto [key_var_idx, key_var_coeff] = cons.key_var; - f_t constraint_rhs = lhs_values[cons.cstr_idx]; // equality constraint + const f_t constraint_rhs = std::round(lhs_values[cons.cstr_idx]); f_t lhs = -constraint_rhs; for (auto [bin_var, coeff] : cons.bin_vars) { cuopt_assert(fixings.count(bin_var), ""); lhs += fixings[bin_var] * coeff; } - fixings[key_var_idx] = std::round(-lhs / key_var_coeff); + const f_t key_val = std::round(-lhs / key_var_coeff); + + // Residual must be exactly 0 after rounding (rejects half-integer / inconsistent carry) + if (!num.isEq(lhs + key_val * key_var_coeff, f_t{0})) { + return papilo::PresolveStatus::kInfeasible; + } + // Dual-role: same var already fixed as a GF(2) binary + if (fixings.count(key_var_idx) && !num.isEq(fixings[key_var_idx], key_val)) { + return papilo::PresolveStatus::kInfeasible; + } + if (!col_flags[key_var_idx].test(papilo::ColFlag::kLbInf) && + key_val < lower_bounds[key_var_idx] - integrality_tolerance) { + return papilo::PresolveStatus::kInfeasible; + } + if (!col_flags[key_var_idx].test(papilo::ColFlag::kUbInf) && + key_val > upper_bounds[key_var_idx] + integrality_tolerance) { + return papilo::PresolveStatus::kInfeasible; + } + + fixings[key_var_idx] = key_val; } // necessary because Papilo asserts on empty TransactionGuard diff --git a/cpp/tests/mip/gf2_presolve_test.cpp b/cpp/tests/mip/gf2_presolve_test.cpp index 2957ee6267..f25b56f91d 100644 --- a/cpp/tests/mip/gf2_presolve_test.cpp +++ b/cpp/tests/mip/gf2_presolve_test.cpp @@ -433,7 +433,8 @@ TEST(gf2_presolve, uses_compact_constraint_indices) problem.set_constraint_upper_bounds(constraint_ub.data(), constraint_ub.size()); auto presolver = std::make_unique>(); - auto result = presolver->apply_presolve_from_op_problem( + presolver->set_reduction_allowlist(std::unordered_set{"gf2presolve"}); + auto result = presolver->apply_presolve_from_op_problem( problem, problem_category_t::MIP, presolver_t::Papilo, false, 1e-6, 1e-12, 20, 1); EXPECT_EQ(result.status, mip::third_party_presolve_status_t::REDUCED); @@ -558,4 +559,43 @@ End EXPECT_EQ(result.status, mip::third_party_presolve_status_t::INFEASIBLE); } +// Dual-role: k is key in c0 and a ±1 bin in c1. GF(2) forces k=1; ℤ key recovery wants k=0. +TEST(gf2_presolve, dual_role_key_bin_conflict_infeasible) +{ + auto result = run_gf2_presolve(R"LP( +Minimize + obj: x0 + k + y +Subject To + c0: x0 + 2 k = 1 + c1: k + 2 y = 1 +Binaries + x0 + k + y +End +)LP"); + EXPECT_EQ(result.status, mip::third_party_presolve_status_t::INFEASIBLE); +} + +// GF(2)-consistent with x0=x1=1, but key recovery gives k=-1 outside [0,1]. +TEST(gf2_presolve, key_out_of_bounds_infeasible) +{ + auto result = run_gf2_presolve(R"LP( +Minimize + obj: x0 + x1 + a + b + k +Subject To + c0: x0 + 2 a = 1 + c1: x1 + 2 b = 1 + c2: x0 + x1 + 2 k = 0 +Binaries + x0 + x1 + a + b + k +End +)LP"); + EXPECT_EQ(result.status, mip::third_party_presolve_status_t::INFEASIBLE); +} + } // namespace cuopt::mathematical_optimization::test