 # Right solver for JuMP to find every solution of a linear system of equations with integer solutions

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.

6 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])

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

Code for benchmarking

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