CppAD: A C++ Algorithmic Differentiation Package 20110419
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_AD_FUN_INCLUDED 00003 # define CPPAD_AD_FUN_INCLUDED 00004 00005 /* -------------------------------------------------------------------------- 00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-11 Bradley M. Bell 00007 00008 CppAD is distributed under multiple licenses. This distribution is under 00009 the terms of the 00010 Common 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 $italic Base$$ 00039 $xref/glossary/Operation/Sequence/operation sequence/1/$$ 00040 is stored in an $code ADFun$$ object by its $xref/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/omp_max_thread.hpp% 00054 cppad/local/optimize.hpp% 00055 omh/fun_deprecated.omh 00056 %$$ 00057 00058 $end 00059 */ 00060 00061 CPPAD_BEGIN_NAMESPACE 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 public: 00078 // ------------------------------------------------------------ 00079 // Private member variables 00080 private: 00081 /// debug checking number of comparision operations that changed 00082 size_t compare_change_; 00083 00084 /// number of taylor_ coefficieint per variable (currently stored) 00085 size_t taylor_per_var_; 00086 00087 /// number of columns currently allocated for taylor_ array 00088 size_t taylor_col_dim_; 00089 00090 /// number of rows (variables) in the recording (play_) 00091 size_t total_num_var_; 00092 00093 /// tape address for the independent variables 00094 CppAD::vector<size_t> ind_taddr_; 00095 00096 /// tape address and parameter flag for the dependent variables 00097 CppAD::vector<size_t> dep_taddr_; 00098 00099 /// which dependent variables are actually parameters 00100 CppAD::vector<bool> dep_parameter_; 00101 00102 /// the operation sequence corresponding to this object 00103 player<Base> play_; 00104 00105 /// results of the forward mode calculations 00106 Base *taylor_; 00107 00108 /// packed results of the forward mode Jacobian sparsity calculations 00109 /// (\c for_jac_sparse_pack_.n_set() != 0 implies 00110 /// for_jac_sparse_set_.n_set() == 0) 00111 sparse_pack for_jac_sparse_pack_; 00112 00113 /// set results of the forward mode Jacobian sparsity calculations 00114 /// (\c for_jac_sparse_set_.n_set() != 0 implies 00115 /// for_jac_sparse_pack_.n_set() == 0) 00116 sparse_set for_jac_sparse_set_; 00117 00118 // ------------------------------------------------------------ 00119 // Private member functions 00120 00121 /// change the operation sequence corresponding to this object 00122 template <typename ADvector> 00123 void Dependent(ADTape<Base> *tape, const ADvector &y); 00124 00125 // ------------------------------------------------------------ 00126 // vector of bool version of ForSparseJac 00127 // (see doxygen in for_sparse_jac.hpp) 00128 template <class VectorSet> 00129 void ForSparseJacCase( 00130 bool set_type , 00131 size_t q , 00132 const VectorSet& r , 00133 VectorSet& s 00134 ); 00135 // vector of std::set<size_t> version of ForSparseJac 00136 // (see doxygen in for_sparse_jac.hpp) 00137 template <class VectorSet> 00138 void ForSparseJacCase( 00139 const std::set<size_t>& set_type , 00140 size_t q , 00141 const VectorSet& r , 00142 VectorSet& s 00143 ); 00144 // ------------------------------------------------------------ 00145 // vector of bool version of RevSparseJac 00146 // (see doxygen in rev_sparse_jac.hpp) 00147 template <class VectorSet> 00148 void RevSparseJacCase( 00149 bool set_type , 00150 size_t p , 00151 const VectorSet& s , 00152 VectorSet& r 00153 ); 00154 // vector of std::set<size_t> version of RevSparseJac 00155 // (see doxygen in rev_sparse_jac.hpp) 00156 template <class VectorSet> 00157 void RevSparseJacCase( 00158 const std::set<size_t>& set_type , 00159 size_t p , 00160 const VectorSet& s , 00161 VectorSet& r 00162 ); 00163 // ------------------------------------------------------------ 00164 // vector of bool version of RevSparseHes 00165 // (see doxygen in rev_sparse_hes.hpp) 00166 template <class VectorSet> 00167 void RevSparseHesCase( 00168 bool set_type , 00169 size_t q , 00170 const VectorSet& s , 00171 VectorSet& h 00172 ); 00173 // vector of std::set<size_t> version of RevSparseHes 00174 // (see doxygen in rev_sparse_hes.hpp) 00175 template <class VectorSet> 00176 void RevSparseHesCase( 00177 const std::set<size_t>& set_type , 00178 size_t q , 00179 const VectorSet& s , 00180 VectorSet& h 00181 ); 00182 // ------------------------------------------------------------ 00183 // vector of bool version of SparseJacobian 00184 // (see doxygen in sparse_jacobian.hpp) 00185 template <class VectorBase, class VectorSet> 00186 void SparseJacobianCase( 00187 bool set_type , 00188 const VectorBase& x , 00189 const VectorSet& p , 00190 VectorBase& jac 00191 ); 00192 // vector of std::set<size_t> version of SparseJacobian 00193 // (see doxygen in sparse_jacobian.hpp) 00194 template <class VectorBase, class VectorSet> 00195 void SparseJacobianCase( 00196 const std::set<size_t>& set_type , 00197 const VectorBase& x , 00198 const VectorSet& p , 00199 VectorBase& jac 00200 ); 00201 // ------------------------------------------------------------ 00202 // vector of bool version of SparseHessian 00203 // (see doxygen in sparse_hessian.hpp) 00204 template <class VectorBase, class VectorSet> 00205 void SparseHessianCase( 00206 bool set_type , 00207 const VectorBase& x , 00208 const VectorBase& w , 00209 const VectorSet& p , 00210 VectorBase& hes 00211 ); 00212 // vector of std::set<size_t> version of SparseHessian 00213 // (see doxygen in sparse_hessian.hpp) 00214 template <class VectorBase, class VectorSet> 00215 void SparseHessianCase( 00216 const std::set<size_t>& set_type , 00217 const VectorBase& x , 00218 const VectorBase& w , 00219 const VectorSet& p , 00220 VectorBase& hes 00221 ); 00222 // ------------------------------------------------------------ 00223 public: 00224 /// copy constructor 00225 ADFun(const ADFun& g) 00226 : total_num_var_(0), taylor_(CPPAD_NULL) 00227 { CppAD::ErrorHandler::Call( 00228 true, 00229 __LINE__, 00230 __FILE__, 00231 "ADFun(const ADFun& g)", 00232 "Attempting to use the ADFun<Base> copy constructor.\n" 00233 "Perhaps you are passing an ADFun<Base> object " 00234 "by value instead of by reference." 00235 ); 00236 } 00237 00238 /// default constructor 00239 ADFun(void); 00240 00241 // assignment operator 00242 // (see doxygen in fun_construct.hpp) 00243 void operator=(const ADFun& f); 00244 00245 /// sequence constructor 00246 template <typename ADvector> 00247 ADFun(const ADvector &x, const ADvector &y); 00248 00249 /// destructor 00250 ~ADFun(void) 00251 { if( taylor_ != CPPAD_NULL ) 00252 CPPAD_TRACK_DEL_VEC(taylor_); 00253 } 00254 00255 /// deprecated: assign a new operation sequence 00256 template <typename ADvector> 00257 void Dependent(const ADvector &y); 00258 00259 /// assign a new operation sequence 00260 template <typename ADvector> 00261 void Dependent(const ADvector &x, const ADvector &y); 00262 00263 /// forward mode sweep 00264 template <typename VectorBase> 00265 VectorBase Forward(size_t p, const VectorBase &u); 00266 00267 /// reverse mode sweep 00268 template <typename VectorBase> 00269 VectorBase Reverse(size_t p, const VectorBase &v); 00270 00271 // forward mode Jacobian sparsity 00272 // (see doxygen documentation in for_sparse_jac.hpp) 00273 template <typename VectorSet> 00274 VectorSet ForSparseJac( 00275 size_t q, const VectorSet &r 00276 ); 00277 // reverse mode Jacobian sparsity 00278 // (see doxygen documentation in rev_sparse_jac.hpp) 00279 template <typename VectorSet> 00280 VectorSet RevSparseJac( 00281 size_t q, const VectorSet &s 00282 ); 00283 // reverse mode Hessian sparsity 00284 // (see doxygen documentation in rev_sparse_hes.hpp) 00285 template <typename VectorSet> 00286 VectorSet RevSparseHes( 00287 size_t q, const VectorSet &s 00288 ); 00289 00290 /// amount of memeory used for Jacobain sparsity pattern 00291 size_t size_forward_bool(void) const 00292 { return for_jac_sparse_pack_.memory(); } 00293 00294 /// free memeory used for Jacobain sparsity pattern 00295 void size_forward_bool(size_t zero) 00296 { CPPAD_ASSERT_KNOWN( 00297 zero == 0, 00298 "size_forward_bool: argument not equal to zero" 00299 ); 00300 for_jac_sparse_pack_.resize(0, 0); 00301 } 00302 00303 /// total number of elements used for Jacobian sparsity pattern 00304 size_t size_forward_set(void) const 00305 { return for_jac_sparse_set_.number_elements(); } 00306 00307 /// free memeory used for Jacobain sparsity pattern 00308 void size_forward_set(size_t zero) 00309 { CPPAD_ASSERT_KNOWN( 00310 zero == 0, 00311 "size_forward_bool: argument not equal to zero" 00312 ); 00313 for_jac_sparse_set_.resize(0, 0); 00314 } 00315 00316 /// number of operators in the operation sequence 00317 size_t size_op(void) const 00318 { return play_.num_rec_op(); } 00319 00320 /// number of operator arguments in the operation sequence 00321 size_t size_op_arg(void) const 00322 { return play_.num_rec_op_arg(); } 00323 00324 /// amount of memory required for the operation sequence 00325 size_t size_op_seq(void) const 00326 { return play_.Memory(); } 00327 00328 /// number of parameters in the operation sequence 00329 size_t size_par(void) const 00330 { return play_.num_rec_par(); } 00331 00332 /// number of taylor_ coefficients currently calculated (per variable) 00333 size_t size_taylor(void) const 00334 { return taylor_per_var_; } 00335 00336 /// number of characters in the operation sequence 00337 size_t size_text(void) const 00338 { return play_.num_rec_text(); } 00339 00340 /// number of variables in opertion sequence 00341 size_t size_var(void) const 00342 { return total_num_var_; } 00343 00344 /// number of VecAD indices in the operation sequence 00345 size_t size_VecAD(void) const 00346 { return play_.num_rec_vecad_ind(); } 00347 00348 /// set number of coefficients currently allocated (per variable) 00349 void capacity_taylor(size_t per_var); 00350 00351 /// number of independent variables 00352 size_t Domain(void) const 00353 { return ind_taddr_.size(); } 00354 00355 /// number of dependent variables 00356 size_t Range(void) const 00357 { return dep_taddr_.size(); } 00358 00359 /// is variable a parameter 00360 bool Parameter(size_t i) 00361 { CPPAD_ASSERT_KNOWN( 00362 i < dep_taddr_.size(), 00363 "Argument to Parameter is >= dimension of range space" 00364 ); 00365 return dep_parameter_[i]; 00366 } 00367 00368 # ifndef NDEBUG 00369 /// in not NDEBUG case, number of comparison operations that change 00370 size_t CompareChange(void) const 00371 { return compare_change_; } 00372 # endif 00373 00374 /// calculate entire Jacobian 00375 template <typename VectorBase> 00376 VectorBase Jacobian(const VectorBase &x); 00377 00378 /// calculate Hessian for one component of f 00379 template <typename VectorBase> 00380 VectorBase Hessian(const VectorBase &x, const VectorBase &w); 00381 template <typename VectorBase> 00382 VectorBase Hessian(const VectorBase &x, size_t i); 00383 00384 /// forward mode calculation of partial w.r.t one domain component 00385 template <typename VectorBase> 00386 VectorBase ForOne( 00387 const VectorBase &x , 00388 size_t j ); 00389 00390 /// reverse mode calculation of derivative of one range component 00391 template <typename VectorBase> 00392 VectorBase RevOne( 00393 const VectorBase &x , 00394 size_t i ); 00395 00396 /// forward mode calculation of a subset of second order partials 00397 template <typename VectorBase, typename VectorSize_t> 00398 VectorBase ForTwo( 00399 const VectorBase &x , 00400 const VectorSize_t &J , 00401 const VectorSize_t &K ); 00402 00403 /// reverse mode calculation of a subset of second order partials 00404 template <typename VectorBase, typename VectorSize_t> 00405 VectorBase RevTwo( 00406 const VectorBase &x , 00407 const VectorSize_t &I , 00408 const VectorSize_t &J ); 00409 00410 /// calculate sparse Jacobians 00411 template <typename VectorBase> 00412 VectorBase SparseJacobian(const VectorBase &x); 00413 00414 /// calculate sparse Jacobians 00415 template <typename VectorBase, typename VectorSet> 00416 VectorBase SparseJacobian(const VectorBase &x, const VectorSet &p); 00417 00418 /// calculate sparse Hessians 00419 template <typename VectorBase> 00420 VectorBase SparseHessian(const VectorBase &x, const VectorBase &w); 00421 00422 /// calculate sparse Hessians 00423 template <typename VectorBase, typename VectorBool> 00424 VectorBase SparseHessian( 00425 const VectorBase &x, const VectorBase &w, const VectorBool &p 00426 ); 00427 00428 // Optimize the tape 00429 // (see doxygen documentation in optimize.hpp) 00430 void optimize(void); 00431 // ------------------- Deprecated ----------------------------- 00432 00433 /// deprecated: number of variables in opertion sequence 00434 size_t Size(void) const 00435 { return total_num_var_; } 00436 00437 /// deprecated: # taylor_ coefficients currently stored (per variable) 00438 size_t Order(void) const 00439 { return taylor_per_var_ - 1; } 00440 00441 /// Deprecated: amount of memory for this object 00442 /// Note that an approximation is used for the std::set<size_t> memory 00443 size_t Memory(void) const 00444 { size_t pervar = taylor_col_dim_ * sizeof(Base) 00445 + for_jac_sparse_pack_.memory() 00446 + 3 * sizeof(size_t) * for_jac_sparse_set_.number_elements(); 00447 size_t total = total_num_var_ * pervar + play_.Memory(); 00448 return total; 00449 } 00450 00451 /// deprecated: # taylor_ coefficients stored (per variable) 00452 size_t taylor_size(void) const 00453 { return taylor_per_var_; } 00454 00455 /// deprecated: Does this AD operation sequence use 00456 //VecAD<Base>::reference operands 00457 bool use_VecAD(void) const 00458 { return play_.num_rec_vecad_ind() > 0; } 00459 }; 00460 // --------------------------------------------------------------------------- 00461 00462 CPPAD_END_NAMESPACE 00463 00464 // non-user interfaces 00465 # include <cppad/local/forward0sweep.hpp> 00466 # include <cppad/local/forward_sweep.hpp> 00467 # include <cppad/local/reverse_sweep.hpp> 00468 # include <cppad/local/for_jac_sweep.hpp> 00469 # include <cppad/local/rev_jac_sweep.hpp> 00470 # include <cppad/local/rev_hes_sweep.hpp> 00471 00472 // user interfaces 00473 # include <cppad/local/independent.hpp> 00474 # include <cppad/local/dependent.hpp> 00475 # include <cppad/local/fun_construct.hpp> 00476 # include <cppad/local/abort_recording.hpp> 00477 # include <cppad/local/fun_eval.hpp> 00478 # include <cppad/local/drivers.hpp> 00479 # include <cppad/local/fun_check.hpp> 00480 # include <cppad/local/omp_max_thread.hpp> 00481 # include <cppad/local/optimize.hpp> 00482 00483 # endif