1- using OrdinaryDiffEq, ForwardDiff, GTPSA, Test
1+ using OrdinaryDiffEq, GTPSA, Test
2+ using DifferentiationInterface
3+ using ADTypes: AutoForwardDiff
4+ using ForwardDiff
5+
6+ # GTPSA is itself an AD engine - these tests compare GTPSA jacobian/hessian
7+ # results against ForwardDiff as a reference implementation
28
39# ODEProblem 1 =======================
410
@@ -20,23 +26,28 @@ sol_GTPSA = solve(prob_GTPSA, Tsit5(), reltol = 1.0e-16, abstol = 1.0e-16)
2026
2127@test sol. u[end ] ≈ scalar .(sol_GTPSA. u[end ]) # scalar gets 0th order part
2228
23- # Compare Jacobian against ForwardDiff
24- J_FD = ForwardDiff . jacobian ([x ... , p ... ]) do t
29+ # Compare Jacobian against AD backends using DifferentiationInterface
30+ function sol_end_problem1 (t)
2531 prob = ODEProblem (f!, t[1 : 3 ], (0.0 , 1.0 ), t[4 : 6 ])
2632 sol = solve (prob, Tsit5 (), reltol = 1.0e-16 , abstol = 1.0e-16 )
27- sol. u[end ]
33+ return sol. u[end ]
2834end
2935
30- @test J_FD ≈ GTPSA. jacobian (sol_GTPSA. u[end ], include_params = true )
36+ @testset " GTPSA Problem 1 Jacobian tests" begin
37+ J_AD = DifferentiationInterface. jacobian (sol_end_problem1, AutoForwardDiff (), [x... , p... ])
38+ @test J_AD ≈ GTPSA. jacobian (sol_GTPSA. u[end ], include_params = true )
39+ end
3140
32- # Compare Hessians against ForwardDiff
33- for i in 1 : 3
34- Hi_FD = ForwardDiff. hessian ([x... , p... ]) do t
35- prob = ODEProblem (f!, t[1 : 3 ], (0.0 , 1.0 ), t[4 : 6 ])
36- sol = solve (prob, Tsit5 (), reltol = 1.0e-16 , abstol = 1.0e-16 )
37- sol. u[end ][i]
41+ @testset " GTPSA Problem 1 Hessian tests" begin
42+ for i in 1 : 3
43+ function sol_end_i_problem1 (t)
44+ prob = ODEProblem (f!, t[1 : 3 ], (0.0 , 1.0 ), t[4 : 6 ])
45+ sol = solve (prob, Tsit5 (), reltol = 1.0e-16 , abstol = 1.0e-16 )
46+ return sol. u[end ][i]
47+ end
48+ Hi_AD = DifferentiationInterface. hessian (sol_end_i_problem1, AutoForwardDiff (), [x... , p... ])
49+ @test Hi_AD ≈ GTPSA. hessian (sol_GTPSA. u[end ][i], include_params = true )
3850 end
39- @test Hi_FD ≈ GTPSA. hessian (sol_GTPSA. u[end ][i], include_params = true )
4051end
4152
4253# ODEProblem 2 =======================
@@ -49,31 +60,36 @@ function qdot!(dp, p, q, params, t)
4960 ]
5061end
5162
52- prob = DynamicalODEProblem (pdot!, qdot!, [0.0 , 0.0 , 0.0 ], [0.0 , 0.0 , 0.0 ], (0.0 , 25.0 ))
53- sol = solve (prob , Yoshida6 (), dt = 1.0 , reltol = 1.0e-16 , abstol = 1.0e-16 )
63+ prob2 = DynamicalODEProblem (pdot!, qdot!, [0.0 , 0.0 , 0.0 ], [0.0 , 0.0 , 0.0 ], (0.0 , 25.0 ))
64+ sol2 = solve (prob2 , Yoshida6 (), dt = 1.0 , reltol = 1.0e-16 , abstol = 1.0e-16 )
5465
55- desc = Descriptor (6 , 2 ) # 6 variables to 2nd order
56- dx = @vars (desc ) # identity map
57- prob_GTPSA = DynamicalODEProblem (pdot!, qdot!, dx [1 : 3 ], dx [4 : 6 ], (0.0 , 25.0 ))
58- sol_GTPSA = solve (prob_GTPSA , Yoshida6 (), dt = 1.0 , reltol = 1.0e-16 , abstol = 1.0e-16 )
66+ desc2 = Descriptor (6 , 2 ) # 6 variables to 2nd order
67+ dx2 = @vars (desc2 ) # identity map
68+ prob_GTPSA2 = DynamicalODEProblem (pdot!, qdot!, dx2 [1 : 3 ], dx2 [4 : 6 ], (0.0 , 25.0 ))
69+ sol_GTPSA2 = solve (prob_GTPSA2 , Yoshida6 (), dt = 1.0 , reltol = 1.0e-16 , abstol = 1.0e-16 )
5970
60- @test sol . u[end ] ≈ scalar .(sol_GTPSA . u[end ]) # scalar gets 0th order part
71+ @test sol2 . u[end ] ≈ scalar .(sol_GTPSA2 . u[end ]) # scalar gets 0th order part
6172
62- # Compare Jacobian against ForwardDiff
63- J_FD = ForwardDiff . jacobian ( zeros ( 6 )) do t
73+ # Compare Jacobian against AD backends using DifferentiationInterface
74+ function sol_end_problem2 (t)
6475 prob = DynamicalODEProblem (pdot!, qdot!, t[1 : 3 ], t[4 : 6 ], (0.0 , 25.0 ))
6576 sol = solve (prob, Yoshida6 (), dt = 1.0 , reltol = 1.0e-16 , abstol = 1.0e-16 )
66- sol. u[end ]
77+ return sol. u[end ]
6778end
6879
69- @test J_FD ≈ GTPSA. jacobian (sol_GTPSA. u[end ], include_params = true )
80+ @testset " GTPSA Problem 2 Jacobian tests" begin
81+ J_AD = DifferentiationInterface. jacobian (sol_end_problem2, AutoForwardDiff (), zeros (6 ))
82+ @test J_AD ≈ GTPSA. jacobian (sol_GTPSA2. u[end ], include_params = true )
83+ end
7084
71- # Compare Hessians against ForwardDiff
72- for i in 1 : 6
73- Hi_FD = ForwardDiff. hessian (zeros (6 )) do t
74- prob = DynamicalODEProblem (pdot!, qdot!, t[1 : 3 ], t[4 : 6 ], (0.0 , 25.0 ))
75- sol = solve (prob, Yoshida6 (), dt = 1.0 , reltol = 1.0e-16 , abstol = 1.0e-16 )
76- sol. u[end ][i]
85+ @testset " GTPSA Problem 2 Hessian tests" begin
86+ for i in 1 : 6
87+ function sol_end_i_problem2 (t)
88+ prob = DynamicalODEProblem (pdot!, qdot!, t[1 : 3 ], t[4 : 6 ], (0.0 , 25.0 ))
89+ sol = solve (prob, Yoshida6 (), dt = 1.0 , reltol = 1.0e-16 , abstol = 1.0e-16 )
90+ return sol. u[end ][i]
91+ end
92+ Hi_AD = DifferentiationInterface. hessian (sol_end_i_problem2, AutoForwardDiff (), zeros (6 ))
93+ @test Hi_AD ≈ GTPSA. hessian (sol_GTPSA2. u[end ][i], include_params = true )
7794 end
78- @test Hi_FD ≈ GTPSA. hessian (sol_GTPSA. u[end ][i], include_params = true )
7995end
0 commit comments