source: trunk/test_more/sparse_hessian.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: 6.7 KB
Line 
1/* $Id: sparse_hessian.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/*
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_TESTVECTOR(AD<double>) a_x(n), a_y(m);
30        CPPAD_TESTVECTOR(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_TESTVECTOR(size_t) r(K), c(K);
73        CPPAD_TESTVECTOR(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_TESTVECTOR(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_TESTVECTOR(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_TESTVECTOR(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_TESTVECTOR(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.