CppAD: A C++ Algorithmic Differentiation Package 20110419
|
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