1- struct GradientCache{CacheType1,CacheType2,CacheType3,CacheType4,fdtype,returntype,inplace}
1+ struct GradientCache{
2+ CacheType1, CacheType2, CacheType3, CacheType4, fdtype, returntype, inplace}
23 fx:: CacheType1
34 c1:: CacheType2
45 c2:: CacheType3
1617Allocating Cache Constructor
1718"""
1819function GradientCache (
19- df,
20- x,
21- fdtype= Val (:central ),
22- returntype= eltype (df),
23- inplace= Val (true ))
24-
20+ df,
21+ x,
22+ fdtype = Val (:central ),
23+ returntype = eltype (df),
24+ inplace = Val (true ))
2525 fdtype isa Type && (fdtype = fdtype ())
2626 inplace isa Type && (inplace = inplace ())
2727 if typeof (x) <: AbstractArray # the vector->scalar case
@@ -59,9 +59,8 @@ function GradientCache(
5959 _c3 = x
6060 end
6161
62- GradientCache{Nothing,typeof (_c1),typeof (_c2),typeof (_c3),fdtype,
63- returntype,inplace}(nothing , _c1, _c2, _c3)
64-
62+ GradientCache{Nothing, typeof (_c1), typeof (_c2), typeof (_c3), fdtype,
63+ returntype, inplace}(nothing , _c1, _c2, _c3)
6564end
6665
6766"""
@@ -103,14 +102,14 @@ GradientCache{Float64, Vector{Float64}, Vector{Float64}, Vector{Float64}, Val{:c
103102```
104103"""
105104function GradientCache (
106- fx:: Fx ,# match order in struct for Setfield
107- c1:: T ,
108- c2:: T ,
109- c3:: T ,
110- fdtype= Val (:central ),
111- returntype= eltype (fx),
112- inplace= Val (true )) where {T,Fx} # Val(false) isn't so important for vector -> scalar, it gets ignored in that case anyway.
113- GradientCache {Fx,T,T,T, fdtype,returntype,inplace} (fx, c1, c2, c3)
105+ fx:: Fx ,# match order in struct for Setfield
106+ c1:: T ,
107+ c2:: T ,
108+ c3:: T ,
109+ fdtype = Val (:central ),
110+ returntype = eltype (fx),
111+ inplace = Val (true )) where {T, Fx} # Val(false) isn't so important for vector -> scalar, it gets ignored in that case anyway.
112+ GradientCache {Fx, T, T, T, fdtype, returntype, inplace} (fx, c1, c2, c3)
114113end
115114
116115"""
@@ -124,23 +123,60 @@ end
124123 absstep=relstep,
125124 dir=true)
126125
127- Gradients are either a vector->scalar map `f(x)`, or a scalar->vector map `f(fx,x)` if `inplace=Val{true}` and `fx=f(x)` if `inplace=Val{false}` .
126+ Compute the gradient of function `f` at point `x` using finite differences .
128127
129- Cache-less.
128+ This is the cache-less version that allocates temporary arrays internally.
129+ Supports both vector→scalar maps `f(x) → scalar` and scalar→vector maps depending
130+ on the `inplace` parameter and function signature.
131+
132+ # Arguments
133+ - `f`: Function to differentiate
134+ - If `typeof(x) <: AbstractArray`: `f(x)` should return a scalar (vector→scalar gradient)
135+ - If `typeof(x) <: Number` and `inplace=Val(true)`: `f(fx, x)` modifies `fx` in-place (scalar→vector gradient)
136+ - If `typeof(x) <: Number` and `inplace=Val(false)`: `f(x)` returns a vector (scalar→vector gradient)
137+ - `x`: Point at which to evaluate the gradient (vector or scalar)
138+ - `fdtype::Type{T1}=Val{:central}`: Finite difference method (`:forward`, `:central`, `:complex`)
139+ - `returntype::Type{T2}=eltype(x)`: Element type of gradient components
140+ - `inplace::Type{Val{T3}}=Val{true}`: Whether to use in-place function evaluation
141+
142+ # Keyword Arguments
143+ - `relstep`: Relative step size (default: method-dependent optimal value)
144+ - `absstep=relstep`: Absolute step size fallback
145+ - `dir=true`: Direction for step size (typically ±1)
146+
147+ # Returns
148+ - Gradient vector `∇f` where `∇f[i] = ∂f/∂x[i]`
149+
150+ # Examples
151+ ```julia
152+ # Vector→scalar gradient
153+ f(x) = x[1]^2 + x[2]^2
154+ x = [1.0, 2.0]
155+ grad = finite_difference_gradient(f, x) # [2.0, 4.0]
156+
157+ # Scalar→vector gradient (out-of-place)
158+ g(t) = [t^2, t^3]
159+ t = 2.0
160+ grad = finite_difference_gradient(g, t, Val(:central), eltype(t), Val(false))
161+ ```
162+
163+ # Notes
164+ - Forward differences: `O(n)` function evaluations, `O(h)` accuracy
165+ - Central differences: `O(2n)` function evaluations, `O(h²)` accuracy
166+ - Complex step: `O(n)` function evaluations, machine precision accuracy
130167"""
131168function finite_difference_gradient (
132- f,
133- x,
134- fdtype= Val (:central ),
135- returntype= eltype (x),
136- inplace= Val (true ),
137- fx= nothing ,
138- c1= nothing ,
139- c2= nothing ;
140- relstep= default_relstep (fdtype, eltype (x)),
141- absstep= relstep,
142- dir= true )
143-
169+ f,
170+ x,
171+ fdtype = Val (:central ),
172+ returntype = eltype (x),
173+ inplace = Val (true ),
174+ fx = nothing ,
175+ c1 = nothing ,
176+ c2 = nothing ;
177+ relstep = default_relstep (fdtype, eltype (x)),
178+ absstep = relstep,
179+ dir = true )
144180 inplace isa Type && (inplace = inplace ())
145181 if typeof (x) <: AbstractArray
146182 df = zero (returntype) .* x
@@ -162,7 +198,8 @@ function finite_difference_gradient(
162198 end
163199 end
164200 cache = GradientCache (df, x, fdtype, returntype, inplace)
165- finite_difference_gradient! (df, f, x, cache, relstep= relstep, absstep= absstep, dir= dir)
201+ finite_difference_gradient! (
202+ df, f, x, cache, relstep = relstep, absstep = absstep, dir = dir)
166203end
167204
168205"""
@@ -181,20 +218,19 @@ Gradients are either a vector->scalar map `f(x)`, or a scalar->vector map `f(fx,
181218Cache-less.
182219"""
183220function finite_difference_gradient! (
184- df,
185- f,
186- x,
187- fdtype= Val (:central ),
188- returntype= eltype (df),
189- inplace= Val (true ),
190- fx= nothing ,
191- c1= nothing ,
192- c2= nothing ;
193- relstep= default_relstep (fdtype, eltype (x)),
194- absstep= relstep)
195-
221+ df,
222+ f,
223+ x,
224+ fdtype = Val (:central ),
225+ returntype = eltype (df),
226+ inplace = Val (true ),
227+ fx = nothing ,
228+ c1 = nothing ,
229+ c2 = nothing ;
230+ relstep = default_relstep (fdtype, eltype (x)),
231+ absstep = relstep)
196232 cache = GradientCache (df, x, fdtype, returntype, inplace)
197- finite_difference_gradient! (df, f, x, cache, relstep= relstep, absstep= absstep)
233+ finite_difference_gradient! (df, f, x, cache, relstep = relstep, absstep = absstep)
198234end
199235
200236"""
@@ -212,32 +248,32 @@ Gradients are either a vector->scalar map `f(x)`, or a scalar->vector map `f(fx,
212248Cached.
213249"""
214250function finite_difference_gradient (
215- f,
216- x,
217- cache:: GradientCache{T1,T2,T3,T4,fdtype,returntype,inplace} ;
218- relstep= default_relstep (fdtype, eltype (x)),
219- absstep= relstep,
220- dir= true ) where {T1,T2,T3,T4,fdtype,returntype,inplace}
221-
251+ f,
252+ x,
253+ cache:: GradientCache{T1, T2, T3, T4, fdtype, returntype, inplace} ;
254+ relstep = default_relstep (fdtype, eltype (x)),
255+ absstep = relstep,
256+ dir = true ) where {T1, T2, T3, T4, fdtype, returntype, inplace}
222257 if typeof (x) <: AbstractArray
223258 df = zero (returntype) .* x
224259 else
225260 df = zero (cache. c1)
226261 end
227- finite_difference_gradient! (df, f, x, cache, relstep= relstep, absstep= absstep, dir= dir)
262+ finite_difference_gradient! (
263+ df, f, x, cache, relstep = relstep, absstep = absstep, dir = dir)
228264 df
229265end
230266
231267# vector of derivatives of a vector->scalar map by each component of a vector x
232268# this ignores the value of "inplace", because it doesn't make much sense
233269function finite_difference_gradient! (
234- df,
235- f,
236- x,
237- cache:: GradientCache{T1,T2,T3,T4,fdtype,returntype,inplace} ;
238- relstep= default_relstep (fdtype, eltype (x)),
239- absstep= relstep,
240- dir= true ) where {T1,T2,T3,T4,fdtype,returntype,inplace}
270+ df,
271+ f,
272+ x,
273+ cache:: GradientCache{T1, T2, T3, T4, fdtype, returntype, inplace} ;
274+ relstep = default_relstep (fdtype, eltype (x)),
275+ absstep = relstep,
276+ dir = true ) where {T1, T2, T3, T4, fdtype, returntype, inplace}
241277
242278 # NOTE: in this case epsilon is a vector, we need two arrays for epsilon and x1
243279 # c1 denotes x1, c2 is epsilon
@@ -247,11 +283,12 @@ function finite_difference_gradient!(
247283 end
248284 copyto! (c1, x)
249285 if fdtype == Val (:forward )
250- @inbounds for i ∈ eachindex (x)
286+ @inbounds for i in eachindex (x)
251287 if ArrayInterface. fast_scalar_indexing (c2)
252288 epsilon = ArrayInterface. allowed_getindex (c2, i) * dir
253289 else
254- epsilon = compute_epsilon (fdtype, one (eltype (x)), relstep, absstep, dir) * dir
290+ epsilon = compute_epsilon (fdtype, one (eltype (x)), relstep, absstep, dir) *
291+ dir
255292 end
256293 c1_old = ArrayInterface. allowed_getindex (c1, i)
257294 ArrayInterface. allowed_setindex! (c1, c1_old + epsilon, i)
@@ -278,11 +315,12 @@ function finite_difference_gradient!(
278315 end
279316 elseif fdtype == Val (:central )
280317 copyto! (c3, x)
281- @inbounds for i ∈ eachindex (x)
318+ @inbounds for i in eachindex (x)
282319 if ArrayInterface. fast_scalar_indexing (c2)
283320 epsilon = ArrayInterface. allowed_getindex (c2, i) * dir
284321 else
285- epsilon = compute_epsilon (fdtype, one (eltype (x)), relstep, absstep, dir) * dir
322+ epsilon = compute_epsilon (fdtype, one (eltype (x)), relstep, absstep, dir) *
323+ dir
286324 end
287325 c1_old = ArrayInterface. allowed_getindex (c1, i)
288326 ArrayInterface. allowed_setindex! (c1, c1_old + epsilon, i)
@@ -303,7 +341,7 @@ function finite_difference_gradient!(
303341 elseif fdtype == Val (:complex ) && returntype <: Real
304342 # we use c1 here to avoid typing issues with x
305343 epsilon_complex = eps (real (eltype (x)))
306- @inbounds for i ∈ eachindex (x)
344+ @inbounds for i in eachindex (x)
307345 c1_old = ArrayInterface. allowed_getindex (c1, i)
308346 ArrayInterface. allowed_setindex! (c1, c1_old + im * epsilon_complex, i)
309347 ArrayInterface. allowed_setindex! (df, imag (f (c1)) / epsilon_complex, i)
@@ -316,13 +354,13 @@ function finite_difference_gradient!(
316354end
317355
318356function finite_difference_gradient! (
319- df:: StridedVector{<:Number} ,
320- f,
321- x:: StridedVector{<:Number} ,
322- cache:: GradientCache{T1,T2,T3,T4,fdtype,returntype,inplace} ;
323- relstep= default_relstep (fdtype, eltype (x)),
324- absstep= relstep,
325- dir= true ) where {T1,T2,T3,T4,fdtype,returntype,inplace}
357+ df:: StridedVector{<:Number} ,
358+ f,
359+ x:: StridedVector{<:Number} ,
360+ cache:: GradientCache{T1, T2, T3, T4, fdtype, returntype, inplace} ;
361+ relstep = default_relstep (fdtype, eltype (x)),
362+ absstep = relstep,
363+ dir = true ) where {T1, T2, T3, T4, fdtype, returntype, inplace}
326364
327365 # c1 is x1 if we need a complex copy of x, otherwise Nothing
328366 # c2 is Nothing
@@ -334,7 +372,7 @@ function finite_difference_gradient!(
334372 end
335373 copyto! (c3, x)
336374 if fdtype == Val (:forward )
337- for i ∈ eachindex (x)
375+ for i in eachindex (x)
338376 epsilon = compute_epsilon (fdtype, x[i], relstep, absstep, dir)
339377 x_old = x[i]
340378 if typeof (fx) != Nothing
@@ -371,7 +409,7 @@ function finite_difference_gradient!(
371409 end
372410 end
373411 elseif fdtype == Val (:central )
374- @inbounds for i ∈ eachindex (x)
412+ @inbounds for i in eachindex (x)
375413 epsilon = compute_epsilon (fdtype, x[i], relstep, absstep, dir)
376414 x_old = x[i]
377415 c3[i] += epsilon
@@ -397,11 +435,12 @@ function finite_difference_gradient!(
397435 df[i] -= im * imag (dfi / (2 * im * epsilon))
398436 end
399437 end
400- elseif fdtype == Val (:complex ) && returntype <: Real && eltype (df) <: Real && eltype (x) <: Real
438+ elseif fdtype == Val (:complex ) && returntype <: Real && eltype (df) <: Real &&
439+ eltype (x) <: Real
401440 copyto! (c1, x)
402441 epsilon_complex = eps (real (eltype (x)))
403442 # we use c1 here to avoid typing issues with x
404- @inbounds for i ∈ eachindex (x)
443+ @inbounds for i in eachindex (x)
405444 c1_old = c1[i]
406445 c1[i] += im * epsilon_complex
407446 df[i] = imag (f (c1)) / epsilon_complex
@@ -416,13 +455,13 @@ end
416455# vector of derivatives of a scalar->vector map
417456# this is effectively a vector of partial derivatives, but we still call it a gradient
418457function finite_difference_gradient! (
419- df,
420- f,
421- x:: Number ,
422- cache:: GradientCache{T1,T2,T3,T4,fdtype,returntype,inplace} ;
423- relstep= default_relstep (fdtype, eltype (x)),
424- absstep= relstep,
425- dir= true ) where {T1,T2,T3,T4,fdtype,returntype,inplace}
458+ df,
459+ f,
460+ x:: Number ,
461+ cache:: GradientCache{T1, T2, T3, T4, fdtype, returntype, inplace} ;
462+ relstep = default_relstep (fdtype, eltype (x)),
463+ absstep = relstep,
464+ dir = true ) where {T1, T2, T3, T4, fdtype, returntype, inplace}
426465
427466 # NOTE: in this case epsilon is a scalar, we need two arrays for fx1 and fx2
428467 # c1 denotes fx1, c2 is fx2, sizes guaranteed by the cache constructor
0 commit comments