CppAD: A C++ Algorithmic Differentiation Package
20130918
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_FOR_TWO_INCLUDED 00003 # define CPPAD_FOR_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 ForTwo$$ 00018 $spell 00019 ddy 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 Forward Mode Second Partial Derivative Driver$$ 00036 00037 $head Syntax$$ 00038 $icode%ddy% = %f%.ForTwo(%x%, %j%, %k%)%$$ 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 ddy [ i * p + \ell ] 00047 = 00048 \DD{ F_i }{ x_{j[ \ell ]} }{ x_{k[ \ell ]} } (x) 00049 \] $$ 00050 for $latex i = 0 , \ldots , m-1$$ 00051 and $latex \ell = 0 , \ldots , p$$, 00052 where $latex p$$ is the size of the vectors $icode j$$ and $icode k$$. 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/ForTwo Uses Forward/ForTwo/ForTwo 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/ForTwo/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 j$$ 00075 The argument $icode j$$ has prototype 00076 $codei% 00077 const %VectorSize_t% &%j% 00078 %$$ 00079 (see $cref/VectorSize_t/ForTwo/VectorSize_t/$$ below) 00080 We use $icode p$$ to denote the size of the vector $icode j$$. 00081 All of the indices in $icode j$$ 00082 must be less than $icode n$$; i.e., 00083 for $latex \ell = 0 , \ldots , p-1$$, $latex j[ \ell ] < n$$. 00084 00085 $head k$$ 00086 The argument $icode k$$ has prototype 00087 $codei% 00088 const %VectorSize_t% &%k% 00089 %$$ 00090 (see $cref/VectorSize_t/ForTwo/VectorSize_t/$$ below) 00091 and its size must be equal to $icode p$$, 00092 the size of the vector $icode j$$. 00093 All of the indices in $icode k$$ 00094 must be less than $icode n$$; i.e., 00095 for $latex \ell = 0 , \ldots , p-1$$, $latex k[ \ell ] < n$$. 00096 00097 $head ddy$$ 00098 The result $icode ddy$$ has prototype 00099 $codei% 00100 %VectorBase% %ddy% 00101 %$$ 00102 (see $cref/VectorBase/ForTwo/VectorBase/$$ below) 00103 and its size is $latex m * p$$. 00104 It contains the requested partial derivatives; to be specific, 00105 for $latex i = 0 , \ldots , m - 1 $$ 00106 and $latex \ell = 0 , \ldots , p - 1$$ 00107 $latex \[ 00108 ddy [ i * p + \ell ] 00109 = 00110 \DD{ F_i }{ x_{j[ \ell ]} }{ x_{k[ \ell ]} } (x) 00111 \] $$ 00112 00113 $head VectorBase$$ 00114 The type $icode VectorBase$$ must be a $cref SimpleVector$$ class with 00115 $cref/elements of type Base/SimpleVector/Elements of Specified Type/$$. 00116 The routine $cref CheckSimpleVector$$ will generate an error message 00117 if this is not the case. 00118 00119 $head VectorSize_t$$ 00120 The type $icode VectorSize_t$$ must be a $cref SimpleVector$$ class with 00121 $cref/elements of type size_t/SimpleVector/Elements of Specified Type/$$. 00122 The routine $cref CheckSimpleVector$$ will generate an error message 00123 if this is not the case. 00124 00125 $head ForTwo Uses Forward$$ 00126 After each call to $cref Forward$$, 00127 the object $icode f$$ contains the corresponding 00128 $cref/Taylor coefficients/glossary/Taylor Coefficient/$$. 00129 After a call to $code ForTwo$$, 00130 the zero order Taylor coefficients correspond to 00131 $icode%f%.Forward(0, %x%)%$$ 00132 and the other coefficients are unspecified. 00133 00134 $head Examples$$ 00135 $children% 00136 example/for_two.cpp 00137 %$$ 00138 The routine 00139 $cref/ForTwo/for_two.cpp/$$ is both an example and test. 00140 It returns $code true$$, if it succeeds and $code false$$ otherwise. 00141 00142 $end 00143 ----------------------------------------------------------------------------- 00144 */ 00145 00146 // BEGIN CppAD namespace 00147 namespace CppAD { 00148 00149 template <typename Base> 00150 template <typename VectorBase, typename VectorSize_t> 00151 VectorBase ADFun<Base>::ForTwo( 00152 const VectorBase &x, 00153 const VectorSize_t &j, 00154 const VectorSize_t &k) 00155 { size_t i; 00156 size_t j1; 00157 size_t k1; 00158 size_t l; 00159 00160 size_t n = Domain(); 00161 size_t m = Range(); 00162 size_t p = j.size(); 00163 00164 // check VectorBase is Simple Vector class with Base type elements 00165 CheckSimpleVector<Base, VectorBase>(); 00166 00167 // check VectorSize_t is Simple Vector class with size_t elements 00168 CheckSimpleVector<size_t, VectorSize_t>(); 00169 00170 CPPAD_ASSERT_KNOWN( 00171 x.size() == n, 00172 "ForTwo: Length of x not equal domain dimension for f." 00173 ); 00174 CPPAD_ASSERT_KNOWN( 00175 j.size() == k.size(), 00176 "ForTwo: Lenght of the j and k vectors are not equal." 00177 ); 00178 // point at which we are evaluating the second partials 00179 Forward(0, x); 00180 00181 00182 // dimension the return value 00183 VectorBase ddy(m * p); 00184 00185 // allocate memory to hold all possible diagonal Taylor coefficients 00186 // (for large sparse cases, this is not efficient) 00187 VectorBase D(m * n); 00188 00189 // boolean flag for which diagonal coefficients are computed 00190 CppAD::vector<bool> c(n); 00191 for(j1 = 0; j1 < n; j1++) 00192 c[j1] = false; 00193 00194 // direction vector in argument space 00195 VectorBase dx(n); 00196 for(j1 = 0; j1 < n; j1++) 00197 dx[j1] = Base(0); 00198 00199 // result vector in range space 00200 VectorBase dy(m); 00201 00202 // compute the diagonal coefficients that are needed 00203 for(l = 0; l < p; l++) 00204 { j1 = j[l]; 00205 k1 = k[l]; 00206 CPPAD_ASSERT_KNOWN( 00207 j1 < n, 00208 "ForTwo: an element of j not less than domain dimension for f." 00209 ); 00210 CPPAD_ASSERT_KNOWN( 00211 k1 < n, 00212 "ForTwo: an element of k not less than domain dimension for f." 00213 ); 00214 size_t count = 2; 00215 while(count) 00216 { count--; 00217 if( ! c[j1] ) 00218 { // diagonal term in j1 direction 00219 c[j1] = true; 00220 dx[j1] = Base(1); 00221 Forward(1, dx); 00222 00223 dx[j1] = Base(0); 00224 dy = Forward(2, dx); 00225 for(i = 0; i < m; i++) 00226 D[i * n + j1 ] = dy[i]; 00227 } 00228 j1 = k1; 00229 } 00230 } 00231 // compute all the requested cross partials 00232 for(l = 0; l < p; l++) 00233 { j1 = j[l]; 00234 k1 = k[l]; 00235 if( j1 == k1 ) 00236 { for(i = 0; i < m; i++) 00237 ddy[i * p + l] = Base(2) * D[i * n + j1]; 00238 } 00239 else 00240 { 00241 // cross term in j1 and k1 directions 00242 dx[j1] = Base(1); 00243 dx[k1] = Base(1); 00244 Forward(1, dx); 00245 00246 dx[j1] = Base(0); 00247 dx[k1] = Base(0); 00248 dy = Forward(2, dx); 00249 00250 // place result in return value 00251 for(i = 0; i < m; i++) 00252 ddy[i * p + l] = dy[i] - D[i*n+j1] - D[i*n+k1]; 00253 00254 } 00255 } 00256 return ddy; 00257 } 00258 00259 } // END CppAD namespace 00260 00261 # endif