FisSal2011: Implement early termination; improve log

src/Cuts/tableau/gmi_dual.jl

@@ -266,13 +266,13 @@ end
 
 function collect_gmi_FisSal2011(
     mps_filename;
-    optimizer,
-    max_cuts_per_round = 1_000_000,
-    time_limit = 3_600,
     interval_print_sec=0.1,
-    silent_solver=true,
+    max_cuts_per_round = 1_000_000,
     max_pool_size_mb = 1024,
-    variant = :fast,
+    optimizer,
+    silent_solver=true,
+    time_limit = 900,
+    variant = :faster,
 )
     variant in [:subg, :hybr, :fast, :faster] || error("unknown variant: $variant")
     if variant == :subg
@@ -292,6 +292,7 @@ function collect_gmi_FisSal2011(
         interval_large_lp = 50
         interval_read_tableau = 1
     end
+    gapcl_best_patience = 2 * interval_large_lp + 1
 
     reset_timer!()
     initial_time = time()
@@ -314,13 +315,14 @@ function collect_gmi_FisSal2011(
         count_backtrack = 0
         count_deterioration = 0
         gapcl_best = 0
+        gapcl_best_history = CircularBuffer{Float64}(gapcl_best_patience)
         gapcl_curr = 0
         last_print_time = 0
         multipliers_best = nothing
         multipliers_curr = nothing
         obj_best = nothing
         obj_curr = nothing
-        obj_hist = []
+        obj_hist = CircularBuffer{Float64}(100)
         obj_initial = nothing
         pool = nothing
         pool_cut_age = nothing
@@ -347,6 +349,11 @@ function collect_gmi_FisSal2011(
         # 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)
+        for (var_idx, var_type) in enumerate(data.var_types)
+            if var_type in ['B', 'I']
+                assert_int(sol_opt[var_idx])
+            end
+        end
 
         # Convert to standard form
         data_s, transforms = convert_to_standard_form(data)
@@ -363,13 +370,19 @@ function collect_gmi_FisSal2011(
         sol_opt_s = forward(transforms, sol_opt)
 
         # Assert converted solution is feasible for standard form problem
+        for (var_idx, var_type) in enumerate(data_s.var_types)
+            if var_type in ['B', 'I']
+                assert_int(sol_opt_s[var_idx])
+            end
+        end
         assert_eq(data_s.constr_lhs * sol_opt_s, data_s.constr_lb)
     end
 
     @info "Standard form model has $(length(data.var_lb)) vars, $(length(data.constr_lb)) constrs"
 
     for round = 1:max_rounds
-        should_print_log = false
+        log_prefix = ' '
+        log_should_print = false
 
         if round > 1
             @timeit "Update objective function" begin
@@ -395,15 +408,17 @@ function collect_gmi_FisSal2011(
 
         @timeit "Optimize LP (lagrangian)" begin
             optimize!(model_s)
+            status = termination_status(model_s)
+            if status != MOI.OPTIMAL
+                error("Non-optimal termination status: $status")
+            end            
             sol_frac = get_x(model_s)
             obj_curr = objective_value(model_s)
             
             push!(obj_hist, obj_curr)
-            if length(obj_hist) > 100
-                popfirst!(obj_hist)
-            end
 
             if obj_best === nothing || obj_curr > obj_best
+                log_prefix = '*'
                 obj_best = obj_curr
                 multipliers_best = multipliers_curr
             end
@@ -412,9 +427,7 @@ function collect_gmi_FisSal2011(
             end
             gapcl_curr = gapcl(obj_curr)
             gapcl_best = gapcl(obj_best)
-            if termination_status(model_s) != MOI.OPTIMAL
+            push!(gapcl_best_history, gapcl_best)
-                error("Non-optimal termination status")
-            end
 
             @timeit "Update μ" begin
                 if variant in [:subg, :hybr]
@@ -508,10 +521,38 @@ function collect_gmi_FisSal2011(
                     end
                 end
             end
+
+            @timeit "Prune the pool" begin
+                pool_size_mb = Base.summarysize(pool) / 1024^2
+                while pool_size_mb >= max_pool_size_mb
+                    @timeit "Identify cuts to remove" begin
+                        σ = sortperm(pool_cut_age, rev=true)
+                        pool_size = length(pool.ub)
+                        n_keep = Int(floor(pool_size * 0.8))
+                        idx_keep = σ[1:n_keep]
+                        idx_remove = σ[(n_keep+1):end]
+                    end
+                    @timeit "Update cut hashes" begin
+                        for idx in idx_remove
+                            cut_data = (pool.lhs[:, idx], pool.lb[idx], pool.ub[idx])
+                            delete!(pool_cut_hashes, hash(cut_data))
+                        end
+                    end
+                    @timeit "Update cut pool" begin
+                        pool.lhs = pool.lhs[:, idx_keep]
+                        pool.ub = pool.ub[idx_keep]
+                        pool.lb = pool.lb[idx_keep]
+                        multipliers_curr = multipliers_curr[idx_keep]
+                        multipliers_best = multipliers_best[idx_keep]
+                        pool_cut_age = pool_cut_age[idx_keep]
+                    end
+                    pool_size_mb = Base.summarysize(pool) / 1024^2
+                end
+            end
         end
         
         if mod(round - 1, interval_large_lp) == 0 || round == max_rounds
-            should_print_log = true
+            log_should_print = true
             @timeit "Update multipliers (large LP)" begin
                 selected_idx = []
                 selected_contrs = []
@@ -519,6 +560,10 @@ function collect_gmi_FisSal2011(
                     @timeit "Optimize LP (extended)" begin
                         set_objective_function(model_s, orig_obj_s)
                         optimize!(model_s)
+                        status = termination_status(model_s)
+                        if status != MOI.OPTIMAL
+                            error("Non-optimal termination status: $status")
+                        end
                         obj_curr = objective_value(model_s)
                         sol_frac = get_x(model_s)
                     end
@@ -568,38 +613,13 @@ function collect_gmi_FisSal2011(
 
                 @timeit "Update best" begin
                     if obj_curr > obj_best
+                        log_prefix = '*'
                         obj_best = obj_curr
                         multipliers_best = multipliers_curr
                     end
                     gapcl_curr = gapcl(obj_curr)
                     gapcl_best = gapcl(obj_best)
-                end
+                    push!(gapcl_best_history, gapcl_best)
-
-                @timeit "Prune cut pool" begin
-                    pool_size_mb = Base.summarysize(pool) / 1024^2
-                    while pool_size_mb >= max_pool_size_mb
-                        @timeit "Identify cuts to remove" begin
-                            σ = sortperm(pool_cut_age, rev=true)
-                            pool_size = length(pool.ub)
-                            n_keep = Int(floor(pool_size * 0.8))
-                            idx_keep = 1:n_keep
-                            idx_remove = σ[(n_keep+1):end]
-                        end
-                        @timeit "Update cut hashes" begin
-                            for idx in idx_remove
-                                cut_data = (pool.lhs[:, idx], pool.lb[idx], pool.ub[idx])
-                                delete!(pool_cut_hashes, hash(cut_data))
-                            end
-                        end
-                        @timeit "Update cut pool" begin
-                            pool.lhs = pool.lhs[:, idx_keep]
-                            pool.ub = pool.ub[idx_keep]
-                            pool.lb = pool.lb[idx_keep]
-                            multipliers_curr = multipliers[idx_keep]
-                            pool_cut_age = pool_cut_age[idx_keep]
-                        end
-                        pool_size_mb = Base.summarysize(pool) / 1024^2
-                    end
                 end
 
                 @timeit "Delete all cut constraints" begin
@@ -609,7 +629,12 @@ function collect_gmi_FisSal2011(
         else
             @timeit "Update multipliers (subgradient)" begin
                 subgrad = (pool.lb' - (sol_frac' * pool.lhs))'
-                λ = μ * (obj_mip - obj_curr) / norm(subgrad)^2
+                subgrad_norm_sq = norm(subgrad)^2
+                if subgrad_norm_sq < 1e-10
+                    λ = 0
+                else
+                    λ = μ * (obj_mip - obj_curr) / subgrad_norm_sq
+                end
                 multipliers_curr = max.(0, multipliers_curr .+ λ * subgrad)
             end
         end
@@ -632,13 +657,14 @@ function collect_gmi_FisSal2011(
         end
 
         if time() - last_print_time > interval_print_sec
-            should_print_log = true
+            log_should_print = true
         end
 
-        if should_print_log
+        if log_should_print
             last_print_time = time()
             @printf(
-                "%8d %10.3e %9.2e %9.2e %9d %9.2f %9d %4d %8.2e %8.2e %8.2e\n",
+                "%c %8d %10.3e %9.2e %9.2e %9d %9.2f %9d %4d %8.2e %8.2e %8.2e\n",
+                log_prefix,
                 round,
                 obj_curr,
                 gapcl_curr,
@@ -653,6 +679,13 @@ function collect_gmi_FisSal2011(
             )
         end
 
+        if length(gapcl_best_history) >= gapcl_best_patience
+            if gapcl_best <= gapcl_best_history[1]
+                @info "No gap closure improvement. Stopping."
+                break
+            end
+        end
+
         elapsed_time = time() - initial_time
         if elapsed_time > time_limit
             @info "Time limit exceeded. Stopping."

src/Cuts/tableau/numerics.jl

@@ -41,6 +41,8 @@ function assert_does_not_cut_off(cuts::ConstraintSet, x::Vector{Float64}; tol =
         ub = cuts.ub[i]
         lb = cuts.lb[i]
         if (val >= ub) || (val <= lb)
+            @show cuts.lhs[i, :]
+            @show x[cuts.lhs[i, :].nzind]
             throw(
                 ErrorException(
                     "inequality $i cuts off integer solution ($lb <= $val <= $ub)",
@@ -49,3 +51,10 @@ function assert_does_not_cut_off(cuts::ConstraintSet, x::Vector{Float64}; tol =
         end
     end
 end
+
+function assert_int(x::Float64, tol=1e-6)
+    fx = frac(x)
+    if min(fx, 1 - fx) >= tol
+        throw(ErrorException("Number must be integer: $x"))
+    end
+end