CppAD: A C++ Algorithmic Differentiation Package 20110419
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_REV_SPARSE_JAC_INCLUDED 00003 # define CPPAD_REV_SPARSE_JAC_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 /* 00017 $begin RevSparseJac$$ 00018 $spell 00019 std 00020 VecAD 00021 var 00022 Jacobian 00023 Jac 00024 const 00025 Bool 00026 Dep 00027 proportional 00028 $$ 00029 00030 $section Jacobian Sparsity Pattern: Reverse Mode$$ 00031 00032 $index RevSparseJac$$ 00033 $index reverse, sparse Jacobian$$ 00034 $index sparse, reverse Jacobian$$ 00035 $index pattern, reverse Jacobian$$ 00036 00037 $head Syntax$$ 00038 $icode%r% = %F%.RevSparseJac(%p%, %s%)%$$ 00039 00040 $head Purpose$$ 00041 We use $latex F : \R^n \rightarrow R^m$$ to denote the 00042 $xref/glossary/AD Function/AD function/$$ corresponding to $icode f$$. 00043 For a fixed $latex p \times m$$ matrix $latex S$$, 00044 the Jacobian of $latex S * F( x )$$ 00045 with respect to $latex x$$ is 00046 $latex \[ 00047 J(x) = S * F^{(1)} ( x ) 00048 \] $$ 00049 Given a 00050 $xref/glossary/Sparsity Pattern/sparsity pattern/$$ 00051 for $latex S$$, 00052 $code RevSparseJac$$ returns a sparsity pattern for the $latex J(x)$$. 00053 00054 $head f$$ 00055 The object $icode f$$ has prototype 00056 $codei% 00057 ADFun<%Base%> %f% 00058 %$$ 00059 00060 $head x$$ 00061 the sparsity pattern is valid for all values of the independent 00062 variables in $latex x \in \R^n$$ 00063 (even if it has $cref/CondExp/$$ or $cref/VecAD/$$ operations). 00064 00065 $head p$$ 00066 The argument $icode p$$ has prototype 00067 $codei% 00068 size_t %p% 00069 %$$ 00070 It specifies the number of rows in 00071 $latex S \in \R^{p \times m}$$ and the 00072 Jacobian $latex J(x) \in \R^{p \times n}$$. 00073 00074 $head s$$ 00075 The argument $icode s$$ has prototype 00076 $codei% 00077 const %VectorSet%& %s% 00078 %$$ 00079 (see $xref/RevSparseJac/VectorSet/VectorSet/$$ below). 00080 If it has elements of type $code bool$$, 00081 its size is $latex p * m$$. 00082 If it has elements of type $code std::set<size_t>$$, 00083 its size is $icode p$$ and all its set elements are between 00084 zero and $latex m - 1$$. 00085 It specifies a 00086 $xref/glossary/Sparsity Pattern/sparsity pattern/$$ 00087 for the matrix $icode S$$. 00088 00089 $head r$$ 00090 The return value $icode r$$ has prototype 00091 $codei% 00092 %VectorSet% %r% 00093 %$$ 00094 (see $xref/RevSparseJac/VectorSet/VectorSet/$$ below). 00095 If it has elements of type $code bool$$, 00096 its size is $latex p * n$$. 00097 If it has elements of type $code std::set<size_t>$$, 00098 its size is $icode p$$. 00099 It specifies a 00100 $xref/glossary/Sparsity Pattern/sparsity pattern/$$ 00101 for the matrix $latex J(x)$$. 00102 00103 $head VectorSet$$ 00104 The type $icode VectorSet$$ must be a $xref/SimpleVector/$$ class with 00105 $xref/SimpleVector/Elements of Specified Type/elements of type/$$ 00106 $code bool$$ or $code std::set<size_t>$$; 00107 see $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for a discussion 00108 of the difference. 00109 00110 $head Entire Sparsity Pattern$$ 00111 Suppose that $latex p = m$$ and 00112 $latex S$$ is the $latex m \times m$$ identity matrix. 00113 In this case, 00114 the corresponding value for $icode r$$ is a 00115 sparsity pattern for the Jacobian $latex J(x) = F^{(1)} ( x )$$. 00116 00117 $head Example$$ 00118 $children% 00119 example/rev_sparse_jac.cpp 00120 %$$ 00121 The file 00122 $xref/RevSparseJac.cpp/$$ 00123 contains an example and test of this operation. 00124 It returns true if it succeeds and false otherwise. 00125 00126 $end 00127 ----------------------------------------------------------------------------- 00128 */ 00129 00130 CPPAD_BEGIN_NAMESPACE 00131 /*! 00132 \file rev_sparse_jac.hpp 00133 Reverse mode Jacobian sparsity patterns. 00134 */ 00135 00136 /*! 00137 Calculate Jacobian sparsity patterns using reverse mode. 00138 00139 The C++ source code corresponding to this operation is 00140 \verbatim 00141 s = f.RevSparseJac(q, r) 00142 \endverbatim 00143 00144 \tparam Base 00145 is the base type for this recording. 00146 00147 \tparam VectorSet 00148 is a simple vector class with elements of type \c bool. 00149 00150 \param p 00151 is the number of rows in the matrix \f$ S \f$. 00152 00153 \param s 00154 is a sparsity pattern for the matrix \f$ S \f$. 00155 00156 \param r 00157 the input value of \a r must be a vector with size \c p*n 00158 where \c n is the number of independent variables 00159 corresponding to the operation sequence stored in \a play. 00160 The input value of the components of \c r does not matter. 00161 On output, \a r is the sparsity pattern for the matrix 00162 \f[ 00163 J(x) = S * F^{(1)} (x) 00164 \f] 00165 where \f$ F \f$ is the function corresponding to the operation sequence 00166 and \a x is any argument value. 00167 00168 \param total_num_var 00169 is the total number of variable in this recording. 00170 00171 \param dep_taddr 00172 maps dependendent variable index 00173 to the corresponding variable in the tape. 00174 00175 \param ind_taddr 00176 maps independent variable index 00177 to the corresponding variable in the tape. 00178 00179 \param play 00180 is the recording that defines the function we are computing the sparsity 00181 pattern for. 00182 */ 00183 00184 template <class Base, class VectorSet> 00185 void RevSparseJacBool( 00186 size_t p , 00187 const VectorSet& s , 00188 VectorSet& r , 00189 size_t total_num_var , 00190 CppAD::vector<size_t>& dep_taddr , 00191 CppAD::vector<size_t>& ind_taddr , 00192 CppAD::player<Base>& play ) 00193 { 00194 // temporary indices 00195 size_t i, j; 00196 00197 // check VectorSet is Simple Vector class with bool elements 00198 CheckSimpleVector<bool, VectorSet>(); 00199 00200 // range and domain dimensions for F 00201 size_t m = dep_taddr.size(); 00202 size_t n = ind_taddr.size(); 00203 00204 CPPAD_ASSERT_KNOWN( 00205 p > 0, 00206 "RevSparseJac: p (first argument) is not greater than zero" 00207 ); 00208 00209 CPPAD_ASSERT_KNOWN( 00210 s.size() == p * m, 00211 "RevSparseJac: s (second argument) length is not equal to\n" 00212 "p (first argument) times range dimension for ADFun object." 00213 ); 00214 00215 // vector of sets that will hold the results 00216 sparse_pack var_sparsity; 00217 var_sparsity.resize(total_num_var, p); 00218 00219 // The sparsity pattern corresponding to the dependent variables 00220 for(i = 0; i < m; i++) 00221 { CPPAD_ASSERT_UNKNOWN( dep_taddr[i] < total_num_var ); 00222 00223 for(j = 0; j < p; j++) if( s[ i * m + j ] ) 00224 var_sparsity.add_element( dep_taddr[i], j ); 00225 } 00226 00227 // evaluate the sparsity patterns 00228 RevJacSweep( 00229 n, 00230 total_num_var, 00231 &play, 00232 var_sparsity 00233 ); 00234 00235 // return values corresponding to dependent variables 00236 CPPAD_ASSERT_UNKNOWN( r.size() == p * n ); 00237 for(j = 0; j < n; j++) 00238 { CPPAD_ASSERT_UNKNOWN( ind_taddr[j] == (j+1) ); 00239 00240 // ind_taddr[j] is operator taddr for j-th independent variable 00241 CPPAD_ASSERT_UNKNOWN( play.GetOp( ind_taddr[j] ) == InvOp ); 00242 00243 // extract the result from var_sparsity 00244 for(i = 0; i < p; i++) 00245 r[ i * n + j ] = false; 00246 CPPAD_ASSERT_UNKNOWN( var_sparsity.end() == p ); 00247 var_sparsity.begin(j+1); 00248 i = var_sparsity.next_element(); 00249 while( i < p ) 00250 { r[ i * n + j ] = true; 00251 i = var_sparsity.next_element(); 00252 } 00253 } 00254 } 00255 00256 /*! 00257 \file rev_sparse_jac.hpp 00258 Reverse mode Jacobian sparsity patterns. 00259 */ 00260 00261 /*! 00262 Calculate Jacobian sparsity patterns using reverse mode. 00263 00264 The C++ source code corresponding to this operation is 00265 \verbatim 00266 s = f.RevSparseJac(q, r) 00267 \endverbatim 00268 00269 \tparam Base 00270 is the base type for this recording. 00271 00272 \tparam VectorSet 00273 is a simple vector class with elements of type \c std::set<size_t>. 00274 00275 \param p 00276 is the number of rows in the matrix \f$ S \f$. 00277 00278 \param s 00279 is a sparsity pattern for the matrix \f$ S \f$. 00280 00281 \param r 00282 the input value of \a r must be a vector with size \a p. 00283 On input, each element of \a r mus be an empty set. 00284 On output, \a r is the sparsity pattern for the matrix 00285 \f[ 00286 J(x) = S * F^{(1)} (x) 00287 \f] 00288 where \f$ F \f$ is the function corresponding to the operation sequence 00289 and \a x is any argument value. 00290 00291 \param total_num_var 00292 is the total number of variable in this recording. 00293 00294 \param dep_taddr 00295 maps dependendent variable index 00296 to the corresponding variable in the tape. 00297 00298 \param ind_taddr 00299 maps independent variable index 00300 to the corresponding variable in the tape. 00301 00302 \param play 00303 is the recording that defines the function we are computing the sparsity 00304 pattern for. 00305 */ 00306 00307 template <class Base, class VectorSet> 00308 void RevSparseJacSet( 00309 size_t p , 00310 const VectorSet& s , 00311 VectorSet& r , 00312 size_t total_num_var , 00313 CppAD::vector<size_t>& dep_taddr , 00314 CppAD::vector<size_t>& ind_taddr , 00315 CppAD::player<Base>& play ) 00316 { 00317 // temporary indices 00318 size_t i, j; 00319 std::set<size_t>::const_iterator itr; 00320 00321 // check VectorSet is Simple Vector class with sets for elements 00322 static std::set<size_t> two, three; 00323 if( two.empty() ) 00324 { two.insert(2); 00325 three.insert(3); 00326 } 00327 CPPAD_ASSERT_UNKNOWN( two.size() == 1 ); 00328 CPPAD_ASSERT_UNKNOWN( three.size() == 1 ); 00329 CheckSimpleVector<std::set<size_t>, VectorSet>(two, three); 00330 00331 // range and domain dimensions for F 00332 size_t m = dep_taddr.size(); 00333 size_t n = ind_taddr.size(); 00334 00335 CPPAD_ASSERT_KNOWN( 00336 p > 0, 00337 "RevSparseJac: p (first argument) is not greater than zero" 00338 ); 00339 00340 CPPAD_ASSERT_KNOWN( 00341 s.size() == p, 00342 "RevSparseJac: s (second argument) length is not equal to " 00343 "p (first argument)." 00344 ); 00345 00346 // vector of sets that will hold the results 00347 sparse_set var_sparsity; 00348 var_sparsity.resize(total_num_var, p); 00349 00350 // The sparsity pattern corresponding to the dependent variables 00351 for(i = 0; i < p; i++) 00352 { itr = s[i].begin(); 00353 while(itr != s[i].end()) 00354 { j = *itr++; 00355 CPPAD_ASSERT_KNOWN( 00356 j < m, 00357 "RevSparseJac: an element of the set s[i] " 00358 "has value greater than or equal m." 00359 ); 00360 CPPAD_ASSERT_UNKNOWN( dep_taddr[j] < total_num_var ); 00361 var_sparsity.add_element( dep_taddr[j], i ); 00362 } 00363 } 00364 00365 // evaluate the sparsity patterns 00366 RevJacSweep( 00367 n, 00368 total_num_var, 00369 &play, 00370 var_sparsity 00371 ); 00372 00373 // return values corresponding to dependent variables 00374 CPPAD_ASSERT_UNKNOWN( r.size() == p ); 00375 for(j = 0; j < n; j++) 00376 { CPPAD_ASSERT_UNKNOWN( ind_taddr[j] == (j+1) ); 00377 00378 // ind_taddr[j] is operator taddr for j-th independent variable 00379 CPPAD_ASSERT_UNKNOWN( play.GetOp( ind_taddr[j] ) == InvOp ); 00380 00381 // extract result from rev_hes_sparsity 00382 // and add corresponding elements to sets in r 00383 CPPAD_ASSERT_UNKNOWN( var_sparsity.end() == p ); 00384 var_sparsity.begin(j+1); 00385 i = var_sparsity.next_element(); 00386 while( i < p ) 00387 { r[i].insert(j); 00388 i = var_sparsity.next_element(); 00389 } 00390 } 00391 } 00392 00393 /*! 00394 User API for Jacobian sparsity patterns using reverse mode. 00395 00396 The C++ source code corresponding to this operation is 00397 \verbatim 00398 s = f.RevSparseJac(q, r) 00399 \endverbatim 00400 00401 \tparam Base 00402 is the base type for this recording. 00403 00404 \tparam VectorSet 00405 is a simple vector with elements of type \c bool. 00406 or \c std::set<size_t>. 00407 00408 \param p 00409 is the number of rows in the matrix \f$ S \f$. 00410 00411 \param s 00412 is a sparsity pattern for the matrix \f$ S \f$. 00413 00414 \return 00415 If \c VectorSet::value_type is \c bool, 00416 the return value \c r is a vector with size \c p*n 00417 where \c n is the number of independent variables 00418 corresponding to the operation sequence stored in \c f. 00419 If \c VectorSet::value_type is \c std::set<size_t>, 00420 the return value \c r is a vector of sets with size \c p 00421 and with all its elements between zero and \c n - 1. 00422 The value of \a r is the sparsity pattern for the matrix 00423 \f[ 00424 J(x) = S * F^{(1)} (x) 00425 \f] 00426 where \f$ F \f$ is the function corresponding to the operation sequence 00427 and \a x is any argument value. 00428 */ 00429 template <class Base> 00430 template <class VectorSet> 00431 VectorSet ADFun<Base>::RevSparseJac( 00432 size_t p , 00433 const VectorSet& s ) 00434 { 00435 VectorSet r; 00436 typedef typename VectorSet::value_type Set_type; 00437 00438 RevSparseJacCase( 00439 Set_type() , 00440 p , 00441 s , 00442 r 00443 ); 00444 return r; 00445 } 00446 00447 /*! 00448 Private helper function for \c RevSparseJac(p, s). 00449 00450 All of the description in the public member function \c RevSparseJac(p, s) 00451 applies. 00452 00453 \param set_type 00454 is a \c bool value. 00455 This arugment is used to dispatch to the proper source code 00456 depending on the value of \c VectorSet::value_type. 00457 00458 \param p 00459 See \c RevSparseJac(p, s). 00460 00461 \param s 00462 See \c RevSparseJac(p, s). 00463 00464 \param r 00465 is the return value for the corresponding call to RevSparseJac(p, s); 00466 */ 00467 00468 template <class Base> 00469 template <class VectorSet> 00470 void ADFun<Base>::RevSparseJacCase( 00471 bool set_type , 00472 size_t p , 00473 const VectorSet& s , 00474 VectorSet& r ) 00475 { size_t n = Domain(); 00476 00477 // dimension of the result vector 00478 r.resize( p * n ); 00479 00480 // store results in r 00481 RevSparseJacBool( 00482 p , 00483 s , 00484 r , 00485 total_num_var_ , 00486 dep_taddr_ , 00487 ind_taddr_ , 00488 play_ 00489 ); 00490 } 00491 00492 /*! 00493 Private helper function for \c RevSparseJac(p, s). 00494 00495 All of the description in the public member function \c RevSparseJac(p, s) 00496 applies. 00497 00498 \param set_type 00499 is a \c std::set<size_t> object. 00500 This arugment is used to dispatch to the proper source code 00501 depending on the value of \c VectorSet::value_type. 00502 00503 \param p 00504 See \c RevSparseJac(p, s). 00505 00506 \param s 00507 See \c RevSparseJac(p, s). 00508 00509 \param r 00510 is the return value for the corresponding call to RevSparseJac(p, s); 00511 */ 00512 00513 template <class Base> 00514 template <class VectorSet> 00515 void ADFun<Base>::RevSparseJacCase( 00516 const std::set<size_t>& set_type , 00517 size_t p , 00518 const VectorSet& s , 00519 VectorSet& r ) 00520 { // dimension of the result vector 00521 r.resize( p ); 00522 00523 // store results in r 00524 RevSparseJacSet( 00525 p , 00526 s , 00527 r , 00528 total_num_var_ , 00529 dep_taddr_ , 00530 ind_taddr_ , 00531 play_ 00532 ); 00533 } 00534 00535 CPPAD_END_NAMESPACE 00536 # endif