diff --git a/src/intervals/arithmetic.jl b/src/intervals/arithmetic.jl index 5efc14d..f3e9166 100644 --- a/src/intervals/arithmetic.jl +++ b/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 diff --git a/src/root_finding/bisect.jl b/src/root_finding/bisect.jl index 104c5c4..93e19a2 100644 --- a/src/root_finding/bisect.jl +++ b/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 diff --git a/test/interval_tests/consistency.jl b/test/interval_tests/consistency.jl index 7e3550b..d541986 100644 --- a/test/interval_tests/consistency.jl +++ b/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 diff --git a/test/root_finding_tests/bisect.jl b/test/root_finding_tests/bisect.jl index 1f1159e..128be08 100644 --- a/test/root_finding_tests/bisect.jl +++ b/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