Generalized Dual numbers
=========================

*#include <audi/gdual.hpp>*

.. cpp:class:: template <typename Cf> gdual

This class represents an element of the algebra :math:`\mathcal P_{n,m}`, (see :ref:`formal_definition`) defined
over the field :math:`\mathbf K` (the field is represented by the template argument **Cf**). We call elements of this algebra
generalized dual numbers as, among other things and when :math:`\mathbf K` is :math:`\mathbb R`, they generalize the dual numbers used for forward automatic differentiation.

Using the multi-index notation, a generalized dual number (for example over the field :math:`\mathbb R` represented by doubles, i.e. :code:`gdual<double>`) 
may be written as:

.. math::

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

.. note::

   A generalized dual number is defined by its truncation order :math:`m` as well as by its expansion point :math:`\mathbf a`. 
   All arithmetic operators +,*,/,- are overloaded so that computing with the type :cpp:class:`gdual` correspond to compute the
   Taylor expansion of the arithmetic computations.

A basic use case, where :math:`m = 2`, :math:`\mathbf a = [1.2, -0.1]` and :math:`f = \frac{x1+x2}{x1-x2}` is thus:

.. code-block:: c++

   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:

.. code-block:: c++

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

.. note::

  The class can be instantiated with any type that is suitable to be a coefficient in an obake polynomial (
  obake::is_cf_v<Cf> must be true). Classical examples would be double, float, std::complex<double>, and
  the audi::vectorized_double type. If obake::is_differentiable_v<Cf> is also true then derivation
  and integration are also availiable and the resulting algebra is a differential algebra.

------------------------------------------------------

Constructors and destructors
------------------------------------------------------

.. cpp:function:: template <typename T> explicit gdual(const T &c)

   .. note::  

      This constructor is only enabled if *T* is not a gdual and if an obake::polynomial can be constructed with coefficients in *T*.

   Constructs a gdual of order :math:`m = 0` with only a constant term :math:`c`. Essentially :math:`T_f = c` in :math:`\mathcal P_{0,0}`

   :param c: The constant term :math:`c`.

------------------------------------------------------

.. cpp:function:: template <typename T> explicit gdual(const T &c, const std::string &symbol, unsigned m)

   .. note::  

      This constructor is only enabled if *T* is not a gdual and if an obake::polynomial can be constructed with coefficients in *T*.

   Constructs a gdual of order :math:`m` representing a variable. Essentially :math:`T_f = c + dx` in :math:`\mathcal P_{1,m}`

   :param c: The constant term :math:`c`.
   :param symbol: The symbol that will represent the variable (e.g. :math:`x`).
   :param m: The truncation order of the algebra.

   :exception: std::invalid_argument if the symbol names starts already with "d" a differential. This avoids symbols like ddx in the obake::polynomial.

------------------------------------------------------

.. cpp:function:: gdual()

   Default constructor. Constructs a gdual of order :math:`m = 0` with only a zero constant term :math:`c`. Essentially :math:`T_f = 0` in :math:`\mathcal P_{0,0}`

------------------------------------------------------

.. cpp:function:: gdual(const gdual &) = default

   Defaulted copy constructor.

------------------------------------------------------

.. cpp:function:: gdual(gdual &&) = default

   Defaulted move constructor.

------------------------------------------------------

.. cpp:function:: ~gdual()

   Destructor. Performs a sanity check on the truncation order and degree of the gdual.

------------------------------------------------------


Methods
-------

Differential algebra operations
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

.. cpp:function:: template<> gdual integrate(const std::string &var_name)

   .. note::
      
      This template is only enabled if *Cf* satisfies obake::is_differentiable, which is
      the case for float, double, std::complex and vectorized_double types. 

   Performs the integration of the gdual with respect to *var_name*. If the *var_name* differential is not in the symbol set
   of the obake::polynomial it is added.
   Note that Information may be lost as the truncation order is preserved.

   :param var_name: Symbol name (cannot start with "d").

   :exception: std::invalid_argument if *var_name* starts with the letter "d".

------------------------------------------------------

.. cpp:function:: template<> gdual partial(const std::string &var_name)

   .. note::
      
      This template is only enabled if *Cf* satisfies obake::is_differentiable, which is
      the case for float, double, std::complex and vectorized_double types. 

   Performs the partial derivative of the gdual with respect to *var_name*. If the *var_name* differential is not in the symbol set
   of the obake::polynomial it is added.

   :param var_name: Symbol name (cannot start with "d").

   :exception: std::invalid_argument if *var_name* starts with the letter "d".

------------------------------------------------------

gdual manipulation
^^^^^^^^^^^^^^^^^^^

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

   Substitute the differential *sym* with *val*. The *Cf* type must be constructable from *val*.

   :param sym: The name of the differential to be substituted.
   :param val: The value to substitute *sym* with.
   :return: A new gdual with the substitution made.

------------------------------------------------------

.. cpp:function:: gdual subs(const std::string &sym, const gdual &val) const

   Substitute the differential *sym* with the gdual *val*

   :param sym: The name of the differential to be substituted.
   :param val: The value to substitute *sym* with.
   :return: A new gdual with the substitution made.

------------------------------------------------------

.. cpp:function:: gdual trim(double epsilon) const

   Sets to zero all coefficients of the gdual with absolute value smaller than *epsilon*.

   :param epsilon: Tolerance for the trim.
   :return: A new gdual without the trimmed coefficients.

------------------------------------------------------

.. cpp:function:: gdual extract_terms(unsigned degree) const

   Extracts all the monomials of a given *degree*.

   :param order: The monomials degree.
   :return:  A new gdual containing only the terms extracted, but preserving the truncation order of the original gdual.

   :exception: std::invalid_argument if the *degree* is higher than the gdual truncation order.

------------------------------------------------------

.. cpp:function:: void extend_symbol_set(const std::vector<std::string> &sym_vars)

   Adds some symbolic variables to the current obake::polynomial
   This is useful in situations where some differential :math:`dx` does not
   appear as its coefficient is zero but we still want to treat the gdual as a function of 
   :math:`x` too (for example when extracting the relative coefficient)

   :param sym_vars: list of symbolic names. It must contain all symbolic names of
     the current obake::polynomial. It may contain more. All symbols must start with the letter "d".

   :exception: std::invalid_argument if any symbol in *sym_vars* does not start with the letter "d"
     or if *sym_vars* does not contain all current symbols too.

------------------------------------------------------

gdual inspection
^^^^^^^^^^^^^^^^^

.. cpp:function:: auto get_symbol_set_size() const

   Returns the size of the symbol set of the obake::polynomial

------------------------------------------------------

.. cpp:function:: std::vector<std::string> get_symbol_set() const

   Returns the symbol set of the obake::polynomial stripping the differentials (i.e. "dx" becomes "x")

------------------------------------------------------

.. cpp:function:: auto evaluate(const std::unordered_map<std::string, double> &dict) const

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

   :param dict: a dictionary (unordered map) containing the values for the differentials.
   :return: the value of the Taylor polynomial.

   .. code-block:: c++

      gdual<double> x1(1., "x1", 2);
      gdual<double> x2(1., "x2", 2);
      auto f = x1*x1 + x2;
      std::cout << f.evaluate({{"dx1", 1.}, {"dx2", 1.}}) << "\n";

------------------------------------------------------

.. cpp:function:: auto degree() const

   Returns the degree of the polynomial. This is necessarily smaller or equal the truncation order.

   :return: the polynomial degree.

------------------------------------------------------

.. cpp:function:: unsigned int get_order() const

   Returns the truncation order of the polynomial.

   :return: the polynomial truncation order.

------------------------------------------------------

.. cpp:function:: template <typename T> auto find_cf(const T &c) const

   Returns the coefficient of the monomial specified by the container *c*.
   The container contains the exponents of the requested monomial. In a three
   variable Taylor expansion with :math:`x, y, z` as symbols, [1, 3, 2] would denote
   the monomial :math:`dx dy^3 dz^2`.

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

   .. warning::
     
     If the monomial requested is not found in the polynomial a zero coefficient will be returned.

   :return: the coefficient of the monomial, if present, zero otherwise.

   :exception: std::invalid_argument if the requested coefficient is beyond the truncation order
   
   .. code-block:: c++

      gdual<double> x1(1.2, "x1", 2);
      gdual<double> x2(-0.2, "x2", 2);
      auto f = sin(x1*x1 + x2);
      std::cout << f.find_cf(std::vector<double>({1,1})) << "\n";

------------------------------------------------------

.. cpp:function:: template <typename T> auto find_cf(std::initializer_list<T> l) const

   Returns the coefficient of the monomial specified by the initializer list *l*.

   .. note::
   
      This method overloads the one above and is provided for convenience.

   :return: the coefficient of the monomial, if present, zero otherwise.

   :exception: std::invalid_argument if the requested coefficient is beyond the truncation order
   
   .. code-block:: c++

      gdual<double> x1(1.2, "x1", 2);
      gdual<double> x2(-0.2, "x2", 2);
      auto f = sin(x1*x1 + x2);
      std::cout << f.find_cf({1,1}) << "\n";

------------------------------------------------------

.. cpp:function:: Cf constant_cf()

   Finds the constant coefficient of the Taylor polynomial. So that if :math:`T_{f} = f_0 + \hat f`, :math:`f_0` is returned

   :return: the constant coefficient.

------------------------------------------------------

.. cpp:function:: template <typename T> auto get_derivative(const T &c) const

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

   .. note:: 
   
      The container describes the derivative requested. In a three
      variable polynomial with :math:`x, y, z` as symbols, [1, 3, 2] would denote
      the sixth order derivative :math:`\frac{d^6}{dxdy^3dz^2}`.

   .. note::

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

   :return: the value of the derivative

   :exception: std::invalid_argument if the requested coefficient is beyond the truncation order

   .. code-block:: c++

      gdual<double> x1(1.2, "x1", 2);
      gdual<double> x2(-0.2, "x2", 2);
      auto f = sin(x1*x1 + x2);
      // This streams the value of df^2/dx1/dx2 in x1=1.2, x2 = -0.2
      std::cout << f.get_derivative(std::vector<double>({1,1})) << "\n";

------------------------------------------------------

.. cpp:function:: template <typename T> auto get_derivative(std::initializer_list<T> l) const

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

   .. note::
   
      This method overloads the one above and is provided for convenience.

   :return: the value of the derivative

   :exception: std::invalid_argument if the requested coefficient is beyond the truncation order

   .. code-block:: c++

      gdual<double> x1(1.2, "x1", 2);
      gdual<double> x2(-0.2, "x2", 2);
      auto f = sin(x1*x1 + x2);
      // This streams the value of df^2/dx1/dx2 in x1=1.2, x2 = -0.2
      std::cout << f.get_derivative({1,1}) << "\n";

------------------------------------------------------

.. cpp:function:: template <typename T> auto get_derivative(const std::unordered_map<std::string, unsigned int> &dict) const

   Returns the (mixed) derivative value of the order specified in *dict*.

   :param dict: a dictionary (unordered map) containing the derivation order (assumes zero for symbols not present).

   :return: the value of the derivative

   :exception: std::invalid_argument if the requested derivative degree is beyond the truncation order

   .. code-block:: c++

      gdual<double> x1(1.2, "x1", 2);
      gdual<double> x2(-0.2, "x2", 2);
      auto f = sin(x1*x1 + x2);
      // This streams the value of df^2/dx1/dx2 in x1=1.2, x2 = -0.2
      std::cout << f.get_derivative({{"dx1", 1u}, {"dx2", 1u}}) << "\n";

------------------------------------------------------

.. cpp:function:: bool is_zero(double tol) const

   Checks all coefficients of the gdual and returns true if all their absolute values are smaller
   than the defined tolerance *tol*.

   :return: whether the gdual is zero within a tolerance.

   .. code-block:: c++

      gdual<double> x(0.1, "x", 6);
      auto f = 1 - sin(x)*sin(x) - cos(x)*cos(x);
      if (f.is_zero(1e-13)) {
        std::cout << "The trigonomoetric identity holds!!" << std::endl;
      }

------------------------------------------------------

Operators
---------

The following operators are implemented: 

  * <<, streaming 
  * ==, equal to 
  * =,  assignement
  * !=, not equal to 
  * +=, addition assignment
  * -=, subtraction assignment
  * \*=, multiplication assignment
  * /=, division assignment
  * -, unary minus
  * +, unary plus
  * +, addition
  * -, subtraction
  * \*, multiplication
  * /, division

They allow to compute with the type gdual as you would operate with a basic type.

.. note:: 

   When relevant, the operators implement order promotion so that, for example, if a gdual of order 2 is added to a
   gdual of order 3 the resulting gdual will have order three.

We specify the documentation of a few operators with non trivial meaning.

------------------------------------------------------

.. cpp:function:: friend std::ostream &operator<<(std::ostream &os, const gdual &d)

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

   .. note::
      
      The print order of the terms will be undefined.

   :param os: target stream.
   :param d: gdual argument.

   :return: reference to *os*

------------------------------------------------------

.. cpp:function:: friend bool operator==(const gdual &d1, const gdual &d2)

   Equality operator. Two gduals are considered equal if all their coefficients are equal.

   .. note:: 
   
      The truncation order of *d1* and *d2* may be different

   :param d1: first argument.
   :param d2: second argument.

   :return: The result of the comparison.

------------------------------------------------------

.. cpp:function:: friend bool operator!=(const gdual &d1, const gdual &d2)

   Non equality operator.

   .. note:: 
   
      The truncation order of *d1* and *d2* may be different

   :param d1: first argument.
   :param d2: second argument.

   :return: The result of the comparison.