source: trunk/test_more/reverse.cpp @ 3008

Last change on this file since 3008 was 2859, checked in by bradbell, 7 years ago

merge in changes from branches/atomic; see bin/svn_merge.sh

  • Property svn:keywords set to Id
File size: 13.4 KB
Line 
1/* $Id: reverse.cpp 2859 2013-05-28 06:03:21Z 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/*
14Old Reverse example now used just for valiadation testing
15*/
16
17# include <cppad/cppad.hpp>
18namespace { // ----------------------------------------------------------
19
20bool reverse_one(void)
21{       bool ok = true;
22
23        using namespace CppAD;
24
25        // independent variable vector
26        CPPAD_TESTVECTOR(AD<double>) U(3);
27        U[0] = 0.; U[1] = 1.; U[2] = 2.;
28        Independent(U);
29
30        // compute sum and product of elements in U
31        AD<double> Sum  = 0.;
32        AD<double> Prod = 1.;
33        size_t i;
34        for(i = 0; i < 3; i++)
35        {       Sum  += U[i];
36                Prod *= U[i];
37        }
38
39        // dependent variable vector
40        CPPAD_TESTVECTOR(AD<double>) V(2);
41        V[0] = Sum;
42        V[1] = Prod;
43
44        // V = f(U)
45        ADFun<double> f(U, V);
46
47        // Evaluate ( v[0] * f_0 + v[1] * f_1 )^(1) [ u0 ] ---------------
48        size_t p  = 1;
49        CPPAD_TESTVECTOR(double) v( f.Range() );
50        CPPAD_TESTVECTOR(double) u0( f.Domain() );
51        CPPAD_TESTVECTOR(double) r1( f.Domain() * p );
52
53        v[0]  = 1.; v[1] = -1.;
54        r1    = f.Reverse(p, v);
55
56        // direct evaluation of gradients of components of f
57        CPPAD_TESTVECTOR(double) g0(3), g1(3);
58        u0[0] = Value(U[0]); u0[1] = Value(U[1]); u0[2] = Value(U[2]);
59        g0[0] =          1.; g0[1] =          1.; g0[2] =          1.;
60        g1[0] = u0[1]*u0[2]; g1[1] = u0[0]*u0[2]; g1[2] = u0[0]*u0[1];
61
62        // compare values
63        for(i = 0; i < 3; i++)
64        {       ok &= NearEqual(r1[i] , 
65                        v[0] * g0[i] + v[1] * g1[i], 1e-10, 1e-10);
66        }
67
68        // -------------------------------------------------------------------
69
70        // Define the function z(t, u0, u1) = f( u0 + u1 * t ) and evaluate
71        // the first order Taylor coefficient column vector z(*, u0, u1)
72        p = 1;
73        CPPAD_TESTVECTOR(double) u1( f.Domain() );
74
75        u1[0] = 2.; u1[1] = -1.; u1[2] = 3.;
76        f.Forward(p, u1);
77
78        // Evaluate the derivaties with respect to u0 of the functions
79        // order 0: v[0] *      z_0 (0, u0, u1) + v[1] *      z_1 (0, u0, u1)
80        // order 1: v[0] * d/dt z_0 (0, u0, u1) + v[1] * d/dt z_1 (0, u0, u1)
81        p    = 2;
82        CPPAD_TESTVECTOR(double) r2( f.Domain() * p );
83        v[0] = -.5; v[1] = .5;
84        r2   = f.Reverse(p, v);
85
86        // check derivative of the zero order term
87        for(i = 0; i < 3; i++)
88        {       ok &= NearEqual(r2[p * i + 0] , 
89                        v[0] * g0[i] + v[1] * g1[i], 1e-10, 1e-10);
90        }
91
92        /*
93        The j-th component of the first order term is
94                d/dt z_j(0, u0, u1) = f_j^{(1)} (u0) * u1
95        We use ei to denote the vector with its i-th component one and all
96        the other components zero. The partial derivative of the j-th
97        component of the first order term with respect u0[i] is
98                ei * f_j^{(2)} ( u0 ) * u1
99        */
100
101
102        // direct evaluation of the Hessian f_1^{(2)} (u0)
103        // (the Hessian f_0^{(2)} is identically zero)
104        CPPAD_TESTVECTOR(double) H1(9);
105        H1[0] =    0.; H1[1] = u0[2]; H1[2] = u0[1];
106        H1[3] = u0[2]; H1[4] =    0.; H1[5] = u0[0];
107        H1[6] = u0[1]; H1[7] = u0[0]; H1[8] =    0.;
108
109
110        size_t j;
111        for(i = 0; i < 3; i++)
112        {       double sum = 0.;
113                for(j = 0; j < 3; j++)
114                        sum += H1[i * 3 + j] * u1[j];
115
116                // note term corresponding to v[0] is zero
117                ok &= NearEqual(r2[p * i + 1], v[1] * sum, 1e-10, 1e-10);
118        }
119
120        return ok;
121}
122
123// define the template function reverse_any_cases<Vector> in empty namespace
124template <typename Vector> 
125bool reverse_any_cases(void)
126{       bool ok = true;
127        using CppAD::AD;
128        using CppAD::NearEqual;
129
130        // domain space vector
131        size_t n = 3;
132        CPPAD_TESTVECTOR(AD<double>) X(n);
133        X[0] = 0.; 
134        X[1] = 1.;
135        X[2] = 2.;
136
137        // declare independent variables and start recording
138        CppAD::Independent(X);
139
140        // range space vector
141        size_t m = 1;
142        CPPAD_TESTVECTOR(AD<double>) Y(m);
143        Y[0] = X[0] * X[1] * X[2];
144
145        // create f : X -> Y and stop recording
146        CppAD::ADFun<double> f(X, Y);
147
148        // define W(t, u) = (u_0 + dx_0*t)*(u_1 + dx_1*t)*(u_2 + dx_2*t)
149        // use zero order forward to evaluate W0(u) = W(0, u)
150        Vector u(n), W0(m);
151        u[0]    = 2.;
152        u[1]    = 3.;
153        u[2]    = 4.;
154        W0      = f.Forward(0, u);
155        double check;
156        check   =  u[0]*u[1]*u[2];
157        ok     &= NearEqual(W0[0] , check, 1e-10, 1e-10);
158
159        // define W_t(t, u) = partial W(t, u) w.r.t t
160        // W_t(t, u)  = (u_0 + dx_0*t)*(u_1 + dx_1*t)*dx_2
161        //            + (u_0 + dx_0*t)*(u_2 + dx_2*t)*dx_1
162        //            + (u_1 + dx_1*t)*(u_2 + dx_2*t)*dx_0
163        // use first order forward mode to evaluate W1(u) = W_t(0, u)
164        Vector dx(n), W1(m);
165        dx[0] = .2;
166        dx[1] = .3;
167        dx[2] = .4;
168        W1    = f.Forward(1, dx);
169        check =  u[0]*u[1]*dx[2] + u[0]*u[2]*dx[1] + u[1]*u[2]*dx[0];
170        ok   &= NearEqual(W1[0], check, 1e-10, 1e-10);
171
172        // define W_tt (t, u) = partial W_t(t, u) w.r.t t
173        // W_tt(t, u) = 2*(u_0 + dx_0*t)*dx_1*dx_2
174        //            + 2*(u_1 + dx_1*t)*dx_0*dx_2
175        //            + 2*(u_3 + dx_3*t)*dx_0*dx_1
176        // use second order forward to evaluate W2(u) = 1/2 * W_tt(0, u)
177        Vector ddx(n), W2(m);
178        ddx[0] = ddx[1] = ddx[2] = 0.;
179        W2     = f.Forward(2, ddx);
180        check  =  u[0]*dx[1]*dx[2] + u[1]*dx[0]*dx[2] + u[2]*dx[0]*dx[1];
181        ok    &= NearEqual(W2[0], check, 1e-10, 1e-10);
182
183        // use third order reverse mode to evaluate derivatives
184        size_t p = 3;
185        Vector w(m), dw(n * p);
186        w[0]   = 1.;
187        dw     = f.Reverse(p, w);
188
189        // check derivative of W0(u) w.r.t. u
190        ok    &= NearEqual(dw[0*p+0], u[1]*u[2], 1e-10, 1e-10);
191        ok    &= NearEqual(dw[1*p+0], u[0]*u[2], 1e-10, 1e-10);
192        ok    &= NearEqual(dw[2*p+0], u[0]*u[1], 1e-10, 1e-10);
193
194        // check derivative of W1(u) w.r.t. u
195        ok    &= NearEqual(dw[0*p+1], u[1]*dx[2] + u[2]*dx[1], 1e-10, 1e-10);
196        ok    &= NearEqual(dw[1*p+1], u[0]*dx[2] + u[2]*dx[0], 1e-10, 1e-10);
197        ok    &= NearEqual(dw[2*p+1], u[0]*dx[1] + u[1]*dx[0], 1e-10, 1e-10);
198
199        // check derivative of W2(u) w.r.t u
200        ok    &= NearEqual(dw[0*p+2], dx[1]*dx[2], 1e-10, 1e-10);
201        ok    &= NearEqual(dw[1*p+2], dx[0]*dx[2], 1e-10, 1e-10);
202        ok    &= NearEqual(dw[2*p+2], dx[0]*dx[1], 1e-10, 1e-10);
203
204        return ok;
205}
206/*
207$comment reverse_any.cpp$$
208$spell
209        Taylor
210$$
211
212$section Reverse Mode General Case: Example and Test$$
213
214$index general, reverse example$$
215$index reverse, general example$$
216$index example, general reverse$$
217$index test, general reverse$$
218
219$index composition, example$$
220$index example, composition$$
221$index test, composition$$
222
223$head Purpose$$
224Break a derivative computation into pieces and only store values at the
225interface of the pieces.
226In actual applications, there may be many functions, but
227for this example there are only two.
228The functions
229$latex F : \B{R}^2 \rightarrow \B{R}^2$$
230and
231$latex G : \B{R}^2 \rightarrow \B{R}^2$$
232defined by
233$latex \[
234        F(x) = \left( \begin{array}{c} x_0 x_1   \\ x_1 - x_0 \end{array} \right)
235        \; , \;
236        G(y) = \left( \begin{array}{c} y_0 - y_1 \\ y_1  y_0   \end{array} \right)
237\] $$
238Another difference is that in actual applications,
239the memory corresponding to function objects not currently being used
240is sometimes returned to the system (see $cref checkpoint.cpp$$).
241
242$head Processing Steps$$
243We apply reverse mode to compute the derivative of
244$latex H : \B{R}^2 \rightarrow \B{R}$$
245is defined by
246$latex \[
247\begin{array}{rcl}
248        H(x)
249        & = & G_0 [ F(x) ] + G_1 [ F(x)  ]
250        \\
251        & = & x_0 x_1 - ( x_1 - x_0 ) + x_0 x_1 ( x_1 - x_0 )
252        \\
253        & = & x_0 x_1 ( 1 - x_0 + x_1 ) - x_1 + x_0
254\end{array}
255\] $$
256Given the zero and first order Taylor coefficients
257$latex x^{(0)} $$ and $latex x^{(1)}$$,
258we use $latex X(t)$$, $latex Y(t)$$ and $latex Z(t)$$
259for the corresponding functions; i.e.,
260$latex \[
261\begin{array}{rcl}
262        X(t) & = & x^{(0)} + x^{(1)} t
263        \\
264        Y(t) & = & F[X(t)] = y^{(0)} + y^{(1)} t  + O(t^2)
265        \\
266        Z(t) & = & G \{ F [ X(t) ] \} = z^{(0)} + z^{(1)} t  + O(t^2)
267        \\
268        h^{(0)} & = & z^{(0)}_0 + z^{(0)}_1
269        \\
270        h^{(1)} & = & z^{(1)}_0 + z^{(1)}_1
271\end{array}
272\] $$
273Here are the processing steps:
274$list number$$
275Use forward mode on $latex F(x)$$ to compute
276$latex y^{(0)}$$ and $latex y^{(1)}$$
277$lnext
278Use forward mode on $latex G(y)$$ to compute
279$latex z^{(0)}$$ and $latex z^{(1)}$$
280$lnext
281Use reverse mode on $latex G(y)$$ to compute the derivative of
282$latex h^{(k)}$$ with respect to
283$latex y^{(0)}$$ and $latex y^{(1)}$$.
284$lnext
285Use reverse mode on $latex F(x)$$ to compute the derivative of
286$latex h^{(k)}$$ with respect to
287$latex x^{(0)}$$ and $latex x^{(1)}$$.
288$lend
289This uses the following relations for $latex k = 0 , 1$$:
290$latex \[
291\begin{array}{rcl}
292        \partial_{x(0)} h^{(k)} [ x^{(0)} , x^{(1)} ]
293        & = &
294        \partial_{y(0)} h^{(k)} [ y^{(0)} , y^{(1)} ]
295        \partial_{x(0)} y^{(0)} [ x^{(0)} , x^{(1)} ]
296        \\
297        & + &
298        \partial_{y(1)} h^{(k)} [ y^{(0)} , y^{(1)} ]
299        \partial_{x(0)} y^{(1)} [ x^{(0)} , x^{(1)} ]
300        \\
301        \partial_{x(1)} h^{(k)} [ x^{(0)} , x^{(1)} ]
302        & = &
303        \partial_{y(0)} h^{(k)} [ y^{(0)} , y^{(1)} ]
304        \partial_{x(1)} y^{(0)} [ x^{(0)} , x^{(1)} ]
305        \\
306        & + &
307        \partial_{y(1)} h^{(k)} [ y^{(0)} , y^{(1)} ]
308        \partial_{x(1)} y^{(1)} [ x^{(0)} , x^{(1)} ]
309\end{array}
310\] $$
311where $latex \partial_{x(0)}$$ denotes the partial with respect
312to $latex x^{(0)}$$.
313
314$code
315$comment%example/reverse_any.cpp%0%// BEGIN C++%// END C++%1%$$
316$$
317$end
318*/
319template <class Vector>
320Vector F_reverse_mul(const Vector& x)
321{       Vector y(2);
322        y[0] = x[0] * x[1];
323        y[1] = x[1] - x[0];
324        return y;
325}
326template <class Vector>
327Vector G_reverse_mul(const Vector& y)
328{       Vector z(2);
329        z[0] = y[0] - y[1];
330        z[1] = y[1] * y[0];
331        return z;
332}
333bool reverse_mul(void)
334{
335        bool ok = true;
336     double eps = 10. * CppAD::numeric_limits<double>::epsilon();
337
338        using CppAD::AD;
339        using CppAD::NearEqual;
340        CppAD::ADFun<double> f, g;
341
342        // Record the function F(x)
343        size_t n    = 2;
344        CPPAD_TESTVECTOR(AD<double>) X(n), Y(n);
345        X[0] = X[1] = 0.;
346        CppAD::Independent(X);
347        Y = F_reverse_mul(X);
348        f.Dependent(X, Y);
349
350        // Record the function G(x)
351        CPPAD_TESTVECTOR(AD<double>) Z(n);
352        Y[0] = Y[1] = 0.;
353        CppAD::Independent(Y);
354        Z = G_reverse_mul(Y);
355        g.Dependent(Y, Z);
356
357        // argument and function values
358        CPPAD_TESTVECTOR(double) x0(n), y0(n), z0(n);
359        x0[0] = 1.;
360        x0[1] = 2.;
361        y0    = f.Forward(0, x0);
362        z0    = g.Forward(0, y0);
363
364        // check function value
365        double check = x0[0] * x0[1] * (1. - x0[0] + x0[1]) - x0[1] + x0[0];
366        double h0    = z0[0] + z0[1];
367        ok          &= NearEqual(h0, check, eps, eps);
368
369        // first order Taylor coefficients
370        CPPAD_TESTVECTOR(double) x1(n), y1(n), z1(n);
371        x1[0] = 3.;
372        x1[1] = 4.;
373        y1    = f.Forward(1, x1);
374        z1    = g.Forward(1, y1);
375
376        // check first order Taylor coefficients
377        check     = x0[0] * x0[1] * (- x1[0] + x1[1]) - x1[1] + x1[0];
378        check    += x1[0] * x0[1] * (1. - x0[0] + x0[1]);
379        check    += x0[0] * x1[1] * (1. - x0[0] + x0[1]);
380        double h1 = z1[0] + z1[1];
381        ok       &= NearEqual(h1, check, eps, eps);
382
383        // ----------------------------------------------------------------
384        // dw^0 (y) = \partial_y^0 h^0 (y)
385        // dw^1 (y) = \partial_y^1 h^0 (y)
386        size_t p = 2;
387        CPPAD_TESTVECTOR(double) w(n*p), dw(n*p);
388        w[0*p+0] = 1.; // coefficient for z^0_0
389        w[1*p+0] = 1.; // coefficient for z^0_1
390        w[0*p+1] = 0.; // coefficient for z^1_0
391        w[1*p+1] = 0.; // coefficient for z^1_1
392        dw       = g.Reverse(p, w);
393
394        // dv^0 = dw^0 * \partial_x^0 y^0 (x) + dw^1 * \partial_x^0 y^1 (x) 
395        // dv^1 = dw^0 * \partial_x^1 y^0 (x) + dw^1 * \partial_x^1 y^1 (x) 
396        CPPAD_TESTVECTOR(double) dv(n*p);
397        dv   = f.Reverse(p, dw); 
398
399        // check partial of h^0 w.r.t x^0_0
400        check  = x0[1] * (1. - x0[0] + x0[1]) + 1.;
401        check -= x0[0] * x0[1];
402        ok    &= NearEqual(dv[0*p+0], check, eps, eps);
403
404        // check partial of h^0 w.r.t x^0_1
405        check  = x0[0] * (1. - x0[0] + x0[1]) - 1.;
406        check += x0[0] * x0[1];
407        ok    &= NearEqual(dv[1*p+0], check, eps, eps);
408
409        // check partial of h^0 w.r.t x^1_0 and x^1_1
410        check  = 0.;
411        ok    &= NearEqual(dv[0*p+1], check, eps, eps);
412        ok    &= NearEqual(dv[1*p+1], check, eps, eps);
413
414        // ----------------------------------------------------------------
415        // dw^0 (y) = \partial_y^0 h^1 (y)
416        // dw^1 (y) = \partial_y^1 h^1 (y)
417        w[0*p+0] = 0.; // coefficient for z^0_0
418        w[1*p+0] = 0.; // coefficient for z^0_1
419        w[0*p+1] = 1.; // coefficient for z^1_0
420        w[1*p+1] = 1.; // coefficient for z^1_1
421        dw       = g.Reverse(p, w);
422
423        // dv^0 = dw^0 * \partial_x^0 y^0 (x) + dw^1 * \partial_x^0 y^1 (x) 
424        // dv^1 = dw^0 * \partial_x^1 y^0 (x) + dw^1 * \partial_x^1 y^1 (x) 
425        dv   = f.Reverse(p, dw); 
426
427        // check partial of h^1 w.r.t x^0_0
428        check  = x0[1] * (- x1[0] + x1[1]);
429        check -= x1[0] * x0[1];
430        check += x1[1] * (1. - x0[0] + x0[1]) - x0[0] * x1[1];
431        ok    &= NearEqual(dv[0*p+0], check, eps, eps);
432
433        // check partial of h^1 w.r.t x^0_1
434        check  = x0[0] * (- x1[0] + x1[1]);
435        check += x1[0] * (1. - x0[0] + x0[1]) + x1[0] * x0[1]; 
436        check += x0[0] * x1[1];
437        ok    &= NearEqual(dv[1*p+0], check, eps, eps);
438
439        // check partial of h^1 w.r.t x^1_0
440        // (by reverse mode identity is equal to partial h^0 w.r.t. x^0_0)
441        check  = 1. - x0[0] * x0[1];
442        check += x0[1] * (1. - x0[0] + x0[1]);
443        ok    &= NearEqual(dv[0*p+1], check, eps, eps);
444
445        // check partial of h^1 w.r.t x^1_1
446        // (by reverse mode identity is equal to partial h^0 w.r.t. x^0_1)
447        check  = x0[0] * x0[1] - 1.;
448        check += x0[0] * (1. - x0[0] + x0[1]);
449        ok    &= NearEqual(dv[1*p+1], check, eps, eps);
450
451        return ok;
452}
453// ----------------------------------------------------------------------------
454} // End empty namespace
455
456# include <vector>
457# include <valarray>
458bool reverse(void)
459{       bool ok = true;
460        ok &= reverse_one();
461        ok &= reverse_mul();
462
463        ok &= reverse_any_cases< CppAD::vector  <double> >();
464        ok &= reverse_any_cases< std::vector    <double> >();
465        ok &= reverse_any_cases< std::valarray  <double> >();
466        return ok;
467}
Note: See TracBrowser for help on using the repository browser.