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