source: trunk/test_more/for_hess.cpp @ 2506

Last change on this file since 2506 was 2506, checked in by bradbell, 8 years ago

Change Licenses: CPL-1.0 -> EPL-1.0, GPL-2.0->GPL-3.0

  • Property svn:keywords set to Id
File size: 3.1 KB
Line 
1/* $Id: for_hess.cpp 2506 2012-10-24 19:36:49Z bradbell $ */
2/* --------------------------------------------------------------------------
3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 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
14/*
15Old ForHess example now used just for valiadation testing
16*/
17
18# include <cppad/cppad.hpp>
19
20
21bool ForHess(void)
22{       bool ok = true;
23
24        using namespace CppAD;
25        using CppAD::exp;
26        using CppAD::sin;
27        using CppAD::cos;
28
29        size_t i;
30
31        // create independent variable vector with assigned values
32        CPPAD_TESTVECTOR(double)      u0(3);
33        CPPAD_TESTVECTOR(AD<double>) U(3);
34        for(i = 0; i < 3; i++)
35                U[i] = u0[i] = double(i+1);
36        Independent( U );
37
38        // define the function
39        CPPAD_TESTVECTOR(AD<double>) Y(2);
40        Y[0] = U[0] * exp( U[1] );
41        Y[1] = U[1] * sin( U[2] ); 
42
43        // create the function y = F(u)
44        ADFun<double> F(U, Y);
45
46        // formulas for the upper triangle of Hessian of F_0
47        CPPAD_TESTVECTOR(double) H0(9);
48        H0[0] = 0.;                    // d^2 y[0] / d_u[0] d_u[0]
49        H0[1] = exp( u0[1] );          // d^2 y[0] / d_u[0] d_u[1]
50        H0[2] = 0.;                    // d^2 y[0] / d_u[0] d_u[2]
51
52        H0[4] = u0[0] * exp( u0[1] );  // d^2 y[0] / d_u[1] d_u[1]
53        H0[5] = 0.;                    // d^2 y[0] / d_u[1] d_u[2]
54
55        H0[8] = 0.;                    // d^2 y[0] / d_u[2] d_u[2]
56
57        // formulas for the upper triangle of Hessian of F_1
58        CPPAD_TESTVECTOR(double) H1(9);
59        H1[0] = 0.;                    // d^2 Y[1] / d_U[0] d_U[0]
60        H1[1] = 0.;                    // d^2 Y[1] / d_U[0] d_U[1]
61        H1[2] = 0.;                    // d^2 Y[1] / d_U[0] d_U[2]
62
63        H1[4] = 0.;                    // d^2 Y[1] / d_U[1] d_U[1]
64        H1[5] = cos( u0[2] );          // d^2 Y[1] / d_U[1] d_U[2]
65
66        H1[8] = - u0[1] * sin( u0[2] );// d^2 Y[1] / d_U[2] d_U[2]
67
68
69        // Define U(t) = u0 + u1 t + u2 t^2 / 2
70        CPPAD_TESTVECTOR(double) u1(3);
71        CPPAD_TESTVECTOR(double) u2(3);
72        for(i = 0; i < 3; i++)
73                u1[i] = u2[i] = 0.;
74
75        size_t j;
76        for(i = 0; i < 3; i++)
77        {       // diagonal of Hessians in i-th coordiante direction
78                u1[i] = 1.;
79                F.Forward(1, u1);
80                CPPAD_TESTVECTOR(double) Di = F.Forward(2, u2);
81                ok &= NearEqual( 2. * Di[0] , H0[ i + 3 * i], 1e-10, 1e-10);
82                ok &= NearEqual( 2. * Di[1] , H1[ i + 3 * i], 1e-10, 1e-10);
83                //
84                for(j = i+1; j < 3; j++)
85                {       // cross term in i and j direction
86                        u1[j] = 1.;
87                        F.Forward(1, u1);
88                        CPPAD_TESTVECTOR(double) Cij = F.Forward(2, u2);
89
90                        // diagonal of Hessian in j-th coordinate direction
91                        u1[i] = 0.;
92                        F.Forward(1, u1);
93                        CPPAD_TESTVECTOR(double) Dj = F.Forward(2, u2);
94
95                        // (i, j) elements of the Hessians
96                        double H0ij = Cij[0] - Di[0] - Dj[0];
97                        ok &= NearEqual( H0ij, H0[j + 3 * i], 1e-10, 1e-10);
98                        double H1ij = Cij[1] - Di[1] - Dj[1];
99                        ok &= NearEqual( H1ij, H1[j + 3 * i], 1e-10, 1e-10);
100
101                        // reset all components of u1 to zero
102                        u1[j] = 0.;
103                }
104        }
105
106        return ok;
107}
Note: See TracBrowser for help on using the repository browser.