CppAD: A C++ Algorithmic Differentiation Package
20130102
|
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-13 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 /* 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 : \B{R}^n \rightarrow R^m$$ to denote the 00042 $cref/AD function/glossary/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 $cref/sparsity pattern/glossary/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 \B{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 \B{R}^{p \times m}$$ and the 00072 Jacobian $latex J(x) \in \B{R}^{p \times n}$$. 00073 00074 $head s$$ 00075 The argument $icode s$$ has prototype 00076 $codei% 00077 const %VectorSet%& %s% 00078 %$$ 00079 (see $cref/VectorSet/RevSparseJac/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 $cref/sparsity pattern/glossary/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 $cref/VectorSet/RevSparseJac/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 $cref/sparsity pattern/glossary/Sparsity Pattern/$$ 00101 for the matrix $latex J(x)$$. 00102 00103 $head VectorSet$$ 00104 The type $icode VectorSet$$ must be a $cref SimpleVector$$ class with 00105 $cref/elements of type/SimpleVector/Elements of Specified 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 $cref rev_sparse_jac.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 # include <cppad/local/std_set.hpp> 00131 00132 CPPAD_BEGIN_NAMESPACE 00133 /*! 00134 \defgroup rev_sparse_jac_hpp rev_sparse_jac.hpp 00135 \{ 00136 \file rev_sparse_jac.hpp 00137 Reverse mode Jacobian sparsity patterns. 00138 */ 00139 00140 // ------------------------------------------------------------------------- 00141 /*! 00142 Calculate Jacobian vector of bools sparsity patterns using reverse mode. 00143 00144 The C++ source code corresponding to this operation is 00145 \verbatim 00146 s = f.RevSparseJac(q, r) 00147 \endverbatim 00148 00149 \tparam Base 00150 is the base type for this recording. 00151 00152 \tparam VectorSet 00153 is a simple vector class with elements of type \c bool. 00154 00155 \param p 00156 is the number of rows in the matrix \f$ S \f$. 00157 00158 \param s 00159 is a sparsity pattern for the matrix \f$ S \f$. 00160 00161 \param r 00162 the input value of \a r must be a vector with size \c p*n 00163 where \c n is the number of independent variables 00164 corresponding to the operation sequence stored in \a play. 00165 The input value of the components of \c r does not matter. 00166 On output, \a r is the sparsity pattern for the matrix 00167 \f[ 00168 J(x) = S * F^{(1)} (x) 00169 \f] 00170 where \f$ F \f$ is the function corresponding to the operation sequence 00171 and \a x is any argument value. 00172 00173 \param total_num_var 00174 is the total number of variable in this recording. 00175 00176 \param dep_taddr 00177 maps dependendent variable index 00178 to the corresponding variable in the tape. 00179 00180 \param ind_taddr 00181 maps independent variable index 00182 to the corresponding variable in the tape. 00183 00184 \param play 00185 is the recording that defines the function we are computing the sparsity 00186 pattern for. 00187 */ 00188 00189 template <class Base, class VectorSet> 00190 void RevSparseJacBool( 00191 size_t p , 00192 const VectorSet& s , 00193 VectorSet& r , 00194 size_t total_num_var , 00195 CppAD::vector<size_t>& dep_taddr , 00196 CppAD::vector<size_t>& ind_taddr , 00197 CppAD::player<Base>& play ) 00198 { 00199 // temporary indices 00200 size_t i, j; 00201 00202 // check VectorSet is Simple Vector class with bool elements 00203 CheckSimpleVector<bool, VectorSet>(); 00204 00205 // range and domain dimensions for F 00206 size_t m = dep_taddr.size(); 00207 size_t n = ind_taddr.size(); 00208 00209 CPPAD_ASSERT_KNOWN( 00210 p > 0, 00211 "RevSparseJac: p (first argument) is not greater than zero" 00212 ); 00213 00214 CPPAD_ASSERT_KNOWN( 00215 size_t(s.size()) == p * m, 00216 "RevSparseJac: s (second argument) length is not equal to\n" 00217 "p (first argument) times range dimension for ADFun object." 00218 ); 00219 00220 // vector of sets that will hold the results 00221 sparse_pack var_sparsity; 00222 var_sparsity.resize(total_num_var, p); 00223 00224 // The sparsity pattern corresponding to the dependent variables 00225 for(i = 0; i < m; i++) 00226 { CPPAD_ASSERT_UNKNOWN( dep_taddr[i] < total_num_var ); 00227 00228 for(j = 0; j < p; j++) if( s[ i * m + j ] ) 00229 var_sparsity.add_element( dep_taddr[i], j ); 00230 } 00231 00232 // evaluate the sparsity patterns 00233 RevJacSweep( 00234 n, 00235 total_num_var, 00236 &play, 00237 var_sparsity 00238 ); 00239 00240 // return values corresponding to dependent variables 00241 CPPAD_ASSERT_UNKNOWN( size_t(r.size()) == p * n ); 00242 for(j = 0; j < n; j++) 00243 { CPPAD_ASSERT_UNKNOWN( ind_taddr[j] == (j+1) ); 00244 00245 // ind_taddr[j] is operator taddr for j-th independent variable 00246 CPPAD_ASSERT_UNKNOWN( play.GetOp( ind_taddr[j] ) == InvOp ); 00247 00248 // extract the result from var_sparsity 00249 for(i = 0; i < p; i++) 00250 r[ i * n + j ] = false; 00251 CPPAD_ASSERT_UNKNOWN( var_sparsity.end() == p ); 00252 var_sparsity.begin(j+1); 00253 i = var_sparsity.next_element(); 00254 while( i < p ) 00255 { r[ i * n + j ] = true; 00256 i = var_sparsity.next_element(); 00257 } 00258 } 00259 } 00260 /*! 00261 Calculate Jacobian vector of sets sparsity patterns using reverse mode. 00262 00263 The C++ source code corresponding to this operation is 00264 \verbatim 00265 s = f.RevSparseJac(q, r) 00266 \endverbatim 00267 00268 \tparam Base 00269 see \c RevSparseJacBool. 00270 00271 \tparam VectorSet 00272 is a simple vector class with elements of type \c std::set<size_t>. 00273 00274 \param p 00275 see \c RevSparseJacBool. 00276 00277 \param s 00278 see \c RevSparseJacBool. 00279 00280 \param r 00281 see \c RevSparseJacBool. 00282 00283 \param total_num_var 00284 see \c RevSparseJacBool. 00285 00286 \param dep_taddr 00287 see \c RevSparseJacBool. 00288 00289 \param ind_taddr 00290 see \c RevSparseJacBool. 00291 00292 \param play 00293 see \c RevSparseJacBool. 00294 */ 00295 template <class Base, class VectorSet> 00296 void RevSparseJacSet( 00297 size_t p , 00298 const VectorSet& s , 00299 VectorSet& r , 00300 size_t total_num_var , 00301 CppAD::vector<size_t>& dep_taddr , 00302 CppAD::vector<size_t>& ind_taddr , 00303 CppAD::player<Base>& play ) 00304 { 00305 // temporary indices 00306 size_t i, j; 00307 std::set<size_t>::const_iterator itr; 00308 00309 // check VectorSet is Simple Vector class with sets for elements 00310 CheckSimpleVector<std::set<size_t>, VectorSet>( 00311 one_element_std_set<size_t>(), two_element_std_set<size_t>() 00312 ); 00313 00314 // domain dimensions for F 00315 size_t n = ind_taddr.size(); 00316 00317 CPPAD_ASSERT_KNOWN( 00318 p > 0, 00319 "RevSparseJac: p (first argument) is not greater than zero" 00320 ); 00321 00322 CPPAD_ASSERT_KNOWN( 00323 size_t(s.size()) == p, 00324 "RevSparseJac: s (second argument) length is not equal to " 00325 "p (first argument)." 00326 ); 00327 00328 // vector of lists that will hold the results 00329 CPPAD_INTERNAL_SPARSE_SET var_sparsity; 00330 var_sparsity.resize(total_num_var, p); 00331 00332 // The sparsity pattern corresponding to the dependent variables 00333 for(i = 0; i < p; i++) 00334 { itr = s[i].begin(); 00335 while(itr != s[i].end()) 00336 { j = *itr++; 00337 CPPAD_ASSERT_KNOWN( 00338 j < dep_taddr.size(), 00339 "RevSparseJac: an element of the set s[i] " 00340 "has value greater than or equal Range dimension." 00341 ); 00342 CPPAD_ASSERT_UNKNOWN( dep_taddr[j] < total_num_var ); 00343 var_sparsity.add_element( dep_taddr[j], i ); 00344 } 00345 } 00346 00347 // evaluate the sparsity patterns 00348 RevJacSweep( 00349 n, 00350 total_num_var, 00351 &play, 00352 var_sparsity 00353 ); 00354 00355 // return values corresponding to dependent variables 00356 CPPAD_ASSERT_UNKNOWN( size_t(r.size()) == p ); 00357 for(j = 0; j < n; j++) 00358 { CPPAD_ASSERT_UNKNOWN( ind_taddr[j] == (j+1) ); 00359 00360 // ind_taddr[j] is operator taddr for j-th independent variable 00361 CPPAD_ASSERT_UNKNOWN( play.GetOp( ind_taddr[j] ) == InvOp ); 00362 00363 // extract result from rev_hes_sparsity 00364 // and add corresponding elements to sets in r 00365 CPPAD_ASSERT_UNKNOWN( var_sparsity.end() == p ); 00366 var_sparsity.begin(j+1); 00367 i = var_sparsity.next_element(); 00368 while( i < p ) 00369 { r[i].insert(j); 00370 i = var_sparsity.next_element(); 00371 } 00372 } 00373 } 00374 // -------------------------------------------------------------------------- 00375 00376 /*! 00377 Private helper function for \c RevSparseJac(p, s). 00378 00379 All of the description in the public member function \c RevSparseJac(p, s) 00380 applies. 00381 00382 \param set_type 00383 is a \c bool value. 00384 This arugment is used to dispatch to the proper source code 00385 depending on the value of \c VectorSet::value_type. 00386 00387 \param p 00388 See \c RevSparseJac(p, s). 00389 00390 \param s 00391 See \c RevSparseJac(p, s). 00392 00393 \param r 00394 is the return value for the corresponding call to RevSparseJac(p, s); 00395 */ 00396 00397 template <class Base> 00398 template <class VectorSet> 00399 void ADFun<Base>::RevSparseJacCase( 00400 bool set_type , 00401 size_t p , 00402 const VectorSet& s , 00403 VectorSet& r ) 00404 { size_t n = Domain(); 00405 00406 // dimension of the result vector 00407 r.resize( p * n ); 00408 00409 // store results in r 00410 RevSparseJacBool( 00411 p , 00412 s , 00413 r , 00414 total_num_var_ , 00415 dep_taddr_ , 00416 ind_taddr_ , 00417 play_ 00418 ); 00419 } 00420 00421 /*! 00422 Private helper function for \c RevSparseJac(p, s). 00423 00424 All of the description in the public member function \c RevSparseJac(p, s) 00425 applies. 00426 00427 \param set_type 00428 is a \c std::set<size_t> object. 00429 This arugment is used to dispatch to the proper source code 00430 depending on the value of \c VectorSet::value_type. 00431 00432 \param p 00433 See \c RevSparseJac(p, s). 00434 00435 \param s 00436 See \c RevSparseJac(p, s). 00437 00438 \param r 00439 is the return value for the corresponding call to RevSparseJac(p, s); 00440 */ 00441 00442 template <class Base> 00443 template <class VectorSet> 00444 void ADFun<Base>::RevSparseJacCase( 00445 const std::set<size_t>& set_type , 00446 size_t p , 00447 const VectorSet& s , 00448 VectorSet& r ) 00449 { // dimension of the result vector 00450 r.resize( p ); 00451 00452 // store results in r 00453 RevSparseJacSet( 00454 p , 00455 s , 00456 r , 00457 total_num_var_ , 00458 dep_taddr_ , 00459 ind_taddr_ , 00460 play_ 00461 ); 00462 } 00463 // -------------------------------------------------------------------------- 00464 /*! 00465 User API for Jacobian sparsity patterns using reverse mode. 00466 00467 The C++ source code corresponding to this operation is 00468 \verbatim 00469 s = f.RevSparseJac(q, r) 00470 \endverbatim 00471 00472 \tparam Base 00473 is the base type for this recording. 00474 00475 \tparam VectorSet 00476 is a simple vector with elements of type \c bool. 00477 or \c std::set<size_t>. 00478 00479 \param p 00480 is the number of rows in the matrix \f$ S \f$. 00481 00482 \param s 00483 is a sparsity pattern for the matrix \f$ S \f$. 00484 00485 \return 00486 If \c VectorSet::value_type is \c bool, 00487 the return value \c r is a vector with size \c p*n 00488 where \c n is the number of independent variables 00489 corresponding to the operation sequence stored in \c f. 00490 If \c VectorSet::value_type is \c std::set<size_t>, 00491 the return value \c r is a vector of sets with size \c p 00492 and with all its elements between zero and \c n - 1. 00493 The value of \a r is the sparsity pattern for the matrix 00494 \f[ 00495 J(x) = S * F^{(1)} (x) 00496 \f] 00497 where \f$ F \f$ is the function corresponding to the operation sequence 00498 and \a x is any argument value. 00499 */ 00500 template <class Base> 00501 template <class VectorSet> 00502 VectorSet ADFun<Base>::RevSparseJac( 00503 size_t p , 00504 const VectorSet& s ) 00505 { 00506 VectorSet r; 00507 typedef typename VectorSet::value_type Set_type; 00508 00509 RevSparseJacCase( 00510 Set_type() , 00511 p , 00512 s , 00513 r 00514 ); 00515 return r; 00516 } 00517 00518 /*! \} */ 00519 CPPAD_END_NAMESPACE 00520 # endif