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