CppAD: A C++ Algorithmic Differentiation Package 20110419
cond_op.hpp
Go to the documentation of this file.
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