Fixed vector version to look more like FADBAD in terms of interface
This commit is contained in:
Jérôme
parent
bd46adc942
commit
864e633552
14 changed files with 220 additions and 113 deletions
4
.gitignore
vendored
Normal file → Executable file
4
.gitignore
vendored
Normal file → Executable file
|
|
@ -1,6 +1,10 @@
|
|||
# executables
|
||||
run
|
||||
|
||||
# catch2 generated files
|
||||
*.gcda
|
||||
*.gcno
|
||||
|
||||
# Prerequisites
|
||||
*.d
|
||||
|
||||
|
|
|
|||
0
AutomaticDifferentiation.hpp
Normal file → Executable file
0
AutomaticDifferentiation.hpp
Normal file → Executable file
280
AutomaticDifferentiationVector.hpp
Normal file → Executable file
280
AutomaticDifferentiationVector.hpp
Normal file → Executable file
|
|
@ -15,63 +15,87 @@ std::ostream & operator<<(std::ostream & out, std::valarray<T> 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<typename Scalar, int N>
|
||||
template<typename Scalar>
|
||||
class DualVector
|
||||
{
|
||||
public:
|
||||
using VectorT = std::valarray<Scalar>;
|
||||
|
||||
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<typename A, typename B, int N>
|
||||
DualVector<B, N> operator+(A const& v, DualVector<B, N> const& x) {
|
||||
return (DualVector<B, N>(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, int N>
|
||||
DualVector<B, N> operator-(A const& v, DualVector<B, N> const& x) {
|
||||
return (DualVector<B, N>(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, int N>
|
||||
DualVector<B, N> operator*(A const& v, DualVector<B, N> const& x) {
|
||||
return (DualVector<B, N>(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, int N>
|
||||
DualVector<B, N> operator/(A const& v, DualVector<B, N> const& x) {
|
||||
return (DualVector<B, N>(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
|
||||
|
|
@ -211,6 +278,7 @@ 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.)));
|
||||
}
|
||||
|
|
@ -231,167 +299,167 @@ template<typename Scalar> Scalar abs(const Scalar & x) {
|
|||
// f(a + b*d) = f(a) + b*f'(a)*d
|
||||
|
||||
// Trigonometric functions
|
||||
template<typename Scalar, int N> DualVector<Scalar, N> cos(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(cos(x.a), -x.b*sin(x.a));
|
||||
template<typename Scalar> DualVector<Scalar> cos(const DualVector<Scalar> & x) {
|
||||
return DualVector<Scalar>(cos(x.a), -x.b*sin(x.a));
|
||||
}
|
||||
|
||||
template<typename Scalar, int N> DualVector<Scalar, N> sin(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(sin(x.a), x.b*cos(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, int N> DualVector<Scalar, N> tan(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(tan(x.a), x.b*sec(x.a)*sec(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, int N> DualVector<Scalar, N> sec(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(sec(x.a), x.b*sec(x.a)*tan(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, int N> DualVector<Scalar, N> cot(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(cot(x.a), x.b*(-csc(x.a)*csc(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, int N> DualVector<Scalar, N> csc(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(csc(x.a), x.b*(-cot(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, int N> DualVector<Scalar, N> acos(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(acos(x.a), x.b*(-Scalar(1.)/sqrt(Scalar(1.)-x.a*x.a)));
|
||||
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, int N> DualVector<Scalar, N> asin(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(asin(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, int N> DualVector<Scalar, N> atan(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(atan(x.a), x.b*(Scalar(1.)/(x.a*x.a + Scalar(1.))));
|
||||
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, int N> DualVector<Scalar, N> asec(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(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> 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, int N> DualVector<Scalar, N> acot(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(acot(x.a), x.b*(-Scalar(1.)/((x.a*x.a)+Scalar(1.))));
|
||||
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, int N> DualVector<Scalar, N> acsc(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(acsc(x.a), x.b*(-Scalar(1.)/(sqrt(Scalar(1.)-Scalar(1.)/(x.a*x.a))*(x.a*x.a))));
|
||||
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, int N> DualVector<Scalar, N> cosh(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(cosh(x.a), x.b*sinh(x.a));
|
||||
template<typename Scalar> DualVector<Scalar> cosh(const DualVector<Scalar> & x) {
|
||||
return DualVector<Scalar>(cosh(x.a), x.b*sinh(x.a));
|
||||
}
|
||||
|
||||
template<typename Scalar, int N> DualVector<Scalar, N> sinh(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(sinh(x.a), x.b*cosh(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, int N> DualVector<Scalar, N> tanh(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(tanh(x.a), x.b*sech(x.a)*sech(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, int N> DualVector<Scalar, N> sech(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(sech(x.a), x.b*(-sech(x.a)*tanh(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, int N> DualVector<Scalar, N> coth(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(coth(x.a), x.b*(-csch(x.a)*csch(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, int N> DualVector<Scalar, N> csch(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(csch(x.a), x.b*(-coth(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, int N> DualVector<Scalar, N> acosh(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(acosh(x.a), x.b*(Scalar(1.)/sqrt((x.a*x.a)-Scalar(1.))));
|
||||
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, int N> DualVector<Scalar, N> asinh(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(asinh(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, int N> DualVector<Scalar, N> atanh(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(atanh(x.a), x.b*(Scalar(1.)/(Scalar(1.)-(x.a*x.a))));
|
||||
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, int N> DualVector<Scalar, N> asech(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(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> 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, int N> DualVector<Scalar, N> acoth(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(acoth(x.a), x.b*(-Scalar(1.)/((x.a*x.a)-Scalar(1.))));
|
||||
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, int N> DualVector<Scalar, N> acsch(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(acsch(x.a), x.b*(-Scalar(1.)/(sqrt(Scalar(1.)/(x.a*x.a)+Scalar(1.))*(x.a*x.a))));
|
||||
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, int N> DualVector<Scalar, N> exp(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(exp(x.a), x.b*exp(x.a));
|
||||
template<typename Scalar> DualVector<Scalar> exp(const DualVector<Scalar> & x) {
|
||||
return DualVector<Scalar>(exp(x.a), x.b*exp(x.a));
|
||||
}
|
||||
|
||||
template<typename Scalar, int N> DualVector<Scalar, N> log(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(log(x.a), x.b/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, int N> DualVector<Scalar, N> exp10(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(exp10(x.a), x.b*(log(Scalar(10.))*exp10(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, int N> DualVector<Scalar, N> log10(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(log10(x.a), x.b/(log(Scalar(10.))*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, int N> DualVector<Scalar, N> exp2(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(exp2(x.a), x.b*(log(Scalar(2.))*exp2(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, int N> DualVector<Scalar, N> log2(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(log2(x.a), x.b/(log(Scalar(2.))*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, int N> DualVector<Scalar, N> pow(const DualVector<Scalar, N> & x, const DualVector<Scalar, N> & n) {
|
||||
template<typename Scalar> DualVector<Scalar> pow(const DualVector<Scalar> & x, const DualVector<Scalar> & n) {
|
||||
return exp(n*log(x));
|
||||
}
|
||||
|
||||
// Other functions
|
||||
template<typename Scalar, int N> DualVector<Scalar, N> sqrt(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(sqrt(x.a), x.b/(Scalar(2.)*sqrt(x.a)));
|
||||
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, int N> DualVector<Scalar, N> sign(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(sign(x.a), DualVector<Scalar, N>::DualVector::__create_VectorT_zeros());
|
||||
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, int N> DualVector<Scalar, N> abs(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(abs(x.a), x.b*sign(x.a));
|
||||
template<typename Scalar> DualVector<Scalar> abs(const DualVector<Scalar> & x) {
|
||||
return DualVector<Scalar>(abs(x.a), x.b*sign(x.a));
|
||||
}
|
||||
template<typename Scalar, int N> DualVector<Scalar, N> fabs(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(fabs(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, int N> DualVector<Scalar, N> heaviside(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(heaviside(x.a), DualVector<Scalar, N>::DualVector::__create_VectorT_zeros());
|
||||
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, int N> DualVector<Scalar, N> floor(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(floor(x.a), DualVector<Scalar, N>::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, int N> DualVector<Scalar, N> ceil(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(ceil(x.a), DualVector<Scalar, N>::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, int N> DualVector<Scalar, N> round(const DualVector<Scalar, N> & x) {
|
||||
return DualVector<Scalar, N>(round(x.a), DualVector<Scalar, N>::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, int N> std::ostream & operator<<(std::ostream & s, const DualVector<Scalar, N> & x)
|
||||
template<typename Scalar> std::ostream & operator<<(std::ostream & s, const DualVector<Scalar> & x)
|
||||
{
|
||||
return (s << x.a);
|
||||
}
|
||||
|
|
|
|||
0
Makefile
Normal file → Executable file
0
Makefile
Normal file → Executable file
0
Makefile.linux
Normal file → Executable file
0
Makefile.linux
Normal file → Executable file
0
Makefile.windows
Normal file → Executable file
0
Makefile.windows
Normal file → Executable file
0
README.md
Normal file → Executable file
0
README.md
Normal file → Executable file
0
coverage_test.txt
Normal file → Executable file
0
coverage_test.txt
Normal file → Executable file
0
examples/firstAndSecondDerivative.cpp
Normal file → Executable file
0
examples/firstAndSecondDerivative.cpp
Normal file → Executable file
0
test_AutomaticDifferentiation.cpp
Normal file → Executable file
0
test_AutomaticDifferentiation.cpp
Normal file → Executable file
0
test_AutomaticDifferentiation_main.cpp
Normal file → Executable file
0
test_AutomaticDifferentiation_main.cpp
Normal file → Executable file
0
test_AutomaticDifferentiation_manual.cpp
Normal file → Executable file
0
test_AutomaticDifferentiation_manual.cpp
Normal file → Executable file
0
test_AutomaticDifferentiation_vector.cpp
Normal file → Executable file
0
test_AutomaticDifferentiation_vector.cpp
Normal file → Executable file
49
test_vector_version_additions_manual.cpp
Normal file → Executable file
49
test_vector_version_additions_manual.cpp
Normal file → Executable file
|
|
@ -53,13 +53,13 @@ int main()
|
|||
{
|
||||
// 1 var function
|
||||
{
|
||||
DualVector<double, 1> x = 1.54;
|
||||
DualVector<double> x = 1.54;
|
||||
|
||||
PRINT_DUAL(x);
|
||||
x.diff(0);
|
||||
PRINT_DUAL(x);
|
||||
|
||||
DualVector<double, 1> y = f1(x);
|
||||
DualVector<double> y = f1(x);
|
||||
|
||||
PRINT_VAR(x);
|
||||
PRINT_DUAL(y);
|
||||
|
|
@ -71,11 +71,11 @@ int main()
|
|||
|
||||
// 3 var function
|
||||
{
|
||||
DualVector<double,3> x(2.), y(3.), z(-1.5);
|
||||
x.diff(0);
|
||||
y.diff(1);
|
||||
z.diff(2);
|
||||
DualVector<double,3> res = f3(x, y, z);
|
||||
DualVector<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);
|
||||
|
||||
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<double> 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<double>::D(0).b.size() == 1);
|
||||
assert(DualVector<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(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(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.);
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue