CppAD: A C++ Algorithmic Differentiation Package  20171217
for_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 for_hes_sweep.hpp
18 Compute Forward 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 ForHesSweep 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 for_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 rev_jac_sparse
67 \b Input:
68 For i = 0, ... , \a numvar - 1
69 the if the function we are computing the Hessian for has a non-zero
70 derivative w.r.t. variable with index i,
71 the set with index i has element zero.
72 Otherwise it has no elements.
73
74 \param for_hes_sparse
75 The forward Hessian sparsity pattern for the variable with index i
76 corresponds to the set with index i in \a for_hes_sparse.
77 The number of rows in this sparsity patter is n+1 and the row
78 with index zero is not used.
79 \n
80 \n
81 \b Input: For i = 1 , ... , \a n
82 the forward Hessian sparsity pattern for the variable with index i is empty.
83 \n
84 \n
85 \b Output: For j = 1 , ... , \a n,
86 the forward Hessian sparsity pattern for the independent dependent variable
87 with index (j-1) is given by the set with index j
88 in \a for_hes_sparse.
89 */
90
91 template <class Base, class Vector_set>
93  const local::player<Base>* play,
94  size_t n,
95  size_t numvar,
96  const Vector_set& for_jac_sparse,
97  const Vector_set& rev_jac_sparse,
98  Vector_set& for_hes_sparse
99 )
100 {
101  OpCode op;
102  size_t i_op;
103  size_t i_var;
104
106
107  // length of the parameter vector (used by CppAD assert macros)
108  const size_t num_par = play->num_par_rec();
109
110  size_t i, j, k;
111
112  // check numvar argument
113  size_t limit = n+1;
114  CPPAD_ASSERT_UNKNOWN( play->num_var_rec() == numvar );
115  CPPAD_ASSERT_UNKNOWN( for_jac_sparse.n_set() == numvar );
116  CPPAD_ASSERT_UNKNOWN( for_hes_sparse.n_set() == limit );
117  CPPAD_ASSERT_UNKNOWN( numvar > 0 );
118
119  // upper limit exclusive for set elements
120  CPPAD_ASSERT_UNKNOWN( for_jac_sparse.end() == limit );
121  CPPAD_ASSERT_UNKNOWN( for_hes_sparse.end() == limit );
122
125  // to the index for the corresponding set in vecad_sparsity.
131  if( num_vecad_vec > 0 )
132  { size_t length;
136  j = 0;
137  for(i = 0; i < num_vecad_vec; i++)
138  { // length of this VecAD
139  length = play->GetVecInd(j);
142  // make all other values for this vector invalid
143  for(k = 1; k <= length; k++)
145  // start of next VecAD
146  j += length + 1;
147  // initialize this vector's reverse jacobian value
149  }
150  CPPAD_ASSERT_UNKNOWN( j == play->num_vec_ind_rec() );
151  }
152  // ------------------------------------------------------------------------
153  // user's atomic op calculator
154  atomic_base<Base>* user_atom = CPPAD_NULL; // user's atomic op calculator
155  //
156  // work space used by UserOp.
157  vector<Base> user_x; // value of parameter arguments to function
158  vector<size_t> user_ix; // variable index (on tape) for each argument
159  vector<size_t> user_iy; // variable index (on tape) for each result
160  //
161  // information set by forward_user (initialization to avoid warnings)
162  size_t user_old=0, user_m=0, user_n=0, user_i=0, user_j=0;
163  // information set by forward_user (necessary initialization)
164  enum_user_state user_state = start_user;
165  // -------------------------------------------------------------------------
166  //
167  // pointer to the beginning of the parameter vector
168  // (used by user atomic functions)
169  const Base* parameter = CPPAD_NULL;
170  if( num_par > 0 )
171  parameter = play->GetPar();
172
173  // Initialize
174  i_op = 0;
175  play->get_op_info(i_op, op, arg, i_var);
176  CPPAD_ASSERT_UNKNOWN( op == BeginOp );
177  bool more_operators = true;
179  vector<size_t> user_usrrp; // parameter index for UsrrpOp operators
180  std::cout << std::endl;
183 # endif
184  bool flag; // temporary for use in switch cases below
185  while(more_operators)
186  {
187  // next op
188  play->get_op_info(++i_op, op, arg, i_var);
189
190  // does the Hessian in question have a non-zero derivative
191  // with respect to this variable
192  bool include = NumRes(op) > 0;
193  if( include )
194  include = rev_jac_sparse.is_element(i_var, 0);
195  //
196  // operators to include even if derivative is zero
197  include |= op == EndOp;
198  include |= op == CSkipOp;
199  include |= op == CSumOp;
200  include |= op == UserOp;
201  include |= op == UsrapOp;
202  include |= op == UsravOp;
203  include |= op == UsrrpOp;
204  include |= op == UsrrvOp;
205  //
206  if( include ) switch( op )
207  { // operators that should not occurr
208  // case BeginOp
209  // -------------------------------------------------
210
211  // operators that do not affect hessian
212  case AbsOp:
215  case CExpOp:
216  case DisOp:
217  case DivvpOp:
218  case InvOp:
219  case LdpOp:
220  case LdvOp:
221  case MulpvOp:
222  case ParOp:
223  case PriOp:
224  case SignOp:
225  case StppOp:
226  case StpvOp:
227  case StvpOp:
228  case StvvOp:
229  case SubvvOp:
230  case SubpvOp:
231  case SubvpOp:
232  case ZmulpvOp:
233  case ZmulvpOp:
234  break;
235  // -------------------------------------------------
236
237  // nonlinear unary operators
238  case AcosOp:
239  case AsinOp:
240  case AtanOp:
241  case CosOp:
242  case CoshOp:
243  case ExpOp:
244  case LogOp:
245  case SinOp:
246  case SinhOp:
247  case SqrtOp:
248  case TanOp:
249  case TanhOp:
251  case AcoshOp:
252  case AsinhOp:
253  case AtanhOp:
254  case Expm1Op:
255  case Log1pOp:
256 # endif
257  CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 )
259  arg[0], for_jac_sparse, for_hes_sparse
260  );
261  break;
262  // -------------------------------------------------
263
264  case CSkipOp:
265  break;
266  // -------------------------------------------------
267
268  case CSumOp:
269  break;
270  // -------------------------------------------------
271
272  case DivvvOp:
275  arg, for_jac_sparse, for_hes_sparse
276  );
277  break;
278  // -------------------------------------------------
279
280  case DivpvOp:
283  arg[1], for_jac_sparse, for_hes_sparse
284  );
285  break;
286  // -------------------------------------------------
287
288  case EndOp:
290  more_operators = false;
291  break;
292  // -------------------------------------------------
293
294  case ErfOp:
295  // arg[1] is always the parameter 0
296  // arg[2] is always the parameter 2 / sqrt(pi)
299  arg[0], for_jac_sparse, for_hes_sparse
300  );
301  break;
302  // -------------------------------------------------
303
304  // -------------------------------------------------
305  // logical comparision operators
306  case EqpvOp:
307  case EqvvOp:
308  case LtpvOp:
309  case LtvpOp:
310  case LtvvOp:
311  case LepvOp:
312  case LevpOp:
313  case LevvOp:
314  case NepvOp:
315  case NevvOp:
317  break;
318  // -------------------------------------------------
319
320  case MulvvOp:
323  arg, for_jac_sparse, for_hes_sparse
324  );
325  break;
326  // -------------------------------------------------
327
328  case PowpvOp:
331  arg[1], for_jac_sparse, for_hes_sparse
332  );
333  break;
334  // -------------------------------------------------
335
336  case PowvpOp:
339  arg[0], for_jac_sparse, for_hes_sparse
340  );
341  break;
342  // -------------------------------------------------
343
344  case PowvvOp:
347  arg, for_jac_sparse, for_hes_sparse
348  );
349  break;
350  // -------------------------------------------------
351
352  case UserOp:
353  // start or end an atomic function call
355  user_state == start_user || user_state == end_user
356  );
357  flag = user_state == start_user;
358  user_atom = play->get_user_info(op, arg, user_old, user_m, user_n);
359  if( flag )
360  { user_state = arg_user;
361  user_i = 0;
362  user_j = 0;
363  //
364  user_x.resize( user_n );
365  user_ix.resize( user_n );
366  user_iy.resize( user_m );
368  user_usrrp.resize( user_m );
369 # endif
370  }
371  else
372  { user_state = start_user;
373  //
374  user_atom->set_old(user_old);
375  user_atom->for_sparse_hes(
376  user_x, user_ix, user_iy,
377  for_jac_sparse, rev_jac_sparse, for_hes_sparse
378  );
379  }
380  break;
381
382  case UsrapOp:
383  // parameter argument for a user atomic function
384  CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
385  CPPAD_ASSERT_UNKNOWN( user_state == arg_user );
386  CPPAD_ASSERT_UNKNOWN( user_i == 0 );
387  CPPAD_ASSERT_UNKNOWN( user_j < user_n );
388  CPPAD_ASSERT_UNKNOWN( size_t( arg[0] ) < num_par );
389  //
390  user_x[user_j] = parameter[arg[0]];
391  user_ix[user_j] = 0; // special variable used for parameters
392  //
393  ++user_j;
394  if( user_j == user_n )
395  user_state = ret_user;
396  break;
397
398  case UsravOp:
399  // variable argument for a user atomic function
400  CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
401  CPPAD_ASSERT_UNKNOWN( user_state == arg_user );
402  CPPAD_ASSERT_UNKNOWN( user_i == 0 );
403  CPPAD_ASSERT_UNKNOWN( user_j < user_n );
404  //
405  // arguemnt variables not avaialbe during sparisty calculations
407  user_ix[user_j] = arg[0]; // variable for this argument
408  //
409  ++user_j;
410  if( user_j == user_n )
411  user_state = ret_user;
412  break;
413
414  case UsrrpOp:
415  // parameter result for a user atomic function
417  CPPAD_ASSERT_UNKNOWN( user_state == ret_user );
418  CPPAD_ASSERT_UNKNOWN( user_i < user_m );
419  CPPAD_ASSERT_UNKNOWN( user_j == user_n );
420  CPPAD_ASSERT_UNKNOWN( size_t( arg[0] ) < num_par );
421  //
422  user_iy[user_i] = 0; // special variable used for parameters
424  // remember argument for delayed tracing
425  user_usrrp[user_i] = arg[0];
426 # endif
427  ++user_i;
428  if( user_i == user_m )
429  user_state = end_user;
430  break;
431
432  case UsrrvOp:
433  // variable result for a user atomic function
435  CPPAD_ASSERT_UNKNOWN( user_state == ret_user );
436  CPPAD_ASSERT_UNKNOWN( user_i < user_m );
437  CPPAD_ASSERT_UNKNOWN( user_j == user_n );
438  //
439  user_iy[user_i] = i_var; // variable index for this result
440  //
441  ++user_i;
442  if( user_i == user_m )
443  user_state = end_user;
444  break;
445  // -------------------------------------------------
446
447  case ZmulvvOp:
450  arg, for_jac_sparse, for_hes_sparse
451  );
452  break;
453
454  // -------------------------------------------------
455
456  default:
458  }
460  typedef typename Vector_set::const_iterator const_iterator;
461  if( op == UserOp && user_state == start_user )
462  { // print operators that have been delayed
463  CPPAD_ASSERT_UNKNOWN( user_m == user_iy.size() );
464  CPPAD_ASSERT_UNKNOWN( i_op > user_m );
468  for(k = 0; k < user_m; k++)
469  { size_t k_var = user_iy[k];
470  // value for this variable
471  for(i = 0; i < limit; i++)
472  { zf_value[i] = false;
473  for(j = 0; j < limit; j++)
474  zh_value[i * limit + j] = false;
475  }
476  const_iterator itr_1(for_jac_sparse, i_var);
477  j = *itr_1;
478  while( j < limit )
479  { zf_value[j] = true;
480  j = *(++itr_1);
481  }
482  for(i = 0; i < limit; i++)
483  { const_iterator itr_2(for_hes_sparse, i);
484  j = *itr_2;
485  while( j < limit )
486  { zh_value[i * limit + j] = true;
487  j = *(++itr_2);
488  }
489  }
490  OpCode op_tmp = UsrrvOp;
491  if( k_var == 0 )
492  { op_tmp = UsrrpOp;
493  arg_tmp[0] = user_usrrp[k];
494  }
495  // k_var is zero when there is no result
496  printOp(
497  std::cout,
498  play,
499  i_op - user_m + k,
500  k_var,
501  op_tmp,
502  arg_tmp
503  );
504  if( k_var > 0 ) printOpResult(
505  std::cout,
506  1,
507  &zf_value,
508  1,
509  &zh_value
510  );
511  std::cout << std::endl;
512  }
513  }
514  for(i = 0; i < limit; i++)
515  { zf_value[i] = false;
516  for(j = 0; j < limit; j++)
517  zh_value[i * limit + j] = false;
518  }
519  const_iterator itr_1(for_jac_sparse, i_var);
520  j = *itr_1;
521  while( j < limit )
522  { zf_value[j] = true;
523  j = *(++itr_1);
524  }
525  for(i = 0; i < limit; i++)
526  { const_iterator itr_2(for_hes_sparse, i);
527  j = *itr_2;
528  while( j < limit )
529  { zh_value[i * limit + j] = true;
530  j = *(++itr_2);
531  }
532  }
533  // must delay print for these cases till after atomic user call
534  bool delay_print = op == UsrrpOp;
535  delay_print |= op == UsrrvOp;
536  if( ! delay_print )
537  { printOp(
538  std::cout,
539  play,
540  i_op,
541  i_var,
542  op,
543  arg
544  );
545  if( NumRes(op) > 0 && (! delay_print) ) printOpResult(
546  std::cout,
547  1,
548  &zf_value,
549  1,
550  &zh_value
551  );
552  std::cout << std::endl;
553  }
554  }
555  std::cout << std::endl;
556 # else
557  }
558 # endif
559
560  return;
561 }
563
564 // preprocessor symbols that are local to this file
566
567 # 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
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
virtual bool for_sparse_hes(const vector< bool > &vx, const vector< bool > &r, const vector< bool > &s, vector< std::set< size_t > > &h, const vector< Base > &x)
Link, after case split, from for_hes_sweep to atomic_base.
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 forward_sparse_hessian_mul_op(const addr_t *arg, const Vector_set &for_jac_sparsity, Vector_set &for_hes_sparsity)
Forward mode Hessian sparsity pattern for multiplication operator.
virtual void set_old(size_t id)
Set value of id (used by deprecated old_atomic class)
size_t size(void) const
number of elements currently in this vector.
Definition: vector.hpp:387
OpCode
Type used to distinguish different AD&lt; Base &gt; atomic operations.
Definition: op_code.hpp:49
void forward_sparse_hessian_pow_op(const addr_t *arg, const Vector_set &for_jac_sparsity, Vector_set &for_hes_sparsity)
Forward mode Hessian sparsity pattern for power operator.
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
Fetch number of VecAD vectors in the recording.
Definition: player.hpp:630
next UsrapOp (UsravOp) is a parameter (variable) argument
Definition: user_state.hpp:22
Check that exp is true, if not terminate execution.
void forward_sparse_hessian_nonlinear_unary_op(size_t i_v, const Vector_set &for_jac_sparsity, Vector_set &for_hes_sparsity)
Forward mode Hessian sparsity pattern for non-linear unary operators.
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