source: trunk/test_more/general/forward.cpp @ 3932

Last change on this file since 3932 was 3932, checked in by bradbell, 2 years ago

merge to branch: trunk
from repository: https://github.com/coin-or/CppAD
start hash code: f4ce6b2601ca057b41100ab6787a8d9cb178e945
end hash code: 2c4a32a99cc6fd35150317d891fc837423e7d37f

commit 2c4a32a99cc6fd35150317d891fc837423e7d37f
Author: Brad Bell <bradbell@…>
Date: Fri May 19 07:36:14 2017 -0700

Move test_more/*.cpp -> test_more/general/*.cpp
(in preparation for other sub-directories of test-more).

commit 4b0be083349be333543036d8e77623285f92c6ec
Author: Brad Bell <bradbell@…>
Date: Fri May 19 05:48:50 2017 -0700

push_git2svn.py: improve detection of moved files.

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