Skip to content

Commit a9d8d01

Browse files
ChrisRackauckasChrisRackauckas
committed
Merge remote-tracking branch 'origin/master'
2 parents 7b856bc + 3194b1f commit a9d8d01

File tree

2 files changed

+137
-59
lines changed

2 files changed

+137
-59
lines changed

‎src/jacobians.jl

Lines changed: 57 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -28,42 +28,92 @@ function finite_difference_jacobian!(J::AbstractMatrix{<:Real}, f, x::AbstractAr
2828
fx::Union{Void,AbstractArray{<:Real}}=nothing, epsilon::Union{Void,AbstractArray{<:Real}}=nothing, returntype=eltype(x),
2929
df::Union{Void,AbstractArray{<:Real}}=nothing)
3030

31-
# TODO: test and rework this
31+
# TODO: test and rework this to support GPUArrays and non-indexable types, if possible
3232
m, n = size(J)
3333
epsilon_elemtype = compute_epsilon_elemtype(epsilon, x)
3434
if fdtype == Val{:forward}
3535
if typeof(fx) == Void
36-
fx = f.(x)
36+
fx = f(x)
3737
end
3838
epsilon_factor = compute_epsilon_factor(Val{:forward}, epsilon_elemtype)
3939
shifted_x = copy(x)
4040
@inbounds for i in 1:n
41-
epsilon = compute_epsilon(t, x[i], epsilon_factor)
41+
epsilon = compute_epsilon(Val{:forward}, x[i], epsilon_factor)
4242
shifted_x[i] += epsilon
43-
J[:, i] .= (f(shifted_x) - f_x) / epsilon
43+
J[:, i] .= (f(shifted_x) - fx) / epsilon
4444
shifted_x[i] = x[i]
4545
end
4646
elseif fdtype == Val{:central}
4747
epsilon_factor = compute_epsilon_factor(Val{:central}, epsilon_elemtype)
4848
shifted_x_plus = copy(x)
4949
shifted_x_minus = copy(x)
5050
@inbounds for i in 1:n
51-
epsilon = compute_epsilon(fdtype, x[i], epsilon_factor)
51+
epsilon = compute_epsilon(Val{:central}, x[i], epsilon_factor)
5252
shifted_x_plus[i] += epsilon
5353
shifted_x_minus[i] -= epsilon
5454
J[:, i] .= (f(shifted_x_plus) - f(shifted_x_minus)) / (epsilon + epsilon)
5555
shifted_x_plus[i] = x[i]
5656
shifted_x_minus[i] = x[i]
5757
end
58+
elseif fdtype == Val{:complex}
59+
x0 = Vector{Complex{eltype(x)}}(x)
60+
epsilon = eps(eltype(x))
61+
@inbounds for i in 1:n
62+
x0[i] += im * epsilon
63+
J[:,i] .= imag.(f(x0)) / epsilon
64+
x0[i] -= im * epsilon
65+
end
5866
else
59-
error("Unrecognized fdtype: must be Val{:forward} or Val{:central}.")
67+
fdtype_error(Val{:Real})
68+
end
69+
if typeof(df) != Void
70+
df .= diag(J)
71+
end
72+
J
73+
end
74+
75+
function finite_difference_jacobian!(J::AbstractMatrix{<:Number}, f, x::AbstractArray{<:Number},
76+
fdtype::DataType, ::Type{Val{:Complex}}, ::Type{Val{:Default}},
77+
fx::Union{Void,AbstractArray{<:Number}}=nothing, epsilon::Union{Void,AbstractArray{<:Real}}=nothing, returntype=eltype(x),
78+
df::Union{Void,AbstractArray{<:Number}}=nothing)
79+
80+
# TODO: test and rework this to support GPUArrays and non-indexable types, if possible
81+
m, n = size(J)
82+
epsilon_elemtype = compute_epsilon_elemtype(epsilon, x)
83+
if fdtype == Val{:forward}
84+
if typeof(fx) == Void
85+
fx = f(x)
86+
end
87+
epsilon_factor = compute_epsilon_factor(Val{:forward}, epsilon_elemtype)
88+
shifted_x = copy(x)
89+
@inbounds for i in 1:n
90+
epsilon = compute_epsilon(Val{:forward}, real(x[i]), epsilon_factor)
91+
shifted_x[i] += epsilon
92+
J[:, i] .= ( real.( f(shifted_x) - fx ) + im*imag.( f(shifted_x) - fx ) ) / epsilon
93+
shifted_x[i] = x[i]
94+
end
95+
elseif fdtype == Val{:central}
96+
epsilon_factor = compute_epsilon_factor(Val{:central}, epsilon_elemtype)
97+
shifted_x_plus = copy(x)
98+
shifted_x_minus = copy(x)
99+
@inbounds for i in 1:n
100+
epsilon = compute_epsilon(Val{:central}, real(x[i]), epsilon_factor)
101+
shifted_x_plus[i] += epsilon
102+
shifted_x_minus[i] -= epsilon
103+
J[:, i] .= ( real(f(shifted_x_plus) - f(shifted_x_minus)) + im*imag(f(shifted_x_plus) - f(shifted_x_minus)) ) / (2 * epsilon)
104+
shifted_x_plus[i] = x[i]
105+
shifted_x_minus[i] = x[i]
106+
end
107+
else
108+
fdtype_error(Val{:Complex})
60109
end
61110
if typeof(df) != Void
62111
df .= diag(J)
63112
end
64113
J
65114
end
66115

116+
#=
67117
function finite_difference_jacobian!(J::StridedMatrix{<:Real}, f, x::StridedArray{<:Real},
68118
fdtype::DataType, ::Type{Val{:Real}}, ::Type{Val{:Default}},
69119
fx::Union{Void,StridedArray{<:Real}}=nothing, epsilon::Union{Void,StridedArray{<:Real}}=nothing, returntype=eltype(x),
@@ -179,3 +229,4 @@ function finite_difference_jacobian!(J::StridedMatrix{<:Number}, f, x::StridedAr
179229
end
180230
J
181231
end
232+
=#

‎test/finitedifftests.jl

Lines changed: 80 additions & 53 deletions
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@ J = zeros(J_ref)
88

99
err_func(a,b) = maximum(abs.(a-b))
1010

11-
# TODO: add tests for GPUArrays, those should work now
11+
# TODO: add tests for GPUArrays
1212
# TODO: add tests for DEDataArrays
1313

1414
# StridedArray tests start here
@@ -39,64 +39,13 @@ err_func(a,b) = maximum(abs.(a-b))
3939
@test err_func(DiffEqDiffTools.finite_difference!(df, sin, x, Val{:complex}, Val{:Real}, Val{:Default}, y, epsilon), df_ref) < 1e-15
4040
end
4141

42-
# Jacobian tests for real-valued callables
43-
@time @testset "Jacobian StridedArray real-valued tests" begin
44-
@test err_func(DiffEqDiffTools.finite_difference_jacobian(sin, x, Val{:forward}), J_ref) < 1e-4
45-
@test err_func(DiffEqDiffTools.finite_difference_jacobian(sin, x, Val{:central}), J_ref) < 1e-8
46-
@test err_func(DiffEqDiffTools.finite_difference_jacobian(sin, x, Val{:complex}), J_ref) < 1e-15
47-
48-
@test err_func(DiffEqDiffTools.finite_difference_jacobian(sin, x, Val{:forward}, Val{:Real}, Val{:Default}, y), J_ref) < 1e-4
49-
@test err_func(DiffEqDiffTools.finite_difference_jacobian(sin, x, Val{:central}, Val{:Real}, Val{:Default}, y), J_ref) < 1e-8
50-
@test err_func(DiffEqDiffTools.finite_difference_jacobian(sin, x, Val{:complex}, Val{:Real}, Val{:Default}, y), J_ref) < 1e-15
51-
52-
@test err_func(DiffEqDiffTools.finite_difference_jacobian(sin, x, Val{:forward}, Val{:Real}, Val{:Default}, y, epsilon), J_ref) < 1e-4
53-
@test err_func(DiffEqDiffTools.finite_difference_jacobian(sin, x, Val{:central}, Val{:Real}, Val{:Default}, y, epsilon), J_ref) < 1e-8
54-
@test err_func(DiffEqDiffTools.finite_difference_jacobian(sin, x, Val{:complex}, Val{:Real}, Val{:Default}, y, epsilon), J_ref) < 1e-15
55-
56-
@test err_func(DiffEqDiffTools.finite_difference_jacobian!(J, sin, x, Val{:forward}), J_ref) < 1e-4
57-
@test err_func(DiffEqDiffTools.finite_difference_jacobian!(J, sin, x, Val{:central}), J_ref) < 1e-8
58-
@test err_func(DiffEqDiffTools.finite_difference_jacobian!(J, sin, x, Val{:complex}), J_ref) < 1e-15
59-
60-
@test err_func(DiffEqDiffTools.finite_difference_jacobian!(J, sin, x, Val{:forward}, Val{:Real}, Val{:Default}, y), J_ref) < 1e-4
61-
@test err_func(DiffEqDiffTools.finite_difference_jacobian!(J, sin, x, Val{:central}, Val{:Real}, Val{:Default}, y), J_ref) < 1e-8
62-
@test err_func(DiffEqDiffTools.finite_difference_jacobian!(J, sin, x, Val{:complex}, Val{:Real}, Val{:Default}, y), J_ref) < 1e-15
63-
64-
@test err_func(DiffEqDiffTools.finite_difference_jacobian!(J, sin, x, Val{:forward}, Val{:Real}, Val{:Default}, y, epsilon), J_ref) < 1e-4
65-
@test err_func(DiffEqDiffTools.finite_difference_jacobian!(J, sin, x, Val{:central}, Val{:Real}, Val{:Default}, y, epsilon), J_ref) < 1e-8
66-
@test err_func(DiffEqDiffTools.finite_difference_jacobian!(J, sin, x, Val{:complex}, Val{:Real}, Val{:Default}, y, epsilon), J_ref) < 1e-15
67-
end
68-
69-
# Jacobian tests w/ derivatives (real-valued callables)
70-
@time @testset "Jacobian StridedArray real-valued derivative tests" begin
71-
DiffEqDiffTools.finite_difference_jacobian!(J, df, sin, x, Val{:forward})
72-
@test err_func(df, df_ref) < 1e-4
73-
DiffEqDiffTools.finite_difference_jacobian!(J, df, sin, x, Val{:central})
74-
@test err_func(df, df_ref) < 1e-8
75-
DiffEqDiffTools.finite_difference_jacobian!(J, df, sin, x, Val{:complex})
76-
@test err_func(df, df_ref) < 1e-15
77-
DiffEqDiffTools.finite_difference_jacobian!(J, df, sin, x, Val{:forward}, Val{:Real}, Val{:Default}, y)
78-
@test err_func(df, df_ref) < 1e-4
79-
DiffEqDiffTools.finite_difference_jacobian!(J, df, sin, x, Val{:central}, Val{:Real}, Val{:Default}, y)
80-
@test err_func(df, df_ref) < 1e-8
81-
DiffEqDiffTools.finite_difference_jacobian!(J, df, sin, x, Val{:complex}, Val{:Real}, Val{:Default}, y)
82-
@test err_func(df, df_ref) < 1e-15
83-
DiffEqDiffTools.finite_difference_jacobian!(J, df, sin, x, Val{:forward}, Val{:Real}, Val{:Default}, y, epsilon)
84-
@test err_func(df, df_ref) < 1e-4
85-
DiffEqDiffTools.finite_difference_jacobian!(J, df, sin, x, Val{:central}, Val{:Real}, Val{:Default}, y, epsilon)
86-
@test err_func(df, df_ref) < 1e-8
87-
DiffEqDiffTools.finite_difference_jacobian!(J, df, sin, x, Val{:complex}, Val{:Real}, Val{:Default}, y, epsilon)
88-
@test err_func(df, df_ref) < 1e-15
89-
end
90-
9142
# derivative tests for complex-valued callables
9243
x = x + im*x
9344
f(x) = cos(real(x)) + im*sin(imag(x))
9445
y = f.(x)
9546
df = zeros(x)
9647
epsilon = zeros(length(x))
9748
df_ref = -sin.(real(x)) + im*cos.(imag(x))
98-
J_ref = diagm(df_ref)
99-
J = zeros(J_ref)
10049

10150
@time @testset "Derivative StridedArray complex-valued tests" begin
10251
@test err_func(DiffEqDiffTools.finite_difference(f, x, Val{:forward}, Val{:Complex}), df_ref) < 1e-4
@@ -118,6 +67,84 @@ J = zeros(J_ref)
11867
@test err_func(DiffEqDiffTools.finite_difference!(df, f, x, Val{:central}, Val{:Complex}, Val{:Default}, y, epsilon), df_ref) < 1e-8
11968
end
12069

70+
function f(x)
71+
fvec = zeros(x)
72+
fvec[1] = (x[1]+3)*(x[2]^3-7)+18
73+
fvec[2] = sin(x[2]*exp(x[1])-1)
74+
fvec
75+
end
76+
x = rand(2)
77+
y = f(x)
78+
J_ref = [[-7+x[2]^3 3*(3+x[1])*x[2]^2]; [exp(x[1])*x[2]*cos(1-exp(x[1])*x[2]) exp(x[1])*cos(1-exp(x[1])*x[2])]]
79+
J = zeros(J_ref)
80+
df = zeros(x)
81+
df_ref = diag(J_ref)
82+
epsilon = zeros(x)
83+
84+
# Jacobian tests for real-valued callables
85+
@time @testset "Jacobian StridedArray real-valued tests" begin
86+
@test err_func(DiffEqDiffTools.finite_difference_jacobian(f, x, Val{:forward}), J_ref) < 1e-4
87+
@test err_func(DiffEqDiffTools.finite_difference_jacobian(f, x, Val{:central}), J_ref) < 1e-8
88+
@test err_func(DiffEqDiffTools.finite_difference_jacobian(f, x, Val{:complex}), J_ref) < 1e-14
89+
90+
@test err_func(DiffEqDiffTools.finite_difference_jacobian(f, x, Val{:forward}, Val{:Real}, Val{:Default}, y), J_ref) < 1e-4
91+
@test err_func(DiffEqDiffTools.finite_difference_jacobian(f, x, Val{:central}, Val{:Real}, Val{:Default}, y), J_ref) < 1e-8
92+
@test err_func(DiffEqDiffTools.finite_difference_jacobian(f, x, Val{:complex}, Val{:Real}, Val{:Default}, y), J_ref) < 1e-14
93+
94+
@test err_func(DiffEqDiffTools.finite_difference_jacobian(f, x, Val{:forward}, Val{:Real}, Val{:Default}, y, epsilon), J_ref) < 1e-4
95+
@test err_func(DiffEqDiffTools.finite_difference_jacobian(f, x, Val{:central}, Val{:Real}, Val{:Default}, y, epsilon), J_ref) < 1e-8
96+
@test err_func(DiffEqDiffTools.finite_difference_jacobian(f, x, Val{:complex}, Val{:Real}, Val{:Default}, y, epsilon), J_ref) < 1e-14
97+
98+
@test err_func(DiffEqDiffTools.finite_difference_jacobian!(J, f, x, Val{:forward}), J_ref) < 1e-4
99+
@test err_func(DiffEqDiffTools.finite_difference_jacobian!(J, f, x, Val{:central}), J_ref) < 1e-8
100+
@test err_func(DiffEqDiffTools.finite_difference_jacobian!(J, f, x, Val{:complex}), J_ref) < 1e-14
101+
102+
@test err_func(DiffEqDiffTools.finite_difference_jacobian!(J, f, x, Val{:forward}, Val{:Real}, Val{:Default}, y), J_ref) < 1e-4
103+
@test err_func(DiffEqDiffTools.finite_difference_jacobian!(J, f, x, Val{:central}, Val{:Real}, Val{:Default}, y), J_ref) < 1e-8
104+
@test err_func(DiffEqDiffTools.finite_difference_jacobian!(J, f, x, Val{:complex}, Val{:Real}, Val{:Default}, y), J_ref) < 1e-14
105+
106+
@test err_func(DiffEqDiffTools.finite_difference_jacobian!(J, f, x, Val{:forward}, Val{:Real}, Val{:Default}, y, epsilon), J_ref) < 1e-4
107+
@test err_func(DiffEqDiffTools.finite_difference_jacobian!(J, f, x, Val{:central}, Val{:Real}, Val{:Default}, y, epsilon), J_ref) < 1e-8
108+
@test err_func(DiffEqDiffTools.finite_difference_jacobian!(J, f, x, Val{:complex}, Val{:Real}, Val{:Default}, y, epsilon), J_ref) < 1e-14
109+
end
110+
111+
# Jacobian tests w/ derivatives (real-valued callables)
112+
@time @testset "Jacobian StridedArray real-valued derivative tests" begin
113+
DiffEqDiffTools.finite_difference_jacobian!(J, df, f, x, Val{:forward})
114+
@show df, df_ref
115+
@test err_func(df, df_ref) < 1e-4
116+
DiffEqDiffTools.finite_difference_jacobian!(J, df, f, x, Val{:central})
117+
@test err_func(df, df_ref) < 1e-8
118+
DiffEqDiffTools.finite_difference_jacobian!(J, df, f, x, Val{:complex})
119+
@test err_func(df, df_ref) < 1e-15
120+
DiffEqDiffTools.finite_difference_jacobian!(J, df, f, x, Val{:forward}, Val{:Real}, Val{:Default}, y)
121+
@test err_func(df, df_ref) < 1e-4
122+
DiffEqDiffTools.finite_difference_jacobian!(J, df, f, x, Val{:central}, Val{:Real}, Val{:Default}, y)
123+
@test err_func(df, df_ref) < 1e-8
124+
DiffEqDiffTools.finite_difference_jacobian!(J, df, f, x, Val{:complex}, Val{:Real}, Val{:Default}, y)
125+
@test err_func(df, df_ref) < 1e-15
126+
DiffEqDiffTools.finite_difference_jacobian!(J, df, f, x, Val{:forward}, Val{:Real}, Val{:Default}, y, epsilon)
127+
@test err_func(df, df_ref) < 1e-4
128+
DiffEqDiffTools.finite_difference_jacobian!(J, df, f, x, Val{:central}, Val{:Real}, Val{:Default}, y, epsilon)
129+
@test err_func(df, df_ref) < 1e-8
130+
DiffEqDiffTools.finite_difference_jacobian!(J, df, f, x, Val{:complex}, Val{:Real}, Val{:Default}, y, epsilon)
131+
@test err_func(df, df_ref) < 1e-15
132+
end
133+
134+
function f(x)
135+
fvec = zeros(x)
136+
fvec[1] = (im*x[1]+3)*(x[2]^3-7)+18
137+
fvec[2] = sin(x[2]*exp(x[1])-1)
138+
fvec
139+
end
140+
x = rand(2) + im*rand(2)
141+
y = f(x)
142+
J_ref = [[im*(-7+x[2]^3) 3*(3+im*x[1])*x[2]^2]; [exp(x[1])*x[2]*cos(1-exp(x[1])*x[2]) exp(x[1])*cos(1-exp(x[1])*x[2])]]
143+
J = zeros(J_ref)
144+
df = zeros(x)
145+
df_ref = diag(J_ref)
146+
epsilon = zeros(real.(x))
147+
121148
# Jacobian tests for complex-valued callables
122149
@time @testset "Jacobian StridedArray complex-valued tests" begin
123150
@test err_func(DiffEqDiffTools.finite_difference_jacobian(f, x, Val{:forward}, Val{:Complex}, Val{:Default}), J_ref) < 1e-4
@@ -139,7 +166,7 @@ end
139166
@test err_func(DiffEqDiffTools.finite_difference_jacobian!(J, f, x, Val{:central}, Val{:Complex}, Val{:Default}, y, epsilon), J_ref) < 1e-8
140167
end
141168

142-
# Jacobian tests w/ derivatives (real-valued callables)
169+
# Jacobian tests w/ derivatives (complex-valued callables)
143170
@time @testset "Jacobian StridedArray complex-valued derivative tests" begin
144171
DiffEqDiffTools.finite_difference_jacobian!(J, df, f, x, Val{:forward}, Val{:Complex}, Val{:Default})
145172
@test err_func(df, df_ref) < 1e-4

0 commit comments

Comments
 (0)