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