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