CppAD: A C++ Algorithmic Differentiation Package
20130918
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_AD_FUN_INCLUDED 00003 # define CPPAD_AD_FUN_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 $begin ADFun$$ 00017 $spell 00018 xk 00019 Ind 00020 bool 00021 taylor_ 00022 sizeof 00023 const 00024 std 00025 ind_taddr_ 00026 dep_taddr_ 00027 $$ 00028 00029 $spell 00030 $$ 00031 00032 $section ADFun Objects$$ 00033 00034 $index ADFun, object$$ 00035 $index object, ADFun$$ 00036 00037 $head Purpose$$ 00038 An AD of $icode Base$$ 00039 $cref/operation sequence/glossary/Operation/Sequence/$$ 00040 is stored in an $code ADFun$$ object by its $cref FunConstruct$$. 00041 The $code ADFun$$ object can then be used to calculate function values, 00042 derivative values, and other values related to the corresponding function. 00043 00044 $childtable% 00045 cppad/local/independent.hpp% 00046 cppad/local/fun_construct.hpp% 00047 cppad/local/dependent.hpp% 00048 cppad/local/abort_recording.hpp% 00049 omh/seq_property.omh% 00050 cppad/local/fun_eval.hpp% 00051 cppad/local/drivers.hpp% 00052 cppad/local/fun_check.hpp% 00053 cppad/local/optimize.hpp% 00054 omh/check_for_nan.omh 00055 %$$ 00056 00057 $end 00058 */ 00059 00060 namespace CppAD { // BEGIN_CPPAD_NAMESPACE 00061 /*! 00062 \file ad_fun.hpp 00063 File used to define the ADFun<Base> class. 00064 */ 00065 00066 /*! 00067 Class used to hold function objects 00068 00069 \tparam Base 00070 A function object has a recording of <tt>AD<Base></tt> operations. 00071 It does it calculations using \c Base operations. 00072 */ 00073 00074 template <class Base> 00075 class ADFun { 00076 // ------------------------------------------------------------ 00077 // Private member variables 00078 private: 00079 /// Check for nan's and report message to user (default value is true). 00080 bool check_for_nan_; 00081 00082 /// debug checking number of comparision operations that changed 00083 size_t compare_change_; 00084 00085 /// number of orders stored in taylor_ 00086 size_t num_order_taylor_; 00087 00088 /// maximum number of orders that will fit in taylor_ 00089 size_t cap_order_taylor_; 00090 00091 /// number of variables in the recording (play_) 00092 size_t num_var_tape_; 00093 00094 /// tape address for the independent variables 00095 CppAD::vector<size_t> ind_taddr_; 00096 00097 /// tape address and parameter flag for the dependent variables 00098 CppAD::vector<size_t> dep_taddr_; 00099 00100 /// which dependent variables are actually parameters 00101 CppAD::vector<bool> dep_parameter_; 00102 00103 /// results of the forward mode calculations 00104 pod_vector<Base> taylor_; 00105 00106 /// which operations can be conditionally skipped 00107 /// Set during forward pass of order zero 00108 pod_vector<bool> cskip_op_; 00109 00110 /// Variable on the tape corresponding to each vecad load operation 00111 /// (if zero, the operation corresponds to a parameter). 00112 pod_vector<addr_t> load_op_; 00113 00114 /// the operation sequence corresponding to this object 00115 player<Base> play_; 00116 00117 /// Packed results of the forward mode Jacobian sparsity calculations. 00118 /// for_jac_sparse_pack_.n_set() != 0 implies other sparsity results 00119 /// are empty 00120 sparse_pack for_jac_sparse_pack_; 00121 00122 /// Set results of the forward mode Jacobian sparsity calculations 00123 /// for_jac_sparse_set_.n_set() != 0 implies for_sparse_pack_ is empty. 00124 CPPAD_INTERNAL_SPARSE_SET for_jac_sparse_set_; 00125 00126 // ------------------------------------------------------------ 00127 // Private member functions 00128 00129 /// change the operation sequence corresponding to this object 00130 template <typename ADvector> 00131 void Dependent(ADTape<Base> *tape, const ADvector &y); 00132 00133 // ------------------------------------------------------------ 00134 // vector of bool version of ForSparseJac 00135 // (see doxygen in for_sparse_jac.hpp) 00136 template <class VectorSet> 00137 void ForSparseJacCase( 00138 bool set_type , 00139 bool transpose , 00140 size_t q , 00141 const VectorSet& r , 00142 VectorSet& s 00143 ); 00144 // vector of std::set<size_t> version of ForSparseJac 00145 // (see doxygen in for_sparse_jac.hpp) 00146 template <class VectorSet> 00147 void ForSparseJacCase( 00148 const std::set<size_t>& set_type , 00149 bool transpose , 00150 size_t q , 00151 const VectorSet& r , 00152 VectorSet& s 00153 ); 00154 // ------------------------------------------------------------ 00155 // vector of bool version of RevSparseJac 00156 // (see doxygen in rev_sparse_jac.hpp) 00157 template <class VectorSet> 00158 void RevSparseJacCase( 00159 bool set_type , 00160 bool transpose , 00161 bool nz_compare, 00162 size_t p , 00163 const VectorSet& s , 00164 VectorSet& r 00165 ); 00166 // vector of std::set<size_t> version of RevSparseJac 00167 // (see doxygen in rev_sparse_jac.hpp) 00168 template <class VectorSet> 00169 void RevSparseJacCase( 00170 const std::set<size_t>& set_type , 00171 bool transpose , 00172 bool nz_compare, 00173 size_t p , 00174 const VectorSet& s , 00175 VectorSet& r 00176 ); 00177 // ------------------------------------------------------------ 00178 // vector of bool version of RevSparseHes 00179 // (see doxygen in rev_sparse_hes.hpp) 00180 template <class VectorSet> 00181 void RevSparseHesCase( 00182 bool set_type , 00183 bool transpose , 00184 size_t q , 00185 const VectorSet& s , 00186 VectorSet& h 00187 ); 00188 // vector of std::set<size_t> version of RevSparseHes 00189 // (see doxygen in rev_sparse_hes.hpp) 00190 template <class VectorSet> 00191 void RevSparseHesCase( 00192 const std::set<size_t>& set_type , 00193 bool transpose , 00194 size_t q , 00195 const VectorSet& s , 00196 VectorSet& h 00197 ); 00198 // ------------------------------------------------------------ 00199 // Forward mode version of SparseJacobian 00200 // (see doxygen in sparse_jacobian.hpp) 00201 template <class VectorBase, class VectorSet, class VectorSize> 00202 size_t SparseJacobianFor( 00203 const VectorBase& x , 00204 VectorSet& p_transpose , 00205 const VectorSize& row , 00206 const VectorSize& col , 00207 VectorBase& jac , 00208 sparse_jacobian_work& work 00209 ); 00210 // Reverse mode version of SparseJacobian 00211 // (see doxygen in sparse_jacobian.hpp) 00212 template <class VectorBase, class VectorSet, class VectorSize> 00213 size_t SparseJacobianRev( 00214 const VectorBase& x , 00215 VectorSet& p , 00216 const VectorSize& row , 00217 const VectorSize& col , 00218 VectorBase& jac , 00219 sparse_jacobian_work& work 00220 ); 00221 // ------------------------------------------------------------ 00222 // combined sparse_set, sparse_list and sparse_pack version of 00223 // SparseHessian (see doxygen in sparse_hessian.hpp) 00224 template <class VectorBase, class VectorSet, class VectorSize> 00225 size_t SparseHessianCompute( 00226 const VectorBase& x , 00227 const VectorBase& w , 00228 VectorSet& sparsity , 00229 const VectorSize& row , 00230 const VectorSize& col , 00231 VectorBase& hes , 00232 sparse_hessian_work& work 00233 ); 00234 // ------------------------------------------------------------ 00235 public: 00236 /// copy constructor 00237 ADFun(const ADFun& g) 00238 : num_var_tape_(0) 00239 { CppAD::ErrorHandler::Call( 00240 true, 00241 __LINE__, 00242 __FILE__, 00243 "ADFun(const ADFun& g)", 00244 "Attempting to use the ADFun<Base> copy constructor.\n" 00245 "Perhaps you are passing an ADFun<Base> object " 00246 "by value instead of by reference." 00247 ); 00248 } 00249 00250 /// default constructor 00251 ADFun(void); 00252 00253 // assignment operator 00254 // (see doxygen in fun_construct.hpp) 00255 void operator=(const ADFun& f); 00256 00257 /// sequence constructor 00258 template <typename ADvector> 00259 ADFun(const ADvector &x, const ADvector &y); 00260 00261 /// destructor 00262 ~ADFun(void) 00263 { } 00264 00265 /// set value of check_for_nan_ 00266 void check_for_nan(bool value) 00267 { check_for_nan_ = value; } 00268 bool check_for_nan(void) const 00269 { return check_for_nan_; } 00270 00271 /// assign a new operation sequence 00272 template <typename ADvector> 00273 void Dependent(const ADvector &x, const ADvector &y); 00274 00275 /// forward mode sweep 00276 template <typename VectorBase> 00277 VectorBase Forward(size_t q, 00278 const VectorBase& x, std::ostream& s = std::cout 00279 ); 00280 00281 /// reverse mode sweep 00282 template <typename VectorBase> 00283 VectorBase Reverse(size_t p, const VectorBase &v); 00284 00285 // forward mode Jacobian sparsity 00286 // (see doxygen documentation in for_sparse_jac.hpp) 00287 template <typename VectorSet> 00288 VectorSet ForSparseJac( 00289 size_t q, const VectorSet &r, bool transpose = false 00290 ); 00291 // reverse mode Jacobian sparsity 00292 // (see doxygen documentation in rev_sparse_jac.hpp) 00293 template <typename VectorSet> 00294 VectorSet RevSparseJac( 00295 size_t q, const VectorSet &s, bool transpose = false, 00296 bool nz_compare = false 00297 ); 00298 // reverse mode Hessian sparsity 00299 // (see doxygen documentation in rev_sparse_hes.hpp) 00300 template <typename VectorSet> 00301 VectorSet RevSparseHes( 00302 size_t q, const VectorSet &s, bool transpose = false 00303 ); 00304 00305 /// amount of memeory used for Jacobain sparsity pattern 00306 size_t size_forward_bool(void) const 00307 { return for_jac_sparse_pack_.memory(); } 00308 00309 /// free memeory used for Jacobain sparsity pattern 00310 void size_forward_bool(size_t zero) 00311 { CPPAD_ASSERT_KNOWN( 00312 zero == 0, 00313 "size_forward_bool: argument not equal to zero" 00314 ); 00315 for_jac_sparse_pack_.resize(0, 0); 00316 } 00317 00318 /// total number of elements used for Jacobian sparsity pattern 00319 size_t size_forward_set(void) const 00320 { return for_jac_sparse_set_.number_elements(); } 00321 00322 /// free memeory used for Jacobain sparsity pattern 00323 void size_forward_set(size_t zero) 00324 { CPPAD_ASSERT_KNOWN( 00325 zero == 0, 00326 "size_forward_bool: argument not equal to zero" 00327 ); 00328 for_jac_sparse_set_.resize(0, 0); 00329 } 00330 00331 /// number of operators in the operation sequence 00332 size_t size_op(void) const 00333 { return play_.num_op_rec(); } 00334 00335 /// number of operator arguments in the operation sequence 00336 size_t size_op_arg(void) const 00337 { return play_.num_op_arg_rec(); } 00338 00339 /// amount of memory required for the operation sequence 00340 size_t size_op_seq(void) const 00341 { return play_.Memory(); } 00342 00343 /// number of parameters in the operation sequence 00344 size_t size_par(void) const 00345 { return play_.num_par_rec(); } 00346 00347 /// number taylor coefficient orders calculated 00348 size_t size_order(void) const 00349 { return num_order_taylor_; } 00350 00351 /// number of characters in the operation sequence 00352 size_t size_text(void) const 00353 { return play_.num_text_rec(); } 00354 00355 /// number of variables in opertion sequence 00356 size_t size_var(void) const 00357 { return num_var_tape_; } 00358 00359 /// number of VecAD indices in the operation sequence 00360 size_t size_VecAD(void) const 00361 { return play_.num_vec_ind_rec(); } 00362 00363 /// set number of orders and directions currently allocated 00364 void capacity_order(size_t c); 00365 00366 /// number of variables in conditional expressions that can be skipped 00367 size_t number_skip(void); 00368 00369 /// number of independent variables 00370 size_t Domain(void) const 00371 { return ind_taddr_.size(); } 00372 00373 /// number of dependent variables 00374 size_t Range(void) const 00375 { return dep_taddr_.size(); } 00376 00377 /// is variable a parameter 00378 bool Parameter(size_t i) 00379 { CPPAD_ASSERT_KNOWN( 00380 i < dep_taddr_.size(), 00381 "Argument to Parameter is >= dimension of range space" 00382 ); 00383 return dep_parameter_[i]; 00384 } 00385 00386 # ifndef NDEBUG 00387 /// in not NDEBUG case, number of comparison operations that change 00388 size_t CompareChange(void) const 00389 { return compare_change_; } 00390 # endif 00391 00392 /// calculate entire Jacobian 00393 template <typename VectorBase> 00394 VectorBase Jacobian(const VectorBase &x); 00395 00396 /// calculate Hessian for one component of f 00397 template <typename VectorBase> 00398 VectorBase Hessian(const VectorBase &x, const VectorBase &w); 00399 template <typename VectorBase> 00400 VectorBase Hessian(const VectorBase &x, size_t i); 00401 00402 /// forward mode calculation of partial w.r.t one domain component 00403 template <typename VectorBase> 00404 VectorBase ForOne( 00405 const VectorBase &x , 00406 size_t j ); 00407 00408 /// reverse mode calculation of derivative of one range component 00409 template <typename VectorBase> 00410 VectorBase RevOne( 00411 const VectorBase &x , 00412 size_t i ); 00413 00414 /// forward mode calculation of a subset of second order partials 00415 template <typename VectorBase, typename VectorSize_t> 00416 VectorBase ForTwo( 00417 const VectorBase &x , 00418 const VectorSize_t &J , 00419 const VectorSize_t &K ); 00420 00421 /// reverse mode calculation of a subset of second order partials 00422 template <typename VectorBase, typename VectorSize_t> 00423 VectorBase RevTwo( 00424 const VectorBase &x , 00425 const VectorSize_t &I , 00426 const VectorSize_t &J ); 00427 00428 /// calculate sparse Jacobians 00429 template <typename VectorBase> 00430 VectorBase SparseJacobian( 00431 const VectorBase &x 00432 ); 00433 template <typename VectorBase, typename VectorSet> 00434 VectorBase SparseJacobian( 00435 const VectorBase &x , 00436 const VectorSet &p 00437 ); 00438 template <class VectorBase, class VectorSet, class VectorSize> 00439 size_t SparseJacobianForward( 00440 const VectorBase& x , 00441 const VectorSet& p , 00442 const VectorSize& r , 00443 const VectorSize& c , 00444 VectorBase& jac , 00445 sparse_jacobian_work& work 00446 ); 00447 template <class VectorBase, class VectorSet, class VectorSize> 00448 size_t SparseJacobianReverse( 00449 const VectorBase& x , 00450 const VectorSet& p , 00451 const VectorSize& r , 00452 const VectorSize& c , 00453 VectorBase& jac , 00454 sparse_jacobian_work& work 00455 ); 00456 00457 /// calculate sparse Hessians 00458 template <typename VectorBase> 00459 VectorBase SparseHessian( 00460 const VectorBase& x , 00461 const VectorBase& w 00462 ); 00463 template <typename VectorBase, typename VectorBool> 00464 VectorBase SparseHessian( 00465 const VectorBase& x , 00466 const VectorBase& w , 00467 const VectorBool& p 00468 ); 00469 template <class VectorBase, class VectorSet, class VectorSize> 00470 size_t SparseHessian( 00471 const VectorBase& x , 00472 const VectorBase& w , 00473 const VectorSet& p , 00474 const VectorSize& r , 00475 const VectorSize& c , 00476 VectorBase& hes , 00477 sparse_hessian_work& work 00478 ); 00479 00480 // Optimize the tape 00481 // (see doxygen documentation in optimize.hpp) 00482 void optimize(void); 00483 // ------------------- Deprecated ----------------------------- 00484 00485 /// deprecated: assign a new operation sequence 00486 template <typename ADvector> 00487 void Dependent(const ADvector &y); 00488 00489 /// Deprecated: number of variables in opertion sequence 00490 size_t Size(void) const 00491 { return num_var_tape_; } 00492 00493 /// Deprecated: # taylor_ coefficients currently stored 00494 /// (per variable,direction) 00495 size_t Order(void) const 00496 { return num_order_taylor_ - 1; } 00497 00498 /// Deprecated: amount of memory for this object 00499 /// Note that an approximation is used for the std::set<size_t> memory 00500 size_t Memory(void) const 00501 { size_t pervar = cap_order_taylor_ * sizeof(Base) 00502 + for_jac_sparse_pack_.memory() 00503 + 3 * sizeof(size_t) * for_jac_sparse_set_.number_elements(); 00504 size_t total = num_var_tape_ * pervar + play_.Memory(); 00505 return total; 00506 } 00507 00508 /// Deprecated: # taylor_ coefficient orderss stored 00509 /// (per variable,direction) 00510 size_t taylor_size(void) const 00511 { return num_order_taylor_; } 00512 00513 /// Deprecated: Does this AD operation sequence use 00514 /// VecAD<Base>::reference operands 00515 bool use_VecAD(void) const 00516 { return play_.num_vec_ind_rec() > 0; } 00517 00518 /// Deprecated: # taylor_ coefficient orders calculated 00519 /// (per variable,direction) 00520 size_t size_taylor(void) const 00521 { return num_order_taylor_; } 00522 00523 /// Deprecated: set number of orders currently allocated 00524 /// (per variable,direction) 00525 void capacity_taylor(size_t per_var); 00526 }; 00527 // --------------------------------------------------------------------------- 00528 00529 } // END_CPPAD_NAMESPACE 00530 00531 // non-user interfaces 00532 # include <cppad/local/forward0sweep.hpp> 00533 # include <cppad/local/forward1sweep.hpp> 00534 # include <cppad/local/reverse_sweep.hpp> 00535 # include <cppad/local/for_jac_sweep.hpp> 00536 # include <cppad/local/rev_jac_sweep.hpp> 00537 # include <cppad/local/rev_hes_sweep.hpp> 00538 00539 // user interfaces 00540 # include <cppad/local/parallel_ad.hpp> 00541 # include <cppad/local/independent.hpp> 00542 # include <cppad/local/dependent.hpp> 00543 # include <cppad/local/fun_construct.hpp> 00544 # include <cppad/local/abort_recording.hpp> 00545 # include <cppad/local/fun_eval.hpp> 00546 # include <cppad/local/drivers.hpp> 00547 # include <cppad/local/fun_check.hpp> 00548 # include <cppad/local/omp_max_thread.hpp> 00549 # include <cppad/local/optimize.hpp> 00550 00551 # endif