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