CppAD: A C++ Algorithmic Differentiation Package 20110419
csum_op.hpp
Go to the documentation of this file.
00001 // $Id$
00002 # ifndef CPPAD_CSUM_OP_INCLUDED
00003 # define CPPAD_CSUM_OP_INCLUDED
00004 
00005 /* --------------------------------------------------------------------------
00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-11 Bradley M. Bell
00007 
00008 CppAD is distributed under multiple licenses. This distribution is under
00009 the terms of the 
00010                     Common Public License Version 1.0.
00011 
00012 A copy of this license is included in the COPYING file of this distribution.
00013 Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
00014 -------------------------------------------------------------------------- */
00015 
00016 CPPAD_BEGIN_NAMESPACE
00017 /*!
00018 \file csum_op.hpp
00019 Forward, reverse and sparsity calculations for cummulative summation.
00020 */
00021 
00022 /*!
00023 Compute forward mode Taylor coefficients for result of op = CsumOp.
00024 
00025 This operation is 
00026 \verbatim
00027         z = p + x(1) + ... + x(m) - y(1) - ... - y(n).
00028 \endverbatim
00029 
00030 \tparam Base
00031 base type for the operator; i.e., this operation was recorded
00032 using AD< \a Base > and computations by this routine are done using type
00033 \a Base.
00034 
00035 \param d
00036 order of the Taylor coefficient that we are computing.
00037 
00038 \param i_z
00039 variable index corresponding to the result for this operation;
00040 i.e. the row index in \a taylor corresponding to z.
00041 
00042 \param arg
00043 \a arg[0] 
00044 is the number of addition variables in this cummulative summation; i.e.,
00045 <tt>m</tt>.
00046 \n
00047 \a arg[1] 
00048 is the number of subtraction variables in this cummulative summation; i.e.,
00049 \c m.
00050 \n
00051 <tt>parameter[ arg[2] ]</tt>
00052 is the parameter value \c p in this cummunative summation.
00053 \n
00054 <tt>arg[2+i]</tt>
00055 for <tt>i = 1 , ... , m</tt> is the value <tt>x(i)</tt>. 
00056 \n
00057 <tt>arg[2+arg[0]+i]</tt>
00058 for <tt>i = 1 , ... , n</tt> is the value <tt>y(i)</tt>. 
00059 
00060 \param num_par
00061 is the number of parameters in \a parameter.
00062 
00063 \param parameter
00064 is the parameter vector for this operation sequence.
00065 
00066 \param nc_taylor
00067 number of colums in the matrix containing all the Taylor coefficients.
00068 
00069 \param taylor
00070 \b Input: <tt>taylor [ arg[2+i] * nc_taylor + k ]</tt>
00071 for <tt>i = 1 , ... , m</tt> 
00072 and <tt>k = 0 , ... , d</tt>
00073 is the k-th order Taylor coefficient corresponding to <tt>x(i)</tt>
00074 \n
00075 \b Input: <tt>taylor [ arg[2+m+i] * nc_taylor + k ]</tt>
00076 for <tt>i = 1 , ... , n</tt> 
00077 and <tt>k = 0 , ... , d</tt>
00078 is the k-th order Taylor coefficient corresponding to <tt>y(i)</tt>
00079 \n
00080 \b Input: <tt>taylor [ i_z * nc_taylor + k ]</tt>
00081 for k = 0 , ... , \a d - 1
00082 is the k-th order Taylor coefficient corresponding to z.
00083 \n
00084 \b Output: <tt>taylor [ i_z * nc_taylor + d ]</tt>
00085 is the \a d-th order Taylor coefficient corresponding to z.
00086 */
00087 template <class Base>
00088 inline void forward_csum_op(
00089         size_t        d           , 
00090         size_t        i_z         ,
00091         const size_t* arg         ,
00092         size_t        num_par     ,
00093         const Base*   parameter   ,
00094         size_t        nc_taylor   ,
00095         Base*         taylor      )
00096 {       static Base zero(0);
00097 
00098         // check assumptions
00099         CPPAD_ASSERT_UNKNOWN( NumRes(CSumOp) == 1 );
00100         CPPAD_ASSERT_UNKNOWN( d < nc_taylor );
00101         CPPAD_ASSERT_UNKNOWN( arg[2] < num_par );
00102         CPPAD_ASSERT_UNKNOWN( 
00103                 arg[0] + arg[1] == arg[ arg[0] + arg[1] + 3 ]
00104         ); 
00105 
00106         // Taylor coefficients corresponding to result
00107         Base* z = taylor + i_z    * nc_taylor;
00108         if( d == 0 )
00109                 z[d] = parameter[ arg[2] ];
00110         else    z[d] = zero;
00111         Base* x;
00112         size_t i, j;
00113         i = arg[0];
00114         j = 2;
00115         while(i--)
00116         {       CPPAD_ASSERT_UNKNOWN( arg[j+1] < i_z );
00117                 x     = taylor + arg[++j] * nc_taylor;
00118                 z[d] += x[d];
00119         }       
00120         i = arg[1];
00121         while(i--)
00122         {       CPPAD_ASSERT_UNKNOWN( arg[j+1] < i_z );
00123                 x     = taylor + arg[++j] * nc_taylor;
00124                 z[d] -= x[d];
00125         }       
00126 }
00127 
00128 /*!
00129 Compute reverse mode Taylor coefficients for result of op = CsumOp.
00130 
00131 This operation is 
00132 \verbatim
00133         z = p + x(1) + ... + x(m) - y(1) - ... - y(n).
00134         H(y, x, w, ...) = G[ z(x, y), y, x, w, ... ] 
00135 \endverbatim
00136 
00137 \tparam Base
00138 base type for the operator; i.e., this operation was recorded
00139 using AD< \a Base > and computations by this routine are done using type
00140 \a Base.
00141 
00142 \param d
00143 order the highest order Taylor coefficient that we are computing
00144 the partial derivatives with respect to.
00145 
00146 \param i_z
00147 variable index corresponding to the result for this operation;
00148 i.e. the row index in \a taylor corresponding to z.
00149 
00150 \param arg
00151 \a arg[0] 
00152 is the number of addition variables in this cummulative summation; i.e.,
00153 <tt>m</tt>.
00154 \n
00155 \a arg[1] 
00156 is the number of subtraction variables in this cummulative summation; i.e.,
00157 \c m.
00158 \n
00159 <tt>parameter[ arg[2] ]</tt>
00160 is the parameter value \c p in this cummunative summation.
00161 \n
00162 <tt>arg[2+i]</tt>
00163 for <tt>i = 1 , ... , m</tt> is the value <tt>x(i)</tt>. 
00164 \n
00165 <tt>arg[2+arg[0]+i]</tt>
00166 for <tt>i = 1 , ... , n</tt> is the value <tt>y(i)</tt>. 
00167 
00168 \param nc_partial
00169 number of colums in the matrix containing all the partial derivatives.
00170 
00171 \param partial
00172 \b Input: <tt>partial [ arg[2+i] * nc_partial + k ]</tt>
00173 for <tt>i = 1 , ... , m</tt> 
00174 and <tt>k = 0 , ... , d</tt>
00175 is the partial derivative of G(z, y, x, w, ...) with respect to the
00176 k-th order Taylor coefficient corresponding to <tt>x(i)</tt>
00177 \n
00178 \b Input: <tt>partial [ arg[2+m+i] * nc_partial + k ]</tt>
00179 for <tt>i = 1 , ... , n</tt> 
00180 and <tt>k = 0 , ... , d</tt>
00181 is the partial derivative of G(z, y, x, w, ...) with respect to the
00182 k-th order Taylor coefficient corresponding to <tt>y(i)</tt>
00183 \n
00184 \b Input: <tt>partial [ i_z * nc_partial + k ]</tt>
00185 for <tt>i = 1 , ... , n</tt> 
00186 and <tt>k = 0 , ... , d</tt>
00187 is the partial derivative of G(z, y, x, w, ...) with respect to the
00188 k-th order Taylor coefficient corresponding to \c z.
00189 \n
00190 \b Output: <tt>partial [ arg[2+i] * nc_partial + k ]</tt>
00191 for <tt>i = 1 , ... , m</tt> 
00192 and <tt>k = 0 , ... , d</tt>
00193 is the partial derivative of H(y, x, w, ...) with respect to the
00194 k-th order Taylor coefficient corresponding to <tt>x(i)</tt>
00195 \n
00196 \b Output: <tt>partial [ arg[2+m+i] * nc_partial + k ]</tt>
00197 for <tt>i = 1 , ... , n</tt> 
00198 and <tt>k = 0 , ... , d</tt>
00199 is the partial derivative of H(y, x, w, ...) with respect to the
00200 k-th order Taylor coefficient corresponding to <tt>y(i)</tt>
00201 */
00202 
00203 template <class Base>
00204 inline void reverse_csum_op(
00205         size_t        d           , 
00206         size_t        i_z         ,
00207         const size_t* arg         ,
00208         size_t        nc_partial  ,
00209         Base*         partial     )
00210 {
00211         // check assumptions
00212         CPPAD_ASSERT_UNKNOWN( NumRes(CSumOp) == 1 );
00213         CPPAD_ASSERT_UNKNOWN( d < nc_partial );
00214 
00215         // Taylor coefficients and partial derivative corresponding to result
00216         Base* pz = partial + i_z * nc_partial;
00217         Base* px;
00218         size_t i, j, k;
00219         size_t d1 = d + 1;
00220         i = arg[0];
00221         j = 2;
00222         while(i--)
00223         {       CPPAD_ASSERT_UNKNOWN( arg[j+1] < i_z );
00224                 px    = partial + arg[++j] * nc_partial;
00225                 k = d1;
00226                 while(k--)
00227                         px[k] += pz[k];
00228         }       
00229         i = arg[1];
00230         while(i--)
00231         {       CPPAD_ASSERT_UNKNOWN( arg[j+1] < i_z );
00232                 px    = partial + arg[++j] * nc_partial;
00233                 k = d1;
00234                 while(k--)
00235                         px[k] -= pz[k];
00236         }       
00237 }
00238 
00239 
00240 /*!
00241 Forward mode Jacobian sparsity pattern for CSumOp operator. 
00242 
00243 This operation is 
00244 \verbatim
00245         z = p + x(1) + ... + x(m) - y(1) - ... - y(n).
00246 \endverbatim
00247 
00248 \tparam Vector_set
00249 is the type used for vectors of sets. It can be either
00250 \c sparse_pack or \c sparse_set.
00251 
00252 \param i_z
00253 variable index corresponding to the result for this operation;
00254 i.e. the index in \a sparsity corresponding to z.
00255 
00256 \param arg
00257 \a arg[0] 
00258 is the number of addition variables in this cummulative summation; i.e.,
00259 <tt>m + n</tt>.
00260 \n
00261 \a arg[1] 
00262 is the number of subtraction variables in this cummulative summation; i.e.,
00263 \c m.
00264 \n
00265 <tt>parameter[ arg[2] ]</tt>
00266 is the parameter value \c p in this cummunative summation.
00267 \n
00268 <tt>arg[2+i]</tt>
00269 for <tt>i = 1 , ... , m</tt> is the value <tt>x(i)</tt>. 
00270 \n
00271 <tt>arg[2+arg[1]+i]</tt>
00272 for <tt>i = 1 , ... , n</tt> is the value <tt>y(i)</tt>. 
00273 
00274 \param sparsity
00275 \b Input: 
00276 For <tt>i = 1 , ... , m</tt>,
00277 the set with index \a arg[2+i] in \a sparsity
00278 is the sparsity bit pattern for <tt>x(i)</tt>.
00279 This identifies which of the independent variables the variable 
00280 <tt>x(i)</tt> depends on. 
00281 \n
00282 \b Input: 
00283 For <tt>i = 1 , ... , n</tt>,
00284 the set with index \a arg[2+arg[0]+i] in \a sparsity
00285 is the sparsity bit pattern for <tt>x(i)</tt>.
00286 This identifies which of the independent variables the variable 
00287 <tt>y(i)</tt> depends on. 
00288 \n
00289 \b Output: 
00290 The set with index \a i_z in \a sparsity
00291 is the sparsity bit pattern for z.
00292 This identifies which of the independent variables the variable z
00293 depends on. 
00294 */
00295 
00296 template <class Vector_set>
00297 inline void forward_sparse_jacobian_csum_op(
00298         size_t           i_z         ,
00299         const size_t*    arg         ,
00300         Vector_set&      sparsity    )
00301 {       sparsity.clear(i_z);
00302 
00303         size_t i, j;
00304         i = arg[0] + arg[1];
00305         j = 2;
00306         while(i--)
00307         {       CPPAD_ASSERT_UNKNOWN( arg[j+1] < i_z );
00308                 sparsity.binary_union(
00309                         i_z        , // index in sparsity for result
00310                         i_z        , // index in sparsity for left operand
00311                         arg[++j]   , // index for right operand 
00312                         sparsity     // sparsity vector for right operand
00313                 );
00314         }       
00315 }
00316 
00317 /*!
00318 Reverse mode Jacobian sparsity pattern for CSumOp operator. 
00319 
00320 This operation is 
00321 \verbatim
00322         z = p + x(1) + ... + x(m) - y(1) - ... - y(n).
00323         H(y, x, w, ...) = G[ z(x, y), y, x, w, ... ] 
00324 \endverbatim
00325 
00326 \tparam Vector_set
00327 is the type used for vectors of sets. It can be either
00328 \c sparse_pack or \c sparse_set.
00329 
00330 \param i_z
00331 variable index corresponding to the result for this operation;
00332 i.e. the index in \a sparsity corresponding to z.
00333 
00334 \param arg
00335 \a arg[0] 
00336 is the number of addition variables in this cummulative summation; i.e.,
00337 <tt>m + n</tt>.
00338 \n
00339 \a arg[1] 
00340 is the number of subtraction variables in this cummulative summation; i.e.,
00341 \c m.
00342 \n
00343 <tt>parameter[ arg[2] ]</tt>
00344 is the parameter value \c p in this cummunative summation.
00345 \n
00346 <tt>arg[2+i]</tt>
00347 for <tt>i = 1 , ... , m</tt> is the value <tt>x(i)</tt>. 
00348 \n
00349 <tt>arg[2+arg[1]+i]</tt>
00350 for <tt>i = 1 , ... , n</tt> is the value <tt>y(i)</tt>. 
00351 
00352 \param sparsity
00353 For <tt>i = 1 , ... , m</tt>,
00354 the set with index \a arg[2+i] in \a sparsity
00355 is the sparsity bit pattern for <tt>x(i)</tt>.
00356 This identifies which of the dependent variables depend on <tt>x(i)</tt>. 
00357 On input, the sparsity patter corresponds to \c G,
00358 and on ouput it corresponds to \c H.
00359 \n
00360 For <tt>i = 1 , ... , m</tt>,
00361 the set with index \a arg[2+arg[0]+i] in \a sparsity
00362 is the sparsity bit pattern for <tt>y(i)</tt>.
00363 This identifies which of the dependent variables depend on <tt>y(i)</tt>. 
00364 On input, the sparsity patter corresponds to \c G,
00365 and on ouput it corresponds to \c H.
00366 \n
00367 \b Input: 
00368 The set with index \a i_z in \a sparsity
00369 is the sparsity bit pattern for z.
00370 On input it corresponds to \c G and on output it is undefined.
00371 */
00372 
00373 template <class Vector_set>
00374 inline void reverse_sparse_jacobian_csum_op(
00375         size_t           i_z         ,
00376         const size_t*    arg         ,
00377         Vector_set&      sparsity    )
00378 {
00379         size_t i, j;
00380         i = arg[0] + arg[1];
00381         j = 2;
00382         while(i--)
00383         {       ++j;
00384                 CPPAD_ASSERT_UNKNOWN( arg[j] < i_z );
00385                 sparsity.binary_union(
00386                         arg[j]     , // index in sparsity for result
00387                         arg[j]     , // index in sparsity for left operand
00388                         i_z        , // index for right operand 
00389                         sparsity     // sparsity vector for right operand
00390                 );
00391         }       
00392 }
00393 /*!
00394 Reverse mode Hessian sparsity pattern for CSumOp operator. 
00395 
00396 This operation is 
00397 \verbatim
00398         z = p + x(1) + ... + x(m) - y(1) - ... - y(n).
00399         H(y, x, w, ...) = G[ z(x, y), y, x, w, ... ] 
00400 \endverbatim
00401 
00402 \tparam Vector_set
00403 is the type used for vectors of sets. It can be either
00404 \c sparse_pack or \c sparse_set.
00405 
00406 \param i_z
00407 variable index corresponding to the result for this operation;
00408 i.e. the index in \a sparsity corresponding to z.
00409 
00410 \param arg
00411 \a arg[0] 
00412 is the number of addition variables in this cummulative summation; i.e.,
00413 <tt>m + n</tt>.
00414 \n
00415 \a arg[1] 
00416 is the number of subtraction variables in this cummulative summation; i.e.,
00417 \c m.
00418 \n
00419 <tt>parameter[ arg[2] ]</tt>
00420 is the parameter value \c p in this cummunative summation.
00421 \n
00422 <tt>arg[2+i]</tt>
00423 for <tt>i = 1 , ... , m</tt> is the value <tt>x(i)</tt>. 
00424 \n
00425 <tt>arg[2+arg[0]+i]</tt>
00426 for <tt>i = 1 , ... , n</tt> is the value <tt>y(i)</tt>. 
00427 
00428 \param rev_jacobian
00429 <tt>rev_jacobian[i_z]</tt>
00430 is all false (true) if the Jabobian of G with respect to z must be zero
00431 (may be non-zero).
00432 \n
00433 \n
00434 For <tt>i = 1 , ... , m</tt> 
00435 <tt>rev_jacobian[ arg[2+i] ]</tt>
00436 is all false (true) if the Jacobian with respect to <tt>x(i)</tt>
00437 is zero (may be non-zero).
00438 On input, it corresponds to the function G,
00439 and on output it corresponds to the function H.
00440 \n
00441 \n
00442 For <tt>i = 1 , ... , n</tt> 
00443 <tt>rev_jacobian[ arg[2+arg[0]+i] ]</tt>
00444 is all false (true) if the Jacobian with respect to <tt>y(i)</tt>
00445 is zero (may be non-zero).
00446 On input, it corresponds to the function G,
00447 and on output it corresponds to the function H.
00448 
00449 \param rev_hes_sparsity
00450 The set with index \a i_z in in \a rev_hes_sparsity
00451 is the Hessian sparsity pattern for the fucntion G
00452 where one of the partials derivative is with respect to z.
00453 \n
00454 \n
00455 For <tt>i = 1 , ... , m</tt> 
00456 The set with index <tt>arg[2+i]</tt> in \a rev_hes_sparsity
00457 is the Hessian sparsity pattern 
00458 where one of the partials derivative is with respect to <tt>x(i)</tt>.
00459 On input, it corresponds to the function G,
00460 and on output it corresponds to the function H.
00461 \n
00462 \n
00463 For <tt>i = 1 , ... , n</tt> 
00464 The set with index <tt>arg[2+arg[0]+i]</tt> in \a rev_hes_sparsity
00465 is the Hessian sparsity pattern 
00466 where one of the partials derivative is with respect to <tt>y(i)</tt>.
00467 On input, it corresponds to the function G,
00468 and on output it corresponds to the function H.
00469 */
00470 
00471 template <class Vector_set>
00472 inline void reverse_sparse_hessian_csum_op(
00473         size_t           i_z                 ,
00474         const size_t*    arg                 ,
00475         bool*            rev_jacobian        ,
00476         Vector_set&      rev_hes_sparsity    )
00477 {
00478         size_t i, j;
00479         i = arg[0] + arg[1];
00480         j = 2;
00481         while(i--)
00482         {       ++j;
00483                 CPPAD_ASSERT_UNKNOWN( arg[j] < i_z );
00484                 rev_hes_sparsity.binary_union(
00485                 arg[j]             , // index in sparsity for result
00486                 arg[j]             , // index in sparsity for left operand
00487                 i_z                , // index for right operand 
00488                 rev_hes_sparsity     // sparsity vector for right operand
00489                 );
00490                 rev_jacobian[arg[j]] |= rev_jacobian[i_z];
00491         }       
00492 }
00493 
00494 CPPAD_END_NAMESPACE
00495 # endif