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

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

Change atomic_base::list to atomic_base::class_object to reflect that it
is a list of all the objects in that class.

vector.hpp: Fix documentation error (see whats_new_13 changes).

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