netgen/libsrc/core/autodiff.hpp

#ifndef FILE_AUTODIFF
#define FILE_AUTODIFF

/**************************************************************************/
/* File:   autodiff.hpp                                                   */
/* Author: Joachim Schoeberl                                              */
/* Date:   24. Oct. 02                                                    */
/**************************************************************************/

namespace ngcore
{
  using ngcore::IfPos;

// Automatic differentiation datatype

  template <int D, typename SCAL = double> class AutoDiffRec;


/**
   Datatype for automatic differentiation.
   Contains function value and D derivatives. Algebraic
   operations are overloaded by using product-rule etc. etc. 
**/
template <int D, typename SCAL = double>
class AutoDiffVec
{
  SCAL val;
  SCAL dval[D?D:1];
public:

  typedef AutoDiffVec<D,SCAL> TELEM;
  typedef SCAL TSCAL;


  /// elements are undefined
  // NETGEN_INLINE AutoDiffVec  () throw() { };
  AutoDiffVec() = default;
  // { val = 0; for (int i = 0; i < D; i++) dval[i] = 0; }  // !

  /// copy constructor
  AutoDiffVec  (const AutoDiffVec & ad2) = default;
  /*
  NETGEN_INLINE AutoDiffVec  (const AutoDiffVec & ad2) throw()
  {
    val = ad2.val;
    for (int i = 0; i < D; i++)
      dval[i] = ad2.dval[i];
  }
  */
  /// initial object with constant value
  NETGEN_INLINE AutoDiffVec  (SCAL aval) throw()
  {
    val = aval;
    for (int i = 0; i < D; i++)
      dval[i] = 0;
  }

  /// init object with (val, e_diffindex)
  NETGEN_INLINE AutoDiffVec  (SCAL aval, int diffindex)  throw()
  {
    val = aval;
    for (int i = 0; i < D; i++)
      dval[i] = 0;
    dval[diffindex] = 1;
  }

  NETGEN_INLINE AutoDiffVec (SCAL aval, const SCAL * grad)
  {
    val = aval;
    LoadGradient (grad);
  }

  /// assign constant value
  NETGEN_INLINE AutoDiffVec & operator= (SCAL aval) throw()
  {
    val = aval;
    for (int i = 0; i < D; i++)
      dval[i] = 0;
    return *this;
  }

  AutoDiffVec & operator= (const AutoDiffVec & ad2) = default;
  
  /// returns value
  NETGEN_INLINE SCAL Value() const throw() { return val; }
  
  /// returns partial derivative
  NETGEN_INLINE SCAL DValue (int i) const throw() { return dval[i]; }

  ///
  NETGEN_INLINE void StoreGradient (SCAL * p) const 
  {
    for (int i = 0; i < D; i++)
      p[i] = dval[i];
  }

  NETGEN_INLINE void LoadGradient (const SCAL * p) 
  {
    for (int i = 0; i < D; i++)
      dval[i] = p[i];
  }

  /// access value
  NETGEN_INLINE SCAL & Value() throw() { return val; }

  /// accesses partial derivative 
  NETGEN_INLINE SCAL & DValue (int i) throw() { return dval[i]; }
};


//@{  AutoDiffVec helper functions.

/// prints AutoDiffVec
template<int D, typename SCAL>
inline ostream & operator<< (ostream & ost, const AutoDiffVec<D,SCAL> & x)
{
  ost << x.Value() << ", D = ";
  for (int i = 0; i < D; i++)
    ost << x.DValue(i) << " ";
  return ost;
}

/// AutoDiffVec plus AutoDiffVec
template<int D, typename SCAL>
NETGEN_INLINE AutoDiffVec<D,SCAL> operator+ (const AutoDiffVec<D,SCAL> & x, const AutoDiffVec<D,SCAL> & y) throw()
{
  AutoDiffVec<D,SCAL> res;
  res.Value () = x.Value()+y.Value();
  // AutoDiffVec<D,SCAL> res(x.Value()+y.Value());
  for (int i = 0; i < D; i++)
    res.DValue(i) = x.DValue(i) + y.DValue(i);
  return res;
}


/// AutoDiffVec minus AutoDiffVec
template<int D, typename SCAL>
NETGEN_INLINE AutoDiffVec<D,SCAL> operator- (const AutoDiffVec<D,SCAL> & x, const AutoDiffVec<D,SCAL> & y) throw()
{
  AutoDiffVec<D,SCAL> res;
  res.Value() = x.Value()-y.Value();
  // AutoDiffVec<D,SCAL> res (x.Value()-y.Value());
  for (int i = 0; i < D; i++)
    res.DValue(i) = x.DValue(i) - y.DValue(i);
  return res;
}

/// double plus AutoDiffVec
  template<int D, typename SCAL, typename SCAL2,
           typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
NETGEN_INLINE AutoDiffVec<D,SCAL> operator+ (SCAL2 x, const AutoDiffVec<D,SCAL> & y) throw()
{
  AutoDiffVec<D,SCAL> res;
  res.Value() = x+y.Value();
  for (int i = 0; i < D; i++)
    res.DValue(i) = y.DValue(i);
  return res;
}

/// AutoDiffVec plus double
  template<int D, typename SCAL, typename SCAL2,
           typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
NETGEN_INLINE AutoDiffVec<D,SCAL> operator+ (const AutoDiffVec<D,SCAL> & y, SCAL2 x) throw()
{
  AutoDiffVec<D,SCAL> res;
  res.Value() = x+y.Value();
  for (int i = 0; i < D; i++)
    res.DValue(i) = y.DValue(i);
  return res;
}


/// minus AutoDiffVec
template<int D, typename SCAL>
NETGEN_INLINE AutoDiffVec<D,SCAL> operator- (const AutoDiffVec<D,SCAL> & x) throw()
{
  AutoDiffVec<D,SCAL> res;
  res.Value() = -x.Value();
  for (int i = 0; i < D; i++)
    res.DValue(i) = -x.DValue(i);
  return res;
}

/// AutoDiffVec minus double
  template<int D, typename SCAL, typename SCAL2,
           typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
NETGEN_INLINE AutoDiffVec<D,SCAL> operator- (const AutoDiffVec<D,SCAL> & x, SCAL2 y) throw()
{
  AutoDiffVec<D,SCAL> res;
  res.Value() = x.Value()-y;
  for (int i = 0; i < D; i++)
    res.DValue(i) = x.DValue(i);
  return res;
}

///
  template<int D, typename SCAL, typename SCAL2,
           typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
NETGEN_INLINE AutoDiffVec<D,SCAL> operator- (SCAL2 x, const AutoDiffVec<D,SCAL> & y) throw()
{
  AutoDiffVec<D,SCAL> res;
  res.Value() = x-y.Value();
  for (int i = 0; i < D; i++)
    res.DValue(i) = -y.DValue(i);
  return res;
}


/// double times AutoDiffVec
  template<int D, typename SCAL, typename SCAL2,
           typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
NETGEN_INLINE AutoDiffVec<D,SCAL> operator* (SCAL2 x, const AutoDiffVec<D,SCAL> & y) throw()
{
  AutoDiffVec<D,SCAL> res;
  res.Value() = x*y.Value();
  for (int i = 0; i < D; i++)
    res.DValue(i) = x*y.DValue(i);
  return res;
}

/// AutoDiffVec times double
  template<int D, typename SCAL, typename SCAL2,
           typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>

  NETGEN_INLINE AutoDiffVec<D,SCAL> operator* (const AutoDiffVec<D,SCAL> & y, SCAL2 x) throw()
{
  AutoDiffVec<D,SCAL> res;
  res.Value() = x*y.Value();
  for (int i = 0; i < D; i++)
    res.DValue(i) = x*y.DValue(i);
  return res;
}

/// AutoDiffVec times AutoDiffVec
template<int D, typename SCAL>
NETGEN_INLINE AutoDiffVec<D,SCAL> operator* (const AutoDiffVec<D,SCAL> & x, const AutoDiffVec<D,SCAL> & y) throw()
{
  AutoDiffVec<D,SCAL> res;
  SCAL hx = x.Value();
  SCAL hy = y.Value();

  res.Value() = hx*hy;
  for (int i = 0; i < D; i++)
    res.DValue(i) = hx*y.DValue(i) + hy*x.DValue(i);

  return res;
}

/// AutoDiffVec times AutoDiffVec
using ngcore::sqr;
template<int D, typename SCAL>
NETGEN_INLINE AutoDiffVec<D,SCAL> sqr (const AutoDiffVec<D,SCAL> & x) throw()
{
  AutoDiffVec<D,SCAL> res;
  SCAL hx = x.Value();
  res.Value() = hx*hx;
  hx *= 2;
  for (int i = 0; i < D; i++)
    res.DValue(i) = hx*x.DValue(i);
  return res;
}

/// Inverse of AutoDiffVec
template<int D, typename SCAL>
NETGEN_INLINE AutoDiffVec<D,SCAL> Inv (const AutoDiffVec<D,SCAL> & x)
{
  AutoDiffVec<D,SCAL> res(1.0 / x.Value());
  for (int i = 0; i < D; i++)
    res.DValue(i) = -x.DValue(i) / (x.Value() * x.Value());
  return res;
}


/// AutoDiffVec div AutoDiffVec
template<int D, typename SCAL>
NETGEN_INLINE AutoDiffVec<D,SCAL> operator/ (const AutoDiffVec<D,SCAL> & x, const AutoDiffVec<D,SCAL> & y)
{
  return x * Inv (y);
}

/// AutoDiffVec div double
template<int D, typename SCAL, typename SCAL2,
         typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
NETGEN_INLINE AutoDiffVec<D,SCAL> operator/ (const AutoDiffVec<D,SCAL> & x, SCAL2 y)
{
  return (1.0/y) * x;
}

  /// double div AutoDiffVec
template<int D, typename SCAL, typename SCAL2,
         typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
  NETGEN_INLINE AutoDiffVec<D,SCAL> operator/ (SCAL2 x, const AutoDiffVec<D,SCAL> & y)
  {
    return x * Inv(y);
  }
  


  
  template <int D, typename SCAL, typename SCAL2>
  NETGEN_INLINE AutoDiffVec<D,SCAL> & operator+= (AutoDiffVec<D,SCAL> & x, SCAL2 y) throw()
  {
    x.Value() += y;
    return x;
  }


  /// 
  template <int D, typename SCAL>
  NETGEN_INLINE AutoDiffVec<D,SCAL> & operator+= (AutoDiffVec<D,SCAL> & x, AutoDiffVec<D,SCAL> y)
  {
    x.Value() += y.Value();
    for (int i = 0; i < D; i++)
      x.DValue(i) += y.DValue(i);
    return x;
  }

  ///
  template <int D, typename SCAL>
  NETGEN_INLINE AutoDiffVec<D,SCAL> & operator-= (AutoDiffVec<D,SCAL> & x, AutoDiffVec<D,SCAL> y)
  {
    x.Value() -= y.Value();
    for (int i = 0; i < D; i++)
      x.DValue(i) -= y.DValue(i);
    return x;

  }

  template <int D, typename SCAL, typename SCAL2>
  NETGEN_INLINE AutoDiffVec<D,SCAL> & operator-= (AutoDiffVec<D,SCAL> & x, SCAL2 y)
  {
    x.Value() -= y;
    return x;
  }

  ///
  template <int D, typename SCAL>
  NETGEN_INLINE AutoDiffVec<D,SCAL> & operator*= (AutoDiffVec<D,SCAL> & x, AutoDiffVec<D,SCAL> y) 
  {
    for (int i = 0; i < D; i++)
      x.DValue(i) = x.DValue(i)*y.Value() + x.Value() * y.DValue(i);
    x.Value() *= y.Value();
    return x;
  }

  ///
  template <int D, typename SCAL, typename SCAL2>
  NETGEN_INLINE AutoDiffVec<D,SCAL> & operator*= (AutoDiffVec<D,SCAL> & x, SCAL2 y) 
  {
    x.Value() *= y;
    for (int i = 0; i < D; i++)
      x.DValue(i) *= y;
    return x;
  }

  ///
  template <int D, typename SCAL>
  NETGEN_INLINE AutoDiffVec<D,SCAL> & operator/= (AutoDiffVec<D,SCAL> & x, SCAL y) 
  {
    SCAL iy = 1.0 / y;
    x.Value() *= iy;
    for (int i = 0; i < D; i++)
      x.DValue(i) *= iy;
    return x;
  }




  /// 
  template <int D, typename SCAL>
  NETGEN_INLINE bool operator== (AutoDiffVec<D,SCAL> x, SCAL val2) 
  {
    return x.Value() == val2;
  }

  ///
  template <int D, typename SCAL>
  NETGEN_INLINE bool operator!= (AutoDiffVec<D,SCAL> x, SCAL val2) throw()
  {
    return x.Value() != val2;
  }

  ///
  template <int D, typename SCAL>
  NETGEN_INLINE bool operator< (AutoDiffVec<D,SCAL> x, SCAL val2) throw()
  {
    return x.Value() < val2;
  }
  
  ///
  template <int D, typename SCAL>
  NETGEN_INLINE bool operator> (AutoDiffVec<D,SCAL> x, SCAL val2) throw()
  {
    return x.Value() > val2;
  }




template<int D, typename SCAL>
NETGEN_INLINE AutoDiffVec<D,SCAL> fabs (const AutoDiffVec<D,SCAL> & x)
{
  double abs = fabs (x.Value());
  AutoDiffVec<D,SCAL> res( abs );
  if (abs != 0.0)
    for (int i = 0; i < D; i++)
      res.DValue(i) = x.Value()*x.DValue(i) / abs;
  else
    for (int i = 0; i < D; i++)
      res.DValue(i) = 0.0;
  return res;
}

using std::sqrt;
template<int D, typename SCAL>
NETGEN_INLINE AutoDiffVec<D,SCAL> sqrt (const AutoDiffVec<D,SCAL> & x)
{
  AutoDiffVec<D,SCAL> res;
  res.Value() = sqrt(x.Value());
  for (int j = 0; j < D; j++)
    res.DValue(j) = 0.5 / res.Value() * x.DValue(j);
  return res;
}

using std::log;
template <int D, typename SCAL>
AutoDiffVec<D,SCAL> log (AutoDiffVec<D,SCAL> x)
{
  AutoDiffVec<D,SCAL> res;
  res.Value() = log(x.Value());
  for (int k = 0; k < D; k++)
    res.DValue(k) = x.DValue(k) / x.Value();
  return res;
}

using std::exp;
template <int D, typename SCAL>
NETGEN_INLINE AutoDiffVec<D,SCAL> exp (AutoDiffVec<D,SCAL> x)
{
  AutoDiffVec<D,SCAL> res;
  res.Value() = exp(x.Value());
  for (int k = 0; k < D; k++)
    res.DValue(k) = x.DValue(k) * res.Value();
  return res;
}

using std::pow;
template <int D, typename SCAL>
NETGEN_INLINE AutoDiffVec<D,SCAL> pow (AutoDiffVec<D,SCAL> x, AutoDiffVec<D,SCAL> y )
{
  return exp(log(x)*y);
}

using std::sin;
/*
template <int D, typename SCAL>
NETGEN_INLINE AutoDiffVec<D,SCAL> sin (AutoDiffVec<D,SCAL> x)
{
  AutoDiffVec<D,SCAL> res;
  res.Value() = sin(x.Value());
  SCAL c = cos(x.Value());
  for (int k = 0; k < D; k++)
    res.DValue(k) = x.DValue(k) * c;
  return res;
}
*/

template <int D, typename SCAL>
NETGEN_INLINE AutoDiffVec<D,SCAL> sin (AutoDiffVec<D,SCAL> x)
{
  return sin(AutoDiffRec<D,SCAL>(x));
}
  
using std::cos;
/*
template <int D, typename SCAL>
NETGEN_INLINE AutoDiffVec<D,SCAL> cos (AutoDiffVec<D,SCAL> x)
{
  AutoDiffVec<D,SCAL> res;
  res.Value() = cos(x.Value());
  SCAL ms = -sin(x.Value());
  for (int k = 0; k < D; k++)
    res.DValue(k) = x.DValue(k) * ms;
  return res;
}
*/
template <int D, typename SCAL>
NETGEN_INLINE AutoDiffVec<D,SCAL> cos (AutoDiffVec<D,SCAL> x)
{
  return cos(AutoDiffRec<D,SCAL>(x));
}

using std::tan;
template <int D, typename SCAL>
NETGEN_INLINE AutoDiffVec<D,SCAL> tan (AutoDiffVec<D,SCAL> x)
{ return sin(x) / cos(x); }

using std::sinh;
template <int D, typename SCAL>
NETGEN_INLINE AutoDiffVec<D,SCAL> sinh (AutoDiffVec<D,SCAL> x)
{
  AutoDiffVec<D,SCAL> res;
  res.Value() = sinh(x.Value());
  SCAL ch = cosh(x.Value());
  for (int k = 0; k < D; k++)
    res.DValue(k) = x.DValue(k) * ch;
  return res;
}

using std::cosh;
template <int D, typename SCAL>
NETGEN_INLINE AutoDiffVec<D,SCAL> cosh (AutoDiffVec<D,SCAL> x)
{
  AutoDiffVec<D,SCAL> res;
  res.Value() = cosh(x.Value());
  SCAL sh = sinh(x.Value());
  for (int k = 0; k < D; k++)
    res.DValue(k) = x.DValue(k) * sh;
  return res;
}

using std::erf;
template <int D, typename SCAL>
NETGEN_INLINE AutoDiffVec<D,SCAL> erf (AutoDiffVec<D,SCAL> x)
{
  return erf(AutoDiffRec<D,SCAL>(x));
}

using std::floor;
template<int D, typename SCAL>
NETGEN_INLINE AutoDiffVec<D,SCAL> floor (const AutoDiffVec<D,SCAL> & x)
{
  AutoDiffVec<D,SCAL> res;
  res.Value() = floor(x.Value());
  for (int j = 0; j < D; j++)
    res.DValue(j) = 0.0;
  return res;
}

using std::ceil;
template<int D, typename SCAL>
NETGEN_INLINE AutoDiffVec<D,SCAL> ceil (const AutoDiffVec<D,SCAL> & x)
{
  AutoDiffVec<D,SCAL> res;
  res.Value() = ceil(x.Value());
  for (int j = 0; j < D; j++)
    res.DValue(j) = 0.0;
  return res;
}


using std::atan;
/*
template <int D, typename SCAL>
NETGEN_INLINE AutoDiffVec<D,SCAL> atan (AutoDiffVec<D,SCAL> x)
{
  AutoDiffVec<D,SCAL> res;
  SCAL a = atan(x.Value());
  res.Value() = a;
  for (int k = 0; k < D; k++)
    res.DValue(k) = x.DValue(k)/(1+x.Value()*x.Value()) ;
  return res;
}
*/
template <int D, typename SCAL>
AutoDiffVec<D,SCAL> atan (AutoDiffVec<D,SCAL> x)
{
  return atan (AutoDiffRec<D,SCAL> (x));
}

using std::atan2;
template <int D, typename SCAL>
NETGEN_INLINE AutoDiffVec<D,SCAL> atan2 (AutoDiffVec<D,SCAL> x, AutoDiffVec<D,SCAL> y)
{
  AutoDiffVec<D,SCAL> res;
  SCAL a = atan2(x.Value(), y.Value());
  res.Value() = a;
  for (int k = 0; k < D; k++)
    res.DValue(k) = (x.Value()*y.DValue(k)-y.Value()*x.DValue(k))/(y.Value()*y.Value()+x.Value()*x.Value());
  return res;
}


using std::acos;
template <int D, typename SCAL>
NETGEN_INLINE AutoDiffVec<D,SCAL> acos (AutoDiffVec<D,SCAL> x)
{
  AutoDiffVec<D,SCAL> res;
  SCAL a = acos(x.Value());
  res.Value() = a;
  SCAL da = -1 / sqrt(1-x.Value()*x.Value());
  for (int k = 0; k < D; k++)
    res.DValue(k) = x.DValue(k)*da;
  return res;
}


using std::asin;
template <int D, typename SCAL>
NETGEN_INLINE AutoDiffVec<D,SCAL> asin (AutoDiffVec<D,SCAL> x)
{
  AutoDiffVec<D,SCAL> res;
  SCAL a = asin(x.Value());
  res.Value() = a;
  SCAL da = 1 / sqrt(1-x.Value()*x.Value());
  for (int k = 0; k < D; k++)
    res.DValue(k) = x.DValue(k)*da;
  return res;
}




  template <int D, typename SCAL, typename TB, typename TC>
  auto IfPos (AutoDiffVec<D,SCAL> a, TB b, TC c) // -> decltype(IfPos (a.Value(), b, c))
  {
    return IfPos (a.Value(), b, c);
  }

  template <int D, typename SCAL>
  NETGEN_INLINE AutoDiffVec<D,SCAL> IfPos (SCAL /* SIMD<double> */ a, AutoDiffVec<D,SCAL> b, AutoDiffVec<D,SCAL> c)
  {
    AutoDiffVec<D,SCAL> res;
    res.Value() = IfPos (a, b.Value(), c.Value());
    for (int j = 0; j < D; j++)
      res.DValue(j) = IfPos (a, b.DValue(j), c.DValue(j));
    return res;
  }

  template <int D, typename SCAL, typename TC>
  NETGEN_INLINE AutoDiffVec<D,SCAL> IfPos (SCAL /* SIMD<double> */ a, AutoDiffVec<D,SCAL> b, TC c)
  {
    return IfPos (a, b, AutoDiffVec<D,SCAL> (c));
  }

//@}


  
  template <int D, typename SCAL>
  class AutoDiffRec
  {
    AutoDiffRec<D-1, SCAL> rec;
    SCAL last;

  public:
    NETGEN_INLINE AutoDiffRec () = default;
    NETGEN_INLINE AutoDiffRec (const AutoDiffRec &) = default;
    NETGEN_INLINE AutoDiffRec (AutoDiffRec<D-1,SCAL> _rec, SCAL _last) : rec(_rec), last(_last) { ; }
    NETGEN_INLINE AutoDiffRec & operator= (const AutoDiffRec &) = default;

    NETGEN_INLINE AutoDiffRec (SCAL aval) : rec(aval), last(0.0) { ; }
    NETGEN_INLINE AutoDiffRec (SCAL aval, int diffindex) : rec(aval, diffindex), last((diffindex==D-1) ? 1.0 : 0.0) { ; }
    NETGEN_INLINE AutoDiffRec (SCAL aval, const SCAL * grad)
      : rec(aval, grad), last(grad[D-1]) { } 

    NETGEN_INLINE AutoDiffRec (const AutoDiffVec<D,SCAL> & ad)
    {
      Value() = ad.Value();
      for (int i = 0; i < D; i++)
        DValue(i) = ad.DValue(i);
    }
    
    NETGEN_INLINE AutoDiffRec & operator= (SCAL aval) { rec = aval; last = 0.0; return *this; }
    NETGEN_INLINE SCAL Value() const { return rec.Value(); }
    NETGEN_INLINE SCAL DValue(int i) const { return (i == D-1) ? last : rec.DValue(i); }
    NETGEN_INLINE SCAL & Value() { return rec.Value(); }
    NETGEN_INLINE SCAL & DValue(int i) { return (i == D-1) ? last : rec.DValue(i); }
    NETGEN_INLINE auto Rec() const { return rec; }
    NETGEN_INLINE auto Last() const { return last; }
    NETGEN_INLINE auto & Rec() { return rec; }
    NETGEN_INLINE auto & Last() { return last; }
    NETGEN_INLINE operator AutoDiffVec<D,SCAL> () const
    {
      AutoDiffVec<D,SCAL> res(Value());
      for (int i = 0; i < D; i++)
        res.DValue(i) = DValue(i);
      return res;
    }
  };

  template<int D, typename SCAL>
  ostream & operator<< (ostream & ost, AutoDiffRec<D,SCAL> ad)
  {
    return ost << AutoDiffVec<D,SCAL> (ad);
  }

  template <typename SCAL>
  class AutoDiffRec<0,SCAL>
  {
    SCAL val;
  public:
    NETGEN_INLINE AutoDiffRec () = default;
    NETGEN_INLINE AutoDiffRec (const AutoDiffRec &) = default;
    NETGEN_INLINE AutoDiffRec (SCAL _val) : val(_val) { ; }
    NETGEN_INLINE AutoDiffRec (SCAL _val, SCAL /* _dummylast */) : val(_val) { ; }
    NETGEN_INLINE AutoDiffRec (SCAL aval, const SCAL * /* grad */)
      : val(aval) { } 

    NETGEN_INLINE AutoDiffRec & operator= (const AutoDiffRec &) = default;
    NETGEN_INLINE AutoDiffRec & operator= (SCAL aval) { val = aval; return *this; }

    NETGEN_INLINE SCAL Value() const { return val; }
    NETGEN_INLINE SCAL DValue(int /* i */) const { return SCAL(0); }
    NETGEN_INLINE SCAL & Value() { return val; }
    // SCAL & DValue(int i) { return val; }
    NETGEN_INLINE auto Rec() const { return val; }
    NETGEN_INLINE auto Last() const { return SCAL(0); }
    NETGEN_INLINE auto & Rec() { return val; }
    NETGEN_INLINE auto & Last() { return val; }
    NETGEN_INLINE operator AutoDiffVec<0,SCAL> () const { return AutoDiffVec<0,SCAL>(); }
  };


  template <typename SCAL>
  class AutoDiffRec<1,SCAL>
  {
    SCAL val;
    SCAL last;
  public:
    NETGEN_INLINE AutoDiffRec () = default;
    NETGEN_INLINE AutoDiffRec (const AutoDiffRec &) = default;
    NETGEN_INLINE AutoDiffRec (SCAL _val) : val(_val), last(0.0) { ; }
    NETGEN_INLINE AutoDiffRec (SCAL _val, SCAL _last) : val(_val), last(_last) { ; }
    NETGEN_INLINE AutoDiffRec (SCAL aval, int diffindex) : val(aval), last((diffindex==0) ? 1.0 : 0.0) { ; }
    NETGEN_INLINE AutoDiffRec (SCAL aval, const SCAL * grad)
      : val(aval), last(grad[0]) { } 
    
    NETGEN_INLINE AutoDiffRec (const AutoDiffVec<1,SCAL> & ad)
    {
      Value() = ad.Value();
      DValue(0) = ad.DValue(0);
    }

    NETGEN_INLINE AutoDiffRec & operator= (const AutoDiffRec &) = default;
    NETGEN_INLINE AutoDiffRec & operator= (SCAL aval) { val = aval; last = 0.0; return *this; }

    NETGEN_INLINE SCAL Value() const { return val; }
    NETGEN_INLINE SCAL DValue(int /* i */) const { return last; }
    NETGEN_INLINE SCAL & Value() { return val; }
    NETGEN_INLINE SCAL & DValue(int /* i */) { return last; }
    NETGEN_INLINE auto Rec() const { return val; }
    NETGEN_INLINE auto Last() const { return last; }
    NETGEN_INLINE auto & Rec() { return val; }
    NETGEN_INLINE auto & Last() { return last; }

    NETGEN_INLINE operator AutoDiffVec<1,SCAL> () const
    {
      AutoDiffVec<1,SCAL> res(Value());
      res.DValue(0) = DValue(0);
      return res;
    }
  };

  template <int D, typename SCAL, typename SCAL2,
            typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
  NETGEN_INLINE AutoDiffRec<D,SCAL> operator+ (SCAL2 a, AutoDiffRec<D,SCAL> b)
  {
    return AutoDiffRec<D,SCAL> (a+b.Rec(), b.Last());
  }
  
  template <int D, typename SCAL, typename SCAL2,
            typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
  NETGEN_INLINE AutoDiffRec<D,SCAL> operator+ (AutoDiffRec<D,SCAL> a, SCAL2 b)
  {
    return AutoDiffRec<D,SCAL> (a.Rec()+b, a.Last());
  }
  
  template <int D, typename SCAL>  
  NETGEN_INLINE AutoDiffRec<D,SCAL> operator+ (AutoDiffRec<D,SCAL> a, AutoDiffRec<D,SCAL> b)
  {
    return AutoDiffRec<D,SCAL> (a.Rec()+b.Rec(), a.Last()+b.Last());
  }

  template <int D, typename SCAL, typename SCAL2,
            typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
  NETGEN_INLINE AutoDiffRec<D,SCAL> operator- (SCAL2 b, AutoDiffRec<D,SCAL> a)
  {
    return AutoDiffRec<D,SCAL> (b-a.Rec(), -a.Last());
  }
  
  template <int D, typename SCAL, typename SCAL2, 
            typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
  NETGEN_INLINE AutoDiffRec<D,SCAL> operator- (AutoDiffRec<D,SCAL> a, SCAL2 b)
  {
    return AutoDiffRec<D,SCAL> (a.Rec()-b, a.Last());
  }
  
  template <int D, typename SCAL>  
  NETGEN_INLINE AutoDiffRec<D,SCAL> operator- (AutoDiffRec<D,SCAL> a, AutoDiffRec<D,SCAL> b)
  {
    return AutoDiffRec<D,SCAL> (a.Rec()-b.Rec(), a.Last()-b.Last());
  }

  /// minus AutoDiff
  template<int D, typename SCAL>
  NETGEN_INLINE AutoDiffRec<D,SCAL> operator- (AutoDiffRec<D,SCAL> a)
  {
    return AutoDiffRec<D,SCAL> (-a.Rec(), -a.Last());
  }

  template <int D, typename SCAL>  
  NETGEN_INLINE AutoDiffRec<D,SCAL> operator* (AutoDiffRec<D,SCAL> a, AutoDiffRec<D,SCAL> b)
  {
    return AutoDiffRec<D,SCAL> (a.Rec()*b.Rec(), a.Value()*b.Last()+b.Value()*a.Last());
  }

  template <int D, typename SCAL, typename SCAL1,
            typename std::enable_if<std::is_convertible<SCAL1,SCAL>::value, int>::type = 0>
  NETGEN_INLINE AutoDiffRec<D,SCAL> operator* (AutoDiffRec<D,SCAL> b, SCAL1 a)
  {
    return AutoDiffRec<D,SCAL> (a*b.Rec(), a*b.Last());
  }

  template <int D, typename SCAL, typename SCAL1,
            typename std::enable_if<std::is_convertible<SCAL1,SCAL>::value, int>::type = 0>
  NETGEN_INLINE AutoDiffRec<D,SCAL> operator* (SCAL1 a, AutoDiffRec<D,SCAL> b)
  {
    return AutoDiffRec<D,SCAL> (a*b.Rec(), a*b.Last());
  }

  template <int D, typename SCAL>
  NETGEN_INLINE AutoDiffRec<D,SCAL> & operator+= (AutoDiffRec<D,SCAL> & a, AutoDiffRec<D,SCAL> b)
  {
    a.Rec() += b.Rec();
    a.Last() += b.Last();
    return a;
  }

  template <int D, typename SCAL>
  NETGEN_INLINE AutoDiffRec<D,SCAL> & operator-= (AutoDiffRec<D,SCAL> & a, double b)
  {
    a.Rec() -= b;
    return a;
  }

  template <int D, typename SCAL>
  NETGEN_INLINE AutoDiffRec<D,SCAL> & operator-= (AutoDiffRec<D,SCAL> & a, AutoDiffRec<D,SCAL> b)
  {
    a.Rec() -= b.Rec();
    a.Last() -= b.Last();
    return a;
  }


  template <int D, typename SCAL>  
  NETGEN_INLINE AutoDiffRec<D,SCAL> & operator*= (AutoDiffRec<D,SCAL> & a, AutoDiffRec<D,SCAL> b)
  {
    a = a*b;
    return a;
  }

  
  template <int D, typename SCAL, typename SCAL2>  
  NETGEN_INLINE AutoDiffRec<D,SCAL> & operator*= (AutoDiffRec<D,SCAL> & b, SCAL2 a)
  {
    b.Rec() *= a;
    b.Last() *= a;
    return b;
  }

  /// Inverse of AutoDiffRec

  template <typename SCAL>
  auto Inv1 (SCAL x)  { return 1.0/x; }

  template<int D, typename SCAL>
  NETGEN_INLINE AutoDiffRec<D,SCAL> Inv1 (AutoDiffRec<D,SCAL> x)
  {
    return AutoDiffRec<D,SCAL> (Inv1(x.Rec()), (-sqr(1.0/x.Value())) * x.Last());
  }

  /// AutoDiffRec div AutoDiffRec
  template<int D, typename SCAL>
  NETGEN_INLINE AutoDiffRec<D,SCAL> operator/ (const AutoDiffRec<D,SCAL> & x, const AutoDiffRec<D,SCAL> & y)
  {
    return x * Inv1 (y);
  }


  /// AutoDiffVec div double
  template<int D, typename SCAL, typename SCAL2,
           typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
  NETGEN_INLINE AutoDiffRec<D,SCAL> operator/ (const AutoDiffRec<D,SCAL> & x, SCAL2 y)
  {
    return (1.0/y) * x;
  }


  /// double div AutoDiffVec
  template<int D, typename SCAL, typename SCAL2,
           typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
  NETGEN_INLINE AutoDiffRec<D,SCAL> operator/ (SCAL2 x, const AutoDiffRec<D,SCAL> & y)
  {
    return x * Inv1(y);
  }







  /// 
  template <int D, typename SCAL>
  NETGEN_INLINE bool operator== (AutoDiffRec<D,SCAL> x, SCAL val2) 
  {
    return x.Value() == val2;
  }

  ///
  template <int D, typename SCAL>
  NETGEN_INLINE bool operator!= (AutoDiffRec<D,SCAL> x, SCAL val2) throw()
  {
    return x.Value() != val2;
  }

  ///
  template <int D, typename SCAL>
  NETGEN_INLINE bool operator< (AutoDiffRec<D,SCAL> x, SCAL val2) throw()
  {
    return x.Value() < val2;
  }
  
  ///
  template <int D, typename SCAL>
  NETGEN_INLINE bool operator> (AutoDiffRec<D,SCAL> x, SCAL val2) throw()
  {
    return x.Value() > val2;
  }

  using std::fabs;
  template<int D, typename SCAL>
  NETGEN_INLINE AutoDiffRec<D,SCAL> fabs (const AutoDiffRec<D,SCAL> & x)
  {
    auto sign = IfPos(x.Value(), SCAL(1.0), IfPos(-x.Value(), SCAL(-1.0), SCAL(0.0)));
    return AutoDiffRec<D,SCAL> (fabs(x.Rec()), sign*x.Last());
    // return fabs (AutoDiffVec<D,SCAL>(x));
    /*
      double abs = fabs (x.Value());
      AutoDiffVec<D,SCAL> res( abs );
      if (abs != 0.0)
      for (int i = 0; i < D; i++)
      res.DValue(i) = x.Value()*x.DValue(i) / abs;
      else
      for (int i = 0; i < D; i++)
      res.DValue(i) = 0.0;
      return res;
    */
  }


  template<int D, typename SCAL>
  NETGEN_INLINE auto sqrt (const AutoDiffRec<D,SCAL> & x)
  {
    return AutoDiffRec<D,SCAL> (sqrt(x.Rec()), (0.5/sqrt(x.Value()))*x.Last());
  }



  template <int D, typename SCAL>
  auto log (AutoDiffRec<D,SCAL> x)
  {
    return AutoDiffRec<D,SCAL> (log(x.Rec()), (1.0/x.Value())*x.Last());
  }

  template <int D, typename SCAL>
  auto exp (AutoDiffRec<D,SCAL> x)
  {
    return AutoDiffRec<D,SCAL> (exp(x.Rec()), exp(x.Value())*x.Last());
  }
  
  template <int D, typename SCAL>
  NETGEN_INLINE AutoDiffRec<D,SCAL> pow (AutoDiffRec<D,SCAL> x, AutoDiffRec<D,SCAL> y )
  {
    return exp(log(x)*y);
  }


  template <int D, typename SCAL>
  auto sin (AutoDiffRec<D,SCAL> x)
  {
    return AutoDiffRec<D,SCAL> (sin(x.Rec()), cos(x.Value())*x.Last());
  }

  template <int D, typename SCAL>
  auto cos (AutoDiffRec<D,SCAL> x)
  {
    return AutoDiffRec<D,SCAL> (cos(x.Rec()), -sin(x.Value())*x.Last());
  }

  template <int D, typename SCAL>
  auto tan (AutoDiffRec<D,SCAL> x)
  {
    return sin(x) / cos(x);
  }

  template <int D, typename SCAL>
  auto sinh (AutoDiffRec<D,SCAL> x)
  {
    return AutoDiffRec<D,SCAL> (sinh(x.Rec()), cosh(x.Value())*x.Last());
  }

  template <int D, typename SCAL>
  auto cosh (AutoDiffRec<D,SCAL> x)
  {
    return AutoDiffRec<D,SCAL> (cosh(x.Rec()), sinh(x.Value())*x.Last());
  }

  template <int D, typename SCAL>
  auto erf (AutoDiffRec<D,SCAL> x)
  {
    return AutoDiffRec<D,SCAL> (erf(x.Rec()), 2. / sqrt(M_PI) * exp(- x.Value() * x.Value())*x.Last());
  }

  template <int D, typename SCAL>
  auto floor (AutoDiffRec<D,SCAL> x)
  {
    return AutoDiffRec<D,SCAL> (floor(x.Rec()), 0.0);
  }

  template <int D, typename SCAL>
  auto ceil (AutoDiffRec<D,SCAL> x)
  {
    return AutoDiffRec<D,SCAL> (ceil(x.Rec()), 0.0);
  }



  template <int D, typename SCAL>
  auto atan (AutoDiffRec<D,SCAL> x)
  {
    return AutoDiffRec<D,SCAL> (atan(x.Rec()), (1./(1.+x.Value()*x.Value()))*x.Last());
  }

  template <int D, typename SCAL>
  auto atan2 (AutoDiffRec<D,SCAL> x, AutoDiffRec<D,SCAL> y)
  {
    return AutoDiffRec<D,SCAL> (atan2(x.Rec(), y.Rec()),
                                (1./(x.Value()*x.Value()+y.Value()*y.Value()))*(x.Value()*y.Last()-y.Value()*x.Last()));
  }

  template <int D, typename SCAL>
  auto acos (AutoDiffRec<D,SCAL> x)
  {
    return AutoDiffRec<D,SCAL> (acos(x.Rec()), (-1./sqrt(1.-x.Value()*x.Value()))*x.Last());
  }

  template <int D, typename SCAL>
  auto asin (AutoDiffRec<D,SCAL> x)
  {
    return AutoDiffRec<D,SCAL> (asin(x.Rec()), (1./sqrt(1.-x.Value()*x.Value()))*x.Last());
  }


  template <int D, typename SCAL, typename TB, typename TC>
  auto IfPos (AutoDiffRec<D,SCAL> a, TB b, TC c) // -> decltype(IfPos (a.Value(), b, c))
  {
    return IfPos (a.Value(), b, c);
  }

  template <int D, typename SCAL>
  NETGEN_INLINE AutoDiffRec<D,SCAL> IfPos (SCAL /* SIMD<double> */ a, AutoDiffRec<D,SCAL> b, AutoDiffRec<D,SCAL> c)
  {
    /*
    AutoDiffRec<D,SCAL> res;
    res.Value() = IfPos (a, b.Value(), c.Value());
    for (int j = 0; j < D; j++)
      res.DValue(j) = IfPos (a, b.DValue(j), c.DValue(j));
    return res;
    */
    return AutoDiffRec<D,SCAL> (IfPos(a, b.Rec(), c.Rec()), IfPos(a, b.Last(), c.Last()));
  }

  template <int D, typename SCAL, typename TC>
  NETGEN_INLINE AutoDiffRec<D,SCAL> IfPos (SCAL /* SIMD<double> */ a, AutoDiffRec<D,SCAL> b, TC c)
  {
    return IfPos (a, b, AutoDiffRec<D,SCAL> (c));
  }



template <int D, typename SCAL = double>
using AutoDiff = AutoDiffRec<D,SCAL>;

}



namespace ngbla
{
  template <typename T> struct is_scalar_type;
  template <int D, typename T>
  struct is_scalar_type<ngcore::AutoDiff<D,T>> { static constexpr bool value = true; };

  // not meaningful for AutoDiff<D,Complex>, since this is
  // not (complex) differentiable anyway
  template<int D, typename SCAL>
  inline auto L2Norm2 (const ngcore::AutoDiff<D,SCAL> & x) 
  {
    return x*x;
  }

  template<int D, typename SCAL>
  inline auto L2Norm (const ngcore::AutoDiff<D,SCAL> & x) throw()
  {
    return IfPos(x.Value(), x, -x);
  }
  

  
  template<int D, typename TAD>
  NETGEN_INLINE auto Conj (const ngcore::AutoDiff<D,TAD> & a) 
  { 
    ngcore::AutoDiff<D,TAD> b; 
    b.Value() = conj(a.Value()); 
    
    for(int i=0;i<D;i++) 
      b.DValue(i) = conj(a.DValue(i)); 
    
    return b;
  }

  
}


#endif