Merge branch 'dev' into feature/replay

Project.toml

@@ -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,4 @@ PyCall = "1"
 Requires = "1"
 Statistics = "1"
 TimerOutputs = "0.5"
+julia = "1"

src/BB/log.jl

@@ -31,9 +31,9 @@ function print_progress(
     node::Node;
     time_elapsed::Float64,
     print_interval::Int,
-    primal_update::Bool,
+    primal_update::Bool
 )::Nothing
-    if (pool.processed % print_interval == 0) || isempty(pool.pending) || primal_update
+    if (pool.processed % print_interval == 0) || isempty(pool.pending)
         if isempty(node.branch_vars)
             branch_var_name = "---"
             branch_lb = "---"

src/BB/nodepool.jl

@@ -11,7 +11,7 @@ function take(
     suggestions::Array{Node}=[],
     time_remaining::Float64,
     gap_limit::Float64,
-    node_limit::Int,
+    node_limit::Int
 )::Union{Symbol,Node}
     t = threadid()
     lock(pool.lock) do
@@ -54,7 +54,7 @@ function offer(
     parent_node::Union{Nothing,Node},
     child_nodes::Vector{Node},
     time_elapsed::Float64=0.0,
-    print_interval::Int = 100,
+    print_interval::Int=100
 )::Nothing
     lock(pool.lock) do
         primal_update = false
@@ -101,6 +101,7 @@ function offer(
             # Update branching variable history
             branch_var = child_nodes[1].branch_vars[end]
             offset = findfirst(isequal(branch_var), parent_node.fractional_variables)
+            if offset !== nothing
                 x = parent_node.fractional_values[offset]
                 obj_change_up = child_nodes[1].obj - parent_node.obj
                 obj_change_down = child_nodes[2].obj - parent_node.obj
@@ -117,6 +118,7 @@ function offer(
                 pool.history.avg_pseudocost_down =
                     mean(vh.pseudocost_down for vh in values(pool.var_history))
             end
+        end
 
         for node in child_nodes
             print_progress(
@@ -136,7 +138,7 @@ function _update_var_history(;
     var::Variable,
     x::Float64,
     obj_change_down::Float64,
-    obj_change_up::Float64,
+    obj_change_up::Float64
 )::Nothing
     # Create new history entry
     if var ∉ keys(pool.var_history)

src/BB/optimize.jl

@@ -17,9 +17,15 @@ function solve!(
     initial_primal_bound::Float64=Inf,
     branch_rule::VariableBranchingRule=ReliabilityBranching(),
     enable_plunging=true,
-)::NodePool
+    replay=nothing
+)::Tuple{NodePool,ReplayInfo}
+
+    if replay === nothing
+        replay = ReplayInfo()
+    end
+
     time_initial = time()
-    pool = NodePool(mip = mip)
+    pool = NodePool(mip=mip, next_index=replay.next_index)
     pool.primal_bound = initial_primal_bound
 
     root_node = _create_node(mip)
@@ -64,10 +70,25 @@ function solve!(
                 @assert status == :Optimal
                 _unset_node_bounds(node)
 
+                if node.index in keys(replay.node_decisions)
+                    decision = replay.node_decisions[node.index]
+                    ids = decision.ids
+                    branch_var = decision.branch_var
+                    var_value = decision.var_value
+                else
                     # Find branching variable
                     ids = generate_indices(pool, 2)
                     branch_var = find_branching_var(branch_rule, node, pool)
 
+                    # Query current fractional value
+                    offset = findfirst(isequal(branch_var), node.fractional_variables)
+                    var_value = node.fractional_values[offset]
+
+                    # Update replay
+                    decision = ReplayNodeDecision(; branch_var, var_value, ids)
+                    replay.node_decisions[node.index] = decision
+                end
+
                 # Find current variable lower and upper bounds
                 offset = findfirst(isequal(branch_var), mip.int_vars)
                 var_lb = mip.int_vars_lb[offset]
@@ -79,10 +100,6 @@ function solve!(
                     end
                 end
 
-                # Query current fractional value
-                offset = findfirst(isequal(branch_var), node.fractional_variables)
-                var_value = node.fractional_values[offset]
-
                 child_zero = _create_node(
                     mip,
                     index=ids[2],
@@ -109,7 +126,8 @@ function solve!(
             end
         end
     end
-    return pool
+    replay.next_index = pool.next_index
+    return pool, replay
 end
 
 function _create_node(
@@ -118,7 +136,7 @@ function _create_node(
     parent::Union{Nothing,Node}=nothing,
     branch_var::Union{Nothing,Variable}=nothing,
     branch_var_lb::Union{Nothing,Float64}=nothing,
-    branch_var_ub::Union{Nothing,Float64} = nothing,
+    branch_var_ub::Union{Nothing,Float64}=nothing
 )::Node
     if parent === nothing
         branch_vars = Variable[]

src/BB/structs.jl

@@ -72,3 +72,14 @@ Base.@kwdef mutable struct NodePool
     history::History = History()
     var_history::Dict{Variable,VariableHistory} = Dict()
 end
+
+Base.@kwdef struct ReplayNodeDecision
+    branch_var
+    var_value
+    ids
+end
+
+Base.@kwdef mutable struct ReplayInfo
+    node_decisions::Dict{Int,ReplayNodeDecision} = Dict()
+    next_index::Int = 1
+end

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,625 @@
+#  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
+
+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,
+)
+    reset_timer!()
+
+    @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
+    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 ExpertDualGmiComponent_before_mip(test_h5, model, _)
+    # Read cuts and optimal solution
+    h5 = H5File(test_h5, "r")
+    sol_opt_dict = Dict(
+        zip(
+            h5.get_array("static_var_names"),
+            convert(Array{Float64}, h5.get_array("mip_var_values")),
+        ),
+    )
+    cut_basis_vars = h5.get_array("cuts_basis_vars")
+    cut_basis_sizes = h5.get_array("cuts_basis_sizes")
+    cut_rows = h5.get_array("cuts_rows")
+    obj_mip = h5.get_scalar("mip_lower_bound")
+    if obj_mip === nothing
+        obj_mip = h5.get_scalar("mip_obj_value")
+    end
+    h5.close()
+
+    # Initialize stats
+    stats_time_convert = 0
+    stats_time_tableau = 0
+    stats_time_gmi = 0
+    all_cuts = nothing
+
+    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, HiGHS.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
+
+    current_basis = nothing
+    for (r, row) in enumerate(cut_rows)
+        stats_time_tableau += @elapsed begin
+            if r == 1 || cut_basis_vars[r, :] != cut_basis_vars[r-1, :]
+                vbb, vnn, cbb, cnn = cut_basis_sizes[r, :]
+                current_basis = Basis(;
+                    var_basic = cut_basis_vars[r, 1:vbb],
+                    var_nonbasic = cut_basis_vars[r, vbb+1:vbb+vnn],
+                    constr_basic = cut_basis_vars[r, vbb+vnn+1:vbb+vnn+cbb],
+                    constr_nonbasic = cut_basis_vars[r, vbb+vnn+cbb+1:vbb+vnn+cbb+cnn],
+                )
+            end
+            tableau = compute_tableau(data_s, current_basis, rows = [row])
+            assert_eq(tableau.lhs * sol_opt_s, tableau.rhs)
+        end
+        stats_time_gmi += @elapsed begin
+            cuts_s = compute_gmi(data_s, tableau)
+            assert_does_not_cut_off(cuts_s, sol_opt_s)
+        end
+        cuts = backwards(transforms, cuts_s)
+        assert_does_not_cut_off(cuts, sol_opt)
+
+        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
+    end
+
+    # Strategy 1: Add all cuts during the first call
+    function cut_callback_1(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
+
+    # Strategy 2: Add violated cuts repeatedly until unable to separate
+    callback_disabled = false
+    function cut_callback_2(cb_data)
+        if callback_disabled
+            return
+        end
+        x = all_variables(model)
+        x_val = callback_value.(cb_data, x)
+        lhs_val = all_cuts.lhs * x_val
+        is_violated = lhs_val .> all_cuts.ub
+        selected_idx = findall(is_violated .== true)
+        selected_cuts = ConstraintSet(
+            lhs=all_cuts.lhs[selected_idx, :],
+            ub=all_cuts.ub[selected_idx],
+            lb=all_cuts.lb[selected_idx],
+        )
+        constrs = build_constraints(model, selected_cuts)
+        if length(constrs) > 0
+            @info "Enforcing $(length(constrs)) cuts..."
+            for c in constrs
+                MOI.submit(model, MOI.UserCut(cb_data), c)
+            end
+        else
+            @info "No violated cuts found. Disabling callback."
+            callback_disabled = true
+        end
+    end
+
+    # Set up cut callback
+    set_attribute(model, MOI.UserCutCallback(), cut_callback_1)
+    # set_attribute(model, MOI.UserCutCallback(), cut_callback_2)
+
+    stats = Dict()
+    stats["ExpertDualGmi: cuts"] = length(all_cuts.lb)
+    stats["ExpertDualGmi: time convert"] = stats_time_convert
+    stats["ExpertDualGmi: time tableau"] = stats_time_tableau
+    stats["ExpertDualGmi: time gmi"] = stats_time_gmi
+    return stats
+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_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
+        vars_to_unique_basis_offset = Dict()
+        unique_basis_vars = nothing
+        unique_basis_sizes = nothing
+        unique_basis_rows = nothing
+    
+        for h5_filename in train_h5
+            h5 = H5File(h5_filename, "r")
+            cut_basis_vars = h5.get_array("cuts_basis_vars")
+            cut_basis_sizes = h5.get_array("cuts_basis_sizes")
+            cut_rows = h5.get_array("cuts_rows")
+            ncuts, nvars = size(cut_basis_vars)
+            if unique_basis_vars === nothing
+                unique_basis_vars = Matrix{Int}(undef, 0, nvars)
+                unique_basis_sizes = Matrix{Int}(undef, 0, 4)
+                unique_basis_rows = Dict{Int,Set{Int}}()                
+            end
+            for i in 1:ncuts
+                vars = cut_basis_vars[i, :]
+                sizes = cut_basis_sizes[i, :]
+                row = cut_rows[i]
+                if vars ∉ keys(vars_to_unique_basis_offset)
+                    offset = size(unique_basis_vars)[1] + 1
+                    vars_to_unique_basis_offset[vars] = offset
+                    unique_basis_vars = [unique_basis_vars; vars']
+                    unique_basis_sizes = [unique_basis_sizes; sizes']
+                    unique_basis_rows[offset] = Set()
+                end
+                offset = vars_to_unique_basis_offset[vars]
+                push!(unique_basis_rows[offset], row)
+            end
+            h5.close()
+        end
+    end
+
+    @timeit "Compute tableaus and cuts" begin
+        all_cuts = nothing
+        for (offset, rows) in unique_basis_rows
+            try
+                vbb, vnn, cbb, cnn = unique_basis_sizes[offset, :]
+                current_basis = Basis(;
+                    var_basic = unique_basis_vars[offset, 1:vbb],
+                    var_nonbasic = unique_basis_vars[offset, vbb+1:vbb+vnn],
+                    constr_basic = unique_basis_vars[offset, vbb+vnn+1:vbb+vnn+cbb],
+                    constr_nonbasic = unique_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 Class1
+        function fit(_, _) end
+        function before_mip(self, test_h5, model, stats)
+            ExpertDualGmiComponent_before_mip(test_h5, model.inner, stats)
+        end
+    end
+    copy!(ExpertDualGmiComponent, Class1)
+
+    @pydef mutable struct Class2
+        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, Class2)
+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,6 +6,8 @@ module MIPLearn
 
 using PyCall
 using SparseArrays
+using PrecompileTools: @setup_workload, @compile_workload
+
 
 include("collectors.jl")
 include("components.jl")
@@ -32,4 +34,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/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
 
 # -----------------------------------------------------------------------------
@@ -71,13 +69,14 @@ function submit(model::JuMP.Model, constr)
 end
 
 function _extract_after_load(model::JuMP.Model, h5)
+    @info "_extract_after_load"
     if JuMP.objective_sense(model) == MOI.MIN_SENSE
         h5.put_scalar("static_sense", "min")
     else
         h5.put_scalar("static_sense", "max")
     end
-    _extract_after_load_vars(model, h5)
-    _extract_after_load_constrs(model, h5)
+    @time _extract_after_load_vars(model, h5)
+    @time _extract_after_load_constrs(model, h5)
 end
 
 function _extract_after_load_vars(model::JuMP.Model, h5)
@@ -155,10 +154,11 @@ function _extract_after_load_constrs(model::JuMP.Model, h5)
 end
 
 function _extract_after_lp(model::JuMP.Model, h5)
+    @info "_extract_after_lp"
     h5.put_scalar("lp_wallclock_time", solve_time(model))
     h5.put_scalar("lp_obj_value", objective_value(model))
-    _extract_after_lp_vars(model, h5)
-    _extract_after_lp_constrs(model, h5)
+    @time _extract_after_lp_vars(model, h5)
+    @time _extract_after_lp_constrs(model, h5)
 end
 
 function _extract_after_lp_vars(model::JuMP.Model, h5)
@@ -184,46 +184,46 @@ function _extract_after_lp_vars(model::JuMP.Model, h5)
     end
     h5.put_array("lp_var_basis_status", to_str_array(basis_status))
 
-    # Sensitivity analysis
-    obj_coeffs = h5.get_array("static_var_obj_coeffs")
-    sensitivity_report = lp_sensitivity_report(model)
-    sa_obj_down, sa_obj_up = Float64[], Float64[]
-    sa_lb_down, sa_lb_up = Float64[], Float64[]
-    sa_ub_down, sa_ub_up = Float64[], Float64[]
-    for (i, v) in enumerate(vars)
-        # Objective function
-        (delta_down, delta_up) = sensitivity_report[v]
-        push!(sa_obj_down, delta_down + obj_coeffs[i])
-        push!(sa_obj_up, delta_up + obj_coeffs[i])
+    # # Sensitivity analysis
+    # obj_coeffs = h5.get_array("static_var_obj_coeffs")
+    # sensitivity_report = lp_sensitivity_report(model)
+    # sa_obj_down, sa_obj_up = Float64[], Float64[]
+    # sa_lb_down, sa_lb_up = Float64[], Float64[]
+    # sa_ub_down, sa_ub_up = Float64[], Float64[]
+    # for (i, v) in enumerate(vars)
+    #     # Objective function
+    #     (delta_down, delta_up) = sensitivity_report[v]
+    #     push!(sa_obj_down, delta_down + obj_coeffs[i])
+    #     push!(sa_obj_up, delta_up + obj_coeffs[i])
 
-        # Lower bound
-        if has_lower_bound(v)
-            constr = LowerBoundRef(v)
-            (delta_down, delta_up) = sensitivity_report[constr]
-            push!(sa_lb_down, lower_bound(v) + delta_down)
-            push!(sa_lb_up, lower_bound(v) + delta_up)
-        else
-            push!(sa_lb_down, -Inf)
-            push!(sa_lb_up, -Inf)
-        end
+    #     # Lower bound
+    #     if has_lower_bound(v)
+    #         constr = LowerBoundRef(v)
+    #         (delta_down, delta_up) = sensitivity_report[constr]
+    #         push!(sa_lb_down, lower_bound(v) + delta_down)
+    #         push!(sa_lb_up, lower_bound(v) + delta_up)
+    #     else
+    #         push!(sa_lb_down, -Inf)
+    #         push!(sa_lb_up, -Inf)
+    #     end
 
-        # Upper bound
-        if has_upper_bound(v)
-            constr = JuMP.UpperBoundRef(v)
-            (delta_down, delta_up) = sensitivity_report[constr]
-            push!(sa_ub_down, upper_bound(v) + delta_down)
-            push!(sa_ub_up, upper_bound(v) + delta_up)
-        else
-            push!(sa_ub_down, Inf)
-            push!(sa_ub_up, Inf)
-        end
-    end
-    h5.put_array("lp_var_sa_obj_up", sa_obj_up)
-    h5.put_array("lp_var_sa_obj_down", sa_obj_down)
-    h5.put_array("lp_var_sa_ub_up", sa_ub_up)
-    h5.put_array("lp_var_sa_ub_down", sa_ub_down)
-    h5.put_array("lp_var_sa_lb_up", sa_lb_up)
-    h5.put_array("lp_var_sa_lb_down", sa_lb_down)
+    #     # Upper bound
+    #     if has_upper_bound(v)
+    #         constr = JuMP.UpperBoundRef(v)
+    #         (delta_down, delta_up) = sensitivity_report[constr]
+    #         push!(sa_ub_down, upper_bound(v) + delta_down)
+    #         push!(sa_ub_up, upper_bound(v) + delta_up)
+    #     else
+    #         push!(sa_ub_down, Inf)
+    #         push!(sa_ub_up, Inf)
+    #     end
+    # end
+    # h5.put_array("lp_var_sa_obj_up", sa_obj_up)
+    # h5.put_array("lp_var_sa_obj_down", sa_obj_down)
+    # h5.put_array("lp_var_sa_ub_up", sa_ub_up)
+    # h5.put_array("lp_var_sa_ub_down", sa_ub_down)
+    # h5.put_array("lp_var_sa_lb_up", sa_lb_up)
+    # h5.put_array("lp_var_sa_lb_down", sa_lb_down)
 end
 
 
@@ -239,7 +239,7 @@ function _extract_after_lp_constrs(model::JuMP.Model, h5)
     duals = Float64[]
     basis_status = []
     constr_idx = 1
-    sensitivity_report = lp_sensitivity_report(model)
+    # sensitivity_report = lp_sensitivity_report(model)
     for (ftype, stype) in JuMP.list_of_constraint_types(model)
         for constr in JuMP.all_constraints(model, ftype, stype)
             length(JuMP.name(constr)) > 0 || continue
@@ -257,21 +257,22 @@ function _extract_after_lp_constrs(model::JuMP.Model, h5)
                 error("Unknown basis status: $b")
             end
 
-            # Sensitivity analysis
-            (delta_down, delta_up) = sensitivity_report[constr]
-            push!(sa_rhs_down, rhs[constr_idx] + delta_down)
-            push!(sa_rhs_up, rhs[constr_idx] + delta_up)
+            # # Sensitivity analysis
+            # (delta_down, delta_up) = sensitivity_report[constr]
+            # push!(sa_rhs_down, rhs[constr_idx] + delta_down)
+            # push!(sa_rhs_up, rhs[constr_idx] + delta_up)
 
             constr_idx += 1
         end
     end
     h5.put_array("lp_constr_dual_values", duals)
     h5.put_array("lp_constr_basis_status", to_str_array(basis_status))
-    h5.put_array("lp_constr_sa_rhs_up", sa_rhs_up)
-    h5.put_array("lp_constr_sa_rhs_down", sa_rhs_down)
+    # h5.put_array("lp_constr_sa_rhs_up", sa_rhs_up)
+    # h5.put_array("lp_constr_sa_rhs_down", sa_rhs_down)
 end
 
 function _extract_after_mip(model::JuMP.Model, h5)
+    @info "_extract_after_mip"
     h5.put_scalar("mip_obj_value", objective_value(model))
     h5.put_scalar("mip_obj_bound", objective_bound(model))
     h5.put_scalar("mip_wallclock_time", solve_time(model))
@@ -297,9 +298,11 @@ 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)
+        if isfinite(var_values[i])
             fix(var, var_values[i], force=true)
         end
     end
+end
 
 function _optimize(model::JuMP.Model)
     # Set up cut callbacks

test/Project.toml

@@ -4,7 +4,9 @@ authors = ["Alinson S. Xavier <git@axavier.org>"]
 version = "0.1.0"
 
 [deps]
+CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 Clp = "e2554f3b-3117-50c0-817c-e040a3ddf72d"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
 Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
@@ -17,6 +19,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

@@ -10,6 +10,9 @@ using Test
 using MIPLearn.BB
 using MIPLearn
 
+using CSV
+using DataFrames
+
 basepath = @__DIR__
 
 function bb_run(optimizer_name, optimizer; large=true)
@@ -86,11 +89,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 +97,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,
@@ -132,3 +131,67 @@ function test_bb()
     @time bb_run("HiGHS", optimizer_with_attributes(HiGHS.Optimizer))
     # @time bb_run("CPLEX", optimizer_with_attributes(CPLEX.Optimizer, "CPXPARAM_Threads" => 1))
 end
+
+function test_bb_replay()
+    rule_sb = BB.StrongBranching()
+    rule_rb = BB.ReliabilityBranching()
+    optimizer = optimizer_with_attributes(HiGHS.Optimizer)
+    filenames = [replace(f, ".h5" => "") for f in glob("test/fixtures/stab/*.h5")]
+    results_filename = "tmp.csv"
+
+    lk = ReentrantLock()
+    results = []
+
+    function push_result(r)
+        lock(lk) do
+            push!(results, r)
+            df = DataFrame()
+            for row in results
+                push!(df, row, cols=:union)
+            end
+            CSV.write(results_filename, df)
+        end
+    end
+
+    function solve(filename; replay=nothing, skip=false, rule)
+        has_replay = (replay !== nothing)
+        h5 = H5File("$filename.h5", "r")
+        mip_obj_bound = h5.get_scalar("mip_obj_bound")
+        @show filename
+        @show has_replay
+        h5.file.close()
+        mip = BB.init(optimizer)
+        BB.read!(mip, "$filename.mps.gz")
+        time_solve = @elapsed begin
+            pool, replay = BB.solve!(
+                mip,
+                initial_primal_bound=mip_obj_bound,
+                print_interval=100,
+                node_limit=1_000,
+                branch_rule=rule,
+                replay=replay,
+            )
+        end
+        if !skip
+            push_result(
+                Dict(
+                    "Filename" => filename,
+                    "Replay?" => has_replay,
+                    "Solve time (s)" => time_solve,
+                    "Relative MIP gap (%)" => round(pool.gap * 100, digits=3)
+                )
+            )
+        end
+        return replay
+    end
+
+    # Solve reference instance
+    replay = solve(filenames[1], skip=true, rule=rule_sb)
+
+    # Solve perturbations
+    for i in 2:6
+        solve(filenames[i], rule=rule_rb, replay=nothing)
+        solve(filenames[i], rule=rule_rb, replay=deepcopy(replay))
+    end
+    return
+end

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,6 +19,8 @@ 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")

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