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