source: branches/opt_cond_exp/test_more/optimize.cpp @ 2977

Last change on this file since 2977 was 2977, checked in by bradbell, 7 years ago

Backout simplification that removed atomic sparsity from optimization.

optimize.cpp: test case that demonstrates the need.

  • Property svn:keywords set to Id
File size: 29.4 KB
Line 
1/* $Id: optimize.cpp 2977 2013-10-20 12:48:12Z bradbell $ */
2/* --------------------------------------------------------------------------
3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 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        // -------------------------------------------------------------------
20        // Test of optimizing out arguments to an atomic function
21        void algo( 
22                const CppAD::vector< CppAD::AD<double> >& ax ,
23                      CppAD::vector< CppAD::AD<double> >& ay )
24        {       ay = ax; }
25
26        bool atomic_arguments(void)
27        {       bool ok = true;
28                using CppAD::AD;
29                using CppAD::vector;
30                vector< AD<double> > au(2), aw(2), ax(2), ay(1);
31
32                // create atomic function corresponding to algo
33                au[0] = 1.0;
34                au[1] = 2.0;
35                CppAD::checkpoint<double> algo_check("algo", algo, au, ax);
36
37                // start recording a new function
38                CppAD::Independent(ax);
39
40                // now use algo_check during the recording
41                au[0] = ax[0] + ax[1]; // this argument requires a new variable
42                au[1] = ax[0] - ax[1]; // this argument also requires a new variable
43                algo_check(au, aw);
44
45                // now create f(x) = x_0 - x_1
46                ay[0] = aw[0];
47                CppAD::ADFun<double> f(ax, ay);
48
49                // number of variables before optimization
50                size_t n_before = f.size_var();
51 
52                // now optimize f so that the calculation of au[1] is removed
53                f.optimize();
54
55                // check difference in number of variables
56                size_t n_after = f.size_var();
57                ok &= n_before == n_after + 1;
58
59                // now compute and check a forward mode calculation
60                vector<double> x(2), y(1);
61                x[0] = 5.0;
62                x[1] = 6.0;
63                y    = f.Forward(0, x);
64                ok  &= (y[0] == x[0] + x[1]); 
65
66                return ok;
67        }
68
69        // -------------------------------------------------------------------
70        // Test the reverse dependency analysis optimization
71        template <class Vector>
72        void depend_fun
73        (const Vector& x, Vector& y, size_t& original, size_t& opt)
74        {       typedef typename Vector::value_type Scalar;
75                Scalar a;
76                Scalar one(1), two(2), three(3), four(4);
77                original = 0;
78                opt      = 0;
79
80                // unary operator where operand is arg[0]
81                a = CppAD::abs(x[0]); 
82                if( a < 1. )
83                        y[0] = sin(x[0]);
84                else    y[0] = cos(x[0]); 
85                original += 3;
86                opt      += 2;
87
88                // binary operator where left operand is a variable
89                // and right operand is a parameter
90                a = x[1] + 2.;
91                if( a < 3. )
92                        y[1] = x[1] * 3.;
93                else    y[1] = x[1] / 2.;
94                original += 2;
95                opt      += 1;
96
97                // binary operator where left operand is a parameter
98                // and right operation is a variable
99                a = 2. - x[2];
100                if( a < 4. )
101                        y[2] = 3. / x[2];
102                else    y[2] = 4. + x[2];
103                original += 2;
104                opt      += 1;
105
106                // binary operator where both operands are variables
107                a = x[3] - x[2];
108                if( a < 4. )
109                        y[3] = x[3] / x[2];
110                else    y[3] = x[3] + x[2];
111                original += 2;
112                opt      += 1;
113
114                // this conditional expression that will be optimized out
115                a = CppAD::CondExpLt(x[0], x[1], x[2], x[3]);
116                // 1 of the following 2 conditional expressions will be kept
117                if( a < 5. )
118                        y[4] = CppAD::CondExpLt(x[4], one, two, three);
119                else    y[4] = CppAD::CondExpLt(x[4], two, three, four);
120                original += 2;
121                opt      += 1;
122
123                // Make sure that a parameter dependent variable
124                // is not optimized out of the operation sequence.
125                // In addition, we do not use the argument x[5], to
126                // make sure it is not optimized out.
127                y[5] = 1.;
128                original += 1;
129                opt      += 1;
130
131                return;
132        }
133
134        bool depend_one(void)
135        {       // Test all except for VecAD operations
136                bool ok = true;
137                using CppAD::AD;
138                size_t original;
139                size_t opt;
140                size_t i, j;
141       
142                // domain space vector
143                size_t n  = 6;
144                CppAD::vector< AD<double> > X(n);
145                for(j = 0; j < n; j++)
146                        X[j] = 1. / double(j + 1); 
147       
148                // declare independent variables and start tape recording
149                CppAD::Independent(X);
150       
151                // range space vector
152                size_t m = n;
153                CppAD::vector< AD<double> > Y(m);
154                depend_fun(X, Y, original, opt);
155       
156                // create f: X -> Y and stop tape recording
157                CppAD::ADFun<double> F;
158                F.Dependent(X, Y); 
159       
160                CppAD::vector<double> x(n), y(m), check(m);
161                for(j = 0; j < n; j++)
162                        x[j] = Value(X[j]);
163                y = F.Forward(0, x);
164                depend_fun(x, check, original, opt);
165                for(i = 0; i < m; i++)
166                        ok &= (y[i] == check[i]);
167       
168                // Check size before optimization
169                ok &= F.size_var() == (n + 1 + original);
170       
171                // Optimize the operation sequence
172                F.optimize();
173       
174                // Check size after optimization
175                ok &= F.size_var() == (n + 1 + opt);
176       
177                // check result now
178                // (should have already been checked if NDEBUG not defined)
179                y = F.Forward(0, x);
180                for(i = 0; i < m; i++)
181                        ok &= (y[i] == check[i]);
182       
183                return ok;
184        }
185
186        bool depend_two(void)
187        {       // Test VecAD operations
188                bool ok = true;
189                using CppAD::AD;
190                size_t i, j;
191
192                // domain space vector
193                size_t n  = 2;
194                CppAD::vector< AD<double> > X(n);
195                for(j = 0; j < n; j++)
196                        X[j] = double(j); 
197
198                // range space vector
199                size_t m = 3;
200                CppAD::vector< AD<double> > Y(m);
201
202                CppAD::VecAD<double> U(m);
203                CppAD::VecAD<double> V(n);
204                for(i = 0; i < m; i++)
205                        U[i] = 0;
206                for(j = 0; j < n; j++)
207                        V[j] = 0;
208       
209                // declare independent variables and start tape recording
210                CppAD::Independent(X);
211
212                // first vecad vector that is a variable
213                U[ X[0] ] = X[1];
214
215                // second vecad vector that is a variable
216                V[ X[0] ] = X[1];
217
218                // Make dependency for vecad vectors different that for
219                // variables because original code used worng dependency info.
220                // Y does not depend on the first variable in the tape; i.e.
221                // the one corresponding to the BeginOp. So make it depend
222                // on the first vecad vector in the tape.
223                for(i = 0; i < m; i++)
224                {       AD<double> I(i);
225                        Y[i] = U[I];
226                }
227       
228                // create f: X -> Y and stop tape recording
229                // Y[ X[0] ] = X[1] and other components of Y are zero.
230                CppAD::ADFun<double> F;
231                F.Dependent(X, Y); 
232
233                // Check number of VecAD vectors plus number of VecAD elements
234                ok &= (F.size_VecAD() == 2 + n + m); 
235       
236                CppAD::vector<double> x(n), y(m);
237                for(j = 0; j < n; j++)
238                        x[j] = double(j);
239
240                y = F.Forward(0, x);
241                for(i = 0; i < m; i++)
242                {       if( i != static_cast<size_t>(x[0]) )
243                                ok &= (y[i] == 0.);
244                        else    ok &= (y[i] == x[1]);
245                }
246
247                F.optimize();
248
249                // Check number of VecAD vectors plus number of VecAD elements
250                ok &= (F.size_VecAD() == 1 + m); 
251                y = F.Forward(0, x);
252                for(i = 0; i < m; i++)
253                {       if( i != static_cast<size_t>(x[0]) )
254                                ok &= (y[i] == 0.);
255                        else    ok &= (y[i] == x[1]);
256                }
257               
258                return ok;
259        }
260        bool depend_three(void)
261        {       // Power function is a special case for optimize
262                bool ok = true;
263                using CppAD::AD;
264                using CppAD::vector;
265
266                size_t n = 3;
267                size_t j;
268       
269                vector< AD<double> >    X(n), Y(n);
270                vector<double>          x(n), y(n); 
271       
272                for(j = 0; j < n; j++)
273                        X[j] = x[j] = double(j+2);
274       
275                CppAD::Independent(X);
276                       
277                Y[0] = pow(X[0], 2.0);
278                Y[1] = pow(2.0, X[1]);
279                Y[2] = pow(X[0], X[1]);
280       
281                CppAD::ADFun<double> F(X, Y);
282                F.optimize();
283                y = F.Forward(0, x);
284
285                // Use identically equal because the result of the operations
286                // have been stored as double and gaurd bits have been dropped.
287                // (This may not be true for some compiler in the future).
288                for(j = 0; j < n; j++)
289                        ok &= ( y[j] == Value(Y[j]) );
290
291                // check reverse mode derivative
292                vector<double>   w(n), dw(n); 
293                w[0] = 0.;
294                w[1] = 0.;
295                w[2] = 1.;
296                dw = F.Reverse(1, w);
297
298                double eps = 20. * std::numeric_limits<double>::epsilon();
299                double check = x[1] * pow( x[0], x[1] - 1. );
300                ok &= CppAD::NearEqual( dw[0], check, eps, eps );
301
302                check = log( x[0] ) * pow( x[0], x[1] );
303                ok &= CppAD::NearEqual( dw[1], check, eps, eps );
304
305                check = 0.;
306                ok &= CppAD::NearEqual( dw[2], check, eps, eps );
307       
308                return ok;
309        }
310        // ===================================================================
311        // Test duplicate operation analysis
312
313        template <class Vector>
314        void duplicate_fun
315        (const Vector& x, Vector& y, size_t& original, size_t& opt)
316        {       typedef typename Vector::value_type Scalar;
317                original = 0;
318                opt      = 0;
319
320                // unary operator where operand is arg[0] and one result
321                Scalar a1 = CppAD::exp(x[0]); 
322                original += 1;
323                opt      += 1;
324
325                // unary operator where operand is arg[0] and two results
326                Scalar b1 = CppAD::sin(x[1]);
327                original += 2;
328                opt      += 2;
329
330                // non-commutative binary operator where left is a variable
331                // and right is a parameter
332                Scalar c1 = x[2] - 3.;
333                original += 1;
334                opt      += 1;
335
336                // non-commutative binary operator where left is a parameter
337                // and right is a variable
338                Scalar d1 = 3. / x[3];
339                original += 1;
340                opt      += 1;
341
342                // non-commutative binary operator where left is a variable
343                // and right is a variable
344                Scalar e1 = pow(x[3], x[4]);
345                original += 3;
346                opt      += 3;
347
348                // commutative binary operator where  left is a variable
349                // and right is a parameter
350                Scalar f1 = x[5] * 5.;
351                original += 1;
352                opt      += 1;
353
354                // commutative binary operator where  left is a variable
355                // and right is a variable
356                Scalar g1 = x[5] + x[6];
357                original += 1;
358                opt      += 1;
359
360                // duplicate variables
361                Scalar a2 = CppAD::exp(x[0]);
362                Scalar b2 = CppAD::sin(x[1]);  // counts for 2 variables
363                Scalar c2 = x[2] - 3.;
364                Scalar d2 = 3. / x[3];
365                Scalar e2 = pow(x[3], x[4]);   // counts for 3 variables
366                Scalar f2 = 5. * x[5];
367                Scalar g2 = x[6] + x[5];
368                original += 10;
369
370                // result vector
371                y[0] = a1 * a2;
372                y[1] = b1 * b2;
373                y[2] = c1 * c2;
374                y[3] = d1 * d2;
375                y[4] = e1 * e2;
376                y[5] = f1 * f2;
377                y[6] = g1 * g2;
378                original += 7;
379                opt      += 7;
380
381                return;
382        }
383        bool duplicate_one(void)
384        {
385                bool ok = true;
386                using CppAD::AD;
387                size_t original;
388                size_t opt;
389                size_t i, j;
390       
391                // domain space vector
392                size_t n  = 7;
393                CppAD::vector< AD<double> > X(n);
394                for(j = 0; j < n; j++)
395                        X[j] = 1. / double(j + 1); 
396       
397                // declare independent variables and start tape recording
398                CppAD::Independent(X);
399       
400                // range space vector
401                size_t m = n;
402                CppAD::vector< AD<double> > Y(m);
403                duplicate_fun(X, Y, original, opt);
404       
405                // create f: X -> Y and stop tape recording
406                CppAD::ADFun<double> F;
407                F.Dependent(X, Y); 
408       
409                CppAD::vector<double> x(n), y(m), check(m);
410                for(j = 0; j < n; j++)
411                        x[j] = Value(X[j]);
412                y = F.Forward(0, x);
413                duplicate_fun(x, check, original, opt);
414                for(i = 0; i < m; i++)
415                        ok &= (y[i] == check[i]);
416       
417                // Check size before optimization
418                ok &= F.size_var() == (n + 1 + original);
419       
420                // Optimize the operation sequence
421                F.optimize();
422       
423                // Check size after optimization
424                ok &= F.size_var() == (n + 1 + opt);
425       
426                // check result now
427                // (should have already been checked if NDEBUG not defined)
428                y = F.Forward(0, x);
429                for(i = 0; i < m; i++)
430                        ok &= (y[i] == check[i]);
431       
432                return ok;
433        }
434        // -------------------------------------------------------------------
435        bool duplicate_two(void)
436        {       // test that duplicate expression removal is relative to
437                // new and not just old argument indices.
438                bool ok = true;
439                using CppAD::AD;
440                size_t i, j;
441
442                // domain space vector
443                size_t n  = 1;
444                CppAD::vector< AD<double> > X(n);
445                for(j = 0; j < n; j++)
446                        X[j] = double(j + 2); 
447
448                // range space vector
449                size_t m = 1;
450                CppAD::vector< AD<double> > Y(m);
451
452                // declare independent variables and start tape recording
453                CppAD::Independent(X);
454
455                // create a new variable
456                AD<double> A1 = X[0] - 2.;
457
458                // create a duplicate variable
459                AD<double> A2 = X[0] - 2.;
460
461                // create a new variable using first version of duplicate
462                AD<double> B1 = A1 / 2.;
463
464                // create a duplicate that can only be dectected using new
465                // argument indices
466                AD<double> B2 = A2 / 2.; 
467
468                // Make a new variable for result
469                // and make it depend on all the variables
470                Y[0] = B1 + B2;
471
472                // create f: X -> Y and stop tape recording
473                CppAD::ADFun<double> F;
474                F.Dependent(X, Y); 
475
476                // check number of variables in original function
477                ok &= (F.size_var() ==  1 + n + m + 4 ); 
478       
479                CppAD::vector<double> x(n), y(m);
480                for(j = 0; j < n; j++)
481                        x[j] = double(j + 2);
482
483                y   = F.Forward(0, x);
484                for(i = 0; i < m; i++)
485                        ok &= ( y[i] == Value( Y[i] ) );
486
487                F.optimize();
488
489                // check number of variables  in optimized version
490                ok &= (F.size_var() == 1 + n + m + 2 ); 
491
492                y   = F.Forward(0, x);
493                for(i = 0; i < m; i++)
494                        ok &= ( y[i] == Value( Y[i] ) );
495
496                return ok;
497        }
498        // -------------------------------------------------------------------
499        bool duplicate_three(void)
500        {       // test that duplicate expression removal is relative to
501                // new and not just old argument indices (commutative case).
502                bool ok = true;
503                using CppAD::AD;
504                size_t i, j;
505
506                // domain space vector
507                size_t n  = 1;
508                CppAD::vector< AD<double> > X(n);
509                for(j = 0; j < n; j++)
510                        X[j] = double(j + 2); 
511
512                // range space vector
513                size_t m = 1;
514                CppAD::vector< AD<double> > Y(m);
515
516                // declare independent variables and start tape recording
517                CppAD::Independent(X);
518
519                // create a new variable
520                AD<double> A1 = X[0] + 2.;
521
522                // create a duplicate variable
523                AD<double> A2 = 2. + X[0];
524
525                // create a new variable using first version of duplicate
526                AD<double> B1 = A1 * 2.;
527
528                // create a duplicate that can only be dectected using new
529                // argument indices
530                AD<double> B2 = 2. * A2; 
531
532                // Make a new variable for result
533                // and make it depend on all the variables
534                Y[0] = B1 + B2;
535
536                // create f: X -> Y and stop tape recording
537                CppAD::ADFun<double> F;
538                F.Dependent(X, Y); 
539
540                // check number of variables in original function
541                ok &= (F.size_var() ==  1 + n + m + 4 ); 
542       
543                CppAD::vector<double> x(n), y(m);
544                for(j = 0; j < n; j++)
545                        x[j] = double(j + 2);
546
547                y   = F.Forward(0, x);
548                for(i = 0; i < m; i++)
549                        ok &= ( y[i] == Value( Y[i] ) );
550
551                F.optimize();
552
553                // check number of variables  in optimized version
554                ok &= (F.size_var() == 1 + n + m + 2 ); 
555
556                y   = F.Forward(0, x);
557                for(i = 0; i < m; i++)
558                        ok &= ( y[i] == Value( Y[i] ) );
559
560                return ok;
561        }
562        // -------------------------------------------------------------------
563        bool duplicate_four(void)
564        {       // Check that unary expression matching not only checks hash code,
565                // and operator, but also operand (old bug that has been fixed).
566                bool ok = true;
567                using CppAD::AD;
568                size_t j;
569
570                // domain space vector
571                size_t n  = 1;
572                CppAD::vector< AD<double> > X(n);
573                X[0] = 1.;
574
575                // range space vector
576                size_t m = 1;
577                CppAD::vector< AD<double> > Y(m);
578
579                // declare independent variables and start tape recording
580                CppAD::Independent(X);
581
582                // check a huge number of same operation with different operands
583                size_t n_operations = std::min(
584                        size_t(CPPAD_HASH_TABLE_SIZE) + 5,
585                        size_t(std::numeric_limits<CPPAD_TAPE_ADDR_TYPE>::max()) - 5
586                );
587                Y[0] = X[0];
588                for(j = 0; j < n_operations; j++)
589                        Y[0] = abs(Y[0]);
590
591                // create f: X -> Y and stop tape recording
592                CppAD::ADFun<double> F;
593                F.Dependent(X, Y); 
594
595                // check number of variables in original function
596                ok &= (F.size_var() ==  1 + n + n_operations ); 
597       
598                CppAD::vector<double> x(n), y(m);
599                x[0] = 1.;
600
601                y   = F.Forward(0, x);
602                ok &= ( y[0] == Value( Y[0] ) );
603
604                F.optimize();
605
606                // check same number of variables in optimized version
607                ok &= (F.size_var() == 1 + n + n_operations ); 
608
609                y   = F.Forward(0, x);
610                ok &= ( y[0] == Value( Y[0] ) );
611
612                return ok;
613        }
614        // ====================================================================
615        bool cummulative_sum(void)
616        {       // test conversion of a sequence of additions and subtraction
617                // to a cummulative summation sequence.
618                bool ok = true;
619                using CppAD::AD;
620                size_t i, j;
621
622                // domain space vector
623                size_t n  = 7;
624                CppAD::vector< AD<double> > X(n);
625                for(j = 0; j < n; j++)
626                        X[j] = double(j + 2); 
627
628                size_t n_original = 1 + n;
629                size_t n_optimize = 1 + n;
630
631                // range space vector
632                size_t m = 2;
633                CppAD::vector< AD<double> > Y(m);
634
635                // declare independent variables and start tape recording
636                CppAD::Independent(X);
637
638                // Operations inside of optimize_cadd
639                Y[0] = 5. + (X[0] + X[1]) + (X[1] - X[2]) // Addvv, Subvv
640                     + (X[2] - 1.) + (2. - X[3])   // Subvp, Subpv
641                     + (X[4] + 3.) + (4. + X[5]);  // Addpv, Addpv (no Addvp)
642                n_original += 12;
643                n_optimize += 1;
644
645
646                // Operations inside of optimize_csub
647                Y[1] = 5. - (X[1] + X[2]) - (X[2] - X[3]) // Addvv, Subvv
648                     - (X[3] - 1.) - (2. - X[4])   // Subvp, Subpv
649                     - (X[5] + 3.) - (4. + X[6]);  // Addpv, Addpv (no Addvp)
650                n_original += 12;
651                n_optimize += 1;
652
653                CppAD::ADFun<double> F;
654                F.Dependent(X, Y); 
655
656                // check number of variables in original function
657                ok &= (F.size_var() ==  n_original ); 
658       
659                CppAD::vector<double> x(n), y(m);
660                for(j = 0; j < n; j++)
661                        x[j] = double(j + 2);
662
663                y   = F.Forward(0, x);
664                for(i = 0; i < m; i++)
665                        ok &= ( y[i] == Value( Y[i] ) );
666
667                F.optimize();
668
669                // check number of variables  in optimized version
670                ok &= (F.size_var() == n_optimize ); 
671
672                y   = F.Forward(0, x);
673                for(i = 0; i < m; i++)
674                        ok &= ( y[i] == Value( Y[i] ) );
675
676                return ok;
677        }
678        // -------------------------------------------------------------------
679        bool forward_csum(void)
680        {       bool ok = true;
681       
682                using namespace CppAD;
683       
684                // independent variable vector
685                CppAD::vector< AD<double> > X(2);
686                X[0] = 0.; 
687                X[1] = 1.;
688                Independent(X);
689       
690                // compute sum of elements in X
691                CppAD::vector< AD<double> > Y(1);
692                Y[0] = X[0] + X[0] + X[1];
693       
694                // create function object F : X -> Y
695                ADFun<double> F(X, Y);
696       
697                // now optimize the operation sequence
698                F.optimize();
699       
700                // use zero order to evaluate F[ (3, 4) ]
701                CppAD::vector<double>  x0( F.Domain() );
702                CppAD::vector<double>  y0( F.Range() );
703                x0[0]    = 3.;
704                x0[1]    = 4.;
705                y0       = F.Forward(0, x0);
706                ok      &= NearEqual(y0[0] , x0[0]+x0[0]+x0[1], 1e-10, 1e-10);
707       
708                // evaluate derivative of F in X[0] direction
709                CppAD::vector<double> x1( F.Domain() );
710                CppAD::vector<double> y1( F.Range() );
711                x1[0]    = 1.;
712                x1[1]    = 0.;
713                y1       = F.Forward(1, x1);
714                ok      &= NearEqual(y1[0] , x1[0]+x1[0]+x1[1], 1e-10, 1e-10);
715       
716                // evaluate second derivative of F in X[0] direction
717                CppAD::vector<double> x2( F.Domain() );
718                CppAD::vector<double> y2( F.Range() );
719                x2[0]       = 0.;
720                x2[1]       = 0.;
721                y2          = F.Forward(2, x2);
722                double F_00 = 2. * y2[0];
723                ok         &= NearEqual(F_00, 0., 1e-10, 1e-10);
724       
725                return ok;
726        }
727        // -------------------------------------------------------------------
728        bool reverse_csum(void)
729        {       bool ok = true;
730       
731                using namespace CppAD;
732       
733                // independent variable vector
734                CppAD::vector< AD<double> > X(2);
735                X[0] = 0.; 
736                X[1] = 1.;
737                Independent(X);
738       
739                // compute sum of elements in X
740                CppAD::vector< AD<double> > Y(1);
741                Y[0] = X[0] - (X[0] - X[1]);
742       
743                // create function object F : X -> Y
744                ADFun<double> F(X, Y);
745       
746                // now optimize the operation sequence
747                F.optimize();
748       
749                // use zero order to evaluate F[ (3, 4) ]
750                CppAD::vector<double>  x0( F.Domain() );
751                CppAD::vector<double>  y0( F.Range() );
752                x0[0]    = 3.;
753                x0[1]    = 4.;
754                y0       = F.Forward(0, x0);
755                ok      &= NearEqual(y0[0] , x0[0]-x0[0]+x0[1], 1e-10, 1e-10);
756       
757                // evaluate derivative of F
758                CppAD::vector<double> dF( F.Domain() );
759                CppAD::vector<double> w( F.Range() );
760                w[0]    = 1.;
761                dF      = F.Reverse(1, w);
762                ok     &= NearEqual(dF[0] , 0., 1e-10, 1e-10);
763                ok     &= NearEqual(dF[1] , 1., 1e-10, 1e-10);
764       
765                return ok;
766        }
767        bool forward_sparse_jacobian_csum()
768        {       bool ok = true;
769                using namespace CppAD;
770       
771                // dimension of the domain space
772                size_t n = 3; 
773       
774                // dimension of the range space
775                size_t m = 2;
776       
777                // independent variable vector
778                CppAD::vector< AD<double> > X(n);
779                X[0] = 2.; 
780                X[1] = 3.;
781                Independent(X);
782       
783                // dependent variable vector
784                CppAD::vector< AD<double> > Y(m);
785       
786                // check results vector
787                CppAD::vector< bool >       Check(m * n);
788       
789                // initialize index into Y
790                size_t index = 0;
791       
792                // Y[0]
793                Y[index]             = X[0] + X[1] + 5.;
794                Check[index * n + 0] = true;
795                Check[index * n + 1] = true;
796                Check[index * n + 2] = false;
797                index++;
798       
799                // Y[1]
800                Y[index]             = Y[0] - (X[1] + X[2]);
801                Check[index * n + 0] = true;
802                Check[index * n + 1] = true;
803                Check[index * n + 2] = true;
804                index++;
805       
806                // check final index
807                assert( index == m );
808       
809                // create function object F : X -> Y
810                ADFun<double> F(X, Y);
811                F.optimize();
812       
813                // ---------------------------------------------------------
814                // dependency matrix for the identity function
815                CppAD::vector< std::set<size_t> > Sx(n);
816                size_t i, j;
817                for(i = 0; i < n; i++)
818                {       assert( Sx[i].empty() );
819                        Sx[i].insert(i);
820                }
821       
822                // evaluate the dependency matrix for F(x)
823                CppAD::vector< std::set<size_t> > Sy(m);
824                Sy = F.ForSparseJac(n, Sx);
825       
826                // check values
827                bool found;
828                for(i = 0; i < m; i++)
829                {       for(j = 0; j < n; j++)
830                        {       found = Sy[i].find(j) != Sy[i].end();
831                                ok &= (found == Check[i * n + j]);
832                        }
833                }       
834       
835                return ok;
836        }
837        bool reverse_sparse_jacobian_csum()
838        {       bool ok = true;
839                using namespace CppAD;
840       
841                // dimension of the domain space
842                size_t n = 3; 
843       
844                // dimension of the range space
845                size_t m = 2;
846       
847                // independent variable vector
848                CppAD::vector< AD<double> > X(n);
849                X[0] = 2.; 
850                X[1] = 3.;
851                Independent(X);
852       
853                // dependent variable vector
854                CppAD::vector< AD<double> > Y(m);
855       
856                // check results vector
857                CppAD::vector< bool >       Check(m * n);
858       
859                // initialize index into Y
860                size_t index = 0;
861       
862                // Y[0]
863                Y[index]             = X[0] + X[1] + 5.;
864                Check[index * n + 0] = true;
865                Check[index * n + 1] = true;
866                Check[index * n + 2] = false;
867                index++;
868       
869                // Y[1]
870                Y[index]             = Y[0] - (X[1] + X[2]);
871                Check[index * n + 0] = true;
872                Check[index * n + 1] = true;
873                Check[index * n + 2] = true;
874                index++;
875       
876                // check final index
877                assert( index == m );
878       
879                // create function object F : X -> Y
880                ADFun<double> F(X, Y);
881                F.optimize();
882       
883                // ----------------------------------------------------------
884                // dependency matrix for the identity function
885                CppAD::vector< bool > Py(m * m);
886                size_t i, j;
887                for(i = 0; i < m; i++)
888                {       for(j = 0; j < m; j++)
889                                Py[ i * m + j ] = (i == j);
890                }
891       
892                // evaluate the dependency matrix for F(x)
893                CppAD::vector< bool > Px(m * n);
894                Px = F.RevSparseJac(m, Py);
895       
896                // check values
897                for(i = 0; i < m; i++)
898                {       for(j = 0; j < n; j++)
899                                ok &= (Px[i * n + j] == Check[i * n + j]);
900                }       
901       
902                return ok;
903        }
904        bool reverse_sparse_hessian_csum(void)
905        {       bool ok = true;
906                using CppAD::AD;
907                size_t i, j;
908       
909                size_t n = 2;
910                CppAD::vector< AD<double> > X(n); 
911                X[0] = 1.;
912                X[1] = 2.;
913                CppAD::Independent(X);
914       
915                size_t m = 1;
916                CppAD::vector< AD<double> > Y(m);
917                Y[0] = (2. + X[0] + X[1] + 3.) * X[0];
918       
919                CppAD::vector<bool> check(n * n);
920                check[0 * n + 0] = true;  // partial w.r.t. x[0], x[0]
921                check[0 * n + 1] = true;  //                x[0], x[1]
922                check[1 * n + 0] = true;  //                x[1], x[0]
923                check[1 * n + 1] = false; //                x[1], x[1]
924       
925                // create function object F : X -> Y
926                CppAD::ADFun<double> F(X, Y);
927                F.optimize();
928       
929                // sparsity pattern for the identity function U(x) = x
930                CppAD::vector<bool> Px(n * n);
931                for(i = 0; i < n; i++)
932                        for(j = 0; j < n; j++)
933                                Px[ i * n + j ] = (i == j);
934       
935                // compute sparsity pattern for Jacobian of F(U(x))
936                CppAD::vector<bool> P_jac(m * n);
937                P_jac = F.ForSparseJac(n, Px);
938       
939                // compute sparsity pattern for Hessian of F_k ( U(x) )
940                CppAD::vector<bool> Py(m);
941                CppAD::vector<bool> Pxx(n * n);
942                Py[0] = true;
943                Pxx = F.RevSparseHes(n, Py);
944                // check values
945                for(i = 0; i < n * n; i++)
946                        ok &= (Pxx[i] == check[i]);
947       
948                return ok;
949        }
950        // check that CondExp properly detects dependencies
951        bool cond_exp_depend(void)
952        {       bool ok = true;
953                using CppAD::AD;
954
955                AD<double> zero(0.), one(1.), two(2.), three(3.);
956       
957                size_t n = 4;
958                CppAD::vector< AD<double> > X(n); 
959                X[0] = zero;
960                X[1] = one;
961                X[2] = two;
962                X[3] = three;
963                CppAD::Independent(X);
964       
965                size_t m = 4;
966                CppAD::vector< AD<double> > Y(m);
967                Y[0] = CondExpLt(X[0] + .5,  one,  two, three);
968                Y[1] = CondExpLt(zero, X[1] + .5,  two, three);
969                Y[2] = CondExpLt(zero,  one, X[2] + .5, three);
970                Y[3] = CondExpLt(zero,  one,  two,  X[3] + .5);
971
972                CppAD::ADFun<double> f(X, Y);
973                f.optimize();
974
975                CppAD::vector<double> x(n), y(m);
976                size_t i;
977                for(i = 0; i < n; i++)
978                        x[i] = double(n - i);
979                y    = f.Forward(0, x);
980
981                if( x[0] + .5 < 1. )
982                        ok  &= y[0] == 2.;
983                else    ok  &= y[0] == 3.;
984                if( 0. < x[1] + .5 )
985                        ok  &= y[1] == 2.;
986                else    ok  &= y[1] == 3.;
987                ok  &= y[2] == x[2] + .5;;
988                ok  &= y[3] == 2.;
989
990                return ok;
991        }
992        // -------------------------------------------------------------------
993        void my_union(
994                std::set<size_t>&         result  ,
995                const std::set<size_t>&   left    ,
996                const std::set<size_t>&   right   )
997        {       std::set<size_t> temp;
998                std::set_union(
999                        left.begin()              ,
1000                        left.end()                ,
1001                        right.begin()             ,
1002                        right.end()               ,
1003                        std::inserter(temp, temp.begin())
1004                );
1005                result.swap(temp);
1006        }
1007
1008        bool old_atomic_forward(
1009                size_t                         id ,
1010                size_t                          k , 
1011                size_t                          n ,
1012                size_t                          m ,
1013                const CppAD::vector<bool>&     vx ,
1014                CppAD::vector<bool>&           vy ,
1015                const CppAD::vector<double>&   tx , 
1016                CppAD::vector<double>&         ty )
1017        {       assert(n == 3 && m == 2);
1018                if( k > 0 ) 
1019                        return false;
1020
1021                // y[0] = x[0] + x[1]
1022                ty[0] = tx[0] + tx[1];
1023
1024                // y[1] = x[1] + x[2]
1025                ty[1] = tx[1] + tx[2];
1026               
1027                if( vy.size() > 0 )
1028                {       vy[0] = (vx[0] | vx[1]);
1029                        vy[1] = (vx[1] | vx[2]);
1030                }
1031                return true; 
1032        }
1033
1034        bool old_atomic_reverse(
1035                size_t                         id ,
1036                size_t                          k , 
1037                size_t                          n , 
1038                size_t                          m , 
1039                const CppAD::vector<double>&   tx , 
1040                const CppAD::vector<double>&   ty ,
1041                CppAD::vector<double>&         px ,
1042                const CppAD::vector<double>&   py )
1043        {       return false; }
1044
1045        bool old_atomic_for_jac_sparse(
1046                size_t                                  id ,
1047                size_t                                   n ,
1048                size_t                                   m ,
1049                size_t                                   q ,
1050                const CppAD::vector< std::set<size_t> >& r ,
1051                CppAD::vector< std::set<size_t>  >&      s )
1052        {       return false; }
1053
1054        bool old_atomic_rev_jac_sparse(
1055                size_t                                  id ,
1056                size_t                                   n ,
1057                size_t                                   m ,
1058                size_t                                   q ,
1059                CppAD::vector< std::set<size_t> >&       r ,
1060                const CppAD::vector< std::set<size_t> >& s )
1061        {       assert(n == 3 && m == 2);
1062                r[0].clear();
1063                r[1].clear();
1064                r[2].clear();
1065                // y[0] = x[0] + x[1]
1066                my_union(r[0], r[0], s[0]);
1067                my_union(r[1], r[1], s[0]);
1068                // y[1] = x[1] + x[2]
1069                my_union(r[1], r[1], s[1]);
1070                my_union(r[2], r[2], s[1]);
1071
1072                return true; 
1073        }
1074
1075        bool old_atomic_rev_hes_sparse(
1076                size_t                                  id ,
1077                size_t                                   n ,
1078                size_t                                   m ,
1079                size_t                                   q ,
1080                const CppAD::vector< std::set<size_t> >& r ,
1081                const CppAD::vector<bool>&               s ,
1082                CppAD::vector<bool>&                     t ,
1083                const CppAD::vector< std::set<size_t> >& u ,
1084                CppAD::vector< std::set<size_t> >&       v )
1085        {       return false; }
1086
1087        CPPAD_USER_ATOMIC(
1088                my_old_atomic             ,
1089                CppAD::vector              ,
1090                double                     ,
1091                old_atomic_forward        ,
1092                old_atomic_reverse        ,
1093                old_atomic_for_jac_sparse ,
1094                old_atomic_rev_jac_sparse ,
1095                old_atomic_rev_hes_sparse
1096        )
1097
1098        bool old_atomic_test(void)
1099        {       bool ok = true;
1100
1101                using CppAD::AD;
1102                size_t j;
1103                size_t n = 3;
1104                size_t m = 2;
1105                CppAD::vector< AD<double> > ax(n), ay(m), az(m);
1106                for(j = 0; j < n; j++)
1107                        ax[j] = AD<double>(j + 1);
1108                CppAD::Independent(ax);
1109
1110                size_t id = 0;
1111                // first call should stay in the tape
1112                my_old_atomic(id++, ax, ay);
1113                // second call will not get used
1114                my_old_atomic(id++, ax, az);
1115                // create function
1116                CppAD::ADFun<double> g(ax, ay);
1117                // should have 1 + n + m + m varaibles
1118                ok &= g.size_var() == (1 + n + m + m);
1119                g.optimize();
1120                // should have 1 + n + m varaibles
1121                ok &= g.size_var() == (1 + n + m);
1122
1123                // now test that the optimized function gives same results
1124                CppAD::vector<double> x(n), y(m);
1125                for(j = 0; j < n; j++)
1126                        x[j] = (j + 1) * (j + 1);
1127                y = g.Forward(0, x);
1128                // y[0] = x[0] + x[1]
1129                ok &= (y[0] == x[0] + x[1]);
1130                // y[1] = x[1] + x[2]
1131                ok &= (y[0] == x[0] + x[1]);
1132
1133                return ok;
1134        }
1135        bool not_identically_equal(void)
1136        {       bool ok = true;
1137                using CppAD::AD;
1138
1139                // independent variable vector
1140                size_t n = 5;
1141                CppAD::vector< AD<double> > ax(n);
1142                size_t j;
1143                for(j = 0; j < n; j++)
1144                        ax[j] = 1. / 3.;
1145                CppAD::Independent(ax);
1146       
1147                // dependent variable vector
1148                size_t m = 1;
1149                CppAD::vector< AD<double> > ay(m);
1150                ay[0]       = 0.;
1151                for(j = 0; j < n; j++)
1152                {       if( j % 2 == 0 )
1153                                ay[0] += ax[j];
1154                        else    ay[0] -= ax[j];
1155                }
1156                CppAD::ADFun<double> f(ax, ay);
1157       
1158                // Used to fail assert in optimize that forward mode results
1159                // are identically equal
1160                f.optimize();
1161       
1162                return ok;
1163        }
1164}
1165
1166bool optimize(void)
1167{       bool ok = true;
1168        // check optimizing out atomic arguments
1169        ok     &= atomic_arguments();
1170        // check reverse dependency analysis optimization
1171        ok     &= depend_one();
1172        ok     &= depend_two();
1173        ok     &= depend_three();
1174        // check removal of duplicate expressions
1175        ok     &= duplicate_one();
1176        ok     &= duplicate_two();
1177        ok     &= duplicate_three();
1178        ok     &= duplicate_four();
1179        // convert sequence of additions to cummulative summation
1180        ok     &= cummulative_sum();
1181        ok     &= forward_csum();
1182        ok     &= reverse_csum();
1183        ok     &= forward_sparse_jacobian_csum();
1184        ok     &= reverse_sparse_jacobian_csum();
1185        ok     &= reverse_sparse_hessian_csum();
1186        // check that CondExp properly detects dependencies
1187        ok     &= cond_exp_depend();
1188        // check old_atomic functions
1189        ok     &= old_atomic_test();
1190        // case where results are not identically equal
1191        ok     &= not_identically_equal();
1192
1193        CppAD::user_atomic<double>::clear();
1194        return ok;
1195}
Note: See TracBrowser for help on using the repository browser.