integer_solutions.jl

using IntervalConstraintProgramming
using IntervalArithmetic
using ModelingToolkit

"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)
    
X = (1.2..2.1) × (3.4..7.2)
integerize(X)

X = (1.2..1.9) × (3.4..7.2)
integerize(X)
isempty(ans)



## Simple example

vars = @variables x, y, u, v

contractors = [X -> C(0..0, X) for C in Contractor.(Ref(vars), [x + y - 3, u + v - 1, x + u - 2, y + v - 2])]
contractors = [contractors; integerize]


X = IntervalBox(0..4, 4)

solutions = pave(contractors, [X], 0.1)

all(diam.(solutions) .== 0)

[Int.(x) for x in mid.(solutions)]


## Complicated example

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

n = length(μ)
vars = @variables y[1:n, 1:n]
vars = vars[1]

vars

contractors = []

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(0..n^2, n^2)

contractors
contractors = [X -> C(0..0, X) for C in contractors]

for C in contractors
    global X = C(X)
end

X

contractors = [contractors; integerize]

solutions = pave(contractors, [X], 1.0)

all(diam.(solutions) .== 0)

[Int.(x) for x in mid.(solutions)]