1 # ifndef CPPAD_CORE_FOR_SPARSE_HES_HPP
2 # define CPPAD_CORE_FOR_SPARSE_HES_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 ForSparseHes$$
17 $spell
18 Andrea Walther
19 std
20 VecAD
21 Jacobian
22 Jac
23 Hessian
24 Hes
25 const
26 Bool
27 Dep
28 proportional
29 var
30 cpp
31 $$
32
33 $section Hessian Sparsity Pattern: Forward Mode$$
34
35 $head Syntax$$
36 $icode%h% = %f%.ForSparseHes(%r%, %s%)
37 %$$
38
39 $head Purpose$$
40 We use $latex F : \B{R}^n \rightarrow \B{R}^m$$ to denote the
41 $cref/AD function/glossary/AD Function/$$ corresponding to $icode f$$.
42 we define
43 $latex \[
44 \begin{array}{rcl}
45 H(x)
46 & = & \partial_x \left[ \partial_u S \cdot F[ x + R \cdot u ] \right]_{u=0}
47 \\
48 & = & R^\R{T} \cdot (S \cdot F)^{(2)} ( x ) \cdot R
49 \end{array}
50 \] $$
51 Where $latex R \in \B{R}^{n \times n}$$ is a diagonal matrix
52 and $latex S \in \B{R}^{1 \times m}$$ is a row vector.
53 Given a
54 $cref/sparsity pattern/glossary/Sparsity Pattern/$$
55 for the diagonal of $latex R$$ and the vector $latex S$$,
56 $code ForSparseHes$$ returns a sparsity pattern for the $latex H(x)$$.
57
58 $head f$$
59 The object $icode f$$ has prototype
60 $codei%
61 const ADFun<%Base%> %f%
62 %$$
63
64 $head x$$
65 If the operation sequence in $icode f$$ is
66 $cref/independent/glossary/Operation/Independent/$$ of
67 the independent variables in $latex x \in B^n$$,
68 the sparsity pattern is valid for all values of
69 (even if it has $cref CondExp$$ or $cref VecAD$$ operations).
70
71 $head r$$
72 The argument $icode r$$ has prototype
73 $codei%
74 const %VectorSet%& %r%
75 %$$
76 (see $cref/VectorSet/ForSparseHes/VectorSet/$$ below)
77 If it has elements of type $code bool$$,
78 its size is $latex n$$.
79 If it has elements of type $code std::set<size_t>$$,
80 its size is one and all the elements of $icode%s%[0]%$$
81 are between zero and $latex n - 1$$.
82 It specifies a
83 $cref/sparsity pattern/glossary/Sparsity Pattern/$$
84 for the diagonal of $latex R$$.
85 The fewer non-zero elements in this sparsity pattern,
86 the faster the calculation should be and the more sparse
87 $latex H(x)$$ should be.
88
89 $head s$$
90 The argument $icode s$$ has prototype
91 $codei%
92 const %VectorSet%& %s%
93 %$$
94 (see $cref/VectorSet/ForSparseHes/VectorSet/$$ below)
95 If it has elements of type $code bool$$,
96 its size is $latex m$$.
97 If it has elements of type $code std::set<size_t>$$,
98 its size is one and all the elements of $icode%s%[0]%$$
99 are between zero and $latex m - 1$$.
100 It specifies a
101 $cref/sparsity pattern/glossary/Sparsity Pattern/$$
102 for the vector $icode S$$.
103 The fewer non-zero elements in this sparsity pattern,
104 the faster the calculation should be and the more sparse
105 $latex H(x)$$ should be.
106
107 $head h$$
108 The result $icode h$$ has prototype
109 $codei%
110 %VectorSet%& %h%
111 %$$
112 (see $cref/VectorSet/ForSparseHes/VectorSet/$$ below).
113 If $icode h$$ has elements of type $code bool$$,
114 its size is $latex n * n$$.
115 If it has elements of type $code std::set<size_t>$$,
116 its size is $latex n$$ and all the set elements are between
117 zero and $icode%n%-1%$$ inclusive.
118 It specifies a
119 $cref/sparsity pattern/glossary/Sparsity Pattern/$$
120 for the matrix $latex H(x)$$.
121
122 $head VectorSet$$
123 The type $icode VectorSet$$ must be a $cref SimpleVector$$ class with
124 $cref/elements of type/SimpleVector/Elements of Specified Type/$$
125 $code bool$$ or $code std::set<size_t>$$;
126 see $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for a discussion
127 of the difference.
128 The type of the elements of
129 $cref/VectorSet/ForSparseHes/VectorSet/$$ must be the
130 same as the type of the elements of $icode r$$.
131
132 $head Algorithm$$
133 See Algorithm II in
134 $italic Computing sparse Hessians with automatic differentiation$$
135 by Andrea Walther.
136 Note that $icode s$$ provides the information so that
137 'dead ends' are not included in the sparsity pattern.
138
139 $head Example$$
140 $children%
141 example/sparse/for_sparse_hes.cpp
142 %$$
143 The file
144 $cref for_sparse_hes.cpp$$
145 contains an example and test of this operation.
146 It returns true if it succeeds and false otherwise.
147
148 $end
149 -----------------------------------------------------------------------------
150 */
151 # include <algorithm>
152 # include <cppad/local/pod_vector.hpp>
153 # include <cppad/local/std_set.hpp>
154
155 namespace CppAD { // BEGIN_CPPAD_NAMESPACE
156 /*!
157 \file for_sparse_hes.hpp
158 Forward mode Hessian sparsity patterns.
159 */
160 // ===========================================================================
161 // ForSparseHesCase
162 /*!
163 Private helper function for ForSparseHes(q, s) bool sparsity.
164
165 All of the description in the public member function ForSparseHes(q, s)
166 applies.
167
168 \param set_type
169 is a \c bool value. This argument is used to dispatch to the proper source
170 code depending on the vlaue of \c VectorSet::value_type.
171
172 \param r
173 See \c ForSparseHes(r, s).
174
175 \param s
176 See \c ForSparseHes(r, s).
177
178 \param h
179 is the return value for the corresponging call to \c ForSparseJac(q, s).
180 */
181 template <class Base>
182 template <class VectorSet>
ForSparseHesCase(bool set_type,const VectorSet & r,const VectorSet & s,VectorSet & h)183 void ADFun<Base>::ForSparseHesCase(
184 bool set_type ,
185 const VectorSet& r ,
186 const VectorSet& s ,
187 VectorSet& h )
188 { size_t n = Domain();
189 size_t m = Range();
190 //
191 // check Vector is Simple VectorSet class with bool elements
192 CheckSimpleVector<bool, VectorSet>();
193 //
194 CPPAD_ASSERT_KNOWN(
195 size_t(r.size()) == n,
196 "ForSparseHes: size of r is not equal to\n"
197 "domain dimension for ADFun object."
198 );
199 CPPAD_ASSERT_KNOWN(
200 size_t(s.size()) == m,
201 "ForSparseHes: size of s is not equal to\n"
202 "range dimension for ADFun object."
203 );
204 //
205 // sparsity pattern corresponding to r
206 local::sparse_pack for_jac_pattern;
207 for_jac_pattern.resize(num_var_tape_, n + 1);
208 for(size_t i = 0; i < n; i++)
209 { CPPAD_ASSERT_UNKNOWN( ind_taddr_[i] < n + 1 );
210 // ind_taddr_[i] is operator taddr for i-th independent variable
211 CPPAD_ASSERT_UNKNOWN( play_.GetOp( ind_taddr_[i] ) == local::InvOp );
212 //
213 // Use add_element when only adding one element per set is added.
214 if( r[i] )
215 for_jac_pattern.add_element( ind_taddr_[i], ind_taddr_[i] );
216 }
217 // compute forward Jacobiain sparsity pattern
218 bool dependency = false;
219 local::for_jac_sweep(
220 &play_,
221 dependency,
222 n,
223 num_var_tape_,
224 for_jac_pattern
225 );
226 // sparsity pattern correspnding to s
227 local::sparse_pack rev_jac_pattern;
228 rev_jac_pattern.resize(num_var_tape_, 1);
229 for(size_t i = 0; i < m; i++)
230 { CPPAD_ASSERT_UNKNOWN( dep_taddr_[i] < num_var_tape_ );
231 //
232 // Use add_element when only adding one element per set is added.
233 if( s[i] )
234 rev_jac_pattern.add_element( dep_taddr_[i], 0);
235 }
236 // compute reverse sparsity pattern for dependency analysis
237 // (note that we are only want non-zero derivatives not true dependency)
238 local::rev_jac_sweep(
239 &play_,
240 dependency,
241 n,
242 num_var_tape_,
243 rev_jac_pattern
244 );
245 // vector of sets that will hold the forward Hessain values
246 local::sparse_pack for_hes_pattern;
247 for_hes_pattern.resize(n+1, n+1);
248 //
249 // compute the Hessian sparsity patterns
250 local::for_hes_sweep(
251 &play_,
252 n,
253 num_var_tape_,
254 for_jac_pattern,
255 rev_jac_pattern,
256 for_hes_pattern
257 );
258 // initialize return values corresponding to independent variables
259 h.resize(n * n);
260 for(size_t i = 0; i < n; i++)
261 { for(size_t j = 0; j < n; j++)
262 h[ i * n + j ] = false;
263 }
264 // copy to result pattern
265 CPPAD_ASSERT_UNKNOWN( for_hes_pattern.end() == n+1 );
266 for(size_t i = 0; i < n; i++)
267 { // ind_taddr_[i] is operator taddr for i-th independent variable
268 CPPAD_ASSERT_UNKNOWN( ind_taddr_[i] == i + 1 );
269 CPPAD_ASSERT_UNKNOWN( play_.GetOp( ind_taddr_[i] ) == local::InvOp );
270
271 // extract the result from for_hes_pattern
272 local::sparse_pack::const_iterator itr(for_hes_pattern, ind_taddr_[i] );
273 size_t j = *itr;
274 while( j < for_hes_pattern.end() )
275 { CPPAD_ASSERT_UNKNOWN( 0 < j )
276 h[ i * n + (j-1) ] = true;
277 j = *(++itr);
278 }
279 }
280 }
281 /*!
282 Private helper function for ForSparseHes(q, s) set sparsity.
283
284 All of the description in the public member function ForSparseHes(q, s)
285 applies.
286
287 \param set_type
288 is a \c std::set<size_t> value.
289 This argument is used to dispatch to the proper source
290 code depending on the vlaue of \c VectorSet::value_type.
291
292 \param r
293 See \c ForSparseHes(r, s).
294
295 \param s
296 See \c ForSparseHes(q, s).
297
298 \param h
299 is the return value for the corresponging call to \c ForSparseJac(q, s).
300 */
301 template <class Base>
302 template <class VectorSet>
ForSparseHesCase(const std::set<size_t> & set_type,const VectorSet & r,const VectorSet & s,VectorSet & h)303 void ADFun<Base>::ForSparseHesCase(
304 const std::set<size_t>& set_type ,
305 const VectorSet& r ,
306 const VectorSet& s ,
307 VectorSet& h )
308 { size_t n = Domain();
309 # ifndef NDEBUG
310 size_t m = Range();
311 # endif
312 std::set<size_t>::const_iterator itr_1;
313 //
314 // check VectorSet is Simple Vector class with sets for elements
315 CheckSimpleVector<std::set<size_t>, VectorSet>(
316 local::one_element_std_set<size_t>(), local::two_element_std_set<size_t>()
317 );
318 CPPAD_ASSERT_KNOWN(
319 r.size() == 1,
320 "ForSparseHes: size of s is not equal to one."
321 );
322 CPPAD_ASSERT_KNOWN(
323 s.size() == 1,
324 "ForSparseHes: size of s is not equal to one."
325 );
326 //
327 // sparsity pattern corresponding to r
328 local::sparse_list for_jac_pattern;
329 for_jac_pattern.resize(num_var_tape_, n + 1);
330 itr_1 = r[0].begin();
331 while( itr_1 != r[0].end() )
332 { size_t i = *itr_1++;
333 CPPAD_ASSERT_UNKNOWN( ind_taddr_[i] < n + 1 );
334 // ind_taddr_[i] is operator taddr for i-th independent variable
335 CPPAD_ASSERT_UNKNOWN( play_.GetOp( ind_taddr_[i] ) == local::InvOp );
336 //
337 // Use add_element when only adding one element per set is added.
338 for_jac_pattern.add_element( ind_taddr_[i], ind_taddr_[i] );
339 }
340 // compute forward Jacobiain sparsity pattern
341 bool dependency = false;
342 local::for_jac_sweep(
343 &play_,
344 dependency,
345 n,
346 num_var_tape_,
347 for_jac_pattern
348 );
349 // sparsity pattern correspnding to s
350 local::sparse_list rev_jac_pattern;
351 rev_jac_pattern.resize(num_var_tape_, 1);
352 itr_1 = s[0].begin();
353 while( itr_1 != s[0].end() )
354 { size_t i = *itr_1++;
355 CPPAD_ASSERT_KNOWN(
356 i < m,
357 "ForSparseHes: an element of the set s[0] has value "
358 "greater than or equal m"
359 );
360 CPPAD_ASSERT_UNKNOWN( dep_taddr_[i] < num_var_tape_ );
361 //
362 // Use add_element when only adding one element per set is added.
363 rev_jac_pattern.add_element( dep_taddr_[i], 0);
364 }
365 //
366 // compute reverse sparsity pattern for dependency analysis
367 // (note that we are only want non-zero derivatives not true dependency)
368 local::rev_jac_sweep(
369 &play_,
370 dependency,
371 n,
372 num_var_tape_,
373 rev_jac_pattern
374 );
375 //
376 // vector of sets that will hold reverse Hessain values
377 local::sparse_list for_hes_pattern;
378 for_hes_pattern.resize(n+1, n+1);
379 //
380 // compute the Hessian sparsity patterns
381 local::for_hes_sweep(
382 &play_,
383 n,
384 num_var_tape_,
385 for_jac_pattern,
386 rev_jac_pattern,
387 for_hes_pattern
388 );
389 // return values corresponding to independent variables
390 // j is index corresponding to reverse mode partial
391 h.resize(n);
392 CPPAD_ASSERT_UNKNOWN( for_hes_pattern.end() == n+1 );
393 for(size_t i = 0; i < n; i++)
394 { CPPAD_ASSERT_UNKNOWN( ind_taddr_[i] == i + 1 );
395 CPPAD_ASSERT_UNKNOWN( play_.GetOp( ind_taddr_[i] ) == local::InvOp );
396
397 // extract the result from for_hes_pattern
398 local::sparse_list::const_iterator itr_2(for_hes_pattern, ind_taddr_[i] );
399 size_t j = *itr_2;
400 while( j < for_hes_pattern.end() )
401 { CPPAD_ASSERT_UNKNOWN( 0 < j )
402 h[i].insert(j-1);
403 j = *(++itr_2);
404 }
405 }
406 }
407
408 // ===========================================================================
409 // ForSparseHes
410
411 /*!
412 User API for Hessian sparsity patterns using reverse mode.
413
414 The C++ source code corresponding to this operation is
415 \verbatim
416 h = f.ForSparseHes(q, r)
417 \endverbatim
418
419 \tparam Base
420 is the base type for this recording.
421
422 \tparam VectorSet
423 is a simple vector with elements of type \c bool
424 or \c std::set<size_t>.
425
426 \param r
427 is a vector with size \c n that specifies the sparsity pattern
428 for the diagonal of the matrix \f$ R \f$,
429 where \c n is the number of independent variables
430 corresponding to the operation sequence stored in \a play.
431
432 \param s
433 is a vector with size \c m that specifies the sparsity pattern
434 for the vector \f$ S \f$,
435 where \c m is the number of dependent variables
436 corresponding to the operation sequence stored in \a play.
437
438 \return
439 The return vector is a sparsity pattern for \f$ H(x) \f$
440 \f[
441 H(x) = R^T ( S * F)^{(2)} (x) R
442 \f]
443 where \f$ F \f$ is the function corresponding to the operation sequence
444 and \a x is any argument value.
445 */
446
447 template <class Base>
448 template <class VectorSet>
ForSparseHes(const VectorSet & r,const VectorSet & s)449 VectorSet ADFun<Base>::ForSparseHes(
450 const VectorSet& r, const VectorSet& s
451 )
452 { VectorSet h;
453 typedef typename VectorSet::value_type Set_type;
454
455 // Should check to make sure q is same as in previous call to
456 // forward sparse Jacobian.
457 ForSparseHesCase(
458 Set_type() ,
459 r ,
460 s ,
461 h
462 );
463
464 return h;
465 }
466 // ===========================================================================
467 // ForSparseHesCheckpoint
468 /*!
469 Hessian sparsity patterns calculation used by checkpoint functions.
470
471 \tparam Base
472 is the base type for this recording.
473
474 \param r
475 is a vector with size n that specifies the sparsity pattern
476 for the diagonal of \f$ R \f$,
477 where n is the number of independent variables
478 corresponding to the operation sequence stored in play_.
479
480 \param s
481 is a vector with size m that specifies the sparsity pattern
482 for the vector \f$ S \f$,
483 where m is the number of dependent variables
484 corresponding to the operation sequence stored in play_.
485
486 \param h
487 The input size and elements of h do not matter.
488 On output, h is the sparsity pattern for the matrix \f$ H(x) R \f$.
489
490 \par Assumptions
491 The forward jacobian sparsity pattern must be currently stored
492 in this ADFUN object.
493 */
494
495 // The checkpoint class is not yet using forward sparse Hessians.
496 # ifdef CPPAD_NOT_DEFINED
497 template <class Base>
ForSparseHesCheckpoint(vector<bool> & r,vector<bool> & s,local::sparse_list & h)498 void ADFun<Base>::ForSparseHesCheckpoint(
499 vector<bool>& r ,
500 vector<bool>& s ,
501 local::sparse_list& h )
502 {
503 size_t n = Domain();
504 size_t m = Range();
505
506 // checkpoint functions should get this right
507 CPPAD_ASSERT_UNKNOWN( for_jac_sparse_pack_.n_set() == 0 );
508 CPPAD_ASSERT_UNKNOWN( for_jac_sparse_set_.n_set() == 0 );
509 CPPAD_ASSERT_UNKNOWN( s.size() == m );
510
511 // Array that holds the reverse Jacobiain dependcy flags.
512 // Initialize as true for dependent variables, flase for others.
513 local::pod_vector<bool> RevJac(num_var_tape_);
514 for(size_t i = 0; i < num_var_tape_; i++)
515 RevJac[i] = false;
516 for(size_t i = 0; i < m; i++)
517 { CPPAD_ASSERT_UNKNOWN( dep_taddr_[i] < num_var_tape_ )
518 RevJac[ dep_taddr_[i] ] = s[i];
519 }
520
521 // holds forward Hessian sparsity pattern for all variables
522 local::sparse_list for_hes_pattern;
523 for_hes_pattern.resize(n+1, n+1);
524
525 // compute Hessian sparsity pattern for all variables
526 local::for_hes_sweep(
527 &play_,
528 n,
529 num_var_tape_,
530 for_jac_sparse_set_,
531 RevJac.data(),
532 for_hes_pattern
533 );
534
535 // dimension the return value
536 if( transpose )
537 h.resize(n, n);
538 else
539 h.resize(n, n);
540
541 // j is index corresponding to reverse mode partial
542 for(size_t j = 0; j < n; j++)
543 { CPPAD_ASSERT_UNKNOWN( ind_taddr_[j] < num_var_tape_ );
544
545 // ind_taddr_[j] is operator taddr for j-th independent variable
546 CPPAD_ASSERT_UNKNOWN( ind_taddr_[j] == j + 1 );
547 CPPAD_ASSERT_UNKNOWN( play_.GetOp( ind_taddr_[j] ) == local::InvOp );
548
549 // extract the result from for_hes_pattern
550 CPPAD_ASSERT_UNKNOWN( for_hes_pattern.end() == q );
551 local::sparse_list::const_iterator itr(for_hes_pattern, .j + 1);
552 size_t i = *itr;
553 while( i < q )
554 { if( transpose )
555 h.post_element(j, i);
556 else h.post_element(i, j);
557 i = *(++itr);
558 }
559 }
560 // process posts
561 for(size_t i = 0; i < n; ++i)
562 h.process_post(i);
563 }
564 # endif
565
566 } // END_CPPAD_NAMESPACE
567 # endif
568