source: branches/cache/cppad/local/ad_fun.hpp @ 3326

Last change on this file since 3326 was 3326, checked in by bradbell, 6 years ago

Implement step 4.1 of plan, computation of last_use.

cache.sh: test of computation of last_use.
optimize.hpp: fix comment.

  • Property svn:keywords set to Id
File size: 16.7 KB
Line 
1/* $Id: ad_fun.hpp 3326 2014-09-13 12:24:13Z bradbell $ */
2# ifndef CPPAD_AD_FUN_INCLUDED
3# define CPPAD_AD_FUN_INCLUDED
4
5/* --------------------------------------------------------------------------
6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-14 Bradley M. Bell
7
8CppAD is distributed under multiple licenses. This distribution is under
9the terms of the
10                    Eclipse Public License Version 1.0.
11
12A copy of this license is included in the COPYING file of this distribution.
13Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
14-------------------------------------------------------------------------- */
15/*
16$begin ADFun$$
17$spell
18        xk
19        Ind
20        bool
21        taylor_
22        sizeof
23        const
24        std
25        ind_taddr_
26        dep_taddr_
27$$
28
29$spell
30$$
31
32$section ADFun Objects$$
33
34$index ADFun, object$$
35$index object, ADFun$$
36
37$head Purpose$$
38An AD of $icode Base$$
39$cref/operation sequence/glossary/Operation/Sequence/$$
40is stored in an $code ADFun$$ object by its $cref FunConstruct$$.
41The $code ADFun$$ object can then be used to calculate function values,
42derivative values, and other values related to the corresponding function.
43
44$childtable%
45        cppad/local/independent.hpp%
46        cppad/local/fun_construct.hpp%
47        cppad/local/dependent.hpp%
48        cppad/local/abort_recording.hpp%
49        omh/seq_property.omh%
50        cppad/local/fun_eval.hpp%
51        cppad/local/drivers.hpp%
52        cppad/local/fun_check.hpp%
53        cppad/local/optimize.hpp%
54        cppad/local/cache.hpp%
55        omh/check_for_nan.omh
56%$$
57
58$end
59*/
60
61namespace CppAD { // BEGIN_CPPAD_NAMESPACE
62/*!
63\file ad_fun.hpp
64File used to define the ADFun<Base> class.
65*/
66
67/*!
68Class used to hold function objects
69
70\tparam Base
71A function object has a recording of <tt>AD<Base></tt> operations.
72It does it calculations using \c Base operations.
73*/
74
75template <class Base>
76class ADFun {
77// ------------------------------------------------------------
78// Private member variables
79private:
80        /// Check for nan's and report message to user (default value is true).
81        bool check_for_nan_;
82
83        /// debug checking number of comparision operations that changed
84        size_t compare_change_;
85
86        /// number of orders stored in taylor_
87        size_t num_order_taylor_;
88
89        /// maximum number of orders that will fit in taylor_
90        size_t cap_order_taylor_;
91
92        /// number of directions stored in taylor_
93        size_t num_direction_taylor_;
94
95        /// number of variables in the recording (play_)
96        size_t num_var_tape_;
97
98        /// tape address for the independent variables
99        CppAD::vector<size_t> ind_taddr_;
100
101        /// tape address and parameter flag for the dependent variables
102        CppAD::vector<size_t> dep_taddr_;
103
104        /// which dependent variables are actually parameters
105        CppAD::vector<bool>   dep_parameter_;
106
107        /// results of the forward mode calculations
108        pod_vector<Base> taylor_;
109
110        /// which operations can be conditionally skipped
111        /// Set during forward pass of order zero
112        pod_vector<bool> cskip_op_;
113
114        /// Variable on the tape corresponding to each vecad load operation
115        /// (if zero, the operation corresponds to a parameter).
116        pod_vector<addr_t> load_op_;
117
118        /// the operation sequence corresponding to this object
119        player<Base> play_;
120
121        /// Packed results of the forward mode Jacobian sparsity calculations.
122        /// for_jac_sparse_pack_.n_set() != 0  implies other sparsity results
123        /// are empty
124        sparse_pack      for_jac_sparse_pack_;
125
126        /// Set results of the forward mode Jacobian sparsity calculations
127        /// for_jac_sparse_set_.n_set() != 0  implies for_sparse_pack_ is empty.
128        CPPAD_INTERNAL_SPARSE_SET  for_jac_sparse_set_;
129
130// ------------------------------------------------------------
131// Private member functions
132
133        /// change the operation sequence corresponding to this object
134        template <typename ADvector>
135        void Dependent(ADTape<Base> *tape, const ADvector &y);
136
137        // ------------------------------------------------------------
138        // vector of bool version of ForSparseJac
139        // (see doxygen in for_sparse_jac.hpp)
140        template <class VectorSet>
141        void ForSparseJacCase(
142                bool               set_type  ,
143                bool               transpose ,
144                size_t             q         ,
145                const VectorSet&   r         , 
146                VectorSet&         s
147        );
148        // vector of std::set<size_t> version of ForSparseJac
149        // (see doxygen in for_sparse_jac.hpp)
150        template <class VectorSet>
151        void ForSparseJacCase(
152                const std::set<size_t>&  set_type  ,
153                bool                     transpose ,
154                size_t                   q         ,
155                const VectorSet&         r         , 
156                VectorSet&               s
157        );
158        // ------------------------------------------------------------
159        // vector of bool version of RevSparseJac
160        // (see doxygen in rev_sparse_jac.hpp)
161        template <class VectorSet>
162        void RevSparseJacCase(
163                bool               set_type  ,
164                bool               transpose ,
165                bool               nz_compare,
166                size_t             p         ,
167                const VectorSet&   s         , 
168                VectorSet&         r
169        );
170        // vector of std::set<size_t> version of RevSparseJac
171        // (see doxygen in rev_sparse_jac.hpp)
172        template <class VectorSet>
173        void RevSparseJacCase(
174                const std::set<size_t>&  set_type  ,
175                bool                     transpose ,
176                bool                     nz_compare,
177                size_t                   p         ,
178                const VectorSet&         s         , 
179                VectorSet&               r
180        );
181        // ------------------------------------------------------------
182        // vector of bool version of RevSparseHes
183        // (see doxygen in rev_sparse_hes.hpp)
184        template <class VectorSet>
185        void RevSparseHesCase(
186                bool               set_type  ,
187                bool               transpose ,
188                size_t             q         ,
189                const VectorSet&   s         , 
190                VectorSet&         h
191        );
192        // vector of std::set<size_t> version of RevSparseHes
193        // (see doxygen in rev_sparse_hes.hpp)
194        template <class VectorSet>
195        void RevSparseHesCase(
196                const std::set<size_t>&  set_type  ,
197                bool                     transpose ,
198                size_t                   q         ,
199                const VectorSet&         s         , 
200                VectorSet&               h
201        );
202        // ------------------------------------------------------------
203        // Forward mode version of SparseJacobian
204        // (see doxygen in sparse_jacobian.hpp)
205        template <class VectorBase, class VectorSet, class VectorSize>
206        size_t SparseJacobianFor(
207                const VectorBase&           x               ,
208                      VectorSet&            p_transpose     ,
209                const VectorSize&           row             ,
210                const VectorSize&           col             ,
211                      VectorBase&           jac             ,
212                      sparse_jacobian_work& work
213        );
214        // Reverse mode version of SparseJacobian
215        // (see doxygen in sparse_jacobian.hpp)
216        template <class VectorBase, class VectorSet, class VectorSize>
217        size_t SparseJacobianRev(
218                const VectorBase&           x               ,
219                      VectorSet&            p               ,
220                const VectorSize&           row             ,
221                const VectorSize&           col             ,
222                      VectorBase&           jac             ,
223                      sparse_jacobian_work& work
224        );
225        // ------------------------------------------------------------
226        // combined sparse_set, sparse_list and sparse_pack version of
227        // SparseHessian (see doxygen in sparse_hessian.hpp)
228        template <class VectorBase, class VectorSet, class VectorSize>
229        size_t SparseHessianCompute(
230                const VectorBase&              x           ,
231                const VectorBase&              w           ,
232                      VectorSet&               sparsity    ,
233                const VectorSize&              row         ,
234                const VectorSize&              col         ,
235                      VectorBase&              hes         ,
236                      sparse_hessian_work&     work
237        );
238// ------------------------------------------------------------
239public:
240        /// copy constructor
241        ADFun(const ADFun& g) 
242        : num_var_tape_(0)
243        {       CppAD::ErrorHandler::Call(
244                true,
245                __LINE__,
246                __FILE__,
247                "ADFun(const ADFun& g)",
248                "Attempting to use the ADFun<Base> copy constructor.\n"
249                "Perhaps you are passing an ADFun<Base> object "
250                "by value instead of by reference."
251                );
252         }
253
254        /// default constructor
255        ADFun(void); 
256
257        // assignment operator
258        // (see doxygen in fun_construct.hpp)
259        void operator=(const ADFun& f);
260
261        /// sequence constructor
262        template <typename ADvector>
263        ADFun(const ADvector &x, const ADvector &y);
264
265        /// destructor
266        ~ADFun(void)
267        { }
268
269        /// set value of check_for_nan_
270        void check_for_nan(bool value)
271        {       check_for_nan_ = value; }
272        bool check_for_nan(void) const
273        {       return check_for_nan_; }
274
275        /// assign a new operation sequence
276        template <typename ADvector>
277        void Dependent(const ADvector &x, const ADvector &y);
278
279        /// forward mode user API, one order multiple directions.
280        template <typename VectorBase>
281        VectorBase Forward(size_t q, size_t r, const VectorBase& x);
282
283        /// forward mode user API, multiple directions one order.
284        template <typename VectorBase>
285        VectorBase Forward(size_t q, 
286                const VectorBase& x, std::ostream& s = std::cout
287        );
288
289        /// reverse mode sweep
290        template <typename VectorBase>
291        VectorBase Reverse(size_t p, const VectorBase &v);
292
293        // forward mode Jacobian sparsity
294        // (see doxygen documentation in for_sparse_jac.hpp)
295        template <typename VectorSet>
296        VectorSet ForSparseJac(
297                size_t q, const VectorSet &r, bool transpose = false
298        );
299        // reverse mode Jacobian sparsity
300        // (see doxygen documentation in rev_sparse_jac.hpp)
301        template <typename VectorSet>
302        VectorSet RevSparseJac(
303                size_t q, const VectorSet &s, bool transpose = false,
304                bool nz_compare = false
305        );
306        // reverse mode Hessian sparsity
307        // (see doxygen documentation in rev_sparse_hes.hpp)
308        template <typename VectorSet>
309        VectorSet RevSparseHes(
310                size_t q, const VectorSet &s, bool transpose = false
311        );
312
313        /// amount of memeory used for Jacobain sparsity pattern
314        size_t size_forward_bool(void) const
315        {       return for_jac_sparse_pack_.memory(); }
316
317        /// free memeory used for Jacobain sparsity pattern
318        void size_forward_bool(size_t zero) 
319        {       CPPAD_ASSERT_KNOWN(
320                        zero == 0,
321                        "size_forward_bool: argument not equal to zero"
322                );
323                for_jac_sparse_pack_.resize(0, 0); 
324        }
325
326        /// total number of elements used for Jacobian sparsity pattern
327        size_t size_forward_set(void) const
328        {       return for_jac_sparse_set_.number_elements(); }
329
330        /// free memeory used for Jacobain sparsity pattern
331        void size_forward_set(size_t zero)
332        {       CPPAD_ASSERT_KNOWN(
333                        zero == 0,
334                        "size_forward_bool: argument not equal to zero"
335                );
336                for_jac_sparse_set_.resize(0, 0); 
337        }
338
339        /// number of operators in the operation sequence
340        size_t size_op(void) const
341        {       return play_.num_op_rec(); }
342
343        /// number of operator arguments in the operation sequence
344        size_t size_op_arg(void) const
345        {       return play_.num_op_arg_rec(); }
346
347        /// amount of memory required for the operation sequence
348        size_t size_op_seq(void) const
349        {       return play_.Memory(); }
350
351        /// number of parameters in the operation sequence
352        size_t size_par(void) const
353        {       return play_.num_par_rec(); }
354
355        /// number taylor coefficient orders calculated
356        size_t size_order(void) const
357        {       return num_order_taylor_; } 
358
359        /// number taylor coefficient directions calculated
360        size_t size_direction(void) const
361        {       return num_direction_taylor_; } 
362
363        /// number of characters in the operation sequence
364        size_t size_text(void) const
365        {       return play_.num_text_rec(); }
366
367        /// number of variables in opertion sequence
368        size_t size_var(void) const
369        {       return num_var_tape_; }
370
371        /// number of VecAD indices in the operation sequence
372        size_t size_VecAD(void) const
373        {       return play_.num_vec_ind_rec(); }
374
375        /// set number of orders currently allocated (user API)
376        void capacity_order(size_t c);
377
378        /// set number of orders and directions currently allocated
379        void capacity_order(size_t c, size_t r);   
380
381        /// number of variables in conditional expressions that can be skipped
382        size_t number_skip(void);   
383
384        /// number of independent variables
385        size_t Domain(void) const
386        {       return ind_taddr_.size(); }
387
388        /// number of dependent variables
389        size_t Range(void) const
390        {       return dep_taddr_.size(); }
391
392        /// is variable a parameter
393        bool Parameter(size_t i)
394        {       CPPAD_ASSERT_KNOWN(
395                        i < dep_taddr_.size(),
396                        "Argument to Parameter is >= dimension of range space"
397                );
398                return dep_parameter_[i]; 
399        }
400
401# ifndef NDEBUG
402        /// in not NDEBUG case, number of comparison operations that change
403        size_t CompareChange(void) const
404        {       return compare_change_; }
405# endif
406
407        /// calculate entire Jacobian
408        template <typename VectorBase>
409        VectorBase Jacobian(const VectorBase &x); 
410
411        /// calculate Hessian for one component of f
412        template <typename VectorBase>
413        VectorBase Hessian(const VectorBase &x, const VectorBase &w); 
414        template <typename VectorBase>
415        VectorBase Hessian(const VectorBase &x, size_t i); 
416
417        /// forward mode calculation of partial w.r.t one domain component
418        template <typename VectorBase>
419        VectorBase ForOne(
420                const VectorBase   &x ,
421                size_t              j );
422
423        /// reverse mode calculation of derivative of one range component
424        template <typename VectorBase>
425        VectorBase RevOne(
426                const VectorBase   &x ,
427                size_t              i );
428
429        /// forward mode calculation of a subset of second order partials
430        template <typename VectorBase, typename VectorSize_t>
431        VectorBase ForTwo(
432                const VectorBase   &x ,
433                const VectorSize_t &J ,
434                const VectorSize_t &K );
435
436        /// reverse mode calculation of a subset of second order partials
437        template <typename VectorBase, typename VectorSize_t>
438        VectorBase RevTwo(
439                const VectorBase   &x ,
440                const VectorSize_t &I ,
441                const VectorSize_t &J );
442
443        /// calculate sparse Jacobians
444        template <typename VectorBase>
445        VectorBase SparseJacobian(
446                const VectorBase &x
447        ); 
448        template <typename VectorBase, typename VectorSet>
449        VectorBase SparseJacobian(
450                const VectorBase &x , 
451                const VectorSet  &p
452        ); 
453        template <class VectorBase, class VectorSet, class VectorSize>
454        size_t SparseJacobianForward(
455                const VectorBase&     x     ,
456                const VectorSet&      p     ,
457                const VectorSize&     r     ,
458                const VectorSize&     c     ,
459                VectorBase&           jac   ,
460                sparse_jacobian_work& work
461        );
462        template <class VectorBase, class VectorSet, class VectorSize>
463        size_t SparseJacobianReverse(
464                const VectorBase&     x    ,
465                const VectorSet&      p    ,
466                const VectorSize&     r    ,
467                const VectorSize&     c    ,
468                VectorBase&           jac  ,
469                sparse_jacobian_work& work
470        );
471
472        /// calculate sparse Hessians
473        template <typename VectorBase>
474        VectorBase SparseHessian(
475                const VectorBase&    x  , 
476                const VectorBase&    w
477        ); 
478        template <typename VectorBase, typename VectorBool>
479        VectorBase SparseHessian(
480                const VectorBase&    x  ,
481                const VectorBase&    w  ,
482                const VectorBool&    p
483        ); 
484        template <class VectorBase, class VectorSet, class VectorSize>
485        size_t SparseHessian(
486                const VectorBase&    x   ,
487                const VectorBase&    w   ,
488                const VectorSet&     p   ,
489                const VectorSize&    r   ,
490                const VectorSize&    c   ,
491                VectorBase&          hes ,
492                sparse_hessian_work& work
493        ); 
494
495        // Optimize the tape
496        // (see doxygen documentation in optimize.hpp)
497        void optimize(void);
498
499        // Create a cache for tape values
500        void cache(void);
501        // ------------------- Deprecated -----------------------------
502
503        /// deprecated: assign a new operation sequence
504        template <typename ADvector>
505        void Dependent(const ADvector &y);
506
507        /// Deprecated: number of variables in opertion sequence
508        size_t Size(void) const
509        {       return num_var_tape_; }
510
511        /// Deprecated: # taylor_ coefficients currently stored
512        /// (per variable,direction)
513        size_t Order(void) const
514        {       return num_order_taylor_ - 1; }
515
516        /// Deprecated: amount of memory for this object
517        /// Note that an approximation is used for the std::set<size_t> memory
518        size_t Memory(void) const
519        {       size_t pervar  = cap_order_taylor_ * sizeof(Base)
520                + for_jac_sparse_pack_.memory()
521                + 3 * sizeof(size_t) * for_jac_sparse_set_.number_elements();
522                size_t total   = num_var_tape_  * pervar + play_.Memory();
523                return total;
524        }
525
526        /// Deprecated: # taylor_ coefficient orderss stored
527        /// (per variable,direction)
528        size_t taylor_size(void) const
529        {       return num_order_taylor_; } 
530
531        /// Deprecated: Does this AD operation sequence use
532        /// VecAD<Base>::reference operands
533        bool use_VecAD(void) const
534        {       return play_.num_vec_ind_rec() > 0; }
535
536        /// Deprecated: # taylor_ coefficient orders calculated
537        /// (per variable,direction)
538        size_t size_taylor(void) const
539        {       return num_order_taylor_; } 
540
541        /// Deprecated: set number of orders currently allocated
542        /// (per variable,direction)
543        void capacity_taylor(size_t per_var);   
544};
545// ---------------------------------------------------------------------------
546
547} // END_CPPAD_NAMESPACE
548
549// non-user interfaces
550# include <cppad/local/forward0sweep.hpp>
551# include <cppad/local/forward1sweep.hpp>
552# include <cppad/local/forward2sweep.hpp>
553# include <cppad/local/reverse_sweep.hpp>
554# include <cppad/local/for_jac_sweep.hpp>
555# include <cppad/local/rev_jac_sweep.hpp>
556# include <cppad/local/rev_hes_sweep.hpp>
557
558// user interfaces
559# include <cppad/local/parallel_ad.hpp>
560# include <cppad/local/independent.hpp>
561# include <cppad/local/dependent.hpp>
562# include <cppad/local/fun_construct.hpp>
563# include <cppad/local/abort_recording.hpp>
564# include <cppad/local/fun_eval.hpp>
565# include <cppad/local/drivers.hpp>
566# include <cppad/local/fun_check.hpp>
567# include <cppad/local/omp_max_thread.hpp>
568# include <cppad/local/optimize.hpp>
569# include <cppad/local/cache.hpp>
570
571# endif
Note: See TracBrowser for help on using the repository browser.