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

#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JPhysics/JPDFTable.hh"


/**
 * Main program to apply TTS to interpolation tables for PDFs.
 */
int main(int argc, char **argv)
{
  using namespace std;

  string              inputFile;
  string              outputFile;
  int                 numberOfPoints;
  double              epsilon;
  double              TTS;
  int                 debug;

  try {

    JParser<> zap;

    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile);
    zap['n'] = make_field(numberOfPoints) = 25;
    zap['e'] = make_field(epsilon)        = 1.0e-10;
    zap['T'] = make_field(TTS)            = 0.0;         // [ns]
    zap['d'] = make_field(debug)          = 0;

    if (zap.read(argc, argv) != 0)
      return 1;
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }


  using namespace JPHYSICS;
  using namespace JTOOLS;

  typedef JSplineFunction1D_t                                     JFunction1D_t;
  typedef JMultipleMap<3, JPolint1FunctionalMap>::type_list       JMapList_t;
  typedef JPDFTable<JFunction1D_t, JMapList_t>                    JPDF_t;
  typedef JArray<3, JFunction1D_t::argument_type>                 JArray_t;

  JPDF_t pdf;


  try {

    NOTICE("loading input from file " << inputFile << "... " << flush);

    pdf.load(inputFile.c_str());

    NOTICE("OK" << endl);
  }
  catch(const JException& error) {
    FATAL(error.what() << endl);
  }

  JFunction1D_t::JSupervisor supervisor(new JFunction1D_t::JDefaultResult(0.0));
    
  for (JPDF_t::super_iterator i = pdf.super_begin(); i != pdf.super_end(); ++i)
    (*i).getValue().setExceptionHandler(supervisor);


  NOTICE("Gauss-Hermite integration... " << flush);

  const JPDF_t::transformer_type& transformer = *pdf.transformer;
  
  JGaussHermite engine(numberOfPoints, epsilon);
  
  for (JPDF_t::super_iterator i = pdf.super_begin(); i != pdf.super_end(); ++i) {
    
    const JArray_t array = (*i).getKey();
    JFunction1D_t& f1    = (*i).getValue();

    //const double   R     = array[0];
    //const double   theta = array[1];
    //const double   phi   = array[2];
    
    const JFunctionGetTransformer<3, JFunction1D_t::value_type> get(transformer, array);
    const JFunctionPutTransformer<3, JFunction1D_t::value_type> put(transformer, array);

    if (!f1.empty()) {
      
      f1.transform(get);
      f1.compile();

      set<double> X;
      
      for (JFunction1D_t::const_iterator j = f1.begin(); j != f1.end(); ++j)
	X.insert(j->getX());
      
      // add points at extrema

      const double tmin = f1. begin()->getX();
      const double tmax = f1.rbegin()->getX();
      
      for (JGaussHermite::const_iterator j = engine.begin(); j != engine.end(); ++j) {
	
	const double& u = j->getX();
	
	if      (u < 0)
	  X.insert(tmin + u*TTS);
	else if (u > 0 && (tmax - tmin) < 20*TTS)
	  X.insert(tmax + u*TTS);
      }

      //const double W = transformer.getWeight(array);

      JFunction1D_t buffer;
	
      for (set<double>::const_iterator x = X.begin(); x != X.end(); ++x) {
	
	double y = 0.0;
	
	for (JGaussHermite::const_iterator j = engine.begin(); j != engine.end(); ++j) {
	  
	  const double u = j->getX();
	  const double v = j->getY() / sqrt(PI);
	  const double w = f1(*x + u*TTS);
	  
	  y += v * w;
	}
	
	buffer[*x] = y;
      }
      
      buffer.transform(put);
      buffer.compile();
	    
      i->second->second->second = buffer;
    }
  }
    
  NOTICE("OK" << endl);


  NOTICE("storing output to file " << outputFile << "... " << flush);

  pdf.store(outputFile.c_str());
    
  NOTICE("OK" << endl);
}