CppAD: A C++ Algorithmic Differentiation Package  20171217
rev_two.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 RevTwo$$17 spell 18 ddw 19 typename 20 Taylor 21 const 22$$ 23 24 25 26 27 28$section Reverse Mode Second Partial Derivative Driver$$29 mindex order easy$$
30
31 $head Syntax$$32 icode%ddw% = %f%.RevTwo(%x%, %i%, %j%)%$$ 33 34 35$head Purpose$$36 We use latex F : B^n \rightarrow B^m$$ to denote the
37 $cref/AD function/glossary/AD Function/$$corresponding to icode f$$. 38 The syntax above sets 39$latex $40 ddw [ k * p + \ell ] 41 = 42 \DD{ F_{i[ \ell ]} }{ x_{j[ \ell ]} }{ x_k } (x) 43$ $$44 for latex k = 0 , \ldots , n-1$$
45 and $latex \ell = 0 , \ldots , p$$, 46 where latex p$$ is the size of the vectors$icode i$$and icode j$$.
47
48 $head f$$49 The object icode f$$ has prototype 50$codei%
52 %$$53 Note that the cref ADFun$$ object $icode f$$is not code const$$ 54 (see$cref/RevTwo Uses Forward/RevTwo/RevTwo Uses Forward/$$below). 55 56 head x$$
57 The argument $icode x$$has prototype 58 codei% 59 const %VectorBase% &%x% 60 %$$ 61 (see$cref/VectorBase/RevTwo/VectorBase/$$below) 62 and its size 63 must be equal to icode n$$, the dimension of the
64 $cref/domain/seq_property/Domain/$$space for icode f$$. 65 It specifies 66 that point at which to evaluate the partial derivatives listed above. 67 68$head i$$69 The argument icode i$$ has prototype
70 $codei% 71 const %VectorSize_t% &%i% 72 %$$73 (see cref/VectorSize_t/RevTwo/VectorSize_t/$$ below) 74 We use$icode p$$to denote the size of the vector icode i$$.
75 All of the indices in $icode i$$76 must be less than icode m$$, the dimension of the 77$cref/range/seq_property/Range/$$space for icode f$$; i.e.,
78 for $latex \ell = 0 , \ldots , p-1$$, latex i[ \ell ] < m$$. 79 80$head j$$81 The argument icode j$$ has prototype
82 $codei% 83 const %VectorSize_t% &%j% 84 %$$85 (see cref/VectorSize_t/RevTwo/VectorSize_t/$$ below) 86 and its size must be equal to$icode p$$, 87 the size of the vector icode i$$.
88 All of the indices in $icode j$$89 must be less than icode n$$; i.e., 90 for$latex \ell = 0 , \ldots , p-1$$, latex j[ \ell ] < n$$.
91
92 $head ddw$$93 The result icode ddw$$ has prototype 94$codei%
95  %VectorBase% %ddw%
96 %$$97 (see cref/VectorBase/RevTwo/VectorBase/$$ below)
98 and its size is $latex n * p$$. 99 It contains the requested partial derivatives; to be specific, 100 for latex k = 0 , \ldots , n - 1$$ 101 and$latex \ell = 0 , \ldots , p - 1$$102 latex $103 ddw [ k * p + \ell ] 104 = 105 \DD{ F_{i[ \ell ]} }{ x_{j[ \ell ]} }{ x_k } (x) 106$$$
107
108 $head VectorBase$$109 The type icode VectorBase$$ must be a$cref SimpleVector$$class with 110 cref/elements of type Base/SimpleVector/Elements of Specified Type/$$.
111 The routine $cref CheckSimpleVector$$will generate an error message 112 if this is not the case. 113 114 head VectorSize_t$$ 115 The type$icode VectorSize_t$$must be a cref SimpleVector$$ class with
116 $cref/elements of type size_t/SimpleVector/Elements of Specified Type/$$. 117 The routine cref CheckSimpleVector$$ will generate an error message 118 if this is not the case. 119 120$head RevTwo Uses Forward$$121 After each call to cref Forward$$,
122 the object $icode f$$contains the corresponding 123 cref/Taylor coefficients/glossary/Taylor Coefficient/$$. 124 After a call to$code RevTwo$$, 125 the zero order Taylor coefficients correspond to 126 icode%f%.Forward(0, %x%)%$$
127 and the other coefficients are unspecified.
128
129 $head Examples$$130 children% 131 example/general/rev_two.cpp 132 %$$ 133 The routine 134$cref/RevTwo/rev_two.cpp/$$is both an example and test. 135 It returns code true$$, if it succeeds and $code false$$otherwise. 136 137$end
138 -----------------------------------------------------------------------------
139 */
140
143
144 template <typename Base>
145 template <typename VectorBase, typename VectorSize_t>
147  const VectorBase &x,
148  const VectorSize_t &i,
149  const VectorSize_t &j)
150 { size_t i1;
151  size_t j1;
152  size_t k;
153  size_t l;
154
155  size_t n = Domain();
156  size_t m = Range();
157  size_t p = i.size();
158
159  // check VectorBase is Simple Vector class with Base elements
160  CheckSimpleVector<Base, VectorBase>();
161
162  // check VectorSize_t is Simple Vector class with size_t elements
163  CheckSimpleVector<size_t, VectorSize_t>();
164
166  x.size() == n,
167  "RevTwo: Length of x not equal domain dimension for f."
168  );
170  i.size() == j.size(),
171  "RevTwo: Lenght of the i and j vectors are not equal."
172  );
173  // point at which we are evaluating the second partials
174  Forward(0, x);
175
176  // dimension the return value
177  VectorBase ddw(n * p);
178
179  // direction vector in argument space
180  VectorBase dx(n);
181  for(j1 = 0; j1 < n; j1++)
182  dx[j1] = Base(0.0);
183
184  // direction vector in range space
185  VectorBase w(m);
186  for(i1 = 0; i1 < m; i1++)
187  w[i1] = Base(0.0);
188
189  // place to hold the results of a reverse calculation
190  VectorBase r(n * 2);
191
192  // check the indices in i and j
193  for(l = 0; l < p; l++)
194  { i1 = i[l];
195  j1 = j[l];
197  i1 < m,
198  "RevTwo: an eleemnt of i not less than range dimension for f."
199  );
201  j1 < n,
202  "RevTwo: an element of j not less than domain dimension for f."
203  );
204  }
205
206  // loop over all forward directions
207  for(j1 = 0; j1 < n; j1++)
208  { // first order forward mode calculation done
209  bool first_done = false;
210  for(l = 0; l < p; l++) if( j[l] == j1 )
211  { if( ! first_done )
212  { first_done = true;
213
214  // first order forward mode in j1 direction
215  dx[j1] = Base(1.0);
216  Forward(1, dx);
217  dx[j1] = Base(0.0);
218  }
219  // execute a reverse in this component direction
220  i1 = i[l];
221  w[i1] = Base(1.0);
222  r = Reverse(2, w);
223  w[i1] = Base(0.0);
224
225  // place the reverse result in return value
226  for(k = 0; k < n; k++)
227  ddw[k * p + l] = r[k * 2 + 1];
228  }
229  }
230  return ddw;
231 }
232
233 } // END CppAD namespace
234
235 # endif