# This file is part of the IntervalArithmetic.jl package; MIT licensed

"""
    setdiff(x::Interval, y::Interval)

Calculate the set difference `x ∖ y`, i.e. the set of values
that are inside the interval `x` but not inside `y`.

Returns an array of intervals.
The array may:

- be empty if `x ⊆ y`;
- contain a single interval, if `y` overlaps `x`
- contain two intervals, if `y` is strictly contained within `x`.
"""
function setdiff(x::Interval{T}, y::Interval{T}) where T
    intersection = x ∩ y

    isempty(intersection) && return [x]
    intersection == x && return [emptyinterval(T)]

    if x.lo == intersection.lo
        if isfinite(intersection.hi)
            return [unsafe_interval(intersection.hi, x.hi)]
        else
            return [emptyinterval(T)]
        end
    end
    if x.hi == intersection.hi
        if isfinite(intersection.lo)
            return [unsafe_interval(x.lo, intersection.lo)]
        else
            return [emptyinterval(T)]
        end
    end

    return [unsafe_interval(x.lo, y.lo),
            unsafe_interval(y.hi, x.hi)]
end
setdiff(x::ValidInterval, ::ValidInterval) = [x]


"""
    distance(a::Interval, b::Interval)

Hausdorff distance between two intervals.
"""
distance(a::Interval, b::Interval) = max(abs(a.lo-b.lo), abs(a.hi-b.hi))

const where_bisect = 0.49609375

"""
    bisect(X::Interval, α=0.49609375)

Split the interval `X` at position α; α=0.5 corresponds to the midpoint.
Returns a tuple of the new intervals.
"""
function bisect(X::Interval, α=where_bisect)
    @assert 0 ≤ α ≤ 1

    m = mid(X, α)

    return (unsafe_interval(X.lo, m), unsafe_interval(m, X.hi))
end

"""
    mid(a::Interval, α=0.5)

Find an intermediate point at a relative position `α`` in the interval `a`.
The default is the true midpoint at `α = 0.5`.

Assumes 0 ≤ α ≤ 1.

Warning: if the parameter `α = 0.5` is explicitly set, the behavior differs
from the default case if the provided `Interval` is not finite, since when
`α` is provided `mid` simply replaces `+Inf` (respectively `-Inf`) by `prevfloat(+Inf)`
(respecively `nextfloat(-Inf)`) for the computation of the intermediate point.
"""
function mid(a::Interval{T}, α) where T

    lo = (a.lo == -Inf ? nextfloat(T(-Inf)) : a.lo)
    hi = (a.hi == +Inf ? prevfloat(T(+Inf)) : a.hi)

    β = T(α)

    midpoint = β * (hi - lo) + lo
    isfinite(midpoint) && return midpoint
    #= Fallback in case of overflow: hi - lo == +Inf.
       This case can not be the default one as it does not pass several
       IEEE1788-2015 tests for small floats.
    =#
    return (1 - β) * lo + β * hi
end

mid(::EmptyInterval{T}, α) where T = T(NaN)
