jerome/AutomaticDifferentiation
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Cleaned up folder structure, moved old stuff to old/

Browse source
This commit is contained in:
Jérôme 2019-03-26 12:40:00 +01:00
parent 864e633552
commit 592f6976ad
8 changed files with 551 additions and 551 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

149
AutomaticDifferentiation.hpp
View file

@ -1,43 +1,127 @@
#ifndef DEF_AUTOMATIC_DIFFERENTIATION
#define DEF_AUTOMATIC_DIFFERENTIATION
#include <assert.h>
#include <cmath>
#include <ostream>
/// Implementation of dual numbers for automatic differentiation
#include <valarray>
template<typename T>
std::ostream & operator<<(std::ostream & out, std::valarray<T> 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<typename Scalar>
class Dual
{
public:
Dual(const Scalar & _a, const Scalar & _b = Scalar(0.0))
using VectorT = std::valarray<Scalar>;
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<typename A, typename B, int N>
Dual<B, N> operator+(A const& v, Dual<B, N> const& x) {
return (Dual<B, N>(v) + x);
template<typename A, typename B>
Dual<B> operator+(A const& v, Dual<B> const& x) {
return (Dual<B>(v) + x);
}
template<typename A, typename B, int N>
Dual<B, N> operator-(A const& v, Dual<B, N> const& x) {
return (Dual<B, N>(v) - x);
template<typename A, typename B>
Dual<B> operator-(A const& v, Dual<B> const& x) {
return (Dual<B>(v) - x);
}
template<typename A, typename B, int N>
Dual<B, N> operator*(A const& v, Dual<B, N> const& x) {
return (Dual<B, N>(v) * x);
template<typename A, typename B>
Dual<B> operator*(A const& v, Dual<B> const& x) {
return (Dual<B>(v) * x);
}
template<typename A, typename B, int N>
Dual<B, N> operator/(A const& v, Dual<B, N> const& x) {
return (Dual<B, N>(v) / x);
template<typename A, typename B>
Dual<B> operator/(A const& v, Dual<B> const& x) {
return (Dual<B>(v) / x);
}
// Basic mathematical functions for Scalar numbers
@ -345,31 +433,30 @@ template<typename Scalar> Dual<Scalar> sqrt(const Dual<Scalar> & x) {
}
template<typename Scalar> Dual<Scalar> sign(const Dual<Scalar> & x) {
return Dual<Scalar>(sign(x.a), Scalar(0.));
return Dual<Scalar>(sign(x.a), Dual<Scalar>::Dual::__create_VectorT_zeros());
}
template<typename Scalar> Dual<Scalar> abs(const Dual<Scalar> & x) {
return Dual<Scalar>(abs(x.a), x.b*sign(x.a));
}
template<typename Scalar> Dual<Scalar> fabs(const Dual<Scalar> & x) {
return Dual<Scalar>(fabs(x.a), x.b*sign(x.a));
}
template<typename Scalar> Dual<Scalar> heaviside(const Dual<Scalar> & x) {
return Dual<Scalar>(heaviside(x.a), Scalar(0.));
return Dual<Scalar>(heaviside(x.a), Dual<Scalar>::Dual::__create_VectorT_zeros());
}
template<typename Scalar> Dual<Scalar> floor(const Dual<Scalar> & x) {
return Dual<Scalar>(floor(x.a), Scalar(0.));
return Dual<Scalar>(floor(x.a), Dual<Scalar>::Dual::__create_VectorT_zeros());
}
template<typename Scalar> Dual<Scalar> ceil(const Dual<Scalar> & x) {
return Dual<Scalar>(ceil(x.a), Scalar(0.));
return Dual<Scalar>(ceil(x.a), Dual<Scalar>::Dual::__create_VectorT_zeros());
}
template<typename Scalar> Dual<Scalar> round(const Dual<Scalar> & x) {
return Dual<Scalar>(round(x.a), Scalar(0.));
return Dual<Scalar>(round(x.a), Dual<Scalar>::Dual::__create_VectorT_zeros());
}
template<typename Scalar> std::ostream & operator<<(std::ostream & s, const Dual<Scalar> & x)

498
AutomaticDifferentiationVector.hpp
View file

@ -1,498 +0,0 @@
#ifndef DEF_AUTOMATIC_DIFFERENTIATION
#define DEF_AUTOMATIC_DIFFERENTIATION
#include <assert.h>
#include <cmath>
#include <ostream>
#include <valarray>
template<typename T>
std::ostream & operator<<(std::ostream & out, std::valarray<T> 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<typename Scalar>
class DualVector
{
public:
using VectorT = std::valarray<Scalar>;
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<typename A, typename B>
DualVector<B> operator+(A const& v, DualVector<B> const& x) {
return (DualVector<B>(v) + x);
}
template<typename A, typename B>
DualVector<B> operator-(A const& v, DualVector<B> const& x) {
return (DualVector<B>(v) - x);
}
template<typename A, typename B>
DualVector<B> operator*(A const& v, DualVector<B> const& x) {
return (DualVector<B>(v) * x);
}
template<typename A, typename B>
DualVector<B> operator/(A const& v, DualVector<B> const& x) {
return (DualVector<B>(v) / x);
}
// Basic mathematical functions for Scalar numbers
// Trigonometric functions
template<typename Scalar> Scalar sec(const Scalar & x) {
return Scalar(1.)/cos(x);
}
template<typename Scalar> Scalar cot(const Scalar & x) {
return cos(x)/sin(x);
}
template<typename Scalar> Scalar csc(const Scalar & x) {
return Scalar(1.)/sin(x);
}
// Inverse trigonometric functions
template<typename Scalar> Scalar asec(const Scalar & x) {
return acos(Scalar(1.)/x);
}
template<typename Scalar> Scalar acot(const Scalar & x) {
return atan(Scalar(1.)/x);
}
template<typename Scalar> Scalar acsc(const Scalar & x) {
return asin(Scalar(1.)/x);
}
// Hyperbolic trigonometric functions
template<typename Scalar> Scalar sech(const Scalar & x) {
return Scalar(1.)/cosh(x);
}
template<typename Scalar> Scalar coth(const Scalar & x) {
return cosh(x)/sinh(x);
}
template<typename Scalar> Scalar csch(const Scalar & x) {
return Scalar(1.)/sinh(x);
}
// Inverse hyperbolic trigonometric functions
template<typename Scalar> Scalar asech(const Scalar & x) {
return log((Scalar(1.) + sqrt(Scalar(1.) - x*x))/x);
}
template<typename Scalar> Scalar acoth(const Scalar & x) {
return Scalar(0.5)*log((x + Scalar(1.))/(x - Scalar(1.)));
}
template<typename Scalar> 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<typename Scalar> Scalar exp10(const Scalar & x) {
return exp(x*log(Scalar(10.)));
}
template<typename Scalar> Scalar sign(const Scalar & x) {
return (x >= Scalar(0.)) ? ((x > Scalar(0.)) ? Scalar(1.) : Scalar(0.)) : Scalar(-1.);
}
template<typename Scalar> Scalar heaviside(const Scalar & x) {
return Scalar(x >= Scalar(0.));
}
template<typename Scalar> 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<typename Scalar> DualVector<Scalar> cos(const DualVector<Scalar> & x) {
return DualVector<Scalar>(cos(x.a), -x.b*sin(x.a));
}
template<typename Scalar> DualVector<Scalar> sin(const DualVector<Scalar> & x) {
return DualVector<Scalar>(sin(x.a), x.b*cos(x.a));
}
template<typename Scalar> DualVector<Scalar> tan(const DualVector<Scalar> & x) {
return DualVector<Scalar>(tan(x.a), x.b*sec(x.a)*sec(x.a));
}
template<typename Scalar> DualVector<Scalar> sec(const DualVector<Scalar> & x) {
return DualVector<Scalar>(sec(x.a), x.b*sec(x.a)*tan(x.a));
}
template<typename Scalar> DualVector<Scalar> cot(const DualVector<Scalar> & x) {
return DualVector<Scalar>(cot(x.a), x.b*(-csc(x.a)*csc(x.a)));
}
template<typename Scalar> DualVector<Scalar> csc(const DualVector<Scalar> & x) {
return DualVector<Scalar>(csc(x.a), x.b*(-cot(x.a)*csc(x.a)));
}
// Inverse trigonometric functions
template<typename Scalar> DualVector<Scalar> acos(const DualVector<Scalar> & x) {
return DualVector<Scalar>(acos(x.a), x.b*(-Scalar(1.)/sqrt(Scalar(1.)-x.a*x.a)));
}
template<typename Scalar> DualVector<Scalar> asin(const DualVector<Scalar> & x) {
return DualVector<Scalar>(asin(x.a), x.b*(Scalar(1.)/sqrt(Scalar(1.)-x.a*x.a)));
}
template<typename Scalar> DualVector<Scalar> atan(const DualVector<Scalar> & x) {
return DualVector<Scalar>(atan(x.a), x.b*(Scalar(1.)/(x.a*x.a + Scalar(1.))));
}
template<typename Scalar> DualVector<Scalar> asec(const DualVector<Scalar> & x) {
return DualVector<Scalar>(asec(x.a), x.b*(Scalar(1.)/(sqrt(Scalar(1.)-Scalar(1.)/(x.a*x.a))*(x.a*x.a))));
}
template<typename Scalar> DualVector<Scalar> acot(const DualVector<Scalar> & x) {
return DualVector<Scalar>(acot(x.a), x.b*(-Scalar(1.)/((x.a*x.a)+Scalar(1.))));
}
template<typename Scalar> DualVector<Scalar> acsc(const DualVector<Scalar> & x) {
return DualVector<Scalar>(acsc(x.a), x.b*(-Scalar(1.)/(sqrt(Scalar(1.)-Scalar(1.)/(x.a*x.a))*(x.a*x.a))));
}
// Hyperbolic trigonometric functions
template<typename Scalar> DualVector<Scalar> cosh(const DualVector<Scalar> & x) {
return DualVector<Scalar>(cosh(x.a), x.b*sinh(x.a));
}
template<typename Scalar> DualVector<Scalar> sinh(const DualVector<Scalar> & x) {
return DualVector<Scalar>(sinh(x.a), x.b*cosh(x.a));
}
template<typename Scalar> DualVector<Scalar> tanh(const DualVector<Scalar> & x) {
return DualVector<Scalar>(tanh(x.a), x.b*sech(x.a)*sech(x.a));
}
template<typename Scalar> DualVector<Scalar> sech(const DualVector<Scalar> & x) {
return DualVector<Scalar>(sech(x.a), x.b*(-sech(x.a)*tanh(x.a)));
}
template<typename Scalar> DualVector<Scalar> coth(const DualVector<Scalar> & x) {
return DualVector<Scalar>(coth(x.a), x.b*(-csch(x.a)*csch(x.a)));
}
template<typename Scalar> DualVector<Scalar> csch(const DualVector<Scalar> & x) {
return DualVector<Scalar>(csch(x.a), x.b*(-coth(x.a)*csch(x.a)));
}
// Inverse hyperbolic trigonometric functions
template<typename Scalar> DualVector<Scalar> acosh(const DualVector<Scalar> & x) {
return DualVector<Scalar>(acosh(x.a), x.b*(Scalar(1.)/sqrt((x.a*x.a)-Scalar(1.))));
}
template<typename Scalar> DualVector<Scalar> asinh(const DualVector<Scalar> & x) {
return DualVector<Scalar>(asinh(x.a), x.b*(Scalar(1.)/sqrt((x.a*x.a)+Scalar(1.))));
}
template<typename Scalar> DualVector<Scalar> atanh(const DualVector<Scalar> & x) {
return DualVector<Scalar>(atanh(x.a), x.b*(Scalar(1.)/(Scalar(1.)-(x.a*x.a))));
}
template<typename Scalar> DualVector<Scalar> asech(const DualVector<Scalar> & x) {
return DualVector<Scalar>(asech(x.a), x.b*(Scalar(-1.)/(sqrt(Scalar(1.)/(x.a*x.a)-Scalar(1.))*(x.a*x.a))));
}
template<typename Scalar> DualVector<Scalar> acoth(const DualVector<Scalar> & x) {
return DualVector<Scalar>(acoth(x.a), x.b*(-Scalar(1.)/((x.a*x.a)-Scalar(1.))));
}
template<typename Scalar> DualVector<Scalar> acsch(const DualVector<Scalar> & x) {
return DualVector<Scalar>(acsch(x.a), x.b*(-Scalar(1.)/(sqrt(Scalar(1.)/(x.a*x.a)+Scalar(1.))*(x.a*x.a))));
}
// Exponential functions
template<typename Scalar> DualVector<Scalar> exp(const DualVector<Scalar> & x) {
return DualVector<Scalar>(exp(x.a), x.b*exp(x.a));
}
template<typename Scalar> DualVector<Scalar> log(const DualVector<Scalar> & x) {
return DualVector<Scalar>(log(x.a), x.b/x.a);
}
template<typename Scalar> DualVector<Scalar> exp10(const DualVector<Scalar> & x) {
return DualVector<Scalar>(exp10(x.a), x.b*(log(Scalar(10.))*exp10(x.a)));
}
template<typename Scalar> DualVector<Scalar> log10(const DualVector<Scalar> & x) {
return DualVector<Scalar>(log10(x.a), x.b/(log(Scalar(10.))*x.a));
}
template<typename Scalar> DualVector<Scalar> exp2(const DualVector<Scalar> & x) {
return DualVector<Scalar>(exp2(x.a), x.b*(log(Scalar(2.))*exp2(x.a)));
}
template<typename Scalar> DualVector<Scalar> log2(const DualVector<Scalar> & x) {
return DualVector<Scalar>(log2(x.a), x.b/(log(Scalar(2.))*x.a));
}
template<typename Scalar> DualVector<Scalar> pow(const DualVector<Scalar> & x, const DualVector<Scalar> & n) {
return exp(n*log(x));
}
// Other functions
template<typename Scalar> DualVector<Scalar> sqrt(const DualVector<Scalar> & x) {
return DualVector<Scalar>(sqrt(x.a), x.b/(Scalar(2.)*sqrt(x.a)));
}
template<typename Scalar> DualVector<Scalar> sign(const DualVector<Scalar> & x) {
return DualVector<Scalar>(sign(x.a), DualVector<Scalar>::DualVector::__create_VectorT_zeros());
}
template<typename Scalar> DualVector<Scalar> abs(const DualVector<Scalar> & x) {
return DualVector<Scalar>(abs(x.a), x.b*sign(x.a));
}
template<typename Scalar> DualVector<Scalar> fabs(const DualVector<Scalar> & x) {
return DualVector<Scalar>(fabs(x.a), x.b*sign(x.a));
}
template<typename Scalar> DualVector<Scalar> heaviside(const DualVector<Scalar> & x) {
return DualVector<Scalar>(heaviside(x.a), DualVector<Scalar>::DualVector::__create_VectorT_zeros());
}
template<typename Scalar> DualVector<Scalar> floor(const DualVector<Scalar> & x) {
return DualVector<Scalar>(floor(x.a), DualVector<Scalar>::DualVector::__create_VectorT_zeros());
}
template<typename Scalar> DualVector<Scalar> ceil(const DualVector<Scalar> & x) {
return DualVector<Scalar>(ceil(x.a), DualVector<Scalar>::DualVector::__create_VectorT_zeros());
}
template<typename Scalar> DualVector<Scalar> round(const DualVector<Scalar> & x) {
return DualVector<Scalar>(round(x.a), DualVector<Scalar>::DualVector::__create_VectorT_zeros());
}
template<typename Scalar> std::ostream & operator<<(std::ostream & s, const DualVector<Scalar> & x)
{
return (s << x.a);
}
#endif

10
examples/firstAndSecondDerivative.cpp
View file

@ -1,4 +1,4 @@
#include "../AutomaticDifferentiationVector.hpp"
#include "../AutomaticDifferentiation.hpp"
#include <iostream>
#include <iomanip>
@ -39,7 +39,7 @@ int main()
// 1st derivative forward
{
using Fd = DualVector<double,1>;
using Fd = Dual<double>;
Fd x = xdbl;
x.diff(0);
Fd y = f(x);
@ -51,10 +51,10 @@ int main()
// 2nd derivative forward
/*
{
using Fdd = DualVector<DualVector<double,2>,2>;
using Fdd = Dual<Dual<double>>;
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;

411
old/AutomaticDifferentiation_scalar.hpp Executable file
View file

@ -0,0 +1,411 @@
#ifndef DEF_AUTOMATIC_DIFFERENTIATION
#define DEF_AUTOMATIC_DIFFERENTIATION
#include <cmath>
#include <ostream>
/// 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<typename Scalar>
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<typename A, typename B, int N>
Dual<B, N> operator+(A const& v, Dual<B, N> const& x) {
return (Dual<B, N>(v) + x);
}
template<typename A, typename B, int N>
Dual<B, N> operator-(A const& v, Dual<B, N> const& x) {
return (Dual<B, N>(v) - x);
}
template<typename A, typename B, int N>
Dual<B, N> operator*(A const& v, Dual<B, N> const& x) {
return (Dual<B, N>(v) * x);
}
template<typename A, typename B, int N>
Dual<B, N> operator/(A const& v, Dual<B, N> const& x) {
return (Dual<B, N>(v) / x);
}
// Basic mathematical functions for Scalar numbers
// Trigonometric functions
template<typename Scalar> Scalar sec(const Scalar & x) {
return Scalar(1.)/cos(x);
}
template<typename Scalar> Scalar cot(const Scalar & x) {
return cos(x)/sin(x);
}
template<typename Scalar> Scalar csc(const Scalar & x) {
return Scalar(1.)/sin(x);
}
// Inverse trigonometric functions
template<typename Scalar> Scalar asec(const Scalar & x) {
return acos(Scalar(1.)/x);
}
template<typename Scalar> Scalar acot(const Scalar & x) {
return atan(Scalar(1.)/x);
}
template<typename Scalar> Scalar acsc(const Scalar & x) {
return asin(Scalar(1.)/x);
}
// Hyperbolic trigonometric functions
template<typename Scalar> Scalar sech(const Scalar & x) {
return Scalar(1.)/cosh(x);
}
template<typename Scalar> Scalar coth(const Scalar & x) {
return cosh(x)/sinh(x);
}
template<typename Scalar> Scalar csch(const Scalar & x) {
return Scalar(1.)/sinh(x);
}
// Inverse hyperbolic trigonometric functions
template<typename Scalar> Scalar asech(const Scalar & x) {
return log((Scalar(1.) + sqrt(Scalar(1.) - x*x))/x);
}
template<typename Scalar> Scalar acoth(const Scalar & x) {
return Scalar(0.5)*log((x + Scalar(1.))/(x - Scalar(1.)));
}
template<typename Scalar> 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<typename Scalar> Scalar exp10(const Scalar & x) {
return exp(x*log(Scalar(10.)));
}
template<typename Scalar> Scalar sign(const Scalar & x) {
return (x >= Scalar(0.)) ? ((x > Scalar(0.)) ? Scalar(1.) : Scalar(0.)) : Scalar(-1.);
}
template<typename Scalar> Scalar heaviside(const Scalar & x) {
return Scalar(x >= Scalar(0.));
}
template<typename Scalar> 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<typename Scalar> Dual<Scalar> cos(const Dual<Scalar> & x) {
return Dual<Scalar>(cos(x.a), -x.b*sin(x.a));
}
template<typename Scalar> Dual<Scalar> sin(const Dual<Scalar> & x) {
return Dual<Scalar>(sin(x.a), x.b*cos(x.a));
}
template<typename Scalar> Dual<Scalar> tan(const Dual<Scalar> & x) {
return Dual<Scalar>(tan(x.a), x.b*sec(x.a)*sec(x.a));
}
template<typename Scalar> Dual<Scalar> sec(const Dual<Scalar> & x) {
return Dual<Scalar>(sec(x.a), x.b*sec(x.a)*tan(x.a));
}
template<typename Scalar> Dual<Scalar> cot(const Dual<Scalar> & x) {
return Dual<Scalar>(cot(x.a), x.b*(-csc(x.a)*csc(x.a)));
}
template<typename Scalar> Dual<Scalar> csc(const Dual<Scalar> & x) {
return Dual<Scalar>(csc(x.a), x.b*(-cot(x.a)*csc(x.a)));
}
// Inverse trigonometric functions
template<typename Scalar> Dual<Scalar> acos(const Dual<Scalar> & x) {
return Dual<Scalar>(acos(x.a), x.b*(-Scalar(1.)/sqrt(Scalar(1.)-x.a*x.a)));
}
template<typename Scalar> Dual<Scalar> asin(const Dual<Scalar> & x) {
return Dual<Scalar>(asin(x.a), x.b*(Scalar(1.)/sqrt(Scalar(1.)-x.a*x.a)));
}
template<typename Scalar> Dual<Scalar> atan(const Dual<Scalar> & x) {
return Dual<Scalar>(atan(x.a), x.b*(Scalar(1.)/(x.a*x.a + Scalar(1.))));
}
template<typename Scalar> Dual<Scalar> asec(const Dual<Scalar> & x) {
return Dual<Scalar>(asec(x.a), x.b*(Scalar(1.)/(sqrt(Scalar(1.)-Scalar(1.)/(x.a*x.a))*(x.a*x.a))));
}
template<typename Scalar> Dual<Scalar> acot(const Dual<Scalar> & x) {
return Dual<Scalar>(acot(x.a), x.b*(-Scalar(1.)/((x.a*x.a)+Scalar(1.))));
}
template<typename Scalar> Dual<Scalar> acsc(const Dual<Scalar> & x) {
return Dual<Scalar>(acsc(x.a), x.b*(-Scalar(1.)/(sqrt(Scalar(1.)-Scalar(1.)/(x.a*x.a))*(x.a*x.a))));
}
// Hyperbolic trigonometric functions
template<typename Scalar> Dual<Scalar> cosh(const Dual<Scalar> & x) {
return Dual<Scalar>(cosh(x.a), x.b*sinh(x.a));
}
template<typename Scalar> Dual<Scalar> sinh(const Dual<Scalar> & x) {
return Dual<Scalar>(sinh(x.a), x.b*cosh(x.a));
}
template<typename Scalar> Dual<Scalar> tanh(const Dual<Scalar> & x) {
return Dual<Scalar>(tanh(x.a), x.b*sech(x.a)*sech(x.a));
}
template<typename Scalar> Dual<Scalar> sech(const Dual<Scalar> & x) {
return Dual<Scalar>(sech(x.a), x.b*(-sech(x.a)*tanh(x.a)));
}
template<typename Scalar> Dual<Scalar> coth(const Dual<Scalar> & x) {
return Dual<Scalar>(coth(x.a), x.b*(-csch(x.a)*csch(x.a)));
}
template<typename Scalar> Dual<Scalar> csch(const Dual<Scalar> & x) {
return Dual<Scalar>(csch(x.a), x.b*(-coth(x.a)*csch(x.a)));
}
// Inverse hyperbolic trigonometric functions
template<typename Scalar> Dual<Scalar> acosh(const Dual<Scalar> & x) {
return Dual<Scalar>(acosh(x.a), x.b*(Scalar(1.)/sqrt((x.a*x.a)-Scalar(1.))));
}
template<typename Scalar> Dual<Scalar> asinh(const Dual<Scalar> & x) {
return Dual<Scalar>(asinh(x.a), x.b*(Scalar(1.)/sqrt((x.a*x.a)+Scalar(1.))));
}
template<typename Scalar> Dual<Scalar> atanh(const Dual<Scalar> & x) {
return Dual<Scalar>(atanh(x.a), x.b*(Scalar(1.)/(Scalar(1.)-(x.a*x.a))));
}
template<typename Scalar> Dual<Scalar> asech(const Dual<Scalar> & x) {
return Dual<Scalar>(asech(x.a), x.b*(Scalar(-1.)/(sqrt(Scalar(1.)/(x.a*x.a)-Scalar(1.))*(x.a*x.a))));
}
template<typename Scalar> Dual<Scalar> acoth(const Dual<Scalar> & x) {
return Dual<Scalar>(acoth(x.a), x.b*(-Scalar(1.)/((x.a*x.a)-Scalar(1.))));
}
template<typename Scalar> Dual<Scalar> acsch(const Dual<Scalar> & x) {
return Dual<Scalar>(acsch(x.a), x.b*(-Scalar(1.)/(sqrt(Scalar(1.)/(x.a*x.a)+Scalar(1.))*(x.a*x.a))));
}
// Exponential functions
template<typename Scalar> Dual<Scalar> exp(const Dual<Scalar> & x) {
return Dual<Scalar>(exp(x.a), x.b*exp(x.a));
}
template<typename Scalar> Dual<Scalar> log(const Dual<Scalar> & x) {
return Dual<Scalar>(log(x.a), x.b/x.a);
}
template<typename Scalar> Dual<Scalar> exp10(const Dual<Scalar> & x) {
return Dual<Scalar>(exp10(x.a), x.b*(log(Scalar(10.))*exp10(x.a)));
}
template<typename Scalar> Dual<Scalar> log10(const Dual<Scalar> & x) {
return Dual<Scalar>(log10(x.a), x.b/(log(Scalar(10.))*x.a));
}
template<typename Scalar> Dual<Scalar> exp2(const Dual<Scalar> & x) {
return Dual<Scalar>(exp2(x.a), x.b*(log(Scalar(2.))*exp2(x.a)));
}
template<typename Scalar> Dual<Scalar> log2(const Dual<Scalar> & x) {
return Dual<Scalar>(log2(x.a), x.b/(log(Scalar(2.))*x.a));
}
template<typename Scalar> Dual<Scalar> pow(const Dual<Scalar> & x, const Dual<Scalar> & n) {
return exp(n*log(x));
}
// Other functions
template<typename Scalar> Dual<Scalar> sqrt(const Dual<Scalar> & x) {
return Dual<Scalar>(sqrt(x.a), x.b/(Scalar(2.)*sqrt(x.a)));
}
template<typename Scalar> Dual<Scalar> sign(const Dual<Scalar> & x) {
return Dual<Scalar>(sign(x.a), Scalar(0.));
}
template<typename Scalar> Dual<Scalar> abs(const Dual<Scalar> & x) {
return Dual<Scalar>(abs(x.a), x.b*sign(x.a));
}
template<typename Scalar> Dual<Scalar> fabs(const Dual<Scalar> & x) {
return Dual<Scalar>(fabs(x.a), x.b*sign(x.a));
}
template<typename Scalar> Dual<Scalar> heaviside(const Dual<Scalar> & x) {
return Dual<Scalar>(heaviside(x.a), Scalar(0.));
}
template<typename Scalar> Dual<Scalar> floor(const Dual<Scalar> & x) {
return Dual<Scalar>(floor(x.a), Scalar(0.));
}
template<typename Scalar> Dual<Scalar> ceil(const Dual<Scalar> & x) {
return Dual<Scalar>(ceil(x.a), Scalar(0.));
}
template<typename Scalar> Dual<Scalar> round(const Dual<Scalar> & x) {
return Dual<Scalar>(round(x.a), Scalar(0.));
}
template<typename Scalar> std::ostream & operator<<(std::ostream & s, const Dual<Scalar> & x)
{
return (s << x.a);
}
#endif

0
test_AutomaticDifferentiation.cpp → old/test_AutomaticDifferentiation.cpp
View file

0
test_AutomaticDifferentiation_manual.cpp → old/test_AutomaticDifferentiation_manual.cpp
View file

0
test_AutomaticDifferentiation_vector.cpp → old/test_AutomaticDifferentiation_vector.cpp
View file

34
test_vector_version_additions_manual.cpp → test_small_manual.cpp
View file

@ -1,4 +1,4 @@
#include "AutomaticDifferentiationVector.hpp"
#include "AutomaticDifferentiation.hpp"
#include <iostream>
#include <iomanip>
@ -53,13 +53,13 @@ int main()
{
// 1 var function
{
DualVector<double> x = 1.54;
Dual<double> x = 1.54;
PRINT_DUAL(x);
x.diff(0);
PRINT_DUAL(x);
DualVector<double> y = f1(x);
Dual<double> y = f1(x);
PRINT_VAR(x);
PRINT_DUAL(y);
@ -71,11 +71,11 @@ int main()
// 3 var function
{
DualVector<double> x(2.), y(3.), z(-1.5);
Dual<double> x(2.), y(3.), z(-1.5);
x.diff(0, 3);
y.diff(1, 3);
z.diff(2, 3);
DualVector<double> res = f3(x, y, z);
Dual<double> 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<double> x;// by default, size(b)=1
Dual<double> 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<double>::D(0).b.size() == 1);
assert(DualVector<double>::D(0).b[0] == 1.);
assert(Dual<double>::D(0).b.size() == 1);
assert(Dual<double>::D(0).b[0] == 1.);
assert(DualVector<double>::D(0, 3).b[0] == 1.);
assert(DualVector<double>::D(0, 3).b[1] == 0.);
assert(DualVector<double>::D(0, 3).b[2] == 0.);
assert(Dual<double>::D(0, 3).b[0] == 1.);
assert(Dual<double>::D(0, 3).b[1] == 0.);
assert(Dual<double>::D(0, 3).b[2] == 0.);
assert(DualVector<double>::D(1, 3).b[0] == 0.);
assert(DualVector<double>::D(1, 3).b[1] == 1.);
assert(DualVector<double>::D(1, 3).b[2] == 0.);
assert(Dual<double>::D(1, 3).b[0] == 0.);
assert(Dual<double>::D(1, 3).b[1] == 1.);
assert(Dual<double>::D(1, 3).b[2] == 0.);
assert(DualVector<double>::D(2, 3).b[0] == 0.);
assert(DualVector<double>::D(2, 3).b[1] == 0.);
assert(DualVector<double>::D(2, 3).b[2] == 1.);
assert(Dual<double>::D(2, 3).b[0] == 0.);
assert(Dual<double>::D(2, 3).b[1] == 0.);
assert(Dual<double>::D(2, 3).b[2] == 1.);
}
return 0;