CppAD: A C++ Algorithmic Differentiation Package 20110419
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-11 Bradley M. Bell
00007 
00008 CppAD is distributed under multiple licenses. This distribution is under
00009 the terms of the 
00010                     Common 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         Taylor
00020         var
00021         CppAD
00022         cppad
00023 $$
00024 
00025 $section Optimize the Tape Corresponding to an ADFun object$$
00026 
00027 $index optimize$$
00028 $index tape, optimize$$
00029 $index sequence, optimize operations$$
00030 $index operations, optimize sequence$$
00031 $index speed, optimize$$
00032 $index memory, optimize$$
00033 
00034 $head Syntax$$
00035 $icode%f%.optimize()%$$
00036 
00037 
00038 $head Purpose$$
00039 The operation sequence corresponding to an $cref/ADFun/$$ object can
00040 be very large and involve many operations; see the
00041 size functions in $cref/seq_property/$$.
00042 This enables one to reduce the number of operations
00043 and thereby reduce the time and the memory required to
00044 compute function and derivative values using an $code ADFun$$ object.
00045 
00046 $head f$$
00047 The object $icode f$$ has prototype
00048 $icode%
00049         ADFun<%Base%> %f%
00050 %$$
00051 
00052 $head Improvements$$
00053 You can see the reduction in number of variables in the operation sequence
00054 by calling the function $cref/f.size_var()/seq_property/size_var/$$
00055 before and after the optimization procedure.
00056 Given that the optimization procedure takes time,
00057 it may be faster to skip this optimize procedure and just compute
00058 derivatives using the original operation sequence.
00059 
00060 $subhead Testing$$
00061 You can run the CppAD $cref/speed/speed_main/$$ tests and see
00062 the corresponding changes in number of variables and execution time.
00063 After the $cref/make test/InstallUnix/make test/$$ command
00064 is run in the $code speed/cppad$$ directory,
00065 you can run the speed tests with optimization 
00066 using the command
00067 $codep
00068         speed/cppad/cppad speed 123 optimize
00069 $$ 
00070 and without optimization using the command
00071 $codep
00072         speed/cppad/cppad speed 123 
00073 $$
00074 Note that $code 123$$ is used for a random number seed and can be
00075 replaced by any integer in the commands above.
00076 
00077 $head Efficiency$$
00078 The $code optimize$$ member function
00079 may greatly reduce the number of variables 
00080 in the operation sequence; see $cref/size_var/seq_property/size_var/$$.
00081 If a $cref/zero order forward/ForwardZero/$$ calculation is done during
00082 the construction of f, it will require more memory
00083 and time that after the optimization procedure.
00084 In addition, it will need to be redone.
00085 For this reason, it is more efficient to use 
00086 $codei%
00087         ADFun<%Base%> %f%;
00088         %f%.Dependent(%x%, %y%);
00089         %f%.optimize();
00090 %$$
00091 instead of
00092 $codei%
00093         ADFun<%Base%> %f%(%x%, %y%)
00094         %f%.optimize();
00095 %$$ 
00096 See the discussion about
00097 $cref/sequence constructors/FunConstruct/Sequence Constructor/$$.
00098 
00099 $head Comparison Operators$$
00100 Any comparison operators that are in the tape are removed by this operation.
00101 Hence the return value of $cref/CompareChange/$$ will always be zero
00102 for an optimized tape (even if $code NDEBUG$$ is not defined).
00103 
00104 $head Checking Optimization$$
00105 $index NDEBUG$$
00106 If $code NDEBUG$$ is not defined
00107 and $cref/f.size_taylor()/size_taylor/$$ is greater than zero,
00108 a $cref/ForwardZero/$$ calculation is done using the optimized version
00109 of $icode f$$ and the results are checked to see that they are
00110 the same as before.
00111 If they are not the same, the
00112 $cref/ErrorHandler/$$ is called with a known error message
00113 related to $icode%f%.optimize()%$$.
00114 
00115 $head Example$$
00116 $children%
00117         example/optimize.cpp
00118 %$$
00119 The file
00120 $xref/optimize.cpp/$$
00121 contains an example and test of this operation.
00122 It returns true if it succeeds and false otherwise.
00123 
00124 $end
00125 -----------------------------------------------------------------------------
00126 */
00127 # include <stack>
00128 
00129 CPPAD_BEGIN_NAMESPACE
00130 /*!
00131 \file optimize.hpp
00132 Routines for optimizing a tape
00133 */
00134 
00135 
00136 /*!
00137 State for this variable set during reverse sweep.
00138 */
00139 enum optimize_connection {
00140         /// There is no operation that connects this variable to the
00141         /// independent variables.
00142         not_connected     ,
00143 
00144         /// There is one or more operations that connects this variable to the
00145         /// independent variables.
00146         yes_connected      ,
00147 
00148         /// There is only one parrent operation that connects this variable to 
00149         /// the independent variables and it is one of the following:
00150         ///  AddvvOp, AddpvOp, SubpvOp, SubvpOp, or SubvvOp.
00151         sum_connected      ,
00152 
00153         /// Satisfies sum_connected and in addition 
00154         /// the parrent operation is one of the following:
00155         ///  AddvvOp, AddpvOp, SubpvOp, SubvpOp, or SubvvOp.
00156         csum_connected
00157 
00158 };
00159 
00160 
00161 /*!
00162 Structure used by \c optimize to hold information about one variable.
00163 in the old operation seqeunce.
00164 */
00165 struct optimize_old_variable {
00166         /// Operator for which this variable is the result, \c NumRes(op) > 0.
00167         /// Set by the reverse sweep at beginning of optimization.
00168         OpCode              op;       
00169 
00170         /// Pointer to first argument (child) for this operator.
00171         /// Set by the reverse sweep at beginning of optimization.
00172         const size_t*       arg;
00173 
00174         /// How is this variable connected to the independent variables
00175         optimize_connection connect; 
00176 
00177         /// Set during forward sweep to the index in the
00178         /// new operation sequence corresponding to this old varable.
00179         size_t new_var;
00180 };
00181 
00182 /*!
00183 Structures used by \c optimize_record_csum
00184 to hold information about one variable.
00185 */
00186 struct optimize_csum_variable {
00187         /// Operator for which this variable is the result, \c NumRes(op) > 0.
00188         OpCode              op;       
00189 
00190         /// Pointer to first argument (child) for this operator.
00191         /// Set by the reverse sweep at beginning of optimization.
00192         const size_t*       arg;
00193 
00194         /// Is this variable added to the summation
00195         /// (if not it is subtracted)
00196         bool                add;
00197 };
00198 
00199 /*!
00200 Structure used to pass work space from \c optimize to \c optimize_record_csum
00201 (so that stacks do not start from zero size every time).
00202 */
00203 struct optimize_csum_stacks {
00204         /// stack of operations in the cummulative summation 
00205         std::stack<struct optimize_csum_variable>   op_stack;
00206         /// stack of variables to be added
00207         std::stack<size_t >                         add_stack;
00208         /// stack of variables to be subtracted
00209         std::stack<size_t >                         sub_stack;
00210 };
00211 
00212 
00213 /*!
00214 Documents arguments that are common to optimization helper functions
00215 (should not be called).
00216 
00217 \param tape
00218 is a vector that maps a variable index, in the old operation sequence,
00219 to an <tt>optimize_old_variable</tt> information record.
00220 Note that the index for this vector must be greater than or equal zero and 
00221 less than <tt>tape.size()</tt>.
00222 
00223 \li <tt>tape[i].op</tt> 
00224 is the operator in the old operation sequence
00225 corresponding to the old variable index \c i.
00226 Assertion: <tt>NumRes(tape[i].op) > 0</tt>.
00227 
00228 \li <tt>tape[i].arg</tt> 
00229 for <tt>j < NumArg( tape[i].op ), tape[i].arg[j]</tt>
00230 is the j-th the argument, in the old operation sequence,
00231 corresponding to the old variable index \c i.
00232 Assertion: <tt>tape[i].arg[j] < i</tt>.
00233 
00234 \li <tt>tape[i].new_var</tt>
00235 Suppose 
00236 <tt>i <= current, j < NumArg( tape[i].op ), and k = tape[i].arg[j]</tt>,
00237 and \c j corresponds to a varialbe for operator <tt>tape[i].op</tt>.
00238 It follows that <tt>tape[k].new_var</tt>
00239 has alread been set to the variable in the new operation sequence 
00240 corresponding to the old variable index \c k.
00241 This means that the \c new_var value has been set
00242 for all the possible arguments that come before \a current.
00243 
00244 \param current
00245 is the index in the old operation sequence for 
00246 the variable corresponding to the result for the current operator.
00247 Assertions: 
00248 <tt>
00249 current < tape.size(),
00250 NumRes( tape[current].op ) > 0.
00251 </tt>
00252 
00253 \param npar
00254 is the number of paraemters corresponding to this operation sequence.
00255 
00256 \param par
00257 is a vector of length \a npar containing the parameters
00258 for this operation sequence; i.e.,
00259 given a parameter index \c i, the corresponding parameter value is
00260 <tt>par[i]</tt>.
00261 */
00262 template <class Base>
00263 void optimize_prototype(
00264         const CppAD::vector<struct optimize_old_variable>& tape           ,
00265         size_t                                             current        ,
00266         size_t                                             npar           ,
00267         const Base*                                        par            )
00268 {       CPPAD_ASSERT_UNKNOWN(false); }
00269 
00270 /*!
00271 Check a unary operator for a complete match with a previous operator.
00272 
00273 A complete match means that the result of the previous operator
00274 can be used inplace of the result for current operator.
00275 
00276 \copydetails optimize_prototype
00277 
00278 \param hash_table_var
00279 is a vector with size CPPAD_HASH_TABLE_SIZE
00280 that maps a hash code to the corresponding 
00281 variable index in the old operation sequence.
00282 All the values in this table must be less than \a current.
00283 
00284 \param code
00285 The input value of code does not matter.
00286 The output value of code is the hash code corresponding to
00287 this operation in the new operation sequence.
00288 
00289 \return
00290 If the return value is zero,
00291 no match was found.
00292 If the return value is greater than zero,
00293 it is the index of a new variable that can be used to replace the 
00294 old variable.
00295 
00296 \par Restrictions:
00297 NumArg( tape[current].op ) == 1
00298 */
00299 template <class Base>
00300 size_t optimize_unary_match(
00301         const CppAD::vector<struct optimize_old_variable>& tape           ,
00302         size_t                                             current        ,
00303         size_t                                             npar           ,
00304         const Base*                                        par            ,
00305         const CppAD::vector<size_t>&                       hash_table_var ,
00306         unsigned short&                                    code           )
00307 {       const size_t *arg = tape[current].arg;
00308         OpCode        op  = tape[current].op;
00309         size_t new_arg[1];
00310         
00311         CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
00312         CPPAD_ASSERT_UNKNOWN( NumRes(op) > 0  );
00313         CPPAD_ASSERT_UNKNOWN( arg[0] < current );
00314         new_arg[0] = tape[arg[0]].new_var;
00315         CPPAD_ASSERT_UNKNOWN( new_arg[0] < current );
00316         code = hash_code(
00317                 op                  , 
00318                 new_arg             ,
00319                 npar                ,
00320                 par
00321         );
00322         size_t  i               = hash_table_var[code];
00323         CPPAD_ASSERT_UNKNOWN( i < current );
00324         if( op == tape[i].op )
00325         {       size_t k = tape[i].arg[0];
00326                 CPPAD_ASSERT_UNKNOWN( k < i );
00327                 if (new_arg[0] == tape[k].new_var )
00328                         return tape[i].new_var;
00329         }
00330         return 0;
00331 } 
00332 
00333 /*!
00334 Check a binary operator for a complete match with a previous operator,
00335 
00336 \copydetails optimize_unary_match
00337 
00338 \par Restrictions:
00339 The binary operator must be an addition, subtraction, multiplication, division
00340 or power operator.
00341 */
00342 template <class Base>
00343 inline size_t optimize_binary_match(
00344         const CppAD::vector<struct optimize_old_variable>& tape           ,
00345         size_t                                             current        ,
00346         size_t                                             npar           ,
00347         const Base*                                        par            ,
00348         const CppAD::vector<size_t>&                       hash_table_var ,
00349         unsigned short&                                    code           )
00350 {       OpCode        op         = tape[current].op;
00351         const size_t* arg        = tape[current].arg;
00352         size_t        new_arg[2];
00353         bool          parameter[2];
00354 
00355         // initialize return value
00356         size_t  match_var = 0;
00357 
00358         CPPAD_ASSERT_UNKNOWN( NumArg(op) == 2 );
00359         CPPAD_ASSERT_UNKNOWN( NumRes(op) >  0 );
00360         switch(op)
00361         {       // parameter op variable ----------------------------------
00362                 case AddpvOp:
00363                 case MulpvOp:
00364                 case DivpvOp:
00365                 case PowpvOp:
00366                 case SubpvOp:
00367                 // arg[0]
00368                 parameter[0] = true;
00369                 new_arg[0]   = arg[0];
00370                 CPPAD_ASSERT_UNKNOWN( arg[0] < npar );
00371                 // arg[1]
00372                 parameter[1] = false;
00373                 new_arg[1]   = tape[arg[1]].new_var;
00374                 CPPAD_ASSERT_UNKNOWN( arg[1] < current );
00375                 break;
00376 
00377                 // variable op parameter -----------------------------------
00378                 case DivvpOp:
00379                 case PowvpOp:
00380                 case SubvpOp:
00381                 // arg[0]
00382                 parameter[0] = false;
00383                 new_arg[0]   = tape[arg[0]].new_var;
00384                 CPPAD_ASSERT_UNKNOWN( arg[0] < current );
00385                 // arg[1]
00386                 parameter[1] = true;
00387                 new_arg[1]   = arg[1];
00388                 CPPAD_ASSERT_UNKNOWN( arg[1] < npar );
00389                 break;
00390 
00391                 // variable op variable -----------------------------------
00392                 case AddvvOp:
00393                 case MulvvOp:
00394                 case DivvvOp:
00395                 case PowvvOp:
00396                 case SubvvOp:
00397                 // arg[0]
00398                 parameter[0] = false;
00399                 new_arg[0]   = tape[arg[0]].new_var;
00400                 CPPAD_ASSERT_UNKNOWN( arg[0] < current );
00401                 // arg[1]
00402                 parameter[1] = false;
00403                 new_arg[1]   = tape[arg[1]].new_var;
00404                 CPPAD_ASSERT_UNKNOWN( arg[1] < current );
00405                 break;
00406 
00407                 // must be one of the cases above
00408                 default:
00409                 CPPAD_ASSERT_UNKNOWN(false);
00410         }
00411         code = hash_code(
00412                 op                  , 
00413                 new_arg             ,
00414                 npar                ,
00415                 par
00416         );
00417         size_t  i  = hash_table_var[code];
00418         CPPAD_ASSERT_UNKNOWN( i < current );
00419         if( op == tape[i].op )
00420         {       bool match = true;
00421                 size_t j;
00422                 for(j = 0; j < 2; j++)
00423                 {       size_t k = tape[i].arg[j];
00424                         if( parameter[j] )
00425                         {       CPPAD_ASSERT_UNKNOWN( k < npar );
00426                                 match &= IdenticalEqualPar(
00427                                         par[ arg[j] ], par[k]
00428                                 );
00429                         }
00430                         else
00431                         {       CPPAD_ASSERT_UNKNOWN( k < i );
00432                                 match &= (new_arg[j] == tape[k].new_var);
00433                         }
00434                 }
00435                 if( match )
00436                         match_var = tape[i].new_var;
00437         }
00438         if( (match_var > 0) | ( (op != AddvvOp) & (op != MulvvOp ) ) )
00439                 return match_var;
00440 
00441         // check for match with argument order switched ----------------------
00442         CPPAD_ASSERT_UNKNOWN( op == AddvvOp || op == MulvvOp );
00443         i          = new_arg[0];
00444         new_arg[0] = new_arg[1];
00445         new_arg[1] = i;
00446         unsigned short code_switch = hash_code(
00447                 op                  , 
00448                 new_arg             ,
00449                 npar                ,
00450                 par
00451         );
00452         i  = hash_table_var[code_switch];
00453         CPPAD_ASSERT_UNKNOWN( i < current );
00454         if( op == tape[i].op )
00455         {       bool match = true;
00456                 size_t j;
00457                 for(j = 0; j < 2; j++)
00458                 {       size_t k = tape[i].arg[j];
00459                         CPPAD_ASSERT_UNKNOWN( k < i );
00460                         match &= (new_arg[j] == tape[k].new_var);
00461                 }
00462                 if( match )
00463                         match_var = tape[i].new_var;
00464         }
00465         return match_var;
00466 } 
00467 
00468 /*!
00469 Record an operation of the form (parameter op variable).
00470 
00471 \copydetails optimize_prototype
00472 
00473 \param rec
00474 is the object that will record the operations.
00475 
00476 \param op
00477 is the operator that we are recording which must be one of the following:
00478 AddpvOp, DivpvOp, MulpvOp, PowvpOp, SubpvOp.
00479  
00480 \param arg
00481 is the vector of arguments for this operator.
00482 
00483 \return
00484 the result value is the index corresponding to the current
00485 operation in the new operation sequence.
00486 */
00487 template <class Base>
00488 size_t optimize_record_pv(
00489         const CppAD::vector<struct optimize_old_variable>& tape           ,
00490         size_t                                             current        ,
00491         size_t                                             npar           ,
00492         const Base*                                        par            ,
00493         recorder<Base>*                                    rec            ,
00494         OpCode                                             op             ,
00495         const size_t*                                      arg            )
00496 {
00497 # ifndef NDEBUG
00498         switch(op)
00499         {       case AddpvOp:
00500                 case DivpvOp:
00501                 case MulpvOp:
00502                 case PowpvOp:
00503                 case SubpvOp:
00504                 break;
00505 
00506                 default:
00507                 CPPAD_ASSERT_UNKNOWN(false);
00508         }
00509 # endif
00510         CPPAD_ASSERT_UNKNOWN( arg[0] < npar    );
00511         CPPAD_ASSERT_UNKNOWN( arg[1] < current );
00512         size_t new_arg[2];
00513         new_arg[0]   = rec->PutPar( par[arg[0]] );
00514         new_arg[1]   = tape[ arg[1] ].new_var;
00515         rec->PutArg( new_arg[0], new_arg[1] );
00516         size_t i     = rec->PutOp(op);
00517         CPPAD_ASSERT_UNKNOWN( new_arg[0] < i );
00518         return i;
00519 }
00520 
00521 
00522 /*!
00523 Record an operation of the form (variable op parameter).
00524 
00525 \copydetails optimize_prototype
00526 
00527 \param rec
00528 is the object that will record the operations.
00529 
00530 \param op
00531 is the operator that we are recording which must be one of the following:
00532 DivvpOp, PowvpOp, SubvpOp.
00533  
00534 \param arg
00535 is the vector of arguments for this operator.
00536 
00537 \return
00538 the result value is the index corresponding to the current
00539 operation in the new operation sequence.
00540 */
00541 template <class Base>
00542 size_t optimize_record_vp(
00543         const CppAD::vector<struct optimize_old_variable>& tape           ,
00544         size_t                                             current        ,
00545         size_t                                             npar           ,
00546         const Base*                                        par            ,
00547         recorder<Base>*                                    rec            ,
00548         OpCode                                             op             ,
00549         const size_t*                                      arg            )
00550 {
00551 # ifndef NDEBUG
00552         switch(op)
00553         {       case DivvpOp:
00554                 case PowvpOp:
00555                 case SubvpOp:
00556                 break;
00557 
00558                 default:
00559                 CPPAD_ASSERT_UNKNOWN(false);
00560         }
00561 # endif
00562         CPPAD_ASSERT_UNKNOWN( arg[0] < current );
00563         CPPAD_ASSERT_UNKNOWN( arg[1] < npar    );
00564         size_t new_arg[2];
00565         new_arg[0]   = tape[ arg[0] ].new_var;
00566         new_arg[1]   = rec->PutPar( par[arg[1]] );
00567         rec->PutArg( new_arg[0], new_arg[1] );
00568         size_t i     = rec->PutOp(op);
00569         CPPAD_ASSERT_UNKNOWN( new_arg[0] < i );
00570         return i;
00571 }
00572 
00573 /*!
00574 Record an operation of the form (variable op variable).
00575 
00576 \copydetails optimize_prototype
00577 
00578 \param rec
00579 is the object that will record the operations.
00580 
00581 \param op
00582 is the operator that we are recording which must be one of the following:
00583 AddvvOp, DivvvOp, MulvvOp, PowvpOp, SubvvOp.
00584  
00585 \param arg
00586 is the vector of arguments for this operator.
00587 
00588 \return
00589 the result value is the index corresponding to the current
00590 operation in the new operation sequence.
00591 */
00592 template <class Base>
00593 size_t optimize_record_vv(
00594         const CppAD::vector<struct optimize_old_variable>& tape           ,
00595         size_t                                             current        ,
00596         size_t                                             npar           ,
00597         const Base*                                        par            ,
00598         recorder<Base>*                                    rec            ,
00599         OpCode                                             op             ,
00600         const size_t*                                      arg            )
00601 {
00602 # ifndef NDEBUG
00603         switch(op)
00604         {       case AddvvOp:
00605                 case DivvvOp:
00606                 case MulvvOp:
00607                 case PowvvOp:
00608                 case SubvvOp:
00609                 break;
00610 
00611                 default:
00612                 CPPAD_ASSERT_UNKNOWN(false);
00613         }
00614 # endif
00615         CPPAD_ASSERT_UNKNOWN( arg[0] < current );
00616         CPPAD_ASSERT_UNKNOWN( arg[1] < current );
00617         size_t new_arg[2];
00618         new_arg[0]   = tape[ arg[0] ].new_var;
00619         new_arg[1]   = tape[ arg[1] ].new_var;
00620         rec->PutArg( new_arg[0], new_arg[1] );
00621         size_t i     = rec->PutOp(op);
00622         CPPAD_ASSERT_UNKNOWN( new_arg[0] < i );
00623         return i;
00624 }
00625 
00626 // ==========================================================================
00627 
00628 /*!
00629 Recording a cummulative cummulative summation starting at its highest parrent.
00630 
00631 \copydetails optimize_prototype
00632 
00633 \param rec
00634 is the object that will record the operations.
00635 
00636 \param work
00637 Is used for computaiton. On input and output,
00638 <tt>work.op_stack.empty()</tt>,
00639 <tt>work.add_stack.empty()</tt>, and
00640 <tt>work.sub_stack.empty()</tt>,
00641 are all true true.
00642 These stacks are passed in so that elements can be allocated once
00643 and then the elements can be reused with calls to \c optimize_record_csum.
00644 
00645 \par Exception
00646 <tt>tape[i].new_var</tt>
00647 is not yet defined for any node \c i that is \c csum_connected
00648 to the \a current node.
00649 For example; suppose that index \c j corresponds to a variable
00650 in the current operator,
00651 <tt>i = tape[current].arg[j]</tt>,
00652 and 
00653 <tt>tape[arg[j]].connect == csum_connected</tt>.
00654 It then follows that
00655 <tt>tape[i].new_var == tape.size()</tt>.
00656 
00657 \par Restriction:
00658 \li <tt>tape[current].op</tt> 
00659 must be one of <tt>AddpvOp, AddvvOp, SubpvOp, SubvpOp, SubvvOp</tt>.
00660 
00661 \li <tt>tape[current].connect</tt> must be \c yes_connected.
00662 */
00663 
00664 
00665 template <class Base>
00666 size_t optimize_record_csum(
00667         const CppAD::vector<struct optimize_old_variable>& tape           ,
00668         size_t                                             current        ,
00669         size_t                                             npar           ,
00670         const Base*                                        par            ,
00671         recorder<Base>*                                    rec            ,
00672         optimize_csum_stacks&                              work           )
00673 {
00674         
00675         CPPAD_ASSERT_UNKNOWN( work.op_stack.empty() );
00676         CPPAD_ASSERT_UNKNOWN( work.add_stack.empty() );
00677         CPPAD_ASSERT_UNKNOWN( work.sub_stack.empty() );
00678         CPPAD_ASSERT_UNKNOWN( tape[current].connect == yes_connected );
00679 
00680         size_t                        i;
00681         OpCode                        op;
00682         const size_t*                 arg;
00683         bool                          add;
00684         struct optimize_csum_variable var;
00685 
00686         var.op  = tape[current].op;
00687         var.arg = tape[current].arg;
00688         var.add = true; 
00689         work.op_stack.push( var );
00690         Base sum_par(0);
00691         while( ! work.op_stack.empty() )
00692         {       var     = work.op_stack.top();
00693                 work.op_stack.pop();
00694                 op      = var.op;
00695                 arg     = var.arg;
00696                 add     = var.add;
00697                 // process first argument to this operator
00698                 switch(op)
00699                 {       case AddpvOp:
00700                         case SubpvOp:
00701                         CPPAD_ASSERT_UNKNOWN( arg[0] < npar );
00702                         if( add )
00703                                 sum_par += par[arg[0]];
00704                         else    sum_par -= par[arg[0]];
00705                         break;
00706 
00707                         case AddvvOp:
00708                         case SubvpOp:
00709                         case SubvvOp:
00710                         if( tape[arg[0]].connect == csum_connected )
00711                         {       CPPAD_ASSERT_UNKNOWN(
00712                                         tape[arg[0]].new_var == tape.size()
00713                                 );
00714                                 var.op  = tape[arg[0]].op;
00715                                 var.arg = tape[arg[0]].arg;
00716                                 var.add = add; 
00717                                 work.op_stack.push( var );
00718                         }
00719                         else if( add )
00720                                 work.add_stack.push(arg[0]);
00721                         else    work.sub_stack.push(arg[0]);
00722                         break;
00723 
00724                         default:
00725                         CPPAD_ASSERT_UNKNOWN(false);
00726                 }
00727                 // process second argument to this operator
00728                 switch(op)
00729                 {
00730                         case SubvpOp:
00731                         CPPAD_ASSERT_UNKNOWN( arg[1] < npar );
00732                         if( add )
00733                                 sum_par -= par[arg[1]];
00734                         else    sum_par += par[arg[1]];
00735                         break;
00736 
00737                         case SubvvOp:
00738                         case SubpvOp:
00739                         add = ! add;
00740 
00741                         case AddvvOp:
00742                         case AddpvOp:
00743                         if( tape[arg[1]].connect == csum_connected )
00744                         {       CPPAD_ASSERT_UNKNOWN(
00745                                         tape[arg[1]].new_var == tape.size()
00746                                 );
00747                                 var.op   = tape[arg[1]].op;
00748                                 var.arg  = tape[arg[1]].arg;
00749                                 var.add  = add;
00750                                 work.op_stack.push( var );
00751                         }
00752                         else if( add )
00753                                 work.add_stack.push(arg[1]);
00754                         else    work.sub_stack.push(arg[1]);
00755                         break;
00756 
00757                         default:
00758                         CPPAD_ASSERT_UNKNOWN(false);
00759                 }
00760         }
00761         // number of variables in this cummulative sum operator
00762         size_t n_add = work.add_stack.size();
00763         size_t n_sub = work.sub_stack.size();
00764         size_t old_arg, new_arg;
00765         rec->PutArg(n_add);                // arg[0]
00766         rec->PutArg(n_sub);                // arg[1]
00767         new_arg = rec->PutPar( sum_par );
00768         rec->PutArg(new_arg);              // arg[2]
00769         for(i = 0; i < n_add; i++)
00770         {       CPPAD_ASSERT_UNKNOWN( ! work.add_stack.empty() );
00771                 old_arg = work.add_stack.top();
00772                 new_arg = tape[old_arg].new_var;
00773                 CPPAD_ASSERT_UNKNOWN( new_arg < tape.size() );
00774                 rec->PutArg(new_arg);      // arg[3+i]
00775                 work.add_stack.pop();
00776         }
00777         for(i = 0; i < n_sub; i++)
00778         {       CPPAD_ASSERT_UNKNOWN( ! work.sub_stack.empty() );
00779                 old_arg = work.sub_stack.top();
00780                 new_arg = tape[old_arg].new_var;
00781                 CPPAD_ASSERT_UNKNOWN( new_arg < tape.size() );
00782                 rec->PutArg(new_arg);      // arg[3 + arg[0] + i]
00783                 work.sub_stack.pop();
00784         }
00785         rec->PutArg(n_add + n_sub);        // arg[3 + arg[0] + arg[1]]
00786         i = rec->PutOp(CSumOp);
00787         CPPAD_ASSERT_UNKNOWN(new_arg < tape.size());
00788 
00789         return i;
00790 }
00791 // ==========================================================================
00792 /*!
00793 Convert a player object to an optimized recorder object
00794 
00795 \tparam Base
00796 base type for the operator; i.e., this operation was recorded
00797 using AD< \a Base > and computations by this routine are done using type 
00798 \a Base.
00799 
00800 \param n
00801 is the number of independent variables on the tape.
00802 
00803 \param dep_taddr
00804 On input this vector contains the indices for each of the dependent
00805 variable values in the operation sequence corresponding to \a play.
00806 Upon return it contains the indices for the same variables but in
00807 the operation sequence corresponding to \a rec.
00808 
00809 \param play
00810 This is the operation sequence that we are optimizing.
00811 It is essentially const, except for play back state which
00812 changes while it plays back the operation seqeunce.
00813 
00814 \param rec
00815 The input contents of this recording does not matter.
00816 Upon return, it contains an optimized verison of the
00817 operation sequence corresponding to \a play.
00818 */
00819 
00820 template <class Base>
00821 void optimize(
00822         size_t                       n         ,
00823         CppAD::vector<size_t>&       dep_taddr ,
00824         player<Base>*                play      ,
00825         recorder<Base>*              rec       ) 
00826 {
00827         // temporary indices
00828         size_t i, j, k;
00829 
00830         // temporary variables
00831         OpCode        op;   // current operator
00832         const size_t *arg;  // operator arguments
00833         size_t        i_var;  // index of first result for current operator
00834 
00835         // range and domain dimensions for F
00836         size_t m = dep_taddr.size();
00837 
00838         // number of variables in the player
00839         const size_t num_var = play->num_rec_var(); 
00840 
00841 # ifndef NDEBUG
00842         // number of paraemters in the player
00843         const size_t num_par = play->num_rec_par();
00844 # endif
00845 
00846         // number of  VecAD indices 
00847         size_t num_vecad_ind   = play->num_rec_vecad_ind();
00848 
00849         // number of VecAD vectors
00850         size_t num_vecad_vec   = play->num_rec_vecad_vec();
00851 
00852         // -------------------------------------------------------------
00853         // data structure that maps variable index in original operation
00854         // sequence to corresponding operator information
00855         CppAD::vector<struct optimize_old_variable> tape(num_var);
00856         // -------------------------------------------------------------
00857         // Determine how each variable is connected to the dependent variables
00858 
00859         // initialize all variables has having no connections
00860         for(i = 0; i < num_var; i++)
00861                 tape[i].connect = not_connected;
00862 
00863         for(j = 0; j < m; j++)
00864         {       // mark dependent variables as having one or more connections
00865                 tape[ dep_taddr[j] ].connect = yes_connected;
00866         }
00867 
00868         // vecad_connect contains a value for each VecAD object.
00869         // vecad maps a VecAD index (which corresponds to the beginning of the
00870         // VecAD object) to the vecad_connect falg for the VecAD object.
00871         CppAD::vector<optimize_connection>   vecad_connect(num_vecad_vec);
00872         CppAD::vector<size_t> vecad(num_vecad_ind);
00873         j = 0;
00874         for(i = 0; i < num_vecad_vec; i++)
00875         {       vecad_connect[i] = not_connected;
00876                 // length of this VecAD
00877                 size_t length = play->GetVecInd(j);
00878                 // set to proper index for this VecAD
00879                 vecad[j] = i; 
00880                 for(k = 1; k <= length; k++)
00881                         vecad[j+k] = num_vecad_vec; // invalid index
00882                 // start of next VecAD
00883                 j       += length + 1;
00884         }
00885         CPPAD_ASSERT_UNKNOWN( j == num_vecad_ind );
00886 
00887         // work space used by UserOp.
00888         typedef std::set<size_t> size_set;
00889         size_t user_q     = 0;       // maximum set element plus one
00890         vector< size_set > user_r;   // sparsity pattern for the argument x
00891         vector< size_set > user_s;   // sparisty pattern for the result y
00892         size_t user_index = 0;       // indentifier for this user_atomic operation
00893         size_t user_id    = 0;       // user identifier for this call to operator
00894         size_t user_i     = 0;       // index in result vector
00895         size_t user_j     = 0;       // index in argument vector
00896         size_t user_m     = 0;       // size of result vector
00897         size_t user_n     = 0;       // size of arugment vector
00898         // next expected operator in a UserOp sequence
00899         enum { user_start, user_arg, user_ret, user_end } user_state;
00900         std::stack<bool> user_keep;
00901 
00902         // Initialize a reverse mode sweep through the operation sequence
00903         size_t i_op;
00904         play->start_reverse(op, arg, i_op, i_var);
00905         CPPAD_ASSERT_UNKNOWN( op == EndOp );
00906         size_t mask;
00907         user_state = user_end;
00908         while(op != BeginOp)
00909         {       // next op
00910                 play->next_reverse(op, arg, i_op, i_var);
00911                 // This if is not necessary becasue last assignment
00912                 // with this value of i_var will have NumRes(op) > 0
00913                 if( NumRes(op) > 0 )
00914                 {       tape[i_var].op = op;
00915                         tape[i_var].arg = arg;
00916                 }
00917 # ifndef NDEBUG
00918                 if( i_op <= n )
00919                 {       CPPAD_ASSERT_UNKNOWN((op == InvOp) | (op == BeginOp));
00920                 }
00921                 else    CPPAD_ASSERT_UNKNOWN((op != InvOp) & (op != BeginOp));
00922 # endif
00923                 switch( op )
00924                 {
00925                         // Unary operator where operand is arg[0]
00926                         case AbsOp:
00927                         case AcosOp:
00928                         case AsinOp:
00929                         case AtanOp:
00930                         case CosOp:
00931                         case CoshOp:
00932                         case DivvpOp:
00933                         case ExpOp:
00934                         case LogOp:
00935                         case PowvpOp:
00936                         case SinOp:
00937                         case SinhOp:
00938                         case SqrtOp:
00939                         if( tape[i_var].connect != not_connected )
00940                                 tape[arg[0]].connect = yes_connected;
00941                         break; // --------------------------------------------
00942 
00943                         // Unary operator where operand is arg[1]
00944                         case DisOp:
00945                         case DivpvOp:
00946                         case MulpvOp:
00947                         case PowpvOp:
00948                         case PrivOp:
00949                         if( tape[i_var].connect != not_connected )
00950                                 tape[arg[1]].connect = yes_connected;
00951                         break; // --------------------------------------------
00952                 
00953                         // Special case for SubvpOp
00954                         case SubvpOp:
00955                         if( tape[i_var].connect != not_connected )
00956                         {
00957                                 if( tape[arg[0]].connect == not_connected )
00958                                         tape[arg[0]].connect = sum_connected;
00959                                 else
00960                                         tape[arg[0]].connect = yes_connected;
00961                                 if( tape[i_var].connect == sum_connected )
00962                                         tape[i_var].connect = csum_connected;
00963                         }
00964                         break; // --------------------------------------------
00965                 
00966                         // Special case for AddpvOp and SubpvOp
00967                         case AddpvOp:
00968                         case SubpvOp:
00969                         if( tape[i_var].connect != not_connected )
00970                         {
00971                                 if( tape[arg[1]].connect == not_connected )
00972                                         tape[arg[1]].connect = sum_connected;
00973                                 else
00974                                         tape[arg[1]].connect = yes_connected;
00975                                 if( tape[i_var].connect == sum_connected )
00976                                         tape[i_var].connect = csum_connected;
00977                         }
00978                         break; // --------------------------------------------
00979 
00980                 
00981                         // Special case for AddvvOp and SubvvOp
00982                         case AddvvOp:
00983                         case SubvvOp:
00984                         if( tape[i_var].connect != not_connected )
00985                         {
00986                                 if( tape[arg[0]].connect == not_connected )
00987                                         tape[arg[0]].connect = sum_connected;
00988                                 else
00989                                         tape[arg[0]].connect = yes_connected;
00990 
00991                                 if( tape[arg[1]].connect == not_connected )
00992                                         tape[arg[1]].connect = sum_connected;
00993                                 else
00994                                         tape[arg[1]].connect = yes_connected;
00995                                 if( tape[i_var].connect == sum_connected )
00996                                         tape[i_var].connect = csum_connected;
00997                         }
00998                         break; // --------------------------------------------
00999 
01000                         // Other binary operators 
01001                         // where operands are arg[0], arg[1]
01002                         case DivvvOp:
01003                         case MulvvOp:
01004                         case PowvvOp:
01005                         if( tape[i_var].connect != not_connected )
01006                         {
01007                                 tape[arg[0]].connect = yes_connected;
01008                                 tape[arg[1]].connect = yes_connected;
01009                         }
01010                         break; // --------------------------------------------
01011 
01012                         // Conditional expression operators
01013                         case CExpOp:
01014                         CPPAD_ASSERT_UNKNOWN( NumArg(CExpOp) == 6 );
01015                         if( tape[i_var].connect != not_connected )
01016                         {
01017                                 mask = 1;
01018                                 for(i = 2; i < 6; i++)
01019                                 {       if( arg[1] & mask )
01020                                         {       CPPAD_ASSERT_UNKNOWN( arg[i] < i_var );
01021                                                 tape[arg[i]].connect = yes_connected;
01022                                         }
01023                                         mask = mask << 1;
01024                                 }
01025                         }
01026                         break;  // --------------------------------------------
01027 
01028                         // Operations where there is noting to do
01029                         case BeginOp:
01030                         case ComOp:
01031                         case EndOp:
01032                         case InvOp:
01033                         case ParOp:
01034                         case PripOp:
01035                         case StppOp:
01036                         break;  // --------------------------------------------
01037 
01038                         // Load using a parameter index
01039                         case LdpOp:
01040                         if( tape[i_var].connect != not_connected )
01041                         {
01042                                 i                = vecad[ arg[0] - 1 ];
01043                                 vecad_connect[i] = yes_connected;
01044                         }
01045                         break; // --------------------------------------------
01046 
01047                         // Load using a variable index
01048                         case LdvOp:
01049                         if( tape[i_var].connect != not_connected )
01050                         {
01051                                 i                    = vecad[ arg[0] - 1 ];
01052                                 vecad_connect[i]     = yes_connected;
01053                                 tape[arg[1]].connect = yes_connected;
01054                         }
01055                         break; // --------------------------------------------
01056 
01057                         // Store a variable using a parameter index
01058                         case StpvOp:
01059                         i = vecad[ arg[0] - 1 ];
01060                         if( vecad_connect[i] != not_connected )
01061                                 tape[arg[2]].connect = yes_connected;
01062                         break; // --------------------------------------------
01063 
01064                         // Store a variable using a variable index
01065                         case StvvOp:
01066                         i = vecad[ arg[0] - 1 ];
01067                         if( vecad_connect[i] )
01068                         {       tape[arg[1]].connect = yes_connected;
01069                                 tape[arg[2]].connect = yes_connected;
01070                         }
01071                         break; 
01072                         // ============================================================
01073                         case UserOp:
01074                         // start or end atomic operation sequence
01075                         CPPAD_ASSERT_UNKNOWN( NumRes( UserOp ) == 0 );
01076                         CPPAD_ASSERT_UNKNOWN( NumArg( UserOp ) == 4 );
01077                         if( user_state == user_end )
01078                         {       user_index = arg[0];
01079                                 user_id    = arg[1];
01080                                 user_n     = arg[2];
01081                                 user_m     = arg[3];
01082                                 user_q     = 1;
01083                                 if(user_r.size() < user_n )
01084                                         user_r.resize(user_n);
01085                                 if(user_s.size() < user_m )
01086                                         user_s.resize(user_m);
01087                                 user_j     = user_n;
01088                                 user_i     = user_m;
01089                                 user_state = user_ret;
01090                                 user_keep.push(false);
01091                         }
01092                         else
01093                         {       CPPAD_ASSERT_UNKNOWN( user_state == user_start );
01094                                 CPPAD_ASSERT_UNKNOWN( user_index == arg[0] );
01095                                 CPPAD_ASSERT_UNKNOWN( user_id    == arg[1] );
01096                                 CPPAD_ASSERT_UNKNOWN( user_n     == arg[2] );
01097                                 CPPAD_ASSERT_UNKNOWN( user_m     == arg[3] );
01098                                 user_state = user_end;
01099                }
01100                         break;
01101 
01102                         case UsrapOp:
01103                         // parameter argument in an atomic operation sequence
01104                         CPPAD_ASSERT_UNKNOWN( user_state == user_arg );
01105                         CPPAD_ASSERT_UNKNOWN( 0 < user_j && user_j <= user_n );
01106                         CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
01107                         CPPAD_ASSERT_UNKNOWN( arg[0] < num_par );
01108                         --user_j;
01109                         if( user_j == 0 )
01110                                 user_state = user_start;
01111                         break;
01112 
01113                         case UsravOp:
01114                         // variable argument in an atomic operation sequence
01115                         CPPAD_ASSERT_UNKNOWN( user_state == user_arg );
01116                         CPPAD_ASSERT_UNKNOWN( 0 < user_j && user_j <= user_n );
01117                         CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
01118                         CPPAD_ASSERT_UNKNOWN( arg[0] <= i_var );
01119                         CPPAD_ASSERT_UNKNOWN( 0 < arg[0] );
01120                         --user_j;
01121                         if( ! user_r[user_j].empty() )
01122                         {       tape[arg[0]].connect = yes_connected;
01123                                 user_keep.top() = true;
01124                         }
01125                         if( user_j == 0 )
01126                                 user_state = user_start;
01127                         break;
01128 
01129                         case UsrrpOp:
01130                         // parameter result in an atomic operation sequence
01131                         CPPAD_ASSERT_UNKNOWN( user_state == user_ret );
01132                         CPPAD_ASSERT_UNKNOWN( 0 < user_i && user_i <= user_m );
01133                         CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
01134                         CPPAD_ASSERT_UNKNOWN( arg[0] < num_par );
01135                         --user_i;
01136                         user_s[user_i].clear();
01137                         if( user_i == 0 )
01138                         {       // call users function for this operation
01139                                 user_atomic<Base>::rev_jac_sparse(user_index, user_id,
01140                                         user_n, user_m, user_q, user_r, user_s
01141                                 );
01142                                 user_state = user_arg;
01143                         }
01144                         break;
01145 
01146                         case UsrrvOp:
01147                         // variable result in an atomic operation sequence
01148                         CPPAD_ASSERT_UNKNOWN( user_state == user_ret );
01149                         CPPAD_ASSERT_UNKNOWN( 0 < user_i && user_i <= user_m );
01150                         --user_i;
01151                         user_s[user_i].clear();
01152                         if( tape[i_var].connect != not_connected )
01153                                 user_s[user_i].insert(0);
01154                         if( user_i == 0 )
01155                         {       // call users function for this operation
01156                                 user_atomic<Base>::rev_jac_sparse(user_index, user_id,
01157                                         user_n, user_m, user_q, user_r, user_s
01158                                 );
01159                                 user_state = user_arg;
01160                         }
01161                         break;
01162                         // ============================================================
01163 
01164                         // all cases should be handled above
01165                         default:
01166                         CPPAD_ASSERT_UNKNOWN(0);
01167                 }
01168         }
01169         // values corresponding to BeginOp
01170         CPPAD_ASSERT_UNKNOWN( i_op == 0 && i_var == 0 && op == BeginOp );
01171         tape[i_var].op = op;
01172         // -------------------------------------------------------------
01173 
01174         // Erase all information in the recording
01175         rec->Erase();
01176 
01177         // Initilaize table mapping hash code to variable index in tape
01178         // as pointing to the BeginOp at the beginning of the tape
01179         CppAD::vector<size_t>  hash_table_var(CPPAD_HASH_TABLE_SIZE);
01180         for(i = 0; i < CPPAD_HASH_TABLE_SIZE; i++)
01181                 hash_table_var[i] = 0;
01182         CPPAD_ASSERT_UNKNOWN( tape[0].op == BeginOp );
01183 
01184         // initialize mapping from old variable index to new variable index
01185         for(i = 0; i < num_var; i++)
01186                 tape[i].new_var = num_var; // invalid index
01187         
01188 
01189         // initialize mapping from old VecAD index to new VecAD index
01190         CppAD::vector<size_t> new_vecad_ind(num_vecad_ind);
01191         for(i = 0; i < num_vecad_ind; i++)
01192                 new_vecad_ind[i] = num_vecad_ind; // invalid index 
01193 
01194         j = 0;     // index into the old set of indices
01195         for(i = 0; i < num_vecad_vec; i++)
01196         {       // length of this VecAD
01197                 size_t length = play->GetVecInd(j);
01198                 if( vecad_connect[i] != not_connected )
01199                 {       // Put this VecAD vector in new recording
01200                         CPPAD_ASSERT_UNKNOWN(length < num_vecad_ind);
01201                         new_vecad_ind[j] = rec->PutVecInd(length);
01202                         for(k = 1; k <= length; k++) new_vecad_ind[j+k] =
01203                                 rec->PutVecInd(
01204                                         rec->PutPar(
01205                                                 play->GetPar( 
01206                                                         play->GetVecInd(j+k)
01207                         ) ) );
01208                 }
01209                 // start of next VecAD
01210                 j       += length + 1;
01211         }
01212         CPPAD_ASSERT_UNKNOWN( j == num_vecad_ind );
01213 
01214         // start playing the operations in the forward direction
01215         play->start_forward(op, arg, i_op, i_var);
01216 
01217         // playing forward skips BeginOp at the beginning, but not EndOp at
01218         // the end.  Put BeginOp at beginning of recording
01219         CPPAD_ASSERT_UNKNOWN( op == BeginOp );
01220         CPPAD_ASSERT_NARG_NRES(BeginOp, 0, 1);
01221         tape[i_var].new_var = rec->PutOp(BeginOp);
01222 
01223         // temporary buffer for new argument values
01224         size_t new_arg[6];
01225 
01226         // temporary work space used by optimize_record_csum
01227         // (decalared here to avoid realloaction of memory)
01228         optimize_csum_stacks csum_work;
01229 
01230         user_state = user_start;
01231         while(op != EndOp)
01232         {       // next op
01233                 play->next_forward(op, arg, i_op, i_var);
01234                 CPPAD_ASSERT_UNKNOWN( (i_op > n)  | (op == InvOp) );
01235                 CPPAD_ASSERT_UNKNOWN( (i_op <= n) | (op != InvOp) );
01236 
01237                 // determine if we should keep this operation in the new
01238                 // operation sequence
01239                 bool keep;
01240                 switch( op )
01241                 {       case ComOp:
01242                         case PripOp:
01243                         case PrivOp:
01244                         keep = false;
01245                         break;
01246 
01247                         case InvOp:
01248                         case EndOp:
01249                         keep = true;
01250                         break;
01251 
01252                         case StppOp:
01253                         case StvpOp:
01254                         case StpvOp:
01255                         case StvvOp:
01256                         CPPAD_ASSERT_UNKNOWN( NumRes(op) == 0 );
01257                         i = vecad[ arg[0] - 1 ];
01258                         keep = vecad_connect[i] != not_connected;
01259                         break;
01260 
01261                         case AddpvOp:
01262                         case AddvvOp:
01263                         case SubpvOp:
01264                         case SubvpOp:
01265                         case SubvvOp:
01266                         keep  = tape[i_var].connect != not_connected;
01267                         keep &= tape[i_var].connect != csum_connected;
01268                         break; 
01269 
01270                         case UserOp:
01271                         case UsrapOp:
01272                         case UsravOp:
01273                         case UsrrpOp:
01274                         case UsrrvOp:
01275                         keep = user_keep.top();
01276                         break;
01277 
01278                         default:
01279                         keep = tape[i_var].connect != not_connected;
01280                         break;
01281                 }
01282 
01283                 size_t         match_var    = 0;
01284                 unsigned short code         = 0;
01285                 bool           replace_hash = false;
01286                 if( keep ) switch( op )
01287                 {
01288                         // Unary operator where operand is arg[0]
01289                         case AbsOp:
01290                         case AcosOp:
01291                         case AsinOp:
01292                         case AtanOp:
01293                         case CosOp:
01294                         case CoshOp:
01295                         case ExpOp:
01296                         case LogOp:
01297                         case SinOp:
01298                         case SinhOp:
01299                         case SqrtOp:
01300                         match_var = optimize_unary_match(
01301                                 tape                ,  // inputs 
01302                                 i_var               ,
01303                                 play->num_rec_par() ,
01304                                 play->GetPar()      ,
01305                                 hash_table_var      ,
01306                                 code                  // outputs
01307                         );
01308                         if( match_var > 0 )
01309                                 tape[i_var].new_var = match_var;
01310                         else
01311                         {
01312                                 replace_hash = true;
01313                                 new_arg[0]   = tape[ arg[0] ].new_var;
01314                                 rec->PutArg( new_arg[0] );
01315                                 i                   = rec->PutOp(op);
01316                                 tape[i_var].new_var = i;
01317                                 CPPAD_ASSERT_UNKNOWN( new_arg[0] < i );
01318                         }
01319                         break;
01320                         // ---------------------------------------------------
01321                         // Binary operators where 
01322                         // left is a variable and right is a parameter
01323                         case SubvpOp:
01324                         if( tape[arg[0]].connect == csum_connected )
01325                         {
01326                                 // convert to a sequence of summation operators
01327                                 tape[i_var].new_var = optimize_record_csum(
01328                                         tape                , // inputs
01329                                         i_var               ,
01330                                         play->num_rec_par() ,
01331                                         play->GetPar()      ,
01332                                         rec                 ,
01333                                         csum_work
01334                                 );
01335                                 // abort rest of this case
01336                                 break;
01337                         }
01338                         case DivvpOp:
01339                         case PowvpOp:
01340                         match_var = optimize_binary_match(
01341                                 tape                ,  // inputs 
01342                                 i_var               ,
01343                                 play->num_rec_par() ,
01344                                 play->GetPar()      ,
01345                                 hash_table_var      ,
01346                                 code                  // outputs
01347                         );
01348                         if( match_var > 0 )
01349                                 tape[i_var].new_var = match_var;
01350                         else
01351                         {       tape[i_var].new_var = optimize_record_vp(
01352                                         tape                , // inputs
01353                                         i_var               ,
01354                                         play->num_rec_par() ,
01355                                         play->GetPar()      ,
01356                                         rec                 ,
01357                                         op                  ,
01358                                         arg
01359                                 );
01360                                 replace_hash = true;
01361                         }
01362                         break;
01363                         // ---------------------------------------------------
01364                         // Binary operators where 
01365                         // left is a parameter and right is a variable
01366                         case SubpvOp:
01367                         case AddpvOp:
01368                         if( tape[arg[1]].connect == csum_connected )
01369                         {
01370                                 // convert to a sequence of summation operators
01371                                 tape[i_var].new_var = optimize_record_csum(
01372                                         tape                , // inputs
01373                                         i_var               ,
01374                                         play->num_rec_par() ,
01375                                         play->GetPar()      ,
01376                                         rec                 ,
01377                                         csum_work
01378                                 );
01379                                 // abort rest of this case
01380                                 break;
01381                         }
01382                         case DivpvOp:
01383                         case MulpvOp:
01384                         case PowpvOp:
01385                         match_var = optimize_binary_match(
01386                                 tape                ,  // inputs 
01387                                 i_var               ,
01388                                 play->num_rec_par() ,
01389                                 play->GetPar()      ,
01390                                 hash_table_var      ,
01391                                 code                  // outputs
01392                         );
01393                         if( match_var > 0 )
01394                                 tape[i_var].new_var = match_var;
01395                         else
01396                         {       tape[i_var].new_var = optimize_record_pv(
01397                                         tape                , // inputs
01398                                         i_var               ,
01399                                         play->num_rec_par() ,
01400                                         play->GetPar()      ,
01401                                         rec                 ,
01402                                         op                  ,
01403                                         arg
01404                                 );
01405                                 replace_hash = true;
01406                         }
01407                         break;
01408                         // ---------------------------------------------------
01409                         // Binary operator where 
01410                         // both operators are variables
01411                         case AddvvOp:
01412                         case SubvvOp:
01413                         if( (tape[arg[0]].connect == csum_connected) |
01414                             (tape[arg[1]].connect == csum_connected)
01415                         )
01416                         {
01417                                 // convert to a sequence of summation operators
01418                                 tape[i_var].new_var = optimize_record_csum(
01419                                         tape                , // inputs
01420                                         i_var               ,
01421                                         play->num_rec_par() ,
01422                                         play->GetPar()      ,
01423                                         rec                 ,
01424                                         csum_work
01425                                 );
01426                                 // abort rest of this case
01427                                 break;
01428                         }
01429                         case DivvvOp:
01430                         case MulvvOp:
01431                         case PowvvOp:
01432                         match_var = optimize_binary_match(
01433                                 tape                ,  // inputs 
01434                                 i_var               ,
01435                                 play->num_rec_par() ,
01436                                 play->GetPar()      ,
01437                                 hash_table_var      ,
01438                                 code                  // outputs
01439                         );
01440                         if( match_var > 0 )
01441                                 tape[i_var].new_var = match_var;
01442                         else
01443                         {       tape[i_var].new_var = optimize_record_vv(
01444                                         tape                , // inputs
01445                                         i_var               ,
01446                                         play->num_rec_par() ,
01447                                         play->GetPar()      ,
01448                                         rec                 ,
01449                                         op                  ,
01450                                         arg
01451                                 );
01452                                 replace_hash = true;
01453                         }
01454                         break;
01455                         // ---------------------------------------------------
01456                         // Conditional expression operators
01457                         case CExpOp:
01458                         CPPAD_ASSERT_NARG_NRES(op, 6, 1);
01459                         new_arg[0] = arg[0];
01460                         new_arg[1] = arg[1];
01461                         mask = 1;
01462                         for(i = 2; i < 6; i++)
01463                         {       if( arg[1] & mask )
01464                                 {       new_arg[i] = tape[arg[i]].new_var;
01465                                         CPPAD_ASSERT_UNKNOWN( 
01466                                                 new_arg[i] < num_var 
01467                                         );
01468                                 }
01469                                 else    new_arg[i] = rec->PutPar( 
01470                                                 play->GetPar( arg[i] )
01471                                 );
01472                                 mask = mask << 1;
01473                         }
01474                         rec->PutArg(
01475                                 new_arg[0] ,
01476                                 new_arg[1] ,
01477                                 new_arg[2] ,
01478                                 new_arg[3] ,
01479                                 new_arg[4] ,
01480                                 new_arg[5] 
01481                         );
01482                         tape[i_var].new_var = rec->PutOp(op);
01483                         break;
01484                         // ---------------------------------------------------
01485                         // Operations with no arguments and no results
01486                         case EndOp:
01487                         CPPAD_ASSERT_NARG_NRES(op, 0, 0);
01488                         rec->PutOp(op);
01489                         break;
01490                         // ---------------------------------------------------
01491                         // Operations with no arguments and one result
01492                         case InvOp:
01493                         CPPAD_ASSERT_NARG_NRES(op, 0, 1);
01494                         tape[i_var].new_var = rec->PutOp(op);
01495                         break;
01496                         // ---------------------------------------------------
01497                         // Operations with one argument that is a parameter
01498                         case ParOp:
01499                         CPPAD_ASSERT_NARG_NRES(op, 1, 1);
01500                         new_arg[0] = rec->PutPar( play->GetPar(arg[0] ) );
01501 
01502                         rec->PutArg( new_arg[0] );
01503                         tape[i_var].new_var = rec->PutOp(op);
01504                         break;
01505                         // ---------------------------------------------------
01506                         // Load using a parameter index
01507                         case LdpOp:
01508                         CPPAD_ASSERT_NARG_NRES(op, 3, 1);
01509                         new_arg[0] = new_vecad_ind[ arg[0] ];
01510                         new_arg[1] = arg[1];
01511                         CPPAD_ASSERT_UNKNOWN( new_arg[0] < num_vecad_ind );
01512                         rec->PutArg( 
01513                                 new_arg[0], 
01514                                 new_arg[1], 
01515                                 0
01516                         );
01517                         tape[i_var].new_var = rec->PutOp(op);
01518                         break;
01519                         // ---------------------------------------------------
01520                         // Load using a variable index
01521                         case LdvOp:
01522                         CPPAD_ASSERT_NARG_NRES(op, 3, 1);
01523                         new_arg[0] = new_vecad_ind[ arg[0] ];
01524                         new_arg[1] = tape[arg[1]].new_var;
01525                         CPPAD_ASSERT_UNKNOWN( new_arg[0] < num_vecad_ind );
01526                         CPPAD_ASSERT_UNKNOWN( new_arg[1] < num_var );
01527                         rec->PutArg( 
01528                                 new_arg[0], 
01529                                 new_arg[1], 
01530                                 0
01531                         );
01532                         tape[i_var].new_var = rec->PutOp(op);
01533                         break;
01534                         // ---------------------------------------------------
01535                         // Store a parameter using a parameter index
01536                         case StppOp:
01537                         CPPAD_ASSERT_NARG_NRES(op, 3, 0);
01538                         new_arg[0] = new_vecad_ind[ arg[0] ];
01539                         new_arg[1] = rec->PutPar( play->GetPar(arg[1]) );
01540                         new_arg[2] = rec->PutPar( play->GetPar(arg[2]) );
01541                         CPPAD_ASSERT_UNKNOWN( new_arg[0] < num_vecad_ind );
01542                         rec->PutArg(
01543                                 new_arg[0], 
01544                                 new_arg[1], 
01545                                 new_arg[2]
01546                         );
01547                         rec->PutOp(op);
01548                         break;
01549                         // ---------------------------------------------------
01550                         // Store a parameter using a variable index
01551                         case StvpOp:
01552                         CPPAD_ASSERT_NARG_NRES(op, 3, 0);
01553                         new_arg[0] = new_vecad_ind[ arg[0] ];
01554                         new_arg[1] = tape[arg[1]].new_var;
01555                         new_arg[2] = rec->PutPar( play->GetPar(arg[2]) );
01556                         CPPAD_ASSERT_UNKNOWN( new_arg[0] < num_vecad_ind );
01557                         CPPAD_ASSERT_UNKNOWN( new_arg[1] < num_var );
01558                         rec->PutArg(
01559                                 new_arg[0], 
01560                                 new_arg[1], 
01561                                 new_arg[2]
01562                         );
01563                         rec->PutOp(op);
01564                         break;
01565                         // ---------------------------------------------------
01566                         // Store a variable using a parameter index
01567                         case StpvOp:
01568                         CPPAD_ASSERT_NARG_NRES(op, 3, 0);
01569                         new_arg[0] = new_vecad_ind[ arg[0] ];
01570                         new_arg[1] = rec->PutPar( play->GetPar(arg[1]) );
01571                         new_arg[2] = tape[arg[2]].new_var;
01572                         CPPAD_ASSERT_UNKNOWN( new_arg[0] < num_vecad_ind );
01573                         CPPAD_ASSERT_UNKNOWN( new_arg[1] < num_var );
01574                         CPPAD_ASSERT_UNKNOWN( new_arg[2] < num_var );
01575                         rec->PutArg(
01576                                 new_arg[0], 
01577                                 new_arg[1], 
01578                                 new_arg[2]
01579                         );
01580                         rec->PutOp(op);
01581                         break;
01582                         // ---------------------------------------------------
01583                         // Store a variable using a variable index
01584                         case StvvOp:
01585                         CPPAD_ASSERT_NARG_NRES(op, 3, 0);
01586                         new_arg[0] = new_vecad_ind[ arg[0] ];
01587                         new_arg[1] = tape[arg[1]].new_var;
01588                         new_arg[2] = tape[arg[2]].new_var;
01589                         CPPAD_ASSERT_UNKNOWN( new_arg[0] < num_vecad_ind );
01590                         CPPAD_ASSERT_UNKNOWN( new_arg[1] < num_var );
01591                         CPPAD_ASSERT_UNKNOWN( new_arg[2] < num_var );
01592                         rec->PutArg(
01593                                 new_arg[0], 
01594                                 new_arg[1], 
01595                                 new_arg[2]
01596                         );
01597                         rec->PutOp(op);
01598                         break;
01599 
01600                         // -----------------------------------------------------------
01601                         case UserOp:
01602                         CPPAD_ASSERT_NARG_NRES(op, 4, 0);
01603                         if( user_state == user_start )
01604                                 user_state = user_arg;
01605                         else
01606                         {       user_state = user_start;
01607                                 user_keep.pop();        
01608                         }
01609                         // user_index, user_id, user_n, user_m
01610                         rec->PutArg(arg[0], arg[1], arg[2], arg[3]);
01611                         rec->PutOp(UserOp);
01612                         break;
01613 
01614                         case UsrapOp:
01615                         CPPAD_ASSERT_NARG_NRES(op, 1, 0);
01616                         new_arg[0] = rec->PutPar( play->GetPar(arg[0]) );
01617                         rec->PutArg(new_arg[0]);
01618                         rec->PutOp(UsrapOp);
01619                         break;
01620 
01621                         case UsravOp:
01622                         CPPAD_ASSERT_NARG_NRES(op, 1, 0);
01623                         new_arg[0] = tape[arg[0]].new_var;
01624                         rec->PutArg(new_arg[0]);
01625                         rec->PutOp(UsravOp);
01626                         break;
01627 
01628                         case UsrrpOp:
01629                         CPPAD_ASSERT_NARG_NRES(op, 1, 0);
01630                         new_arg[0] = rec->PutPar( play->GetPar(arg[0]) );
01631                         rec->PutArg(new_arg[0]);
01632                         rec->PutOp(UsrrpOp);
01633                         break;
01634                         
01635                         case UsrrvOp:
01636                         CPPAD_ASSERT_NARG_NRES(op, 0, 1);
01637                         tape[i_var].new_var = rec->PutOp(UsrrvOp);
01638                         break;
01639                         // ---------------------------------------------------
01640 
01641                         // all cases should be handled above
01642                         default:
01643                         CPPAD_ASSERT_UNKNOWN(false);
01644 
01645                 }
01646                 if( replace_hash )
01647                 {       // The old variable index i_var corresponds to the 
01648                         // new variable index tape[i_var].new_var. In addition
01649                         // this is the most recent variable that has this code.
01650                         hash_table_var[code] = i_var;
01651                 }
01652 
01653         }
01654         // modify the dependent variable vector to new indices
01655         for(i = 0; i < dep_taddr.size(); i++ )
01656         {       CPPAD_ASSERT_UNKNOWN( tape[ dep_taddr[i] ].new_var < num_var );
01657                 dep_taddr[i] = tape[ dep_taddr[i] ].new_var;
01658         }
01659 }
01660 
01661 /*!
01662 Optimize a player object operation sequence
01663 
01664 The operation sequence for this object is replaced by one with fewer operations
01665 but the same funcition and derivative values.
01666 
01667 \tparam Base
01668 base type for the operator; i.e., this operation was recorded
01669 using AD< \a Base > and computations by this routine are done using type 
01670 \a Base.
01671 
01672 */
01673 template <class Base>
01674 void ADFun<Base>::optimize(void)
01675 {       // place to store the optimized version of the recording
01676         recorder<Base> rec;
01677 
01678         // number of independent variables
01679         size_t n = ind_taddr_.size();
01680 
01681 # ifndef NDEBUG
01682         size_t i, j, m = dep_taddr_.size();
01683         CppAD::vector<Base> x(n), y(m), check(m);
01684         bool check_zero_order = taylor_per_var_ > 0;
01685         if( check_zero_order )
01686         {       // zero order coefficients for independent vars
01687                 for(j = 0; j < n; j++)
01688                 {       CPPAD_ASSERT_UNKNOWN( play_.GetOp(j+1) == InvOp );
01689                         CPPAD_ASSERT_UNKNOWN( ind_taddr_[j]    == j+1   );
01690                         x[j] = taylor_[ ind_taddr_[j] * taylor_col_dim_ + 0];
01691                 }
01692                 // zero order coefficients for dependent vars
01693                 for(i = 0; i < m; i++)
01694                 {       CPPAD_ASSERT_UNKNOWN( dep_taddr_[i] < total_num_var_ );
01695                         y[i] = taylor_[ dep_taddr_[i] * taylor_col_dim_ + 0];
01696                 }
01697         }
01698 # endif
01699 
01700         // create the optimized recording
01701         CppAD::optimize<Base>(n, dep_taddr_, &play_, &rec);
01702 
01703         // now replace the recording
01704         play_ = rec;
01705 
01706         // number of variables in the recording
01707         total_num_var_ = rec.num_rec_var();
01708 
01709         // free memory allocated for sparse Jacobian calculation
01710         // (the results are no longer valid)
01711         for_jac_sparse_pack_.resize(0, 0);
01712         for_jac_sparse_set_.resize(0,0);
01713 
01714         // free old Taylor coefficient memory
01715         if( taylor_ != CPPAD_NULL )
01716                 CPPAD_TRACK_DEL_VEC(taylor_);
01717         taylor_         = CPPAD_NULL;
01718         taylor_per_var_ = 0;
01719         taylor_col_dim_ = 0;
01720 
01721 # ifndef NDEBUG
01722         if( check_zero_order )
01723         {
01724                 // zero order forward calculation using new operation sequence
01725                 check = Forward(0, x);
01726 
01727                 // check results
01728                 for(i = 0; i < m; i++) CPPAD_ASSERT_KNOWN( 
01729                         check[i] == y[i] ,
01730                         "Error during check of f.optimize()."
01731                 );
01732 
01733                 // Erase memory that this calculation was done so NDEBUG gives 
01734                 // same final state for this object (from users perspective)
01735                 taylor_per_var_ = 0;
01736         }
01737 # endif
01738 }
01739 
01740 CPPAD_END_NAMESPACE
01741 # endif