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