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

So that’s in 16 dimensions? That is indeed rather high. I can’t look at this right now but will try later

Here’s a gist with a better approach:

Basically integerize needs to be mixed in throughout the process, and we use Contractors instead of Separators.

(Separators are used when you want to find inner and outer approximations. We want only an outer approximation here.)

This new version produces the following results in a few seconds:

24-element Array{StaticArrays.SArray{Tuple{16},Int64,1,16},1}:
[0, 2, 0, 0, 0, 0, 2, 0, 0, 0, 1, 0, 0, 0, 0, 1]
[0, 2, 0, 0, 0, 0, 2, 0, 0, 0, 0, 1, 0, 0, 1, 0]
[0, 2, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0]
[0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1]
[0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0]
⋮
[0, 0, 1, 1, 0, 2, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0]
[0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0]
[0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0]
[0, 0, 2, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0]
[0, 0, 1, 1, 0, 0, 2, 0, 0, 1, 0, 0, 0, 1, 0, 0]


Do you know what the correct solutions are in this case?

1 Like

By the way, after applying each contractor once to the initial box [0..16]^16, it already reduces to

[0, 0] × [0, 2] × [0, 2] × [0, 1] × [0, 0] × [0, 2] × [0, 2] × [0, 1] × [0, 0] × [0, 1] × [0, 1] × [0, 1] × [0, 0] × [0, 1] × [0, 1] × [0, 1]


This really shows the power of constraint propagation!

I forgot to say that the world-age error comes from code generation inside Separator / Contractor.

If you want to get around this while still being inside a function, you need something like

Base.invokelatest(pave, C, X)

There are the 24 solutions (returned by Mathematica):

0,0,0,0,0,0,1,1,1,2,0,0,1,0,0,0
0,0,0,0,0,0,1,1,2,1,0,0,0,1,0,0
0,0,0,0,0,1,0,1,1,1,1,0,1,0,0,0
0,0,0,0,0,1,0,1,2,0,1,0,0,1,0,0
0,0,0,0,0,1,0,1,2,1,0,0,0,0,1,0
0,0,0,0,0,1,1,0,1,1,0,1,1,0,0,0
0,0,0,0,0,1,1,0,2,0,0,1,0,1,0,0
0,0,0,0,0,1,1,0,2,1,0,0,0,0,0,1
0,0,0,0,0,2,0,0,1,0,1,1,1,0,0,0
0,0,0,0,0,2,0,0,2,0,0,1,0,0,1,0
0,0,0,0,0,2,0,0,2,0,1,0,0,0,0,1
0,0,0,0,1,0,0,1,0,2,1,0,1,0,0,0
0,0,0,0,1,0,0,1,1,1,1,0,0,1,0,0
0,0,0,0,1,0,0,1,1,2,0,0,0,0,1,0
0,0,0,0,1,0,1,0,0,2,0,1,1,0,0,0
0,0,0,0,1,0,1,0,1,1,0,1,0,1,0,0
0,0,0,0,1,0,1,0,1,2,0,0,0,0,0,1
0,0,0,0,1,1,0,0,0,1,1,1,1,0,0,0
0,0,0,0,1,1,0,0,1,0,1,1,0,1,0,0
0,0,0,0,1,1,0,0,1,1,0,1,0,0,1,0
0,0,0,0,1,1,0,0,1,1,1,0,0,0,0,1
0,0,0,0,2,0,0,0,0,1,1,1,0,1,0,0
0,0,0,0,2,0,0,0,0,2,0,1,0,0,1,0
0,0,0,0,2,0,0,0,0,2,1,0,0,0,0,1


I tried checking both solutions (yours and Mathematica) using sort with no avail. I’ll send the results soon, but I’m confident both are equal.

It looks like the variables are stored in a different order, which is making things difficult to compare?

Both solutions are equal.
I exported your solution into an csv and compared it against Mathematica solution. Unfortunately, I have been unable to compare it using Julia.
I tried the following
sortslices(mma_output,dims=1)-sortslices(hcat(collect(output)...)',dims=1), where mma_output was obtained by mma_output = readdlm("/home/david/example.csv",',',Int).
I was expecting to find a zero matrix but, as you can check, it is not the case. example.csv is a file containing the list in post 14.

What order is Mathematica using for the variables?
I tried just transposing the matrix of variables in the Julia code, but that didn’t help.

I would just apply a permutation matrix that maps the order to the correct one.
I’m not sure that just sorting is enough.

The orden I gave to Mathematica is

# non Julia code
{x[1, 1], x[2, 1], x[3, 1], x[4, 1], x[1, 2], x[2, 2], x[3, 2],
x[4, 2], x[1, 3], x[2, 3], x[3, 3], x[4, 3], x[1, 4], x[2, 4],
x[3, 4], x[4, 4]}


Inspecting vars, it looks like this is the same order as in Julia. Still, I am unable to match both solutions in Julia.

Edit. Silly mistake. I interchanged \nu and \mu. Now running sortslices(mma_output,dims=1)-sortslices(hcat(collect(output)...)',dims=1) results in the zero matrix.

1 Like

Here’s a way to force Julia to use exactly that order of the variables. Then the output looks the same:

vars = @variables y[1:n, 1:n]
vars = vars[1]
vars = permutedims(vars, (1, 2))

contractors = []

for i in 1:n
push!(contractors, Contractor(vec(vars), sum(vars[i, 1:n]) - μ[i]))
end

for i in 1:n
push!(contractors, Contractor(vec(vars), sum(vars[1:n, i]) - ν[i]))
end

2 Likes

Hi @damaro, here is how you do this in Gurobi (I’m using Gurobi 0.9.12, JuMP 0.21.7)

## Packages
using JuMP
using Gurobi

## This function will return an array of named tuples in the variables with objective value.
function example_little_sum(; verbose::Bool = true)
## Model
model = Model(Gurobi.Optimizer)

@variable(model, 0 <= x <= 4, Int)
@variable(model, 0 <= y <= 4, Int)
@variable(model, 0 <= u <= 4, Int)
@variable(model, 0 <= v <= 4, Int)
@constraint(model, x + y == 3)
@constraint(model, u + v == 1)
@constraint(model, x + u == 2)
@constraint(model, y + v == 2)
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)
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("Num of results: ", num_results); end

Results = Array{NamedTuple{(:x,:y,:u,:v,:obj),NTuple{5,Int64}},1}() ## Note the conversion to Int64
for i in 1:num_results
sol = (x=JuMP.value(x; result=i), y=JuMP.value(y; result=i),
u=JuMP.value(u; result=i), v=JuMP.value(v; result=i),
obj=JuMP.objective_value(model; result = i))
push!(Results,sol)
if verbose
println(Results[i])
end
end

return Results
end

3 Likes

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