CppAD: A C++ Algorithmic Differentiation Package  20130918
optimize.hpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 # ifndef CPPAD_OPTIMIZE_INCLUDED
00003 # define CPPAD_OPTIMIZE_INCLUDED
00004 
00005 /* --------------------------------------------------------------------------
00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-14 Bradley M. Bell
00007 
00008 CppAD is distributed under multiple licenses. This distribution is under
00009 the terms of the 
00010                     Eclipse Public License Version 1.0.
00011 
00012 A copy of this license is included in the COPYING file of this distribution.
00013 Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
00014 -------------------------------------------------------------------------- */
00015 
00016 /*
00017 $begin optimize$$
00018 $spell
00019      jac
00020      bool
00021      Taylor
00022      var
00023      CppAD
00024      cppad
00025      std
00026      CondExpEq
00027 $$
00028 
00029 $section Optimize an ADFun Object Tape$$
00030 
00031 $index optimize$$
00032 $index tape, optimize$$
00033 $index sequence, optimize operations$$
00034 $index operations, optimize sequence$$
00035 $index speed, optimize$$
00036 $index memory, optimize$$
00037 
00038 $head Syntax$$
00039 $icode%f%.optimize()%$$
00040 
00041 
00042 $head Purpose$$
00043 The operation sequence corresponding to an $cref ADFun$$ object can
00044 be very large and involve many operations; see the
00045 size functions in $cref seq_property$$.
00046 The $icode%f%.optimize%$$ procedure reduces the number of operations,
00047 and thereby the time and the memory, required to
00048 compute function and derivative values. 
00049 
00050 $head f$$
00051 The object $icode f$$ has prototype
00052 $codei%
00053      ADFun<%Base%> %f%
00054 %$$
00055 
00056 $head Improvements$$
00057 You can see the reduction in number of variables in the operation sequence
00058 by calling the function $cref/f.size_var()/seq_property/size_var/$$
00059 before and after the optimization procedure.
00060 Given that the optimization procedure takes time,
00061 it may be faster to skip this optimize procedure and just compute
00062 derivatives using the original operation sequence.
00063 
00064 $subhead Testing$$
00065 You can run the CppAD $cref/speed/speed_main/$$ tests and see
00066 the corresponding changes in number of variables and execution time; 
00067 see $cref cmake_check$$.
00068 
00069 $head Efficiency$$
00070 The $code optimize$$ member function
00071 may greatly reduce the number of variables 
00072 in the operation sequence; see $cref/size_var/seq_property/size_var/$$.
00073 If a $cref/zero order forward/forward_zero/$$ calculation is done during
00074 the construction of $icode f$$, it will require more memory
00075 and time than required after the optimization procedure.
00076 In addition, it will need to be redone.
00077 For this reason, it is more efficient to use 
00078 $codei%
00079      ADFun<%Base%> %f%;
00080      %f%.Dependent(%x%, %y%);
00081      %f%.optimize();
00082 %$$
00083 instead of
00084 $codei%
00085      ADFun<%Base%> %f%(%x%, %y%)
00086      %f%.optimize();
00087 %$$ 
00088 See the discussion about
00089 $cref/sequence constructors/FunConstruct/Sequence Constructor/$$.
00090 
00091 $head Comparison Operators$$
00092 Any comparison operators that are in the tape are removed by this operation.
00093 Hence the return value of $cref CompareChange$$ will always be zero
00094 for an optimized tape (even if $code NDEBUG$$ is not defined).
00095 
00096 $head Atomic Functions$$
00097 There are some subtitle issue with optimized $cref atomic$$ functions
00098 $latex v = g(u)$$:
00099 
00100 $subhead rev_sparse_jac$$
00101 The $cref atomic_rev_sparse_jac$$ function is be used to determine
00102 which components of $icode u$$ affect the dependent variables of $icode f$$.
00103 The current setting of the 
00104 $cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for each
00105 atomic function is used to determine if the $code bool$$ or 
00106 $code std::set<size_t>$$ version of $cref atomic_rev_sparse_jac$$ is used.
00107 
00108 $subhead nan$$
00109 If $icode%u%[%i%]%$$ does not affect the value of
00110 the dependent variables for $icode f$$,
00111 the value of $icode%u%[%i%]%$$ is set to $cref nan$$.
00112 
00113 
00114 $head Checking Optimization$$
00115 $index NDEBUG$$
00116 If $cref/NDEBUG/Faq/Speed/NDEBUG/$$ is not defined,
00117 and $cref/f.size_order()/size_order/$$ is greater than zero,
00118 a $cref forward_zero$$ calculation is done using the optimized version
00119 of $icode f$$ and the results are checked to see that they are
00120 the same as before.
00121 If they are not the same, the
00122 $cref ErrorHandler$$ is called with a known error message
00123 related to $icode%f%.optimize()%$$.
00124 
00125 $head Example$$
00126 $children%
00127      example/optimize.cpp
00128 %$$
00129 The file
00130 $cref optimize.cpp$$
00131 contains an example and test of this operation.
00132 It returns true if it succeeds and false otherwise.
00133 
00134 $end
00135 -----------------------------------------------------------------------------
00136 */
00137 # include <stack>
00138 
00139 namespace CppAD { // BEGIN_CPPAD_NAMESPACE
00140 namespace optimize { // BEGIN_CPPAD_OPTIMIZE_NAMESPACE
00141 /*!
00142 \file optimize.hpp
00143 Routines for optimizing a tape
00144 */
00145 
00146 
00147 /*!
00148 State for this variable set during reverse sweep.
00149 */
00150 enum enum_connect_type {
00151      /// There is no operation that connects this variable to the
00152      /// independent variables.
00153      not_connected        ,
00154 
00155      /// There is one or more operations that connects this variable to the
00156      /// independent variables.
00157      yes_connected        ,
00158 
00159      /// There is only one parrent that connects this variable to the 
00160      /// independent variables and the parent is a summation operation; i.e.,
00161      /// AddvvOp, AddpvOp, SubpvOp, SubvpOp, or SubvvOp.
00162      sum_connected        ,
00163 
00164      /// Satisfies the sum_connected assumptions above and in addition 
00165      /// this variable is the result of summation operator.
00166      csum_connected       ,
00167 
00168      /// This node is only connected in the case where the comparision is 
00169      /// true for the conditional expression with index \c connect_index.
00170      cexp_connected
00171 
00172 };
00173 
00174 /*!
00175 Class used to hold information about one conditional expression.
00176 */
00177 class class_cexp_pair {
00178 public:
00179      /// If this is true (false) this connection is only for the case where
00180      /// the comparision in the conditional expression is true (false)
00181      bool compare_; 
00182      /// This is the index of the conditional expression (in cksip_info)
00183      /// for this connection
00184      size_t index_;
00185      /// constructor
00186      class_cexp_pair(const bool& compare, const size_t& index)
00187      : compare_(compare), index_(index)
00188      { }
00189      /// assignment operator
00190      void operator=(const class_cexp_pair& right)
00191      {    index_   = right.index_;
00192           compare_ = right.compare_;
00193      }
00194      /// not equal operator
00195      bool operator!=(const class_cexp_pair& right)
00196      {    return (index_ != right.index_) | (compare_ != right.compare_);
00197      }
00198      /// Less than operator 
00199      /// (required for intersection of two sets of class_cexp_pair elements).
00200      bool operator<(const class_cexp_pair& right) const
00201      {    if( index_ == right.index_ )
00202                return size_t(compare_) < size_t(right.compare_); 
00203           return index_ < right.index_;
00204      }
00205 };
00206 /*!
00207 Compute intersection of two sets of class_cexp_pair elements.
00208 
00209 \param left
00210 first operand of the intersection
00211 
00212 \param right 
00213 second operand of the intersection
00214 
00215 \result
00216 the intersection of left and right
00217 */
00218 inline std::set<class_cexp_pair> intersection(
00219      std::set<class_cexp_pair>& left  ,
00220      std::set<class_cexp_pair>& right )
00221 {    std::set<class_cexp_pair> result;
00222      std::set_intersection(
00223           left.begin()  ,
00224           left.end()    ,
00225           right.begin() ,
00226           right.end()   ,
00227           std::inserter(result, result.begin())
00228      );
00229      return result; 
00230 }
00231 
00232 /*!
00233 Structure used by \c optimize to hold information about one variable.
00234 in the old operation seqeunce.
00235 */
00236 struct struct_old_variable {
00237      /// Operator for which this variable is the result, \c NumRes(op) > 0.
00238      /// Set by the reverse sweep at beginning of optimization.
00239      OpCode              op;       
00240 
00241      /// Pointer to first argument (child) for this operator.
00242      /// Set by the reverse sweep at beginning of optimization.
00243      const addr_t*       arg;
00244 
00245      /// How is this variable connected to the independent variables
00246      enum_connect_type connect_type; 
00247 
00248      /*!
00249      If \c connect_type is \c cexp_connected,
00250      this is the corresponding infromation for the conditional connections.
00251      */
00252      std::set<class_cexp_pair> cexp_set;
00253 
00254      /// New operation sequence corresponding to this old varable.
00255      /// Set during forward sweep to the index in the new tape
00256      addr_t new_var;
00257 
00258      /// New operator index for this varable.
00259      /// Set during forward sweep to the index in the new tape
00260      size_t new_op;
00261 };
00262 
00263 struct struct_size_pair {
00264      size_t i_op;  // an operator index
00265      size_t i_var; // a variable index
00266 };
00267 
00268 /*!
00269 Structures used by \c record_csum
00270 to hold information about one variable.
00271 */
00272 struct struct_csum_variable {
00273      /// Operator for which this variable is the result, \c NumRes(op) > 0.
00274      OpCode              op;       
00275 
00276      /// Pointer to first argument (child) for this operator.
00277      /// Set by the reverse sweep at beginning of optimization.
00278      const addr_t*       arg;
00279 
00280      /// Is this variable added to the summation
00281      /// (if not it is subtracted)
00282      bool                add;
00283 };
00284 
00285 /*!
00286 Structure used to pass work space from \c optimize to \c record_csum
00287 (so that stacks do not start from zero size every time).
00288 */
00289 struct struct_csum_stacks {
00290      /// stack of operations in the cummulative summation 
00291      std::stack<struct struct_csum_variable>   op_stack;
00292      /// stack of variables to be added
00293      std::stack<size_t >                         add_stack;
00294      /// stack of variables to be subtracted
00295      std::stack<size_t >                         sub_stack;
00296 };
00297 
00298 /*!
00299 CExpOp information that is copied to corresponding CSkipOp
00300 */
00301 struct struct_cskip_info {
00302      /// comparision operator
00303      CompareOp cop;
00304      /// (flag & 1) is true if and only if left is a variable
00305      /// (flag & 2) is true if and only if right is a variable
00306      size_t flag;
00307      /// index for left comparison operand
00308      size_t left; 
00309      /// index for right comparison operand
00310      size_t right; 
00311      /// maximum variable index between left and right
00312      size_t max_left_right;
00313      /// set of variables to skip on true
00314      CppAD::vector<size_t> skip_var_true;
00315      /// set of variables to skip on false
00316      CppAD::vector<size_t> skip_var_false;
00317      /// set of operations to skip on true
00318      CppAD::vector<size_t> skip_op_true;
00319      /// set of operations to skip on false
00320      CppAD::vector<size_t> skip_op_false;
00321      /// size of skip_op_true
00322      size_t n_op_true;
00323      /// size of skip_op_false
00324      size_t n_op_false;
00325      /// index in the argument recording of first argument for this CSkipOp
00326      size_t i_arg;
00327 };
00328 /*!
00329 Connection information for a user atomic function
00330 */
00331 struct struct_user_info {
00332      /// type of connection for this atomic function
00333      enum_connect_type connect_type;
00334      /// If this is an conditional connection, this is the information
00335      /// of the correpsonding CondExpOp operators
00336      std::set<class_cexp_pair> cexp_set;
00337      /// If this is a conditional connection, this is the operator
00338      /// index of the beginning of the atomic call sequence; i.e.,
00339      /// the first UserOp.
00340      size_t op_begin;
00341      /// If this is a conditional connection, this is one more than the
00342      ///  operator index of the ending of the atomic call sequence; i.e.,
00343      /// the second UserOp.
00344      size_t op_end;
00345 };
00346 
00347 /*!
00348 Shared documentation for optimization helper functions (not called).
00349 
00350 <!-- define prototype -->
00351 \param tape
00352 is a vector that maps a variable index, in the old operation sequence,
00353 to an <tt>struct_old_variable</tt> information record.
00354 Note that the index for this vector must be greater than or equal zero and 
00355 less than <tt>tape.size()</tt>.
00356 
00357 \li <tt>tape[i].op</tt> 
00358 is the operator in the old operation sequence
00359 corresponding to the old variable index \c i.
00360 Assertion: <tt>NumRes(tape[i].op) > 0</tt>.
00361 
00362 \li <tt>tape[i].arg</tt> 
00363 for <tt>j < NumArg( tape[i].op ), tape[i].arg[j]</tt>
00364 is the j-th the argument, in the old operation sequence,
00365 corresponding to the old variable index \c i.
00366 Assertion: <tt>tape[i].arg[j] < i</tt>.
00367 
00368 \li <tt>tape[i].new_var</tt>
00369 Suppose 
00370 <tt>i <= current, j < NumArg( tape[i].op ), and k = tape[i].arg[j]</tt>,
00371 and \c j corresponds to a variable for operator <tt>tape[i].op</tt>.
00372 It follows that <tt>tape[k].new_var</tt>
00373 has alread been set to the variable in the new operation sequence 
00374 corresponding to the old variable index \c k.
00375 This means that the \c new_var value has been set
00376 for all the possible arguments that come before \a current.
00377 
00378 \param current
00379 is the index in the old operation sequence for 
00380 the variable corresponding to the result for the current operator.
00381 Assertions: 
00382 <tt>
00383 current < tape.size(),
00384 NumRes( tape[current].op ) > 0.
00385 </tt>
00386 
00387 \param npar
00388 is the number of parameters corresponding to this operation sequence.
00389 
00390 \param par
00391 is a vector of length \a npar containing the parameters
00392 for this operation sequence; i.e.,
00393 given a parameter index \c i, the corresponding parameter value is
00394 <tt>par[i]</tt>.
00395 <!-- end prototype -->
00396 */
00397 template <class Base>
00398 void prototype(
00399      const CppAD::vector<struct struct_old_variable>& tape           ,
00400      size_t                                             current        ,
00401      size_t                                             npar           ,
00402      const Base*                                        par            )
00403 {    CPPAD_ASSERT_UNKNOWN(false); }
00404 
00405 /*!
00406 Check a unary operator for a complete match with a previous operator.
00407 
00408 A complete match means that the result of the previous operator
00409 can be used inplace of the result for current operator.
00410 
00411 <!-- replace prototype -->
00412 \param tape
00413 is a vector that maps a variable index, in the old operation sequence,
00414 to an <tt>struct_old_variable</tt> information record.
00415 Note that the index for this vector must be greater than or equal zero and 
00416 less than <tt>tape.size()</tt>.
00417 
00418 \li <tt>tape[i].op</tt> 
00419 is the operator in the old operation sequence
00420 corresponding to the old variable index \c i.
00421 Assertion: <tt>NumRes(tape[i].op) > 0</tt>.
00422 
00423 \li <tt>tape[i].arg</tt> 
00424 for <tt>j < NumArg( tape[i].op ), tape[i].arg[j]</tt>
00425 is the j-th the argument, in the old operation sequence,
00426 corresponding to the old variable index \c i.
00427 Assertion: <tt>tape[i].arg[j] < i</tt>.
00428 
00429 \li <tt>tape[i].new_var</tt>
00430 Suppose 
00431 <tt>i <= current, j < NumArg( tape[i].op ), and k = tape[i].arg[j]</tt>,
00432 and \c j corresponds to a variable for operator <tt>tape[i].op</tt>.
00433 It follows that <tt>tape[k].new_var</tt>
00434 has alread been set to the variable in the new operation sequence 
00435 corresponding to the old variable index \c k.
00436 This means that the \c new_var value has been set
00437 for all the possible arguments that come before \a current.
00438 
00439 \param current
00440 is the index in the old operation sequence for 
00441 the variable corresponding to the result for the current operator.
00442 Assertions: 
00443 <tt>
00444 current < tape.size(),
00445 NumRes( tape[current].op ) > 0.
00446 </tt>
00447 
00448 \param npar
00449 is the number of parameters corresponding to this operation sequence.
00450 
00451 \param par
00452 is a vector of length \a npar containing the parameters
00453 for this operation sequence; i.e.,
00454 given a parameter index \c i, the corresponding parameter value is
00455 <tt>par[i]</tt>.
00456 <!-- end prototype -->
00457 
00458 \param hash_table_var
00459 is a vector with size CPPAD_HASH_TABLE_SIZE
00460 that maps a hash code to the corresponding 
00461 variable index in the old operation sequence.
00462 All the values in this table must be less than \a current.
00463 
00464 \param code
00465 The input value of code does not matter.
00466 The output value of code is the hash code corresponding to
00467 this operation in the new operation sequence.
00468 
00469 \return
00470 If the return value is zero,
00471 no match was found.
00472 If the return value is greater than zero,
00473 it is the index of a new variable that can be used to replace the 
00474 old variable.
00475 
00476 \par Restrictions:
00477 NumArg( tape[current].op ) == 1
00478 */
00479 template <class Base>
00480 size_t unary_match(
00481      const CppAD::vector<struct struct_old_variable>& tape           ,
00482      size_t                                             current        ,
00483      size_t                                             npar           ,
00484      const Base*                                        par            ,
00485      const CppAD::vector<size_t>&                       hash_table_var ,
00486      unsigned short&                                    code           )
00487 {    const addr_t* arg = tape[current].arg;
00488      OpCode        op  = tape[current].op;
00489      addr_t new_arg[1];
00490      
00491      CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
00492      CPPAD_ASSERT_UNKNOWN( NumRes(op) > 0  );
00493      CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < current );
00494      new_arg[0] = tape[arg[0]].new_var;
00495      CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < current );
00496      code = hash_code(
00497           op                  , 
00498           new_arg             ,
00499           npar                ,
00500           par
00501      );
00502      size_t  i               = hash_table_var[code];
00503      CPPAD_ASSERT_UNKNOWN( i < current );
00504      if( op == tape[i].op )
00505      {    size_t k = tape[i].arg[0];
00506           CPPAD_ASSERT_UNKNOWN( k < i );
00507           if (new_arg[0] == tape[k].new_var )
00508                return tape[i].new_var;
00509      }
00510      return 0;
00511 } 
00512 
00513 /*!
00514 Check a binary operator for a complete match with a previous operator,
00515 
00516 <!-- replace prototype -->
00517 \param tape
00518 is a vector that maps a variable index, in the old operation sequence,
00519 to an <tt>struct_old_variable</tt> information record.
00520 Note that the index for this vector must be greater than or equal zero and 
00521 less than <tt>tape.size()</tt>.
00522 
00523 \li <tt>tape[i].op</tt> 
00524 is the operator in the old operation sequence
00525 corresponding to the old variable index \c i.
00526 Assertion: <tt>NumRes(tape[i].op) > 0</tt>.
00527 
00528 \li <tt>tape[i].arg</tt> 
00529 for <tt>j < NumArg( tape[i].op ), tape[i].arg[j]</tt>
00530 is the j-th the argument, in the old operation sequence,
00531 corresponding to the old variable index \c i.
00532 Assertion: <tt>tape[i].arg[j] < i</tt>.
00533 
00534 \li <tt>tape[i].new_var</tt>
00535 Suppose 
00536 <tt>i <= current, j < NumArg( tape[i].op ), and k = tape[i].arg[j]</tt>,
00537 and \c j corresponds to a variable for operator <tt>tape[i].op</tt>.
00538 It follows that <tt>tape[k].new_var</tt>
00539 has alread been set to the variable in the new operation sequence 
00540 corresponding to the old variable index \c k.
00541 This means that the \c new_var value has been set
00542 for all the possible arguments that come before \a current.
00543 
00544 \param current
00545 is the index in the old operation sequence for 
00546 the variable corresponding to the result for the current operator.
00547 Assertions: 
00548 <tt>
00549 current < tape.size(),
00550 NumRes( tape[current].op ) > 0.
00551 </tt>
00552 
00553 \param npar
00554 is the number of parameters corresponding to this operation sequence.
00555 
00556 \param par
00557 is a vector of length \a npar containing the parameters
00558 for this operation sequence; i.e.,
00559 given a parameter index \c i, the corresponding parameter value is
00560 <tt>par[i]</tt>.
00561 <!-- end prototype -->
00562 
00563 \param hash_table_var
00564 is a vector with size CPPAD_HASH_TABLE_SIZE
00565 that maps a hash code to the corresponding 
00566 variable index in the old operation sequence.
00567 All the values in this table must be less than \a current.
00568 
00569 \param code
00570 The input value of code does not matter.
00571 The output value of code is the hash code corresponding to
00572 this operation in the new operation sequence.
00573 
00574 \return
00575 If the return value is zero,
00576 no match was found.
00577 If the return value is greater than zero,
00578 it is the index of a new variable that can be used to replace the 
00579 old variable.
00580 
00581 
00582 \par Restrictions:
00583 The binary operator must be an addition, subtraction, multiplication, division
00584 or power operator.  NumArg( tape[current].op ) == 1.
00585 */
00586 template <class Base>
00587 inline size_t binary_match(
00588      const CppAD::vector<struct struct_old_variable>& tape           ,
00589      size_t                                             current        ,
00590      size_t                                             npar           ,
00591      const Base*                                        par            ,
00592      const CppAD::vector<size_t>&                       hash_table_var ,
00593      unsigned short&                                    code           )
00594 {    OpCode        op         = tape[current].op;
00595      const addr_t* arg        = tape[current].arg;
00596      addr_t        new_arg[2];
00597      bool          parameter[2];
00598 
00599      // initialize return value
00600      size_t  match_var = 0;
00601 
00602      CPPAD_ASSERT_UNKNOWN( NumArg(op) == 2 );
00603      CPPAD_ASSERT_UNKNOWN( NumRes(op) >  0 );
00604      switch(op)
00605      {    // parameter op variable ----------------------------------
00606           case AddpvOp:
00607           case MulpvOp:
00608           case DivpvOp:
00609           case PowpvOp:
00610           case SubpvOp:
00611           // arg[0]
00612           parameter[0] = true;
00613           new_arg[0]   = arg[0];
00614           CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < npar );
00615           // arg[1]
00616           parameter[1] = false;
00617           new_arg[1]   = tape[arg[1]].new_var;
00618           CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < current );
00619           break;
00620 
00621           // variable op parameter -----------------------------------
00622           case DivvpOp:
00623           case PowvpOp:
00624           case SubvpOp:
00625           // arg[0]
00626           parameter[0] = false;
00627           new_arg[0]   = tape[arg[0]].new_var;
00628           CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < current );
00629           // arg[1]
00630           parameter[1] = true;
00631           new_arg[1]   = arg[1];
00632           CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < npar );
00633           break;
00634 
00635           // variable op variable -----------------------------------
00636           case AddvvOp:
00637           case MulvvOp:
00638           case DivvvOp:
00639           case PowvvOp:
00640           case SubvvOp:
00641           // arg[0]
00642           parameter[0] = false;
00643           new_arg[0]   = tape[arg[0]].new_var;
00644           CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < current );
00645           // arg[1]
00646           parameter[1] = false;
00647           new_arg[1]   = tape[arg[1]].new_var;
00648           CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < current );
00649           break;
00650 
00651           // must be one of the cases above
00652           default:
00653           CPPAD_ASSERT_UNKNOWN(false);
00654      }
00655      code = hash_code(
00656           op                  , 
00657           new_arg             ,
00658           npar                ,
00659           par
00660      );
00661      size_t  i  = hash_table_var[code];
00662      CPPAD_ASSERT_UNKNOWN( i < current );
00663      if( op == tape[i].op )
00664      {    bool match = true;
00665           size_t j;
00666           for(j = 0; j < 2; j++)
00667           {    size_t k = tape[i].arg[j];
00668                if( parameter[j] )
00669                {    CPPAD_ASSERT_UNKNOWN( k < npar );
00670                     match &= IdenticalEqualPar(
00671                          par[ arg[j] ], par[k]
00672                     );
00673                }
00674                else
00675                {    CPPAD_ASSERT_UNKNOWN( k < i );
00676                     match &= (new_arg[j] == tape[k].new_var);
00677                }
00678           }
00679           if( match )
00680                match_var = tape[i].new_var;
00681      }
00682      if( (match_var > 0) | ( (op != AddvvOp) & (op != MulvvOp ) ) )
00683           return match_var;
00684 
00685      // check for match with argument order switched ----------------------
00686      CPPAD_ASSERT_UNKNOWN( op == AddvvOp || op == MulvvOp );
00687      i          = new_arg[0];
00688      new_arg[0] = new_arg[1];
00689      new_arg[1] = i;
00690      unsigned short code_switch = hash_code(
00691           op                  , 
00692           new_arg             ,
00693           npar                ,
00694           par
00695      );
00696      i  = hash_table_var[code_switch];
00697      CPPAD_ASSERT_UNKNOWN( i < current );
00698      if( op == tape[i].op )
00699      {    bool match = true;
00700           size_t j;
00701           for(j = 0; j < 2; j++)
00702           {    size_t k = tape[i].arg[j];
00703                CPPAD_ASSERT_UNKNOWN( k < i );
00704                match &= (new_arg[j] == tape[k].new_var);
00705           }
00706           if( match )
00707                match_var = tape[i].new_var;
00708      }
00709      return match_var;
00710 } 
00711 
00712 /*!
00713 Record an operation of the form (parameter op variable).
00714 
00715 <!-- replace prototype -->
00716 \param tape
00717 is a vector that maps a variable index, in the old operation sequence,
00718 to an <tt>struct_old_variable</tt> information record.
00719 Note that the index for this vector must be greater than or equal zero and 
00720 less than <tt>tape.size()</tt>.
00721 
00722 \li <tt>tape[i].op</tt> 
00723 is the operator in the old operation sequence
00724 corresponding to the old variable index \c i.
00725 Assertion: <tt>NumRes(tape[i].op) > 0</tt>.
00726 
00727 \li <tt>tape[i].arg</tt> 
00728 for <tt>j < NumArg( tape[i].op ), tape[i].arg[j]</tt>
00729 is the j-th the argument, in the old operation sequence,
00730 corresponding to the old variable index \c i.
00731 Assertion: <tt>tape[i].arg[j] < i</tt>.
00732 
00733 \li <tt>tape[i].new_var</tt>
00734 Suppose 
00735 <tt>i <= current, j < NumArg( tape[i].op ), and k = tape[i].arg[j]</tt>,
00736 and \c j corresponds to a variable for operator <tt>tape[i].op</tt>.
00737 It follows that <tt>tape[k].new_var</tt>
00738 has alread been set to the variable in the new operation sequence 
00739 corresponding to the old variable index \c k.
00740 This means that the \c new_var value has been set
00741 for all the possible arguments that come before \a current.
00742 
00743 \param current
00744 is the index in the old operation sequence for 
00745 the variable corresponding to the result for the current operator.
00746 Assertions: 
00747 <tt>
00748 current < tape.size(),
00749 NumRes( tape[current].op ) > 0.
00750 </tt>
00751 
00752 \param npar
00753 is the number of parameters corresponding to this operation sequence.
00754 
00755 \param par
00756 is a vector of length \a npar containing the parameters
00757 for this operation sequence; i.e.,
00758 given a parameter index \c i, the corresponding parameter value is
00759 <tt>par[i]</tt>.
00760 <!-- end prototype -->
00761 
00762 \param rec
00763 is the object that will record the operations.
00764 
00765 \param op
00766 is the operator that we are recording which must be one of the following:
00767 AddpvOp, DivpvOp, MulpvOp, PowvpOp, SubpvOp.
00768  
00769 \param arg
00770 is the vector of arguments for this operator.
00771 
00772 \return
00773 the result is the operaiton and variable index corresponding to the current
00774 operation in the new operation sequence.
00775 */
00776 template <class Base>
00777 struct_size_pair record_pv(
00778      const CppAD::vector<struct struct_old_variable>& tape           ,
00779      size_t                                             current        ,
00780      size_t                                             npar           ,
00781      const Base*                                        par            ,
00782      recorder<Base>*                                    rec            ,
00783      OpCode                                             op             ,
00784      const addr_t*                                      arg            )
00785 {
00786 # ifndef NDEBUG
00787      switch(op)
00788      {    case AddpvOp:
00789           case DivpvOp:
00790           case MulpvOp:
00791           case PowpvOp:
00792           case SubpvOp:
00793           break;
00794 
00795           default:
00796           CPPAD_ASSERT_UNKNOWN(false);
00797      }
00798 # endif
00799      CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < npar    );
00800      CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < current );
00801      addr_t new_arg[2];
00802      new_arg[0]   = rec->PutPar( par[arg[0]] );
00803      new_arg[1]   = tape[ arg[1] ].new_var;
00804      rec->PutArg( new_arg[0], new_arg[1] );
00805 
00806      struct_size_pair ret;
00807      ret.i_op  = rec->num_op_rec();
00808      ret.i_var = rec->PutOp(op);
00809      CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < ret.i_var );
00810      return ret;
00811 }
00812 
00813 
00814 /*!
00815 Record an operation of the form (variable op parameter).
00816 
00817 <!-- replace prototype -->
00818 \param tape
00819 is a vector that maps a variable index, in the old operation sequence,
00820 to an <tt>struct_old_variable</tt> information record.
00821 Note that the index for this vector must be greater than or equal zero and 
00822 less than <tt>tape.size()</tt>.
00823 
00824 \li <tt>tape[i].op</tt> 
00825 is the operator in the old operation sequence
00826 corresponding to the old variable index \c i.
00827 Assertion: <tt>NumRes(tape[i].op) > 0</tt>.
00828 
00829 \li <tt>tape[i].arg</tt> 
00830 for <tt>j < NumArg( tape[i].op ), tape[i].arg[j]</tt>
00831 is the j-th the argument, in the old operation sequence,
00832 corresponding to the old variable index \c i.
00833 Assertion: <tt>tape[i].arg[j] < i</tt>.
00834 
00835 \li <tt>tape[i].new_var</tt>
00836 Suppose 
00837 <tt>i <= current, j < NumArg( tape[i].op ), and k = tape[i].arg[j]</tt>,
00838 and \c j corresponds to a variable for operator <tt>tape[i].op</tt>.
00839 It follows that <tt>tape[k].new_var</tt>
00840 has alread been set to the variable in the new operation sequence 
00841 corresponding to the old variable index \c k.
00842 This means that the \c new_var value has been set
00843 for all the possible arguments that come before \a current.
00844 
00845 \param current
00846 is the index in the old operation sequence for 
00847 the variable corresponding to the result for the current operator.
00848 Assertions: 
00849 <tt>
00850 current < tape.size(),
00851 NumRes( tape[current].op ) > 0.
00852 </tt>
00853 
00854 \param npar
00855 is the number of parameters corresponding to this operation sequence.
00856 
00857 \param par
00858 is a vector of length \a npar containing the parameters
00859 for this operation sequence; i.e.,
00860 given a parameter index \c i, the corresponding parameter value is
00861 <tt>par[i]</tt>.
00862 <!-- end prototype -->
00863 
00864 \param rec
00865 is the object that will record the operations.
00866 
00867 \param op
00868 is the operator that we are recording which must be one of the following:
00869 DivvpOp, PowvpOp, SubvpOp.
00870  
00871 \param arg
00872 is the vector of arguments for this operator.
00873 
00874 \return
00875 the result operation and variable index corresponding to the current
00876 operation in the new operation sequence.
00877 */
00878 template <class Base>
00879 struct_size_pair record_vp(
00880      const CppAD::vector<struct struct_old_variable>& tape           ,
00881      size_t                                             current        ,
00882      size_t                                             npar           ,
00883      const Base*                                        par            ,
00884      recorder<Base>*                                    rec            ,
00885      OpCode                                             op             ,
00886      const addr_t*                                      arg            )
00887 {
00888 # ifndef NDEBUG
00889      switch(op)
00890      {    case DivvpOp:
00891           case PowvpOp:
00892           case SubvpOp:
00893           break;
00894 
00895           default:
00896           CPPAD_ASSERT_UNKNOWN(false);
00897      }
00898 # endif
00899      CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < current );
00900      CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < npar    );
00901      addr_t new_arg[2];
00902      new_arg[0]   = tape[ arg[0] ].new_var;
00903      new_arg[1]   = rec->PutPar( par[arg[1]] );
00904      rec->PutArg( new_arg[0], new_arg[1] );
00905 
00906      struct_size_pair ret;
00907      ret.i_op  = rec->num_op_rec();
00908      ret.i_var = rec->PutOp(op);
00909      CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < ret.i_var );
00910      return ret;
00911 }
00912 
00913 /*!
00914 Record an operation of the form (variable op variable).
00915 
00916 <!-- replace prototype -->
00917 \param tape
00918 is a vector that maps a variable index, in the old operation sequence,
00919 to an <tt>struct_old_variable</tt> information record.
00920 Note that the index for this vector must be greater than or equal zero and 
00921 less than <tt>tape.size()</tt>.
00922 
00923 \li <tt>tape[i].op</tt> 
00924 is the operator in the old operation sequence
00925 corresponding to the old variable index \c i.
00926 Assertion: <tt>NumRes(tape[i].op) > 0</tt>.
00927 
00928 \li <tt>tape[i].arg</tt> 
00929 for <tt>j < NumArg( tape[i].op ), tape[i].arg[j]</tt>
00930 is the j-th the argument, in the old operation sequence,
00931 corresponding to the old variable index \c i.
00932 Assertion: <tt>tape[i].arg[j] < i</tt>.
00933 
00934 \li <tt>tape[i].new_var</tt>
00935 Suppose 
00936 <tt>i <= current, j < NumArg( tape[i].op ), and k = tape[i].arg[j]</tt>,
00937 and \c j corresponds to a variable for operator <tt>tape[i].op</tt>.
00938 It follows that <tt>tape[k].new_var</tt>
00939 has alread been set to the variable in the new operation sequence 
00940 corresponding to the old variable index \c k.
00941 This means that the \c new_var value has been set
00942 for all the possible arguments that come before \a current.
00943 
00944 \param current
00945 is the index in the old operation sequence for 
00946 the variable corresponding to the result for the current operator.
00947 Assertions: 
00948 <tt>
00949 current < tape.size(),
00950 NumRes( tape[current].op ) > 0.
00951 </tt>
00952 
00953 \param npar
00954 is the number of parameters corresponding to this operation sequence.
00955 
00956 \param par
00957 is a vector of length \a npar containing the parameters
00958 for this operation sequence; i.e.,
00959 given a parameter index \c i, the corresponding parameter value is
00960 <tt>par[i]</tt>.
00961 <!-- end prototype -->
00962 
00963 \param rec
00964 is the object that will record the operations.
00965 
00966 \param op
00967 is the operator that we are recording which must be one of the following:
00968 AddvvOp, DivvvOp, MulvvOp, PowvpOp, SubvvOp.
00969  
00970 \param arg
00971 is the vector of arguments for this operator.
00972 
00973 \return
00974 the result is the operation and variable index corresponding to the current
00975 operation in the new operation sequence.
00976 */
00977 template <class Base>
00978 struct_size_pair record_vv(
00979      const CppAD::vector<struct struct_old_variable>& tape           ,
00980      size_t                                             current        ,
00981      size_t                                             npar           ,
00982      const Base*                                        par            ,
00983      recorder<Base>*                                    rec            ,
00984      OpCode                                             op             ,
00985      const addr_t*                                      arg            )
00986 {
00987 # ifndef NDEBUG
00988      switch(op)
00989      {    case AddvvOp:
00990           case DivvvOp:
00991           case MulvvOp:
00992           case PowvvOp:
00993           case SubvvOp:
00994           break;
00995 
00996           default:
00997           CPPAD_ASSERT_UNKNOWN(false);
00998      }
00999 # endif
01000      CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < current );
01001      CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < current );
01002      addr_t new_arg[2];
01003      new_arg[0]   = tape[ arg[0] ].new_var;
01004      new_arg[1]   = tape[ arg[1] ].new_var;
01005      rec->PutArg( new_arg[0], new_arg[1] );
01006 
01007      struct_size_pair ret;
01008      ret.i_op  = rec->num_op_rec();
01009      ret.i_var = rec->PutOp(op);
01010      CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < ret.i_var );
01011      CPPAD_ASSERT_UNKNOWN( size_t(new_arg[1]) < ret.i_var );
01012      return ret;
01013 }
01014 
01015 // ==========================================================================
01016 
01017 /*!
01018 Recording a cummulative cummulative summation starting at its highest parrent.
01019 
01020 <!-- replace prototype -->
01021 \param tape
01022 is a vector that maps a variable index, in the old operation sequence,
01023 to an <tt>struct_old_variable</tt> information record.
01024 Note that the index for this vector must be greater than or equal zero and 
01025 less than <tt>tape.size()</tt>.
01026 
01027 \li <tt>tape[i].op</tt> 
01028 is the operator in the old operation sequence
01029 corresponding to the old variable index \c i.
01030 Assertion: <tt>NumRes(tape[i].op) > 0</tt>.
01031 
01032 \li <tt>tape[i].arg</tt> 
01033 for <tt>j < NumArg( tape[i].op ), tape[i].arg[j]</tt>
01034 is the j-th the argument, in the old operation sequence,
01035 corresponding to the old variable index \c i.
01036 Assertion: <tt>tape[i].arg[j] < i</tt>.
01037 
01038 \li <tt>tape[i].new_var</tt>
01039 Suppose 
01040 <tt>i <= current, j < NumArg( tape[i].op ), and k = tape[i].arg[j]</tt>,
01041 and \c j corresponds to a variable for operator <tt>tape[i].op</tt>.
01042 It follows that <tt>tape[k].new_var</tt>
01043 has alread been set to the variable in the new operation sequence 
01044 corresponding to the old variable index \c k.
01045 This means that the \c new_var value has been set
01046 for all the possible arguments that come before \a current.
01047 
01048 \param current
01049 is the index in the old operation sequence for 
01050 the variable corresponding to the result for the current operator.
01051 Assertions: 
01052 <tt>
01053 current < tape.size(),
01054 NumRes( tape[current].op ) > 0.
01055 </tt>
01056 
01057 \param npar
01058 is the number of parameters corresponding to this operation sequence.
01059 
01060 \param par
01061 is a vector of length \a npar containing the parameters
01062 for this operation sequence; i.e.,
01063 given a parameter index \c i, the corresponding parameter value is
01064 <tt>par[i]</tt>.
01065 <!-- end prototype -->
01066 
01067 \param rec
01068 is the object that will record the operations.
01069 
01070 \param work
01071 Is used for computation. On input and output,
01072 <tt>work.op_stack.empty()</tt>,
01073 <tt>work.add_stack.empty()</tt>, and
01074 <tt>work.sub_stack.empty()</tt>,
01075 are all true true.
01076 These stacks are passed in so that elements can be allocated once
01077 and then the elements can be reused with calls to \c record_csum.
01078 
01079 \par Exception
01080 <tt>tape[i].new_var</tt>
01081 is not yet defined for any node \c i that is \c csum_connected
01082 to the \a current node
01083 (or that is \c sum_connected to a node that is \c csum_connected).
01084 For example; suppose that index \c j corresponds to a variable
01085 in the current operator,
01086 <tt>i = tape[current].arg[j]</tt>,
01087 and 
01088 <tt>tape[arg[j]].connect_type == csum_connected</tt>.
01089 It then follows that
01090 <tt>tape[i].new_var == tape.size()</tt>.
01091 
01092 \par Restriction:
01093 \li <tt>tape[current].op</tt> 
01094 must be one of <tt>AddpvOp, AddvvOp, SubpvOp, SubvpOp, SubvvOp</tt>.
01095 
01096 \li <tt>tape[current].connect_type</tt> must be \c yes_connected.
01097 
01098 \li <tt>tape[j].connect_type == csum_connected</tt> for some index
01099 j that is a variable operand for the current operation.
01100 */
01101 
01102 
01103 template <class Base>
01104 struct_size_pair record_csum(
01105      const CppAD::vector<struct struct_old_variable>& tape           ,
01106      size_t                                             current        ,
01107      size_t                                             npar           ,
01108      const Base*                                        par            ,
01109      recorder<Base>*                                    rec            ,
01110      struct_csum_stacks&                              work           )
01111 {
01112      
01113      CPPAD_ASSERT_UNKNOWN( work.op_stack.empty() );
01114      CPPAD_ASSERT_UNKNOWN( work.add_stack.empty() );
01115      CPPAD_ASSERT_UNKNOWN( work.sub_stack.empty() );
01116      CPPAD_ASSERT_UNKNOWN( tape[current].connect_type == yes_connected );
01117 
01118      size_t                        i;
01119      OpCode                        op;
01120      const addr_t*                 arg;
01121      bool                          add;
01122      struct struct_csum_variable var;
01123 
01124      var.op  = tape[current].op;
01125      var.arg = tape[current].arg;
01126      var.add = true; 
01127      work.op_stack.push( var );
01128      Base sum_par(0);
01129 
01130 # ifndef NDEBUG
01131      bool ok = false;
01132      if( var.op == SubvpOp )
01133           ok = tape[ tape[current].arg[0] ].connect_type == csum_connected;
01134      if( var.op == AddpvOp || var.op == SubpvOp )
01135           ok = tape[ tape[current].arg[1] ].connect_type == csum_connected;
01136      if( var.op == AddvvOp || var.op == SubvvOp )
01137      {    ok  = tape[ tape[current].arg[0] ].connect_type == csum_connected;
01138           ok |= tape[ tape[current].arg[1] ].connect_type == csum_connected;
01139      }
01140      CPPAD_ASSERT_UNKNOWN( ok );
01141 # endif
01142      while( ! work.op_stack.empty() )
01143      {    var     = work.op_stack.top();
01144           work.op_stack.pop();
01145           op      = var.op;
01146           arg     = var.arg;
01147           add     = var.add;
01148           // process first argument to this operator
01149           switch(op)
01150           {    case AddpvOp:
01151                case SubpvOp:
01152                CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < npar );
01153                if( add )
01154                     sum_par += par[arg[0]];
01155                else sum_par -= par[arg[0]];
01156                break;
01157 
01158                case AddvvOp:
01159                case SubvpOp:
01160                case SubvvOp:
01161                if( tape[arg[0]].connect_type == csum_connected )
01162                {    CPPAD_ASSERT_UNKNOWN(
01163                          size_t(tape[arg[0]].new_var) == tape.size()
01164                     );
01165                     var.op  = tape[arg[0]].op;
01166                     var.arg = tape[arg[0]].arg;
01167                     var.add = add; 
01168                     work.op_stack.push( var );
01169                }
01170                else if( add )
01171                     work.add_stack.push(arg[0]);
01172                else work.sub_stack.push(arg[0]);
01173                break;
01174 
01175                default:
01176                CPPAD_ASSERT_UNKNOWN(false);
01177           }
01178           // process second argument to this operator
01179           switch(op)
01180           {
01181                case SubvpOp:
01182                CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < npar );
01183                if( add )
01184                     sum_par -= par[arg[1]];
01185                else sum_par += par[arg[1]];
01186                break;
01187 
01188                case SubvvOp:
01189                case SubpvOp:
01190                add = ! add;
01191 
01192                case AddvvOp:
01193                case AddpvOp:
01194                if( tape[arg[1]].connect_type == csum_connected )
01195                {    CPPAD_ASSERT_UNKNOWN(
01196                          size_t(tape[arg[1]].new_var) == tape.size()
01197                     );
01198                     var.op   = tape[arg[1]].op;
01199                     var.arg  = tape[arg[1]].arg;
01200                     var.add  = add;
01201                     work.op_stack.push( var );
01202                }
01203                else if( add )
01204                     work.add_stack.push(arg[1]);
01205                else work.sub_stack.push(arg[1]);
01206                break;
01207 
01208                default:
01209                CPPAD_ASSERT_UNKNOWN(false);
01210           }
01211      }
01212      // number of variables in this cummulative sum operator
01213      size_t n_add = work.add_stack.size();
01214      size_t n_sub = work.sub_stack.size();
01215      size_t old_arg, new_arg;
01216      rec->PutArg(n_add);                // arg[0]
01217      rec->PutArg(n_sub);                // arg[1]
01218      new_arg = rec->PutPar( sum_par );
01219      rec->PutArg(new_arg);              // arg[2]
01220      for(i = 0; i < n_add; i++)
01221      {    CPPAD_ASSERT_UNKNOWN( ! work.add_stack.empty() );
01222           old_arg = work.add_stack.top();
01223           new_arg = tape[old_arg].new_var;
01224           CPPAD_ASSERT_UNKNOWN( new_arg < tape.size() );
01225           rec->PutArg(new_arg);      // arg[3+i]
01226           work.add_stack.pop();
01227      }
01228      for(i = 0; i < n_sub; i++)
01229      {    CPPAD_ASSERT_UNKNOWN( ! work.sub_stack.empty() );
01230           old_arg = work.sub_stack.top();
01231           new_arg = tape[old_arg].new_var;
01232           CPPAD_ASSERT_UNKNOWN( new_arg < tape.size() );
01233           rec->PutArg(new_arg);      // arg[3 + arg[0] + i]
01234           work.sub_stack.pop();
01235      }
01236      rec->PutArg(n_add + n_sub);        // arg[3 + arg[0] + arg[1]]
01237 
01238 
01239      struct_size_pair ret;
01240      ret.i_op  = rec->num_op_rec();
01241      ret.i_var = rec->PutOp(CSumOp);
01242      CPPAD_ASSERT_UNKNOWN( new_arg < ret.i_var );
01243      return ret;
01244 }
01245 // ==========================================================================
01246 /*!
01247 Convert a player object to an optimized recorder object
01248 
01249 \tparam Base
01250 base type for the operator; i.e., this operation was recorded
01251 using AD< \a Base > and computations by this routine are done using type 
01252 \a Base.
01253 
01254 \param n
01255 is the number of independent variables on the tape.
01256 
01257 \param dep_taddr
01258 On input this vector contains the indices for each of the dependent
01259 variable values in the operation sequence corresponding to \a play.
01260 Upon return it contains the indices for the same variables but in
01261 the operation sequence corresponding to \a rec.
01262 
01263 \param play
01264 This is the operation sequence that we are optimizing.
01265 It is essentially const, except for play back state which
01266 changes while it plays back the operation seqeunce.
01267 
01268 \param rec
01269 The input contents of this recording does not matter.
01270 Upon return, it contains an optimized verison of the
01271 operation sequence corresponding to \a play.
01272 */
01273 
01274 template <class Base>
01275 void optimize_run(
01276      size_t                       n         ,
01277      CppAD::vector<size_t>&       dep_taddr ,
01278      player<Base>*                play      ,
01279      recorder<Base>*              rec       ) 
01280 {
01281      // temporary indices
01282      size_t i, j, k;
01283 
01284      // temporary variables
01285      OpCode        op;   // current operator
01286      const addr_t* arg;  // operator arguments
01287      size_t        i_var;  // index of first result for current operator
01288 
01289      // range and domain dimensions for F
01290      size_t m = dep_taddr.size();
01291 
01292      // number of variables in the player
01293      const size_t num_var = play->num_var_rec(); 
01294 
01295 # ifndef NDEBUG
01296      // number of parameters in the player
01297      const size_t num_par = play->num_par_rec();
01298 # endif
01299 
01300      // number of  VecAD indices 
01301      size_t num_vecad_ind   = play->num_vec_ind_rec();
01302 
01303      // number of VecAD vectors
01304      size_t num_vecad_vec   = play->num_vecad_vec_rec();
01305 
01306      // -------------------------------------------------------------
01307      // data structure that maps variable index in original operation
01308      // sequence to corresponding operator information
01309      CppAD::vector<struct struct_old_variable> tape(num_var);
01310      // -------------------------------------------------------------
01311      // Determine how each variable is connected to the dependent variables
01312 
01313      // initialize all variables has having no connections
01314      for(i = 0; i < num_var; i++)
01315           tape[i].connect_type = not_connected;
01316 
01317      for(j = 0; j < m; j++)
01318      {    // mark dependent variables as having one or more connections
01319           tape[ dep_taddr[j] ].connect_type = yes_connected;
01320      }
01321 
01322      // vecad_connect contains a value for each VecAD object.
01323      // vecad maps a VecAD index (which corresponds to the beginning of the
01324      // VecAD object) to the vecad_connect falg for the VecAD object.
01325      CppAD::vector<enum_connect_type>   vecad_connect(num_vecad_vec);
01326      CppAD::vector<size_t> vecad(num_vecad_ind);
01327      j = 0;
01328      for(i = 0; i < num_vecad_vec; i++)
01329      {    vecad_connect[i] = not_connected;
01330           // length of this VecAD
01331           size_t length = play->GetVecInd(j);
01332           // set to proper index for this VecAD
01333           vecad[j] = i; 
01334           for(k = 1; k <= length; k++)
01335                vecad[j+k] = num_vecad_vec; // invalid index
01336           // start of next VecAD
01337           j       += length + 1;
01338      }
01339      CPPAD_ASSERT_UNKNOWN( j == num_vecad_ind );
01340 
01341      // work space used by UserOp.
01342      typedef std::set<size_t> size_set;
01343      vector<size_set> user_r_set;   // set sparsity pattern for result
01344      vector<size_set> user_s_set;   // set sparisty pattern for argument
01345      vector<bool>     user_r_bool;  // bool sparsity pattern for result 
01346      vector<bool>     user_s_bool;  // bool sparisty pattern for argument
01347      //
01348      size_t user_q     = 0;       // column dimension for sparsity patterns
01349      size_t user_index = 0;       // indentifier for this user_atomic operation
01350      size_t user_id    = 0;       // user identifier for this call to operator
01351      size_t user_i     = 0;       // index in result vector
01352      size_t user_j     = 0;       // index in argument vector
01353      size_t user_m     = 0;       // size of result vector
01354      size_t user_n     = 0;       // size of arugment vector
01355      //
01356      atomic_base<Base>* user_atom = CPPAD_NULL; // current user atomic function 
01357      bool               user_set  = true;       // use set sparsity (or bool)
01358 
01359      // next expected operator in a UserOp sequence
01360      enum { user_start, user_arg, user_ret, user_end } user_state;
01361 
01362      // During reverse mode, compute type of connection for each call to
01363      // a user atomic function.
01364      CppAD::vector<struct_user_info>    user_info;
01365      size_t                             user_curr = 0;
01366 
01367      /// During reverse mode, information for each CSkip operation
01368      CppAD::vector<struct_cskip_info>   cskip_info;
01369 
01370      // Initialize a reverse mode sweep through the operation sequence
01371      size_t i_op;
01372      play->reverse_start(op, arg, i_op, i_var);
01373      CPPAD_ASSERT_UNKNOWN( op == EndOp );
01374      size_t mask;
01375      user_state = user_end;
01376      while(op != BeginOp)
01377      {    // next op
01378           play->reverse_next(op, arg, i_op, i_var);
01379 
01380           // Store the operator corresponding to each variable
01381           if( NumRes(op) > 0 )
01382           {    tape[i_var].op = op;
01383                tape[i_var].arg = arg;
01384           }
01385 # ifndef NDEBUG
01386           if( i_op <= n )
01387           {    CPPAD_ASSERT_UNKNOWN((op == InvOp) | (op == BeginOp));
01388           }
01389           else CPPAD_ASSERT_UNKNOWN((op != InvOp) & (op != BeginOp));
01390 # endif
01391           enum_connect_type connect_type      = tape[i_var].connect_type;
01392           std::set<class_cexp_pair>& cexp_set = tape[i_var].cexp_set;
01393           switch( op )
01394           {
01395                // Unary operator where operand is arg[0]
01396                case AbsOp:
01397                case AcosOp:
01398                case AsinOp:
01399                case AtanOp:
01400                case CosOp:
01401                case CoshOp:
01402                case DivvpOp:
01403                case ExpOp:
01404                case LogOp:
01405                case PowvpOp:
01406                case SignOp:
01407                case SinOp:
01408                case SinhOp:
01409                case SqrtOp:
01410                case TanOp:
01411                case TanhOp:
01412                switch( connect_type )
01413                {    case not_connected:
01414                     break;
01415      
01416                     case yes_connected:
01417                     case sum_connected:
01418                     case csum_connected: 
01419                     tape[arg[0]].connect_type = yes_connected;
01420                     break;
01421 
01422                     case cexp_connected:
01423                     if( tape[arg[0]].connect_type == not_connected )
01424                     {    tape[arg[0]].connect_type = cexp_connected;
01425                          tape[arg[0]].cexp_set     = cexp_set;
01426                     }
01427                     else if( tape[arg[0]].connect_type == cexp_connected )
01428                     {    tape[arg[0]].cexp_set = intersection(
01429                               tape[arg[0]].cexp_set, cexp_set
01430                          );
01431                          if( tape[arg[0]].cexp_set.empty() )
01432                               tape[arg[0]].connect_type = yes_connected;
01433                     }
01434                     else tape[arg[0]].connect_type = yes_connected;
01435                     break;
01436 
01437                     default:
01438                     CPPAD_ASSERT_UNKNOWN(false);
01439                }
01440                break; // --------------------------------------------
01441 
01442                // Unary operator where operand is arg[1]
01443                case DisOp:
01444                case DivpvOp:
01445                case MulpvOp:
01446                case PowpvOp:
01447                switch( connect_type )
01448                {    case not_connected:
01449                     break;
01450      
01451                     case yes_connected:
01452                     case sum_connected:
01453                     case csum_connected: 
01454                     tape[arg[1]].connect_type = yes_connected;
01455                     break;
01456 
01457                     case cexp_connected:
01458                     if( tape[arg[1]].connect_type == not_connected )
01459                     {    tape[arg[1]].connect_type = cexp_connected;
01460                          tape[arg[1]].cexp_set     = cexp_set;
01461                     }
01462                     else if( tape[arg[1]].connect_type == cexp_connected )
01463                     {    tape[arg[1]].cexp_set = intersection(
01464                               tape[arg[1]].cexp_set, cexp_set
01465                          );
01466                          if( tape[arg[1]].cexp_set.empty() )
01467                               tape[arg[1]].connect_type = yes_connected;
01468                     }
01469                     else tape[arg[1]].connect_type = yes_connected;
01470                     break;
01471 
01472                     default:
01473                     CPPAD_ASSERT_UNKNOWN(false);
01474                }
01475                break; // --------------------------------------------
01476           
01477                // Special case for SubvpOp
01478                case SubvpOp:
01479                switch( connect_type )
01480                {    case not_connected:
01481                     break;
01482      
01483                     case yes_connected:
01484                     case sum_connected:
01485                     case csum_connected: 
01486                     if( tape[arg[0]].connect_type == not_connected )
01487                          tape[arg[0]].connect_type = sum_connected;
01488                     else tape[arg[0]].connect_type = yes_connected;
01489                     break;
01490 
01491                     case cexp_connected:
01492                     if( tape[arg[0]].connect_type == not_connected )
01493                     {    tape[arg[0]].connect_type = cexp_connected;
01494                          tape[arg[0]].cexp_set     = cexp_set;
01495                     }
01496                     else if( tape[arg[0]].connect_type == cexp_connected )
01497                     {    tape[arg[0]].cexp_set = intersection(
01498                               tape[arg[0]].cexp_set, cexp_set
01499                          );
01500                          if( tape[arg[0]].cexp_set.empty() )
01501                               tape[arg[0]].connect_type = yes_connected;
01502                     }
01503                     else tape[arg[0]].connect_type = yes_connected;
01504                     break;
01505 
01506                     default:
01507                     CPPAD_ASSERT_UNKNOWN(false);
01508                }
01509                if( connect_type == sum_connected )
01510                {    // convert sum to csum connection for this variable
01511                     tape[i_var].connect_type = connect_type = csum_connected;
01512                }
01513                break; // --------------------------------------------
01514           
01515                // Special case for AddpvOp and SubpvOp
01516                case AddpvOp:
01517                case SubpvOp:
01518                switch( connect_type )
01519                {    case not_connected:
01520                     break;
01521      
01522                     case yes_connected:
01523                     case sum_connected:
01524                     case csum_connected:
01525                     if( tape[arg[1]].connect_type == not_connected )
01526                          tape[arg[1]].connect_type = sum_connected;
01527                     else tape[arg[1]].connect_type = yes_connected;
01528                     break;
01529 
01530                     case cexp_connected:
01531                     if( tape[arg[1]].connect_type == not_connected )
01532                     {    tape[arg[1]].connect_type = cexp_connected;
01533                          tape[arg[1]].cexp_set     = cexp_set;
01534                     }
01535                     else if( tape[arg[1]].connect_type == cexp_connected )
01536                     {    tape[arg[1]].cexp_set = intersection(
01537                               tape[arg[1]].cexp_set, cexp_set
01538                          );
01539                          if( tape[arg[1]].cexp_set.empty() )
01540                               tape[arg[1]].connect_type = yes_connected;
01541                     }
01542                     else tape[arg[1]].connect_type = yes_connected;
01543                     break;
01544 
01545                     default:
01546                     CPPAD_ASSERT_UNKNOWN(false);
01547                }
01548                if( connect_type == sum_connected )
01549                {    // convert sum to csum connection for this variable
01550                     tape[i_var].connect_type = connect_type = csum_connected;
01551                }
01552                break; // --------------------------------------------
01553 
01554           
01555                // Special case for AddvvOp and SubvvOp
01556                case AddvvOp:
01557                case SubvvOp:
01558                for(i = 0; i < 2; i++) switch( connect_type )
01559                {    case not_connected:
01560                     break;
01561 
01562                     case yes_connected:
01563                     case sum_connected:
01564                     case csum_connected:
01565                     if( tape[arg[i]].connect_type == not_connected )
01566                          tape[arg[i]].connect_type = sum_connected;
01567                     else tape[arg[i]].connect_type = yes_connected;
01568                     break;
01569 
01570                     case cexp_connected:
01571                     if( tape[arg[i]].connect_type == not_connected )
01572                     {    tape[arg[i]].connect_type = cexp_connected;
01573                          tape[arg[i]].cexp_set     = cexp_set;
01574                     }
01575                     else if( tape[arg[i]].connect_type == cexp_connected )
01576                     {    tape[arg[i]].cexp_set = intersection(
01577                               tape[arg[i]].cexp_set, cexp_set
01578                          );
01579                          if( tape[arg[i]].cexp_set.empty() )
01580                               tape[arg[i]].connect_type = yes_connected;
01581                     }
01582                     else tape[arg[i]].connect_type = yes_connected;
01583                     break;
01584 
01585                     default:
01586                     CPPAD_ASSERT_UNKNOWN(false);
01587                }
01588                if( connect_type == sum_connected )
01589                {    // convert sum to csum connection for this variable
01590                     tape[i_var].connect_type = connect_type = csum_connected;
01591                }
01592                break; // --------------------------------------------
01593 
01594                // Other binary operators 
01595                // where operands are arg[0], arg[1]
01596                case DivvvOp:
01597                case MulvvOp:
01598                case PowvvOp:
01599                for(i = 0; i < 2; i++) switch( connect_type )
01600                {    case not_connected:
01601                     break;
01602 
01603                     case yes_connected:
01604                     case sum_connected:
01605                     case csum_connected:
01606                     tape[arg[i]].connect_type = yes_connected;
01607                     break;
01608 
01609                     case cexp_connected:
01610                     if( tape[arg[i]].connect_type == not_connected )
01611                     {    tape[arg[i]].connect_type = cexp_connected;
01612                          tape[arg[i]].cexp_set     = cexp_set;
01613                     }
01614                     else if( tape[arg[i]].connect_type == cexp_connected )
01615                     {    tape[arg[i]].cexp_set = intersection(
01616                               tape[arg[i]].cexp_set, cexp_set
01617                          );
01618                          if( tape[arg[i]].cexp_set.empty() )
01619                               tape[arg[i]].connect_type = yes_connected;
01620                     }
01621                     else tape[arg[i]].connect_type = yes_connected;
01622                     break;
01623 
01624                     default:
01625                     CPPAD_ASSERT_UNKNOWN(false);
01626                }
01627                break; // --------------------------------------------
01628 
01629                // Conditional expression operators
01630                case CExpOp:
01631                CPPAD_ASSERT_UNKNOWN( NumArg(CExpOp) == 6 );
01632                if( connect_type != not_connected )
01633                {    struct_cskip_info info;
01634                     info.cop        = CompareOp( arg[0] );
01635                     info.flag       = arg[1];
01636                     info.left       = arg[2];
01637                     info.right      = arg[3];
01638                     info.n_op_true  = 0;
01639                     info.n_op_false = 0;
01640                     //
01641                     size_t index    = 0;
01642                     if( arg[1] & 1 )
01643                     {    index = std::max(index, info.left);
01644                          tape[info.left].connect_type = yes_connected;
01645                     }
01646                     if( arg[1] & 2 )
01647                     {    index = std::max(index, info.right);
01648                          tape[info.right].connect_type = yes_connected;
01649                     }
01650                     CPPAD_ASSERT_UNKNOWN( index > 0 );
01651                     info.max_left_right = index;
01652                     //
01653                     index = cskip_info.size();
01654                     cskip_info.push_back(info);
01655                     //
01656                     if( arg[1] & 4 )
01657                     {    tape[arg[4]].connect_type = cexp_connected;
01658                          tape[arg[4]].cexp_set     = cexp_set;
01659                          tape[arg[4]].cexp_set.insert(
01660                               class_cexp_pair(true, index)
01661                          );
01662                     }
01663                     if( arg[1] & 8 )
01664                     {    tape[arg[5]].connect_type = cexp_connected;
01665                          tape[arg[5]].cexp_set     = cexp_set;
01666                          tape[arg[5]].cexp_set.insert(
01667                               class_cexp_pair(false, index)
01668                          );
01669                     }
01670                }
01671                break;  // --------------------------------------------
01672 
01673                // Operations where there is noting to do
01674                case ComOp:
01675                case EndOp:
01676                case ParOp:
01677                case PriOp:
01678                break;  // --------------------------------------------
01679 
01680                // Operators that never get removed
01681                case BeginOp:
01682                case InvOp:
01683                tape[i_var].connect_type = yes_connected;
01684                break;
01685 
01686                // Load using a parameter index
01687                case LdpOp:
01688                if( tape[i_var].connect_type != not_connected )
01689                {
01690                     i                = vecad[ arg[0] - 1 ];
01691                     vecad_connect[i] = yes_connected;
01692                }
01693                break; // --------------------------------------------
01694 
01695                // Load using a variable index
01696                case LdvOp:
01697                if( tape[i_var].connect_type != not_connected )
01698                {
01699                     i                    = vecad[ arg[0] - 1 ];
01700                     vecad_connect[i]     = yes_connected;
01701                     tape[arg[1]].connect_type = yes_connected;
01702                }
01703                break; // --------------------------------------------
01704 
01705                // Store a variable using a parameter index
01706                case StpvOp:
01707                i = vecad[ arg[0] - 1 ];
01708                if( vecad_connect[i] != not_connected )
01709                     tape[arg[2]].connect_type = yes_connected;
01710                break; // --------------------------------------------
01711 
01712                // Store a variable using a variable index
01713                case StvvOp:
01714                i = vecad[ arg[0] - 1 ];
01715                if( vecad_connect[i] )
01716                {    tape[arg[1]].connect_type = yes_connected;
01717                     tape[arg[2]].connect_type = yes_connected;
01718                }
01719                break; 
01720                // ============================================================
01721                case UserOp:
01722                // start or end atomic operation sequence
01723                CPPAD_ASSERT_UNKNOWN( NumRes( UserOp ) == 0 );
01724                CPPAD_ASSERT_UNKNOWN( NumArg( UserOp ) == 4 );
01725                if( user_state == user_end )
01726                {    user_index = arg[0];
01727                     user_id    = arg[1];
01728                     user_n     = arg[2];
01729                     user_m     = arg[3];
01730                     user_q     = 1;
01731                     user_atom  = atomic_base<Base>::class_object(user_index);
01732                     user_set   = user_atom->sparsity() ==
01733                          atomic_base<Base>::set_sparsity_enum;
01734                     //
01735                     if(user_s_set.size() != user_n )
01736                          user_s_set.resize(user_n);
01737                     if(user_r_set.size() != user_m )
01738                          user_r_set.resize(user_m);
01739                     //
01740                     // Note user_q is 1, but use it for clarity of code
01741                     if(user_s_bool.size() != user_n * user_q )
01742                          user_s_bool.resize(user_n * user_q);
01743                     if(user_r_bool.size() != user_m * user_q )
01744                          user_r_bool.resize(user_m * user_q);
01745                     //
01746                     user_j     = user_n;
01747                     user_i     = user_m;
01748                     user_state = user_ret;
01749                     //
01750                     struct_user_info info;
01751                     info.connect_type = not_connected;
01752                     info.op_end       = i_op + 1;
01753                     user_info.push_back(info);
01754                     
01755                }
01756                else
01757                {    CPPAD_ASSERT_UNKNOWN( user_state == user_start );
01758                     CPPAD_ASSERT_UNKNOWN( user_index == size_t(arg[0]) );
01759                     CPPAD_ASSERT_UNKNOWN( user_id    == size_t(arg[1]) );
01760                     CPPAD_ASSERT_UNKNOWN( user_n     == size_t(arg[2]) );
01761                     CPPAD_ASSERT_UNKNOWN( user_m     == size_t(arg[3]) );
01762                     user_state = user_end;
01763                     //
01764                     CPPAD_ASSERT_UNKNOWN( user_curr + 1 == user_info.size() );
01765                     user_info[user_curr].op_begin = i_op;
01766                     user_curr                     = user_info.size();
01767                }
01768                break;
01769 
01770                case UsrapOp:
01771                // parameter argument in an atomic operation sequence
01772                CPPAD_ASSERT_UNKNOWN( user_state == user_arg );
01773                CPPAD_ASSERT_UNKNOWN( 0 < user_j && user_j <= user_n );
01774                CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
01775                CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
01776                --user_j;
01777                if( user_j == 0 )
01778                     user_state = user_start;
01779                break;
01780 
01781                case UsravOp:
01782                // variable argument in an atomic operation sequence
01783                CPPAD_ASSERT_UNKNOWN( user_state == user_arg );
01784                CPPAD_ASSERT_UNKNOWN( 0 < user_j && user_j <= user_n );
01785                CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
01786                CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) <= i_var );
01787                CPPAD_ASSERT_UNKNOWN( 0 < arg[0] );
01788                --user_j;
01789                if( user_set )
01790                {    if( ! user_s_set[user_j].empty() )
01791                          tape[arg[0]].connect_type = 
01792                               user_info[user_curr].connect_type;
01793                }
01794                else
01795                {    if( user_s_bool[user_j] )
01796                          tape[arg[0]].connect_type = 
01797                               user_info[user_curr].connect_type;
01798                }
01799                if( user_j == 0 )
01800                     user_state = user_start;
01801                break;
01802 
01803                case UsrrvOp:
01804                // variable result in an atomic operation sequence
01805                CPPAD_ASSERT_UNKNOWN( user_state == user_ret );
01806                CPPAD_ASSERT_UNKNOWN( 0 < user_i && user_i <= user_m );
01807                --user_i;
01808                user_r_set[user_i].clear();
01809                user_r_bool[user_i] = false;
01810                switch( connect_type )
01811                {    case not_connected:
01812                     break;
01813 
01814                     case yes_connected:
01815                     case sum_connected:
01816                     case csum_connected:
01817                     user_info[user_curr].connect_type = yes_connected;
01818                     user_r_set[user_i].insert(0);
01819                     user_r_bool[user_i] = true;
01820                     break;
01821 
01822                     case cexp_connected:
01823                     if( user_info[user_curr].connect_type == not_connected )
01824                     {    user_info[user_curr].connect_type  = connect_type;
01825                          user_info[user_curr].cexp_set      = cexp_set;
01826                     }
01827                     else if(user_info[user_curr].connect_type==cexp_connected)
01828                     {    user_info[user_curr].cexp_set = intersection(
01829                               user_info[user_curr].cexp_set, cexp_set
01830                          );
01831                          if( user_info[user_curr].cexp_set.empty() )
01832                               user_info[user_curr].connect_type = yes_connected;
01833                     }
01834                     else user_info[user_curr].connect_type = yes_connected;
01835                     break;
01836 
01837                     default:
01838                     CPPAD_ASSERT_UNKNOWN(false);
01839                }
01840                // drop into op = UsrrpOp code to handle case where user_i == 0
01841                // for both UsrrvOp and UsrrpOp together.
01842 
01843                case UsrrpOp:
01844                if( op == UsrrpOp )
01845                {    // parameter result in an atomic operation sequence
01846                     CPPAD_ASSERT_UNKNOWN( user_state == user_ret );
01847                     CPPAD_ASSERT_UNKNOWN( 0 < user_i && user_i <= user_m );
01848                     CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
01849                     CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
01850                     --user_i;
01851                     user_r_set[user_i].clear();
01852                     user_r_bool[user_i] = false;
01853                }
01854                if( user_i == 0 )
01855                {    // call users function for this operation
01856                     user_atom->set_id(user_id);
01857 # ifdef NDEBUG
01858                     if( user_set )
01859                     {    user_atom->
01860                               rev_sparse_jac(user_q, user_r_set, user_s_set);
01861                     }
01862                     else
01863                     {    user_atom->
01864                               rev_sparse_jac(user_q, user_r_bool, user_s_bool);
01865                     }
01866 # else
01867                     bool flag;
01868                     if( user_set )
01869                     {    flag = user_atom->
01870                               rev_sparse_jac(user_q, user_r_set, user_s_set);
01871                     }
01872                     else
01873                     {    flag = user_atom->
01874                               rev_sparse_jac(user_q, user_r_bool, user_s_bool);
01875                     }
01876                     if( ! flag )
01877                     {    std::string s =
01878                               "Optimizing an ADFun object"
01879                               " that contains the atomic function\n\t";
01880                          s += user_atom->afun_name();
01881                          s += "\nCurrent atomoic_sparsity is set to";
01882                          //
01883                          if( user_set )
01884                               s += " std::set\nand std::set";
01885                          else s += " bool\nand bool";
01886                          //
01887                          s += " version of rev_sparse_jac returned false";
01888                          CPPAD_ASSERT_KNOWN(false, s.c_str() );
01889                     }
01890 # endif
01891                     user_state = user_arg;
01892                }
01893                break;
01894                // ============================================================
01895 
01896                // all cases should be handled above
01897                default:
01898                CPPAD_ASSERT_UNKNOWN(0);
01899           }
01900      }
01901      // values corresponding to BeginOp
01902      CPPAD_ASSERT_UNKNOWN( i_op == 0 && i_var == 0 && op == BeginOp );
01903      tape[i_var].op = op;
01904      // -------------------------------------------------------------
01905 
01906      // Determine which variables can be conditionally skipped
01907      for(i = 0; i < num_var; i++)
01908      {    if( tape[i].connect_type == cexp_connected )
01909           {    std::set<class_cexp_pair>::const_iterator itr = 
01910                     tape[i].cexp_set.begin();
01911                while( itr != tape[i].cexp_set.end() )
01912                {    j = itr->index_;
01913                     if( itr->compare_ == true )
01914                          cskip_info[j].skip_var_false.push_back(i);
01915                     else cskip_info[j].skip_var_true.push_back(i);
01916                     itr++;
01917                }
01918           }
01919      }
01920      // Determine size of skip information in user_info 
01921      for(i = 0; i < user_info.size(); i++)
01922      {    if( user_info[i].connect_type == cexp_connected )
01923           {    std::set<class_cexp_pair>::const_iterator itr = 
01924                     user_info[i].cexp_set.begin();
01925                while( itr != user_info[i].cexp_set.end() )
01926                {    j = itr->index_;
01927                     if( itr->compare_ == true )
01928                          cskip_info[j].n_op_false = 
01929                               user_info[i].op_end - user_info[i].op_begin;
01930                     else
01931                          cskip_info[j].n_op_true = 
01932                               user_info[i].op_end - user_info[i].op_begin;
01933                     itr++;
01934                }
01935           }
01936      }
01937 
01938      // Sort the conditional skip information by the maximum of the
01939      // index for the left and right comparision operands
01940      CppAD::vector<size_t> cskip_info_order( cskip_info.size() );
01941      {    CppAD::vector<size_t> keys( cskip_info.size() );
01942           for(i = 0; i < cskip_info.size(); i++)
01943                keys[i] = std::max( cskip_info[i].left, cskip_info[i].right );
01944           CppAD::index_sort(keys, cskip_info_order);
01945      }
01946      size_t cskip_info_next = 0;
01947 
01948 
01949      // Initilaize table mapping hash code to variable index in tape
01950      // as pointing to the BeginOp at the beginning of the tape
01951      CppAD::vector<size_t>  hash_table_var(CPPAD_HASH_TABLE_SIZE);
01952      for(i = 0; i < CPPAD_HASH_TABLE_SIZE; i++)
01953           hash_table_var[i] = 0;
01954      CPPAD_ASSERT_UNKNOWN( tape[0].op == BeginOp );
01955 
01956      // initialize mapping from old variable index to new 
01957      // operator and variable index
01958      for(i = 0; i < num_var; i++)
01959      {    tape[i].new_op  = 0;       // invalid index (except for BeginOp)
01960           tape[i].new_var = num_var; // invalid index
01961      }
01962 
01963      // Erase all information in the old recording
01964      rec->free();
01965 
01966      // initialize mapping from old VecAD index to new VecAD index
01967      CppAD::vector<size_t> new_vecad_ind(num_vecad_ind);
01968      for(i = 0; i < num_vecad_ind; i++)
01969           new_vecad_ind[i] = num_vecad_ind; // invalid index 
01970 
01971      j = 0;     // index into the old set of indices
01972      for(i = 0; i < num_vecad_vec; i++)
01973      {    // length of this VecAD
01974           size_t length = play->GetVecInd(j);
01975           if( vecad_connect[i] != not_connected )
01976           {    // Put this VecAD vector in new recording
01977                CPPAD_ASSERT_UNKNOWN(length < num_vecad_ind);
01978                new_vecad_ind[j] = rec->PutVecInd(length);
01979                for(k = 1; k <= length; k++) new_vecad_ind[j+k] =
01980                     rec->PutVecInd(
01981                          rec->PutPar(
01982                               play->GetPar( 
01983                                    play->GetVecInd(j+k)
01984                ) ) );
01985           }
01986           // start of next VecAD
01987           j       += length + 1;
01988      }
01989      CPPAD_ASSERT_UNKNOWN( j == num_vecad_ind );
01990 
01991      // start playing the operations in the forward direction
01992      play->forward_start(op, arg, i_op, i_var);
01993      CPPAD_ASSERT_UNKNOWN( user_curr == user_info.size() );
01994 
01995      // playing forward skips BeginOp at the beginning, but not EndOp at
01996      // the end.  Put BeginOp at beginning of recording
01997      CPPAD_ASSERT_UNKNOWN( op == BeginOp );
01998      CPPAD_ASSERT_NARG_NRES(BeginOp, 1, 1);
01999      tape[i_var].new_op  = rec->num_op_rec();
02000      tape[i_var].new_var = rec->PutOp(BeginOp);
02001      rec->PutArg(0);
02002 
02003 
02004      // temporary buffer for new argument values
02005      addr_t new_arg[6];
02006 
02007      // temporary work space used by record_csum
02008      // (decalared here to avoid realloaction of memory)
02009      struct_csum_stacks csum_work;
02010 
02011      // tempory used to hold a size_pair
02012      struct_size_pair size_pair;
02013 
02014      user_state = user_start;
02015      while(op != EndOp)
02016      {    // next op
02017           play->forward_next(op, arg, i_op, i_var);
02018           CPPAD_ASSERT_UNKNOWN( (i_op > n)  | (op == InvOp) );
02019           CPPAD_ASSERT_UNKNOWN( (i_op <= n) | (op != InvOp) );
02020 
02021           // determine if we should insert a conditional skip here
02022           bool skip = cskip_info_next < cskip_info.size();
02023           skip     &= (op != BeginOp) & (op != InvOp);
02024           if( skip )
02025           {    j     = cskip_info_order[cskip_info_next];
02026                if( NumRes(op) > 0 )
02027                     skip &= cskip_info[j].max_left_right < i_var;
02028                else
02029                     skip &= cskip_info[j].max_left_right <= i_var;
02030           }
02031           if( skip )
02032           {    cskip_info_next++;
02033                skip &= cskip_info[j].skip_var_true.size() > 0 ||
02034                          cskip_info[j].skip_var_false.size() > 0;
02035                if( skip )
02036                {    struct_cskip_info info = cskip_info[j];
02037                     CPPAD_ASSERT_UNKNOWN( NumRes(CSkipOp) == 0 );
02038                     size_t n_true  = 
02039                          info.skip_var_true.size() + info.n_op_true;
02040                     size_t n_false = 
02041                          info.skip_var_false.size() + info.n_op_false;
02042                     size_t n_arg   = 7 + n_true + n_false; 
02043                     // reserve space for the arguments to this operator but 
02044                     // delay setting them until we have all the new addresses
02045                     cskip_info[j].i_arg = rec->ReserveArg(n_arg);
02046                     CPPAD_ASSERT_UNKNOWN( cskip_info[j].i_arg > 0 );
02047                     rec->PutOp(CSkipOp);
02048                }
02049                else cskip_info[j].i_arg = 0;
02050           }
02051 
02052           // determine if we should keep this operation in the new
02053           // operation sequence
02054           bool keep;
02055           switch( op )
02056           {    case ComOp:
02057                case PriOp:
02058                keep = false;
02059                break;
02060 
02061                case InvOp:
02062                case EndOp:
02063                keep = true;
02064                break;
02065 
02066                case StppOp:
02067                case StvpOp:
02068                case StpvOp:
02069                case StvvOp:
02070                CPPAD_ASSERT_UNKNOWN( NumRes(op) == 0 );
02071                i = vecad[ arg[0] - 1 ];
02072                keep = vecad_connect[i] != not_connected;
02073                break;
02074 
02075                case AddpvOp:
02076                case AddvvOp:
02077                case SubpvOp:
02078                case SubvpOp:
02079                case SubvvOp:
02080                keep  = tape[i_var].connect_type != not_connected;
02081                keep &= tape[i_var].connect_type != csum_connected;
02082                break; 
02083 
02084                case UserOp:
02085                case UsrapOp:
02086                case UsravOp:
02087                case UsrrpOp:
02088                case UsrrvOp:
02089                keep = true;
02090                break;
02091 
02092                default:
02093                keep = tape[i_var].connect_type != not_connected;
02094                break;
02095           }
02096 
02097           size_t         match_var    = 0;
02098           unsigned short code         = 0;
02099           bool           replace_hash = false;
02100           if( keep ) switch( op )
02101           {
02102                // Unary operator where operand is arg[0]
02103                case AbsOp:
02104                case AcosOp:
02105                case AsinOp:
02106                case AtanOp:
02107                case CosOp:
02108                case CoshOp:
02109                case ExpOp:
02110                case LogOp:
02111                case SignOp:
02112                case SinOp:
02113                case SinhOp:
02114                case SqrtOp:
02115                case TanOp:
02116                case TanhOp:
02117                match_var = unary_match(
02118                     tape                ,  // inputs 
02119                     i_var               ,
02120                     play->num_par_rec() ,
02121                     play->GetPar()      ,
02122                     hash_table_var      ,
02123                     code                  // outputs
02124                );
02125                if( match_var > 0 )
02126                     tape[i_var].new_var = match_var;
02127                else
02128                {
02129                     replace_hash = true;
02130                     new_arg[0]   = tape[ arg[0] ].new_var;
02131                     rec->PutArg( new_arg[0] );
02132                     tape[i_var].new_op  = rec->num_op_rec();
02133                     tape[i_var].new_var = i = rec->PutOp(op);
02134                     CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < i );
02135                }
02136                break;
02137                // ---------------------------------------------------
02138                // Binary operators where 
02139                // left is a variable and right is a parameter
02140                case SubvpOp:
02141                if( tape[arg[0]].connect_type == csum_connected )
02142                {
02143                     // convert to a sequence of summation operators
02144                     size_pair = record_csum(
02145                          tape                , // inputs
02146                          i_var               ,
02147                          play->num_par_rec() ,
02148                          play->GetPar()      ,
02149                          rec                 ,
02150                          csum_work
02151                     );
02152                     tape[i_var].new_op  = size_pair.i_op;
02153                     tape[i_var].new_var = size_pair.i_var;
02154                     // abort rest of this case
02155                     break;
02156                }
02157                case DivvpOp:
02158                case PowvpOp:
02159                match_var = binary_match(
02160                     tape                ,  // inputs 
02161                     i_var               ,
02162                     play->num_par_rec() ,
02163                     play->GetPar()      ,
02164                     hash_table_var      ,
02165                     code                  // outputs
02166                );
02167                if( match_var > 0 )
02168                     tape[i_var].new_var = match_var;
02169                else
02170                {    size_pair = record_vp(
02171                          tape                , // inputs
02172                          i_var               ,
02173                          play->num_par_rec() ,
02174                          play->GetPar()      ,
02175                          rec                 ,
02176                          op                  ,
02177                          arg
02178                     );
02179                     tape[i_var].new_op  = size_pair.i_op;
02180                     tape[i_var].new_var = size_pair.i_var;
02181                     replace_hash = true;
02182                }
02183                break;
02184                // ---------------------------------------------------
02185                // Binary operators where 
02186                // left is a parameter and right is a variable
02187                case SubpvOp:
02188                case AddpvOp:
02189                if( tape[arg[1]].connect_type == csum_connected )
02190                {
02191                     // convert to a sequence of summation operators
02192                     size_pair = record_csum(
02193                          tape                , // inputs
02194                          i_var               ,
02195                          play->num_par_rec() ,
02196                          play->GetPar()      ,
02197                          rec                 ,
02198                          csum_work
02199                     );
02200                     tape[i_var].new_op  = size_pair.i_op;
02201                     tape[i_var].new_var = size_pair.i_var;
02202                     // abort rest of this case
02203                     break;
02204                }
02205                case DivpvOp:
02206                case MulpvOp:
02207                case PowpvOp:
02208                match_var = binary_match(
02209                     tape                ,  // inputs 
02210                     i_var               ,
02211                     play->num_par_rec() ,
02212                     play->GetPar()      ,
02213                     hash_table_var      ,
02214                     code                  // outputs
02215                );
02216                if( match_var > 0 )
02217                     tape[i_var].new_var = match_var;
02218                else
02219                {    size_pair = record_pv(
02220                          tape                , // inputs
02221                          i_var               ,
02222                          play->num_par_rec() ,
02223                          play->GetPar()      ,
02224                          rec                 ,
02225                          op                  ,
02226                          arg
02227                     );
02228                     tape[i_var].new_op  = size_pair.i_op;
02229                     tape[i_var].new_var = size_pair.i_var;
02230                     replace_hash = true;
02231                }
02232                break;
02233                // ---------------------------------------------------
02234                // Binary operator where 
02235                // both operators are variables
02236                case AddvvOp:
02237                case SubvvOp:
02238                if( (tape[arg[0]].connect_type == csum_connected) |
02239                    (tape[arg[1]].connect_type == csum_connected)
02240                )
02241                {
02242                     // convert to a sequence of summation operators
02243                     size_pair = record_csum(
02244                          tape                , // inputs
02245                          i_var               ,
02246                          play->num_par_rec() ,
02247                          play->GetPar()      ,
02248                          rec                 ,
02249                          csum_work
02250                     );
02251                     tape[i_var].new_op  = size_pair.i_op;
02252                     tape[i_var].new_var = size_pair.i_var;
02253                     // abort rest of this case
02254                     break;
02255                }
02256                case DivvvOp:
02257                case MulvvOp:
02258                case PowvvOp:
02259                match_var = binary_match(
02260                     tape                ,  // inputs 
02261                     i_var               ,
02262                     play->num_par_rec() ,
02263                     play->GetPar()      ,
02264                     hash_table_var      ,
02265                     code                  // outputs
02266                );
02267                if( match_var > 0 )
02268                     tape[i_var].new_var = match_var;
02269                else
02270                {    size_pair = record_vv(
02271                          tape                , // inputs
02272                          i_var               ,
02273                          play->num_par_rec() ,
02274                          play->GetPar()      ,
02275                          rec                 ,
02276                          op                  ,
02277                          arg
02278                     );
02279                     tape[i_var].new_op  = size_pair.i_op;
02280                     tape[i_var].new_var = size_pair.i_var;
02281                     replace_hash = true;
02282                }
02283                break;
02284                // ---------------------------------------------------
02285                // Conditional expression operators
02286                case CExpOp:
02287                CPPAD_ASSERT_NARG_NRES(op, 6, 1);
02288                new_arg[0] = arg[0];
02289                new_arg[1] = arg[1];
02290                mask = 1;
02291                for(i = 2; i < 6; i++)
02292                {    if( arg[1] & mask )
02293                     {    new_arg[i] = tape[arg[i]].new_var;
02294                          CPPAD_ASSERT_UNKNOWN( 
02295                               size_t(new_arg[i]) < num_var 
02296                          );
02297                     }
02298                     else new_arg[i] = rec->PutPar( 
02299                               play->GetPar( arg[i] )
02300                     );
02301                     mask = mask << 1;
02302                }
02303                rec->PutArg(
02304                     new_arg[0] ,
02305                     new_arg[1] ,
02306                     new_arg[2] ,
02307                     new_arg[3] ,
02308                     new_arg[4] ,
02309                     new_arg[5] 
02310                );
02311                tape[i_var].new_op  = rec->num_op_rec();
02312                tape[i_var].new_var = rec->PutOp(op);
02313                break;
02314                // ---------------------------------------------------
02315                // Operations with no arguments and no results
02316                case EndOp:
02317                CPPAD_ASSERT_NARG_NRES(op, 0, 0);
02318                rec->PutOp(op);
02319                break;
02320                // ---------------------------------------------------
02321                // Operations with no arguments and one result
02322                case InvOp:
02323                CPPAD_ASSERT_NARG_NRES(op, 0, 1);
02324                tape[i_var].new_op  = rec->num_op_rec();
02325                tape[i_var].new_var = rec->PutOp(op);
02326                break;
02327                // ---------------------------------------------------
02328                // Operations with one argument that is a parameter
02329                case ParOp:
02330                CPPAD_ASSERT_NARG_NRES(op, 1, 1);
02331                new_arg[0] = rec->PutPar( play->GetPar(arg[0] ) );
02332 
02333                rec->PutArg( new_arg[0] );
02334                tape[i_var].new_op  = rec->num_op_rec();
02335                tape[i_var].new_var = rec->PutOp(op);
02336                break;
02337                // ---------------------------------------------------
02338                // Load using a parameter index
02339                case LdpOp:
02340                CPPAD_ASSERT_NARG_NRES(op, 3, 1);
02341                new_arg[0] = new_vecad_ind[ arg[0] ];
02342                new_arg[1] = arg[1];
02343                new_arg[2] = rec->num_load_op_rec();
02344                CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
02345                rec->PutArg( 
02346                     new_arg[0], 
02347                     new_arg[1], 
02348                     new_arg[2]
02349                );
02350                tape[i_var].new_op  = rec->num_op_rec();
02351                tape[i_var].new_var = rec->PutLoadOp(op);
02352                break;
02353                // ---------------------------------------------------
02354                // Load using a variable index
02355                case LdvOp:
02356                CPPAD_ASSERT_NARG_NRES(op, 3, 1);
02357                new_arg[0] = new_vecad_ind[ arg[0] ];
02358                new_arg[1] = tape[arg[1]].new_var;
02359                new_arg[2] = rec->num_load_op_rec();
02360                CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
02361                CPPAD_ASSERT_UNKNOWN( size_t(new_arg[1]) < num_var );
02362                rec->PutArg( 
02363                     new_arg[0], 
02364                     new_arg[1], 
02365                     new_arg[2]
02366                );
02367                tape[i_var].new_var = rec->num_op_rec();
02368                tape[i_var].new_var = rec->PutLoadOp(op);
02369                break;
02370                // ---------------------------------------------------
02371                // Store a parameter using a parameter index
02372                case StppOp:
02373                CPPAD_ASSERT_NARG_NRES(op, 3, 0);
02374                new_arg[0] = new_vecad_ind[ arg[0] ];
02375                new_arg[1] = rec->PutPar( play->GetPar(arg[1]) );
02376                new_arg[2] = rec->PutPar( play->GetPar(arg[2]) );
02377                CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
02378                rec->PutArg(
02379                     new_arg[0], 
02380                     new_arg[1], 
02381                     new_arg[2]
02382                );
02383                rec->PutOp(op);
02384                break;
02385                // ---------------------------------------------------
02386                // Store a parameter using a variable index
02387                case StvpOp:
02388                CPPAD_ASSERT_NARG_NRES(op, 3, 0);
02389                new_arg[0] = new_vecad_ind[ arg[0] ];
02390                new_arg[1] = tape[arg[1]].new_var;
02391                new_arg[2] = rec->PutPar( play->GetPar(arg[2]) );
02392                CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
02393                CPPAD_ASSERT_UNKNOWN( size_t(new_arg[1]) < num_var );
02394                rec->PutArg(
02395                     new_arg[0], 
02396                     new_arg[1], 
02397                     new_arg[2]
02398                );
02399                rec->PutOp(op);
02400                break;
02401                // ---------------------------------------------------
02402                // Store a variable using a parameter index
02403                case StpvOp:
02404                CPPAD_ASSERT_NARG_NRES(op, 3, 0);
02405                new_arg[0] = new_vecad_ind[ arg[0] ];
02406                new_arg[1] = rec->PutPar( play->GetPar(arg[1]) );
02407                new_arg[2] = tape[arg[2]].new_var;
02408                CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
02409                CPPAD_ASSERT_UNKNOWN( size_t(new_arg[1]) < num_var );
02410                CPPAD_ASSERT_UNKNOWN( size_t(new_arg[2]) < num_var );
02411                rec->PutArg(
02412                     new_arg[0], 
02413                     new_arg[1], 
02414                     new_arg[2]
02415                );
02416                rec->PutOp(op);
02417                break;
02418                // ---------------------------------------------------
02419                // Store a variable using a variable index
02420                case StvvOp:
02421                CPPAD_ASSERT_NARG_NRES(op, 3, 0);
02422                new_arg[0] = new_vecad_ind[ arg[0] ];
02423                new_arg[1] = tape[arg[1]].new_var;
02424                new_arg[2] = tape[arg[2]].new_var;
02425                CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
02426                CPPAD_ASSERT_UNKNOWN( size_t(new_arg[1]) < num_var );
02427                CPPAD_ASSERT_UNKNOWN( size_t(new_arg[2]) < num_var );
02428                rec->PutArg(
02429                     new_arg[0], 
02430                     new_arg[1], 
02431                     new_arg[2]
02432                );
02433                rec->PutOp(op);
02434                break;
02435 
02436                // -----------------------------------------------------------
02437                case UserOp:
02438                CPPAD_ASSERT_NARG_NRES(op, 4, 0);
02439                if( user_state == user_start )
02440                {    user_state = user_arg;
02441                     CPPAD_ASSERT_UNKNOWN( user_curr > 0 );
02442                     user_curr--;
02443                     user_info[user_curr].op_begin = rec->num_op_rec();
02444                }
02445                else
02446                {    user_state = user_start;
02447                     user_info[user_curr].op_end = rec->num_op_rec() + 1;
02448                }
02449                // user_index, user_id, user_n, user_m
02450                if( user_info[user_curr].connect_type != not_connected )
02451                {    rec->PutArg(arg[0], arg[1], arg[2], arg[3]);
02452                     rec->PutOp(UserOp);
02453                }
02454                break;
02455 
02456                case UsrapOp:
02457                CPPAD_ASSERT_NARG_NRES(op, 1, 0);
02458                if( user_info[user_curr].connect_type != not_connected )
02459                {    new_arg[0] = rec->PutPar( play->GetPar(arg[0]) );
02460                     rec->PutArg(new_arg[0]);
02461                     rec->PutOp(UsrapOp);
02462                }
02463                break;
02464 
02465                case UsravOp:
02466                CPPAD_ASSERT_NARG_NRES(op, 1, 0);
02467                if( user_info[user_curr].connect_type != not_connected )
02468                {    new_arg[0] = tape[arg[0]].new_var;
02469                     if( size_t(new_arg[0]) < num_var )
02470                     {    rec->PutArg(new_arg[0]);
02471                          rec->PutOp(UsravOp);
02472                     }
02473                     else
02474                     {    // This argument does not affect the result and
02475                          // has been optimized out so use nan in its place.
02476                          new_arg[0] = rec->PutPar( nan(Base(0)) );
02477                          rec->PutArg(new_arg[0]);
02478                          rec->PutOp(UsrapOp);
02479                     }
02480                }
02481                break;
02482 
02483                case UsrrpOp:
02484                CPPAD_ASSERT_NARG_NRES(op, 1, 0);
02485                if( user_info[user_curr].connect_type != not_connected )
02486                {    new_arg[0] = rec->PutPar( play->GetPar(arg[0]) );
02487                     rec->PutArg(new_arg[0]);
02488                     rec->PutOp(UsrrpOp);
02489                }
02490                break;
02491                
02492                case UsrrvOp:
02493                CPPAD_ASSERT_NARG_NRES(op, 0, 1);
02494                if( user_info[user_curr].connect_type != not_connected )
02495                {    tape[i_var].new_op  = rec->num_op_rec();
02496                     tape[i_var].new_var = rec->PutOp(UsrrvOp);
02497                }
02498                break;
02499                // ---------------------------------------------------
02500 
02501                // all cases should be handled above
02502                default:
02503                CPPAD_ASSERT_UNKNOWN(false);
02504 
02505           }
02506           if( replace_hash )
02507           {    // The old variable index i_var corresponds to the 
02508                // new variable index tape[i_var].new_var. In addition
02509                // this is the most recent variable that has this code.
02510                hash_table_var[code] = i_var;
02511           }
02512 
02513      }
02514      // modify the dependent variable vector to new indices
02515      for(i = 0; i < dep_taddr.size(); i++ )
02516      {    CPPAD_ASSERT_UNKNOWN( size_t(tape[dep_taddr[i]].new_var) < num_var );
02517           dep_taddr[i] = tape[ dep_taddr[i] ].new_var;
02518      }
02519 
02520 # ifndef NDEBUG
02521      size_t num_new_op = rec->num_op_rec();
02522      for(i_var = 0; i_var < tape.size(); i_var++)
02523           CPPAD_ASSERT_UNKNOWN( tape[i_var].new_op < num_new_op );
02524 # endif
02525 
02526      // Move skip information from user_info to cskip_info
02527      for(i = 0; i < user_info.size(); i++)
02528      {    if( user_info[i].connect_type == cexp_connected )
02529           {    std::set<class_cexp_pair>::const_iterator itr =
02530                     user_info[i].cexp_set.begin();
02531                while( itr != user_info[i].cexp_set.end() )
02532                {    j = itr->index_;
02533                     k = user_info[i].op_begin;
02534                     while(k < user_info[i].op_end)
02535                     {    if( itr->compare_ == true )
02536                               cskip_info[j].skip_op_false.push_back(k++);
02537                          else cskip_info[j].skip_op_true.push_back(k++);
02538                     }
02539                     itr++;
02540                }
02541           }
02542      }
02543 
02544      // fill in the arguments for the CSkip operations
02545      CPPAD_ASSERT_UNKNOWN( cskip_info_next == cskip_info.size() );
02546      for(i = 0; i < cskip_info.size(); i++)
02547      {    struct_cskip_info info = cskip_info[i];
02548           if( info.i_arg > 0 )
02549           {    CPPAD_ASSERT_UNKNOWN( info.n_op_true==info.skip_op_true.size() );
02550                CPPAD_ASSERT_UNKNOWN(info.n_op_false==info.skip_op_false.size());
02551                size_t n_true  = 
02552                     info.skip_var_true.size() + info.skip_op_true.size();
02553                size_t n_false = 
02554                     info.skip_var_false.size() + info.skip_op_false.size();
02555                size_t i_arg   = info.i_arg;
02556                rec->ReplaceArg(i_arg++, info.cop   );
02557                rec->ReplaceArg(i_arg++, info.flag  );
02558                rec->ReplaceArg(i_arg++, info.left  );
02559                rec->ReplaceArg(i_arg++, info.right );
02560                rec->ReplaceArg(i_arg++, n_true     );
02561                rec->ReplaceArg(i_arg++, n_false    );
02562                for(j = 0; j < info.skip_var_true.size(); j++)
02563                {    i_var = info.skip_var_true[j];
02564                     CPPAD_ASSERT_UNKNOWN( tape[i_var].new_op > 0 );
02565                     rec->ReplaceArg(i_arg++, tape[i_var].new_op );
02566                } 
02567                for(j = 0; j < info.skip_op_true.size(); j++)
02568                {    i_op = info.skip_op_true[j];
02569                     rec->ReplaceArg(i_arg++, i_op);
02570                } 
02571                for(j = 0; j < info.skip_var_false.size(); j++)
02572                {    i_var = info.skip_var_false[j];
02573                     CPPAD_ASSERT_UNKNOWN( tape[i_var].new_op > 0 );
02574                     rec->ReplaceArg(i_arg++, tape[i_var].new_op );
02575                } 
02576                for(j = 0; j < info.skip_op_false.size(); j++)
02577                {    i_op = info.skip_op_false[j];
02578                     rec->ReplaceArg(i_arg++, i_op);
02579                } 
02580                rec->ReplaceArg(i_arg++, n_true + n_false);
02581 # ifndef NDEBUG
02582                size_t n_arg   = 7 + n_true + n_false; 
02583                CPPAD_ASSERT_UNKNOWN( info.i_arg + n_arg == i_arg );
02584 # endif
02585           }
02586      }
02587 }
02588 
02589 } // END_CPPAD_OPTIMIZE_NAMESPACE
02590 
02591 /*!
02592 Optimize a player object operation sequence
02593 
02594 The operation sequence for this object is replaced by one with fewer operations
02595 but the same funcition and derivative values.
02596 
02597 \tparam Base
02598 base type for the operator; i.e., this operation was recorded
02599 using AD< \a Base > and computations by this routine are done using type 
02600 \a Base.
02601 
02602 */
02603 template <class Base>
02604 void ADFun<Base>::optimize(void)
02605 {    // place to store the optimized version of the recording
02606      recorder<Base> rec;
02607 
02608      // number of independent variables
02609      size_t n = ind_taddr_.size();
02610 
02611 # ifndef NDEBUG
02612      size_t i, j, m = dep_taddr_.size();
02613      CppAD::vector<Base> x(n), y(m), check(m);
02614      Base max_taylor(0);
02615      bool check_zero_order = num_order_taylor_ > 0;
02616      if( check_zero_order )
02617      {    // zero order coefficients for independent vars
02618           for(j = 0; j < n; j++)
02619           {    CPPAD_ASSERT_UNKNOWN( play_.GetOp(j+1) == InvOp );
02620                CPPAD_ASSERT_UNKNOWN( ind_taddr_[j]    == j+1   );
02621                x[j] = taylor_[ ind_taddr_[j] * cap_order_taylor_ + 0];
02622           }
02623           // zero order coefficients for dependent vars
02624           for(i = 0; i < m; i++)
02625           {    CPPAD_ASSERT_UNKNOWN( dep_taddr_[i] < num_var_tape_  );
02626                y[i] = taylor_[ dep_taddr_[i] * cap_order_taylor_ + 0];
02627           }
02628           // maximum zero order coefficient not counting BeginOp at beginning
02629           // (which is correpsonds to uninitialized memory).
02630           for(i = 1; i < num_var_tape_; i++)
02631           {    if(  abs_geq(taylor_[i*cap_order_taylor_+0] , max_taylor) )
02632                     max_taylor = taylor_[i*cap_order_taylor_+0];
02633           }
02634      }
02635 # endif
02636 
02637      // create the optimized recording
02638      CppAD::optimize::optimize_run<Base>(n, dep_taddr_, &play_, &rec);
02639 
02640      // number of variables in the recording
02641      num_var_tape_  = rec.num_var_rec();
02642 
02643      // now replace the recording
02644      play_.get(rec);
02645 
02646      // free memory allocated for sparse Jacobian calculation
02647      // (the results are no longer valid)
02648      for_jac_sparse_pack_.resize(0, 0);
02649      for_jac_sparse_set_.resize(0,0);
02650 
02651      // free old Taylor coefficient memory
02652      taylor_.free();
02653      num_order_taylor_     = 0;
02654      cap_order_taylor_     = 0;
02655 
02656      // resize and initilaize conditional skip vector
02657      // (must use player size because it now has the recoreder information)
02658      cskip_op_.erase();
02659      cskip_op_.extend( play_.num_op_rec() );
02660 
02661 # ifndef NDEBUG
02662      if( check_zero_order )
02663      {
02664           // zero order forward calculation using new operation sequence
02665           check = Forward(0, x);
02666 
02667           // check results
02668           Base eps = 10. * epsilon<Base>();
02669           for(i = 0; i < m; i++) CPPAD_ASSERT_KNOWN( 
02670                abs_geq( eps * max_taylor , check[i] - y[i] ) ,
02671                "Error during check of f.optimize()."
02672           );
02673 
02674           // Erase memory that this calculation was done so NDEBUG gives 
02675           // same final state for this object (from users perspective)
02676           num_order_taylor_     = 0;
02677      }
02678 # endif
02679 }
02680 
02681 } // END_CPPAD_NAMESPACE
02682 # endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines