CppAD: A C++ Algorithmic Differentiation Package 20110419
|
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-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 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 $xref/glossary/AD Function/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 $xref/glossary/Sparsity Pattern/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 $xref/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 $xref/ForSparseJac/VectorSet/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 $xref/ForSparseJac/VectorSet/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 $xref/glossary/Sparsity Pattern/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 $xref/ForSparseJac/VectorSet/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 $xref/glossary/Sparsity Pattern/sparsity pattern/$$ 00147 for the matrix $latex S(x)$$. 00148 00149 $head VectorSet$$ 00150 The type $icode VectorSet$$ must be a $xref/SimpleVector/$$ class with 00151 $xref/SimpleVector/Elements of Specified Type/elements of 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 $xref/ForSparseJac.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 CPPAD_BEGIN_NAMESPACE 00177 /*! 00178 \file for_sparse_jac.hpp 00179 Forward mode Jacobian sparsity patterns. 00180 */ 00181 00182 /*! 00183 Calculate Jacobian sparsity patterns using forward mode. 00184 00185 The C++ source code corresponding to this operation is 00186 \verbatim 00187 s = f.ForSparseJac(q, r) 00188 \endverbatim 00189 00190 \tparam Base 00191 is the base type for this recording. 00192 00193 \tparam VectorSet 00194 is a simple vector class with elements of type \c bool. 00195 00196 \tparam Sparsity 00197 is either \c sparse_pack or \c sparse_set. 00198 00199 \param q 00200 is the number of columns in the matrix \f$ R \f$. 00201 00202 \param r 00203 is a sparsity pattern for the matrix \f$ R \f$. 00204 00205 \param s 00206 The input value of \a s must be a vector with size \c m*q 00207 where \c m is the number of dependent variables 00208 corresponding to the operation sequence stored in \a play. 00209 The input value of the components of \c s does not matter. 00210 On output, \a s is the sparsity pattern for the matrix 00211 \f[ 00212 S(x) = F^{(1)} (x) * R 00213 \f] 00214 where \f$ F \f$ is the function corresponding to the operation sequence 00215 and \a x is any argument value. 00216 00217 \param total_num_var 00218 is the total number of variable in this recording. 00219 00220 \param dep_taddr 00221 maps dependendent variable index 00222 to the corresponding variable in the tape. 00223 00224 \param ind_taddr 00225 maps independent variable index 00226 to the corresponding variable in the tape. 00227 00228 \param play 00229 is the recording that defines the function we are computing the sparsity 00230 pattern for. 00231 00232 \param for_jac_sparsity 00233 the input value of \a for_jac_sparsity does not matter. 00234 On output, \a for_jac_sparsity.n_set() == \a total_num_var 00235 and \a for_jac_sparsity.end() == \a q. 00236 It contains the forward sparsity pattern for all of the variables on the 00237 tape (given the sparsity pattern for the independent variables is \f$ R \f$). 00238 */ 00239 template <class Base, class VectorSet, class Sparsity> 00240 void ForSparseJacBool( 00241 size_t q , 00242 const VectorSet& r , 00243 VectorSet& s , 00244 size_t total_num_var , 00245 CppAD::vector<size_t>& dep_taddr , 00246 CppAD::vector<size_t>& ind_taddr , 00247 CppAD::player<Base>& play , 00248 Sparsity& for_jac_sparsity ) 00249 { 00250 // temporary indices 00251 size_t i, j; 00252 00253 // check VectorSet is Simple Vector class with bool elements 00254 CheckSimpleVector<bool, VectorSet>(); 00255 00256 // range and domain dimensions for F 00257 size_t m = dep_taddr.size(); 00258 size_t n = ind_taddr.size(); 00259 00260 CPPAD_ASSERT_KNOWN( 00261 q > 0, 00262 "ForSparseJac: q (first arugment) is not greater than zero" 00263 ); 00264 00265 CPPAD_ASSERT_KNOWN( 00266 r.size() == n * q, 00267 "ForSparseJac: r (second argument) length is not equal to\n" 00268 "q (first argument) times domain dimension for ADFun object." 00269 ); 00270 00271 // allocate memory for the requested sparsity calculation 00272 for_jac_sparsity.resize(total_num_var, q); 00273 00274 // set values corresponding to independent variables 00275 for(i = 0; i < n; i++) 00276 { CPPAD_ASSERT_UNKNOWN( ind_taddr[i] < total_num_var ); 00277 // ind_taddr[i] is operator taddr for i-th independent variable 00278 CPPAD_ASSERT_UNKNOWN( play.GetOp( ind_taddr[i] ) == InvOp ); 00279 00280 // set bits that are true 00281 for(j = 0; j < q; j++) if( r[ i * q + j ] ) 00282 for_jac_sparsity.add_element( ind_taddr[i], j); 00283 } 00284 00285 // evaluate the sparsity patterns 00286 ForJacSweep( 00287 n, 00288 total_num_var, 00289 &play, 00290 for_jac_sparsity 00291 ); 00292 00293 // return values corresponding to dependent variables 00294 CPPAD_ASSERT_UNKNOWN( s.size() == m * q ); 00295 for(i = 0; i < m; i++) 00296 { CPPAD_ASSERT_UNKNOWN( dep_taddr[i] < total_num_var ); 00297 00298 // extract the result from for_jac_sparsity 00299 for(j = 0; j < q; j++) 00300 s[ i * q + j ] = false; 00301 CPPAD_ASSERT_UNKNOWN( for_jac_sparsity.end() == q ); 00302 for_jac_sparsity.begin( dep_taddr[i] ); 00303 j = for_jac_sparsity.next_element(); 00304 while( j < q ) 00305 { s[i * q + j ] = true; 00306 j = for_jac_sparsity.next_element(); 00307 } 00308 } 00309 } 00310 00311 /*! 00312 Calculate Jacobian sparsity patterns using forward mode. 00313 00314 The C++ source code corresponding to this operation is 00315 \verbatim 00316 s = f.ForSparseJac(q, r) 00317 \endverbatim 00318 00319 \tparam Base 00320 is the base type for this recording. 00321 00322 \tparam VectorSet 00323 is a simple vector class with elements of the type \c std::set<size_t>. 00324 00325 \tparam Sparsity 00326 is either \c sparse_pack or \c sparse_set. 00327 00328 \param q 00329 is the number of columns in the matrix \f$ R \f$. 00330 00331 \param r 00332 is a sparsity pattern for the matrix \f$ R \f$. 00333 00334 \param s 00335 The input value of \a s must be a vector with size \c m 00336 where \c m is the number of dependent variables 00337 corresponding to the operation sequence stored in \a play. 00338 On input, each element of \a s must be an empty set. 00339 On output, \a s is the sparsity pattern for the matrix 00340 \f[ 00341 S(x) = F^{(1)} (x) * R 00342 \f] 00343 where \f$ F \f$ is the function corresponding to the operation sequence 00344 and \a x is any argument value. 00345 00346 \param total_num_var 00347 is the total number of variable in this recording. 00348 00349 \param dep_taddr 00350 maps dependendent variable index 00351 to the corresponding variable in the tape. 00352 00353 \param ind_taddr 00354 maps independent variable index 00355 to the corresponding variable in the tape. 00356 00357 \param play 00358 is the recording that defines the function we are computing the sparsity 00359 pattern for. 00360 00361 \param for_jac_sparsity 00362 the input value of \a for_jac_sparsity does not matter. 00363 On output, \a for_jac_sparsity.n_set() == \a total_num_var 00364 and \a for_jac_sparsity.end() == \a q. 00365 It contains the forward sparsity pattern for all of the variables on the 00366 tape (given the sparsity pattern for the independent variables is \f$ R \f$). 00367 */ 00368 00369 template <class Base, class VectorSet, class Sparsity> 00370 void ForSparseJacSet( 00371 size_t q , 00372 const VectorSet& r , 00373 VectorSet& s , 00374 size_t total_num_var , 00375 CppAD::vector<size_t>& dep_taddr , 00376 CppAD::vector<size_t>& ind_taddr , 00377 CppAD::player<Base>& play , 00378 Sparsity& for_jac_sparsity ) 00379 { 00380 // temporary indices 00381 size_t i, j; 00382 std::set<size_t>::const_iterator itr; 00383 00384 // check VectorSet is Simple Vector class with sets for elements 00385 static std::set<size_t> two, three; 00386 if( two.empty() ) 00387 { two.insert(2); 00388 three.insert(3); 00389 } 00390 CPPAD_ASSERT_UNKNOWN( two.size() == 1 ); 00391 CPPAD_ASSERT_UNKNOWN( three.size() == 1 ); 00392 CheckSimpleVector<std::set<size_t>, VectorSet>(two, three); 00393 00394 // range and domain dimensions for F 00395 size_t m = dep_taddr.size(); 00396 size_t n = ind_taddr.size(); 00397 00398 CPPAD_ASSERT_KNOWN( 00399 q > 0, 00400 "ForSparseJac: q (first arugment) is not greater than zero" 00401 ); 00402 00403 CPPAD_ASSERT_KNOWN( 00404 r.size() == n, 00405 "ForSparseJac: r (second argument) length is not equal to\n" 00406 "the domain dimension for ADFun object." 00407 ); 00408 00409 // allocate memory for the requested sparsity calculation 00410 for_jac_sparsity.resize(total_num_var, q); 00411 00412 // set values corresponding to independent variables 00413 for(i = 0; i < n; i++) 00414 { CPPAD_ASSERT_UNKNOWN( ind_taddr[i] < total_num_var ); 00415 // ind_taddr[i] is operator taddr for i-th independent variable 00416 CPPAD_ASSERT_UNKNOWN( play.GetOp( ind_taddr[i] ) == InvOp ); 00417 00418 // add the elements that are present 00419 itr = r[i].begin(); 00420 while( itr != r[i].end() ) 00421 { j = *itr++; 00422 CPPAD_ASSERT_KNOWN( 00423 j < q, 00424 "ForSparseJac: an element of the set r[i] " 00425 "has value greater than or equal q." 00426 ); 00427 for_jac_sparsity.add_element( ind_taddr[i], j); 00428 } 00429 } 00430 00431 // evaluate the sparsity patterns 00432 ForJacSweep( 00433 n, 00434 total_num_var, 00435 &play, 00436 for_jac_sparsity 00437 ); 00438 00439 // return values corresponding to dependent variables 00440 CPPAD_ASSERT_UNKNOWN( s.size() == m ); 00441 for(i = 0; i < m; i++) 00442 { CPPAD_ASSERT_UNKNOWN( dep_taddr[i] < total_num_var ); 00443 CPPAD_ASSERT_UNKNOWN( s[i].empty() ); 00444 00445 // extract results from for_jac_sparsity 00446 // and add corresponding elements to sets in s 00447 CPPAD_ASSERT_UNKNOWN( for_jac_sparsity.end() == q ); 00448 for_jac_sparsity.begin( dep_taddr[i] ); 00449 j = for_jac_sparsity.next_element(); 00450 while( j < q ) 00451 { s[i].insert(j); 00452 j = for_jac_sparsity.next_element(); 00453 } 00454 } 00455 } 00456 00457 00458 /*! 00459 User API for Jacobian sparsity patterns using forward mode. 00460 00461 The C++ source code corresponding to this operation is 00462 \verbatim 00463 s = f.ForSparseJac(q, r) 00464 \endverbatim 00465 00466 \tparam Base 00467 is the base type for this recording. 00468 00469 \tparam VectorSet 00470 is a simple vector with elements of type \c bool 00471 or \c std::set<size_t>. 00472 00473 \param q 00474 is the number of columns in the matrix \f$ R \f$. 00475 00476 \param r 00477 is a sparsity pattern for the matrix \f$ R \f$. 00478 00479 \return 00480 If \c VectorSet::value_type is \c bool, 00481 the return value \c s is a vector with size \c m*q 00482 where \c m is the number of dependent variables 00483 corresponding to the operation sequence stored in \c f. 00484 If \c VectorSet::value_type is \c std::set<size_t>, 00485 the return value \c s is a vector of sets with size \c m 00486 and with all its elements between zero and \a q - 1. 00487 The value of \a s is the sparsity pattern for the matrix 00488 \f[ 00489 S(x) = F^{(1)} (x) * R 00490 \f] 00491 where \f$ F \f$ is the function corresponding to the operation sequence 00492 and \a x is any argument value. 00493 00494 \par Side Effects 00495 If \c VectorSet::value_type is \c bool, 00496 the forward sparsity pattern for all of the variables on the 00497 tape is stored in \c for_jac_sparse_pack__. 00498 In this case 00499 \verbatim 00500 for_jac_sparse_pack_.n_set() == total_num_var_ 00501 for_jac_sparse_pack_.end() == q 00502 for_jac_sparse_set_.n_set() == 0 00503 for_jac_sparse_set_.end() == 0 00504 \endverbatim 00505 \n 00506 \n 00507 If \c VectorSet::value_type is \c std::set<size_t>, 00508 the forward sparsity pattern for all of the variables on the 00509 tape is stored in \c for_jac_sparse_set__. 00510 In this case 00511 \verbatim 00512 for_jac_sparse_set_.n_set() == total_num_var_ 00513 for_jac_sparse_set_.end() == q 00514 for_jac_sparse_pack_.n_set() == 0 00515 for_jac_sparse_pack_.end() == 0 00516 \endverbatim 00517 */ 00518 template <class Base> 00519 template <class VectorSet> 00520 VectorSet ADFun<Base>::ForSparseJac( 00521 size_t q , 00522 const VectorSet& r ) 00523 { VectorSet s; 00524 typedef typename VectorSet::value_type Set_type; 00525 00526 ForSparseJacCase( 00527 Set_type() , 00528 q , 00529 r , 00530 s 00531 ); 00532 00533 return s; 00534 } 00535 00536 /*! 00537 Private helper function for ForSparseJac(q, r). 00538 00539 All of the description in the public member function ForSparseJac(q, r) 00540 applies. 00541 00542 \param set_type 00543 is a \c bool value. This argument is used to dispatch to the proper source 00544 code depending on the value of \c VectorSet::value_type. 00545 00546 \param q 00547 See \c ForSparseJac(q, r). 00548 00549 \param r 00550 See \c ForSparseJac(q, r). 00551 00552 \param s 00553 is the return value for the corresponding call to \c ForSparseJac(q, r). 00554 */ 00555 00556 template <class Base> 00557 template <class VectorSet> 00558 void ADFun<Base>::ForSparseJacCase( 00559 bool set_type , 00560 size_t q , 00561 const VectorSet& r , 00562 VectorSet& s ) 00563 { size_t m = Range(); 00564 00565 // dimension size of result vector 00566 s.resize( m * q ); 00567 00568 // free any results stored in for_jac_sparse_set_ 00569 for_jac_sparse_set_.resize(0, 0); 00570 00571 // store results in r and for_jac_sparse_pack_ 00572 ForSparseJacBool( 00573 q , 00574 r , 00575 s , 00576 total_num_var_ , 00577 dep_taddr_ , 00578 ind_taddr_ , 00579 play_ , 00580 for_jac_sparse_pack_ 00581 ); 00582 } 00583 00584 00585 /*! 00586 Private helper function for \c ForSparseJac(q, r). 00587 00588 All of the description in the public member function \c ForSparseJac(q, r) 00589 applies. 00590 00591 \param set_type 00592 is a \c std::set<size_t> object. 00593 This argument is used to dispatch to the proper source 00594 code depending on the value of \c VectorSet::value_type. 00595 00596 \param q 00597 See \c ForSparseJac(q, r). 00598 00599 \param r 00600 See \c ForSparseJac(q, r). 00601 00602 \param s 00603 is the return value for the corresponding call to \c ForSparseJac(q, r). 00604 */ 00605 template <class Base> 00606 template <class VectorSet> 00607 void ADFun<Base>::ForSparseJacCase( 00608 const std::set<size_t>& set_type , 00609 size_t q , 00610 const VectorSet& r , 00611 VectorSet& s ) 00612 { size_t m = Range(); 00613 00614 // dimension size of result vector 00615 s.resize( m ); 00616 00617 // free any results stored in for_jac_sparse_set_ 00618 for_jac_sparse_pack_.resize(0, 0); 00619 00620 // store results in r and for_jac_sparse_pack_ 00621 CppAD::ForSparseJacSet( 00622 q , 00623 r , 00624 s , 00625 total_num_var_ , 00626 dep_taddr_ , 00627 ind_taddr_ , 00628 play_ , 00629 for_jac_sparse_set_ 00630 ); 00631 } 00632 00633 00634 CPPAD_END_NAMESPACE 00635 # endif