diff --git a/src/Cuts/tableau/gmi_dual.jl b/src/Cuts/tableau/gmi_dual.jl index 06fe13c..4250ecf 100644 --- a/src/Cuts/tableau/gmi_dual.jl +++ b/src/Cuts/tableau/gmi_dual.jl @@ -269,8 +269,8 @@ function collect_gmi_FisSal2011( optimizer, max_rounds = 10_000, max_cuts_per_round = 1_000_000, - time_limit = 1_000_000, - interval_print=1, + time_limit = 900, + print_interval_sec=1, silent_solver=true, max_pool_size_mb = 1024, ) @@ -291,41 +291,28 @@ function collect_gmi_FisSal2011( end @timeit "Initialize" begin - backtrack_count = 0 - deterioration_count = 0 - # ε = 0.01 - basis_curr = nothing - basis_prev = nothing basis_seen = Set{UInt64}() + count_backtrack = 0 + count_deterioration = 0 gapcl_best = 0 gapcl_curr = 0 - multipliers_curr = nothing - multipliers_best = nothing last_print_time = 0 - obj_curr = 0 - obj_initial = nothing + multipliers_best = nothing + multipliers_curr = nothing obj_best = nothing + obj_curr = nothing obj_hist = [] + obj_initial = nothing 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 + μ = 10 end - gap(v) = 100 * abs(obj_mip - v) / abs(obj_mip) 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 model = read_from_file(mps_filename) set_optimizer(model, optimizer) @@ -361,13 +348,11 @@ function collect_gmi_FisSal2011( 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(multipliers, 0.1) v = sparse(pool.lhs * multipliers_curr) lagr_term = AffExpr(dot(multipliers_curr, pool.lb)) for offset in 1:nnz(v) @@ -391,7 +376,6 @@ function collect_gmi_FisSal2011( 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)") push!(obj_hist, obj_curr) if length(obj_hist) > 100 @@ -407,28 +391,21 @@ function collect_gmi_FisSal2011( 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 - push!(stats_ncuts, 0) - else - push!(stats_ncuts, length(pool.lb)) - end if termination_status(model_s) != MOI.OPTIMAL error("Non-optimal termination status") end Δ = obj_mip - obj_best if obj_curr < obj_best - Δ - deterioration_count += 1 + count_deterioration += 1 else - deterioration_count = 0 + count_deterioration = 0 end - if deterioration_count >= 10 + if count_deterioration >= 10 μ *= 0.5 multipliers_curr = multipliers_best - deterioration_count = 0 - backtrack_count += 1 + count_deterioration = 0 + count_backtrack += 1 elseif length(obj_hist) >= 100 obj_hist_avg = mean(obj_hist) improv = obj_best - obj_hist[1] @@ -446,18 +423,15 @@ function collect_gmi_FisSal2011( 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) + basis = get_basis(model_s) + basis_hash = hash(basis) push!(basis_seen, basis_hash) 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 @timeit "Compute tableau rows" begin - tableau = compute_tableau(data_s, basis_curr, x = sol_frac, rows = selected_rows) - - # Assert tableau rows have been computed correctly + tableau = compute_tableau(data_s, basis, x = sol_frac, rows = selected_rows) assert_eq(tableau.lhs * sol_frac, tableau.rhs, atol=1e-3) assert_eq(tableau.lhs * sol_opt_s, tableau.rhs, atol=1e-3) end @@ -507,14 +481,15 @@ function collect_gmi_FisSal2011( end end - selected_idx = [] - selected_contrs = [] - if round == 1 || round == max_rounds + if mod(round, 1000) == 1 || round == max_rounds @timeit "Update multipliers (large LP)" begin + selected_idx = [] + selected_contrs = [] while true @timeit "Optimize LP (extended)" begin set_objective_function(model_s, orig_obj_s) optimize!(model_s) + obj_curr = objective_value(model_s) sol_frac = get_x(model_s) end @@ -550,8 +525,6 @@ function collect_gmi_FisSal2011( end end - pool_active = length(selected_idx) - @timeit "Find dual values for all selected cuts" begin multipliers_curr .= 0 pool_cut_age .+= 1 @@ -563,6 +536,15 @@ function collect_gmi_FisSal2011( 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 pool_size_mb = Base.summarysize(pool) / 1024^2 while pool_size_mb >= max_pool_size_mb @@ -602,20 +584,13 @@ function collect_gmi_FisSal2011( end end - elapsed_time = time() - initial_time - if elapsed_time > time_limit - @info "Time limit exceeded. Stopping." - break - end - if round == 1 @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", "obj", "cl_curr", "cl_best", - "active", "pool_cuts", "pool_mb", "bases", @@ -626,55 +601,30 @@ function collect_gmi_FisSal2011( ) end - if time() - last_print_time > interval_print + if time() - last_print_time > print_interval_sec last_print_time = time() @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, obj_curr, gapcl_curr, gapcl_best, - length(selected_idx), length(pool.ub), pool_size_mb, length(basis_seen), - backtrack_count, + count_backtrack, Δ, μ, λ, ) 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 + elapsed_time = time() - initial_time + if elapsed_time > time_limit + @info "Time limit exceeded. Stopping." + break + end + end @info "Best gap closure: $(gapcl_best)" @@ -689,7 +639,6 @@ function collect_gmi_FisSal2011( "obj_final" => obj_curr, "obj_initial" => obj_initial, "obj_mip" => obj_mip, - "pool_active" => pool_active, "pool_size_mb" => pool_size_mb, "pool_total" => length(pool.lb), "time" => stats_time,