multidual.h

Go to the documentation of this file.

#ifndef THEORETICA_MULTIDUAL_H
#define THEORETICA_MULTIDUAL_H
 
 
 
#ifndef THEORETICA_NO_PRINT
#include <sstream>
#include <ostream>
#endif
 
#include "../algebra/vec.h"
#include "../algebra/mat.h"
 
 
namespace theoretica {
 
    template<unsigned int N = 0>
    class multidual {
        
        public:
            
            real a;
 
            vec<real, N> v;
 
            // The template argument of the vector type used
            static constexpr size_t vector_argument = N;
 
 
            multidual() : a(0) {}
 
 
            multidual(real r, vec<real, N> u) : a(r), v(u) {}
 
 
            multidual(real r) : a(r), v(vec<real, N>()) {}
 
            ~multidual() = default;
 
 
            inline multidual& operator=(real x) {
                a = x;
                v = vec<real, N>();
                return *this;
            }
 
 
            inline real Re() const {
                return a;
            }
 
 
            inline real& Re() {
                return a;
            }
 
 
            inline vec<real, N> Dual() const {
                return v;
            }
 
 
            inline vec<real, N>& Dual() {
                return v;
            }
 
 
            inline multidual conjugate() const {
                return multidual(a, -v);
            }
 
 
            inline multidual inverse() const {
 
                if(a == 0) {
                    TH_MATH_ERROR("multidual::inverse", 0, DIV_BY_ZERO);
                    return multidual(nan(), vec<real, N>(nan(), N));
                }
 
                return multidual(1.0 / a, v * (-1 / (a * a)));
            }
 
 
            inline multidual operator+() const {
                return multidual(a, v);
            }
 
 
            inline multidual operator+(const multidual& other) const {
                return multidual(a + other.a, v + other.v);
            }
 
 
            inline multidual operator+(real r) const {
                return multidual(a + r, v);
            }
 
 
            inline multidual operator-() const {
                return multidual(-a, -v);
            }
 
 
            inline multidual operator-(const multidual& other) const {
                return multidual(a - other.a, v - other.v);
            }
 
 
            inline multidual operator-(real r) const {
                return multidual(a - r, v);
            }
 
 
            inline multidual operator*(const multidual& other) const {
                return multidual(a * other.a, other.v * a + v * other.a);
            }
 
 
            inline multidual operator*(real r) const {
                return multidual(a * r, v * r);
            }
 
 
            inline multidual operator/(const multidual& other) const {
 
                if(a == 0) {
                    TH_MATH_ERROR("multidual::operator/", 0, DIV_BY_ZERO);
                    return multidual(nan(), vec<real, N>(nan(), N));
                }
 
                return multidual(a / other.a,
                    (v * other.a - other.v * a) / (other.a * other.a));
            }
 
            inline multidual operator/(real r) const {
                return multidual(a / r, v / r);
            }
 
 
            inline multidual& operator+=(const multidual& other) {
 
                a += other.a;
                v += other.v;
                return *this;
            }
 
 
            inline multidual& operator+=(real r) {
 
                a += r;
                return *this;
            }
 
            inline multidual& operator-=(const multidual& other) {
 
                a -= other.a;
                v -= other.v;
                return *this;
            }
 
            inline multidual& operator-=(real r) {
 
                a -= r;
                return *this;
            }
 
            inline multidual& operator*=(const multidual& other) {
 
                a = (a * other.a);
                v = (other.v * a) + (v * other.a);
                return *this;
            }
 
            inline multidual& operator*=(real r) {
 
                a *= r;
                v *= r;
                return *this;
            }
 
            inline multidual& operator/=(const multidual& other) {
 
                a = (a / other.a);
                v = (v * other.a - other.v * a) / (other.a * other.a);
                return *this;
            }
 
            inline multidual& operator/=(real r) {
 
                if(r == 0) {
                    TH_MATH_ERROR("multidual::operator/=", 0, DIV_BY_ZERO);
                    a = nan();
                    v = vec<real, N>(nan(), N);
                    return *this;
                }
 
                a /= r;
                v /= r;
 
                return *this;
            }
 
 
            inline bool operator==(const multidual& other) {
                return (a == other.a) && (v == other.v);
            }
 
 
            inline static vec<multidual<N>, N> make_argument(
                const vec<real, N>& x) {
 
                vec<multidual<N>, N> arg;
                arg.resize(x.size());
 
                for (unsigned int i = 0; i < x.size(); ++i)
                    arg[i] = multidual<N>(
                        x[i], vec<real, N>::euclidean_base(i, x.size())
                    );
 
                return arg;
            }
 
 
            inline static vec<real, N> extract_real(
                const vec<multidual<N>, N>& v) {
 
                vec<real, N> x;
                x.resize(v.size());
 
                for (unsigned int i = 0; i < N; ++i)
                    x[i] = v.get(i).Re();
 
                return x;
            }
 
 
            inline static mat<real, N, N> extract_dual(
                const vec<multidual<N>, N>& v) {
 
                mat<real, N, N> J;
                J.resize(v.size(), v.size());
 
                for (unsigned int i = 0; i < N; ++i)
                    for (unsigned int j = 0; j < N; ++j)
                        J.at(j, i) = v.get(j).Dual().get(i);
                
                return J;
            }
 
 
            inline static void extract(
                const vec<multidual<N>, N>& v,
                vec<real, N>& x,
                mat<real, N, N>& J) {
 
                for (unsigned int i = 0; i < N; ++i) {
                    
                    for (unsigned int j = 0; j < N; ++j)
                        J.at(j, i) = v.get(j).Dual().get(i);
                    
                    x[i] = v.get(i).Re();
                }
            }
 
 
            unsigned int size() const {
                return v.size();
            }
 
 
            inline void resize(unsigned int size) {
                v.resize(size);
            }
 
 
            // Friend operators to enable equations of the form
            // (real) op. (multidual)
            
            inline friend multidual operator+(real a, const multidual& d) {
                return d + a;
            }
 
            inline friend multidual operator-(real a, const multidual& d) {
                return -d + a;
            }
 
            inline friend multidual operator*(real a, const multidual& d) {
                return d * a;
            }
 
            inline friend multidual operator/(real a, const multidual& d) {
                return a * d.inverse();
            }
 
 
#ifndef THEORETICA_NO_PRINT
 
            inline std::string to_string(const std::string& epsilon = "e") const {
 
                std::stringstream res;
                res << a << " + " << v << epsilon;
 
                return res.str();
            }
 
 
            inline operator std::string() {
                return to_string();
            }
 
 
            inline friend std::ostream& operator<<(std::ostream& out, const multidual& obj) {
                return out << obj.to_string();
            }
 
#endif
        
    };
 
}
 
#endif