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

#include "TRandom3.h"

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

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

namespace {

  /**
   * Test method for interpolation with derivatives.
   *
   * Accuracy of  nth  derivative of  N  degree polynome
   * is equal to n-jth derivative of N-j degree polynome.
   *
   * \param  argc    number of command line arguments
   * \param  argv    list   of command line arguments
   * \return         status
   */
  template<int N>
  int do_main(int argc, char **argv)
  {
    using namespace std;
    using namespace JPP;

    unsigned int  numberOfEvents;
    unsigned int  numberOfBins;
    int           debug;

    try {

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

      zap['n'] = make_field(numberOfEvents);
      zap['N'] = make_field(numberOfBins)     = 21;
      zap['d'] = make_field(debug)            =  3;

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

    if (N < 1) {
      FATAL("Number of degrees " << N << " < 1." << endl);
    }

    
    using namespace JPP;

    // Analytical function and its derivatives.
    
    vector<JTrigonometric> f1(1, JTrigonometric(sin));

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


    typedef JResultPolynome<N, double>                                                JResult_t;
    typedef JPolintFunction1D<N, JElement2D<double, double>, JCollection, JResult_t>  JFunction1D_t;

    JFunction1D_t  polint;

    const double xmin = 0.0;
    const double xmax = PI;
 
    //const double dx   = ((N+1)/2) * (xmax - xmin) / (double) (numberOfBins - 2*((N+1)/2));
    const double dx   = 0.0 * PI;
    
    const JGrid<double> grid(numberOfBins + 1, xmin - dx, xmax + dx);

    polint.configure(grid, f1[0]);
  
    polint.compile();
  
    polint.setExceptionHandler(new typename JFunction1D_t::JDefaultResult(JMATH::zero));


    if (numberOfEvents != 0) {

      JQuantile Q[N+1];

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

	ostringstream os;

	os << "f^" << i;

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

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

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

	for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
	  Q[i].put(f1[i](x) - result.y[i]);
	}
      }

      cout << "\nInterpolation with " << N << " degrees polynomial:" << endl;

      for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
	Q[i].print(cout, i == 0);
      }
    }
    
    return 0;
  }
}

/**
 * \file
 *
 * Example program to test 1D interpolation of a trigonometric function.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  if (do_main<2>(argc, argv) != 0) { return 1; }
  if (do_main<3>(argc, argv) != 0) { return 1; }
  if (do_main<4>(argc, argv) != 0) { return 1; }
  if (do_main<5>(argc, argv) != 0) { return 1; }
}