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