CppAD: A C++ Algorithmic Differentiation Package  20171217
forward.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 // documened after Forward but included here so easy to see
19
21 /*!
22 \file forward.hpp
23 User interface to forward mode computations.
24 */
25
26 /*!
27 Multiple orders, one direction, forward mode Taylor coefficieints.
28
29 \tparam Base
30 The type used during the forward mode computations; i.e., the corresponding
31 recording of operations used the type AD<Base>.
32
33 \tparam VectorBase
34 is a Simple Vector class with eleements of type Base.
35
36 \param q
37 is the hightest order for this forward mode computation; i.e.,
38 after this calculation there will be <code>q+1</code>
39 Taylor coefficients per variable.
40
41 \param xq
42 contains Taylor coefficients for the independent variables.
43 The size of xq must either be n or <code>(q+1)*n</code>,
44 We define <code>p = q + 1 - xq.size()/n</code>.
45 For <code>j = 0 , ... , n-1</code>,
46 <code>k = p, ... , q</code>, are
47 <code>xq[ (q+1-p)*j + k - p ]</code>
48 is the k-th order coefficient for the j-th independent variable.
49
50 \param s
51 Is the stream where output corresponding to PriOp operations will written.
52
53 \return
54 contains Taylor coefficients for the dependent variables.
55 The size of the return value y is <code>m*(q+1-p)</code>.
56 For <code>i = 0, ... , m-1</code>,
57 <code>k = p, ..., q</code>,
58 <code>y[(q+1-p)*i + (k-p)]</code>
59 is the k-th order coefficient for the i-th dependent variable.
60
61 \par taylor_
62 The Taylor coefficients up to order p-1 are inputs
63 and the coefficents from order p through q are outputs.
64 Let <code>N = num_var_tape_</code>, and
65 <code>C = cap_order_taylor_</code>.
66 Note that for
67 <code>i = 1 , ..., N-1</code>,
68 <code>k = 0 , ..., q</code>,
69 <code>taylor_[ C*i + k ]</code>
70 is the k-th order cofficent,
71 for the i-th varaible on the tape.
72 (The first independent variable has index one on the tape
73 and there is no variable with index zero.)
74 */
75
76 template <typename Base>
77 template <typename VectorBase>
79  size_t q ,
80  const VectorBase& xq ,
81  std::ostream& s )
82 { // temporary indices
83  size_t i, j, k;
84
85  // number of independent variables
87
88  // number of dependent variables
90
91  // check Vector is Simple Vector class with Base type elements
92  CheckSimpleVector<Base, VectorBase>();
93
94
96  size_t(xq.size()) == n || size_t(xq.size()) == n*(q+1),
97  "Forward(q, xq): xq.size() is not equal n or n*(q+1)"
98  );
99
100  // lowest order we are computing
101  size_t p = q + 1 - size_t(xq.size()) / n;
102  CPPAD_ASSERT_UNKNOWN( p == 0 || p == q );
104  q <= num_order_taylor_ || p == 0,
105  "Forward(q, xq): Number of Taylor coefficient orders stored in this"
106  " ADFun\nis less than q and xq.size() != n*(q+1)."
107  );
109  p <= 1 || num_direction_taylor_ == 1,
110  "Forward(q, xq): computing order q >= 2"
111  " and number of directions is not one."
112  "\nMust use Forward(q, r, xq) for this case"
113  );
114  // does taylor_ need more orders or fewer directions
115  if( (cap_order_taylor_ <= q) | (num_direction_taylor_ != 1) )
116  { if( p == 0 )
117  { // no need to copy old values during capacity_order
118  num_order_taylor_ = 0;
119  }
120  else num_order_taylor_ = q;
121  size_t c = std::max(q + 1, cap_order_taylor_);
122  size_t r = 1;
123  capacity_order(c, r);
124  }
125  CPPAD_ASSERT_UNKNOWN( cap_order_taylor_ > q );
126  CPPAD_ASSERT_UNKNOWN( num_direction_taylor_ == 1 );
127
128  // short hand notation for order capacity
129  size_t C = cap_order_taylor_;
130
131  // The optimizer may skip a step that does not affect dependent variables.
132  // Initilaizing zero order coefficients avoids following valgrind warning:
133  // "Conditional jump or move depends on uninitialised value(s)".
134  for(j = 0; j < num_var_tape_; j++)
135  { for(k = p; k <= q; k++)
136  taylor_[C * j + k] = CppAD::numeric_limits<Base>::quiet_NaN();
137  }
138
139  // set Taylor coefficients for independent variables
140  for(j = 0; j < n; j++)
142
145
146  if( p == q )
147  taylor_[ C * ind_taddr_[j] + q] = xq[j];
148  else
149  { for(k = 0; k <= q; k++)
150  taylor_[ C * ind_taddr_[j] + k] = xq[ (q+1)*j + k];
151  }
152  }
153
154  // evaluate the derivatives
155  CPPAD_ASSERT_UNKNOWN( cskip_op_.size() == play_.num_op_rec() );
157  if( q == 0 )
158  { local::forward0sweep(&play_, s, true,
159  n, num_var_tape_, C,
161  compare_change_count_,
162  compare_change_number_,
163  compare_change_op_index_
164  );
165  }
166  else
167  { local::forward1sweep(&play_, s, true, p, q,
168  n, num_var_tape_, C,
170  compare_change_count_,
171  compare_change_number_,
172  compare_change_op_index_
173  );
174  }
175
176  // return Taylor coefficients for dependent variables
177  VectorBase yq;
178  if( p == q )
179  { yq.resize(m);
180  for(i = 0; i < m; i++)
182  yq[i] = taylor_[ C * dep_taddr_[i] + q];
183  }
184  }
185  else
186  { yq.resize(m * (q+1) );
187  for(i = 0; i < m; i++)
188  { for(k = 0; k <= q; k++)
189  yq[ (q+1) * i + k] =
190  taylor_[ C * dep_taddr_[i] + k ];
191  }
192  }
193 # ifndef NDEBUG
194  if( check_for_nan_ )
195  { bool ok = true;
196  size_t index = m;
197  if( p == 0 )
198  { for(i = 0; i < m; i++)
199  { // Visual Studio 2012, CppAD required in front of isnan ?
200  if( CppAD::isnan( yq[ (q+1) * i + 0 ] ) )
201  { ok = false;
202  if( index == m )
203  index = i;
204  }
205  }
206  }
207  if( ! ok )
208  { CPPAD_ASSERT_UNKNOWN( index < m );
209  //
211  for(j = 0; j < n; j++)
212  x0[j] = taylor_[ C * ind_taddr_[j] + 0 ];
213  std::string file_name;
214  put_check_for_nan(x0, file_name);
215  std::stringstream ss;
216  ss <<
217  "yq = f.Forward(q, xq): a zero order Taylor coefficient is nan.\n"
218  "Corresponding independent variables vector was written "
219  "to binary a file.\n"
220  "vector_size = " << n << "\n" <<
221  "file_name = " << file_name << "\n" <<
222  "index = " << index << "\n";
223  // ss.str() returns a string object with a copy of the current
224  // contents in the stream buffer.
225  std::string msg_str = ss.str();
226  // msg_str.c_str() returns a pointer to the c-string
227  // representation of the string object's value.
228  const char* msg_char_star = msg_str.c_str();
230  true,
231  __LINE__,
232  __FILE__,
233  "if( CppAD::isnan( yq[ (q+1) * index + 0 ] )",
234  msg_char_star
235  );
236  }
238  "with the value nan."
239  );
240  if( 0 < q )
241  { for(i = 0; i < m; i++)
242  { for(k = p; k <= q; k++)
243  { // Studio 2012, CppAD required in front of isnan ?
244  ok &= ! CppAD::isnan( yq[ (q+1-p)*i + k-p ] );
245  }
246  }
247  }
249  "yq = f.Forward(q, xq): has a non-zero order Taylor coefficient\n"
250  "with the value nan (but zero order coefficients are not nan)."
251  );
252  }
253 # endif
254
255  // now we have q + 1 taylor_ coefficient orders per variable
256  num_order_taylor_ = q + 1;
257
258  return yq;
259 }
260
261 /*!
262 One order, multiple directions, forward mode Taylor coefficieints.
263
264 \tparam Base
265 The type used during the forward mode computations; i.e., the corresponding
266 recording of operations used the type AD<Base>.
267
268 \tparam VectorBase
269 is a Simple Vector class with eleements of type Base.
270
271 \param q
272 is the order for this forward mode computation,
273 <code>q > 0</code>.
274 There must be at least <code>q</code> Taylor coefficients
275 per variable before this call.
276 After this call there will be <code>q+1</code>
277 Taylor coefficients per variable.
278
279 \param r
280 is the number of directions for this calculation.
281 If <code>q != 1</code>, \c r must be the same as in the previous
282 call to Forward where \c q was equal to one.
283
284 \param xq
285 contains Taylor coefficients for the independent variables.
286 The size of xq must either be <code>r*n</code>,
287 For <code>j = 0 , ... , n-1</code>,
288 <code>ell = 0, ... , r-1</code>,
289 <code>xq[ ( r*j + ell ]</code>
290 is the q-th order coefficient for the j-th independent variable
291 and the ell-th direction.
292
293 \return
294 contains Taylor coefficients for the dependent variables.
295 The size of the return value \c y is <code>r*m</code>.
296 For <code>i = 0, ... , m-1</code>,
297 <code>ell = 0, ... , r-1</code>,
298 <code>y[ r*i + ell ]</code>
299 is the q-th order coefficient for the i-th dependent variable
300 and the ell-th direction.
301
302 \par taylor_
303 The Taylor coefficients up to order <code>q-1</code> are inputs
304 and the coefficents of order \c q are outputs.
305 Let <code>N = num_var_tape_</code>, and
306 <code>C = cap_order_taylor_</code>.
307 Note that for
308 <code>i = 1 , ..., N-1</code>,
309 <code>taylor_[ (C-1)*r*i + i + 0 ]</code>
310 is the zero order cofficent,
311 for the i-th varaible, and all directions.
312 For <code>i = 1 , ..., N-1</code>,
313 <code>k = 1 , ..., q</code>,
314 <code>ell = 0 , ..., r-1</code>,
315 <code>taylor_[ (C-1)*r*i + i + (k-1)*r + ell + 1 ]</code>
316 is the k-th order cofficent,
317 for the i-th varaible, and ell-th direction.
318 (The first independent variable has index one on the tape
319 and there is no variable with index zero.)
320 */
321
322 template <typename Base>
323 template <typename VectorBase>
325  size_t q ,
326  size_t r ,
327  const VectorBase& xq )
328 { // temporary indices
329  size_t i, j, ell;
330
331  // number of independent variables
333
334  // number of dependent variables
336
337  // check Vector is Simple Vector class with Base type elements
338  CheckSimpleVector<Base, VectorBase>();
339
340  CPPAD_ASSERT_KNOWN( q > 0, "Forward(q, r, xq): q == 0" );
342  size_t(xq.size()) == r * n,
343  "Forward(q, r, xq): xq.size() is not equal r * n"
344  );
346  q <= num_order_taylor_ ,
347  "Forward(q, r, xq): Number of Taylor coefficient orders stored in"
348  " this ADFun is less than q"
349  );
351  q == 1 || num_direction_taylor_ == r ,
352  "Forward(q, r, xq): q > 1 and number of Taylor directions r"
353  " is not same as previous Forward(1, r, xq)"
354  );
355
356  // does taylor_ need more orders or new number of directions
357  if( cap_order_taylor_ <= q || num_direction_taylor_ != r )
358  { if( num_direction_taylor_ != r )
359  num_order_taylor_ = 1;
360
361  size_t c = std::max(q + 1, cap_order_taylor_);
362  capacity_order(c, r);
363  }
364  CPPAD_ASSERT_UNKNOWN( cap_order_taylor_ > q );
365  CPPAD_ASSERT_UNKNOWN( num_direction_taylor_ == r )
366
367  // short hand notation for order capacity
368  size_t c = cap_order_taylor_;
369
370  // set Taylor coefficients for independent variables
371  for(j = 0; j < n; j++)
373
376
377  for(ell = 0; ell < r; ell++)
378  { size_t index = ((c-1)*r + 1)*ind_taddr_[j] + (q-1)*r + ell + 1;
379  taylor_[ index ] = xq[ r * j + ell ];
380  }
381  }
382
383  // evaluate the derivatives
384  CPPAD_ASSERT_UNKNOWN( cskip_op_.size() == play_.num_op_rec() );
387  &play_,
388  q,
389  r,
390  n,
391  num_var_tape_,
392  c,
393  taylor_.data(),
394  cskip_op_.data(),
396  );
397
398  // return Taylor coefficients for dependent variables
399  VectorBase yq;
400  yq.resize(r * m);
401  for(i = 0; i < m; i++)
403  for(ell = 0; ell < r; ell++)
404  { size_t index = ((c-1)*r + 1)*dep_taddr_[i] + (q-1)*r + ell + 1;
405  yq[ r * i + ell ] = taylor_[ index ];
406  }
407  }
408 # ifndef NDEBUG
409  if( check_for_nan_ )
410  { bool ok = true;
411  for(i = 0; i < m; i++)
412  { for(ell = 0; ell < r; ell++)
413  { // Studio 2012, CppAD required in front of isnan ?
414  ok &= ! CppAD::isnan( yq[ r * i + ell ] );
415  }
416  }
418  "yq = f.Forward(q, r, xq): has a non-zero order Taylor coefficient\n"
419  "with the value nan (but zero order coefficients are not nan)."
420  );
421  }
422 # endif
423
424  // now we have q + 1 taylor_ coefficient orders per variable
425  num_order_taylor_ = q + 1;
426
427  return yq;
428 }
429
430
432 # endif
Check that exp is true, if not print msg and terminate execution.
void put_check_for_nan(const CppAD::vector< Base > &vec, std::string &file_name)
Stores a vector in a file when nans occur.
void forward0sweep(const local::player< Base > *play, std::ostream &s_out, bool print, size_t n, size_t numvar, size_t J, Base *taylor, bool *cskip_op, pod_vector< addr_t > &var_by_load_op, size_t compare_change_count, size_t &compare_change_number, size_t &compare_change_op_index)
Compute zero order forward mode Taylor coefficients.
bool isnan(const Scalar &s)
Definition: nan.hpp:169
static void Call(bool known, int line, const char *file, const char *exp, const char *msg)
All tthese defaults correspond to errors.