#include <string>
#include <iostream>
#include <sstream>
#include <iomanip>
#include <vector>

#include "TRandom3.h"

#include "JLang/JException.hh"
#include "JMath/JPolynome.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JGrid.hh"
#include "JTools/JQuantile.hh"

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


/**
 * \file
 *
 * Example program to test 1D 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        x0;
  int           debug;

  try {

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

    zap['n'] = make_field(numberOfEvents)   =  0;
    zap['N'] = make_field(numberOfBins)     = 21;
    zap['P'] = make_field(f1);
    zap['x'] = make_field(x0)               =  0.0;
    zap['d'] = make_field(debug)            =  1;

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

  if (f1.empty()) {
    FATAL("Invalid polynomial.");
  }


  const int N = 3;                      // degree of polynomial interpolation
  
  typedef JResultPolynome<N, double>                                                JResult_t;
  typedef JPolintFunction1D<N, JElement2D<double, double>, JCollection, JResult_t>  JFunction1D_t;


  JFunction1D_t  polint;


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

  polint.configure(JGrid<double>(numberOfBins, xmin, xmax), f1);
  
  polint.compile();
  
  polint.setExceptionHandler(new JFunction1D_t::JDefaultResult(JMATH::zero));
    

  vector<JPolynome> fp(1, f1);

  for (int i = 1; i != (int) f1.size(); ++i) {
    fp.push_back(fp.rbegin()->getDerivative());
  }

  DEBUG("polynome: ");

  for (int i = 0; i != (int) f1.size(); ++i) {
    DEBUG(' ' << fp[i](x0));
  }
  DEBUG(endl);

  DEBUG("result:   ");

  JResult_t result = polint(x0);

  for (int i = 0; i != N+1; ++i) {
    DEBUG(' ' << result.y[i]);
  }
  DEBUG(endl);


  if (numberOfEvents != 0) {

    JQuantile Q[N+1];

    for (int j = 0; j != sizeof(Q)/sizeof(Q[0]); ++j) {

      ostringstream os;

      os << "f^" << j;

      Q[j].setTitle(os.str());
    }

    const int M = min(fp.size(), sizeof(Q)/sizeof(Q[0]));

    for (unsigned int i = 0; i != numberOfEvents; ++i) {

      const double x   = gRandom->Uniform(xmin, xmax);
      
      JResult_t result = polint(x);

      for (int j = 0; j != M; ++j) {
	Q[j].put(fp[j](x) - result.y[j]);
      }
    }

    for (int j = 0; j != M; ++j) {
      Q[j].print(cout, j == 0);
    }
  }
}