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