source: trunk/test_more/pow.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: 10.1 KB
Line 
1/* $Id: pow.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 just used for validation testing.
15*/
16# include <cppad/cppad.hpp>
17
18namespace { // BEGIN empty namespace
19
20bool PowTestOne(void)
21{       bool ok = true;
22
23        using CppAD::AD;
24        using CppAD::NearEqual;
25
26        // domain space vector
27        size_t n  = 2;
28        double x = 0.5;
29        double y = 2.;
30        CPPAD_TESTVECTOR(AD<double>) XY(n);
31        XY[0]      = x;
32        XY[1]      = y;
33
34        // declare independent variables and start tape recording
35        CppAD::Independent(XY);
36
37        // range space vector
38        size_t m = 3;
39        CPPAD_TESTVECTOR(AD<double>) Z(m);
40        Z[0] = CppAD::pow(XY[0], XY[1]);  // pow(variable, variable)
41        Z[1] = CppAD::pow(XY[0], y);      // pow(variable, parameter)
42        Z[2] = CppAD::pow(x,     XY[1]);  // pow(parameter, variable)
43
44        // create f: XY -> Z and stop tape recording
45        CppAD::ADFun<double> f(XY, Z); 
46
47        // check value
48        double check = std::pow(x, y);
49        size_t i;
50        for(i = 0; i < m; i++)
51                ok &= NearEqual(Z[i] , check,  1e-10 , 1e-10);
52
53        // forward computation of first partial w.r.t. x
54        CPPAD_TESTVECTOR(double) dxy(n);
55        CPPAD_TESTVECTOR(double) dz(m);
56        dxy[0] = 1.;
57        dxy[1] = 0.;
58        dz    = f.Forward(1, dxy);
59        check = y * std::pow(x, y-1.);
60        ok   &= NearEqual(dz[0], check, 1e-10, 1e-10);
61        ok   &= NearEqual(dz[1], check, 1e-10, 1e-10);
62        ok   &= NearEqual(dz[2],    0., 1e-10, 1e-10);
63
64        // forward computation of first partial w.r.t. y
65        dxy[0] = 0.;
66        dxy[1] = 1.;
67        dz    = f.Forward(1, dxy);
68        check = std::log(x) * std::pow(x, y);
69        ok   &= NearEqual(dz[0], check, 1e-10, 1e-10);
70        ok   &= NearEqual(dz[1],    0., 1e-10, 1e-10);
71        ok   &= NearEqual(dz[2], check, 1e-10, 1e-10);
72
73        // reverse computation of derivative of z[0] + z[1] + z[2]
74        CPPAD_TESTVECTOR(double)  w(m);
75        CPPAD_TESTVECTOR(double) dw(n);
76        w[0]  = 1.;
77        w[1]  = 1.;
78        w[2]  = 1.;
79        dw    = f.Reverse(1, w);
80        check = y * std::pow(x, y-1.);
81        ok   &= NearEqual(dw[0], 2. * check, 1e-10, 1e-10);
82        check = std::log(x) * std::pow(x, y);
83        ok   &= NearEqual(dw[1], 2. * check, 1e-10, 1e-10);
84
85        // use a VecAD<Base>::reference object with pow
86        CppAD::VecAD<double> v(2);
87        AD<double> zero(0);
88        AD<double> one(1);
89        v[zero]           = XY[0];
90        v[one]            = XY[1];
91        AD<double> result = CppAD::pow(v[zero], v[one]);
92        ok               &= NearEqual(result, Z[0], 1e-10, 1e-10);
93
94        return ok;
95}
96
97bool PowTestTwo(void)
98{       bool ok = true;
99
100        using CppAD::pow;
101        using CppAD::exp;
102        using namespace CppAD;
103
104
105        // independent variable vector, indices, values, and declaration
106        CPPAD_TESTVECTOR(AD<double>) U(2);
107        size_t s = 0;
108        size_t t = 1;
109        U[s]     = 2.;
110        U[t]     = 3.;
111        Independent(U);
112
113        // dependent variable vector and indices
114        CPPAD_TESTVECTOR(AD<double>) Z(2);
115        size_t x = 0;
116        size_t y = 1;
117
118
119        // dependent variable values
120        AD<double> u = exp(U[s]);        // u = exp(s)
121        Z[x]         = pow(u, U[t]);     // x = exp(s * t)
122        Z[y]         = pow(Z[x], u);     // y = exp( s * t * exp(s) )
123
124        // create f: U -> Z and vectors used for derivative calculations
125        ADFun<double> f(U, Z);
126        CPPAD_TESTVECTOR(double) v( f.Domain() );
127        CPPAD_TESTVECTOR(double) w( f.Range() );
128
129        /*
130        u_s  (s, t) = u
131        u_t  (s, t) = 0
132        y_s  (s, t) = (1 + s) t * u * y
133        y_t  (s, t) = s * u * y
134        y_st (s, t) = ( u + s * u ) * y
135                    + ( t * u + s * t * u ) * s * u * y
136        */
137
138        // check values
139        ok &= NearEqual(Z[x] , exp(2. * 3.),              1e-10 , 1e-10);
140        ok &= NearEqual(Z[y] , exp( 2. * 3. * exp(2.) ),  1e-10 , 1e-10);
141
142        // forward computation of partials w.r.t. s
143        v[s] = 1.;
144        v[t] = 0.;
145        w = f.Forward(1, v);
146        ok &= ( w[x] == U[t] * Z[x] );                   // dx/ds
147        ok &= ( w[y] == (1. + U[s]) * U[t] * u * Z[y] ); // dy/ds
148
149        // forward computation of partials w.r.t. t
150        v[s] = 0.;
151        v[t] = 1.;
152        w = f.Forward(1, v);
153        ok &= ( w[y] == U[s] * u * Z[y] );               // dy/dt
154
155        // forward computation of second Taylor coefficient w.r.t. t
156        v[t] = 1.;
157        w    = f.Forward(1, v);
158        v[t] = 0.;
159        CPPAD_TESTVECTOR(double) f_tt = f.Forward(2, v);
160
161        // forward computation of second Taylor coefficient w.r.t. s
162        v[s] = 1.;
163        w    = f.Forward(1, v);
164        v[s] = 0.;
165        CPPAD_TESTVECTOR(double) f_ss = f.Forward(2, v);
166
167        // second Taylor coefficient w.r.t. direction r = (s,t)
168        v[s] = 1.;
169        v[t] = 1.;
170        w    = f.Forward(1, v);
171        v[s] = 0.;
172        v[t] = 0.;
173        CPPAD_TESTVECTOR(double) f_rr = f.Forward(2, v);
174
175        // check second order partial of y
176        ok &= NearEqual(
177                f_rr[y] - f_ss[y] - f_tt[y], 
178                (1. + U[s]) * u * Z[y] + 
179                        (1. + U[s]) * U[t] * u * U[s] * u * Z[y],
180                1e-10 ,
181                1e-10 
182        ); 
183
184        return ok;
185}
186
187bool PowTestThree(void)
188{       bool ok = true;
189
190        using CppAD::AD;
191        using CppAD::NearEqual;
192
193        // domain space vector
194        size_t n  = 1;
195        CPPAD_TESTVECTOR(AD<double>) x(n);
196        x[0]      = 2.;
197
198        // declare independent variables and start tape recording
199        CppAD::Independent(x);
200
201        // range space vector
202        size_t m = 4;
203        CPPAD_TESTVECTOR(AD<double>) y(m);
204
205        // some special cases
206        y[0] = pow(x[0], 0.);
207        y[1] = pow(0., x[0]);
208        y[2] = pow(x[0], 1.);
209        y[3] = pow(1., x[0]);
210
211        // create f: x -> y and stop tape recording
212        CppAD::ADFun<double> f(x, y); 
213
214        // check function values
215        ok  &= (Value(y[0]) == 1.);
216        ok  &= (Value(y[1]) == 0.);
217        ok  &= (Value(y[2]) == Value(x[0]));
218        ok  &= (Value(y[3]) == 1.);
219
220        // forward computation of first partial w.r.t. x
221        CPPAD_TESTVECTOR(double) dx(n);
222        CPPAD_TESTVECTOR(double) dy(m);
223        dx[0] = 1.;
224        dy    = f.Forward(1, dx);
225        ok   &= (dy[0] == 0.);
226        ok   &= (dy[1] == 0.);
227        ok   &= NearEqual(dy[2], 1., 1e-10, 1e-10);
228        ok   &= (dy[3] == 0.);
229
230        // reverse mode computation of derivative of y[0]+y[1]+y[2]+y[3]
231        CPPAD_TESTVECTOR(double)  w(m);
232        CPPAD_TESTVECTOR(double) dw(n);
233        w[0] = 1.;
234        w[1] = 1.;
235        w[2] = 1.;
236        w[3] = 1.;
237        dw   = f.Reverse(1, w);
238        ok  &= NearEqual(dw[0], 1., 1e-10, 1e-10);
239
240        return ok;     
241}
242
243bool PowTestFour(void)
244{       bool ok = true;
245
246        using CppAD::AD;
247        using CppAD::NearEqual;
248
249        // domain space vector
250        size_t n  = 1;
251        double x0 = -2;
252        CPPAD_TESTVECTOR(AD<double>) x(n);
253        x[0]      = x0;
254
255        // declare independent variables and start tape recording
256        CppAD::Independent(x);
257
258        // range space vector
259        size_t m = 5;
260        CPPAD_TESTVECTOR(AD<double>) y(m);
261
262        // some special cases (skip zero raised to a negative power)
263        y[0] = pow(1., x[0]);
264        size_t i;
265        for(i = 1; i < m; i++) 
266                y[i] = pow(x[0], i-1);   // pow(AD<double>, int)
267
268        // create f: x -> y and stop tape recording
269        CppAD::ADFun<double> f(x, y); 
270
271        ok  &= (Value(y[0]) == 1.);
272        double check;
273        for(i = 1; i < m; i++)
274        {       check = std::pow(x0, double(i-1));
275                ok   &= NearEqual(y[i], check, 1e-10, 1e-10);
276        }
277
278        // forward computation of first partial w.r.t. x
279        CPPAD_TESTVECTOR(double) dx(n);
280        CPPAD_TESTVECTOR(double) dy(m);
281        dx[0] = 1.;
282        dy    = f.Forward(1, dx);
283        ok   &= (dy[0] == 0.);
284        double sum = 0;
285        for(i = 1; i < m; i++)
286        {       if( i == 1 )
287                        check = 0.;
288                else    check = double(i-1) * std::pow(x0, double(i-2));
289                ok   &= NearEqual(dy[i], check, 1e-10, 1e-10);
290                sum  += check;
291        }
292
293        // reverse mode computation of derivative of y[0] + .. y[m-1];
294        CPPAD_TESTVECTOR(double)  w(m);
295        CPPAD_TESTVECTOR(double) dw(n);
296        for(i = 0; i < m; i++)
297                w[i] = 1.;
298        dw   = f.Reverse(1, w);
299        ok  &= NearEqual(dw[0], sum, 1e-10, 1e-10);
300
301        return ok;     
302}
303bool PowTestFive(void)
304{       bool ok = true;
305
306        using CppAD::AD;
307        using CppAD::NearEqual;
308
309        // domain space vector
310        size_t n  = 1;
311        double x0 = -1.;
312        CPPAD_TESTVECTOR(AD<double>) x(n);
313        x[0]      = x0;
314
315        // declare independent variables and start tape recording
316        CppAD::Independent(x);
317
318        // range space vector
319        size_t m = 1;
320        CPPAD_TESTVECTOR(AD<double>) y(m);
321
322        // case of zero raised to a positive integer power
323        double e = 2.;
324        y[0] = pow(x[0], int(e)); // use pow(AD<double>, int)
325
326        // create f: x -> y and stop tape recording
327        CppAD::ADFun<double> f(x, y); 
328
329        // check function value
330        ok  &= (Value(y[0]) == pow(x0, e) );
331
332        // forward computation of first partial w.r.t. x[1]
333        double d1 = e * pow(x0, (e-1));
334        CPPAD_TESTVECTOR(double) dx(n);
335        CPPAD_TESTVECTOR(double) dy(m);
336        dx[0] = 1.;
337        dy    = f.Forward(1, dx);
338        ok   &= NearEqual(dy[0], d1, 1e-10, 1e-10);
339
340        // reverse mode computation of second partials
341        // x.r.t. x[1],x[0]  and x[1], x[1]
342        double d2 = e * (e-1) * pow(x0, (e-2));
343        CPPAD_TESTVECTOR(double)   w(m);
344        CPPAD_TESTVECTOR(double) ddw(2*n);
345        w[0] = 1.;
346        ddw  = f.Reverse(2, w);
347        ok  &= NearEqual(ddw[0], d1, 1e-10, 1e-10);
348        ok  &= NearEqual(ddw[1], d2, 1e-10, 1e-10);
349
350        return ok;     
351}
352bool PowTestSix(void)
353{       bool ok = true;
354
355        using CppAD::AD;
356        using CppAD::NearEqual;
357
358        // domain space vector
359        size_t n  = 1;
360        double x0 = 1.5;
361        CPPAD_TESTVECTOR(AD<double>) x(n);
362        x[0]      = x0;
363
364        // domain space vector
365        CPPAD_TESTVECTOR(AD< AD<double> >) X(n);
366        X[0]      = x[0];
367
368        // declare independent variables and start tape recording
369        CppAD::Independent(X);
370
371        // range space vector
372        size_t m = 1;
373        CPPAD_TESTVECTOR(AD< AD<double> >) Y(m);
374
375        // case of AD< AD<double> > raised to a double power
376        double e = 2.5;
377        Y[0] = pow(X[0], e); 
378
379        // create F: X -> Y and stop tape recording
380        CppAD::ADFun< AD<double> > F(X, Y); 
381
382        // check function value
383        ok  &= (Value( Value(Y[0]) ) == pow(x0, e) );
384
385        // forward computation of first partial w.r.t. x[1]
386        double d1 = e * pow(x0, (e-1));
387        CPPAD_TESTVECTOR(AD<double>) dx(n);
388        CPPAD_TESTVECTOR(AD<double>) dy(m);
389        dx[0] = 1.;
390        dy    = F.Forward(1, dx);
391        ok   &= NearEqual(dy[0], d1, 1e-10, 1e-10);
392
393        // reverse mode computation of second partials
394        // x.r.t. x[1],x[0]  and x[1], x[1]
395        double d2 = e * (e-1) * pow(x0, (e-2));
396        CPPAD_TESTVECTOR(AD<double>)   w(m);
397        CPPAD_TESTVECTOR(AD<double>) ddw(2*n);
398        w[0] = 1.;
399        ddw  = F.Reverse(2, w);
400        ok  &= NearEqual(ddw[0], d1, 1e-10, 1e-10);
401        ok  &= NearEqual(ddw[1], d2, 1e-10, 1e-10);
402
403        return ok;     
404}
405
406} // END empty namespace
407 
408bool Pow(void)
409{       bool ok = true;
410        ok     &= PowTestOne();
411        ok     &= PowTestTwo();
412        ok     &= PowTestThree();
413        ok     &= PowTestFour();
414        ok     &= PowTestFive();
415        ok     &= PowTestSix();
416        return ok;
417}
Note: See TracBrowser for help on using the repository browser.