From 592f6976ad3fbbecdf90e2ee1adf5790906007dc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A9r=C3=B4me?= Date: Tue, 26 Mar 2019 12:40:00 +0100 Subject: [PATCH] Cleaned up folder structure, moved old stuff to old/ --- AutomaticDifferentiation.hpp | 149 ++++-- AutomaticDifferentiationVector.hpp | 498 ------------------ examples/firstAndSecondDerivative.cpp | 10 +- old/AutomaticDifferentiation_scalar.hpp | 411 +++++++++++++++ .../test_AutomaticDifferentiation.cpp | 0 .../test_AutomaticDifferentiation_manual.cpp | 0 .../test_AutomaticDifferentiation_vector.cpp | 0 ...itions_manual.cpp => test_small_manual.cpp | 34 +- 8 files changed, 551 insertions(+), 551 deletions(-) delete mode 100755 AutomaticDifferentiationVector.hpp create mode 100755 old/AutomaticDifferentiation_scalar.hpp rename test_AutomaticDifferentiation.cpp => old/test_AutomaticDifferentiation.cpp (100%) rename test_AutomaticDifferentiation_manual.cpp => old/test_AutomaticDifferentiation_manual.cpp (100%) rename test_AutomaticDifferentiation_vector.cpp => old/test_AutomaticDifferentiation_vector.cpp (100%) rename test_vector_version_additions_manual.cpp => test_small_manual.cpp (75%) diff --git a/AutomaticDifferentiation.hpp b/AutomaticDifferentiation.hpp index 153964f..2460d9d 100755 --- a/AutomaticDifferentiation.hpp +++ b/AutomaticDifferentiation.hpp @@ -1,43 +1,127 @@ #ifndef DEF_AUTOMATIC_DIFFERENTIATION #define DEF_AUTOMATIC_DIFFERENTIATION +#include #include #include -/// Implementation of dual numbers for automatic differentiation +#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: - Dual(const Scalar & _a, const Scalar & _b = Scalar(0.0)) + 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) {} - static Dual d() { - return Dual(Scalar(0.), Scalar(1.)); + 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 < MINAB(b.size(), b_old.size()) ; j++) + b[j] = b_old[j]; + } + b[i] = Scalar(1.); + return *this; } - Dual & operator+=(const Dual & x) { + /// 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) { + Dual & operator-=(const Dual & x) + { a -= x.a; b -= x.b; return *this; } - Dual & operator*=(const Dual & x) { + Dual & operator*=(const Dual & x) + { b = a*x.b + b*x.a; a *= x.a; return *this; } - Dual & operator/=(const Dual & x) { + Dual & operator/=(const Dual & x) + { b = (x.a*b - a*x.b)/(x.a*x.a); a /= x.a; return *this; @@ -68,7 +152,8 @@ class Dual return (res += x); } - Dual operator+(void) const { // +x + Dual operator+(void) const // +x + { return (*this); } @@ -77,16 +162,19 @@ class Dual return (res -= x); } - Dual operator-(void) const { // -x + Dual operator-(void) const // -x + { return Dual(-a, -b); } - Dual operator*(const Dual & x) const { + Dual operator*(const Dual & x) const + { Dual res(*this); return (res *= x); } - Dual operator/(const Dual & x) const { + Dual operator/(const Dual & x) const + { Dual res(*this); return (res /= x); } @@ -116,24 +204,24 @@ class Dual } Scalar a; /// Real part - Scalar b; /// Infinitesimal part + VectorT b; /// Infinitesimal parts }; -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); } -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 @@ -345,31 +433,30 @@ template Dual sqrt(const Dual & x) { } template Dual sign(const Dual & x) { - return Dual(sign(x.a), Scalar(0.)); + 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), Scalar(0.)); + return Dual(heaviside(x.a), Dual::Dual::__create_VectorT_zeros()); } template Dual floor(const Dual & x) { - return Dual(floor(x.a), Scalar(0.)); + return Dual(floor(x.a), Dual::Dual::__create_VectorT_zeros()); } template Dual ceil(const Dual & x) { - return Dual(ceil(x.a), Scalar(0.)); + return Dual(ceil(x.a), Dual::Dual::__create_VectorT_zeros()); } template Dual round(const Dual & x) { - return Dual(round(x.a), Scalar(0.)); + return Dual(round(x.a), Dual::Dual::__create_VectorT_zeros()); } template std::ostream & operator<<(std::ostream & s, const Dual & x) diff --git a/AutomaticDifferentiationVector.hpp b/AutomaticDifferentiationVector.hpp deleted file mode 100755 index 31fe706..0000000 --- a/AutomaticDifferentiationVector.hpp +++ /dev/null @@ -1,498 +0,0 @@ -#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+DualVector::d(0), y+DualVector::d(1), z+DualVector::d(2), ...) -/// reference : http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.89.7749&rep=rep1&type=pdf -template -class DualVector -{ - public: - using VectorT = std::valarray; - - static VectorT __create_VectorT_zeros(int N = 1) - { - assert(N >= 0); - VectorT res(Scalar(0.), N); - return res; - } - - DualVector(const Scalar & _a = Scalar(0.), const VectorT & _b = DualVector::__create_VectorT_zeros()) - : a(_a), - b(_b) - {} - - DualVector const& operator=(Scalar const& _a) - { - *this = DualVector(_a); - } - - /// Use this function to set what variable is to be derived : x + DualVector::d(i) - static DualVector D(int i = 0, int N = 1) - { - assert(i >= 0); - assert(i < N); - VectorT _d = DualVector::__create_VectorT_zeros(N); - _d[i] = Scalar(1.); - return DualVector(Scalar(0.), _d); - } - - /// Use this function to set what variable is to be derived. - DualVector 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 < MINAB(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]; - } - - DualVector & operator+=(const DualVector & x) - { - a += x.a; - b += x.b; - return *this; - } - - DualVector & operator-=(const DualVector & x) - { - a -= x.a; - b -= x.b; - return *this; - } - - DualVector & operator*=(const DualVector & x) - { - b = a*x.b + b*x.a; - a *= x.a; - return *this; - } - - DualVector & operator/=(const DualVector & x) - { - b = (x.a*b - a*x.b)/(x.a*x.a); - a /= x.a; - return *this; - } - - DualVector & operator++() { // ++x - return ((*this) += Scalar(1.)); - } - - DualVector & operator--() { // --x - return ((*this) -= Scalar(1.)); - } - - DualVector operator++(int) { // x++ - DualVector copy = *this; - (*this) += Scalar(1.); - return copy; - } - - DualVector operator--(int) { // x-- - DualVector copy = *this; - (*this) -= Scalar(1.); - return copy; - } - - DualVector operator+(const DualVector & x) const { - DualVector res(*this); - return (res += x); - } - - DualVector operator+(void) const // +x - { - return (*this); - } - - DualVector operator-(const DualVector & x) const { - DualVector res(*this); - return (res -= x); - } - - DualVector operator-(void) const // -x - { - return DualVector(-a, -b); - } - - DualVector operator*(const DualVector & x) const - { - DualVector res(*this); - return (res *= x); - } - - DualVector operator/(const DualVector & x) const - { - DualVector res(*this); - return (res /= x); - } - - bool operator==(const DualVector & x) const { - return (a == x.a); - } - - bool operator!=(const DualVector & x) const { - return (a != x.a); - } - - bool operator<(const DualVector & x) const { - return (a < x.a); - } - - bool operator<=(const DualVector & x) const { - return (a <= x.a); - } - - bool operator>(const DualVector & x) const { - return (a > x.a); - } - - bool operator>=(const DualVector & x) const { - return (a >= x.a); - } - - Scalar a; /// Real part - VectorT b; /// Infinitesimal parts -}; - -template -DualVector operator+(A const& v, DualVector const& x) { - return (DualVector(v) + x); -} -template -DualVector operator-(A const& v, DualVector const& x) { - return (DualVector(v) - x); -} -template -DualVector operator*(A const& v, DualVector const& x) { - return (DualVector(v) * x); -} -template -DualVector operator/(A const& v, DualVector const& x) { - return (DualVector(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 DualVector numbers -// f(a + b*d) = f(a) + b*f'(a)*d - -// Trigonometric functions -template DualVector cos(const DualVector & x) { - return DualVector(cos(x.a), -x.b*sin(x.a)); -} - -template DualVector sin(const DualVector & x) { - return DualVector(sin(x.a), x.b*cos(x.a)); -} - -template DualVector tan(const DualVector & x) { - return DualVector(tan(x.a), x.b*sec(x.a)*sec(x.a)); -} - -template DualVector sec(const DualVector & x) { - return DualVector(sec(x.a), x.b*sec(x.a)*tan(x.a)); -} - -template DualVector cot(const DualVector & x) { - return DualVector(cot(x.a), x.b*(-csc(x.a)*csc(x.a))); -} - -template DualVector csc(const DualVector & x) { - return DualVector(csc(x.a), x.b*(-cot(x.a)*csc(x.a))); -} - -// Inverse trigonometric functions -template DualVector acos(const DualVector & x) { - return DualVector(acos(x.a), x.b*(-Scalar(1.)/sqrt(Scalar(1.)-x.a*x.a))); -} - -template DualVector asin(const DualVector & x) { - return DualVector(asin(x.a), x.b*(Scalar(1.)/sqrt(Scalar(1.)-x.a*x.a))); -} - -template DualVector atan(const DualVector & x) { - return DualVector(atan(x.a), x.b*(Scalar(1.)/(x.a*x.a + Scalar(1.)))); -} - -template DualVector asec(const DualVector & x) { - return DualVector(asec(x.a), x.b*(Scalar(1.)/(sqrt(Scalar(1.)-Scalar(1.)/(x.a*x.a))*(x.a*x.a)))); -} - -template DualVector acot(const DualVector & x) { - return DualVector(acot(x.a), x.b*(-Scalar(1.)/((x.a*x.a)+Scalar(1.)))); -} - -template DualVector acsc(const DualVector & x) { - return DualVector(acsc(x.a), x.b*(-Scalar(1.)/(sqrt(Scalar(1.)-Scalar(1.)/(x.a*x.a))*(x.a*x.a)))); -} - -// Hyperbolic trigonometric functions -template DualVector cosh(const DualVector & x) { - return DualVector(cosh(x.a), x.b*sinh(x.a)); -} - -template DualVector sinh(const DualVector & x) { - return DualVector(sinh(x.a), x.b*cosh(x.a)); -} - -template DualVector tanh(const DualVector & x) { - return DualVector(tanh(x.a), x.b*sech(x.a)*sech(x.a)); -} - -template DualVector sech(const DualVector & x) { - return DualVector(sech(x.a), x.b*(-sech(x.a)*tanh(x.a))); -} - -template DualVector coth(const DualVector & x) { - return DualVector(coth(x.a), x.b*(-csch(x.a)*csch(x.a))); -} - -template DualVector csch(const DualVector & x) { - return DualVector(csch(x.a), x.b*(-coth(x.a)*csch(x.a))); -} - -// Inverse hyperbolic trigonometric functions -template DualVector acosh(const DualVector & x) { - return DualVector(acosh(x.a), x.b*(Scalar(1.)/sqrt((x.a*x.a)-Scalar(1.)))); -} - -template DualVector asinh(const DualVector & x) { - return DualVector(asinh(x.a), x.b*(Scalar(1.)/sqrt((x.a*x.a)+Scalar(1.)))); -} - -template DualVector atanh(const DualVector & x) { - return DualVector(atanh(x.a), x.b*(Scalar(1.)/(Scalar(1.)-(x.a*x.a)))); -} - -template DualVector asech(const DualVector & x) { - return DualVector(asech(x.a), x.b*(Scalar(-1.)/(sqrt(Scalar(1.)/(x.a*x.a)-Scalar(1.))*(x.a*x.a)))); -} - -template DualVector acoth(const DualVector & x) { - return DualVector(acoth(x.a), x.b*(-Scalar(1.)/((x.a*x.a)-Scalar(1.)))); -} - -template DualVector acsch(const DualVector & x) { - return DualVector(acsch(x.a), x.b*(-Scalar(1.)/(sqrt(Scalar(1.)/(x.a*x.a)+Scalar(1.))*(x.a*x.a)))); -} - -// Exponential functions -template DualVector exp(const DualVector & x) { - return DualVector(exp(x.a), x.b*exp(x.a)); -} - -template DualVector log(const DualVector & x) { - return DualVector(log(x.a), x.b/x.a); -} - -template DualVector exp10(const DualVector & x) { - return DualVector(exp10(x.a), x.b*(log(Scalar(10.))*exp10(x.a))); -} - -template DualVector log10(const DualVector & x) { - return DualVector(log10(x.a), x.b/(log(Scalar(10.))*x.a)); -} - -template DualVector exp2(const DualVector & x) { - return DualVector(exp2(x.a), x.b*(log(Scalar(2.))*exp2(x.a))); -} - -template DualVector log2(const DualVector & x) { - return DualVector(log2(x.a), x.b/(log(Scalar(2.))*x.a)); -} - -template DualVector pow(const DualVector & x, const DualVector & n) { - return exp(n*log(x)); -} - -// Other functions -template DualVector sqrt(const DualVector & x) { - return DualVector(sqrt(x.a), x.b/(Scalar(2.)*sqrt(x.a))); -} - -template DualVector sign(const DualVector & x) { - return DualVector(sign(x.a), DualVector::DualVector::__create_VectorT_zeros()); -} - -template DualVector abs(const DualVector & x) { - return DualVector(abs(x.a), x.b*sign(x.a)); -} -template DualVector fabs(const DualVector & x) { - return DualVector(fabs(x.a), x.b*sign(x.a)); -} - -template DualVector heaviside(const DualVector & x) { - return DualVector(heaviside(x.a), DualVector::DualVector::__create_VectorT_zeros()); -} - -template DualVector floor(const DualVector & x) { - return DualVector(floor(x.a), DualVector::DualVector::__create_VectorT_zeros()); -} - -template DualVector ceil(const DualVector & x) { - return DualVector(ceil(x.a), DualVector::DualVector::__create_VectorT_zeros()); -} - -template DualVector round(const DualVector & x) { - return DualVector(round(x.a), DualVector::DualVector::__create_VectorT_zeros()); -} - -template std::ostream & operator<<(std::ostream & s, const DualVector & x) -{ - return (s << x.a); -} - -#endif - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/examples/firstAndSecondDerivative.cpp b/examples/firstAndSecondDerivative.cpp index af4aae5..29aa8c8 100755 --- a/examples/firstAndSecondDerivative.cpp +++ b/examples/firstAndSecondDerivative.cpp @@ -1,4 +1,4 @@ -#include "../AutomaticDifferentiationVector.hpp" +#include "../AutomaticDifferentiation.hpp" #include #include @@ -39,7 +39,7 @@ int main() // 1st derivative forward { - using Fd = DualVector; + using Fd = Dual; Fd x = xdbl; x.diff(0); Fd y = f(x); @@ -51,10 +51,10 @@ int main() // 2nd derivative forward /* { - using Fdd = DualVector,2>; + using Fdd = Dual>; Fdd x(xdbl); - x.diff(0); - x.a.diff(1); + x.diff(0,2); + x.x().diff(1,2); Fdd y = f(x); cout << "\nForward 2nd der\n"; cout << "f(x) = " << y.a.a << endl; diff --git a/old/AutomaticDifferentiation_scalar.hpp b/old/AutomaticDifferentiation_scalar.hpp new file mode 100755 index 0000000..153964f --- /dev/null +++ b/old/AutomaticDifferentiation_scalar.hpp @@ -0,0 +1,411 @@ +#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/test_AutomaticDifferentiation.cpp b/old/test_AutomaticDifferentiation.cpp similarity index 100% rename from test_AutomaticDifferentiation.cpp rename to old/test_AutomaticDifferentiation.cpp diff --git a/test_AutomaticDifferentiation_manual.cpp b/old/test_AutomaticDifferentiation_manual.cpp similarity index 100% rename from test_AutomaticDifferentiation_manual.cpp rename to old/test_AutomaticDifferentiation_manual.cpp diff --git a/test_AutomaticDifferentiation_vector.cpp b/old/test_AutomaticDifferentiation_vector.cpp similarity index 100% rename from test_AutomaticDifferentiation_vector.cpp rename to old/test_AutomaticDifferentiation_vector.cpp diff --git a/test_vector_version_additions_manual.cpp b/test_small_manual.cpp similarity index 75% rename from test_vector_version_additions_manual.cpp rename to test_small_manual.cpp index a6bfbfd..3a4946b 100755 --- a/test_vector_version_additions_manual.cpp +++ b/test_small_manual.cpp @@ -1,4 +1,4 @@ -#include "AutomaticDifferentiationVector.hpp" +#include "AutomaticDifferentiation.hpp" #include #include @@ -53,13 +53,13 @@ int main() { // 1 var function { - DualVector x = 1.54; + Dual x = 1.54; PRINT_DUAL(x); x.diff(0); PRINT_DUAL(x); - DualVector y = f1(x); + Dual y = f1(x); PRINT_VAR(x); PRINT_DUAL(y); @@ -71,11 +71,11 @@ int main() // 3 var function { - DualVector x(2.), y(3.), z(-1.5); + Dual x(2.), y(3.), z(-1.5); x.diff(0, 3); y.diff(1, 3); z.diff(2, 3); - DualVector res = f3(x, y, z); + Dual res = f3(x, y, z); PRINT_DUAL(res); PRINT_VAR(res.b[0]); @@ -92,7 +92,7 @@ int main() // test d() diff and D { - DualVector x;// by default, size(b)=1 + Dual x;// by default, size(b)=1 assert(x.b.size() == 1); // x.d(3);// assertion thrown x.d(0);// good @@ -109,20 +109,20 @@ int main() assert(x.b[2] == 1.); assert(x.b[2] == 1.); - assert(DualVector::D(0).b.size() == 1); - assert(DualVector::D(0).b[0] == 1.); + assert(Dual::D(0).b.size() == 1); + assert(Dual::D(0).b[0] == 1.); - assert(DualVector::D(0, 3).b[0] == 1.); - assert(DualVector::D(0, 3).b[1] == 0.); - assert(DualVector::D(0, 3).b[2] == 0.); + assert(Dual::D(0, 3).b[0] == 1.); + assert(Dual::D(0, 3).b[1] == 0.); + assert(Dual::D(0, 3).b[2] == 0.); - assert(DualVector::D(1, 3).b[0] == 0.); - assert(DualVector::D(1, 3).b[1] == 1.); - assert(DualVector::D(1, 3).b[2] == 0.); + assert(Dual::D(1, 3).b[0] == 0.); + assert(Dual::D(1, 3).b[1] == 1.); + assert(Dual::D(1, 3).b[2] == 0.); - assert(DualVector::D(2, 3).b[0] == 0.); - assert(DualVector::D(2, 3).b[1] == 0.); - assert(DualVector::D(2, 3).b[2] == 1.); + assert(Dual::D(2, 3).b[0] == 0.); + assert(Dual::D(2, 3).b[1] == 0.); + assert(Dual::D(2, 3).b[2] == 1.); } return 0;