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

#include "TRandom3.h"

#include "JLang/JException.hh"
#include "JMath/JPolynome.hh"
#include "JTools/JPolfit.hh"
#include "JTools/JPolint.hh"
#include "JTools/JElement.hh"
#include "JTools/JGridCollection.hh"
#include "JTools/JGrid.hh"
#include "JTools/JQuantile.hh"
#include "JTools/JToolsToolkit.hh"

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


/**
 * \file
 *
 * Example program to test 1D Legendre interpolation of a polynome.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  unsigned int  numberOfEvents;
  unsigned int  numberOfBins;
  JPolynome     f1;
  double        sigma;
  double        P;
  int           debug;

  try {

    JParser<> zap("Example program to test 1D Legendre interpolation of a polynome.");

    zap['n'] = make_field(numberOfEvents)   =  0;
    zap['N'] = make_field(numberOfBins)     = 21;
    zap['P'] = make_field(f1);
    zap['s'] = make_field(sigma)            = 0.0;
    zap['p'] = make_field(P)                = 0.0;
    zap['d'] = make_field(debug)            =  3;

    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }


  const int N = 7;                        // number of points for Legendre interpolation
  const int M = 3;                        // degree of Legendre polynome

  typedef JPolfitFunction1D<N, 
			    M, 
			    JElement2D<double, double>, 
			    JGridCollection, 
			    double>       JFunction1D_t;

  JFunction1D_t g1;

  JPolintFunction1D<M, 
		    JElement2D<double, double>, 
		    JGridCollection, 
		    double> h1;

  const double xmin = -10.0;
  const double xmax = +10.0;  

  const JGrid<double> grid(numberOfBins, xmin, xmax);

  for (int i = 0; i != grid.getSize(); ++i) {

    const double x = grid.getX(i);

    g1[x]= f1(x) + (sigma > 0.0 && gRandom->Rndm() < P ? gRandom->Gaus(0.0, sigma) : 0.0);
  }

  copy(g1, h1);

  g1.compile();
  h1.compile();

  for (JFunction1D_t::const_iterator i = g1.begin(); i != g1.end(); ++i) {
    DEBUG(FIXED(12,5) << i->getX() << ' ' << FIXED(12,5) << i->getY() << endl);
  }
  
    
  if (numberOfEvents != 0) {

    JQuantile Q[2];
    
    for (unsigned int i = 0; i != numberOfEvents; ++i) {
      
      const double x = gRandom->Uniform(xmin, xmax);
      const double y = f1(x);
      const double u = g1(x);
      const double v = h1(x);

      DEBUG(FIXED(12,5) << x << ' ' << FIXED(12,5) << y << ' ' << FIXED(12,5) << u << ' ' << FIXED(12,5) << v << endl);

      Q[0].put(u - y);
      Q[1].put(v - y);
    }

    Q[0].print(cout);
    Q[1].print(cout, false);
  }

  return 0;
}