source: trunk/test_more/reverse.cpp @ 2354

Last change on this file since 2354 was 1771, checked in by bradbell, 9 years ago

Fix warnings about variables shadowing other identifiers (gcc -Wshadow).
rosen_34.cpp, bender_quad.cpp, ode_grea.cpp, ode_gear_control.cpp,
reverse.cpp, k_gt_one.cpp, multiple_solution.cpp, get_started.cpp,
ode_simple.hpp, ode_fast.hpp, sparse_set.hpp, sparse_pack.hpp, det_minor.cpp

Update version number:
AUTHORS, configure, configure.ac, config.h, configure.hpp

Add -Wshadow to compliation command:
build.sh, */test_one.in

makefile.in: fix location of cppad.pc.
whats_new_10.omh: user's view of the changes.

  • Property svn:keywords set to Id
File size: 6.3 KB
Line 
1/* $Id: reverse.cpp 1771 2011-01-01 15:41:51Z 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/*
14Old Reverse example now used just for valiadation testing
15*/
16
17# include <cppad/cppad.hpp>
18namespace { // ----------------------------------------------------------
19
20bool Reverse(void)
21{       bool ok = true;
22
23        using namespace CppAD;
24
25        // independent variable vector
26        CPPAD_TEST_VECTOR< 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_TEST_VECTOR< 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_TEST_VECTOR<double> v( f.Range() );
50        CPPAD_TEST_VECTOR<double> u0( f.Domain() );
51        CPPAD_TEST_VECTOR<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_TEST_VECTOR<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_TEST_VECTOR<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_TEST_VECTOR<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_TEST_VECTOR<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_TEST_VECTOR< 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_TEST_VECTOR< 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} // End empty namespace
207
208# include <vector>
209# include <valarray>
210bool reverse(void)
211{       bool ok = true;
212        ok &= Reverse();
213
214        ok &= reverse_any_cases< CppAD::vector  <double> >();
215        ok &= reverse_any_cases< std::vector    <double> >();
216        ok &= reverse_any_cases< std::valarray  <double> >();
217        return ok;
218}
Note: See TracBrowser for help on using the repository browser.