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

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

atomic_base.hpp: remove documentaiton text used in old atomic functions.
get_started.cpp: remove executable property (set by mistake).

  • Property svn:keywords set to Id
File size: 36.1 KB
Line 
1/* $Id: atomic_base.hpp 2902 2013-09-19 12:19:11Z 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%$$
567If $icode%vx%.size() == 0%$$,
568then $icode%vy%.size() == 0%$$ and neither of these vectors
569should be used.
570
571$head vy$$
572The $code forward$$ argument $icode vy$$ has prototype
573$codei%
574        CppAD::vector<bool>& %vy%
575%$$
576If $icode%vy%.size() == 0%$$, it should not be used.
577Otherwise,
578$icode%p% == 0%$$ and $icode%vy%.size() == %m%$$.
579The input values of the elements of $icode vy$$
580are not specified (must not matter).
581Upon return, for $latex j = 0 , \ldots , m-1$$,
582$icode%vy%[%i%]%$$ is true if and only if
583$icode%ay%[%i%]%$$ is a variable
584(CppAD uses $icode vy$$ to reduce the necessary computations).
585
586$head tx$$
587The argument $icode tx$$ has prototype
588$codei%
589        const CppAD::vector<%Base%>& %tx%
590%$$
591and $icode%tx%.size() == (%p%+1)*%n%$$.
592For $latex j = 0 , \ldots , n-1$$ and $latex k = 0 , \ldots , p$$,
593we use the Taylor coefficient notation
594$latex \[
595\begin{array}{rcl}
596        x_j^k    & = & tx [ j * ( p + 1 ) + k ]
597        \\
598        X_j (t)  & = & x_j^0 + x_j^1 t^1 + \cdots + x_j^p t^p
599\end{array}
600\] $$
601Note that superscripts represent an index for $latex x_j^k$$
602and an exponent for $latex t^k$$.
603Also note that the Taylor coefficients for $latex X(t)$$ correspond
604to the derivatives of $latex X(t)$$ at $latex t = 0$$ in the following way:
605$latex \[
606        x_j^k = \frac{1}{ k ! } X_j^{(k)} (0)
607\] $$
608
609$head ty$$
610The argument $icode ty$$ has prototype
611$codei%
612        CppAD::vector<%Base%>& %ty%
613%$$
614and $icode%tx%.size() == (%p%+1)*%m%$$.
615Upon return,
616For $latex i = 0 , \ldots , m-1$$ and $latex k = 0 , \ldots , p$$,
617$latex \[
618\begin{array}{rcl}
619        Y_i (t)  & = & f_i [ X(t) ]
620        \\
621        Y_i (t)  & = & y_i^0 + y_i^1 t^1 + \cdots + y_i^p t^p + o ( t^p )
622        \\
623        ty [ i * ( p + 1 ) + k ] & = & y_i^k
624\end{array}
625\] $$
626where $latex o( t^p ) / t^p \rightarrow 0$$ as $latex t \rightarrow 0$$.
627Note that superscripts represent an index for $latex y_j^k$$
628and an exponent for $latex t^k$$.
629Also note that the Taylor coefficients for $latex Y(t)$$ correspond
630to the derivatives of $latex Y(t)$$ at $latex t = 0$$ in the following way:
631$latex \[
632        y_j^k = \frac{1}{ k ! } Y_j^{(k)} (0)
633\] $$
634If $latex q > 0$$,
635for $latex i = 0 , \ldots , m-1$$ and $latex k = 0 , \ldots , q-1$$,
636the input of $icode ty$$ satisfies
637$latex \[
638        ty [ i * ( p + 1 ) + k ] = y_i^k
639\]$$
640and hence the corresponding elements need not be recalculated.
641
642$head ok$$
643If the required results are calculated, $icode ok$$ should be true.
644Otherwise, it should be false.
645
646$head Example$$
647For example, suppose that $icode%p% == 2%$$,
648and you know how to compute the function $latex f(x)$$,
649its first derivative $latex f^{(1)} (x)$$,
650and it component wise Hessian $latex f_i^{(2)} (x)$$.
651Then you can compute $icode ty$$ using the following formulas:
652$latex \[
653\begin{array}{rcl}
654y_i^0 & = & Y(0)
655        = f_i ( x^0 )
656\\
657y_i^1 & = & Y^{(1)} ( 0 )
658        = f_i^{(1)} ( x^0 ) X^{(1)} ( 0 )
659        = f_i^{(1)} ( x^0 ) x^1
660\\
661y_i^2
662& = & \frac{1}{2 !} Y^{(2)} (0)
663\\
664& = & \frac{1}{2} X^{(1)} (0)^\R{T} f_i^{(2)} ( x^0 ) X^{(1)} ( 0 )
665  +   \frac{1}{2} f_i^{(1)} ( x^0 ) X^{(2)} ( 0 )
666\\
667& = & \frac{1}{2} (x^1)^\R{T} f_i^{(2)} ( x^0 ) x^1
668  +    f_i^{(1)} ( x^0 ) x^2
669\end{array}
670\] $$
671For $latex i = 0 , \ldots , m-1$$, and $latex k = 0 , 1 , 2$$,
672$latex \[
673        ty [ i * (p + 1) + k ] = y_i^k
674\] $$
675 
676$end
677-----------------------------------------------------------------------------
678*/
679/*!
680Link from atomic_base to forward mode
681
682\param q [in]
683lowerest order for this forward mode calculation.
684
685\param p [in]
686highest order for this forward mode calculation.
687
688\param vx [in]
689if size not zero, which components of \c x are variables
690
691\param vy [out]
692if size not zero, which components of \c y are variables
693
694\param tx [in]
695Taylor coefficients corresponding to \c x for this calculation.
696
697\param ty [out]
698Taylor coefficient corresponding to \c y for this calculation
699
700See the forward mode in user's documentation for base_atomic
701*/
702virtual bool forward(
703        size_t                    q  ,
704        size_t                    p  ,
705        const vector<bool>&       vx ,
706              vector<bool>&       vy ,
707        const vector<Base>&       tx ,
708              vector<Base>&       ty )
709{       return false; }
710/*
711-----------------------------------------------------------------------------
712$begin atomic_reverse$$
713$spell
714        afun
715        ty
716        px
717        py
718        Taylor
719        const
720        CppAD
721$$
722
723$section Atomic Reverse Mode$$
724$index atomic, reverse callback$$
725$index reverse, atomic callback$$
726$index reverse, atomic virtual$$
727$spell
728        bool
729$$
730
731$head Syntax$$
732$icode%ok% = %afun%.reverse(%p%, %tx%, %ty%, %px%, %py%)%$$
733
734$head Purpose$$
735This function is used by $cref/reverse/Reverse/$$
736to compute derivatives.
737
738$head Implementation$$
739If you are using
740$cref/reverse/Reverse/$$ mode,
741this virtual function must be defined by the
742$cref/atomic_user/atomic_ctor/atomic_user/$$ class.
743It can just return $icode%ok% == false%$$
744(and not compute anything) for values
745of $icode p$$ that are greater than those used by your
746$cref/reverse/Reverse/$$ mode calculations.
747
748$head p$$
749The argument $icode p$$ has prototype
750$codei%
751        size_t %p%
752%$$
753It specifies the highest order Taylor coefficient that
754computing the derivative of.
755
756$head tx$$
757The argument $icode tx$$ has prototype
758$codei%
759        const CppAD::vector<%Base%>& %tx%
760%$$
761and $icode%tx%.size() == (%p%+1)*%n%$$.
762For $latex j = 0 , \ldots , n-1$$ and $latex k = 0 , \ldots , p$$,
763we use the Taylor coefficient notation
764$latex \[
765\begin{array}{rcl}
766        x_j^k    & = & tx [ j * ( p + 1 ) + k ]
767        \\
768        X_j (t)  & = & x_j^0 + x_j^1 t^1 + \cdots + x_j^p t^p
769\end{array}
770\] $$
771Note that superscripts represent an index for $latex x_j^k$$
772and an exponent for $latex t^k$$.
773Also note that the Taylor coefficients for $latex X(t)$$ correspond
774to the derivatives of $latex X(t)$$ at $latex t = 0$$ in the following way:
775$latex \[
776        x_j^k = \frac{1}{ k ! } X_j^{(k)} (0)
777\] $$
778
779$head ty$$
780The argument $icode ty$$ has prototype
781$codei%
782        const CppAD::vector<%Base%>& %ty%
783%$$
784and $icode%tx%.size() == (%p%+1)*%m%$$.
785For $latex i = 0 , \ldots , m-1$$ and $latex k = 0 , \ldots , p$$,
786we use the Taylor coefficient notation
787$latex \[
788\begin{array}{rcl}
789        Y_i (t)  & = & f_i [ X(t) ]
790        \\
791        Y_i (t)  & = & y_i^0 + y_i^1 t^1 + \cdots + y_i^p t^p + o ( t^p )
792        \\
793        y_i^k    & = & ty [ i * ( p + 1 ) + k ]
794\end{array}
795\] $$
796where $latex o( t^p ) / t^p \rightarrow 0$$ as $latex t \rightarrow 0$$.
797Note that superscripts represent an index for $latex y_j^k$$
798and an exponent for $latex t^k$$.
799Also note that the Taylor coefficients for $latex Y(t)$$ correspond
800to the derivatives of $latex Y(t)$$ at $latex t = 0$$ in the following way:
801$latex \[
802        y_j^k = \frac{1}{ k ! } Y_j^{(k)} (0)
803\] $$
804
805
806$head F, G, H$$
807We use the notation $latex \{ x_j^k \} \in B^{n \times (p+1)}$$ for
808$latex \[
809        \{ x_j^k \W{:} j = 0 , \ldots , n-1, k = 0 , \ldots , p \}
810\]$$
811We use the notation $latex \{ y_i^k \} \in B^{m \times (p+1)}$$ for
812$latex \[
813        \{ y_i^k \W{:} i = 0 , \ldots , m-1, k = 0 , \ldots , p \}
814\]$$
815We define the function
816$latex F : B^{n \times (p+1)} \rightarrow B^{m \times (p+1)}$$ by
817$latex \[
818        y_i^k = F_i^k [ \{ x_j^k \} ]
819\] $$
820We use $latex G : B^{m \times (p+1)} \rightarrow B$$
821to denote an arbitrary scalar valued function of $latex \{ y_i^k \}$$.
822We use $latex H : B^{n \times (p+1)} \rightarrow B$$
823defined by
824$latex \[
825        H ( \{ x_j^k \} ) = G[ F( \{ x_j^k \} ) ]
826\] $$
827
828$head py$$
829The argument $icode py$$ has prototype
830$codei%
831        const CppAD::vector<%Base%>& %py%
832%$$
833and $icode%py%.size() == m * (%p%+1)%$$.
834For $latex i = 0 , \ldots , m-1$$, $latex k = 0 , \ldots , p$$,
835$latex \[
836        py[ i * (p + 1 ) + k ] = \partial G / \partial y_i^k
837\] $$
838
839$subhead px$$
840The $icode px$$ has prototype
841$codei%
842        CppAD::vector<%Base%>& %px%
843%$$
844and $icode%px%.size() == n * (%p%+1)%$$.
845The input values of the elements of $icode px$$
846are not specified (must not matter).
847Upon return,
848for $latex j = 0 , \ldots , n-1$$ and $latex \ell = 0 , \ldots , p$$,
849$latex \[
850\begin{array}{rcl}
851px [ j * (p + 1) + \ell ] & = & \partial H / \partial x_j^\ell
852\\
853& = &
854( \partial G / \partial \{ y_i^k \} )
855        ( \partial \{ y_i^k \} / \partial x_j^\ell )
856\\
857& = &
858\sum_{i=0}^{m-1} \sum_{k=0}^p
859( \partial G / \partial y_i^k ) ( \partial y_i^k / \partial x_j^\ell )
860\\
861& = &
862\sum_{i=0}^{m-1} \sum_{k=\ell}^p
863py[ i * (p + 1 ) + k ] ( \partial F_i^k / \partial x_j^\ell )
864\end{array}
865\] $$
866Note that we have used the fact that for $latex k < \ell$$,
867$latex \partial F_i^k / \partial x_j^\ell = 0$$.
868
869$head ok$$
870The return value $icode ok$$ has prototype
871$codei%
872        bool %ok%
873%$$
874If it is $code true$$, the corresponding evaluation succeeded,
875otherwise it failed.
876
877$end
878-----------------------------------------------------------------------------
879*/
880/*!
881Link from reverse mode sweep to users routine.
882
883\param p [in]
884highest order for this reverse mode calculation.
885
886\param tx [in]
887Taylor coefficients corresponding to \c x for this calculation.
888
889\param ty [in]
890Taylor coefficient corresponding to \c y for this calculation
891
892\param px [out]
893Partials w.r.t. the \c x Taylor coefficients.
894
895\param py [in]
896Partials w.r.t. the \c y Taylor coefficients.
897
898See atomic_reverse mode use documentation
899*/
900virtual bool reverse(
901        size_t                    p  ,
902        const vector<Base>&       tx ,
903        const vector<Base>&       ty ,
904              vector<Base>&       px ,
905        const vector<Base>&       py )
906{       return false; }
907/*
908-------------------------------------- ---------------------------------------
909$begin atomic_for_sparse_jac$$
910$spell
911        afun
912        Jacobian
913        jac
914        const
915        CppAD
916        std
917        bool
918        std
919$$
920
921$section Atomic Forward Jacobian Sparsity Patterns$$
922$index atomic, for_sparse_jac callback$$
923$index for_sparse_jac, atomic callback$$
924$index for_sparse_jac, atomic virtual$$
925
926$head Syntax$$
927$icode%ok% = %afun%.for_sparse_jac(%q%, %r%, %s%)%$$
928
929$head Purpose$$
930This function is used by $cref ForSparseJac$$ to compute
931Jacobian sparsity patterns.
932For a fixed matrix $latex R \in B^{n \times q}$$,
933the Jacobian of $latex f( x + R * u)$$ with respect to $latex u \in B^q$$ is
934$latex \[
935        S(x) = f^{(1)} (x) * R
936\] $$
937Given a $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for $latex R$$,
938$code for_sparse_jac$$ computes a sparsity pattern for $latex S(x)$$.
939
940$head Implementation$$
941If you are using $cref ForSparseJac$$,
942this virtual function must be defined by the
943$cref/atomic_user/atomic_ctor/atomic_user/$$ class.
944
945$subhead q$$
946The argument $icode q$$ has prototype
947$codei%
948     size_t %q%
949%$$
950It specifies the number of columns in
951$latex R \in B^{n \times q}$$ and the Jacobian
952$latex S(x) \in B^{m \times q}$$.
953
954$subhead r$$
955This argument has prototype
956$codei%
957     const %atomic_sparsity%& %r%
958%$$
959and is a $cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
960$latex R \in B^{n \times q}$$.
961
962$subhead s$$
963This argument has prototype
964$codei%
965        %atomic_sparsity%& %s%
966%$$
967The input values of its elements
968are not specified (must not matter).
969Upon return, $icode s$$ is a
970$cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
971$latex S(x) \in B^{m \times q}$$.
972
973$head ok$$
974The return value $icode ok$$ has prototype
975$codei%
976        bool %ok%
977%$$
978If it is $code true$$, the corresponding evaluation succeeded,
979otherwise it failed.
980
981$end
982-----------------------------------------------------------------------------
983*/
984/*!
985Link from forward Jacobian sparsity sweep to atomic_base
986
987\param q
988is the column dimension for the Jacobian sparsity partterns.
989
990\param r
991is the Jacobian sparsity pattern for the argument vector x
992
993\param s
994is the Jacobian sparsity pattern for the result vector y
995*/
996virtual bool for_sparse_jac(
997        size_t                                  q  ,
998        const vector< std::set<size_t> >&       r  ,
999              vector< std::set<size_t> >&       s  )
1000{       return false; }
1001virtual bool for_sparse_jac(
1002        size_t                                  q  ,
1003        const vector<bool>&                     r  ,
1004              vector<bool>&                     s  )
1005{       return false; }
1006/*
1007-------------------------------------- ---------------------------------------
1008$begin atomic_rev_sparse_jac$$
1009$spell
1010        rt
1011        afun
1012        Jacobian
1013        jac
1014        CppAD
1015        std
1016        bool
1017        const
1018        hes
1019$$
1020
1021$section Atomic Reverse Jacobian Sparsity Patterns$$
1022$index atomic, rev_sparse_jac callback$$
1023$index rev_sparse_jac, atomic callback$$
1024$index rev_sparse_jac, atomic virtual$$
1025
1026$head Syntax$$
1027$icode%ok% = %afun%.rev_sparse_jac(%q%, %rt%, %st%)%$$
1028
1029$head Purpose$$
1030This function is used by $cref RevSparseJac$$ to compute
1031Jacobian sparsity patterns.
1032For a fixed matrix $latex R \in B^{q \times m}$$,
1033the Jacobian of $latex R * f( x )$$ with respect to $latex x \in B^q$$ is
1034$latex \[
1035        S(x) = R * f^{(1)} (x)
1036\] $$
1037Given a $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for $latex R$$,
1038$code rev_sparse_jac$$ computes a sparsity pattern for $latex S(x)$$.
1039
1040$head Implementation$$
1041If you are using $cref RevSparseJac$$,
1042this virtual function must be defined by the
1043$cref/atomic_user/atomic_ctor/atomic_user/$$ class.
1044
1045$subhead q$$
1046The argument $icode q$$ has prototype
1047$codei%
1048     size_t %q%
1049%$$
1050It specifies the number of rows in
1051$latex R \in B^{q \times m}$$ and the Jacobian
1052$latex S(x) \in B^{q \times n}$$.
1053
1054$subhead rt$$
1055This argument has prototype
1056$codei%
1057     const %atomic_sparsity%& %rt%
1058%$$
1059and is a
1060$cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
1061$latex R^\R{T} \in B^{m \times q}$$.
1062
1063$subhead st$$
1064This argument has prototype
1065$codei%
1066        %atomic_sparsity%& %st%
1067%$$
1068The input value of its elements
1069are not specified (must not matter).
1070Upon return, $icode s$$ is a
1071$cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
1072$latex S(x)^\R{T} \in B^{n \times q}$$.
1073
1074$head ok$$
1075The return value $icode ok$$ has prototype
1076$codei%
1077        bool %ok%
1078%$$
1079If it is $code true$$, the corresponding evaluation succeeded,
1080otherwise it failed.
1081
1082$end
1083-----------------------------------------------------------------------------
1084*/
1085/*!
1086Link from reverse Jacobian sparsity sweep to atomic_base
1087
1088\param q [in]
1089is the row dimension for the Jacobian sparsity partterns
1090
1091\param rt [out]
1092is the tansposed Jacobian sparsity pattern w.r.t to range variables y
1093
1094\param st [in]
1095is the tansposed Jacobian sparsity pattern for the argument variables x
1096*/
1097virtual bool rev_sparse_jac(
1098        size_t                                  q  ,
1099        const vector< std::set<size_t> >&       rt ,
1100              vector< std::set<size_t> >&       st )
1101{       return false; }
1102virtual bool rev_sparse_jac(
1103        size_t                                  q  ,
1104        const vector<bool>&                     rt ,
1105              vector<bool>&                     st )
1106{       return false; }
1107/*
1108-------------------------------------- ---------------------------------------
1109$begin atomic_rev_sparse_hes$$
1110$spell
1111        vx
1112        afun
1113        Jacobian
1114        jac
1115        CppAD
1116        std
1117        bool
1118        hes
1119        const
1120$$
1121
1122$section Atomic Reverse Hessian Sparsity Patterns$$
1123$index atomic, rev_sparse_hes callback$$
1124$index rev_sparse_hes, atomic callback$$
1125$index rev_sparse_hes, atomic virtual$$
1126
1127$head Syntax$$
1128$icode%ok% = %afun%.rev_sparse_hes(%vx%, %s%, %t%, %q%, %r%, %u%, %v%)%$$
1129
1130$head Purpose$$
1131This function is used by $cref RevSparseHes$$ to compute
1132Hessian sparsity patterns.
1133There is an unspecified scalar valued function
1134$latex g : B^m \rightarrow B$$.
1135Given a $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for
1136$latex R \in B^{n \times q}$$,
1137and information about the function $latex z = g(y)$$,
1138this routine computes the sparsity pattern for
1139$latex \[
1140        V(x) = (g \circ f)^{(2)}( x ) R
1141\] $$
1142
1143$head Implementation$$
1144If you are using and $cref RevSparseHes$$,
1145this virtual function must be defined by the
1146$cref/atomic_user/atomic_ctor/atomic_user/$$ class.
1147
1148$subhead vx$$
1149The argument $icode vx$$ has prototype
1150$codei%
1151     const CppAD:vector<bool>& %vx%
1152%$$
1153$icode%vx%.size() == %n%$$, and
1154for $latex j = 0 , \ldots , n-1$$,
1155$icode%vx%[%j%]%$$ is true if and only if
1156$icode%ax%[%j%]%$$ is a $cref/variable/glossary/Variable/$$
1157in the corresponding call to
1158$codei%
1159        %afun%(%ax%, %ay%, %id%)
1160%$$
1161
1162$subhead s$$
1163The argument $icode s$$ has prototype
1164$codei%
1165     const CppAD:vector<bool>& %s%
1166%$$
1167and its size is $icode m$$.
1168It is a sparsity pattern for
1169$latex S(x) = g^{(1)} (y) \in B^{1 \times m}$$.
1170
1171$subhead t$$
1172This argument has prototype
1173$codei%
1174     CppAD:vector<bool>& %t%
1175%$$
1176and its size is $icode m$$.
1177The input values of its elements
1178are not specified (must not matter).
1179Upon return, $icode t$$ is a
1180sparsity pattern for
1181$latex T(x) \in B^{1 \times n}$$ where
1182$latex \[
1183        T(x) = (g \circ f)^{(1)} (x) = S(x) * f^{(1)} (x)
1184\]$$
1185
1186$subhead q$$
1187The argument $icode q$$ has prototype
1188$codei%
1189     size_t %q%
1190%$$
1191It specifies the number of columns in
1192$latex R \in B^{n \times q}$$,
1193$latex U(x) \in B^{m \times q}$$, and
1194$latex V(x) \in B^{n \times q}$$.
1195
1196$subhead r$$
1197This argument has prototype
1198$codei%
1199     const %atomic_sparsity%& %r%
1200%$$
1201and is a $cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
1202$latex R \in B^{n \times q}$$.
1203
1204$head u$$
1205This argument has prototype
1206$codei%
1207     const %atomic_sparsity%& %u%
1208%$$
1209and is a $cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
1210$latex U(x) \in B^{m \times q}$$ which is defined by
1211$latex \[
1212\begin{array}{rcl}
1213U(x)
1214& = &
1215\partial_u \{ \partial_y g[ y + f^{(1)} (x) R u ] \}_{u=0}
1216\\
1217& = &
1218\partial_u \{ g^{(1)} [ y + f^{(1)} (x) R u ] \}_{u=0}
1219\\
1220& = &
1221g^{(2)} (y) f^{(1)} (x) R
1222\end{array}
1223\] $$
1224
1225$subhead v$$
1226This argument has prototype
1227$codei%
1228     %atomic_sparsity%& %v%
1229%$$
1230The input value of its elements
1231are not specified (must not matter).
1232Upon return, $icode v$$ is a
1233$cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
1234$latex V(x) \in B^{n \times q}$$ which is defined by
1235$latex \[
1236\begin{array}{rcl}
1237V(x)
1238& = &
1239\partial_u [ \partial_x (g \circ f) ( x + R u )  ]_{u=0}
1240\\
1241& = &
1242\partial_u [ (g \circ f)^{(1)}( x + R u )  ]_{u=0}
1243\\
1244& = &
1245(g \circ f)^{(2)}( x ) R
1246\\
1247& = &
1248f^{(1)} (x)^\R{T} g^{(2)} ( y ) f^{(1)} (x)  R
1249+
1250\sum_{i=1}^m g_i^{(1)} (y) \; f_i^{(2)} (x) R
1251\\
1252& = &
1253f^{(1)} (x)^\R{T} U(x)
1254+
1255\sum_{i=1}^m S_i (x) \; f_i^{(2)} (x) R
1256\end{array}
1257\] $$
1258
1259$end
1260-----------------------------------------------------------------------------
1261*/
1262/*!
1263Link from reverse Hessian sparsity sweep to base_atomic
1264
1265\param vx [in]
1266which componens of x are variables.
1267
1268\param s [in]
1269is the reverse Jacobian sparsity pattern w.r.t the result vector y.
1270
1271\param t [out]
1272is the reverse Jacobian sparsity pattern w.r.t the argument vector x.
1273
1274\param q [in]
1275is the column dimension for the sparsity partterns.
1276
1277\param r [in]
1278is the forward Jacobian sparsity pattern w.r.t the argument vector x
1279
1280\param u [in]
1281is the Hessian sparsity pattern w.r.t the result vector y.
1282
1283\param v [out]
1284is the Hessian sparsity pattern w.r.t the argument vector x.
1285*/
1286virtual bool rev_sparse_hes(
1287        const vector<bool>&                     vx ,
1288        const vector<bool>&                     s  ,
1289              vector<bool>&                     t  ,
1290        size_t                                  q  ,
1291        const vector< std::set<size_t> >&       r  ,
1292        const vector< std::set<size_t> >&       u  ,
1293              vector< std::set<size_t> >&       v  )
1294{       return false; }
1295virtual bool rev_sparse_hes(
1296        const vector<bool>&                     vx ,
1297        const vector<bool>&                     s  ,
1298              vector<bool>&                     t  ,
1299        size_t                                  q  ,
1300        const vector<bool>&                     r  ,
1301        const vector<bool>&                     u  ,
1302              vector<bool>&                     v  )
1303{       return false; }
1304/*
1305------------------------------------------------------------------------------
1306$begin atomic_base_clear$$
1307
1308$section Free Static Variables$$
1309$index free, atomic static$$
1310$index atomic, free static$$
1311$index static, free atomic$$
1312$index clear, atomic static$$
1313
1314$head Syntax$$
1315$codei%atomic_base<%Base%>::clear()%$$
1316
1317$head Purpose$$
1318The $code atomic_base$$ class holds onto static work space in order to
1319increase speed by avoiding system memory allocation calls.
1320This call makes to work space $cref/available/ta_available/$$ to
1321for other uses by the same thread.
1322This should be called when you are done using the
1323user atomic functions for a specific value of $icode Base$$.
1324
1325$head Restriction$$
1326This routine cannot be called
1327while in $cref/parallel/ta_in_parallel/$$ execution mode.
1328
1329$end
1330------------------------------------------------------------------------------
1331*/
1332/*!
1333Free all thread_alloc static memory held by atomic_base (avoids reallocations).
1334(This does not include class_object() which is an std::vector.)
1335*/
1336/// Free vector memory used by this class (work space)
1337static void clear(void)
1338{       CPPAD_ASSERT_KNOWN(
1339                ! thread_alloc::in_parallel() ,
1340                "cannot use atomic_base clear during parallel execution"
1341        );
1342        size_t i = class_object().size();
1343        while(i--)
1344        {       size_t thread = CPPAD_MAX_NUM_THREADS;
1345                while(thread--)
1346                {
1347                        atomic_base* op = class_object()[i];
1348                        if( op != CPPAD_NULL )
1349                        {       op->afun_vx_[thread].clear();
1350                                op->afun_vy_[thread].clear();
1351                                op->afun_tx_[thread].clear();
1352                                op->afun_ty_[thread].clear();
1353                        }
1354                }
1355        }
1356        return;
1357}
1358// -------------------------------------------------------------------------
1359/*!
1360Set value of id (used by deprecated old_atomic class)
1361
1362This function is called just before calling any of the virtual funcitons
1363and has the corresponding id of the corresponding virtual call.
1364*/
1365virtual void set_id(size_t id) 
1366{ }
1367// ---------------------------------------------------------------------------
1368};
1369/*! \} */
1370CPPAD_END_NAMESPACE
1371# endif
Note: See TracBrowser for help on using the repository browser.