CppAD: A C++ Algorithmic Differentiation Package
20130918
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_CHECKPOINT_INCLUDED 00003 # define CPPAD_CHECKPOINT_INCLUDED 00004 00005 /* -------------------------------------------------------------------------- 00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-14 Bradley M. Bell 00007 00008 CppAD is distributed under multiple licenses. This distribution is under 00009 the terms of the 00010 Eclipse Public License Version 1.0. 00011 00012 A copy of this license is included in the COPYING file of this distribution. 00013 Please visit http://www.coin-or.org/CppAD/ for information on other licenses. 00014 -------------------------------------------------------------------------- */ 00015 00016 namespace CppAD { // BEGIN_CPPAD_NAMESPACE 00017 /*! 00018 \file checkpoint.hpp 00019 defining checkpoint functions. 00020 */ 00021 00022 /* 00023 $begin checkpoint$$ 00024 $spell 00025 cppad.hpp 00026 CppAD 00027 checkpoint 00028 checkpointing 00029 algo 00030 afun 00031 const 00032 $$ 00033 00034 $section Checkpointing Functions$$ 00035 $index function, checkpoint$$ 00036 $index checkpoint, function$$ 00037 00038 $head Syntax$$ 00039 $codei%checkpoint<%Base%> %afun%(%name%, %algo%, %ax%, %ay%) 00040 %afun%.option(%option_value%) 00041 %algo%(%ax%, %ay%) 00042 %afun%(%ax%, %ay%) 00043 checkpoint<%Base%>::clear()%$$ 00044 00045 $head Purpose$$ 00046 You can reduce the size of the tape and memory required for AD by 00047 checkpointing functions of the form $latex y = f(x)$$ where 00048 $latex f : B^n \rightarrow B^m$$. 00049 00050 $head Method$$ 00051 The $code checkpoint$$ class is derived from $code atomic_base$$ 00052 and makes this easy. 00053 It implements all the $code atomic_base$$ 00054 $cref/virtual functions/atomic_base/Virtual Functions/$$ 00055 and hence its source code $code cppad/local/checkpoint.hpp$$ 00056 provides an example implementation of $cref atomic_base$$. 00057 The difference is that $code checkpoint.hpp$$ uses AD 00058 instead of user provided derivatives. 00059 00060 $head constructor$$ 00061 The constructor 00062 $codei% 00063 checkpoint<%Base%> %afun%(%name%, %algo%, %ax%, %ay%) 00064 %$$ 00065 cannot be called in $cref/parallel/ta_in_parallel/$$ mode. 00066 In addition, you cannot currently be recording 00067 $codei%AD<%Base%>%$$ operations when the constructor is called. 00068 This class is implemented as a derived class of 00069 $cref/atomic_base/atomic_ctor/atomic_base/$$ and hence 00070 some of its error message will refer to $code atomic_base$$. 00071 00072 $head Base$$ 00073 The type $icode Base$$ specifies the base type for AD operations. 00074 00075 $head ADVector$$ 00076 The type $icode ADVector$$ must be a 00077 $cref/simple vector class/SimpleVector/$$ with elements of type 00078 $codei%AD<%Base%>%$$. 00079 00080 $head name$$ 00081 This $icode checkpoint$$ constructor argument has prototype 00082 $codei% 00083 const char* %name% 00084 %$$ 00085 It is the name used for error reporting. 00086 The suggested value for $icode name$$ is $icode afun$$; i.e., 00087 the same name as used for the function. 00088 00089 $head ax$$ 00090 This argument has prototype 00091 $codei% 00092 const %ADVector%& %ax% 00093 %$$ 00094 and size must be equal to $icode n$$. 00095 It specifies vector $latex x \in B^n$$ 00096 at which an $codei%AD<%Base%>%$$ version of 00097 $latex y = f(x)$$ is to be evaluated. 00098 00099 $head ay$$ 00100 This argument has prototype 00101 $codei% 00102 %ADVector%& %ay% 00103 %$$ 00104 Its input size must be equal to $icode m$$ and does not change. 00105 The input values of its elements do not matter. 00106 Upon return, it is an $codei%AD<%Base%>%$$ version of 00107 $latex y = f(x)$$. 00108 00109 $head option$$ 00110 The $code option$$ syntax can be used to set the type of sparsity 00111 pattern used by $icode afun$$. 00112 This is an $codei%atomic_base<%Base%>%$$ function and its documentation 00113 can be found at $cref atomic_option$$. 00114 00115 $head algo$$ 00116 The type of $icode algo$$ is arbitrary, except for the fact that 00117 the syntax 00118 $codei% 00119 %algo%(%ax%, %ay%) 00120 %$$ 00121 must evaluate the function $latex y = f(x)$$ using 00122 $codei%AD<%Base%>%$$ operations. 00123 In addition, we assume that the 00124 $cref/operation sequence/glossary/Operation/Sequence/$$ 00125 does not depend on the value of $icode ax$$. 00126 00127 $head afun$$ 00128 Given $icode ax$$ it computes the corresponding value of $icode ay$$ 00129 using the operation sequence corresponding to $icode algo$$. 00130 If $codei%AD<%Base%>%$$ operations are being recorded, 00131 it enters the computation as single operation in the recording 00132 see $cref/start recording/Independent/Start Recording/$$. 00133 (Currently each use of $icode afun$$ actually corresponds to 00134 $icode%m%+%n%+2%$$ operations and creates $icode m$$ new variables, 00135 but this is not part of the CppAD specifications and my change.) 00136 00137 $head clear$$ 00138 The $code atomic_base$$ class holds onto static work space in order to 00139 increase speed by avoiding system memory allocation calls. 00140 This call makes to work space $cref/available/ta_available/$$ to 00141 for other uses by the same thread. 00142 This should be called when you are done using the 00143 user atomic functions for a specific value of $icode Base$$. 00144 00145 $subhead Restriction$$ 00146 The $code clear$$ routine cannot be called 00147 while in $cref/parallel/ta_in_parallel/$$ execution mode. 00148 00149 $children% 00150 example/atomic/checkpoint.cpp 00151 %$$ 00152 $head Example$$ 00153 The file $cref checkpoint.cpp$$ contains an example and test 00154 of these operations. 00155 It returns true if it succeeds and false if it fails. 00156 00157 $end 00158 */ 00159 template <class Base> 00160 class checkpoint : public atomic_base<Base> { 00161 private: 00162 ADFun<Base> f_; 00163 public: 00164 /*! 00165 Constructor of a checkpoint object 00166 00167 \param name [in] 00168 is the user's name for the AD version of this atomic operation. 00169 00170 \param algo [in/out] 00171 user routine that compute AD function values 00172 (not const because state may change during evaluation). 00173 00174 \param ax [in] 00175 argument value where algo operation sequence is taped. 00176 00177 \param ay [out] 00178 function value at specified argument value. 00179 */ 00180 template <class Algo, class ADVector> 00181 checkpoint(const char* name, 00182 Algo& algo, const ADVector& ax, ADVector& ay) 00183 : atomic_base<Base>(name) 00184 { CheckSimpleVector< CppAD::AD<Base> , ADVector>(); 00185 00186 // make a copy of ax because Independent modifies AD information 00187 ADVector x_tmp(ax); 00188 // delcare x_tmp as the independent variables 00189 Independent(x_tmp); 00190 // record mapping from x_tmp to ay 00191 algo(x_tmp, ay); 00192 // create function f_ : x -> y 00193 f_.Dependent(ay); 00194 // suppress checking for nan in f_ results 00195 // (see optimize documentation for atomic functions) 00196 f_.check_for_nan(false); 00197 // now optimize (we expect to use this function many times). 00198 f_.optimize(); 00199 } 00200 /*! 00201 Implement the user call to <tt>afun(ax, ay)</tt>. 00202 00203 \tparam ADVector 00204 A simple vector class with elements of type <code>AD<Base></code>. 00205 00206 \param id 00207 optional parameter which must be zero if present. 00208 00209 \param ax 00210 is the argument vector for this call, 00211 <tt>ax.size()</tt> determines the number of arguments. 00212 00213 \param ay 00214 is the result vector for this call, 00215 <tt>ay.size()</tt> determines the number of results. 00216 */ 00217 template <class ADVector> 00218 void operator()(const ADVector& ax, ADVector& ay, size_t id = 0) 00219 { CPPAD_ASSERT_KNOWN( 00220 id == 0, 00221 "checkpoint: id is non-zero in afun(ax, ay, id)" 00222 ); 00223 this->atomic_base<Base>::operator()(ax, ay, id); 00224 } 00225 /*! 00226 Link from user_atomic to forward mode 00227 00228 \copydetails atomic_base::forward 00229 */ 00230 virtual bool forward( 00231 size_t p , 00232 size_t q , 00233 const vector<bool>& vx , 00234 vector<bool>& vy , 00235 const vector<Base>& tx , 00236 vector<Base>& ty ) 00237 { 00238 CPPAD_ASSERT_UNKNOWN( f_.size_var() > 0 ); 00239 CPPAD_ASSERT_UNKNOWN( tx.size() % (q+1) == 0 ); 00240 CPPAD_ASSERT_UNKNOWN( ty.size() % (q+1) == 0 ); 00241 size_t n = tx.size() / (q+1); 00242 size_t m = ty.size() / (q+1); 00243 bool ok = true; 00244 size_t i, j; 00245 00246 // 2DO: test both forward and reverse vy information 00247 if( vx.size() > 0 ) 00248 { //Compute Jacobian sparsity pattern. 00249 vector< std::set<size_t> > s(m); 00250 if( n <= m ) 00251 { vector< std::set<size_t> > r(n); 00252 for(j = 0; j < n; j++) 00253 r[j].insert(j); 00254 s = f_.ForSparseJac(n, r); 00255 } 00256 else 00257 { vector< std::set<size_t> > r(m); 00258 for(i = 0; i < m; i++) 00259 r[i].insert(i); 00260 s = f_.RevSparseJac(m, r); 00261 } 00262 std::set<size_t>::const_iterator itr; 00263 for(i = 0; i < m; i++) 00264 { vy[i] = false; 00265 for(itr = s[i].begin(); itr != s[i].end(); itr++) 00266 { j = *itr; 00267 assert( j < n ); 00268 // y[i] depends on the value of x[j] 00269 vy[i] |= vx[j]; 00270 } 00271 } 00272 } 00273 ty = f_.Forward(q, tx); 00274 00275 // no longer need the Taylor coefficients in f_ 00276 // (have to reconstruct them every time) 00277 f_.capacity_order(0); 00278 return ok; 00279 } 00280 /*! 00281 Link from user_atomic to reverse mode 00282 00283 \copydetails atomic_base::reverse 00284 */ 00285 virtual bool reverse( 00286 size_t q , 00287 const vector<Base>& tx , 00288 const vector<Base>& ty , 00289 vector<Base>& px , 00290 const vector<Base>& py ) 00291 { 00292 CPPAD_ASSERT_UNKNOWN( f_.size_var() > 0 ); 00293 CPPAD_ASSERT_UNKNOWN( tx.size() % (q+1) == 0 ); 00294 CPPAD_ASSERT_UNKNOWN( ty.size() % (q+1) == 0 ); 00295 bool ok = true; 00296 00297 // put proper forward mode coefficients in f_ 00298 # ifdef NDEBUG 00299 f_.Forward(q, tx); 00300 # else 00301 size_t n = tx.size() / (q+1); 00302 size_t m = ty.size() / (q+1); 00303 CPPAD_ASSERT_UNKNOWN( px.size() == n * (q+1) ); 00304 CPPAD_ASSERT_UNKNOWN( py.size() == m * (q+1) ); 00305 size_t i, j, k; 00306 // 00307 vector<Base> check_ty = f_.Forward(q, tx); 00308 for(i = 0; i < m; i++) 00309 { for(k = 0; k <= q; k++) 00310 { j = i * (q+1) + k; 00311 CPPAD_ASSERT_UNKNOWN( check_ty[j] == ty[j] ); 00312 } 00313 } 00314 # endif 00315 // now can run reverse mode 00316 px = f_.Reverse(q+1, py); 00317 00318 // no longer need the Taylor coefficients in f_ 00319 // (have to reconstruct them every time) 00320 f_.capacity_order(0); 00321 return ok; 00322 } 00323 /*! 00324 Link from user_atomic to forward sparse Jacobian 00325 00326 \copydetails atomic_base::for_sparse_jac 00327 */ 00328 virtual bool for_sparse_jac( 00329 size_t q , 00330 const vector< std::set<size_t> >& r , 00331 vector< std::set<size_t> >& s ) 00332 { 00333 bool ok = true; 00334 s = f_.ForSparseJac(q, r); 00335 00336 // no longer need the forward mode sparsity pattern 00337 // (have to reconstruct them every time) 00338 f_.size_forward_set(0); 00339 00340 return ok; 00341 } 00342 /*! 00343 Link from user_atomic to forward sparse Jacobian 00344 00345 \copydetails atomic_base::for_sparse_jac 00346 */ 00347 virtual bool for_sparse_jac( 00348 size_t q , 00349 const vector<bool>& r , 00350 vector<bool>& s ) 00351 { 00352 bool ok = true; 00353 s = f_.ForSparseJac(q, r); 00354 00355 // no longer need the forward mode sparsity pattern 00356 // (have to reconstruct them every time) 00357 f_.size_forward_bool(0); 00358 00359 return ok; 00360 } 00361 /*! 00362 Link from user_atomic to forward sparse Jacobian 00363 00364 \copydetails atomic_base::rev_sparse_jac 00365 */ 00366 virtual bool rev_sparse_jac( 00367 size_t q , 00368 const vector< std::set<size_t> >& rt , 00369 vector< std::set<size_t> >& st ) 00370 { 00371 bool ok = true; 00372 00373 // compute rt 00374 // 2DO: remove need for nz_compare all the time. It is only really 00375 // necessary when optimizer calls this member function. 00376 bool transpose = true; 00377 bool nz_compare = true; 00378 st = f_.RevSparseJac(q, rt, transpose, nz_compare); 00379 00380 return ok; 00381 } 00382 /*! 00383 Link from user_atomic to forward sparse Jacobian 00384 00385 \copydetails atomic_base::rev_sparse_jac 00386 */ 00387 virtual bool rev_sparse_jac( 00388 size_t q , 00389 const vector<bool>& rt , 00390 vector<bool>& st ) 00391 { 00392 bool ok = true; 00393 00394 // compute rt 00395 bool transpose = true; 00396 bool nz_compare = true; 00397 // 2DO: remove need for nz_compare all the time. It is only really 00398 // necessary when optimizer calls this member function. 00399 st = f_.RevSparseJac(q, rt, transpose, nz_compare); 00400 00401 return ok; 00402 } 00403 /*! 00404 Link from user_atomic to forward sparse Jacobian 00405 00406 \copydetails atomic_base::rev_sparse_hes 00407 */ 00408 virtual bool rev_sparse_hes( 00409 const vector<bool>& vx , 00410 const vector<bool>& s , 00411 vector<bool>& t , 00412 size_t q , 00413 const vector< std::set<size_t> >& r , 00414 const vector< std::set<size_t> >& u , 00415 vector< std::set<size_t> >& v ) 00416 { size_t n = v.size(); 00417 size_t m = u.size(); 00418 CPPAD_ASSERT_UNKNOWN( r.size() == v.size() ); 00419 CPPAD_ASSERT_UNKNOWN( s.size() == m ); 00420 CPPAD_ASSERT_UNKNOWN( t.size() == n ); 00421 bool ok = true; 00422 bool transpose = true; 00423 std::set<size_t>::const_iterator itr; 00424 00425 // compute sparsity pattern for T(x) = S(x) * f'(x) 00426 t = f_.RevSparseJac(1, s); 00427 # ifndef NDEBUG 00428 for(size_t j = 0; j < n; j++) 00429 CPPAD_ASSERT_UNKNOWN( vx[j] || ! t[j] ) 00430 # endif 00431 00432 // V(x) = f'(x)^T * g''(y) * f'(x) * R + g'(y) * f''(x) * R 00433 // U(x) = g''(y) * f'(x) * R 00434 // S(x) = g'(y) 00435 00436 // compute sparsity pattern for A(x) = f'(x)^T * U(x) 00437 vector< std::set<size_t> > a(n); 00438 a = f_.RevSparseJac(q, u, transpose); 00439 00440 // set version of s 00441 vector< std::set<size_t> > set_s(1); 00442 CPPAD_ASSERT_UNKNOWN( set_s[0].empty() ); 00443 size_t i; 00444 for(i = 0; i < m; i++) 00445 if( s[i] ) 00446 set_s[0].insert(i); 00447 00448 // compute sparsity pattern for H(x) = (S(x) * F)''(x) * R 00449 // (store it in v) 00450 f_.ForSparseJac(q, r); 00451 v = f_.RevSparseHes(q, set_s, transpose); 00452 00453 // compute sparsity pattern for V(x) = A(x) + H(x) 00454 for(i = 0; i < n; i++) 00455 { for(itr = a[i].begin(); itr != a[i].end(); itr++) 00456 { size_t j = *itr; 00457 CPPAD_ASSERT_UNKNOWN( j < q ); 00458 v[i].insert(j); 00459 } 00460 } 00461 00462 // no longer need the forward mode sparsity pattern 00463 // (have to reconstruct them every time) 00464 f_.size_forward_set(0); 00465 00466 return ok; 00467 } 00468 /*! 00469 Link from user_atomic to forward sparse Jacobian 00470 00471 \copydetails atomic_base::rev_sparse_hes 00472 */ 00473 virtual bool rev_sparse_hes( 00474 const vector<bool>& vx , 00475 const vector<bool>& s , 00476 vector<bool>& t , 00477 size_t q , 00478 const vector<bool>& r , 00479 const vector<bool>& u , 00480 vector<bool>& v ) 00481 { 00482 CPPAD_ASSERT_UNKNOWN( r.size() == v.size() ); 00483 CPPAD_ASSERT_UNKNOWN( s.size() == u.size() / q ); 00484 CPPAD_ASSERT_UNKNOWN( t.size() == v.size() / q ); 00485 size_t n = t.size(); 00486 bool ok = true; 00487 bool transpose = true; 00488 std::set<size_t>::const_iterator itr; 00489 size_t i, j; 00490 00491 // compute sparsity pattern for T(x) = S(x) * f'(x) 00492 t = f_.RevSparseJac(1, s); 00493 # ifndef NDEBUG 00494 for(j = 0; j < n; j++) 00495 CPPAD_ASSERT_UNKNOWN( vx[j] || ! t[j] ) 00496 # endif 00497 00498 // V(x) = f'(x)^T * g''(y) * f'(x) * R + g'(y) * f''(x) * R 00499 // U(x) = g''(y) * f'(x) * R 00500 // S(x) = g'(y) 00501 00502 // compute sparsity pattern for A(x) = f'(x)^T * U(x) 00503 vector<bool> a(n * q); 00504 a = f_.RevSparseJac(q, u, transpose); 00505 00506 // compute sparsity pattern for H(x) =(S(x) * F)''(x) * R 00507 // (store it in v) 00508 f_.ForSparseJac(q, r); 00509 v = f_.RevSparseHes(q, s, transpose); 00510 00511 // compute sparsity pattern for V(x) = A(x) + H(x) 00512 for(i = 0; i < n; i++) 00513 { for(j = 0; j < q; j++) 00514 v[ i * q + j ] |= a[ i * q + j]; 00515 } 00516 00517 // no longer need the forward mode sparsity pattern 00518 // (have to reconstruct them every time) 00519 f_.size_forward_set(0); 00520 00521 return ok; 00522 } 00523 }; 00524 00525 } // END_CPPAD_NAMESPACE 00526 # endif