From 864e633552d591f0ea9c893359530bb6e2fad1cf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A9r=C3=B4me?= Date: Tue, 26 Mar 2019 12:30:17 +0100 Subject: [PATCH] Fixed vector version to look more like FADBAD in terms of interface --- .gitignore | 4 + AutomaticDifferentiation.hpp | 0 AutomaticDifferentiationVector.hpp | 280 ++++++++++++++--------- Makefile | 0 Makefile.linux | 0 Makefile.windows | 0 README.md | 0 coverage_test.txt | 0 examples/firstAndSecondDerivative.cpp | 0 test_AutomaticDifferentiation.cpp | 0 test_AutomaticDifferentiation_main.cpp | 0 test_AutomaticDifferentiation_manual.cpp | 0 test_AutomaticDifferentiation_vector.cpp | 0 test_vector_version_additions_manual.cpp | 49 +++- 14 files changed, 220 insertions(+), 113 deletions(-) mode change 100644 => 100755 .gitignore mode change 100644 => 100755 AutomaticDifferentiation.hpp mode change 100644 => 100755 AutomaticDifferentiationVector.hpp mode change 100644 => 100755 Makefile mode change 100644 => 100755 Makefile.linux mode change 100644 => 100755 Makefile.windows mode change 100644 => 100755 README.md mode change 100644 => 100755 coverage_test.txt mode change 100644 => 100755 examples/firstAndSecondDerivative.cpp mode change 100644 => 100755 test_AutomaticDifferentiation.cpp mode change 100644 => 100755 test_AutomaticDifferentiation_main.cpp mode change 100644 => 100755 test_AutomaticDifferentiation_manual.cpp mode change 100644 => 100755 test_AutomaticDifferentiation_vector.cpp mode change 100644 => 100755 test_vector_version_additions_manual.cpp diff --git a/.gitignore b/.gitignore old mode 100644 new mode 100755 index 175b764..762afd8 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,10 @@ # executables run +# catch2 generated files +*.gcda +*.gcno + # Prerequisites *.d diff --git a/AutomaticDifferentiation.hpp b/AutomaticDifferentiation.hpp old mode 100644 new mode 100755 diff --git a/AutomaticDifferentiationVector.hpp b/AutomaticDifferentiationVector.hpp old mode 100644 new mode 100755 index c95d78c..31fe706 --- a/AutomaticDifferentiationVector.hpp +++ b/AutomaticDifferentiationVector.hpp @@ -15,63 +15,87 @@ std::ostream & operator<<(std::ostream & out, std::valarray const& v) 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 +template class DualVector { public: using VectorT = std::valarray; - static VectorT __create_VectorT_zeros() + static VectorT __create_VectorT_zeros(int N = 1) { + assert(N >= 0); VectorT res(Scalar(0.), N); return res; } - DualVector(const Scalar & _a, const VectorT & _b = DualVector::__create_VectorT_zeros()) + DualVector(const Scalar & _a = Scalar(0.), const VectorT & _b = DualVector::__create_VectorT_zeros()) : a(_a), b(_b) {} - DualVector(const Scalar & _a, const Scalar & _b) - : a(_a) - { - b = VectorT(_b, N); - } - 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) + static DualVector D(int i = 0, int N = 1) { assert(i >= 0); assert(i < N); - VectorT _d = DualVector::__create_VectorT_zeros(); + 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) + DualVector const& diff(int i = 0, int N = 1) { assert(i >= 0); assert(i < N); - b = DualVector::__create_VectorT_zeros(); + 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 < N); + 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]; } @@ -103,8 +127,27 @@ class DualVector return *this; } - DualVector operator+(const DualVector & x) const - { + 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); } @@ -136,25 +179,49 @@ class DualVector 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); } -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 @@ -211,6 +278,7 @@ 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.))); } @@ -231,167 +299,167 @@ template Scalar abs(const Scalar & x) { // 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 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 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 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 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 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))); +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 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 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 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 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 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)))); +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 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 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 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 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 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))); +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 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 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 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 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 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)))); +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 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 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 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 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 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 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) { +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 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 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 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 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 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 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 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 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) +template std::ostream & operator<<(std::ostream & s, const DualVector & x) { return (s << x.a); } diff --git a/Makefile b/Makefile old mode 100644 new mode 100755 diff --git a/Makefile.linux b/Makefile.linux old mode 100644 new mode 100755 diff --git a/Makefile.windows b/Makefile.windows old mode 100644 new mode 100755 diff --git a/README.md b/README.md old mode 100644 new mode 100755 diff --git a/coverage_test.txt b/coverage_test.txt old mode 100644 new mode 100755 diff --git a/examples/firstAndSecondDerivative.cpp b/examples/firstAndSecondDerivative.cpp old mode 100644 new mode 100755 diff --git a/test_AutomaticDifferentiation.cpp b/test_AutomaticDifferentiation.cpp old mode 100644 new mode 100755 diff --git a/test_AutomaticDifferentiation_main.cpp b/test_AutomaticDifferentiation_main.cpp old mode 100644 new mode 100755 diff --git a/test_AutomaticDifferentiation_manual.cpp b/test_AutomaticDifferentiation_manual.cpp old mode 100644 new mode 100755 diff --git a/test_AutomaticDifferentiation_vector.cpp b/test_AutomaticDifferentiation_vector.cpp old mode 100644 new mode 100755 diff --git a/test_vector_version_additions_manual.cpp b/test_vector_version_additions_manual.cpp old mode 100644 new mode 100755 index f0c25e8..a6bfbfd --- a/test_vector_version_additions_manual.cpp +++ b/test_vector_version_additions_manual.cpp @@ -53,13 +53,13 @@ int main() { // 1 var function { - DualVector x = 1.54; + DualVector x = 1.54; PRINT_DUAL(x); x.diff(0); PRINT_DUAL(x); - DualVector y = f1(x); + DualVector 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); - x.diff(0); - y.diff(1); - z.diff(2); - DualVector res = f3(x, y, z); + DualVector 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); PRINT_DUAL(res); PRINT_VAR(res.b[0]); @@ -90,5 +90,40 @@ int main() TEST_EQ_DOUBLE(res.b[2], df3z(x.a, y.a, z.a)); } + // test d() diff and D + { + DualVector x;// by default, size(b)=1 + assert(x.b.size() == 1); + // x.d(3);// assertion thrown + x.d(0);// good + assert(x.d(0) == x.b[0]); + + // assignment through d() + x.d(0) = 1.;// good + assert(x.d(0) == 1.); + + // b resizing through d + x.diff(2, 3);// good + assert(x.b[0] == 1.); + assert(x.b[1] == 0.); + 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(DualVector::D(0, 3).b[0] == 1.); + assert(DualVector::D(0, 3).b[1] == 0.); + assert(DualVector::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(DualVector::D(2, 3).b[0] == 0.); + assert(DualVector::D(2, 3).b[1] == 0.); + assert(DualVector::D(2, 3).b[2] == 1.); + } + return 0; }