FisSal2011: Adjust constants

Project.toml

@@ -1,11 +1,12 @@
 name = "MIPLearn"
 uuid = "2b1277c3-b477-4c49-a15e-7ba350325c68"
 authors = ["Alinson S Xavier <git@axavier.org>"]
-version = "0.4.0"
+version = "0.4.2"
 
 [deps]
 Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+Gurobi = "2e9cd046-0924-5485-92f1-d5272153d98b"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
@@ -15,10 +16,12 @@ KLU = "ef3ab10e-7fda-4108-b977-705223b18434"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
+SCIP = "82193955-e24f-5292-bf16-6f2c5261a85f"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
@@ -26,16 +29,19 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 [compat]
 Conda = "1"
 DataStructures = "0.18"
+Gurobi = "1.7.5"
 HDF5 = "0.16"
 HiGHS = "1"
 JLD2 = "0.4"
 JSON = "0.21"
-julia = "1"
 JuMP = "1"
 KLU = "0.4"
 MathOptInterface = "1"
 OrderedCollections = "1"
+PrecompileTools = "1"
 PyCall = "1"
 Requires = "1"
+SCIP = "0.12"
 Statistics = "1"
 TimerOutputs = "0.5"
+julia = "1"

deps/build.jl

@@ -5,7 +5,7 @@ function install_miplearn()
     Conda.update()
     pip = joinpath(dirname(pyimport("sys").executable), "pip")
     isfile(pip) || error("$pip: invalid path")
-    run(`$pip install miplearn==0.4.0`)
+    run(`$pip install miplearn==0.4.4`)
 end
 
 install_miplearn()

src/BB/log.jl

@@ -6,7 +6,7 @@ using Printf
 
 function print_progress_header()
     @printf(
-        "%8s %9s %9s %13s %13s %9s %6s %13s %6s %-24s %9s %9s %6s %6s",
+        "%8s %9s %9s %13s %13s %9s %9s %13s %9s %-24s %9s %9s %6s %6s",
         "time",
         "processed",
         "pending",
@@ -46,7 +46,7 @@ function print_progress(
             branch_ub = @sprintf("%9.2f", last(node.branch_ub))
         end
         @printf(
-            "%8.2f %9d %9d %13.6e %13.6e %9.2e %6d %13.6e %6s %-24s %9s %9s %6d %6d",
+            "%8.2f %9d %9d %13.6e %13.6e %9.2e %9d %13.6e %9s %-24s %9s %9s %6d %6d",
             time_elapsed,
             pool.processed,
             length(pool.processing) + length(pool.pending),

src/BB/lp.jl

@@ -134,7 +134,11 @@ function _get_int_variables(
                 var_ub = constr.upper
                 MOI.delete(optimizer, _upper_bound_index(var))
             end
-            MOI.add_constraint(optimizer, var, MOI.Interval(var_lb, var_ub))
+            MOI.add_constraint(
+                optimizer,
+                MOI.VariableIndex(var.index),
+                MOI.Interval(var_lb, var_ub),
+            )
         end
         push!(vars, var)
         push!(lb, var_lb)

src/Cuts/Cuts.jl

@@ -4,15 +4,22 @@
 
 module Cuts
 
+using PyCall
+
 import ..to_str_array
 
 include("tableau/structs.jl")
 
 # include("blackbox/cplex.jl")
-include("tableau/collect.jl")
+include("tableau/numerics.jl")
 include("tableau/gmi.jl")
+include("tableau/gmi_dual.jl")
 include("tableau/moi.jl")
 include("tableau/tableau.jl")
 include("tableau/transform.jl")
 
+function __init__()
+    __init_gmi_dual__()
+end
+
 end # module

src/Cuts/tableau/collect.jl

@@ -1,184 +0,0 @@
-#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
-#  Copyright (C) 2020-2023, UChicago Argonne, LLC. All rights reserved.
-#  Released under the modified BSD license. See COPYING.md for more details.
-
-import ..H5File
-
-using OrderedCollections
-
-function collect_gmi(mps_filename; optimizer, max_rounds = 10, max_cuts_per_round = 100)
-    @info mps_filename
-    reset_timer!()
-
-    # Open HDF5 file
-    h5_filename = replace(mps_filename, ".mps.gz" => ".h5")
-    h5 = H5File(h5_filename)
-
-    # Read optimal solution
-    sol_opt_dict = Dict(
-        zip(
-            h5.get_array("static_var_names"),
-            convert(Array{Float64}, h5.get_array("mip_var_values")),
-        ),
-    )
-
-    # Read optimal value
-    obj_mip = h5.get_scalar("mip_lower_bound")
-    if obj_mip === nothing
-        obj_mip = h5.get_scalar("mip_obj_value")
-    end
-    obj_lp = nothing
-    h5.file.close()
-
-    # Define relative MIP gap
-    gap(v) = 100 * abs(obj_mip - v) / abs(v)
-
-    # Initialize stats
-    stats_obj = []
-    stats_gap = []
-    stats_ncuts = []
-    stats_time_convert = 0
-    stats_time_solve = 0
-    stats_time_select = 0
-    stats_time_tableau = 0
-    stats_time_gmi = 0
-    all_cuts = nothing
-
-    # Read problem
-    model = read_from_file(mps_filename)
-
-    for round = 1:max_rounds
-        @info "Round $(round)..."
-
-        stats_time_convert = @elapsed 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 all(data.constr_lb .- 1e-3 .<= data.constr_lhs * sol_opt)
-            @assert all(data.constr_lhs * sol_opt .<= data.constr_ub .+ 1e-3)
-
-            # Convert to standard form
-            data_s, transforms = convert_to_standard_form(data)
-            model_s = to_model(data_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 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
-            push!(stats_obj, obj)
-            push!(stats_gap, gap(obj))
-            push!(stats_ncuts, 0)
-        end
-        if termination_status(model_s) != MOI.OPTIMAL
-            return
-        end
-
-        # Select tableau rows
-        basis = get_basis(model_s)
-        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)
-        end
-
-        # Compute selected tableau rows
-        stats_time_tableau += @elapsed begin
-            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
-        end
-
-        # Compute GMI cuts
-        stats_time_gmi += @elapsed begin
-            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
-
-            # Abort if no cuts are left
-            if length(cuts_s.lb) == 0
-                @info "No cuts generated. Aborting."
-                break
-            end
-        end
-
-        # Add GMI cuts to original problem
-        cuts = backwards(transforms, cuts_s)
-        assert_does_not_cut_off(cuts, sol_opt)
-        constrs = add_constraint_set(model, cuts)
-
-        # Optimize original form
-        set_optimizer(model, optimizer)
-        undo_relax = relax_integrality(model)
-        optimize!(model)
-        obj = objective_value(model)
-        push!(stats_obj, obj)
-        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]
-        drop = findall(useful .== false)
-        keep = findall(useful .== true)
-        delete.(model, constrs[drop])
-        if all_cuts === nothing
-            all_cuts = cuts
-        else
-            all_cuts.lhs = [all_cuts.lhs; cuts.lhs[keep, :]]
-            all_cuts.lb = [all_cuts.lb; cuts.lb[keep]]
-            all_cuts.lb = [all_cuts.lb; cuts.lb[keep]]
-        end
-        push!(stats_ncuts, length(all_cuts.lb))
-
-        undo_relax()
-    end
-
-    # Store cuts
-    if all_cuts !== nothing
-        @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.file.close()
-    end
-
-    return OrderedDict(
-        "instance" => mps_filename,
-        "max_rounds" => max_rounds,
-        "rounds" => length(stats_obj) - 1,
-        "time_convert" => stats_time_convert,
-        "time_solve" => stats_time_solve,
-        "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,
-    )
-end
-
-export collect_gmi

src/Cuts/tableau/gmi.jl

@@ -2,16 +2,196 @@
 #  Copyright (C) 2020-2023, UChicago Argonne, LLC. All rights reserved.
 #  Released under the modified BSD license. See COPYING.md for more details.
 
+import ..H5File
+
+using OrderedCollections
 using SparseArrays
+using Statistics
 using TimerOutputs
 
-@inline frac(x::Float64) = x - floor(x)
+function collect_gmi(
+    mps_filename;
+    optimizer,
+    max_rounds = 10,
+    max_cuts_per_round = 100,
+    atol = 1e-4,
+)
+    reset_timer!()
+
+    # Open HDF5 file
+    h5_filename = replace(mps_filename, ".mps.gz" => ".h5")
+    h5 = H5File(h5_filename)
+
+    # Read optimal solution
+    sol_opt_dict = Dict(
+        zip(
+            h5.get_array("static_var_names"),
+            convert(Array{Float64}, h5.get_array("mip_var_values")),
+        ),
+    )
+
+    # Read optimal value
+    obj_mip = h5.get_scalar("mip_lower_bound")
+    if obj_mip === nothing
+        obj_mip = h5.get_scalar("mip_obj_value")
+    end
+    obj_lp = h5.get_scalar("lp_obj_value")
+    h5.file.close()
+
+    # Define relative MIP gap
+    gap(v) = 100 * abs(obj_mip - v) / abs(v)
+
+    # Initialize stats
+    stats_obj = []
+    stats_gap = []
+    stats_ncuts = []
+    stats_time_convert = 0
+    stats_time_solve = 0
+    stats_time_select = 0
+    stats_time_tableau = 0
+    stats_time_gmi = 0
+    all_cuts = nothing
+
+    # Read problem
+    model = read_from_file(mps_filename)
+
+    for round = 1:max_rounds
+        @info "Round $(round)..."
+
+        stats_time_convert = @elapsed 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)
+            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
+
+        # Optimize standard form
+        optimize!(model_s)
+        stats_time_solve += solve_time(model_s)
+        obj = objective_value(model_s)
+
+        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)
+        end
+        if termination_status(model_s) != MOI.OPTIMAL
+            return
+        end
+
+        # Select tableau rows
+        basis = get_basis(model_s)
+        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)
+        end
+
+        # Compute selected tableau rows
+        stats_time_tableau += @elapsed begin
+            tableau = compute_tableau(data_s, basis, sol_frac, rows = selected_rows)
+
+            # Assert tableau rows have been computed correctly
+            assert_eq(tableau.lhs * sol_frac, tableau.rhs)
+            assert_eq(tableau.lhs * sol_opt_s, tableau.rhs)
+        end
+
+        # Compute GMI cuts
+        stats_time_gmi += @elapsed begin
+            cuts_s = compute_gmi(data_s, tableau)
+
+            # Assert cuts have been generated correctly
+            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. Stopping."
+                break
+            end
+        end
+
+        # Add GMI cuts to original problem
+        cuts = backwards(transforms, cuts_s)
+        assert_does_not_cut_off(cuts, sol_opt)
+        constrs = add_constraint_set(model, cuts)
+
+        # Optimize original form
+        set_optimizer(model, optimizer)
+        undo_relax = relax_integrality(model)
+        optimize!(model)
+        obj = objective_value(model)
+        push!(stats_obj, obj)
+        push!(stats_gap, gap(obj))
+
+        # Store useful cuts; drop useless ones from the problem
+        useful = [abs(shadow_price(c)) > atol for c in constrs]
+        drop = findall(useful .== false)
+        keep = findall(useful .== true)
+        delete.(model, constrs[drop])
+        if all_cuts === nothing
+            all_cuts = cuts
+        else
+            all_cuts.lhs = [all_cuts.lhs; cuts.lhs[keep, :]]
+            all_cuts.lb = [all_cuts.lb; cuts.lb[keep]]
+            all_cuts.ub = [all_cuts.ub; cuts.ub[keep]]
+        end
+        push!(stats_ncuts, length(all_cuts.lb))
+
+        undo_relax()
+    end
+
+    # Store cuts
+    if all_cuts !== nothing
+        @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.file.close()
+    end
+
+    return OrderedDict(
+        "instance" => mps_filename,
+        "max_rounds" => max_rounds,
+        "rounds" => length(stats_obj) - 1,
+        "time_convert" => stats_time_convert,
+        "time_solve" => stats_time_solve,
+        "time_tableau" => stats_time_tableau,
+        "time_gmi" => stats_time_gmi,
+        "obj_mip" => obj_mip,
+        "stats_obj" => stats_obj,
+        "stats_gap" => stats_gap,
+        "stats_ncuts" => stats_ncuts,
+    )
+end
 
 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') &&
-        (frac(x[basis.var_basic[r]]) > atol)
+        r for r = 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,66 +199,122 @@ 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
-    nrows, ncols = size(tableau.lhs)
-    ub = Float64[Inf 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[]
-    lhs_V = Float64[]
+# function compute_gmi(data::ProblemData, tableau::Tableau)::ConstraintSet
+#     @timeit "Initialization" begin
+#         nrows, ncols = size(tableau.lhs)
+#         ub = Float64[Inf 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[]
+#         lhs_V = Float64[]
+#     end
+#     @timeit "Compute coefficients" begin
+#         for k = 1:nnz(tableau.lhs)
+#             i::Int = tableau_I[k]
+#             j::Int = tableau_J[k]
+#             v::Float64 = 0.0
+#             frac_alpha_j = frac(tableau_V[k])
+#             alpha_j = tableau_V[k]
+#             beta = frac(tableau.rhs[i])
+#             if data.var_types[j] == 'C'
+#                 if alpha_j >= 0
+#                     v = alpha_j / beta
+#                 else
+#                     v = -alpha_j / (1 - beta)
+#                 end
+#             else
+#                 if frac_alpha_j < beta
+#                     v = frac_alpha_j / beta
+#                 else
+#                     v = (1 - frac_alpha_j) / (1 - beta)
+#                 end
+#             end
+#             if abs(v) > 1e-8
+#                 push!(lhs_I, i)
+#                 push!(lhs_J, tableau_J[k])
+#                 push!(lhs_V, v)
+#             end
+#         end
+#     end
+#     @timeit "Convert to ConstraintSet" begin
+#         lhs = sparse(lhs_I, lhs_J, lhs_V, nrows, ncols)
+#         cs = ConstraintSet(; lhs, ub, lb)
+#     end
+#     return cs
+# end
+
+function compute_gmi(data::ProblemData, tableau::Tableau)::ConstraintSet
+    @timeit "Initialization" begin
+        nrows::Int, ncols::Int = size(tableau.lhs)
+        var_types::Vector{Char} = data.var_types
+        tableau_rhs::Vector{Float64} = tableau.rhs
+        tableau_I::Vector{Int}, tableau_J::Vector{Int}, tableau_V::Vector{Float64} = findnz(tableau.lhs)
+    end
+
+    @timeit "Pre-allocation" begin
+        cut_ub::Vector{Float64} = fill(Inf, nrows)
+        cut_lb::Vector{Float64} = fill(0.999, nrows)
+        nnz_tableau::Int = length(tableau_I)
+        cut_lhs_I::Vector{Int} = Vector{Int}(undef, nnz_tableau)
+        cut_lhs_J::Vector{Int} = Vector{Int}(undef, nnz_tableau)
+        cut_lhs_V::Vector{Float64} = Vector{Float64}(undef, nnz_tableau)
+        cut_hash::Vector{UInt64} = zeros(UInt64, nrows)
+        nnz_count::Int = 0
+    end
+    
     @timeit "Compute coefficients" begin
-        for k = 1:nnz(tableau.lhs)
+        @inbounds for k = 1:nnz_tableau
             i::Int = tableau_I[k]
-            v::Float64 = 0.0
-            alpha_j = frac(tableau_V[k])
-            beta = frac(tableau.rhs[i])
-            if data.var_types[i] == "C"
+            j::Int = tableau_J[k]
+            alpha_j::Float64 = tableau_V[k]
+            frac_alpha_j::Float64 = alpha_j - floor(alpha_j)
+            beta_i::Float64 = tableau_rhs[i]
+            beta::Float64 = beta_i - floor(beta_i)
+            v::Float64 = 0
+            
+            # Compute coefficient
+            if 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
-                push!(lhs_I, i)
-                push!(lhs_J, tableau_J[k])
-                push!(lhs_V, v)
+            
+            # Store if significant
+            if abs(v) > 1e-8
+                nnz_count += 1
+                cut_lhs_I[nnz_count] = i
+                cut_lhs_J[nnz_count] = j
+                cut_lhs_V[nnz_count] = v
+                cut_hash[i] = hash(j, cut_hash[i])
+                cut_hash[i] = hash(v, cut_hash[i])
             end
         end
-        lhs = sparse(lhs_I, lhs_J, lhs_V, nrows, ncols)
     end
-    return ConstraintSet(; lhs, ub, lb)
+
+    @timeit "Resize arrays to actual size" begin
+        resize!(cut_lhs_I, nnz_count)
+        resize!(cut_lhs_J, nnz_count)
+        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)
+    end
+    
+    return cs
 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
+export compute_gmi,
+    frac, select_gmi_rows, assert_cuts_off, assert_does_not_cut_off, collect_gmi

src/Cuts/tableau/moi.jl

@@ -9,6 +9,8 @@ function ProblemData(model::Model)::ProblemData
 
     # Objective function
     obj = objective_function(model)
+    obj_offset = obj.constant
+    obj_sense = objective_sense(model)
     obj = [v ∈ keys(obj.terms) ? obj.terms[v] : 0.0 for v in vars]
 
     # Variable types, lower bounds and upper bounds
@@ -86,8 +88,9 @@ function ProblemData(model::Model)::ProblemData
     @assert length(constr_ub) == m
 
     return ProblemData(
-        obj_offset = 0.0;
         obj,
+        obj_offset,
+        obj_sense,
         constr_lb,
         constr_ub,
         constr_lhs,
@@ -102,6 +105,7 @@ function to_model(data::ProblemData, tol = 1e-6)::Model
     model = Model()
 
     # Variables
+    obj_expr = AffExpr(data.obj_offset)
     nvars = length(data.obj)
     @variable(model, x[1:nvars])
     for i = 1:nvars
@@ -117,8 +121,9 @@ function to_model(data::ProblemData, tol = 1e-6)::Model
         if isfinite(data.var_ub[i])
             set_upper_bound(x[i], data.var_ub[i])
         end
-        set_objective_coefficient(model, x[i], data.obj[i])
+        add_to_expression!(obj_expr, x[i], data.obj[i])
     end
+    @objective(model, data.obj_sense, obj_expr)
 
     # Constraints
     lhs = data.constr_lhs * x
@@ -140,19 +145,9 @@ function to_model(data::ProblemData, tol = 1e-6)::Model
 end
 
 function add_constraint_set(model::JuMP.Model, cs::ConstraintSet)
-    vars = all_variables(model)
-    nrows, _ = size(cs.lhs)
-    constrs = []
-    for i = 1:nrows
-        c = nothing
-        if isinf(cs.ub[i])
-            c = @constraint(model, cs.lb[i] <= dot(cs.lhs[i, :], vars))
-        elseif isinf(cs.lb[i])
-            c = @constraint(model, dot(cs.lhs[i, :], vars) <= cs.ub[i])
-        else
-            c = @constraint(model, cs.lb[i] <= dot(cs.lhs[i, :], vars) <= cs.ub[i])
-        end
-        push!(constrs, c)
+    constrs = build_constraints(model, cs)
+    for c in constrs
+        add_constraint(model, c)
     end
     return constrs
 end
@@ -164,4 +159,30 @@ function set_warm_start(model::JuMP.Model, x::Vector{Float64})
     end
 end
 
-export to_model, ProblemData, add_constraint_set, set_warm_start
+function build_constraints(model::JuMP.Model, cs::ConstraintSet)
+    vars = all_variables(model)
+    nrows, _ = size(cs.lhs)
+    constrs = []
+    for i = 1:nrows
+        # Build LHS expression
+        row = cs.lhs[i, :]
+        lhs_expr = AffExpr()
+        for (offset, val) in enumerate(row.nzval)
+            add_to_expression!(lhs_expr, vars[row.nzind[offset]], val)
+        end
+
+        # Build JuMP constraint
+        c = nothing
+        if isinf(cs.ub[i])
+            c = @build_constraint(cs.lb[i] <= lhs_expr)
+        elseif isinf(cs.lb[i])
+            c = @build_constraint(lhs_expr <= cs.ub[i])
+        else
+            c = @build_constraint(cs.lb[i] <= lhs_expr <= cs.ub[i])
+        end
+        push!(constrs, c)
+    end
+    return constrs
+end
+
+export to_model, ProblemData, add_constraint_set, set_warm_start, build_constraints

src/Cuts/tableau/numerics.jl

@@ -0,0 +1,52 @@
+@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)
+    vals = cuts.lhs * x
+    for i = 1:length(cuts.lb)
+        if (vals[i] <= cuts.ub[i] - tol) && (vals[i] >= cuts.lb[i] + tol)
+            throw(ErrorException("inequality $i fails to cut off fractional solution: $(cuts.lb[i]) <= $(vals[i]) <= $(cuts.ub[i])"))
+        end
+    end
+end
+
+function assert_does_not_cut_off(cuts::ConstraintSet, x::Vector{Float64}; tol = 1e-6)
+    vals = cuts.lhs * x
+    for i = 1:length(cuts.lb)
+        if (vals[i] >= cuts.ub[i]) || (vals[i] <= cuts.lb[i])
+            throw(ErrorException("inequality $i cuts off integer solution: $(cuts.lb[i]) <= $(vals[i]) <= $(cuts.ub[i])"))
+        end
+    end
+end
+
+function assert_int(x::Float64, tol=1e-5)
+    fx = frac(x)
+    if min(fx, 1 - fx) >= tol
+        throw(ErrorException("Number must be integer: $x"))
+    end
+end

src/Cuts/tableau/structs.jl

@@ -7,6 +7,7 @@ using SparseArrays
 Base.@kwdef mutable struct ProblemData
     obj::Vector{Float64}
     obj_offset::Float64
+    obj_sense::Any
     constr_lb::Vector{Float64}
     constr_ub::Vector{Float64}
     constr_lhs::SparseMatrixCSC
@@ -17,10 +18,10 @@ Base.@kwdef mutable struct ProblemData
 end
 
 Base.@kwdef mutable struct Tableau
-    obj::Any
-    lhs::Any
-    rhs::Any
-    z::Any
+    obj::Vector{Float64}
+    lhs::SparseMatrixCSC
+    rhs::Vector{Float64}
+    z::Float64
 end
 
 Base.@kwdef mutable struct Basis
@@ -34,6 +35,7 @@ Base.@kwdef mutable struct ConstraintSet
     lhs::SparseMatrixCSC
     ub::Vector{Float64}
     lb::Vector{Float64}
+    hash::Union{Nothing,Vector{UInt64}} = nothing
 end
 
 export ProblemData, Tableau, Basis, ConstraintSet

src/Cuts/tableau/tableau.jl

@@ -4,48 +4,161 @@
 
 using KLU
 using TimerOutputs
+using Gurobi
 
 function get_basis(model::JuMP.Model)::Basis
-    var_basic = Int[]
-    var_nonbasic = Int[]
-    constr_basic = Int[]
-    constr_nonbasic = Int[]
-
-    # Variables
-    for (i, var) in enumerate(all_variables(model))
-        bstatus = MOI.get(model, MOI.VariableBasisStatus(), var)
-        if bstatus == MOI.BASIC
-            push!(var_basic, i)
-        elseif bstatus == MOI.NONBASIC_AT_LOWER
-            push!(var_nonbasic, i)
-        else
-            error("Unknown basis status: $bstatus")
-        end
+    if isa(unsafe_backend(model), Gurobi.Optimizer)
+        return get_basis_gurobi(model)
     end
 
-    # Constraints
-    constr_index = 1
-    for (ftype, stype) in list_of_constraint_types(model)
-        for constr in all_constraints(model, ftype, stype)
-            if ftype == VariableRef
-                # nop
-            elseif ftype == AffExpr
-                bstatus = MOI.get(model, MOI.ConstraintBasisStatus(), constr)
-                if bstatus == MOI.BASIC
-                    push!(constr_basic, constr_index)
-                elseif bstatus == MOI.NONBASIC
-                    push!(constr_nonbasic, constr_index)
-                else
-                    error("Unknown basis status: $bstatus")
-                end
-                constr_index += 1
+    @timeit "Initialization" begin
+        var_basic = Int[]
+        var_nonbasic = Int[]
+        constr_basic = Int[]
+        constr_nonbasic = Int[]
+        nvars = num_variables(model)
+        sizehint!(var_basic, nvars)
+        sizehint!(var_nonbasic, nvars)
+    end
+
+    @timeit "Query variables" begin
+        for (i, var) in enumerate(all_variables(model))
+            bstatus = MOI.get(model, MOI.VariableBasisStatus(), var)
+            if bstatus == MOI.BASIC
+                push!(var_basic, i)
+            elseif bstatus == MOI.NONBASIC_AT_LOWER
+                push!(var_nonbasic, i)
             else
-                error("Unsupported constraint type: ($ftype, $stype)")
+                error("Unknown basis status: $bstatus")
             end
         end
     end
 
-    return Basis(; var_basic, var_nonbasic, constr_basic, constr_nonbasic)
+    @timeit "Query constraints" begin
+        constr_index = 1
+        for (ftype, stype) in list_of_constraint_types(model)
+            for constr in all_constraints(model, ftype, stype)
+                if ftype == VariableRef
+                    # nop
+                elseif ftype == AffExpr
+                    bstatus = MOI.get(model, MOI.ConstraintBasisStatus(), constr)
+                    if bstatus == MOI.BASIC
+                        push!(constr_basic, constr_index)
+                    elseif bstatus == MOI.NONBASIC
+                        push!(constr_nonbasic, constr_index)
+                    else
+                        error("Unknown basis status: $bstatus")
+                    end
+                    constr_index += 1
+                else
+                    error("Unsupported constraint type: ($ftype, $stype)")
+                end
+            end
+        end
+    end
+
+    @timeit "Build basis struct" begin
+        basis = Basis(; var_basic, var_nonbasic, constr_basic, constr_nonbasic)
+    end
+
+    return basis
+end
+
+function set_basis(model::JuMP.Model, basis::Basis)
+    if isa(unsafe_backend(model), Gurobi.Optimizer)
+        # NOP
+        return
+    end
+
+    @timeit "Initialization" begin
+        nvars = num_variables(model)
+        gurobi_model = unsafe_backend(model).inner
+    end
+
+    @timeit "Set variable basis" begin
+        var_basis_statuses = Vector{Cint}(undef, nvars)
+        fill!(var_basis_statuses, -1)  # Default to GRB_NONBASIC_LOWER
+        
+        for var_idx in basis.var_basic
+            var_basis_statuses[var_idx] = 0  # GRB_BASIC
+        end
+        
+        ret = GRBsetintattrarray(gurobi_model, "VBasis", 0, nvars, var_basis_statuses)
+        if ret != 0
+            error("Failed to set variable basis statuses in Gurobi: error code $ret")
+        end
+    end
+
+    @timeit "Set constraint basis" begin
+        nconstr = num_constraints(model, AffExpr, MOI.EqualTo{Float64})
+        constr_basis_statuses = Vector{Cint}(undef, nconstr)
+        fill!(constr_basis_statuses, -1)  # Default to GRB_NONBASIC
+        
+        for constr_idx in basis.constr_basic
+            constr_basis_statuses[constr_idx] = 0  # GRB_BASIC
+        end
+        
+        ret = GRBsetintattrarray(gurobi_model, "CBasis", 0, nconstr, constr_basis_statuses)
+        if ret != 0
+            error("Failed to set constraint basis statuses in Gurobi: error code $ret")
+        end
+    end
+
+    return nothing
+end
+
+function get_basis_gurobi(model::JuMP.Model)::Basis
+    @timeit "Initialization" begin
+        var_basic = Int[]
+        var_nonbasic = Int[]
+        constr_basic = Int[]
+        constr_nonbasic = Int[]
+        nvars = num_variables(model)
+        sizehint!(var_basic, nvars)
+        sizehint!(var_nonbasic, nvars)
+        gurobi_model = unsafe_backend(model).inner
+    end
+
+    @timeit "Query variables" begin
+        var_basis_statuses = Vector{Cint}(undef, nvars)
+        ret = GRBgetintattrarray(gurobi_model, "VBasis", 0, nvars, var_basis_statuses)
+        if ret != 0
+            error("Failed to get variable basis statuses from Gurobi: error code $ret")
+        end
+        for i in 1:nvars
+            if var_basis_statuses[i] == 0  # GRB_BASIC
+                push!(var_basic, i)
+            elseif var_basis_statuses[i] == -1  # GRB_NONBASIC_LOWER
+                push!(var_nonbasic, i)
+            else
+                error("Unknown variable basis status: $(var_basis_statuses[i])")
+            end
+        end
+    end
+
+    @timeit "Query constraints" begin
+        nconstr = num_constraints(model, AffExpr, MOI.EqualTo{Float64})
+        constr_basis_statuses = Vector{Cint}(undef, nconstr)
+        ret = GRBgetintattrarray(gurobi_model, "CBasis", 0, nconstr, constr_basis_statuses)
+        if ret != 0
+            error("Failed to get constraint basis statuses from Gurobi: error code $ret")
+        end
+        for i in 1:nconstr
+            if constr_basis_statuses[i] == 0  # GRB_BASIC
+                push!(constr_basic, i)
+            elseif constr_basis_statuses[i] == -1  # GRB_NONBASIC
+                push!(constr_nonbasic, i)
+            else
+                error("Unknown constraint basis status: $(constr_basis_statuses[i])")
+            end
+        end
+    end
+
+    @timeit "Build basis struct" begin
+        basis = Basis(; var_basic, var_nonbasic, constr_basic, constr_nonbasic)
+    end
+
+    return basis
 end
 
 function get_x(model::JuMP.Model)
@@ -54,11 +167,16 @@ end
 
 function compute_tableau(
     data::ProblemData,
-    basis::Basis,
-    x::Vector{Float64};
+    basis::Basis;
+    x::Union{Nothing,Vector{Float64}} = nothing,
     rows::Union{Vector{Int},Nothing} = nothing,
     tol = 1e-8,
+    estimated_density = 0.10,
 )::Tableau
+    if isnan(estimated_density) || estimated_density <= 0
+        estimated_density = 0.10
+    end
+    
     @timeit "Split data" begin
         nrows, ncols = size(data.constr_lhs)
         lhs_slacks = sparse(I, nrows, nrows)
@@ -73,53 +191,80 @@ function compute_tableau(
         factor = klu(sparse(lhs_b'))
     end
 
-    @timeit "Compute tableau LHS" begin
-        tableau_lhs_I = Int[]
-        tableau_lhs_J = Int[]
-        tableau_lhs_V = Float64[]
-        for k = 1:length(rows)
-            @timeit "Prepare inputs" begin
-                i = rows[k]
-                e = zeros(nrows)
-                e[i] = 1.0
-            end
-            @timeit "Solve" begin
-                sol = factor \ e
-            end
-            @timeit "Multiply" begin
-                row = sol' * data.constr_lhs
-            end
-            @timeit "Sparsify & copy" begin
-                for (j, v) in enumerate(row)
-                    if abs(v) < tol
-                        continue
-                    end
-                    push!(tableau_lhs_I, k)
-                    push!(tableau_lhs_J, j)
-                    push!(tableau_lhs_V, v)
-                end
-            end
-        end
-        tableau_lhs =
-            sparse(tableau_lhs_I, tableau_lhs_J, tableau_lhs_V, length(rows), ncols)
+    @timeit "Initialize arrays" begin
+        num_rows = length(rows)
+        tableau_rhs::Array{Float64} = zeros(num_rows)
+        tableau_rowptr::Array{Int} = zeros(Int, num_rows + 1)
+        tableau_colval::Array{Int} = Int[]
+        tableau_nzval::Array{Float64} = Float64[]
+        estimated_nnz::Int = round(num_rows * ncols * estimated_density)
+        sizehint!(tableau_colval, estimated_nnz)
+        sizehint!(tableau_nzval, estimated_nnz)
+        e::Array{Float64} = zeros(nrows)
+        sol::Array{Float64} = zeros(nrows)
+        tableau_row::Array{Float64} = zeros(ncols)
     end
 
-    @timeit "Compute tableau RHS" begin
-        tableau_rhs = [x[basis.var_basic]; zeros(length(basis.constr_basic))][rows]
+    A = data.constr_lhs'
+    b = data.constr_ub
+    tableau_rowptr[1] = 1
+
+    @timeit "Process rows" begin
+        for k in eachindex(rows)
+            @timeit "Solve" begin
+                fill!(e, 0.0)
+                e[rows[k]] = 1.0
+                ldiv!(sol, factor, e)
+            end
+            @timeit "Compute row" begin
+                mul!(tableau_row, A, sol)
+                tableau_rhs[k] = dot(sol, b)
+            end
+            needed_space = length(tableau_colval) + ncols
+            if needed_space > estimated_nnz
+                @timeit "Grow arrays" begin
+                    estimated_nnz *= 2
+                    sizehint!(tableau_colval, estimated_nnz)
+                    sizehint!(tableau_nzval, estimated_nnz)
+                end
+            end
+            @timeit "Collect nonzeros for row" begin
+                for j in 1:ncols
+                    val = tableau_row[j]
+                    if abs(val) > tol
+                        push!(tableau_colval, j)
+                        push!(tableau_nzval, val)
+                    end
+                end
+            end
+            tableau_rowptr[k + 1] = length(tableau_colval) + 1
+        end
+    end
+
+    @timeit "Shrink arrays" begin
+        sizehint!(tableau_colval, length(tableau_colval))
+        sizehint!(tableau_nzval, length(tableau_nzval))
+    end
+
+    @timeit "Build sparse matrix" begin
+        tableau_lhs_transposed = SparseMatrixCSC(ncols, num_rows, tableau_rowptr, tableau_colval, tableau_nzval)
+        tableau_lhs = transpose(tableau_lhs_transposed)
     end
 
     @timeit "Compute tableau objective row" begin
         sol = factor \ obj_b
         tableau_obj = -data.obj' + sol' * data.constr_lhs
         tableau_obj[abs.(tableau_obj).<tol] .= 0
+        tableau_obj = Array(tableau_obj')
     end
 
-    return Tableau(
-        obj = tableau_obj,
-        lhs = tableau_lhs,
-        rhs = tableau_rhs,
-        z = dot(data.obj, x),
-    )
+    # Compute z if solution is provided
+    z = 0
+    if x !== nothing
+        z = dot(data.obj, x)
+    end
+
+    return Tableau(obj = tableau_obj, lhs = tableau_lhs, rhs = tableau_rhs, z = z)
 end
 
-export get_basis, get_x, compute_tableau
+export get_basis, get_basis_gurobi, set_basis, get_x, compute_tableau

src/Cuts/tableau/transform.jl

@@ -96,46 +96,70 @@ Base.@kwdef mutable struct AddSlackVariables <: Transform
 end
 
 function forward!(t::AddSlackVariables, data::ProblemData)
-    nrows, ncols = size(data.constr_lhs)
-    isequality = abs.(data.constr_ub .- data.constr_lb) .< 1e-6
-    eq = [i for i = 1:nrows if isequality[i]]
-    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)
-
-    t.M1 = [
-        I spzeros(ncols, GE + LE)
-        data.constr_lhs[ge, :] spzeros(GE, GE + LE)
-        -data.constr_lhs[le, :] spzeros(LE, GE + LE)
-    ]
-    t.M2 = [
-        zeros(ncols)
-        data.constr_lb[ge]
-        -data.constr_ub[le]
-    ]
-    t.ncols_orig = ncols
-    t.GE, t.LE = GE, LE
-    t.lhs_ge = data.constr_lhs[ge, :]
-    t.lhs_le = data.constr_lhs[le, :]
-    t.rhs_ge = data.constr_lb[ge]
-    t.rhs_le = data.constr_ub[le]
-
-    data.constr_lhs = [
-        data.constr_lhs[eq, :] spzeros(EQ, GE) spzeros(EQ, LE)
-        data.constr_lhs[ge, :] -I spzeros(GE, LE)
-        data.constr_lhs[le, :] spzeros(LE, GE) I
-    ]
-    data.obj = [data.obj; zeros(GE + LE)]
-    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.constr_lb = [
-        data.constr_lb[eq]
-        data.constr_lb[ge]
-        data.constr_ub[le]
-    ]
-    data.constr_ub = copy(data.constr_lb)
+    @timeit "Identify constraint type" begin
+        nrows, ncols = size(data.constr_lhs)
+        isequality = abs.(data.constr_ub .- data.constr_lb) .< 1e-6
+        eq = [i for i = 1:nrows if isequality[i]]
+        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)
+    end
+    @timeit "Identify slack type" begin
+        constr_lhs_t = sparse(data.constr_lhs')
+        function is_integral(row_idx, rhs)
+            rhs_is_integer = abs(rhs - round(rhs)) <= 1e-6
+            cols, coeffs = findnz(constr_lhs_t[:, row_idx])[1:2]
+            vars_are_integer = all(j -> data.var_types[j] ∈ ['I', 'B'], cols)
+            coeffs_are_integer = all(v -> abs(v - round(v)) <= 1e-6, coeffs)
+            return rhs_is_integer && vars_are_integer && coeffs_are_integer
+        end
+        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]
+        ]
+    end
+    @timeit "Build M1" begin
+        t.M1 = [
+            I spzeros(ncols, GE + LE)
+            data.constr_lhs[ge, :] spzeros(GE, GE + LE)
+            -data.constr_lhs[le, :] spzeros(LE, GE + LE)
+        ]
+    end
+    @timeit "Build M2" begin
+        t.M2 = [
+            zeros(ncols)
+            data.constr_lb[ge]
+            -data.constr_ub[le]
+        ]
+    end
+    @timeit "Build t.lhs, t.rhs" begin
+        t.ncols_orig = ncols
+        t.GE, t.LE = GE, LE
+        t.lhs_ge = data.constr_lhs[ge, :]
+        t.lhs_le = data.constr_lhs[le, :]
+        t.rhs_ge = data.constr_lb[ge]
+        t.rhs_le = data.constr_ub[le]
+    end
+    @timeit "Build data.constr_lhs" begin
+        data.constr_lhs = [
+            data.constr_lhs[eq, :] spzeros(EQ, GE) spzeros(EQ, LE)
+            data.constr_lhs[ge, :] -I spzeros(GE, LE)
+            data.constr_lhs[le, :] spzeros(LE, GE) I
+        ]
+    end
+    @timeit "Build other data fields" begin
+        data.obj = [data.obj; zeros(GE + LE)]
+        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; slack_types]
+        data.constr_lb = [
+            data.constr_lb[eq]
+            data.constr_lb[ge]
+            data.constr_ub[le]
+        ]
+        data.constr_ub = copy(data.constr_lb)
+    end
 end
 
 function backwards!(t::AddSlackVariables, c::ConstraintSet)
@@ -155,71 +179,55 @@ end
 # -----------------------------------------------------------------------------
 
 Base.@kwdef mutable struct SplitFreeVars <: Transform
-    F::Int = 0
-    B::Int = 0
-    free::Vector{Int} = []
-    others::Vector{Int} = []
+    ncols::Int = 0
+    is_var_free::Vector{Bool} = []
 end
 
 function forward!(t::SplitFreeVars, data::ProblemData)
     lhs = data.constr_lhs
     _, ncols = size(lhs)
-    free = [i for i = 1:ncols if !isfinite(data.var_lb[i]) && !isfinite(data.var_ub[i])]
-    others = [i for i = 1:ncols if isfinite(data.var_lb[i]) || isfinite(data.var_ub[i])]
-    t.F = length(free)
-    t.B = length(others)
-    t.free, t.others = free, others
+    is_var_free = [!isfinite(data.var_lb[i]) && !isfinite(data.var_ub[i]) for i = 1:ncols]
+    free_idx = findall(is_var_free)
     data.obj = [
-        data.obj[others]
-        data.obj[free]
-        -data.obj[free]
+        data.obj
+        [-data.obj[i] for i in free_idx]
     ]
-    data.constr_lhs = [lhs[:, others] lhs[:, free] -lhs[:, free]]
     data.var_lb = [
-        data.var_lb[others]
-        [0.0 for _ in free]
-        [0.0 for _ in free]
+        [is_var_free[i] ? 0.0 : data.var_lb[i] for i in 1:ncols]
+        [0 for _ in free_idx]
     ]
     data.var_ub = [
-        data.var_ub[others]
-        [Inf for _ in free]
-        [Inf for _ in free]
+        [is_var_free[i] ? Inf : data.var_ub[i] for i in 1:ncols]
+        [Inf for _ in free_idx]
     ]
     data.var_types = [
-        data.var_types[others]
-        data.var_types[free]
-        data.var_types[free]
+        data.var_types
+        [data.var_types[i] for i in free_idx]
     ]
     data.var_names = [
-        data.var_names[others]
-        ["$(v)_p" for v in data.var_names[free]]
-        ["$(v)_m" for v in data.var_names[free]]
+        data.var_names
+        ["$(data.var_names[i])_neg" for i in free_idx]
     ]
+    data.constr_lhs = [lhs -lhs[:, free_idx]]
+    t.is_var_free, t.ncols = is_var_free, ncols
 end
 
 function backwards!(t::SplitFreeVars, c::ConstraintSet)
-    # Convert GE constraints into LE
-    nrows, _ = size(c.lhs)
-    ge = [i for i = 1:nrows if isfinite(c.lb[i])]
-    c.ub[ge], c.lb[ge] = -c.lb[ge], -c.ub[ge]
-    c.lhs[ge, :] *= -1
+    ncols, is_var_free = t.ncols, t.is_var_free
+    free_idx = findall(is_var_free)
 
-    # Assert only LE constraints are left (EQ constraints are not supported)
-    @assert all(c.lb .== -Inf)
-
-    # Take minimum (weakest) coefficient
-    B, F = t.B, t.F
-    for i = 1:F
-        c.lhs[:, B+i] = min.(c.lhs[:, B+i], -c.lhs[:, B+F+i])
+    for (offset, var_idx) in enumerate(free_idx)
+        @assert c.lhs[:, var_idx] == -c.lhs[:, ncols+offset]
     end
-    c.lhs = c.lhs[:, 1:(B+F)]
+    c.lhs = c.lhs[:, 1:ncols]
 end
 
 function forward(t::SplitFreeVars, p::Vector{Float64})::Vector{Float64}
+    ncols, is_var_free = t.ncols, t.is_var_free
+    free_idx = findall(is_var_free)
     return [
-        p[t.others]
-        max.(p[t.free], 0)
-        max.(-p[t.free], 0)
+        [is_var_free[i] ? max(0, p[i]) : p[i] for i in 1:ncols]
+        [max(0, -p[i]) for i in free_idx]
     ]
 end
 

src/MIPLearn.jl

@@ -6,11 +6,14 @@ module MIPLearn
 
 using PyCall
 using SparseArrays
+using PrecompileTools: @setup_workload, @compile_workload
+
 
 include("collectors.jl")
 include("components.jl")
 include("extractors.jl")
 include("io.jl")
+include("problems/maxcut.jl")
 include("problems/setcover.jl")
 include("problems/stab.jl")
 include("problems/tsp.jl")
@@ -22,6 +25,7 @@ function __init__()
     __init_components__()
     __init_extractors__()
     __init_io__()
+    __init_problems_maxcut__()
     __init_problems_setcover__()
     __init_problems_stab__()
     __init_problems_tsp__()
@@ -32,4 +36,51 @@ end
 include("BB/BB.jl")
 include("Cuts/Cuts.jl")
 
+# Precompilation
+# =============================================================================
+
+# function __precompile_cuts__()
+#     function build_model(mps_filename)
+#         model = read_from_file(mps_filename)
+#         set_optimizer(model, SCIP.Optimizer)
+#         return JumpModel(model)
+#     end
+#     BASEDIR = dirname(@__FILE__)
+#     mps_filename = "$BASEDIR/../test/fixtures/bell5.mps.gz"
+#     h5_filename = "$BASEDIR/../test/fixtures/bell5.h5"
+#     collect_gmi_dual(
+#         mps_filename;
+#         optimizer=HiGHS.Optimizer,
+#         max_rounds = 10,
+#         max_cuts_per_round = 500,
+#     )
+#     knn = KnnDualGmiComponent(
+#         extractor = H5FieldsExtractor(instance_fields = ["static_var_obj_coeffs"]),
+#         k = 2,
+#     )
+#     knn.fit([h5_filename, h5_filename])
+#     solver = LearningSolver(
+#         components = [
+#             ExpertPrimalComponent(action = SetWarmStart()),
+#             knn,
+#         ],
+#         skip_lp = true,
+#     )
+#     solver.optimize(mps_filename, build_model)
+# end
+
+# @setup_workload begin
+#     using SCIP
+#     using HiGHS
+#     using MIPLearn.Cuts
+#     using PrecompileTools: @setup_workload, @compile_workload
+
+#     __init__()
+#     Cuts.__init__()
+
+#     @compile_workload begin
+#         __precompile_cuts__()
+#     end
+# end
+
 end # module

src/components.jl

@@ -53,8 +53,14 @@ function __init_components__()
     )
     copy!(SelectTopSolutions, pyimport("miplearn.components.primal.mem").SelectTopSolutions)
     copy!(MergeTopSolutions, pyimport("miplearn.components.primal.mem").MergeTopSolutions)
-    copy!(MemorizingCutsComponent, pyimport("miplearn.components.cuts.mem").MemorizingCutsComponent)
-    copy!(MemorizingLazyComponent, pyimport("miplearn.components.lazy.mem").MemorizingLazyComponent)
+    copy!(
+        MemorizingCutsComponent,
+        pyimport("miplearn.components.cuts.mem").MemorizingCutsComponent,
+    )
+    copy!(
+        MemorizingLazyComponent,
+        pyimport("miplearn.components.lazy.mem").MemorizingLazyComponent,
+    )
 end
 
 export MinProbabilityClassifier,

src/io.jl

@@ -39,14 +39,14 @@ end
 function PyObject(m::SparseMatrixCSC)
     pyimport("scipy.sparse").csc_matrix(
         (m.nzval, m.rowval .- 1, m.colptr .- 1),
-        shape=size(m),
+        shape = size(m),
     ).tocoo()
 end
 
 function write_jld2(
     objs::Vector,
     dirname::AbstractString;
-    prefix::AbstractString=""
+    prefix::AbstractString = "",
 )::Vector{String}
     mkpath(dirname)
     filenames = [@sprintf("%s/%s%05d.jld2", dirname, prefix, i) for i = 1:length(objs)]

src/problems/maxcut.jl

@@ -0,0 +1,31 @@
+#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
+#  Copyright (C) 2020-2025, UChicago Argonne, LLC. All rights reserved.
+#  Released under the modified BSD license. See COPYING.md for more details.
+
+using JuMP
+
+global MaxCutData = PyNULL()
+global MaxCutGenerator = PyNULL()
+
+function __init_problems_maxcut__()
+    copy!(MaxCutData, pyimport("miplearn.problems.maxcut").MaxCutData)
+    copy!(MaxCutGenerator, pyimport("miplearn.problems.maxcut").MaxCutGenerator)
+end
+
+function build_maxcut_model_jump(data::Any; optimizer)
+    if data isa String
+        data = read_pkl_gz(data)
+    end
+    nodes = collect(data.graph.nodes())
+    edges = collect(data.graph.edges())
+    model = Model(optimizer)
+    @variable(model, x[nodes], Bin)
+    @objective(
+        model,
+        Min,
+        sum(-data.weights[i] * x[e[1]] * (1 - x[e[2]]) for (i, e) in enumerate(edges))
+    )
+    return JumpModel(model)
+end
+
+export MaxCutData, MaxCutGenerator, build_maxcut_model_jump

src/problems/setcover.jl

@@ -13,7 +13,7 @@ function __init_problems_setcover__()
     copy!(SetCoverGenerator, pyimport("miplearn.problems.setcover").SetCoverGenerator)
 end
 
-function build_setcover_model_jump(data::Any; optimizer=HiGHS.Optimizer)
+function build_setcover_model_jump(data::Any; optimizer = HiGHS.Optimizer)
     if data isa String
         data = read_pkl_gz(data)
     end

src/problems/stab.jl

@@ -10,10 +10,13 @@ global MaxWeightStableSetGenerator = PyNULL()
 
 function __init_problems_stab__()
     copy!(MaxWeightStableSetData, pyimport("miplearn.problems.stab").MaxWeightStableSetData)
-    copy!(MaxWeightStableSetGenerator, pyimport("miplearn.problems.stab").MaxWeightStableSetGenerator)
+    copy!(
+        MaxWeightStableSetGenerator,
+        pyimport("miplearn.problems.stab").MaxWeightStableSetGenerator,
+    )
 end
 
-function build_stab_model_jump(data::Any; optimizer=HiGHS.Optimizer)
+function build_stab_model_jump(data::Any; optimizer = HiGHS.Optimizer)
     nx = pyimport("networkx")
 
     if data isa String
@@ -50,11 +53,7 @@ function build_stab_model_jump(data::Any; optimizer=HiGHS.Optimizer)
         end
     end
 
-    return JumpModel(
-        model,
-        cuts_separate=cuts_separate,
-        cuts_enforce=cuts_enforce,
-    )
+    return JumpModel(model, cuts_separate = cuts_separate, cuts_enforce = cuts_enforce)
 end
 
 export MaxWeightStableSetData, MaxWeightStableSetGenerator, build_stab_model_jump

src/problems/tsp.jl

@@ -9,7 +9,10 @@ global TravelingSalesmanGenerator = PyNULL()
 
 function __init_problems_tsp__()
     copy!(TravelingSalesmanData, pyimport("miplearn.problems.tsp").TravelingSalesmanData)
-    copy!(TravelingSalesmanGenerator, pyimport("miplearn.problems.tsp").TravelingSalesmanGenerator)
+    copy!(
+        TravelingSalesmanGenerator,
+        pyimport("miplearn.problems.tsp").TravelingSalesmanGenerator,
+    )
 end
 
 function build_tsp_model_jump(data::Any; optimizer)
@@ -19,17 +22,15 @@ function build_tsp_model_jump(data::Any; optimizer)
         data = read_pkl_gz(data)
     end
     model = Model(optimizer)
-    edges = [(i, j) for i in 1:data.n_cities for j in (i+1):data.n_cities]
+    edges = [(i, j) for i = 1:data.n_cities for j = (i+1):data.n_cities]
     x = @variable(model, x[edges], Bin)
-    @objective(model, Min, sum(
-        x[(i, j)] * data.distances[i, j] for (i, j) in edges
-    ))
+    @objective(model, Min, sum(x[(i, j)] * data.distances[i, j] for (i, j) in edges))
 
     # Eq: Must choose two edges adjacent to each node
     @constraint(
         model,
         eq_degree[i in 1:data.n_cities],
-        sum(x[(min(i, j), max(i, j))] for j in 1:data.n_cities if i != j) == 2
+        sum(x[(min(i, j), max(i, j))] for j = 1:data.n_cities if i != j) == 2
     )
 
     function lazy_separate(cb_data)
@@ -41,10 +42,8 @@ function build_tsp_model_jump(data::Any; optimizer)
         for component in nx.connected_components(graph)
             if length(component) < data.n_cities
                 cut_edges = [
-                    [e[1], e[2]]
-                    for e in edges
-                    if (e[1] ∈ component && e[2] ∉ component)
-                    ||
+                    [e[1], e[2]] for
+                    e in edges if (e[1] ∈ component && e[2] ∉ component) ||
                     (e[1] ∉ component && e[2] ∈ component)
                 ]
                 push!(violations, cut_edges)
@@ -63,9 +62,9 @@ function build_tsp_model_jump(data::Any; optimizer)
 
     return JumpModel(
         model,
-        lazy_enforce=lazy_enforce,
-        lazy_separate=lazy_separate,
-        lp_optimizer=optimizer,
+        lazy_enforce = lazy_enforce,
+        lazy_separate = lazy_separate,
+        lp_optimizer = optimizer,
     )
 end
 

src/solvers/jump.jl

@@ -18,7 +18,7 @@ Base.@kwdef mutable struct _JumpModelExtData
     cuts_separate::Union{Function,Nothing} = nothing
     lazy_enforce::Union{Function,Nothing} = nothing
     lazy_separate::Union{Function,Nothing} = nothing
-    lp_optimizer
+    lp_optimizer::Any
 end
 
 function JuMP.copy_extension_data(
@@ -26,9 +26,7 @@ function JuMP.copy_extension_data(
     new_model::AbstractModel,
     ::AbstractModel,
 )
-    new_model.ext[:miplearn] = _JumpModelExtData(
-        lp_optimizer=old_ext.lp_optimizer
-    )
+    new_model.ext[:miplearn] = _JumpModelExtData(lp_optimizer = old_ext.lp_optimizer)
 end
 
 # -----------------------------------------------------------------------------
@@ -91,14 +89,27 @@ function _extract_after_load_vars(model::JuMP.Model, h5)
         for v in vars
     ]
     types = [JuMP.is_binary(v) ? "B" : JuMP.is_integer(v) ? "I" : "C" for v in vars]
-    obj = objective_function(model, AffExpr)
-    obj_coeffs = [v ∈ keys(obj.terms) ? obj.terms[v] : 0.0 for v in vars]
+
+    # Linear obj terms
+    obj = objective_function(model, QuadExpr)
+    obj_coeffs_linear = [v ∈ keys(obj.aff.terms) ? obj.aff.terms[v] : 0.0 for v in vars]
+
+    # Quadratic obj terms
+    if length(obj.terms) > 0
+        nvars = length(vars)
+        obj_coeffs_quad = zeros(nvars, nvars)
+        for (pair, coeff) in obj.terms
+            obj_coeffs_quad[pair.a.index.value, pair.b.index.value] = coeff
+        end
+        h5.put_array("static_var_obj_coeffs_quad", obj_coeffs_quad)
+    end
+
     h5.put_array("static_var_names", to_str_array(JuMP.name.(vars)))
     h5.put_array("static_var_types", to_str_array(types))
     h5.put_array("static_var_lower_bounds", lb)
     h5.put_array("static_var_upper_bounds", ub)
-    h5.put_array("static_var_obj_coeffs", obj_coeffs)
-    h5.put_scalar("static_obj_offset", obj.constant)
+    h5.put_array("static_var_obj_coeffs", obj_coeffs_linear)
+    h5.put_scalar("static_obj_offset", obj.aff.constant)
 end
 
 function _extract_after_load_constrs(model::JuMP.Model, h5)
@@ -145,7 +156,7 @@ function _extract_after_load_constrs(model::JuMP.Model, h5)
         end
     end
     if isempty(names)
-        error("no model constraints found; note that MIPLearn ignores unnamed constraints")
+        return
     end
     lhs = sparse(lhs_rows, lhs_cols, lhs_values, length(rhs), JuMP.num_variables(model))
     h5.put_sparse("static_constr_lhs", lhs)
@@ -284,9 +295,11 @@ function _extract_after_mip(model::JuMP.Model, h5)
 
     # Slacks
     lhs = h5.get_sparse("static_constr_lhs")
-    rhs = h5.get_array("static_constr_rhs")
-    slacks = abs.(lhs * x - rhs)
-    h5.put_array("mip_constr_slacks", slacks)
+    if lhs !== nothing
+        rhs = h5.get_array("static_constr_rhs")
+        slacks = abs.(lhs * x - rhs)
+        h5.put_array("mip_constr_slacks", slacks)
+    end
 
     # Cuts and lazy constraints
     ext = model.ext[:miplearn]
@@ -297,7 +310,7 @@ end
 function _fix_variables(model::JuMP.Model, var_names, var_values, stats)
     vars = [variable_by_name(model, v) for v in var_names]
     for (i, var) in enumerate(vars)
-        fix(var, var_values[i], force=true)
+        fix(var, var_values[i], force = true)
     end
 end
 
@@ -392,19 +405,19 @@ function __init_solvers_jump__()
         function __init__(
             self,
             inner;
-            cuts_enforce::Union{Function,Nothing}=nothing,
-            cuts_separate::Union{Function,Nothing}=nothing,
-            lazy_enforce::Union{Function,Nothing}=nothing,
-            lazy_separate::Union{Function,Nothing}=nothing,
-            lp_optimizer=HiGHS.Optimizer,
+            cuts_enforce::Union{Function,Nothing} = nothing,
+            cuts_separate::Union{Function,Nothing} = nothing,
+            lazy_enforce::Union{Function,Nothing} = nothing,
+            lazy_separate::Union{Function,Nothing} = nothing,
+            lp_optimizer = HiGHS.Optimizer,
         )
             self.inner = inner
             self.inner.ext[:miplearn] = _JumpModelExtData(
-                cuts_enforce=cuts_enforce,
-                cuts_separate=cuts_separate,
-                lazy_enforce=lazy_enforce,
-                lazy_separate=lazy_separate,
-                lp_optimizer=lp_optimizer,
+                cuts_enforce = cuts_enforce,
+                cuts_separate = cuts_separate,
+                lazy_enforce = lazy_enforce,
+                lazy_separate = lazy_separate,
+                lp_optimizer = lp_optimizer,
             )
         end
 
@@ -414,7 +427,7 @@ function __init_solvers_jump__()
             constrs_lhs,
             constrs_sense,
             constrs_rhs,
-            stats=nothing,
+            stats = nothing,
         ) = _add_constrs(
             self.inner,
             from_str_array(var_names),
@@ -430,14 +443,14 @@ function __init_solvers_jump__()
 
         extract_after_mip(self, h5) = _extract_after_mip(self.inner, h5)
 
-        fix_variables(self, var_names, var_values, stats=nothing) =
+        fix_variables(self, var_names, var_values, stats = nothing) =
             _fix_variables(self.inner, from_str_array(var_names), var_values, stats)
 
         optimize(self) = _optimize(self.inner)
 
         relax(self) = Class(_relax(self.inner))
 
-        set_warm_starts(self, var_names, var_values, stats=nothing) =
+        set_warm_starts(self, var_names, var_values, stats = nothing) =
             _set_warm_starts(self.inner, from_str_array(var_names), var_values, stats)
 
         write(self, filename) = _write(self.inner, filename)

test/Project.toml

@@ -17,6 +17,7 @@ MIPLearn = "2b1277c3-b477-4c49-a15e-7ba350325c68"
 PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 SCIP = "82193955-e24f-5292-bf16-6f2c5261a85f"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]

test/src/BB/test_bb.jl

@@ -86,11 +86,7 @@ function bb_run(optimizer_name, optimizer; large = true)
                     BB.ReliabilityBranching(aggregation = :min, collect = true),
                 ]
                     h5 = H5File("$FIXTURES/$instance.h5")
-                    mip_lower_bound = h5.get_scalar("mip_lower_bound")
-                    mip_upper_bound = h5.get_scalar("mip_upper_bound")
-                    mip_sense = h5.get_scalar("mip_sense")
-                    mip_primal_bound =
-                        mip_sense == "min" ? mip_upper_bound : mip_lower_bound
+                    mip_obj_bound = h5.get_scalar("mip_obj_bound")
                     h5.file.close()
 
                     mip = BB.init(optimizer)
@@ -98,7 +94,7 @@ function bb_run(optimizer_name, optimizer; large = true)
                     @info optimizer_name, branch_rule, instance
                     @time BB.solve!(
                         mip,
-                        initial_primal_bound = mip_primal_bound,
+                        initial_primal_bound = mip_obj_bound,
                         print_interval = 1,
                         node_limit = 25,
                         branch_rule = branch_rule,

test/src/Cuts/tableau/test_gmi.jl

@@ -0,0 +1,23 @@
+#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
+#  Copyright (C) 2020-2024, UChicago Argonne, LLC. All rights reserved.
+#  Released under the modified BSD license. See COPYING.md for more details.
+
+using HiGHS
+
+function test_cuts_tableau_gmi()
+    mps_filename = "$BASEDIR/../fixtures/bell5.mps.gz"
+    h5_filename = "$BASEDIR/../fixtures/bell5.h5"
+    collect_gmi(mps_filename, optimizer = HiGHS.Optimizer)
+    h5 = H5File(h5_filename, "r")
+    try
+        cuts_lb = h5.get_array("cuts_lb")
+        cuts_ub = h5.get_array("cuts_ub")
+        cuts_lhs = h5.get_sparse("cuts_lhs")
+        n_cuts = length(cuts_lb)
+        @test n_cuts > 0
+        @test n_cuts == length(cuts_ub)
+        @test cuts_lhs.shape[1] == n_cuts
+    finally
+        h5.close()
+    end
+end

test/src/Cuts/tableau/test_gmi_dual.jl

@@ -0,0 +1,70 @@
+#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
+#  Copyright (C) 2020-2024, UChicago Argonne, LLC. All rights reserved.
+#  Released under the modified BSD license. See COPYING.md for more details.
+
+using SCIP
+using HiGHS
+using MIPLearn.Cuts
+
+function test_cuts_tableau_gmi_dual_collect()
+    mps_filename = "$BASEDIR/../fixtures/bell5.mps.gz"
+    h5_filename = "$BASEDIR/../fixtures/bell5.h5"
+    stats = collect_gmi_dual(mps_filename, optimizer = HiGHS.Optimizer)
+    h5 = H5File(h5_filename, "r")
+    try
+        cuts_basis_vars = h5.get_array("cuts_basis_vars")
+        cuts_basis_sizes = h5.get_array("cuts_basis_sizes")
+        cuts_rows = h5.get_array("cuts_rows")
+        @test size(cuts_basis_vars) == (15, 402)
+        @test size(cuts_basis_sizes) == (15, 4)
+        @test size(cuts_rows) == (15,)
+    finally
+        h5.close()
+    end
+end
+
+function test_cuts_tableau_gmi_dual_usage()
+    function build_model(mps_filename)
+        model = read_from_file(mps_filename)
+        set_optimizer(model, SCIP.Optimizer)
+        return JumpModel(model)
+    end
+
+    mps_filename = "$BASEDIR/../fixtures/bell5.mps.gz"
+    h5_filename = "$BASEDIR/../fixtures/bell5.h5"
+    rm(h5_filename, force=true)
+
+    # Run basic collector
+    bc = BasicCollector(write_mps = false, skip_lp = true)
+    bc.collect([mps_filename], build_model)
+
+    # Run dual GMI collector
+    @info "Running dual GMI collector..."
+    collect_gmi_dual(mps_filename, optimizer = HiGHS.Optimizer)
+
+    # # Test expert component
+    # solver = LearningSolver(
+    #     components = [
+    #         ExpertPrimalComponent(action = SetWarmStart()),
+    #         ExpertDualGmiComponent(),
+    #     ],
+    #     skip_lp = true,
+    # )
+    # solver.optimize(mps_filename, build_model)
+
+    # Test kNN component
+    knn = KnnDualGmiComponent(
+        extractor = H5FieldsExtractor(instance_fields = ["static_var_obj_coeffs"]),
+        k = 2,
+    )
+    knn.fit([h5_filename, h5_filename])
+    solver = LearningSolver(
+        components = [
+            ExpertPrimalComponent(action = SetWarmStart()),
+            knn,
+        ],
+        skip_lp = true,
+    )
+    solver.optimize(mps_filename, build_model)
+    return
+end

test/src/MIPLearnT.jl

@@ -19,9 +19,12 @@ include("BB/test_bb.jl")
 include("components/test_cuts.jl")
 include("components/test_lazy.jl")
 include("Cuts/BlackBox/test_cplex.jl")
+include("Cuts/tableau/test_gmi.jl")
+include("Cuts/tableau/test_gmi_dual.jl")
 include("problems/test_setcover.jl")
 include("problems/test_stab.jl")
 include("problems/test_tsp.jl")
+include("problems/test_maxcut.jl")
 include("solvers/test_jump.jl")
 include("test_io.jl")
 include("test_usage.jl")
@@ -35,6 +38,7 @@ function runtests()
         test_problems_setcover()
         test_problems_stab()
         test_problems_tsp()
+        test_problems_maxcut()
         test_solvers_jump()
         test_usage()
         test_cuts()
@@ -43,8 +47,8 @@ function runtests()
 end
 
 function format()
-    JuliaFormatter.format(BASEDIR, verbose=true)
-    JuliaFormatter.format("$BASEDIR/../../src", verbose=true)
+    JuliaFormatter.format(BASEDIR, verbose = true)
+    JuliaFormatter.format("$BASEDIR/../../src", verbose = true)
     return
 end
 

test/src/components/test_cuts.jl

@@ -10,34 +10,32 @@ function gen_stab()
     randint = pyimport("scipy.stats").randint
     np.random.seed(42)
     gen = MaxWeightStableSetGenerator(
-        w=uniform(10.0, scale=1.0),
-        n=randint(low=50, high=51),
-        p=uniform(loc=0.5, scale=0.0),
-        fix_graph=true,
+        w = uniform(10.0, scale = 1.0),
+        n = randint(low = 50, high = 51),
+        p = uniform(loc = 0.5, scale = 0.0),
+        fix_graph = true,
     )
     data = gen.generate(1)
-    data_filenames = write_pkl_gz(data, "$BASEDIR/../fixtures", prefix="stab-n50-")
+    data_filenames = write_pkl_gz(data, "$BASEDIR/../fixtures", prefix = "stab-n50-")
     collector = BasicCollector()
     collector.collect(
         data_filenames,
-        data -> build_stab_model_jump(data, optimizer=SCIP.Optimizer),
-        progress=true,
-        verbose=true,
+        data -> build_stab_model_jump(data, optimizer = SCIP.Optimizer),
+        progress = true,
+        verbose = true,
     )
 end
 
 function test_cuts()
     data_filenames = ["$BASEDIR/../fixtures/stab-n50-00000.pkl.gz"]
     clf = pyimport("sklearn.dummy").DummyClassifier()
-    extractor = H5FieldsExtractor(
-        instance_fields=["static_var_obj_coeffs"],
-    )
-    comp = MemorizingCutsComponent(clf=clf, extractor=extractor)
-    solver = LearningSolver(components=[comp])
+    extractor = H5FieldsExtractor(instance_fields = ["static_var_obj_coeffs"])
+    comp = MemorizingCutsComponent(clf = clf, extractor = extractor)
+    solver = LearningSolver(components = [comp])
     solver.fit(data_filenames)
-    stats = solver.optimize(
+    model, stats = solver.optimize(
         data_filenames[1],
-        data -> build_stab_model_jump(data, optimizer=SCIP.Optimizer),
+        data -> build_stab_model_jump(data, optimizer = SCIP.Optimizer),
     )
     @test stats["Cuts: AOT"] > 0
 end

test/src/components/test_lazy.jl

@@ -11,36 +11,34 @@ function gen_tsp()
     np.random.seed(42)
 
     gen = TravelingSalesmanGenerator(
-        x=uniform(loc=0.0, scale=1000.0),
-        y=uniform(loc=0.0, scale=1000.0),
-        n=randint(low=20, high=21),
-        gamma=uniform(loc=1.0, scale=0.25),
-        fix_cities=true,
-        round=true,
+        x = uniform(loc = 0.0, scale = 1000.0),
+        y = uniform(loc = 0.0, scale = 1000.0),
+        n = randint(low = 20, high = 21),
+        gamma = uniform(loc = 1.0, scale = 0.25),
+        fix_cities = true,
+        round = true,
     )
     data = gen.generate(1)
-    data_filenames = write_pkl_gz(data, "$BASEDIR/../fixtures", prefix="tsp-n20-")
+    data_filenames = write_pkl_gz(data, "$BASEDIR/../fixtures", prefix = "tsp-n20-")
     collector = BasicCollector()
     collector.collect(
         data_filenames,
-        data -> build_tsp_model_jump(data, optimizer=GLPK.Optimizer),
-        progress=true,
-        verbose=true,
+        data -> build_tsp_model_jump(data, optimizer = GLPK.Optimizer),
+        progress = true,
+        verbose = true,
     )
 end
 
 function test_lazy()
     data_filenames = ["$BASEDIR/../fixtures/tsp-n20-00000.pkl.gz"]
     clf = pyimport("sklearn.dummy").DummyClassifier()
-    extractor = H5FieldsExtractor(
-        instance_fields=["static_var_obj_coeffs"],
-    )
-    comp = MemorizingLazyComponent(clf=clf, extractor=extractor)
-    solver = LearningSolver(components=[comp])
+    extractor = H5FieldsExtractor(instance_fields = ["static_var_obj_coeffs"])
+    comp = MemorizingLazyComponent(clf = clf, extractor = extractor)
+    solver = LearningSolver(components = [comp])
     solver.fit(data_filenames)
-    stats = solver.optimize(
+    model, stats = solver.optimize(
         data_filenames[1],
-        data -> build_tsp_model_jump(data, optimizer=GLPK.Optimizer),
+        data -> build_tsp_model_jump(data, optimizer = GLPK.Optimizer),
     )
     @test stats["Lazy Constraints: AOT"] > 0
 end

test/src/problems/test_maxcut.jl

@@ -0,0 +1,54 @@
+#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
+#  Copyright (C) 2020-2025, UChicago Argonne, LLC. All rights reserved.
+#  Released under the modified BSD license. See COPYING.md for more details.
+
+using PyCall
+
+function test_problems_maxcut()
+    np = pyimport("numpy")
+    random = pyimport("random")
+    scipy_stats = pyimport("scipy.stats")
+    randint = scipy_stats.randint
+    uniform = scipy_stats.uniform
+
+    # Set random seed
+    random.seed(42)
+    np.random.seed(42)
+
+    # Build random instance
+    data = MaxCutGenerator(
+        n = randint(low = 10, high = 11),
+        p = uniform(loc = 0.5, scale = 0.0),
+        fix_graph = false,
+    ).generate(
+        1,
+    )[1]
+
+    # Build model
+    model = build_maxcut_model_jump(data, optimizer = SCIP.Optimizer)
+
+    # Check static features
+    h5 = H5File(tempname(), "w")
+    model.extract_after_load(h5)
+    obj_linear = h5.get_array("static_var_obj_coeffs")
+    obj_quad = h5.get_array("static_var_obj_coeffs_quad")
+    @test obj_linear == [3.0, 1.0, 3.0, 1.0, -1.0, 0.0, -1.0, 0.0, -1.0, 0.0]
+    @test obj_quad == [
+        0.0 0.0 -1.0 1.0 -1.0 0.0 0.0 0.0 -1.0 -1.0
+        0.0 0.0 1.0 -1.0 0.0 -1.0 -1.0 0.0 0.0 1.0
+        0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 -1.0 -1.0
+        0.0 0.0 0.0 0.0 0.0 -1.0 1.0 -1.0 0.0 0.0
+        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
+        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
+        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 -1.0
+        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
+        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+    ]
+
+    # Check optimal solution
+    model.optimize()
+    model.extract_after_mip(h5)
+    @test h5.get_scalar("mip_obj_value") == -4
+    h5.close()
+end

test/src/problems/test_stab.jl

@@ -8,11 +8,11 @@ using SCIP
 function test_problems_stab()
     nx = pyimport("networkx")
     data = MaxWeightStableSetData(
-        graph=nx.gnp_random_graph(25, 0.5, seed=42),
-        weights=repeat([1.0], 25),
+        graph = nx.gnp_random_graph(25, 0.5, seed = 42),
+        weights = repeat([1.0], 25),
     )
     h5 = H5File(tempname(), "w")
-    model = build_stab_model_jump(data, optimizer=SCIP.Optimizer)
+    model = build_stab_model_jump(data, optimizer = SCIP.Optimizer)
     model.extract_after_load(h5)
     model.optimize()
     model.extract_after_mip(h5)

test/src/problems/test_tsp.jl

@@ -10,17 +10,12 @@ function test_problems_tsp()
     squareform = pyimport("scipy.spatial.distance").squareform
 
     data = TravelingSalesmanData(
-        n_cities=6,
-        distances=squareform(pdist([
-            [0.0, 0.0],
-            [1.0, 0.0],
-            [2.0, 0.0],
-            [3.0, 0.0],
-            [0.0, 1.0],
-            [3.0, 1.0],
-        ])),
+        n_cities = 6,
+        distances = squareform(
+            pdist([[0.0, 0.0], [1.0, 0.0], [2.0, 0.0], [3.0, 0.0], [0.0, 1.0], [3.0, 1.0]]),
+        ),
     )
-    model = build_tsp_model_jump(data, optimizer=GLPK.Optimizer)
+    model = build_tsp_model_jump(data, optimizer = GLPK.Optimizer)
     model.optimize()
     @test objective_value(model.inner) == 8.0
     return

test/src/test_io.jl

@@ -4,6 +4,7 @@
 
 using MIPLearn
 using JLD2
+using SparseArrays
 
 struct _TestStruct
     n::Int
@@ -35,6 +36,8 @@ function test_h5()
     _test_roundtrip_array(h5, [1, 2, 3])
     _test_roundtrip_array(h5, [1.0, 2.0, 3.0])
     _test_roundtrip_str_array(h5, ["A", "BB", "CCC"])
+    _test_roundtrip_sparse(h5, sparse([1; 2; 3], [1; 2; 3], [1; 2; 3]))
+    # _test_roundtrip_sparse(h5, sparse([1; 2; 3], [1; 2; 3], [1; 2; 3], 4, 4))
     @test h5.get_array("unknown-key") === nothing
     h5.close()
 end
@@ -46,7 +49,7 @@ function test_jld2()
         _TestStruct(2, [1.0, 2.0, 3.0]),
         _TestStruct(3, [3.0, 3.0, 3.0]),
     ]
-    filenames = write_jld2(data, dirname, prefix="obj")
+    filenames = write_jld2(data, dirname, prefix = "obj")
     @test all(
         filenames .==
         ["$dirname/obj00001.jld2", "$dirname/obj00002.jld2", "$dirname/obj00003.jld2"],
@@ -79,3 +82,11 @@ function _test_roundtrip_str_array(h5, original)
     @test recovered !== nothing
     @test all(original .== recovered)
 end
+
+function _test_roundtrip_sparse(h5, original)
+    h5.put_sparse("key", original)
+    recovered = MIPLearn.convert(SparseMatrixCSC, h5.get_sparse("key"))
+    @test recovered !== nothing
+    @test size(original) == size(recovered)
+    @test all(original .== recovered)
+end

test/src/test_usage.jl

@@ -13,16 +13,16 @@ function test_usage()
 
     @debug "Setting up LearningSolver..."
     solver = LearningSolver(
-        components=[
+        components = [
             IndependentVarsPrimalComponent(
-                base_clf=SingleClassFix(
+                base_clf = SingleClassFix(
                     MinProbabilityClassifier(
-                        base_clf=LogisticRegression(),
-                        thresholds=[0.95, 0.95],
+                        base_clf = LogisticRegression(),
+                        thresholds = [0.95, 0.95],
                     ),
                 ),
-                extractor=AlvLouWeh2017Extractor(),
-                action=SetWarmStart(),
+                extractor = AlvLouWeh2017Extractor(),
+                action = SetWarmStart(),
             ),
         ],
     )