From 0b5ec4740ed4237683ee963ac59327209926b580 Mon Sep 17 00:00:00 2001 From: "Alinson S. Xavier" Date: Mon, 4 Aug 2025 15:17:17 -0500 Subject: [PATCH] FisSal2011: partial implementation --- src/Cuts/tableau/gmi.jl | 109 +++++++++-- src/Cuts/tableau/gmi_dual.jl | 368 ++++++++++++++++++++++------------- src/Cuts/tableau/structs.jl | 8 +- src/Cuts/tableau/tableau.jl | 1 + 4 files changed, 328 insertions(+), 158 deletions(-) diff --git a/src/Cuts/tableau/gmi.jl b/src/Cuts/tableau/gmi.jl index bdde291..600b912 100644 --- a/src/Cuts/tableau/gmi.jl +++ b/src/Cuts/tableau/gmi.jl @@ -199,23 +199,81 @@ function select_gmi_rows(data, basis, x; max_rows = 10, atol = 0.001) return [candidate_rows[perm[i]] for i = 1:min(length(perm), max_rows)] end +# function compute_gmi(data::ProblemData, tableau::Tableau)::ConstraintSet +# @timeit "Initialization" begin +# nrows, ncols = size(tableau.lhs) +# ub = Float64[Inf for _ = 1:nrows] +# lb = Float64[0.999 for _ = 1:nrows] +# tableau_I, tableau_J, tableau_V = findnz(tableau.lhs) +# lhs_I = Int[] +# lhs_J = Int[] +# lhs_V = Float64[] +# end +# @timeit "Compute coefficients" begin +# for k = 1:nnz(tableau.lhs) +# i::Int = tableau_I[k] +# j::Int = tableau_J[k] +# v::Float64 = 0.0 +# frac_alpha_j = frac(tableau_V[k]) +# alpha_j = tableau_V[k] +# beta = frac(tableau.rhs[i]) +# if data.var_types[j] == 'C' +# if alpha_j >= 0 +# v = alpha_j / beta +# else +# v = -alpha_j / (1 - beta) +# end +# else +# if frac_alpha_j < beta +# v = frac_alpha_j / beta +# else +# v = (1 - frac_alpha_j) / (1 - beta) +# end +# end +# if abs(v) > 1e-8 +# push!(lhs_I, i) +# push!(lhs_J, tableau_J[k]) +# push!(lhs_V, v) +# end +# end +# end +# @timeit "Convert to ConstraintSet" begin +# lhs = sparse(lhs_I, lhs_J, lhs_V, nrows, ncols) +# cs = ConstraintSet(; lhs, ub, lb) +# end +# return cs +# end + function compute_gmi(data::ProblemData, tableau::Tableau)::ConstraintSet - nrows, ncols = size(tableau.lhs) - ub = Float64[Inf for _ = 1:nrows] - lb = Float64[0.999 for _ = 1:nrows] - tableau_I, tableau_J, tableau_V = findnz(tableau.lhs) - lhs_I = Int[] - lhs_J = Int[] - lhs_V = Float64[] + @timeit "Initialization" begin + nrows::Int, ncols::Int = size(tableau.lhs) + var_types::Vector{Char} = data.var_types + tableau_rhs::Vector{Float64} = tableau.rhs + tableau_I::Vector{Int}, tableau_J::Vector{Int}, tableau_V::Vector{Float64} = findnz(tableau.lhs) + end + + @timeit "Pre-allocation" begin + ub::Vector{Float64} = fill(Inf, nrows) + lb::Vector{Float64} = fill(0.999, nrows) + nnz_tableau::Int = length(tableau_I) + lhs_I::Vector{Int} = Vector{Int}(undef, nnz_tableau) + lhs_J::Vector{Int} = Vector{Int}(undef, nnz_tableau) + lhs_V::Vector{Float64} = Vector{Float64}(undef, nnz_tableau) + nnz_count::Int = 0 + end + @timeit "Compute coefficients" begin - for k = 1:nnz(tableau.lhs) + @inbounds for k = 1:nnz_tableau i::Int = tableau_I[k] j::Int = tableau_J[k] - v::Float64 = 0.0 - frac_alpha_j = frac(tableau_V[k]) - alpha_j = tableau_V[k] - beta = frac(tableau.rhs[i]) - if data.var_types[j] == 'C' + alpha_j::Float64 = tableau_V[k] + frac_alpha_j::Float64 = alpha_j - floor(alpha_j) + beta_i::Float64 = tableau_rhs[i] + beta::Float64 = beta_i - floor(beta_i) + v::Float64 = 0 + + # Compute coefficient + if var_types[j] == 'C' if alpha_j >= 0 v = alpha_j / beta else @@ -228,16 +286,31 @@ function compute_gmi(data::ProblemData, tableau::Tableau)::ConstraintSet v = (1 - frac_alpha_j) / (1 - beta) end end + + # Store if significant if abs(v) > 1e-8 - push!(lhs_I, i) - push!(lhs_J, tableau_J[k]) - push!(lhs_V, v) + nnz_count += 1 + lhs_I[nnz_count] = i + lhs_J[nnz_count] = j + lhs_V[nnz_count] = v end end - lhs = sparse(lhs_I, lhs_J, lhs_V, nrows, ncols) end - return ConstraintSet(; lhs, ub, lb) + + @timeit "Resize arrays to actual size" begin + resize!(lhs_I, nnz_count) + resize!(lhs_J, nnz_count) + resize!(lhs_V, nnz_count) + end + + @timeit "Convert to ConstraintSet" begin + lhs::SparseMatrixCSC = sparse(lhs_I, lhs_J, lhs_V, nrows, ncols) + cs::ConstraintSet = ConstraintSet(; lhs, ub, lb) + end + + return cs end + export compute_gmi, frac, select_gmi_rows, assert_cuts_off, assert_does_not_cut_off, collect_gmi diff --git a/src/Cuts/tableau/gmi_dual.jl b/src/Cuts/tableau/gmi_dual.jl index d6535eb..06fe13c 100644 --- a/src/Cuts/tableau/gmi_dual.jl +++ b/src/Cuts/tableau/gmi_dual.jl @@ -269,9 +269,10 @@ function collect_gmi_FisSal2011( optimizer, max_rounds = 10_000, max_cuts_per_round = 1_000_000, - time_limit = 30, + time_limit = 1_000_000, interval_print=1, - max_cut_age=10, + silent_solver=true, + max_pool_size_mb = 1024, ) reset_timer!() initial_time = time() @@ -290,26 +291,37 @@ function collect_gmi_FisSal2011( end @timeit "Initialize" begin - stats_obj = [] - stats_gap = [] - stats_ncuts = [] - pool = nothing - cut_age = nothing - lambda_curr = nothing - lambda_best = nothing - last_print_time = 0 - obj_initial = nothing - obj_curr = 0 - obj_best = 0 - noimprov_count = 0 backtrack_count = 0 + deterioration_count = 0 + # ε = 0.01 + basis_curr = nothing + basis_prev = nothing + basis_seen = Set{UInt64}() gapcl_best = 0 + gapcl_curr = 0 + multipliers_curr = nothing + multipliers_best = nothing + last_print_time = 0 + obj_curr = 0 + obj_initial = nothing + obj_best = nothing + obj_hist = [] + pool = nothing + pool_active = 0 + pool_cut_age = nothing + pool_cut_hashes = Set{UInt64}() + pool_size_mb = 0 + stats_gap = [] + stats_ncuts = [] + stats_obj = [] + μ = 10 + λ = 0 end gap(v) = 100 * abs(obj_mip - v) / abs(obj_mip) gapcl(v) = 100 * (v - obj_initial) / (obj_mip - obj_initial) - function perturb(v, ε=0.001) + function perturb(v, ε) p = (1 - ε) .+ 2 * ε * rand(length(v)) return v .* p end @@ -333,6 +345,9 @@ function collect_gmi_FisSal2011( # Convert to standard form data_s, transforms = convert_to_standard_form(data) model_s = to_model(data_s) + if silent_solver + set_silent(model_s) + end vars_s = all_variables(model_s) orig_obj_s = objective_function(model_s) set_optimizer(model_s, optimizer) @@ -345,13 +360,16 @@ function collect_gmi_FisSal2011( assert_eq(data_s.constr_lhs * sol_opt_s, data_s.constr_lb) end + @info "Standard form model has $(length(data.var_lb)) vars, $(length(data.constr_lb)) constrs" + @show obj_mip + for round = 1:max_rounds if round > 1 @timeit "Update objective function" begin # Build Lagrangian term - lambda_perturbed = perturb(lambda_curr) - v = sparse(pool.lhs' * lambda_perturbed) - lagr_term = AffExpr(dot(lambda_perturbed, pool.lb)) + # lambda_perturbed = perturb(multipliers, 0.1) + v = sparse(pool.lhs * multipliers_curr) + lagr_term = AffExpr(dot(multipliers_curr, pool.lb)) for offset in 1:nnz(v) var_idx = v.nzind[offset] add_to_expression!( @@ -370,32 +388,25 @@ function collect_gmi_FisSal2011( end @timeit "Optimize LP (lagrangian)" begin - set_silent(model_s) optimize!(model_s) sol_frac = get_x(model_s) obj_curr = objective_value(model_s) obj_curr <= obj_mip || error("LP value higher than MIP value: $(obj_curr) > $(obj_mip)") - - if round == 1 - obj_initial = obj_curr + + push!(obj_hist, obj_curr) + if length(obj_hist) > 100 + popfirst!(obj_hist) end - if obj_curr >= obj_best + if obj_best === nothing || obj_curr > obj_best obj_best = obj_curr - gapcl_best = gapcl(obj_best) - lambda_best = lambda_curr - noimprov_count = 0 - else - noimprov_count += 1 + multipliers_best = multipliers_curr end - - if noimprov_count > 10 - lambda_curr = lambda_best - backtrack_count += 1 - noimprov_count = 0 - continue + if round == 1 + obj_initial = obj_curr end - + gapcl_curr = gapcl(obj_curr) + gapcl_best = gapcl(obj_best) push!(stats_obj, obj_curr) push!(stats_gap, gap(obj_curr)) if round == 1 @@ -406,121 +417,191 @@ function collect_gmi_FisSal2011( if termination_status(model_s) != MOI.OPTIMAL error("Non-optimal termination status") end - end - @timeit "Select tableau rows" begin - basis = get_basis(model_s) - if round == 1 - original_basis = basis + Δ = obj_mip - obj_best + if obj_curr < obj_best - Δ + deterioration_count += 1 + else + deterioration_count = 0 + end + if deterioration_count >= 10 + μ *= 0.5 + multipliers_curr = multipliers_best + deterioration_count = 0 + backtrack_count += 1 + elseif length(obj_hist) >= 100 + obj_hist_avg = mean(obj_hist) + improv = obj_best - obj_hist[1] + if improv < 0.01 * Δ + if obj_best - obj_hist_avg < 0.001 * Δ + μ = 10 * μ + elseif obj_best - obj_hist_avg < 0.01 * Δ + μ = 2 * μ + else + μ = 0.5 * μ + end + end end - selected_rows = - select_gmi_rows(data_s, basis, sol_frac, max_rows = max_cuts_per_round) end - @timeit "Compute tableau rows" begin - tableau = compute_tableau(data_s, basis, x = sol_frac, rows = selected_rows) + if mod(round, 10) == 1 + @timeit "Select tableau rows" begin + basis_prev = basis_curr + basis_curr = get_basis(model_s) + basis_hash = hash(basis_curr) + push!(basis_seen, basis_hash) + selected_rows = + select_gmi_rows(data_s, basis_curr, sol_frac, max_rows = max_cuts_per_round) + end - # Assert tableau rows have been computed correctly - assert_eq(tableau.lhs * sol_frac, tableau.rhs, atol=1e-3) - assert_eq(tableau.lhs * sol_opt_s, tableau.rhs, atol=1e-3) - end + @timeit "Compute tableau rows" begin + tableau = compute_tableau(data_s, basis_curr, x = sol_frac, rows = selected_rows) - @timeit "Compute GMI cuts" begin - cuts_s = compute_gmi(data_s, tableau) - assert_cuts_off(cuts_s, sol_frac) - assert_does_not_cut_off(cuts_s, sol_opt_s) - ncuts = length(cuts_s.lb) - end + # Assert tableau rows have been computed correctly + assert_eq(tableau.lhs * sol_frac, tableau.rhs, atol=1e-3) + assert_eq(tableau.lhs * sol_opt_s, tableau.rhs, atol=1e-3) + end - @timeit "Add new cuts to the pool" begin - if round == 1 - pool = cuts_s - lambda_curr = zeros(ncuts) - lambda_best = zeros(ncuts) - cut_age = zeros(ncuts) - else - pool.lhs = [pool.lhs; cuts_s.lhs] - pool.lb = [pool.lb; cuts_s.lb] - pool.ub = [pool.ub; cuts_s.ub] - lambda_curr = [lambda_curr; zeros(ncuts)] - lambda_best = [lambda_best; zeros(ncuts)] - cut_age = [cut_age; zeros(ncuts)] + @timeit "Compute GMI cuts" begin + cuts_s = compute_gmi(data_s, tableau) + assert_cuts_off(cuts_s, sol_frac) + assert_does_not_cut_off(cuts_s, sol_opt_s) + ncuts = length(cuts_s.lb) + end + + @timeit "Add new cuts to the pool" begin + @timeit "Compute cut hashses" begin + unique_indices = Int[] + for i in 1:ncuts + cut_data = (cuts_s.lhs[i, :], cuts_s.lb[i], cuts_s.ub[i]) + cut_hash = hash(cut_data) + if !(cut_hash in pool_cut_hashes) + push!(pool_cut_hashes, cut_hash) + push!(unique_indices, i) + end + end + end + @timeit "Append unique cuts" begin + if round == 1 + pool = ConstraintSet( + lhs = sparse(cuts_s.lhs[unique_indices, :]'), + lb = cuts_s.lb[unique_indices], + ub = cuts_s.ub[unique_indices] + ) + ncuts_unique = length(unique_indices) + multipliers_curr = zeros(ncuts_unique) + multipliers_best = zeros(ncuts_unique) + pool_cut_age = zeros(ncuts_unique) + else + if !isempty(unique_indices) + ncuts_unique = length(unique_indices) + pool.lhs = [pool.lhs sparse(cuts_s.lhs[unique_indices, :]')] + pool.lb = [pool.lb; cuts_s.lb[unique_indices]] + pool.ub = [pool.ub; cuts_s.ub[unique_indices]] + multipliers_curr = [multipliers_curr; zeros(ncuts_unique)] + multipliers_best = [multipliers_best; zeros(ncuts_unique)] + pool_cut_age = [pool_cut_age; zeros(ncuts_unique)] + end + end + end end end - # @timeit "Update multipliers (subgradient)" begin - # subgrad = pool.lb .- pool.lhs * sol_frac - # lambda = max.(0, lambda .+ 0.01 * subgrad) - # end - selected_idx = [] selected_contrs = [] + if round == 1 || round == max_rounds + @timeit "Update multipliers (large LP)" begin + while true + @timeit "Optimize LP (extended)" begin + set_objective_function(model_s, orig_obj_s) + optimize!(model_s) + sol_frac = get_x(model_s) + end - @timeit "Update multipliers (large LP)" begin - while true - @timeit "Optimize LP (extended)" begin - set_objective_function(model_s, orig_obj_s) - optimize!(model_s) - sol_frac = get_x(model_s) - end + @timeit "Computing cut violations" begin + violations = (pool.lb' - (sol_frac' * pool.lhs))' + end - @timeit "Find most violated cut" begin - violations = pool.lb .- pool.lhs * sol_frac - σ = sortperm(violations, rev=true) - end + @timeit "Sorting cut violations" begin + σ = sortperm(violations, rev=true) + end - # Stop if all cuts are satisfied - if violations[σ[1]] <= 1e-6 - break - end + if violations[σ[1]] <= 1e-6 + break + end - @timeit "Add constraint to the model" begin - push!(selected_idx, σ[1]) - cut_lhs = pool.lhs[σ[1], :] - cut_lhs_value = 0.0 - cut_lb = pool.lb[σ[1]] - cut_expr = AffExpr() - for offset in 1:nnz(cut_lhs) - var_idx = cut_lhs.nzind[offset] - add_to_expression!( - cut_expr, - vars_s[var_idx], - cut_lhs.nzval[offset], - ) - cut_lhs_value += sol_frac[var_idx] * cut_lhs.nzval[offset] + @timeit "Add constraint to the model" begin + push!(selected_idx, σ[1]) + cut_lhs = pool.lhs[:, σ[1]] + cut_lhs_value = 0.0 + cut_lb = pool.lb[σ[1]] + cut_expr = AffExpr() + for offset in 1:nnz(cut_lhs) + var_idx = cut_lhs.nzind[offset] + add_to_expression!( + cut_expr, + vars_s[var_idx], + cut_lhs.nzval[offset], + ) + cut_lhs_value += sol_frac[var_idx] * cut_lhs.nzval[offset] + end + cut_constr = @constraint(model_s, cut_expr >= cut_lb) + push!(selected_contrs, cut_constr) end - cut_constr = @constraint(model_s, cut_expr >= cut_lb) - push!(selected_contrs, cut_constr) end - end - @timeit "Find dual values for all selected cuts" begin - lambda_curr .= 0 - cut_age .+= 1 - for (offset, idx) in enumerate(selected_idx) - lambda_curr[idx] = -shadow_price(selected_contrs[offset]) - if lambda_curr[idx] > 1e-5 - cut_age[idx] = 0 + pool_active = length(selected_idx) + + @timeit "Find dual values for all selected cuts" begin + multipliers_curr .= 0 + pool_cut_age .+= 1 + for (offset, idx) in enumerate(selected_idx) + multipliers_curr[idx] = -shadow_price(selected_contrs[offset]) + if multipliers_curr[idx] > 1e-5 + pool_cut_age[idx] = 0 + end end end - end - # Filter cut pool - keep = findall(cut_age .< max_cut_age) - pool.lhs = pool.lhs[keep, :] - pool.ub = pool.ub[keep] - pool.lb = pool.lb[keep] - lambda_curr = lambda_curr[keep] - lambda_best = lambda_best[keep] - cut_age = cut_age[keep] + @timeit "Prune cut pool" begin + pool_size_mb = Base.summarysize(pool) / 1024^2 + while pool_size_mb >= max_pool_size_mb + @timeit "Identify cuts to remove" begin + σ = sortperm(pool_cut_age, rev=true) + pool_size = length(pool.ub) + n_keep = Int(floor(pool_size * 0.8)) + idx_keep = 1:n_keep + idx_remove = σ[(n_keep+1):end] + end + @timeit "Update cut hashes" begin + for idx in idx_remove + cut_data = (pool.lhs[:, idx], pool.lb[idx], pool.ub[idx]) + delete!(pool_cut_hashes, hash(cut_data)) + end + end + @timeit "Update cut pool" begin + pool.lhs = pool.lhs[:, idx_keep] + pool.ub = pool.ub[idx_keep] + pool.lb = pool.lb[idx_keep] + multipliers_curr = multipliers[idx_keep] + pool_cut_age = pool_cut_age[idx_keep] + end + pool_size_mb = Base.summarysize(pool) / 1024^2 + end + end - @timeit "Delete all cut constraints" begin - delete.(model_s, selected_contrs) + @timeit "Delete all cut constraints" begin + delete.(model_s, selected_contrs) + end + end + else + @timeit "Update multipliers (subgradient)" begin + subgrad = (pool.lb' - (sol_frac' * pool.lhs))' + λ = μ * (obj_mip - obj_curr) / norm(subgrad)^2 + multipliers_curr = max.(0, multipliers_curr .+ λ * subgrad) end end - push!(stats_ncuts, length(pool.lb)) - elapsed_time = time() - initial_time if elapsed_time > time_limit @info "Time limit exceeded. Stopping." @@ -529,28 +610,38 @@ function collect_gmi_FisSal2011( if round == 1 @printf( - "%8s %12s %12s %12s %8s %8s %9s\n", + "%8s %9s %7s %7s %8s %9s %9s %9s %4s %8s %8s %8s\n", "round", "obj", - "gapcl_curr", - "gapcl_best", + "cl_curr", + "cl_best", "active", - "pool", - "backtrack", + "pool_cuts", + "pool_mb", + "bases", + "bktk", + "Δ", + "μ", + "λ", ) end if time() - last_print_time > interval_print last_print_time = time() @printf( - "%8d %12.6e %12.2f %12.2f %8d %8d %9d\n", + "%8d %9.3e %7.2f %7.2f %8d %9d %9.2f %9d %4d %8.2e %8.2e %8.2e\n", round, obj_curr, - gapcl(obj_curr), + gapcl_curr, gapcl_best, length(selected_idx), length(pool.ub), + pool_size_mb, + length(basis_seen), backtrack_count, + Δ, + μ, + λ, ) end end @@ -585,18 +676,23 @@ function collect_gmi_FisSal2011( # end # end + @info "Best gap closure: $(gapcl_best)" + to = TimerOutputs.get_defaulttimer() stats_time = TimerOutputs.tottime(to) / 1e9 print_timer() return OrderedDict( + "gapcl_best" => gapcl_best, + "gapcl_curr" => gapcl_curr, "instance" => mps_filename, - "max_rounds" => max_rounds, - "rounds" => length(stats_obj) - 1, + "obj_final" => obj_curr, + "obj_initial" => obj_initial, "obj_mip" => obj_mip, - "stats_obj" => stats_obj, - "stats_ncuts" => stats_ncuts, - "stats_time" => stats_time, + "pool_active" => pool_active, + "pool_size_mb" => pool_size_mb, + "pool_total" => length(pool.lb), + "time" => stats_time, ) end diff --git a/src/Cuts/tableau/structs.jl b/src/Cuts/tableau/structs.jl index 7a2bb07..b6bca2e 100644 --- a/src/Cuts/tableau/structs.jl +++ b/src/Cuts/tableau/structs.jl @@ -18,10 +18,10 @@ Base.@kwdef mutable struct ProblemData end Base.@kwdef mutable struct Tableau - obj::Any - lhs::Any - rhs::Any - z::Any + obj::Vector{Float64} + lhs::SparseMatrixCSC + rhs::Vector{Float64} + z::Float64 end Base.@kwdef mutable struct Basis diff --git a/src/Cuts/tableau/tableau.jl b/src/Cuts/tableau/tableau.jl index 06bf3bd..270ec16 100644 --- a/src/Cuts/tableau/tableau.jl +++ b/src/Cuts/tableau/tableau.jl @@ -102,6 +102,7 @@ function compute_tableau( sol = factor \ obj_b tableau_obj = -data.obj' + sol' * data.constr_lhs tableau_obj[abs.(tableau_obj).