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