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

#include "TMath.h"

#include "JTools/JCollection.hh"
#include "JTools/JGrid.hh"
#include "JTools/JQuantiles.hh"
#include "JTools/JFunction1D_t.hh"

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


/**
 * Function.
 *
 * \param           x
 * \return          function value
 */
inline Double_t g1(const Double_t x)
{
  return TMath::Gaus(x, 0.0, 1.0, kTRUE);
}


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


/**
 * \file
 * Example program to test JTOOLS::JQuantiles calculation of a function.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;

  int     numberOfBins;
  double  Q;
  double  precision;
  int     debug;

  try {

    JParser<> zap("Example program to test quantiles calculation of a function.");

    zap['N'] = make_field(numberOfBins)   = 13;
    zap['Q'] = make_field(Q)              = 0.5;
    zap['e'] = make_field(precision)      = 1.0e-3;
    zap['d'] = make_field(debug)          = 3;

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

  if (numberOfBins < 2) {
    FATAL("Fatal error: number of bins " << numberOfBins << endl); 
  }

  if (Q < 0.0 || Q > 1.0) {
    FATAL("Fatal error: quantile " << Q << endl); 
  }

  using namespace JPP;

  const Double_t xmin = -5.0;
  const Double_t xmax = +5.0;

  JQuantiles quantiles(make_grid(numberOfBins, xmin, xmax), g1, Q);

  NOTICE("quantity  " << setw(10)    << "calculated"                      << " / " << setw(10)    << "true"                 << endl);
  NOTICE("X         " << FIXED(10,6) << quantiles.getX()                  << " / " << FIXED(10,6) << 0.0                    << endl);
  NOTICE("FWHM      " << FIXED(10,6) << quantiles.getFWHM()               << " / " << FIXED(10,6) << 2.0*sqrt(2.0*log(2.0)) << endl);
  NOTICE("integral  " << FIXED(10,6) << quantiles.getIntegral()           << " / " << FIXED(10,6) << G1(xmax) - G1(xmin)    << endl);
  NOTICE("quantile  " << FIXED(10,6) << (G1(quantiles.getUpperLimit()) - 
					 G1(quantiles.getLowerLimit()))   << " / " << FIXED(10,6) << Q                      << endl);

  ASSERT(fabs(quantiles.getX()        - 0.0)                    < precision);
  ASSERT(fabs(quantiles.getFWHM()     - 2.0*sqrt(2.0*log(2.0))) < precision);
  ASSERT(fabs(quantiles.getIntegral() - (G1(xmax) - G1(xmin)))  < precision);
  ASSERT(fabs((G1(quantiles.getUpperLimit()) -
	       G1(quantiles.getLowerLimit())) - Q)              < precision);
}