Generalized Dual numbers

#include <audi/gdual.hpp>

template<typename Cf, typename Monomial = obake::d_packed_monomial<std::uint64_t, 8>>
class gdual

Generalized dual number class.

This class represents a generalized dual number, in a nutshell, a truncated multivariate Taylor polynomial. Using the multi-index notation, a generalized dual number may be written as:

\[T_f(\mathbf x) = \sum_{|\alpha| = 0}^m \frac{(\mathbf x-\mathbf a)^\alpha}{\alpha!}(\partial^\alpha f)(\mathbf a) \]

and thus depends on the order \(m\) as well as on the expansion point \(\mathbf a\). All arithmetic operators +,*,/,- are overloaded so that the Taylor expansion of arithmetic computations is obtained. A basic use case, where \(m = 2\), \(\mathbf a = [1.2, -0.1]\) and \(f = \frac{x1+x2}{x1-x2}\) is thus:

gdual<double> x1(1.2, "x1", 2);
gdual<double> x2(-0.1, "x2", 2);
std::cout << (x1+x2) / (x1-x2) << "\n";

resulting in the output:

0.118343*dx1+0.846154+1.42012*dx2+1.0924*dx2**2-0.0910332*dx1**2-1.00137*dx1*dx2

Integration and differentiation are also implemented so that the generalized dual computations are formally made in a differential algebra.

Note

The class can be instantiated with any type that is suitable to be a coefficient in a obake polynomial.

Low-level interface

friend class boost::serialization::access
inline decltype(m_p._get_s_table()) _container()

Get a mutable reference to the container of terms.

Returns:

a reference to the internal container of terms.

inline decltype(m_p._get_s_table()) _container() const

Get a const reference to the container of terms.

Returns:

a const reference to the internal container of terms.

inline const p_type &_poly() const

Get a const reference to the polynomial.

Returns:

a const reference to the internal polynomial.

inline p_type &_poly()

Get a mutable reference to the polynomial.

Returns:

a reference to the internal polynomial.

Comparison operators

inline friend bool operator==(const gdual &d1, const gdual &d2)

Overloaded Equality operator.

Compares the single polynomial coefficients of two audi::gdual objects and returns true if equal.

/note The truncation order of d1 and d2 may be different

Parameters:
Returns:

The result of the comparison

inline friend bool operator!=(const gdual &d1, const gdual &d2)

Overloaded Non equality operator.

Compares the single polynomial coefficients of two audi::gdual objects and returns false if equal.

Parameters:
Returns:

The result of the comparison

inline friend bool operator<(const gdual &d1, const gdual &d2)

Overloaded less than operator.

Compares the constant coefficients of two gduals

Parameters:
Returns:

The result of the comparison

inline friend bool operator>(const gdual &d1, const gdual &d2)

Overloaded more than operator.

Compares the constant coefficients of two gduals

Parameters:
Returns:

The result of the comparison

Algebraic operators

template<typename T, typename U>
inline friend gdual_if_enabled<T, U> operator+(const T &d1, const U &d2)

Overloaded addition operator.

Implements the sum operation between truncated Taylor polynomials.

Note

In order for this overload to be active (SFINAE rules), at least one of the arguments must be an audi::gdual, while the second argument may only be a double or int.

Parameters:
Returns:

the sum between d1 and d2

template<typename T, typename U>
inline friend gdual_if_enabled<T, U> operator-(const T &d1, const U &d2)

Overloaded difference operator.

Implements the difference operation between truncated Taylor polynomials.

Note

In order for this overload to be active (SFINAE rules), at least one of the arguments must be an audi::gdual, while the second argument may only be a double or int.

Parameters:
Returns:

the difference between d1 and d2

template<typename T, typename U>
inline friend gdual_if_enabled<T, U> operator*(const T &d1, const U &d2)

Overloaded multiplication operator.

Implements the multiplication operation between truncated Taylor polynomials.

Note

In order for this overload to be active (SFINAE rules), at least one of the arguments must be an audi::gdual, while the second argument may only be a double or int.

Note

The truncated polynomial multiplication operator is at the very heart of AuDi and its details / performances are those of the obake multiplication algorithm which is, essentially, used.

Parameters:
Returns:

the (truncated) multiplication between d1 and d2

template<typename T, typename U>
inline friend gdual_if_enabled<T, U> operator/(const T &d1, const U &d2)

Overloaded division operator.

Implements the division operation between truncated Taylor polynomials. Essentially, defined (in case of two audi::gdual) by a multiplication and the reciprocal rule:

\[T_g = \frac 1{T_f} = \frac 1{f_0} \left(1 +\sum_{k=1}^m (-1)^k (\hat f / f_0)^k\right) \]

where \(T_f = f_0 + \hat f\).

Note

In order for this overload to be active (SFINAE rules), at least one of the arguments must be an audi::gdual, while the second argument may only be a double or int.

Parameters:
Returns:

an audi:gdual containing the (truncated) multiplication between d1 and d2

template<typename T>
inline decltype(*this = *this + d1) operator+=(const T &d1)

Add and assignment operator.

template<typename T>
inline decltype(*this = *this - d1) operator-=(const T &d1)

Subtract and assignment operator.

template<typename T>
inline decltype(*this = *this * d1) operator*=(const T &d1)

Multiply and assignment operator.

template<typename T>
inline decltype(*this = *this / d1) operator/=(const T &d1)

Divide and assignment operator.

inline gdual operator-() const

Negate operator.

inline gdual operator+() const

Identity operator.

Public Types

using cf_type = Cf
using key_type = Monomial
using p_type = obake::polynomial<key_type, cf_type>

Public Functions

gdual(const gdual&) = default

Defaulted copy constructor.

gdual(gdual&&) = default

Defaulted move constructor.

inline gdual()

Default constuctor.

inline ~gdual()

Destructor (contains a sanity check)

template<typename T, generic_ctor_enabler<T> = 0>
inline explicit gdual(const T &value)
template<typename T>
inline explicit gdual(const std::initializer_list<T> &value)
template<typename T, generic_ctor_enabler<T> = 0>
inline explicit gdual(const T &value, const std::string &symbol, unsigned order)
template<typename T>
inline explicit gdual(const std::initializer_list<T> &value, const std::string &symbol, unsigned order)
gdual &operator=(const gdual&) = default

Defaulted assignment operator.

gdual &operator=(gdual&&) = default

Defaulted assignment operator.

inline auto get_symbol_set_size() const

Gets the symbol set size.

Returns the size of the symbol set.

Returns:

the size of the symbol set.

inline std::vector<std::string> get_symbol_set() const

Returns the symbol set.

Returns the symbol set (not the differentials)

Returns:

the symbol set

inline void extend_symbol_set(const std::vector<std::string> &sym_vars)

Extends the symbol set.

Adds some symbolic variables to the current polynomial This is useful in situations where some variable \( dx\) does not appear in the final polynomial but we still want to treat the Taylor expansion as a function of \( dx\) too (for example when extracting the relative coefficient)

Parameters:

sym_vars – list of symbolic names. It must contain all symbolic names of the current polynomial. It may contain more. All symbols must start with the letter “d”.

Throws:

std::invalid_argument –

  • if any symbol in sym_vars does not start with the letter “d”

template<enable_if_t<obake::is_differentiable<Cf>::value, int> = 0>
inline gdual integrate(const std::string &var_name)

Integration.

Performs the integration of the gdual with respect to the symbol.

Note

If the var_name differential is not in the symbol set, then it is added.

Note

Information may be lost as the truncation order is preserved.

Parameters:

var_name – Symbol name (cannot start with “d”).

Throws:

  • std::invalid_argument –

    • if var_name starts with the letter “d” (this avoid creating confusing names for symbol’s differentials)

  • unspecified – any exception thrown by:

    • obake::truncate_degree,

    • obake::integrate

template<enable_if_t<obake::is_differentiable<Cf>::value, int> = 0>
inline gdual partial(const std::string &var_name)

Partial derivative.

Performs the partial derivative of the gdual with respect to the symbol

Note

If the var_name differential is not in the symbol set, then it is added.

Parameters:

var_name – Symbol name (cannot start with “d”).

Throws:

  • std::invalid_argument –

    • if symbol starts with the letter “d” (this avoid creating confusing names for symbol’s differentials)

  • unspecified – any exception thrown by:

    • obake::diff,

template<typename T>
inline gdual subs(const std::string &sym, const T &val) const

Substitute symbol with value.

Substitute the symbol sym with the value val

Throws:

unspecified – any exception thrown by:

  • obake::subs,

inline gdual subs(const std::string &sym, const gdual &val) const

Substitute a symbol with a gdual.

Substitute the symbol sym with the gdual val The returned gdual has the same order, its symbol set will be expanded with the symbol set of val and then trimmed of the substituded symbol and truncated to the order of the gdual.

Parameters:

  • sym – Symbol to be substituded (i.e. “dx”)

  • val – The gdual to substitute sym with.

Throws:

unspecified – any exception thrown by:

  • obake::subs,

  • obake::trim

  • obake::truncate_degree

inline gdual trim(double epsilon) const

Trims small coefficients.

All coefficients smaller than a set magnitude epsilon are set to zero

Parameters:

epsilon – Tolerance used to trim coefficients.

Returns:

The new trimmed gdual.

inline gdual extract_terms(unsigned order) const

Extracts all terms of some order.

Extracts all the terms with the given order into a new gdual

Parameters:

order – Order of the terms to be extracted

Throws:

std::invalid_argument –

  • if the order is higher than the gdual order

Returns:

A new gdual containing the terms extracted, but preserving the order of the original gdual.

inline auto evaluate(const obake::symbol_map<double> &dict) const

Evaluates the Taylor polynomial.

Evaluates the Taylor polynomial using the values in dict for all the differentials (variables variations)

Throws:

unspecified – any exception thrown by:

  • obake::evaluate,

inline auto degree() const

Current degree.

Returns the current degree of the polynomial represented as an audi::gdual. This may be different from the truncation order \(m\) and, in particular will be smaller or equal.

Returns:

the current degree of the polynomial

inline unsigned int get_order() const

Getter for the truncation order.

Returns the truncation order of the underlying \(\mathcal P_{n,m}\) algebra.

Returns:

the truncation order

template<typename T>
inline Cf find_cf(const T &c) const

Finds the coefficient of a particular monomial.

Returns the coefficient of the monomial specified in the container c

This method will first construct a term with zero coefficient and key initialised from the begin/end iterators of c and the symbol set of this, and it will then try to locate the term inside this. If the term is found, its coefficient will be returned. Otherwise, a coefficient initialised from 0 will be returned.

Note

The container contains the exponents of the requested monomial. In a three variable polynomial with “x”, “y” and “z” as symbols, [1, 3, 2] would denote the term x y^3 z^2.

Note

Alphabetical order is used to order the symbol set and thus specify the coefficients.

Throws:

std::invalid_argument –

  • if the requested coefficient is beyond the truncation order

Returns:

the coefficient

template<typename T>
inline Cf find_cf(std::initializer_list<T> l) const

Finds the coefficient of a particular monomial.

Returns the coefficient of the monomial specified in the initializer list l

Note

This method is identical to the other overload with the same name, and it is provided for convenience.

Throws:

std::invalid_argument –

  • if the requested coefficient is beyond the truncation order

Returns:

the coefficient

template<typename T>
inline Cf get_derivative(const T &c) const

Gets the derivative value.

Returns the (mixed) derivative value of order specified by the container c

Note

The container contains the order requested. In a three variable polynomial with “x”, “y” and “z” as symbols, [1, 3, 2] would denote the sixth order derivative \( \frac{d^6}{dxdy^3dz^2}\).

Note

No computations are made at this points as all derivatives are already contained in the Taylor expansion

Throws:

std::invalid_argument –

  • if the requested coefficient is beyond the truncation order

Returns:

the value of the derivative

template<typename T>
inline Cf get_derivative(std::initializer_list<T> l) const

Gets the derivative value.

Returns the (mixed) derivative value of order specified by the initializer list l

Note

This method is identical to the other overload with the same name, and it is provided for convenience.

Note

No computations are made at this points as all derivatives are already contained in the Taylor expansion

Throws:

std::invalid_argument –

  • if the requested coefficient is beyond the truncation order

Returns:

the value of the derivative

inline auto get_derivative(const std::unordered_map<std::string, unsigned> &dict) const

Gets the derivative value.

Returns the (mixed) derivative value of order specified by the unordered_map dict

Note

To get the following derivative: \( \frac{d^6}{dxdy^3dz^2}\) the input should be {{“dx”, 1u},{“dy”,3u},{“dz”,2u}}

Note

The current implementation call internally the other templated implementations. When obake will implement the sparse monomial this will change and be more efficient.

Throws:

  • unspecified – all exceptions thrown by the templated version call.

  • std::invalid_argument – if one of the symbols is not found in the expression

Returns:

the value of the derivative

inline Cf constant_cf() const

Finds the constant coefficient.

Returns the coefficient of the of the constant part of the polynomial so that if \(T_{f} = f_0 + \hat f\), \(f_0\) is returned.

Returns:

the coefficient

inline bool is_zero(double tol) const

Determines if a gdual is zero within tolerance.

Returns:

true if all coefficients of the gdual are zero within a tolerance tol

inline std::string info() const

Info on the gdual.

Will assemble a string containing information on the gdual. It uses the obake streaming operator for the type obake::series. Refer to that documentation for further details

Friends

inline friend std::ostream &operator<<(std::ostream &os, const gdual &d)

Overloaded stream operator.

Will direct to stream a human-readable representation of the generalized dual number.

Note

The print order of the terms will be undefined.

Parameters:

  • os – [inout] target stream.

  • d – audi::gdual argument.

Returns:

reference to os.