#ifndef _utl_QuadraticFitter_h_
#define _utl_QuadraticFitter_h_

#include <cmath>

#include <utl/Math.h>
#include <utl/Histogram.h>
#include <utl/AugerException.h>
#include <utl/QuadraticFitData.h>


namespace utl {


  class PoissonianErrors {
  public:
    static double GetError(const double& value)
    { return std::sqrt(value); }

    static double GetError2(const double& value)
    { return value; }
  };


  class FlatErrors {
  public:
    static double GetError(const double& /*value*/) { return 1; }
    static double GetError2(const double& /*value*/) { return 1; }
  };


  /**
    \class QuadraticFitter QuadraticFitter.h "utl/QuadraticFitter.h"

    \author Darko Veberic
    \date 27 August 2007
    \version $Id: QuadraticFitter.h 14717 2009-09-17 20:24:36Z lukas $
  */

  template<class Histogram, class ErrorPolicy /*= utl::PoissonianErrors*/>
  class QuadraticFitter {
  public:
    QuadraticFitter(const Histogram& h)
      : fStartBin(0), fStopBin(h.GetNBins()), fHistogram(h) { }

    QuadraticFitter(const Histogram& h,
                    const int startBin, const int stopBin)
      : fStartBin(startBin), fStopBin(stopBin), fHistogram(h)
    {
      if (startBin < 0 || stopBin - startBin < 3 || stopBin >= int(h.GetNBins()))
        throw utl::OutOfBoundException("start, stop bin out of bounds of histogram");
    }

    QuadraticFitter(const Histogram& h,
                    const double xStart, const double xStop)
      : fHistogram(h)
    {
      int start = h.GetBinIndex(xStart);
      if (start == Histogram::eUnderflow)
        start = 0;
      fStartBin = start;
      int stop = h.GetBinIndex(xStop);
      if (stop == Histogram::eOverflow)
        stop = h.GetNBins() - 1;
      fStopBin = stop;
      if (start < 0 || stop < 0 || fStopBin - fStartBin < 3)
        throw utl::OutOfBoundException("start, stop bin out of bounds of histogram");
    }

    double
    GetExtremePosition()
    {
      MakeSums();
      MakeCoeficientTimesDeterminant(1);
      MakeCoeficientTimesDeterminant(2);
      return - fK[1] / (2 * fK[2]);
    }

    double
    GetExtremePosition(double& posError, int& ndof)
    {
      const double pos = GetExtremePosition();

      const double sigma2k1 = Sqr(fx2Sum) - fx0Sum * fx4Sum; // /det
      const double sigma2k2 = Sqr(fx1Sum) - fx0Sum * fx2Sum; // /det
      const double covk1k2  = fx0Sum * fx3Sum - fx1Sum * fx2Sum; // /det
      const double k22 = Sqr(fK[2]); // /det^2

      MakeDeterminant();

      posError =
        0.5 * sqrt(fDet*(sigma2k1 + pos*covk1k2 + sigma2k2*Sqr(fK[1])/k22) / k22);

      ndof = fN - 3;

      return pos;
    }

    double
    GetExtremePosition(double& posError, double& chi2, int& ndof)
    {
      const double pos = GetExtremePosition(posError, ndof);

      MakeCoeficientTimesDeterminant(0);

      chi2 = ((Sqr(fK[0]) * fx0Sum + 2*fK[0]*fK[1] * fx1Sum + (Sqr(fK[1]) + 2*fK[0]*fK[2]) * fx2Sum
               + 2*fK[1]*fK[2] * fx3Sum + Sqr(fK[2]) * fx4Sum) / Sqr(fDet)
              - 2 * (fK[0] * fx0ySum + fK[1] * fx1ySum + fK[2] * fx2ySum) / fDet + fx0y2Sum);

      return pos;
    }

    double
    GetExtremePosition(double& posError, double& chi2, int& ndof, double polynomialCoeff[3])
    {
      const double pos = GetExtremePosition(posError, chi2, ndof);
      for (int i = 0; i < 3; ++i)
        polynomialCoeff[i] = fK[i] / fDet;
      return pos;
    }

    void
    GetFitData(QuadraticFitData& fit)
    {
      fit.fStart = fHistogram.GetBinCenter(fStartBin);
      fit.fStop = fHistogram.GetBinCenter(fStopBin);
      fit.fExtremePosition =
        GetExtremePosition(fit.fExtremePositionError, fit.fChi2, fit.fNdof, fit.fCoefficients);
    }

  private:
    void
    MakeSums()
    {
      fx0Sum = fx1Sum = fx2Sum = fx3Sum = fx4Sum =
        fx0ySum = fx1ySum = fx2ySum = fx0y2Sum = 0;

      for (unsigned int i = fStartBin; i < fStopBin; ++i) {
        const double x = fHistogram.GetBinCenter(i);
        const double raw = fHistogram.GetBin(i);
        const double invSigma2 = Sqr(fHistogram.GetBinSize(i)) / ErrorPolicy::GetError2(raw);
        const double y = fHistogram.GetBinAverage(i);
        fx0Sum += invSigma2;
        fx1Sum += x * invSigma2;
        fx0ySum += y * invSigma2;
        fx0y2Sum += Sqr(y) * invSigma2;
        const double xy = x*y;
        fx1ySum += xy * invSigma2;
        const double x2 = Sqr(x);
        fx2Sum += x2 * invSigma2;
        fx2ySum += x*xy * invSigma2;
        const double x3 = x2*x;
        fx3Sum += x3 * invSigma2;
        fx4Sum += x*x3 * invSigma2;
      }

      fN = fStopBin - fStartBin;
    }

    void
    MakeDeterminant()
    {
      fDet =   Sqr(fx2Sum) * fx2Sum + fx0Sum * Sqr(fx3Sum)
             + Sqr(fx1Sum) * fx4Sum
             - fx2Sum * (2 * fx1Sum * fx3Sum + fx0Sum * fx4Sum);
    }

    void
    MakeCoeficientTimesDeterminant(const int i)
    {
      switch (i) {
      case 0:
        fK[0] =   Sqr(fx3Sum) * fx0ySum + fx4Sum * (fx1Sum * fx1ySum - fx2Sum * fx0ySum)
                + Sqr(fx2Sum) * fx2ySum - fx3Sum * (fx2Sum * fx1ySum + fx1Sum * fx2ySum); // /det
      case 1:
        fK[1] =   fx4Sum * (fx1Sum * fx0ySum - fx0Sum * fx1ySum) + Sqr(fx2Sum) * fx1ySum
                + fx0Sum * fx3Sum * fx2ySum - fx2Sum * (fx3Sum * fx0ySum + fx1Sum * fx2ySum); // /det
      case 2:
        fK[2] =   Sqr(fx2Sum) * fx0ySum + fx3Sum * (fx0Sum * fx1ySum - fx1Sum * fx0ySum)
                + Sqr(fx1Sum) * fx2ySum - fx2Sum * (fx1Sum * fx1ySum + fx0Sum * fx2ySum); // /det
      }
    }

    unsigned int fStartBin;
    unsigned int fStopBin;
    const Histogram& fHistogram;

    double fx0Sum;
    double fx1Sum;
    double fx2Sum;
    double fx3Sum;
    double fx4Sum;
    double fx0ySum;
    double fx1ySum;
    double fx2ySum;
    double fx0y2Sum;
    int fN;
    double fDet;
    double fK[3];
  };


  // convinience makers
  template<class Histogram, class ErrorPolicy>
  inline
  QuadraticFitter<Histogram, ErrorPolicy>
  MakeQuadraticFitter(const Histogram& h, const ErrorPolicy /*ep*/)
  {
    return QuadraticFitter<Histogram, ErrorPolicy>(h);
  }

  template<class Histogram>
  inline
  QuadraticFitter<Histogram>
  MakeQuadraticFitter(const Histogram& h)
  {
    return QuadraticFitter<Histogram>(h);
  }


  template<class Histogram, typename Type, class ErrorPolicy>
  inline
  QuadraticFitter<Histogram, ErrorPolicy>
  MakeQuadraticFitter(const Histogram& h,
                      const Type xStart, const Type xStop,
                      const ErrorPolicy ep)
  {
    return QuadraticFitter<Histogram, ErrorPolicy>(h, xStart, xStop);
  }

  template<class Histogram, typename Type>
  inline
  QuadraticFitter<Histogram>
  MakeQuadraticFitter(const Histogram& h,
                      const Type xStart, const Type xStop)
  {
    return QuadraticFitter<Histogram>(h, xStart, xStop);
  }

}

#endif