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