CppAD: A C++ Algorithmic Differentiation Package
20130102
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_REVERSE_SWEEP_INCLUDED 00003 # define CPPAD_REVERSE_SWEEP_INCLUDED 00004 00005 /* -------------------------------------------------------------------------- 00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 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 CPPAD_BEGIN_NAMESPACE 00018 /*! 00019 \defgroup reverse_sweep_hpp reverse_sweep.hpp 00020 \{ 00021 \file reverse_sweep.hpp 00022 Compute derivatives of arbitrary order Taylor coefficients. 00023 */ 00024 00025 /*! 00026 \def CPPAD_REVERSE_SWEEP_TRACE 00027 This value is either zero or one. 00028 Zero is the normal operational value. 00029 If it is one, a trace of every reverse_sweep computation is printed. 00030 */ 00031 # define CPPAD_REVERSE_SWEEP_TRACE 0 00032 00033 /*! 00034 Compute derivative of arbitrary order forward mode Taylor coefficients. 00035 00036 \tparam Base 00037 base type for the operator; i.e., this operation sequence was recorded 00038 using AD< \a Base > and computations by this routine are done using type 00039 \a Base. 00040 00041 \param d 00042 is the highest order Taylor coefficients that 00043 we are computing the derivative of. 00044 00045 \param n 00046 is the number of independent variables on the tape. 00047 00048 \param numvar 00049 is the total number of variables on the tape. 00050 This is also equal to the number of rows in the matrix \a Taylor; i.e., 00051 Rec->num_rec_var(). 00052 00053 \param Rec 00054 The information stored in \a Rec 00055 is a recording of the operations corresponding to the function 00056 \f[ 00057 F : {\bf R}^n \rightarrow {\bf R}^m 00058 \f] 00059 where \f$ n \f$ is the number of independent variables and 00060 \f$ m \f$ is the number of dependent variables. 00061 We define the function 00062 \f$ G : {\bf R}^{n \times d} \rightarrow {\bf R} \f$ by 00063 \f[ 00064 G( u ) = \frac{1}{d !} \frac{ \partial^d }{ \partial t^d } 00065 \left[ 00066 \sum_{i=1}^m w_i F_i ( u^{(0)} + u^{(1)} t + \cdots + u^{(d)} t^d ) 00067 \right]_{t=0} 00068 \f] 00069 Note that the scale factor 1 / a d converts 00070 the \a d-th partial derivative to the \a d-th order Taylor coefficient. 00071 This routine computes the derivative of \f$ G(u) \f$ 00072 with respect to all the Taylor coefficients 00073 \f$ u^{(k)} \f$ for \f$ k = 0 , ... , d \f$. 00074 The vector \f$ w \in {\bf R}^m \f$, and 00075 value of \f$ u \in {\bf R}^{n \times d} \f$ 00076 at which the derivative is computed, 00077 are defined below. 00078 \n 00079 \n 00080 The object \a Rec is effectly constant. 00081 There is an exception to this, 00082 while palying back the tape 00083 the object \a Rec holds information about the current location 00084 with in the tape and this changes during palyback. 00085 00086 \param J 00087 Is the number of columns in the coefficient matrix \a Taylor. 00088 This must be greater than or equal \a d + 1. 00089 00090 \param Taylor 00091 For i = 1 , ... , \a numvar, and for k = 0 , ... , \a d, 00092 \a Taylor [ i * J + k ] 00093 is the k-th order Taylor coefficient corresponding to 00094 variable with index i on the tape. 00095 The value \f$ u \in {\bf R}^{n \times d} \f$, 00096 at which the derivative is computed, 00097 is defined by 00098 \f$ u_j^{(k)} \f$ = \a Taylor [ j * J + k ] 00099 for j = 1 , ... , \a n, and for k = 0 , ... , \a d. 00100 00101 \param K 00102 Is the number of columns in the partial derivative matrix \a Partial. 00103 It must be greater than or equal \a d + 1. 00104 00105 \param Partial 00106 \b Input: 00107 The last \f$ m \f$ rows of \a Partial are inputs. 00108 The vector \f$ v \f$, used to define \f$ G(u) \f$, 00109 is specified by these rows. 00110 For i = 0 , ... , m - 1, \a Partial [ ( \a numvar - m + i ) * K + d ] = v_i. 00111 For i = 0 , ... , m - 1 and for k = 0 , ... , d - 1, 00112 \a Partial [ ( \a numvar - m + i ) * K + k ] = 0. 00113 \n 00114 \n 00115 \b Temporary: 00116 For i = n+1 , ... , \a numvar - 1 and for k = 0 , ... , d, 00117 the value of \a Partial [ i * K + k ] is used for temporary work space 00118 and its output value is not defined. 00119 \n 00120 \n 00121 \b Output: 00122 For j = 1 , ... , n and for k = 0 , ... , d, 00123 \a Partial [ j * K + k ] 00124 is the partial derivative of \f$ G( u ) \f$ with 00125 respect to \f$ u_j^{(k)} \f$. 00126 00127 \par Assumptions 00128 The first operator on the tape is a BeginOp, 00129 and the next \a n operators are InvOp operations for the 00130 corresponding independent variables. 00131 */ 00132 template <class Base> 00133 void ReverseSweep( 00134 size_t d, 00135 size_t n, 00136 size_t numvar, 00137 player<Base>* Rec, 00138 size_t J, 00139 const Base* Taylor, 00140 size_t K, 00141 Base* Partial 00142 ) 00143 { 00144 OpCode op; 00145 size_t i_op; 00146 size_t i_var; 00147 00148 const addr_t* arg = 0; 00149 00150 // check numvar argument 00151 CPPAD_ASSERT_UNKNOWN( Rec->num_rec_var() == numvar ); 00152 CPPAD_ASSERT_UNKNOWN( numvar > 0 ); 00153 00154 // length of the parameter vector (used by CppAD assert macros) 00155 const size_t num_par = Rec->num_rec_par(); 00156 00157 // pointer to the beginning of the parameter vector 00158 const Base* parameter = 0; 00159 if( num_par > 0 ) 00160 parameter = Rec->GetPar(); 00161 00162 // work space used by UserOp. 00163 const size_t user_k = d; // order of this forward mode calculation 00164 const size_t user_k1 = d+1; // number of orders for this calculation 00165 vector<size_t> user_ix; // variable indices for argument vector 00166 vector<Base> user_tx; // argument vector Taylor coefficients 00167 vector<Base> user_ty; // result vector Taylor coefficients 00168 vector<Base> user_px; // partials w.r.t argument vector 00169 vector<Base> user_py; // partials w.r.t. result vector 00170 size_t user_index = 0; // indentifier for this user_atomic operation 00171 size_t user_id = 0; // user identifier for this call to operator 00172 size_t user_i = 0; // index in result vector 00173 size_t user_j = 0; // index in argument vector 00174 size_t user_m = 0; // size of result vector 00175 size_t user_n = 0; // size of arugment vector 00176 // next expected operator in a UserOp sequence 00177 enum { user_start, user_arg, user_ret, user_end } user_state = user_end; 00178 00179 // temporary indices 00180 size_t j, ell; 00181 00182 // Initialize 00183 Rec->start_reverse(op, arg, i_op, i_var); 00184 CPPAD_ASSERT_UNKNOWN( op == EndOp ); 00185 # if CPPAD_REVERSE_SWEEP_TRACE 00186 std::cout << std::endl; 00187 # endif 00188 while(op != BeginOp ) 00189 { // next op 00190 Rec->next_reverse(op, arg, i_op, i_var); 00191 # ifndef NDEBUG 00192 if( i_op <= n ) 00193 { CPPAD_ASSERT_UNKNOWN((op == InvOp) | (op == BeginOp)); 00194 } 00195 else CPPAD_ASSERT_UNKNOWN((op != InvOp) & (op != BeginOp)); 00196 # endif 00197 00198 // rest of informaiton depends on the case 00199 # if CPPAD_REVERSE_SWEEP_TRACE 00200 size_t i_tmp = i_var; 00201 const Base* Z_tmp = Taylor + i_var * J; 00202 const Base* pZ_tmp = Partial + i_var * K; 00203 00204 printOp( 00205 std::cout, 00206 Rec, 00207 i_tmp, 00208 op, 00209 arg, 00210 d + 1, 00211 Z_tmp, 00212 d + 1, 00213 pZ_tmp 00214 ); 00215 # endif 00216 00217 switch( op ) 00218 { 00219 00220 case AbsOp: 00221 reverse_abs_op( 00222 d, i_var, arg[0], J, Taylor, K, Partial 00223 ); 00224 break; 00225 // -------------------------------------------------- 00226 00227 case AddvvOp: 00228 reverse_addvv_op( 00229 d, i_var, arg, parameter, J, Taylor, K, Partial 00230 ); 00231 break; 00232 // -------------------------------------------------- 00233 00234 case AddpvOp: 00235 CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par ); 00236 reverse_addpv_op( 00237 d, i_var, arg, parameter, J, Taylor, K, Partial 00238 ); 00239 break; 00240 // -------------------------------------------------- 00241 00242 case AcosOp: 00243 // sqrt(1 - x * x), acos(x) 00244 CPPAD_ASSERT_UNKNOWN( i_var < numvar ); 00245 reverse_acos_op( 00246 d, i_var, arg[0], J, Taylor, K, Partial 00247 ); 00248 break; 00249 // -------------------------------------------------- 00250 00251 case AsinOp: 00252 // sqrt(1 - x * x), asin(x) 00253 CPPAD_ASSERT_UNKNOWN( i_var < numvar ); 00254 reverse_asin_op( 00255 d, i_var, arg[0], J, Taylor, K, Partial 00256 ); 00257 break; 00258 // -------------------------------------------------- 00259 00260 case AtanOp: 00261 // 1 + x * x, atan(x) 00262 CPPAD_ASSERT_UNKNOWN( i_var < numvar ); 00263 reverse_atan_op( 00264 d, i_var, arg[0], J, Taylor, K, Partial 00265 ); 00266 break; 00267 // ------------------------------------------------- 00268 00269 case BeginOp: 00270 CPPAD_ASSERT_NARG_NRES(op, 0, 1); 00271 break; 00272 // -------------------------------------------------- 00273 00274 case CSumOp: 00275 // CSumOp has a variable number of arguments and 00276 // next_reverse thinks it one has one argument. 00277 // We must inform next_reverse of this special case. 00278 Rec->reverse_csum(op, arg, i_op, i_var); 00279 reverse_csum_op( 00280 d, i_var, arg, K, Partial 00281 ); 00282 // end of a cummulative summation 00283 break; 00284 // ------------------------------------------------- 00285 00286 case CExpOp: 00287 reverse_cond_op( 00288 d, 00289 i_var, 00290 arg, 00291 num_par, 00292 parameter, 00293 J, 00294 Taylor, 00295 K, 00296 Partial 00297 ); 00298 break; 00299 // -------------------------------------------------- 00300 00301 case ComOp: 00302 break; 00303 // -------------------------------------------------- 00304 00305 case CosOp: 00306 CPPAD_ASSERT_UNKNOWN( i_var < numvar ); 00307 reverse_cos_op( 00308 d, i_var, arg[0], J, Taylor, K, Partial 00309 ); 00310 break; 00311 // -------------------------------------------------- 00312 00313 case CoshOp: 00314 CPPAD_ASSERT_UNKNOWN( i_var < numvar ); 00315 reverse_cosh_op( 00316 d, i_var, arg[0], J, Taylor, K, Partial 00317 ); 00318 break; 00319 // -------------------------------------------------- 00320 00321 case DisOp: 00322 // Derivative of discrete operation is zero so no 00323 // contribution passes through this operation. 00324 break; 00325 // -------------------------------------------------- 00326 00327 case DivvvOp: 00328 reverse_divvv_op( 00329 d, i_var, arg, parameter, J, Taylor, K, Partial 00330 ); 00331 break; 00332 // -------------------------------------------------- 00333 00334 case DivpvOp: 00335 CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par ); 00336 reverse_divpv_op( 00337 d, i_var, arg, parameter, J, Taylor, K, Partial 00338 ); 00339 break; 00340 // -------------------------------------------------- 00341 00342 case DivvpOp: 00343 CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < num_par ); 00344 reverse_divvp_op( 00345 d, i_var, arg, parameter, J, Taylor, K, Partial 00346 ); 00347 break; 00348 // -------------------------------------------------- 00349 00350 case ExpOp: 00351 reverse_exp_op( 00352 d, i_var, arg[0], J, Taylor, K, Partial 00353 ); 00354 break; 00355 // -------------------------------------------------- 00356 case LdpOp: 00357 reverse_load_op( 00358 op, d, i_var, arg, J, Taylor, K, Partial 00359 ); 00360 break; 00361 // ------------------------------------------------- 00362 00363 case LdvOp: 00364 reverse_load_op( 00365 op, d, i_var, arg, J, Taylor, K, Partial 00366 ); 00367 break; 00368 // ------------------------------------------------- 00369 00370 case InvOp: 00371 break; 00372 // -------------------------------------------------- 00373 00374 case LogOp: 00375 reverse_log_op( 00376 d, i_var, arg[0], J, Taylor, K, Partial 00377 ); 00378 break; 00379 // -------------------------------------------------- 00380 00381 case MulvvOp: 00382 reverse_mulvv_op( 00383 d, i_var, arg, parameter, J, Taylor, K, Partial 00384 ); 00385 break; 00386 // -------------------------------------------------- 00387 00388 case MulpvOp: 00389 CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par ); 00390 reverse_mulpv_op( 00391 d, i_var, arg, parameter, J, Taylor, K, Partial 00392 ); 00393 break; 00394 // -------------------------------------------------- 00395 00396 case ParOp: 00397 break; 00398 // -------------------------------------------------- 00399 00400 case PowvpOp: 00401 CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < num_par ); 00402 reverse_powvp_op( 00403 d, i_var, arg, parameter, J, Taylor, K, Partial 00404 ); 00405 break; 00406 // ------------------------------------------------- 00407 00408 case PowpvOp: 00409 CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par ); 00410 reverse_powpv_op( 00411 d, i_var, arg, parameter, J, Taylor, K, Partial 00412 ); 00413 break; 00414 // ------------------------------------------------- 00415 00416 case PowvvOp: 00417 reverse_powvv_op( 00418 d, i_var, arg, parameter, J, Taylor, K, Partial 00419 ); 00420 break; 00421 // -------------------------------------------------- 00422 00423 case PriOp: 00424 // no result so nothing to do 00425 break; 00426 // -------------------------------------------------- 00427 00428 case SignOp: 00429 CPPAD_ASSERT_UNKNOWN( i_var < numvar ); 00430 reverse_sign_op( 00431 d, i_var, arg[0], J, Taylor, K, Partial 00432 ); 00433 break; 00434 // ------------------------------------------------- 00435 00436 case SinOp: 00437 CPPAD_ASSERT_UNKNOWN( i_var < numvar ); 00438 reverse_sin_op( 00439 d, i_var, arg[0], J, Taylor, K, Partial 00440 ); 00441 break; 00442 // ------------------------------------------------- 00443 00444 case SinhOp: 00445 CPPAD_ASSERT_UNKNOWN( i_var < numvar ); 00446 reverse_sinh_op( 00447 d, i_var, arg[0], J, Taylor, K, Partial 00448 ); 00449 break; 00450 // -------------------------------------------------- 00451 00452 case SqrtOp: 00453 reverse_sqrt_op( 00454 d, i_var, arg[0], J, Taylor, K, Partial 00455 ); 00456 break; 00457 // -------------------------------------------------- 00458 00459 case StppOp: 00460 break; 00461 // -------------------------------------------------- 00462 00463 case StpvOp: 00464 break; 00465 // ------------------------------------------------- 00466 00467 case StvpOp: 00468 break; 00469 // ------------------------------------------------- 00470 00471 case StvvOp: 00472 break; 00473 // -------------------------------------------------- 00474 00475 case SubvvOp: 00476 reverse_subvv_op( 00477 d, i_var, arg, parameter, J, Taylor, K, Partial 00478 ); 00479 break; 00480 // -------------------------------------------------- 00481 00482 case SubpvOp: 00483 CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par ); 00484 reverse_subpv_op( 00485 d, i_var, arg, parameter, J, Taylor, K, Partial 00486 ); 00487 break; 00488 // -------------------------------------------------- 00489 00490 case SubvpOp: 00491 CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < num_par ); 00492 reverse_subvp_op( 00493 d, i_var, arg, parameter, J, Taylor, K, Partial 00494 ); 00495 break; 00496 // ------------------------------------------------- 00497 00498 case TanOp: 00499 CPPAD_ASSERT_UNKNOWN( i_var < numvar ); 00500 reverse_tan_op( 00501 d, i_var, arg[0], J, Taylor, K, Partial 00502 ); 00503 break; 00504 // ------------------------------------------------- 00505 00506 case TanhOp: 00507 CPPAD_ASSERT_UNKNOWN( i_var < numvar ); 00508 reverse_tanh_op( 00509 d, i_var, arg[0], J, Taylor, K, Partial 00510 ); 00511 break; 00512 // -------------------------------------------------- 00513 00514 case UserOp: 00515 // start or end an atomic operation sequence 00516 CPPAD_ASSERT_UNKNOWN( NumRes( UserOp ) == 0 ); 00517 CPPAD_ASSERT_UNKNOWN( NumArg( UserOp ) == 4 ); 00518 if( user_state == user_end ) 00519 { user_index = arg[0]; 00520 user_id = arg[1]; 00521 user_n = arg[2]; 00522 user_m = arg[3]; 00523 if(user_ix.size() < user_n) 00524 user_ix.resize(user_n); 00525 if(user_tx.size() < user_n * user_k1) 00526 { user_tx.resize(user_n * user_k1); 00527 user_px.resize(user_n * user_k1); 00528 } 00529 if(user_ty.size() < user_m * user_k1) 00530 { user_ty.resize(user_m * user_k1); 00531 user_py.resize(user_m * user_k1); 00532 } 00533 user_j = user_n; 00534 user_i = user_m; 00535 user_state = user_ret; 00536 } 00537 else 00538 { CPPAD_ASSERT_UNKNOWN( user_state == user_start ); 00539 CPPAD_ASSERT_UNKNOWN( user_index == size_t(arg[0]) ); 00540 CPPAD_ASSERT_UNKNOWN( user_id == size_t(arg[1]) ); 00541 CPPAD_ASSERT_UNKNOWN( user_n == size_t(arg[2]) ); 00542 CPPAD_ASSERT_UNKNOWN( user_m == size_t(arg[3]) ); 00543 user_state = user_end; 00544 00545 // call users function for this operation 00546 user_atomic<Base>::reverse(user_index, user_id, 00547 user_k, user_n, user_m, user_tx, user_ty, 00548 user_px, user_py 00549 ); 00550 for(j = 0; j < user_n; j++) if( user_ix[j] > 0 ) 00551 { for(ell = 0; ell < user_k1; ell++) 00552 Partial[user_ix[j] * K + ell] += 00553 user_px[j * user_k1 + ell]; 00554 } 00555 } 00556 break; 00557 00558 case UsrapOp: 00559 // parameter argument in an atomic operation sequence 00560 CPPAD_ASSERT_UNKNOWN( user_state == user_arg ); 00561 CPPAD_ASSERT_UNKNOWN( 0 < user_j && user_j <= user_n ); 00562 CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 ); 00563 CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par ); 00564 --user_j; 00565 user_ix[user_j] = 0; 00566 user_tx[user_j * user_k1 + 0] = parameter[ arg[0]]; 00567 for(ell = 1; ell < user_k1; ell++) 00568 user_tx[user_j * user_k1 + ell] = Base(0.); 00569 00570 if( user_j == 0 ) 00571 user_state = user_start; 00572 break; 00573 00574 case UsravOp: 00575 // variable argument in an atomic operation sequence 00576 CPPAD_ASSERT_UNKNOWN( user_state == user_arg ); 00577 CPPAD_ASSERT_UNKNOWN( 0 < user_j && user_j <= user_n ); 00578 CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 ); 00579 CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) <= i_var ); 00580 CPPAD_ASSERT_UNKNOWN( 0 < arg[0] ); 00581 --user_j; 00582 user_ix[user_j] = arg[0]; 00583 for(ell = 0; ell < user_k1; ell++) 00584 user_tx[user_j*user_k1 + ell] = Taylor[ arg[0] * J + ell]; 00585 if( user_j == 0 ) 00586 user_state = user_start; 00587 break; 00588 00589 case UsrrpOp: 00590 // parameter result in an atomic operation sequence 00591 CPPAD_ASSERT_UNKNOWN( user_state == user_ret ); 00592 CPPAD_ASSERT_UNKNOWN( 0 < user_i && user_i <= user_m ); 00593 CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 ); 00594 CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par ); 00595 --user_i; 00596 for(ell = 0; ell < user_k1; ell++) 00597 { user_py[user_i * user_k1 + ell] = Base(0.); 00598 user_ty[user_i * user_k1 + ell] = Base(0.); 00599 } 00600 user_ty[user_i * user_k1 + 0] = parameter[ arg[0] ]; 00601 if( user_i == 0 ) 00602 user_state = user_arg; 00603 break; 00604 00605 case UsrrvOp: 00606 // variable result in an atomic operation sequence 00607 CPPAD_ASSERT_UNKNOWN( user_state == user_ret ); 00608 CPPAD_ASSERT_UNKNOWN( 0 < user_i && user_i <= user_m ); 00609 --user_i; 00610 for(ell = 0; ell < user_k1; ell++) 00611 { user_py[user_i * user_k1 + ell] = 00612 Partial[i_var * K + ell]; 00613 user_ty[user_i * user_k1 + ell] = 00614 Taylor[i_var * J + ell]; 00615 } 00616 if( user_i == 0 ) 00617 user_state = user_arg; 00618 break; 00619 // ------------------------------------------------------------ 00620 00621 default: 00622 CPPAD_ASSERT_UNKNOWN(false); 00623 } 00624 } 00625 # if CPPAD_REVERSE_SWEEP_TRACE 00626 std::cout << std::endl; 00627 # endif 00628 // values corresponding to BeginOp 00629 CPPAD_ASSERT_UNKNOWN( i_op == 0 ); 00630 CPPAD_ASSERT_UNKNOWN( i_var == 0 ); 00631 } 00632 00633 /*! \} */ 00634 CPPAD_END_NAMESPACE 00635 00636 // preprocessor symbols that are local to this file 00637 # undef CPPAD_REVERSE_SWEEP_TRACE 00638 00639 # endif