CppAD: A C++ Algorithmic Differentiation Package 20110419
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_COND_OP_INCLUDED 00003 # define CPPAD_COND_OP_INCLUDED 00004 /* -------------------------------------------------------------------------- 00005 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-10 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 00015 CPPAD_BEGIN_NAMESPACE 00016 /*! 00017 \file cond_op.hpp 00018 Forward, reverse, and sparse operations for conditional expressions. 00019 */ 00020 00021 /*! 00022 Compute forward mode Taylor coefficients for op = CExpOp. 00023 00024 \copydetails conditional_exp_op 00025 00026 \param d 00027 is the order of the Taylor coefficient of z that we are computing. 00028 00029 \param taylor 00030 \b Input: 00031 For j = 0, 1, 2, 3 and k = 0 , ... , \a d, 00032 if y_j is a variable then 00033 \a taylor [ \a arg[2+j] * nc_taylor + k ] 00034 is the k-th order Taylor coefficient corresponding to y_j. 00035 \n 00036 \b Input: \a taylor [ \a i_z * \a nc_taylor + k ] 00037 for k = 0 , ... , \a d - 1 00038 is the k-th order Taylor coefficient corresponding to z. 00039 \n 00040 \b Output: \a taylor [ \a i_z * \a nc_taylor + \a d ] 00041 is the \a d-th order Taylor coefficient corresponding to z. 00042 00043 */ 00044 template <class Base> 00045 inline void forward_cond_op( 00046 size_t d , 00047 size_t i_z , 00048 const size_t* arg , 00049 size_t num_par , 00050 const Base* parameter , 00051 size_t nc_taylor , 00052 Base* taylor ) 00053 { Base y_0, y_1, y_2, y_3; 00054 Base zero(0); 00055 Base* z; 00056 00057 CPPAD_ASSERT_UNKNOWN( arg[0] < static_cast<size_t> (CompareNe) ); 00058 CPPAD_ASSERT_UNKNOWN( NumArg(CExpOp) == 6 ); 00059 CPPAD_ASSERT_UNKNOWN( NumRes(CExpOp) == 1 ); 00060 CPPAD_ASSERT_UNKNOWN( arg[1] != 0 ); 00061 00062 if( arg[1] & 1 ) 00063 { CPPAD_ASSERT_UNKNOWN( arg[2] < i_z ); 00064 y_0 = taylor[ arg[2] * nc_taylor + 0 ]; 00065 } 00066 else 00067 { CPPAD_ASSERT_UNKNOWN( arg[2] < num_par ); 00068 y_0 = parameter[ arg[2] ]; 00069 } 00070 if( arg[1] & 2 ) 00071 { CPPAD_ASSERT_UNKNOWN( arg[3] < i_z ); 00072 y_1 = taylor[ arg[3] * nc_taylor + 0 ]; 00073 } 00074 else 00075 { CPPAD_ASSERT_UNKNOWN( arg[3] < num_par ); 00076 y_1 = parameter[ arg[3] ]; 00077 } 00078 # if CPPAD_USE_FORWARD0SWEEP 00079 CPPAD_ASSERT_UNKNOWN( d > 0 ); 00080 # else 00081 if( d == 0 ) 00082 { if( arg[1] & 4 ) 00083 { CPPAD_ASSERT_UNKNOWN( arg[4] < i_z ); 00084 y_2 = taylor[ arg[4] * nc_taylor + 0 ]; 00085 } 00086 else 00087 { CPPAD_ASSERT_UNKNOWN( arg[4] < num_par ); 00088 y_2 = parameter[ arg[4] ]; 00089 } 00090 if( arg[1] & 8 ) 00091 { CPPAD_ASSERT_UNKNOWN( arg[5] < i_z ); 00092 y_3 = taylor[ arg[5] * nc_taylor + 0 ]; 00093 } 00094 else 00095 { CPPAD_ASSERT_UNKNOWN( arg[5] < num_par ); 00096 y_3 = parameter[ arg[5] ]; 00097 } 00098 } 00099 else 00100 # endif 00101 { if( arg[1] & 4 ) 00102 { CPPAD_ASSERT_UNKNOWN( arg[4] < i_z ); 00103 y_2 = taylor[ arg[4] * nc_taylor + d]; 00104 } 00105 else y_2 = zero; 00106 if( arg[1] & 8 ) 00107 { CPPAD_ASSERT_UNKNOWN( arg[5] < i_z ); 00108 y_3 = taylor[ arg[5] * nc_taylor + d]; 00109 } 00110 else y_3 = zero; 00111 } 00112 z = taylor + i_z * nc_taylor; 00113 z[d] = CondExpOp( 00114 CompareOp( arg[0] ), 00115 y_0, 00116 y_1, 00117 y_2, 00118 y_3 00119 ); 00120 return; 00121 } 00122 00123 /*! 00124 Compute zero order forward mode Taylor coefficients for op = CExpOp. 00125 00126 \copydetails conditional_exp_op 00127 00128 \param taylor 00129 \b Input: 00130 For j = 0, 1, 2, 3, 00131 if y_j is a variable then 00132 \a taylor [ \a arg[2+j] * nc_taylor + 0 ] 00133 is the zero order Taylor coefficient corresponding to y_j. 00134 \n 00135 \b Output: \a taylor [ \a i_z * \a nc_taylor + 0 ] 00136 is the zero order Taylor coefficient corresponding to z. 00137 */ 00138 template <class Base> 00139 inline void forward_cond_op_0( 00140 size_t i_z , 00141 const size_t* arg , 00142 size_t num_par , 00143 const Base* parameter , 00144 size_t nc_taylor , 00145 Base* taylor ) 00146 { Base y_0, y_1, y_2, y_3; 00147 Base* z; 00148 00149 CPPAD_ASSERT_UNKNOWN( arg[0] < static_cast<size_t> (CompareNe) ); 00150 CPPAD_ASSERT_UNKNOWN( NumArg(CExpOp) == 6 ); 00151 CPPAD_ASSERT_UNKNOWN( NumRes(CExpOp) == 1 ); 00152 CPPAD_ASSERT_UNKNOWN( arg[1] != 0 ); 00153 00154 if( arg[1] & 1 ) 00155 { CPPAD_ASSERT_UNKNOWN( arg[2] < i_z ); 00156 y_0 = taylor[ arg[2] * nc_taylor + 0 ]; 00157 } 00158 else 00159 { CPPAD_ASSERT_UNKNOWN( arg[2] < num_par ); 00160 y_0 = parameter[ arg[2] ]; 00161 } 00162 if( arg[1] & 2 ) 00163 { CPPAD_ASSERT_UNKNOWN( arg[3] < i_z ); 00164 y_1 = taylor[ arg[3] * nc_taylor + 0 ]; 00165 } 00166 else 00167 { CPPAD_ASSERT_UNKNOWN( arg[3] < num_par ); 00168 y_1 = parameter[ arg[3] ]; 00169 } 00170 if( arg[1] & 4 ) 00171 { CPPAD_ASSERT_UNKNOWN( arg[4] < i_z ); 00172 y_2 = taylor[ arg[4] * nc_taylor + 0 ]; 00173 } 00174 else 00175 { CPPAD_ASSERT_UNKNOWN( arg[4] < num_par ); 00176 y_2 = parameter[ arg[4] ]; 00177 } 00178 if( arg[1] & 8 ) 00179 { CPPAD_ASSERT_UNKNOWN( arg[5] < i_z ); 00180 y_3 = taylor[ arg[5] * nc_taylor + 0 ]; 00181 } 00182 else 00183 { CPPAD_ASSERT_UNKNOWN( arg[5] < num_par ); 00184 y_3 = parameter[ arg[5] ]; 00185 } 00186 z = taylor + i_z * nc_taylor; 00187 z[0] = CondExpOp( 00188 CompareOp( arg[0] ), 00189 y_0, 00190 y_1, 00191 y_2, 00192 y_3 00193 ); 00194 return; 00195 } 00196 00197 /*! 00198 Compute reverse mode Taylor coefficients for op = CExpOp. 00199 00200 \copydetails conditional_exp_op 00201 00202 This routine is given the partial derivatives of a function 00203 G( z , y , x , w , ... ) 00204 and it uses them to compute the partial derivatives of 00205 \verbatim 00206 H( y , x , w , u , ... ) = G[ z(y) , y , x , w , u , ... ] 00207 \endverbatim 00208 where y above represents y_0, y_1, y_2, y_3. 00209 00210 \param d 00211 is the order of the Taylor coefficient of z that we are computing. 00212 00213 \param taylor 00214 \b Input: 00215 For j = 0, 1, 2, 3 and k = 0 , ... , \a d, 00216 if y_j is a variable then 00217 \a taylor [ \a arg[2+j] * nc_taylor + k ] 00218 is the k-th order Taylor coefficient corresponding to y_j. 00219 \n 00220 \a taylor [ \a i_z * \a nc_taylor + k ] 00221 for k = 0 , ... , \a d 00222 is the k-th order Taylor coefficient corresponding to z. 00223 00224 \param nc_partial 00225 number of columns in the matrix containing the Taylor coefficients. 00226 00227 \param partial 00228 \b Input: 00229 For j = 0, 1, 2, 3 and k = 0 , ... , \a d, 00230 if y_j is a variable then 00231 \a partial [ \a arg[2+j] * nc_partial + k ] 00232 is the partial derivative of G( z , y , x , w , u , ... ) 00233 with respect to the k-th order Taylor coefficient corresponding to y_j. 00234 \n 00235 \b Input: \a partial [ \a i_z * \a nc_taylor + k ] 00236 for k = 0 , ... , \a d 00237 is the partial derivative of G( z , y , x , w , u , ... ) 00238 with respect to the k-th order Taylor coefficient corresponding to z. 00239 \n 00240 \b Output: 00241 For j = 0, 1, 2, 3 and k = 0 , ... , \a d, 00242 if y_j is a variable then 00243 \a partial [ \a arg[2+j] * nc_partial + k ] 00244 is the partial derivative of H( y , x , w , u , ... ) 00245 with respect to the k-th order Taylor coefficient corresponding to y_j. 00246 00247 */ 00248 template <class Base> 00249 inline void reverse_cond_op( 00250 size_t d , 00251 size_t i_z , 00252 const size_t* arg , 00253 size_t num_par , 00254 const Base* parameter , 00255 size_t nc_taylor , 00256 const Base* taylor , 00257 size_t nc_partial , 00258 Base* partial ) 00259 { Base y_0, y_1; 00260 Base zero(0); 00261 Base* pz; 00262 Base* py_2; 00263 Base* py_3; 00264 00265 CPPAD_ASSERT_UNKNOWN( arg[0] < static_cast<size_t> (CompareNe) ); 00266 CPPAD_ASSERT_UNKNOWN( NumArg(CExpOp) == 6 ); 00267 CPPAD_ASSERT_UNKNOWN( NumRes(CExpOp) == 1 ); 00268 CPPAD_ASSERT_UNKNOWN( arg[1] != 0 ); 00269 00270 pz = partial + i_z * nc_partial + 0; 00271 if( arg[1] & 1 ) 00272 { CPPAD_ASSERT_UNKNOWN( arg[2] < i_z ); 00273 y_0 = taylor[ arg[2] * nc_taylor + 0 ]; 00274 } 00275 else 00276 { CPPAD_ASSERT_UNKNOWN( arg[2] < num_par ); 00277 y_0 = parameter[ arg[2] ]; 00278 } 00279 if( arg[1] & 2 ) 00280 { CPPAD_ASSERT_UNKNOWN( arg[3] < i_z ); 00281 y_1 = taylor[ arg[3] * nc_taylor + 0 ]; 00282 } 00283 else 00284 { CPPAD_ASSERT_UNKNOWN( arg[3] < num_par ); 00285 y_1 = parameter[ arg[3] ]; 00286 } 00287 if( arg[1] & 4 ) 00288 { CPPAD_ASSERT_UNKNOWN( arg[4] < i_z ); 00289 py_2 = partial + arg[4] * nc_partial; 00290 size_t j = d + 1; 00291 while(j--) 00292 { py_2[j] += CondExpOp( 00293 CompareOp( arg[0] ), 00294 y_0, 00295 y_1, 00296 pz[j], 00297 zero 00298 ); 00299 } 00300 } 00301 if( arg[1] & 8 ) 00302 { CPPAD_ASSERT_UNKNOWN( arg[5] < i_z ); 00303 py_3 = partial + arg[5] * nc_partial; 00304 size_t j = d + 1; 00305 while(j--) 00306 { py_3[j] += CondExpOp( 00307 CompareOp( arg[0] ), 00308 y_0, 00309 y_1, 00310 zero, 00311 pz[j] 00312 ); 00313 } 00314 } 00315 return; 00316 } 00317 00318 /*! 00319 Compute forward Jacobian sparsity patterns for op = CExpOp. 00320 00321 \copydetails sparse_conditional_exp_op 00322 00323 \param sparsity 00324 \b Input: 00325 if y_2 is a variable, the set with index t is 00326 the sparsity pattern corresponding to y_2. 00327 This identifies which of the independent variables the variable y_2 00328 depends on. 00329 \n 00330 \b Input: 00331 if y_3 is a variable, the set with index t is 00332 the sparsity pattern corresponding to y_3. 00333 This identifies which of the independent variables the variable y_3 00334 depends on. 00335 \n 00336 \b Output: 00337 The set with index T is 00338 the sparsity pattern corresponding to z. 00339 This identifies which of the independent variables the variable z 00340 depends on. 00341 */ 00342 template <class Vector_set> 00343 inline void forward_sparse_jacobian_cond_op( 00344 size_t i_z , 00345 const size_t* arg , 00346 size_t num_par , 00347 Vector_set& sparsity ) 00348 { 00349 CPPAD_ASSERT_UNKNOWN( arg[0] < static_cast<size_t> (CompareNe) ); 00350 CPPAD_ASSERT_UNKNOWN( NumArg(CExpOp) == 6 ); 00351 CPPAD_ASSERT_UNKNOWN( NumRes(CExpOp) == 1 ); 00352 CPPAD_ASSERT_UNKNOWN( arg[1] != 0 ); 00353 00354 # ifndef NDEBUG 00355 if( arg[1] & 1 ) 00356 { CPPAD_ASSERT_UNKNOWN( arg[2] < i_z ); 00357 } 00358 else 00359 { CPPAD_ASSERT_UNKNOWN( arg[2] < num_par ); 00360 } 00361 if( arg[1] & 2 ) 00362 { CPPAD_ASSERT_UNKNOWN( arg[3] < i_z ); 00363 } 00364 else 00365 { CPPAD_ASSERT_UNKNOWN( arg[3] < num_par ); 00366 } 00367 # endif 00368 if( arg[1] & 4 ) 00369 { CPPAD_ASSERT_UNKNOWN( arg[4] < i_z ); 00370 if( arg[1] & 8 ) 00371 { CPPAD_ASSERT_UNKNOWN( arg[5] < i_z ); 00372 sparsity.binary_union(i_z, arg[4], arg[5], sparsity); 00373 } 00374 else 00375 { CPPAD_ASSERT_UNKNOWN( arg[5] < num_par ); 00376 sparsity.assignment(i_z, arg[4], sparsity); 00377 } 00378 } 00379 else 00380 { CPPAD_ASSERT_UNKNOWN( arg[4] < num_par ); 00381 if( arg[1] & 8 ) 00382 { CPPAD_ASSERT_UNKNOWN( arg[5] < i_z ); 00383 sparsity.assignment(i_z, arg[5], sparsity); 00384 } 00385 else 00386 { CPPAD_ASSERT_UNKNOWN( arg[5] < num_par ); 00387 sparsity.clear(i_z); 00388 } 00389 } 00390 return; 00391 } 00392 00393 /*! 00394 Compute reverse Jacobian sparsity patterns for op = CExpOp. 00395 00396 This routine is given the sparsity patterns 00397 for a function G(z, y, x, ... ) 00398 and it uses them to compute the sparsity patterns for 00399 \verbatim 00400 H( y, x, w , u , ... ) = G[ z(x,y) , y , x , w , u , ... ] 00401 \endverbatim 00402 where y represents the combination of y_0, y_1, y_2, and y_3. 00403 00404 \copydetails sparse_conditional_exp_op 00405 00406 00407 \param sparsity 00408 if y_2 is a variable, the set with index t is 00409 the sparsity pattern corresponding to y_2. 00410 This identifies which of the dependent variables depend on the variable y_2. 00411 On input, this pattern corresponds to the function G. 00412 On ouput, it corresponds to the function H. 00413 \n 00414 \n 00415 if y_3 is a variable, the set with index t is 00416 the sparsity pattern corresponding to y_3. 00417 This identifies which of the dependent variables depeond on the variable y_3. 00418 On input, this pattern corresponds to the function G. 00419 On ouput, it corresponds to the function H. 00420 \n 00421 \b Output: 00422 The set with index T is 00423 the sparsity pattern corresponding to z. 00424 This identifies which of the dependent variables depend on the variable z. 00425 On input and output, this pattern corresponds to the function G. 00426 */ 00427 template <class Vector_set> 00428 inline void reverse_sparse_jacobian_cond_op( 00429 size_t i_z , 00430 const size_t* arg , 00431 size_t num_par , 00432 Vector_set& sparsity ) 00433 { 00434 CPPAD_ASSERT_UNKNOWN( arg[0] < static_cast<size_t> (CompareNe) ); 00435 CPPAD_ASSERT_UNKNOWN( NumArg(CExpOp) == 6 ); 00436 CPPAD_ASSERT_UNKNOWN( NumRes(CExpOp) == 1 ); 00437 CPPAD_ASSERT_UNKNOWN( arg[1] != 0 ); 00438 00439 # ifndef NDEBUG 00440 if( arg[1] & 1 ) 00441 { CPPAD_ASSERT_UNKNOWN( arg[2] < i_z ); 00442 } 00443 else 00444 { CPPAD_ASSERT_UNKNOWN( arg[2] < num_par ); 00445 } 00446 if( arg[1] & 2 ) 00447 { CPPAD_ASSERT_UNKNOWN( arg[3] < i_z ); 00448 } 00449 else 00450 { CPPAD_ASSERT_UNKNOWN( arg[3] < num_par ); 00451 } 00452 if( ! ( arg[1] & 4 ) ) 00453 { CPPAD_ASSERT_UNKNOWN( arg[4] < num_par ); 00454 } 00455 if( ! ( arg[1] & 8 ) ) 00456 { CPPAD_ASSERT_UNKNOWN( arg[5] < num_par ); 00457 } 00458 # endif 00459 if( arg[1] & 4 ) 00460 { CPPAD_ASSERT_UNKNOWN( arg[4] < i_z ); 00461 sparsity.binary_union(arg[4], arg[4], i_z, sparsity); 00462 } 00463 if( arg[1] & 8 ) 00464 { CPPAD_ASSERT_UNKNOWN( arg[5] < i_z ); 00465 sparsity.binary_union(arg[5], arg[5], i_z, sparsity); 00466 } 00467 return; 00468 } 00469 00470 /*! 00471 Compute reverse Hessian sparsity patterns for op = CExpOp. 00472 00473 This routine is given the sparsity patterns 00474 for a function G(z, y, x, ... ) 00475 and it uses them to compute the sparsity patterns for 00476 \verbatim 00477 H( y, x, w , u , ... ) = G[ z(x,y) , y , x , w , u , ... ] 00478 \endverbatim 00479 where y represents the combination of y_0, y_1, y_2, and y_3. 00480 00481 \copydetails sparse_conditional_exp_op 00482 00483 00484 \param jac_reverse 00485 \a jac_reverse[i_z] 00486 is false (true) if the Jacobian of G with respect to z is always zero 00487 (may be non-zero). 00488 \n 00489 \n 00490 \a jac_reverse[ arg[4] ] 00491 If y_2 is a variable, 00492 \a jac_reverse[ arg[4] ] 00493 is false (true) if the Jacobian with respect to y_2 is always zero 00494 (may be non-zero). 00495 On input, it corresponds to the function G, 00496 and on output it corresponds to the function H. 00497 \n 00498 \n 00499 \a jac_reverse[ arg[5] ] 00500 If y_3 is a variable, 00501 \a jac_reverse[ arg[5] ] 00502 is false (true) if the Jacobian with respect to y_3 is always zero 00503 (may be non-zero). 00504 On input, it corresponds to the function G, 00505 and on output it corresponds to the function H. 00506 00507 \param hes_sparsity 00508 The set with index \a i_z in \a hes_sparsity 00509 is the Hessian sparsity pattern for the function G 00510 where one of the partials is with respect to z. 00511 \n 00512 \n 00513 If y_2 is a variable, 00514 the set with index \a arg[4] in \a hes_sparsity 00515 is the Hessian sparsity pattern 00516 where one of the partials is with respect to y_2. 00517 On input, this pattern corresponds to the function G. 00518 On output, this pattern corresponds to the function H. 00519 \n 00520 \n 00521 If y_3 is a variable, 00522 the set with index \a arg[5] in \a hes_sparsity 00523 is the Hessian sparsity pattern 00524 where one of the partials is with respect to y_3. 00525 On input, this pattern corresponds to the function G. 00526 On output, this pattern corresponds to the function H. 00527 */ 00528 template <class Vector_set> 00529 inline void reverse_sparse_hessian_cond_op( 00530 size_t i_z , 00531 const size_t* arg , 00532 size_t num_par , 00533 bool* jac_reverse , 00534 Vector_set& hes_sparsity ) 00535 { 00536 00537 CPPAD_ASSERT_UNKNOWN( arg[0] < static_cast<size_t> (CompareNe) ); 00538 CPPAD_ASSERT_UNKNOWN( NumArg(CExpOp) == 6 ); 00539 CPPAD_ASSERT_UNKNOWN( NumRes(CExpOp) == 1 ); 00540 CPPAD_ASSERT_UNKNOWN( arg[1] != 0 ); 00541 00542 # ifndef NDEBUG 00543 if( arg[1] & 1 ) 00544 { CPPAD_ASSERT_UNKNOWN( arg[2] < i_z ); 00545 } 00546 else 00547 { CPPAD_ASSERT_UNKNOWN( arg[2] < num_par ); 00548 } 00549 if( arg[1] & 2 ) 00550 { CPPAD_ASSERT_UNKNOWN( arg[3] < i_z ); 00551 } 00552 else 00553 { CPPAD_ASSERT_UNKNOWN( arg[3] < num_par ); 00554 } 00555 if( ! ( arg[1] & 4 ) ) 00556 { CPPAD_ASSERT_UNKNOWN( arg[4] < num_par ); 00557 } 00558 if( ! ( arg[1] & 8 ) ) 00559 { CPPAD_ASSERT_UNKNOWN( arg[5] < num_par ); 00560 } 00561 # endif 00562 if( arg[1] & 4 ) 00563 { CPPAD_ASSERT_UNKNOWN( arg[4] < i_z ); 00564 00565 hes_sparsity.binary_union(arg[4], arg[4], i_z, hes_sparsity); 00566 jac_reverse[ arg[4] ] |= jac_reverse[i_z]; 00567 } 00568 if( arg[1] & 8 ) 00569 { CPPAD_ASSERT_UNKNOWN( arg[5] < i_z ); 00570 00571 hes_sparsity.binary_union(arg[5], arg[5], i_z, hes_sparsity); 00572 jac_reverse[ arg[5] ] |= jac_reverse[i_z]; 00573 } 00574 return; 00575 } 00576 00577 CPPAD_END_NAMESPACE 00578 # endif