Add methods for IntervalBox previously in IntervalRootFinding (#160)

* Add methods for IntervalBox previously in IntervalRootFinding

* Remove a reference to IntervalRootFinding, and fix complex interval test

* Implement suggestions of the review

* Incorporate a missing suggestion of the review

Author: lbenet

src/IntervalArithmetic.jl

@@ -18,13 +18,13 @@ import Base:
     asin, acos, atan, atan2,
     sinh, cosh, tanh, asinh, acosh, atanh,
     union, intersect, isempty,
-    convert, promote_rule, eltype,
+    convert, promote_rule, eltype, size,
     BigFloat, float, widen, big,
     ∩, ∪, ⊆, eps,
     floor, ceil, trunc, sign, round,
     expm1, log1p,
     precision,
-    isfinite, isnan,
+    isfinite, isnan, isinf, iszero,
     show, showall,
     isinteger, setdiff,
     parse
@@ -40,9 +40,9 @@ export
     @interval, @biginterval, @floatinterval, @make_interval,
     diam, radius, mid, mag, mig, hull,
     emptyinterval, ∅, ∞, isempty, isinterior, isdisjoint, ⪽,
-    precedes, strictprecedes, ≺, ⊂, ⊃, ⊇,
+    precedes, strictprecedes, ≺, ⊂, ⊃, ⊇, contains_zero,
     entireinterval, isentire, nai, isnai, isthin, iscommon, isatomic,
-    widen, inf, sup,
+    widen, inf, sup, bisect, 
     parameters, eps, dist,
     pi_interval,
     midpoint_radius, interval_from_midpoint_radius,
@@ -88,6 +88,7 @@ end
 
 include("intervals/intervals.jl")
 include("multidim/multidim.jl")
+include("bisect.jl")
 include("decorations/decorations.jl")
 
 include("parsing.jl")

src/bisect.jl

@@ -0,0 +1,42 @@
+
+const where_bisect = 0.49609375
+
+doc"""
+    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 (Interval(X.lo, m), Interval(m, X.hi))
+end
+
+doc"""
+    bisect(X::IntervalBox, α=0.49609375)
+
+Bisect the `IntervalBox` `X` at position α ∈ [0,1] along its longest side.
+"""
+function bisect(X::IntervalBox, α=where_bisect)
+    i = indmax(diam.(X))  # find longest side
+
+    return bisect(X, i, α)
+end
+
+doc"""
+    bisect(X::IntervalBox, i::Integer, α=0.49609375)
+
+Bisect the `IntervalBox` in side number `i`.
+"""
+function bisect(X::IntervalBox, i::Integer, α=where_bisect)
+
+    x1, x2 = bisect(X[i], α)
+
+    X1 = setindex(X, x1, i)
+    X2 = setindex(X, x2, i)
+
+    return (X1, X2)
+end

src/intervals/intervals.jl

@@ -58,6 +58,9 @@ Interval{T}(x) where T = Interval(convert(T, x))
 
 Interval{T}(x::Interval) where T = atomic(Interval{T}, x)
 
+size(x::Interval) = (1,)
+
+
 """
     is_valid_interval(a::Real, b::Real)
 

src/intervals/special.jl

@@ -22,6 +22,7 @@ entireinterval(x::Interval{T}) where T<:Real = entireinterval(T)
 entireinterval() = entireinterval(precision(Interval)[1])
 
 isentire(x::Interval) = x.lo == -Inf && x.hi == Inf
+isinf(x::Interval) = isentire(x)
 isunbounded(x::Interval) = x.lo == -Inf || x.hi == Inf
 
 
@@ -75,7 +76,10 @@ This occurs when the interval is empty, or when the upper bound equals the lower
 """
 isatomic(x::Interval) = isempty(x) || (x.hi == x.lo) || (x.hi == nextfloat(x.lo))
 
-Base.iszero(x::Interval) = iszero(x.lo) && iszero(x.hi)
+iszero(x::Interval) = iszero(x.lo) && iszero(x.hi)
+
+contains_zero(X::Interval{T}) where {T} = zero(T) ∈ X
+
 
 # """
 #     widen(x)

src/multidim/intervalbox.jl

@@ -73,6 +73,12 @@ diam(X::IntervalBox) = maximum(diam.(X.v))
 
 emptyinterval(X::IntervalBox{N,T}) where {N,T} = IntervalBox(emptyinterval.(X.v))
 
+isinf(X::IntervalBox) = any(isinf(X))
+
+isinterior(X::IntervalBox{N,T}, Y::IntervalBox{N,T}) where {N,T} = all(isinterior.(X, Y))
+
+contains_zero(X::SVector) = all(contains_zero(X))
+contains_zero(X::IntervalBox) = all(contains_zero(X))
 
 import Base.×
 ×(a::Interval...) = IntervalBox(a...)

test/interval_tests/bisect.jl

@@ -0,0 +1,30 @@
+using IntervalArithmetic
+using Base.Test
+
+
+@testset "`bisect` function" begin
+    X = 0..1
+    @test bisect(X, 0.5) == (0..0.5, 0.5..1)
+    @test bisect(X, 0.25) == (0..0.25, 0.25..1)
+
+    @test bisect(X) == (interval(0.0, 0.49609375), interval(0.49609375, 1.0))
+
+    X = -∞..∞
+    @test bisect(X, 0.5) == (-∞..0, 0..∞)
+    @test bisect(X, 0.75) == (-∞..0, 0..∞)
+
+    X = 1..∞
+    @test bisect(X) == (Interval(1, prevfloat(∞)), Interval(prevfloat(∞), ∞))
+
+    X = (0..1) × (0..2)
+    @test bisect(X, 0.5) == ( (0..1) × (0..1), (0..1) × (1..2) )
+    @test bisect(X, 0.25) == ( (0..1) × (0..0.5), (0..1) × (0.5..2) )
+    @test bisect(X, 1, 0.5) == ( (0..0.5) × (0..2), (0.5..1) × (0..2) )
+    @test bisect(X, 1, 0.25) == ( (0..0.25) × (0..2), (0.25..1) × (0..2) )
+
+    @test bisect(X) ==  (IntervalBox(0..1, interval(0.0, 0.9921875)),
+                         IntervalBox(0..1, Interval(0.9921875, 2.0)))
+
+    X = (-∞..∞) × (-∞..∞)
+    @test bisect(X) == ( (-∞..0) × (-∞..∞), (0..∞) × (-∞..∞))
+end

test/interval_tests/complex.jl

@@ -14,7 +14,6 @@ end
 
 @testset "Complex functions" begin
     Z = (3 ± 1e-7) + (4 ± 1e-7)*im
-    @test sin(Z) == Interval(3.853734949309744, 3.8537411265295543) - Interval(27.016810169394066, 27.016816346613904)*im
-
+    @test sin(Z) == complex(sin(real(Z))*cosh(imag(Z)),sinh(imag(Z))*cos(real(Z)))
     @test exp(-im * Interval(π)) == Interval(-1.0, -0.9999999999999999) - Interval(1.224646799147353e-16, 1.2246467991473532e-16)*im
 end

test/interval_tests/construction.jl

@@ -22,6 +22,7 @@ using Base.Test
 
     # Naive constructors, with no conversion involved
     @test Interval(1) == Interval(1.0, 1.0)
+    @test size(Interval(1)) == (1,)
     @test Interval(big(1)) == Interval(1.0, 1.0)
     @test Interval(eu) == Interval(1.0*eu)
     @test Interval(1//10) == Interval{Rational{Int}}(1//10, 1//10)

test/interval_tests/intervals.jl

@@ -11,5 +11,5 @@ include("loops.jl")
 include("parsing.jl")
 include("rounding_macros.jl")
 include("rounding.jl")
-
-# include("complex.jl")
+include("bisect.jl")
+include("complex.jl")