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

collect_gmi_dual: profile, do not filter at the end

Browse Source
feature/replay^2
Alinson S. Xavier 1 year ago
parent beab75a16d
commit f89903cf68
Signed by: isoron
GPG Key ID: 0DA8E4B9E1109DCA
1 changed files with 105 additions and 175 deletions
Show Stats Download Patch File Download Diff File

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

@ -32,57 +32,43 @@ function collect_gmi_dual(
)
reset_timer!()
# Open HDF5 file
h5_filename = replace(mps_filename, ".mps.gz" => ".h5")
h5 = H5File(h5_filename)
# Read optimal solution
sol_opt_dict = Dict(
zip(
h5.get_array("static_var_names"),
convert(Array{Float64}, h5.get_array("mip_var_values")),
),
)
# Read optimal value
obj_mip = h5.get_scalar("mip_lower_bound")
if obj_mip === nothing
@timeit "Read H5" begin
h5_filename = replace(mps_filename, ".mps.gz" => ".h5")
h5 = H5File(h5_filename)
sol_opt_dict = Dict(
zip(
h5.get_array("static_var_names"),
convert(Array{Float64}, h5.get_array("mip_var_values")),
),
)
obj_mip = h5.get_scalar("mip_obj_value")
h5.file.close()
end
obj_lp = h5.get_scalar("lp_obj_value")
h5.file.close()
# Define relative MIP gap
gap(v) = 100 * abs(obj_mip - v) / abs(v)
# Initialize stats
stats_obj = []
stats_gap = []
stats_ncuts = []
stats_time_convert = 0
stats_time_solve = 0
stats_time_select = 0
stats_time_tableau = 0
stats_time_gmi = 0
stats_time_dual = 0
stats_time_dual_2 = 0
all_cuts = nothing
all_cuts_v2 = nothing
cuts_all = nothing
cuts_all_v2 = nothing
original_basis = nothing
# Read problem
model = read_from_file(mps_filename)
@timeit "Initialize" begin
stats_obj = []
stats_gap = []
stats_ncuts = []
all_cuts = nothing
all_cuts_v2 = nothing
cuts_all = nothing
cuts_all_v2 = nothing
original_basis = nothing
end
# Read original objective function
or_obj_f = objective_function(model)
revised_obj = objective_function(model)
@timeit "Read problem" begin
model = read_from_file(mps_filename)
or_obj_f = objective_function(model)
revised_obj = objective_function(model)
end
for round = 1:max_rounds
@info "Round $(round)..."
stats_time_convert = @elapsed begin
@timeit "Convert to standard form" begin
# Update objective function
set_objective_function(model, revised_obj)
@ -107,51 +93,40 @@ function collect_gmi_dual(
# Assert converted solution is feasible for standard form problem
assert_eq(data_s.constr_lhs * sol_opt_s, data_s.constr_lb)
end
# Optimize standard form
optimize!(model_s)
stats_time_solve += solve_time(model_s)
obj = objective_value(model_s) + data_s.obj_offset
if round == 1
# Assert standard form problem has same value as original
if obj_lp !== nothing
assert_eq(obj, obj_lp)
@timeit "Optimize standard form" begin
optimize!(model_s)
if round == 1
obj = objective_value(model_s) + data_s.obj_offset
push!(stats_obj, obj)
push!(stats_gap, gap(obj))
push!(stats_ncuts, 0)
end
if termination_status(model_s) != MOI.OPTIMAL
error("Non-optimal termination status")
end
obj_lp = obj
push!(stats_obj, obj)
push!(stats_gap, gap(obj))
push!(stats_ncuts, 0)
end
if termination_status(model_s) != MOI.OPTIMAL
return
end
# Store original basis and select tableau rows
basis = get_basis(model_s)
if round == 1
original_basis = basis
end
sol_frac = get_x(model_s)
stats_time_select += @elapsed begin
@timeit "Select tableau rows" begin
basis = get_basis(model_s)
if round == 1
original_basis = basis
end
sol_frac = get_x(model_s)
selected_rows =
select_gmi_rows(data_s, basis, sol_frac, max_rows = max_cuts_per_round)
end
# Compute selected tableau rows
stats_time_tableau += @elapsed begin
@timeit "Compute tableau rows" begin
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)
assert_eq(tableau.lhs * sol_opt_s, tableau.rhs)
end
# Compute GMI cuts
stats_time_gmi += @elapsed begin
@timeit "Compute GMI cuts" begin
cuts_s = compute_gmi(data_s, tableau)
# Assert cuts have been generated correctly
@ -165,39 +140,40 @@ function collect_gmi_dual(
end
end
# Add GMI cuts to original problem
cuts = backwards(transforms, cuts_s)
if round == 1
cuts_all = cuts
basis_vec = repeat([basis], length(selected_rows))
cuts_all_v2 =
ConstraintSet_v2(cuts.lhs, cuts.ub, cuts.lb, basis_vec, selected_rows)
else
# v1 struct
cuts_all.lb = [cuts_all.lb; cuts.lb]
cuts_all.ub = [cuts_all.ub; cuts.ub]
cuts_all.lhs = [cuts_all.lhs; cuts.lhs]
# v2 struct
cuts_all_v2.lb = [cuts_all_v2.lb; cuts.lb]
cuts_all_v2.ub = [cuts_all_v2.ub; cuts.ub]
cuts_all_v2.lhs = [cuts_all_v2.lhs; cuts.lhs]
cuts_all_v2.Bss = [cuts_all_v2.Bss; repeat([basis], length(selected_rows))]
cuts_all_v2.Bv = [cuts_all_v2.Bv; selected_rows]
@timeit "Add GMI cuts to original problem" begin
cuts = backwards(transforms, cuts_s)
if round == 1
cuts_all = cuts
basis_vec = repeat([basis], length(selected_rows))
cuts_all_v2 =
ConstraintSet_v2(cuts.lhs, cuts.ub, cuts.lb, basis_vec, selected_rows)
else
# v1 struct
cuts_all.lb = [cuts_all.lb; cuts.lb]
cuts_all.ub = [cuts_all.ub; cuts.ub]
cuts_all.lhs = [cuts_all.lhs; cuts.lhs]
# v2 struct
cuts_all_v2.lb = [cuts_all_v2.lb; cuts.lb]
cuts_all_v2.ub = [cuts_all_v2.ub; cuts.ub]
cuts_all_v2.lhs = [cuts_all_v2.lhs; cuts.lhs]
cuts_all_v2.Bss = [cuts_all_v2.Bss; repeat([basis], length(selected_rows))]
cuts_all_v2.Bv = [cuts_all_v2.Bv; selected_rows]
end
constrs, gmi_exps = add_constraint_set_dual_v2(model, cuts_all)
end
constrs, gmi_exps = add_constraint_set_dual_v2(model, cuts_all)
# Optimize original form
set_objective_function(model, or_obj_f)
set_optimizer(model, optimizer)
undo_relax = relax_integrality(model)
optimize!(model)
obj = objective_value(model)
push!(stats_obj, obj)
push!(stats_gap, gap(obj))
# Reoptimize with updated obj function
stats_time_dual += @elapsed begin
@timeit "Optimize original form" begin
set_objective_function(model, or_obj_f)
set_optimizer(model, optimizer)
undo_relax = relax_integrality(model)
optimize!(model)
obj = objective_value(model)
push!(stats_obj, obj)
push!(stats_gap, gap(obj))
end
@timeit "Reoptimize with updated obj function" begin
revised_obj = (
or_obj_f - sum(
shadow_price(c) * gmi_exps[iz] for (iz, c) in enumerate(constrs)
@ -213,93 +189,47 @@ function collect_gmi_dual(
undo_relax()
end
# Filter out useless cuts
stats_time_dual_2 += @elapsed begin
set_objective_function(model, or_obj_f)
keep = []
obj_gmi = obj_lp
if (cuts_all !== nothing)
constrs, gmi_exps = add_constraint_set_dual_v2(model, cuts_all)
for (i, c) in enumerate(constrs)
set_name(c, @sprintf("gomory_%05d", i))
@timeit "Store cuts" begin
if all_cuts !== nothing
cut_sizezz = length(all_cuts_v2.Bv)
var_totall =
length(original_basis.var_basic) +
length(original_basis.var_nonbasic) +
length(original_basis.constr_basic) +
length(original_basis.constr_nonbasic)
bm_size = Array{Int64,2}(undef, cut_sizezz, 4)
basis_matrix = Array{Int64,2}(undef, cut_sizezz, var_totall)
for ii = 1:cut_sizezz
vb = all_cuts_v2.Bss[ii].var_basic
vn = all_cuts_v2.Bss[ii].var_nonbasic
cb = all_cuts_v2.Bss[ii].constr_basic
cn = all_cuts_v2.Bss[ii].constr_nonbasic
bm_size[ii, :] = [length(vb) length(vn) length(cb) length(cn)]
basis_matrix[ii, :] = [vb' vn' cb' cn']
end
set_optimizer(model, optimizer)
undo_relax = relax_integrality(model)
optimize!(model)
obj = objective_value(model)
obj_gmi = obj
push!(stats_obj, obj)
push!(stats_gap, gap(obj))
# Store useful cuts; drop useless ones from the problem
useful = [-shadow_price(c) > 1e-3 for c in constrs]
drop = findall(useful .== false)
keep = findall(useful .== true)
all_cuts = ConstraintSet(;
lhs = cuts_all.lhs[keep, :],
lb = cuts_all.lb[keep],
ub = cuts_all.ub[keep],
)
all_cuts_v2 = ConstraintSet_v2(;
lhs = cuts_all_v2.lhs[keep, :],
lb = cuts_all_v2.lb[keep],
ub = cuts_all_v2.ub[keep],
Bss = cuts_all_v2.Bss[keep],
Bv = cuts_all_v2.Bv[keep],
)
delete.(model, constrs[drop])
undo_relax()
@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", basis_matrix)
h5.put_array("cuts_basis_sizes", bm_size)
h5.put_array("cuts_rows", all_cuts_v2.Bv)
h5.file.close()
end
end
basis = original_basis
if all_cuts !== nothing
cut_sizezz = length(all_cuts_v2.Bv)
var_totall =
length(basis.var_basic) +
length(basis.var_nonbasic) +
length(basis.constr_basic) +
length(basis.constr_nonbasic)
bm_size = Array{Int64,2}(undef, cut_sizezz, 4)
basis_matrix = Array{Int64,2}(undef, cut_sizezz, var_totall)
for ii = 1:cut_sizezz
vb = all_cuts_v2.Bss[ii].var_basic
vn = all_cuts_v2.Bss[ii].var_nonbasic
cb = all_cuts_v2.Bss[ii].constr_basic
cn = all_cuts_v2.Bss[ii].constr_nonbasic
bm_size[ii, :] = [length(vb) length(vn) length(cb) length(cn)]
basis_matrix[ii, :] = [vb' vn' cb' cn']
end
# Store cuts
@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", basis_matrix)
h5.put_array("cuts_basis_sizes", bm_size)
h5.put_array("cuts_rows", all_cuts_v2.Bv)
h5.file.close()
end
print_timer()
return OrderedDict(
"instance" => mps_filename,
"max_rounds" => max_rounds,
"rounds" => length(stats_obj) - 1,
"time_convert" => stats_time_convert,
"time_solve" => stats_time_solve,
"time_tableau" => stats_time_tableau,
"time_gmi" => stats_time_gmi,
"time_dual" => stats_time_dual,
"time_dual_2" => stats_time_dual_2,
"obj_mip" => obj_mip,
"obj_lp" => obj_lp,
"stats_obj" => stats_obj,
"stats_gap" => stats_gap,
"stats_ncuts" => length(keep),
)
end

Write Preview
Loading…
Cancel
Save