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