FisSal2011: Implement miplearn variant; minor fixes

src/Cuts/tableau/gmi.jl

@@ -306,6 +306,7 @@ function compute_gmi(data::ProblemData, tableau::Tableau)::ConstraintSet
         resize!(cut_lhs_V, nnz_count)
     end
 
+    # TODO: Build cut in compressed row format instead of converting
     @timeit "Convert to ConstraintSet" begin
         cut_lhs::SparseMatrixCSC = sparse(cut_lhs_I, cut_lhs_J, cut_lhs_V, nrows, ncols)
         cs::ConstraintSet = ConstraintSet(; lhs=cut_lhs, ub=cut_ub, lb=cut_lb, hash=cut_hash)

src/Cuts/tableau/gmi_dual.jl

@@ -266,15 +266,15 @@ end
 
 function collect_gmi_FisSal2011(
     mps_filename;
-    interval_print_sec=0.1,
+    interval_print_sec = 1,
     max_cuts_per_round = 1_000_000,
     max_pool_size_mb = 1024,
     optimizer,
-    silent_solver=true,
+    silent_solver = true,
-    time_limit = 3_600,
+    time_limit = 14_400,
-    variant = :fast,
+    variant = :miplearn,
 )
-    variant in [:subg, :hybr, :fast, :faster] || error("unknown variant: $variant")
+    variant in [:subg, :hybr, :fast, :faster, :miplearn] || error("unknown variant: $variant")
     if variant == :subg
         max_rounds = 10_000
         interval_large_lp = 10_000
@@ -291,8 +291,12 @@ function collect_gmi_FisSal2011(
         max_rounds = 500
         interval_large_lp = 50
         interval_read_tableau = 1
+    elseif variant == :miplearn
+        max_rounds = 1_000_000
+        interval_large_lp = 100
+        interval_read_tableau = 1
     end
-    gapcl_best_patience = 2 * interval_large_lp + 1
+    gapcl_best_patience = 2 * interval_large_lp + 5
 
     reset_timer!()
     initial_time = time()
@@ -407,6 +411,7 @@ function collect_gmi_FisSal2011(
         end
 
         @timeit "Optimize LP (lagrangian)" begin
+            set_silent(model_s)
             optimize!(model_s)
             status = termination_status(model_s)
             if status != MOI.OPTIMAL
@@ -414,9 +419,10 @@ function collect_gmi_FisSal2011(
             end            
             sol_frac = get_x(model_s)
             obj_curr = objective_value(model_s)
+        end
 
+        @timeit "Update history and μ" begin
             push!(obj_hist, obj_curr)
-
             if obj_best === nothing || obj_curr > obj_best
                 log_prefix = '*'
                 obj_best = obj_curr
@@ -428,38 +434,35 @@ function collect_gmi_FisSal2011(
             gapcl_curr = gapcl(obj_curr)
             gapcl_best = gapcl(obj_best)
             push!(gapcl_best_history, gapcl_best)
-
+            if variant in [:subg, :hybr]
-            @timeit "Update μ" begin
+                Δ = obj_mip - obj_best
-                if variant in [:subg, :hybr]
+                if obj_curr < obj_best - Δ
-                    Δ = obj_mip - obj_best
+                    count_deterioration += 1
-                    if obj_curr < obj_best - Δ
+                else
-                        count_deterioration += 1
+                    count_deterioration = 0
-                    else
+                end
-                        count_deterioration = 0
+                if count_deterioration >= 10
-                    end
+                    μ *= 0.5
-                    if count_deterioration >= 10
+                    multipliers_curr = multipliers_best
-                        μ *= 0.5
+                    count_deterioration = 0
-                        multipliers_curr = multipliers_best
+                    count_backtrack += 1
-                        count_deterioration = 0
+                elseif length(obj_hist) >= 100
-                        count_backtrack += 1
+                    obj_hist_avg = mean(obj_hist)
-                    elseif length(obj_hist) >= 100
+                    improv = obj_best - obj_hist[1]
-                        obj_hist_avg = mean(obj_hist)
+                    if improv < 0.01 * Δ
-                        improv = obj_best - obj_hist[1]
+                        if obj_best - obj_hist_avg < 0.001 * Δ
-                        if improv < 0.01 * Δ
+                            μ = 10 * μ
-                            if obj_best - obj_hist_avg < 0.001 * Δ
+                        elseif obj_best - obj_hist_avg < 0.01 * Δ
-                                μ = 10 * μ
+                            μ = 2 * μ
-                            elseif obj_best - obj_hist_avg < 0.01 * Δ
+                        else
-                                μ = 2 * μ
+                            μ = 0.5 * μ
-                            else
-                                μ = 0.5 * μ
-                            end
                         end
                     end
-                elseif variant in [:fast, :faster]
-                    μ = 0.01
-                else
-                    error("not implemented")
                 end
+            elseif variant in [:fast, :faster, :miplearn]
+                μ = 0.01
+            else
+                error("not implemented")
             end
         end
 
@@ -479,6 +482,9 @@ function collect_gmi_FisSal2011(
 
             @timeit "Compute GMI cuts" begin
                 cuts_s = compute_gmi(data_s, tableau)
+            end
+
+            @timeit "Check cut validity" begin
                 assert_cuts_off(cuts_s, sol_frac)
                 assert_does_not_cut_off(cuts_s, sol_opt_s)
                 ncuts = length(cuts_s.lb)
@@ -495,6 +501,7 @@ function collect_gmi_FisSal2011(
                         end
                     end
                 end
+                # TODO: Reduce allocations and improve performance
                 @timeit "Append unique cuts" begin
                     if round == 1
                         pool = ConstraintSet(
@@ -559,6 +566,7 @@ function collect_gmi_FisSal2011(
                 selected_contrs = []
                 while true
                     @timeit "Optimize LP (extended)" begin
+                        set_silent(model_s)
                         set_objective_function(model_s, orig_obj_s)
                         optimize!(model_s)
                         status = termination_status(model_s)
@@ -620,7 +628,6 @@ function collect_gmi_FisSal2011(
                     end
                     gapcl_curr = gapcl(obj_curr)
                     gapcl_best = gapcl(obj_best)
-                    push!(gapcl_best_history, gapcl_best)
                 end
 
                 @timeit "Delete all cut constraints" begin
@@ -680,6 +687,7 @@ function collect_gmi_FisSal2011(
             )
         end
         
+        push!(gapcl_best_history, gapcl_best)
         if length(gapcl_best_history) >= gapcl_best_patience
             if gapcl_best <= gapcl_best_history[1]
                 @info "No gap closure improvement. Stopping."