interval_arithmetic_solution_split

interval_arithmetic_solution.jl

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 = (ModelingToolkit.@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 setup_solve(μ::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]

    return contractors, X
end

function find_solution(contractors::Vector{Function},X::IntervalBox{N, Float64}) where N
    solutions = Base.invokelatest(pave, contractors,[X], 1.0)
    # return [Int.(x) for x in mid.(solutions)]
end

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

contractors, X = setup_solve(μ,ν)
solutions = @time find_solution(contractors,X)
@btime find_solution(contractors,X) samples=10