[/          Copyright Matthew Pulver 2018 - 2019.
// Distributed under the Boost Software License, Version 1.0.
//    (See accompanying file LICENSE_1_0.txt or copy at
//          https://www.boost.org/LICENSE_1_0.txt)]

[section:autodiff Automatic Differentiation]
[template autodiff_equation[name]  '''<inlinemediaobject><imageobject><imagedata fileref="../equations/autodiff/'''[name]'''"></imagedata></imageobject></inlinemediaobject>''']

[h1:synopsis Synopsis]

    #include <boost/math/differentiation/autodiff.hpp>
    
    namespace boost {
    namespace math {
    namespace differentiation {
    
    // Function returning a single variable of differentiation. Recommended: Use auto for type.
    template <typename RealType, size_t Order, size_t... Orders>
    autodiff_fvar<RealType, Order, Orders...> make_fvar(RealType const& ca);
    
    // Function returning multiple independent variables of differentiation in a std::tuple.
    template<typename RealType, size_t... Orders, typename... RealTypes>
    auto make_ftuple(RealTypes const&... ca);
    
    // Type of combined autodiff types. Recommended: Use auto for return type (C++14).
    template <typename RealType, typename... RealTypes>
    using promote = typename detail::promote_args_n<RealType, RealTypes...>::type;
    
    namespace detail {
    
    // Single autodiff variable. Use make_fvar() or make_ftuple() to instantiate.
    template <typename RealType, size_t Order>
    class fvar {
     public:
      // Query return value of function to get the derivatives.
      template <typename... Orders>
      get_type_at<RealType, sizeof...(Orders) - 1> derivative(Orders... orders) const;
    
      // All of the arithmetic and comparison operators are overloaded.
      template <typename RealType2, size_t Order2>
      fvar& operator+=(fvar<RealType2, Order2> const&);
    
      fvar& operator+=(root_type const&);
    
      // ...
    };
    
    // Standard math functions are overloaded and called via argument-dependent lookup (ADL).
    template <typename RealType, size_t Order>
    fvar<RealType, Order> floor(fvar<RealType, Order> const&);
    
    template <typename RealType, size_t Order>
    fvar<RealType, Order> exp(fvar<RealType, Order> const&);
    
    // ...
    
    }  // namespace detail
    
    }  // namespace differentiation
    }  // namespace math
    }  // namespace boost

[h1:description Description]

Autodiff is a header-only C++ library that facilitates the [@https://en.wikipedia.org/wiki/Automatic_differentiation
automatic differentiation] (forward mode) of mathematical functions of single and multiple variables.

This implementation is based upon the [@https://en.wikipedia.org/wiki/Taylor_series Taylor series] expansion of
an analytic function /f/ at the point ['x[sub 0]]:

[/ Thanks to http://www.tlhiv.org/ltxpreview/ for LaTeX-to-SVG conversions. ]
[/ \Large\begin{align*}
f(x_0+\varepsilon) &= f(x_0) + f'(x_0)\varepsilon + \frac{f''(x_0)}{2!}\varepsilon^2 + \frac{f'''(x_0)}{3!}\varepsilon^3 + \cdots \\
  &= \sum_{n=0}^N\frac{f^{(n)}(x_0)}{n!}\varepsilon^n + O\left(\varepsilon^{N+1}\right).
\end{align*} ]
[:[:[autodiff_equation taylor_series.svg]]]

The essential idea of autodiff is the substitution of numbers with polynomials in the evaluation of ['f(x[sub 0])].
By substituting the number ['x[sub 0]] with the first-order polynomial ['x[sub 0]+\u03b5], and using the same
algorithm to compute ['f(x[sub 0]+\u03b5)], the resulting polynomial in ['\u03b5] contains the function's derivatives
['f'(x[sub 0])], ['f''(x[sub 0])], ['f\'\'\'(x[sub 0])], ... within the coefficients. Each coefficient is equal to the
derivative of its respective order, divided by the factorial of the order.

In greater detail, assume one is interested in calculating the first /N/ derivatives of /f/ at ['x[sub 0]]. Without
loss of precision to the calculation of the derivatives, all terms ['O(\u03b5[super N+1])] that include powers
of ['\u03b5] greater than /N/ can be discarded. (This is due to the fact that each term in a polynomial depends
only upon equal and lower-order terms under arithmetic operations.) Under these truncation rules, /f/ provides a
polynomial-to-polynomial transformation:

[/ \Large$$f \qquad : \qquad x_0+\varepsilon \qquad \mapsto \qquad
    \sum_{n=0}^Ny_n\varepsilon^n=\sum_{n=0}^N\frac{f^{(n)}(x_0)}{n!}\varepsilon^n.$$ ]
[:[:[autodiff_equation polynomial_transform.svg]]]

C++'s ability to overload operators and functions allows for the creation of a class `fvar` ([_f]orward-mode autodiff
[_var]iable) that represents polynomials in ['\u03b5]. Thus the same algorithm /f/ that calculates the numeric value
of ['y[sub 0]=f(x[sub 0])], when written to accept and return variables of a generic (template) type, is also used
to calculate the polynomial ['\u03a3[sub n]y[sub n]\u03b5[super n]=f(x[sub 0]+\u03b5)].  The derivatives
['f[super (n)](x[sub 0])] are then found from the product of the respective factorial ['n!] and coefficient
['y[sub n]]:

[/ \Large$$\frac{d^nf}{dx^n}(x_0)=n!y_n.$$ ]
[:[:[autodiff_equation derivative_formula.svg]]]


[h1:examples Examples]

[h2:example-single-variable Example 1: Single-variable derivatives]

[h3 Calculate derivatives of ['f(x)=x[super 4]] at /x/=2.]

In this example, `make_fvar<double, Order>(2.0)` instantiates the polynomial 2+['\u03b5]. The `Order=5`
means that enough space is allocated (on the stack) to hold a polynomial of up to degree 5 during the proceeding
computation.

Internally, this is modeled by a `std::array<double,6>` whose elements `{2, 1, 0, 0, 0, 0}` correspond
to the 6 coefficients of the polynomial upon initialization. Its fourth power, at the end of the computation, is
a polynomial with coefficients `y = {16, 32, 24, 8, 1, 0}`.  The derivatives are obtained using the formula
['f[super (n)](2)=n!*y[n]].

    #include <boost/math/differentiation/autodiff.hpp>
    #include <iostream>
    
    template <typename T>
    T fourth_power(T const& x) {
      T x4 = x * x;  // retval in operator*() uses x4's memory via NRVO.
      x4 *= x4;      // No copies of x4 are made within operator*=() even when squaring.
      return x4;     // x4 uses y's memory in main() via NRVO.
    }
    
    int main() {
      using namespace boost::math::differentiation;
    
      constexpr unsigned Order = 5;                  // Highest order derivative to be calculated.
      auto const x = make_fvar<double, Order>(2.0);  // Find derivatives at x=2.
      auto const y = fourth_power(x);
      for (unsigned i = 0; i <= Order; ++i)
        std::cout << "y.derivative(" << i << ") = " << y.derivative(i) << std::endl;
      return 0;
    }
    /*
    Output:
    y.derivative(0) = 16
    y.derivative(1) = 32
    y.derivative(2) = 48
    y.derivative(3) = 48
    y.derivative(4) = 24
    y.derivative(5) = 0
    */

The above calculates

[/ \Large\begin{alignat*}{3}
{\tt y.derivative(0)} &=& f(2) =&& \left.x^4\right|_{x=2} &= 16\\
{\tt y.derivative(1)} &=& f'(2) =&& \left.4\cdot x^3\right|_{x=2} &= 32\\
{\tt y.derivative(2)} &=& f''(2) =&& \left.4\cdot 3\cdot x^2\right|_{x=2} &= 48\\
{\tt y.derivative(3)} &=& f'''(2) =&& \left.4\cdot 3\cdot2\cdot x\right|_{x=2} &= 48\\
{\tt y.derivative(4)} &=& f^{(4)}(2) =&& 4\cdot 3\cdot2\cdot1 &= 24\\
{\tt y.derivative(5)} &=& f^{(5)}(2) =&& 0 &
\end{alignat*} ]
[:[:[autodiff_equation example1.svg]]]

[h2:example-multiprecision
Example 2: Multi-variable mixed partial derivatives with multi-precision data type]

[/ \Large$\frac{\partial^{12}f}{\partial w^{3}\partial x^{2}\partial y^{4}\partial z^{3}}(11,12,13,14)$]
[/ \Large$f(w,x,y,z)=\exp\left(w\sin\left(\frac{x\log(y)}{z}\right)+\sqrt{\frac{wz}{xy}}\right)+\frac{w^2}{\tan(z)}$]
[h3 Calculate [autodiff_equation mixed12.svg] with a precision of about 50 decimal digits,
where [autodiff_equation example2f.svg].]

In this example, `make_ftuple<float50, Nw, Nx, Ny, Nz>(11, 12, 13, 14)` returns a `std::tuple` of 4
independent `fvar` variables, with values of 11, 12, 13, and 14, for which the maximum order derivative to
be calculated for each are 3, 2, 4, 3, respectively. The order of the variables is important, as it is the same
order used when calling `v.derivative(Nw, Nx, Ny, Nz)` in the example below.

    #include <boost/math/differentiation/autodiff.hpp>
    #include <boost/multiprecision/cpp_bin_float.hpp>
    #include <iostream>
    
    using namespace boost::math::differentiation;
    
    template <typename W, typename X, typename Y, typename Z>
    promote<W, X, Y, Z> f(const W& w, const X& x, const Y& y, const Z& z) {
      using namespace std;
      return exp(w * sin(x * log(y) / z) + sqrt(w * z / (x * y))) + w * w / tan(z);
    }
    
    int main() {
      using float50 = boost::multiprecision::cpp_bin_float_50;
    
      constexpr unsigned Nw = 3;  // Max order of derivative to calculate for w
      constexpr unsigned Nx = 2;  // Max order of derivative to calculate for x
      constexpr unsigned Ny = 4;  // Max order of derivative to calculate for y
      constexpr unsigned Nz = 3;  // Max order of derivative to calculate for z
      // Declare 4 independent variables together into a std::tuple.
      auto const variables = make_ftuple<float50, Nw, Nx, Ny, Nz>(11, 12, 13, 14);
      auto const& w = std::get<0>(variables);  // Up to Nw derivatives at w=11
      auto const& x = std::get<1>(variables);  // Up to Nx derivatives at x=12
      auto const& y = std::get<2>(variables);  // Up to Ny derivatives at y=13
      auto const& z = std::get<3>(variables);  // Up to Nz derivatives at z=14
      auto const v = f(w, x, y, z);
      // Calculated from Mathematica symbolic differentiation.
      float50 const answer("1976.319600747797717779881875290418720908121189218755");
      std::cout << std::setprecision(std::numeric_limits<float50>::digits10)
                << "mathematica   : " << answer << '\n'
                << "autodiff      : " << v.derivative(Nw, Nx, Ny, Nz) << '\n'
                << std::setprecision(3)
                << "relative error: " << (v.derivative(Nw, Nx, Ny, Nz) / answer - 1) << '\n';
      return 0;
    }
    /*
    Output:
    mathematica   : 1976.3196007477977177798818752904187209081211892188
    autodiff      : 1976.3196007477977177798818752904187209081211892188
    relative error: 2.67e-50
    */

[h2:example-black_scholes
Example 3: Black-Scholes Option Pricing with Greeks Automatically Calculated]
[h3 Calculate greeks directly from the Black-Scholes pricing function.]

Below is the standard Black-Scholes pricing function written as a function template, where the price, volatility
(sigma), time to expiration (tau) and interest rate are template parameters. This means that any greek based on
these 4 variables can be calculated using autodiff. The below example calculates delta and gamma where the variable
of differentiation is only the price. For examples of more exotic greeks, see `example/black_scholes.cpp`.

```
#include <boost/math/differentiation/autodiff.hpp>
#include <iostream>

using namespace boost::math::constants;
using namespace boost::math::differentiation;

// Equations and function/variable names are from
// https://en.wikipedia.org/wiki/Greeks_(finance)#Formulas_for_European_option_Greeks

// Standard normal cumulative distribution function
template <typename X>
X Phi(X const& x) {
  return 0.5 * erfc(-one_div_root_two<X>() * x);
}

enum class CP { call, put };

// Assume zero annual dividend yield (q=0).
template <typename Price, typename Sigma, typename Tau, typename Rate>
promote<Price, Sigma, Tau, Rate> black_scholes_option_price(CP cp,
                                                            double K,
                                                            Price const& S,
                                                            Sigma const& sigma,
                                                            Tau const& tau,
                                                            Rate const& r) {
  using namespace std;
  auto const d1 = (log(S / K) + (r + sigma * sigma / 2) * tau) / (sigma * sqrt(tau));
  auto const d2 = (log(S / K) + (r - sigma * sigma / 2) * tau) / (sigma * sqrt(tau));
  switch (cp) {
    case CP::call:
      return S * Phi(d1) - exp(-r * tau) * K * Phi(d2);
    case CP::put:
      return exp(-r * tau) * K * Phi(-d2) - S * Phi(-d1);
  }
}

int main() {
  double const K = 100.0;                    // Strike price.
  auto const S = make_fvar<double, 2>(105);  // Stock price.
  double const sigma = 5;                    // Volatility.
  double const tau = 30.0 / 365;             // Time to expiration in years. (30 days).
  double const r = 1.25 / 100;               // Interest rate.
  auto const call_price = black_scholes_option_price(CP::call, K, S, sigma, tau, r);
  auto const put_price = black_scholes_option_price(CP::put, K, S, sigma, tau, r);

  std::cout << "black-scholes call price = " << call_price.derivative(0) << '\n'
            << "black-scholes put  price = " << put_price.derivative(0) << '\n'
            << "call delta = " << call_price.derivative(1) << '\n'
            << "put  delta = " << put_price.derivative(1) << '\n'
            << "call gamma = " << call_price.derivative(2) << '\n'
            << "put  gamma = " << put_price.derivative(2) << '\n';
  return 0;
}
/*
Output:
black-scholes call price = 56.5136
black-scholes put  price = 51.4109
call delta = 0.773818
put  delta = -0.226182
call gamma = 0.00199852
put  gamma = 0.00199852
*/
```

[h1 Advantages of Automatic Differentiation]
The above examples illustrate some of the advantages of using autodiff:

* Elimination of code redundancy. The existence of /N/ separate functions to calculate derivatives is a form
  of code redundancy, with all the liabilities that come with it:
    * Changes to one function require /N/ additional changes to other functions. In the 3rd example above,
      consider how much larger and inter-dependent the above code base would be if a separate function were
      written for [@https://en.wikipedia.org/wiki/Greeks_(finance)#Formulas_for_European_option_Greeks
      each Greek] value.
    * Dependencies upon a derivative function for a different purpose will break when changes are made to
      the original function. What doesn't need to exist cannot break.
    * Code bloat, reducing conceptual integrity. Control over the evolution of code is easier/safer when
      the code base is smaller and able to be intuitively grasped.
* Accuracy of derivatives over finite difference methods. Single-iteration finite difference methods always include
  a ['\u0394x] free variable that must be carefully chosen for each application. If ['\u0394x] is too small, then
  numerical errors become large. If ['\u0394x] is too large, then mathematical errors become large.  With autodiff,
  there are no free variables to set and the accuracy of the answer is generally superior to finite difference
  methods even with the best choice of ['\u0394x].

[h1 Manual]
Additional details are in the [@../differentiation/autodiff.pdf autodiff manual].

[endsect]