CppAD: A C++ Algorithmic Differentiation Package  20171217
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
for_two.hpp
Go to the documentation of this file.
1 # ifndef CPPAD_CORE_FOR_TWO_HPP
2 # define CPPAD_CORE_FOR_TWO_HPP
3 
4 /* --------------------------------------------------------------------------
5 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
6 
7 CppAD is distributed under multiple licenses. This distribution is under
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.
12 Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
13 -------------------------------------------------------------------------- */
14 
15 /*
16 $begin ForTwo$$
17 $spell
18  ddy
19  typename
20  Taylor
21  const
22 $$
23 
24 
25 
26 
27 
28 $section Forward Mode Second Partial Derivative Driver$$
29 $mindex order easy$$
30 
31 $head Syntax$$
32 $icode%ddy% = %f%.ForTwo(%x%, %j%, %k%)%$$
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  ddy [ i * p + \ell ]
41  =
42  \DD{ F_i }{ x_{j[ \ell ]} }{ x_{k[ \ell ]} } (x)
43 \] $$
44 for $latex i = 0 , \ldots , m-1$$
45 and $latex \ell = 0 , \ldots , p$$,
46 where $latex p$$ is the size of the vectors $icode j$$ and $icode k$$.
47 
48 $head f$$
49 The object $icode f$$ has prototype
50 $codei%
51  ADFun<%Base%> %f%
52 %$$
53 Note that the $cref ADFun$$ object $icode f$$ is not $code const$$
54 (see $cref/ForTwo Uses Forward/ForTwo/ForTwo 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/ForTwo/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 j$$
69 The argument $icode j$$ has prototype
70 $codei%
71  const %VectorSize_t% &%j%
72 %$$
73 (see $cref/VectorSize_t/ForTwo/VectorSize_t/$$ below)
74 We use $icode p$$ to denote the size of the vector $icode j$$.
75 All of the indices in $icode j$$
76 must be less than $icode n$$; i.e.,
77 for $latex \ell = 0 , \ldots , p-1$$, $latex j[ \ell ] < n$$.
78 
79 $head k$$
80 The argument $icode k$$ has prototype
81 $codei%
82  const %VectorSize_t% &%k%
83 %$$
84 (see $cref/VectorSize_t/ForTwo/VectorSize_t/$$ below)
85 and its size must be equal to $icode p$$,
86 the size of the vector $icode j$$.
87 All of the indices in $icode k$$
88 must be less than $icode n$$; i.e.,
89 for $latex \ell = 0 , \ldots , p-1$$, $latex k[ \ell ] < n$$.
90 
91 $head ddy$$
92 The result $icode ddy$$ has prototype
93 $codei%
94  %VectorBase% %ddy%
95 %$$
96 (see $cref/VectorBase/ForTwo/VectorBase/$$ below)
97 and its size is $latex m * p$$.
98 It contains the requested partial derivatives; to be specific,
99 for $latex i = 0 , \ldots , m - 1 $$
100 and $latex \ell = 0 , \ldots , p - 1$$
101 $latex \[
102  ddy [ i * p + \ell ]
103  =
104  \DD{ F_i }{ x_{j[ \ell ]} }{ x_{k[ \ell ]} } (x)
105 \] $$
106 
107 $head VectorBase$$
108 The type $icode VectorBase$$ must be a $cref SimpleVector$$ class with
109 $cref/elements of type Base/SimpleVector/Elements of Specified Type/$$.
110 The routine $cref CheckSimpleVector$$ will generate an error message
111 if this is not the case.
112 
113 $head VectorSize_t$$
114 The type $icode VectorSize_t$$ must be a $cref SimpleVector$$ class with
115 $cref/elements of type size_t/SimpleVector/Elements of Specified Type/$$.
116 The routine $cref CheckSimpleVector$$ will generate an error message
117 if this is not the case.
118 
119 $head ForTwo Uses Forward$$
120 After each call to $cref Forward$$,
121 the object $icode f$$ contains the corresponding
122 $cref/Taylor coefficients/glossary/Taylor Coefficient/$$.
123 After a call to $code ForTwo$$,
124 the zero order Taylor coefficients correspond to
125 $icode%f%.Forward(0, %x%)%$$
126 and the other coefficients are unspecified.
127 
128 $head Examples$$
129 $children%
130  example/general/for_two.cpp
131 %$$
132 The routine
133 $cref/ForTwo/for_two.cpp/$$ is both an example and test.
134 It returns $code true$$, if it succeeds and $code false$$ otherwise.
135 
136 $end
137 -----------------------------------------------------------------------------
138 */
139 
140 // BEGIN CppAD namespace
141 namespace CppAD {
142 
143 template <typename Base>
144 template <typename VectorBase, typename VectorSize_t>
146  const VectorBase &x,
147  const VectorSize_t &j,
148  const VectorSize_t &k)
149 { size_t i;
150  size_t j1;
151  size_t k1;
152  size_t l;
153 
154  size_t n = Domain();
155  size_t m = Range();
156  size_t p = j.size();
157 
158  // check VectorBase is Simple Vector class with Base type elements
159  CheckSimpleVector<Base, VectorBase>();
160 
161  // check VectorSize_t is Simple Vector class with size_t elements
162  CheckSimpleVector<size_t, VectorSize_t>();
163 
165  x.size() == n,
166  "ForTwo: Length of x not equal domain dimension for f."
167  );
169  j.size() == k.size(),
170  "ForTwo: Lenght of the j and k vectors are not equal."
171  );
172  // point at which we are evaluating the second partials
173  Forward(0, x);
174 
175 
176  // dimension the return value
177  VectorBase ddy(m * p);
178 
179  // allocate memory to hold all possible diagonal Taylor coefficients
180  // (for large sparse cases, this is not efficient)
181  VectorBase D(m * n);
182 
183  // boolean flag for which diagonal coefficients are computed
184  CppAD::vector<bool> c(n);
185  for(j1 = 0; j1 < n; j1++)
186  c[j1] = false;
187 
188  // direction vector in argument space
189  VectorBase dx(n);
190  for(j1 = 0; j1 < n; j1++)
191  dx[j1] = Base(0.0);
192 
193  // result vector in range space
194  VectorBase dy(m);
195 
196  // compute the diagonal coefficients that are needed
197  for(l = 0; l < p; l++)
198  { j1 = j[l];
199  k1 = k[l];
201  j1 < n,
202  "ForTwo: an element of j not less than domain dimension for f."
203  );
205  k1 < n,
206  "ForTwo: an element of k not less than domain dimension for f."
207  );
208  size_t count = 2;
209  while(count)
210  { count--;
211  if( ! c[j1] )
212  { // diagonal term in j1 direction
213  c[j1] = true;
214  dx[j1] = Base(1.0);
215  Forward(1, dx);
216 
217  dx[j1] = Base(0.0);
218  dy = Forward(2, dx);
219  for(i = 0; i < m; i++)
220  D[i * n + j1 ] = dy[i];
221  }
222  j1 = k1;
223  }
224  }
225  // compute all the requested cross partials
226  for(l = 0; l < p; l++)
227  { j1 = j[l];
228  k1 = k[l];
229  if( j1 == k1 )
230  { for(i = 0; i < m; i++)
231  ddy[i * p + l] = Base(2.0) * D[i * n + j1];
232  }
233  else
234  {
235  // cross term in j1 and k1 directions
236  dx[j1] = Base(1.0);
237  dx[k1] = Base(1.0);
238  Forward(1, dx);
239 
240  dx[j1] = Base(0.0);
241  dx[k1] = Base(0.0);
242  dy = Forward(2, dx);
243 
244  // place result in return value
245  for(i = 0; i < m; i++)
246  ddy[i * p + l] = dy[i] - D[i*n+j1] - D[i*n+k1];
247 
248  }
249  }
250  return ddy;
251 }
252 
253 } // END CppAD namespace
254 
255 # endif
#define CPPAD_ASSERT_KNOWN(exp, msg)
Check that exp is true, if not print msg and terminate execution.
VectorBase ForTwo(const VectorBase &x, const VectorSize_t &J, const VectorSize_t &K)
forward mode calculation of a subset of second order partials
Definition: for_two.hpp:145