From 305153158268aed04bafdb98b94423005d07a9f4 Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Tue, 17 Mar 2026 16:33:24 -0700 Subject: [PATCH 01/14] Complement variables in order to find knapsack cuts --- cpp/src/branch_and_bound/branch_and_bound.cpp | 2 +- cpp/src/cuts/cuts.cpp | 60 ++++++++++++++----- cpp/src/cuts/cuts.hpp | 3 +- 3 files changed, 49 insertions(+), 16 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index ba2b63983..9113f5c4b 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -2220,7 +2220,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut i_t num_cuts = cut_pool.get_best_cuts(cuts_to_add, cut_rhs, cut_types); if (num_cuts == 0) { break; } cut_info.record_cut_types(cut_types); -#ifdef PRINT_CUT_POOL_TYPES +#if 1 cut_pool.print_cutpool_types(); print_cut_types("In LP ", cut_types, settings_); printf("Cut pool size: %d\n", cut_pool.pool_size()); diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index c74a12b36..7a9b44a60 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -697,6 +697,7 @@ knapsack_generation_t::knapsack_generation_t( if (row_len < 3) { continue; } bool is_knapsack = true; f_t sum_pos = 0.0; + f_t sum_neg = 0.0; for (i_t p = row_start; p < row_end; p++) { const i_t j = Arow.j[p]; if (is_slack_[j]) { continue; } @@ -710,22 +711,24 @@ knapsack_generation_t::knapsack_generation_t( break; } if (aj < 0.0) { - is_knapsack = false; - break; + sum_pos += -aj; + sum_neg += -aj; + } else { + sum_pos += aj; } - sum_pos += aj; } if (is_knapsack) { - const f_t beta = lp.rhs[i]; + const f_t beta = lp.rhs[i] + sum_neg; if (std::abs(beta - std::round(beta)) <= settings.integer_tol) { if (beta > 0.0 && beta <= sum_pos && std::abs(sum_pos / (row_len - 1) - beta) > 1e-3) { - if (verbose) { + if (1) { settings.log.printf( - "Knapsack constraint %d row len %d beta %e sum_pos %e sum_pos / (row_len - 1) %e\n", + "Knapsack constraint %d row len %d beta %e sum_neg %e sum_pos %e sum_pos / (row_len - 1) %e\n", i, row_len, beta, + sum_neg, sum_pos, sum_pos / (row_len - 1)); } @@ -735,14 +738,14 @@ knapsack_generation_t::knapsack_generation_t( } } -#ifdef PRINT_KNAPSACK_INFO +#if 1 i_t num_knapsack_constraints = knapsack_constraints_.size(); settings.log.printf("Number of knapsack constraints %d\n", num_knapsack_constraints); #endif } template -i_t knapsack_generation_t::generate_knapsack_cuts( +i_t knapsack_generation_t::generate_knapsack_cut( const lp_problem_t& lp, const simplex_solver_settings_t& settings, csr_matrix_t& Arow, @@ -750,6 +753,7 @@ i_t knapsack_generation_t::generate_knapsack_cuts( const std::vector& var_types, const std::vector& xstar, i_t knapsack_row, + std::vector& is_complemented, inequality_t& cut) { const bool verbose = false; @@ -760,11 +764,20 @@ i_t knapsack_generation_t::generate_knapsack_cuts( // Remove the slacks from the inequality f_t seperation_rhs = 0.0; if (verbose) { settings.log.printf(" Knapsack : "); } + std::vector complemented_variables; + complemented_variables.reserve(knapsack_inequality.i.size()); for (i_t k = 0; k < knapsack_inequality.i.size(); k++) { const i_t j = knapsack_inequality.i[k]; if (is_slack_[j]) { knapsack_inequality.x[k] = 0.0; } else { + const f_t aj = knapsack_inequality.x[k]; + if (aj < 0.0) { + knapsack_rhs -= aj; + knapsack_inequality.x[k] *= -1.0; + complemented_variables.push_back(j); + is_complemented[j] = 1; + } if (verbose) { settings.log.printf(" %g x%d +", knapsack_inequality.x[k], j); } seperation_rhs += knapsack_inequality.x[k]; } @@ -784,7 +797,15 @@ i_t knapsack_generation_t::generate_knapsack_cuts( settings.log.printf("seperation_rhs %g\n", seperation_rhs); } - if (seperation_rhs <= 0.0) { return -1; } + auto restore_complemented = [&complemented_variables, &is_complemented]() { + for (i_t j : complemented_variables) { + is_complemented[j] = 0; + } + }; + if (seperation_rhs <= 0.0) { + restore_complemented(); + return -1; + } std::vector values; values.resize(knapsack_inequality.i.size() - 1); @@ -795,7 +816,8 @@ i_t knapsack_generation_t::generate_knapsack_cuts( for (i_t k = 0; k < knapsack_inequality.i.size(); k++) { const i_t j = knapsack_inequality.i[k]; if (!is_slack_[j]) { - const f_t vj = std::min(1.0, std::max(0.0, 1.0 - xstar[j])); + const f_t xstar_j = is_complemented[j] ? 1.0 - xstar[j] : xstar[j]; + const f_t vj = std::min(1.0, std::max(0.0, 1.0 - xstar_j)); objective_constant += vj; values[h] = vj; weights[h] = knapsack_inequality.x[k]; @@ -807,14 +829,14 @@ i_t knapsack_generation_t::generate_knapsack_cuts( if (verbose) { settings.log.printf("Calling solve_knapsack_problem\n"); } f_t objective = solve_knapsack_problem(values, weights, seperation_rhs, solution); - if (std::isnan(objective)) { return -1; } + if (std::isnan(objective)) { restore_complemented(); return -1; } if (verbose) { settings.log.printf("objective %e objective_constant %e\n", objective, objective_constant); } f_t seperation_value = -objective + objective_constant; if (verbose) { settings.log.printf("seperation_value %e\n", seperation_value); } const f_t tol = 1e-6; - if (seperation_value >= 1.0 - tol) { return -1; } + if (seperation_value >= 1.0 - tol) { restore_complemented(); return -1; } i_t cover_size = 0; for (i_t k = 0; k < solution.size(); k++) { @@ -833,6 +855,14 @@ i_t knapsack_generation_t::generate_knapsack_cuts( } } cut.rhs = -cover_size + 1; + + for (i_t k = 0; k < cut.size(); k++) { + const i_t j = cut.index(k); + if (is_complemented[j]) { + cut.vector.x[k] *= -1.0; + cut.rhs += 1.0; + } + } cut.sort(); // The cut is in the form: - sum_{j in cover} x_j >= -cover_size + 1 @@ -845,6 +875,7 @@ i_t knapsack_generation_t::generate_knapsack_cuts( settings.log.printf("Knapsack cut %d violation %e < 0\n", knapsack_row, violation); } + restore_complemented(); if (violation >= -tol) { return -1; } return 0; } @@ -1086,10 +1117,11 @@ void cut_generation_t::generate_knapsack_cuts( const std::vector& xstar) { if (knapsack_generation_.num_knapsack_constraints() > 0) { + std::vector is_complemented(lp.num_cols, 0); for (i_t knapsack_row : knapsack_generation_.get_knapsack_constraints()) { inequality_t cut(lp.num_cols); - i_t knapsack_status = knapsack_generation_.generate_knapsack_cuts( - lp, settings, Arow, new_slacks, var_types, xstar, knapsack_row, cut); + i_t knapsack_status = knapsack_generation_.generate_knapsack_cut( + lp, settings, Arow, new_slacks, var_types, xstar, knapsack_row, is_complemented, cut); if (knapsack_status == 0) { cut_pool_.add_cut(cut_type_t::KNAPSACK, cut); } } } diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index 91806d81a..e4f9c1932 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -288,13 +288,14 @@ class knapsack_generation_t { const std::vector& new_slacks, const std::vector& var_types); - i_t generate_knapsack_cuts(const lp_problem_t& lp, + i_t generate_knapsack_cut(const lp_problem_t& lp, const simplex_solver_settings_t& settings, csr_matrix_t& Arow, const std::vector& new_slacks, const std::vector& var_types, const std::vector& xstar, i_t knapsack_row, + std::vector& is_complemented, inequality_t& cut); i_t num_knapsack_constraints() const { return knapsack_constraints_.size(); } From c2a7331030a49fa74b07b37d0d0f0d54b153c3ec Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Tue, 17 Mar 2026 20:20:07 -0700 Subject: [PATCH 02/14] Convert coefficients to integers in rows when looking for knapsack constraints. Fix variables to reduce the size of seperation problem --- cpp/src/cuts/cuts.cpp | 159 ++++++++++++++++++++++++++++-------------- cpp/src/cuts/cuts.hpp | 28 ++++---- 2 files changed, 120 insertions(+), 67 deletions(-) diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index 7a9b44a60..f978e8a4a 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -691,21 +691,22 @@ knapsack_generation_t::knapsack_generation_t( } for (i_t i = 0; i < lp.num_rows; i++) { - const i_t row_start = Arow.row_start[i]; - const i_t row_end = Arow.row_start[i + 1]; - const i_t row_len = row_end - row_start; + inequality_t inequality(Arow, i, lp.rhs[i]); + inequality_t rational_inequality = inequality; + if (!rational_coefficients(var_types, inequality, rational_inequality)) { + continue; + } + inequality = rational_inequality; + + const i_t row_len = rational_inequality.size(); if (row_len < 3) { continue; } bool is_knapsack = true; f_t sum_pos = 0.0; f_t sum_neg = 0.0; - for (i_t p = row_start; p < row_end; p++) { - const i_t j = Arow.j[p]; + for (i_t p = 0; p < row_len; p++) { + const i_t j = inequality.index(p); if (is_slack_[j]) { continue; } - const f_t aj = Arow.x[p]; - if (std::abs(aj - std::round(aj)) > settings.integer_tol) { - is_knapsack = false; - break; - } + const f_t aj = inequality.coeff(p); if (var_types[j] != variable_type_t::INTEGER || lp.lower[j] != 0.0 || lp.upper[j] != 1.0) { is_knapsack = false; break; @@ -719,7 +720,7 @@ knapsack_generation_t::knapsack_generation_t( } if (is_knapsack) { - const f_t beta = lp.rhs[i] + sum_neg; + const f_t beta = inequality.rhs + sum_neg; if (std::abs(beta - std::round(beta)) <= settings.integer_tol) { if (beta > 0.0 && beta <= sum_pos && std::abs(sum_pos / (row_len - 1) - beta) > 1e-3) { if (1) { @@ -758,37 +759,61 @@ i_t knapsack_generation_t::generate_knapsack_cut( { const bool verbose = false; // Get the row associated with the knapsack constraint - sparse_vector_t knapsack_inequality(Arow, knapsack_row); - f_t knapsack_rhs = lp.rhs[knapsack_row]; + inequality_t knapsack_inequality(Arow, knapsack_row, lp.rhs[knapsack_row]); + inequality_t rational_knapsack_inequality = knapsack_inequality; + if (!rational_coefficients(var_types, knapsack_inequality, rational_knapsack_inequality)) { + return -1; + } + knapsack_inequality = rational_knapsack_inequality; + + // Given the following knapsack constraint: + // sum_j a_j x_j <= beta + // + // We solve the following separation problem: + // minimize sum_j (1 - xstar_j) z_j + // subject to sum_j a_j z_j > beta + // z_j in {0, 1} + // When z_j = 1, then j is in the cover. + // Let phi_star be the optimal objective of this problem. + // We have a violated cover when phi_star < 1.0 + // + // We convert this problem into a 0-1 knapsack problem + // maximize sum_j (1 - xstar_j) zbar_j + // subject to sum_j a_j zbar_j <= sum_j a_j - (beta + 1) + // zbar_j in {0, 1} + // where zbar_j = 1 - z_j + // This problem is in the form of a 0-1 knapsack problem + // which we can solve with dynamic programming or generate + // a heuristic solution with a greedy algorithm. // Remove the slacks from the inequality f_t seperation_rhs = 0.0; if (verbose) { settings.log.printf(" Knapsack : "); } std::vector complemented_variables; - complemented_variables.reserve(knapsack_inequality.i.size()); - for (i_t k = 0; k < knapsack_inequality.i.size(); k++) { - const i_t j = knapsack_inequality.i[k]; + complemented_variables.reserve(knapsack_inequality.size()); + for (i_t k = 0; k < knapsack_inequality.size(); k++) { + const i_t j = knapsack_inequality.index(k); if (is_slack_[j]) { - knapsack_inequality.x[k] = 0.0; + knapsack_inequality.vector.x[k] = 0.0; } else { - const f_t aj = knapsack_inequality.x[k]; + const f_t aj = knapsack_inequality.vector.x[k]; if (aj < 0.0) { - knapsack_rhs -= aj; - knapsack_inequality.x[k] *= -1.0; + knapsack_inequality.rhs -= aj; + knapsack_inequality.vector.x[k] *= -1.0; complemented_variables.push_back(j); is_complemented[j] = 1; } - if (verbose) { settings.log.printf(" %g x%d +", knapsack_inequality.x[k], j); } - seperation_rhs += knapsack_inequality.x[k]; + if (verbose) { settings.log.printf(" %g x%d +", knapsack_inequality.vector.x[k], j); } + seperation_rhs += knapsack_inequality.vector.x[k]; } } - if (verbose) { settings.log.printf(" <= %g\n", knapsack_rhs); } - seperation_rhs -= (knapsack_rhs + 1); + if (verbose) { settings.log.printf(" <= %g\n", knapsack_inequality.rhs); } + seperation_rhs -= (knapsack_inequality.rhs + 1); if (verbose) { settings.log.printf("\t"); - for (i_t k = 0; k < knapsack_inequality.i.size(); k++) { - const i_t j = knapsack_inequality.i[k]; + for (i_t k = 0; k < knapsack_inequality.size(); k++) { + const i_t j = knapsack_inequality.index(k); if (!is_slack_[j]) { if (std::abs(xstar[j]) > 1e-3) { settings.log.printf("x_relax[%d]= %g ", j, xstar[j]); } } @@ -808,24 +833,49 @@ i_t knapsack_generation_t::generate_knapsack_cut( } std::vector values; - values.resize(knapsack_inequality.i.size() - 1); + values.reserve(knapsack_inequality.size() - 1); std::vector weights; - weights.resize(knapsack_inequality.i.size() - 1); - i_t h = 0; + weights.reserve(knapsack_inequality.size() - 1); + std::vector indices; + indices.reserve(knapsack_inequality.size() - 1); f_t objective_constant = 0.0; - for (i_t k = 0; k < knapsack_inequality.i.size(); k++) { - const i_t j = knapsack_inequality.i[k]; + std::vector fixed_variables; + std::vector fixed_values; + const f_t x_tol = 1e-5; + for (i_t k = 0; k < knapsack_inequality.size(); k++) { + const i_t j = knapsack_inequality.index(k); if (!is_slack_[j]) { const f_t xstar_j = is_complemented[j] ? 1.0 - xstar[j] : xstar[j]; const f_t vj = std::min(1.0, std::max(0.0, 1.0 - xstar_j)); + if (xstar_j < x_tol) { + // if xstar_j is close to 0, then we can fix z to zero + fixed_variables.push_back(j); + fixed_values.push_back(0.0); + seperation_rhs -= knapsack_inequality.vector.x[k]; + // No need to adjust the objective constant + continue; + } + if (xstar_j > 1.0 - x_tol) { + // if xstar_j is close to 1, then we can fix z to 1 + fixed_variables.push_back(j); + fixed_values.push_back(1.0); + // Note seperation rhs is unchanged + objective_constant += vj; + continue; + } objective_constant += vj; - values[h] = vj; - weights[h] = knapsack_inequality.x[k]; - h++; + values.push_back(vj); + weights.push_back(knapsack_inequality.vector.x[k]); + indices.push_back(j); } } std::vector solution; - solution.resize(knapsack_inequality.i.size() - 1); + solution.resize(values.size()); + + if (seperation_rhs <= 0.0) { + restore_complemented(); + return -1; + } if (verbose) { settings.log.printf("Calling solve_knapsack_problem\n"); } f_t objective = solve_knapsack_problem(values, weights, seperation_rhs, solution); @@ -842,17 +892,20 @@ i_t knapsack_generation_t::generate_knapsack_cut( for (i_t k = 0; k < solution.size(); k++) { if (solution[k] == 0.0) { cover_size++; } } + for (i_t k = 0; k < fixed_values.size(); k++) { + if (fixed_values[k] == 1.0) { cover_size++; } + } cut.reserve(cover_size); cut.clear(); - h = 0; - for (i_t k = 0; k < knapsack_inequality.i.size(); k++) { - const i_t j = knapsack_inequality.i[k]; - if (!is_slack_[j]) { - if (solution[h] == 0.0) { cut.push_back(j, -1.0); } - h++; - } + for (i_t k = 0; k < solution.size(); k++) { + const i_t j = indices[k]; + if (solution[k] == 0.0) { cut.push_back(j, -1.0); } + } + for (i_t k = 0; k < fixed_variables.size(); k++) { + const i_t j = fixed_variables[k]; + if (fixed_values[k] == 1.0) { cut.push_back(j, -1.0); } } cut.rhs = -cover_size + 1; @@ -1740,7 +1793,7 @@ void cut_generation_t::generate_gomory_cuts( complemented_mir.remove_small_coefficients(lp.lower, lp.upper, cut_A_float); inequality_t cut_A(lp.num_cols); - if (cut_ok) { cut_ok = gomory_cut.rational_coefficients(var_types, cut_A_float, cut_A); } + if (cut_ok) { cut_ok = rational_coefficients(var_types, cut_A_float, cut_A); } // See if the inequality is violated by the original relaxation solution f_t cut_A_violation = complemented_mir.compute_violation(cut_A, xstar); @@ -1777,7 +1830,7 @@ void cut_generation_t::generate_gomory_cuts( complemented_mir.remove_small_coefficients(lp.lower, lp.upper, cut_B_float); inequality_t cut_B(lp.num_cols); - if (cut_ok) { cut_ok = gomory_cut.rational_coefficients(var_types, cut_B_float, cut_B); } + if (cut_ok) { cut_ok = rational_coefficients(var_types, cut_B_float, cut_B); } bool B_valid = false; f_t cut_B_distance = 0.0; @@ -1962,11 +2015,11 @@ i_t tableau_equality_t::generate_base_equality( return 0; } -template -bool mixed_integer_gomory_cut_t::rational_approximation(f_t x, - int64_t max_denominator, - int64_t& numerator, - int64_t& denominator) +template +bool rational_approximation(f_t x, + int64_t max_denominator, + int64_t& numerator, + int64_t& denominator) { int64_t a, p0 = 0, q0 = 1, p1 = 1, q1 = 0; f_t val = x; @@ -2002,7 +2055,7 @@ bool mixed_integer_gomory_cut_t::rational_approximation(f_t x, } template -bool mixed_integer_gomory_cut_t::rational_coefficients( +bool rational_coefficients( const std::vector& var_types, const inequality_t& input_inequality, inequality_t& rational_inequality) @@ -2039,8 +2092,7 @@ bool mixed_integer_gomory_cut_t::rational_coefficients( return true; } -template -int64_t mixed_integer_gomory_cut_t::gcd(const std::vector& integers) +int64_t gcd(const std::vector& integers) { if (integers.empty()) { return 0; } @@ -2051,8 +2103,7 @@ int64_t mixed_integer_gomory_cut_t::gcd(const std::vector& in return result; } -template -int64_t mixed_integer_gomory_cut_t::lcm(const std::vector& integers) +int64_t lcm(const std::vector& integers) { if (integers.empty()) { return 0; } int64_t result = diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index e4f9c1932..b79ff1e8b 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -166,6 +166,21 @@ f_t fractional_part(f_t a) return a - std::floor(a); } +template +bool rational_coefficients(const std::vector& var_types, + const inequality_t& inequality, + inequality_t& rational_inequality); + +template +bool rational_approximation(f_t x, + int64_t max_denominator, + int64_t& numerator, + int64_t& denominator); + +int64_t gcd(const std::vector& integers); + +int64_t lcm(const std::vector& integers); + template bool add_work_estimate(f_t accesses, f_t* work_estimate, @@ -465,19 +480,6 @@ template class mixed_integer_gomory_cut_t { public: mixed_integer_gomory_cut_t() {} - - bool rational_coefficients(const std::vector& var_types, - const inequality_t& input_inequality, - inequality_t& rational_inequality); - - private: - bool rational_approximation(f_t x, - int64_t max_denominator, - int64_t& numerator, - int64_t& denominator); - - int64_t gcd(const std::vector& integers); - int64_t lcm(const std::vector& integers); }; template From b07feed81ab3099c369e134293aa66bb62b0b401 Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Wed, 18 Mar 2026 15:43:50 -0700 Subject: [PATCH 03/14] Lift knapsack cuts --- cpp/src/branch_and_bound/branch_and_bound.cpp | 6 - cpp/src/cuts/cuts.cpp | 402 ++++++++++++++++-- cpp/src/cuts/cuts.hpp | 27 +- 3 files changed, 399 insertions(+), 36 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 9113f5c4b..2133f156f 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -724,12 +724,6 @@ void branch_and_bound_t::set_final_solution(mip_solution_t& obj, is_maximization ? "Upper" : "Lower", user_bound); - { - const f_t root_lp_obj = root_lp_current_lower_bound_.load(); - if (std::isfinite(root_lp_obj)) { - settings_.log.printf("Root LP dual objective (last): %.16e\n", root_lp_obj); - } - } if (gap <= settings_.absolute_mip_gap_tol || gap_rel <= settings_.relative_mip_gap_tol) { solver_status_ = mip_status_t::OPTIMAL; diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index f978e8a4a..115fa82e6 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -680,12 +680,16 @@ knapsack_generation_t::knapsack_generation_t( csr_matrix_t& Arow, const std::vector& new_slacks, const std::vector& var_types) - : settings_(settings) + : is_slack_(lp.num_cols, 0), + is_complemented_(lp.num_cols, 0), + is_marked_(lp.num_cols, 0), + workspace_(lp.num_cols, 0.0), + complemented_xstar_(lp.num_cols, 0.0), + settings_(settings) { const bool verbose = false; knapsack_constraints_.reserve(lp.num_rows); - is_slack_.resize(lp.num_cols, 0); for (i_t j : new_slacks) { is_slack_[j] = 1; } @@ -754,7 +758,6 @@ i_t knapsack_generation_t::generate_knapsack_cut( const std::vector& var_types, const std::vector& xstar, i_t knapsack_row, - std::vector& is_complemented, inequality_t& cut) { const bool verbose = false; @@ -801,7 +804,7 @@ i_t knapsack_generation_t::generate_knapsack_cut( knapsack_inequality.rhs -= aj; knapsack_inequality.vector.x[k] *= -1.0; complemented_variables.push_back(j); - is_complemented[j] = 1; + is_complemented_[j] = 1; } if (verbose) { settings.log.printf(" %g x%d +", knapsack_inequality.vector.x[k], j); } seperation_rhs += knapsack_inequality.vector.x[k]; @@ -822,13 +825,9 @@ i_t knapsack_generation_t::generate_knapsack_cut( settings.log.printf("seperation_rhs %g\n", seperation_rhs); } - auto restore_complemented = [&complemented_variables, &is_complemented]() { - for (i_t j : complemented_variables) { - is_complemented[j] = 0; - } - }; + if (seperation_rhs <= 0.0) { - restore_complemented(); + restore_complemented(complemented_variables); return -1; } @@ -845,7 +844,8 @@ i_t knapsack_generation_t::generate_knapsack_cut( for (i_t k = 0; k < knapsack_inequality.size(); k++) { const i_t j = knapsack_inequality.index(k); if (!is_slack_[j]) { - const f_t xstar_j = is_complemented[j] ? 1.0 - xstar[j] : xstar[j]; + const f_t xstar_j = is_complemented_[j] ? 1.0 - xstar[j] : xstar[j]; + complemented_xstar_[j] = xstar_j; const f_t vj = std::min(1.0, std::max(0.0, 1.0 - xstar_j)); if (xstar_j < x_tol) { // if xstar_j is close to 0, then we can fix z to zero @@ -873,20 +873,20 @@ i_t knapsack_generation_t::generate_knapsack_cut( solution.resize(values.size()); if (seperation_rhs <= 0.0) { - restore_complemented(); + restore_complemented(complemented_variables); return -1; } if (verbose) { settings.log.printf("Calling solve_knapsack_problem\n"); } f_t objective = solve_knapsack_problem(values, weights, seperation_rhs, solution); - if (std::isnan(objective)) { restore_complemented(); return -1; } + if (std::isnan(objective)) { restore_complemented(complemented_variables); return -1; } if (verbose) { settings.log.printf("objective %e objective_constant %e\n", objective, objective_constant); } f_t seperation_value = -objective + objective_constant; if (verbose) { settings.log.printf("seperation_value %e\n", seperation_value); } const f_t tol = 1e-6; - if (seperation_value >= 1.0 - tol) { restore_complemented(); return -1; } + if (seperation_value >= 1.0 - tol) { restore_complemented(complemented_variables); return -1; } i_t cover_size = 0; for (i_t k = 0; k < solution.size(); k++) { @@ -909,30 +909,373 @@ i_t knapsack_generation_t::generate_knapsack_cut( } cut.rhs = -cover_size + 1; - for (i_t k = 0; k < cut.size(); k++) { - const i_t j = cut.index(k); - if (is_complemented[j]) { - cut.vector.x[k] *= -1.0; - cut.rhs += 1.0; - } - } - cut.sort(); - // The cut is in the form: - sum_{j in cover} x_j >= -cover_size + 1 // Which is equivalent to: sum_{j in cover} x_j <= cover_size - 1 // Verify the cut is violated - f_t dot = cut.vector.dot(xstar); + f_t dot = cut.vector.dot(complemented_xstar_); f_t violation = dot - cut.rhs; - if (verbose) { + if (1) { settings.log.printf("Knapsack cut %d violation %e < 0\n", knapsack_row, violation); } - restore_complemented(); - if (violation >= -tol) { return -1; } + if (violation >= -tol) { restore_complemented(complemented_variables); return -1; } + + // Compute the minimal cover and partition the variables into C1 and C2 + inequality_t minimal_cover_cut(lp.num_cols); + std::vector c1_partition; + std::vector c2_partition; + minimial_cover_and_partition(knapsack_inequality, + cut, + complemented_xstar_, + minimal_cover_cut, + c1_partition, + c2_partition); + + // Lift the cut + inequality_t lifted_cut(lp.num_cols); + lift_knapsack_cut(knapsack_inequality, + minimal_cover_cut, + c1_partition, + c2_partition, + lifted_cut); + + lifted_cut.negate(); + + + for (i_t k = 0; k < lifted_cut.size(); k++) { + const i_t j = lifted_cut.index(k); + if (is_complemented_[j]) { + lifted_cut.vector.x[k] *= -1.0; + lifted_cut.rhs += 1.0; + } + } + lifted_cut.sort(); + + // Verify the cut is violated + f_t lifted_dot = lifted_cut.vector.dot(xstar); + f_t lifted_violation = lifted_dot - lifted_cut.rhs; + if (1) { + settings.log.printf("Knapsack cut %d lifted violation %e < 0\n", knapsack_row, lifted_violation); + } + + f_t best_violation; + if (lifted_violation < violation) { + cut = lifted_cut; + best_violation = lifted_violation; + } else { + for (i_t k = 0; k < cut.size(); k++) { + const i_t j = cut.index(k); + if (is_complemented_[j]) { + cut.vector.x[k] *= -1.0; + cut.rhs += 1.0; + } + } + cut.sort(); + best_violation = violation; + } + + if (best_violation >= -tol) { + restore_complemented(complemented_variables); + return -1; + } + + restore_complemented(complemented_variables); return 0; } +template +bool knapsack_generation_t::is_minimal_cover(f_t cover_sum, + f_t beta, + const std::vector& cover_coefficients) +{ + // Check if the cover is minimial + // A set C is a cover if + // sum_{j in C} a_j > beta + // A set C is a minimal cover if + // sum_{k in C \ {j}} a_k <= beta for all j in C + bool minimal = true; + + // cover_sum = sum_{j in C} a_j + + // A cover is minimal if cover_sum - a_j <= beta for all j in C + + for (i_t k = 0; k < cover_coefficients.size(); k++) { + const f_t a_j = cover_coefficients[k]; + if (a_j == 0.0) { continue; } + if (cover_sum - a_j > beta) { + minimal = false; + break; + } + } + return minimal; + +} + +template +void knapsack_generation_t::minimial_cover_and_partition( + const inequality_t& knapsack_inequality, + const inequality_t& negated_base_cut, + const std::vector& xstar, + inequality_t& minimal_cover_cut, + std::vector& c1_partition, + std::vector& c2_partition) +{ + // Compute the minimal cover cut + inequality_t base_cut = negated_base_cut; + base_cut.negate(); + + std::vector cover_indicies; + cover_indicies.reserve(base_cut.size()); + + std::vector cover_coefficients; + cover_coefficients.reserve(base_cut.size()); + + std::vector score; + score.reserve(base_cut.size()); + + for (i_t k = 0; k < knapsack_inequality.size(); k++) { + const i_t j = knapsack_inequality.index(k); + workspace_[j] = knapsack_inequality.coeff(k); + } + + for (i_t k = 0; k < base_cut.size(); k++) { + const i_t j = base_cut.index(k); + const f_t xstar_j = xstar[j]; + score.push_back((1.0 - xstar_j) / workspace_[j]); + cover_indicies.push_back(j); + cover_coefficients.push_back(workspace_[j]); + } + + f_t cover_sum = std::accumulate(cover_coefficients.begin(), cover_coefficients.end(), 0.0); + + bool is_minimal = is_minimal_cover(cover_sum, knapsack_inequality.rhs, cover_coefficients); + + if (is_minimal) { + minimal_cover_cut = base_cut; + return; + } + + // We don't have a minimal cover. So sort the score from smallest to largest breaking ties by largest to smallest a_j + std::vector permuation(cover_indicies.size()); + std::iota(permuation.begin(), permuation.end(), 0); + std::sort(permuation.begin(), permuation.end(), [&](i_t a, i_t b) { + if (score[a] < score[b]) { + return true; + } else if (score[a] == score[b]) { + return cover_coefficients[a] > cover_coefficients[b]; + } else { + return false; + } + }); + + const f_t beta = knapsack_inequality.rhs; + for (i_t k = 0; k < permuation.size(); k++) { + const i_t h = permuation[k]; + const f_t a_j = cover_coefficients[h]; + if (cover_sum - a_j > beta) { + // sum_{k in C} a_k - a_j > beta + // so sum_{k in C \ {k}} a_k > beta and C \ {k} remains a cover + + cover_sum -= a_j; + // Set the coefficient to 0 to remove it from the cover + cover_coefficients[h] = 0.0; + + is_minimal = is_minimal_cover(cover_sum, beta, cover_coefficients); + if (is_minimal) { + break; + } + } else { + // C \ {j} is no longer a cover. + continue; + } + } + + // Go through and correct cover_indicies and cover_coefficients + for (i_t k = 0; k < cover_coefficients.size(); ) { + if (cover_coefficients[k] == 0.0) { + cover_indicies[k] = cover_indicies.back(); + cover_indicies.pop_back(); + cover_coefficients[k] = cover_coefficients.back(); + cover_coefficients.pop_back(); + } else { + k++; + } + } + + // We now have a minimal cover cut + // sum_{j in C} x_j <= |C| - 1 + minimal_cover_cut.vector.i = cover_indicies; + minimal_cover_cut.vector.x.resize(cover_indicies.size(), 1.0); + minimal_cover_cut.rhs = cover_coefficients.size() - 1; + + // Now we need to partition the variables into C1 and C2 + // C2 = {j in C | x_j = 1} + // C1 = C / C2 + + const f_t x_tol = 1e-5; + for (i_t j : cover_indicies) { + if (xstar[j] > 1.0 - x_tol) { + c2_partition.push_back(j); + } else { + c1_partition.push_back(j); + } + } +} + +template +void knapsack_generation_t::lift_knapsack_cut( + const inequality_t& knapsack_inequality, + const inequality_t& base_cut, + const std::vector& c1_partition, + const std::vector& c2_partition, + inequality_t& lifted_cut) +{ + // The base cut is in the form: sum_{j in cover} x_j <= |cover| - 1 + + // We will attempt to lift the cut by adding a new variable x_k with k not in C to the base cut + // so that the cut becomes + // sum_{j in cover} x_j + alpha_k * x_k <= |cover| - 1 + // + // We can do this for multiple variables so that in the end the cut becomes + // + // sum_{j in cover} x_j + sum_{k in F} alpha_k * x_k <= |cover| - 1 + + // Determine which variables are in the knapsack constraint and not in the cover + std::vector marked_variables; + marked_variables.reserve(knapsack_inequality.size()); + for (i_t k = 0; k < knapsack_inequality.size(); k++) { + const i_t j = knapsack_inequality.index(k); + if (!is_marked_[j]) { + is_marked_[j] = 1; // is_marked_[j] = 1 for all j in N + marked_variables.push_back(j); + } + } + for (i_t k = 0; k < base_cut.size(); k++) { + const i_t j = base_cut.index(k); + if (is_marked_[j]) { + is_marked_[j] = 0; // is_marked_[j] = 1 for all j in N \ C + // OK to leave marked_variables unchanged as marked_variables will be a superset of all dirty is_marked + } + } + std::vector remaining_variables; + std::vector remaining_indices; + std::vector remaining_coefficients; + remaining_variables.reserve(knapsack_inequality.size()); + remaining_indices.reserve(knapsack_inequality.size()); + remaining_coefficients.reserve(knapsack_inequality.size()); + + for (i_t k = 0; k < knapsack_inequality.size(); k++) { + const i_t j = knapsack_inequality.index(k); + if (is_marked_[j]) { + remaining_variables.push_back(j); + remaining_indices.push_back(k); + remaining_coefficients.push_back(knapsack_inequality.coeff(k)); + } + } + + // We start with F = {} and lift remaining variables one by one + // For a variable k not in C union F, the inequality + // + // alpha_k * x_k + sum_{j in C} x_j <= |C| - 1 + // is trivially satisfied when x_k = 0. + // If x_k = 1, then the inequality will be valid for all alpha_k + // such that + // alpha_k + maximize sum_{j in C} x_j <= |C| - 1 + // subject to a_k + sum_{j in C} a_j x_j <= beta + // + // where here we require a_k + sum_{j in C} a_j x_j <= beta so that the inequality + // is valid for the set { x_j in {0, 1}^(|C| + 1) | sum_{j in C union k} a_j x_j <= beta} + // + // Let phi^star_k denote the optimal objective value of the problem + // + // maximize sum_{j in C} x_j + // subject to a_k + sum_{j in C} a_j x_j <= beta + // x_j in {0, 1} for all j in C + // Then alpha_k <= |C| - 1 - phi^star_k + // and we can set alpha_k = |C| - 1 - phi^star_k + // + // We can continue this process for each variable k not in C union F + // + // Assume the valid inequality + // sum_{j in C} x_j + sum_{j in F} alpha_j * x_j <= |C| - 1 + // has been obtained so far. We now add the variable x_k with k not in C union F to the inequality. + // So we have + // alpha_k * x_k + sum_{j in C} x_j + sum_{j in F} alpha_j * x_j <= |C| - 1 + // + // Again, this is trivially satisfied when x_k = 0. And we can determine the max value of alpha_k + // by solving the 0-1 knapsack problem: + // + // maximize sum_{j in C} x_j + sum_{j in F} alpha_j * x_j + // subject to sum_{j in C} a_j x_j + sum_{j in F} a_j * x_j <= beta - a_k + // x_j in {0, 1} for all j in C union F + // + // Let phi^star_k denote the optimal objective value of the knapsack problem. + // The lifted coefficient alpha_k = |C| - 1 - phi^star_k + + + // Construct weight and values for C + std::vector values; + values.reserve(knapsack_inequality.size()); + + std::vector weights; + weights.reserve(knapsack_inequality.size()); + + for (i_t k = 0; k < knapsack_inequality.size(); k++) { + const i_t j = knapsack_inequality.index(k); + if (!is_marked_[j]) { + // j is in C + weights.push_back(knapsack_inequality.coeff(k)); + values.push_back(1.0); + } + } + + std::vector F; + std::vector alpha; + + std::vector solution; + + f_t cover_size = base_cut.size(); + + lifted_cut = base_cut; + + // TODO sort the remaining variables / coefficients according to some ordering + + while (remaining_variables.size() > 0) { + const i_t k = remaining_variables.back(); + const f_t a_k = remaining_coefficients.back(); + + f_t capacity = knapsack_inequality.rhs - a_k; + + f_t objective = solve_knapsack_problem(values, weights, capacity, solution); + if (std::isnan(objective)) { break; } + + f_t alpha_k = std::max(0.0, cover_size - 1.0 - objective); + + if (alpha_k > 0.0) { + printf("Lifted variable %d with alpha %g\n", k, alpha_k); + F.push_back(k); + alpha.push_back(alpha_k); + values.push_back(alpha_k); + weights.push_back(a_k); + + lifted_cut.vector.i.push_back(k); + lifted_cut.vector.x.push_back(alpha_k); + } + + // Remove the variable from the remaining variables and coefficients + remaining_variables.pop_back(); + remaining_coefficients.pop_back(); + } + printf("Lifted %ld variables\n", F.size()); + + // Restore is_marked_ + for (i_t j : marked_variables) { + is_marked_[j] = 0; + } + +} + template f_t knapsack_generation_t::greedy_knapsack_problem(const std::vector& values, const std::vector& weights, @@ -1054,6 +1397,8 @@ f_t knapsack_generation_t::solve_knapsack_problem(const std::vector dp(n + 1, sum_value + 1, INT_INF); dense_matrix_t take(n + 1, sum_value + 1, 0); @@ -1170,11 +1515,10 @@ void cut_generation_t::generate_knapsack_cuts( const std::vector& xstar) { if (knapsack_generation_.num_knapsack_constraints() > 0) { - std::vector is_complemented(lp.num_cols, 0); for (i_t knapsack_row : knapsack_generation_.get_knapsack_constraints()) { inequality_t cut(lp.num_cols); i_t knapsack_status = knapsack_generation_.generate_knapsack_cut( - lp, settings, Arow, new_slacks, var_types, xstar, knapsack_row, is_complemented, cut); + lp, settings, Arow, new_slacks, var_types, xstar, knapsack_row, cut); if (knapsack_status == 0) { cut_pool_.add_cut(cut_type_t::KNAPSACK, cut); } } } diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index b79ff1e8b..aad1a51a2 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -310,13 +310,34 @@ class knapsack_generation_t { const std::vector& var_types, const std::vector& xstar, i_t knapsack_row, - std::vector& is_complemented, inequality_t& cut); i_t num_knapsack_constraints() const { return knapsack_constraints_.size(); } const std::vector& get_knapsack_constraints() const { return knapsack_constraints_; } private: + void restore_complemented(const std::vector& complemented_variables) { + for (i_t j : complemented_variables) { + is_complemented_[j] = 0; + } + } + bool is_minimal_cover(f_t cover_sum, + f_t beta, + const std::vector& cover_coefficients); + + void minimial_cover_and_partition(const inequality_t& knapsack_inequality, + const inequality_t& negated_base_cut, + const std::vector& xstar, + inequality_t& minimal_cover_cut, + std::vector& c1_partition, + std::vector& c2_partition); + + void lift_knapsack_cut(const inequality_t& knapsack_inequality, + const inequality_t& base_cut, + const std::vector& c1_partition, + const std::vector& c2_partition, + inequality_t& lifted_cut); + // Generate a heuristic solution to the 0-1 knapsack problem f_t greedy_knapsack_problem(const std::vector& values, const std::vector& weights, @@ -331,6 +352,10 @@ class knapsack_generation_t { std::vector is_slack_; std::vector knapsack_constraints_; + std::vector is_complemented_; + std::vector is_marked_; + std::vector workspace_; + std::vector complemented_xstar_; const simplex_solver_settings_t& settings_; }; From 9b6c6153a7d8db053adc1fc039ff4651a54399e1 Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Wed, 18 Mar 2026 16:54:00 -0700 Subject: [PATCH 04/14] Minor tweaks --- cpp/src/cuts/cuts.cpp | 63 +++++++++++++------------------------------ 1 file changed, 19 insertions(+), 44 deletions(-) diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index 115fa82e6..3f49b5f2f 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -725,20 +725,19 @@ knapsack_generation_t::knapsack_generation_t( if (is_knapsack) { const f_t beta = inequality.rhs + sum_neg; - if (std::abs(beta - std::round(beta)) <= settings.integer_tol) { - if (beta > 0.0 && beta <= sum_pos && std::abs(sum_pos / (row_len - 1) - beta) > 1e-3) { - if (1) { - settings.log.printf( - "Knapsack constraint %d row len %d beta %e sum_neg %e sum_pos %e sum_pos / (row_len - 1) %e\n", - i, - row_len, - beta, - sum_neg, - sum_pos, - sum_pos / (row_len - 1)); - } - knapsack_constraints_.push_back(i); + if (beta > 0.0 && beta <= sum_pos && std::abs(sum_pos / (row_len - 1) - beta) > 1e-3) { + if (verbose) { + settings.log.printf( + "Knapsack constraint %d row len %d beta %e sum_neg %e sum_pos %e sum_pos / (row_len - " + "1) %e\n", + i, + row_len, + beta, + sum_neg, + sum_pos, + sum_pos / (row_len - 1)); } + knapsack_constraints_.push_back(i); } } } @@ -912,14 +911,6 @@ i_t knapsack_generation_t::generate_knapsack_cut( // The cut is in the form: - sum_{j in cover} x_j >= -cover_size + 1 // Which is equivalent to: sum_{j in cover} x_j <= cover_size - 1 - // Verify the cut is violated - f_t dot = cut.vector.dot(complemented_xstar_); - f_t violation = dot - cut.rhs; - if (1) { - settings.log.printf("Knapsack cut %d violation %e < 0\n", knapsack_row, violation); - } - - if (violation >= -tol) { restore_complemented(complemented_variables); return -1; } // Compute the minimal cover and partition the variables into C1 and C2 inequality_t minimal_cover_cut(lp.num_cols); @@ -939,7 +930,6 @@ i_t knapsack_generation_t::generate_knapsack_cut( c1_partition, c2_partition, lifted_cut); - lifted_cut.negate(); @@ -955,31 +945,16 @@ i_t knapsack_generation_t::generate_knapsack_cut( // Verify the cut is violated f_t lifted_dot = lifted_cut.vector.dot(xstar); f_t lifted_violation = lifted_dot - lifted_cut.rhs; - if (1) { + if (verbose) { settings.log.printf("Knapsack cut %d lifted violation %e < 0\n", knapsack_row, lifted_violation); } - f_t best_violation; - if (lifted_violation < violation) { - cut = lifted_cut; - best_violation = lifted_violation; - } else { - for (i_t k = 0; k < cut.size(); k++) { - const i_t j = cut.index(k); - if (is_complemented_[j]) { - cut.vector.x[k] *= -1.0; - cut.rhs += 1.0; - } - } - cut.sort(); - best_violation = violation; - } - - if (best_violation >= -tol) { + if (lifted_violation >= -tol) { restore_complemented(complemented_variables); return -1; } + cut = lifted_cut; restore_complemented(complemented_variables); return 0; } @@ -1253,7 +1228,7 @@ void knapsack_generation_t::lift_knapsack_cut( f_t alpha_k = std::max(0.0, cover_size - 1.0 - objective); if (alpha_k > 0.0) { - printf("Lifted variable %d with alpha %g\n", k, alpha_k); + settings_.log.printf("Lifted variable %d with alpha %g\n", k, alpha_k); F.push_back(k); alpha.push_back(alpha_k); values.push_back(alpha_k); @@ -1267,7 +1242,7 @@ void knapsack_generation_t::lift_knapsack_cut( remaining_variables.pop_back(); remaining_coefficients.pop_back(); } - printf("Lifted %ld variables\n", F.size()); + settings_.log.printf("Lifted %ld variables\n", F.size()); // Restore is_marked_ for (i_t j : marked_variables) { @@ -1474,8 +1449,8 @@ bool cut_generation_t::generate_cuts(const lp_problem_t& lp, f_t cut_start_time = tic(); generate_knapsack_cuts(lp, settings, Arow, new_slacks, var_types, xstar); f_t cut_generation_time = toc(cut_start_time); - if (cut_generation_time > 1.0) { - settings.log.debug("Knapsack cut generation time %.2f seconds\n", cut_generation_time); + if (1 || cut_generation_time > 1.0) { + settings.log.printf("Knapsack cut generation time %.2f seconds\n", cut_generation_time); } } From 3709bf9e2656ba2fcfdaf3f1fbb9dd8ef04b0008 Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Thu, 19 Mar 2026 10:09:47 -0700 Subject: [PATCH 05/14] Fix bugs found by coderabbit --- cpp/src/cuts/cuts.cpp | 34 +++++++++++++++++++++++++++++++--- cpp/src/cuts/cuts.hpp | 15 --------------- 2 files changed, 31 insertions(+), 18 deletions(-) diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index 3f49b5f2f..22392484e 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -445,6 +445,21 @@ void extend_clique_vertices(std::vector& clique_vertices, } // namespace +template +bool rational_coefficients(const std::vector& var_types, + const inequality_t& inequality, + inequality_t& rational_inequality); + +template +bool rational_approximation(f_t x, + int64_t max_denominator, + int64_t& numerator, + int64_t& denominator); + +int64_t gcd(const std::vector& integers); + +int64_t lcm(const std::vector& integers); + // This function is only used in tests std::vector> find_maximal_cliques_for_test( const std::vector>& adjacency_list, @@ -876,8 +891,14 @@ i_t knapsack_generation_t::generate_knapsack_cut( return -1; } - if (verbose) { settings.log.printf("Calling solve_knapsack_problem\n"); } - f_t objective = solve_knapsack_problem(values, weights, seperation_rhs, solution); + f_t objective = 0.0; + if (!values.empty()) { + if (verbose) { settings.log.printf("Calling solve_knapsack_problem\n"); } + + objective = solve_knapsack_problem(values, weights, seperation_rhs, solution); + } else { + solution.clear(); + } if (std::isnan(objective)) { restore_complemented(complemented_variables); return -1; } if (verbose) { settings.log.printf("objective %e objective_constant %e\n", objective, objective_constant); @@ -933,11 +954,18 @@ i_t knapsack_generation_t::generate_knapsack_cut( lifted_cut.negate(); + // The cut is now in the form: + // -\sum_{j in C} x_j - \sum_{j in F} alpha_j x_j >= -cover_size + 1 for (i_t k = 0; k < lifted_cut.size(); k++) { const i_t j = lifted_cut.index(k); + // \sum_{k!=j} d_k x_k + d_j xbar_j >= gamma + // xbar_j = 1 - x_j + // \sum_{k!=j} d_k x_k + d_j (1 - x_j) >= gamma + // \sum_{k!=j} d_k x_k + d_j - d_j x_j >= gamma + // \sum_{k!=j} d_k x_k + d_j x_j >= gamma - d_j if (is_complemented_[j]) { + lifted_cut.rhs -= lifted_cut.vector.x[k]; lifted_cut.vector.x[k] *= -1.0; - lifted_cut.rhs += 1.0; } } lifted_cut.sort(); diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index aad1a51a2..81a2bc2e2 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -166,21 +166,6 @@ f_t fractional_part(f_t a) return a - std::floor(a); } -template -bool rational_coefficients(const std::vector& var_types, - const inequality_t& inequality, - inequality_t& rational_inequality); - -template -bool rational_approximation(f_t x, - int64_t max_denominator, - int64_t& numerator, - int64_t& denominator); - -int64_t gcd(const std::vector& integers); - -int64_t lcm(const std::vector& integers); - template bool add_work_estimate(f_t accesses, f_t* work_estimate, From 9c54e01605c2d36e4574c0368ce313d16b282e3c Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Thu, 19 Mar 2026 11:10:51 -0700 Subject: [PATCH 06/14] Use greedy approximation in lifting. Sort lifting variables --- cpp/src/cuts/cuts.cpp | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index 22392484e..a7469ecb5 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -1242,16 +1242,20 @@ void knapsack_generation_t::lift_knapsack_cut( lifted_cut = base_cut; - // TODO sort the remaining variables / coefficients according to some ordering + // Sort the coefficients such that the largest coefficients are lifted first + // We will pop the largest coefficients from the back of the permutation + std::vector permutation; + best_score_last_permutation(remaining_coefficients, permutation); - while (remaining_variables.size() > 0) { - const i_t k = remaining_variables.back(); - const f_t a_k = remaining_coefficients.back(); + while (permutation.size() > 0) { + const i_t h = permutation.back(); + const i_t k = remaining_variables[h]; + const f_t a_k = remaining_coefficients[h]; f_t capacity = knapsack_inequality.rhs - a_k; - f_t objective = solve_knapsack_problem(values, weights, capacity, solution); - if (std::isnan(objective)) { break; } + f_t objective = greedy_knapsack_problem(values, weights, capacity, solution); + if (std::isnan(objective)) { settings_.log.printf("lifting knapsack problem failed\n"); break; } f_t alpha_k = std::max(0.0, cover_size - 1.0 - objective); @@ -1266,9 +1270,8 @@ void knapsack_generation_t::lift_knapsack_cut( lifted_cut.vector.x.push_back(alpha_k); } - // Remove the variable from the remaining variables and coefficients - remaining_variables.pop_back(); - remaining_coefficients.pop_back(); + // Remove the variable from the permutation + permutation.pop_back(); } settings_.log.printf("Lifted %ld variables\n", F.size()); From 7b60c084e06aaa649fad09ceccfc60a6157901a8 Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Thu, 19 Mar 2026 12:49:06 -0700 Subject: [PATCH 07/14] Fix typo --- cpp/src/cuts/cuts.cpp | 10 +++------- cpp/src/cuts/cuts.hpp | 12 ++++++------ 2 files changed, 9 insertions(+), 13 deletions(-) diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index a7469ecb5..3602e3832 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -937,12 +937,8 @@ i_t knapsack_generation_t::generate_knapsack_cut( inequality_t minimal_cover_cut(lp.num_cols); std::vector c1_partition; std::vector c2_partition; - minimial_cover_and_partition(knapsack_inequality, - cut, - complemented_xstar_, - minimal_cover_cut, - c1_partition, - c2_partition); + minimal_cover_and_partition( + knapsack_inequality, cut, complemented_xstar_, minimal_cover_cut, c1_partition, c2_partition); // Lift the cut inequality_t lifted_cut(lp.num_cols); @@ -1016,7 +1012,7 @@ bool knapsack_generation_t::is_minimal_cover(f_t cover_sum, } template -void knapsack_generation_t::minimial_cover_and_partition( +void knapsack_generation_t::minimal_cover_and_partition( const inequality_t& knapsack_inequality, const inequality_t& negated_base_cut, const std::vector& xstar, diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index 81a2bc2e2..bf9aad172 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -310,12 +310,12 @@ class knapsack_generation_t { f_t beta, const std::vector& cover_coefficients); - void minimial_cover_and_partition(const inequality_t& knapsack_inequality, - const inequality_t& negated_base_cut, - const std::vector& xstar, - inequality_t& minimal_cover_cut, - std::vector& c1_partition, - std::vector& c2_partition); + void minimal_cover_and_partition(const inequality_t& knapsack_inequality, + const inequality_t& negated_base_cut, + const std::vector& xstar, + inequality_t& minimal_cover_cut, + std::vector& c1_partition, + std::vector& c2_partition); void lift_knapsack_cut(const inequality_t& knapsack_inequality, const inequality_t& base_cut, From 6d3b5c42b03bc0243d4798267291de328d29aac0 Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Thu, 19 Mar 2026 19:11:30 -0700 Subject: [PATCH 08/14] Don't zero out score of rows we have tried generating an MIR cut from. This allows rows we have examined to be choosen in aggregation. Helps on dcmulti and a1c1s1. --- cpp/src/cuts/cuts.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index 3602e3832..e45d2a0dd 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -1836,6 +1836,8 @@ void cut_generation_t::generate_mir_cuts( // Transform the relaxation solution std::vector transformed_xstar; complemented_mir.bound_substitution(lp, variable_bounds, var_types, xstar, transformed_xstar); + std::vector initial_scores = scores; + std::vector starting_count(lp.num_rows, 0); const i_t max_cuts = std::min(lp.num_rows, 100000); f_t work_estimate = 0.0; @@ -2069,8 +2071,8 @@ void cut_generation_t::generate_mir_cuts( // Clear the aggregated rows aggregated_rows.clear(); - // Set the score of the current row to zero - scores[i] = 0.0; + starting_count[i]++; + scores[i] = initial_scores[i] * std::pow(0.9, starting_count[i]); score_queue.push(std::make_pair(scores[i], i)); work_estimate += std::log2(std::max(1, static_cast(score_queue.size()))); } From fec6d17db654cdd5ca2d84683ec5ee21664dc11b Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Thu, 19 Mar 2026 22:00:42 -0700 Subject: [PATCH 09/14] Add implied bound cuts --- cpp/src/branch_and_bound/branch_and_bound.cpp | 3 + cpp/src/branch_and_bound/branch_and_bound.hpp | 2 + cpp/src/cuts/cuts.cpp | 125 ++++++++++++++++++ cpp/src/cuts/cuts.hpp | 56 +++++++- .../dual_simplex/simplex_solver_settings.hpp | 2 + cpp/src/dual_simplex/solve.cpp | 6 +- cpp/src/mip_heuristics/diversity/lns/rins.cu | 5 +- .../diversity/recombiners/sub_mip.cuh | 4 +- cpp/src/mip_heuristics/solver.cu | 74 +++++++++++ 9 files changed, 271 insertions(+), 6 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 2133f156f..1345cb3e7 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -243,9 +243,11 @@ branch_and_bound_t::branch_and_bound_t( const user_problem_t& user_problem, const simplex_solver_settings_t& solver_settings, f_t start_time, + const probing_implied_bounds_t& probing_implied_bounds, std::shared_ptr> clique_table) : original_problem_(user_problem), settings_(solver_settings), + probing_implied_bounds_(probing_implied_bounds), clique_table_(std::move(clique_table)), original_lp_(user_problem.handle_ptr, 1, 1, 1), Arow_(1, 1, 0), @@ -2144,6 +2146,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut new_slacks_, var_types_, original_problem_, + probing_implied_bounds_, clique_table_, &clique_table_future_, &signal_extend_cliques_); diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index 98674b7f9..3c5ee6538 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -77,6 +77,7 @@ class branch_and_bound_t { branch_and_bound_t(const user_problem_t& user_problem, const simplex_solver_settings_t& solver_settings, f_t start_time, + const probing_implied_bounds_t& probing_implied_bounds, std::shared_ptr> clique_table = nullptr); // Set an initial guess based on the user_problem. This should be called before solve. @@ -148,6 +149,7 @@ class branch_and_bound_t { private: const user_problem_t& original_problem_; const simplex_solver_settings_t settings_; + const probing_implied_bounds_t& probing_implied_bounds_; std::shared_ptr> clique_table_; std::future>> clique_table_future_; std::atomic signal_extend_cliques_{false}; diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index e45d2a0dd..c00341ed3 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -1445,6 +1445,121 @@ f_t knapsack_generation_t::solve_knapsack_problem(const std::vector +void cut_generation_t::generate_implied_bounds_cuts( + const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + const std::vector& var_types, + const std::vector& xstar, + variable_bounds_t& variable_bounds, + f_t start_time) +{ + if (probing_implied_bounds_.zero_offsets.empty()) { return; } + + const f_t tol = 1e-4; + i_t num_cuts = 0; + const i_t pib_cols = static_cast(probing_implied_bounds_.zero_offsets.size()) - 1; + const i_t n_cols = std::min(lp.num_cols, pib_cols); + + for (i_t j = 0; j < n_cols; j++) { + if (var_types[j] == variable_type_t::CONTINUOUS) { continue; } + const f_t xstar_j = xstar[j]; + + // x_j = 0 implications + const i_t zero_begin = probing_implied_bounds_.zero_offsets[j]; + const i_t zero_end = probing_implied_bounds_.zero_offsets[j + 1]; + for (i_t p = zero_begin; p < zero_end; p++) { + const i_t i = probing_implied_bounds_.zero_variables[p]; + const f_t l_i = lp.lower[i]; + const f_t u_i = lp.upper[i]; + + // Tightened upper bound: x_j = 0 implies y_i <= b, where b < u_i + // Valid inequality: y_i <= b + (u_i - b)*x_j or -y_i + (u_i - b)*x_j >= -b + const f_t b_ub = probing_implied_bounds_.zero_upper_bound[p]; + if (b_ub < u_i - tol) { + const f_t coeff_j = u_i - b_ub; + const f_t y_i = xstar[i]; + const f_t lhs_val = -y_i + coeff_j * xstar_j; + const f_t rhs_val = -b_ub; + if (lhs_val < rhs_val - tol) { + inequality_t cut; + cut.push_back(i, -1.0); + cut.push_back(j, coeff_j); + cut.rhs = -b_ub; + cut_pool_.add_cut(cut_type_t::IMPLIED_BOUNDS, cut); + num_cuts++; + } + } + + // Tightened lower bound: x_j = 0 implies y_i >= b, where b > l_i + // Valid inequality: y_i >= b - (b - l_i)*x_j or y_i + (b - l_i)*x_j >= b + const f_t b_lb = probing_implied_bounds_.zero_lower_bound[p]; + if (b_lb > l_i + tol) { + const f_t coeff_j = b_lb - l_i; + const f_t y_i = xstar[i]; + const f_t lhs_val = y_i + coeff_j * xstar_j; + const f_t rhs_val = b_lb; + if (lhs_val < rhs_val - tol) { + inequality_t cut; + cut.push_back(i, 1.0); + cut.push_back(j, coeff_j); + cut.rhs = b_lb; + cut_pool_.add_cut(cut_type_t::IMPLIED_BOUNDS, cut); + num_cuts++; + } + } + } + + // x_j = 1 implications + const i_t one_begin = probing_implied_bounds_.one_offsets[j]; + const i_t one_end = probing_implied_bounds_.one_offsets[j + 1]; + for (i_t p = one_begin; p < one_end; p++) { + const i_t i = probing_implied_bounds_.one_variables[p]; + const f_t l_i = lp.lower[i]; + const f_t u_i = lp.upper[i]; + + // Tightened upper bound: x_j = 1 implies y_i <= b, where b < u_i + // Valid inequality: y_i <= u_i - (u_i - b)*x_j or -y_i - (u_i - b)*x_j >= -u_i + const f_t b_ub = probing_implied_bounds_.one_upper_bound[p]; + if (b_ub < u_i - tol) { + const f_t coeff_j = -(u_i - b_ub); + const f_t y_i = xstar[i]; + const f_t lhs_val = -y_i + coeff_j * xstar_j; + const f_t rhs_val = -u_i; + if (lhs_val < rhs_val - tol) { + inequality_t cut; + cut.push_back(i, -1.0); + cut.push_back(j, coeff_j); + cut.rhs = -u_i; + cut_pool_.add_cut(cut_type_t::IMPLIED_BOUNDS, cut); + num_cuts++; + } + } + + // Tightened lower bound: x_j = 1 implies y_i >= b, where b > l_i + // Valid inequality: y_i >= l_i + (b - l_i)*x_j or y_i - (b - l_i)*x_j >= l_i + const f_t b_lb = probing_implied_bounds_.one_lower_bound[p]; + if (b_lb > l_i + tol) { + const f_t coeff_j = -(b_lb - l_i); + const f_t lhs_val = xstar[i] + coeff_j * xstar_j; + const f_t rhs_val = l_i; + if (lhs_val < rhs_val - tol) { + inequality_t cut; + cut.push_back(i, 1.0); + cut.push_back(j, coeff_j); + cut.rhs = rhs_val; + cut_pool_.add_cut(cut_type_t::IMPLIED_BOUNDS, cut); + num_cuts++; + } + } + } + } + + if (num_cuts > 0) { + settings.log.debug("Generated %d implied bounds cuts from probing\n", num_cuts); + } +} + template bool cut_generation_t::generate_cuts(const lp_problem_t& lp, const simplex_solver_settings_t& settings, @@ -1504,6 +1619,16 @@ bool cut_generation_t::generate_cuts(const lp_problem_t& lp, settings.log.debug("Clique cut generation time %.2f seconds\n", cut_generation_time); } } + + // Generate implied bounds cuts + if (settings.implied_bounds_cuts != 0) { + f_t cut_start_time = tic(); + generate_implied_bounds_cuts(lp, settings, var_types, xstar, variable_bounds, start_time); + f_t cut_generation_time = toc(cut_start_time); + if (cut_generation_time > 1.0) { + settings.log.debug("Implied bounds cut generation time %.2f seconds\n", cut_generation_time); + } + } return true; } diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index bf9aad172..2653ee793 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -38,7 +38,8 @@ enum cut_type_t : int8_t { KNAPSACK = 2, CHVATAL_GOMORY = 3, CLIQUE = 4, - MAX_CUT_TYPE = 5 + IMPLIED_BOUNDS = 5, + MAX_CUT_TYPE = 6 }; template @@ -62,6 +63,46 @@ cut_gap_closure_t compute_cut_gap_closure(f_t objective_reference, return {initial_gap, final_gap, gap_closed, gap_closed_ratio}; } +template +struct probing_implied_bounds_t { + // Probing implications stored in CSR format, indexed by binary variable x_j. + // + // "zero" = implications discovered when probing x_j = 0. + // "one" = implications discovered when probing x_j = 1. + // + // For a binary variable x_j, the range + // zero_offsets[j] .. zero_offsets[j+1] + // indexes into the flat arrays zero_variables, zero_lower_bound, zero_upper_bound. + // + // For each position p in that range: + // zero_variables[p] = i if variable y_i bounds were tightened + // when x_j was fixed to 0 and constraints were propagated. + // zero_lower_bound[p] = tightened lower bound on y_i (i.e., x_j = 0 => y_i >= zero_lower_bound[p]). + // zero_upper_bound[p] = tightened upper bound on y_i (i.e., x_j = 0 => y_i <= zero_upper_bound[p]). + // + // The one arrays are analogous for probing x_j = 1. + // + // Non-binary variables have empty ranges (zero_offsets[j] == zero_offsets[j+1]). + // Offsets vectors have size num_cols + 1. + + probing_implied_bounds_t() = default; + + probing_implied_bounds_t(i_t num_cols) + : zero_offsets(num_cols + 1, 0), one_offsets(num_cols + 1, 0) + { + } + + std::vector zero_offsets; + std::vector zero_variables; + std::vector zero_lower_bound; + std::vector zero_upper_bound; + + std::vector one_offsets; + std::vector one_variables; + std::vector one_lower_bound; + std::vector one_upper_bound; +}; + template struct inequality_t { inequality_t() : vector(), rhs(0.0) {} @@ -132,7 +173,7 @@ struct cut_info_t { } } const char* cut_type_names[MAX_CUT_TYPE] = { - "Gomory ", "MIR ", "Knapsack ", "Strong CG", "Clique "}; + "Gomory ", "MIR ", "Knapsack ", "Strong CG", "Clique ", "Implied Bounds"}; std::array num_cuts = {0}; }; @@ -362,12 +403,14 @@ class cut_generation_t { const std::vector& new_slacks, const std::vector& var_types, const user_problem_t& user_problem, + const probing_implied_bounds_t& probing_implied_bounds, std::shared_ptr> clique_table = nullptr, std::future>>* clique_table_future = nullptr, std::atomic* signal_extend = nullptr) : cut_pool_(cut_pool), knapsack_generation_(lp, settings, Arow, new_slacks, var_types), user_problem_(user_problem), + probing_implied_bounds_(probing_implied_bounds), clique_table_(std::move(clique_table)), clique_table_future_(clique_table_future), signal_extend_(signal_extend) @@ -426,9 +469,18 @@ class cut_generation_t { const std::vector& reduced_costs, f_t start_time); + // Generate implied bounds cuts from probing implications + void generate_implied_bounds_cuts(const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + const std::vector& var_types, + const std::vector& xstar, + variable_bounds_t& variable_bounds, + f_t start_time); + cut_pool_t& cut_pool_; knapsack_generation_t knapsack_generation_; const user_problem_t& user_problem_; + const probing_implied_bounds_t& probing_implied_bounds_; std::shared_ptr> clique_table_; std::future>>* clique_table_future_{nullptr}; std::atomic* signal_extend_{nullptr}; diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index eadd93040..67e52f9f7 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -100,6 +100,7 @@ struct simplex_solver_settings_t { mir_cuts(-1), mixed_integer_gomory_cuts(-1), knapsack_cuts(-1), + implied_bounds_cuts(-1), clique_cuts(-1), strong_chvatal_gomory_cuts(-1), reduced_cost_strengthening(-1), @@ -180,6 +181,7 @@ struct simplex_solver_settings_t { i_t mixed_integer_gomory_cuts; // -1 automatic, 0 to disable, >0 to enable mixed integer Gomory // cuts i_t knapsack_cuts; // -1 automatic, 0 to disable, >0 to enable knapsack cuts + i_t implied_bounds_cuts; // -1 automatic, 0 to disable, >0 to enable implied bounds cuts i_t clique_cuts; // -1 automatic, 0 to disable, >0 to enable clique cuts i_t strong_chvatal_gomory_cuts; // -1 automatic, 0 to disable, >0 to enable strong Chvatal Gomory // cuts diff --git a/cpp/src/dual_simplex/solve.cpp b/cpp/src/dual_simplex/solve.cpp index d300d6011..c0f919020 100644 --- a/cpp/src/dual_simplex/solve.cpp +++ b/cpp/src/dual_simplex/solve.cpp @@ -706,7 +706,8 @@ i_t solve(const user_problem_t& problem, { i_t status; if (is_mip(problem) && !settings.relaxation) { - branch_and_bound_t branch_and_bound(problem, settings, tic()); + probing_implied_bounds_t empty_probing(problem.num_cols); + branch_and_bound_t branch_and_bound(problem, settings, tic(), empty_probing); mip_solution_t mip_solution(problem.num_cols); mip_status_t mip_status = branch_and_bound.solve(mip_solution); if (mip_status == mip_status_t::OPTIMAL) { @@ -745,7 +746,8 @@ i_t solve_mip_with_guess(const user_problem_t& problem, { i_t status; if (is_mip(problem)) { - branch_and_bound_t branch_and_bound(problem, settings, tic()); + probing_implied_bounds_t empty_probing(problem.num_cols); + branch_and_bound_t branch_and_bound(problem, settings, tic(), empty_probing); branch_and_bound.set_initial_guess(guess); mip_status_t mip_status = branch_and_bound.solve(solution); if (mip_status == mip_status_t::OPTIMAL) { diff --git a/cpp/src/mip_heuristics/diversity/lns/rins.cu b/cpp/src/mip_heuristics/diversity/lns/rins.cu index d7d760101..de33d5487 100644 --- a/cpp/src/mip_heuristics/diversity/lns/rins.cu +++ b/cpp/src/mip_heuristics/diversity/lns/rins.cu @@ -22,6 +22,7 @@ #include #include +#include #include namespace cuopt::linear_programming::detail { @@ -274,8 +275,10 @@ void rins_t::run_rins() f_t objective) { rins_solution_queue.push_back(solution); }; + dual_simplex::probing_implied_bounds_t empty_probing( + branch_and_bound_problem.num_cols); dual_simplex::branch_and_bound_t branch_and_bound( - branch_and_bound_problem, branch_and_bound_settings, dual_simplex::tic()); + branch_and_bound_problem, branch_and_bound_settings, dual_simplex::tic(), empty_probing); branch_and_bound.set_initial_guess(cuopt::host_copy(fixed_assignment, rins_handle.get_stream())); branch_and_bound_status = branch_and_bound.solve(branch_and_bound_solution); diff --git a/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh b/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh index a867141d0..42e33d214 100644 --- a/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh +++ b/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh @@ -117,8 +117,10 @@ class sub_mip_recombiner_t : public recombiner_t { // disable B&B logs, so that it is not interfering with the main B&B thread branch_and_bound_settings.log.log = false; + dual_simplex::probing_implied_bounds_t empty_probing( + branch_and_bound_problem.num_cols); dual_simplex::branch_and_bound_t branch_and_bound( - branch_and_bound_problem, branch_and_bound_settings, dual_simplex::tic()); + branch_and_bound_problem, branch_and_bound_settings, dual_simplex::tic(), empty_probing); branch_and_bound_status = branch_and_bound.solve(branch_and_bound_solution); if (solution_vector.size() > 0) { cuopt_assert(fixed_assignment.size() == branch_and_bound_solution.x.size(), diff --git a/cpp/src/mip_heuristics/solver.cu b/cpp/src/mip_heuristics/solver.cu index f25c093af..0e03a529f 100644 --- a/cpp/src/mip_heuristics/solver.cu +++ b/cpp/src/mip_heuristics/solver.cu @@ -192,6 +192,8 @@ solution_t mip_solver_t::run_solver() branch_and_bound_solution_helper_t solution_helper(&dm, branch_and_bound_settings); dual_simplex::mip_solution_t branch_and_bound_solution(1); + dual_simplex::probing_implied_bounds_t probing_implied_bounds; + bool run_bb = !context.settings.heuristics_only; if (run_bb) { // Convert the presolved problem to dual_simplex::user_problem_t @@ -199,6 +201,77 @@ solution_t mip_solver_t::run_solver() // Resize the solution now that we know the number of columns/variables branch_and_bound_solution.resize(branch_and_bound_problem.num_cols); + // Extract probing cache into CPU-only CSR struct for implied bounds cuts + { + const i_t num_cols = branch_and_bound_problem.num_cols; + probing_implied_bounds = dual_simplex::probing_implied_bounds_t(num_cols); + auto& pc = dm.ls.constraint_prop.bounds_update.probing_cache.probing_cache; + auto& rev_ids = context.problem_ptr->reverse_original_ids; + + // First pass: count entries per binary variable + for (auto& [var_idx, entries] : pc) { + if (entries[0].val_interval.interval_type != interval_type_t::EQUALS) { continue; } + i_t j = (var_idx < static_cast(rev_ids.size())) ? rev_ids[var_idx] : -1; + if (j < 0 || j >= num_cols) { continue; } + + for (auto& [imp_var, bound] : entries[0].var_to_cached_bound_map) { + i_t i = (imp_var < static_cast(rev_ids.size())) ? rev_ids[imp_var] : -1; + if (i < 0 || i >= num_cols) { continue; } + probing_implied_bounds.zero_offsets[j + 1]++; + } + for (auto& [imp_var, bound] : entries[1].var_to_cached_bound_map) { + i_t i = (imp_var < static_cast(rev_ids.size())) ? rev_ids[imp_var] : -1; + if (i < 0 || i >= num_cols) { continue; } + probing_implied_bounds.one_offsets[j + 1]++; + } + } + + // Prefix sum + for (i_t j = 0; j < num_cols; j++) { + probing_implied_bounds.zero_offsets[j + 1] += probing_implied_bounds.zero_offsets[j]; + probing_implied_bounds.one_offsets[j + 1] += probing_implied_bounds.one_offsets[j]; + } + + // Allocate flat arrays + i_t zero_nnz = probing_implied_bounds.zero_offsets[num_cols]; + i_t one_nnz = probing_implied_bounds.one_offsets[num_cols]; + probing_implied_bounds.zero_variables.resize(zero_nnz); + probing_implied_bounds.zero_lower_bound.resize(zero_nnz); + probing_implied_bounds.zero_upper_bound.resize(zero_nnz); + probing_implied_bounds.one_variables.resize(one_nnz); + probing_implied_bounds.one_lower_bound.resize(one_nnz); + probing_implied_bounds.one_upper_bound.resize(one_nnz); + + // Second pass: fill flat arrays using write cursors + std::vector zero_cursor(probing_implied_bounds.zero_offsets); + std::vector one_cursor(probing_implied_bounds.one_offsets); + + for (auto& [var_idx, entries] : pc) { + if (entries[0].val_interval.interval_type != interval_type_t::EQUALS) { continue; } + i_t j = (var_idx < static_cast(rev_ids.size())) ? rev_ids[var_idx] : -1; + if (j < 0 || j >= num_cols) { continue; } + + for (auto& [imp_var, bound] : entries[0].var_to_cached_bound_map) { + i_t i = (imp_var < static_cast(rev_ids.size())) ? rev_ids[imp_var] : -1; + if (i < 0 || i >= num_cols) { continue; } + i_t p = zero_cursor[j]++; + probing_implied_bounds.zero_variables[p] = i; + probing_implied_bounds.zero_lower_bound[p] = bound.lb; + probing_implied_bounds.zero_upper_bound[p] = bound.ub; + } + for (auto& [imp_var, bound] : entries[1].var_to_cached_bound_map) { + i_t i = (imp_var < static_cast(rev_ids.size())) ? rev_ids[imp_var] : -1; + if (i < 0 || i >= num_cols) { continue; } + i_t p = one_cursor[j]++; + probing_implied_bounds.one_variables[p] = i; + probing_implied_bounds.one_lower_bound[p] = bound.lb; + probing_implied_bounds.one_upper_bound[p] = bound.ub; + } + } + + CUOPT_LOG_INFO("Probing implied bounds: %d zero entries, %d one entries", zero_nnz, one_nnz); + } + // Fill in the settings for branch and bound branch_and_bound_settings.time_limit = timer_.get_time_limit(); branch_and_bound_settings.node_limit = context.settings.node_limit; @@ -269,6 +342,7 @@ solution_t mip_solver_t::run_solver() branch_and_bound_problem, branch_and_bound_settings, timer_.get_tic_start(), + probing_implied_bounds, context.problem_ptr->clique_table); context.branch_and_bound_ptr = branch_and_bound.get(); auto* stats_ptr = &context.stats; From e7e690b28abc1cf567d639a9b88bebcf380072f9 Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Fri, 20 Mar 2026 19:31:35 -0700 Subject: [PATCH 10/14] Negative slack bug strikes again --- cpp/src/cuts/cuts.cpp | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index c00341ed3..162b2aba4 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -724,7 +724,13 @@ knapsack_generation_t::knapsack_generation_t( f_t sum_neg = 0.0; for (i_t p = 0; p < row_len; p++) { const i_t j = inequality.index(p); - if (is_slack_[j]) { continue; } + if (is_slack_[j]) { + if (inequality.coeff(p) < 0.0) { + is_knapsack = false; + break; + } + continue; + } const f_t aj = inequality.coeff(p); if (var_types[j] != variable_type_t::INTEGER || lp.lower[j] != 0.0 || lp.upper[j] != 1.0) { is_knapsack = false; @@ -1167,6 +1173,7 @@ void knapsack_generation_t::lift_knapsack_cut( for (i_t k = 0; k < knapsack_inequality.size(); k++) { const i_t j = knapsack_inequality.index(k); if (is_marked_[j]) { + if (is_slack_[j]) { continue; } remaining_variables.push_back(j); remaining_indices.push_back(k); remaining_coefficients.push_back(knapsack_inequality.coeff(k)); @@ -1250,7 +1257,7 @@ void knapsack_generation_t::lift_knapsack_cut( f_t capacity = knapsack_inequality.rhs - a_k; - f_t objective = greedy_knapsack_problem(values, weights, capacity, solution); + f_t objective = solve_knapsack_problem(values, weights, capacity, solution); if (std::isnan(objective)) { settings_.log.printf("lifting knapsack problem failed\n"); break; } f_t alpha_k = std::max(0.0, cover_size - 1.0 - objective); From ecce002092a5f2709cbffa8d3a48a31470cf2755 Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Fri, 20 Mar 2026 19:58:10 -0700 Subject: [PATCH 11/14] Use exact knapsack with integer values and fractional weights in lifting --- cpp/src/cuts/cuts.cpp | 80 ++++++++++++++++++++++++++++++++++++++++--- cpp/src/cuts/cuts.hpp | 5 +++ 2 files changed, 81 insertions(+), 4 deletions(-) diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index 162b2aba4..6115739a8 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -1221,7 +1221,7 @@ void knapsack_generation_t::lift_knapsack_cut( // Construct weight and values for C - std::vector values; + std::vector values; values.reserve(knapsack_inequality.size()); std::vector weights; @@ -1232,7 +1232,7 @@ void knapsack_generation_t::lift_knapsack_cut( if (!is_marked_[j]) { // j is in C weights.push_back(knapsack_inequality.coeff(k)); - values.push_back(1.0); + values.push_back(1); } } @@ -1257,7 +1257,7 @@ void knapsack_generation_t::lift_knapsack_cut( f_t capacity = knapsack_inequality.rhs - a_k; - f_t objective = solve_knapsack_problem(values, weights, capacity, solution); + f_t objective = exact_knapsack_problem_integer_values_fraction_values(values, weights, capacity, solution); if (std::isnan(objective)) { settings_.log.printf("lifting knapsack problem failed\n"); break; } f_t alpha_k = std::max(0.0, cover_size - 1.0 - objective); @@ -1266,7 +1266,7 @@ void knapsack_generation_t::lift_knapsack_cut( settings_.log.printf("Lifted variable %d with alpha %g\n", k, alpha_k); F.push_back(k); alpha.push_back(alpha_k); - values.push_back(alpha_k); + values.push_back(static_cast(alpha_k)); weights.push_back(a_k); lifted_cut.vector.i.push_back(k); @@ -1452,6 +1452,78 @@ f_t knapsack_generation_t::solve_knapsack_problem(const std::vector +f_t knapsack_generation_t::exact_knapsack_problem_integer_values_fraction_values( + const std::vector& values, + const std::vector& weights, + f_t rhs, + std::vector& solution) +{ + // Solve the knapsack problem + // maximize sum_{j=0}^n values[j] * solution[j] + // subject to sum_{j=0}^n weights[j] * solution[j] <= rhs + // values: values of the items + // weights: weights of the items + // return the value of the solution + + const i_t n = weights.size(); + + const bool verbose = false; + i_t sum_value = std::accumulate(values.begin(), values.end(), 0); + if (verbose) { settings_.log.printf("sum value %d\n", sum_value); } + const i_t max_size = 10000; + if (sum_value <= 0.0 || sum_value >= max_size) { + if (verbose) { + settings_.log.printf("sum value %d is negative or too large\n", + sum_value); + } + return std::numeric_limits::quiet_NaN(); + } + + solution.assign(n, 0.0); + + // dp(j, v) = minimum weight using first j items to get value v + dense_matrix_t dp(n + 1, sum_value + 1, inf); + dense_matrix_t take(n + 1, sum_value + 1, 0); + dp(0, 0) = 0; + + // 4. Dynamic programming + for (i_t j = 1; j <= n; ++j) { + for (i_t v = 0; v <= sum_value; ++v) { + // Do not take item i-1 + dp(j, v) = dp(j - 1, v); + + // Take item j-1 if possible + if (v >= values[j - 1]) { + f_t candidate = dp(j - 1, v - values[j - 1]) + weights[j - 1]; + if (candidate < dp(j, v)) { + dp(j, v) = candidate; + take(j, v) = 1; + } + } + } + } + + // 5. Find best achievable value within capacity + i_t best_value = 0; + for (i_t v = 0; v <= sum_value; ++v) { + if (dp(n, v) <= rhs) { best_value = v; } + } + + // 6. Backtrack to recover solution + i_t v = best_value; + for (i_t j = n; j >= 1; --j) { + if (take(j, v)) { + solution[j - 1] = 1.0; + v -= values[j - 1]; + } else { + solution[j - 1] = 0.0; + } + } + + return f_t(best_value); +} + template void cut_generation_t::generate_implied_bounds_cuts( const lp_problem_t& lp, diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index 2653ee793..c05ccd41b 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -376,6 +376,11 @@ class knapsack_generation_t { f_t rhs, std::vector& solution); + f_t exact_knapsack_problem_integer_values_fraction_values(const std::vector& values, + const std::vector& weights, + f_t rhs, + std::vector& solution); + std::vector is_slack_; std::vector knapsack_constraints_; std::vector is_complemented_; From dba801a1c9cf78ebd0e5e20810b6e621c102dc9c Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Fri, 20 Mar 2026 20:07:59 -0700 Subject: [PATCH 12/14] Add parameter for implied bound cuts --- cpp/include/cuopt/linear_programming/constants.h | 1 + cpp/include/cuopt/linear_programming/mip/solver_settings.hpp | 1 + cpp/src/math_optimization/solver_settings.cu | 1 + cpp/src/mip_heuristics/solver.cu | 1 + 4 files changed, 4 insertions(+) diff --git a/cpp/include/cuopt/linear_programming/constants.h b/cpp/include/cuopt/linear_programming/constants.h index 86becfe06..930b793ea 100644 --- a/cpp/include/cuopt/linear_programming/constants.h +++ b/cpp/include/cuopt/linear_programming/constants.h @@ -64,6 +64,7 @@ #define CUOPT_MIP_MIXED_INTEGER_ROUNDING_CUTS "mip_mixed_integer_rounding_cuts" #define CUOPT_MIP_MIXED_INTEGER_GOMORY_CUTS "mip_mixed_integer_gomory_cuts" #define CUOPT_MIP_KNAPSACK_CUTS "mip_knapsack_cuts" +#define CUOPT_MIP_IMPLIED_BOUNDS_CUTS "mip_implied_bounds_cuts" #define CUOPT_MIP_CLIQUE_CUTS "mip_clique_cuts" #define CUOPT_MIP_STRONG_CHVATAL_GOMORY_CUTS "mip_strong_chvatal_gomory_cuts" #define CUOPT_MIP_REDUCED_COST_STRENGTHENING "mip_reduced_cost_strengthening" diff --git a/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp b/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp index 2c92d2623..0a3dd6d25 100644 --- a/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp +++ b/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp @@ -94,6 +94,7 @@ class mip_solver_settings_t { i_t mixed_integer_gomory_cuts = -1; i_t knapsack_cuts = -1; i_t clique_cuts = -1; + i_t implied_bounds_cuts = -1; i_t strong_chvatal_gomory_cuts = -1; i_t reduced_cost_strengthening = -1; f_t cut_change_threshold = -1.0; diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index a60d508fa..7f23449fd 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -95,6 +95,7 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_MIP_MIXED_INTEGER_GOMORY_CUTS, &mip_settings.mixed_integer_gomory_cuts, -1, 1, -1}, {CUOPT_MIP_KNAPSACK_CUTS, &mip_settings.knapsack_cuts, -1, 1, -1}, {CUOPT_MIP_CLIQUE_CUTS, &mip_settings.clique_cuts, -1, 1, -1}, + {CUOPT_MIP_IMPLIED_BOUNDS_CUTS, &mip_settings.implied_bounds_cuts, -1, 1, -1}, {CUOPT_MIP_STRONG_CHVATAL_GOMORY_CUTS, &mip_settings.strong_chvatal_gomory_cuts, -1, 1, -1}, {CUOPT_MIP_REDUCED_COST_STRENGTHENING, &mip_settings.reduced_cost_strengthening, -1, std::numeric_limits::max(), -1}, {CUOPT_NUM_GPUS, &pdlp_settings.num_gpus, 1, 2, 1}, diff --git a/cpp/src/mip_heuristics/solver.cu b/cpp/src/mip_heuristics/solver.cu index 0e03a529f..6f479a45c 100644 --- a/cpp/src/mip_heuristics/solver.cu +++ b/cpp/src/mip_heuristics/solver.cu @@ -293,6 +293,7 @@ solution_t mip_solver_t::run_solver() branch_and_bound_settings.mixed_integer_gomory_cuts = context.settings.mixed_integer_gomory_cuts; branch_and_bound_settings.knapsack_cuts = context.settings.knapsack_cuts; + branch_and_bound_settings.implied_bounds_cuts = context.settings.implied_bounds_cuts; branch_and_bound_settings.clique_cuts = context.settings.clique_cuts; branch_and_bound_settings.strong_chvatal_gomory_cuts = context.settings.strong_chvatal_gomory_cuts; From f98a1f5e4712585722b11b68a2fac70e91025609 Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Fri, 20 Mar 2026 20:44:46 -0700 Subject: [PATCH 13/14] implied bounds -> implied bound --- .../cuopt/linear_programming/constants.h | 2 +- .../mip/solver_settings.hpp | 2 +- cpp/src/branch_and_bound/branch_and_bound.cpp | 6 +-- cpp/src/branch_and_bound/branch_and_bound.hpp | 4 +- cpp/src/cuts/cuts.cpp | 41 ++++++++-------- cpp/src/cuts/cuts.hpp | 25 +++++----- .../dual_simplex/simplex_solver_settings.hpp | 4 +- cpp/src/dual_simplex/solve.cpp | 4 +- cpp/src/math_optimization/solver_settings.cu | 2 +- cpp/src/mip_heuristics/diversity/lns/rins.cu | 2 +- .../diversity/recombiners/sub_mip.cuh | 2 +- cpp/src/mip_heuristics/solver.cu | 48 +++++++++---------- 12 files changed, 70 insertions(+), 72 deletions(-) diff --git a/cpp/include/cuopt/linear_programming/constants.h b/cpp/include/cuopt/linear_programming/constants.h index 930b793ea..ec0db2e7e 100644 --- a/cpp/include/cuopt/linear_programming/constants.h +++ b/cpp/include/cuopt/linear_programming/constants.h @@ -64,7 +64,7 @@ #define CUOPT_MIP_MIXED_INTEGER_ROUNDING_CUTS "mip_mixed_integer_rounding_cuts" #define CUOPT_MIP_MIXED_INTEGER_GOMORY_CUTS "mip_mixed_integer_gomory_cuts" #define CUOPT_MIP_KNAPSACK_CUTS "mip_knapsack_cuts" -#define CUOPT_MIP_IMPLIED_BOUNDS_CUTS "mip_implied_bounds_cuts" +#define CUOPT_MIP_IMPLIED_BOUND_CUTS "mip_implied_bound_cuts" #define CUOPT_MIP_CLIQUE_CUTS "mip_clique_cuts" #define CUOPT_MIP_STRONG_CHVATAL_GOMORY_CUTS "mip_strong_chvatal_gomory_cuts" #define CUOPT_MIP_REDUCED_COST_STRENGTHENING "mip_reduced_cost_strengthening" diff --git a/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp b/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp index 0a3dd6d25..8a1dc5232 100644 --- a/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp +++ b/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp @@ -94,7 +94,7 @@ class mip_solver_settings_t { i_t mixed_integer_gomory_cuts = -1; i_t knapsack_cuts = -1; i_t clique_cuts = -1; - i_t implied_bounds_cuts = -1; + i_t implied_bound_cuts = -1; i_t strong_chvatal_gomory_cuts = -1; i_t reduced_cost_strengthening = -1; f_t cut_change_threshold = -1.0; diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 1345cb3e7..418d493b4 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -243,11 +243,11 @@ branch_and_bound_t::branch_and_bound_t( const user_problem_t& user_problem, const simplex_solver_settings_t& solver_settings, f_t start_time, - const probing_implied_bounds_t& probing_implied_bounds, + const probing_implied_bound_t& probing_implied_bound, std::shared_ptr> clique_table) : original_problem_(user_problem), settings_(solver_settings), - probing_implied_bounds_(probing_implied_bounds), + probing_implied_bound_(probing_implied_bound), clique_table_(std::move(clique_table)), original_lp_(user_problem.handle_ptr, 1, 1, 1), Arow_(1, 1, 0), @@ -2146,7 +2146,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut new_slacks_, var_types_, original_problem_, - probing_implied_bounds_, + probing_implied_bound_, clique_table_, &clique_table_future_, &signal_extend_cliques_); diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index 3c5ee6538..bbb8e04bc 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -77,7 +77,7 @@ class branch_and_bound_t { branch_and_bound_t(const user_problem_t& user_problem, const simplex_solver_settings_t& solver_settings, f_t start_time, - const probing_implied_bounds_t& probing_implied_bounds, + const probing_implied_bound_t& probing_implied_bound, std::shared_ptr> clique_table = nullptr); // Set an initial guess based on the user_problem. This should be called before solve. @@ -149,7 +149,7 @@ class branch_and_bound_t { private: const user_problem_t& original_problem_; const simplex_solver_settings_t settings_; - const probing_implied_bounds_t& probing_implied_bounds_; + const probing_implied_bound_t& probing_implied_bound_; std::shared_ptr> clique_table_; std::future>> clique_table_future_; std::atomic signal_extend_cliques_{false}; diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index 6115739a8..098d9b8d1 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -1525,19 +1525,18 @@ f_t knapsack_generation_t::exact_knapsack_problem_integer_values_fract } template -void cut_generation_t::generate_implied_bounds_cuts( +void cut_generation_t::generate_implied_bound_cuts( const lp_problem_t& lp, const simplex_solver_settings_t& settings, const std::vector& var_types, const std::vector& xstar, - variable_bounds_t& variable_bounds, f_t start_time) { - if (probing_implied_bounds_.zero_offsets.empty()) { return; } + if (probing_implied_bound_.zero_offsets.empty()) { return; } const f_t tol = 1e-4; i_t num_cuts = 0; - const i_t pib_cols = static_cast(probing_implied_bounds_.zero_offsets.size()) - 1; + const i_t pib_cols = static_cast(probing_implied_bound_.zero_offsets.size()) - 1; const i_t n_cols = std::min(lp.num_cols, pib_cols); for (i_t j = 0; j < n_cols; j++) { @@ -1545,16 +1544,16 @@ void cut_generation_t::generate_implied_bounds_cuts( const f_t xstar_j = xstar[j]; // x_j = 0 implications - const i_t zero_begin = probing_implied_bounds_.zero_offsets[j]; - const i_t zero_end = probing_implied_bounds_.zero_offsets[j + 1]; + const i_t zero_begin = probing_implied_bound_.zero_offsets[j]; + const i_t zero_end = probing_implied_bound_.zero_offsets[j + 1]; for (i_t p = zero_begin; p < zero_end; p++) { - const i_t i = probing_implied_bounds_.zero_variables[p]; + const i_t i = probing_implied_bound_.zero_variables[p]; const f_t l_i = lp.lower[i]; const f_t u_i = lp.upper[i]; // Tightened upper bound: x_j = 0 implies y_i <= b, where b < u_i // Valid inequality: y_i <= b + (u_i - b)*x_j or -y_i + (u_i - b)*x_j >= -b - const f_t b_ub = probing_implied_bounds_.zero_upper_bound[p]; + const f_t b_ub = probing_implied_bound_.zero_upper_bound[p]; if (b_ub < u_i - tol) { const f_t coeff_j = u_i - b_ub; const f_t y_i = xstar[i]; @@ -1565,14 +1564,14 @@ void cut_generation_t::generate_implied_bounds_cuts( cut.push_back(i, -1.0); cut.push_back(j, coeff_j); cut.rhs = -b_ub; - cut_pool_.add_cut(cut_type_t::IMPLIED_BOUNDS, cut); + cut_pool_.add_cut(cut_type_t::IMPLIED_BOUND, cut); num_cuts++; } } // Tightened lower bound: x_j = 0 implies y_i >= b, where b > l_i // Valid inequality: y_i >= b - (b - l_i)*x_j or y_i + (b - l_i)*x_j >= b - const f_t b_lb = probing_implied_bounds_.zero_lower_bound[p]; + const f_t b_lb = probing_implied_bound_.zero_lower_bound[p]; if (b_lb > l_i + tol) { const f_t coeff_j = b_lb - l_i; const f_t y_i = xstar[i]; @@ -1583,23 +1582,23 @@ void cut_generation_t::generate_implied_bounds_cuts( cut.push_back(i, 1.0); cut.push_back(j, coeff_j); cut.rhs = b_lb; - cut_pool_.add_cut(cut_type_t::IMPLIED_BOUNDS, cut); + cut_pool_.add_cut(cut_type_t::IMPLIED_BOUND, cut); num_cuts++; } } } // x_j = 1 implications - const i_t one_begin = probing_implied_bounds_.one_offsets[j]; - const i_t one_end = probing_implied_bounds_.one_offsets[j + 1]; + const i_t one_begin = probing_implied_bound_.one_offsets[j]; + const i_t one_end = probing_implied_bound_.one_offsets[j + 1]; for (i_t p = one_begin; p < one_end; p++) { - const i_t i = probing_implied_bounds_.one_variables[p]; + const i_t i = probing_implied_bound_.one_variables[p]; const f_t l_i = lp.lower[i]; const f_t u_i = lp.upper[i]; // Tightened upper bound: x_j = 1 implies y_i <= b, where b < u_i // Valid inequality: y_i <= u_i - (u_i - b)*x_j or -y_i - (u_i - b)*x_j >= -u_i - const f_t b_ub = probing_implied_bounds_.one_upper_bound[p]; + const f_t b_ub = probing_implied_bound_.one_upper_bound[p]; if (b_ub < u_i - tol) { const f_t coeff_j = -(u_i - b_ub); const f_t y_i = xstar[i]; @@ -1610,14 +1609,14 @@ void cut_generation_t::generate_implied_bounds_cuts( cut.push_back(i, -1.0); cut.push_back(j, coeff_j); cut.rhs = -u_i; - cut_pool_.add_cut(cut_type_t::IMPLIED_BOUNDS, cut); + cut_pool_.add_cut(cut_type_t::IMPLIED_BOUND, cut); num_cuts++; } } // Tightened lower bound: x_j = 1 implies y_i >= b, where b > l_i // Valid inequality: y_i >= l_i + (b - l_i)*x_j or y_i - (b - l_i)*x_j >= l_i - const f_t b_lb = probing_implied_bounds_.one_lower_bound[p]; + const f_t b_lb = probing_implied_bound_.one_lower_bound[p]; if (b_lb > l_i + tol) { const f_t coeff_j = -(b_lb - l_i); const f_t lhs_val = xstar[i] + coeff_j * xstar_j; @@ -1627,7 +1626,7 @@ void cut_generation_t::generate_implied_bounds_cuts( cut.push_back(i, 1.0); cut.push_back(j, coeff_j); cut.rhs = rhs_val; - cut_pool_.add_cut(cut_type_t::IMPLIED_BOUNDS, cut); + cut_pool_.add_cut(cut_type_t::IMPLIED_BOUND, cut); num_cuts++; } } @@ -1699,10 +1698,10 @@ bool cut_generation_t::generate_cuts(const lp_problem_t& lp, } } - // Generate implied bounds cuts - if (settings.implied_bounds_cuts != 0) { + // Generate implied bound cuts + if (settings.implied_bound_cuts != 0) { f_t cut_start_time = tic(); - generate_implied_bounds_cuts(lp, settings, var_types, xstar, variable_bounds, start_time); + generate_implied_bound_cuts(lp, settings, var_types, xstar, start_time); f_t cut_generation_time = toc(cut_start_time); if (cut_generation_time > 1.0) { settings.log.debug("Implied bounds cut generation time %.2f seconds\n", cut_generation_time); diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index c05ccd41b..61767ecc8 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -38,7 +38,7 @@ enum cut_type_t : int8_t { KNAPSACK = 2, CHVATAL_GOMORY = 3, CLIQUE = 4, - IMPLIED_BOUNDS = 5, + IMPLIED_BOUND = 5, MAX_CUT_TYPE = 6 }; @@ -64,7 +64,7 @@ cut_gap_closure_t compute_cut_gap_closure(f_t objective_reference, } template -struct probing_implied_bounds_t { +struct probing_implied_bound_t { // Probing implications stored in CSR format, indexed by binary variable x_j. // // "zero" = implications discovered when probing x_j = 0. @@ -85,9 +85,9 @@ struct probing_implied_bounds_t { // Non-binary variables have empty ranges (zero_offsets[j] == zero_offsets[j+1]). // Offsets vectors have size num_cols + 1. - probing_implied_bounds_t() = default; + probing_implied_bound_t() = default; - probing_implied_bounds_t(i_t num_cols) + probing_implied_bound_t(i_t num_cols) : zero_offsets(num_cols + 1, 0), one_offsets(num_cols + 1, 0) { } @@ -408,14 +408,14 @@ class cut_generation_t { const std::vector& new_slacks, const std::vector& var_types, const user_problem_t& user_problem, - const probing_implied_bounds_t& probing_implied_bounds, + const probing_implied_bound_t& probing_implied_bound, std::shared_ptr> clique_table = nullptr, std::future>>* clique_table_future = nullptr, std::atomic* signal_extend = nullptr) : cut_pool_(cut_pool), knapsack_generation_(lp, settings, Arow, new_slacks, var_types), user_problem_(user_problem), - probing_implied_bounds_(probing_implied_bounds), + probing_implied_bound_(probing_implied_bound), clique_table_(std::move(clique_table)), clique_table_future_(clique_table_future), signal_extend_(signal_extend) @@ -475,17 +475,16 @@ class cut_generation_t { f_t start_time); // Generate implied bounds cuts from probing implications - void generate_implied_bounds_cuts(const lp_problem_t& lp, - const simplex_solver_settings_t& settings, - const std::vector& var_types, - const std::vector& xstar, - variable_bounds_t& variable_bounds, - f_t start_time); + void generate_implied_bound_cuts(const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + const std::vector& var_types, + const std::vector& xstar, + f_t start_time); cut_pool_t& cut_pool_; knapsack_generation_t knapsack_generation_; const user_problem_t& user_problem_; - const probing_implied_bounds_t& probing_implied_bounds_; + const probing_implied_bound_t& probing_implied_bound_; std::shared_ptr> clique_table_; std::future>>* clique_table_future_{nullptr}; std::atomic* signal_extend_{nullptr}; diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index 67e52f9f7..170bdb1be 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -100,7 +100,7 @@ struct simplex_solver_settings_t { mir_cuts(-1), mixed_integer_gomory_cuts(-1), knapsack_cuts(-1), - implied_bounds_cuts(-1), + implied_bound_cuts(-1), clique_cuts(-1), strong_chvatal_gomory_cuts(-1), reduced_cost_strengthening(-1), @@ -181,7 +181,7 @@ struct simplex_solver_settings_t { i_t mixed_integer_gomory_cuts; // -1 automatic, 0 to disable, >0 to enable mixed integer Gomory // cuts i_t knapsack_cuts; // -1 automatic, 0 to disable, >0 to enable knapsack cuts - i_t implied_bounds_cuts; // -1 automatic, 0 to disable, >0 to enable implied bounds cuts + i_t implied_bound_cuts; // -1 automatic, 0 to disable, >0 to enable implied bound cuts i_t clique_cuts; // -1 automatic, 0 to disable, >0 to enable clique cuts i_t strong_chvatal_gomory_cuts; // -1 automatic, 0 to disable, >0 to enable strong Chvatal Gomory // cuts diff --git a/cpp/src/dual_simplex/solve.cpp b/cpp/src/dual_simplex/solve.cpp index c0f919020..5fb2f6c6e 100644 --- a/cpp/src/dual_simplex/solve.cpp +++ b/cpp/src/dual_simplex/solve.cpp @@ -706,7 +706,7 @@ i_t solve(const user_problem_t& problem, { i_t status; if (is_mip(problem) && !settings.relaxation) { - probing_implied_bounds_t empty_probing(problem.num_cols); + probing_implied_bound_t empty_probing(problem.num_cols); branch_and_bound_t branch_and_bound(problem, settings, tic(), empty_probing); mip_solution_t mip_solution(problem.num_cols); mip_status_t mip_status = branch_and_bound.solve(mip_solution); @@ -746,7 +746,7 @@ i_t solve_mip_with_guess(const user_problem_t& problem, { i_t status; if (is_mip(problem)) { - probing_implied_bounds_t empty_probing(problem.num_cols); + probing_implied_bound_t empty_probing(problem.num_cols); branch_and_bound_t branch_and_bound(problem, settings, tic(), empty_probing); branch_and_bound.set_initial_guess(guess); mip_status_t mip_status = branch_and_bound.solve(solution); diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index 7f23449fd..0ac5633d2 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -95,7 +95,7 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_MIP_MIXED_INTEGER_GOMORY_CUTS, &mip_settings.mixed_integer_gomory_cuts, -1, 1, -1}, {CUOPT_MIP_KNAPSACK_CUTS, &mip_settings.knapsack_cuts, -1, 1, -1}, {CUOPT_MIP_CLIQUE_CUTS, &mip_settings.clique_cuts, -1, 1, -1}, - {CUOPT_MIP_IMPLIED_BOUNDS_CUTS, &mip_settings.implied_bounds_cuts, -1, 1, -1}, + {CUOPT_MIP_IMPLIED_BOUND_CUTS, &mip_settings.implied_bound_cuts, -1, 1, -1}, {CUOPT_MIP_STRONG_CHVATAL_GOMORY_CUTS, &mip_settings.strong_chvatal_gomory_cuts, -1, 1, -1}, {CUOPT_MIP_REDUCED_COST_STRENGTHENING, &mip_settings.reduced_cost_strengthening, -1, std::numeric_limits::max(), -1}, {CUOPT_NUM_GPUS, &pdlp_settings.num_gpus, 1, 2, 1}, diff --git a/cpp/src/mip_heuristics/diversity/lns/rins.cu b/cpp/src/mip_heuristics/diversity/lns/rins.cu index de33d5487..25b9d14cb 100644 --- a/cpp/src/mip_heuristics/diversity/lns/rins.cu +++ b/cpp/src/mip_heuristics/diversity/lns/rins.cu @@ -275,7 +275,7 @@ void rins_t::run_rins() f_t objective) { rins_solution_queue.push_back(solution); }; - dual_simplex::probing_implied_bounds_t empty_probing( + dual_simplex::probing_implied_bound_t empty_probing( branch_and_bound_problem.num_cols); dual_simplex::branch_and_bound_t branch_and_bound( branch_and_bound_problem, branch_and_bound_settings, dual_simplex::tic(), empty_probing); diff --git a/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh b/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh index 42e33d214..f6e4f172c 100644 --- a/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh +++ b/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh @@ -117,7 +117,7 @@ class sub_mip_recombiner_t : public recombiner_t { // disable B&B logs, so that it is not interfering with the main B&B thread branch_and_bound_settings.log.log = false; - dual_simplex::probing_implied_bounds_t empty_probing( + dual_simplex::probing_implied_bound_t empty_probing( branch_and_bound_problem.num_cols); dual_simplex::branch_and_bound_t branch_and_bound( branch_and_bound_problem, branch_and_bound_settings, dual_simplex::tic(), empty_probing); diff --git a/cpp/src/mip_heuristics/solver.cu b/cpp/src/mip_heuristics/solver.cu index 6f479a45c..f477613f7 100644 --- a/cpp/src/mip_heuristics/solver.cu +++ b/cpp/src/mip_heuristics/solver.cu @@ -192,7 +192,7 @@ solution_t mip_solver_t::run_solver() branch_and_bound_solution_helper_t solution_helper(&dm, branch_and_bound_settings); dual_simplex::mip_solution_t branch_and_bound_solution(1); - dual_simplex::probing_implied_bounds_t probing_implied_bounds; + dual_simplex::probing_implied_bound_t probing_implied_bound; bool run_bb = !context.settings.heuristics_only; if (run_bb) { @@ -204,7 +204,7 @@ solution_t mip_solver_t::run_solver() // Extract probing cache into CPU-only CSR struct for implied bounds cuts { const i_t num_cols = branch_and_bound_problem.num_cols; - probing_implied_bounds = dual_simplex::probing_implied_bounds_t(num_cols); + probing_implied_bound = dual_simplex::probing_implied_bound_t(num_cols); auto& pc = dm.ls.constraint_prop.bounds_update.probing_cache.probing_cache; auto& rev_ids = context.problem_ptr->reverse_original_ids; @@ -217,34 +217,34 @@ solution_t mip_solver_t::run_solver() for (auto& [imp_var, bound] : entries[0].var_to_cached_bound_map) { i_t i = (imp_var < static_cast(rev_ids.size())) ? rev_ids[imp_var] : -1; if (i < 0 || i >= num_cols) { continue; } - probing_implied_bounds.zero_offsets[j + 1]++; + probing_implied_bound.zero_offsets[j + 1]++; } for (auto& [imp_var, bound] : entries[1].var_to_cached_bound_map) { i_t i = (imp_var < static_cast(rev_ids.size())) ? rev_ids[imp_var] : -1; if (i < 0 || i >= num_cols) { continue; } - probing_implied_bounds.one_offsets[j + 1]++; + probing_implied_bound.one_offsets[j + 1]++; } } // Prefix sum for (i_t j = 0; j < num_cols; j++) { - probing_implied_bounds.zero_offsets[j + 1] += probing_implied_bounds.zero_offsets[j]; - probing_implied_bounds.one_offsets[j + 1] += probing_implied_bounds.one_offsets[j]; + probing_implied_bound.zero_offsets[j + 1] += probing_implied_bound.zero_offsets[j]; + probing_implied_bound.one_offsets[j + 1] += probing_implied_bound.one_offsets[j]; } // Allocate flat arrays - i_t zero_nnz = probing_implied_bounds.zero_offsets[num_cols]; - i_t one_nnz = probing_implied_bounds.one_offsets[num_cols]; - probing_implied_bounds.zero_variables.resize(zero_nnz); - probing_implied_bounds.zero_lower_bound.resize(zero_nnz); - probing_implied_bounds.zero_upper_bound.resize(zero_nnz); - probing_implied_bounds.one_variables.resize(one_nnz); - probing_implied_bounds.one_lower_bound.resize(one_nnz); - probing_implied_bounds.one_upper_bound.resize(one_nnz); + i_t zero_nnz = probing_implied_bound.zero_offsets[num_cols]; + i_t one_nnz = probing_implied_bound.one_offsets[num_cols]; + probing_implied_bound.zero_variables.resize(zero_nnz); + probing_implied_bound.zero_lower_bound.resize(zero_nnz); + probing_implied_bound.zero_upper_bound.resize(zero_nnz); + probing_implied_bound.one_variables.resize(one_nnz); + probing_implied_bound.one_lower_bound.resize(one_nnz); + probing_implied_bound.one_upper_bound.resize(one_nnz); // Second pass: fill flat arrays using write cursors - std::vector zero_cursor(probing_implied_bounds.zero_offsets); - std::vector one_cursor(probing_implied_bounds.one_offsets); + std::vector zero_cursor(probing_implied_bound.zero_offsets); + std::vector one_cursor(probing_implied_bound.one_offsets); for (auto& [var_idx, entries] : pc) { if (entries[0].val_interval.interval_type != interval_type_t::EQUALS) { continue; } @@ -255,17 +255,17 @@ solution_t mip_solver_t::run_solver() i_t i = (imp_var < static_cast(rev_ids.size())) ? rev_ids[imp_var] : -1; if (i < 0 || i >= num_cols) { continue; } i_t p = zero_cursor[j]++; - probing_implied_bounds.zero_variables[p] = i; - probing_implied_bounds.zero_lower_bound[p] = bound.lb; - probing_implied_bounds.zero_upper_bound[p] = bound.ub; + probing_implied_bound.zero_variables[p] = i; + probing_implied_bound.zero_lower_bound[p] = bound.lb; + probing_implied_bound.zero_upper_bound[p] = bound.ub; } for (auto& [imp_var, bound] : entries[1].var_to_cached_bound_map) { i_t i = (imp_var < static_cast(rev_ids.size())) ? rev_ids[imp_var] : -1; if (i < 0 || i >= num_cols) { continue; } i_t p = one_cursor[j]++; - probing_implied_bounds.one_variables[p] = i; - probing_implied_bounds.one_lower_bound[p] = bound.lb; - probing_implied_bounds.one_upper_bound[p] = bound.ub; + probing_implied_bound.one_variables[p] = i; + probing_implied_bound.one_lower_bound[p] = bound.lb; + probing_implied_bound.one_upper_bound[p] = bound.ub; } } @@ -293,7 +293,7 @@ solution_t mip_solver_t::run_solver() branch_and_bound_settings.mixed_integer_gomory_cuts = context.settings.mixed_integer_gomory_cuts; branch_and_bound_settings.knapsack_cuts = context.settings.knapsack_cuts; - branch_and_bound_settings.implied_bounds_cuts = context.settings.implied_bounds_cuts; + branch_and_bound_settings.implied_bound_cuts = context.settings.implied_bound_cuts; branch_and_bound_settings.clique_cuts = context.settings.clique_cuts; branch_and_bound_settings.strong_chvatal_gomory_cuts = context.settings.strong_chvatal_gomory_cuts; @@ -343,7 +343,7 @@ solution_t mip_solver_t::run_solver() branch_and_bound_problem, branch_and_bound_settings, timer_.get_tic_start(), - probing_implied_bounds, + probing_implied_bound, context.problem_ptr->clique_table); context.branch_and_bound_ptr = branch_and_bound.get(); auto* stats_ptr = &context.stats; From 23d82954daa38bd1d5a41aeb0180abd4b407b24f Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Fri, 20 Mar 2026 20:52:05 -0700 Subject: [PATCH 14/14] Don't generate MIR cut off the same base row --- cpp/src/cuts/cuts.cpp | 2 +- cpp/src/cuts/cuts.hpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index 098d9b8d1..b9b0aebd1 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -2276,7 +2276,7 @@ void cut_generation_t::generate_mir_cuts( starting_count[i]++; scores[i] = initial_scores[i] * std::pow(0.9, starting_count[i]); - score_queue.push(std::make_pair(scores[i], i)); + //score_queue.push(std::make_pair(scores[i], i)); work_estimate += std::log2(std::max(1, static_cast(score_queue.size()))); } } diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index 61767ecc8..e71c2557c 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -173,7 +173,7 @@ struct cut_info_t { } } const char* cut_type_names[MAX_CUT_TYPE] = { - "Gomory ", "MIR ", "Knapsack ", "Strong CG", "Clique ", "Implied Bounds"}; + "Gomory ", "MIR ", "Knapsack ", "Strong CG ", "Clique ", "Implied Bounds"}; std::array num_cuts = {0}; };