#ifndef __JSPLINE__
#define __JSPLINE__

#include <utility>
#include <limits>

#include "JLang/JException.hh"
#include "JLang/JClass.hh"
#include "JMath/JMath.hh"
#include "JTools/JCollection.hh"
#include "JTools/JFunctional.hh"
#include "JTools/JFunctionalMap.hh"
#include "JTools/JElement.hh"


namespace JTOOLS {

  namespace {
    using JLANG::JClass;
    using JLANG::JNoValue;
    using JLANG::JDivisionByZero;
    using JLANG::JEmptyCollection;
    using JLANG::JValueOutOfRange;
    using JMATH::JMath;
  }


  /**
   * Customized data structure for extended result of spline interpolation.
   * This data structure contains the interpolated values of the function, the first derivative and the integral(s).
   */
  template<class T>
  class JFunctionValuePDF : public JMath< JFunctionValuePDF<T> > {
  public:
    /**
     * Default constructor.
     */
    JFunctionValuePDF() :
      y(0.0),
      z(0.0),
      v(0.0),
      V(0.0)
    {}


    /**
     * Constructor.
     *
     * \param  __y   function value
     */
    JFunctionValuePDF(const T& __y) :
      y(__y),
      z(0.0),
      v(0.0),
      V(0.0)
    {}


    /**
     * Constructor.
     *
     * \param  __y   function value
     * \param  __z   first derivative
     * \param  __v   integral <xmin,x]
     * \param  __V   integral <xmin,xmax>
     */
    JFunctionValuePDF(typename JClass<T>::argument_type __y,
		      typename JClass<T>::argument_type __z,
		      typename JClass<T>::argument_type __v,
		      typename JClass<T>::argument_type __V) :
      y(__y),
      z(__z),
      v(__v),
      V(__V)
    {}


    /**
     * Prefix unary minus for JFunctionValuePDF.
     *
     * \return         this JFunctionValuePDF
     */
    const JFunctionValuePDF& operator-()
    {
      y = -y;
      z = -z;
      v = -v;
      V = -V;

      return *this;
    }


    /**
     * Addition operator for JFunctionValuePDF.
     *
     * \param  b       that JFunctionValuePDF
     * \return         this JFunctionValuePDF
     */
    const JFunctionValuePDF& operator+=(const JFunctionValuePDF& b)
    {
      y += b.y;
      z += b.z;
      v += b.v;
      V += b.V;

      return *this;
    }


    /**
     * Subtraction operator for JFunctionValuePDF.
     *
     * \param  b       that JFunctionValuePDF
     * \return         this JFunctionValuePDF
     */
    const JFunctionValuePDF& operator-=(const JFunctionValuePDF& b)
    {
      y -= b.y;
      z -= b.z;
      v -= b.v;
      V -= b.V;

      return *this;
    }


    /**
     * Multiplication operator for JFunctionValuePDF.
     *
     * \param  b       double
     * \return         this JFunctionValuePDF
     */
    const JFunctionValuePDF& operator*=(const double& b)
    {
      y *= b;
      z *= b;
      v *= b;
      V *= b;

      return *this;
    }


    /**
     * Division operator for JFunctionValuePDF.
     *
     * \param  b       value
     * \return         this JFunctionValuePDF
     */
    const JFunctionValuePDF& operator/=(const double& b)
    {
      y /= b;
      z /= b;
      v /= b;
      V /= b;

      return *this;
    }


    T y;   //!< function value
    T z;   //!< first derivative
    T v;   //!< integral <xmin,x]
    T V;   //!< integral <xmin,xmax>
  };


  /**
   * Customized data structure for extended result of spline interpolation.
   * This data structure contains the interpolated values of the function, the first and second derivatives.
   */
  template<class T>
  class JFunctionValueHesse : public JMath< JFunctionValueHesse<T> > {
  public:
    /**
     * Default constructor.
     */
    JFunctionValueHesse() :
      f  (0.0),
      fp (0.0),
      fpp(0.0)
    {}


    /**
     * Constructor.
     *
     * \param  __f   function value
     */
    JFunctionValueHesse(typename JClass<T>::argument_type __f) :
      f  (__f),
      fp (0.0),
      fpp(0.0)
    {}


    /**
     * Constructor.
     *
     * \param  __f   function value
     * \param  __fp  first  derivative
     * \param  __fpp second derivative
     */
    JFunctionValueHesse(typename JClass<T>::argument_type __f,
			typename JClass<T>::argument_type __fp,
			typename JClass<T>::argument_type __fpp) :
      f  (__f  ),
      fp (__fp ),
      fpp(__fpp)
    {}


    /**
     * Prefix unary minus for JFunctionValueHesse.
     *
     * \return         this JFunctionValueHesse
     */
    const JFunctionValueHesse& operator-()
    {
      f   = -f;
      fp  = -fp;
      fpp = -fpp;

      return *this;
    }


    /**
     * Addition operator for JFunctionValueHesse.
     *
     * \param  b       that JFunctionValueHesse
     * \return         this JFunctionValueHesse
     */
    const JFunctionValueHesse& operator+=(const JFunctionValueHesse& b)
    {
      f   += b.f;
      fp  += b.fp;
      fpp += b.fpp;

      return *this;
    }


    /**
     * Subtraction operator for JFunctionValueHesse.
     *
     * \param  b       that JFunctionValueHesse
     * \return         this JFunctionValueHesse
     */
    const JFunctionValueHesse& operator-=(const JFunctionValueHesse& b)
    {
      f   -= b.f;
      fp  -= b.fp;
      fpp -= b.fpp;

      return *this;
    }


    /**
     * Multiplication operator for JFunctionValueHesse.
     *
     * \param  b       double
     * \return         this JFunctionValueHesse
     */
    const JFunctionValueHesse& operator*=(const double& b)
    {
      f   *= b;
      fp  *= b;
      fpp *= b;

      return *this;
    }


    /**
     * Division operator for JFunctionValueHesse.
     *
     * \param  b       double
     * \return         this JFunctionValueHesse
     */
    const JFunctionValueHesse& operator/=(const double& b)
    {
      f   /= b;
      fp  /= b;
      fpp /= b;

      return *this;
    }


    T f;      //!< function value
    T fp;     //!< first  derivative
    T fpp;    //!< second derivative
  };


  /**
   * Auxiliary class to define first derivates of the spline function at two extrema.
   */
  template<class T>
  class JSplineBounds {
  public:
    /** 
     * Default constructor.
     */
    JSplineBounds() :
      fp_at_xmin(false, 0.0),
      fp_at_xmax(false, 0.0)
    {}


    /** 
     * Constructor.
     *
     * \param fpAtXmin    first derivative at smallest x
     * \param fpAtXmax    first derivative at largest  x
     */
    JSplineBounds(typename JClass<T>::argument_type fpAtXmin,
		  typename JClass<T>::argument_type fpAtXmax) :
      fp_at_xmin(true, fpAtXmin),
      fp_at_xmax(true, fpAtXmax)
    {}


    /**
     * Set first derivative of function at smallest x.
     * 
     * \param fpAtXmin    first derivative at smallest x
     */
    void setFirstDerivativeAtXmin(typename JClass<T>::argument_type fpAtXmin) 
    {
      fp_at_xmin.first  = true;
      fp_at_xmin.second = fpAtXmin;
    }


    /**
     * Set first derivative of function at largest  x.
     * 
     * \param fpAtXmax    first derivative at largest x
     */
    void setFirstDerivativeAtXmax(typename JClass<T>::argument_type fpAtXmax) 
    {
      fp_at_xmax.first  = true;
      fp_at_xmax.second = fpAtXmax;
    }


    /**
     * Has first derivative of function at smallest x.
     * 
     * \return        true if first derivative of function is set, else false
     */
    const bool& hasFirstDerivativeAtXmin() const 
    { 
      return fp_at_xmin.first;
    }


    /**
     * Has first derivative of function at largest  x.
     * 
     * \return        true if first derivative of function is set, else false
     */
    const bool& hasFirstDerivativeAtXmax() const 
    { 
      return fp_at_xmax.first;
    }


    /**
     * Get first derivative of function at smallest x.
     * 
     * \return        first derivative of function
     */
    const T& getFirstDerivativeAtXmin() const 
    { 
      if (fp_at_xmin.first)
	return fp_at_xmin.second;
      else
	throw JNoValue("JSplineBounds: missing first derivative.");
    }


    /**
     * Get first derivative of function at largest  x.
     * 
     * \return        first derivative of function
     */
    const T& getFirstDerivativeAtXmax() const 
    { 
      if (fp_at_xmax.first)
	return fp_at_xmax.second;
      else
	throw JNoValue("JSplineBounds: missing first derivative.");
    }

  protected:
    std::pair<bool,T> fp_at_xmin;
    std::pair<bool,T> fp_at_xmax;
  };


  /**
   * Template interface definition for compilation of spline object.
   *
   * Spline interpolation code is taken from reference:
   * Numerical Recipes in C++, W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery,
   * Cambridge University Press.
   */
  template<class JElement_t, template<class> class JCollection_t>
  class JSpline :
    public JCollection_t<JElement_t>,
    public virtual JFunction<>
  {
  public:

    typedef typename JElement_t::key_type                                        key_type;
    typedef typename JElement_t::mapped_type                                     mapped_type;

    typedef JCollection_t<JElement_t>                                            JContainer_t;

    typedef typename JContainer_t::const_iterator                                const_iterator;
    typedef typename JContainer_t::const_reverse_iterator                        const_reverse_iterator;

    typedef typename JContainer_t::iterator                                      iterator;
    typedef typename JContainer_t::reverse_iterator                              reverse_iterator;
 

    /**
     * Determination of second derivatives with no bounds.
     */
    virtual void compile() 
    {
      compile(JSplineBounds<mapped_type>());
    }


    /**
     * Determination of second derivatives with specified bounds.
     *
     * \param bounds    first derivatives at two extrema.
     */
    virtual void compile(const JSplineBounds<mapped_type>& bounds) 
    {
      const int numberOfElements = this->size();

      using namespace std;

      if (numberOfElements > 2) {
	
	std::vector<double> buffer(numberOfElements);

	if (bounds.hasFirstDerivativeAtXmin()) {
	  
	  iterator j = this->begin();
	  iterator i = j++;
	   
	  const double      dx = i->getX() - j->getX();
	  const mapped_type dy = i->getY() - j->getY();
	   
	  buffer[0] = -0.5;

	  i->set2ndDerivative((3.0/dx) * (dy/dx - bounds.getFirstDerivativeAtXmin()));

	} else {

	  buffer[0] = 0.0;

	  this->begin()->set2ndDerivative(0.0);
	}

	int index = 1;
       
	for (iterator k = this->begin(), i = k++, j = k++; k != this->end(); ++i, ++j, ++k, ++index) {

	  const key_type&    x1 = i->getX();
	  const key_type&    x2 = j->getX();
	  const key_type&    x3 = k->getX();

	  const mapped_type& y1 = i->getY();
	  const mapped_type& y2 = j->getY();
	  const mapped_type& y3 = k->getY();

	  const double sig = (x2 - x1) / (x3 - x1);
	  const double h   = sig * buffer[index-1] + 2.0;
	 
	  buffer[index]    = (sig - 1.0) / h;
	 
	  j->set2ndDerivative((y3-y2) / (x3-x2)  -  (y2-y1) / (x2-x1));

	  j->set2ndDerivative((6.0 * j->get2ndDerivative() / (x3-x1)  -  sig * i->get2ndDerivative()) / h);
	}
       

	if (bounds.hasFirstDerivativeAtXmax()) {

	  reverse_iterator j = this->rbegin();
	  reverse_iterator i = j++;

	  index = numberOfElements - 2;
	  
	  const double      dx = i->getX() - j->getX();
	  const mapped_type dy = i->getY() - j->getY();
	   
	  i->set2ndDerivative((3.0/dx) * (bounds.getFirstDerivativeAtXmax() - dy/dx));

	  i->set2ndDerivative((i->get2ndDerivative() - 0.5*j->get2ndDerivative()) / (0.5*buffer[index] + 1.0));

	} else {

	  this->rbegin()->set2ndDerivative(0.0);
	}

	{
	  reverse_iterator j = this->rbegin();
	  reverse_iterator i = j++;
	  
	  index = numberOfElements - 2;
	  
	  for ( ; j != this->rend(); ++i, ++j, --index)
	    j->set2ndDerivative(j->get2ndDerivative() + i->get2ndDerivative() * buffer[index]);
	}   

      } else {

	for (iterator i = this->begin(); i != this->end(); ++i)
	  i->set2ndDerivative(0.0);
      }
    }
  };


  /**
   * Template definition for spline interpolation method.
   */
  template<class JElement_t, template<class> class JCollection_t, class JResult_t = typename JElement_t::mapped_type>
  class JSplineFunction1D;


  /**
   * Template specialisation for spline interpolation method with returning mapped value type of container element.
   * This class implements the JFunction1D<> interface.
   */
  template<class JElement_t, template<class> class JCollection_t>
  class JSplineFunction1D<JElement_t, JCollection_t, typename JElement_t::mapped_type> :
    public JSpline<JElement_t, JCollection_t>,
    public JFunction1D<typename JElement_t::key_type, typename JElement_t::mapped_type>
  {
  public:

    typedef typename JElement_t::key_type                                        key_type;
    typedef typename JElement_t::mapped_type                                     mapped_type;
    typedef JElement_t                                                           value_type;

    typedef JCollection_t<JElement_t>                                            JContainer_t;

    typedef typename JContainer_t::const_iterator                                const_iterator;
    typedef typename JContainer_t::const_reverse_iterator                        const_reverse_iterator;

    typedef typename JContainer_t::iterator                                      iterator;
    typedef typename JContainer_t::reverse_iterator                              reverse_iterator;

    typedef JSpline<JElement_t, JCollection_t>                                   JSpline_t;

    typedef JFunction1D<key_type, mapped_type>                                   JFunction1D_t;           
    typedef typename JFunction1D_t::argument_type                                argument_type;
    typedef typename JFunction1D_t::result_type                                  result_type;
    typedef typename JFunction<key_type, result_type>::JExceptionHandler         JExceptionHandler_t;


    /**
     * Default constructor.			
     */
    JSplineFunction1D() :
      JSpline_t(),
      JFunction1D_t()
    {}


    /**
     * Cubic spline interpolation method.
     *
     * \param  x      argument value
     * \return        function value
     */
    mapped_type operator()(const key_type& x) const 
    {
      static const key_type eps = std::numeric_limits<key_type>::min();

      if (this->empty())
	return this->getExceptionHandler().action(JEmptyCollection("JSplineFunction1D<>::operator() no data."));

      const_iterator i = this->lower_bound(x);

      if ((i == this->begin() && this->getComparator()(x + eps, *(i++))) ||
	  (i == this->end()   && this->getComparator()(*(--i), x - eps)) ) {

	return this->getExceptionHandler().action(JValueOutOfRange("JSplineFunction1D<>::operator() x out of range."));
      }

      const_iterator j = i--;

      const double dx =  j->getX() - i->getX();
      const double a  = (j->getX() - x) / dx;
      const double b  = 1.0 - a;

      mapped_type y = a * i->getY() + b * j->getY()
        - a*b * ((a+1.0)*i->get2ndDerivative() + (b+1.0)*j->get2ndDerivative()) * dx*dx/6;

      return y;
    }
  };

  
  /**
   * Template specialisation for spline interpolation method with returning JFunctionValuePDF data structure.
   * This class implements the JFunction1D<> interface.
   */
  template<class JElement_t, template<class> class JCollection_t>
  class JSplineFunction1D<JElement_t, JCollection_t, JFunctionValuePDF<typename JElement_t::mapped_type> > :
    public JSpline<JElement_t, JCollection_t>, 
    public JFunction1D<typename JElement_t::key_type, JFunctionValuePDF<typename JElement_t::mapped_type> >
  {
  public:

    typedef typename JElement_t::key_type                                        key_type;
    typedef typename JElement_t::mapped_type                                     mapped_type;
    typedef JElement_t                                                           value_type;

    typedef JCollection_t<JElement_t>                                            JContainer_t;

    typedef typename JContainer_t::const_iterator                                const_iterator;
    typedef typename JContainer_t::const_reverse_iterator                        const_reverse_iterator;

    typedef typename JContainer_t::iterator                                      iterator;
    typedef typename JContainer_t::reverse_iterator                              reverse_iterator;

    typedef JSpline<JElement_t, JCollection_t>                                   JSpline_t;

    typedef JFunction1D<key_type, JFunctionValuePDF<mapped_type> >               JFunction1D_t;           
    typedef typename JFunction1D_t::argument_type                                argument_type;
    typedef typename JFunction1D_t::result_type                                  result_type;
    typedef typename JFunction<key_type, result_type>::JExceptionHandler         JExceptionHandler_t;


    /**
     * Default constructor.			
     */
    JSplineFunction1D() :
      JSpline_t(),
      JFunction1D_t()
    {}


    /**
     * Determination of second derivatives with no bounds.
     */
    virtual void compile()
    {
      compile(JSplineBounds<mapped_type>());
    }


    /**
     * Determination of second derivatives and integral values with specified bounds.
     *
     * \param bounds    first derivatives at two extrema.
     */
    void compile(const JSplineBounds<mapped_type>& bounds) 
    {
      if (!this->empty()) {

	JSpline_t::compile(bounds);
 
	this->begin()->setIntegral(0.0);

	for (iterator j = this->begin(), i = j++; j != this->end(); ++i, ++j) {
	  
	  const double dx = j->getX() - i->getX();
	  const double y  = i->getY() + j->getY();
	  //const double z  = i->get2ndDerivative() + j->get2ndDerivative();
	
	  const double v = dx * 0.50 * y;
	  /*
	  const double w = dx * 0.25 * z*dx*dx/6;

	  if ( (v > 0.0 && w < v) ||
	       (v < 0.0 && w > v) )
	    j->setIntegral(i->getIntegral() + v - w);
	  else
	    j->setIntegral(i->getIntegral() + v);
	  */
	  j->setIntegral(i->getIntegral() + v);
	}
      }
    }


    /**
     * Cubic spline interpolation method.
     *
     * \param  x      argument value
     * \return        function value
     */
    result_type operator()(const key_type& x) const 
    {
      static const key_type eps = std::numeric_limits<key_type>::min();

      if (this->empty())
	return this->getExceptionHandler().action(JEmptyCollection("JSplineFunction1D<>::operator() no data."));

      result_type result;

      const_iterator i = this->lower_bound(x);

      if        (i == this->begin() && this->getComparator()(x + eps, *(i++))) {

	result   = this->getExceptionHandler().action(JValueOutOfRange("JSplineFunction1D<>::operator() x out of range."));
	result.z = 0.0;
	result.v = 0.0;
	result.V = this->rbegin()->getIntegral();

	return result;

      } else if (i == this->end()   && this->getComparator()(*(--i), x - eps)) {

	result   = this->getExceptionHandler().action(JValueOutOfRange("JSplineFunction1D<>::operator() x out of range."));
	result.z = 0.0;
	result.v = this->rbegin()->getIntegral();
	result.V = this->rbegin()->getIntegral();

	return result;
      }

      const_iterator j = i--;

      const double dx =  j->getX() - i->getX();
      const double a  = (j->getX() - x) / dx;
      const double b  = 1.0 - a;

      result.y  = a * i->getY() + b * j->getY()
	- a*b * ((a+1.0)*i->get2ndDerivative() + (b+1.0)*j->get2ndDerivative()) * dx*dx/6;

      result.z  = (j->getY() - i->getY() + (i->get2ndDerivative()*(1.0 - 3*a*a) -
					    j->get2ndDerivative()*(1.0 - 3*b*b)) * dx*dx/6) / dx;
      
      result.v  = i->getIntegral()
	+ 0.5*dx * (i->getY() - 0.5*i->get2ndDerivative()*dx*dx/6)
	- 0.5*dx * ((a*a*i->getY() - b*b*j->getY()) +
		    (i->get2ndDerivative() * a*a*(0.5*a*a - 1.0) -
		     j->get2ndDerivative() * b*b*(0.5*b*b - 1.0)) * dx*dx/6);

      result.V  = this->rbegin()->getIntegral();

      return result;
    }
  };


  /**
   * Template specialisation for spline interpolation method with result type JFunctionValueHesse.
   * This class implements the JFunction1D<> interface.
   */
  template<class JElement_t, template<class> class JCollection_t>
  class JSplineFunction1D<JElement_t, JCollection_t, JFunctionValueHesse<typename JElement_t::mapped_type> > :
    public JSpline<JElement_t, JCollection_t>,
    public JFunction1D<typename JElement_t::key_type, JFunctionValueHesse<typename JElement_t::mapped_type> >
  {
  public:

    typedef typename JElement_t::key_type                                        key_type;
    typedef typename JElement_t::mapped_type                                     mapped_type;
    typedef JElement_t                                                           value_type;

    typedef JCollection_t<JElement_t>                                            JContainer_t;

    typedef typename JContainer_t::const_iterator                                const_iterator;
    typedef typename JContainer_t::const_reverse_iterator                        const_reverse_iterator;

    typedef typename JContainer_t::iterator                                      iterator;
    typedef typename JContainer_t::reverse_iterator                              reverse_iterator;

    typedef JSpline<JElement_t, JCollection_t>                                   JSpline_t;

    typedef JFunction1D<key_type, JFunctionValueHesse<mapped_type> >             JFunction1D_t;
    typedef typename JFunction1D_t::argument_type                                argument_type;
    typedef typename JFunction1D_t::result_type                                  result_type;
    typedef typename JFunction<key_type, result_type>::JExceptionHandler         JExceptionHandler_t;


    /**
     * Default constructor.
     */
    JSplineFunction1D() :
      JSpline_t(),
      JFunction1D_t()
    {}


    /**
     * Cubic spline interpolation method.
     *
     * \param  x      argument value
     * \return        function value
     */
    result_type operator()(const key_type& x) const
    {
      static const key_type eps = std::numeric_limits<key_type>::min();

      if (this->empty())
        return this->getExceptionHandler().action(JEmptyCollection("JSplineFunction1D<>::operator() no data."));

      const_iterator i = this->lower_bound(x);

      result_type result;

      if        (i == this->begin() && this->getComparator()(x + eps, *(i++))) {

        result     = this->getExceptionHandler().action(JValueOutOfRange("JSplineFunction1D<>::operator() x out of range."));
        result.fp  = 0.0;
        result.fpp = 0.0;

        return result;

      } else if (i == this->end()   && this->getComparator()(*(--i), x - eps)) {

        result     = this->getExceptionHandler().action(JValueOutOfRange("JSplineFunction1D<>::operator() x out of range."));
        result.fp  = 0.0;
        result.fpp = 0.0;

        return result;
      }

      const_iterator j = i--;

      const double dx =  j->getX() - i->getX();
      const double a  = (j->getX() - x) / dx;
      const double b  = 1.0 - a;

      result.f   = a * i->getY() + b * j->getY()
	- a*b * ((a+1.0)*i->get2ndDerivative() + (b+1.0)*j->get2ndDerivative()) * dx*dx/6;

      result.fp  = (j->getY() - i->getY() + (i->get2ndDerivative()*(1.0 - 3*a*a) -
					     j->get2ndDerivative()*(1.0 - 3*b*b)) * dx*dx/6) / dx;

      result.fpp = a * i->get2ndDerivative() + b * j->get2ndDerivative();

      return result;
    }
  };


  /**
   * Template definition of a functional map.
   */
  template<class JKey_t, class JValue_t, template<class, class> class JMap_t, class JResult_t>
  class JSplineMap;


  /**
   * Template specialisation of a functional map.
   * This class extends the map class and implements the JFunction<> interface.
   * The function evaluation is based on a spline interpolation.
   */
  template<class JFunction_t, template<class, class> class JMap_t, class JResult_t>
  class JSplineMap<typename JFunction_t::argument_type, JFunction_t, JMap_t, JResult_t> :
    public JFunctionalMap<typename JFunction_t::argument_type, JFunction_t, JMap_t, JResult_t>
  {
  public:

    typedef JFunctionalMap<typename JFunction_t::argument_type, 
			   JFunction_t, 
			   JMap_t, 
			   JResult_t>                            JFunctionalMap_t;

    typedef typename JFunctionalMap_t::argument_type             argument_type;
    typedef typename JFunctionalMap_t::result_type               result_type;

    typedef typename JFunctionalMap_t::const_iterator            const_iterator;
    typedef typename JFunctionalMap_t::const_reverse_iterator    const_reverse_iterator;

    typedef typename JFunctionalMap_t::iterator                  iterator;
    typedef typename JFunctionalMap_t::reverse_iterator          reverse_iterator;

    typedef typename JFunctionalMap_t::key_type                  key_type;
    typedef typename JFunctionalMap_t::mapped_type               mapped_type;

    typedef JSplineElement2D<argument_type, typename JFunction_t::result_type>  JElement_t;
    typedef JSplineFunction1D<JElement_t, JVector, result_type>                 JSpline_t;


    /**
     * Default constructor.
     */
    JSplineMap() : 
      JFunctionalMap_t()
    {}


    /**
     * Function compilation.
     */
    void compile()
    {
      for (iterator i = this->begin(); i != this->end(); ++i) {

        i->second.compile();
	
	buffer.insert(JElement_t(i->first));
      }
    }


    /**
     * Recursive interpolation method implementation.
     *
     * \param  pX     pointer to argument list
     * \return        function value
     */
    result_type evaluate(const argument_type* pX) const 
    {
      const double x = *pX;

      ++pX;  // next argument value

      const_iterator               p = this->begin();    // input
      typename JSpline_t::iterator q = buffer.begin();   // output

      for ( ; q != buffer.end(); ++q, ++p)
	q->setValue(p->second.evaluate(pX));

      buffer.compile();

      return buffer(x);
    }

  protected:
    mutable JSpline_t buffer;
  };


  /**
   * Helper methods to convert spline result to function value.
   */
  template<class T>
  inline double get_value(const JFunctionValuePDF  <T>& value) { return value.y; }

  template<class T>
  inline double get_value(const JFunctionValueHesse<T>& value) { return value.f; }
}

#endif