fix jacobian in case where target functions's returned value is alias…

…ed with input (JuliaMath#106)
johanbluecreek · Mar 16, 2017 · f7b5556 · f7b5556
1 parent 2ec7168
commit f7b5556
Show file tree

Hide file tree

Showing 2 changed files with 18 additions and 12 deletions.
diff --git 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
@@ -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()
 #