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