CppAD: A C++ Algorithmic Differentiation Package
20130102
|
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