CppAD: A C++ Algorithmic Differentiation Package  20130918
link_sparse_hessian.cpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 /* --------------------------------------------------------------------------
00003 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-14 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 \{
00146 \file link_sparse_hessian.cpp
00147 Defines and implement sparse Hessian speed link to package specific code.
00148 */
00149 namespace {
00150      using CppAD::vector;
00151 
00152      /*!
00153      Class used by choose_row_col to determin order of row and column indices
00154      */
00155      class Key {
00156      public:
00157           /// row index
00158           size_t row_;
00159           /// column index
00160           size_t col_;
00161           /// default constructor
00162           Key(void)
00163           { }
00164           /*!
00165           Construct from a value for row and col
00166 
00167           \param row
00168           row value for this key
00169 
00170           \param col
00171           column value for this key
00172           */
00173           Key(size_t row, size_t col)
00174           : row_(row), col_(col)
00175           { }
00176           /*!
00177           Compare this key with another key using < operator
00178 
00179           \param other
00180           the other key.
00181           */
00182           bool operator<(const Key& other) const
00183           {    if( row_ == other.row_ )
00184                     return col_ < other.col_;
00185                return row_ < other.row_;
00186           }
00187      };
00188 
00189      /*!
00190      Function that randomly choose the row and column indices
00191 
00192      \param n [in]
00193      is the dimension of the argument space for the function f(x).
00194 
00195      \param row [out]
00196      the input size and elements of \c row do not matter.
00197      Upon return it is the chosen row indices.
00198 
00199      \param col [out]
00200      the input size and elements of \c col do not matter.
00201      Upon return it is the chosen column indices.
00202      */
00203      void choose_row_col(
00204           size_t          n   ,
00205           vector<size_t>& row , 
00206           vector<size_t>& col )
00207      {    size_t r, c, k, K = 5 * n;
00208 
00209           // get the random indices
00210           vector<double>  random(2 * K);
00211           CppAD::uniform_01(2 * K, random);
00212 
00213           // sort the temporary row and colunn choices
00214           vector<Key>   keys(K + n);
00215           vector<size_t> ind(K + n);
00216           for(k = 0; k < K; k++)
00217           {    r = size_t( n * random[k] );
00218                r = std::min(n-1, r);
00219                //
00220                c = size_t( n * random[k + K] );
00221                c = std::min(n-1, c);
00222                //
00223                // force to lower triangle
00224                if( c > r )
00225                     std::swap(r, c);
00226                //
00227                keys[k] = Key(r, c);
00228           }
00229           // include the diagonal
00230           for(k = 0; k < n; k++)
00231                keys[k + K] = Key(k, k);  
00232           CppAD::index_sort(keys, ind);
00233 
00234           // remove duplicates while setting the return value for row and col
00235           row.resize(0);
00236           col.resize(0);
00237           size_t r_previous = keys[ ind[0] ].row_;
00238           size_t c_previous = keys[ ind[0] ].col_;
00239           row.push_back(r_previous);
00240           col.push_back(c_previous);
00241           for(k = 1; k < K; k++)
00242           {    r = keys[ ind[k] ].row_;
00243                c = keys[ ind[k] ].row_;
00244                if( r != r_previous || c != c_previous)
00245                {    row.push_back(r);
00246                     col.push_back(c);
00247                }
00248                r_previous = r;
00249                c_previous = c;
00250           }
00251      }
00252 }
00253 
00254 /*!
00255 Package specific implementation of a sparse Hessian claculation.
00256 
00257 \param size [in]
00258 is the size of the domain space; i.e. specifies \c n.
00259 
00260 \param repeat [in]
00261 number of times tha the test is repeated.
00262 
00263 \param x [out]
00264 is a vector of size \c n containing
00265 the argument at which the Hessian was evaluated during the last repetition.
00266 
00267 \param row [in]
00268 is the row indices correpsonding to non-zero Hessian entries.
00269 
00270 \param col [in]
00271 is the column indices corresponding to non-zero Hessian entries.
00272 
00273 \param hessian [out]
00274 is a vector with size <code>n * n</code> 
00275 containing the value of the Hessian of f(x) 
00276 corresponding to the last repetition.
00277 
00278 \return
00279 is true, if the sparse Hessian speed test is implemented for this package,
00280 and false otherwise.
00281 */
00282 extern bool link_sparse_hessian(
00283      size_t                           size      ,
00284      size_t                           repeat    ,
00285      const CppAD::vector<size_t>&     row       ,
00286      const CppAD::vector<size_t>&     col       , 
00287            CppAD::vector<double>&     x         ,
00288            CppAD::vector<double>&     hessian
00289 );
00290 
00291 /*!
00292 Is sparse Hessian test avaialable.
00293 
00294 \return
00295 true, if spare Hessian available for this package, and false otherwise.
00296 */
00297 bool available_sparse_hessian(void)
00298 {    
00299      size_t n      = 2;
00300      size_t repeat = 1;
00301      vector<double> x(n);
00302      vector<double> hessian(n * n);
00303      vector<size_t> row, col; 
00304      choose_row_col(n, row, col);
00305 
00306      return link_sparse_hessian(n, repeat, row, col, x, hessian);
00307      exit(0);
00308 }
00309 /*!
00310 Does final sparse Hessian value pass correctness test.
00311 
00312 \param is_package_double
00313 if true, we are checking function values instead of derivatives.
00314 
00315 \return
00316 true, if correctness test passes, and false otherwise.
00317 */
00318 bool correct_sparse_hessian(bool is_package_double)
00319 {    
00320      size_t n      = 10;
00321      size_t repeat = 1;
00322      vector<double> x(n);
00323      vector<double> hessian(n * n);
00324      vector<size_t> row, col;
00325      choose_row_col(n, row, col);
00326 
00327      link_sparse_hessian(n, repeat, row, col, x, hessian);
00328 
00329      size_t order, size;
00330      if( is_package_double)
00331      {    order = 0;  // check function value
00332           size  = 1;
00333      }
00334      else
00335      {    order = 2;     // check hessian value
00336           size  = n * n;
00337      }
00338      CppAD::vector<double> check(size);
00339      CppAD::sparse_hes_fun<double>(n, x, row, col, order, check);
00340      bool ok = true;
00341      size_t k;
00342      for( k = 0; k < size; k++)
00343           ok &= CppAD::NearEqual(check[k], hessian[k], 1e-10, 1e-10);
00344 
00345      return ok;
00346 }
00347 
00348 /*!
00349 Sparse Hessian speed test.
00350 
00351 \param size
00352 is the dimension of the argument space for this speed test.
00353 
00354 \param repeat
00355 is the number of times to repeate the speed test.
00356 */
00357 void speed_sparse_hessian(size_t size, size_t repeat)
00358 {    size_t n = size;    
00359      vector<double> x(n);
00360      vector<double> hessian(n * n);
00361      vector<size_t> row, col;
00362      choose_row_col(n, row, col);
00363 
00364      // note that cppad/sparse_hessian.cpp assumes that x.size() == size
00365      link_sparse_hessian(n, repeat, row, col, x, hessian);
00366      return;
00367 }
00368 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines