CppAD: A C++ Algorithmic Differentiation Package 20110419
|
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