CppAD: A C++ Algorithmic Differentiation Package  20171217
for_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 ForSparseJac$$17 spell 18 std 19 var 20 Jacobian 21 Jac 22 const 23 Bool 24 proportional 25 VecAD 26 CondExpRel 27 optimizer 28 cpp 29$$ 30 31$section Jacobian Sparsity Pattern: Forward Mode$$32 33 head Syntax$$
34 $icode%s% = %f%.ForSparseJac(%q%, %r%) 35 %$$36 icode%s% = %f%.ForSparseJac(%q%, %r%, %transpose%, %dependency%)%$$ 37 38$head Purpose$$39 We use latex F : B^n \rightarrow B^m$$ to denote the
40 $cref/AD function/glossary/AD Function/$$corresponding to icode f$$. 41 For a fixed$latex n \times q$$matrix latex R$$,
42 the Jacobian of $latex F[ x + R * u ]$$43 with respect to latex u$$ at$latex u = 0$$is 44 latex $45 S(x) = F^{(1)} ( x ) * R 46$$$
47 Given a
48 $cref/sparsity pattern/glossary/Sparsity Pattern/$$49 for latex R$$, 50$code ForSparseJac$$returns a sparsity pattern for the latex S(x)$$.
51
52 $head f$$53 The object icode f$$ has prototype 54$codei%
56 %$$57 Note that the cref ADFun$$ object $icode f$$is not code const$$. 58 After a call to$code ForSparseJac$$, the sparsity pattern 59 for each of the variables in the operation sequence 60 is held in icode f$$ (for possible later use by $cref RevSparseHes$$). 61 These sparsity patterns are stored with elements of type code bool$$ 62 or elements of type$code std::set<size_t>$$63 (see cref/VectorSet/ForSparseJac/VectorSet/$$ below).
64
65 $subhead size_forward_bool$$66 After code ForSparseJac$$, if$icode k$$is a code size_t$$ object,
67 $codei% 68 %k% = %f%.size_forward_bool() 69 %$$70 sets icode k$$ to the amount of memory (in unsigned character units) 71 used to store the sparsity pattern with elements of type$code bool$$72 in the function object icode f$$.
73 If the sparsity patterns for the previous $code ForSparseJac$$used 74 elements of type code bool$$, 75 the return value for$code size_forward_bool$$will be non-zero. 76 Otherwise, its return value will be zero. 77 This sparsity pattern is stored for use by cref RevSparseHes$$ and
78 when it is not longer needed, it can be deleted
79 (and the corresponding memory freed) using
80 $codei% 81 %f%.size_forward_bool(0) 82 %$$83 After this call, icode%f%.size_forward_bool()%$$ will return zero. 84 85$subhead size_forward_set$$86 After code ForSparseJac$$, if $icode k$$is a code size_t$$ object, 87$codei%
88  %k% = %f%.size_forward_set()
89 %$$90 sets icode k$$ to the amount of memory (in unsigned character units)
91 used to store the
92 $cref/vector of sets/glossary/Sparsity Pattern/Vector of Sets/$$93 sparsity patterns. 94 If the sparsity patterns for this operation use elements of type code bool$$, 95 the return value for$code size_forward_set$$will be zero. 96 Otherwise, its return value will be non-zero. 97 This sparsity pattern is stored for use by cref RevSparseHes$$ and
98 when it is not longer needed, it can be deleted
99 (and the corresponding memory freed) using
100 $codei% 101 %f%.size_forward_set(0) 102 %$$103 After this call, icode%f%.size_forward_set()%$$ will return zero. 104 105$head x$$106 If the operation sequence in icode f$$ is
107 $cref/independent/glossary/Operation/Independent/$$of 108 the independent variables in latex x \in B^n$$, 109 the sparsity pattern is valid for all values of 110 (even if it has$cref CondExp$$or cref VecAD$$ operations).
111
112 $head q$$113 The argument icode q$$ has prototype 114$codei%
115  size_t %q%
116 %$$117 It specifies the number of columns in 118 latex R \in B^{n \times q}$$ and the Jacobian
119 $latex S(x) \in B^{m \times q}$$. 120 121 head transpose$$ 122 The argument$icode transpose$$has prototype 123 codei% 124 bool %transpose% 125 %$$
126 The default value $code false$$is used when icode transpose$$ is not present. 127 128$head dependency$$129 The argument icode dependency$$ has prototype
130 $codei% 131 bool %dependency% 132 %$$133 If icode dependency$$ is true, 134 the$cref/dependency pattern/dependency.cpp/Dependency Pattern/$$135 (instead of sparsity pattern) is computed. 136 137 head r$$
138 The argument $icode r$$has prototype 139 codei% 140 const %VectorSet%& %r% 141 %$$ 142 see$cref/VectorSet/ForSparseJac/VectorSet/$$below. 143 144 subhead transpose false$$
145 If $icode r$$has elements of type code bool$$, 146 its size is$latex n * q$$. 147 If it has elements of type code std::set<size_t>$$,
148 its size is $latex n$$and all the set elements must be between 149 zero and icode%q%-1%$$ inclusive. 150 It specifies a 151$cref/sparsity pattern/glossary/Sparsity Pattern/$$152 for the matrix latex R \in B^{n \times q}$$.
153
154 $subhead transpose true$$155 If icode r$$ has elements of type$code bool$$, 156 its size is latex q * n$$.
157 If it has elements of type $code std::set<size_t>$$, 158 its size is latex q$$ and all the set elements must be between 159 zero and$icode%n%-1%$$inclusive. 160 It specifies a 161 cref/sparsity pattern/glossary/Sparsity Pattern/$$
162 for the matrix $latex R^\R{T} \in B^{q \times n}$$. 163 164 head s$$ 165 The return value$icode s$$has prototype 166 codei% 167 %VectorSet% %s% 168 %$$
169 see $cref/VectorSet/ForSparseJac/VectorSet/$$below. 170 171 subhead transpose false$$ 172 If$icode s$$has elements of type code bool$$,
173 its size is $latex m * q$$. 174 If it has elements of type code std::set<size_t>$$, 175 its size is$latex m$$and all its set elements are between 176 zero and icode%q%-1%$$ inclusive.
177 It specifies a
178 $cref/sparsity pattern/glossary/Sparsity Pattern/$$179 for the matrix latex S(x) \in B^{m \times q}$$. 180 181$subhead transpose true$$182 If icode s$$ has elements of type $code bool$$, 183 its size is latex q * m$$. 184 If it has elements of type$code std::set<size_t>$$, 185 its size is latex q$$ and all its set elements are between
186 zero and $icode%m%-1%$$inclusive. 187 It specifies a 188 cref/sparsity pattern/glossary/Sparsity Pattern/$$ 189 for the matrix$latex S(x)^\R{T} \in B^{q \times m}$$. 190 191 head VectorSet$$
192 The type $icode VectorSet$$must be a cref SimpleVector$$ class with 193$cref/elements of type/SimpleVector/Elements of Specified Type/$$194 code bool$$ or $code std::set<size_t>$$; 195 see cref/sparsity pattern/glossary/Sparsity Pattern/$$ for a discussion 196 of the difference. 197 198$head Entire Sparsity Pattern$$199 Suppose that latex q = n$$ and
200 $latex R$$is the latex n \times n$$ identity matrix. 201 In this case, 202 the corresponding value for$icode s$$is a 203 sparsity pattern for the Jacobian latex S(x) = F^{(1)} ( x )$$.
204
205 $head Example$$206 children% 207 example/sparse/for_sparse_jac.cpp 208 %$$ 209 The file 210$cref for_sparse_jac.cpp$$211 contains an example and test of this operation. 212 It returns true if it succeeds and false otherwise. 213 The file 214 cref/sparsity_sub.cpp/sparsity_sub.cpp/ForSparseJac/$$
215 contains an example and test of using $code ForSparseJac$$216 to compute the sparsity pattern for a subset of the Jacobian. 217 218$end
219 -----------------------------------------------------------------------------
220 */
221
223
225 /*!
226 \file for_sparse_jac.hpp
227 Forward mode Jacobian sparsity patterns.
228 */
229 // ---------------------------------------------------------------------------
230 /*!
231 Private helper function for ForSparseJac(q, r) boolean sparsity patterns.
232
233 All of the description in the public member function ForSparseJac(q, r)
234 applies.
235
236 \param set_type
237 is a \c bool value. This argument is used to dispatch to the proper source
238 code depending on the value of \c VectorSet::value_type.
239
240 \param transpose
241 See \c ForSparseJac(q, r, transpose, dependency).
242
243 \param dependency
244 See \c ForSparseJac(q, r, transpose, dependency).
245
246 \param q
247 See \c ForSparseJac(q, r, transpose, dependency).
248
249 \param r
250 See \c ForSparseJac(q, r, transpose, dependency).
251
252 \param s
253 is the return value for the corresponding call to \c ForSparseJac(q, r).
254 */
255
256 template <class Base>
257 template <class VectorSet>
259  bool set_type ,
260  bool transpose ,
261  bool dependency ,
262  size_t q ,
263  const VectorSet& r ,
264  VectorSet& s )
265 { size_t m = Range();
266  size_t n = Domain();
267
268  // check VectorSet is Simple Vector class with bool elements
269  CheckSimpleVector<bool, VectorSet>();
270
271  // dimension size of result vector
272  s.resize( m * q );
273
275  q > 0,
276  "ForSparseJac: q is not greater than zero"
277  );
279  size_t(r.size()) == n * q,
280  "ForSparseJac: size of r is not equal to\n"
281  "q times domain dimension for ADFun object."
282  );
283  //
284  // allocate memory for the requested sparsity calculation result
285  for_jac_sparse_pack_.resize(num_var_tape_, q);
286
287  // set values corresponding to independent variables
288  for(size_t i = 0; i < n; i++)
292
293  // set bits that are true
294  if( transpose )
295  { for(size_t j = 0; j < q; j++) if( r[ j * n + i ] )
297  }
298  else
299  { for(size_t j = 0; j < q; j++) if( r[ i * q + j ] )
301  }
302  }
303  // process posts
304  for(size_t j = 0; j < n; j++)
306
307  // evaluate the sparsity patterns
309  &play_,
310  dependency,
311  n,
312  num_var_tape_,
313  for_jac_sparse_pack_
314  );
315
316  // return values corresponding to dependent variables
317  CPPAD_ASSERT_UNKNOWN( size_t(s.size()) == m * q );
318  for(size_t i = 0; i < m; i++)
320
321  // extract the result from for_jac_sparse_pack_
322  if( transpose )
323  { for(size_t j = 0; j < q; j++)
324  s[ j * m + i ] = false;
325  }
326  else
327  { for(size_t j = 0; j < q; j++)
328  s[ i * q + j ] = false;
329  }
330  CPPAD_ASSERT_UNKNOWN( for_jac_sparse_pack_.end() == q );
333  size_t j = *itr;
334  while( j < q )
335  { if( transpose )
336  s[j * m + i] = true;
337  else s[i * q + j] = true;
338  j = *(++itr);
339  }
340  }
341 }
342 // ---------------------------------------------------------------------------
343 /*!
344 Private helper function for \c ForSparseJac(q, r) set sparsity.
345
346 All of the description in the public member function \c ForSparseJac(q, r)
347 applies.
348
349 \param set_type
350 is a \c std::set<size_t> object.
351 This argument is used to dispatch to the proper source
352 code depending on the value of \c VectorSet::value_type.
353
354 \param transpose
355 See \c ForSparseJac(q, r, transpose, dependency).
356
357 \param dependency
358 See \c ForSparseJac(q, r, transpose, dependency).
359
360 \param q
361 See \c ForSparseJac(q, r, transpose, dependency).
362
363 \param r
364 See \c ForSparseJac(q, r, transpose, dependency).
365
366 \param s
367 is the return value for the corresponding call to \c ForSparseJac(q, r).
368 */
369 template <class Base>
370 template <class VectorSet>
372  const std::set<size_t>& set_type ,
373  bool transpose ,
374  bool dependency ,
375  size_t q ,
376  const VectorSet& r ,
377  VectorSet& s )
378 { size_t m = Range();
379  size_t n = Domain();
380
381  // check VectorSet is Simple Vector class with sets for elements
382  CheckSimpleVector<std::set<size_t>, VectorSet>(
383  local::one_element_std_set<size_t>(), local::two_element_std_set<size_t>()
384  );
385
386  // dimension size of result vector
387  if( transpose )
388  s.resize(q);
389  else s.resize( m );
390
391  // temporary iterator
392  std::set<size_t>::const_iterator itr_1;
393
395  q > 0,
396  "ForSparseJac: q is not greater than zero"
397  );
399  size_t(r.size()) == n || transpose,
400  "ForSparseJac: size of r is not equal to n and transpose is false."
401  );
403  size_t(r.size()) == q || ! transpose,
404  "ForSparseJac: size of r is not equal to q and transpose is true."
405  );
406  //
407  // allocate memory for the requested sparsity calculation
408  for_jac_sparse_set_.resize(num_var_tape_, q);
409
410  // set values corresponding to independent variables
411  if( transpose )
412  { for(size_t i = 0; i < q; i++)
413  { // add the elements that are present
414  itr_1 = r[i].begin();
415  while( itr_1 != r[i].end() )
416  { size_t j = *itr_1++;
418  j < n,
419  "ForSparseJac: transpose is true and element of the set\n"
420  "r[j] has value greater than or equal n."
421  );
423  // operator for j-th independent variable
425  play_.GetOp( ind_taddr_[j] ) == local::InvOp
426  );
428  }
429  }
430  }
431  else
432  { for(size_t i = 0; i < n; i++)
436
437  // add the elements that are present
438  itr_1 = r[i].begin();
439  while( itr_1 != r[i].end() )
440  { size_t j = *itr_1++;
442  j < q,
443  "ForSparseJac: an element of the set r[i] "
444  "has value greater than or equal q."
445  );
447  }
448  }
449  }
450  // process posts
451  for(size_t j = 0; j < n; j++)
453
454  // evaluate the sparsity patterns
456  &play_,
457  dependency,
458  n,
459  num_var_tape_,
460  for_jac_sparse_set_
461  );
462
463  // return values corresponding to dependent variables
464  CPPAD_ASSERT_UNKNOWN( size_t(s.size()) == m || transpose );
465  CPPAD_ASSERT_UNKNOWN( size_t(s.size()) == q || ! transpose );
466  for(size_t i = 0; i < m; i++)
468
469  // extract results from for_jac_sparse_set_
470  // and add corresponding elements to sets in s
471  CPPAD_ASSERT_UNKNOWN( for_jac_sparse_set_.end() == q );
474  size_t j = *itr_2;
475  while( j < q )
476  { if( transpose )
477  s[j].insert(i);
478  else s[i].insert(j);
479  j = *(++itr_2);
480  }
481  }
482 }
483 // ---------------------------------------------------------------------------
484
485 /*!
486 User API for Jacobian sparsity patterns using forward mode.
487
488 The C++ source code corresponding to this operation is
489 \verbatim
490  s = f.ForSparseJac(q, r, transpose, dependency)
491 \endverbatim
492
493 \tparam Base
494 is the base type for this recording.
495
496 \tparam VectorSet
497 is a simple vector with elements of type \c bool
498 or \c std::set<size_t>.
499
500 \param q
501 is the number of columns in the matrix \f$R \f$.
502
503 \param r
504 is a sparsity pattern for the matrix \f$R \f$.
505
506 \param transpose
507 are sparsity patterns for \f$R \f$ and \f$S(x) \f$ transposed.
508
509 \param dependency
510 Are the derivatives with respect to left and right of the expression below
511 considered to be non-zero:
512 \code
513  CondExpRel(left, right, if_true, if_false)
514 \endcode
515 This is used by the optimizer to obtain the correct dependency relations.
516
517 \return
518 The value of \c transpose is false (true),
519 the return value is a sparsity pattern for \f$S(x) \f$ (\f$S(x)^T \f$) where
520 \f[
521  S(x) = F^{(1)} (x) * R
522 \f]
523 where \f$F \f$ is the function corresponding to the operation sequence
524 and \a x is any argument value.
525 If \c VectorSet::value_type is \c bool,
526 the return value has size \f$m * q \f$ (\f$q * m \f$).
527 where \c m is the number of dependent variables
528 corresponding to the operation sequence stored in \c f.
529 If \c VectorSet::value_type is \c std::set<size_t>,
530 the return value has size \f$m \f$ ( \f$q \f$ )
531 and with all its elements between zero and
532 \f$q - 1 \f$ ( \f$m - 1 \f$).
533
534 \par Side Effects
535 If \c VectorSet::value_type is \c bool,
536 the forward sparsity pattern for all of the variables on the
537 tape is stored in \c for_jac_sparse_pack__.
538 In this case
539 \verbatim
540  for_jac_sparse_pack_.n_set() == num_var_tape_
541  for_jac_sparse_pack_.end() == q
542  for_jac_sparse_set_.n_set() == 0
543  for_jac_sparse_set_.end() == 0
544 \endverbatim
545 \n
546 \n
547 If \c VectorSet::value_type is \c std::set<size_t>,
548 the forward sparsity pattern for all of the variables on the
549 tape is stored in \c for_jac_sparse_set__.
550 In this case
551 \verbatim
552  for_jac_sparse_set_.n_set() == num_var_tape_
553  for_jac_sparse_set_.end() == q
554  for_jac_sparse_pack_.n_set() == 0
555  for_jac_sparse_pack_.end() == 0
556 \endverbatim
557 */
558 template <class Base>
559 template <class VectorSet>
561  size_t q ,
562  const VectorSet& r ,
563  bool transpose ,
564  bool dependency )
565 { VectorSet s;
566  typedef typename VectorSet::value_type Set_type;
567
568  // free all memory currently in sparsity patterns
569  for_jac_sparse_pack_.resize(0, 0);
570  for_jac_sparse_set_.resize(0, 0);
571
572  ForSparseJacCase(
573  Set_type() ,
574  transpose ,
575  dependency ,
576  q ,
577  r ,
578  s
579  );
580
581  return s;
582 }
583 // ===========================================================================
584 // ForSparseJacCheckpoint
585 /*!
586 Forward mode Jacobian sparsity calculation used by checkpoint functions.
587
588 \tparam Base
589 is the base type for this recording.
590
591 \param transpose
592 is true (false) s is equal to \f$S(x) \f$ (\f$S(x)^T \f$)
593 where
594 \f[
595  S(x) = F^{(1)} (x) * R
596 \f]
597 where \f$F \f$ is the function corresponding to the operation sequence
598 and \f$x \f$ is any argument value.
599
600 \param q
601 is the number of columns in the matrix \f$R \f$.
602
603 \param r
604 is a sparsity pattern for the matrix \f$R \f$.
605
606 \param transpose
607 are the sparsity patterns for \f$R \f$ and \f$S(x) \f$ transposed.
608
609 \param dependency
610 Are the derivatives with respect to left and right of the expression below
611 considered to be non-zero:
612 \code
613  CondExpRel(left, right, if_true, if_false)
614 \endcode
615 This is used by the optimizer to obtain the correct dependency relations.
616
617 \param s
618 The input size and elements of s do not matter.
619 On output, s is the sparsity pattern for the matrix \f$S(x) \f$
620 or \f$S(x)^T \f$ depending on transpose.
621
622 \par Side Effects
623 If \c VectorSet::value_type is \c bool,
624 the forward sparsity pattern for all of the variables on the
625 tape is stored in \c for_jac_sparse_pack__.
626 In this case
627 \verbatim
628  for_jac_sparse_pack_.n_set() == num_var_tape_
629  for_jac_sparse_pack_.end() == q
630  for_jac_sparse_set_.n_set() == 0
631  for_jac_sparse_set_.end() == 0
632 \endverbatim
633 \n
634 \n
635 If \c VectorSet::value_type is \c std::set<size_t>,
636 the forward sparsity pattern for all of the variables on the
637 tape is stored in \c for_jac_sparse_set__.
638 In this case
639 \verbatim
640  for_jac_sparse_set_.n_set() == num_var_tape_
641  for_jac_sparse_set_.end() == q
642  for_jac_sparse_pack_.n_set() == 0
643  for_jac_sparse_pack_.end() == 0
644 \endverbatim
645 */
646 template <class Base>
648  size_t q ,
649  const local::sparse_list& r ,
650  bool transpose ,
651  bool dependency ,
652  local::sparse_list& s )
653 { size_t n = Domain();
654  size_t m = Range();
655
656 # ifndef NDEBUG
657  if( transpose )
658  { CPPAD_ASSERT_UNKNOWN( r.n_set() == q );
659  CPPAD_ASSERT_UNKNOWN( r.end() == n );
660  }
661  else
662  { CPPAD_ASSERT_UNKNOWN( r.n_set() == n );
663  CPPAD_ASSERT_UNKNOWN( r.end() == q );
664  }
665  for(size_t j = 0; j < n; j++)
668  }
669 # endif
670
671  // free all memory currently in sparsity patterns
672  for_jac_sparse_pack_.resize(0, 0);
673  for_jac_sparse_set_.resize(0, 0);
674
675  // allocate new sparsity pattern
676  for_jac_sparse_set_.resize(num_var_tape_, q);
677
678  // set sparsity pattern for dependent variables
679  if( transpose )
680  { for(size_t i = 0; i < q; i++)
682  size_t j = *itr;
683  while( j < n )
684  { for_jac_sparse_set_.post_element( ind_taddr_[j], i );
685  j = *(++itr);
686  }
687  }
688  }
689  else
690  { for(size_t j = 0; j < n; j++)
692  size_t i = *itr;
693  while( i < q )
694  { for_jac_sparse_set_.post_element( ind_taddr_[j], i );
695  i = *(++itr);
696  }
697  }
698  }
699  // process posts
700  for(size_t j = 0; j < n; j++)
702
703  // evaluate the sparsity pattern for all variables
705  &play_,
706  dependency,
707  n,
708  num_var_tape_,
709  for_jac_sparse_set_
710  );
711
712  // dimension the return value
713  if( transpose )
714  s.resize(q, m);
715  else
716  s.resize(m, q);
717
718  // return values corresponding to dependent variables
719  for(size_t i = 0; i < m; i++)
721
722  // extract the result from for_jac_sparse_set_
723  CPPAD_ASSERT_UNKNOWN( for_jac_sparse_set_.end() == q );
726  size_t j = *itr;
727  while( j < q )
728  { if( transpose )
729  s.post_element(j, i);
730  else
731  s.post_element(i, j);
732  j = *(++itr);
733  }
734  }
735  // process posts
736  for(size_t i = 0; i < s.n_set(); ++i)
737  s.process_post(i);
738
739 }
740
741
743 # 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 ForSparseJacCheckpoint(size_t q, const local::sparse_list &r, bool transpose, bool dependency, local::sparse_list &s)
Forward mode Jacobian sparsity calculation used by checkpoint functions.
void ForSparseJacCase(bool set_type, bool transpose, bool dependency, size_t q, const VectorSet &r, VectorSet &s)
Private helper function for ForSparseJac(q, r) boolean sparsity patterns.
void post_element(size_t i, size_t element)
Post an element for delayed addition to a set.
Two constant standard sets (currently used for concept checking).
Check that exp is true, if not terminate execution.
VectorSet ForSparseJac(size_t q, const VectorSet &r, bool transpose=false, bool dependency=false)
User API for Jacobian sparsity patterns using forward mode.
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.
Scalar value_type
void for_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 independent variables, ForJacSweep computes the sparsity pattern f...