diff --git a/AutomaticDifferentiation.hpp b/AutomaticDifferentiation.hpp new file mode 100644 index 0000000..a8af999 --- /dev/null +++ b/AutomaticDifferentiation.hpp @@ -0,0 +1,390 @@ +#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 +}; + +// 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.)); +} + +// 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(fabs(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/AutomaticDifferentiationVector.hpp b/AutomaticDifferentiationVector.hpp new file mode 100644 index 0000000..e6455df --- /dev/null +++ b/AutomaticDifferentiationVector.hpp @@ -0,0 +1,367 @@ +#ifndef DEF_AUTOMATIC_DIFFERENTIATION +#define DEF_AUTOMATIC_DIFFERENTIATION + +#include +#include + +#include + +/// 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 = Eigen::Matrix; + + DualVector(const Scalar & _a, const VectorT & _b = VectorT::Zero()) + : a(_a), + b(_b) + {} + + DualVector(const Scalar & _a, const Scalar & _b) + : a(_a) + { + b.fill(_b); + } + + static DualVector d(int i = 0) + { + assert(i >= 0); + assert(i < N); + VectorT _d = VectorT::Zero(); + _d[i] = Scalar(1.); + return DualVector(Scalar(0.), _d); + } + + 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+(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); + } + + Scalar a; /// Real part + VectorT b; /// Infinitesimal parts +}; + +// 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); +} + +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.)); +} + +// 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::VectorT::Zero()); +} + +template DualVector abs(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::VectorT::Zero()); +} + +template DualVector floor(const DualVector & x) { + return DualVector(floor(x.a), DualVector::VectorT::Zero()); +} + +template DualVector ceil(const DualVector & x) { + return DualVector(ceil(x.a), DualVector::VectorT::Zero()); +} + +template DualVector round(const DualVector & x) { + return DualVector(round(x.a), DualVector::VectorT::Zero()); +} + +template std::ostream & operator<<(std::ostream & s, const DualVector & x) +{ + return (s << x.a); +} + +#endif + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/Makefile.linux b/Makefile.linux new file mode 100644 index 0000000..6b64939 --- /dev/null +++ b/Makefile.linux @@ -0,0 +1,36 @@ +CXX = g++ +INCLUDE = -I../third_party_libs/Catch2-master/single_include +CFLAGS = -Wall -std=c++11 -O0 -fprofile-arcs -ftest-coverage $(INCLUDE) +LDFLAGS = -lgcov +LDFLAGS_ALL = $(LDFLAGS) +OUTPUT_NAME = test_AutomaticDifferentiation + +# scalar version +test_AutomaticDifferentiation:test_AutomaticDifferentiation.o test_AutomaticDifferentiation_main.o + $(CXX) -o $(OUTPUT_NAME) test_AutomaticDifferentiation.o test_AutomaticDifferentiation_main.o $(LDFLAGS_ALL) + +test_AutomaticDifferentiation_manual:test_AutomaticDifferentiation_manual.o + $(CXX) -o $(OUTPUT_NAME) test_AutomaticDifferentiation_manual.o $(LDFLAGS_ALL) + +test_AutomaticDifferentiation_main.o: test_AutomaticDifferentiation_main.cpp + $(CXX) $(CFLAGS) $(LDFLAGS) -c test_AutomaticDifferentiation_main.cpp + +test_AutomaticDifferentiation.o: test_AutomaticDifferentiation.cpp + $(CXX) $(CFLAGS) $(LDFLAGS) -c test_AutomaticDifferentiation.cpp + +test_AutomaticDifferentiation_manual.o: test_AutomaticDifferentiation_manual.cpp + $(CXX) $(CFLAGS) $(LDFLAGS) -c test_AutomaticDifferentiation_manual.cpp + +# Vector version +test_AutomaticDifferentiation_vector:test_AutomaticDifferentiation_vector.o + $(CXX) -o test_AutomaticDifferentiation_vector.exe test_AutomaticDifferentiation_vector.o $(LDFLAGS_ALL) + +test_AutomaticDifferentiation_vector.o: test_AutomaticDifferentiation_vector.cpp + $(CXX) $(CFLAGS) $(LDFLAGS) -c test_AutomaticDifferentiation_vector.cpp + +clean: + rm *.o + +cleaner: + rm *.o + rm $(OUTPUT_NAME) diff --git a/Makefile.windows b/Makefile.windows new file mode 100644 index 0000000..4a11c55 --- /dev/null +++ b/Makefile.windows @@ -0,0 +1,36 @@ +CXX = g++ +INCLUDE = -I D:/Users/jerome/Documents/Programmation/third_party_libs/Eigen3 -ID:/Users/jerome/Documents/Programmation/third_party_libs/Catch2-master/single_include +CFLAGS = -Wall -std=c++11 $(INCLUDE) +LDFLAGS = +LDFLAGS_ALL = $(LDFLAGS) +OUTPUT_NAME = test_AutomaticDifferentiation + +# scalar version +test_AutomaticDifferentiation:test_AutomaticDifferentiation.o test_AutomaticDifferentiation_main.o + $(CXX) -o $(OUTPUT_NAME) test_AutomaticDifferentiation.o test_AutomaticDifferentiation_main.o $(LDFLAGS_ALL) + +test_AutomaticDifferentiation_manual:test_AutomaticDifferentiation_manual.o + $(CXX) -o $(OUTPUT_NAME) test_AutomaticDifferentiation_manual.o $(LDFLAGS_ALL) + +test_AutomaticDifferentiation_main.o: test_AutomaticDifferentiation_main.cpp + $(CXX) $(CFLAGS) $(LDFLAGS) -c test_AutomaticDifferentiation_main.cpp + +test_AutomaticDifferentiation.o: test_AutomaticDifferentiation.cpp + $(CXX) $(CFLAGS) $(LDFLAGS) -c test_AutomaticDifferentiation.cpp + +test_AutomaticDifferentiation_manual.o: test_AutomaticDifferentiation_manual.cpp + $(CXX) $(CFLAGS) $(LDFLAGS) -c test_AutomaticDifferentiation_manual.cpp + +# Vector version +test_AutomaticDifferentiation_vector:test_AutomaticDifferentiation_vector.o + $(CXX) -o test_AutomaticDifferentiation_vector.exe test_AutomaticDifferentiation_vector.o $(LDFLAGS_ALL) + +test_AutomaticDifferentiation_vector.o: test_AutomaticDifferentiation_vector.cpp + $(CXX) $(CFLAGS) $(LDFLAGS) -c test_AutomaticDifferentiation_vector.cpp + +clean: + rm *.o + +cleaner: + rm *.o + rm $(OUTPUT_NAME) diff --git a/coverage_test.txt b/coverage_test.txt new file mode 100644 index 0000000..40c9207 --- /dev/null +++ b/coverage_test.txt @@ -0,0 +1,11 @@ +To help the next person with a similar question, here's what I ended up doing (thanks to /u/flyingcaribou): +1- I'm using cmake to build my project and Catch for unit tests. +2- In the root CMakeLists.txt I create multiple targets, including a test executable (tests.exe)which will run all unit tests. The compile flags are -O0 -g --coverage -std=c++14. +3- ./tests.exe +4- In the build directory under ./CMakeFiles/, each target has its own subdirectory. In there (e.g. under build/CMakeFiles/targetName/src) a file with gcno extension will be created as the result of running tests. I used lcov to parse these coverage files. +5- brew install lcov +6- lcov -c -d CMakeFiles -o cov.info +7- genhtml cov.info -o out +8- open out/index.html + +https://medium.com/@naveen.maltesh/generating-code-coverage-report-using-gnu-gcov-lcov-ee54a4de3f11 diff --git a/test_AutomaticDifferentiation.cpp b/test_AutomaticDifferentiation.cpp new file mode 100644 index 0000000..f3ddcce --- /dev/null +++ b/test_AutomaticDifferentiation.cpp @@ -0,0 +1,353 @@ +#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/test_AutomaticDifferentiation_main.cpp b/test_AutomaticDifferentiation_main.cpp new file mode 100644 index 0000000..0465d9b --- /dev/null +++ b/test_AutomaticDifferentiation_main.cpp @@ -0,0 +1,4 @@ +// Let Catch provide main() +#define CATCH_CONFIG_MAIN + +#include diff --git a/test_AutomaticDifferentiation_manual.cpp b/test_AutomaticDifferentiation_manual.cpp new file mode 100644 index 0000000..1f4fa93 --- /dev/null +++ b/test_AutomaticDifferentiation_manual.cpp @@ -0,0 +1,503 @@ +#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/test_AutomaticDifferentiation_vector.cpp b/test_AutomaticDifferentiation_vector.cpp new file mode 100644 index 0000000..e967687 --- /dev/null +++ b/test_AutomaticDifferentiation_vector.cpp @@ -0,0 +1,461 @@ +#include "AutomaticDifferentiationVector.hpp" + +#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_DUALVECTOR(x) std::cout << #x << "\t= " << std::fixed << std::setprecision(4) << std::setw(10) << (x).a << ", " << std::setw(10) << (x).b.transpose() << 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(abs((a) - (b)) <= TEST_EQ_TOL) + +#define TEST_FUNCTION_ON_DualVector_DOUBLE(fct, x) \ + if(abs((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); + + double x = 2; + double y = 5; + + D X(x), Y(y); + + assert(X.a == x); + assert(X.b == D::VectorT::Zero()); + assert(Y.a == y); + assert(Y.b == D::VectorT::Zero()); + 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 = ((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 = 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.transpose()); + + 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)); +} + + + + + + + + + + + + + + + + + + + + + + + + + + +