1+ module ForwardDiffStaticArraysExt
2+
3+ using ForwardDiff
4+ import DiffResults
5+ using LinearAlgebra
6+ import DiffResults: DiffResult, ImmutableDiffResult, MutableDiffResult
7+ import ForwardDiff: Dual, Chunk, extract_jacobian!, extract_value!, extract_gradient!, extract_jacobian!,
8+ GradientConfig, JacobianConfig, HessianConfig, ImmutableDiffResult
9+
10+ isdefined (Base, :get_extension ) ? (using StaticArrays) : (using .. StaticArrays)
11+
12+ @generated function dualize (:: Type{T} , x:: StaticArray ) where T
13+ N = length (x)
14+ dx = Expr (:tuple , [:(Dual {T} (x[$ i], chunk, Val {$i} ())) for i in 1 : N]. .. )
15+ V = StaticArrays. similar_type (x, Dual{T,eltype (x),N})
16+ return quote
17+ chunk = Chunk {$N} ()
18+ $ (Expr (:meta , :inline ))
19+ return $ V ($ (dx))
20+ end
21+ end
22+
23+ @inline static_dual_eval (:: Type{T} , f, x:: StaticArray ) where T = f (dualize (T, x))
24+
25+ @inline ForwardDiff. gradient (f, x:: StaticArray ) = vector_mode_gradient (f, x)
26+ @inline ForwardDiff. gradient (f, x:: StaticArray , cfg:: GradientConfig ) = gradient (f, x)
27+ @inline ForwardDiff. gradient (f, x:: StaticArray , cfg:: GradientConfig , :: Val ) = gradient (f, x)
28+
29+ @inline ForwardDiff. gradient! (result:: Union{AbstractArray,DiffResult} , f, x:: StaticArray ) = vector_mode_gradient! (result, f, x)
30+ @inline ForwardDiff. gradient! (result:: Union{AbstractArray,DiffResult} , f, x:: StaticArray , cfg:: GradientConfig ) = gradient! (result, f, x)
31+ @inline ForwardDiff. gradient! (result:: Union{AbstractArray,DiffResult} , f, x:: StaticArray , cfg:: GradientConfig , :: Val ) = gradient! (result, f, x)
32+
33+ @generated function extract_gradient (:: Type{T} , y:: Real , x:: S ) where {T,S<: StaticArray }
34+ result = Expr (:tuple , [:(partials (T, y, $ i)) for i in 1 : length (x)]. .. )
35+ return quote
36+ $ (Expr (:meta , :inline ))
37+ V = StaticArrays. similar_type (S, valtype ($ y))
38+ return V ($ result)
39+ end
40+ end
41+
42+ @inline function ForwardDiff. vector_mode_gradient (f, x:: StaticArray )
43+ T = typeof (Tag (f, eltype (x)))
44+ return extract_gradient (T, static_dual_eval (T, f, x), x)
45+ end
46+
47+ @inline function ForwardDiff. vector_mode_gradient! (result, f, x:: StaticArray )
48+ T = typeof (Tag (f, eltype (x)))
49+ return extract_gradient! (T, result, static_dual_eval (T, f, x))
50+ end
51+
52+ @inline ForwardDiff. jacobian (f, x:: StaticArray ) = vector_mode_jacobian (f, x)
53+ @inline ForwardDiff. jacobian (f, x:: StaticArray , cfg:: JacobianConfig ) = jacobian (f, x)
54+ @inline ForwardDiff. jacobian (f, x:: StaticArray , cfg:: JacobianConfig , :: Val ) = jacobian (f, x)
55+
56+ @inline ForwardDiff. jacobian! (result:: Union{AbstractArray,DiffResult} , f, x:: StaticArray ) = vector_mode_jacobian! (result, f, x)
57+ @inline ForwardDiff. jacobian! (result:: Union{AbstractArray,DiffResult} , f, x:: StaticArray , cfg:: JacobianConfig ) = jacobian! (result, f, x)
58+ @inline ForwardDiff. jacobian! (result:: Union{AbstractArray,DiffResult} , f, x:: StaticArray , cfg:: JacobianConfig , :: Val ) = jacobian! (result, f, x)
59+
60+ @generated function extract_jacobian (:: Type{T} , ydual:: StaticArray , x:: S ) where {T,S<: StaticArray }
61+ M, N = length (ydual), length (x)
62+ result = Expr (:tuple , [:(partials (T, ydual[$ i], $ j)) for i in 1 : M, j in 1 : N]. .. )
63+ return quote
64+ $ (Expr (:meta , :inline ))
65+ V = StaticArrays. similar_type (S, valtype (eltype ($ ydual)), Size ($ M, $ N))
66+ return V ($ result)
67+ end
68+ end
69+
70+ function extract_jacobian (:: Type{T} , ydual:: AbstractArray , x:: StaticArray ) where T
71+ result = similar (ydual, valtype (eltype (ydual)), length (ydual), length (x))
72+ return extract_jacobian! (T, result, ydual, length (x))
73+ end
74+
75+ @inline function ForwardDiff. vector_mode_jacobian (f, x:: StaticArray )
76+ T = typeof (Tag (f, eltype (x)))
77+ return extract_jacobian (T, static_dual_eval (T, f, x), x)
78+ end
79+
80+ @inline function ForwardDiff. vector_mode_jacobian! (result, f, x:: StaticArray )
81+ T = typeof (Tag (f, eltype (x)))
82+ ydual = static_dual_eval (T, f, x)
83+ result = extract_jacobian! (T, result, ydual, length (x))
84+ result = extract_value! (T, result, ydual)
85+ return result
86+ end
87+
88+ @inline function ForwardDiff. vector_mode_jacobian! (result:: ImmutableDiffResult , f, x:: StaticArray )
89+ T = typeof (Tag (f, eltype (x)))
90+ ydual = static_dual_eval (T, f, x)
91+ result = DiffResults. jacobian! (result, extract_jacobian (T, ydual, x))
92+ result = DiffResults. value! (d -> value (T,d), result, ydual)
93+ return result
94+ end
95+
96+ ForwardDiff. hessian (f, x:: StaticArray ) = jacobian (y -> gradient (f, y), x)
97+ ForwardDiff. hessian (f, x:: StaticArray , cfg:: HessianConfig ) = hessian (f, x)
98+ ForwardDiff. hessian (f, x:: StaticArray , cfg:: HessianConfig , :: Val ) = hessian (f, x)
99+
100+ ForwardDiff. hessian! (result:: AbstractArray , f, x:: StaticArray ) = jacobian! (result, y -> gradient (f, y), x)
101+
102+ ForwardDiff. hessian! (result:: MutableDiffResult , f, x:: StaticArray ) = hessian! (result, f, x, HessianConfig (f, result, x))
103+
104+ ForwardDiff. hessian! (result:: ImmutableDiffResult , f, x:: StaticArray , cfg:: HessianConfig ) = hessian! (result, f, x)
105+ ForwardDiff. hessian! (result:: ImmutableDiffResult , f, x:: StaticArray , cfg:: HessianConfig , :: Val ) = hessian! (result, f, x)
106+
107+ function ForwardDiff. hessian! (result:: ImmutableDiffResult , f, x:: StaticArray )
108+ T = typeof (Tag (f, eltype (x)))
109+ d1 = dualize (T, x)
110+ d2 = dualize (T, d1)
111+ fd2 = f (d2)
112+ val = value (T,value (T,fd2))
113+ grad = extract_gradient (T,value (T,fd2), x)
114+ hess = extract_jacobian (T,partials (T,fd2), x)
115+ result = DiffResults. hessian! (result, hess)
116+ result = DiffResults. gradient! (result, grad)
117+ result = DiffResults. value! (result, val)
118+ return result
119+ end
120+
121+ function LinearAlgebra. eigvals (A:: Symmetric{<:Dual{Tg,T,N}, <:StaticArrays.StaticMatrix} ) where {Tg,T<: Real ,N}
122+ λ,Q = eigen (Symmetric (value .(parent (A))))
123+ parts = ntuple (j -> diag (Q' * getindex .(partials .(A), j) * Q), N)
124+ Dual {Tg} .(λ, tuple .(parts... ))
125+ end
126+
127+ function LinearAlgebra. eigen (A:: Symmetric{<:Dual{Tg,T,N}, <:StaticArrays.StaticMatrix} ) where {Tg,T<: Real ,N}
128+ λ = eigvals (A)
129+ _,Q = eigen (Symmetric (value .(parent (A))))
130+ parts = ntuple (j -> Q* _lyap_div! (Q' * getindex .(partials .(A), j) * Q - Diagonal (getindex .(partials .(λ), j)), value .(λ)), N)
131+ Eigen (λ,Dual {Tg} .(Q, tuple .(parts... )))
132+ end
133+
134+ end
0 commit comments