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