From a0fffe6aa143485afd22cb362ecc2e44cc059771 Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Fri, 19 Sep 2025 15:03:34 +0000 Subject: [PATCH 1/9] gf2 xor constraint detection --- .../linear_programming/cuopt/run_mip.cpp | 3 +- cpp/src/mip/CMakeLists.txt | 3 + .../mip/presolve/bounds_update_helpers.cuh | 4 +- cpp/src/mip/presolve/gf2_presolve.cu | 146 ++++++++++++++++++ cpp/src/mip/presolve/gf2_presolve.cuh | 32 ++++ cpp/src/mip/presolve/trivial_presolve.cuh | 5 + 6 files changed, 190 insertions(+), 3 deletions(-) create mode 100644 cpp/src/mip/presolve/gf2_presolve.cu create mode 100644 cpp/src/mip/presolve/gf2_presolve.cuh diff --git a/benchmarks/linear_programming/cuopt/run_mip.cpp b/benchmarks/linear_programming/cuopt/run_mip.cpp index 64b1bb4bb..904274736 100644 --- a/benchmarks/linear_programming/cuopt/run_mip.cpp +++ b/benchmarks/linear_programming/cuopt/run_mip.cpp @@ -210,7 +210,8 @@ int run_single_file(std::string file_path, settings.log_to_console = log_to_console; settings.tolerances.relative_tolerance = 1e-12; settings.tolerances.absolute_tolerance = 1e-6; - settings.presolve = true; + settings.presolve = false; + settings.mip_scaling = false; cuopt::linear_programming::benchmark_info_t benchmark_info; settings.benchmark_info_ptr = &benchmark_info; auto start_run_solver = std::chrono::high_resolution_clock::now(); diff --git a/cpp/src/mip/CMakeLists.txt b/cpp/src/mip/CMakeLists.txt index 8e72ca70e..02ce5944f 100644 --- a/cpp/src/mip/CMakeLists.txt +++ b/cpp/src/mip/CMakeLists.txt @@ -45,10 +45,13 @@ set(MIP_NON_LP_FILES ${CMAKE_CURRENT_SOURCE_DIR}/presolve/multi_probe.cu ${CMAKE_CURRENT_SOURCE_DIR}/presolve/probing_cache.cu ${CMAKE_CURRENT_SOURCE_DIR}/presolve/trivial_presolve.cu + ${CMAKE_CURRENT_SOURCE_DIR}/presolve/gf2_presolve.cu ${CMAKE_CURRENT_SOURCE_DIR}/feasibility_jump/feasibility_jump.cu ${CMAKE_CURRENT_SOURCE_DIR}/feasibility_jump/feasibility_jump_kernels.cu ) +set_source_files_properties(${CMAKE_CURRENT_SOURCE_DIR}/presolve/gf2_presolve.cu DIRECTORY ${CMAKE_SOURCE_DIR} PROPERTIES COMPILE_OPTIONS "-g;-O0") + # Choose which files to include based on build mode if(BUILD_LP_ONLY) set(MIP_SRC_FILES ${MIP_LP_NECESSARY_FILES}) diff --git a/cpp/src/mip/presolve/bounds_update_helpers.cuh b/cpp/src/mip/presolve/bounds_update_helpers.cuh index 77574d527..3de6b5537 100644 --- a/cpp/src/mip/presolve/bounds_update_helpers.cuh +++ b/cpp/src/mip/presolve/bounds_update_helpers.cuh @@ -24,13 +24,13 @@ namespace cuopt::linear_programming::detail { // Activity calculation template -inline __device__ f_t min_act_of_var(f_t coeff, f_t var_lb, f_t var_ub) +inline __host__ __device__ f_t min_act_of_var(f_t coeff, f_t var_lb, f_t var_ub) { return (coeff < 0.) ? coeff * var_ub : coeff * var_lb; } template -inline __device__ f_t max_act_of_var(f_t coeff, f_t var_lb, f_t var_ub) +inline __host__ __device__ f_t max_act_of_var(f_t coeff, f_t var_lb, f_t var_ub) { return (coeff < 0.) ? coeff * var_lb : coeff * var_ub; } diff --git a/cpp/src/mip/presolve/gf2_presolve.cu b/cpp/src/mip/presolve/gf2_presolve.cu new file mode 100644 index 000000000..62ff83542 --- /dev/null +++ b/cpp/src/mip/presolve/gf2_presolve.cu @@ -0,0 +1,146 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights + * reserved. SPDX-License-Identifier: Apache-2.0 + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include + +#include +#include + +namespace cuopt::linear_programming::detail { + +// Inspired from +// https://github.com/snowberryfield/printemps/blob/0dc1f9d2b78a667451c11d8ab156861ce54e9fb7/printemps/model_component/constraint.h +// (MIT License) + +#define NOT_GF2(reason, ...) \ + { \ + printf("NO : Cons %d is not gf2: " reason "\n", cstr_idx, ##__VA_ARGS__); \ + goto not_valid; \ + } + +// TODO: implement fully on GPU + cuDSS? +template +void gf2_presolve(problem_t& problem) +{ + auto h_constraint_lower_bounds = cuopt::host_copy(problem.constraint_lower_bounds); + auto h_constraint_upper_bounds = cuopt::host_copy(problem.constraint_upper_bounds); + auto h_coefficients = cuopt::host_copy(problem.coefficients); + auto h_variables = cuopt::host_copy(problem.variables); + auto h_is_binary_variable = cuopt::host_copy(problem.is_binary_variable); + auto h_variable_bounds = cuopt::host_copy(problem.variable_bounds); + auto h_offsets = cuopt::host_copy(problem.offsets); + auto h_variable_types = cuopt::host_copy(problem.variable_types); + + for (i_t cstr_idx = 0; cstr_idx < problem.n_constraints; ++cstr_idx) { + i_t key_var_idx = -1; + f_t key_var_coeff = 0.0; + f_t key_var_lb = 0.0; + f_t key_var_ub = 0.0; + f_t min_activity_without_key = 0.0; + f_t max_activity_without_key = 0.0; + + //__asm__ __volatile__("int3"); + + f_t rhs = round(h_constraint_lower_bounds[cstr_idx]); + + // needs to be an equality constraint + if (h_constraint_lower_bounds[cstr_idx] != h_constraint_upper_bounds[cstr_idx]) + NOT_GF2("not an equality constraint"); + if (!isfinite(h_constraint_lower_bounds[cstr_idx])) NOT_GF2("not finite"); + if (!is_integer(h_constraint_lower_bounds[cstr_idx], problem.tolerances.integrality_tolerance)) + NOT_GF2("rhs not integer %f", h_constraint_lower_bounds[cstr_idx]); + + // only accept 0, 1, -1 as rhs + if (rhs != 0.0 && rhs != 1.0 && rhs != -1.0) NOT_GF2("invalid rhs %f", rhs); + + // check if nnzs match the pattern of a gf2 constraint + for (i_t j = h_offsets[cstr_idx]; j < h_offsets[cstr_idx + 1]; ++j) { + if (!is_integer(h_coefficients[j], problem.tolerances.integrality_tolerance)) + NOT_GF2("coeff not integer (%d,%f)", j, h_coefficients[j]); + i_t var_idx = h_variables[j]; + f_t coeff = round(h_coefficients[j]); + if (h_variable_types[var_idx] != var_t::INTEGER) + NOT_GF2("var not integer (%d,%f)", var_idx, h_coefficients[j]); + bool is_binary = h_is_binary_variable[var_idx]; + + // only the key variable can have a coefficient of 2 + if (is_binary && (abs(coeff) != 1.0 && abs(coeff) != 2.0)) + NOT_GF2("invalid coef bin (%d,%f)", var_idx, coeff); + if (!is_binary && (abs(coeff) != 2.0)) + NOT_GF2("invalid coef non-bin (%d,%f)", var_idx, coeff); + + // key var + if (abs(coeff) == 2.0) { + // can only be one + if (key_var_idx != -1) NOT_GF2("multiple key variables"); + key_var_idx = var_idx; + key_var_coeff = coeff; + } else { + min_activity_without_key += min_act_of_var( + coeff, get_lower(h_variable_bounds[var_idx]), get_upper(h_variable_bounds[var_idx])); + max_activity_without_key += max_act_of_var( + coeff, get_lower(h_variable_bounds[var_idx]), get_upper(h_variable_bounds[var_idx])); + } + } + + // no key var found + if (key_var_idx == -1) NOT_GF2("no key variable"); + + if (key_var_coeff > 0) { + std::swap(min_activity_without_key, max_activity_without_key); + min_activity_without_key *= -1; + max_activity_without_key *= -1; + } + cuopt_assert(min_activity_without_key <= max_activity_without_key, "invalid activities?"); + + key_var_lb = get_lower(h_variable_bounds[key_var_idx]); + key_var_ub = get_upper(h_variable_bounds[key_var_idx]); + if (isfinite(key_var_lb) && key_var_lb > ceil(min_activity_without_key / 2.0)) + NOT_GF2("invalid key variable bounds lower ((%f,%f), %f)", + key_var_lb, + key_var_ub, + min_activity_without_key); + if (isfinite(key_var_ub) && key_var_ub < floor(max_activity_without_key / 2.0)) + NOT_GF2("invalid key variable bounds upper ((%f,%f), %f)", + key_var_lb, + key_var_ub, + max_activity_without_key); + + printf("YES: Cons %d is gf2\n", cstr_idx); + continue; + + not_valid: + continue; + } + + exit(0); +} + +#define INSTANTIATE(F_TYPE) \ + template void gf2_presolve(problem_t & problem); + +#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/presolve/gf2_presolve.cuh b/cpp/src/mip/presolve/gf2_presolve.cuh new file mode 100644 index 000000000..9e1cccc1b --- /dev/null +++ b/cpp/src/mip/presolve/gf2_presolve.cuh @@ -0,0 +1,32 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2025 NVIDIA CORPORATION & AFFILIATES. All rights + * reserved. SPDX-License-Identifier: Apache-2.0 + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#pragma once + +#include +#include +#include +#include + +#include "utils.cuh" + +namespace cuopt::linear_programming::detail { + +template +void gf2_presolve(problem_t& problem); + +} // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip/presolve/trivial_presolve.cuh b/cpp/src/mip/presolve/trivial_presolve.cuh index d7fe14233..081c1297f 100644 --- a/cpp/src/mip/presolve/trivial_presolve.cuh +++ b/cpp/src/mip/presolve/trivial_presolve.cuh @@ -19,6 +19,7 @@ #include #include +#include #include #include #include @@ -356,6 +357,10 @@ void trivial_presolve(problem_t& problem) // The problem has been solved by presolve. Mark its empty status as valid if (problem.n_variables == 0) { problem.empty = true; } problem.check_problem_representation(true); + + gf2_presolve(problem); + problem.recompute_auxilliary_data( + false); // check problem representation later once cstr bounds are computed } } // namespace cuopt::linear_programming::detail From e16e6081c9d71c602a5642a963035e7c622e33e5 Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Sat, 20 Sep 2025 09:06:34 +0000 Subject: [PATCH 2/9] working, b4 moving to papilo presolver --- cpp/src/mip/presolve/gf2_presolve.cu | 291 ++++++++++++++++++++++ cpp/src/mip/presolve/trivial_presolve.cuh | 3 +- 2 files changed, 293 insertions(+), 1 deletion(-) diff --git a/cpp/src/mip/presolve/gf2_presolve.cu b/cpp/src/mip/presolve/gf2_presolve.cu index 62ff83542..8a31a5e82 100644 --- a/cpp/src/mip/presolve/gf2_presolve.cu +++ b/cpp/src/mip/presolve/gf2_presolve.cu @@ -20,6 +20,8 @@ #include #include +#include + namespace cuopt::linear_programming::detail { // Inspired from @@ -32,6 +34,208 @@ namespace cuopt::linear_programming::detail { goto not_valid; \ } +class BinaryMatrix { + private: + std::vector> m_rows; + + public: + /*************************************************************************/ + BinaryMatrix(void) { this->initialize(); } + + /*************************************************************************/ + BinaryMatrix(const int a_NUMBER_OF_ROWS, const int a_NUMBER_OF_COLUMNS) + { + this->setup(a_NUMBER_OF_ROWS, a_NUMBER_OF_COLUMNS); + } + + /*************************************************************************/ + inline void initialize(void) { m_rows.clear(); } + + /*************************************************************************/ + inline void setup(const int a_NUMBER_OF_ROWS, const int a_NUMBER_OF_COLUMNS) + { + m_rows.resize(a_NUMBER_OF_ROWS, std::vector(a_NUMBER_OF_COLUMNS, 0)); + } + + /*************************************************************************/ + inline int number_of_rows(void) const { return m_rows.size(); } + + /*************************************************************************/ + inline int number_of_columns(void) const { return m_rows.front().size(); } + + /*************************************************************************/ + inline std::vector& operator[](const int a_ROW) { return m_rows[a_ROW]; } + + /*************************************************************************/ + inline std::vector const& operator[](const int a_ROW) const { return m_rows[a_ROW]; } + + /*************************************************************************/ + inline void print(void) const + { + const int NUMBER_OF_ROWS = this->number_of_rows(); + const int NUMBER_OF_COLUMNS = this->number_of_columns(); + for (auto i = 0; i < NUMBER_OF_ROWS; i++) { + for (auto j = 0; j < NUMBER_OF_COLUMNS; j++) { + std::cout << m_rows[i][j] << " "; + } + std::cout << std::endl; + } + } + + /*************************************************************************/ + inline std::pair inverse_and_rank(void) const + { + const int SIZE = this->number_of_rows(); + BinaryMatrix A = *this; + BinaryMatrix B(SIZE, SIZE); + int rank = 0; + + for (auto i = 0; i < SIZE; i++) { + B[i][i] = 1; + } + + for (auto j = 0; j < SIZE; j++) { + int row = -1; + for (auto i = rank; i < SIZE; i++) { + if (A[i][j] == 1) { + row = i; + break; + } + } + if (row == -1) { continue; } + + if (row != rank) { + swap(A[row], A[rank]); + swap(B[row], B[rank]); + } + + for (auto i = rank + 1; i < SIZE; i++) { + if (A[i][j] == 1) { + for (auto k = 0; k < SIZE; k++) { + A[i][k] = (A[i][k] + A[rank][k]) & 1; + B[i][k] = (B[i][k] + B[rank][k]) & 1; + } + } + } + rank++; + } + + if (rank == SIZE) { + for (auto j = SIZE - 1; j > 0; j--) { + for (auto i = j - 1; i >= 0; i--) { + if (A[i][j] == 1) { + for (auto k = 0; k < SIZE; k++) { + A[i][k] = (A[j][k] + A[i][k]) & 1; + B[i][k] = (B[j][k] + B[i][k]) & 1; + } + } + } + } + } + + return {B, rank}; + } + + /*************************************************************************/ + inline std::vector dot(const std::vector& a_VECTOR) const + { + std::vector result(m_rows.size(), 0); + const int NUMBER_OF_ROWS = this->number_of_rows(); + const int NUMBER_OF_COLUMNS = this->number_of_columns(); + for (auto i = 0; i < NUMBER_OF_ROWS; i++) { + for (auto j = 0; j < NUMBER_OF_COLUMNS; j++) { + result[i] += m_rows[i][j] * a_VECTOR[j]; + } + result[i] &= 1; + } + return result; + } + + /*************************************************************************/ + inline BinaryMatrix dot(const BinaryMatrix& a_MATRIX) const + { + const int NUMBER_OF_ROWS = this->number_of_rows(); + const int NUMBER_OF_COLUMNS = this->number_of_columns(); + const int RESULT_NUMBER_OF_COLUMNS = a_MATRIX.number_of_columns(); + + BinaryMatrix result(NUMBER_OF_ROWS, RESULT_NUMBER_OF_COLUMNS); + + for (auto i = 0; i < NUMBER_OF_ROWS; i++) { + for (auto j = 0; j < RESULT_NUMBER_OF_COLUMNS; j++) { + for (auto k = 0; k < NUMBER_OF_COLUMNS; k++) { + result[i][j] += m_rows[i][k] * a_MATRIX[k][j]; + } + result[i][j] &= 1; + } + } + return result; + } + + /*************************************************************************/ + inline BinaryMatrix reachability(void) const + { + auto reachability = *this; + const int SIZE = m_rows.size(); + + std::vector> nonzeros(SIZE); + for (auto i = 0; i < SIZE; i++) { + for (auto j = 0; j < SIZE; j++) { + if (reachability[i][j] > 0) { nonzeros[i].insert(j); } + } + } + + for (auto l = 0; l < SIZE; l++) { + bool is_updated = false; + for (auto i = 0; i < SIZE; i++) { + for (auto j = 0; j < SIZE; j++) { + if (reachability[i][j] > 0) { continue; } + + for (auto&& k : nonzeros[i]) { + if (reachability[k][j]) { + reachability[i][j] = 1; + is_updated = true; + break; + } + } + } + } + if (!is_updated) { break; } + } + + return reachability; + } + + /*************************************************************************/ + inline static BinaryMatrix identity(const int a_SIZE) + { + auto identity = BinaryMatrix(a_SIZE, a_SIZE); + for (auto i = 0; i < a_SIZE; i++) { + identity[i][i] = 1; + } + return identity; + } +}; + +struct gf2_constraint_t { + size_t cstr_idx; + std::vector> bin_vars; + std::pair key_var; + size_t rhs; // 0 or 1 + + gf2_constraint_t() = default; + gf2_constraint_t(size_t cstr_idx, + std::vector> bin_vars, + std::pair key_var, + size_t rhs) + : cstr_idx(cstr_idx), bin_vars(std::move(bin_vars)), key_var(key_var), rhs(rhs) + { + } + gf2_constraint_t(const gf2_constraint_t& other) = default; + gf2_constraint_t(gf2_constraint_t&& other) noexcept = default; + gf2_constraint_t& operator=(const gf2_constraint_t& other) = default; + gf2_constraint_t& operator=(gf2_constraint_t&& other) noexcept = default; +}; + // TODO: implement fully on GPU + cuDSS? template void gf2_presolve(problem_t& problem) @@ -45,6 +249,12 @@ void gf2_presolve(problem_t& problem) auto h_offsets = cuopt::host_copy(problem.offsets); auto h_variable_types = cuopt::host_copy(problem.variable_types); + // maps problem:var_idx -> gf2matrix:bin/key_idx + std::unordered_map gf2_bin_vars; + std::unordered_map gf2_key_vars; + + std::vector gf2_constraints; + for (i_t cstr_idx = 0; cstr_idx < problem.n_constraints; ++cstr_idx) { i_t key_var_idx = -1; f_t key_var_coeff = 0.0; @@ -53,6 +263,8 @@ void gf2_presolve(problem_t& problem) f_t min_activity_without_key = 0.0; f_t max_activity_without_key = 0.0; + std::vector> constraint_bin_vars; + //__asm__ __volatile__("int3"); f_t rhs = round(h_constraint_lower_bounds[cstr_idx]); @@ -89,11 +301,16 @@ void gf2_presolve(problem_t& problem) if (key_var_idx != -1) NOT_GF2("multiple key variables"); key_var_idx = var_idx; key_var_coeff = coeff; + + gf2_key_vars.insert({var_idx, gf2_key_vars.size()}); } else { min_activity_without_key += min_act_of_var( coeff, get_lower(h_variable_bounds[var_idx]), get_upper(h_variable_bounds[var_idx])); max_activity_without_key += max_act_of_var( coeff, get_lower(h_variable_bounds[var_idx]), get_upper(h_variable_bounds[var_idx])); + + constraint_bin_vars.push_back({var_idx, (double)coeff}); + gf2_bin_vars.insert({var_idx, gf2_bin_vars.size()}); } } @@ -121,12 +338,86 @@ void gf2_presolve(problem_t& problem) max_activity_without_key); printf("YES: Cons %d is gf2\n", cstr_idx); + gf2_constraints.emplace_back((size_t)cstr_idx, + std::move(constraint_bin_vars), + std::pair{key_var_idx, (double)key_var_coeff}, + ((size_t)rhs) % 2); continue; not_valid: continue; } + if (gf2_key_vars.size() != gf2_constraints.size()) { + printf( + "invalid key variable count: %d/%d", (int)gf2_key_vars.size(), (int)gf2_constraints.size()); + exit(0); + } + if (gf2_bin_vars.size() != gf2_constraints.size()) { + printf("invalid binary variable count: %d/%d", + (int)gf2_bin_vars.size(), + (int)gf2_constraints.size()); + exit(0); + } + + printf("Valid GF2 structure!\n"); + + // maps gf2:bin/key_idx -> problem:var_idx + std::unordered_map gf2_bin_vars_invmap; + std::unordered_map gf2_key_vars_invmap; + for (const auto& [var_idx, gf2_idx] : gf2_bin_vars) { + cuopt_assert(gf2_bin_vars_invmap.count(gf2_idx) == 0, "not a bijection??"); + gf2_bin_vars_invmap.insert({gf2_idx, var_idx}); + } + for (const auto& [var_idx, gf2_idx] : gf2_key_vars) { + cuopt_assert(gf2_key_vars_invmap.count(gf2_idx) == 0, "not a bijection??"); + gf2_key_vars_invmap.insert({gf2_idx, var_idx}); + } + + BinaryMatrix A(gf2_constraints.size(), gf2_constraints.size()); + std::vector b(gf2_constraints.size()); + for (const auto& cons : gf2_constraints) { + for (auto [bin_var, _] : cons.bin_vars) { + cuopt_assert(gf2_bin_vars.count(bin_var) == 1, ""); + A[cons.cstr_idx][gf2_bin_vars[bin_var]] = 1; + } + b[cons.cstr_idx] = cons.rhs; + } + + auto [inverse, rank] = A.inverse_and_rank(); + if (rank != (int)gf2_constraints.size()) { + printf("non invertible\n"); + exit(0); + } + + auto solution = inverse.dot(b); + + std::unordered_map fixings; + + for (size_t sol_idx = 0; sol_idx < gf2_constraints.size(); ++sol_idx) { + fixings[gf2_bin_vars_invmap[sol_idx]] = solution[sol_idx]; + } + + // compute fixings for the key variables by solving the corresponding constraint + for (const auto& cons : gf2_constraints) { + auto [key_var_idx, key_var_coeff] = cons.key_var; + f_t rhs = h_constraint_lower_bounds[cons.cstr_idx]; // eq constraint + f_t lhs = 0.0; + for (auto [bin_var, coeff] : cons.bin_vars) { + lhs += fixings[bin_var] * coeff; + } + lhs -= rhs; + + cuopt_assert(fixings.count(key_var_idx) == 0, "key var unexpectedly already fixed"); + cuopt_assert(key_var_coeff != 0, "key var coeff is 0"); + fixings[key_var_idx] = round(-lhs / key_var_coeff); + } + + printf("Fixings:\n"); + for (const auto& [var_idx, fixing] : fixings) { + if (fixing != 0) printf("%s %d\n", problem.var_names[var_idx].c_str(), (int)round(fixing)); + } + exit(0); } diff --git a/cpp/src/mip/presolve/trivial_presolve.cuh b/cpp/src/mip/presolve/trivial_presolve.cuh index 081c1297f..f2cb61050 100644 --- a/cpp/src/mip/presolve/trivial_presolve.cuh +++ b/cpp/src/mip/presolve/trivial_presolve.cuh @@ -346,6 +346,8 @@ void test_reverse_matches(const problem_t& pb) template void trivial_presolve(problem_t& problem) { + gf2_presolve(problem); + cuopt_expects(problem.preprocess_called, error_type_t::RuntimeError, "preprocess_problem should be called before running the solver"); @@ -358,7 +360,6 @@ void trivial_presolve(problem_t& problem) if (problem.n_variables == 0) { problem.empty = true; } problem.check_problem_representation(true); - gf2_presolve(problem); problem.recompute_auxilliary_data( false); // check problem representation later once cstr bounds are computed } From 9d9501461b97cd3845e27e30a850839ef5ea2e19 Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Sat, 20 Sep 2025 10:31:56 +0000 Subject: [PATCH 3/9] added as a papilo presolver --- .../linear_programming/cuopt/run_mip.cpp | 2 +- cpp/src/mip/CMakeLists.txt | 4 +- cpp/src/mip/presolve/gf2_presolve.cpp | 216 +++++++++ cpp/src/mip/presolve/gf2_presolve.cu | 437 ------------------ cpp/src/mip/presolve/gf2_presolve.cuh | 32 -- cpp/src/mip/presolve/gf2_presolve.hpp | 176 +++++++ cpp/src/mip/presolve/third_party_presolve.cpp | 4 + cpp/src/mip/presolve/trivial_presolve.cuh | 3 - 8 files changed, 398 insertions(+), 476 deletions(-) create mode 100644 cpp/src/mip/presolve/gf2_presolve.cpp delete mode 100644 cpp/src/mip/presolve/gf2_presolve.cu delete mode 100644 cpp/src/mip/presolve/gf2_presolve.cuh create mode 100644 cpp/src/mip/presolve/gf2_presolve.hpp diff --git a/benchmarks/linear_programming/cuopt/run_mip.cpp b/benchmarks/linear_programming/cuopt/run_mip.cpp index 904274736..8ce9fe5e7 100644 --- a/benchmarks/linear_programming/cuopt/run_mip.cpp +++ b/benchmarks/linear_programming/cuopt/run_mip.cpp @@ -210,7 +210,7 @@ int run_single_file(std::string file_path, settings.log_to_console = log_to_console; settings.tolerances.relative_tolerance = 1e-12; settings.tolerances.absolute_tolerance = 1e-6; - settings.presolve = false; + settings.presolve = true; settings.mip_scaling = false; cuopt::linear_programming::benchmark_info_t benchmark_info; settings.benchmark_info_ptr = &benchmark_info; diff --git a/cpp/src/mip/CMakeLists.txt b/cpp/src/mip/CMakeLists.txt index 02ce5944f..e694a9f6f 100644 --- a/cpp/src/mip/CMakeLists.txt +++ b/cpp/src/mip/CMakeLists.txt @@ -21,6 +21,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/gf2_presolve.cpp ${CMAKE_CURRENT_SOURCE_DIR}/solution/solution.cu ) @@ -45,13 +46,10 @@ set(MIP_NON_LP_FILES ${CMAKE_CURRENT_SOURCE_DIR}/presolve/multi_probe.cu ${CMAKE_CURRENT_SOURCE_DIR}/presolve/probing_cache.cu ${CMAKE_CURRENT_SOURCE_DIR}/presolve/trivial_presolve.cu - ${CMAKE_CURRENT_SOURCE_DIR}/presolve/gf2_presolve.cu ${CMAKE_CURRENT_SOURCE_DIR}/feasibility_jump/feasibility_jump.cu ${CMAKE_CURRENT_SOURCE_DIR}/feasibility_jump/feasibility_jump_kernels.cu ) -set_source_files_properties(${CMAKE_CURRENT_SOURCE_DIR}/presolve/gf2_presolve.cu DIRECTORY ${CMAKE_SOURCE_DIR} PROPERTIES COMPILE_OPTIONS "-g;-O0") - # Choose which files to include based on build mode if(BUILD_LP_ONLY) set(MIP_SRC_FILES ${MIP_LP_NECESSARY_FILES}) diff --git a/cpp/src/mip/presolve/gf2_presolve.cpp b/cpp/src/mip/presolve/gf2_presolve.cpp new file mode 100644 index 000000000..5fbc81fc1 --- /dev/null +++ b/cpp/src/mip/presolve/gf2_presolve.cpp @@ -0,0 +1,216 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights + * reserved. SPDX-License-Identifier: Apache-2.0 + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "gf2_presolve.hpp" +#include +#include +#include + +#define NOT_GF2(reason, ...) \ + { \ + printf("NO : Cons %d is not gf2: " reason "\n", cstr_idx, ##__VA_ARGS__); \ + goto not_valid; \ + } + +namespace cuopt::linear_programming::detail { + +template +papilo::PresolveStatus GF2Presolve::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 int num_rows = constraint_matrix.getNRows(); + + std::unordered_map gf2_bin_vars; + std::unordered_map gf2_key_vars; + std::vector gf2_constraints; + + const f_t integrality_tolerance = num.getFeasTol(); + + for (int cstr_idx = 0; cstr_idx < num_rows; ++cstr_idx) { + int key_var_idx = -1; + f_t key_var_coeff = 0.0; + + std::vector> constraint_bin_vars; + + // Check constraint coefficients + auto row_coeff = constraint_matrix.getRowCoefficients(cstr_idx); + const int* row_indices = row_coeff.getIndices(); + const f_t* row_values = row_coeff.getValues(); + const int row_length = row_coeff.getLength(); + f_t rhs = std::round(lhs_values[cstr_idx]); + + // Check if this is an equality constraint + if (!num.isEq(lhs_values[cstr_idx], rhs_values[cstr_idx])) + NOT_GF2("not eq", lhs_values[cstr_idx], rhs_values[cstr_idx]); + if (!std::isfinite(lhs_values[cstr_idx])) NOT_GF2("not finite", lhs_values[cstr_idx]); + if (!is_integer(lhs_values[cstr_idx], integrality_tolerance)) + NOT_GF2("not integer", lhs_values[cstr_idx]); + + // Only accept 0, 1, -1 as rhs + if (rhs != 0.0 && rhs != 1.0 && rhs != -1.0) NOT_GF2("not 0, 1, -1", rhs); + + for (int j = 0; j < row_length; ++j) { + if (!is_integer(row_values[j], integrality_tolerance)) { + NOT_GF2("coeff not integer", row_values[j]); + } + + int var_idx = row_indices[j]; + f_t coeff = std::round(row_values[j]); + + // Check if variable is integer + if (!col_flags[var_idx].test(papilo::ColFlag::kIntegral)) { + NOT_GF2("not integral", var_idx); + } + + bool is_binary = col_flags[var_idx].test(papilo::ColFlag::kLbInf) ? false + : col_flags[var_idx].test(papilo::ColFlag::kUbInf) + ? false + : (lower_bounds[var_idx] == 0.0 && upper_bounds[var_idx] == 1.0); + + // Check coefficient constraints + if (is_binary && (std::abs(coeff) != 1.0 && std::abs(coeff) != 2.0)) { + NOT_GF2("not binary", var_idx); + } + if (!is_binary && (std::abs(coeff) != 2.0)) { NOT_GF2("not binary", var_idx); } + + // Key variable (coefficient of 2) + if (std::abs(coeff) == 2.0) { + 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 + f_t var_lb = col_flags[var_idx].test(papilo::ColFlag::kLbInf) + ? -std::numeric_limits::infinity() + : lower_bounds[var_idx]; + f_t var_ub = col_flags[var_idx].test(papilo::ColFlag::kUbInf) + ? std::numeric_limits::infinity() + : upper_bounds[var_idx]; + + 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"); + + gf2_constraints.emplace_back((size_t)cstr_idx, + std::move(constraint_bin_vars), + std::pair{key_var_idx, key_var_coeff}, + ((size_t)rhs) % 2); + continue; + not_valid: + continue; + } + + // If no GF2 constraints found, return unchanged + if (gf2_constraints.empty()) { + printf("No GF2 constraints found\n"); + exit(0); + return papilo::PresolveStatus::kUnchanged; + } + + // Validate structure + if (gf2_key_vars.size() != gf2_constraints.size() || + gf2_bin_vars.size() != gf2_constraints.size()) { + printf("GF2 constraints: %d, %d, %d\n", + gf2_key_vars.size(), + gf2_bin_vars.size(), + gf2_constraints.size()); + exit(0); + return papilo::PresolveStatus::kUnchanged; + } + + // Create inverse mappings + std::unordered_map gf2_bin_vars_invmap; + for (const auto& [var_idx, gf2_idx] : gf2_bin_vars) { + gf2_bin_vars_invmap.insert({gf2_idx, var_idx}); + } + + // Build binary matrix + BinaryMatrix A(gf2_constraints.size(), gf2_constraints.size()); + std::vector b(gf2_constraints.size()); + for (const auto& cons : gf2_constraints) { + for (auto [bin_var, _] : cons.bin_vars) { + A[cons.cstr_idx][gf2_bin_vars[bin_var]] = 1; + } + b[cons.cstr_idx] = cons.rhs; + } + + auto [inverse, rank] = A.inverse_and_rank(); + if (rank != (int)gf2_constraints.size()) { + printf("Non-invertible, rank: %d/%d\n", rank, gf2_constraints.size()); + exit(0); + return papilo::PresolveStatus::kUnchanged; // Non-invertible + } + + auto solution = inverse.dot(b); + + 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]; + } + + // Compute fixings for key variables + for (const auto& cons : gf2_constraints) { + auto [key_var_idx, key_var_coeff] = cons.key_var; + f_t constraint_rhs = lhs_values[cons.cstr_idx]; // equality constraint + f_t lhs = 0.0; + for (auto [bin_var, coeff] : cons.bin_vars) { + lhs += fixings[bin_var] * coeff; + } + lhs -= constraint_rhs; + + fixings[key_var_idx] = std::round(-lhs / key_var_coeff); + } + + papilo::PresolveStatus status = papilo::PresolveStatus::kUnchanged; + for (const auto& [var_idx, fixing] : fixings) { + if (num.isZero(fixing)) { + reductions.fixCol(var_idx, 0); + } else { + reductions.fixCol(var_idx, fixing); + } + status = papilo::PresolveStatus::kReduced; + } + + printf("GF2 presolved!!\n"); + + return status; +} + +// Explicit template instantiations +template class GF2Presolve; +template class GF2Presolve; + +} // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip/presolve/gf2_presolve.cu b/cpp/src/mip/presolve/gf2_presolve.cu deleted file mode 100644 index 8a31a5e82..000000000 --- a/cpp/src/mip/presolve/gf2_presolve.cu +++ /dev/null @@ -1,437 +0,0 @@ -/* - * SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights - * reserved. SPDX-License-Identifier: Apache-2.0 - * - * Licensed under the Apache License, Version 2.0 (the "License"); - * you may not use this file except in compliance with the License. - * You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -#include - -#include -#include - -#include - -namespace cuopt::linear_programming::detail { - -// Inspired from -// https://github.com/snowberryfield/printemps/blob/0dc1f9d2b78a667451c11d8ab156861ce54e9fb7/printemps/model_component/constraint.h -// (MIT License) - -#define NOT_GF2(reason, ...) \ - { \ - printf("NO : Cons %d is not gf2: " reason "\n", cstr_idx, ##__VA_ARGS__); \ - goto not_valid; \ - } - -class BinaryMatrix { - private: - std::vector> m_rows; - - public: - /*************************************************************************/ - BinaryMatrix(void) { this->initialize(); } - - /*************************************************************************/ - BinaryMatrix(const int a_NUMBER_OF_ROWS, const int a_NUMBER_OF_COLUMNS) - { - this->setup(a_NUMBER_OF_ROWS, a_NUMBER_OF_COLUMNS); - } - - /*************************************************************************/ - inline void initialize(void) { m_rows.clear(); } - - /*************************************************************************/ - inline void setup(const int a_NUMBER_OF_ROWS, const int a_NUMBER_OF_COLUMNS) - { - m_rows.resize(a_NUMBER_OF_ROWS, std::vector(a_NUMBER_OF_COLUMNS, 0)); - } - - /*************************************************************************/ - inline int number_of_rows(void) const { return m_rows.size(); } - - /*************************************************************************/ - inline int number_of_columns(void) const { return m_rows.front().size(); } - - /*************************************************************************/ - inline std::vector& operator[](const int a_ROW) { return m_rows[a_ROW]; } - - /*************************************************************************/ - inline std::vector const& operator[](const int a_ROW) const { return m_rows[a_ROW]; } - - /*************************************************************************/ - inline void print(void) const - { - const int NUMBER_OF_ROWS = this->number_of_rows(); - const int NUMBER_OF_COLUMNS = this->number_of_columns(); - for (auto i = 0; i < NUMBER_OF_ROWS; i++) { - for (auto j = 0; j < NUMBER_OF_COLUMNS; j++) { - std::cout << m_rows[i][j] << " "; - } - std::cout << std::endl; - } - } - - /*************************************************************************/ - inline std::pair inverse_and_rank(void) const - { - const int SIZE = this->number_of_rows(); - BinaryMatrix A = *this; - BinaryMatrix B(SIZE, SIZE); - int rank = 0; - - for (auto i = 0; i < SIZE; i++) { - B[i][i] = 1; - } - - for (auto j = 0; j < SIZE; j++) { - int row = -1; - for (auto i = rank; i < SIZE; i++) { - if (A[i][j] == 1) { - row = i; - break; - } - } - if (row == -1) { continue; } - - if (row != rank) { - swap(A[row], A[rank]); - swap(B[row], B[rank]); - } - - for (auto i = rank + 1; i < SIZE; i++) { - if (A[i][j] == 1) { - for (auto k = 0; k < SIZE; k++) { - A[i][k] = (A[i][k] + A[rank][k]) & 1; - B[i][k] = (B[i][k] + B[rank][k]) & 1; - } - } - } - rank++; - } - - if (rank == SIZE) { - for (auto j = SIZE - 1; j > 0; j--) { - for (auto i = j - 1; i >= 0; i--) { - if (A[i][j] == 1) { - for (auto k = 0; k < SIZE; k++) { - A[i][k] = (A[j][k] + A[i][k]) & 1; - B[i][k] = (B[j][k] + B[i][k]) & 1; - } - } - } - } - } - - return {B, rank}; - } - - /*************************************************************************/ - inline std::vector dot(const std::vector& a_VECTOR) const - { - std::vector result(m_rows.size(), 0); - const int NUMBER_OF_ROWS = this->number_of_rows(); - const int NUMBER_OF_COLUMNS = this->number_of_columns(); - for (auto i = 0; i < NUMBER_OF_ROWS; i++) { - for (auto j = 0; j < NUMBER_OF_COLUMNS; j++) { - result[i] += m_rows[i][j] * a_VECTOR[j]; - } - result[i] &= 1; - } - return result; - } - - /*************************************************************************/ - inline BinaryMatrix dot(const BinaryMatrix& a_MATRIX) const - { - const int NUMBER_OF_ROWS = this->number_of_rows(); - const int NUMBER_OF_COLUMNS = this->number_of_columns(); - const int RESULT_NUMBER_OF_COLUMNS = a_MATRIX.number_of_columns(); - - BinaryMatrix result(NUMBER_OF_ROWS, RESULT_NUMBER_OF_COLUMNS); - - for (auto i = 0; i < NUMBER_OF_ROWS; i++) { - for (auto j = 0; j < RESULT_NUMBER_OF_COLUMNS; j++) { - for (auto k = 0; k < NUMBER_OF_COLUMNS; k++) { - result[i][j] += m_rows[i][k] * a_MATRIX[k][j]; - } - result[i][j] &= 1; - } - } - return result; - } - - /*************************************************************************/ - inline BinaryMatrix reachability(void) const - { - auto reachability = *this; - const int SIZE = m_rows.size(); - - std::vector> nonzeros(SIZE); - for (auto i = 0; i < SIZE; i++) { - for (auto j = 0; j < SIZE; j++) { - if (reachability[i][j] > 0) { nonzeros[i].insert(j); } - } - } - - for (auto l = 0; l < SIZE; l++) { - bool is_updated = false; - for (auto i = 0; i < SIZE; i++) { - for (auto j = 0; j < SIZE; j++) { - if (reachability[i][j] > 0) { continue; } - - for (auto&& k : nonzeros[i]) { - if (reachability[k][j]) { - reachability[i][j] = 1; - is_updated = true; - break; - } - } - } - } - if (!is_updated) { break; } - } - - return reachability; - } - - /*************************************************************************/ - inline static BinaryMatrix identity(const int a_SIZE) - { - auto identity = BinaryMatrix(a_SIZE, a_SIZE); - for (auto i = 0; i < a_SIZE; i++) { - identity[i][i] = 1; - } - return identity; - } -}; - -struct gf2_constraint_t { - size_t cstr_idx; - std::vector> bin_vars; - std::pair key_var; - size_t rhs; // 0 or 1 - - gf2_constraint_t() = default; - gf2_constraint_t(size_t cstr_idx, - std::vector> bin_vars, - std::pair key_var, - size_t rhs) - : cstr_idx(cstr_idx), bin_vars(std::move(bin_vars)), key_var(key_var), rhs(rhs) - { - } - gf2_constraint_t(const gf2_constraint_t& other) = default; - gf2_constraint_t(gf2_constraint_t&& other) noexcept = default; - gf2_constraint_t& operator=(const gf2_constraint_t& other) = default; - gf2_constraint_t& operator=(gf2_constraint_t&& other) noexcept = default; -}; - -// TODO: implement fully on GPU + cuDSS? -template -void gf2_presolve(problem_t& problem) -{ - auto h_constraint_lower_bounds = cuopt::host_copy(problem.constraint_lower_bounds); - auto h_constraint_upper_bounds = cuopt::host_copy(problem.constraint_upper_bounds); - auto h_coefficients = cuopt::host_copy(problem.coefficients); - auto h_variables = cuopt::host_copy(problem.variables); - auto h_is_binary_variable = cuopt::host_copy(problem.is_binary_variable); - auto h_variable_bounds = cuopt::host_copy(problem.variable_bounds); - auto h_offsets = cuopt::host_copy(problem.offsets); - auto h_variable_types = cuopt::host_copy(problem.variable_types); - - // maps problem:var_idx -> gf2matrix:bin/key_idx - std::unordered_map gf2_bin_vars; - std::unordered_map gf2_key_vars; - - std::vector gf2_constraints; - - for (i_t cstr_idx = 0; cstr_idx < problem.n_constraints; ++cstr_idx) { - i_t key_var_idx = -1; - f_t key_var_coeff = 0.0; - f_t key_var_lb = 0.0; - f_t key_var_ub = 0.0; - f_t min_activity_without_key = 0.0; - f_t max_activity_without_key = 0.0; - - std::vector> constraint_bin_vars; - - //__asm__ __volatile__("int3"); - - f_t rhs = round(h_constraint_lower_bounds[cstr_idx]); - - // needs to be an equality constraint - if (h_constraint_lower_bounds[cstr_idx] != h_constraint_upper_bounds[cstr_idx]) - NOT_GF2("not an equality constraint"); - if (!isfinite(h_constraint_lower_bounds[cstr_idx])) NOT_GF2("not finite"); - if (!is_integer(h_constraint_lower_bounds[cstr_idx], problem.tolerances.integrality_tolerance)) - NOT_GF2("rhs not integer %f", h_constraint_lower_bounds[cstr_idx]); - - // only accept 0, 1, -1 as rhs - if (rhs != 0.0 && rhs != 1.0 && rhs != -1.0) NOT_GF2("invalid rhs %f", rhs); - - // check if nnzs match the pattern of a gf2 constraint - for (i_t j = h_offsets[cstr_idx]; j < h_offsets[cstr_idx + 1]; ++j) { - if (!is_integer(h_coefficients[j], problem.tolerances.integrality_tolerance)) - NOT_GF2("coeff not integer (%d,%f)", j, h_coefficients[j]); - i_t var_idx = h_variables[j]; - f_t coeff = round(h_coefficients[j]); - if (h_variable_types[var_idx] != var_t::INTEGER) - NOT_GF2("var not integer (%d,%f)", var_idx, h_coefficients[j]); - bool is_binary = h_is_binary_variable[var_idx]; - - // only the key variable can have a coefficient of 2 - if (is_binary && (abs(coeff) != 1.0 && abs(coeff) != 2.0)) - NOT_GF2("invalid coef bin (%d,%f)", var_idx, coeff); - if (!is_binary && (abs(coeff) != 2.0)) - NOT_GF2("invalid coef non-bin (%d,%f)", var_idx, coeff); - - // key var - if (abs(coeff) == 2.0) { - // can only be one - if (key_var_idx != -1) NOT_GF2("multiple key variables"); - key_var_idx = var_idx; - key_var_coeff = coeff; - - gf2_key_vars.insert({var_idx, gf2_key_vars.size()}); - } else { - min_activity_without_key += min_act_of_var( - coeff, get_lower(h_variable_bounds[var_idx]), get_upper(h_variable_bounds[var_idx])); - max_activity_without_key += max_act_of_var( - coeff, get_lower(h_variable_bounds[var_idx]), get_upper(h_variable_bounds[var_idx])); - - constraint_bin_vars.push_back({var_idx, (double)coeff}); - gf2_bin_vars.insert({var_idx, gf2_bin_vars.size()}); - } - } - - // no key var found - if (key_var_idx == -1) NOT_GF2("no key variable"); - - if (key_var_coeff > 0) { - std::swap(min_activity_without_key, max_activity_without_key); - min_activity_without_key *= -1; - max_activity_without_key *= -1; - } - cuopt_assert(min_activity_without_key <= max_activity_without_key, "invalid activities?"); - - key_var_lb = get_lower(h_variable_bounds[key_var_idx]); - key_var_ub = get_upper(h_variable_bounds[key_var_idx]); - if (isfinite(key_var_lb) && key_var_lb > ceil(min_activity_without_key / 2.0)) - NOT_GF2("invalid key variable bounds lower ((%f,%f), %f)", - key_var_lb, - key_var_ub, - min_activity_without_key); - if (isfinite(key_var_ub) && key_var_ub < floor(max_activity_without_key / 2.0)) - NOT_GF2("invalid key variable bounds upper ((%f,%f), %f)", - key_var_lb, - key_var_ub, - max_activity_without_key); - - printf("YES: Cons %d is gf2\n", cstr_idx); - gf2_constraints.emplace_back((size_t)cstr_idx, - std::move(constraint_bin_vars), - std::pair{key_var_idx, (double)key_var_coeff}, - ((size_t)rhs) % 2); - continue; - - not_valid: - continue; - } - - if (gf2_key_vars.size() != gf2_constraints.size()) { - printf( - "invalid key variable count: %d/%d", (int)gf2_key_vars.size(), (int)gf2_constraints.size()); - exit(0); - } - if (gf2_bin_vars.size() != gf2_constraints.size()) { - printf("invalid binary variable count: %d/%d", - (int)gf2_bin_vars.size(), - (int)gf2_constraints.size()); - exit(0); - } - - printf("Valid GF2 structure!\n"); - - // maps gf2:bin/key_idx -> problem:var_idx - std::unordered_map gf2_bin_vars_invmap; - std::unordered_map gf2_key_vars_invmap; - for (const auto& [var_idx, gf2_idx] : gf2_bin_vars) { - cuopt_assert(gf2_bin_vars_invmap.count(gf2_idx) == 0, "not a bijection??"); - gf2_bin_vars_invmap.insert({gf2_idx, var_idx}); - } - for (const auto& [var_idx, gf2_idx] : gf2_key_vars) { - cuopt_assert(gf2_key_vars_invmap.count(gf2_idx) == 0, "not a bijection??"); - gf2_key_vars_invmap.insert({gf2_idx, var_idx}); - } - - BinaryMatrix A(gf2_constraints.size(), gf2_constraints.size()); - std::vector b(gf2_constraints.size()); - for (const auto& cons : gf2_constraints) { - for (auto [bin_var, _] : cons.bin_vars) { - cuopt_assert(gf2_bin_vars.count(bin_var) == 1, ""); - A[cons.cstr_idx][gf2_bin_vars[bin_var]] = 1; - } - b[cons.cstr_idx] = cons.rhs; - } - - auto [inverse, rank] = A.inverse_and_rank(); - if (rank != (int)gf2_constraints.size()) { - printf("non invertible\n"); - exit(0); - } - - auto solution = inverse.dot(b); - - std::unordered_map fixings; - - for (size_t sol_idx = 0; sol_idx < gf2_constraints.size(); ++sol_idx) { - fixings[gf2_bin_vars_invmap[sol_idx]] = solution[sol_idx]; - } - - // compute fixings for the key variables by solving the corresponding constraint - for (const auto& cons : gf2_constraints) { - auto [key_var_idx, key_var_coeff] = cons.key_var; - f_t rhs = h_constraint_lower_bounds[cons.cstr_idx]; // eq constraint - f_t lhs = 0.0; - for (auto [bin_var, coeff] : cons.bin_vars) { - lhs += fixings[bin_var] * coeff; - } - lhs -= rhs; - - cuopt_assert(fixings.count(key_var_idx) == 0, "key var unexpectedly already fixed"); - cuopt_assert(key_var_coeff != 0, "key var coeff is 0"); - fixings[key_var_idx] = round(-lhs / key_var_coeff); - } - - printf("Fixings:\n"); - for (const auto& [var_idx, fixing] : fixings) { - if (fixing != 0) printf("%s %d\n", problem.var_names[var_idx].c_str(), (int)round(fixing)); - } - - exit(0); -} - -#define INSTANTIATE(F_TYPE) \ - template void gf2_presolve(problem_t & problem); - -#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/presolve/gf2_presolve.cuh b/cpp/src/mip/presolve/gf2_presolve.cuh deleted file mode 100644 index 9e1cccc1b..000000000 --- a/cpp/src/mip/presolve/gf2_presolve.cuh +++ /dev/null @@ -1,32 +0,0 @@ -/* - * SPDX-FileCopyrightText: Copyright (c) 2025 NVIDIA CORPORATION & AFFILIATES. All rights - * reserved. SPDX-License-Identifier: Apache-2.0 - * - * Licensed under the Apache License, Version 2.0 (the "License"); - * you may not use this file except in compliance with the License. - * You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -#pragma once - -#include -#include -#include -#include - -#include "utils.cuh" - -namespace cuopt::linear_programming::detail { - -template -void gf2_presolve(problem_t& problem); - -} // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip/presolve/gf2_presolve.hpp b/cpp/src/mip/presolve/gf2_presolve.hpp new file mode 100644 index 000000000..6d389b377 --- /dev/null +++ b/cpp/src/mip/presolve/gf2_presolve.hpp @@ -0,0 +1,176 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights + * reserved. SPDX-License-Identifier: Apache-2.0 + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#pragma once + +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wstringop-overflow" // ignore boost error for pip wheel build +#include +#include +#include +#include +#pragma GCC diagnostic pop + +namespace cuopt::linear_programming::detail { + +template +class GF2Presolve : public papilo::PresolveMethod { + public: + GF2Presolve() : papilo::PresolveMethod() + { + this->setName("gf2presolve"); + 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: + struct gf2_constraint_t { + size_t cstr_idx; + std::vector> bin_vars; + std::pair key_var; + size_t rhs; // 0 or 1 + + gf2_constraint_t() = default; + gf2_constraint_t(size_t cstr_idx, + std::vector> bin_vars, + std::pair key_var, + size_t rhs) + : cstr_idx(cstr_idx), bin_vars(std::move(bin_vars)), key_var(key_var), rhs(rhs) + { + } + gf2_constraint_t(const gf2_constraint_t& other) = default; + gf2_constraint_t(gf2_constraint_t&& other) noexcept = default; + gf2_constraint_t& operator=(const gf2_constraint_t& other) = default; + gf2_constraint_t& operator=(gf2_constraint_t&& other) noexcept = default; + }; + + class BinaryMatrix { + private: + std::vector> m_rows; + + public: + BinaryMatrix(void) { this->initialize(); } + BinaryMatrix(const int a_NUMBER_OF_ROWS, const int a_NUMBER_OF_COLUMNS) + { + this->setup(a_NUMBER_OF_ROWS, a_NUMBER_OF_COLUMNS); + } + + inline void initialize(void) { m_rows.clear(); } + inline void setup(const int a_NUMBER_OF_ROWS, const int a_NUMBER_OF_COLUMNS) + { + m_rows.resize(a_NUMBER_OF_ROWS, std::vector(a_NUMBER_OF_COLUMNS, 0)); + } + + inline int number_of_rows(void) const { return m_rows.size(); } + inline int number_of_columns(void) const { return m_rows.front().size(); } + inline std::vector& operator[](const int a_ROW) { return m_rows[a_ROW]; } + inline std::vector const& operator[](const int a_ROW) const { return m_rows[a_ROW]; } + + inline std::pair inverse_and_rank(void) const + { + const int SIZE = this->number_of_rows(); + BinaryMatrix A = *this; + BinaryMatrix B(SIZE, SIZE); + int rank = 0; + + for (auto i = 0; i < SIZE; i++) { + B[i][i] = 1; + } + + for (auto j = 0; j < SIZE; j++) { + int row = -1; + for (auto i = rank; i < SIZE; i++) { + if (A[i][j] == 1) { + row = i; + break; + } + } + if (row == -1) { continue; } + + if (row != rank) { + swap(A[row], A[rank]); + swap(B[row], B[rank]); + } + + for (auto i = rank + 1; i < SIZE; i++) { + if (A[i][j] == 1) { + for (auto k = 0; k < SIZE; k++) { + A[i][k] = (A[i][k] + A[rank][k]) & 1; + B[i][k] = (B[i][k] + B[rank][k]) & 1; + } + } + } + rank++; + } + + if (rank == SIZE) { + for (auto j = SIZE - 1; j > 0; j--) { + for (auto i = j - 1; i >= 0; i--) { + if (A[i][j] == 1) { + for (auto k = 0; k < SIZE; k++) { + A[i][k] = (A[j][k] + A[i][k]) & 1; + B[i][k] = (B[j][k] + B[i][k]) & 1; + } + } + } + } + } + + return {B, rank}; + } + + inline std::vector dot(const std::vector& a_VECTOR) const + { + std::vector result(m_rows.size(), 0); + const int NUMBER_OF_ROWS = this->number_of_rows(); + const int NUMBER_OF_COLUMNS = this->number_of_columns(); + for (auto i = 0; i < NUMBER_OF_ROWS; i++) { + for (auto j = 0; j < NUMBER_OF_COLUMNS; j++) { + result[i] += m_rows[i][j] * a_VECTOR[j]; + } + result[i] &= 1; + } + return result; + } + }; + + bool is_integer(f_t value, f_t tolerance) const + { + return std::abs(value - std::round(value)) <= tolerance; + } + + f_t min_act_of_var(f_t coeff, f_t lb, f_t ub) const + { + if (coeff >= 0) return coeff * lb; + return coeff * ub; + } + + f_t max_act_of_var(f_t coeff, f_t lb, f_t ub) const + { + if (coeff >= 0) return coeff * ub; + return coeff * lb; + } +}; + +} // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip/presolve/third_party_presolve.cpp b/cpp/src/mip/presolve/third_party_presolve.cpp index 211db4605..46c2ea4bc 100644 --- a/cpp/src/mip/presolve/third_party_presolve.cpp +++ b/cpp/src/mip/presolve/third_party_presolve.cpp @@ -18,6 +18,7 @@ #include #include #include +#include #include #include @@ -308,6 +309,9 @@ void set_presolve_methods(papilo::Presolve& presolver, problem_category_t c { using uptr = std::unique_ptr>; + // cuopt custom presolvers + presolver.addPresolveMethod(uptr(new cuopt::linear_programming::detail::GF2Presolve())); + // fast presolvers presolver.addPresolveMethod(uptr(new papilo::SingletonCols())); presolver.addPresolveMethod(uptr(new papilo::CoefficientStrengthening())); diff --git a/cpp/src/mip/presolve/trivial_presolve.cuh b/cpp/src/mip/presolve/trivial_presolve.cuh index f2cb61050..f16a36924 100644 --- a/cpp/src/mip/presolve/trivial_presolve.cuh +++ b/cpp/src/mip/presolve/trivial_presolve.cuh @@ -19,7 +19,6 @@ #include #include -#include #include #include #include @@ -346,8 +345,6 @@ void test_reverse_matches(const problem_t& pb) template void trivial_presolve(problem_t& problem) { - gf2_presolve(problem); - cuopt_expects(problem.preprocess_called, error_type_t::RuntimeError, "preprocess_problem should be called before running the solver"); From 182ac189a45118de6b606228ca28ce68c51642c1 Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Sat, 20 Sep 2025 11:12:13 +0000 Subject: [PATCH 4/9] cleanup --- .../linear_programming/cuopt/run_mip.cpp | 1 - .../mip/presolve/bounds_update_helpers.cuh | 4 +- cpp/src/mip/presolve/gf2_presolve.cpp | 118 +++++++++++++----- cpp/src/mip/presolve/gf2_presolve.hpp | 104 +-------------- cpp/src/mip/presolve/trivial_presolve.cuh | 3 - 5 files changed, 89 insertions(+), 141 deletions(-) diff --git a/benchmarks/linear_programming/cuopt/run_mip.cpp b/benchmarks/linear_programming/cuopt/run_mip.cpp index 8ce9fe5e7..64b1bb4bb 100644 --- a/benchmarks/linear_programming/cuopt/run_mip.cpp +++ b/benchmarks/linear_programming/cuopt/run_mip.cpp @@ -211,7 +211,6 @@ int run_single_file(std::string file_path, settings.tolerances.relative_tolerance = 1e-12; settings.tolerances.absolute_tolerance = 1e-6; settings.presolve = true; - settings.mip_scaling = false; cuopt::linear_programming::benchmark_info_t benchmark_info; settings.benchmark_info_ptr = &benchmark_info; auto start_run_solver = std::chrono::high_resolution_clock::now(); diff --git a/cpp/src/mip/presolve/bounds_update_helpers.cuh b/cpp/src/mip/presolve/bounds_update_helpers.cuh index 3de6b5537..77574d527 100644 --- a/cpp/src/mip/presolve/bounds_update_helpers.cuh +++ b/cpp/src/mip/presolve/bounds_update_helpers.cuh @@ -24,13 +24,13 @@ namespace cuopt::linear_programming::detail { // Activity calculation template -inline __host__ __device__ f_t min_act_of_var(f_t coeff, f_t var_lb, f_t var_ub) +inline __device__ f_t min_act_of_var(f_t coeff, f_t var_lb, f_t var_ub) { return (coeff < 0.) ? coeff * var_ub : coeff * var_lb; } template -inline __host__ __device__ f_t max_act_of_var(f_t coeff, f_t var_lb, f_t var_ub) +inline __device__ f_t max_act_of_var(f_t coeff, f_t var_lb, f_t var_ub) { return (coeff < 0.) ? coeff * var_lb : coeff * var_ub; } diff --git a/cpp/src/mip/presolve/gf2_presolve.cpp b/cpp/src/mip/presolve/gf2_presolve.cpp index 5fbc81fc1..244349dff 100644 --- a/cpp/src/mip/presolve/gf2_presolve.cpp +++ b/cpp/src/mip/presolve/gf2_presolve.cpp @@ -16,18 +16,78 @@ */ #include "gf2_presolve.hpp" + +#include + #include #include -#include +#if GF2_PRESOLVE_DEBUG #define NOT_GF2(reason, ...) \ - { \ + do { \ printf("NO : Cons %d is not gf2: " reason "\n", cstr_idx, ##__VA_ARGS__); \ goto not_valid; \ - } + } while (0) +#else +#define NOT_GF2(reason, ...) \ + do { \ + goto not_valid; \ + } while (0) +#endif namespace cuopt::linear_programming::detail { +// 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) +{ + int i, j, k; + const int N = A.size(); + for (i = 0; i < N; i++) { + // Find pivot + int pivot = -1; + for (j = i; j < N; j++) { + if (A[j][i]) { + pivot = j; + 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; + } + + // 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]; + } + } + } + + // 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 + } + return true; // Success +} + template papilo::PresolveStatus GF2Presolve::execute(const papilo::Problem& problem, const papilo::ProblemUpdate& problemUpdate, @@ -132,20 +192,14 @@ papilo::PresolveStatus GF2Presolve::execute(const papilo::Problem& pro } // If no GF2 constraints found, return unchanged - if (gf2_constraints.empty()) { - printf("No GF2 constraints found\n"); - exit(0); - return papilo::PresolveStatus::kUnchanged; - } + 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; } // Validate structure if (gf2_key_vars.size() != gf2_constraints.size() || gf2_bin_vars.size() != gf2_constraints.size()) { - printf("GF2 constraints: %d, %d, %d\n", - gf2_key_vars.size(), - gf2_bin_vars.size(), - gf2_constraints.size()); - exit(0); return papilo::PresolveStatus::kUnchanged; } @@ -156,7 +210,9 @@ papilo::PresolveStatus GF2Presolve::execute(const papilo::Problem& pro } // Build binary matrix - BinaryMatrix A(gf2_constraints.size(), gf2_constraints.size()); + // 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 (const auto& cons : gf2_constraints) { for (auto [bin_var, _] : cons.bin_vars) { @@ -165,32 +221,24 @@ papilo::PresolveStatus GF2Presolve::execute(const papilo::Problem& pro b[cons.cstr_idx] = cons.rhs; } - auto [inverse, rank] = A.inverse_and_rank(); - if (rank != (int)gf2_constraints.size()) { - printf("Non-invertible, rank: %d/%d\n", rank, gf2_constraints.size()); - exit(0); - return papilo::PresolveStatus::kUnchanged; // Non-invertible - } - - auto solution = inverse.dot(b); + std::vector solution(gf2_constraints.size()); + bool feasible = gf2_solve(A, b, solution); + if (!feasible) { 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]; } - // Compute fixings for key variables + // Compute fixings for key variables by solving for the constraint for (const auto& cons : gf2_constraints) { auto [key_var_idx, key_var_coeff] = cons.key_var; f_t constraint_rhs = lhs_values[cons.cstr_idx]; // equality constraint - f_t lhs = 0.0; + f_t lhs = -constraint_rhs; for (auto [bin_var, coeff] : cons.bin_vars) { lhs += fixings[bin_var] * coeff; } - lhs -= constraint_rhs; - fixings[key_var_idx] = std::round(-lhs / key_var_coeff); } @@ -204,13 +252,19 @@ papilo::PresolveStatus GF2Presolve::execute(const papilo::Problem& pro status = papilo::PresolveStatus::kReduced; } - printf("GF2 presolved!!\n"); - return status; } -// Explicit template instantiations -template class GF2Presolve; -template class GF2Presolve; +#define INSTANTIATE(F_TYPE) template class GF2Presolve; + +#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/presolve/gf2_presolve.hpp b/cpp/src/mip/presolve/gf2_presolve.hpp index 6d389b377..b7a6633aa 100644 --- a/cpp/src/mip/presolve/gf2_presolve.hpp +++ b/cpp/src/mip/presolve/gf2_presolve.hpp @@ -65,112 +65,10 @@ class GF2Presolve : public papilo::PresolveMethod { gf2_constraint_t& operator=(gf2_constraint_t&& other) noexcept = default; }; - class BinaryMatrix { - private: - std::vector> m_rows; - - public: - BinaryMatrix(void) { this->initialize(); } - BinaryMatrix(const int a_NUMBER_OF_ROWS, const int a_NUMBER_OF_COLUMNS) - { - this->setup(a_NUMBER_OF_ROWS, a_NUMBER_OF_COLUMNS); - } - - inline void initialize(void) { m_rows.clear(); } - inline void setup(const int a_NUMBER_OF_ROWS, const int a_NUMBER_OF_COLUMNS) - { - m_rows.resize(a_NUMBER_OF_ROWS, std::vector(a_NUMBER_OF_COLUMNS, 0)); - } - - inline int number_of_rows(void) const { return m_rows.size(); } - inline int number_of_columns(void) const { return m_rows.front().size(); } - inline std::vector& operator[](const int a_ROW) { return m_rows[a_ROW]; } - inline std::vector const& operator[](const int a_ROW) const { return m_rows[a_ROW]; } - - inline std::pair inverse_and_rank(void) const - { - const int SIZE = this->number_of_rows(); - BinaryMatrix A = *this; - BinaryMatrix B(SIZE, SIZE); - int rank = 0; - - for (auto i = 0; i < SIZE; i++) { - B[i][i] = 1; - } - - for (auto j = 0; j < SIZE; j++) { - int row = -1; - for (auto i = rank; i < SIZE; i++) { - if (A[i][j] == 1) { - row = i; - break; - } - } - if (row == -1) { continue; } - - if (row != rank) { - swap(A[row], A[rank]); - swap(B[row], B[rank]); - } - - for (auto i = rank + 1; i < SIZE; i++) { - if (A[i][j] == 1) { - for (auto k = 0; k < SIZE; k++) { - A[i][k] = (A[i][k] + A[rank][k]) & 1; - B[i][k] = (B[i][k] + B[rank][k]) & 1; - } - } - } - rank++; - } - - if (rank == SIZE) { - for (auto j = SIZE - 1; j > 0; j--) { - for (auto i = j - 1; i >= 0; i--) { - if (A[i][j] == 1) { - for (auto k = 0; k < SIZE; k++) { - A[i][k] = (A[j][k] + A[i][k]) & 1; - B[i][k] = (B[j][k] + B[i][k]) & 1; - } - } - } - } - } - - return {B, rank}; - } - - inline std::vector dot(const std::vector& a_VECTOR) const - { - std::vector result(m_rows.size(), 0); - const int NUMBER_OF_ROWS = this->number_of_rows(); - const int NUMBER_OF_COLUMNS = this->number_of_columns(); - for (auto i = 0; i < NUMBER_OF_ROWS; i++) { - for (auto j = 0; j < NUMBER_OF_COLUMNS; j++) { - result[i] += m_rows[i][j] * a_VECTOR[j]; - } - result[i] &= 1; - } - return result; - } - }; - - bool is_integer(f_t value, f_t tolerance) const + inline bool is_integer(f_t value, f_t tolerance) const { return std::abs(value - std::round(value)) <= tolerance; } - - f_t min_act_of_var(f_t coeff, f_t lb, f_t ub) const - { - if (coeff >= 0) return coeff * lb; - return coeff * ub; - } - - f_t max_act_of_var(f_t coeff, f_t lb, f_t ub) const - { - if (coeff >= 0) return coeff * ub; - return coeff * lb; - } }; } // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip/presolve/trivial_presolve.cuh b/cpp/src/mip/presolve/trivial_presolve.cuh index f16a36924..d7fe14233 100644 --- a/cpp/src/mip/presolve/trivial_presolve.cuh +++ b/cpp/src/mip/presolve/trivial_presolve.cuh @@ -356,9 +356,6 @@ void trivial_presolve(problem_t& problem) // The problem has been solved by presolve. Mark its empty status as valid if (problem.n_variables == 0) { problem.empty = true; } problem.check_problem_representation(true); - - problem.recompute_auxilliary_data( - false); // check problem representation later once cstr bounds are computed } } // namespace cuopt::linear_programming::detail From dc61271f0d373f836ae3f77688c155c7ed832c42 Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Mon, 22 Sep 2025 12:41:21 +0000 Subject: [PATCH 5/9] switch to an actually hard instance :) --- cpp/tests/linear_programming/c_api_tests/c_api_tests.cpp | 2 +- datasets/mip/download_miplib_test_dataset.sh | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/cpp/tests/linear_programming/c_api_tests/c_api_tests.cpp b/cpp/tests/linear_programming/c_api_tests/c_api_tests.cpp index d01f4319b..10d339748 100644 --- a/cpp/tests/linear_programming/c_api_tests/c_api_tests.cpp +++ b/cpp/tests/linear_programming/c_api_tests/c_api_tests.cpp @@ -71,7 +71,7 @@ INSTANTIATE_TEST_SUITE_P( 5, CUOPT_METHOD_DUAL_SIMPLEX), // LP, Dual Simplex std::make_tuple("/linear_programming/square41/square41.mps", 5, CUOPT_METHOD_PDLP), // LP, PDLP - std::make_tuple("/mip/enlight_hard.mps", 5, CUOPT_METHOD_DUAL_SIMPLEX) // MIP + std::make_tuple("/mip/supportcase22.mps", 5, CUOPT_METHOD_DUAL_SIMPLEX) // MIP )); TEST(c_api, iteration_limit) diff --git a/datasets/mip/download_miplib_test_dataset.sh b/datasets/mip/download_miplib_test_dataset.sh index d0eb0d2b9..1adc3a784 100755 --- a/datasets/mip/download_miplib_test_dataset.sh +++ b/datasets/mip/download_miplib_test_dataset.sh @@ -34,6 +34,7 @@ INSTANCES=( "neos5" "swath1" "enlight_hard" + "supportcase22" ) BASE_URL="https://miplib.zib.de/WebData/instances" From 3e1fe0c7c2991d70d3f0ca196d0815fac17f3a16 Mon Sep 17 00:00:00 2001 From: Alice Boucher <160623740+aliceb-nv@users.noreply.github.com> Date: Mon, 22 Sep 2025 17:52:50 +0200 Subject: [PATCH 6/9] Update cpp/src/mip/presolve/gf2_presolve.cpp Co-authored-by: Rajesh Gandham --- cpp/src/mip/presolve/gf2_presolve.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/cpp/src/mip/presolve/gf2_presolve.cpp b/cpp/src/mip/presolve/gf2_presolve.cpp index 244349dff..e612fc404 100644 --- a/cpp/src/mip/presolve/gf2_presolve.cpp +++ b/cpp/src/mip/presolve/gf2_presolve.cpp @@ -243,6 +243,7 @@ papilo::PresolveStatus GF2Presolve::execute(const papilo::Problem& pro } papilo::PresolveStatus status = papilo::PresolveStatus::kUnchanged; + TransactionGuard rg{reductions}; for (const auto& [var_idx, fixing] : fixings) { if (num.isZero(fixing)) { reductions.fixCol(var_idx, 0); From 6a955ea0bec5a2123a76cccf1a018f0f68d65243 Mon Sep 17 00:00:00 2001 From: Alice Boucher <160623740+aliceb-nv@users.noreply.github.com> Date: Mon, 22 Sep 2025 17:53:00 +0200 Subject: [PATCH 7/9] Update cpp/src/mip/presolve/third_party_presolve.cpp Co-authored-by: Rajesh Gandham --- cpp/src/mip/presolve/third_party_presolve.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/cpp/src/mip/presolve/third_party_presolve.cpp b/cpp/src/mip/presolve/third_party_presolve.cpp index 46c2ea4bc..a6c8f8659 100644 --- a/cpp/src/mip/presolve/third_party_presolve.cpp +++ b/cpp/src/mip/presolve/third_party_presolve.cpp @@ -310,6 +310,7 @@ void set_presolve_methods(papilo::Presolve& presolver, problem_category_t c using uptr = std::unique_ptr>; // cuopt custom presolvers +if (category == problem_category_t::MIP) presolver.addPresolveMethod(uptr(new cuopt::linear_programming::detail::GF2Presolve())); // fast presolvers From 7b28b242d3f1c969886082d604e7b6512be6a2c2 Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Mon, 22 Sep 2025 16:22:00 +0000 Subject: [PATCH 8/9] transaction guard, fix build --- cpp/src/mip/presolve/gf2_presolve.cpp | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/cpp/src/mip/presolve/gf2_presolve.cpp b/cpp/src/mip/presolve/gf2_presolve.cpp index e612fc404..c396a12c1 100644 --- a/cpp/src/mip/presolve/gf2_presolve.cpp +++ b/cpp/src/mip/presolve/gf2_presolve.cpp @@ -168,13 +168,6 @@ papilo::PresolveStatus GF2Presolve::execute(const papilo::Problem& pro gf2_key_vars.insert({var_idx, gf2_key_vars.size()}); } else { // Binary variable - f_t var_lb = col_flags[var_idx].test(papilo::ColFlag::kLbInf) - ? -std::numeric_limits::infinity() - : lower_bounds[var_idx]; - f_t var_ub = col_flags[var_idx].test(papilo::ColFlag::kUbInf) - ? std::numeric_limits::infinity() - : upper_bounds[var_idx]; - constraint_bin_vars.push_back({var_idx, coeff}); gf2_bin_vars.insert({var_idx, gf2_bin_vars.size()}); } @@ -243,7 +236,7 @@ papilo::PresolveStatus GF2Presolve::execute(const papilo::Problem& pro } papilo::PresolveStatus status = papilo::PresolveStatus::kUnchanged; - TransactionGuard rg{reductions}; + papilo::TransactionGuard rg{reductions}; for (const auto& [var_idx, fixing] : fixings) { if (num.isZero(fixing)) { reductions.fixCol(var_idx, 0); From 1588c156964b508b1a61bafb536e0a19a7130eff Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Mon, 22 Sep 2025 16:24:24 +0000 Subject: [PATCH 9/9] style fix --- cpp/src/mip/presolve/third_party_presolve.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cpp/src/mip/presolve/third_party_presolve.cpp b/cpp/src/mip/presolve/third_party_presolve.cpp index b067484de..dc2d4b00e 100644 --- a/cpp/src/mip/presolve/third_party_presolve.cpp +++ b/cpp/src/mip/presolve/third_party_presolve.cpp @@ -309,8 +309,8 @@ void set_presolve_methods(papilo::Presolve& presolver, problem_category_t c using uptr = std::unique_ptr>; // cuopt custom presolvers -if (category == problem_category_t::MIP) - presolver.addPresolveMethod(uptr(new cuopt::linear_programming::detail::GF2Presolve())); + if (category == problem_category_t::MIP) + presolver.addPresolveMethod(uptr(new cuopt::linear_programming::detail::GF2Presolve())); // fast presolvers presolver.addPresolveMethod(uptr(new papilo::SingletonCols()));