DualGMI: Skip empty H5 files

Project.toml

@@ -1,7 +1,7 @@
 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"
@@ -15,10 +15,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"
@@ -30,7 +32,6 @@ HDF5 = "0.16"
 HiGHS = "1"
 JLD2 = "0.4"
 JSON = "0.21"
-julia = "1"
 JuMP = "1"
 KLU = "0.4"
 MathOptInterface = "1"
@@ -39,3 +40,6 @@ PyCall = "1"
 Requires = "1"
 Statistics = "1"
 TimerOutputs = "0.5"
+julia = "1"
+PrecompileTools = "1"
+SCIP = "0.12"

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!()
 
-function select_gmi_rows(data, basis, x; max_rows = 10, atol = 0.001)
+    # 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 = 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 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,10 +199,10 @@ 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[]
@@ -30,23 +210,25 @@ function compute_gmi(data::ProblemData, tableau::Tableau, tol = 1e-8)::Constrain
     @timeit "Compute coefficients" begin
         for k = 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 +239,5 @@ 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
+export compute_gmi,
+    frac, select_gmi_rows, assert_cuts_off, assert_does_not_cut_off, collect_gmi

src/Cuts/tableau/gmi_dual.jl

@@ -0,0 +1,580 @@
+#  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.
+
+using Printf
+using JuMP
+using HiGHS
+using Random
+using DataStructures
+
+import ..H5FieldsExtractor
+
+global ExpertDualGmiComponent = PyNULL()
+global KnnDualGmiComponent = PyNULL()
+
+Base.@kwdef mutable struct _KnnDualGmiData
+    k = nothing
+    extractor = nothing
+    train_h5 = nothing
+    model = nothing
+    strategy = nothing
+end
+
+function collect_gmi_dual(
+    mps_filename;
+    optimizer,
+    max_rounds = 10,
+    max_cuts_per_round = 500,
+    time_limit = 3_600,
+)
+    reset_timer!()
+    initial_time = time()
+
+    @timeit "Read H5" begin
+        h5_filename = replace(mps_filename, ".mps.gz" => ".h5")
+        h5 = H5File(h5_filename, "r")
+        sol_opt_dict = Dict(
+            zip(
+                h5.get_array("static_var_names"),
+                convert(Array{Float64}, h5.get_array("mip_var_values")),
+            ),
+        )
+        obj_mip = h5.get_scalar("mip_obj_value")
+        h5.file.close()
+    end
+
+    # Define relative MIP gap
+    gap(v) = 100 * abs(obj_mip - v) / abs(obj_mip)
+
+    @timeit "Initialize" begin
+        stats_obj = []
+        stats_gap = []
+        stats_ncuts = []
+        original_basis = nothing
+        all_cuts = nothing
+        all_cuts_bases = nothing
+        all_cuts_rows = nothing
+        last_round_obj = nothing
+    end
+
+    @timeit "Read problem" begin
+        model = read_from_file(mps_filename)
+        set_optimizer(model, optimizer)
+        obj_original = objective_function(model)
+    end
+
+    for round = 1:max_rounds
+        @info "Round $(round)..."
+
+        @timeit "Convert model to standard form" begin
+            # Extract problem data
+            data = ProblemData(model)
+
+            # Construct optimal solution vector (with correct variable sequence)
+            sol_opt = [sol_opt_dict[n] for n in data.var_names]
+
+            # Assert optimal solution is feasible for the original problem
+            assert_leq(data.constr_lb, data.constr_lhs * sol_opt)
+            assert_leq(data.constr_lhs * sol_opt, data.constr_ub)
+
+            # Convert to standard form
+            data_s, transforms = convert_to_standard_form(data)
+            model_s = to_model(data_s)
+            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
+
+        @timeit "Optimize standard model" begin
+            @info "Optimizing standard model..."
+            optimize!(model_s)
+            obj = objective_value(model_s)
+            if round == 1
+                push!(stats_obj, obj)
+                push!(stats_gap, gap(obj))
+                push!(stats_ncuts, 0)
+            else
+                if obj ≈ last_round_obj
+                    @info ("No improvement in obj value. Aborting.")
+                    break
+                end
+            end
+            if termination_status(model_s) != MOI.OPTIMAL
+                error("Non-optimal termination status")
+            end
+            last_round_obj = obj
+        end
+
+        @timeit "Select tableau rows" begin
+            basis = get_basis(model_s)
+            if round == 1
+                original_basis = basis
+            end
+            sol_frac = get_x(model_s)
+            selected_rows =
+                select_gmi_rows(data_s, basis, sol_frac, max_rows = max_cuts_per_round)
+        end
+
+        @timeit "Compute tableau rows" begin
+            tableau = compute_tableau(data_s, basis, x = sol_frac, rows = selected_rows)
+
+            # Assert tableau rows have been computed correctly
+            assert_eq(tableau.lhs * sol_frac, tableau.rhs, atol=1e-3)
+            assert_eq(tableau.lhs * sol_opt_s, tableau.rhs, atol=1e-3)
+        end
+
+        @timeit "Compute GMI cuts" begin
+            cuts_s = compute_gmi(data_s, tableau)
+
+            # Assert cuts 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. Aborting."
+                break
+            else
+                @info "Generated $(length(cuts_s.lb)) cuts"
+            end
+        end
+
+        @timeit "Add GMI cuts to original model" begin
+            @timeit "Convert to original form" begin
+                cuts = backwards(transforms, cuts_s)
+            end
+
+            @timeit "Prepare bv" begin
+                bv = repeat([basis], length(selected_rows))
+            end
+
+            @timeit "Append matrices" begin
+                if round == 1
+                    all_cuts = cuts
+                    all_cuts_bases = bv
+                    all_cuts_rows = selected_rows
+                else
+                    all_cuts.lhs = [all_cuts.lhs; cuts.lhs]
+                    all_cuts.lb = [all_cuts.lb; cuts.lb]
+                    all_cuts.ub = [all_cuts.ub; cuts.ub]
+                    all_cuts_bases = [all_cuts_bases; bv]
+                    all_cuts_rows = [all_cuts_rows; selected_rows]
+                end
+            end
+
+            @timeit "Add to model" begin
+                @info "Adding $(length(all_cuts.lb)) constraints to original model"
+                constrs, gmi_exps = add_constraint_set_dual_v2(model, all_cuts)
+            end
+        end
+
+        @timeit "Optimize original model" begin
+            set_objective_function(model, obj_original)
+            undo_relax = relax_integrality(model)
+            @info "Optimizing original model (constr)..."
+            optimize!(model)
+            obj = objective_value(model)
+            push!(stats_obj, obj)
+            push!(stats_gap, gap(obj))
+            sp = [shadow_price(c) for c in constrs]
+            undo_relax()
+            useful = [abs(sp[i]) > 1e-6 for (i, _) in enumerate(constrs)]
+            keep = findall(useful .== true)
+        end
+
+        @timeit "Filter out useless cuts" begin
+            @info "Keeping $(length(keep)) useful cuts"
+            all_cuts.lhs = all_cuts.lhs[keep, :]
+            all_cuts.lb = all_cuts.lb[keep]
+            all_cuts.ub = all_cuts.ub[keep]
+            all_cuts_bases = all_cuts_bases[keep, :]
+            all_cuts_rows = all_cuts_rows[keep, :]
+            push!(stats_ncuts, length(all_cuts_rows))
+            if isempty(keep)
+                break
+            end
+        end
+
+        @timeit "Update obj function of original model" begin
+            delete.(model, constrs)
+            set_objective_function(
+                model,
+                obj_original -
+                sum(sp[i] * gmi_exps[i] for (i, c) in enumerate(constrs) if useful[i]),
+            )
+        end
+
+        elapsed_time = time() - initial_time
+        if elapsed_time > time_limit
+            @info "Time limit exceeded. Stopping."
+            break
+        end
+    end
+
+    @timeit "Store cuts in H5 file" begin
+        if all_cuts !== nothing
+            ncuts = length(all_cuts_rows)
+            total =
+                length(original_basis.var_basic) +
+                length(original_basis.var_nonbasic) +
+                length(original_basis.constr_basic) +
+                length(original_basis.constr_nonbasic)
+            all_cuts_basis_sizes = Array{Int64,2}(undef, ncuts, 4)
+            all_cuts_basis_vars = Array{Int64,2}(undef, ncuts, total)
+            for i = 1:ncuts
+                vb = all_cuts_bases[i].var_basic
+                vn = all_cuts_bases[i].var_nonbasic
+                cb = all_cuts_bases[i].constr_basic
+                cn = all_cuts_bases[i].constr_nonbasic
+                all_cuts_basis_sizes[i, :] = [length(vb) length(vn) length(cb) length(cn)]
+                all_cuts_basis_vars[i, :] = [vb' vn' cb' cn']
+            end
+            @info "Storing $(length(all_cuts.ub)) GMI cuts..."
+            h5 = H5File(h5_filename)
+            h5.put_sparse("cuts_lhs", all_cuts.lhs)
+            h5.put_array("cuts_lb", all_cuts.lb)
+            h5.put_array("cuts_ub", all_cuts.ub)
+            h5.put_array("cuts_basis_vars", all_cuts_basis_vars)
+            h5.put_array("cuts_basis_sizes", all_cuts_basis_sizes)
+            h5.put_array("cuts_rows", all_cuts_rows)
+            h5.file.close()
+        end
+    end
+
+    to = TimerOutputs.get_defaulttimer()
+    stats_time = TimerOutputs.tottime(to) / 1e9
+    print_timer()
+
+    return OrderedDict(
+        "instance" => mps_filename,
+        "max_rounds" => max_rounds,
+        "rounds" => length(stats_obj) - 1,
+        "obj_mip" => obj_mip,
+        "stats_obj" => stats_obj,
+        "stats_gap" => stats_gap,
+        "stats_ncuts" => stats_ncuts,
+        "stats_time" => stats_time,
+    )
+end
+
+function add_constraint_set_dual_v2(model::JuMP.Model, cs::ConstraintSet)
+    vars = all_variables(model)
+    nrows, ncols = size(cs.lhs)
+
+    @timeit "Transpose LHS" begin
+        lhs_t = spzeros(ncols, nrows)
+        ftranspose!(lhs_t, cs.lhs, x -> x)
+        lhs_t_rows = rowvals(lhs_t)
+        lhs_t_vals = nonzeros(lhs_t)
+    end
+
+    constrs = []
+    gmi_exps = []
+    for i = 1:nrows
+        c = nothing
+        gmi_exp = nothing
+        gmi_exp2 = nothing
+        @timeit "Build expr" begin
+            expr = AffExpr()
+            for k in nzrange(lhs_t, i)
+                add_to_expression!(expr, lhs_t_vals[k], vars[lhs_t_rows[k]])
+            end
+        end
+        @timeit "Add constraints" begin
+            if isinf(cs.ub[i])
+                c = @constraint(model, cs.lb[i] <= expr)
+                gmi_exp = cs.lb[i] - expr
+            elseif isinf(cs.lb[i])
+                c = @constraint(model, expr <= cs.ub[i])
+                gmi_exp = expr - cs.ub[i]
+            else
+                c = @constraint(model, cs.lb[i] <= expr <= cs.ub[i])
+                gmi_exp = cs.lb[i] - expr
+                gmi_exp2 = expr - cs.ub[i]
+            end
+        end
+        @timeit "Update structs" begin
+            push!(constrs, c)
+            push!(gmi_exps, gmi_exp)
+            if !isnothing(gmi_exp2)
+                push!(gmi_exps, gmi_exp2)
+            end
+        end
+    end
+    return constrs, gmi_exps
+end
+
+function _dualgmi_features(h5_filename, extractor)
+    h5 = H5File(h5_filename, "r")
+    try
+        return extractor.get_instance_features(h5)
+    finally
+        h5.close()
+    end
+end
+
+function _dualgmi_compress_h5(h5_filename)
+    vars_to_basis_offset = Dict()
+    basis_vars = []
+    basis_sizes = []
+    cut_basis::Array{Int} = []
+    cut_row::Array{Int} = []
+
+    h5 = H5File(h5_filename, "r")
+    orig_cut_basis_vars = h5.get_array("cuts_basis_vars")
+    orig_cut_basis_sizes = h5.get_array("cuts_basis_sizes")
+    orig_cut_rows = h5.get_array("cuts_rows")
+    h5.close()
+    if orig_cut_basis_vars === nothing
+        @warn "orig_cut_basis_vars is null; skipping file"
+        return
+    end
+    ncuts, _ = size(orig_cut_basis_vars)
+    if ncuts == 0
+        return
+    end
+
+    for i in 1:ncuts
+        vars = orig_cut_basis_vars[i, :]
+        sizes = orig_cut_basis_sizes[i, :]
+        row = orig_cut_rows[i]
+        if vars ∉ keys(vars_to_basis_offset)
+            offset = size(basis_vars)[1] + 1
+            vars_to_basis_offset[vars] = offset
+            push!(basis_vars, vars)
+            push!(basis_sizes, sizes)
+        end
+        offset = vars_to_basis_offset[vars]
+        push!(cut_basis, offset)
+        push!(cut_row, row)
+    end
+
+    basis_vars = hcat(basis_vars...)'
+    basis_sizes = hcat(basis_sizes...)'
+    _, n_vars = size(basis_vars)
+    if n_vars == 0
+        @warn "n_vars is zero; skipping file"
+        return
+    end
+
+    h5 = H5File(h5_filename, "r+")
+    h5.put_array("gmi_basis_vars", basis_vars)
+    h5.put_array("gmi_basis_sizes", basis_sizes)
+    h5.put_array("gmi_cut_basis", cut_basis)
+    h5.put_array("gmi_cut_row", cut_row)
+    h5.file.close()
+end
+
+function _dualgmi_generate(train_h5, model)
+    @timeit "Read problem data" begin
+        data = ProblemData(model)
+    end
+    @timeit "Convert to standard form" begin
+        data_s, transforms = convert_to_standard_form(data)
+    end
+    @timeit "Collect cuts from H5 files" begin
+        basis_vars_to_basis_offset = Dict()
+        combined_basis_sizes = nothing
+        combined_basis_sizes_list = Any[]
+        combined_basis_vars = nothing
+        combined_basis_vars_list = Any[]
+        combined_cut_rows = Any[]
+        for h5_filename in train_h5
+            @timeit "get_array (new)" begin
+                h5 = H5File(h5_filename, "r")
+                gmi_basis_vars = h5.get_array("gmi_basis_vars")
+                if gmi_basis_vars === nothing
+                    @warn "$(h5_filename) does not contain gmi_basis_vars; skipping"
+                    continue
+                end
+                gmi_basis_sizes = h5.get_array("gmi_basis_sizes")
+                gmi_cut_basis = h5.get_array("gmi_cut_basis")
+                gmi_cut_row = h5.get_array("gmi_cut_row")
+                h5.close()
+            end
+            @timeit "combine basis" begin
+                nbasis, _ = size(gmi_basis_vars)
+                local_to_combined_offset = Dict()
+                for local_offset in 1:nbasis
+                    vars = gmi_basis_vars[local_offset, :]
+                    sizes = gmi_basis_sizes[local_offset, :]
+                    if vars ∉ keys(basis_vars_to_basis_offset)
+                        combined_offset = length(combined_basis_vars_list) + 1
+                        basis_vars_to_basis_offset[vars] = combined_offset
+                        push!(combined_basis_vars_list, vars)
+                        push!(combined_basis_sizes_list, sizes)
+                        push!(combined_cut_rows, Set{Int}())
+                    end
+                    combined_offset = basis_vars_to_basis_offset[vars]
+                    local_to_combined_offset[local_offset] = combined_offset
+                end
+            end
+            @timeit "combine rows" begin
+                ncuts = length(gmi_cut_row)
+                for i in 1:ncuts
+                    local_offset = gmi_cut_basis[i]
+                    combined_offset = local_to_combined_offset[local_offset]
+                    row = gmi_cut_row[i]
+                    push!(combined_cut_rows[combined_offset], row)
+                end
+            end
+            @timeit "convert lists to matrices" begin
+                combined_basis_vars = hcat(combined_basis_vars_list...)'
+                combined_basis_sizes = hcat(combined_basis_sizes_list...)'
+            end
+        end
+    end
+    @timeit "Compute tableaus and cuts" begin
+        all_cuts = nothing
+        nbasis = length(combined_cut_rows)
+        for offset in 1:nbasis
+            rows = combined_cut_rows[offset]
+            try
+                vbb, vnn, cbb, cnn = combined_basis_sizes[offset, :]
+                current_basis = Basis(;
+                    var_basic = combined_basis_vars[offset, 1:vbb],
+                    var_nonbasic = combined_basis_vars[offset, vbb+1:vbb+vnn],
+                    constr_basic = combined_basis_vars[offset, vbb+vnn+1:vbb+vnn+cbb],
+                    constr_nonbasic = combined_basis_vars[offset, vbb+vnn+cbb+1:vbb+vnn+cbb+cnn],
+                )
+                tableau = compute_tableau(data_s, current_basis; rows=collect(rows))
+                cuts_s = compute_gmi(data_s, tableau)
+                cuts = backwards(transforms, cuts_s)
+                if all_cuts === nothing
+                    all_cuts = cuts
+                else
+                    all_cuts.lhs = [all_cuts.lhs; cuts.lhs]
+                    all_cuts.lb = [all_cuts.lb; cuts.lb]
+                    all_cuts.ub = [all_cuts.ub; cuts.ub]
+                end
+            catch e
+                if isa(e, AssertionError)
+                    @warn "Numerical error detected. Skipping cuts from current tableau."
+                    continue
+                else
+                    rethrow(e)
+                end
+            end
+        end
+    end
+    return all_cuts
+end
+
+function _dualgmi_set_callback(model, all_cuts)
+    function cut_callback(cb_data)
+        if all_cuts !== nothing
+            constrs = build_constraints(model, all_cuts)
+            @info "Enforcing $(length(constrs)) cuts..."
+            for c in constrs
+                MOI.submit(model, MOI.UserCut(cb_data), c)
+            end
+            all_cuts = nothing
+        end
+    end
+    set_attribute(model, MOI.UserCutCallback(), cut_callback)
+end
+
+function KnnDualGmiComponent_fit(data::_KnnDualGmiData, train_h5)
+    x = hcat([_dualgmi_features(filename, data.extractor) for filename in train_h5]...)'
+    model = pyimport("sklearn.neighbors").NearestNeighbors(n_neighbors = length(train_h5))
+    model.fit(x)
+    data.model = model
+    data.train_h5 = train_h5
+end
+
+
+function KnnDualGmiComponent_before_mip(data::_KnnDualGmiData, test_h5, model, _)
+    reset_timer!()
+
+    @timeit "Extract features" begin
+        x = _dualgmi_features(test_h5, data.extractor)
+        x = reshape(x, 1, length(x))
+    end
+
+    @timeit "Find neighbors" begin
+        neigh_dist, neigh_ind = data.model.kneighbors(x, return_distance = true)
+        neigh_ind = neigh_ind .+ 1
+        N = length(neigh_dist)
+        k = min(N, data.k)
+
+        if data.strategy == "near"
+            selected = collect(1:k)
+        elseif data.strategy == "far"
+            selected = collect((N - k + 1) : N)
+        elseif data.strategy == "rand"
+            selected = shuffle(collect(1:N))[1:k]
+        else
+            error("unknown strategy: $(data.strategy)")
+        end
+
+        @info "Dual GMI: Selected neighbors ($(data.strategy)):"
+        neigh_dist = neigh_dist[selected]
+        neigh_ind = neigh_ind[selected]
+        for i in 1:k
+            h5_filename = data.train_h5[neigh_ind[i]]
+            dist = neigh_dist[i]
+            @info "    $(h5_filename) dist=$(dist)"
+        end
+    end
+
+    @info "Dual GMI: Generating cuts..."
+    @timeit "Generate cuts" begin
+        time_generate = @elapsed begin
+            cuts = _dualgmi_generate(data.train_h5[neigh_ind], model)
+        end
+        @info "Dual GMI: Generated $(length(cuts.lb)) unique cuts in $(time_generate) seconds"
+    end
+
+    @timeit "Set callback" begin
+        _dualgmi_set_callback(model, cuts)
+    end
+
+    print_timer()
+
+    stats = Dict()
+    stats["KnnDualGmi: k"] = k
+    stats["KnnDualGmi: strategy"] = data.strategy
+    stats["KnnDualGmi: cuts"] = length(cuts.lb)
+    stats["KnnDualGmi: time generate"] = time_generate
+    return stats
+end
+
+function __init_gmi_dual__()
+    @pydef mutable struct KnnDualGmiComponentPy
+        function __init__(self; extractor, k = 3, strategy = "near")
+            self.data = _KnnDualGmiData(; extractor, k, strategy)
+        end
+        function fit(self, train_h5)
+            KnnDualGmiComponent_fit(self.data, train_h5)
+        end
+        function before_mip(self, test_h5, model, stats)
+            return @time KnnDualGmiComponent_before_mip(self.data, test_h5, model.inner, stats)
+        end
+    end
+    copy!(KnnDualGmiComponent, KnnDualGmiComponentPy)
+
+    @pydef mutable struct ExpertDualGmiComponentPy
+        function __init__(self)
+            self.inner = KnnDualGmiComponentPy(
+                extractor=H5FieldsExtractor(instance_fields=["static_var_obj_coeffs"]),
+                k=1,
+            )
+        end
+        function fit(self, train_h5)
+        end
+        function before_mip(self, test_h5, model, stats)
+            self.inner.fit([test_h5])
+            return self.inner.before_mip(test_h5, model, stats)
+        end
+    end
+    copy!(ExpertDualGmiComponent, ExpertDualGmiComponentPy)
+end
+
+export collect_gmi_dual, expert_gmi_dual, ExpertDualGmiComponent, KnnDualGmiComponent
+

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

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

src/Cuts/tableau/tableau.jl

@@ -54,8 +54,8 @@ 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,
 )::Tableau
@@ -73,11 +73,12 @@ 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 "Compute tableau" begin
+        @timeit "Initialize" begin
+            tableau_rhs = zeros(length(rows))
+            tableau_lhs = zeros(length(rows), ncols)
+        end
+        for k in eachindex(1:length(rows))
             @timeit "Prepare inputs" begin
                 i = rows[k]
                 e = zeros(nrows)
@@ -87,25 +88,14 @@ function compute_tableau(
                 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)
+                tableau_lhs[k, :] = sol' * data.constr_lhs
+                tableau_rhs[k] = sol' * data.constr_ub
             end
         end
+        @timeit "Sparsify" begin
+            tableau_lhs[abs.(tableau_lhs) .<= tol] .= 0
+            tableau_lhs = sparse(tableau_lhs)
         end
-        tableau_lhs =
-            sparse(tableau_lhs_I, tableau_lhs_J, tableau_lhs_V, length(rows), ncols)
-    end
-
-    @timeit "Compute tableau RHS" begin
-        tableau_rhs = [x[basis.var_basic]; zeros(length(basis.constr_basic))][rows]
     end
 
     @timeit "Compute tableau objective row" begin
@@ -114,12 +104,13 @@ function compute_tableau(
         tableau_obj[abs.(tableau_obj).<tol] .= 0
     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

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

@@ -46,7 +46,7 @@ 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/stab.jl

@@ -10,7 +10,10 @@ 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)
@@ -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)

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")
+    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]

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()

test/src/components/test_cuts.jl

@@ -29,13 +29,11 @@ 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"],
-    )
+    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),
     )

test/src/components/test_lazy.jl

@@ -32,13 +32,11 @@ 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"],
-    )
+    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),
     )

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_tsp.jl

@@ -11,14 +11,9 @@ function test_problems_tsp()
 
     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],
-        ])),
+        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.optimize()

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