source: trunk/test_more/vec_ad.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: 6.4 KB
Line 
1/* $Id: vec_ad.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/*
14Old examples now only used for validation testing
15*/
16// BEGIN C++
17
18# include <cppad/cppad.hpp>
19# include <cmath>
20
21
22namespace { // Begin empty namespace
23
24bool VecADTestOne(void)
25{       bool ok = true;
26
27        using namespace CppAD;
28        using CppAD::sin;
29        using CppAD::cos;
30
31
32        size_t      n = 3;
33        AD<double>  N(n);
34
35        AD<double>  a;
36        size_t      i;
37
38        // create the array
39        CppAD::VecAD<double> V(n);
40
41        // check assignment from double (while not taping)
42        for(i = 0; i < n; i++)
43                V[i] = double(n - i);
44
45        // check assignment from an AD<double> (while not taping)
46        for(i = 0; i < n; i++)
47                V[i] = 2. * V[i];
48
49        // check array values (while not taping)
50        for(i = 0; i < n; i++)
51                ok &= ( V[i] == 2. * double(n - i) ); 
52
53        // independent variable
54        CPPAD_TESTVECTOR(AD<double>) X(1);
55        X[0] = double(n - 1);
56        Independent(X);
57
58        // check assignment from double during taping
59        a = -1.;
60        for(i = 0; i < n; i++)
61        {       a += 1.;
62                V[a] = double(n - i);
63        }
64
65        // check assignment from AD<double> during taping
66        a = -1.;
67        for(i = 0; i < n; i++)
68        {       a += 1.;
69                V[a] = sin( X[0] ) * V[a];
70        }
71
72        // dependent variable
73        CPPAD_TESTVECTOR(AD<double>) Z(1);
74        Z[0] = V[ X[0] ];
75
76        // create f: X -> Z
77        ADFun<double> f(X, Z);
78        CPPAD_TESTVECTOR(double)  x( f.Domain() );
79        CPPAD_TESTVECTOR(double) dx( f.Domain() );
80        CPPAD_TESTVECTOR(double)  z( f.Range() );
81        CPPAD_TESTVECTOR(double) dz( f.Range() );
82
83        double vx;
84        for(i = 0; i < n; i++)
85        {       // check that the indexing operation was taped
86                x[0] = double(i);
87                z    = f.Forward(0, x);
88                vx   = double(n - i);
89                ok  &= NearEqual(z[0], sin(x[0]) * vx, 1e-10, 1e-10); 
90
91                // note that derivative of v[x] w.r.t. x is zero
92                dx[0] = 1.;
93                dz    = f.Forward(1, dx);
94                ok   &= NearEqual(dz[0], cos(x[0]) * vx, 1e-10, 1e-10); 
95
96                // reverse mode calculation of same value
97                dz[0] = 1.;
98                dx    = f.Reverse(1, dz);
99                ok   &= NearEqual(dx[0], cos(x[0]) * vx, 1e-10, 1e-10); 
100        }
101
102
103        return ok;
104}
105
106// create the discrete function AD<double> Floor(const AD<double> &X)
107double Floor(const double &x)
108{       return std::floor(x); } 
109CPPAD_DISCRETE_FUNCTION(double, Floor)
110
111bool VecADTestTwo(void)
112{       bool ok = true;
113        using namespace CppAD;
114       
115        double pi    = 4. * CppAD::atan(1.);
116        size_t nx    = 10;                       // number of x grid point
117        double xLow  = 0;                        // minimum value for x
118        double xUp   = 2 * pi;                   // maximum value for x
119        double xStep = (xUp - xLow) / (nx - 1);  // step size in x
120        double xCur;                             // current value for x
121
122        // fill in the data vector on a uniform grid
123        VecAD<double> Data(nx);
124        size_t i;
125        for(i = 0; i < nx; i++)
126        {       xCur = xLow + double(i) * xStep; 
127                // use size_t indexing of Data while not taping
128                Data[i] = CppAD::sin(xCur); 
129        }
130
131        // declare independent variable
132        CPPAD_TESTVECTOR(AD<double>) X(1);
133        X[0] = 2.;
134        Independent(X);
135
136        // compute index corresponding to current value of X[0]
137        AD<double> I = X[0] / xStep;
138        AD<double> Ifloor = Floor(I);
139
140        // make sure Ifloor >= 0  (during taping)
141        AD<double> Zero(0);
142        Ifloor = CondExpLt(Ifloor, Zero, Zero, Ifloor);
143
144        // make sure Ifloor <= nx - 2 (during taping)
145        AD<double> Nxminus2(nx - 2);
146        Ifloor = CondExpGt(Ifloor, Nxminus2, Nxminus2, Ifloor);
147
148        // Iceil is Ifloor + 1
149        AD<double> Iceil  = Ifloor + 1.;
150
151        // linear interpolate Data
152        CPPAD_TESTVECTOR(AD<double>) Y(1);
153        Y[0] = Data[Ifloor] + (I - Ifloor) * (Data[Iceil] - Data[Ifloor]);
154
155        // create f: X -> Y that linearly interpolates the data vector
156        ADFun<double> f(X, Y);
157
158        // evaluate the linear interpolant at the mid point for first interval
159        CPPAD_TESTVECTOR(double)  x(1);
160        CPPAD_TESTVECTOR(double)  y(1);
161        x[0] = xStep / 2.;
162        y    = f.Forward(0, x);
163        ok  &= NearEqual(y[0], (Data[0] + Data[1])/2., 1e-10, 1e-10);
164
165        // evalute the derivative with respect to x
166        CPPAD_TESTVECTOR(double) dx(1);
167        CPPAD_TESTVECTOR(double) dy(1);
168        dx[0] = 1.;
169        dy    = f.Forward(1, dx);
170        ok   &= NearEqual(dy[0], (Data[1] - Data[0]) / xStep, 1e-10, 1e-10);
171
172        return ok;
173}
174
175# include <limits>
176bool SecondOrderReverse(void)
177{       // Bradley M. Bell 2009-07-06
178        // Reverse mode for LdpOp was only modifying the highest order partial
179        // This test demonstrated the bug
180        bool ok = true;
181        using CppAD::AD;
182        using CppAD::NearEqual;
183        double eps = 10. * std::numeric_limits<double>::epsilon();
184
185        size_t n = 1;
186        CPPAD_TESTVECTOR(AD<double>) X(n);
187        X[0] = 2.;
188        CppAD::Independent(X);
189
190        size_t m = 2;
191        CPPAD_TESTVECTOR(AD<double>) Y(m);
192
193        // The LdpOp instruction corresponds to operations with VecAD vectors.
194        CppAD::VecAD<double> Z(2);
195        AD<double> zero = 0;
196        Z[zero] = X[0] + 1;
197
198        // The LdvOp instruction corresponds to the index being a variable.
199        AD<double> one = X[0] - 1; // one in a variable here
200        Z[one]  = X[0] + 1.;
201       
202
203        // Compute a function where the second order partial for y
204        // depends on the first order partials for z
205        // This will use the LdpOp instruction because the index
206        // access to z is the parameter zero.
207        Y[0] = Z[zero] * Z[zero];
208        Y[1] = Z[one]  * Z[one];
209
210        CppAD::ADFun<double> f(X, Y);
211
212        // first order forward
213        CPPAD_TESTVECTOR(double) dx(n);
214        size_t p = 1;
215        dx[0]    = 1.;
216        f.Forward(p, dx);
217
218        // Test LdpOp
219        // second order reverse (test exp_if_true case)
220        CPPAD_TESTVECTOR(double) w(m), dw(2 * n);
221        w[0] = 1.;
222        w[1] = 0.;
223        p    = 2;
224        dw = f.Reverse(p, w);
225
226        // check first derivative in dw
227        double check = 2. * (Value( X[0] ) + 1.);
228        ok &= NearEqual(dw[0], check, eps, eps); 
229
230        // check second derivative in dw
231        check = 2.;
232        ok &= NearEqual(dw[1], check, eps, eps); 
233
234        // Test LdvOp
235        // second order reverse (test exp_if_true case)
236        w[0] = 0.;
237        w[1] = 1.;
238        p    = 2;
239        dw = f.Reverse(p, w);
240
241        // check first derivative in dw
242        check = 2. * (Value( X[0] ) + 1.);
243        ok &= NearEqual(dw[0], check, eps, eps); 
244
245        // check second derivative in dw
246        check = 2.;
247        ok &= NearEqual(dw[1], check, eps, eps); 
248
249        return ok;
250}
251
252} // END empty namespace
253
254bool VecAD(void)
255{       bool ok = true;
256        ok &= VecADTestOne();
257        ok &= VecADTestTwo(); 
258        ok &= SecondOrderReverse();
259        return ok;
260}
Note: See TracBrowser for help on using the repository browser.