TMB Documentation  v1.9.11
mask.hpp
1 template<int base, int n>
2 struct pow_t {
3  static const long int value = base * pow_t<base, n-1>::value;
4 };
5 template<int base>
6 struct pow_t<base, 0> {
7  static const long int value = 1;
8 };
9 template<long int mask>
10 struct mask_t {
11  template<int length, int i=0>
12  struct set_length {
13  static const int base = 8;
14  static const long int power = pow_t<base, i>::value;
15  typedef set_length<length, i+1> next_index_t;
16  static const int flag = ( (mask / power) % base ) != 0;
17  static const int count = flag + next_index_t::count;
18  static const int Id = count - 1;
19  static const int Index = length - i - 1;
20  next_index_t next_index;
21  template<class S, class T>
22  void copy(S &dest, const T &orig) {
23  dest[Index] = (flag ? orig[Id] : 0);
24  next_index.copy(dest, orig);
25  }
26  template<class S, class T>
27  void activate_derivs(S &var, T &value) {
28  var[Index] = value[Index];
29  if (flag) var[Index].setid(Id);
30  next_index.activate_derivs(var, value);
31  }
32  };
33  template<int length>
34  struct set_length<length, length> {
35  static const int count = 0;
36  void trace() { }
37  template<class S, class T>
38  void copy(S &dest, const T &orig) { }
39  template<class S, class T>
40  void activate_derivs(S &var, T &value) { }
41  };
42 };
43 
44 template<int nvar>
45 struct ADTypes {
46  typedef tiny_ad::variable<1, nvar> order1;
47  typedef tiny_ad::variable<2, nvar> order2;
48  typedef tiny_ad::variable<3, nvar> order3;
49 };
50 
51 // 'TMB_BIND_ATOMIC' depends on these:
52 #define NCHAR(x) sizeof(#x)-1
53 #define OCTAL(x) 0 ## x
54 
73 #define TMB_BIND_ATOMIC(NAME,MASK,CALL) \
74 TMB_ATOMIC_VECTOR_FUNCTION( \
75  NAME \
76  , \
77  (size_t) \
78  pow((double) \
79  atomic::mask_t<OCTAL(MASK)>::set_length< NCHAR(MASK) >::count, \
80  CppAD::Integer(tx[NCHAR(MASK)])) \
81  , \
82  int order = CppAD::Integer(tx[NCHAR(MASK)]); \
83  typedef \
84  atomic::mask_t<OCTAL(MASK)>::set_length<NCHAR(MASK)> mask_type; \
85  mask_type mask; \
86  static const int nvar = mask_type::count; \
87  atomic::tiny_vec_ref<double> tyref(&ty[0], ty.size()); \
88  if(order==0) { \
89  typedef double Float; \
90  CppAD::vector<Float> x(tx); \
91  ty[0] = CALL; \
92  } \
93  else if (order==1) { \
94  typedef typename atomic::ADTypes<nvar>::order1 Float; \
95  Float x[NCHAR(MASK)]; \
96  mask.activate_derivs(x, tx); \
97  tyref = CALL.getDeriv(); \
98  } \
99  else if (order==2) { \
100  typedef typename atomic::ADTypes<nvar>::order2 Float; \
101  Float x[NCHAR(MASK)]; \
102  mask.activate_derivs(x, tx); \
103  tyref = CALL.getDeriv(); \
104  } \
105  else if (order==3) { \
106  typedef typename atomic::ADTypes<nvar>::order3 Float; \
107  Float x[NCHAR(MASK)]; \
108  mask.activate_derivs(x, tx); \
109  tyref = CALL.getDeriv(); \
110  } \
111  else \
112  Rf_error("Order not implemented"); \
113  , \
114  typedef \
115  atomic::mask_t<OCTAL(MASK)>::set_length<NCHAR(MASK)> mask_type; \
116  mask_type mask; \
117  static const int nvar = mask_type::count; \
118  CppAD::vector<Type> tx_(tx); \
119  tx_[NCHAR(MASK)] = tx_[NCHAR(MASK)] + Type(1.0); \
120  vector<Type> tmp = NAME(tx_); \
121  matrix<Type> m = tmp.matrix(); \
122  m.resize(nvar, m.size() / nvar); \
123  vector<Type> w = py; \
124  vector<Type> px_ = m * w.matrix(); \
125  mask.copy(px, px_); \
126  px[NCHAR(MASK)] = 0; \
127  )
128 
129 // ======================================================================
130 
131 
132 #ifdef TMBAD_FRAMEWORK
133 
134 #undef TMB_BIND_ATOMIC
135 #ifndef TMB_MAX_ORDER
136 #define TMB_MAX_ORDER 3
137 #endif
138 
139 #define TMB_BIND_ATOMIC(NAME,MASK,CALL) \
140 template<int order, int ninput, int noutput, long int mask> \
141 struct NAME ## Eval { \
142  typedef typename \
143  atomic::mask_t<mask>::template set_length<ninput> mask_type; \
144  mask_type mask_; \
145  static const int nvar = mask_type::count; \
146  template <class S, class T> \
147  void operator()(S* tx, T* ty) { \
148  typedef atomic::tiny_ad::variable<order, nvar> Float; \
149  atomic::tiny_vec_ref<double> tyref(&(ty[0]), noutput); \
150  Float x[ninput]; \
151  mask_.activate_derivs(x, tx); \
152  tyref = (CALL).getDeriv(); \
153  } \
154 }; \
155 template<int ninput, int noutput, long int mask> \
156 struct NAME ## Eval<0, ninput, noutput, mask> { \
157  template <class S, class T> \
158  void operator()(S* tx, T* ty) { \
159  S* x = tx; \
160  ty[0] = (CALL); \
161  } \
162 }; \
163 template<int order, int ninput, int noutput, long int mask> \
164 struct NAME ## Op : TMBad::global::Operator<ninput, noutput> { \
165  static const bool add_forward_replay_copy = true; \
166  typedef typename \
167  atomic::mask_t<mask>::template set_length<ninput> mask_type; \
168  mask_type mask_; \
169  static const int nvar = mask_type::count; \
170  template <class S, class T> \
171  void eval(S* tx, T* ty) const { \
172  NAME ## Eval<order, ninput, noutput, mask>()(tx, ty); \
173  } \
174  std::vector<TMBad::ad_plain> \
175  add_to_tape(const std::vector<TMBad::ad_plain> &x) { \
176  TMBad::OperatorPure* pOp = TMBad::get_glob()->getOperator<NAME ## Op>(); \
177  return \
178  TMBad::get_glob()->add_to_stack<NAME ## Op>(pOp, x); \
179  } \
180  std::vector<TMBad::ad_plain> \
181  operator()(const std::vector<TMBad::ad_plain> &x) { \
182  return add_to_tape(x); \
183  } \
184  Eigen::Matrix<TMBad::ad_aug, nvar, noutput / nvar> \
185  operator()(const Eigen::Array<TMBad::ad_aug, ninput, 1> &x) { \
186  std::vector<TMBad::ad_plain> x_(&(x(0)), &(x(0)) + x.size()); \
187  Eigen::Matrix<TMBad::ad_aug, nvar, noutput / nvar> ans; \
188  std::vector<TMBad::ad_plain> y = add_to_tape(x_); \
189  for (size_t i=0; i<y.size(); i++) ans(i) = y[i]; \
190  return ans; \
191  } \
192  Eigen::Matrix<double, nvar, noutput / nvar> \
193  operator()(const Eigen::Array<double, ninput, 1> &x) { \
194  Eigen::Matrix<double, nvar, noutput / nvar> ans; \
195  eval(&(x(0)), &(ans(0))); \
196  return ans; \
197  } \
198  template<class Type> \
199  void forward(TMBad::ForwardArgs<Type> &args) { \
200  Rf_error("Un-implemented method request"); \
201  } \
202  void forward(TMBad::ForwardArgs<double> &args) { \
203  double x[ninput]; \
204  for (size_t i=0; i<ninput; i++) x[i] = args.x(i); \
205  eval(x, &(args.y(0))); \
206  } \
207  template<class Type> \
208  void reverse(TMBad::ReverseArgs<Type> &args) { \
209  Eigen::Array<Type, ninput, 1> tx; \
210  for (size_t i=0; i<ninput; i++) tx(i) = args.x(i); \
211  Eigen::Matrix<Type, noutput, 1> w; \
212  for (size_t i=0; i<noutput; i++) w(i) = args.dy(i); \
213  NAME ## Op<order+1, ninput, noutput * nvar, mask> foo; \
214  Eigen::Matrix<Type, nvar, noutput> ty; \
215  ty = foo(tx); \
216  Eigen::Matrix<Type, nvar, 1> tyw = ty * w; \
217  Type tmp[ninput]; \
218  mask_.copy(tmp, &(tyw[0])); \
219  for (size_t i=0; i<ninput; i++) args.dx(i) += tmp[i]; \
220  } \
221  void reverse(TMBad::ReverseArgs<TMBad::Writer> &args) { \
222  Rf_error("Un-implemented method request"); \
223  } \
224  const char* op_name() { return #NAME ; } \
225 }; \
226 template<int ninput, int noutput, long int mask> \
227 struct NAME ## Op<TMB_MAX_ORDER+1, ninput, noutput, mask> { \
228  typedef typename \
229  atomic::mask_t<mask>::template set_length<ninput> mask_type; \
230  mask_type mask_; \
231  static const int nvar = mask_type::count; \
232  Eigen::Matrix<TMBad::ad_aug, nvar, noutput / nvar> \
233  operator()(const Eigen::Array<TMBad::ad_aug, ninput, 1> &x) { \
234  Eigen::Matrix<TMBad::ad_aug, nvar, noutput / nvar> ans; \
235  Rf_error("Order not implemented. Please increase TMB_MAX_ORDER"); \
236  return ans; \
237  } \
238  Eigen::Matrix<double, nvar, noutput / nvar> \
239  operator()(const Eigen::Array<double, ninput, 1> &x) { \
240  Eigen::Matrix<double, nvar, noutput / nvar> ans; \
241  Rf_error("Order not implemented. Please increase TMB_MAX_ORDER"); \
242  return ans; \
243  } \
244 }; \
245 template<class dummy=void> \
246 CppAD::vector<double> \
247 NAME (const CppAD::vector<double> &x) CSKIP_ATOMIC({ \
248  int n = x.size() - 1; \
249  int order = CppAD::Integer(x[n]); \
250  typedef NAME ## Op<0, NCHAR(MASK), 1, OCTAL(MASK)> Foo0; \
251  static const int nvar = Foo0::nvar; \
252  typedef NAME ## Op<1, NCHAR(MASK), nvar, OCTAL(MASK)> Foo1; \
253  if (order==0) { \
254  CppAD::vector<double> y(1); \
255  y[0] = CALL; \
256  return y; \
257  } \
258  else if (order==1) { \
259  Foo1 foo1; \
260  CppAD::vector<double> y(nvar); \
261  foo1.eval(&x[0], &y[0]); \
262  return y; \
263  } \
264  else { \
265  Rf_error("This interface is limited to 0th and 1st deriv order"); \
266  } \
267 }) \
268 template<class dummy=void> \
269 CppAD::vector<TMBad::ad_aug> \
270 NAME (const CppAD::vector<TMBad::ad_aug> &x) CSKIP_ATOMIC({ \
271  bool all_constant = true; \
272  for (size_t i = 0; i<x.size(); i++) \
273  all_constant &= x[i].constant(); \
274  if (all_constant) { \
275  CppAD::vector<double> xd(x.size()); \
276  for (size_t i=0; i<xd.size(); i++) xd[i] = x[i].Value(); \
277  CppAD::vector<double> yd = NAME(xd); \
278  CppAD::vector<TMBad::ad_aug> y(yd.size()); \
279  for (size_t i=0; i<yd.size(); i++) y[i] = yd[i]; \
280  return y; \
281  } \
282  int n = x.size() - 1; \
283  int order = CppAD::Integer(x[n]); \
284  std::vector<TMBad::ad_plain> x_(&(x[0]), &(x[0]) + n); \
285  std::vector<TMBad::ad_plain> y_; \
286  typedef NAME ## Op<0, NCHAR(MASK), 1, OCTAL(MASK)> Foo0; \
287  static const int nvar = Foo0::nvar; \
288  typedef NAME ## Op<1, NCHAR(MASK), nvar, OCTAL(MASK)> Foo1; \
289  if (order==0) { \
290  Foo0 foo0; \
291  y_ = foo0(x_); \
292  } \
293  else if (order==1) { \
294  Foo1 foo1; \
295  y_ = foo1(x_); \
296  } \
297  else { \
298  Rf_error("This interface is limited to 0th and 1st deriv order"); \
299  } \
300  CppAD::vector<TMBad::ad_aug> y(y_.size()); \
301  for (size_t i=0; i<y.size(); i++) y[i] = y_[i]; \
302  return y; \
303 }) \
304 IF_TMB_PRECOMPILE_ATOMICS( \
305 template \
306 CppAD::vector<TMBad::ad_aug> \
307 NAME<> (const CppAD::vector<TMBad::ad_aug> &x); \
308 template \
309 CppAD::vector<double> \
310 NAME<> (const CppAD::vector<double> &x); \
311 )
312 
313 #endif // TMBAD_FRAMEWORK
License: GPL v2