CppAD: A C++ Algorithmic Differentiation Package  20171217
abs_normal_fun.hpp
Go to the documentation of this file.
3 /* --------------------------------------------------------------------------
5
7 the terms of the
8  Eclipse Public License Version 1.0.
9
10 A copy of this license is included in the COPYING file of this distribution.
12 -------------------------------------------------------------------------- */
13 /*
14 $begin abs_normal_fun$$15 spell 16 const 17$$ 18 19 20$section Create An Abs-normal Representation of a Function$$21 22 head Syntax$$
23 $icode%f%.abs_normal_fun(%g%, %a%)%$$24 25 head f$$ 26 The object$icode f$$has prototype 27 codei% 28 const ADFun<%Base%>& %f% 29 %$$
30 It represents a function $latex f : \B{R}^n \rightarrow \B{R}^m$$. 31 We assume that the only non-smooth terms in the representation are 32 absolute value functions and use latex s \in \B{Z}_+$$ 33 to represent the number of these terms. 34 35$subhead n$$36 We use icode n$$ to denote the dimension of the domain space for $icode f$$. 37 38 subhead m$$ 39 We use$icode m$$to denote the dimension of the range space for icode f$$.
40
41 $subhead s$$42 We use icode s$$ to denote the number of absolute value terms in$icode f$$. 43 44 45 head a$$
46 The object $icode a$$has prototype 47 codei% 48 ADFun<%Base%> %a% 49 %$$ 50 The initial function representation in$icode a$$is lost. 51 Upon return it represents the result of the absolute terms 52 latex a : \B{R}^n \rightarrow \B{R}^s$$; see $latex a(x)$$defined below. 53 Note that icode a$$ is constructed by copying$icode f$$54 and then changing the dependent variables. There may 55 be many calculations in this representation that are not necessary 56 and can be removed using 57 codei% 58 %a%.optimize() 59 %$$
60 This optimization is not done automatically by $code abs_normal_fun$$61 because it may take a significant amount of time. 62 63 subhead zeta$$ 64 Let$latex \zeta_0 ( x )$$65 denote the argument for the first absolute value term in latex f(x)$$,
66 $latex \zeta_1 ( x , |\zeta_0 (x)| )$$for the second term, and so on. 67 68 subhead a(x)$$ 69 For$latex i = 0 , \ldots , {s-1}$$define 70 latex $71 a_i (x) 72 = 73 | \zeta_i ( x , a_0 (x) , \ldots , a_{i-1} (x ) ) | 74$$$
75 This defines $latex a : \B{R}^n \rightarrow \B{R}^s$$. 76 77 head g$$ 78 The object$icode g$$has prototype 79 codei% 80 ADFun<%Base%> %g% 81 %$$
82 The initial function representation in $icode g$$is lost. 83 Upon return it represents the smooth function 84 latex g : \B{R}^{n + s} \rightarrow \B{R}^{m + s}$$ is defined by 85$latex $86 g( x , u ) 87 = 88 \left[ \begin{array}{c} y(x, u) \\ z(x, u) \end{array} \right] 89$ $$90 were latex y(x, u)$$ and $latex z(x, u)$$are defined below. 91 92 subhead z(x, u)$$ 93 Define the smooth function 94$latex z : \B{R}^{n + s} \rightarrow \B{R}^s$$by 95 latex $96 z_i ( x , u ) = \zeta_i ( x , u_0 , \ldots , u_{i-1} ) 97$$$
98 Note that the partial of $latex z_i$$with respect to latex u_j$$ is zero 99 for$latex j \geq i$$. 100 101 subhead y(x, u)$$
102 There is a smooth function
103 $latex y : \B{R}^{n + s} \rightarrow \B{R}^m$$104 such that latex y( x , u ) = f(x)$$ whenever$latex u = a(x)$$. 105 106 head Affine Approximation$$
107 We define the affine approximations
108 $latex $109 \begin{array}{rcl} 110 y[ \hat{x} ]( x , u ) 111 & = & 112 y ( \hat{x}, a( \hat{x} ) ) 113 + \partial_x y ( \hat{x}, a( \hat{x} ) ) ( x - \hat{x} ) 114 + \partial_u y ( \hat{x}, a( \hat{x} ) ) ( u - a( \hat{x} ) ) 115 \\ 116 z[ \hat{x} ]( x , u ) 117 & = & 118 z ( \hat{x}, a( \hat{x} ) ) 119 + \partial_x z ( \hat{x}, a( \hat{x} ) ) ( x - \hat{x} ) 120 + \partial_u z ( \hat{x}, a( \hat{x} ) ) ( u - a( \hat{x} ) ) 121 \end{array} 122$ $$123 It follows that 124 latex $125 \begin{array}{rcl} 126 y( x , u ) 127 & = & 128 y[ \hat{x} ]( x , u ) + o ( x - \hat{x}, u - a( \hat{x} ) ) 129 \\ 130 z( x , u ) 131 & = & 132 z[ \hat{x} ]( x , u ) + o ( x - \hat{x}, u - a( \hat{x} ) ) 133 \end{array} 134$$$ 135 136$head Abs-normal Approximation$$137 138 subhead Approximating a(x)$$
139 The function $latex a(x)$$is not smooth, but it is equal to 140 latex | z(x, u) |$$ when$latex u = a(x)$$. 141 Furthermore 142 latex $143 z[ \hat{x} ]( x , u ) 144 = 145 z ( \hat{x}, a( \hat{x} ) ) 146 + \partial_x z ( \hat{x}, a( \hat{x} ) ) ( x - \hat{x} ) 147 + \partial_u z ( \hat{x}, a( \hat{x} ) ) ( u - a( \hat{x} ) ) 148$$$
149 The partial of $latex z_i$$with respect to latex u_j$$ is zero 150 for$latex j \geq i$$. It follows that 151 latex $152 z_i [ \hat{x} ]( x , u ) 153 = 154 z_i ( \hat{x}, a( \hat{x} ) ) 155 + \partial_x z_i ( \hat{x}, a( \hat{x} ) ) ( x - \hat{x} ) 156 + \sum_{j < i} \partial_{u(j)} 157 z_i ( \hat{x}, a( \hat{x} ) ) ( u_j - a_j ( \hat{x} ) ) 158$$$
159 Considering the case $latex i = 0$$we define 160 latex $161 a_0 [ \hat{x} ]( x ) 162 = 163 | z_0 [ \hat{x} ]( x , u ) | 164 = 165 \left| 166 z_0 ( \hat{x}, a( \hat{x} ) ) 167 + \partial_x z_0 ( \hat{x}, a( \hat{x} ) ) ( x - \hat{x} ) 168 \right| 169$$$ 170 It follows that 171$latex $172 a_0 (x) = a_0 [ \hat{x} ]( x ) + o ( x - \hat{x} ) 173$ $$174 In general, we define latex a_i [ \hat{x} ]$$ using
175 $latex a_j [ \hat{x} ]$$for latex j < i$$ as follows: 176$latex $177 a_i [ \hat{x} ]( x ) 178 = 179 \left | 180 z_i ( \hat{x}, a( \hat{x} ) ) 181 + \partial_x z_i ( \hat{x}, a( \hat{x} ) ) ( x - \hat{x} ) 182 + \sum_{j < i} \partial_{u(j)} 183 z_i ( \hat{x}, a( \hat{x} ) ) 184 ( a_j [ \hat{x} ] ( x ) - a_j ( \hat{x} ) ) 185 \right| 186$ $$187 It follows that 188 latex $189 a (x) = a[ \hat{x} ]( x ) + o ( x - \hat{x} ) 190$$$
191 Note that in the case where $latex z(x, u)$$and latex y(x, u)$$ are 192 affine, 193$latex $194 a[ \hat{x} ]( x ) = a( x ) 195$ $$196 197 198 subhead Approximating f(x)$$
199 $latex $200 f(x) 201 = 202 y ( x , a(x ) ) 203 = 204 y [ \hat{x} ] ( x , a[ \hat{x} ] ( x ) ) 205 + o( \Delta x ) 206$ $$207 208 head Correspondence to Literature$$ 209 Using the notation 210$latex Z = \partial_x z(\hat{x}, \hat{u})$$, 211 latex L = \partial_u z(\hat{x}, \hat{u})$$,
212 $latex J = \partial_x y(\hat{x}, \hat{u})$$, 213 latex Y = \partial_u y(\hat{x}, \hat{u})$$, 214 the approximation for$latex z$$and latex y$$ are
215 $latex $216 \begin{array}{rcl} 217 z[ \hat{x} ]( x , u ) 218 & = & 219 z ( \hat{x}, a( \hat{x} ) ) + Z ( x - \hat{x} ) + L ( u - a( \hat{x} ) ) 220 \\ 221 y[ \hat{x} ]( x , u ) 222 & = & 223 y ( \hat{x}, a( \hat{x} ) ) + J ( x - \hat{x} ) + Y ( u - a( \hat{x} ) ) 224 \end{array} 225$ $$226 Moving the terms with latex \hat{x}$$ together, we have 227$latex $228 \begin{array}{rcl} 229 z[ \hat{x} ]( x , u ) 230 & = & 231 z ( \hat{x}, a( \hat{x} ) ) - Z \hat{x} - L a( \hat{x} ) + Z x + L u 232 \\ 233 y[ \hat{x} ]( x , u ) 234 & = & 235 y ( \hat{x}, a( \hat{x} ) ) - J \hat{x} - Y a( \hat{x} ) + J x + Y u 236 \end{array} 237$ $$238 Using the notation 239 latex c = z ( \hat{x}, \hat{u} ) - Z \hat{x} - L \hat{u}$$,
240 $latex b = y ( \hat{x}, \hat{u} ) - J \hat{x} - Y \hat{u}$$, 241 we have 242 latex $243 \begin{array}{rcl} 244 z[ \hat{x} ]( x , u ) & = & c + Z x + L u 245 \\ 246 y[ \hat{x} ]( x , u ) & = & b + J x + Y u 247 \end{array} 248$$$ 249 Considering the affine case, where the approximations are exact, 250 and choosing$latex u = a(x) = |z(x, u)|$$, we obtain 251 latex $252 \begin{array}{rcl} 253 z( x , a(x ) ) & = & c + Z x + L |z( x , a(x ) )| 254 \\ 255 y( x , a(x ) ) & = & b + J x + Y |z( x , a(x ) )| 256 \end{array} 257$$$
258 This is Equation (2) of the
259 $cref/reference/abs_normal/Reference/$$. 260 261 children%example/abs_normal/get_started.cpp 262 %$$ 263$head Example$$264 The file cref abs_get_started.cpp$$ contains
265 an example and test using this operation.
266
267 \$end
268 -------------------------------------------------------------------------------
269 */
270 /*!
271 file abs_normal_fun.hpp
272 Create an abs-normal representation of a function
273 */
274
276 /*!
277 Create an abs-normal representation of an ADFun object.
278
279 \tparam Base
280 base type for this abs-normal form and for the function beging represented;
281 i.e., f.
282
283 \param f
284 is the function that this object will represent in abs-normal form.
285 This is effectively const except that the play back state play_
286 is used.
287 */
288
289 # define NOT_YET_COMPILING 0
290
291 template <class Base>
293 { using namespace local;
294
295  // -----------------------------------------------------------------------
296  // Forward sweep to determine number of absolute value operations in f
297  // -----------------------------------------------------------------------
298  // The argument and result index in f for each absolute value operator
301  //
302  OpCode op; // this operator
304  size_t i_op; // index of this operator
305  size_t i_var; // variable index for this operator
306  i_op = 0;
307  play_.get_op_info(i_op, op, arg, i_var);
308  CPPAD_ASSERT_UNKNOWN( op == BeginOp );
309  //
310  bool more_operators = true;
311  while( more_operators )
312  {
313  // next op
314  play_.get_op_info(++i_op, op, arg, i_var);
315  switch( op )
316  { // absolute value operator
317  case AbsOp:
319  f_abs_arg.push_back( arg[0] );
320  f_abs_res.push_back( i_var );
321  break;
322
323  case CSumOp:
324  break;
325
326  case CSkipOp:
327  break;
328
329  case EndOp:
330  more_operators = false;
331  break;
332
333  default:
334  break;
335  }
336  }
337  // ------------------------------------------------------------------------
338  // Forward sweep to create new recording
339  // ------------------------------------------------------------------------
340  // recorder for new operation sequence
341  recorder<Base> rec;
342  //
343  // number of variables in both operation sequences
344  // (the AbsOp operators are replace by InvOp operators)
345  const size_t num_var = play_.num_var_rec();
346  //
347  // mapping from old variable index to new variable index
350  );
352  for(i_var = 0; i_var < num_var; i_var++)
353  f2g_var[i_var] = addr_t( num_var ); // invalid (should not be used)
354  //
355  // record the independent variables in f
356  i_op = 0;
357  play_.get_op_info(i_op, op, arg, i_var);
358  CPPAD_ASSERT_UNKNOWN( op == BeginOp );
359  more_operators = true;
360  while( more_operators )
361  { switch( op )
362  {
363  // phantom variable
364  case BeginOp:
366  CPPAD_ASSERT_UNKNOWN( arg[0] == 0 );
367  rec.PutArg(0);
368  f2g_var[i_var] = rec.PutOp(op);
369  break;
370
371  // independent variables
372  case InvOp:
374  f2g_var[i_var] = rec.PutOp(op);
375  break;
376
377  // end of independent variables
378  default:
379  more_operators = false;
380  break;
381  }
382  if( more_operators )
383  play_.get_op_info(++i_op, op, arg, i_var);
384  }
385  // add one for the phantom variable
386  CPPAD_ASSERT_UNKNOWN( 1 + Domain() == i_var );
387  //
388  // record the independent variables corresponding AbsOp results
389  size_t index_abs;
390  for(index_abs = 0; index_abs < f_abs_res.size(); index_abs++)
391  f2g_var[ f_abs_res[index_abs] ] = rec.PutOp(InvOp);
392  //
393  // used to hold new argument vector
395  //
396  // Parameters in recording of f
397  const Base* f_parameter = play_.GetPar();
398  //
399  // now loop through the rest of the
400  more_operators = true;
401  index_abs = 0;
402  while( more_operators )
404  switch( op )
405  {
406  // check setting of f_abs_arg and f_abs_res;
407  case AbsOp:
409  CPPAD_ASSERT_UNKNOWN( f_abs_arg[index_abs] == arg[0] );
410  CPPAD_ASSERT_UNKNOWN( f_abs_res[index_abs] == i_var );
411  CPPAD_ASSERT_UNKNOWN( f2g_var[i_var] > 0 );
412  ++index_abs;
413  break;
414
415  // These operators come at beginning of take and are handled above
416  case InvOp:
418  break;
419
420  // ---------------------------------------------------------------
421  // Unary operators, argument a parameter, one result
422  case ParOp:
424  new_arg[0] = rec.PutPar( f_parameter[ arg[0] ] );
425  rec.PutArg( new_arg[0] );
426  f2g_var[i_var] = rec.PutOp(op);
427  break;
428
429  // --------------------------------------------------------------
430  // Unary operators, argument a variable, one result
431  // (excluding the absolute value operator AbsOp)
432  case AcosOp:
433  case AcoshOp:
434  case AsinOp:
435  case AsinhOp:
436  case AtanOp:
437  case AtanhOp:
438  case CosOp:
439  case CoshOp:
440  case ExpOp:
441  case Expm1Op:
442  case LogOp:
443  case Log1pOp:
444  case SignOp:
445  case SinOp:
446  case SinhOp:
447  case SqrtOp:
448  case TanOp:
449  case TanhOp:
450  // some of these operators have an auxillary result; e.g.,
451  // sine and cosine are computed togeather.
452  CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
453  CPPAD_ASSERT_UNKNOWN( NumRes(op) == 1 || NumRes(op) == 2 );
454  CPPAD_ASSERT_UNKNOWN( size_t( f2g_var[ arg[0] ] ) < num_var );
455  new_arg[0] = f2g_var[ arg[0] ];
456  rec.PutArg( new_arg[0] );
457  f2g_var[i_var] = rec.PutOp( op );
458  break;
459
460  case ErfOp:
462  CPPAD_ASSERT_UNKNOWN( size_t( f2g_var[ arg[0] ] ) < num_var );
463  // Error function is a special case
464  // second argument is always the parameter 0
465  // third argument is always the parameter 2 / sqrt(pi)
466  rec.PutArg( rec.PutPar( Base(0.0) ) );
467  rec.PutArg( rec.PutPar(
468  Base( 1.0 / std::sqrt( std::atan(1.0) ) )
469  ) );
470  f2g_var[i_var] = rec.PutOp(op);
471  break;
472  // --------------------------------------------------------------
473  // Binary operators, left variable, right parameter, one result
474  case SubvpOp:
475  case DivvpOp:
476  case PowvpOp:
477  case ZmulvpOp:
479  CPPAD_ASSERT_UNKNOWN( size_t( f2g_var[ arg[0] ] ) < num_var );
480  new_arg[0] = f2g_var[ arg[0] ];
481  new_arg[1] = rec.PutPar( f_parameter[ arg[1] ] );
482  rec.PutArg( new_arg[0], new_arg[1] );
483  f2g_var[i_var] = rec.PutOp(op);
484  break;
485  // ---------------------------------------------------
486  // Binary operators, left index, right variable, one result
487  case DisOp:
488  CPPAD_ASSERT_UNKNOWN( size_t( f2g_var[ arg[1] ] ) < num_var );
489  new_arg[0] = arg[0];
490  new_arg[1] = f2g_var[ arg[1] ];
491  rec.PutArg( new_arg[0], new_arg[1] );
492  f2g_var[i_var] = rec.PutOp(op);
493  break;
494
495  // --------------------------------------------------------------
496  // Binary operators, left parameter, right variable, one result
498  case SubpvOp:
499  case MulpvOp:
500  case DivpvOp:
501  case PowpvOp:
502  case ZmulpvOp:
504  CPPAD_ASSERT_UNKNOWN( size_t( f2g_var[ arg[1] ] ) < num_var );
505  new_arg[0] = rec.PutPar( f_parameter[ arg[0] ] );
506  new_arg[1] = f2g_var[ arg[1] ];
507  rec.PutArg( new_arg[0], new_arg[1] );
508  f2g_var[i_var] = rec.PutOp(op);
509  break;
510  // --------------------------------------------------------------
511  // Binary operators, left and right variables, one result
513  case SubvvOp:
514  case MulvvOp:
515  case DivvvOp:
516  case ZmulvvOp:
518  CPPAD_ASSERT_UNKNOWN( size_t( f2g_var[ arg[0] ] ) < num_var );
519  CPPAD_ASSERT_UNKNOWN( size_t( f2g_var[ arg[1] ] ) < num_var );
520  new_arg[0] = f2g_var[ arg[0] ];
521  new_arg[1] = f2g_var[ arg[1] ];
522  rec.PutArg( new_arg[0], new_arg[1] );
523  f2g_var[i_var] = rec.PutOp(op);
524  break;
525  // ---------------------------------------------------
526  // Conditional expression operators
527  case CExpOp:
529  new_arg[0] = arg[0];
530  new_arg[1] = arg[1];
532  for(size_t i = 2; i < 6; i++)
533  { if( arg[1] & mask )
534  { CPPAD_ASSERT_UNKNOWN( size_t(f2g_var[arg[i]]) < num_var );
535  new_arg[i] = f2g_var[ arg[i] ];
536  }
537  else
538  new_arg[i] = rec.PutPar( f_parameter[ arg[i] ] );
540  }
541  rec.PutArg(
542  new_arg[0] ,
543  new_arg[1] ,
544  new_arg[2] ,
545  new_arg[3] ,
546  new_arg[4] ,
547  new_arg[5]
548  );
549  f2g_var[i_var] = rec.PutOp(op);
550  break;
551
552  // --------------------------------------------------
554  case EndOp:
556  rec.PutOp(op);
557  more_operators = false;
558  break;
559
560  // ---------------------------------------------------
562  case LepvOp:
563  case LtpvOp:
564  case EqpvOp:
565  case NepvOp:
567  new_arg[0] = rec.PutPar( f_parameter[ arg[0] ] );
568  new_arg[1] = f2g_var[ arg[1] ];
569  rec.PutArg(new_arg[0], new_arg[1]);
570  rec.PutOp(op);
571  break;
572  //
573  case LevpOp:
574  case LtvpOp:
576  new_arg[0] = f2g_var[ arg[0] ];
577  new_arg[1] = rec.PutPar( f_parameter[ arg[1] ] );
578  rec.PutArg(new_arg[0], new_arg[1]);
579  rec.PutOp(op);
580  break;
581  //
582  case LevvOp:
583  case LtvvOp:
584  case EqvvOp:
585  case NevvOp:
587  new_arg[0] = f2g_var[ arg[0] ];
588  new_arg[1] = f2g_var[ arg[1] ];
589  rec.PutArg(new_arg[0], new_arg[1]);
590  rec.PutOp(op);
591  break;
592
593  // ---------------------------------------------------
594  // print forward operator
595  case PriOp:
597  //
598  // arg[0]
599  new_arg[0] = arg[0];
600  //
601  // arg[1]
602  if( arg[0] & 1 )
603  {
604  CPPAD_ASSERT_UNKNOWN( size_t( f2g_var[ arg[1] ] ) < num_var );
605  new_arg[1] = f2g_var[ arg[1] ];
606  }
607  else
608  { new_arg[1] = rec.PutPar( f_parameter[ arg[1] ] );
609  }
610  //
611  // arg[3]
612  if( arg[0] & 2 )
613  {
614  CPPAD_ASSERT_UNKNOWN( size_t( f2g_var[ arg[3] ] ) < num_var );
615  new_arg[3] = f2g_var[ arg[3] ];
616  }
617  else
618  { new_arg[3] = rec.PutPar( f_parameter[ arg[3] ] );
619  }
620  new_arg[2] = rec.PutTxt( play_.GetTxt( arg[2] ) );
621  new_arg[4] = rec.PutTxt( play_.GetTxt( arg[4] ) );
622  //
623  rec.PutArg(
624  new_arg[0] ,
625  new_arg[1] ,
626  new_arg[2] ,
627  new_arg[3] ,
628  new_arg[4]
629  );
630  // no result
631  rec.PutOp(op);
632  break;
633
634  // ---------------------------------------------------
636
637  // Load using a parameter index
638  case LdpOp:
640  new_arg[0] = arg[0];
641  new_arg[1] = arg[1];
642  new_arg[2] = arg[2];
643  rec.PutArg(
644  new_arg[0],
645  new_arg[1],
646  new_arg[2]
647  );
649  break;
650
651  // Load using a variable index
652  case LdvOp:
654  CPPAD_ASSERT_UNKNOWN( size_t( f2g_var[ arg[1] ] ) < num_var );
655  new_arg[0] = arg[0];
656  new_arg[1] = f2g_var[ arg[1] ];
657  new_arg[2] = arg[2];
658  rec.PutArg(
659  new_arg[0],
660  new_arg[1],
661  new_arg[2]
662  );
664  break;
665
666  // Store a parameter using a parameter index
667  case StppOp:
669  new_arg[0] = arg[0];
670  new_arg[1] = rec.PutPar( f_parameter[ arg[1] ] );
671  new_arg[2] = rec.PutPar( f_parameter[ arg[2] ] );
672  rec.PutArg(
673  new_arg[0],
674  new_arg[1],
675  new_arg[2]
676  );
677  rec.PutOp(op);
678  break;
679
680  // Store a parameter using a variable index
681  case StvpOp:
683  CPPAD_ASSERT_UNKNOWN( size_t( f2g_var[ arg[1] ] ) < num_var );
684  new_arg[0] = arg[0];
685  new_arg[1] = f2g_var[ arg[1] ];
686  new_arg[2] = rec.PutPar( f_parameter[ arg[2] ] );
687  rec.PutArg(
688  new_arg[0],
689  new_arg[1],
690  new_arg[2]
691  );
692  rec.PutOp(op);
693  break;
694
695  // Store a variable using a parameter index
696  case StpvOp:
698  CPPAD_ASSERT_UNKNOWN( size_t( f2g_var[ arg[2] ] ) < num_var );
699  new_arg[0] = arg[0];
700  new_arg[1] = rec.PutPar( f_parameter[ arg[1] ] );
701  new_arg[2] = f2g_var[ arg[2] ];
702  rec.PutArg(
703  new_arg[0],
704  new_arg[1],
705  new_arg[2]
706  );
707  rec.PutOp(op);
708  break;
709
710  // Store a variable using a variable index
711  case StvvOp:
713  CPPAD_ASSERT_UNKNOWN( size_t( f2g_var[ arg[1] ] ) < num_var );
714  CPPAD_ASSERT_UNKNOWN( size_t( f2g_var[ arg[2] ] ) < num_var );
715  new_arg[0] = arg[0];
716  new_arg[1] = f2g_var[ arg[1] ];
717  new_arg[2] = f2g_var[ arg[2] ];
718  rec.PutArg(
719  new_arg[0],
720  new_arg[1],
721  new_arg[2]
722  );
723  break;
724
725  // -----------------------------------------------------------
726  // user atomic function call operators
727
728  case UserOp:
730  // atomic_index, user_old, user_n, user_m
731  rec.PutArg(arg[0], arg[1], arg[2], arg[3]);
732  rec.PutOp(UserOp);
733  break;
734
735  case UsrapOp:
737  new_arg[0] = rec.PutPar( f_parameter[ arg[0] ] );
738  rec.PutArg(new_arg[0]);
739  rec.PutOp(UsrapOp);
740  break;
741
742  case UsravOp:
744  CPPAD_ASSERT_UNKNOWN( size_t( f2g_var[arg[0]] ) < num_var );
745  new_arg[0] = f2g_var[ arg[0] ];
746  rec.PutArg(new_arg[0]);
747  rec.PutOp(UsravOp);
748  break;
749
750  case UsrrpOp:
752  new_arg[0] = rec.PutPar( f_parameter[ arg[0] ] );
753  rec.PutArg(new_arg[0]);
754  rec.PutOp(UsrrpOp);
755  break;
756
757  case UsrrvOp:
759  f2g_var[i_var] = rec.PutOp(UsrrvOp);
760  break;
761  // ---------------------------------------------------
762
763  // all cases should be handled above
764  default:
766  }
767  if( more_operators )
768  play_.get_op_info(++i_op, op, arg, i_var);
769  }
770  // Check a few expected results
771  CPPAD_ASSERT_UNKNOWN( rec.num_op_rec() == play_.num_op_rec() );
772  CPPAD_ASSERT_UNKNOWN( rec.num_var_rec() == play_.num_var_rec() );
774
775  // -----------------------------------------------------------------------
776  // Use rec to create the function g
777  // -----------------------------------------------------------------------
778
779  // number of variables in the recording
780  g.num_var_tape_ = rec.num_var_rec();
781
782  // dimension cskip_op vector to number of operators
783  g.cskip_op_.resize( rec.num_op_rec() );
784
785  // independent variables in g: (x, u)
786  size_t s = f_abs_res.size();
787  size_t n = Domain();
789  // (x, u)
790  for(size_t j = 0; j < n; j++)
793  }
794  for(size_t j = 0; j < s; j++)
795  { g.ind_taddr_[n + j] = f2g_var[ f_abs_res[j] ];
796  CPPAD_ASSERT_UNKNOWN( g.ind_taddr_[n + j] == n + j + 1 );
797  }
798
799  // dependent variable in g: (y, z)
800  CPPAD_ASSERT_UNKNOWN( s == f_abs_arg.size() );
801  size_t m = Range();
803  for(size_t i = 0; i < m; i++)
806  }
807  for(size_t i = 0; i < s; i++)
808  { g.dep_taddr_[m + i] = f2g_var[ f_abs_arg[i] ];
810  }
811
812  // which dependent variables are parameters
813  g.dep_parameter_.resize(m + s);
814  for(size_t i = 0; i < m; i++)
815  g.dep_parameter_[i] = dep_parameter_[i];
816  for(size_t i = 0; i < s; i++)
817  g.dep_parameter_[m + i] = false;
818
819  // free memory allocated for sparse Jacobian calculation
820  // (the resutls are no longer valid)
821  g.for_jac_sparse_pack_.resize(0, 0);
822  g.for_jac_sparse_set_.resize(0, 0);
823
824  // free taylor coefficient memory
825  g.taylor_.clear();
826  g.num_order_taylor_ = 0;
827  g.cap_order_taylor_ = 0;
828
829  // Transferring the recording swaps its vectors so do this last
830  // replace the recording in g (this ADFun object)
831  g.play_.get(rec, n + s);
832
833  // resize subgraph_info_
834  g.subgraph_info_.resize(
837  g.play_.num_op_rec(), // n_op
838  g.play_.num_var_rec() // n_var
839  );
840
841  // ------------------------------------------------------------------------
842  // Create the function a
843  // ------------------------------------------------------------------------
844
846  a = *this;
847
848  // dependent variables in a(x)
849  CPPAD_ASSERT_UNKNOWN( s == f_abs_arg.size() );
851  for(size_t i = 0; i < s; i++)
854  }
855
856  // free memory allocated for sparse Jacobian calculation
857  // (the resutls are no longer valid)
858  a.for_jac_sparse_pack_.resize(0, 0);
859  a.for_jac_sparse_set_.resize(0, 0);
860
861  // free taylor coefficient memory
862  a.taylor_.clear();
863  a.num_order_taylor_ = 0;
864  a.cap_order_taylor_ = 0;
865 }
866
868
869 # endif
local::sparse_list for_jac_sparse_set_
Set results of the forward mode Jacobian sparsity calculations for_jac_sparse_set_.n_set() != 0 implies for_sparse_pack_ is empty.
Class used to hold function objects.
size_t num_order_taylor_
number of orders stored in taylor_
The CppAD Simple Vector template class.
Definition: vector.hpp:338
size_t cap_order_taylor_
maximum number of orders that will fit in taylor_
local::sparse_pack for_jac_sparse_pack_
Packed results of the forward mode Jacobian sparsity calculations. for_jac_sparse_pack_.n_set() != 0 implies other sparsity results are empty.
std::complex< double > atan(const std::complex< double > &x)
size_t NumArg(OpCode op)
Number of arguments for a specified operator.
Definition: op_code.hpp:175
tape address for the independent variables
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
local::player< Base > play_
the operation sequence corresponding to this object
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
size_t num_var_tape_
number of variables in the recording (play_)
local::pod_vector< bool > cskip_op_
which operations can be conditionally skipped Set during forward pass of order zero ...
Check that exp is true, if not terminate execution.