ANL-CEEESA
/
MIPLearn.jl
mirror of https://github.com/ANL-CEEESA/MIPLearn.jl.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

FisSal2011: clean up, improve gap closure on MIPLIB 3 (65.5%)

Browse Source
fs11_01
Alinson S. Xavier 2 months ago
parent 0b5ec4740e
commit 0a0d133161
1 changed files with 41 additions and 92 deletions
Show Stats Download Patch File Download Diff File

133
src/Cuts/tableau/gmi_dual.jl
Unescape View File

@ -269,8 +269,8 @@ function collect_gmi_FisSal2011(
optimizer, optimizer,
max_rounds = 10_000, max_rounds = 10_000,
max_cuts_per_round = 1_000_000, max_cuts_per_round = 1_000_000,
time_limit = 1_000_000, time_limit = 900,
interval_print=1, print_interval_sec=1,
silent_solver=true, silent_solver=true,
max_pool_size_mb = 1024, max_pool_size_mb = 1024,
) )
@ -291,41 +291,28 @@ function collect_gmi_FisSal2011(
end end
@timeit "Initialize" begin @timeit "Initialize" begin
backtrack_count = 0
deterioration_count = 0
# ε = 0.01
basis_curr = nothing
basis_prev = nothing
basis_seen = Set{UInt64}() basis_seen = Set{UInt64}()
count_backtrack = 0
count_deterioration = 0
gapcl_best = 0 gapcl_best = 0
gapcl_curr = 0 gapcl_curr = 0
multipliers_curr = nothing
multipliers_best = nothing
last_print_time = 0 last_print_time = 0
obj_curr = 0 multipliers_best = nothing
obj_initial = nothing multipliers_curr = nothing
obj_best = nothing obj_best = nothing
obj_curr = nothing
obj_hist = [] obj_hist = []
obj_initial = nothing
pool = nothing pool = nothing
pool_active = 0
pool_cut_age = nothing pool_cut_age = nothing
pool_cut_hashes = Set{UInt64}() pool_cut_hashes = Set{UInt64}()
pool_size_mb = 0 pool_size_mb = 0
stats_gap = []
stats_ncuts = []
stats_obj = []
μ = 10
λ = 0 λ = 0
μ = 10
end end
gap(v) = 100 * abs(obj_mip - v) / abs(obj_mip)
gapcl(v) = 100 * (v - obj_initial) / (obj_mip - obj_initial) gapcl(v) = 100 * (v - obj_initial) / (obj_mip - obj_initial)
function perturb(v, ε)
p = (1 - ε) .+ 2 * ε * rand(length(v))
return v .* p
end
@timeit "Read problem" begin @timeit "Read problem" begin
model = read_from_file(mps_filename) model = read_from_file(mps_filename)
set_optimizer(model, optimizer) set_optimizer(model, optimizer)
@ -361,13 +348,11 @@ function collect_gmi_FisSal2011(
end end
@info "Standard form model has $(length(data.var_lb)) vars, $(length(data.constr_lb)) constrs" @info "Standard form model has $(length(data.var_lb)) vars, $(length(data.constr_lb)) constrs"
@show obj_mip
for round = 1:max_rounds for round = 1:max_rounds
if round > 1 if round > 1
@timeit "Update objective function" begin @timeit "Update objective function" begin
# Build Lagrangian term # Build Lagrangian term
# lambda_perturbed = perturb(multipliers, 0.1)
v = sparse(pool.lhs * multipliers_curr) v = sparse(pool.lhs * multipliers_curr)
lagr_term = AffExpr(dot(multipliers_curr, pool.lb)) lagr_term = AffExpr(dot(multipliers_curr, pool.lb))
for offset in 1:nnz(v) for offset in 1:nnz(v)
@ -391,7 +376,6 @@ function collect_gmi_FisSal2011(
optimize!(model_s) optimize!(model_s)
sol_frac = get_x(model_s) sol_frac = get_x(model_s)
obj_curr = objective_value(model_s) obj_curr = objective_value(model_s)
obj_curr <= obj_mip || error("LP value higher than MIP value: $(obj_curr) > $(obj_mip)")
push!(obj_hist, obj_curr) push!(obj_hist, obj_curr)
if length(obj_hist) > 100 if length(obj_hist) > 100
@ -407,28 +391,21 @@ function collect_gmi_FisSal2011(
end end
gapcl_curr = gapcl(obj_curr) gapcl_curr = gapcl(obj_curr)
gapcl_best = gapcl(obj_best) gapcl_best = gapcl(obj_best)
push!(stats_obj, obj_curr)
push!(stats_gap, gap(obj_curr))
if round == 1
push!(stats_ncuts, 0)
else
push!(stats_ncuts, length(pool.lb))
end
if termination_status(model_s) != MOI.OPTIMAL if termination_status(model_s) != MOI.OPTIMAL
error("Non-optimal termination status") error("Non-optimal termination status")
end end
Δ = obj_mip - obj_best Δ = obj_mip - obj_best
if obj_curr < obj_best - Δ if obj_curr < obj_best - Δ
deterioration_count += 1 count_deterioration += 1
else else
deterioration_count = 0 count_deterioration = 0
end end
if deterioration_count >= 10 if count_deterioration >= 10
μ *= 0.5 μ *= 0.5
multipliers_curr = multipliers_best multipliers_curr = multipliers_best
deterioration_count = 0 count_deterioration = 0
backtrack_count += 1 count_backtrack += 1
elseif length(obj_hist) >= 100 elseif length(obj_hist) >= 100
obj_hist_avg = mean(obj_hist) obj_hist_avg = mean(obj_hist)
improv = obj_best - obj_hist[1] improv = obj_best - obj_hist[1]
@ -446,18 +423,15 @@ function collect_gmi_FisSal2011(
if mod(round, 10) == 1 if mod(round, 10) == 1
@timeit "Select tableau rows" begin @timeit "Select tableau rows" begin
basis_prev = basis_curr basis = get_basis(model_s)
basis_curr = get_basis(model_s) basis_hash = hash(basis)
basis_hash = hash(basis_curr)
push!(basis_seen, basis_hash) push!(basis_seen, basis_hash)
selected_rows = selected_rows =
select_gmi_rows(data_s, basis_curr, sol_frac, max_rows = max_cuts_per_round) select_gmi_rows(data_s, basis, sol_frac, max_rows = max_cuts_per_round)
end end
@timeit "Compute tableau rows" begin @timeit "Compute tableau rows" begin
tableau = compute_tableau(data_s, basis_curr, x = sol_frac, rows = selected_rows) tableau = compute_tableau(data_s, basis, x = sol_frac, rows = selected_rows)
# Assert tableau rows have been computed correctly
assert_eq(tableau.lhs * sol_frac, tableau.rhs, atol=1e-3) assert_eq(tableau.lhs * sol_frac, tableau.rhs, atol=1e-3)
assert_eq(tableau.lhs * sol_opt_s, tableau.rhs, atol=1e-3) assert_eq(tableau.lhs * sol_opt_s, tableau.rhs, atol=1e-3)
end end
@ -507,14 +481,15 @@ function collect_gmi_FisSal2011(
end end
end end
selected_idx = [] if mod(round, 1000) == 1 || round == max_rounds
selected_contrs = []
if round == 1 || round == max_rounds
@timeit "Update multipliers (large LP)" begin @timeit "Update multipliers (large LP)" begin
selected_idx = []
selected_contrs = []
while true while true
@timeit "Optimize LP (extended)" begin @timeit "Optimize LP (extended)" begin
set_objective_function(model_s, orig_obj_s) set_objective_function(model_s, orig_obj_s)
optimize!(model_s) optimize!(model_s)
obj_curr = objective_value(model_s)
sol_frac = get_x(model_s) sol_frac = get_x(model_s)
end end
@ -550,8 +525,6 @@ function collect_gmi_FisSal2011(
end end
end end
pool_active = length(selected_idx)
@timeit "Find dual values for all selected cuts" begin @timeit "Find dual values for all selected cuts" begin
multipliers_curr .= 0 multipliers_curr .= 0
pool_cut_age .+= 1 pool_cut_age .+= 1
@ -563,6 +536,15 @@ function collect_gmi_FisSal2011(
end end
end end
@timeit "Update best" begin
if obj_curr > obj_best
obj_best = obj_curr
multipliers_best = multipliers_curr
end
gapcl_curr = gapcl(obj_curr)
gapcl_best = gapcl(obj_best)
end
@timeit "Prune cut pool" begin @timeit "Prune cut pool" begin
pool_size_mb = Base.summarysize(pool) / 1024^2 pool_size_mb = Base.summarysize(pool) / 1024^2
while pool_size_mb >= max_pool_size_mb while pool_size_mb >= max_pool_size_mb
@ -602,20 +584,13 @@ function collect_gmi_FisSal2011(
end end
end end
elapsed_time = time() - initial_time
if elapsed_time > time_limit
@info "Time limit exceeded. Stopping."
break
end
if round == 1 if round == 1
@printf( @printf(
"%8s %9s %7s %7s %8s %9s %9s %9s %4s %8s %8s %8s\n", "%8s %10s %9s %9s %9s %9s %9s %4s %8s %8s %8s\n",
"round", "round",
"obj", "obj",
"cl_curr", "cl_curr",
"cl_best", "cl_best",
"active",
"pool_cuts", "pool_cuts",
"pool_mb", "pool_mb",
"bases", "bases",
@ -626,55 +601,30 @@ function collect_gmi_FisSal2011(
) )
end end
if time() - last_print_time > interval_print if time() - last_print_time > print_interval_sec
last_print_time = time() last_print_time = time()
@printf( @printf(
"%8d %9.3e %7.2f %7.2f %8d %9d %9.2f %9d %4d %8.2e %8.2e %8.2e\n", "%8d %10.3e %9.2e %9.2e %9d %9.2f %9d %4d %8.2e %8.2e %8.2e\n",
round, round,
obj_curr, obj_curr,
gapcl_curr, gapcl_curr,
gapcl_best, gapcl_best,
length(selected_idx),
length(pool.ub), length(pool.ub),
pool_size_mb, pool_size_mb,
length(basis_seen), length(basis_seen),
backtrack_count, count_backtrack,
Δ, Δ,
μ, μ,
λ, λ,
) )
end end
end
# @timeit "Store cuts in H5 file" begin elapsed_time = time() - initial_time
# if all_cuts !== nothing if elapsed_time > time_limit
# ncuts = length(all_cuts_rows) @info "Time limit exceeded. Stopping."
# total = break
# length(original_basis.var_basic) + end
# length(original_basis.var_nonbasic) + end
# length(original_basis.constr_basic) +
# length(original_basis.constr_nonbasic)
# all_cuts_basis_sizes = Array{Int64,2}(undef, ncuts, 4)
# all_cuts_basis_vars = Array{Int64,2}(undef, ncuts, total)
# for i = 1:ncuts
# vb = all_cuts_bases[i].var_basic
# vn = all_cuts_bases[i].var_nonbasic
# cb = all_cuts_bases[i].constr_basic
# cn = all_cuts_bases[i].constr_nonbasic
# all_cuts_basis_sizes[i, :] = [length(vb) length(vn) length(cb) length(cn)]
# all_cuts_basis_vars[i, :] = [vb' vn' cb' cn']
# end
# @info "Storing $(length(all_cuts.ub)) GMI cuts..."
# h5 = H5File(h5_filename)
# h5.put_sparse("cuts_lhs", all_cuts.lhs)
# h5.put_array("cuts_lb", all_cuts.lb)
# h5.put_array("cuts_ub", all_cuts.ub)
# h5.put_array("cuts_basis_vars", all_cuts_basis_vars)
# h5.put_array("cuts_basis_sizes", all_cuts_basis_sizes)
# h5.put_array("cuts_rows", all_cuts_rows)
# h5.file.close()
# end
# end
@info "Best gap closure: $(gapcl_best)" @info "Best gap closure: $(gapcl_best)"
@ -689,7 +639,6 @@ function collect_gmi_FisSal2011(
"obj_final" => obj_curr, "obj_final" => obj_curr,
"obj_initial" => obj_initial, "obj_initial" => obj_initial,
"obj_mip" => obj_mip, "obj_mip" => obj_mip,
"pool_active" => pool_active,
"pool_size_mb" => pool_size_mb, "pool_size_mb" => pool_size_mb,
"pool_total" => length(pool.lb), "pool_total" => length(pool.lb),
"time" => stats_time, "time" => stats_time,

Write Preview
Loading…
Cancel
Save