CppAD: A C++ Algorithmic Differentiation Package  20171217
erf_op.hpp
Go to the documentation of this file.
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
19
20
22 /*!
23 \file erf_op.hpp
24 Forward and reverse mode calculations for z = erf(x).
25 */
26
27 /*!
28 Forward mode Taylor coefficient for result of op = ErfOp.
29
30 The C++ source code corresponding to this operation is
31 \verbatim
32  z = erf(x)
33 \endverbatim
34
35 \tparam Base
36 base type for the operator; i.e., this operation was recorded
37 using AD< \a Base > and computations by this routine are done using type
38 \a Base.
39
40 \param p
41 lowest order of the Taylor coefficients that we are computing.
42
43 \param q
44 highest order of the Taylor coefficients that we are computing.
45
46 \param i_z
47 variable index corresponding to the last (primary) result for this operation;
48 i.e. the row index in \a taylor corresponding to z.
49 The auxillary results are called y_j have index \a i_z - j.
50
51 \param arg
52 arg[0]: is the variable index corresponding to x.
53 \n
54 arg[1]: is the parameter index corresponding to the value zero.
55 \n
56 \arg[2]: is the parameter index correspodning to the value 2 / sqrt(pi).
57
58 \param parameter
59 parameter[ arg[1] ] is the value zero,
60 and parameter[ arg[2] ] is the value 2 / sqrt(pi).
61
62 \param cap_order
63 maximum number of orders that will fit in the \c taylor array.
64
65 \param taylor
66 \b Input:
67 taylor [ arg[0] * cap_order + k ]
68 for k = 0 , ... , q,
69 is the k-th order Taylor coefficient corresponding to x.
70 \n
71 \b Input:
72 taylor [ i_z * cap_order + k ]
73 for k = 0 , ... , p - 1,
74 is the k-th order Taylor coefficient corresponding to z.
75 \n
76 \b Input:
77 taylor [ ( i_z - j) * cap_order + k ]
78 for k = 0 , ... , p-1,
79 and j = 0 , ... , 4,
80 is the k-th order Taylor coefficient corresponding to the j-th result for z.
81 \n
82 \b Output:
83 taylor [ (i_z-j) * cap_order + k ],
84 for k = p , ... , q,
85 and j = 0 , ... , 4,
86 is the k-th order Taylor coefficient corresponding to the j-th result for z.
87
88 \par Checked Assertions
89 \li NumArg(op) == 3
90 \li NumRes(op) == 5
91 \li q < cap_order
92 \li p <= q
93 \li std::numeric_limits<addr_t>::max() >= i_z + 2
94 */
95 template <class Base>
96 inline void forward_erf_op(
97  size_t p ,
98  size_t q ,
99  size_t i_z ,
101  const Base* parameter ,
102  size_t cap_order ,
103  Base* taylor )
104 {
105  // check assumptions
108  CPPAD_ASSERT_UNKNOWN( q < cap_order );
109  CPPAD_ASSERT_UNKNOWN( p <= q );
111
112  // array used to pass parameter values for sub-operations
114
115  // convert from final result to first result
116  i_z -= 4; // 4 = NumRes(ErfOp) - 1;
117
118  // z_0 = x * x
119  addr[0] = arg[0]; // x
120  addr[1] = arg[0]; // x
121  forward_mulvv_op(p, q, i_z+0, addr, parameter, cap_order, taylor);
122
123  // z_1 = - x * x
124  addr[0] = arg[1]; // zero
126  forward_subpv_op(p, q, i_z+1, addr, parameter, cap_order, taylor);
127
128  // z_2 = exp( - x * x )
129  forward_exp_op(p, q, i_z+2, i_z+1, cap_order, taylor);
130
131  // z_3 = (2 / sqrt(pi)) * exp( - x * x )
132  addr[0] = arg[2]; // 2 / sqrt(pi)
134  forward_mulpv_op(p, q, i_z+3, addr, parameter, cap_order, taylor);
135
136  // pointers to taylor coefficients for x , z_3, and z_4
137  Base* x = taylor + arg[0] * cap_order;
138  Base* z_3 = taylor + (i_z+3) * cap_order;
139  Base* z_4 = taylor + (i_z+4) * cap_order;
140
141  // calculte z_4 coefficients
142  if( p == 0 )
143  { // z4 (t) = erf[x(t)]
144  z_4[0] = erf(x[0]);
145  p++;
146  }
147  for(size_t j = p; j <= q; j++)
148  { // z_4' (t) = erf'[x(t)] * x'(t) = z3(t) * x'(t)
149  // z_4[1] + 2 * z_4[2] * t + ... =
150  // (z_3[0] + z_3[1] * t + ...) * (x[1] + 2 * x[2] * t + ...)
151  Base base_j = static_cast<Base>(double(j));
152  z_4[j] = static_cast<Base>(0);
153  for(size_t k = 1; k <= j; k++)
154  z_4[j] += (Base(double(k)) / base_j) * x[k] * z_3[j-k];
155  }
156 }
157
158 /*!
159 Zero order Forward mode Taylor coefficient for result of op = ErfOp.
160
161 The C++ source code corresponding to this operation is
162 \verbatim
163  z = erf(x)
164 \endverbatim
165
166 \tparam Base
167 base type for the operator; i.e., this operation was recorded
168 using AD< \a Base > and computations by this routine are done using type
169 \a Base.
170
171 \param i_z
172 variable index corresponding to the last (primary) result for this operation;
173 i.e. the row index in \a taylor corresponding to z.
174 The auxillary results are called y_j have index \a i_z - j.
175
176 \param arg
177 arg[0]: is the variable index corresponding to x.
178 \n
179 arg[1]: is the parameter index corresponding to the value zero.
180 \n
181 \arg[2]: is the parameter index correspodning to the value 2 / sqrt(pi).
182
183 \param parameter
184 parameter[ arg[1] ] is the value zero,
185 and parameter[ arg[2] ] is the value 2 / sqrt(pi).
186
187 \param cap_order
188 maximum number of orders that will fit in the \c taylor array.
189
190 \param taylor
191 \b Input:
192 taylor [ arg[0] * cap_order + 0 ]
193 is the zero order Taylor coefficient corresponding to x.
194 \n
195 \b Input:
196 taylor [ i_z * cap_order + 0 ]
197 is the zero order Taylor coefficient corresponding to z.
198 \n
199 \b Output:
200 taylor [ (i_z-j) * cap_order + 0 ],
201 for j = 0 , ... , 4,
202 is the zero order Taylor coefficient for j-th result corresponding to z.
203
204 \par Checked Assertions
205 \li NumArg(op) == 3
206 \li NumRes(op) == 5
207 \li q < cap_order
208 \li p <= q
209 \li std::numeric_limits<addr_t>::max() >= i_z + 2
210 */
211 template <class Base>
212 inline void forward_erf_op_0(
213  size_t i_z ,
215  const Base* parameter ,
216  size_t cap_order ,
217  Base* taylor )
218 {
219  // check assumptions
222  CPPAD_ASSERT_UNKNOWN( 0 < cap_order );
224
225  // array used to pass parameter values for sub-operations
227
228  // convert from final result to first result
229  i_z -= 4; // 4 = NumRes(ErfOp) - 1;
230
231  // z_0 = x * x
232  addr[0] = arg[0]; // x
233  addr[1] = arg[0]; // x
234  forward_mulvv_op_0(i_z+0, addr, parameter, cap_order, taylor);
235
236  // z_1 = - x * x
237  addr[0] = arg[1]; // zero
239  forward_subpv_op_0(i_z+1, addr, parameter, cap_order, taylor);
240
241  // z_2 = exp( - x * x )
242  forward_exp_op_0(i_z+2, i_z+1, cap_order, taylor);
243
244  // z_3 = (2 / sqrt(pi)) * exp( - x * x )
245  addr[0] = arg[2]; // 2 / sqrt(pi)
247  forward_mulpv_op_0(i_z+3, addr, parameter, cap_order, taylor);
248
249  // zero order Taylor coefficient for z_4
250  Base* x = taylor + arg[0] * cap_order;
251  Base* z_4 = taylor + (i_z + 4) * cap_order;
252  z_4[0] = erf(x[0]);
253 }
254 /*!
255 Forward mode Taylor coefficient for result of op = ErfOp.
256
257 The C++ source code corresponding to this operation is
258 \verbatim
259  z = erf(x)
260 \endverbatim
261
262 \tparam Base
263 base type for the operator; i.e., this operation was recorded
264 using AD< \a Base > and computations by this routine are done using type
265 \a Base.
266
267 \param q
268 order of the Taylor coefficients that we are computing.
269
270 \param r
271 number of directions for the Taylor coefficients that we afre computing.
272
273 \param i_z
274 variable index corresponding to the last (primary) result for this operation;
275 i.e. the row index in \a taylor corresponding to z.
276 The auxillary results have index i_z - j for j = 0 , ... , 4
277 (and include z).
278
279 \param arg
280 arg[0]: is the variable index corresponding to x.
281 \n
282 arg[1]: is the parameter index corresponding to the value zero.
283 \n
284 \arg[2]: is the parameter index correspodning to the value 2 / sqrt(pi).
285
286 \param parameter
287 parameter[ arg[1] ] is the value zero,
288 and parameter[ arg[2] ] is the value 2 / sqrt(pi).
289
290 \param cap_order
291 maximum number of orders that will fit in the \c taylor array.
292
293 \par tpv
294 We use the notation
295 <code>tpv = (cap_order-1) * r + 1</code>
296 which is the number of Taylor coefficients per variable
297
298 \param taylor
299 \b Input: If x is a variable,
300 <code>taylor [ arg[0] * tpv + 0 ]</code>,
301 is the zero order Taylor coefficient for all directions and
302 <code>taylor [ arg[0] * tpv + (k-1)*r + ell + 1 ]</code>,
303 for k = 1 , ... , q,
304 ell = 0, ..., r-1,
305 is the k-th order Taylor coefficient
306 corresponding to x and the ell-th direction.
307 \n
308 \b Input:
309 taylor [ (i_z - j) * tpv + 0 ]
310 is the zero order Taylor coefficient for all directions and the
311 j-th result for z.
312 for k = 1 , ... , q-1,
313 ell = 0, ... , r-1,
314 <code>
315 taylor[ (i_z - j) * tpv + (k-1)*r + ell + 1]
316 </code>
317 is the Taylor coefficient for the k-th order, ell-th direction,
318 and j-th auzillary result.
319 \n
320 \b Output:
321 taylor [ (i_z-j) * tpv + (q-1)*r + ell + 1 ],
322 for ell = 0 , ... , r-1,
323 is the Taylor coefficient for the q-th order, ell-th direction,
324 and j-th auzillary result.
325
326 \par Checked Assertions
327 \li NumArg(op) == 3
328 \li NumRes(op) == 5
329 \li 0 < q < cap_order
330 */
331 template <class Base>
332 inline void forward_erf_op_dir(
333  size_t q ,
334  size_t r ,
335  size_t i_z ,
337  const Base* parameter ,
338  size_t cap_order ,
339  Base* taylor )
340 {
341  // check assumptions
344  CPPAD_ASSERT_UNKNOWN( q < cap_order );
345  CPPAD_ASSERT_UNKNOWN( 0 < q );
347
348  // array used to pass parameter values for sub-operations
350
351  // convert from final result to first result
352  i_z -= 4; // 4 = NumRes(ErfOp) - 1;
353
354  // z_0 = x * x
355  addr[0] = arg[0]; // x
356  addr[1] = arg[0]; // x
357  forward_mulvv_op_dir(q, r, i_z+0, addr, parameter, cap_order, taylor);
358
359  // z_1 = - x * x
360  addr[0] = arg[1]; // zero
362  forward_subpv_op_dir(q, r, i_z+1, addr, parameter, cap_order, taylor);
363
364  // z_2 = exp( - x * x )
365  forward_exp_op_dir(q, r, i_z+2, i_z+1, cap_order, taylor);
366
367  // z_3 = (2 / sqrt(pi)) * exp( - x * x )
368  addr[0] = arg[2]; // 2 / sqrt(pi)
370  forward_mulpv_op_dir(q, r, i_z+3, addr, parameter, cap_order, taylor);
371
372  // pointers to taylor coefficients for x , z_3, and z_4
373  size_t num_taylor_per_var = (cap_order - 1) * r + 1;
374  Base* x = taylor + arg[0] * num_taylor_per_var;
375  Base* z_3 = taylor + (i_z+3) * num_taylor_per_var;
376  Base* z_4 = taylor + (i_z+4) * num_taylor_per_var;
377
378  // z_4' (t) = erf'[x(t)] * x'(t) = z3(t) * x'(t)
379  // z_4[1] + 2 * z_4[2] * t + ... =
380  // (z_3[0] + z_3[1] * t + ...) * (x[1] + 2 * x[2] * t + ...)
381  Base base_q = static_cast<Base>(double(q));
382  for(size_t ell = 0; ell < r; ell++)
383  { // index in z_4 and x for q-th order term
384  size_t m = (q-1)*r + ell + 1;
385  // initialize q-th order term summation
386  z_4[m] = z_3[0] * x[m];
387  for(size_t k = 1; k < q; k++)
388  { size_t x_index = (k-1)*r + ell + 1;
389  size_t z3_index = (q-k-1)*r + ell + 1;
390  z_4[m] += (Base(double(k)) / base_q) * x[x_index] * z_3[z3_index];
391  }
392  }
393 }
394
395 /*!
396 Compute reverse mode partial derivatives for result of op = ErfOp.
397
398 The C++ source code corresponding to this operation is
399 \verbatim
400  z = erf(x)
401 \endverbatim
402
403 \tparam Base
404 base type for the operator; i.e., this operation was recorded
405 using AD< \a Base > and computations by this routine are done using type
406 \a Base.
407
408 \param d
409 highest order Taylor of the Taylor coefficients that we are computing
410 the partial derivatives with respect to.
411
412 \param i_z
413 variable index corresponding to the last (primary) result for this operation;
414 i.e. the row index in \a taylor corresponding to z.
415 The auxillary results are called y_j have index \a i_z - j.
416
417 \param arg
418 arg[0]: is the variable index corresponding to x.
419 \n
420 arg[1]: is the parameter index corresponding to the value zero.
421 \n
422 \arg[2]: is the parameter index correspodning to the value 2 / sqrt(pi).
423
424 \param parameter
425 parameter[ arg[1] ] is the value zero,
426 and parameter[ arg[2] ] is the value 2 / sqrt(pi).
427
428 \param cap_order
429 maximum number of orders that will fit in the \c taylor array.
430
431 \param taylor
432 \b Input:
433 taylor [ arg[0] * cap_order + k ]
434 for k = 0 , ... , d,
435 is the k-th order Taylor coefficient corresponding to x.
436 \n
437 taylor [ (i_z - j) * cap_order + k ]
438 for k = 0 , ... , d,
439 and for j = 0 , ... , 4,
440 is the k-th order Taylor coefficient corresponding to the j-th result
441 for this operation.
442
443 \param nc_partial
444 number of columns in the matrix containing all the partial derivatives
445
446 \param partial
447 \b Input:
448 partial [ arg[0] * nc_partial + k ]
449 for k = 0 , ... , d,
450 is the partial derivative of G( z , x , w , u , ... ) with respect to
451 the k-th order Taylor coefficient for x.
452 \n
453 \b Input:
454 partial [ (i_z - j) * nc_partial + k ]
455 for k = 0 , ... , d,
456 and for j = 0 , ... , 4,
457 is the partial derivative of G( z , x , w , u , ... ) with respect to
458 the k-th order Taylor coefficient for the j-th result of this operation.
459 \n
460 \b Output:
461 partial [ arg[0] * nc_partial + k ]
462 for k = 0 , ... , d,
463 is the partial derivative of H( x , w , u , ... ) with respect to
464 the k-th order Taylor coefficient for x.
465 \n
466 \b Output:
467 partial [ (i_z-j) * nc_partial + k ]
468 for k = 0 , ... , d,
469 and for j = 0 , ... , 4,
470 may be used as work space; i.e., may change in an unspecified manner.
471
472 \par Checked Assertions
473 \li NumArg(op) == 3
474 \li NumRes(op) == 5
475 \li q < cap_order
476 \li p <= q
477 */
478 template <class Base>
479 inline void reverse_erf_op(
480  size_t d ,
481  size_t i_z ,
483  const Base* parameter ,
484  size_t cap_order ,
485  const Base* taylor ,
486  size_t nc_partial ,
487  Base* partial )
488 {
489  // check assumptions
492  CPPAD_ASSERT_UNKNOWN( d < cap_order );
494
495  // array used to pass parameter values for sub-operations
497
498  // If pz is zero, make sure this operation has no effect
499  // (zero times infinity or nan would be non-zero).
500  Base* pz = partial + i_z * nc_partial;
501  bool skip(true);
502  for(size_t i_d = 0; i_d <= d; i_d++)
503  skip &= IdenticalZero(pz[i_d]);
504  if( skip )
505  return;
506
507  // convert from final result to first result
508  i_z -= 4; // 4 = NumRes(ErfOp) - 1;
509
510  // Taylor coefficients and partials corresponding to x
511  const Base* x = taylor + arg[0] * cap_order;
512  Base* px = partial + arg[0] * nc_partial;
513
514  // Taylor coefficients and partials corresponding to z_3
515  const Base* z_3 = taylor + (i_z+3) * cap_order;
516  Base* pz_3 = partial + (i_z+3) * nc_partial;
517
518  // Taylor coefficients and partials corresponding to z_4
519  Base* pz_4 = partial + (i_z+4) * nc_partial;
520
521  // Reverse z_4
522  size_t j = d;
523  while(j)
524  { pz_4[j] /= Base(double(j));
525  for(size_t k = 1; k <= j; k++)
526  { px[k] += azmul(pz_4[j], z_3[j-k]) * Base(double(k));
527  pz_3[j-k] += azmul(pz_4[j], x[k]) * Base(double(k));
528  }
529  j--;
530  }
531  px[0] += azmul(pz_4[0], z_3[0]);
532
533  // z_3 = (2 / sqrt(pi)) * exp( - x * x )
534  addr[0] = arg[2]; // 2 / sqrt(pi)
537  d, i_z+3, addr, parameter, cap_order, taylor, nc_partial, partial
538  );
539
540  // z_2 = exp( - x * x )
542  d, i_z+2, i_z+1, cap_order, taylor, nc_partial, partial
543  );
544
545  // z_1 = - x * x
546  addr[0] = arg[1]; // zero
549  d, i_z+1, addr, parameter, cap_order, taylor, nc_partial, partial
550  );
551
552  // z_0 = x * x
553  addr[0] = arg[0]; // x
554  addr[1] = arg[0]; // x
556  d, i_z+0, addr, parameter, cap_order, taylor, nc_partial, partial
557  );
558
559 }
560
561
std::complex< double > erf(const std::complex< double > &x)
Definition: azmul.hpp:94
void reverse_mulvv_op(size_t d, size_t i_z, const addr_t *arg, const Base *parameter, size_t cap_order, const Base *taylor, size_t nc_partial, Base *partial)
Compute reverse mode partial derivatives for result of op = MulvvOp.
Definition: mul_op.hpp:158
void forward_subpv_op_dir(size_t q, size_t r, size_t i_z, const addr_t *arg, const Base *parameter, size_t cap_order, Base *taylor)
Multiple directions forward mode Taylor coefficients for op = SubpvOp.
Definition: sub_op.hpp:238
void forward_erf_op_0(size_t i_z, const addr_t *arg, const Base *parameter, size_t cap_order, Base *taylor)
Zero order Forward mode Taylor coefficient for result of op = ErfOp.
Definition: erf_op.hpp:212
size_t NumArg(OpCode op)
Number of arguments for a specified operator.
Definition: op_code.hpp:175
void reverse_mulpv_op(size_t d, size_t i_z, const addr_t *arg, const Base *parameter, size_t cap_order, const Base *taylor, size_t nc_partial, Base *partial)
Compute reverse mode partial derivative for result of op = MulpvOp.
Definition: mul_op.hpp:326
void reverse_exp_op(size_t d, size_t i_z, size_t i_x, size_t cap_order, const Base *taylor, size_t nc_partial, Base *partial)
Reverse mode partial derivatives for result of op = ExpOp.
Definition: exp_op.hpp:146
void forward_erf_op(size_t p, size_t q, size_t i_z, const addr_t *arg, const Base *parameter, size_t cap_order, Base *taylor)
Forward mode Taylor coefficient for result of op = ErfOp.
Definition: erf_op.hpp:96
size_t NumRes(OpCode op)
Number of variables resulting from the specified operation.
Definition: op_code.hpp:281
void reverse_erf_op(size_t d, size_t i_z, const addr_t *arg, const Base *parameter, size_t cap_order, const Base *taylor, size_t nc_partial, Base *partial)
Compute reverse mode partial derivatives for result of op = ErfOp.
Definition: erf_op.hpp:479
void forward_mulpv_op_0(size_t i_z, const addr_t *arg, const Base *parameter, size_t cap_order, Base *taylor)
Compute zero order forward mode Taylor coefficient for result of op = MulpvOp.
Definition: mul_op.hpp:291
void forward_erf_op_dir(size_t q, size_t r, size_t i_z, const addr_t *arg, const Base *parameter, size_t cap_order, Base *taylor)
Forward mode Taylor coefficient for result of op = ErfOp.
Definition: erf_op.hpp:332
Forward and reverse mode calculations for z = x - y.
Forward and reverse mode calculations for z = exp(x).
bool IdenticalZero(const std::complex< double > &x)
void forward_exp_op_dir(size_t q, size_t r, size_t i_z, size_t i_x, size_t cap_order, Base *taylor)
Multiple direction forward mode Taylor coefficient for op = ExpOp.
Definition: exp_op.hpp:78
void forward_subpv_op(size_t p, size_t q, size_t i_z, const addr_t *arg, const Base *parameter, size_t cap_order, Base *taylor)
Compute forward mode Taylor coefficients for result of op = SubpvOp.
Definition: sub_op.hpp:196
void forward_mulvv_op_dir(size_t q, size_t r, size_t i_z, const addr_t *arg, const Base *parameter, size_t cap_order, Base *taylor)
Multiple directions forward mode Taylor coefficients for op = MulvvOp.
Definition: mul_op.hpp:79
void forward_mulvv_op(size_t p, size_t q, size_t i_z, const addr_t *arg, const Base *parameter, size_t cap_order, Base *taylor)
Compute forward mode Taylor coefficients for result of op = MulvvOp.
Definition: mul_op.hpp:37
void forward_mulpv_op(size_t p, size_t q, size_t i_z, const addr_t *arg, const Base *parameter, size_t cap_order, Base *taylor)
Compute forward mode Taylor coefficients for result of op = MulpvOp.
Definition: mul_op.hpp:211
void forward_exp_op_0(size_t i_z, size_t i_x, size_t cap_order, Base *taylor)
Zero order forward mode Taylor coefficient for result of op = ExpOp.
Definition: exp_op.hpp:117
void reverse_subpv_op(size_t d, size_t i_z, const addr_t *arg, const Base *parameter, size_t cap_order, const Base *taylor, size_t nc_partial, Base *partial)
Compute reverse mode partial derivative for result of op = SubpvOp.
Definition: sub_op.hpp:312
Check that exp is true, if not terminate execution.
void forward_mulpv_op_dir(size_t q, size_t r, size_t i_z, const addr_t *arg, const Base *parameter, size_t cap_order, Base *taylor)
Multiple directions forward mode Taylor coefficients for op = MulpvOp.
Definition: mul_op.hpp:250
void forward_mulvv_op_0(size_t i_z, const addr_t *arg, const Base *parameter, size_t cap_order, Base *taylor)
Compute zero order forward mode Taylor coefficients for result of op = MulvvOp.
Definition: mul_op.hpp:124
void forward_exp_op(size_t p, size_t q, size_t i_z, size_t i_x, size_t cap_order, Base *taylor)
Forward mode Taylor coefficient for result of op = ExpOp.
Definition: exp_op.hpp:34
Forward and reverse mode calculations for z = x * y.
void forward_subpv_op_0(size_t i_z, const addr_t *arg, const Base *parameter, size_t cap_order, Base *taylor)
Compute zero order forward mode Taylor coefficient for result of op = SubpvOp.
Definition: sub_op.hpp:277