FisSal2011: Adjust constants

Project.toml

@@ -1,11 +1,12 @@
 name = "MIPLearn"
 uuid = "2b1277c3-b477-4c49-a15e-7ba350325c68"
 authors = ["Alinson S Xavier <git@axavier.org>"]
-version = "0.4.0"
+version = "0.4.2"
 
 [deps]
 Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+Gurobi = "2e9cd046-0924-5485-92f1-d5272153d98b"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
@@ -28,6 +29,7 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 [compat]
 Conda = "1"
 DataStructures = "0.18"
+Gurobi = "1.7.5"
 HDF5 = "0.16"
 HiGHS = "1"
 JLD2 = "0.4"
@@ -36,8 +38,10 @@ JuMP = "1"
 KLU = "0.4"
 MathOptInterface = "1"
 OrderedCollections = "1"
+PrecompileTools = "1"
 PyCall = "1"
 Requires = "1"
+SCIP = "0.12"
 Statistics = "1"
 TimerOutputs = "0.5"
 julia = "1"

deps/build.jl

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

src/BB/log.jl

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

src/BB/lp.jl

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

src/Cuts/tableau/gmi.jl

@@ -185,7 +185,7 @@ function collect_gmi(
     )
 end
 
-function select_gmi_rows(data, basis, x; max_rows = 10, atol = 1e-4)
+function select_gmi_rows(data, basis, x; max_rows = 10, atol = 0.001)
     candidate_rows = [
         r for r = 1:length(basis.var_basic) if (
             (data.var_types[basis.var_basic[r]] != 'C') &&
@@ -199,23 +199,82 @@ function select_gmi_rows(data, basis, x; max_rows = 10, atol = 1e-4)
     return [candidate_rows[perm[i]] for i = 1:min(length(perm), max_rows)]
 end
 
+# function compute_gmi(data::ProblemData, tableau::Tableau)::ConstraintSet
+#     @timeit "Initialization" begin
+#         nrows, ncols = size(tableau.lhs)
+#         ub = Float64[Inf for _ = 1:nrows]
+#         lb = Float64[0.999 for _ = 1:nrows]
+#         tableau_I, tableau_J, tableau_V = findnz(tableau.lhs)
+#         lhs_I = Int[]
+#         lhs_J = Int[]
+#         lhs_V = Float64[]
+#     end
+#     @timeit "Compute coefficients" begin
+#         for k = 1:nnz(tableau.lhs)
+#             i::Int = tableau_I[k]
+#             j::Int = tableau_J[k]
+#             v::Float64 = 0.0
+#             frac_alpha_j = frac(tableau_V[k])
+#             alpha_j = tableau_V[k]
+#             beta = frac(tableau.rhs[i])
+#             if data.var_types[j] == 'C'
+#                 if alpha_j >= 0
+#                     v = alpha_j / beta
+#                 else
+#                     v = -alpha_j / (1 - beta)
+#                 end
+#             else
+#                 if frac_alpha_j < beta
+#                     v = frac_alpha_j / beta
+#                 else
+#                     v = (1 - frac_alpha_j) / (1 - beta)
+#                 end
+#             end
+#             if abs(v) > 1e-8
+#                 push!(lhs_I, i)
+#                 push!(lhs_J, tableau_J[k])
+#                 push!(lhs_V, v)
+#             end
+#         end
+#     end
+#     @timeit "Convert to ConstraintSet" begin
+#         lhs = sparse(lhs_I, lhs_J, lhs_V, nrows, ncols)
+#         cs = ConstraintSet(; lhs, ub, lb)
+#     end
+#     return cs
+# end
+
 function compute_gmi(data::ProblemData, tableau::Tableau)::ConstraintSet
-    nrows, ncols = size(tableau.lhs)
-    ub = Float64[Inf for _ = 1:nrows]
-    lb = Float64[0.9999 for _ = 1:nrows]
-    tableau_I, tableau_J, tableau_V = findnz(tableau.lhs)
-    lhs_I = Int[]
-    lhs_J = Int[]
-    lhs_V = Float64[]
+    @timeit "Initialization" begin
+        nrows::Int, ncols::Int = size(tableau.lhs)
+        var_types::Vector{Char} = data.var_types
+        tableau_rhs::Vector{Float64} = tableau.rhs
+        tableau_I::Vector{Int}, tableau_J::Vector{Int}, tableau_V::Vector{Float64} = findnz(tableau.lhs)
+    end
+
+    @timeit "Pre-allocation" begin
+        cut_ub::Vector{Float64} = fill(Inf, nrows)
+        cut_lb::Vector{Float64} = fill(0.999, nrows)
+        nnz_tableau::Int = length(tableau_I)
+        cut_lhs_I::Vector{Int} = Vector{Int}(undef, nnz_tableau)
+        cut_lhs_J::Vector{Int} = Vector{Int}(undef, nnz_tableau)
+        cut_lhs_V::Vector{Float64} = Vector{Float64}(undef, nnz_tableau)
+        cut_hash::Vector{UInt64} = zeros(UInt64, nrows)
+        nnz_count::Int = 0
+    end
+    
     @timeit "Compute coefficients" begin
-        for k = 1:nnz(tableau.lhs)
+        @inbounds for k = 1:nnz_tableau
             i::Int = tableau_I[k]
             j::Int = tableau_J[k]
-            v::Float64 = 0.0
-            frac_alpha_j = frac(tableau_V[k])
-            alpha_j = tableau_V[k]
-            beta = frac(tableau.rhs[i])
-            if data.var_types[j] == 'C'
+            alpha_j::Float64 = tableau_V[k]
+            frac_alpha_j::Float64 = alpha_j - floor(alpha_j)
+            beta_i::Float64 = tableau_rhs[i]
+            beta::Float64 = beta_i - floor(beta_i)
+            v::Float64 = 0
+            
+            # Compute coefficient
+            if var_types[j] == 'C'
                 if alpha_j >= 0
                     v = alpha_j / beta
                 else
@@ -228,16 +287,34 @@ function compute_gmi(data::ProblemData, tableau::Tableau)::ConstraintSet
                     v = (1 - frac_alpha_j) / (1 - beta)
                 end
             end
+            
+            # Store if significant
             if abs(v) > 1e-8
-                push!(lhs_I, i)
-                push!(lhs_J, tableau_J[k])
-                push!(lhs_V, v)
+                nnz_count += 1
+                cut_lhs_I[nnz_count] = i
+                cut_lhs_J[nnz_count] = j
+                cut_lhs_V[nnz_count] = v
+                cut_hash[i] = hash(j, cut_hash[i])
+                cut_hash[i] = hash(v, cut_hash[i])
             end
         end
-        lhs = sparse(lhs_I, lhs_J, lhs_V, nrows, ncols)
     end
-    return ConstraintSet(; lhs, ub, lb)
+
+    @timeit "Resize arrays to actual size" begin
+        resize!(cut_lhs_I, nnz_count)
+        resize!(cut_lhs_J, nnz_count)
+        resize!(cut_lhs_V, nnz_count)
+    end
+
+    # TODO: Build cut in compressed row format instead of converting
+    @timeit "Convert to ConstraintSet" begin
+        cut_lhs::SparseMatrixCSC = sparse(cut_lhs_I, cut_lhs_J, cut_lhs_V, nrows, ncols)
+        cs::ConstraintSet = ConstraintSet(; lhs=cut_lhs, ub=cut_ub, lb=cut_lb, hash=cut_hash)
+    end
+    
+    return cs
 end
 
+
 export compute_gmi,
     frac, select_gmi_rows, assert_cuts_off, assert_does_not_cut_off, collect_gmi

src/Cuts/tableau/numerics.jl

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

src/Cuts/tableau/structs.jl

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

src/Cuts/tableau/tableau.jl

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

src/Cuts/tableau/transform.jl

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

src/MIPLearn.jl

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

@@ -89,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)
@@ -143,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)
@@ -282,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]

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