Refactor contractors (#188)

* Refactor contractors

* Rename Contractor to AbstractContractor

docs/src/internals.md

@@ -39,7 +39,7 @@ Newton{typeof(sin),typeof(cos)}(sin, cos)
 julia> C(Root(pi ± 0.001, :unknown), 1e-10)
 Root([3.14159, 3.1416], :unique)
 
-julia> C(Root(2 ± 0.001, :unkown), 1e-10)
+julia> C(Root(2 ± 0.001, :unknown), 1e-10)
 Root([1.99899, 2.00101], :empty)
 ```
 

src/contractors.jl

@@ -1,126 +1,99 @@
-export Contractor
 export Bisection, Newton, Krawczyk
 
 """
-    𝒩(f, f′, X, α)
+    AbstractContractor{F}
 
-Single-variable Newton operator.
-
-The symbol for the operator is accessed with `\\scrN<tab>`.
-"""
-function 𝒩(f, f′, X::Interval{T}, α) where {T}
-    m = Interval(mid(X, α))
-
-    return m - (f(m) / f′(X))
-end
-
-"""
-    𝒩(f, jacobian, X, α)
-
-Multi-variable Newton operator.
+Abstract type for contractors.
 """
-function 𝒩(f::Function, jacobian::Function, X::IntervalBox, α)  # multidimensional Newton operator
-    m = Interval.(mid(X, α))
-    J = jacobian(X)
-    y = gauss_elimination_interval(J, f(m)) # J \ f(m)
-    return IntervalBox(m .- y)
-end
+abstract type AbstractContractor{F} end
 
 """
-    𝒦(f, f′, X, α)
-
-Single-variable Krawczyk operator.
+    Bisection{F} <: AbstractContractor{F}
 
-The symbol for the operator is accessed with `\\scrK<tab>`.
+AbstractContractor type for the bisection method.
 """
-function 𝒦(f, f′, X::Interval{T}, α) where {T}
-    m = Interval(mid(X, α))
-    Y = 1 / f′(m)
-
-    return m - Y*f(m) + (1 - Y*f′(X)) * (X - m)
+struct Bisection{F} <: AbstractContractor{F}
+    f::F
 end
 
 """
-    𝒦(f, jacobian, X, α)
-
-Multi-variable Krawczyk operator.
-"""
-function 𝒦(f, jacobian, X::IntervalBox{T}, α) where {T}
-    m = mid(X, α)
-    mm = IntervalBox(m)
-    J = jacobian(X)
-    Y = mid.(inv(jacobian(mm)))
+    Newton{F, FP} <: AbstractContractor{F}
 
-    return m - Y*f(mm) + (I - Y*J) * (X.v - m)
-end
+AbstractContractor type for the interval Newton method.
 
+# Fields
+    - `f::F`: function whose roots are searched
+    - `f::FP`: derivative or jacobian of `f`
 
-"""
-    Contractor{F}
+-----
 
-Abstract type for contractors.
-"""
-abstract type Contractor{F} end
+    (N::Newton)(X, α=where_bisect)
 
-"""
-    Bisection{F} <: Contractor{F}
+Contract an interval `X` using Newton operator and return the
+contracted interval together with its status.
 
-Contractor type for the bisection method.
+# Inputs
+    - `R`: Root object containing the interval to contract.
+    - `α`: Point of bisection of intervals.
 """
-struct Bisection{F} <: Contractor{F}
+struct Newton{F, FP} <: AbstractContractor{F}
     f::F
+    f′::FP   # use \prime<TAB> for ′
 end
 
-function (contractor::Bisection)(r, tol)
-    X = interval(r)
-    image = (contractor.f)(X)
-
-    if root_status(r) == :empty || !(contains_zero(image))
-        return Root(X, :empty)
-    end
-
-    return Root(X, :unknown)
+function (N::Newton)(X::Interval ; α=where_bisect)
+    m = Interval(mid(X, α))
+    return m - (N.f(m) / N.f′(X))
 end
 
-for (Method, 𝒪) in ((:Newton, 𝒩), (:Krawczyk, 𝒦))
-    doc = """
-        $Method{F, FP} <: Contractor{F}
+function (N::Newton)(X::IntervalBox ; α=where_bisect)
+    m = Interval.(mid(X, α))
+    J = N.f′(X)
+    y = gauss_elimination_interval(J, N.f(m))  # J \ f(m)
+    return IntervalBox(m .- y)
+end
 
-    Contractor type for the $Method method.
+"""
+    Krawczyk{F, FP} <: AbstractContractor{F}
 
-    # Fields
-      - `f::F`: function whose roots are searched
-      - `f::FP`: derivative or jacobian of `f`
+AbstractContractor type for the interval Krawczyk method.
 
-    -----
+# Fields
+    - `f::F`: function whose roots are searched
+    - `f::FP`: derivative or jacobian of `f`
 
-        (C::$Method)(X, tol, α=where_bisect)
+-----
 
-    Contract an interval `X` using $Method operator and return the
-    contracted interval together with its status.
+    (K::Krawczyk)(X ; α=where_bisect)
 
-    # Inputs
-      - `R`: Root object containing the interval to contract.
-      - `tol`: Precision to which unique solutions are refined.
-      - `α`: Point of bisection of intervals.
-    """
+Contract an interval `X` using Krawczyk operator and return the
+contracted interval together with its status.
 
-    @eval begin
-        struct $Method{F, FP} <: Contractor{F}
-            f::F
-            f′::FP   # use \prime<TAB> for ′
-        end
+# Inputs
+    - `R`: Root object containing the interval to contract.
+    - `α`: Point of bisection of intervals.
+"""
+struct Krawczyk{F, FP} <: AbstractContractor{F}
+    f::F
+    f′::FP   # use \prime<TAB> for ′
+end
 
-        function (C::$Method)(R, tol, α=where_bisect)
-            op = x -> $𝒪(C.f, C.f′, x, α)
-            rt = determine_region_status(op, C.f, R)
-            return refine(op, rt, tol)
-        end
+function (K::Krawczyk)(X::Interval ; α=where_bisect)
+    m = Interval(mid(X, α))
+    Y = 1 / K.f′(m)
 
-        @doc $doc $Method
-    end
+    return m - Y*K.f(m) + (1 - Y*K.f′(X)) * (X - m)
 end
 
+function (K::Krawczyk)(X::IntervalBox ; α=where_bisect)
+    jacobian = K.f′
+    m = mid(X, α)
+    mm = IntervalBox(m)
+    J = jacobian(X)
+    Y = mid.(inv(jacobian(mm)))
+
+    return m - Y*K.f(mm) + (I - Y*J) * (X.v - m)
+end
 
 """
     safe_isempty(X)
@@ -131,52 +104,57 @@ of `Interval`.
 safe_isempty(X) = isempty(IntervalBox(X))
 
 """
-    determine_region_status(contract, f, R)
+    contract(contractor, R)
 
-Contraction operation for contractors using the first derivative of the
-function.
-
-Currently `Newton` and `Krawczyk` contractors use this.
+Contract the region R using the given contractor.
 """
-function determine_region_status(op, f, R)
+function contract(B::Bisection, R::Root)
     X = interval(R)
-    former_status = root_status(R)
+    R.status == :empty && return Root(X, :empty)
 
-    imX = f(X)
+    imX = B.f(X)
 
-    if former_status == :empty || !(contains_zero(imX))
+    if !(contains_zero(imX)) || safe_isempty(imX)
         return Root(X, :empty)
     end
 
-    safe_isempty(imX) && return Root(X, :empty)  # X is fully outside of the domain of f
+    return Root(X, :unknown)
+end
+
+function contract(C::Union{Newton, Krawczyk}, R::Root)
+    # We first check with the simple bisection method
+    # If we can prove it is empty at this point, we don't go further
+    R2 = contract(Bisection(C.f), R)
+    R2.status == :empty && return R2
 
-    contracted_X = op(X)
+    X = interval(R)
+    contracted_X = C(X)
 
     # Only happens if X is partially out of the domain of f
     safe_isempty(contracted_X) && return Root(X, :unknown)  # force bisection
 
-    # given that have the Jacobian, can also do mean value form
     NX = contracted_X ∩ X
 
     isinf(X) && return Root(NX, :unknown)  # force bisection
     safe_isempty(NX) && return Root(X, :empty)
 
-    if former_status == :unique || NX ⪽ X  # isinterior; know there's a unique root inside
+    if R.status == :unique || NX ⪽ X  # isinterior, we know there's a unique root inside
         return Root(NX, :unique)
     end
 
     return Root(NX, :unknown)
 end
 
 """
-    refine(op, X::Region, tol)
+    refine(C, X::Region, tol)
+
+Refine a interval known to contain a solution.
 
-Generic refine operation for Krawczyk and Newton.
 This function assumes that it is already known that `X` contains a unique root.
 """
-function refine(op, X::Region, tol)
-    while diam(X) > tol  # avoid problem with tiny floating-point numbers if 0 is a root
-        NX = op(X) ∩ X
+function refine(C::Union{Newton, Krawczyk}, X::Region, root_problem)
+    while diam(X) > root_problem.abstol
+        NX = C(X) ∩ X
         NX == X && break  # reached limit of precision
         X = NX
     end
@@ -190,7 +168,7 @@ end
 Wrap the refine method to leave unchanged intervals that are not guaranteed to
 contain an unique solution.
 """
-function refine(op, R::Root, tol)
+function refine(C::AbstractContractor, R::Root, root_problem)
     root_status(R) != :unique && return R
-    return Root(refine(op, interval(R), tol), :unique)
+    return Root(refine(C, interval(R), root_problem), :unique)
 end

src/root_object.jl

@@ -12,7 +12,7 @@ by the `roots` function.
 # Fields
   - `interval`: a region (either `Interval` of `IntervalBox`) searched for
         roots.
-  - `status`: the status of the region, valid values are `:empty`, `unkown` and
+  - `status`: the status of the region, valid values are `:empty`, `unknown` and
         `:unique`.
 """
 struct Root{T}

src/roots.jl

@@ -10,23 +10,25 @@ isnan(r::Root) = isnan(interval(r))
 struct RootProblem{T}
     abstol::T
 end
-
+   
 function bisect(r::Root)
     Y1, Y2 = bisect(interval(r))
     return Root(Y1, :unknown), Root(Y2, :unknown)
 end
 
 function process(contractor, root_problem, r::Root)
-    contracted_root = contractor(r, root_problem.abstol)
-    status = root_status(contracted_root)
+    contracted_root = contract(contractor, r)
+    refined_root = refine(contractor, contracted_root, root_problem)
+
+    status = root_status(refined_root)
 
-    status == :unique && return :store, contracted_root
-    status == :empty && return :prune, contracted_root
+    status == :unique && return :store, refined_root
+    status == :empty && return :prune, refined_root
 
     if status == :unknown
         # Avoid infinite division of intervals with singularity
-        isnan(contracted_root) && diam(r) < root_problem.abstol && return :store, r
-        diam(contracted_root) < root_problem.abstol && return :store, contracted_root
+        isnan(refined_root) && diam(r) < root_problem.abstol && return :store, r
+        diam(refined_root) < root_problem.abstol && return :store, refined_root
 
         return :branch, r
     else
@@ -84,26 +86,26 @@ Inputs:
 """
 function roots(f::Function, X, contractor::Type{C}=default_contractor,
                search_order::Type{S}=default_search_order,
-               tol::Float64=default_tolerance) where {C <: Contractor, S <: SearchOrder}
+               tol::Float64=default_tolerance) where {C <: AbstractContractor, S <: SearchOrder}
 
     _roots(f, X, contractor, search_order, tol)
 end
 
 function roots(f::Function, deriv::Function, X, contractor::Type{C}=default_contractor,
                search_order::Type{S}=default_search_order,
-               tol::Float64=default_tolerance) where {C <: Contractor, S <: SearchOrder}
+               tol::Float64=default_tolerance) where {C <: AbstractContractor, S <: SearchOrder}
 
     _roots(f, deriv, X, contractor, search_order, tol)
 end
 
 function roots(f::Function, X, contractor::Type{C},
-               tol::Float64) where {C <: Contractor}
+               tol::Float64) where {C <: AbstractContractor}
 
     _roots(f, X, contractor, default_search_order, tol)
 end
 
 function roots(f::Function, deriv::Function, X, contractor::Type{C},
-               tol::Float64) where {C <: Contractor}
+               tol::Float64) where {C <: AbstractContractor}
 
     _roots(f, deriv, X, contractor, default_search_order, tol)
 end
@@ -153,35 +155,35 @@ end
 
 # Acting on `Interval`
 function _roots(f, X::Region, contractor::Type{C},
-               search_order::Type{S}, tol::Float64) where {C <: Contractor, S <: SearchOrder}
+               search_order::Type{S}, tol::Float64) where {C <: AbstractContractor, S <: SearchOrder}
 
     _roots(f, Root(X, :unknown), contractor, search_order, tol)
 end
 
 function _roots(f, deriv, X::Region, contractor::Type{C},
-               search_order::Type{S}, tol::Float64) where {C <: Contractor, S <: SearchOrder}
+               search_order::Type{S}, tol::Float64) where {C <: AbstractContractor, S <: SearchOrder}
 
     _roots(f, deriv, Root(X, :unknown), contractor, search_order, tol)
 end
 
 
 # Acting on `Vector` of `Root`
 function _roots(f, V::Vector{Root{T}}, contractor::Type{C},
-               search_order::Type{S}, tol::Float64) where {T, C <: Contractor, S <: SearchOrder}
+               search_order::Type{S}, tol::Float64) where {T, C <: AbstractContractor, S <: SearchOrder}
 
     vcat(_roots.(f, V, contractor, search_order, tol)...)
 end
 
 function _roots(f, deriv, V::Vector{Root{T}}, contractor::Type{C},
-               search_order::Type{S}, tol::Float64) where {T, C <: Contractor, S <: SearchOrder}
+               search_order::Type{S}, tol::Float64) where {T, C <: AbstractContractor, S <: SearchOrder}
 
     vcat(_roots.(f, deriv, V, contractor, search_order, tol)...)
 end
 
 
 # Acting on complex `Interval`
 function _roots(f, Xc::Complex{Interval{T}}, contractor::Type{C},
-               search_order::Type{S}, tol::Float64) where {T, C <: Contractor, S <: SearchOrder}
+               search_order::Type{S}, tol::Float64) where {T, C <: AbstractContractor, S <: SearchOrder}
 
     g = realify(f)
     Y = IntervalBox(reim(Xc)...)