1+ #define BOOST_TEST_DYN_LINK
2+ #include < boost/test/unit_test.hpp>
3+
4+ namespace tt = boost::test_tools;
5+
6+ #include " ../const.h"
7+ #include < adolc/adolc.h>
8+ #include < array>
9+ #include < numeric>
10+ #include < vector>
11+
12+ // Euler step (double version)
13+ int euler_step (size_t n, double *y) {
14+ y[0 ] = y[0 ] + 0.01 * y[0 ];
15+ y[1 ] = y[1 ] + 0.01 * 2 * y[1 ];
16+ return 1 ;
17+ }
18+
19+ // Euler step (adouble version)
20+ int euler_step_act (size_t n, adouble *y) {
21+ y[0 ] = y[0 ] + 0.01 * y[0 ];
22+ y[1 ] = y[1 ] + 0.01 * 2 * y[1 ];
23+ return 1 ;
24+ }
25+
26+ BOOST_AUTO_TEST_SUITE (test_checkpoint_example)
27+ BOOST_AUTO_TEST_CASE(Checkpointing_Gradient_Comparison) {
28+ const int16_t tag_full = 1 ; // Tag for full taping
29+ const int16_t tag_part = 2 ; // Tag for partial taping with checkpointing
30+ const int16_t tag_check = 3 ; // Tag for checkpointing
31+ const size_t n = 2 ; // Number of state variables
32+ const int steps = 100 ; // Number of time steps
33+ const double t0 = 0.0 ; // Initial time
34+ const double tf = 1.0 ; // Final time
35+
36+ // State variables (double and adouble versions)
37+ std::vector<double > y_double (n);
38+ ensureContiguousLocations (2 * n);
39+ std::vector<adouble> y_adouble_1 (n);
40+ std::vector<adouble> y_adouble_2 (n);
41+
42+ // Control variables (double and adouble versions)
43+ std::vector<double > conp{1.0 , 1.0 }; // Initial control values
44+ std::vector<adouble> con (n);
45+
46+ // Target value and gradient
47+ std::vector<double > out (2 );
48+ std::vector<double > grad_full (n); // Gradient from full taping
49+ std::vector<double > grad_part (n); // Gradient from checkpointing
50+
51+ // Full taping of the time step loop
52+ trace_on (tag_full);
53+ con[0 ] <<= conp[0 ];
54+ con[1 ] <<= conp[1 ];
55+ y_adouble_1[0 ] = con[0 ];
56+ y_adouble_1[1 ] = con[1 ];
57+
58+ for (int i = 0 ; i < steps; i++) {
59+ euler_step_act (n, y_adouble_1.data ());
60+ }
61+ y_adouble_1[0 ] + y_adouble_1[1 ] >>= out[0 ];
62+ trace_off ();
63+
64+ // Compute gradient using full taping
65+ gradient (tag_full, n, conp.data (), grad_full.data ());
66+
67+ // Checkpointing setup
68+ CP_Context cpc (euler_step_act); // Checkpointing context
69+ cpc.setDoubleFct (euler_step); // Double version of the time step function
70+ cpc.setNumberOfSteps (steps); // Number of time steps
71+ cpc.setNumberOfCheckpoints (5 ); // Number of checkpoints
72+ cpc.setDimensionXY (n); // Dimension of input/output
73+ cpc.setInput (y_adouble_2.data ()); // Input vector
74+ cpc.setOutput (y_adouble_2.data ()); // Output vector
75+ cpc.setTapeNumber (tag_check); // Tape number for checkpointing
76+ cpc.setAlwaysRetaping (false ); // Do not always retape
77+
78+ // Partial taping with checkpointing
79+ trace_on (tag_part);
80+ con[0 ] <<= conp[0 ];
81+ con[1 ] <<= conp[1 ];
82+ y_adouble_2[0 ] = con[0 ];
83+ y_adouble_2[1 ] = con[1 ];
84+
85+ cpc.checkpointing (); // Perform checkpointing
86+
87+ y_adouble_2[0 ] + y_adouble_2[1 ] >>= out[1 ];
88+ trace_off ();
89+
90+ // test if both taping results are equal
91+ BOOST_TEST (out[0 ] == out[1 ], tt::tolerance (tol));
92+ // Compute gradient using checkpointing
93+ gradient (tag_part, n, conp.data (), grad_part.data ());
94+ // Compare gradients from full taping and checkpointing
95+ for (size_t i = 0 ; i < n; i++) {
96+ BOOST_TEST (grad_full[i] == grad_part[i], tt::tolerance (tol));
97+ }
98+ }
99+ BOOST_AUTO_TEST_CASE (Checkpointing_fov_reverse) {
100+ const int16_t tag_full = 1 ; // Tag for full taping
101+ const int16_t tag_part = 2 ; // Tag for partial taping with checkpointing
102+ const int16_t tag_check = 3 ; // Tag for checkpointing
103+ const size_t n = 2 ; // Number of state variables
104+ const int steps = 100 ; // Number of time steps
105+ const double t0 = 0.0 ; // Initial time
106+ const double tf = 1.0 ; // Final time
107+
108+ // State variables (double and adouble versions)
109+ std::vector<double > y_double (n);
110+ ensureContiguousLocations (2 * n);
111+ std::vector<adouble> y_adouble_1 (n);
112+ std::vector<adouble> y_adouble_2 (n);
113+
114+ // Control variables (double and adouble versions)
115+ std::vector<double > conp{1.0 , 1.0 }; // Initial control values
116+ std::vector<adouble> con (n);
117+
118+ // Target value and gradient
119+ std::vector<double > out (2 );
120+ std::vector<double > grad_full (n); // Gradient from full taping
121+ std::vector<double > grad_part (n); // Gradient from checkpointing
122+
123+ // Full taping of the time step loop
124+ trace_on (tag_full, 1 );
125+ con[0 ] <<= conp[0 ];
126+ con[1 ] <<= conp[1 ];
127+ y_adouble_1[0 ] = con[0 ];
128+ y_adouble_1[1 ] = con[1 ];
129+
130+ for (int i = 0 ; i < steps; i++) {
131+ euler_step_act (n, y_adouble_1.data ());
132+ }
133+ y_adouble_1[0 ] + y_adouble_1[1 ] >>= out[0 ];
134+ trace_off ();
135+
136+ // weights
137+ double **U = myalloc2 (2 , 1 );
138+ U[0 ][0 ] = 1.0 ;
139+ U[1 ][0 ] = -1.0 ;
140+
141+ // outputs
142+ double **Z_full = myalloc2 (2 , 2 );
143+ double **Z_part = myalloc2 (2 , 2 );
144+
145+ // Compute vector-mode reverse
146+ fov_reverse (tag_full, 1 , 2 , 2 , U, Z_full);
147+
148+ // Checkpointing setup
149+ CP_Context cpc (euler_step_act); // Checkpointing context
150+ cpc.setDoubleFct (euler_step); // Double version of the time step function
151+ cpc.setNumberOfSteps (steps); // Number of time steps
152+ cpc.setNumberOfCheckpoints (5 ); // Number of checkpoints
153+ cpc.setDimensionXY (n); // Dimension of input/output
154+ cpc.setInput (y_adouble_2.data ()); // Input vector
155+ cpc.setOutput (y_adouble_2.data ()); // Output vector
156+ cpc.setTapeNumber (tag_check); // Tape number for checkpointing
157+ cpc.setAlwaysRetaping (true ); // Do always retape
158+
159+ // Partial taping with checkpointing
160+ trace_on (tag_part, 1 );
161+ con[0 ] <<= conp[0 ];
162+ con[1 ] <<= conp[1 ];
163+ y_adouble_2[0 ] = con[0 ];
164+ y_adouble_2[1 ] = con[1 ];
165+
166+ cpc.checkpointing (); // Perform checkpointing
167+
168+ y_adouble_2[0 ] + y_adouble_2[1 ] >>= out[1 ];
169+ trace_off ();
170+
171+ // test if both taping results are equal
172+ BOOST_TEST (out[0 ] == out[1 ], tt::tolerance (tol));
173+
174+ // Compute gradient using checkpointing
175+ fov_reverse (tag_part, 1 , 2 , 2 , U, Z_part);
176+ // Compare gradients from full taping and checkpointing
177+ for (size_t i = 0 ; i < 2 ; ++i) {
178+ for (size_t j = 0 ; j < 2 ; ++j)
179+ BOOST_TEST (Z_full[i][j] == Z_part[i][j], tt::tolerance (tol));
180+ }
181+ myfree2 (U);
182+ myfree2 (Z_full);
183+ myfree2 (Z_part);
184+ }
185+ BOOST_AUTO_TEST_SUITE_END ()
0 commit comments