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

Change comment so end of group is included in doxygen processing
(do not know why works without this command at end).

• Property svn:keywords set to Id
File size: 17.3 KB
Line
1/* $Id: reverse_sweep.hpp 2625 2012-12-23 14:34:12Z bradbell$ */
4
5/* --------------------------------------------------------------------------
7
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.
14-------------------------------------------------------------------------- */
15
16
18/*!
19\defgroup reverse_sweep_hpp reverse_sweep.hpp
20\{
21\file reverse_sweep.hpp
22Compute derivatives of arbitrary order Taylor coefficients.
23*/
24
25/*!
27This value is either zero or one.
28Zero is the normal operational value.
29If it is one, a trace of every reverse_sweep computation is printed.
30*/
32
33/*!
34Compute derivative of arbitrary order forward mode Taylor coefficients.
35
36\tparam Base
37base type for the operator; i.e., this operation sequence was recorded
38using AD< \a Base > and computations by this routine are done using type
39\a Base.
40
41\param d
42is the highest order Taylor coefficients that
43we are computing the derivative of.
44
45\param n
46is the number of independent variables on the tape.
47
48\param numvar
49is the total number of variables on the tape.
50This is also equal to the number of rows in the matrix \a Taylor; i.e.,
51Rec->num_rec_var().
52
53\param Rec
54The information stored in \a Rec
55is a recording of the operations corresponding to the function
56\f[
57        F : {\bf R}^n \rightarrow {\bf R}^m
58\f]
59where \f$n \f$ is the number of independent variables and
60\f$m \f$ is the number of dependent variables.
61We define the function
62\f$G : {\bf R}^{n \times d} \rightarrow {\bf R} \f$ by
63\f[
64G( u ) = \frac{1}{d !} \frac{ \partial^d }{ \partial t^d }
65\left[
66        \sum_{i=1}^m w_i  F_i ( u^{(0)} + u^{(1)} t + \cdots + u^{(d)} t^d )
67\right]_{t=0}
68\f]
69Note that the scale factor  1 / a d  converts
70the \a d-th partial derivative to the \a d-th order Taylor coefficient.
71This routine computes the derivative of \f$G(u) \f$
72with respect to all the Taylor coefficients
73\f$u^{(k)} \f$ for \f$k = 0 , ... , d \f$.
74The vector \f$w \in {\bf R}^m \f$, and
75value of \f$u \in {\bf R}^{n \times d} \f$
76at which the derivative is computed,
77are defined below.
78\n
79\n
80The object \a Rec is effectly constant.
81There is an exception to this,
82while palying back the tape
83the object \a Rec holds information about the current location
84with in the tape and this changes during palyback.
85
86\param J
87Is the number of columns in the coefficient matrix \a Taylor.
88This must be greater than or equal \a d + 1.
89
90\param Taylor
91For i = 1 , ... , \a numvar, and for k = 0 , ... , \a d,
92\a Taylor [ i * J + k ]
93is the k-th order Taylor coefficient corresponding to
94variable with index i on the tape.
95The value \f$u \in {\bf R}^{n \times d} \f$,
96at which the derivative is computed,
97is defined by
98\f$u_j^{(k)} \f$ = \a Taylor [ j * J + k ]
99for j = 1 , ... , \a n, and for k = 0 , ... , \a d.
100
101\param K
102Is the number of columns in the partial derivative matrix \a Partial.
103It must be greater than or equal \a d + 1.
104
105\param Partial
106\b Input:
107The last \f$m \f$ rows of \a Partial are inputs.
108The vector \f$v \f$, used to define \f$G(u) \f$,
109is specified by these rows.
110For i = 0 , ... , m - 1, \a Partial [ ( \a numvar - m + i ) * K + d ] = v_i.
111For i = 0 , ... , m - 1 and for k = 0 , ... , d - 1,
112\a Partial [ ( \a numvar - m + i ) * K + k ] = 0.
113\n
114\n
115\b Temporary:
116For i = n+1 , ... , \a numvar - 1 and for k = 0 , ... , d,
117the value of \a Partial [ i * K + k ] is used for temporary work space
118and its output value is not defined.
119\n
120\n
121\b Output:
122For j = 1 , ... , n and for k = 0 , ... , d,
123\a Partial [ j * K + k ]
124is the partial derivative of \f$G( u ) \f$ with
125respect to \f$u_j^{(k)} \f$.
126
127\par Assumptions
128The first operator on the tape is a BeginOp,
129and the next \a n operators are InvOp operations for the
130corresponding independent variables.
131*/
132template <class Base>
133void ReverseSweep(
134        size_t                d,
135        size_t                n,
136        size_t                numvar,
137        player<Base>*         Rec,
138        size_t                J,
139        const Base*           Taylor,
140        size_t                K,
141        Base*                 Partial
142)
143{
144        OpCode           op;
145        size_t         i_op;
146        size_t        i_var;
147
148        const addr_t*   arg = 0;
149
150        // check numvar argument
151        CPPAD_ASSERT_UNKNOWN( Rec->num_rec_var() == numvar );
152        CPPAD_ASSERT_UNKNOWN( numvar > 0 );
153
154        // length of the parameter vector (used by CppAD assert macros)
155        const size_t num_par = Rec->num_rec_par();
156
157        // pointer to the beginning of the parameter vector
158        const Base* parameter = 0;
159        if( num_par > 0 )
160                parameter = Rec->GetPar();
161
162        // work space used by UserOp.
163        const size_t user_k  = d;    // order of this forward mode calculation
164        const size_t user_k1 = d+1;  // number of orders for this calculation
165        vector<size_t> user_ix;      // variable indices for argument vector
166        vector<Base> user_tx;        // argument vector Taylor coefficients
167        vector<Base> user_ty;        // result vector Taylor coefficients
168        vector<Base> user_px;        // partials w.r.t argument vector
169        vector<Base> user_py;        // partials w.r.t. result vector
170        size_t user_index = 0;       // indentifier for this user_atomic operation
171        size_t user_id    = 0;       // user identifier for this call to operator
172        size_t user_i     = 0;       // index in result vector
173        size_t user_j     = 0;       // index in argument vector
174        size_t user_m     = 0;       // size of result vector
175        size_t user_n     = 0;       // size of arugment vector
176        // next expected operator in a UserOp sequence
177        enum { user_start, user_arg, user_ret, user_end } user_state = user_end;
178
179        // temporary indices
180        size_t j, ell;
181
182        // Initialize
183        Rec->start_reverse(op, arg, i_op, i_var);
184        CPPAD_ASSERT_UNKNOWN( op == EndOp );
186        std::cout << std::endl;
187# endif
188        while(op != BeginOp )
189        {       // next op
190                Rec->next_reverse(op, arg, i_op, i_var);
191# ifndef NDEBUG
192                if( i_op <= n )
193                {       CPPAD_ASSERT_UNKNOWN((op == InvOp) | (op == BeginOp));
194                }
195                else    CPPAD_ASSERT_UNKNOWN((op != InvOp) & (op != BeginOp));
196# endif
197
198                // rest of informaiton depends on the case
200                size_t       i_tmp  = i_var;
201                const Base*  Z_tmp  = Taylor + i_var * J;
202                const Base*  pZ_tmp = Partial + i_var * K;
203
204                printOp(
205                        std::cout,
206                        Rec,
207                        i_tmp,
208                        op,
209                        arg,
210                        d + 1,
211                        Z_tmp,
212                        d + 1,
213                        pZ_tmp
214                );
215# endif
216
217                switch( op )
218                {
219
220                        case AbsOp:
221                        reverse_abs_op(
222                                d, i_var, arg[0], J, Taylor, K, Partial
223                        );
224                        break;
225                        // --------------------------------------------------
226
229                                d, i_var, arg, parameter, J, Taylor, K, Partial
230                        );
231                        break;
232                        // --------------------------------------------------
233
235                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
237                                d, i_var, arg, parameter, J, Taylor, K, Partial
238                        );
239                        break;
240                        // --------------------------------------------------
241
242                        case AcosOp:
243                        // sqrt(1 - x * x), acos(x)
244                        CPPAD_ASSERT_UNKNOWN( i_var < numvar );
245                        reverse_acos_op(
246                                d, i_var, arg[0], J, Taylor, K, Partial
247                        );
248                        break;
249                        // --------------------------------------------------
250
251                        case AsinOp:
252                        // sqrt(1 - x * x), asin(x)
253                        CPPAD_ASSERT_UNKNOWN( i_var < numvar );
254                        reverse_asin_op(
255                                d, i_var, arg[0], J, Taylor, K, Partial
256                        );
257                        break;
258                        // --------------------------------------------------
259
260                        case AtanOp:
261                        // 1 + x * x, atan(x)
262                        CPPAD_ASSERT_UNKNOWN( i_var < numvar );
263                        reverse_atan_op(
264                                d, i_var, arg[0], J, Taylor, K, Partial
265                        );
266                        break;
267                        // -------------------------------------------------
268
269                        case BeginOp:
271                        break;
272                        // --------------------------------------------------
273
274                        case CSumOp:
275                        // CSumOp has a variable number of arguments and
276                        // next_reverse thinks it one has one argument.
277                        // We must inform next_reverse of this special case.
278                        Rec->reverse_csum(op, arg, i_op, i_var);
279                        reverse_csum_op(
280                                d, i_var, arg, K, Partial
281                        );
282                        // end of a cummulative summation
283                        break;
284                        // -------------------------------------------------
285
286                        case CExpOp:
287                        reverse_cond_op(
288                                d,
289                                i_var,
290                                arg,
291                                num_par,
292                                parameter,
293                                J,
294                                Taylor,
295                                K,
296                                Partial
297                        );
298                        break;
299                        // --------------------------------------------------
300
301                        case ComOp:
302                        break;
303                        // --------------------------------------------------
304
305                        case CosOp:
306                        CPPAD_ASSERT_UNKNOWN( i_var < numvar );
307                        reverse_cos_op(
308                                d, i_var, arg[0], J, Taylor, K, Partial
309                        );
310                        break;
311                        // --------------------------------------------------
312
313                        case CoshOp:
314                        CPPAD_ASSERT_UNKNOWN( i_var < numvar );
315                        reverse_cosh_op(
316                                d, i_var, arg[0], J, Taylor, K, Partial
317                        );
318                        break;
319                        // --------------------------------------------------
320
321                        case DisOp:
322                        // Derivative of discrete operation is zero so no
323                        // contribution passes through this operation.
324                        break;
325                        // --------------------------------------------------
326
327                        case DivvvOp:
328                        reverse_divvv_op(
329                                d, i_var, arg, parameter, J, Taylor, K, Partial
330                        );
331                        break;
332                        // --------------------------------------------------
333
334                        case DivpvOp:
335                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
336                        reverse_divpv_op(
337                                d, i_var, arg, parameter, J, Taylor, K, Partial
338                        );
339                        break;
340                        // --------------------------------------------------
341
342                        case DivvpOp:
343                        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < num_par );
344                        reverse_divvp_op(
345                                d, i_var, arg, parameter, J, Taylor, K, Partial
346                        );
347                        break;
348                        // --------------------------------------------------
349
350                        case ExpOp:
351                        reverse_exp_op(
352                                d, i_var, arg[0], J, Taylor, K, Partial
353                        );
354                        break;
355                        // --------------------------------------------------
356                        case LdpOp:
358                                op, d, i_var, arg, J, Taylor, K, Partial
359                        );
360                        break;
361                        // -------------------------------------------------
362
363                        case LdvOp:
365                                op, d, i_var, arg, J, Taylor, K, Partial
366                        );
367                        break;
368                        // -------------------------------------------------
369
370                        case InvOp:
371                        break;
372                        // --------------------------------------------------
373
374                        case LogOp:
375                        reverse_log_op(
376                                d, i_var, arg[0], J, Taylor, K, Partial
377                        );
378                        break;
379                        // --------------------------------------------------
380
381                        case MulvvOp:
382                        reverse_mulvv_op(
383                                d, i_var, arg, parameter, J, Taylor, K, Partial
384                        );
385                        break;
386                        // --------------------------------------------------
387
388                        case MulpvOp:
389                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
390                        reverse_mulpv_op(
391                                d, i_var, arg, parameter, J, Taylor, K, Partial
392                        );
393                        break;
394                        // --------------------------------------------------
395
396                        case ParOp:
397                        break;
398                        // --------------------------------------------------
399
400                        case PowvpOp:
401                        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < num_par );
402                        reverse_powvp_op(
403                                d, i_var, arg, parameter, J, Taylor, K, Partial
404                        );
405                        break;
406                        // -------------------------------------------------
407
408                        case PowpvOp:
409                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
410                        reverse_powpv_op(
411                                d, i_var, arg, parameter, J, Taylor, K, Partial
412                        );
413                        break;
414                        // -------------------------------------------------
415
416                        case PowvvOp:
417                        reverse_powvv_op(
418                                d, i_var, arg, parameter, J, Taylor, K, Partial
419                        );
420                        break;
421                        // --------------------------------------------------
422
423                        case PriOp:
424                        // no result so nothing to do
425                        break;
426                        // --------------------------------------------------
427
428                        case SignOp:
429                        CPPAD_ASSERT_UNKNOWN( i_var < numvar );
430                        reverse_sign_op(
431                                d, i_var, arg[0], J, Taylor, K, Partial
432                        );
433                        break;
434                        // -------------------------------------------------
435
436                        case SinOp:
437                        CPPAD_ASSERT_UNKNOWN( i_var < numvar );
438                        reverse_sin_op(
439                                d, i_var, arg[0], J, Taylor, K, Partial
440                        );
441                        break;
442                        // -------------------------------------------------
443
444                        case SinhOp:
445                        CPPAD_ASSERT_UNKNOWN( i_var < numvar );
446                        reverse_sinh_op(
447                                d, i_var, arg[0], J, Taylor, K, Partial
448                        );
449                        break;
450                        // --------------------------------------------------
451
452                        case SqrtOp:
453                        reverse_sqrt_op(
454                                d, i_var, arg[0], J, Taylor, K, Partial
455                        );
456                        break;
457                        // --------------------------------------------------
458
459                        case StppOp:
460                        break;
461                        // --------------------------------------------------
462
463                        case StpvOp:
464                        break;
465                        // -------------------------------------------------
466
467                        case StvpOp:
468                        break;
469                        // -------------------------------------------------
470
471                        case StvvOp:
472                        break;
473                        // --------------------------------------------------
474
475                        case SubvvOp:
476                        reverse_subvv_op(
477                                d, i_var, arg, parameter, J, Taylor, K, Partial
478                        );
479                        break;
480                        // --------------------------------------------------
481
482                        case SubpvOp:
483                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
484                        reverse_subpv_op(
485                                d, i_var, arg, parameter, J, Taylor, K, Partial
486                        );
487                        break;
488                        // --------------------------------------------------
489
490                        case SubvpOp:
491                        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < num_par );
492                        reverse_subvp_op(
493                                d, i_var, arg, parameter, J, Taylor, K, Partial
494                        );
495                        break;
496                        // -------------------------------------------------
497
498                        case TanOp:
499                        CPPAD_ASSERT_UNKNOWN( i_var < numvar );
500                        reverse_tan_op(
501                                d, i_var, arg[0], J, Taylor, K, Partial
502                        );
503                        break;
504                        // -------------------------------------------------
505
506                        case TanhOp:
507                        CPPAD_ASSERT_UNKNOWN( i_var < numvar );
508                        reverse_tanh_op(
509                                d, i_var, arg[0], J, Taylor, K, Partial
510                        );
511                        break;
512                        // --------------------------------------------------
513
514                        case UserOp:
515                        // start or end an atomic operation sequence
516                        CPPAD_ASSERT_UNKNOWN( NumRes( UserOp ) == 0 );
517                        CPPAD_ASSERT_UNKNOWN( NumArg( UserOp ) == 4 );
518                        if( user_state == user_end )
519                        {       user_index = arg[0];
520                                user_id    = arg[1];
521                                user_n     = arg[2];
522                                user_m     = arg[3];
523                                if(user_ix.size() < user_n)
524                                        user_ix.resize(user_n);
525                                if(user_tx.size() < user_n * user_k1)
526                                {       user_tx.resize(user_n * user_k1);
527                                        user_px.resize(user_n * user_k1);
528                                }
529                                if(user_ty.size() < user_m * user_k1)
530                                {       user_ty.resize(user_m * user_k1);
531                                        user_py.resize(user_m * user_k1);
532                                }
533                                user_j     = user_n;
534                                user_i     = user_m;
535                                user_state = user_ret;
536                        }
537                        else
538                        {       CPPAD_ASSERT_UNKNOWN( user_state == user_start );
539                                CPPAD_ASSERT_UNKNOWN( user_index == size_t(arg[0]) );
540                                CPPAD_ASSERT_UNKNOWN( user_id    == size_t(arg[1]) );
541                                CPPAD_ASSERT_UNKNOWN( user_n     == size_t(arg[2]) );
542                                CPPAD_ASSERT_UNKNOWN( user_m     == size_t(arg[3]) );
543                                user_state = user_end;
544
545                                // call users function for this operation
546                                user_atomic<Base>::reverse(user_index, user_id,
547                                        user_k, user_n, user_m, user_tx, user_ty,
548                                        user_px, user_py
549                                );
550                                for(j = 0; j < user_n; j++) if( user_ix[j] > 0 )
551                                {       for(ell = 0; ell < user_k1; ell++)
552                                                Partial[user_ix[j] * K + ell] +=
553                                                        user_px[j * user_k1 + ell];
554                                }
555                        }
556                        break;
557
558                        case UsrapOp:
559                        // parameter argument in an atomic operation sequence
560                        CPPAD_ASSERT_UNKNOWN( user_state == user_arg );
561                        CPPAD_ASSERT_UNKNOWN( 0 < user_j && user_j <= user_n );
562                        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
563                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
564                        --user_j;
565                        user_ix[user_j] = 0;
566                        user_tx[user_j * user_k1 + 0] = parameter[ arg[0]];
567                        for(ell = 1; ell < user_k1; ell++)
568                                user_tx[user_j * user_k1 + ell] = Base(0.);
569
570                        if( user_j == 0 )
571                                user_state = user_start;
572                        break;
573
574                        case UsravOp:
575                        // variable argument in an atomic operation sequence
576                        CPPAD_ASSERT_UNKNOWN( user_state == user_arg );
577                        CPPAD_ASSERT_UNKNOWN( 0 < user_j && user_j <= user_n );
578                        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
579                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) <= i_var );
580                        CPPAD_ASSERT_UNKNOWN( 0 < arg[0] );
581                        --user_j;
582                        user_ix[user_j] = arg[0];
583                        for(ell = 0; ell < user_k1; ell++)
584                                user_tx[user_j*user_k1 + ell] = Taylor[ arg[0] * J + ell];
585                        if( user_j == 0 )
586                                user_state = user_start;
587                        break;
588
589                        case UsrrpOp:
590                        // parameter result in an atomic operation sequence
591                        CPPAD_ASSERT_UNKNOWN( user_state == user_ret );
592                        CPPAD_ASSERT_UNKNOWN( 0 < user_i && user_i <= user_m );
593                        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
594                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
595                        --user_i;
596                        for(ell = 0; ell < user_k1; ell++)
597                        {       user_py[user_i * user_k1 + ell] = Base(0.);
598                                user_ty[user_i * user_k1 + ell] = Base(0.);
599                        }
600                        user_ty[user_i * user_k1 + 0] = parameter[ arg[0] ];
601                        if( user_i == 0 )
602                                user_state = user_arg;
603                        break;
604
605                        case UsrrvOp:
606                        // variable result in an atomic operation sequence
607                        CPPAD_ASSERT_UNKNOWN( user_state == user_ret );
608                        CPPAD_ASSERT_UNKNOWN( 0 < user_i && user_i <= user_m );
609                        --user_i;
610                        for(ell = 0; ell < user_k1; ell++)
611                        {       user_py[user_i * user_k1 + ell] =
612                                                Partial[i_var * K + ell];
613                                user_ty[user_i * user_k1 + ell] =
614                                                Taylor[i_var * J + ell];
615                        }
616                        if( user_i == 0 )
617                                user_state = user_arg;
618                        break;
619                        // ------------------------------------------------------------
620
621                        default:
623                }
624        }
626        std::cout << std::endl;
627# endif
628        // values corresponding to BeginOp
629        CPPAD_ASSERT_UNKNOWN( i_op == 0 );
630        CPPAD_ASSERT_UNKNOWN( i_var == 0 );
631}
632
633/*! \} */