diff --git a/AutomaticDifferentiation_base.hpp b/AutomaticDifferentiation_base.hpp index 6b801a6..77121eb 100755 --- a/AutomaticDifferentiation_base.hpp +++ b/AutomaticDifferentiation_base.hpp @@ -14,15 +14,58 @@ #define __tplDualBase_t template /// Implementation of dual numbers for automatic differentiation. -/// This implementation uses vectors for b so that function gradients can be computed in one function call. -/// Set the index of every variable with the ::D(int i) function and call the function to be computed : f(x+__Dual_DualBase::D(0), y+__Dual_DualBase::D(1), z+__Dual_DualBase::D(2), ...) -/// Or use x.diff(int i) : __Dual_DualBase x(xf), y(yf), z(zf); x.diff(0); y.diff(1); z.diff(2); auto gradF = f(x,y,z).b; cout << gradF << endl; /// -/// - N is the number of derivatives that you want to compute simultaneously. Make sure to not mix __Dual_DualBase numbers with different N values. -/// - if bdynamic = false, the vectors are allocated on the stack. -/// - if bdynamic = true, the vectors are allocated dynamically on the heap. +/// Description +/// ----------- +/// Dual numbers serve to compute arbitrary function gradients efficiently and easily, with machine precision. +/// They can be used as a drop-in replacement for any of the base arithmetic types (double, float, etc). +/// The dual numbers are used in the so called forward automatic differentiation. /// -/// reference : http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.89.7749&rep=rep1&type=pdf +/// Here are some of the advantages of forward automatic differentiation compared to the backward method : +/// +/// - Contrary to the backward method, no graph must be computed, and the memory footprint of the forward method is greatly reduced. +/// - Contrary to popular belief, there *IS* a way to compute the whole gradient of a function using only a *single* function call (see examples below). +/// - The foward method is the *only one* that can be used to compute gradients of complicated numerical functions, such as the result of a numerical integration. +/// The backward method in these cases explodes the memory limit and crawls to a stop as it tries to record *all* the operations involved in the evaluation of the function. +/// +/// Template parameters and sub-classes +/// ----------------------------------- +/// There are 3 things to consider when using the Dual class to compute gradients : +/// +/// - What arithmetic type is going to be used. +/// - The number of variables with respect to which the gradient will be computed. +/// - Whether the vectors should be allocated dynamically on the heap or statically on the stack. +/// +/// Since the class defined as a template, any arithmetic type will do. Typically (but not limited to) : +/// +/// - float +/// - double +/// - long double +/// - quad float (quadmath) +/// - arbitrary precision (boost, gmp, mpfr, ...) +/// - fixed precision +/// - ... +/// +/// The number of variables with respect to which the gradient will be computed depends entirely on the problem to be solved. +/// +/// Finally, the type of memory management depends on how many variables will be part of the gradient computation : since the vector type used is provided by Eigen, +/// their recommendation should be followed : +/// +/// - Static for N < ~15 -> DualS<> +/// - Dynamic for N > ~15 -> DualD<> +/// +/// Typical use cases +/// ----------------- +/// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp} +/// Scalar xf = ..., yf = ..., zf = ...; +/// Dual x(xf), y(yf), z(zf), fx; +/// x.diff(0); y.diff(1); z.diff(2); +/// fx = f(x, y, z); // <--- A single function call ! +/// Scalar dfdx = fx.d(0), +/// dfdy = fx.d(1), +/// dfdz = fx.d(2); +/// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +/// Reference for the underlaying mathematics : http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.89.7749&rep=rep1&type=pdf template struct __Dual_DualBase { @@ -40,159 +83,183 @@ struct __Dual_DualBase #endif } - __Dual_DualBase(const Scalar & _a = Scalar()) - : a(_a) - { SetBToZero(); } + __Dual_DualBase(const Scalar & _a = Scalar()) + : a(_a) + { SetBToZero(); } - __Dual_DualBase(const Scalar & _a, const VectorT & _b) - : a(_a) - { - assert(_b.size() == N); - b = _b; - } + __Dual_DualBase(const Scalar & _a, const VectorT & _b) + : a(_a) + { + assert(_b.size() == N); + b = _b; + } - // __Dual_DualBase const& operator=(Scalar const& _a) - // { - // return *this; - // } + /// Use this function to set what variable is to be derived. + /// + /// The two following statements are *not exactly* equivalent, but produce the same effect (the two last cases are equivalent) : + /// + /// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp} + /// // Using Dual::D(int i) + /// Dual x(xf), y(yf), z(zf), fx; + /// fx = f(x+Dual::D(0), y+Dual::D(1), z+Dual::D(2)); + /// Scalar dfdx = fx.d(0); + /// Scalar dfdy = fx.d(1); + /// Scalar dfdz = fx.d(2); + /// + /// // Using diff(int i) before the function call + /// Dual x(xf), y(yf), z(zf), fx; + /// x.diff(0); y.diff(1); z.diff(2); + /// fx = f(x, y, z); + /// Scalar dfdx = fx.d(0); + /// Scalar dfdy = fx.d(1); + /// Scalar dfdz = fx.d(2); + /// + /// // Using diff(int i) directly during the function call + /// Dual x(xf), y(yf), z(zf), fx; + /// fx = f(x.diff(0), y.diff(1), z.diff(2)); + /// Scalar dfdx = fx.d(0); + /// Scalar dfdy = fx.d(1); + /// Scalar dfdz = fx.d(2); + /// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + /// + static __Dual_DualBase D(int i = 0) + { + assert(i >= 0); + assert(i < N); + __Dual_DualBase res(Scalar(0)); + res.b[i] = Scalar(1); + return res; + } - /// Use this function to set what variable is to be derived : x + __Dual_DualBase::d(i) - static __Dual_DualBase D(int i = 0) - { - assert(i >= 0); - assert(i < N); - __Dual_DualBase res(Scalar(0)); - res.b[i] = Scalar(1); - return res; - } + /// Use this function to set what variable is to be derived. Only one derivative can be toggled at once using this function. + /// For example, If x.b = {1 1 0}, after transformation y = f(x), y.b = {dy/dx, dy/dx, 0} + /// + /// Only one derivative should be selected per variable. + /// + /// In order to compute the gradient of a function, the following code can be used : + /// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp} + /// Dual x(xd), y(yd), z(zd), fxyz; + /// x.diff(0,3); // Set the first derivative to be that of x. + /// y.diff(1,3); // Set the first derivative to be that of y. + /// z.diff(2,3); // Set the first derivative to be that of z. + /// fxyz = f(x, y, z); // Evaluate the function to differentiate. + /// S dfdx = fxyz.d(0); // df/dx + /// S dfdy = fxyz.d(1); // df/dy + /// S dfdz = fxyz.d(2); // df/dz + /// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + __Dual_DualBase const& diff(int i = 0) + { + assert(i >= 0); + assert(i < N); + SetBToZero(); + b[i] = Scalar(1); + return *this; + } - /// Use this function to set what variable is to be derived. Only one derivative can be toggled at once using this function. - /// For example, If x.b = {1 1 0}, after transformation y = f(x), y.b = {dy/dx, dy/dx, 0} - /// Only one derivative should be selected per variable. - /// In order to compute the gradient of a function, the following code can be used : - /// __Dual_DualBase x(xd), y(yd), z(zd), fxyz; - /// x.diff(0,3); // Set the first derivative to be that of x. - /// y.diff(1,3); // Set the first derivative to be that of y. - /// z.diff(2,3); // Set the first derivative to be that of z. - /// fxyz = f(x, y, z); // Evaluate the function to differentiate. - /// S dfdx = fxyz.d(0); // df/dx - /// S dfdy = fxyz.d(1); // df/dy - /// S dfdz = fxyz.d(2); // df/dz - __Dual_DualBase const& diff(int i = 0) - { - assert(i >= 0); - assert(i < N); - SetBToZero(); - b[i] = Scalar(1); - return *this; - } + /// Returns a reference to the value + Scalar const& x() const { return a; } + Scalar & x() { return a; } - /// Returns a reference to the value - Scalar const& x() const { return a; } - Scalar & x() { return a; } + /// Returns a reference to the vector of infinitesimal parts + VectorT const& B() const { return b; } + VectorT & B() { return b; } - /// Returns a reference to the vector of infinitesimal parts - VectorT const& B() const { return b; } - VectorT & B() { return b; } + /// Returns the derivative value at index i. An assertion protects against i >= N. + Scalar const& d(int i) const + { + assert(i >= 0); + assert(i < N); + return b[i]; + } - /// Returns the derivative value at index i - Scalar const& d(int i) const - { - assert(i >= 0); - assert(i < N); - return b[i]; - } + /// Returns the derivative value at index i. An assertion protects against i >= N. + Scalar & d(int i) + { + assert(i >= 0); + assert(i < N); + return b[i]; + } - /// Returns the derivative value at index i - Scalar & d(int i) - { - assert(i >= 0); - assert(i < N); - return b[i]; - } + __Dual_DualBase & operator+=(const __Dual_DualBase & x) + { + a += x.a; + b += x.b; + return *this; + } - __Dual_DualBase & operator+=(const __Dual_DualBase & x) - { - a += x.a; - b += x.b; - // x() += x.x(); - // B() += x.B(); - return *this; - } + __Dual_DualBase & operator-=(const __Dual_DualBase & x) + { + a -= x.a; + b -= x.b; + return *this; + } - __Dual_DualBase & operator-=(const __Dual_DualBase & x) - { - a -= x.a; - b -= x.b; - return *this; - } + __Dual_DualBase & operator*=(const __Dual_DualBase & x) + { + b = a*x.b + b*x.a; + a *= x.a; + return *this; + } - __Dual_DualBase & operator*=(const __Dual_DualBase & x) - { - b = a*x.b + b*x.a; - a *= x.a; - return *this; - } + __Dual_DualBase & operator/=(const __Dual_DualBase & x) + { + b = (x.a*b - a*x.b)/(x.a*x.a); + a /= x.a; + return *this; + } - __Dual_DualBase & operator/=(const __Dual_DualBase & x) - { - b = (x.a*b - a*x.b)/(x.a*x.a); - a /= x.a; - return *this; - } + __Dual_DualBase & operator++() { return ((*this) += Scalar(1.)); }// ++x + __Dual_DualBase & operator--() { return ((*this) -= Scalar(1.)); }// --x - __Dual_DualBase & operator++() { return ((*this) += Scalar(1.)); }// ++x - __Dual_DualBase & operator--() { return ((*this) -= Scalar(1.)); }// --x + __Dual_DualBase operator++(int) { // x++ + __Dual_DualBase copy = *this; + (*this) += Scalar(1.); + return copy; + } - __Dual_DualBase operator++(int) { // x++ - __Dual_DualBase copy = *this; - (*this) += Scalar(1.); - return copy; - } + __Dual_DualBase operator--(int) { // x-- + __Dual_DualBase copy = *this; + (*this) -= Scalar(1.); + return copy; + } - __Dual_DualBase operator--(int) { // x-- - __Dual_DualBase copy = *this; - (*this) -= Scalar(1.); - return copy; - } + __Dual_DualBase operator+(const __Dual_DualBase & x) const { + __Dual_DualBase res(*this); + return (res += x); + } - __Dual_DualBase operator+(const __Dual_DualBase & x) const { - __Dual_DualBase res(*this); - return (res += x); - } + __Dual_DualBase operator-(const __Dual_DualBase & x) const { + __Dual_DualBase res(*this); + return (res -= x); + } + __Dual_DualBase operator*(const __Dual_DualBase & x) const + { + __Dual_DualBase res(*this); + return (res *= x); + } - __Dual_DualBase operator-(const __Dual_DualBase & x) const { - __Dual_DualBase res(*this); - return (res -= x); - } + __Dual_DualBase operator/(const __Dual_DualBase & x) const + { + __Dual_DualBase res(*this); + return (res /= x); + } - __Dual_DualBase operator*(const __Dual_DualBase & x) const - { - __Dual_DualBase res(*this); - return (res *= x); - } + __Dual_DualBase operator+(void) const { return (*this); }// +x + __Dual_DualBase operator-(void) const { return __Dual_DualBase(-a, -b); }// -x - __Dual_DualBase operator/(const __Dual_DualBase & x) const - { - __Dual_DualBase res(*this); - return (res /= x); - } + bool operator==(const __Dual_DualBase & x) const { return (a == x.a); } + bool operator!=(const __Dual_DualBase & x) const { return (a != x.a); } + bool operator<(const __Dual_DualBase & x) const { return (a < x.a); } + bool operator<=(const __Dual_DualBase & x) const { return (a <= x.a); } + bool operator>(const __Dual_DualBase & x) const { return (a > x.a); } + bool operator>=(const __Dual_DualBase & x) const { return (a >= x.a); } - __Dual_DualBase operator+(void) const { return (*this); }// +x - __Dual_DualBase operator-(void) const { return __Dual_DualBase(-a, -b); }// -x + /// Explicit conversion of the dual number to *ANY* type. Clearely, not every conversion makes sense. Use at your own risk. + template explicit operator T() const { return static_cast(a); } - bool operator==(const __Dual_DualBase & x) const { return (a == x.a); } - bool operator!=(const __Dual_DualBase & x) const { return (a != x.a); } - bool operator<(const __Dual_DualBase & x) const { return (a < x.a); } - bool operator<=(const __Dual_DualBase & x) const { return (a <= x.a); } - bool operator>(const __Dual_DualBase & x) const { return (a > x.a); } - bool operator>=(const __Dual_DualBase & x) const { return (a >= x.a); } - - template explicit operator T() const { return static_cast(a); } - - Scalar a; /// Real part - VectorT b; /// Infinitesimal parts + Scalar a; ///< Real part + VectorT b; ///< Infinitesimal parts }; template diff --git a/AutomaticDifferentiation_dynamic.hpp b/AutomaticDifferentiation_dynamic.hpp index 8668781..c256fe8 100755 --- a/AutomaticDifferentiation_dynamic.hpp +++ b/AutomaticDifferentiation_dynamic.hpp @@ -10,4 +10,35 @@ #include "AutomaticDifferentiation_base.hpp" +/// Macro to create a function object that returns the gradient of the function at X. +/// Suitable for functions that take an N-dimensional vector and return a scalar. +/// This version uses dynamically allocated arrays for everything (Dual and vectors to and from the function). +/// Designed to work with functions, lambdas, etc. +#define CREATE_GRAD_FUNCTION_OBJECT_N_1_D(Func, GradFuncName, N) \ +template \ +struct GradFuncName { \ + using ArrayOfScalar = Eigen::Array; \ + using ArrayOfDual = Eigen::Array, -1, 1>; \ + \ + ArrayOfScalar operator()(ArrayOfScalar const& x) { \ + Scalar fx; \ + ArrayOfScalar gradfx; \ + get_f_grad(x, fx, gradfx); \ + return gradfx; \ + } \ + void get_f_grad(ArrayOfScalar const& x, Scalar & fx, ArrayOfScalar & gradfx) { \ + ArrayOfDual X(N); \ + for(int i = 0 ; i < N ; i++) \ + { \ + X[i] = x[i]; \ + X[i].diff(i); \ + } \ + auto Y = Func>(X); \ + gradfx.resize(N); \ + for (int i = 0; i < N; i++) \ + gradfx[i] = Y.d(i); \ + fx = Y.x(); \ + } \ +} + #endif diff --git a/AutomaticDifferentiation_static.hpp b/AutomaticDifferentiation_static.hpp index b9e0987..c07c357 100755 --- a/AutomaticDifferentiation_static.hpp +++ b/AutomaticDifferentiation_static.hpp @@ -22,7 +22,7 @@ /// /// // Alternatively : /// Scalar fx, dfdx; -/// gradFct.get_f_grad(x, fx, dfdx); // Evaluation of the gradient of fct at x. +/// gradFct.get_f_grad(x, fx, dfdx); // Evaluation of value AND the gradient of fct at x. template struct GradFunc1 { @@ -35,11 +35,9 @@ struct GradFunc1 // Function that returns the 1st derivative of the function f at x. Scalar operator()(Scalar const& x) { - // differentiate using the dual number and return the .b component. - DualS X(x); - X.diff(0); - DualS Y = f(X); - return Y.d(0); + Scalar fx, gradfx; + get_f_grad(x, fx, gradfx); + return gradfx; } /// Function that returns both the function value and the gradient of f at x. Use this preferably over separate calls to f and to gradf. @@ -54,4 +52,55 @@ struct GradFunc1 } }; +/// Macro to create a function object that returns the gradient of the function at X. +/// Suitable for functions that take an scalar and return a scalar. +/// Designed to work with functions, lambdas, etc. +#define CREATE_GRAD_FUNCTION_OBJECT_1_1(Func, GradFuncName) \ +struct GradFuncName { \ + template \ + Scalar operator()(Scalar const& x) { \ + Scalar fx, gradfx; \ + get_f_grad(x, fx, gradfx); \ + return gradfx; \ + } \ + template \ + void get_f_grad(Scalar const& x, Scalar & fx, Scalar & gradfx) { \ + DualS X(x); \ + X.diff(0); \ + DualS Y = Func>(X); \ + fx = Y.x(); \ + gradfx = Y.d(0); \ + } \ +} + +/// Macro to create a function object that returns the gradient of the function at X. +/// Suitable for functions that take an N-dimensional vector and return a scalar. +/// This version uses statically allocated arrays for everything (Dual and vectors to and from the function). +/// Designed to work with functions, lambdas, etc. +#define CREATE_GRAD_FUNCTION_OBJECT_N_1_S(Func, GradFuncName, N) \ +template \ +struct GradFuncName { \ + using ArrayOfScalar = Eigen::Array; \ + using ArrayOfDual = Eigen::Array, N, 1>; \ + \ + ArrayOfScalar operator()(ArrayOfScalar const& x) { \ + Scalar fx; \ + ArrayOfScalar gradfx; \ + get_f_grad(x, fx, gradfx); \ + return gradfx; \ + } \ + void get_f_grad(ArrayOfScalar const& x, Scalar & fx, ArrayOfScalar & gradfx) { \ + ArrayOfDual X(N); \ + for(int i = 0 ; i < N ; i++) \ + { \ + X[i] = x[i]; \ + X[i].diff(i); \ + } \ + auto Y = Func>(X); \ + for (int i = 0; i < N; i++) \ + gradfx[i] = Y.d(i); \ + fx = Y.x(); \ + } \ +} + #endif diff --git a/Doxyfile b/Doxyfile index 012f224..a345eef 100755 --- a/Doxyfile +++ b/Doxyfile @@ -32,7 +32,7 @@ DOXYFILE_ENCODING = UTF-8 # title of most generated pages and in a few other places. # The default value is: My Project. -PROJECT_NAME = "C++ Quick Start Project" +PROJECT_NAME = "Automatic Differentiation" # The PROJECT_NUMBER tag can be used to enter a project or revision number. This # could be handy for archiving the generated documentation or if some version diff --git a/old/AutomaticDifferentiation_scalar.hpp b/old/AutomaticDifferentiation_scalar.hpp deleted file mode 100755 index 153964f..0000000 --- a/old/AutomaticDifferentiation_scalar.hpp +++ /dev/null @@ -1,411 +0,0 @@ -#ifndef DEF_AUTOMATIC_DIFFERENTIATION -#define DEF_AUTOMATIC_DIFFERENTIATION - -#include -#include - -/// Implementation of dual numbers for automatic differentiation -/// reference : http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.89.7749&rep=rep1&type=pdf -template -class Dual -{ - public: - Dual(const Scalar & _a, const Scalar & _b = Scalar(0.0)) - : a(_a), - b(_b) - {} - - static Dual d() { - return Dual(Scalar(0.), Scalar(1.)); - } - - Dual & operator+=(const Dual & x) { - a += x.a; - b += x.b; - return *this; - } - - Dual & operator-=(const Dual & x) { - a -= x.a; - b -= x.b; - return *this; - } - - Dual & operator*=(const Dual & x) { - b = a*x.b + b*x.a; - a *= x.a; - return *this; - } - - Dual & operator/=(const Dual & x) { - b = (x.a*b - a*x.b)/(x.a*x.a); - a /= x.a; - return *this; - } - - Dual & operator++() { // ++x - return ((*this) += Scalar(1.)); - } - - Dual & operator--() { // --x - return ((*this) -= Scalar(1.)); - } - - Dual operator++(int) { // x++ - Dual copy = *this; - (*this) += Scalar(1.); - return copy; - } - - Dual operator--(int) { // x-- - Dual copy = *this; - (*this) -= Scalar(1.); - return copy; - } - - Dual operator+(const Dual & x) const { - Dual res(*this); - return (res += x); - } - - Dual operator+(void) const { // +x - return (*this); - } - - Dual operator-(const Dual & x) const { - Dual res(*this); - return (res -= x); - } - - Dual operator-(void) const { // -x - return Dual(-a, -b); - } - - Dual operator*(const Dual & x) const { - Dual res(*this); - return (res *= x); - } - - Dual operator/(const Dual & x) const { - Dual res(*this); - return (res /= x); - } - - bool operator==(const Dual & x) const { - return (a == x.a); - } - - bool operator!=(const Dual & x) const { - return (a != x.a); - } - - bool operator<(const Dual & x) const { - return (a < x.a); - } - - bool operator<=(const Dual & x) const { - return (a <= x.a); - } - - bool operator>(const Dual & x) const { - return (a > x.a); - } - - bool operator>=(const Dual & x) const { - return (a >= x.a); - } - - Scalar a; /// Real part - Scalar b; /// Infinitesimal part -}; - -template -Dual operator+(A const& v, Dual const& x) { - return (Dual(v) + x); -} -template -Dual operator-(A const& v, Dual const& x) { - return (Dual(v) - x); -} -template -Dual operator*(A const& v, Dual const& x) { - return (Dual(v) * x); -} -template -Dual operator/(A const& v, Dual const& x) { - return (Dual(v) / x); -} - -// Basic mathematical functions for Scalar numbers - -// Trigonometric functions -template Scalar sec(const Scalar & x) { - return Scalar(1.)/cos(x); -} - -template Scalar cot(const Scalar & x) { - return cos(x)/sin(x); -} - -template Scalar csc(const Scalar & x) { - return Scalar(1.)/sin(x); -} - -// Inverse trigonometric functions -template Scalar asec(const Scalar & x) { - return acos(Scalar(1.)/x); -} - -template Scalar acot(const Scalar & x) { - return atan(Scalar(1.)/x); -} - -template Scalar acsc(const Scalar & x) { - return asin(Scalar(1.)/x); -} - -// Hyperbolic trigonometric functions -template Scalar sech(const Scalar & x) { - return Scalar(1.)/cosh(x); -} - -template Scalar coth(const Scalar & x) { - return cosh(x)/sinh(x); -} - -template Scalar csch(const Scalar & x) { - return Scalar(1.)/sinh(x); -} - -// Inverse hyperbolic trigonometric functions -template Scalar asech(const Scalar & x) { - return log((Scalar(1.) + sqrt(Scalar(1.) - x*x))/x); -} - -template Scalar acoth(const Scalar & x) { - return Scalar(0.5)*log((x + Scalar(1.))/(x - Scalar(1.))); -} - -template Scalar acsch(const Scalar & x) { - return (x >= Scalar(0.)) ? log((Scalar(1.) + sqrt(Scalar(1.) + x*x))/x) : log((Scalar(1.) - sqrt(Scalar(1.) + x*x))/x); -} - -// Other functions -template Scalar exp10(const Scalar & x) { - return exp(x*log(Scalar(10.))); -} - -template Scalar sign(const Scalar & x) { - return (x >= Scalar(0.)) ? ((x > Scalar(0.)) ? Scalar(1.) : Scalar(0.)) : Scalar(-1.); -} - -template Scalar heaviside(const Scalar & x) { - return Scalar(x >= Scalar(0.)); -} - -template Scalar abs(const Scalar & x) { - return (x >= Scalar(0.)) ? x : -x; -} - -// Basic mathematical functions for Dual numbers -// f(a + b*d) = f(a) + b*f'(a)*d - -// Trigonometric functions -template Dual cos(const Dual & x) { - return Dual(cos(x.a), -x.b*sin(x.a)); -} - -template Dual sin(const Dual & x) { - return Dual(sin(x.a), x.b*cos(x.a)); -} - -template Dual tan(const Dual & x) { - return Dual(tan(x.a), x.b*sec(x.a)*sec(x.a)); -} - -template Dual sec(const Dual & x) { - return Dual(sec(x.a), x.b*sec(x.a)*tan(x.a)); -} - -template Dual cot(const Dual & x) { - return Dual(cot(x.a), x.b*(-csc(x.a)*csc(x.a))); -} - -template Dual csc(const Dual & x) { - return Dual(csc(x.a), x.b*(-cot(x.a)*csc(x.a))); -} - -// Inverse trigonometric functions -template Dual acos(const Dual & x) { - return Dual(acos(x.a), x.b*(-Scalar(1.)/sqrt(Scalar(1.)-x.a*x.a))); -} - -template Dual asin(const Dual & x) { - return Dual(asin(x.a), x.b*(Scalar(1.)/sqrt(Scalar(1.)-x.a*x.a))); -} - -template Dual atan(const Dual & x) { - return Dual(atan(x.a), x.b*(Scalar(1.)/(x.a*x.a + Scalar(1.)))); -} - -template Dual asec(const Dual & x) { - return Dual(asec(x.a), x.b*(Scalar(1.)/(sqrt(Scalar(1.)-Scalar(1.)/(x.a*x.a))*(x.a*x.a)))); -} - -template Dual acot(const Dual & x) { - return Dual(acot(x.a), x.b*(-Scalar(1.)/((x.a*x.a)+Scalar(1.)))); -} - -template Dual acsc(const Dual & x) { - return Dual(acsc(x.a), x.b*(-Scalar(1.)/(sqrt(Scalar(1.)-Scalar(1.)/(x.a*x.a))*(x.a*x.a)))); -} - -// Hyperbolic trigonometric functions -template Dual cosh(const Dual & x) { - return Dual(cosh(x.a), x.b*sinh(x.a)); -} - -template Dual sinh(const Dual & x) { - return Dual(sinh(x.a), x.b*cosh(x.a)); -} - -template Dual tanh(const Dual & x) { - return Dual(tanh(x.a), x.b*sech(x.a)*sech(x.a)); -} - -template Dual sech(const Dual & x) { - return Dual(sech(x.a), x.b*(-sech(x.a)*tanh(x.a))); -} - -template Dual coth(const Dual & x) { - return Dual(coth(x.a), x.b*(-csch(x.a)*csch(x.a))); -} - -template Dual csch(const Dual & x) { - return Dual(csch(x.a), x.b*(-coth(x.a)*csch(x.a))); -} - -// Inverse hyperbolic trigonometric functions -template Dual acosh(const Dual & x) { - return Dual(acosh(x.a), x.b*(Scalar(1.)/sqrt((x.a*x.a)-Scalar(1.)))); -} - -template Dual asinh(const Dual & x) { - return Dual(asinh(x.a), x.b*(Scalar(1.)/sqrt((x.a*x.a)+Scalar(1.)))); -} - -template Dual atanh(const Dual & x) { - return Dual(atanh(x.a), x.b*(Scalar(1.)/(Scalar(1.)-(x.a*x.a)))); -} - -template Dual asech(const Dual & x) { - return Dual(asech(x.a), x.b*(Scalar(-1.)/(sqrt(Scalar(1.)/(x.a*x.a)-Scalar(1.))*(x.a*x.a)))); -} - -template Dual acoth(const Dual & x) { - return Dual(acoth(x.a), x.b*(-Scalar(1.)/((x.a*x.a)-Scalar(1.)))); -} - -template Dual acsch(const Dual & x) { - return Dual(acsch(x.a), x.b*(-Scalar(1.)/(sqrt(Scalar(1.)/(x.a*x.a)+Scalar(1.))*(x.a*x.a)))); -} - -// Exponential functions -template Dual exp(const Dual & x) { - return Dual(exp(x.a), x.b*exp(x.a)); -} - -template Dual log(const Dual & x) { - return Dual(log(x.a), x.b/x.a); -} - -template Dual exp10(const Dual & x) { - return Dual(exp10(x.a), x.b*(log(Scalar(10.))*exp10(x.a))); -} - -template Dual log10(const Dual & x) { - return Dual(log10(x.a), x.b/(log(Scalar(10.))*x.a)); -} - -template Dual exp2(const Dual & x) { - return Dual(exp2(x.a), x.b*(log(Scalar(2.))*exp2(x.a))); -} - -template Dual log2(const Dual & x) { - return Dual(log2(x.a), x.b/(log(Scalar(2.))*x.a)); -} - -template Dual pow(const Dual & x, const Dual & n) { - return exp(n*log(x)); -} - -// Other functions -template Dual sqrt(const Dual & x) { - return Dual(sqrt(x.a), x.b/(Scalar(2.)*sqrt(x.a))); -} - -template Dual sign(const Dual & x) { - return Dual(sign(x.a), Scalar(0.)); -} - -template Dual abs(const Dual & x) { - return Dual(abs(x.a), x.b*sign(x.a)); -} - -template Dual fabs(const Dual & x) { - return Dual(fabs(x.a), x.b*sign(x.a)); -} - -template Dual heaviside(const Dual & x) { - return Dual(heaviside(x.a), Scalar(0.)); -} - -template Dual floor(const Dual & x) { - return Dual(floor(x.a), Scalar(0.)); -} - -template Dual ceil(const Dual & x) { - return Dual(ceil(x.a), Scalar(0.)); -} - -template Dual round(const Dual & x) { - return Dual(round(x.a), Scalar(0.)); -} - -template std::ostream & operator<<(std::ostream & s, const Dual & x) -{ - return (s << x.a); -} - -#endif - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/old/test_AutomaticDifferentiation.cpp b/old/test_AutomaticDifferentiation.cpp deleted file mode 100755 index f3ddcce..0000000 --- a/old/test_AutomaticDifferentiation.cpp +++ /dev/null @@ -1,353 +0,0 @@ -#include -#include "AutomaticDifferentiation.hpp" - -#define TEST_EQ_TOL 1e-9 -#define EQ_DOUBLE(a, b) (fabs((a) - (b)) <= TEST_EQ_TOL) -#define EQ_DOUBLE_TOL(a, b, tol) (fabs((a) - (b)) <= (tol)) -#define CHECK_FUNCTION_ON_DUAL_DOUBLE(fct, x) CHECK( EQ_DOUBLE((fct(Dual((x))).a), fct((x))) ) - -#define TEST_DERIVATIVE_NUM(fct, x, DX_NUM_DIFF, tol) \ - CHECK( fabs(fct((Dual((x))+Dual::d())).b - (((fct((x)+DX_NUM_DIFF)) - (fct((x)-DX_NUM_DIFF)))/(2*DX_NUM_DIFF))) <= tol ) - -// Test function for nth derivative -template -Scalar testFunction1(const Scalar & x) -{ - return Scalar(1.)/atan(Scalar(1.) - pow(x, Scalar(2.))); -} - -template -Scalar dtestFunction1(const Scalar & x) -{ - return testFunction1(Dual(x) + Dual::d()).b; -} - -template -Scalar ddtestFunction1(const Scalar & x) -{ - return dtestFunction1(Dual(x) + Dual::d()).b; -} - -template -Scalar dtestFunction1_sym(const Scalar & x) -{ - return (Scalar(2.)*x)/((Scalar(1.) + pow(Scalar(1.) - pow(x,Scalar(2.)),Scalar(2.)))*pow(atan(Scalar(1.) - pow(x,Scalar(2.))),Scalar(2.))); -} - -template -Scalar ddtestFunction1_sym(const Scalar & x) -{ - return (Scalar(8.)*pow(x,Scalar(2.)))/(pow(Scalar(1.) + pow(Scalar(1.) - pow(x,Scalar(2.)),Scalar(2.)),Scalar(2.))* pow(atan(Scalar(1.) - pow(x,Scalar(2.))),Scalar(3.))) + (Scalar(8.)*pow(x,Scalar(2.))*(Scalar(1.) - pow(x,Scalar(2.))))/(pow(Scalar(1.) + pow(Scalar(1.) - pow(x,Scalar(2.)),Scalar(2.)),Scalar(2.))* pow(atan(Scalar(1.) - pow(x,Scalar(2.))),Scalar(2.))) + Scalar(2.)/((Scalar(1.) + pow(Scalar(1.) - pow(x,Scalar(2.)),Scalar(2.)))*pow(atan(Scalar(1.) - pow(x,Scalar(2.))),Scalar(2.))); -} - -// Length of vector from coordinates -template -Scalar f3(const Scalar & x, const Scalar & y, const Scalar & z) -{ - return sqrt(z*z+y*y+x*x); -} - -template -Scalar df3x(const Scalar & x, const Scalar & y, const Scalar & z) -{ - return x/sqrt(z*z+y*y+x*x); -} - -template -Scalar df3y(const Scalar & x, const Scalar & y, const Scalar & z) -{ - return y/sqrt(z*z+y*y+x*x); -} - -template -Scalar df3z(const Scalar & x, const Scalar & y, const Scalar & z) -{ - return z/sqrt(z*z+y*y+x*x); -} - -TEST_CASE( "Basic Dual class tests", "[basic]" ) { - - // For each section, these variables are anew - double x1 = 1.62, x2 = 0.62, x3 = 3.14, x4 = 2.71; - using D = Dual; - - SECTION( "Constructors" ) { - // REQUIRE( D().a == 0. ); - // REQUIRE( D().b == 0. ); - REQUIRE( D(x1).a == x1 ); - REQUIRE( D(x1).b == 0. ); - REQUIRE( D(x1, x2).a == x1 ); - REQUIRE( D(x1, x2).b == x2 ); - - // copy constructor - D X(x1, x2); - D Y(X); - REQUIRE( X.a == Y.a ); - REQUIRE( X.b == Y.b ); - - // d() function - REQUIRE( D::d().a == 0. ); - REQUIRE( D::d().b == 1. ); - } - - SECTION( "Comparison operators" ) { - D X(x1, x2); - D Y(x3, x4); - - // equal - REQUIRE( (X == X) ); - REQUIRE_FALSE( (X == Y) ); - - // different - REQUIRE_FALSE( (X != X) ); - REQUIRE( (X != Y) ); - - // lower than - REQUIRE( (X < Y) == (X.a < Y.a) ); - REQUIRE( (X < X) == (X.a < X.a) ); - REQUIRE( (X <= Y) == (X.a <= Y.a) ); - REQUIRE( (X <= X) == (X.a <= X.a) ); - - // greater than - REQUIRE( (X > Y) == (X.a > Y.a) ); - REQUIRE( (X > X) == (X.a > X.a) ); - REQUIRE( (X >= Y) == (X.a >= Y.a) ); - REQUIRE( (X >= X) == (X.a >= X.a) ); - } - - SECTION( "Operators for operations" ) { - D X(x1, x2); - D Y(x3, x4); - - REQUIRE( (X+Y).a == X.a+Y.a ); - REQUIRE( (X-Y).a == X.a-Y.a ); - REQUIRE( (X*Y).a == X.a*Y.a ); - REQUIRE( (X/Y).a == X.a/Y.a ); - - REQUIRE( (X+Y).b == X.b+Y.b ); - REQUIRE( (X-Y).b == X.b-Y.b ); - REQUIRE( (X*Y).b == X.a*Y.b+X.b*Y.a ); - REQUIRE( (X/Y).b == (Y.a*X.b - X.a*Y.b)/(Y.a*Y.a) ); - - // increment and decrement operators - double aValue = X.a; - REQUIRE( (++X).a == ++aValue); - REQUIRE( (--X).a == --aValue); - REQUIRE( (X++).a == aValue++); - REQUIRE( (X).a == aValue);// check that it was incremented properly - REQUIRE( (X--).a == aValue--); - REQUIRE( (X).a == aValue);// check that it was decremented properly - } -} - -TEST_CASE( "Scalar functions tests", "[scalarFunctions]" ) { - SECTION( "Scalar functions" ) { - REQUIRE( fabs(-5.) == 5. ); - REQUIRE( fabs(5.) == 5. ); - - CHECK( EQ_DOUBLE(cos(1.62), -0.04918382191417056) ); - CHECK( EQ_DOUBLE(sin(1.62), 0.998789743470524) ); - CHECK( EQ_DOUBLE(tan(1.62), -20.30728204110463) ); - CHECK( EQ_DOUBLE(sec(1.62), -20.33188884233264) ); - CHECK( EQ_DOUBLE(cot(1.62), -0.04924341908365026) ); - CHECK( EQ_DOUBLE(csc(1.62), 1.001211723025179) ); - - // inverse trigonometric functions - CHECK( EQ_DOUBLE(acos(0.62), 0.902053623592525) ); - CHECK( EQ_DOUBLE(asin(0.62), 0.6687427032023717) ); - CHECK( EQ_DOUBLE(atan(0.62), 0.5549957273385867) ); - CHECK( EQ_DOUBLE(asec(1.62), 0.905510600165641) ); - CHECK( EQ_DOUBLE(acot(0.62), 1.01580059945631) ); - CHECK( EQ_DOUBLE(acsc(1.62), 0.6652857266292561) ); - - // hyperbolic trigonometric functions - CHECK( EQ_DOUBLE(cosh(1.62), 2.625494507823741) ); - CHECK( EQ_DOUBLE(sinh(1.62), 2.427595808740127) ); - CHECK( EQ_DOUBLE(tanh(1.62), 0.924624218982788) ); - CHECK( EQ_DOUBLE(sech(1.62), 0.3808806291615117) ); - CHECK( EQ_DOUBLE(coth(1.62), 1.081520448491102) ); - CHECK( EQ_DOUBLE(csch(1.62), 0.4119301888723312) ); - - // inverse hyperbolic trigonometric functions - CHECK( EQ_DOUBLE(acosh(1.62), 1.062819127408777) ); - CHECK( EQ_DOUBLE(asinh(1.62), 1.259535895278778) ); - CHECK( EQ_DOUBLE(atanh(0.62), 0.7250050877529992) ); - CHECK( EQ_DOUBLE(asech(0.62), 1.057231115568124) ); - CHECK( EQ_DOUBLE(acoth(1.62), 0.7206050593580027) ); - CHECK( EQ_DOUBLE(acsch(1.62), 0.5835891509960214) ); - - // other functions - CHECK( EQ_DOUBLE(exp10(1.62), pow(10., 1.62)) ); - CHECK( EQ_DOUBLE(sign(1.62), 1.) ); - CHECK( EQ_DOUBLE(sign(-1.62), -1.) ); - CHECK( EQ_DOUBLE(sign(0.), 0.) ); - CHECK( EQ_DOUBLE(heaviside(-1.62), 0.) ); - CHECK( EQ_DOUBLE(heaviside(0.), 1.) ); - CHECK( EQ_DOUBLE(heaviside(1.62), 1.) ); - } -} - -TEST_CASE( "Functions on dual numbers", "[FunctionsOnDualNumbers]" ) { - SECTION( "Functions on Dual numbers" ) { - CHECK_FUNCTION_ON_DUAL_DOUBLE(cos, 1.62); - CHECK_FUNCTION_ON_DUAL_DOUBLE(sin, 1.62); - CHECK_FUNCTION_ON_DUAL_DOUBLE(tan, 1.62); - CHECK_FUNCTION_ON_DUAL_DOUBLE(sec, 1.62); - CHECK_FUNCTION_ON_DUAL_DOUBLE(cot, 1.62); - CHECK_FUNCTION_ON_DUAL_DOUBLE(csc, 1.62); - - // inverse trigonometric functions - CHECK_FUNCTION_ON_DUAL_DOUBLE(acos, 0.62); - CHECK_FUNCTION_ON_DUAL_DOUBLE(asin, 0.62); - CHECK_FUNCTION_ON_DUAL_DOUBLE(atan, 0.62); - CHECK_FUNCTION_ON_DUAL_DOUBLE(asec, 1.62); - CHECK_FUNCTION_ON_DUAL_DOUBLE(acot, 0.62); - CHECK_FUNCTION_ON_DUAL_DOUBLE(acsc, 1.62); - - // exponential functions - CHECK_FUNCTION_ON_DUAL_DOUBLE(exp, 1.62); - CHECK_FUNCTION_ON_DUAL_DOUBLE(log, 1.62); - CHECK_FUNCTION_ON_DUAL_DOUBLE(exp10, 1.62); - CHECK_FUNCTION_ON_DUAL_DOUBLE(log10, 1.62); - CHECK_FUNCTION_ON_DUAL_DOUBLE(exp2, 1.62); - CHECK_FUNCTION_ON_DUAL_DOUBLE(log2, 1.62); - CHECK_FUNCTION_ON_DUAL_DOUBLE(sqrt, 1.62); - - // hyperbolic trigonometric functions - CHECK_FUNCTION_ON_DUAL_DOUBLE(cosh, 1.62); - CHECK_FUNCTION_ON_DUAL_DOUBLE(sinh, 1.62); - CHECK_FUNCTION_ON_DUAL_DOUBLE(tanh, 1.62); - CHECK_FUNCTION_ON_DUAL_DOUBLE(sech, 1.62); - CHECK_FUNCTION_ON_DUAL_DOUBLE(coth, 1.62); - CHECK_FUNCTION_ON_DUAL_DOUBLE(csch, 1.62); - - // inverse hyperbolic trigonometric functions - CHECK_FUNCTION_ON_DUAL_DOUBLE(acosh, 1.62); - CHECK_FUNCTION_ON_DUAL_DOUBLE(asinh, 1.62); - CHECK_FUNCTION_ON_DUAL_DOUBLE(atanh, 0.62); - CHECK_FUNCTION_ON_DUAL_DOUBLE(asech, 0.62); - CHECK_FUNCTION_ON_DUAL_DOUBLE(acoth, 1.62); - CHECK_FUNCTION_ON_DUAL_DOUBLE(acsch, 1.62); - - // other functions - CHECK_FUNCTION_ON_DUAL_DOUBLE(sign, -1.62); - CHECK_FUNCTION_ON_DUAL_DOUBLE(sign, 0.00); - CHECK_FUNCTION_ON_DUAL_DOUBLE(sign, 1.62); - - CHECK( abs(Dual(-1.62)) == fabs(-1.62) ); - CHECK( fabs(Dual(-1.62)) == fabs(-1.62) ); - CHECK_FUNCTION_ON_DUAL_DOUBLE(fabs, -1.62); - CHECK_FUNCTION_ON_DUAL_DOUBLE(fabs, 1.62); - - CHECK_FUNCTION_ON_DUAL_DOUBLE(heaviside, -1.62); - CHECK_FUNCTION_ON_DUAL_DOUBLE(heaviside, 0.00); - CHECK_FUNCTION_ON_DUAL_DOUBLE(heaviside, 1.62); - - CHECK_FUNCTION_ON_DUAL_DOUBLE(floor, 1.2); - CHECK_FUNCTION_ON_DUAL_DOUBLE(ceil, 1.2); - CHECK_FUNCTION_ON_DUAL_DOUBLE(round, 1.2); - } - - SECTION( "Function derivatives checked numerically" ) { - double x = 0.5, x2 = 1.5, x3 = -x, dx = 1e-6, tol = 1e-9; - TEST_DERIVATIVE_NUM(cos, x, dx, tol); - TEST_DERIVATIVE_NUM(sin, x, dx, tol); - TEST_DERIVATIVE_NUM(tan, x, dx, tol); - TEST_DERIVATIVE_NUM(sec, x, dx, tol); - TEST_DERIVATIVE_NUM(cot, x, dx, tol); - TEST_DERIVATIVE_NUM(csc, x, dx, tol); - - TEST_DERIVATIVE_NUM(acos, x, dx, tol); - TEST_DERIVATIVE_NUM(asin, x, dx, tol); - TEST_DERIVATIVE_NUM(atan, x, dx, tol); - TEST_DERIVATIVE_NUM(asec, x2, dx, tol); - TEST_DERIVATIVE_NUM(acot, x, dx, tol); - TEST_DERIVATIVE_NUM(acsc, x2, dx, tol); - - TEST_DERIVATIVE_NUM(cosh, x, dx, tol); - TEST_DERIVATIVE_NUM(sinh, x, dx, tol); - TEST_DERIVATIVE_NUM(tanh, x, dx, tol); - TEST_DERIVATIVE_NUM(sech, x, dx, tol); - TEST_DERIVATIVE_NUM(coth, x, dx, tol); - TEST_DERIVATIVE_NUM(csch, x, dx, tol); - - TEST_DERIVATIVE_NUM(acosh, x2, dx, tol); - TEST_DERIVATIVE_NUM(asinh, x, dx, tol); - TEST_DERIVATIVE_NUM(atanh, x, dx, tol); - TEST_DERIVATIVE_NUM(asech, x, dx, tol); - TEST_DERIVATIVE_NUM(acoth, x2, dx, tol); - TEST_DERIVATIVE_NUM(acsch, x, dx, tol); - - TEST_DERIVATIVE_NUM(exp, x2, dx, tol); - TEST_DERIVATIVE_NUM(log, x, dx, tol); - TEST_DERIVATIVE_NUM(exp10, x, dx, tol); - TEST_DERIVATIVE_NUM(log10, x, dx, tol); - TEST_DERIVATIVE_NUM(exp2, x2, dx, tol); - TEST_DERIVATIVE_NUM(log2, x, dx, tol); - - TEST_DERIVATIVE_NUM(sqrt, x2, dx, tol); - TEST_DERIVATIVE_NUM(sign, x3, dx, tol); - TEST_DERIVATIVE_NUM(sign, x, dx, tol); - TEST_DERIVATIVE_NUM(fabs, x3, dx, tol); - TEST_DERIVATIVE_NUM(fabs, x, dx, tol); - TEST_DERIVATIVE_NUM(heaviside, x3, dx, tol); - TEST_DERIVATIVE_NUM(heaviside, x, dx, tol); - TEST_DERIVATIVE_NUM(floor, 1.6, dx, tol); - TEST_DERIVATIVE_NUM(ceil, 1.6, dx, tol); - TEST_DERIVATIVE_NUM(round, 1.6, dx, tol); - } -} - -TEST_CASE( "Nth derivative", "[NthDerivative]" ) { - double x = 0.7; - double dx = 1e-6; - double fx = testFunction1(x); - double dfdx = (testFunction1(x+dx) - testFunction1(x-dx))/(2*dx); - double d2fdx2 = (testFunction1(x-dx) - 2*fx + testFunction1(x+dx))/(dx*dx); - - CHECK( testFunction1(Dual(x)).a == fx ); - CHECK( EQ_DOUBLE_TOL(dtestFunction1_sym(x), dfdx, 1e-6) ); - CHECK( EQ_DOUBLE_TOL(ddtestFunction1_sym(x), d2fdx2, 1e-4) ); - - CHECK( testFunction1(Dual(x) + Dual::d()).b == dtestFunction1_sym(x) ); - CHECK( dtestFunction1(x) == dtestFunction1_sym(x) ); - CHECK( EQ_DOUBLE(ddtestFunction1(x), ddtestFunction1_sym(x)) ); -} - -TEST_CASE( "Partial derivatives", "[PartialDerivatives]" ) { - double x = 0.7, y = -2., z = 1.5; - - CHECK( f3(Dual(x), Dual(y), Dual(z)) == f3(x,y,z) ); - CHECK( f3(Dual(x)+Dual::d(), Dual(y), Dual(z)).b == df3x(x,y,z) ); - CHECK( f3(Dual(x), Dual(y)+Dual::d(), Dual(z)).b == df3y(x,y,z) ); - CHECK( f3(Dual(x), Dual(y), Dual(z)+Dual::d()).b == df3z(x,y,z) ); -} - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/old/test_AutomaticDifferentiation_manual.cpp b/old/test_AutomaticDifferentiation_manual.cpp deleted file mode 100755 index 1f4fa93..0000000 --- a/old/test_AutomaticDifferentiation_manual.cpp +++ /dev/null @@ -1,503 +0,0 @@ -#include "AutomaticDifferentiation.hpp" - -#include -#include -#include -#include -using std::cout; -using std::cout; -using std::setw; -#define PRINT_VAR(x) std::cout << #x << "\t= " << std::setprecision(16) << (x) << std::endl -#define PRINT_DUAL(x) std::cout << #x << "\t= " << std::fixed << std::setprecision(4) << std::setw(10) << (x).a << ", " << std::setw(10) << (x).b << std::endl - -#define TEST_EQ_TOL 1e-9 -// #define TEST_FUNCTION_ON_DUAL_DOUBLE(fct, x) assert(fct(Dual(x)).a == fct(x)) -// #define TEST_FUNCTION_ON_DUAL_DOUBLE(fct, x) assert(abs((fct(Dual((x))).a) - (fct((x)))) <= TEST_EQ_TOL) -#define TEST_EQ_DOUBLE(a, b) assert(abs((a) - (b)) <= TEST_EQ_TOL) - -#define TEST_FUNCTION_ON_DUAL_DOUBLE(fct, x) \ - if(abs((fct(Dual((x))).a) - fct((x))) > TEST_EQ_TOL) { \ - std::cerr << "Assertion failed at " << __FILE__ << ":" << __LINE__ << " : " << #fct << ">(" << (x) << ") != " << #fct << "(" << (x) << ")" << "\n";\ - std::cerr << "Got " << (fct(Dual(x)).a) << " ; expected " << (fct((x))) << "\n";\ - exit(1);\ - } - -template void print_T() { std::cout << __PRETTY_FUNCTION__ << '\n'; } - -template -Scalar f1(const Scalar & x) -{ - return Scalar(5.)*x*x*x + Scalar(3.)*x*x - Scalar(2.)*x + Scalar(4.); -} - -template -Scalar df1(const Scalar & x) -{ - return Scalar(15.)*x*x + Scalar(6.)*x - Scalar(2.); -} - -template -Scalar ddf1(const Scalar & x) -{ - return Scalar(30.)*x + Scalar(6.); -} - -template -Scalar g1(Scalar x) { - return f1(Dual(x) + Dual::d()).b; -} - -template -Scalar h1(Scalar x) { - return g1(Dual(x) + Dual::d()).b; -} - -template D f2(D x) { - return (x + D(2.0)) * (x + D(1.0)); -} - -template -Scalar f3(const Scalar & x, const Scalar & y, const Scalar & z) -{ - return sqrt(z*z+y*y+x*x); -} - -template -Scalar df3x(const Scalar & x, const Scalar & y, const Scalar & z) -{ - return x/sqrt(z*z+y*y+x*x); -} - -template -Scalar df3y(const Scalar & x, const Scalar & y, const Scalar & z) -{ - return y/sqrt(z*z+y*y+x*x); -} - -template -Scalar df3z(const Scalar & x, const Scalar & y, const Scalar & z) -{ - return z/sqrt(z*z+y*y+x*x); -} - -void test_basic(); -void test_scalar_functions(); -void test_derivative_all(); -void test_derivative_pow(); -void test_derivative_simple(); -void test_derivative_simple_2(); -void test_derivative_simple_3(); -void test_derivative_nested(); - -int main() -{ - test_scalar_functions(); - test_basic(); - test_derivative_all(); - test_derivative_pow(); - // test_derivative_simple(); - // test_derivative_simple_2(); - // test_derivative_simple_3(); - // test_derivative_nested(); - - return 0; -} - -void test_basic() -{ - cout << "\ntest_basic()\n"; - using D = Dual; - cout.precision(16); - - double x = 2; - double y = 5; - - D X(x), Y(y), Z(x, y); - - D W = Z; - - assert(X.a == x); - assert(X.b == 0.); - assert(Y.a == y); - assert(Y.b == 0.); - assert((X+Y).a == (x+y)); - assert((X-Y).a == (x-y)); - assert((X*Y).a == (x*y)); - assert((X/Y).a == (x/y)); - assert(-X.a == -x); - - assert(D(1., 2.).a == 1.); - assert(D(1., 2.).b == 2.); - assert(W.a == Z.a); - assert(W.b == Z.b); - - assert((D(1., 2.)+D(4., 7.)).a == 5.); - assert((D(1., 2.)+D(4., 7.)).b == 9.); - - // test all the value returned by non-linear functions - assert(heaviside(-1.) == 0.); - assert(heaviside(0.) == 1.); - assert(heaviside(1.) == 1.); - - assert(sign(-10.) == -1.); - assert(sign(0.) == 0.); - assert(sign(10.) == 1.); - - PRINT_VAR(abs(-10.)); - PRINT_VAR(abs(10.)); - PRINT_VAR(abs(Dual(-10.)).a); - PRINT_VAR(abs(Dual(10.)).a); - PRINT_VAR(abs(Dual(-1.62)).a); - PRINT_VAR(abs(-1.62)); - - PRINT_VAR(exp10(-3.)); - PRINT_VAR(exp10(-2.)); - PRINT_VAR(exp10(-1.)); - PRINT_VAR(exp10(0.)); - PRINT_VAR(exp10(1.)); - PRINT_VAR(exp10(2.)); - PRINT_VAR(exp10(3.)); - - x = -1.5; - PRINT_VAR((x >= (0.)) ? pow((10.), x) : (1.)/pow((10.), -x)); - - PRINT_VAR(pow(10., -3.)); - PRINT_VAR(pow(10., 3.)); - PRINT_VAR(1./pow(10., 3.)); - - PRINT_VAR(exp2(-3.)); - PRINT_VAR(exp2(3.)); - - TEST_EQ_DOUBLE(exp10(-3.), 1./exp10(3.)); - TEST_EQ_DOUBLE(exp10(-3.), 0.001); - TEST_EQ_DOUBLE(exp10( 3.), 1000.); - - PRINT_VAR(atanh(0.62)); - PRINT_VAR(atanh(Dual(0.62)).a); - - PRINT_VAR(acsc(Dual(1.62)).a); - PRINT_VAR(acsc(1.62)); - - TEST_EQ_DOUBLE(pow(Dual(1.62), Dual(1.5)).a, pow(1.62, 1.5)); - - // trigonometric functions - TEST_FUNCTION_ON_DUAL_DOUBLE(cos, 1.62); - TEST_FUNCTION_ON_DUAL_DOUBLE(sin, 1.62); - TEST_FUNCTION_ON_DUAL_DOUBLE(tan, 1.62); - TEST_FUNCTION_ON_DUAL_DOUBLE(sec, 1.62); - TEST_FUNCTION_ON_DUAL_DOUBLE(cot, 1.62); - TEST_FUNCTION_ON_DUAL_DOUBLE(csc, 1.62); - - // inverse trigonometric functions - TEST_FUNCTION_ON_DUAL_DOUBLE(acos, 0.62); - TEST_FUNCTION_ON_DUAL_DOUBLE(asin, 0.62); - TEST_FUNCTION_ON_DUAL_DOUBLE(atan, 0.62); - TEST_FUNCTION_ON_DUAL_DOUBLE(asec, 1.62); - TEST_FUNCTION_ON_DUAL_DOUBLE(acot, 0.62); - TEST_FUNCTION_ON_DUAL_DOUBLE(acsc, 1.62); - - // exponential functions - TEST_FUNCTION_ON_DUAL_DOUBLE(exp, 1.62); - TEST_FUNCTION_ON_DUAL_DOUBLE(log, 1.62); - TEST_FUNCTION_ON_DUAL_DOUBLE(exp10, 1.62); - TEST_FUNCTION_ON_DUAL_DOUBLE(log10, 1.62); - TEST_FUNCTION_ON_DUAL_DOUBLE(exp2, 1.62); - TEST_FUNCTION_ON_DUAL_DOUBLE(log2, 1.62); - TEST_FUNCTION_ON_DUAL_DOUBLE(sqrt, 1.62); - - // hyperbolic trigonometric functions - TEST_FUNCTION_ON_DUAL_DOUBLE(cosh, 1.62); - TEST_FUNCTION_ON_DUAL_DOUBLE(sinh, 1.62); - TEST_FUNCTION_ON_DUAL_DOUBLE(tanh, 1.62); - TEST_FUNCTION_ON_DUAL_DOUBLE(sech, 1.62); - TEST_FUNCTION_ON_DUAL_DOUBLE(coth, 1.62); - TEST_FUNCTION_ON_DUAL_DOUBLE(csch, 1.62); - - // inverse hyperbolic trigonometric functions - TEST_FUNCTION_ON_DUAL_DOUBLE(acosh, 1.62); - TEST_FUNCTION_ON_DUAL_DOUBLE(asinh, 1.62); - TEST_FUNCTION_ON_DUAL_DOUBLE(atanh, 0.62); - TEST_FUNCTION_ON_DUAL_DOUBLE(asech, 1.62); - TEST_FUNCTION_ON_DUAL_DOUBLE(acoth, 1.62); - TEST_FUNCTION_ON_DUAL_DOUBLE(acsch, 1.62); - - // other functions - TEST_FUNCTION_ON_DUAL_DOUBLE(sign, -1.62); - TEST_FUNCTION_ON_DUAL_DOUBLE(sign, 0.00); - TEST_FUNCTION_ON_DUAL_DOUBLE(sign, 1.62); - - TEST_FUNCTION_ON_DUAL_DOUBLE(abs, -1.62); - TEST_FUNCTION_ON_DUAL_DOUBLE(abs, 1.62); - - TEST_FUNCTION_ON_DUAL_DOUBLE(heaviside, -1.62); - TEST_FUNCTION_ON_DUAL_DOUBLE(heaviside, 0.00); - TEST_FUNCTION_ON_DUAL_DOUBLE(heaviside, 1.62); - - TEST_FUNCTION_ON_DUAL_DOUBLE(floor, 1.2); - TEST_FUNCTION_ON_DUAL_DOUBLE(ceil, 1.2); - TEST_FUNCTION_ON_DUAL_DOUBLE(round, 1.2); -} - -void test_scalar_functions() -{ - // test basic scalar functions numerically with values from mathematica - cout << "\ntest_scalar_functions()\n"; - - // trigonometric functions - //{-0.04918382191417056,0.998789743470524,-20.30728204110463,-20.33188884233264,-0.04924341908365026,1.001211723025179} - TEST_EQ_DOUBLE(cos(1.62), -0.04918382191417056); - TEST_EQ_DOUBLE(sin(1.62), 0.998789743470524); - TEST_EQ_DOUBLE(tan(1.62), -20.30728204110463); - TEST_EQ_DOUBLE(sec(1.62), -20.33188884233264); - TEST_EQ_DOUBLE(cot(1.62), -0.04924341908365026); - TEST_EQ_DOUBLE(csc(1.62), 1.001211723025179); - - // inverse trigonometric functions - // {0.902053623592525,0.6687427032023717,0.5549957273385867,0.905510600165641,1.01580059945631,0.6652857266292561} - TEST_EQ_DOUBLE(acos(0.62), 0.902053623592525); - TEST_EQ_DOUBLE(asin(0.62), 0.6687427032023717); - TEST_EQ_DOUBLE(atan(0.62), 0.5549957273385867); - TEST_EQ_DOUBLE(asec(1.62), 0.905510600165641); - TEST_EQ_DOUBLE(acot(0.62), 1.01580059945631); - TEST_EQ_DOUBLE(acsc(1.62), 0.6652857266292561); - - // hyperbolic trigonometric functions - // {2.625494507823741,2.427595808740127,0.924624218982788,0.3808806291615117,1.081520448491102,0.4119301888723312} - TEST_EQ_DOUBLE(cosh(1.62), 2.625494507823741); - TEST_EQ_DOUBLE(sinh(1.62), 2.427595808740127); - TEST_EQ_DOUBLE(tanh(1.62), 0.924624218982788); - TEST_EQ_DOUBLE(sech(1.62), 0.3808806291615117); - TEST_EQ_DOUBLE(coth(1.62), 1.081520448491102); - TEST_EQ_DOUBLE(csch(1.62), 0.4119301888723312); - - // inverse hyperbolic trigonometric functions - // {1.062819127408777,1.259535895278778,0.7250050877529992,1.057231115568124,0.7206050593580027,0.5835891509960214} - TEST_EQ_DOUBLE(acosh(1.62), 1.062819127408777); - TEST_EQ_DOUBLE(asinh(1.62), 1.259535895278778); - TEST_EQ_DOUBLE(atanh(0.62), 0.7250050877529992); - TEST_EQ_DOUBLE(asech(0.62), 1.057231115568124); - TEST_EQ_DOUBLE(acoth(1.62), 0.7206050593580027); - TEST_EQ_DOUBLE(acsch(1.62), 0.5835891509960214); -} - -#define TEST_DERIVATIVE_NUM(fct, x, DX_NUM_DIFF, tol); \ -{\ - double dfdx = fct((Dual((x))+Dual::d())).b;\ - double dfdx_num = ((fct((x)+DX_NUM_DIFF)) - (fct((x)-DX_NUM_DIFF)))/(2*DX_NUM_DIFF);\ - bool ok = (abs(dfdx - dfdx_num) <= tol) ? true : false;\ - cout << setw(10) << #fct << "(" << (x) << ") : " << setw(25) << dfdx << " " << setw(25) << dfdx_num << " " << setw(25) << (dfdx-dfdx_num) << "\t" << ok << "\n";\ -} - -void test_derivative_all() -{ - cout << "\ntest_derivative_all()\n"; - // test all derivatives numerically and check that they add up - double x = 0.5, x2 = 1.5, x3 = -x, dx = 1e-6, tol = 1e-9; - TEST_DERIVATIVE_NUM(cos, x, dx, tol); - TEST_DERIVATIVE_NUM(sin, x, dx, tol); - TEST_DERIVATIVE_NUM(tan, x, dx, tol); - TEST_DERIVATIVE_NUM(sec, x, dx, tol); - TEST_DERIVATIVE_NUM(cot, x, dx, tol); - TEST_DERIVATIVE_NUM(csc, x, dx, tol); - - TEST_DERIVATIVE_NUM(acos, x, dx, tol); - TEST_DERIVATIVE_NUM(asin, x, dx, tol); - TEST_DERIVATIVE_NUM(atan, x, dx, tol); - TEST_DERIVATIVE_NUM(asec, x2, dx, tol); - TEST_DERIVATIVE_NUM(acot, x, dx, tol); - TEST_DERIVATIVE_NUM(acsc, x2, dx, tol); - - TEST_DERIVATIVE_NUM(cosh, x, dx, tol); - TEST_DERIVATIVE_NUM(sinh, x, dx, tol); - TEST_DERIVATIVE_NUM(tanh, x, dx, tol); - TEST_DERIVATIVE_NUM(sech, x, dx, tol); - TEST_DERIVATIVE_NUM(coth, x, dx, tol); - TEST_DERIVATIVE_NUM(csch, x, dx, tol); - - TEST_DERIVATIVE_NUM(acosh, x2, dx, tol); - TEST_DERIVATIVE_NUM(asinh, x, dx, tol); - TEST_DERIVATIVE_NUM(atanh, x, dx, tol); - TEST_DERIVATIVE_NUM(asech, x, dx, tol); - TEST_DERIVATIVE_NUM(acoth, x2, dx, tol); - TEST_DERIVATIVE_NUM(acsch, x, dx, tol); - - TEST_DERIVATIVE_NUM(exp, x2, dx, tol); - TEST_DERIVATIVE_NUM(log, x, dx, tol); - TEST_DERIVATIVE_NUM(exp10, x, dx, tol); - TEST_DERIVATIVE_NUM(log10, x, dx, tol); - TEST_DERIVATIVE_NUM(exp2, x2, dx, tol); - TEST_DERIVATIVE_NUM(log2, x, dx, tol); - - TEST_DERIVATIVE_NUM(sqrt, x2, dx, tol); - TEST_DERIVATIVE_NUM(sign, x3, dx, tol); - TEST_DERIVATIVE_NUM(sign, x, dx, tol); - TEST_DERIVATIVE_NUM(abs, x3, dx, tol); - TEST_DERIVATIVE_NUM(abs, x, dx, tol); - TEST_DERIVATIVE_NUM(heaviside, x3, dx, tol); - TEST_DERIVATIVE_NUM(heaviside, x, dx, tol); - TEST_DERIVATIVE_NUM(floor, 1.6, dx, tol); - TEST_DERIVATIVE_NUM(ceil, 1.6, dx, tol); - TEST_DERIVATIVE_NUM(round, 1.6, dx, tol); -} -void test_derivative_pow() -{ - // test the derivatives of the power function - cout << "\ntest_derivative_pow()\n"; - using D = Dual; - - double a = 1.5, b = 5.4; - double c = pow(a, b); - double dcda = b*pow(a, b-1); - double dcdb = pow(a, b)*log(a); - D A(a), B(b); - - PRINT_VAR(a); - PRINT_VAR(b); - PRINT_VAR(c); - PRINT_VAR(dcda); - PRINT_VAR(dcdb); - - PRINT_DUAL(pow(A, B)); - - PRINT_DUAL(pow(A+D::d(), B)); - PRINT_DUAL(pow(A, B+D::d())); - - double dcda_AD = pow(A+D::d(), B).b; - double dcdb_AD = pow(A, B+D::d()).b; - - TEST_EQ_DOUBLE(pow(A, B).a, c); - TEST_EQ_DOUBLE(dcda_AD, dcda); - TEST_EQ_DOUBLE(dcdb_AD, dcdb); -} - -void test_derivative_simple() -{ - cout << "\ntest_derivative_simple()\n"; - using D = Dual; - - D d(0., 1.); - D x = 3.5; - D y = f1(x+d); - D dy = df1(x); - D ddy = ddf1(x); - - PRINT_DUAL(x); - PRINT_DUAL(x+d); - PRINT_DUAL(y); - PRINT_DUAL(dy); - PRINT_VAR(g1(x)); - PRINT_VAR(h1(x.a)); - PRINT_VAR(ddy); - - assert(y.b == dy.a); -} - -void test_derivative_simple_2() -{ - cout << "\ntest_derivative_simple_2()\n"; - using D = Dual; - - D d(0., 1.); - D x = 3.; - D x2 = x+d; - D y = f2(x); - D y2 = f2(x2); - D y3 = f2(D(3.0)+d); - - PRINT_DUAL(x); - PRINT_DUAL(x+d); - PRINT_DUAL(y); - PRINT_DUAL(y2); - PRINT_DUAL(y3); - - assert(y.a == 20.); - assert(y2.a == 20.); - assert(y3.a == 20.); - assert(y2.b == 9.); - assert(y3.b == 9.); -} - -void test_derivative_simple_3() -{ - // partial derivatives using scalar implementation - cout << "\ntest_derivative_simple_3()\n"; - using D = Dual; - - D x(2.), y(4.), z(-7.); - D L = f3(x, y, z); - - PRINT_DUAL(L); - PRINT_VAR(f3(x+D::d(), y, z).b); - PRINT_VAR(f3(x, y+D::d(), z).b); - PRINT_VAR(f3(x, y, z+D::d()).b); - - PRINT_VAR(df3x(x.a, y.a, z.a)); - PRINT_VAR(df3y(x.a, y.a, z.a)); - PRINT_VAR(df3z(x.a, y.a, z.a)); -} - -template -Scalar testFunctionNesting(const Scalar & x) -{ - return Scalar(1.)/atan(Scalar(1.) - pow(x, Scalar(2.))); -} - -template -struct TestFunctionNestingFunctor -{ - Scalar operator()(const Scalar & x) - { - return Scalar(1.)/atan(Scalar(1.) - pow(x, Scalar(2.))); - } -}; - -/* template -std::vector computeNfirstDerivatives(FunctorType & functor, const Scalar & x, const int & n, int depth = 0) -{ - // compute the n first derivatives using the recursive, nested approach - std::vector derivatives; - - if(depth < n) - derivatives.push_back(computeNfirstDerivatives(functor, Dual(x) + Dual::d(), n, ++depth)[0].b); - - return derivatives; -} - -void test_derivative_nested() -{ - // nested function to compute the first nth derivatives using scalar implementation - cout << "\ntest_derivative_nested()\n"; - using D = Dual; - - double x = 0.7; - TestFunctionNestingFunctor ffun; - - std::vector derivatives = computeNfirstDerivatives >(ffun, x, 1); - - -} */ - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/old/test_AutomaticDifferentiation_vector.cpp b/old/test_AutomaticDifferentiation_vector.cpp deleted file mode 100755 index f39dcae..0000000 --- a/old/test_AutomaticDifferentiation_vector.cpp +++ /dev/null @@ -1,473 +0,0 @@ -#include "AutomaticDifferentiationVector.hpp" - -#include -#include -#include -using std::cout; -using std::cout; -using std::setw; -using std::abs; -#define PRINT_VAR(x) std::cout << #x << "\t= " << std::setprecision(16) << (x) << std::endl -#define PRINT_DUALVECTOR(x) std::cout << #x << "\t= " << std::fixed << std::setprecision(4) << std::setw(10) << (x).a << ", " << std::setw(10) << (x).b << std::endl - -#define TEST_EQ_TOL 1e-9 -// #define TEST_FUNCTION_ON_DualVector_DOUBLE(fct, x) assert(fct(DualVector(x)).a == fct(x)) -// #define TEST_FUNCTION_ON_DualVector_DOUBLE(fct, x) assert(abs((fct(DualVector((x))).a) - (fct((x)))) <= TEST_EQ_TOL) -#define TEST_EQ_DOUBLE(a, b) assert(std::abs((a) - (b)) <= TEST_EQ_TOL) - -#define TEST_FUNCTION_ON_DualVector_DOUBLE(fct, x) \ - if(fabs((fct(DualVector((x))).a) - fct((x))) > TEST_EQ_TOL) { \ - std::cerr << "Assertion failed at " << __FILE__ << ":" << __LINE__ << " : " << #fct << ">(" << (x) << ") != " << #fct << "(" << (x) << ")" << "\n";\ - std::cerr << "Got " << (fct(DualVector(x)).a) << " ; expected " << (fct((x))) << "\n";\ - exit(1);\ - } - -template void print_T() { std::cout << __PRETTY_FUNCTION__ << '\n'; } - -template -Scalar f1(const Scalar & x) -{ - return Scalar(5.)*x*x*x + Scalar(3.)*x*x - Scalar(2.)*x + Scalar(4.); -} - -template -Scalar df1(const Scalar & x) -{ - return Scalar(15.)*x*x + Scalar(6.)*x - Scalar(2.); -} - -template -Scalar ddf1(const Scalar & x) -{ - return Scalar(30.)*x + Scalar(6.); -} - -template -Vector g1(Scalar x) { - return f1(DualVector(x) + DualVector::d()).b; -} - -template -Vector h1(Scalar x) { - return g1(DualVector(x) + DualVector::d()).b; -} - -template D f2(D x) { - return (x + D(2.0)) * (x + D(1.0)); -} - -template -Scalar f3(const Scalar & x, const Scalar & y, const Scalar & z) -{ - return sqrt(z*z+y*y+x*x); -} - -template -Scalar df3x(const Scalar & x, const Scalar & y, const Scalar & z) -{ - return x/sqrt(z*z+y*y+x*x); -} - -template -Scalar df3y(const Scalar & x, const Scalar & y, const Scalar & z) -{ - return y/sqrt(z*z+y*y+x*x); -} - -template -Scalar df3z(const Scalar & x, const Scalar & y, const Scalar & z) -{ - return z/sqrt(z*z+y*y+x*x); -} - -void test_basic(); -void test_scalar_functions(); -void test_derivative_all(); -void test_derivative_pow(); -void test_derivative_simple(); -void test_derivative_simple_2(); -void test_derivative_simple_3(); - -int main() -{ - test_scalar_functions(); - test_basic(); - test_derivative_all(); - test_derivative_pow(); - // test_derivative_simple(); - test_derivative_simple_2(); - test_derivative_simple_3(); - - return 0; -} - -void test_basic() -{ - cout << "\ntest_basic()\n"; - using D = DualVector; - cout.precision(16); - - DualVector::VectorT zeroVec(0., 3); - - double x = 2; - double y = 5; - - D X(x), Y(y); - - assert(X.a == x); - for(size_t i = 0 ; i < 3 ; i++) - assert(X.b[i] == 0.); - assert(Y.a == y); - for(size_t i = 0 ; i < 3 ; i++) - assert(Y.b[i] == 0.); - assert((X+Y).a == (x+y)); - assert((X-Y).a == (x-y)); - assert((X*Y).a == (x*y)); - assert((X/Y).a == (x/y)); - assert(-X.a == -x); - - PRINT_DUALVECTOR(X+Y); - PRINT_DUALVECTOR(X-Y); - PRINT_DUALVECTOR(X*Y); - PRINT_DUALVECTOR(X/Y); - - assert(D(1., 2.).a == 1.); - // assert(D(1., 2.).b == 2.); - - assert((D(1., 2.)+D(4., 7.)).a == 5.); - // assert((D(1., 2.)+D(4., 7.)).b == 9.); - - // test all the value returned by non-linear functions - assert(heaviside(-1.) == 0.); - assert(heaviside(0.) == 1.); - assert(heaviside(1.) == 1.); - - assert(sign(-10.) == -1.); - assert(sign(0.) == 0.); - assert(sign(10.) == 1.); - - PRINT_VAR(abs(-10.)); - PRINT_VAR(abs(10.)); - PRINT_VAR(abs(DualVector(-10.)).a); - PRINT_VAR(abs(DualVector(10.)).a); - PRINT_VAR(abs(DualVector(-1.62)).a); - PRINT_VAR(abs(-1.62)); - - PRINT_VAR(exp10(-3.)); - PRINT_VAR(exp10(-2.)); - PRINT_VAR(exp10(-1.)); - PRINT_VAR(exp10(0.)); - PRINT_VAR(exp10(1.)); - PRINT_VAR(exp10(2.)); - PRINT_VAR(exp10(3.)); - - x = -1.5; - PRINT_VAR((x >= (0.)) ? pow((10.), x) : (1.)/pow((10.), -x)); - - PRINT_VAR(pow(10., -3.)); - PRINT_VAR(pow(10., 3.)); - PRINT_VAR(1./pow(10., 3.)); - - PRINT_VAR(exp2(-3.)); - PRINT_VAR(exp2(3.)); - - TEST_EQ_DOUBLE(exp10(-3.), 1./exp10(3.)); - TEST_EQ_DOUBLE(exp10(-3.), 0.001); - TEST_EQ_DOUBLE(exp10( 3.), 1000.); - - PRINT_VAR(atanh(0.62)); - PRINT_VAR(atanh(DualVector(0.62)).a); - - PRINT_VAR(acsc(DualVector(1.62)).a); - PRINT_VAR(acsc(1.62)); - - TEST_EQ_DOUBLE(pow(DualVector(1.62), DualVector(1.5)).a, pow(1.62, 1.5)); - - // trigonometric functions - TEST_FUNCTION_ON_DualVector_DOUBLE(cos, 1.62); - TEST_FUNCTION_ON_DualVector_DOUBLE(sin, 1.62); - TEST_FUNCTION_ON_DualVector_DOUBLE(tan, 1.62); - TEST_FUNCTION_ON_DualVector_DOUBLE(sec, 1.62); - TEST_FUNCTION_ON_DualVector_DOUBLE(cot, 1.62); - TEST_FUNCTION_ON_DualVector_DOUBLE(csc, 1.62); - - // inverse trigonometric functions - TEST_FUNCTION_ON_DualVector_DOUBLE(acos, 0.62); - TEST_FUNCTION_ON_DualVector_DOUBLE(asin, 0.62); - TEST_FUNCTION_ON_DualVector_DOUBLE(atan, 0.62); - TEST_FUNCTION_ON_DualVector_DOUBLE(asec, 1.62); - TEST_FUNCTION_ON_DualVector_DOUBLE(acot, 0.62); - TEST_FUNCTION_ON_DualVector_DOUBLE(acsc, 1.62); - - // exponential functions - TEST_FUNCTION_ON_DualVector_DOUBLE(exp, 1.62); - TEST_FUNCTION_ON_DualVector_DOUBLE(log, 1.62); - TEST_FUNCTION_ON_DualVector_DOUBLE(exp10, 1.62); - TEST_FUNCTION_ON_DualVector_DOUBLE(log10, 1.62); - TEST_FUNCTION_ON_DualVector_DOUBLE(exp2, 1.62); - TEST_FUNCTION_ON_DualVector_DOUBLE(log2, 1.62); - TEST_FUNCTION_ON_DualVector_DOUBLE(sqrt, 1.62); - - // hyperbolic trigonometric functions - TEST_FUNCTION_ON_DualVector_DOUBLE(cosh, 1.62); - TEST_FUNCTION_ON_DualVector_DOUBLE(sinh, 1.62); - TEST_FUNCTION_ON_DualVector_DOUBLE(tanh, 1.62); - TEST_FUNCTION_ON_DualVector_DOUBLE(sech, 1.62); - TEST_FUNCTION_ON_DualVector_DOUBLE(coth, 1.62); - TEST_FUNCTION_ON_DualVector_DOUBLE(csch, 1.62); - - // inverse hyperbolic trigonometric functions - TEST_FUNCTION_ON_DualVector_DOUBLE(acosh, 1.62); - TEST_FUNCTION_ON_DualVector_DOUBLE(asinh, 1.62); - TEST_FUNCTION_ON_DualVector_DOUBLE(atanh, 0.62); - TEST_FUNCTION_ON_DualVector_DOUBLE(asech, 1.62); - TEST_FUNCTION_ON_DualVector_DOUBLE(acoth, 1.62); - TEST_FUNCTION_ON_DualVector_DOUBLE(acsch, 1.62); - - // other functions - TEST_FUNCTION_ON_DualVector_DOUBLE(sign, -1.62); - TEST_FUNCTION_ON_DualVector_DOUBLE(sign, 0.00); - TEST_FUNCTION_ON_DualVector_DOUBLE(sign, 1.62); - - TEST_FUNCTION_ON_DualVector_DOUBLE(abs, -1.62); - TEST_FUNCTION_ON_DualVector_DOUBLE(abs, 1.62); - - TEST_FUNCTION_ON_DualVector_DOUBLE(heaviside, -1.62); - TEST_FUNCTION_ON_DualVector_DOUBLE(heaviside, 0.00); - TEST_FUNCTION_ON_DualVector_DOUBLE(heaviside, 1.62); - - TEST_FUNCTION_ON_DualVector_DOUBLE(floor, 1.2); - TEST_FUNCTION_ON_DualVector_DOUBLE(ceil, 1.2); - TEST_FUNCTION_ON_DualVector_DOUBLE(round, 1.2); -} - -void test_scalar_functions() -{ - // test basic scalar functions numerically with values from mathematica - cout << "\ntest_scalar_functions()\n"; - - // trigonometric functions - //{-0.04918382191417056,0.998789743470524,-20.30728204110463,-20.33188884233264,-0.04924341908365026,1.001211723025179} - TEST_EQ_DOUBLE(cos(1.62), -0.04918382191417056); - TEST_EQ_DOUBLE(sin(1.62), 0.998789743470524); - TEST_EQ_DOUBLE(tan(1.62), -20.30728204110463); - TEST_EQ_DOUBLE(sec(1.62), -20.33188884233264); - TEST_EQ_DOUBLE(cot(1.62), -0.04924341908365026); - TEST_EQ_DOUBLE(csc(1.62), 1.001211723025179); - - // inverse trigonometric functions - // {0.902053623592525,0.6687427032023717,0.5549957273385867,0.905510600165641,1.01580059945631,0.6652857266292561} - TEST_EQ_DOUBLE(acos(0.62), 0.902053623592525); - TEST_EQ_DOUBLE(asin(0.62), 0.6687427032023717); - TEST_EQ_DOUBLE(atan(0.62), 0.5549957273385867); - TEST_EQ_DOUBLE(asec(1.62), 0.905510600165641); - TEST_EQ_DOUBLE(acot(0.62), 1.01580059945631); - TEST_EQ_DOUBLE(acsc(1.62), 0.6652857266292561); - - // hyperbolic trigonometric functions - // {2.625494507823741,2.427595808740127,0.924624218982788,0.3808806291615117,1.081520448491102,0.4119301888723312} - TEST_EQ_DOUBLE(cosh(1.62), 2.625494507823741); - TEST_EQ_DOUBLE(sinh(1.62), 2.427595808740127); - TEST_EQ_DOUBLE(tanh(1.62), 0.924624218982788); - TEST_EQ_DOUBLE(sech(1.62), 0.3808806291615117); - TEST_EQ_DOUBLE(coth(1.62), 1.081520448491102); - TEST_EQ_DOUBLE(csch(1.62), 0.4119301888723312); - - // inverse hyperbolic trigonometric functions - // {1.062819127408777,1.259535895278778,0.7250050877529992,1.057231115568124,0.7206050593580027,0.5835891509960214} - TEST_EQ_DOUBLE(acosh(1.62), 1.062819127408777); - TEST_EQ_DOUBLE(asinh(1.62), 1.259535895278778); - TEST_EQ_DOUBLE(atanh(0.62), 0.7250050877529992); - TEST_EQ_DOUBLE(asech(0.62), 1.057231115568124); - TEST_EQ_DOUBLE(acoth(1.62), 0.7206050593580027); - TEST_EQ_DOUBLE(acsch(1.62), 0.5835891509960214); -} - -#define TEST_DERIVATIVE_NUM(fct, x, DX_NUM_DIFF, tol) \ -{\ - double dfdx = fct((DualVector((x))+DualVector::d())).b[0];\ - double dfdx_num = double(-1.)/double(60.) * fct(x-3*dx)\ - + double( 3.)/double(20.) * fct(x-2*dx)\ - + double(-3.)/double(4. ) * fct(x-1*dx)\ - + double( 3.)/double(4. ) * fct(x+1*dx)\ - + double(-3.)/double(20.) * fct(x+2*dx)\ - + double( 1.)/double(60.) * fct(x+3*dx);\ - dfdx_num /= dx;\ - bool ok = (std::abs(dfdx - dfdx_num) <= tol) ? true : false;\ - cout << setw(10) << #fct << "(" << (x) << ") : " << setw(25) << dfdx << " " << setw(25) << dfdx_num << " " << setw(25) << (dfdx-dfdx_num) << "\t" << ok << "\n";\ - assert(ok);\ -} - -void test_derivative_all() -{ - cout << "\ntest_derivative_all()\n"; - // test all derivatives numerically and check that they add up - double x = 0.5, x2 = 1.5, x3 = -x, dx = 1e-3, tol = 1e-9; - TEST_DERIVATIVE_NUM(cos, x, dx, tol); - TEST_DERIVATIVE_NUM(sin, x, dx, tol); - TEST_DERIVATIVE_NUM(tan, x, dx, tol); - TEST_DERIVATIVE_NUM(sec, x, dx, tol); - TEST_DERIVATIVE_NUM(cot, x, dx, tol); - TEST_DERIVATIVE_NUM(csc, x, dx, tol); - - TEST_DERIVATIVE_NUM(acos, x, dx, tol); - TEST_DERIVATIVE_NUM(asin, x, dx, tol); - TEST_DERIVATIVE_NUM(atan, x, dx, tol); - TEST_DERIVATIVE_NUM(asec, x2, dx, tol); - TEST_DERIVATIVE_NUM(acot, x, dx, tol); - TEST_DERIVATIVE_NUM(acsc, x2, dx, tol); - - TEST_DERIVATIVE_NUM(cosh, x, dx, tol); - TEST_DERIVATIVE_NUM(sinh, x, dx, tol); - TEST_DERIVATIVE_NUM(tanh, x, dx, tol); - TEST_DERIVATIVE_NUM(sech, x, dx, tol); - TEST_DERIVATIVE_NUM(coth, x, dx, tol); - TEST_DERIVATIVE_NUM(csch, x, dx, tol); - - TEST_DERIVATIVE_NUM(acosh, x2, dx, tol); - TEST_DERIVATIVE_NUM(asinh, x, dx, tol); - TEST_DERIVATIVE_NUM(atanh, x, dx, tol); - TEST_DERIVATIVE_NUM(asech, x, dx, tol); - TEST_DERIVATIVE_NUM(acoth, x2, dx, tol); - TEST_DERIVATIVE_NUM(acsch, x, dx, tol); - - TEST_DERIVATIVE_NUM(exp, x2, dx, tol); - TEST_DERIVATIVE_NUM(log, x, dx, tol); - TEST_DERIVATIVE_NUM(exp10, x, dx, tol); - TEST_DERIVATIVE_NUM(log10, x, dx, tol); - TEST_DERIVATIVE_NUM(exp2, x2, dx, tol); - TEST_DERIVATIVE_NUM(log2, x, dx, tol); - - TEST_DERIVATIVE_NUM(sqrt, x2, dx, tol); - TEST_DERIVATIVE_NUM(sign, x3, dx, tol); - TEST_DERIVATIVE_NUM(sign, x, dx, tol); - TEST_DERIVATIVE_NUM(abs, x3, dx, tol); - TEST_DERIVATIVE_NUM(abs, x, dx, tol); - TEST_DERIVATIVE_NUM(heaviside, x3, dx, tol); - TEST_DERIVATIVE_NUM(heaviside, x, dx, tol); - TEST_DERIVATIVE_NUM(floor, 1.6, dx, tol); - TEST_DERIVATIVE_NUM(ceil, 1.6, dx, tol); - TEST_DERIVATIVE_NUM(round, 1.6, dx, tol); -} - -void test_derivative_pow() -{ - // test the derivatives of the power function - cout << "\ntest_derivative_pow()\n"; - using D = DualVector; - - double a = 1.5, b = 5.4; - double c = pow(a, b); - double dcda = b*pow(a, b-1); - double dcdb = pow(a, b)*log(a); - D A(a), B(b); - - PRINT_VAR(a); - PRINT_VAR(b); - PRINT_VAR(c); - PRINT_VAR(dcda); - PRINT_VAR(dcdb); - - PRINT_DUALVECTOR(pow(A, B)); - - PRINT_DUALVECTOR(pow(A+D::d(0), B+D::d(1))); - - D A1(A+D::d(0)), B1(B+D::d(1)), C(0.); - PRINT_DUALVECTOR(exp(B1*log(A1))); - - C = pow(A1, B1); - TEST_EQ_DOUBLE(C.a , c); - TEST_EQ_DOUBLE(C.b[0], dcda); - TEST_EQ_DOUBLE(C.b[1], dcdb); -} - -void test_derivative_simple() -{ - cout << "\ntest_derivative_simple()\n"; - /*using D = DualVector; - - D d = D::d(); - D x = 3.5; - D y = f1(x+D::d()); - D dy = df1(x); - D ddy = ddf1(x); - - PRINT_DUALVECTOR(x); - PRINT_DUALVECTOR(x+d); - PRINT_DUALVECTOR(y); - PRINT_DUALVECTOR(dy); - PRINT_VAR((g1(x))); - // PRINT_VAR(h1(x.a)); - PRINT_VAR(ddy); - - assert(y.b[0] == dy.a);//*/ -} - -void test_derivative_simple_2() -{ - cout << "\ntest_derivative_simple_2()\n"; - using D = DualVector; - - D d(0., 1.); - D x = 3.; - D x2 = x+d; - D y = f2(x); - D y2 = f2(x2); - D y3 = f2(D(3.0)+d); - - PRINT_DUALVECTOR(x); - PRINT_DUALVECTOR(x+d); - PRINT_DUALVECTOR(y); - PRINT_DUALVECTOR(y2); - PRINT_DUALVECTOR(y3); - - assert(y.a == 20.); - assert(y2.a == 20.); - assert(y3.a == 20.); - assert(y2.b[0] == 9.); - assert(y3.b[0] == 9.); -} - -void test_derivative_simple_3() -{ - // partial derivatives using scalar implementation - cout << "\ntest_derivative_simple_3()\n"; - using D = DualVector; - - D x(2.), y(4.), z(-7.); - D L = f3(x, y, z); - - PRINT_DUALVECTOR(L); - PRINT_VAR(f3(x+D::d(0), y+D::d(1), z+D::d(2)).b); - - PRINT_VAR(df3x(x.a, y.a, z.a)); - PRINT_VAR(df3y(x.a, y.a, z.a)); - PRINT_VAR(df3z(x.a, y.a, z.a)); -} - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/test/dual_test.cpp b/test/dual_test.cpp index 1fe8668..20602e7 100755 --- a/test/dual_test.cpp +++ b/test/dual_test.cpp @@ -547,6 +547,3 @@ TEST_CASE( "Dual numbers : Jacobian using Eigen arrays and Dual directly", "[dua for(int j = 0 ; j < 3 ; j++) CHECK(fabs(xcart[i].d(j) - Jxcart(i,j)) < tol); } - -// CHECK() -// REQUIRE() diff --git a/tests/functor_tests/Makefile b/tests/functor_tests/Makefile new file mode 100755 index 0000000..5171130 --- /dev/null +++ b/tests/functor_tests/Makefile @@ -0,0 +1,37 @@ +# Declaration of variables +C = clang +COMMON_FLAGS = -Wall -MMD +C_FLAGS = $(COMMON_FLAGS) +CC = clang++ +CC_FLAGS = $(COMMON_FLAGS) -std=c++17 -O3#-O0 -g +LD_FLAGS = +INCLUDES = -I/usr/include/eigen3 -I../.. + +# File names +EXEC = run +CSOURCES = $(wildcard *.c) +COBJECTS = $(CSOURCES:.c=.o) +SOURCES = $(wildcard *.cpp) +OBJECTS = $(SOURCES:.cpp=.o) + +# Main target +$(EXEC): $(COBJECTS) $(OBJECTS) + $(CC) $(COBJECTS) $(OBJECTS) -o $(EXEC) $(LD_FLAGS) + +# To obtain object files +%.o: %.cpp + $(CC) $(INCLUDES) $(CC_FLAGS) -o $@ -c $< + +# To obtain object files +%.o: %.c + $(C) $(INCLUDES) $(C_FLAGS) -o $@ -c $< + +-include $(SOURCES:%.cpp=%.d) +-include $(CSOURCES:%.c=%.d) + +# To remove generated files +clean: + rm -f $(COBJECTS) $(OBJECTS) $(SOURCES:%.cpp=%.d) $(CSOURCES:%.c=%.d) + +cleaner: clean + rm -f $(EXEC) diff --git a/tests/functor_tests/functor_tests.cpp b/tests/functor_tests/functor_tests.cpp new file mode 100755 index 0000000..9d235fb --- /dev/null +++ b/tests/functor_tests/functor_tests.cpp @@ -0,0 +1,156 @@ +#include +#include + +#include +#include + +using std::cout; +using std::endl; +using namespace Eigen; + +template +bool check_almost_equal(T a, T b, T tol) +{ + return std::abs(a - b) < tol; +} + +template +bool check_almost_equalV(const V & v1, const V & v2, T tol) +{ + bool ok = true; + for(size_t i = 0 ; i < v1.size() ; i++) + if(fabs(v1[i] - v2[i]) > tol) + { + ok = false; + break; + } + return ok; +} + +// --------------------------------------------------------------------------------------------------------------------------------------------------- + +template +T fct1(const T & x) +{ + return cos(x) + sqrt(x) + cbrt(sin(x)) + exp(-x*x)/log(x) + atanh(5*x) - 2*csc(x); +} + +template +T dfct1(const T & x) +{ + return -2*x*exp(-pow(x, 2))/log(x) - sin(x) + 2*cot(x)*csc(x) + (1.0/3.0)*cos(x)/pow(sin(x), 2.0/3.0) + 5/(-25*pow(x, 2) + 1) - exp(-pow(x, 2))/(x*pow(log(x), 2)) + (1.0/2.0)/sqrt(x); +} + +// --------------------------------------------------------------------------------------------------------------------------------------------------- + +/// Test function : length of vector +template +T fctNto1(const V & x) +{ + T res = T(0); + for(int i = 0 ; i < x.size() ; i++) + res += x[i]*x[i]; + return sqrt(res); +} + +/// Test function gradient : length of vector +template +V dfctNto1(const V & x) +{ + T length = fctNto1(x); + V res = x; + for(int i = 0 ; i < x.size() ; i++) + res[i] = x[i]/length; + return res; +} + +// --------------------------------------------------------------------------------------------------------------------------------------------------- + +constexpr int BIG_GLOBAL_N_STATIC = 10; +constexpr int BIG_GLOBAL_N_DYNAMIC = 100; + +CREATE_GRAD_FUNCTION_OBJECT_1_1(fct1, Grad_fct1); +CREATE_GRAD_FUNCTION_OBJECT_N_1_S(fctNto1, Grad_fctNto1S, BIG_GLOBAL_N_STATIC); +CREATE_GRAD_FUNCTION_OBJECT_N_1_D(fctNto1, Grad_fctNto1D, BIG_GLOBAL_N_DYNAMIC); + +int main() +{ + cout.precision(16); + using S = double; + constexpr S tol = 1e-12; + + { PRINT_STR("\nGradFunc1\n"); + S x = 0.1, fx, gradfx; + auto grad_fct1 = GradFunc1(fct1>, x); + grad_fct1.get_f_grad(x, fx, gradfx); + PRINT_VAR(x); + PRINT_VAR(fx); + PRINT_VAR(gradfx); + PRINT_VAR(dfct1(x)); + assert(check_almost_equal(gradfx, dfct1(x), tol)); + } + + { PRINT_STR("\nCREATE_GRAD_FUNCTION_OBJECT_1_1\n"); + S x = 0.1, fx, gradfx; + auto grad_fct1 = Grad_fct1(); + grad_fct1.get_f_grad(x, fx, gradfx); + PRINT_VAR(x); + PRINT_VAR(fx); + PRINT_VAR(gradfx); + PRINT_VAR(dfct1(x)); + assert(check_almost_equal(gradfx, dfct1(x), tol)); + } + + { PRINT_STR("\nCREATE_GRAD_FUNCTION_OBJECT_N_1_S\n"); + constexpr int N = BIG_GLOBAL_N_STATIC; + // using D = DualD; + using V = Eigen::Array; + V x; + for(int i = 0 ; i < x.size() ; i++) + x[i] = S(i); + S fx = fctNto1(x); + V dfx = dfctNto1(x); + PRINT_VAR(fx); + PRINT_VEC(dfx); + + Grad_fctNto1S grad_fctNto1; + S fx2; V gradfx2; + grad_fctNto1.get_f_grad(x, fx2, gradfx2); + PRINT_VAR(fx2); + PRINT_VEC(gradfx2); + PRINT_VEC(grad_fctNto1(x)); + + for (size_t i = 0; i < N; i++) + cout << dfx[i] << "\t" << gradfx2[i] << "\t" << (dfx[i] - gradfx2[i]) << endl; + + assert(check_almost_equal(fx, fx2, tol)); + assert(check_almost_equalV(dfx, gradfx2, tol)); + } + + //* + { PRINT_STR("\nCREATE_GRAD_FUNCTION_OBJECT_N_1_D\n"); + constexpr int N = BIG_GLOBAL_N_DYNAMIC; + // using D = DualD; + using V = Eigen::Array; + V x(N); + for(int i = 0 ; i < x.size() ; i++) + x[i] = S(i); + S fx = fctNto1(x); + V dfx = dfctNto1(x); + PRINT_VAR(fx); + PRINT_VEC(dfx); + + Grad_fctNto1D grad_fctNto1; + S fx2; V gradfx2; + grad_fctNto1.get_f_grad(x, fx2, gradfx2); + PRINT_VAR(fx2); + PRINT_VEC(gradfx2); + PRINT_VEC(grad_fctNto1(x)); + + for (size_t i = 0; i < N; i++) + cout << dfx[i] << "\t" << gradfx2[i] << "\t" << (dfx[i] - gradfx2[i]) << endl; + + assert(check_almost_equal(fx, fx2, tol)); + assert(check_almost_equalV(dfx, gradfx2, tol)); + }//*/ +} diff --git a/utils.hpp b/utils.hpp index 3f5243e..384e9cb 100755 --- a/utils.hpp +++ b/utils.hpp @@ -22,8 +22,8 @@ std::ostream& operator<<(std::ostream& os, const C& objs) return os; }//*/ -/// Convenience template class to do StdPairValueCatcher(a, b) = someFunctionThatReturnsAnStdPair() -/// Instead of doing a = std::get<0>(result_from_function), b = std::get<1>(result_from_function) +/// Convenience template class to do `StdPairValueCatcher(a, b) = someFunctionThatReturnsAnStdPair()` +/// Instead of doing `a = std::get<0>(result_from_function), b = std::get<1>(result_from_function)` template struct StdPairValueCatcher {