@@ -6,15 +6,15 @@ Very heavily inspired by Calculus.jl, but with an emphasis on performance and Di
6
6
Compute the finite difference interval epsilon.
7
7
Reference: Numerical Recipes, chapter 5.7.
8
8
=#
9
- @inline function compute_epsilon {T<:Real} (:: Type{Val{:forward}} , x:: T , eps_sqrt:: T = sqrt (eps (T)))
9
+ @inline function compute_epsilon (:: Type{Val{:forward}} , x:: T , eps_sqrt:: T = sqrt (eps (T))) where T <: Real
10
10
eps_sqrt * max (one (T), abs (x))
11
11
end
12
12
13
- @inline function compute_epsilon {T<:Real} (:: Type{Val{:central}} , x:: T , eps_cbrt:: T = cbrt (eps (T)))
13
+ @inline function compute_epsilon (:: Type{Val{:central}} , x:: T , eps_cbrt:: T = cbrt (eps (T))) where T <: Real
14
14
eps_cbrt * max (one (T), abs (x))
15
15
end
16
16
17
- @inline function compute_epsilon_factor {T<:Real} (fdtype:: DataType , :: Type{T} )
17
+ @inline function compute_epsilon_factor (fdtype:: DataType , :: Type{T} ) where T <: Real
18
18
if fdtype== Val{:forward }
19
19
return sqrt (eps (T))
20
20
elseif fdtype== Val{:central }
@@ -30,12 +30,12 @@ Compute the derivative df of a real-valued callable f on a collection of points
30
30
Generic fallbacks for AbstractArrays that are not StridedArrays.
31
31
# TODO : test the fallbacks
32
32
=#
33
- function finite_difference {T<:Real} (f, x:: AbstractArray{T} , fdtype:: DataType , fx:: Union{Void,AbstractArray{T}} = nothing , funtype:: DataType = Val{:Default })
33
+ function finite_difference (f, x:: AbstractArray{T} , fdtype:: DataType , fx:: Union{Void,AbstractArray{T}} = nothing , funtype:: DataType = Val{:Default }) where T <: Real
34
34
df = zeros (T, size (x))
35
35
finite_difference! (df, f, x, fdtype, fx, funtype)
36
36
end
37
37
38
- function finite_difference! {T<:Real} (df:: AbstractArray{T} , f, x:: AbstractArray{T} , fdtype:: DataType , fx:: Union{Void,AbstractArray{T}} , :: Type{Val{:Default}} )
38
+ function finite_difference! (df:: AbstractArray{T} , f, x:: AbstractArray{T} , fdtype:: DataType , fx:: Union{Void,AbstractArray{T}} , :: Type{Val{:Default}} ) where T <: Real
39
39
if fdtype == Val{:forward }
40
40
epsilon_factor = compute_epsilon_factor (fdtype, T)
41
41
@. epsilon = compute_epsilon (fdtype, x, epsilon_factor)
@@ -55,7 +55,7 @@ function finite_difference!{T<:Real}(df::AbstractArray{T}, f, x::AbstractArray{T
55
55
df
56
56
end
57
57
58
- function finite_difference! {T<:Real} (df:: AbstractArray{T} , f, x:: T , fdtype:: DataType , fx:: AbstractArray{T} , :: Type{Val{:DiffEqDerivativeWrapper}} )
58
+ function finite_difference! (df:: AbstractArray{T} , f, x:: T , fdtype:: DataType , fx:: AbstractArray{T} , :: Type{Val{:DiffEqDerivativeWrapper}} ) where T <: Real
59
59
fx1 = f. fx1
60
60
if fdtype == Val{:forward }
61
61
epsilon = compute_epsilon (fdtype, x)
79
79
Compute the derivative df of a real-valued callable f on a collection of points x.
80
80
Optimized implementations for StridedArrays.
81
81
=#
82
- function finite_difference! {T<:Real} (df:: StridedArray{T} , f, x:: StridedArray{T} , :: Type{Val{:central}} , :: Union{Void,StridedArray{T}} , :: Type{Val{:Default}} )
82
+ function finite_difference! (df:: StridedArray{T} , f, x:: StridedArray{T} , :: Type{Val{:central}} , :: Union{Void,StridedArray{T}} , :: Type{Val{:Default}} ) where T <: Real
83
83
epsilon_factor = compute_epsilon_factor (Val{:central }, T)
84
84
@inbounds for i in 1 : length (x)
85
85
epsilon = compute_epsilon (Val{:central }, x[i], epsilon_factor)
@@ -90,7 +90,7 @@ function finite_difference!{T<:Real}(df::StridedArray{T}, f, x::StridedArray{T},
90
90
df
91
91
end
92
92
93
- function finite_difference! {T<:Real} (df:: StridedArray{T} , f, x:: StridedArray{T} , :: Type{Val{:forward}} , :: Void , :: Type{Val{:Default}} )
93
+ function finite_difference! (df:: StridedArray{T} , f, x:: StridedArray{T} , :: Type{Val{:forward}} , :: Void , :: Type{Val{:Default}} ) where T <: Real
94
94
epsilon_factor = compute_epsilon_factor (Val{:forward }, T)
95
95
@inbounds for i in 1 : length (x)
96
96
epsilon = compute_epsilon (Val{:forward }, x[i], epsilon_factor)
@@ -101,7 +101,7 @@ function finite_difference!{T<:Real}(df::StridedArray{T}, f, x::StridedArray{T},
101
101
df
102
102
end
103
103
104
- function finite_difference! {T<:Real} (df:: StridedArray{T} , f, x:: StridedArray{T} , :: Type{Val{:forward}} , fx:: StridedArray{T} , :: Type{Val{:Default}} )
104
+ function finite_difference! (df:: StridedArray{T} , f, x:: StridedArray{T} , :: Type{Val{:forward}} , fx:: StridedArray{T} , :: Type{Val{:Default}} ) where T <: Real
105
105
epsilon_factor = compute_epsilon_factor (Val{:forward }, T)
106
106
@inbounds for i in 1 : length (x)
107
107
epsilon = compute_epsilon (Val{:forward }, x[i], epsilon_factor)
116
116
Compute the derivative df of a real-valued callable f on a collection of points x.
117
117
Single point implementations.
118
118
=#
119
- function finite_difference {T<:Real} (f, x:: T , fdtype:: DataType , f_x:: Union{Void,T} = nothing )
119
+ function finite_difference (f, x:: T , fdtype:: DataType , f_x:: Union{Void,T} = nothing ) where T <: Real
120
120
if fdtype == Val{:complex }
121
121
epsilon = eps (T)
122
122
return imag (f (x+ im* epsilon)) / epsilon
@@ -126,15 +126,15 @@ function finite_difference{T<:Real}(f, x::T, fdtype::DataType, f_x::Union{Void,T
126
126
end
127
127
end
128
128
129
- @inline function finite_difference_kernel {T<:Real} (f, x:: T , :: Type{Val{:forward}} , epsilon:: T , fx:: Union{Void,T} )
129
+ @inline function finite_difference_kernel (f, x:: T , :: Type{Val{:forward}} , epsilon:: T , fx:: Union{Void,T} ) where T <: Real
130
130
if typeof (fx) == Void
131
131
return (f (x+ epsilon) - f (x)) / epsilon
132
132
else
133
133
return (f (x+ epsilon) - fx) / epsilon
134
134
end
135
135
end
136
136
137
- @inline function finite_difference_kernel {T<:Real} (f, x:: T , :: Type{Val{:central}} , epsilon:: T , :: Union{Void,T} = nothing )
137
+ @inline function finite_difference_kernel (f, x:: T , :: Type{Val{:central}} , epsilon:: T , :: Union{Void,T} = nothing ) where T <: Real
138
138
(f (x+ epsilon) - f (x- epsilon)) / (2 * epsilon)
139
139
end
140
140
144
144
#=
145
145
Compute the Jacobian matrix of a real-valued callable f: R^n -> R^m.
146
146
=#
147
- function finite_difference_jacobian {T<:Real} (f, x:: AbstractArray{T} , fdtype:: DataType = Val{:central }, funtype:: DataType = Val{:Default })
147
+ function finite_difference_jacobian (f, x:: AbstractArray{T} , fdtype:: DataType = Val{:central }, funtype:: DataType = Val{:Default }) where T <: Real
148
148
if funtype== Val{:Default }
149
149
fx = f .(x)
150
150
elseif funtype== Val{:DiffEqJacobianWrapper }
@@ -156,7 +156,7 @@ function finite_difference_jacobian{T<:Real}(f, x::AbstractArray{T}, fdtype::Dat
156
156
finite_difference_jacobian! (J, f, x, fdtype, fx, funtype)
157
157
end
158
158
159
- function finite_difference_jacobian! {T<:Real} (J:: AbstractArray{T} , f, x:: AbstractArray{T} , fdtype:: DataType , fx:: AbstractArray{T} , :: DataType )
159
+ function finite_difference_jacobian! (J:: AbstractArray{T} , f, x:: AbstractArray{T} , fdtype:: DataType , fx:: AbstractArray{T} , :: DataType ) where T <: Real
160
160
# This is an inefficient fallback that only makes sense if setindex/getindex are unavailable, e.g. GPUArrays etc.
161
161
m, n = size (J)
162
162
epsilon_factor = compute_epsilon_factor (fdtype, T)
@@ -185,7 +185,7 @@ function finite_difference_jacobian!{T<:Real}(J::AbstractArray{T}, f, x::Abstrac
185
185
J
186
186
end
187
187
188
- function finite_difference_jacobian! {T<:Real} (J:: StridedArray{T} , f, x:: StridedArray{T} , :: Type{Val{:central}} , fx:: StridedArray{T} , :: Type{Val{:Default}} )
188
+ function finite_difference_jacobian! (J:: StridedArray{T} , f, x:: StridedArray{T} , :: Type{Val{:central}} , fx:: StridedArray{T} , :: Type{Val{:Default}} ) where T <: Real
189
189
m, n = size (J)
190
190
epsilon_factor = compute_epsilon_factor (Val{:central }, T)
191
191
@inbounds for i in 1 : n
@@ -202,7 +202,7 @@ function finite_difference_jacobian!{T<:Real}(J::StridedArray{T}, f, x::StridedA
202
202
J
203
203
end
204
204
205
- function finite_difference_jacobian! {T<:Real} (J:: StridedArray{T} , f, x:: StridedArray{T} , :: Type{Val{:forward}} , fx:: StridedArray{T} , :: Type{Val{:Default}} )
205
+ function finite_difference_jacobian! (J:: StridedArray{T} , f, x:: StridedArray{T} , :: Type{Val{:forward}} , fx:: StridedArray{T} , :: Type{Val{:Default}} ) where T <: Real
206
206
m, n = size (J)
207
207
epsilon_factor = compute_epsilon_factor (Val{:forward }, T)
208
208
@inbounds for i in 1 : n
@@ -219,7 +219,7 @@ function finite_difference_jacobian!{T<:Real}(J::StridedArray{T}, f, x::StridedA
219
219
J
220
220
end
221
221
222
- function finite_difference_jacobian! {T<:Real} (J:: StridedArray{T} , f, x:: StridedArray{T} , :: Type{Val{:complex}} , fx:: StridedArray{T} , :: Type{Val{:Default}} )
222
+ function finite_difference_jacobian! (J:: StridedArray{T} , f, x:: StridedArray{T} , :: Type{Val{:complex}} , fx:: StridedArray{T} , :: Type{Val{:Default}} ) where T <: Real
223
223
m, n = size (J)
224
224
epsilon = eps (T)
225
225
epsilon_inv = one (T) / epsilon
@@ -236,7 +236,7 @@ function finite_difference_jacobian!{T<:Real}(J::StridedArray{T}, f, x::StridedA
236
236
end
237
237
238
238
# efficient implementations for OrdinaryDiffEq Jacobian wrappers, assuming the system function supplies StridedArrays
239
- function finite_difference_jacobian! {T<:Real} (J:: StridedArray{T} , f, x:: StridedArray{T} , :: Type{Val{:forward}} , fx:: StridedArray{T} , :: Type{Val{:JacobianWrapper}} )
239
+ function finite_difference_jacobian! (J:: StridedArray{T} , f, x:: StridedArray{T} , :: Type{Val{:forward}} , fx:: StridedArray{T} , :: Type{Val{:JacobianWrapper}} ) where T <: Real
240
240
m, n = size (J)
241
241
epsilon_factor = compute_epsilon_factor (Val{:forward }, T)
242
242
x1, fx1 = f. x1, f. fx1
@@ -254,7 +254,7 @@ function finite_difference_jacobian!{T<:Real}(J::StridedArray{T}, f, x::StridedA
254
254
J
255
255
end
256
256
257
- function finite_difference_jacobian! {T<:Real} (J:: StridedArray{T} , f, x:: StridedArray{T} , :: Type{Val{:central}} , fx:: StridedArray{T} , :: Type{Val{:JacobianWrapper}} )
257
+ function finite_difference_jacobian! (J:: StridedArray{T} , f, x:: StridedArray{T} , :: Type{Val{:central}} , fx:: StridedArray{T} , :: Type{Val{:JacobianWrapper}} ) where T <: Real
258
258
m, n = size (J)
259
259
epsilon_factor = compute_epsilon_factor (Val{:central }, T)
260
260
x1, fx1 = f. x1, f. fx1
0 commit comments