source: trunk/test_more/sparse_hessian.cpp @ 2401

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

Merge in changes from branches/sparse.

  1. Update version number to today.
  2. Add sparse return user API interfaces to SparseJacobian? and SparseHessian?.
  3. Include new interface and imporve sparse Hessian and Jacobian examples.
  4. More testing of sparse Hessian and Jacobians.
  5. Add is_element to sparse_set and sparse_pack.

.: merge information.
search.sh: search cppad/speed directory, drop makefile.in from results.
svn_merge.sh: parameters for this merge.
parallel_ad.hpp: initilize is_element statics.
sparse_evaluate.hpp: fix typo in documentation.
glossary.omh: convert to modern omhelp, improve AD Function entry.
link_sparse_hessian.cpp: fix typo in documentation.

  • Property svn:keywords set to Id
File size: 6.7 KB
Line 
1/* $Id: sparse_hessian.cpp 2401 2012-05-24 16:30:25Z 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                    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 { // ---------------------------------------------------------
19
20bool rc_tridiagonal(void)
21{       bool ok = true;
22        using CppAD::AD;
23        using CppAD::NearEqual;
24        size_t i, j, k, ell;
25        double eps = 10. * CppAD::epsilon<double>();
26
27        size_t n = 12; // must be greater than or equal 3; see n_sweep.
28        size_t m = n - 1;
29        CPPAD_TEST_VECTOR< AD<double> > a_x(n), a_y(m);
30        CPPAD_TEST_VECTOR<double> x(n), check(n * n), w(m);
31        for(j = 0; j < n; j++)
32                a_x[j] = x[j] = double(j+1);
33
34        // declare independent variables and start taping
35        CppAD::Independent(a_x);
36
37        for(ell = 0; ell < n * n; ell++)
38                check[ell] = 0.0;
39
40        for(i = 0; i < m; i++)
41        {       AD<double> diff = a_x[i+1] - a_x[i];   
42                a_y[i] = 0.5 * diff * diff / double(i+2); 
43                w[i] = double(i+1);
44                ell         = i * n + i;
45                check[ell] += w[i] / double(i+2);
46                ell         = (i+1) * n + i+1;
47                check[ell] += w[i] / double(i+2);
48                ell         = (i+1) * n + i;
49                check[ell] -= w[i] / double(i+2);
50                ell         = i * n + i+1;
51                check[ell] -= w[i] / double(i+2);
52        }
53
54        // create f: x -> y
55        CppAD::ADFun<double> f(a_x, a_y);
56                 
57        // determine the sparsity pattern p for Hessian of w^T f
58        typedef CppAD::vector< std::set<size_t> > VectorSet;
59        VectorSet p_r(n);
60        for(j = 0; j < n; j++)
61                p_r[j].insert(j);
62        f.ForSparseJac(n, p_r);
63        //
64        VectorSet p_s(1);
65        for(i = 0; i < m; i++)
66                if( w[i] != 0 )
67                        p_s[0].insert(i);
68        VectorSet p_h = f.RevSparseHes(n, p_s);
69
70        // requres the upper triangle of the Hessian
71        size_t K = 2 * n - 1;
72        CPPAD_TEST_VECTOR<size_t> r(K), c(K);
73        CPPAD_TEST_VECTOR<double> hes(K);
74        k = 0;
75        for(i = 0; i < n; i++)
76        {       r[k] = i;
77                c[k] = i;
78                k++;
79                if( i < n-1 )
80                {       r[k] = i;
81                        c[k] = i+1;
82                        k++;
83                }
84        } 
85        ok &= k == K;
86
87        // test computing sparse Hessian
88        CppAD::sparse_hessian_work work;
89        size_t n_sweep = f.SparseHessian(x, w, p_h, r, c, hes, work);
90        ok &= n_sweep == 3;
91        for(k = 0; k < K; k++)
92        {       ell = r[k] * n + c[k];
93                ok &=  NearEqual(check[ell], hes[k], eps, eps);
94                ell = c[k] * n + r[k];
95                ok &=  NearEqual(check[ell], hes[k], eps, eps);
96        }
97
98        return ok;
99}
100
101
102template <class VectorBase, class VectorBool> 
103bool bool_case()
104{       bool ok = true;
105        using CppAD::AD;
106        using CppAD::NearEqual;
107        size_t i, j, k;
108
109        // domain space vector
110        size_t n = 3;
111        CPPAD_TEST_VECTOR< AD<double> >  X(n);
112        for(i = 0; i < n; i++)
113                X[i] = AD<double> (0);
114
115        // declare independent variables and starting recording
116        CppAD::Independent(X);
117
118        size_t m = 1;
119        CPPAD_TEST_VECTOR< AD<double> >  Y(m);
120        Y[0] = X[0] * X[0] + X[0] * X[1] + X[1] * X[1] + X[2] * X[2];
121
122        // create f: X -> Y and stop tape recording
123        CppAD::ADFun<double> f(X, Y);
124
125        // new value for the independent variable vector
126        VectorBase x(n);
127        for(i = 0; i < n; i++)
128                x[i] = double(i);
129
130        // second derivative of y[1]
131        VectorBase w(m);
132        w[0] = 1.;
133        VectorBase h( n * n );
134        h = f.SparseHessian(x, w);
135        /*
136            [ 2 1 0 ]
137        h = [ 1 2 0 ]
138            [ 0 0 2 ]
139        */
140        VectorBase check(n * n);
141        check[0] = 2.; check[1] = 1.; check[2] = 0.;
142        check[3] = 1.; check[4] = 2.; check[5] = 0.;
143        check[6] = 0.; check[7] = 0.; check[8] = 2.;
144        for(k = 0; k < n * n; k++)
145                ok &=  NearEqual(check[k], h[k], 1e-10, 1e-10 );
146
147        // determine the sparsity pattern p for Hessian of w^T F
148        VectorBool r(n * n);
149        for(j = 0; j < n; j++)
150        {       for(k = 0; k < n; k++)
151                r[j * n + k] = false;
152                r[j * n + j] = true;
153        }
154        f.ForSparseJac(n, r);
155        //
156        VectorBool s(m);
157        for(i = 0; i < m; i++)
158                s[i] = w[i] != 0;
159        VectorBool p = f.RevSparseHes(n, s);
160
161        // test passing sparsity pattern
162        h = f.SparseHessian(x, w, p);
163        for(k = 0; k < n * n; k++)
164                ok &=  NearEqual(check[k], h[k], 1e-10, 1e-10 );
165
166        return ok;
167}
168template <class VectorBase, class VectorSet> 
169bool set_case()
170{       bool ok = true;
171        using CppAD::AD;
172        using CppAD::NearEqual;
173        size_t i, j, k;
174
175        // domain space vector
176        size_t n = 3;
177        CPPAD_TEST_VECTOR< AD<double> >  X(n);
178        for(i = 0; i < n; i++)
179                X[i] = AD<double> (0);
180
181        // declare independent variables and starting recording
182        CppAD::Independent(X);
183
184        size_t m = 1;
185        CPPAD_TEST_VECTOR< AD<double> >  Y(m);
186        Y[0] = X[0] * X[0] + X[0] * X[1] + X[1] * X[1] + X[2] * X[2];
187
188        // create f: X -> Y and stop tape recording
189        CppAD::ADFun<double> f(X, Y);
190
191        // new value for the independent variable vector
192        VectorBase x(n);
193        for(i = 0; i < n; i++)
194                x[i] = double(i);
195
196        // second derivative of y[1]
197        VectorBase w(m);
198        w[0] = 1.;
199        VectorBase h( n * n );
200        h = f.SparseHessian(x, w);
201        /*
202            [ 2 1 0 ]
203        h = [ 1 2 0 ]
204            [ 0 0 2 ]
205        */
206        VectorBase check(n * n);
207        check[0] = 2.; check[1] = 1.; check[2] = 0.;
208        check[3] = 1.; check[4] = 2.; check[5] = 0.;
209        check[6] = 0.; check[7] = 0.; check[8] = 2.;
210        for(k = 0; k < n * n; k++)
211                ok &=  NearEqual(check[k], h[k], 1e-10, 1e-10 );
212
213        // determine the sparsity pattern p for Hessian of w^T F
214        VectorSet r(n);
215        for(j = 0; j < n; j++)
216                r[j].insert(j);
217        f.ForSparseJac(n, r);
218        //
219        VectorSet s(1);
220        for(i = 0; i < m; i++)
221                if( w[i] != 0 )
222                        s[0].insert(i);
223        VectorSet p = f.RevSparseHes(n, s);
224
225        // test passing sparsity pattern
226        h = f.SparseHessian(x, w, p);
227        for(k = 0; k < n * n; k++)
228                ok &=  NearEqual(check[k], h[k], 1e-10, 1e-10 );
229
230        return ok;
231}
232} // End empty namespace
233# include <vector>
234# include <valarray>
235bool sparse_hessian(void)
236{       bool ok = true;
237
238        ok &= rc_tridiagonal();
239        // ---------------------------------------------------------------
240        // vector of bool cases
241        ok &= bool_case< CppAD::vector  <double>, CppAD::vectorBool   >();
242        ok &= bool_case< std::vector    <double>, CppAD::vector<bool> >();
243        ok &= bool_case< std::valarray  <double>, std::vector<bool>   >();
244        // ---------------------------------------------------------------
245        // vector of set cases
246        typedef std::vector< std::set<size_t> >   std_vector_set;
247        typedef CppAD::vector< std::set<size_t> > cppad_vector_set;
248        //
249        ok &= set_case< CppAD::vector<double>, std_vector_set   >();
250        ok &= set_case< std::valarray<double>, std_vector_set   >();
251        ok &= set_case< std::vector<double>,   cppad_vector_set >();
252        ok &= set_case< CppAD::vector<double>, cppad_vector_set >();
253        //
254        // According to section 26.3.2.3 of the 1998 C++ standard
255        // a const valarray does not return references to its elements.
256        // so do not include it in the testing for sets.
257        // ---------------------------------------------------------------
258        return ok;
259}
Note: See TracBrowser for help on using the repository browser.