Fix state before refactoring everything.

Author: dextorious

src/derivatives.jl

@@ -1,29 +1,26 @@
 #=
-Derivative of f : R -> R or f : C -> C at a single point x.
+Sinple-point derivatives of scalar->scalar maps.
 =#
-function finite_difference_derivative(f, x::T, fdtype::DataType, funtype::DataType=Val{:Real},
-    f_x::Union{Void,T}=nothing) where T<:Number
+function finite_difference_derivative(f, x::T, fdtype::DataType=Val{:central},
+    returntype::DataType=eltype(x), f_x::Union{Void,T}=nothing) where T<:Number
 
-    if funtype == Val{:Real}
-        if fdtype == Val{:complex}
-            epsilon = eps(T)
-        else
-            epsilon = compute_epsilon(fdtype, x)
-        end
-    elseif funtype == Val{:Complex}
-        epsilon = compute_epsilon(fdtype, real(x))
-    else
-        fdtype_error(funtype)
+    epsilon = compute_epsilon(fdtype, real(x))
+    if fdtype == Val{:forward}
+        return (f(x+epsilon) - f(x)) / epsilon
+    elseif fdtype == Val{:central}
+        return (f(x+epsilon) - f(x-epsilon)) / (2*epsilon)
+    elseif fdtype == Val{:complex} && returntype == Val{:Real}
+        return imag(f(x+im*epsilon)) / epsilon
     end
-
-    _finite_difference_kernel(f, x, fdtype, funtype, epsilon, f_x)
+    fdtype_error(returntype)
 end
 
 #=
-Finite difference kernels for single point derivatives of f : R -> R.
-These are currently underused because of inlining / broadcast issues.
+Finite difference kernels for single point derivatives.
+These are currently unused because of inlining / broadcast issues.
 Revisit this in Julia v0.7 / 1.0.
 =#
+#=
 @inline function _finite_difference_kernel(f, x::T, ::Type{Val{:forward}}, ::Type{Val{:Real}},
     epsilon::T, fx::Union{Void,T}=nothing) where T<:Real
 
@@ -61,6 +58,7 @@ end
 
     real(f(x+epsilon) - f(x-epsilon)) / (2 * epsilon) + im*imag(f(x+im*epsilon) - f(x-im*epsilon)) / (2 * epsilon)
 end
+=#
 # Single point derivative implementations end here.
 
 
@@ -179,9 +177,9 @@ function _finite_difference_derivative!(df::AbstractArray{<:Number}, f, x::Abstr
         if typeof(fx) == Void
             fx = f.(x)
         end
-        @. df = real((f(x+epsilon) - fx)) / epsilon + im*imag((f(x+im*epsilon) - fx)) / epsilon
+        @. df = real((f(x+epsilon) - fx)) / epsilon + im*imag((f(x+epsilon) - fx)) / epsilon
     elseif fdtype == Val{:central}
-        @. df = real(f(x+epsilon) - f(x-epsilon)) / (2 * epsilon) + im*imag(f(x+im*epsilon) - f(x-im*epsilon)) / (2 * epsilon)
+        @. df = real(f(x+epsilon) - f(x-epsilon)) / (2 * epsilon) + im*imag(f(x+epsilon) - f(x-epsilon)) / (2 * epsilon)
     else
         fdtype_error(Val{:Complex})
     end
@@ -242,13 +240,13 @@ function _finite_difference_derivative!(df::StridedArray{<:Number}, f, x::Stride
             else
                 fxi = fx[i]
             end
-            df[i] = ( real( f(x[i]+epsilon) - fxi ) + im*imag( f(x[i]+im*epsilon) - fxi ) ) / epsilon
+            df[i] = ( real( f(x[i]+epsilon) - fxi ) + im*imag( f(x[i]+epsilon) - fxi ) ) / epsilon
         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}, 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)
+            df[i] = ( real( f(x[i]+epsilon) - f(x[i]-epsilon) ) + im*imag( f(x[i]+epsilon) - f(x[i]-epsilon) ) ) / (2 * epsilon)
         end
     else
         fdtype_error(Val{:Complex})

src/finitediff.jl

@@ -14,13 +14,17 @@ end
     eps_cbrt * max(one(T), abs(x))
 end
 
+@inline function compute_epsilon(::Type{Val{:complex}}, x::T) where T<:Real
+    eps(T)
+end
+
 @inline function compute_epsilon_factor(fdtype::DataType, ::Type{T}) where T<:Number
     if fdtype==Val{:forward}
         return sqrt(eps(T))
     elseif fdtype==Val{:central}
         return cbrt(eps(T))
     else
-        error("Unrecognized fdtype $fdtype: must be Val{:forward} or Val{:central}.")
+        return one(T)
     end
 end
 

src/gradients.jl

@@ -195,43 +195,25 @@ function finite_difference_gradient!(df::AbstractArray{<:Number}, f, x::Abstract
             c2[i] += epsilon
             if typeof(fx) == Void
                 df[i]  = real(f(c2) - f(x)) / epsilon
-            else
-                df[i]  = real(f(c2) - fx) / epsilon
-            end
-            c2[i] -= epsilon
-            c2[i] += im*epsilon
-            if typeof(fx) == Void
                 df[i] += im*imag(f(c2) - f(x)) / epsilon
             else
+                df[i]  = real(f(c2) - fx) / epsilon
                 df[i] += im*imag(f(c2) - fx) / epsilon
             end
-            c2[i] -= im*epsilon
+            c2[i] -= epsilon
         end
     elseif fdtype == Val{:central}
         @inbounds for i ∈ eachindex(x)
             epsilon = c1[i]
             c2[i]  += epsilon
             x[i]   -= epsilon
             df[i]   = real(f(c2) - f(x)) / (2*epsilon)
-            c2[i]  -= c1[i]
-            x[i]   += c1[i]
-            c2[i]  += im*epsilon
-            x[i]   -= im*epsilon
             df[i]  += im*imag(f(c2) - f(x)) / (2*epsilon)
-            c2[i]  -= im*epsilon
-            x[i]   += im*epsilon
-        end
-    elseif fdtype == Val{:complex}
-        epsilon_elemtype = compute_epsilon_elemtype(nothing, x)
-        epsilon_complex = eps(epsilon_elemtype)
-        # we use c1 here to avoid typing issues with x
-        @inbounds for i ∈ eachindex(x)
-            c1[i] += im*epsilon_complex
-            df[i]  = imag(f(c1)) / epsilon_complex
-            c1[i] -= im*epsilon_complex
+            c2[i]  -= c1[i]
+            x[i]   += c1[i]            
         end
     else
-        fdtype_error(Val{:Complex})
+        fdtype_error(Val{:complex})
     end
     df
 end
@@ -241,6 +223,34 @@ end
 function finite_difference_gradient!(df::AbstractArray{<:Number}, f, x::Number,
     cache::GradientCache{T1,T2,T3,fdtype,Val{:Complex}}) where {T1,T2,T3,fdtype}
 
-    error("Not implemented yet.")
+    # NOTE: in this case epsilon is a scalar, we need two arrays for fx1 and fx2
+    # c1 denotes fx1, c2 is fx2, sizes guaranteed by the cache constructor
+    fx, c1, c2 = cache.fx, cache.c1, cache.c2
+
+    epsilon_elemtype = compute_epsilon_elemtype(nothing, x)
+    if fdtype == Val{:forward}
+        epsilon_factor = compute_epsilon_factor(fdtype, real(eltype(x)))
+        epsilon = compute_epsilon(Val{:forward}, real(x), epsilon_factor)
+        c1 .= f(x+epsilon)
+        if typeof(fx) != Void
+            @. df  = real(c1 - fx) / epsilon
+            @. df += im*imag(c1 - fx) / epsilon
+        else
+            c2 .= f(x)
+            @. df  = real(c1 - c2) / epsilon
+            @. df += im*imag(c1 - c2) / epsilon
+        end
+    elseif fdtype == Val{:central}
+        epsilon_factor = compute_epsilon_factor(fdtype, real(eltype(x)))
+        epsilon = compute_epsilon(Val{:central}, real(x), epsilon_factor)
+        c1 .= f(x+epsilon)
+        c2 .= f(x-epsilon)
+        @. df  = real(c1 - c2) / (2*epsilon)
+        #c1 .= im*imag(f(x+im*epsilon))
+        #c2 .= im*imag(f(x-im*epsilon))
+        @. df += im*imag(c1 - c2) / (2*epsilon)
+    elseif fdtype == Val{:complex}
+        fdtype_error(Val{:complex})
+    end
     df
 end

test/finitedifftests.jl

@@ -52,17 +52,19 @@ end
 end
 
 x = x + im*x
-f(x) = cos(real(x)) + im*sin(imag(x))
+#f(x) = cos(real(x)) + im*sin(imag(x))
+f(x) = sin(x) + cos(x)
 y = f.(x)
 df = zeros(x)
 epsilon = zeros(length(x))
-df_ref = -sin.(real(x)) + im*cos.(imag(x))
+#df_ref = -sin.(real(x)) + im*cos.(imag(x))
+df_ref = cos.(x) - sin.(x)
 forward_cache = DiffEqDiffTools.DerivativeCache(x, y, epsilon, Val{:forward})
 central_cache = DiffEqDiffTools.DerivativeCache(x, nothing, epsilon, Val{:central})
 
 @time @testset "Derivative single point complex-valued tests" begin
-    @test err_func(DiffEqDiffTools.finite_difference_derivative(f, π/4+im*π/4, Val{:forward}, Val{:Complex}), -√2/2 + im*√2/2) < 1e-4
-    @test err_func(DiffEqDiffTools.finite_difference_derivative(f, π/4+im*π/4, Val{:central}, Val{:Complex}), -√2/2 + im*√2/2) < 1e-8
+    @test err_func(DiffEqDiffTools.finite_difference_derivative(f, π/4+im*π/4, Val{:forward}, Val{:Complex}), cos(π/4+im*π/4)-sin(π/4+im*π/4)) < 1e-4
+    @test err_func(DiffEqDiffTools.finite_difference_derivative(f, π/4+im*π/4, Val{:central}, Val{:Complex}), cos(π/4+im*π/4)-sin(π/4+im*π/4)) < 1e-8
 end
 
 @time @testset "Derivative StridedArray complex-valued tests" begin
@@ -157,8 +159,26 @@ complex_cache = DiffEqDiffTools.GradientCache(df,x,fx,nothing,nothing,Val{:compl
     @test err_func(DiffEqDiffTools.finite_difference_gradient!(df, f, x, complex_cache), df_ref) < 1e-15
 end
 
+f(x) = [sin(x), cos(x)]
+@show x = (2π * rand()) * (1 + im)
+@show fx = f(x)
+df = zeros(fx)
+@show df_ref = [cos(x), -sin(x)]
+forward_cache = DiffEqDiffTools.GradientCache(df,x,fx,nothing,nothing,Val{:forward},Val{:Real})
+central_cache = DiffEqDiffTools.GradientCache(df,x,fx,nothing,nothing,Val{:central},Val{:Real})
+complex_cache = DiffEqDiffTools.GradientCache(df,x,fx,nothing,nothing,Val{:complex},Val{:Real})
+@show DiffEqDiffTools.finite_difference_gradient(f, x, Val{:forward})
+@show DiffEqDiffTools.finite_difference_gradient(f, x, Val{:central})
+
 @time @testset "Gradient of f:R->R^n complex-valued tests" begin
-    # TODO
+    @test err_func(DiffEqDiffTools.finite_difference_gradient(f, x, Val{:forward}), df_ref) < 1e-4
+    @test err_func(DiffEqDiffTools.finite_difference_gradient(f, x, Val{:central}), df_ref) < 1e-7
+
+    @test err_func(DiffEqDiffTools.finite_difference_gradient!(df, f, x, Val{:forward}), df_ref) < 1e-4
+    @test err_func(DiffEqDiffTools.finite_difference_gradient!(df, f, x, Val{:central}), df_ref) < 1e-7
+
+    @test err_func(DiffEqDiffTools.finite_difference_gradient!(df, f, x, forward_cache), df_ref) < 1e-4
+    @test err_func(DiffEqDiffTools.finite_difference_gradient!(df, f, x, central_cache), df_ref) < 1e-7
 end
 
 # Jacobian tests