@@ -140,15 +140,63 @@ end
140140
141141# # Array to scalar
142142
143- arr_to_num (x:: AbstractArray ):: Number = sum (sin, x)
143+ const DEFAULT_α = 4
144+ const DEFAULT_β = 6
144145
145- arr_to_num_gradient (x) = cos .(x)
146- arr_to_num_hvp (x, v) = - sin .(x) .* v
146+ arr_to_num_aux_linalg (x; α, β) = sum (vec (x .^ α) .* transpose (vec (x .^ β)))
147+
148+ function arr_to_num_aux_no_linalg (x; α, β)
149+ n = length (x)
150+ s = zero (eltype (x))
151+ for i in 1 : n, j in 1 : n
152+ s += x[i]^ α * x[j]^ β
153+ end
154+ return s
155+ end
156+
157+ function arr_to_num_aux_gradient (x; α, β)
158+ x = Array (x) # GPU arrays don't like indexing
159+ g = similar (x)
160+ for k in eachindex (g, x)
161+ g[k] = (
162+ α * x[k]^ (α - 1 ) * sum (x[j]^ β for j in eachindex (x) if j != k) +
163+ β * x[k]^ (β - 1 ) * sum (x[i]^ α for i in eachindex (x) if i != k) +
164+ (α + β) * x[k]^ (α + β - 1 )
165+ )
166+ end
167+ return g
168+ end
169+
170+ function arr_to_num_aux_hessian (x; α, β)
171+ x = Array (x) # GPU arrays don't like indexing
172+ H = similar (x, length (x), length (x))
173+ for k in axes (H, 1 ), l in axes (H, 2 )
174+ if k == l
175+ H[k, k] = (
176+ α * (α - 1 ) * x[k]^ (α - 2 ) * sum (x[j]^ β for j in eachindex (x) if j != k) +
177+ β * (β - 1 ) * x[k]^ (β - 2 ) * sum (x[i]^ α for i in eachindex (x) if i != k) +
178+ (α + β) * (α + β - 1 ) * x[k]^ (α + β - 2 )
179+ )
180+ else
181+ H[k, l] = α * β * (x[k]^ (α - 1 ) * x[l]^ (β - 1 ) + x[k]^ (β - 1 ) * x[l]^ (α - 1 ))
182+ end
183+ end
184+ return H
185+ end
186+
187+ arr_to_num_linalg (x:: AbstractArray ):: Number =
188+ arr_to_num_aux_linalg (x; α= DEFAULT_α, β= DEFAULT_β)
189+ arr_to_num_no_linalg (x:: AbstractArray ):: Number =
190+ arr_to_num_aux_no_linalg (x; α= DEFAULT_α, β= DEFAULT_β)
191+
192+ arr_to_num_gradient (x) = arr_to_num_aux_gradient (x; α= DEFAULT_α, β= DEFAULT_β)
193+ arr_to_num_hessian (x) = arr_to_num_aux_hessian (x; α= DEFAULT_α, β= DEFAULT_β)
147194arr_to_num_pushforward (x, dx) = dot (arr_to_num_gradient (x), dx)
148195arr_to_num_pullback (x, dy) = arr_to_num_gradient (x) .* dy
149- arr_to_num_hessian (x ) = Matrix ( Diagonal ( - sin .( vec (x)) ))
196+ arr_to_num_hvp (x, v ) = reshape ( arr_to_num_hessian (x) * vec (v), size (x ))
150197
151- function arr_to_num_scenarios_onearg (x:: AbstractArray )
198+ function arr_to_num_scenarios_onearg (x:: AbstractArray ; linalg= true )
199+ arr_to_num = linalg ? arr_to_num_linalg : arr_to_num_no_linalg
152200 # pushforward stays out of place
153201 scens = AbstractScenario[]
154202 for place in (:outofplace , :inplace )
@@ -448,24 +496,24 @@ const IMAT = Matrix((1:2) .* transpose(1:3))
448496
449497Create a vector of [`AbstractScenario`](@ref)s with standard array types.
450498"""
451- function default_scenarios ()
499+ function default_scenarios (; linalg = true )
452500 return vcat (
453501 # one argument
454- num_to_num_scenarios_onearg (randn ()),
455- num_to_arr_scenarios_onearg (randn (), IVEC),
456- num_to_arr_scenarios_onearg (randn (), IMAT),
457- arr_to_num_scenarios_onearg (randn (6 )),
458- arr_to_num_scenarios_onearg (randn (2 , 3 )),
459- vec_to_vec_scenarios_onearg (randn (6 )),
460- vec_to_mat_scenarios_onearg (randn (6 )),
461- mat_to_vec_scenarios_onearg (randn (2 , 3 )),
462- mat_to_mat_scenarios_onearg (randn (2 , 3 )),
502+ num_to_num_scenarios_onearg (rand ()),
503+ num_to_arr_scenarios_onearg (rand (), IVEC),
504+ num_to_arr_scenarios_onearg (rand (), IMAT),
505+ arr_to_num_scenarios_onearg (rand (6 ); linalg ),
506+ arr_to_num_scenarios_onearg (rand (2 , 3 ); linalg ),
507+ vec_to_vec_scenarios_onearg (rand (6 )),
508+ vec_to_mat_scenarios_onearg (rand (6 )),
509+ mat_to_vec_scenarios_onearg (rand (2 , 3 )),
510+ mat_to_mat_scenarios_onearg (rand (2 , 3 )),
463511 # two arguments
464- num_to_arr_scenarios_twoarg (randn (), IVEC),
465- num_to_arr_scenarios_twoarg (randn (), IMAT),
466- vec_to_vec_scenarios_twoarg (randn (6 )),
467- vec_to_mat_scenarios_twoarg (randn (6 )),
468- mat_to_vec_scenarios_twoarg (randn (2 , 3 )),
469- mat_to_mat_scenarios_twoarg (randn (2 , 3 )),
512+ num_to_arr_scenarios_twoarg (rand (), IVEC),
513+ num_to_arr_scenarios_twoarg (rand (), IMAT),
514+ vec_to_vec_scenarios_twoarg (rand (6 )),
515+ vec_to_mat_scenarios_twoarg (rand (6 )),
516+ mat_to_vec_scenarios_twoarg (rand (2 , 3 )),
517+ mat_to_mat_scenarios_twoarg (rand (2 , 3 )),
470518 )
471519end
0 commit comments