source: trunk/test_more/optimize.cpp @ 3638

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

merge to branch: trunk
from repository: https://github.com/coin-or/CppAD
start hash code: 3249a4c6f8216de001805d3e22bcdec63fed2d25
end hash code: 22c1453e46e3eb16e9a96d679a0a1fd459c849c0

commit 22c1453e46e3eb16e9a96d679a0a1fd459c849c0
Author: Brad Bell <bradbell@…>
Date: Tue Feb 10 07:19:54 2015 -0700

Fix bug in optimization and sparsity calculations that include the C++11 erf function.

commit df8602cb2129722c0b373fa467a76289769c026e
Author: Brad Bell <bradbell@…>
Date: Mon Feb 9 20:25:19 2015 -0700

erf.sh: Bug report by Micheal Braun.

  • Property svn:keywords set to Id
File size: 43.0 KB
Line 
1/* $Id: optimize.cpp 3638 2015-02-10 15:04:04Z 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// 2DO: Test that optimize.hpp use of base_atomic<Base>::rev_sparse_jac works.
13
14# include <limits>
15# include <cppad/cppad.hpp>
16
17
18namespace {
19        // accuracy for almost equal checks
20        double eps = 10. * std::numeric_limits<double>::epsilon();
21        using CppAD::NearEqual;
22
23        // note this enum type is not part of the API (but its values are)
24        CppAD::atomic_base<double>::option_enum atomic_sparsity_option;
25        //
26        // ----------------------------------------------------------------
27        // Test nested conditional expressions.
28        int nested_cond_exp(void)
29        {       bool ok = true;
30                using CppAD::AD;
31                using CppAD::vector;
32       
33                // independent variable vector
34                vector< AD<double> > ax(2), ay(1);
35                ax[0] = 1.0;
36                ax[1] = 2.0;
37                Independent(ax);
38       
39                // first conditional expression
40                AD<double> ac1 = CondExpLe(ax[0], ax[1], 2.0 * ax[0], 3.0 * ax[1] );
41       
42                // second conditional expression
43                AD<double> ac2 = CondExpGe(ax[0], ax[1], 4.0 * ax[0], 5.0 * ax[1] );
44       
45                // third conditional expression
46                AD<double> ac3 = CondExpLt(ax[0], ax[1], 6.0 * ac1, 7.0 * ac2 );
47       
48                // create function object f : ax -> ay
49                ay[0] = ac3;
50                CppAD::ADFun<double> f(ax, ay);
51       
52                // now optimize the operation sequence
53                f.optimize();
54       
55                // now zero order forward
56                vector<double> x(2), y(1);
57                for(size_t i = 0; i < 3; i++)
58                {       x[0] = 1.0 - double(i);
59                        x[1] = - x[0];
60                        y    = f.Forward(0, x);
61                        //
62                        // first conditional expression
63                        double c1;
64                        if( x[0] <= x[1] )
65                                c1 = 2.0 * x[0];
66                        else    c1 = 3.0 * x[1];
67                        //
68                        // second conditional expression
69                        double c2;
70                        if( x[0] >= x[1] )
71                                c2 = 4.0 * x[0];
72                        else    c2 = 5.0 * x[1];
73       
74                        // third conditional expression
75                        double c3;
76                        if( x[0] < x[1] )
77                                c3 = 6.0 * c1;
78                        else    c3 = 7.0 * c2;
79       
80                        ok &= y[0] == c3;
81                }
82                return ok;
83        }
84        // ----------------------------------------------------------------
85        // Test for bug where checkpoint function did not depend on
86        // the operands in the logical comparison because of the CondExp
87        // sparsity pattern.
88        void j_algo( 
89                const CppAD::vector< CppAD::AD<double> >& ax ,
90                      CppAD::vector< CppAD::AD<double> >& ay )
91        {       ay[0] = CondExpGt(ax[0], ax[1], ax[2], ax[3]); }
92       
93        bool atomic_cond_exp_sparsity(void)
94        {       bool ok = true;
95                using CppAD::AD;
96                using CppAD::vector;
97       
98                // Create a checkpoint version of the function g
99                vector< AD<double> > au(4), av(1);
100                for(size_t i = 0; i < 4; i++)
101                        au[i] = AD<double>(i);
102                CppAD::checkpoint<double> j_check("j_check", j_algo, au, av);
103       
104                // independent variable vector
105                vector< AD<double> > ax(2), ay(1);
106                ax[0] = 1.;
107                ax[1] = 1.;
108                Independent(ax);
109       
110                // call atomic function that does not get used
111                for(size_t i = 0; i < 4; i++) 
112                        au[i] = ax[0] + AD<double>(i + 1) * ax[1];
113                j_check(au, ay);
114       
115                // create function object f : ax -> ay
116                CppAD::ADFun<double> f(ax, ay);
117
118       
119                // now optimize the operation sequence
120                j_check.option( atomic_sparsity_option );
121                f.optimize();
122       
123                // check result where true case is used; i.e., au[0] > au[1]
124                vector<double> x(2), y(1);
125                x[0] = 1.;
126                x[1] = -1;
127                y    = f.Forward(0, x);
128                ok  &= y[0] == x[0] + double(3) * x[1];
129               
130       
131                // check result where false case is used; i.e., au[0] <= au[1]
132                x[0] = 1.;
133                x[1] = 1;
134                y    = f.Forward(0, x);
135                ok  &= y[0] == x[0] + double(4) * x[1];
136               
137                return ok;
138        }
139        // -------------------------------------------------------------------
140        // Test conditional optimizing out call to an atomic function call
141        void k_algo(
142                const CppAD::vector< CppAD::AD<double> >& x ,
143                      CppAD::vector< CppAD::AD<double> >& y )
144        {       y[0] = x[0] + x[1]; }
145
146        void h_algo(
147                const CppAD::vector< CppAD::AD<double> >& x ,
148                      CppAD::vector< CppAD::AD<double> >& y )
149        {       y[0] = x[0] - x[1]; }
150
151        bool atomic_cond_exp(void)
152        {       bool ok = true;
153                typedef CppAD::vector< CppAD::AD<double> > ADVector;
154
155                // Create a checkpoint version of the function g
156                ADVector ax(2), ag(1), ah(1), ay(1);
157                ax[0] = 0.;
158                ax[1] = 1.;
159                CppAD::checkpoint<double> k_check("k_check", k_algo, ax, ag);
160                CppAD::checkpoint<double> h_check("h_check", h_algo, ax, ah);
161
162                // independent variable vector
163                Independent(ax);
164
165                // atomic function calls that get conditionally used
166                k_check(ax, ag);
167                h_check(ax, ah);
168
169                // conditional expression
170                ay[0] = CondExpLt(ax[0], ax[1], ag[0], ah[0]); 
171       
172                // create function object f : ax -> ay
173                CppAD::ADFun<double> f;
174                f.Dependent(ax, ay);
175       
176                // use zero order to evaluate when condition is true
177                CppAD::vector<double>  x(2), dx(2);
178                CppAD::vector<double>  y(1), dy(1), w(1);
179                x[0] = 3.;
180                x[1] = 4.;
181                y    = f.Forward(0, x);
182                ok  &= y[0] == x[0] + x[1];
183
184                // before optimize
185                k_check.option( atomic_sparsity_option );
186                h_check.option( atomic_sparsity_option );
187                ok  &= f.number_skip() == 0;
188
189                // now optimize the operation sequence
190                f.optimize();
191
192                // optimized zero order forward when condition is false
193                x[0] = 4.;
194                x[1] = 3.;
195                y    = f.Forward(0, x);
196                ok  &= y[0] == x[0] - x[1];
197
198                // after optimize can skip either call to g or call to h
199                ok  &= f.number_skip() == 1;
200
201                // optimized first order forward
202                dx[0] = 2.;
203                dx[1] = 1.;
204                dy    = f.Forward(1, dx);
205                ok   &= dy[0] == dx[0] - dx[1];
206
207                // optimized first order reverse
208                w[0]  = 1.;
209                dx    = f.Reverse(1, w);
210                ok   &= dx[0] == 1.;
211                ok   &= dx[1] == -1.;
212       
213                return ok;
214        }
215        // -------------------------------------------------------------------
216        // Test of optimizing out arguments to an atomic function
217        void g_algo( 
218                const CppAD::vector< CppAD::AD<double> >& ax ,
219                      CppAD::vector< CppAD::AD<double> >& ay )
220        {       ay = ax; }
221
222        bool atomic_no_used(void)
223        {       bool ok = true;
224                using CppAD::AD;
225                using CppAD::vector;
226       
227                // Create a checkpoint version of the function g
228                vector< AD<double> > ax(2), ay(2), az(1);
229                ax[0] = 0.;
230                ax[1] = 1.;
231                CppAD::checkpoint<double> g_check("g_check", g_algo, ax, ay);
232       
233                // independent variable vector
234                Independent(ax);
235       
236                // call atomic function that does not get used
237                g_check(ax, ay);
238       
239                // conditional expression
240                az[0] = CondExpLt(ax[0], ax[1], ax[0] + ax[1], ax[0] - ax[1]); 
241               
242                // create function object f : ax -> az
243                CppAD::ADFun<double> f(ax, az);
244
245                // number of variables before optimization
246                // (include ay[0] and ay[1])
247                size_t n_before = f.size_var();
248               
249                // now optimize the operation sequence
250                g_check.option( atomic_sparsity_option );
251                f.optimize();
252
253                // number of variables after optimization
254                // (does not include ay[0] and ay[1])
255                size_t n_after = f.size_var();
256                ok            &= n_after + 2 == n_before;
257       
258                // check optimization works ok
259                vector<double> x(2), z(1);
260                x[0] = 4.;
261                x[1] = 3.;
262                z    = f.Forward(0, x);
263                ok  &= z[0] == x[0] - x[1];
264               
265                return ok;
266        }
267        bool atomic_arguments(void)
268        {       bool ok = true;
269                using CppAD::AD;
270                using CppAD::vector;
271                vector< AD<double> > au(2), aw(2), ax(2), ay(1);
272
273                // create atomic function corresponding to g_algo
274                au[0] = 1.0;
275                au[1] = 2.0;
276                CppAD::checkpoint<double> g_check("g_algo", g_algo, au, ax);
277
278                // start recording a new function
279                CppAD::Independent(ax);
280
281                // now use g_check during the recording
282                au[0] = ax[0] + ax[1]; // this argument requires a new variable
283                au[1] = ax[0] - ax[1]; // this argument also requires a new variable
284                g_check(au, aw);
285
286                // now create f(x) = x_0 - x_1
287                ay[0] = aw[0];
288                CppAD::ADFun<double> f(ax, ay);
289
290                // number of variables before optimization
291                size_t n_before = f.size_var();
292 
293                // now optimize f so that the calculation of au[1] is removed
294                g_check.option( atomic_sparsity_option );
295                f.optimize();
296
297                // check difference in number of variables
298                size_t n_after = f.size_var();
299                ok &= n_before == n_after + 1;
300
301                // now compute and check a forward mode calculation
302                vector<double> x(2), y(1);
303                x[0] = 5.0;
304                x[1] = 6.0;
305                y    = f.Forward(0, x);
306                ok  &= (y[0] == x[0] + x[1]); 
307
308                return ok;
309        }
310
311        // -------------------------------------------------------------------
312        // Test the reverse dependency analysis optimization
313        template <class Vector>
314        void depend_fun
315        (const Vector& x, Vector& y, size_t& original, size_t& opt)
316        {       typedef typename Vector::value_type Scalar;
317                Scalar not_used;
318                Scalar one(1), two(2), three(3), four(4);
319
320                // independent variable and phantom at beginning
321                original = 1 + x.size();
322                opt      = 1 + x.size();
323
324                // unary operator where operand is arg[0]
325                // (note that sin corresponds to two tape variables)
326                not_used = CppAD::abs(x[0]); 
327                y[0]     = sin(x[0]);
328                original += 3;
329                opt      += 2;
330
331                // binary operator where left operand is a variable
332                // and right operand is a parameter
333                not_used = not_used + 2.;
334                y[1]     = x[1] * 3.;
335                original += 2;
336                opt      += 1;
337
338                // binary operator where left operand is a parameter
339                // and right operation is a variable
340                not_used = 2. - not_used;
341                y[2]     = 3. / x[2];
342                original += 2;
343                opt      += 1;
344
345                // binary operator where both operands are variables
346                not_used  = x[3] - not_used;
347                y[3]      = x[3] / x[2];
348                original += 2;
349                opt      += 1;
350
351                // conditional expression that will be optimized out
352                not_used = CppAD::CondExpLt(x[0], x[1], x[2], x[3]) + not_used;
353                y[4]     = CppAD::CondExpLt(x[4], one, two, three);
354                original += 3;
355                opt      += 1;
356
357                // y[5] does not depend on the value of not_used.
358                // Make sure a parameter, corresponding to a dependent variable,
359                // is not optimized out of the operation sequence.
360                y[5]      = 0.0 * not_used;
361                original += 1;
362                opt      += 1;
363
364                // Wwe do not use the argument x[5], to
365                // make sure it is not optimized out.
366
367                return;
368        }
369
370        bool depend_one(void)
371        {       // Test all except for VecAD operations
372                bool ok = true;
373                using CppAD::AD;
374                size_t original;
375                size_t opt;
376                size_t i, j;
377       
378                // domain space vector
379                size_t n  = 6;
380                CppAD::vector< AD<double> > X(n);
381                for(j = 0; j < n; j++)
382                        X[j] = 1. / double(j + 1); 
383       
384                // declare independent variables and start tape recording
385                CppAD::Independent(X);
386       
387                // range space vector
388                size_t m = n;
389                CppAD::vector< AD<double> > Y(m);
390                depend_fun(X, Y, original, opt);
391       
392                // create f: X -> Y and stop tape recording
393                CppAD::ADFun<double> F;
394                F.Dependent(X, Y); 
395       
396                CppAD::vector<double> x(n), y(m), check(m);
397                for(j = 0; j < n; j++)
398                        x[j] = Value(X[j]);
399                y = F.Forward(0, x);
400                depend_fun(x, check, original, opt);
401                for(i = 0; i < m; i++)
402                        ok &= NearEqual(y[i], check[i], eps, eps);
403       
404                // Check size before optimization
405                ok &= F.size_var() == original;
406       
407                // Optimize the operation sequence
408                F.optimize();
409       
410                // Check size after optimization
411                ok &= F.size_var() == opt;
412       
413                // check result now
414                // (should have already been checked if NDEBUG not defined)
415                y = F.Forward(0, x);
416                for(i = 0; i < m; i++)
417                        ok &= NearEqual(y[i], check[i], eps, eps);
418       
419                return ok;
420        }
421
422        bool depend_two(void)
423        {       // Test VecAD operations
424                bool ok = true;
425                using CppAD::AD;
426                size_t i, j;
427
428                // domain space vector
429                size_t n  = 2;
430                CppAD::vector< AD<double> > X(n);
431                for(j = 0; j < n; j++)
432                        X[j] = double(j); 
433
434                // range space vector
435                size_t m = 3;
436                CppAD::vector< AD<double> > Y(m);
437
438                CppAD::VecAD<double> U(m);
439                CppAD::VecAD<double> V(n);
440                for(i = 0; i < m; i++)
441                        U[i] = 0;
442                for(j = 0; j < n; j++)
443                        V[j] = 0;
444       
445                // declare independent variables and start tape recording
446                CppAD::Independent(X);
447
448                // first vecad vector that is a variable
449                U[ X[0] ] = X[1];
450
451                // second vecad vector that is a variable
452                V[ X[0] ] = X[1];
453
454                // Make dependency for vecad vectors different that for
455                // variables because original code used worng dependency info.
456                // Y does not depend on the first variable in the tape; i.e.
457                // the one corresponding to the BeginOp. So make it depend
458                // on the first vecad vector in the tape.
459                for(i = 0; i < m; i++)
460                {       AD<double> I(i);
461                        Y[i] = U[I];
462                }
463       
464                // create f: X -> Y and stop tape recording
465                // Y[ X[0] ] = X[1] and other components of Y are zero.
466                CppAD::ADFun<double> F;
467                F.Dependent(X, Y); 
468
469                // Check number of VecAD vectors plus number of VecAD elements
470                ok &= (F.size_VecAD() == 2 + n + m); 
471       
472                CppAD::vector<double> x(n), y(m);
473                for(j = 0; j < n; j++)
474                        x[j] = double(j);
475
476                y = F.Forward(0, x);
477                for(i = 0; i < m; i++)
478                {       if( i != static_cast<size_t>(x[0]) )
479                                ok &= (y[i] == 0.);
480                        else    ok &= (y[i] == x[1]);
481                }
482
483                F.optimize();
484
485                // Check number of VecAD vectors plus number of VecAD elements
486                ok &= (F.size_VecAD() == 1 + m); 
487                y = F.Forward(0, x);
488                for(i = 0; i < m; i++)
489                {       if( i != static_cast<size_t>(x[0]) )
490                                ok &= (y[i] == 0.);
491                        else    ok &= (y[i] == x[1]);
492                }
493               
494                return ok;
495        }
496        bool depend_three(void)
497        {       // Power function is a special case for optimize
498                bool ok = true;
499                using CppAD::AD;
500                using CppAD::vector;
501
502                size_t n = 3;
503                size_t j;
504       
505                vector< AD<double> >    X(n), Y(n);
506                vector<double>          x(n), y(n); 
507       
508                for(j = 0; j < n; j++)
509                        X[j] = x[j] = double(j+2);
510       
511                CppAD::Independent(X);
512                       
513                Y[0] = pow(X[0], 2.0);
514                Y[1] = pow(2.0, X[1]);
515                Y[2] = pow(X[0], X[1]);
516       
517                CppAD::ADFun<double> F(X, Y);
518                F.optimize();
519                y = F.Forward(0, x);
520
521                // Use identically equal because the result of the operations
522                // have been stored as double and gaurd bits have been dropped.
523                // (This may not be true for some compiler in the future).
524                for(j = 0; j < n; j++)
525                        ok &= ( y[j] == Value(Y[j]) );
526
527                // check reverse mode derivative
528                vector<double>   w(n), dw(n); 
529                w[0] = 0.;
530                w[1] = 0.;
531                w[2] = 1.;
532                dw = F.Reverse(1, w);
533
534                double check = x[1] * pow( x[0], x[1] - 1. );
535                ok &= NearEqual( dw[0], check, eps, eps );
536
537                check = log( x[0] ) * pow( x[0], x[1] );
538                ok &= NearEqual( dw[1], check, eps, eps );
539
540                check = 0.;
541                ok &= NearEqual( dw[2], check, eps, eps );
542       
543                return ok;
544        }
545        bool depend_four(void)
546        {       // erf function is a special case for optimize
547                bool ok = true;
548# if CPPAD_COMPILER_HAS_ERF
549                using CppAD::AD;
550                using CppAD::vector;
551
552                size_t n = 1;
553                size_t m = 1;
554                vector< AD<double> > X(n), Y(m);
555                vector<double>       x(n);
556                X[0] = x[0] = double(0.5);
557
558                CppAD::Independent(X);
559
560                Y[0] = erf(X[0]) + erf(X[0]);
561
562                CppAD::ADFun<double> F(X, Y);
563
564                vector<double> y_original     = F.Forward(0, x);
565                size_t         size_original  = F.size_var();
566                F.optimize();
567                ok &= F.size_var() + 5 == size_original;
568                vector<double> y = F.Forward(0, x);
569                ok &=  NearEqual(y[0], y_original[0], eps, eps);
570# endif
571                return ok;
572        }
573        // ===================================================================
574        // Test duplicate operation analysis
575
576        template <class Vector>
577        void duplicate_fun
578        (const Vector& x, Vector& y, size_t& original, size_t& opt)
579        {       typedef typename Vector::value_type Scalar;
580                original = 0;
581                opt      = 0;
582
583                // unary operator where operand is arg[0] and one result
584                Scalar a1 = CppAD::exp(x[0]); 
585                original += 1;
586                opt      += 1;
587
588                // unary operator where operand is arg[0] and two results
589                Scalar b1 = CppAD::sin(x[1]);
590                original += 2;
591                opt      += 2;
592
593                // non-commutative binary operator where left is a variable
594                // and right is a parameter
595                Scalar c1 = x[2] - 3.;
596                original += 1;
597                opt      += 1;
598
599                // non-commutative binary operator where left is a parameter
600                // and right is a variable
601                Scalar d1 = 3. / x[3];
602                original += 1;
603                opt      += 1;
604
605                // non-commutative binary operator where left is a variable
606                // and right is a variable
607                Scalar e1 = pow(x[3], x[4]);
608                original += 3;
609                opt      += 3;
610
611                // commutative binary operator where  left is a variable
612                // and right is a parameter
613                Scalar f1 = x[5] * 5.;
614                original += 1;
615                opt      += 1;
616
617                // commutative binary operator where  left is a variable
618                // and right is a variable
619                Scalar g1 = x[5] + x[6];
620                original += 1;
621                opt      += 1;
622
623                // duplicate variables
624                Scalar a2 = CppAD::exp(x[0]);
625                Scalar b2 = CppAD::sin(x[1]);  // counts for 2 variables
626                Scalar c2 = x[2] - 3.;
627                Scalar d2 = 3. / x[3];
628                Scalar e2 = pow(x[3], x[4]);   // counts for 3 variables
629                Scalar f2 = 5. * x[5];
630                Scalar g2 = x[6] + x[5];
631                original += 10;
632
633                // result vector
634                y[0] = a1 * a2;
635                y[1] = b1 * b2;
636                y[2] = c1 * c2;
637                y[3] = d1 * d2;
638                y[4] = e1 * e2;
639                y[5] = f1 * f2;
640                y[6] = g1 * g2;
641                original += 7;
642                opt      += 7;
643
644                return;
645        }
646        bool duplicate_one(void)
647        {
648                bool ok = true;
649                using CppAD::AD;
650                size_t original;
651                size_t opt;
652                size_t i, j;
653       
654                // domain space vector
655                size_t n  = 7;
656                CppAD::vector< AD<double> > X(n);
657                for(j = 0; j < n; j++)
658                        X[j] = 1. / double(j + 1); 
659       
660                // declare independent variables and start tape recording
661                CppAD::Independent(X);
662       
663                // range space vector
664                size_t m = n;
665                CppAD::vector< AD<double> > Y(m);
666                duplicate_fun(X, Y, original, opt);
667       
668                // create f: X -> Y and stop tape recording
669                CppAD::ADFun<double> F;
670                F.Dependent(X, Y); 
671       
672                CppAD::vector<double> x(n), y(m), check(m);
673                for(j = 0; j < n; j++)
674                        x[j] = Value(X[j]);
675                y = F.Forward(0, x);
676                duplicate_fun(x, check, original, opt);
677                for(i = 0; i < m; i++)
678                        ok &= NearEqual(y[i], check[i], eps, eps);
679       
680                // Check size before optimization
681                ok &= F.size_var() == (n + 1 + original);
682       
683                // Optimize the operation sequence
684                F.optimize();
685       
686                // Check size after optimization
687                ok &= F.size_var() == (n + 1 + opt);
688       
689                // check result now
690                // (should have already been checked if NDEBUG not defined)
691                y = F.Forward(0, x);
692                for(i = 0; i < m; i++)
693                        ok &= NearEqual(y[i], check[i], eps, eps);
694       
695                return ok;
696        }
697        // -------------------------------------------------------------------
698        bool duplicate_two(void)
699        {       // test that duplicate expression removal is relative to
700                // new and not just old argument indices.
701                bool ok = true;
702                using CppAD::AD;
703                size_t i, j;
704
705                // domain space vector
706                size_t n  = 1;
707                CppAD::vector< AD<double> > X(n);
708                for(j = 0; j < n; j++)
709                        X[j] = double(j + 2); 
710
711                // range space vector
712                size_t m = 1;
713                CppAD::vector< AD<double> > Y(m);
714
715                // declare independent variables and start tape recording
716                CppAD::Independent(X);
717
718                // create a new variable
719                AD<double> A1 = X[0] - 2.;
720
721                // create a duplicate variable
722                AD<double> A2 = X[0] - 2.;
723
724                // create a new variable using first version of duplicate
725                AD<double> B1 = A1 / 2.;
726
727                // create a duplicate that can only be dectected using new
728                // argument indices
729                AD<double> B2 = A2 / 2.; 
730
731                // Make a new variable for result
732                // and make it depend on all the variables
733                Y[0] = B1 + B2;
734
735                // create f: X -> Y and stop tape recording
736                CppAD::ADFun<double> F;
737                F.Dependent(X, Y); 
738
739                // check number of variables in original function
740                ok &= (F.size_var() ==  1 + n + m + 4 ); 
741       
742                CppAD::vector<double> x(n), y(m);
743                for(j = 0; j < n; j++)
744                        x[j] = double(j + 2);
745
746                y   = F.Forward(0, x);
747                for(i = 0; i < m; i++)
748                        ok &= ( y[i] == Value( Y[i] ) );
749
750                F.optimize();
751
752                // check number of variables  in optimized version
753                ok &= (F.size_var() == 1 + n + m + 2 ); 
754
755                y   = F.Forward(0, x);
756                for(i = 0; i < m; i++)
757                        ok &= ( y[i] == Value( Y[i] ) );
758
759                return ok;
760        }
761        // -------------------------------------------------------------------
762        bool duplicate_three(void)
763        {       // test that duplicate expression removal is relative to
764                // new and not just old argument indices (commutative case).
765                bool ok = true;
766                using CppAD::AD;
767                size_t i, j;
768
769                // domain space vector
770                size_t n  = 1;
771                CppAD::vector< AD<double> > X(n);
772                for(j = 0; j < n; j++)
773                        X[j] = double(j + 2); 
774
775                // range space vector
776                size_t m = 1;
777                CppAD::vector< AD<double> > Y(m);
778
779                // declare independent variables and start tape recording
780                CppAD::Independent(X);
781
782                // create a new variable
783                AD<double> A1 = X[0] + 2.;
784
785                // create a duplicate variable
786                AD<double> A2 = 2. + X[0];
787
788                // create a new variable using first version of duplicate
789                AD<double> B1 = A1 * 2.;
790
791                // create a duplicate that can only be dectected using new
792                // argument indices
793                AD<double> B2 = 2. * A2; 
794
795                // Make a new variable for result
796                // and make it depend on all the variables
797                Y[0] = B1 + B2;
798
799                // create f: X -> Y and stop tape recording
800                CppAD::ADFun<double> F;
801                F.Dependent(X, Y); 
802
803                // check number of variables in original function
804                ok &= (F.size_var() ==  1 + n + m + 4 ); 
805       
806                CppAD::vector<double> x(n), y(m);
807                for(j = 0; j < n; j++)
808                        x[j] = double(j + 2);
809
810                y   = F.Forward(0, x);
811                for(i = 0; i < m; i++)
812                        ok &= ( y[i] == Value( Y[i] ) );
813
814                F.optimize();
815
816                // check number of variables  in optimized version
817                ok &= (F.size_var() == 1 + n + m + 2 ); 
818
819                y   = F.Forward(0, x);
820                for(i = 0; i < m; i++)
821                        ok &= ( y[i] == Value( Y[i] ) );
822
823                return ok;
824        }
825        // -------------------------------------------------------------------
826        bool duplicate_four(void)
827        {       // Check that unary expression matching not only checks hash code,
828                // and operator, but also operand (old bug that has been fixed).
829                bool ok = true;
830                using CppAD::AD;
831                size_t j;
832
833                // domain space vector
834                size_t n  = 1;
835                CppAD::vector< AD<double> > X(n);
836                X[0] = 1.;
837
838                // range space vector
839                size_t m = 1;
840                CppAD::vector< AD<double> > Y(m);
841
842                // declare independent variables and start tape recording
843                CppAD::Independent(X);
844
845                // check a huge number of same operation with different operands
846                size_t n_operations = std::min(
847                        size_t(CPPAD_HASH_TABLE_SIZE) + 5,
848                        size_t(std::numeric_limits<CPPAD_TAPE_ADDR_TYPE>::max()) - 5
849                );
850                Y[0] = X[0];
851                for(j = 0; j < n_operations; j++)
852                        Y[0] = abs(Y[0]);
853
854                // create f: X -> Y and stop tape recording
855                CppAD::ADFun<double> F;
856                F.Dependent(X, Y); 
857
858                // check number of variables in original function
859                ok &= (F.size_var() ==  1 + n + n_operations ); 
860       
861                CppAD::vector<double> x(n), y(m);
862                x[0] = 1.;
863
864                y   = F.Forward(0, x);
865                ok &= ( y[0] == Value( Y[0] ) );
866
867                F.optimize();
868
869                // check same number of variables in optimized version
870                ok &= (F.size_var() == 1 + n + n_operations ); 
871
872                y   = F.Forward(0, x);
873                ok &= ( y[0] == Value( Y[0] ) );
874
875                return ok;
876        }
877        // ====================================================================
878        bool cummulative_sum(void)
879        {       // test conversion of a sequence of additions and subtraction
880                // to a cummulative summation sequence.
881                bool ok = true;
882                using CppAD::AD;
883                size_t i, j;
884
885                // domain space vector
886                size_t n  = 7;
887                CppAD::vector< AD<double> > X(n);
888                for(j = 0; j < n; j++)
889                        X[j] = double(j + 2); 
890
891                size_t n_original = 1 + n;
892                size_t n_optimize = 1 + n;
893
894                // range space vector
895                size_t m = 2;
896                CppAD::vector< AD<double> > Y(m);
897
898                // declare independent variables and start tape recording
899                CppAD::Independent(X);
900
901                // Operations inside of optimize_cadd
902                Y[0] = 5. + (X[0] + X[1]) + (X[1] - X[2]) // Addvv, Subvv
903                     + (X[2] - 1.) + (2. - X[3])   // Subvp, Subpv
904                     + (X[4] + 3.) + (4. + X[5]);  // Addpv, Addpv (no Addvp)
905                n_original += 12;
906                n_optimize += 1;
907
908
909                // Operations inside of optimize_csub
910                Y[1] = 5. - (X[1] + X[2]) - (X[2] - X[3]) // Addvv, Subvv
911                     - (X[3] - 1.) - (2. - X[4])   // Subvp, Subpv
912                     - (X[5] + 3.) - (4. + X[6]);  // Addpv, Addpv (no Addvp)
913                n_original += 12;
914                n_optimize += 1;
915
916                CppAD::ADFun<double> F;
917                F.Dependent(X, Y); 
918
919                // check number of variables in original function
920                ok &= (F.size_var() ==  n_original ); 
921       
922                CppAD::vector<double> x(n), y(m);
923                for(j = 0; j < n; j++)
924                        x[j] = double(j + 2);
925
926                y   = F.Forward(0, x);
927                for(i = 0; i < m; i++)
928                        ok &= ( y[i] == Value( Y[i] ) );
929
930                F.optimize();
931
932                // check number of variables  in optimized version
933                ok &= (F.size_var() == n_optimize ); 
934
935                y   = F.Forward(0, x);
936                for(i = 0; i < m; i++)
937                        ok &= ( y[i] == Value( Y[i] ) );
938
939                return ok;
940        }
941        // -------------------------------------------------------------------
942        bool forward_csum(void)
943        {       bool ok = true;
944       
945                using namespace CppAD;
946       
947                // independent variable vector
948                CppAD::vector< AD<double> > X(2);
949                X[0] = 0.; 
950                X[1] = 1.;
951                Independent(X);
952       
953                // compute sum of elements in X
954                CppAD::vector< AD<double> > Y(1);
955                Y[0] = X[0] + X[0] + X[1];
956       
957                // create function object F : X -> Y
958                ADFun<double> F(X, Y);
959       
960                // now optimize the operation sequence
961                F.optimize();
962       
963                // use zero order to evaluate F[ (3, 4) ]
964                CppAD::vector<double>  x0( F.Domain() );
965                CppAD::vector<double>  y0( F.Range() );
966                x0[0]    = 3.;
967                x0[1]    = 4.;
968                y0       = F.Forward(0, x0);
969                ok      &= NearEqual(y0[0] , x0[0]+x0[0]+x0[1], eps, eps);
970       
971                // evaluate derivative of F in X[0] direction
972                CppAD::vector<double> x1( F.Domain() );
973                CppAD::vector<double> y1( F.Range() );
974                x1[0]    = 1.;
975                x1[1]    = 0.;
976                y1       = F.Forward(1, x1);
977                ok      &= NearEqual(y1[0] , x1[0]+x1[0]+x1[1], eps, eps);
978       
979                // evaluate second derivative of F in X[0] direction
980                CppAD::vector<double> x2( F.Domain() );
981                CppAD::vector<double> y2( F.Range() );
982                x2[0]       = 0.;
983                x2[1]       = 0.;
984                y2          = F.Forward(2, x2);
985                double F_00 = 2. * y2[0];
986                ok         &= NearEqual(F_00, 0., eps, eps);
987       
988                return ok;
989        }
990        // -------------------------------------------------------------------
991        bool reverse_csum(void)
992        {       bool ok = true;
993       
994                using namespace CppAD;
995       
996                // independent variable vector
997                CppAD::vector< AD<double> > X(2);
998                X[0] = 0.; 
999                X[1] = 1.;
1000                Independent(X);
1001       
1002                // compute sum of elements in X
1003                CppAD::vector< AD<double> > Y(1);
1004                Y[0] = X[0] - (X[0] - X[1]);
1005       
1006                // create function object F : X -> Y
1007                ADFun<double> F(X, Y);
1008       
1009                // now optimize the operation sequence
1010                F.optimize();
1011       
1012                // use zero order to evaluate F[ (3, 4) ]
1013                CppAD::vector<double>  x0( F.Domain() );
1014                CppAD::vector<double>  y0( F.Range() );
1015                x0[0]    = 3.;
1016                x0[1]    = 4.;
1017                y0       = F.Forward(0, x0);
1018                ok      &= NearEqual(y0[0] , x0[0]-x0[0]+x0[1], eps, eps);
1019       
1020                // evaluate derivative of F
1021                CppAD::vector<double> dF( F.Domain() );
1022                CppAD::vector<double> w( F.Range() );
1023                w[0]    = 1.;
1024                dF      = F.Reverse(1, w);
1025                ok     &= NearEqual(dF[0] , 0., eps, eps);
1026                ok     &= NearEqual(dF[1] , 1., eps, eps);
1027       
1028                return ok;
1029        }
1030        bool forward_sparse_jacobian()
1031        {       bool ok = true;
1032                using namespace CppAD;
1033       
1034                // dimension of the domain space
1035                size_t n = 3; 
1036       
1037                // dimension of the range space
1038                size_t m = 3;
1039       
1040                // independent variable vector
1041                CppAD::vector< AD<double> > X(n);
1042                X[0] = 2.; 
1043                X[1] = 3.;
1044                X[2] = 4.;
1045                Independent(X);
1046       
1047                // dependent variable vector
1048                CppAD::vector< AD<double> > Y(m);
1049       
1050                // check results vector
1051                CppAD::vector< bool >       Check(m * n);
1052       
1053                // initialize index into Y
1054                size_t index = 0;
1055       
1056                // Y[0]
1057                Y[index]             = X[0] + X[1] + 5.;
1058                Check[index * n + 0] = true;
1059                Check[index * n + 1] = true;
1060                Check[index * n + 2] = false;
1061                index++;
1062       
1063                // Y[1]
1064                Y[index]             = Y[0] - (X[1] + X[2]);
1065                Check[index * n + 0] = true;
1066                Check[index * n + 1] = true;
1067                Check[index * n + 2] = true;
1068                index++;
1069
1070                // Y[2]
1071                // 2DO: There is a subtitle issue that has to do with using reverse
1072                // jacobian sparsity patterns during the optimization process.
1073                // We need an option to include X[0] in the sparsity pattern
1074                // so the optimizer can know it affects the results.
1075                Y[index]             = CondExpLe(X[0], X[1], X[1]+X[1], X[2]-X[2]);
1076                Check[index * n + 0] = false;
1077                Check[index * n + 1] = true;
1078                Check[index * n + 2] = true;
1079                index++;
1080       
1081                // check final index
1082                assert( index == m );
1083       
1084                // create function object F : X -> Y
1085                ADFun<double> F(X, Y);
1086                F.optimize();
1087       
1088                // ---------------------------------------------------------
1089                // dependency matrix for the identity function
1090                CppAD::vector< std::set<size_t> > Sx(n);
1091                size_t i, j;
1092                for(i = 0; i < n; i++)
1093                {       assert( Sx[i].empty() );
1094                        Sx[i].insert(i);
1095                }
1096       
1097                // evaluate the dependency matrix for F(x)
1098                CppAD::vector< std::set<size_t> > Sy(m);
1099                Sy = F.ForSparseJac(n, Sx);
1100       
1101                // check values
1102                bool found;
1103                for(i = 0; i < m; i++)
1104                {       for(j = 0; j < n; j++)
1105                        {       found = Sy[i].find(j) != Sy[i].end();
1106                                ok &= (found == Check[i * n + j]);
1107                        }
1108                }       
1109       
1110                return ok;
1111        }
1112        bool reverse_sparse_jacobian()
1113        {       bool ok = true;
1114                using namespace CppAD;
1115       
1116                // dimension of the domain space
1117                size_t n = 3; 
1118       
1119                // dimension of the range space
1120                size_t m = 3;
1121       
1122                // independent variable vector
1123                CppAD::vector< AD<double> > X(n);
1124                X[0] = 2.; 
1125                X[1] = 3.;
1126                X[2] = 4.;
1127                Independent(X);
1128       
1129                // dependent variable vector
1130                CppAD::vector< AD<double> > Y(m);
1131       
1132                // check results vector
1133                CppAD::vector< bool >       Check(m * n);
1134       
1135                // initialize index into Y
1136                size_t index = 0;
1137       
1138                // Y[0]
1139                Y[index]             = X[0] + X[1] + 5.;
1140                Check[index * n + 0] = true;
1141                Check[index * n + 1] = true;
1142                Check[index * n + 2] = false;
1143                index++;
1144       
1145                // Y[1]
1146                Y[index]             = Y[0] - (X[1] + X[2]);
1147                Check[index * n + 0] = true;
1148                Check[index * n + 1] = true;
1149                Check[index * n + 2] = true;
1150                index++;
1151
1152                // Y[2]
1153                Y[index]             = CondExpLe(X[0], X[1], X[1]+X[1], X[2]-X[2]);
1154                Check[index * n + 0] = false;
1155                Check[index * n + 1] = true;
1156                Check[index * n + 2] = true;
1157                index++;
1158       
1159                // check final index
1160                assert( index == m );
1161       
1162                // create function object F : X -> Y
1163                ADFun<double> F(X, Y);
1164                F.optimize();
1165       
1166                // ----------------------------------------------------------
1167                // dependency matrix for the identity function
1168                CppAD::vector< bool > Py(m * m);
1169                size_t i, j;
1170                for(i = 0; i < m; i++)
1171                {       for(j = 0; j < m; j++)
1172                                Py[ i * m + j ] = (i == j);
1173                }
1174       
1175                // evaluate the dependency matrix for F(x)
1176                CppAD::vector< bool > Px(m * n);
1177                Px = F.RevSparseJac(m, Py);
1178       
1179                // check values
1180                for(i = 0; i < m; i++)
1181                {       for(j = 0; j < n; j++)
1182                                ok &= (Px[i * n + j] == Check[i * n + j]);
1183                }       
1184       
1185                return ok;
1186        }
1187        bool reverse_sparse_hessian(void)
1188        {       bool ok = true;
1189                using CppAD::AD;
1190                size_t i, j;
1191       
1192                size_t n = 3;
1193                CppAD::vector< AD<double> > X(n); 
1194                X[0] = 1.;
1195                X[1] = 2.;
1196                X[2] = 3.;
1197                CppAD::Independent(X);
1198       
1199                size_t m = 1;
1200                CppAD::vector< AD<double> > Y(m);
1201                Y[0] = CondExpGe( X[0], X[1], 
1202                        X[0] + (2. + X[1] + 3.) * X[1], 
1203                        X[0] + (2. + X[2] + 3.) * X[1]
1204                );
1205       
1206                CppAD::vector<bool> check(n * n);
1207                check[0 * n + 0] = false; // partial w.r.t. x[0], x[0]
1208                check[0 * n + 1] = false; //                x[0], x[1]
1209                check[0 * n + 2] = false; //                x[0], x[2]
1210
1211                check[1 * n + 0] = false; // partial w.r.t. x[1], x[0]
1212                check[1 * n + 1] = true;  //                x[1], x[1]
1213                check[1 * n + 2] = true;  //                x[1], x[2]
1214
1215                check[2 * n + 0] = false; // partial w.r.t. x[2], x[0]
1216                check[2 * n + 1] = true;  //                x[2], x[1]
1217                check[2 * n + 2] = false; //                x[2], x[2]
1218       
1219                // create function object F : X -> Y
1220                CppAD::ADFun<double> F(X, Y);
1221                F.optimize();
1222       
1223                // sparsity pattern for the identity function U(x) = x
1224                CppAD::vector<bool> Px(n * n);
1225                for(i = 0; i < n; i++)
1226                        for(j = 0; j < n; j++)
1227                                Px[ i * n + j ] = (i == j);
1228       
1229                // compute sparsity pattern for Jacobian of F(U(x))
1230                CppAD::vector<bool> P_jac(m * n);
1231                P_jac = F.ForSparseJac(n, Px);
1232       
1233                // compute sparsity pattern for Hessian of F_k ( U(x) )
1234                CppAD::vector<bool> Py(m);
1235                CppAD::vector<bool> Pxx(n * n);
1236                Py[0] = true;
1237                Pxx = F.RevSparseHes(n, Py);
1238                // check values
1239                for(i = 0; i < n * n; i++)
1240                        ok &= (Pxx[i] == check[i]);
1241       
1242                return ok;
1243        }
1244        // check that CondExp properly detects dependencies
1245        bool cond_exp_depend(void)
1246        {       bool ok = true;
1247                using CppAD::AD;
1248
1249                AD<double> zero(0.), one(1.), two(2.), three(3.);
1250       
1251                size_t n = 4;
1252                CppAD::vector< AD<double> > X(n); 
1253                X[0] = zero;
1254                X[1] = one;
1255                X[2] = two;
1256                X[3] = three;
1257                CppAD::Independent(X);
1258       
1259                size_t m = 4;
1260                CppAD::vector< AD<double> > Y(m);
1261                Y[0] = CondExpLt(X[0] + .5,  one,  two, three);
1262                Y[1] = CondExpLt(zero, X[1] + .5,  two, three);
1263                Y[2] = CondExpLt(zero,  one, X[2] + .5, three);
1264                Y[3] = CondExpLt(zero,  one,  two,  X[3] + .5);
1265
1266                CppAD::ADFun<double> f(X, Y);
1267                f.optimize();
1268
1269                CppAD::vector<double> x(n), y(m);
1270                size_t i;
1271                for(i = 0; i < n; i++)
1272                        x[i] = double(n - i);
1273                y    = f.Forward(0, x);
1274
1275                if( x[0] + .5 < 1. )
1276                        ok  &= y[0] == 2.;
1277                else    ok  &= y[0] == 3.;
1278                if( 0. < x[1] + .5 )
1279                        ok  &= y[1] == 2.;
1280                else    ok  &= y[1] == 3.;
1281                ok  &= y[2] == x[2] + .5;;
1282                ok  &= y[3] == 2.;
1283
1284                return ok;
1285        }
1286        // -------------------------------------------------------------------
1287        void my_union(
1288                std::set<size_t>&         result  ,
1289                const std::set<size_t>&   left    ,
1290                const std::set<size_t>&   right   )
1291        {       std::set<size_t> temp;
1292                std::set_union(
1293                        left.begin()              ,
1294                        left.end()                ,
1295                        right.begin()             ,
1296                        right.end()               ,
1297                        std::inserter(temp, temp.begin())
1298                );
1299                result.swap(temp);
1300        }
1301
1302        bool old_atomic_forward(
1303                size_t                         id ,
1304                size_t                          k , 
1305                size_t                          n ,
1306                size_t                          m ,
1307                const CppAD::vector<bool>&     vx ,
1308                CppAD::vector<bool>&           vy ,
1309                const CppAD::vector<double>&   tx , 
1310                CppAD::vector<double>&         ty )
1311        {       assert(n == 3 && m == 2);
1312                if( k > 0 ) 
1313                        return false;
1314
1315                // y[0] = x[0] + x[1]
1316                ty[0] = tx[0] + tx[1];
1317
1318                // y[1] = x[1] + x[2]
1319                ty[1] = tx[1] + tx[2];
1320               
1321                if( vy.size() > 0 )
1322                {       vy[0] = (vx[0] | vx[1]);
1323                        vy[1] = (vx[1] | vx[2]);
1324                }
1325                return true; 
1326        }
1327
1328        bool old_atomic_reverse(
1329                size_t                         id ,
1330                size_t                          k , 
1331                size_t                          n , 
1332                size_t                          m , 
1333                const CppAD::vector<double>&   tx , 
1334                const CppAD::vector<double>&   ty ,
1335                CppAD::vector<double>&         px ,
1336                const CppAD::vector<double>&   py )
1337        {       return false; }
1338
1339        bool old_atomic_for_jac_sparse(
1340                size_t                                  id ,
1341                size_t                                   n ,
1342                size_t                                   m ,
1343                size_t                                   q ,
1344                const CppAD::vector< std::set<size_t> >& r ,
1345                CppAD::vector< std::set<size_t>  >&      s )
1346        {       return false; }
1347
1348        bool old_atomic_rev_jac_sparse(
1349                size_t                                  id ,
1350                size_t                                   n ,
1351                size_t                                   m ,
1352                size_t                                   q ,
1353                CppAD::vector< std::set<size_t> >&       r ,
1354                const CppAD::vector< std::set<size_t> >& s )
1355        {       assert(n == 3 && m == 2);
1356                r[0].clear();
1357                r[1].clear();
1358                r[2].clear();
1359                // y[0] = x[0] + x[1]
1360                my_union(r[0], r[0], s[0]);
1361                my_union(r[1], r[1], s[0]);
1362                // y[1] = x[1] + x[2]
1363                my_union(r[1], r[1], s[1]);
1364                my_union(r[2], r[2], s[1]);
1365
1366                return true; 
1367        }
1368
1369        bool old_atomic_rev_hes_sparse(
1370                size_t                                  id ,
1371                size_t                                   n ,
1372                size_t                                   m ,
1373                size_t                                   q ,
1374                const CppAD::vector< std::set<size_t> >& r ,
1375                const CppAD::vector<bool>&               s ,
1376                CppAD::vector<bool>&                     t ,
1377                const CppAD::vector< std::set<size_t> >& u ,
1378                CppAD::vector< std::set<size_t> >&       v )
1379        {       return false; }
1380
1381        CPPAD_USER_ATOMIC(
1382                my_old_atomic             ,
1383                CppAD::vector              ,
1384                double                     ,
1385                old_atomic_forward        ,
1386                old_atomic_reverse        ,
1387                old_atomic_for_jac_sparse ,
1388                old_atomic_rev_jac_sparse ,
1389                old_atomic_rev_hes_sparse
1390        )
1391
1392        bool old_atomic_test(void)
1393        {       bool ok = true;
1394
1395                using CppAD::AD;
1396                size_t j;
1397                size_t n = 3;
1398                size_t m = 2;
1399                CppAD::vector< AD<double> > ax(n), ay(m), az(m);
1400                for(j = 0; j < n; j++)
1401                        ax[j] = AD<double>(j + 1);
1402                CppAD::Independent(ax);
1403
1404                size_t id = 0;
1405                // first call should stay in the tape
1406                my_old_atomic(id++, ax, ay);
1407                // second call will not get used
1408                my_old_atomic(id++, ax, az);
1409                // create function
1410                CppAD::ADFun<double> g(ax, ay);
1411                // should have 1 + n + m + m varaibles
1412                ok &= g.size_var() == (1 + n + m + m);
1413                g.optimize();
1414                // should have 1 + n + m varaibles
1415                ok &= g.size_var() == (1 + n + m);
1416
1417                // now test that the optimized function gives same results
1418                CppAD::vector<double> x(n), y(m);
1419                for(j = 0; j < n; j++)
1420                        x[j] = (j + 1) * (j + 1);
1421                y = g.Forward(0, x);
1422                // y[0] = x[0] + x[1]
1423                ok &= (y[0] == x[0] + x[1]);
1424                // y[1] = x[1] + x[2]
1425                ok &= (y[0] == x[0] + x[1]);
1426
1427                return ok;
1428        }
1429        bool not_identically_equal(void)
1430        {       bool ok = true;
1431                using CppAD::AD;
1432
1433                // independent variable vector
1434                size_t n = 5;
1435                CppAD::vector< AD<double> > ax(n);
1436                size_t j;
1437                for(j = 0; j < n; j++)
1438                        ax[j] = 1. / 3.;
1439                CppAD::Independent(ax);
1440       
1441                // dependent variable vector
1442                size_t m = 1;
1443                CppAD::vector< AD<double> > ay(m);
1444                ay[0]       = 0.;
1445                for(j = 0; j < n; j++)
1446                {       if( j % 2 == 0 )
1447                                ay[0] += ax[j];
1448                        else    ay[0] -= ax[j];
1449                }
1450                CppAD::ADFun<double> f(ax, ay);
1451       
1452                // Used to fail assert in optimize that forward mode results
1453                // are identically equal
1454                f.optimize();
1455       
1456                return ok;
1457        }
1458        // -----------------------------------------------------------------------
1459        double floor(const double& x)
1460        {       return std::floor(x); }
1461        CPPAD_DISCRETE_FUNCTION(double, floor)
1462        bool discrete_function(void)
1463        {       bool ok = true;
1464                using CppAD::vector;
1465       
1466                vector< CppAD::AD<double> > ax(1), ay(1);
1467                ax[0] = 0.0; 
1468                CppAD::Independent(ax);
1469                ay[0] =  floor(ax[0]) + floor(ax[0]); 
1470                CppAD::ADFun<double> f(ax, ay);
1471       
1472                size_t size_before = f.size_var();
1473                f.optimize(); 
1474                size_t size_after = f.size_var();
1475                ok &= size_after + 1 == size_before;
1476       
1477                vector<double> x(1), y(1);
1478                x[0] = -2.2;
1479                y    = f.Forward(0, x);
1480                ok &= y[0] == -6.0;
1481       
1482                return ok;
1483        }
1484        // ----------------------------------------------------------------
1485        void i_algo( 
1486                const CppAD::vector< CppAD::AD<double> >& ax ,
1487                      CppAD::vector< CppAD::AD<double> >& ay )
1488        {       ay[0] = 1.0 / ax[0]; }
1489        //
1490        // Test bug where atomic functions were not properly conditionally skipped.
1491        bool cond_exp_skip_atomic(void)
1492        {       bool ok = true;
1493                using CppAD::AD;
1494                using CppAD::vector;
1495       
1496                // Create a checkpoint version of the function i_algo
1497                vector< AD<double> > au(1), av(1), aw(1);
1498                au[0] = 1.0;
1499                CppAD::checkpoint<double> i_check("i_check", i_algo, au, av);
1500       
1501                // independent variable vector
1502                vector< AD<double> > ax(2), ay(1);
1503                ax[0] = 1.0;
1504                ax[1] = 2.0;
1505                Independent(ax);
1506       
1507                // call atomic function that does not get used
1508                au[0] = ax[0];
1509                i_check(au, av);
1510                au[0] = ax[1];
1511                i_check(au, aw);
1512                AD<double> zero = 0.0;
1513                ay[0] = CondExpGt(av[0], zero, av[0], aw[0]);
1514
1515                // create function object f : ax -> ay
1516                CppAD::ADFun<double> f(ax, ay);
1517
1518                // run case that skips the second call to afun
1519                // (can use trace in forward0sweep.hpp to see this).
1520                vector<double> x(2), y_before(1), y_after(1);
1521                x[0]      = 1.0;
1522                x[1]      = 2.0;
1523                y_before  = f.Forward(0, x);
1524                f.optimize();
1525                y_after   = f.Forward(0, x);
1526       
1527                ok &= y_before[0] == y_after[0];
1528               
1529                return ok;
1530        }
1531        //
1532        // Test bug where conditional dependence did not pass through
1533        // atomic functions
1534        bool cond_exp_atomic_dependence(void)
1535        {       bool ok = true;
1536                using CppAD::AD;
1537                using CppAD::vector;
1538       
1539                // Create a checkpoint version of the function i_algo
1540                vector< AD<double> > au(1), av(1), aw(1);
1541                au[0] = 1.0;
1542                CppAD::checkpoint<double> i_check("i_check", i_algo, au, av);
1543       
1544                vector< AD<double> > ax(2), ay(1);
1545                AD<double> zero = 0.0; 
1546                ax[0] = 1.0;
1547                ax[1] = 1.0;
1548                Independent(ax);
1549                av[0] = ax[0] + ax[1];
1550                i_check(av, aw);
1551                ay[0] = CondExpGt(aw[0], zero, zero, aw[0]);
1552                CppAD::ADFun<double> f;
1553                f.Dependent(ax, ay);
1554
1555                // run case that skips the second call to afun
1556                // (but not for order zero)
1557                vector<double> x(2), y_before(1), y_after(1);
1558                vector<double> dx(2), dy_before(1), dy_after(1);
1559                x[0]      = 1.0;
1560                x[1]      = 1.0;
1561                y_before  = f.Forward(0, x);
1562                dx[0]     = 2.0;
1563                dx[1]     = 2.0;
1564                dy_before = f.Forward(1, dx);
1565                f.optimize();
1566                y_after   = f.Forward(0, x);
1567                dy_after  = f.Forward(1, dx);
1568
1569                ok &= y_before[0]  == y_after[0];
1570                ok &= dy_before[0] == dy_after[0];
1571
1572                return ok;
1573        }
1574        // -----------------------------------------------------------------------
1575        // Test reverse mode conditionalay skipping commands.
1576        template <class Type>
1577        Type my_max(const CppAD::vector<Type>& arg)
1578        {       Type res = arg[0];
1579                for(size_t j = 0;j < arg.size(); j++)
1580                res = CondExpGt(res, arg[j], res, arg[j]);
1581                return res;
1582        }
1583        bool cond_exp_reverse(void)
1584        {       bool ok = true;
1585                size_t n = 3;
1586                using CppAD::vector;
1587                using CppAD::AD;
1588       
1589                vector< AD<double> > ax(n), ay(1);
1590                for(size_t j = 0; j < n; j++)
1591                        ax[j] = 1.0;
1592                Independent(ax);
1593                ay[0] = my_max(ax) + my_max(ax);
1594                CppAD::ADFun<double> f(ax, ay);
1595       
1596                f.optimize();
1597       
1598                vector<double> x(n), w(1), dx(n);
1599                for(size_t j = 0;j < n; j++)
1600                        x[j] = double(j);
1601                f.Forward(0, x);
1602                w[0] = 1.0;
1603                dx = f.Reverse(1, w);
1604                for(size_t j = 0; j < n; j++)
1605                {       if( j == n-1 )
1606                                ok &= dx[j] == 2.0;
1607                        else
1608                                ok &= dx[j] == 0.0;
1609                }
1610                return ok;
1611        }
1612        // Test case where an expression depends on both the true
1613        // and false cases (bug fixed 2014-12-22)
1614        bool cond_exp_both_true_and_false(void)
1615        {       bool ok = true;
1616                using CppAD::vector;
1617                using CppAD::AD;
1618
1619                // f(x) = x[0] + x[0] if x[0] >= 3
1620                //      = x[0] + x[1] otherwise
1621                vector< AD<double> > ax(2), ay(3);
1622                ax[0] = 1.0;
1623                ax[1] = 2.0;
1624                Independent(ax);
1625                AD<double> three(3);
1626                AD<double> value = ax[0] + ax[1];
1627                // a simple value
1628                ay[0]  = CppAD::CondExpGe(ax[0], three, value, value);
1629                // a  binary exprpression
1630                ay[1]  = CppAD::CondExpGe(ax[0], three, ax[0]-ax[1], ax[0]-ax[1]);
1631                // a unary expression
1632                ay[2]  = CppAD::CondExpGe(ax[0], three, exp(ax[0]), exp(ax[0]) );
1633                CppAD::ADFun<double> f(ax, ay);
1634                f.optimize();
1635
1636                // check case where x[0] >= 3
1637                vector<double> x(2), y(3);
1638                x[0] = 4.0;
1639                x[1] = 2.0;
1640                y    = f.Forward(0, x);
1641                ok  &= y[0] == x[0] + x[1];
1642                ok  &= y[1] == x[0] - x[1];
1643                ok  &= y[2] == exp(x[0]);
1644
1645                // check case where x[0] < 3
1646                x[0] = 1.0;
1647                x[1] = 2.0;
1648                y    = f.Forward(0, x);
1649                ok  &= y[0] == x[0] + x[1];
1650                ok  &= y[1] == x[0] - x[1];
1651                ok  &= y[2] == exp(x[0]);
1652
1653                return ok;
1654        }
1655}
1656
1657bool optimize(void)
1658{       bool ok = true;
1659
1660        atomic_sparsity_option = CppAD::atomic_base<double>::bool_sparsity_enum;
1661        for(size_t i = 0; i < 2; i++)
1662        {       // check conditional expression sparsity pattern
1663                // (used to optimize calls to atomic functions).
1664                ok     &= atomic_cond_exp_sparsity();
1665                // check optimizing out entire atomic function
1666                ok     &= atomic_cond_exp();
1667                // check optimizing out atomic arguments
1668                ok     &= atomic_no_used();
1669                ok     &= atomic_arguments();
1670                atomic_sparsity_option = 
1671                        CppAD::atomic_base<double>::set_sparsity_enum;
1672        }
1673        // check nested conditional expressions
1674        ok     &= nested_cond_exp();
1675        // check reverse dependency analysis optimization
1676        ok     &= depend_one();
1677        ok     &= depend_two();
1678        ok     &= depend_three();
1679        ok     &= depend_four();
1680        // check removal of duplicate expressions
1681        ok     &= duplicate_one();
1682        ok     &= duplicate_two();
1683        ok     &= duplicate_three();
1684        ok     &= duplicate_four();
1685        // convert sequence of additions to cummulative summation
1686        ok     &= cummulative_sum();
1687        ok     &= forward_csum();
1688        ok     &= reverse_csum();
1689        // sparsity patterns
1690        ok     &= forward_sparse_jacobian();
1691        ok     &= reverse_sparse_jacobian();
1692        ok     &= reverse_sparse_hessian();
1693        // check that CondExp properly detects dependencies
1694        ok     &= cond_exp_depend();
1695        // check old_atomic functions
1696        ok     &= old_atomic_test();
1697        // case where results are not identically equal
1698        ok     &= not_identically_equal();
1699        // case where a discrete function is used
1700        ok     &= discrete_function();
1701        // check conditional skip of an atomic function
1702        ok     &= cond_exp_skip_atomic();
1703        // check conditional dependence through atomic function
1704        ok     &= cond_exp_atomic_dependence();
1705        // check reverse mode conditional skipping
1706        ok     &= cond_exp_reverse();
1707        // check case where an expresion needed by both true and false case
1708        ok     &=  cond_exp_both_true_and_false();
1709        //
1710        CppAD::user_atomic<double>::clear();
1711        return ok;
1712}
Note: See TracBrowser for help on using the repository browser.