From 26a268cbb3b6d509d4b8d896b213eb09b8fa59ff Mon Sep 17 00:00:00 2001 From: Hugo Linsenmaier Date: Mon, 6 Apr 2026 16:05:59 -0700 Subject: [PATCH 01/17] Implement flow cover cuts --- .../cuopt/linear_programming/constants.h | 1 + .../mip/solver_settings.hpp | 1 + cpp/src/cuts/cuts.cpp | 332 ++++++++++++++++++ cpp/src/cuts/cuts.hpp | 28 +- .../dual_simplex/simplex_solver_settings.hpp | 2 + cpp/src/math_optimization/solver_settings.cu | 1 + cpp/src/mip_heuristics/solver.cu | 1 + cpp/tests/mip/cuts_test.cu | 134 +++++++ 8 files changed, 498 insertions(+), 2 deletions(-) diff --git a/cpp/include/cuopt/linear_programming/constants.h b/cpp/include/cuopt/linear_programming/constants.h index 06eacb3408..69b1f223e9 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_FLOW_COVER_CUTS "mip_flow_cover_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" diff --git a/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp b/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp index 14c4d227bc..2b4ee28eb2 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 mir_cuts = -1; i_t mixed_integer_gomory_cuts = -1; i_t knapsack_cuts = -1; + i_t flow_cover_cuts = -1; i_t clique_cuts = -1; i_t implied_bound_cuts = -1; i_t strong_chvatal_gomory_cuts = -1; diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index 6d7d97ef0a..5986bbca6b 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -898,6 +898,307 @@ knapsack_generation_t::knapsack_generation_t( i_t num_knapsack_constraints = knapsack_constraints_.size(); settings.log.printf("Number of knapsack constraints %d\n", num_knapsack_constraints); #endif + + // Identify flow cover constraints: rows with binary variables having both + // positive and negative coefficients (mixed 0-1 structure). + // The constraint has the form: sum_{j in N1} a_j x_j <= b + sum_{j in N2} a_j x_j + // which is stored as: sum_{j in N1} a_j x_j - sum_{j in N2} a_j x_j <= b + flow_cover_constraints_.reserve(lp.num_rows); + for (i_t i = 0; i < lp.num_rows; i++) { + 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 valid = true; + bool has_pos = false; + bool has_neg = false; + i_t num_binaries = 0; + for (i_t p = 0; p < row_len; p++) { + const i_t j = inequality.index(p); + if (is_slack_[j]) { + if (inequality.coeff(p) < 0.0) { + valid = false; + break; + } + continue; + } + if (var_types[j] != variable_type_t::INTEGER || lp.lower[j] != 0.0 || lp.upper[j] != 1.0) { + valid = false; + break; + } + const f_t aj = inequality.coeff(p); + if (aj > 0.0) { + has_pos = true; + } else if (aj < 0.0) { + has_neg = true; + } + num_binaries++; + } + + if (valid && has_pos && has_neg && num_binaries >= 3) { flow_cover_constraints_.push_back(i); } + } +} + +template +i_t knapsack_generation_t::generate_flow_cover_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 flow_cover_row, + inequality_t& cut) +{ + const bool verbose = false; + // Get the row: sum_{j in N1} a_j x_j - sum_{j in N2} a_j x_j <= b + // which represents: sum_{j in N1} y_j <= b + sum_{j in N2} y_j + // where y_j = a_j * x_j for binary x_j. + inequality_t inequality(Arow, flow_cover_row, lp.rhs[flow_cover_row]); + inequality_t rational_inequality = inequality; + if (!rational_coefficients(var_types, inequality, rational_inequality)) { return -1; } + inequality = rational_inequality; + + // Separate into N1 (positive coefficients) and N2 (negative coefficients) + // N1: j with a_j > 0 (left-hand side flows) + // N2: j with a_j < 0 (right-hand side flows, stored as negative in the constraint) + struct flow_var_t { + i_t j; // variable index + f_t a_j; // absolute value of coefficient + }; + std::vector N1, N2; + N1.reserve(inequality.size()); + N2.reserve(inequality.size()); + + f_t b = inequality.rhs; // RHS of the constraint + + for (i_t k = 0; k < inequality.size(); k++) { + const i_t j = inequality.index(k); + const f_t aj = inequality.coeff(k); + if (is_slack_[j]) { continue; } + if (aj > 0.0) { + N1.push_back({j, aj}); + } else if (aj < 0.0) { + N2.push_back({j, -aj}); // store absolute value + } + } + + if (N1.empty() || N2.empty()) { return -1; } + + // Solve the separation knapsack problem to find cover C = (C1, C2) + // minimize sum_{j in N1} (1 - x_j*) z_j + sum_{j in N2} (-x_j*) z_j (note: N2 values negative) + // subject to sum_{j in N1} a_j z_j - sum_{j in N2} a_j z_j > b + // z in {0,1}^n + // + // Convert to maximization knapsack (zbar_j = 1 - z_j): + // maximize sum_{j in N1} (1 - x_j*) zbar_j + sum_{j in N2} x_j* zbar_j + // subject to sum_{j in N1} a_j zbar_j + sum_{j in N2} a_j zbar_j <= sum_all_a - b - epsilon + // + // Wait, let's be more careful. The separation problem objective is: + // zeta = sum_{j in N1} (1 - x_j*) z_j - sum_{j in N2} x_j* z_j + // We want to minimize this. A violated flow cover exists when zeta < 1 - sum_{j in N2} x_j*. + // + // Convert to maximize by flipping with zbar = 1 - z: + // maximize sum_{j in N1} (1 - x_j*) zbar_j + sum_{j in N2} x_j* zbar_j + // subject to sum_{j in N1} a_j (1 - zbar_j) - sum_{j in N2} a_j (1 - zbar_j) > b + // which is: sum_N1 a_j - sum_N1 a_j zbar_j - sum_N2 a_j + sum_N2 a_j zbar_j > b + // => -sum_N1 a_j zbar_j + sum_N2 a_j zbar_j > b - sum_N1 a_j + sum_N2 a_j + // => sum_N1 a_j zbar_j - sum_N2 a_j zbar_j < sum_N1 a_j - sum_N2 a_j - b + // + // This is a knapsack: maximize value s.t. sum weights <= capacity + // where weights for N1 items = a_j, weights for N2 items = -a_j (negative weights!) + // + // Negative weights complicate this. Instead, let's use the greedy approach directly + // on the separation problem. We solve: + // minimize sum_{j in N1} (1 - x_j*) z_j - sum_{j in N2} x_j* z_j + // s.t. sum_{j in N1} a_j z_j - sum_{j in N2} a_j z_j > b + // + // Greedy: sort N1 by (1-x_j*)/a_j ascending, N2 by x_j*/a_j descending + // Include all of N2 first (they help the constraint and reduce objective), + // then greedily add N1 items. + + // Greedy separation heuristic + f_t constraint_lhs = 0.0; + f_t objective = 0.0; + + // Start by including all N2 variables in the cover (z_j = 1 for j in N2) + // This contributes -a_j to constraint and -x_j* to objective + std::vector in_cover_n2(N2.size(), true); + for (i_t k = 0; k < static_cast(N2.size()); k++) { + constraint_lhs -= N2[k].a_j; + objective -= xstar[N2[k].j]; + } + + // Sort N1 by (1 - x_j*) / a_j ratio (ascending = cheapest first) + std::vector n1_perm(N1.size()); + std::iota(n1_perm.begin(), n1_perm.end(), 0); + std::sort(n1_perm.begin(), n1_perm.end(), [&](i_t a, i_t idx_b) { + f_t ratio_a = (1.0 - xstar[N1[a].j]) / N1[a].a_j; + f_t ratio_b = (1.0 - xstar[N1[idx_b].j]) / N1[idx_b].a_j; + return ratio_a < ratio_b; + }); + + // Add N1 variables greedily until constraint > b + std::vector in_cover_n1(N1.size(), false); + for (i_t idx : n1_perm) { + constraint_lhs += N1[idx].a_j; + objective += (1.0 - xstar[N1[idx].j]); + in_cover_n1[idx] = true; + if (constraint_lhs > b) { break; } + } + + if (constraint_lhs <= b) { return -1; } // Could not find a cover + + // Remove N2 variables that are expensive (high x_j*) and not needed + // Sort N2 by x_j*/a_j descending (most expensive first) + std::vector n2_perm(N2.size()); + std::iota(n2_perm.begin(), n2_perm.end(), 0); + std::sort(n2_perm.begin(), n2_perm.end(), [&](i_t a, i_t idx_b) { + f_t ratio_a = xstar[N2[a].j] / N2[a].a_j; + f_t ratio_b = xstar[N2[idx_b].j] / N2[idx_b].a_j; + return ratio_a > ratio_b; + }); + + for (i_t idx : n2_perm) { + // Try removing this N2 variable from the cover + if (constraint_lhs + N2[idx].a_j > b) { + constraint_lhs += N2[idx].a_j; + objective += xstar[N2[idx].j]; + in_cover_n2[idx] = false; + } + } + + // Build C1 (subset of N1 in cover) and C2 (subset of N2 in cover) + std::vector C1_indices, C2_indices; + std::vector C1_coeffs, C2_coeffs; + for (i_t k = 0; k < static_cast(N1.size()); k++) { + if (in_cover_n1[k]) { + C1_indices.push_back(N1[k].j); + C1_coeffs.push_back(N1[k].a_j); + } + } + for (i_t k = 0; k < static_cast(N2.size()); k++) { + if (in_cover_n2[k]) { + C2_indices.push_back(N2[k].j); + C2_coeffs.push_back(N2[k].a_j); + } + } + + if (C1_indices.empty()) { return -1; } + + // Compute lambda = sum_{j in C1} a_j - sum_{j in C2} a_j - b + f_t sum_C1 = std::accumulate(C1_coeffs.begin(), C1_coeffs.end(), 0.0); + f_t sum_C2 = std::accumulate(C2_coeffs.begin(), C2_coeffs.end(), 0.0); + f_t lambda = sum_C1 - sum_C2 - b; + + if (lambda <= 1e-10) { return -1; } + + // Determine L2 = {j in N2 \ C2 : lambda * x_j* < a_j * x_j*} + // = {j in N2 \ C2 : lambda < a_j} (for x_j* > 0) + // Actually from the theory: L2 = {j in N2 \ C2 : lambda * x_j* < y_j*} + // where y_j* = a_j * x_j* for binary variables. + // So: lambda * x_j* < a_j * x_j*, which simplifies to lambda < a_j when x_j* > 0. + // When x_j* = 0, both sides are 0, so j is not in L2. + std::vector L2_indices; + std::vector L2_coeffs; + std::vector rest_N2_indices; // N2 \ (C2 union L2) + std::vector rest_N2_coeffs; + + for (i_t k = 0; k < static_cast(N2.size()); k++) { + if (in_cover_n2[k]) { continue; } // skip C2 + const i_t j = N2[k].j; + const f_t aj = N2[k].a_j; + if (xstar[j] > 1e-10 && lambda < aj) { + L2_indices.push_back(j); + L2_coeffs.push_back(aj); + } else { + rest_N2_indices.push_back(j); + rest_N2_coeffs.push_back(aj); + } + } + + // Evaluate the flow cover inequality violation: + // LHS = sum_{j in C1} y_j* + sum_{j in C1} (a_j - lambda)^+ (1 - x_j*) + // RHS = b + sum_{j in C2} a_j + lambda * sum_{j in L2} x_j* + sum_{j in rest_N2} y_j* + f_t lhs = 0.0; + for (i_t k = 0; k < static_cast(C1_indices.size()); k++) { + const i_t j = C1_indices[k]; + const f_t aj = C1_coeffs[k]; + lhs += aj * xstar[j]; // y_j* = a_j * x_j* + lhs += std::max(0.0, aj - lambda) * (1.0 - xstar[j]); // (a_j - lambda)^+ (1 - x_j*) + } + + f_t rhs = b + sum_C2; + for (i_t k = 0; k < static_cast(L2_indices.size()); k++) { + rhs += lambda * xstar[L2_indices[k]]; + } + for (i_t k = 0; k < static_cast(rest_N2_indices.size()); k++) { + rhs += rest_N2_coeffs[k] * xstar[rest_N2_indices[k]]; // y_j* = a_j * x_j* + } + + f_t violation = lhs - rhs; + const f_t tol = 1e-6; + if (violation <= tol) { return -1; } + + if (verbose) { + settings.log.printf( + "Flow cover cut row %d violation %e lambda %e\n", flow_cover_row, violation, lambda); + } + + // Build the cut inequality: LHS <= RHS + // sum_{j in C1} a_j x_j + sum_{j in C1} (a_j - lambda)^+ (1 - x_j) <= b + sum_C2 + lambda * + // sum_{j in L2} x_j + sum_{j in rest_N2} a_j x_j + // + // Rearranging to >= form for cut pool (cut'*x >= rhs): + // -sum_{j in C1} a_j x_j - sum_{j in C1} (a_j - lambda)^+ (1 - x_j) + lambda sum_{j in L2} x_j + // + sum_{j in rest_N2} a_j x_j >= -(b + sum_C2) + // + // Expanding (a_j - lambda)^+ (1 - x_j): + // For j with a_j > lambda: (a_j - lambda)(1 - x_j) = (a_j - lambda) - (a_j - lambda) x_j + // So the x_j coefficient for C1 becomes: -a_j + (a_j - lambda)^+ = -min(a_j, lambda) + // And there's a constant: -sum_{j in C1} (a_j - lambda)^+ + // + // In >= form: + // sum_{j in C1} min(a_j, lambda) x_j + lambda sum_{j in L2} x_j + sum_{j in rest_N2} a_j x_j + // >= sum_{j in C1} (a_j - lambda)^+ - (b + sum_C2) + // = sum_{j in C1} a_j - |C1^+| lambda ... hmm let me just build it directly. + + // Build cut in the form: cut'*x >= cut.rhs + // From: sum_{j in C1} y_j + sum_{j in C1} (a_j-lambda)^+(1-x_j) <= b + sum_C2 + lambda*sum_L2 + // x_j + sum_rest a_j x_j Negate for >= form. + cut.clear(); + cut.reserve(C1_indices.size() + L2_indices.size() + rest_N2_indices.size()); + + f_t constant = 0.0; + for (i_t k = 0; k < static_cast(C1_indices.size()); k++) { + const i_t j = C1_indices[k]; + const f_t aj = C1_coeffs[k]; + const f_t lift = std::max(0.0, aj - lambda); + // Coefficient of x_j: -a_j + lift = -(a_j - lift) = -min(a_j, lambda) + // But also the constant from lift*(1-x_j): -lift + f_t coeff = -aj + lift; // = -min(a_j, lambda) + cut.push_back(j, coeff); + constant -= lift; + } + for (i_t k = 0; k < static_cast(L2_indices.size()); k++) { + cut.push_back(L2_indices[k], lambda); + } + for (i_t k = 0; k < static_cast(rest_N2_indices.size()); k++) { + cut.push_back(rest_N2_indices[k], rest_N2_coeffs[k]); + } + cut.rhs = -(b + sum_C2) + constant; + + // Verify violation + f_t dot = cut.vector.dot(xstar); + f_t cut_violation = dot - cut.rhs; + if (cut_violation >= -tol) { return -1; } + + cut.sort(); + return 0; } template @@ -1803,6 +2104,16 @@ bool cut_generation_t::generate_cuts(const lp_problem_t& lp, } } + // Generate Flow Cover cuts + if (settings.flow_cover_cuts != 0) { + f_t cut_start_time = tic(); + generate_flow_cover_cuts(lp, settings, Arow, new_slacks, var_types, xstar, start_time); + f_t cut_generation_time = toc(cut_start_time); + if (cut_generation_time > 1.0) { + settings.log.debug("Flow cover cut generation time %.2f seconds\n", cut_generation_time); + } + } + // Generate MIR and CG cuts if (settings.mir_cuts != 0 || settings.strong_chvatal_gomory_cuts != 0) { f_t cut_start_time = tic(); @@ -1860,6 +2171,27 @@ void cut_generation_t::generate_knapsack_cuts( } } +template +void cut_generation_t::generate_flow_cover_cuts( + 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, + f_t start_time) +{ + if (knapsack_generation_.num_flow_cover_constraints() > 0) { + for (i_t flow_cover_row : knapsack_generation_.get_flow_cover_constraints()) { + if (toc(start_time) >= settings.time_limit) { return; } + inequality_t cut(lp.num_cols); + i_t status = knapsack_generation_.generate_flow_cover_cut( + lp, settings, Arow, new_slacks, var_types, xstar, flow_cover_row, cut); + if (status == 0) { cut_pool_.add_cut(cut_type_t::FLOW_COVER, cut); } + } + } +} + template bool cut_generation_t::generate_clique_cuts( const lp_problem_t& lp, diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index 2da9760e27..464746edc4 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -39,7 +39,8 @@ enum cut_type_t : int8_t { CHVATAL_GOMORY = 3, CLIQUE = 4, IMPLIED_BOUND = 5, - MAX_CUT_TYPE = 6 + FLOW_COVER = 6, + MAX_CUT_TYPE = 7 }; template @@ -178,7 +179,8 @@ struct cut_info_t { "Knapsack ", "Strong CG ", "Clique ", - "Implied Bounds"}; + "Implied Bounds", + "Flow Cover "}; std::array num_cuts = {0}; }; @@ -345,9 +347,21 @@ class knapsack_generation_t { i_t knapsack_row, inequality_t& cut); + i_t generate_flow_cover_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 flow_cover_row, + inequality_t& cut); + i_t num_knapsack_constraints() const { return knapsack_constraints_.size(); } const std::vector& get_knapsack_constraints() const { return knapsack_constraints_; } + i_t num_flow_cover_constraints() const { return flow_cover_constraints_.size(); } + const std::vector& get_flow_cover_constraints() const { return flow_cover_constraints_; } + private: void restore_complemented(const std::vector& complemented_variables) { @@ -389,6 +403,7 @@ class knapsack_generation_t { std::vector is_slack_; std::vector knapsack_constraints_; + std::vector flow_cover_constraints_; std::vector is_complemented_; std::vector is_marked_; std::vector workspace_; @@ -473,6 +488,15 @@ class cut_generation_t { const std::vector& xstar, f_t start_time); + // Generate all flow cover cuts + void generate_flow_cover_cuts(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, + f_t start_time); + // Generate clique cuts from conflict graph cliques bool generate_clique_cuts(const lp_problem_t& lp, const simplex_solver_settings_t& settings, diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index cfc120e477..dcfff0d280 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), + flow_cover_cuts(-1), implied_bound_cuts(-1), clique_cuts(-1), strong_chvatal_gomory_cuts(-1), @@ -184,6 +185,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 flow_cover_cuts; // -1 automatic, 0 to disable, >0 to enable flow cover 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 diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index c23b1d27ca..8afe1b461e 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -131,6 +131,7 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_MIP_MIXED_INTEGER_ROUNDING_CUTS, &mip_settings.mir_cuts, -1, 1, -1}, {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_FLOW_COVER_CUTS, &mip_settings.flow_cover_cuts, -1, 1, -1}, {CUOPT_MIP_CLIQUE_CUTS, &mip_settings.clique_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}, diff --git a/cpp/src/mip_heuristics/solver.cu b/cpp/src/mip_heuristics/solver.cu index 0bbf48d95e..8ed5a60ada 100644 --- a/cpp/src/mip_heuristics/solver.cu +++ b/cpp/src/mip_heuristics/solver.cu @@ -342,6 +342,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.flow_cover_cuts = context.settings.flow_cover_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 = diff --git a/cpp/tests/mip/cuts_test.cu b/cpp/tests/mip/cuts_test.cu index 33b5457dfe..fb559e8b0c 100644 --- a/cpp/tests/mip/cuts_test.cu +++ b/cpp/tests/mip/cuts_test.cu @@ -1434,4 +1434,138 @@ TEST(cuts, clique_neos8_phase4_lp_infeasibility_binary_search) EXPECT_EQ(first_infeasible.value(), injected_index); } +// Problem data for Example 9.8 from Wolsey: flow cover cut test. +// +// The mixed 0-1 set X has the knapsack constraint: +// 3*x0 + 3*x1 + 6*x2 - 3*x3 - 5*x4 - 1*x5 <= 4 +// with all binary variables (0-indexed: x0..x5 correspond to 1-indexed x1..x6). +// +// The LP relaxation solution is x* = (1, 0, 2/3, 1, 0, 0). +// With C1={0,2}, C2={3}, lambda=2, L2={4}, the flow cover inequality simplifies to: +// x0 + x2 - x4 <= 1 +// which is violated by x*: 1 + 2/3 - 0 = 5/3 > 1. +// +// We construct an objective that drives the LP relaxation toward this fractional solution. +mps_parser::mps_data_model_t create_flow_cover_problem_example_9_8() +{ + mps_parser::mps_data_model_t problem; + + // 6 binary variables x0..x5 + // Constraint: 3*x0 + 3*x1 + 6*x2 - 3*x3 - 5*x4 - 1*x5 <= 4 + // We also add x3 <= 1 explicitly (it is binary, but coefficients are negative in the knapsack + // row so the solver needs it to identify flow cover structure). + // + // CSR format (1 row): + std::vector offsets = {0, 6}; + std::vector indices = {0, 1, 2, 3, 4, 5}; + std::vector coefficients = {3.0, 3.0, 6.0, -3.0, -5.0, -1.0}; + problem.set_csr_constraint_matrix(coefficients.data(), + coefficients.size(), + indices.data(), + indices.size(), + offsets.data(), + offsets.size()); + + // Constraint bounds: -inf <= sum <= 4 + std::vector lower_bounds = {-std::numeric_limits::infinity()}; + std::vector upper_bounds = {4.0}; + problem.set_constraint_lower_bounds(lower_bounds.data(), lower_bounds.size()); + problem.set_constraint_upper_bounds(upper_bounds.data(), upper_bounds.size()); + + // Variable bounds: all binary [0, 1] + std::vector var_lower_bounds = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; + std::vector var_upper_bounds = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0}; + problem.set_variable_lower_bounds(var_lower_bounds.data(), var_lower_bounds.size()); + problem.set_variable_upper_bounds(var_upper_bounds.data(), var_upper_bounds.size()); + + // Objective: maximize 3*x0 + 6*x2 + 3*x3 (minimize negated) + // This pushes x0=1, x2 as high as possible, x3=1 (opening the RHS). + // LP relaxation gives x0=1, x2=2/3, x3=1 => obj = 3 + 4 + 3 = 10 + // Integer optimum would be x0=1, x2=0, x3=1 or similar. + std::vector objective_coefficients = {-3.0, 0.0, -6.0, -3.0, 0.0, 0.0}; + problem.set_objective_coefficients(objective_coefficients.data(), objective_coefficients.size()); + + // All integer (binary) + std::vector variable_types = {'I', 'I', 'I', 'I', 'I', 'I'}; + problem.set_variable_types(variable_types); + + problem.set_maximize(false); + return problem; +} + +// End-to-end test: solve the Example 9.8 problem with flow cover cuts enabled. +// The flow cover cut x0 + x2 - x4 <= 1 should cut off the fractional LP solution. +TEST(cuts, flow_cover_example_9_8_end_to_end) +{ + const raft::handle_t handle_{}; + mip_solver_settings_t settings; + + auto problem = create_flow_cover_problem_example_9_8(); + + settings.time_limit = 5.0; + settings.max_cut_passes = 10; + settings.flow_cover_cuts = 1; // explicitly enable flow cover cuts + settings.presolver = presolver_t::None; + + mip_solution_t solution = solve_mip(&handle_, problem, settings); + EXPECT_EQ(solution.get_termination_status(), mip_termination_status_t::Optimal); + + // The integer optimal: x0=1, x1=0, x2=0, x3=1, x4=0, x5=0 + // gives 3(1) + 3(0) + 6(0) - 3(1) - 5(0) - 1(0) = 0 <= 4, feasible. + // Objective = -3(1) + 0 + -6(0) + -3(1) + 0 + 0 = -6 + // + // With x0=1, x2=1, x3=1, x4=1: + // 3(1)+3(0)+6(1)-3(1)-5(1)-1(0) = 3+6-3-5 = 1 <= 4, feasible. + // Objective = -3-6-3 = -12. Even better. + // + // With x0=1, x1=1, x2=1, x3=1, x4=1, x5=1: + // 3+3+6-3-5-1 = 3 <= 4, feasible. + // Objective = -3+0-6-3+0+0 = -12 (x1, x4, x5 have zero obj coeff). + // Actually x0=1, x2=1, x3=1, x4=1: obj = -3-6-3 = -12 + double obj_val = solution.get_objective_value(); + EXPECT_LE(obj_val, -6.0 + 1e-3); // at least as good as -6 +} + +// Test that the flow cover cut from Example 9.8 violates the LP relaxation. +// The fractional point is x* = (1, 0, 2/3, 1, 0, 0). +// The flow cover inequality (simplified) is: x0 + x2 - x4 <= 1. +// Violation: 1 + 2/3 - 0 = 5/3 > 1. +TEST(cuts, flow_cover_example_9_8_violation_check) +{ + // Verify the cut violation algebraically + // LP relaxation: x* = (1, 0, 2/3, 1, 0, 0) + std::vector xstar = {1.0, 0.0, 2.0 / 3.0, 1.0, 0.0, 0.0}; + + // Flow cover inequality from Example 9.8: + // y1 + y3 + 1*(1-x1) + 4*(1-x3) <= 7 + 2*x5 + y6 + // With y_j = a_j * x_j: + // 3*x0 + 6*x2 + 1*(1-x0) + 4*(1-x2) <= 7 + 2*x4 + 1*x5 + // 3*x0 + 6*x2 + 1 - x0 + 4 - 4*x2 <= 7 + 2*x4 + x5 + // 2*x0 + 2*x2 + 5 <= 7 + 2*x4 + x5 + // 2*x0 + 2*x2 - 2*x4 - x5 <= 2 + // Dividing is not needed for the violation check. + // + // LHS = 2(1) + 2(2/3) - 2(0) - 0 = 2 + 4/3 = 10/3 + // RHS = 2 + // Violation = 10/3 - 2 = 4/3 + + double lhs = 2.0 * xstar[0] + 2.0 * xstar[2] - 2.0 * xstar[4] - 1.0 * xstar[5]; + double rhs = 2.0; + double violation = lhs - rhs; + + EXPECT_GT(violation, 1.0); // violation is 4/3 ~ 1.333 + EXPECT_NEAR(violation, 4.0 / 3.0, 1e-10); + + // Also verify that any integer feasible solution satisfies the cut. + // x0=1, x2=1, x3=1, x4=1 => lhs = 2+2-2-0 = 2 <= 2. Tight but feasible. + std::vector x_int = {1.0, 0.0, 1.0, 1.0, 1.0, 0.0}; + double lhs_int = 2.0 * x_int[0] + 2.0 * x_int[2] - 2.0 * x_int[4] - 1.0 * x_int[5]; + EXPECT_LE(lhs_int, rhs + 1e-10); + + // x0=1, x2=0, x3=1 => lhs = 2+0-0-0 = 2 <= 2. Feasible. + std::vector x_int2 = {1.0, 0.0, 0.0, 1.0, 0.0, 0.0}; + double lhs_int2 = 2.0 * x_int2[0] + 2.0 * x_int2[2] - 2.0 * x_int2[4] - 1.0 * x_int2[5]; + EXPECT_LE(lhs_int2, rhs + 1e-10); +} + } // namespace cuopt::linear_programming::test From b37aa78f8869d1d7c9640e1b9d4a4bab691c1115 Mon Sep 17 00:00:00 2001 From: Hugo Linsenmaier Date: Mon, 20 Apr 2026 14:03:15 -0700 Subject: [PATCH 02/17] Reuse existing knapsack code --- cpp/src/cuts/cuts.cpp | 307 +++++++++++++++++-------------------- cpp/tests/mip/cuts_test.cu | 159 +++++++++++-------- 2 files changed, 236 insertions(+), 230 deletions(-) diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index 5986bbca6b..fa6c24b29a 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -955,181 +955,180 @@ i_t knapsack_generation_t::generate_flow_cover_cut( inequality_t& cut) { const bool verbose = false; - // Get the row: sum_{j in N1} a_j x_j - sum_{j in N2} a_j x_j <= b - // which represents: sum_{j in N1} y_j <= b + sum_{j in N2} y_j - // where y_j = a_j * x_j for binary x_j. + static_cast(new_slacks); inequality_t inequality(Arow, flow_cover_row, lp.rhs[flow_cover_row]); inequality_t rational_inequality = inequality; if (!rational_coefficients(var_types, inequality, rational_inequality)) { return -1; } inequality = rational_inequality; - // Separate into N1 (positive coefficients) and N2 (negative coefficients) - // N1: j with a_j > 0 (left-hand side flows) - // N2: j with a_j < 0 (right-hand side flows, stored as negative in the constraint) struct flow_var_t { - i_t j; // variable index - f_t a_j; // absolute value of coefficient + i_t j; + f_t a_j; + bool in_n2; }; - std::vector N1, N2; - N1.reserve(inequality.size()); - N2.reserve(inequality.size()); + std::vector flow_vars; + flow_vars.reserve(inequality.size()); - f_t b = inequality.rhs; // RHS of the constraint + std::vector complemented_variables; + complemented_variables.reserve(inequality.size()); + + f_t b = inequality.rhs; + f_t sum_all_a = 0.0; + f_t sum_n2_a = 0.0; for (i_t k = 0; k < inequality.size(); k++) { const i_t j = inequality.index(k); const f_t aj = inequality.coeff(k); if (is_slack_[j]) { continue; } if (aj > 0.0) { - N1.push_back({j, aj}); + flow_vars.push_back({j, aj, false}); + sum_all_a += aj; } else if (aj < 0.0) { - N2.push_back({j, -aj}); // store absolute value + const f_t weight = -aj; + flow_vars.push_back({j, weight, true}); + complemented_variables.push_back(j); + is_complemented_[j] = 1; + sum_all_a += weight; + sum_n2_a += weight; } } - if (N1.empty() || N2.empty()) { return -1; } + if (flow_vars.empty() || sum_n2_a <= 0.0 || sum_n2_a >= sum_all_a) { + restore_complemented(complemented_variables); + return -1; + } - // Solve the separation knapsack problem to find cover C = (C1, C2) - // minimize sum_{j in N1} (1 - x_j*) z_j + sum_{j in N2} (-x_j*) z_j (note: N2 values negative) - // subject to sum_{j in N1} a_j z_j - sum_{j in N2} a_j z_j > b - // z in {0,1}^n - // - // Convert to maximization knapsack (zbar_j = 1 - z_j): - // maximize sum_{j in N1} (1 - x_j*) zbar_j + sum_{j in N2} x_j* zbar_j - // subject to sum_{j in N1} a_j zbar_j + sum_{j in N2} a_j zbar_j <= sum_all_a - b - epsilon - // - // Wait, let's be more careful. The separation problem objective is: - // zeta = sum_{j in N1} (1 - x_j*) z_j - sum_{j in N2} x_j* z_j - // We want to minimize this. A violated flow cover exists when zeta < 1 - sum_{j in N2} x_j*. - // - // Convert to maximize by flipping with zbar = 1 - z: - // maximize sum_{j in N1} (1 - x_j*) zbar_j + sum_{j in N2} x_j* zbar_j - // subject to sum_{j in N1} a_j (1 - zbar_j) - sum_{j in N2} a_j (1 - zbar_j) > b - // which is: sum_N1 a_j - sum_N1 a_j zbar_j - sum_N2 a_j + sum_N2 a_j zbar_j > b - // => -sum_N1 a_j zbar_j + sum_N2 a_j zbar_j > b - sum_N1 a_j + sum_N2 a_j - // => sum_N1 a_j zbar_j - sum_N2 a_j zbar_j < sum_N1 a_j - sum_N2 a_j - b - // - // This is a knapsack: maximize value s.t. sum weights <= capacity - // where weights for N1 items = a_j, weights for N2 items = -a_j (negative weights!) - // - // Negative weights complicate this. Instead, let's use the greedy approach directly - // on the separation problem. We solve: - // minimize sum_{j in N1} (1 - x_j*) z_j - sum_{j in N2} x_j* z_j - // s.t. sum_{j in N1} a_j z_j - sum_{j in N2} a_j z_j > b - // - // Greedy: sort N1 by (1-x_j*)/a_j ascending, N2 by x_j*/a_j descending - // Include all of N2 first (they help the constraint and reduce objective), - // then greedily add N1 items. - - // Greedy separation heuristic - f_t constraint_lhs = 0.0; - f_t objective = 0.0; - - // Start by including all N2 variables in the cover (z_j = 1 for j in N2) - // This contributes -a_j to constraint and -x_j* to objective - std::vector in_cover_n2(N2.size(), true); - for (i_t k = 0; k < static_cast(N2.size()); k++) { - constraint_lhs -= N2[k].a_j; - objective -= xstar[N2[k].j]; - } - - // Sort N1 by (1 - x_j*) / a_j ratio (ascending = cheapest first) - std::vector n1_perm(N1.size()); - std::iota(n1_perm.begin(), n1_perm.end(), 0); - std::sort(n1_perm.begin(), n1_perm.end(), [&](i_t a, i_t idx_b) { - f_t ratio_a = (1.0 - xstar[N1[a].j]) / N1[a].a_j; - f_t ratio_b = (1.0 - xstar[N1[idx_b].j]) / N1[idx_b].a_j; - return ratio_a < ratio_b; - }); + const f_t transformed_beta = b + sum_n2_a; + f_t separation_rhs = sum_all_a - (transformed_beta + 1.0); + if (separation_rhs <= 0.0) { + restore_complemented(complemented_variables); + return -1; + } - // Add N1 variables greedily until constraint > b - std::vector in_cover_n1(N1.size(), false); - for (i_t idx : n1_perm) { - constraint_lhs += N1[idx].a_j; - objective += (1.0 - xstar[N1[idx].j]); - in_cover_n1[idx] = true; - if (constraint_lhs > b) { break; } - } - - if (constraint_lhs <= b) { return -1; } // Could not find a cover - - // Remove N2 variables that are expensive (high x_j*) and not needed - // Sort N2 by x_j*/a_j descending (most expensive first) - std::vector n2_perm(N2.size()); - std::iota(n2_perm.begin(), n2_perm.end(), 0); - std::sort(n2_perm.begin(), n2_perm.end(), [&](i_t a, i_t idx_b) { - f_t ratio_a = xstar[N2[a].j] / N2[a].a_j; - f_t ratio_b = xstar[N2[idx_b].j] / N2[idx_b].a_j; - return ratio_a > ratio_b; - }); + std::vector values; + std::vector weights; + std::vector indices; + std::vector fixed_values; + std::vector fixed_variables; + values.reserve(flow_vars.size()); + weights.reserve(flow_vars.size()); + indices.reserve(flow_vars.size()); + fixed_values.reserve(flow_vars.size()); + fixed_variables.reserve(flow_vars.size()); - for (i_t idx : n2_perm) { - // Try removing this N2 variable from the cover - if (constraint_lhs + N2[idx].a_j > b) { - constraint_lhs += N2[idx].a_j; - objective += xstar[N2[idx].j]; - in_cover_n2[idx] = false; + f_t objective_constant = 0.0; + const f_t x_tol = 1e-5; + for (const auto& flow_var : flow_vars) { + const i_t j = flow_var.j; + const f_t transformed_xstar = flow_var.in_n2 ? 1.0 - xstar[j] : xstar[j]; + complemented_xstar_[j] = transformed_xstar; + const f_t value = std::min(1.0, std::max(0.0, 1.0 - transformed_xstar)); + if (transformed_xstar < x_tol) { + fixed_variables.push_back(j); + fixed_values.push_back(0.0); + separation_rhs -= flow_var.a_j; + continue; } + if (transformed_xstar > 1.0 - x_tol) { + fixed_variables.push_back(j); + fixed_values.push_back(1.0); + objective_constant += value; + continue; + } + objective_constant += value; + values.push_back(value); + weights.push_back(flow_var.a_j); + indices.push_back(j); + } + + if (separation_rhs <= 0.0) { + restore_complemented(complemented_variables); + return -1; + } + + std::vector solution(values.size(), 0.0); + f_t objective = 0.0; + if (!values.empty()) { + objective = solve_knapsack_problem(values, weights, separation_rhs, solution); + } + if (std::isnan(objective)) { + restore_complemented(complemented_variables); + return -1; + } + + const f_t separation_value = -objective + objective_constant; + const f_t tol = 1e-6; + if (separation_value >= 1.0 - tol) { + restore_complemented(complemented_variables); + return -1; } - // Build C1 (subset of N1 in cover) and C2 (subset of N2 in cover) std::vector C1_indices, C2_indices; std::vector C1_coeffs, C2_coeffs; - for (i_t k = 0; k < static_cast(N1.size()); k++) { - if (in_cover_n1[k]) { - C1_indices.push_back(N1[k].j); - C1_coeffs.push_back(N1[k].a_j); - } + std::vector transformed_cover_n2_indices; + std::vector transformed_cover_n2_coeffs; + std::fill(is_marked_.begin(), is_marked_.end(), 0); + + for (i_t k = 0; k < static_cast(solution.size()); k++) { + const i_t j = indices[k]; + if (solution[k] != 0.0) { continue; } + is_marked_[j] = 1; } - for (i_t k = 0; k < static_cast(N2.size()); k++) { - if (in_cover_n2[k]) { - C2_indices.push_back(N2[k].j); - C2_coeffs.push_back(N2[k].a_j); + for (i_t k = 0; k < static_cast(fixed_variables.size()); k++) { + if (fixed_values[k] != 1.0) { continue; } + const i_t j = fixed_variables[k]; + is_marked_[j] = 1; + } + + for (const auto& flow_var : flow_vars) { + const i_t j = flow_var.j; + if (!flow_var.in_n2) { + if (is_marked_[j]) { + C1_indices.push_back(j); + C1_coeffs.push_back(flow_var.a_j); + } + } else if (is_marked_[j]) { + transformed_cover_n2_indices.push_back(j); + transformed_cover_n2_coeffs.push_back(flow_var.a_j); + } else { + C2_indices.push_back(j); + C2_coeffs.push_back(flow_var.a_j); } } - if (C1_indices.empty()) { return -1; } + if (C1_indices.empty()) { + restore_complemented(complemented_variables); + return -1; + } - // Compute lambda = sum_{j in C1} a_j - sum_{j in C2} a_j - b f_t sum_C1 = std::accumulate(C1_coeffs.begin(), C1_coeffs.end(), 0.0); f_t sum_C2 = std::accumulate(C2_coeffs.begin(), C2_coeffs.end(), 0.0); f_t lambda = sum_C1 - sum_C2 - b; + if (lambda <= tol) { + restore_complemented(complemented_variables); + return -1; + } - if (lambda <= 1e-10) { return -1; } - - // Determine L2 = {j in N2 \ C2 : lambda * x_j* < a_j * x_j*} - // = {j in N2 \ C2 : lambda < a_j} (for x_j* > 0) - // Actually from the theory: L2 = {j in N2 \ C2 : lambda * x_j* < y_j*} - // where y_j* = a_j * x_j* for binary variables. - // So: lambda * x_j* < a_j * x_j*, which simplifies to lambda < a_j when x_j* > 0. - // When x_j* = 0, both sides are 0, so j is not in L2. std::vector L2_indices; - std::vector L2_coeffs; - std::vector rest_N2_indices; // N2 \ (C2 union L2) + std::vector rest_N2_indices; std::vector rest_N2_coeffs; - - for (i_t k = 0; k < static_cast(N2.size()); k++) { - if (in_cover_n2[k]) { continue; } // skip C2 - const i_t j = N2[k].j; - const f_t aj = N2[k].a_j; - if (xstar[j] > 1e-10 && lambda < aj) { + for (i_t k = 0; k < static_cast(transformed_cover_n2_indices.size()); k++) { + const i_t j = transformed_cover_n2_indices[k]; + const f_t aj = transformed_cover_n2_coeffs[k]; + if (lambda * xstar[j] < aj * xstar[j] - tol) { L2_indices.push_back(j); - L2_coeffs.push_back(aj); } else { rest_N2_indices.push_back(j); rest_N2_coeffs.push_back(aj); } } - // Evaluate the flow cover inequality violation: - // LHS = sum_{j in C1} y_j* + sum_{j in C1} (a_j - lambda)^+ (1 - x_j*) - // RHS = b + sum_{j in C2} a_j + lambda * sum_{j in L2} x_j* + sum_{j in rest_N2} y_j* f_t lhs = 0.0; for (i_t k = 0; k < static_cast(C1_indices.size()); k++) { const i_t j = C1_indices[k]; const f_t aj = C1_coeffs[k]; - lhs += aj * xstar[j]; // y_j* = a_j * x_j* - lhs += std::max(0.0, aj - lambda) * (1.0 - xstar[j]); // (a_j - lambda)^+ (1 - x_j*) + lhs += aj * xstar[j]; + lhs += std::max(0.0, aj - lambda) * (1.0 - xstar[j]); } f_t rhs = b + sum_C2; @@ -1137,52 +1136,30 @@ i_t knapsack_generation_t::generate_flow_cover_cut( rhs += lambda * xstar[L2_indices[k]]; } for (i_t k = 0; k < static_cast(rest_N2_indices.size()); k++) { - rhs += rest_N2_coeffs[k] * xstar[rest_N2_indices[k]]; // y_j* = a_j * x_j* + rhs += rest_N2_coeffs[k] * xstar[rest_N2_indices[k]]; } - f_t violation = lhs - rhs; - const f_t tol = 1e-6; - if (violation <= tol) { return -1; } + const f_t violation = lhs - rhs; + if (violation <= tol) { + restore_complemented(complemented_variables); + return -1; + } if (verbose) { settings.log.printf( "Flow cover cut row %d violation %e lambda %e\n", flow_cover_row, violation, lambda); } - // Build the cut inequality: LHS <= RHS - // sum_{j in C1} a_j x_j + sum_{j in C1} (a_j - lambda)^+ (1 - x_j) <= b + sum_C2 + lambda * - // sum_{j in L2} x_j + sum_{j in rest_N2} a_j x_j - // - // Rearranging to >= form for cut pool (cut'*x >= rhs): - // -sum_{j in C1} a_j x_j - sum_{j in C1} (a_j - lambda)^+ (1 - x_j) + lambda sum_{j in L2} x_j - // + sum_{j in rest_N2} a_j x_j >= -(b + sum_C2) - // - // Expanding (a_j - lambda)^+ (1 - x_j): - // For j with a_j > lambda: (a_j - lambda)(1 - x_j) = (a_j - lambda) - (a_j - lambda) x_j - // So the x_j coefficient for C1 becomes: -a_j + (a_j - lambda)^+ = -min(a_j, lambda) - // And there's a constant: -sum_{j in C1} (a_j - lambda)^+ - // - // In >= form: - // sum_{j in C1} min(a_j, lambda) x_j + lambda sum_{j in L2} x_j + sum_{j in rest_N2} a_j x_j - // >= sum_{j in C1} (a_j - lambda)^+ - (b + sum_C2) - // = sum_{j in C1} a_j - |C1^+| lambda ... hmm let me just build it directly. - - // Build cut in the form: cut'*x >= cut.rhs - // From: sum_{j in C1} y_j + sum_{j in C1} (a_j-lambda)^+(1-x_j) <= b + sum_C2 + lambda*sum_L2 - // x_j + sum_rest a_j x_j Negate for >= form. cut.clear(); cut.reserve(C1_indices.size() + L2_indices.size() + rest_N2_indices.size()); - f_t constant = 0.0; + f_t lift_sum = 0.0; for (i_t k = 0; k < static_cast(C1_indices.size()); k++) { const i_t j = C1_indices[k]; const f_t aj = C1_coeffs[k]; const f_t lift = std::max(0.0, aj - lambda); - // Coefficient of x_j: -a_j + lift = -(a_j - lift) = -min(a_j, lambda) - // But also the constant from lift*(1-x_j): -lift - f_t coeff = -aj + lift; // = -min(a_j, lambda) - cut.push_back(j, coeff); - constant -= lift; + cut.push_back(j, -aj + lift); + lift_sum += lift; } for (i_t k = 0; k < static_cast(L2_indices.size()); k++) { cut.push_back(L2_indices[k], lambda); @@ -1190,14 +1167,16 @@ i_t knapsack_generation_t::generate_flow_cover_cut( for (i_t k = 0; k < static_cast(rest_N2_indices.size()); k++) { cut.push_back(rest_N2_indices[k], rest_N2_coeffs[k]); } - cut.rhs = -(b + sum_C2) + constant; - - // Verify violation - f_t dot = cut.vector.dot(xstar); - f_t cut_violation = dot - cut.rhs; - if (cut_violation >= -tol) { return -1; } + cut.rhs = -(b + sum_C2) + lift_sum; cut.sort(); + const f_t dot = cut.vector.dot(xstar); + if (dot >= cut.rhs - tol) { + restore_complemented(complemented_variables); + return -1; + } + + restore_complemented(complemented_variables); return 0; } diff --git a/cpp/tests/mip/cuts_test.cu b/cpp/tests/mip/cuts_test.cu index fb559e8b0c..cb3a9598f8 100644 --- a/cpp/tests/mip/cuts_test.cu +++ b/cpp/tests/mip/cuts_test.cu @@ -1493,79 +1493,106 @@ mps_parser::mps_data_model_t create_flow_cover_problem_example_9_8( return problem; } -// End-to-end test: solve the Example 9.8 problem with flow cover cuts enabled. -// The flow cover cut x0 + x2 - x4 <= 1 should cut off the fractional LP solution. -TEST(cuts, flow_cover_example_9_8_end_to_end) -{ - const raft::handle_t handle_{}; - mip_solver_settings_t settings; +struct flow_cover_test_problem_t { + raft::handle_t handle; + dual_simplex::simplex_solver_settings_t settings; + dual_simplex::lp_problem_t lp; + dual_simplex::csr_matrix_t Arow; + std::vector new_slacks; + std::vector var_types; - auto problem = create_flow_cover_problem_example_9_8(); + flow_cover_test_problem_t() : handle(), settings(), lp(&handle, 1, 1, 1), Arow(0, 0, 0) {} +}; - settings.time_limit = 5.0; - settings.max_cut_passes = 10; - settings.flow_cover_cuts = 1; // explicitly enable flow cover cuts - settings.presolver = presolver_t::None; +flow_cover_test_problem_t build_flow_cover_test_problem() +{ + flow_cover_test_problem_t test_problem; + auto model = create_flow_cover_problem_example_9_8(); + auto op_problem = mps_data_model_to_optimization_problem(&test_problem.handle, model); + detail::problem_t mip_problem(op_problem); + dual_simplex::user_problem_t host_problem(op_problem.get_handle_ptr()); + mip_problem.get_host_user_problem(host_problem); - mip_solution_t solution = solve_mip(&handle_, problem, settings); - EXPECT_EQ(solution.get_termination_status(), mip_termination_status_t::Optimal); + dual_simplex::dualize_info_t dualize_info; + dual_simplex::convert_user_problem( + host_problem, test_problem.settings, test_problem.lp, test_problem.new_slacks, dualize_info); + test_problem.var_types = host_problem.var_types; + if (test_problem.lp.num_cols > static_cast(test_problem.var_types.size())) { + test_problem.var_types.resize(test_problem.lp.num_cols, + dual_simplex::variable_type_t::CONTINUOUS); + } + test_problem.lp.A.to_compressed_row(test_problem.Arow); + return test_problem; +} - // The integer optimal: x0=1, x1=0, x2=0, x3=1, x4=0, x5=0 - // gives 3(1) + 3(0) + 6(0) - 3(1) - 5(0) - 1(0) = 0 <= 4, feasible. - // Objective = -3(1) + 0 + -6(0) + -3(1) + 0 + 0 = -6 - // - // With x0=1, x2=1, x3=1, x4=1: - // 3(1)+3(0)+6(1)-3(1)-5(1)-1(0) = 3+6-3-5 = 1 <= 4, feasible. - // Objective = -3-6-3 = -12. Even better. - // - // With x0=1, x1=1, x2=1, x3=1, x4=1, x5=1: - // 3+3+6-3-5-1 = 3 <= 4, feasible. - // Objective = -3+0-6-3+0+0 = -12 (x1, x4, x5 have zero obj coeff). - // Actually x0=1, x2=1, x3=1, x4=1: obj = -3-6-3 = -12 - double obj_val = solution.get_objective_value(); - EXPECT_LE(obj_val, -6.0 + 1e-3); // at least as good as -6 +TEST(cuts, flow_cover_example_9_8_generates_expected_cut) +{ + auto test_problem = build_flow_cover_test_problem(); + std::vector xstar(test_problem.lp.num_cols, 0.0); + xstar[0] = 1.0; + xstar[2] = 2.0 / 3.0; + xstar[3] = 1.0; + + dual_simplex::knapsack_generation_t generator(test_problem.lp, + test_problem.settings, + test_problem.Arow, + test_problem.new_slacks, + test_problem.var_types); + ASSERT_EQ(generator.num_flow_cover_constraints(), 1); + + dual_simplex::inequality_t cut(test_problem.lp.num_cols); + const int status = generator.generate_flow_cover_cut(test_problem.lp, + test_problem.settings, + test_problem.Arow, + test_problem.new_slacks, + test_problem.var_types, + xstar, + generator.get_flow_cover_constraints()[0], + cut); + ASSERT_EQ(status, 0); + + ASSERT_EQ(cut.size(), 4); + EXPECT_EQ(cut.index(0), 0); + EXPECT_EQ(cut.index(1), 2); + EXPECT_EQ(cut.index(2), 4); + EXPECT_EQ(cut.index(3), 5); + EXPECT_NEAR(cut.coeff(0), -2.0, 1e-10); + EXPECT_NEAR(cut.coeff(1), -2.0, 1e-10); + EXPECT_NEAR(cut.coeff(2), 5.0, 1e-10); + EXPECT_NEAR(cut.coeff(3), 1.0, 1e-10); + EXPECT_NEAR(cut.rhs, -2.0, 1e-10); + + const double dot = cut.vector.dot(xstar); + EXPECT_LT(dot, cut.rhs - 1e-6); + EXPECT_NEAR(dot - cut.rhs, -4.0 / 3.0, 1e-10); } -// Test that the flow cover cut from Example 9.8 violates the LP relaxation. -// The fractional point is x* = (1, 0, 2/3, 1, 0, 0). -// The flow cover inequality (simplified) is: x0 + x2 - x4 <= 1. -// Violation: 1 + 2/3 - 0 = 5/3 > 1. -TEST(cuts, flow_cover_example_9_8_violation_check) +TEST(cuts, flow_cover_example_9_8_skips_nonviolated_point) { - // Verify the cut violation algebraically - // LP relaxation: x* = (1, 0, 2/3, 1, 0, 0) - std::vector xstar = {1.0, 0.0, 2.0 / 3.0, 1.0, 0.0, 0.0}; - - // Flow cover inequality from Example 9.8: - // y1 + y3 + 1*(1-x1) + 4*(1-x3) <= 7 + 2*x5 + y6 - // With y_j = a_j * x_j: - // 3*x0 + 6*x2 + 1*(1-x0) + 4*(1-x2) <= 7 + 2*x4 + 1*x5 - // 3*x0 + 6*x2 + 1 - x0 + 4 - 4*x2 <= 7 + 2*x4 + x5 - // 2*x0 + 2*x2 + 5 <= 7 + 2*x4 + x5 - // 2*x0 + 2*x2 - 2*x4 - x5 <= 2 - // Dividing is not needed for the violation check. - // - // LHS = 2(1) + 2(2/3) - 2(0) - 0 = 2 + 4/3 = 10/3 - // RHS = 2 - // Violation = 10/3 - 2 = 4/3 - - double lhs = 2.0 * xstar[0] + 2.0 * xstar[2] - 2.0 * xstar[4] - 1.0 * xstar[5]; - double rhs = 2.0; - double violation = lhs - rhs; - - EXPECT_GT(violation, 1.0); // violation is 4/3 ~ 1.333 - EXPECT_NEAR(violation, 4.0 / 3.0, 1e-10); - - // Also verify that any integer feasible solution satisfies the cut. - // x0=1, x2=1, x3=1, x4=1 => lhs = 2+2-2-0 = 2 <= 2. Tight but feasible. - std::vector x_int = {1.0, 0.0, 1.0, 1.0, 1.0, 0.0}; - double lhs_int = 2.0 * x_int[0] + 2.0 * x_int[2] - 2.0 * x_int[4] - 1.0 * x_int[5]; - EXPECT_LE(lhs_int, rhs + 1e-10); - - // x0=1, x2=0, x3=1 => lhs = 2+0-0-0 = 2 <= 2. Feasible. - std::vector x_int2 = {1.0, 0.0, 0.0, 1.0, 0.0, 0.0}; - double lhs_int2 = 2.0 * x_int2[0] + 2.0 * x_int2[2] - 2.0 * x_int2[4] - 1.0 * x_int2[5]; - EXPECT_LE(lhs_int2, rhs + 1e-10); + auto test_problem = build_flow_cover_test_problem(); + std::vector xstar(test_problem.lp.num_cols, 0.0); + xstar[0] = 1.0; + xstar[2] = 1.0; + xstar[3] = 1.0; + xstar[4] = 1.0; + + dual_simplex::knapsack_generation_t generator(test_problem.lp, + test_problem.settings, + test_problem.Arow, + test_problem.new_slacks, + test_problem.var_types); + ASSERT_EQ(generator.num_flow_cover_constraints(), 1); + + dual_simplex::inequality_t cut(test_problem.lp.num_cols); + const int status = generator.generate_flow_cover_cut(test_problem.lp, + test_problem.settings, + test_problem.Arow, + test_problem.new_slacks, + test_problem.var_types, + xstar, + generator.get_flow_cover_constraints()[0], + cut); + EXPECT_EQ(status, -1); } } // namespace cuopt::linear_programming::test From 6404cf4e77ba17baf29d838c0bd5d62ece904201 Mon Sep 17 00:00:00 2001 From: Hugo Linsenmaier Date: Wed, 22 Apr 2026 16:28:44 +0000 Subject: [PATCH 03/17] Fix compile errors --- cpp/CMakeLists.txt | 17 ++++++++++++++++- .../mip_heuristics/feasibility_jump/fj_cpu.cu | 6 ------ .../mip_heuristics/local_search/local_search.cu | 2 +- 3 files changed, 17 insertions(+), 8 deletions(-) diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index abc6a19ab2..8eecf49066 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -175,7 +175,7 @@ if(${CMAKE_CUDA_COMPILER_VERSION} VERSION_GREATER_EQUAL 12.9) endif() list(APPEND CUOPT_CUDA_FLAGS -Werror=cross-execution-space-call -Wno-deprecated-declarations -Xcompiler=-Werror --default-stream=per-thread) if("${CMAKE_CUDA_HOST_COMPILER}" MATCHES "clang" OR "${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang") - list(APPEND CUOPT_CUDA_FLAGS -Xcompiler=-Wall) + list(APPEND CUOPT_CUDA_FLAGS -Xcompiler=-Wall -Xcompiler=-Wno-unused-function) else() list(APPEND CUOPT_CUDA_FLAGS -Xcompiler=-Wall -Wno-error=non-template-friend) endif() @@ -200,6 +200,21 @@ if(NOT DISABLE_OPENMP) if(OPENMP_FOUND) message(VERBOSE "cuOpt: OpenMP found in ${OpenMP_CXX_INCLUDE_DIRS}") + + # Mixed OpenMP runtimes can yield different omp_lock_t layouts across translation units + # (e.g. clang/libomp for CXX and nvcc-hosted CUDA/libgomp), which later surfaces as + # new/delete type mismatches in inline wrappers such as utilities/omp_helpers.hpp. + if(DEFINED OpenMP_CXX_LIB_NAMES AND DEFINED OpenMP_CUDA_LIB_NAMES AND + NOT "${OpenMP_CXX_LIB_NAMES}" STREQUAL "${OpenMP_CUDA_LIB_NAMES}") + message(FATAL_ERROR + "cuOpt requires a consistent OpenMP runtime across CXX and CUDA compilation units. " + "Detected OpenMP_CXX_LIB_NAMES='${OpenMP_CXX_LIB_NAMES}' but " + "OpenMP_CUDA_LIB_NAMES='${OpenMP_CUDA_LIB_NAMES}'. " + "This configuration can mix libomp and libgomp, producing ABI mismatches such as " + "omp_lock_t size differences and ASAN new-delete-type-mismatch failures during teardown. " + "Configure both toolchains to use the same OpenMP runtime, typically by aligning the " + "CUDA host compiler with the CXX compiler/OpenMP implementation.") + endif() endif() endif() diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu index 4eaa5b6a21..bb07a89577 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu @@ -355,12 +355,6 @@ static void log_regression_features(fj_cpu_climber_t& fj_cpu, // Dynamic runtime features double violated_ratio = (double)fj_cpu.violated_constraints.size() / n_cstrs; - // Compute per-iteration metrics - double nnz_per_move = 0.0; - i_t total_moves = - fj_cpu.n_lift_moves_window + fj_cpu.n_mtm_viol_moves_window + fj_cpu.n_mtm_sat_moves_window; - if (total_moves > 0) { nnz_per_move = (double)fj_cpu.nnz_processed_window / total_moves; } - double eval_intensity = (double)fj_cpu.nnz_processed_window / 1000.0; // Cache and locality metrics diff --git a/cpp/src/mip_heuristics/local_search/local_search.cu b/cpp/src/mip_heuristics/local_search/local_search.cu index da29511d70..b96b48a413 100644 --- a/cpp/src/mip_heuristics/local_search/local_search.cu +++ b/cpp/src/mip_heuristics/local_search/local_search.cu @@ -125,7 +125,7 @@ void local_search_t::start_cpufj_lptopt_scratch_threads( solution_lp, default_weights, default_weights, 0., context.preempt_heuristic_solver_); scratch_cpu_fj_on_lp_opt.fj_cpu->log_prefix = "******* scratch on LP optimal: "; scratch_cpu_fj_on_lp_opt.fj_cpu->improvement_callback = - [this, &population](f_t obj, const std::vector& h_vec, double /*work_units*/) { + [&population](f_t obj, const std::vector& h_vec, double /*work_units*/) { population.add_external_solution(h_vec, obj, solution_origin_t::CPUFJ); if (obj < local_search_best_obj) { CUOPT_LOG_DEBUG("******* New local search best obj %g, best overall %g", From 27159468132f18190d3973b653d3f2341724c606 Mon Sep 17 00:00:00 2001 From: Hugo Linsenmaier Date: Tue, 28 Apr 2026 04:04:43 +0000 Subject: [PATCH 04/17] Implement 0-1 single node flow relaxation --- .../linear_programming/run_mps_files.sh | 38 +- cpp/src/cuts/cuts.cpp | 817 +++++++++++++----- cpp/src/cuts/cuts.hpp | 8 +- cpp/tests/mip/cuts_test.cu | 297 ++++++- 4 files changed, 925 insertions(+), 235 deletions(-) diff --git a/benchmarks/linear_programming/run_mps_files.sh b/benchmarks/linear_programming/run_mps_files.sh index b731eb4697..20eb3af4e2 100755 --- a/benchmarks/linear_programming/run_mps_files.sh +++ b/benchmarks/linear_programming/run_mps_files.sh @@ -24,6 +24,8 @@ # --batch-num : Batch number. This allows to split the work across multiple batches uniformly when resources are limited. # --n-batches : Number of batches # --log-to-console : Log to console +# --recursive : Recursively search for .mps/.MPS/.SIF files under --path +# --cut-mode : Cut family configuration: default, no-cuts, flow-cover-only # # Examples: # # Run all MPS files in /data/lp using 2 GPUs with 1 hour time limit @@ -75,6 +77,8 @@ Optional Arguments: --log-to-console Log to console --model-list File containing a list of models to run --pdlp-tolerances Tolerances for PDLP solver (default: 1e-4) + --recursive Recursively search for .mps/.MPS/.SIF files under --path + --cut-mode MODE Cut family configuration: default, no-cuts, flow-cover-only -h, --help Show this help message and exit Examples: @@ -193,6 +197,16 @@ while [[ $# -gt 0 ]]; do PDLP_TOLERANCES="$2" shift 2 ;; + --recursive) + echo "RECURSIVE: true" + RECURSIVE=true + shift + ;; + --cut-mode) + echo "CUT_MODE: $2" + CUT_MODE="$2" + shift 2 + ;; *) echo "Unknown argument: $1" print_help @@ -221,6 +235,8 @@ N_BATCHES=${N_BATCHES:-1} LOG_TO_CONSOLE=${LOG_TO_CONSOLE:-true} MODEL_LIST=${MODEL_LIST:-} PDLP_TOLERANCES=${PDLP_TOLERANCES:-1e-4} +RECURSIVE=${RECURSIVE:-false} +CUT_MODE=${CUT_MODE:-default} # Validate GPUS_PER_INSTANCE if [[ "$GPUS_PER_INSTANCE" != "1" && "$GPUS_PER_INSTANCE" != "2" ]]; then @@ -228,6 +244,15 @@ if [[ "$GPUS_PER_INSTANCE" != "1" && "$GPUS_PER_INSTANCE" != "2" ]]; then exit 1 fi +case "$CUT_MODE" in + default|no-cuts|flow-cover-only) + ;; + *) + echo "Error: --cut-mode must be one of: default, no-cuts, flow-cover-only" + exit 1 + ;; +esac + # Determine GPU list if [[ -n "$CUDA_VISIBLE_DEVICES" ]]; then IFS=',' read -ra GPU_LIST <<< "$CUDA_VISIBLE_DEVICES" @@ -331,8 +356,12 @@ if [[ -n "$MODEL_LIST" ]]; then exit 1 fi else - # Gather both .mps and .SIF files in the directory - mapfile -t mps_files < <(ls "$MPS_DIR"/*.mps "$MPS_DIR"/*.SIF 2>/dev/null) + if [ "$RECURSIVE" = true ]; then + mapfile -t mps_files < <(find "$MPS_DIR" -type f \( -name "*.mps" -o -name "*.MPS" -o -name "*.SIF" \) | sort) + else + # Gather .mps/.MPS and .SIF files in the directory + mapfile -t mps_files < <(ls "$MPS_DIR"/*.mps "$MPS_DIR"/*.MPS "$MPS_DIR"/*.SIF 2>/dev/null) + fi echo "Found ${#mps_files[@]} .mps and .SIF files in $MPS_DIR" fi @@ -420,6 +449,11 @@ worker() { if [ -n "$METHOD" ]; then args="$args --method $METHOD" fi + if [ "$CUT_MODE" = "no-cuts" ]; then + args="$args --mip-mixed-integer-rounding-cuts 0 --mip-mixed-integer-gomory-cuts 0 --mip-knapsack-cuts 0 --mip-flow-cover-cuts 0 --mip-clique-cuts 0 --mip-implied-bound-cuts 0 --mip-strong-chvatal-gomory-cuts 0 --mip-reduced-cost-strengthening 0" + elif [ "$CUT_MODE" = "flow-cover-only" ]; then + args="$args --mip-mixed-integer-rounding-cuts 0 --mip-mixed-integer-gomory-cuts 0 --mip-knapsack-cuts 0 --mip-flow-cover-cuts 1 --mip-clique-cuts 0 --mip-implied-bound-cuts 0 --mip-strong-chvatal-gomory-cuts 0 --mip-reduced-cost-strengthening 0" + fi if [ -n "$PDLP_TOLERANCES" ]; then args="$args --absolute-primal-tolerance $PDLP_TOLERANCES --absolute-dual-tolerance $PDLP_TOLERANCES --relative-primal-tolerance $PDLP_TOLERANCES --relative-dual-tolerance $PDLP_TOLERANCES --absolute-gap-tolerance $PDLP_TOLERANCES --relative-gap-tolerance $PDLP_TOLERANCES" fi diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index fa6c24b29a..25bbfb8d87 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -12,9 +12,12 @@ #include #include +#include #include #include #include +#include +#include #include #include @@ -899,47 +902,27 @@ knapsack_generation_t::knapsack_generation_t( settings.log.printf("Number of knapsack constraints %d\n", num_knapsack_constraints); #endif - // Identify flow cover constraints: rows with binary variables having both - // positive and negative coefficients (mixed 0-1 structure). - // The constraint has the form: sum_{j in N1} a_j x_j <= b + sum_{j in N2} a_j x_j - // which is stored as: sum_{j in N1} a_j x_j - sum_{j in N2} a_j x_j <= b + // Wolter's separator is applied to every row after attempting to construct a 0-1 + // single-node-flow relaxation. Rows that do not admit the relaxation are rejected + // inside generate_flow_cover_cut. flow_cover_constraints_.reserve(lp.num_rows); for (i_t i = 0; i < lp.num_rows; i++) { - 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 valid = true; - bool has_pos = false; - bool has_neg = false; - i_t num_binaries = 0; - for (i_t p = 0; p < row_len; p++) { - const i_t j = inequality.index(p); - if (is_slack_[j]) { - if (inequality.coeff(p) < 0.0) { - valid = false; - break; - } - continue; - } - if (var_types[j] != variable_type_t::INTEGER || lp.lower[j] != 0.0 || lp.upper[j] != 1.0) { - valid = false; + if (Arow.row_start[i + 1] <= Arow.row_start[i]) { continue; } + flow_cover_constraints_.push_back(i); + bool has_slack = false; + i_t slack_col = -1; + for (i_t p = Arow.row_start[i]; p < Arow.row_start[i + 1]; p++) { + if (is_slack_[Arow.j[p]]) { + has_slack = true; + slack_col = Arow.j[p]; break; } - const f_t aj = inequality.coeff(p); - if (aj > 0.0) { - has_pos = true; - } else if (aj < 0.0) { - has_neg = true; - } - num_binaries++; } - - if (valid && has_pos && has_neg && num_binaries >= 3) { flow_cover_constraints_.push_back(i); } + const bool equality_row = + !has_slack || + (slack_col >= 0 && std::abs(lp.lower[slack_col]) <= 1e-9 && + std::abs(lp.upper[slack_col]) <= 1e-9); + if (equality_row) { flow_cover_constraints_.push_back(-(i + 1)); } } } @@ -949,6 +932,7 @@ i_t knapsack_generation_t::generate_flow_cover_cut( const simplex_solver_settings_t& settings, csr_matrix_t& Arow, const std::vector& new_slacks, + const variable_bounds_t& variable_bounds, const std::vector& var_types, const std::vector& xstar, i_t flow_cover_row, @@ -956,227 +940,632 @@ i_t knapsack_generation_t::generate_flow_cover_cut( { const bool verbose = false; static_cast(new_slacks); - inequality_t inequality(Arow, flow_cover_row, lp.rhs[flow_cover_row]); - inequality_t rational_inequality = inequality; - if (!rational_coefficients(var_types, inequality, rational_inequality)) { return -1; } - inequality = rational_inequality; - - struct flow_var_t { - i_t j; - f_t a_j; + const f_t tol = 1e-6; + const f_t bound_tol = 1e-8; + + struct snf_arc_t { + f_t u; bool in_n2; + i_t x_col; + f_t x_value; + f_t y_const; + i_t y_col; + f_t y_coeff; + f_t y_x_coeff; + f_t y_value; }; - std::vector flow_vars; - flow_vars.reserve(inequality.size()); - std::vector complemented_variables; - complemented_variables.reserve(inequality.size()); + struct snf_candidate_t { + snf_arc_t arc; + f_t b_shift; + f_t distance; + bool absorbs_binary_coeff; + }; + + auto is_binary_variable = [&](i_t j) { + return var_types[j] == variable_type_t::INTEGER && std::abs(lp.lower[j]) <= bound_tol && + std::abs(lp.upper[j] - 1.0) <= bound_tol; + }; - f_t b = inequality.rhs; - f_t sum_all_a = 0.0; - f_t sum_n2_a = 0.0; + auto clamp01 = [&](f_t x) { return std::min(1.0, std::max(0.0, x)); }; - for (i_t k = 0; k < inequality.size(); k++) { - const i_t j = inequality.index(k); - const f_t aj = inequality.coeff(k); + auto arc_x_value = [&](i_t x_col, f_t fixed_value) { + return x_col >= 0 ? clamp01(xstar[x_col]) : fixed_value; + }; + + auto build_arc = [&](f_t u, + bool in_n2, + i_t x_col, + f_t fixed_x, + f_t y_const, + i_t y_col, + f_t y_coeff, + f_t y_x_coeff) { + const f_t x_value = arc_x_value(x_col, fixed_x); + f_t y_value = y_const; + if (y_col >= 0) { y_value += y_coeff * xstar[y_col]; } + if (x_col >= 0) { + y_value += y_x_coeff * xstar[x_col]; + } else { + y_value += y_x_coeff * fixed_x; + } + return snf_arc_t{u, in_n2, x_col, x_value, y_const, y_col, y_coeff, y_x_coeff, y_value}; + }; + + auto solve_strict_kpsnf = [&](const std::vector& values, + const std::vector& weights, + f_t strict_capacity, + std::vector& solution) { + const i_t n = static_cast(weights.size()); + solution.assign(n, 0.0); + if (n == 0) { return static_cast(0.0); } + + const f_t relaxed_capacity = std::max(0.0, strict_capacity - 1e-9); + return greedy_knapsack_problem(values, weights, relaxed_capacity, solution); + }; + + const bool reverse_row_requested = flow_cover_row < 0; + const i_t actual_flow_cover_row = reverse_row_requested ? -flow_cover_row - 1 : flow_cover_row; + inequality_t row(Arow, actual_flow_cover_row, lp.rhs[actual_flow_cover_row]); + i_t slack_count = 0; + f_t slack_coeff = 0.0; + bool negate_row = reverse_row_requested; + const i_t row_len = row.size(); + for (i_t k = 0; k < row_len; k++) { + const i_t j = row.index(k); + if (!is_slack_[j]) { continue; } + slack_count++; + slack_coeff = row.coeff(k); + } + if (slack_count > 1) { return -1; } + if (slack_count == 1) { + if (std::abs(slack_coeff) <= tol) { return -1; } + negate_row = reverse_row_requested ? slack_coeff > 0.0 : slack_coeff < 0.0; + } + + std::vector> continuous_terms; + continuous_terms.reserve(row_len); + std::vector binary_columns; + std::unordered_map binary_coefficients; + f_t b = negate_row ? -row.rhs : row.rhs; + if (!std::isfinite(b)) { return -1; } + + for (i_t k = 0; k < row_len; k++) { + const i_t j = row.index(k); if (is_slack_[j]) { continue; } - if (aj > 0.0) { - flow_vars.push_back({j, aj, false}); - sum_all_a += aj; - } else if (aj < 0.0) { - const f_t weight = -aj; - flow_vars.push_back({j, weight, true}); - complemented_variables.push_back(j); - is_complemented_[j] = 1; - sum_all_a += weight; - sum_n2_a += weight; + const f_t coeff = negate_row ? -row.coeff(k) : row.coeff(k); + if (std::abs(coeff) <= tol) { continue; } + if (is_binary_variable(j)) { + if (binary_coefficients.emplace(j, coeff).second) { + binary_columns.push_back(j); + } else { + binary_coefficients[j] += coeff; + } + } else { + continuous_terms.push_back({j, coeff}); + } + } + + std::vector arcs; + arcs.reserve(row_len); + f_t b_shift = 0.0; + + auto push_candidate = [&](std::vector& candidates, + f_t u, + bool in_n2, + i_t x_col, + f_t fixed_x, + f_t y_const, + i_t y_col, + f_t y_coeff, + f_t y_x_coeff, + f_t shift, + f_t bound_value, + f_t y_star, + bool absorbs_binary_coeff) { + if (u < -bound_tol) { return; } + const f_t capacity = std::max(0.0, u); + if (capacity <= tol) { return; } + const f_t distance = std::abs(bound_value - y_star); + candidates.push_back( + {build_arc(capacity, in_n2, x_col, fixed_x, y_const, y_col, y_coeff, y_x_coeff), + shift, + distance, + absorbs_binary_coeff}); + }; + + for (const auto& [j, c] : continuous_terms) { + const f_t lower_j = lp.lower[j]; + const f_t upper_j = lp.upper[j]; + const f_t y_star = xstar[j]; + std::vector candidates; + + auto add_variable_upper_candidates = [&]() { + const i_t start = variable_bounds.upper_offsets[j]; + const i_t end = variable_bounds.upper_offsets[j + 1]; + for (i_t p = start; p < end; p++) { + const i_t x_col = variable_bounds.upper_variables[p]; + if (!is_binary_variable(x_col)) { continue; } + const f_t gamma = variable_bounds.upper_weights[p]; + const f_t alpha = variable_bounds.upper_biases[p]; + if (!std::isfinite(gamma) || !std::isfinite(alpha) || lower_j <= -inf) { continue; } + const f_t direct_coeff = + binary_coefficients.count(x_col) != 0 ? binary_coefficients[x_col] : 0.0; + const std::array a_values = {direct_coeff, 0.0}; + const i_t num_a_values = std::abs(direct_coeff) > tol ? 2 : 1; + for (i_t h = 0; h < num_a_values; h++) { + const f_t a = a_values[h]; + if (c > 0.0) { + if (c * (lower_j - alpha) >= -bound_tol && + c * (lower_j - alpha) + a >= -bound_tol && c * gamma + a >= -bound_tol) { + push_candidate(candidates, + c * gamma + a, + false, + x_col, + 0.0, + -c * alpha, + j, + c, + a, + c * alpha, + gamma * xstar[x_col] + alpha, + y_star, + std::abs(a) > tol); + } + } else { + if (c * (lower_j - alpha) <= bound_tol && + c * (lower_j - alpha) + a <= bound_tol && c * gamma + a <= bound_tol) { + push_candidate(candidates, + -(c * gamma + a), + true, + x_col, + 0.0, + c * alpha, + j, + -c, + -a, + c * alpha, + gamma * xstar[x_col] + alpha, + y_star, + std::abs(a) > tol); + } + } + } + } + }; + + auto add_variable_lower_candidates = [&]() { + const i_t start = variable_bounds.lower_offsets[j]; + const i_t end = variable_bounds.lower_offsets[j + 1]; + for (i_t p = start; p < end; p++) { + const i_t x_col = variable_bounds.lower_variables[p]; + if (!is_binary_variable(x_col)) { continue; } + const f_t gamma = variable_bounds.lower_weights[p]; + const f_t alpha = variable_bounds.lower_biases[p]; + if (!std::isfinite(gamma) || !std::isfinite(alpha) || upper_j >= inf) { continue; } + const f_t direct_coeff = + binary_coefficients.count(x_col) != 0 ? binary_coefficients[x_col] : 0.0; + const std::array a_values = {direct_coeff, 0.0}; + const i_t num_a_values = std::abs(direct_coeff) > tol ? 2 : 1; + for (i_t h = 0; h < num_a_values; h++) { + const f_t a = a_values[h]; + if (c > 0.0) { + if (c * (upper_j - alpha) <= bound_tol && + c * (upper_j - alpha) + a <= bound_tol && c * gamma + a <= bound_tol) { + push_candidate(candidates, + -(c * gamma + a), + true, + x_col, + 0.0, + c * alpha, + j, + -c, + -a, + c * alpha, + gamma * xstar[x_col] + alpha, + y_star, + std::abs(a) > tol); + } + } else { + if (c * (upper_j - alpha) >= -bound_tol && + c * (upper_j - alpha) + a >= -bound_tol && c * gamma + a >= -bound_tol) { + push_candidate(candidates, + c * gamma + a, + false, + x_col, + 0.0, + -c * alpha, + j, + c, + a, + c * alpha, + gamma * xstar[x_col] + alpha, + y_star, + std::abs(a) > tol); + } + } + } + } + }; + + auto add_simple_bound_candidates = [&]() { + if (lower_j <= -inf || upper_j >= inf) { return; } + const f_t capacity = std::abs(c) * (upper_j - lower_j); + if (capacity <= tol) { return; } + if (c > 0.0) { + push_candidate(candidates, + capacity, + false, + -1, + 1.0, + -c * lower_j, + j, + c, + 0.0, + c * lower_j, + upper_j, + y_star, + false); + push_candidate(candidates, + capacity, + true, + -1, + 1.0, + c * upper_j, + j, + -c, + 0.0, + c * upper_j, + lower_j, + y_star, + false); + } else { + push_candidate(candidates, + capacity, + true, + -1, + 1.0, + c * lower_j, + j, + -c, + 0.0, + c * lower_j, + upper_j, + y_star, + false); + push_candidate(candidates, + capacity, + false, + -1, + 1.0, + -c * upper_j, + j, + c, + 0.0, + c * upper_j, + lower_j, + y_star, + false); + } + }; + + add_variable_upper_candidates(); + add_variable_lower_candidates(); + add_simple_bound_candidates(); + + if (candidates.empty()) { return -1; } + auto best = std::min_element( + candidates.begin(), candidates.end(), [](const snf_candidate_t& a, const snf_candidate_t& b) { + return a.distance < b.distance; + }); + if (best->arc.y_value < -1e-5 || + best->arc.y_value - best->arc.u * best->arc.x_value > 1e-4 * std::max(1.0, best->arc.u)) { + return -1; } + if (best->arc.y_value < 0.0) { best->arc.y_value = 0.0; } + arcs.push_back(best->arc); + b_shift += best->b_shift; + if (best->absorbs_binary_coeff) { binary_coefficients[best->arc.x_col] = 0.0; } } - if (flow_vars.empty() || sum_n2_a <= 0.0 || sum_n2_a >= sum_all_a) { - restore_complemented(complemented_variables); - return -1; + for (i_t j : binary_columns) { + const f_t coeff = binary_coefficients[j]; + if (std::abs(coeff) <= tol) { continue; } + const f_t u = std::abs(coeff); + arcs.push_back(build_arc( + u, coeff < 0.0, j, 0.0, 0.0, -1, 0.0, u)); } - const f_t transformed_beta = b + sum_n2_a; - f_t separation_rhs = sum_all_a - (transformed_beta + 1.0); - if (separation_rhs <= 0.0) { - restore_complemented(complemented_variables); - return -1; - } + if (arcs.empty()) { return -1; } + f_t snf_b = b - b_shift; + + auto compute_structure = [&]() { + bool has_binary_controlled_arc = false; + bool has_n1 = false; + f_t sum_n1_u = 0.0; + for (const auto& arc : arcs) { + if (arc.x_col >= 0) { has_binary_controlled_arc = true; } + if (!arc.in_n2) { + has_n1 = true; + sum_n1_u += arc.u; + } + } + const f_t cover_capacity = -snf_b + sum_n1_u; + return std::make_tuple(has_binary_controlled_arc, has_n1, cover_capacity); + }; + + auto [has_binary_controlled_arc, has_n1, cover_capacity] = compute_structure(); + if (!has_binary_controlled_arc || !has_n1 || cover_capacity <= tol) { return -1; } std::vector values; std::vector weights; - std::vector indices; - std::vector fixed_values; - std::vector fixed_variables; - values.reserve(flow_vars.size()); - weights.reserve(flow_vars.size()); - indices.reserve(flow_vars.size()); - fixed_values.reserve(flow_vars.size()); - fixed_variables.reserve(flow_vars.size()); + std::vector item_to_arc; + values.reserve(arcs.size()); + weights.reserve(arcs.size()); + item_to_arc.reserve(arcs.size()); + for (i_t k = 0; k < static_cast(arcs.size()); k++) { + const auto& arc = arcs[k]; + if (arc.u <= tol) { continue; } + const f_t value = arc.in_n2 ? arc.x_value : 1.0 - arc.x_value; + values.push_back(std::max(0.0, value)); + weights.push_back(arc.u); + item_to_arc.push_back(k); + } + if (values.empty()) { return -1; } - f_t objective_constant = 0.0; - const f_t x_tol = 1e-5; - for (const auto& flow_var : flow_vars) { - const i_t j = flow_var.j; - const f_t transformed_xstar = flow_var.in_n2 ? 1.0 - xstar[j] : xstar[j]; - complemented_xstar_[j] = transformed_xstar; - const f_t value = std::min(1.0, std::max(0.0, 1.0 - transformed_xstar)); - if (transformed_xstar < x_tol) { - fixed_variables.push_back(j); - fixed_values.push_back(0.0); - separation_rhs -= flow_var.a_j; - continue; - } - if (transformed_xstar > 1.0 - x_tol) { - fixed_variables.push_back(j); - fixed_values.push_back(1.0); - objective_constant += value; - continue; + std::vector solution; + const f_t objective = solve_strict_kpsnf(values, weights, cover_capacity, solution); + if (std::isnan(objective)) { return -1; } + static_cast(objective); + + std::vector in_c1(arcs.size(), 0); + std::vector in_c2(arcs.size(), 0); + for (i_t item = 0; item < static_cast(solution.size()); item++) { + const i_t arc_index = item_to_arc[item]; + if (arcs[arc_index].in_n2) { + if (solution[item] > 0.5) { in_c2[arc_index] = 1; } + } else { + if (solution[item] <= 0.5) { in_c1[arc_index] = 1; } } - objective_constant += value; - values.push_back(value); - weights.push_back(flow_var.a_j); - indices.push_back(j); } - if (separation_rhs <= 0.0) { - restore_complemented(complemented_variables); - return -1; + f_t sum_c1_u = 0.0; + f_t sum_c2_u = 0.0; + bool c1_nonempty = false; + for (i_t k = 0; k < static_cast(arcs.size()); k++) { + if (in_c1[k]) { + sum_c1_u += arcs[k].u; + c1_nonempty = true; + } + if (in_c2[k]) { sum_c2_u += arcs[k].u; } } + const f_t lambda = -snf_b + sum_c1_u - sum_c2_u; + if (!c1_nonempty || lambda <= tol) { return -1; } - std::vector solution(values.size(), 0.0); - f_t objective = 0.0; - if (!values.empty()) { - objective = solve_knapsack_problem(values, weights, separation_rhs, solution); - } - if (std::isnan(objective)) { - restore_complemented(complemented_variables); - return -1; - } + auto fractional_part_local = [](f_t value) { + const f_t floor_value = std::floor(value); + return value - floor_value; + }; - const f_t separation_value = -objective + objective_constant; - const f_t tol = 1e-6; - if (separation_value >= 1.0 - tol) { - restore_complemented(complemented_variables); - return -1; - } + auto mir_F = [&](f_t alpha, f_t value) { + const f_t floor_value = std::floor(value); + const f_t frac_value = value - floor_value; + return floor_value + std::max(0.0, frac_value - alpha) / (1.0 - alpha); + }; - std::vector C1_indices, C2_indices; - std::vector C1_coeffs, C2_coeffs; - std::vector transformed_cover_n2_indices; - std::vector transformed_cover_n2_coeffs; - std::fill(is_marked_.begin(), is_marked_.end(), 0); + std::vector ubar_candidates; + auto add_ubar_candidate = [&](f_t candidate) { + if (candidate <= lambda + tol || !std::isfinite(candidate)) { return; } + for (f_t existing : ubar_candidates) { + if (std::abs(existing - candidate) <= 1e-8 * std::max(1.0, std::abs(candidate))) { + return; + } + } + ubar_candidates.push_back(candidate); + }; - for (i_t k = 0; k < static_cast(solution.size()); k++) { - const i_t j = indices[k]; - if (solution[k] != 0.0) { continue; } - is_marked_[j] = 1; - } - for (i_t k = 0; k < static_cast(fixed_variables.size()); k++) { - if (fixed_values[k] != 1.0) { continue; } - const i_t j = fixed_variables[k]; - is_marked_[j] = 1; + f_t max_c1_ltilde2 = -inf; + f_t max_c1 = -inf; + f_t max_u_ge_lambda = -inf; + for (i_t k = 0; k < static_cast(arcs.size()); k++) { + const auto& arc = arcs[k]; + if (arc.u > lambda + tol) { add_ubar_candidate(arc.u); } + if (arc.u >= lambda - tol) { max_u_ge_lambda = std::max(max_u_ge_lambda, arc.u); } + if (!arc.in_n2 && in_c1[k] && arc.u > lambda + tol) { + max_c1 = std::max(max_c1, arc.u); + max_c1_ltilde2 = std::max(max_c1_ltilde2, arc.u); + } + if (arc.in_n2 && !in_c2[k] && arc.u > lambda + tol && + std::min(arc.u, lambda) * arc.x_value <= arc.y_value + tol) { + max_c1_ltilde2 = std::max(max_c1_ltilde2, arc.u); + } } + add_ubar_candidate(max_c1_ltilde2); + add_ubar_candidate(max_c1); + add_ubar_candidate(max_u_ge_lambda + 1.0); + add_ubar_candidate(lambda + 1.0); + if (ubar_candidates.empty()) { return -1; } + + f_t best_violation = 0.0; + f_t best_ubar = 0.0; + std::vector best_in_l1(arcs.size(), 0); + std::vector best_in_l2(arcs.size(), 0); + + auto contribution = [&](const snf_arc_t& arc, + bool arc_in_c1, + bool arc_in_c2, + bool arc_in_l1, + bool arc_in_l2, + f_t ubar) { + const f_t f_pos = mir_F(fractional_part_local(-lambda / ubar), arc.u / ubar); + const f_t f_neg = mir_F(fractional_part_local(-lambda / ubar), -arc.u / ubar); + if (arc_in_c1) { + return arc.y_value + (arc.u + lambda * f_neg) * (1.0 - arc.x_value); + } + if (!arc.in_n2) { + return arc_in_l1 ? arc.y_value - (arc.u - lambda * f_pos) * arc.x_value + : static_cast(0.0); + } + if (arc_in_c2) { + return -arc.u + lambda * f_pos * (1.0 - arc.x_value); + } + return arc_in_l2 ? lambda * f_neg * arc.x_value : -arc.y_value; + }; - for (const auto& flow_var : flow_vars) { - const i_t j = flow_var.j; - if (!flow_var.in_n2) { - if (is_marked_[j]) { - C1_indices.push_back(j); - C1_coeffs.push_back(flow_var.a_j); + for (f_t ubar : ubar_candidates) { + const f_t beta = -lambda / ubar; + const f_t f_beta = fractional_part_local(beta); + if (f_beta < 0.01 || f_beta > 0.95) { continue; } + + std::vector in_l1(arcs.size(), 0); + std::vector in_l2(arcs.size(), 0); + for (i_t k = 0; k < static_cast(arcs.size()); k++) { + const auto& arc = arcs[k]; + const f_t f_pos = mir_F(f_beta, arc.u / ubar); + const f_t f_neg = mir_F(f_beta, -arc.u / ubar); + if (!arc.in_n2 && !in_c1[k] && + arc.y_value - (arc.u - lambda * f_pos) * arc.x_value >= -tol) { + in_l1[k] = 1; + } + if (arc.in_n2 && !in_c2[k] && lambda * f_neg * arc.x_value >= -arc.y_value - tol) { + in_l2[k] = 1; } - } else if (is_marked_[j]) { - transformed_cover_n2_indices.push_back(j); - transformed_cover_n2_coeffs.push_back(flow_var.a_j); - } else { - C2_indices.push_back(j); - C2_coeffs.push_back(flow_var.a_j); } - } - if (C1_indices.empty()) { - restore_complemented(complemented_variables); - return -1; + f_t lhs_value = 0.0; + for (i_t k = 0; k < static_cast(arcs.size()); k++) { + lhs_value += contribution(arcs[k], in_c1[k], in_c2[k], in_l1[k], in_l2[k], ubar); + } + const f_t violation = lhs_value - snf_b; + if (violation > best_violation + tol) { + best_violation = violation; + best_ubar = ubar; + best_in_l1 = std::move(in_l1); + best_in_l2 = std::move(in_l2); + } } - f_t sum_C1 = std::accumulate(C1_coeffs.begin(), C1_coeffs.end(), 0.0); - f_t sum_C2 = std::accumulate(C2_coeffs.begin(), C2_coeffs.end(), 0.0); - f_t lambda = sum_C1 - sum_C2 - b; - if (lambda <= tol) { - restore_complemented(complemented_variables); - return -1; + std::vector sgfci_in_l2(arcs.size(), 0); + f_t sgfci_lhs_value = 0.0; + for (i_t k = 0; k < static_cast(arcs.size()); k++) { + const auto& arc = arcs[k]; + if (in_c1[k]) { + sgfci_lhs_value += arc.y_value + std::max(0.0, arc.u - lambda) * (1.0 - arc.x_value); + } else if (!arc.in_n2) { + continue; + } else if (in_c2[k]) { + sgfci_lhs_value -= arc.u; + } else if (std::min(arc.u, lambda) * arc.x_value <= arc.y_value + tol) { + sgfci_in_l2[k] = 1; + sgfci_lhs_value -= std::min(arc.u, lambda) * arc.x_value; + } else { + sgfci_lhs_value -= arc.y_value; + } } + const f_t sgfci_violation = sgfci_lhs_value - snf_b; + if (best_violation <= tol && sgfci_violation <= tol) { return -1; } + const bool use_cmirfci = best_violation >= sgfci_violation; - std::vector L2_indices; - std::vector rest_N2_indices; - std::vector rest_N2_coeffs; - for (i_t k = 0; k < static_cast(transformed_cover_n2_indices.size()); k++) { - const i_t j = transformed_cover_n2_indices[k]; - const f_t aj = transformed_cover_n2_coeffs[k]; - if (lambda * xstar[j] < aj * xstar[j] - tol) { - L2_indices.push_back(j); + f_t lhs_constant = 0.0; + std::vector lhs_indices; + std::unordered_map lhs_coefficients; + lhs_indices.reserve(arcs.size() * 2); + + auto add_lhs_coeff = [&](i_t j, f_t value) { + if (j < 0 || std::abs(value) <= 1e-12) { return; } + if (lhs_coefficients.emplace(j, value).second) { + lhs_indices.push_back(j); } else { - rest_N2_indices.push_back(j); - rest_N2_coeffs.push_back(aj); + lhs_coefficients[j] += value; } - } + }; - f_t lhs = 0.0; - for (i_t k = 0; k < static_cast(C1_indices.size()); k++) { - const i_t j = C1_indices[k]; - const f_t aj = C1_coeffs[k]; - lhs += aj * xstar[j]; - lhs += std::max(0.0, aj - lambda) * (1.0 - xstar[j]); - } + auto add_y = [&](const snf_arc_t& arc, f_t multiplier) { + lhs_constant += multiplier * arc.y_const; + if (arc.y_col >= 0) { add_lhs_coeff(arc.y_col, multiplier * arc.y_coeff); } + if (arc.x_col >= 0) { + add_lhs_coeff(arc.x_col, multiplier * arc.y_x_coeff); + } else { + lhs_constant += multiplier * arc.y_x_coeff * arc.x_value; + } + }; - f_t rhs = b + sum_C2; - for (i_t k = 0; k < static_cast(L2_indices.size()); k++) { - rhs += lambda * xstar[L2_indices[k]]; - } - for (i_t k = 0; k < static_cast(rest_N2_indices.size()); k++) { - rhs += rest_N2_coeffs[k] * xstar[rest_N2_indices[k]]; - } + auto add_x = [&](const snf_arc_t& arc, f_t multiplier) { + if (arc.x_col >= 0) { + add_lhs_coeff(arc.x_col, multiplier); + } else { + lhs_constant += multiplier * arc.x_value; + } + }; - const f_t violation = lhs - rhs; - if (violation <= tol) { - restore_complemented(complemented_variables); - return -1; - } + auto add_one_minus_x = [&](const snf_arc_t& arc, f_t multiplier) { + lhs_constant += multiplier; + add_x(arc, -multiplier); + }; - if (verbose) { - settings.log.printf( - "Flow cover cut row %d violation %e lambda %e\n", flow_cover_row, violation, lambda); + if (use_cmirfci) { + const f_t f_beta_best = fractional_part_local(-lambda / best_ubar); + for (i_t k = 0; k < static_cast(arcs.size()); k++) { + const auto& arc = arcs[k]; + const f_t f_pos = mir_F(f_beta_best, arc.u / best_ubar); + const f_t f_neg = mir_F(f_beta_best, -arc.u / best_ubar); + if (in_c1[k]) { + add_y(arc, 1.0); + add_one_minus_x(arc, arc.u + lambda * f_neg); + } else if (!arc.in_n2) { + if (best_in_l1[k]) { + add_y(arc, 1.0); + add_x(arc, -(arc.u - lambda * f_pos)); + } + } else if (in_c2[k]) { + lhs_constant -= arc.u; + add_one_minus_x(arc, lambda * f_pos); + } else if (best_in_l2[k]) { + add_x(arc, lambda * f_neg); + } else { + add_y(arc, -1.0); + } + } + } else { + for (i_t k = 0; k < static_cast(arcs.size()); k++) { + const auto& arc = arcs[k]; + if (in_c1[k]) { + add_y(arc, 1.0); + add_one_minus_x(arc, std::max(0.0, arc.u - lambda)); + } else if (!arc.in_n2) { + continue; + } else if (in_c2[k]) { + lhs_constant -= arc.u; + } else if (sgfci_in_l2[k]) { + add_x(arc, -std::min(arc.u, lambda)); + } else { + add_y(arc, -1.0); + } + } } cut.clear(); - cut.reserve(C1_indices.size() + L2_indices.size() + rest_N2_indices.size()); - - f_t lift_sum = 0.0; - for (i_t k = 0; k < static_cast(C1_indices.size()); k++) { - const i_t j = C1_indices[k]; - const f_t aj = C1_coeffs[k]; - const f_t lift = std::max(0.0, aj - lambda); - cut.push_back(j, -aj + lift); - lift_sum += lift; - } - for (i_t k = 0; k < static_cast(L2_indices.size()); k++) { - cut.push_back(L2_indices[k], lambda); - } - for (i_t k = 0; k < static_cast(rest_N2_indices.size()); k++) { - cut.push_back(rest_N2_indices[k], rest_N2_coeffs[k]); + cut.reserve(lhs_indices.size()); + for (i_t j : lhs_indices) { + const f_t coeff = lhs_coefficients[j]; + if (std::abs(coeff) > 1e-9) { cut.push_back(j, -coeff); } } - cut.rhs = -(b + sum_C2) + lift_sum; - + if (cut.size() == 0) { return -1; } + cut.rhs = lhs_constant - snf_b; cut.sort(); + const f_t dot = cut.vector.dot(xstar); - if (dot >= cut.rhs - tol) { - restore_complemented(complemented_variables); - return -1; + if (dot >= cut.rhs - tol) { return -1; } + + if (verbose) { + settings.log.printf("Flow cover row %d violation %e lambda %e ubar %e type %s\n", + flow_cover_row, + use_cmirfci ? best_violation : sgfci_violation, + lambda, + best_ubar, + use_cmirfci ? "c-MIRFCI" : "SGFCI"); } - restore_complemented(complemented_variables); return 0; } @@ -2086,7 +2475,8 @@ bool cut_generation_t::generate_cuts(const lp_problem_t& lp, // Generate Flow Cover cuts if (settings.flow_cover_cuts != 0) { f_t cut_start_time = tic(); - generate_flow_cover_cuts(lp, settings, Arow, new_slacks, var_types, xstar, start_time); + generate_flow_cover_cuts( + lp, settings, Arow, new_slacks, 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("Flow cover cut generation time %.2f seconds\n", cut_generation_time); @@ -2158,6 +2548,7 @@ void cut_generation_t::generate_flow_cover_cuts( const std::vector& new_slacks, const std::vector& var_types, const std::vector& xstar, + variable_bounds_t& variable_bounds, f_t start_time) { if (knapsack_generation_.num_flow_cover_constraints() > 0) { @@ -2165,7 +2556,7 @@ void cut_generation_t::generate_flow_cover_cuts( if (toc(start_time) >= settings.time_limit) { return; } inequality_t cut(lp.num_cols); i_t status = knapsack_generation_.generate_flow_cover_cut( - lp, settings, Arow, new_slacks, var_types, xstar, flow_cover_row, cut); + lp, settings, Arow, new_slacks, variable_bounds, var_types, xstar, flow_cover_row, cut); if (status == 0) { cut_pool_.add_cut(cut_type_t::FLOW_COVER, cut); } } } diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index 464746edc4..0aa74ec534 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -329,6 +329,9 @@ class cut_pool_t { const f_t min_cut_distance_{1e-4}; }; +template +class variable_bounds_t; + template class knapsack_generation_t { public: @@ -351,6 +354,7 @@ class knapsack_generation_t { const simplex_solver_settings_t& settings, csr_matrix_t& Arow, const std::vector& new_slacks, + const variable_bounds_t& variable_bounds, const std::vector& var_types, const std::vector& xstar, i_t flow_cover_row, @@ -415,9 +419,6 @@ class knapsack_generation_t { template class mixed_integer_rounding_cut_t; -template -class variable_bounds_t; - template class cut_generation_t { public: @@ -495,6 +496,7 @@ class cut_generation_t { const std::vector& new_slacks, const std::vector& var_types, const std::vector& xstar, + variable_bounds_t& variable_bounds, f_t start_time); // Generate clique cuts from conflict graph cliques diff --git a/cpp/tests/mip/cuts_test.cu b/cpp/tests/mip/cuts_test.cu index cb3a9598f8..223f6cea6a 100644 --- a/cpp/tests/mip/cuts_test.cu +++ b/cpp/tests/mip/cuts_test.cu @@ -1441,9 +1441,11 @@ TEST(cuts, clique_neos8_phase4_lp_infeasibility_binary_search) // with all binary variables (0-indexed: x0..x5 correspond to 1-indexed x1..x6). // // The LP relaxation solution is x* = (1, 0, 2/3, 1, 0, 0). -// With C1={0,2}, C2={3}, lambda=2, L2={4}, the flow cover inequality simplifies to: -// x0 + x2 - x4 <= 1 -// which is violated by x*: 1 + 2/3 - 0 = 5/3 > 1. +// Wolter's c-MIRFCI separator can generate, for example: +// 2*x0 + 2*x2 - 2*x4 - x5 <= 2 +// which is violated by x*: 2 + 4/3 > 2. +// The always-approximate KPSNF flow-cover selection may choose a different +// violated valid cut, so the tests below check validity rather than exact coefficients. // // We construct an objective that drives the LP relaxation toward this fractional solution. mps_parser::mps_data_model_t create_flow_cover_problem_example_9_8() @@ -1493,6 +1495,63 @@ mps_parser::mps_data_model_t create_flow_cover_problem_example_9_8( return problem; } +mps_parser::mps_data_model_t create_mixed_flow_cover_problem_example_9_8() +{ + mps_parser::mps_data_model_t problem; + + // Variables 0..5 are binary x variables. Variables 6..11 are continuous y variables. + // The rows model: + // y0 + y1 + y2 <= 4 + y3 + y4 + y5 + // yj <= u_j * x_j + // for u = (3, 3, 6, 3, 5, 1). + std::vector offsets = {0, 6, 8, 10, 12, 14, 16, 18}; + std::vector indices = {6, 7, 8, 9, 10, 11, 0, 6, 1, 7, 2, 8, 3, 9, 4, 10, 5, 11}; + std::vector coefficients = {1.0, + 1.0, + 1.0, + -1.0, + -1.0, + -1.0, + -3.0, + 1.0, + -3.0, + 1.0, + -6.0, + 1.0, + -3.0, + 1.0, + -5.0, + 1.0, + -1.0, + 1.0}; + problem.set_csr_constraint_matrix(coefficients.data(), + coefficients.size(), + indices.data(), + indices.size(), + offsets.data(), + offsets.size()); + + std::vector lower_bounds(7, -std::numeric_limits::infinity()); + std::vector upper_bounds = {4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; + problem.set_constraint_lower_bounds(lower_bounds.data(), lower_bounds.size()); + problem.set_constraint_upper_bounds(upper_bounds.data(), upper_bounds.size()); + + std::vector var_lower_bounds(12, 0.0); + std::vector var_upper_bounds = { + 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 3.0, 3.0, 6.0, 3.0, 5.0, 1.0}; + problem.set_variable_lower_bounds(var_lower_bounds.data(), var_lower_bounds.size()); + problem.set_variable_upper_bounds(var_upper_bounds.data(), var_upper_bounds.size()); + + std::vector objective_coefficients(12, 0.0); + problem.set_objective_coefficients(objective_coefficients.data(), objective_coefficients.size()); + + std::vector variable_types = {'I', 'I', 'I', 'I', 'I', 'I', 'C', 'C', 'C', 'C', 'C', 'C'}; + problem.set_variable_types(variable_types); + + problem.set_maximize(false); + return problem; +} + struct flow_cover_test_problem_t { raft::handle_t handle; dual_simplex::simplex_solver_settings_t settings; @@ -1504,10 +1563,10 @@ struct flow_cover_test_problem_t { flow_cover_test_problem_t() : handle(), settings(), lp(&handle, 1, 1, 1), Arow(0, 0, 0) {} }; -flow_cover_test_problem_t build_flow_cover_test_problem() +flow_cover_test_problem_t build_flow_cover_test_problem( + const mps_parser::mps_data_model_t& model) { flow_cover_test_problem_t test_problem; - auto model = create_flow_cover_problem_example_9_8(); auto op_problem = mps_data_model_to_optimization_problem(&test_problem.handle, model); detail::problem_t mip_problem(op_problem); dual_simplex::user_problem_t host_problem(op_problem.get_handle_ptr()); @@ -1525,6 +1584,112 @@ flow_cover_test_problem_t build_flow_cover_test_problem() return test_problem; } +flow_cover_test_problem_t build_flow_cover_test_problem() +{ + return build_flow_cover_test_problem(create_flow_cover_problem_example_9_8()); +} + +bool is_flow_cover_example_9_8_binary_feasible(const std::vector& x) +{ + const double activity = 3.0 * x[0] + 3.0 * x[1] + 6.0 * x[2] - 3.0 * x[3] - 5.0 * x[4] - x[5]; + return activity <= 4.0 + 1e-9; +} + +std::vector mixed_flow_cover_fractional_solution(int num_cols) +{ + std::vector xstar(num_cols, 0.0); + xstar[0] = 1.0; + xstar[2] = 2.0 / 3.0; + xstar[3] = 1.0; + xstar[6] = 3.0; + xstar[8] = 4.0; + xstar[9] = 3.0; + return xstar; +} + +bool mixed_flow_cover_y_feasible(const std::vector& y) +{ + const double activity = y[0] + y[1] + y[2] - y[3] - y[4] - y[5]; + return activity <= 4.0 + 1e-8; +} + +void expect_mixed_flow_cover_cut_valid_at_point(const dual_simplex::inequality_t& cut, + const std::vector& point, + const std::string& label) +{ + EXPECT_GE(cut.vector.dot(point), cut.rhs - 1e-7) << label; +} + +void expect_mixed_flow_cover_cut_valid_at_simple_extreme_points( + const dual_simplex::inequality_t& cut, int num_cols) +{ + const std::vector capacities = {3.0, 3.0, 6.0, 3.0, 5.0, 1.0}; + const std::vector flow_signs = {1.0, 1.0, 1.0, -1.0, -1.0, -1.0}; + int checked_points = 0; + + for (int x_mask = 0; x_mask < 64; x_mask++) { + std::vector y_upper(6, 0.0); + for (int j = 0; j < 6; j++) { + if (((x_mask >> j) & 1) != 0) { y_upper[j] = capacities[j]; } + } + + for (int y_mask = 0; y_mask < 64; y_mask++) { + std::vector y(6, 0.0); + for (int j = 0; j < 6; j++) { + if (((y_mask >> j) & 1) != 0) { y[j] = y_upper[j]; } + } + if (!mixed_flow_cover_y_feasible(y)) { continue; } + + std::vector point(num_cols, 0.0); + for (int j = 0; j < 6; j++) { + point[j] = ((x_mask >> j) & 1) != 0 ? 1.0 : 0.0; + point[6 + j] = y[j]; + } + expect_mixed_flow_cover_cut_valid_at_point( + cut, + point, + "box vertex x_mask=" + std::to_string(x_mask) + " y_mask=" + std::to_string(y_mask)); + checked_points++; + } + + for (int free_j = 0; free_j < 6; free_j++) { + for (int bound_mask = 0; bound_mask < 32; bound_mask++) { + std::vector y(6, 0.0); + int bit = 0; + for (int j = 0; j < 6; j++) { + if (j == free_j) { continue; } + if (((bound_mask >> bit) & 1) != 0) { y[j] = y_upper[j]; } + bit++; + } + + double fixed_activity = 0.0; + for (int j = 0; j < 6; j++) { + if (j != free_j) { fixed_activity += flow_signs[j] * y[j]; } + } + + const double y_free = (4.0 - fixed_activity) / flow_signs[free_j]; + if (y_free < -1e-8 || y_free > y_upper[free_j] + 1e-8) { continue; } + y[free_j] = std::max(0.0, std::min(y_upper[free_j], y_free)); + if (!mixed_flow_cover_y_feasible(y)) { continue; } + + std::vector point(num_cols, 0.0); + for (int j = 0; j < 6; j++) { + point[j] = ((x_mask >> j) & 1) != 0 ? 1.0 : 0.0; + point[6 + j] = y[j]; + } + expect_mixed_flow_cover_cut_valid_at_point( + cut, + point, + "flow-tight vertex x_mask=" + std::to_string(x_mask) + + " free_j=" + std::to_string(free_j) + " bound_mask=" + std::to_string(bound_mask)); + checked_points++; + } + } + } + + EXPECT_GT(checked_points, 0); +} + TEST(cuts, flow_cover_example_9_8_generates_expected_cut) { auto test_problem = build_flow_cover_test_problem(); @@ -1538,6 +1703,11 @@ TEST(cuts, flow_cover_example_9_8_generates_expected_cut) test_problem.Arow, test_problem.new_slacks, test_problem.var_types); + dual_simplex::variable_bounds_t variable_bounds(test_problem.lp, + test_problem.settings, + test_problem.var_types, + test_problem.Arow, + test_problem.new_slacks); ASSERT_EQ(generator.num_flow_cover_constraints(), 1); dual_simplex::inequality_t cut(test_problem.lp.num_cols); @@ -1545,26 +1715,74 @@ TEST(cuts, flow_cover_example_9_8_generates_expected_cut) test_problem.settings, test_problem.Arow, test_problem.new_slacks, + variable_bounds, test_problem.var_types, xstar, generator.get_flow_cover_constraints()[0], cut); ASSERT_EQ(status, 0); - ASSERT_EQ(cut.size(), 4); - EXPECT_EQ(cut.index(0), 0); - EXPECT_EQ(cut.index(1), 2); - EXPECT_EQ(cut.index(2), 4); - EXPECT_EQ(cut.index(3), 5); - EXPECT_NEAR(cut.coeff(0), -2.0, 1e-10); - EXPECT_NEAR(cut.coeff(1), -2.0, 1e-10); - EXPECT_NEAR(cut.coeff(2), 5.0, 1e-10); - EXPECT_NEAR(cut.coeff(3), 1.0, 1e-10); - EXPECT_NEAR(cut.rhs, -2.0, 1e-10); - const double dot = cut.vector.dot(xstar); EXPECT_LT(dot, cut.rhs - 1e-6); - EXPECT_NEAR(dot - cut.rhs, -4.0 / 3.0, 1e-10); + + int feasible_assignments = 0; + for (int mask = 0; mask < 64; mask++) { + std::vector point(test_problem.lp.num_cols, 0.0); + for (int j = 0; j < 6; j++) { + point[j] = ((mask >> j) & 1) != 0 ? 1.0 : 0.0; + } + if (!is_flow_cover_example_9_8_binary_feasible(point)) { continue; } + EXPECT_GE(cut.vector.dot(point), cut.rhs - 1e-7) << "mask=" << mask; + feasible_assignments++; + } + EXPECT_GT(feasible_assignments, 0); +} + +TEST(cuts, flow_cover_example_9_8_generated_cut_valid_for_binary_feasible_points) +{ + auto test_problem = build_flow_cover_test_problem(); + std::vector xstar(test_problem.lp.num_cols, 0.0); + xstar[0] = 1.0; + xstar[2] = 2.0 / 3.0; + xstar[3] = 1.0; + + dual_simplex::knapsack_generation_t generator(test_problem.lp, + test_problem.settings, + test_problem.Arow, + test_problem.new_slacks, + test_problem.var_types); + dual_simplex::variable_bounds_t variable_bounds(test_problem.lp, + test_problem.settings, + test_problem.var_types, + test_problem.Arow, + test_problem.new_slacks); + ASSERT_EQ(generator.num_flow_cover_constraints(), 1); + + dual_simplex::inequality_t cut(test_problem.lp.num_cols); + const int status = generator.generate_flow_cover_cut(test_problem.lp, + test_problem.settings, + test_problem.Arow, + test_problem.new_slacks, + variable_bounds, + test_problem.var_types, + xstar, + generator.get_flow_cover_constraints()[0], + cut); + ASSERT_EQ(status, 0); + EXPECT_LT(cut.vector.dot(xstar), cut.rhs - 1e-6); + + int feasible_assignments = 0; + for (int mask = 0; mask < 64; mask++) { + std::vector point(test_problem.lp.num_cols, 0.0); + for (int j = 0; j < 6; j++) { + point[j] = ((mask >> j) & 1) != 0 ? 1.0 : 0.0; + } + if (!is_flow_cover_example_9_8_binary_feasible(point)) { continue; } + EXPECT_GE(cut.vector.dot(point), cut.rhs - 1e-7) << "mask=" << mask; + feasible_assignments++; + } + + EXPECT_GT(feasible_assignments, 0); } TEST(cuts, flow_cover_example_9_8_skips_nonviolated_point) @@ -1581,6 +1799,11 @@ TEST(cuts, flow_cover_example_9_8_skips_nonviolated_point) test_problem.Arow, test_problem.new_slacks, test_problem.var_types); + dual_simplex::variable_bounds_t variable_bounds(test_problem.lp, + test_problem.settings, + test_problem.var_types, + test_problem.Arow, + test_problem.new_slacks); ASSERT_EQ(generator.num_flow_cover_constraints(), 1); dual_simplex::inequality_t cut(test_problem.lp.num_cols); @@ -1588,6 +1811,7 @@ TEST(cuts, flow_cover_example_9_8_skips_nonviolated_point) test_problem.settings, test_problem.Arow, test_problem.new_slacks, + variable_bounds, test_problem.var_types, xstar, generator.get_flow_cover_constraints()[0], @@ -1595,4 +1819,43 @@ TEST(cuts, flow_cover_example_9_8_skips_nonviolated_point) EXPECT_EQ(status, -1); } +TEST(cuts, flow_cover_mixed_single_node_flow_generates_valid_cut) +{ + auto test_problem = build_flow_cover_test_problem(create_mixed_flow_cover_problem_example_9_8()); + const std::vector xstar = mixed_flow_cover_fractional_solution(test_problem.lp.num_cols); + + dual_simplex::knapsack_generation_t generator(test_problem.lp, + test_problem.settings, + test_problem.Arow, + test_problem.new_slacks, + test_problem.var_types); + dual_simplex::variable_bounds_t variable_bounds(test_problem.lp, + test_problem.settings, + test_problem.var_types, + test_problem.Arow, + test_problem.new_slacks); + ASSERT_GT(generator.num_flow_cover_constraints(), 0); + + int generated_cuts = 0; + for (const int flow_cover_row : generator.get_flow_cover_constraints()) { + dual_simplex::inequality_t cut(test_problem.lp.num_cols); + const int status = generator.generate_flow_cover_cut(test_problem.lp, + test_problem.settings, + test_problem.Arow, + test_problem.new_slacks, + variable_bounds, + test_problem.var_types, + xstar, + flow_cover_row, + cut); + if (status != 0) { continue; } + + EXPECT_LT(cut.vector.dot(xstar), cut.rhs - 1e-6) << "row=" << flow_cover_row; + expect_mixed_flow_cover_cut_valid_at_simple_extreme_points(cut, test_problem.lp.num_cols); + generated_cuts++; + } + + EXPECT_GT(generated_cuts, 0); +} + } // namespace cuopt::linear_programming::test From a2771c2254a15fdc07b1cc7fcd84cea1eb5f6e6d Mon Sep 17 00:00:00 2001 From: Hugo Linsenmaier Date: Wed, 29 Apr 2026 18:40:56 +0000 Subject: [PATCH 05/17] Implement variable bound extension --- cpp/src/cuts/cuts.cpp | 481 +++++++++++++++--- .../presolve/third_party_presolve.cpp | 1 - cpp/tests/mip/cuts_test.cu | 63 ++- 3 files changed, 472 insertions(+), 73 deletions(-) diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index 25bbfb8d87..ae47921358 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -15,6 +15,7 @@ #include #include #include +#include #include #include #include @@ -34,6 +35,44 @@ namespace { enum class clique_cut_build_status_t : int8_t { NO_CUT = 0, CUT_ADDED = 1, INFEASIBLE = 2 }; +enum flow_cover_cut_status_t : int { + FLOW_COVER_CUT_FOUND = 0, + FLOW_COVER_REJECT_SLACK = -1, + FLOW_COVER_REJECT_RHS = -2, + FLOW_COVER_REJECT_CONTINUOUS_TERM = -3, + FLOW_COVER_REJECT_RELAXATION = -4, + FLOW_COVER_REJECT_NO_ARCS = -5, + FLOW_COVER_REJECT_NO_COVER = -6, + FLOW_COVER_REJECT_NO_ITEMS = -7, + FLOW_COVER_REJECT_KPSNF = -8, + FLOW_COVER_REJECT_LAMBDA = -9, + FLOW_COVER_REJECT_UBAR = -10, + FLOW_COVER_REJECT_NO_VIOLATION = -11, + FLOW_COVER_REJECT_EMPTY_CUT = -12, + FLOW_COVER_REJECT_MAPPED_CUT = -13, +}; + +const char* flow_cover_status_name(int status) +{ + switch (status) { + case FLOW_COVER_CUT_FOUND: return "cut_found"; + case FLOW_COVER_REJECT_SLACK: return "bad_slack"; + case FLOW_COVER_REJECT_RHS: return "bad_rhs"; + case FLOW_COVER_REJECT_CONTINUOUS_TERM: return "continuous_term_no_candidate"; + case FLOW_COVER_REJECT_RELAXATION: return "invalid_snf_relaxation"; + case FLOW_COVER_REJECT_NO_ARCS: return "no_arcs"; + case FLOW_COVER_REJECT_NO_COVER: return "no_cover_capacity"; + case FLOW_COVER_REJECT_NO_ITEMS: return "no_kpsnf_items"; + case FLOW_COVER_REJECT_KPSNF: return "kpsnf_failed"; + case FLOW_COVER_REJECT_LAMBDA: return "no_positive_lambda"; + case FLOW_COVER_REJECT_UBAR: return "no_ubar_candidate"; + case FLOW_COVER_REJECT_NO_VIOLATION: return "no_cmir_or_sgfci_violation"; + case FLOW_COVER_REJECT_EMPTY_CUT: return "empty_mapped_cut"; + case FLOW_COVER_REJECT_MAPPED_CUT: return "mapped_cut_not_violated"; + default: return "unknown"; + } +} + #if DEBUG_CLIQUE_CUTS #define CLIQUE_CUTS_DEBUG(...) \ do { \ @@ -914,14 +953,13 @@ knapsack_generation_t::knapsack_generation_t( for (i_t p = Arow.row_start[i]; p < Arow.row_start[i + 1]; p++) { if (is_slack_[Arow.j[p]]) { has_slack = true; - slack_col = Arow.j[p]; + slack_col = Arow.j[p]; break; } } const bool equality_row = - !has_slack || - (slack_col >= 0 && std::abs(lp.lower[slack_col]) <= 1e-9 && - std::abs(lp.upper[slack_col]) <= 1e-9); + !has_slack || (slack_col >= 0 && std::abs(lp.lower[slack_col]) <= 1e-9 && + std::abs(lp.upper[slack_col]) <= 1e-9); if (equality_row) { flow_cover_constraints_.push_back(-(i + 1)); } } } @@ -938,7 +976,11 @@ i_t knapsack_generation_t::generate_flow_cover_cut( i_t flow_cover_row, inequality_t& cut) { - const bool verbose = false; + const bool verbose = false; + const bool flow_cover_diag_verbose = std::getenv("CUOPT_FLOW_COVER_DIAG_VERBOSE") != nullptr; + const char* flow_cover_diag_row_env = std::getenv("CUOPT_FLOW_COVER_DIAG_ROW"); + const i_t flow_cover_diag_row = + flow_cover_diag_row_env != nullptr ? static_cast(std::atoi(flow_cover_diag_row_env)) : -1; static_cast(new_slacks); const f_t tol = 1e-6; const f_t bound_tol = 1e-8; @@ -1005,7 +1047,7 @@ i_t knapsack_generation_t::generate_flow_cover_cut( }; const bool reverse_row_requested = flow_cover_row < 0; - const i_t actual_flow_cover_row = reverse_row_requested ? -flow_cover_row - 1 : flow_cover_row; + const i_t actual_flow_cover_row = reverse_row_requested ? -flow_cover_row - 1 : flow_cover_row; inequality_t row(Arow, actual_flow_cover_row, lp.rhs[actual_flow_cover_row]); i_t slack_count = 0; f_t slack_coeff = 0.0; @@ -1017,9 +1059,9 @@ i_t knapsack_generation_t::generate_flow_cover_cut( slack_count++; slack_coeff = row.coeff(k); } - if (slack_count > 1) { return -1; } + if (slack_count > 1) { return FLOW_COVER_REJECT_SLACK; } if (slack_count == 1) { - if (std::abs(slack_coeff) <= tol) { return -1; } + if (std::abs(slack_coeff) <= tol) { return FLOW_COVER_REJECT_SLACK; } negate_row = reverse_row_requested ? slack_coeff > 0.0 : slack_coeff < 0.0; } @@ -1028,7 +1070,7 @@ i_t knapsack_generation_t::generate_flow_cover_cut( std::vector binary_columns; std::unordered_map binary_coefficients; f_t b = negate_row ? -row.rhs : row.rhs; - if (!std::isfinite(b)) { return -1; } + if (!std::isfinite(b)) { return FLOW_COVER_REJECT_RHS; } for (i_t k = 0; k < row_len; k++) { const i_t j = row.index(k); @@ -1096,8 +1138,8 @@ i_t knapsack_generation_t::generate_flow_cover_cut( for (i_t h = 0; h < num_a_values; h++) { const f_t a = a_values[h]; if (c > 0.0) { - if (c * (lower_j - alpha) >= -bound_tol && - c * (lower_j - alpha) + a >= -bound_tol && c * gamma + a >= -bound_tol) { + if (c * (lower_j - alpha) >= -bound_tol && c * (lower_j - alpha) + a >= -bound_tol && + c * gamma + a >= -bound_tol) { push_candidate(candidates, c * gamma + a, false, @@ -1113,8 +1155,8 @@ i_t knapsack_generation_t::generate_flow_cover_cut( std::abs(a) > tol); } } else { - if (c * (lower_j - alpha) <= bound_tol && - c * (lower_j - alpha) + a <= bound_tol && c * gamma + a <= bound_tol) { + if (c * (lower_j - alpha) <= bound_tol && c * (lower_j - alpha) + a <= bound_tol && + c * gamma + a <= bound_tol) { push_candidate(candidates, -(c * gamma + a), true, @@ -1150,8 +1192,8 @@ i_t knapsack_generation_t::generate_flow_cover_cut( for (i_t h = 0; h < num_a_values; h++) { const f_t a = a_values[h]; if (c > 0.0) { - if (c * (upper_j - alpha) <= bound_tol && - c * (upper_j - alpha) + a <= bound_tol && c * gamma + a <= bound_tol) { + if (c * (upper_j - alpha) <= bound_tol && c * (upper_j - alpha) + a <= bound_tol && + c * gamma + a <= bound_tol) { push_candidate(candidates, -(c * gamma + a), true, @@ -1167,8 +1209,8 @@ i_t knapsack_generation_t::generate_flow_cover_cut( std::abs(a) > tol); } } else { - if (c * (upper_j - alpha) >= -bound_tol && - c * (upper_j - alpha) + a >= -bound_tol && c * gamma + a >= -bound_tol) { + if (c * (upper_j - alpha) >= -bound_tol && c * (upper_j - alpha) + a >= -bound_tol && + c * gamma + a >= -bound_tol) { push_candidate(candidates, c * gamma + a, false, @@ -1253,14 +1295,63 @@ i_t knapsack_generation_t::generate_flow_cover_cut( add_variable_lower_candidates(); add_simple_bound_candidates(); - if (candidates.empty()) { return -1; } + if (candidates.empty()) { + if (flow_cover_diag_verbose) { + const i_t upper_start = variable_bounds.upper_offsets[j]; + const i_t upper_end = variable_bounds.upper_offsets[j + 1]; + const i_t lower_start = variable_bounds.lower_offsets[j]; + const i_t lower_end = variable_bounds.lower_offsets[j + 1]; + settings.log.printf( + "Flow cover row %d rejected: continuous col=%d coeff=%e lower=%e upper=%e " + "type=%d upper_vb=%d lower_vb=%d\n", + flow_cover_row, + j, + c, + lower_j, + upper_j, + static_cast(var_types[j]), + upper_end - upper_start, + lower_end - lower_start); + if (flow_cover_diag_row == flow_cover_row) { + for (i_t p = upper_start; p < upper_end; p++) { + const i_t x_col = variable_bounds.upper_variables[p]; + const f_t direct_coeff = + binary_coefficients.count(x_col) != 0 ? binary_coefficients[x_col] : 0.0; + settings.log.printf( + " upper candidate bound_var=%d gamma=%e alpha=%e direct_coeff=%e " + "binary=%d xstar=%e\n", + x_col, + variable_bounds.upper_weights[p], + variable_bounds.upper_biases[p], + direct_coeff, + static_cast(is_binary_variable(x_col)), + x_col >= 0 ? xstar[x_col] : 0.0); + } + for (i_t p = lower_start; p < lower_end; p++) { + const i_t x_col = variable_bounds.lower_variables[p]; + const f_t direct_coeff = + binary_coefficients.count(x_col) != 0 ? binary_coefficients[x_col] : 0.0; + settings.log.printf( + " lower candidate bound_var=%d gamma=%e alpha=%e direct_coeff=%e " + "binary=%d xstar=%e\n", + x_col, + variable_bounds.lower_weights[p], + variable_bounds.lower_biases[p], + direct_coeff, + static_cast(is_binary_variable(x_col)), + x_col >= 0 ? xstar[x_col] : 0.0); + } + } + } + return FLOW_COVER_REJECT_CONTINUOUS_TERM; + } auto best = std::min_element( candidates.begin(), candidates.end(), [](const snf_candidate_t& a, const snf_candidate_t& b) { return a.distance < b.distance; }); - if (best->arc.y_value < -1e-5 || - best->arc.y_value - best->arc.u * best->arc.x_value > 1e-4 * std::max(1.0, best->arc.u)) { - return -1; + if (best->arc.y_value < -1e-5 || best->arc.y_value - best->arc.u * best->arc.x_value > + 1e-4 * std::max(1.0, best->arc.u)) { + return FLOW_COVER_REJECT_RELAXATION; } if (best->arc.y_value < 0.0) { best->arc.y_value = 0.0; } arcs.push_back(best->arc); @@ -1272,17 +1363,16 @@ i_t knapsack_generation_t::generate_flow_cover_cut( const f_t coeff = binary_coefficients[j]; if (std::abs(coeff) <= tol) { continue; } const f_t u = std::abs(coeff); - arcs.push_back(build_arc( - u, coeff < 0.0, j, 0.0, 0.0, -1, 0.0, u)); + arcs.push_back(build_arc(u, coeff < 0.0, j, 0.0, 0.0, -1, 0.0, u)); } - if (arcs.empty()) { return -1; } + if (arcs.empty()) { return FLOW_COVER_REJECT_NO_ARCS; } f_t snf_b = b - b_shift; auto compute_structure = [&]() { bool has_binary_controlled_arc = false; - bool has_n1 = false; - f_t sum_n1_u = 0.0; + bool has_n1 = false; + f_t sum_n1_u = 0.0; for (const auto& arc : arcs) { if (arc.x_col >= 0) { has_binary_controlled_arc = true; } if (!arc.in_n2) { @@ -1295,7 +1385,9 @@ i_t knapsack_generation_t::generate_flow_cover_cut( }; auto [has_binary_controlled_arc, has_n1, cover_capacity] = compute_structure(); - if (!has_binary_controlled_arc || !has_n1 || cover_capacity <= tol) { return -1; } + if (!has_binary_controlled_arc || !has_n1 || cover_capacity <= tol) { + return FLOW_COVER_REJECT_NO_COVER; + } std::vector values; std::vector weights; @@ -1311,11 +1403,11 @@ i_t knapsack_generation_t::generate_flow_cover_cut( weights.push_back(arc.u); item_to_arc.push_back(k); } - if (values.empty()) { return -1; } + if (values.empty()) { return FLOW_COVER_REJECT_NO_ITEMS; } std::vector solution; const f_t objective = solve_strict_kpsnf(values, weights, cover_capacity, solution); - if (std::isnan(objective)) { return -1; } + if (std::isnan(objective)) { return FLOW_COVER_REJECT_KPSNF; } static_cast(objective); std::vector in_c1(arcs.size(), 0); @@ -1329,8 +1421,8 @@ i_t knapsack_generation_t::generate_flow_cover_cut( } } - f_t sum_c1_u = 0.0; - f_t sum_c2_u = 0.0; + f_t sum_c1_u = 0.0; + f_t sum_c2_u = 0.0; bool c1_nonempty = false; for (i_t k = 0; k < static_cast(arcs.size()); k++) { if (in_c1[k]) { @@ -1340,7 +1432,18 @@ i_t knapsack_generation_t::generate_flow_cover_cut( if (in_c2[k]) { sum_c2_u += arcs[k].u; } } const f_t lambda = -snf_b + sum_c1_u - sum_c2_u; - if (!c1_nonempty || lambda <= tol) { return -1; } + if (!c1_nonempty || lambda <= tol) { + if (flow_cover_diag_verbose) { + settings.log.printf( + "Flow cover row %d rejected: lambda=%e snf_b=%e cover_capacity=%e arcs=%d\n", + flow_cover_row, + lambda, + snf_b, + cover_capacity, + static_cast(arcs.size())); + } + return FLOW_COVER_REJECT_LAMBDA; + } auto fractional_part_local = [](f_t value) { const f_t floor_value = std::floor(value); @@ -1364,8 +1467,8 @@ i_t knapsack_generation_t::generate_flow_cover_cut( ubar_candidates.push_back(candidate); }; - f_t max_c1_ltilde2 = -inf; - f_t max_c1 = -inf; + f_t max_c1_ltilde2 = -inf; + f_t max_c1 = -inf; f_t max_u_ge_lambda = -inf; for (i_t k = 0; k < static_cast(arcs.size()); k++) { const auto& arc = arcs[k]; @@ -1384,7 +1487,16 @@ i_t knapsack_generation_t::generate_flow_cover_cut( add_ubar_candidate(max_c1); add_ubar_candidate(max_u_ge_lambda + 1.0); add_ubar_candidate(lambda + 1.0); - if (ubar_candidates.empty()) { return -1; } + if (ubar_candidates.empty()) { + if (flow_cover_diag_verbose) { + settings.log.printf("Flow cover row %d rejected: no ubar lambda=%e snf_b=%e arcs=%d\n", + flow_cover_row, + lambda, + snf_b, + static_cast(arcs.size())); + } + return FLOW_COVER_REJECT_UBAR; + } f_t best_violation = 0.0; f_t best_ubar = 0.0; @@ -1399,16 +1511,12 @@ i_t knapsack_generation_t::generate_flow_cover_cut( f_t ubar) { const f_t f_pos = mir_F(fractional_part_local(-lambda / ubar), arc.u / ubar); const f_t f_neg = mir_F(fractional_part_local(-lambda / ubar), -arc.u / ubar); - if (arc_in_c1) { - return arc.y_value + (arc.u + lambda * f_neg) * (1.0 - arc.x_value); - } + if (arc_in_c1) { return arc.y_value + (arc.u + lambda * f_neg) * (1.0 - arc.x_value); } if (!arc.in_n2) { return arc_in_l1 ? arc.y_value - (arc.u - lambda * f_pos) * arc.x_value : static_cast(0.0); } - if (arc_in_c2) { - return -arc.u + lambda * f_pos * (1.0 - arc.x_value); - } + if (arc_in_c2) { return -arc.u + lambda * f_pos * (1.0 - arc.x_value); } return arc_in_l2 ? lambda * f_neg * arc.x_value : -arc.y_value; }; @@ -1423,8 +1531,7 @@ i_t knapsack_generation_t::generate_flow_cover_cut( const auto& arc = arcs[k]; const f_t f_pos = mir_F(f_beta, arc.u / ubar); const f_t f_neg = mir_F(f_beta, -arc.u / ubar); - if (!arc.in_n2 && !in_c1[k] && - arc.y_value - (arc.u - lambda * f_pos) * arc.x_value >= -tol) { + if (!arc.in_n2 && !in_c1[k] && arc.y_value - (arc.u - lambda * f_pos) * arc.x_value >= -tol) { in_l1[k] = 1; } if (arc.in_n2 && !in_c2[k] && lambda * f_neg * arc.x_value >= -arc.y_value - tol) { @@ -1463,7 +1570,36 @@ i_t knapsack_generation_t::generate_flow_cover_cut( } } const f_t sgfci_violation = sgfci_lhs_value - snf_b; - if (best_violation <= tol && sgfci_violation <= tol) { return -1; } + if (best_violation <= tol && sgfci_violation <= tol) { + if (flow_cover_diag_verbose) { + i_t n1_count = 0; + i_t n2_count = 0; + i_t c1_count = 0; + i_t c2_count = 0; + for (i_t k = 0; k < static_cast(arcs.size()); k++) { + n1_count += arcs[k].in_n2 ? 0 : 1; + n2_count += arcs[k].in_n2 ? 1 : 0; + c1_count += in_c1[k] ? 1 : 0; + c2_count += in_c2[k] ? 1 : 0; + } + settings.log.printf( + "Flow cover row %d rejected: no violation cmir=%e sgfci=%e lambda=%e ubar_count=%d " + "snf_b=%e cover_capacity=%e arcs=%d n1=%d n2=%d c1=%d c2=%d\n", + flow_cover_row, + best_violation, + sgfci_violation, + lambda, + static_cast(ubar_candidates.size()), + snf_b, + cover_capacity, + static_cast(arcs.size()), + n1_count, + n2_count, + c1_count, + c2_count); + } + return FLOW_COVER_REJECT_NO_VIOLATION; + } const bool use_cmirfci = best_violation >= sgfci_violation; f_t lhs_constant = 0.0; @@ -1550,12 +1686,12 @@ i_t knapsack_generation_t::generate_flow_cover_cut( const f_t coeff = lhs_coefficients[j]; if (std::abs(coeff) > 1e-9) { cut.push_back(j, -coeff); } } - if (cut.size() == 0) { return -1; } + if (cut.size() == 0) { return FLOW_COVER_REJECT_EMPTY_CUT; } cut.rhs = lhs_constant - snf_b; cut.sort(); const f_t dot = cut.vector.dot(xstar); - if (dot >= cut.rhs - tol) { return -1; } + if (dot >= cut.rhs - tol) { return FLOW_COVER_REJECT_MAPPED_CUT; } if (verbose) { settings.log.printf("Flow cover row %d violation %e lambda %e ubar %e type %s\n", @@ -1566,7 +1702,7 @@ i_t knapsack_generation_t::generate_flow_cover_cut( use_cmirfci ? "c-MIRFCI" : "SGFCI"); } - return 0; + return FLOW_COVER_CUT_FOUND; } template @@ -2552,13 +2688,40 @@ void cut_generation_t::generate_flow_cover_cuts( f_t start_time) { if (knapsack_generation_.num_flow_cover_constraints() > 0) { + const bool flow_cover_diag = std::getenv("CUOPT_FLOW_COVER_DIAG") != nullptr; + std::unordered_map status_counts; for (i_t flow_cover_row : knapsack_generation_.get_flow_cover_constraints()) { if (toc(start_time) >= settings.time_limit) { return; } inequality_t cut(lp.num_cols); i_t status = knapsack_generation_.generate_flow_cover_cut( lp, settings, Arow, new_slacks, variable_bounds, var_types, xstar, flow_cover_row, cut); + if (flow_cover_diag) { status_counts[status]++; } if (status == 0) { cut_pool_.add_cut(cut_type_t::FLOW_COVER, cut); } } + if (flow_cover_diag) { + settings.log.printf("Flow cover diagnostics: rows=%d\n", + knapsack_generation_.num_flow_cover_constraints()); + const std::array statuses = {FLOW_COVER_CUT_FOUND, + FLOW_COVER_REJECT_SLACK, + FLOW_COVER_REJECT_RHS, + FLOW_COVER_REJECT_CONTINUOUS_TERM, + FLOW_COVER_REJECT_RELAXATION, + FLOW_COVER_REJECT_NO_ARCS, + FLOW_COVER_REJECT_NO_COVER, + FLOW_COVER_REJECT_NO_ITEMS, + FLOW_COVER_REJECT_KPSNF, + FLOW_COVER_REJECT_LAMBDA, + FLOW_COVER_REJECT_UBAR, + FLOW_COVER_REJECT_NO_VIOLATION, + FLOW_COVER_REJECT_EMPTY_CUT, + FLOW_COVER_REJECT_MAPPED_CUT}; + for (int status : statuses) { + auto count_it = status_counts.find(status); + if (count_it != status_counts.end()) { + settings.log.printf(" %s: %d\n", flow_cover_status_name(status), count_it->second); + } + } + } } } @@ -3509,7 +3672,47 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, if (settings.sub_mip) { return; // Don't compute the variable upper/lower bounds inside sub-MIP } - f_t start_time = tic(); + f_t start_time = tic(); + const f_t bound_tol = 1e-9; + constexpr i_t max_propagation_rounds = 3; + constexpr i_t max_propagated_bounds_per_variable = 128; + + struct variable_bound_edge_t { + i_t variable; + f_t weight; + f_t bias; + }; + + auto is_binary_variable = [&](i_t j) { + return var_types[j] == variable_type_t::INTEGER && std::abs(lp.lower[j]) <= bound_tol && + std::abs(lp.upper[j] - 1.0) <= bound_tol; + }; + + auto valid_bound_edge = [&](const variable_bound_edge_t& edge) { + return edge.variable >= 0 && edge.variable < lp.num_cols && std::isfinite(edge.weight) && + std::isfinite(edge.bias); + }; + + auto already_has_bound = [](const std::vector& bounds, + const variable_bound_edge_t& candidate) { + constexpr f_t duplicate_tol = 1e-9; + for (const auto& bound : bounds) { + if (bound.variable != candidate.variable) { continue; } + const f_t scale = std::max(1.0, + std::max({std::abs(bound.weight), + std::abs(bound.bias), + std::abs(candidate.weight), + std::abs(candidate.bias)})); + if (std::abs(bound.weight - candidate.weight) <= duplicate_tol * scale && + std::abs(bound.bias - candidate.bias) <= duplicate_tol * scale) { + return true; + } + } + return false; + }; + + std::vector> continuous_upper_bounds(lp.num_cols); + std::vector> continuous_lower_bounds(lp.num_cols); // Construct the slack map slack_map_.resize(lp.num_rows, -1); @@ -3598,12 +3801,13 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, const i_t col_end = lp.A.col_start[j + 1]; for (i_t p = col_start; p < col_end; p++) { const i_t i = lp.A.i[p]; - if (num_integer_in_row[i] < 1) { continue; } if (num_neg_inf_[i] > 2 && num_pos_inf_[i] > 2) { continue; } 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; if (row_len < 2) { continue; } + const bool small_continuous_linking_row = num_integer_in_row[i] == 0 && row_len <= 3; + if (num_integer_in_row[i] < 1 && !small_continuous_linking_row) { continue; } const f_t a_ij = lp.A.x[p]; const f_t slack_lower = lp.lower[slack_map_[i]]; const f_t slack_upper = lp.upper[slack_map_[i]]; @@ -3627,7 +3831,7 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, // we can't derive a bound for (i_t q = row_start; q < row_end; q++) { const i_t l = Arow.j[q]; - if (var_types[l] == variable_type_t::CONTINUOUS) { continue; } + if (l == j || l == slack_map_[i]) { continue; } // sum_{k != l, k != j} a_ik x_k + a_ij x_j + a_il x_l <= beta // a_ij x_j <= -a_il x_l + beta - sum_{k != l, k != j} a_ik x_k const f_t a_il = Arow.x[q]; @@ -3638,10 +3842,15 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, // We have a valid variable upper bound // x_j <= -a_il/a_ij * x_l + beta/a_ij - 1/a_ij * sum_{k != l, k != j} a_ik * // bound(x_k) - upper_variables.push_back(l); - upper_weights.push_back(-a_il / a_ij); - upper_biases.push_back(beta / a_ij - (1.0 / a_ij) * sum); - upper_edges++; + const variable_bound_edge_t edge{l, -a_il / a_ij, beta / a_ij - (1.0 / a_ij) * sum}; + if (var_types[l] == variable_type_t::CONTINUOUS) { + continuous_upper_bounds[j].push_back(edge); + } else { + upper_variables.push_back(edge.variable); + upper_weights.push_back(edge.weight); + upper_biases.push_back(edge.bias); + upper_edges++; + } } } } @@ -3658,7 +3867,7 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, for (i_t q = row_start; q < row_end; q++) { const i_t l = Arow.j[q]; - if (var_types[l] == variable_type_t::CONTINUOUS) { continue; } + if (l == j || l == slack_map_[i]) { continue; } // sum_{k != l, k != j} a_ik x_k + a_ij x_j + a_il x_l >= beta // a_ij x_j >= -a_il x_l + beta - sum_{k != l, k != j} a_ik x_k const f_t a_il = Arow.x[q]; @@ -3669,10 +3878,15 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, // We have a valid variable upper bound // x_j <= -a_il/a_ij * x_l + beta/a_ij - 1/a_ij * sum_{k != l, k != j} a_ik * // bound(x_k) - upper_variables.push_back(l); - upper_weights.push_back(-a_il / a_ij); - upper_biases.push_back(beta / a_ij - (1.0 / a_ij) * sum); - upper_edges++; + const variable_bound_edge_t edge{l, -a_il / a_ij, beta / a_ij - (1.0 / a_ij) * sum}; + if (var_types[l] == variable_type_t::CONTINUOUS) { + continuous_upper_bounds[j].push_back(edge); + } else { + upper_variables.push_back(edge.variable); + upper_weights.push_back(edge.weight); + upper_biases.push_back(edge.bias); + upper_edges++; + } } } } @@ -3690,12 +3904,13 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, const i_t col_start = lp.A.col_start[j]; const i_t col_end = lp.A.col_start[j + 1]; for (i_t p = col_start; p < col_end; p++) { - const i_t i = lp.A.i[p]; - if (num_integer_in_row[i] < 1) { continue; } + const i_t i = lp.A.i[p]; 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; if (row_len < 2) { continue; } + const bool small_continuous_linking_row = num_integer_in_row[i] == 0 && row_len <= 3; + if (num_integer_in_row[i] < 1 && !small_continuous_linking_row) { continue; } const f_t a_ij = lp.A.x[p]; const f_t slack_lower = lp.lower[slack_map_[i]]; const f_t slack_upper = lp.upper[slack_map_[i]]; @@ -3714,7 +3929,7 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, for (i_t q = row_start; q < row_end; q++) { const i_t l = Arow.j[q]; - if (var_types[l] == variable_type_t::CONTINUOUS) { continue; } + if (l == j || l == slack_map_[i]) { continue; } // sum_{k != l, k != j} a_ik x_k + a_ij x_j + a_il x_l <= beta // a_ij x_j <= -a_il x_l + beta - sum_{k != l, k != j} a_ik x_k // x_j >= -a_il/a_ij * x_l + beta/a_ij - 1/a_ij * sum_{k != l, k != j} a_ik * bound(x_k) @@ -3726,10 +3941,15 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, // We have a valid variable lower bound // x_j >= -a_il/a_ij * x_l + beta/a_ij - 1/a_ij * sum_{k != l, k != j} a_ik * // bound(x_k) - lower_variables.push_back(l); - lower_weights.push_back(-a_il / a_ij); - lower_biases.push_back(beta / a_ij - (1.0 / a_ij) * sum); - lower_edges++; + const variable_bound_edge_t edge{l, -a_il / a_ij, beta / a_ij - (1.0 / a_ij) * sum}; + if (var_types[l] == variable_type_t::CONTINUOUS) { + continuous_lower_bounds[j].push_back(edge); + } else { + lower_variables.push_back(edge.variable); + lower_weights.push_back(edge.weight); + lower_biases.push_back(edge.bias); + lower_edges++; + } } } } @@ -3746,7 +3966,7 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, for (i_t q = row_start; q < row_end; q++) { const i_t l = Arow.j[q]; - if (var_types[l] == variable_type_t::CONTINUOUS) { continue; } + if (l == j || l == slack_map_[i]) { continue; } // sum_{k != l, k != j} a_ik x_k + a_ij x_j + a_il x_l >= beta // a_ij x_j >= -a_il x_l + beta - sum_{k != l, k != j} a_ik x_k const f_t a_il = Arow.x[q]; @@ -3757,10 +3977,15 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, // We have a valid variable lower bound // x_j >= -a_il/a_ij * x_l + beta/a_ij - 1/a_ij * sum_{k != l, k != j} a_ik * // bound(x_k) - lower_variables.push_back(l); - lower_weights.push_back(-a_il / a_ij); - lower_biases.push_back(beta / a_ij - (1.0 / a_ij) * sum); - lower_edges++; + const variable_bound_edge_t edge{l, -a_il / a_ij, beta / a_ij - (1.0 / a_ij) * sum}; + if (var_types[l] == variable_type_t::CONTINUOUS) { + continuous_lower_bounds[j].push_back(edge); + } else { + lower_variables.push_back(edge.variable); + lower_weights.push_back(edge.weight); + lower_biases.push_back(edge.bias); + lower_edges++; + } } } } @@ -3769,6 +3994,120 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, } lower_offsets[lp.num_cols] = lower_edges; settings.log.printf("%d variable lower bounds in %.2f seconds\n", lower_edges, toc(start_time)); + + std::vector> upper_bounds(lp.num_cols); + std::vector> lower_bounds(lp.num_cols); + std::vector> binary_upper_bounds(lp.num_cols); + std::vector> binary_lower_bounds(lp.num_cols); + + for (i_t j = 0; j < lp.num_cols; j++) { + upper_bounds[j].reserve(upper_offsets[j + 1] - upper_offsets[j]); + for (i_t p = upper_offsets[j]; p < upper_offsets[j + 1]; p++) { + variable_bound_edge_t edge{upper_variables[p], upper_weights[p], upper_biases[p]}; + upper_bounds[j].push_back(edge); + if (is_binary_variable(edge.variable)) { binary_upper_bounds[j].push_back(edge); } + } + lower_bounds[j].reserve(lower_offsets[j + 1] - lower_offsets[j]); + for (i_t p = lower_offsets[j]; p < lower_offsets[j + 1]; p++) { + variable_bound_edge_t edge{lower_variables[p], lower_weights[p], lower_biases[p]}; + lower_bounds[j].push_back(edge); + if (is_binary_variable(edge.variable)) { binary_lower_bounds[j].push_back(edge); } + } + } + + std::vector propagated_upper_count(lp.num_cols, 0); + std::vector propagated_lower_count(lp.num_cols, 0); + i_t propagated_upper_edges = 0; + i_t propagated_lower_edges = 0; + + auto add_propagated_bound = [&](i_t j, + variable_bound_edge_t edge, + std::vector>& all_bounds, + std::vector>& binary_bounds, + std::vector& propagated_count, + i_t& propagated_edges) { + if (!valid_bound_edge(edge) || !is_binary_variable(edge.variable) || edge.variable == j) { + return false; + } + if (std::abs(edge.weight) <= bound_tol) { return false; } + if (propagated_count[j] >= max_propagated_bounds_per_variable) { return false; } + if (already_has_bound(binary_bounds[j], edge)) { return false; } + binary_bounds[j].push_back(edge); + all_bounds[j].push_back(edge); + propagated_count[j]++; + propagated_edges++; + return true; + }; + + auto compose_bound = [](const variable_bound_edge_t& outer, const variable_bound_edge_t& inner) { + return variable_bound_edge_t{ + inner.variable, outer.weight * inner.weight, outer.bias + outer.weight * inner.bias}; + }; + + for (i_t round = 0; round < max_propagation_rounds; round++) { + bool changed = false; + for (i_t j = 0; j < lp.num_cols; j++) { + if (var_types[j] != variable_type_t::CONTINUOUS) { continue; } + + for (const auto& edge : continuous_upper_bounds[j]) { + if (!valid_bound_edge(edge) || std::abs(edge.weight) <= bound_tol) { continue; } + const auto& source_bounds = edge.weight > 0.0 ? binary_upper_bounds[edge.variable] + : binary_lower_bounds[edge.variable]; + for (const auto& source_bound : source_bounds) { + changed |= add_propagated_bound(j, + compose_bound(edge, source_bound), + upper_bounds, + binary_upper_bounds, + propagated_upper_count, + propagated_upper_edges); + } + } + + for (const auto& edge : continuous_lower_bounds[j]) { + if (!valid_bound_edge(edge) || std::abs(edge.weight) <= bound_tol) { continue; } + const auto& source_bounds = edge.weight > 0.0 ? binary_lower_bounds[edge.variable] + : binary_upper_bounds[edge.variable]; + for (const auto& source_bound : source_bounds) { + changed |= add_propagated_bound(j, + compose_bound(edge, source_bound), + lower_bounds, + binary_lower_bounds, + propagated_lower_count, + propagated_lower_edges); + } + } + } + if (!changed) { break; } + } + + auto flatten_bounds = [&](const std::vector>& bounds, + std::vector& offsets, + std::vector& variables, + std::vector& weights, + std::vector& biases) { + offsets.assign(lp.num_cols + 1, 0); + variables.clear(); + weights.clear(); + biases.clear(); + for (i_t j = 0; j < lp.num_cols; j++) { + offsets[j] = static_cast(variables.size()); + for (const auto& edge : bounds[j]) { + variables.push_back(edge.variable); + weights.push_back(edge.weight); + biases.push_back(edge.bias); + } + } + offsets[lp.num_cols] = static_cast(variables.size()); + }; + + if (propagated_upper_edges > 0 || propagated_lower_edges > 0) { + flatten_bounds(upper_bounds, upper_offsets, upper_variables, upper_weights, upper_biases); + flatten_bounds(lower_bounds, lower_offsets, lower_variables, lower_weights, lower_biases); + settings.log.printf( + "%d propagated variable upper bounds and %d propagated variable lower bounds\n", + propagated_upper_edges, + propagated_lower_edges); + } } template diff --git a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp index d94cf5aa67..f25dfb762c 100644 --- a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp +++ b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp @@ -577,7 +577,6 @@ void set_presolve_methods(papilo::Presolve& presolver, presolver.addPresolveMethod(uptr(new papilo::DualInfer())); presolver.addPresolveMethod(uptr(new papilo::SimpleSubstitution())); presolver.addPresolveMethod(uptr(new papilo::Sparsify())); - presolver.addPresolveMethod(uptr(new papilo::Substitution())); } else { CUOPT_LOG_INFO("Disabling the presolver methods that do not support dual postsolve"); } diff --git a/cpp/tests/mip/cuts_test.cu b/cpp/tests/mip/cuts_test.cu index 223f6cea6a..158c595116 100644 --- a/cpp/tests/mip/cuts_test.cu +++ b/cpp/tests/mip/cuts_test.cu @@ -1552,6 +1552,43 @@ mps_parser::mps_data_model_t create_mixed_flow_cover_problem_exampl return problem; } +mps_parser::mps_data_model_t create_variable_bound_chain_problem() +{ + mps_parser::mps_data_model_t problem; + + // x0 is binary, x1 and x2 are continuous, and the rows imply + // x2 <= 2 * x1 <= 2 * x0. + std::vector offsets = {0, 2, 4}; + std::vector indices = {2, 1, 1, 0}; + std::vector coefficients = {1.0, -2.0, 1.0, -1.0}; + problem.set_csr_constraint_matrix(coefficients.data(), + coefficients.size(), + indices.data(), + indices.size(), + offsets.data(), + offsets.size()); + + std::vector lower_bounds(2, -std::numeric_limits::infinity()); + std::vector upper_bounds = {0.0, 0.0}; + problem.set_constraint_lower_bounds(lower_bounds.data(), lower_bounds.size()); + problem.set_constraint_upper_bounds(upper_bounds.data(), upper_bounds.size()); + + std::vector var_lower_bounds = {0.0, 0.0, 0.0}; + std::vector var_upper_bounds = { + 1.0, std::numeric_limits::infinity(), std::numeric_limits::infinity()}; + problem.set_variable_lower_bounds(var_lower_bounds.data(), var_lower_bounds.size()); + problem.set_variable_upper_bounds(var_upper_bounds.data(), var_upper_bounds.size()); + + std::vector objective_coefficients(3, 0.0); + problem.set_objective_coefficients(objective_coefficients.data(), objective_coefficients.size()); + + std::vector variable_types = {'I', 'C', 'C'}; + problem.set_variable_types(variable_types); + + problem.set_maximize(false); + return problem; +} + struct flow_cover_test_problem_t { raft::handle_t handle; dual_simplex::simplex_solver_settings_t settings; @@ -1816,7 +1853,7 @@ TEST(cuts, flow_cover_example_9_8_skips_nonviolated_point) xstar, generator.get_flow_cover_constraints()[0], cut); - EXPECT_EQ(status, -1); + EXPECT_LT(status, 0); } TEST(cuts, flow_cover_mixed_single_node_flow_generates_valid_cut) @@ -1858,4 +1895,28 @@ TEST(cuts, flow_cover_mixed_single_node_flow_generates_valid_cut) EXPECT_GT(generated_cuts, 0); } +TEST(cuts, variable_bounds_propagates_continuous_link_to_binary_bound) +{ + auto test_problem = build_flow_cover_test_problem(create_variable_bound_chain_problem()); + + dual_simplex::variable_bounds_t variable_bounds(test_problem.lp, + test_problem.settings, + test_problem.var_types, + test_problem.Arow, + test_problem.new_slacks); + + bool found_propagated_bound = false; + const int binary_col = 0; + const int continuous_col = 2; + for (int p = variable_bounds.upper_offsets[continuous_col]; + p < variable_bounds.upper_offsets[continuous_col + 1]; + p++) { + if (variable_bounds.upper_variables[p] != binary_col) { continue; } + found_propagated_bound |= std::abs(variable_bounds.upper_weights[p] - 2.0) <= 1e-9 && + std::abs(variable_bounds.upper_biases[p]) <= 1e-9; + } + + EXPECT_TRUE(found_propagated_bound); +} + } // namespace cuopt::linear_programming::test From 79265a2d5a59686d09d8b16b1be12111de5d8f04 Mon Sep 17 00:00:00 2001 From: Hugo Linsenmaier Date: Wed, 29 Apr 2026 23:37:26 +0000 Subject: [PATCH 06/17] Use wolter 5.1 approximation for knapsack --- cpp/src/cuts/cuts.cpp | 52 ++++++++++++++++++++++++++++++++----------- cpp/src/cuts/cuts.hpp | 6 +++++ 2 files changed, 45 insertions(+), 13 deletions(-) diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index ae47921358..89738dd4c4 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -1034,18 +1034,6 @@ i_t knapsack_generation_t::generate_flow_cover_cut( return snf_arc_t{u, in_n2, x_col, x_value, y_const, y_col, y_coeff, y_x_coeff, y_value}; }; - auto solve_strict_kpsnf = [&](const std::vector& values, - const std::vector& weights, - f_t strict_capacity, - std::vector& solution) { - const i_t n = static_cast(weights.size()); - solution.assign(n, 0.0); - if (n == 0) { return static_cast(0.0); } - - const f_t relaxed_capacity = std::max(0.0, strict_capacity - 1e-9); - return greedy_knapsack_problem(values, weights, relaxed_capacity, solution); - }; - const bool reverse_row_requested = flow_cover_row < 0; const i_t actual_flow_cover_row = reverse_row_requested ? -flow_cover_row - 1 : flow_cover_row; inequality_t row(Arow, actual_flow_cover_row, lp.rhs[actual_flow_cover_row]); @@ -1406,7 +1394,7 @@ i_t knapsack_generation_t::generate_flow_cover_cut( if (values.empty()) { return FLOW_COVER_REJECT_NO_ITEMS; } std::vector solution; - const f_t objective = solve_strict_kpsnf(values, weights, cover_capacity, solution); + const f_t objective = prefix_ratio_knapsack_problem(values, weights, cover_capacity, solution); if (std::isnan(objective)) { return FLOW_COVER_REJECT_KPSNF; } static_cast(objective); @@ -2275,6 +2263,44 @@ f_t knapsack_generation_t::greedy_knapsack_problem(const std::vector +f_t knapsack_generation_t::prefix_ratio_knapsack_problem( + const std::vector& values, + const std::vector& weights, + f_t strict_rhs, + std::vector& solution) +{ + const i_t n = static_cast(weights.size()); + solution.assign(n, 0.0); + if (n == 0) { return static_cast(0.0); } + + // Wolter Algorithm 5.1 for KPSNF^rat: sort by p_j / u_j, add items while + // the strict capacity remains satisfied, and stop at the first item that + // does not fit. This is intentionally not the generic greedy knapsack + // heuristic, which keeps scanning for smaller later items. + std::vector permutation(n); + std::iota(permutation.begin(), permutation.end(), 0); + std::sort(permutation.begin(), permutation.end(), [&](i_t a, i_t b) { + const f_t ratio_a = values[a] / weights[a]; + const f_t ratio_b = values[b] / weights[b]; + if (ratio_a != ratio_b) { return ratio_a > ratio_b; } + return weights[a] > weights[b]; + }); + + f_t total_weight = 0.0; + f_t total_value = 0.0; + for (i_t item : permutation) { + if (total_weight + weights[item] < strict_rhs) { + solution[item] = 1.0; + total_weight += weights[item]; + total_value += values[item]; + } else { + break; + } + } + return total_value; +} + template f_t knapsack_generation_t::solve_knapsack_problem(const std::vector& values, const std::vector& weights, diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index 0aa74ec534..4741c4b330 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -394,6 +394,12 @@ class knapsack_generation_t { f_t rhs, std::vector& solution); + // Generate a ratio-sorted prefix solution to a strict 0-1 knapsack problem. + f_t prefix_ratio_knapsack_problem(const std::vector& values, + const std::vector& weights, + f_t strict_rhs, + std::vector& solution); + // Solve a 0-1 knapsack problem using dynamic programming f_t solve_knapsack_problem(const std::vector& values, const std::vector& weights, From f490e33f27cc565129d5a8ab7d337fddacc96bd8 Mon Sep 17 00:00:00 2001 From: Hugo Linsenmaier Date: Thu, 30 Apr 2026 00:24:20 +0000 Subject: [PATCH 07/17] Fix invalid cut generation --- cpp/src/cuts/cuts.cpp | 99 +++++++++++++++++++++++++++---------------- 1 file changed, 62 insertions(+), 37 deletions(-) diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index 89738dd4c4..fb4bed846b 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -941,26 +941,35 @@ knapsack_generation_t::knapsack_generation_t( settings.log.printf("Number of knapsack constraints %d\n", num_knapsack_constraints); #endif - // Wolter's separator is applied to every row after attempting to construct a 0-1 - // single-node-flow relaxation. Rows that do not admit the relaxation are rejected - // inside generate_flow_cover_cut. - flow_cover_constraints_.reserve(lp.num_rows); + // Wolter's separator is applied to every finite row side after attempting to construct a 0-1 + // single-node-flow relaxation. Rows that do not admit the relaxation are rejected inside + // generate_flow_cover_cut. + flow_cover_constraints_.reserve(2 * lp.num_rows); for (i_t i = 0; i < lp.num_rows; i++) { if (Arow.row_start[i + 1] <= Arow.row_start[i]) { continue; } - flow_cover_constraints_.push_back(i); - bool has_slack = false; - i_t slack_col = -1; + i_t slack_col = -1; + f_t slack_coeff = 0.0; + i_t slack_count = 0; for (i_t p = Arow.row_start[i]; p < Arow.row_start[i + 1]; p++) { if (is_slack_[Arow.j[p]]) { - has_slack = true; - slack_col = Arow.j[p]; - break; + slack_col = Arow.j[p]; + slack_coeff = Arow.x[p]; + slack_count++; } } - const bool equality_row = - !has_slack || (slack_col >= 0 && std::abs(lp.lower[slack_col]) <= 1e-9 && - std::abs(lp.upper[slack_col]) <= 1e-9); - if (equality_row) { flow_cover_constraints_.push_back(-(i + 1)); } + if (slack_count > 1 || (slack_count == 1 && std::abs(slack_coeff) <= 1e-9)) { continue; } + if (slack_count == 0) { + flow_cover_constraints_.push_back(i); + flow_cover_constraints_.push_back(-(i + 1)); + continue; + } + + const f_t sigma_slack_lower = + slack_coeff > 0.0 ? slack_coeff * lp.lower[slack_col] : slack_coeff * lp.upper[slack_col]; + const f_t sigma_slack_upper = + slack_coeff > 0.0 ? slack_coeff * lp.upper[slack_col] : slack_coeff * lp.lower[slack_col]; + if (sigma_slack_lower > -inf) { flow_cover_constraints_.push_back(i); } + if (sigma_slack_upper < inf) { flow_cover_constraints_.push_back(-(i + 1)); } } } @@ -1038,26 +1047,40 @@ i_t knapsack_generation_t::generate_flow_cover_cut( const i_t actual_flow_cover_row = reverse_row_requested ? -flow_cover_row - 1 : flow_cover_row; inequality_t row(Arow, actual_flow_cover_row, lp.rhs[actual_flow_cover_row]); i_t slack_count = 0; + i_t slack_col = -1; f_t slack_coeff = 0.0; - bool negate_row = reverse_row_requested; const i_t row_len = row.size(); for (i_t k = 0; k < row_len; k++) { const i_t j = row.index(k); if (!is_slack_[j]) { continue; } slack_count++; + slack_col = j; slack_coeff = row.coeff(k); } if (slack_count > 1) { return FLOW_COVER_REJECT_SLACK; } + bool negate_row = reverse_row_requested; + f_t b = negate_row ? -row.rhs : row.rhs; if (slack_count == 1) { if (std::abs(slack_coeff) <= tol) { return FLOW_COVER_REJECT_SLACK; } - negate_row = reverse_row_requested ? slack_coeff > 0.0 : slack_coeff < 0.0; + const f_t sigma_slack_lower = + slack_coeff > 0.0 ? slack_coeff * lp.lower[slack_col] : slack_coeff * lp.upper[slack_col]; + const f_t sigma_slack_upper = + slack_coeff > 0.0 ? slack_coeff * lp.upper[slack_col] : slack_coeff * lp.lower[slack_col]; + if (!reverse_row_requested) { + if (sigma_slack_lower <= -inf) { return FLOW_COVER_REJECT_RHS; } + negate_row = false; + b = row.rhs - sigma_slack_lower; + } else { + if (sigma_slack_upper >= inf) { return FLOW_COVER_REJECT_RHS; } + negate_row = true; + b = -row.rhs + sigma_slack_upper; + } } std::vector> continuous_terms; continuous_terms.reserve(row_len); std::vector binary_columns; std::unordered_map binary_coefficients; - f_t b = negate_row ? -row.rhs : row.rhs; if (!std::isfinite(b)) { return FLOW_COVER_REJECT_RHS; } for (i_t k = 0; k < row_len; k++) { @@ -2264,11 +2287,10 @@ f_t knapsack_generation_t::greedy_knapsack_problem(const std::vector -f_t knapsack_generation_t::prefix_ratio_knapsack_problem( - const std::vector& values, - const std::vector& weights, - f_t strict_rhs, - std::vector& solution) +f_t knapsack_generation_t::prefix_ratio_knapsack_problem(const std::vector& values, + const std::vector& weights, + f_t strict_rhs, + std::vector& solution) { const i_t n = static_cast(weights.size()); solution.assign(n, 0.0); @@ -2280,11 +2302,10 @@ f_t knapsack_generation_t::prefix_ratio_knapsack_problem( // heuristic, which keeps scanning for smaller later items. std::vector permutation(n); std::iota(permutation.begin(), permutation.end(), 0); - std::sort(permutation.begin(), permutation.end(), [&](i_t a, i_t b) { + std::stable_sort(permutation.begin(), permutation.end(), [&](i_t a, i_t b) { const f_t ratio_a = values[a] / weights[a]; const f_t ratio_b = values[b] / weights[b]; - if (ratio_a != ratio_b) { return ratio_a > ratio_b; } - return weights[a] > weights[b]; + return ratio_a > ratio_b; }); f_t total_weight = 0.0; @@ -3834,12 +3855,14 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, if (row_len < 2) { continue; } const bool small_continuous_linking_row = num_integer_in_row[i] == 0 && row_len <= 3; if (num_integer_in_row[i] < 1 && !small_continuous_linking_row) { continue; } - const f_t a_ij = lp.A.x[p]; - const f_t slack_lower = lp.lower[slack_map_[i]]; - const f_t slack_upper = lp.upper[slack_map_[i]]; - const f_t slack_coeff_i = slack_coeff[i]; - const f_t sigma_slack_lower = slack_coeff_i == 1.0 ? slack_lower : -slack_upper; - const f_t sigma_slack_upper = slack_coeff_i == 1.0 ? slack_upper : -slack_lower; + const f_t a_ij = lp.A.x[p]; + const f_t slack_lower = lp.lower[slack_map_[i]]; + const f_t slack_upper = lp.upper[slack_map_[i]]; + const f_t slack_coeff_i = slack_coeff[i]; + const f_t sigma_slack_lower = + slack_coeff_i > 0.0 ? slack_coeff_i * slack_lower : slack_coeff_i * slack_upper; + const f_t sigma_slack_upper = + slack_coeff_i > 0.0 ? slack_coeff_i * slack_upper : slack_coeff_i * slack_lower; if (sigma_slack_lower > -inf) { const f_t beta = lp.rhs[i] - sigma_slack_lower; @@ -3937,12 +3960,14 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, if (row_len < 2) { continue; } const bool small_continuous_linking_row = num_integer_in_row[i] == 0 && row_len <= 3; if (num_integer_in_row[i] < 1 && !small_continuous_linking_row) { continue; } - const f_t a_ij = lp.A.x[p]; - const f_t slack_lower = lp.lower[slack_map_[i]]; - const f_t slack_upper = lp.upper[slack_map_[i]]; - const f_t slack_coeff_i = slack_coeff[i]; - const f_t sigma_slack_lower = slack_coeff_i == 1.0 ? slack_lower : -slack_upper; - const f_t sigma_slack_upper = slack_coeff_i == 1.0 ? slack_upper : -slack_lower; + const f_t a_ij = lp.A.x[p]; + const f_t slack_lower = lp.lower[slack_map_[i]]; + const f_t slack_upper = lp.upper[slack_map_[i]]; + const f_t slack_coeff_i = slack_coeff[i]; + const f_t sigma_slack_lower = + slack_coeff_i > 0.0 ? slack_coeff_i * slack_lower : slack_coeff_i * slack_upper; + const f_t sigma_slack_upper = + slack_coeff_i > 0.0 ? slack_coeff_i * slack_upper : slack_coeff_i * slack_lower; if (sigma_slack_lower > -inf) { const f_t beta = lp.rhs[i] - sigma_slack_lower; From 6e7481b8bfd614bfbe7b1d0cf55332bd36ea8e2f Mon Sep 17 00:00:00 2001 From: Hugo Linsenmaier Date: Thu, 30 Apr 2026 00:29:31 +0000 Subject: [PATCH 08/17] Add asserts --- cpp/src/cuts/cuts.cpp | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index fb4bed846b..4019fd106f 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -1083,6 +1083,21 @@ i_t knapsack_generation_t::generate_flow_cover_cut( std::unordered_map binary_coefficients; if (!std::isfinite(b)) { return FLOW_COVER_REJECT_RHS; } + cuopt_assert( + [&]() { + f_t row_side_activity = 0.0; + f_t row_side_scale = std::max(1.0, std::abs(b)); + for (i_t k = 0; k < row_len; k++) { + const i_t j = row.index(k); + if (is_slack_[j]) { continue; } + const f_t coeff = negate_row ? -row.coeff(k) : row.coeff(k); + row_side_activity += coeff * xstar[j]; + row_side_scale += std::abs(coeff * xstar[j]); + } + return row_side_activity <= b + 1e-5 * row_side_scale; + }(), + "Flow cover normalized row side excludes LP solution"); + for (i_t k = 0; k < row_len; k++) { const i_t j = row.index(k); if (is_slack_[j]) { continue; } @@ -1379,6 +1394,21 @@ i_t knapsack_generation_t::generate_flow_cover_cut( if (arcs.empty()) { return FLOW_COVER_REJECT_NO_ARCS; } f_t snf_b = b - b_shift; + cuopt_assert( + [&]() { + f_t snf_activity = 0.0; + f_t snf_scale = std::max(1.0, std::abs(snf_b)); + for (const auto& arc : arcs) { + const f_t arc_tol = 1e-5 * std::max(1.0, arc.u); + if (arc.y_value < -arc_tol) { return false; } + if (arc.y_value > arc.u * arc.x_value + arc_tol) { return false; } + const f_t signed_y = arc.in_n2 ? -arc.y_value : arc.y_value; + snf_activity += signed_y; + snf_scale += std::abs(signed_y); + } + return snf_activity <= snf_b + 1e-5 * snf_scale; + }(), + "Flow cover SNF relaxation excludes LP solution"); auto compute_structure = [&]() { bool has_binary_controlled_arc = false; From c44b30c35386eb8d18556991f1a63bce9631448a Mon Sep 17 00:00:00 2001 From: Hugo Linsenmaier Date: Thu, 30 Apr 2026 02:32:44 +0000 Subject: [PATCH 09/17] Clean up --- cpp/src/cuts/cuts.cpp | 213 ++++-------------------------------------- 1 file changed, 18 insertions(+), 195 deletions(-) diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index 4019fd106f..004b1d3f01 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -15,7 +15,6 @@ #include #include #include -#include #include #include #include @@ -35,44 +34,6 @@ namespace { enum class clique_cut_build_status_t : int8_t { NO_CUT = 0, CUT_ADDED = 1, INFEASIBLE = 2 }; -enum flow_cover_cut_status_t : int { - FLOW_COVER_CUT_FOUND = 0, - FLOW_COVER_REJECT_SLACK = -1, - FLOW_COVER_REJECT_RHS = -2, - FLOW_COVER_REJECT_CONTINUOUS_TERM = -3, - FLOW_COVER_REJECT_RELAXATION = -4, - FLOW_COVER_REJECT_NO_ARCS = -5, - FLOW_COVER_REJECT_NO_COVER = -6, - FLOW_COVER_REJECT_NO_ITEMS = -7, - FLOW_COVER_REJECT_KPSNF = -8, - FLOW_COVER_REJECT_LAMBDA = -9, - FLOW_COVER_REJECT_UBAR = -10, - FLOW_COVER_REJECT_NO_VIOLATION = -11, - FLOW_COVER_REJECT_EMPTY_CUT = -12, - FLOW_COVER_REJECT_MAPPED_CUT = -13, -}; - -const char* flow_cover_status_name(int status) -{ - switch (status) { - case FLOW_COVER_CUT_FOUND: return "cut_found"; - case FLOW_COVER_REJECT_SLACK: return "bad_slack"; - case FLOW_COVER_REJECT_RHS: return "bad_rhs"; - case FLOW_COVER_REJECT_CONTINUOUS_TERM: return "continuous_term_no_candidate"; - case FLOW_COVER_REJECT_RELAXATION: return "invalid_snf_relaxation"; - case FLOW_COVER_REJECT_NO_ARCS: return "no_arcs"; - case FLOW_COVER_REJECT_NO_COVER: return "no_cover_capacity"; - case FLOW_COVER_REJECT_NO_ITEMS: return "no_kpsnf_items"; - case FLOW_COVER_REJECT_KPSNF: return "kpsnf_failed"; - case FLOW_COVER_REJECT_LAMBDA: return "no_positive_lambda"; - case FLOW_COVER_REJECT_UBAR: return "no_ubar_candidate"; - case FLOW_COVER_REJECT_NO_VIOLATION: return "no_cmir_or_sgfci_violation"; - case FLOW_COVER_REJECT_EMPTY_CUT: return "empty_mapped_cut"; - case FLOW_COVER_REJECT_MAPPED_CUT: return "mapped_cut_not_violated"; - default: return "unknown"; - } -} - #if DEBUG_CLIQUE_CUTS #define CLIQUE_CUTS_DEBUG(...) \ do { \ @@ -985,11 +946,7 @@ i_t knapsack_generation_t::generate_flow_cover_cut( i_t flow_cover_row, inequality_t& cut) { - const bool verbose = false; - const bool flow_cover_diag_verbose = std::getenv("CUOPT_FLOW_COVER_DIAG_VERBOSE") != nullptr; - const char* flow_cover_diag_row_env = std::getenv("CUOPT_FLOW_COVER_DIAG_ROW"); - const i_t flow_cover_diag_row = - flow_cover_diag_row_env != nullptr ? static_cast(std::atoi(flow_cover_diag_row_env)) : -1; + static_cast(settings); static_cast(new_slacks); const f_t tol = 1e-6; const f_t bound_tol = 1e-8; @@ -1057,21 +1014,21 @@ i_t knapsack_generation_t::generate_flow_cover_cut( slack_col = j; slack_coeff = row.coeff(k); } - if (slack_count > 1) { return FLOW_COVER_REJECT_SLACK; } + if (slack_count > 1) { return -1; } bool negate_row = reverse_row_requested; f_t b = negate_row ? -row.rhs : row.rhs; if (slack_count == 1) { - if (std::abs(slack_coeff) <= tol) { return FLOW_COVER_REJECT_SLACK; } + if (std::abs(slack_coeff) <= tol) { return -1; } const f_t sigma_slack_lower = slack_coeff > 0.0 ? slack_coeff * lp.lower[slack_col] : slack_coeff * lp.upper[slack_col]; const f_t sigma_slack_upper = slack_coeff > 0.0 ? slack_coeff * lp.upper[slack_col] : slack_coeff * lp.lower[slack_col]; if (!reverse_row_requested) { - if (sigma_slack_lower <= -inf) { return FLOW_COVER_REJECT_RHS; } + if (sigma_slack_lower <= -inf) { return -1; } negate_row = false; b = row.rhs - sigma_slack_lower; } else { - if (sigma_slack_upper >= inf) { return FLOW_COVER_REJECT_RHS; } + if (sigma_slack_upper >= inf) { return -1; } negate_row = true; b = -row.rhs + sigma_slack_upper; } @@ -1081,7 +1038,7 @@ i_t knapsack_generation_t::generate_flow_cover_cut( continuous_terms.reserve(row_len); std::vector binary_columns; std::unordered_map binary_coefficients; - if (!std::isfinite(b)) { return FLOW_COVER_REJECT_RHS; } + if (!std::isfinite(b)) { return -1; } cuopt_assert( [&]() { @@ -1321,63 +1278,14 @@ i_t knapsack_generation_t::generate_flow_cover_cut( add_variable_lower_candidates(); add_simple_bound_candidates(); - if (candidates.empty()) { - if (flow_cover_diag_verbose) { - const i_t upper_start = variable_bounds.upper_offsets[j]; - const i_t upper_end = variable_bounds.upper_offsets[j + 1]; - const i_t lower_start = variable_bounds.lower_offsets[j]; - const i_t lower_end = variable_bounds.lower_offsets[j + 1]; - settings.log.printf( - "Flow cover row %d rejected: continuous col=%d coeff=%e lower=%e upper=%e " - "type=%d upper_vb=%d lower_vb=%d\n", - flow_cover_row, - j, - c, - lower_j, - upper_j, - static_cast(var_types[j]), - upper_end - upper_start, - lower_end - lower_start); - if (flow_cover_diag_row == flow_cover_row) { - for (i_t p = upper_start; p < upper_end; p++) { - const i_t x_col = variable_bounds.upper_variables[p]; - const f_t direct_coeff = - binary_coefficients.count(x_col) != 0 ? binary_coefficients[x_col] : 0.0; - settings.log.printf( - " upper candidate bound_var=%d gamma=%e alpha=%e direct_coeff=%e " - "binary=%d xstar=%e\n", - x_col, - variable_bounds.upper_weights[p], - variable_bounds.upper_biases[p], - direct_coeff, - static_cast(is_binary_variable(x_col)), - x_col >= 0 ? xstar[x_col] : 0.0); - } - for (i_t p = lower_start; p < lower_end; p++) { - const i_t x_col = variable_bounds.lower_variables[p]; - const f_t direct_coeff = - binary_coefficients.count(x_col) != 0 ? binary_coefficients[x_col] : 0.0; - settings.log.printf( - " lower candidate bound_var=%d gamma=%e alpha=%e direct_coeff=%e " - "binary=%d xstar=%e\n", - x_col, - variable_bounds.lower_weights[p], - variable_bounds.lower_biases[p], - direct_coeff, - static_cast(is_binary_variable(x_col)), - x_col >= 0 ? xstar[x_col] : 0.0); - } - } - } - return FLOW_COVER_REJECT_CONTINUOUS_TERM; - } + if (candidates.empty()) { return -1; } auto best = std::min_element( candidates.begin(), candidates.end(), [](const snf_candidate_t& a, const snf_candidate_t& b) { return a.distance < b.distance; }); if (best->arc.y_value < -1e-5 || best->arc.y_value - best->arc.u * best->arc.x_value > 1e-4 * std::max(1.0, best->arc.u)) { - return FLOW_COVER_REJECT_RELAXATION; + return -1; } if (best->arc.y_value < 0.0) { best->arc.y_value = 0.0; } arcs.push_back(best->arc); @@ -1392,7 +1300,7 @@ i_t knapsack_generation_t::generate_flow_cover_cut( arcs.push_back(build_arc(u, coeff < 0.0, j, 0.0, 0.0, -1, 0.0, u)); } - if (arcs.empty()) { return FLOW_COVER_REJECT_NO_ARCS; } + if (arcs.empty()) { return -1; } f_t snf_b = b - b_shift; cuopt_assert( [&]() { @@ -1427,7 +1335,7 @@ i_t knapsack_generation_t::generate_flow_cover_cut( auto [has_binary_controlled_arc, has_n1, cover_capacity] = compute_structure(); if (!has_binary_controlled_arc || !has_n1 || cover_capacity <= tol) { - return FLOW_COVER_REJECT_NO_COVER; + return -1; } std::vector values; @@ -1444,11 +1352,11 @@ i_t knapsack_generation_t::generate_flow_cover_cut( weights.push_back(arc.u); item_to_arc.push_back(k); } - if (values.empty()) { return FLOW_COVER_REJECT_NO_ITEMS; } + if (values.empty()) { return -1; } std::vector solution; const f_t objective = prefix_ratio_knapsack_problem(values, weights, cover_capacity, solution); - if (std::isnan(objective)) { return FLOW_COVER_REJECT_KPSNF; } + if (std::isnan(objective)) { return -1; } static_cast(objective); std::vector in_c1(arcs.size(), 0); @@ -1473,18 +1381,7 @@ i_t knapsack_generation_t::generate_flow_cover_cut( if (in_c2[k]) { sum_c2_u += arcs[k].u; } } const f_t lambda = -snf_b + sum_c1_u - sum_c2_u; - if (!c1_nonempty || lambda <= tol) { - if (flow_cover_diag_verbose) { - settings.log.printf( - "Flow cover row %d rejected: lambda=%e snf_b=%e cover_capacity=%e arcs=%d\n", - flow_cover_row, - lambda, - snf_b, - cover_capacity, - static_cast(arcs.size())); - } - return FLOW_COVER_REJECT_LAMBDA; - } + if (!c1_nonempty || lambda <= tol) { return -1; } auto fractional_part_local = [](f_t value) { const f_t floor_value = std::floor(value); @@ -1528,16 +1425,7 @@ i_t knapsack_generation_t::generate_flow_cover_cut( add_ubar_candidate(max_c1); add_ubar_candidate(max_u_ge_lambda + 1.0); add_ubar_candidate(lambda + 1.0); - if (ubar_candidates.empty()) { - if (flow_cover_diag_verbose) { - settings.log.printf("Flow cover row %d rejected: no ubar lambda=%e snf_b=%e arcs=%d\n", - flow_cover_row, - lambda, - snf_b, - static_cast(arcs.size())); - } - return FLOW_COVER_REJECT_UBAR; - } + if (ubar_candidates.empty()) { return -1; } f_t best_violation = 0.0; f_t best_ubar = 0.0; @@ -1611,36 +1499,7 @@ i_t knapsack_generation_t::generate_flow_cover_cut( } } const f_t sgfci_violation = sgfci_lhs_value - snf_b; - if (best_violation <= tol && sgfci_violation <= tol) { - if (flow_cover_diag_verbose) { - i_t n1_count = 0; - i_t n2_count = 0; - i_t c1_count = 0; - i_t c2_count = 0; - for (i_t k = 0; k < static_cast(arcs.size()); k++) { - n1_count += arcs[k].in_n2 ? 0 : 1; - n2_count += arcs[k].in_n2 ? 1 : 0; - c1_count += in_c1[k] ? 1 : 0; - c2_count += in_c2[k] ? 1 : 0; - } - settings.log.printf( - "Flow cover row %d rejected: no violation cmir=%e sgfci=%e lambda=%e ubar_count=%d " - "snf_b=%e cover_capacity=%e arcs=%d n1=%d n2=%d c1=%d c2=%d\n", - flow_cover_row, - best_violation, - sgfci_violation, - lambda, - static_cast(ubar_candidates.size()), - snf_b, - cover_capacity, - static_cast(arcs.size()), - n1_count, - n2_count, - c1_count, - c2_count); - } - return FLOW_COVER_REJECT_NO_VIOLATION; - } + if (best_violation <= tol && sgfci_violation <= tol) { return -1; } const bool use_cmirfci = best_violation >= sgfci_violation; f_t lhs_constant = 0.0; @@ -1727,23 +1586,14 @@ i_t knapsack_generation_t::generate_flow_cover_cut( const f_t coeff = lhs_coefficients[j]; if (std::abs(coeff) > 1e-9) { cut.push_back(j, -coeff); } } - if (cut.size() == 0) { return FLOW_COVER_REJECT_EMPTY_CUT; } + if (cut.size() == 0) { return -1; } cut.rhs = lhs_constant - snf_b; cut.sort(); const f_t dot = cut.vector.dot(xstar); - if (dot >= cut.rhs - tol) { return FLOW_COVER_REJECT_MAPPED_CUT; } + if (dot >= cut.rhs - tol) { return -1; } - if (verbose) { - settings.log.printf("Flow cover row %d violation %e lambda %e ubar %e type %s\n", - flow_cover_row, - use_cmirfci ? best_violation : sgfci_violation, - lambda, - best_ubar, - use_cmirfci ? "c-MIRFCI" : "SGFCI"); - } - - return FLOW_COVER_CUT_FOUND; + return 0; } template @@ -2765,40 +2615,13 @@ void cut_generation_t::generate_flow_cover_cuts( f_t start_time) { if (knapsack_generation_.num_flow_cover_constraints() > 0) { - const bool flow_cover_diag = std::getenv("CUOPT_FLOW_COVER_DIAG") != nullptr; - std::unordered_map status_counts; for (i_t flow_cover_row : knapsack_generation_.get_flow_cover_constraints()) { if (toc(start_time) >= settings.time_limit) { return; } inequality_t cut(lp.num_cols); i_t status = knapsack_generation_.generate_flow_cover_cut( lp, settings, Arow, new_slacks, variable_bounds, var_types, xstar, flow_cover_row, cut); - if (flow_cover_diag) { status_counts[status]++; } if (status == 0) { cut_pool_.add_cut(cut_type_t::FLOW_COVER, cut); } } - if (flow_cover_diag) { - settings.log.printf("Flow cover diagnostics: rows=%d\n", - knapsack_generation_.num_flow_cover_constraints()); - const std::array statuses = {FLOW_COVER_CUT_FOUND, - FLOW_COVER_REJECT_SLACK, - FLOW_COVER_REJECT_RHS, - FLOW_COVER_REJECT_CONTINUOUS_TERM, - FLOW_COVER_REJECT_RELAXATION, - FLOW_COVER_REJECT_NO_ARCS, - FLOW_COVER_REJECT_NO_COVER, - FLOW_COVER_REJECT_NO_ITEMS, - FLOW_COVER_REJECT_KPSNF, - FLOW_COVER_REJECT_LAMBDA, - FLOW_COVER_REJECT_UBAR, - FLOW_COVER_REJECT_NO_VIOLATION, - FLOW_COVER_REJECT_EMPTY_CUT, - FLOW_COVER_REJECT_MAPPED_CUT}; - for (int status : statuses) { - auto count_it = status_counts.find(status); - if (count_it != status_counts.end()) { - settings.log.printf(" %s: %d\n", flow_cover_status_name(status), count_it->second); - } - } - } } } From ff56b94792f8a24383dcf48f996d90874507422c Mon Sep 17 00:00:00 2001 From: Hugo Linsenmaier Date: Tue, 5 May 2026 04:36:27 +0000 Subject: [PATCH 10/17] Re add substitution --- cpp/src/mip_heuristics/presolve/third_party_presolve.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp index f25dfb762c..d94cf5aa67 100644 --- a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp +++ b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp @@ -577,6 +577,7 @@ void set_presolve_methods(papilo::Presolve& presolver, presolver.addPresolveMethod(uptr(new papilo::DualInfer())); presolver.addPresolveMethod(uptr(new papilo::SimpleSubstitution())); presolver.addPresolveMethod(uptr(new papilo::Sparsify())); + presolver.addPresolveMethod(uptr(new papilo::Substitution())); } else { CUOPT_LOG_INFO("Disabling the presolver methods that do not support dual postsolve"); } From 070874a108e90d88a862bdebf8ebd30f24d975f0 Mon Sep 17 00:00:00 2001 From: Hugo Linsenmaier Date: Tue, 5 May 2026 04:58:38 +0000 Subject: [PATCH 11/17] Add comments --- cpp/src/cuts/cuts.cpp | 185 +++++++++++++++++++++++++++++++++++++++--- 1 file changed, 173 insertions(+), 12 deletions(-) diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index 004b1d3f01..16b6b652d6 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -951,23 +951,52 @@ i_t knapsack_generation_t::generate_flow_cover_cut( const f_t tol = 1e-6; const f_t bound_tol = 1e-8; + // Implementation outline: + // + // 1. Pick one finite side of one MIP row and normalize it as a'x <= b. + // 2. Build an explicit 0-1 single-node-flow (SNF) relaxation for that row side: + // + // sum_{j in N1} y_j - sum_{j in N2} y_j <= b_snf, + // 0 <= y_j <= u_j x_j. + // + // Wolter's separator starts from this SNF form. cuOpt starts from general + // MIP rows, so this step is a mapping layer: every synthetic SNF y_j must + // remember how to substitute back to the original variables. + // 3. Run the KPSNF cover-selection heuristic to choose C1 and C2. + // 4. Test c-MIRFCI candidates and the simpler SGFCI fallback at the current LP point. + // 5. Substitute each synthetic y_j back into original variables and emit a cuOpt cut. + // + // A row side is converted to a 0-1 single-node-flow relaxation + // + // sum_{j in N1} y_j - sum_{j in N2} y_j <= b, + // 0 <= y_j <= u_j x_j, x_j in {0, 1}. + // + // Each arc stores both the SNF data (u, N1/N2 membership, binary controller x_j) + // and the affine expression used to map y_j back to the original row variables. struct snf_arc_t { - f_t u; - bool in_n2; - i_t x_col; - f_t x_value; + f_t u; // Capacity u_j in 0 <= y_j <= u_j x_j. + bool in_n2; // false: y_j appears with + sign (N1), true: with - sign (N2). + i_t x_col; // Binary controller column, or -1 for a fixed-control simple-bound arc. + f_t x_value; // Current LP value of x_col, or the fixed-control value. + + // cuOpt mapping layer: synthetic y_j equals + // + // y_const + y_coeff * original_y_col + y_x_coeff * x_col. + // + // This is not thesis notation; it is how we later substitute the SNF cut back + // into the original MIP variables. f_t y_const; i_t y_col; f_t y_coeff; f_t y_x_coeff; - f_t y_value; + f_t y_value; // Current LP value of the synthetic y_j. }; struct snf_candidate_t { - snf_arc_t arc; - f_t b_shift; - f_t distance; - bool absorbs_binary_coeff; + snf_arc_t arc; // One possible SNF rewrite of an original row term. + f_t b_shift; // Constant moved from the row activity into the SNF RHS. + f_t distance; // How close the active bound used by this rewrite is to z*. + bool absorbs_binary_coeff; // True if a direct x coefficient was folded into this arc. }; auto is_binary_variable = [&](i_t j) { @@ -981,6 +1010,9 @@ i_t knapsack_generation_t::generate_flow_cover_cut( return x_col >= 0 ? clamp01(xstar[x_col]) : fixed_value; }; + // Materialize the synthetic y_j mapping and evaluate it at the current LP point. + // Pure simple-bound arcs have no binary controller, so x_col = -1 and fixed_x is + // used only to evaluate constants like u * 1. auto build_arc = [&](f_t u, bool in_n2, i_t x_col, @@ -1015,6 +1047,21 @@ i_t knapsack_generation_t::generate_flow_cover_cut( slack_coeff = row.coeff(k); } if (slack_count > 1) { return -1; } + // A row can contribute either finite side. Internally the row is equality-like: + // + // a x + sigma * slack = rhs. + // + // A finite lower slack activity gives: + // + // a x <= rhs - sigma_slack_lower. + // + // A finite upper slack activity gives: + // + // a x >= rhs - sigma_slack_upper, + // + // which we negate to keep the normalized form a'x <= b. + // + // b is the RHS after moving the selected finite slack activity into the constant. bool negate_row = reverse_row_requested; f_t b = negate_row ? -row.rhs : row.rhs; if (slack_count == 1) { @@ -1034,6 +1081,11 @@ i_t knapsack_generation_t::generate_flow_cover_cut( } } + // Split the normalized row a'x <= b into two groups: + // + // - Binary coefficients a_j x_j can directly become SNF arcs with y_j = |a_j| x_j. + // - Continuous terms c z need a bound z <= gamma x + alpha or z >= gamma x + alpha + // before they can become arcs with 0 <= y <= u x. std::vector> continuous_terms; continuous_terms.reserve(row_len); std::vector binary_columns; @@ -1075,6 +1127,15 @@ i_t knapsack_generation_t::generate_flow_cover_cut( arcs.reserve(row_len); f_t b_shift = 0.0; + // Candidate arcs are alternative algebraic rewrites of one original row term c z. + // For example, with z <= gamma x + alpha and a positive row coefficient c: + // + // c z = c alpha + y, y = c (z - alpha), 0 <= y <= c gamma x. + // + // The c alpha constant is accumulated in b_shift, so the final SNF RHS is + // snf_b = b - b_shift. If the same binary controller x also has a direct row + // coefficient a x, we try two candidates: absorb a x into this arc capacity or + // leave a x for a later pure-binary arc. auto push_candidate = [&](std::vector& candidates, f_t u, bool in_n2, @@ -1105,6 +1166,13 @@ i_t knapsack_generation_t::generate_flow_cover_cut( const f_t y_star = xstar[j]; std::vector candidates; + // Use upper variable bounds z <= gamma x + alpha. With c > 0 this naturally + // creates an N1 arc: + // + // c z = c alpha + y, y = c(z - alpha) + a x, u = c gamma + a. + // + // With c < 0 the signs flip and the same rewrite creates an N2 arc. The guards + // below check that y is nonnegative at the valid endpoint and that u is nonnegative. auto add_variable_upper_candidates = [&]() { const i_t start = variable_bounds.upper_offsets[j]; const i_t end = variable_bounds.upper_offsets[j + 1]; @@ -1116,6 +1184,9 @@ i_t knapsack_generation_t::generate_flow_cover_cut( if (!std::isfinite(gamma) || !std::isfinite(alpha) || lower_j <= -inf) { continue; } const f_t direct_coeff = binary_coefficients.count(x_col) != 0 ? binary_coefficients[x_col] : 0.0; + // First try folding the direct x coefficient into this arc, then try the + // pure variable-bound arc. If there is no direct coefficient, only the pure + // candidate is generated. const std::array a_values = {direct_coeff, 0.0}; const i_t num_a_values = std::abs(direct_coeff) > tol ? 2 : 1; for (i_t h = 0; h < num_a_values; h++) { @@ -1159,6 +1230,9 @@ i_t knapsack_generation_t::generate_flow_cover_cut( } }; + // Use lower variable bounds z >= gamma x + alpha. This is the symmetric case: + // with c > 0 we represent the term through an N2 arc, and with c < 0 through + // an N1 arc. The endpoint checks mirror the upper-bound case. auto add_variable_lower_candidates = [&]() { const i_t start = variable_bounds.lower_offsets[j]; const i_t end = variable_bounds.lower_offsets[j + 1]; @@ -1170,6 +1244,8 @@ i_t knapsack_generation_t::generate_flow_cover_cut( if (!std::isfinite(gamma) || !std::isfinite(alpha) || upper_j >= inf) { continue; } const f_t direct_coeff = binary_coefficients.count(x_col) != 0 ? binary_coefficients[x_col] : 0.0; + // As above, test both variants: absorb a direct x coefficient into the + // arc or keep it as a separate pure-binary arc. const std::array a_values = {direct_coeff, 0.0}; const i_t num_a_values = std::abs(direct_coeff) > tol ? 2 : 1; for (i_t h = 0; h < num_a_values; h++) { @@ -1213,6 +1289,9 @@ i_t knapsack_generation_t::generate_flow_cover_cut( } }; + // If no binary-controlled variable bound is available, finite simple bounds + // l <= z <= u still give a valid fixed-control arc. These arcs can preserve + // the SNF row algebra, but they do not create useful binary cover structure. auto add_simple_bound_candidates = [&]() { if (lower_j <= -inf || upper_j >= inf) { return; } const f_t capacity = std::abs(c) * (upper_j - lower_j); @@ -1279,6 +1358,10 @@ i_t knapsack_generation_t::generate_flow_cover_cut( add_simple_bound_candidates(); if (candidates.empty()) { return -1; } + // Pick the rewrite whose active bound is closest to z*. Several rewrites can + // be algebraically valid for the integer set, but the separator evaluates cuts + // at the current LP point. We therefore keep the rewrite that is least likely + // to exclude (x*, y*) from the constructed SNF relaxation. auto best = std::min_element( candidates.begin(), candidates.end(), [](const snf_candidate_t& a, const snf_candidate_t& b) { return a.distance < b.distance; @@ -1293,6 +1376,10 @@ i_t knapsack_generation_t::generate_flow_cover_cut( if (best->absorbs_binary_coeff) { binary_coefficients[best->arc.x_col] = 0.0; } } + // Remaining pure binary coefficients a x become arcs with y = u x and u = |a|: + // + // a > 0: +u x contributes as an N1 arc, + // a < 0: -u x contributes as an N2 arc. for (i_t j : binary_columns) { const f_t coeff = binary_coefficients[j]; if (std::abs(coeff) <= tol) { continue; } @@ -1318,6 +1405,13 @@ i_t knapsack_generation_t::generate_flow_cover_cut( }(), "Flow cover SNF relaxation excludes LP solution"); + // Coarse structural test before separation. A flow cover exists only if the N1 + // capacities can exceed the RHS: + // + // lambda(C = N1, empty C2) = -b + sum_{j in N1} u_j > 0. + // + // This is the capacity surplus the KPSNF separator tries to keep positive while + // choosing a cover that is violated by the current LP point. auto compute_structure = [&]() { bool has_binary_controlled_arc = false; bool has_n1 = false; @@ -1334,9 +1428,7 @@ i_t knapsack_generation_t::generate_flow_cover_cut( }; auto [has_binary_controlled_arc, has_n1, cover_capacity] = compute_structure(); - if (!has_binary_controlled_arc || !has_n1 || cover_capacity <= tol) { - return -1; - } + if (!has_binary_controlled_arc || !has_n1 || cover_capacity <= tol) { return -1; } std::vector values; std::vector weights; @@ -1344,6 +1436,18 @@ i_t knapsack_generation_t::generate_flow_cover_cut( values.reserve(arcs.size()); weights.reserve(arcs.size()); item_to_arc.reserve(arcs.size()); + // KPSNF separation chooses which arcs to keep out of the cover by solving: + // + // max sum_{N1} (1 - x*_j) zbar_j + sum_{N2} x*_j z_j + // s.t. sum u_j z_j < -b + sum_{N1} u_j. + // + // Interpretation after solving: + // + // N1: zbar_j = 0 means j enters C1. + // N2: z_j = 1 means j enters C2. + // + // The current implementation uses the ratio-prefix approximation from the + // rational KPSNF path instead of the exact integer DP path. for (i_t k = 0; k < static_cast(arcs.size()); k++) { const auto& arc = arcs[k]; if (arc.u <= tol) { continue; } @@ -1359,6 +1463,9 @@ i_t knapsack_generation_t::generate_flow_cover_cut( if (std::isnan(objective)) { return -1; } static_cast(objective); + // Convert the KPSNF solution into the flow-cover sets: + // + // C1 = {j in N1 : zbar_j = 0}, C2 = {j in N2 : z_j = 1}. std::vector in_c1(arcs.size(), 0); std::vector in_c2(arcs.size(), 0); for (i_t item = 0; item < static_cast(solution.size()); item++) { @@ -1380,6 +1487,11 @@ i_t knapsack_generation_t::generate_flow_cover_cut( } if (in_c2[k]) { sum_c2_u += arcs[k].u; } } + // lambda is the cover excess: + // + // lambda = -b + sum_{j in C1} u_j - sum_{j in C2} u_j. + // + // A positive lambda is required for the flow-cover inequality to cut anything. const f_t lambda = -snf_b + sum_c1_u - sum_c2_u; if (!c1_nonempty || lambda <= tol) { return -1; } @@ -1388,12 +1500,18 @@ i_t knapsack_generation_t::generate_flow_cover_cut( return value - floor_value; }; + // MIR rounding function used by the c-MIR flow-cover inequality: + // + // F_alpha(t) = floor(t) + max(0, frac(t) - alpha) / (1 - alpha). auto mir_F = [&](f_t alpha, f_t value) { const f_t floor_value = std::floor(value); const f_t frac_value = value - floor_value; return floor_value + std::max(0.0, frac_value - alpha) / (1.0 - alpha); }; + // ubar controls the c-MIR scaling beta = -lambda / ubar. The rounded coefficients + // only change at capacity-related breakpoints, so we test a small candidate set + // derived from capacities in and around C1, C2, L1, and L2. std::vector ubar_candidates; auto add_ubar_candidate = [&](f_t candidate) { if (candidate <= lambda + tol || !std::isfinite(candidate)) { return; } @@ -1432,6 +1550,9 @@ i_t knapsack_generation_t::generate_flow_cover_cut( std::vector best_in_l1(arcs.size(), 0); std::vector best_in_l2(arcs.size(), 0); + // Evaluate one candidate c-MIRFCI at the current LP point. This mirrors the + // inequality terms but does not write the cut yet; it is only used to choose the + // best ubar and the L1/L2 sets. auto contribution = [&](const snf_arc_t& arc, bool arc_in_c1, bool arc_in_c2, @@ -1454,6 +1575,12 @@ i_t knapsack_generation_t::generate_flow_cover_cut( const f_t f_beta = fractional_part_local(beta); if (f_beta < 0.01 || f_beta > 0.95) { continue; } + // L1 and L2 are chosen by local comparison at (x*, y*): + // + // - For N1 \ C1, use the strengthened expression only if it contributes at + // least as much violation as dropping the term. + // - For N2 \ C2, use the x-term only if it contributes at least as much + // violation as the baseline -y term. std::vector in_l1(arcs.size(), 0); std::vector in_l2(arcs.size(), 0); for (i_t k = 0; k < static_cast(arcs.size()); k++) { @@ -1481,6 +1608,9 @@ i_t knapsack_generation_t::generate_flow_cover_cut( } } + // SGFCI fallback: this is the non-c-MIR strengthened flow-cover inequality. We + // evaluate it too because it can be violated when no tested ubar gives a useful + // c-MIRFCI, and we keep whichever of the two cuts is more violated. std::vector sgfci_in_l2(arcs.size(), 0); f_t sgfci_lhs_value = 0.0; for (i_t k = 0; k < static_cast(arcs.size()); k++) { @@ -1516,6 +1646,16 @@ i_t knapsack_generation_t::generate_flow_cover_cut( } }; + // Map the selected SNF inequality + // + // lhs_snf(x, y) <= b + // + // back through each arc's affine y-expression, accumulating + // + // lhs_constant + sum_j lhs_coefficients[j] * original_x_j <= snf_b. + // + // add_y substitutes synthetic y_j. add_x and add_one_minus_x handle the binary + // controller terms that appear in the c-MIRFCI/SGFCI formula. auto add_y = [&](const snf_arc_t& arc, f_t multiplier) { lhs_constant += multiplier * arc.y_const; if (arc.y_col >= 0) { add_lhs_coeff(arc.y_col, multiplier * arc.y_coeff); } @@ -1580,6 +1720,13 @@ i_t knapsack_generation_t::generate_flow_cover_cut( } } + // cuOpt stores cuts in the opposite >= form. Therefore + // + // lhs_constant + lhs_coefficients * x <= snf_b + // + // becomes + // + // -lhs_coefficients * x >= lhs_constant - snf_b. cut.clear(); cut.reserve(lhs_indices.size()); for (i_t j : lhs_indices) { @@ -3611,6 +3758,11 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, return false; }; + // This is not the flow-cover inequality itself. It is row-to-SNF preprocessing: + // the separator needs bounds like z <= gamma x + alpha or z >= gamma x + alpha + // to turn continuous row terms into arcs with 0 <= y <= u x. Edges whose + // controller is already binary can be used directly; edges through continuous + // variables are kept so short linking chains can be composed into binary bounds. std::vector> continuous_upper_bounds(lp.num_cols); std::vector> continuous_lower_bounds(lp.num_cols); @@ -3943,6 +4095,15 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, return true; }; + // Compose affine bounds during row-to-SNF preprocessing. For example, + // + // z <= a w + b, w <= c x + d + // + // gives + // + // z <= a c x + (b + a d). + // + // If a is negative, the caller swaps upper/lower source bounds. auto compose_bound = [](const variable_bound_edge_t& outer, const variable_bound_edge_t& inner) { return variable_bound_edge_t{ inner.variable, outer.weight * inner.weight, outer.bias + outer.weight * inner.bias}; From b240a996a8d767a03d5c1c89d61cc394925d346d Mon Sep 17 00:00:00 2001 From: Hugo Linsenmaier Date: Tue, 5 May 2026 05:21:59 +0000 Subject: [PATCH 12/17] Add example --- cpp/src/cuts/cuts.cpp | 55 ++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 54 insertions(+), 1 deletion(-) diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index 16b6b652d6..86fd1e705a 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -954,6 +954,14 @@ i_t knapsack_generation_t::generate_flow_cover_cut( // Implementation outline: // // 1. Pick one finite side of one MIP row and normalize it as a'x <= b. + // + // Example: + // + // 4 z + 3 x2 - 2 x3 <= 11, + // 1 <= z <= 5 x1 + 1, + // + // where x1, x2, and x3 are binary and z is continuous. + // // 2. Build an explicit 0-1 single-node-flow (SNF) relaxation for that row side: // // sum_{j in N1} y_j - sum_{j in N2} y_j <= b_snf, @@ -962,10 +970,33 @@ i_t knapsack_generation_t::generate_flow_cover_cut( // Wolter's separator starts from this SNF form. cuOpt starts from general // MIP rows, so this step is a mapping layer: every synthetic SNF y_j must // remember how to substitute back to the original variables. + // + // In the example: + // + // 4 z -> y1 = 4(z - 1), 0 <= y1 <= 20 x1, y1 in N1 + // 3 x2 -> y2 = 3 x2, 0 <= y2 <= 3 x2, y2 in N1 + // -2 x3 -> y3 = 2 x3, 0 <= y3 <= 2 x3, y3 in N2 + // + // The constant 4 from 4 z = 4 + y1 shifts the RHS, giving + // + // y1 + y2 - y3 <= 7. + // // 3. Run the KPSNF cover-selection heuristic to choose C1 and C2. + // + // For example, if KPSNF chooses C1 = {y1, y2} and C2 = {y3}, then + // + // lambda = -7 + 20 + 3 - 2 = 14. + // // 4. Test c-MIRFCI candidates and the simpler SGFCI fallback at the current LP point. + // + // Continuing the example, ubar candidates are capacities such as 20 or + // lambda + 1; L1/L2 are chosen by comparing candidate terms at (x*, y*). + // // 5. Substitute each synthetic y_j back into original variables and emit a cuOpt cut. // + // In the example, any final SNF inequality is mapped back with + // y1 = 4z - 4, y2 = 3x2, and y3 = 2x3. + // // A row side is converted to a 0-1 single-node-flow relaxation // // sum_{j in N1} y_j - sum_{j in N2} y_j <= b, @@ -984,7 +1015,8 @@ i_t knapsack_generation_t::generate_flow_cover_cut( // y_const + y_coeff * original_y_col + y_x_coeff * x_col. // // This is not thesis notation; it is how we later substitute the SNF cut back - // into the original MIP variables. + // into the original MIP variables. In the running example, the flow arc stores + // y1 = -4 + 4z, while the direct binary arcs store y2 = 3x2 and y3 = 2x3. f_t y_const; i_t y_col; f_t y_coeff; @@ -1086,6 +1118,9 @@ i_t knapsack_generation_t::generate_flow_cover_cut( // - Binary coefficients a_j x_j can directly become SNF arcs with y_j = |a_j| x_j. // - Continuous terms c z need a bound z <= gamma x + alpha or z >= gamma x + alpha // before they can become arcs with 0 <= y <= u x. + // + // In the running example, z is the only continuous term; x2 and x3 stay in + // binary_coefficients until they are converted to direct arcs. std::vector> continuous_terms; continuous_terms.reserve(row_len); std::vector binary_columns; @@ -1136,6 +1171,9 @@ i_t knapsack_generation_t::generate_flow_cover_cut( // snf_b = b - b_shift. If the same binary controller x also has a direct row // coefficient a x, we try two candidates: absorb a x into this arc capacity or // leave a x for a later pure-binary arc. + // + // In the running example, c = 4, gamma = 5, alpha = 1, and a = 0, so the + // candidate is y1 = 4z - 4 with capacity u = 20 and b_shift = 4. auto push_candidate = [&](std::vector& candidates, f_t u, bool in_n2, @@ -1173,6 +1211,7 @@ i_t knapsack_generation_t::generate_flow_cover_cut( // // With c < 0 the signs flip and the same rewrite creates an N2 arc. The guards // below check that y is nonnegative at the valid endpoint and that u is nonnegative. + // The running example uses this branch for 4z with z <= 5x1 + 1. auto add_variable_upper_candidates = [&]() { const i_t start = variable_bounds.upper_offsets[j]; const i_t end = variable_bounds.upper_offsets[j + 1]; @@ -1380,6 +1419,9 @@ i_t knapsack_generation_t::generate_flow_cover_cut( // // a > 0: +u x contributes as an N1 arc, // a < 0: -u x contributes as an N2 arc. + // + // In the running example, 3x2 becomes y2 = 3x2 in N1 and -2x3 becomes + // -y3 with y3 = 2x3 in N2. for (i_t j : binary_columns) { const f_t coeff = binary_coefficients[j]; if (std::abs(coeff) <= tol) { continue; } @@ -1412,6 +1454,9 @@ i_t knapsack_generation_t::generate_flow_cover_cut( // // This is the capacity surplus the KPSNF separator tries to keep positive while // choosing a cover that is violated by the current LP point. + // + // In the running example, N1 capacities are 20 and 3, the SNF RHS is 7, and + // the coarse surplus is -7 + 20 + 3 = 16 before any N2 arc is put in C2. auto compute_structure = [&]() { bool has_binary_controlled_arc = false; bool has_n1 = false; @@ -1448,6 +1493,9 @@ i_t knapsack_generation_t::generate_flow_cover_cut( // // The current implementation uses the ratio-prefix approximation from the // rational KPSNF path instead of the exact integer DP path. + // + // In the running example, the KPSNF items have weights 20, 3, and 2. The + // values are 1 - x1*, 1 - x2*, and x3* because y1/y2 are in N1 and y3 is in N2. for (i_t k = 0; k < static_cast(arcs.size()); k++) { const auto& arc = arcs[k]; if (arc.u <= tol) { continue; } @@ -1492,6 +1540,8 @@ i_t knapsack_generation_t::generate_flow_cover_cut( // lambda = -b + sum_{j in C1} u_j - sum_{j in C2} u_j. // // A positive lambda is required for the flow-cover inequality to cut anything. + // In the running example with C1 = {y1, y2} and C2 = {y3}, lambda is + // -7 + 20 + 3 - 2 = 14. const f_t lambda = -snf_b + sum_c1_u - sum_c2_u; if (!c1_nonempty || lambda <= tol) { return -1; } @@ -1512,6 +1562,7 @@ i_t knapsack_generation_t::generate_flow_cover_cut( // ubar controls the c-MIR scaling beta = -lambda / ubar. The rounded coefficients // only change at capacity-related breakpoints, so we test a small candidate set // derived from capacities in and around C1, C2, L1, and L2. + // In the running example, one natural ubar candidate is 20, the largest C1 capacity. std::vector ubar_candidates; auto add_ubar_candidate = [&](f_t candidate) { if (candidate <= lambda + tol || !std::isfinite(candidate)) { return; } @@ -1656,6 +1707,8 @@ i_t knapsack_generation_t::generate_flow_cover_cut( // // add_y substitutes synthetic y_j. add_x and add_one_minus_x handle the binary // controller terms that appear in the c-MIRFCI/SGFCI formula. + // In the running example, add_y(y1, m) contributes -4m to the constant and + // 4m to the coefficient of z; add_y(y2, m) contributes 3m to x2. auto add_y = [&](const snf_arc_t& arc, f_t multiplier) { lhs_constant += multiplier * arc.y_const; if (arc.y_col >= 0) { add_lhs_coeff(arc.y_col, multiplier * arc.y_coeff); } From e62eecd2415750402df51e4918995e9314cd0e13 Mon Sep 17 00:00:00 2001 From: Hugo Linsenmaier Date: Tue, 5 May 2026 05:34:17 +0000 Subject: [PATCH 13/17] Generate minimalist tests --- cpp/tests/mip/cuts_test.cu | 319 ++++++------------------------------- 1 file changed, 52 insertions(+), 267 deletions(-) diff --git a/cpp/tests/mip/cuts_test.cu b/cpp/tests/mip/cuts_test.cu index 4ca252faa3..8e1271db60 100644 --- a/cpp/tests/mip/cuts_test.cu +++ b/cpp/tests/mip/cuts_test.cu @@ -1446,33 +1446,22 @@ TEST(cuts, clique_neos8_phase4_lp_infeasibility_binary_search) EXPECT_EQ(first_infeasible.value(), injected_index); } -// Problem data for Example 9.8 from Wolsey: flow cover cut test. +// Minimal 0-1 single-node-flow relaxation for the flow-cover separator. // -// The mixed 0-1 set X has the knapsack constraint: -// 3*x0 + 3*x1 + 6*x2 - 3*x3 - 5*x4 - 1*x5 <= 4 -// with all binary variables (0-indexed: x0..x5 correspond to 1-indexed x1..x6). +// y0 + y1 - y2 <= 4 +// 0 <= y0 <= 3*x0, 0 <= y1 <= 6*x1, 0 <= y2 <= 3*x2 // -// The LP relaxation solution is x* = (1, 0, 2/3, 1, 0, 0). -// Wolter's c-MIRFCI separator can generate, for example: -// 2*x0 + 2*x2 - 2*x4 - x5 <= 2 -// which is violated by x*: 2 + 4/3 > 2. -// The always-approximate KPSNF flow-cover selection may choose a different -// violated valid cut, so the tests below check validity rather than exact coefficients. -// -// We construct an objective that drives the LP relaxation toward this fractional solution. -mps_parser::mps_data_model_t create_flow_cover_problem_example_9_8() +// The fractional point x* = (1, 2/3, 1), y* = (3, 4, 3) satisfies the relaxation +// but violates the generated c-MIR flow-cover cut. This is a reduced version of a +// standard flow-cover example; the test checks validity instead of exact coefficients +// because the approximate KPSNF selection may choose a different valid cut. +mps_parser::mps_data_model_t create_small_single_node_flow_problem() { mps_parser::mps_data_model_t problem; - // 6 binary variables x0..x5 - // Constraint: 3*x0 + 3*x1 + 6*x2 - 3*x3 - 5*x4 - 1*x5 <= 4 - // We also add x3 <= 1 explicitly (it is binary, but coefficients are negative in the knapsack - // row so the solver needs it to identify flow cover structure). - // - // CSR format (1 row): - std::vector offsets = {0, 6}; - std::vector indices = {0, 1, 2, 3, 4, 5}; - std::vector coefficients = {3.0, 3.0, 6.0, -3.0, -5.0, -1.0}; + std::vector offsets = {0, 3, 5, 7, 9}; + std::vector indices = {3, 4, 5, 0, 3, 1, 4, 2, 5}; + std::vector coefficients = {1.0, 1.0, -1.0, -3.0, 1.0, -6.0, 1.0, -3.0, 1.0}; problem.set_csr_constraint_matrix(coefficients.data(), coefficients.size(), indices.data(), @@ -1480,84 +1469,20 @@ mps_parser::mps_data_model_t create_flow_cover_problem_example_9_8( offsets.data(), offsets.size()); - // Constraint bounds: -inf <= sum <= 4 - std::vector lower_bounds = {-std::numeric_limits::infinity()}; - std::vector upper_bounds = {4.0}; + std::vector lower_bounds(4, -std::numeric_limits::infinity()); + std::vector upper_bounds = {4.0, 0.0, 0.0, 0.0}; problem.set_constraint_lower_bounds(lower_bounds.data(), lower_bounds.size()); problem.set_constraint_upper_bounds(upper_bounds.data(), upper_bounds.size()); - // Variable bounds: all binary [0, 1] - std::vector var_lower_bounds = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; - std::vector var_upper_bounds = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0}; + std::vector var_lower_bounds(6, 0.0); + std::vector var_upper_bounds = {1.0, 1.0, 1.0, 3.0, 6.0, 3.0}; problem.set_variable_lower_bounds(var_lower_bounds.data(), var_lower_bounds.size()); problem.set_variable_upper_bounds(var_upper_bounds.data(), var_upper_bounds.size()); - // Objective: maximize 3*x0 + 6*x2 + 3*x3 (minimize negated) - // This pushes x0=1, x2 as high as possible, x3=1 (opening the RHS). - // LP relaxation gives x0=1, x2=2/3, x3=1 => obj = 3 + 4 + 3 = 10 - // Integer optimum would be x0=1, x2=0, x3=1 or similar. - std::vector objective_coefficients = {-3.0, 0.0, -6.0, -3.0, 0.0, 0.0}; + std::vector objective_coefficients(6, 0.0); problem.set_objective_coefficients(objective_coefficients.data(), objective_coefficients.size()); - // All integer (binary) - std::vector variable_types = {'I', 'I', 'I', 'I', 'I', 'I'}; - problem.set_variable_types(variable_types); - - problem.set_maximize(false); - return problem; -} - -mps_parser::mps_data_model_t create_mixed_flow_cover_problem_example_9_8() -{ - mps_parser::mps_data_model_t problem; - - // Variables 0..5 are binary x variables. Variables 6..11 are continuous y variables. - // The rows model: - // y0 + y1 + y2 <= 4 + y3 + y4 + y5 - // yj <= u_j * x_j - // for u = (3, 3, 6, 3, 5, 1). - std::vector offsets = {0, 6, 8, 10, 12, 14, 16, 18}; - std::vector indices = {6, 7, 8, 9, 10, 11, 0, 6, 1, 7, 2, 8, 3, 9, 4, 10, 5, 11}; - std::vector coefficients = {1.0, - 1.0, - 1.0, - -1.0, - -1.0, - -1.0, - -3.0, - 1.0, - -3.0, - 1.0, - -6.0, - 1.0, - -3.0, - 1.0, - -5.0, - 1.0, - -1.0, - 1.0}; - problem.set_csr_constraint_matrix(coefficients.data(), - coefficients.size(), - indices.data(), - indices.size(), - offsets.data(), - offsets.size()); - - std::vector lower_bounds(7, -std::numeric_limits::infinity()); - std::vector upper_bounds = {4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; - problem.set_constraint_lower_bounds(lower_bounds.data(), lower_bounds.size()); - problem.set_constraint_upper_bounds(upper_bounds.data(), upper_bounds.size()); - - std::vector var_lower_bounds(12, 0.0); - std::vector var_upper_bounds = { - 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 3.0, 3.0, 6.0, 3.0, 5.0, 1.0}; - problem.set_variable_lower_bounds(var_lower_bounds.data(), var_lower_bounds.size()); - problem.set_variable_upper_bounds(var_upper_bounds.data(), var_upper_bounds.size()); - - std::vector objective_coefficients(12, 0.0); - problem.set_objective_coefficients(objective_coefficients.data(), objective_coefficients.size()); - - std::vector variable_types = {'I', 'I', 'I', 'I', 'I', 'I', 'C', 'C', 'C', 'C', 'C', 'C'}; + std::vector variable_types = {'I', 'I', 'I', 'C', 'C', 'C'}; problem.set_variable_types(variable_types); problem.set_maximize(false); @@ -1633,100 +1558,89 @@ flow_cover_test_problem_t build_flow_cover_test_problem( return test_problem; } -flow_cover_test_problem_t build_flow_cover_test_problem() -{ - return build_flow_cover_test_problem(create_flow_cover_problem_example_9_8()); -} - -bool is_flow_cover_example_9_8_binary_feasible(const std::vector& x) -{ - const double activity = 3.0 * x[0] + 3.0 * x[1] + 6.0 * x[2] - 3.0 * x[3] - 5.0 * x[4] - x[5]; - return activity <= 4.0 + 1e-9; -} - -std::vector mixed_flow_cover_fractional_solution(int num_cols) +std::vector single_node_flow_fractional_solution(int num_cols) { std::vector xstar(num_cols, 0.0); xstar[0] = 1.0; - xstar[2] = 2.0 / 3.0; - xstar[3] = 1.0; - xstar[6] = 3.0; - xstar[8] = 4.0; - xstar[9] = 3.0; + xstar[1] = 2.0 / 3.0; + xstar[2] = 1.0; + xstar[3] = 3.0; + xstar[4] = 4.0; + xstar[5] = 3.0; return xstar; } -bool mixed_flow_cover_y_feasible(const std::vector& y) +bool single_node_flow_y_feasible(const std::vector& y) { - const double activity = y[0] + y[1] + y[2] - y[3] - y[4] - y[5]; + const double activity = y[0] + y[1] - y[2]; return activity <= 4.0 + 1e-8; } -void expect_mixed_flow_cover_cut_valid_at_point(const dual_simplex::inequality_t& cut, +void expect_single_node_flow_cut_valid_at_point(const dual_simplex::inequality_t& cut, const std::vector& point, const std::string& label) { EXPECT_GE(cut.vector.dot(point), cut.rhs - 1e-7) << label; } -void expect_mixed_flow_cover_cut_valid_at_simple_extreme_points( +void expect_single_node_flow_cut_valid_at_extreme_points( const dual_simplex::inequality_t& cut, int num_cols) { - const std::vector capacities = {3.0, 3.0, 6.0, 3.0, 5.0, 1.0}; - const std::vector flow_signs = {1.0, 1.0, 1.0, -1.0, -1.0, -1.0}; + const std::vector capacities = {3.0, 6.0, 3.0}; + const std::vector flow_signs = {1.0, 1.0, -1.0}; int checked_points = 0; - for (int x_mask = 0; x_mask < 64; x_mask++) { - std::vector y_upper(6, 0.0); - for (int j = 0; j < 6; j++) { + for (int x_mask = 0; x_mask < 8; x_mask++) { + std::vector y_upper(3, 0.0); + for (int j = 0; j < 3; j++) { if (((x_mask >> j) & 1) != 0) { y_upper[j] = capacities[j]; } } - for (int y_mask = 0; y_mask < 64; y_mask++) { - std::vector y(6, 0.0); - for (int j = 0; j < 6; j++) { + for (int y_mask = 0; y_mask < 8; y_mask++) { + std::vector y(3, 0.0); + for (int j = 0; j < 3; j++) { if (((y_mask >> j) & 1) != 0) { y[j] = y_upper[j]; } } - if (!mixed_flow_cover_y_feasible(y)) { continue; } + if (!single_node_flow_y_feasible(y)) { continue; } std::vector point(num_cols, 0.0); - for (int j = 0; j < 6; j++) { + for (int j = 0; j < 3; j++) { point[j] = ((x_mask >> j) & 1) != 0 ? 1.0 : 0.0; - point[6 + j] = y[j]; + point[3 + j] = y[j]; } - expect_mixed_flow_cover_cut_valid_at_point( + expect_single_node_flow_cut_valid_at_point( cut, point, "box vertex x_mask=" + std::to_string(x_mask) + " y_mask=" + std::to_string(y_mask)); checked_points++; } - for (int free_j = 0; free_j < 6; free_j++) { - for (int bound_mask = 0; bound_mask < 32; bound_mask++) { - std::vector y(6, 0.0); + for (int free_j = 0; free_j < 3; free_j++) { + for (int bound_mask = 0; bound_mask < 4; bound_mask++) { + std::vector y(3, 0.0); int bit = 0; - for (int j = 0; j < 6; j++) { + for (int j = 0; j < 3; j++) { if (j == free_j) { continue; } if (((bound_mask >> bit) & 1) != 0) { y[j] = y_upper[j]; } bit++; } double fixed_activity = 0.0; - for (int j = 0; j < 6; j++) { + for (int j = 0; j < 3; j++) { if (j != free_j) { fixed_activity += flow_signs[j] * y[j]; } } const double y_free = (4.0 - fixed_activity) / flow_signs[free_j]; if (y_free < -1e-8 || y_free > y_upper[free_j] + 1e-8) { continue; } y[free_j] = std::max(0.0, std::min(y_upper[free_j], y_free)); - if (!mixed_flow_cover_y_feasible(y)) { continue; } + if (!single_node_flow_y_feasible(y)) { continue; } std::vector point(num_cols, 0.0); - for (int j = 0; j < 6; j++) { + for (int j = 0; j < 3; j++) { point[j] = ((x_mask >> j) & 1) != 0 ? 1.0 : 0.0; - point[6 + j] = y[j]; + point[3 + j] = y[j]; } - expect_mixed_flow_cover_cut_valid_at_point( + expect_single_node_flow_cut_valid_at_point( cut, point, "flow-tight vertex x_mask=" + std::to_string(x_mask) + @@ -1739,139 +1653,10 @@ void expect_mixed_flow_cover_cut_valid_at_simple_extreme_points( EXPECT_GT(checked_points, 0); } -TEST(cuts, flow_cover_example_9_8_generates_expected_cut) -{ - auto test_problem = build_flow_cover_test_problem(); - std::vector xstar(test_problem.lp.num_cols, 0.0); - xstar[0] = 1.0; - xstar[2] = 2.0 / 3.0; - xstar[3] = 1.0; - - dual_simplex::knapsack_generation_t generator(test_problem.lp, - test_problem.settings, - test_problem.Arow, - test_problem.new_slacks, - test_problem.var_types); - dual_simplex::variable_bounds_t variable_bounds(test_problem.lp, - test_problem.settings, - test_problem.var_types, - test_problem.Arow, - test_problem.new_slacks); - ASSERT_EQ(generator.num_flow_cover_constraints(), 1); - - dual_simplex::inequality_t cut(test_problem.lp.num_cols); - const int status = generator.generate_flow_cover_cut(test_problem.lp, - test_problem.settings, - test_problem.Arow, - test_problem.new_slacks, - variable_bounds, - test_problem.var_types, - xstar, - generator.get_flow_cover_constraints()[0], - cut); - ASSERT_EQ(status, 0); - - const double dot = cut.vector.dot(xstar); - EXPECT_LT(dot, cut.rhs - 1e-6); - - int feasible_assignments = 0; - for (int mask = 0; mask < 64; mask++) { - std::vector point(test_problem.lp.num_cols, 0.0); - for (int j = 0; j < 6; j++) { - point[j] = ((mask >> j) & 1) != 0 ? 1.0 : 0.0; - } - if (!is_flow_cover_example_9_8_binary_feasible(point)) { continue; } - EXPECT_GE(cut.vector.dot(point), cut.rhs - 1e-7) << "mask=" << mask; - feasible_assignments++; - } - EXPECT_GT(feasible_assignments, 0); -} - -TEST(cuts, flow_cover_example_9_8_generated_cut_valid_for_binary_feasible_points) -{ - auto test_problem = build_flow_cover_test_problem(); - std::vector xstar(test_problem.lp.num_cols, 0.0); - xstar[0] = 1.0; - xstar[2] = 2.0 / 3.0; - xstar[3] = 1.0; - - dual_simplex::knapsack_generation_t generator(test_problem.lp, - test_problem.settings, - test_problem.Arow, - test_problem.new_slacks, - test_problem.var_types); - dual_simplex::variable_bounds_t variable_bounds(test_problem.lp, - test_problem.settings, - test_problem.var_types, - test_problem.Arow, - test_problem.new_slacks); - ASSERT_EQ(generator.num_flow_cover_constraints(), 1); - - dual_simplex::inequality_t cut(test_problem.lp.num_cols); - const int status = generator.generate_flow_cover_cut(test_problem.lp, - test_problem.settings, - test_problem.Arow, - test_problem.new_slacks, - variable_bounds, - test_problem.var_types, - xstar, - generator.get_flow_cover_constraints()[0], - cut); - ASSERT_EQ(status, 0); - EXPECT_LT(cut.vector.dot(xstar), cut.rhs - 1e-6); - - int feasible_assignments = 0; - for (int mask = 0; mask < 64; mask++) { - std::vector point(test_problem.lp.num_cols, 0.0); - for (int j = 0; j < 6; j++) { - point[j] = ((mask >> j) & 1) != 0 ? 1.0 : 0.0; - } - if (!is_flow_cover_example_9_8_binary_feasible(point)) { continue; } - EXPECT_GE(cut.vector.dot(point), cut.rhs - 1e-7) << "mask=" << mask; - feasible_assignments++; - } - - EXPECT_GT(feasible_assignments, 0); -} - -TEST(cuts, flow_cover_example_9_8_skips_nonviolated_point) -{ - auto test_problem = build_flow_cover_test_problem(); - std::vector xstar(test_problem.lp.num_cols, 0.0); - xstar[0] = 1.0; - xstar[2] = 1.0; - xstar[3] = 1.0; - xstar[4] = 1.0; - - dual_simplex::knapsack_generation_t generator(test_problem.lp, - test_problem.settings, - test_problem.Arow, - test_problem.new_slacks, - test_problem.var_types); - dual_simplex::variable_bounds_t variable_bounds(test_problem.lp, - test_problem.settings, - test_problem.var_types, - test_problem.Arow, - test_problem.new_slacks); - ASSERT_EQ(generator.num_flow_cover_constraints(), 1); - - dual_simplex::inequality_t cut(test_problem.lp.num_cols); - const int status = generator.generate_flow_cover_cut(test_problem.lp, - test_problem.settings, - test_problem.Arow, - test_problem.new_slacks, - variable_bounds, - test_problem.var_types, - xstar, - generator.get_flow_cover_constraints()[0], - cut); - EXPECT_LT(status, 0); -} - -TEST(cuts, flow_cover_mixed_single_node_flow_generates_valid_cut) +TEST(cuts, flow_cover_generates_valid_single_node_flow_cut) { - auto test_problem = build_flow_cover_test_problem(create_mixed_flow_cover_problem_example_9_8()); - const std::vector xstar = mixed_flow_cover_fractional_solution(test_problem.lp.num_cols); + auto test_problem = build_flow_cover_test_problem(create_small_single_node_flow_problem()); + const std::vector xstar = single_node_flow_fractional_solution(test_problem.lp.num_cols); dual_simplex::knapsack_generation_t generator(test_problem.lp, test_problem.settings, @@ -1900,7 +1685,7 @@ TEST(cuts, flow_cover_mixed_single_node_flow_generates_valid_cut) if (status != 0) { continue; } EXPECT_LT(cut.vector.dot(xstar), cut.rhs - 1e-6) << "row=" << flow_cover_row; - expect_mixed_flow_cover_cut_valid_at_simple_extreme_points(cut, test_problem.lp.num_cols); + expect_single_node_flow_cut_valid_at_extreme_points(cut, test_problem.lp.num_cols); generated_cuts++; } From 45bc40549b4994a51dccbbb1270a73eeb4089717 Mon Sep 17 00:00:00 2001 From: Hugo Linsenmaier Date: Tue, 5 May 2026 06:27:04 +0000 Subject: [PATCH 14/17] Disable variable bound prop --- cpp/src/cuts/cuts.cpp | 187 ++----------------------------------- cpp/tests/mip/cuts_test.cu | 61 ------------ 2 files changed, 10 insertions(+), 238 deletions(-) diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index 86fd1e705a..5a4463b50c 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -3772,10 +3772,7 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, if (settings.sub_mip) { return; // Don't compute the variable upper/lower bounds inside sub-MIP } - f_t start_time = tic(); - const f_t bound_tol = 1e-9; - constexpr i_t max_propagation_rounds = 3; - constexpr i_t max_propagated_bounds_per_variable = 128; + f_t start_time = tic(); struct variable_bound_edge_t { i_t variable; @@ -3783,41 +3780,10 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, f_t bias; }; - auto is_binary_variable = [&](i_t j) { - return var_types[j] == variable_type_t::INTEGER && std::abs(lp.lower[j]) <= bound_tol && - std::abs(lp.upper[j] - 1.0) <= bound_tol; - }; - - auto valid_bound_edge = [&](const variable_bound_edge_t& edge) { - return edge.variable >= 0 && edge.variable < lp.num_cols && std::isfinite(edge.weight) && - std::isfinite(edge.bias); - }; - - auto already_has_bound = [](const std::vector& bounds, - const variable_bound_edge_t& candidate) { - constexpr f_t duplicate_tol = 1e-9; - for (const auto& bound : bounds) { - if (bound.variable != candidate.variable) { continue; } - const f_t scale = std::max(1.0, - std::max({std::abs(bound.weight), - std::abs(bound.bias), - std::abs(candidate.weight), - std::abs(candidate.bias)})); - if (std::abs(bound.weight - candidate.weight) <= duplicate_tol * scale && - std::abs(bound.bias - candidate.bias) <= duplicate_tol * scale) { - return true; - } - } - return false; - }; - // This is not the flow-cover inequality itself. It is row-to-SNF preprocessing: - // the separator needs bounds like z <= gamma x + alpha or z >= gamma x + alpha - // to turn continuous row terms into arcs with 0 <= y <= u x. Edges whose - // controller is already binary can be used directly; edges through continuous - // variables are kept so short linking chains can be composed into binary bounds. - std::vector> continuous_upper_bounds(lp.num_cols); - std::vector> continuous_lower_bounds(lp.num_cols); + // the separator needs direct bounds like z <= gamma x + alpha or + // z >= gamma x + alpha to turn continuous row terms into arcs with 0 <= y <= u x. + // Bounds through continuous linking variables are intentionally not propagated here. // Construct the slack map slack_map_.resize(lp.num_rows, -1); @@ -3911,8 +3877,7 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, const i_t row_end = Arow.row_start[i + 1]; const i_t row_len = row_end - row_start; if (row_len < 2) { continue; } - const bool small_continuous_linking_row = num_integer_in_row[i] == 0 && row_len <= 3; - if (num_integer_in_row[i] < 1 && !small_continuous_linking_row) { continue; } + if (num_integer_in_row[i] < 1) { continue; } const f_t a_ij = lp.A.x[p]; const f_t slack_lower = lp.lower[slack_map_[i]]; const f_t slack_upper = lp.upper[slack_map_[i]]; @@ -3950,9 +3915,7 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, // x_j <= -a_il/a_ij * x_l + beta/a_ij - 1/a_ij * sum_{k != l, k != j} a_ik * // bound(x_k) const variable_bound_edge_t edge{l, -a_il / a_ij, beta / a_ij - (1.0 / a_ij) * sum}; - if (var_types[l] == variable_type_t::CONTINUOUS) { - continuous_upper_bounds[j].push_back(edge); - } else { + if (var_types[l] != variable_type_t::CONTINUOUS) { upper_variables.push_back(edge.variable); upper_weights.push_back(edge.weight); upper_biases.push_back(edge.bias); @@ -3986,9 +3949,7 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, // x_j <= -a_il/a_ij * x_l + beta/a_ij - 1/a_ij * sum_{k != l, k != j} a_ik * // bound(x_k) const variable_bound_edge_t edge{l, -a_il / a_ij, beta / a_ij - (1.0 / a_ij) * sum}; - if (var_types[l] == variable_type_t::CONTINUOUS) { - continuous_upper_bounds[j].push_back(edge); - } else { + if (var_types[l] != variable_type_t::CONTINUOUS) { upper_variables.push_back(edge.variable); upper_weights.push_back(edge.weight); upper_biases.push_back(edge.bias); @@ -4016,8 +3977,7 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, const i_t row_end = Arow.row_start[i + 1]; const i_t row_len = row_end - row_start; if (row_len < 2) { continue; } - const bool small_continuous_linking_row = num_integer_in_row[i] == 0 && row_len <= 3; - if (num_integer_in_row[i] < 1 && !small_continuous_linking_row) { continue; } + if (num_integer_in_row[i] < 1) { continue; } const f_t a_ij = lp.A.x[p]; const f_t slack_lower = lp.lower[slack_map_[i]]; const f_t slack_upper = lp.upper[slack_map_[i]]; @@ -4051,9 +4011,7 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, // x_j >= -a_il/a_ij * x_l + beta/a_ij - 1/a_ij * sum_{k != l, k != j} a_ik * // bound(x_k) const variable_bound_edge_t edge{l, -a_il / a_ij, beta / a_ij - (1.0 / a_ij) * sum}; - if (var_types[l] == variable_type_t::CONTINUOUS) { - continuous_lower_bounds[j].push_back(edge); - } else { + if (var_types[l] != variable_type_t::CONTINUOUS) { lower_variables.push_back(edge.variable); lower_weights.push_back(edge.weight); lower_biases.push_back(edge.bias); @@ -4087,9 +4045,7 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, // x_j >= -a_il/a_ij * x_l + beta/a_ij - 1/a_ij * sum_{k != l, k != j} a_ik * // bound(x_k) const variable_bound_edge_t edge{l, -a_il / a_ij, beta / a_ij - (1.0 / a_ij) * sum}; - if (var_types[l] == variable_type_t::CONTINUOUS) { - continuous_lower_bounds[j].push_back(edge); - } else { + if (var_types[l] != variable_type_t::CONTINUOUS) { lower_variables.push_back(edge.variable); lower_weights.push_back(edge.weight); lower_biases.push_back(edge.bias); @@ -4103,129 +4059,6 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, } lower_offsets[lp.num_cols] = lower_edges; settings.log.printf("%d variable lower bounds in %.2f seconds\n", lower_edges, toc(start_time)); - - std::vector> upper_bounds(lp.num_cols); - std::vector> lower_bounds(lp.num_cols); - std::vector> binary_upper_bounds(lp.num_cols); - std::vector> binary_lower_bounds(lp.num_cols); - - for (i_t j = 0; j < lp.num_cols; j++) { - upper_bounds[j].reserve(upper_offsets[j + 1] - upper_offsets[j]); - for (i_t p = upper_offsets[j]; p < upper_offsets[j + 1]; p++) { - variable_bound_edge_t edge{upper_variables[p], upper_weights[p], upper_biases[p]}; - upper_bounds[j].push_back(edge); - if (is_binary_variable(edge.variable)) { binary_upper_bounds[j].push_back(edge); } - } - lower_bounds[j].reserve(lower_offsets[j + 1] - lower_offsets[j]); - for (i_t p = lower_offsets[j]; p < lower_offsets[j + 1]; p++) { - variable_bound_edge_t edge{lower_variables[p], lower_weights[p], lower_biases[p]}; - lower_bounds[j].push_back(edge); - if (is_binary_variable(edge.variable)) { binary_lower_bounds[j].push_back(edge); } - } - } - - std::vector propagated_upper_count(lp.num_cols, 0); - std::vector propagated_lower_count(lp.num_cols, 0); - i_t propagated_upper_edges = 0; - i_t propagated_lower_edges = 0; - - auto add_propagated_bound = [&](i_t j, - variable_bound_edge_t edge, - std::vector>& all_bounds, - std::vector>& binary_bounds, - std::vector& propagated_count, - i_t& propagated_edges) { - if (!valid_bound_edge(edge) || !is_binary_variable(edge.variable) || edge.variable == j) { - return false; - } - if (std::abs(edge.weight) <= bound_tol) { return false; } - if (propagated_count[j] >= max_propagated_bounds_per_variable) { return false; } - if (already_has_bound(binary_bounds[j], edge)) { return false; } - binary_bounds[j].push_back(edge); - all_bounds[j].push_back(edge); - propagated_count[j]++; - propagated_edges++; - return true; - }; - - // Compose affine bounds during row-to-SNF preprocessing. For example, - // - // z <= a w + b, w <= c x + d - // - // gives - // - // z <= a c x + (b + a d). - // - // If a is negative, the caller swaps upper/lower source bounds. - auto compose_bound = [](const variable_bound_edge_t& outer, const variable_bound_edge_t& inner) { - return variable_bound_edge_t{ - inner.variable, outer.weight * inner.weight, outer.bias + outer.weight * inner.bias}; - }; - - for (i_t round = 0; round < max_propagation_rounds; round++) { - bool changed = false; - for (i_t j = 0; j < lp.num_cols; j++) { - if (var_types[j] != variable_type_t::CONTINUOUS) { continue; } - - for (const auto& edge : continuous_upper_bounds[j]) { - if (!valid_bound_edge(edge) || std::abs(edge.weight) <= bound_tol) { continue; } - const auto& source_bounds = edge.weight > 0.0 ? binary_upper_bounds[edge.variable] - : binary_lower_bounds[edge.variable]; - for (const auto& source_bound : source_bounds) { - changed |= add_propagated_bound(j, - compose_bound(edge, source_bound), - upper_bounds, - binary_upper_bounds, - propagated_upper_count, - propagated_upper_edges); - } - } - - for (const auto& edge : continuous_lower_bounds[j]) { - if (!valid_bound_edge(edge) || std::abs(edge.weight) <= bound_tol) { continue; } - const auto& source_bounds = edge.weight > 0.0 ? binary_lower_bounds[edge.variable] - : binary_upper_bounds[edge.variable]; - for (const auto& source_bound : source_bounds) { - changed |= add_propagated_bound(j, - compose_bound(edge, source_bound), - lower_bounds, - binary_lower_bounds, - propagated_lower_count, - propagated_lower_edges); - } - } - } - if (!changed) { break; } - } - - auto flatten_bounds = [&](const std::vector>& bounds, - std::vector& offsets, - std::vector& variables, - std::vector& weights, - std::vector& biases) { - offsets.assign(lp.num_cols + 1, 0); - variables.clear(); - weights.clear(); - biases.clear(); - for (i_t j = 0; j < lp.num_cols; j++) { - offsets[j] = static_cast(variables.size()); - for (const auto& edge : bounds[j]) { - variables.push_back(edge.variable); - weights.push_back(edge.weight); - biases.push_back(edge.bias); - } - } - offsets[lp.num_cols] = static_cast(variables.size()); - }; - - if (propagated_upper_edges > 0 || propagated_lower_edges > 0) { - flatten_bounds(upper_bounds, upper_offsets, upper_variables, upper_weights, upper_biases); - flatten_bounds(lower_bounds, lower_offsets, lower_variables, lower_weights, lower_biases); - settings.log.printf( - "%d propagated variable upper bounds and %d propagated variable lower bounds\n", - propagated_upper_edges, - propagated_lower_edges); - } } template diff --git a/cpp/tests/mip/cuts_test.cu b/cpp/tests/mip/cuts_test.cu index 8e1271db60..d9ab0ad066 100644 --- a/cpp/tests/mip/cuts_test.cu +++ b/cpp/tests/mip/cuts_test.cu @@ -1489,43 +1489,6 @@ mps_parser::mps_data_model_t create_small_single_node_flow_problem( return problem; } -mps_parser::mps_data_model_t create_variable_bound_chain_problem() -{ - mps_parser::mps_data_model_t problem; - - // x0 is binary, x1 and x2 are continuous, and the rows imply - // x2 <= 2 * x1 <= 2 * x0. - std::vector offsets = {0, 2, 4}; - std::vector indices = {2, 1, 1, 0}; - std::vector coefficients = {1.0, -2.0, 1.0, -1.0}; - problem.set_csr_constraint_matrix(coefficients.data(), - coefficients.size(), - indices.data(), - indices.size(), - offsets.data(), - offsets.size()); - - std::vector lower_bounds(2, -std::numeric_limits::infinity()); - std::vector upper_bounds = {0.0, 0.0}; - problem.set_constraint_lower_bounds(lower_bounds.data(), lower_bounds.size()); - problem.set_constraint_upper_bounds(upper_bounds.data(), upper_bounds.size()); - - std::vector var_lower_bounds = {0.0, 0.0, 0.0}; - std::vector var_upper_bounds = { - 1.0, std::numeric_limits::infinity(), std::numeric_limits::infinity()}; - problem.set_variable_lower_bounds(var_lower_bounds.data(), var_lower_bounds.size()); - problem.set_variable_upper_bounds(var_upper_bounds.data(), var_upper_bounds.size()); - - std::vector objective_coefficients(3, 0.0); - problem.set_objective_coefficients(objective_coefficients.data(), objective_coefficients.size()); - - std::vector variable_types = {'I', 'C', 'C'}; - problem.set_variable_types(variable_types); - - problem.set_maximize(false); - return problem; -} - struct flow_cover_test_problem_t { raft::handle_t handle; dual_simplex::simplex_solver_settings_t settings; @@ -1692,28 +1655,4 @@ TEST(cuts, flow_cover_generates_valid_single_node_flow_cut) EXPECT_GT(generated_cuts, 0); } -TEST(cuts, variable_bounds_propagates_continuous_link_to_binary_bound) -{ - auto test_problem = build_flow_cover_test_problem(create_variable_bound_chain_problem()); - - dual_simplex::variable_bounds_t variable_bounds(test_problem.lp, - test_problem.settings, - test_problem.var_types, - test_problem.Arow, - test_problem.new_slacks); - - bool found_propagated_bound = false; - const int binary_col = 0; - const int continuous_col = 2; - for (int p = variable_bounds.upper_offsets[continuous_col]; - p < variable_bounds.upper_offsets[continuous_col + 1]; - p++) { - if (variable_bounds.upper_variables[p] != binary_col) { continue; } - found_propagated_bound |= std::abs(variable_bounds.upper_weights[p] - 2.0) <= 1e-9 && - std::abs(variable_bounds.upper_biases[p]) <= 1e-9; - } - - EXPECT_TRUE(found_propagated_bound); -} - } // namespace cuopt::linear_programming::test From 2908d725c44197492ac723c753ba27dd31ef640e Mon Sep 17 00:00:00 2001 From: Hugo Linsenmaier Date: Wed, 6 May 2026 13:08:48 -0700 Subject: [PATCH 15/17] Fix compile errors --- cpp/tests/mip/cuts_test.cu | 20 +++++++++----------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/cpp/tests/mip/cuts_test.cu b/cpp/tests/mip/cuts_test.cu index db99d78160..d68f7e84e5 100644 --- a/cpp/tests/mip/cuts_test.cu +++ b/cpp/tests/mip/cuts_test.cu @@ -32,6 +32,7 @@ #include #include #include +#include #include #include #include @@ -1420,25 +1421,22 @@ mps_parser::mps_data_model_t create_small_single_node_flow_problem( std::vector offsets = {0, 3, 5, 7, 9}; std::vector indices = {3, 4, 5, 0, 3, 1, 4, 2, 5}; std::vector coefficients = {1.0, 1.0, -1.0, -3.0, 1.0, -6.0, 1.0, -3.0, 1.0}; - problem.set_csr_constraint_matrix(coefficients.data(), - coefficients.size(), - indices.data(), - indices.size(), - offsets.data(), - offsets.size()); + problem.set_csr_constraint_matrix(std::span{coefficients}, + std::span{indices}, + std::span{offsets}); std::vector lower_bounds(4, -std::numeric_limits::infinity()); std::vector upper_bounds = {4.0, 0.0, 0.0, 0.0}; - problem.set_constraint_lower_bounds(lower_bounds.data(), lower_bounds.size()); - problem.set_constraint_upper_bounds(upper_bounds.data(), upper_bounds.size()); + problem.set_constraint_lower_bounds(std::span{lower_bounds}); + problem.set_constraint_upper_bounds(std::span{upper_bounds}); std::vector var_lower_bounds(6, 0.0); std::vector var_upper_bounds = {1.0, 1.0, 1.0, 3.0, 6.0, 3.0}; - problem.set_variable_lower_bounds(var_lower_bounds.data(), var_lower_bounds.size()); - problem.set_variable_upper_bounds(var_upper_bounds.data(), var_upper_bounds.size()); + problem.set_variable_lower_bounds(std::span{var_lower_bounds}); + problem.set_variable_upper_bounds(std::span{var_upper_bounds}); std::vector objective_coefficients(6, 0.0); - problem.set_objective_coefficients(objective_coefficients.data(), objective_coefficients.size()); + problem.set_objective_coefficients(std::span{objective_coefficients}); std::vector variable_types = {'I', 'I', 'I', 'C', 'C', 'C'}; problem.set_variable_types(variable_types); From 393cff1b4d30f89dcdb7b8b07d919e56afaad42a Mon Sep 17 00:00:00 2001 From: Hugo Linsenmaier Date: Mon, 11 May 2026 17:05:19 -0700 Subject: [PATCH 16/17] Fix indexing bug --- cpp/src/cuts/cuts.cpp | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index aacdeab3e9..617b7df00f 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -1031,7 +1031,7 @@ i_t knapsack_generation_t::generate_flow_cover_cut( bool absorbs_binary_coeff; // True if a direct x coefficient was folded into this arc. }; - auto is_binary_variable = [&](i_t j) { + auto is_zero_one_integer_variable = [&](i_t j) { return var_types[j] == variable_type_t::INTEGER && std::abs(lp.lower[j]) <= bound_tol && std::abs(lp.upper[j] - 1.0) <= bound_tol; }; @@ -1147,7 +1147,7 @@ i_t knapsack_generation_t::generate_flow_cover_cut( if (is_slack_[j]) { continue; } const f_t coeff = negate_row ? -row.coeff(k) : row.coeff(k); if (std::abs(coeff) <= tol) { continue; } - if (is_binary_variable(j)) { + if (is_zero_one_integer_variable(j)) { if (binary_coefficients.emplace(j, coeff).second) { binary_columns.push_back(j); } else { @@ -1217,7 +1217,7 @@ i_t knapsack_generation_t::generate_flow_cover_cut( const i_t end = variable_bounds.upper_offsets[j + 1]; for (i_t p = start; p < end; p++) { const i_t x_col = variable_bounds.upper_variables[p]; - if (!is_binary_variable(x_col)) { continue; } + if (!is_zero_one_integer_variable(x_col)) { continue; } const f_t gamma = variable_bounds.upper_weights[p]; const f_t alpha = variable_bounds.upper_biases[p]; if (!std::isfinite(gamma) || !std::isfinite(alpha) || lower_j <= -inf) { continue; } @@ -1277,7 +1277,7 @@ i_t knapsack_generation_t::generate_flow_cover_cut( const i_t end = variable_bounds.lower_offsets[j + 1]; for (i_t p = start; p < end; p++) { const i_t x_col = variable_bounds.lower_variables[p]; - if (!is_binary_variable(x_col)) { continue; } + if (!is_zero_one_integer_variable(x_col)) { continue; } const f_t gamma = variable_bounds.lower_weights[p]; const f_t alpha = variable_bounds.lower_biases[p]; if (!std::isfinite(gamma) || !std::isfinite(alpha) || upper_j >= inf) { continue; } @@ -3877,8 +3877,9 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, if (row_len < 2) { continue; } if (num_integer_in_row[i] < 1) { continue; } const f_t a_ij = lp.A.x[p]; - const f_t slack_lower = lp.lower[slack_map_[i]]; - const f_t slack_upper = lp.upper[slack_map_[i]]; + const i_t slack_index = slack_map_[i]; + const f_t slack_lower = slack_index >= 0 ? lp.lower[slack_index] : 0.0; + const f_t slack_upper = slack_index >= 0 ? lp.upper[slack_index] : 0.0; const f_t slack_coeff_i = slack_coeff[i]; const f_t sigma_slack_lower = slack_coeff_i > 0.0 ? slack_coeff_i * slack_lower : slack_coeff_i * slack_upper; @@ -3977,8 +3978,9 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, if (row_len < 2) { continue; } if (num_integer_in_row[i] < 1) { continue; } const f_t a_ij = lp.A.x[p]; - const f_t slack_lower = lp.lower[slack_map_[i]]; - const f_t slack_upper = lp.upper[slack_map_[i]]; + const i_t slack_index = slack_map_[i]; + const f_t slack_lower = slack_index >= 0 ? lp.lower[slack_index] : 0.0; + const f_t slack_upper = slack_index >= 0 ? lp.upper[slack_index] : 0.0; const f_t slack_coeff_i = slack_coeff[i]; const f_t sigma_slack_lower = slack_coeff_i > 0.0 ? slack_coeff_i * slack_lower : slack_coeff_i * slack_upper; From e0db214c07f95a0b0d721a747147da34b1c7f135 Mon Sep 17 00:00:00 2001 From: Hugo Linsenmaier Date: Mon, 11 May 2026 17:17:19 -0700 Subject: [PATCH 17/17] Merge knapsack paths --- cpp/src/cuts/cuts.cpp | 60 ++++++++++++++++--------------------------- cpp/src/cuts/cuts.hpp | 20 +++++++-------- 2 files changed, 31 insertions(+), 49 deletions(-) diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index 617b7df00f..d2c02cdb04 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -1507,7 +1507,8 @@ i_t knapsack_generation_t::generate_flow_cover_cut( if (values.empty()) { return -1; } std::vector solution; - const f_t objective = prefix_ratio_knapsack_problem(values, weights, cover_capacity, solution); + const f_t objective = greedy_knapsack_problem( + values, weights, cover_capacity, solution, greedy_knapsack_mode_t::STRICT_RATIO_PREFIX); if (std::isnan(objective)) { return -1; } static_cast(objective); @@ -2315,7 +2316,8 @@ template f_t knapsack_generation_t::greedy_knapsack_problem(const std::vector& values, const std::vector& weights, f_t rhs, - std::vector& solution) + std::vector& solution, + greedy_knapsack_mode_t mode) { i_t n = weights.size(); solution.assign(n, 0.0); @@ -2333,6 +2335,24 @@ f_t knapsack_generation_t::greedy_knapsack_problem(const std::vector ratios[j]; }); + if (mode == greedy_knapsack_mode_t::STRICT_RATIO_PREFIX) { + // Wolter Algorithm 5.1 for KPSNF^rat: take the ratio-sorted prefix while + // the strict capacity remains satisfied and stop at the first item that + // does not fit. + f_t total_weight = 0.0; + f_t total_value = 0.0; + for (i_t item : perm) { + if (total_weight + weights[item] < rhs) { + solution[item] = 1.0; + total_weight += weights[item]; + total_value += values[item]; + } else { + break; + } + } + return total_value; + } + // Greedy select items with the best value / weight ratio until the remaining capacity is // exhausted f_t remaining = rhs; @@ -2366,42 +2386,6 @@ f_t knapsack_generation_t::greedy_knapsack_problem(const std::vector -f_t knapsack_generation_t::prefix_ratio_knapsack_problem(const std::vector& values, - const std::vector& weights, - f_t strict_rhs, - std::vector& solution) -{ - const i_t n = static_cast(weights.size()); - solution.assign(n, 0.0); - if (n == 0) { return static_cast(0.0); } - - // Wolter Algorithm 5.1 for KPSNF^rat: sort by p_j / u_j, add items while - // the strict capacity remains satisfied, and stop at the first item that - // does not fit. This is intentionally not the generic greedy knapsack - // heuristic, which keeps scanning for smaller later items. - std::vector permutation(n); - std::iota(permutation.begin(), permutation.end(), 0); - std::stable_sort(permutation.begin(), permutation.end(), [&](i_t a, i_t b) { - const f_t ratio_a = values[a] / weights[a]; - const f_t ratio_b = values[b] / weights[b]; - return ratio_a > ratio_b; - }); - - f_t total_weight = 0.0; - f_t total_value = 0.0; - for (i_t item : permutation) { - if (total_weight + weights[item] < strict_rhs) { - solution[item] = 1.0; - total_weight += weights[item]; - total_value += values[item]; - } else { - break; - } - } - return total_value; -} - template f_t knapsack_generation_t::solve_knapsack_problem(const std::vector& values, const std::vector& weights, diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index 41c682e64f..44ba48dadf 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -388,17 +388,15 @@ class knapsack_generation_t { 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, - f_t rhs, - std::vector& solution); - - // Generate a ratio-sorted prefix solution to a strict 0-1 knapsack problem. - f_t prefix_ratio_knapsack_problem(const std::vector& values, - const std::vector& weights, - f_t strict_rhs, - std::vector& solution); + enum class greedy_knapsack_mode_t { SCAN_ALL_WITH_BEST_SINGLE, STRICT_RATIO_PREFIX }; + + // Generate a heuristic solution to the 0-1 knapsack problem. + f_t greedy_knapsack_problem( + const std::vector& values, + const std::vector& weights, + f_t rhs, + std::vector& solution, + greedy_knapsack_mode_t mode = greedy_knapsack_mode_t::SCAN_ALL_WITH_BEST_SINGLE); // Solve a 0-1 knapsack problem using dynamic programming f_t solve_knapsack_problem(const std::vector& values,