|
1 | 1 |
|
2 |
| -if VERSION ≥ v"1.4" |
3 |
| - function frule((_, ẋ, ṗ), ::typeof(evalpoly), x, p::Union{Tuple,AbstractVector}) |
4 |
| - Δx, Δp = ẋ, unthunk(ṗ) |
5 |
| - N = length(p) |
6 |
| - @inbounds y = p[N] |
7 |
| - Δy = Δp[N] |
8 |
| - @inbounds for i in (N - 1):-1:1 |
9 |
| - Δy = muladd(Δx, y, muladd(x, Δy, Δp[i])) |
10 |
| - y = muladd(x, y, p[i]) |
11 |
| - end |
12 |
| - return y, Δy |
| 2 | +function frule((_, ẋ, ṗ), ::typeof(evalpoly), x, p::Union{Tuple,AbstractVector}) |
| 3 | + Δx, Δp = ẋ, unthunk(ṗ) |
| 4 | + N = length(p) |
| 5 | + @inbounds y = p[N] |
| 6 | + Δy = Δp[N] |
| 7 | + @inbounds for i in (N - 1):-1:1 |
| 8 | + Δy = muladd(Δx, y, muladd(x, Δy, Δp[i])) |
| 9 | + y = muladd(x, y, p[i]) |
13 | 10 | end
|
| 11 | + return y, Δy |
| 12 | +end |
14 | 13 |
|
15 |
| - function rrule(::typeof(evalpoly), x, p::Union{Tuple,AbstractVector}) |
16 |
| - y, ys = _evalpoly_intermediates(x, p) |
17 |
| - project_x = ProjectTo(x) |
18 |
| - project_p = p isa Tuple ? identity : ProjectTo(p) |
19 |
| - function evalpoly_pullback(Δy) |
20 |
| - ∂x, ∂p = _evalpoly_back(x, p, ys, Δy) |
21 |
| - return NoTangent(), project_x(∂x), project_p(∂p) |
22 |
| - end |
23 |
| - return y, evalpoly_pullback |
| 14 | +function rrule(::typeof(evalpoly), x, p::Union{Tuple,AbstractVector}) |
| 15 | + y, ys = _evalpoly_intermediates(x, p) |
| 16 | + project_x = ProjectTo(x) |
| 17 | + project_p = p isa Tuple ? identity : ProjectTo(p) |
| 18 | + function evalpoly_pullback(Δy) |
| 19 | + ∂x, ∂p = _evalpoly_back(x, p, ys, Δy) |
| 20 | + return NoTangent(), project_x(∂x), project_p(∂p) |
24 | 21 | end
|
| 22 | + return y, evalpoly_pullback |
| 23 | +end |
25 | 24 |
|
26 |
| - function rrule(::typeof(evalpoly), x, p::Vector{<:Matrix}) # does not type infer with ProjectTo |
27 |
| - y, ys = _evalpoly_intermediates(x, p) |
28 |
| - function evalpoly_pullback(Δy) |
29 |
| - ∂x, ∂p = _evalpoly_back(x, p, ys, Δy) |
30 |
| - return NoTangent(), ∂x, ∂p |
31 |
| - end |
32 |
| - return y, evalpoly_pullback |
| 25 | +function rrule(::typeof(evalpoly), x, p::Vector{<:Matrix}) # does not type infer with ProjectTo |
| 26 | + y, ys = _evalpoly_intermediates(x, p) |
| 27 | + function evalpoly_pullback(Δy) |
| 28 | + ∂x, ∂p = _evalpoly_back(x, p, ys, Δy) |
| 29 | + return NoTangent(), ∂x, ∂p |
33 | 30 | end
|
| 31 | + return y, evalpoly_pullback |
| 32 | +end |
34 | 33 |
|
35 |
| - # evalpoly but storing intermediates |
36 |
| - function _evalpoly_intermediates(x, p::Tuple) |
37 |
| - return if @generated |
38 |
| - N = length(p.parameters) |
39 |
| - exs = [] |
40 |
| - vars = [] |
41 |
| - ex = :(p[$N]) |
42 |
| - for i in 1:(N - 1) |
43 |
| - yi = Symbol("y", i) |
44 |
| - push!(vars, yi) |
45 |
| - push!(exs, :($yi = $ex)) |
46 |
| - ex = :(muladd(x, $yi, p[$(N - i)])) |
47 |
| - end |
48 |
| - push!(exs, :(y = $ex)) |
49 |
| - Expr(:block, exs..., :(y, ($(vars...),))) |
50 |
| - else |
51 |
| - _evalpoly_intermediates_fallback(x, p) |
| 34 | +# evalpoly but storing intermediates |
| 35 | +function _evalpoly_intermediates(x, p::Tuple) |
| 36 | + return if @generated |
| 37 | + N = length(p.parameters) |
| 38 | + exs = [] |
| 39 | + vars = [] |
| 40 | + ex = :(p[$N]) |
| 41 | + for i in 1:(N - 1) |
| 42 | + yi = Symbol("y", i) |
| 43 | + push!(vars, yi) |
| 44 | + push!(exs, :($yi = $ex)) |
| 45 | + ex = :(muladd(x, $yi, p[$(N - i)])) |
52 | 46 | end
|
| 47 | + push!(exs, :(y = $ex)) |
| 48 | + Expr(:block, exs..., :(y, ($(vars...),))) |
| 49 | + else |
| 50 | + _evalpoly_intermediates_fallback(x, p) |
53 | 51 | end
|
54 |
| - function _evalpoly_intermediates_fallback(x, p::Tuple) |
55 |
| - N = length(p) |
56 |
| - y = p[N] |
57 |
| - ys = (y, ntuple(N - 2) do i |
58 |
| - return y = muladd(x, y, p[N - i]) |
59 |
| - end...) |
60 |
| - y = muladd(x, y, p[1]) |
61 |
| - return y, ys |
62 |
| - end |
63 |
| - function _evalpoly_intermediates(x, p) |
64 |
| - N = length(p) |
65 |
| - @inbounds yn = one(x) * p[N] |
66 |
| - ys = similar(p, typeof(yn), N - 1) |
67 |
| - @inbounds ys[1] = yn |
68 |
| - @inbounds for i in 2:(N - 1) |
69 |
| - ys[i] = muladd(x, ys[i - 1], p[N - i + 1]) |
70 |
| - end |
71 |
| - @inbounds y = muladd(x, ys[N - 1], p[1]) |
72 |
| - return y, ys |
| 52 | +end |
| 53 | +function _evalpoly_intermediates_fallback(x, p::Tuple) |
| 54 | + N = length(p) |
| 55 | + y = p[N] |
| 56 | + ys = (y, ntuple(N - 2) do i |
| 57 | + return y = muladd(x, y, p[N - i]) |
| 58 | + end...) |
| 59 | + y = muladd(x, y, p[1]) |
| 60 | + return y, ys |
| 61 | +end |
| 62 | +function _evalpoly_intermediates(x, p) |
| 63 | + N = length(p) |
| 64 | + @inbounds yn = one(x) * p[N] |
| 65 | + ys = similar(p, typeof(yn), N - 1) |
| 66 | + @inbounds ys[1] = yn |
| 67 | + @inbounds for i in 2:(N - 1) |
| 68 | + ys[i] = muladd(x, ys[i - 1], p[N - i + 1]) |
73 | 69 | end
|
| 70 | + @inbounds y = muladd(x, ys[N - 1], p[1]) |
| 71 | + return y, ys |
| 72 | +end |
74 | 73 |
|
75 |
| - # TODO: Handle following cases |
76 |
| - # 1) x is a UniformScaling, pᵢ is a matrix |
77 |
| - # 2) x is a matrix, pᵢ is a UniformScaling |
78 |
| - @inline _evalpoly_backx(x, yi, ∂yi) = ∂yi * yi' |
79 |
| - @inline _evalpoly_backx(x, yi, ∂x, ∂yi) = muladd(∂yi, yi', ∂x) |
80 |
| - @inline _evalpoly_backx(x::Number, yi, ∂yi) = conj(dot(∂yi, yi)) |
81 |
| - @inline _evalpoly_backx(x::Number, yi, ∂x, ∂yi) = _evalpoly_backx(x, yi, ∂yi) + ∂x |
| 74 | +# TODO: Handle following cases |
| 75 | +# 1) x is a UniformScaling, pᵢ is a matrix |
| 76 | +# 2) x is a matrix, pᵢ is a UniformScaling |
| 77 | +@inline _evalpoly_backx(x, yi, ∂yi) = ∂yi * yi' |
| 78 | +@inline _evalpoly_backx(x, yi, ∂x, ∂yi) = muladd(∂yi, yi', ∂x) |
| 79 | +@inline _evalpoly_backx(x::Number, yi, ∂yi) = conj(dot(∂yi, yi)) |
| 80 | +@inline _evalpoly_backx(x::Number, yi, ∂x, ∂yi) = _evalpoly_backx(x, yi, ∂yi) + ∂x |
82 | 81 |
|
83 |
| - @inline _evalpoly_backp(pi, ∂yi) = ∂yi |
| 82 | +@inline _evalpoly_backp(pi, ∂yi) = ∂yi |
84 | 83 |
|
85 |
| - function _evalpoly_back(x, p::Tuple, ys, Δy) |
86 |
| - return if @generated |
87 |
| - exs = [] |
88 |
| - vars = [] |
89 |
| - N = length(p.parameters) |
90 |
| - for i in 2:(N - 1) |
91 |
| - ∂pi = Symbol("∂p", i) |
92 |
| - push!(vars, ∂pi) |
93 |
| - push!(exs, :(∂x = _evalpoly_backx(x, ys[$(N - i)], ∂x, ∂yi))) |
94 |
| - push!(exs, :($∂pi = _evalpoly_backp(p[$i], ∂yi))) |
95 |
| - push!(exs, :(∂yi = x′ * ∂yi)) |
96 |
| - end |
97 |
| - push!(vars, :(_evalpoly_backp(p[$N], ∂yi))) # ∂pN |
98 |
| - Expr( |
99 |
| - :block, |
100 |
| - :(x′ = x'), |
101 |
| - :(∂yi = Δy), |
102 |
| - :(∂p1 = _evalpoly_backp(p[1], ∂yi)), |
103 |
| - :(∂x = _evalpoly_backx(x, ys[$(N - 1)], ∂yi)), |
104 |
| - :(∂yi = x′ * ∂yi), |
105 |
| - exs..., |
106 |
| - :(∂p = (∂p1, $(vars...))), |
107 |
| - :(∂x, Tangent{typeof(p),typeof(∂p)}(∂p)), |
108 |
| - ) |
109 |
| - else |
110 |
| - _evalpoly_back_fallback(x, p, ys, Δy) |
| 84 | +function _evalpoly_back(x, p::Tuple, ys, Δy) |
| 85 | + return if @generated |
| 86 | + exs = [] |
| 87 | + vars = [] |
| 88 | + N = length(p.parameters) |
| 89 | + for i in 2:(N - 1) |
| 90 | + ∂pi = Symbol("∂p", i) |
| 91 | + push!(vars, ∂pi) |
| 92 | + push!(exs, :(∂x = _evalpoly_backx(x, ys[$(N - i)], ∂x, ∂yi))) |
| 93 | + push!(exs, :($∂pi = _evalpoly_backp(p[$i], ∂yi))) |
| 94 | + push!(exs, :(∂yi = x′ * ∂yi)) |
111 | 95 | end
|
| 96 | + push!(vars, :(_evalpoly_backp(p[$N], ∂yi))) # ∂pN |
| 97 | + Expr( |
| 98 | + :block, |
| 99 | + :(x′ = x'), |
| 100 | + :(∂yi = Δy), |
| 101 | + :(∂p1 = _evalpoly_backp(p[1], ∂yi)), |
| 102 | + :(∂x = _evalpoly_backx(x, ys[$(N - 1)], ∂yi)), |
| 103 | + :(∂yi = x′ * ∂yi), |
| 104 | + exs..., |
| 105 | + :(∂p = (∂p1, $(vars...))), |
| 106 | + :(∂x, Tangent{typeof(p),typeof(∂p)}(∂p)), |
| 107 | + ) |
| 108 | + else |
| 109 | + _evalpoly_back_fallback(x, p, ys, Δy) |
112 | 110 | end
|
113 |
| - function _evalpoly_back_fallback(x, p::Tuple, ys, Δy) |
114 |
| - x′ = x' |
115 |
| - ∂yi = unthunk(Δy) |
116 |
| - N = length(p) |
117 |
| - ∂p1 = _evalpoly_backp(p[1], ∂yi) |
| 111 | +end |
| 112 | +function _evalpoly_back_fallback(x, p::Tuple, ys, Δy) |
| 113 | + x′ = x' |
| 114 | + ∂yi = unthunk(Δy) |
| 115 | + N = length(p) |
| 116 | + ∂p1 = _evalpoly_backp(p[1], ∂yi) |
| 117 | + ∂x = _evalpoly_backx(x, ys[N - 1], ∂yi) |
| 118 | + ∂yi = x′ * ∂yi |
| 119 | + ∂p = ( |
| 120 | + ∂p1, |
| 121 | + ntuple(N - 2) do i |
| 122 | + ∂x = _evalpoly_backx(x, ys[N-i-1], ∂x, ∂yi) |
| 123 | + ∂pi = _evalpoly_backp(p[i+1], ∂yi) |
| 124 | + ∂yi = x′ * ∂yi |
| 125 | + return ∂pi |
| 126 | + end..., |
| 127 | + _evalpoly_backp(p[N], ∂yi), # ∂pN |
| 128 | + ) |
| 129 | + return ∂x, Tangent{typeof(p),typeof(∂p)}(∂p) |
| 130 | +end |
| 131 | +function _evalpoly_back(x, p, ys, Δy) |
| 132 | + x′ = x' |
| 133 | + ∂yi = one(x′) * Δy |
| 134 | + N = length(p) |
| 135 | + @inbounds ∂p1 = _evalpoly_backp(p[1], ∂yi) |
| 136 | + ∂p = similar(p, typeof(∂p1)) |
| 137 | + @inbounds begin |
118 | 138 | ∂x = _evalpoly_backx(x, ys[N - 1], ∂yi)
|
119 | 139 | ∂yi = x′ * ∂yi
|
120 |
| - ∂p = ( |
121 |
| - ∂p1, |
122 |
| - ntuple(N - 2) do i |
123 |
| - ∂x = _evalpoly_backx(x, ys[N-i-1], ∂x, ∂yi) |
124 |
| - ∂pi = _evalpoly_backp(p[i+1], ∂yi) |
125 |
| - ∂yi = x′ * ∂yi |
126 |
| - return ∂pi |
127 |
| - end..., |
128 |
| - _evalpoly_backp(p[N], ∂yi), # ∂pN |
129 |
| - ) |
130 |
| - return ∂x, Tangent{typeof(p),typeof(∂p)}(∂p) |
131 |
| - end |
132 |
| - function _evalpoly_back(x, p, ys, Δy) |
133 |
| - x′ = x' |
134 |
| - ∂yi = one(x′) * Δy |
135 |
| - N = length(p) |
136 |
| - @inbounds ∂p1 = _evalpoly_backp(p[1], ∂yi) |
137 |
| - ∂p = similar(p, typeof(∂p1)) |
138 |
| - @inbounds begin |
139 |
| - ∂x = _evalpoly_backx(x, ys[N - 1], ∂yi) |
| 140 | + ∂p[1] = ∂p1 |
| 141 | + for i in 2:(N - 1) |
| 142 | + ∂x = _evalpoly_backx(x, ys[N - i], ∂x, ∂yi) |
| 143 | + ∂p[i] = _evalpoly_backp(p[i], ∂yi) |
140 | 144 | ∂yi = x′ * ∂yi
|
141 |
| - ∂p[1] = ∂p1 |
142 |
| - for i in 2:(N - 1) |
143 |
| - ∂x = _evalpoly_backx(x, ys[N - i], ∂x, ∂yi) |
144 |
| - ∂p[i] = _evalpoly_backp(p[i], ∂yi) |
145 |
| - ∂yi = x′ * ∂yi |
146 |
| - end |
147 |
| - ∂p[N] = _evalpoly_backp(p[N], ∂yi) |
148 | 145 | end
|
149 |
| - return ∂x, ∂p |
| 146 | + ∂p[N] = _evalpoly_backp(p[N], ∂yi) |
150 | 147 | end
|
| 148 | + return ∂x, ∂p |
151 | 149 | end
|
0 commit comments