CppAD: A C++ Algorithmic Differentiation Package 20110419
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_REV_TWO_INCLUDED 00003 # define CPPAD_REV_TWO_INCLUDED 00004 00005 /* -------------------------------------------------------------------------- 00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-09 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 RevTwo$$ 00018 $spell 00019 ddw 00020 typename 00021 Taylor 00022 const 00023 $$ 00024 00025 00026 $index partial, second order driver$$ 00027 $index second, order partial driver$$ 00028 $index driver, second order partial$$ 00029 00030 $index easy, partial$$ 00031 $index driver, easy partial$$ 00032 $index partial, easy$$ 00033 00034 00035 $section Reverse Mode Second Partial Derivative Driver$$ 00036 00037 $head Syntax$$ 00038 $syntax%%ddw% = %f%.RevTwo(%x%, %i%, %j%)%$$ 00039 00040 00041 $head Purpose$$ 00042 We use $latex F : B^n \rightarrow B^m$$ to denote the 00043 $xref/glossary/AD Function/AD function/$$ corresponding to $italic f$$. 00044 The syntax above sets 00045 $latex \[ 00046 ddw [ k * p + \ell ] 00047 = 00048 \DD{ F_{i[ \ell ]} }{ x_{j[ \ell ]} }{ x_k } (x) 00049 \] $$ 00050 for $latex k = 0 , \ldots , n-1$$ 00051 and $latex \ell = 0 , \ldots , p$$, 00052 where $latex p$$ is the size of the vectors $italic i$$ and $italic j$$. 00053 00054 $head f$$ 00055 The object $italic f$$ has prototype 00056 $syntax% 00057 ADFun<%Base%> %f% 00058 %$$ 00059 Note that the $xref/ADFun/$$ object $italic f$$ is not $code const$$ 00060 (see $xref/RevTwo/RevTwo Uses Forward/RevTwo Uses Forward/$$ below). 00061 00062 $head x$$ 00063 The argument $italic x$$ has prototype 00064 $syntax% 00065 const %VectorBase% &%x% 00066 %$$ 00067 (see $xref/RevTwo/VectorBase/VectorBase/$$ below) 00068 and its size 00069 must be equal to $italic n$$, the dimension of the 00070 $xref/seq_property/Domain/domain/$$ space for $italic f$$. 00071 It specifies 00072 that point at which to evaluate the partial derivatives listed above. 00073 00074 $head i$$ 00075 The argument $italic i$$ has prototype 00076 $syntax% 00077 const %VectorSize_t% &%i% 00078 %$$ 00079 (see $xref/RevTwo/VectorSize_t/VectorSize_t/$$ below) 00080 We use $italic p$$ to denote the size of the vector $italic i$$. 00081 All of the indices in $italic i$$ 00082 must be less than $italic m$$, the dimension of the 00083 $xref/seq_property/Range/range/$$ space for $italic f$$; i.e., 00084 for $latex \ell = 0 , \ldots , p-1$$, $latex i[ \ell ] < m$$. 00085 00086 $head j$$ 00087 The argument $italic j$$ has prototype 00088 $syntax% 00089 const %VectorSize_t% &%j% 00090 %$$ 00091 (see $xref/RevTwo/VectorSize_t/VectorSize_t/$$ below) 00092 and its size must be equal to $italic p$$, 00093 the size of the vector $italic i$$. 00094 All of the indices in $italic j$$ 00095 must be less than $italic n$$; i.e., 00096 for $latex \ell = 0 , \ldots , p-1$$, $latex j[ \ell ] < n$$. 00097 00098 $head ddw$$ 00099 The result $italic ddw$$ has prototype 00100 $syntax% 00101 %VectorBase% %ddw% 00102 %$$ 00103 (see $xref/RevTwo/VectorBase/VectorBase/$$ below) 00104 and its size is $latex n * p$$. 00105 It contains the requested partial derivatives; to be specific, 00106 for $latex k = 0 , \ldots , n - 1 $$ 00107 and $latex \ell = 0 , \ldots , p - 1$$ 00108 $latex \[ 00109 ddw [ k * p + \ell ] 00110 = 00111 \DD{ F_{i[ \ell ]} }{ x_{j[ \ell ]} }{ x_k } (x) 00112 \] $$ 00113 00114 $head VectorBase$$ 00115 The type $italic VectorBase$$ must be a $xref/SimpleVector/$$ class with 00116 $xref/SimpleVector/Elements of Specified Type/elements of type Base/$$. 00117 The routine $xref/CheckSimpleVector/$$ will generate an error message 00118 if this is not the case. 00119 00120 $head VectorSize_t$$ 00121 The type $italic VectorSize_t$$ must be a $xref/SimpleVector/$$ class with 00122 $xref/SimpleVector/Elements of Specified Type/elements of type size_t/$$. 00123 The routine $xref/CheckSimpleVector/$$ will generate an error message 00124 if this is not the case. 00125 00126 $head RevTwo Uses Forward$$ 00127 After each call to $xref/Forward/$$, 00128 the object $italic f$$ contains the corresponding 00129 $xref/glossary/Taylor Coefficient/Taylor coefficients/$$. 00130 After $code RevTwo$$, 00131 the previous calls to $xref/Forward/$$ are undefined. 00132 00133 $head Examples$$ 00134 $children% 00135 example/rev_two.cpp 00136 %$$ 00137 The routine 00138 $xref/RevTwo.cpp//RevTwo/$$ is both an example and test. 00139 It returns $code true$$, if it succeeds and $code false$$ otherwise. 00140 00141 $end 00142 ----------------------------------------------------------------------------- 00143 */ 00144 00145 // BEGIN CppAD namespace 00146 namespace CppAD { 00147 00148 template <typename Base> 00149 template <typename VectorBase, typename VectorSize_t> 00150 VectorBase ADFun<Base>::RevTwo( 00151 const VectorBase &x, 00152 const VectorSize_t &i, 00153 const VectorSize_t &j) 00154 { size_t i1; 00155 size_t j1; 00156 size_t k; 00157 size_t l; 00158 00159 size_t n = Domain(); 00160 size_t m = Range(); 00161 size_t p = i.size(); 00162 00163 // check VectorBase is Simple Vector class with Base elements 00164 CheckSimpleVector<Base, VectorBase>(); 00165 00166 // check VectorSize_t is Simple Vector class with size_t elements 00167 CheckSimpleVector<size_t, VectorSize_t>(); 00168 00169 CPPAD_ASSERT_KNOWN( 00170 x.size() == n, 00171 "RevTwo: Length of x not equal domain dimension for f." 00172 ); 00173 CPPAD_ASSERT_KNOWN( 00174 i.size() == j.size(), 00175 "RevTwo: Lenght of the i and j vectors are not equal." 00176 ); 00177 // point at which we are evaluating the second partials 00178 Forward(0, x); 00179 00180 // dimension the return value 00181 VectorBase ddw(n * p); 00182 00183 // direction vector in argument space 00184 VectorBase dx(n); 00185 for(j1 = 0; j1 < n; j1++) 00186 dx[j1] = Base(0); 00187 00188 // direction vector in range space 00189 VectorBase w(m); 00190 for(i1 = 0; i1 < m; i1++) 00191 w[i1] = Base(0); 00192 00193 // place to hold the results of a reverse calculation 00194 VectorBase r(n * 2); 00195 00196 // check the indices in i and j 00197 for(l = 0; l < p; l++) 00198 { i1 = i[l]; 00199 j1 = j[l]; 00200 CPPAD_ASSERT_KNOWN( 00201 i1 < m, 00202 "RevTwo: an eleemnt of i not less than range dimension for f." 00203 ); 00204 CPPAD_ASSERT_KNOWN( 00205 j1 < n, 00206 "RevTwo: an element of j not less than domain dimension for f." 00207 ); 00208 } 00209 00210 // loop over all forward directions 00211 for(j1 = 0; j1 < n; j1++) 00212 { // first order forward mode calculation done 00213 bool first_done = false; 00214 for(l = 0; l < p; l++) if( j[l] == j1 ) 00215 { if( ! first_done ) 00216 { first_done = true; 00217 00218 // first order forward mode in j1 direction 00219 dx[j1] = Base(1); 00220 Forward(1, dx); 00221 dx[j1] = Base(0); 00222 } 00223 // execute a reverse in this component direction 00224 i1 = i[l]; 00225 w[i1] = Base(1); 00226 r = Reverse(2, w); 00227 w[i1] = Base(0); 00228 00229 // place the reverse result in return value 00230 for(k = 0; k < n; k++) 00231 ddw[k * p + l] = r[k * 2 + 1]; 00232 } 00233 } 00234 return ddw; 00235 } 00236 00237 } // END CppAD namespace 00238 00239 # endif