JLegendreEstimator.hh

Go to the documentation of this file.

#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