Make dual GMI cuts stronger

src/Cuts/tableau/gmi.jl

@@ -185,7 +185,7 @@ function collect_gmi(
     )
 end
 
-function select_gmi_rows(data, basis, x; max_rows = 10, atol = 1e-4)
+function select_gmi_rows(data, basis, x; max_rows = 10, atol = 0.001)
     candidate_rows = [
         r for r = 1:length(basis.var_basic) if (
             (data.var_types[basis.var_basic[r]] != 'C') &&
@@ -202,7 +202,7 @@ end
 function compute_gmi(data::ProblemData, tableau::Tableau)::ConstraintSet
     nrows, ncols = size(tableau.lhs)
     ub = Float64[Inf for _ = 1:nrows]
-    lb = Float64[0.9999 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[]

src/Cuts/tableau/gmi_dual.jl

@@ -7,6 +7,7 @@ using JuMP
 using HiGHS
 using Random
 using DataStructures
+using Statistics
 
 import ..H5FieldsExtractor
 
@@ -25,7 +26,7 @@ function collect_gmi_dual(
     mps_filename;
     optimizer,
     max_rounds = 10,
-    max_cuts_per_round = 500,
+    max_cuts_per_round = 1_000_000,
     time_limit = 3_600,
 )
     reset_timer!()
@@ -263,6 +264,342 @@ function collect_gmi_dual(
     )
 end
 
+function collect_gmi_FisSal2011(
+    mps_filename;
+    optimizer,
+    max_rounds = 10_000,
+    max_cuts_per_round = 1_000_000,
+    time_limit = 30,
+    interval_print=1,
+    max_cut_age=10,
+)
+    reset_timer!()
+    initial_time = time()
+
+    @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
+    
+    @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
+        gapcl_best = 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)
+        p = (1 - ε) .+ 2 * ε * rand(length(v))
+        return v .* p
+    end
+
+    @timeit "Read problem" begin
+        model = read_from_file(mps_filename)
+        set_optimizer(model, optimizer)
+    end
+
+    @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)
+        vars_s = all_variables(model_s)
+        orig_obj_s = objective_function(model_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
+
+    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))
+                for offset in 1:nnz(v)
+                    var_idx = v.nzind[offset]
+                    add_to_expression!(
+                        lagr_term,
+                        vars_s[var_idx],
+                        - v.nzval[offset],
+                    )
+                end
+
+                # Update objective
+                set_objective_function(
+                    model_s,
+                    orig_obj_s + lagr_term,
+                )
+            end
+        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
+            end
+
+            if obj_curr >= obj_best
+                obj_best = obj_curr
+                gapcl_best = gapcl(obj_best)
+                lambda_best = lambda_curr
+                noimprov_count = 0
+            else
+                noimprov_count += 1
+            end
+
+            if noimprov_count > 10
+                lambda_curr = lambda_best
+                backtrack_count += 1
+                noimprov_count = 0
+                continue
+            end
+
+            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
+        end
+
+        @timeit "Select tableau rows" begin
+            basis = get_basis(model_s)
+            if round == 1
+                original_basis = basis
+            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)
+
+            # 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_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
+            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)]
+            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 = []
+
+        @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 "Find most violated cut" begin
+                    violations = pool.lb .- pool.lhs * sol_frac
+                    σ = sortperm(violations, rev=true)
+                end
+
+                # Stop if all cuts are satisfied
+                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]
+                    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
+                    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 "Delete all cut constraints" begin
+                delete.(model_s, selected_contrs)
+            end
+        end
+
+        push!(stats_ncuts, length(pool.lb))
+
+        elapsed_time = time() - initial_time
+        if elapsed_time > time_limit
+            @info "Time limit exceeded. Stopping."
+            break
+        end
+
+        if round == 1
+            @printf(
+                "%8s %12s %12s %12s %8s %8s %9s\n",
+                "round",
+                "obj",
+                "gapcl_curr",
+                "gapcl_best",
+                "active",
+                "pool",
+                "backtrack",
+            )
+        end
+
+        if time() - last_print_time > interval_print
+            last_print_time = time()
+            @printf(
+                "%8d %12.6e %12.2f %12.2f %8d %8d %9d\n",
+                round,
+                obj_curr,
+                gapcl(obj_curr),
+                gapcl_best,
+                length(selected_idx),
+                length(pool.ub),
+                backtrack_count,
+            )
+        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_ncuts" => stats_ncuts,
+        "stats_time" => stats_time,
+    )
+end
+
 function add_constraint_set_dual_v2(model::JuMP.Model, cs::ConstraintSet)
     vars = all_variables(model)
     nrows, ncols = size(cs.lhs)
@@ -576,5 +913,5 @@ function __init_gmi_dual__()
     copy!(ExpertDualGmiComponent, ExpertDualGmiComponentPy)
 end
 
-export collect_gmi_dual, expert_gmi_dual, ExpertDualGmiComponent, KnnDualGmiComponent
+export collect_gmi_dual, expert_gmi_dual, ExpertDualGmiComponent, KnnDualGmiComponent, collect_gmi_FisSal2011
 

src/Cuts/tableau/transform.jl

@@ -102,7 +102,14 @@ function forward!(t::AddSlackVariables, data::ProblemData)
     ge = [i for i = 1:nrows if isfinite(data.constr_lb[i]) && !isequality[i]]
     le = [i for i = 1:nrows if isfinite(data.constr_ub[i]) && !isequality[i]]
     EQ, GE, LE = length(eq), length(ge), length(le)
-
+    is_integral(row_idx, rhs) = (
+        abs(rhs - round(rhs)) <= 1e-6 && 
+        all(j -> data.var_types[j] ∈ ['I', 'B'], findnz(data.constr_lhs[row_idx, :])[1])
+    )
+    slack_types = [
+        [is_integral(ge[i], data.constr_lb[ge[i]]) ? 'I' : 'C' for i = 1:GE];
+        [is_integral(le[i], data.constr_ub[le[i]]) ? 'I' : 'C' for i = 1:LE]
+    ]
     t.M1 = [
         I spzeros(ncols, GE + LE)
         data.constr_lhs[ge, :] spzeros(GE, GE + LE)
@@ -129,7 +136,7 @@ function forward!(t::AddSlackVariables, data::ProblemData)
     data.var_lb = [data.var_lb; zeros(GE + LE)]
     data.var_ub = [data.var_ub; [Inf for _ = 1:(GE+LE)]]
     data.var_names = [data.var_names; ["__s$i" for i = 1:(GE+LE)]]
-    data.var_types = [data.var_types; ['C' for _ = 1:(GE+LE)]]
+    data.var_types = [data.var_types; slack_types]
     data.constr_lb = [
         data.constr_lb[eq]
         data.constr_lb[ge]