CppAD: A C++ Algorithmic Differentiation Package 20110419
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_REV_JAC_SWEEP_INCLUDED 00003 # define CPPAD_REV_JAC_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 CPPAD_BEGIN_NAMESPACE 00017 /*! 00018 \file rev_jac_sweep.hpp 00019 Compute Reverse mode Jacobian sparsity patterns. 00020 */ 00021 00022 /*! 00023 \def CPPAD_REV_JAC_SWEEP_TRACE 00024 This value is either zero or one. 00025 Zero is the normal operational value. 00026 If it is one, a trace of every rev_jac_sweep computation is printed. 00027 */ 00028 # define CPPAD_REV_JAC_SWEEP_TRACE 0 00029 00030 /*! 00031 Given the sparsity pattern for the dependent variables, 00032 RevJacSweep computes the sparsity pattern for all the independent variables. 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 \tparam Vector_set 00040 is the type used for vectors of sets. It can be either 00041 \c sparse_pack or \c sparse_set. 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; i.e., 00048 \a play->num_rec_var(). 00049 This is also the number of rows in the entire sparsity pattern \a RevJac. 00050 00051 \param play 00052 The information stored in \a play 00053 is a recording of the operations corresponding to a 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 00058 and \f$ m \f$ is the number of dependent variables. 00059 The object \a play is effectly constant. 00060 It is not declared const because while playing back the tape 00061 the object \a play holds information about the currentl location 00062 with in the tape and this changes during playback. 00063 00064 \param var_sparsity 00065 For i = 0 , ... , \a numvar - 1, 00066 (all the variables on the tape) 00067 the forward Jacobian sparsity pattern for variable i 00068 corresponds to the set with index i in \a var_sparsity. 00069 \b 00070 \b 00071 \b Input: 00072 For i = 0 , ... , \a numvar - 1, 00073 the forward Jacobian sparsity pattern for variable i is an input 00074 if i corresponds to a dependent variable. 00075 Otherwise the sparsity patten is empty. 00076 \n 00077 \n 00078 \b Output: For j = 1 , ... , \a n, 00079 the sparsity pattern for the dependent variable with index (j-1) 00080 is given by the set with index index j in \a var_sparsity. 00081 */ 00082 00083 template <class Base, class Vector_set> 00084 void RevJacSweep( 00085 size_t n, 00086 size_t numvar, 00087 player<Base> *play, 00088 Vector_set& var_sparsity 00089 ) 00090 { 00091 OpCode op; 00092 size_t i_op; 00093 size_t i_var; 00094 00095 const size_t *arg = 0; 00096 00097 size_t i, j, k; 00098 00099 // length of the parameter vector (used by CppAD assert macros) 00100 const size_t num_par = play->num_rec_par(); 00101 00102 // check numvar argument 00103 CPPAD_ASSERT_UNKNOWN( numvar > 0 ); 00104 CPPAD_ASSERT_UNKNOWN( play->num_rec_var() == numvar ); 00105 CPPAD_ASSERT_UNKNOWN( var_sparsity.n_set() == numvar ); 00106 00107 // upper limit (exclusive) for elements in the set 00108 size_t limit = var_sparsity.end(); 00109 00110 // vecad_sparsity contains a sparsity pattern for each VecAD object. 00111 // vecad_ind maps a VecAD index (beginning of the VecAD object) 00112 // to the index of the corresponding set in vecad_sparsity. 00113 size_t num_vecad_ind = play->num_rec_vecad_ind(); 00114 size_t num_vecad_vec = play->num_rec_vecad_vec(); 00115 Vector_set vecad_sparsity; 00116 vecad_sparsity.resize(num_vecad_vec, limit); 00117 size_t* vecad_ind = CPPAD_NULL; 00118 if( num_vecad_vec > 0 ) 00119 { size_t length; 00120 vecad_ind = CPPAD_TRACK_NEW_VEC(num_vecad_ind, vecad_ind); 00121 j = 0; 00122 for(i = 0; i < num_vecad_vec; i++) 00123 { // length of this VecAD 00124 length = play->GetVecInd(j); 00125 // set to proper index for this VecAD 00126 vecad_ind[j] = i; 00127 for(k = 1; k <= length; k++) 00128 vecad_ind[j+k] = num_vecad_vec; // invalid index 00129 // start of next VecAD 00130 j += length + 1; 00131 } 00132 CPPAD_ASSERT_UNKNOWN( j == play->num_rec_vecad_ind() ); 00133 } 00134 00135 // work space used by UserOp. 00136 typedef std::set<size_t> size_set; 00137 const size_t user_q = limit; // maximum element plus one 00138 size_set::iterator set_itr; // iterator for a standard set 00139 size_set::iterator set_end; // end of iterator sequence 00140 vector< size_set > user_r; // sparsity pattern for the argument x 00141 vector< size_set > user_s; // sparisty pattern for the result y 00142 size_t user_index = 0; // indentifier for this user_atomic operation 00143 size_t user_id = 0; // user identifier for this call to operator 00144 size_t user_i = 0; // index in result vector 00145 size_t user_j = 0; // index in argument vector 00146 size_t user_m = 0; // size of result vector 00147 size_t user_n = 0; // size of arugment vector 00148 // next expected operator in a UserOp sequence 00149 enum { user_start, user_arg, user_ret, user_end } user_state = user_end; 00150 00151 // Initialize 00152 play->start_reverse(op, arg, i_op, i_var); 00153 CPPAD_ASSERT_UNKNOWN( op == EndOp ); 00154 # if CPPAD_REV_JAC_SWEEP_TRACE 00155 std::cout << std::endl; 00156 CppAD::vector<bool> z_value(limit); 00157 # endif 00158 while(op != BeginOp ) 00159 { 00160 // next op 00161 play->next_reverse(op, arg, i_op, i_var); 00162 # ifndef NDEBUG 00163 if( i_op <= n ) 00164 { CPPAD_ASSERT_UNKNOWN((op == InvOp) | (op == BeginOp)); 00165 } 00166 else CPPAD_ASSERT_UNKNOWN((op != InvOp) & (op != BeginOp)); 00167 # endif 00168 00169 // rest of information depends on the case 00170 switch( op ) 00171 { 00172 case AbsOp: 00173 CPPAD_ASSERT_NARG_NRES(op, 1, 1); 00174 reverse_sparse_jacobian_unary_op( 00175 i_var, arg[0], var_sparsity 00176 ); 00177 break; 00178 // ------------------------------------------------- 00179 00180 case AddvvOp: 00181 CPPAD_ASSERT_NARG_NRES(op, 2, 1); 00182 reverse_sparse_jacobian_binary_op( 00183 i_var, arg, var_sparsity 00184 ); 00185 break; 00186 // ------------------------------------------------- 00187 00188 case AddpvOp: 00189 CPPAD_ASSERT_NARG_NRES(op, 2, 1); 00190 reverse_sparse_jacobian_unary_op( 00191 i_var, arg[1], var_sparsity 00192 ); 00193 break; 00194 // ------------------------------------------------- 00195 00196 case AcosOp: 00197 // acos(x) and sqrt(1 - x * x) are computed in pairs 00198 // but i_var + 1 should only be used here 00199 CPPAD_ASSERT_NARG_NRES(op, 1, 2); 00200 reverse_sparse_jacobian_unary_op( 00201 i_var, arg[0], var_sparsity 00202 ); 00203 break; 00204 // ------------------------------------------------- 00205 00206 case AsinOp: 00207 // asin(x) and sqrt(1 - x * x) are computed in pairs 00208 // but i_var + 1 should only be used here 00209 CPPAD_ASSERT_NARG_NRES(op, 1, 2); 00210 reverse_sparse_jacobian_unary_op( 00211 i_var, arg[0], var_sparsity 00212 ); 00213 break; 00214 // ------------------------------------------------- 00215 00216 case AtanOp: 00217 // atan(x) and 1 + x * x must be computed in pairs 00218 // but i_var + 1 should only be used here 00219 CPPAD_ASSERT_NARG_NRES(op, 1, 2); 00220 reverse_sparse_jacobian_unary_op( 00221 i_var, arg[0], var_sparsity 00222 ); 00223 break; 00224 // ------------------------------------------------- 00225 00226 case BeginOp: 00227 CPPAD_ASSERT_NARG_NRES(op, 0, 1); 00228 break; 00229 // ------------------------------------------------- 00230 00231 case CSumOp: 00232 // CSumOp has a variable number of arguments and 00233 // next_reverse thinks it one has one argument. 00234 // We must inform next_reverse of this special case. 00235 play->reverse_csum(op, arg, i_op, i_var); 00236 reverse_sparse_jacobian_csum_op( 00237 i_var, arg, var_sparsity 00238 ); 00239 break; 00240 // ------------------------------------------------- 00241 00242 case CExpOp: 00243 reverse_sparse_jacobian_cond_op( 00244 i_var, arg, num_par, var_sparsity 00245 ); 00246 break; 00247 // --------------------------------------------------- 00248 00249 case ComOp: 00250 CPPAD_ASSERT_NARG_NRES(op, 4, 0); 00251 CPPAD_ASSERT_UNKNOWN( arg[1] > 1 ); 00252 break; 00253 // -------------------------------------------------- 00254 00255 case CosOp: 00256 // cosine and sine must come in pairs 00257 // but i_var should only be used here 00258 CPPAD_ASSERT_NARG_NRES(op, 1, 2); 00259 reverse_sparse_jacobian_unary_op( 00260 i_var, arg[0], var_sparsity 00261 ); 00262 break; 00263 // --------------------------------------------------- 00264 00265 case CoshOp: 00266 // hyperbolic cosine and sine must come in pairs 00267 // but i_var should only be used here 00268 CPPAD_ASSERT_NARG_NRES(op, 1, 2); 00269 reverse_sparse_jacobian_unary_op( 00270 i_var, arg[0], var_sparsity 00271 ); 00272 break; 00273 // ------------------------------------------------- 00274 00275 case DisOp: 00276 CPPAD_ASSERT_NARG_NRES(op, 2, 1); 00277 00278 break; 00279 // ------------------------------------------------- 00280 00281 case DivvvOp: 00282 CPPAD_ASSERT_NARG_NRES(op, 2, 1); 00283 reverse_sparse_jacobian_binary_op( 00284 i_var, arg, var_sparsity 00285 ); 00286 break; 00287 // ------------------------------------------------- 00288 00289 case DivpvOp: 00290 CPPAD_ASSERT_NARG_NRES(op, 2, 1); 00291 reverse_sparse_jacobian_unary_op( 00292 i_var, arg[1], var_sparsity 00293 ); 00294 break; 00295 // ------------------------------------------------- 00296 00297 case DivvpOp: 00298 CPPAD_ASSERT_NARG_NRES(op, 2, 1); 00299 reverse_sparse_jacobian_unary_op( 00300 i_var, arg[0], var_sparsity 00301 ); 00302 break; 00303 // ------------------------------------------------- 00304 00305 case ExpOp: 00306 CPPAD_ASSERT_NARG_NRES(op, 1, 1); 00307 reverse_sparse_jacobian_unary_op( 00308 i_var, arg[0], var_sparsity 00309 ); 00310 break; 00311 // ------------------------------------------------- 00312 00313 case InvOp: 00314 CPPAD_ASSERT_NARG_NRES(op, 0, 1); 00315 break; 00316 // ------------------------------------------------- 00317 00318 case LdpOp: 00319 reverse_sparse_jacobian_load_op( 00320 op, 00321 i_var, 00322 arg, 00323 num_vecad_ind, 00324 vecad_ind, 00325 var_sparsity, 00326 vecad_sparsity 00327 ); 00328 break; 00329 // ------------------------------------------------- 00330 00331 case LdvOp: 00332 reverse_sparse_jacobian_load_op( 00333 op, 00334 i_var, 00335 arg, 00336 num_vecad_ind, 00337 vecad_ind, 00338 var_sparsity, 00339 vecad_sparsity 00340 ); 00341 break; 00342 // ------------------------------------------------- 00343 00344 case LogOp: 00345 CPPAD_ASSERT_NARG_NRES(op, 1, 1); 00346 reverse_sparse_jacobian_unary_op( 00347 i_var, arg[0], var_sparsity 00348 ); 00349 break; 00350 // ------------------------------------------------- 00351 00352 case MulvvOp: 00353 CPPAD_ASSERT_NARG_NRES(op, 2, 1); 00354 reverse_sparse_jacobian_binary_op( 00355 i_var, arg, var_sparsity 00356 ); 00357 break; 00358 // ------------------------------------------------- 00359 00360 case MulpvOp: 00361 CPPAD_ASSERT_NARG_NRES(op, 2, 1); 00362 reverse_sparse_jacobian_unary_op( 00363 i_var, arg[1], var_sparsity 00364 ); 00365 break; 00366 // ------------------------------------------------- 00367 00368 case ParOp: 00369 CPPAD_ASSERT_NARG_NRES(op, 1, 1); 00370 00371 break; 00372 // ------------------------------------------------- 00373 00374 case PowvpOp: 00375 reverse_sparse_jacobian_unary_op( 00376 i_var, arg[0], var_sparsity 00377 ); 00378 break; 00379 // ------------------------------------------------- 00380 00381 case PowpvOp: 00382 CPPAD_ASSERT_NARG_NRES(op, 2, 3); 00383 reverse_sparse_jacobian_unary_op( 00384 i_var, arg[1], var_sparsity 00385 ); 00386 break; 00387 // ------------------------------------------------- 00388 00389 case PowvvOp: 00390 CPPAD_ASSERT_NARG_NRES(op, 2, 3); 00391 reverse_sparse_jacobian_binary_op( 00392 i_var, arg, var_sparsity 00393 ); 00394 break; 00395 // ------------------------------------------------- 00396 00397 case PripOp: 00398 CPPAD_ASSERT_NARG_NRES(op, 2, 0); 00399 00400 break; 00401 // ------------------------------------------------- 00402 00403 case PrivOp: 00404 CPPAD_ASSERT_NARG_NRES(op, 2, 0); 00405 break; 00406 // ------------------------------------------------- 00407 00408 case SinOp: 00409 // sine and cosine must come in pairs 00410 // but i_var should only be used here 00411 CPPAD_ASSERT_NARG_NRES(op, 1, 2); 00412 reverse_sparse_jacobian_unary_op( 00413 i_var, arg[0], var_sparsity 00414 ); 00415 break; 00416 // ------------------------------------------------- 00417 00418 case SinhOp: 00419 // hyperbolic sine and cosine must come in pairs 00420 // but i_var should only be used here 00421 CPPAD_ASSERT_NARG_NRES(op, 1, 2); 00422 reverse_sparse_jacobian_unary_op( 00423 i_var, arg[0], var_sparsity 00424 ); 00425 break; 00426 // ------------------------------------------------- 00427 00428 case SqrtOp: 00429 CPPAD_ASSERT_NARG_NRES(op, 1, 1); 00430 reverse_sparse_jacobian_unary_op( 00431 i_var, arg[0], var_sparsity 00432 ); 00433 break; 00434 // ------------------------------------------------- 00435 00436 case StppOp: 00437 // sparsity cannot proagate through a parameter 00438 CPPAD_ASSERT_NARG_NRES(op, 3, 0); 00439 break; 00440 // ------------------------------------------------- 00441 00442 case StpvOp: 00443 reverse_sparse_jacobian_store_op( 00444 op, 00445 arg, 00446 num_vecad_ind, 00447 vecad_ind, 00448 var_sparsity, 00449 vecad_sparsity 00450 ); 00451 break; 00452 // ------------------------------------------------- 00453 00454 case StvpOp: 00455 CPPAD_ASSERT_NARG_NRES(op, 3, 0); 00456 break; 00457 // ------------------------------------------------- 00458 00459 case StvvOp: 00460 reverse_sparse_jacobian_store_op( 00461 op, 00462 arg, 00463 num_vecad_ind, 00464 vecad_ind, 00465 var_sparsity, 00466 vecad_sparsity 00467 ); 00468 break; 00469 // ------------------------------------------------- 00470 00471 case SubvvOp: 00472 CPPAD_ASSERT_NARG_NRES(op, 2, 1); 00473 reverse_sparse_jacobian_binary_op( 00474 i_var, arg, var_sparsity 00475 ); 00476 break; 00477 // ------------------------------------------------- 00478 00479 case SubpvOp: 00480 CPPAD_ASSERT_NARG_NRES(op, 2, 1); 00481 reverse_sparse_jacobian_unary_op( 00482 i_var, arg[1], var_sparsity 00483 ); 00484 break; 00485 // ------------------------------------------------- 00486 00487 case SubvpOp: 00488 CPPAD_ASSERT_NARG_NRES(op, 2, 1); 00489 reverse_sparse_jacobian_unary_op( 00490 i_var, arg[0], var_sparsity 00491 ); 00492 break; 00493 // ------------------------------------------------- 00494 00495 case UserOp: 00496 // start or end atomic operation sequence 00497 CPPAD_ASSERT_UNKNOWN( NumRes( UserOp ) == 0 ); 00498 CPPAD_ASSERT_UNKNOWN( NumArg( UserOp ) == 4 ); 00499 if( user_state == user_end ) 00500 { user_index = arg[0]; 00501 user_id = arg[1]; 00502 user_n = arg[2]; 00503 user_m = arg[3]; 00504 if(user_r.size() < user_n ) 00505 user_r.resize(user_n); 00506 if(user_s.size() < user_m ) 00507 user_s.resize(user_m); 00508 user_j = user_n; 00509 user_i = user_m; 00510 user_state = user_ret; 00511 } 00512 else 00513 { CPPAD_ASSERT_UNKNOWN( user_state == user_start ); 00514 CPPAD_ASSERT_UNKNOWN( user_index == arg[0] ); 00515 CPPAD_ASSERT_UNKNOWN( user_id == arg[1] ); 00516 CPPAD_ASSERT_UNKNOWN( user_n == arg[2] ); 00517 CPPAD_ASSERT_UNKNOWN( user_m == arg[3] ); 00518 user_state = user_end; 00519 } 00520 break; 00521 00522 case UsrapOp: 00523 // parameter argument in an atomic operation sequence 00524 CPPAD_ASSERT_UNKNOWN( user_state == user_arg ); 00525 CPPAD_ASSERT_UNKNOWN( 0 < user_j && user_j <= user_n ); 00526 CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 ); 00527 CPPAD_ASSERT_UNKNOWN( arg[0] < num_par ); 00528 --user_j; 00529 if( user_j == 0 ) 00530 user_state = user_start; 00531 break; 00532 00533 case UsravOp: 00534 // variable argument in an atomic operation sequence 00535 CPPAD_ASSERT_UNKNOWN( user_state == user_arg ); 00536 CPPAD_ASSERT_UNKNOWN( 0 < user_j && user_j <= user_n ); 00537 CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 ); 00538 CPPAD_ASSERT_UNKNOWN( arg[0] <= i_var ); 00539 CPPAD_ASSERT_UNKNOWN( 0 < arg[0] ); 00540 --user_j; 00541 // It might be faster if we add set union to var_sparsity 00542 // where one of the sets is not in var_sparsity. 00543 set_itr = user_r[user_j].begin(); 00544 set_end = user_r[user_j].end(); 00545 while( set_itr != set_end ) 00546 var_sparsity.add_element(arg[0], *set_itr++); 00547 if( user_j == 0 ) 00548 user_state = user_start; 00549 break; 00550 00551 case UsrrpOp: 00552 // parameter result in an atomic operation sequence 00553 CPPAD_ASSERT_UNKNOWN( user_state == user_ret ); 00554 CPPAD_ASSERT_UNKNOWN( 0 < user_i && user_i <= user_m ); 00555 CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 ); 00556 CPPAD_ASSERT_UNKNOWN( arg[0] < num_par ); 00557 --user_i; 00558 user_s[user_i].clear(); 00559 if( user_i == 0 ) 00560 { // call users function for this operation 00561 user_atomic<Base>::rev_jac_sparse(user_index, user_id, 00562 user_n, user_m, user_q, user_r, user_s 00563 ); 00564 user_state = user_arg; 00565 } 00566 break; 00567 00568 case UsrrvOp: 00569 // variable result in an atomic operation sequence 00570 CPPAD_ASSERT_UNKNOWN( user_state == user_ret ); 00571 CPPAD_ASSERT_UNKNOWN( 0 < user_i && user_i <= user_m ); 00572 --user_i; 00573 user_s[user_i].clear(); 00574 var_sparsity.begin(i_var); 00575 i = var_sparsity.next_element(); 00576 while( i < user_q ) 00577 { user_s[user_i].insert(i); 00578 i = var_sparsity.next_element(); 00579 } 00580 if( user_i == 0 ) 00581 { // call users function for this operation 00582 user_atomic<Base>::rev_jac_sparse(user_index, user_id, 00583 user_n, user_m, user_q, user_r, user_s 00584 ); 00585 user_state = user_arg; 00586 } 00587 break; 00588 // ------------------------------------------------- 00589 00590 default: 00591 CPPAD_ASSERT_UNKNOWN(0); 00592 } 00593 # if CPPAD_REV_JAC_SWEEP_TRACE 00594 for(j = 0; j < limit; j++) 00595 z_value[j] = false; 00596 var_sparsity.begin(i_var); 00597 j = var_sparsity.next_element(); 00598 while( j < limit ) 00599 { z_value[j] = true; 00600 j = var_sparsity.next_element(); 00601 } 00602 printOp( 00603 std::cout, 00604 play, 00605 i_var, 00606 op, 00607 arg, 00608 0, 00609 (CppAD::vector<bool> *) CPPAD_NULL, 00610 1, 00611 &z_value 00612 ); 00613 # endif 00614 } 00615 // values corresponding to BeginOp 00616 CPPAD_ASSERT_UNKNOWN( i_op == 0 ); 00617 CPPAD_ASSERT_UNKNOWN( i_var == 0 ); 00618 00619 if( vecad_ind != CPPAD_NULL ) 00620 CPPAD_TRACK_DEL_VEC( vecad_ind); 00621 00622 return; 00623 } 00624 00625 # undef CPPAD_REV_JAC_SWEEP_TRACE 00626 00627 CPPAD_END_NAMESPACE 00628 # endif