CppAD: A C++ Algorithmic Differentiation Package 20110419
jacobian.hpp
Go to the documentation of this file.
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