# source:trunk/test_more/forward.cpp@3008

Last change on this file since 3008 was 2570, checked in by bradbell, 7 years ago
2. Fix error in the NDEBUG case.
• Property svn:keywords set to `Id`
File size: 6.2 KB
Line
1/* \$Id: forward.cpp 2570 2012-11-14 18:25:24Z bradbell \$ */
2/* --------------------------------------------------------------------------
4
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.
11-------------------------------------------------------------------------- */
12
13/*
14Two old Forward example now used just for valiadation testing
15*/
16
18
19namespace { // Begin empty namespace
20
21template <typename VectorDouble> // vector class, elements of type double
22bool ForwardCases(void)
23{       bool ok = true;
24
26
27        // independent variable vector
29        X[0] = 0.;
30        X[1] = 1.;
31        Independent(X);
32
33        // compute product of elements in X
35        Y[0] = X[0] * X[0] * X[1];
36
37        // create function object F : X -> Y
39
40        // use zero order to evaluate F[ (3, 4) ]
41        VectorDouble x0( F.Domain() );
42        VectorDouble y0( F.Range() );
43        x0[0]    = 3.;
44        x0[1]    = 4.;
45        y0       = F.Forward(0, x0);
46        ok      &= NearEqual(y0[0] , x0[0]*x0[0]*x0[1], 1e-10, 1e-10);
47
48        // evaluate derivative of F in X[0] direction
49        VectorDouble x1( F.Domain() );
50        VectorDouble y1( F.Range() );
51        x1[0]    = 1.;
52        x1[1]    = 0.;
53        y1       = F.Forward(1, x1);
54        ok      &= NearEqual(y1[0] , 2.*x0[0]*x0[1], 1e-10, 1e-10);
55
56        // evaluate second derivative of F in X[0] direction
57        VectorDouble x2( F.Domain() );
58        VectorDouble y2( F.Range() );
59        x2[0]       = 0.;
60        x2[1]       = 0.;
61        y2          = F.Forward(2, x2);
62        double F_00 = 2. * y2[0];
63        ok         &= NearEqual(F_00, 2.*x0[1], 1e-10, 1e-10);
64
65        // evalute derivative of F in X[1] direction
66        x1[0]    = 0.;
67        x1[1]    = 1.;
68        y1       = F.Forward(1, x1);
69        ok      &= NearEqual(y1[0] , x0[0]*x0[0], 1e-10, 1e-10);
70
71        // evaluate second derivative of F in X[1] direction
72        y2          = F.Forward(2, x2);
73        double F_11 = 2. * y2[0];
74        ok         &= NearEqual(F_11, 0., 1e-10, 1e-10);
75
76        // evalute derivative of F in X[0] + X[1] direction
77        x1[0]    = 1.;
78        x1[1]    = 1.;
79        y1       = F.Forward(1, x1);
80        ok      &= NearEqual(y1[0], 2.*x0[0]*x0[1] + x0[0]*x0[0], 1e-10, 1e-10);
81
82        // use second derivative of F in X[0] direction to
83        // compute second partial of F w.r.t X[1] w.r.t X[2]
84        y2          = F.Forward(2, x2);
85        double F_01 = y2[0] - F_00 / 2. - F_11 / 2.;
86        ok         &= NearEqual(F_01 , 2.*x0[0], 1e-10, 1e-10);
87
88        return ok;
89}
90
91bool ForwardOlder(void)
92{       bool ok = true;
93
95
96        // independent variable vector
98        U[0] = 0.; U[1] = 1.; U[2] = 2.;
99        Independent(U);
100
101        // compute sum and product of elements in U
104        size_t i;
105        for(i = 0; i < 3; i++)
106        {       sum  += U[i];
107                prod *= U[i];
108        }
109
110        // dependent variable vector
112        V[0] = sum;
113        V[1] = prod;
114
115        // V = f(U)
117
118        // use ADFun object to evaluate f[ (1, 2, 3)^T ] -----------------
121        size_t p;
122        p     = 0;
123        u0[0] = 1.; u0[1] = 2.; u0[2] = 3.;
124        v0    = f.Forward(p, u0);
125
126        // direct evaluation of f[ u0 ]
128        f0[0] = u0[0] + u0[1] + u0[2];
129        f0[1] = u0[0] * u0[1] * u0[2];
130
131        // compare values
132        ok &= NearEqual(v0[0] , f0[0], 1e-10, 1e-10);
133        ok &= NearEqual(v0[1] , f0[1], 1e-10, 1e-10);
134
135        // use ADFun object to evaluate f^(1) [ u0 ] * u1 -----------------
138        p     = 1;
139        u1[0] = 1.; u1[1] = 1.; u1[2] = 1.;
140        v1    = f.Forward(p, u1);
141
142        // direct evaluation of gradients of components of f
144        g0[0] =          1.; g0[1] =          1.; g0[2] =          1.;
145        g1[0] = u0[1]*u0[2]; g1[1] = u0[0]*u0[2]; g1[2] = u0[0]*u0[1];
146
147        // compare values
148        ok &= NearEqual(v1[0] ,
149                g0[0]*u1[0] + g0[1]*u1[1] + g0[2]*u1[2] , 1e-10, 1e-10);
150        ok &= NearEqual(v1[1] ,
151                g1[0]*u1[0] + g1[1]*u1[1] + g1[2]*u1[2] , 1e-10, 1e-10);
152
153        // use ADFun object to evaluate ------------------------------------
154        // (1/2) * { f^(1)[ u0 ] * u2 + u1^T * f^(2)[ u0 ] * u1 }
157        p     = 2;
158        u2[0] = .5; u2[1] = .4; u2[2] = .3;
159        v2    = f.Forward(p, u2);
160
161        // direct evaluation of Hessian of second components of f
162        // (the Hessian of the first component is zero)
164        H1[0] =    0.; H1[1] = u0[2]; H1[2] = u0[1];
165        H1[3] = u0[2]; H1[4] =    0.; H1[5] = u0[0];
166        H1[6] = u0[1]; H1[7] = u0[0]; H1[8] =    0.;
167
168        // compare values
169        ok &= NearEqual(v2[0] ,
170                g0[0]*u2[0] + g0[1]*u2[1] + g0[2]*u2[2] , 1e-10, 1e-10);
171
172        size_t j;
173        double v2_1 = 0.;
174        for(i = 0; i < 3; i++)
175        {       v2_1 += g1[i] * u2[i];
176                for(j = 0; j < 3; j++)
177                        v2_1 += .5 * u1[i] * H1[i * 3 + j] * u1[j];
178        }
179        ok &= NearEqual(v2[1], v2_1, 1e-10, 1e-10);
180
181
182        return ok;
183}
184
185void my_error_handler(
186        bool known           ,
187        int  line            ,
188        const char *file     ,
189        const char *exp      ,
190        const char *msg      )
191{       // error hander must not return, so throw an exception
192        throw std::string(msg);
193}
194
195bool forward_nan(void)
196{
197
200
201        size_t n = 2, m = 1;
202        vector< AD<double> > a_x(n), a_y(m);
203        a_x[0] = 1.;
204        a_x[1] = 2.;
205        Independent(a_x);
206        a_y[0] = a_x[0] / a_x[1];
208        //
209        vector<double> x(n), y(m);
210        x[0] = 0.;
211        x[1] = 0.;
212
213        // replace the default CppAD error handler
215
216        bool ok = false;
217        try {
218                y    = f.Forward(0, x);
219        }
220        catch( std::string msg )
221        {       // check that the message contains "nan"
222                ok = msg.find("nan") != std::string::npos;
223        }
224
225        return ok;
226}
227} // END empty namespace
228
229# include <vector>
230# include <valarray>
231bool Forward(void)
232{       bool ok = true;
233        ok &= ForwardCases< CppAD::vector  <double> >();
234        ok &= ForwardCases< std::vector    <double> >();
235        ok &= ForwardCases< std::valarray  <double> >();
236        ok &= ForwardOlder();
237# ifndef NDEBUG
238        // CppAD does not check for nan when NDEBUG is defined
239        ok &= forward_nan();
240# endif
241        return ok;
242}
Note: See TracBrowser for help on using the repository browser.