source: trunk/test_more/mul_level.cpp @ 2354

Last change on this file since 2354 was 2237, checked in by bradbell, 8 years ago

abs.sh: Demonstrate bug in multi-level abs function.
build.sh: copy configure.hpp from work directory to source.
mul_level.cpp: add abs test case but ifdef it out (does not work yet).
test_one.sh.in: proper handleing of lib64 (copied from example dir).

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