CppAD: A C++ Algorithmic Differentiation Package 20110419
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_SPARSE_JACOBIAN_INCLUDED 00003 # define CPPAD_SPARSE_JACOBIAN_INCLUDED 00004 00005 /* -------------------------------------------------------------------------- 00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-10 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 sparse_jacobian$$ 00018 $spell 00019 valarray 00020 std 00021 CppAD 00022 Bool 00023 jac 00024 Jacobian 00025 const 00026 Taylor 00027 $$ 00028 00029 $section Sparse Jacobian: Easy Driver$$ 00030 $index SparseJacobian$$ 00031 $index jacobian, sparse$$ 00032 00033 $head Syntax$$ 00034 $codei%%jac% = %f%.SparseJacobian(%x%) 00035 %$$ 00036 $codei%%jac% = %f%.SparseJacobian(%x%, %p%)%$$ 00037 00038 $head Purpose$$ 00039 We use $latex F : \R^n \rightarrow \R^m$$ do denote the 00040 $cref/AD function/glossary/AD Function/$$ 00041 corresponding to $icode f$$. 00042 The syntax above sets $icode jac$$ to the Jacobian 00043 $latex \[ 00044 jac = F^{(1)} (x) 00045 \] $$ 00046 This routine assumes 00047 that the matrix $latex F^{(1)} (x) \in \R^{m \times n}$$ is sparse 00048 and uses this assumption to reduce the amount of computation necessary. 00049 One should use speed tests (e.g. $cref/speed_test/$$) 00050 to verify that results are computed faster 00051 than when using the routine $cref/Jacobian/$$. 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 (see $cref/Uses Forward/sparse_jacobian/Uses Forward/$$ below). 00060 00061 $head x$$ 00062 The argument $icode x$$ has prototype 00063 $codei% 00064 const %VectorBase%& %x% 00065 %$$ 00066 (see $cref/VectorBase/sparse_jacobian/VectorBase/$$ below) 00067 and its size 00068 must be equal to $icode n$$, the dimension of the 00069 $cref/domain/seq_property/Domain/$$ space for $icode f$$. 00070 It specifies 00071 that point at which to evaluate the Jacobian. 00072 00073 $head p$$ 00074 The argument $icode p$$ is optional and has prototype 00075 $syntax% 00076 const %VectorSet%& %p% 00077 %$$ 00078 (see $cref/VectorSet/sparse_jacobian/VectorSet/$$ below). 00079 If it has elements of type $code bool$$, 00080 its size is $latex m * n$$. 00081 If it has elements of type $code std::set<size_t>$$, 00082 its size is $latex m$$ and all its set elements are between 00083 zero and $latex n - 1$$. 00084 It specifies a 00085 $cref/sparsity pattern/glossary/Sparsity Pattern/$$ 00086 for the Jacobian $latex F^{(1)} (x)$$. 00087 $pre 00088 00089 $$ 00090 If this sparsity pattern does not change between calls to 00091 $codei SparseJacobian$$, it should be faster to calculate $icode p$$ once 00092 (using $cref/ForSparseJac/$$ or $cref/RevSparseJac/$$) 00093 and then pass $icode p$$ to $codei SparseJacobian$$. 00094 In addition, 00095 if you specify $icode p$$, CppAD will use the same 00096 type of sparsity representation 00097 (vectors of $code bool$$ or vectors of $code std::set<size_t>$$) 00098 for its internal calculations. 00099 Otherwise, the representation 00100 for the internal calculations is unspecified. 00101 00102 $head jac$$ 00103 The result $icode jac$$ has prototype 00104 $codei% 00105 %VectorBase% %jac% 00106 %$$ 00107 and its size is $latex m * n$$. 00108 For $latex i = 0 , \ldots , m - 1$$, 00109 and $latex j = 0 , \ldots , n - 1 $$ 00110 $latex \[ 00111 jac [ i * n + j ] = \D{ F_i }{ x_j } (x) 00112 \] $$ 00113 00114 $head VectorBase$$ 00115 The type $icode VectorBase$$ must be a $cref/SimpleVector/$$ class with 00116 $cref/elements of type/SimpleVector/Elements of Specified Type/$$ 00117 $icode Base$$. 00118 The routine $cref/CheckSimpleVector/$$ will generate an error message 00119 if this is not the case. 00120 00121 $head VectorSet$$ 00122 The type $icode VectorSet$$ must be a $xref/SimpleVector/$$ class with 00123 $xref/SimpleVector/Elements of Specified Type/elements of type/$$ 00124 $code bool$$ or $code std::set<size_t>$$; 00125 see $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for a discussion 00126 of the difference. 00127 The routine $cref/CheckSimpleVector/$$ will generate an error message 00128 if this is not the case. 00129 00130 $subhead Restrictions$$ 00131 If $icode VectorSet$$ has elements of $code std::set<size_t>$$, 00132 then $icode%p%[%i%]%$$ must return a reference (not a copy) to the 00133 corresponding set. 00134 According to section 26.3.2.3 of the 1998 C++ standard, 00135 $code std::valarray< std::set<size_t> >$$ does not satisfy 00136 this condition. 00137 00138 $head Uses Forward$$ 00139 After each call to $cref/Forward/$$, 00140 the object $icode f$$ contains the corresponding 00141 $cref/Taylor coefficients/glossary/Taylor Coefficient/$$. 00142 After $code SparseJacobian$$, 00143 the previous calls to $xref/Forward/$$ are undefined. 00144 00145 $head Example$$ 00146 $children% 00147 example/sparse_jacobian.cpp 00148 %$$ 00149 The routine 00150 $cref/sparse_jacobian.cpp/$$ 00151 is examples and tests of $code sparse_jacobian$$. 00152 It return $code true$$, if it succeeds and $code false$$ otherwise. 00153 00154 $end 00155 ----------------------------------------------------------------------------- 00156 */ 00157 00158 CPPAD_BEGIN_NAMESPACE 00159 /*! 00160 \file sparse_jacobian.hpp 00161 Sparse Jacobian driver routine and helper functions. 00162 */ 00163 00164 /*! 00165 Private helper function for SparseJacobian(x, p). 00166 00167 All of the description in the public member function SparseJacobian(x, p) 00168 applies. 00169 00170 \param set_type 00171 is a \c bool value. This argument is used to dispatch to the proper souce 00172 code depending on the value of \c VectorSet::value_type. 00173 00174 \param x 00175 See \c SparseJacobian(x, p). 00176 00177 \param p 00178 See \c SparseJacobian(x, p). 00179 00180 \param jac 00181 is the return value for the corresponding call to \c SparseJacobian(x, p). 00182 On input, it must have size equal to the domain times range dimension 00183 for this ADFun<Base> object. 00184 On return, it will contain the Jacobian. 00185 */ 00186 template <class Base> 00187 template <class VectorBase, class VectorSet> 00188 void ADFun<Base>::SparseJacobianCase( 00189 bool set_type , 00190 const VectorBase& x , 00191 const VectorSet& p , 00192 VectorBase& jac ) 00193 { 00194 typedef CppAD::vector<size_t> SizeVector; 00195 typedef CppAD::vectorBool VectorBool; 00196 size_t i, j, k; 00197 00198 size_t m = Range(); 00199 size_t n = Domain(); 00200 00201 // some values 00202 const Base zero(0); 00203 const Base one(1); 00204 00205 // check VectorSet is Simple Vector class with bool elements 00206 CheckSimpleVector<bool, VectorSet>(); 00207 00208 // check VectorBase is Simple Vector class with Base type elements 00209 CheckSimpleVector<Base, VectorBase>(); 00210 00211 CPPAD_ASSERT_KNOWN( 00212 x.size() == n, 00213 "SparseJacobian: size of x not equal domain dimension for f" 00214 ); 00215 CPPAD_ASSERT_KNOWN( 00216 p.size() == m * n, 00217 "SparseJacobian: using bool values and size of p " 00218 " not equal range dimension times domain dimension for f" 00219 ); 00220 CPPAD_ASSERT_UNKNOWN(jac.size() == m * n); 00221 00222 // point at which we are evaluating the Jacobian 00223 Forward(0, x); 00224 00225 // initialize the return value 00226 for(i = 0; i < m; i++) 00227 for(j = 0; j < n; j++) 00228 jac[i * n + j] = zero; 00229 00230 if( n <= m ) 00231 { // use forward mode ---------------------------------------- 00232 00233 // initial coloring 00234 SizeVector color(n); 00235 for(j = 0; j < n; j++) 00236 color[j] = j; 00237 00238 // See GreedyPartialD2Coloring Algorithm Section 3.6.2 of 00239 // Graph Coloring in Optimization Revisited by 00240 // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen 00241 VectorBool forbidden(n); 00242 for(j = 0; j < n; j++) 00243 { // initial all colors as ok for this column 00244 for(k = 0; k < n; k++) 00245 forbidden[k] = false; 00246 // for each row that is connected to column j 00247 for(i = 0; i < m; i++) if( p[i * n + j] ) 00248 { // for each column that is connected to row i 00249 for(k = 0; k < n; k++) 00250 if( p[i * n + k] & (j != k) ) 00251 forbidden[ color[k] ] = true; 00252 } 00253 k = 0; 00254 while( forbidden[k] && k < n ) 00255 { k++; 00256 CPPAD_ASSERT_UNKNOWN( k < n ); 00257 } 00258 color[j] = k; 00259 } 00260 size_t n_color = 1; 00261 for(k = 0; k < n; k++) 00262 n_color = std::max(n_color, color[k] + 1); 00263 00264 // direction vector for calls to forward 00265 VectorBase dx(n); 00266 00267 // location for return values from Reverse 00268 VectorBase dy(m); 00269 00270 // loop over colors 00271 size_t c; 00272 for(c = 0; c < n_color; c++) 00273 { // determine all the colums with this color 00274 for(j = 0; j < n; j++) 00275 { if( color[j] == c ) 00276 dx[j] = one; 00277 else dx[j] = zero; 00278 } 00279 // call forward mode for all these columns at once 00280 dy = Forward(1, dx); 00281 00282 // set the corresponding components of the result 00283 for(j = 0; j < n; j++) if( color[j] == c ) 00284 { for(i = 0; i < m; i++) 00285 if( p[ i * n + j ] ) 00286 jac[i * n + j] = dy[i]; 00287 } 00288 } 00289 } 00290 else 00291 { // use reverse mode ---------------------------------------- 00292 00293 // initial coloring 00294 SizeVector color(m); 00295 for(i = 0; i < m; i++) 00296 color[i] = i; 00297 00298 // See GreedyPartialD2Coloring Algorithm Section 3.6.2 of 00299 // Graph Coloring in Optimization Revisited by 00300 // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen 00301 VectorBool forbidden(m); 00302 for(i = 0; i < m; i++) 00303 { // initial all colors as ok for this row 00304 for(k = 0; k < m; k++) 00305 forbidden[k] = false; 00306 // for each column that is connected to row i 00307 for(j = 0; j < n; j++) if( p[i * n + j] ) 00308 { // for each row that is connected to column j 00309 for(k = 0; k < m; k++) 00310 if( p[k * n + j] & (i != k) ) 00311 forbidden[ color[k] ] = true; 00312 } 00313 k = 0; 00314 while( forbidden[k] && k < m ) 00315 { k++; 00316 CPPAD_ASSERT_UNKNOWN( k < n ); 00317 } 00318 color[i] = k; 00319 } 00320 size_t n_color = 1; 00321 for(k = 0; k < m; k++) 00322 n_color = std::max(n_color, color[k] + 1); 00323 00324 // weight vector for calls to reverse 00325 VectorBase w(m); 00326 00327 // location for return values from Reverse 00328 VectorBase dw(n); 00329 00330 // loop over colors 00331 size_t c; 00332 for(c = 0; c < n_color; c++) 00333 { // determine all the rows with this color 00334 for(i = 0; i < m; i++) 00335 { if( color[i] == c ) 00336 w[i] = one; 00337 else w[i] = zero; 00338 } 00339 // call reverse mode for all these rows at once 00340 dw = Reverse(1, w); 00341 00342 // set the corresponding components of the result 00343 for(i = 0; i < m; i++) if( color[i] == c ) 00344 { for(j = 0; j < n; j++) 00345 if( p[ i * n + j ] ) 00346 jac[i * n + j] = dw[j]; 00347 } 00348 } 00349 } 00350 } 00351 00352 /*! 00353 Private helper function for SparseJacobian(x, p). 00354 00355 All of the description in the public member function SparseJacobian(x, p) 00356 applies. 00357 00358 \param set_type 00359 is a \c std::set<size_t> value. 00360 This argument is used to dispatch to the proper souce 00361 code depending on the vlaue of \c VectorSet::value_type. 00362 00363 \param x 00364 See \c SparseJacobian(x, p). 00365 00366 \param p 00367 See \c SparseJacobian(x, p). 00368 00369 \param jac 00370 is the return value for the corresponding call to \c SparseJacobian(x, p). 00371 On input it must have size equalt to the domain times range dimension 00372 for this ADFun<Base> object. 00373 On return, it will contain the Jacobian. 00374 */ 00375 template <class Base> 00376 template <class VectorBase, class VectorSet> 00377 void ADFun<Base>::SparseJacobianCase( 00378 const std::set<size_t>& set_type , 00379 const VectorBase& x , 00380 const VectorSet& p , 00381 VectorBase& jac ) 00382 { 00383 typedef CppAD::vector<size_t> SizeVector; 00384 typedef CppAD::vectorBool VectorBool; 00385 size_t i, j, k; 00386 00387 size_t m = Range(); 00388 size_t n = Domain(); 00389 00390 // some values 00391 const Base zero(0); 00392 const Base one(1); 00393 00394 // check VectorSet is Simple Vector class with sets for elements 00395 static std::set<size_t> two, three; 00396 if( two.empty() ) 00397 { two.insert(2); 00398 three.insert(3); 00399 } 00400 CPPAD_ASSERT_UNKNOWN( two.size() == 1 ); 00401 CPPAD_ASSERT_UNKNOWN( three.size() == 1 ); 00402 CheckSimpleVector<std::set<size_t>, VectorSet>(two, three); 00403 00404 // check VectorBase is Simple Vector class with Base type elements 00405 CheckSimpleVector<Base, VectorBase>(); 00406 00407 CPPAD_ASSERT_KNOWN( 00408 x.size() == n, 00409 "SparseJacobian: size of x not equal domain dimension for f" 00410 ); 00411 CPPAD_ASSERT_KNOWN( 00412 p.size() == m, 00413 "SparseJacobian: using sets and size of p " 00414 "not equal range dimension for f" 00415 ); 00416 CPPAD_ASSERT_UNKNOWN(jac.size() == m * n); 00417 00418 // point at which we are evaluating the Jacobian 00419 Forward(0, x); 00420 00421 // initialize the return value 00422 for(i = 0; i < m; i++) 00423 for(j = 0; j < n; j++) 00424 jac[i * n + j] = zero; 00425 00426 // create a copy of the transpose sparsity pattern 00427 VectorSet q(n); 00428 std::set<size_t>::const_iterator itr_i, itr_j; 00429 for(i = 0; i < m; i++) 00430 { itr_j = p[i].begin(); 00431 while( itr_j != p[i].end() ) 00432 { j = *itr_j++; 00433 q[j].insert(i); 00434 } 00435 } 00436 00437 if( n <= m ) 00438 { // use forward mode ---------------------------------------- 00439 00440 // initial coloring 00441 SizeVector color(n); 00442 for(j = 0; j < n; j++) 00443 color[j] = j; 00444 00445 // See GreedyPartialD2Coloring Algorithm Section 3.6.2 of 00446 // Graph Coloring in Optimization Revisited by 00447 // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen 00448 VectorBool forbidden(n); 00449 for(j = 0; j < n; j++) 00450 { // initial all colors as ok for this column 00451 for(k = 0; k < n; k++) 00452 forbidden[k] = false; 00453 00454 // for each row connected to column j 00455 itr_i = q[j].begin(); 00456 while( itr_i != q[j].end() ) 00457 { i = *itr_i++; 00458 // for each column connected to row i 00459 itr_j = p[i].begin(); 00460 while( itr_j != p[i].end() ) 00461 { // if this is not j, forbid it 00462 k = *itr_j++; 00463 forbidden[ color[k] ] = (k != j); 00464 } 00465 } 00466 k = 0; 00467 while( forbidden[k] && k < n ) 00468 { k++; 00469 CPPAD_ASSERT_UNKNOWN( k < n ); 00470 } 00471 color[j] = k; 00472 } 00473 size_t n_color = 1; 00474 for(k = 0; k < n; k++) 00475 n_color = std::max(n_color, color[k] + 1); 00476 00477 // direction vector for calls to forward 00478 VectorBase dx(n); 00479 00480 // location for return values from Reverse 00481 VectorBase dy(m); 00482 00483 // loop over colors 00484 size_t c; 00485 for(c = 0; c < n_color; c++) 00486 { // determine all the colums with this color 00487 for(j = 0; j < n; j++) 00488 { if( color[j] == c ) 00489 dx[j] = one; 00490 else dx[j] = zero; 00491 } 00492 // call forward mode for all these columns at once 00493 dy = Forward(1, dx); 00494 00495 // set the corresponding components of the result 00496 for(j = 0; j < n; j++) if( color[j] == c ) 00497 { itr_i = q[j].begin(); 00498 while( itr_i != q[j].end() ) 00499 { i = *itr_i++; 00500 jac[i * n + j] = dy[i]; 00501 } 00502 } 00503 } 00504 } 00505 else 00506 { // use reverse mode ---------------------------------------- 00507 00508 // initial coloring 00509 SizeVector color(m); 00510 for(i = 0; i < m; i++) 00511 color[i] = i; 00512 00513 // See GreedyPartialD2Coloring Algorithm Section 3.6.2 of 00514 // Graph Coloring in Optimization Revisited by 00515 // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen 00516 VectorBool forbidden(m); 00517 for(i = 0; i < m; i++) 00518 { // initial all colors as ok for this row 00519 for(k = 0; k < m; k++) 00520 forbidden[k] = false; 00521 00522 // for each column connected to row i 00523 itr_j = p[i].begin(); 00524 while( itr_j != p[i].end() ) 00525 { j = *itr_j++; 00526 // for each row connected to column j 00527 itr_i = q[j].begin(); 00528 while( itr_i != q[j].end() ) 00529 { // if this is not i, forbid it 00530 k = *itr_i++; 00531 forbidden[ color[k] ] = (k != i); 00532 } 00533 } 00534 k = 0; 00535 while( forbidden[k] && k < m ) 00536 { k++; 00537 CPPAD_ASSERT_UNKNOWN( k < n ); 00538 } 00539 color[i] = k; 00540 } 00541 size_t n_color = 1; 00542 for(k = 0; k < m; k++) 00543 n_color = std::max(n_color, color[k] + 1); 00544 00545 // weight vector for calls to reverse 00546 VectorBase w(m); 00547 00548 // location for return values from Reverse 00549 VectorBase dw(n); 00550 00551 // loop over colors 00552 size_t c; 00553 for(c = 0; c < n_color; c++) 00554 { // determine all the rows with this color 00555 for(i = 0; i < m; i++) 00556 { if( color[i] == c ) 00557 w[i] = one; 00558 else w[i] = zero; 00559 } 00560 // call reverse mode for all these rows at once 00561 dw = Reverse(1, w); 00562 00563 // set the corresponding components of the result 00564 for(i = 0; i < m; i++) if( color[i] == c ) 00565 { itr_j = p[i].begin(); 00566 while( itr_j != p[i].end() ) 00567 { j = *itr_j++; 00568 jac[i * n + j] = dw[j]; 00569 } 00570 } 00571 } 00572 } 00573 } 00574 00575 /*! 00576 Compute a sparse Jacobian. 00577 00578 The C++ source code corresponding to this operation is 00579 \verbatim 00580 jac = SparseJacobian(x, p) 00581 \endverbatim 00582 00583 \tparam Base 00584 is the base type for the recording that is stored in this 00585 ADFun<Base object. 00586 00587 \tparam VectorBase 00588 is a simple vector class with elements of type \a Base. 00589 00590 \tparam VectorSet 00591 is a simple vector class with elements of type 00592 \c bool or \c std::set<size_t>. 00593 00594 \param x 00595 is a vector specifing the point at which to compute the Jacobian. 00596 00597 \param p 00598 is the sparsity pattern for the Jacobian that we are calculating. 00599 00600 \return 00601 Will be a vector if size \c m * n containing the Jacobian at the 00602 specified point (in row major order). 00603 */ 00604 template <class Base> 00605 template <class VectorBase, class VectorSet> 00606 VectorBase ADFun<Base>::SparseJacobian( 00607 const VectorBase& x, const VectorSet& p 00608 ) 00609 { size_t m = Range(); 00610 size_t n = Domain(); 00611 VectorBase jac(m * n); 00612 00613 typedef typename VectorSet::value_type Set_type; 00614 SparseJacobianCase( Set_type(), x, p, jac); 00615 00616 return jac; 00617 } 00618 00619 /*! 00620 Compute a sparse Jacobian. 00621 00622 The C++ source code corresponding to this operation is 00623 \verbatim 00624 jac = SparseJacobian(x) 00625 \endverbatim 00626 00627 \tparam Base 00628 is the base type for the recording that is stored in this 00629 ADFun<Base object. 00630 00631 \tparam VectorBase 00632 is a simple vector class with elements of the \a Base. 00633 00634 \param x 00635 is a vector specifing the point at which to compute the Jacobian. 00636 00637 \return 00638 Will be a vector of size \c m * n containing the Jacobian at the 00639 specified point (in row major order). 00640 */ 00641 template <class Base> 00642 template <class VectorBase> 00643 VectorBase ADFun<Base>::SparseJacobian( const VectorBase& x ) 00644 { typedef CppAD::vectorBool VectorBool; 00645 00646 size_t m = Range(); 00647 size_t n = Domain(); 00648 00649 // sparsity pattern for Jacobian 00650 VectorBool p(m * n); 00651 00652 // return result 00653 VectorBase jac(m * n); 00654 00655 if( n <= m ) 00656 { size_t j, k; 00657 00658 // use forward mode 00659 VectorBool r(n * n); 00660 for(j = 0; j < n; j++) 00661 { for(k = 0; k < n; k++) 00662 r[j * n + k] = false; 00663 r[j * n + j] = true; 00664 } 00665 p = ForSparseJac(n, r); 00666 } 00667 else 00668 { size_t i, k; 00669 00670 // use reverse mode 00671 VectorBool s(m * m); 00672 for(i = 0; i < m; i++) 00673 { for(k = 0; k < m; k++) 00674 s[i * m + k] = false; 00675 s[i * m + i] = true; 00676 } 00677 p = RevSparseJac(m, s); 00678 } 00679 bool set_type = true; // only used to dispatch compiler to proper case 00680 SparseJacobianCase(set_type, x, p, jac); 00681 return jac; 00682 } 00683 00684 00685 CPPAD_END_NAMESPACE 00686 # endif