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