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"
@@ -41,3 +41,5 @@ 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",
@@ -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)
+    if (pool.processed % print_interval == 0) || isempty(pool.pending) || primal_update
         if isempty(node.branch_vars)
             branch_var_name = "---"
             branch_lb = "---"
@@ -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/BB/nodepool.jl

@@ -8,10 +8,10 @@ import Base.Threads: threadid
 
 function take(
     pool::NodePool;
-    suggestions::Array{Node}=[],
+    suggestions::Array{Node} = [],
     time_remaining::Float64,
     gap_limit::Float64,
-    node_limit::Int
+    node_limit::Int,
 )::Union{Symbol,Node}
     t = threadid()
     lock(pool.lock) do
@@ -53,8 +53,8 @@ function offer(
     pool::NodePool;
     parent_node::Union{Nothing,Node},
     child_nodes::Vector{Node},
-    time_elapsed::Float64=0.0,
-    print_interval::Int=100
+    time_elapsed::Float64 = 0.0,
+    print_interval::Int = 100,
 )::Nothing
     lock(pool.lock) do
         primal_update = false
@@ -101,32 +101,30 @@ 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
-                _update_var_history(
-                    pool=pool,
-                    var=branch_var,
-                    x=x,
-                    obj_change_down=obj_change_down,
-                    obj_change_up=obj_change_up,
-                )
-                # Update global history
-                pool.history.avg_pseudocost_up =
-                    mean(vh.pseudocost_up for vh in values(pool.var_history))
-                pool.history.avg_pseudocost_down =
-                    mean(vh.pseudocost_down for vh in values(pool.var_history))
-            end
+            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
+            _update_var_history(
+                pool = pool,
+                var = branch_var,
+                x = x,
+                obj_change_down = obj_change_down,
+                obj_change_up = obj_change_up,
+            )
+            # Update global history
+            pool.history.avg_pseudocost_up =
+                mean(vh.pseudocost_up for vh in values(pool.var_history))
+            pool.history.avg_pseudocost_down =
+                mean(vh.pseudocost_down for vh in values(pool.var_history))
         end
 
         for node in child_nodes
             print_progress(
                 pool,
                 node,
-                time_elapsed=time_elapsed,
-                print_interval=print_interval,
-                primal_update=isfinite(node.obj) && isempty(node.fractional_variables),
+                time_elapsed = time_elapsed,
+                print_interval = print_interval,
+                primal_update = isfinite(node.obj) && isempty(node.fractional_variables),
             )
         end
     end
@@ -138,7 +136,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

@@ -10,22 +10,16 @@ import ..H5File
 
 function solve!(
     mip::MIP;
-    time_limit::Float64=Inf,
-    node_limit::Int=typemax(Int),
-    gap_limit::Float64=1e-4,
-    print_interval::Int=5,
-    initial_primal_bound::Float64=Inf,
-    branch_rule::VariableBranchingRule=ReliabilityBranching(),
-    enable_plunging=true,
-    replay=nothing
-)::Tuple{NodePool,ReplayInfo}
-
-    if replay === nothing
-        replay = ReplayInfo()
-    end
-
+    time_limit::Float64 = Inf,
+    node_limit::Int = typemax(Int),
+    gap_limit::Float64 = 1e-4,
+    print_interval::Int = 5,
+    initial_primal_bound::Float64 = Inf,
+    branch_rule::VariableBranchingRule = ReliabilityBranching(),
+    enable_plunging = true,
+)::NodePool
     time_initial = time()
-    pool = NodePool(mip=mip, next_index=replay.next_index)
+    pool = NodePool(mip = mip)
     pool.primal_bound = initial_primal_bound
 
     root_node = _create_node(mip)
@@ -40,9 +34,9 @@ function solve!(
 
     offer(
         pool,
-        parent_node=nothing,
-        child_nodes=[root_node],
-        print_interval=print_interval,
+        parent_node = nothing,
+        child_nodes = [root_node],
+        print_interval = print_interval,
     )
     @threads for t = 1:nthreads()
         child_one, child_zero, suggestions = nothing, nothing, Node[]
@@ -53,10 +47,10 @@ function solve!(
             end
             node = take(
                 pool,
-                suggestions=suggestions,
-                time_remaining=time_limit - time_elapsed,
-                node_limit=node_limit,
-                gap_limit=gap_limit,
+                suggestions = suggestions,
+                time_remaining = time_limit - time_elapsed,
+                node_limit = node_limit,
+                gap_limit = gap_limit,
             )
             if node == :END
                 break
@@ -70,24 +64,9 @@ 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 branching variable
+                ids = generate_indices(pool, 2)
+                branch_var = find_branching_var(branch_rule, node, pool)
 
                 # Find current variable lower and upper bounds
                 offset = findfirst(isequal(branch_var), mip.int_vars)
@@ -100,43 +79,46 @@ 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],
-                    parent=node,
-                    branch_var=branch_var,
-                    branch_var_lb=var_lb,
-                    branch_var_ub=floor(var_value),
+                    index = ids[2],
+                    parent = node,
+                    branch_var = branch_var,
+                    branch_var_lb = var_lb,
+                    branch_var_ub = floor(var_value),
                 )
                 child_one = _create_node(
                     mip,
-                    index=ids[1],
-                    parent=node,
-                    branch_var=branch_var,
-                    branch_var_lb=ceil(var_value),
-                    branch_var_ub=var_ub,
+                    index = ids[1],
+                    parent = node,
+                    branch_var = branch_var,
+                    branch_var_lb = ceil(var_value),
+                    branch_var_ub = var_ub,
                 )
                 offer(
                     pool,
-                    parent_node=node,
-                    child_nodes=[child_one, child_zero],
-                    time_elapsed=time_elapsed,
-                    print_interval=print_interval,
+                    parent_node = node,
+                    child_nodes = [child_one, child_zero],
+                    time_elapsed = time_elapsed,
+                    print_interval = print_interval,
                 )
             end
         end
     end
-    replay.next_index = pool.next_index
-    return pool, replay
+    return pool
 end
 
 function _create_node(
     mip;
-    index::Int=0,
-    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
+    index::Int = 0,
+    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,
 )::Node
     if parent === nothing
         branch_vars = Variable[]

src/BB/structs.jl

@@ -72,14 +72,3 @@ 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/tableau/gmi_dual.jl

@@ -8,6 +8,8 @@ using HiGHS
 using Random
 using DataStructures
 
+import ..H5FieldsExtractor
+
 global ExpertDualGmiComponent = PyNULL()
 global KnnDualGmiComponent = PyNULL()
 
@@ -24,8 +26,10 @@ function collect_gmi_dual(
     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")
@@ -205,6 +209,12 @@ function collect_gmi_dual(
                 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
@@ -253,138 +263,6 @@ function collect_gmi_dual(
     )
 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)
@@ -441,6 +319,58 @@ function _dualgmi_features(h5_filename, extractor)
     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)
@@ -448,54 +378,71 @@ function _dualgmi_generate(train_h5, model)
     @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
-    
+        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
-            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()
+            @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
-                offset = vars_to_unique_basis_offset[vars]
-                push!(unique_basis_rows[offset], row)
+                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
-            h5.close()
         end
     end
-
     @timeit "Compute tableaus and cuts" begin
         all_cuts = nothing
-        for (offset, rows) in unique_basis_rows
+        nbasis = length(combined_cut_rows)
+        for offset in 1:nbasis
+            rows = combined_cut_rows[offset]
             try
-                vbb, vnn, cbb, cnn = unique_basis_sizes[offset, :]
+                vbb, vnn, cbb, cnn = combined_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],
+                    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)
@@ -599,15 +546,7 @@ function KnnDualGmiComponent_before_mip(data::_KnnDualGmiData, test_h5, model, _
 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
+    @pydef mutable struct KnnDualGmiComponentPy
         function __init__(self; extractor, k = 3, strategy = "near")
             self.data = _KnnDualGmiData(; extractor, k, strategy)
         end
@@ -618,7 +557,23 @@ function __init_gmi_dual__()
             return @time KnnDualGmiComponent_before_mip(self.data, test_h5, model.inner, stats)
         end
     end
-    copy!(KnnDualGmiComponent, Class2)
+    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/MIPLearn.jl

@@ -13,6 +13,7 @@ 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")
@@ -24,6 +25,7 @@ function __init__()
     __init_components__()
     __init_extractors__()
     __init_io__()
+    __init_problems_maxcut__()
     __init_problems_setcover__()
     __init_problems_stab__()
     __init_problems_tsp__()
@@ -37,48 +39,48 @@ 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
+# 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
+# @setup_workload begin
+#     using SCIP
+#     using HiGHS
+#     using MIPLearn.Cuts
+#     using PrecompileTools: @setup_workload, @compile_workload
 
-    __init__()
-    Cuts.__init__()
+#     __init__()
+#     Cuts.__init__()
 
-    @compile_workload begin
-        __precompile_cuts__()
-    end
-end
+#     @compile_workload begin
+#         __precompile_cuts__()
+#     end
+# end
 
 end # module

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/solvers/jump.jl

@@ -69,14 +69,13 @@ 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
-    @time _extract_after_load_vars(model, h5)
-    @time _extract_after_load_constrs(model, h5)
+    _extract_after_load_vars(model, h5)
+    _extract_after_load_constrs(model, h5)
 end
 
 function _extract_after_load_vars(model::JuMP.Model, h5)
@@ -90,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)
@@ -144,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)
@@ -154,11 +166,10 @@ 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))
-    @time _extract_after_lp_vars(model, h5)
-    @time _extract_after_lp_constrs(model, h5)
+    _extract_after_lp_vars(model, h5)
+    _extract_after_lp_constrs(model, h5)
 end
 
 function _extract_after_lp_vars(model::JuMP.Model, h5)
@@ -184,46 +195,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 +250,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,22 +268,21 @@ 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))
@@ -285,9 +295,11 @@ function _extract_after_mip(model::JuMP.Model, h5)
 
     # Slacks
     lhs = h5.get_sparse("static_constr_lhs")
-    rhs = h5.get_array("static_constr_rhs")
-    slacks = abs.(lhs * x - rhs)
-    h5.put_array("mip_constr_slacks", slacks)
+    if lhs !== nothing
+        rhs = h5.get_array("static_constr_rhs")
+        slacks = abs.(lhs * x - rhs)
+        h5.put_array("mip_constr_slacks", slacks)
+    end
 
     # Cuts and lazy constraints
     ext = model.ext[:miplearn]
@@ -298,9 +310,7 @@ end
 function _fix_variables(model::JuMP.Model, var_names, var_values, stats)
     vars = [variable_by_name(model, v) for v in var_names]
     for (i, var) in enumerate(vars)
-        if isfinite(var_values[i])
-            fix(var, var_values[i], force=true)
-        end
+        fix(var, var_values[i], force = true)
     end
 end
 
@@ -417,7 +427,7 @@ function __init_solvers_jump__()
             constrs_lhs,
             constrs_sense,
             constrs_rhs,
-            stats=nothing,
+            stats = nothing,
         ) = _add_constrs(
             self.inner,
             from_str_array(var_names),
@@ -433,14 +443,14 @@ function __init_solvers_jump__()
 
         extract_after_mip(self, h5) = _extract_after_mip(self.inner, h5)
 
-        fix_variables(self, var_names, var_values, stats=nothing) =
+        fix_variables(self, var_names, var_values, stats = nothing) =
             _fix_variables(self.inner, from_str_array(var_names), var_values, stats)
 
         optimize(self) = _optimize(self.inner)
 
         relax(self) = Class(_relax(self.inner))
 
-        set_warm_starts(self, var_names, var_values, stats=nothing) =
+        set_warm_starts(self, var_names, var_values, stats = nothing) =
             _set_warm_starts(self.inner, from_str_array(var_names), var_values, stats)
 
         write(self, filename) = _write(self.inner, filename)

test/Project.toml

@@ -4,9 +4,7 @@ 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"

test/src/BB/test_bb.jl

@@ -10,12 +10,9 @@ using Test
 using MIPLearn.BB
 using MIPLearn
 
-using CSV
-using DataFrames
-
 basepath = @__DIR__
 
-function bb_run(optimizer_name, optimizer; large=true)
+function bb_run(optimizer_name, optimizer; large = true)
     @testset "Solve ($optimizer_name)" begin
         @testset "interface" begin
             filename = "$FIXTURES/danoint.mps.gz"
@@ -28,7 +25,7 @@ function bb_run(optimizer_name, optimizer; large=true)
 
             status, obj = BB.solve_relaxation!(mip)
             @test status == :Optimal
-            @test round(obj, digits=6) == 62.637280
+            @test round(obj, digits = 6) == 62.637280
 
             @test BB.name(mip, mip.int_vars[1]) == "xab"
             @test BB.name(mip, mip.int_vars[2]) == "xac"
@@ -38,26 +35,26 @@ function bb_run(optimizer_name, optimizer; large=true)
             @test mip.int_vars_ub[1] == 1.0
 
             vals = BB.values(mip, mip.int_vars)
-            @test round(vals[1], digits=6) == 0.046933
-            @test round(vals[2], digits=6) == 0.000841
-            @test round(vals[3], digits=6) == 0.248696
+            @test round(vals[1], digits = 6) == 0.046933
+            @test round(vals[2], digits = 6) == 0.000841
+            @test round(vals[3], digits = 6) == 0.248696
 
             # Probe (up and down are feasible)
             probe_up, probe_down = BB.probe(mip, mip.int_vars[1], 0.5, 0.0, 1.0, 1_000_000)
-            @test round(probe_down, digits=6) == 62.690000
-            @test round(probe_up, digits=6) == 62.714100
+            @test round(probe_down, digits = 6) == 62.690000
+            @test round(probe_up, digits = 6) == 62.714100
 
             # Fix one variable to zero
             BB.set_bounds!(mip, mip.int_vars[1:1], [0.0], [0.0])
             status, obj = BB.solve_relaxation!(mip)
             @test status == :Optimal
-            @test round(obj, digits=6) == 62.690000
+            @test round(obj, digits = 6) == 62.690000
 
             # Fix one variable to one and another variable variable to zero
             BB.set_bounds!(mip, mip.int_vars[1:2], [1.0, 0.0], [1.0, 0.0])
             status, obj = BB.solve_relaxation!(mip)
             @test status == :Optimal
-            @test round(obj, digits=6) == 62.714777
+            @test round(obj, digits = 6) == 62.714777
 
             # Fix all binary variables to one, making problem infeasible
             N = length(mip.int_vars)
@@ -71,7 +68,7 @@ function bb_run(optimizer_name, optimizer; large=true)
             BB.set_bounds!(mip, mip.int_vars, zeros(N), ones(N))
             status, obj = BB.solve_relaxation!(mip)
             @test status == :Optimal
-            @test round(obj, digits=6) == 62.637280
+            @test round(obj, digits = 6) == 62.637280
         end
 
         @testset "varbranch" begin
@@ -85,8 +82,8 @@ function bb_run(optimizer_name, optimizer; large=true)
                     BB.StrongBranching(),
                     BB.ReliabilityBranching(),
                     BB.HybridBranching(),
-                    BB.StrongBranching(aggregation=:min),
-                    BB.ReliabilityBranching(aggregation=:min, collect=true),
+                    BB.StrongBranching(aggregation = :min),
+                    BB.ReliabilityBranching(aggregation = :min, collect = true),
                 ]
                     h5 = H5File("$FIXTURES/$instance.h5")
                     mip_obj_bound = h5.get_scalar("mip_obj_bound")
@@ -107,13 +104,13 @@ function bb_run(optimizer_name, optimizer; large=true)
         end
 
         @testset "collect" begin
-            rule = BB.ReliabilityBranching(collect=true)
+            rule = BB.ReliabilityBranching(collect = true)
             BB.collect!(
                 optimizer,
                 "$FIXTURES/bell5.mps.gz",
-                node_limit=100,
-                print_interval=10,
-                branch_rule=rule,
+                node_limit = 100,
+                print_interval = 10,
+                branch_rule = rule,
             )
             n_sb = rule.stats.num_strong_branch_calls
             h5 = H5File("$FIXTURES/bell5.h5")
@@ -131,67 +128,3 @@ 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/MIPLearnT.jl

@@ -24,6 +24,7 @@ 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")
@@ -37,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/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