CppAD: A C++ Algorithmic Differentiation Package  20130918
ad_fun.hpp
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines