Improve mid and bisect (#261)

* Add a parameter \alpha to mid, and add some docs of mid, diam, rad

* Add special case for radius of rational interval

* Add tests of mid with parameter

* Add bisect tests for infinite intervals

* Added specialized mid for rational intervals

* Fix test

src/intervals/arithmetic.jl

@@ -307,26 +307,54 @@ end
 # mid, diam, radius
 
 # Compare pg. 64 of the IEEE 1788-2015 standard:
-function mid{T}(a::Interval{T})
+doc"""
+    mid(a::Interval, α=0.5)
+
+Find the midpoint (or, in general, an intermediate point) at a distance α along the interval `a`. The default is the true midpoint at α=0.5.
+"""
+function mid{T}(a::Interval{T}, α=0.5)
+
     isempty(a) && return convert(T, NaN)
     isentire(a) && return zero(a.lo)
 
     a.lo == -∞ && return nextfloat(a.lo)
     a.hi == +∞ && return prevfloat(a.hi)
 
-    (a.lo + a.hi) / 2  # rounds to nearest
+    @assert 0 ≤ α ≤ 1
+
+    return (1-α) * a.lo + α * a.hi  # rounds to nearest
 end
 
+mid{T}(a::Interval{Rational{T}}) = (1//2) * (a.lo + a.hi)
+
+
+doc"""
+    diam(a::Interval)
+
+Return the diameter (length) of the `Interval` `a`.
+"""
 function diam{T<:Real}(a::Interval{T})
     isempty(a) && return convert(T, NaN)
+
     @round_up(a.hi - a.lo) # cf page 64 of IEEE1788
 end
 
+doc"""
+    radius(a::Interval)
+
+Return the radius of the `Interval` `a`, such that
+`a ⊆ m ± radius`, where `m = mid(a)` is the midpoint.
+"""
 # Should `radius` this yield diam(a)/2? This affects other functions!
 function radius(a::Interval)
     isempty(a) && return convert(eltype(a), NaN)
     m = mid(a)
-    max(m - a.lo, a.hi - m)
+    return max(m - a.lo, a.hi - m)
+end
+
+function radius{T}(a::Interval{Rational{T}})
+    m = (a.lo + a.hi) / 2
+    return max(m - a.lo, a.hi - m)
 end
 
 # cancelplus and cancelminus

src/root_finding/bisect.jl

@@ -6,17 +6,20 @@ Split the interval `X` at position α; α=0.5 corresponds to the midpoint.
 Returns a tuple of the new intervals.
 """
 function bisect(X::Interval, α=0.5)
-    m = (1-α) * X.lo + α * X.hi
-    return ( Interval(X.lo, m), Interval(m, X.hi) )
+    @assert 0 ≤ α ≤ 1
+
+    m = mid(X, α)
+
+    return (Interval(X.lo, m), Interval(m, X.hi))
 end
 
 doc"""
     bisect(X::IntervalBox, α=0.5)
 
-Bisect the `IntervalBox` in its longest side.
+Bisect the `IntervalBox` `X` at position α ∈ [0,1] along its longest side.
 """
 function bisect(X::IntervalBox, α=0.5)
-    i = findmax([diam(x) for x in X])[2]  # find longest side
+    i = indmax(diam.(X))  # find longest side
 
     return bisect(X, i, α)
 end

test/interval_tests/consistency.jl

@@ -165,11 +165,17 @@ c = @interval(0.25, 4.0)
         @test mid(1..2) == 1.5
         @test mid(0.1..0.3) == 0.2
         @test mid(-10..5) == -2.5
-        @test mid(-∞..1) == -1.7976931348623157e308
-        @test mid(1..∞) == 1.7976931348623157e308
+        @test mid(-∞..1) == nextfloat(-∞)
+        @test mid(1..∞) == prevfloat(∞)
         @test isnan(mid(emptyinterval()))
     end
 
+    @testset "mid with parameter" begin
+        @test mid(0..1, 0.75) == 0.75
+        @test mid(1..∞, 0.75) == prevfloat(∞)
+        @test mid(-∞..∞, 0.75) == 0
+    end
+
     @testset "diam" begin
 
         @test diam( Interval(1//2) ) == 0//1

test/root_finding_tests/bisect.jl

@@ -3,15 +3,23 @@ using Base.Test
 
 
 @testset "Bisection tests" begin
-    x = Interval(0, 1)
+    X = 0..1
+    @test bisect(X) == (0..0.5, 0.5..1)
+    @test bisect(X, 0.25) == (0..0.25, 0.25..1)
 
-    @test bisect(x) == (0..0.5, 0.5..1)
-    @test bisect(x, 0.25) == (0..0.25, 0.25..1)
+    X = -∞..∞
+    @test bisect(X) == (-∞..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..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) )
+
+    X = (-∞..∞) × (-∞..∞)
+    @test bisect(X) == ( (-∞..0) × (-∞..∞), (0..∞) × (-∞..∞))
 end