ad3.h
Go to the documentation of this file.
#ifndef _RHEO_AD_POINT_H
#define _RHEO_AD_POINT_H
/*Class:ad3
NAME: @code{ad3} - automatic differentiation with respect to variables in R^3
@clindex ad3
DESCRIPTION:
@noindent
The @code{ad3} class defines a forward automatic differentiation
in R^3 for computing automatically the gradient of a function from R^3 to R.
The implementation uses a simple forward automatic differentiation method.
EXAMPLE:
@example
template<class T>
T f (const point_basic<T>& x) @{ return sqr(x[0]) + x[0]*x[1] + 1/x[2]; @}
...
point_basic<ad3> x_ad = ad3::point (x);
ad3 y_ad = f(x_ad);
cout << "f ="<<y_ad <<endl
<< "gradq(f)="<<y_ad.grad() <<endl;
@end example
AUTHOR: Pierre.Saramito@imag.fr
DATE: 26 september 2017
End:
*/
#include "rheolef/point.h"
namespace rheolef {
template <class T>
class ad3_basic {
public:
// typedefs:
typedef typename point_basic<T>::size_type size_type;
typedef T value_type;
// allocators:
ad3_basic ();
ad3_basic (const ad3_basic&);
ad3_basic<T>& operator= (const ad3_basic<T>&);
template <class U>
ad3_basic(const U& x);
// helper:
static point_basic<ad3_basic<T> >
point (const point_basic<T>& x0);
// accessors:
const T& value() const;
const point_basic<T>& grad() const;
// algebra:
ad3_basic<T> operator+ () const;
ad3_basic<T> operator- () const;
ad3_basic<T> operator+ (const ad3_basic<T>& b) const;
ad3_basic<T> operator- (const ad3_basic<T>& b) const;
ad3_basic<T> operator* (const ad3_basic<T>& b) const;
ad3_basic<T> operator/ (const ad3_basic<T>& b) const;
ad3_basic<T>& operator += (const ad3_basic<T>& b);
ad3_basic<T>& operator -= (const ad3_basic<T>& b);
ad3_basic<T>& operator *= (const ad3_basic<T>& b);
ad3_basic<T>& operator /= (const ad3_basic<T>& b);
// data:
public:
T _v;
point_basic<T> _g;
};
typedef ad3_basic<Float> ad3;
// output:
template <class T>
std::ostream& operator<< (std::ostream& os, const ad3_basic<T>& a);
// get the underlying float
template<class T> struct float_traits<ad3_basic<T> > { typedef typename float_traits<T>::type type; };
// ----------------------------------------------------------------------------
// inlined
// ----------------------------------------------------------------------------
template <class T>
inline
ad3_basic<T>::ad3_basic()
: _v(0),
_g()
{
}
template <class T>
inline
ad3_basic<T>::ad3_basic(const ad3_basic<T>& a)
: _v(a._v),
_g(a._g)
{
}
template <class T>
inline
ad3_basic<T>&
ad3_basic<T>::operator= (const ad3_basic<T>& a)
{
_v = a._v;
_g = a._g;
return *this;
}
template <class T>
template <class U>
inline
ad3_basic<T>::ad3_basic(const U& v)
: _v(v),
_g()
{
}
template <class T>
inline
point_basic<ad3_basic<T> >
ad3_basic<T>::point (const point_basic<T>& x0)
{
point_basic<ad3_basic<T> > x;
for (size_type i = 0; i < 3; ++i) {
x[i]._v = x0[i];
x[i]._g[i] = 1;
}
return x;
}
template <class T>
inline
const T&
ad3_basic<T>::value () const
{
return _v;
}
template <class T>
inline
const point_basic<T>&
ad3_basic<T>::grad () const
{
return _g;
}
template <class T>
inline
std::ostream&
operator<< (std::ostream& os, const ad3_basic<T>& a)
{
os << a.value();
return os;
}
// ----------------------------------------------------------------------------
// op+-
// ----------------------------------------------------------------------------
template <class T>
inline
ad3_basic<T>
ad3_basic<T>::operator+ () const
{
ad3_basic<T> c;
c._v = _v;
c._g = _g;
return c;
}
template <class T>
inline
ad3_basic<T>
ad3_basic<T>::operator- () const
{
ad3_basic<T> c;
c._v = - _v;
c._g = - _g;
return c;
}
template <class T>
inline
ad3_basic<T>
ad3_basic<T>::operator+ (const ad3_basic<T>& b) const
{
ad3_basic<T> c;
c._v = _v + b._v;
c._g = _g + b._g;
return c;
}
template <class T>
inline
ad3_basic<T>
ad3_basic<T>::operator- (const ad3_basic<T>& b) const
{
ad3_basic<T> c;
c._v = _v - b._v;
c._g = _g - b._g;
return c;
}
template <class T, class U>
inline
typename
std::enable_if<
details::is_rheolef_arithmetic<U>::value
,ad3_basic<T>
>::type
operator+ (const U& a, const ad3_basic<T>& b)
{
ad3_basic<T> c;
c._v = a + b._v;
c._g = b._g;
return c;
}
template <class T, class U>
inline
typename
std::enable_if<
details::is_rheolef_arithmetic<U>::value
,ad3_basic<T>
>::type
operator+ (const ad3_basic<T>& a, const U& b)
{
ad3_basic<T> c;
c._v = a._v + b;
c._g = a._g;
return c;
}
template <class T, class U>
inline
typename
std::enable_if<
details::is_rheolef_arithmetic<U>::value
,ad3_basic<T>
>::type
operator- (const U& a, const ad3_basic<T>& b)
{
ad3_basic<T> c;
c._v = a - b._v;
c._g = - b._g;
return c;
}
template <class T, class U>
inline
typename
std::enable_if<
details::is_rheolef_arithmetic<U>::value
,ad3_basic<T>
>::type
operator- (const ad3_basic<T>& a, const U& b)
{
ad3_basic<T> c;
c._v = a._v - b;
c._g = a._g;
return c;
}
template <class T>
inline
ad3_basic<T>&
ad3_basic<T>::operator+= (const ad3_basic<T>& b)
{
*this = (*this) + b;
return *this;
}
template <class T, class U>
inline
typename
std::enable_if<
details::is_rheolef_arithmetic<U>::value
,ad3_basic<T>&
>::type
operator+= (ad3_basic<T>& a, const U& b)
{
a = a + b;
return a;
}
template <class T>
inline
ad3_basic<T>&
ad3_basic<T>::operator-= (const ad3_basic<T>& b)
{
*this = (*this) - b;
return *this;
}
template <class T, class U>
inline
typename
std::enable_if<
details::is_rheolef_arithmetic<U>::value
,ad3_basic<T>&
>::type
operator-= (ad3_basic<T>& a, const U& b)
{
a = a - b;
return a;
}
// ----------------------------------------------------------------------------
// op*
// ----------------------------------------------------------------------------
template <class T>
inline
ad3_basic<T>
ad3_basic<T>::operator* (const ad3_basic<T>& b) const
{
ad3_basic<T> c;
c._v = _v*b._v;
c._g = _v*b._g + _g*b._v;
return c;
}
template <class T, class U>
inline
typename
std::enable_if<
details::is_rheolef_arithmetic<U>::value
,ad3_basic<T>
>::type
operator* (const U& a, const ad3_basic<T>& b)
{
ad3_basic<T> c;
c._v = a*b._v;
c._g = a*b._g;
return c;
}
template <class T, class U>
inline
typename
std::enable_if<
details::is_rheolef_arithmetic<U>::value
,ad3_basic<T>
>::type
operator* (const ad3_basic<T>& a, const U& b)
{
ad3_basic<T> c;
c._v = a._v*b;
c._g = a._g*b;
return c;
}
template <class T>
inline
ad3_basic<T>&
ad3_basic<T>::operator*= (const ad3_basic<T>& b)
{
*this = (*this) * b;
return *this;
}
template <class T, class U>
inline
typename
std::enable_if<
details::is_rheolef_arithmetic<U>::value
,ad3_basic<T>&
>::type
operator*= (ad3_basic<T>& a, const U& b)
{
a = a * b;
return a;
}
// ----------------------------------------------------------------------------
// op/
// ----------------------------------------------------------------------------
template <class T>
inline
ad3_basic<T>
ad3_basic<T>::operator/ (const ad3_basic<T>& b) const
{
ad3_basic<T> c;
c._v = _v/b._v;
c._g = _g/b._v - _v*b._g/b._v;
return c;
}
template <class T, class U>
inline
typename
std::enable_if<
details::is_rheolef_arithmetic<U>::value
,ad3_basic<T>
>::type
operator/ (const U& a, const ad3_basic<T>& b)
{
ad3_basic<T> c;
c._v = a/b._v;
c._g = - a*b._g/b._v;
return c;
}
template <class T, class U>
inline
typename
std::enable_if<
details::is_rheolef_arithmetic<U>::value
,ad3_basic<T>
>::type
operator/ (const ad3_basic<T>& a, const U& b)
{
ad3_basic<T> c;
c._v = a._v/b;
c._g = a._g/b;
return c;
}
template <class T>
inline
ad3_basic<T>&
ad3_basic<T>::operator/= (const ad3_basic<T>& b)
{
*this = (*this) / b;
return *this;
}
template <class T, class U>
inline
typename
std::enable_if<
details::is_rheolef_arithmetic<U>::value
,ad3_basic<T>&
>::type
operator/= (ad3_basic<T>& a, const U& b)
{
a = a / b;
return a;
}
// -----------------------------------------------------------------------------
// comparators
// -----------------------------------------------------------------------------
#define _RHEOLEF_ad3_comparator(OP) \
template <class T> \
inline \
bool \
operator OP (const ad3_basic<T> &a, const ad3_basic<T> &b) \
{ \
return a._v OP b._v; \
} \
template <class T, class U> \
inline \
typename \
std::enable_if< \
details::is_rheolef_arithmetic<U>::value \
,bool \
>::type \
operator OP (const ad3_basic<T> &a, const U &b) \
{ \
return a._v OP b; \
} \
template <class T, class U> \
inline \
typename \
std::enable_if< \
details::is_rheolef_arithmetic<U>::value \
,bool \
>::type \
operator OP (const U &a, const ad3_basic<T> &b) \
{ \
return a OP b._v; \
}
_RHEOLEF_ad3_comparator(==)
_RHEOLEF_ad3_comparator(!=)
_RHEOLEF_ad3_comparator(>)
_RHEOLEF_ad3_comparator(>=)
_RHEOLEF_ad3_comparator(<)
_RHEOLEF_ad3_comparator(<=)
#undef _RHEOLEF_ad3_comparator
// ----------------------------------------------------------------------------
// std math functions: sqr, pow, exp...
// ----------------------------------------------------------------------------
// unary functions
// ----------------------------------------------------------------------------
template <class T>
inline
ad3_basic<T>
sqr (const ad3_basic<T>& a)
{
ad3_basic<T> c;
c._v = sqr(a._v);
c._g = (2*a._v)*a._g;
return c;
}
template <class T>
inline
ad3_basic<T>
fabs (const ad3_basic<T>& a)
{
ad3_basic<T> c;
c._v = fabs(a._v);
c._g = (a._v >= 0 ? 1 : -1)*a._g;
return c;
}
#ifdef TODO
// ----------------------------------------------------------------------------
// unary functions (cont.)
// ----------------------------------------------------------------------------
template <class T>
INLINE2 ad3_basic<T> exp (const ad3_basic<T>& a)
{
ad3_basic<T> c(exp(a.v));
c.touchg(a.gsize);
for (int i=0;i<a.gsize;i++) c.g[i]=a.g[i]*c.v;
return c;
}
template <class T>
INLINE2 ad3_basic<T> log (const ad3_basic<T>& a)
{
ad3_basic<T> c(log(a.v));
c.touchg(a.gsize);
for (int i=0;i<a.gsize;i++) c.g[i]=a.g[i]/a.v;
return c;
}
template <class T>
INLINE2 ad3_basic<T> sqrt (const ad3_basic<T>& a)
{
ad3_basic<T> c(sqrt(a.v));
T tmp(c.v*FADBAD_TWO);
c.touchg(a.gsize);
for (int i=0;i<a.gsize;i++) c.g[i]=a.g[i]/tmp;
return c;
}
template <class T>
INLINE2 ad3_basic<T> sin (const ad3_basic<T>& a)
{
ad3_basic<T> c(sin(a.v));
T tmp(cos(a.v));
c.touchg(a.gsize);
for (int i=0;i<a.gsize;i++) c.g[i]=a.g[i]*tmp;
return c;
}
template <class T>
INLINE2 ad3_basic<T> cos (const ad3_basic<T>& a)
{
ad3_basic<T> c(cos(a.v));
T tmp(-sin(a.v));
c.touchg(a.gsize);
for (int i=0;i<a.gsize;i++) c.g[i]=a.g[i]*tmp;
return c;
}
template <class T>
INLINE2 ad3_basic<T> tan (const ad3_basic<T>& a)
{
ad3_basic<T> c(tan(a.v));
c.touchg(a.gsize);
T tmp(FADBAD_ONE+_sqr(c.v));
for (int i=0;i<a.gsize;i++) c.g[i]=a.g[i]*tmp;
return c;
}
template <class T>
INLINE2 ad3_basic<T> asin (const ad3_basic<T>& a)
{
ad3_basic<T> c(asin(a.v));
c.touchg(a.gsize);
T tmp(FADBAD_ONE/sqrt(FADBAD_ONE-_sqr(a.v)));
for (int i=0;i<a.gsize;i++) c.g[i]=a.g[i]*tmp;
return c;
}
template <class T>
INLINE2 ad3_basic<T> acos (const ad3_basic<T>& a)
{
ad3_basic<T> c(acos(a.v));
c.touchg(a.gsize);
T tmp(-FADBAD_ONE/sqrt(FADBAD_ONE-_sqr(a.v)));
for (int i=0;i<a.gsize;i++) c.g[i]=a.g[i]*tmp;
return c;
}
template <class T>
INLINE2 ad3_basic<T> atan (const ad3_basic<T>& a)
{
ad3_basic<T> c(atan(a.v));
c.touchg(a.gsize);
T tmp(FADBAD_ONE/(FADBAD_ONE+_sqr(a.v)));
for (int i=0;i<a.gsize;i++) c.g[i]=a.g[i]*tmp;
return c;
}
// ----------------------------------------------------------------------------
// binary functions
// ----------------------------------------------------------------------------
template <class T>
INLINE2 ad3_basic<T> pow (const BaseType& a, const ad3_basic<T>& b)
{
ad3_basic<T> c(pow(a,b.v));
T tmp(c.v*log(a));
c.touchg(b.gsize);
for (int i=0;i<b.gsize;i++) c.g[i]=tmp*b.g[i];
return c;
}
template <class T>
INLINE2 ad3_basic<T> pow (const ad3_basic<T>& a,const BaseType& b)
{
ad3_basic<T> c(pow(a.v,b));
T tmp(b*pow(a.v,b-FADBAD_ONE));
c.touchg(a.gsize);
for (int i=0;i<a.gsize;i++) c.g[i]=tmp*a.g[i];
return c;
}
template <class T>
INLINE2 ad3_basic<T> pow (const ad3_basic<T>& a, const ad3_basic<T>& b)
{
if (a.gsize==0) return pow1(a.v,b);
if (b.gsize==0) return pow2(a,b.v);
ad3_basic<T> c(pow(a.v,b.v));
USER_ASSERT(a.gsize==b.gsize,"derivative vectors not of same size in pow");
T tmp(b.v*pow(a.v,b.v-FADBAD_ONE)),tmp1(c.v*log(a.v));
c.touchg(a.gsize);
for (int i=0;i<a.gsize;i++)
c.g[i]=tmp*a.g[i]+tmp1*b.g[i];
return c;
}
template <class T>
INLINE2 ad3_basic<T> pow (const ad3_basic<T>& a,int b)
{
ad3_basic<T> c(pow(a.v,b));
c.touchg(a.gsize);
T tmp(b*pow(a.v,b-1));
for (int i=0;i<a.gsize;i++) c.g[i]=a.g[i]*tmp;
return c;
}
#endif // TODO
} // namespace rheolef
#endif // _RHEO_AD_POINT_H
rheolef::ad3_basic::operator+=
ad3_basic< T > & operator+=(const ad3_basic< T > &b)
Definition: ad3.h:274
rheolef::ad3
ad3_basic< Float > ad3
Definition: ad3.h:95
rheolef::operator+=
std::enable_if< details::is_rheolef_arithmetic< U >::value,ad3_basic< T > & >::type operator+=(ad3_basic< T > &a, const U &b)
Definition: ad3.h:286
rheolef::float_traits
helper for std::complex<T>: get basic T type
Definition: Float.h:93
rheolef::ad3_basic::operator*=
ad3_basic< T > & operator*=(const ad3_basic< T > &b)
Definition: ad3.h:355
rheolef::ad3_basic::point
static point_basic< ad3_basic< T > > point(const point_basic< T > &x0)
Definition: ad3.h:141
rheolef::point_basic
Definition: point.h:87
rheolef::ad3_basic::operator-=
ad3_basic< T > & operator-=(const ad3_basic< T > &b)
Definition: ad3.h:294
rheolef::ad3_basic::operator*
ad3_basic< T > operator*(const ad3_basic< T > &b) const
Definition: ad3.h:317
rheolef::point_basic::size_type
size_t size_type
Definition: point.h:92
rheolef::operator*
csr< T, sequential > operator*(const T &lambda, const csr< T, sequential > &a)
Definition: csr.h:437
rheolef::ad3_basic::grad
const point_basic< T > & grad() const
Definition: ad3.h:160
rheolef::pow
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
Definition: space_mult.h:120
mkgeo_ball.c
c
Definition: mkgeo_ball.sh:153
rheolef::ad3_basic::value_type
T value_type
Definition: ad3.h:56
rheolef::ad3_basic
Definition: ad3.h:50
rheolef::operator/
dia< T, M > operator/(const T &lambda, const dia< T, M > &d)
Definition: dia.h:145
rheolef::ad3_basic::operator=
ad3_basic< T > & operator=(const ad3_basic< T > &)
Definition: ad3.h:124
rheolef::ad3_basic::size_type
point_basic< T >::size_type size_type
Definition: ad3.h:55
rheolef::float_traits< ad3_basic< T > >::type
float_traits< T >::type type
Definition: ad3.h:102
a
Definition: diffusion_isotropic.h:25
rheolef::type
rheolef::std type
rheolef::exp
tensor_basic< T > exp(const tensor_basic< T > &a, size_t d)
Definition: tensor-exp.cc:92
rheolef::operator-
csr< T, sequential > operator-(const csr< T, sequential > &a)
Definition: csr.h:447
rheolef::ad3_basic::operator/
ad3_basic< T > operator/(const ad3_basic< T > &b) const
Definition: ad3.h:378
rheolef::operator/=
std::enable_if< details::is_rheolef_arithmetic< U >::value,ad3_basic< T > & >::type operator/=(ad3_basic< T > &a, const U &b)
Definition: ad3.h:428
rheolef::ad3_basic::_g
point_basic< T > _g
Definition: ad3.h:93
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::float_traits::type
T type
Definition: Float.h:94
rheolef::ad3_basic::operator/=
ad3_basic< T > & operator/=(const ad3_basic< T > &b)
Definition: ad3.h:416
mkgeo_ball.b
b
Definition: mkgeo_ball.sh:152
mkgeo_ball.a
a
Definition: mkgeo_ball.sh:151
rheolef::ad3_basic::value
const T & value() const
Definition: ad3.h:153
rheolef::operator+
std::enable_if< details::is_rheolef_arithmetic< U >::value,ad3_basic< T >>::type operator+(const U &a, const ad3_basic< T > &b)
Definition: ad3.h:222
rheolef::operator-=
std::enable_if< details::is_rheolef_arithmetic< U >::value,ad3_basic< T > & >::type operator-=(ad3_basic< T > &a, const U &b)
Definition: ad3.h:306
rheolef::details::is_rheolef_arithmetic
Definition: Float.h:150
mkgeo_ball.tmp
tmp
Definition: mkgeo_ball.sh:380
rheolef::operator<<
std::ostream & operator<<(std::ostream &os, const catchmark &m)
Definition: catchmark.h:99
rheolef::operator*=
std::enable_if< details::is_rheolef_arithmetic< U >::value,ad3_basic< T > & >::type operator*=(ad3_basic< T > &a, const U &b)
Definition: ad3.h:367
rheolef::ad3_basic::operator-
ad3_basic< T > operator-() const
Definition: ad3.h:188
T
Expr1::float_type T
Definition: field_expr.h:218
rheolef::ad3_basic::ad3_basic
ad3_basic()
Definition: ad3.h:109
rheolef::ad3_basic::_v
T _v
Definition: ad3.h:92
rheolef::ad3_basic::operator+
ad3_basic< T > operator+() const
Definition: ad3.h:178