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

#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JPhysics/JPDFTable.hh"

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


/**
 * \file
 *
 * Auxiliary program to compare PDF tables of the arrival time of the Cherenkov light from a shower.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;

  string      inputFileA;
  string      inputFileB;
  double      precision;
  int         debug;

  try { 

    JParser<> zap("Auxiliary program to compare PDF tables of the arrival time of the Cherenkov light from a shower.");

    zap['a'] = make_field(inputFileA);
    zap['b'] = make_field(inputFileB);
    zap['p'] = make_field(precision)        = 1.0e-4;
    zap['d'] = make_field(debug)            = 0;

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


  using namespace JPP;

  typedef JSplineFunction1D_t                                     JFunction1D_t;
  typedef JMAPLIST<JPolint1FunctionalMap,
		   JPolint1FunctionalMap,
		   JPolint1FunctionalGridMap,
		   JPolint1FunctionalGridMap>::maplist            JMapList_t;
  typedef JPDFTable<JFunction1D_t, JMapList_t>                    JPDF_t;

  typedef JFunction1D_t::argument_type                            argument_type;
  typedef JArray<4, argument_type>                                JArray_t;


  JPDF_t pdfA;
  JPDF_t pdfB;


  const JFunction1D_t::JSupervisor supervisor(new JFunction1D_t::JDefaultResult(0.0));


  struct JEntry_t {

    const char* file_name;
    JPDF_t&     pdf;

    void load() 
    { 
      pdf.load(file_name); 
    }

  } buffer[] = { 
    { inputFileA.c_str(), pdfA },
    { inputFileB.c_str(), pdfB }
  };


  for (JEntry_t* p = buffer; p != buffer + sizeof(buffer)/sizeof(buffer[0]); ++p) {

    try {
      
      NOTICE("loading input from file " << p->file_name << "... " << flush);
      
      p->load();
      
      NOTICE("OK" << endl);
    }
    catch(const JException& error) {
      FATAL(error.what() << endl);
    }

    p->pdf.compile();
    p->pdf.setExceptionHandler(supervisor);
  }    

  bool ok = true;

  JPDF_t::super_const_iterator i = pdfA.super_begin();
  JPDF_t::super_const_iterator j = pdfB.super_begin();

  for ( ; i != pdfA.super_end() && 
	  j != pdfB.super_end(); ++i, ++j) {

    const double Wa = pdfA.transformer->getWeight(JArray_t((*i).getKey()));
    const double Wb = pdfB.transformer->getWeight(JArray_t((*j).getKey()));

    if (fabs(i->first                 - j->first)                                  >  precision ||
	fabs(i->second->first         - j->second->first)                          >  precision ||
	fabs(i->second->second->first - j->second->second->first)                  >  precision ||
	fabs(i->second->second->second->first - j->second->second->second->first)  >  precision ||
	fabs(Wa                       - Wb)                                        >  precision) {

      DEBUG("a[] " 
	    << i->first                         << ' '
	    << i->second->first                 << ' '
	    << i->second->second->first         << ' '
	    << i->second->second->second->first << ' '
	    << Wa                               << endl);

      DEBUG("b[] " 
	    << j->first                         << ' '
	    << j->second->first                 << ' '
	    << j->second->second->first         << ' '
	    << j->second->second->second->first << ' '
	    << Wb                               << endl);

      ok = false;
    }


    const JFunction1D_t& fa = (*i).getValue();
    const JFunction1D_t& fb = (*j).getValue();

    JFunction1D_t::const_iterator p = fa.begin();
    JFunction1D_t::const_iterator q = fb.begin();

    while (p != fa.end() && 
	   q != fb.end()) {
      
      if (fabs(p->getX() - q->getX()) > precision) {

	DEBUG("a.f[] " << p->getX() << endl);
	DEBUG("b.f[] " << q->getX() << endl);

	ok = false;
      }

      if (fabs(p->getY() - q->getY()) > precision) {

	DEBUG("a.f() " << p->getX() << ' ' << p->getY() << endl);
	DEBUG("b.f() " << q->getX() << ' ' << q->getY() << endl);

	ok = false;
      }

      if      (p->getX() < q->getX())
	++p;
      else if (q->getX() < p->getX())
	++q;
      else {
	++p;
	++q;
      }
    }

    for ( ; p != fa.end(); ++p) {

      DEBUG("a.f() " << p->getX() << endl);

      ok = false;
    }

    for ( ; q != fb.end(); ++q) {

      DEBUG("b.f() " << q->getX() << endl);

      ok = false;
    }
  }

  for ( ; i != pdfA.super_end(); ++i) {

    DEBUG("a[] " 
	  << i->first                         << ' '
	  << i->second->first                 << ' '
	  << i->second->second->first         << ' '
	  << i->second->second->second->first << endl); 
    
    ok = false;
  }

  for ( ; j != pdfB.super_end(); ++j) {

    DEBUG("b[] " 
	  << j->first                         << ' '
	  << j->second->first                 << ' '
	  << j->second->second->first         << ' '
	  << j->second->second->second->first << endl);

    ok = false;
  }

  if (!ok) {
    ERROR(inputFileA << " and " << inputFileB << " differ." << endl);
  }
}