Gradient bugfixes.

Author: dextorious

src/DiffEqDiffTools.jl

@@ -5,6 +5,7 @@ module DiffEqDiffTools
 include("function_wrappers.jl")
 include("finitediff.jl")
 include("derivatives.jl")
+include("gradients.jl")
 include("jacobians.jl")
 
 end # module

src/finitediff.jl

@@ -20,7 +20,7 @@ end
     elseif fdtype==Val{:central}
         return cbrt(eps(T))
     else
-        error("Unrecognized fdtype: must be Val{:forward} or Val{:central}.")
+        error("Unrecognized fdtype $fdtype: must be Val{:forward} or Val{:central}.")
     end
 end
 

src/gradients.jl

@@ -1,4 +1,4 @@
-struct GradientCache{CacheType, CacheType2, CacheType3, fdtype, RealOrComplex}
+struct GradientCache{CacheType1, CacheType2, CacheType3, fdtype, RealOrComplex}
     fx :: CacheType1
     c1 :: CacheType2
     c2 :: CacheType3
@@ -7,19 +7,15 @@ end
 function GradientCache(
     df     :: AbstractArray{<:Number},
     x      :: Union{<:Number, AbstractArray{<:Number}},
-    fx     :: Union{Void,AbstractArray{<:Number}} = nothing,
-    c1     :: AbstractArray{<:Number} = nothing,
-    c2     :: AbstractArray{<:Number} = nothing,
+    fx     :: Union{Void,<:Number,AbstractArray{<:Number}} = nothing,
+    c1     :: Union{Void,AbstractArray{<:Number}} = nothing,
+    c2     :: Union{Void,AbstractArray{<:Number}} = nothing,
     fdtype :: DataType = Val{:central},
     RealOrComplex :: DataType =
-        fdtype==Val{:complex} ? Val{:Real} : eltype(x1) <: Complex ? Val{:Complex} : Val{:Real}
+        fdtype==Val{:complex} ? Val{:Real} : eltype(x) <: Complex ? Val{:Complex} : Val{:Real}
     )
 
-    if fdtype == Val{:complex} && RealOrComplex == Val{:Complex}
-        fdtype_error(Val{:Complex})
-    end
-
-    if fdtype != Val{:forward}
+    if fdtype != Val{:forward} && typeof(fx) != Void
         warn("Pre-computed function values are only useful for fdtype == Val{:forward}.")
         _fx = nothing
     else
@@ -28,8 +24,8 @@ function GradientCache(
     end
 
     if typeof(x) <: AbstractArray # the f:R^n->R case
-        # need cache arrays for epsilon and x1
-        epsilon_elemtype = compute_epsilon_elemtype(epsilon, x)
+        # 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))
         else
@@ -56,7 +52,7 @@ function GradientCache(
             _c2 = c2
         end
     end
-    GradientCache{typeof(fx),typeof(c1),typeof(c2),fdtype,RealOrComplex}(fx,c1,c2)
+    GradientCache{typeof(_fx),typeof(_c1),typeof(_c2),fdtype,RealOrComplex}(_fx,_c1,_c2)
 end
 
 function finite_difference_gradient(f, x, fdtype::DataType=Val{:central},
@@ -107,20 +103,18 @@ function finite_difference_gradient!(df::AbstractArray{<:Number}, f, x::Abstract
     # c2 is x1, pre-set to the values of x by the cache constructor
     fx, c1, c2 = cache.fx, cache.c1, cache.c2
     if fdtype == Val{:forward}
-        # TODO: do we even need c2 for the forward case?
         @inbounds for i ∈ eachindex(x)
             c2[i] += c1[i]
             df[i]  = (f(c2) - f(x)) / c1[i]
             c2[i] -= c1[i]
         end
     elseif fdtype == Val{:central}
         @inbounds for i ∈ eachindex(x)
-            #copy!(x1,x)
             c2[i] += c1[i]
             x[i]  -= c1[i]
             df[i]  = (f(c2) - f(x)) / (2*c1[i])
             c2[i] -= c1[i]
-            x[i]  += c1[i] # revert any changes to x
+            x[i]  += c1[i]
         end
     elseif fdtype == Val{:complex}
         # TODO

test/finitedifftests.jl

@@ -85,6 +85,30 @@ end
     @test err_func(DiffEqDiffTools.finite_difference_derivative!(df, f, x, central_cache), df_ref) < 1e-8
 end
 
+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})
+
+@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!(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, 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
+end
+
+
 function f(fvec,x)
     fvec[1] = (x[1]+3)*(x[2]^3-7)+18
     fvec[2] = sin(x[2]*exp(x[1])-1)