1
+ module ForwardDiffStaticArraysExt
2
+
3
+ using ForwardDiff
4
+ using ForwardDiff. DiffResults: DiffResults, DiffResult, ImmutableDiffResult, MutableDiffResult
5
+ using ForwardDiff. LinearAlgebra
6
+ import ForwardDiff: Dual, Chunk, value, extract_jacobian!, extract_value!, extract_gradient!, extract_jacobian!,
7
+ GradientConfig, JacobianConfig, HessianConfig, vector_mode_gradient, vector_mode_gradient!,
8
+ Tag, valtype, partials, gradient, gradient!, jacobian, jacobian!, hessian, hessian!, vector_mode_jacobian,
9
+ vector_mode_jacobian!
10
+ 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* ForwardDiff. _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