CppAD: A C++ Algorithmic Differentiation Package  20171217
csum_op.hpp
Go to the documentation of this file.
1 // \$Id: csum_op.hpp 3845 2016-11-19 01:50:47Z bradbell \$
4
5 /* --------------------------------------------------------------------------
7
9 the terms of the
10  Eclipse Public License Version 1.0.
11
12 A copy of this license is included in the COPYING file of this distribution.
14 -------------------------------------------------------------------------- */
15
17 /*!
18 \file csum_op.hpp
19 Forward, reverse and sparsity calculations for cummulative summation.
20 */
21
22 /*!
23 Compute forward mode Taylor coefficients for result of op = CsumOp.
24
25 This operation is
26 \verbatim
27  z = s + x(1) + ... + x(m) - y(1) - ... - y(n).
28 \endverbatim
29
30 \tparam Base
31 base type for the operator; i.e., this operation was recorded
32 using AD< \a Base > and computations by this routine are done using type
33 \a Base.
34
35 \param p
36 lowest order of the Taylor coefficient that we are computing.
37
38 \param q
39 highest order of the Taylor coefficient that we are computing.
40
41 \param i_z
42 variable index corresponding to the result for this operation;
43 i.e. the row index in \a taylor corresponding to z.
44
45 \param arg
46 \a arg[0]
47 is the number of addition variables in this cummulative summation; i.e.,
48 <tt>m</tt>.
49 \n
50 \a arg[1]
51 is the number of subtraction variables in this cummulative summation; i.e.,
52 \c m.
53 \n
54 <tt>parameter[ arg[2] ]</tt>
55 is the parameter value \c s in this cummunative summation.
56 \n
57 <tt>arg[2+i]</tt>
58 for <tt>i = 1 , ... , m</tt> is the variable index of <tt>x(i)</tt>.
59 \n
60 <tt>arg[2+arg[0]+i]</tt>
61 for <tt>i = 1 , ... , n</tt> is the variable index of <tt>y(i)</tt>.
62
63 \param num_par
64 is the number of parameters in \a parameter.
65
66 \param parameter
67 is the parameter vector for this operation sequence.
68
69 \param cap_order
70 number of colums in the matrix containing all the Taylor coefficients.
71
72 \param taylor
73 \b Input: <tt>taylor [ arg[2+i] * cap_order + k ]</tt>
74 for <tt>i = 1 , ... , m</tt>
75 and <tt>k = 0 , ... , q</tt>
76 is the k-th order Taylor coefficient corresponding to <tt>x(i)</tt>
77 \n
78 \b Input: <tt>taylor [ arg[2+m+i] * cap_order + k ]</tt>
79 for <tt>i = 1 , ... , n</tt>
80 and <tt>k = 0 , ... , q</tt>
81 is the k-th order Taylor coefficient corresponding to <tt>y(i)</tt>
82 \n
83 \b Input: <tt>taylor [ i_z * cap_order + k ]</tt>
84 for k = 0 , ... , p,
85 is the k-th order Taylor coefficient corresponding to z.
86 \n
87 \b Output: <tt>taylor [ i_z * cap_order + k ]</tt>
88 for k = p , ... , q,
89 is the \a k-th order Taylor coefficient corresponding to z.
90 */
91 template <class Base>
92 inline void forward_csum_op(
93  size_t p ,
94  size_t q ,
95  size_t i_z ,
97  size_t num_par ,
98  const Base* parameter ,
99  size_t cap_order ,
100  Base* taylor )
101 { Base zero(0);
102  size_t i, j, k;
103
104  // check assumptions
106  CPPAD_ASSERT_UNKNOWN( q < cap_order );
107  CPPAD_ASSERT_UNKNOWN( p <= q );
108  CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < num_par );
110  arg[0] + arg[1] == arg[ arg[0] + arg[1] + 3 ]
111  );
112
113  // Taylor coefficients corresponding to result
114  Base* z = taylor + i_z * cap_order;
115  for(k = p; k <= q; k++)
116  z[k] = zero;
117  if( p == 0 )
118  z[p] = parameter[ arg[2] ];
119  Base* x;
120  i = arg[0];
121  j = 2;
122  while(i--)
123  { CPPAD_ASSERT_UNKNOWN( size_t(arg[j+1]) < i_z );
124  x = taylor + arg[++j] * cap_order;
125  for(k = p; k <= q; k++)
126  z[k] += x[k];
127  }
128  i = arg[1];
129  while(i--)
130  { CPPAD_ASSERT_UNKNOWN( size_t(arg[j+1]) < i_z );
131  x = taylor + arg[++j] * cap_order;
132  for(k = p; k <= q; k++)
133  z[k] -= x[k];
134  }
135 }
136
137 /*!
138 Multiple direction forward mode Taylor coefficients for op = CsumOp.
139
140 This operation is
141 \verbatim
142  z = s + x(1) + ... + x(m) - y(1) - ... - y(n).
143 \endverbatim
144
145 \tparam Base
146 base type for the operator; i.e., this operation was recorded
147 using AD<Base> and computations by this routine are done using type
148 \a Base.
149
150 \param q
151 order ot the Taylor coefficients that we are computing.
152
153 \param r
154 number of directions for Taylor coefficients that we are computing.
155
156 \param i_z
157 variable index corresponding to the result for this operation;
158 i.e. the row index in \a taylor corresponding to z.
159
160 \param arg
161 \a arg[0]
162 is the number of addition variables in this cummulative summation; i.e.,
163 <tt>m</tt>.
164 \n
165 \a arg[1]
166 is the number of subtraction variables in this cummulative summation; i.e.,
167 \c m.
168 \n
169 <tt>parameter[ arg[2] ]</tt>
170 is the parameter value \c s in this cummunative summation.
171 \n
172 <tt>arg[2+i]</tt>
173 for <tt>i = 1 , ... , m</tt> is the variable index of <tt>x(i)</tt>.
174 \n
175 <tt>arg[2+arg[0]+i]</tt>
176 for <tt>i = 1 , ... , n</tt> is the variable index of <tt>y(i)</tt>.
177
178 \param num_par
179 is the number of parameters in \a parameter.
180
181 \param parameter
182 is the parameter vector for this operation sequence.
183
184 \param cap_order
185 number of colums in the matrix containing all the Taylor coefficients.
186
187 \param taylor
188 \b Input: <tt>taylor [ arg[2+i]*((cap_order-1)*r + 1) + 0 ]</tt>
189 for <tt>i = 1 , ... , m</tt>
190 is the 0-th order Taylor coefficient corresponding to <tt>x(i)</tt> and
191 <tt>taylor [ arg[2+i]*((cap_order-1)*r + 1) + (q-1)*r + ell + 1 ]</tt>
192 for <tt>i = 1 , ... , m</tt>,
193 <tt>ell = 0 , ... , r-1</tt>
194 is the q-th order Taylor coefficient corresponding to <tt>x(i)</tt>
195 and direction ell.
196 \n
197 \b Input: <tt>taylor [ arg[2+m+i]*((cap_order-1)*r + 1) + 0 ]</tt>
198 for <tt>i = 1 , ... , n</tt>
199 is the 0-th order Taylor coefficient corresponding to <tt>y(i)</tt> and
200 <tt>taylor [ arg[2+m+i]*((cap_order-1)*r + 1) + (q-1)*r + ell + 1 ]</tt>
201 for <tt>i = 1 , ... , n</tt>,
202 <tt>ell = 0 , ... , r-1</tt>
203 is the q-th order Taylor coefficient corresponding to <tt>y(i)</tt>
204 and direction ell.
205 \n
206 \b Output: <tt>taylor [ i_z*((cap_order-1)*r+1) + (q-1)*r + ell + 1 ]</tt>
207 is the \a q-th order Taylor coefficient corresponding to z
208 for direction <tt>ell = 0 , ... , r-1</tt>.
209 */
210 template <class Base>
212  size_t q ,
213  size_t r ,
214  size_t i_z ,
216  size_t num_par ,
217  const Base* parameter ,
218  size_t cap_order ,
219  Base* taylor )
220 { Base zero(0);
221  size_t i, j, ell;
222
223  // check assumptions
225  CPPAD_ASSERT_UNKNOWN( q < cap_order );
226  CPPAD_ASSERT_UNKNOWN( 0 < q );
227  CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < num_par );
229  arg[0] + arg[1] == arg[ arg[0] + arg[1] + 3 ]
230  );
231
232  // Taylor coefficients corresponding to result
233  size_t num_taylor_per_var = (cap_order-1) * r + 1;
234  size_t m = (q-1)*r + 1;
235  Base* z = taylor + i_z * num_taylor_per_var + m;
236  for(ell = 0; ell < r; ell++)
237  z[ell] = zero;
238  Base* x;
239  i = arg[0];
240  j = 2;
241  while(i--)
242  { CPPAD_ASSERT_UNKNOWN( size_t(arg[j+1]) < i_z );
243  x = taylor + arg[++j] * num_taylor_per_var + m;
244  for(ell = 0; ell < r; ell++)
245  z[ell] += x[ell];
246  }
247  i = arg[1];
248  while(i--)
249  { CPPAD_ASSERT_UNKNOWN( size_t(arg[j+1]) < i_z );
250  x = taylor + arg[++j] * num_taylor_per_var + m;
251  for(ell = 0; ell < r; ell++)
252  z[ell] -= x[ell];
253  }
254 }
255
256 /*!
257 Compute reverse mode Taylor coefficients for result of op = CsumOp.
258
259 This operation is
260 \verbatim
261  z = q + x(1) + ... + x(m) - y(1) - ... - y(n).
262  H(y, x, w, ...) = G[ z(x, y), y, x, w, ... ]
263 \endverbatim
264
265 \tparam Base
266 base type for the operator; i.e., this operation was recorded
267 using AD< \a Base > and computations by this routine are done using type
268 \a Base.
269
270 \param d
271 order the highest order Taylor coefficient that we are computing
272 the partial derivatives with respect to.
273
274 \param i_z
275 variable index corresponding to the result for this operation;
276 i.e. the row index in \a taylor corresponding to z.
277
278 \param arg
279 \a arg[0]
280 is the number of addition variables in this cummulative summation; i.e.,
281 <tt>m</tt>.
282 \n
283 \a arg[1]
284 is the number of subtraction variables in this cummulative summation; i.e.,
285 \c m.
286 \n
287 <tt>parameter[ arg[2] ]</tt>
288 is the parameter value \c q in this cummunative summation.
289 \n
290 <tt>arg[2+i]</tt>
291 for <tt>i = 1 , ... , m</tt> is the value <tt>x(i)</tt>.
292 \n
293 <tt>arg[2+arg[0]+i]</tt>
294 for <tt>i = 1 , ... , n</tt> is the value <tt>y(i)</tt>.
295
296 \param nc_partial
297 number of colums in the matrix containing all the partial derivatives.
298
299 \param partial
300 \b Input: <tt>partial [ arg[2+i] * nc_partial + k ]</tt>
301 for <tt>i = 1 , ... , m</tt>
302 and <tt>k = 0 , ... , d</tt>
303 is the partial derivative of G(z, y, x, w, ...) with respect to the
304 k-th order Taylor coefficient corresponding to <tt>x(i)</tt>
305 \n
306 \b Input: <tt>partial [ arg[2+m+i] * nc_partial + k ]</tt>
307 for <tt>i = 1 , ... , n</tt>
308 and <tt>k = 0 , ... , d</tt>
309 is the partial derivative of G(z, y, x, w, ...) with respect to the
310 k-th order Taylor coefficient corresponding to <tt>y(i)</tt>
311 \n
312 \b Input: <tt>partial [ i_z * nc_partial + k ]</tt>
313 for <tt>i = 1 , ... , n</tt>
314 and <tt>k = 0 , ... , d</tt>
315 is the partial derivative of G(z, y, x, w, ...) with respect to the
316 k-th order Taylor coefficient corresponding to \c z.
317 \n
318 \b Output: <tt>partial [ arg[2+i] * nc_partial + k ]</tt>
319 for <tt>i = 1 , ... , m</tt>
320 and <tt>k = 0 , ... , d</tt>
321 is the partial derivative of H(y, x, w, ...) with respect to the
322 k-th order Taylor coefficient corresponding to <tt>x(i)</tt>
323 \n
324 \b Output: <tt>partial [ arg[2+m+i] * nc_partial + k ]</tt>
325 for <tt>i = 1 , ... , n</tt>
326 and <tt>k = 0 , ... , d</tt>
327 is the partial derivative of H(y, x, w, ...) with respect to the
328 k-th order Taylor coefficient corresponding to <tt>y(i)</tt>
329 */
330
331 template <class Base>
332 inline void reverse_csum_op(
333  size_t d ,
334  size_t i_z ,
336  size_t nc_partial ,
337  Base* partial )
338 {
339  // check assumptions
341  CPPAD_ASSERT_UNKNOWN( d < nc_partial );
342
343  // Taylor coefficients and partial derivative corresponding to result
344  Base* pz = partial + i_z * nc_partial;
345  Base* px;
346  size_t i, j, k;
347  size_t d1 = d + 1;
348  i = arg[0];
349  j = 2;
350  while(i--)
351  { CPPAD_ASSERT_UNKNOWN( size_t(arg[j+1]) < i_z );
352  px = partial + arg[++j] * nc_partial;
353  k = d1;
354  while(k--)
355  px[k] += pz[k];
356  }
357  i = arg[1];
358  while(i--)
359  { CPPAD_ASSERT_UNKNOWN( size_t(arg[j+1]) < i_z );
360  px = partial + arg[++j] * nc_partial;
361  k = d1;
362  while(k--)
363  px[k] -= pz[k];
364  }
365 }
366
367
368 /*!
369 Forward mode Jacobian sparsity pattern for CSumOp operator.
370
371 This operation is
372 \verbatim
373  z = q + x(1) + ... + x(m) - y(1) - ... - y(n).
374 \endverbatim
375
376 \tparam Vector_set
377 is the type used for vectors of sets. It can be either
378 sparse_pack or sparse_list.
379
380 \param i_z
381 variable index corresponding to the result for this operation;
382 i.e. the index in \a sparsity corresponding to z.
383
384 \param arg
385 \a arg[0]
386 is the number of addition variables in this cummulative summation; i.e.,
387 <tt>m + n</tt>.
388 \n
389 \a arg[1]
390 is the number of subtraction variables in this cummulative summation; i.e.,
391 \c m.
392 \n
393 <tt>parameter[ arg[2] ]</tt>
394 is the parameter value \c q in this cummunative summation.
395 \n
396 <tt>arg[2+i]</tt>
397 for <tt>i = 1 , ... , m</tt> is the value <tt>x(i)</tt>.
398 \n
399 <tt>arg[2+arg[1]+i]</tt>
400 for <tt>i = 1 , ... , n</tt> is the value <tt>y(i)</tt>.
401
402 \param sparsity
403 \b Input:
404 For <tt>i = 1 , ... , m</tt>,
405 the set with index \a arg[2+i] in \a sparsity
406 is the sparsity bit pattern for <tt>x(i)</tt>.
407 This identifies which of the independent variables the variable
408 <tt>x(i)</tt> depends on.
409 \n
410 \b Input:
411 For <tt>i = 1 , ... , n</tt>,
412 the set with index \a arg[2+arg[0]+i] in \a sparsity
413 is the sparsity bit pattern for <tt>x(i)</tt>.
414 This identifies which of the independent variables the variable
415 <tt>y(i)</tt> depends on.
416 \n
417 \b Output:
418 The set with index \a i_z in \a sparsity
419 is the sparsity bit pattern for z.
420 This identifies which of the independent variables the variable z
421 depends on.
422 */
423
424 template <class Vector_set>
426  size_t i_z ,
428  Vector_set& sparsity )
429 { sparsity.clear(i_z);
430
431  size_t i, j;
432  i = arg[0] + arg[1];
433  j = 2;
434  while(i--)
435  { CPPAD_ASSERT_UNKNOWN( size_t(arg[j+1]) < i_z );
436  sparsity.binary_union(
437  i_z , // index in sparsity for result
438  i_z , // index in sparsity for left operand
439  arg[++j] , // index for right operand
440  sparsity // sparsity vector for right operand
441  );
442  }
443 }
444
445 /*!
446 Reverse mode Jacobian sparsity pattern for CSumOp operator.
447
448 This operation is
449 \verbatim
450  z = q + x(1) + ... + x(m) - y(1) - ... - y(n).
451  H(y, x, w, ...) = G[ z(x, y), y, x, w, ... ]
452 \endverbatim
453
454 \tparam Vector_set
455 is the type used for vectors of sets. It can be either
456 sparse_pack or sparse_list.
457
458 \param i_z
459 variable index corresponding to the result for this operation;
460 i.e. the index in \a sparsity corresponding to z.
461
462 \param arg
463 \a arg[0]
464 is the number of addition variables in this cummulative summation; i.e.,
465 <tt>m + n</tt>.
466 \n
467 \a arg[1]
468 is the number of subtraction variables in this cummulative summation; i.e.,
469 \c m.
470 \n
471 <tt>parameter[ arg[2] ]</tt>
472 is the parameter value \c q in this cummunative summation.
473 \n
474 <tt>arg[2+i]</tt>
475 for <tt>i = 1 , ... , m</tt> is the value <tt>x(i)</tt>.
476 \n
477 <tt>arg[2+arg[1]+i]</tt>
478 for <tt>i = 1 , ... , n</tt> is the value <tt>y(i)</tt>.
479
480 \param sparsity
481 For <tt>i = 1 , ... , m</tt>,
482 the set with index \a arg[2+i] in \a sparsity
483 is the sparsity bit pattern for <tt>x(i)</tt>.
484 This identifies which of the dependent variables depend on <tt>x(i)</tt>.
485 On input, the sparsity patter corresponds to \c G,
486 and on ouput it corresponds to \c H.
487 \n
488 For <tt>i = 1 , ... , m</tt>,
489 the set with index \a arg[2+arg[0]+i] in \a sparsity
490 is the sparsity bit pattern for <tt>y(i)</tt>.
491 This identifies which of the dependent variables depend on <tt>y(i)</tt>.
492 On input, the sparsity patter corresponds to \c G,
493 and on ouput it corresponds to \c H.
494 \n
495 \b Input:
496 The set with index \a i_z in \a sparsity
497 is the sparsity bit pattern for z.
498 On input it corresponds to \c G and on output it is undefined.
499 */
500
501 template <class Vector_set>
503  size_t i_z ,
505  Vector_set& sparsity )
506 {
507  size_t i, j;
508  i = arg[0] + arg[1];
509  j = 2;
510  while(i--)
511  { ++j;
512  CPPAD_ASSERT_UNKNOWN( size_t(arg[j]) < i_z );
513  sparsity.binary_union(
514  arg[j] , // index in sparsity for result
515  arg[j] , // index in sparsity for left operand
516  i_z , // index for right operand
517  sparsity // sparsity vector for right operand
518  );
519  }
520 }
521 /*!
522 Reverse mode Hessian sparsity pattern for CSumOp operator.
523
524 This operation is
525 \verbatim
526  z = q + x(1) + ... + x(m) - y(1) - ... - y(n).
527  H(y, x, w, ...) = G[ z(x, y), y, x, w, ... ]
528 \endverbatim
529
530 \tparam Vector_set
531 is the type used for vectors of sets. It can be either
532 sparse_pack or sparse_list.
533
534 \param i_z
535 variable index corresponding to the result for this operation;
536 i.e. the index in \a sparsity corresponding to z.
537
538 \param arg
539 \a arg[0]
540 is the number of addition variables in this cummulative summation; i.e.,
541 <tt>m + n</tt>.
542 \n
543 \a arg[1]
544 is the number of subtraction variables in this cummulative summation; i.e.,
545 \c m.
546 \n
547 <tt>parameter[ arg[2] ]</tt>
548 is the parameter value \c q in this cummunative summation.
549 \n
550 <tt>arg[2+i]</tt>
551 for <tt>i = 1 , ... , m</tt> is the value <tt>x(i)</tt>.
552 \n
553 <tt>arg[2+arg[0]+i]</tt>
554 for <tt>i = 1 , ... , n</tt> is the value <tt>y(i)</tt>.
555
556 \param rev_jacobian
557 <tt>rev_jacobian[i_z]</tt>
558 is all false (true) if the Jabobian of G with respect to z must be zero
559 (may be non-zero).
560 \n
561 \n
562 For <tt>i = 1 , ... , m</tt>
563 <tt>rev_jacobian[ arg[2+i] ]</tt>
564 is all false (true) if the Jacobian with respect to <tt>x(i)</tt>
565 is zero (may be non-zero).
566 On input, it corresponds to the function G,
567 and on output it corresponds to the function H.
568 \n
569 \n
570 For <tt>i = 1 , ... , n</tt>
571 <tt>rev_jacobian[ arg[2+arg[0]+i] ]</tt>
572 is all false (true) if the Jacobian with respect to <tt>y(i)</tt>
573 is zero (may be non-zero).
574 On input, it corresponds to the function G,
575 and on output it corresponds to the function H.
576
577 \param rev_hes_sparsity
578 The set with index \a i_z in in \a rev_hes_sparsity
579 is the Hessian sparsity pattern for the fucntion G
580 where one of the partials derivative is with respect to z.
581 \n
582 \n
583 For <tt>i = 1 , ... , m</tt>
584 The set with index <tt>arg[2+i]</tt> in \a rev_hes_sparsity
585 is the Hessian sparsity pattern
586 where one of the partials derivative is with respect to <tt>x(i)</tt>.
587 On input, it corresponds to the function G,
588 and on output it corresponds to the function H.
589 \n
590 \n
591 For <tt>i = 1 , ... , n</tt>
592 The set with index <tt>arg[2+arg[0]+i]</tt> in \a rev_hes_sparsity
593 is the Hessian sparsity pattern
594 where one of the partials derivative is with respect to <tt>y(i)</tt>.
595 On input, it corresponds to the function G,
596 and on output it corresponds to the function H.
597 */
598
599 template <class Vector_set>
601  size_t i_z ,
603  bool* rev_jacobian ,
604  Vector_set& rev_hes_sparsity )
605 {
606  size_t i, j;
607  i = arg[0] + arg[1];
608  j = 2;
609  while(i--)
610  { ++j;
611  CPPAD_ASSERT_UNKNOWN( size_t(arg[j]) < i_z );
612  rev_hes_sparsity.binary_union(
613  arg[j] , // index in sparsity for result
614  arg[j] , // index in sparsity for left operand
615  i_z , // index for right operand
616  rev_hes_sparsity // sparsity vector for right operand
617  );
618  rev_jacobian[arg[j]] |= rev_jacobian[i_z];
619  }
620 }
621
623 # endif
void forward_sparse_jacobian_csum_op(size_t i_z, const addr_t *arg, Vector_set &sparsity)
Forward mode Jacobian sparsity pattern for CSumOp operator.
Definition: csum_op.hpp:425
size_t NumRes(OpCode op)
Number of variables resulting from the specified operation.
Definition: op_code.hpp:281
void forward_csum_op_dir(size_t q, size_t r, size_t i_z, const addr_t *arg, size_t num_par, const Base *parameter, size_t cap_order, Base *taylor)
Multiple direction forward mode Taylor coefficients for op = CsumOp.
Definition: csum_op.hpp:211
void forward_csum_op(size_t p, size_t q, size_t i_z, const addr_t *arg, size_t num_par, const Base *parameter, size_t cap_order, Base *taylor)
Compute forward mode Taylor coefficients for result of op = CsumOp.
Definition: csum_op.hpp:92