CppAD: A C++ Algorithmic Differentiation Package 20110419
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_JACOBIAN_INCLUDED 00003 # define CPPAD_JACOBIAN_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 Jacobian$$ 00018 $spell 00019 jac 00020 typename 00021 Taylor 00022 Jacobian 00023 DetLu 00024 const 00025 $$ 00026 00027 $index Jacobian, driver$$ 00028 $index first, derivative$$ 00029 $index driver, Jacobian$$ 00030 00031 $section Jacobian: Driver Routine$$ 00032 00033 $head Syntax$$ 00034 $syntax%%jac% = %f%.Jacobian(%x%)%$$ 00035 00036 00037 $head Purpose$$ 00038 We use $latex F : B^n \rightarrow B^m$$ to denote the 00039 $xref/glossary/AD Function/AD function/$$ corresponding to $italic f$$. 00040 The syntax above sets $italic jac$$ to the 00041 Jacobian of $italic F$$ evaluated at $italic x$$; i.e., 00042 $latex \[ 00043 jac = F^{(1)} (x) 00044 \] $$ 00045 00046 $head f$$ 00047 The object $italic f$$ has prototype 00048 $syntax% 00049 ADFun<%Base%> %f% 00050 %$$ 00051 Note that the $xref/ADFun/$$ object $italic f$$ is not $code const$$ 00052 (see $cref/Forward or Reverse/Jacobian/Forward or Reverse/$$ below). 00053 00054 $head x$$ 00055 The argument $italic x$$ has prototype 00056 $syntax% 00057 const %Vector% &%x% 00058 %$$ 00059 (see $xref/Jacobian/Vector/Vector/$$ below) 00060 and its size 00061 must be equal to $italic n$$, the dimension of the 00062 $xref/seq_property/Domain/domain/$$ space for $italic f$$. 00063 It specifies 00064 that point at which to evaluate the Jacobian. 00065 00066 $head jac$$ 00067 The result $italic jac$$ has prototype 00068 $syntax% 00069 %Vector% %jac% 00070 %$$ 00071 (see $xref/Jacobian/Vector/Vector/$$ below) 00072 and its size is $latex m * n$$; i.e., the product of the 00073 $xref/seq_property/Domain/domain/$$ 00074 and 00075 $xref/seq_property/Range/range/$$ 00076 dimensions for $italic f$$. 00077 For $latex i = 0 , \ldots , m - 1 $$ 00078 and $latex j = 0 , \ldots , n - 1$$ 00079 $latex \[. 00080 jac[ i * n + j ] = \D{ F_i }{ x_j } ( x ) 00081 \] $$ 00082 00083 00084 $head Vector$$ 00085 The type $italic Vector$$ must be a $xref/SimpleVector/$$ class with 00086 $xref/SimpleVector/Elements of Specified Type/elements of type/$$ 00087 $italic Base$$. 00088 The routine $xref/CheckSimpleVector/$$ will generate an error message 00089 if this is not the case. 00090 00091 $head Forward or Reverse$$ 00092 This will use order zero Forward mode and either 00093 order one Forward or order one Reverse to compute the Jacobian 00094 (depending on which it estimates will require less work). 00095 After each call to $xref/Forward/$$, 00096 the object $italic f$$ contains the corresponding 00097 $xref/glossary/Taylor Coefficient/Taylor coefficients/$$. 00098 After each call to $code Jacobian$$, 00099 the previous calls to $xref/Forward/$$ are unspecified. 00100 00101 $head Example$$ 00102 $children% 00103 example/jacobian.cpp 00104 %$$ 00105 The routine 00106 $xref/Jacobian.cpp//Jacobian/$$ is both an example and test. 00107 It returns $code true$$, if it succeeds and $code false$$ otherwise. 00108 00109 $end 00110 ----------------------------------------------------------------------------- 00111 */ 00112 00113 // BEGIN CppAD namespace 00114 namespace CppAD { 00115 00116 template <typename Base, typename Vector> 00117 void JacobianFor(ADFun<Base> &f, const Vector &x, Vector &jac) 00118 { size_t i; 00119 size_t j; 00120 00121 size_t m = f.Domain(); 00122 size_t n = f.Range(); 00123 00124 // check Vector is Simple Vector class with Base type elements 00125 CheckSimpleVector<Base, Vector>(); 00126 00127 CPPAD_ASSERT_UNKNOWN( x.size() == f.Domain() ); 00128 CPPAD_ASSERT_UNKNOWN( jac.size() == f.Range() * f.Domain() ); 00129 00130 // argument and result for forward mode calculations 00131 Vector u(m); 00132 Vector v(n); 00133 00134 // initialize all the components 00135 for(j = 0; j < m; j++) 00136 u[j] = Base(0); 00137 00138 // loop through the different coordinate directions 00139 for(j = 0; j < m; j++) 00140 { // set u to the j-th coordinate direction 00141 u[j] = Base(1); 00142 00143 // compute the partial of f w.r.t. this coordinate direction 00144 v = f.Forward(1, u); 00145 00146 // reset u to vector of all zeros 00147 u[j] = Base(0); 00148 00149 // return the result 00150 for(i = 0; i < n; i++) 00151 jac[ i * m + j ] = v[i]; 00152 } 00153 } 00154 template <typename Base, typename Vector> 00155 void JacobianRev(ADFun<Base> &f, const Vector &x, Vector &jac) 00156 { size_t i; 00157 size_t j; 00158 00159 size_t m = f.Domain(); 00160 size_t n = f.Range(); 00161 00162 CPPAD_ASSERT_UNKNOWN( x.size() == f.Domain() ); 00163 CPPAD_ASSERT_UNKNOWN( jac.size() == f.Range() * f.Domain() ); 00164 00165 // argument and result for reverse mode calculations 00166 Vector u(m); 00167 Vector v(n); 00168 00169 // initialize all the components 00170 for(i = 0; i < n; i++) 00171 v[i] = Base(0); 00172 00173 // loop through the different coordinate directions 00174 for(i = 0; i < n; i++) 00175 { if( f.Parameter(i) ) 00176 { // return zero for this component of f 00177 for(j = 0; j < m; j++) 00178 jac[ i * m + j ] = Base(0); 00179 } 00180 else 00181 { 00182 // set v to the i-th coordinate direction 00183 v[i] = Base(1); 00184 00185 // compute the derivative of this component of f 00186 u = f.Reverse(1, v); 00187 00188 // reset v to vector of all zeros 00189 v[i] = Base(0); 00190 00191 // return the result 00192 for(j = 0; j < m; j++) 00193 jac[ i * m + j ] = u[j]; 00194 } 00195 } 00196 } 00197 00198 template <typename Base> 00199 template <typename Vector> 00200 Vector ADFun<Base>::Jacobian(const Vector &x) 00201 { size_t i; 00202 size_t m = Domain(); 00203 size_t n = Range(); 00204 00205 CPPAD_ASSERT_KNOWN( 00206 x.size() == m, 00207 "Jacobian: length of x not equal domain dimension for F" 00208 ); 00209 00210 // point at which we are evaluating the Jacobian 00211 Forward(0, x); 00212 00213 // work factor for forward mode 00214 size_t workForward = m; 00215 00216 // work factor for reverse mode 00217 size_t workReverse = 0; 00218 for(i = 0; i < n; i++) 00219 { if( ! Parameter(i) ) 00220 ++workReverse; 00221 } 00222 00223 // choose the method with the least work 00224 Vector jac( n * m ); 00225 if( workForward <= workReverse ) 00226 JacobianFor(*this, x, jac); 00227 else JacobianRev(*this, x, jac); 00228 00229 return jac; 00230 } 00231 00232 } // END CppAD namespace 00233 00234 # endif