Preliminary work on complex derivatives complete, everything should at least work.

dextorious · dextorious · commit 4803b91b8147 · 2017-11-08T01:11:24.000+02:00
diff --git a/src/derivatives.jl b/src/derivatives.jl
@@ -42,7 +42,7 @@ function finite_difference!(df::AbstractArray{<:Real}, f, x::AbstractArray{<:Rea
         epsilon_complex = eps(epsilon_elemtype)
         @. df = imag(f(x+im*epsilon_complex)) / epsilon_complex
     else
-        error("Unrecognized fdtype: valid values are Val{:forward}, Val{:central} and Val{:complex}.")
+        fdtype_error(Val{:Real})
     end
     df
 end
@@ -61,7 +61,6 @@ function finite_difference!(df::AbstractArray{<:Number}, f, x::AbstractArray{<:N
         end
     end
     if fdtype == Val{:forward}
-        @show typeof(x)
         epsilon_factor = compute_epsilon_factor(Val{:forward}, eltype(epsilon))
         @. epsilon = compute_epsilon(Val{:forward}, real(x), epsilon_factor)
         if typeof(fx) == Void
@@ -72,8 +71,8 @@ function finite_difference!(df::AbstractArray{<:Number}, f, x::AbstractArray{<:N
         epsilon_factor = compute_epsilon_factor(Val{:central}, eltype(epsilon))
         @. epsilon = compute_epsilon(Val{:central}, real(x), epsilon_factor)
         @. df = real(f(x+epsilon) - f(x-epsilon)) / (2 * epsilon) + im*imag(f(x+im*epsilon) - f(x-epsilon)) / (2 * epsilon)
-    elseif fdtype == Val{:complex}
-        error("Invalid fdtype value, Val{:complex} not implemented for complex-valued functions.")
+    else
+        fdtype_error(Val{:Complex})
     end
     df
 end
@@ -82,50 +81,90 @@ end
 #=
 Optimized implementations for StridedArrays.
 =#
-function finite_difference!(df::StridedArray{<:Real}, f, x::StridedArray{<:Real},
-    ::Type{Val{:central}}, ::Type{Val{:Real}}, ::Type{Val{:Default}},
+# for R -> R^n
+function finite_difference!(df::StridedArray{<:Real}, f, x::Real,
+    fdtype::DataType, ::Type{Val{:Real}}, ::Type{Val{:Default}},
     fx::Union{Void,StridedArray{<:Real}}=nothing, epsilon::Union{Void,StridedArray{<:Real}}=nothing, return_type::DataType=eltype(x))
 
     epsilon_elemtype = compute_epsilon_elemtype(epsilon, x)
-    epsilon_factor = compute_epsilon_factor(Val{:central}, epsilon_elemtype)
-    @inbounds for i in 1 : length(x)
-        epsilon = compute_epsilon(Val{:central}, x[i], epsilon_factor)
-        epsilon_double_inv = one(typeof(epsilon)) / (2*epsilon)
-        x_plus, x_minus = x[i]+epsilon, x[i]-epsilon
-        df[i] = (f(x_plus) - f(x_minus)) * epsilon_double_inv
+    if fdtype == Val{:forward}
+        epsilon = compute_epsilon(Val{:forward}, x)
+        if typeof(fx) == Void
+            df .= (f(x+epsilon) - f(x)) / epsilon
+        else
+            df .= (f(x+epsilon) - fx) / epsilon
+        end
+    elseif fdtype == Val{:central}
+        epsilon = compute_epsilon(Val{:central}, x)
+        df .= (f(x+epsilon) - f(x-epsilon)) / (2*epsilon)
+    elseif fdtype == Val{:complex}
+        epsilon = eps(eltype(x))
+        df .= imag(f(x+im*epsilon)) / epsilon
+    else
+        fdtype_error(Val{:Real})
     end
     df
 end
 
+# for R^n -> R^n
 function finite_difference!(df::StridedArray{<:Real}, f, x::StridedArray{<:Real},
-    ::Type{Val{:forward}}, ::Type{Val{:Real}}, ::Type{Val{:Default}},
+    fdtype::DataType, ::Type{Val{:Real}}, ::Type{Val{:Default}},
     fx::Union{Void,StridedArray{<:Real}}=nothing, epsilon::Union{Void,StridedArray{<:Real}}=nothing, return_type::DataType=eltype(x))
 
     epsilon_elemtype = compute_epsilon_elemtype(epsilon, x)
-    epsilon_factor = compute_epsilon_factor(Val{:forward}, epsilon_elemtype)
-    @inbounds for i in 1 : length(x)
-        epsilon = compute_epsilon(Val{:forward}, x[i], epsilon_factor)
-        x_plus = x[i] + epsilon
-        if typeof(fx) == Void
-            df[i] = (f(x_plus) - f(x[i])) / epsilon
-        else
-            df[i] = (f(x_plus) - fx[i]) / epsilon
+    if fdtype == Val{:forward}
+        epsilon_factor = compute_epsilon_factor(Val{:forward}, epsilon_elemtype)
+        @inbounds for i in 1 : length(x)
+            epsilon = compute_epsilon(Val{:forward}, x[i], epsilon_factor)
+            x_plus = x[i] + epsilon
+            if typeof(fx) == Void
+                df[i] = (f(x_plus) - f(x[i])) / epsilon
+            else
+                df[i] = (f(x_plus) - fx[i]) / epsilon
+            end
+        end
+    elseif fdtype == Val{:central}
+        epsilon_factor = compute_epsilon_factor(Val{:central}, epsilon_elemtype)
+        @inbounds for i in 1 : length(x)
+            epsilon = compute_epsilon(Val{:central}, x[i], epsilon_factor)
+            epsilon_double_inv = one(typeof(epsilon)) / (2*epsilon)
+            x_plus, x_minus = x[i]+epsilon, x[i]-epsilon
+            df[i] = (f(x_plus) - f(x_minus)) * epsilon_double_inv
+        end
+    elseif fdtype == Val{:complex}
+        epsilon_complex = eps(eltype(x))
+        @inbounds for i in 1 : length(x)
+            df[i] = imag(f(x[i]+im*epsilon_complex)) / epsilon_complex
         end
+    else
+        fdtype_error(Val{:Real})
     end
     df
 end
 
-function finite_difference!(df::StridedArray{<:Real}, f, x::StridedArray{<:Real},
-    ::Type{Val{:complex}}, ::Type{Val{:Real}}, ::Type{Val{:Default}},
-    fx::Union{Void,StridedArray{<:Real}}=nothing, epsilon::Union{Void,StridedArray{<:Real}}=nothing, return_type::DataType=eltype(x))
+# C -> C^n
+function finite_difference!(df::StridedArray{<:Number}, f, x::Number,
+    fdtype::DataType, ::Type{Val{:Complex}}, ::Type{Val{:Default}},
+    fx::Union{Void,StridedArray{<:Number}}=nothing, epsilon::Union{Void,StridedArray{<:Real}}=nothing, return_type::DataType=eltype(x))
 
-    epsilon_complex = eps(eltype(x))
-    @inbounds for i in 1 : length(x)
-        df[i] = imag(f(x[i]+im*epsilon_complex)) / epsilon_complex
+    epsilon_elemtype = compute_epsilon_elemtype(epsilon, x)
+    if fdtype == Val{:forward}
+        epsilon = compute_epsilon(Val{:forward}, real(x[i]))
+        if typeof(fx) == Void
+            df .= ( real( f(x+epsilon) - f(x) ) + im*imag( f(x+im*epsilon) - f(x) ) ) / epsilon
+        else
+            df .= ( real( f(x+epsilon) - fx ) + im*imag( f(x+im*epsilon) - fx )) / epsilon
+        end
+    elseif fdtype == Val{:central}
+        epsilon = compute_epsilon(Val{:central}, real(x[i]))
+        df .= (real(f(x+epsilon) - f(x-epsilon)) + im*imag(f(x+im*epsilon) - f(x-im*epsilon))) / (2 * epsilon)
+    else
+        fdtype_error(Val{:Complex})
     end
     df
 end
 
+# C^n -> C^n
 function finite_difference!(df::StridedArray{<:Number}, f, x::StridedArray{<:Number},
     fdtype::DataType, ::Type{Val{:Complex}}, ::Type{Val{:Default}},
     fx::Union{Void,StridedArray{<:Number}}=nothing, epsilon::Union{Void,StridedArray{<:Real}}=nothing, return_type::DataType=eltype(x))
@@ -147,8 +186,8 @@ function finite_difference!(df::StridedArray{<:Number}, f, x::StridedArray{<:Num
             epsilon = compute_epsilon(Val{:central}, real(x[i]), epsilon_factor)
             df[i] = (real(f(x[i]+epsilon) - f(x[i]-epsilon)) + im*imag(f(x[i]+im*epsilon) - f(x[i]-im*epsilon))) / (2 * epsilon)
         end
-    elseif fdtype == Val{:complex}
-        error("Invalid fdtype value, Val{:complex} not implemented for complex-valued functions.")
+    else
+        fdtype_error(Val{:Complex})
     end
     df
 end
@@ -169,6 +208,8 @@ function finite_difference(f, x::T, fdtype::DataType, funtype::DataType=Val{:Rea
     elseif funtype == Val{:Complex}
         epsilon = compute_epsilon(fdtype, real(x))
         return finite_difference_kernel(f, x, fdtype, funtype, epsilon, f_x)
+    else
+        fdtype_error(funtype)
     end
 end
 
@@ -186,7 +227,7 @@ end
 
 @inline function finite_difference_kernel(f, x::Number, ::Type{Val{:forward}}, ::Type{Val{:Complex}}, epsilon::Real, fx::Union{Void,<:Number}=nothing)
     if typeof(fx) == Void
-        return real((f(x[i]+epsilon) - f(x[i]))) / epsilon + im*imag((f(x[i]+im*epsilon) - fx[i])) / epsilon
+        return real((f(x[i]+epsilon) - f(x[i]))) / epsilon + im*imag((f(x[i]+im*epsilon) - f(x[i]))) / epsilon
     else
         return real((f(x[i]+epsilon) - fx[i])) / epsilon + im*imag((f(x[i]+im*epsilon) - fx[i])) / epsilon
     end
diff --git a/src/diffeqwrappers.jl b/src/diffeqwrappers.jl
@@ -1,62 +1,90 @@
-function finite_difference!(df::AbstractArray{<:Real}, f, x::AbstractArray{<:Real},
-    fdtype::DataType, ::Type{Val{:Real}}, ::Type{Val{:DiffEqDerivativeWrapper}},
-    fx::Union{Void,AbstractArray{<:Real}}=nothing, epsilon::Union{Void,AbstractArray{<:Real}}=nothing, return_type::DataType=eltype(x))
+function finite_difference!(df::AbstractArray{<:Number}, f, x::Union{Number,AbstractArray{<:Number}},
+    fdtype::DataType, funtype::DataType, ::Type{Val{:DiffEqDerivativeWrapper}},
+    fx::Union{Void,AbstractArray{<:Number}}=nothing, epsilon::Union{Void,AbstractArray{<:Real}}=nothing, return_type::DataType=eltype(x))
 
-    # TODO: test this one, and figure out what happens with epsilon
-    fx1 = f.fx1
-    if fdtype == Val{:forward}
-        epsilon = compute_epsilon(Val{:forward}, x)
-        f(fx, x)
-        f(fx1, x+epsilon)
-        @. df = (fx1 - fx) / epsilon
-    elseif fdtype == Val{:central}
-        epsilon = compute_epsilon(Val{:central}, x)
-        f(fx, x-epsilon)
-        f(fx1, x+epsilon)
-        @. df = (fx1 - fx) / (2 * epsilon)
-    elseif fdtype == Val{:complex}
-        epsilon = eps(eltype(x))
-        f(fx, f(x+im*epsilon))
-        @. df = imag(fx) / epsilon
-    end
+    # TODO: optimized implementations for specific wrappers using the added DiffEq caching where appopriate
+
+    finite_difference!(df, f, x, fdtype, funtype, Val{:Default}, fx, epsilon, return_type)
     df
 end
 
-# AbstractArray{T} should be OK if JacobianWrapper is provided
-function finite_difference_jacobian!(J::AbstractArray{T}, f, x::StridedArray{T}, ::Type{Val{:forward}}, fx::StridedArray{T}, ::Type{Val{:JacobianWrapper}}) where T<:Real
+function finite_difference_jacobian!(J::AbstractMatrix{<:Real}, f, x::AbstractArray{<:Real},
+    fdtype::DataType, ::Type{Val{:Real}}, ::Type{Val{:JacobianWrapper}},
+    fx::AbstractArray{<:Real}, epsilon::Union{Void,AbstractArray{<:Real}}=nothing, return_type::DataType=eltype(x))
+
     m, n = size(J)
-    epsilon_factor = compute_epsilon_factor(Val{:forward}, T)
+    epsilon_elemtype = compute_epsilon_elemtype(epsilon, x)
     x1, fx1 = f.x1, f.fx1
     copy!(x1, x)
-    copy!(fx1, fx)
-    @inbounds for i in 1:n
-        epsilon = compute_epsilon(Val{:forward}, x[i], epsilon_factor)
-        epsilon_inv = one(T) / epsilon
-        x1[i] += epsilon
-        f(fx, x)
-        f(fx1, x1)
-        @. J[:,i] = (fx-fx1) * epsilon_inv
-        x1[i] -= epsilon
+    if fdtype == Val{:forward}
+        epsilon_factor = compute_epsilon_factor(Val{:forward}, epsilon_elemtype)
+        @inbounds for i ∈ 1:n
+            epsilon = compute_epsilon(Val{:forward}, x[i], epsilon_factor)
+            x1[i] += epsilon
+            f(fx1, x1)
+            f(fx, x)
+            @. J[:,i] = (fx1 - fx) / epsilon
+            x1[i] -= epsilon
+        end
+    elseif fdtype == Val{:central}
+        epsilon_factor = compute_epsilon_factor(Val{:central}, epsilon_elemtype)
+        @inbounds for i ∈ 1:n
+            epsilon = compute_epsilon(Val{:central}, x[i], epsilon_factor)
+            x1[i] += epsilon
+            x[i] -= epsilon
+            f(fx1, x1)
+            f(fx, x)
+            @. J[:,i] = (fx1 - fx) / (2*epsilon)
+            x1[i] -= epsilon
+            x[i] += epsilon
+        end
+    elseif fdtype == Val{:complex}
+        x0 = Complex{eltype(x)}(x)
+        epsilon = eps(eltype(x))
+        @inbounds for i ∈ 1:n
+            x0[i] += im * epsilon
+            @. J[:,i] = imag(f(x0)) / epsilon
+            x0[i] -= im * epsilon
+        end
+    else
+        fdtype_error(Val{:Real})
     end
     J
 end
 
-function finite_difference_jacobian!(J::AbstractArray{T}, f, x::StridedArray{T}, ::Type{Val{:central}}, fx::StridedArray{T}, ::Type{Val{:JacobianWrapper}}) where T<:Real
+function finite_difference_jacobian!(J::AbstractMatrix{<:Number}, f, x::AbstractArray{<:Number},
+    fdtype::DataType, ::Type{Val{:Complex}}, ::Type{Val{:JacobianWrapper}},
+    fx::AbstractArray{<:Number}, epsilon::Union{Void,AbstractArray{<:Real}}=nothing, return_type::DataType=eltype(x))
+
+    # TODO: test this
     m, n = size(J)
-    epsilon_factor = compute_epsilon_factor(Val{:central}, T)
+    epsilon_elemtype = compute_epsilon_elemtype(epsilon, x)
     x1, fx1 = f.x1, f.fx1
     copy!(x1, x)
-    copy!(fx1, fx)
-    @inbounds for i in 1:n
-        epsilon = compute_epsilon(Val{:central}, x[i], epsilon_factor)
-        epsilon_double_inv = one(T) / (2 * epsilon)
-        x[i] += epsilon
-        x1[i] -= epsilon
-        f(fx, x)
-        f(fx1, x1)
-        @. J[:,i] = (fx-fx1) * epsilon_double_inv
-        x[i] -= epsilon
-        x1[i] += epsilon
+    if fdtype == Val{:forward}
+        epsilon_factor = compute_epsilon_factor(Val{:forward}, epsilon_elemtype)
+        @inbounds for i ∈ 1:n
+            epsilon = compute_epsilon(Val{:forward}, real(x[i]), epsilon_factor)
+            x1[i] += epsilon
+            f(fx1, x1)
+            f(fx, x)
+            @. J[:,i] = ( real( (fx1 - fx) ) + im*imag( (fx1 - fx) ) ) / epsilon
+            x1[i] -= epsilon
+        end
+    elseif fdtype == Val{:central}
+        epsilon_factor = compute_epsilon_factor(Val{:central}, epsilon_elemtype)
+        @inbounds for i ∈ 1:n
+            epsilon = compute_epsilon(Val{:central}, real(x[i]), epsilon_factor)
+            x1[i] += epsilon
+            x[i] -= epsilon
+            f(fx1, x1)
+            f(fx, x)
+            @. J[:,i] = ( real( (fx1 - fx) ) + im*imag( fx1 - fx ) ) / (2*epsilon)
+            x1[i] -= epsilon
+            x[i] += epsilon
+        end
+    else
+        fdtype_error(Val{:Complex})
     end
     J
 end
diff --git a/src/finitediff.jl b/src/finitediff.jl
@@ -22,6 +22,7 @@ end
     else
         error("Unrecognized fdtype: must be Val{:forward} or Val{:central}.")
     end
+    nothing
 end
 
 function compute_epsilon_elemtype(epsilon, x)
@@ -34,6 +35,18 @@ function compute_epsilon_elemtype(epsilon, x)
     else
         error("Could not compute epsilon type.")
     end
+    nothing
+end
+
+function fdtype_error(funtype::DataType=Val{:Real})
+    if funtype == Val{:Real}
+        error("Unrecognized fdtype: valid values are Val{:forward}, Val{:central} and Val{:complex}.")
+    elseif funtype == Val{:Complex}
+        error("Unrecognized fdtype: valid values are Val{:forward} or Val{:central}.")
+    else
+        error("Unrecognized funtype: valid values are Val{:Real} or Val{:Complex}.")
+    end
+    nothing
 end
 
 include("derivatives.jl")
diff --git a/src/jacobians.jl b/src/jacobians.jl
@@ -15,21 +15,15 @@ function finite_difference_jacobian!(J::AbstractMatrix{<:Number}, f, x::Abstract
 end
 
 function finite_difference_jacobian!(J::AbstractMatrix{<:Real}, f, x::AbstractArray{<:Real},
-    fdtype::DataType, ::Type{Val{:Real}}, wrappertype::DataType=Val{:Default},
+    fdtype::DataType, ::Type{Val{:Real}}, ::Type{Val{:Default}},
     fx::Union{Void,AbstractArray{<:Real}}=nothing, epsilon::Union{Void,AbstractArray{<:Real}}=nothing, returntype=eltype(x))
 
     # TODO: test and rework this
     m, n = size(J)
     epsilon_elemtype = compute_epsilon_elemtype(epsilon, x)
     if fdtype == Val{:forward}
         if typeof(fx) == Void
-            if wrappertype==Val{:Default}
-                fx = f.(x)
-            elseif wrappertype==Val{:DiffEqJacobianWrapper}
-                fx = f(x)
-            else
-                error("Unrecognized wrappertype: must be Val{:Default} or Val{:DiffEqJacobianWrapper}.")
-            end
+            fx = f.(x)
         end
         epsilon_factor = compute_epsilon_factor(Val{:forward}, epsilon_elemtype)
         shifted_x = copy(x)
