CppAD: A C++ Algorithmic Differentiation Package  20130102
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-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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines