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/include/cuopt/linear_programming/constants.h b/cpp/include/cuopt/linear_programming/constants.h index b251b3eaba..5fb0779e90 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 ae0187e454..e548f2d500 100644 --- a/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp +++ b/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp @@ -104,6 +104,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 0b93ece0c7..25615fcfdf 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -12,9 +12,11 @@ #include #include +#include #include #include #include +#include #include #include @@ -821,6 +823,89 @@ void cut_pool_t::drop_cuts() // TODO: Implement this } +enum class flow_cover_bound_side_t : int8_t { UPPER, LOWER }; + +template +f_t flow_cover_scaled_primal_tol(const flow_cover_context_t& context, f_t scale) +{ + return context.settings.primal_tol * std::max(1.0, std::abs(scale)); +} + +template +bool flow_cover_is_zero_one_integer_variable(const flow_cover_context_t& context, i_t j) +{ + const f_t bound_tol = context.settings.primal_tol; + return context.var_types[j] == variable_type_t::INTEGER && + std::abs(context.lp.lower[j]) <= bound_tol && + std::abs(context.lp.upper[j] - 1.0) <= bound_tol; +} + +template +f_t flow_cover_arc_tol(const flow_cover_context_t& context, + const snf_arc_t& arc) +{ + return flow_cover_scaled_primal_tol(context, arc.u); +} + +template +snf_arc_t flow_cover_build_arc(const flow_cover_context_t& context, + 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) +{ + auto clamp01 = [](f_t x) { return std::min(1.0, std::max(0.0, x)); }; + const f_t x_value = x_col >= 0 ? clamp01(context.xstar[x_col]) : fixed_x; + f_t y_value = y_const; + if (y_col >= 0) { y_value += y_coeff * context.xstar[y_col]; } + if (x_col >= 0) { + y_value += y_x_coeff * context.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}; +} + +template +void flow_cover_try_add_candidate(const flow_cover_context_t& context, + const flow_cover_arc_spec_t& spec, + f_t y_star, + std::vector>& candidates) +{ + const f_t bound_tol = context.settings.primal_tol; + const f_t feasibility_tol = context.settings.primal_tol; + if (spec.u < -bound_tol) { return; } + const f_t capacity = std::max(0.0, spec.u); + if (capacity <= feasibility_tol) { return; } + + snf_candidate_t candidate; + candidate.arc = flow_cover_build_arc(context, + capacity, + spec.in_n2, + spec.x_col, + spec.fixed_x, + spec.y_const, + spec.y_col, + spec.y_coeff, + spec.y_x_coeff); + candidate.b_shift = spec.b_shift; + candidate.distance = std::abs(spec.active_bound - y_star); + candidate.absorbs_binary_coeff = spec.absorbs_binary_coeff; + candidates.push_back(candidate); +} + +template +f_t flow_cover_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); +} + template knapsack_generation_t::knapsack_generation_t( const lp_problem_t& lp, @@ -898,6 +983,626 @@ 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 + + // 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; } + 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]]) { + slack_col = Arow.j[p]; + slack_coeff = Arow.x[p]; + slack_count++; + } + } + if (slack_count > 1 || (slack_count == 1 && std::abs(slack_coeff) <= settings.zero_tol)) { + continue; + } + if (slack_count == 0) { + flow_cover_constraints_.push_back({i, false}); + flow_cover_constraints_.push_back({i, true}); + 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, false}); } + if (sigma_slack_upper < inf) { flow_cover_constraints_.push_back({i, true}); } + } +} + +template +bool knapsack_generation_t::normalize_row_side( + const flow_cover_context_t& context, + const flow_cover_row_t& flow_cover_row, + inequality_t& row, + f_t& b, + bool& negate_row) +{ + const f_t coefficient_tol = static_cast(1e-6); + auto& scratch = flow_cover_scratch_; + row = + inequality_t(context.Arow, flow_cover_row.row, context.lp.rhs[flow_cover_row.row]); + + i_t slack_count = 0; + i_t slack_col = -1; + f_t slack_coeff = 0.0; + 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 false; } + + negate_row = flow_cover_row.reverse; + b = negate_row ? -row.rhs : row.rhs; + if (slack_count == 1) { + if (std::abs(slack_coeff) <= coefficient_tol) { return false; } + const f_t sigma_slack_lower = slack_coeff > 0.0 ? slack_coeff * context.lp.lower[slack_col] + : slack_coeff * context.lp.upper[slack_col]; + const f_t sigma_slack_upper = slack_coeff > 0.0 ? slack_coeff * context.lp.upper[slack_col] + : slack_coeff * context.lp.lower[slack_col]; + if (!flow_cover_row.reverse) { + if (sigma_slack_lower <= -inf) { return false; } + negate_row = false; + b = row.rhs - sigma_slack_lower; + } else { + if (sigma_slack_upper >= inf) { return false; } + negate_row = true; + b = -row.rhs + sigma_slack_upper; + } + } + if (!std::isfinite(b)) { return false; } + + 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 * context.xstar[j]; + row_side_scale += std::abs(coeff * context.xstar[j]); + } + return row_side_activity <= b + flow_cover_scaled_primal_tol(context, row_side_scale); + }(), + "Flow cover normalized row side excludes LP solution"); + + scratch.continuous_terms.reserve(row_len); + scratch.binary_columns.reserve(row_len); + auto add_binary_coeff = [&](i_t j, f_t coeff) { + if (!scratch.binary_coefficients_touched[j]) { + scratch.binary_coefficients_touched[j] = 1; + scratch.touched_binary_coefficients.push_back(j); + scratch.binary_columns.push_back(j); + } + scratch.binary_coefficients[j] += coeff; + }; + + 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); + if (std::abs(coeff) <= coefficient_tol) { continue; } + if (flow_cover_is_zero_one_integer_variable(context, j)) { + add_binary_coeff(j, coeff); + } else { + scratch.continuous_terms.push_back({j, coeff}); + } + } + + return true; +} + +template +bool knapsack_generation_t::build_snf_relaxation( + const flow_cover_context_t& context, f_t b, f_t& snf_b) +{ + auto& scratch = flow_cover_scratch_; + const f_t coefficient_tol = static_cast(1e-6); + const f_t feasibility_tol = context.settings.primal_tol; + const f_t bound_tol = context.settings.primal_tol; + f_t b_shift = 0.0; + + scratch.arcs.reserve(scratch.continuous_terms.size() + scratch.binary_columns.size()); + + auto add_variable_bound_candidates = [&](i_t j, f_t c, flow_cover_bound_side_t side) { + const bool use_upper_bound = side == flow_cover_bound_side_t::UPPER; + const f_t lower_j = context.lp.lower[j]; + const f_t upper_j = context.lp.upper[j]; + if (use_upper_bound && lower_j <= -inf) { return; } + if (!use_upper_bound && upper_j >= inf) { return; } + + const i_t start = use_upper_bound ? context.variable_bounds.upper_offsets[j] + : context.variable_bounds.lower_offsets[j]; + const i_t end = use_upper_bound ? context.variable_bounds.upper_offsets[j + 1] + : context.variable_bounds.lower_offsets[j + 1]; + const f_t endpoint = use_upper_bound ? lower_j : upper_j; + + for (i_t p = start; p < end; p++) { + const i_t x_col = use_upper_bound ? context.variable_bounds.upper_variables[p] + : context.variable_bounds.lower_variables[p]; + if (!flow_cover_is_zero_one_integer_variable(context, x_col)) { continue; } + const f_t gamma = use_upper_bound ? context.variable_bounds.upper_weights[p] + : context.variable_bounds.lower_weights[p]; + const f_t alpha = use_upper_bound ? context.variable_bounds.upper_biases[p] + : context.variable_bounds.lower_biases[p]; + if (!std::isfinite(gamma) || !std::isfinite(alpha)) { continue; } + + const f_t direct_coeff = scratch.binary_coefficients_touched[x_col] + ? scratch.binary_coefficients[x_col] + : static_cast(0.0); + const std::array a_values = {direct_coeff, 0.0}; + const i_t num_a_values = std::abs(direct_coeff) > coefficient_tol ? 2 : 1; + for (i_t h = 0; h < num_a_values; h++) { + const f_t a = a_values[h]; + const bool in_n2 = use_upper_bound ? c < 0.0 : c > 0.0; + const f_t signed_capacity = c * gamma + a; + const f_t endpoint_term = c * (endpoint - alpha); + const bool valid_endpoint = + in_n2 ? (endpoint_term <= bound_tol && endpoint_term + a <= bound_tol && + signed_capacity <= bound_tol) + : (endpoint_term >= -bound_tol && endpoint_term + a >= -bound_tol && + signed_capacity >= -bound_tol); + if (!valid_endpoint) { continue; } + + flow_cover_arc_spec_t spec; + spec.u = in_n2 ? -signed_capacity : signed_capacity; + spec.in_n2 = in_n2; + spec.x_col = x_col; + spec.fixed_x = 0.0; + spec.y_const = in_n2 ? c * alpha : -c * alpha; + spec.y_col = j; + spec.y_coeff = in_n2 ? -c : c; + spec.y_x_coeff = in_n2 ? -a : a; + spec.b_shift = c * alpha; + spec.active_bound = gamma * context.xstar[x_col] + alpha; + spec.absorbs_binary_coeff = std::abs(a) > coefficient_tol; + flow_cover_try_add_candidate(context, spec, context.xstar[j], scratch.candidates); + } + } + }; + + auto add_simple_bound_candidates = [&](i_t j, f_t c) { + const f_t lower_j = context.lp.lower[j]; + const f_t upper_j = context.lp.upper[j]; + if (lower_j <= -inf || upper_j >= inf) { return; } + const f_t capacity = std::abs(c) * (upper_j - lower_j); + if (capacity <= feasibility_tol) { return; } + + auto add_simple_candidate = + [&](bool in_n2, f_t y_const, f_t y_coeff, f_t shift, f_t active_bound) { + flow_cover_arc_spec_t spec; + spec.u = capacity; + spec.in_n2 = in_n2; + spec.x_col = -1; + spec.fixed_x = 1.0; + spec.y_const = y_const; + spec.y_col = j; + spec.y_coeff = y_coeff; + spec.y_x_coeff = 0.0; + spec.b_shift = shift; + spec.active_bound = active_bound; + spec.absorbs_binary_coeff = false; + flow_cover_try_add_candidate(context, spec, context.xstar[j], scratch.candidates); + }; + + if (c > 0.0) { + add_simple_candidate(false, -c * lower_j, c, c * lower_j, upper_j); + add_simple_candidate(true, c * upper_j, -c, c * upper_j, lower_j); + } else { + add_simple_candidate(true, c * lower_j, -c, c * lower_j, upper_j); + add_simple_candidate(false, -c * upper_j, c, c * upper_j, lower_j); + } + }; + + for (const auto& [j, c] : scratch.continuous_terms) { + scratch.candidates.clear(); + add_variable_bound_candidates(j, c, flow_cover_bound_side_t::UPPER); + add_variable_bound_candidates(j, c, flow_cover_bound_side_t::LOWER); + add_simple_bound_candidates(j, c); + + if (scratch.candidates.empty()) { return false; } + auto best = + std::min_element(scratch.candidates.begin(), + scratch.candidates.end(), + [](const snf_candidate_t& a, const snf_candidate_t& b) { + return a.distance < b.distance; + }); + const f_t arc_lower_tol = + std::max(10 * context.settings.primal_tol, static_cast(1e-5)); + const f_t arc_upper_tol = + std::max(100 * context.settings.primal_tol, static_cast(1e-4)) * + std::max(1.0, best->arc.u); + if (best->arc.y_value < -arc_lower_tol || + best->arc.y_value > best->arc.u * best->arc.x_value + arc_upper_tol) { + return false; + } + if (best->arc.y_value < 0.0) { best->arc.y_value = 0.0; } + scratch.arcs.push_back(best->arc); + b_shift += best->b_shift; + if (best->absorbs_binary_coeff && best->arc.x_col >= 0) { + scratch.binary_coefficients[best->arc.x_col] = 0.0; + } + } + + for (i_t j : scratch.binary_columns) { + const f_t coeff = scratch.binary_coefficients[j]; + if (std::abs(coeff) <= coefficient_tol) { continue; } + const f_t u = std::abs(coeff); + scratch.arcs.push_back(flow_cover_build_arc(context, u, coeff < 0.0, j, 0.0, 0.0, -1, 0.0, u)); + } + + if (scratch.arcs.empty()) { return false; } + 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 : scratch.arcs) { + const f_t arc_tol = flow_cover_arc_tol(context, arc); + 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 + context.settings.primal_tol * snf_scale; + }(), + "Flow cover SNF relaxation excludes LP solution"); + + return true; +} + +template +bool knapsack_generation_t::select_flow_cover( + const flow_cover_context_t& context, f_t snf_b, flow_cover_selection_t& selection) +{ + auto& scratch = flow_cover_scratch_; + const f_t feasibility_tol = context.settings.primal_tol; + + bool has_binary_controlled_arc = false; + bool has_n1 = false; + f_t sum_n1_u = 0.0; + for (const auto& arc : scratch.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; + if (!has_binary_controlled_arc || !has_n1 || cover_capacity <= feasibility_tol) { return false; } + + scratch.values.reserve(scratch.arcs.size()); + scratch.weights.reserve(scratch.arcs.size()); + scratch.item_to_arc.reserve(scratch.arcs.size()); + for (i_t k = 0; k < static_cast(scratch.arcs.size()); k++) { + const auto& arc = scratch.arcs[k]; + if (arc.u <= feasibility_tol) { continue; } + const f_t value = arc.in_n2 ? arc.x_value : 1.0 - arc.x_value; + scratch.values.push_back(std::max(0.0, value)); + scratch.weights.push_back(arc.u); + scratch.item_to_arc.push_back(k); + } + if (scratch.values.empty()) { return false; } + + const f_t objective = greedy_knapsack_problem(scratch.values, + scratch.weights, + cover_capacity, + scratch.solution, + greedy_knapsack_mode_t::STRICT_RATIO_PREFIX); + if (std::isnan(objective)) { return false; } + + scratch.in_c1.assign(scratch.arcs.size(), 0); + scratch.in_c2.assign(scratch.arcs.size(), 0); + for (i_t item = 0; item < static_cast(scratch.solution.size()); item++) { + const i_t arc_index = scratch.item_to_arc[item]; + if (scratch.arcs[arc_index].in_n2) { + if (scratch.solution[item] > 0.5) { scratch.in_c2[arc_index] = 1; } + } else { + if (scratch.solution[item] <= 0.5) { scratch.in_c1[arc_index] = 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(scratch.arcs.size()); k++) { + if (scratch.in_c1[k]) { + sum_c1_u += scratch.arcs[k].u; + c1_nonempty = true; + } + if (scratch.in_c2[k]) { sum_c2_u += scratch.arcs[k].u; } + } + + selection.lambda = -snf_b + sum_c1_u - sum_c2_u; + return c1_nonempty && selection.lambda > feasibility_tol; +} + +template +flow_cover_evaluation_t knapsack_generation_t::evaluate_cmirfci( + const flow_cover_context_t& context, f_t snf_b, f_t lambda) +{ + auto& scratch = flow_cover_scratch_; + constexpr f_t min_mir_beta_fraction = 0.01; + constexpr f_t max_mir_beta_fraction = 0.95; + const f_t min_violation = static_cast(1e-6); + const f_t lambda_tol = min_violation; + f_t max_c1_ltilde2 = -inf; + f_t max_c1 = -inf; + f_t max_u_ge_lambda = -inf; + + auto add_ubar_candidate = [&](f_t candidate) { + if (candidate <= lambda + lambda_tol || !std::isfinite(candidate)) { return; } + for (f_t existing : scratch.ubar_candidates) { + const f_t duplicate_tol = static_cast(1e-8) * std::max(1.0, std::abs(candidate)); + if (std::abs(existing - candidate) <= duplicate_tol) { return; } + } + scratch.ubar_candidates.push_back(candidate); + }; + + for (i_t k = 0; k < static_cast(scratch.arcs.size()); k++) { + const auto& arc = scratch.arcs[k]; + if (arc.u > lambda + lambda_tol) { add_ubar_candidate(arc.u); } + if (arc.u >= lambda - lambda_tol) { max_u_ge_lambda = std::max(max_u_ge_lambda, arc.u); } + if (!arc.in_n2 && scratch.in_c1[k] && arc.u > lambda + lambda_tol) { + max_c1 = std::max(max_c1, arc.u); + max_c1_ltilde2 = std::max(max_c1_ltilde2, arc.u); + } + if (arc.in_n2 && !scratch.in_c2[k] && arc.u > lambda + lambda_tol && + std::min(arc.u, lambda) * arc.x_value <= arc.y_value + min_violation) { + 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); + + flow_cover_evaluation_t best{0.0, 0.0}; + if (scratch.ubar_candidates.empty()) { return best; } + + scratch.best_in_l1.assign(scratch.arcs.size(), 0); + scratch.best_in_l2.assign(scratch.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_beta = fractional_part(-lambda / ubar); + const f_t f_pos = flow_cover_mir_f(f_beta, arc.u / ubar); + const f_t f_neg = flow_cover_mir_f(f_beta, -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 (f_t ubar : scratch.ubar_candidates) { + const f_t beta = -lambda / ubar; + const f_t f_beta = fractional_part(beta); + if (f_beta < min_mir_beta_fraction || f_beta > max_mir_beta_fraction) { continue; } + + scratch.in_l1.assign(scratch.arcs.size(), 0); + scratch.in_l2.assign(scratch.arcs.size(), 0); + for (i_t k = 0; k < static_cast(scratch.arcs.size()); k++) { + const auto& arc = scratch.arcs[k]; + const f_t f_pos = flow_cover_mir_f(f_beta, arc.u / ubar); + const f_t f_neg = flow_cover_mir_f(f_beta, -arc.u / ubar); + if (!arc.in_n2 && !scratch.in_c1[k] && + arc.y_value - (arc.u - lambda * f_pos) * arc.x_value >= -min_violation) { + scratch.in_l1[k] = 1; + } + if (arc.in_n2 && !scratch.in_c2[k] && + lambda * f_neg * arc.x_value >= -arc.y_value - min_violation) { + scratch.in_l2[k] = 1; + } + } + + f_t lhs_value = 0.0; + for (i_t k = 0; k < static_cast(scratch.arcs.size()); k++) { + lhs_value += contribution(scratch.arcs[k], + scratch.in_c1[k], + scratch.in_c2[k], + scratch.in_l1[k], + scratch.in_l2[k], + ubar); + } + const f_t violation = lhs_value - snf_b; + if (violation > best.violation + min_violation) { + best.violation = violation; + best.ubar = ubar; + scratch.best_in_l1 = scratch.in_l1; + scratch.best_in_l2 = scratch.in_l2; + } + } + + return best; +} + +template +flow_cover_evaluation_t knapsack_generation_t::evaluate_sgfci( + const flow_cover_context_t& context, f_t snf_b, f_t lambda) +{ + auto& scratch = flow_cover_scratch_; + const f_t min_violation = static_cast(1e-6); + scratch.sgfci_in_l2.assign(scratch.arcs.size(), 0); + + f_t lhs_value = 0.0; + for (i_t k = 0; k < static_cast(scratch.arcs.size()); k++) { + const auto& arc = scratch.arcs[k]; + if (scratch.in_c1[k]) { + 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 (scratch.in_c2[k]) { + lhs_value -= arc.u; + } else if (std::min(arc.u, lambda) * arc.x_value <= arc.y_value + min_violation) { + scratch.sgfci_in_l2[k] = 1; + lhs_value -= std::min(arc.u, lambda) * arc.x_value; + } else { + lhs_value -= arc.y_value; + } + } + + return {lhs_value - snf_b, 0.0}; +} + +template +bool knapsack_generation_t::emit_flow_cover_cut( + const flow_cover_context_t& context, + f_t snf_b, + f_t lambda, + const flow_cover_evaluation_t& cmirfci, + const flow_cover_evaluation_t& sgfci, + inequality_t& cut) +{ + auto& scratch = flow_cover_scratch_; + const f_t accumulator_tol = static_cast(1e-12); + const f_t output_drop_tol = static_cast(1e-9); + const bool use_cmirfci = cmirfci.violation >= sgfci.violation; + f_t lhs_constant = 0.0; + + scratch.lhs_indices.reserve(scratch.arcs.size() * 2); + auto add_lhs_coeff = [&](i_t j, f_t value) { + if (j < 0 || std::abs(value) <= accumulator_tol) { return; } + if (!scratch.lhs_coefficients_touched[j]) { + scratch.lhs_coefficients_touched[j] = 1; + scratch.touched_lhs_coefficients.push_back(j); + scratch.lhs_indices.push_back(j); + } + scratch.lhs_coefficients[j] += value; + }; + + 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; + } + }; + + 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; + } + }; + + auto add_one_minus_x = [&](const snf_arc_t& arc, f_t multiplier) { + lhs_constant += multiplier; + add_x(arc, -multiplier); + }; + + if (use_cmirfci) { + const f_t f_beta_best = fractional_part(-lambda / cmirfci.ubar); + for (i_t k = 0; k < static_cast(scratch.arcs.size()); k++) { + const auto& arc = scratch.arcs[k]; + const f_t f_pos = flow_cover_mir_f(f_beta_best, arc.u / cmirfci.ubar); + const f_t f_neg = flow_cover_mir_f(f_beta_best, -arc.u / cmirfci.ubar); + if (scratch.in_c1[k]) { + add_y(arc, 1.0); + add_one_minus_x(arc, arc.u + lambda * f_neg); + } else if (!arc.in_n2) { + if (scratch.best_in_l1[k]) { + add_y(arc, 1.0); + add_x(arc, -(arc.u - lambda * f_pos)); + } + } else if (scratch.in_c2[k]) { + lhs_constant -= arc.u; + add_one_minus_x(arc, lambda * f_pos); + } else if (scratch.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(scratch.arcs.size()); k++) { + const auto& arc = scratch.arcs[k]; + if (scratch.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 (scratch.in_c2[k]) { + lhs_constant -= arc.u; + } else if (scratch.sgfci_in_l2[k]) { + add_x(arc, -std::min(arc.u, lambda)); + } else { + add_y(arc, -1.0); + } + } + } + + cut.clear(); + cut.reserve(scratch.lhs_indices.size()); + for (i_t j : scratch.lhs_indices) { + const f_t coeff = scratch.lhs_coefficients[j]; + if (std::abs(coeff) > output_drop_tol) { cut.push_back(j, -coeff); } + } + if (cut.size() == 0) { return false; } + cut.rhs = lhs_constant - snf_b; + cut.sort(); + + const f_t dot = cut.vector.dot(context.xstar); + const f_t violation_tol = + std::max(context.settings.primal_tol, static_cast(1e-6)); + return dot < cut.rhs - violation_tol; +} + +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 variable_bounds_t& variable_bounds, + const std::vector& var_types, + const std::vector& xstar, + const flow_cover_row_t& flow_cover_row, + inequality_t& cut) +{ + flow_cover_context_t context{lp, settings, Arow, variable_bounds, var_types, xstar}; + flow_cover_scratch_.clear(lp.num_cols); + + inequality_t row; + f_t b = 0.0; + bool negate_row = false; + if (!normalize_row_side(context, flow_cover_row, row, b, negate_row)) { return -1; } + + f_t snf_b = 0.0; + if (!build_snf_relaxation(context, b, snf_b)) { return -1; } + + flow_cover_selection_t selection{0.0}; + if (!select_flow_cover(context, snf_b, selection)) { return -1; } + + const auto cmirfci = evaluate_cmirfci(context, snf_b, selection.lambda); + const auto sgfci = evaluate_sgfci(context, snf_b, selection.lambda); + const f_t violation_tol = std::max(settings.primal_tol, static_cast(1e-6)); + if (cmirfci.violation <= violation_tol && sgfci.violation <= violation_tol) { return -1; } + + if (!emit_flow_cover_cut(context, snf_b, selection.lambda, cmirfci, sgfci, cut)) { return -1; } + return 0; } template @@ -1419,7 +2124,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); @@ -1437,6 +2143,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; @@ -1803,6 +2527,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, 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); + } + } + // Generate MIR and CG cuts if (settings.mir_cuts != 0 || settings.strong_chvatal_gomory_cuts != 0) { f_t cut_start_time = tic(); @@ -1860,6 +2594,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& var_types, + const std::vector& xstar, + variable_bounds_t& variable_bounds, + f_t start_time) +{ + if (knapsack_generation_.num_flow_cover_constraints() > 0) { + for (const auto& 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, variable_bounds, 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, @@ -2807,6 +3562,19 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, } f_t start_time = tic(); + struct variable_bound_edge_t { + i_t variable; + f_t weight; + f_t bias; + }; + + // This is shared variable-bound preprocessing used by multiple cut separators. + // Flow-cover applies its own 0-1 controller filter when it consumes these bounds. + auto is_variable_bound_controller = [&](i_t j) { + return var_types[j] != variable_type_t::CONTINUOUS; + }; + std::vector num_noncontinuous_in_row(lp.num_rows, 0); + // Construct the slack map slack_map_.resize(lp.num_rows, -1); std::vector slack_coeff(lp.num_rows, 0.0); @@ -2823,7 +3591,6 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, // The constraints are in the form: // sum_j a_j x_j + sigma * slack = beta - std::vector num_integer_in_row(lp.num_rows, 0); // Compute the upper activities of the constraints for (i_t i = 0; i < lp.num_rows; i++) { const i_t row_start = Arow.row_start[i]; @@ -2833,6 +3600,7 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, for (i_t p = row_start; p < row_end; p++) { const i_t j = Arow.j[p]; if (j == slack_index) { continue; } + if (is_variable_bound_controller(j)) { num_noncontinuous_in_row[i]++; } const f_t aj = Arow.x[p]; const f_t uj = lp.upper[j]; const f_t lj = lp.lower[j]; @@ -2850,8 +3618,6 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, num_pos_inf_[i]++; } } - - if (var_types[j] == variable_type_t::INTEGER) { num_integer_in_row[i]++; } } upper_activities_[i] = activity; } @@ -2894,18 +3660,21 @@ 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_noncontinuous_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 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 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; + 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; @@ -2916,14 +3685,9 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, if (a_ij > 0.0 && num_neg_inf_[i] <= 2) { const f_t lower_activity_j = lower_activity(lp.lower[j], lp.upper[j], a_ij); - // This is inefficient if num_neg_inf_[i] > 0 - // If num_neg_inf_[i] == 1 and var_types[s] != INTEGER, we can't derive a bound - // If num_neg_inf_[i] == 2 and var_types[s ^ j] != INTEGER, we can't derive a bound - // If num_neg_inf_[i] == 2 and var_types[s ^ j] == INTEGER, and lower_activity_j != -inf, - // 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 (!is_variable_bound_controller(l)) { 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]; @@ -2934,9 +3698,10 @@ 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); + const variable_bound_edge_t edge{l, -a_il / a_ij, beta / a_ij - (1.0 / a_ij) * sum}; + upper_variables.push_back(edge.variable); + upper_weights.push_back(edge.weight); + upper_biases.push_back(edge.bias); upper_edges++; } } @@ -2954,7 +3719,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 (!is_variable_bound_controller(l)) { 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]; @@ -2965,9 +3730,10 @@ 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); + const variable_bound_edge_t edge{l, -a_il / a_ij, beta / a_ij - (1.0 / a_ij) * sum}; + upper_variables.push_back(edge.variable); + upper_weights.push_back(edge.weight); + upper_biases.push_back(edge.bias); upper_edges++; } } @@ -2987,17 +3753,20 @@ 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_noncontinuous_in_row[i] < 1) { 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 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 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; + 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; @@ -3010,7 +3779,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 (!is_variable_bound_controller(l)) { 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) @@ -3022,9 +3791,10 @@ 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); + const variable_bound_edge_t edge{l, -a_il / a_ij, beta / a_ij - (1.0 / a_ij) * sum}; + lower_variables.push_back(edge.variable); + lower_weights.push_back(edge.weight); + lower_biases.push_back(edge.bias); lower_edges++; } } @@ -3042,7 +3812,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 (!is_variable_bound_controller(l)) { 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]; @@ -3053,9 +3823,10 @@ 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); + const variable_bound_edge_t edge{l, -a_il / a_ij, beta / a_ij - (1.0 / a_ij) * sum}; + lower_variables.push_back(edge.variable); + lower_weights.push_back(edge.weight); + lower_biases.push_back(edge.bias); lower_edges++; } } @@ -3639,7 +4410,7 @@ void complemented_mixed_integer_rounding_cut_t::transform_inequality( scratch_pad_.get_pad(inequality.vector.i, inequality.vector.x); // At this point we have converted all the continuous variables to be nonnegative // Note that since continuous variables had VUB or VLB, they modified - // the integer variables. + // the controller variables. // We clear the scratch pad. As it is no longer needed. scratch_pad_.clear_pad(); diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index 2d2a2dcd21..f29dd5c026 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -20,6 +20,7 @@ #include #include #include +#include #include #include @@ -39,7 +40,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 +180,8 @@ struct cut_info_t { "Knapsack ", "Strong CG ", "Clique ", - "Implied Bounds"}; + "Implied Bounds", + "Flow Cover "}; std::array num_cuts = {0}; }; @@ -327,6 +330,142 @@ class cut_pool_t { const f_t min_cut_distance_{1e-4}; }; +template +class variable_bounds_t; + +template +struct flow_cover_row_t { + i_t row; + bool reverse; +}; + +template +struct snf_arc_t { + 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. + f_t y_const; + i_t y_col; + f_t y_coeff; + f_t y_x_coeff; + f_t y_value; // Current LP value of the synthetic y_j. +}; + +template +struct snf_candidate_t { + snf_arc_t arc; + f_t b_shift; + f_t distance; + bool absorbs_binary_coeff; +}; + +template +struct flow_cover_arc_spec_t { + 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 b_shift; + f_t active_bound; + bool absorbs_binary_coeff; +}; + +template +struct flow_cover_context_t { + const lp_problem_t& lp; + const simplex_solver_settings_t& settings; + csr_matrix_t& Arow; + const variable_bounds_t& variable_bounds; + const std::vector& var_types; + const std::vector& xstar; +}; + +template +struct flow_cover_scratch_t { + std::vector> continuous_terms; + std::vector binary_columns; + std::vector binary_coefficients; + std::vector binary_coefficients_touched; + std::vector touched_binary_coefficients; + std::vector> arcs; + std::vector> candidates; + std::vector values; + std::vector weights; + std::vector item_to_arc; + std::vector solution; + std::vector in_c1; + std::vector in_c2; + std::vector ubar_candidates; + std::vector in_l1; + std::vector in_l2; + std::vector best_in_l1; + std::vector best_in_l2; + std::vector sgfci_in_l2; + std::vector lhs_indices; + std::vector lhs_coefficients; + std::vector lhs_coefficients_touched; + std::vector touched_lhs_coefficients; + + void clear(i_t num_cols) + { + continuous_terms.clear(); + binary_columns.clear(); + arcs.clear(); + candidates.clear(); + values.clear(); + weights.clear(); + item_to_arc.clear(); + solution.clear(); + in_c1.clear(); + in_c2.clear(); + ubar_candidates.clear(); + in_l1.clear(); + in_l2.clear(); + best_in_l1.clear(); + best_in_l2.clear(); + sgfci_in_l2.clear(); + lhs_indices.clear(); + + if (static_cast(binary_coefficients.size()) < num_cols) { + binary_coefficients.assign(num_cols, 0.0); + binary_coefficients_touched.assign(num_cols, 0); + } else { + for (i_t j : touched_binary_coefficients) { + binary_coefficients[j] = 0.0; + binary_coefficients_touched[j] = 0; + } + } + touched_binary_coefficients.clear(); + + if (static_cast(lhs_coefficients.size()) < num_cols) { + lhs_coefficients.assign(num_cols, 0.0); + lhs_coefficients_touched.assign(num_cols, 0); + } else { + for (i_t j : touched_lhs_coefficients) { + lhs_coefficients[j] = 0.0; + lhs_coefficients_touched[j] = 0; + } + } + touched_lhs_coefficients.clear(); + } +}; + +template +struct flow_cover_selection_t { + f_t lambda; +}; + +template +struct flow_cover_evaluation_t { + f_t violation; + f_t ubar; +}; + template class knapsack_generation_t { public: @@ -345,9 +484,24 @@ 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 variable_bounds_t& variable_bounds, + const std::vector& var_types, + const std::vector& xstar, + const flow_cover_row_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) { @@ -370,11 +524,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); + 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, @@ -387,12 +545,41 @@ class knapsack_generation_t { f_t rhs, std::vector& solution); + bool normalize_row_side(const flow_cover_context_t& context, + const flow_cover_row_t& flow_cover_row, + inequality_t& row, + f_t& b, + bool& negate_row); + + bool build_snf_relaxation(const flow_cover_context_t& context, f_t b, f_t& snf_b); + + bool select_flow_cover(const flow_cover_context_t& context, + f_t snf_b, + flow_cover_selection_t& selection); + + flow_cover_evaluation_t evaluate_cmirfci(const flow_cover_context_t& context, + f_t snf_b, + f_t lambda); + + flow_cover_evaluation_t evaluate_sgfci(const flow_cover_context_t& context, + f_t snf_b, + f_t lambda); + + bool emit_flow_cover_cut(const flow_cover_context_t& context, + f_t snf_b, + f_t lambda, + const flow_cover_evaluation_t& cmirfci, + const flow_cover_evaluation_t& sgfci, + inequality_t& cut); + std::vector is_slack_; std::vector knapsack_constraints_; + std::vector> flow_cover_constraints_; std::vector is_complemented_; std::vector is_marked_; std::vector workspace_; std::vector complemented_xstar_; + flow_cover_scratch_t flow_cover_scratch_; const simplex_solver_settings_t& settings_; }; @@ -400,9 +587,6 @@ class knapsack_generation_t { template class mixed_integer_rounding_cut_t; -template -class variable_bounds_t; - template class cut_generation_t { public: @@ -470,6 +654,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& var_types, + const std::vector& xstar, + variable_bounds_t& variable_bounds, + 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 b968ad18ea..b8c213f427 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -132,6 +132,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 540e31800b..a0d96b2528 100644 --- a/cpp/src/mip_heuristics/solver.cu +++ b/cpp/src/mip_heuristics/solver.cu @@ -337,6 +337,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 1348d7e7e4..ef0634dcbb 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 @@ -1404,4 +1405,210 @@ TEST(cuts, clique_neos8_phase4_lp_infeasibility_binary_search) EXPECT_EQ(first_infeasible.value(), injected_index); } +// Minimal 0-1 single-node-flow relaxation for the flow-cover separator. +// +// y0 + y1 - y2 <= 4 +// 0 <= y0 <= 3*x0, 0 <= y1 <= 6*x1, 0 <= y2 <= 3*x2 +// +// 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; + + 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(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(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(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(std::span{objective_coefficients}); + + std::vector variable_types = {'I', 'I', 'I', '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; + dual_simplex::lp_problem_t lp; + dual_simplex::csr_matrix_t Arow; + std::vector new_slacks; + std::vector var_types; + + 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( + const mps_parser::mps_data_model_t& model) +{ + flow_cover_test_problem_t test_problem; + 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); + + 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; +} + +std::vector single_node_flow_fractional_solution(int num_cols) +{ + std::vector xstar(num_cols, 0.0); + xstar[0] = 1.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 single_node_flow_y_feasible(const std::vector& y) +{ + const double activity = y[0] + y[1] - y[2]; + return activity <= 4.0 + 1e-8; +} + +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_single_node_flow_cut_valid_at_extreme_points( + const dual_simplex::inequality_t& cut, int num_cols) +{ + 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 < 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 < 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 (!single_node_flow_y_feasible(y)) { continue; } + + std::vector point(num_cols, 0.0); + for (int j = 0; j < 3; j++) { + point[j] = ((x_mask >> j) & 1) != 0 ? 1.0 : 0.0; + point[3 + j] = y[j]; + } + 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 < 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 < 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 < 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 (!single_node_flow_y_feasible(y)) { continue; } + + std::vector point(num_cols, 0.0); + for (int j = 0; j < 3; j++) { + point[j] = ((x_mask >> j) & 1) != 0 ? 1.0 : 0.0; + point[3 + j] = y[j]; + } + expect_single_node_flow_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_generates_valid_single_node_flow_cut) +{ + 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, + 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 auto& 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, + 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.row << " reverse=" << flow_cover_row.reverse; + expect_single_node_flow_cut_valid_at_extreme_points(cut, test_problem.lp.num_cols); + generated_cuts++; + } + + EXPECT_GT(generated_cuts, 0); +} + } // namespace cuopt::linear_programming::test