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