source: trunk/test_more/forward_dir.cpp @ 3682

Last change on this file since 3682 was 3682, checked in by bradbell, 5 years ago

merge to branch: trunk
from repository: https://github.com/coin-or/CppAD
start hash code: c85ed888ce93510f0df134d9fec6b7e185552f41
end hash code: febe930d34888cf37df862a1bf118229b9bf37a5

commit febe930d34888cf37df862a1bf118229b9bf37a5
Author: Brad Bell <bradbell@…>
Date: Fri May 8 14:21:47 2015 -0700

  1. Automatically generated changes for deprecated auto-tools install.
  2. Fix bug in acosh_op multiple direction forward mode.


configure.ac: Add new symbols for acosh and atanh detection (this install never detects them).
acosh_op.hpp: fix bug in multiple direction forward mode.
whats_new_15.omh: mention bug fix.
forward_dir.cpp: test asinh, acosh, and atanh.

commit 5d007c969a1e8bb0f983629f8134022a47ff2d75
Author: Brad Bell <bradbell@…>
Date: Fri May 8 12:27:18 2015 -0700

Remove trialing white space.

commit d88357106a95ee156ecf2101334dc51b22162950
Author: Brad Bell <bradbell@…>
Date: Fri May 8 12:26:58 2015 -0700

Add atanh operator.

commit 03890672ea394862373e4b1aa19a79cf04a71e73
Author: Brad Bell <bradbell@…>
Date: Fri May 8 11:57:45 2015 -0700

Change order to erf asinh acosh.

commit 9657e6dae7b5ef30ea1ea79d344f655e4db4498e
Author: Brad Bell <bradbell@…>
Date: Fri May 8 05:56:33 2015 -0700

Remove trialing white space.

commit cca1d8b8f344582f4c27c3c9c16745625e7e395e
Author: Brad Bell <bradbell@…>
Date: Fri May 8 05:56:06 2015 -0700

Change Atan -> atan.

commit ffa1894b1f97ca342aa84b65b171baad03568f58
Author: Brad Bell <bradbell@…>
Date: Fri May 8 05:37:49 2015 -0700

  1. Extend atan to include atanh.
  2. Remove index commands for words in headings.
  3. Change AtanForward? -> atan_forward.
  4. Remove trailing white space.


acos_forward.omh: add heading as in tan_forward.
asin_forward.omh: add heading as in tan_forward.
atan_forward.omh: change AtanForward? -> atan_forward.

  • Property svn:keywords set to Id
File size: 42.8 KB
Line 
1/* $Id: forward_dir.cpp 3682 2015-05-08 21:59:05Z bradbell $ */
2/* --------------------------------------------------------------------------
3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-15 Bradley M. Bell
4
5CppAD is distributed under multiple licenses. This distribution is under
6the terms of the
7                    Eclipse Public License Version 1.0.
8
9A copy of this license is included in the COPYING file of this distribution.
10Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
11-------------------------------------------------------------------------- */
12// test multiple directions operators
13// MulvvOp is tested by example/forward_dir.cpp
14
15# include <limits>
16# include <cmath>
17# include <cppad/cppad.hpp>
18
19namespace {
20        using CppAD::AD;
21        using CppAD::NearEqual;
22        // ---------------------------------------------------------------------
23        // Used the check that fun is an idenity funciton
24        typedef AD<double> (*adfun)(const AD<double>&);
25        bool check_identity(adfun identity, double argument)
26        {
27                bool ok = true;
28                double eps = 10. * std::numeric_limits<double>::epsilon();
29                size_t j;
30
31
32                // domain space vector
33                size_t n = 1;
34                CPPAD_TESTVECTOR(AD<double>) ax(n);
35                ax[0] = argument;
36
37                // declare independent variables and starting recording
38                CppAD::Independent(ax);
39
40                // range space vector
41                size_t m = 1;
42                CPPAD_TESTVECTOR(AD<double>) ay(m);
43                ay[0] = identity(ax[0]);
44
45                // create f: x -> y and stop tape recording
46                CppAD::ADFun<double> f(ax, ay);
47
48                // first order Taylor coefficients
49                size_t r = 2, ell;
50                CPPAD_TESTVECTOR(double) x1(r*n), y1;
51                for(ell = 0; ell < r; ell++)
52                {       for(j = 0; j < n; j++)
53                                x1[ r * j + ell ] = double(j + ell + 1);
54                }
55                y1  = f.Forward(1, r, x1);
56                ok &= y1.size() == r*m;
57
58                // secondorder Taylor coefficients
59                CPPAD_TESTVECTOR(double) x2(r*n), y2;
60                for(ell = 0; ell < r; ell++)
61                {       for(j = 0; j < n; j++)
62                                x2[ r * j + ell ] = double(j + ell + 2);
63                }
64                y2  = f.Forward(2, r, x2);
65                ok &= y2.size() == r*m;
66                //
67                // Y_0  (t)    = F[X_0(t)] = X_0(t)
68                //             =  0.5 + 1t + 2t^2
69                double y_1_0   = 1.0;
70                double y_2_0   = 2.0;
71                //
72                // Y_1  (t)    = F[X_1(t)] = X_1(t)
73                //             =  0.5 + 2t + 3t^2
74                double y_1_1   = 2.0;
75                double y_2_1   = 3.0;
76                //
77                ok  &= NearEqual(y1[0] , y_1_0, eps, eps);
78                ok  &= NearEqual(y1[1] , y_1_1, eps, eps);
79                ok  &= NearEqual(y2[0] , y_2_0, eps, eps);
80                ok  &= NearEqual(y2[1] , y_2_1, eps, eps);
81                //
82                return ok;
83        }
84        // ---------------------------------------------------------------------
85        // AbsOp
86        bool abs_op(void)
87        {       bool ok = true;
88                double eps = 10. * std::numeric_limits<double>::epsilon();
89                size_t j;
90
91                // domain space vector
92                size_t n = 2;
93                CPPAD_TESTVECTOR(AD<double>) ax(n);
94                ax[0] = 0.5;
95                ax[1] = -1.0;
96
97                // declare independent variables and starting recording
98                CppAD::Independent(ax);
99
100                // range space vector
101                size_t m = 1;
102                CPPAD_TESTVECTOR(AD<double>) ay(m);
103                ay[0] = abs( ax[0] ) + abs( 2.0 * ax[1] );
104
105                // create f: x -> y and stop tape recording
106                CppAD::ADFun<double> f(ax, ay);
107
108                // first order Taylor coefficients
109                size_t r = 2, ell;
110                CPPAD_TESTVECTOR(double) x1(r*n), y1;
111                for(ell = 0; ell < r; ell++)
112                {       for(j = 0; j < n; j++)
113                                x1[ r * j + ell ] = double(j + ell + 1);
114                }
115                y1  = f.Forward(1, r, x1);
116                ok &= y1.size() == r*m;
117
118                // secondorder Taylor coefficients
119                CPPAD_TESTVECTOR(double) x2(r*n), y2;
120                for(ell = 0; ell < r; ell++)
121                {       for(j = 0; j < n; j++)
122                                x2[ r * j + ell ] = double(j + ell + 2);
123                }
124                y2  = f.Forward(2, r, x2);
125                ok &= y2.size() == r*m;
126                //
127                // Y_0  (t)    = F[X_0(t)]
128                //             = abs(0.5 + 1t + 2t^2) + abs( 2*(-1.0 + 2t + 3t^2 ) )
129                double y_1_0   = -3.0;
130                double y_2_0   = -4.0;
131                //
132                // Y_1  (t)    = F[X_1(t)]
133                //             = abs(0.5 + 2t + 3t^2) + abs( 2*(-1.0 + 3t + 4t^2 ) )
134                double y_1_1   = -4.0;
135                double y_2_1   = -5.0;
136                //
137                ok  &= NearEqual(y1[0] , y_1_0, eps, eps);
138                ok  &= NearEqual(y1[1] , y_1_1, eps, eps);
139                ok  &= NearEqual(y2[0] , y_2_0, eps, eps);
140                ok  &= NearEqual(y2[1] , y_2_1, eps, eps);
141                //
142                return ok;
143        }
144        // ---------------------------------------------------------------------
145        // AddpvOp
146        bool addpv_op(void)
147        {       bool ok = true;
148                double eps = 10. * std::numeric_limits<double>::epsilon();
149                size_t j;
150
151                // domain space vector
152                size_t n = 1;
153                CPPAD_TESTVECTOR(AD<double>) ax(n);
154                ax[0] = 0.5;
155
156                // declare independent variables and starting recording
157                CppAD::Independent(ax);
158
159                // range space vector
160                size_t m = 1;
161                CPPAD_TESTVECTOR(AD<double>) ay(m);
162                ay[0] = 2.0 + ax[0];
163
164                // create f: x -> y and stop tape recording
165                CppAD::ADFun<double> f(ax, ay);
166
167                // first order Taylor coefficients
168                size_t r = 2, ell;
169                CPPAD_TESTVECTOR(double) x1(r*n), y1;
170                for(ell = 0; ell < r; ell++)
171                {       for(j = 0; j < n; j++)
172                                x1[ r * j + ell ] = double(j + ell + 1);
173                }
174                y1  = f.Forward(1, r, x1);
175                ok &= y1.size() == r*m;
176
177                // secondorder Taylor coefficients
178                CPPAD_TESTVECTOR(double) x2(r*n), y2;
179                for(ell = 0; ell < r; ell++)
180                {       for(j = 0; j < n; j++)
181                                x2[ r * j + ell ] = double(j + ell + 3);
182                }
183                y2  = f.Forward(2, r, x2);
184                ok &= y2.size() == r*m;
185                //
186                // Y_0 (t)     = F[X_0(t)]
187                //             =  2.0 + (0.5 + 1t + 3t^2)
188                double y_1_0   = 1.0;
189                double y_2_0   = 3.0;
190                //
191                // Y_1 (t)     = F[X_1(t)]
192                //             =  2.0 + (0.5 + 2t + 4t^2)
193                double y_1_1   = 2.0;
194                double y_2_1   = 4.0;
195                //
196                ok  &= NearEqual(y1[0] , y_1_0, eps, eps);
197                ok  &= NearEqual(y1[1] , y_1_1, eps, eps);
198                ok  &= NearEqual(y2[0] , y_2_0, eps, eps);
199                ok  &= NearEqual(y2[1] , y_2_1, eps, eps);
200                //
201                return ok;
202        }
203        // ---------------------------------------------------------------------
204        // AddvvOp
205        bool addvv_op(void)
206        {       bool ok = true;
207                double eps = 10. * std::numeric_limits<double>::epsilon();
208                size_t j;
209
210                // domain space vector
211                size_t n = 2;
212                CPPAD_TESTVECTOR(AD<double>) ax(n);
213                ax[0] = 0.5;
214                ax[1] = 2.0;
215
216                // declare independent variables and starting recording
217                CppAD::Independent(ax);
218
219                // range space vector
220                size_t m = 1;
221                CPPAD_TESTVECTOR(AD<double>) ay(m);
222                ay[0] = ax[0] + ax[1];
223
224                // create f: x -> y and stop tape recording
225                CppAD::ADFun<double> f(ax, ay);
226
227                // first order Taylor coefficients
228                size_t r = 2, ell;
229                CPPAD_TESTVECTOR(double) x1(r*n), y1;
230                for(ell = 0; ell < r; ell++)
231                {       for(j = 0; j < n; j++)
232                                x1[ r * j + ell ] = double(j + ell + 1);
233                }
234                y1  = f.Forward(1, r, x1);
235                ok &= y1.size() == r*m;
236
237                // secondorder Taylor coefficients
238                CPPAD_TESTVECTOR(double) x2(r*n), y2;
239                for(ell = 0; ell < r; ell++)
240                {       for(j = 0; j < n; j++)
241                                x2[ r * j + ell ] = double(j + ell + 2);
242                }
243                y2  = f.Forward(2, r, x2);
244                ok &= y2.size() == r*m;
245                //
246                // Y_0 (t)     = F[X_0(t)]
247                //             =  (0.5 + 1t + 2t^2) + (2.0 + 2t + 3t^2)
248                double y_1_0   = 1.0 + 2.0;
249                double y_2_0   = 2.0 + 3.0;
250                //
251                // Y_1 (t)     = F[X_1(t)]
252                //             =  (2.0 + 2t + 3t^2) + (2.0 + 3t + 4t^2)
253                double y_1_1   = 2.0 + 3.0;
254                double y_2_1   = 3.0 + 4.0;
255                //
256                ok  &= NearEqual(y1[0] , y_1_0, eps, eps);
257                ok  &= NearEqual(y1[1] , y_1_1, eps, eps);
258                ok  &= NearEqual(y2[0] , y_2_0, eps, eps);
259                ok  &= NearEqual(y2[1] , y_2_1, eps, eps);
260                //
261                return ok;
262        }
263        // ---------------------------------------------------------------------
264        // CosOp
265        bool cos_op(void)
266        {       bool ok = true;
267                double eps = 10. * std::numeric_limits<double>::epsilon();
268                size_t j;
269
270
271                // domain space vector
272                size_t n = 1;
273                CPPAD_TESTVECTOR(AD<double>) ax(n);
274                ax[0] = 0.5;
275
276                // declare independent variables and starting recording
277                CppAD::Independent(ax);
278
279                // range space vector
280                size_t m = 1;
281                CPPAD_TESTVECTOR(AD<double>) ay(m);
282                ay[0] = cos( ax[0] );
283
284                // create f: x -> y and stop tape recording
285                CppAD::ADFun<double> f(ax, ay);
286
287                // first order Taylor coefficients
288                size_t r = 2, ell;
289                CPPAD_TESTVECTOR(double) x1(r*n), y1;
290                for(ell = 0; ell < r; ell++)
291                {       for(j = 0; j < n; j++)
292                                x1[ r * j + ell ] = double(j + ell + 1);
293                }
294                y1  = f.Forward(1, r, x1);
295                ok &= y1.size() == r*m;
296
297                // secondorder Taylor coefficients
298                CPPAD_TESTVECTOR(double) x2(r*n), y2;
299                for(ell = 0; ell < r; ell++)
300                {       for(j = 0; j < n; j++)
301                                x2[ r * j + ell ] = double(j + ell + 2);
302                }
303                y2  = f.Forward(2, r, x2);
304                ok &= y2.size() == r*m;
305                //
306                // Y_0  (t)    = F[X_0(t)]
307                //             =  cos( 0.5 + 1t + 2t^2 )
308                // Y_0' (t)    = -sin( 0.5 + 1t + 2t^2) * (1 + 4t)
309                double y_1_0   = - sin(0.5);
310                double y_2_0   = - ( cos(0.5) + 4.0 * sin(0.5) ) / 2.0;
311                //
312                // Y_1  (t)    = F[X_1(t)]
313                //             =  cos( 0.5 + 2t + 3t^2 )
314                // Y_1' (t)    = -sin( 0.5 + 2t + 3t^2) * (2 + 6t)
315                double y_1_1   = - sin(0.5) * 2.0;
316                double y_2_1   = - ( cos(0.5) * 4.0 + 6.0 * sin(0.5) ) / 2.0;
317                //
318                ok  &= NearEqual(y1[0] , y_1_0, eps, eps);
319                ok  &= NearEqual(y1[1] , y_1_1, eps, eps);
320                ok  &= NearEqual(y2[0] , y_2_0, eps, eps);
321                ok  &= NearEqual(y2[1] , y_2_1, eps, eps);
322                //
323                return ok;
324        }
325        // ---------------------------------------------------------------------
326        // CoshOp
327        bool cosh_op(void)
328        {       bool ok = true;
329                double eps = 10. * std::numeric_limits<double>::epsilon();
330                size_t j;
331
332
333                // domain space vector
334                size_t n = 1;
335                CPPAD_TESTVECTOR(AD<double>) ax(n);
336                ax[0] = 0.5;
337
338                // declare independent variables and starting recording
339                CppAD::Independent(ax);
340
341                // range space vector
342                size_t m = 1;
343                CPPAD_TESTVECTOR(AD<double>) ay(m);
344                ay[0] = cosh( ax[0] );
345
346                // create f: x -> y and stop tape recording
347                CppAD::ADFun<double> f(ax, ay);
348
349                // first order Taylor coefficients
350                size_t r = 2, ell;
351                CPPAD_TESTVECTOR(double) x1(r*n), y1;
352                for(ell = 0; ell < r; ell++)
353                {       for(j = 0; j < n; j++)
354                                x1[ r * j + ell ] = double(j + ell + 1);
355                }
356                y1  = f.Forward(1, r, x1);
357                ok &= y1.size() == r*m;
358
359                // secondorder Taylor coefficients
360                CPPAD_TESTVECTOR(double) x2(r*n), y2;
361                for(ell = 0; ell < r; ell++)
362                {       for(j = 0; j < n; j++)
363                                x2[ r * j + ell ] = double(j + ell + 2);
364                }
365                y2  = f.Forward(2, r, x2);
366                ok &= y2.size() == r*m;
367                //
368                // Y_0  (t)    = F[X_0(t)]
369                //             = cosh( 0.5 + 1t + 2t^2 )
370                // Y_0' (t)    = sinh( 0.5 + 1t + 2t^2) * (1 + 4t)
371                double y_1_0   = sinh(0.5);
372                double y_2_0   = ( sinh(0.5) * 4.0 + cosh(0.5) ) / 2.0;
373                //
374                // Y_1  (t)    = F[X_1(t)]
375                //             = cosh( 0.5 + 2t + 3t^2 )
376                // Y_1' (t)    = sinh( 0.5 + 2t + 3t^2) * (2 + 6t)
377                double y_1_1   = sinh(0.5) * 2.0;
378                double y_2_1   = ( sinh(0.5) * 6.0 + cosh(0.5) * 4.0 ) / 2.0;
379                //
380                ok  &= NearEqual(y1[0] , y_1_0, eps, eps);
381                ok  &= NearEqual(y1[1] , y_1_1, eps, eps);
382                ok  &= NearEqual(y2[0] , y_2_0, eps, eps);
383                ok  &= NearEqual(y2[1] , y_2_1, eps, eps);
384                //
385                return ok;
386        }
387        // ---------------------------------------------------------------------
388        // CExpOp
389        bool cexp_op(void)
390        {       bool ok = true;
391                double eps = 10. * std::numeric_limits<double>::epsilon();
392                size_t j;
393
394                // domain space vector
395                size_t n = 4;
396                CPPAD_TESTVECTOR(AD<double>) ax(n);
397                for(j = 0; j < n; j++)
398                        ax[j] = double(j);
399
400                // declare independent variables and starting recording
401                CppAD::Independent(ax);
402
403                // range space vector
404                size_t m = 2;
405                CPPAD_TESTVECTOR(AD<double>) ay(m);
406                ay[0] = CondExpLt(ax[0], ax[1], ax[2], ax[3]);
407                ay[1] = CondExpGt(ax[0], ax[1], ax[2], ax[3]);
408
409                // create f: x -> y and stop tape recording
410                CppAD::ADFun<double> f(ax, ay);
411
412                // first order Taylor coefficients
413                size_t r = 2, ell;
414                CPPAD_TESTVECTOR(double) x1(r*n), y1;
415                for(ell = 0; ell < r; ell++)
416                {       for(j = 0; j < n; j++)
417                                x1[ r * j + ell ] = double(j + ell + 1);
418                }
419                y1  = f.Forward(1, r, x1);
420                ok &= y1.size() == r*m;
421
422                // secondorder Taylor coefficients
423                CPPAD_TESTVECTOR(double) x2(r*n), y2;
424                for(ell = 0; ell < r; ell++)
425                {       for(j = 0; j < n; j++)
426                                x2[ r * j + ell ] = double(j + ell + 2);
427                }
428                y2  = f.Forward(2, r, x2);
429                ok &= y2.size() == r*m;
430                //
431                // Y0_0 (t)     = X2_0(t)
432                //             =  2.0 + 3t + 4t^2
433                double y0_1_0  = 3.0;
434                double y0_2_0  = 4.0;
435                //
436                // Y1_0 (t)     = X3_0(t)
437                //             =  3.0 + 4t + 5t^2
438                double y1_1_0  = 4.0;
439                double y1_2_0  = 5.0;
440                //
441                // Y0_1 (t)     = X2_1(t)
442                //             =  2.0 + 4t + 5t^2
443                double y0_1_1  = 4.0;
444                double y0_2_1  = 5.0;
445                //
446                // Y1_1 (t)     = X3_0(t)
447                //             =  3.0 + 5t + 6t^2
448                double y1_1_1  = 5.0;
449                double y1_2_1  = 6.0;
450                //
451                ok  &= NearEqual(y1[0*r+0] , y0_1_0, eps, eps);
452                ok  &= NearEqual(y1[1*r+0] , y1_1_0, eps, eps);
453                ok  &= NearEqual(y1[0*r+1] , y0_1_1, eps, eps);
454                ok  &= NearEqual(y1[1*r+1] , y1_1_1, eps, eps);
455                //
456                ok  &= NearEqual(y2[0*r+0] , y0_2_0, eps, eps);
457                ok  &= NearEqual(y2[1*r+0] , y1_2_0, eps, eps);
458                ok  &= NearEqual(y2[0*r+1] , y0_2_1, eps, eps);
459                ok  &= NearEqual(y2[1*r+1] , y1_2_1, eps, eps);
460                //
461                return ok;
462        }
463
464        // ---------------------------------------------------------------------
465        // CSumOp
466        bool csum_op(void)
467        {       bool ok = true;
468                double eps = 10. * std::numeric_limits<double>::epsilon();
469                size_t j;
470
471                // domain space vector
472                size_t n = 3;
473                CPPAD_TESTVECTOR(AD<double>) ax(n);
474                for(j = 0; j < n; j++)
475                        ax[j] = double(j);
476
477                // declare independent variables and starting recording
478                CppAD::Independent(ax);
479
480                // range space vector
481                size_t m = 1;
482                CPPAD_TESTVECTOR(AD<double>) ay(m);
483                ay[0] = 0.0;
484                for(j = 0; j < n; j++)
485                        ay[0] += ax[j];
486
487                // create f: x -> y and stop tape recording
488                CppAD::ADFun<double> f(ax, ay);
489
490                // optmize the tape so converts summation to on CSumOp operator
491                f.optimize();
492
493                // zero order
494                CPPAD_TESTVECTOR(double) x0(n);
495                for(j = 0; j < n; j++)
496                        x0[j] = double(j);
497                f.Forward(0, x0);
498
499                // first order Taylor coefficients
500                size_t r = 2, ell;
501                CPPAD_TESTVECTOR(double) x1(r*n), y1;
502                for(ell = 0; ell < r; ell++)
503                {       for(j = 0; j < n; j++)
504                                x1[ r * j + ell ] = double(j + ell + 1);
505                }
506                y1  = f.Forward(1, r, x1);
507                ok &= y1.size() == r*m;
508
509                // secondorder Taylor coefficients
510                CPPAD_TESTVECTOR(double) x2(r*n), y2;
511                for(ell = 0; ell < r; ell++)
512                {       for(j = 0; j < n; j++)
513                                x2[ r * j + ell ] = double(j + ell + 3);
514                }
515                y2  = f.Forward(2, r, x2);
516                ok &= y2.size() == r*m;
517                //
518                double check = 0.0;
519                for(j = 0; j < n; j++)
520                        check += x1[ r * j + 0];
521                ok  &= NearEqual(y1[0] , check, eps, eps);
522                //
523                check = 0.0;
524                for(j = 0; j < n; j++)
525                        check += x1[ r * j + 1];
526                ok  &= NearEqual(y1[1] , check, eps, eps);
527                //
528                check = 0.0;
529                for(j = 0; j < n; j++)
530                        check += x2[ r * j + 0];
531                ok  &= NearEqual(y2[0] , check, eps, eps);
532                //
533                check = 0.0;
534                for(j = 0; j < n; j++)
535                        check += x2[ r * j + 1];
536                ok  &= NearEqual(y2[1] , check, eps, eps);
537                //
538                return ok;
539        }
540        // ---------------------------------------------------------------------
541        // DisOp (test assuming that AddvvOp is correct)
542        double round_off(const double& x)
543        {       // std::round(x); is C++11, so we avoid using it
544                return std::floor( x + 0.5 );
545        }
546        CPPAD_DISCRETE_FUNCTION(double, round_off)
547        bool dis_op(void)
548        {       bool ok = true;
549                double eps = 10. * std::numeric_limits<double>::epsilon();
550                size_t j;
551
552                // domain space vector
553                size_t n = 1;
554                CPPAD_TESTVECTOR(AD<double>) ax(n);
555                ax[0] = 0.5;
556
557                // declare independent variables and starting recording
558                CppAD::Independent(ax);
559
560                // range space vector
561                size_t m = 1;
562                CPPAD_TESTVECTOR(AD<double>) ay(m);
563                ay[0] = round_off( ax[0] ) + ax[0];
564
565                // create f: x -> y and stop tape recording
566                CppAD::ADFun<double> f(ax, ay);
567
568                // zero order
569                CPPAD_TESTVECTOR(double) x0(n), y0;
570                x0[0] = 2.2;
571                y0  = f.Forward(0, x0);
572                ok &= y0.size() == m;
573                ok &= NearEqual(y0[0], round_off(x0[0]) + x0[0], eps, eps);
574
575                // first order Taylor coefficients
576                size_t r = 2, ell;
577                CPPAD_TESTVECTOR(double) x1(r*n), y1;
578                for(ell = 0; ell < r; ell++)
579                {       for(j = 0; j < n; j++)
580                                x1[ r * j + ell ] = double(j + ell + 1);
581                }
582                y1  = f.Forward(1, r, x1);
583                ok &= y1.size() == r*m;
584
585                // secondorder Taylor coefficients
586                CPPAD_TESTVECTOR(double) x2(r*n), y2;
587                for(ell = 0; ell < r; ell++)
588                {       for(j = 0; j < n; j++)
589                                x2[ r * j + ell ] = double(j + ell + 2);
590                }
591                y2  = f.Forward(2, r, x2);
592                ok &= y2.size() == r*m;
593                //
594                //
595                // Y_0 (t)     = F[X_0(t)]
596                //             =  2.0 + (2.2 + 1t + 2t^2)
597                double y_1_0   = 1.0;
598                double y_2_0   = 2.0;
599                //
600                // Y_1 (t)     = F[X_1(t)]
601                //             =  2.0 + (2.2 + 2t + 3t^2)
602                double y_1_1   = 2.0;
603                double y_2_1   = 3.0;
604                //
605                ok  &= NearEqual(y1[0] , y_1_0, eps, eps);
606                ok  &= NearEqual(y1[1] , y_1_1, eps, eps);
607                ok  &= NearEqual(y2[0] , y_2_0, eps, eps);
608                ok  &= NearEqual(y2[1] , y_2_1, eps, eps);
609                //
610                return ok;
611        }
612        // ---------------------------------------------------------------------
613        // DivpvOp (testing assumping MulpvOp is correct)
614        bool divpv_op(void)
615        {       bool ok = true;
616                double eps = 10. * std::numeric_limits<double>::epsilon();
617                size_t j;
618
619                // domain space vector
620                size_t n = 1;
621                CPPAD_TESTVECTOR(AD<double>) ax(n);
622                ax[0] = 0.5;
623
624                // declare independent variables and starting recording
625                CppAD::Independent(ax);
626
627                // range space vector
628                size_t m = 1;
629                CPPAD_TESTVECTOR(AD<double>) ay(m);
630                ay[0] = (2.0 / ax[0]) * (ax[0] * ax[0]);
631
632                // create f: x -> y and stop tape recording
633                CppAD::ADFun<double> f(ax, ay);
634
635                // first order Taylor coefficients
636                size_t r = 2, ell;
637                CPPAD_TESTVECTOR(double) x1(r*n), y1;
638                for(ell = 0; ell < r; ell++)
639                {       for(j = 0; j < n; j++)
640                                x1[ r * j + ell ] = double(j + ell + 1);
641                }
642                y1  = f.Forward(1, r, x1);
643                ok &= y1.size() == r*m;
644
645                // secondorder Taylor coefficients
646                CPPAD_TESTVECTOR(double) x2(r*n), y2;
647                for(ell = 0; ell < r; ell++)
648                {       for(j = 0; j < n; j++)
649                                x2[ r * j + ell ] = double(j + ell + 2);
650                }
651                y2  = f.Forward(2, r, x2);
652                ok &= y2.size() == r*m;
653                //
654                // Y_0 (t)     = F[X_0(t)]
655                //             = 2.0 * (0.5 + 1t + 2t^2)
656                double y_1_0   = 2.0;
657                double y_2_0   = 4.0;
658                //
659                // Y_1 (t)     = F[X_1(t)]
660                //             = 2.0 * (0.5 + 2t + 3t^2)/2.0
661                double y_1_1   = 4.0;
662                double y_2_1   = 6.0;
663                //
664                ok  &= NearEqual(y1[0] , y_1_0, eps, eps);
665                ok  &= NearEqual(y1[1] , y_1_1, eps, eps);
666                ok  &= NearEqual(y2[0] , y_2_0, eps, eps);
667                ok  &= NearEqual(y2[1] , y_2_1, eps, eps);
668                //
669                return ok;
670        }
671        // ---------------------------------------------------------------------
672        // DivvpOp
673        bool divvp_op(void)
674        {       bool ok = true;
675                double eps = 10. * std::numeric_limits<double>::epsilon();
676                size_t j;
677
678                // domain space vector
679                size_t n = 1;
680                CPPAD_TESTVECTOR(AD<double>) ax(n);
681                ax[0] = 0.5;
682
683                // declare independent variables and starting recording
684                CppAD::Independent(ax);
685
686                // range space vector
687                size_t m = 1;
688                CPPAD_TESTVECTOR(AD<double>) ay(m);
689                ay[0] = ax[0] / 2.0;
690
691                // create f: x -> y and stop tape recording
692                CppAD::ADFun<double> f(ax, ay);
693
694                // first order Taylor coefficients
695                size_t r = 2, ell;
696                CPPAD_TESTVECTOR(double) x1(r*n), y1;
697                for(ell = 0; ell < r; ell++)
698                {       for(j = 0; j < n; j++)
699                                x1[ r * j + ell ] = double(j + ell + 1);
700                }
701                y1  = f.Forward(1, r, x1);
702                ok &= y1.size() == r*m;
703
704                // secondorder Taylor coefficients
705                CPPAD_TESTVECTOR(double) x2(r*n), y2;
706                for(ell = 0; ell < r; ell++)
707                {       for(j = 0; j < n; j++)
708                                x2[ r * j + ell ] = double(j + ell + 3);
709                }
710                y2  = f.Forward(2, r, x2);
711                ok &= y2.size() == r*m;
712                //
713                // Y_0 (t)     = F[X_0(t)]
714                //             =  (0.5 + 1t + 3t^2)/2.0
715                double y_1_0   = 1.0 / 2.0;
716                double y_2_0   = 3.0 / 2.0;
717                //
718                // Y_1 (t)     = F[X_1(t)]
719                //             =  (0.5 + 2t + 4t^2)/2.0
720                double y_1_1   = 2.0 / 2.0;
721                double y_2_1   = 4.0 / 2.0;
722                //
723                ok  &= NearEqual(y1[0] , y_1_0, eps, eps);
724                ok  &= NearEqual(y1[1] , y_1_1, eps, eps);
725                ok  &= NearEqual(y2[0] , y_2_0, eps, eps);
726                ok  &= NearEqual(y2[1] , y_2_1, eps, eps);
727                //
728                return ok;
729        }
730        // ---------------------------------------------------------------------
731        // ExpOp
732        bool exp_op(void)
733        {       bool ok = true;
734                double eps = 10. * std::numeric_limits<double>::epsilon();
735                size_t j;
736
737                // domain space vector
738                size_t n = 1;
739                CPPAD_TESTVECTOR(AD<double>) ax(n);
740                ax[0] = 0.5;
741
742                // declare independent variables and starting recording
743                CppAD::Independent(ax);
744
745                // range space vector
746                size_t m = 1;
747                CPPAD_TESTVECTOR(AD<double>) ay(m);
748                ay[0] = exp( ax[0] );
749
750                // create f: x -> y and stop tape recording
751                CppAD::ADFun<double> f(ax, ay);
752
753                // first order Taylor coefficients
754                size_t r = 2, ell;
755                CPPAD_TESTVECTOR(double) x1(r*n), y1;
756                for(ell = 0; ell < r; ell++)
757                {       for(j = 0; j < n; j++)
758                                x1[ r * j + ell ] = double(j + ell + 1);
759                }
760                y1  = f.Forward(1, r, x1);
761                ok &= y1.size() == r*m;
762
763                // secondorder Taylor coefficients
764                CPPAD_TESTVECTOR(double) x2(r*n), y2;
765                for(ell = 0; ell < r; ell++)
766                {       for(j = 0; j < n; j++)
767                                x2[ r * j + ell ] = double(j + ell + 2);
768                }
769                y2  = f.Forward(2, r, x2);
770                ok &= y2.size() == r*m;
771                //
772                // Y_0  (t)    = F[X_0(t)]
773                //             =  exp(0.5 + 1t + 2t^2)
774                // Y_0' (t)    =  exp(0.5 + 1t + 2t^2)*(1 + 4t)
775                double y_1_0   = exp(0.5);
776                double y_2_0   = ( exp(0.5)*4.0 + exp(0.5) ) / 2.0;
777                //
778                // Y_1  (t)    = F[X_1(t)]
779                //             =  exp(0.5 + 2t + 3t^2)
780                // Y_1' (t)    =  exp(0.5 + 2t + 3t^2)*(2 + 6t)
781                double y_1_1   = exp(0.5)*2.0;
782                double y_2_1   = ( exp(0.5)*6.0 + exp(0.5)*2.0*2.0 ) / 2.0;
783                //
784                ok  &= NearEqual(y1[0] , y_1_0, eps, eps);
785                ok  &= NearEqual(y1[1] , y_1_1, eps, eps);
786                ok  &= NearEqual(y2[0] , y_2_0, eps, eps);
787                ok  &= NearEqual(y2[1] , y_2_1, eps, eps);
788                //
789                return ok;
790        }
791        // ---------------------------------------------------------------------
792        // LdpOp and LdvOp (test assuming AdvvOp is correct)
793        bool load_op(void)
794        {       bool ok = true;
795                double eps = 10. * std::numeric_limits<double>::epsilon();
796                size_t j;
797
798                // domain space vector
799                size_t n = 2;
800                CPPAD_TESTVECTOR(AD<double>) ax(n);
801                ax[0] = 0.0;
802                ax[1] = 1.0;
803
804                // declare independent variables and starting recording
805                CppAD::Independent(ax);
806
807                // Store operations
808                CppAD::VecAD<double> avec(3);
809                avec[ AD<double>(0) ]    = ax[0];  // store a variable
810                avec[ AD<double>(1) ]    = ax[1];  // store a variable
811                avec[ AD<double>(2) ]    = 5.0;    // store a parameter
812
813                // range space vector
814                size_t m = 2;
815                CPPAD_TESTVECTOR(AD<double>) ay(m);
816                ay[0] = avec[ AD<double>(0) ];     // load using parameter index
817                ay[1] = avec[ ax[0] ];             // load using variable index
818
819                // create f: x -> y and stop tape recording
820                CppAD::ADFun<double> f(ax, ay);
821
822                // zero order Taylor coefficients
823                CPPAD_TESTVECTOR(double) x0(n), y0;
824                x0[0] = 2;
825                x0[1] = 3;
826                y0  = f.Forward(0, x0);
827                ok &= y0.size() == m;
828                // y[0] = avec[0] = x[0]
829                ok &= y0[0] == x0[0];
830                // y[1] = avec[ x[0] ] = avec[2] = 5.0
831                ok &= y0[1] == 5.0;
832
833                // first order Taylor coefficients
834                size_t r = 2, ell;
835                CPPAD_TESTVECTOR(double) x1(r*n), y1;
836                for(ell = 0; ell < r; ell++)
837                {       for(j = 0; j < n; j++)
838                                x1[ r * j + ell ] = double(j + ell + 1);
839                }
840                y1  = f.Forward(1, r, x1);
841                ok &= y1.size() == r*m;
842
843                // secondorder Taylor coefficients
844                CPPAD_TESTVECTOR(double) x2(r*n), y2;
845                for(ell = 0; ell < r; ell++)
846                {       for(j = 0; j < n; j++)
847                                x2[ r * j + ell ] = double(j + ell + 2);
848                }
849                y2  = f.Forward(2, r, x2);
850                ok &= y2.size() == r*m;
851                //
852                // Y0_0 (t)    = 2.0 + 1t + 2t^2
853                double y0_1_0  = 1.0;
854                double y0_2_0  = 2.0;
855                //
856                // Y1_0 (t)    = 5.0
857                double y1_1_0  = 0.0;
858                double y1_2_0  = 0.0;
859                //
860                // Y0_1 (t)    = 2.0 + 2t + 3t^2
861                double y0_1_1  = 2.0;
862                double y0_2_1  = 3.0;
863                //
864                // Y1_1 (t)    = 5.0
865                double y1_1_1  = 0.0;
866                double y1_2_1  = 0.0;
867                //
868                ok  &= NearEqual(y1[0*r+0] , y0_1_0, eps, eps);
869                ok  &= NearEqual(y1[1*r+0] , y1_1_0, eps, eps);
870                ok  &= NearEqual(y1[0*r+1] , y0_1_1, eps, eps);
871                ok  &= NearEqual(y1[1*r+1] , y1_1_1, eps, eps);
872                //
873                ok  &= NearEqual(y2[0*r+0] , y0_2_0, eps, eps);
874                ok  &= NearEqual(y2[1*r+0] , y1_2_0, eps, eps);
875                ok  &= NearEqual(y2[0*r+1] , y0_2_1, eps, eps);
876                ok  &= NearEqual(y2[1*r+1] , y1_2_1, eps, eps);
877                //
878                return ok;
879        }
880
881        // ---------------------------------------------------------------------
882        // MulpvOp
883        bool mulpv_op(void)
884        {       bool ok = true;
885                double eps = 10. * std::numeric_limits<double>::epsilon();
886                size_t j;
887
888                // domain space vector
889                size_t n = 1;
890                CPPAD_TESTVECTOR(AD<double>) ax(n);
891                ax[0] = 0.5;
892
893                // declare independent variables and starting recording
894                CppAD::Independent(ax);
895
896                // range space vector
897                size_t m = 1;
898                CPPAD_TESTVECTOR(AD<double>) ay(m);
899                ay[0] = 2.0 * ax[0];
900
901                // create f: x -> y and stop tape recording
902                CppAD::ADFun<double> f(ax, ay);
903
904                // first order Taylor coefficients
905                size_t r = 2, ell;
906                CPPAD_TESTVECTOR(double) x1(r*n), y1;
907                for(ell = 0; ell < r; ell++)
908                {       for(j = 0; j < n; j++)
909                                x1[ r * j + ell ] = double(j + ell + 1);
910                }
911                y1  = f.Forward(1, r, x1);
912                ok &= y1.size() == r*m;
913
914                // secondorder Taylor coefficients
915                CPPAD_TESTVECTOR(double) x2(r*n), y2;
916                for(ell = 0; ell < r; ell++)
917                {       for(j = 0; j < n; j++)
918                                x2[ r * j + ell ] = double(j + ell + 3);
919                }
920                y2  = f.Forward(2, r, x2);
921                ok &= y2.size() == r*m;
922                //
923                // Y_0 (t)     = F[X_0(t)]
924                //             =  2.0 * (0.5 + 1t + 3t^2)
925                double y_1_0   = 2.0 * 1.0;
926                double y_2_0   = 2.0 * 3.0;
927                //
928                // Y_1 (t)     = F[X_1(t)]
929                //             =  2.0 * (0.5 + 2t + 4t^2)
930                double y_1_1   = 2.0 * 2.0;
931                double y_2_1   = 2.0 * 4.0;
932                //
933                ok  &= NearEqual(y1[0] , y_1_0, eps, eps);
934                ok  &= NearEqual(y1[1] , y_1_1, eps, eps);
935                ok  &= NearEqual(y2[0] , y_2_0, eps, eps);
936                ok  &= NearEqual(y2[1] , y_2_1, eps, eps);
937                //
938                return ok;
939        }
940        // ---------------------------------------------------------------------
941        // ParOp
942        bool par_op(void)
943        {       bool ok = true;
944                size_t j;
945
946                // domain space vector
947                size_t n = 1;
948                CPPAD_TESTVECTOR(AD<double>) ax(n);
949                ax[0] = 0.5;
950
951                // declare independent variables and starting recording
952                CppAD::Independent(ax);
953
954                // range space vector
955                size_t m = 1;
956                CPPAD_TESTVECTOR(AD<double>) ay(m);
957                ay[0] = 0.0 * ax[0];
958
959                // create f: x -> y and stop tape recording
960                CppAD::ADFun<double> f(ax, ay);
961
962                // first order Taylor coefficients
963                size_t r = 2, ell;
964                CPPAD_TESTVECTOR(double) x1(r*n), y1;
965                for(ell = 0; ell < r; ell++)
966                {       for(j = 0; j < n; j++)
967                                x1[ r * j + ell ] = double(j + ell + 1);
968                }
969                y1  = f.Forward(1, r, x1);
970                ok &= y1.size() == r*m;
971
972                // secondorder Taylor coefficients
973                CPPAD_TESTVECTOR(double) x2(r*n), y2;
974                for(ell = 0; ell < r; ell++)
975                {       for(j = 0; j < n; j++)
976                                x2[ r * j + ell ] = double(j + ell + 2);
977                }
978                y2  = f.Forward(2, r, x2);
979                ok &= y2.size() == r*m;
980                //
981                // Y_0  (t)    = 0.0
982                for(ell = 0; ell < r; ell++)
983                {       ok &= y1[ell] == 0.0;
984                        ok &= y2[ell] == 0.0;
985                }
986                return ok;
987        }
988        // ---------------------------------------------------------------------
989        // PowvvOp  (test assuming LogOp, ExpOp and DivvvOp are correct)
990        bool powvv_op(void)
991        {       bool ok = true;
992                double eps = 10. * std::numeric_limits<double>::epsilon();
993                size_t j;
994
995                // domain space vector
996                size_t n = 2;
997                CPPAD_TESTVECTOR(AD<double>) ax(n);
998                ax[0] = 0.5;
999                ax[1] = 2.0;
1000
1001                // declare independent variables and starting recording
1002                CppAD::Independent(ax);
1003
1004                // range space vector
1005                size_t m = 2;
1006                CPPAD_TESTVECTOR(AD<double>) ay(m);
1007                ay[0] = log( pow( exp(ax[0]) , ax[1] ) ) / ax[1] ;
1008                ay[1] = log( pow( exp(ax[0]) , ax[1] ) ) / ax[0] ;
1009
1010                // create f: x -> y and stop tape recording
1011                CppAD::ADFun<double> f(ax, ay);
1012
1013                // first order Taylor coefficients
1014                size_t r = 2, ell;
1015                CPPAD_TESTVECTOR(double) x1(r*n), y1;
1016                for(ell = 0; ell < r; ell++)
1017                {       for(j = 0; j < n; j++)
1018                                x1[ r * j + ell ] = double(j + ell + 1);
1019                }
1020                y1  = f.Forward(1, r, x1);
1021                ok &= y1.size() == r*m;
1022
1023                // secondorder Taylor coefficients
1024                CPPAD_TESTVECTOR(double) x2(r*n), y2;
1025                for(ell = 0; ell < r; ell++)
1026                {       for(j = 0; j < n; j++)
1027                                x2[ r * j + ell ] = double(j + ell + 2);
1028                }
1029                y2  = f.Forward(2, r, x2);
1030                ok &= y2.size() == r*m;
1031                //
1032                // Y0_0 (t)    = 0.5 + 1t + 2t^2
1033                double y0_1_0  = 1.0;
1034                double y0_2_0  = 2.0;
1035                //
1036                // Y0_1 (t)    = 0.5 + 2t + 3t^2
1037                double y0_1_1  = 2.0;
1038                double y0_2_1  = 3.0;
1039                //
1040                // Y1_0 (t)    = 2.0 + 2t + 3t^2
1041                double y1_1_0  = 2.0;
1042                double y1_2_0  = 3.0;
1043                //
1044                // Y1_1 (t)    = 2.0 + 3t + 4t^2
1045                double y1_1_1  = 3.0;
1046                double y1_2_1  = 4.0;
1047                //
1048                ok  &= NearEqual(y1[0*r+0] , y0_1_0, eps, eps);
1049                ok  &= NearEqual(y1[1*r+0] , y1_1_0, eps, eps);
1050                ok  &= NearEqual(y1[0*r+1] , y0_1_1, eps, eps);
1051                ok  &= NearEqual(y1[1*r+1] , y1_1_1, eps, eps);
1052                //
1053                ok  &= NearEqual(y2[0*r+0] , y0_2_0, eps, eps);
1054                ok  &= NearEqual(y2[1*r+0] , y1_2_0, eps, eps);
1055                ok  &= NearEqual(y2[0*r+1] , y0_2_1, eps, eps);
1056                ok  &= NearEqual(y2[1*r+1] , y1_2_1, eps, eps);
1057                //
1058                return ok;
1059        }
1060        // ---------------------------------------------------------------------
1061        // SignOp (test assuming that MulvvOp is correct)
1062        bool sign_op(void)
1063        {       bool ok = true;
1064                double eps = 10. * std::numeric_limits<double>::epsilon();
1065                size_t j;
1066
1067                // domain space vector
1068                size_t n = 1;
1069                CPPAD_TESTVECTOR(AD<double>) ax(n);
1070                ax[0] = 0.5;
1071
1072                // declare independent variables and starting recording
1073                CppAD::Independent(ax);
1074
1075                // range space vector
1076                size_t m = 1;
1077                CPPAD_TESTVECTOR(AD<double>) ay(m);
1078                ay[0] = sign( ax[0] ) * ax[0];
1079
1080                // create f: x -> y and stop tape recording
1081                CppAD::ADFun<double> f(ax, ay);
1082
1083                // zero order
1084                CPPAD_TESTVECTOR(double) x0(n), y0;
1085                x0[0] = -3.0;
1086                y0  = f.Forward(0, x0);
1087                ok &= y0.size() == m;
1088                ok &= NearEqual(y0[0], CppAD::abs(x0[0]), eps, eps);
1089
1090                // first order Taylor coefficients
1091                size_t r = 2, ell;
1092                CPPAD_TESTVECTOR(double) x1(r*n), y1;
1093                for(ell = 0; ell < r; ell++)
1094                {       for(j = 0; j < n; j++)
1095                                x1[ r * j + ell ] = double(j + ell + 1);
1096                }
1097                y1  = f.Forward(1, r, x1);
1098                ok &= y1.size() == r*m;
1099
1100                // secondorder Taylor coefficients
1101                CPPAD_TESTVECTOR(double) x2(r*n), y2;
1102                for(ell = 0; ell < r; ell++)
1103                {       for(j = 0; j < n; j++)
1104                                x2[ r * j + ell ] = double(j + ell + 2);
1105                }
1106                y2  = f.Forward(2, r, x2);
1107                ok &= y2.size() == r*m;
1108                //
1109                //
1110                // Y_0 (t)     = F[X_0(t)]
1111                //             =  -(-3.0 + 1t + 2t^2)
1112                double y_1_0   = -1.0;
1113                double y_2_0   = -2.0;
1114                //
1115                // Y_1 (t)     = F[X_1(t)]
1116                //             =  -(-3.0 + 2t + 3t^2)
1117                double y_1_1   = -2.0;
1118                double y_2_1   = -3.0;
1119                //
1120                ok  &= NearEqual(y1[0] , y_1_0, eps, eps);
1121                ok  &= NearEqual(y1[1] , y_1_1, eps, eps);
1122                ok  &= NearEqual(y2[0] , y_2_0, eps, eps);
1123                ok  &= NearEqual(y2[1] , y_2_1, eps, eps);
1124                //
1125                return ok;
1126        }
1127
1128        // ---------------------------------------------------------------------
1129        // SinOp
1130        bool sin_op(void)
1131        {       bool ok = true;
1132                double eps = 10. * std::numeric_limits<double>::epsilon();
1133                size_t j;
1134
1135
1136                // domain space vector
1137                size_t n = 1;
1138                CPPAD_TESTVECTOR(AD<double>) ax(n);
1139                ax[0] = 0.5;
1140
1141                // declare independent variables and starting recording
1142                CppAD::Independent(ax);
1143
1144                // range space vector
1145                size_t m = 1;
1146                CPPAD_TESTVECTOR(AD<double>) ay(m);
1147                ay[0] = sin( ax[0] );
1148
1149                // create f: x -> y and stop tape recording
1150                CppAD::ADFun<double> f(ax, ay);
1151
1152                // first order Taylor coefficients
1153                size_t r = 2, ell;
1154                CPPAD_TESTVECTOR(double) x1(r*n), y1;
1155                for(ell = 0; ell < r; ell++)
1156                {       for(j = 0; j < n; j++)
1157                                x1[ r * j + ell ] = double(j + ell + 1);
1158                }
1159                y1  = f.Forward(1, r, x1);
1160                ok &= y1.size() == r*m;
1161
1162                // secondorder Taylor coefficients
1163                CPPAD_TESTVECTOR(double) x2(r*n), y2;
1164                for(ell = 0; ell < r; ell++)
1165                {       for(j = 0; j < n; j++)
1166                                x2[ r * j + ell ] = double(j + ell + 2);
1167                }
1168                y2  = f.Forward(2, r, x2);
1169                ok &= y2.size() == r*m;
1170                //
1171                // Y_0  (t)    = F[X_0(t)]
1172                //             = sin( 0.5 + 1t + 2t^2 )
1173                // Y_0' (t)    = cos( 0.5 + 1t + 2t^2) * (1 + 4t)
1174                double y_1_0   = cos(0.5);
1175                double y_2_0   = ( cos(0.5) * 4.0 - sin(0.5) ) / 2.0;
1176                //
1177                // Y_1  (t)    = F[X_1(t)]
1178                //             = sin( 0.5 + 2t + 3t^2 )
1179                // Y_1' (t)    = cos( 0.5 + 2t + 3t^2) * (2 + 6t)
1180                double y_1_1   = cos(0.5) * 2.0;
1181                double y_2_1   = ( cos(0.5) * 6.0 - sin(0.5) * 4.0 ) / 2.0;
1182                //
1183                ok  &= NearEqual(y1[0] , y_1_0, eps, eps);
1184                ok  &= NearEqual(y1[1] , y_1_1, eps, eps);
1185                ok  &= NearEqual(y2[0] , y_2_0, eps, eps);
1186                ok  &= NearEqual(y2[1] , y_2_1, eps, eps);
1187                //
1188                return ok;
1189        }
1190        // ---------------------------------------------------------------------
1191        // SinhOp
1192        bool sinh_op(void)
1193        {       bool ok = true;
1194                double eps = 10. * std::numeric_limits<double>::epsilon();
1195                size_t j;
1196
1197
1198                // domain space vector
1199                size_t n = 1;
1200                CPPAD_TESTVECTOR(AD<double>) ax(n);
1201                ax[0] = 0.5;
1202
1203                // declare independent variables and starting recording
1204                CppAD::Independent(ax);
1205
1206                // range space vector
1207                size_t m = 1;
1208                CPPAD_TESTVECTOR(AD<double>) ay(m);
1209                ay[0] = sinh( ax[0] );
1210
1211                // create f: x -> y and stop tape recording
1212                CppAD::ADFun<double> f(ax, ay);
1213
1214                // first order Taylor coefficients
1215                size_t r = 2, ell;
1216                CPPAD_TESTVECTOR(double) x1(r*n), y1;
1217                for(ell = 0; ell < r; ell++)
1218                {       for(j = 0; j < n; j++)
1219                                x1[ r * j + ell ] = double(j + ell + 1);
1220                }
1221                y1  = f.Forward(1, r, x1);
1222                ok &= y1.size() == r*m;
1223
1224                // secondorder Taylor coefficients
1225                CPPAD_TESTVECTOR(double) x2(r*n), y2;
1226                for(ell = 0; ell < r; ell++)
1227                {       for(j = 0; j < n; j++)
1228                                x2[ r * j + ell ] = double(j + ell + 2);
1229                }
1230                y2  = f.Forward(2, r, x2);
1231                ok &= y2.size() == r*m;
1232                //
1233                // Y_0  (t)    = F[X_0(t)]
1234                //             = sinh( 0.5 + 1t + 2t^2 )
1235                // Y_0' (t)    = cosh( 0.5 + 1t + 2t^2) * (1 + 4t)
1236                double y_1_0   = cosh(0.5);
1237                double y_2_0   = ( cosh(0.5) * 4.0 + sinh(0.5) ) / 2.0;
1238                //
1239                // Y_1  (t)    = F[X_1(t)]
1240                //             = sinh( 0.5 + 2t + 3t^2 )
1241                // Y_1' (t)    = cosh( 0.5 + 2t + 3t^2) * (2 + 6t)
1242                double y_1_1   = cosh(0.5) * 2.0;
1243                double y_2_1   = ( cosh(0.5) * 6.0 + sinh(0.5) * 4.0 ) / 2.0;
1244                //
1245                ok  &= NearEqual(y1[0] , y_1_0, eps, eps);
1246                ok  &= NearEqual(y1[1] , y_1_1, eps, eps);
1247                ok  &= NearEqual(y2[0] , y_2_0, eps, eps);
1248                ok  &= NearEqual(y2[1] , y_2_1, eps, eps);
1249                //
1250                return ok;
1251        }
1252        // ---------------------------------------------------------------------
1253        // SubpvOp
1254        bool subpv_op(void)
1255        {       bool ok = true;
1256                double eps = 10. * std::numeric_limits<double>::epsilon();
1257                size_t j;
1258
1259                // domain space vector
1260                size_t n = 1;
1261                CPPAD_TESTVECTOR(AD<double>) ax(n);
1262                ax[0] = 0.5;
1263
1264                // declare independent variables and starting recording
1265                CppAD::Independent(ax);
1266
1267                // range space vector
1268                size_t m = 1;
1269                CPPAD_TESTVECTOR(AD<double>) ay(m);
1270                ay[0] = 2.0 - ax[0];
1271
1272                // create f: x -> y and stop tape recording
1273                CppAD::ADFun<double> f(ax, ay);
1274
1275                // first order Taylor coefficients
1276                size_t r = 2, ell;
1277                CPPAD_TESTVECTOR(double) x1(r*n), y1;
1278                for(ell = 0; ell < r; ell++)
1279                {       for(j = 0; j < n; j++)
1280                                x1[ r * j + ell ] = double(j + ell + 1);
1281                }
1282                y1  = f.Forward(1, r, x1);
1283                ok &= y1.size() == r*m;
1284
1285                // secondorder Taylor coefficients
1286                CPPAD_TESTVECTOR(double) x2(r*n), y2;
1287                for(ell = 0; ell < r; ell++)
1288                {       for(j = 0; j < n; j++)
1289                                x2[ r * j + ell ] = double(j + ell + 3);
1290                }
1291                y2  = f.Forward(2, r, x2);
1292                ok &= y2.size() == r*m;
1293                //
1294                // Y_0 (t)     = F[X_0(t)]
1295                //             =  2.0 - (0.5 + 1t + 3t^2)/2.0
1296                double y_1_0   = - 1.0;
1297                double y_2_0   = - 3.0;
1298                //
1299                // Y_1 (t)     = F[X_1(t)]
1300                //             =  3.0 - (0.5 + 2t + 4t^2)/2.0
1301                double y_1_1   = - 2.0;
1302                double y_2_1   = - 4.0;
1303                //
1304                ok  &= NearEqual(y1[0] , y_1_0, eps, eps);
1305                ok  &= NearEqual(y1[1] , y_1_1, eps, eps);
1306                ok  &= NearEqual(y2[0] , y_2_0, eps, eps);
1307                ok  &= NearEqual(y2[1] , y_2_1, eps, eps);
1308                //
1309                return ok;
1310        }
1311        // ---------------------------------------------------------------------
1312        // SubvvOp
1313        bool subvv_op(void)
1314        {       bool ok = true;
1315                double eps = 10. * std::numeric_limits<double>::epsilon();
1316                size_t j;
1317
1318                // domain space vector
1319                size_t n = 2;
1320                CPPAD_TESTVECTOR(AD<double>) ax(n);
1321                ax[0] = 0.5;
1322                ax[1] = 2.0;
1323
1324                // declare independent variables and starting recording
1325                CppAD::Independent(ax);
1326
1327                // range space vector
1328                size_t m = 1;
1329                CPPAD_TESTVECTOR(AD<double>) ay(m);
1330                ay[0] = ax[0] - 2.0 * ax[1];
1331
1332                // create f: x -> y and stop tape recording
1333                CppAD::ADFun<double> f(ax, ay);
1334
1335                // first order Taylor coefficients
1336                size_t r = 2, ell;
1337                CPPAD_TESTVECTOR(double) x1(r*n), y1;
1338                for(ell = 0; ell < r; ell++)
1339                {       for(j = 0; j < n; j++)
1340                                x1[ r * j + ell ] = double(j + ell + 1);
1341                }
1342                y1  = f.Forward(1, r, x1);
1343                ok &= y1.size() == r*m;
1344
1345                // secondorder Taylor coefficients
1346                CPPAD_TESTVECTOR(double) x2(r*n), y2;
1347                for(ell = 0; ell < r; ell++)
1348                {       for(j = 0; j < n; j++)
1349                                x2[ r * j + ell ] = double(j + ell + 2);
1350                }
1351                y2  = f.Forward(2, r, x2);
1352                ok &= y2.size() == r*m;
1353                //
1354                // Y_0 (t)     = F[X_0(t)]
1355                //             =  (0.5 + 1t + 2t^2) - 2.0 * (2.0 + 2t + 3t^2)
1356                double y_1_0   = 1.0 - 4.0;
1357                double y_2_0   = 2.0 - 6.0;
1358                //
1359                // Y_1 (t)     = F[X_1(t)]
1360                //             =  (2.0 + 2t + 3t^2) - 2.0 * (2.0 + 3t + 4t^2)
1361                double y_1_1   = 2.0 - 6.0;
1362                double y_2_1   = 3.0 - 8.0;
1363                //
1364                ok  &= NearEqual(y1[0] , y_1_0, eps, eps);
1365                ok  &= NearEqual(y1[1] , y_1_1, eps, eps);
1366                ok  &= NearEqual(y2[0] , y_2_0, eps, eps);
1367                ok  &= NearEqual(y2[1] , y_2_1, eps, eps);
1368                //
1369                return ok;
1370        }
1371        // ---------------------------------------------------------------------
1372        // TanOp
1373        bool tan_op(void)
1374        {       bool ok = true;
1375                double eps = 10. * std::numeric_limits<double>::epsilon();
1376                size_t j;
1377
1378                // domain space vector
1379                size_t n = 1;
1380                CPPAD_TESTVECTOR(AD<double>) ax(n);
1381                ax[0] = 0.5;
1382
1383                // declare independent variables and starting recording
1384                CppAD::Independent(ax);
1385
1386                // range space vector
1387                size_t m = 1;
1388                CPPAD_TESTVECTOR(AD<double>) ay(m);
1389                ay[0] = tan( ax[0] );
1390
1391                // create f: x -> y and stop tape recording
1392                CppAD::ADFun<double> f(ax, ay);
1393
1394                // first order Taylor coefficients
1395                size_t r = 2, ell;
1396                CPPAD_TESTVECTOR(double) x1(r*n), y1;
1397                for(ell = 0; ell < r; ell++)
1398                {       for(j = 0; j < n; j++)
1399                                x1[ r * j + ell ] = double(j + ell + 1);
1400                }
1401                y1  = f.Forward(1, r, x1);
1402                ok &= y1.size() == r*m;
1403
1404                // secondorder Taylor coefficients
1405                CPPAD_TESTVECTOR(double) x2(r*n), y2;
1406                for(ell = 0; ell < r; ell++)
1407                {       for(j = 0; j < n; j++)
1408                                x2[ r * j + ell ] = double(j + ell + 2);
1409                }
1410                y2  = f.Forward(2, r, x2);
1411                ok &= y2.size() == r*m;
1412                //
1413                // Y_0  (t)    = F[X_0(t)]
1414                //             =  tan(0.5 + 1t + 2t^2)
1415                // Y_0' (t)    =  cos(0.5 + 1t + 2t^2)^(-2)*(1 + 4t)
1416                // Y_0''(0)    = 2*cos(0.5)^(-3)*sin(0.5) + 4*cos(0.5)^(-2)
1417                double sec_sq  = 1.0 / ( cos(0.5) * cos(0.5) );
1418                double y_1_0   = sec_sq;
1419                double y_2_0   = (2.0 * tan(0.5) + 4.0) * sec_sq / 2.0;
1420                //
1421                // Y_1  (t)    = F[X_1(t)]
1422                //             = tan(0.5 + 2t + 3t^2)
1423                // Y_1' (t)    = cos(0.5 + 2t + 3t^2)^(-2)*(2 + 6t)
1424                // Y_1''(0)    = 2*cos(0.5)^(-3)*sin(0.5)*2*2 + 6*cos(0.5)^(-2)
1425                double y_1_1   = 2.0 * sec_sq;
1426                double y_2_1   = (8.0 * tan(0.5) + 6.0) * sec_sq / 2.0;
1427                //
1428                ok  &= NearEqual(y1[0] , y_1_0, eps, eps);
1429                ok  &= NearEqual(y1[1] , y_1_1, eps, eps);
1430                ok  &= NearEqual(y2[0] , y_2_0, eps, eps);
1431                ok  &= NearEqual(y2[1] , y_2_1, eps, eps);
1432                //
1433                return ok;
1434        }
1435        // ---------------------------------------------------------------------
1436        // Usr*Op
1437        typedef CPPAD_TESTVECTOR(AD<double>) avector;
1438        void usr_algo(const avector& x, avector& z)
1439        {       z[0] = ( x[0] + x[1] ) / 2.0;
1440                z[1] = x[0] * x[1];
1441                z[2] = ( x[0] - x[1] ) / 2.0;
1442                return;
1443        }
1444        bool usr_op(void)
1445        {       bool ok = true;
1446                double eps = 10. * std::numeric_limits<double>::epsilon();
1447                size_t j;
1448
1449                // define checkpoint function
1450                size_t n = 2;
1451                avector ax(n), az(3);
1452                ax[0] = 0.5;
1453                ax[1] = 2.0;
1454                CppAD::checkpoint<double> usr_check("usr_check", usr_algo, ax, az);
1455
1456                // declare independent variables and starting recording
1457                CppAD::Independent(ax);
1458
1459                // record checkpoint function
1460                usr_check(ax, az);
1461
1462                // range space vector
1463                size_t m = 2;
1464                avector ay(m);
1465                ay[0] = az[0] + az[2]; // = ax[0]
1466                ay[1] = az[0] - az[2]; // = ax[1]
1467
1468                // create f: x -> y and stop tape recording
1469                CppAD::ADFun<double> f(ax, ay);
1470
1471                // first order Taylor coefficients
1472                size_t r = 2, ell;
1473                CPPAD_TESTVECTOR(double) x1(r*n), y1;
1474                for(ell = 0; ell < r; ell++)
1475                {       for(j = 0; j < n; j++)
1476                                x1[ r * j + ell ] = double(j + ell + 1);
1477                }
1478                y1  = f.Forward(1, r, x1);
1479                ok &= y1.size() == r*m;
1480
1481                // secondorder Taylor coefficients
1482                CPPAD_TESTVECTOR(double) x2(r*n), y2;
1483                for(ell = 0; ell < r; ell++)
1484                {       for(j = 0; j < n; j++)
1485                                x2[ r * j + ell ] = double(j + ell + 2);
1486                }
1487                y2  = f.Forward(2, r, x2);
1488                ok &= y2.size() == r*m;
1489                //
1490                // Y0_0 (t)    = 0.5 + 1t + 2t^2
1491                double y0_1_0  = 1.0;
1492                double y0_2_0  = 2.0;
1493                //
1494                // Y0_1 (t)    = 0.5 + 2t + 3t^2
1495                double y0_1_1  = 2.0;
1496                double y0_2_1  = 3.0;
1497                //
1498                // Y1_0 (t)    = 2.0 + 2t + 3t^2
1499                double y1_1_0  = 2.0;
1500                double y1_2_0  = 3.0;
1501                //
1502                // Y1_1 (t)    = 2.0 + 3t + 4t^2
1503                double y1_1_1  = 3.0;
1504                double y1_2_1  = 4.0;
1505                //
1506                ok  &= NearEqual(y1[0*r+0] , y0_1_0, eps, eps);
1507                ok  &= NearEqual(y1[1*r+0] , y1_1_0, eps, eps);
1508                ok  &= NearEqual(y1[0*r+1] , y0_1_1, eps, eps);
1509                ok  &= NearEqual(y1[1*r+1] , y1_1_1, eps, eps);
1510                //
1511                ok  &= NearEqual(y2[0*r+0] , y0_2_0, eps, eps);
1512                ok  &= NearEqual(y2[1*r+0] , y1_2_0, eps, eps);
1513                ok  &= NearEqual(y2[0*r+1] , y0_2_1, eps, eps);
1514                ok  &= NearEqual(y2[1*r+1] , y1_2_1, eps, eps);
1515                //
1516                return ok;
1517        }
1518        // ---------------------------------------------------------------------
1519        // Inverse functions assume the following already tested:
1520        // CosOp, SinOp, TanOp, ExpOp, MulvvOp, DivvpOp, AddpvOp
1521        //
1522        // AcosOp
1523        AD<double> acos_fun(const AD<double>& x)
1524        {       return acos( cos(x) ); }
1525        bool acos_op(void)
1526        {       return check_identity(acos_fun, 0.5); }
1527        //
1528        // AcoshOp
1529        AD<double> acosh_fun(const AD<double>& x)
1530        {       return acosh( cosh(x) ); }
1531        bool acosh_op(void)
1532        {       return check_identity(acosh_fun, 0.5); }
1533        //
1534        // AsinOp
1535        AD<double> asin_fun(const AD<double>& x)
1536        {       return asin( sin(x) ); }
1537        bool asin_op(void)
1538        {       return check_identity(asin_fun, 0.5); }
1539        //
1540        // AsinhOp
1541        AD<double> asinh_fun(const AD<double>& x)
1542        {       return asinh( sinh(x) ); }
1543        bool asinh_op(void)
1544        {       return check_identity(asinh_fun, 0.5); }
1545        //
1546        // AtanOp
1547        AD<double> atan_fun(const AD<double>& x)
1548        {       return atan( tan(x) ); }
1549        bool atan_op(void)
1550        {       return check_identity(atan_fun, 0.5); }
1551        //
1552        // AtanhOp
1553        AD<double> atanh_fun(const AD<double>& x)
1554        {       return atanh( tanh(x) ); }
1555        bool atanh_op(void)
1556        {       return check_identity(atanh_fun, 0.5); }
1557        //
1558        // LogOp
1559        AD<double> log_fun(const AD<double>& x)
1560        {       return log( exp(x) ); }
1561        bool log_op(void)
1562        {       return check_identity(log_fun, 0.5); }
1563        //
1564        // DivvvOp
1565        AD<double> divvv_fun(const AD<double>& x)
1566        {       return (x * x) / x; }
1567        bool divvv_op(void)
1568        {       return check_identity(divvv_fun, 0.5); }
1569        //
1570        // PowpvOp
1571        AD<double> powpv_fun(const AD<double>& x )
1572        {       return log( pow( exp(3.0) , x ) ) / 3.0; }
1573        bool powpv_op(void)
1574        {       return check_identity(powpv_fun, 0.5); }
1575        //
1576        // PowvpOp
1577        AD<double> powvp_fun(const AD<double>& x )
1578        {       return log( pow( exp(x) , 3.0 ) ) / 3.0; }
1579        bool powvp_op(void)
1580        {       return check_identity(powvp_fun, 0.5); }
1581        //
1582        // SqrtOp
1583        AD<double> sqrt_fun(const AD<double>& x )
1584        {       return sqrt( x * x ); }
1585        bool sqrt_op(void)
1586        {       return check_identity(sqrt_fun, 0.5); }
1587        //
1588        // SubvpOp
1589        AD<double> subvp_fun(const AD<double>& x )
1590        {       return 3.0 + ( x - 3.0 ); }
1591        bool subvp_op(void)
1592        {       return check_identity(subvp_fun, 0.5); }
1593        //
1594        // TanhOp
1595        AD<double> tanh_fun(const AD<double>& x )
1596        {       AD<double> z = tanh(x);
1597                return log( (1.0 + z) / (1.0 - z) ) / 2.0;
1598        }
1599        bool tanh_op(void)
1600        {       return check_identity(tanh_fun, 0.5); }
1601}
1602
1603bool forward_dir(void)
1604{       bool ok = true;
1605        //
1606        ok     &= abs_op();
1607        ok     &= acos_op();
1608        ok     &= acosh_op();
1609        ok     &= asin_op();
1610        ok     &= asinh_op();
1611        ok     &= atan_op();
1612        ok     &= atanh_op();
1613        ok     &= addpv_op();
1614        ok     &= addvv_op();
1615        ok     &= cexp_op();
1616        ok     &= cosh_op();
1617        ok     &= cos_op();
1618        ok     &= csum_op();
1619        ok     &= dis_op();
1620        ok     &= divpv_op();
1621        ok     &= divvp_op();
1622        ok     &= divvv_op();
1623        ok     &= exp_op();
1624        ok     &= load_op();
1625        ok     &= log_op();
1626        ok     &= mulpv_op();
1627        ok     &= par_op();
1628        ok     &= powpv_op();
1629        ok     &= powvp_op();
1630        ok     &= powvv_op();
1631        ok     &= sign_op();
1632        ok     &= sin_op();
1633        ok     &= sinh_op();
1634        ok     &= subpv_op();
1635        ok     &= subvp_op();
1636        ok     &= subvv_op();
1637        ok     &= sqrt_op();
1638        ok     &= tan_op();
1639        ok     &= tanh_op();
1640        ok     &= usr_op();
1641        //
1642        return ok;
1643}
Note: See TracBrowser for help on using the repository browser.