Add the option to get the derivatives from the Jacobians without extr…

…a computation.
JuliaDiff · Nov 8, 2017 · f788176 · f788176
1 parent 4803b91
commit f788176
Show file tree

Hide file tree

Showing 2 changed files with 75 additions and 7 deletions.
diff --git a/src/jacobians.jl b/src/jacobians.jl
@@ -1,22 +1,32 @@
 # Compute the Jacobian matrix of a real-valued callable f.
 function finite_difference_jacobian(f, x::AbstractArray{<:Number},
     fdtype::DataType=Val{:central}, funtype::DataType=Val{:Real}, wrappertype::DataType=Val{:Default},
-    fx::Union{Void,AbstractArray{<:Number}}=nothing, epsilon::Union{Void,AbstractArray{<:Real}}=nothing, returntype=eltype(x))
+    fx::Union{Void,AbstractArray{<:Number}}=nothing, epsilon::Union{Void,AbstractArray{<:Real}}=nothing, returntype=eltype(x),
+    df::Union{Void,AbstractArray{<:Number}}=nothing)
 
     J = zeros(returntype, length(x), length(x))
-    finite_difference_jacobian!(J, f, x, fdtype, funtype, wrappertype, fx, epsilon, returntype)
+    finite_difference_jacobian!(J, f, x, fdtype, funtype, wrappertype, fx, epsilon, returntype, df)
 end
 
-function finite_difference_jacobian!(J::AbstractMatrix{<:Number}, f, x::AbstractArray{<:Number},
+function finite_difference_jacobian!(J::AbstractMatrix{<:Number}, df::AbstractVector, f, x::AbstractArray{<:Number},
     fdtype::DataType=Val{:central}, funtype::DataType=Val{:Real}, wrappertype::DataType=Val{:Default},
     fx::Union{Void,AbstractArray{<:Number}}=nothing, epsilon::Union{Void,AbstractArray{<:Number}}=nothing, returntype=eltype(x))
 
-    finite_difference_jacobian!(J, f, x, fdtype, funtype, wrappertype, fx, epsilon, returntype)
+    finite_difference_jacobian!(J, f, x, fdtype, funtype, wrappertype, fx, epsilon, returntype, df)
+end
+
+function finite_difference_jacobian!(J::AbstractMatrix{<:Number}, f, x::AbstractArray{<:Number},
+    fdtype::DataType=Val{:central}, funtype::DataType=Val{:Real}, wrappertype::DataType=Val{:Default},
+    fx::Union{Void,AbstractArray{<:Number}}=nothing, epsilon::Union{Void,AbstractArray{<:Number}}=nothing, returntype=eltype(x),
+    df::Union{Void,AbstractArray{<:Number}}=nothing)
+
+    finite_difference_jacobian!(J, f, x, fdtype, funtype, wrappertype, fx, epsilon, returntype, df)
 end
 
 function finite_difference_jacobian!(J::AbstractMatrix{<:Real}, f, x::AbstractArray{<:Real},
     fdtype::DataType, ::Type{Val{:Real}}, ::Type{Val{:Default}},
-    fx::Union{Void,AbstractArray{<:Real}}=nothing, epsilon::Union{Void,AbstractArray{<:Real}}=nothing, returntype=eltype(x))
+    fx::Union{Void,AbstractArray{<:Real}}=nothing, epsilon::Union{Void,AbstractArray{<:Real}}=nothing, returntype=eltype(x),
+    df::Union{Void,AbstractArray{<:Real}}=nothing)
 
     # TODO: test and rework this
     m, n = size(J)
@@ -48,12 +58,16 @@ function finite_difference_jacobian!(J::AbstractMatrix{<:Real}, f, x::AbstractAr
     else
         error("Unrecognized fdtype: must be Val{:forward} or Val{:central}.")
     end
+    if typeof(df) != Void
+        df .= diag(J)
+    end
     J
 end
 
 function finite_difference_jacobian!(J::StridedMatrix{<:Real}, f, x::StridedArray{<:Real},
     fdtype::DataType, ::Type{Val{:Real}}, ::Type{Val{:Default}},
-    fx::Union{Void,StridedArray{<:Real}}=nothing, epsilon::Union{Void,StridedArray{<:Real}}=nothing, returntype=eltype(x))
+    fx::Union{Void,StridedArray{<:Real}}=nothing, epsilon::Union{Void,StridedArray{<:Real}}=nothing, returntype=eltype(x),
+    df::Union{Void,StridedArray{<:Real}}=nothing)
 
     m, n = size(J)
     epsilon_elemtype = compute_epsilon_elemtype(epsilon, x)
@@ -73,6 +87,9 @@ function finite_difference_jacobian!(J::StridedMatrix{<:Real}, f, x::StridedArra
                             J[j,i] = (f(x[j]+epsilon) - fx[j]) * epsilon_inv
                         end
                     end
+                    if typeof(df) != Void
+                        df[i] = J[j,i]
+                    end
                 else
                     J[j,i] = zero(returntype)
                 end
@@ -86,6 +103,9 @@ function finite_difference_jacobian!(J::StridedMatrix{<:Real}, f, x::StridedArra
             for j in 1:m
                 if i==j
                     J[j,i] = (f(x[j]+epsilon) - f(x[j]-epsilon)) * epsilon_double_inv
+                    if typeof(df) != Void
+                        df[i] = J[j,i]
+                    end
                 else
                     J[j,i] = zero(returntype)
                 end
@@ -98,6 +118,9 @@ function finite_difference_jacobian!(J::StridedMatrix{<:Real}, f, x::StridedArra
             for j in 1:m
                 if i==j
                     J[j,i] = imag(f(x[j]+im*epsilon)) * epsilon_inv
+                    if typeof(df) != Void
+                        df[i] = J[j,i]
+                    end
                 else
                     J[j,i] = zero(returntype)
                 end
@@ -109,7 +132,8 @@ end
 
 function finite_difference_jacobian!(J::StridedMatrix{<:Number}, f, x::StridedArray{<:Number},
     fdtype::DataType, ::Type{Val{:Complex}}, ::Type{Val{:Default}},
-    fx::Union{Void,StridedArray{<:Number}}=nothing, epsilon::Union{Void,StridedArray{<:Real}}=nothing, returntype=eltype(x))
+    fx::Union{Void,StridedArray{<:Number}}=nothing, epsilon::Union{Void,StridedArray{<:Real}}=nothing, returntype=eltype(x),
+    df::Union{Void,StridedArray{<:Number}}=nothing)
 
     # TODO: finish this
     m, n = size(J)
@@ -126,6 +150,9 @@ function finite_difference_jacobian!(J::StridedMatrix{<:Number}, f, x::StridedAr
                     else
                         J[j,i] = ( real( f(x[j]+epsilon) - fx[j] ) + im*imag( f(x[j]+im*epsilon) - fx[j] ) ) * epsilon_inv
                     end
+                    if typeof(df) != Void
+                        df[i] = J[j,i]
+                    end
                 else
                     J[j,i] = zero(returntype)
                 end
@@ -139,6 +166,9 @@ function finite_difference_jacobian!(J::StridedMatrix{<:Number}, f, x::StridedAr
             for j in 1:m
                 if i==j
                     J[j,i] = ( real( f(x[j]+epsilon)-f(x[j]-epsilon) ) + im*imag( f(x[j]+im*epsilon) - f(x[j]-im*epsilon) ) ) * epsilon_double_inv
+                    if typeof(df) != Void
+                        df[i] = J[j,i]
+                    end
                 else
                     J[j,i] = zero(returntype)
                 end

diff --git a/test/finitedifftests.jl b/test/finitedifftests.jl
@@ -66,6 +66,28 @@ end
     @test err_func(DiffEqDiffTools.finite_difference_jacobian!(J, sin, x, Val{:complex}, Val{:Real}, Val{:Default}, y, epsilon), J_ref) < 1e-15
 end
 
+# Jacobian tests w/ derivatives (real-valued callables)
+@time @testset "Jacobian StridedArray real-valued derivative tests" begin
+    DiffEqDiffTools.finite_difference_jacobian!(J, df, sin, x, Val{:forward})
+    @test err_func(df, df_ref) < 1e-4
+    DiffEqDiffTools.finite_difference_jacobian!(J, df, sin, x, Val{:central})
+    @test err_func(df, df_ref) < 1e-8
+    DiffEqDiffTools.finite_difference_jacobian!(J, df, sin, x, Val{:complex})
+    @test err_func(df, df_ref) < 1e-15
+    DiffEqDiffTools.finite_difference_jacobian!(J, df, sin, x, Val{:forward}, Val{:Real}, Val{:Default}, y)
+    @test err_func(df, df_ref) < 1e-4
+    DiffEqDiffTools.finite_difference_jacobian!(J, df, sin, x, Val{:central}, Val{:Real}, Val{:Default}, y)
+    @test err_func(df, df_ref) < 1e-8
+    DiffEqDiffTools.finite_difference_jacobian!(J, df, sin, x, Val{:complex}, Val{:Real}, Val{:Default}, y)
+    @test err_func(df, df_ref) < 1e-15
+    DiffEqDiffTools.finite_difference_jacobian!(J, df, sin, x, Val{:forward}, Val{:Real}, Val{:Default}, y, epsilon)
+    @test err_func(df, df_ref) < 1e-4
+    DiffEqDiffTools.finite_difference_jacobian!(J, df, sin, x, Val{:central}, Val{:Real}, Val{:Default}, y, epsilon)
+    @test err_func(df, df_ref) < 1e-8
+    DiffEqDiffTools.finite_difference_jacobian!(J, df, sin, x, Val{:complex}, Val{:Real}, Val{:Default}, y, epsilon)
+    @test err_func(df, df_ref) < 1e-15
+end
+
 # derivative tests for complex-valued callables
 x = x + im*x
 f(x) = cos(real(x)) + im*sin(imag(x))
@@ -116,4 +138,20 @@ end
     @test err_func(DiffEqDiffTools.finite_difference_jacobian!(J, f, x, Val{:forward}, Val{:Complex}, Val{:Default}, y, epsilon), J_ref) < 1e-4
     @test err_func(DiffEqDiffTools.finite_difference_jacobian!(J, f, x, Val{:central}, Val{:Complex}, Val{:Default}, y, epsilon), J_ref) < 1e-8
 end
+
+# Jacobian tests w/ derivatives (real-valued callables)
+@time @testset "Jacobian StridedArray complex-valued derivative tests" begin
+    DiffEqDiffTools.finite_difference_jacobian!(J, df, f, x, Val{:forward}, Val{:Complex}, Val{:Default})
+    @test err_func(df, df_ref) < 1e-4
+    DiffEqDiffTools.finite_difference_jacobian!(J, df, f, x, Val{:central}, Val{:Complex}, Val{:Default})
+    @test err_func(df, df_ref) < 1e-8
+    DiffEqDiffTools.finite_difference_jacobian!(J, df, f, x, Val{:forward}, Val{:Complex}, Val{:Default}, y)
+    @test err_func(df, df_ref) < 1e-4
+    DiffEqDiffTools.finite_difference_jacobian!(J, df, f, x, Val{:central}, Val{:Complex}, Val{:Default}, y)
+    @test err_func(df, df_ref) < 1e-8
+    DiffEqDiffTools.finite_difference_jacobian!(J, df, f, x, Val{:forward}, Val{:Complex}, Val{:Default}, y, epsilon)
+    @test err_func(df, df_ref) < 1e-4
+    DiffEqDiffTools.finite_difference_jacobian!(J, df, f, x, Val{:central}, Val{:Complex}, Val{:Default}, y, epsilon)
+    @test err_func(df, df_ref) < 1e-8
+end
 # StridedArray tests end here