source: trunk/test_more/optimize.cpp @ 1533

Last change on this file since 1533 was 1497, checked in by bradbell, 11 years ago

trunk: Merge in branches/optimize with minor chages to following files:

cppad/local/optimize.hpp
cppad/local/cap_taylor.hpp
omh/forward.omh
omh/seq_property.omh
omh/whats_new_09.omh

  • Property svn:keywords set to Id
File size: 4.8 KB
Line 
1/* $Id: optimize.cpp 1497 2009-08-13 16:57:34Z bradbell $ */
2/* --------------------------------------------------------------------------
3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-09 Bradley M. Bell
4
5CppAD is distributed under multiple licenses. This distribution is under
6the terms of the
7                    Common 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 <cppad/cppad.hpp>
14
15namespace {
16
17        template <class Vector>
18        void fun(const Vector& x, Vector& y, size_t& original, size_t& opt)
19        {       typedef typename Vector::value_type Scalar;
20                Scalar a;
21                Scalar one(1), two(2), three(3), four(4);
22                original = 0;
23                opt      = 0;
24
25                // unary operator where operand is arg[0]
26                a = abs(x[0]); 
27                if( a < 1. )
28                        y[0] = sin(x[0]);
29                else    y[0] = cos(x[0]); 
30                original += 3;
31                opt      += 2;
32
33                // binary operator where left operand is a variable
34                // and right operand is a parameter
35                a = x[1] + 2.;
36                if( a < 3. )
37                        y[1] = x[1] * 3.;
38                else    y[1] = x[1] / 2.;
39                original += 2;
40                opt      += 1;
41
42                // binary operator where left operand is a parameter
43                // and right operation is a variable
44                a = 2. - x[2];
45                if( a < 4. )
46                        y[2] = 3. / x[2];
47                else    y[2] = 4. + x[2];
48                original += 2;
49                opt      += 1;
50
51                // binary operator where both operands are variables
52                a = x[3] - x[2];
53                if( a < 4. )
54                        y[3] = x[3] / x[2];
55                else    y[3] = x[3] + x[2];
56                original += 2;
57                opt      += 1;
58
59                // a conditional expression that will be optimized out
60                a = CppAD::CondExpLt(x[0], x[1], x[2], x[3]);
61                if( a < 5. )
62                        y[4] = CppAD::CondExpLt(x[4], one, two, three);
63                else    y[4] = CppAD::CondExpLt(x[4], two, three, four);
64                original += 2;
65                opt      += 1;
66
67                // Make sure that a parameter dependent variable
68                // is not optimized out of the operation sequence.
69                // In addition, we do not use the argument x[5], to
70                // make sure it is not optimized out.
71                y[5] = 1.;
72                original += 1;
73                opt      += 1;
74
75                return;
76        }
77
78        bool CaseOne(void)
79        {       // Test all except for VecAD operations
80                bool ok = true;
81                using CppAD::AD;
82                size_t original;
83                size_t opt;
84                size_t i, j;
85       
86                // domain space vector
87                size_t n  = 6;
88                CPPAD_TEST_VECTOR< AD<double> > X(n);
89                for(j = 0; j < n; j++)
90                        X[j] = 1. / double(j + 1); 
91       
92                // declare independent variables and start tape recording
93                CppAD::Independent(X);
94       
95                // range space vector
96                size_t m = n;
97                CPPAD_TEST_VECTOR< AD<double> > Y(m);
98                fun(X, Y, original, opt);
99       
100                // create f: X -> Y and stop tape recording
101                CppAD::ADFun<double> F;
102                F.Dependent(X, Y); 
103       
104                CPPAD_TEST_VECTOR<double> x(n), y(m), check(m);
105                for(j = 0; j < n; j++)
106                        x[j] = Value(X[j]);
107                y = F.Forward(0, x);
108                fun(x, check, original, opt);
109                for(i = 0; i < m; i++)
110                        ok &= (y[i] == check[i]);
111       
112                // Check size before optimization
113                ok &= F.size_var() == (n + 1 + original);
114       
115                // Optimize the operation sequence
116                F.optimize();
117       
118                // Check size after optimization
119                ok &= F.size_var() == (n + 1 + opt);
120       
121                // check result now
122                // (should have already been checked if NDEBUG not defined)
123                y = F.Forward(0, x);
124                for(i = 0; i < m; i++)
125                        ok &= (y[i] == check[i]);
126       
127                return ok;
128        }
129
130        bool CaseTwo(void)
131        {       // Test VecAD operations
132                bool ok = true;
133                using CppAD::AD;
134                size_t i, j;
135
136                // domain space vector
137                size_t n  = 2;
138                CPPAD_TEST_VECTOR< AD<double> > X(n);
139                for(j = 0; j < n; j++)
140                        X[j] = double(j); 
141
142                // range space vector
143                size_t m = 3;
144                CPPAD_TEST_VECTOR< AD<double> > Y(m);
145
146                CppAD::VecAD<double> U(n);
147                CppAD::VecAD<double> V(m);
148                for(j = 0; j < n; j++)
149                        U[j] = 0;
150                for(i = 0; i < m; i++)
151                        V[i] = 0;
152       
153                // declare independent variables and start tape recording
154                CppAD::Independent(X);
155
156                // make U a variable
157                U[ X[0] ] = X[1];
158
159                // make V a variable
160                V[ X[0] ] = X[1];
161
162                // Y only depend on V (and not on U)
163                for(i = 0; i < m; i++)
164                {       AD<double> I(i);
165                        Y[i] = V[I];
166                }
167       
168                // create f: X -> Y and stop tape recording
169                // Y[ X[0] ] = X[1] and other components of Y are zero.
170                CppAD::ADFun<double> F;
171                F.Dependent(X, Y); 
172
173                // Check number of VecAD vectors plus number of VecAD elements
174                ok &= (F.size_VecAD() == 2 + n + m); 
175       
176                CPPAD_TEST_VECTOR<double> x(n), y(m);
177                for(j = 0; j < n; j++)
178                        x[j] = double(j);
179
180                y = F.Forward(0, x);
181                for(i = 0; i < m; i++)
182                {       if( i != static_cast<size_t>(x[0]) )
183                                ok &= (y[i] == 0.);
184                        else    ok &= (y[i] == x[1]);
185                }
186
187                F.optimize();
188
189                // Check number of VecAD vectors plus number of VecAD elements
190                ok &= (F.size_VecAD() == 1 + m); 
191                y = F.Forward(0, x);
192                for(i = 0; i < m; i++)
193                {       if( i != static_cast<size_t>(x[0]) )
194                                ok &= (y[i] == 0.);
195                        else    ok &= (y[i] == x[1]);
196                }
197               
198                return ok;
199        }
200}
201
202bool optimize(void)
203{       bool ok = true;
204        ok     &= CaseOne();
205        ok     &= CaseTwo();
206
207        return ok;
208}
Note: See TracBrowser for help on using the repository browser.