gurobi_solution_split

gurobi_solution_split.jl

## Packages
using JuMP
using Gurobi

using BenchmarkTools

# function example_big_sum(rowsum, colsum; verbose::Bool = true)
function create_model(rowsum, colsum; verbose::Bool = true)
    ## Model:
    model = Model(Gurobi.Optimizer)
    @assert length(rowsum) == length(colsum)

    N = length(rowsum)
    INDEX = 1:N
    @variable(model, y[INDEX,INDEX] >= 0, Int)

    JuMP.@constraint(model, rowCons[i=INDEX], sum(y[i,j] for j in INDEX) == rowsum[i])
    JuMP.@constraint(model, colCons[j=INDEX], sum(y[i,j] for i in INDEX) == colsum[j])

    if verbose; print(model); end

    ## Gurobi parameters - see https://www.gurobi.com/documentation/9.0/refman/finding_multiple_solutions.html
    JuMP.set_optimizer_attribute(model, "PoolSearchMode", 2)   # exhaustive search mode
    JuMP.set_optimizer_attribute(model, "PoolSolutions", 100)  # 100 is an arbitrary (large enough) whole number

    if !verbose; JuMP.set_optimizer_attribute(model, "OutputFlag", 0); end

    return model
end

function optimize_model(model::JuMP.AbstractModel)
    ## Optimize:
    JuMP.optimize!(model)
    return nothing
end

function print_results(model::JuMP.AbstractModel,N::Int64; verbose::Bool = true)
    INDEX = 1:N
    y = model[:y]
    ## Results:
    num_results = result_count(model)
    if verbose; println("Number of results: ", num_results, "\n"); end

    Results = Array{Array{Int64},1}()  ## Note the conversion to Int64
    for n in 1:num_results
        sol = Array{Int64}(undef, N, N)
        for i in INDEX
            for j in INDEX
                sol[i,j] = JuMP.value(y[i,j]; result=n)
            end
        end

        push!(Results,sol)
        if verbose
            println(Results[n])
        end
    end

    return Results
end

μ = [1,1,0]
ν = [1,0,1]

model = create_model(μ, ν; verbose=false)
@time optimize_model(model)
@btime optimize_model(model) samples=10
print_results(model,length(μ),verbose=false)

# @btime example_big_sum($μ, $ν; verbose=false) samples=10