commit
cbedb02a9f
@ -1,184 +0,0 @@
|
|||||||
# MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
|
|
||||||
# Copyright (C) 2020-2023, UChicago Argonne, LLC. All rights reserved.
|
|
||||||
# Released under the modified BSD license. See COPYING.md for more details.
|
|
||||||
|
|
||||||
import ..H5File
|
|
||||||
|
|
||||||
using OrderedCollections
|
|
||||||
|
|
||||||
function collect_gmi(mps_filename; optimizer, max_rounds = 10, max_cuts_per_round = 100)
|
|
||||||
@info mps_filename
|
|
||||||
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
|
|
||||||
obj_mip = h5.get_scalar("mip_obj_value")
|
|
||||||
end
|
|
||||||
obj_lp = nothing
|
|
||||||
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
|
|
||||||
all_cuts = nothing
|
|
||||||
|
|
||||||
# Read problem
|
|
||||||
model = read_from_file(mps_filename)
|
|
||||||
|
|
||||||
for round = 1:max_rounds
|
|
||||||
@info "Round $(round)..."
|
|
||||||
|
|
||||||
stats_time_convert = @elapsed begin
|
|
||||||
# Extract problem data
|
|
||||||
data = ProblemData(model)
|
|
||||||
|
|
||||||
# Construct optimal solution vector (with correct variable sequence)
|
|
||||||
sol_opt = [sol_opt_dict[n] for n in data.var_names]
|
|
||||||
|
|
||||||
# Assert optimal solution is feasible for the original problem
|
|
||||||
@assert all(data.constr_lb .- 1e-3 .<= data.constr_lhs * sol_opt)
|
|
||||||
@assert all(data.constr_lhs * sol_opt .<= data.constr_ub .+ 1e-3)
|
|
||||||
|
|
||||||
# Convert to standard form
|
|
||||||
data_s, transforms = convert_to_standard_form(data)
|
|
||||||
model_s = to_model(data_s)
|
|
||||||
set_optimizer(model_s, optimizer)
|
|
||||||
relax_integrality(model_s)
|
|
||||||
|
|
||||||
# Convert optimal solution to standard form
|
|
||||||
sol_opt_s = forward(transforms, sol_opt)
|
|
||||||
|
|
||||||
# Assert converted solution is feasible for standard form problem
|
|
||||||
@assert 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 obj_lp === nothing
|
|
||||||
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
|
|
||||||
|
|
||||||
# Select tableau rows
|
|
||||||
basis = get_basis(model_s)
|
|
||||||
sol_frac = get_x(model_s)
|
|
||||||
stats_time_select += @elapsed begin
|
|
||||||
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
|
|
||||||
tableau = compute_tableau(data_s, basis, sol_frac, rows = selected_rows)
|
|
||||||
|
|
||||||
# Assert tableau rows have been computed correctly
|
|
||||||
@assert tableau.lhs * sol_frac ≈ tableau.rhs
|
|
||||||
@assert tableau.lhs * sol_opt_s ≈ tableau.rhs
|
|
||||||
end
|
|
||||||
|
|
||||||
# Compute GMI cuts
|
|
||||||
stats_time_gmi += @elapsed begin
|
|
||||||
cuts_s = compute_gmi(data_s, tableau)
|
|
||||||
|
|
||||||
# Assert cuts have been generated correctly
|
|
||||||
try
|
|
||||||
assert_cuts_off(cuts_s, sol_frac)
|
|
||||||
assert_does_not_cut_off(cuts_s, sol_opt_s)
|
|
||||||
catch
|
|
||||||
@warn "Invalid cuts detected. Discarding round $round cuts and aborting."
|
|
||||||
break
|
|
||||||
end
|
|
||||||
|
|
||||||
# Abort if no cuts are left
|
|
||||||
if length(cuts_s.lb) == 0
|
|
||||||
@info "No cuts generated. Aborting."
|
|
||||||
break
|
|
||||||
end
|
|
||||||
end
|
|
||||||
|
|
||||||
# Add GMI cuts to original problem
|
|
||||||
cuts = backwards(transforms, cuts_s)
|
|
||||||
assert_does_not_cut_off(cuts, sol_opt)
|
|
||||||
constrs = add_constraint_set(model, cuts)
|
|
||||||
|
|
||||||
# Optimize original form
|
|
||||||
set_optimizer(model, optimizer)
|
|
||||||
undo_relax = relax_integrality(model)
|
|
||||||
optimize!(model)
|
|
||||||
obj = objective_value(model)
|
|
||||||
push!(stats_obj, obj)
|
|
||||||
push!(stats_gap, gap(obj))
|
|
||||||
|
|
||||||
# Store useful cuts; drop useless ones from the problem
|
|
||||||
useful = [abs(shadow_price(c)) > 1e-3 for c in constrs]
|
|
||||||
drop = findall(useful .== false)
|
|
||||||
keep = findall(useful .== true)
|
|
||||||
delete.(model, constrs[drop])
|
|
||||||
if all_cuts === nothing
|
|
||||||
all_cuts = cuts
|
|
||||||
else
|
|
||||||
all_cuts.lhs = [all_cuts.lhs; cuts.lhs[keep, :]]
|
|
||||||
all_cuts.lb = [all_cuts.lb; cuts.lb[keep]]
|
|
||||||
all_cuts.lb = [all_cuts.lb; cuts.lb[keep]]
|
|
||||||
end
|
|
||||||
push!(stats_ncuts, length(all_cuts.lb))
|
|
||||||
|
|
||||||
undo_relax()
|
|
||||||
end
|
|
||||||
|
|
||||||
# Store cuts
|
|
||||||
if all_cuts !== nothing
|
|
||||||
@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.file.close()
|
|
||||||
end
|
|
||||||
|
|
||||||
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,
|
|
||||||
"obj_mip" => obj_mip,
|
|
||||||
"obj_lp" => obj_lp,
|
|
||||||
"stats_obj" => stats_obj,
|
|
||||||
"stats_gap" => stats_gap,
|
|
||||||
"stats_ncuts" => stats_ncuts,
|
|
||||||
)
|
|
||||||
end
|
|
||||||
|
|
||||||
export collect_gmi
|
|
||||||
@ -0,0 +1,625 @@
|
|||||||
|
# MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
|
||||||
|
# Copyright (C) 2020-2023, UChicago Argonne, LLC. All rights reserved.
|
||||||
|
# Released under the modified BSD license. See COPYING.md for more details.
|
||||||
|
|
||||||
|
using Printf
|
||||||
|
using JuMP
|
||||||
|
using HiGHS
|
||||||
|
using Random
|
||||||
|
using DataStructures
|
||||||
|
|
||||||
|
global ExpertDualGmiComponent = PyNULL()
|
||||||
|
global KnnDualGmiComponent = PyNULL()
|
||||||
|
|
||||||
|
Base.@kwdef mutable struct _KnnDualGmiData
|
||||||
|
k = nothing
|
||||||
|
extractor = nothing
|
||||||
|
train_h5 = nothing
|
||||||
|
model = nothing
|
||||||
|
strategy = nothing
|
||||||
|
end
|
||||||
|
|
||||||
|
function collect_gmi_dual(
|
||||||
|
mps_filename;
|
||||||
|
optimizer,
|
||||||
|
max_rounds = 10,
|
||||||
|
max_cuts_per_round = 500,
|
||||||
|
)
|
||||||
|
reset_timer!()
|
||||||
|
|
||||||
|
@timeit "Read H5" begin
|
||||||
|
h5_filename = replace(mps_filename, ".mps.gz" => ".h5")
|
||||||
|
h5 = H5File(h5_filename, "r")
|
||||||
|
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
|
||||||
|
|
||||||
|
# Define relative MIP gap
|
||||||
|
gap(v) = 100 * abs(obj_mip - v) / abs(obj_mip)
|
||||||
|
|
||||||
|
@timeit "Initialize" begin
|
||||||
|
stats_obj = []
|
||||||
|
stats_gap = []
|
||||||
|
stats_ncuts = []
|
||||||
|
original_basis = nothing
|
||||||
|
all_cuts = nothing
|
||||||
|
all_cuts_bases = nothing
|
||||||
|
all_cuts_rows = nothing
|
||||||
|
last_round_obj = nothing
|
||||||
|
end
|
||||||
|
|
||||||
|
@timeit "Read problem" begin
|
||||||
|
model = read_from_file(mps_filename)
|
||||||
|
set_optimizer(model, optimizer)
|
||||||
|
obj_original = objective_function(model)
|
||||||
|
end
|
||||||
|
|
||||||
|
for round = 1:max_rounds
|
||||||
|
@info "Round $(round)..."
|
||||||
|
|
||||||
|
@timeit "Convert model to standard form" begin
|
||||||
|
# Extract problem data
|
||||||
|
data = ProblemData(model)
|
||||||
|
|
||||||
|
# Construct optimal solution vector (with correct variable sequence)
|
||||||
|
sol_opt = [sol_opt_dict[n] for n in data.var_names]
|
||||||
|
|
||||||
|
# Assert optimal solution is feasible for the original problem
|
||||||
|
assert_leq(data.constr_lb, data.constr_lhs * sol_opt)
|
||||||
|
assert_leq(data.constr_lhs * sol_opt, data.constr_ub)
|
||||||
|
|
||||||
|
# Convert to standard form
|
||||||
|
data_s, transforms = convert_to_standard_form(data)
|
||||||
|
model_s = to_model(data_s)
|
||||||
|
set_optimizer(model_s, optimizer)
|
||||||
|
relax_integrality(model_s)
|
||||||
|
|
||||||
|
# Convert optimal solution to standard form
|
||||||
|
sol_opt_s = forward(transforms, sol_opt)
|
||||||
|
|
||||||
|
# Assert converted solution is feasible for standard form problem
|
||||||
|
assert_eq(data_s.constr_lhs * sol_opt_s, data_s.constr_lb)
|
||||||
|
end
|
||||||
|
|
||||||
|
@timeit "Optimize standard model" begin
|
||||||
|
@info "Optimizing standard model..."
|
||||||
|
optimize!(model_s)
|
||||||
|
obj = objective_value(model_s)
|
||||||
|
if round == 1
|
||||||
|
push!(stats_obj, obj)
|
||||||
|
push!(stats_gap, gap(obj))
|
||||||
|
push!(stats_ncuts, 0)
|
||||||
|
else
|
||||||
|
if obj ≈ last_round_obj
|
||||||
|
@info ("No improvement in obj value. Aborting.")
|
||||||
|
break
|
||||||
|
end
|
||||||
|
end
|
||||||
|
if termination_status(model_s) != MOI.OPTIMAL
|
||||||
|
error("Non-optimal termination status")
|
||||||
|
end
|
||||||
|
last_round_obj = obj
|
||||||
|
end
|
||||||
|
|
||||||
|
@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
|
||||||
|
|
||||||
|
@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, atol=1e-3)
|
||||||
|
assert_eq(tableau.lhs * sol_opt_s, tableau.rhs, atol=1e-3)
|
||||||
|
end
|
||||||
|
|
||||||
|
@timeit "Compute GMI cuts" begin
|
||||||
|
cuts_s = compute_gmi(data_s, tableau)
|
||||||
|
|
||||||
|
# Assert cuts have been generated correctly
|
||||||
|
assert_cuts_off(cuts_s, sol_frac)
|
||||||
|
assert_does_not_cut_off(cuts_s, sol_opt_s)
|
||||||
|
|
||||||
|
# Abort if no cuts are left
|
||||||
|
if length(cuts_s.lb) == 0
|
||||||
|
@info "No cuts generated. Aborting."
|
||||||
|
break
|
||||||
|
else
|
||||||
|
@info "Generated $(length(cuts_s.lb)) cuts"
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
@timeit "Add GMI cuts to original model" begin
|
||||||
|
@timeit "Convert to original form" begin
|
||||||
|
cuts = backwards(transforms, cuts_s)
|
||||||
|
end
|
||||||
|
|
||||||
|
@timeit "Prepare bv" begin
|
||||||
|
bv = repeat([basis], length(selected_rows))
|
||||||
|
end
|
||||||
|
|
||||||
|
@timeit "Append matrices" begin
|
||||||
|
if round == 1
|
||||||
|
all_cuts = cuts
|
||||||
|
all_cuts_bases = bv
|
||||||
|
all_cuts_rows = selected_rows
|
||||||
|
else
|
||||||
|
all_cuts.lhs = [all_cuts.lhs; cuts.lhs]
|
||||||
|
all_cuts.lb = [all_cuts.lb; cuts.lb]
|
||||||
|
all_cuts.ub = [all_cuts.ub; cuts.ub]
|
||||||
|
all_cuts_bases = [all_cuts_bases; bv]
|
||||||
|
all_cuts_rows = [all_cuts_rows; selected_rows]
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
@timeit "Add to model" begin
|
||||||
|
@info "Adding $(length(all_cuts.lb)) constraints to original model"
|
||||||
|
constrs, gmi_exps = add_constraint_set_dual_v2(model, all_cuts)
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
@timeit "Optimize original model" begin
|
||||||
|
set_objective_function(model, obj_original)
|
||||||
|
undo_relax = relax_integrality(model)
|
||||||
|
@info "Optimizing original model (constr)..."
|
||||||
|
optimize!(model)
|
||||||
|
obj = objective_value(model)
|
||||||
|
push!(stats_obj, obj)
|
||||||
|
push!(stats_gap, gap(obj))
|
||||||
|
sp = [shadow_price(c) for c in constrs]
|
||||||
|
undo_relax()
|
||||||
|
useful = [abs(sp[i]) > 1e-6 for (i, _) in enumerate(constrs)]
|
||||||
|
keep = findall(useful .== true)
|
||||||
|
end
|
||||||
|
|
||||||
|
@timeit "Filter out useless cuts" begin
|
||||||
|
@info "Keeping $(length(keep)) useful cuts"
|
||||||
|
all_cuts.lhs = all_cuts.lhs[keep, :]
|
||||||
|
all_cuts.lb = all_cuts.lb[keep]
|
||||||
|
all_cuts.ub = all_cuts.ub[keep]
|
||||||
|
all_cuts_bases = all_cuts_bases[keep, :]
|
||||||
|
all_cuts_rows = all_cuts_rows[keep, :]
|
||||||
|
push!(stats_ncuts, length(all_cuts_rows))
|
||||||
|
if isempty(keep)
|
||||||
|
break
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
@timeit "Update obj function of original model" begin
|
||||||
|
delete.(model, constrs)
|
||||||
|
set_objective_function(
|
||||||
|
model,
|
||||||
|
obj_original -
|
||||||
|
sum(sp[i] * gmi_exps[i] for (i, c) in enumerate(constrs) if useful[i]),
|
||||||
|
)
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
@timeit "Store cuts in H5 file" begin
|
||||||
|
if all_cuts !== nothing
|
||||||
|
ncuts = length(all_cuts_rows)
|
||||||
|
total =
|
||||||
|
length(original_basis.var_basic) +
|
||||||
|
length(original_basis.var_nonbasic) +
|
||||||
|
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
|
||||||
|
|
||||||
|
to = TimerOutputs.get_defaulttimer()
|
||||||
|
stats_time = TimerOutputs.tottime(to) / 1e9
|
||||||
|
print_timer()
|
||||||
|
|
||||||
|
return OrderedDict(
|
||||||
|
"instance" => mps_filename,
|
||||||
|
"max_rounds" => max_rounds,
|
||||||
|
"rounds" => length(stats_obj) - 1,
|
||||||
|
"obj_mip" => obj_mip,
|
||||||
|
"stats_obj" => stats_obj,
|
||||||
|
"stats_gap" => stats_gap,
|
||||||
|
"stats_ncuts" => stats_ncuts,
|
||||||
|
"stats_time" => stats_time,
|
||||||
|
)
|
||||||
|
end
|
||||||
|
|
||||||
|
function ExpertDualGmiComponent_before_mip(test_h5, model, _)
|
||||||
|
# Read cuts and optimal solution
|
||||||
|
h5 = H5File(test_h5, "r")
|
||||||
|
sol_opt_dict = Dict(
|
||||||
|
zip(
|
||||||
|
h5.get_array("static_var_names"),
|
||||||
|
convert(Array{Float64}, h5.get_array("mip_var_values")),
|
||||||
|
),
|
||||||
|
)
|
||||||
|
cut_basis_vars = h5.get_array("cuts_basis_vars")
|
||||||
|
cut_basis_sizes = h5.get_array("cuts_basis_sizes")
|
||||||
|
cut_rows = h5.get_array("cuts_rows")
|
||||||
|
obj_mip = h5.get_scalar("mip_lower_bound")
|
||||||
|
if obj_mip === nothing
|
||||||
|
obj_mip = h5.get_scalar("mip_obj_value")
|
||||||
|
end
|
||||||
|
h5.close()
|
||||||
|
|
||||||
|
# Initialize stats
|
||||||
|
stats_time_convert = 0
|
||||||
|
stats_time_tableau = 0
|
||||||
|
stats_time_gmi = 0
|
||||||
|
all_cuts = nothing
|
||||||
|
|
||||||
|
stats_time_convert = @elapsed begin
|
||||||
|
# Extract problem data
|
||||||
|
data = ProblemData(model)
|
||||||
|
|
||||||
|
# Construct optimal solution vector (with correct variable sequence)
|
||||||
|
sol_opt = [sol_opt_dict[n] for n in data.var_names]
|
||||||
|
|
||||||
|
# Assert optimal solution is feasible for the original problem
|
||||||
|
assert_leq(data.constr_lb, data.constr_lhs * sol_opt)
|
||||||
|
assert_leq(data.constr_lhs * sol_opt, data.constr_ub)
|
||||||
|
|
||||||
|
# Convert to standard form
|
||||||
|
data_s, transforms = convert_to_standard_form(data)
|
||||||
|
model_s = to_model(data_s)
|
||||||
|
set_optimizer(model_s, HiGHS.Optimizer)
|
||||||
|
relax_integrality(model_s)
|
||||||
|
|
||||||
|
# Convert optimal solution to standard form
|
||||||
|
sol_opt_s = forward(transforms, sol_opt)
|
||||||
|
|
||||||
|
# Assert converted solution is feasible for standard form problem
|
||||||
|
assert_eq(data_s.constr_lhs * sol_opt_s, data_s.constr_lb)
|
||||||
|
|
||||||
|
end
|
||||||
|
|
||||||
|
current_basis = nothing
|
||||||
|
for (r, row) in enumerate(cut_rows)
|
||||||
|
stats_time_tableau += @elapsed begin
|
||||||
|
if r == 1 || cut_basis_vars[r, :] != cut_basis_vars[r-1, :]
|
||||||
|
vbb, vnn, cbb, cnn = cut_basis_sizes[r, :]
|
||||||
|
current_basis = Basis(;
|
||||||
|
var_basic = cut_basis_vars[r, 1:vbb],
|
||||||
|
var_nonbasic = cut_basis_vars[r, vbb+1:vbb+vnn],
|
||||||
|
constr_basic = cut_basis_vars[r, vbb+vnn+1:vbb+vnn+cbb],
|
||||||
|
constr_nonbasic = cut_basis_vars[r, vbb+vnn+cbb+1:vbb+vnn+cbb+cnn],
|
||||||
|
)
|
||||||
|
end
|
||||||
|
tableau = compute_tableau(data_s, current_basis, rows = [row])
|
||||||
|
assert_eq(tableau.lhs * sol_opt_s, tableau.rhs)
|
||||||
|
end
|
||||||
|
stats_time_gmi += @elapsed begin
|
||||||
|
cuts_s = compute_gmi(data_s, tableau)
|
||||||
|
assert_does_not_cut_off(cuts_s, sol_opt_s)
|
||||||
|
end
|
||||||
|
cuts = backwards(transforms, cuts_s)
|
||||||
|
assert_does_not_cut_off(cuts, sol_opt)
|
||||||
|
|
||||||
|
if all_cuts === nothing
|
||||||
|
all_cuts = cuts
|
||||||
|
else
|
||||||
|
all_cuts.lhs = [all_cuts.lhs; cuts.lhs]
|
||||||
|
all_cuts.lb = [all_cuts.lb; cuts.lb]
|
||||||
|
all_cuts.ub = [all_cuts.ub; cuts.ub]
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
# Strategy 1: Add all cuts during the first call
|
||||||
|
function cut_callback_1(cb_data)
|
||||||
|
if all_cuts !== nothing
|
||||||
|
constrs = build_constraints(model, all_cuts)
|
||||||
|
@info "Enforcing $(length(constrs)) cuts..."
|
||||||
|
for c in constrs
|
||||||
|
MOI.submit(model, MOI.UserCut(cb_data), c)
|
||||||
|
end
|
||||||
|
all_cuts = nothing
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
# Strategy 2: Add violated cuts repeatedly until unable to separate
|
||||||
|
callback_disabled = false
|
||||||
|
function cut_callback_2(cb_data)
|
||||||
|
if callback_disabled
|
||||||
|
return
|
||||||
|
end
|
||||||
|
x = all_variables(model)
|
||||||
|
x_val = callback_value.(cb_data, x)
|
||||||
|
lhs_val = all_cuts.lhs * x_val
|
||||||
|
is_violated = lhs_val .> all_cuts.ub
|
||||||
|
selected_idx = findall(is_violated .== true)
|
||||||
|
selected_cuts = ConstraintSet(
|
||||||
|
lhs=all_cuts.lhs[selected_idx, :],
|
||||||
|
ub=all_cuts.ub[selected_idx],
|
||||||
|
lb=all_cuts.lb[selected_idx],
|
||||||
|
)
|
||||||
|
constrs = build_constraints(model, selected_cuts)
|
||||||
|
if length(constrs) > 0
|
||||||
|
@info "Enforcing $(length(constrs)) cuts..."
|
||||||
|
for c in constrs
|
||||||
|
MOI.submit(model, MOI.UserCut(cb_data), c)
|
||||||
|
end
|
||||||
|
else
|
||||||
|
@info "No violated cuts found. Disabling callback."
|
||||||
|
callback_disabled = true
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
# Set up cut callback
|
||||||
|
set_attribute(model, MOI.UserCutCallback(), cut_callback_1)
|
||||||
|
# set_attribute(model, MOI.UserCutCallback(), cut_callback_2)
|
||||||
|
|
||||||
|
stats = Dict()
|
||||||
|
stats["ExpertDualGmi: cuts"] = length(all_cuts.lb)
|
||||||
|
stats["ExpertDualGmi: time convert"] = stats_time_convert
|
||||||
|
stats["ExpertDualGmi: time tableau"] = stats_time_tableau
|
||||||
|
stats["ExpertDualGmi: time gmi"] = stats_time_gmi
|
||||||
|
return stats
|
||||||
|
end
|
||||||
|
|
||||||
|
function add_constraint_set_dual_v2(model::JuMP.Model, cs::ConstraintSet)
|
||||||
|
vars = all_variables(model)
|
||||||
|
nrows, ncols = size(cs.lhs)
|
||||||
|
|
||||||
|
@timeit "Transpose LHS" begin
|
||||||
|
lhs_t = spzeros(ncols, nrows)
|
||||||
|
ftranspose!(lhs_t, cs.lhs, x -> x)
|
||||||
|
lhs_t_rows = rowvals(lhs_t)
|
||||||
|
lhs_t_vals = nonzeros(lhs_t)
|
||||||
|
end
|
||||||
|
|
||||||
|
constrs = []
|
||||||
|
gmi_exps = []
|
||||||
|
for i = 1:nrows
|
||||||
|
c = nothing
|
||||||
|
gmi_exp = nothing
|
||||||
|
gmi_exp2 = nothing
|
||||||
|
@timeit "Build expr" begin
|
||||||
|
expr = AffExpr()
|
||||||
|
for k in nzrange(lhs_t, i)
|
||||||
|
add_to_expression!(expr, lhs_t_vals[k], vars[lhs_t_rows[k]])
|
||||||
|
end
|
||||||
|
end
|
||||||
|
@timeit "Add constraints" begin
|
||||||
|
if isinf(cs.ub[i])
|
||||||
|
c = @constraint(model, cs.lb[i] <= expr)
|
||||||
|
gmi_exp = cs.lb[i] - expr
|
||||||
|
elseif isinf(cs.lb[i])
|
||||||
|
c = @constraint(model, expr <= cs.ub[i])
|
||||||
|
gmi_exp = expr - cs.ub[i]
|
||||||
|
else
|
||||||
|
c = @constraint(model, cs.lb[i] <= expr <= cs.ub[i])
|
||||||
|
gmi_exp = cs.lb[i] - expr
|
||||||
|
gmi_exp2 = expr - cs.ub[i]
|
||||||
|
end
|
||||||
|
end
|
||||||
|
@timeit "Update structs" begin
|
||||||
|
push!(constrs, c)
|
||||||
|
push!(gmi_exps, gmi_exp)
|
||||||
|
if !isnothing(gmi_exp2)
|
||||||
|
push!(gmi_exps, gmi_exp2)
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
return constrs, gmi_exps
|
||||||
|
end
|
||||||
|
|
||||||
|
function _dualgmi_features(h5_filename, extractor)
|
||||||
|
h5 = H5File(h5_filename, "r")
|
||||||
|
try
|
||||||
|
return extractor.get_instance_features(h5)
|
||||||
|
finally
|
||||||
|
h5.close()
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
function _dualgmi_generate(train_h5, model)
|
||||||
|
@timeit "Read problem data" begin
|
||||||
|
data = ProblemData(model)
|
||||||
|
end
|
||||||
|
@timeit "Convert to standard form" begin
|
||||||
|
data_s, transforms = convert_to_standard_form(data)
|
||||||
|
end
|
||||||
|
|
||||||
|
@timeit "Collect cuts from H5 files" begin
|
||||||
|
vars_to_unique_basis_offset = Dict()
|
||||||
|
unique_basis_vars = nothing
|
||||||
|
unique_basis_sizes = nothing
|
||||||
|
unique_basis_rows = nothing
|
||||||
|
|
||||||
|
for h5_filename in train_h5
|
||||||
|
h5 = H5File(h5_filename, "r")
|
||||||
|
cut_basis_vars = h5.get_array("cuts_basis_vars")
|
||||||
|
cut_basis_sizes = h5.get_array("cuts_basis_sizes")
|
||||||
|
cut_rows = h5.get_array("cuts_rows")
|
||||||
|
ncuts, nvars = size(cut_basis_vars)
|
||||||
|
if unique_basis_vars === nothing
|
||||||
|
unique_basis_vars = Matrix{Int}(undef, 0, nvars)
|
||||||
|
unique_basis_sizes = Matrix{Int}(undef, 0, 4)
|
||||||
|
unique_basis_rows = Dict{Int,Set{Int}}()
|
||||||
|
end
|
||||||
|
for i in 1:ncuts
|
||||||
|
vars = cut_basis_vars[i, :]
|
||||||
|
sizes = cut_basis_sizes[i, :]
|
||||||
|
row = cut_rows[i]
|
||||||
|
if vars ∉ keys(vars_to_unique_basis_offset)
|
||||||
|
offset = size(unique_basis_vars)[1] + 1
|
||||||
|
vars_to_unique_basis_offset[vars] = offset
|
||||||
|
unique_basis_vars = [unique_basis_vars; vars']
|
||||||
|
unique_basis_sizes = [unique_basis_sizes; sizes']
|
||||||
|
unique_basis_rows[offset] = Set()
|
||||||
|
end
|
||||||
|
offset = vars_to_unique_basis_offset[vars]
|
||||||
|
push!(unique_basis_rows[offset], row)
|
||||||
|
end
|
||||||
|
h5.close()
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
@timeit "Compute tableaus and cuts" begin
|
||||||
|
all_cuts = nothing
|
||||||
|
for (offset, rows) in unique_basis_rows
|
||||||
|
try
|
||||||
|
vbb, vnn, cbb, cnn = unique_basis_sizes[offset, :]
|
||||||
|
current_basis = Basis(;
|
||||||
|
var_basic = unique_basis_vars[offset, 1:vbb],
|
||||||
|
var_nonbasic = unique_basis_vars[offset, vbb+1:vbb+vnn],
|
||||||
|
constr_basic = unique_basis_vars[offset, vbb+vnn+1:vbb+vnn+cbb],
|
||||||
|
constr_nonbasic = unique_basis_vars[offset, vbb+vnn+cbb+1:vbb+vnn+cbb+cnn],
|
||||||
|
)
|
||||||
|
|
||||||
|
tableau = compute_tableau(data_s, current_basis; rows=collect(rows))
|
||||||
|
cuts_s = compute_gmi(data_s, tableau)
|
||||||
|
cuts = backwards(transforms, cuts_s)
|
||||||
|
if all_cuts === nothing
|
||||||
|
all_cuts = cuts
|
||||||
|
else
|
||||||
|
all_cuts.lhs = [all_cuts.lhs; cuts.lhs]
|
||||||
|
all_cuts.lb = [all_cuts.lb; cuts.lb]
|
||||||
|
all_cuts.ub = [all_cuts.ub; cuts.ub]
|
||||||
|
end
|
||||||
|
catch e
|
||||||
|
if isa(e, AssertionError)
|
||||||
|
@warn "Numerical error detected. Skipping cuts from current tableau."
|
||||||
|
continue
|
||||||
|
else
|
||||||
|
rethrow(e)
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
return all_cuts
|
||||||
|
end
|
||||||
|
|
||||||
|
function _dualgmi_set_callback(model, all_cuts)
|
||||||
|
function cut_callback(cb_data)
|
||||||
|
if all_cuts !== nothing
|
||||||
|
constrs = build_constraints(model, all_cuts)
|
||||||
|
@info "Enforcing $(length(constrs)) cuts..."
|
||||||
|
for c in constrs
|
||||||
|
MOI.submit(model, MOI.UserCut(cb_data), c)
|
||||||
|
end
|
||||||
|
all_cuts = nothing
|
||||||
|
end
|
||||||
|
end
|
||||||
|
set_attribute(model, MOI.UserCutCallback(), cut_callback)
|
||||||
|
end
|
||||||
|
|
||||||
|
function KnnDualGmiComponent_fit(data::_KnnDualGmiData, train_h5)
|
||||||
|
x = hcat([_dualgmi_features(filename, data.extractor) for filename in train_h5]...)'
|
||||||
|
model = pyimport("sklearn.neighbors").NearestNeighbors(n_neighbors = length(train_h5))
|
||||||
|
model.fit(x)
|
||||||
|
data.model = model
|
||||||
|
data.train_h5 = train_h5
|
||||||
|
end
|
||||||
|
|
||||||
|
|
||||||
|
function KnnDualGmiComponent_before_mip(data::_KnnDualGmiData, test_h5, model, _)
|
||||||
|
reset_timer!()
|
||||||
|
|
||||||
|
@timeit "Extract features" begin
|
||||||
|
x = _dualgmi_features(test_h5, data.extractor)
|
||||||
|
x = reshape(x, 1, length(x))
|
||||||
|
end
|
||||||
|
|
||||||
|
@timeit "Find neighbors" begin
|
||||||
|
neigh_dist, neigh_ind = data.model.kneighbors(x, return_distance = true)
|
||||||
|
neigh_ind = neigh_ind .+ 1
|
||||||
|
N = length(neigh_dist)
|
||||||
|
k = min(N, data.k)
|
||||||
|
|
||||||
|
if data.strategy == "near"
|
||||||
|
selected = collect(1:k)
|
||||||
|
elseif data.strategy == "far"
|
||||||
|
selected = collect((N - k + 1) : N)
|
||||||
|
elseif data.strategy == "rand"
|
||||||
|
selected = shuffle(collect(1:N))[1:k]
|
||||||
|
else
|
||||||
|
error("unknown strategy: $(data.strategy)")
|
||||||
|
end
|
||||||
|
|
||||||
|
@info "Dual GMI: Selected neighbors ($(data.strategy)):"
|
||||||
|
neigh_dist = neigh_dist[selected]
|
||||||
|
neigh_ind = neigh_ind[selected]
|
||||||
|
for i in 1:k
|
||||||
|
h5_filename = data.train_h5[neigh_ind[i]]
|
||||||
|
dist = neigh_dist[i]
|
||||||
|
@info " $(h5_filename) dist=$(dist)"
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
@info "Dual GMI: Generating cuts..."
|
||||||
|
@timeit "Generate cuts" begin
|
||||||
|
time_generate = @elapsed begin
|
||||||
|
cuts = _dualgmi_generate(data.train_h5[neigh_ind], model)
|
||||||
|
end
|
||||||
|
@info "Dual GMI: Generated $(length(cuts.lb)) unique cuts in $(time_generate) seconds"
|
||||||
|
end
|
||||||
|
|
||||||
|
@timeit "Set callback" begin
|
||||||
|
_dualgmi_set_callback(model, cuts)
|
||||||
|
end
|
||||||
|
|
||||||
|
print_timer()
|
||||||
|
|
||||||
|
stats = Dict()
|
||||||
|
stats["KnnDualGmi: k"] = k
|
||||||
|
stats["KnnDualGmi: strategy"] = data.strategy
|
||||||
|
stats["KnnDualGmi: cuts"] = length(cuts.lb)
|
||||||
|
stats["KnnDualGmi: time generate"] = time_generate
|
||||||
|
return stats
|
||||||
|
end
|
||||||
|
|
||||||
|
function __init_gmi_dual__()
|
||||||
|
@pydef mutable struct Class1
|
||||||
|
function fit(_, _) end
|
||||||
|
function before_mip(self, test_h5, model, stats)
|
||||||
|
ExpertDualGmiComponent_before_mip(test_h5, model.inner, stats)
|
||||||
|
end
|
||||||
|
end
|
||||||
|
copy!(ExpertDualGmiComponent, Class1)
|
||||||
|
|
||||||
|
@pydef mutable struct Class2
|
||||||
|
function __init__(self; extractor, k = 3, strategy = "near")
|
||||||
|
self.data = _KnnDualGmiData(; extractor, k, strategy)
|
||||||
|
end
|
||||||
|
function fit(self, train_h5)
|
||||||
|
KnnDualGmiComponent_fit(self.data, train_h5)
|
||||||
|
end
|
||||||
|
function before_mip(self, test_h5, model, stats)
|
||||||
|
return @time KnnDualGmiComponent_before_mip(self.data, test_h5, model.inner, stats)
|
||||||
|
end
|
||||||
|
end
|
||||||
|
copy!(KnnDualGmiComponent, Class2)
|
||||||
|
end
|
||||||
|
|
||||||
|
export collect_gmi_dual, expert_gmi_dual, ExpertDualGmiComponent, KnnDualGmiComponent
|
||||||
|
|
||||||
@ -0,0 +1,51 @@
|
|||||||
|
@inline frac(x::Float64) = x - floor(x)
|
||||||
|
|
||||||
|
@inline frac2(x::Float64) = ceil(x) - x
|
||||||
|
|
||||||
|
function assert_leq(a, b; atol = 0.01)
|
||||||
|
if !all(a .<= b .+ atol)
|
||||||
|
delta = a .- b
|
||||||
|
for i in eachindex(delta)
|
||||||
|
if delta[i] > atol
|
||||||
|
@info "Assertion failed: a[$i] = $(a[i]) <= $(b[i]) = b[$i]"
|
||||||
|
end
|
||||||
|
end
|
||||||
|
error("assert_leq failed")
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
function assert_eq(a, b; atol = 1e-4)
|
||||||
|
if !all(abs.(a .- b) .<= atol)
|
||||||
|
delta = abs.(a .- b)
|
||||||
|
for i in eachindex(delta)
|
||||||
|
if delta[i] > atol
|
||||||
|
@info "Assertion failed: a[$i] = $(a[i]) == $(b[i]) = b[$i]"
|
||||||
|
end
|
||||||
|
end
|
||||||
|
error("assert_eq failed")
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
function assert_cuts_off(cuts::ConstraintSet, x::Vector{Float64}, tol = 1e-6)
|
||||||
|
for i = 1:length(cuts.lb)
|
||||||
|
val = cuts.lhs[i, :]' * x
|
||||||
|
if (val <= cuts.ub[i] - tol) && (val >= cuts.lb[i] + tol)
|
||||||
|
throw(ErrorException("inequality fails to cut off fractional solution"))
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
function assert_does_not_cut_off(cuts::ConstraintSet, x::Vector{Float64}; tol = 1e-6)
|
||||||
|
for i = 1:length(cuts.lb)
|
||||||
|
val = cuts.lhs[i, :]' * x
|
||||||
|
ub = cuts.ub[i]
|
||||||
|
lb = cuts.lb[i]
|
||||||
|
if (val >= ub) || (val <= lb)
|
||||||
|
throw(
|
||||||
|
ErrorException(
|
||||||
|
"inequality $i cuts off integer solution ($lb <= $val <= $ub)",
|
||||||
|
),
|
||||||
|
)
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
@ -0,0 +1,59 @@
|
|||||||
|
# MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
|
||||||
|
# Copyright (C) 2020-2024, UChicago Argonne, LLC. All rights reserved.
|
||||||
|
# Released under the modified BSD license. See COPYING.md for more details.
|
||||||
|
|
||||||
|
using JuMP
|
||||||
|
using HiGHS
|
||||||
|
|
||||||
|
global MaxWeightStableSetData = PyNULL()
|
||||||
|
global MaxWeightStableSetGenerator = PyNULL()
|
||||||
|
|
||||||
|
function __init_problems_stab__()
|
||||||
|
copy!(MaxWeightStableSetData, pyimport("miplearn.problems.stab").MaxWeightStableSetData)
|
||||||
|
copy!(
|
||||||
|
MaxWeightStableSetGenerator,
|
||||||
|
pyimport("miplearn.problems.stab").MaxWeightStableSetGenerator,
|
||||||
|
)
|
||||||
|
end
|
||||||
|
|
||||||
|
function build_stab_model_jump(data::Any; optimizer = HiGHS.Optimizer)
|
||||||
|
nx = pyimport("networkx")
|
||||||
|
|
||||||
|
if data isa String
|
||||||
|
data = read_pkl_gz(data)
|
||||||
|
end
|
||||||
|
model = Model(optimizer)
|
||||||
|
|
||||||
|
# Variables and objective function
|
||||||
|
nodes = data.graph.nodes
|
||||||
|
x = @variable(model, x[nodes], Bin)
|
||||||
|
@objective(model, Min, sum(-data.weights[i+1] * x[i] for i in nodes))
|
||||||
|
|
||||||
|
# Edge inequalities
|
||||||
|
for (i1, i2) in data.graph.edges
|
||||||
|
@constraint(model, x[i1] + x[i2] <= 1, base_name = "eq_edge[$i1,$i2]")
|
||||||
|
end
|
||||||
|
|
||||||
|
function cuts_separate(cb_data)
|
||||||
|
x_val = callback_value.(Ref(cb_data), x)
|
||||||
|
violations = []
|
||||||
|
for clique in nx.find_cliques(data.graph)
|
||||||
|
if sum(x_val[i] for i in clique) > 1.0001
|
||||||
|
push!(violations, sort(clique))
|
||||||
|
end
|
||||||
|
end
|
||||||
|
return violations
|
||||||
|
end
|
||||||
|
|
||||||
|
function cuts_enforce(violations)
|
||||||
|
@info "Adding $(length(violations)) clique cuts..."
|
||||||
|
for clique in violations
|
||||||
|
constr = @build_constraint(sum(x[i] for i in clique) <= 1)
|
||||||
|
submit(model, constr)
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
return JumpModel(model, cuts_separate = cuts_separate, cuts_enforce = cuts_enforce)
|
||||||
|
end
|
||||||
|
|
||||||
|
export MaxWeightStableSetData, MaxWeightStableSetGenerator, build_stab_model_jump
|
||||||
@ -0,0 +1,71 @@
|
|||||||
|
# MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
|
||||||
|
# Copyright (C) 2020-2024, UChicago Argonne, LLC. All rights reserved.
|
||||||
|
# Released under the modified BSD license. See COPYING.md for more details.
|
||||||
|
|
||||||
|
using JuMP
|
||||||
|
|
||||||
|
global TravelingSalesmanData = PyNULL()
|
||||||
|
global TravelingSalesmanGenerator = PyNULL()
|
||||||
|
|
||||||
|
function __init_problems_tsp__()
|
||||||
|
copy!(TravelingSalesmanData, pyimport("miplearn.problems.tsp").TravelingSalesmanData)
|
||||||
|
copy!(
|
||||||
|
TravelingSalesmanGenerator,
|
||||||
|
pyimport("miplearn.problems.tsp").TravelingSalesmanGenerator,
|
||||||
|
)
|
||||||
|
end
|
||||||
|
|
||||||
|
function build_tsp_model_jump(data::Any; optimizer)
|
||||||
|
nx = pyimport("networkx")
|
||||||
|
|
||||||
|
if data isa String
|
||||||
|
data = read_pkl_gz(data)
|
||||||
|
end
|
||||||
|
model = Model(optimizer)
|
||||||
|
edges = [(i, j) for i = 1:data.n_cities for j = (i+1):data.n_cities]
|
||||||
|
x = @variable(model, x[edges], Bin)
|
||||||
|
@objective(model, Min, sum(x[(i, j)] * data.distances[i, j] for (i, j) in edges))
|
||||||
|
|
||||||
|
# Eq: Must choose two edges adjacent to each node
|
||||||
|
@constraint(
|
||||||
|
model,
|
||||||
|
eq_degree[i in 1:data.n_cities],
|
||||||
|
sum(x[(min(i, j), max(i, j))] for j = 1:data.n_cities if i != j) == 2
|
||||||
|
)
|
||||||
|
|
||||||
|
function lazy_separate(cb_data)
|
||||||
|
x_val = callback_value.(Ref(cb_data), x)
|
||||||
|
violations = []
|
||||||
|
selected_edges = [e for e in edges if x_val[e] > 0.5]
|
||||||
|
graph = nx.Graph()
|
||||||
|
graph.add_edges_from(selected_edges)
|
||||||
|
for component in nx.connected_components(graph)
|
||||||
|
if length(component) < data.n_cities
|
||||||
|
cut_edges = [
|
||||||
|
[e[1], e[2]] for
|
||||||
|
e in edges if (e[1] ∈ component && e[2] ∉ component) ||
|
||||||
|
(e[1] ∉ component && e[2] ∈ component)
|
||||||
|
]
|
||||||
|
push!(violations, cut_edges)
|
||||||
|
end
|
||||||
|
end
|
||||||
|
return violations
|
||||||
|
end
|
||||||
|
|
||||||
|
function lazy_enforce(violations)
|
||||||
|
@info "Adding $(length(violations)) subtour elimination eqs..."
|
||||||
|
for violation in violations
|
||||||
|
constr = @build_constraint(sum(x[(e[1], e[2])] for e in violation) >= 2)
|
||||||
|
submit(model, constr)
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
return JumpModel(
|
||||||
|
model,
|
||||||
|
lazy_enforce = lazy_enforce,
|
||||||
|
lazy_separate = lazy_separate,
|
||||||
|
lp_optimizer = optimizer,
|
||||||
|
)
|
||||||
|
end
|
||||||
|
|
||||||
|
export TravelingSalesmanData, TravelingSalesmanGenerator, build_tsp_model_jump
|
||||||
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
@ -0,0 +1,23 @@
|
|||||||
|
# MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
|
||||||
|
# Copyright (C) 2020-2024, UChicago Argonne, LLC. All rights reserved.
|
||||||
|
# Released under the modified BSD license. See COPYING.md for more details.
|
||||||
|
|
||||||
|
using HiGHS
|
||||||
|
|
||||||
|
function test_cuts_tableau_gmi()
|
||||||
|
mps_filename = "$BASEDIR/../fixtures/bell5.mps.gz"
|
||||||
|
h5_filename = "$BASEDIR/../fixtures/bell5.h5"
|
||||||
|
collect_gmi(mps_filename, optimizer = HiGHS.Optimizer)
|
||||||
|
h5 = H5File(h5_filename, "r")
|
||||||
|
try
|
||||||
|
cuts_lb = h5.get_array("cuts_lb")
|
||||||
|
cuts_ub = h5.get_array("cuts_ub")
|
||||||
|
cuts_lhs = h5.get_sparse("cuts_lhs")
|
||||||
|
n_cuts = length(cuts_lb)
|
||||||
|
@test n_cuts > 0
|
||||||
|
@test n_cuts == length(cuts_ub)
|
||||||
|
@test cuts_lhs.shape[1] == n_cuts
|
||||||
|
finally
|
||||||
|
h5.close()
|
||||||
|
end
|
||||||
|
end
|
||||||
@ -0,0 +1,70 @@
|
|||||||
|
# MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
|
||||||
|
# Copyright (C) 2020-2024, UChicago Argonne, LLC. All rights reserved.
|
||||||
|
# Released under the modified BSD license. See COPYING.md for more details.
|
||||||
|
|
||||||
|
using SCIP
|
||||||
|
using HiGHS
|
||||||
|
using MIPLearn.Cuts
|
||||||
|
|
||||||
|
function test_cuts_tableau_gmi_dual_collect()
|
||||||
|
mps_filename = "$BASEDIR/../fixtures/bell5.mps.gz"
|
||||||
|
h5_filename = "$BASEDIR/../fixtures/bell5.h5"
|
||||||
|
stats = collect_gmi_dual(mps_filename, optimizer = HiGHS.Optimizer)
|
||||||
|
h5 = H5File(h5_filename, "r")
|
||||||
|
try
|
||||||
|
cuts_basis_vars = h5.get_array("cuts_basis_vars")
|
||||||
|
cuts_basis_sizes = h5.get_array("cuts_basis_sizes")
|
||||||
|
cuts_rows = h5.get_array("cuts_rows")
|
||||||
|
@test size(cuts_basis_vars) == (15, 402)
|
||||||
|
@test size(cuts_basis_sizes) == (15, 4)
|
||||||
|
@test size(cuts_rows) == (15,)
|
||||||
|
finally
|
||||||
|
h5.close()
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
function test_cuts_tableau_gmi_dual_usage()
|
||||||
|
function build_model(mps_filename)
|
||||||
|
model = read_from_file(mps_filename)
|
||||||
|
set_optimizer(model, SCIP.Optimizer)
|
||||||
|
return JumpModel(model)
|
||||||
|
end
|
||||||
|
|
||||||
|
mps_filename = "$BASEDIR/../fixtures/bell5.mps.gz"
|
||||||
|
h5_filename = "$BASEDIR/../fixtures/bell5.h5"
|
||||||
|
rm(h5_filename, force=true)
|
||||||
|
|
||||||
|
# Run basic collector
|
||||||
|
bc = BasicCollector(write_mps = false, skip_lp = true)
|
||||||
|
bc.collect([mps_filename], build_model)
|
||||||
|
|
||||||
|
# Run dual GMI collector
|
||||||
|
@info "Running dual GMI collector..."
|
||||||
|
collect_gmi_dual(mps_filename, optimizer = HiGHS.Optimizer)
|
||||||
|
|
||||||
|
# # Test expert component
|
||||||
|
# solver = LearningSolver(
|
||||||
|
# components = [
|
||||||
|
# ExpertPrimalComponent(action = SetWarmStart()),
|
||||||
|
# ExpertDualGmiComponent(),
|
||||||
|
# ],
|
||||||
|
# skip_lp = true,
|
||||||
|
# )
|
||||||
|
# solver.optimize(mps_filename, build_model)
|
||||||
|
|
||||||
|
# Test kNN component
|
||||||
|
knn = KnnDualGmiComponent(
|
||||||
|
extractor = H5FieldsExtractor(instance_fields = ["static_var_obj_coeffs"]),
|
||||||
|
k = 2,
|
||||||
|
)
|
||||||
|
knn.fit([h5_filename, h5_filename])
|
||||||
|
solver = LearningSolver(
|
||||||
|
components = [
|
||||||
|
ExpertPrimalComponent(action = SetWarmStart()),
|
||||||
|
knn,
|
||||||
|
],
|
||||||
|
skip_lp = true,
|
||||||
|
)
|
||||||
|
solver.optimize(mps_filename, build_model)
|
||||||
|
return
|
||||||
|
end
|
||||||
@ -0,0 +1,41 @@
|
|||||||
|
# MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
|
||||||
|
# Copyright (C) 2020-2024, UChicago Argonne, LLC. All rights reserved.
|
||||||
|
# Released under the modified BSD license. See COPYING.md for more details.
|
||||||
|
|
||||||
|
using SCIP
|
||||||
|
|
||||||
|
function gen_stab()
|
||||||
|
np = pyimport("numpy")
|
||||||
|
uniform = pyimport("scipy.stats").uniform
|
||||||
|
randint = pyimport("scipy.stats").randint
|
||||||
|
np.random.seed(42)
|
||||||
|
gen = MaxWeightStableSetGenerator(
|
||||||
|
w = uniform(10.0, scale = 1.0),
|
||||||
|
n = randint(low = 50, high = 51),
|
||||||
|
p = uniform(loc = 0.5, scale = 0.0),
|
||||||
|
fix_graph = true,
|
||||||
|
)
|
||||||
|
data = gen.generate(1)
|
||||||
|
data_filenames = write_pkl_gz(data, "$BASEDIR/../fixtures", prefix = "stab-n50-")
|
||||||
|
collector = BasicCollector()
|
||||||
|
collector.collect(
|
||||||
|
data_filenames,
|
||||||
|
data -> build_stab_model_jump(data, optimizer = SCIP.Optimizer),
|
||||||
|
progress = true,
|
||||||
|
verbose = true,
|
||||||
|
)
|
||||||
|
end
|
||||||
|
|
||||||
|
function test_cuts()
|
||||||
|
data_filenames = ["$BASEDIR/../fixtures/stab-n50-00000.pkl.gz"]
|
||||||
|
clf = pyimport("sklearn.dummy").DummyClassifier()
|
||||||
|
extractor = H5FieldsExtractor(instance_fields = ["static_var_obj_coeffs"])
|
||||||
|
comp = MemorizingCutsComponent(clf = clf, extractor = extractor)
|
||||||
|
solver = LearningSolver(components = [comp])
|
||||||
|
solver.fit(data_filenames)
|
||||||
|
model, stats = solver.optimize(
|
||||||
|
data_filenames[1],
|
||||||
|
data -> build_stab_model_jump(data, optimizer = SCIP.Optimizer),
|
||||||
|
)
|
||||||
|
@test stats["Cuts: AOT"] > 0
|
||||||
|
end
|
||||||
@ -0,0 +1,44 @@
|
|||||||
|
# MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
|
||||||
|
# Copyright (C) 2020-2024, UChicago Argonne, LLC. All rights reserved.
|
||||||
|
# Released under the modified BSD license. See COPYING.md for more details.
|
||||||
|
|
||||||
|
using GLPK
|
||||||
|
|
||||||
|
function gen_tsp()
|
||||||
|
np = pyimport("numpy")
|
||||||
|
uniform = pyimport("scipy.stats").uniform
|
||||||
|
randint = pyimport("scipy.stats").randint
|
||||||
|
np.random.seed(42)
|
||||||
|
|
||||||
|
gen = TravelingSalesmanGenerator(
|
||||||
|
x = uniform(loc = 0.0, scale = 1000.0),
|
||||||
|
y = uniform(loc = 0.0, scale = 1000.0),
|
||||||
|
n = randint(low = 20, high = 21),
|
||||||
|
gamma = uniform(loc = 1.0, scale = 0.25),
|
||||||
|
fix_cities = true,
|
||||||
|
round = true,
|
||||||
|
)
|
||||||
|
data = gen.generate(1)
|
||||||
|
data_filenames = write_pkl_gz(data, "$BASEDIR/../fixtures", prefix = "tsp-n20-")
|
||||||
|
collector = BasicCollector()
|
||||||
|
collector.collect(
|
||||||
|
data_filenames,
|
||||||
|
data -> build_tsp_model_jump(data, optimizer = GLPK.Optimizer),
|
||||||
|
progress = true,
|
||||||
|
verbose = true,
|
||||||
|
)
|
||||||
|
end
|
||||||
|
|
||||||
|
function test_lazy()
|
||||||
|
data_filenames = ["$BASEDIR/../fixtures/tsp-n20-00000.pkl.gz"]
|
||||||
|
clf = pyimport("sklearn.dummy").DummyClassifier()
|
||||||
|
extractor = H5FieldsExtractor(instance_fields = ["static_var_obj_coeffs"])
|
||||||
|
comp = MemorizingLazyComponent(clf = clf, extractor = extractor)
|
||||||
|
solver = LearningSolver(components = [comp])
|
||||||
|
solver.fit(data_filenames)
|
||||||
|
model, stats = solver.optimize(
|
||||||
|
data_filenames[1],
|
||||||
|
data -> build_tsp_model_jump(data, optimizer = GLPK.Optimizer),
|
||||||
|
)
|
||||||
|
@test stats["Lazy Constraints: AOT"] > 0
|
||||||
|
end
|
||||||
@ -0,0 +1,22 @@
|
|||||||
|
# MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
|
||||||
|
# Copyright (C) 2020-2024, UChicago Argonne, LLC. All rights reserved.
|
||||||
|
# Released under the modified BSD license. See COPYING.md for more details.
|
||||||
|
|
||||||
|
using PyCall
|
||||||
|
using SCIP
|
||||||
|
|
||||||
|
function test_problems_stab()
|
||||||
|
nx = pyimport("networkx")
|
||||||
|
data = MaxWeightStableSetData(
|
||||||
|
graph = nx.gnp_random_graph(25, 0.5, seed = 42),
|
||||||
|
weights = repeat([1.0], 25),
|
||||||
|
)
|
||||||
|
h5 = H5File(tempname(), "w")
|
||||||
|
model = build_stab_model_jump(data, optimizer = SCIP.Optimizer)
|
||||||
|
model.extract_after_load(h5)
|
||||||
|
model.optimize()
|
||||||
|
model.extract_after_mip(h5)
|
||||||
|
@test h5.get_scalar("mip_obj_value") == -6
|
||||||
|
@test h5.get_scalar("mip_cuts")[1:20] == "[[0,8,11,13],[0,8,13"
|
||||||
|
h5.close()
|
||||||
|
end
|
||||||
@ -0,0 +1,22 @@
|
|||||||
|
# MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
|
||||||
|
# Copyright (C) 2020-2024, UChicago Argonne, LLC. All rights reserved.
|
||||||
|
# Released under the modified BSD license. See COPYING.md for more details.
|
||||||
|
|
||||||
|
using GLPK
|
||||||
|
using JuMP
|
||||||
|
|
||||||
|
function test_problems_tsp()
|
||||||
|
pdist = pyimport("scipy.spatial.distance").pdist
|
||||||
|
squareform = pyimport("scipy.spatial.distance").squareform
|
||||||
|
|
||||||
|
data = TravelingSalesmanData(
|
||||||
|
n_cities = 6,
|
||||||
|
distances = squareform(
|
||||||
|
pdist([[0.0, 0.0], [1.0, 0.0], [2.0, 0.0], [3.0, 0.0], [0.0, 1.0], [3.0, 1.0]]),
|
||||||
|
),
|
||||||
|
)
|
||||||
|
model = build_tsp_model_jump(data, optimizer = GLPK.Optimizer)
|
||||||
|
model.optimize()
|
||||||
|
@test objective_value(model.inner) == 8.0
|
||||||
|
return
|
||||||
|
end
|
||||||
Loading…
Reference in new issue