CppAD: A C++ Algorithmic Differentiation Package
20130102
|
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 /*! \} */