source: trunk/test_more/general/sparse_sub_hes.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.

File size: 5.7 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/*
13$begin sparse_sub_hes.cpp$$
14$spell
15$$
16
17$section Sparse Hessian on Subset of Variables: Example and Test$$
18
19$head Purpose$$
20This example uses a
21$cref/column subset/sparse_hessian/p/Column Subset/$$ of the sparsity pattern
22to compute the Hessian for a subset of the variables.
23The values in the rest of the sparsity pattern do not matter.
24
25$head See Also$$
26$cref sub_sparse_hes.cpp$$
27
28$code
29$comment%example/sparse/sparse_sub_hes.cpp%0%// BEGIN C++%// END C++%1%$$
30$$
31
32$end
33*/
34// BEGIN C++
35# include <cppad/cppad.hpp>
36namespace { // BEGIN_EMPTY_NAMESPACE
37
38// --------------------------------------------------------------------------
39CppAD::ADFun<double> record_function(size_t n)
40{       // must be greater than or equal 3; see n_sweep below
41        assert( n >= 3 );
42        //
43        using CppAD::AD;
44        typedef CppAD::vector< AD<double> >     a_vector;
45        //
46        // domain space vector
47        a_vector a_x(n);
48        for(size_t j = 0; j < n; j++)
49                a_x[j] = AD<double> (0);
50
51        // declare independent variables and starting recording
52        CppAD::Independent(a_x);
53
54        // range space vector
55        size_t m = 1;
56        a_vector a_y(m);
57        a_y[0] = 0.0;
58        for(size_t j = 1; j < n; j++)
59                a_y[0] += a_x[j-1] * a_x[j] * a_x[j];
60
61        // create f: x -> y and stop tape recording
62        // (without executing zero order forward calculation)
63        CppAD::ADFun<double> f;
64        f.Dependent(a_x, a_y);
65        //
66        return f;
67}
68// --------------------------------------------------------------------------
69bool test_set(const char* color_method)
70{       bool ok = true;
71        //
72        typedef CppAD::vector< double >                   d_vector;
73        typedef CppAD::vector<size_t>                     i_vector;
74        typedef CppAD::vector< std::set<size_t> >         s_vector;
75        //
76        size_t n = 12;
77        CppAD::ADFun<double> f = record_function(n);
78        //
79        // sparsity patteren for the sub-set of variables we are computing
80        // the hessian w.r.t.
81        size_t n_sub = 4;
82        s_vector r(n);
83        for(size_t j = 0; j < n_sub; j++)
84        {       assert(  r[j].empty() );
85                r[j].insert(j);
86        }
87
88        // store forward sparsity for J(x) = F^{(1)} (x) * R
89        f.ForSparseJac(n_sub, r);
90
91        // compute sparsity pattern for H(x) = (S * F)^{(2)} ( x ) * R
92        s_vector s(1);
93        assert(  s[0].empty() );
94        s[0].insert(0);
95        bool transpose = true;
96        s_vector h = f.RevSparseHes(n_sub, s, transpose);
97
98        // set the row and column indices that correspond to lower triangle
99        i_vector row, col;
100        for(size_t i = 0; i < n_sub; i++)
101        {       if( i > 0 )
102                {       // diagonal element
103                        row.push_back(i);
104                        col.push_back(i);
105                        // lower diagonal element
106                        row.push_back(i);
107                        col.push_back(i-1);
108                }
109        }
110
111        // weighting for the Hessian
112        d_vector w(1);
113        w[0] = 1.0;
114
115        // compute Hessian
116        CppAD::sparse_hessian_work work;
117        work.color_method = color_method;
118        d_vector x(n), hes( row.size() );
119        for(size_t j = 0; j < n; j++)
120                x[j] = double(j+1);
121        f.SparseHessian(x, w, h, row, col, hes, work);
122
123        // check the values in the sparse hessian
124        for(size_t ell = 0; ell < row.size(); ell++)
125        {       size_t i = row[ell];
126                size_t j = col[ell];
127                if( i == j )
128                        ok &= hes[ell] == 2.0 * x[i-1];
129                else
130                {       ok &= j+1 == i;
131                        ok &= hes[ell] == 2.0 * x[i];
132                }
133                ell++;
134        }
135        return ok;
136}
137// --------------------------------------------------------------------------
138bool test_bool(const char* color_method)
139{       bool ok = true;
140        //
141        typedef CppAD::vector< double >    d_vector;
142        typedef CppAD::vector<size_t>      i_vector;
143        typedef CppAD::vector<bool>        s_vector;
144        //
145        size_t n = 12;
146        CppAD::ADFun<double> f = record_function(n);
147        //
148        // sparsity patteren for the sub-set of variables we are computing
149        // the hessian w.r.t.
150        size_t n_sub = 4;
151        s_vector r(n * n_sub);
152        for(size_t i = 0; i < n; i++)
153        {       for(size_t j = 0; j < n_sub; j++)
154                        r[ i * n_sub + j ] = (i == j);
155        }
156
157        // store forward sparsity for J(x) = F^{(1)} (x) * R
158        f.ForSparseJac(n_sub, r);
159
160        // compute sparsity pattern for H(x) = (S * F)^{(2)} ( x ) * R
161        s_vector s(1);
162        s[0] = true;
163        bool transpose = true;
164        s_vector h = f.RevSparseHes(n_sub, s, transpose);
165
166        // set the row and column indices that correspond to lower triangle
167        i_vector row, col;
168        for(size_t i = 0; i < n_sub; i++)
169        {       if( i > 0 )
170                {       // diagonal element
171                        row.push_back(i);
172                        col.push_back(i);
173                        // lower diagonal element
174                        row.push_back(i);
175                        col.push_back(i-1);
176                }
177        }
178
179        // weighting for the Hessian
180        d_vector w(1);
181        w[0] = 1.0;
182
183        // extend sparsity pattern (values in extended columns do not matter)
184        s_vector h_extended(n * n);
185        for(size_t i = 0; i < n; i++)
186        {       for(size_t j = 0; j < n_sub; j++)
187                        h_extended[ i * n + j ] = h[ i * n_sub + j ];
188                for(size_t j = n_sub; j < n; j++)
189                        h_extended[ i * n + j ] = false;
190        }
191        // compute Hessian
192        CppAD::sparse_hessian_work work;
193        work.color_method = color_method;
194        d_vector x(n), hes( row.size() );
195        for(size_t j = 0; j < n; j++)
196                x[j] = double(j+1);
197        f.SparseHessian(x, w, h_extended, row, col, hes, work);
198
199        // check the values in the sparse hessian
200        for(size_t ell = 0; ell < row.size(); ell++)
201        {       size_t i = row[ell];
202                size_t j = col[ell];
203                if( i == j )
204                        ok &= hes[ell] == 2.0 * x[i-1];
205                else
206                {       ok &= j+1 == i;
207                        ok &= hes[ell] == 2.0 * x[i];
208                }
209                ell++;
210        }
211        return ok;
212}
213} // END_EMPTY_NAMESPACE
214
215bool sparse_sub_hes(void)
216{       bool ok = true;
217        ok &= test_set("cppad.symmetric");
218        ok &= test_set("cppad.general");
219        //
220        ok &= test_bool("cppad.symmetric");
221        ok &= test_bool("cppad.general");
222        return ok;
223}
224// END C++
Note: See TracBrowser for help on using the repository browser.