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

**URL:** <https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709>\
**Category:** Optimization (Mathematical)\
**Tags:** question, jump\
**Created:** [August 10, 2020, 9:44pm UTC](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709 "2020-08-10T21:44:52Z")\
**Posts on this page:** 9\
**Page:** 2

<div class="post-metadata">

**Author:** ![damaro](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/damaro/32/17002_2.png) [@damaro](https://discourse.julialang.org/u/damaro)\
**Post date:** [May 22, 2021, 1:45am UTC](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709/21 "2021-05-22T01:45:04Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![jd-foster](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jd-foster/32/35824_2.png) [@jd-foster](https://discourse.julialang.org/u/jd-foster)\
**Post date:** [May 24, 2021, 12:58pm UTC](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709/22 "2021-05-24T12:58:36Z")

</div>

No problem.  
While it’s fresh in my mind, here is a more general solution for future reference:

```julia
## 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.

---

<div class="post-metadata">

**Author:** ![damaro](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/damaro/32/17002_2.png) [@damaro](https://discourse.julialang.org/u/damaro)\
**Post date:** [June 3, 2021, 7:44pm UTC](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709/23 "2021-06-03T19:44:37Z")

</div>

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

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

```julia
pkg > remove ModelingToolkit
pkg > add ModelingToolkit@3.16.0
pkg > pin ModelingToolkit@3.16.0

```

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [June 3, 2021, 10:23pm UTC](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709/25 "2021-06-03T22:23:20Z")

</div>

> [@damaro](#):
>
> @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.

What are the respective timings and what is the benchmarking code? Can you post a complete code that can be copied and pasted?

---

<div class="post-metadata">

**Author:** ![damaro](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/damaro/32/17002_2.png) [@damaro](https://discourse.julialang.org/u/damaro)\
**Post date:** [June 3, 2021, 10:31pm UTC](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709/26 "2021-06-03T22:31:09Z")

</div>

I am far from being a expert benchmarking, but here is what I did using `@btime` from `BenchmarkingTools`:

```julia
μ = [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:

```julia
@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](https://gist.github.com/davidamaro/b2ae3a106f7e8385d0e28f3dc6bffb9a)

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [June 3, 2021, 11:51pm UTC](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709/27 "2021-06-03T23:51:08Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![jd-foster](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jd-foster/32/35824_2.png) [@jd-foster](https://discourse.julialang.org/u/jd-foster)\
**Post date:** [June 4, 2021, 1:38am UTC](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709/28 "2021-06-04T01:38:08Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![damaro](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/damaro/32/17002_2.png) [@damaro](https://discourse.julialang.org/u/damaro)\
**Post date:** [June 4, 2021, 2:39am UTC](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709/29 "2021-06-04T02:39:10Z")

</div>

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 ‘]’)

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

---

<div class="post-metadata">

**Author:** ![jd-foster](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jd-foster/32/35824_2.png) [@jd-foster](https://discourse.julialang.org/u/jd-foster)\
**Post date:** [June 4, 2021, 5:03am UTC](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709/30 "2021-06-04T05:03:25Z")

</div>

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:

- [Interval method gist](https://gist.github.com/jd-foster/40ed9bcb18a5d919945bed3850e99bcb)
- [Gurobi method gist](https://gist.github.com/jd-foster/018d45f4f5ef7c1e1de870804ff7dc14)

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

```julia
  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**

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

[Previous page](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709.md?page=1)