CppAD: A C++ Algorithmic Differentiation Package
20130918
|
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