CppAD: A C++ Algorithmic Differentiation Package 20110419
reverse_sweep.hpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 # ifndef CPPAD_REVERSE_SWEEP_INCLUDED
00003 # define CPPAD_REVERSE_SWEEP_INCLUDED
00004 
00005 /* --------------------------------------------------------------------------
00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-11 Bradley M. Bell
00007 
00008 CppAD is distributed under multiple licenses. This distribution is under
00009 the terms of the 
00010                     Common Public License Version 1.0.
00011 
00012 A copy of this license is included in the COPYING file of this distribution.
00013 Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
00014 -------------------------------------------------------------------------- */
00015 
00016 
00017 CPPAD_BEGIN_NAMESPACE
00018 /*!
00019 \file reverse_sweep.hpp
00020 Compute derivatives of arbitrary order Taylor coefficients.
00021 */
00022 
00023 /*!
00024 \def CPPAD_REVERSE_SWEEP_TRACE
00025 This value is either zero or one. 
00026 Zero is the normal operational value.
00027 If it is one, a trace of every reverse_sweep computation is printed.
00028 */
00029 # define CPPAD_REVERSE_SWEEP_TRACE 0
00030 
00031 /*!
00032 Compute derivative of arbitrary order forward mode Taylor coefficients.
00033 
00034 \tparam Base
00035 base type for the operator; i.e., this operation sequence was recorded
00036 using AD< \a Base > and computations by this routine are done using type 
00037 \a Base.
00038 
00039 \param d
00040 is the highest order Taylor coefficients that 
00041 we are computing the derivative of.
00042 
00043 \param n
00044 is the number of independent variables on the tape.
00045 
00046 \param numvar
00047 is the total number of variables on the tape.
00048 This is also equal to the number of rows in the matrix \a Taylor; i.e.,
00049 Rec->num_rec_var().
00050 
00051 \param Rec
00052 The information stored in \a Rec
00053 is a recording of the operations corresponding to the function
00054 \f[
00055         F : {\bf R}^n \rightarrow {\bf R}^m
00056 \f]
00057 where \f$ n \f$ is the number of independent variables and
00058 \f$ m \f$ is the number of dependent variables.
00059 We define the function 
00060 \f$ G : {\bf R}^{n \times d} \rightarrow {\bf R} \f$ by
00061 \f[
00062 G( u ) = \frac{1}{d !} \frac{ \partial^d }{ \partial t^d } 
00063 \left[ 
00064         \sum_{i=1}^m w_i  F_i ( u^{(0)} + u^{(1)} t + \cdots + u^{(d)} t^d )
00065 \right]_{t=0}
00066 \f]
00067 Note that the scale factor  1 / a d  converts 
00068 the \a d-th partial derivative to the \a d-th order Taylor coefficient.
00069 This routine computes the derivative of \f$ G(u) \f$
00070 with respect to all the Taylor coefficients
00071 \f$ u^{(k)} \f$ for \f$ k = 0 , ... , d \f$.
00072 The vector \f$ w \in {\bf R}^m \f$, and
00073 value of \f$ u \in {\bf R}^{n \times d} \f$
00074 at which the derivative is computed,
00075 are defined below.
00076 \n
00077 \n
00078 The object \a Rec is effectly constant.
00079 There is an exception to this,
00080 while palying back the tape
00081 the object \a Rec holds information about the current location
00082 with in the tape and this changes during palyback. 
00083 
00084 \param J
00085 Is the number of columns in the coefficient matrix \a Taylor.
00086 This must be greater than or equal \a d + 1.
00087 
00088 \param Taylor
00089 For i = 1 , ... , \a numvar, and for k = 0 , ... , \a d,
00090 \a Taylor [ i * J + k ]
00091 is the k-th order Taylor coefficient corresponding to 
00092 variable with index i on the tape.
00093 The value \f$ u \in {\bf R}^{n \times d} \f$,
00094 at which the derivative is computed,
00095 is defined by
00096 \f$ u_j^{(k)} \f$ = \a Taylor [ j * J + k ]
00097 for j = 1 , ... , \a n, and for k = 0 , ... , \a d.
00098 
00099 \param K
00100 Is the number of columns in the partial derivative matrix \a Partial.
00101 It must be greater than or equal \a d + 1.
00102 
00103 \param Partial
00104 \b Input:
00105 The last \f$ m \f$ rows of \a Partial are inputs.
00106 The vector \f$ v \f$, used to define \f$ G(u) \f$,
00107 is specified by these rows. 
00108 For i = 0 , ... , m - 1, \a Partial [ ( \a numvar - m + i ) * K + d ] = v_i.
00109 For i = 0 , ... , m - 1 and for k = 0 , ... , d - 1, 
00110 \a Partial [ ( \a numvar - m + i ) * K + k ] = 0.
00111 \n
00112 \n
00113 \b Temporary:
00114 For i = n+1 , ... , \a numvar - 1 and for k = 0 , ... , d, 
00115 the value of \a Partial [ i * K + k ] is used for temporary work space
00116 and its output value is not defined. 
00117 \n
00118 \n
00119 \b Output:
00120 For j = 1 , ... , n and for k = 0 , ... , d, 
00121 \a Partial [ j * K + k ] 
00122 is the partial derivative of \f$ G( u ) \f$ with 
00123 respect to \f$ u_j^{(k)} \f$.
00124 
00125 \par Assumptions
00126 The first operator on the tape is a BeginOp,
00127 and the next \a n operators are InvOp operations for the 
00128 corresponding independent variables.
00129 */
00130 template <class Base>
00131 void ReverseSweep(
00132         size_t                d,
00133         size_t                n,
00134         size_t                numvar,
00135         player<Base>*         Rec,
00136         size_t                J,
00137         const Base*           Taylor,
00138         size_t                K,
00139         Base*                 Partial
00140 )
00141 {
00142         OpCode           op;
00143         size_t         i_op;
00144         size_t        i_var;
00145 
00146         const size_t   *arg = 0;
00147 
00148         // check numvar argument
00149         CPPAD_ASSERT_UNKNOWN( Rec->num_rec_var() == numvar );
00150         CPPAD_ASSERT_UNKNOWN( numvar > 0 );
00151 
00152         // length of the parameter vector (used by CppAD assert macros)
00153         const size_t num_par = Rec->num_rec_par();
00154 
00155         // pointer to the beginning of the parameter vector
00156         const Base* parameter = 0;
00157         if( num_par > 0 )
00158                 parameter = Rec->GetPar();
00159 
00160         // work space used by UserOp.
00161         const size_t user_k  = d;    // order of this forward mode calculation
00162         const size_t user_k1 = d+1;  // number of orders for this calculation
00163         vector<size_t> user_ix;      // variable indices for argument vector
00164         vector<Base> user_tx;        // argument vector Taylor coefficients
00165         vector<Base> user_ty;        // result vector Taylor coefficients
00166         vector<Base> user_px;        // partials w.r.t argument vector
00167         vector<Base> user_py;        // partials w.r.t. result vector
00168         size_t user_index = 0;       // indentifier for this user_atomic operation
00169         size_t user_id    = 0;       // user identifier for this call to operator
00170         size_t user_i     = 0;       // index in result vector
00171         size_t user_j     = 0;       // index in argument vector
00172         size_t user_m     = 0;       // size of result vector
00173         size_t user_n     = 0;       // size of arugment vector
00174         // next expected operator in a UserOp sequence
00175         enum { user_start, user_arg, user_ret, user_end } user_state = user_end;
00176 
00177         // temporary indices
00178         size_t j, ell;
00179 
00180         // Initialize
00181         Rec->start_reverse(op, arg, i_op, i_var);
00182         CPPAD_ASSERT_UNKNOWN( op == EndOp );
00183 # if CPPAD_REVERSE_SWEEP_TRACE
00184         std::cout << std::endl;
00185 # endif
00186         while(op != BeginOp )
00187         {       // next op
00188                 Rec->next_reverse(op, arg, i_op, i_var);
00189 # ifndef NDEBUG
00190                 if( i_op <= n )
00191                 {       CPPAD_ASSERT_UNKNOWN((op == InvOp) | (op == BeginOp));
00192                 }
00193                 else    CPPAD_ASSERT_UNKNOWN((op != InvOp) & (op != BeginOp));
00194 # endif
00195 
00196                 // rest of informaiton depends on the case
00197 # if CPPAD_REVERSE_SWEEP_TRACE
00198                 size_t       i_tmp  = i_var;
00199                 const Base*  Z_tmp  = Taylor + i_var * J;
00200                 const Base*  pZ_tmp = Partial + i_var * K;
00201 
00202                 printOp(
00203                         std::cout, 
00204                         Rec,
00205                         i_tmp,
00206                         op, 
00207                         arg,
00208                         d + 1, 
00209                         Z_tmp, 
00210                         d + 1, 
00211                         pZ_tmp 
00212                 );
00213 # endif
00214 
00215                 switch( op )
00216                 {
00217 
00218                         case AbsOp:
00219                         reverse_abs_op(
00220                                 d, i_var, arg[0], J, Taylor, K, Partial
00221                         );
00222                         break;
00223                         // --------------------------------------------------
00224 
00225                         case AddvvOp:
00226                         reverse_addvv_op(
00227                                 d, i_var, arg, parameter, J, Taylor, K, Partial
00228                         );
00229                         break;
00230                         // --------------------------------------------------
00231 
00232                         case AddpvOp:
00233                         CPPAD_ASSERT_UNKNOWN( arg[0] < num_par );
00234                         reverse_addpv_op(
00235                                 d, i_var, arg, parameter, J, Taylor, K, Partial
00236                         );
00237                         break;
00238                         // --------------------------------------------------
00239 
00240                         case AcosOp:
00241                         // results: sqrt(1 - x * x), acos(x) 
00242                         CPPAD_ASSERT_UNKNOWN( i_var < numvar );
00243                         reverse_acos_op(
00244                                 d, i_var, arg[0], J, Taylor, K, Partial
00245                         );
00246                         break;
00247                         // --------------------------------------------------
00248 
00249                         case AsinOp:
00250                         // results: sqrt(1 - x * x), sin(x) 
00251                         CPPAD_ASSERT_UNKNOWN( i_var < numvar );
00252                         reverse_asin_op(
00253                                 d, i_var, arg[0], J, Taylor, K, Partial
00254                         );
00255                         break;
00256                         // --------------------------------------------------
00257 
00258                         case AtanOp:
00259                         // results: 1 + x * x, atan(x) 
00260                         CPPAD_ASSERT_UNKNOWN( i_var < numvar );
00261                         reverse_atan_op(
00262                                 d, i_var, arg[0], J, Taylor, K, Partial
00263                         );
00264                         break;
00265                         // -------------------------------------------------
00266 
00267                         case BeginOp:
00268                         CPPAD_ASSERT_NARG_NRES(op, 0, 1);
00269                         break;
00270                         // --------------------------------------------------
00271 
00272                         case CSumOp:
00273                         // CSumOp has a variable number of arguments and
00274                         // next_reverse thinks it one has one argument.
00275                         // We must inform next_reverse of this special case.
00276                         Rec->reverse_csum(op, arg, i_op, i_var);
00277                         reverse_csum_op(
00278                                 d, i_var, arg, K, Partial
00279                         );
00280                         // end of a cummulative summation
00281                         break;
00282                         // -------------------------------------------------
00283 
00284                         case CExpOp:
00285                         reverse_cond_op(
00286                                 d, 
00287                                 i_var, 
00288                                 arg, 
00289                                 num_par, 
00290                                 parameter, 
00291                                 J, 
00292                                 Taylor,
00293                                 K, 
00294                                 Partial
00295                         );
00296                         break;
00297                         // --------------------------------------------------
00298 
00299                         case ComOp:
00300                         break;
00301                         // --------------------------------------------------
00302 
00303                         case CosOp:
00304                         CPPAD_ASSERT_UNKNOWN( i_var < numvar );
00305                         reverse_cos_op(
00306                                 d, i_var, arg[0], J, Taylor, K, Partial
00307                         );
00308                         break;
00309                         // --------------------------------------------------
00310 
00311                         case CoshOp:
00312                         CPPAD_ASSERT_UNKNOWN( i_var < numvar );
00313                         reverse_cosh_op(
00314                                 d, i_var, arg[0], J, Taylor, K, Partial
00315                         );
00316                         break;
00317                         // --------------------------------------------------
00318 
00319                         case DisOp:
00320                         // Derivative of discrete operation is zero so no
00321                         // contribution passes through this operation. 
00322                         break;
00323                         // --------------------------------------------------
00324 
00325                         case DivvvOp:
00326                         reverse_divvv_op(
00327                                 d, i_var, arg, parameter, J, Taylor, K, Partial
00328                         );
00329                         break;
00330                         // --------------------------------------------------
00331 
00332                         case DivpvOp:
00333                         CPPAD_ASSERT_UNKNOWN( arg[0] < num_par );
00334                         reverse_divpv_op(
00335                                 d, i_var, arg, parameter, J, Taylor, K, Partial
00336                         );
00337                         break;
00338                         // --------------------------------------------------
00339 
00340                         case DivvpOp:
00341                         CPPAD_ASSERT_UNKNOWN( arg[1] < num_par );
00342                         reverse_divvp_op(
00343                                 d, i_var, arg, parameter, J, Taylor, K, Partial
00344                         );
00345                         break;
00346                         // --------------------------------------------------
00347 
00348                         case ExpOp:
00349                         reverse_exp_op(
00350                                 d, i_var, arg[0], J, Taylor, K, Partial
00351                         );
00352                         break;
00353                         // --------------------------------------------------
00354                         case LdpOp:
00355                         reverse_load_op(
00356                                 op, d, i_var, arg, J, Taylor, K, Partial
00357                         );
00358                         break;
00359                         // -------------------------------------------------
00360 
00361                         case LdvOp:
00362                         reverse_load_op(
00363                                 op, d, i_var, arg, J, Taylor, K, Partial
00364                         );
00365                         break;
00366                         // -------------------------------------------------
00367 
00368                         case InvOp:
00369                         break;
00370                         // --------------------------------------------------
00371 
00372                         case LogOp:
00373                         reverse_log_op(
00374                                 d, i_var, arg[0], J, Taylor, K, Partial
00375                         );
00376                         break;
00377                         // --------------------------------------------------
00378 
00379                         case MulvvOp:
00380                         reverse_mulvv_op(
00381                                 d, i_var, arg, parameter, J, Taylor, K, Partial
00382                         );
00383                         break;
00384                         // --------------------------------------------------
00385 
00386                         case MulpvOp:
00387                         CPPAD_ASSERT_UNKNOWN( arg[0] < num_par );
00388                         reverse_mulpv_op(
00389                                 d, i_var, arg, parameter, J, Taylor, K, Partial
00390                         );
00391                         break;
00392                         // --------------------------------------------------
00393 
00394                         case ParOp:
00395                         break;
00396                         // --------------------------------------------------
00397 
00398                         case PowvpOp:
00399                         CPPAD_ASSERT_UNKNOWN( arg[1] < num_par );
00400                         reverse_powvp_op(
00401                                 d, i_var, arg, parameter, J, Taylor, K, Partial
00402                         );
00403                         break;
00404                         // -------------------------------------------------
00405 
00406                         case PowpvOp:
00407                         CPPAD_ASSERT_UNKNOWN( arg[0] < num_par );
00408                         reverse_powpv_op(
00409                                 d, i_var, arg, parameter, J, Taylor, K, Partial
00410                         );
00411                         break;
00412                         // -------------------------------------------------
00413 
00414                         case PowvvOp:
00415                         reverse_powvv_op(
00416                                 d, i_var, arg, parameter, J, Taylor, K, Partial
00417                         );
00418                         break;
00419                         // --------------------------------------------------
00420 
00421                         case PripOp:
00422                         // no result so nothing to do
00423                         break;
00424                         // --------------------------------------------------
00425 
00426                         case PrivOp:
00427                         // no result so nothing to do
00428                         break;
00429 
00430                         // -------------------------------------------------
00431 
00432                         case SinOp:
00433                         CPPAD_ASSERT_UNKNOWN( i_var < numvar );
00434                         reverse_sin_op(
00435                                 d, i_var, arg[0], J, Taylor, K, Partial
00436                         );
00437                         break;
00438                         // -------------------------------------------------
00439 
00440                         case SinhOp:
00441                         CPPAD_ASSERT_UNKNOWN( i_var < numvar );
00442                         reverse_sinh_op(
00443                                 d, i_var, arg[0], J, Taylor, K, Partial
00444                         );
00445                         break;
00446                         // --------------------------------------------------
00447 
00448                         case SqrtOp:
00449                         reverse_sqrt_op(
00450                                 d, i_var, arg[0], J, Taylor, K, Partial
00451                         );
00452                         break;
00453                         // --------------------------------------------------
00454 
00455                         case StppOp:
00456                         break;
00457                         // --------------------------------------------------
00458 
00459                         case StpvOp:
00460                         break;
00461                         // -------------------------------------------------
00462 
00463                         case StvpOp:
00464                         break;
00465                         // -------------------------------------------------
00466 
00467                         case StvvOp:
00468                         break;
00469                         // --------------------------------------------------
00470 
00471                         case SubvvOp:
00472                         reverse_subvv_op(
00473                                 d, i_var, arg, parameter, J, Taylor, K, Partial
00474                         );
00475                         break;
00476                         // --------------------------------------------------
00477 
00478                         case SubpvOp:
00479                         CPPAD_ASSERT_UNKNOWN( arg[0] < num_par );
00480                         reverse_subpv_op(
00481                                 d, i_var, arg, parameter, J, Taylor, K, Partial
00482                         );
00483                         break;
00484                         // --------------------------------------------------
00485 
00486                         case SubvpOp:
00487                         CPPAD_ASSERT_UNKNOWN( arg[1] < num_par );
00488                         reverse_subvp_op(
00489                                 d, i_var, arg, parameter, J, Taylor, K, Partial
00490                         );
00491                         break;
00492                         // --------------------------------------------------
00493 
00494                         case UserOp:
00495                         // start or end an atomic operation sequence
00496                         CPPAD_ASSERT_UNKNOWN( NumRes( UserOp ) == 0 );
00497                         CPPAD_ASSERT_UNKNOWN( NumArg( UserOp ) == 4 );
00498                         if( user_state == user_end )
00499                         {       user_index = arg[0];
00500                                 user_id    = arg[1];
00501                                 user_n     = arg[2];
00502                                 user_m     = arg[3];
00503                                 if(user_ix.size() < user_n)
00504                                         user_ix.resize(user_n);
00505                                 if(user_tx.size() < user_n * user_k1)
00506                                 {       user_tx.resize(user_n * user_k1);
00507                                         user_px.resize(user_n * user_k1);
00508                                 }
00509                                 if(user_ty.size() < user_m * user_k1)
00510                                 {       user_ty.resize(user_m * user_k1);
00511                                         user_py.resize(user_m * user_k1);
00512                                 }
00513                                 user_j     = user_n;
00514                                 user_i     = user_m;
00515                                 user_state = user_ret;
00516                         }
00517                         else
00518                         {       CPPAD_ASSERT_UNKNOWN( user_state == user_start );
00519                                 CPPAD_ASSERT_UNKNOWN( user_index == arg[0] );
00520                                 CPPAD_ASSERT_UNKNOWN( user_id    == arg[1] );
00521                                 CPPAD_ASSERT_UNKNOWN( user_n     == arg[2] );
00522                                 CPPAD_ASSERT_UNKNOWN( user_m     == arg[3] );
00523                                 user_state = user_start;
00524                                 user_state = user_end;
00525 
00526                                 // call users function for this operation
00527                                 user_atomic<Base>::reverse(user_index, user_id,
00528                                         user_k, user_n, user_m, user_tx, user_ty,
00529                                         user_px, user_py
00530                                 );
00531                                 for(j = 0; j < user_n; j++) if( user_ix[j] > 0 )
00532                                 {       for(ell = 0; ell < user_k1; ell++)
00533                                                 Partial[user_ix[j] * K + ell] +=
00534                                                         user_px[j * user_k1 + ell];
00535                                 }
00536                         }
00537                         break;
00538 
00539                         case UsrapOp:
00540                         // parameter argument in an atomic operation sequence
00541                         CPPAD_ASSERT_UNKNOWN( user_state == user_arg );
00542                         CPPAD_ASSERT_UNKNOWN( 0 < user_j && user_j <= user_n );
00543                         CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
00544                         CPPAD_ASSERT_UNKNOWN( arg[0] < num_par );
00545                         --user_j;
00546                         user_ix[user_j] = 0;
00547                         user_tx[user_j * user_k1 + 0] = parameter[ arg[0]];
00548                         for(ell = 1; ell < user_k1; ell++)
00549                                 user_tx[user_j * user_k1 + ell] = Base(0.);
00550                         
00551                         if( user_j == 0 )
00552                                 user_state = user_start;
00553                         break;
00554 
00555                         case UsravOp:
00556                         // variable argument in an atomic operation sequence
00557                         CPPAD_ASSERT_UNKNOWN( user_state == user_arg );
00558                         CPPAD_ASSERT_UNKNOWN( 0 < user_j && user_j <= user_n );
00559                         CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
00560                         CPPAD_ASSERT_UNKNOWN( arg[0] <= i_var );
00561                         CPPAD_ASSERT_UNKNOWN( 0 < arg[0] );
00562                         --user_j;
00563                         user_ix[user_j] = arg[0];
00564                         for(ell = 0; ell < user_k1; ell++)
00565                                 user_tx[user_j*user_k1 + ell] = Taylor[ arg[0] * J + ell];
00566                         if( user_j == 0 )
00567                                 user_state = user_start;
00568                         break;
00569 
00570                         case UsrrpOp:
00571                         // parameter result in an atomic operation sequence
00572                         CPPAD_ASSERT_UNKNOWN( user_state == user_ret );
00573                         CPPAD_ASSERT_UNKNOWN( 0 < user_i && user_i <= user_m );
00574                         CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
00575                         CPPAD_ASSERT_UNKNOWN( arg[0] < num_par );
00576                         --user_i;
00577                         for(ell = 0; ell < user_k1; ell++)
00578                         {       user_py[user_i * user_k1 + ell] = Base(0.);
00579                                 user_ty[user_i * user_k1 + ell] = Base(0.);
00580                         }
00581                         user_ty[user_i * user_k1 + 0] = parameter[ arg[0] ];
00582                         if( user_i == 0 )
00583                                 user_state = user_arg;
00584                         break;
00585 
00586                         case UsrrvOp:
00587                         // variable result in an atomic operation sequence
00588                         CPPAD_ASSERT_UNKNOWN( user_state == user_ret );
00589                         CPPAD_ASSERT_UNKNOWN( 0 < user_i && user_i <= user_m );
00590                         --user_i;
00591                         for(ell = 0; ell < user_k1; ell++)
00592                         {       user_py[user_i * user_k1 + ell] =
00593                                                 Partial[i_var * K + ell];
00594                                 user_ty[user_i * user_k1 + ell] =
00595                                                 Taylor[i_var * J + ell];
00596                         }
00597                         if( user_i == 0 )
00598                                 user_state = user_arg;
00599                         break;
00600                         // ------------------------------------------------------------
00601 
00602                         default:
00603                         CPPAD_ASSERT_UNKNOWN(false);
00604                 }
00605         }
00606 # if CPPAD_REVERSE_SWEEP_TRACE
00607         std::cout << std::endl;
00608 # endif
00609         // values corresponding to BeginOp
00610         CPPAD_ASSERT_UNKNOWN( i_op == 0 );
00611         CPPAD_ASSERT_UNKNOWN( i_var == 0 );
00612 }
00613 
00614 # undef CPPAD_REVERSE_SWEEP_TRACE
00615 
00616 CPPAD_END_NAMESPACE
00617 # endif