CppAD: A C++ Algorithmic Differentiation Package  20130102
link_sparse_hessian.cpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 /* --------------------------------------------------------------------------
00003 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
00004 
00005 CppAD is distributed under multiple licenses. This distribution is under
00006 the terms of the 
00007                     Eclipse Public License Version 1.0.
00008 
00009 A copy of this license is included in the COPYING file of this distribution.
00010 Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
00011 -------------------------------------------------------------------------- */
00012 
00013 /*
00014 $begin link_sparse_hessian$$
00015 $spell
00016      const
00017      bool
00018      CppAD
00019 $$
00020 
00021 $index link_sparse_hessian$$
00022 $index sparse, speed test$$
00023 $index speed, test sparse$$
00024 $index test, sparse speed$$
00025 
00026 $section Speed Testing Sparse Hessian$$
00027 
00028 $head Prototype$$
00029 $codei%extern bool link_sparse_hessian(
00030      size_t                        %size%      ,
00031      size_t                        %repeat%    ,
00032      CppAD::vector<double>&        %x%         ,
00033      const CppAD::vector<size_t>&  %row%       ,
00034      const CppAD::vector<size_t>&  %col%       , 
00035      CppAD::vector<double>&        %hessian%
00036 );
00037 %$$
00038 
00039 $head Method$$
00040 Given a row index vector $latex row$$
00041 and a second column vector $latex col$$,
00042 the corresponding function 
00043 $latex f : \B{R}^n \rightarrow \B{R} $$ 
00044 is defined by $cref sparse_hes_fun$$.
00045 The non-zero entries in the Hessian of this function have 
00046 one of the following forms:
00047 $latex \[
00048      \DD{f}{x[row[k]]}{x[row[k]]}
00049      \; , \;
00050      \DD{f}{x[row[k]]}{x[col[k]]}
00051      \; , \;
00052      \DD{f}{x[col[k]]}{x[row[k]]}
00053      \; , \;
00054      \DD{f}{x[col[k]]}{x[col[k]]}
00055 \] $$
00056 for some $latex k $$ between zero and $latex K-1 $$.
00057 All the other terms of the Hessian are zero.
00058 
00059 $head size$$
00060 The argument $icode size$$, referred to as $latex n$$ below,
00061 is the dimension of the domain space for $latex f(x)$$.
00062 
00063 $head repeat$$
00064 The argument $icode repeat$$ is the number of times
00065 to repeat the test
00066 (with a different value for $icode x$$ corresponding to
00067 each repetition).
00068 
00069 $head x$$
00070 The argument $icode x$$ has prototype
00071 $codei%
00072         CppAD::vector<double> &%x%
00073 %$$
00074 and its size is $latex n$$; i.e., $icode%x%.size() == %size%$$.
00075 The input value of the elements of $icode x$$ does not matter.
00076 On output, it has been set to the
00077 argument value for which the function,
00078 or its derivative, is being evaluated.
00079 The value of this vector need not change with each repetition.
00080 
00081 $head row$$
00082 The argument $icode row$$ has prototype
00083 $codei%
00084      const CppAD::vector<size_t> %row%
00085 %$$
00086 Its size defines the value $latex K$$.
00087 It contains the row indices for the corresponding function $latex f(x)$$.
00088 All the elements of $icode row$$ are between zero and $latex n-1$$.
00089 
00090 $head col$$
00091 The argument $icode col$$ has prototype
00092 $codei%
00093      const CppAD::vector<size_t> %col%
00094 %$$
00095 Its size must be the same as $icode row$$; i.e., $latex K$$.
00096 It contains the column indices for the corresponding function 
00097 $latex f(x)$$.
00098 All the elements of $icode col$$ are between zero and $latex n-1$$.
00099 $pre
00100 
00101 $$
00102 There are no duplicate row and column entires; i.e., if $icode%j% != %k%$$,
00103 $codei%
00104      %row%[%j%] != %row%[%k%] || %col%[%j%] != %col%[%k%]
00105 %$$
00106 Only the lower triangle of the Hessian is included in the indices; i.e.,
00107 $codei%row%[%k%] >= %col%[%k%]%$$.
00108 Furthermore, for all the non-zero entries in the lower triangle
00109 are included; i.e., if $latex i \geq j$$ and
00110 $latex \DD{f}{x[i]}{x[j]} \neq 0$$,
00111 there is an index $icode k$$ such that
00112 $icode%i% = %row%[%k%]%$$ and
00113 $icode%j% = %col%[%k%]%$$.
00114 
00115 $head hessian$$
00116 The argument $icode hessian$$ has prototype
00117 $codei%
00118      CppAD::vector<double>  &hessian
00119 %$$
00120 and its size is $latex n \times n$$.
00121 The input value of its elements does not matter. 
00122 The output value of its elements is the Hessian of the function $latex f(x)$$.
00123 To be more specific, for
00124 $latex i = 0 , \ldots , n-1$$,
00125 $latex j = 0 , \ldots , n-1$$,
00126 $latex \[
00127      \DD{f}{x[i]}{x[j]} (x) = hessian [ i * n + j ]
00128 \] $$
00129 
00130 $subhead double$$
00131 In the case where $icode package$$ is $code double$$,
00132 only the first element of $icode hessian$$ is used and it is actually 
00133 the value of $latex f(x)$$ (derivatives are not computed).
00134 
00135 $end 
00136 -----------------------------------------------------------------------------
00137 */
00138 # include <cppad/vector.hpp>
00139 # include <cppad/near_equal.hpp>
00140 # include <cppad/speed/sparse_hes_fun.hpp>
00141 # include <cppad/speed/uniform_01.hpp>
00142 # include <cppad/index_sort.hpp>
00143 
00144 /*!
00145 \defgroup link_sparse_hessian_cpp link_sparse_hessian.cpp
00146 \{
00147 \file link_sparse_hessian.cpp
00148 Defines and implement sparse Hessian speed link to package specific code.
00149 */
00150 namespace {
00151      using CppAD::vector;
00152 
00153      /*!
00154      Class used by choose_row_col to determin order of row and column indices
00155      */
00156      class Key {
00157      public:
00158           /// row index
00159           size_t row_;
00160           /// column index
00161           size_t col_;
00162           /// default constructor
00163           Key(void)
00164           { }
00165           /*!
00166           Construct from a value for row and col
00167 
00168           \param row
00169           row value for this key
00170 
00171           \param col
00172           column value for this key
00173           */
00174           Key(size_t row, size_t col)
00175           : row_(row), col_(col)
00176           { }
00177           /*!
00178           Compare this key with another key using < operator
00179 
00180           \param other
00181           the other key.
00182           */
00183           bool operator<(const Key& other) const
00184           {    if( row_ == other.row_ )
00185                     return col_ < other.col_;
00186                return row_ < other.row_;
00187           }
00188      };
00189 
00190      /*!
00191      Function that randomly choose the row and column indices
00192 
00193      \param n [in]
00194      is the dimension of the argument space for the function f(x).
00195 
00196      \param row [out]
00197      the input size and elements of \c row do not matter.
00198      Upon return it is the chosen row indices.
00199 
00200      \param col [out]
00201      the input size and elements of \c col do not matter.
00202      Upon return it is the chosen column indices.
00203      */
00204      void choose_row_col(
00205           size_t          n   ,
00206           vector<size_t>& row , 
00207           vector<size_t>& col )
00208      {    size_t r, c, k, K = 5 * n;
00209 
00210           // get the random indices
00211           vector<double>  random(2 * K);
00212           CppAD::uniform_01(2 * K, random);
00213 
00214           // sort the temporary row and colunn choices
00215           vector<Key>   keys(K + n);
00216           vector<size_t> ind(K + n);
00217           for(k = 0; k < K; k++)
00218           {    r = size_t( n * random[k] );
00219                r = std::min(n-1, r);
00220                //
00221                c = size_t( n * random[k + K] );
00222                c = std::min(n-1, c);
00223                //
00224                // force to lower triangle
00225                if( c > r )
00226                     std::swap(r, c);
00227                //
00228                keys[k] = Key(r, c);
00229           }
00230           // include the diagonal
00231           for(k = 0; k < n; k++)
00232                keys[k + K] = Key(k, k);  
00233           CppAD::index_sort(keys, ind);
00234 
00235           // remove duplicates while setting the return value for row and col
00236           row.resize(0);
00237           col.resize(0);
00238           size_t r_previous = keys[ ind[0] ].row_;
00239           size_t c_previous = keys[ ind[0] ].col_;
00240           row.push_back(r_previous);
00241           col.push_back(c_previous);
00242           for(k = 1; k < K; k++)
00243           {    r = keys[ ind[k] ].row_;
00244                c = keys[ ind[k] ].row_;
00245                if( r != r_previous || c != c_previous)
00246                {    row.push_back(r);
00247                     col.push_back(c);
00248                }
00249                r_previous = r;
00250                c_previous = c;
00251           }
00252      }
00253 }
00254 
00255 /*!
00256 Package specific implementation of a sparse Hessian claculation.
00257 
00258 \param size [in]
00259 is the size of the domain space; i.e. specifies \c n.
00260 
00261 \param repeat [in]
00262 number of times tha the test is repeated.
00263 
00264 \param x [out]
00265 is a vector of size \c n containing
00266 the argument at which the Hessian was evaluated during the last repetition.
00267 
00268 \param row [in]
00269 is the row indices correpsonding to non-zero Hessian entries.
00270 
00271 \param col [in]
00272 is the column indices corresponding to non-zero Hessian entries.
00273 
00274 \param hessian [out]
00275 is a vector with size <code>n * n</code> 
00276 containing the value of the Hessian of f(x) 
00277 corresponding to the last repetition.
00278 
00279 \return
00280 is true, if the sparse Hessian speed test is implemented for this package,
00281 and false otherwise.
00282 */
00283 extern bool link_sparse_hessian(
00284      size_t                           size       ,
00285      size_t                           repeat     ,
00286      CppAD::vector<double>            &x         ,
00287      const CppAD::vector<size_t>      &row       ,
00288      const CppAD::vector<size_t>      &col       , 
00289      CppAD::vector<double>            &hessian
00290 );
00291 
00292 /*!
00293 Is sparse Hessian test avaialable.
00294 
00295 \return
00296 true, if spare Hessian available for this package, and false otherwise.
00297 */
00298 bool available_sparse_hessian(void)
00299 {    
00300      size_t n      = 2;
00301      size_t repeat = 1;
00302      vector<double> x(n);
00303      vector<double> hessian(n * n);
00304      vector<size_t> row, col; 
00305      choose_row_col(n, row, col);
00306 
00307      return link_sparse_hessian(n, repeat, x, row, col, hessian);
00308      exit(0);
00309 }
00310 /*!
00311 Does final sparse Hessian value pass correctness test.
00312 
00313 \return
00314 true, if correctness test passes, and false otherwise.
00315 */
00316 bool correct_sparse_hessian(bool is_package_double)
00317 {    
00318      size_t n      = 10;
00319      size_t repeat = 1;
00320      vector<double> x(n);
00321      vector<double> hessian(n * n);
00322      vector<size_t> row, col;
00323      choose_row_col(n, row, col);
00324 
00325      link_sparse_hessian(n, repeat, x, row, col, hessian);
00326 
00327      size_t order, size;
00328      if( is_package_double)
00329      {    order = 0;  // check function value
00330           size  = 1;
00331      }
00332      else
00333      {    order = 2;     // check hessian value
00334           size  = n * n;
00335      }
00336      CppAD::vector<double> check(size);
00337      CppAD::sparse_hes_fun<double>(n, x, row, col, order, check);
00338      bool ok = true;
00339      size_t k;
00340      for( k = 0; k < size; k++)
00341           ok &= CppAD::NearEqual(check[k], hessian[k], 1e-10, 1e-10);
00342 
00343      return ok;
00344 }
00345 
00346 /*!
00347 Sparse Hessian speed test.
00348 
00349 \param size
00350 is the dimension of the argument space for this speed test.
00351 
00352 \param repeat
00353 is the number of times to repeate the speed test.
00354 */
00355 void speed_sparse_hessian(size_t size, size_t repeat)
00356 {    size_t n = size;    
00357      vector<double> x(n);
00358      vector<double> hessian(n * n);
00359      vector<size_t> row, col;
00360      choose_row_col(n, row, col);
00361 
00362      // note that cppad/sparse_hessian.cpp assumes that x.size() == size
00363      link_sparse_hessian(n, repeat, x, row, col, hessian);
00364      return;
00365 }
00366 
00367 /*! \} */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines