/
diagnostic_functions.jl
334 lines (273 loc) · 10.8 KB
/
diagnostic_functions.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
"""
est_dΣdθ(glo, kernel_hyperparameters; return_est=true, return_anal=false, return_dif=false, return_bool=false, dif=1e-6, print_stuff=true)
Estimate the covariance derivatives of `glo` (a GLOM model) at
`kernel_hyperparameters`with forward differences
# Keyword Arguments
- `return_est::Bool=true`: Add the numerically estimated dΣdθs to the output
- `return_anal::Bool=false`: Add the analytical dΣdθs to the output
- `return_dif::Bool=false`: Add the difference between the dΣdθs to the output
- `return_bool::Bool=false`: Return just a similarity Boolean
- `dif::Real=1e-6`: The step size used in the forward difference method
- `print_stuff::Bool=true`: Whether to print things
# Output
- If `return_bool==true`, returns a Boolean for whether the analytical and
numerical estimates for the dΣdθs are approximately the same.
- Else returns a vector with some combination of the numerically estimated
dΣdθs, analytical dΣdθs, and differences between them
"""
function est_dΣdθ(glo::GLO, kernel_hyperparameters::Vector{T}; return_est::Bool=true, return_anal::Bool=false, return_dif::Bool=false, return_bool::Bool=false, dif::Real=1e-6, print_stuff::Bool=true, kwargs...) where {T<:Real}
total_hyperparameters = append!(collect(Iterators.flatten(glo.a)), kernel_hyperparameters)
if print_stuff
println()
println("est_dΣdθ: Check that our analytical dΣdθ are close to numerical estimates")
println("hyperparameters: ", total_hyperparameters)
end
x = glo.x_obs
return_vec = []
coeff_orders = coefficient_orders(glo.n_out, glo.n_dif, a=glo.a)
# construct estimated dΣdθs
if return_est || return_dif || return_bool
val = covariance(glo, total_hyperparameters)
est_dΣdθs = zeros(length(total_hyperparameters), glo.n_out * length(x), glo.n_out * length(x))
for i in 1:length(total_hyperparameters)
if total_hyperparameters[i]!=0
hold = copy(total_hyperparameters)
hold[i] += dif
est_dΣdθs[i, :, :] = (covariance(glo, hold) - val) / dif
end
end
if return_est; append!(return_vec, [est_dΣdθs]) end
end
# construct analytical dΣdθs
if return_anal || return_dif || return_bool
anal_dΣdθs = zeros(length(total_hyperparameters), glo.n_out * length(x), glo.n_out * length(x))
for i in 1:length(total_hyperparameters)
anal_dΣdθs[i, :, :] = covariance(glo, total_hyperparameters; dΣdθs_total=[i], kwargs...)
end
if return_anal; append!(return_vec, [anal_dΣdθs]) end
end
if return_dif || return_bool
difs = est_dΣdθs - anal_dΣdθs
append!(return_vec, [difs])
end
# find whether the analytical and estimated dΣdθs are approximately the same
if return_bool
no_differences = true
min_thres = 5e-3
max_val = 0
for ind in 1:length(total_hyperparameters)
est = est_dΣdθs[ind,:,:]
est[est .== 0] .= 1e-8
val = mean(abs.(difs[ind,:,:] ./ est))
max_val = maximum([max_val, val])
if val>min_thres
no_differences = false
if print_stuff; println("dΣ/dθ$ind has a average ratioed difference of: ", val) end
end
end
if print_stuff; println("Maximum average dΣ/dθ ratioed difference of: ", max_val) end
return no_differences
end
return return_vec
end
"""
est_∇(f::Function, inputs::Vector{<:Real}; dif=1e-7, ignore_0_inputs=false)
Estimate the gradient of `f` at `inputs` with forward differences
# Examples
```jldoctest
julia> f(x) = x[1] + x[2]^2;
julia> isapprox([1,4], GPLinearODEMaker.est_∇(f, [2, 2.]); rtol=1e-5)
true
```
"""
function est_∇(f::Function, inputs::Vector{<:Real}; dif::Real=1e-7, ignore_0_inputs::Bool=false)
# original value
val = f(inputs)
#estimate gradient
j = 1
if ignore_0_inputs
grad = zeros(length(remove_zeros(inputs)))
else
grad = zeros(length(inputs))
end
for i in 1:length(inputs)
if !ignore_0_inputs || inputs[i]!=0
hold = copy(inputs)
hold[i] += dif
grad[j] = (f(hold) - val) / dif
j += 1
end
end
return grad
end
"""
est_∇nlogL_GLOM(glo, total_hyperparameters; dif=1e-7)
Estimate the gradient of `nlogL_GLOM(glo, total_hyperparameters)` with forward
differences
"""
function est_∇nlogL_GLOM(glo::GLO, total_hyperparameters::Vector{T}; dif::Real=1e-7) where {T<:Real}
f(inputs) = nlogL_GLOM(glo, inputs)
return est_∇(f, total_hyperparameters; dif=dif, ignore_0_inputs=true)
end
"""
test_∇(est_G, G; print_stuff=true, function_name="function")
Check if two gradient vectors (`est_G` and `G`) are approximately the same
# Examples
```jldoctest
julia> GPLinearODEMaker.test_∇([1., 2, 3], [1.0001, 2, 3]; print_stuff=false)
true
julia> GPLinearODEMaker.test_∇([1., 2, 3], [0., 2, 3]; print_stuff=false)
false
```
"""
function test_∇(est_G::Vector{T}, G::Vector{T}; print_stuff::Bool=true, function_name::String="function") where {T<:Real}
if print_stuff
println()
println("test_∇: Check that our analytical $function_name is close to numerical estimates")
# println("only values for non-zero hyperparameters are shown!")
# println("hyperparameters: ", total_hyperparameters)
println("analytical: ", G)
println("numerical : ", est_G)
end
no_mismatch = true
for i in 1:length(G)
if !isapprox(G[i], est_G[i], rtol=5e-2)
if print_stuff
println("mismatch df/dθ" * string(i))
end
no_mismatch = false
end
end
return no_mismatch
end
"""
test_∇nlogL_GLOM(glo, kernel_hyperparameters; dif=1e-4, print_stuff=true)
Check if `∇nlogL_GLOM(glo, total_hyperparameters)` is close to numerical
estimates provided by `est_∇nlogL_GLOM(glo, total_hyperparameters)`
"""
function test_∇nlogL_GLOM(glo::GLO, kernel_hyperparameters::Vector{T}; dif::Real=1e-4, print_stuff::Bool=true) where {T<:Real}
total_hyperparameters = append!(collect(Iterators.flatten(glo.a)), kernel_hyperparameters)
return test_∇(est_∇nlogL_GLOM(glo, total_hyperparameters; dif=dif),
∇nlogL_GLOM(glo, total_hyperparameters);
print_stuff=print_stuff, function_name="∇nlogL_GLOM")
end
"""
est_∇∇(g::Function, inputs::Vector{<:Real}; dif=1e-7, ignore_0_inputs=false)
Estimate the Hessian of a function whose gradients are provided by `g` at
`inputs` with forward differences
# Examples
```jldoctest
julia> g(x) = [(x[2] ^ 2) / 2, x[1] * x[2]]; # h(x) = [0 x[2], x[2] x[1]]
julia> isapprox([0. 9; 9 4], GPLinearODEMaker.est_∇∇(g, [4., 9]); rtol=1e-5)
true
```
"""
function est_∇∇(g::Function, inputs::Vector{<:Real}; dif::Real=1e-7, ignore_0_inputs::Bool=false)
val = g(inputs)
#estimate hessian
j = 1
dim = length(remove_zeros(inputs))
hess = zeros(dim, dim)
hold = copy(inputs)
for i in 1:length(inputs)
if !ignore_0_inputs || inputs[i]!=0
hold[i] += dif
hess[j, :] = (g(hold) - val) / dif
j += 1
hold[i] = inputs[i]
end
end
return symmetric_A(hess)
end
"""
est_∇∇_from_f(f::Function, inputs::Vector{<:Real}; dif=1e-7, ignore_0_inputs=false)
Estimate the Hessian of `f` at `inputs` with forward differences. WARNING: The
result is very sensitive to `dif`
# Examples
```jldoctest
julia> f(x) = (x[1] * x[2] ^ 2) / 2; # h(x) = [0 x[2], x[2] x[1]]
julia> isapprox([0. 9; 9 4], GPLinearODEMaker.est_∇∇_from_f(f, [4., 9]; dif=1e-4); rtol=1e-3)
true
```
"""
function est_∇∇_from_f(f::Function, inputs::Vector{<:Real}; dif::Real=1e-7, ignore_0_inputs::Bool=false)
val = est_∇(f, inputs; dif=dif, ignore_0_inputs=ignore_0_inputs)
#estimate hessian
j = 1
hess = zeros(length(inputs), length(inputs))
hold = copy(inputs)
for i in 1:length(inputs)
if !ignore_0_inputs || inputs[i]!=0
hold[i] += dif
hess[j, :] = (est_∇(f, hold; dif=dif, ignore_0_inputs=ignore_0_inputs) - val) / dif
j += 1
hold[i] = inputs[i]
end
end
return symmetric_A(hess)
end
"""
est_∇∇nlogL_GLOM(glo, total_hyperparameters; dif=1e-4)
Estimate the Hessian of `nlogL_GLOM(glo, total_hyperparameters)` with forward
differences
"""
function est_∇∇nlogL_GLOM(glo::GLO, total_hyperparameters::Vector{T}; dif::Real=1e-4) where {T<:Real}
g(inputs) = ∇nlogL_GLOM(glo, inputs)
return est_∇∇(g, total_hyperparameters; ignore_0_inputs=true)
end
"""
test_∇∇(est_H, H; print_stuff=true, function_name="function", rtol=1e-3)
Check if two Hessian matrices (`est_H` and `H`) are approximately the same
# Examples
```jldoctest
julia> GPLinearODEMaker.test_∇∇([1. 2; 3 4], [1.0001 2; 3 4]; print_stuff=false)
true
julia> GPLinearODEMaker.test_∇∇([1. 2; 3 4], [0. 2; 3 4]; print_stuff=false)
false
```
"""
function test_∇∇(est_H::Union{Symmetric,Matrix{T}}, H::Union{Symmetric,Matrix{T}}; print_stuff::Bool=true, function_name::String="function", rtol::Real=1e-3) where {T<:Real}
if print_stuff
println()
println("test_∇∇: Check that our analytical $function_name is close to numerical estimates")
# println("only values for non-zero hyperparameters are shown!")
# println("hyperparameters: ", total_hyperparameters)
println("analytical:")
for i in 1:size(H, 1)
println(H[i, :])
end
println("numerical:")
for i in 1:size(est_H, 1)
println(est_H[i, :])
end
end
no_mismatch = true
matches = fill(1, size(H))
for i in 1:size(H, 1)
for j in 1:size(H, 2)
if !(isapprox(H[i, j], est_H[i, j], rtol=rtol))
# println("mismatch d2nlogL/dθ" * string(i) * "dθ" * string(j))
matches[i, j] = 0
no_mismatch = false
end
end
end
if !no_mismatch && print_stuff
println("mismatches at 0s")
for i in 1:size(H, 1)
println(matches[i, :])
end
end
return no_mismatch
end
"""
test_∇∇nlogL_GLOM(glo, kernel_hyperparameters; dif=1e-4, print_stuff=true)
Check if `∇∇nlogL_GLOM(glo, total_hyperparameters)` is close to numerical
estimates provided by `est_∇∇nlogL_GLOM(glo, total_hyperparameters)`
"""
function test_∇∇nlogL_GLOM(glo::GLO, kernel_hyperparameters::Vector{T}; dif::Real=1e-4, print_stuff::Bool=true) where {T<:Real}
total_hyperparameters = append!(collect(Iterators.flatten(glo.a)), kernel_hyperparameters)
H = ∇∇nlogL_GLOM(glo, total_hyperparameters)
est_H = est_∇∇nlogL_GLOM(glo, total_hyperparameters; dif=dif)
return test_∇∇(est_H, H; print_stuff=print_stuff, function_name="∇∇nlogL_GLOM")
end