CppAD: A C++ Algorithmic Differentiation Package  20171217
rev_hes_sweep.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
16 /*!
17 \file rev_hes_sweep.hpp
18 Compute Reverse mode Hessian sparsity patterns.
19 */
20
21 /*!
23 This value is either zero or one.
24 Zero is the normal operational value.
25 If it is one, a trace of every rev_hes_sweep computation is printed.
26 */
28
29 /*!
30 Given the forward Jacobian sparsity pattern for all the variables,
31 and the reverse Jacobian sparsity pattern for the dependent variables,
32 RevHesSweep computes the Hessian sparsity pattern for all the independent
33 variables.
34
35 \tparam Base
36 this operation sequence was recorded using AD<Base>.
37
38 \tparam Vector_set
39 is the type used for vectors of sets. It can be either
40 sparse_pack or sparse_list.
41
42 \param n
43 is the number of independent variables on the tape.
44
45 \param numvar
46 is the total number of variables on the tape; i.e.,
47 \a play->num_var_rec().
48 This is also the number of rows in the entire sparsity pattern
49 \a rev_hes_sparse.
50
51 \param play
52 The information stored in \a play
53 is a recording of the operations corresponding to a function
54 \f[
55  F : {\bf R}^n \rightarrow {\bf R}^m
56 \f]
57 where \f$n \f$ is the number of independent variables
58 and \f$m \f$ is the number of dependent variables.
59
60 \param for_jac_sparse
61 For i = 0 , ... , \a numvar - 1,
62 (for all the variables on the tape),
63 the forward Jacobian sparsity pattern for the variable with index i
64 corresponds to the set with index i in \a for_jac_sparse.
65
66 \param RevJac
67 \b Input:
68 For i = 0, ... , \a numvar - 1
69 the if the variable with index i on the tape is an dependent variable and
70 included in the Hessian, \a RevJac[ i ] is equal to true,
71 otherwise it is equal to false.
72 \n
73 \n
74 \b Output: The values in \a RevJac upon return are not specified; i.e.,
75 it is used for temporary work space.
76
77 \param rev_hes_sparse
78 The reverse Hessian sparsity pattern for the variable with index i
79 corresponds to the set with index i in \a rev_hes_sparse.
80 \n
81 \n
82 \b Input: For i = 0 , ... , \a numvar - 1
83 the reverse Hessian sparsity pattern for the variable with index i is empty.
84 \n
85 \n
86 \b Output: For j = 1 , ... , \a n,
87 the reverse Hessian sparsity pattern for the independent dependent variable
88 with index (j-1) is given by the set with index j
89 in \a rev_hes_sparse.
90 The values in the rest of \a rev_hes_sparse are not specified; i.e.,
91 they are used for temporary work space.
92 */
93
94 template <class Base, class Vector_set>
96  const local::player<Base>* play,
97  size_t n,
98  size_t numvar,
99  const Vector_set& for_jac_sparse,
100  bool* RevJac,
101  Vector_set& rev_hes_sparse
102 )
103 {
104  OpCode op;
105  size_t i_op;
106  size_t i_var;
107
109
110  // length of the parameter vector (used by CppAD assert macros)
111  const size_t num_par = play->num_par_rec();
112
113  size_t i, j, k;
114
115  // check numvar argument
116  CPPAD_ASSERT_UNKNOWN( play->num_var_rec() == numvar );
117  CPPAD_ASSERT_UNKNOWN( for_jac_sparse.n_set() == numvar );
118  CPPAD_ASSERT_UNKNOWN( rev_hes_sparse.n_set() == numvar );
119  CPPAD_ASSERT_UNKNOWN( numvar > 0 );
120
121  // upper limit exclusive for set elements
122  size_t limit = rev_hes_sparse.end();
123  CPPAD_ASSERT_UNKNOWN( for_jac_sparse.end() == limit );
124
125  // check number of sets match
127  for_jac_sparse.n_set() == rev_hes_sparse.n_set()
128  );
129
132  // to the index for the corresponding set in vecad_sparsity.
138  if( num_vecad_vec > 0 )
139  { size_t length;
143  j = 0;
144  for(i = 0; i < num_vecad_vec; i++)
145  { // length of this VecAD
146  length = play->GetVecInd(j);
149  // make all other values for this vector invalid
150  for(k = 1; k <= length; k++)
152  // start of next VecAD
153  j += length + 1;
154  // initialize this vector's reverse jacobian value
156  }
157  CPPAD_ASSERT_UNKNOWN( j == play->num_vec_ind_rec() );
158  }
159
160  // ----------------------------------------------------------------------
161  // user's atomic op calculator
162  atomic_base<Base>* user_atom = CPPAD_NULL; // user's atomic op calculator
163  //
164  // work space used by UserOp.
165  vector<Base> user_x; // parameters in x as integers
166  vector<size_t> user_ix; // variable indices for argument vector
167  vector<size_t> user_iy; // variable indices for result vector
168  //
169  // information set by forward_user (initialization to avoid warnings)
170  size_t user_old=0, user_m=0, user_n=0, user_i=0, user_j=0;
171  // information set by forward_user (necessary initialization)
172  enum_user_state user_state = end_user; // proper initialization
173  // ----------------------------------------------------------------------
174  //
175  // pointer to the beginning of the parameter vector
176  // (used by atomic functions
177  const Base* parameter = CPPAD_NULL;
178  if( num_par > 0 )
179  parameter = play->GetPar();
180  //
181  // Initialize
182  i_op = play->num_op_rec();
183  play->get_op_info(--i_op, op, arg, i_var);
184  CPPAD_ASSERT_UNKNOWN( op == EndOp );
186  std::cout << std::endl;
189 # endif
190  bool more_operators = true;
191  while(more_operators)
192  { bool flag; // temporary for use in switch cases
193  //
194  // next op
195  play->get_op_info(--i_op, op, arg, i_var);
196
197  // rest of information depends on the case
198  switch( op )
199  {
200  case AbsOp:
203  i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
204  );
205  break;
206  // -------------------------------------------------
207
211  i_var, arg, RevJac, for_jac_sparse, rev_hes_sparse
212  );
213  break;
214  // -------------------------------------------------
215
219  i_var, arg[1], RevJac, for_jac_sparse, rev_hes_sparse
220  );
221  break;
222  // -------------------------------------------------
223
224  case AcosOp:
225  // sqrt(1 - x * x), acos(x)
228  i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
229  );
230  break;
231  // -------------------------------------------------
232
234  case AcoshOp:
235  // sqrt(x * x - 1), acosh(x)
238  i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
239  );
240  break;
241 # endif
242  // -------------------------------------------------
243
244  case AsinOp:
245  // sqrt(1 - x * x), asin(x)
248  i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
249  );
250  break;
251  // -------------------------------------------------
252
254  case AsinhOp:
255  // sqrt(1 + x * x), asinh(x)
258  i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
259  );
260  break;
261 # endif
262  // -------------------------------------------------
263
264  case AtanOp:
265  // 1 + x * x, atan(x)
268  i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
269  );
270  break;
271  // -------------------------------------------------
272
274  case AtanhOp:
275  // 1 - x * x, atanh(x)
278  i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
279  );
280  break;
281 # endif
282  // -------------------------------------------------
283
284  case BeginOp:
286  more_operators = false;
287  break;
288  // -------------------------------------------------
289
290  case CSkipOp:
291  break;
292  // -------------------------------------------------
293
294  case CSumOp:
296  i_var, arg, RevJac, rev_hes_sparse
297  );
298  break;
299  // -------------------------------------------------
300
301  case CExpOp:
303  i_var, arg, num_par, RevJac, rev_hes_sparse
304  );
305  break;
306  // ---------------------------------------------------
307
308  case CosOp:
309  // sin(x), cos(x)
312  i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
313  );
314  break;
315  // ---------------------------------------------------
316
317  case CoshOp:
318  // sinh(x), cosh(x)
321  i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
322  );
323  break;
324  // -------------------------------------------------
325
326  case DisOp:
328  // derivativve is identically zero
329  break;
330  // -------------------------------------------------
331
332  case DivvvOp:
335  i_var, arg, RevJac, for_jac_sparse, rev_hes_sparse
336  );
337  break;
338  // -------------------------------------------------
339
340  case DivpvOp:
343  i_var, arg[1], RevJac, for_jac_sparse, rev_hes_sparse
344  );
345  break;
346  // -------------------------------------------------
347
348  case DivvpOp:
351  i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
352  );
353  break;
354  // -------------------------------------------------
355
356  case ErfOp:
357  // arg[1] is always the parameter 0
358  // arg[2] is always the parameter 2 / sqrt(pi)
361  i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
362  );
363  break;
364  // -------------------------------------------------
365
366  case ExpOp:
369  i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
370  );
371  break;
372  // -------------------------------------------------
373
375  case Expm1Op:
378  i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
379  );
380  break;
381 # endif
382  // -------------------------------------------------
383
384  case InvOp:
386  // Z is already defined
387  break;
388  // -------------------------------------------------
389
390  case LdpOp:
392  op,
393  i_var,
394  arg,
397  rev_hes_sparse,
399  RevJac,
401  );
402  break;
403  // -------------------------------------------------
404
405  case LdvOp:
407  op,
408  i_var,
409  arg,
412  rev_hes_sparse,
414  RevJac,
416  );
417  break;
418  // -------------------------------------------------
419
420  case EqpvOp:
421  case EqvvOp:
422  case LtpvOp:
423  case LtvpOp:
424  case LtvvOp:
425  case LepvOp:
426  case LevpOp:
427  case LevvOp:
428  case NepvOp:
429  case NevvOp:
431  break;
432  // -------------------------------------------------
433
434  case LogOp:
437  i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
438  );
439  break;
440  // -------------------------------------------------
441
443  case Log1pOp:
446  i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
447  );
448  break;
449 # endif
450  // -------------------------------------------------
451
452  case MulpvOp:
455  i_var, arg[1], RevJac, for_jac_sparse, rev_hes_sparse
456  );
457  break;
458  // -------------------------------------------------
459
460  case MulvvOp:
463  i_var, arg, RevJac, for_jac_sparse, rev_hes_sparse
464  );
465  break;
466  // -------------------------------------------------
467
468  case ParOp:
470
471  break;
472  // -------------------------------------------------
473
474  case PowpvOp:
477  i_var, arg[1], RevJac, for_jac_sparse, rev_hes_sparse
478  );
479  break;
480  // -------------------------------------------------
481
482  case PowvpOp:
485  i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
486  );
487  break;
488  // -------------------------------------------------
489
490  case PowvvOp:
493  i_var, arg, RevJac, for_jac_sparse, rev_hes_sparse
494  );
495  break;
496  // -------------------------------------------------
497
498  case PriOp:
500  break;
501  // -------------------------------------------------
502
503  case SignOp:
505  // Derivative is identiaclly zero
506  break;
507  // -------------------------------------------------
508
509  case SinOp:
510  // cos(x), sin(x)
513  i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
514  );
515  break;
516  // -------------------------------------------------
517
518  case SinhOp:
519  // cosh(x), sinh(x)
522  i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
523  );
524  break;
525  // -------------------------------------------------
526
527  case SqrtOp:
530  i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
531  );
532  break;
533  // -------------------------------------------------
534
535  case StppOp:
536  // sparsity cannot propagate through a parameter
538  break;
539  // -------------------------------------------------
540
541  case StpvOp:
543  op,
544  arg,
547  rev_hes_sparse,
549  RevJac,
551  );
552  break;
553  // -------------------------------------------------
554
555  case StvpOp:
556  // sparsity cannot propagate through a parameter
558  break;
559  // -------------------------------------------------
560
561  case StvvOp:
563  op,
564  arg,
567  rev_hes_sparse,
569  RevJac,
571  );
572  break;
573  // -------------------------------------------------
574
575  case SubvvOp:
578  i_var, arg, RevJac, for_jac_sparse, rev_hes_sparse
579  );
580  break;
581  // -------------------------------------------------
582
583  case SubpvOp:
586  i_var, arg[1], RevJac, for_jac_sparse, rev_hes_sparse
587  );
588  break;
589  // -------------------------------------------------
590
591  case SubvpOp:
594  i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
595  );
596  break;
597  // -------------------------------------------------
598
599  case TanOp:
600  // tan(x)^2, tan(x)
603  i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
604  );
605  break;
606  // -------------------------------------------------
607
608  case TanhOp:
609  // tanh(x)^2, tanh(x)
612  i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
613  );
614  break;
615  // -------------------------------------------------
616
617  case UserOp:
618  // start or end an atomic function call
620  user_state == start_user || user_state == end_user
621  );
622  flag = user_state == end_user;
623  user_atom = play->get_user_info(op, arg, user_old, user_m, user_n);
624  if( flag )
625  { user_state = ret_user;
626  user_i = user_m;
627  user_j = user_n;
628  //
629  user_x.resize(user_n);
630  user_ix.resize(user_n);
631  user_iy.resize(user_m);
632  }
633  else
634  { user_state = end_user;
635  //
636  // call users function for this operation
637  user_atom->set_old(user_old);
638  user_atom->rev_sparse_hes(
639  user_x, user_ix, user_iy,
640  for_jac_sparse, RevJac, rev_hes_sparse
641  );
642  }
643  break;
644
645  case UsrapOp:
646  // parameter argument in an atomic operation sequence
647  CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
648  CPPAD_ASSERT_UNKNOWN( user_state == arg_user );
649  CPPAD_ASSERT_UNKNOWN( user_i == 0 );
650  CPPAD_ASSERT_UNKNOWN( user_j <= user_n );
651  CPPAD_ASSERT_UNKNOWN( size_t( arg[0] ) < num_par );
652  //
653  --user_j;
654  // argument parameter value
655  user_x[user_j] = parameter[arg[0]];
656  // special variable index used for parameters
657  user_ix[user_j] = 0;
658  //
659  if( user_j == 0 )
660  user_state = start_user;
661  break;
662
663  case UsravOp:
664  // variable argument in an atomic operation sequence
665  CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
666  CPPAD_ASSERT_UNKNOWN( user_state == arg_user );
667  CPPAD_ASSERT_UNKNOWN( user_i == 0 );
668  CPPAD_ASSERT_UNKNOWN( user_j <= user_n );
669  //
670  --user_j;
671  // argument variables not available during sparsity calculations
673  // variable index for this argument
674  user_ix[user_j] = arg[0];
675  //
676  if( user_j == 0 )
677  user_state = start_user;
678  break;
679
680  case UsrrpOp:
681  // parameter result for a user atomic function
683  CPPAD_ASSERT_UNKNOWN( user_state == ret_user );
684  CPPAD_ASSERT_UNKNOWN( user_i <= user_m );
685  CPPAD_ASSERT_UNKNOWN( user_j == user_n );
686  CPPAD_ASSERT_UNKNOWN( size_t( arg[0] ) < num_par );
687  //
688  --user_i;
689  user_iy[user_i] = 0; // special variable used for parameters
690  //
691  if( user_i == 0 )
692  user_state = arg_user;
693  break;
694
695  case UsrrvOp:
696  // variable result for a user atomic function
698  CPPAD_ASSERT_UNKNOWN( user_state == ret_user );
699  CPPAD_ASSERT_UNKNOWN( user_i <= user_m );
700  CPPAD_ASSERT_UNKNOWN( user_j == user_n );
701  //
702  --user_i;
703  user_iy[user_i] = i_var; // variable for this result
704  //
705  if( user_i == 0 )
706  user_state = arg_user;
707  break;
708  // -------------------------------------------------
709
710  case ZmulpvOp:
713  i_var, arg[1], RevJac, for_jac_sparse, rev_hes_sparse
714  );
715  break;
716  // -------------------------------------------------
717
718  case ZmulvpOp:
721  i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse
722  );
723  break;
724  // -------------------------------------------------
725
726  case ZmulvvOp:
729  i_var, arg, RevJac, for_jac_sparse, rev_hes_sparse
730  );
731  break;
732
733  // -------------------------------------------------
734
735  default:
737  }
739  for(j = 0; j < limit; j++)
740  { zf_value[j] = false;
741  zh_value[j] = false;
742  }
743  typename Vector_set::const_iterator itr_jac(for_jac_sparse, i_var);
744  j = *itr_jac;
745  while( j < limit )
746  { zf_value[j] = true;
747  j = *(++itr_jac);
748  }
749  typename Vector_set::const_iterator itr_hes(rev_hes_sparse, i_var);
750  j = *itr_hes;
751  while( j < limit )
752  { zh_value[j] = true;
753  j = *(++itr_hes);
754  }
755  printOp(
756  std::cout,
757  play,
758  i_op,
759  i_var,
760  op,
761  arg
762  );
763  // should also print RevJac[i_var], but printOpResult does not
764  // yet allow for this
765  if( NumRes(op) > 0 && op != BeginOp ) printOpResult(
766  std::cout,
767  1,
768  &zf_value,
769  1,
770  &zh_value
771  );
772  std::cout << std::endl;
773  }
774  std::cout << std::endl;
775 # else
776  }
777 # endif
778  // values corresponding to BeginOp
779  CPPAD_ASSERT_UNKNOWN( i_op == 0 );
780  CPPAD_ASSERT_UNKNOWN( i_var == 0 );
781
782  return;
783 }
785
786 // preprocessor symbols that are local to this file
788
789 # endif
void resize(size_t n)
resize the vector (existing elements preserved when n &lt;= capacity_).
Definition: pod_vector.hpp:180
size_t num_par_rec(void) const
Fetch number of parameters in the recording.
Definition: player.hpp:638
size_t extend(size_t n)
Increase the number of elements the end of this vector (existing elements are always preserved)...
Definition: pod_vector.hpp:115
void printOpResult(std::ostream &os, size_t nfz, const Value *fz, size_t nrz, const Value *rz)
Prints the result values correspnding to an operator.
Definition: op_code.hpp:852
size_t num_vec_ind_rec(void) const
Fetch number of VecAD indices in the recording.
Definition: player.hpp:626
static Float quiet_NaN(void)
not a number
next UsrrpOp (UsrrvOp) is a parameter (variable) result
Definition: user_state.hpp:25
void reverse_sparse_hessian_addsub_op(size_t i_z, const addr_t *arg, bool *jac_reverse, const Vector_set &for_jac_sparsity, Vector_set &rev_hes_sparsity)
Reverse mode Hessian sparsity pattern for add and subtract operators.
size_t NumArg(OpCode op)
Number of arguments for a specified operator.
Definition: op_code.hpp:175
Class used to store and play back an operation sequence recording.
size_t GetVecInd(size_t i) const
Fetch a VecAD index from the recording.
Definition: player.hpp:571
void resize(size_t n)
change the number of elements in this vector.
Definition: vector.hpp:399
size_t NumRes(OpCode op)
Number of variables resulting from the specified operation.
Definition: op_code.hpp:281
next UserOp marks end of a user atomic call
Definition: user_state.hpp:28
atomic_base< Base > * get_user_info(const OpCode op, const addr_t *op_arg, size_t &user_old, size_t &user_m, size_t &user_n) const
unpack extra information corresponding to a UserOp
Definition: player.hpp:517
void reverse_sparse_hessian_nonlinear_unary_op(size_t i_z, size_t i_x, bool *rev_jacobian, const Vector_set &for_jac_sparsity, Vector_set &rev_hes_sparsity)
Reverse mode Hessian sparsity pattern for non-linear unary operators.
void rev_hes_sweep(const local::player< Base > *play, size_t n, size_t numvar, const Vector_set &for_jac_sparse, bool *RevJac, Vector_set &rev_hes_sparse)
Given the forward Jacobian sparsity pattern for all the variables, and the reverse Jacobian sparsity ...
virtual void set_old(size_t id)
Set value of id (used by deprecated old_atomic class)
OpCode
Type used to distinguish different AD&lt; Base &gt; atomic operations.
Definition: op_code.hpp:49
next UserOp marks beginning of a user atomic call
Definition: user_state.hpp:19
void get_op_info(size_t op_index, OpCode &op, const addr_t *&op_arg, size_t &var_index) const
fetch the information corresponding to an operator
Definition: player.hpp:485
void reverse_sparse_hessian_cond_op(size_t i_z, const addr_t *arg, size_t num_par, bool *jac_reverse, Vector_set &hes_sparsity)
Compute reverse Hessian sparsity patterns for op = CExpOp.
Definition: cond_op.hpp:1279
Fetch number of VecAD vectors in the recording.
Definition: player.hpp:630
Type * data(void)
current data pointer is no longer valid after any of the following: extend, erase, operator=, and ~pod_vector. Take extreem care when using this function.
Definition: pod_vector.hpp:89
void reverse_sparse_hessian_linear_unary_op(size_t i_z, size_t i_x, bool *rev_jacobian, const Vector_set &for_jac_sparsity, Vector_set &rev_hes_sparsity)
Reverse mode Hessian sparsity pattern for linear unary operators.
next UsrapOp (UsravOp) is a parameter (variable) argument
Definition: user_state.hpp:22
void reverse_sparse_hessian_pow_op(size_t i_z, const addr_t *arg, bool *jac_reverse, const Vector_set &for_jac_sparsity, Vector_set &rev_hes_sparsity)
Reverse mode Hessian sparsity pattern for power function.
virtual bool rev_sparse_hes(const vector< bool > &vx, const vector< bool > &s, vector< bool > &t, size_t q, const vector< std::set< size_t > > &r, const vector< std::set< size_t > > &u, vector< std::set< size_t > > &v, const vector< Base > &x)
Link from reverse Hessian sparsity sweep to base_atomic.
Reverse mode Hessian sparsity operations for LdpOp and LdvOp.
void reverse_sparse_hessian_mul_op(size_t i_z, const addr_t *arg, bool *jac_reverse, const Vector_set &for_jac_sparsity, Vector_set &rev_hes_sparsity)
Reverse mode Hessian sparsity pattern for multiplication operator.
Check that exp is true, if not terminate execution.
void reverse_sparse_hessian_csum_op(size_t i_z, const addr_t *arg, bool *rev_jacobian, Vector_set &rev_hes_sparsity)
Reverse mode Hessian sparsity pattern for CSumOp operator.
Definition: csum_op.hpp:600
void reverse_sparse_hessian_div_op(size_t i_z, const addr_t *arg, bool *jac_reverse, const Vector_set &for_jac_sparsity, Vector_set &rev_hes_sparsity)
Reverse mode Hessian sparsity pattern for division operator.
size_t num_op_rec(void) const
Fetch number of operators in the recording.
Definition: player.hpp:622
void printOp(std::ostream &os, const local::player< Base > *play, size_t i_op, size_t i_var, OpCode op, const addr_t *ind)
Prints a single operator and its operands.
Definition: op_code.hpp:547