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.
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.
@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
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:
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.
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
:
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)
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.