CppAD: A C++ Algorithmic Differentiation Package  20171217
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
abs_normal_fun.hpp
Go to the documentation of this file.
1 # ifndef CPPAD_CORE_ABS_NORMAL_FUN_HPP
2 # define CPPAD_CORE_ABS_NORMAL_FUN_HPP
3 /* --------------------------------------------------------------------------
4 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
5 
6 CppAD is distributed under multiple licenses. This distribution is under
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.
11 Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
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 
275 namespace CppAD { // BEGIN_CPPAD_NAMESPACE
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
299  CppAD::vector<addr_t> f_abs_arg;
300  CppAD::vector<size_t> f_abs_res;
301  //
302  OpCode op; // this operator
303  const addr_t* arg = CPPAD_NULL; // arguments for 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:
318  CPPAD_ASSERT_NARG_NRES(op, 1, 1);
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
349  std::numeric_limits<addr_t>::max() >= num_var
350  );
351  CppAD::vector<addr_t> f2g_var(num_var);
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:
365  CPPAD_ASSERT_NARG_NRES(op, 1, 1);
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:
373  CPPAD_ASSERT_NARG_NRES(op, 0, 1);
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
394  addr_t new_arg[6];
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 )
403  { addr_t mask; // temporary used in some switch cases
404  switch( op )
405  {
406  // check setting of f_abs_arg and f_abs_res;
407  case AbsOp:
408  CPPAD_ASSERT_NARG_NRES(op, 1, 1);
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:
417  CPPAD_ASSERT_UNKNOWN(false);
418  break;
419 
420  // ---------------------------------------------------------------
421  // Unary operators, argument a parameter, one result
422  case ParOp:
423  CPPAD_ASSERT_NARG_NRES(op, 1, 1);
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:
461  CPPAD_ASSERT_NARG_NRES(op, 3, 5);
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:
478  CPPAD_ASSERT_NARG_NRES(op, 2, 1);
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
497  case AddpvOp:
498  case SubpvOp:
499  case MulpvOp:
500  case DivpvOp:
501  case PowpvOp:
502  case ZmulpvOp:
503  CPPAD_ASSERT_NARG_NRES(op, 2, 1);
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
512  case AddvvOp:
513  case SubvvOp:
514  case MulvvOp:
515  case DivvvOp:
516  case ZmulvvOp:
517  CPPAD_ASSERT_NARG_NRES(op, 2, 1);
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:
528  CPPAD_ASSERT_NARG_NRES(op, 6, 1);
529  new_arg[0] = arg[0];
530  new_arg[1] = arg[1];
531  mask = 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] ] );
539  mask = mask << 1;
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  // --------------------------------------------------
553  // Operators with no arguments and no results
554  case EndOp:
555  CPPAD_ASSERT_NARG_NRES(op, 0, 0);
556  rec.PutOp(op);
557  more_operators = false;
558  break;
559 
560  // ---------------------------------------------------
561  // Operations with two arguments and no results
562  case LepvOp:
563  case LtpvOp:
564  case EqpvOp:
565  case NepvOp:
566  CPPAD_ASSERT_NARG_NRES(op, 2, 0);
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:
575  CPPAD_ASSERT_NARG_NRES(op, 2, 0);
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:
586  CPPAD_ASSERT_NARG_NRES(op, 2, 0);
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:
596  CPPAD_ASSERT_NARG_NRES(op, 5, 0);
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  // ---------------------------------------------------
635  // VecAD operators
636 
637  // Load using a parameter index
638  case LdpOp:
639  CPPAD_ASSERT_NARG_NRES(op, 3, 1);
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  );
648  f2g_var[i_var] = rec.PutLoadOp(op);
649  break;
650 
651  // Load using a variable index
652  case LdvOp:
653  CPPAD_ASSERT_NARG_NRES(op, 3, 1);
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  );
663  f2g_var[i_var] = rec.PutLoadOp(op);
664  break;
665 
666  // Store a parameter using a parameter index
667  case StppOp:
668  CPPAD_ASSERT_NARG_NRES(op, 3, 0);
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:
682  CPPAD_ASSERT_NARG_NRES(op, 3, 0);
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:
697  CPPAD_ASSERT_NARG_NRES(op, 3, 0);
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:
712  CPPAD_ASSERT_NARG_NRES(op, 3, 0);
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:
729  CPPAD_ASSERT_NARG_NRES(op, 4, 0);
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:
736  CPPAD_ASSERT_NARG_NRES(op, 1, 0);
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:
743  CPPAD_ASSERT_NARG_NRES(op, 1, 0);
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:
751  CPPAD_ASSERT_NARG_NRES(op, 1, 0);
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:
758  CPPAD_ASSERT_NARG_NRES(op, 0, 1);
759  f2g_var[i_var] = rec.PutOp(UsrrvOp);
760  break;
761  // ---------------------------------------------------
762 
763  // all cases should be handled above
764  default:
765  CPPAD_ASSERT_UNKNOWN(false);
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() );
773  CPPAD_ASSERT_UNKNOWN( rec.num_load_op_rec() == play_.num_load_op_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();
788  g.ind_taddr_.resize(n + s);
789  // (x, u)
790  for(size_t j = 0; j < n; j++)
791  { g.ind_taddr_[j] = f2g_var[ ind_taddr_[j] ];
792  CPPAD_ASSERT_UNKNOWN( g.ind_taddr_[j] == j + 1 );
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();
802  g.dep_taddr_.resize(m + s);
803  for(size_t i = 0; i < m; i++)
804  { g.dep_taddr_[i] = f2g_var[ dep_taddr_[i] ];
805  CPPAD_ASSERT_UNKNOWN( g.dep_taddr_[i] < num_var );
806  }
807  for(size_t i = 0; i < s; i++)
808  { g.dep_taddr_[m + i] = f2g_var[ f_abs_arg[i] ];
809  CPPAD_ASSERT_UNKNOWN( g.dep_taddr_[m + i] < num_var );
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(
835  g.ind_taddr_.size(), // n_ind
836  g.dep_taddr_.size(), // n_dep
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 
845  // start with a copy of f
846  a = *this;
847 
848  // dependent variables in a(x)
849  CPPAD_ASSERT_UNKNOWN( s == f_abs_arg.size() );
850  a.dep_taddr_.resize(s);
851  for(size_t i = 0; i < s; i++)
852  { a.dep_taddr_[i] = f_abs_res[i];
853  CPPAD_ASSERT_UNKNOWN( a.dep_taddr_[i] < num_var );
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 
867 } // END_CPPAD_NAMESPACE
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.
Definition: ad_fun.hpp:140
Class used to hold function objects.
Definition: ad_fun.hpp:69
size_t num_order_taylor_
number of orders stored in taylor_
Definition: ad_fun.hpp:94
The CppAD Simple Vector template class.
Definition: vector.hpp:338
size_t cap_order_taylor_
maximum number of orders that will fit in taylor_
Definition: ad_fun.hpp:97
CPPAD_TAPE_ADDR_TYPE addr_t
Definition: declare_ad.hpp:44
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.
Definition: ad_fun.hpp:136
void abs_normal_fun(ADFun &g, ADFun &a) const
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
CppAD::vector< size_t > ind_taddr_
tape address for the independent variables
Definition: ad_fun.hpp:106
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
Definition: ad_fun.hpp:131
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_)
Definition: ad_fun.hpp:103
AD< Base > sqrt(const AD< Base > &x)
local::pod_vector< bool > cskip_op_
which operations can be conditionally skipped Set during forward pass of order zero ...
Definition: ad_fun.hpp:124
#define CPPAD_ASSERT_UNKNOWN(exp)
Check that exp is true, if not terminate execution.
#define CPPAD_ASSERT_NARG_NRES(op, n_arg, n_res)
Check that operator op has the specified number of of arguments and results.
CppAD::vector< bool > dep_parameter_
which dependent variables are actually parameters
Definition: ad_fun.hpp:112
CppAD::vector< size_t > dep_taddr_
tape address and parameter flag for the dependent variables
Definition: ad_fun.hpp:109
local::pod_vector< Base > taylor_
results of the forward mode calculations
Definition: ad_fun.hpp:115
local::subgraph::subgraph_info subgraph_info_
subgraph information for this object
Definition: ad_fun.hpp:143