-
-
Notifications
You must be signed in to change notification settings - Fork 27
/
ad.jl
112 lines (94 loc) · 4.18 KB
/
ad.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
using DelayDiffEq
import FiniteDiff
import ForwardDiff
using Test
# Hutchinson's equation
f(u, h, p, t) = p[2] * u * (1 - h(p, t - p[1]) / p[3])
h(p, t) = p[4]
@testset "Gradient" begin
function test(p)
prob = DDEProblem(f, p[5], h, eltype(p).((0.0, 10.0)), copy(p);
constant_lags = (p[1],))
sol = solve(prob, MethodOfSteps(Tsit5()); abstol = 1e-14, reltol = 1e-14,
saveat = 1.0)
sum(sol)
end
# without delay length estimation
p = [1.5, 1.0, 0.5]
findiff = FiniteDiff.finite_difference_gradient(p -> test(vcat(1, p, p[end])), p)
fordiff = ForwardDiff.gradient(p -> test(vcat(1, p, p[end])), p)
@test maximum(abs.(findiff .- fordiff)) < 1e-6
# with delay length estimation and without discontinuity
p = [1.0, 1.5, 1.0, 0.5]
findiff2 = FiniteDiff.finite_difference_gradient(p -> test(vcat(p, p[end])), p)
fordiff2 = ForwardDiff.gradient(p -> test(vcat(p, p[end])), p)
@test_broken maximum(abs.(findiff2 .- fordiff2)) < 2
# with discontinuity and without delay length estimation
p = [1.5, 1.0, 0.5, 1.0]
findiff3 = FiniteDiff.finite_difference_gradient(p -> test(vcat(1, p)), p)
fordiff3 = ForwardDiff.gradient(p -> test(vcat(1, p)), p)
@test maximum(abs.(findiff3 .- fordiff3)) < 9
# consistency checks
@test findiff2[2:end] ≈ findiff
@test fordiff2[2:end] ≈ fordiff
@test_broken findiff3[1:(end - 1)] ≈ findiff
@test_broken fordiff3[1:(end - 1)] ≈ fordiff
end
@testset "Jacobian" begin
function test(p)
prob = DDEProblem(f, p[5], h, eltype(p).((0.0, 10.0)), copy(p);
constant_lags = (p[1],))
sol = solve(prob, MethodOfSteps(Tsit5()); abstol = 1e-14, reltol = 1e-14,
saveat = 1.0)
sol.u
end
# without delay length estimation
p = [1.5, 1.0, 0.5]
findiff = FiniteDiff.finite_difference_jacobian(p -> test(vcat(1, p, p[end])), p)
fordiff = ForwardDiff.jacobian(p -> test(vcat(1, p, p[end])), p)
@test maximum(abs.(findiff .- fordiff)) < 2e-6
# with delay length estimation and without discontinuity
p = [1.0, 1.5, 1.0, 0.5]
findiff2 = FiniteDiff.finite_difference_jacobian(p -> test(vcat(p, p[end])), p)
fordiff2 = ForwardDiff.jacobian(p -> test(vcat(p, p[end])), p)
@test_broken maximum(abs.(findiff2 .- fordiff2)) < 3
# with discontinuity and without delay length estimation
p = [1.5, 1.0, 0.5, 1.0]
findiff3 = FiniteDiff.finite_difference_jacobian(p -> test(vcat(1, p)), p)
fordiff3 = ForwardDiff.jacobian(p -> test(vcat(1, p)), p)
@test maximum(abs.(findiff3 .- fordiff3)) < 1
# consistency checks
@test findiff2[:, 2:end] ≈ findiff
@test fordiff2[:, 2:end] ≈ fordiff
@test_broken findiff3[:, 1:(end - 1)] ≈ findiff
@test_broken fordiff3[:, 1:(end - 1)] ≈ fordiff
end
@testset "Hessian" begin
function test(p)
prob = DDEProblem(f, p[5], h, eltype(p).((0.0, 10.0)), copy(p);
constant_lags = (p[1],))
sol = solve(prob, MethodOfSteps(Tsit5()); abstol = 1e-14, reltol = 1e-14,
saveat = 1.0)
sum(sol)
end
# without delay length estimation and without discontinuity
p = [1.5, 1.0, 0.5]
findiff = FiniteDiff.finite_difference_hessian(p -> test(vcat(1, p, p[end])), p)
fordiff = ForwardDiff.hessian(p -> test(vcat(1, p, p[end])), p)
@test maximum(abs.(findiff .- fordiff)) < 1e-5
# with delay length estimation and without discontinuity
p = [1.0, 1.5, 1.0, 0.5]
findiff2 = FiniteDiff.finite_difference_hessian(p -> test(vcat(p, p[end])), p)
fordiff2 = ForwardDiff.hessian(p -> test(vcat(p, p[end])), p)
@test_broken maximum(abs.(findiff2 .- fordiff2)) < 25
# with discontinuity and without delay length estimation
p = [1.5, 1.0, 0.5, 1.0]
findiff3 = FiniteDiff.finite_difference_hessian(p -> test(vcat(1, p)), p)
fordiff3 = ForwardDiff.hessian(p -> test(vcat(1, p)), p)
@test maximum(abs.(findiff3 .- fordiff3)) < 1
# consistency checks
@test findiff2[2:end, 2:end] ≈ findiff
@test fordiff2[2:end, 2:end] ≈ fordiff
@test_broken findiff3[1:(end - 1), 1:(end - 1)] ≈ findiff
@test_broken fordiff3[1:(end - 1), 1:(end - 1)] ≈ fordiff
end