#ifndef __JFIT__JLEGENDREESTIMATOR__
#define __JFIT__JLEGENDREESTIMATOR__

#include "JMath/JLegendre.hh"
#include "JMath/JZero.hh"
#include "JMath/JMatrixNS.hh"
#include "JFit/JEstimator.hh"


/**
 * \file
 * Linear fit of JMATH::JLegendre.
 * \author mdejong
 */
namespace JFIT {}
namespace JPP { using namespace JFIT; }

namespace JFIT {

  using JMATH::JLegendre;

  /**
   * Linear fit of Legendre polynomial.
   */
  template<class JOrdinate_t, size_t N>
  class JEstimator< JLegendre<JOrdinate_t, N> > :
    public JLegendre<JOrdinate_t, N>
  {
  public:
    typedef JOrdinate_t                    ordinate_type;


    /**
     * Constructor.
     *
     * The template argument <tt>T</tt> refers to an iterator of a data structure which should have the following member methods:
     *   - double        %getX();       //
     *   - ordinate_type %getY();       //
     * where <tt>ordinate_type</tt> refers to the first template atgument of class JLegendre.
     *
     * \param  __begin          begin of data
     * \param  __end            end   of data
     */
    template<class T>
    JEstimator(T __begin, T __end)
    {
      using namespace std;
      using namespace JPP;

      if (distance(__begin, __end) >= NUMBER_OF_PARAMETERS) {

	vector<ordinate_type> Y(NUMBER_OF_PARAMETERS, getZero<ordinate_type>());
	vector<double>        h(NUMBER_OF_PARAMETERS, getZero<double>());

	V.resize(NUMBER_OF_PARAMETERS);
	V.reset();

	T i = __begin;

	this->xmin = i->getX();  advance(i, distance(__begin, __end) - 1);
	this->xmax = i->getX();

	for (i = __begin; i != __end; ++i) {

	  const double z = this->getX(i->getX());

	  for (size_t n = 0; n <= N; ++n) {
	    h[n] = legendre(n, z);
	  }

	  for (size_t row = 0; row <= N; ++row) {

	    V(row, row) += h[row] * h[row];
	    Y[row]      += h[row] * i->getY();

	    for (size_t col = 0; col != row; ++col) {
	      V(row, col) += h[row] * h[col];
	    }
	  }
	}

	for (size_t row = 0; row <= N; ++row) {
	  for (size_t col = 0; col != row; ++col) {
	    V(col, row) = V(row, col);
	  }
	}

	V.invert();

	for (size_t row = 0; row <= N; ++row) {

	  (*this)[row] = zero;

	  for (size_t col = 0; col <= N; ++col) {
	    (*this)[row] += V(row, col) * Y[col];
	  }
	}

      } else {
        throw JValueOutOfRange("JEstimator<JLegendre>::JEstimator(): Not enough data points.");
      }
    }


    static const int NUMBER_OF_PARAMETERS  =  N + 1;  //!< number of parameters of fit
    JMATH::JMatrixNS V;                               //!< co-variance matrix of fit parameters
  };
}

#endif