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