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