Hi @jd-foster .
Thanks for your help. Right now I am not working on this, so I will edit this comment when I am able to test your response.
No problem.
While it’s fresh in my mind, here is a more general solution for future reference:
## Packages
using JuMP
using Gurobi
function example_big_sum(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)
@constraint(model, rowCons[i=INDEX], sum(y[i,j] for j in INDEX) == rowsum[i])
@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
## Optimize:
JuMP.optimize!(model)
## 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
μ = [0,2,3,1]
ν = [2,2,1,1]
solutions = example_big_sum(μ,ν);
FYI: 24 solutions found as other posts show.
7 Likes
@jd-foster I was able to compare your solution with that provided by @dpsanders. Your solution is faster (no surprise here, I suppose) and produce the same results. Below is the function (and its helper function) that I used to compare against your solution.
n::Int64 = length(μ)
vars = (@variables y[1:n, 1:n])[1]
contractors::Array{Contractor,1} = Contractor[]
for i in 1:n
push!(contractors, Contractor(vars, sum(vars[i, 1:n]) - μ[i]))
end
for i in 1:n
push!(contractors, Contractor(vars, sum(vars[1:n, i]) - ν[i]))
end
X::IntervalBox = IntervalBox(0..n^2, n^2)
X, contractors
end
function find_solutions(μ::Array{Int64,1}, ν::Array{Int64,1})
n::Int64 = length(μ)
x::IntervalBox{n^2, Float64}, c::Array{Contractor,1} = generate_contractors(μ,ν)
contractors::Array{Function,1} = [X -> C(0..0, X) for C in c]
helper = pop!(contractors)
X::IntervalBox{n^2, Float64} = Base.invokelatest(helper, x)
for C in contractors
X = Base.invokelatest(C, X)
end
contractors = [contractors; integerize]
solutions::Array{IntervalBox{n^2,Float64},1} = Base.invokelatest(pave, contractors,[X], 1.0)
[Int.(x) for x in mid.(solutions)]
end
BTW, the method using IntervalArithmetic.jl no longer works due to an incompatibility between this and ModelingToolkit.jl. You need to downgrade it. This is how I did it using Pkg REPL
pkg > remove ModelingToolkit
pkg > add ModelingToolkit@3.16.0
pkg > pin ModelingToolkit@3.16.0
1 Like
What are the respective timings and what is the benchmarking code? Can you post a complete code that can be copied and pasted?
I am far from being a expert benchmarking, but here is what I did using @btime from BenchmarkingTools:
μ = [1,1,0]
ν = [1,0,1]
@btime find_solutions($μ,$ν) samples=10
and I obtained 481.762 ms (856141 allocations: 52.75 MiB). Meanwhile, for the Gurobi solution:
@btime example_big_sum($μ, $ν; verbose=false) samples=10
I obtain 5.243 ms (3111 allocations: 188.68 KiB).
I will post the complete code later.
Here’s my code for the solution using interval arithmetic:
1 Like
You are benchmarking everything, including the construction of the model. While that is certainly interesting, I think it would also be interesting to benchmark just the solving step, after the model is already constructed.
1 Like
Interesting to hear this @damaro
I am unable to run the Interval… based code as is, I’m getting UndefVarError: Contractor not defined as an error, and not being familar with the packages, I’m not even sure where Contractor is defined; best guess is IntervalConstraintProgramming ?
Maybe you did not downgrade ModelingToolkit.jl to a v3 version. If you don’t commonly use this package, the easiest option to downgrade is to run the following in the Pkg REPL (pressing ‘]’)
pkg > remove ModelingToolkit
pkg > add ModelingToolkit@3.16.0
pkg > pin ModelingToolkit@3.16.0
Please, let us know if this works.
BTW, both JuMP and ModelingToolkit define the @variables macro. So, I use one Julia session for the IntervalArithmetic method and other one for your suggestion using Gurobi.
Thanks, just needed to update again after the pin.
If you split out the code and just benchmark the solving step, things look fairly similar (at least for this small test case).
See:
Shown below are the times for re-running the scripts through 4 times; the first line is the initial compilation call with @time, the second line with @btime ... samples=10:
- Interval method:
0.742040 seconds (1.79 M allocations: 102.037 MiB, 2.99% gc time, 99.97% compilation time)
13.408 μs (64 allocations: 9.91 KiB)
0.044696 seconds (59.92 k allocations: 3.183 MiB, 99.77% compilation time)
9.921 μs (64 allocations: 9.91 KiB)
0.045063 seconds (59.93 k allocations: 3.184 MiB, 99.72% compilation time)
10.261 μs (64 allocations: 9.91 KiB)
0.045872 seconds (59.94 k allocations: 3.184 MiB, 99.76% compilation time)
12.664 μs (64 allocations: 9.91 KiB)
- Gurobi method
6.713763 seconds (16.86 M allocations: 992.390 MiB, 7.46% gc time, 1.77% compilation time)
25.157 μs (26 allocations: 800 bytes)
0.003869 seconds (1.69 k allocations: 121.092 KiB, 60.37% compilation time)
24.683 μs (26 allocations: 800 bytes)
0.003868 seconds (1.69 k allocations: 121.170 KiB, 60.21% compilation time)
25.832 μs (26 allocations: 800 bytes)
0.004770 seconds (1.69 k allocations: 121.092 KiB, 55.82% compilation time)
24.244 μs (26 allocations: 800 bytes)
Qualifier: I have no idea how each code does caching, so that might be wrinkle in interpreting this.
2 Likes