source: trunk/test_more/optimize.cpp @ 3642

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

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

commit ee39862aa18207f6085df3cfacf8353d7a0682c3
Author: Brad Bell <bradbell@…>
Date: Wed Feb 11 07:15:13 2015 -0700

optimize.hpp: fixed bug in optimization of conditional expressions.
optimize.cpp: test that demonstrates bug.

commit 2578d9745b6261d802337bd14b38efedacf680b0
Author: Brad Bell <bradbell@…>
Date: Wed Feb 11 04:45:33 2015 -0700

cond_exp.sh: a new bug report.

commit 0ab482aef633270ebfc3db2e2533d5ad8d696aaa
Author: Brad Bell <bradbell@…>
Date: Wed Feb 11 04:32:05 2015 -0700

erf.sh: This bug was fixed; see whats new 2015-02-10.

  • Property svn:keywords set to Id
File size: 43.9 KB
Line 
1/* $Id: optimize.cpp 3642 2015-02-11 14:37:10Z 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        // check that CondExp properly handels expressions that get
1287        // removed during opitmization
1288        bool cond_exp_removed(void)
1289        {       bool ok = true;
1290                using CppAD::AD;
1291                AD<double> zero(0.);
1292       
1293                size_t n = 1;
1294                CppAD::vector< AD<double> > X(n); 
1295                X[0] = 1.0;
1296                CppAD::Independent(X);
1297       
1298                size_t m = 1;
1299                CppAD::vector< AD<double> > Y(m);
1300
1301                AD<double> true_case  = sin(X[0]) + sin(X[0]);
1302                AD<double> false_case = cos(X[0]) + cos(X[0]);
1303                Y[0] = CondExpLt(X[0],  zero,  true_case, false_case);
1304
1305                CppAD::ADFun<double> f(X, Y);
1306                f.optimize();
1307
1308                CppAD::vector<double> x(n), y(m), w(m), dw(n);
1309                x[0] = 1.0;
1310                y    = f.Forward(0, x);
1311                ok &= NearEqual(y[0], false_case, eps, eps);
1312
1313                w[0] = 1.0;
1314                dw   = f.Reverse(1, w);
1315                // derivative of cos is minus sin
1316                ok &= NearEqual(dw[0], - true_case, eps, eps);
1317
1318                return ok;
1319        }
1320        // -------------------------------------------------------------------
1321        void my_union(
1322                std::set<size_t>&         result  ,
1323                const std::set<size_t>&   left    ,
1324                const std::set<size_t>&   right   )
1325        {       std::set<size_t> temp;
1326                std::set_union(
1327                        left.begin()              ,
1328                        left.end()                ,
1329                        right.begin()             ,
1330                        right.end()               ,
1331                        std::inserter(temp, temp.begin())
1332                );
1333                result.swap(temp);
1334        }
1335
1336        bool old_atomic_forward(
1337                size_t                         id ,
1338                size_t                          k , 
1339                size_t                          n ,
1340                size_t                          m ,
1341                const CppAD::vector<bool>&     vx ,
1342                CppAD::vector<bool>&           vy ,
1343                const CppAD::vector<double>&   tx , 
1344                CppAD::vector<double>&         ty )
1345        {       assert(n == 3 && m == 2);
1346                if( k > 0 ) 
1347                        return false;
1348
1349                // y[0] = x[0] + x[1]
1350                ty[0] = tx[0] + tx[1];
1351
1352                // y[1] = x[1] + x[2]
1353                ty[1] = tx[1] + tx[2];
1354               
1355                if( vy.size() > 0 )
1356                {       vy[0] = (vx[0] | vx[1]);
1357                        vy[1] = (vx[1] | vx[2]);
1358                }
1359                return true; 
1360        }
1361
1362        bool old_atomic_reverse(
1363                size_t                         id ,
1364                size_t                          k , 
1365                size_t                          n , 
1366                size_t                          m , 
1367                const CppAD::vector<double>&   tx , 
1368                const CppAD::vector<double>&   ty ,
1369                CppAD::vector<double>&         px ,
1370                const CppAD::vector<double>&   py )
1371        {       return false; }
1372
1373        bool old_atomic_for_jac_sparse(
1374                size_t                                  id ,
1375                size_t                                   n ,
1376                size_t                                   m ,
1377                size_t                                   q ,
1378                const CppAD::vector< std::set<size_t> >& r ,
1379                CppAD::vector< std::set<size_t>  >&      s )
1380        {       return false; }
1381
1382        bool old_atomic_rev_jac_sparse(
1383                size_t                                  id ,
1384                size_t                                   n ,
1385                size_t                                   m ,
1386                size_t                                   q ,
1387                CppAD::vector< std::set<size_t> >&       r ,
1388                const CppAD::vector< std::set<size_t> >& s )
1389        {       assert(n == 3 && m == 2);
1390                r[0].clear();
1391                r[1].clear();
1392                r[2].clear();
1393                // y[0] = x[0] + x[1]
1394                my_union(r[0], r[0], s[0]);
1395                my_union(r[1], r[1], s[0]);
1396                // y[1] = x[1] + x[2]
1397                my_union(r[1], r[1], s[1]);
1398                my_union(r[2], r[2], s[1]);
1399
1400                return true; 
1401        }
1402
1403        bool old_atomic_rev_hes_sparse(
1404                size_t                                  id ,
1405                size_t                                   n ,
1406                size_t                                   m ,
1407                size_t                                   q ,
1408                const CppAD::vector< std::set<size_t> >& r ,
1409                const CppAD::vector<bool>&               s ,
1410                CppAD::vector<bool>&                     t ,
1411                const CppAD::vector< std::set<size_t> >& u ,
1412                CppAD::vector< std::set<size_t> >&       v )
1413        {       return false; }
1414
1415        CPPAD_USER_ATOMIC(
1416                my_old_atomic             ,
1417                CppAD::vector              ,
1418                double                     ,
1419                old_atomic_forward        ,
1420                old_atomic_reverse        ,
1421                old_atomic_for_jac_sparse ,
1422                old_atomic_rev_jac_sparse ,
1423                old_atomic_rev_hes_sparse
1424        )
1425
1426        bool old_atomic_test(void)
1427        {       bool ok = true;
1428
1429                using CppAD::AD;
1430                size_t j;
1431                size_t n = 3;
1432                size_t m = 2;
1433                CppAD::vector< AD<double> > ax(n), ay(m), az(m);
1434                for(j = 0; j < n; j++)
1435                        ax[j] = AD<double>(j + 1);
1436                CppAD::Independent(ax);
1437
1438                size_t id = 0;
1439                // first call should stay in the tape
1440                my_old_atomic(id++, ax, ay);
1441                // second call will not get used
1442                my_old_atomic(id++, ax, az);
1443                // create function
1444                CppAD::ADFun<double> g(ax, ay);
1445                // should have 1 + n + m + m varaibles
1446                ok &= g.size_var() == (1 + n + m + m);
1447                g.optimize();
1448                // should have 1 + n + m varaibles
1449                ok &= g.size_var() == (1 + n + m);
1450
1451                // now test that the optimized function gives same results
1452                CppAD::vector<double> x(n), y(m);
1453                for(j = 0; j < n; j++)
1454                        x[j] = (j + 1) * (j + 1);
1455                y = g.Forward(0, x);
1456                // y[0] = x[0] + x[1]
1457                ok &= (y[0] == x[0] + x[1]);
1458                // y[1] = x[1] + x[2]
1459                ok &= (y[0] == x[0] + x[1]);
1460
1461                return ok;
1462        }
1463        bool not_identically_equal(void)
1464        {       bool ok = true;
1465                using CppAD::AD;
1466
1467                // independent variable vector
1468                size_t n = 5;
1469                CppAD::vector< AD<double> > ax(n);
1470                size_t j;
1471                for(j = 0; j < n; j++)
1472                        ax[j] = 1. / 3.;
1473                CppAD::Independent(ax);
1474       
1475                // dependent variable vector
1476                size_t m = 1;
1477                CppAD::vector< AD<double> > ay(m);
1478                ay[0]       = 0.;
1479                for(j = 0; j < n; j++)
1480                {       if( j % 2 == 0 )
1481                                ay[0] += ax[j];
1482                        else    ay[0] -= ax[j];
1483                }
1484                CppAD::ADFun<double> f(ax, ay);
1485       
1486                // Used to fail assert in optimize that forward mode results
1487                // are identically equal
1488                f.optimize();
1489       
1490                return ok;
1491        }
1492        // -----------------------------------------------------------------------
1493        double floor(const double& x)
1494        {       return std::floor(x); }
1495        CPPAD_DISCRETE_FUNCTION(double, floor)
1496        bool discrete_function(void)
1497        {       bool ok = true;
1498                using CppAD::vector;
1499       
1500                vector< CppAD::AD<double> > ax(1), ay(1);
1501                ax[0] = 0.0; 
1502                CppAD::Independent(ax);
1503                ay[0] =  floor(ax[0]) + floor(ax[0]); 
1504                CppAD::ADFun<double> f(ax, ay);
1505       
1506                size_t size_before = f.size_var();
1507                f.optimize(); 
1508                size_t size_after = f.size_var();
1509                ok &= size_after + 1 == size_before;
1510       
1511                vector<double> x(1), y(1);
1512                x[0] = -2.2;
1513                y    = f.Forward(0, x);
1514                ok &= y[0] == -6.0;
1515       
1516                return ok;
1517        }
1518        // ----------------------------------------------------------------
1519        void i_algo( 
1520                const CppAD::vector< CppAD::AD<double> >& ax ,
1521                      CppAD::vector< CppAD::AD<double> >& ay )
1522        {       ay[0] = 1.0 / ax[0]; }
1523        //
1524        // Test bug where atomic functions were not properly conditionally skipped.
1525        bool cond_exp_skip_atomic(void)
1526        {       bool ok = true;
1527                using CppAD::AD;
1528                using CppAD::vector;
1529       
1530                // Create a checkpoint version of the function i_algo
1531                vector< AD<double> > au(1), av(1), aw(1);
1532                au[0] = 1.0;
1533                CppAD::checkpoint<double> i_check("i_check", i_algo, au, av);
1534       
1535                // independent variable vector
1536                vector< AD<double> > ax(2), ay(1);
1537                ax[0] = 1.0;
1538                ax[1] = 2.0;
1539                Independent(ax);
1540       
1541                // call atomic function that does not get used
1542                au[0] = ax[0];
1543                i_check(au, av);
1544                au[0] = ax[1];
1545                i_check(au, aw);
1546                AD<double> zero = 0.0;
1547                ay[0] = CondExpGt(av[0], zero, av[0], aw[0]);
1548
1549                // create function object f : ax -> ay
1550                CppAD::ADFun<double> f(ax, ay);
1551
1552                // run case that skips the second call to afun
1553                // (can use trace in forward0sweep.hpp to see this).
1554                vector<double> x(2), y_before(1), y_after(1);
1555                x[0]      = 1.0;
1556                x[1]      = 2.0;
1557                y_before  = f.Forward(0, x);
1558                f.optimize();
1559                y_after   = f.Forward(0, x);
1560       
1561                ok &= y_before[0] == y_after[0];
1562               
1563                return ok;
1564        }
1565        //
1566        // Test bug where conditional dependence did not pass through
1567        // atomic functions
1568        bool cond_exp_atomic_dependence(void)
1569        {       bool ok = true;
1570                using CppAD::AD;
1571                using CppAD::vector;
1572       
1573                // Create a checkpoint version of the function i_algo
1574                vector< AD<double> > au(1), av(1), aw(1);
1575                au[0] = 1.0;
1576                CppAD::checkpoint<double> i_check("i_check", i_algo, au, av);
1577       
1578                vector< AD<double> > ax(2), ay(1);
1579                AD<double> zero = 0.0; 
1580                ax[0] = 1.0;
1581                ax[1] = 1.0;
1582                Independent(ax);
1583                av[0] = ax[0] + ax[1];
1584                i_check(av, aw);
1585                ay[0] = CondExpGt(aw[0], zero, zero, aw[0]);
1586                CppAD::ADFun<double> f;
1587                f.Dependent(ax, ay);
1588
1589                // run case that skips the second call to afun
1590                // (but not for order zero)
1591                vector<double> x(2), y_before(1), y_after(1);
1592                vector<double> dx(2), dy_before(1), dy_after(1);
1593                x[0]      = 1.0;
1594                x[1]      = 1.0;
1595                y_before  = f.Forward(0, x);
1596                dx[0]     = 2.0;
1597                dx[1]     = 2.0;
1598                dy_before = f.Forward(1, dx);
1599                f.optimize();
1600                y_after   = f.Forward(0, x);
1601                dy_after  = f.Forward(1, dx);
1602
1603                ok &= y_before[0]  == y_after[0];
1604                ok &= dy_before[0] == dy_after[0];
1605
1606                return ok;
1607        }
1608        // -----------------------------------------------------------------------
1609        // Test reverse mode conditionalay skipping commands.
1610        template <class Type>
1611        Type my_max(const CppAD::vector<Type>& arg)
1612        {       Type res = arg[0];
1613                for(size_t j = 0;j < arg.size(); j++)
1614                res = CondExpGt(res, arg[j], res, arg[j]);
1615                return res;
1616        }
1617        bool cond_exp_reverse(void)
1618        {       bool ok = true;
1619                size_t n = 3;
1620                using CppAD::vector;
1621                using CppAD::AD;
1622       
1623                vector< AD<double> > ax(n), ay(1);
1624                for(size_t j = 0; j < n; j++)
1625                        ax[j] = 1.0;
1626                Independent(ax);
1627                ay[0] = my_max(ax) + my_max(ax);
1628                CppAD::ADFun<double> f(ax, ay);
1629       
1630                f.optimize();
1631       
1632                vector<double> x(n), w(1), dx(n);
1633                for(size_t j = 0;j < n; j++)
1634                        x[j] = double(j);
1635                f.Forward(0, x);
1636                w[0] = 1.0;
1637                dx = f.Reverse(1, w);
1638                for(size_t j = 0; j < n; j++)
1639                {       if( j == n-1 )
1640                                ok &= dx[j] == 2.0;
1641                        else
1642                                ok &= dx[j] == 0.0;
1643                }
1644                return ok;
1645        }
1646        // Test case where an expression depends on both the true
1647        // and false cases (bug fixed 2014-12-22)
1648        bool cond_exp_both_true_and_false(void)
1649        {       bool ok = true;
1650                using CppAD::vector;
1651                using CppAD::AD;
1652
1653                // f(x) = x[0] + x[0] if x[0] >= 3
1654                //      = x[0] + x[1] otherwise
1655                vector< AD<double> > ax(2), ay(3);
1656                ax[0] = 1.0;
1657                ax[1] = 2.0;
1658                Independent(ax);
1659                AD<double> three(3);
1660                AD<double> value = ax[0] + ax[1];
1661                // a simple value
1662                ay[0]  = CppAD::CondExpGe(ax[0], three, value, value);
1663                // a  binary exprpression
1664                ay[1]  = CppAD::CondExpGe(ax[0], three, ax[0]-ax[1], ax[0]-ax[1]);
1665                // a unary expression
1666                ay[2]  = CppAD::CondExpGe(ax[0], three, exp(ax[0]), exp(ax[0]) );
1667                CppAD::ADFun<double> f(ax, ay);
1668                f.optimize();
1669
1670                // check case where x[0] >= 3
1671                vector<double> x(2), y(3);
1672                x[0] = 4.0;
1673                x[1] = 2.0;
1674                y    = f.Forward(0, x);
1675                ok  &= y[0] == x[0] + x[1];
1676                ok  &= y[1] == x[0] - x[1];
1677                ok  &= y[2] == exp(x[0]);
1678
1679                // check case where x[0] < 3
1680                x[0] = 1.0;
1681                x[1] = 2.0;
1682                y    = f.Forward(0, x);
1683                ok  &= y[0] == x[0] + x[1];
1684                ok  &= y[1] == x[0] - x[1];
1685                ok  &= y[2] == exp(x[0]);
1686
1687                return ok;
1688        }
1689}
1690
1691bool optimize(void)
1692{       bool ok = true;
1693        atomic_sparsity_option = CppAD::atomic_base<double>::bool_sparsity_enum;
1694        for(size_t i = 0; i < 2; i++)
1695        {       // check conditional expression sparsity pattern
1696                // (used to optimize calls to atomic functions).
1697                ok     &= atomic_cond_exp_sparsity();
1698                // check optimizing out entire atomic function
1699                ok     &= atomic_cond_exp();
1700                // check optimizing out atomic arguments
1701                ok     &= atomic_no_used();
1702                ok     &= atomic_arguments();
1703                atomic_sparsity_option = 
1704                        CppAD::atomic_base<double>::set_sparsity_enum;
1705        }
1706        // check nested conditional expressions
1707        ok     &= nested_cond_exp();
1708        // check reverse dependency analysis optimization
1709        ok     &= depend_one();
1710        ok     &= depend_two();
1711        ok     &= depend_three();
1712        ok     &= depend_four();
1713        // check removal of duplicate expressions
1714        ok     &= duplicate_one();
1715        ok     &= duplicate_two();
1716        ok     &= duplicate_three();
1717        ok     &= duplicate_four();
1718        // convert sequence of additions to cummulative summation
1719        ok     &= cummulative_sum();
1720        ok     &= forward_csum();
1721        ok     &= reverse_csum();
1722        // sparsity patterns
1723        ok     &= forward_sparse_jacobian();
1724        ok     &= reverse_sparse_jacobian();
1725        ok     &= reverse_sparse_hessian();
1726        // check that CondExp properly detects dependencies
1727        ok     &= cond_exp_depend();
1728        // check that it properly handles expressions that have been removed
1729        ok     &= cond_exp_removed();
1730        // check old_atomic functions
1731        ok     &= old_atomic_test();
1732        // case where results are not identically equal
1733        ok     &= not_identically_equal();
1734        // case where a discrete function is used
1735        ok     &= discrete_function();
1736        // check conditional skip of an atomic function
1737        ok     &= cond_exp_skip_atomic();
1738        // check conditional dependence through atomic function
1739        ok     &= cond_exp_atomic_dependence();
1740        // check reverse mode conditional skipping
1741        ok     &= cond_exp_reverse();
1742        // check case where an expresion needed by both true and false case
1743        ok     &=  cond_exp_both_true_and_false();
1744        //
1745        CppAD::user_atomic<double>::clear();
1746        return ok;
1747}
Note: See TracBrowser for help on using the repository browser.