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

#include "TMath.h"
#include "TRandom.h"

#include "JTools/JFunction1D_t.hh"
#include "JTools/JTransformableMultiFunction.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JTools/JQuantile.hh"

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


namespace {

  using namespace JPP;

  /**
   * Test function.
   *
   * \param  x       x value
   * \param  y       y value
   * \return         function value
   */
  inline Double_t f2(const Double_t x,
		     const Double_t y)
  {
    return TMath::Gaus(y, 0.0, x, kTRUE);
  }


  /**
   * Derivative of test function.
   *
   * \param  x       x value
   * \param  y       y value
   * \return         derivative value
   */
  inline Double_t g2(const Double_t x,
		     const Double_t y)
  {
    return TMath::Gaus(y, 0.0, x, kTRUE) * -y/x;
  }


  /**
   * Integral of test function.
   *
   * \param  x       x value
   * \param  y       y value
   * \return         integral value
   */
  inline Double_t G2(const Double_t x,
		     const Double_t y)
  {
    return 0.5 * (1.0 + TMath::Erf(sqrt(0.5)*y/x));
  }


  typedef JTransformableMultiFunction<JSplineFunction1S_t, JMAPLIST<JPolint1FunctionalMap>::maplist>   function_type;
  
  typedef typename function_type::transformer_type   transformer_type;
  typedef typename function_type::result_type        result_type;


  /**
   * Transformer.
   */
  struct JTransformer_t :
    public transformer_type
  {
    typedef typename transformer_type::clone_type              clone_type;
    typedef typename transformer_type::argument_type           argument_type;
    typedef typename transformer_type::const_array_type        const_array_type;


    /**
     * Evaluate y as a function of x.
     *
     * \param  buffer               x value
     * \param  xn                   y value
     * \return                      y value
     */
    virtual argument_type putXn(const_array_type& buffer, const argument_type xn) const 
    {
      return xn * buffer[0];
    }


    /**
     * Evaluate y as a function of x.
     *
     * \param  buffer               x value
     * \param  xn                   y value
     * \return                      y value
     */
    virtual argument_type getXn(const_array_type& buffer, const argument_type xn) const 
    {
      return xn / buffer[0];
    }


    /**
     * Weight function.
     *
     * \param  buffer               x0 - xn-1 values
     * \return                      weight
     */
    virtual double getWeight(const_array_type& buffer) const 
    {
      return 1.0;
    }


    /**
     * Read from input.
     *
     * \param  in    reader
     * \return       reader
     */
    virtual JReader& read(JReader& in) 
    {
      return in;
    }


    /**
     * Write to output.
     *
     * \param  out   writer
     * \return       writer
     */
    virtual JWriter& write(JWriter& out) const 
    {
      return out;
    }


    /**
     * Get clone of this object.
     *
     * \return              pointer to newly created object
     */
    virtual clone_type clone() const 
    {
      return new JTransformer_t(*this);
    }
  };
}

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

  int     numberOfEvents;
  int     numberOfBins;
  bool    transform;
  int     debug;

  try {

    JParser<> zap("Example program to test transformable function.");

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

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

  if (numberOfEvents <= 0) {
    FATAL("No events." << endl);
  }


  const double xmin = 0.1;
  const double xmax = 1.0;
  const double dx   = (xmax - xmin) / 3;

  function_type h2;

  for (double x = xmin; x < xmax + 0.5*dx; x += dx) {

    const double ymin = -5.0 * xmax;
    const double ymax = +5.0 * xmax;
    const double dy   = (ymax - ymin) / (numberOfBins - 1);

    for (double y = ymin; y < ymax + 0.5*dy; y += dy) {
      h2[x][y] = f2(x,y);
    }
  }

  h2.setExceptionHandler(new typename function_type::JDefaultResult(JMATH::zero));

  if (transform) {
    h2.transform(JTransformer_t());
  }

  h2.compile();


  JQuantile Q[] = {
    JQuantile("function"),
    JQuantile("derivative"),
    JQuantile("integral")
  };
  
  const double ymin = -5.0 * xmax;
  const double ymax = +5.0 * xmax;

  for (int i = 0; i != numberOfEvents; ++i) {
      
    const double x = gRandom->Uniform(xmin, xmax);
    const double y = gRandom->Uniform(ymin, ymax);
    
    try {
      
      const result_type value = h2(x,y);

      const double f  = f2(x,y);
      const double fp = g2(x,y);
      const double v  = G2(x,y);

      Q[0].put(value.f  - f);
      Q[1].put(value.fp - fp);
      Q[2].put(value.v  - v);
    }
    catch(const std::exception& error) {}
  }

    
  if (debug >= notice_t) {

    cout << "     ";
    for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
      cout << "  " << setw(10) << Q[i].getTitle();
    }
    cout << endl;

    cout << "mean ";
    for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
      cout << "  " << SCIENTIFIC(10,2) << Q[i].getMean();
    }
    cout << endl;

    cout << "RMS  ";
    for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
      cout << "  " << SCIENTIFIC(10,2) << Q[i].getSTDev();
    }
    cout << endl;
  }
}