/*! \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
 */
#pragma once

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

namespace kwantix{

/// \defgroup BVP

/*! \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
 * kwantix::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 kwantix::point and returns double
  double (*myfunc)(const point &p);
  /// Exception builder
  kwantix::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;
};

/// \}

} //namespace kwantix