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