1 /* --------------------------------------------------------------------------
2 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
3
4 CppAD is distributed under the terms of the
5 Eclipse Public License Version 2.0.
6
7 This Source Code may also be made available under the following
8 Secondary License when the conditions for such availability set forth
9 in the Eclipse Public License, Version 2.0 are satisfied:
10 GNU General Public License, Version 2.0 or later.
11 ---------------------------------------------------------------------------- */
12
13 /*
14 Old SparseHessian example
15 */
16
17 # include <cppad/cppad.hpp>
18 namespace { // ---------------------------------------------------------
19
rc_tridiagonal(void)20 bool rc_tridiagonal(void)
21 { bool ok = true;
22 using CppAD::AD;
23 using CppAD::NearEqual;
24 size_t i, j, k, ell;
25 double eps10 = 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> > SetVector;
59 SetVector p_r(n);
60 for(j = 0; j < n; j++)
61 p_r[j].insert(j);
62 f.ForSparseJac(n, p_r);
63 //
64 SetVector p_s(1);
65 for(i = 0; i < m; i++)
66 if( w[i] != 0 )
67 p_s[0].insert(i);
68 SetVector 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], eps10, eps10);
94 ell = c[k] * n + r[k];
95 ok &= NearEqual(check[ell], hes[k], eps10, eps10);
96 }
97
98 return ok;
99 }
100
101
102 template <class BaseVector, class BoolVector>
bool_case()103 bool bool_case()
104 { bool ok = true;
105 using CppAD::AD;
106 using CppAD::NearEqual;
107 size_t i, j, k;
108 double eps99 = 99.0 * std::numeric_limits<double>::epsilon();
109
110 // domain space vector
111 size_t n = 3;
112 CPPAD_TESTVECTOR(AD<double>) X(n);
113 for(i = 0; i < n; i++)
114 X[i] = AD<double> (0);
115
116 // declare independent variables and starting recording
117 CppAD::Independent(X);
118
119 size_t m = 1;
120 CPPAD_TESTVECTOR(AD<double>) Y(m);
121 Y[0] = X[0] * X[0] + X[0] * X[1] + X[1] * X[1] + X[2] * X[2];
122
123 // create f: X -> Y and stop tape recording
124 CppAD::ADFun<double> f(X, Y);
125
126 // new value for the independent variable vector
127 BaseVector x(n);
128 for(i = 0; i < n; i++)
129 x[i] = double(i);
130
131 // second derivative of y[1]
132 BaseVector w(m);
133 w[0] = 1.;
134 BaseVector h( n * n );
135 h = f.SparseHessian(x, w);
136 /*
137 [ 2 1 0 ]
138 h = [ 1 2 0 ]
139 [ 0 0 2 ]
140 */
141 BaseVector check(n * n);
142 check[0] = 2.; check[1] = 1.; check[2] = 0.;
143 check[3] = 1.; check[4] = 2.; check[5] = 0.;
144 check[6] = 0.; check[7] = 0.; check[8] = 2.;
145 for(k = 0; k < n * n; k++)
146 ok &= NearEqual(check[k], h[k], eps99, eps99 );
147
148 // determine the sparsity pattern p for Hessian of w^T F
149 BoolVector r(n * n);
150 for(j = 0; j < n; j++)
151 { for(k = 0; k < n; k++)
152 r[j * n + k] = false;
153 r[j * n + j] = true;
154 }
155 f.ForSparseJac(n, r);
156 //
157 BoolVector s(m);
158 for(i = 0; i < m; i++)
159 s[i] = w[i] != 0;
160 BoolVector p = f.RevSparseHes(n, s);
161
162 // test passing sparsity pattern
163 h = f.SparseHessian(x, w, p);
164 for(k = 0; k < n * n; k++)
165 ok &= NearEqual(check[k], h[k], eps99, eps99 );
166
167 return ok;
168 }
169 template <class BaseVector, class SetVector>
set_case()170 bool set_case()
171 { bool ok = true;
172 using CppAD::AD;
173 using CppAD::NearEqual;
174 size_t i, j, k;
175 double eps99 = 99.0 * std::numeric_limits<double>::epsilon();
176
177 // domain space vector
178 size_t n = 3;
179 CPPAD_TESTVECTOR(AD<double>) X(n);
180 for(i = 0; i < n; i++)
181 X[i] = AD<double> (0);
182
183 // declare independent variables and starting recording
184 CppAD::Independent(X);
185
186 size_t m = 1;
187 CPPAD_TESTVECTOR(AD<double>) Y(m);
188 Y[0] = X[0] * X[0] + X[0] * X[1] + X[1] * X[1] + X[2] * X[2];
189
190 // create f: X -> Y and stop tape recording
191 CppAD::ADFun<double> f(X, Y);
192
193 // new value for the independent variable vector
194 BaseVector x(n);
195 for(i = 0; i < n; i++)
196 x[i] = double(i);
197
198 // second derivative of y[1]
199 BaseVector w(m);
200 w[0] = 1.;
201 BaseVector h( n * n );
202 h = f.SparseHessian(x, w);
203 /*
204 [ 2 1 0 ]
205 h = [ 1 2 0 ]
206 [ 0 0 2 ]
207 */
208 BaseVector check(n * n);
209 check[0] = 2.; check[1] = 1.; check[2] = 0.;
210 check[3] = 1.; check[4] = 2.; check[5] = 0.;
211 check[6] = 0.; check[7] = 0.; check[8] = 2.;
212 for(k = 0; k < n * n; k++)
213 ok &= NearEqual(check[k], h[k], eps99, eps99 );
214
215 // determine the sparsity pattern p for Hessian of w^T F
216 SetVector r(n);
217 for(j = 0; j < n; j++)
218 r[j].insert(j);
219 f.ForSparseJac(n, r);
220 //
221 SetVector s(1);
222 for(i = 0; i < m; i++)
223 if( w[i] != 0 )
224 s[0].insert(i);
225 SetVector p = f.RevSparseHes(n, s);
226
227 // test passing sparsity pattern
228 h = f.SparseHessian(x, w, p);
229 for(k = 0; k < n * n; k++)
230 ok &= NearEqual(check[k], h[k], eps99, eps99 );
231
232 return ok;
233 }
234 } // End empty namespace
235 # include <vector>
236 # include <valarray>
sparse_hessian(void)237 bool sparse_hessian(void)
238 { bool ok = true;
239
240 ok &= rc_tridiagonal();
241 // ---------------------------------------------------------------
242 // vector of bool cases
243 ok &= bool_case< CppAD::vector <double>, CppAD::vectorBool >();
244 ok &= bool_case< std::vector <double>, CppAD::vector<bool> >();
245 ok &= bool_case< std::valarray <double>, std::vector<bool> >();
246 // ---------------------------------------------------------------
247 // vector of set cases
248 typedef std::vector< std::set<size_t> > std_vector_set;
249 typedef CppAD::vector< std::set<size_t> > cppad_vector_set;
250 //
251 ok &= set_case< CppAD::vector<double>, std_vector_set >();
252 ok &= set_case< std::valarray<double>, std_vector_set >();
253 ok &= set_case< std::vector<double>, cppad_vector_set >();
254 ok &= set_case< CppAD::vector<double>, cppad_vector_set >();
255 //
256 // According to section 26.3.2.3 of the 1998 C++ standard
257 // a const valarray does not return references to its elements.
258 // so do not include it in the testing for sets.
259 // ---------------------------------------------------------------
260 return ok;
261 }
262