CppAD: A C++ Algorithmic Differentiation Package  20171217
jacobian.hpp
Go to the documentation of this file.
3
4 /* --------------------------------------------------------------------------
6
8 the terms of the
9  Eclipse Public License Version 1.0.
10
11 A copy of this license is included in the COPYING file of this distribution.
13 -------------------------------------------------------------------------- */
14
15 /*
16 $begin Jacobian$$17 spell 18 jac 19 typename 20 Taylor 21 Jacobian 22 DetLu 23 const 24$$ 25 26 27$section Jacobian: Driver Routine$$28 mindex Jacobian first derivative$$
29
30 $head Syntax$$31 icode%jac% = %f%.Jacobian(%x%)%$$ 32 33 34$head Purpose$$35 We use latex F : B^n \rightarrow B^m$$ to denote the
36 $cref/AD function/glossary/AD Function/$$corresponding to icode f$$. 37 The syntax above sets$icode jac$$to the 38 Jacobian of icode F$$ evaluated at $icode x$$; i.e., 39 latex $40 jac = F^{(1)} (x) 41$$$ 42 43$head f$$44 The object icode f$$ has prototype
45 $codei% 46 ADFun<%Base%> %f% 47 %$$48 Note that the cref ADFun$$ object$icode f$$is not code const$$
49 (see $cref/Forward or Reverse/Jacobian/Forward or Reverse/$$below). 50 51 head x$$ 52 The argument$icode x$$has prototype 53 codei% 54 const %Vector% &%x% 55 %$$
56 (see $cref/Vector/Jacobian/Vector/$$below) 57 and its size 58 must be equal to icode n$$, the dimension of the 59$cref/domain/seq_property/Domain/$$space for icode f$$.
60 It specifies
61 that point at which to evaluate the Jacobian.
62
63 $head jac$$64 The result icode jac$$ has prototype 65$codei%
66  %Vector% %jac%
67 %$$68 (see cref/Vector/Jacobian/Vector/$$ below)
69 and its size is $latex m * n$$; i.e., the product of the 70 cref/domain/seq_property/Domain/$$ 71 and 72$cref/range/seq_property/Range/$$73 dimensions for icode f$$.
74 For $latex i = 0 , \ldots , m - 1 $$75 and latex j = 0 , \ldots , n - 1$$ 76$latex $. 77 jac[ i * n + j ] = \D{ F_i }{ x_j } ( x ) 78$ $$79 80 81 head Vector$$
82 The type $icode Vector$$must be a cref SimpleVector$$ class with 83$cref/elements of type/SimpleVector/Elements of Specified Type/$$84 icode Base$$.
85 The routine $cref CheckSimpleVector$$will generate an error message 86 if this is not the case. 87 88 head Forward or Reverse$$ 89 This will use order zero Forward mode and either 90 order one Forward or order one Reverse to compute the Jacobian 91 (depending on which it estimates will require less work). 92 After each call to$cref Forward$$, 93 the object icode f$$ contains the corresponding
94 $cref/Taylor coefficients/glossary/Taylor Coefficient/$$. 95 After a call to code Jacobian$$, 96 the zero order Taylor coefficients correspond to 97$icode%f%.Forward(0, %x%)%$$98 and the other coefficients are unspecified. 99 100 head Example$$
101 $children% 102 example/general/jacobian.cpp 103 %$$104 The routine 105 cref/Jacobian/jacobian.cpp/$$ is both an example and test. 106 It returns$code true$$, if it succeeds and code false$$ otherwise.
107
108 \$end
109 -----------------------------------------------------------------------------
110 */
111
114
115 template <typename Base, typename Vector>
116 void JacobianFor(ADFun<Base> &f, const Vector &x, Vector &jac)
117 { size_t i;
118  size_t j;
119
120  size_t n = f.Domain();
121  size_t m = f.Range();
122
123  // check Vector is Simple Vector class with Base type elements
124  CheckSimpleVector<Base, Vector>();
125
126  CPPAD_ASSERT_UNKNOWN( size_t(x.size()) == f.Domain() );
127  CPPAD_ASSERT_UNKNOWN( size_t(jac.size()) == f.Range() * f.Domain() );
128
129  // argument and result for forward mode calculations
130  Vector u(n);
131  Vector v(m);
132
133  // initialize all the components
134  for(j = 0; j < n; j++)
135  u[j] = Base(0.0);
136
137  // loop through the different coordinate directions
138  for(j = 0; j < n; j++)
139  { // set u to the j-th coordinate direction
140  u[j] = Base(1.0);
141
142  // compute the partial of f w.r.t. this coordinate direction
143  v = f.Forward(1, u);
144
145  // reset u to vector of all zeros
146  u[j] = Base(0.0);
147
148  // return the result
149  for(i = 0; i < m; i++)
150  jac[ i * n + j ] = v[i];
151  }
152 }
153 template <typename Base, typename Vector>
154 void JacobianRev(ADFun<Base> &f, const Vector &x, Vector &jac)
155 { size_t i;
156  size_t j;
157
158  size_t n = f.Domain();
159  size_t m = f.Range();
160
161  CPPAD_ASSERT_UNKNOWN( size_t(x.size()) == f.Domain() );
162  CPPAD_ASSERT_UNKNOWN( size_t(jac.size()) == f.Range() * f.Domain() );
163
164  // argument and result for reverse mode calculations
165  Vector u(n);
166  Vector v(m);
167
168  // initialize all the components
169  for(i = 0; i < m; i++)
170  v[i] = Base(0.0);
171
172  // loop through the different coordinate directions
173  for(i = 0; i < m; i++)
174  { if( f.Parameter(i) )
175  { // return zero for this component of f
176  for(j = 0; j < n; j++)
177  jac[ i * n + j ] = Base(0.0);
178  }
179  else
180  {
181  // set v to the i-th coordinate direction
182  v[i] = Base(1.0);
183
184  // compute the derivative of this component of f
185  u = f.Reverse(1, v);
186
187  // reset v to vector of all zeros
188  v[i] = Base(0.0);
189
190  // return the result
191  for(j = 0; j < n; j++)
192  jac[ i * n + j ] = u[j];
193  }
194  }
195 }
196
197 template <typename Base>
198 template <typename Vector>
200 { size_t i;
201  size_t n = Domain();
202  size_t m = Range();
203
205  size_t(x.size()) == n,
206  "Jacobian: length of x not equal domain dimension for F"
207  );
208
209  // point at which we are evaluating the Jacobian
210  Forward(0, x);
211
212  // work factor for forward mode
213  size_t workForward = n;
214
215  // work factor for reverse mode
216  size_t workReverse = 0;
217  for(i = 0; i < m; i++)
218  { if( ! Parameter(i) )
219  ++workReverse;
220  }
221
222  // choose the method with the least work
223  Vector jac( n * m );
225  if( workForward < workReverse )
226 # else
227  if( workForward <= workReverse )
228 # endif
229  JacobianFor(*this, x, jac);
230  else JacobianRev(*this, x, jac);
231
232  return jac;
233 }
234
235 } // END CppAD namespace
236
237 # endif
Check that exp is true, if not print msg and terminate execution.
size_t Range(void) const
number of dependent variables
Class used to hold function objects.
void JacobianRev(ADFun< Base > &f, const Vector &x, Vector &jac)
Definition: jacobian.hpp:154
VectorBase Jacobian(const VectorBase &x)
calculate entire Jacobian
Check that exp is true, if not terminate execution.
bool Parameter(size_t i)
is variable a parameter