source: trunk/test_more/optimize.cpp @ 2506

Last change on this file since 2506 was 2506, checked in by bradbell, 8 years ago

Change Licenses: CPL-1.0 -> EPL-1.0, GPL-2.0->GPL-3.0

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