source: trunk/test_more/forward_dir.cpp @ 3520

Last change on this file since 3520 was 3305, checked in by bradbell, 6 years ago

forward_dir.cpp: remove another use of std::round.

  • Property svn:keywords set to Id
File size: 41.4 KB
Line 
1/* $Id: forward_dir.cpp 3305 2014-05-24 16:37:42Z bradbell $ */
2/* --------------------------------------------------------------------------
3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-14 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        // PowvvOp  (test assuming LogOp, ExpOp and DivvvOp are correct)
942        bool powvv_op(void)
943        {       bool ok = true;
944                double eps = 10. * std::numeric_limits<double>::epsilon();
945                size_t j;
946       
947                // domain space vector
948                size_t n = 2;
949                CPPAD_TESTVECTOR(AD<double>) ax(n);
950                ax[0] = 0.5; 
951                ax[1] = 2.0;
952       
953                // declare independent variables and starting recording
954                CppAD::Independent(ax);
955       
956                // range space vector
957                size_t m = 2;
958                CPPAD_TESTVECTOR(AD<double>) ay(m);
959                ay[0] = log( pow( exp(ax[0]) , ax[1] ) ) / ax[1] ;
960                ay[1] = log( pow( exp(ax[0]) , ax[1] ) ) / ax[0] ;
961       
962                // create f: x -> y and stop tape recording
963                CppAD::ADFun<double> f(ax, ay);
964       
965                // first order Taylor coefficients
966                size_t r = 2, ell;
967                CPPAD_TESTVECTOR(double) x1(r*n), y1;
968                for(ell = 0; ell < r; ell++)
969                {       for(j = 0; j < n; j++)
970                                x1[ r * j + ell ] = double(j + ell + 1);
971                }
972                y1  = f.Forward(1, r, x1);
973                ok &= y1.size() == r*m;
974               
975                // secondorder Taylor coefficients
976                CPPAD_TESTVECTOR(double) x2(r*n), y2;
977                for(ell = 0; ell < r; ell++)
978                {       for(j = 0; j < n; j++)
979                                x2[ r * j + ell ] = double(j + ell + 2);
980                }
981                y2  = f.Forward(2, r, x2);
982                ok &= y2.size() == r*m;
983                //
984                // Y0_0 (t)    = 0.5 + 1t + 2t^2
985                double y0_1_0  = 1.0;
986                double y0_2_0  = 2.0;
987                //
988                // Y0_1 (t)    = 0.5 + 2t + 3t^2
989                double y0_1_1  = 2.0;
990                double y0_2_1  = 3.0;
991                //
992                // Y1_0 (t)    = 2.0 + 2t + 3t^2
993                double y1_1_0  = 2.0;
994                double y1_2_0  = 3.0;
995                //
996                // Y1_1 (t)    = 2.0 + 3t + 4t^2
997                double y1_1_1  = 3.0;
998                double y1_2_1  = 4.0;
999                //
1000                ok  &= NearEqual(y1[0*r+0] , y0_1_0, eps, eps);
1001                ok  &= NearEqual(y1[1*r+0] , y1_1_0, eps, eps);
1002                ok  &= NearEqual(y1[0*r+1] , y0_1_1, eps, eps);
1003                ok  &= NearEqual(y1[1*r+1] , y1_1_1, eps, eps);
1004                //
1005                ok  &= NearEqual(y2[0*r+0] , y0_2_0, eps, eps);
1006                ok  &= NearEqual(y2[1*r+0] , y1_2_0, eps, eps);
1007                ok  &= NearEqual(y2[0*r+1] , y0_2_1, eps, eps);
1008                ok  &= NearEqual(y2[1*r+1] , y1_2_1, eps, eps);
1009                //
1010                return ok;
1011        }
1012        // ---------------------------------------------------------------------
1013        // SignOp (test assuming that MulvvOp is correct)
1014        bool sign_op(void)
1015        {       bool ok = true;
1016                double eps = 10. * std::numeric_limits<double>::epsilon();
1017                size_t j;
1018       
1019                // domain space vector
1020                size_t n = 1;
1021                CPPAD_TESTVECTOR(AD<double>) ax(n);
1022                ax[0] = 0.5; 
1023       
1024                // declare independent variables and starting recording
1025                CppAD::Independent(ax);
1026       
1027                // range space vector
1028                size_t m = 1;
1029                CPPAD_TESTVECTOR(AD<double>) ay(m);
1030                ay[0] = sign( ax[0] ) * ax[0];
1031       
1032                // create f: x -> y and stop tape recording
1033                CppAD::ADFun<double> f(ax, ay);
1034
1035                // zero order
1036                CPPAD_TESTVECTOR(double) x0(n), y0;
1037                x0[0] = -3.0;
1038                y0  = f.Forward(0, x0);
1039                ok &= y0.size() == m;
1040                ok &= NearEqual(y0[0], CppAD::abs(x0[0]), eps, eps); 
1041       
1042                // first order Taylor coefficients
1043                size_t r = 2, ell;
1044                CPPAD_TESTVECTOR(double) x1(r*n), y1;
1045                for(ell = 0; ell < r; ell++)
1046                {       for(j = 0; j < n; j++)
1047                                x1[ r * j + ell ] = double(j + ell + 1);
1048                }
1049                y1  = f.Forward(1, r, x1);
1050                ok &= y1.size() == r*m;
1051               
1052                // secondorder Taylor coefficients
1053                CPPAD_TESTVECTOR(double) x2(r*n), y2;
1054                for(ell = 0; ell < r; ell++)
1055                {       for(j = 0; j < n; j++)
1056                                x2[ r * j + ell ] = double(j + ell + 2);
1057                }
1058                y2  = f.Forward(2, r, x2);
1059                ok &= y2.size() == r*m;
1060                //
1061                //
1062                // Y_0 (t)     = F[X_0(t)]
1063                //             =  -(-3.0 + 1t + 2t^2)
1064                double y_1_0   = -1.0;
1065                double y_2_0   = -2.0;
1066                //
1067                // Y_1 (t)     = F[X_1(t)]
1068                //             =  -(-3.0 + 2t + 3t^2)
1069                double y_1_1   = -2.0;
1070                double y_2_1   = -3.0;
1071                //
1072                ok  &= NearEqual(y1[0] , y_1_0, eps, eps);
1073                ok  &= NearEqual(y1[1] , y_1_1, eps, eps);
1074                ok  &= NearEqual(y2[0] , y_2_0, eps, eps);
1075                ok  &= NearEqual(y2[1] , y_2_1, eps, eps);
1076                //
1077                return ok;
1078        }
1079
1080        // ---------------------------------------------------------------------
1081        // SinOp
1082        bool sin_op(void)
1083        {       bool ok = true;
1084                double eps = 10. * std::numeric_limits<double>::epsilon();
1085                size_t j;
1086
1087       
1088                // domain space vector
1089                size_t n = 1;
1090                CPPAD_TESTVECTOR(AD<double>) ax(n);
1091                ax[0] = 0.5; 
1092       
1093                // declare independent variables and starting recording
1094                CppAD::Independent(ax);
1095       
1096                // range space vector
1097                size_t m = 1;
1098                CPPAD_TESTVECTOR(AD<double>) ay(m);
1099                ay[0] = sin( ax[0] );
1100       
1101                // create f: x -> y and stop tape recording
1102                CppAD::ADFun<double> f(ax, ay);
1103       
1104                // first order Taylor coefficients
1105                size_t r = 2, ell;
1106                CPPAD_TESTVECTOR(double) x1(r*n), y1;
1107                for(ell = 0; ell < r; ell++)
1108                {       for(j = 0; j < n; j++)
1109                                x1[ r * j + ell ] = double(j + ell + 1);
1110                }
1111                y1  = f.Forward(1, r, x1);
1112                ok &= y1.size() == r*m;
1113               
1114                // secondorder Taylor coefficients
1115                CPPAD_TESTVECTOR(double) x2(r*n), y2;
1116                for(ell = 0; ell < r; ell++)
1117                {       for(j = 0; j < n; j++)
1118                                x2[ r * j + ell ] = double(j + ell + 2);
1119                }
1120                y2  = f.Forward(2, r, x2);
1121                ok &= y2.size() == r*m;
1122                //
1123                // Y_0  (t)    = F[X_0(t)]
1124                //             = sin( 0.5 + 1t + 2t^2 )
1125                // Y_0' (t)    = cos( 0.5 + 1t + 2t^2) * (1 + 4t)
1126                double y_1_0   = cos(0.5);
1127                double y_2_0   = ( cos(0.5) * 4.0 - sin(0.5) ) / 2.0;
1128                //
1129                // Y_1  (t)    = F[X_1(t)]
1130                //             = sin( 0.5 + 2t + 3t^2 )
1131                // Y_1' (t)    = cos( 0.5 + 2t + 3t^2) * (2 + 6t)
1132                double y_1_1   = cos(0.5) * 2.0;
1133                double y_2_1   = ( cos(0.5) * 6.0 - sin(0.5) * 4.0 ) / 2.0;
1134                //
1135                ok  &= NearEqual(y1[0] , y_1_0, eps, eps);
1136                ok  &= NearEqual(y1[1] , y_1_1, eps, eps);
1137                ok  &= NearEqual(y2[0] , y_2_0, eps, eps);
1138                ok  &= NearEqual(y2[1] , y_2_1, eps, eps);
1139                //
1140                return ok;
1141        }
1142        // ---------------------------------------------------------------------
1143        // SinhOp
1144        bool sinh_op(void)
1145        {       bool ok = true;
1146                double eps = 10. * std::numeric_limits<double>::epsilon();
1147                size_t j;
1148
1149       
1150                // domain space vector
1151                size_t n = 1;
1152                CPPAD_TESTVECTOR(AD<double>) ax(n);
1153                ax[0] = 0.5; 
1154       
1155                // declare independent variables and starting recording
1156                CppAD::Independent(ax);
1157       
1158                // range space vector
1159                size_t m = 1;
1160                CPPAD_TESTVECTOR(AD<double>) ay(m);
1161                ay[0] = sinh( ax[0] );
1162       
1163                // create f: x -> y and stop tape recording
1164                CppAD::ADFun<double> f(ax, ay);
1165       
1166                // first order Taylor coefficients
1167                size_t r = 2, ell;
1168                CPPAD_TESTVECTOR(double) x1(r*n), y1;
1169                for(ell = 0; ell < r; ell++)
1170                {       for(j = 0; j < n; j++)
1171                                x1[ r * j + ell ] = double(j + ell + 1);
1172                }
1173                y1  = f.Forward(1, r, x1);
1174                ok &= y1.size() == r*m;
1175               
1176                // secondorder Taylor coefficients
1177                CPPAD_TESTVECTOR(double) x2(r*n), y2;
1178                for(ell = 0; ell < r; ell++)
1179                {       for(j = 0; j < n; j++)
1180                                x2[ r * j + ell ] = double(j + ell + 2);
1181                }
1182                y2  = f.Forward(2, r, x2);
1183                ok &= y2.size() == r*m;
1184                //
1185                // Y_0  (t)    = F[X_0(t)]
1186                //             = sinh( 0.5 + 1t + 2t^2 )
1187                // Y_0' (t)    = cosh( 0.5 + 1t + 2t^2) * (1 + 4t)
1188                double y_1_0   = cosh(0.5);
1189                double y_2_0   = ( cosh(0.5) * 4.0 + sinh(0.5) ) / 2.0;
1190                //
1191                // Y_1  (t)    = F[X_1(t)]
1192                //             = sinh( 0.5 + 2t + 3t^2 )
1193                // Y_1' (t)    = cosh( 0.5 + 2t + 3t^2) * (2 + 6t)
1194                double y_1_1   = cosh(0.5) * 2.0;
1195                double y_2_1   = ( cosh(0.5) * 6.0 + sinh(0.5) * 4.0 ) / 2.0;
1196                //
1197                ok  &= NearEqual(y1[0] , y_1_0, eps, eps);
1198                ok  &= NearEqual(y1[1] , y_1_1, eps, eps);
1199                ok  &= NearEqual(y2[0] , y_2_0, eps, eps);
1200                ok  &= NearEqual(y2[1] , y_2_1, eps, eps);
1201                //
1202                return ok;
1203        }
1204        // ---------------------------------------------------------------------
1205        // SubpvOp
1206        bool subpv_op(void)
1207        {       bool ok = true;
1208                double eps = 10. * std::numeric_limits<double>::epsilon();
1209                size_t j;
1210       
1211                // domain space vector
1212                size_t n = 1;
1213                CPPAD_TESTVECTOR(AD<double>) ax(n);
1214                ax[0] = 0.5; 
1215       
1216                // declare independent variables and starting recording
1217                CppAD::Independent(ax);
1218       
1219                // range space vector
1220                size_t m = 1;
1221                CPPAD_TESTVECTOR(AD<double>) ay(m);
1222                ay[0] = 2.0 - ax[0];
1223       
1224                // create f: x -> y and stop tape recording
1225                CppAD::ADFun<double> f(ax, ay);
1226       
1227                // first order Taylor coefficients
1228                size_t r = 2, ell;
1229                CPPAD_TESTVECTOR(double) x1(r*n), y1;
1230                for(ell = 0; ell < r; ell++)
1231                {       for(j = 0; j < n; j++)
1232                                x1[ r * j + ell ] = double(j + ell + 1);
1233                }
1234                y1  = f.Forward(1, r, x1);
1235                ok &= y1.size() == r*m;
1236               
1237                // secondorder Taylor coefficients
1238                CPPAD_TESTVECTOR(double) x2(r*n), y2;
1239                for(ell = 0; ell < r; ell++)
1240                {       for(j = 0; j < n; j++)
1241                                x2[ r * j + ell ] = double(j + ell + 3);
1242                }
1243                y2  = f.Forward(2, r, x2);
1244                ok &= y2.size() == r*m;
1245                //
1246                // Y_0 (t)     = F[X_0(t)]
1247                //             =  2.0 - (0.5 + 1t + 3t^2)/2.0
1248                double y_1_0   = - 1.0;
1249                double y_2_0   = - 3.0;
1250                //
1251                // Y_1 (t)     = F[X_1(t)]
1252                //             =  3.0 - (0.5 + 2t + 4t^2)/2.0
1253                double y_1_1   = - 2.0;
1254                double y_2_1   = - 4.0;
1255                //
1256                ok  &= NearEqual(y1[0] , y_1_0, eps, eps);
1257                ok  &= NearEqual(y1[1] , y_1_1, eps, eps);
1258                ok  &= NearEqual(y2[0] , y_2_0, eps, eps);
1259                ok  &= NearEqual(y2[1] , y_2_1, eps, eps);
1260                //
1261                return ok;
1262        }
1263        // ---------------------------------------------------------------------
1264        // SubvvOp
1265        bool subvv_op(void)
1266        {       bool ok = true;
1267                double eps = 10. * std::numeric_limits<double>::epsilon();
1268                size_t j;
1269       
1270                // domain space vector
1271                size_t n = 2;
1272                CPPAD_TESTVECTOR(AD<double>) ax(n);
1273                ax[0] = 0.5; 
1274                ax[1] = 2.0;
1275       
1276                // declare independent variables and starting recording
1277                CppAD::Independent(ax);
1278       
1279                // range space vector
1280                size_t m = 1;
1281                CPPAD_TESTVECTOR(AD<double>) ay(m);
1282                ay[0] = ax[0] - 2.0 * ax[1];
1283       
1284                // create f: x -> y and stop tape recording
1285                CppAD::ADFun<double> f(ax, ay);
1286       
1287                // first order Taylor coefficients
1288                size_t r = 2, ell;
1289                CPPAD_TESTVECTOR(double) x1(r*n), y1;
1290                for(ell = 0; ell < r; ell++)
1291                {       for(j = 0; j < n; j++)
1292                                x1[ r * j + ell ] = double(j + ell + 1);
1293                }
1294                y1  = f.Forward(1, r, x1);
1295                ok &= y1.size() == r*m;
1296               
1297                // secondorder Taylor coefficients
1298                CPPAD_TESTVECTOR(double) x2(r*n), y2;
1299                for(ell = 0; ell < r; ell++)
1300                {       for(j = 0; j < n; j++)
1301                                x2[ r * j + ell ] = double(j + ell + 2);
1302                }
1303                y2  = f.Forward(2, r, x2);
1304                ok &= y2.size() == r*m;
1305                //
1306                // Y_0 (t)     = F[X_0(t)]
1307                //             =  (0.5 + 1t + 2t^2) - 2.0 * (2.0 + 2t + 3t^2)
1308                double y_1_0   = 1.0 - 4.0;
1309                double y_2_0   = 2.0 - 6.0;
1310                //
1311                // Y_1 (t)     = F[X_1(t)]
1312                //             =  (2.0 + 2t + 3t^2) - 2.0 * (2.0 + 3t + 4t^2)
1313                double y_1_1   = 2.0 - 6.0;
1314                double y_2_1   = 3.0 - 8.0;
1315                //
1316                ok  &= NearEqual(y1[0] , y_1_0, eps, eps);
1317                ok  &= NearEqual(y1[1] , y_1_1, eps, eps);
1318                ok  &= NearEqual(y2[0] , y_2_0, eps, eps);
1319                ok  &= NearEqual(y2[1] , y_2_1, eps, eps);
1320                //
1321                return ok;
1322        }
1323        // ---------------------------------------------------------------------
1324        // TanOp
1325        bool tan_op(void)
1326        {       bool ok = true;
1327                double eps = 10. * std::numeric_limits<double>::epsilon();
1328                size_t j;
1329       
1330                // domain space vector
1331                size_t n = 1;
1332                CPPAD_TESTVECTOR(AD<double>) ax(n);
1333                ax[0] = 0.5; 
1334       
1335                // declare independent variables and starting recording
1336                CppAD::Independent(ax);
1337       
1338                // range space vector
1339                size_t m = 1;
1340                CPPAD_TESTVECTOR(AD<double>) ay(m);
1341                ay[0] = tan( ax[0] );
1342       
1343                // create f: x -> y and stop tape recording
1344                CppAD::ADFun<double> f(ax, ay);
1345       
1346                // first order Taylor coefficients
1347                size_t r = 2, ell;
1348                CPPAD_TESTVECTOR(double) x1(r*n), y1;
1349                for(ell = 0; ell < r; ell++)
1350                {       for(j = 0; j < n; j++)
1351                                x1[ r * j + ell ] = double(j + ell + 1);
1352                }
1353                y1  = f.Forward(1, r, x1);
1354                ok &= y1.size() == r*m;
1355               
1356                // secondorder Taylor coefficients
1357                CPPAD_TESTVECTOR(double) x2(r*n), y2;
1358                for(ell = 0; ell < r; ell++)
1359                {       for(j = 0; j < n; j++)
1360                                x2[ r * j + ell ] = double(j + ell + 2);
1361                }
1362                y2  = f.Forward(2, r, x2);
1363                ok &= y2.size() == r*m;
1364                //
1365                // Y_0  (t)    = F[X_0(t)]
1366                //             =  tan(0.5 + 1t + 2t^2)
1367                // Y_0' (t)    =  cos(0.5 + 1t + 2t^2)^(-2)*(1 + 4t)
1368                // Y_0''(0)    = 2*cos(0.5)^(-3)*sin(0.5) + 4*cos(0.5)^(-2)
1369                double sec_sq  = 1.0 / ( cos(0.5) * cos(0.5) );
1370                double y_1_0   = sec_sq;
1371                double y_2_0   = (2.0 * tan(0.5) + 4.0) * sec_sq / 2.0;
1372                //
1373                // Y_1  (t)    = F[X_1(t)]
1374                //             = tan(0.5 + 2t + 3t^2)
1375                // Y_1' (t)    = cos(0.5 + 2t + 3t^2)^(-2)*(2 + 6t)
1376                // Y_1''(0)    = 2*cos(0.5)^(-3)*sin(0.5)*2*2 + 6*cos(0.5)^(-2)
1377                double y_1_1   = 2.0 * sec_sq;
1378                double y_2_1   = (8.0 * tan(0.5) + 6.0) * sec_sq / 2.0;
1379                //
1380                ok  &= NearEqual(y1[0] , y_1_0, eps, eps);
1381                ok  &= NearEqual(y1[1] , y_1_1, eps, eps);
1382                ok  &= NearEqual(y2[0] , y_2_0, eps, eps);
1383                ok  &= NearEqual(y2[1] , y_2_1, eps, eps);
1384                //
1385                return ok;
1386        }
1387        // ---------------------------------------------------------------------
1388        // Usr*Op 
1389        typedef CPPAD_TESTVECTOR(AD<double>) avector;
1390        void usr_algo(const avector& x, avector& z)
1391        {       z[0] = ( x[0] + x[1] ) / 2.0;
1392                z[1] = x[0] * x[1];
1393                z[2] = ( x[0] - x[1] ) / 2.0;
1394                return;
1395        }
1396        bool usr_op(void)
1397        {       bool ok = true;
1398                double eps = 10. * std::numeric_limits<double>::epsilon();
1399                size_t j;
1400
1401                // define checkpoint function
1402                size_t n = 2;
1403                avector ax(n), az(3);
1404                ax[0] = 0.5; 
1405                ax[1] = 2.0;
1406                CppAD::checkpoint<double> usr_check("usr_check", usr_algo, ax, az);
1407
1408                // declare independent variables and starting recording
1409                CppAD::Independent(ax);
1410
1411                // record checkpoint function
1412                usr_check(ax, az);
1413       
1414                // range space vector
1415                size_t m = 2;
1416                avector ay(m);
1417                ay[0] = az[0] + az[2]; // = ax[0]
1418                ay[1] = az[0] - az[2]; // = ax[1]
1419
1420                // create f: x -> y and stop tape recording
1421                CppAD::ADFun<double> f(ax, ay);
1422       
1423                // first order Taylor coefficients
1424                size_t r = 2, ell;
1425                CPPAD_TESTVECTOR(double) x1(r*n), y1;
1426                for(ell = 0; ell < r; ell++)
1427                {       for(j = 0; j < n; j++)
1428                                x1[ r * j + ell ] = double(j + ell + 1);
1429                }
1430                y1  = f.Forward(1, r, x1);
1431                ok &= y1.size() == r*m;
1432               
1433                // secondorder Taylor coefficients
1434                CPPAD_TESTVECTOR(double) x2(r*n), y2;
1435                for(ell = 0; ell < r; ell++)
1436                {       for(j = 0; j < n; j++)
1437                                x2[ r * j + ell ] = double(j + ell + 2);
1438                }
1439                y2  = f.Forward(2, r, x2);
1440                ok &= y2.size() == r*m;
1441                //
1442                // Y0_0 (t)    = 0.5 + 1t + 2t^2
1443                double y0_1_0  = 1.0;
1444                double y0_2_0  = 2.0;
1445                //
1446                // Y0_1 (t)    = 0.5 + 2t + 3t^2
1447                double y0_1_1  = 2.0;
1448                double y0_2_1  = 3.0;
1449                //
1450                // Y1_0 (t)    = 2.0 + 2t + 3t^2
1451                double y1_1_0  = 2.0;
1452                double y1_2_0  = 3.0;
1453                //
1454                // Y1_1 (t)    = 2.0 + 3t + 4t^2
1455                double y1_1_1  = 3.0;
1456                double y1_2_1  = 4.0;
1457                //
1458                ok  &= NearEqual(y1[0*r+0] , y0_1_0, eps, eps);
1459                ok  &= NearEqual(y1[1*r+0] , y1_1_0, eps, eps);
1460                ok  &= NearEqual(y1[0*r+1] , y0_1_1, eps, eps);
1461                ok  &= NearEqual(y1[1*r+1] , y1_1_1, eps, eps);
1462                //
1463                ok  &= NearEqual(y2[0*r+0] , y0_2_0, eps, eps);
1464                ok  &= NearEqual(y2[1*r+0] , y1_2_0, eps, eps);
1465                ok  &= NearEqual(y2[0*r+1] , y0_2_1, eps, eps);
1466                ok  &= NearEqual(y2[1*r+1] , y1_2_1, eps, eps);
1467                //
1468                return ok;
1469        }
1470        // ---------------------------------------------------------------------
1471        // Inverse functions assume the following already tested:
1472        // CosOp, SinOp, TanOp, ExpOp, MulvvOp, DivvpOp, AddpvOp
1473        //
1474        // AcosOp
1475        AD<double> acos_fun(const AD<double>& x)
1476        {       return acos( cos(x) ); }
1477        bool acos_op(void)
1478        {       return check_identity(acos_fun, 0.5); }
1479        //
1480        // AsinOp
1481        AD<double> asin_fun(const AD<double>& x)
1482        {       return asin( sin(x) ); }
1483        bool asin_op(void)
1484        {       return check_identity(asin_fun, 0.5); }
1485        //
1486        // AtanOp
1487        AD<double> atan_fun(const AD<double>& x)
1488        {       return atan( tan(x) ); }
1489        bool atan_op(void)
1490        {       return check_identity(atan_fun, 0.5); }
1491        //
1492        // LogOp
1493        AD<double> log_fun(const AD<double>& x)
1494        {       return log( exp(x) ); }
1495        bool log_op(void)
1496        {       return check_identity(log_fun, 0.5); }
1497        //
1498        // DivvvOp
1499        AD<double> divvv_fun(const AD<double>& x)
1500        {       return (x * x) / x; }
1501        bool divvv_op(void)
1502        {       return check_identity(divvv_fun, 0.5); }
1503        //
1504        // PowpvOp
1505        AD<double> powpv_fun(const AD<double>& x )
1506        {       return log( pow( exp(3.0) , x ) ) / 3.0; }
1507        bool powpv_op(void)
1508        {       return check_identity(powpv_fun, 0.5); }
1509        //
1510        // PowvpOp
1511        AD<double> powvp_fun(const AD<double>& x )
1512        {       return log( pow( exp(x) , 3.0 ) ) / 3.0; }
1513        bool powvp_op(void)
1514        {       return check_identity(powvp_fun, 0.5); }
1515        //
1516        // SqrtOp
1517        AD<double> sqrt_fun(const AD<double>& x )
1518        {       return sqrt( x * x ); }
1519        bool sqrt_op(void)
1520        {       return check_identity(sqrt_fun, 0.5); }
1521        //
1522        // SubvpOp
1523        AD<double> subvp_fun(const AD<double>& x )
1524        {       return 3.0 + ( x - 3.0 ); }
1525        bool subvp_op(void)
1526        {       return check_identity(subvp_fun, 0.5); }
1527        //
1528        // TanhOp
1529        AD<double> tanh_fun(const AD<double>& x )
1530        {       AD<double> z = tanh(x);
1531                return log( (1.0 + z) / (1.0 - z) ) / 2.0;
1532        }
1533        bool tanh_op(void)
1534        {       return check_identity(tanh_fun, 0.5); }
1535}
1536
1537bool forward_dir(void)
1538{       bool ok = true;
1539        //
1540        ok     &= abs_op();
1541        ok     &= acos_op();
1542        ok     &= asin_op();
1543        ok     &= atan_op();
1544        ok     &= addpv_op();
1545        ok     &= addvv_op();
1546        ok     &= cexp_op();
1547        ok     &= cosh_op();
1548        ok     &= cos_op();
1549        ok     &= csum_op();
1550        ok     &= dis_op();
1551        ok     &= divpv_op();
1552        ok     &= divvp_op();
1553        ok     &= divvv_op();
1554        ok     &= exp_op();
1555        ok     &= load_op();
1556        ok     &= log_op();
1557        ok     &= mulpv_op();
1558        ok     &= powpv_op();
1559        ok     &= powvp_op();
1560        ok     &= powvv_op();
1561        ok     &= sign_op();
1562        ok     &= sin_op();
1563        ok     &= sinh_op();
1564        ok     &= subpv_op();
1565        ok     &= subvp_op();
1566        ok     &= subvv_op();
1567        ok     &= sqrt_op();
1568        ok     &= tan_op();
1569        ok     &= tanh_op();
1570        ok     &= usr_op();
1571        //
1572        return ok;
1573}
Note: See TracBrowser for help on using the repository browser.