source: trunk/cppad/local/atomic_base.hpp @ 2901

Last change on this file since 2901 was 2901, checked in by bradbell, 7 years ago

Detect and report atomic function name when user
attempts to use an atomic function that has been deleted.

  • Property svn:keywords set to Id
File size: 36.3 KB
Line 
1/* $Id: atomic_base.hpp 2901 2013-09-19 12:07:15Z bradbell $ */
2# ifndef CPPAD_ATOMIC_BASE_INCLUDED
3# define CPPAD_ATOMIC_BASE_INCLUDED
4
5/* --------------------------------------------------------------------------
6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 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# include <set>
17# include <cppad/local/cppad_assert.hpp>
18// needed before one can use CPPAD_ASSERT_FIRST_CALL_NOT_PARALLEL
19# include <cppad/thread_alloc.hpp>
20
21CPPAD_BEGIN_NAMESPACE
22/*!
23\defgroup atomic_base.hpp atomic_base.hpp
24\{
25\file atomic_base.hpp
26Base class for atomic user operations.
27*/
28
29template <class Base>
30class atomic_base {
31// ===================================================================
32public:
33        enum option_enum { bool_sparsity_enum, set_sparsity_enum};
34private:
35        // ------------------------------------------------------
36        // constants
37        //
38        /// index of this object in class_object
39        const size_t index_;
40
41        // -----------------------------------------------------
42        // variables
43        //
44        /// sparsity pattern this object is currently using
45        /// (set by constructor and option member functions)
46        option_enum sparsity_;
47
48        /// temporary work space used afun, declared here to avoid memory
49        /// allocation/deallocation for each call to afun
50        vector<bool>  afun_vx_[CPPAD_MAX_NUM_THREADS];
51        vector<bool>  afun_vy_[CPPAD_MAX_NUM_THREADS];
52        vector<Base>  afun_tx_[CPPAD_MAX_NUM_THREADS];
53        vector<Base>  afun_ty_[CPPAD_MAX_NUM_THREADS];
54
55        // -----------------------------------------------------
56        // static member functions
57        //
58        /// List of all the object in this class
59        static std::vector<atomic_base *>& class_object(void)
60        {       CPPAD_ASSERT_FIRST_CALL_NOT_PARALLEL;
61                static std::vector<atomic_base *> list_;
62                return list_;
63        }
64        /// List of names for each object in this class
65        static std::vector<std::string>& class_name(void)
66        {       CPPAD_ASSERT_FIRST_CALL_NOT_PARALLEL;
67                static std::vector<std::string> list_;
68                return list_;
69        }
70        // =====================================================================
71public:
72        // -----------------------------------------------------
73        // member functions not in user API
74        //
75        /// current sparsity setting
76        option_enum sparsity(void) const
77        {       return sparsity_; }
78
79        /// Name corresponding to a base_atomic object
80        const std::string& afun_name(void) const
81        {       return class_name()[index_]; }
82/*
83$begin atomic_ctor$$
84$spell
85        std
86        afun
87        arg
88        CppAD
89        bool
90        ctor
91        const
92$$
93
94$section Atomic Function Constructor$$
95$index constructor, atomic function$$
96$index atomic, function constructor$$
97$index function, atomic constructor$$
98
99$head Syntax$$
100$icode%atomic_user afun%(%ctor_arg_list%)
101%$$
102$codei%atomic_base<%Base%>(%name%)
103%$$
104
105$head atomic_user$$
106
107$subhead ctor_arg_list$$
108Is a list of arguments for the $icode atomic_user$$ constructor.
109
110$subhead afun$$
111The object $icode afun$$ must stay in scope for as long
112as the corresponding atomic function is used.
113This includes use by any $cref/ADFun<Base>/ADFun/$$ that
114has this $icode atomic_user$$ operation in its
115$cref/operation sequence/glossary/Operation/Sequence/$$.
116
117$subhead Implementation$$
118The user defined $icode atomic_user$$ class is a publicly derived class of
119$codei%atomic_base<%Base%>%$$.
120It should be declared as follows:
121$codei%
122        class %atomic_user% : public CppAD::atomic_base<%Base%> {
123        public:
124                %atomic_user%(%ctor_arg_list%) : atomic_base<%Base%>(%name%)
125        %...%
126        };
127%$$
128where $icode ...$$ 
129denotes the rest of the implementation of the derived class.
130This includes completing the constructor and
131all the virtual functions that have their
132$code atomic_base$$ implementations replaced by
133$icode atomic_user$$ implementations.
134
135$head atomic_base$$
136
137$subhead Restrictions$$
138The $code atomic_base$$ constructor cannot be called in
139$cref/parallel/ta_in_parallel/$$ mode.
140
141$subhead Base$$
142The template parameter determines the
143$icode Base$$ type for this $codei%AD<%Base%>%$$ atomic operation.
144
145$subhead name$$
146This $icode atomic_base$$ constructor argument has either of the
147following prototypes
148$codei%
149        const char*        %name%
150        const std::string& %name%
151%$$
152It is the name for this atomic function and is used for error reporting.
153The suggested value for $icode name$$ is $icode afun$$ or $icode atomic_user$$,
154i.e., the name of the corresponding atomic object or class.
155
156$end
157*/
158/*!
159Base class for atomic_user functions.
160
161\tparam Base
162This class is used for defining an AD<Base> atomic operation y = f(x).
163*/
164/// make sure user does not invoke the default constructor
165atomic_base(void)
166{       CPPAD_ASSERT_KNOWN(false,
167                "Attempt to use the atomic_base default constructor"
168        );
169}
170/*!
171Constructor
172
173\param name
174name used for error reporting
175*/
176atomic_base( const std::string&  name) :
177index_( class_object().size() )     ,
178sparsity_( set_sparsity_enum )
179{       CPPAD_ASSERT_KNOWN(
180                ! thread_alloc::in_parallel() ,
181                "atomic_base: constructor cannot be called in parallel mode."
182        );
183        class_object().push_back(this);
184        class_name().push_back(name);
185        CPPAD_ASSERT_UNKNOWN( class_object().size() == class_name().size() );
186}
187/// destructor informs CppAD that this atomic function with this index
188/// has dropped out of scope by setting its pointer to null
189virtual ~atomic_base(void)
190{       CPPAD_ASSERT_UNKNOWN( class_object().size() > index_ );
191        // change object pointer to null, but leave name for error reporting
192        class_object()[index_] = CPPAD_NULL;
193}
194/// atomic_base function object corresponding to a certain index
195static atomic_base* class_object(size_t index)
196{       CPPAD_ASSERT_UNKNOWN( class_object().size() > index );
197        return class_object()[index];
198}
199/// atomic_base function name corresponding to a certain index
200static const std::string& class_name(size_t index)
201{       CPPAD_ASSERT_UNKNOWN( class_name().size() > index );
202        return class_name()[index];
203}
204/*
205$begin atomic_option$$
206$spell
207        enum
208        afun
209        bool
210        CppAD
211        std
212        typedef
213$$
214
215$section Set Atomic Function Options$$
216$index atomic, options$$
217$index options, atomic$$
218
219$head Syntax$$
220$icode%afun%.option(%option_value%)%$$
221
222$head atomic_sparsity$$
223$index atomic_sparsity$$
224$index sparsity, atomic$$
225You can used this option to set to type used for
226$icode afun$$ sparsity patterns.
227This does not apply individual calls to $icode afun$$,
228but rather all its uses between when the sparsity pattern is set and when
229it is changed.
230If neither the $code set_sparsity_enum$$ or
231$code bool_sparsity_enum$$ option is set,
232the type for $icode atomic_sparsity$$ is one of the two choices below
233(and otherwise unspecified).
234
235$subhead bool_sparsity_enum$$
236$index bool_sparsity_enum$$
237If $icode option_value$$ is $code atomic_base::bool_sparsity_enum$$,
238then the type used by $icode afun$$ for
239$cref/sparsity patterns/glossary/Sparsity Pattern/$$,
240(after the option is set) will be
241$codei%
242        typedef CppAD::vector<bool> %atomic_sparsity%
243%$$
244If $icode r$$ is a sparsity pattern
245for a matrix $latex R \in B^{p \times q}$$:
246$icode%r%.size() == %p% * %q%$$.
247
248$subhead set_sparsity_enum$$
249$index set_sparsity_enum$$
250If $icode option_value$$ is $code atomic_base::set_sparsity_enum$$,
251then the type used by $icode afun$$ for
252$cref/sparsity patterns/glossary/Sparsity Pattern/$$,
253(after the option is set) will be
254$codei%
255        typedef CppAD::vector< std::set<size_t> > %atomic_sparsity%
256%$$
257If $icode r$$ is a sparsity pattern
258for a matrix $code R \in B^{p \times q}$$:
259$icode%r%.size() == %p%$$, and for $latex i = 0 , \ldots , p-1$$,
260the elements of $icode%r%[%i%]%$$ are between zero and $latex q-1$$ inclusive.
261
262$end
263*/
264void option(enum option_enum option_value)
265{       switch( option_value )
266        {       case bool_sparsity_enum:
267                case set_sparsity_enum:
268                sparsity_ = option_value;
269                break;
270
271                default:
272                CPPAD_ASSERT_KNOWN(
273                        false,
274                        "atoic_base::option: option_value is not valid"
275                );
276        }
277        return;
278}
279/*
280-----------------------------------------------------------------------------
281$begin atomic_afun$$
282
283$spell
284        afun
285        const
286        CppAD
287$$
288
289$section Using an Atomic Function$$
290$index atomic, use function$$
291
292$head Syntax$$
293$icode%afun%(%ax%, %ay%)%$$
294
295$head Purpose$$
296Given $icode ax$$,
297this call computes the corresponding value of $icode ay$$.
298If $codei%AD<%Base%>%$$ operations are being recorded,
299it enters the computation as an atomic operation in the recording;
300see $cref/start recording/Independent/Start Recording/$$.
301
302$head ADVector$$
303The type $icode ADVector$$ must be a
304$cref/simple vector class/SimpleVector/$$ with elements of type
305$codei%AD<%Base%>%$$; see $cref/Base/atomic_ctor/atomic_base/Base/$$.
306
307$head afun$$
308is a $cref/atomic_user/atomic_ctor/atomic_user/$$ object
309and this $icode afun$$ function call is implemented by the
310$cref/atomic_base/atomic_ctor/atomic_base/$$ class.
311
312$head ax$$
313This argument has prototype
314$codei%
315        const %ADVector%& %ax%
316%$$
317and size must be equal to $icode n$$.
318It specifies vector $latex x \in B^n$$
319at which an $codei%AD<%Base%>%$$ version of
320$latex y = f(x)$$ is to be evaluated; see
321$cref/Base/atomic_ctor/atomic_base/Base/$$.
322
323$head ay$$
324This argument has prototype
325$codei%
326        %ADVector%& %ay%
327%$$
328and size must be equal to $icode m$$.
329The input values of its elements
330are not specified (must not matter).
331Upon return, it is an $codei%AD<%Base%>%$$ version of
332$latex y = f(x)$$.
333
334$end
335-----------------------------------------------------------------------------
336*/
337/*!
338Implement the user call to <tt>afun(ax, ay)</tt> and old_atomic call to
339<tt>afun(ax, ay, id)</tt>.
340
341\tparam ADVector
342A simple vector class with elements of type <code>AD<Base></code>.
343
344\param id
345optional extra information vector that is just passed through by CppAD,
346and used by old_atomic derived class (not other derived classes).
347This is an extra parameter to the virtual callbacks for old_atomic;
348see the set_id member function.
349
350\param ax
351is the argument vector for this call,
352<tt>ax.size()</tt> determines the number of arguments.
353
354\param ay
355is the result vector for this call,
356<tt>ay.size()</tt> determines the number of results.
357*/
358template <class ADVector>
359void operator()(
360        const ADVector&  ax     ,
361              ADVector&  ay     ,
362        size_t           id = 0 )
363{       size_t i, j;
364        size_t n = ax.size();
365        size_t m = ay.size();
366# ifndef NDEBUG
367        bool ok;
368        std::string msg = "atomic_base: " + afun_name() + ".eval: ";
369        if( (n == 0) | (m == 0) )
370        {       msg += "ax.size() or ay.size() is zero";
371                CPPAD_ASSERT_KNOWN(false, msg.c_str() );
372        }
373# endif
374        size_t thread = thread_alloc::thread_num();
375        vector <Base>& tx  = afun_tx_[thread];
376        vector <Base>& ty  = afun_ty_[thread];
377        vector <bool>& vx  = afun_vx_[thread];
378        vector <bool>& vy  = afun_vy_[thread];
379        //
380        if( vx.size() != n )
381        {       vx.resize(n);
382                tx.resize(n);
383        }
384        if( vy.size() != m )
385        {       vy.resize(m);
386                ty.resize(m);
387        }
388        //
389        // Determine tape corresponding to variables in ax
390        tape_id_t     tape_id  = 0;
391        ADTape<Base>* tape     = CPPAD_NULL;
392        for(j = 0; j < n; j++)
393        {       tx[j]  = ax[j].value_;
394                vx[j]  = Variable( ax[j] );
395                if( vx[j] )
396                {
397                        if( tape_id == 0 )
398                        {       tape    = ax[j].tape_this();
399                                tape_id = ax[j].tape_id_;
400                                CPPAD_ASSERT_UNKNOWN( tape != CPPAD_NULL );
401                        }
402# ifndef NDEBUG
403                        if( tape_id != ax[j].tape_id_ )
404                        {       msg += afun_name() + 
405                                ": ax contains variables from different threads.";
406                                CPPAD_ASSERT_KNOWN(false, msg.c_str());
407                        }
408# endif
409                }
410        }
411        // Use zero order forward mode to compute values
412        size_t q = 0, p = 0;
413        set_id(id);
414# ifdef NDEBUG
415        forward(q, p, vx, vy, tx, ty); 
416# else
417        ok = forward(q, p, vx, vy, tx, ty); 
418        if( ! ok )
419        {       msg += afun_name() + ": ok is false for "
420                        "zero order forward mode calculation.";
421                CPPAD_ASSERT_KNOWN(false, msg.c_str());
422        }
423# endif
424        bool record_operation = false;
425        for(i = 0; i < m; i++)
426        {
427                // pass back values
428                ay[i].value_ = ty[i];
429
430                // initialize entire vector parameters (not in tape)
431                ay[i].tape_id_ = 0;
432                ay[i].taddr_   = 0;
433
434                // we need to record this operation if
435                // any of the eleemnts of ay are variables,
436                record_operation |= vy[i];
437        }
438# ifndef NDEBUG
439        if( record_operation & (tape == CPPAD_NULL) )
440        {       msg += 
441                "all elements of vx are false but vy contains a true element";
442                CPPAD_ASSERT_KNOWN(false, msg.c_str() );
443        }
444# endif
445        // if tape is not null, ay is on the tape
446        if( record_operation )
447        {
448                // Operator that marks beginning of this atomic operation
449                CPPAD_ASSERT_UNKNOWN( NumRes(UserOp) == 0 );
450                CPPAD_ASSERT_UNKNOWN( NumArg(UserOp) == 4 );
451                tape->Rec_.PutArg(index_, id, n, m);
452                tape->Rec_.PutOp(UserOp);
453
454                // Now put n operators, one for each element of arugment vector
455                CPPAD_ASSERT_UNKNOWN( NumRes(UsravOp) == 0 );
456                CPPAD_ASSERT_UNKNOWN( NumRes(UsrapOp) == 0 );
457                CPPAD_ASSERT_UNKNOWN( NumArg(UsravOp) == 1 );
458                CPPAD_ASSERT_UNKNOWN( NumArg(UsrapOp) == 1 );
459                for(j = 0; j < n; j++)
460                {       if( vx[j] )
461                        {       // information for an argument that is a variable
462                                tape->Rec_.PutArg(ax[j].taddr_);
463                                tape->Rec_.PutOp(UsravOp);
464                        }
465                        else
466                        {       // information for an arugment that is parameter
467                                addr_t par = tape->Rec_.PutPar(ax[j].value_);
468                                tape->Rec_.PutArg(par);
469                                tape->Rec_.PutOp(UsrapOp);
470                        }
471                }
472
473                // Now put m operators, one for each element of result vector
474                CPPAD_ASSERT_UNKNOWN( NumArg(UsrrpOp) == 1 );
475                CPPAD_ASSERT_UNKNOWN( NumRes(UsrrpOp) == 0 );
476                CPPAD_ASSERT_UNKNOWN( NumArg(UsrrvOp) == 0 );
477                CPPAD_ASSERT_UNKNOWN( NumRes(UsrrvOp) == 1 );
478                for(i = 0; i < m; i++)
479                {       if( vy[i] )
480                        {       ay[i].taddr_    = tape->Rec_.PutOp(UsrrvOp);
481                                ay[i].tape_id_  = tape_id;
482                        }
483                        else
484                        {       addr_t par = tape->Rec_.PutPar(ay[i].value_);
485                                tape->Rec_.PutArg(par);
486                                tape->Rec_.PutOp(UsrrpOp);
487                        }
488                }
489
490                // Put a duplicate UserOp at end of UserOp sequence
491                tape->Rec_.PutArg(index_, id, n, m);
492                tape->Rec_.PutOp(UserOp);
493        } 
494        return;
495}
496/*
497-----------------------------------------------------------------------------
498$begin atomic_forward$$
499$spell
500        hes
501        afun
502        vx
503        vy
504        ty
505        Taylor
506        const
507        CppAD
508        bool
509$$
510
511$section Atomic Forward Mode$$
512$index atomic, forward callback$$
513$index forward, atomic callback$$
514$index forward, atomic virtual$$
515
516
517$head Syntax$$
518$icode%ok% = %afun%.forward(%q%, %p%, %vx%, %vy%, %tx%, %ty%)%$$
519
520$head Purpose$$
521This virtual function is used by $cref atomic_afun$$
522to evaluate function values.
523It is also used buy $cref/forward/Forward/$$
524to compute function vales and derivatives.
525
526$head Implementation$$
527This virtual function must be defined by the
528$cref/atomic_user/atomic_ctor/atomic_user/$$ class.
529It can just return $icode%ok% == false%$$
530(and not compute anything) for values
531of $icode%p% > 0%$$ that are greater than those used by your
532$cref/forward/Forward/$$ mode calculations.
533
534$head q$$
535The argument $icode q$$ has prototype
536$codei%
537        size_t %q%
538%$$
539It specifies the lowest order Taylor coefficient that we are evaluating.
540During calls to $cref atomic_afun$$, $icode%q% == 0%$$.
541
542$head p$$
543The argument $icode p$$ has prototype
544$codei%
545        size_t %p%
546%$$
547It specifies the highest order Taylor coefficient that we are evaluating.
548During calls to $cref atomic_afun$$, $icode%p% == 0%$$.
549
550$head vx$$
551The $code forward$$ argument $icode vx$$ has prototype
552$codei%
553        const CppAD::vector<bool>& %vx%
554%$$
555The case $icode%vx%.size() > 0%$$ only occurs while evaluating a call to
556$cref atomic_afun$$.
557In this case,
558$icode%q% == %p% == 0%$$,
559$icode%vx%.size() == %n%$$, and
560for $latex j = 0 , \ldots , n-1$$,
561$icode%vx%[%j%]%$$ is true if and only if
562$icode%ax%[%j%]%$$ is a $cref/variable/glossary/Variable/$$
563in the corresponding call to
564$codei%
565        %afun%(%ax%, %ay%, %id%)
566%$$
567This is only necessary to use for second order cross terms because
568CppAD knows which components are variables; e.g., see
569$code rev_sparse_hes$$ in the $cref atomic_matrix_mul.hpp$$ example.
570$pre
571
572$$
573If $icode%vx%.size() == 0%$$,
574then $icode%vy%.size() == 0%$$ and neither of these vectors
575should be used.
576
577$head vy$$
578The $code forward$$ argument $icode vy$$ has prototype
579$codei%
580        CppAD::vector<bool>& %vy%
581%$$
582If $icode%vy%.size() == 0%$$, it should not be used.
583Otherwise,
584$icode%p% == 0%$$ and $icode%vy%.size() == %m%$$.
585The input values of the elements of $icode vy$$
586are not specified (must not matter).
587Upon return, for $latex j = 0 , \ldots , m-1$$,
588$icode%vy%[%i%]%$$ is true if and only if
589$icode%ay%[%i%]%$$ is a variable
590(CppAD uses $icode vy$$ to reduce the necessary computations).
591
592$head tx$$
593The argument $icode tx$$ has prototype
594$codei%
595        const CppAD::vector<%Base%>& %tx%
596%$$
597and $icode%tx%.size() == (%p%+1)*%n%$$.
598For $latex j = 0 , \ldots , n-1$$ and $latex k = 0 , \ldots , p$$,
599we use the Taylor coefficient notation
600$latex \[
601\begin{array}{rcl}
602        x_j^k    & = & tx [ j * ( p + 1 ) + k ]
603        \\
604        X_j (t)  & = & x_j^0 + x_j^1 t^1 + \cdots + x_j^p t^p
605\end{array}
606\] $$
607Note that superscripts represent an index for $latex x_j^k$$
608and an exponent for $latex t^k$$.
609Also note that the Taylor coefficients for $latex X(t)$$ correspond
610to the derivatives of $latex X(t)$$ at $latex t = 0$$ in the following way:
611$latex \[
612        x_j^k = \frac{1}{ k ! } X_j^{(k)} (0)
613\] $$
614
615$head ty$$
616The argument $icode ty$$ has prototype
617$codei%
618        CppAD::vector<%Base%>& %ty%
619%$$
620and $icode%tx%.size() == (%p%+1)*%m%$$.
621Upon return,
622For $latex i = 0 , \ldots , m-1$$ and $latex k = 0 , \ldots , p$$,
623$latex \[
624\begin{array}{rcl}
625        Y_i (t)  & = & f_i [ X(t) ]
626        \\
627        Y_i (t)  & = & y_i^0 + y_i^1 t^1 + \cdots + y_i^p t^p + o ( t^p )
628        \\
629        ty [ i * ( p + 1 ) + k ] & = & y_i^k
630\end{array}
631\] $$
632where $latex o( t^p ) / t^p \rightarrow 0$$ as $latex t \rightarrow 0$$.
633Note that superscripts represent an index for $latex y_j^k$$
634and an exponent for $latex t^k$$.
635Also note that the Taylor coefficients for $latex Y(t)$$ correspond
636to the derivatives of $latex Y(t)$$ at $latex t = 0$$ in the following way:
637$latex \[
638        y_j^k = \frac{1}{ k ! } Y_j^{(k)} (0)
639\] $$
640If $latex q > 0$$,
641for $latex i = 0 , \ldots , m-1$$ and $latex k = 0 , \ldots , q-1$$,
642the input of $icode ty$$ satisfies
643$latex \[
644        ty [ i * ( p + 1 ) + k ] = y_i^k
645\]$$
646and hence the corresponding elements need not be recalculated.
647
648$head ok$$
649If the required results are calculated, $icode ok$$ should be true.
650Otherwise, it should be false.
651
652$head Example$$
653For example, suppose that $icode%p% == 2%$$,
654and you know how to compute the function $latex f(x)$$,
655its first derivative $latex f^{(1)} (x)$$,
656and it component wise Hessian $latex f_i^{(2)} (x)$$.
657Then you can compute $icode ty$$ using the following formulas:
658$latex \[
659\begin{array}{rcl}
660y_i^0 & = & Y(0)
661        = f_i ( x^0 )
662\\
663y_i^1 & = & Y^{(1)} ( 0 )
664        = f_i^{(1)} ( x^0 ) X^{(1)} ( 0 )
665        = f_i^{(1)} ( x^0 ) x^1
666\\
667y_i^2
668& = & \frac{1}{2 !} Y^{(2)} (0)
669\\
670& = & \frac{1}{2} X^{(1)} (0)^\R{T} f_i^{(2)} ( x^0 ) X^{(1)} ( 0 )
671  +   \frac{1}{2} f_i^{(1)} ( x^0 ) X^{(2)} ( 0 )
672\\
673& = & \frac{1}{2} (x^1)^\R{T} f_i^{(2)} ( x^0 ) x^1
674  +    f_i^{(1)} ( x^0 ) x^2
675\end{array}
676\] $$
677For $latex i = 0 , \ldots , m-1$$, and $latex k = 0 , 1 , 2$$,
678$latex \[
679        ty [ i * (p + 1) + k ] = y_i^k
680\] $$
681 
682$end
683-----------------------------------------------------------------------------
684*/
685/*!
686Link from atomic_base to forward mode
687
688\param q [in]
689lowerest order for this forward mode calculation.
690
691\param p [in]
692highest order for this forward mode calculation.
693
694\param vx [in]
695if size not zero, which components of \c x are variables
696
697\param vy [out]
698if size not zero, which components of \c y are variables
699
700\param tx [in]
701Taylor coefficients corresponding to \c x for this calculation.
702
703\param ty [out]
704Taylor coefficient corresponding to \c y for this calculation
705
706See the forward mode in user's documentation for base_atomic
707*/
708virtual bool forward(
709        size_t                    q  ,
710        size_t                    p  ,
711        const vector<bool>&       vx ,
712              vector<bool>&       vy ,
713        const vector<Base>&       tx ,
714              vector<Base>&       ty )
715{       return false; }
716/*
717-----------------------------------------------------------------------------
718$begin atomic_reverse$$
719$spell
720        afun
721        ty
722        px
723        py
724        Taylor
725        const
726        CppAD
727$$
728
729$section Atomic Reverse Mode$$
730$index atomic, reverse callback$$
731$index reverse, atomic callback$$
732$index reverse, atomic virtual$$
733$spell
734        bool
735$$
736
737$head Syntax$$
738$icode%ok% = %afun%.reverse(%p%, %tx%, %ty%, %px%, %py%)%$$
739
740$head Purpose$$
741This function is used by $cref/reverse/Reverse/$$
742to compute derivatives.
743
744$head Implementation$$
745If you are using
746$cref/reverse/Reverse/$$ mode,
747this virtual function must be defined by the
748$cref/atomic_user/atomic_ctor/atomic_user/$$ class.
749It can just return $icode%ok% == false%$$
750(and not compute anything) for values
751of $icode p$$ that are greater than those used by your
752$cref/reverse/Reverse/$$ mode calculations.
753
754$head p$$
755The argument $icode p$$ has prototype
756$codei%
757        size_t %p%
758%$$
759It specifies the highest order Taylor coefficient that
760computing the derivative of.
761
762$head tx$$
763The argument $icode tx$$ has prototype
764$codei%
765        const CppAD::vector<%Base%>& %tx%
766%$$
767and $icode%tx%.size() == (%p%+1)*%n%$$.
768For $latex j = 0 , \ldots , n-1$$ and $latex k = 0 , \ldots , p$$,
769we use the Taylor coefficient notation
770$latex \[
771\begin{array}{rcl}
772        x_j^k    & = & tx [ j * ( p + 1 ) + k ]
773        \\
774        X_j (t)  & = & x_j^0 + x_j^1 t^1 + \cdots + x_j^p t^p
775\end{array}
776\] $$
777Note that superscripts represent an index for $latex x_j^k$$
778and an exponent for $latex t^k$$.
779Also note that the Taylor coefficients for $latex X(t)$$ correspond
780to the derivatives of $latex X(t)$$ at $latex t = 0$$ in the following way:
781$latex \[
782        x_j^k = \frac{1}{ k ! } X_j^{(k)} (0)
783\] $$
784
785$head ty$$
786The argument $icode ty$$ has prototype
787$codei%
788        const CppAD::vector<%Base%>& %ty%
789%$$
790and $icode%tx%.size() == (%p%+1)*%m%$$.
791For $latex i = 0 , \ldots , m-1$$ and $latex k = 0 , \ldots , p$$,
792we use the Taylor coefficient notation
793$latex \[
794\begin{array}{rcl}
795        Y_i (t)  & = & f_i [ X(t) ]
796        \\
797        Y_i (t)  & = & y_i^0 + y_i^1 t^1 + \cdots + y_i^p t^p + o ( t^p )
798        \\
799        y_i^k    & = & ty [ i * ( p + 1 ) + k ]
800\end{array}
801\] $$
802where $latex o( t^p ) / t^p \rightarrow 0$$ as $latex t \rightarrow 0$$.
803Note that superscripts represent an index for $latex y_j^k$$
804and an exponent for $latex t^k$$.
805Also note that the Taylor coefficients for $latex Y(t)$$ correspond
806to the derivatives of $latex Y(t)$$ at $latex t = 0$$ in the following way:
807$latex \[
808        y_j^k = \frac{1}{ k ! } Y_j^{(k)} (0)
809\] $$
810
811
812$head F, G, H$$
813We use the notation $latex \{ x_j^k \} \in B^{n \times (p+1)}$$ for
814$latex \[
815        \{ x_j^k \W{:} j = 0 , \ldots , n-1, k = 0 , \ldots , p \}
816\]$$
817We use the notation $latex \{ y_i^k \} \in B^{m \times (p+1)}$$ for
818$latex \[
819        \{ y_i^k \W{:} i = 0 , \ldots , m-1, k = 0 , \ldots , p \}
820\]$$
821We define the function
822$latex F : B^{n \times (p+1)} \rightarrow B^{m \times (p+1)}$$ by
823$latex \[
824        y_i^k = F_i^k [ \{ x_j^k \} ]
825\] $$
826We use $latex G : B^{m \times (p+1)} \rightarrow B$$
827to denote an arbitrary scalar valued function of $latex \{ y_i^k \}$$.
828We use $latex H : B^{n \times (p+1)} \rightarrow B$$
829defined by
830$latex \[
831        H ( \{ x_j^k \} ) = G[ F( \{ x_j^k \} ) ]
832\] $$
833
834$head py$$
835The argument $icode py$$ has prototype
836$codei%
837        const CppAD::vector<%Base%>& %py%
838%$$
839and $icode%py%.size() == m * (%p%+1)%$$.
840For $latex i = 0 , \ldots , m-1$$, $latex k = 0 , \ldots , p$$,
841$latex \[
842        py[ i * (p + 1 ) + k ] = \partial G / \partial y_i^k
843\] $$
844
845$subhead px$$
846The $icode px$$ has prototype
847$codei%
848        CppAD::vector<%Base%>& %px%
849%$$
850and $icode%px%.size() == n * (%p%+1)%$$.
851The input values of the elements of $icode px$$
852are not specified (must not matter).
853Upon return,
854for $latex j = 0 , \ldots , n-1$$ and $latex \ell = 0 , \ldots , p$$,
855$latex \[
856\begin{array}{rcl}
857px [ j * (p + 1) + \ell ] & = & \partial H / \partial x_j^\ell
858\\
859& = &
860( \partial G / \partial \{ y_i^k \} )
861        ( \partial \{ y_i^k \} / \partial x_j^\ell )
862\\
863& = &
864\sum_{i=0}^{m-1} \sum_{k=0}^p
865( \partial G / \partial y_i^k ) ( \partial y_i^k / \partial x_j^\ell )
866\\
867& = &
868\sum_{i=0}^{m-1} \sum_{k=\ell}^p
869py[ i * (p + 1 ) + k ] ( \partial F_i^k / \partial x_j^\ell )
870\end{array}
871\] $$
872Note that we have used the fact that for $latex k < \ell$$,
873$latex \partial F_i^k / \partial x_j^\ell = 0$$.
874
875$head ok$$
876The return value $icode ok$$ has prototype
877$codei%
878        bool %ok%
879%$$
880If it is $code true$$, the corresponding evaluation succeeded,
881otherwise it failed.
882
883$end
884-----------------------------------------------------------------------------
885*/
886/*!
887Link from reverse mode sweep to users routine.
888
889\param p [in]
890highest order for this reverse mode calculation.
891
892\param tx [in]
893Taylor coefficients corresponding to \c x for this calculation.
894
895\param ty [in]
896Taylor coefficient corresponding to \c y for this calculation
897
898\param px [out]
899Partials w.r.t. the \c x Taylor coefficients.
900
901\param py [in]
902Partials w.r.t. the \c y Taylor coefficients.
903
904See atomic_reverse mode use documentation
905*/
906virtual bool reverse(
907        size_t                    p  ,
908        const vector<Base>&       tx ,
909        const vector<Base>&       ty ,
910              vector<Base>&       px ,
911        const vector<Base>&       py )
912{       return false; }
913/*
914-------------------------------------- ---------------------------------------
915$begin atomic_for_sparse_jac$$
916$spell
917        afun
918        Jacobian
919        jac
920        const
921        CppAD
922        std
923        bool
924        std
925$$
926
927$section Atomic Forward Jacobian Sparsity Patterns$$
928$index atomic, for_sparse_jac callback$$
929$index for_sparse_jac, atomic callback$$
930$index for_sparse_jac, atomic virtual$$
931
932$head Syntax$$
933$icode%ok% = %afun%.for_sparse_jac(%q%, %r%, %s%)%$$
934
935$head Purpose$$
936This function is used by $cref ForSparseJac$$ to compute
937Jacobian sparsity patterns.
938For a fixed matrix $latex R \in B^{n \times q}$$,
939the Jacobian of $latex f( x + R * u)$$ with respect to $latex u \in B^q$$ is
940$latex \[
941        S(x) = f^{(1)} (x) * R
942\] $$
943Given a $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for $latex R$$,
944$code for_sparse_jac$$ computes a sparsity pattern for $latex S(x)$$.
945
946$head Implementation$$
947If you are using $cref ForSparseJac$$,
948this virtual function must be defined by the
949$cref/atomic_user/atomic_ctor/atomic_user/$$ class.
950
951$subhead q$$
952The argument $icode q$$ has prototype
953$codei%
954     size_t %q%
955%$$
956It specifies the number of columns in
957$latex R \in B^{n \times q}$$ and the Jacobian
958$latex S(x) \in B^{m \times q}$$.
959
960$subhead r$$
961This argument has prototype
962$codei%
963     const %atomic_sparsity%& %r%
964%$$
965and is a $cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
966$latex R \in B^{n \times q}$$.
967
968$subhead s$$
969This argument has prototype
970$codei%
971        %atomic_sparsity%& %s%
972%$$
973The input values of its elements
974are not specified (must not matter).
975Upon return, $icode s$$ is a
976$cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
977$latex S(x) \in B^{m \times q}$$.
978
979$head ok$$
980The return value $icode ok$$ has prototype
981$codei%
982        bool %ok%
983%$$
984If it is $code true$$, the corresponding evaluation succeeded,
985otherwise it failed.
986
987$end
988-----------------------------------------------------------------------------
989*/
990/*!
991Link from forward Jacobian sparsity sweep to atomic_base
992
993\param q
994is the column dimension for the Jacobian sparsity partterns.
995
996\param r
997is the Jacobian sparsity pattern for the argument vector x
998
999\param s
1000is the Jacobian sparsity pattern for the result vector y
1001*/
1002virtual bool for_sparse_jac(
1003        size_t                                  q  ,
1004        const vector< std::set<size_t> >&       r  ,
1005              vector< std::set<size_t> >&       s  )
1006{       return false; }
1007virtual bool for_sparse_jac(
1008        size_t                                  q  ,
1009        const vector<bool>&                     r  ,
1010              vector<bool>&                     s  )
1011{       return false; }
1012/*
1013-------------------------------------- ---------------------------------------
1014$begin atomic_rev_sparse_jac$$
1015$spell
1016        rt
1017        afun
1018        Jacobian
1019        jac
1020        CppAD
1021        std
1022        bool
1023        const
1024        hes
1025$$
1026
1027$section Atomic Reverse Jacobian Sparsity Patterns$$
1028$index atomic, rev_sparse_jac callback$$
1029$index rev_sparse_jac, atomic callback$$
1030$index rev_sparse_jac, atomic virtual$$
1031
1032$head Syntax$$
1033$icode%ok% = %afun%.rev_sparse_jac(%q%, %rt%, %st%)%$$
1034
1035$head Purpose$$
1036This function is used by $cref RevSparseJac$$ to compute
1037Jacobian sparsity patterns.
1038For a fixed matrix $latex R \in B^{q \times m}$$,
1039the Jacobian of $latex R * f( x )$$ with respect to $latex x \in B^q$$ is
1040$latex \[
1041        S(x) = R * f^{(1)} (x)
1042\] $$
1043Given a $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for $latex R$$,
1044$code rev_sparse_jac$$ computes a sparsity pattern for $latex S(x)$$.
1045
1046$head Implementation$$
1047If you are using $cref RevSparseJac$$,
1048this virtual function must be defined by the
1049$cref/atomic_user/atomic_ctor/atomic_user/$$ class.
1050
1051$subhead q$$
1052The argument $icode q$$ has prototype
1053$codei%
1054     size_t %q%
1055%$$
1056It specifies the number of rows in
1057$latex R \in B^{q \times m}$$ and the Jacobian
1058$latex S(x) \in B^{q \times n}$$.
1059
1060$subhead rt$$
1061This argument has prototype
1062$codei%
1063     const %atomic_sparsity%& %rt%
1064%$$
1065and is a
1066$cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
1067$latex R^\R{T} \in B^{m \times q}$$.
1068
1069$subhead st$$
1070This argument has prototype
1071$codei%
1072        %atomic_sparsity%& %st%
1073%$$
1074The input value of its elements
1075are not specified (must not matter).
1076Upon return, $icode s$$ is a
1077$cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
1078$latex S(x)^\R{T} \in B^{n \times q}$$.
1079
1080$head ok$$
1081The return value $icode ok$$ has prototype
1082$codei%
1083        bool %ok%
1084%$$
1085If it is $code true$$, the corresponding evaluation succeeded,
1086otherwise it failed.
1087
1088$end
1089-----------------------------------------------------------------------------
1090*/
1091/*!
1092Link from reverse Jacobian sparsity sweep to atomic_base
1093
1094\param q [in]
1095is the row dimension for the Jacobian sparsity partterns
1096
1097\param rt [out]
1098is the tansposed Jacobian sparsity pattern w.r.t to range variables y
1099
1100\param st [in]
1101is the tansposed Jacobian sparsity pattern for the argument variables x
1102*/
1103virtual bool rev_sparse_jac(
1104        size_t                                  q  ,
1105        const vector< std::set<size_t> >&       rt ,
1106              vector< std::set<size_t> >&       st )
1107{       return false; }
1108virtual bool rev_sparse_jac(
1109        size_t                                  q  ,
1110        const vector<bool>&                     rt ,
1111              vector<bool>&                     st )
1112{       return false; }
1113/*
1114-------------------------------------- ---------------------------------------
1115$begin atomic_rev_sparse_hes$$
1116$spell
1117        vx
1118        afun
1119        Jacobian
1120        jac
1121        CppAD
1122        std
1123        bool
1124        hes
1125        const
1126$$
1127
1128$section Atomic Reverse Hessian Sparsity Patterns$$
1129$index atomic, rev_sparse_hes callback$$
1130$index rev_sparse_hes, atomic callback$$
1131$index rev_sparse_hes, atomic virtual$$
1132
1133$head Syntax$$
1134$icode%ok% = %afun%.rev_sparse_hes(%vx%, %s%, %t%, %q%, %r%, %u%, %v%)%$$
1135
1136$head Purpose$$
1137This function is used by $cref RevSparseHes$$ to compute
1138Hessian sparsity patterns.
1139There is an unspecified scalar valued function
1140$latex g : B^m \rightarrow B$$.
1141Given a $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for
1142$latex R \in B^{n \times q}$$,
1143and information about the function $latex z = g(y)$$,
1144this routine computes the sparsity pattern for
1145$latex \[
1146        V(x) = (g \circ f)^{(2)}( x ) R
1147\] $$
1148
1149$head Implementation$$
1150If you are using and $cref RevSparseHes$$,
1151this virtual function must be defined by the
1152$cref/atomic_user/atomic_ctor/atomic_user/$$ class.
1153
1154$subhead vx$$
1155The argument $icode vx$$ has prototype
1156$codei%
1157     const CppAD:vector<bool>& %vx%
1158%$$
1159$icode%vx%.size() == %n%$$, and
1160for $latex j = 0 , \ldots , n-1$$,
1161$icode%vx%[%j%]%$$ is true if and only if
1162$icode%ax%[%j%]%$$ is a $cref/variable/glossary/Variable/$$
1163in the corresponding call to
1164$codei%
1165        %afun%(%ax%, %ay%, %id%)
1166%$$
1167
1168$subhead s$$
1169The argument $icode s$$ has prototype
1170$codei%
1171     const CppAD:vector<bool>& %s%
1172%$$
1173and its size is $icode m$$.
1174It is a sparsity pattern for
1175$latex S(x) = g^{(1)} (y) \in B^{1 \times m}$$.
1176
1177$subhead t$$
1178This argument has prototype
1179$codei%
1180     CppAD:vector<bool>& %t%
1181%$$
1182and its size is $icode m$$.
1183The input values of its elements
1184are not specified (must not matter).
1185Upon return, $icode t$$ is a
1186sparsity pattern for
1187$latex T(x) \in B^{1 \times n}$$ where
1188$latex \[
1189        T(x) = (g \circ f)^{(1)} (x) = S(x) * f^{(1)} (x)
1190\]$$
1191
1192$subhead q$$
1193The argument $icode q$$ has prototype
1194$codei%
1195     size_t %q%
1196%$$
1197It specifies the number of columns in
1198$latex R \in B^{n \times q}$$,
1199$latex U(x) \in B^{m \times q}$$, and
1200$latex V(x) \in B^{n \times q}$$.
1201
1202$subhead r$$
1203This argument has prototype
1204$codei%
1205     const %atomic_sparsity%& %r%
1206%$$
1207and is a $cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
1208$latex R \in B^{n \times q}$$.
1209
1210$head u$$
1211This argument has prototype
1212$codei%
1213     const %atomic_sparsity%& %u%
1214%$$
1215and is a $cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
1216$latex U(x) \in B^{m \times q}$$ which is defined by
1217$latex \[
1218\begin{array}{rcl}
1219U(x)
1220& = &
1221\partial_u \{ \partial_y g[ y + f^{(1)} (x) R u ] \}_{u=0}
1222\\
1223& = &
1224\partial_u \{ g^{(1)} [ y + f^{(1)} (x) R u ] \}_{u=0}
1225\\
1226& = &
1227g^{(2)} (y) f^{(1)} (x) R
1228\end{array}
1229\] $$
1230
1231$subhead v$$
1232This argument has prototype
1233$codei%
1234     %atomic_sparsity%& %v%
1235%$$
1236The input value of its elements
1237are not specified (must not matter).
1238Upon return, $icode v$$ is a
1239$cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
1240$latex V(x) \in B^{n \times q}$$ which is defined by
1241$latex \[
1242\begin{array}{rcl}
1243V(x)
1244& = &
1245\partial_u [ \partial_x (g \circ f) ( x + R u )  ]_{u=0}
1246\\
1247& = &
1248\partial_u [ (g \circ f)^{(1)}( x + R u )  ]_{u=0}
1249\\
1250& = &
1251(g \circ f)^{(2)}( x ) R
1252\\
1253& = &
1254f^{(1)} (x)^\R{T} g^{(2)} ( y ) f^{(1)} (x)  R
1255+
1256\sum_{i=1}^m g_i^{(1)} (y) \; f_i^{(2)} (x) R
1257\\
1258& = &
1259f^{(1)} (x)^\R{T} U(x)
1260+
1261\sum_{i=1}^m S_i (x) \; f_i^{(2)} (x) R
1262\end{array}
1263\] $$
1264
1265$end
1266-----------------------------------------------------------------------------
1267*/
1268/*!
1269Link from reverse Hessian sparsity sweep to base_atomic
1270
1271\param vx [in]
1272which componens of x are variables.
1273
1274\param s [in]
1275is the reverse Jacobian sparsity pattern w.r.t the result vector y.
1276
1277\param t [out]
1278is the reverse Jacobian sparsity pattern w.r.t the argument vector x.
1279
1280\param q [in]
1281is the column dimension for the sparsity partterns.
1282
1283\param r [in]
1284is the forward Jacobian sparsity pattern w.r.t the argument vector x
1285
1286\param u [in]
1287is the Hessian sparsity pattern w.r.t the result vector y.
1288
1289\param v [out]
1290is the Hessian sparsity pattern w.r.t the argument vector x.
1291*/
1292virtual bool rev_sparse_hes(
1293        const vector<bool>&                     vx ,
1294        const vector<bool>&                     s  ,
1295              vector<bool>&                     t  ,
1296        size_t                                  q  ,
1297        const vector< std::set<size_t> >&       r  ,
1298        const vector< std::set<size_t> >&       u  ,
1299              vector< std::set<size_t> >&       v  )
1300{       return false; }
1301virtual bool rev_sparse_hes(
1302        const vector<bool>&                     vx ,
1303        const vector<bool>&                     s  ,
1304              vector<bool>&                     t  ,
1305        size_t                                  q  ,
1306        const vector<bool>&                     r  ,
1307        const vector<bool>&                     u  ,
1308              vector<bool>&                     v  )
1309{       return false; }
1310/*
1311------------------------------------------------------------------------------
1312$begin atomic_base_clear$$
1313
1314$section Free Static Variables$$
1315$index free, atomic static$$
1316$index atomic, free static$$
1317$index static, free atomic$$
1318$index clear, atomic static$$
1319
1320$head Syntax$$
1321$codei%atomic_base<%Base%>::clear()%$$
1322
1323$head Purpose$$
1324The $code atomic_base$$ class holds onto static work space in order to
1325increase speed by avoiding system memory allocation calls.
1326This call makes to work space $cref/available/ta_available/$$ to
1327for other uses by the same thread.
1328This should be called when you are done using the
1329user atomic functions for a specific value of $icode Base$$.
1330
1331$head Restriction$$
1332This routine cannot be called
1333while in $cref/parallel/ta_in_parallel/$$ execution mode.
1334
1335$end
1336------------------------------------------------------------------------------
1337*/
1338/*!
1339Free all thread_alloc static memory held by atomic_base (avoids reallocations).
1340(This does not include class_object() which is an std::vector.)
1341*/
1342/// Free vector memory used by this class (work space)
1343static void clear(void)
1344{       CPPAD_ASSERT_KNOWN(
1345                ! thread_alloc::in_parallel() ,
1346                "cannot use atomic_base clear during parallel execution"
1347        );
1348        size_t i = class_object().size();
1349        while(i--)
1350        {       size_t thread = CPPAD_MAX_NUM_THREADS;
1351                while(thread--)
1352                {
1353                        atomic_base* op = class_object()[i];
1354                        if( op != CPPAD_NULL )
1355                        {       op->afun_vx_[thread].clear();
1356                                op->afun_vy_[thread].clear();
1357                                op->afun_tx_[thread].clear();
1358                                op->afun_ty_[thread].clear();
1359                        }
1360                }
1361        }
1362        return;
1363}
1364// -------------------------------------------------------------------------
1365/*!
1366Set value of id (used by deprecated old_atomic class)
1367
1368This function is called just before calling any of the virtual funcitons
1369and has the corresponding id of the corresponding virtual call.
1370*/
1371virtual void set_id(size_t id) 
1372{ }
1373// ---------------------------------------------------------------------------
1374};
1375/*! \} */
1376CPPAD_END_NAMESPACE
1377# endif
Note: See TracBrowser for help on using the repository browser.