source: trunk/test_more/forward.cpp @ 3793

Last change on this file since 3793 was 3744, checked in by bradbell, 4 years ago

merge to branch: trunk
from repository: https://github.com/coin-or/CppAD
start hash code: 5c92f4e222d7edf44cc2633c5671418dd05473d0
end hash code: 079578bc3de0b44d7c647216b46090f45bbeca9d

commit 079578bc3de0b44d7c647216b46090f45bbeca9d
Author: Brad Bell <bradbell@…>
Date: Fri Nov 6 07:02:07 2015 -0800

Store and recover independent variables when nan occurs during Forward(0, x).


forward.hpp: remove invisible white space.
forward.cpp: remove invisible white space.

  • Property svn:keywords set to Id
File size: 6.3 KB
Line 
1/* $Id: forward.cpp 3744 2015-11-06 15:33:15Z bradbell $ */
2/* --------------------------------------------------------------------------
3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-15 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/*
14Two old Forward example now used just for valiadation testing
15*/
16
17# include <cppad/cppad.hpp>
18
19namespace { // Begin empty namespace
20
21template <typename VectorDouble> // vector class, elements of type double
22bool ForwardCases(void)
23{       bool ok = true;
24
25        using namespace CppAD;
26
27        // independent variable vector
28        CPPAD_TESTVECTOR(AD<double>) X(2);
29        X[0] = 0.;
30        X[1] = 1.;
31        Independent(X);
32
33        // compute product of elements in X
34        CPPAD_TESTVECTOR(AD<double>) Y(1);
35        Y[0] = X[0] * X[0] * X[1];
36
37        // create function object F : X -> Y
38        ADFun<double> 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
94        using namespace CppAD;
95
96        // independent variable vector
97        CPPAD_TESTVECTOR(AD<double>) U(3);
98        U[0] = 0.; U[1] = 1.; U[2] = 2.;
99        Independent(U);
100
101        // compute sum and product of elements in U
102        AD<double> sum  = 0.;
103        AD<double> prod = 1.;
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
111        CPPAD_TESTVECTOR(AD<double>) V(2);
112        V[0] = sum;
113        V[1] = prod;
114
115        // V = f(U)
116        ADFun<double> f(U, V);
117
118        // use ADFun object to evaluate f[ (1, 2, 3)^T ] -----------------
119        CPPAD_TESTVECTOR(double) u0( f.Domain() );
120        CPPAD_TESTVECTOR(double) v0( f.Range() );
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 ]
127        CPPAD_TESTVECTOR(double) f0(2);
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 -----------------
136        CPPAD_TESTVECTOR(double) u1( f.Domain() );
137        CPPAD_TESTVECTOR(double) v1( f.Range() );
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
143        CPPAD_TESTVECTOR(double) g0(3), g1(3);
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 }
155        CPPAD_TESTVECTOR(double) u2( f.Domain() );
156        CPPAD_TESTVECTOR(double) v2( f.Range() );
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)
163        CPPAD_TESTVECTOR(double) H1(9);
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        std::string message = msg;
193        throw message;
194}
195
196bool forward_nan(void)
197{
198
199        using CppAD::vector;
200        using CppAD::AD;
201
202        size_t n = 2, m = 1;
203        vector< AD<double> > a_x(n), a_y(m);
204        a_x[0] = 1.;
205        a_x[1] = 2.;
206        Independent(a_x);
207        a_y[0] = a_x[0] / a_x[1];
208        CppAD::ADFun<double> f(a_x, a_y);
209        //
210        vector<double> x(n), y(m);
211        x[0] = 0.;
212        x[1] = 0.;
213
214        // replace the default CppAD error handler
215        CppAD::ErrorHandler info(my_error_handler);
216
217        bool ok = false;
218        try {
219                y    = f.Forward(0, x);
220        }
221        catch( std::string msg )
222        {       // check that the message contains
223                // "vector_size = " and "file_name = "
224                ok = msg.find("vector_size = ") != std::string::npos;
225                ok = msg.find("file_name = ") != std::string::npos;
226        }
227
228        return ok;
229}
230} // END empty namespace
231
232# include <vector>
233# include <valarray>
234bool Forward(void)
235{       bool ok = true;
236        ok &= ForwardCases< CppAD::vector  <double> >();
237        ok &= ForwardCases< std::vector    <double> >();
238        ok &= ForwardCases< std::valarray  <double> >();
239        ok &= ForwardOlder();
240# ifndef NDEBUG
241        // CppAD does not check for nan when NDEBUG is defined
242        ok &= forward_nan();
243# endif
244        return ok;
245}
Note: See TracBrowser for help on using the repository browser.