CppAD: A C++ Algorithmic Differentiation Package  20130102
sparse_hessian.hpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 # ifndef CPPAD_SPARSE_HESSIAN_INCLUDED
00003 # define CPPAD_SPARSE_HESSIAN_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 sparse_hessian$$
00018 $spell
00019      recomputed
00020      CppAD
00021      valarray
00022      std
00023      Bool
00024      hes
00025      const
00026      Taylor
00027 $$
00028 
00029 $section Sparse Hessian: Easy Driver$$
00030 $index SparseHessian$$
00031 $index hessian, sparse$$
00032 
00033 $head Syntax$$
00034 $icode%hes% = %f%.SparseHessian(%x%, %w%)
00035 %hes% = %f%.SparseHessian(%x%, %w%, %p%)
00036 %n_sweep% = %f%.SparseHessian(%x%, %w%, %p%, %row%, %col%, %hes%, %work%)
00037 %$$
00038 
00039 $head Purpose$$
00040 We use $latex n$$ for the $cref/domain/seq_property/Domain/$$ size,
00041 and $latex m$$ for the $cref/range/seq_property/Range/$$ size of $icode f$$.
00042 We use $latex F : \B{R}^n \rightarrow \B{R}^m$$ do denote the
00043 $cref/AD function/glossary/AD Function/$$
00044 corresponding to $icode f$$. 
00045 The syntax above sets $icode hes$$ to the Hessian 
00046 $latex \[
00047      H(x) = \dpow{2}{x} \sum_{i=1}^m w_i F_i (x) 
00048 \] $$
00049 This routine takes advantage of the sparsity of the Hessian
00050 in order to reduce the amount of computation necessary.
00051 If $icode row$$ and $icode col$$ are present, it also takes
00052 advantage of the reduced set of elements of the Hessian that
00053 need to be computed.
00054 One can use speed tests (e.g. $cref speed_test$$)
00055 to verify that results are computed faster
00056 than when using the routine $cref Hessian$$.
00057 
00058 $head f$$
00059 The object $icode f$$ has prototype
00060 $codei%
00061      ADFun<%Base%> %f%
00062 %$$
00063 Note that the $cref ADFun$$ object $icode f$$ is not $code const$$
00064 (see $cref/Uses Forward/sparse_hessian/Uses Forward/$$ below).
00065 
00066 $head x$$
00067 The argument $icode x$$ has prototype
00068 $codei%
00069      const %VectorBase%& %x%
00070 %$$
00071 (see $cref/VectorBase/sparse_hessian/VectorBase/$$ below)
00072 and its size 
00073 must be equal to $icode n$$, the dimension of the
00074 $cref/domain/seq_property/Domain/$$ space for $icode f$$.
00075 It specifies
00076 that point at which to evaluate the Hessian.
00077 
00078 $head w$$
00079 The argument $icode w$$ has prototype
00080 $codei%
00081      const %VectorBase%& %w%
00082 %$$
00083 and size $latex m$$.
00084 It specifies the value of $latex w_i$$ in the expression 
00085 for $icode hes$$.
00086 The more components of $latex w$$ that are identically zero,
00087 the more sparse the resulting Hessian may be (and hence the more efficient
00088 the calculation of $icode hes$$ may be).
00089 
00090 $head p$$
00091 The argument $icode p$$ is optional and has prototype
00092 $codei%
00093      const %VectorSet%& %p%
00094 %$$
00095 (see $cref/VectorSet/sparse_hessian/VectorSet/$$ below)
00096 If it has elements of type $code bool$$,
00097 its size is $latex n * n$$.
00098 If it has elements of type $code std::set<size_t>$$,
00099 its size is $latex n$$ and all its set elements are between
00100 zero and $latex n - 1$$.
00101 It specifies a
00102 $cref/sparsity pattern/glossary/Sparsity Pattern/$$ 
00103 for the Hessian $latex H(x)$$.
00104 $pre
00105 
00106 $$
00107 If this sparsity pattern does not change between calls to 
00108 $codei SparseHessian$$, it should be faster to calculate $icode p$$ once and
00109 pass this argument to $codei SparseHessian$$.
00110 In addition,
00111 if you specify $icode p$$, CppAD will use the same
00112 type of sparsity representation
00113 (vectors of $code bool$$ or vectors of $code std::set<size_t>$$)
00114 for its internal calculations.
00115 Otherwise, the representation
00116 for the internal calculations is unspecified.
00117 
00118 $head row, col$$
00119 The arguments $icode row$$ and $icode col$$ are optional and have prototype
00120 $codei%
00121      const %VectorSize%& %row%
00122      const %VectorSize%& %col%
00123 %$$
00124 (see $cref/VectorSize/sparse_hessian/VectorSize/$$ below).
00125 They specify which rows and columns of $latex H (x)$$ are
00126 returned and in what order.
00127 We use $latex K$$ to denote the value $icode%hes%.size()%$$
00128 which must also equal the size of $icode row$$ and $icode col$$.
00129 Furthermore,
00130 for $latex k = 0 , \ldots , K-1$$, it must hold that
00131 $latex row[k] < n$$ and $latex col[k] < n$$.
00132 In addition, 
00133 all of the $latex (row[k], col[k])$$ pairs must correspond to a true value
00134 in the sparsity pattern $icode p$$.
00135 
00136 $head hes$$
00137 The result $icode hes$$ has prototype
00138 $codei%
00139      %VectorBase% %hes%
00140 %$$
00141 In the case where $icode row$$ and $icode col$$ are not present,
00142 the size of $icode hes$$ is $latex n * n$$ and
00143 its size is $latex n * n$$.
00144 In this case, for $latex i = 0 , \ldots , n - 1 $$ 
00145 and $latex ell = 0 , \ldots , n - 1$$
00146 $latex \[
00147      hes [ j * n + \ell ] = \DD{ w^{\rm T} F }{ x_j }{ x_\ell } ( x )
00148 \] $$
00149 $pre
00150 
00151 $$
00152 In the case where the arguments $icode row$$ and $icode col$$ are present,
00153 we use $latex K$$ to denote the size of $icode hes$$. 
00154 The input value of its elements does not matter.
00155 Upon return, for $latex k = 0 , \ldots , K - 1$$,
00156 $latex \[
00157      hes [ k ] = \DD{ w^{\rm T} F }{ x_j }{ x_\ell } (x)
00158      \; , \;
00159      \; {\rm where} \;
00160      j = row[k]
00161      \; {\rm and } \;
00162      \ell = col[k]
00163 \] $$
00164 
00165 $head work$$
00166 If this argument is present, it has prototype
00167 $codei%
00168      sparse_hessian_work& %work%
00169 %$$
00170 This object can only be used with the routines $code SparseHessian$$.
00171 During its the first use, information is stored in $icode work$$.
00172 This is used to reduce the work done by future calls to $code SparseHessian$$ 
00173 with the same $icode f$$, $icode p$$, $icode row$$, and $icode col$$.
00174 If a future call is make where any of these values have changed,
00175 you must first call $icode%work%.clear()%$$
00176 to inform CppAD that this information needs to be recomputed.
00177 
00178 $head n_sweep$$
00179 The return value $icode n_sweep$$ has prototype
00180 $codei%
00181      size_t %n_sweep%
00182 %$$
00183 It is the number of first order forward sweeps 
00184 used to compute the requested Hessian values.
00185 Each first forward sweep is followed by a second order reverse sweep 
00186 so it is also the number of reverse sweeps.
00187 This is proportional to the total work that $code SparseHessian$$ does, 
00188 not counting the zero order forward sweep, 
00189 or the work to combine multiple columns into a single 
00190 forward-reverse sweep pair.
00191 
00192 $head VectorBase$$
00193 The type $icode VectorBase$$ must be a $cref SimpleVector$$ class with
00194 $cref/elements of type/SimpleVector/Elements of Specified Type/$$
00195 $icode Base$$.
00196 The routine $cref CheckSimpleVector$$ will generate an error message
00197 if this is not the case.
00198 
00199 $head VectorSet$$
00200 The type $icode VectorSet$$ must be a $cref SimpleVector$$ class with
00201 $cref/elements of type/SimpleVector/Elements of Specified Type/$$
00202 $code bool$$ or $code std::set<size_t>$$;
00203 see $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for a discussion
00204 of the difference.
00205 The routine $cref CheckSimpleVector$$ will generate an error message
00206 if this is not the case.
00207 
00208 $subhead Restrictions$$
00209 If $icode VectorSet$$ has elements of $code std::set<size_t>$$,
00210 then $icode%p%[%i%]%$$ must return a reference (not a copy) to the 
00211 corresponding set.
00212 According to section 26.3.2.3 of the 1998 C++ standard,
00213 $code std::valarray< std::set<size_t> >$$ does not satisfy
00214 this condition. 
00215 
00216 $head VectorSize$$
00217 The type $icode VectorSize$$ must be a $cref SimpleVector$$ class with
00218 $cref/elements of type/SimpleVector/Elements of Specified Type/$$
00219 $code size_t$$.
00220 The routine $cref CheckSimpleVector$$ will generate an error message
00221 if this is not the case.
00222 
00223 $head Uses Forward$$
00224 After each call to $cref Forward$$,
00225 the object $icode f$$ contains the corresponding 
00226 $cref/Taylor coefficients/glossary/Taylor Coefficient/$$.
00227 After a call to any of the sparse Hessian routines,
00228 the zero order Taylor coefficients correspond to
00229 $icode%f%.Forward(0, %x%)%$$
00230 and the other coefficients are unspecified.
00231 
00232 $head Example$$
00233 $children%
00234      example/sparse_hessian.cpp
00235 %$$
00236 The routine
00237 $cref sparse_hessian.cpp$$ 
00238 is examples and tests of $code sparse_hessian$$.
00239 It return $code true$$, if it succeeds and $code false$$ otherwise.
00240 
00241 
00242 $end
00243 -----------------------------------------------------------------------------
00244 */
00245 # include <cppad/local/std_set.hpp>
00246 
00247 CPPAD_BEGIN_NAMESPACE
00248 /*!
00249 \defgroup sparse_hessian_hpp sparse_hessian.hpp
00250 \{
00251 \file sparse_hessian.hpp
00252 Sparse Hessian driver routine and helper functions.
00253 */
00254 // ===========================================================================
00255 /*!
00256 class used by SparseHessian to hold information 
00257 so it does not need to be recomputed.
00258 */
00259 class sparse_hessian_work {
00260      public:
00261           /// version of user r array sorted by row or column
00262           CppAD::vector<size_t> r_sort;
00263           /// version of user c array sorted by row or column
00264           CppAD::vector<size_t> c_sort;
00265           /// mapping from sorted array indices to user array indices
00266           CppAD::vector<size_t> k_sort;
00267           /// results of the coloring algorithm
00268           CppAD::vector<size_t> color;
00269           /// inform CppAD that this information needs to be recomputed
00270           void clear(void)
00271           {    r_sort.clear();
00272                c_sort.clear();
00273                k_sort.clear();
00274                color.clear();
00275           }
00276 };
00277 // ===========================================================================
00278 /*!
00279 Private helper function that does computation for all Sparse Hessian cases.
00280 
00281 \tparam Base
00282 See \c SparseHessian(x, w, p, row, col, hes, work).
00283 
00284 \tparam VectorBase
00285 See \c SparseHessian(x, w, p, row, col, hes, work).
00286 
00287 \tparam VectorSet
00288 is either \c sparse_pack, \c sparse_set or \c sparse_list.
00289 
00290 \param x
00291 See \c SparseHessian(x, w, p, row, col, hes, work).
00292 
00293 \param w
00294 See \c SparseHessian(x, w, p, row, col, hes, work).
00295 
00296 \param sparsity
00297 If <code>work.color.size() != 0</code>, then \c sparsity is not used.
00298 Otherwise, it is a 
00299 sparsity pattern for the Hessian of this ADFun<Base> object.
00300 
00301 \param hes
00302 See \c SparseHessian(x, w, p, row, col, hes, work).
00303 
00304 \param work
00305 See \c SparseHessian(x, w, p, row, col, hes, work).
00306 
00307 \return
00308 See \c SparseHessian(x, w, p, row, col, hes, work).
00309 */
00310 template<class Base>
00311 template <class VectorBase, class VectorSet>
00312 size_t ADFun<Base>::SparseHessianCompute(
00313      const VectorBase&     x           ,
00314      const VectorBase&     w           ,
00315      VectorSet&            sparsity    ,
00316      VectorBase&           hes         ,
00317      sparse_hessian_work& work         )
00318 {
00319      using   CppAD::vectorBool;
00320      size_t i, j, k, ell;
00321 
00322      CppAD::vector<size_t>& row(work.r_sort);
00323      CppAD::vector<size_t>& col(work.c_sort);
00324      CppAD::vector<size_t>& user_k(work.k_sort);
00325      CppAD::vector<size_t>& color(work.color);
00326 
00327      size_t n = Domain();
00328 
00329      // some values
00330      const Base zero(0);
00331      const Base one(1);
00332 
00333      // check VectorBase is Simple Vector class with Base type elements
00334      CheckSimpleVector<Base, VectorBase>();
00335 
00336      CPPAD_ASSERT_UNKNOWN( size_t(x.size()) == n );
00337      CPPAD_ASSERT_UNKNOWN( color.size() == 0 || color.size() == n );
00338 
00339      // number of components of Hessian that are required
00340      size_t K = hes.size();
00341      CPPAD_ASSERT_UNKNOWN( row.size() == K+1 );
00342      CPPAD_ASSERT_UNKNOWN( col.size() == K );
00343      CPPAD_ASSERT_UNKNOWN( row[K] == n );
00344 
00345      // Point at which we are evaluating the Hessian
00346      Forward(0, x);
00347 
00348      // Rows of the Hessian (i below) correspond to the forward mode index
00349      // and columns (j below) correspond to the reverse mode index.
00350      if( color.size() == 0 )
00351      {
00352           CPPAD_ASSERT_UNKNOWN( sparsity.n_set() ==  n );
00353           CPPAD_ASSERT_UNKNOWN( sparsity.end() ==  n );
00354 
00355           // rows and columns that are in the returned hessian
00356           VectorSet r_used, c_used;
00357           r_used.resize(n, n);
00358           c_used.resize(n, n);
00359           k = 0;
00360           while( k < K )
00361           {    CPPAD_ASSERT_UNKNOWN( row[k] < n && col[k] < n );
00362                CPPAD_ASSERT_UNKNOWN( k == 0 || row[k-1] <= row[k] );
00363                CPPAD_ASSERT_KNOWN(
00364                     sparsity.is_element(row[k], col[k]) ,
00365                     "SparseHessian: an (row, col) pair is not in sparsity pattern."
00366                );
00367                r_used.add_element(col[k], row[k]);
00368                c_used.add_element(row[k], col[k]);
00369                k++;
00370           }
00371      
00372           // given a column index, which rows are non-zero and not used
00373           VectorSet not_used;
00374           not_used.resize(n, n);
00375           for(i = 0; i < n; i++)
00376           {    sparsity.begin(i);
00377                j = sparsity.next_element();
00378                while( j != sparsity.end() )
00379                {    if( ! r_used.is_element(j , i) )
00380                          not_used.add_element(j, i);
00381                     j = sparsity.next_element();
00382                }
00383           }
00384 
00385           // initial coloring
00386           color.resize(n);
00387           for(i = 0; i < n; i++)
00388                color[i] = i;
00389      
00390           // See GreedyPartialD2Coloring Algorithm Section 3.6.2 of
00391           // Graph Coloring in Optimization Revisited by
00392           // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen
00393           vectorBool forbidden(n);
00394           for(i = 1; i < n; i++)
00395           {
00396                // initial all colors as ok for this row
00397                // (value of forbidden for ell > i does not matter)
00398                for(ell = 0; ell <= i; ell++)
00399                     forbidden[ell] = false;
00400      
00401                // -----------------------------------------------------
00402                // Forbid colors that this row would destroy results for.
00403                // for each column that is non-zero for this row
00404                sparsity.begin(i);
00405                j = sparsity.next_element();
00406                while( j != sparsity.end() )
00407                {    // for each row that this column uses
00408                     r_used.begin(j);
00409                     ell = r_used.next_element();
00410                     while( ell != r_used.end() )
00411                     {    // if this is not the same row, forbid its color
00412                          if( ell < i )
00413                               forbidden[ color[ell] ] = true;
00414                          ell = r_used.next_element();
00415                     }
00416                     j = sparsity.next_element();
00417                }
00418      
00419                // -------------------------------------------------------
00420                // Forbid colors that would destroy the results for this row.
00421                // for each column that this row used
00422                c_used.begin(i);
00423                j = c_used.next_element();
00424                while( j != c_used.end() )
00425                {    // For each row that is non-zero for this column
00426                     // (the used rows have already been checked above).
00427                     not_used.begin(j);
00428                     ell = not_used.next_element();
00429                     while( ell != not_used.end() )
00430                     {    // if this is not the same row, forbid its color
00431                          if( ell < i )
00432                               forbidden[ color[ell] ] = true;
00433                          ell = not_used.next_element();
00434                     }
00435                     j = c_used.next_element();
00436                }
00437 
00438                // pick the color with the smallest index
00439                ell = 0;
00440                while( forbidden[ell] )
00441                {    ell++;
00442                     CPPAD_ASSERT_UNKNOWN( ell <= i );
00443                }
00444                color[i] = ell;
00445           }
00446      }
00447      size_t n_color = 1;
00448      for(ell = 0; ell < n; ell++) 
00449           n_color = std::max(n_color, color[ell] + 1);
00450 
00451      // direction vector for calls to forward (rows of the Hessian)
00452      VectorBase u(n);
00453 
00454      // location for return values from reverse (columns of the Hessian)
00455      VectorBase ddw(2 * n);
00456 
00457      // initialize the return value
00458      for(k = 0; k < K; k++)
00459           hes[k] = zero;
00460 
00461      // loop over colors
00462      size_t n_sweep = 0;
00463      for(ell = 0; ell < n_color; ell++)
00464      {    bool any = false;
00465           k = 0;
00466           for(i = 0; i < n; i++) if( color[i] == ell )
00467           {    // find first k such that row[k] has color ell
00468                if( ! any )
00469                {    while( row[k] < i )
00470                          k++;
00471                     any = row[k] == i;
00472                }
00473           }
00474           if( any )
00475           {    n_sweep++;
00476                // combine all rows with this color
00477                for(i = 0; i < n; i++)
00478                {    u[i] = zero;
00479                     if( color[i] == ell )
00480                          u[i] = one;
00481                }
00482                // call forward mode for all these rows at once
00483                Forward(1, u);
00484 
00485                // evaluate derivative of w^T * F'(x) * u
00486                ddw = Reverse(2, w);
00487 
00488                // set the corresponding components of the result
00489                for(i = 0; i < n; i++) if( color[i] == ell )
00490                {    // find first index in c for this column
00491                     while( row[k] < i )
00492                          k++;
00493                     // extract the results for this row
00494                     while( row[k] == i ) 
00495                     {    hes[ user_k[k] ] = ddw[ col[k] * 2 + 1 ];
00496                          k++;
00497                     }
00498                }
00499           }
00500      }
00501      return n_sweep;
00502 }
00503 // ===========================================================================
00504 /*!
00505 Private helper function for vector of \c bool sparsity pattern cases.
00506 
00507 \tparam Base
00508 See \c SparseHessian(x, w, p, row, col, hes, work).
00509 
00510 \tparam VectorBase
00511 See \c SparseHessian(x, w, p, row, col, hes, work).
00512 
00513 \tparam VectorSet
00514 is a simple vector with elements of type \c bool.
00515 
00516 \param set_type
00517 has element type for vector representing the sparsity sets.
00518 
00519 \param x
00520 See \c SparseHessian(x, w, p, row, col, hes, work).
00521 
00522 \param w
00523 See \c SparseHessian(x, w, p, row, col, hes, work).
00524 
00525 \param p
00526 Sparsity pattern for the Hessian of this ADFun<Base> object.
00527 
00528 \param hes
00529 See \c SparseHessian(x, w, p, row, col, hes, work).
00530 
00531 \param work
00532 See \c SparseHessian(x, w, p, row, col, hes, work).
00533 
00534 \return
00535 See \c SparseHessian(x, w, p, row, col, hes, work).
00536 */
00537 template <class Base>
00538 template <class VectorBase, class VectorSet>
00539 size_t ADFun<Base>::SparseHessianCase(
00540      bool                  set_type  ,
00541      const VectorBase&     x         ,
00542      const VectorBase&     w         ,
00543      const VectorSet&      p         ,
00544      VectorBase&           hes       ,
00545      sparse_hessian_work&  work      )
00546 {
00547      size_t n = Domain();
00548 
00549      // check VectorSet is Simple Vector class with bool elements
00550      CheckSimpleVector<bool, VectorSet>();
00551 
00552      // check VectorBase is Simple Vector class with Base type elements
00553      CheckSimpleVector<Base, VectorBase>();
00554 
00555      CPPAD_ASSERT_UNKNOWN( size_t(x.size()) == n );
00556      CPPAD_ASSERT_UNKNOWN( size_t(w.size()) == Range() );
00557      CPPAD_ASSERT_KNOWN(
00558           size_t(p.size()) == n * n,
00559           "SparseHessian: using bool values for sparsity and p.size() "
00560           "not equal square of domain dimension for f"
00561      );
00562  
00563      sparse_pack sparsity;
00564      if( work.color.size() == 0 )
00565      {    bool transpose = false;
00566           typedef typename VectorSet::value_type Set_type;
00567           sparsity_user2internal(sparsity, p, n, n, transpose);
00568      }
00569      
00570      // compute the Hessian
00571      size_t n_sweep = SparseHessianCompute(x, w, sparsity, hes, work);
00572 
00573      return n_sweep;
00574 }
00575 /*!
00576 Private helper function for vector of std::set<size_t> sparsity pattern cases.
00577 
00578 \tparam Base
00579 See \c SparseHessian(x, w, p, row, col, hes, work).
00580 
00581 \tparam VectorBase
00582 See \c SparseHessian(x, w, p, row, col, hes, work).
00583 
00584 \tparam VectorSet
00585 is a simple vector with elements of type <code>std::set<size_t></code>.
00586 
00587 \param set_type
00588 has element type for vector representing the sparsity sets.
00589 
00590 \param x
00591 See \c SparseHessian(x, w, p, row, col, hes, work).
00592 
00593 \param w
00594 See \c SparseHessian(x, w, p, row, col, hes, work).
00595 
00596 \param p
00597 Sparsity pattern for the Hessian of this ADFun<Base> object.
00598 
00599 \param hes
00600 See \c SparseHessian(x, w, p, row, col, hes, work).
00601 
00602 \param work
00603 See \c SparseHessian(x, w, p, row, col, hes, work).
00604 
00605 \return
00606 See \c SparseHessian(x, w, p, row, col, hes, work).
00607 */
00608 template <class Base>
00609 template <class VectorBase, class VectorSet>
00610 size_t ADFun<Base>::SparseHessianCase(
00611      const std::set<size_t>&     set_type  ,
00612      const VectorBase&           x         ,
00613      const VectorBase&           w         ,
00614      const VectorSet&            p         ,
00615      VectorBase&                 hes       ,
00616      sparse_hessian_work&        work      )
00617 {
00618      size_t n = Domain();
00619 
00620      // check VectorSet is Simple Vector class with std::set<size_t> elements
00621      CheckSimpleVector<std::set<size_t>, VectorSet>(
00622           one_element_std_set<size_t>(), two_element_std_set<size_t>()
00623      );
00624 
00625      // check VectorBase is Simple Vector class with Base type elements
00626      CheckSimpleVector<Base, VectorBase>();
00627 
00628      CPPAD_ASSERT_UNKNOWN( size_t(x.size()) == n );
00629      CPPAD_ASSERT_UNKNOWN( size_t(w.size()) == Range() );
00630      CPPAD_ASSERT_KNOWN(
00631           size_t(p.size()) == n,
00632           "SparseHessian: using std::set<size_t> for sparsity and p.size() "
00633           "not equal domain dimension for f"
00634      );
00635  
00636      CPPAD_INTERNAL_SPARSE_SET sparsity;
00637      if( work.color.size() == 0 )
00638      {    bool transpose = false;
00639           typedef typename VectorSet::value_type Set_type;
00640           sparsity_user2internal(sparsity, p, n, n, transpose);
00641      }
00642      
00643      // compute the Hessian
00644      size_t n_sweep = SparseHessianCompute(x, w, sparsity, hes, work);
00645 
00646      return n_sweep;
00647 }
00648 // ===========================================================================
00649 /*!
00650 Private helper function for SparseHessian(x, w, p).
00651 
00652 All of the description in the public member function SparseHessian(x, w, p)
00653 applies.
00654 
00655 \param set_type
00656 is a \c bool value. This argument is used to dispatch to the proper source
00657 code depending on the value of \c VectorSet::value_type.
00658 
00659 \param x
00660 See \c SparseHessian(x, w, p).
00661 
00662 \param w
00663 See \c SparseHessian(x, w, p).
00664 
00665 \param p
00666 See \c SparseHessian(x, w, p).
00667 
00668 \param hes
00669 is the return value for the corresponding call to \c SparseHessian(x, w, p).
00670 On input, it must have size equal to the domain times range dimension
00671 for this ADFun<Base> object.
00672 On return, it will contain the Hessian.
00673 */
00674 
00675 template <class Base>
00676 template <class VectorBase, class VectorSet>
00677 void ADFun<Base>::SparseHessianCase(
00678      bool               set_type         ,
00679      const VectorBase&  x                ,
00680      const VectorBase&  w                ,
00681      const VectorSet&   p                ,
00682      VectorBase&        hes              )
00683 {
00684      size_t n = Domain();
00685      size_t i, j, k;
00686 
00687      // check Vector is Simple Vector class with bool elements
00688      CheckSimpleVector<bool, VectorSet>();
00689 
00690      // check Vector is Simple Vector class with Base type elements
00691      CheckSimpleVector<Base, VectorBase>();
00692 
00693      CPPAD_ASSERT_KNOWN(
00694           size_t(x.size()) == n,
00695           "SparseHessian: size of x not equal domain dimension for f"
00696      ); 
00697      CPPAD_ASSERT_KNOWN(
00698           size_t(w.size()) == Range(),
00699           "SparseHessian: size of w not equal range dimension for f"
00700      );
00701      CPPAD_ASSERT_KNOWN(
00702           size_t(p.size()) == n * n,
00703           "SparseHessian: using bool values and size of p "
00704           "not equal square of domain dimension for f"
00705      );
00706      sparse_hessian_work work;
00707      CppAD::vector<size_t>& row(work.r_sort);
00708      CppAD::vector<size_t>& col(work.c_sort);
00709      CppAD::vector<size_t>& user_k(work.k_sort);
00710 
00711      k = 0;
00712      for(i = 0; i < n; i++)
00713      {    for(j = 0; j < n; j++)
00714                if( p[i * n + j] )
00715                     k++;
00716      }
00717      size_t K = k;
00718      VectorBase H(K);
00719      row.resize(K+1);
00720      col.resize(K);
00721      user_k.resize(K);
00722      k = 0;
00723      for(i = 0; i < n; i++)
00724      {    for(j = 0; j < n; j++)
00725                if( p[i * n + j] )
00726                {    row[k]      = i;
00727                     col[k]      = j;
00728                     user_k[k] = k;  
00729                     k++;
00730                }
00731      }
00732      row[K] = n;
00733  
00734      SparseHessianCase(set_type, x, w, p, H, work);
00735 
00736      Base zero(0);
00737      for(i = 0; i < n; i++)
00738      {    for(j = 0; j < n; j++)
00739                hes[i * n + j] = zero; 
00740      }
00741      for(k = 0; k < K; k++)
00742           hes[row[k] * n + col[k]] = H[k]; 
00743 }
00744 /*!
00745 Private helper function for SparseHessian(x, w, p).
00746 
00747 All of the description in the public member function SparseHessian(x, w, p)
00748 applies.
00749 
00750 \param set_type
00751 is a \c std::set<size_t> value.
00752 This argument is used to dispatch to the proper source
00753 code depending on the value of \c VectorSet::value_type.
00754 
00755 \param x
00756 See \c SparseHessian(x, w, p).
00757 
00758 \param w
00759 See \c SparseHessian(x, w, p).
00760 
00761 \param p
00762 See \c SparseHessian(x, w, p).
00763 
00764 \param hes
00765 is the return value for the corresponding call to \c SparseHessian(x, w, p).
00766 On input, it must have size equal to the domain times range dimension
00767 for this ADFun<Base> object.
00768 On return, it will contain the Hessian.
00769 */
00770 
00771 template <class Base>
00772 template <class VectorBase, class VectorSet>
00773 void ADFun<Base>::SparseHessianCase(
00774      const std::set<size_t>&  set_type         ,
00775      const VectorBase&        x                ,
00776      const VectorBase&        w                ,
00777      const VectorSet&         p                ,
00778      VectorBase&              hes              )
00779 {
00780      size_t n = Domain();
00781      size_t i, j, k;
00782 
00783      // check VectorSet is Simple Vector class with sets for elements
00784      CheckSimpleVector<std::set<size_t>, VectorSet>(
00785           one_element_std_set<size_t>(), two_element_std_set<size_t>()
00786      );
00787 
00788      // check Vector is Simple Vector class with Base type elements
00789      CheckSimpleVector<Base, VectorBase>();
00790 
00791      CPPAD_ASSERT_KNOWN(
00792           size_t(x.size()) == n,
00793           "SparseHessian: size of x not equal domain dimension for f"
00794      ); 
00795      CPPAD_ASSERT_KNOWN(
00796           size_t(w.size()) == Range(),
00797           "SparseHessian: size of w not equal range dimension for f"
00798      );
00799      CPPAD_ASSERT_KNOWN(
00800           size_t(p.size()) == n,
00801           "SparseHessian: using size_t sets for sparsity pattern and p.size() "
00802           "not equal domain dimension for f"
00803      );
00804      sparse_hessian_work work;
00805      CppAD::vector<size_t>& row(work.r_sort);
00806      CppAD::vector<size_t>& col(work.c_sort);
00807      CppAD::vector<size_t>& user_k(work.k_sort);
00808 
00809      k = 0;
00810      std::set<size_t>::const_iterator itr;
00811      for(i = 0; i < n; i++)
00812      {    itr = p[i].begin();
00813           while( itr != p[i].end() )
00814           {    itr++;
00815                k++;
00816           }
00817      }
00818      size_t K = k;
00819      VectorBase H(K);
00820      row.resize(K+1);
00821      col.resize(K);
00822      user_k.resize(K);
00823      k = 0;
00824      for(i = 0; i < n; i++)
00825      {    itr = p[i].begin();
00826           while( itr != p[i].end() )
00827           {    row[k]      = i;
00828                col[k]      = *itr;
00829                user_k[k] = k;  
00830                itr++;
00831                k++;
00832           }
00833      }
00834      row[K] = n;
00835  
00836      SparseHessianCase(set_type, x, w, p, H, work);
00837 
00838      Base zero(0);
00839      for(i = 0; i < n; i++)
00840      {    for(j = 0; j < n; j++)
00841                hes[i * n + j] = zero; 
00842      }
00843      for(k = 0; k < K; k++)
00844           hes[row[k] * n + col[k]] = H[k]; 
00845 }
00846 // ===========================================================================
00847 // Public Member Functions
00848 // ===========================================================================
00849 /*!
00850 Compute user specified subset of a sparse Hessian.
00851 
00852 The C++ source code corresponding to this operation is
00853 \verbatim
00854      SparceHessian(x, w, p, row, col, hes, work)
00855 \endverbatim
00856 
00857 \tparam Base
00858 is the base type for the recording that is stored in this ADFun<Base object.
00859 
00860 \tparam VectorBase
00861 is a simple vector class with elements of type \a Base.
00862 
00863 \tparam VectorSet
00864 is a simple vector class with elements of type 
00865 \c bool or \c std::set<size_t>.
00866 
00867 \tparam VectorSize
00868 is a simple vector class with elements of type \c size_t.
00869 
00870 \param x
00871 is a vector specifing the point at which to compute the Hessian.
00872 
00873 \param w
00874 is the weighting vector that defines a scalar valued function by
00875 a weighted sum of the components of the vector valued function
00876 $latex F(x)$$.
00877 
00878 \param p
00879 is the sparsity pattern for the Hessian that we are calculating.
00880 
00881 \param row
00882 is the vector of row indices for the returned Hessian values.
00883 
00884 \param col
00885 is the vector of columns indices for the returned Hessian values.
00886 It must have the same size are r.
00887 
00888 \param hes
00889 is the vector of Hessian values.
00890 It must have the same size are r. 
00891 The return value <code>hes[k]</code> is the second partial of 
00892 \f$ w^{\rm T} F(x)\f$ with respect to the
00893 <code>row[k]</code> and <code>col[k]</code> component of \f$ x\f$.
00894 
00895 \param work
00896 contains information that depends on the function object, sparsity pattern,
00897 \c row, and \c col vector.
00898 If these values are the same, \c work does not need to be recomputed.
00899 To be more specific, 
00900 \c r_sort is sorted copy of \c row ,
00901 \c c_sort is sorted copy of \c col ,
00902 <code>k_sort[k]</code> is the original index corresponding to the
00903 values <code>r_sort[k]</code> and <code>c_sort[k]</code>.
00904 The order for the sort is by columns.
00905 Let \c n the domain dimension,
00906 and \c K the size of \c row , \c col , and \c hes.
00907 There is one extra entry 
00908 in the sorted row array and it has value <code>r_sort[K]=n</code>.
00909 The \c color vector is set and used by \c SparseHessianCompute.
00910 
00911 \return
00912 Is the number of first order forward sweeps used to compute the
00913 requested Hessian values.
00914 (This is also equal to the number of second order reverse sweeps.)
00915 The total work, not counting the zero order
00916 forward sweep, or the time to combine computations, is proportional to this
00917 return value.
00918 */
00919 template<class Base>
00920 template <class VectorBase, class VectorSet, class VectorSize>
00921 size_t ADFun<Base>::SparseHessian(
00922      const VectorBase&     x    ,
00923      const VectorBase&     w    ,
00924      const VectorSet&      p    ,
00925      const VectorSize&     row  ,
00926      const VectorSize&     col  ,
00927      VectorBase&           hes  ,
00928      sparse_hessian_work&  work )
00929 {
00930      size_t n = Domain();
00931      size_t k, K = hes.size();
00932      if( work.r_sort.size() == 0 )
00933      {    // create version of (row, col, k) sorted by row value
00934           work.c_sort.resize(K);
00935           work.r_sort.resize(K+1);
00936           work.k_sort.resize(K);
00937 
00938           // put sorting indices in k_sort
00939           index_sort(row, work.k_sort);
00940 
00941           for(k = 0; k < K; k++)
00942           {    work.r_sort[k] = row[ work.k_sort[k] ];
00943                work.c_sort[k] = col[ work.k_sort[k] ];
00944           }
00945           work.r_sort[K] = n;
00946      }
00947 # ifndef NDEBUG
00948      CPPAD_ASSERT_KNOWN(
00949           size_t(x.size()) == n ,
00950           "SparseHessian: size of x not equal domain dimension for f."
00951      ); 
00952      CPPAD_ASSERT_KNOWN(
00953           size_t(row.size()) == K && size_t(col.size()) == K ,
00954           "SparseHessian: either r or c does not have the same size as ehs."
00955      ); 
00956      CPPAD_ASSERT_KNOWN(
00957           work.r_sort.size() == K+1 &&
00958           work.c_sort.size() == K   &&
00959           work.k_sort.size() == K   ,
00960           "SparseHessian: invalid value in work."
00961      );
00962      CPPAD_ASSERT_KNOWN(
00963           work.color.size() == 0 || work.color.size() == n,
00964           "SparseHessian: invalid value in work."
00965      );
00966      for(k = 0; k < K; k++)
00967      {    CPPAD_ASSERT_KNOWN(
00968                row[k] < n,
00969                "SparseHessian: invalid value in r."
00970           );
00971           CPPAD_ASSERT_KNOWN(
00972                col[k] < n,
00973                "SparseHessian: invalid value in c."
00974           );
00975           CPPAD_ASSERT_KNOWN(
00976                work.k_sort[k] < K,
00977                "SparseHessian: invalid value in work."
00978           );
00979           CPPAD_ASSERT_KNOWN(
00980                work.r_sort[k] == row[ work.k_sort[k] ] ,
00981                "SparseHessian: invalid value in work."
00982           );
00983           CPPAD_ASSERT_KNOWN(
00984                work.c_sort[k] == col[ work.k_sort[k] ],
00985                "SparseHessian: invalid value in work."
00986           );
00987      }
00988      if( work.color.size() != 0 )
00989           for(size_t j = 0; j < n; j++) CPPAD_ASSERT_KNOWN(
00990                work.color[j] < n,
00991                "SparseHessian: invalid value in work."
00992      );
00993 # endif
00994  
00995      typedef typename VectorSet::value_type Set_type;
00996      size_t n_sweep = SparseHessianCase(Set_type(), x, w, p, hes, work);
00997      return n_sweep;
00998 }
00999 /*!
01000 Compute a sparse Hessian.
01001 
01002 The C++ source code coresponding to this operation is
01003 \verbatim
01004      hes = SparseHessian(x, w, p)
01005 \endverbatim
01006 
01007 
01008 \tparam Base
01009 is the base type for the recording that is stored in this
01010 ADFun<Base object.
01011 
01012 \tparam VectorBase
01013 is a simple vector class with elements of the \a Base.
01014 
01015 \tparam VectorSet
01016 is a simple vector class with elements of type
01017 \c bool or \c std::set<size_t>.
01018 
01019 \param x
01020 is a vector specifing the point at which to compute the Hessian.
01021 
01022 \param w
01023 The Hessian is computed for a weighted sum of the components
01024 of the function corresponding to this ADFun<Base> object.
01025 The argument \a w specifies the weights for each component.
01026 It must have size equal to the range dimension for this ADFun<Base> object.
01027 
01028 \param p
01029 is a sparsity pattern for the Hessian.
01030 
01031 \return
01032 Will be a vector of size \c n * n containing the Hessian of 
01033 at the point specified by \a x
01034 (where \c n is the domain dimension for this ADFun<Base> object).
01035 */
01036 template <class Base>
01037 template <class VectorBase, class VectorSet>
01038 VectorBase ADFun<Base>::SparseHessian(
01039      const VectorBase& x, const VectorBase& w, const VectorSet& p
01040 )
01041 {    size_t n = Domain();
01042      VectorBase hes(n * n);
01043 
01044      typedef typename VectorSet::value_type Set_type;
01045      SparseHessianCase(Set_type(), x, w, p, hes);
01046 
01047      return hes;
01048 }
01049 /*!
01050 Compute a sparse Hessian
01051 
01052 The C++ source code coresponding to this operation is
01053 \verbatim
01054      hes = SparseHessian(x, w)
01055 \endverbatim
01056 
01057 
01058 \tparam Base
01059 is the base type for the recording that is stored in this
01060 ADFun<Base object.
01061 
01062 \tparam VectorBase
01063 is a simple vector class with elements of the \a Base.
01064 
01065 \param x
01066 is a vector specifing the point at which to compute the Hessian.
01067 
01068 \param w
01069 The Hessian is computed for a weighted sum of the components
01070 of the function corresponding to this ADFun<Base> object.
01071 The argument \a w specifies the weights for each component.
01072 It must have size equal to the range dimension for this ADFun<Base> object.
01073 
01074 \return
01075 Will be a vector of size \c n * n containing the Hessian of 
01076 at the point specified by \a x
01077 (where \c n is the domain dimension for this ADFun<Base> object).
01078 */
01079 template <class Base>
01080 template <class VectorBase>
01081 VectorBase ADFun<Base>::SparseHessian(const VectorBase &x, const VectorBase &w)
01082 {    size_t i, j, k;
01083      typedef CppAD::vectorBool VectorBool;
01084 
01085      size_t m = Range();
01086      size_t n = Domain();
01087 
01088      // determine the sparsity pattern p for Hessian of w^T F
01089      VectorBool r(n * n);
01090      for(j = 0; j < n; j++)
01091      {    for(k = 0; k < n; k++)
01092                r[j * n + k] = false;
01093           r[j * n + j] = true;
01094      }
01095      ForSparseJac(n, r);
01096      //
01097      VectorBool s(m);
01098      for(i = 0; i < m; i++)
01099           s[i] = w[i] != 0;
01100      VectorBool p = RevSparseHes(n, s);
01101 
01102      // compute sparse Hessian
01103      VectorBase hes(n * n);
01104      bool set_type = true; // only used to dispatch compiler to proper case
01105      SparseHessianCase(set_type, x, w, p, hes);
01106 
01107      return hes;
01108 }
01109 
01110 /*! \} */
01111 CPPAD_END_NAMESPACE
01112 # endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines