/*****************************************************************************
* File: Derivative.hh
* Author: Keith Schwarz (htiek@cs.stanford.edu)
*
* A function that produces a function object containing an approximation of
* the derivative of a real-valued function.  Internally, the function works
* by handing back a function object that uses a secant-line approximation of
* the derivative, as shown here:
*
*            f(x + e) - f(x - e)
* f'(x) ~=   -------------------
*                     2e
*
* This can be shown to produce the derivative with a very good approximation
* for small e.
*/

#ifndef Derivative_Included
#define Derivative_Included

/* Forward declarations to make TakeDerivative make sense. */
namespace derivative_detail {
template <typename UnaryFunction> class Derivative;
}

/**
* Function: TakeDerivative(UnaryFunction function, double epsilon = 0.000001);
* Usage: transform(v.begin(), v.end(), v.begin(), TakeDerivative(sinc));
* --------------------------------------------------------------------------
* Given a unary function (and optionally an epsilon), returns a function object
* containing an approximation of the first derivative of that function.
*/

template <typename UnaryFunction>
derivative_detail::Derivative<UnaryFunction>
TakeDerivative(UnaryFunction function, double epsilon = 0.000001);

/* * * * * Implementation Below This Point * * * * */
namespace derivative_detail {
/* Functor class responsible for approximating a derivative. */
template <typename UnaryFunction> class Derivative {
public:
/* Constructor stores the function and epsilon to use. */
Derivative(UnaryFunction function, double epsilon)
: mFunction(function), mEpsilon(epsilon) {
// Handled in initializer list.
}

/* Functor operator produces a derivative approximation. */
double operator() (double x) const {
return (mFunction(x + mEpsilon) - mFunction(x - mEpsilon)) / (2.0 * mEpsilon);
}

private:
/* The function to call and the obligatory epsilon. */
UnaryFunction mFunction;
double mEpsilon;
};
}

/* Actual implementation of TakeDerivative. */
template <typename UnaryFunction>
derivative_detail::Derivative<UnaryFunction>
TakeDerivative(UnaryFunction function, double epsilon) {
/* Construct a Derivative of the correct type, then hand it back. */
return derivative_detail::Derivative<UnaryFunction>(function, epsilon);
}

#endif