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