Fix thread safety issue in finite_difference_jacobian! when using

central differencing.

Fix initialization of HessianCache

Test thread safety for gradients, Jacobians, and Hessians

Author: andreasnoack

src/hessians.jl

@@ -11,9 +11,10 @@ function HessianCache(xpp,xpm,xmp,xmm,
     HessianCache{typeof(xpp),fdtype,inplace}(xpp,xpm,xmp,xmm)
 end
 
-function HessianCache(x,fdtype=Val{:hcentral},
-                        inplace = x isa StaticArray ? Val{false} : Val{true})
-    HessianCache{typeof(x),fdtype,inplace}(copy(x),copy(x),copy(x),copy(x))
+function HessianCache(x, fdtype=Val{:hcentral},
+                      inplace = x isa StaticArray ? Val{false} : Val{true})
+    cx = copy(x)
+    HessianCache{typeof(cx),fdtype,inplace}(cx, copy(x), copy(x), copy(x))
 end
 
 function finite_difference_hessian(f, x,

src/jacobians.jl

@@ -356,13 +356,13 @@ function finite_difference_jacobian!(
             if sparsity isa Nothing
                 x1_save = ArrayInterface.allowed_getindex(x1,color_i)
                 epsilon = compute_epsilon(Val{:forward}, x1_save, relstep, absstep, dir)
-                ArrayInterface.allowed_setindex!(x1,x1_save + epsilon,color_i)
+                ArrayInterface.allowed_setindex!(x1, x1_save + epsilon, color_i)
                 f(fx1, x1)
                 # J is dense, so either it is truly dense or this is the
                 # compressed form of the coloring, so write into it.
                 @. J[:,color_i] = (vfx1 - vfx) / epsilon
                 # Now return x1 back to its original value
-                ArrayInterface.allowed_setindex!(x1,x1_save,color_i)
+                ArrayInterface.allowed_setindex!(x1, x1_save, color_i)
             else # Perturb along the colorvec vector
                 @. fx1 = x1 * (_color == color_i)
                 tmp = norm(fx1)
@@ -381,9 +381,9 @@ function finite_difference_jacobian!(
                     += means requires a zero'd out start
                     =#
                     if J isa SparseMatrixCSC
-                        @. void_setindex!((J.nzval,),getindex((J.nzval,),rows_index) + (getindex((_color,),cols_index) == color_i) * getindex((vfx1,),rows_index),rows_index)
+                        @. void_setindex!((J.nzval,), getindex((J.nzval,), rows_index) + (getindex((_color,), cols_index) == color_i) * getindex((vfx1,), rows_index), rows_index)
                     else
-                        @. void_setindex!((J,),getindex((J,),rows_index, cols_index) + (getindex((_color,),cols_index) == color_i) * getindex((vfx1,),rows_index),rows_index, cols_index)
+                        @. void_setindex!((J,), getindex((J,), rows_index, cols_index) + (getindex((_color,), cols_index) == color_i) * getindex((vfx1,), rows_index), rows_index, cols_index)
                     end
                 end
                 # Now return x1 back to its original value
@@ -394,16 +394,14 @@ function finite_difference_jacobian!(
         vfx1 = _vec(fx1)
         @inbounds for color_i ∈ 1:maximum(colorvec)
             if sparsity isa Nothing
-                x_save = ArrayInterface.allowed_getindex(x,color_i)
-                x1_save = ArrayInterface.allowed_getindex(x1,color_i)
+                x_save = ArrayInterface.allowed_getindex(x, color_i)
                 epsilon = compute_epsilon(Val{:central}, x_save, relstep, absstep, dir)
-                ArrayInterface.allowed_setindex!(x1,x1_save+epsilon,color_i)
-                ArrayInterface.allowed_setindex!(x,x_save-epsilon,color_i)
+                ArrayInterface.allowed_setindex!(x1, x_save + epsilon, color_i)
                 f(fx1, x1)
-                f(fx, x)
+                ArrayInterface.allowed_setindex!(x1, x_save - epsilon, color_i)
+                f(fx, x1)
                 @. J[:,color_i] = (vfx1 - vfx) / 2epsilon
-                ArrayInterface.allowed_setindex!(x1,x1_save,color_i)
-                ArrayInterface.allowed_setindex!(x,x_save,color_i)
+                ArrayInterface.allowed_setindex!(x1, x_save, color_i)
             else # Perturb along the colorvec vector
                 @. fx1 = x1 * (_color == color_i)
                 tmp = norm(fx1)
@@ -417,9 +415,9 @@ function finite_difference_jacobian!(
                     _colorediteration!(J,sparsity,rows_index,cols_index,vfx1,colorvec,color_i,n)
                 else
                     if J isa SparseMatrixCSC
-                        @. void_setindex!((J.nzval,),getindex((J.nzval,),rows_index) + (getindex((_color,),cols_index) == color_i) * getindex((vfx1,),rows_index),rows_index)
+                        @. void_setindex!((J.nzval,), getindex((J.nzval,), rows_index) + (getindex((_color,), cols_index) == color_i) * getindex((vfx1,), rows_index), rows_index)
                     else
-                        @. void_setindex!((J,),getindex((J,),rows_index, cols_index) + (getindex((_color,),cols_index) == color_i) * getindex((vfx1,),rows_index),rows_index, cols_index)
+                        @. void_setindex!((J,), getindex((J,), rows_index, cols_index) + (getindex((_color,), cols_index) == color_i) * getindex((vfx1,), rows_index), rows_index, cols_index)
                     end
                 end
                 @. x1 = x1 - epsilon * (_color == color_i)
@@ -430,11 +428,11 @@ function finite_difference_jacobian!(
         epsilon = eps(eltype(x))
         @inbounds for color_i ∈ 1:maximum(colorvec)
             if sparsity isa Nothing
-                x1_save = ArrayInterface.allowed_getindex(x1,color_i)
-                ArrayInterface.allowed_setindex!(x1,x1_save + im*epsilon, color_i)
+                x1_save = ArrayInterface.allowed_getindex(x1, color_i)
+                ArrayInterface.allowed_setindex!(x1, x1_save + im*epsilon, color_i)
                 f(fx,x1)
                 @. J[:,color_i] = imag(vfx) / epsilon
-                ArrayInterface.allowed_setindex!(x1,x1_save,color_i)
+                ArrayInterface.allowed_setindex!(x1, x1_save,color_i)
             else # Perturb along the colorvec vector
                 @. x1 = x1 + im * epsilon * (_color == color_i)
                 f(fx,x1)
@@ -443,9 +441,9 @@ function finite_difference_jacobian!(
                     _colorediteration!(J,sparsity,rows_index,cols_index,vfx,colorvec,color_i,n)
                 else
                    if J isa SparseMatrixCSC
-                        @. void_setindex!((J.nzval,),getindex((J.nzval,),rows_index) + (getindex((_color,),cols_index) == color_i) * getindex((vfx,),rows_index),rows_index)
+                        @. void_setindex!((J.nzval,), getindex((J.nzval,), rows_index) + (getindex((_color,), cols_index) == color_i) * getindex((vfx,),rows_index), rows_index)
                     else
-                        @. void_setindex!((J,),getindex((J,),rows_index, cols_index) + (getindex((_color,),cols_index) == color_i) * getindex((vfx,),rows_index),rows_index, cols_index)
+                        @. void_setindex!((J,), getindex((J,), rows_index, cols_index) + (getindex((_color,), cols_index) == color_i) * getindex((vfx,), rows_index), rows_index, cols_index)
                     end
                 end
                 @. x1 = x1 - im * epsilon * (_color == color_i)

test/finitedifftests.jl

@@ -425,3 +425,33 @@ H = similar(H_ref)
     @test err_func(FiniteDiff.finite_difference_hessian!(H, f, x), H_ref) < 1e-4
     @test err_func(FiniteDiff.finite_difference_hessian!(H, f, x, hcache), H_ref) < 1e-4
 end
+
+# Thread safety
+# create an an abstract array type that doesn't allow setindex
+struct ImmutableVector <: DenseVector{Float64}
+    x::Vector{Float64}
+end
+Base.size(x::ImmutableVector) = size(x.x)
+Base.getindex(x::ImmutableVector, i::Integer) = x.x[i]
+@testset "thread safety" begin
+    @testset "Gradients with diff type $difftype"  for difftype in (Val{:forward}, Val{:central}, Val{:complex})
+        g = FiniteDiff.finite_difference_gradient(sum, ImmutableVector(ones(2)), difftype)
+        @test g ≈ ones(2)
+        FiniteDiff.finite_difference_gradient!(g, sum, ImmutableVector(ones(2)), difftype)
+        @test g ≈ ones(2)
+    end
+
+    @testset "Hessians (only supported diff type is :hcentral)" begin
+        H = FiniteDiff.finite_difference_hessian(t -> sum(abs2, t)/2, ImmutableVector(ones(2)))
+        @test H ≈ Matrix(I, 2, 2)
+        FiniteDiff.finite_difference_hessian!(parent(H), t -> sum(abs2, t)/2, ImmutableVector(ones(2)))
+        @test H ≈ Matrix(I, 2, 2)
+    end
+
+    @testset "Jacobians with diff type $difftype"  for difftype in (Val{:forward}, Val{:central}, Val{:complex})
+        J = FiniteDiff.finite_difference_jacobian(identity, ImmutableVector(ones(2)), difftype)
+        @test J ≈ Matrix(I, 2, 2)
+        FiniteDiff.finite_difference_jacobian!(J, (out, in) -> out .= in, ImmutableVector(ones(2)), difftype)
+        @test J ≈ Matrix(I, 2, 2)
+    end
+end