Complex-mode f:R^n->R gradient fixes and tests.

Author: dextorious

src/gradients.jl

@@ -25,20 +25,26 @@ function GradientCache(
 
     if typeof(x) <: AbstractArray # the f:R^n->R case
         # need cache arrays for epsilon (c1) and x1 (c2)
-        epsilon_elemtype = compute_epsilon_elemtype(nothing, x)
-        if typeof(c1) == Void || eltype(c1) != epsilon_elemtype
-            _c1 = zeros(epsilon_elemtype, size(x))
+        if fdtype != Val{:complex} # complex-mode FD only needs one cache, for x+eps*im
+            epsilon_elemtype = compute_epsilon_elemtype(nothing, x)
+            if typeof(c1) == Void || eltype(c1) != epsilon_elemtype
+                _c1 = zeros(epsilon_elemtype, size(x))
+            else
+                _c1 = c1
+            end
+            epsilon_factor = compute_epsilon_factor(fdtype, real(eltype(x)))
+            @. _c1 = compute_epsilon(fdtype, real(x), epsilon_factor)
+
+            if typeof(c2) != typeof(x) || size(c2) != size(x)
+                _c2 = copy(x)
+            else
+                copy!(_c2, x)
+            end
         else
-            _c1 = c1
+            _c1 = x + 0*im
+            _c2 = nothing
         end
-        epsilon_factor = compute_epsilon_factor(fdtype, real(eltype(x)))
-        @. _c1 = compute_epsilon(fdtype, real(x), epsilon_factor)
 
-        if typeof(c2) != typeof(x) || size(c2) != size(x)
-            _c2 = copy(x)
-        else
-            copy!(_c2, x)
-        end
     else # the f:R->R^n case
         # need cache arrays for fx1 and fx2
         if typeof(c1) != typeof(df) || size(c1) != size(df)
@@ -117,7 +123,14 @@ function finite_difference_gradient!(df::AbstractArray{<:Number}, f, x::Abstract
             x[i]  += c1[i]
         end
     elseif fdtype == Val{:complex}
-        # TODO
+        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
+        end
     end
     df
 end

test/finitedifftests.jl

@@ -1,3 +1,9 @@
+
+# TODO: add tests for GPUArrays
+# TODO: add tests for DEDataArrays
+
+
+# Derivative tests
 x = collect(linspace(-2π, 2π, 100))
 y = sin.(x)
 df = zeros(100)
@@ -9,9 +15,6 @@ complex_cache = DiffEqDiffTools.DerivativeCache(x, nothing, nothing, Val{:comple
 
 err_func(a,b) = maximum(abs.(a-b))
 
-# TODO: add tests for GPUArrays
-# TODO: add tests for DEDataArrays
-
 @time @testset "Derivative single point real-valued tests" begin
     @test err_func(DiffEqDiffTools.finite_difference_derivative(sin, π/4, Val{:forward}), √2/2) < 1e-4
     @test err_func(DiffEqDiffTools.finite_difference_derivative(sin, π/4, Val{:central}), √2/2) < 1e-8
@@ -85,30 +88,43 @@ end
     @test err_func(DiffEqDiffTools.finite_difference_derivative!(df, f, x, central_cache), df_ref) < 1e-8
 end
 
+# Gradient tests
 f(x) = 2x[1] + x[2]^2
 x = rand(2)
 fx = f(x)
 df = zeros(2)
 df_ref = [2., 2*x[2]]
 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})
+complex_cache = DiffEqDiffTools.GradientCache(df,x,fx,nothing,nothing,Val{:complex},Val{:Real})
 
 @time @testset "Gradient of f:R^n->R real-valued tests" begin
     @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-8
-    #@test err_func(DiffEqDiffTools.finite_difference_gradient(f, x, Val{:complex}), df_ref) < 1e-15
+    @test err_func(DiffEqDiffTools.finite_difference_gradient(f, x, Val{:complex}), df_ref) < 1e-15
 
     @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-8
-    #@test err_func(DiffEqDiffTools.finite_difference_gradient!(df, f, x, Val{:complex}), df_ref) < 1e-15
+    @test err_func(DiffEqDiffTools.finite_difference_gradient!(df, f, x, Val{:complex}), df_ref) < 1e-15
 
     @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-8
-    #@test err_func(DiffEqDiffTools.finite_difference_gradient!(df, f, x, complex_cache), df_ref) < 1e-15
+    @test err_func(DiffEqDiffTools.finite_difference_gradient!(df, f, x, complex_cache), df_ref) < 1e-15
+end
+
+@time @testset "Gradient of f:R^n->R complex-valued tests" begin
+    # TODO
 end
 
+@time @testset "Gradient of f:R->R^n real-valued tests" begin
+    # TODO
+end
+
+@time @testset "Gradient of f:R->R^n complex-valued tests" begin
+    # TODO
+end
 
+# Jacobian tests
 function f(fvec,x)
     fvec[1] = (x[1]+3)*(x[2]^3-7)+18
     fvec[2] = sin(x[2]*exp(x[1])-1)