From d9124bfb4bdfedede760b1d6ad335667e6c87bc7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A9r=C3=B4me?= Date: Sun, 31 Mar 2019 15:46:52 +0200 Subject: [PATCH] MAJOR MILESTONE : rewrote everything so that 2 template classes are created : static and dynamic. --- .gitignore | 1 + AutomaticDifferentiation.hpp | 557 +----------------- AutomaticDifferentiation_base.hpp | 443 ++++++++++++++ AutomaticDifferentiation_dynamic.hpp | 13 + AutomaticDifferentiation_math.hpp | 77 +++ AutomaticDifferentiation_static.hpp | 57 ++ Makefile | 3 + test/Makefile | 10 +- test/dual_test.cpp | 524 +++++++++++++++- test/py_dual_test_generate_function_tests.py | 65 ++ tests/minimalistic_test_autodiff/Makefile | 37 ++ .../minimalistic_test_autodiff.cpp | 257 ++++++++ utils.hpp | 3 +- 13 files changed, 1485 insertions(+), 562 deletions(-) create mode 100755 AutomaticDifferentiation_base.hpp create mode 100755 AutomaticDifferentiation_dynamic.hpp create mode 100755 AutomaticDifferentiation_math.hpp create mode 100755 AutomaticDifferentiation_static.hpp create mode 100755 test/py_dual_test_generate_function_tests.py create mode 100755 tests/minimalistic_test_autodiff/Makefile create mode 100755 tests/minimalistic_test_autodiff/minimalistic_test_autodiff.cpp diff --git a/.gitignore b/.gitignore index 4f79056..8ea2d5c 100755 --- a/.gitignore +++ b/.gitignore @@ -7,6 +7,7 @@ html/ # catch2 generated files *.gcda *.gcno +*.gcov # Prerequisites *.d diff --git a/AutomaticDifferentiation.hpp b/AutomaticDifferentiation.hpp index c104968..fa7a931 100755 --- a/AutomaticDifferentiation.hpp +++ b/AutomaticDifferentiation.hpp @@ -1,550 +1,7 @@ -#ifndef DEF_AUTOMATIC_DIFFERENTIATION -#define DEF_AUTOMATIC_DIFFERENTIATION - -#include -#include -#include - -#include - -// template -// std::ostream & operator<<(std::ostream & out, std::valarray const& v) -// { -// for(size_t i = 0 ; i < v.size() ; i++) -// out << v[i] << " "; -// return out; -// } - -#define MINAB(a, b) (((a) < (b)) ? (a) : (b)) - -/// 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::d(0), y+Dual::d(1), z+Dual::d(2), ...) -/// reference : http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.89.7749&rep=rep1&type=pdf -template -class Dual -{ - public: - using VectorT = std::valarray; - - static VectorT __create_VectorT_zeros(int N = 1) - { - assert(N >= 0); - VectorT res(Scalar(0.), N); - return res; - } - - Dual(const Scalar & _a = Scalar(0.), const VectorT & _b = Dual::__create_VectorT_zeros()) - : a(_a), - b(_b) - {} - - Dual const& operator=(Scalar const& _a) - { - *this = Dual(_a); - } - - /// Use this function to set what variable is to be derived : x + Dual::d(i) - static Dual D(int i = 0, int N = 1) - { - assert(i >= 0); - assert(i < N); - VectorT _d = Dual::__create_VectorT_zeros(N); - _d[i] = Scalar(1.); - return Dual(Scalar(0.), _d); - } - - /// Use this function to set what variable is to be derived. - Dual const& diff(int i = 0, int N = 1) - { - assert(i >= 0); - assert(i < N); - if(N != b.size()) - { - // copy old data into new b vector - VectorT b_old = b; - b.resize(N); - for(size_t j = 0 ; j < min(b.size(), b_old.size()) ; j++) - b[j] = b_old[j]; - } - b[i] = Scalar(1.); - return *this; - } - - /// Returns the value - Scalar const& x() const - { - return a; - } - - /// Returns the value - Scalar & x() - { - return a; - } - - /// Returns the derivative value at index i - Scalar const& d(int i) const - { - assert(i >= 0); - assert(i < b.size()); - return b[i]; - } - - /// Returns the derivative value at index i - Scalar & d(int i) - { - assert(i >= 0); - assert(i < b.size()); - return b[i]; - } - - 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 - VectorT b; /// Infinitesimal parts -}; - -//* -template -Dual operator+(A const& v, Dual const& x) { - return (Dual(static_cast(v)) + x); -} -template -Dual operator-(A const& v, Dual const& x) { - return (Dual(static_cast(v)) - x); -} -template -Dual operator*(A const& v, Dual const& x) { - return (Dual(static_cast(v)) * x); -} -template -Dual operator/(A const& v, Dual const& x) { - return (Dual(static_cast(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 cbrt(const Dual & x) { - return Dual(cbrt(x.a), x.b/(Scalar(3.)*pow(x.a, Scalar(2.)/Scalar(3.)))); -} - -template Dual sign(const Dual & x) { - return Dual(sign(x.a), Dual::Dual::__create_VectorT_zeros()); -} - -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), Dual::Dual::__create_VectorT_zeros()); -} - -template Dual floor(const Dual & x) { - return Dual(floor(x.a), Dual::Dual::__create_VectorT_zeros()); -} - -template Dual ceil(const Dual & x) { - return Dual(ceil(x.a), Dual::Dual::__create_VectorT_zeros()); -} - -template Dual round(const Dual & x) { - return Dual(round(x.a), Dual::Dual::__create_VectorT_zeros()); -} - -template std::ostream & operator<<(std::ostream & s, const Dual & x) -{ - return (s << x.a); -} - -/// Function object to evaluate the derivative of a function anywhere without explicitely using the Dual numbers. -//* -template -struct GradFunc -{ - Func f; - - GradFunc(Func f_, Scalar) : f(f_) {} - - // 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. - Dual X(x); - X.diff(0,1); - Dual Y = f(X); - return Y.d(0); - } - - /// 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. - void get_f_grad(Scalar const& x, Scalar & fx, Scalar & gradfx) - { - // differentiate using the Dual number and return the .b component. - Dual X(x); - X.diff(0,1); - Dual Y = f(X); - fx = Y.x(); - gradfx = Y.d(0); - } -};//*/ - -//* -/// Macro to create a function object that returns the gradient of the function at X. -/// Designed to work with functions, lambdas, etc. -#define CREATE_GRAD_FUNCTION_OBJECT(Func, GradFuncName) \ -struct GradFuncName { \ - template \ - Scalar operator()(Scalar const& x) { \ - Dual X(x); \ - X.diff(0,1); \ - Dual Y = Func>(X); \ - return Y.d(0); \ - } \ - template \ - void get_f_grad(Scalar const& x, Scalar & fx, Scalar & gradfx) { \ - Dual X(x); \ - X.diff(0,1); \ - Dual 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. -/// Designed to work with function objects. -#define CREATE_GRAD_FUNCTION_OBJECT_FUNCTOR(Func, GradFuncName) \ -struct GradFuncName { \ - template \ - Scalar operator()(Scalar const& x) { \ - Dual X(x); \ - X.diff(0,1); \ - Func f; \ - Dual Y = f(X); \ - return Y.d(0); \ - } \ - template \ - void get_f_grad(Scalar const& x, Scalar & fx, Scalar & gradfx) { \ - Dual X(x); \ - X.diff(0,1); \ - Func f; \ - Dual Y = f(X); \ - fx = Y.x(); \ - gradfx = Y.d(0); \ - } \ -} -//*/ - -#endif +#ifndef DEF_AUTOMATIC_DIFFERENTIATION +#define DEF_AUTOMATIC_DIFFERENTIATION + +#include "AutomaticDifferentiation_static.hpp" +#include "AutomaticDifferentiation_dynamic.hpp" + +#endif diff --git a/AutomaticDifferentiation_base.hpp b/AutomaticDifferentiation_base.hpp new file mode 100755 index 0000000..6b801a6 --- /dev/null +++ b/AutomaticDifferentiation_base.hpp @@ -0,0 +1,443 @@ +#include +#include +#include + +#include +#include "AutomaticDifferentiation_math.hpp" + +// Convenience defines to make the code easier to read. + +/// Type of the DualBase object. +#define __DualBase_t __Dual_DualBase + +/// Template declaration of DualBase object. +#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. +/// +/// reference : http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.89.7749&rep=rep1&type=pdf +template +struct __Dual_DualBase +{ + static_assert(N > 0, "N must be > 0."); + + using VectorT = __Dual_VectorT; + + /// Sets all the elements of the b vector to 0. + void SetBToZero() + { + #if __Dual_bdynamic == 1 + b = VectorT::Zero(N); + #else + b = VectorT::Zero(); + #endif + } + + __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& operator=(Scalar const& _a) + // { + // return *this; + // } + + /// 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 : + /// __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 vector of infinitesimal parts + VectorT const& B() const { return b; } + VectorT & B() { return b; } + + /// 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 + 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; + // 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) + { + 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++() { 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+(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 + + 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 +}; + +template +__Dual_DualBase operator+(A const& v, __Dual_DualBase const& x) { + return (__Dual_DualBase(v) + x); +} +template +__Dual_DualBase operator-(A const& v, __Dual_DualBase const& x) { + return (__Dual_DualBase(v) - x); +} +template +__Dual_DualBase operator*(A const& v, __Dual_DualBase const& x) { + return (__Dual_DualBase(v) * x); +} +template +__Dual_DualBase operator/(A const& v, __Dual_DualBase const& x) { + return (__Dual_DualBase(v) / x); +} + +// Basic mathematical functions for __Dual_DualBase numbers +// f(a + b*d) = f(a) + b*f'(a)*d + +// Trigonometric functions +__tplDualBase_t __DualBase_t cos(const __DualBase_t & x) { + return __DualBase_t(cos(x.a), -x.b*sin(x.a)); +} + +__tplDualBase_t __DualBase_t sin(const __DualBase_t & x) { + return __DualBase_t(sin(x.a), x.b*cos(x.a)); +} + +__tplDualBase_t __DualBase_t tan(const __DualBase_t & x) { + return __DualBase_t(tan(x.a), x.b*sec(x.a)*sec(x.a)); +} + +__tplDualBase_t __DualBase_t sec(const __DualBase_t & x) { + return __DualBase_t(sec(x.a), x.b*sec(x.a)*tan(x.a)); +} + +__tplDualBase_t __DualBase_t cot(const __DualBase_t & x) { + return __DualBase_t(cot(x.a), x.b*(-csc(x.a)*csc(x.a))); +} + +__tplDualBase_t __DualBase_t csc(const __DualBase_t & x) { + return __DualBase_t(csc(x.a), x.b*(-cot(x.a)*csc(x.a))); +} + +// Inverse trigonometric functions +__tplDualBase_t __DualBase_t acos(const __DualBase_t & x) { + return __DualBase_t(acos(x.a), x.b*(-Scalar(1.)/sqrt(Scalar(1.)-x.a*x.a))); +} + +__tplDualBase_t __DualBase_t asin(const __DualBase_t & x) { + return __DualBase_t(asin(x.a), x.b*(Scalar(1.)/sqrt(Scalar(1.)-x.a*x.a))); +} + +__tplDualBase_t __DualBase_t atan(const __DualBase_t & x) { + return __DualBase_t(atan(x.a), x.b*(Scalar(1.)/(x.a*x.a + Scalar(1.)))); +} + +__tplDualBase_t __DualBase_t asec(const __DualBase_t & x) { + return __DualBase_t(asec(x.a), x.b*(Scalar(1.)/(sqrt(Scalar(1.)-Scalar(1.)/(x.a*x.a))*(x.a*x.a)))); +} + +__tplDualBase_t __DualBase_t acot(const __DualBase_t & x) { + return __DualBase_t(acot(x.a), x.b*(-Scalar(1.)/((x.a*x.a)+Scalar(1.)))); +} + +__tplDualBase_t __DualBase_t acsc(const __DualBase_t & x) { + return __DualBase_t(acsc(x.a), x.b*(-Scalar(1.)/(sqrt(Scalar(1.)-Scalar(1.)/(x.a*x.a))*(x.a*x.a)))); +} + +// Hyperbolic trigonometric functions +__tplDualBase_t __DualBase_t cosh(const __DualBase_t & x) { + return __DualBase_t(cosh(x.a), x.b*sinh(x.a)); +} + +__tplDualBase_t __DualBase_t sinh(const __DualBase_t & x) { + return __DualBase_t(sinh(x.a), x.b*cosh(x.a)); +} + +__tplDualBase_t __DualBase_t tanh(const __DualBase_t & x) { + return __DualBase_t(tanh(x.a), x.b*sech(x.a)*sech(x.a)); +} + +__tplDualBase_t __DualBase_t sech(const __DualBase_t & x) { + return __DualBase_t(sech(x.a), x.b*(-sech(x.a)*tanh(x.a))); +} + +__tplDualBase_t __DualBase_t coth(const __DualBase_t & x) { + return __DualBase_t(coth(x.a), x.b*(-csch(x.a)*csch(x.a))); +} + +__tplDualBase_t __DualBase_t csch(const __DualBase_t & x) { + return __DualBase_t(csch(x.a), x.b*(-coth(x.a)*csch(x.a))); +} + +// Inverse hyperbolic trigonometric functions +__tplDualBase_t __DualBase_t acosh(const __DualBase_t & x) { + return __DualBase_t(acosh(x.a), x.b*(Scalar(1.)/sqrt((x.a*x.a)-Scalar(1.)))); +} + +__tplDualBase_t __DualBase_t asinh(const __DualBase_t & x) { + return __DualBase_t(asinh(x.a), x.b*(Scalar(1.)/sqrt((x.a*x.a)+Scalar(1.)))); +} + +__tplDualBase_t __DualBase_t atanh(const __DualBase_t & x) { + return __DualBase_t(atanh(x.a), x.b*(Scalar(1.)/(Scalar(1.)-(x.a*x.a)))); +} + +__tplDualBase_t __DualBase_t asech(const __DualBase_t & x) { + return __DualBase_t(asech(x.a), x.b*(Scalar(-1.)/(sqrt(Scalar(1.)/(x.a*x.a)-Scalar(1.))*(x.a*x.a)))); +} + +__tplDualBase_t __DualBase_t acoth(const __DualBase_t & x) { + return __DualBase_t(acoth(x.a), x.b*(-Scalar(1.)/((x.a*x.a)-Scalar(1.)))); +} + +__tplDualBase_t __DualBase_t acsch(const __DualBase_t & x) { + return __DualBase_t(acsch(x.a), x.b*(-Scalar(1.)/(sqrt(Scalar(1.)/(x.a*x.a)+Scalar(1.))*(x.a*x.a)))); +} + +// Exponential functions +__tplDualBase_t __DualBase_t exp(const __DualBase_t & x) { + return __DualBase_t(exp(x.a), x.b*exp(x.a)); +} + +__tplDualBase_t __DualBase_t log(const __DualBase_t & x) { + return __DualBase_t(log(x.a), x.b/x.a); +} + +__tplDualBase_t __DualBase_t exp10(const __DualBase_t & x) { + return __DualBase_t(exp10(x.a), x.b*(log(Scalar(10.))*exp10(x.a))); +} + +__tplDualBase_t __DualBase_t log10(const __DualBase_t & x) { + return __DualBase_t(log10(x.a), x.b/(log(Scalar(10.))*x.a)); +} + +__tplDualBase_t __DualBase_t exp2(const __DualBase_t & x) { + return __DualBase_t(exp2(x.a), x.b*(log(Scalar(2.))*exp2(x.a))); +} + +__tplDualBase_t __DualBase_t log2(const __DualBase_t & x) { + return __DualBase_t(log2(x.a), x.b/(log(Scalar(2.))*x.a)); +} + +__tplDualBase_t __DualBase_t pow(const __DualBase_t & x, const __DualBase_t & n) { + return exp(n*log(x)); +} + +template __DualBase_t pow(const __DualBase_t & x, const Scalar2 & n) { + return exp(__DualBase_t(static_cast(n))*log(x)); +} + +template __DualBase_t pow(const Scalar2 & x, const __DualBase_t & n) { + return exp(__DualBase_t(n)*log(static_cast(x))); +} + +// Other functions +__tplDualBase_t __DualBase_t sqrt(const __DualBase_t & x) { + return __DualBase_t(sqrt(x.a), x.b/(Scalar(2.)*sqrt(x.a))); +} + +__tplDualBase_t __DualBase_t cbrt(const __DualBase_t & x) { + return __DualBase_t(cbrt(x.a), x.b/(Scalar(3.)*pow(x.a, Scalar(2.)/Scalar(3.)))); +} + +__tplDualBase_t __DualBase_t sign(const __DualBase_t & x) { + return __DualBase_t(sign(x.a)); +} + +__tplDualBase_t __DualBase_t abs(const __DualBase_t & x) { + return __DualBase_t(abs(x.a), x.b*sign(x.a)); +} +__tplDualBase_t __DualBase_t fabs(const __DualBase_t & x) { + return __DualBase_t(fabs(x.a), x.b*sign(x.a)); +} + +__tplDualBase_t __DualBase_t heaviside(const __DualBase_t & x) { + return __DualBase_t(heaviside(x.a)); +} + +__tplDualBase_t __DualBase_t floor(const __DualBase_t & x) { + return __DualBase_t(floor(x.a)); +} + +__tplDualBase_t __DualBase_t ceil(const __DualBase_t & x) { + return __DualBase_t(ceil(x.a)); +} + +__tplDualBase_t __DualBase_t round(const __DualBase_t & x) { + return __DualBase_t(round(x.a)); +} + +__tplDualBase_t std::ostream & operator<<(std::ostream & s, const __Dual_DualBase & x) +{ + return (s << x.a); +} + +// ----------------------------------------------------------------------------------------------------------- + +/* +/// Macro to create a function object that returns the gradient of the function at X. +/// Designed to work with functions, lambdas, etc. +#define CREATE_GRAD_FUNCTION_OBJECT(Func, GradFuncName) \ +struct GradFuncName { \ + template \ + Scalar operator()(Scalar const& x) { \ + __Dual_DualBase X(x); \ + X.diff(0,1); \ + __Dual_DualBase Y = Func<__Dual_DualBase>(X); \ + return Y.d(0); \ + } \ + template \ + void get_f_grad(Scalar const& x, Scalar & fx, Scalar & gradfx) { \ + __Dual_DualBase X(x); \ + X.diff(0,1); \ + __Dual_DualBase Y = Func<__Dual_DualBase>(X); \ + fx = Y.x(); \ + gradfx = Y.d(0); \ + } \ +} +//*/ +/* +/// Macro to create a function object that returns the gradient of the function at X. +/// Designed to work with function objects. +#define CREATE_GRAD_FUNCTION_OBJECT_FUNCTOR(Func, GradFuncName) \ +struct GradFuncName { \ + template \ + Scalar operator()(Scalar const& x) { \ + __Dual_DualBase X(x); \ + X.diff(0,1); \ + Func f; \ + __Dual_DualBase Y = f(X); \ + return Y.d(0); \ + } \ + template \ + void get_f_grad(Scalar const& x, Scalar & fx, Scalar & gradfx) { \ + __Dual_DualBase X(x); \ + X.diff(0,1); \ + Func f; \ + __Dual_DualBase Y = f(X); \ + fx = Y.x(); \ + gradfx = Y.d(0); \ + } \ +} +//*/ diff --git a/AutomaticDifferentiation_dynamic.hpp b/AutomaticDifferentiation_dynamic.hpp new file mode 100755 index 0000000..8668781 --- /dev/null +++ b/AutomaticDifferentiation_dynamic.hpp @@ -0,0 +1,13 @@ +#ifndef DEF_AUTOMATIC_DIFFERENTIATION_DYNAMIC +#define DEF_AUTOMATIC_DIFFERENTIATION_DYNAMIC + +#undef __Dual_DualBase +#undef __Dual_VectorT +#undef __Dual_bdynamic +#define __Dual_DualBase DualD +#define __Dual_VectorT Eigen::Array +#define __Dual_bdynamic 1 + +#include "AutomaticDifferentiation_base.hpp" + +#endif diff --git a/AutomaticDifferentiation_math.hpp b/AutomaticDifferentiation_math.hpp new file mode 100755 index 0000000..0f950bc --- /dev/null +++ b/AutomaticDifferentiation_math.hpp @@ -0,0 +1,77 @@ +#ifndef DEF_AUTOMATIC_DIFFERENTIATION_MATH +#define DEF_AUTOMATIC_DIFFERENTIATION_MATH + +#include + +// 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; +} + +#endif diff --git a/AutomaticDifferentiation_static.hpp b/AutomaticDifferentiation_static.hpp new file mode 100755 index 0000000..b9e0987 --- /dev/null +++ b/AutomaticDifferentiation_static.hpp @@ -0,0 +1,57 @@ +#ifndef DEF_AUTOMATIC_DIFFERENTIATION_STATIC +#define DEF_AUTOMATIC_DIFFERENTIATION_STATIC + +#undef __Dual_DualBase +#undef __Dual_VectorT +#undef __Dual_bdynamic +#define __Dual_DualBase DualS +#define __Dual_VectorT Eigen::Array +#define __Dual_bdynamic 0 + +#include "AutomaticDifferentiation_base.hpp" + +/// Function object to evaluate the derivative of a univariate function anywhere without explicitely using the dual numbers. +/// +/// Here is an example of use : +/// +/// template fct(T x) { return exp(-x*x); } +/// +/// Scalar x = 3.14; +/// auto gradFct = GradFunc1(fct, x); // x is only passed to the function so that the full type of GradFunc1 does not need to be written manually. +/// Scalar dfdx = gradFct(x); // Evaluation of the gradient of fct at x. +/// +/// // Alternatively : +/// Scalar fx, dfdx; +/// gradFct.get_f_grad(x, fx, dfdx); // Evaluation of the gradient of fct at x. +template +struct GradFunc1 +{ + Func f; + + /// Constructor of the gradient function object. The second parameter IS NOT USED FOR THE COMPUTATION OF THE GRADIENT, but only as a convenience for the + /// automatic detection of the Scalar variable type. + GradFunc1(Func f_, Scalar) : f(f_) {} + + // 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); + } + + /// 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. + void get_f_grad(Scalar const& x, Scalar & fx, Scalar & gradfx) + { + // differentiate using the dual number and return the .b component. + DualS X(x); + X.diff(0); + DualS Y = f(X); + fx = Y.x(); + gradfx = Y.d(0); + } +}; + +#endif diff --git a/Makefile b/Makefile index 01e733e..1a720a2 100755 --- a/Makefile +++ b/Makefile @@ -29,6 +29,9 @@ test_AutomaticDifferentiation_vector: $(OBJ)/test_AutomaticDifferentiation_vecto $(OBJ)/test_AutomaticDifferentiation_vector.o: test_AutomaticDifferentiation_vector.cpp $(CXX) $(CFLAGS) $(LDFLAGS) -c test_AutomaticDifferentiation_vector.cpp +doc: + doxygen Doxyfile + clean: rm *.o diff --git a/test/Makefile b/test/Makefile index a1efaa2..5fa1871 100755 --- a/test/Makefile +++ b/test/Makefile @@ -1,11 +1,11 @@ # Declaration of variables C = clang -COMMON_FLAGS = -Wall -MMD +COMMON_FLAGS = -Wall -MMD -fprofile-arcs -ftest-coverage C_FLAGS = $(COMMON_FLAGS) -CC = clang++ +CC = g++ CC_FLAGS = $(COMMON_FLAGS) -std=c++17 -g -O0#-O2 -LD_FLAGS = -lgmp -lmpfr -INCLUDES = -I../../Catch2-master/single_include +LD_FLAGS = -lgmp -lmpfr -lgcov +INCLUDES = -I../../Catch2-master/single_include -I/usr/include/eigen3 # File names EXEC = run @@ -31,7 +31,7 @@ $(EXEC): $(COBJECTS) $(OBJECTS) # To remove generated files clean: - rm -f $(COBJECTS) $(OBJECTS) $(SOURCES:%.cpp=%.d) $(CSOURCES:%.c=%.d) + rm -f $(COBJECTS) $(OBJECTS) $(SOURCES:%.cpp=%.d) $(CSOURCES:%.c=%.d) *.gcov *.gcda *.gcno cleaner: clean rm -f $(EXEC) diff --git a/test/dual_test.cpp b/test/dual_test.cpp index f73803a..1fe8668 100755 --- a/test/dual_test.cpp +++ b/test/dual_test.cpp @@ -20,20 +20,532 @@ bool check_almost_equal(const std::vector & v1, const std::vector & v2, T return ok; } +using S = double; +constexpr S tol = 1e-12; + +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); +} + +template +struct Fct1 +{ + T operator()(const T & x) { return fct1(x); } +}; + +template +T fct2(const T & x, const T & y, const T & z) +{ + T two(2); + return 2*sqrt(x)*sqrt(y)/log(sin(z)) + sqrt(y) + exp(-pow(x, two) - pow(y, two) - pow(z, two))/(log(x)*log(y)*log(z)) + cbrt(sin(z))*tan(y) + cos(x) + 2*atanh(5*x)*csc(y); +} +template +T dfct2dx(const T & x, const T & y, const T & z) +{ + T two(2); + return -2*x*exp(-pow(x, two) - pow(y, two) - pow(z, two))/(log(x)*log(y)*log(z)) - sin(x) + 10*csc(y)/(-25*pow(x, two) + 1) - exp(-pow(x, two) - pow(y, two) - pow(z, two))/(x*pow(log(x), two)*log(y)*log(z)) + sqrt(y)/(sqrt(x)*log(sin(z))); +} +template +T dfct2dy(const T & x, const T & y, const T & z) +{ + T two(2); + return sqrt(x)/(sqrt(y)*log(sin(z))) - 2*y*exp(-pow(x, two) - pow(y, two) - pow(z, two))/(log(x)*log(y)*log(z)) + (pow(tan(y), two) + 1)*cbrt(sin(z)) - 2*cot(y)*atanh(5*x)*csc(y) - exp(-pow(x, two) - pow(y, two) - pow(z, two))/(y*log(x)*pow(log(y), two)*log(z)) + (1.0/2.0)/sqrt(y); +} +template +T dfct2dz(const T & x, const T & y, const T & z) +{ + T two(2); + return -2*sqrt(x)*sqrt(y)*cos(z)/(pow(log(sin(z)), two)*sin(z)) - 2*z*exp(-pow(x, two) - pow(y, two) - pow(z, two))/(log(x)*log(y)*log(z)) + (1.0/3.0)*cos(z)*tan(y)/pow(sin(z), 2.0/3.0) - exp(-pow(x, two) - pow(y, two) - pow(z, two))/(z*log(x)*log(y)*pow(log(z), two)); +} + +/// Spherical coordinates [r, lat, lon] to cartesian [x, y, z] transformation. +template +Eigen::Array sph2cart(const Eigen::Array & xsph) +{ + Eigen::Array xcart; + xcart[0] = xsph[0]*cos(xsph[1])*cos(xsph[2]); + xcart[1] = xsph[0]*cos(xsph[1])*sin(xsph[2]); + xcart[2] = xsph[0]*sin(xsph[1]); + return xcart; +} + +/// Gradient of spherical coordinates [r, lat, lon] to cartesian [x, y, z] transformation. +template +Eigen::Matrix sph2cartJ(const Eigen::Array & xsph) +{ + // Jij = Dfi/dxj + Eigen::Matrix Jcart; + Jcart << cos(xsph[2])*cos(xsph[1]), -xsph[0]*sin(xsph[1])*cos(xsph[2]), -xsph[0]*sin(xsph[2])*cos(xsph[1]), + sin(xsph[2])*cos(xsph[1]), -xsph[0]*sin(xsph[2])*sin(xsph[1]), xsph[0]*cos(xsph[2])*cos(xsph[1]), + sin(xsph[1]), xsph[0]*cos(xsph[1]), T(0); + return Jcart; +} + +// dxdr = cos(xsph[2])*cos(xsph[1]), +// dxth = -xsph[0]*sin(xsph[1])*cos(xsph[2]), +// dxdphi = -xsph[0]*sin(xsph[2])*cos(xsph[1]) +// +// dydr = sin(xsph[2])*cos(xsph[1]), +// dydth = -xsph[0]*sin(xsph[2])*sin(xsph[1]), +// dydxsph[2] = xsph[0]*cos(xsph[2])*cos(xsph[1]) +// +// dzdr = sin(xsph[1]), +// dzth = xsph[0]*cos(xsph[1]), +// dzphi = 0 + TEST_CASE( "Dual numbers : Constructors", "[dual]" ) { std::cout.precision(16); - using S = double; - using D = Dual; - using V = std::vector; + constexpr int N = 2; + using D = DualS; + using V = D::VectorT; SECTION( "Default Constructor" ) { - + D a; + REQUIRE(a.a == 0); + REQUIRE(a.b.size() == N); + for(int i = 0 ; i < N ; i++) + REQUIRE(a.b[i] == 0); } - SECTION( "Constructor with int" ) { - + SECTION( "Constructor with a" ) { + S x = 3.14159; + D a(x); + REQUIRE(a.a == x); + REQUIRE(a.b.size() == N); + for(int i = 0 ; i < N ; i++) + REQUIRE(a.b[i] == 0); } + + SECTION( "Constructor with a and b" ) { + S x = 3.14159; + V y = {1.5, 3.2}; + D a(x, y); + REQUIRE(a.a == x); + REQUIRE(a.b.size() == 2); + REQUIRE(a.b[0] == y[0]); + REQUIRE(a.b[1] == y[1]); + } + + SECTION( "operator= Dual" ) { + S x = 3.14159; + V y = {1.5, 3.2}; + D a, b(x, y); + a = b; + REQUIRE(a.a == x); + REQUIRE(a.b.size() == 2); + REQUIRE(a.b[0] == y[0]); + REQUIRE(a.b[1] == y[1]); + } + + SECTION( "operator= Dual" ) { + S x = 3.14159; + D a; + a = x; + REQUIRE(a.a == x); + REQUIRE(a.b.size() == N); + for(int i = 0 ; i < N ; i++) + REQUIRE(a.b[i] == 0); + } +} + +TEST_CASE( "Dual numbers : access parts", "[dual]" ) +{ + std::cout.precision(16); + constexpr int N = 2; + using D = DualS; + using V = D::VectorT; + + S x = 3.14159; + V v = {1.5, 3.2}; + + SECTION( "x()" ) { + D a(x, v); + REQUIRE(a.x() == x); + } + + SECTION( "d()" ) { + D a(x, v); + REQUIRE(a.d(0) == v[0]); + REQUIRE(a.d(1) == v[1]); + } +} + +TEST_CASE( "Dual numbers : D() and diff()", "[dual]" ) +{ + std::cout.precision(16); + constexpr int N = 3; + using D = DualS; + + S x = 3.14159; + + SECTION( "D(0)" ) { + D a = D::D(0); + REQUIRE(a.a == 0); + REQUIRE(a.b.size() == N); + REQUIRE(a.d(0) == 1); + } + + SECTION( "D(0, 3)" ) { + D a = D::D(0); + REQUIRE(a.a == 0); + REQUIRE(a.b.size() == N); + REQUIRE(a.d(0) == 1); + REQUIRE(a.d(1) == 0); + REQUIRE(a.d(2) == 0); + } + + SECTION( "D(1, 3)" ) { + D a = D::D(1); + REQUIRE(a.a == 0); + REQUIRE(a.b.size() == N); + REQUIRE(a.d(0) == 0); + REQUIRE(a.d(1) == 1); + REQUIRE(a.d(2) == 0); + } + + SECTION( "D(2, 3)" ) { + D a = D::D(2); + REQUIRE(a.b.size() == N); + REQUIRE(a.d(0) == 0); + REQUIRE(a.d(1) == 0); + REQUIRE(a.d(2) == 1); + } + + SECTION( "diff(0)" ) { + D a(x); + a.diff(0); + REQUIRE(a.b.size() == N); + REQUIRE(a.d(0) == 1); + REQUIRE(a.d(1) == 0); + REQUIRE(a.d(2) == 0); + } + + SECTION( "diff(i,N)" ) { + D a(x); + for(int i = 0 ; i < N ; i++) + { + a.diff(i); + REQUIRE(a.a == x); + REQUIRE(a.b.size() == N); + for(int j = 0 ; j < N ; j++) + REQUIRE(a.d(j) == (j==i)); + } + } +} + +#define MAKE_TEST_SECTION_COMPARISON_OPERATOR(Op, A, B)\ +SECTION( "operator"#Op ) { \ + D a(A), b(B); \ + bool resD = a Op b; \ + bool resS = A Op B; \ + REQUIRE(resD == resS); \ +} + +TEST_CASE( "Dual numbers : comparison operators", "[dual]" ) +{ + std::cout.precision(16); + constexpr int N = 3; + using D = DualS; + + MAKE_TEST_SECTION_COMPARISON_OPERATOR(==, 1., 1.) + MAKE_TEST_SECTION_COMPARISON_OPERATOR(==, 1., 2.) + + MAKE_TEST_SECTION_COMPARISON_OPERATOR(!=, 1., 1.) + MAKE_TEST_SECTION_COMPARISON_OPERATOR(!=, 1., 2.) + + MAKE_TEST_SECTION_COMPARISON_OPERATOR(<, 1., 1.) + MAKE_TEST_SECTION_COMPARISON_OPERATOR(<, 2., 1.) + MAKE_TEST_SECTION_COMPARISON_OPERATOR(<, 1., 2.) + + MAKE_TEST_SECTION_COMPARISON_OPERATOR(<=, 1., 1.) + MAKE_TEST_SECTION_COMPARISON_OPERATOR(<=, 2., 1.) + MAKE_TEST_SECTION_COMPARISON_OPERATOR(<=, 1., 2.) + + MAKE_TEST_SECTION_COMPARISON_OPERATOR(>, 1., 1.) + MAKE_TEST_SECTION_COMPARISON_OPERATOR(>, 2., 1.) + MAKE_TEST_SECTION_COMPARISON_OPERATOR(>, 1., 2.) + + MAKE_TEST_SECTION_COMPARISON_OPERATOR(>=, 1., 1.) + MAKE_TEST_SECTION_COMPARISON_OPERATOR(>=, 2., 1.) + MAKE_TEST_SECTION_COMPARISON_OPERATOR(>=, 1., 2.) +} + +#define MAKE_TEST_SECTION_BASIC_BINARY_OPERATOR(Op, A, B)\ +SECTION( "operator"#Op ) { \ + D a(A), b(B); \ + S resD = (a Op b).a; \ + S resS = A Op B; \ + REQUIRE(resD == resS); \ +} + +TEST_CASE( "Dual numbers : basic binary operators", "[dual]" ) +{ + std::cout.precision(16); + constexpr int N = 3; + using D = DualS; + + MAKE_TEST_SECTION_BASIC_BINARY_OPERATOR(+, 3.14, 2.71) + MAKE_TEST_SECTION_BASIC_BINARY_OPERATOR(-, 3.14, 2.71) + MAKE_TEST_SECTION_BASIC_BINARY_OPERATOR(*, 3.14, 2.71) + MAKE_TEST_SECTION_BASIC_BINARY_OPERATOR(/, 3.14, 2.71) +} + +#define MAKE_TEST_SECTION_BASIC_BINARY_OPERATOR_TYPE_MISMATCH_POST(Op, A, B)\ +SECTION( "operator"#Op ) { \ + D a = D(A) Op (B); \ + S resS = A Op B; \ + REQUIRE(a.a == resS); \ +} + +TEST_CASE( "Dual numbers : basic binary operators with type mismatch POST", "[dual]" ) +{ + std::cout.precision(16); + constexpr int N = 3; + using D = DualS; + + MAKE_TEST_SECTION_BASIC_BINARY_OPERATOR_TYPE_MISMATCH_POST(+, 3.14, (int)2) + MAKE_TEST_SECTION_BASIC_BINARY_OPERATOR_TYPE_MISMATCH_POST(-, 3.14, (int)2) + MAKE_TEST_SECTION_BASIC_BINARY_OPERATOR_TYPE_MISMATCH_POST(*, 3.14, (int)2) + MAKE_TEST_SECTION_BASIC_BINARY_OPERATOR_TYPE_MISMATCH_POST(/, 3.14, (int)2) +} + +#define MAKE_TEST_SECTION_BASIC_BINARY_OPERATOR_TYPE_MISMATCH_PRE(Op, A, B)\ +SECTION( "operator"#Op ) { \ + D a = (A) Op D(B); \ + S resS = A Op B; \ + REQUIRE(a.a == resS); \ +} + +TEST_CASE( "Dual numbers : basic binary operators with type mismatch PRE", "[dual]" ) +{ + std::cout.precision(16); + constexpr int N = 3; + using D = DualS; + + MAKE_TEST_SECTION_BASIC_BINARY_OPERATOR_TYPE_MISMATCH_PRE(+, (int)2, 3.14) + MAKE_TEST_SECTION_BASIC_BINARY_OPERATOR_TYPE_MISMATCH_PRE(-, (int)2, 3.14) + MAKE_TEST_SECTION_BASIC_BINARY_OPERATOR_TYPE_MISMATCH_PRE(*, (int)2, 3.14) + MAKE_TEST_SECTION_BASIC_BINARY_OPERATOR_TYPE_MISMATCH_PRE(/, (int)2, 3.14) +} + +TEST_CASE( "Dual numbers : basic unary operators", "[dual]" ) +{ + std::cout.precision(16); + constexpr int N = 3; + using D = DualS; + + SECTION( "operator++" ) { + S sa = 3.14; + D a(sa); + a++; + sa++; + REQUIRE(a.a == sa); + } + + SECTION( "operator--" ) { + S sa = 3.14; + D a(sa); + a--; + sa--; + REQUIRE(a.a == sa); + } + + SECTION( "operator++" ) { + S sa = 3.14; + D a(sa); + ++a; + ++sa; + REQUIRE(a.a == sa); + } + + SECTION( "operator--" ) { + S sa = 3.14; + D a(sa); + --a; + --sa; + REQUIRE(a.a == sa); + } + + SECTION( "operator+" ) { + S sa = 3.14; + D a(sa); + a = +a; + sa = +sa; + REQUIRE(a.a == sa); + } + + SECTION( "operator-" ) { + S sa = 3.14; + D a(sa); + a = -a; + sa = -sa; + REQUIRE(a.a == sa); + } +} + +TEST_CASE( "Basic numbers : extended math library", "[dual]" ) +{ + S x = 0.500000000000000000000000000000; + REQUIRE(fabs(sec(x) - 1.139493927324549016333321560523) < tol); + REQUIRE(fabs(cot(x) - 1.830487721712451998357096272230) < tol); + REQUIRE(fabs(csc(x) - 2.085829642933488159428634389769) < tol); + REQUIRE(fabs(asec(3*x) - 0.841068670567930221082519892661) < tol); + REQUIRE(fabs(acot(x) - 1.107148717794090408972351724515) < tol); + REQUIRE(fabs(acsc(3*x) - 0.729727656226966336916461841611) < tol); + REQUIRE(fabs(sech(x) - 0.886818883970073912337284127716) < tol); + REQUIRE(fabs(coth(x) - 2.163953413738652908904214200447) < tol); + REQUIRE(fabs(csch(x) - 1.919034751334943722511638952710) < tol); + REQUIRE(fabs(asech(x) - 1.316957896924816795447554795828) < tol); + REQUIRE(fabs(acoth(3*x) - 0.804718956217050140899971211184) < tol); + REQUIRE(fabs(acsch(x) - 1.443635475178810301244425318146) < tol); + REQUIRE(fabs(exp10(x) - 3.162277660168379522787063251599) < tol); + REQUIRE(fabs(sign(x) - 1.000000000000000000000000000000) < tol); + REQUIRE(fabs(sign(-x) - -1.000000000000000000000000000000) < tol); + REQUIRE(fabs(heaviside(x) - 1.000000000000000000000000000000) < tol); + REQUIRE(fabs(heaviside(0.) - 1.000000000000000000000000000000) < tol); + REQUIRE(fabs(heaviside(-x) - 0.000000000000000000000000000000) < tol); +} + +TEST_CASE( "Dual numbers : Elementary function derivatives", "[dual]" ) +{ + S xf = 0.500000000000000000000000000000; + constexpr int N = 3; + using D = DualS; + + D x(xf); x.diff(0); + CHECK(fabs((cos(x)).d(0) - (-sin(xf))) < tol); + CHECK(fabs((sin(x)).d(0) - (cos(xf))) < tol); + CHECK(fabs((tan(x)).d(0) - (pow(tan(xf), 2) + 1)) < tol); + CHECK(fabs((acos(x)).d(0) - (-1/sqrt(-pow(xf, 2) + 1))) < tol); + CHECK(fabs((asin(x)).d(0) - (pow(-pow(xf, 2) + 1, -1.0/2.0))) < tol); + CHECK(fabs((atan(x)).d(0) - (1.0/(pow(xf, 2) + 1))) < tol); + CHECK(fabs((sec(x)).d(0) - (tan(xf)*sec(xf))) < tol); + CHECK(fabs((cot(x)).d(0) - (-pow(cot(xf), 2) - 1)) < tol); + CHECK(fabs((csc(x)).d(0) - (-cot(xf)*csc(xf))) < tol); + CHECK(fabs((asec(3*x)).d(0) - ((1.0/3.0)/(pow(xf, 2)*sqrt(1 - 1.0/9.0/pow(xf, 2))))) < tol); + CHECK(fabs((acot(x)).d(0) - (-1/(pow(xf, 2) + 1))) < tol); + CHECK(fabs((acsc(3*x)).d(0) - (-1.0/3.0/(pow(xf, 2)*sqrt(1 - 1.0/9.0/pow(xf, 2))))) < tol); + CHECK(fabs((cosh(x)).d(0) - (sinh(xf))) < tol); + CHECK(fabs((sinh(x)).d(0) - (cosh(xf))) < tol); + CHECK(fabs((atanh(x)).d(0) - (1.0/(-pow(xf, 2) + 1))) < tol); + CHECK(fabs((asinh(x)).d(0) - (pow(pow(xf, 2) + 1, -1.0/2.0))) < tol); + CHECK(fabs((acosh(3*x)).d(0) - (3/sqrt(9*pow(xf, 2) - 1))) < tol); + CHECK(fabs((atanh(x)).d(0) - (1.0/(-pow(xf, 2) + 1))) < tol); + CHECK(fabs((sech(x)).d(0) - (-tanh(xf)*sech(xf))) < tol); + CHECK(fabs((coth(x)).d(0) - (-1/pow(sinh(xf), 2))) < tol); + CHECK(fabs((csch(x)).d(0) - (-coth(xf)*csch(xf))) < tol); + CHECK(fabs((asech(x)).d(0) - (-1/(xf*sqrt(-pow(xf, 2) + 1)))) < tol); + CHECK(fabs((acoth(3*x)).d(0) - (3/(-9*pow(xf, 2) + 1))) < tol); + CHECK(fabs((acsch(x)).d(0) - (-1/(pow(xf, 2)*sqrt(1 + pow(xf, -2))))) < tol); + CHECK(fabs(exp(x).d(0) - (exp(xf))) < tol); + CHECK(fabs(exp2(x).d(0) - (pow(2, xf)*M_LN2)) < tol); + CHECK(fabs(exp10(x).d(0) - (pow(10, xf)*M_LN10)) < tol); + CHECK(fabs((log(x)).d(0) - (1.0/xf)) < tol); + CHECK(fabs((log2(x)).d(0) - (1/(xf*M_LN2))) < tol); + CHECK(fabs((log10(x)).d(0) - (1/(xf*M_LN10))) < tol); + CHECK(fabs((sqrt(x)).d(0) - ((1.0/2.0)/sqrt(xf))) < tol); + CHECK(fabs(cbrt(x).d(0) - ((1.0/3.0)/pow(xf, 2.0/3.0))) < tol); + CHECK(fabs((abs(x)).d(0) - ((((xf) > 0) - ((xf) < 0)))) < tol); + CHECK(fabs((fabs(x)).d(0) - ((((xf) > 0) - ((xf) < 0)))) < tol); + CHECK(fabs((pow(x, D(1.0/7.0))).d(0) - ((1.0/7.0)/pow(xf, 6.0/7.0))) < tol); + CHECK(fabs((pow(x,(1.0/7.0))).d(0) - ((1.0/7.0)/pow(xf, 6.0/7.0))) < tol); + CHECK(fabs((pow(D(1.0/7.0), x)).d(0) - (-pow(7, -xf)*log(7))) < tol); + CHECK(fabs((pow(1.0/7.0, x)).d(0) - (-pow(7, -xf)*log(7))) < tol); +} + +TEST_CASE( "Dual numbers : derivatives of compositions of functions ", "[dual]" ) +{ + S xf = 0.500000000000000000000000000000; + constexpr int N = 3; + using D = DualS; + + D x(xf); x.diff(0); + CHECK(fabs((pow(x, pow(x, x))).d(0) - (pow(xf, pow(xf, xf))*(pow(xf, xf)*(log(xf) + 1)*log(xf) + pow(xf, xf)/xf))) < tol); + CHECK(fabs((log(sin(x))).d(0) - (cos(xf)/sin(xf))) < tol); + CHECK(fabs((pow(x, 2)/log(x)).d(0) - (2*xf/log(xf) - xf/pow(log(xf), 2))) < tol); + CHECK(fabs((cbrt(sin(x))).d(0) - ((1.0/3.0)*cos(xf)/pow(sin(xf), 2.0/3.0))) < tol); +} + +TEST_CASE( "Dual numbers : univariate function derivative", "[dual]" ) +{ + S xf = 0.1; + constexpr int N = 3; + using D = DualS; + D x(xf); x.diff(0); + CHECK(fabs(fct1(x).x() - fct1(xf)) < tol); + CHECK(fabs(fct1(x).d(0) - dfct1(xf)) < tol); +} + +TEST_CASE( "Dual numbers : univariate function derivative using GradFunc", "[dual]" ) +{ + S xf = 0.1; + auto gradfct1 = GradFunc1(fct1>, xf); + CHECK(fabs(gradfct1(xf) - dfct1(xf)) < tol); + + S fx, dfx; + gradfct1.get_f_grad(xf, fx, dfx); + CHECK(fabs(fx - fct1(xf)) < tol); + CHECK(fabs(dfx - dfct1(xf)) < tol); +} + +TEST_CASE( "Dual numbers : multivariate function derivative", "[dual]" ) +{ + constexpr int N = 3; + using D = DualS; + + S xf = 0.11, yf = .12, zf = .13; + D x(xf); x.diff(0); + D y(yf); y.diff(1); + D z(zf); z.diff(2); + D fxyz = fct2(x, y, z); + REQUIRE(x.a == xf); + REQUIRE(y.a == yf); + REQUIRE(z.a == zf); + print(fxyz) + print(fxyz.d(0)) + print(fxyz.d(1)) + print(fxyz.d(2)) + + print(dfct2dx(xf, yf, zf)) + print(dfct2dy(xf, yf, zf)) + print(dfct2dz(xf, yf, zf)) + + CHECK(fabs(fxyz.d(0) - dfct2dx(xf, yf, zf)) < tol); + CHECK(fabs(fxyz.d(1) - dfct2dy(xf, yf, zf)) < tol); + CHECK(fabs(fxyz.d(2) - dfct2dz(xf, yf, zf)) < tol); +} + +TEST_CASE( "Dual numbers : Jacobian using Eigen arrays and Dual directly", "[dual]" ) +{ + constexpr int N = 3; + using D = DualS; + + constexpr S deg = M_PI/180.; + S r = 3, theta = 30*deg, phi = 120*deg; + Eigen::Array xsph(r, theta, phi); + xsph[0].diff(0); + xsph[1].diff(1); + xsph[2].diff(2); + Eigen::Array xcart = sph2cart(xsph); + print(xcart.transpose()) + print(xcart[0].b.transpose()) + print(xcart[1].b.transpose()) + print(xcart[2].b.transpose()) + auto Jxcart = sph2cartJ(Eigen::Array(r, theta, phi)); + std::cout << "Jxcart =\n" << (Jxcart) << std::endl; + + for(int i = 0 ; i < 3 ; i++) + for(int j = 0 ; j < 3 ; j++) + CHECK(fabs(xcart[i].d(j) - Jxcart(i,j)) < tol); } // CHECK() diff --git a/test/py_dual_test_generate_function_tests.py b/test/py_dual_test_generate_function_tests.py new file mode 100755 index 0000000..cd5dc87 --- /dev/null +++ b/test/py_dual_test_generate_function_tests.py @@ -0,0 +1,65 @@ +from sympy import * + +#%% create a list of test cases for the extended math library for automatic testing +x = S(1)/2 +print("S x = %.30f;" % x) +print("S tol = 1e-12;") +print("REQUIRE(fabs(sec(x) - %.30f) < tol);" % sec(x)) +print("REQUIRE(fabs(cot(x) - %.30f) < tol);" % cot(x)) +print("REQUIRE(fabs(csc(x) - %.30f) < tol);" % csc(x)) +print("REQUIRE(fabs(asec(3*x) - %.30f) < tol);" % asec(3*x)) +print("REQUIRE(fabs(acot(x) - %.30f) < tol);" % acot(x)) +print("REQUIRE(fabs(acsc(3*x) - %.30f) < tol);" % acsc(3*x)) +print("REQUIRE(fabs(sech(x) - %.30f) < tol);" % sech(x)) +print("REQUIRE(fabs(coth(x) - %.30f) < tol);" % coth(x)) +print("REQUIRE(fabs(csch(x) - %.30f) < tol);" % csch(x)) +print("REQUIRE(fabs(asech(x) - %.30f) < tol);" % asech(x)) +print("REQUIRE(fabs(acoth(3*x) - %.30f) < tol);" % acoth(3*x)) +print("REQUIRE(fabs(acsch(x) - %.30f) < tol);" % acsch(x)) +print("REQUIRE(fabs(exp10(x) - %.30f) < tol);" % 10**x) +print("REQUIRE(fabs(sign(x) - %.30f) < tol);" % sign(x)) +print("REQUIRE(fabs(sign(-x) - %.30f) < tol);" % sign(-x)) +print("REQUIRE(fabs(heaviside(x) - %.30f) < tol);" % 1) +print("REQUIRE(fabs(heaviside(0.) - %.30f) < tol);" % 1) +print("REQUIRE(fabs(heaviside(-x) - %.30f) < tol);" % 0) + +#%% Create list +x, xf = symbols("x xf", real=True) +# xf = S(1)/2 +print("S xf = 0.500000000000000000000000000000;") +print("D x(xf); x.diff(0);") +print("S tol = 1e-12;") +expr = (cos(x), sin(x), tan(x), acos(x), asin(x), atan(x), sec(x), cot(x), csc(x), asec(3*x), acot(x), acsc(3*x), cosh(x), sinh(x), atanh(x), asinh(x), acosh(3*x), atanh(x), sech(x), coth(x), csch(x), asech(x), acoth(3*x), acsch(x), exp(x), pow(2, x), pow(10, x), log(x), log(x)/log(2), log(x)/log(10), sqrt(x), cbrt(x), abs(x), pow(x, S(1)/7), pow(S(1)/7, x)) + +for i in range(len(expr)): + f = expr[i] + df = f.diff(x) + print("CHECK(fabs((%s).d(0) - (%s)) < tol);" % (ccode(f), ccode(df.subs(x, xf)))) + +#%% Generate small test functions +expr2 = (x**(x**x), log(sin(x)), pow(x, 2)/(log(x)), cbrt(sin(x))) +for i in range(len(expr2)): + f = expr2[i] + df = f.diff(x) + print("CHECK(fabs((%s).d(0) - (%s)) < tol);" % (ccode(f), ccode(df.subs(x, xf)))) + +#%% Generate test functions +x, y, z = symbols("x y z", real=True) +print(" ") +print("fct 1 :") +fct1 = (cos(x) + sqrt(x) + cbrt(sin(x)) + exp(-x*x)/log(x) + atanh(5*x) - 2*csc(x)) +print(ccode(fct1)) +print(ccode(fct1.diff(x))) + +print(" ") +print("fct 2 :") +fct2 = (cos(x) + sqrt(y) + cbrt(sin(z))*tan(y) + 2*sqrt(x)*sqrt(y)/log(sin(z)) + exp(-x*x - y*y - z*z)/(log(x)*log(y)*log(z)) + atanh(5*x)*2*csc(y)) +print("fct2(x,y,z) =", ccode(fct2)) +print("dfct2/dx =", ccode(fct2.diff(x))) +print("dfct2/dy =", ccode(fct2.diff(y))) +print("dfct2/dz =", ccode(fct2.diff(z))) + +fct2.subs(((x,.11),(y,.12),(z,.13))) +fct2.subs(((x,.11),(y,.12),(z,.14))) +(fct2.subs(((x,.11),(y,.12),(z,.14)))-fct2.subs(((x,.11),(y,.12),(z,.13))))/.01 +fct2.diff(z).subs(((x,.11),(y,.12),(z,.13))) diff --git a/tests/minimalistic_test_autodiff/Makefile b/tests/minimalistic_test_autodiff/Makefile new file mode 100755 index 0000000..967bbd2 --- /dev/null +++ b/tests/minimalistic_test_autodiff/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 -O0 +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/minimalistic_test_autodiff/minimalistic_test_autodiff.cpp b/tests/minimalistic_test_autodiff/minimalistic_test_autodiff.cpp new file mode 100755 index 0000000..58982ed --- /dev/null +++ b/tests/minimalistic_test_autodiff/minimalistic_test_autodiff.cpp @@ -0,0 +1,257 @@ +#include +#include + +#include +#include + +using std::cout; +using std::endl; +using namespace Eigen; + +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); +} + +template +struct Fct1 +{ + T operator()(const T & x) { return fct1(x); } +}; + +template +T fct2(const T & x, const T & y, const T & z) +{ + T two(2); + return 2*sqrt(x)*sqrt(y)/log(sin(z)) + sqrt(y) + exp(-pow(x, two) - pow(y, two) - pow(z, two))/(log(x)*log(y)*log(z)) + cbrt(sin(z))*tan(y) + cos(x) + 2*atanh(5*x)*csc(y); +} +template +T dfct2dx(const T & x, const T & y, const T & z) +{ + T two(2); + return -2*x*exp(-pow(x, two) - pow(y, two) - pow(z, two))/(log(x)*log(y)*log(z)) - sin(x) + 10*csc(y)/(-25*pow(x, two) + 1) - exp(-pow(x, two) - pow(y, two) - pow(z, two))/(x*pow(log(x), two)*log(y)*log(z)) + sqrt(y)/(sqrt(x)*log(sin(z))); +} +template +T dfct2dy(const T & x, const T & y, const T & z) +{ + T two(2); + return sqrt(x)/(sqrt(y)*log(sin(z))) - 2*y*exp(-pow(x, two) - pow(y, two) - pow(z, two))/(log(x)*log(y)*log(z)) + (pow(tan(y), two) + 1)*cbrt(sin(z)) - 2*cot(y)*atanh(5*x)*csc(y) - exp(-pow(x, two) - pow(y, two) - pow(z, two))/(y*log(x)*pow(log(y), two)*log(z)) + (1.0/2.0)/sqrt(y); +} +template +T dfct2dz(const T & x, const T & y, const T & z) +{ + T two(2); + return -2*sqrt(x)*sqrt(y)*cos(z)/(pow(log(sin(z)), two)*sin(z)) - 2*z*exp(-pow(x, two) - pow(y, two) - pow(z, two))/(log(x)*log(y)*log(z)) + (1.0/3.0)*cos(z)*tan(y)/pow(sin(z), 2.0/3.0) - exp(-pow(x, two) - pow(y, two) - pow(z, two))/(z*log(x)*log(y)*pow(log(z), two)); +} + +void test_derivatives(); +int main() +{ + {// C++14 version + { + PRINT_STR("\nStatic\n"); + { + DualS a; + PRINT_VAR(a.a); + PRINT_VAR(a.b.transpose()); + } + { + DualS a(3.14); + PRINT_VAR(a.a); + PRINT_VAR(a.b.transpose()); + } + { + DualS a(3.14, {1.,2.,3.}); + PRINT_VAR(a.a); + PRINT_VAR(a.b.transpose()); + } + { + DualS a(3.14, {1.,2.,3.}), b(2.71, {-1., 3., 4.}); + auto c = a + b; + PRINT_VAR(a.a); + PRINT_VAR(a.b.transpose()); + PRINT_VAR(b.a); + PRINT_VAR(b.b.transpose()); + PRINT_VAR(c.a); + PRINT_VAR(c.b.transpose()); + } + { + double af = 3.14, bf = 2.71; + DualS a(af, {1.,0.}), b(bf, {0., 1.}); + PRINT_VAR((a+b).a); + PRINT_VAR((a+b).b.transpose()); + PRINT_VAR((a-b).a); + PRINT_VAR((a-b).b.transpose()); + PRINT_VAR((a*b).a); + PRINT_VAR((a*b).b.transpose()); + PRINT_VAR((a/b).a); + PRINT_VAR((a/b).b.transpose()); + PRINT_VAR(((a+b)*(a+b)).a); + PRINT_VAR(((a+b)*(a+b)).b.transpose()); + PRINT_VAR((sqrt(a)*sqrt(b)).b.transpose()); + PRINT_VAR((1.0/2.0)*sqrt(bf)/sqrt(af)); + PRINT_VAR((1.0/2.0)*sqrt(af)/sqrt(bf)); + } + } + { + PRINT_STR("\nDynamic\n"); + { + DualD a; + PRINT_VAR(a.a); + PRINT_VAR(a.b.transpose()); + } + { + DualD a(3.14); + PRINT_VAR(a.a); + PRINT_VAR(a.b.transpose()); + } + { + DualD::VectorT ab(3); ab << 1.,2.,3.; + DualD a(3.14, ab); + PRINT_VAR(a.a); + PRINT_VAR(a.b.transpose()); + } + { + DualD::VectorT ab(3); ab << 1.,2.,3.; + DualD::VectorT bb(3); bb << -1.,3.,4.; + DualD a(3.14, ab), b(2.71, bb); + auto c = a + b; + PRINT_VAR(a.a); + PRINT_VAR(a.b.transpose()); + PRINT_VAR(b.a); + PRINT_VAR(b.b.transpose()); + PRINT_VAR(c.a); + PRINT_VAR(c.b.transpose()); + } + { + double af = 3.14, bf = 2.71; + DualD::VectorT ab(3); ab << 1.,2.,3.; + DualD::VectorT bb(3); bb << -1.,3.,4.; + DualD a(af, ab), b(bf, bb); + auto c = (1 + 3*a) + b*2 + 1/b + (3 - a); + PRINT_VAR(a.a); + PRINT_VAR(a.b.transpose()); + PRINT_VAR(b.a); + PRINT_VAR(b.b.transpose()); + PRINT_VAR(c.a); + PRINT_VAR(c.b.transpose()); + assert(c.a == ((1 + 3*af) + bf*2 + 1/bf + (3 - af))); + } + { + double af = 3.14, bf = 2.71; + DualD a(af), b(bf); a.diff(0); a.diff(1); + PRINT_VAR((a+b).a); + PRINT_VAR((a+b).b.transpose()); + PRINT_VAR((a-b).a); + PRINT_VAR((a-b).b.transpose()); + PRINT_VAR((a*b).a); + PRINT_VAR((a*b).b.transpose()); + PRINT_VAR((a/b).a); + PRINT_VAR((a/b).b.transpose()); + PRINT_VAR(((a+b)*(a+b)).a); + PRINT_VAR(((a+b)*(a+b)).b.transpose()); + PRINT_VAR((sqrt(a)*sqrt(b)).b.transpose()); + PRINT_VAR((1.0/2.0)*sqrt(bf)/sqrt(af)); + PRINT_VAR((1.0/2.0)*sqrt(af)/sqrt(bf)); + } + } + + #if 0 + {// should throw an compile-time assertion + Dual a; + } + #endif + #if 0 + { + DualD ad; + DualS as; + ad.b.resize(5); + as.b.resize(5);// should throw an compile-time error + } + #endif + } + + test_derivatives(); + + return 0; +} + +#define PRINT_F_DF(Term) cout << #Term << "\t: " << (Term).x() << "\t" << (Term).d(0) << "\n" + +void test_derivatives() +{ + cout.precision(16); + + using S = double; + constexpr S tol = 1e-12; + + { + PRINT_STR("--------------------------------"); + DualS x(.1); x.diff(0); + PRINT_VAR(x); + PRINT_F_DF(cos(x)); + PRINT_F_DF(sin(x)); + PRINT_F_DF(tan(x)); + PRINT_F_DF(sec(x)); + PRINT_F_DF(csc(x)); + PRINT_F_DF(cot(x)); + + PRINT_F_DF(sqrt(x)); + PRINT_F_DF(cbrt(x)); + PRINT_F_DF(exp(x)); + PRINT_F_DF(2*exp(x)); + PRINT_F_DF(exp(2*x)); + } + + #if 1 + {// univariate + PRINT_STR("--------------------------------"); + using D = DualS; + S xf = 0.1; + D x(xf); x.diff(0); + + PRINT_VAR(fct1(x)); + PRINT_VAR(fct1(xf)); + + PRINT_VAR(fct1(x).d(0)); + PRINT_VAR(dfct1(xf)); + + assert(fabs(fct1(x).x() - fct1(xf)) < tol); + assert(fabs(fct1(x).d(0) - dfct1(xf)) < tol); + } + #endif + + #if 1 + { + using D = DualS; + S xf = 0.11, yf = .12, zf = .13; + D x(xf); x.diff(0); + D y(yf); y.diff(1); + D z(zf); z.diff(2); + D fxyz = fct2(x, y, z); + assert(x.a == xf); + assert(y.a == yf); + assert(z.a == zf); + PRINT_VAR(fxyz); + PRINT_VAR(fxyz.d(0)); + PRINT_VAR(fxyz.d(1)); + PRINT_VAR(fxyz.d(2)); + + PRINT_VAR(dfct2dx(xf, yf, zf)); + PRINT_VAR(dfct2dy(xf, yf, zf)); + PRINT_VAR(dfct2dz(xf, yf, zf)); + + // CHECK(fabs(fxyz.d(0) - dfct2dx(xf, yf, zf)) < tol); + // CHECK(fabs(fxyz.d(1) - dfct2dy(xf, yf, zf)) < tol); + // CHECK(fabs(fxyz.d(2) - dfct2dz(xf, yf, zf)) < tol); + } + #endif + +} diff --git a/utils.hpp b/utils.hpp index b419ae6..3f5243e 100755 --- a/utils.hpp +++ b/utils.hpp @@ -13,13 +13,14 @@ template auto min(const T & a, const T2 & b) -> decltype(a | b) { return (a < b) ? a : b; } template auto max(const T & a, const T2 & b) -> decltype(a | b) { return (a < b) ? b : a; } +/* template class C, class... Args> std::ostream& operator<<(std::ostream& os, const C& objs) { for (auto const& obj : objs) os << obj << ' '; 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)