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