CppAD: A C++ Algorithmic Differentiation Package 20110419
reverse.hpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 # ifndef CPPAD_REVERSE_INCLUDED
00003 # define CPPAD_REVERSE_INCLUDED
00004 
00005 /* --------------------------------------------------------------------------
00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-10 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 $begin Reverse$$
00018 $spell
00019         typename
00020         xk
00021         xp
00022         dw
00023         Ind
00024         uj
00025         std
00026         arg
00027         const
00028         taylor_
00029 $$
00030 
00031 $section Reverse Mode$$ 
00032 
00033 $childtable%
00034         omh/reverse.omh
00035 %$$
00036 
00037 $end
00038 -----------------------------------------------------------------------------
00039 */
00040 # include <algorithm>
00041 
00042 CPPAD_BEGIN_NAMESPACE
00043 /*!
00044 \file reverse.hpp
00045 Compute derivatives using reverse mode.
00046 */
00047 
00048 
00049 /*!
00050 Use reverse mode to compute derivative of forward mode Taylor coefficients.
00051 
00052 The function 
00053 \f$ X : {\rm R} \times {\rm R}^{n \times p} \rightarrow {\rm R} \f$ 
00054 is defined by
00055 \f[
00056 X(t , u) = \sum_{k=0}^{p-1} u^{(k)} t^k
00057 \f]
00058 The function 
00059 \f$ Y : {\rm R} \times {\rm R}^{n \times p} \rightarrow {\rm R} \f$ 
00060 is defined by
00061 \f[
00062 Y(t , u) = F[ X(t, u) ]
00063 \f]
00064 The function 
00065 \f$ W : {\rm R}^{n \times p} \rightarrow {\rm R} \f$ is defined by
00066 \f[
00067 W(u) = \sum_{k=0}^{p-1} ( w^{(k)} )^{\rm T} 
00068 \frac{1}{k !} \frac{ \partial^k } { t^k } Y(0, u)
00069 \f]
00070 
00071 \tparam Base
00072 base type for the operator; i.e., this operation sequence was recorded
00073 using AD< \a Base > and computations by this routine are done using type 
00074 \a Base.
00075 
00076 \tparam VectorBase
00077 is a Simple Vector class with elements of type \a Base.
00078 
00079 \param p
00080 is the number of the number of Taylor coefficients that are being
00081 differentiated (per variable).
00082 
00083 \param w
00084 is the weighting for each of the Taylor coefficients corresponding
00085 to dependent variables.
00086 If the argument \a w has size <tt>m * p </tt>,
00087 for \f$ k = 0 , \ldots , p-1 \f$ and \f$ i = 0, \ldots , m-1 \f$,
00088 \f[
00089         w_i^{(k)} = w [ i * p + k ]
00090 \f]
00091 If the argument \a w has size \c m ,
00092 for \f$ k = 0 , \ldots , p-1 \f$ and \f$ i = 0, \ldots , m-1 \f$,
00093 \f[
00094 w_i^{(k)} = \left\{ \begin{array}{ll}
00095         w [ i ] & {\rm if} \; k = p-1
00096         \\
00097         0       & {\rm otherwise}
00098 \end{array} \right.
00099 \f]
00100 
00101 \return
00102 Is a vector \f$ dw \f$ such that
00103 for \f$ j = 0 , \ldots , n-1 \f$ and
00104 \f$ k = 0 , \ldots , p-1 \f$
00105 \f[ 
00106         dw[ j * p + k ] = W^{(1)} ( x )_{j,k}
00107 \f]
00108 where the matrix \f$ x \f$ is the value for \f$ u \f$
00109 that corresponding to the forward mode Taylor coefficients
00110 for the independent variables as specified by previous calls to Forward.
00111 
00112 */
00113 template <typename Base>
00114 template <typename VectorBase>
00115 VectorBase ADFun<Base>::Reverse(size_t p, const VectorBase &w) 
00116 {       // temporary indices
00117         size_t i, j, k;
00118 
00119         // number of independent variables
00120         size_t n = ind_taddr_.size();
00121 
00122         // number of dependent variables
00123         size_t m = dep_taddr_.size();
00124 
00125         Base *Partial = CPPAD_NULL;
00126         Partial       = CPPAD_TRACK_NEW_VEC(total_num_var_ * p, Partial);
00127 
00128         // update maximum memory requirement
00129         // memoryMax = std::max( memoryMax, 
00130         //      Memory() + total_num_var_ * p * sizeof(Base)
00131         // );
00132 
00133         // check VectorBase is Simple Vector class with Base type elements
00134         CheckSimpleVector<Base, VectorBase>();
00135 
00136         CPPAD_ASSERT_KNOWN(
00137                 w.size() == m || w.size() == (m * p),
00138                 "Argument w to Reverse does not have length equal to\n"
00139                 "the dimension of the range for the corresponding ADFun."
00140         );
00141         CPPAD_ASSERT_KNOWN(
00142                 p > 0,
00143                 "The first argument to Reverse must be greater than zero."
00144         );  
00145         CPPAD_ASSERT_KNOWN(
00146                 taylor_per_var_ >= p,
00147                 "Less that p taylor_ coefficients are currently stored"
00148                 " in this ADFun object."
00149         );  
00150 
00151         // initialize entire Partial matrix to zero
00152         for(i = 0; i < total_num_var_; i++)
00153                 for(j = 0; j < p; j++)
00154                         Partial[i * p + j] = Base(0);
00155 
00156         // set the dependent variable direction
00157         // (use += because two dependent variables can point to same location)
00158         for(i = 0; i < m; i++)
00159         {       CPPAD_ASSERT_UNKNOWN( dep_taddr_[i] < total_num_var_ );
00160                 if( w.size() == m )
00161                         Partial[dep_taddr_[i] * p + p - 1] += w[i];
00162                 else
00163                 {       for(k = 0; k < p; k++)
00164                                 Partial[ dep_taddr_[i] * p + k ] = w[i * p + k ];
00165                 }
00166         }
00167 
00168         // evaluate the derivatives
00169         ReverseSweep(
00170                 p - 1,
00171                 n,
00172                 total_num_var_,
00173                 &play_,
00174                 taylor_col_dim_,
00175                 taylor_,
00176                 p,
00177                 Partial
00178         );
00179 
00180         // return the derivative values
00181         VectorBase value(n * p);
00182         for(j = 0; j < n; j++)
00183         {       CPPAD_ASSERT_UNKNOWN( ind_taddr_[j] < total_num_var_ );
00184 
00185                 // independent variable taddr equals its operator taddr 
00186                 CPPAD_ASSERT_UNKNOWN( play_.GetOp( ind_taddr_[j] ) == InvOp );
00187 
00188                 // by the Reverse Identity Theorem 
00189                 // partial of y^{(k)} w.r.t. u^{(0)} is equal to
00190                 // partial of y^{(p-1)} w.r.t. u^{(p - 1 - k)}
00191                 if( w.size() == m )
00192                 {       for(k = 0; k < p; k++)
00193                                 value[j * p + k ] = 
00194                                         Partial[ind_taddr_[j] * p + p - 1 - k];
00195                 }
00196                 else
00197                 {       for(k = 0; k < p; k++)
00198                                 value[j * p + k ] =
00199                                         Partial[ind_taddr_[j] * p + k];
00200                 }
00201         }
00202 
00203         // done with the Partial array
00204         CPPAD_TRACK_DEL_VEC(Partial);
00205 
00206         return value;
00207 }
00208         
00209 
00210 } // END CppAD namespace
00211         
00212 
00213 # endif