CppAD: A C++ Algorithmic Differentiation Package
20130102
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_FORWARD0SWEEP_INCLUDED 00003 # define CPPAD_FORWARD0SWEEP_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 CPPAD_BEGIN_NAMESPACE 00017 /*! 00018 \defgroup forward0sweep_hpp forward0sweep.hpp 00019 \{ 00020 \file forward0sweep.hpp 00021 Compute zero order forward mode Taylor coefficients. 00022 */ 00023 00024 /*! 00025 \def CPPAD_FORWARD0SWEEP_TRACE 00026 This value is either zero or one. 00027 Zero is the normal operational value. 00028 If it is one, a trace of every forward0sweep computation is printed. 00029 (Note that forward0sweep is not used if CPPAD_USE_FORWARD0SWEEP is zero). 00030 */ 00031 # define CPPAD_FORWARD0SWEEP_TRACE 0 00032 00033 /*! 00034 Compute zero order forward mode Taylor coefficients. 00035 00036 \tparam Base 00037 The type used during the forward mode computations; i.e., the corresponding 00038 recording of operations used the type \c AD<Base>. 00039 00040 \param s_out 00041 Is the stream where output corresponding to \c PriOp operations will 00042 be written. 00043 00044 \param print 00045 If \a print is false, 00046 suppress the output that is otherwise generated by the c PriOp instructions. 00047 00048 \param n 00049 is the number of independent variables on the tape. 00050 00051 \param numvar 00052 is the total number of variables on the tape. 00053 This is also equal to the number of rows in the matrix \a Taylor; i.e., 00054 \a Rec->num_rec_var(). 00055 00056 \param Rec 00057 The information stored in \a Rec 00058 is a recording of the operations corresponding to the function 00059 \f[ 00060 F : {\bf R}^n \rightarrow {\bf R}^m 00061 \f] 00062 where \f$ n \f$ is the number of independent variables and 00063 \f$ m \f$ is the number of dependent variables. 00064 \n 00065 \n 00066 The object \a Rec is effectly constant. 00067 There are two exceptions to this. 00068 The first exception is that while palying back the tape 00069 the object \a Rec holds information about the current location 00070 with in the tape and this changes during palyback. 00071 The second exception is the fact that the 00072 zero order ( \a d = 0 ) versions of the VecAD operators LdpOp and LdvOp 00073 modify the corresponding \a op_arg values returned by 00074 \ref player::next_forward and \ref player::next_reverse; see the 00075 \link load_op.hpp LdpOp and LdvOp \endlink operations. 00076 00077 \param J 00078 Is the number of columns in the coefficient matrix \a Taylor. 00079 This must be greater than or equal one. 00080 00081 \param Taylor 00082 \b Input: For j = 1 , ... , \a n, 00083 \a Taylor [ j * J + 0 ] 00084 variable with index i on the tape 00085 (independent variable with index (j-1) in the independent variable vector). 00086 \n 00087 \n 00088 \b Output: For i = \a n + 1, ... , \a numvar - 1, 00089 \a Taylor [ i * J + 0 ] 00090 is the zero order Taylor coefficient for the variable with 00091 index i on the tape. 00092 00093 \a return 00094 The return value is equal to the number of ComOp operations 00095 that have a different result from when the information in 00096 \a Rec was recorded. 00097 (Note that if NDEBUG is true, there are no ComOp operations 00098 in Rec and hence this return value is always zero.) 00099 */ 00100 00101 template <class Base> 00102 size_t forward0sweep( 00103 std::ostream& s_out, 00104 bool print, 00105 size_t n, 00106 size_t numvar, 00107 player<Base> *Rec, 00108 size_t J, 00109 Base *Taylor 00110 ) 00111 { CPPAD_ASSERT_UNKNOWN( J >= 1 ); 00112 00113 // op code for current instruction 00114 OpCode op; 00115 00116 // index for current instruction 00117 size_t i_op; 00118 00119 // next variables 00120 size_t i_var; 00121 00122 // constant and non-constant version of the operation argument indices 00123 addr_t* non_const_arg; 00124 const addr_t* arg = 0; 00125 00126 // initialize the comparision operator (ComOp) counter 00127 size_t compareCount = 0; 00128 00129 // This is an order zero calculation, initialize vector indices 00130 pod_vector<size_t> VectorInd; // address for each element 00131 pod_vector<bool> VectorVar; // is element a variable 00132 size_t i = Rec->num_rec_vecad_ind(); 00133 if( i > 0 ) 00134 { VectorInd.extend(i); 00135 VectorVar.extend(i); 00136 while(i--) 00137 { VectorInd[i] = Rec->GetVecInd(i); 00138 VectorVar[i] = false; 00139 } 00140 } 00141 00142 // work space used by UserOp. 00143 const size_t user_k = 0; // order of this forward mode calculation 00144 vector<Base> user_tx; // argument vector Taylor coefficients 00145 vector<Base> user_ty; // result vector Taylor coefficients 00146 size_t user_index = 0; // indentifier for this user_atomic operation 00147 size_t user_id = 0; // user identifier for this call to operator 00148 size_t user_i = 0; // index in result vector 00149 size_t user_j = 0; // index in argument vector 00150 size_t user_m = 0; // size of result vector 00151 size_t user_n = 0; // size of arugment vector 00152 // next expected operator in a UserOp sequence 00153 enum { user_start, user_arg, user_ret, user_end } user_state = user_start; 00154 00155 // check numvar argument 00156 CPPAD_ASSERT_UNKNOWN( Rec->num_rec_var() == numvar ); 00157 00158 // length of the parameter vector (used by CppAD assert macros) 00159 const size_t num_par = Rec->num_rec_par(); 00160 00161 // length of the text vector (used by CppAD assert macros) 00162 const size_t num_text = Rec->num_rec_text(); 00163 00164 // pointer to the beginning of the parameter vector 00165 const Base* parameter = CPPAD_NULL; 00166 if( num_par > 0 ) 00167 parameter = Rec->GetPar(); 00168 00169 // pointer to the beginning of the text vector 00170 const char* text = CPPAD_NULL; 00171 if( num_text > 0 ) 00172 text = Rec->GetTxt(0); 00173 00174 // skip the BeginOp at the beginning of the recording 00175 Rec->start_forward(op, arg, i_op, i_var); 00176 CPPAD_ASSERT_UNKNOWN( op == BeginOp ); 00177 # if CPPAD_FORWARD0SWEEP_TRACE 00178 std::cout << std::endl; 00179 # endif 00180 while(op != EndOp) 00181 { 00182 // this op 00183 Rec->next_forward(op, arg, i_op, i_var); 00184 # ifndef NDEBUG 00185 if( i_op <= n ) 00186 { CPPAD_ASSERT_UNKNOWN((op == InvOp) | (op == BeginOp)); 00187 } 00188 else CPPAD_ASSERT_UNKNOWN((op != InvOp) & (op != BeginOp)); 00189 # endif 00190 00191 // action to take depends on the case 00192 switch( op ) 00193 { 00194 case AbsOp: 00195 forward_abs_op_0(i_var, arg[0], J, Taylor); 00196 break; 00197 // ------------------------------------------------- 00198 00199 case AddvvOp: 00200 forward_addvv_op_0(i_var, arg, parameter, J, Taylor); 00201 break; 00202 // ------------------------------------------------- 00203 00204 case AddpvOp: 00205 CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par ); 00206 forward_addpv_op_0(i_var, arg, parameter, J, Taylor); 00207 break; 00208 // ------------------------------------------------- 00209 00210 case AcosOp: 00211 // sqrt(1 - x * x), acos(x) 00212 CPPAD_ASSERT_UNKNOWN( i_var < numvar ); 00213 forward_acos_op_0(i_var, arg[0], J, Taylor); 00214 break; 00215 // ------------------------------------------------- 00216 00217 case AsinOp: 00218 // sqrt(1 - x * x), asin(x) 00219 CPPAD_ASSERT_UNKNOWN( i_var < numvar ); 00220 forward_asin_op_0(i_var, arg[0], J, Taylor); 00221 break; 00222 // ------------------------------------------------- 00223 00224 case AtanOp: 00225 // 1 + x * x, atan(x) 00226 CPPAD_ASSERT_UNKNOWN( i_var < numvar ); 00227 forward_atan_op_0(i_var, arg[0], J, Taylor); 00228 break; 00229 // ------------------------------------------------- 00230 00231 case CSumOp: 00232 // CSumOp has a variable number of arguments and 00233 // next_forward thinks it one has one argument. 00234 // we must inform next_forward of this special case. 00235 Rec->forward_csum(op, arg, i_op, i_var); 00236 forward_csum_op( 00237 0, i_var, arg, num_par, parameter, J, Taylor 00238 ); 00239 break; 00240 00241 // ------------------------------------------------- 00242 case CExpOp: 00243 // Use the general case with d == 0 00244 // (could create an optimzied verison for this case) 00245 forward_cond_op_0( 00246 i_var, arg, num_par, parameter, J, Taylor 00247 ); 00248 break; 00249 // --------------------------------------------------- 00250 case ComOp: 00251 forward_comp_op_0( 00252 compareCount, arg, num_par, parameter, J, Taylor 00253 ); 00254 break; 00255 // --------------------------------------------------- 00256 00257 case CosOp: 00258 // sin(x), cos(x) 00259 CPPAD_ASSERT_UNKNOWN( i_var < numvar ); 00260 forward_cos_op_0(i_var, arg[0], J, Taylor); 00261 break; 00262 // --------------------------------------------------- 00263 00264 case CoshOp: 00265 // sinh(x), cosh(x) 00266 CPPAD_ASSERT_UNKNOWN( i_var < numvar ); 00267 forward_cosh_op_0(i_var, arg[0], J, Taylor); 00268 break; 00269 // ------------------------------------------------- 00270 00271 case DisOp: 00272 forward_dis_op_0(i_var, arg, J, Taylor); 00273 break; 00274 // ------------------------------------------------- 00275 00276 case DivvvOp: 00277 forward_divvv_op_0(i_var, arg, parameter, J, Taylor); 00278 break; 00279 // ------------------------------------------------- 00280 00281 case DivpvOp: 00282 CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par ); 00283 forward_divpv_op_0(i_var, arg, parameter, J, Taylor); 00284 break; 00285 // ------------------------------------------------- 00286 00287 case DivvpOp: 00288 CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < num_par ); 00289 forward_divvp_op_0(i_var, arg, parameter, J, Taylor); 00290 break; 00291 // ------------------------------------------------- 00292 00293 case EndOp: 00294 CPPAD_ASSERT_NARG_NRES(op, 0, 0); 00295 break; 00296 // ------------------------------------------------- 00297 00298 case ExpOp: 00299 forward_exp_op_0(i_var, arg[0], J, Taylor); 00300 break; 00301 // ------------------------------------------------- 00302 00303 case InvOp: 00304 break; 00305 // ------------------------------------------------- 00306 00307 case LdpOp: 00308 CPPAD_ASSERT_UNKNOWN( VectorInd.size() != 0 ); 00309 CPPAD_ASSERT_UNKNOWN( VectorVar.size() != 0 ); 00310 non_const_arg = Rec->forward_non_const_arg(); 00311 forward_load_p_op_0( 00312 i_var, 00313 non_const_arg, 00314 num_par, 00315 parameter, 00316 J, 00317 Taylor, 00318 Rec->num_rec_vecad_ind(), 00319 VectorVar.data(), 00320 VectorInd.data() 00321 ); 00322 break; 00323 // ------------------------------------------------- 00324 00325 case LdvOp: 00326 CPPAD_ASSERT_UNKNOWN( VectorInd.size() != 0 ); 00327 CPPAD_ASSERT_UNKNOWN( VectorVar.size() != 0 ); 00328 non_const_arg = Rec->forward_non_const_arg(); 00329 forward_load_v_op_0( 00330 i_var, 00331 non_const_arg, 00332 num_par, 00333 parameter, 00334 J, 00335 Taylor, 00336 Rec->num_rec_vecad_ind(), 00337 VectorVar.data(), 00338 VectorInd.data() 00339 ); 00340 break; 00341 // ------------------------------------------------- 00342 00343 case LogOp: 00344 forward_log_op_0(i_var, arg[0], J, Taylor); 00345 break; 00346 // ------------------------------------------------- 00347 00348 case MulvvOp: 00349 forward_mulvv_op_0(i_var, arg, parameter, J, Taylor); 00350 break; 00351 // ------------------------------------------------- 00352 00353 case MulpvOp: 00354 CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par ); 00355 forward_mulpv_op_0(i_var, arg, parameter, J, Taylor); 00356 break; 00357 // ------------------------------------------------- 00358 00359 case ParOp: 00360 forward_par_op_0( 00361 i_var, arg, num_par, parameter, J, Taylor 00362 ); 00363 break; 00364 // ------------------------------------------------- 00365 00366 case PowvpOp: 00367 CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < num_par ); 00368 forward_powvp_op_0(i_var, arg, parameter, J, Taylor); 00369 break; 00370 // ------------------------------------------------- 00371 00372 case PowpvOp: 00373 CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par ); 00374 forward_powpv_op_0(i_var, arg, parameter, J, Taylor); 00375 break; 00376 // ------------------------------------------------- 00377 00378 case PowvvOp: 00379 forward_powvv_op_0(i_var, arg, parameter, J, Taylor); 00380 break; 00381 // ------------------------------------------------- 00382 00383 case PriOp: 00384 if( print ) forward_pri_0(s_out, 00385 i_var, arg, num_text, text, num_par, parameter, J, Taylor 00386 ); 00387 break; 00388 // ------------------------------------------------- 00389 00390 case SignOp: 00391 // cos(x), sin(x) 00392 CPPAD_ASSERT_UNKNOWN( i_var < numvar ); 00393 forward_sign_op_0(i_var, arg[0], J, Taylor); 00394 break; 00395 // ------------------------------------------------- 00396 00397 case SinOp: 00398 // cos(x), sin(x) 00399 CPPAD_ASSERT_UNKNOWN( i_var < numvar ); 00400 forward_sin_op_0(i_var, arg[0], J, Taylor); 00401 break; 00402 // ------------------------------------------------- 00403 00404 case SinhOp: 00405 // cosh(x), sinh(x) 00406 CPPAD_ASSERT_UNKNOWN( i_var < numvar ); 00407 forward_sinh_op_0(i_var, arg[0], J, Taylor); 00408 break; 00409 // ------------------------------------------------- 00410 00411 case SqrtOp: 00412 forward_sqrt_op_0(i_var, arg[0], J, Taylor); 00413 break; 00414 // ------------------------------------------------- 00415 00416 case StppOp: 00417 forward_store_pp_op_0( 00418 i_var, 00419 arg, 00420 num_par, 00421 J, 00422 Taylor, 00423 Rec->num_rec_vecad_ind(), 00424 VectorVar.data(), 00425 VectorInd.data() 00426 ); 00427 break; 00428 // ------------------------------------------------- 00429 00430 case StpvOp: 00431 forward_store_pv_op_0( 00432 i_var, 00433 arg, 00434 num_par, 00435 J, 00436 Taylor, 00437 Rec->num_rec_vecad_ind(), 00438 VectorVar.data(), 00439 VectorInd.data() 00440 ); 00441 break; 00442 // ------------------------------------------------- 00443 00444 case StvpOp: 00445 forward_store_vp_op_0( 00446 i_var, 00447 arg, 00448 num_par, 00449 J, 00450 Taylor, 00451 Rec->num_rec_vecad_ind(), 00452 VectorVar.data(), 00453 VectorInd.data() 00454 ); 00455 break; 00456 // ------------------------------------------------- 00457 00458 case StvvOp: 00459 forward_store_vv_op_0( 00460 i_var, 00461 arg, 00462 num_par, 00463 J, 00464 Taylor, 00465 Rec->num_rec_vecad_ind(), 00466 VectorVar.data(), 00467 VectorInd.data() 00468 ); 00469 break; 00470 // ------------------------------------------------- 00471 00472 case SubvvOp: 00473 forward_subvv_op_0(i_var, arg, parameter, J, Taylor); 00474 break; 00475 // ------------------------------------------------- 00476 00477 case SubpvOp: 00478 CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par ); 00479 forward_subpv_op_0(i_var, arg, parameter, J, Taylor); 00480 break; 00481 // ------------------------------------------------- 00482 00483 case SubvpOp: 00484 CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < num_par ); 00485 forward_subvp_op_0(i_var, arg, parameter, J, Taylor); 00486 break; 00487 // ------------------------------------------------- 00488 00489 case TanOp: 00490 // tan(x)^2, tan(x) 00491 CPPAD_ASSERT_UNKNOWN( i_var < numvar ); 00492 forward_tan_op_0(i_var, arg[0], J, Taylor); 00493 break; 00494 // ------------------------------------------------- 00495 00496 case TanhOp: 00497 // tanh(x)^2, tanh(x) 00498 CPPAD_ASSERT_UNKNOWN( i_var < numvar ); 00499 forward_tanh_op_0(i_var, arg[0], J, Taylor); 00500 break; 00501 // ------------------------------------------------- 00502 00503 case UserOp: 00504 // start or end an atomic operation sequence 00505 CPPAD_ASSERT_UNKNOWN( NumRes( UserOp ) == 0 ); 00506 CPPAD_ASSERT_UNKNOWN( NumArg( UserOp ) == 4 ); 00507 if( user_state == user_start ) 00508 { user_index = arg[0]; 00509 user_id = arg[1]; 00510 user_n = arg[2]; 00511 user_m = arg[3]; 00512 if(user_tx.size() < user_n) 00513 user_tx.resize(user_n); 00514 if(user_ty.size() < user_m) 00515 user_ty.resize(user_m); 00516 user_j = 0; 00517 user_i = 0; 00518 user_state = user_arg; 00519 } 00520 else 00521 { CPPAD_ASSERT_UNKNOWN( user_state == user_end ); 00522 CPPAD_ASSERT_UNKNOWN( user_index == size_t(arg[0]) ); 00523 CPPAD_ASSERT_UNKNOWN( user_id == size_t(arg[1]) ); 00524 CPPAD_ASSERT_UNKNOWN( user_n == size_t(arg[2]) ); 00525 CPPAD_ASSERT_UNKNOWN( user_m == size_t(arg[3]) ); 00526 user_state = user_start; 00527 } 00528 break; 00529 00530 case UsrapOp: 00531 // parameter argument in an atomic operation sequence 00532 CPPAD_ASSERT_UNKNOWN( user_state == user_arg ); 00533 CPPAD_ASSERT_UNKNOWN( user_j < user_n ); 00534 CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par ); 00535 user_tx[user_j++] = parameter[ arg[0] ]; 00536 if( user_j == user_n ) 00537 { // call users function for this operation 00538 user_atomic<Base>::forward(user_index, user_id, 00539 user_k, user_n, user_m, user_tx, user_ty 00540 ); 00541 user_state = user_ret; 00542 } 00543 break; 00544 00545 case UsravOp: 00546 // variable argument in an atomic operation sequence 00547 CPPAD_ASSERT_UNKNOWN( user_state == user_arg ); 00548 CPPAD_ASSERT_UNKNOWN( user_j < user_n ); 00549 CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) <= i_var ); 00550 user_tx[user_j++] = Taylor[ arg[0] * J + 0 ]; 00551 if( user_j == user_n ) 00552 { // call users function for this operation 00553 user_atomic<Base>::forward(user_index, user_id, 00554 user_k, user_n, user_m, user_tx, user_ty 00555 ); 00556 user_state = user_ret; 00557 } 00558 break; 00559 00560 case UsrrpOp: 00561 // parameter result in an atomic operation sequence 00562 CPPAD_ASSERT_UNKNOWN( user_state == user_ret ); 00563 CPPAD_ASSERT_UNKNOWN( user_i < user_m ); 00564 user_i++; 00565 if( user_i == user_m ) 00566 user_state = user_end; 00567 break; 00568 00569 case UsrrvOp: 00570 // variable result in an atomic operation sequence 00571 CPPAD_ASSERT_UNKNOWN( user_state == user_ret ); 00572 CPPAD_ASSERT_UNKNOWN( user_i < user_m ); 00573 Taylor[ i_var * J + 0 ] = user_ty[user_i++]; 00574 if( user_i == user_m ) 00575 user_state = user_end; 00576 break; 00577 // ------------------------------------------------- 00578 00579 default: 00580 CPPAD_ASSERT_UNKNOWN(0); 00581 } 00582 # if CPPAD_FORWARD0SWEEP_TRACE 00583 size_t d = 0; 00584 size_t i_tmp = i_var; 00585 Base* Z_tmp = Taylor + i_var * J; 00586 00587 printOp( 00588 std::cout, 00589 Rec, 00590 i_tmp, 00591 op, 00592 arg, 00593 d + 1, 00594 Z_tmp, 00595 0, 00596 (Base *) CPPAD_NULL 00597 ); 00598 } 00599 std::cout << std::endl; 00600 # else 00601 } 00602 # endif 00603 CPPAD_ASSERT_UNKNOWN( user_state == user_start ); 00604 CPPAD_ASSERT_UNKNOWN( i_var + 1 == Rec->num_rec_var() ); 00605 00606 return compareCount; 00607 } 00608 00609 /*! \} */ 00610 CPPAD_END_NAMESPACE 00611 00612 // preprocessor symbols that are local to this file 00613 # undef CPPAD_FORWARD0SWEEP_TRACE 00614 00615 # endif