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