NGSolve  5.3
autodiff.hpp
1 #ifndef FILE_AUTODIFF
2 #define FILE_AUTODIFF
3 
4 /**************************************************************************/
5 /* File: autodiff.hpp */
6 /* Author: Joachim Schoeberl */
7 /* Date: 24. Oct. 02 */
8 /**************************************************************************/
9 
10 namespace ngstd
11 {
12 
13 // Automatic differentiation datatype
14 
15  template <int D, typename SCAL = double> class AutoDiffRec;
16 
17 
23 template <int D, typename SCAL = double>
24 class AutoDiffVec // : public AlignedAlloc<AutoDiff<D,SCAL>>
25 {
26  SCAL val;
27  SCAL dval[D?D:1];
28 public:
29 
30  typedef AutoDiffVec<D,SCAL> TELEM;
31  typedef SCAL TSCAL;
32 
33 
35  // INLINE AutoDiffVec () throw() { };
36  AutoDiffVec() = default;
37  // { val = 0; for (int i = 0; i < D; i++) dval[i] = 0; } // !
38 
40  AutoDiffVec (const AutoDiffVec & ad2) = default;
41  /*
42  INLINE AutoDiffVec (const AutoDiffVec & ad2) throw()
43  {
44  val = ad2.val;
45  for (int i = 0; i < D; i++)
46  dval[i] = ad2.dval[i];
47  }
48  */
50  INLINE AutoDiffVec (SCAL aval) throw()
51  {
52  val = aval;
53  for (int i = 0; i < D; i++)
54  dval[i] = 0;
55  }
56 
58  INLINE AutoDiffVec (SCAL aval, int diffindex) throw()
59  {
60  val = aval;
61  for (int i = 0; i < D; i++)
62  dval[i] = 0;
63  dval[diffindex] = 1;
64  }
65 
66  INLINE AutoDiffVec (SCAL aval, const SCAL * grad)
67  {
68  val = aval;
69  LoadGradient (grad);
70  }
71 
73  INLINE AutoDiffVec & operator= (SCAL aval) throw()
74  {
75  val = aval;
76  for (int i = 0; i < D; i++)
77  dval[i] = 0;
78  return *this;
79  }
80 
81  AutoDiffVec & operator= (const AutoDiffVec & ad2) = default;
82 
84  INLINE SCAL Value() const throw() { return val; }
85 
87  INLINE SCAL DValue (int i) const throw() { return dval[i]; }
88 
90  INLINE void StoreGradient (SCAL * p) const
91  {
92  for (int i = 0; i < D; i++)
93  p[i] = dval[i];
94  }
95 
96  INLINE void LoadGradient (const SCAL * p)
97  {
98  for (int i = 0; i < D; i++)
99  dval[i] = p[i];
100  }
101 
103  INLINE SCAL & Value() throw() { return val; }
104 
106  INLINE SCAL & DValue (int i) throw() { return dval[i]; }
107 };
108 
109 
111 
113 template<int D, typename SCAL>
114 inline ostream & operator<< (ostream & ost, const AutoDiffVec<D,SCAL> & x)
115 {
116  ost << x.Value() << ", D = ";
117  for (int i = 0; i < D; i++)
118  ost << x.DValue(i) << " ";
119  return ost;
120 }
121 
123 template<int D, typename SCAL>
124 INLINE AutoDiffVec<D,SCAL> operator+ (const AutoDiffVec<D,SCAL> & x, const AutoDiffVec<D,SCAL> & y) throw()
125 {
126  AutoDiffVec<D,SCAL> res;
127  res.Value () = x.Value()+y.Value();
128  // AutoDiffVec<D,SCAL> res(x.Value()+y.Value());
129  for (int i = 0; i < D; i++)
130  res.DValue(i) = x.DValue(i) + y.DValue(i);
131  return res;
132 }
133 
134 
136 template<int D, typename SCAL>
137 INLINE AutoDiffVec<D,SCAL> operator- (const AutoDiffVec<D,SCAL> & x, const AutoDiffVec<D,SCAL> & y) throw()
138 {
139  AutoDiffVec<D,SCAL> res;
140  res.Value() = x.Value()-y.Value();
141  // AutoDiffVec<D,SCAL> res (x.Value()-y.Value());
142  for (int i = 0; i < D; i++)
143  res.DValue(i) = x.DValue(i) - y.DValue(i);
144  return res;
145 }
146 
148  template<int D, typename SCAL, typename SCAL2,
149  typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
150 INLINE AutoDiffVec<D,SCAL> operator+ (SCAL2 x, const AutoDiffVec<D,SCAL> & y) throw()
151 {
152  AutoDiffVec<D,SCAL> res;
153  res.Value() = x+y.Value();
154  for (int i = 0; i < D; i++)
155  res.DValue(i) = y.DValue(i);
156  return res;
157 }
158 
160  template<int D, typename SCAL, typename SCAL2,
161  typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
162 INLINE AutoDiffVec<D,SCAL> operator+ (const AutoDiffVec<D,SCAL> & y, SCAL2 x) throw()
163 {
164  AutoDiffVec<D,SCAL> res;
165  res.Value() = x+y.Value();
166  for (int i = 0; i < D; i++)
167  res.DValue(i) = y.DValue(i);
168  return res;
169 }
170 
171 
173 template<int D, typename SCAL>
174 INLINE AutoDiffVec<D,SCAL> operator- (const AutoDiffVec<D,SCAL> & x) throw()
175 {
176  AutoDiffVec<D,SCAL> res;
177  res.Value() = -x.Value();
178  for (int i = 0; i < D; i++)
179  res.DValue(i) = -x.DValue(i);
180  return res;
181 }
182 
184  template<int D, typename SCAL, typename SCAL2,
185  typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
186 INLINE AutoDiffVec<D,SCAL> operator- (const AutoDiffVec<D,SCAL> & x, SCAL2 y) throw()
187 {
188  AutoDiffVec<D,SCAL> res;
189  res.Value() = x.Value()-y;
190  for (int i = 0; i < D; i++)
191  res.DValue(i) = x.DValue(i);
192  return res;
193 }
194 
196  template<int D, typename SCAL, typename SCAL2,
197  typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
198 INLINE AutoDiffVec<D,SCAL> operator- (SCAL2 x, const AutoDiffVec<D,SCAL> & y) throw()
199 {
200  AutoDiffVec<D,SCAL> res;
201  res.Value() = x-y.Value();
202  for (int i = 0; i < D; i++)
203  res.DValue(i) = -y.DValue(i);
204  return res;
205 }
206 
207 
209  template<int D, typename SCAL, typename SCAL2,
210  typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
211 INLINE AutoDiffVec<D,SCAL> operator* (SCAL2 x, const AutoDiffVec<D,SCAL> & y) throw()
212 {
213  AutoDiffVec<D,SCAL> res;
214  res.Value() = x*y.Value();
215  for (int i = 0; i < D; i++)
216  res.DValue(i) = x*y.DValue(i);
217  return res;
218 }
219 
221  template<int D, typename SCAL, typename SCAL2,
222  typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
223 
224  INLINE AutoDiffVec<D,SCAL> operator* (const AutoDiffVec<D,SCAL> & y, SCAL2 x) throw()
225 {
226  AutoDiffVec<D,SCAL> res;
227  res.Value() = x*y.Value();
228  for (int i = 0; i < D; i++)
229  res.DValue(i) = x*y.DValue(i);
230  return res;
231 }
232 
234 template<int D, typename SCAL>
235 INLINE AutoDiffVec<D,SCAL> operator* (const AutoDiffVec<D,SCAL> & x, const AutoDiffVec<D,SCAL> & y) throw()
236 {
237  AutoDiffVec<D,SCAL> res;
238  SCAL hx = x.Value();
239  SCAL hy = y.Value();
240 
241  res.Value() = hx*hy;
242  for (int i = 0; i < D; i++)
243  res.DValue(i) = hx*y.DValue(i) + hy*x.DValue(i);
244 
245  return res;
246 }
247 
249 template<int D, typename SCAL>
250 INLINE AutoDiffVec<D,SCAL> sqr (const AutoDiffVec<D,SCAL> & x) throw()
251 {
252  AutoDiffVec<D,SCAL> res;
253  SCAL hx = x.Value();
254  res.Value() = hx*hx;
255  hx *= 2;
256  for (int i = 0; i < D; i++)
257  res.DValue(i) = hx*x.DValue(i);
258  return res;
259 }
260 
262 template<int D, typename SCAL>
263 INLINE AutoDiffVec<D,SCAL> Inv (const AutoDiffVec<D,SCAL> & x)
264 {
265  AutoDiffVec<D,SCAL> res(1.0 / x.Value());
266  for (int i = 0; i < D; i++)
267  res.DValue(i) = -x.DValue(i) / (x.Value() * x.Value());
268  return res;
269 }
270 
271 
273 template<int D, typename SCAL>
274 INLINE AutoDiffVec<D,SCAL> operator/ (const AutoDiffVec<D,SCAL> & x, const AutoDiffVec<D,SCAL> & y)
275 {
276  return x * Inv (y);
277 }
278 
280 template<int D, typename SCAL, typename SCAL2,
281  typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
282 INLINE AutoDiffVec<D,SCAL> operator/ (const AutoDiffVec<D,SCAL> & x, SCAL2 y)
283 {
284  return (1.0/y) * x;
285 }
286 
288 template<int D, typename SCAL, typename SCAL2,
289  typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
290  INLINE AutoDiffVec<D,SCAL> operator/ (SCAL2 x, const AutoDiffVec<D,SCAL> & y)
291  {
292  return x * Inv(y);
293  }
294 
295 
296 
297 
298  template <int D, typename SCAL, typename SCAL2>
299  INLINE AutoDiffVec<D,SCAL> & operator+= (AutoDiffVec<D,SCAL> & x, SCAL2 y) throw()
300  {
301  x.Value() += y;
302  return x;
303  }
304 
305 
307  template <int D, typename SCAL>
308  INLINE AutoDiffVec<D,SCAL> & operator+= (AutoDiffVec<D,SCAL> & x, AutoDiffVec<D,SCAL> y)
309  {
310  x.Value() += y.Value();
311  for (int i = 0; i < D; i++)
312  x.DValue(i) += y.DValue(i);
313  return x;
314  }
315 
317  template <int D, typename SCAL>
318  INLINE AutoDiffVec<D,SCAL> & operator-= (AutoDiffVec<D,SCAL> & x, AutoDiffVec<D,SCAL> y)
319  {
320  x.Value() -= y.Value();
321  for (int i = 0; i < D; i++)
322  x.DValue(i) -= y.DValue(i);
323  return x;
324 
325  }
326 
327  template <int D, typename SCAL, typename SCAL2>
328  INLINE AutoDiffVec<D,SCAL> & operator-= (AutoDiffVec<D,SCAL> & x, SCAL2 y)
329  {
330  x.Value() -= y;
331  return x;
332  }
333 
335  template <int D, typename SCAL>
336  INLINE AutoDiffVec<D,SCAL> & operator*= (AutoDiffVec<D,SCAL> & x, AutoDiffVec<D,SCAL> y)
337  {
338  for (int i = 0; i < D; i++)
339  x.DValue(i) = x.DValue(i)*y.Value() + x.Value() * y.DValue(i);
340  x.Value() *= y.Value();
341  return x;
342  }
343 
345  template <int D, typename SCAL, typename SCAL2>
346  INLINE AutoDiffVec<D,SCAL> & operator*= (AutoDiffVec<D,SCAL> & x, SCAL2 y)
347  {
348  x.Value() *= y;
349  for (int i = 0; i < D; i++)
350  x.DValue(i) *= y;
351  return x;
352  }
353 
355  template <int D, typename SCAL>
356  INLINE AutoDiffVec<D,SCAL> & operator/= (AutoDiffVec<D,SCAL> & x, SCAL y)
357  {
358  SCAL iy = 1.0 / y;
359  x.Value() *= iy;
360  for (int i = 0; i < D; i++)
361  x.DValue(i) *= iy;
362  return x;
363  }
364 
365 
366 
367 
369  template <int D, typename SCAL>
370  INLINE bool operator== (AutoDiffVec<D,SCAL> x, SCAL val2)
371  {
372  return x.Value() == val2;
373  }
374 
376  template <int D, typename SCAL>
377  INLINE bool operator!= (AutoDiffVec<D,SCAL> x, SCAL val2) throw()
378  {
379  return x.Value() != val2;
380  }
381 
383  template <int D, typename SCAL>
384  INLINE bool operator< (AutoDiffVec<D,SCAL> x, SCAL val2) throw()
385  {
386  return x.Value() < val2;
387  }
388 
390  template <int D, typename SCAL>
391  INLINE bool operator> (AutoDiffVec<D,SCAL> x, SCAL val2) throw()
392  {
393  return x.Value() > val2;
394  }
395 
396 
397 
398 
399 template<int D, typename SCAL>
400 INLINE AutoDiffVec<D,SCAL> fabs (const AutoDiffVec<D,SCAL> & x)
401 {
402  double abs = fabs (x.Value());
403  AutoDiffVec<D,SCAL> res( abs );
404  if (abs != 0.0)
405  for (int i = 0; i < D; i++)
406  res.DValue(i) = x.Value()*x.DValue(i) / abs;
407  else
408  for (int i = 0; i < D; i++)
409  res.DValue(i) = 0.0;
410  return res;
411 }
412 
413 using std::sqrt;
414 template<int D, typename SCAL>
415 INLINE AutoDiffVec<D,SCAL> sqrt (const AutoDiffVec<D,SCAL> & x)
416 {
417  AutoDiffVec<D,SCAL> res;
418  res.Value() = sqrt(x.Value());
419  for (int j = 0; j < D; j++)
420  res.DValue(j) = 0.5 / res.Value() * x.DValue(j);
421  return res;
422 }
423 
424 using std::log;
425 template <int D, typename SCAL>
426 AutoDiffVec<D,SCAL> log (AutoDiffVec<D,SCAL> x)
427 {
428  AutoDiffVec<D,SCAL> res;
429  res.Value() = log(x.Value());
430  for (int k = 0; k < D; k++)
431  res.DValue(k) = x.DValue(k) / x.Value();
432  return res;
433 }
434 
435 using std::exp;
436 template <int D, typename SCAL>
437 INLINE AutoDiffVec<D,SCAL> exp (AutoDiffVec<D,SCAL> x)
438 {
439  AutoDiffVec<D,SCAL> res;
440  res.Value() = exp(x.Value());
441  for (int k = 0; k < D; k++)
442  res.DValue(k) = x.DValue(k) * res.Value();
443  return res;
444 }
445 
446 using std::pow;
447 template <int D, typename SCAL>
448 INLINE AutoDiffVec<D,SCAL> pow (AutoDiffVec<D,SCAL> x, AutoDiffVec<D,SCAL> y )
449 {
450  return exp(log(x)*y);
451 }
452 
453 using std::sin;
454 /*
455 template <int D, typename SCAL>
456 INLINE AutoDiffVec<D,SCAL> sin (AutoDiffVec<D,SCAL> x)
457 {
458  AutoDiffVec<D,SCAL> res;
459  res.Value() = sin(x.Value());
460  SCAL c = cos(x.Value());
461  for (int k = 0; k < D; k++)
462  res.DValue(k) = x.DValue(k) * c;
463  return res;
464 }
465 */
466 
467 template <int D, typename SCAL>
468 INLINE AutoDiffVec<D,SCAL> sin (AutoDiffVec<D,SCAL> x)
469 {
470  return sin(AutoDiffRec<D,SCAL>(x));
471 }
472 
473 using std::cos;
474 /*
475 template <int D, typename SCAL>
476 INLINE AutoDiffVec<D,SCAL> cos (AutoDiffVec<D,SCAL> x)
477 {
478  AutoDiffVec<D,SCAL> res;
479  res.Value() = cos(x.Value());
480  SCAL ms = -sin(x.Value());
481  for (int k = 0; k < D; k++)
482  res.DValue(k) = x.DValue(k) * ms;
483  return res;
484 }
485 */
486 template <int D, typename SCAL>
487 INLINE AutoDiffVec<D,SCAL> cos (AutoDiffVec<D,SCAL> x)
488 {
489  return cos(AutoDiffRec<D,SCAL>(x));
490 }
491 
492 using std::tan;
493 template <int D, typename SCAL>
494 INLINE AutoDiffVec<D,SCAL> tan (AutoDiffVec<D,SCAL> x)
495 { return sin(x) / cos(x); }
496 
497 using std::sinh;
498 template <int D, typename SCAL>
499 INLINE AutoDiffVec<D,SCAL> sinh (AutoDiffVec<D,SCAL> x)
500 {
501  AutoDiffVec<D,SCAL> res;
502  res.Value() = sinh(x.Value());
503  SCAL ch = cosh(x.Value());
504  for (int k = 0; k < D; k++)
505  res.DValue(k) = x.DValue(k) * ch;
506  return res;
507 }
508 
509 using std::cosh;
510 template <int D, typename SCAL>
511 INLINE AutoDiffVec<D,SCAL> cosh (AutoDiffVec<D,SCAL> x)
512 {
513  AutoDiffVec<D,SCAL> res;
514  res.Value() = cosh(x.Value());
515  SCAL sh = sinh(x.Value());
516  for (int k = 0; k < D; k++)
517  res.DValue(k) = x.DValue(k) * sh;
518  return res;
519 }
520 
521 using std::floor;
522 template<int D, typename SCAL>
523 INLINE AutoDiffVec<D,SCAL> floor (const AutoDiffVec<D,SCAL> & x)
524 {
525  AutoDiffVec<D,SCAL> res;
526  res.Value() = floor(x.Value());
527  for (int j = 0; j < D; j++)
528  res.DValue(j) = 0.0;
529  return res;
530 }
531 
532 using std::ceil;
533 template<int D, typename SCAL>
534 INLINE AutoDiffVec<D,SCAL> ceil (const AutoDiffVec<D,SCAL> & x)
535 {
536  AutoDiffVec<D,SCAL> res;
537  res.Value() = ceil(x.Value());
538  for (int j = 0; j < D; j++)
539  res.DValue(j) = 0.0;
540  return res;
541 }
542 
543 
544 using std::atan;
545 /*
546 template <int D, typename SCAL>
547 INLINE AutoDiffVec<D,SCAL> atan (AutoDiffVec<D,SCAL> x)
548 {
549  AutoDiffVec<D,SCAL> res;
550  SCAL a = atan(x.Value());
551  res.Value() = a;
552  for (int k = 0; k < D; k++)
553  res.DValue(k) = x.DValue(k)/(1+x.Value()*x.Value()) ;
554  return res;
555 }
556 */
557 template <int D, typename SCAL>
558 AutoDiffVec<D,SCAL> atan (AutoDiffVec<D,SCAL> x)
559 {
560  return atan (AutoDiffRec<D,SCAL> (x));
561 }
562 
563 using std::atan2;
564 template <int D, typename SCAL>
565 INLINE AutoDiffVec<D,SCAL> atan2 (AutoDiffVec<D,SCAL> x, AutoDiffVec<D,SCAL> y)
566 {
567  AutoDiffVec<D,SCAL> res;
568  SCAL a = atan2(x.Value(), y.Value());
569  res.Value() = a;
570  for (int k = 0; k < D; k++)
571  res.DValue(k) = (x.Value()*y.DValue(k)-y.Value()*x.DValue(k))/(y.Value()*y.Value()+x.Value()*x.Value());
572  return res;
573 }
574 
575 
576 using std::acos;
577 template <int D, typename SCAL>
578 INLINE AutoDiffVec<D,SCAL> acos (AutoDiffVec<D,SCAL> x)
579 {
580  AutoDiffVec<D,SCAL> res;
581  SCAL a = acos(x.Value());
582  res.Value() = a;
583  SCAL da = -1 / sqrt(1-x.Value()*x.Value());
584  for (int k = 0; k < D; k++)
585  res.DValue(k) = x.DValue(k)*da;
586  return res;
587 }
588 
589 
590 using std::asin;
591 template <int D, typename SCAL>
592 INLINE AutoDiffVec<D,SCAL> asin (AutoDiffVec<D,SCAL> x)
593 {
594  AutoDiffVec<D,SCAL> res;
595  SCAL a = asin(x.Value());
596  res.Value() = a;
597  SCAL da = 1 / sqrt(1-x.Value()*x.Value());
598  for (int k = 0; k < D; k++)
599  res.DValue(k) = x.DValue(k)*da;
600  return res;
601 }
602 
603 
604 
605 
606  template <int D, typename SCAL, typename TB, typename TC>
607  auto IfPos (AutoDiffVec<D,SCAL> a, TB b, TC c) // -> decltype(IfPos (a.Value(), b, c))
608  {
609  return IfPos (a.Value(), b, c);
610  }
611 
612  template <int D, typename SCAL>
613  INLINE AutoDiffVec<D,SCAL> IfPos (SCAL /* SIMD<double> */ a, AutoDiffVec<D,SCAL> b, AutoDiffVec<D,SCAL> c)
614  {
615  AutoDiffVec<D,SCAL> res;
616  res.Value() = IfPos (a, b.Value(), c.Value());
617  for (int j = 0; j < D; j++)
618  res.DValue(j) = IfPos (a, b.DValue(j), c.DValue(j));
619  return res;
620  }
621 
622  template <int D, typename SCAL, typename TC>
623  INLINE AutoDiffVec<D,SCAL> IfPos (SCAL /* SIMD<double> */ a, AutoDiffVec<D,SCAL> b, TC c)
624  {
625  return IfPos (a, b, AutoDiffVec<D,SCAL> (c));
626  }
627 
629 
630 
631 
632  template <int D, typename SCAL>
633  class AutoDiffRec // : public AlignedAlloc<AutoDiffRec<D,SCAL>>
634  {
635  AutoDiffRec<D-1, SCAL> rec;
636  SCAL last;
637 
638  public:
639  INLINE AutoDiffRec () = default;
640  INLINE AutoDiffRec (const AutoDiffRec &) = default;
641  INLINE AutoDiffRec (AutoDiffRec<D-1,SCAL> _rec, SCAL _last) : rec(_rec), last(_last) { ; }
642  INLINE AutoDiffRec & operator= (const AutoDiffRec &) = default;
643 
644  INLINE AutoDiffRec (SCAL aval) : rec(aval), last(0.0) { ; }
645  INLINE AutoDiffRec (SCAL aval, int diffindex) : rec(aval, diffindex), last((diffindex==D-1) ? 1.0 : 0.0) { ; }
646  INLINE AutoDiffRec (SCAL aval, const SCAL * grad)
647  : rec(aval, grad), last(grad[D-1]) { }
648 
649  INLINE AutoDiffRec (const AutoDiffVec<D,SCAL> & ad)
650  {
651  Value() = ad.Value();
652  for (int i = 0; i < D; i++)
653  DValue(i) = ad.DValue(i);
654  }
655 
656  INLINE AutoDiffRec & operator= (SCAL aval) { rec = aval; last = 0.0; return *this; }
657  INLINE SCAL Value() const { return rec.Value(); }
658  INLINE SCAL DValue(int i) const { return (i == D-1) ? last : rec.DValue(i); }
659  INLINE SCAL & Value() { return rec.Value(); }
660  INLINE SCAL & DValue(int i) { return (i == D-1) ? last : rec.DValue(i); }
661  INLINE auto Rec() const { return rec; }
662  INLINE auto Last() const { return last; }
663  INLINE auto & Rec() { return rec; }
664  INLINE auto & Last() { return last; }
665  INLINE operator AutoDiffVec<D,SCAL> () const
666  {
667  AutoDiffVec<D,SCAL> res(Value());
668  for (int i = 0; i < D; i++)
669  res.DValue(i) = DValue(i);
670  return res;
671  }
672  };
673 
674  template<int D, typename SCAL>
675  ostream & operator<< (ostream & ost, AutoDiffRec<D,SCAL> ad)
676  {
677  return ost << AutoDiffVec<D,SCAL> (ad);
678  }
679 
680  template <typename SCAL>
681  class AutoDiffRec<0,SCAL> : public AlignedAlloc<AutoDiffRec<0,SCAL>>
682  {
683  SCAL val;
684  public:
685  INLINE AutoDiffRec () = default;
686  INLINE AutoDiffRec (const AutoDiffRec &) = default;
687  INLINE AutoDiffRec (SCAL _val) : val(_val) { ; }
688  INLINE AutoDiffRec (SCAL _val, SCAL /* _dummylast */) : val(_val) { ; }
689  INLINE AutoDiffRec (SCAL aval, const SCAL * /* grad */)
690  : val(aval) { }
691 
692  INLINE AutoDiffRec & operator= (const AutoDiffRec &) = default;
693  INLINE AutoDiffRec & operator= (SCAL aval) { val = aval; return *this; }
694 
695  INLINE SCAL Value() const { return val; }
696  INLINE SCAL DValue(int i) const { return SCAL(0); }
697  INLINE SCAL & Value() { return val; }
698  // SCAL & DValue(int i) { return val; }
699  INLINE auto Rec() const { return val; }
700  INLINE auto Last() const { return SCAL(0); }
701  INLINE auto & Rec() { return val; }
702  INLINE auto & Last() { return val; }
703  INLINE operator AutoDiffVec<0,SCAL> () const { return AutoDiffVec<0,SCAL>(); }
704  };
705 
706 
707  template <typename SCAL>
708  class AutoDiffRec<1,SCAL> : public AlignedAlloc<AutoDiffRec<1,SCAL>>
709  {
710  SCAL val;
711  SCAL last;
712  public:
713  INLINE AutoDiffRec () = default;
714  INLINE AutoDiffRec (const AutoDiffRec &) = default;
715  INLINE AutoDiffRec (SCAL _val) : val(_val), last(0.0) { ; }
716  INLINE AutoDiffRec (SCAL _val, SCAL _last) : val(_val), last(_last) { ; }
717  INLINE AutoDiffRec (SCAL aval, int diffindex) : val(aval), last((diffindex==0) ? 1.0 : 0.0) { ; }
718  INLINE AutoDiffRec (SCAL aval, const SCAL * grad)
719  : val(aval), last(grad[0]) { }
720 
721  INLINE AutoDiffRec (const AutoDiffVec<1,SCAL> & ad)
722  {
723  Value() = ad.Value();
724  DValue(0) = ad.DValue(0);
725  }
726 
727  INLINE AutoDiffRec & operator= (const AutoDiffRec &) = default;
728  INLINE AutoDiffRec & operator= (SCAL aval) { val = aval; last = 0.0; return *this; }
729 
730  INLINE SCAL Value() const { return val; }
731  INLINE SCAL DValue(int i) const { return last; }
732  INLINE SCAL & Value() { return val; }
733  INLINE SCAL & DValue(int i) { return last; }
734  INLINE auto Rec() const { return val; }
735  INLINE auto Last() const { return last; }
736  INLINE auto & Rec() { return val; }
737  INLINE auto & Last() { return last; }
738 
739  INLINE operator AutoDiffVec<1,SCAL> () const
740  {
741  AutoDiffVec<1,SCAL> res(Value());
742  res.DValue(0) = DValue(0);
743  return res;
744  }
745  };
746 
747  template <int D, typename SCAL, typename SCAL2,
748  typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
749  INLINE AutoDiffRec<D,SCAL> operator+ (SCAL2 a, AutoDiffRec<D,SCAL> b)
750  {
751  return AutoDiffRec<D,SCAL> (a+b.Rec(), b.Last());
752  }
753 
754  template <int D, typename SCAL, typename SCAL2,
755  typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
756  INLINE AutoDiffRec<D,SCAL> operator+ (AutoDiffRec<D,SCAL> a, SCAL2 b)
757  {
758  return AutoDiffRec<D,SCAL> (a.Rec()+b, a.Last());
759  }
760 
761  template <int D, typename SCAL>
762  INLINE AutoDiffRec<D,SCAL> operator+ (AutoDiffRec<D,SCAL> a, AutoDiffRec<D,SCAL> b)
763  {
764  return AutoDiffRec<D,SCAL> (a.Rec()+b.Rec(), a.Last()+b.Last());
765  }
766 
767  template <int D, typename SCAL, typename SCAL2,
768  typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
769  INLINE AutoDiffRec<D,SCAL> operator- (SCAL2 b, AutoDiffRec<D,SCAL> a)
770  {
771  return AutoDiffRec<D,SCAL> (b-a.Rec(), -a.Last());
772  }
773 
774  template <int D, typename SCAL, typename SCAL2,
775  typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
776  INLINE AutoDiffRec<D,SCAL> operator- (AutoDiffRec<D,SCAL> a, SCAL2 b)
777  {
778  return AutoDiffRec<D,SCAL> (a.Rec()-b, a.Last());
779  }
780 
781  template <int D, typename SCAL>
782  INLINE AutoDiffRec<D,SCAL> operator- (AutoDiffRec<D,SCAL> a, AutoDiffRec<D,SCAL> b)
783  {
784  return AutoDiffRec<D,SCAL> (a.Rec()-b.Rec(), a.Last()-b.Last());
785  }
786 
788  template<int D, typename SCAL>
789  INLINE AutoDiffRec<D,SCAL> operator- (AutoDiffRec<D,SCAL> a)
790  {
791  return AutoDiffRec<D,SCAL> (-a.Rec(), -a.Last());
792  }
793 
794  template <int D, typename SCAL>
795  INLINE AutoDiffRec<D,SCAL> operator* (AutoDiffRec<D,SCAL> a, AutoDiffRec<D,SCAL> b)
796  {
797  return AutoDiffRec<D,SCAL> (a.Rec()*b.Rec(), a.Value()*b.Last()+b.Value()*a.Last());
798  }
799 
800  template <int D, typename SCAL, typename SCAL1,
801  typename std::enable_if<std::is_convertible<SCAL1,SCAL>::value, int>::type = 0>
802  INLINE AutoDiffRec<D,SCAL> operator* (AutoDiffRec<D,SCAL> b, SCAL1 a)
803  {
804  return AutoDiffRec<D,SCAL> (a*b.Rec(), a*b.Last());
805  }
806 
807  template <int D, typename SCAL, typename SCAL1,
808  typename std::enable_if<std::is_convertible<SCAL1,SCAL>::value, int>::type = 0>
809  INLINE AutoDiffRec<D,SCAL> operator* (SCAL1 a, AutoDiffRec<D,SCAL> b)
810  {
811  return AutoDiffRec<D,SCAL> (a*b.Rec(), a*b.Last());
812  }
813 
814  template <int D, typename SCAL>
815  INLINE AutoDiffRec<D,SCAL> & operator+= (AutoDiffRec<D,SCAL> & a, AutoDiffRec<D,SCAL> b)
816  {
817  a.Rec() += b.Rec();
818  a.Last() += b.Last();
819  return a;
820  }
821 
822  template <int D, typename SCAL>
823  INLINE AutoDiffRec<D,SCAL> & operator-= (AutoDiffRec<D,SCAL> & a, double b)
824  {
825  a.Rec() -= b;
826  return a;
827  }
828 
829  template <int D, typename SCAL>
830  INLINE AutoDiffRec<D,SCAL> & operator-= (AutoDiffRec<D,SCAL> & a, AutoDiffRec<D,SCAL> b)
831  {
832  a.Rec() -= b.Rec();
833  a.Last() -= b.Last();
834  return a;
835  }
836 
837 
838  template <int D, typename SCAL>
839  INLINE AutoDiffRec<D,SCAL> & operator*= (AutoDiffRec<D,SCAL> & a, AutoDiffRec<D,SCAL> b)
840  {
841  a = a*b;
842  return a;
843  }
844 
845 
846  template <int D, typename SCAL, typename SCAL2>
847  INLINE AutoDiffRec<D,SCAL> & operator*= (AutoDiffRec<D,SCAL> & b, SCAL2 a)
848  {
849  b.Rec() *= a;
850  b.Last() *= a;
851  return b;
852  }
853 
855 
856  template <typename SCAL>
857  auto Inv1 (SCAL x) { return 1.0/x; }
858 
859  template<int D, typename SCAL>
860  INLINE AutoDiffRec<D,SCAL> Inv1 (AutoDiffRec<D,SCAL> x)
861  {
862  return AutoDiffRec<D,SCAL> (Inv1(x.Rec()), (-sqr(1.0/x.Value())) * x.Last());
863  }
864 
866  template<int D, typename SCAL>
867  INLINE AutoDiffRec<D,SCAL> operator/ (const AutoDiffRec<D,SCAL> & x, const AutoDiffRec<D,SCAL> & y)
868  {
869  return x * Inv1 (y);
870  }
871 
872 
874  template<int D, typename SCAL, typename SCAL2,
875  typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
876  INLINE AutoDiffRec<D,SCAL> operator/ (const AutoDiffRec<D,SCAL> & x, SCAL2 y)
877  {
878  return (1.0/y) * x;
879  }
880 
881 
883  template<int D, typename SCAL, typename SCAL2,
884  typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
885  INLINE AutoDiffRec<D,SCAL> operator/ (SCAL2 x, const AutoDiffRec<D,SCAL> & y)
886  {
887  return x * Inv1(y);
888  }
889 
890 
891 
892 
893 
894 
895 
897  template <int D, typename SCAL>
898  INLINE bool operator== (AutoDiffRec<D,SCAL> x, SCAL val2)
899  {
900  return x.Value() == val2;
901  }
902 
904  template <int D, typename SCAL>
905  INLINE bool operator!= (AutoDiffRec<D,SCAL> x, SCAL val2) throw()
906  {
907  return x.Value() != val2;
908  }
909 
911  template <int D, typename SCAL>
912  INLINE bool operator< (AutoDiffRec<D,SCAL> x, SCAL val2) throw()
913  {
914  return x.Value() < val2;
915  }
916 
918  template <int D, typename SCAL>
919  INLINE bool operator> (AutoDiffRec<D,SCAL> x, SCAL val2) throw()
920  {
921  return x.Value() > val2;
922  }
923 
924  using std::fabs;
925  template<int D, typename SCAL>
926  INLINE AutoDiffRec<D,SCAL> fabs (const AutoDiffRec<D,SCAL> & x)
927  {
928  auto sign = IfPos(x.Value(), SCAL(1.0), IfPos(-x.Value(), SCAL(-1.0), SCAL(0.0)));
929  return AutoDiffRec<D,SCAL> (fabs(x.Rec()), sign*x.Last());
930  // return fabs (AutoDiffVec<D,SCAL>(x));
931  /*
932  double abs = fabs (x.Value());
933  AutoDiffVec<D,SCAL> res( abs );
934  if (abs != 0.0)
935  for (int i = 0; i < D; i++)
936  res.DValue(i) = x.Value()*x.DValue(i) / abs;
937  else
938  for (int i = 0; i < D; i++)
939  res.DValue(i) = 0.0;
940  return res;
941  */
942  }
943 
944 
945  template<int D, typename SCAL>
946  INLINE auto sqrt (const AutoDiffRec<D,SCAL> & x)
947  {
948  return AutoDiffRec<D,SCAL> (sqrt(x.Rec()), (0.5/sqrt(x.Value()))*x.Last());
949  }
950 
951 
952 
953  template <int D, typename SCAL>
954  auto log (AutoDiffRec<D,SCAL> x)
955  {
956  return AutoDiffRec<D,SCAL> (log(x.Rec()), (1.0/x.Value())*x.Last());
957  }
958 
959  template <int D, typename SCAL>
960  auto exp (AutoDiffRec<D,SCAL> x)
961  {
962  return AutoDiffRec<D,SCAL> (exp(x.Rec()), exp(x.Value())*x.Last());
963  }
964 
965  template <int D, typename SCAL>
966  INLINE AutoDiffRec<D,SCAL> pow (AutoDiffRec<D,SCAL> x, AutoDiffRec<D,SCAL> y )
967  {
968  return exp(log(x)*y);
969  }
970 
971 
972  template <int D, typename SCAL>
973  auto sin (AutoDiffRec<D,SCAL> x)
974  {
975  return AutoDiffRec<D,SCAL> (sin(x.Rec()), cos(x.Value())*x.Last());
976  }
977 
978  template <int D, typename SCAL>
979  auto cos (AutoDiffRec<D,SCAL> x)
980  {
981  return AutoDiffRec<D,SCAL> (cos(x.Rec()), -sin(x.Value())*x.Last());
982  }
983 
984  template <int D, typename SCAL>
985  auto tan (AutoDiffRec<D,SCAL> x)
986  {
987  return sin(x) / cos(x);
988  }
989 
990  template <int D, typename SCAL>
991  auto sinh (AutoDiffRec<D,SCAL> x)
992  {
993  return AutoDiffRec<D,SCAL> (sinh(x.Rec()), cosh(x.Value())*x.Last());
994  }
995 
996  template <int D, typename SCAL>
997  auto cosh (AutoDiffRec<D,SCAL> x)
998  {
999  return AutoDiffRec<D,SCAL> (cosh(x.Rec()), sinh(x.Value())*x.Last());
1000  }
1001 
1002  template <int D, typename SCAL>
1003  auto floor (AutoDiffRec<D,SCAL> x)
1004  {
1005  return AutoDiffRec<D,SCAL> (floor(x.Rec()), 0.0);
1006  }
1007 
1008  template <int D, typename SCAL>
1009  auto ceil (AutoDiffRec<D,SCAL> x)
1010  {
1011  return AutoDiffRec<D,SCAL> (ceil(x.Rec()), 0.0);
1012  }
1013 
1014 
1015 
1016  template <int D, typename SCAL>
1017  auto atan (AutoDiffRec<D,SCAL> x)
1018  {
1019  return AutoDiffRec<D,SCAL> (atan(x.Rec()), (1./(1.+x.Value()*x.Value()))*x.Last());
1020  }
1021 
1022  template <int D, typename SCAL>
1023  auto atan2 (AutoDiffRec<D,SCAL> x, AutoDiffRec<D,SCAL> y)
1024  {
1025  return AutoDiffRec<D,SCAL> (atan2(x.Rec(), y.Rec()),
1026  (1./(x.Value()*x.Value()+y.Value()*y.Value()))*(x.Value()*y.Last()-y.Value()*x.Last()));
1027  }
1028 
1029  template <int D, typename SCAL>
1030  auto acos (AutoDiffRec<D,SCAL> x)
1031  {
1032  return AutoDiffRec<D,SCAL> (acos(x.Rec()), (-1./sqrt(1.-x.Value()*x.Value()))*x.Last());
1033  }
1034 
1035  template <int D, typename SCAL>
1036  auto asin (AutoDiffRec<D,SCAL> x)
1037  {
1038  return AutoDiffRec<D,SCAL> (asin(x.Rec()), (1./sqrt(1.-x.Value()*x.Value()))*x.Last());
1039  }
1040 
1041 
1042  template <int D, typename SCAL, typename TB, typename TC>
1043  auto IfPos (AutoDiffRec<D,SCAL> a, TB b, TC c) // -> decltype(IfPos (a.Value(), b, c))
1044  {
1045  return IfPos (a.Value(), b, c);
1046  }
1047 
1048  template <int D, typename SCAL>
1049  INLINE AutoDiffRec<D,SCAL> IfPos (SCAL /* SIMD<double> */ a, AutoDiffRec<D,SCAL> b, AutoDiffRec<D,SCAL> c)
1050  {
1051  /*
1052  AutoDiffRec<D,SCAL> res;
1053  res.Value() = IfPos (a, b.Value(), c.Value());
1054  for (int j = 0; j < D; j++)
1055  res.DValue(j) = IfPos (a, b.DValue(j), c.DValue(j));
1056  return res;
1057  */
1058  return AutoDiffRec<D,SCAL> (IfPos(a, b.Rec(), c.Rec()), IfPos(a, b.Last(), c.Last()));
1059  }
1060 
1061  template <int D, typename SCAL, typename TC>
1062  INLINE AutoDiffRec<D,SCAL> IfPos (SCAL /* SIMD<double> */ a, AutoDiffRec<D,SCAL> b, TC c)
1063  {
1064  return IfPos (a, b, AutoDiffRec<D,SCAL> (c));
1065  }
1066 
1067 
1068 
1069 template <int D, typename SCAL = double>
1070 using AutoDiff = AutoDiffRec<D,SCAL>;
1071 
1072 }
1073 
1074 #endif
1075 
1076 
ngstd::AutoDiffVec::operator=
INLINE AutoDiffVec & operator=(SCAL aval)
assign constant value
Definition: autodiff.hpp:73
ngstd::AutoDiffVec::AutoDiffVec
AutoDiffVec()=default
elements are undefined
ngstd::AutoDiffVec::AutoDiffVec
INLINE AutoDiffVec(SCAL aval)
initial object with constant value
Definition: autodiff.hpp:50
ngstd::operator<<
ostream & operator<<(ostream &ost, const AutoDiffDiff< D, SCAL > &x)
Prints AudoDiffDiff.
Definition: autodiffdiff.hpp:255
ngstd::AutoDiffVec::DValue
INLINE SCAL DValue(int i) const
returns partial derivative
Definition: autodiff.hpp:87
ngstd::AutoDiffVec::AutoDiffVec
INLINE AutoDiffVec(SCAL aval, int diffindex)
init object with (val, e_diffindex)
Definition: autodiff.hpp:58
ngstd::sqr
INLINE AutoDiffVec< D, SCAL > sqr(const AutoDiffVec< D, SCAL > &x)
AutoDiffVec times AutoDiffVec.
Definition: autodiff.hpp:250
ngstd::AutoDiffVec::Value
INLINE SCAL & Value()
access value
Definition: autodiff.hpp:103
ngstd
namespace for standard data types and algorithms.
Definition: ngstd.hpp:81
ngstd::AutoDiffRec
Definition: autodiff.hpp:15
ngstd::AutoDiffVec
Datatype for automatic differentiation.
Definition: autodiff.hpp:24
ngstd::AutoDiffVec::Value
INLINE SCAL Value() const
returns value
Definition: autodiff.hpp:84
ngstd::AutoDiffVec::DValue
INLINE SCAL & DValue(int i)
accesses partial derivative
Definition: autodiff.hpp:106
ngstd::Inv1
auto Inv1(SCAL x)
Inverse of AutoDiffRec.
Definition: autodiff.hpp:857
Generated by   doxygen 1.8.15