CppAD: A C++ Algorithmic Differentiation Package 20110419
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-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 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 $xref/glossary/AD Function/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 $xref/glossary/Sparsity Pattern/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 $xref/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 $xref/ForSparseJac/VectorSet/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 $xref/ForSparseJac/VectorSet/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 $xref/glossary/Sparsity Pattern/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 $xref/ForSparseJac/VectorSet/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 $xref/glossary/Sparsity Pattern/sparsity pattern/$$ 
00147 for the matrix $latex S(x)$$.
00148 
00149 $head VectorSet$$
00150 The type $icode VectorSet$$ must be a $xref/SimpleVector/$$ class with
00151 $xref/SimpleVector/Elements of Specified Type/elements of 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 $xref/ForSparseJac.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 CPPAD_BEGIN_NAMESPACE
00177 /*!
00178 \file for_sparse_jac.hpp
00179 Forward mode Jacobian sparsity patterns.
00180 */
00181 
00182 /*!
00183 Calculate Jacobian sparsity patterns using forward mode.
00184 
00185 The C++ source code corresponding to this operation is
00186 \verbatim
00187         s = f.ForSparseJac(q, r)
00188 \endverbatim
00189 
00190 \tparam Base
00191 is the base type for this recording.
00192 
00193 \tparam VectorSet
00194 is a simple vector class with elements of type \c bool.
00195 
00196 \tparam Sparsity
00197 is either \c sparse_pack or \c sparse_set. 
00198 
00199 \param q
00200 is the number of columns in the matrix \f$ R \f$.
00201 
00202 \param r
00203 is a sparsity pattern for the matrix \f$ R \f$.
00204 
00205 \param s
00206 The input value of \a s must be a vector with size \c m*q
00207 where \c m is the number of dependent variables
00208 corresponding to the operation sequence stored in \a play. 
00209 The input value of the components of \c s does not matter.
00210 On output, \a s is the sparsity pattern for the matrix
00211 \f[
00212         S(x) = F^{(1)} (x) * R
00213 \f]
00214 where \f$ F \f$ is the function corresponding to the operation sequence
00215 and \a x is any argument value.
00216 
00217 \param total_num_var
00218 is the total number of variable in this recording.
00219 
00220 \param dep_taddr
00221 maps dependendent variable index
00222 to the corresponding variable in the tape.
00223 
00224 \param ind_taddr
00225 maps independent variable index
00226 to the corresponding variable in the tape.
00227 
00228 \param play
00229 is the recording that defines the function we are computing the sparsity 
00230 pattern for.
00231 
00232 \param for_jac_sparsity
00233 the input value of \a for_jac_sparsity does not matter.
00234 On output, \a for_jac_sparsity.n_set() == \a total_num_var
00235 and \a for_jac_sparsity.end() == \a q.
00236 It contains the forward sparsity pattern for all of the variables on the
00237 tape (given the sparsity pattern for the independent variables is \f$ R \f$).
00238 */
00239 template <class Base, class VectorSet, class Sparsity> 
00240 void ForSparseJacBool(
00241         size_t                 q                , 
00242         const VectorSet&       r                ,
00243         VectorSet&             s                ,
00244         size_t                 total_num_var    ,
00245         CppAD::vector<size_t>& dep_taddr        ,
00246         CppAD::vector<size_t>& ind_taddr        ,
00247         CppAD::player<Base>&   play             ,
00248         Sparsity&              for_jac_sparsity )
00249 {
00250         // temporary indices
00251         size_t i, j;
00252 
00253         // check VectorSet is Simple Vector class with bool elements
00254         CheckSimpleVector<bool, VectorSet>();
00255 
00256         // range and domain dimensions for F
00257         size_t m = dep_taddr.size();
00258         size_t n = ind_taddr.size();
00259 
00260         CPPAD_ASSERT_KNOWN(
00261                 q > 0,
00262                 "ForSparseJac: q (first arugment) is not greater than zero"
00263         );
00264 
00265         CPPAD_ASSERT_KNOWN(
00266                 r.size() == n * q,
00267                 "ForSparseJac: r (second argument) length is not equal to\n"
00268                 "q (first argument) times domain dimension for ADFun object."
00269         );
00270 
00271         // allocate memory for the requested sparsity calculation
00272         for_jac_sparsity.resize(total_num_var, q);
00273 
00274         // set values corresponding to independent variables
00275         for(i = 0; i < n; i++)
00276         {       CPPAD_ASSERT_UNKNOWN( ind_taddr[i] < total_num_var );
00277                 // ind_taddr[i] is operator taddr for i-th independent variable
00278                 CPPAD_ASSERT_UNKNOWN( play.GetOp( ind_taddr[i] ) == InvOp );
00279 
00280                 // set bits that are true
00281                 for(j = 0; j < q; j++) if( r[ i * q + j ] )
00282                         for_jac_sparsity.add_element( ind_taddr[i], j);
00283         }
00284 
00285         // evaluate the sparsity patterns
00286         ForJacSweep(
00287                 n,
00288                 total_num_var,
00289                 &play,
00290                 for_jac_sparsity
00291         );
00292 
00293         // return values corresponding to dependent variables
00294         CPPAD_ASSERT_UNKNOWN( s.size() == m * q );
00295         for(i = 0; i < m; i++)
00296         {       CPPAD_ASSERT_UNKNOWN( dep_taddr[i] < total_num_var );
00297 
00298                 // extract the result from for_jac_sparsity
00299                 for(j = 0; j < q; j++)
00300                         s[ i * q + j ] = false;
00301                 CPPAD_ASSERT_UNKNOWN( for_jac_sparsity.end() == q );
00302                 for_jac_sparsity.begin( dep_taddr[i] );
00303                 j = for_jac_sparsity.next_element();
00304                 while( j < q )
00305                 {       s[i * q + j ] = true;
00306                         j = for_jac_sparsity.next_element();
00307                 }
00308         }
00309 }
00310 
00311 /*!
00312 Calculate Jacobian sparsity patterns using forward mode.
00313 
00314 The C++ source code corresponding to this operation is
00315 \verbatim
00316         s = f.ForSparseJac(q, r)
00317 \endverbatim
00318 
00319 \tparam Base
00320 is the base type for this recording.
00321 
00322 \tparam VectorSet
00323 is a simple vector class with elements of the type \c std::set<size_t>.
00324 
00325 \tparam Sparsity
00326 is either \c sparse_pack or \c sparse_set. 
00327 
00328 \param q
00329 is the number of columns in the matrix \f$ R \f$.
00330 
00331 \param r
00332 is a sparsity pattern for the matrix \f$ R \f$.
00333 
00334 \param s
00335 The input value of \a s must be a vector with size \c m
00336 where \c m is the number of dependent variables
00337 corresponding to the operation sequence stored in \a play. 
00338 On input, each element of \a s must be an empty set.
00339 On output, \a s is the sparsity pattern for the matrix
00340 \f[
00341         S(x) = F^{(1)} (x) * R
00342 \f]
00343 where \f$ F \f$ is the function corresponding to the operation sequence
00344 and \a x is any argument value.
00345 
00346 \param total_num_var
00347 is the total number of variable in this recording.
00348 
00349 \param dep_taddr
00350 maps dependendent variable index
00351 to the corresponding variable in the tape.
00352 
00353 \param ind_taddr
00354 maps independent variable index
00355 to the corresponding variable in the tape.
00356 
00357 \param play
00358 is the recording that defines the function we are computing the sparsity 
00359 pattern for.
00360 
00361 \param for_jac_sparsity
00362 the input value of \a for_jac_sparsity does not matter.
00363 On output, \a for_jac_sparsity.n_set() == \a total_num_var
00364 and \a for_jac_sparsity.end() == \a q.
00365 It contains the forward sparsity pattern for all of the variables on the
00366 tape (given the sparsity pattern for the independent variables is \f$ R \f$).
00367 */
00368 
00369 template <class Base, class VectorSet, class Sparsity> 
00370 void ForSparseJacSet(
00371         size_t                  q                , 
00372         const VectorSet&        r                ,
00373         VectorSet&              s                ,
00374         size_t                  total_num_var    ,
00375         CppAD::vector<size_t>&  dep_taddr        ,
00376         CppAD::vector<size_t>&  ind_taddr        ,
00377         CppAD::player<Base>&    play             ,
00378         Sparsity&               for_jac_sparsity )
00379 {
00380         // temporary indices
00381         size_t i, j;
00382         std::set<size_t>::const_iterator itr;
00383 
00384         // check VectorSet is Simple Vector class with sets for elements
00385         static std::set<size_t> two, three;
00386         if( two.empty() )
00387         {       two.insert(2);
00388                 three.insert(3);
00389         }
00390         CPPAD_ASSERT_UNKNOWN( two.size() == 1 );
00391         CPPAD_ASSERT_UNKNOWN( three.size() == 1 );
00392         CheckSimpleVector<std::set<size_t>, VectorSet>(two, three);
00393 
00394         // range and domain dimensions for F
00395         size_t m = dep_taddr.size();
00396         size_t n = ind_taddr.size();
00397 
00398         CPPAD_ASSERT_KNOWN(
00399                 q > 0,
00400                 "ForSparseJac: q (first arugment) is not greater than zero"
00401         );
00402 
00403         CPPAD_ASSERT_KNOWN(
00404                 r.size() == n,
00405                 "ForSparseJac: r (second argument) length is not equal to\n"
00406                 "the domain dimension for ADFun object."
00407         );
00408 
00409         // allocate memory for the requested sparsity calculation
00410         for_jac_sparsity.resize(total_num_var, q);
00411 
00412         // set values corresponding to independent variables
00413         for(i = 0; i < n; i++)
00414         {       CPPAD_ASSERT_UNKNOWN( ind_taddr[i] < total_num_var );
00415                 // ind_taddr[i] is operator taddr for i-th independent variable
00416                 CPPAD_ASSERT_UNKNOWN( play.GetOp( ind_taddr[i] ) == InvOp );
00417 
00418                 // add the elements that are present
00419                 itr = r[i].begin();
00420                 while( itr != r[i].end() )
00421                 {       j = *itr++;
00422                         CPPAD_ASSERT_KNOWN(
00423                                 j < q,
00424                                 "ForSparseJac: an element of the set r[i] "
00425                                 "has value greater than or equal q."
00426                         );
00427                         for_jac_sparsity.add_element( ind_taddr[i], j);
00428                 }
00429         }
00430 
00431         // evaluate the sparsity patterns
00432         ForJacSweep(
00433                 n,
00434                 total_num_var,
00435                 &play,
00436                 for_jac_sparsity
00437         );
00438 
00439         // return values corresponding to dependent variables
00440         CPPAD_ASSERT_UNKNOWN( s.size() == m );
00441         for(i = 0; i < m; i++)
00442         {       CPPAD_ASSERT_UNKNOWN( dep_taddr[i] < total_num_var );
00443                 CPPAD_ASSERT_UNKNOWN( s[i].empty() );
00444 
00445                 // extract results from for_jac_sparsity
00446                 // and add corresponding elements to sets in s
00447                 CPPAD_ASSERT_UNKNOWN( for_jac_sparsity.end() == q );
00448                 for_jac_sparsity.begin( dep_taddr[i] );
00449                 j = for_jac_sparsity.next_element();
00450                 while( j < q )
00451                 {       s[i].insert(j);
00452                         j = for_jac_sparsity.next_element();
00453                 }
00454         }
00455 }
00456 
00457 
00458 /*!
00459 User API for Jacobian sparsity patterns using forward mode.
00460 
00461 The C++ source code corresponding to this operation is
00462 \verbatim
00463         s = f.ForSparseJac(q, r)
00464 \endverbatim
00465 
00466 \tparam Base
00467 is the base type for this recording.
00468 
00469 \tparam VectorSet
00470 is a simple vector with elements of type \c bool
00471 or \c std::set<size_t>.
00472 
00473 \param q
00474 is the number of columns in the matrix \f$ R \f$.
00475 
00476 \param r
00477 is a sparsity pattern for the matrix \f$ R \f$.
00478 
00479 \return
00480 If \c VectorSet::value_type is \c bool,
00481 the return value \c s is a vector with size \c m*q
00482 where \c m is the number of dependent variables
00483 corresponding to the operation sequence stored in \c f. 
00484 If \c VectorSet::value_type is \c std::set<size_t>,
00485 the return value \c s is a vector of sets with size \c m
00486 and with all its elements between zero and \a q - 1.
00487 The value of \a s is the sparsity pattern for the matrix
00488 \f[
00489         S(x) = F^{(1)} (x) * R
00490 \f]
00491 where \f$ F \f$ is the function corresponding to the operation sequence
00492 and \a x is any argument value.
00493 
00494 \par Side Effects
00495 If \c VectorSet::value_type is \c bool,
00496 the forward sparsity pattern for all of the variables on the
00497 tape is stored in \c for_jac_sparse_pack__.
00498 In this case 
00499 \verbatim
00500         for_jac_sparse_pack_.n_set() == total_num_var_
00501         for_jac_sparse_pack_.end() == q
00502         for_jac_sparse_set_.n_set()  == 0
00503         for_jac_sparse_set_.end()  == 0
00504 \endverbatim
00505 \n
00506 \n
00507 If \c VectorSet::value_type is \c std::set<size_t>,
00508 the forward sparsity pattern for all of the variables on the
00509 tape is stored in \c for_jac_sparse_set__.
00510 In this case 
00511 \verbatim
00512         for_jac_sparse_set_.n_set()   == total_num_var_
00513         for_jac_sparse_set_.end()   == q
00514         for_jac_sparse_pack_.n_set()  == 0
00515         for_jac_sparse_pack_.end()  == 0
00516 \endverbatim
00517 */
00518 template <class Base>
00519 template <class VectorSet>
00520 VectorSet ADFun<Base>::ForSparseJac(
00521         size_t             q             , 
00522         const VectorSet&   r             )
00523 {       VectorSet s;
00524         typedef typename VectorSet::value_type Set_type; 
00525 
00526         ForSparseJacCase(
00527                 Set_type()  ,
00528                 q           ,
00529                 r           ,
00530                 s
00531         );
00532 
00533         return s;
00534 }
00535 
00536 /*! 
00537 Private helper function for ForSparseJac(q, r).
00538 
00539 All of the description in the public member function ForSparseJac(q, r) 
00540 applies. 
00541 
00542 \param set_type
00543 is a \c bool value. This argument is used to dispatch to the proper source
00544 code depending on the value of \c VectorSet::value_type.
00545 
00546 \param q
00547 See \c ForSparseJac(q, r).
00548 
00549 \param r
00550 See \c ForSparseJac(q, r).
00551 
00552 \param s
00553 is the return value for the corresponding call to \c ForSparseJac(q, r).
00554 */
00555 
00556 template <class Base>
00557 template <class VectorSet>
00558 void ADFun<Base>::ForSparseJacCase(
00559         bool                set_type      ,
00560         size_t              q             ,
00561         const VectorSet&    r             ,
00562         VectorSet&          s             )
00563 {       size_t m = Range();
00564 
00565         // dimension size of result vector
00566         s.resize( m * q );
00567 
00568         // free any results stored in for_jac_sparse_set_       
00569         for_jac_sparse_set_.resize(0, 0);
00570 
00571         // store results in r and for_jac_sparse_pack_
00572         ForSparseJacBool(
00573                 q                , 
00574                 r                ,
00575                 s                ,
00576                 total_num_var_   ,
00577                 dep_taddr_       ,
00578                 ind_taddr_       ,
00579                 play_            ,
00580                 for_jac_sparse_pack_ 
00581         );
00582 }
00583 
00584 
00585 /*! 
00586 Private helper function for \c ForSparseJac(q, r).
00587 
00588 All of the description in the public member function \c ForSparseJac(q, r) 
00589 applies. 
00590 
00591 \param set_type
00592 is a \c std::set<size_t> object. 
00593 This argument is used to dispatch to the proper source
00594 code depending on the value of \c VectorSet::value_type.
00595 
00596 \param q
00597 See \c ForSparseJac(q, r).
00598 
00599 \param r
00600 See \c ForSparseJac(q, r).
00601 
00602 \param s
00603 is the return value for the corresponding call to \c ForSparseJac(q, r).
00604 */
00605 template <class Base>
00606 template <class VectorSet>
00607 void ADFun<Base>::ForSparseJacCase(
00608         const std::set<size_t>&    set_type      ,
00609         size_t                     q             ,
00610         const VectorSet&           r             ,
00611         VectorSet&                 s             )
00612 {       size_t m = Range();
00613 
00614         // dimension size of result vector
00615         s.resize( m );
00616 
00617         // free any results stored in for_jac_sparse_set_       
00618         for_jac_sparse_pack_.resize(0, 0);
00619 
00620         // store results in r and for_jac_sparse_pack_
00621         CppAD::ForSparseJacSet(
00622                 q                , 
00623                 r                ,
00624                 s                ,
00625                 total_num_var_   ,
00626                 dep_taddr_       ,
00627                 ind_taddr_       ,
00628                 play_            ,
00629                 for_jac_sparse_set_ 
00630         );
00631 }
00632 
00633 
00634 CPPAD_END_NAMESPACE
00635 # endif