/*! \file func.hpp
 * \brief Real-valued functions defined here
 * 
 * A general class for a real-valued functions from \f$\mathbb{R}^n\f$
 * to \f$\mathbb{R}\f$ instead of using naked function pointers, meant
 * to be a class that can be derived from. Wrappers for GSL-style
 * function pointers are also provided.
 *
 *  \file func.cpp
 *  \brief Implementations and instantiations of the functions and
 *  classes declared in func.hpp
 */
#ifndef __FUNC_HPP__
#define __FUNC_HPP__

#include "linalg.hpp"
#include "error.hpp"
#include <gsl/gsl_deriv.h>

namespace bvp{
  using namespace linalg;

  /*! \brief  A real-valued function from \f$\mathbb{R}^n\f$ to
   *  \f$\mathbb{R}\f$ .
   *
   * That is, a realfunc is a real-valued function that takes in
   * linalg::point and returns doubles. This class is meant to be
   * generic and any other real-valued function (e.g. a
   * radial_basis_function) should derive from this one.
   *
   * Note that the base class provides derivatives using
   * finite-difference approximations from the GSL. Derived classes
   * are intended to provide their own exact derivatives instead of
   * using the default finite-difference approximations.
   * 
   */
  class realfunc{
  public:
    /// Create a useless realfunc, not assigned to anything.
    realfunc(); 
    /// Create a realfunc from a function pointer.
    realfunc( double(*f)(const point&));
    /// Destroys nothing, but declared virtual for derived classes.
    virtual ~realfunc(){};

    /// Give a function pointer to this realfunc.
    void set_function_ptr(double (f_in)(const point &p));
    
    /// Evaluate the function at point
    double operator()(const point& p) const;
    /// Evaluate the function at point
    virtual double at(const point& p) const;

    /// First derivative at point \f$x\f$ and in the \f$k\f$th direction.
    virtual double d(const point& x, size_t k) const; 
    /// Second derivative (FIXME: actually implement this)
    virtual double d2(const point& x, size_t k1, size_t k2) const ; 
  protected:
    /// Machine epsilon
    static double eps; 
    /// Square root of epsilon
    static double sqrteps;
    /// Cube root of epsilon
    static double root3eps;
    /// Have already initialised epsilon above?
    static bool initialised ;
  private:
    /// Pointer to a function that takes linalg::point and returns double
    double (*myfunc)(const point &p);
    /// Exception builder
    error_handling::badArgument
    no_init(int line, string file) const; 
  };


  ///A function wrapper for calling GSL derivatives. 
  class gsl_function_wrapper{
  public:
    /*! \brief Wraps a realfunc for the GSL
     *
     * Turns a realfunc into a one-dimensional GSL function from
     * \f$\mathbb{R}\f$ to \f$\mathbb{R}\f$ by localising the function
     * at a point along a given direction.
     *
     * \param f - The realfunc to wrap.
     * \param p - The point on which to localise.
     * \param idx - The direction along with to localise.
     */
    gsl_function_wrapper(const realfunc &f, point p, size_t idx);
    /// Same as the above constructor
    void set_params( const realfunc &f, point p, size_t idx);
    /// Pointer to the GSL-style function
    gsl_function* get_gsl_function() const;

    /// The actual function that's passed to the GSL
    static double takemyaddress(double xi, void* nothing);
  private:
    /// No nontrivial construction!
    gsl_function_wrapper();
    
    /// The point at which the function is localised
    static point x;
    /// The direction along which the function is localised
    static size_t index;    
    /// The function being localised
    static realfunc myfunc;
    
    /// The GSL function formed with the above data.
    static gsl_function* f;
  };
}

#endif //__FUNC_HPP__