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