Implement Krawczyk method (#57)

* Krawczyk contractor implementation

* Add tests and fix complex roots for Krawczyk

* Fix requested changes

- Uniformized tests
- Uniformized contraction for NewtonLike contractors
- Add a bit of documentation

Author: Kolaru

src/contractors.jl

@@ -1,14 +1,20 @@
-# contractors:
-"""
+export Bisection, Newton, Krawczyk
+
+Base.isinf(X::IntervalBox) = any(isinf.(X))
+IntervalArithmetic.mid(X::IntervalBox, α) = mid.(X, α)
+
+doc"""
     Contractor{F}
 
     Abstract type for contractors.
 """
 abstract type Contractor{F} end
 
-export Bisection, Newton
+doc"""
+    Bisection{F} <: Contractor{F}
 
-# bisection
+    Contractor type for the bisection method.
+"""
 struct Bisection{F} <: Contractor{F}
     f::F
 end
@@ -23,38 +29,54 @@ function (contractor::Bisection)(X, tol)
     return :unknown, X
 end
 
-# Newton
-struct Newton{F,FP} <: Contractor{F}
-    f::F
-    f′::FP   # use \prime<TAB> for ′
-end
+doc"""
+    newtonlike_contract(op, X, tol)
 
-Base.isinf(X::IntervalBox) = any(isinf.(X))
+    Contraction operation for contractors using the first derivative of the
+    function. This contraction use a bisection scheme to refine the intervals
+    with `:unkown` status.
 
-function (C::Newton)(X, tol)
+    Currently `Newton` and `Krawczyk` contractors uses this.
+"""
+function newtonlike_contract(op, C, X, tol)
     # use Bisection contractor for this:
     if !(contains_zero(C.f(X)))
         return :empty, X
     end
 
     # given that have the Jacobian, can also do mean value form
 
-    NX = 𝒩(C.f, C.f′, X) ∩ X
+    NX = op(C.f, C.f′, X) ∩ X
 
     isempty(NX) && return :empty, X
-
-    if isinf(X)
-        return :unknown, NX  # force bisection
-    end
+    isinf(X) && return :unknown, NX  # force bisection
 
     if NX ⪽ X  # isinterior; know there's a unique root inside
-        NX =  refine(X -> 𝒩(C.f, C.f′, X), NX, tol)
+        NX =  refine(X -> op(C.f, C.f′, X), NX, tol)
         return :unique, NX
     end
 
     return :unknown, NX
 end
 
+doc"""
+    Newton{F, FP} <: Contractor{F}
+
+    Contractor type for the Newton method.
+
+    # Fields
+        - `f::F`: function whose roots are searched
+        - `f::FP`: derivative or jacobian of `f`
+"""
+struct Newton{F,FP} <: Contractor{F}
+    f::F
+    f′::FP   # use \prime<TAB> for ′
+end
+
+function (C::Newton)(X, tol)
+    newtonlike_contract(𝒩, C, X, tol)
+end
+
 
 doc"""
 Single-variable Newton operator
@@ -77,20 +99,57 @@ function 𝒩{T}(f, X::Interval{T}, dX::Interval{T})
     m - (f(m) / dX)
 end
 
-IntervalArithmetic.mid(X::IntervalBox, α) = mid.(X, α)
-
 doc"""
 Multi-variable Newton operator.
 """
 function 𝒩(f::Function, jacobian::Function, X::IntervalBox)  # multidimensional Newton operator
-
     m = IntervalBox(Interval.(mid(X, where_bisect)))
     J = jacobian(SVector(X))
 
     return IntervalBox(m - (J \ f(m)))
 end
 
 
+doc"""
+    Krawczyk{F, FP} <: Contractor{F}
+
+    Contractor type for the Krawczyk method.
+
+    # Fields
+        - `f::F`: function whose roots are searched
+        - `f::FP`: derivative or jacobian of `f`
+"""
+struct Krawczyk{F, FP} <: Contractor{F}
+    f::F
+    f′::FP   # use \prime<TAB> for ′
+end
+
+function (C::Krawczyk)(X, tol)
+    newtonlike_contract(𝒦, C, X, tol)
+end
+
+
+doc"""
+Single-variable Krawczyk operator
+"""
+function 𝒦(f, f′, X::Interval{T}) where {T}
+    m = Interval(mid(X))
+    Y = 1/f′(m)
+
+    m - Y*f(m) + (1 - Y*f′(X))*(X - m)
+end
+
+doc"""
+Multi-variable Krawczyk operator
+"""
+function 𝒦(f, jacobian, X::IntervalBox{T}) where {T}
+    m = mid(X)
+    J = jacobian(X)
+    Y = inv(jacobian(m))
+    m = IntervalBox(Interval.(m))
+
+    IntervalBox(m - Y*f(m) + (I - Y*J)*(X - m))
+end
 
 """
 Generic refine operation for Krawczyk and Newton.

src/roots.jl

@@ -82,7 +82,7 @@ function as the contractor.
 
 Inputs:
 - `X`: `Interval` or `IntervalBox`
-- `contractor`: function that determines the status of a given box `X`. It
+- `contract`: function that determines the status of a given box `X`. It
     returns the new box and a symbol indicating the status. Current possible
     values are of type `Bisection` or `Newton`
 
@@ -116,9 +116,10 @@ contains_zero(X::IntervalBox) = all(contains_zero.(X))
 
 
 IntervalLike{T} = Union{Interval{T}, IntervalBox{T}}
+NewtonLike = Union{Type{Newton}, Type{Krawczyk}}
 
 """
-    roots(f, X, contractor, tol=1e-3)
+    roots(f, X, contractor, tol=1e-3 ; deriv=nothing)
 
 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.
@@ -128,53 +129,57 @@ Inputs:
 - `X`: `Interval` or `IntervalBox`
 - `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` and `Newton`.
+    the status. Current possible values are `Bisection`, `Newton` and `Krawczyk`
+- `deriv` ; explicit derivative of `f` for `Newton` and `Krawczyk`
 
 """
 # Contractor specific `roots` functions
 function roots(f, X::IntervalLike{T}, ::Type{Bisection}, tol::Float64=1e-3; deriv = nothing) where {T}
     branch_and_prune(X, Bisection(f), tol)
 end
 
-function roots(f, X::Interval{T}, ::Type{Newton}, tol::Float64=1e-3; deriv = nothing) where {T}
+function roots(f, X::Interval{T}, C::NewtonLike, tol::Float64=1e-3; deriv = nothing) where {T}
 
     if deriv == nothing
         deriv = x -> ForwardDiff.derivative(f, x)
     end
 
-    branch_and_prune(X, Newton(f, deriv), tol)
+    branch_and_prune(X, C(f, deriv), tol)
 end
 
-function roots(f, X::IntervalBox{T}, ::Type{Newton}, tol::Float64=1e-3; deriv = nothing) where {T}
+function roots(f, X::IntervalBox{T}, C::NewtonLike, tol::Float64=1e-3; deriv = nothing) where {T}
 
     if deriv == nothing
         deriv = x -> ForwardDiff.jacobian(f, x)
     end
 
-    branch_and_prune(X, Newton(f, deriv), tol)
+    branch_and_prune(X, C(f, deriv), tol)
 end
 
-
 roots(f, r::Root, contractor::Type{C}, tol::Float64=1e-3; deriv = nothing) where {C<:Contractor}  = roots(f, r.interval, contractor, tol; deriv = deriv)
 
+
 # Acting on a Vector:
 
 # TODO: Use previous status information about roots:
-roots(f, V::Vector{Root{T}}, contractor::Type{C}, tol::Float64=1e-3; deriv = nothing) where {T, C<:Contractor} = vcat(roots.(f, V, contractor, tol; deriv = deriv)...)
+roots(f, V::Vector{Root{T}}, contractor::Type{C}, tol::Float64=1e-3;
+    deriv = nothing) where {T, C<:Contractor} = vcat(roots.(f, V, contractor, tol; deriv = deriv)...)
 
 
 
 # Complex:
 
-function roots(f, Xc::Complex{Interval{T}}, contractor::Type{C}, tol::Float64=1e-3) where {T, C<:Contractor}
+function roots(f, Xc::Complex{Interval{T}}, contractor::Type{C},
+        tol::Float64=1e-3) where {T, C<:Contractor}
+
     g = realify(f)
     Y = IntervalBox(reim(Xc))
     rts = roots(g, Y, contractor, tol)
 
     return [Root(Complex(root.interval...), root.status) for root in rts]
 end
 
-function roots(f, Xc::Complex{Interval{T}}, ::Type{Newton}, tol::Float64=1e-15;
+function roots(f, Xc::Complex{Interval{T}}, C::NewtonLike, tol::Float64=1e-3;
         deriv = nothing) where {T}
 
     g = realify(f)
@@ -186,7 +191,7 @@ function roots(f, Xc::Complex{Interval{T}}, ::Type{Newton}, tol::Float64=1e-15;
     end
 
     Y = IntervalBox(reim(Xc))
-    rts = roots(g, Y, Newton, tol; deriv=g_prime)
+    rts = roots(g, Y, C, tol; deriv=g_prime)
 
     return [Root(Complex(root.interval...), root.status) for root in rts]
 end

test/roots.jl

@@ -2,19 +2,38 @@
 using IntervalArithmetic, IntervalRootFinding, StaticArrays
 using Base.Test
 
+function all_unique(rts)
+    all(root_status.(rts) .== :unique)
+end
+
+function test_newtonlike(f, X, method, nsol, tol=1e-3; deriv=nothing)
+    rts = roots(f, X, method)
+    @test length(rts) == nsol
+    @test all_unique(rts)
+    @test rts == roots(f, X, method; deriv = deriv)
+end
+
+newtonlike_methods = [Newton, Krawczyk]
+
 @testset "1D roots" begin
+    # Default
     rts = roots(sin, -5..5)
     @test length(rts) == 3
-    @test length(find(x->x==:unique, [root.status for root in rts])) == 3
+    @test all_unique(rts)
 
+    # Bisection
     rts = roots(sin, -5..6, Bisection)
     @test length(rts) == 3
+
+    # Refinement
     rts = roots(sin, rts, Newton)
-    @test all(root.status == :unique for root in rts)
+    @test all_unique(rts)
 
-    rts = roots(sin, -5..5, Newton)
-    @test rts == roots(sin, -5..5, Newton; deriv = cos)
+    for method in newtonlike_methods
+        test_newtonlike(sin, -5..5, method, 3; deriv=cos)
+    end
 
+    # Infinite interval
     rts = roots(x -> x^2 - 2, -∞..∞)
     @test length(rts) == 2
 end
@@ -25,19 +44,18 @@ end
     f(X) = f(X...)
     X = (-6..6) × (-6..6)
 
+    # Bisection
     rts = roots(f, X, Bisection, 1e-3)
     @test length(rts) == 4
 
-    rts = roots(f, rts, Newton)
-    @test length(rts) == 2
-
-    rts = roots(f, X, Newton)
-    @test rts == roots(f, X, Newton; deriv = xx -> ForwardDiff.jacobian(f, xx))
+    for method in newtonlike_methods
+        test_newtonlike(f, X, method, 2; deriv = xx -> ForwardDiff.jacobian(f, xx))
+    end
 
+    # Infinite interval
     X = IntervalBox(-∞..∞, 2)
     rts = roots(f, X, Newton)
     @test length(rts) == 2
-
 end
 
 
@@ -55,9 +73,13 @@ end
 
     X = (-5..5)
     XX = IntervalBox(X, 3)
-    rts = roots(g, XX, Newton)
-    @test length(rts) == 4
-    @test rts == roots(g, XX, Newton; deriv = xx -> ForwardDiff.jacobian(g, xx))
+
+    for method in newtonlike_methods
+        rts = roots(g, XX, method)
+        @test length(rts) == 4
+        @test all_unique(rts)
+        @test rts == roots(g, XX, method; deriv = xx -> ForwardDiff.jacobian(g, xx))
+    end
 end
 
 @testset "Stationary points" begin
@@ -66,34 +88,27 @@ end
     XX = IntervalBox(-5..6, 2)
     tol = 1e-5
 
-    rts = roots(gradf, XX, Newton, tol)
-    @test length(rts) == 25
-    @test rts == roots(gradf, XX, Newton, tol; deriv = xx -> ForwardDiff.jacobian(gradf, xx))
+    for method in newtonlike_methods
+        test_newtonlike(gradf, XX, method, 25, tol; deriv = xx -> ForwardDiff.jacobian(gradf, xx))
+    end
 end
 
 @testset "Complex roots" begin
     X = -5..5
     Xc = Complex(X, X)
     f(z) = z^3 - 1
 
-    rts = roots(f, Xc, Bisection, 1e-3)
-    @test length(rts) == 7
-    rts = roots(f, rts, Newton)
-    @test length(rts) == 3
+    # Default
     rts = roots(f, Xc)
     @test length(rts) == 3
 
-    rts = roots(f, Xc, Newton)
-    rts2 = roots(f, Xc, Newton; deriv = z -> 3*z^2)
-
-    intervals = [reim(rt.interval) for rt in rts]
-    intervals2 = [reim(rt.interval) for rt in rts2]
+    # Bisection
+    rts = roots(f, Xc, Bisection, 1e-3)
+    @test length(rts) == 7
 
-    d = []
-    for (I, I2) in zip(sort(intervals), sort(intervals2))
-        append!(d, abs.(I .- I2))
+    for method in newtonlike_methods
+        test_newtonlike(f, Xc, method, 3; deriv = z -> 3*z^2)
     end
-    @test all(d .< 1e-15)
 end
 
 @testset "RootSearch interface" begin