gurobi_soution

## Packages
using JuMP
using Gurobi

using BenchmarkTools

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

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

@btime example_big_sum($μ, $ν; verbose=false) samples=10

interval_arithmetic_solution

using IntervalConstraintProgramming
using IntervalArithmetic
using ModelingToolkit

using BenchmarkTools

"Cs is a list of contractors"
function pave(Cs, working::Vector{IntervalBox{N,T}}, ϵ, bisection_point=nothing) where {N,T}

    boundary = SubPaving{N,T}()

    while !isempty(working)

        X = pop!(working)

        # apply each contractor in a round-robin fashion:
        for C in Cs
            X = C(X)
            
            if isempty(X)
                break
            end
        end
        
        if isempty(X)
            continue
        end

        if diam(X) < ϵ
            push!(boundary, X)

        else
            if bisection_point == nothing
                push!(working, bisect(X)...)
            else
                push!(working, bisect(X, bisection_point)...)
            end

        end

    end

    return boundary

end

function integerize(X::Interval)
    a = ceil(X.lo)
    b = floor(X.hi)
    
    if a > b
        return emptyinterval(X)
    end

    return Interval(a, b)
end 

integerize(X::IntervalBox) = integerize.(X)

function generate_contractors(μ::Array{Int64,1}, ν::Array{Int64,1})
    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_solution(μ::Array{Int64,1}, ν::Array{Int64,1})

    # intervalo y contractors
    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

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

@btime find_solution($μ,$ν) samples=10