AutomaticDifferentiation/AutomaticDifferentiation_static.hpp

#ifndef DEF_AUTOMATIC_DIFFERENTIATION_STATIC
#define DEF_AUTOMATIC_DIFFERENTIATION_STATIC

#undef __Dual_DualBase
#undef __Dual_VectorT
#undef __Dual_bdynamic
#define __Dual_DualBase DualS
#define __Dual_VectorT Eigen::Array<Scalar, N, 1>
#define __Dual_bdynamic 0

#include "AutomaticDifferentiation_base.hpp"

/// Function object to evaluate the derivative of a univariate function anywhere without explicitely using the dual numbers.
///
/// Here is an example of use :
///
///     template<typename T> fct(T x) { return exp(-x*x); }
///
///     Scalar x = 3.14;
///     auto gradFct = GradFunc1(fct, x);    // x is only passed to the function so that the full type of GradFunc1 does not need to be written manually.
///     Scalar dfdx = gradFct(x);            // Evaluation of the gradient of fct at x.
///
///     // Alternatively :
///     Scalar fx, dfdx;
///     gradFct.get_f_grad(x, fx, dfdx);     // Evaluation of value AND the gradient of fct at x.
template<typename Func, typename Scalar>
struct GradFunc1
{
    Func f;

    /// Constructor of the gradient function object. The second parameter IS NOT USED FOR THE COMPUTATION OF THE GRADIENT, but only as a convenience for the
    /// automatic detection of the Scalar variable type.
    GradFunc1(Func f_, Scalar) : f(f_) {}

	// Function that returns the 1st derivative of the function f at x.
    Scalar operator()(Scalar const& x)
    {
        Scalar fx, gradfx;
        get_f_grad(x, fx, gradfx);
        return gradfx;
    }

	/// Function that returns both the function value and the gradient of f at x. Use this preferably over separate calls to f and to gradf.
    void get_f_grad(Scalar const& x, Scalar & fx, Scalar & gradfx)
    {
        // differentiate using the dual number and return the .b component.
        DualS<Scalar, 1> X(x);
		X.diff(0);
		DualS<Scalar, 1> Y = f(X);
		fx = Y.x();
		gradfx = Y.d(0);
    }
};

/// Macro to create a function object that returns the gradient of the function at X.
/// Suitable for functions that take an scalar and return a scalar.
/// Designed to work with functions, lambdas, etc.
#define CREATE_GRAD_FUNCTION_OBJECT_1_1(Func, GradFuncName)                     \
struct GradFuncName {	                                                        \
    template<typename Scalar>											        \
    Scalar operator()(Scalar const& x) { 									    \
        Scalar fx, gradfx;                                                      \
        get_f_grad(x, fx, gradfx);                                              \
        return gradfx;														    \
    }																			\
    template<typename Scalar>													\
    void get_f_grad(Scalar const& x, Scalar & fx, Scalar & gradfx) { 			\
        DualS<Scalar,1> X(x);													\
	    X.diff(0);														        \
	    DualS<Scalar,1> Y = Func<DualS<Scalar,1>>(X);						    \
	    fx = Y.x();															    \
	    gradfx = Y.d(0);													    \
    }																			\
}

/// Macro to create a function object that returns the gradient of the function at X.
/// Suitable for functions that take an N-dimensional vector and return a scalar.
/// This version uses statically allocated arrays for everything (Dual and vectors to and from the function).
/// Designed to work with functions, lambdas, etc.
#define CREATE_GRAD_FUNCTION_OBJECT_N_1_S(Func, GradFuncName, N)                        \
template<typename Scalar>                                                               \
struct GradFuncName {	                                                                \
    using ArrayOfScalar = Eigen::Array<Scalar, N, 1>;                                   \
    using ArrayOfDual = Eigen::Array<DualS<Scalar,N>, N, 1>;                            \
                                                                                        \
    ArrayOfScalar operator()(ArrayOfScalar const& x) {                                  \
        Scalar fx;                                                                      \
        ArrayOfScalar gradfx;                                                           \
        get_f_grad(x, fx, gradfx);                                                      \
        return gradfx;														            \
    }																			        \
    void get_f_grad(ArrayOfScalar const& x, Scalar & fx, ArrayOfScalar & gradfx) {      \
        ArrayOfDual X(N);                                                               \
        for(int i = 0 ; i < N ; i++)                                                    \
        {                                                                               \
            X[i] = x[i];                                                                \
            X[i].diff(i);                                                               \
        }                                                                               \
        auto Y = Func<DualS<Scalar,N>>(X);                                              \
        for (int i = 0; i < N; i++)                                                     \
            gradfx[i] = Y.d(i);                                                         \
	    fx = Y.x();															            \
    }																			        \
}

#endif