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

I’m not positive now.

I tried @dpsanders’s solution on a larger —but still on topic—
problem and found a huge use of resources. Example below.

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

components = 4 #  μ |> length
vars = @variables y[1:components,1:components]

constriction_list = Separator[]
for i in 1:components
    push!(constriction_list, Separator(vars..., sum(vars[1][i,1:components]) == μ[i]))
end
for i in 1:components
    push!(constriction_list, Separator(vars..., sum(vars[1][1:components,i]) == ν[i]))
end

X = IntervalBox(0..(components*components), (components*components))  # domain
for c1 in constriction_list
    X = c1(X)[1]
end

ini = ∩(constriction_list[1], constriction_list[2])
for ci in constriction_list[3:end]
    ini = ∩(ini, ci)
end
p = pave(ini, X, 1e-2)

unique(integerize.(p.boundary))

So, is this (large memory usage) an issue with the way I extended the solution or indeed IntervalConstraintProgramming
is not apt to this new problem?

ps1. I was unable to use a function for the previous code. I had a problem similar to The applicable method may be too new: running in world age 28047, while current world is 28055.
ps2. In case someone wants to copy-paste code from this post, you need to add to @dpsanders’s code using IntervalArithmetic

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:

https://gist.github.com/dpsanders/5cf6448b3ef3d280b6d70a8fefe34963

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

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