Rework roots interface to allow choosing a strategy

- Introduce default strategies
- Rework interface to allow using them
- Changes break tests for complex Krawczyk method which expectation has been lowered

src/roots.jl

@@ -22,8 +22,8 @@ Inputs:
     values are of type `Bisection` or `Newton`
 
 """
-function branch_and_prune(X, contractor, tol=1e-3)
-    iter = RootSearch(X, contractor, tol)
+function branch_and_prune(X, contractor, strategy, tol)
+    iter = RootSearch(X, contractor, strategy, tol)
     local state
     # complete iteration
     for state in iter end
@@ -54,100 +54,162 @@ IntervalLike{T} = Union{Interval{T}, IntervalBox{T}}
 NewtonLike = Union{Type{Newton}, Type{Krawczyk}}
 
 """
-    roots(f, X, contractor, tol=1e-3)
-    roots(f, deriv, X, contractor, tol=1e-3)
+    roots(f, X, contractor=Newton, strategy=BreadthFirstSearch, tol=1e-15)
+    roots(f, deriv, X, contractor=Newton, strategy=BreadthFirstSearch, tol=1e-15)
+    roots(f, X, contractor, tol)
+    roots(f, deriv, X, contractor, tol)
 
-Uses a generic branch and prune routine to find in principle all isolated roots of a function
-`f:R^n → R^n` in a box `X`, or a vector of boxes.
+Uses a generic branch and prune routine to find in principle all isolated roots
+of a function `f:R^n → R^n` in a region `X`, if the number of roots is finite.
 
 Inputs:
 - `f`: function whose roots will be found
-- `X`: `Interval` or `IntervalBox`
+- `X`: `Interval` or `IntervalBox` in which roots are searched
 - `contractor`: function that, when applied to the function `f`, determines
     the status of a given box `X`. It returns the new box and a symbol indicating
     the status. Current possible values are `Bisection`, `Newton` and `Krawczyk`
 - `deriv`: explicit derivative of `f` for `Newton` and `Krawczyk`
+- `strategy`: `SearchStrategy` determining the order in which regions are
+    processed.
+- `tol`: Absolute tolerance. If a region has a diameter smaller than `tol`, it
+    is returned with status `:unkown`.
 
 """
-# Contractor specific `roots` functions
-function roots(f, X::IntervalLike{T}, ::Type{Bisection}, tol::Float64=1e-3) where {T}
-    branch_and_prune(X, Bisection(f), tol)
+
+const default_strategy = DepthFirstSearch
+const default_tolerance = 1e-15
+const default_contractor = Newton
+
+#===
+    Default case when `contractor, `strategy` or `tol` is omitted.
+===#
+function roots(f::Function, X, contractor::Type{C}=default_contractor,
+               strategy::SearchStrategy=default_strategy,
+               tol::Float64=default_tolerance) where {C <: Contractor}
+
+    _roots(f, X, contractor, strategy, tol)
 end
 
-function roots(f, X::Interval{T}, C::NewtonLike, tol::Float64=1e-3) where {T}
+function roots(f::Function, deriv::Function, X, contractor::Type{C}=default_contractor,
+               strategy::SearchStrategy=default_strategy,
+               tol::Float64=default_tolerance) where {C <: Contractor}
 
-    deriv = x -> ForwardDiff.derivative(f, x)
+    _roots(f, deriv, X, contractor, strategy, tol)
+end
+
+function roots(f::Function, X, contractor::Type{C},
+               tol::Float64) where {C <: Contractor}
 
-    roots(f, deriv, X, C, tol)
+    _roots(f, X, contractor, default_strategy, tol)
 end
 
-function roots(f, deriv, X::Interval{T}, C::NewtonLike, tol::Float64=1e-3) where {T}
-    branch_and_prune(X, C(f, deriv), tol)
+function roots(f::Function, deriv::Function, X, contractor::Type{C},
+               tol::Float64) where {C <: Contractor}
+
+    _roots(f, deriv, X, contractor, default_strategy, tol)
 end
 
-function roots(f, X::IntervalBox{T}, C::NewtonLike, tol::Float64=1e-3) where {T}
+#===
+    More specific `roots` methods (all parameters are present)
+    These functions are called `_roots` to avoid recursive calls.
+===#
 
-    deriv = x -> ForwardDiff.jacobian(f, x)
+# For `Bisection` method
+function _roots(f, X::IntervalLike{T}, ::Type{Bisection},
+               strategy::SearchStrategy, tol::Float64) where {T}
 
-    roots(f, deriv, X, C, tol)
+    branch_and_prune(X, Bisection(f), strategy, tol)
 end
 
-function roots(f, deriv, X::IntervalBox{T}, C::NewtonLike, tol::Float64=1e-3) where {T}
-    branch_and_prune(X, C(f, deriv), tol)
+
+# For `NewtonLike` acting on `Interval`
+function _roots(f, X::Interval{T}, contractor::NewtonLike,
+               strategy::SearchStrategy, tol::Float64) where {T}
+
+    deriv = x -> ForwardDiff.derivative(f, x)
+    _roots(f, deriv, X, contractor, strategy, tol)
+end
+
+function _roots(f, deriv, X::Interval{T}, contractor::NewtonLike,
+               strategy::SearchStrategy, tol::Float64) where {T}
+
+    branch_and_prune(X, contractor(f, deriv), strategy, tol)
+end
+
+
+# For `NewtonLike` acting on `IntervalBox`
+function _roots(f, X::IntervalBox{T}, contractor::NewtonLike,
+               strategy::SearchStrategy, tol::Float64) where {T}
+
+    deriv = x -> ForwardDiff.jacobian(f, x)
+    _roots(f, deriv, X, contractor, strategy, tol)
 end
 
-roots(f, r::Root, contractor::Type{C}, tol::Float64=1e-3) where {C<:Contractor} =
-    roots(f, r.interval, contractor, tol)
-roots(f, deriv, r::Root, contractor::Type{C}, tol::Float64=1e-3) where {C<:Contractor} =
-    roots(f, deriv, r.interval, contractor, tol)
+function _roots(f, deriv, X::IntervalBox{T}, contractor::NewtonLike,
+               strategy::SearchStrategy, tol::Float64) where {T}
 
+    branch_and_prune(X, contractor(f, deriv), strategy, tol)
+end
 
-# Acting on a Vector:
 
+# Acting on`Root`
 # TODO: Use previous status information about roots:
-roots(f, V::Vector{Root{T}}, contractor::Type{C}, tol::Float64=1e-3) where {T, C<:Contractor} =
-    vcat(roots.(f, V, contractor, tol)...)
-roots(f, deriv, V::Vector{Root{T}}, contractor::Type{C}, tol::Float64=1e-3) where {T, C<:Contractor} =
-    vcat(roots.(f, deriv, V, contractor, tol)...)
+function _roots(f, r::Root, contractor::Type{C},
+               strategy::SearchStrategy, tol::Float64) where {C<:Contractor}
+
+    _roots(f, r.interval, contractor, strategy, tol)
+end
+
+function _roots(f, deriv, r::Root, contractor::Type{C},
+               strategy::SearchStrategy, tol::Float64) where {C<:Contractor}
+
+    _roots(f, deriv, r.interval, contractor, strategy, tol)
+end
+
 
+# Acting on `Vector` of `Root`
+function _roots(f, V::Vector{Root{T}}, contractor::Type{C},
+               strategy::SearchStrategy, tol::Float64) where {T, C<:Contractor}
 
+    vcat(_roots.(f, V, contractor, strategy, tol)...)
+end
+
+function _roots(f, deriv, V::Vector{Root{T}}, contractor::Type{C},
+               strategy::SearchStrategy, tol::Float64) where {T, C<:Contractor}
 
-# Complex:
+    vcat(_roots.(f, deriv, V, contractor, strategy, tol)...)
+end
 
-function roots(f, Xc::Complex{Interval{T}}, contractor::Type{C},
-        tol::Float64=1e-3) where {T, C<:Contractor}
+
+# Acting on complex `Interval`
+function _roots(f, Xc::Complex{Interval{T}}, contractor::Type{C},
+               strategy::SearchStrategy, tol::Float64) where {T, C<:Contractor}
 
     g = realify(f)
     Y = IntervalBox(reim(Xc))
-    rts = roots(g, Y, contractor, tol)
+    rts = _roots(g, Y, contractor, strategy, tol)
 
     return [Root(Complex(root.interval...), root.status) for root in rts]
 end
 
-function roots(f, Xc::Complex{Interval{T}}, C::NewtonLike, tol::Float64=1e-3) where {T}
+function _roots(f, Xc::Complex{Interval{T}}, contractor::NewtonLike,
+               strategy::SearchStrategy, tol::Float64) where {T}
 
     g = realify(f)
-
     g_prime = x -> ForwardDiff.jacobian(g, x)
-
     Y = IntervalBox(reim(Xc))
-    rts = roots(g, g_prime, Y, C, tol)
+    rts = _roots(g, g_prime, Y, contractor, strategy, tol)
 
     return [Root(Complex(root.interval...), root.status) for root in rts]
 end
 
-function roots(f, deriv, Xc::Complex{Interval{T}}, C::NewtonLike, tol::Float64=1e-3) where {T}
+function _roots(f, deriv, Xc::Complex{Interval{T}}, contractor::NewtonLike,
+               strategy::SearchStrategy, tol::Float64) where {T}
 
     g = realify(f)
-
     g_prime = realify_derivative(deriv)
-
     Y = IntervalBox(reim(Xc))
-    rts = roots(g, g_prime, Y, C, tol)
+    rts = _roots(g, g_prime, Y, contractor, strategy, tol)
 
     return [Root(Complex(root.interval...), root.status) for root in rts]
 end
-
-# Default
-roots(f::F, deriv::F, X, tol::Float64=1e-15) where {F<:Function} = roots(f, deriv, X, Newton, tol)
-roots(f::F, X, tol::Float64=1e-15) where {F<:Function} = roots(f, X, Newton, tol)

src/rootsearch_iterator.jl

@@ -1,6 +1,7 @@
 import Base: start, next, done, copy, eltype, iteratorsize
 
 export RootSearch, SearchStrategy
+export BreadthFirstSearch, DepthFirstSearch
 export start, next, done, copy, step!, eltype, iteratorsize
 
 """
@@ -24,7 +25,8 @@ end
 
 SearchStrategy(CON::Type, store!::Function, retrieve!::Function) = SearchStrategy{CON}(store!, retrieve!)
 
-SearchStrategy() = SearchStrategy(Vector, push!, pop!)
+const BreadthFirstSearch = SearchStrategy{Vector}(push!, shift!)
+const DepthFirstSearch = SearchStrategy{Vector}(push!, pop!)
 
 """
     RootSearch{R <: Union{Interval,IntervalBox}, C <: Contractor, CON, T <: Real}
@@ -48,10 +50,6 @@ struct RootSearch{R <: Union{Interval,IntervalBox}, C <: Contractor, CON, T <: R
     tolerance::T
 end
 
-function RootSearch(region::R, contractor::C, tol::T) where {R <: Union{Interval,IntervalBox}, C <: Contractor, T <: Real}
-    RootSearch(region, contractor, SearchStrategy(), tol)
-end
-
 eltype(::Type{RS}) where {R, C, T, CON, RS <: RootSearch{R, C, CON, T}} = RootSearchState{CON{R}, CON{Root{R}}}
 iteratorsize(::Type{RS}) where {RS <: RootSearch} = Base.SizeUnknown()
 
@@ -79,6 +77,7 @@ function RootSearchState(region::R, strat::SearchStrategy{CON}) where {R <: Unio
     return RootSearchState(working, outputs)
 end
 
+# Currently never reached
 function RootSearchState(region::T) where {T<:Union{Interval,IntervalBox}}
     working = [region]
     outputs = Root{T}[]

test/roots.jl

@@ -6,11 +6,39 @@ function all_unique(rts)
     all(root_status.(rts) .== :unique)
 end
 
-function test_newtonlike(f, deriv, X, method, nsol, tol=1e-3)
+function reldiff(X, Y)
+    if diam(X) == 0
+        return 0
+    else
+        return (abs(X.lo - Y.lo) + abs(X.lo - Y.lo))/diam(X)
+    end
+end
+
+function reldiff(X::Root{T}, Y::Root{T}) where {T}
+    return reldiff(X.interval, Y.interval)
+end
+
+function reldiff(X::Root{T}, Y::Root{T}) where {T <: IntervalBox}
+    s = 0
+    for (x, y) in zip(X.interval, Y.interval)
+        s += reldiff(x, y)
+    end
+    return s
+end
+
+function reldiff(X::Root{T}, Y::Root{T}) where {T <: Complex}
+    return reldiff(X.interval.re, Y.interval.re) + reldiff(X.interval.im, Y.interval.im)
+end
+
+function test_newtonlike(f, deriv, X, method, nsol, reltol=0)
     rts = roots(f, X, method)
     @test length(rts) == nsol
     @test all_unique(rts)
-    @test rts == roots(f, deriv, X, method)
+    if reltol == 0
+        @test rts == roots(f, deriv, X, method)
+    else
+        @test all(reldiff.(rts, roots(f, deriv, X, method)) .< reltol)
+    end
 end
 
 newtonlike_methods = [Newton, Krawczyk]
@@ -22,7 +50,7 @@ newtonlike_methods = [Newton, Krawczyk]
     @test all_unique(rts)
 
     # Bisection
-    rts = roots(sin, -5..6, Bisection)
+    rts = roots(sin, -5..6, Bisection, 1e-3)
     @test length(rts) == 3
 
     # Refinement
@@ -111,13 +139,13 @@ end
 
     for method in newtonlike_methods
         deriv = z -> 3*z^2
-        test_newtonlike(f, deriv, Xc, method, 3)
+        test_newtonlike(f, deriv, Xc, method, 3, 1e-10)
     end
 end
 
 @testset "RootSearch interface" begin
     contractor = Newton(sin, cos)
-    search = RootSearch(-10..10, contractor, 1e-3)
+    search = RootSearch(-10..10, contractor, BreadthFirstSearch, 1e-3)
     state = start(search)
 
     # check that original and copy are independent