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