#include <string>
#include <iostream>
#include <iomanip>

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TMath.h"
#include "TString.h"

#include "JTools/JFunction1D_t.hh"
#include "JTools/JQuadrature.hh"
#include "JTools/JGrid.hh"

#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

  
// command line options
  
int     numberOfPoints;
bool    quadrature;
int     debug;


namespace {
  
  /**
   * Test function f1().
   *
   * \param  x       abscissa value
   * \return         ordinate value
   */
  inline Double_t f1(const Double_t x)
  {
    return TMath::Gaus(x, 0.0, 1.0, kTRUE);
  }


  /**
   * Derivative of f1().
   *
   * \param  x       abscissa value
   * \return         derivative value
   */
  inline Double_t g1(const Double_t x)
  {
    return f1(x) * -x;
  }


  /**
   * Integral of f1().
   *
   * \param  x       abscissa value
   * \return         integral value
   */
  inline Double_t G1(const Double_t x)
  {
    return 0.5 * (1.0 + TMath::Erf(sqrt(0.5)*x));
  }

  
  /**
   * Method to test interpolations for PDF.
   */
  template<class JFunction1D_t>
  void test(const char* title)
  {
    using namespace JPP;

    JFunction1D_t  buffer;

    const int    nx   = 1000;
    const double xmin = -3.5;
    const double xmax = +3.5;

    if (!quadrature)
      buffer.configure(make_grid(numberOfPoints, xmin, xmax), f1);
    else
      buffer.configure(JGaussHermite(numberOfPoints), f1);

    buffer.compile();
    
    TH1D f0(TString("f0[") +  title + "]", NULL, nx,  xmin, xmax);
    TH1D g0(TString("g0[") +  title + "]", NULL, nx,  xmin, xmax);
    TH1D G0(TString("G0[") +  title + "]", NULL, nx,  xmin, xmax);

    for (int i = 1; i <= f0.GetNbinsX(); ++i) {

      const Double_t x = f0.GetBinCenter(i);

      try {

	const typename JFunction1D_t::result_type result = buffer(x);
    
	f0.SetBinContent(i, result.f);
	g0.SetBinContent(i, result.fp);
	G0.SetBinContent(i, result.v);
      }
      catch(const std::exception& error) {
	//ERROR(error.what() << endl);
      }
    }

    f0.Write();
    g0.Write();
    G0.Write();
  }
}

/**
 * \file
 *
 * Example program to test interpolations for PDF.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  string  outputFile;
  
  try {
    
    JParser<> zap("Example program to test interpolations for PDF.");
    
    zap['o'] = make_field(outputFile);
    zap['N'] = make_field(numberOfPoints) = 15;
    zap['Q'] = make_field(quadrature);
    zap['d'] = make_field(debug)          = 3;
    
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }  

  
  TFile out(outputFile.c_str(), "recreate");

  const int    nx   = 1000;
  const double xmin = -3.5;
  const double xmax = +3.5;

  TH1D f0(TString("f0[") + "true" + "]", NULL, nx,  xmin, xmax);
  TH1D g0(TString("g0[") + "true" + "]", NULL, nx,  xmin, xmax);
  TH1D G0(TString("G0[") + "true" + "]", NULL, nx,  xmin, xmax);

  for (int i = 1; i <= f0.GetNbinsX(); ++i) {
      
    const Double_t x = f0.GetBinCenter(i);
    
    f0.SetBinContent(i, f1(x));
    g0.SetBinContent(i, g1(x));
    G0.SetBinContent(i, G1(x));
  }
  
  {
    const int N = 4;
    
    typedef JPolintFunction1D<N, JPolintElement2S<double, double>, JCollection, JResultPDF<double> >   JFunction1D_t;

    test<JFunction1D_t>("Polint");
  }

  {
    typedef JSplineFunction1S_t JFunction1D_t;
    
    test<JFunction1D_t>("Spline");
  }

  out.Write();
  out.Close();
}