CppAD: A C++ Algorithmic Differentiation Package  20171217
rev_sparse_jac.hpp
Go to the documentation of this file.
3
4 /* --------------------------------------------------------------------------
6
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.
13 -------------------------------------------------------------------------- */
14
15 /*
16 $begin RevSparseJac$$17 spell 18 optimizer 19 nz 20 CondExpRel 21 std 22 VecAD 23 var 24 Jacobian 25 Jac 26 const 27 Bool 28 Dep 29 proportional 30$$ 31 32$section Jacobian Sparsity Pattern: Reverse Mode$$33 mindex RevSparseJac sparse$$
34
35 $head Syntax$$36 icode%s% = %f%.RevSparseJac(%q%, %r%) 37 %$$ 38$icode%s% = %f%.RevSparseJac(%q%, %r%, %transpose%, %dependency%)%$$39 40 head Purpose$$
41 We use $latex F : B^n \rightarrow B^m$$to denote the 42 cref/AD function/glossary/AD Function/$$ corresponding to$icode f$$. 43 For a fixed matrix latex R \in B^{q \times m}$$,
44 the Jacobian of $latex R * F( x )$$45 with respect to latex x$$ is 46$latex $47 S(x) = R * F^{(1)} ( x ) 48$ $$49 Given a 50 cref/sparsity pattern/glossary/Sparsity Pattern/$$
51 for $latex R$$, 52 code RevSparseJac$$ returns a sparsity pattern for the$latex S(x)$$. 53 54 head f$$
55 The object $icode f$$has prototype 56 codei% 57 ADFun<%Base%> %f% 58 %$$ 59 60$head x$$61 If the operation sequence in icode f$$ is
62 $cref/independent/glossary/Operation/Independent/$$of 63 the independent variables in latex x \in B^n$$, 64 the sparsity pattern is valid for all values of 65 (even if it has$cref CondExp$$or cref VecAD$$ operations).
66
67 $head q$$68 The argument icode q$$ has prototype 69$codei%
70  size_t %q%
71 %$$72 It specifies the number of rows in 73 latex R \in B^{q \times m}$$ and the
74 Jacobian $latex S(x) \in B^{q \times n}$$. 75 76 head transpose$$ 77 The argument$icode transpose$$has prototype 78 codei% 79 bool %transpose% 80 %$$
81 The default value $code false$$is used when icode transpose$$ is not present. 82 83$head dependency$$84 The argument icode dependency$$ has prototype
85 $codei% 86 bool %dependency% 87 %$$88 If icode dependency$$ is true, 89 the$cref/dependency pattern/dependency.cpp/Dependency Pattern/$$90 (instead of sparsity pattern) is computed. 91 92 head r$$
93 The argument $icode s$$has prototype 94 codei% 95 const %VectorSet%& %r% 96 %$$ 97 see$cref/VectorSet/RevSparseJac/VectorSet/$$below. 98 99 subhead transpose false$$
100 If $icode r$$has elements of type code bool$$, 101 its size is$latex q * m$$. 102 If it has elements of type code std::set<size_t>$$,
103 its size is $icode q$$and all its set elements are between 104 zero and latex m - 1$$. 105 It specifies a 106$cref/sparsity pattern/glossary/Sparsity Pattern/$$107 for the matrix latex R \in B^{q \times m}$$.
108
109 $subhead transpose true$$110 If icode r$$ has elements of type$code bool$$, 111 its size is latex m * q$$.
112 If it has elements of type $code std::set<size_t>$$, 113 its size is icode m$$ and all its set elements are between 114 zero and$latex q - 1$$. 115 It specifies a 116 cref/sparsity pattern/glossary/Sparsity Pattern/$$
117 for the matrix $latex R^\R{T} \in B^{m \times q}$$. 118 119 head s$$ 120 The return value$icode s$$has prototype 121 codei% 122 %VectorSet% %s% 123 %$$
124 see $cref/VectorSet/RevSparseJac/VectorSet/$$below. 125 126 subhead transpose false$$ 127 If it has elements of type$code bool$$, 128 its size is latex q * n$$.
129 If it has elements of type $code std::set<size_t>$$, 130 its size is icode q$$ and all its set elements are between 131 zero and$latex n - 1$$. 132 It specifies a 133 cref/sparsity pattern/glossary/Sparsity Pattern/$$
134 for the matrix $latex S(x) \in {q \times n}$$. 135 136 subhead transpose true$$ 137 If it has elements of type$code bool$$, 138 its size is latex n * q$$.
139 If it has elements of type $code std::set<size_t>$$, 140 its size is icode n$$ and all its set elements are between 141 zero and$latex q - 1$$. 142 It specifies a 143 cref/sparsity pattern/glossary/Sparsity Pattern/$$
144 for the matrix $latex S(x)^\R{T} \in {n \times q}$$. 145 146 head VectorSet$$ 147 The type$icode VectorSet$$must be a cref SimpleVector$$ class with
148 $cref/elements of type/SimpleVector/Elements of Specified Type/$$149 code bool$$ or$code std::set<size_t>$$; 150 see cref/sparsity pattern/glossary/Sparsity Pattern/$$ for a discussion
151 of the difference.
152
153 $head Entire Sparsity Pattern$$154 Suppose that latex q = m$$ and 155$latex R$$is the latex m \times m$$ identity matrix.
156 In this case,
157 the corresponding value for $icode s$$is a 158 sparsity pattern for the Jacobian latex S(x) = F^{(1)} ( x )$$. 159 160$head Example$$161 children% 162 example/sparse/rev_sparse_jac.cpp 163 %$$
164 The file
165 $cref rev_sparse_jac.cpp$$166 contains an example and test of this operation. 167 It returns true if it succeeds and false otherwise. 168 169$end
170 -----------------------------------------------------------------------------
171 */
172
174
176 /*!
177 \file rev_sparse_jac.hpp
178 Reverse mode Jacobian sparsity patterns.
179 */
180 // =========================================================================
181 // RevSparseJacCase
182
183 /*!
184 Private helper function for RevSparseJac(q, r, transpose) boolean sparsity.
185
186 All of the description in the public member function
187 \c RevSparseJac(q, r, transpose) apply.
188
189 \param set_type
190 is a \c bool value.
191 This argument is used to dispatch to the proper source code
192 depending on the value of \c VectorSet::value_type.
193
194 \param transpose
195 See \c RevSparseJac(q, r, transpose, dependency)
196
197 \param dependency
198 See \c RevSparseJac(q, r, transpose, dependency)
199
200 \param q
201 See \c RevSparseJac(q, r, transpose, dependency)
202
203 \param r
204 See \c RevSparseJac(q, r, transpose, dependency)
205
206 \param s
207 is the return value for the corresponding call to
208 RevSparseJac(q, r, transpose).
209 */
210
211 template <class Base>
212 template <class VectorSet>
214  bool set_type ,
215  bool transpose ,
216  bool dependency ,
217  size_t q ,
218  const VectorSet& r ,
219  VectorSet& s )
220 { size_t n = Domain();
221  size_t m = Range();
222
223  // dimension of the result vector
224  s.resize( q * n );
225
226  // check VectorSet is Simple Vector class with bool elements
227  CheckSimpleVector<bool, VectorSet>();
228  //
230  q > 0,
231  "RevSparseJac: q is not greater than zero"
232  );
234  size_t(r.size()) == q * m,
235  "RevSparseJac: size of r is not equal to\n"
236  "q times range dimension for ADFun object."
237  );
238  //
239  // vector of sets that will hold the results
240  local::sparse_pack var_sparsity;
241  var_sparsity.resize(num_var_tape_, q);
242
243  // The sparsity pattern corresponding to the dependent variables
244  for(size_t i = 0; i < m; i++)
246  if( transpose )
247  { for(size_t j = 0; j < q; j++) if( r[ i * q + j ] )
249  }
250  else
251  { for(size_t j = 0; j < q; j++) if( r[ j * m + i ] )
253  }
254  }
255  // process posts
256  for(size_t i = 0; i < m; i++)
258
259  // evaluate the sparsity patterns
261  &play_,
262  dependency,
263  n,
264  num_var_tape_,
265  var_sparsity
266  );
267
268  // return values corresponding to dependent variables
269  CPPAD_ASSERT_UNKNOWN( size_t(s.size()) == q * n );
270  for(size_t j = 0; j < n; j++)
272
275
276  // extract the result from var_sparsity
277  if( transpose )
278  { for(size_t i = 0; i < q; i++)
279  s[ j * q + i ] = false;
280  }
281  else
282  { for(size_t i = 0; i < q; i++)
283  s[ i * n + j ] = false;
284  }
285  CPPAD_ASSERT_UNKNOWN( var_sparsity.end() == q );
286  local::sparse_pack::const_iterator itr(var_sparsity, j+1);
287  size_t i = *itr;
288  while( i < q )
289  { if( transpose )
290  s[ j * q + i ] = true;
291  else s[ i * n + j ] = true;
292  i = *(++itr);
293  }
294  }
295 }
296
297 /*!
298 Private helper function for RevSparseJac(q, r, transpose) set sparsity
299
300 All of the description in the public member function
301 \c RevSparseJac(q, r, transpose) apply.
302
303 \param set_type
304 is a \c std::set<size_t> object.
305 This argument is used to dispatch to the proper source code
306 depending on the value of \c VectorSet::value_type.
307
308 \param transpose
309 See \c RevSparseJac(q, r, transpose, dependency)
310
311 \param dependency
312 See \c RevSparseJac(q, r, transpose, dependency)
313
314 \param q
315 See \c RevSparseJac(q, r, transpose, dependency)
316
317 \param r
318 See \c RevSparseJac(q, r, transpose, dependency)
319
320 \param s
321 is the return value for the corresponding call to RevSparseJac(q, r, transpose)
322 */
323
324 template <class Base>
325 template <class VectorSet>
327  const std::set<size_t>& set_type ,
328  bool transpose ,
329  bool dependency ,
330  size_t q ,
331  const VectorSet& r ,
332  VectorSet& s )
333 { // dimension of the result vector
334  if( transpose )
335  s.resize( Domain() );
336  else s.resize( q );
337
338  // temporary indices
339  std::set<size_t>::const_iterator itr_1;
340
341  // check VectorSet is Simple Vector class with sets for elements
342  CheckSimpleVector<std::set<size_t>, VectorSet>(
343  local::one_element_std_set<size_t>(), local::two_element_std_set<size_t>()
344  );
345
346  // domain dimensions for F
349
351  q > 0,
352  "RevSparseJac: q is not greater than zero"
353  );
355  size_t(r.size()) == q || transpose,
356  "RevSparseJac: size of r is not equal to q and transpose is false."
357  );
359  size_t(r.size()) == m || ! transpose,
360  "RevSparseJac: size of r is not equal to m and transpose is true."
361  );
362
363  // vector of lists that will hold the results
364  local::sparse_list var_sparsity;
365  var_sparsity.resize(num_var_tape_, q);
366
367  // The sparsity pattern corresponding to the dependent variables
368  if( transpose )
369  { for(size_t i = 0; i < m; i++)
370  { itr_1 = r[i].begin();
371  while(itr_1 != r[i].end())
372  { size_t j = *itr_1++;
374  j < q,
375  "RevSparseJac: transpose is true and element of the set\n"
376  "r[i] has value greater than or equal q."
377  );
380  }
381  }
382  }
383  else
384  { for(size_t i = 0; i < q; i++)
385  { itr_1 = r[i].begin();
386  while(itr_1 != r[i].end())
387  { size_t j = *itr_1++;
389  j < m,
390  "RevSparseJac: transpose is false and element of the set\n"
391  "r[i] has value greater than or equal range dimension."
392  );
395  }
396  }
397  }
398  // process posts
399  for(size_t i = 0; i < m; i++)
401
402  // evaluate the sparsity patterns
404  &play_,
405  dependency,
406  n,
407  num_var_tape_,
408  var_sparsity
409  );
410
411  // return values corresponding to dependent variables
412  CPPAD_ASSERT_UNKNOWN( size_t(s.size()) == q || transpose );
413  CPPAD_ASSERT_UNKNOWN( size_t(s.size()) == n || ! transpose );
414  for(size_t j = 0; j < n; j++)
416
419
420  CPPAD_ASSERT_UNKNOWN( var_sparsity.end() == q );
421  local::sparse_list::const_iterator itr_2(var_sparsity, j+1);
422  size_t i = *itr_2;
423  while( i < q )
424  { if( transpose )
425  s[j].insert(i);
426  else s[i].insert(j);
427  i = *(++itr_2);
428  }
429  }
430 }
431
432 // =========================================================================
433 // RevSparseJac
434 /*!
435 User API for Jacobian sparsity patterns using reverse mode.
436
437 The C++ source code corresponding to this operation is
438 \verbatim
439  s = f.RevSparseJac(q, r, transpose, dependency)
440 \endverbatim
441
442 \tparam Base
443 is the base type for this recording.
444
445 \tparam VectorSet
446 is a simple vector with elements of type \c bool.
447 or \c std::set<size_t>.
448
449 \param q
450 is the number of rows in the matrix \f$R \f$.
451
452 \param r
453 is a sparsity pattern for the matrix \f$R \f$.
454
455 \param transpose
456 are the sparsity patterns for \f$R \f$ and \f$S(x) \f$ transposed.
457
458 \param dependency
459 Are the derivatives with respect to left and right of the expression below
460 considered to be non-zero:
461 \code
462  CondExpRel(left, right, if_true, if_false)
463 \endcode
464 This is used by the optimizer to obtain the correct dependency relations.
465
466
467 \return
468 If \c transpose is false (true), the return value is a sparsity pattern
469 for \f$S(x) \f$ (\f$S(x)^T \f$) where
470 \f[
471  S(x) = R * F^{(1)} (x)
472 \f]
473 and \f$F \f$ is the function corresponding to the operation sequence
474 and \a x is any argument value.
475 If \c VectorSet::value_type is \c bool,
476 the return value has size \f$q * n \f$ ( \f$n * q \f$).
477 If \c VectorSet::value_type is \c std::set<size_t>,
478 the return value has size \f$q \f$ ( \f$n \f$)
479 and with all its elements between zero and \f$n - 1 \f$ (\f$q - 1 \f$).
480 */
481 template <class Base>
482 template <class VectorSet>
484  size_t q ,
485  const VectorSet& r ,
486  bool transpose ,
487  bool dependency )
488 {
489  VectorSet s;
490  typedef typename VectorSet::value_type Set_type;
491
492  RevSparseJacCase(
493  Set_type() ,
494  transpose ,
495  dependency ,
496  q ,
497  r ,
498  s
499  );
500  return s;
501 }
502 // ===========================================================================
503 // RevSparseJacCheckpoint
504 /*!
505 Reverse mode Jacobian sparsity calculation used by checkpoint functions.
506
507 \tparam Base
508 is the base type for this recording.
509
510 \param transpose
511 is true (false) s is equal to \f$S(x) \f$ (\f$S(x)^T \f$)
512 where
513 \f[
514  S(x) = R * F^{(1)} (x)
515 \f]
516 where \f$F \f$ is the function corresponding to the operation sequence
517 and \f$x \f$ is any argument value.
518
519 \param q
520 is the number of rows in the matrix \f$R \f$.
521
522 \param r
523 is a sparsity pattern for the matrix \f$R \f$.
524
525 \param transpose
526 are the sparsity patterns for \f$R \f$ and \f$S(x) \f$ transposed.
527
528 \param dependency
529 Are the derivatives with respect to left and right of the expression below
530 considered to be non-zero:
531 \code
532  CondExpRel(left, right, if_true, if_false)
533 \endcode
534 This is used by the optimizer to obtain the correct dependency relations.
535
536 \param s
537 The input size and elements of s do not matter.
538 On output, s is the sparsity pattern for the matrix \f$S(x) \f$
539 or \f$S(x)^T \f$ depending on transpose.
540
541 */
542 template <class Base>
544  size_t q ,
545  const local::sparse_list& r ,
546  bool transpose ,
547  bool dependency ,
548  local::sparse_list& s )
549 { size_t n = Domain();
550  size_t m = Range();
551
552 # ifndef NDEBUG
553  if( transpose )
554  { CPPAD_ASSERT_UNKNOWN( r.n_set() == m );
555  CPPAD_ASSERT_UNKNOWN( r.end() == q );
556  }
557  else
558  { CPPAD_ASSERT_UNKNOWN( r.n_set() == q );
559  CPPAD_ASSERT_UNKNOWN( r.end() == m );
560  }
561  for(size_t i = 0; i < m; i++)
563 # endif
564
565  // holds reverse Jacobian sparsity pattern for all variables
566  local::sparse_list var_sparsity;
567  var_sparsity.resize(num_var_tape_, q);
568
569  // set sparsity pattern for dependent variables
570  if( transpose )
571  { for(size_t i = 0; i < m; i++)
573  size_t j = *itr;
574  while( j < q )
575  { var_sparsity.post_element( dep_taddr_[i], j );
576  j = *(++itr);
577  }
578  }
579  }
580  else
581  { for(size_t j = 0; j < q; j++)
583  size_t i = *itr;
584  while( i < m )
585  { var_sparsity.post_element( dep_taddr_[i], j );
586  i = *(++itr);
587  }
588  }
589  }
590  // process posts
591  for(size_t i = 0; i < m; i++)
593
594  // evaluate the sparsity pattern for all variables
596  &play_,
597  dependency,
598  n,
599  num_var_tape_,
600  var_sparsity
601  );
602
603  // dimension the return value
604  if( transpose )
605  s.resize(n, m);
606  else
607  s.resize(m, n);
608
609  // return values corresponding to independent variables
610  for(size_t j = 0; j < n; j++)
612
615
616  // extract the result from var_sparsity
617  CPPAD_ASSERT_UNKNOWN( var_sparsity.end() == q );
618  local::sparse_list::const_iterator itr(var_sparsity, j+1);
619  size_t i = *itr;
620  while( i < q )
621  { if( transpose )
622  s.post_element(j, i);
623  else
624  s.post_element(i, j);
625  i = *(++itr);
626  }
627  }
628  // process posts
629  for(size_t i = 0; i < s.n_set(); i++)
630  s.process_post(i);
631
632 }
633
635 # endif
Check that exp is true, if not print msg and terminate execution.
size_t end(void) const
Fetch end for this vector of sets object.
Vector of sets of positive integers, each set stored as a singly linked list.
Definition: sparse_list.hpp:35
void process_post(size_t i)
process post entries for a specific set.
size_t end(void) const
Fetch end for this vector of sets object.
void post_element(size_t i, size_t element)
Post an element for delayed addition to a set.
void RevSparseJacCase(bool set_type, bool transpose, bool dependency, size_t p, const VectorSet &s, VectorSet &r)
Private helper function for RevSparseJac(q, r, transpose) boolean sparsity.
VectorSet RevSparseJac(size_t q, const VectorSet &s, bool transpose=false, bool dependency=false)
User API for Jacobian sparsity patterns using reverse mode.
Vector of sets of postivie integers, each set stored as a packed boolean array.
Definition: sparse_pack.hpp:33
void post_element(size_t i, size_t element)
Post an element for delayed addition to a set.
void RevSparseJacCheckpoint(size_t q, const local::sparse_list &r, bool transpose, bool dependency, local::sparse_list &s)
Reverse mode Jacobian sparsity calculation used by checkpoint functions.
Two constant standard sets (currently used for concept checking).
Check that exp is true, if not terminate execution.
void resize(size_t n_set, size_t end)
Start a new vector of sets.
cons_iterator for one set of positive integers in a sparse_pack object.
void process_post(size_t i)
process post entries for a specific set.
size_t n_set(void) const
Fetch n_set for vector of sets object.
cons_iterator for one set of positive integers in a sparse_list object.
void rev_jac_sweep(const local::player< Base > *play, bool dependency, size_t n, size_t numvar, Vector_set &var_sparsity)
Given the sparsity pattern for the dependent variables, RevJacSweep computes the sparsity pattern for...
Scalar value_type
void resize(size_t n_set, size_t end)
Change number of sets, set end, and initialize all sets as empty.