Explicitely give derivative to Newton contractor. (#51)

* Explicitely give derivative to Newton contractor.

* Refactor and  methods.

* Refactor roots and branch_and_prune methods.

* Fix notation and simplify Newton type

Author: Kolaru

src/branch_and_prune.jl

@@ -5,42 +5,89 @@ export branch_and_prune, Bisection, Newton
 
 diam(x::Root) = diam(x.interval)
 
-
 Base.size(x::Interval) = (1,)
 
 isinterior{N}(X::IntervalBox{N}, Y::IntervalBox{N}) = all(isinterior.(X, Y))
 
+# contractors:
+"""
+    Contractor{F}
 
+    Abstract type for contractors.
 """
-    branch_and_prune(f, X, contractor, tol=1e-3)
+abstract type Contractor{F} end
 
-Generic branch and prune routine for finding isolated roots of a function
-`f:R^n → R^n` in a box.
+export Bisection, Newton
 
-Inputs:
-- `f`: function whose roots will be found
-- `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`.
+struct Bisection{F} <: Contractor{F}
+    f::F
+end
 
-"""
-function branch_and_prune(f, X, contractor, tol=1e-3)
+function (contractor::Bisection)(X)
+    image = contractor.f(X)
 
-    input_dim = length(X)
-    output_dim = length(X)
+    if !(contains_zero(image))
+        return :empty, X
+    end
+
+    return :unknown, X
+end
 
-    # @show input_dim
-    # @show output_dim
 
-    # if !(input_dim == output_dim)
-    #     throw(ArgumentError("Input dimension ($input_dim) and output dimension ($output_dim) must be the same."))
-    # end
+struct Newton{F,FP,O} <: Contractor{F}
+    f::F
+    f_prime::FP
+    op::O
+end
 
-    contract = contractor(Val{input_dim}, f)
+function Newton(f::Function, f_prime::Function)
+    Newton(f, f_prime, N)
+end
 
-    # main algorithm:
+function Newton(f_prime::Function)
+    return NewtonConstructor(f_prime)
+end
 
+function (C::Newton)(X)
+    # use Bisection contractor for this:
+    if !(contains_zero(IntervalBox(C.f(X))))
+        return :empty, X
+    end
+
+    # given that have the Jacobian, can also do mean value form
+
+    NX = C.op(C.f, C.f_prime, X) ∩ X
+
+    isempty(NX) && return :empty, X
+
+    if NX ⪽ X  # isinterior; know there's a unique root inside
+        NX =  refine(X -> C.op(C.f, C.f_prime, X), NX)
+        return :unique, NX
+    end
+
+    return :unknown, NX
+end
+
+
+struct NewtonConstructor{FP}
+    f_prime::FP
+end
+
+
+"""
+    branch_and_prune(X, contract, tol=1e-3)
+
+Generic branch and prune routine for finding isolated roots using the `contract`
+function as the contractor.
+
+Inputs:
+- `X`: `Interval` or `IntervalBox`
+- `contractor`: 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`
+
+"""
+function branch_and_prune(X, contractor, tol=1e-3)
     working = [X]
     outputs = Root{typeof(X)}[]
 
@@ -51,7 +98,7 @@ function branch_and_prune(f, X, contractor, tol=1e-3)
         # @show working
         X = pop!(working)
 
-        status, output = contract(X)
+        status, output = contractor(X)
 
         if status == :empty
             continue
@@ -67,23 +114,11 @@ function branch_and_prune(f, X, contractor, tol=1e-3)
 
             push!(working, X1, X2)
         end
-
     end
 
     return outputs
 end
 
-branch_and_prune(f, X::Root, contractor, tol=1e-3) =
-    branch_and_prune(f, X.interval, contractor, tol)
-
-
-function branch_and_prune(f, V::Vector{Root{T}}, contractor, tol=1e-3) where {T}
-    reduce(append!, Root{T}[], [branch_and_prune(f, X.interval, contractor, tol) for X in V])
-end
-
-function branch_and_prune(f, V::Vector{T}, contractor, tol=1e-3) where {T}
-    reduce(append!, Root{T}[], [branch_and_prune(f, X, contractor, tol) for X in V])
-end
 
 export recursively_branch_and_prune
 
@@ -99,54 +134,12 @@ function recursively_branch_and_prune(h, X, contractor=BisectionContractor, fina
     return roots
 end
 
-"""
-If the input interval is complex, treat `f` as a complex function, currently of one complex variable `z`.
-"""
-function branch_and_prune{T}(f, Xc::Complex{Interval{T}}, contractor, tol=1e-3)
-
-    g = realify(f)
-    Y = IntervalBox(reim(Xc))
-
-    roots = branch_and_prune(g, Y, contractor, tol)
-
-    # @show roots
-
-    return [Root(Complex(root.interval...), root.status) for root in roots]
-end
-
-
 
 contains_zero{T}(X::Interval{T}) = zero(T) ∈ X
 contains_zero(X::SVector) = all(contains_zero.(X))
 contains_zero(X::IntervalBox) = all(contains_zero(X[i]) for i in 1:length(X))
 
 
-# contractors:
-
-abstract type Contractor{F} end
-
-export Bisection, Newton
-
-struct Bisection{F} <: Contractor{F}
-    dimension::Int
-    f::F
-end
-
-Bisection{n}(::Type{Val{n}}, f) = Bisection(n, f)
-
-
-function (contractor::Bisection)(X)
-    image = contractor.f(X)
-
-    if !(contains_zero(image))
-        return :empty, X
-    end
-
-    return :unknown, X
-end
-
-
-
 """
 Generic refine operation for Krawczyk and Newton.
 This assumes that it is already known that `X` contains a unique root.
@@ -166,49 +159,7 @@ function refine(op, X)
 end
 
 
-
-struct Newton{F,FP,O} <: Contractor{F}
-    dimension::Int
-    f::F
-    fp::FP
-    op::O
-end
-
-function Newton(::Type{Val{1}}, f::Function)
-    f_prime = x -> ForwardDiff.derivative(f, x)
-    Newton(1, f, f_prime, N)
-end
-
-function Newton{n}(::Type{Val{n}}, f::Function)
-    f_prime = x -> ForwardDiff.jacobian(f, x)
-    Newton(n, f, f_prime, N)
-end
-
-
-function (C::Newton)(X)
-
-    # use Bisection contractor for this:
-    if !(contains_zero(IntervalBox(C.f(X))))
-        return :empty, X
-    end
-
-    # given that have the Jacobian, can also do mean value form
-
-
-    NX = C.op(C.f, C.fp, X) ∩ X
-
-    isempty(NX) && return :empty, X
-
-
-    if NX ⪽ X  # isinterior; know there's a unique root inside
-        NX =  refine(X -> C.op(C.f, C.fp, X), NX)
-        return :unique, NX
-    end
-
-
-    return :unknown, NX
-end
-
+IntervalLike{T} = Union{Interval{T}, IntervalBox{T}}
 
 """
     roots(f, X, contractor, tol=1e-3)
@@ -221,9 +172,52 @@ 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
+    `Newton(f_prime)` where `f_prime` is the derivative or jacobian of `f`.
 
 """
-roots{C<:Contractor}(f, X, contractor::Type{C}, tol::Float64=1e-3) = branch_and_prune(f, X, contractor, tol)
+# 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)
+end
+
+function roots(f, X::Interval{T}, ::Type{Newton}, tol::Float64=1e-3) where {T}
+    branch_and_prune(X, Newton(f, x -> ForwardDiff.derivative(f, x)), tol)
+end
+
+function roots(f, X::IntervalBox{T}, ::Type{Newton}, tol::Float64=1e-3) where {T}
+    branch_and_prune(X, Newton(f, x -> ForwardDiff.jacobian(f, x)), tol)
+end
+
+function roots{T}(f, X::IntervalLike{T}, nc::NewtonConstructor, tol::Float64=1e-3)
+    branch_and_prune(X, Newton(f, nc.f_prime), tol)
+end
+
+# `roots` function for cases where `X` is not an `Interval` or `IntervalBox`
+function roots(f, V::Vector{Root{T}}, contractor::Type{C}, tol::Float64=1e-3) where {T, C<:Contractor}
+    reduce(append!, Root{T}[], [roots(f, X.interval, contractor, tol) for X in V])
+end
+
+function roots(f, V::Vector{T}, contractor::Type{C}, tol::Float64=1e-3) where {T, C<:Contractor}
+    reduce(append!, Root{T}[], [roots(f, X, contractor, tol) for X in V])
+end
+
+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}}, nc::NewtonConstructor, tol::Float64=1e-3) where {T}
+    g = realify(f)
+    g_prime = realify_derivative(nc.f_prime)
+    Y = IntervalBox(reim(Xc))
+    rts = roots(g, Y, Newton(g_prime), tol)
+
+    return [Root(Complex(root.interval...), root.status) for root in rts]
+end
 
-roots(f, X, tol::Float64=1e-3) = branch_and_prune(f, X, Newton, tol)
+# Default
+roots(f, X, tol::Float64=1e-3) = roots(f, X, Newton, tol)

src/complex.jl

@@ -38,3 +38,16 @@ function realify(f)
     return g
 
 end
+
+"""
+Takes the derivative of a complex function and returns the real jacobian
+that implements it.
+"""
+function realify_derivative(fp)
+    function g_jac(x)
+        z = Compl(x[1], x[2])
+        fpz = fp(z)
+        SMatrix{2, 2}(fpz.re, fpz.im, -fpz.im, fpz.re)
+    end
+    return g_jac
+end

test/roots.jl

@@ -11,18 +11,25 @@ using Base.Test
     @test length(rts) == 3
     rts = roots(sin, rts, Newton)
     @test all(root.status == :unique for root in rts)
+
+    rts = roots(sin, -5..5, Newton)
+    @test rts == roots(sin, -5..5, Newton(x -> ForwardDiff.derivative(sin, x)))
 end
 
 
 @testset "2D roots" begin
     f(x, y) = SVector(x^2 + y^2 - 1, y - 2x)
     f(X) = f(X...)
+    X = (-6..6) × (-6..6)
 
-    rts = roots(f, (-6..6) × (-6..6), Bisection, 1e-3)
+    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(xx -> ForwardDiff.jacobian(f, xx)))
 end
 
 @testset "Complex roots" begin
@@ -36,6 +43,16 @@ end
     @test length(rts) == 3
     rts = roots(f, Xc)
     @test length(rts) == 3
+
+    rts = roots(f, Xc, Newton)
+    rts2 = roots(f, Xc, Newton(z -> 3*z^2))
+    intervals = [reim(rt.interval) for rt in rts]
+    intervals2 = [reim(rt.interval) for rt in rts2]
+    d = []
+    for (I, I2) in zip(sort(intervals), sort(intervals2))
+        append!(d, abs.(I .- I2))
+    end
+    @test all(d .< 1e-15)
 end
 
 # From R docs:
@@ -51,12 +68,19 @@ end
     end
 
     X = (-5..5)
-    rts = roots(g, IntervalBox(X, 3))
+    XX = IntervalBox(X, 3)
+    rts = roots(g, XX, Newton)
     @test length(rts) == 4
+    @test rts == roots(g, XX, Newton(xx -> ForwardDiff.jacobian(g, xx)))
 end
 
 @testset "Stationary points" begin
     f(xx) = ( (x, y) = xx; sin(x) * sin(y) )
-    rts = roots(∇(f), IntervalBox(-5..6, 2), Newton, 1e-5)
+    gradf = ∇(f)
+    XX = IntervalBox(-5..6, 2)
+    tol = 1e-5
+
+    rts = roots(gradf, XX, Newton, tol)
     @test length(rts) == 25
+    @test rts == roots(gradf, XX, Newton(xx -> ForwardDiff.jacobian(gradf, xx)), tol)
 end