1 | # ifndef CPPAD_CORE_SPARSE_HESSIAN_HPP |
---|
2 | # define CPPAD_CORE_SPARSE_HESSIAN_HPP |
---|
3 | |
---|
4 | /* -------------------------------------------------------------------------- |
---|
5 | CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell |
---|
6 | |
---|
7 | CppAD is distributed under multiple licenses. This distribution is under |
---|
8 | the terms of the |
---|
9 | Eclipse Public License Version 1.0. |
---|
10 | |
---|
11 | A copy of this license is included in the COPYING file of this distribution. |
---|
12 | Please visit http://www.coin-or.org/CppAD/ for information on other licenses. |
---|
13 | -------------------------------------------------------------------------- */ |
---|
14 | |
---|
15 | /* |
---|
16 | $begin sparse_hessian$$ |
---|
17 | $spell |
---|
18 | jacobian |
---|
19 | recomputed |
---|
20 | CppAD |
---|
21 | valarray |
---|
22 | std |
---|
23 | Bool |
---|
24 | hes |
---|
25 | const |
---|
26 | Taylor |
---|
27 | cppad |
---|
28 | cmake |
---|
29 | colpack |
---|
30 | $$ |
---|
31 | |
---|
32 | $section Sparse Hessian$$ |
---|
33 | |
---|
34 | $head Syntax$$ |
---|
35 | $icode%hes% = %f%.SparseHessian(%x%, %w%) |
---|
36 | %hes% = %f%.SparseHessian(%x%, %w%, %p%) |
---|
37 | %n_sweep% = %f%.SparseHessian(%x%, %w%, %p%, %row%, %col%, %hes%, %work%) |
---|
38 | %$$ |
---|
39 | |
---|
40 | $head Purpose$$ |
---|
41 | We use $latex n$$ for the $cref/domain/seq_property/Domain/$$ size, |
---|
42 | and $latex m$$ for the $cref/range/seq_property/Range/$$ size of $icode f$$. |
---|
43 | We use $latex F : \B{R}^n \rightarrow \B{R}^m$$ do denote the |
---|
44 | $cref/AD function/glossary/AD Function/$$ |
---|
45 | corresponding to $icode f$$. |
---|
46 | The syntax above sets $icode hes$$ to the Hessian |
---|
47 | $latex \[ |
---|
48 | H(x) = \dpow{2}{x} \sum_{i=1}^m w_i F_i (x) |
---|
49 | \] $$ |
---|
50 | This routine takes advantage of the sparsity of the Hessian |
---|
51 | in order to reduce the amount of computation necessary. |
---|
52 | If $icode row$$ and $icode col$$ are present, it also takes |
---|
53 | advantage of the reduced set of elements of the Hessian that |
---|
54 | need to be computed. |
---|
55 | One can use speed tests (e.g. $cref speed_test$$) |
---|
56 | to verify that results are computed faster |
---|
57 | than when using the routine $cref Hessian$$. |
---|
58 | |
---|
59 | $head f$$ |
---|
60 | The object $icode f$$ has prototype |
---|
61 | $codei% |
---|
62 | ADFun<%Base%> %f% |
---|
63 | %$$ |
---|
64 | Note that the $cref ADFun$$ object $icode f$$ is not $code const$$ |
---|
65 | (see $cref/Uses Forward/sparse_hessian/Uses Forward/$$ below). |
---|
66 | |
---|
67 | $head x$$ |
---|
68 | The argument $icode x$$ has prototype |
---|
69 | $codei% |
---|
70 | const %VectorBase%& %x% |
---|
71 | %$$ |
---|
72 | (see $cref/VectorBase/sparse_hessian/VectorBase/$$ below) |
---|
73 | and its size |
---|
74 | must be equal to $icode n$$, the dimension of the |
---|
75 | $cref/domain/seq_property/Domain/$$ space for $icode f$$. |
---|
76 | It specifies |
---|
77 | that point at which to evaluate the Hessian. |
---|
78 | |
---|
79 | $head w$$ |
---|
80 | The argument $icode w$$ has prototype |
---|
81 | $codei% |
---|
82 | const %VectorBase%& %w% |
---|
83 | %$$ |
---|
84 | and size $latex m$$. |
---|
85 | It specifies the value of $latex w_i$$ in the expression |
---|
86 | for $icode hes$$. |
---|
87 | The more components of $latex w$$ that are identically zero, |
---|
88 | the more sparse the resulting Hessian may be (and hence the more efficient |
---|
89 | the calculation of $icode hes$$ may be). |
---|
90 | |
---|
91 | $head p$$ |
---|
92 | The argument $icode p$$ is optional and has prototype |
---|
93 | $codei% |
---|
94 | const %VectorSet%& %p% |
---|
95 | %$$ |
---|
96 | (see $cref/VectorSet/sparse_hessian/VectorSet/$$ below) |
---|
97 | If it has elements of type $code bool$$, |
---|
98 | its size is $latex n * n$$. |
---|
99 | If it has elements of type $code std::set<size_t>$$, |
---|
100 | its size is $latex n$$ and all its set elements are between |
---|
101 | zero and $latex n - 1$$. |
---|
102 | It specifies a |
---|
103 | $cref/sparsity pattern/glossary/Sparsity Pattern/$$ |
---|
104 | for the Hessian $latex H(x)$$. |
---|
105 | |
---|
106 | $subhead Purpose$$ |
---|
107 | If this sparsity pattern does not change between calls to |
---|
108 | $codei SparseHessian$$, it should be faster to calculate $icode p$$ once and |
---|
109 | pass this argument to $codei SparseHessian$$. |
---|
110 | If you specify $icode p$$, CppAD will use the same |
---|
111 | type of sparsity representation |
---|
112 | (vectors of $code bool$$ or vectors of $code std::set<size_t>$$) |
---|
113 | for its internal calculations. |
---|
114 | Otherwise, the representation |
---|
115 | for the internal calculations is unspecified. |
---|
116 | |
---|
117 | $subhead work$$ |
---|
118 | If you specify $icode work$$ in the calling sequence, |
---|
119 | it is not necessary to keep the sparsity pattern; see the heading |
---|
120 | $cref/p/sparse_hessian/work/p/$$ under the $icode work$$ description. |
---|
121 | |
---|
122 | $subhead Column Subset$$ |
---|
123 | If the arguments $icode row$$ and $icode col$$ are present, |
---|
124 | and $cref/color_method/sparse_hessian/work/color_method/$$ is |
---|
125 | $code cppad.general$$ or $code cppad.symmetric$$, |
---|
126 | it is not necessary to compute the entire sparsity pattern. |
---|
127 | Only the following subset of column values will matter: |
---|
128 | $codei% |
---|
129 | { %col%[%k%] : %k% = 0 , %...% , %K%-1 } |
---|
130 | %$$. |
---|
131 | |
---|
132 | |
---|
133 | $head row, col$$ |
---|
134 | The arguments $icode row$$ and $icode col$$ are optional and have prototype |
---|
135 | $codei% |
---|
136 | const %VectorSize%& %row% |
---|
137 | const %VectorSize%& %col% |
---|
138 | %$$ |
---|
139 | (see $cref/VectorSize/sparse_hessian/VectorSize/$$ below). |
---|
140 | They specify which rows and columns of $latex H (x)$$ are |
---|
141 | returned and in what order. |
---|
142 | We use $latex K$$ to denote the value $icode%hes%.size()%$$ |
---|
143 | which must also equal the size of $icode row$$ and $icode col$$. |
---|
144 | Furthermore, |
---|
145 | for $latex k = 0 , \ldots , K-1$$, it must hold that |
---|
146 | $latex row[k] < n$$ and $latex col[k] < n$$. |
---|
147 | In addition, |
---|
148 | all of the $latex (row[k], col[k])$$ pairs must correspond to a true value |
---|
149 | in the sparsity pattern $icode p$$. |
---|
150 | |
---|
151 | $head hes$$ |
---|
152 | The result $icode hes$$ has prototype |
---|
153 | $codei% |
---|
154 | %VectorBase% %hes% |
---|
155 | %$$ |
---|
156 | In the case where $icode row$$ and $icode col$$ are not present, |
---|
157 | the size of $icode hes$$ is $latex n * n$$ and |
---|
158 | its size is $latex n * n$$. |
---|
159 | In this case, for $latex i = 0 , \ldots , n - 1 $$ |
---|
160 | and $latex ell = 0 , \ldots , n - 1$$ |
---|
161 | $latex \[ |
---|
162 | hes [ j * n + \ell ] = \DD{ w^{\rm T} F }{ x_j }{ x_\ell } ( x ) |
---|
163 | \] $$ |
---|
164 | $pre |
---|
165 | |
---|
166 | $$ |
---|
167 | In the case where the arguments $icode row$$ and $icode col$$ are present, |
---|
168 | we use $latex K$$ to denote the size of $icode hes$$. |
---|
169 | The input value of its elements does not matter. |
---|
170 | Upon return, for $latex k = 0 , \ldots , K - 1$$, |
---|
171 | $latex \[ |
---|
172 | hes [ k ] = \DD{ w^{\rm T} F }{ x_j }{ x_\ell } (x) |
---|
173 | \; , \; |
---|
174 | \; {\rm where} \; |
---|
175 | j = row[k] |
---|
176 | \; {\rm and } \; |
---|
177 | \ell = col[k] |
---|
178 | \] $$ |
---|
179 | |
---|
180 | $head work$$ |
---|
181 | If this argument is present, it has prototype |
---|
182 | $codei% |
---|
183 | sparse_hessian_work& %work% |
---|
184 | %$$ |
---|
185 | This object can only be used with the routines $code SparseHessian$$. |
---|
186 | During its the first use, information is stored in $icode work$$. |
---|
187 | This is used to reduce the work done by future calls to $code SparseHessian$$ |
---|
188 | with the same $icode f$$, $icode p$$, $icode row$$, and $icode col$$. |
---|
189 | If a future call is made where any of these values have changed, |
---|
190 | you must first call $icode%work%.clear()%$$ |
---|
191 | to inform CppAD that this information needs to be recomputed. |
---|
192 | |
---|
193 | $subhead color_method$$ |
---|
194 | The coloring algorithm determines which rows and columns |
---|
195 | can be computed during the same sweep. |
---|
196 | This field has prototype |
---|
197 | $codei% |
---|
198 | std::string %work%.color_method |
---|
199 | %$$ |
---|
200 | This value only matters on the first call to $code sparse_hessian$$ that |
---|
201 | follows the $icode work$$ constructor or a call to |
---|
202 | $icode%work%.clear()%$$. |
---|
203 | $codei% |
---|
204 | |
---|
205 | "cppad.symmetric" |
---|
206 | %$$ |
---|
207 | This is the default coloring method (after a constructor or $code clear()$$). |
---|
208 | It takes advantage of the fact that the Hessian matrix |
---|
209 | is symmetric to find a coloring that requires fewer |
---|
210 | $cref/sweeps/sparse_hessian/n_sweep/$$. |
---|
211 | $codei% |
---|
212 | |
---|
213 | "cppad.general" |
---|
214 | %$$ |
---|
215 | This is the same as the $code "cppad"$$ method for the |
---|
216 | $cref/sparse_jacobian/sparse_jacobian/work/color_method/$$ calculation. |
---|
217 | $codei% |
---|
218 | |
---|
219 | "colpack.symmetric" |
---|
220 | %$$ |
---|
221 | This method requires that |
---|
222 | $cref colpack_prefix$$ was specified on the |
---|
223 | $cref/cmake command/cmake/CMake Command/$$ line. |
---|
224 | It also takes advantage of the fact that the Hessian matrix is symmetric. |
---|
225 | $codei% |
---|
226 | |
---|
227 | "colpack.general" |
---|
228 | %$$ |
---|
229 | This is the same as the $code "colpack"$$ method for the |
---|
230 | $cref/sparse_jacobian/sparse_jacobian/work/color_method/$$ calculation. |
---|
231 | |
---|
232 | $subhead colpack.star Deprecated 2017-06-01$$ |
---|
233 | The $code colpack.star$$ method is deprecated. |
---|
234 | It is the same as the $code colpack.symmetric$$ |
---|
235 | which should be used instead. |
---|
236 | |
---|
237 | $subhead p$$ |
---|
238 | If $icode work$$ is present, and it is not the first call after |
---|
239 | its construction or a clear, |
---|
240 | the sparsity pattern $icode p$$ is not used. |
---|
241 | This enables one to free the sparsity pattern |
---|
242 | and still compute corresponding sparse Hessians. |
---|
243 | |
---|
244 | $head n_sweep$$ |
---|
245 | The return value $icode n_sweep$$ has prototype |
---|
246 | $codei% |
---|
247 | size_t %n_sweep% |
---|
248 | %$$ |
---|
249 | It is the number of first order forward sweeps |
---|
250 | used to compute the requested Hessian values. |
---|
251 | Each first forward sweep is followed by a second order reverse sweep |
---|
252 | so it is also the number of reverse sweeps. |
---|
253 | This is proportional to the total work that $code SparseHessian$$ does, |
---|
254 | not counting the zero order forward sweep, |
---|
255 | or the work to combine multiple columns into a single |
---|
256 | forward-reverse sweep pair. |
---|
257 | |
---|
258 | $head VectorBase$$ |
---|
259 | The type $icode VectorBase$$ must be a $cref SimpleVector$$ class with |
---|
260 | $cref/elements of type/SimpleVector/Elements of Specified Type/$$ |
---|
261 | $icode Base$$. |
---|
262 | The routine $cref CheckSimpleVector$$ will generate an error message |
---|
263 | if this is not the case. |
---|
264 | |
---|
265 | $head VectorSet$$ |
---|
266 | The type $icode VectorSet$$ must be a $cref SimpleVector$$ class with |
---|
267 | $cref/elements of type/SimpleVector/Elements of Specified Type/$$ |
---|
268 | $code bool$$ or $code std::set<size_t>$$; |
---|
269 | see $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for a discussion |
---|
270 | of the difference. |
---|
271 | The routine $cref CheckSimpleVector$$ will generate an error message |
---|
272 | if this is not the case. |
---|
273 | |
---|
274 | $subhead Restrictions$$ |
---|
275 | If $icode VectorSet$$ has elements of $code std::set<size_t>$$, |
---|
276 | then $icode%p%[%i%]%$$ must return a reference (not a copy) to the |
---|
277 | corresponding set. |
---|
278 | According to section 26.3.2.3 of the 1998 C++ standard, |
---|
279 | $code std::valarray< std::set<size_t> >$$ does not satisfy |
---|
280 | this condition. |
---|
281 | |
---|
282 | $head VectorSize$$ |
---|
283 | The type $icode VectorSize$$ must be a $cref SimpleVector$$ class with |
---|
284 | $cref/elements of type/SimpleVector/Elements of Specified Type/$$ |
---|
285 | $code size_t$$. |
---|
286 | The routine $cref CheckSimpleVector$$ will generate an error message |
---|
287 | if this is not the case. |
---|
288 | |
---|
289 | $head Uses Forward$$ |
---|
290 | After each call to $cref Forward$$, |
---|
291 | the object $icode f$$ contains the corresponding |
---|
292 | $cref/Taylor coefficients/glossary/Taylor Coefficient/$$. |
---|
293 | After a call to any of the sparse Hessian routines, |
---|
294 | the zero order Taylor coefficients correspond to |
---|
295 | $icode%f%.Forward(0, %x%)%$$ |
---|
296 | and the other coefficients are unspecified. |
---|
297 | |
---|
298 | $children% |
---|
299 | example/sparse/sparse_hessian.cpp% |
---|
300 | example/sparse/sub_sparse_hes.cpp% |
---|
301 | example/sparse/sparse_sub_hes.cpp |
---|
302 | %$$ |
---|
303 | |
---|
304 | $head Example$$ |
---|
305 | The routine |
---|
306 | $cref sparse_hessian.cpp$$ |
---|
307 | is examples and tests of $code sparse_hessian$$. |
---|
308 | It return $code true$$, if it succeeds and $code false$$ otherwise. |
---|
309 | |
---|
310 | $head Subset Hessian$$ |
---|
311 | The routine |
---|
312 | $cref sub_sparse_hes.cpp$$ |
---|
313 | is an example and test that compute a sparse Hessian |
---|
314 | for a subset of the variables. |
---|
315 | It returns $code true$$, for success, and $code false$$ otherwise. |
---|
316 | |
---|
317 | $end |
---|
318 | ----------------------------------------------------------------------------- |
---|
319 | */ |
---|
320 | # include <cppad/local/std_set.hpp> |
---|
321 | # include <cppad/local/color_general.hpp> |
---|
322 | # include <cppad/local/color_symmetric.hpp> |
---|
323 | |
---|
324 | namespace CppAD { // BEGIN_CPPAD_NAMESPACE |
---|
325 | /*! |
---|
326 | \file sparse_hessian.hpp |
---|
327 | Sparse Hessian driver routine and helper functions. |
---|
328 | */ |
---|
329 | // =========================================================================== |
---|
330 | /*! |
---|
331 | class used by SparseHessian to hold information |
---|
332 | so it does not need to be recomputed. |
---|
333 | */ |
---|
334 | class sparse_hessian_work { |
---|
335 | public: |
---|
336 | /// Coloring method: "cppad", or "colpack" |
---|
337 | /// (this field is set by user) |
---|
338 | std::string color_method; |
---|
339 | /// row and column indicies for return values |
---|
340 | /// (some may be reflected by star coloring algorithm) |
---|
341 | CppAD::vector<size_t> row; |
---|
342 | CppAD::vector<size_t> col; |
---|
343 | /// indices that sort the user row and col arrays by color |
---|
344 | CppAD::vector<size_t> order; |
---|
345 | /// results of the coloring algorithm |
---|
346 | CppAD::vector<size_t> color; |
---|
347 | |
---|
348 | /// constructor |
---|
349 | sparse_hessian_work(void) : color_method("cppad.symmetric") |
---|
350 | { } |
---|
351 | /// inform CppAD that this information needs to be recomputed |
---|
352 | void clear(void) |
---|
353 | { color_method = "cppad.symmetric"; |
---|
354 | row.clear(); |
---|
355 | col.clear(); |
---|
356 | order.clear(); |
---|
357 | color.clear(); |
---|
358 | } |
---|
359 | }; |
---|
360 | // =========================================================================== |
---|
361 | /*! |
---|
362 | Private helper function that does computation for all Sparse Hessian cases. |
---|
363 | |
---|
364 | \tparam Base |
---|
365 | is the base type for the recording that is stored in this ADFun<Base object. |
---|
366 | |
---|
367 | \tparam VectorBase |
---|
368 | is a simple vector class with elements of type \a Base. |
---|
369 | |
---|
370 | \tparam VectorSet |
---|
371 | is a simple vector class with elements of type |
---|
372 | \c bool or \c std::set<size_t>. |
---|
373 | |
---|
374 | \tparam VectorSize |
---|
375 | is sparse_pack or sparse_list. |
---|
376 | |
---|
377 | \param x [in] |
---|
378 | is a vector specifing the point at which to compute the Hessian. |
---|
379 | |
---|
380 | \param w [in] |
---|
381 | is the weighting vector that defines a scalar valued function by |
---|
382 | a weighted sum of the components of the vector valued function |
---|
383 | $latex F(x)$$. |
---|
384 | |
---|
385 | \param sparsity [in] |
---|
386 | is the sparsity pattern for the Hessian that we are calculating. |
---|
387 | |
---|
388 | \param user_row [in] |
---|
389 | is the vector of row indices for the returned Hessian values. |
---|
390 | |
---|
391 | \param user_col [in] |
---|
392 | is the vector of columns indices for the returned Hessian values. |
---|
393 | It must have the same size as user_row. |
---|
394 | |
---|
395 | \param hes [out] |
---|
396 | is the vector of Hessian values. |
---|
397 | It must have the same size as user_row. |
---|
398 | The return value <code>hes[k]</code> is the second partial of |
---|
399 | \f$ w^{\rm T} F(x)\f$ with respect to the |
---|
400 | <code>row[k]</code> and <code>col[k]</code> component of \f$ x\f$. |
---|
401 | |
---|
402 | \param work |
---|
403 | This structure contains information that is computed by \c SparseHessianCompute. |
---|
404 | If the sparsity pattern, \c row vector, or \c col vectors |
---|
405 | are not the same between calls to \c SparseHessianCompute, |
---|
406 | \c work.clear() must be called to reinitialize \c work. |
---|
407 | |
---|
408 | \return |
---|
409 | Is the number of first order forward sweeps used to compute the |
---|
410 | requested Hessian values. |
---|
411 | (This is also equal to the number of second order reverse sweeps.) |
---|
412 | The total work, not counting the zero order |
---|
413 | forward sweep, or the time to combine computations, is proportional to this |
---|
414 | return value. |
---|
415 | */ |
---|
416 | template<class Base> |
---|
417 | template <class VectorBase, class VectorSet, class VectorSize> |
---|
418 | size_t ADFun<Base>::SparseHessianCompute( |
---|
419 | const VectorBase& x , |
---|
420 | const VectorBase& w , |
---|
421 | VectorSet& sparsity , |
---|
422 | const VectorSize& user_row , |
---|
423 | const VectorSize& user_col , |
---|
424 | VectorBase& hes , |
---|
425 | sparse_hessian_work& work ) |
---|
426 | { |
---|
427 | using CppAD::vectorBool; |
---|
428 | size_t i, k, ell; |
---|
429 | |
---|
430 | CppAD::vector<size_t>& row(work.row); |
---|
431 | CppAD::vector<size_t>& col(work.col); |
---|
432 | CppAD::vector<size_t>& color(work.color); |
---|
433 | CppAD::vector<size_t>& order(work.order); |
---|
434 | |
---|
435 | size_t n = Domain(); |
---|
436 | |
---|
437 | // some values |
---|
438 | const Base zero(0); |
---|
439 | const Base one(1); |
---|
440 | |
---|
441 | // check VectorBase is Simple Vector class with Base type elements |
---|
442 | CheckSimpleVector<Base, VectorBase>(); |
---|
443 | |
---|
444 | // number of components of Hessian that are required |
---|
445 | size_t K = hes.size(); |
---|
446 | CPPAD_ASSERT_UNKNOWN( size_t( user_row.size() ) == K ); |
---|
447 | CPPAD_ASSERT_UNKNOWN( size_t( user_col.size() ) == K ); |
---|
448 | |
---|
449 | CPPAD_ASSERT_UNKNOWN( size_t(x.size()) == n ); |
---|
450 | CPPAD_ASSERT_UNKNOWN( color.size() == 0 || color.size() == n ); |
---|
451 | CPPAD_ASSERT_UNKNOWN( row.size() == 0 || row.size() == K ); |
---|
452 | CPPAD_ASSERT_UNKNOWN( col.size() == 0 || col.size() == K ); |
---|
453 | |
---|
454 | |
---|
455 | // Point at which we are evaluating the Hessian |
---|
456 | Forward(0, x); |
---|
457 | |
---|
458 | // check for case where nothing (except Forward above) to do |
---|
459 | if( K == 0 ) |
---|
460 | return 0; |
---|
461 | |
---|
462 | // Rows of the Hessian (i below) correspond to the forward mode index |
---|
463 | // and columns (j below) correspond to the reverse mode index. |
---|
464 | if( color.size() == 0 ) |
---|
465 | { |
---|
466 | CPPAD_ASSERT_UNKNOWN( sparsity.n_set() == n ); |
---|
467 | CPPAD_ASSERT_UNKNOWN( sparsity.end() == n ); |
---|
468 | |
---|
469 | // copy user rwo and col to work space |
---|
470 | row.resize(K); |
---|
471 | col.resize(K); |
---|
472 | for(k = 0; k < K; k++) |
---|
473 | { row[k] = user_row[k]; |
---|
474 | col[k] = user_col[k]; |
---|
475 | } |
---|
476 | |
---|
477 | // execute coloring algorithm |
---|
478 | color.resize(n); |
---|
479 | if( work.color_method == "cppad.general" ) |
---|
480 | local::color_general_cppad(sparsity, row, col, color); |
---|
481 | else if( work.color_method == "cppad.symmetric" ) |
---|
482 | local::color_symmetric_cppad(sparsity, row, col, color); |
---|
483 | else if( work.color_method == "colpack.general" ) |
---|
484 | { |
---|
485 | # if CPPAD_HAS_COLPACK |
---|
486 | local::color_general_colpack(sparsity, row, col, color); |
---|
487 | # else |
---|
488 | CPPAD_ASSERT_KNOWN( |
---|
489 | false, |
---|
490 | "SparseHessian: work.color_method = colpack.general " |
---|
491 | "and colpack_prefix missing from cmake command line." |
---|
492 | ); |
---|
493 | # endif |
---|
494 | } |
---|
495 | else if( |
---|
496 | work.color_method == "colpack.symmetric" || |
---|
497 | work.color_method == "colpack.star" |
---|
498 | ) |
---|
499 | { |
---|
500 | # if CPPAD_HAS_COLPACK |
---|
501 | local::color_symmetric_colpack(sparsity, row, col, color); |
---|
502 | # else |
---|
503 | CPPAD_ASSERT_KNOWN( |
---|
504 | false, |
---|
505 | "SparseHessian: work.color_method is " |
---|
506 | "colpack.symmetric or colpack.star\n" |
---|
507 | "and colpack_prefix missing from cmake command line." |
---|
508 | ); |
---|
509 | # endif |
---|
510 | } |
---|
511 | else |
---|
512 | { CPPAD_ASSERT_KNOWN( |
---|
513 | false, |
---|
514 | "SparseHessian: work.color_method is not valid." |
---|
515 | ); |
---|
516 | } |
---|
517 | |
---|
518 | // put sorting indices in color order |
---|
519 | VectorSize key(K); |
---|
520 | order.resize(K); |
---|
521 | for(k = 0; k < K; k++) |
---|
522 | key[k] = color[ row[k] ]; |
---|
523 | index_sort(key, order); |
---|
524 | |
---|
525 | } |
---|
526 | size_t n_color = 1; |
---|
527 | for(ell = 0; ell < n; ell++) if( color[ell] < n ) |
---|
528 | n_color = std::max(n_color, color[ell] + 1); |
---|
529 | |
---|
530 | // direction vector for calls to forward (rows of the Hessian) |
---|
531 | VectorBase u(n); |
---|
532 | |
---|
533 | // location for return values from reverse (columns of the Hessian) |
---|
534 | VectorBase ddw(2 * n); |
---|
535 | |
---|
536 | // initialize the return value |
---|
537 | for(k = 0; k < K; k++) |
---|
538 | hes[k] = zero; |
---|
539 | |
---|
540 | // loop over colors |
---|
541 | # ifndef NDEBUG |
---|
542 | const std::string& coloring = work.color_method; |
---|
543 | # endif |
---|
544 | k = 0; |
---|
545 | for(ell = 0; ell < n_color; ell++) |
---|
546 | if( k == K ) |
---|
547 | { // kludge because colpack returns colors that are not used |
---|
548 | // (it does not know about the subset corresponding to row, col) |
---|
549 | CPPAD_ASSERT_UNKNOWN( |
---|
550 | coloring == "colpack.general" || |
---|
551 | coloring == "colpack.symmetic" || |
---|
552 | coloring == "colpack.star" |
---|
553 | ); |
---|
554 | } |
---|
555 | else if( color[ row[ order[k] ] ] != ell ) |
---|
556 | { // kludge because colpack returns colors that are not used |
---|
557 | // (it does not know about the subset corresponding to row, col) |
---|
558 | CPPAD_ASSERT_UNKNOWN( |
---|
559 | coloring == "colpack.general" || |
---|
560 | coloring == "colpack.symmetic" || |
---|
561 | coloring == "colpack.star" |
---|
562 | ); |
---|
563 | } |
---|
564 | else |
---|
565 | { CPPAD_ASSERT_UNKNOWN( color[ row[ order[k] ] ] == ell ); |
---|
566 | |
---|
567 | // combine all rows with this color |
---|
568 | for(i = 0; i < n; i++) |
---|
569 | { u[i] = zero; |
---|
570 | if( color[i] == ell ) |
---|
571 | u[i] = one; |
---|
572 | } |
---|
573 | // call forward mode for all these rows at once |
---|
574 | Forward(1, u); |
---|
575 | |
---|
576 | // evaluate derivative of w^T * F'(x) * u |
---|
577 | ddw = Reverse(2, w); |
---|
578 | |
---|
579 | // set the corresponding components of the result |
---|
580 | while( k < K && color[ row[ order[k] ] ] == ell ) |
---|
581 | { hes[ order[k] ] = ddw[ col[ order[k] ] * 2 + 1 ]; |
---|
582 | k++; |
---|
583 | } |
---|
584 | } |
---|
585 | return n_color; |
---|
586 | } |
---|
587 | // =========================================================================== |
---|
588 | // Public Member Functions |
---|
589 | // =========================================================================== |
---|
590 | /*! |
---|
591 | Compute user specified subset of a sparse Hessian. |
---|
592 | |
---|
593 | The C++ source code corresponding to this operation is |
---|
594 | \verbatim |
---|
595 | SparceHessian(x, w, p, row, col, hes, work) |
---|
596 | \endverbatim |
---|
597 | |
---|
598 | \tparam Base |
---|
599 | is the base type for the recording that is stored in this ADFun<Base object. |
---|
600 | |
---|
601 | \tparam VectorBase |
---|
602 | is a simple vector class with elements of type \a Base. |
---|
603 | |
---|
604 | \tparam VectorSet |
---|
605 | is a simple vector class with elements of type |
---|
606 | \c bool or \c std::set<size_t>. |
---|
607 | |
---|
608 | \tparam VectorSize |
---|
609 | is a simple vector class with elements of type \c size_t. |
---|
610 | |
---|
611 | \param x [in] |
---|
612 | is a vector specifing the point at which to compute the Hessian. |
---|
613 | |
---|
614 | \param w [in] |
---|
615 | is the weighting vector that defines a scalar valued function by |
---|
616 | a weighted sum of the components of the vector valued function |
---|
617 | $latex F(x)$$. |
---|
618 | |
---|
619 | \param p [in] |
---|
620 | is the sparsity pattern for the Hessian that we are calculating. |
---|
621 | |
---|
622 | \param row [in] |
---|
623 | is the vector of row indices for the returned Hessian values. |
---|
624 | |
---|
625 | \param col [in] |
---|
626 | is the vector of columns indices for the returned Hessian values. |
---|
627 | It must have the same size are r. |
---|
628 | |
---|
629 | \param hes [out] |
---|
630 | is the vector of Hessian values. |
---|
631 | It must have the same size are r. |
---|
632 | The return value <code>hes[k]</code> is the second partial of |
---|
633 | \f$ w^{\rm T} F(x)\f$ with respect to the |
---|
634 | <code>row[k]</code> and <code>col[k]</code> component of \f$ x\f$. |
---|
635 | |
---|
636 | \param work |
---|
637 | This structure contains information that is computed by \c SparseHessianCompute. |
---|
638 | If the sparsity pattern, \c row vector, or \c col vectors |
---|
639 | are not the same between calls to \c SparseHessian, |
---|
640 | \c work.clear() must be called to reinitialize \c work. |
---|
641 | |
---|
642 | \return |
---|
643 | Is the number of first order forward sweeps used to compute the |
---|
644 | requested Hessian values. |
---|
645 | (This is also equal to the number of second order reverse sweeps.) |
---|
646 | The total work, not counting the zero order |
---|
647 | forward sweep, or the time to combine computations, is proportional to this |
---|
648 | return value. |
---|
649 | */ |
---|
650 | template<class Base> |
---|
651 | template <class VectorBase, class VectorSet, class VectorSize> |
---|
652 | size_t ADFun<Base>::SparseHessian( |
---|
653 | const VectorBase& x , |
---|
654 | const VectorBase& w , |
---|
655 | const VectorSet& p , |
---|
656 | const VectorSize& row , |
---|
657 | const VectorSize& col , |
---|
658 | VectorBase& hes , |
---|
659 | sparse_hessian_work& work ) |
---|
660 | { |
---|
661 | size_t n = Domain(); |
---|
662 | size_t K = hes.size(); |
---|
663 | # ifndef NDEBUG |
---|
664 | size_t k; |
---|
665 | CPPAD_ASSERT_KNOWN( |
---|
666 | size_t(x.size()) == n , |
---|
667 | "SparseHessian: size of x not equal domain dimension for f." |
---|
668 | ); |
---|
669 | CPPAD_ASSERT_KNOWN( |
---|
670 | size_t(row.size()) == K && size_t(col.size()) == K , |
---|
671 | "SparseHessian: either r or c does not have the same size as ehs." |
---|
672 | ); |
---|
673 | CPPAD_ASSERT_KNOWN( |
---|
674 | work.color.size() == 0 || work.color.size() == n, |
---|
675 | "SparseHessian: invalid value in work." |
---|
676 | ); |
---|
677 | for(k = 0; k < K; k++) |
---|
678 | { CPPAD_ASSERT_KNOWN( |
---|
679 | row[k] < n, |
---|
680 | "SparseHessian: invalid value in r." |
---|
681 | ); |
---|
682 | CPPAD_ASSERT_KNOWN( |
---|
683 | col[k] < n, |
---|
684 | "SparseHessian: invalid value in c." |
---|
685 | ); |
---|
686 | } |
---|
687 | if( work.color.size() != 0 ) |
---|
688 | for(size_t j = 0; j < n; j++) CPPAD_ASSERT_KNOWN( |
---|
689 | work.color[j] <= n, |
---|
690 | "SparseHessian: invalid value in work." |
---|
691 | ); |
---|
692 | # endif |
---|
693 | // check for case where there is nothing to compute |
---|
694 | size_t n_sweep = 0; |
---|
695 | if( K == 0 ) |
---|
696 | return n_sweep; |
---|
697 | |
---|
698 | typedef typename VectorSet::value_type Set_type; |
---|
699 | typedef typename local::internal_sparsity<Set_type>::pattern_type Pattern_type; |
---|
700 | Pattern_type s; |
---|
701 | if( work.color.size() == 0 ) |
---|
702 | { bool transpose = false; |
---|
703 | const char* error_msg = "SparseHessian: sparsity pattern" |
---|
704 | " does not have proper row or column dimension"; |
---|
705 | sparsity_user2internal(s, p, n, n, transpose, error_msg); |
---|
706 | } |
---|
707 | n_sweep = SparseHessianCompute(x, w, s, row, col, hes, work); |
---|
708 | return n_sweep; |
---|
709 | } |
---|
710 | /*! |
---|
711 | Compute a sparse Hessian. |
---|
712 | |
---|
713 | The C++ source code coresponding to this operation is |
---|
714 | \verbatim |
---|
715 | hes = SparseHessian(x, w, p) |
---|
716 | \endverbatim |
---|
717 | |
---|
718 | |
---|
719 | \tparam Base |
---|
720 | is the base type for the recording that is stored in this |
---|
721 | ADFun<Base object. |
---|
722 | |
---|
723 | \tparam VectorBase |
---|
724 | is a simple vector class with elements of the \a Base. |
---|
725 | |
---|
726 | \tparam VectorSet |
---|
727 | is a simple vector class with elements of type |
---|
728 | \c bool or \c std::set<size_t>. |
---|
729 | |
---|
730 | \param x [in] |
---|
731 | is a vector specifing the point at which to compute the Hessian. |
---|
732 | |
---|
733 | \param w [in] |
---|
734 | The Hessian is computed for a weighted sum of the components |
---|
735 | of the function corresponding to this ADFun<Base> object. |
---|
736 | The argument \a w specifies the weights for each component. |
---|
737 | It must have size equal to the range dimension for this ADFun<Base> object. |
---|
738 | |
---|
739 | \param p [in] |
---|
740 | is a sparsity pattern for the Hessian. |
---|
741 | |
---|
742 | \return |
---|
743 | Will be a vector of size \c n * n containing the Hessian of |
---|
744 | at the point specified by \a x |
---|
745 | (where \c n is the domain dimension for this ADFun<Base> object). |
---|
746 | */ |
---|
747 | template <class Base> |
---|
748 | template <class VectorBase, class VectorSet> |
---|
749 | VectorBase ADFun<Base>::SparseHessian( |
---|
750 | const VectorBase& x, const VectorBase& w, const VectorSet& p |
---|
751 | ) |
---|
752 | { size_t i, j, k; |
---|
753 | |
---|
754 | size_t n = Domain(); |
---|
755 | VectorBase hes(n * n); |
---|
756 | |
---|
757 | CPPAD_ASSERT_KNOWN( |
---|
758 | size_t(x.size()) == n, |
---|
759 | "SparseHessian: size of x not equal domain size for f." |
---|
760 | ); |
---|
761 | |
---|
762 | typedef typename VectorSet::value_type Set_type; |
---|
763 | typedef typename local::internal_sparsity<Set_type>::pattern_type Pattern_type; |
---|
764 | |
---|
765 | // initialize the return value as zero |
---|
766 | Base zero(0); |
---|
767 | for(i = 0; i < n; i++) |
---|
768 | for(j = 0; j < n; j++) |
---|
769 | hes[i * n + j] = zero; |
---|
770 | |
---|
771 | // arguments to SparseHessianCompute |
---|
772 | Pattern_type s; |
---|
773 | CppAD::vector<size_t> row; |
---|
774 | CppAD::vector<size_t> col; |
---|
775 | sparse_hessian_work work; |
---|
776 | bool transpose = false; |
---|
777 | const char* error_msg = "SparseHessian: sparsity pattern" |
---|
778 | " does not have proper row or column dimension"; |
---|
779 | sparsity_user2internal(s, p, n, n, transpose, error_msg); |
---|
780 | k = 0; |
---|
781 | for(i = 0; i < n; i++) |
---|
782 | { typename Pattern_type::const_iterator itr(s, i); |
---|
783 | j = *itr; |
---|
784 | while( j != s.end() ) |
---|
785 | { row.push_back(i); |
---|
786 | col.push_back(j); |
---|
787 | k++; |
---|
788 | j = *(++itr); |
---|
789 | } |
---|
790 | } |
---|
791 | size_t K = k; |
---|
792 | VectorBase H(K); |
---|
793 | |
---|
794 | // now we have folded this into the following case |
---|
795 | SparseHessianCompute(x, w, s, row, col, H, work); |
---|
796 | |
---|
797 | // now set the non-zero return values |
---|
798 | for(k = 0; k < K; k++) |
---|
799 | hes[ row[k] * n + col[k] ] = H[k]; |
---|
800 | |
---|
801 | return hes; |
---|
802 | } |
---|
803 | /*! |
---|
804 | Compute a sparse Hessian |
---|
805 | |
---|
806 | The C++ source code coresponding to this operation is |
---|
807 | \verbatim |
---|
808 | hes = SparseHessian(x, w) |
---|
809 | \endverbatim |
---|
810 | |
---|
811 | |
---|
812 | \tparam Base |
---|
813 | is the base type for the recording that is stored in this |
---|
814 | ADFun<Base object. |
---|
815 | |
---|
816 | \tparam VectorBase |
---|
817 | is a simple vector class with elements of the \a Base. |
---|
818 | |
---|
819 | \param x [in] |
---|
820 | is a vector specifing the point at which to compute the Hessian. |
---|
821 | |
---|
822 | \param w [in] |
---|
823 | The Hessian is computed for a weighted sum of the components |
---|
824 | of the function corresponding to this ADFun<Base> object. |
---|
825 | The argument \a w specifies the weights for each component. |
---|
826 | It must have size equal to the range dimension for this ADFun<Base> object. |
---|
827 | |
---|
828 | \return |
---|
829 | Will be a vector of size \c n * n containing the Hessian of |
---|
830 | at the point specified by \a x |
---|
831 | (where \c n is the domain dimension for this ADFun<Base> object). |
---|
832 | */ |
---|
833 | template <class Base> |
---|
834 | template <class VectorBase> |
---|
835 | VectorBase ADFun<Base>::SparseHessian(const VectorBase &x, const VectorBase &w) |
---|
836 | { size_t i, j, k; |
---|
837 | typedef CppAD::vectorBool VectorBool; |
---|
838 | |
---|
839 | size_t m = Range(); |
---|
840 | size_t n = Domain(); |
---|
841 | |
---|
842 | // determine the sparsity pattern p for Hessian of w^T F |
---|
843 | VectorBool r(n * n); |
---|
844 | for(j = 0; j < n; j++) |
---|
845 | { for(k = 0; k < n; k++) |
---|
846 | r[j * n + k] = false; |
---|
847 | r[j * n + j] = true; |
---|
848 | } |
---|
849 | ForSparseJac(n, r); |
---|
850 | // |
---|
851 | VectorBool s(m); |
---|
852 | for(i = 0; i < m; i++) |
---|
853 | s[i] = w[i] != 0; |
---|
854 | VectorBool p = RevSparseHes(n, s); |
---|
855 | |
---|
856 | // compute sparse Hessian |
---|
857 | return SparseHessian(x, w, p); |
---|
858 | } |
---|
859 | |
---|
860 | } // END_CPPAD_NAMESPACE |
---|
861 | # endif |
---|