1 /* --------------------------------------------------------------------------
2 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-20 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 $begin colpack_hessian.cpp$$
15 $spell
16 colpack_hessian
17 jacobian
18 $$
19
20 $section ColPack: Sparse Hessian Example and Test$$
21
22
23 $srcthisfile%0%// BEGIN C++%// END C++%1%$$
24
25 $end
26 */
27 // BEGIN C++
28
29 # include <cppad/cppad.hpp>
colpack_hessian(void)30 bool colpack_hessian(void)
31 { bool ok = true;
32 using CppAD::AD;
33 using CppAD::NearEqual;
34 typedef CPPAD_TESTVECTOR(AD<double>) a_vector;
35 typedef CPPAD_TESTVECTOR(double) d_vector;
36 typedef CppAD::vector<size_t> i_vector;
37 size_t i, j, k, ell;
38 double eps = 10. * CppAD::numeric_limits<double>::epsilon();
39
40 // domain space vector
41 size_t n = 5;
42 a_vector a_x(n);
43 for(j = 0; j < n; j++)
44 a_x[j] = AD<double> (0);
45
46 // declare independent variables and starting recording
47 CppAD::Independent(a_x);
48
49 // colpack example case where hessian is a spear head
50 // i.e, H(i, j) non zero implies i = 0, j = 0, or i = j
51 AD<double> sum = 0.0;
52 // partial_0 partial_j = x[j]
53 // partial_j partial_j = x[0]
54 for(j = 1; j < n; j++)
55 sum += a_x[0] * a_x[j] * a_x[j] / 2.0;
56 //
57 // partial_i partial_i = 2 * x[i]
58 for(i = 0; i < n; i++)
59 sum += a_x[i] * a_x[i] * a_x[i] / 3.0;
60
61 // declare dependent variables
62 size_t m = 1;
63 a_vector a_y(m);
64 a_y[0] = sum;
65
66 // create f: x -> y and stop tape recording
67 CppAD::ADFun<double> f(a_x, a_y);
68
69 // new value for the independent variable vector
70 d_vector x(n);
71 for(j = 0; j < n; j++)
72 x[j] = double(j + 1);
73
74 /*
75 [ 2 2 3 4 5 ]
76 hes = [ 2 5 0 0 0 ]
77 [ 3 0 7 0 0 ]
78 [ 4 0 0 9 0 ]
79 [ 5 0 0 0 11 ]
80 */
81 d_vector check(n * n);
82 for(i = 0; i < n; i++)
83 { for(j = 0; j < n; j++)
84 { size_t index = i * n + j;
85 check[index] = 0.0;
86 if( i == 0 && 1 <= j )
87 check[index] += x[j];
88 if( 1 <= i && j == 0 )
89 check[index] += x[i];
90 if( i == j )
91 { check[index] += 2.0 * x[i];
92 if( i != 0 )
93 check[index] += x[0];
94 }
95 }
96 }
97 // Normally one would use f.RevSparseHes to compute
98 // sparsity pattern, but for this example we extract it from check.
99 std::vector< std::set<size_t> > p(n);
100 i_vector row, col;
101 for(i = 0; i < n; i++)
102 { for(j = 0; j < n; j++)
103 { ell = i * n + j;
104 if( check[ell] != 0. )
105 { // insert this non-zero entry in sparsity pattern
106 p[i].insert(j);
107
108 // the Hessian is symmetric, so only lower triangle
109 if( j <= i )
110 { row.push_back(i);
111 col.push_back(j);
112 }
113 }
114 }
115 }
116 size_t K = row.size();
117 d_vector hes(K);
118
119 // default coloring method is cppad.symmetric
120 CppAD::sparse_hessian_work work;
121 ok &= work.color_method == "cppad.symmetric";
122
123 // contrast and check results for both CppAD and Colpack
124 for(size_t i_method = 0; i_method < 4; i_method++)
125 { // empty work structure
126 switch(i_method)
127 { case 0:
128 work.color_method = "cppad.symmetric";
129 break;
130
131 case 1:
132 work.color_method = "cppad.general";
133 break;
134
135 case 2:
136 work.color_method = "colpack.symmetric";
137 break;
138
139 case 3:
140 work.color_method = "colpack.general";
141 break;
142 }
143
144 // compute Hessian
145 d_vector w(m);
146 w[0] = 1.0;
147 size_t n_sweep = f.SparseHessian(x, w, p, row, col, hes, work);
148 //
149 // check result
150 for(k = 0; k < K; k++)
151 { ell = row[k] * n + col[k];
152 ok &= NearEqual(check[ell], hes[k], eps, eps);
153 }
154 if(
155 work.color_method == "cppad.symmetric"
156 || work.color_method == "colpack.symmetric"
157 )
158 ok &= n_sweep == 2;
159 else
160 ok &= n_sweep == 5;
161 //
162 // check that clear resets color_method to cppad.symmetric
163 work.clear();
164 ok &= work.color_method == "cppad.symmetric";
165 }
166
167 return ok;
168 }
169 // END C++
170