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