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

Last change on this file since 2979 was 2979, checked in by bradbell, 7 years ago
  1. Add operator index to trace output.
  2. Check that new operator values have been set and are valid.
  3. Reinitialize cskip_op_ after optimization.

optimize.hpp: Add test of entirely removing an atomic call.

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