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