source: trunk/test_more/optimize.cpp @ 3732

Last change on this file since 3732 was 3732, checked in by bradbell, 4 years ago

merge to branch: trunk
from repository: https://github.com/coin-or/CppAD
start hash code: 401c99935c98da4b2086f3777f7d4bee3ebebb22
end hash code: 37afaa0d224e420a8cd052c18ceadf1c684ffd34

commit 37afaa0d224e420a8cd052c18ceadf1c684ffd34
Author: Brad Bell <bradbell@…>
Date: Thu Sep 24 08:05:47 2015 -0700

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