CppAD: A C++ Algorithmic Differentiation Package  20171217
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
reverse.hpp
Go to the documentation of this file.
1 # ifndef CPPAD_CORE_REVERSE_HPP
2 # define CPPAD_CORE_REVERSE_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 # include <algorithm>
17 
18 namespace CppAD { // BEGIN_CPPAD_NAMESPACE
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
99  size_t n = ind_taddr_.size();
100 
101  // number of dependent variables
102  size_t m = dep_taddr_.size();
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++)
144  { CPPAD_ASSERT_UNKNOWN( dep_taddr_[i] < num_var_tape_ );
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() );
156  CPPAD_ASSERT_UNKNOWN( load_op_.size() == play_.num_load_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(),
167  load_op_,
168  subgraph_info_.entire_graph()
169  );
170 
171  // return the derivative values
172  VectorBase value(n * q);
173  for(j = 0; j < n; j++)
174  { CPPAD_ASSERT_UNKNOWN( ind_taddr_[j] < num_var_tape_ );
175 
176  // independent variable taddr equals its operator taddr
177  CPPAD_ASSERT_UNKNOWN( play_.GetOp( ind_taddr_[j] ) == local::InvOp );
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 
202 } // END_CPPAD_NAMESPACE
203 # endif
#define CPPAD_ASSERT_KNOWN(exp, msg)
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
#define CPPAD_ASSERT_UNKNOWN(exp)
Check that exp is true, if not terminate execution.
void reverse_sweep(size_t d, size_t n, size_t numvar, const local::player< Base > *play, size_t J, const Base *Taylor, size_t K, Base *Partial, bool *cskip_op, const pod_vector< addr_t > &var_by_load_op, const pod_vector< addr_t > &subgraph)
Compute derivative of arbitrary order forward mode Taylor coefficients.
VectorBase Reverse(size_t p, const VectorBase &v)
reverse mode sweep
Definition: reverse.hpp:91