CppAD: A C++ Algorithmic Differentiation Package  20130918
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-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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines