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