From 1e47561f9c03d673d00b12e9da91209e888de83f Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Tue, 7 Oct 2025 08:10:03 +0000 Subject: [PATCH 1/8] more robust feasibility checking --- .../feasibility_jump_impl_common.cuh | 18 +++++ .../feasibility_jump_kernels.cu | 23 +----- cpp/src/mip/feasibility_jump/fj_cpu.cu | 70 +++++++++++++------ 3 files changed, 71 insertions(+), 40 deletions(-) diff --git a/cpp/src/mip/feasibility_jump/feasibility_jump_impl_common.cuh b/cpp/src/mip/feasibility_jump/feasibility_jump_impl_common.cuh index 71513a2ec9..013aedbf79 100644 --- a/cpp/src/mip/feasibility_jump/feasibility_jump_impl_common.cuh +++ b/cpp/src/mip/feasibility_jump/feasibility_jump_impl_common.cuh @@ -19,6 +19,24 @@ namespace cuopt::linear_programming::detail { +template +HDI f_t fj_kahan_babushka_neumaier_sum(Iterator begin, Iterator end) +{ + f_t sum = 0; + f_t c = 0; + for (Iterator it = begin; it != end; ++it) { + f_t delta = *it; + f_t t = sum + delta; + if (fabs(sum) > fabs(delta)) { + c += (sum - t) + delta; + } else { + c += (delta - t) + sum; + } + sum = t; + } + return sum + c; +} + // Returns the current slack, and the variable delta that would nullify this slack ("tighten" it) template HDI thrust::tuple get_mtm_for_bound( diff --git a/cpp/src/mip/feasibility_jump/feasibility_jump_kernels.cu b/cpp/src/mip/feasibility_jump/feasibility_jump_kernels.cu index dc63d94cd9..4cb8bbd1b1 100644 --- a/cpp/src/mip/feasibility_jump/feasibility_jump_kernels.cu +++ b/cpp/src/mip/feasibility_jump/feasibility_jump_kernels.cu @@ -35,24 +35,6 @@ namespace cg = cooperative_groups; namespace cuopt::linear_programming::detail { -template -static DI f_t KahanBabushkaNeumaierSum(Iterator begin, Iterator end) -{ - f_t sum = 0; - f_t c = 0; - for (Iterator it = begin; it != end; ++it) { - f_t delta = *it; - f_t t = sum + delta; - if (fabs(sum) > fabs(delta)) { - c += (sum - t) + delta; - } else { - c += (delta - t) + sum; - } - sum = t; - } - return sum + c; -} - template DI thrust::pair move_objective_score( const typename fj_t::climber_data_t::view_t& fj, i_t var_idx, f_t delta) @@ -176,7 +158,7 @@ __global__ void init_lhs_and_violation(typename fj_t::climber_data_t:: return fj.pb.coefficients[j] * fj.incumbent_assignment[fj.pb.variables[j]]; }); fj.incumbent_lhs[cstr_idx] = - KahanBabushkaNeumaierSum(delta_it + offset_begin, delta_it + offset_end); + fj_kahan_babushka_neumaier_sum(delta_it + offset_begin, delta_it + offset_end); fj.incumbent_lhs_sumcomp[cstr_idx] = 0; f_t th_violation = fj.excess_score(cstr_idx, fj.incumbent_lhs[cstr_idx]); @@ -457,7 +439,8 @@ DI bool check_feasibility(const typename fj_t::climber_data_t::view_t& return fj.pb.coefficients[j] * fj.incumbent_assignment[fj.pb.variables[j]]; }); - f_t lhs = KahanBabushkaNeumaierSum(delta_it + offset_begin, delta_it + offset_end); + f_t lhs = + fj_kahan_babushka_neumaier_sum(delta_it + offset_begin, delta_it + offset_end); cuopt_assert(fj.cstr_satisfied(cIdx, lhs), "constraint violated"); } cuopt_func_call((check_variable_feasibility(fj, check_integer))); diff --git a/cpp/src/mip/feasibility_jump/fj_cpu.cu b/cpp/src/mip/feasibility_jump/fj_cpu.cu index ff4f1a1eb1..86033a8dfa 100644 --- a/cpp/src/mip/feasibility_jump/fj_cpu.cu +++ b/cpp/src/mip/feasibility_jump/fj_cpu.cu @@ -124,6 +124,20 @@ static inline bool tabu_check(fj_cpu_climber_t& fj_cpu, } } +template +static bool check_variable_feasibility(const typename fj_t::climber_data_t::view_t& fj, + bool check_integer = true) +{ + for (i_t var_idx = 0; var_idx < fj.pb.n_variables; var_idx += 1) { + auto val = fj.incumbent_assignment[var_idx]; + bool feasible = fj.pb.check_variable_within_bounds(var_idx, val); + + if (!feasible) return false; + if (check_integer && fj.pb.is_integer_var(var_idx)) return false; + } + return true; +} + template static inline fj_staged_score_t compute_score(fj_cpu_climber_t& fj_cpu, i_t var_idx, @@ -321,17 +335,22 @@ static void apply_move(fj_cpu_climber_t& fj_cpu, fj_cpu.h_incumbent_objective += fj_cpu.h_obj_coeffs[var_idx] * delta; if (fj_cpu.h_incumbent_objective < fj_cpu.h_best_objective && fj_cpu.violated_constraints.empty()) { - cuopt_assert(fj_cpu.satisfied_constraints.size() == fj_cpu.view.pb.n_constraints, ""); - fj_cpu.h_best_objective = - fj_cpu.h_incumbent_objective - fj_cpu.settings.parameters.breakthrough_move_epsilon; - fj_cpu.h_best_assignment = fj_cpu.h_assignment; - CUOPT_LOG_TRACE("%sCPUFJ: new best objective: %g\n", - fj_cpu.log_prefix.c_str(), - fj_cpu.pb_ptr->get_user_obj_from_solver_obj(fj_cpu.h_best_objective)); - if (fj_cpu.improvement_callback) { - fj_cpu.improvement_callback(fj_cpu.h_best_objective, fj_cpu.h_assignment); + // recompute the LHS values to cancel out accumulation errors, then check if feasibility remains + recompute_lhs(fj_cpu); + + if (!fj_cpu.violated_constraints.empty() && check_variable_feasibility(fj_cpu.view)) { + cuopt_assert(fj_cpu.satisfied_constraints.size() == fj_cpu.view.pb.n_constraints, ""); + fj_cpu.h_best_objective = + fj_cpu.h_incumbent_objective - fj_cpu.settings.parameters.breakthrough_move_epsilon; + fj_cpu.h_best_assignment = fj_cpu.h_assignment; + CUOPT_LOG_TRACE("%sCPUFJ: new best objective: %g\n", + fj_cpu.log_prefix.c_str(), + fj_cpu.pb_ptr->get_user_obj_from_solver_obj(fj_cpu.h_best_objective)); + if (fj_cpu.improvement_callback) { + fj_cpu.improvement_callback(fj_cpu.h_best_objective, fj_cpu.h_assignment); + } + fj_cpu.feasible_found = true; } - fj_cpu.feasible_found = true; } i_t tabu_tenure = fj_cpu.settings.parameters.tabu_tenure_min + @@ -537,7 +556,7 @@ static thrust::tuple find_mtm_move_sat( } template -static void init_lhs(fj_cpu_climber_t& fj_cpu) +static void recompute_lhs(fj_cpu_climber_t& fj_cpu) { cuopt_assert(fj_cpu.h_lhs.size() == fj_cpu.view.pb.n_constraints, "h_lhs size mismatch"); @@ -545,15 +564,16 @@ static void init_lhs(fj_cpu_climber_t& fj_cpu) fj_cpu.satisfied_constraints.clear(); for (i_t cstr_idx = 0; cstr_idx < fj_cpu.view.pb.n_constraints; ++cstr_idx) { auto [offset_begin, offset_end] = fj_cpu.view.pb.range_for_constraint(cstr_idx); - f_t lhs = 0; - for (i_t i = offset_begin; i < offset_end; ++i) { - lhs += fj_cpu.h_coefficients[i] * fj_cpu.h_assignment[fj_cpu.h_variables[i]]; - } - - fj_cpu.h_lhs[cstr_idx] = lhs; + auto delta_it = + thrust::make_transform_iterator(thrust::make_counting_iterator(0), [fj = fj_cpu.view](i_t j) { + return fj.pb.coefficients[j] * fj.incumbent_assignment[fj.pb.variables[j]]; + }); + fj_cpu.h_lhs[cstr_idx] = + fj_kahan_babushka_neumaier_sum(delta_it + offset_begin, delta_it + offset_end); + fj_cpu.h_lhs_sumcomp[cstr_idx] = 0; f_t cstr_tolerance = fj_cpu.view.get_corrected_tolerance(cstr_idx); - f_t new_cost = fj_cpu.view.excess_score(cstr_idx, lhs); + f_t new_cost = fj_cpu.view.excess_score(cstr_idx, fj_cpu.h_lhs[cstr_idx]); if (new_cost < -cstr_tolerance) { fj_cpu.violated_constraints.insert(cstr_idx); } else { @@ -695,7 +715,7 @@ static void perturb(fj_cpu_climber_t& fj_cpu) fj_cpu.h_assignment[var_idx] = val; } - init_lhs(fj_cpu); + recompute_lhs(fj_cpu); } template @@ -759,6 +779,8 @@ static void init_fj_cpu(fj_cpu_climber_t& fj_cpu, fj_cpu.view.incumbent_assignment = raft::device_span(fj_cpu.h_assignment.data(), fj_cpu.h_assignment.size()); fj_cpu.view.incumbent_lhs = raft::device_span(fj_cpu.h_lhs.data(), fj_cpu.h_lhs.size()); + fj_cpu.view.incumbent_lhs_sumcomp = + raft::device_span(fj_cpu.h_lhs_sumcomp.data(), fj_cpu.h_lhs_sumcomp.size()); fj_cpu.view.tabu_nodec_until = raft::device_span(fj_cpu.h_tabu_nodec_until.data(), fj_cpu.h_tabu_nodec_until.size()); fj_cpu.view.tabu_noinc_until = @@ -824,7 +846,7 @@ static void init_fj_cpu(fj_cpu_climber_t& fj_cpu, fj_cpu.var_bitmap.resize(fj_cpu.view.pb.n_variables, false); fj_cpu.iter_mtm_vars.reserve(fj_cpu.view.pb.n_variables); - init_lhs(fj_cpu); + recompute_lhs(fj_cpu); } template @@ -914,6 +936,14 @@ bool fj_t::cpu_solve(fj_cpu_climber_t& fj_cpu, f_t in_time_l break; } + // periodically recompute the LHS and violation scores + // to correct any accumulated numerical errors + cuopt_assert(fj_cpu.settings.parameters.lhs_refresh_period > 0, + "lhs_refresh_period should be positive"); + if (fj_cpu.iterations % fj_cpu.settings.parameters.lhs_refresh_period == 0) { + recompute_lhs(fj_cpu); + } + fj_move_t move = fj_move_t{-1, 0}; fj_staged_score_t score = fj_staged_score_t::invalid(); // Perform lift moves From 2a186658791ff1d1097cca7332240a9aa1dd4a09 Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Tue, 7 Oct 2025 08:12:37 +0000 Subject: [PATCH 2/8] bump 1 --- benchmarks/linear_programming/cuopt/run_mip.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/linear_programming/cuopt/run_mip.cpp b/benchmarks/linear_programming/cuopt/run_mip.cpp index a4f52cb4ed..52168277c6 100644 --- a/benchmarks/linear_programming/cuopt/run_mip.cpp +++ b/benchmarks/linear_programming/cuopt/run_mip.cpp @@ -299,7 +299,7 @@ int main(int argc, char* argv[]) { argparse::ArgumentParser program("solve_MIP"); - // Define all arguments with appropriate defaults and help messages + // Define all arguments with appropriate defaults and help messages 1 program.add_argument("--path").help("input path").required(); program.add_argument("--run-dir") From 9924b6d830e46fc9558a59021728d8246abd90a9 Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Tue, 7 Oct 2025 08:14:05 +0000 Subject: [PATCH 3/8] bump 2 --- benchmarks/linear_programming/cuopt/run_mip.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/linear_programming/cuopt/run_mip.cpp b/benchmarks/linear_programming/cuopt/run_mip.cpp index 52168277c6..a4f52cb4ed 100644 --- a/benchmarks/linear_programming/cuopt/run_mip.cpp +++ b/benchmarks/linear_programming/cuopt/run_mip.cpp @@ -299,7 +299,7 @@ int main(int argc, char* argv[]) { argparse::ArgumentParser program("solve_MIP"); - // Define all arguments with appropriate defaults and help messages 1 + // Define all arguments with appropriate defaults and help messages program.add_argument("--path").help("input path").required(); program.add_argument("--run-dir") From 3e83a345f5426cbba3cda64a8b8d103907b2064d Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Tue, 7 Oct 2025 08:20:17 +0000 Subject: [PATCH 4/8] fix incorrect variable feasibility check --- cpp/src/mip/feasibility_jump/fj_cpu.cu | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/cpp/src/mip/feasibility_jump/fj_cpu.cu b/cpp/src/mip/feasibility_jump/fj_cpu.cu index 86033a8dfa..783e020a53 100644 --- a/cpp/src/mip/feasibility_jump/fj_cpu.cu +++ b/cpp/src/mip/feasibility_jump/fj_cpu.cu @@ -133,7 +133,9 @@ static bool check_variable_feasibility(const typename fj_t::climber_da bool feasible = fj.pb.check_variable_within_bounds(var_idx, val); if (!feasible) return false; - if (check_integer && fj.pb.is_integer_var(var_idx)) return false; + if (check_integer && fj.pb.is_integer_var(var_idx) && + !fj.pb.is_integer(fj.incumbent_assignment[var_idx])) + return false; // 1 } return true; } From 8dfbbef522f942450820c66a46a1ff374085ba4e Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Tue, 7 Oct 2025 08:20:33 +0000 Subject: [PATCH 5/8] bump --- cpp/src/mip/feasibility_jump/fj_cpu.cu | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/src/mip/feasibility_jump/fj_cpu.cu b/cpp/src/mip/feasibility_jump/fj_cpu.cu index 783e020a53..ad88c8d059 100644 --- a/cpp/src/mip/feasibility_jump/fj_cpu.cu +++ b/cpp/src/mip/feasibility_jump/fj_cpu.cu @@ -135,7 +135,7 @@ static bool check_variable_feasibility(const typename fj_t::climber_da if (!feasible) return false; if (check_integer && fj.pb.is_integer_var(var_idx) && !fj.pb.is_integer(fj.incumbent_assignment[var_idx])) - return false; // 1 + return false; } return true; } From 931441fdf3a803b5765aab248b96370d853a62ef Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Tue, 7 Oct 2025 08:26:07 +0000 Subject: [PATCH 6/8] fix typo --- cpp/src/mip/feasibility_jump/fj_cpu.cu | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/cpp/src/mip/feasibility_jump/fj_cpu.cu b/cpp/src/mip/feasibility_jump/fj_cpu.cu index ad88c8d059..c77f65287e 100644 --- a/cpp/src/mip/feasibility_jump/fj_cpu.cu +++ b/cpp/src/mip/feasibility_jump/fj_cpu.cu @@ -338,9 +338,10 @@ static void apply_move(fj_cpu_climber_t& fj_cpu, if (fj_cpu.h_incumbent_objective < fj_cpu.h_best_objective && fj_cpu.violated_constraints.empty()) { // recompute the LHS values to cancel out accumulation errors, then check if feasibility remains + // 1 recompute_lhs(fj_cpu); - if (!fj_cpu.violated_constraints.empty() && check_variable_feasibility(fj_cpu.view)) { + if (fj_cpu.violated_constraints.empty() && check_variable_feasibility(fj_cpu.view)) { cuopt_assert(fj_cpu.satisfied_constraints.size() == fj_cpu.view.pb.n_constraints, ""); fj_cpu.h_best_objective = fj_cpu.h_incumbent_objective - fj_cpu.settings.parameters.breakthrough_move_epsilon; From ae1ac963332d0fb146be9411cb53a2f37b1f6b5c Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Tue, 7 Oct 2025 08:26:23 +0000 Subject: [PATCH 7/8] bump --- cpp/src/mip/feasibility_jump/fj_cpu.cu | 1 - 1 file changed, 1 deletion(-) diff --git a/cpp/src/mip/feasibility_jump/fj_cpu.cu b/cpp/src/mip/feasibility_jump/fj_cpu.cu index c77f65287e..df839ab667 100644 --- a/cpp/src/mip/feasibility_jump/fj_cpu.cu +++ b/cpp/src/mip/feasibility_jump/fj_cpu.cu @@ -338,7 +338,6 @@ static void apply_move(fj_cpu_climber_t& fj_cpu, if (fj_cpu.h_incumbent_objective < fj_cpu.h_best_objective && fj_cpu.violated_constraints.empty()) { // recompute the LHS values to cancel out accumulation errors, then check if feasibility remains - // 1 recompute_lhs(fj_cpu); if (fj_cpu.violated_constraints.empty() && check_variable_feasibility(fj_cpu.view)) { From 022da68b84f8da92fe967f0317d0b74e99e8bad5 Mon Sep 17 00:00:00 2001 From: Alice Boucher Date: Tue, 7 Oct 2025 14:33:28 +0000 Subject: [PATCH 8/8] fix fickle CI test --- cpp/tests/linear_programming/c_api_tests/c_api_tests.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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 6a10193079..d72f53fb26 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 @@ -60,7 +60,7 @@ TEST_P(TimeLimitTestFixture, time_limit) // Dual simplex is spending some time for factorizing the basis, and this computation does not // check for time limit - double excess_allowed_time = 2.0; + double excess_allowed_time = 3.0; EXPECT_NEAR(solve_time, target_solve_time, excess_allowed_time); } INSTANTIATE_TEST_SUITE_P( @@ -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/supportcase22.mps", 5, CUOPT_METHOD_DUAL_SIMPLEX) // MIP + std::make_tuple("/mip/supportcase22.mps", 15, CUOPT_METHOD_DUAL_SIMPLEX) // MIP )); TEST(c_api, iteration_limit)