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