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