source: trunk/test_more/mul_level.cpp @ 3071

Last change on this file since 3071 was 3071, checked in by bradbell, 6 years ago

test_one.sh.in: update choices when compiling eigen versus not.
whats_new_13.omh: improve discussion of eigen.
mul_level.cpp: test conversion from double to AD<AD<double>>.
test_one.sh.in: update choices when compiling eigen versus not.

  • Property svn:keywords set to Id
File size: 9.5 KB
Line 
1/* $Id: mul_level.cpp 3071 2014-01-01 04:02:27Z 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# ifdef CPPAD_ADOLC_TEST
15# include <adolc/adouble.h>
16# include <adolc/taping.h>
17# include <adolc/interfaces.h>
18// adouble definitions not in Adolc distribution and
19// required in order to use CppAD::AD<adouble>
20# include <cppad/example/base_adolc.hpp>
21# endif
22
23# include <cppad/cppad.hpp>
24# include <limits>
25
26namespace { // BEGIN empty namespace
27
28bool One(void) 
29{       bool ok = true;                          // initialize test result
30        using CppAD::NearEqual;
31        double eps = 10. * std::numeric_limits<double>::epsilon();
32
33
34        typedef CppAD::AD<double>   ADdouble;    // for one level of taping
35        typedef CppAD::AD<ADdouble> ADDdouble;   // for two levels of taping
36        size_t n = 2;                            // dimension for example
37
38        // value of the independent variables
39        CPPAD_TESTVECTOR(ADDdouble) aa_x(n);
40        aa_x[0] = 1.; aa_x[1] = 3.; // test conversion double to AD< AD<double> >
41        aa_x[0] = 2. * aa_x[0];     // test double * AD< AD<double> >
42        CppAD::Independent(aa_x);
43
44        // compute the function f(x) = 2 * x[0] * x[1]
45        CPPAD_TESTVECTOR(ADDdouble) aa_f(1);
46        aa_f[0] = 2. * aa_x[0] * aa_x[1];
47        CppAD::ADFun<ADdouble> F(aa_x, aa_f);
48
49        // value of the independent variables
50        CPPAD_TESTVECTOR(ADdouble)   a_x(n);
51        a_x[0] = 2.; a_x[1] = 3.;
52        Independent(a_x);
53
54        // re-evaluate f(2, 3) (must get deepedence on a_x).
55        size_t p = 0;
56        CPPAD_TESTVECTOR(ADdouble) a_fp(1);
57        a_fp    = F.Forward(p, a_x);
58        ok     &= NearEqual(a_fp[0], 2. * a_x[0] * a_x[1], eps, eps);
59
60        // compute the function g(x) = 2 * partial_x[0] f(x) = 4 * x[1]
61        p = 1;
62        CPPAD_TESTVECTOR(ADdouble) a_dx(n), a_g(1);
63        a_dx[0] = 1.; a_dx[1] = 0.;
64        a_fp    = F.Forward(p, a_dx);
65        a_g[0]  = 2. * a_fp[0];
66        CppAD::ADFun<double> G(a_x, a_g);
67
68        // compute partial_x[1] g(x)
69        CPPAD_TESTVECTOR(double)  xp(n), gp(1);
70        p = 0;
71        xp[0] = 4.; xp[1] = 5.;
72        gp    = G.Forward(p, xp);
73        ok   &= NearEqual(gp[0], 4. * xp[1], eps, eps);
74
75        p = 1;
76        xp[0] = 0.; xp[1] = 1.;
77        gp    = G.Forward(p, xp);
78        ok   &= NearEqual(gp[0], 4., eps, eps);
79
80        return ok;
81}
82
83// f(x) = |x|^2 = .5 * ( x[0]^2 + ... + x[n-1]^2 )
84template <class Type>
85Type f_Two(CPPAD_TESTVECTOR(Type) &x)
86{       Type sum;
87
88        // check assignment of AD< AD<double> > = double
89        sum  = .5;
90        sum += .5;
91
92        size_t i = x.size();
93        while(i--)
94                sum += x[i] * x[i];
95
96        // check computed assignment AD< AD<double> > -= int
97        sum -= 1; 
98
99        // check double * AD< AD<double> >
100        return .5 * sum;
101} 
102
103bool Two(void) 
104{       bool ok = true;                          // initialize test result
105
106        typedef CppAD::AD<double>   ADdouble;    // for one level of taping
107        typedef CppAD::AD<ADdouble> ADDdouble;   // for two levels of taping
108        size_t n = 5;                            // dimension for example
109        size_t j;                                // a temporary index variable
110
111        CPPAD_TESTVECTOR(double)       x(n);
112        CPPAD_TESTVECTOR(ADdouble)   a_x(n);
113        CPPAD_TESTVECTOR(ADDdouble) aa_x(n);
114
115        // value of the independent variables
116        for(j = 0; j < n; j++)
117                a_x[j] = x[j] = double(j); // x[j] = j
118        Independent(a_x);                  // a_x is indedendent for ADdouble
119        for(j = 0; j < n; j++)
120                aa_x[j] = a_x[j];          // track how aa_x depends on a_x
121        CppAD::Independent(aa_x);          // aa_x is independent for ADDdouble
122
123        // compute function
124        CPPAD_TESTVECTOR(ADDdouble) aa_f(1);    // scalar valued function
125        aa_f[0] = f_Two(aa_x);                   // has only one component
126
127        // declare inner function (corresponding to ADDdouble calculation)
128        CppAD::ADFun<ADdouble> a_F(aa_x, aa_f);
129
130        // compute f'(x)
131        size_t p = 1;                        // order of derivative of a_F
132        CPPAD_TESTVECTOR(ADdouble) a_w(1);  // weight vector for a_F
133        CPPAD_TESTVECTOR(ADdouble) a_df(n); // value of derivative
134        a_w[0] = 1;                          // weighted function same as a_F
135        a_df   = a_F.Reverse(p, a_w);        // gradient of f
136
137        // declare outter function (corresponding to ADdouble calculation)
138        CppAD::ADFun<double> df(a_x, a_df);
139
140        // compute the d/dx of f'(x) * v = f''(x) * v
141        CPPAD_TESTVECTOR(double) v(n);
142        CPPAD_TESTVECTOR(double) ddf_v(n);
143        for(j = 0; j < n; j++)
144                v[j] = double(n - j);
145        ddf_v = df.Reverse(p, v);
146
147        // f(x)       = .5 * ( x[0]^2 + x[1]^2 + ... + x[n-1]^2 )
148        // f'(x)      = (x[0], x[1], ... , x[n-1])
149        // f''(x) * v = ( v[0], v[1],  ... , x[n-1] )
150        for(j = 0; j < n; j++)
151                ok &= CppAD::NearEqual(ddf_v[j], v[j], 1e-10, 1e-10);
152
153        return ok;
154}
155
156# ifdef CPPAD_ADOLC_TEST
157
158bool adolc(void) 
159{       bool ok = true;                   // initialize test result
160
161        typedef adouble      ADdouble;         // for first level of taping
162        typedef CppAD::AD<ADdouble> ADDdouble; // for second level of taping
163        size_t n = 5;                          // number independent variables
164
165        CPPAD_TESTVECTOR(double)       x(n);
166        CPPAD_TESTVECTOR(ADdouble)   a_x(n);
167        CPPAD_TESTVECTOR(ADDdouble) aa_x(n);
168
169        // value of the independent variables
170        int tag = 0;                         // Adolc setup
171        int keep = 1;
172        trace_on(tag, keep);
173        size_t j;
174        for(j = 0; j < n; j++)
175        {       x[j] = double(j);           // x[j] = j
176                a_x[j] <<= x[j];            // a_x is independent for ADdouble
177        }
178        for(j = 0; j < n; j++)
179                aa_x[j] = a_x[j];          // track how aa_x depends on a_x
180        CppAD::Independent(aa_x);          // aa_x is independent for ADDdouble
181
182        // compute function
183        CPPAD_TESTVECTOR(ADDdouble) aa_f(1);    // scalar valued function
184        aa_f[0] = f_Two(aa_x);                   // has only one component
185
186        // declare inner function (corresponding to ADDdouble calculation)
187        CppAD::ADFun<ADdouble> a_F(aa_x, aa_f);
188
189        // compute f'(x)
190        size_t p = 1;                        // order of derivative of a_F
191        CPPAD_TESTVECTOR(ADdouble) a_w(1);  // weight vector for a_F
192        CPPAD_TESTVECTOR(ADdouble) a_df(n); // value of derivative
193        a_w[0] = 1;                          // weighted function same as a_F
194        a_df   = a_F.Reverse(p, a_w);        // gradient of f
195
196        // declare outter function
197        // (corresponding to the tape of adouble operations)
198        double df_j;
199        for(j = 0; j < n; j++)
200                a_df[j] >>= df_j;
201        trace_off();
202
203        // compute the d/dx of f'(x) * v = f''(x) * v
204        size_t m      = n;                     // # dependent in f'(x)
205        double *v = CPPAD_NULL, *ddf_v = CPPAD_NULL;
206        v     = CPPAD_TRACK_NEW_VEC(m, v);     // track v = new double[m]
207        ddf_v = CPPAD_TRACK_NEW_VEC(n, ddf_v); // track ddf_v = new double[n]
208        for(j = 0; j < n; j++)
209                v[j] = double(n - j);
210        fos_reverse(tag, int(m), int(n), v, ddf_v);
211
212        // f(x)       = .5 * ( x[0]^2 + x[1]^2 + ... + x[n-1]^2 )
213        // f'(x)      = (x[0], x[1], ... , x[n-1])
214        // f''(x) * v = ( v[0], v[1],  ... , x[n-1] )
215        for(j = 0; j < n; j++)
216                ok &= CppAD::NearEqual(ddf_v[j], v[j], 1e-10, 1e-10);
217
218        CPPAD_TRACK_DEL_VEC(v);                 // check usage of delete
219        CPPAD_TRACK_DEL_VEC(ddf_v);
220        return ok;
221}
222
223# endif // CPPAD_ADOLC_TEST
224
225bool std_math(void)
226{       bool ok = true;
227        using CppAD::AD;
228        using CppAD::Independent;
229        using CppAD::ADFun;
230        double eps = std::numeric_limits<double>::epsilon();
231
232
233        typedef AD<double>      ADdouble; // for first level of taping
234        typedef AD<ADdouble>   ADDdouble; // for second level of taping
235        size_t n = 1;         // number independent variables
236        size_t m = 1;         // number dependent and independent variables
237
238        CPPAD_TESTVECTOR(double)       x(n),   y(m);
239        CPPAD_TESTVECTOR(ADdouble)    ax(n),  ay(m);
240        CPPAD_TESTVECTOR(ADDdouble)  aax(n), aay(m);
241
242        // create af(x) = tanh(x)
243        aax[0] = 1.;
244        Independent( aax );
245        aay[0] = tanh(aax[0]);
246        ADFun<ADdouble> af(aax, aay);
247
248        // create g(x) = af(x)
249        ax[0] = 1.;
250        Independent( ax );
251        ay = af.Forward(0, ax);
252        ADFun<double> g(ax, ay);
253
254        // evaluate h(x) = g(x)
255        x[0] = 1.;
256        y = g.Forward(0, x);
257
258        // check result
259        double check  = tanh(x[0]);
260        ok &= CppAD::NearEqual(y[0], check, eps, eps);
261
262        return ok;
263}
264
265// supress this test until ../bug/abs.sh is fixed
266# ifdef CPPAD_NOT_DEFINED
267bool abs(void)
268{       bool ok = true;
269        using CppAD::AD;
270        using CppAD::Independent;
271        using CppAD::ADFun;
272        double eps = std::numeric_limits<double>::epsilon();
273
274
275        typedef AD<double>      ADdouble; // for first level of taping
276        typedef AD<ADdouble>   ADDdouble; // for second level of taping
277        size_t n = 1;         // number independent variables
278        size_t m = 1;         // number dependent and independent variables
279
280        CPPAD_TESTVECTOR(double)       x(n),   y(m);
281        CPPAD_TESTVECTOR(ADdouble)    ax(n),  ay(m);
282        CPPAD_TESTVECTOR(ADDdouble)  aax(n), aay(m);
283
284        // create af(x) = abs(x)
285        aax[0] = 1.;
286        Independent( aax );
287        aay[0] = abs(aax[0]);
288        ADFun<ADdouble> af(aax, aay);
289
290        // create g(x) = af'(x)
291        ax[0] = 1.;
292        Independent( ax );
293        ay = af.Jacobian(ax);
294        ADFun<double> g(ax, ay);
295
296        // evaluate g(x) at same x as recording
297        x[0] = 1.;
298        y = g.Forward(0, x);
299
300        // check result
301        double check  = 1.;
302        ok &= CppAD::NearEqual(y[0], check, eps, eps);
303
304        // evaluate g(x) at different x from recording
305        // (but abs is an atomic operation so derivative should work)
306        x[0] = -1.;
307        y = g.Forward(0, x);
308
309        // check result
310        check  = -1.;
311        ok &= CppAD::NearEqual(y[0], check, eps, eps);
312
313        return ok;
314}
315# endif
316
317} // END empty namespace
318
319bool mul_level(void)
320{       bool ok = true;
321        ok     &= One();
322        ok     &= Two();
323# ifdef CPPAD_ADOLC_TEST
324        ok     &= adolc();
325# endif
326        ok     &= std_math();
327# ifdef CPPAD_NOT_DEFINED
328        ok     &= abs();
329# endif
330
331        return ok;
332}
Note: See TracBrowser for help on using the repository browser.