source: trunk/test_more/sparse_hessian.cpp @ 1601

Last change on this file since 1601 was 1558, checked in by bradbell, 11 years ago

trunk: Extend sparse Hessian driver to optionally use set representation of sparsity.

check_doxygen.sh: add sparse_hessian.hpp to list of documented files.
*/sparse_hessian.cpp: example and test of using sets for sparse calculations.
makefile.am: add sparse_hessian to list of tests.
makefile.in: changes automatically generated by change in makefile.am.
test_more.cpp: add sparse_hessian to list of tests.
sparse_jacobina.cpp: minor edits.
whats_new_09.omh: user's view of the changes.
ad_fun.hpp: add new sparse Hessian private helper functions to class.
sparse_jacobian.hpp: minor edits.
sparse_hessian.cpp: change to allow for vectors of sets sparsity patterns.

  • Property svn:keywords set to Id
File size: 4.9 KB
Line 
1/* $Id: sparse_hessian.cpp 1558 2009-10-23 15:46:58Z bradbell $ */
2/* --------------------------------------------------------------------------
3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-09 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 SparseHessian example
15*/
16
17# include <cppad/cppad.hpp>
18namespace { // ---------------------------------------------------------
19template <class VectorBase, class VectorBool> 
20bool bool_case()
21{       bool ok = true;
22        using CppAD::AD;
23        using CppAD::NearEqual;
24        size_t i, j, k;
25
26        // domain space vector
27        size_t n = 3;
28        CPPAD_TEST_VECTOR< AD<double> >  X(n);
29        for(i = 0; i < n; i++)
30                X[i] = AD<double> (0);
31
32        // declare independent variables and starting recording
33        CppAD::Independent(X);
34
35        size_t m = 1;
36        CPPAD_TEST_VECTOR< AD<double> >  Y(m);
37        Y[0] = X[0] * X[0] + X[0] * X[1] + X[1] * X[1] + X[2] * X[2];
38
39        // create f: X -> Y and stop tape recording
40        CppAD::ADFun<double> f(X, Y);
41
42        // new value for the independent variable vector
43        VectorBase x(n);
44        for(i = 0; i < n; i++)
45                x[i] = double(i);
46
47        // second derivative of y[1]
48        VectorBase w(m);
49        w[0] = 1.;
50        VectorBase h( n * n );
51        h = f.SparseHessian(x, w);
52        /*
53            [ 2 1 0 ]
54        h = [ 1 2 0 ]
55            [ 0 0 2 ]
56        */
57        VectorBase check(n * n);
58        check[0] = 2.; check[1] = 1.; check[2] = 0.;
59        check[3] = 1.; check[4] = 2.; check[5] = 0.;
60        check[6] = 0.; check[7] = 0.; check[8] = 2.;
61        for(k = 0; k < n * n; k++)
62                ok &=  NearEqual(check[k], h[k], 1e-10, 1e-10 );
63
64        // determine the sparsity pattern p for Hessian of w^T F
65        VectorBool r(n * n);
66        for(j = 0; j < n; j++)
67        {       for(k = 0; k < n; k++)
68                        r[j * n + k] = false;
69                r[j * n + j] = true;
70        }
71        f.ForSparseJac(n, r);
72        //
73        VectorBool s(m);
74        for(i = 0; i < m; i++)
75                s[i] = w[i] != 0;
76        VectorBool p = f.RevSparseHes(n, s);
77
78        // test passing sparsity pattern
79        h = f.SparseHessian(x, w, p);
80        for(k = 0; k < n * n; k++)
81                ok &=  NearEqual(check[k], h[k], 1e-10, 1e-10 );
82
83        return ok;
84}
85template <class VectorBase, class VectorSet> 
86bool set_case()
87{       bool ok = true;
88        using CppAD::AD;
89        using CppAD::NearEqual;
90        size_t i, j, k;
91
92        // domain space vector
93        size_t n = 3;
94        CPPAD_TEST_VECTOR< AD<double> >  X(n);
95        for(i = 0; i < n; i++)
96                X[i] = AD<double> (0);
97
98        // declare independent variables and starting recording
99        CppAD::Independent(X);
100
101        size_t m = 1;
102        CPPAD_TEST_VECTOR< AD<double> >  Y(m);
103        Y[0] = X[0] * X[0] + X[0] * X[1] + X[1] * X[1] + X[2] * X[2];
104
105        // create f: X -> Y and stop tape recording
106        CppAD::ADFun<double> f(X, Y);
107
108        // new value for the independent variable vector
109        VectorBase x(n);
110        for(i = 0; i < n; i++)
111                x[i] = double(i);
112
113        // second derivative of y[1]
114        VectorBase w(m);
115        w[0] = 1.;
116        VectorBase h( n * n );
117        h = f.SparseHessian(x, w);
118        /*
119            [ 2 1 0 ]
120        h = [ 1 2 0 ]
121            [ 0 0 2 ]
122        */
123        VectorBase check(n * n);
124        check[0] = 2.; check[1] = 1.; check[2] = 0.;
125        check[3] = 1.; check[4] = 2.; check[5] = 0.;
126        check[6] = 0.; check[7] = 0.; check[8] = 2.;
127        for(k = 0; k < n * n; k++)
128                ok &=  NearEqual(check[k], h[k], 1e-10, 1e-10 );
129
130        // determine the sparsity pattern p for Hessian of w^T F
131        VectorSet r(n);
132        for(j = 0; j < n; j++)
133                r[j].insert(j);
134        f.ForSparseJac(n, r);
135        //
136        VectorSet s(1);
137        for(i = 0; i < m; i++)
138                if( w[i] != 0 )
139                        s[0].insert(i);
140        VectorSet p = f.RevSparseHes(n, s);
141
142        // test passing sparsity pattern
143        h = f.SparseHessian(x, w, p);
144        for(k = 0; k < n * n; k++)
145                ok &=  NearEqual(check[k], h[k], 1e-10, 1e-10 );
146
147        return ok;
148}
149} // End empty namespace
150# include <vector>
151# include <valarray>
152bool sparse_hessian(void)
153{       bool ok = true;
154        // vector of bool cases
155        ok &= bool_case< CppAD::vector  <double>, CppAD::vectorBool   >();
156        ok &= bool_case< std::vector    <double>, CppAD::vector<bool> >();
157        ok &= bool_case< std::valarray  <double>, std::vector<bool>   >();
158        // ---------------------------------------------------------------
159        // vector of set cases
160        typedef std::vector< std::set<size_t> >   std_vector_set;
161        typedef CppAD::vector< std::set<size_t> > cppad_vector_set;
162        //
163        ok &= set_case< CppAD::vector<double>, std_vector_set   >();
164        ok &= set_case< std::valarray<double>, std_vector_set   >();
165        ok &= set_case< std::vector<double>,   cppad_vector_set >();
166        ok &= set_case< CppAD::vector<double>, cppad_vector_set >();
167        //
168        // According to section 26.3.2.3 of the 1998 C++ standard
169        // a const valarray does not return references to its elements.
170        // so do not include it in the testing for sets.
171        // ---------------------------------------------------------------
172        return ok;
173}
Note: See TracBrowser for help on using the repository browser.