#ifndef __JTOOLS__JPOLFIT__
#define __JTOOLS__JPOLFIT__

#include <cmath>
#include <iterator>
#include <algorithm>
#include <vector>
#include <utility>
#include <sstream>

#include "JLang/JException.hh"
#include "JLang/JAssert.hh"
#include "JLang/JStreamAvailable.hh"
#include "JMath/JMathSupportkit.hh"
#include "JMath/JLegendre.hh"
#include "JMath/JZero.hh"
#include "JTools/JFunctional.hh"
#include "JTools/JDistance.hh"


/**
 * \author mdejong
 */

namespace JTOOLS {}
namespace JPP { using namespace JTOOLS; }

namespace JTOOLS {

  using JLANG::JException;
  using JLANG::JFunctionalException;
  using JLANG::JNumericalPrecision;
  using JLANG::JValueOutOfRange;


  /**
   * Functional collection with Legendre polynomial interpolation.
   */
  template<unsigned int N,      // number of points to interpolate
	   unsigned int M,      // degree of polynomial
	   class JElement_t,
	   template<class, class> class JCollection_t,
	   class JResult_t,
	   class JDistance_t>
  class JPolfitFunction :
    public JCollection_t<JElement_t, JDistance_t>,
    public JFunction<typename JElement_t::abscissa_type,
		     typename JResultType<typename JElement_t::ordinate_type>::result_type>
  {
  public:

    typedef JCollection_t<JElement_t, JDistance_t>                               collection_type;

    typedef typename collection_type::abscissa_type                              abscissa_type;
    typedef typename collection_type::ordinate_type                              ordinate_type;
    typedef typename collection_type::value_type                                 value_type;
    typedef typename collection_type::distance_type                              distance_type;

    typedef typename collection_type::const_iterator                             const_iterator;
    typedef typename collection_type::const_reverse_iterator                     const_reverse_iterator;
    typedef typename collection_type::iterator                                   iterator;
    typedef typename collection_type::reverse_iterator                           reverse_iterator;

    typedef typename JResultType<ordinate_type>::result_type                     data_type;
    typedef JFunction<abscissa_type, data_type>                                  function_type;           

    typedef typename function_type::argument_type                                argument_type;
    typedef typename function_type::result_type                                  result_type;
    typedef typename function_type::JExceptionHandler                            exceptionhandler_type;


    /**
     * Default constructor.			
     */
    JPolfitFunction()
    {
      STATIC_CHECK((N > M));
    }


    /**
     * Recursive interpolation method implementation.
     *
     * \param  pX              pointer to abscissa values
     * \return                 function value
     */
    virtual result_type evaluate(const argument_type* pX) const override 
    {
      using namespace std;
      using namespace JPP;

      const argument_type x = *pX;

      if (this->size() <= 1u) {

	try {
	  return this->getExceptionHandler().action();
	}
	catch (const JException& error) {

	  std::ostringstream os;

	  os << __FILE__ << ':' << __LINE__ << " not enough data " << STREAM("?") << x;
	
	  throw JFunctionalException(os.str());
	}
      }

      const_iterator p = this->lower_bound(x);

      if ((p == this->begin() && this->getDistance(x, (p++)->getX()) > distance_type::precision) ||
          (p == this->end()   && this->getDistance((--p)->getX(), x) > distance_type::precision)) {

	try {
	  return this->getExceptionHandler().action();
	}
	catch (const JException& error) {

	  std::ostringstream os;

	  os << __FILE__ << ':' << __LINE__ << " abscissa out of range " 
	     << STREAM("?") << x                      << " <> " 
	     << STREAM("?") << this->begin() ->getX() << ' '
	     << STREAM("?") << this->rbegin()->getX();
	  
	  throw JValueOutOfRange(os.str());
	}
      }

      ++pX;  // next argument value


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

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

      buffer.clear();

      for (const_iterator q = p; n != 0; ++q, --n) {
	buffer.push_back(make_pair(this->getDistance(p->getX(), q->getX()), function_type::getValue(q->getY(), pX)));
      }

      const JLegendre<result_type, M> g1(buffer.begin(), buffer.end());

      return g1(this->getDistance(p->getX(), x));
    }

  private: 
    /**
     * Function compilation.
     */
    virtual void do_compile() override 
    {}

    mutable std::vector< std::pair<abscissa_type, result_type> > buffer;
  };


  /**
   * Template class for polynomial interpolation in 1D
   *
   * This class implements the JFunction1D interface.
   */
  template<unsigned int N,
	   unsigned int M,
	   class JElement_t,
	   template<class, class> class JCollection_t,
	   class JResult_t   = typename JElement_t::ordinate_type,
	   class JDistance_t = JDistance<typename JElement_t::abscissa_type> >
  class JPolfitFunction1D :
    public JPolfitFunction<N, M, JElement_t, JCollection_t, JResult_t, JDistance_t>,
    public JFunction1D<typename JElement_t::abscissa_type, JResult_t> 
  {
  public:

    typedef JCollection_t<JElement_t, JDistance_t>                               collection_type;

    typedef typename collection_type::abscissa_type                              abscissa_type;
    typedef typename collection_type::ordinate_type                              ordinate_type;
    typedef typename collection_type::value_type                                 value_type;
    typedef typename collection_type::distance_type                              distance_type;

    typedef typename collection_type::const_iterator                             const_iterator;
    typedef typename collection_type::const_reverse_iterator                     const_reverse_iterator;
    typedef typename collection_type::iterator                                   iterator;
    typedef typename collection_type::reverse_iterator                           reverse_iterator;

    typedef JFunction1D<abscissa_type, JResult_t>                                function_type;           

    typedef typename function_type::argument_type                                argument_type;
    typedef typename function_type::result_type                                  result_type;
    typedef typename function_type::JExceptionHandler                            exceptionhandler_type;


    /**
     * Default contructor.
     */
    JPolfitFunction1D()
    {}
  };
}

#endif