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