#ifndef __JPOLINT__
#define __JPOLINT__

#include <cmath>

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


namespace JTOOLS {

  namespace {
    using JLANG::JNumericalPrecision;
    using JLANG::JEmptyCollection;
    using JLANG::JValueOutOfRange;
  }


  /**
   * Template implementation for polynomial interpolation method.
   *
   * Polynomial 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<unsigned int N, class JElement_t, template<class> class JCollection_t>
  class JPolintFunction1D :
    public JCollection_t<JElement_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 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<argument_type, result_type>::JExceptionHandler    JExceptionHandler_t;



    /**
     * Default constructor.			
     */
    JPolintFunction1D() :
      JCollection_t<JElement_t>(),
      JFunction1D_t()
    {}


    /**
     * Function compilation.
     */
    virtual void compile() 
    {}


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

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

      const_iterator p = this->lower_bound(x);

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

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


      int n = std::min(N + 1, (unsigned int) this->size());  // number of points to interpolate
      

      for (int i = n/2; i != 0 && p != this->end();   --i, ++p) {}
      for (int i = n  ; i != 0 && p != this->begin(); --i, --p) {}


      argument_type d = std::numeric_limits<argument_type>::max();
      int           j = 0;

      for (int i = 0; i != n; ++p, ++i) {
	
	u[i] = p->getX() - x;
	v[i] = p->getY();
	w[i] = p->getY();

	if (fabs(u[i]) < d) {
	  d = u[i];
	  j = i;
	}
      }

      
      result_type y = v[j--];

      for (int m = 1; m != n; ++m) {

	for (int i = 0; i != n-m; ++i) {

	  const argument_type ho = u[ i ];
	  const argument_type hp = u[i+m];

	  d = ho - hp;

	  if (d == 0.0)
	    throw JNumericalPrecision("JPolint<>::operator()");

	  d = (v[i+1] - w[i]) / d;

	  v[i] = ho*d;
	  w[i] = hp*d;
	}

	if (2 *(j+1) < n - m)
	  y += v[j+1];
	else
	  y += w[j--];
      }

      return y;
    }

  private: 
    mutable argument_type u[N+1];
    mutable result_type   v[N+1];
    mutable result_type   w[N+1];
  };


  /**
   * Template specialisation for zero degree polynomial interpolation method.
   * The interpolation is based on a simple lookup table.
   */
  template<class JElement_t, template<class> class JCollection_t>
  class JPolintFunction1D<0, JElement_t, JCollection_t> :
    public JCollection_t<JElement_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 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<argument_type, result_type>::JExceptionHandler    JExceptionHandler_t;


    /**
     * Default constructor.			
     */
    JPolintFunction1D() :
      JCollection_t<JElement_t>(),
      JFunction1D_t()
    {}


    /**
     * Function compilation.
     */
    virtual void compile()
    {}


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

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

      const_iterator p = this->lower_bound(x);

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

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

      const_iterator q = p--;

      if (q == this->begin() || q->getX() - x < x - p->getX()) 
	return q->getY();
      else
	return p->getY();
    }
  };


  /**
   * Template specialisation for first degree polynomial interpolation method.
   * The linear interpolation method is identical to to the general template but is optimised in terms of speed.
   */
  template<class JElement_t, template<class> class JCollection_t>
  class JPolintFunction1D<1, JElement_t, JCollection_t> :
    public JCollection_t<JElement_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 JFunction1D<key_type, mapped_type>                                   JFunction1D_t;           
    typedef typename JFunction1D<key_type, mapped_type>::JExceptionHandler       JExceptionHandler_t;
    typedef typename JFunction1D_t::argument_type                                argument_type;
    typedef typename JFunction1D_t::result_type                                  result_type;


    /**
     * Default constructor.			
     */
    JPolintFunction1D() :
      JCollection_t<JElement_t>(),
      JFunction1D_t()
    {}


    /**
     * Function compilation.
     */
    virtual void compile()
    {}


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

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

      const_iterator p = this->lower_bound(x);

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

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

      const_iterator q = p--;

      const argument_type dx =  q->getX() - p->getX();
      const argument_type a  = (q->getX() - x) / dx;
      const argument_type b  = 1.0 - a;

      return a * p->getY() + b * q->getY();
    }
  };
 

  /**
   * Template definition of a functional map with polynomial interpolation.
   */
  template<unsigned int N, class JKey_t, class JValue_t, template<class, class> class JMap_t>
  class JPolintMap;


  /**
   * Template specialisation of a functional map with polynomial interpolation.
   * This class extends the map class and implements the JFunction<> interface.
   * The function evaluation is based on a polynomial interpolation.
   */
  template<unsigned int N, class JFunction_t, template<class, class> class JMap_t>
  class JPolintMap<N, typename JFunction_t::argument_type, JFunction_t, JMap_t> :
    public JFunctionalMap<typename JFunction_t::argument_type, JFunction_t, JMap_t, typename JFunction_t::result_type> 
  {
  public:
    
    typedef JFunctionalMap<typename JFunction_t::argument_type, 
			   JFunction_t, 
			   JMap_t, 
			   typename JFunction_t::result_type>    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 typename JFunctionalMap_t::value_type                value_type;


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


    /**
     * Recursive interpolation method implementation.
     *
     * \param  pX     pointer to argument list
     * \return        function value
     */
    result_type evaluate(const argument_type* pX) const 
    {
      static const key_type eps = std::numeric_limits<key_type>::min();

      const double x = *pX;

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

      const_iterator p = this->lower_bound(x);

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

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

      ++pX;  // next argument value


      int n = std::min(N + 1, (unsigned int) this->size());  // number of points to interpolate
      

      for (int i = n/2; i != 0 && p != this->end();   --i, ++p) {}
      for (int i = n  ; i != 0 && p != this->begin(); --i, --p) {}


      argument_type d = std::numeric_limits<argument_type>::max();
      int           j = 0;

      for (int i = 0; i != n; ++p, ++i) {
	
	u[i] = p->first - x;
	v[i] = p->second.evaluate(pX);
	w[i] = v[i];

	if (fabs(u[i]) < d) {
	  d = u[i];
	  j = i;
	}
      }


      result_type y = v[j--];

      for (int m = 1; m != n; ++m) {

	for (int i = 0; i != n-m; ++i) {

	  const argument_type ho = u[ i ];
	  const argument_type hp = u[i+m];

	  d = ho - hp;

	  if (d == 0.0)
	    throw JNumericalPrecision("JPolintMap<>::evaluate()");

	  const result_type z = (v[i+1] - w[i]) / d;

	  v[i] = ho*z;
	  w[i] = hp*z;
	}

	if (2*(j+1) <  n-m)
	  y += v[j+1];
	else
	  y += w[j--];
      }

      return y;
    }


  private: 
    mutable argument_type u[N+1];
    mutable result_type   v[N+1];
    mutable result_type   w[N+1];
  };


  /**
   * Template specialisation of a functional map with zero degree polynomial interpolation.
   * This class extends the map class and implements the JFunction<> interface.
   * The interpolation is based on a simple lookup table.
   */
  template<class JFunction_t, template<class, class> class JMap_t>
  class JPolintMap<0, typename JFunction_t::argument_type, JFunction_t, JMap_t> :
    public JFunctionalMap<typename JFunction_t::argument_type, JFunction_t, JMap_t, typename JFunction_t::result_type> 
  {
  public:
    
    typedef JFunctionalMap<typename JFunction_t::argument_type, 
			   JFunction_t, 
			   JMap_t, 
			   typename JFunction_t::result_type>    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 typename JFunctionalMap_t::value_type                value_type;


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


    /**
     * Recursive interpolation method implementation.
     *
     * \param  pX     pointer to argument list
     * \return        function value
     */
    result_type evaluate(const argument_type* pX) const 
    {
      static const key_type eps = std::numeric_limits<key_type>::min();

      const double x = *pX;

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

      const_iterator p = this->lower_bound(x);

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

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

      ++pX;  // next argument value

      const_iterator q = p--;

      if (q == this->begin() || q->first - x < x - p->first)
        return q->second.evaluate(pX);
      else
        return p->second.evaluate(pX);
    }
  };


  /**
   * Template specialisation of a functional map with first degree polynomial interpolation.
   * This class extends the map class and implements the JFunction<> interface.
   * The linear interpolation method is identical to to the general template but is optimised in terms of speed.
   */
  template<class JFunction_t, template<class, class> class JMap_t>
  class JPolintMap<1, typename JFunction_t::argument_type, JFunction_t, JMap_t> :
    public JFunctionalMap<typename JFunction_t::argument_type, JFunction_t, JMap_t, typename JFunction_t::result_type> 
  {
  public:
    
    typedef JFunctionalMap<typename JFunction_t::argument_type, 
			   JFunction_t, 
			   JMap_t, 
			   typename JFunction_t::result_type>    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 typename JFunctionalMap_t::value_type                value_type;


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


    /**
     * Recursive interpolation method implementation.
     *
     * \param  pX     pointer to argument list
     * \return        function value
     */
    result_type evaluate(const argument_type* pX) const 
    {
      static const key_type eps = std::numeric_limits<key_type>::min();

      const double x = *pX;

      if (this->empty() || this->size() == 1)
        return this->getExceptionHandler().action(JEmptyCollection("JFunctionalMap<>::evaluate() no data."));

      const_iterator p = this->lower_bound(x);
      
      if ((p == this->begin() && this->getComparator()(x + eps, *(p++))) ||
	  (p == this->end()   && this->getComparator()(*(--p), x - eps)) ) {
	
	return this->getExceptionHandler().action(JValueOutOfRange("JFunctionalMap<>::evaluate() x out of range."));
      }
      
      ++pX;  // next argument value

      const_iterator q = p--;
      
      const argument_type dx =  q->first - p->first;
      const argument_type a  = (q->first - x) / dx;
      const argument_type b  = 1.0 - a;
      
      return a * p->second.evaluate(pX) + b * q->second.evaluate(pX);
    }
  };
}

#endif