Make GMI cuts more stable

src/Cuts/Cuts.jl

@@ -9,6 +9,7 @@ import ..to_str_array
 include("tableau/structs.jl")
 
 # include("blackbox/cplex.jl")
+include("tableau/numerics.jl")
 include("tableau/collect.jl")
 include("tableau/gmi.jl")
 include("tableau/moi.jl")

src/Cuts/tableau/collect.jl

@@ -5,8 +5,10 @@
 import ..H5File
 
 using OrderedCollections
+using Statistics
 
-function collect_gmi(mps_filename; optimizer, max_rounds = 10, max_cuts_per_round = 100)
+
+function collect_gmi(mps_filename; optimizer, max_rounds=10, max_cuts_per_round=100, atol=1e-4)
     @info mps_filename
     reset_timer!()
 
@@ -27,7 +29,7 @@ function collect_gmi(mps_filename; optimizer, max_rounds = 10, max_cuts_per_roun
     if obj_mip === nothing
         obj_mip = h5.get_scalar("mip_obj_value")
     end
-    obj_lp = nothing
+    obj_lp = h5.get_scalar("lp_obj_value")
     h5.file.close()
 
     # Define relative MIP gap
@@ -58,8 +60,8 @@ function collect_gmi(mps_filename; optimizer, max_rounds = 10, max_cuts_per_roun
             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)
+            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)
@@ -71,15 +73,17 @@ function collect_gmi(mps_filename; optimizer, max_rounds = 10, max_cuts_per_roun
             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
+            assert_eq(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
+
+        if round == 1
+            # Assert standard form problem has same value as original
+            assert_eq(obj, obj_lp)
             push!(stats_obj, obj)
             push!(stats_gap, gap(obj))
             push!(stats_ncuts, 0)
@@ -93,16 +97,16 @@ function collect_gmi(mps_filename; optimizer, max_rounds = 10, max_cuts_per_roun
         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)
+                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)
+            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
+            assert_eq(tableau.lhs * sol_frac, tableau.rhs)
+            assert_eq(tableau.lhs * sol_opt_s, tableau.rhs)
         end
 
         # Compute GMI cuts
@@ -110,17 +114,12 @@ function collect_gmi(mps_filename; optimizer, max_rounds = 10, max_cuts_per_roun
             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
+            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."
+                @info "No cuts generated. Stopping."
                 break
             end
         end
@@ -139,7 +138,7 @@ function collect_gmi(mps_filename; optimizer, max_rounds = 10, max_cuts_per_roun
         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]
+        useful = [abs(shadow_price(c)) > atol for c in constrs]
         drop = findall(useful .== false)
         keep = findall(useful .== true)
         delete.(model, constrs[drop])
@@ -174,7 +173,6 @@ function collect_gmi(mps_filename; optimizer, max_rounds = 10, max_cuts_per_roun
         "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,

src/Cuts/tableau/gmi.jl

@@ -5,13 +5,14 @@
 using SparseArrays
 using TimerOutputs
 
-@inline frac(x::Float64) = x - floor(x)
-
-function select_gmi_rows(data, basis, x; max_rows = 10, atol = 0.001)
+function select_gmi_rows(data, basis, x; max_rows=10, atol=1e-4)
     candidate_rows = [
         r for
-        r = 1:length(basis.var_basic) if (data.var_types[basis.var_basic[r]] != 'C') &&
-        (frac(x[basis.var_basic[r]]) > atol)
+        r in 1:length(basis.var_basic) if (
+            (data.var_types[basis.var_basic[r]] != 'C') &&
+            (frac(x[basis.var_basic[r]]) > atol) &&
+            (frac2(x[basis.var_basic[r]]) > atol)
+        )
     ]
     candidate_vals = frac.(x[basis.var_basic[candidate_rows]])
     score = abs.(candidate_vals .- 0.5)
@@ -19,34 +20,36 @@ function select_gmi_rows(data, basis, x; max_rows = 10, atol = 0.001)
     return [candidate_rows[perm[i]] for i = 1:min(length(perm), max_rows)]
 end
 
-function compute_gmi(data::ProblemData, tableau::Tableau, tol = 1e-8)::ConstraintSet
+function compute_gmi(data::ProblemData, tableau::Tableau)::ConstraintSet
     nrows, ncols = size(tableau.lhs)
     ub = Float64[Inf for _ = 1:nrows]
-    lb = Float64[0.999 for _ = 1:nrows]
+    lb = Float64[0.9999 for _ = 1:nrows]
     tableau_I, tableau_J, tableau_V = findnz(tableau.lhs)
     lhs_I = Int[]
     lhs_J = Int[]
     lhs_V = Float64[]
     @timeit "Compute coefficients" begin
-        for k = 1:nnz(tableau.lhs)
+        for k in 1:nnz(tableau.lhs)
             i::Int = tableau_I[k]
+            j::Int = tableau_J[k]
             v::Float64 = 0.0
-            alpha_j = frac(tableau_V[k])
+            frac_alpha_j = frac(tableau_V[k])
+            alpha_j = tableau_V[k]
             beta = frac(tableau.rhs[i])
-            if data.var_types[i] == "C"
+            if data.var_types[j] == 'C'
                 if alpha_j >= 0
                     v = alpha_j / beta
                 else
-                    v = alpha_j / (1 - beta)
+                    v = -alpha_j / (1 - beta)
                 end
             else
-                if alpha_j <= beta
-                    v = alpha_j / beta
+                if frac_alpha_j < beta
+                    v = frac_alpha_j / beta
                 else
-                    v = (1 - alpha_j) / (1 - beta)
+                    v = (1 - frac_alpha_j) / (1 - beta)
                 end
             end
-            if abs(v) > tol
+            if abs(v) > 1e-8
                 push!(lhs_I, i)
                 push!(lhs_J, tableau_J[k])
                 push!(lhs_V, v)
@@ -57,28 +60,4 @@ function compute_gmi(data::ProblemData, tableau::Tableau, tol = 1e-8)::Constrain
     return ConstraintSet(; lhs, ub, lb)
 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
-
 export compute_gmi, frac, select_gmi_rows, assert_cuts_off, assert_does_not_cut_off

src/Cuts/tableau/numerics.jl

@@ -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