source: trunk/test_more/pow.cpp @ 992

Last change on this file since 992 was 992, checked in by bradbell, 13 years ago

trunk: Extend pow function to work for double, AD<Base> with any Base.

whats_new_07.omh: user's view of the changes.
pow.cpp: test of use of pow( AD< AD<double> > , double ).
base_complex.hpp: correction of minor typo.
pow.hpp: extend to work with AD<Base> and double.

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