CppAD: A C++ Algorithmic Differentiation Package  20171217
reverse.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 # include <algorithm>
17
19 /*!
20 \file reverse.hpp
21 Compute derivatives using reverse mode.
22 */
23
24
25 /*!
26 Use reverse mode to compute derivative of forward mode Taylor coefficients.
27
28 The function
29 \f$X : {\bf R} \times {\bf R}^{n \times q} \rightarrow {\bf R} \f$
30 is defined by
31 \f[
32 X(t , u) = \sum_{k=0}^{q-1} u^{(k)} t^k
33 \f]
34 The function
35 \f$Y : {\bf R} \times {\bf R}^{n \times q} \rightarrow {\bf R} \f$
36 is defined by
37 \f[
38 Y(t , u) = F[ X(t, u) ]
39 \f]
40 The function
41 \f$W : {\bf R}^{n \times q} \rightarrow {\bf R} \f$ is defined by
42 \f[
43 W(u) = \sum_{k=0}^{q-1} ( w^{(k)} )^{\rm T}
44 \frac{1}{k !} \frac{ \partial^k } { t^k } Y(0, u)
45 \f]
46
47 \tparam Base
48 base type for the operator; i.e., this operation sequence was recorded
49 using AD< \a Base > and computations by this routine are done using type
50 \a Base.
51
52 \tparam VectorBase
53 is a Simple Vector class with elements of type \a Base.
54
55 \param q
56 is the number of the number of Taylor coefficients that are being
57 differentiated (per variable).
58
59 \param w
60 is the weighting for each of the Taylor coefficients corresponding
61 to dependent variables.
62 If the argument \a w has size <tt>m * q </tt>,
63 for \f$k = 0 , \ldots , q-1 \f$ and \f$i = 0, \ldots , m-1 \f$,
64 \f[
65  w_i^{(k)} = w [ i * q + k ]
66 \f]
67 If the argument \a w has size \c m ,
68 for \f$k = 0 , \ldots , q-1 \f$ and \f$i = 0, \ldots , m-1 \f$,
69 \f[
70 w_i^{(k)} = \left\{ \begin{array}{ll}
71  w [ i ] & {\rm if} \; k = q-1
72  \\
73  0 & {\rm otherwise}
74 \end{array} \right.
75 \f]
76
77 \return
78 Is a vector \f$dw \f$ such that
79 for \f$j = 0 , \ldots , n-1 \f$ and
80 \f$k = 0 , \ldots , q-1 \f$
81 \f[
82  dw[ j * q + k ] = W^{(1)} ( x )_{j,k}
83 \f]
84 where the matrix \f$x \f$ is the value for \f$u \f$
85 that corresponding to the forward mode Taylor coefficients
86 for the independent variables as specified by previous calls to Forward.
87
88 */
89 template <typename Base>
90 template <typename VectorBase>
91 VectorBase ADFun<Base>::Reverse(size_t q, const VectorBase &w)
92 { // constants
93  const Base zero(0);
94
95  // temporary indices
96  size_t i, j, k;
97
98  // number of independent variables
100
101  // number of dependent variables
103
104  // check VectorBase is Simple Vector class with Base type elements
105  CheckSimpleVector<Base, VectorBase>();
106
108  size_t(w.size()) == m || size_t(w.size()) == (m * q),
109  "Argument w to Reverse does not have length equal to\n"
110  "the dimension of the range or dimension of range times q."
111  );
113  q > 0,
114  "The first argument to Reverse must be greater than zero."
115  );
117  num_order_taylor_ >= q,
118  "Less than q Taylor coefficients are currently stored"
119  " in this ADFun object."
120  );
121  // special case where multiple forward directions have been computed,
122  // but we are only using the one direction zero order results
123  if( (q == 1) & (num_direction_taylor_ > 1) )
124  { num_order_taylor_ = 1; // number of orders to copy
125  size_t c = cap_order_taylor_; // keep the same capacity setting
126  size_t r = 1; // only keep one direction
127  capacity_order(c, r);
128  }
130  num_direction_taylor_ == 1,
131  "Reverse mode for Forward(q, r, xq) with more than one direction"
132  "\n(r > 1) is not yet supported for q > 1."
133  );
134
135  // initialize entire Partial matrix to zero
136  local::pod_vector<Base> Partial(num_var_tape_ * q);
137  for(i = 0; i < num_var_tape_; i++)
138  for(j = 0; j < q; j++)
139  Partial[i * q + j] = zero;
140
141  // set the dependent variable direction
142  // (use += because two dependent variables can point to same location)
143  for(i = 0; i < m; i++)
145  if( size_t(w.size()) == m )
146  Partial[dep_taddr_[i] * q + q - 1] += w[i];
147  else
148  { for(k = 0; k < q; k++)
149  // ? should use += here, first make test to demonstrate bug
150  Partial[ dep_taddr_[i] * q + k ] = w[i * q + k ];
151  }
152  }
153
154  // evaluate the derivatives
155  CPPAD_ASSERT_UNKNOWN( cskip_op_.size() == play_.num_op_rec() );
158  q - 1,
159  n,
160  num_var_tape_,
161  &play_,
162  cap_order_taylor_,
163  taylor_.data(),
164  q,
165  Partial.data(),
166  cskip_op_.data(),
168  subgraph_info_.entire_graph()
169  );
170
171  // return the derivative values
172  VectorBase value(n * q);
173  for(j = 0; j < n; j++)
175
178
179  // by the Reverse Identity Theorem
180  // partial of y^{(k)} w.r.t. u^{(0)} is equal to
181  // partial of y^{(q-1)} w.r.t. u^{(q - 1 - k)}
182  if( size_t(w.size()) == m )
183  { for(k = 0; k < q; k++)
184  value[j * q + k ] =
185  Partial[ind_taddr_[j] * q + q - 1 - k];
186  }
187  else
188  { for(k = 0; k < q; k++)
189  value[j * q + k ] =
190  Partial[ind_taddr_[j] * q + k];
191  }
192  }
193  CPPAD_ASSERT_KNOWN( ! ( hasnan(value) && check_for_nan_ ) ,
194  "dw = f.Reverse(q, w): has a nan,\n"
195  "but none of its Taylor coefficents are nan."
196  );
197
198  return value;
199 }
200
201
203 # endif
Check that exp is true, if not print msg and terminate execution.
File used to define pod_vector class.
bool hasnan(const Vector &v)
Definition: nan.hpp:174
Type * data(void)
current data pointer is no longer valid after any of the following: extend, erase, operator=, and ~pod_vector. Take extreem care when using this function.
Definition: pod_vector.hpp:89