From 2ae3604a055923e727188bf31a9720e3f97771f9 Mon Sep 17 00:00:00 2001 From: Jarrett Revels Date: Wed, 15 Mar 2017 18:39:05 -0400 Subject: [PATCH] fix jacobian in case where target functions's returned value is aliased with input --- src/finite_difference.jl | 23 +++++++++++------------ test/derivative.jl | 7 +++++++ 2 files changed, 18 insertions(+), 12 deletions(-) diff --git a/src/finite_difference.jl b/src/finite_difference.jl index db42fec..e7a9081 100644 --- a/src/finite_difference.jl +++ b/src/finite_difference.jl @@ -166,24 +166,23 @@ function finite_difference_jacobian!{R <: Number, # Iterate over each dimension of the gradient separately. if dtype == :forward + shifted_x = copy(x) for i = 1:n @forwardrule x[i] epsilon - oldx = x[i] - x[i] = oldx + epsilon - f_xplusdx = f(x) - x[i] = oldx - J[:, i] = (f_xplusdx - f_x) / epsilon + shifted_x[i] += epsilon + J[:, i] = (f(shifted_x) - f_x) / epsilon + shifted_x[i] = x[i] end elseif dtype == :central + shifted_x_plus = copy(x) + shifted_x_minus = copy(x) for i = 1:n @centralrule x[i] epsilon - oldx = x[i] - x[i] = oldx + epsilon - f_xplusdx = f(x) - x[i] = oldx - epsilon - f_xminusdx = f(x) - x[i] = oldx - J[:, i] = (f_xplusdx - f_xminusdx) / (epsilon + epsilon) + shifted_x_plus[i] += epsilon + shifted_x_minus[i] -= epsilon + J[:, i] = (f(shifted_x_plus) - f(shifted_x_minus)) / (epsilon + epsilon) + shifted_x_plus[i] = x[i] + shifted_x_minus[i] = x[i] end else error("dtype must :forward or :central") diff --git a/test/derivative.jl b/test/derivative.jl index c482cc1..02117d0 100644 --- a/test/derivative.jl +++ b/test/derivative.jl @@ -31,6 +31,13 @@ f4(x::Vector) = (100.0 - x[1])^2 + (50.0 - x[2])^2 @test norm(Calculus.gradient(f4, :central)([100.0, 50.0]) - [0.0, 0.0]) < 10e-4 @test norm(Calculus.gradient(f4)([100.0, 50.0]) - [0.0, 0.0]) < 10e-4 +# +# jacobian() +# + +@test norm(Calculus.jacobian(identity, rand(3), :forward) - eye(3)) < 10e-4 +@test norm(Calculus.jacobian(identity, rand(3), :central) - eye(3)) < 10e-4 + # # second_derivative() #