source: trunk/test_more/forward_mul.cpp @ 3008

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

merge in changes from branches/atomic; see bin/svn_merge.sh

  • Property svn:keywords set to Id
File size: 5.6 KB
Line 
1/* $Id: forward_mul.cpp 2859 2013-05-28 06:03:21Z 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
13/*
14$begin forward_mul.cpp$$
15$spell
16        Cpp
17$$
18
19$section Forward Mode: Example and Test of Multiple Orders$$
20$index forward, multiple orders$$
21$index multiple, forward orders$$
22$index order, multiple forward$$
23$index example, forward multiple orders$$
24$index test, forward multiple orders$$
25
26$code
27$verbatim%example/forward_mul.cpp%0%// BEGIN C++%// END C++%1%$$
28$$
29
30$end
31*/
32// BEGIN C++
33# include <cppad/cppad.hpp>
34
35namespace {
36
37        double my_discrete(const double& x)
38        {       return static_cast<int> ( x ); }
39        CPPAD_DISCRETE_FUNCTION(double, my_discrete)
40
41}
42bool forward_mul(void)
43{       bool ok = true;
44        using CppAD::AD;
45        using CppAD::NearEqual;
46        size_t j, k;
47        double eps = 10. * CppAD::numeric_limits<double>::epsilon();
48
49        // domain space vector
50        size_t n = 23, m = n;
51        CPPAD_TESTVECTOR(AD<double>) X(n), Y(m);
52        for(j = 0; j < n; j++)
53                X[j] = 0.0;
54
55        // declare independent variables and starting recording
56        CppAD::Independent(X);
57
58        // identity function values
59        size_t i = 0;
60        size_t identity_begin = i;
61        Y[i] = cos( acos( X[i] ) );                   i++; // AcosOp,  CosOp
62        Y[i] = sin( asin( X[i] ) );                   i++; // AsinOp,  SinOp
63        Y[i] = tan( atan( X[i] ) );                   i++; // AtanOp,  TanOp
64        Y[i] = CondExpGt(X[i], X[i-1], X[i], X[i-2]); i++; // CExpOp
65        Y[i] = X[i-1] * X[i] / X[i-1];                i++; // DivvvOp, MulvvOp
66        Y[i] = X[i] * X[i] * 1.0 / X[i];              i++; // DivpvOp
67        Y[i] = 5.0 * X[i] / 5.0;                      i++; // DivvpOp, MulpvOp
68        Y[i] = exp( log( X[i] ) );                    i++; // ExpOp,   LogOp
69        Y[i] = pow( sqrt( X[i] ), 2.0);               i++; // PowvpOp, SqrtOp
70        Y[i] = log( pow( std::exp(1.), X[i] ) );      i++; // PowpvOp
71        Y[i] = log( pow( X[i], X[i] ) ) / log( X[i]); i++; // PowvvOp
72        Y[i] = -2. - ((X[i-1] - X[i]) - 2.) + X[i-1]; i++; // Sub*Op: pv, vv, vp
73        size_t identity_end = i;
74
75        // other functions
76        Y[i] = abs( X[i] );         i++;   // AbsOp
77        Y[i] = X[i-1] + X[i] + 2.0; i++;   // AddvvOp, AddvpOp
78        Y[i] = cosh( X[i] );        i++;   // CoshOp
79        Y[i] = my_discrete( X[i] ); i++;   // DisOp
80        Y[i] = 4.0;                 i++;   // ParOp
81        Y[i] = sign( X[i] );        i++;   // SignOp
82        Y[i] = sinh( X[i] );        i++;   // SinhOp
83        Y[i] = tanh(X[i]);          i++;   // TanhOp
84
85        // VecAD operations
86        CppAD::VecAD<double> V(n);
87        AD<double> index = 1.;
88        V[index] = 3.0;
89        Y[i]     = V[index];            i++;   // StppOp, LdpOp
90        V[index] = X[0];
91        Y[i]     = V[index];            i++;   // StpvOp, LdpOp
92        index    = double(n) * X[3];
93        V[index] = X[1];
94        Y[i]     = V[index];            i++;   // StvvOp, LdvOp
95
96        // create f: X -> Y and stop tape recording
97        assert( i == m );
98        CppAD::ADFun<double> f;
99        f.Dependent(X, Y);
100
101        // initially, no values stored in f
102        ok &= f.size_taylor() == 0;
103
104        // Set X_j (t) = x + t
105        size_t p = 2, p1 = p+1;
106        CPPAD_TESTVECTOR(double) x(n), x_p(n * p1), y_p(m * p1);
107        for(j = 0; j < n; j++)
108        {       x[j]            = double(j) / double(n); 
109                x_p[j * p1 + 0] = x[j]; // order 0
110                x_p[j * p1 + 1] = 1.;   // order 1
111                x_p[j * p1 + 2] = 0.;   // order 2
112        }
113        // compute orders 0, 1, 2
114        y_p  = f.Forward(p, x_p);
115
116        // identity functions
117        CPPAD_TESTVECTOR(double) y(p1);
118        i = 0;
119        for(j = identity_begin; j != identity_end; j++)
120        {       y[0] = x[j];
121                y[1] = 1.0;
122                y[2] = 0.0;
123                for(k = 0; k < p1; k++)
124                        ok  &= NearEqual(y[k] , y_p[i * p1 + k], eps, eps);
125                i++;
126        }
127
128        // y_i = abs( x_i )
129        y[0] = CppAD::abs( x[i] );
130        y[1] = CppAD::sign( x[i] );
131        y[2] = 0.0;
132        for(k = 0; k < p1; k++)
133                ok  &= NearEqual(y[k] , y_p[i * p1 + k], eps, eps);
134
135        // y_i = x_[i-1] + x_i + 2
136        i++;
137        y[0] = x[i-1] + x[i] + 2.0;
138        y[1] = 2.0;
139        y[2] = 0.0;
140        for(k = 0; k < p1; k++)
141                ok  &= NearEqual(y[k] , y_p[i * p1 + k], eps, eps);
142
143        // y_i = cosh( x_i )
144        i++;
145        y[0] = CppAD::cosh( x[i] );
146        y[1] = CppAD::sinh( x[i] );
147        y[2] = CppAD::cosh( x[i] ) / 2.0;
148        for(k = 0; k < p1; k++)
149                ok  &= NearEqual(y[k] , y_p[i * p1 + k], eps, eps);
150
151        // y_i = my_discrete( x_i )
152        i++;
153        y[0] = my_discrete( x[i] );
154        y[1] = 0.0;
155        y[2] = 0.0;
156        for(k = 0; k < p1; k++)
157                ok  &= NearEqual(y[k] , y_p[i * p1 + k], eps, eps);
158
159        // y_i = 4
160        i++;
161        y[0] = 4.0;
162        y[1] = 0.0;
163        y[2] = 0.0;
164        for(k = 0; k < p1; k++)
165                ok  &= NearEqual(y[k] , y_p[i * p1 + k], eps, eps);
166
167        // y_i = sign( x_i )
168        i++;
169        y[0] = CppAD::sign( x[i] );
170        y[1] = 0.0;
171        y[2] = 0.0;
172        for(k = 0; k < p1; k++)
173                ok  &= NearEqual(y[k] , y_p[i * p1 + k], eps, eps);
174
175        // y_i = sinh( x_i )
176        i++;
177        y[0] = CppAD::sinh( x[i] );
178        y[1] = CppAD::cosh( x[i] );
179        y[2] = CppAD::sinh( x[i] ) / 2.0;
180        for(k = 0; k < p1; k++)
181                ok  &= NearEqual(y[k] , y_p[i * p1 + k], eps, eps);
182
183        // y_i = tanh( x_i )
184        i++;
185        y[0] = CppAD::tanh( x[i] );
186        y[1] = 1.0 - y[0] * y[0];
187        y[2] = - 2.0 * y[0] * y[1] / 2.0;
188        for(k = 0; k < p1; k++)
189                ok  &= NearEqual(y[k] , y_p[i * p1 + k], eps, eps);
190
191        // y_i = 3.0;
192        i++;
193        y[0] = 3.0;
194        y[1] = 0.0;
195        y[2] = 0.0;
196        for(k = 0; k < p1; k++)
197                ok  &= NearEqual(y[k] , y_p[i * p1 + k], eps, eps);
198
199        // y_i = x_0
200        i++;
201        y[0] = x[0];
202        y[1] = 1.0;
203        y[2] = 0.0;
204        for(k = 0; k < p1; k++)
205                ok  &= NearEqual(y[k] , y_p[i * p1 + k], eps, eps);
206
207        // y_i = x_1
208        i++;
209        y[0] = x[1];
210        y[1] = 1.0;
211        y[2] = 0.0;
212        for(k = 0; k < p1; k++)
213                ok  &= NearEqual(y[k] , y_p[i * p1 + k], eps, eps);
214
215        return ok;
216}
217// END C++
Note: See TracBrowser for help on using the repository browser.