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

#include "TRandom3.h"

#include "JMath/JPolynome.hh"
#include "JTools/JMapList.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JTools/JMultiFunction.hh"
#include "JTools/JQuantile.hh"

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


namespace {

  using namespace JPP;


  /**
   * 3D polynomial function.
   */
  struct JFunction3D {

    /**
     * Result.
     */
    struct result_type {
      /**
       * Constructor.
       *
       * \param  __f         function value  f
       * \param  __fx        derivative     df/dx
       * \param  __fy        derivative     df/dy
       * \param  __fz        derivative     df/dz
       */
      result_type(const double __f,
		  const double __fx,
		  const double __fy,
		  const double __fz) :
	f (__f),
	fx(__fx),
	fy(__fy),
	fz(__fz)
      {}


      /**
       * Type conversion operator.
       *
       * \return             function value
       */
      operator double() const
      { 
	return f;
      }
    

      const double f;
      const double fx;
      const double fy;
      const double fz;
    };
  

    /**
     * Constructor.
     *
     * In each dimension, the function value and the derivative are given by: 
     * <pre>
     *      f (x) = f1.getValue(x);
     *      f'(x) = f1.getDerivative(x);
     * </pre>
     *
     * \param  f1          polynome
     */
    JFunction3D(const JPolynome& f1) :
      __f1(f1)
    {}


    /**
     * Evaluation of function value and derivatives.
     *
     * \param  x           x-value
     * \param  y           y-value
     * \param  z           z-value
     * \return             function and derivatives
     */
    inline result_type operator()(const double x,
				  const double y,
				  const double z) const
    {
      return result_type(f (x) * f (y) * f (z),
			 fp(x) * f (y) * f (z),
			 f (x) * fp(y) * f (z),
			 f (x) * f (y) * fp(z));
    }


  protected:
    /**
     * Evaluation of function value.
     *
     * \param  x           x-value
     * \return             function value
     */
    inline double f(const double x) const
    { 
      return __f1.getValue(x);
    }


    /**
     * Evaluation of derivative.
     *
     * \param  x           x-value
     * \return             derivative
     */
    inline double fp(const double x) const
    { 
      return __f1.getDerivative(x);
    }

    JPolynome __f1;
  };


  /**
   * Test method for multi-dimensional interpolation with derivatives.
   *
   * \param  argc    number of command line arguments
   * \param  argv    list   of command line arguments
   * \param  title   title of this template process
   * \return         status
   */
  template<class JMultiFunction_t>
  int do_main(int argc, char **argv, const char* title)
  {
    using namespace std;

    int         numberOfEvents;
    int         numberOfBins;
    JPolynome   f1;
    int         debug;

    try {

      JParser<> zap("Example program to test multi-dimensional interpolation with derivatives.");

      zap['n'] = make_field(numberOfEvents)   = 1000;
      zap['N'] = make_field(numberOfBins)     = 11;
      zap['P'] = make_field(f1)               = JPolynome(1.0, 1.0, 1.0);
      zap['d'] = make_field(debug)            = 3;

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


    if (numberOfEvents <= 0) { FATAL("Number of events " << numberOfEvents << endl); }
    if (numberOfBins   <= 2) { FATAL("Number of bins   " << numberOfBins   << endl); }

    using namespace JPP;


    const JFunction3D f3(f1);

    const double xmin = -1.0;
    const double xmax = +1.0;
    const double dx   = (xmax - xmin) / (numberOfBins - 1);
   

    JMultiFunction_t g3;

    for (double x = xmin; x < xmax + 0.5*dx; x += dx) {
      for (double y = xmin; y < xmax + 0.5*dx; y += dx) {
	for (double z = xmin; z < xmax + 0.5*dx; z += dx) {
	  g3[x][y][z] = f3(x,y,z);
	}
      }
    }

    g3.compile();
    //g3.setExceptionHandler(new JFunction1D_t::JDefaultResult(JMATH::zero));


    JQuantile f ("f    ");
    JQuantile fx("df/dx");
    JQuantile fy("df/dy");
    JQuantile fz("df/dz");
    
    for (int i = 0; i != numberOfEvents; ++i) {
      
      const double x  = gRandom->Uniform(xmin, xmax);
      const double y  = gRandom->Uniform(xmin, xmax);
      const double z  = gRandom->Uniform(xmin, xmax);
      
      const          JFunction3D     ::result_type u = f3(x,y,z);
      const typename JMultiFunction_t::result_type v = g3(x,y,z);
      
      f .put(u.f   -  v.f .f .f);
      fx.put(u.fx  -  v.fp.f .f);
      fy.put(u.fy  -  v.f .fp.f);
      fz.put(u.fz  -  v.f .f .fp);
    }
    
    cout << endl << title << ":" << endl;

    for (const JQuantile* p : { &f, &fx, &fy, &fz }) {
      p->print(cout, p == &f);
    }

    return 0;
  }
}


/**
 * \file
 * Example program to test multi-dimensional interpolation with derivatives.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  {
    typedef JGridPolint3Function1H_t                          JFunction1D_t;
    typedef JMAPLIST<JPolint3FunctionalGridMapH,
		     JPolint3FunctionalGridMapH>::maplist     JMaplist_t;
    typedef JMultiFunction<JFunction1D_t, JMaplist_t>         JMultiFunction_t;

    if (do_main<JMultiFunction_t>(argc, argv, "Polint") != 0) return 1;
  }

  {
    typedef JGridSplineFunction1H_t                           JFunction1D_t;
    typedef JMAPLIST<JSplineFunctionalGridMapH,
		     JSplineFunctionalGridMapH>::maplist      JMaplist_t;
    typedef JMultiFunction<JFunction1D_t, JMaplist_t>         JMultiFunction_t;

    if (do_main<JMultiFunction_t>(argc, argv, "Spline") != 0) return 1;
  }
}