FisSal2011: Adjust constants

Project.toml

@@ -6,6 +6,7 @@ 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,10 +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"
-PrecompileTools = "1"
-SCIP = "0.12"

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)
+    @timeit "Initialization" begin
-    ub = Float64[Inf for _ = 1:nrows]
+        nrows::Int, ncols::Int = size(tableau.lhs)
-    lb = Float64[0.9999 for _ = 1:nrows]
+        var_types::Vector{Char} = data.var_types
-    tableau_I, tableau_J, tableau_V = findnz(tableau.lhs)
+        tableau_rhs::Vector{Float64} = tableau.rhs
-    lhs_I = Int[]
+        tableau_I::Vector{Int}, tableau_J::Vector{Int}, tableau_V::Vector{Float64} = findnz(tableau.lhs)
-    lhs_J = Int[]
+    end
-    lhs_V = Float64[]
+
+    @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
+            alpha_j::Float64 = tableau_V[k]
-            frac_alpha_j = frac(tableau_V[k])
+            frac_alpha_j::Float64 = alpha_j - floor(alpha_j)
-            alpha_j = tableau_V[k]
+            beta_i::Float64 = tableau_rhs[i]
-            beta = frac(tableau.rhs[i])
+            beta::Float64 = beta_i - floor(beta_i)
-            if data.var_types[j] == 'C'
+            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)
+                nnz_count += 1
-                push!(lhs_J, tableau_J[k])
+                cut_lhs_I[nnz_count] = i
-                push!(lhs_V, v)
+                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/gmi_dual.jl

@@ -7,6 +7,7 @@ using JuMP
 using HiGHS
 using Random
 using DataStructures
+using Statistics
 
 import ..H5FieldsExtractor
 
@@ -25,7 +26,7 @@ function collect_gmi_dual(
     mps_filename;
     optimizer,
     max_rounds = 10,
-    max_cuts_per_round = 500,
+    max_cuts_per_round = 1_000_000,
     time_limit = 3_600,
 )
     reset_timer!()
@@ -263,6 +264,689 @@ function collect_gmi_dual(
     )
 end
 
+function collect_gmi_FisSal2011(
+    mps_filename;
+    interval_print_sec = 1,
+    max_cuts_per_round = 1_000,
+    max_pool_size_mb = 1024,
+    optimizer,
+    silent_solver = true,
+    time_limit = 300,
+    variant = :miplearn,
+    verify_cuts = true,
+)
+    variant in [:subg, :hybr, :fast, :faster, :miplearn] || error("unknown variant: $variant")
+    if variant == :subg
+        max_rounds = 10_000
+        interval_large_lp = 10_000
+        interval_read_tableau = 10
+    elseif variant == :hybr
+        max_rounds = 10_000
+        interval_large_lp = 1_000
+        interval_read_tableau = 10
+    elseif variant == :fast
+        max_rounds = 1_000
+        interval_large_lp = 100
+        interval_read_tableau = 1
+    elseif variant == :faster
+        max_rounds = 500
+        interval_large_lp = 50
+        interval_read_tableau = 1
+    elseif variant == :miplearn
+        max_rounds = 1_000_000
+        interval_large_lp = 100
+        interval_read_tableau = 1
+    end
+    gapcl_best_patience = 2 * interval_large_lp + 5
+
+    reset_timer!()
+    initial_time = time()
+
+    @timeit "Read H5" begin
+        h5_filename = replace(mps_filename, ".mps.gz" => ".h5")
+        h5 = H5File(h5_filename, "r")
+        sol_opt_dict = Dict(
+            zip(
+                h5.get_array("static_var_names"),
+                convert(Array{Float64}, h5.get_array("mip_var_values")),
+            ),
+        )
+        obj_mip = h5.get_scalar("mip_obj_value")
+        h5.file.close()
+    end
+    
+    @timeit "Initialize" begin
+        count_backtrack = 0
+        count_deterioration = 0
+        gapcl_best = 0
+        gapcl_best_history = CircularBuffer{Float64}(gapcl_best_patience)
+        gapcl_curr = 0
+        last_print_time = 0
+        multipliers_best = Float64[]
+        multipliers_curr = Float64[]
+        obj_best = nothing
+        obj_curr = nothing
+        obj_hist = CircularBuffer{Float64}(100)
+        obj_initial = nothing
+        pool = nothing
+        pool_cut_age = nothing
+        pool_cut_hashes = Set{UInt64}()
+        pool_size_mb = 0
+        tableau_density::Float32 = 0.05
+        basis_cache = nothing
+        λ, Δ = 0, 0
+        μ = 10
+        
+        basis_vars_to_id = Dict()
+        basis_id_to_vars = Dict{Int, Vector{Int}}()
+        basis_id_to_sizes = Dict{Int, Vector{Int}}()
+        next_basis_id = 1
+        cut_basis_id = Int[]
+        cut_row = Int[]
+    end
+
+    gapcl(v) = 100 * (v - obj_initial) / (obj_mip - obj_initial)
+
+    @timeit "Read problem" begin
+        model = read_from_file(mps_filename)
+        set_optimizer(model, optimizer)
+    end
+
+    @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)
+        for (var_idx, var_type) in enumerate(data.var_types)
+            if var_type in ['B', 'I']
+                assert_int(sol_opt[var_idx])
+            end
+        end
+
+        # Convert to standard form
+        data_s, transforms = convert_to_standard_form(data)
+        model_s = to_model(data_s)
+        if silent_solver
+            set_silent(model_s)
+        end
+        vars_s = all_variables(model_s)
+        orig_obj_s = objective_function(model_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
+        for (var_idx, var_type) in enumerate(data_s.var_types)
+            if var_type in ['B', 'I']
+                assert_int(sol_opt_s[var_idx])
+            end
+        end
+        assert_eq(data_s.constr_lhs * sol_opt_s, data_s.constr_lb)
+    end
+
+    @info "Standard form model has $(length(data_s.var_lb)) vars, $(length(data_s.constr_lb)) constrs"
+
+    for round = 1:max_rounds
+        log_prefix = ' '
+        log_should_print = false
+        is_last_iteration = false
+        if round == max_rounds
+            is_last_iteration = true
+        end
+
+        elapsed_time = time() - initial_time
+        if elapsed_time > time_limit
+            @info "Time limit exceeded. Stopping after current iteration."
+            is_last_iteration = true
+        end
+
+        if round > 1
+            @timeit "Build Lagrangian term" begin
+                @timeit "mul" begin
+                    active_idx = findall(multipliers_curr .> 1e-6)
+                    v = sparse(pool.lhs[:, active_idx] * multipliers_curr[active_idx])
+                end
+                @timeit "dot" begin
+                    lagr_term = AffExpr(dot(multipliers_curr, pool.lb))
+                end
+                @timeit "add_to_expression!" begin
+                    for offset in 1:nnz(v)
+                        var_idx = v.nzind[offset]
+                        add_to_expression!(
+                            lagr_term,
+                            vars_s[var_idx],
+                            - v.nzval[offset],
+                        )
+                    end
+                end
+            end
+            @timeit "Update objective" begin
+                set_objective_function(
+                    model_s,
+                    orig_obj_s + lagr_term,
+                )
+            end
+        end
+
+        @timeit "Optimize LP (lagrangian)" begin
+            basis_cache === nothing || set_basis(model_s, basis_cache)
+            set_silent(model_s)
+            optimize!(model_s)
+            basis_cache = get_basis(model_s)
+            status = termination_status(model_s)
+            if status == MOI.DUAL_INFEASIBLE
+                @warn "LP is unbounded (dual infeasible). Resetting to best known multipliers."
+                copy!(multipliers_curr, multipliers_best)
+                obj_curr = obj_best
+                continue
+            elseif status != MOI.OPTIMAL
+                error("Non-optimal termination status: $status")
+            end            
+            sol_frac = get_x(model_s)
+            obj_curr = objective_value(model_s)
+        end
+
+        @timeit "Update history and μ" begin
+            push!(obj_hist, obj_curr)
+            if obj_best === nothing || obj_curr > obj_best
+                log_prefix = '*'
+                obj_best = obj_curr
+                copy!(multipliers_best, multipliers_curr)
+            end
+            if round == 1
+                obj_initial = obj_curr
+            end
+            gapcl_curr = gapcl(obj_curr)
+            gapcl_best = gapcl(obj_best)
+            push!(gapcl_best_history, gapcl_best)
+            if variant in [:subg, :hybr]
+                Δ = obj_mip - obj_best
+                if obj_curr < obj_best - Δ
+                    count_deterioration += 1
+                else
+                    count_deterioration = 0
+                end
+                if count_deterioration >= 10
+                    μ *= 0.5
+                    copy!(multipliers_curr, multipliers_best)
+                    count_deterioration = 0
+                    count_backtrack += 1
+                elseif length(obj_hist) >= 100
+                    obj_hist_avg = mean(obj_hist)
+                    improv = obj_best - obj_hist[1]
+                    if improv < 0.01 * Δ
+                        if obj_best - obj_hist_avg < 0.001 * Δ
+                            μ = 10 * μ
+                        elseif obj_best - obj_hist_avg < 0.01 * Δ
+                            μ = 2 * μ
+                        else
+                            μ = 0.5 * μ
+                        end
+                    end
+                end
+            elseif variant in [:fast, :faster, :miplearn]
+                μ = 0.01
+            else
+                error("not implemented")
+            end
+        end
+
+        if mod(round - 1, interval_read_tableau) == 0
+            @timeit "Get basis" begin
+                basis = get_basis(model_s)
+            end
+            @timeit "Select tableau rows" begin
+                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, estimated_density=tableau_density * 1.05)
+                tableau_density = nnz(tableau.lhs) / length(tableau.lhs)
+                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)
+            end
+
+            @timeit "Check cut validity" begin
+                assert_cuts_off(cuts_s, sol_frac)
+                assert_does_not_cut_off(cuts_s, sol_opt_s)
+                ncuts = length(cuts_s.lb)
+            end
+
+            @timeit "Add new cuts to the pool" begin
+                @timeit "Compute cut hashses" begin
+                    unique_indices = Int[]
+                    for i in 1:ncuts
+                        cut_hash = cuts_s.hash[i]
+                        if !(cut_hash in pool_cut_hashes)
+                            push!(pool_cut_hashes, cut_hash)
+                            push!(unique_indices, i)
+                        end
+                    end
+                end
+                @timeit "Append unique cuts" begin
+                    @timeit "Track basis" begin
+                        vb = basis.var_basic
+                        vn = basis.var_nonbasic
+                        cb = basis.constr_basic
+                        cn = basis.constr_nonbasic
+                        basis_vars = [vb; vn; cb; cn]
+                        basis_sizes = [length(vb), length(vn), length(cb), length(cn)]
+                        
+                        if basis_vars ∉ keys(basis_vars_to_id)
+                            basis_id = next_basis_id
+                            basis_vars_to_id[basis_vars] = basis_id
+                            basis_id_to_vars[basis_id] = basis_vars
+                            basis_id_to_sizes[basis_id] = basis_sizes
+                            next_basis_id += 1
+                        else
+                            basis_id = basis_vars_to_id[basis_vars]
+                        end
+                    end
+                    
+                    if round == 1
+                        pool = ConstraintSet(
+                            lhs = sparse(cuts_s.lhs[unique_indices, :]'),
+                            lb = cuts_s.lb[unique_indices],
+                            ub = cuts_s.ub[unique_indices],
+                            hash = cuts_s.hash[unique_indices],
+                        )
+                        ncuts_unique = length(unique_indices)
+                        multipliers_curr = zeros(ncuts_unique)
+                        multipliers_best = zeros(ncuts_unique)
+                        pool_cut_age = zeros(ncuts_unique)
+                        for i in unique_indices
+                            push!(cut_basis_id, basis_id)
+                            push!(cut_row, selected_rows[i])
+                        end
+                    else
+                        if !isempty(unique_indices)
+                            @timeit "Append LHS" begin
+                                # Transpose cuts matrix for better performance
+                                new_cuts_lhs = sparse(cuts_s.lhs[unique_indices, :]')
+                                
+                                # Resize existing matrix in-place to accommodate new columns
+                                old_cols = pool.lhs.n
+                                new_cols = new_cuts_lhs.n
+                                total_cols = old_cols + new_cols
+                                resize!(pool.lhs.colptr, total_cols + 1)
+                                
+                                # Append new column pointers with offset
+                                old_nnz = nnz(pool.lhs)
+                                for i in 1:new_cols
+                                    pool.lhs.colptr[old_cols + i + 1] = old_nnz + new_cuts_lhs.colptr[i + 1]
+                                end
+                                
+                                # Expand rowval and nzval arrays
+                                append!(pool.lhs.rowval, new_cuts_lhs.rowval)
+                                append!(pool.lhs.nzval, new_cuts_lhs.nzval)
+                                
+                                # Update matrix dimensions
+                                pool.lhs = SparseMatrixCSC(pool.lhs.m, total_cols, pool.lhs.colptr, pool.lhs.rowval, pool.lhs.nzval)
+                            end
+                            @timeit "Append others" begin
+                                ncuts_unique = length(unique_indices)
+                                append!(pool.lb, cuts_s.lb[unique_indices])
+                                append!(pool.ub, cuts_s.ub[unique_indices])
+                                append!(pool.hash, cuts_s.hash[unique_indices])
+                                append!(multipliers_curr, zeros(ncuts_unique))
+                                append!(multipliers_best, zeros(ncuts_unique))
+                                append!(pool_cut_age, zeros(ncuts_unique))
+                                for i in unique_indices
+                                    push!(cut_basis_id, basis_id)
+                                    push!(cut_row, selected_rows[i])
+                                end
+                            end
+                        end
+                    end
+                end
+            end
+
+
+            @timeit "Prune the pool" begin
+                pool_size_mb = Base.summarysize(pool) / 1024^2
+                while pool_size_mb >= max_pool_size_mb
+                    @timeit "Identify cuts to remove" begin
+                        scores = collect(zip(multipliers_best .> 1e-6, -pool_cut_age))
+                        σ = sortperm(scores, rev=true)
+                        pool_size = length(pool.ub)
+                        n_keep = Int(floor(pool_size * 0.8))
+                        idx_keep = σ[1:n_keep]
+                        idx_remove = σ[(n_keep+1):end]
+                        
+                        positive_multipliers_dropped = sum(multipliers_best[idx_remove] .> 1e-6)
+                        @info "Dropping $(length(idx_remove)) cuts ($(positive_multipliers_dropped) with multipliers_best)"
+                    end
+                    @timeit "Update cut hashes" begin
+                        for idx in idx_remove
+                            cut_hash = pool.hash[idx]
+                            delete!(pool_cut_hashes, cut_hash)
+                        end
+                    end
+                    @timeit "Update cut pool" begin
+                        pool.lhs = pool.lhs[:, idx_keep]
+                        pool.ub = pool.ub[idx_keep]
+                        pool.lb = pool.lb[idx_keep]
+                        pool.hash = pool.hash[idx_keep]
+                        multipliers_curr = multipliers_curr[idx_keep]
+                        multipliers_best = multipliers_best[idx_keep]
+                        pool_cut_age = pool_cut_age[idx_keep]
+                        cut_basis_id = cut_basis_id[idx_keep]
+                        cut_row = cut_row[idx_keep]
+                    end
+                    @timeit "Update known bases" begin
+                        used_basis_ids = Set(cut_basis_id)
+                        for basis_id in collect(keys(basis_id_to_vars))
+                            if basis_id ∉ used_basis_ids
+                                basis_vars = basis_id_to_vars[basis_id]
+                                delete!(basis_vars_to_id, basis_vars)
+                                delete!(basis_id_to_vars, basis_id)
+                                delete!(basis_id_to_sizes, basis_id)
+                            end
+                        end
+                    end
+                    pool_size_mb = Base.summarysize(pool) / 1024^2
+                end
+            end
+        end
+        
+        if mod(round - 1, interval_large_lp) == 0 || is_last_iteration
+            log_should_print = true
+            @timeit "Update multipliers (large LP)" begin
+                selected_idx = []
+                selected_contrs = []
+                while true
+                    @timeit "Optimize LP (extended)" begin
+                        set_silent(model_s)
+                        set_objective_function(model_s, orig_obj_s)
+                        optimize!(model_s)
+                        status = termination_status(model_s)
+                        if status != MOI.OPTIMAL
+                            error("Non-optimal termination status: $status")
+                        end
+                        obj_curr = objective_value(model_s)
+                        sol_frac = get_x(model_s)
+                    end
+
+                    @timeit "Computing cut violations" begin
+                        violations = pool.lb - pool.lhs' * sol_frac
+                    end
+
+                    @timeit "Sorting cut violations" begin
+                        σ = sortperm(violations, rev=true)
+                    end
+
+                    if violations[σ[1]] <= 1e-6
+                        break
+                    end
+
+                    @timeit "Add constraints to the model" begin
+                        ncuts = min(max(1, sum(violations .> 1e-6) ÷ 10), length(σ))
+                        for i in 1:ncuts
+                            if violations[σ[i]] <= 1e-6
+                                break
+                            end
+                            cut_lhs = pool.lhs[:, σ[i]]
+                            cut_lhs_value = 0.0
+                            cut_lb = pool.lb[σ[i]]
+                            cut_expr = AffExpr()
+                            for offset in 1:nnz(cut_lhs)
+                                var_idx = cut_lhs.nzind[offset]
+                                add_to_expression!(
+                                    cut_expr,
+                                    vars_s[var_idx],
+                                    cut_lhs.nzval[offset],
+                                )
+                                cut_lhs_value += sol_frac[var_idx] * cut_lhs.nzval[offset]
+                            end
+                            cut_constr = @constraint(model_s, cut_expr >= cut_lb)
+                            push!(selected_idx, σ[i])
+                            push!(selected_contrs, cut_constr)
+                        end
+                    end
+                end
+
+                @timeit "Find dual values for all selected cuts" begin
+                    multipliers_curr .= 0
+                    pool_cut_age .+= 1
+                    for (offset, idx) in enumerate(selected_idx)
+                        multipliers_curr[idx] = -shadow_price(selected_contrs[offset])
+                        if multipliers_curr[idx] > 1e-5
+                            pool_cut_age[idx] = 0
+                        end
+                    end
+                end
+
+                @timeit "Update best" begin
+                    if obj_curr > obj_best
+                        log_prefix = '*'
+                        obj_best = obj_curr
+                        copy!(multipliers_best, multipliers_curr)
+                    end
+                    gapcl_curr = gapcl(obj_curr)
+                    gapcl_best = gapcl(obj_best)
+                end
+
+                @timeit "Delete all cut constraints" begin
+                    delete.(model_s, selected_contrs)
+                end
+            end
+        else
+            @timeit "Update multipliers (subgradient)" begin
+                subgrad = (pool.lb' - (sol_frac' * pool.lhs))'
+                subgrad_norm_sq = norm(subgrad)^2
+                if subgrad_norm_sq < 1e-10
+                    λ = 0
+                else
+                    λ = μ * (obj_mip - obj_curr) / subgrad_norm_sq
+                end
+                multipliers_curr = max.(0, multipliers_curr .+ λ * subgrad)
+            end
+        end
+
+        if round == 1
+            @printf(
+                "  %8s %8s %10s %9s %9s %9s %9s %4s %8s %8s %8s\n",
+                "time",
+                "round",
+                "obj",
+                "cl_curr",
+                "cl_best",
+                "pool_cuts",
+                "pool_mb",
+                "bktk",
+                "Δ",
+                "μ",
+                "λ",
+            )
+        end
+
+        if time() - last_print_time > interval_print_sec
+            log_should_print = true
+        end
+
+        if is_last_iteration
+            log_should_print = true
+        end
+
+        if log_should_print
+            last_print_time = time()
+            @printf(
+                "%c %8.2f %8d %10.3e %9.2e %9.2e %9d %9.2f %4d %8.2e %8.2e %8.2e\n",
+                log_prefix,
+                time() - initial_time,
+                round,
+                obj_curr,
+                gapcl_curr,
+                gapcl_best,
+                length(pool.ub),
+                pool_size_mb,
+                count_backtrack,
+                Δ,
+                μ,
+                λ,
+            )
+        end
+        
+        push!(gapcl_best_history, gapcl_best)
+        if length(gapcl_best_history) >= gapcl_best_patience
+            if gapcl_best <= gapcl_best_history[1]
+                @info "No gap closure improvement. Stopping."
+                break
+            end
+        end
+
+        if is_last_iteration
+            break
+        end
+    end
+
+    @info "Best gap closure: $(gapcl_best)"
+    
+    @timeit "Keep only active cuts" begin
+        positive_idx = findall(multipliers_best .> 1e-6)
+        if length(positive_idx) == 0 && gapcl_best > 1e-6
+            error("gap closure with zero cuts")
+        end
+
+        @timeit "Clean up cut pool" begin
+            pool.lhs = pool.lhs[:, positive_idx]
+            pool.lb = pool.lb[positive_idx]
+            pool.ub = pool.ub[positive_idx]
+            pool.hash = pool.hash[positive_idx]
+            multipliers_best = multipliers_best[positive_idx]
+            multipliers_curr = multipliers_curr[positive_idx]
+            cut_basis_id = cut_basis_id[positive_idx]
+            cut_row = cut_row[positive_idx]
+        end
+
+        @timeit "Clean up known bases" begin
+            used_basis_ids = Set(cut_basis_id)
+            for basis_id in collect(keys(basis_id_to_vars))
+                if basis_id ∉ used_basis_ids
+                    basis_vars = basis_id_to_vars[basis_id]
+                    delete!(basis_vars_to_id, basis_vars)
+                    delete!(basis_id_to_vars, basis_id)
+                    delete!(basis_id_to_sizes, basis_id)
+                end
+            end
+        end
+
+        @info "Keeping $(length(positive_idx)) cuts from $(length(used_basis_ids)) unique bases"
+    end
+
+    to = TimerOutputs.get_defaulttimer()
+    stats_time = TimerOutputs.tottime(to) / 1e9
+    print_timer()
+
+    if length(positive_idx) > 0
+        @timeit "Write cuts to H5" begin
+            if !isempty(cut_basis_id)
+                @timeit "Convert IDs to offsets" begin
+                    id_to_offset = Dict{Int, Int}()
+                    gmi_basis_vars = []
+                    gmi_basis_sizes = []
+                    for (offset, basis_id) in enumerate(sort(collect(keys(basis_id_to_vars))))
+                        id_to_offset[basis_id] = offset
+                        push!(gmi_basis_vars, basis_id_to_vars[basis_id])
+                        push!(gmi_basis_sizes, basis_id_to_sizes[basis_id])
+                    end
+                    gmi_cut_basis = [id_to_offset[basis_id] for basis_id in cut_basis_id]
+                    gmi_cut_row = cut_row
+                end
+                
+                @timeit "Convert to matrices" begin
+                    gmi_basis_vars_matrix = hcat(gmi_basis_vars...)'
+                    gmi_basis_sizes_matrix = hcat(gmi_basis_sizes...)'
+                end
+
+                @timeit "Write H5" begin
+                    h5 = H5File(h5_filename, "r+")
+                    h5.put_array("gmi_basis_vars", gmi_basis_vars_matrix)
+                    h5.put_array("gmi_basis_sizes", gmi_basis_sizes_matrix)
+                    h5.put_array("gmi_cut_basis", gmi_cut_basis)
+                    h5.put_array("gmi_cut_row", gmi_cut_row)
+                    h5.file.close()
+                end
+            end
+        end
+
+        if verify_cuts
+            @timeit "Verify cuts in current model" begin
+                @info "Verifying cuts in current standard form model using pool..."
+                if !isempty(cut_basis_id)
+                    @info "Adding $(length(pool.lb)) cuts from pool to current model"
+                    pool.lhs = sparse(pool.lhs')
+                    constrs = build_constraints(model_s, pool)
+                    add_constraint.(model_s, constrs)
+                    set_objective_function(model_s, orig_obj_s)
+                    optimize!(model_s)
+                    status = termination_status(model_s)
+                    if status != MOI.OPTIMAL
+                        error("Non-optimal termination status: $status")
+                    end
+                    obj_verify_s = objective_value(model_s)
+                    gapcl_verify_s = gapcl(obj_verify_s)
+                    @show gapcl_verify_s
+                    @show gapcl_best
+                    if abs(gapcl_best - gapcl_verify_s) > 0.01
+                        error("Gap closures differ: $(gapcl_best) ≠ $(gapcl_verify_s)")
+                    end
+                    @info "Current model gap closure matches: $(gapcl_best) ≈ $(gapcl_verify_s)"
+                else
+                    @warn "No cuts in pool to verify"
+                end
+            end
+
+            @timeit "Verify stored cuts" begin
+                @info "Verifying stored cuts..."
+                model_verify = read_from_file(mps_filename)
+                set_optimizer(model_verify, optimizer)
+                verification_cuts = _dualgmi_generate([h5_filename], model_verify; test_h5=h5_filename)
+                constrs = build_constraints(model_verify, verification_cuts)
+                add_constraint.(model_verify, constrs)
+                relax_integrality(model_verify)
+                optimize!(model_verify)
+                status = termination_status(model_verify)
+                if status != MOI.OPTIMAL
+                    error("Non-optimal termination status: $status")
+                end
+                obj_verify = objective_value(model_verify)
+                gapcl_verify = gapcl(obj_verify)
+                @show gapcl_verify
+                @show gapcl_best
+                if abs(gapcl_best - gapcl_verify) > 0.01
+                    error("Gap closures differ: $(gapcl_best) ≠ $(gapcl_verify)")
+                end
+                @info "Gap closure matches gapcl_best: $(gapcl_best) ≈ $(gapcl_verify)"
+            end
+        end
+    end
+
+    return OrderedDict(
+        "gapcl_best" => gapcl_best,
+        "gapcl_curr" => gapcl_curr,
+        "instance" => mps_filename,
+        "obj_final" => obj_curr,
+        "obj_initial" => obj_initial,
+        "obj_mip" => obj_mip,
+        "pool_size_mb" => pool_size_mb,
+        "pool_total" => length(pool.lb),
+        "time" => stats_time,
+    )
+end
+
 function add_constraint_set_dual_v2(model::JuMP.Model, cs::ConstraintSet)
     vars = all_variables(model)
     nrows, ncols = size(cs.lhs)
@@ -371,13 +1055,34 @@ function _dualgmi_compress_h5(h5_filename)
     h5.file.close()
 end
 
-function _dualgmi_generate(train_h5, model)
+function _dualgmi_generate(train_h5, model; test_h5=nothing)
     @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 "Read optimal solution from test H5" begin
+        sol_opt_dict = nothing
+        sol_opt = nothing
+        sol_opt_s = nothing
+        if test_h5 !== nothing
+            try
+                h5 = H5File(test_h5, "r")
+                var_names = h5.get_array("static_var_names")
+                var_values = h5.get_array("mip_var_values")
+                h5.close()
+                if var_names !== nothing && var_values !== nothing
+                    sol_opt_dict = Dict(zip(var_names, convert(Array{Float64}, var_values)))
+                    sol_opt = [sol_opt_dict[n] for n in data.var_names]
+                    sol_opt_s = forward(transforms, sol_opt)
+                    @info "Loaded optimal solution for cut validation"
+                end
+            catch e
+                @warn "Could not read optimal solution from test H5 file: $e"
+            end
+        end
+    end
     @timeit "Collect cuts from H5 files" begin
         basis_vars_to_basis_offset = Dict()
         combined_basis_sizes = nothing
@@ -446,6 +1151,10 @@ function _dualgmi_generate(train_h5, model)
                 tableau = compute_tableau(data_s, current_basis; rows=collect(rows))
                 cuts_s = compute_gmi(data_s, tableau)
                 cuts = backwards(transforms, cuts_s)
+                if sol_opt_s !== nothing && sol_opt !== nothing
+                    assert_does_not_cut_off(cuts_s, sol_opt_s)
+                    assert_does_not_cut_off(cuts, sol_opt)
+                end
                 if all_cuts === nothing
                     all_cuts = cuts
                 else
@@ -576,5 +1285,5 @@ function __init_gmi_dual__()
     copy!(ExpertDualGmiComponent, ExpertDualGmiComponentPy)
 end
 
-export collect_gmi_dual, expert_gmi_dual, ExpertDualGmiComponent, KnnDualGmiComponent
+export collect_gmi_dual, expert_gmi_dual, ExpertDualGmiComponent, KnnDualGmiComponent, collect_gmi_FisSal2011
 

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 (vals[i] <= cuts.ub[i] - tol) && (vals[i] >= cuts.lb[i] + tol)
-        if (val <= cuts.ub[i] - tol) && (val >= cuts.lb[i] + tol)
+            throw(ErrorException("inequality $i fails to cut off fractional solution: $(cuts.lb[i]) <= $(vals[i]) <= $(cuts.ub[i])"))
-            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)
+    vals = cuts.lhs * x
     for i = 1:length(cuts.lb)
-        val = cuts.lhs[i, :]' * x
+        if (vals[i] >= cuts.ub[i]) || (vals[i] <= cuts.lb[i])
-        ub = cuts.ub[i]
+            throw(ErrorException("inequality $i cuts off integer solution: $(cuts.lb[i]) <= $(vals[i]) <= $(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
+
+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
+    obj::Vector{Float64}
-    lhs::Any
+    lhs::SparseMatrixCSC
-    rhs::Any
+    rhs::Vector{Float64}
-    z::Any
+    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[]
+    if isa(unsafe_backend(model), Gurobi.Optimizer)
-    var_nonbasic = Int[]
+        return get_basis_gurobi(model)
-    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
     end
 
-    # Constraints
+    @timeit "Initialization" begin
-    constr_index = 1
+        var_basic = Int[]
-    for (ftype, stype) in list_of_constraint_types(model)
+        var_nonbasic = Int[]
-        for constr in all_constraints(model, ftype, stype)
+        constr_basic = Int[]
-            if ftype == VariableRef
+        constr_nonbasic = Int[]
-                # nop
+        nvars = num_variables(model)
-            elseif ftype == AffExpr
+        sizehint!(var_basic, nvars)
-                bstatus = MOI.get(model, MOI.ConstraintBasisStatus(), constr)
+        sizehint!(var_nonbasic, nvars)
-                if bstatus == MOI.BASIC
+    end
-                    push!(constr_basic, constr_index)
+
-                elseif bstatus == MOI.NONBASIC
+    @timeit "Query variables" begin
-                    push!(constr_nonbasic, constr_index)
+        for (i, var) in enumerate(all_variables(model))
-                else
+            bstatus = MOI.get(model, MOI.VariableBasisStatus(), var)
-                    error("Unknown basis status: $bstatus")
+            if bstatus == MOI.BASIC
-                end
+                push!(var_basic, i)
-                constr_index += 1
+            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 arrays" begin
-        @timeit "Initialize" begin
+        num_rows = length(rows)
-            tableau_rhs = zeros(length(rows))
+        tableau_rhs::Array{Float64} = zeros(num_rows)
-            tableau_lhs = zeros(length(rows), ncols)
+        tableau_rowptr::Array{Int} = zeros(Int, num_rows + 1)
-        end
+        tableau_colval::Array{Int} = Int[]
-        for k in eachindex(1:length(rows))
+        tableau_nzval::Array{Float64} = Float64[]
-            @timeit "Prepare inputs" begin
+        estimated_nnz::Int = round(num_rows * ncols * estimated_density)
-                i = rows[k]
+        sizehint!(tableau_colval, estimated_nnz)
-                e = zeros(nrows)
+        sizehint!(tableau_nzval, estimated_nnz)
-                e[i] = 1.0
+        e::Array{Float64} = zeros(nrows)
-            end
+        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
+            @timeit "Compute row" begin
-                tableau_lhs[k, :] = sol' * data.constr_lhs
+                mul!(tableau_row, A, sol)
-                tableau_rhs[k] = sol' * data.constr_ub
+                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
+    end
-            tableau_lhs[abs.(tableau_lhs) .<= tol] .= 0
+
-            tableau_lhs = sparse(tableau_lhs)
+    @timeit "Shrink arrays" begin
-        end
+        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)
+    @timeit "Identify constraint type" begin
-    isequality = abs.(data.constr_ub .- data.constr_lb) .< 1e-6
+        nrows, ncols = size(data.constr_lhs)
-    eq = [i for i = 1:nrows if isequality[i]]
+        isequality = abs.(data.constr_ub .- data.constr_lb) .< 1e-6
-    ge = [i for i = 1:nrows if isfinite(data.constr_lb[i]) && !isequality[i]]
+        eq = [i for i = 1:nrows if isequality[i]]
-    le = [i for i = 1:nrows if isfinite(data.constr_ub[i]) && !isequality[i]]
+        ge = [i for i = 1:nrows if isfinite(data.constr_lb[i]) && !isequality[i]]
-    EQ, GE, LE = length(eq), length(ge), length(le)
+        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 = [
+    end
-        I spzeros(ncols, GE + LE)
+    @timeit "Identify slack type" begin
-        data.constr_lhs[ge, :] spzeros(GE, GE + LE)
+        constr_lhs_t = sparse(data.constr_lhs')
-        -data.constr_lhs[le, :] spzeros(LE, GE + LE)
+        function is_integral(row_idx, rhs)
-    ]
+            rhs_is_integer = abs(rhs - round(rhs)) <= 1e-6
-    t.M2 = [
+            cols, coeffs = findnz(constr_lhs_t[:, row_idx])[1:2]
-        zeros(ncols)
+            vars_are_integer = all(j -> data.var_types[j] ∈ ['I', 'B'], cols)
-        data.constr_lb[ge]
+            coeffs_are_integer = all(v -> abs(v - round(v)) <= 1e-6, coeffs)
-        -data.constr_ub[le]
+            return rhs_is_integer && vars_are_integer && coeffs_are_integer
-    ]
+        end
-    t.ncols_orig = ncols
+        slack_types = [
-    t.GE, t.LE = GE, LE
+            [is_integral(ge[i], data.constr_lb[ge[i]]) ? 'I' : 'C' for i = 1:GE];
-    t.lhs_ge = data.constr_lhs[ge, :]
+            [is_integral(le[i], data.constr_ub[le[i]]) ? 'I' : 'C' for i = 1:LE]
-    t.lhs_le = data.constr_lhs[le, :]
+        ]
-    t.rhs_ge = data.constr_lb[ge]
+    end
-    t.rhs_le = data.constr_ub[le]
+    @timeit "Build M1" begin
-
+        t.M1 = [
-    data.constr_lhs = [
+            I spzeros(ncols, GE + LE)
-        data.constr_lhs[eq, :] spzeros(EQ, GE) spzeros(EQ, LE)
+            data.constr_lhs[ge, :] spzeros(GE, GE + LE)
-        data.constr_lhs[ge, :] -I spzeros(GE, LE)
+            -data.constr_lhs[le, :] spzeros(LE, GE + LE)
-        data.constr_lhs[le, :] spzeros(LE, GE) I
+        ]
-    ]
+    end
-    data.obj = [data.obj; zeros(GE + LE)]
+    @timeit "Build M2" begin
-    data.var_lb = [data.var_lb; zeros(GE + LE)]
+        t.M2 = [
-    data.var_ub = [data.var_ub; [Inf for _ = 1:(GE+LE)]]
+            zeros(ncols)
-    data.var_names = [data.var_names; ["__s$i" for i = 1:(GE+LE)]]
+            data.constr_lb[ge]
-    data.var_types = [data.var_types; ['C' for _ = 1:(GE+LE)]]
+            -data.constr_ub[le]
-    data.constr_lb = [
+        ]
-        data.constr_lb[eq]
+    end
-        data.constr_lb[ge]
+    @timeit "Build t.lhs, t.rhs" begin
-        data.constr_ub[le]
+        t.ncols_orig = ncols
-    ]
+        t.GE, t.LE = GE, LE
-    data.constr_ub = copy(data.constr_lb)
+        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
+    ncols::Int = 0
-    B::Int = 0
+    is_var_free::Vector{Bool} = []
-    free::Vector{Int} = []
-    others::Vector{Int} = []
 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])]
+    is_var_free = [!isfinite(data.var_lb[i]) && !isfinite(data.var_ub[i]) for i = 1:ncols]
-    others = [i for i = 1:ncols if isfinite(data.var_lb[i]) || isfinite(data.var_ub[i])]
+    free_idx = findall(is_var_free)
-    t.F = length(free)
-    t.B = length(others)
-    t.free, t.others = free, others
     data.obj = [
-        data.obj[others]
+        data.obj
-        data.obj[free]
+        [-data.obj[i] for i in free_idx]
-        -data.obj[free]
     ]
-    data.constr_lhs = [lhs[:, others] lhs[:, free] -lhs[:, free]]
     data.var_lb = [
-        data.var_lb[others]
+        [is_var_free[i] ? 0.0 : data.var_lb[i] for i in 1:ncols]
-        [0.0 for _ in free]
+        [0 for _ in free_idx]
-        [0.0 for _ in free]
     ]
     data.var_ub = [
-        data.var_ub[others]
+        [is_var_free[i] ? Inf : data.var_ub[i] for i in 1:ncols]
-        [Inf for _ in free]
+        [Inf for _ in free_idx]
-        [Inf for _ in free]
     ]
     data.var_types = [
-        data.var_types[others]
+        data.var_types
-        data.var_types[free]
+        [data.var_types[i] for i in free_idx]
-        data.var_types[free]
     ]
     data.var_names = [
-        data.var_names[others]
+        data.var_names
-        ["$(v)_p" for v in data.var_names[free]]
+        ["$(data.var_names[i])_neg" for i in free_idx]
-        ["$(v)_m" for v in data.var_names[free]]
     ]
+    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
+    ncols, is_var_free = t.ncols, t.is_var_free
-    nrows, _ = size(c.lhs)
+    free_idx = findall(is_var_free)
-    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
 
-    # Assert only LE constraints are left (EQ constraints are not supported)
+    for (offset, var_idx) in enumerate(free_idx)
-    @assert all(c.lb .== -Inf)
+        @assert c.lhs[:, var_idx] == -c.lhs[:, ncols+offset]
-
-    # 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])
     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]
+        [is_var_free[i] ? max(0, p[i]) : p[i] for i in 1:ncols]
-        max.(p[t.free], 0)
+        [max(0, -p[i]) for i in free_idx]
-        max.(-p[t.free], 0)
     ]
 end