#ifndef __JOSCPROB__JOSCPROBINTERPOLATOR__
#define __JOSCPROB__JOSCPROBINTERPOLATOR__

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

#include "JIO/JSerialisable.hh"
#include "JIO/JFileStreamIO.hh"

#include "JLang/JClonable.hh"
#include "JLang/JObjectIO.hh"
#include "JLang/JException.hh"

#include "JTools/JPolint.hh"
#include "JTools/JMapList.hh"
#include "JTools/JCollection.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JMultiFunction.hh"
#include "JTools/JFunctionalMap_t.hh"

#include "JOscProb/JOscChannel.hh"
#include "JOscProb/JOscParameters.hh"
#include "JOscProb/JOscProbToolkit.hh"
#include "JOscProb/JBaselineCalculator.hh"
#include "JOscProb/JOscProbInterpolatorInterface.hh"


/**
 * \author bjung, mdejong
 */

namespace JOSCPROB {}
namespace JPP { using namespace JOSCPROB; }

namespace JOSCPROB {

  using JEEP::JMessage;

  using JLANG::JClonable;

  using JTOOLS::JMultiFunction;
  

  /** 
   * Template definition of a multi-dimensional oscillation probability interpolation table.
   */
  template<template<class, class> class JCollection_t = JTOOLS::JCollection,
	   class JFunction1D_t                        = JTOOLS::JPolintFunction1D <1,
										   JTOOLS::JElement2D<double, JTOOLS::JArray<NUMBER_OF_OSCCHANNELS, double> >,
										   JCollection_t,
										   JTOOLS::JArray<NUMBER_OF_OSCCHANNELS, double> >,
	   class JFunctionalMaplist_t                 = JTOOLS::JMAPLIST          <JTOOLS::JPolint1FunctionalMap,
										   JTOOLS::JPolint1FunctionalMap,
										   JTOOLS::JPolint1FunctionalMap,
										   JTOOLS::JPolint1FunctionalMap,
										   JTOOLS::JPolint1FunctionalMap,
										   JTOOLS::JPolint1FunctionalMap,
										   JTOOLS::JPolint2FunctionalMap>::maplist >
  class JOscProbInterpolator :
    public JMultiFunction                                                               <JFunction1D_t, JFunctionalMaplist_t>,
    public JClonable<JOscProbInterpolatorInterface, JOscProbInterpolator <JCollection_t, JFunction1D_t, JFunctionalMaplist_t> >,
    public JMessage                                <JOscProbInterpolator <JCollection_t, JFunction1D_t, JFunctionalMaplist_t> >
  {
  public:

    typedef JOscProbInterpolator<JCollection_t, JFunction1D_t, JFunctionalMaplist_t>                interpolator_type;
    typedef JMultiFunction<JFunction1D_t, JFunctionalMaplist_t>                                    multifunction_type;

    enum { NUMBER_OF_DIMENSIONS = multifunction_type::NUMBER_OF_DIMENSIONS };

    typedef typename multifunction_type::abscissa_type                                                  abscissa_type;    
    typedef typename multifunction_type::argument_type                                                  argument_type;
    typedef typename multifunction_type::result_type                                                      result_type;
    typedef typename multifunction_type::value_type                                                        value_type;

    typedef typename multifunction_type::multimap_type                                                  multimap_type;

    typedef typename multifunction_type::super_const_iterator                                    super_const_iterator;
    typedef typename multifunction_type::super_iterator                                                super_iterator;
    typedef typename multifunction_type::function_type                                                  function_type;

    using JMessage<interpolator_type>::debug;

    
    /**
     * Default constructor.
     */
    JOscProbInterpolator() :
      multifunction_type(),
      parameters(),
      getBaseline()
    {
      this->set(JOscParameters(false)); // Initialize buffer with NuFIT NO best fit parameters
    }


    /**
     * Constructor.
     *
     * \param  fileName             oscillation probability table filename
     */
    JOscProbInterpolator(const char* fileName) :
      multifunction_type(),
      parameters(),
      getBaseline()
    {
      this->load(fileName);
      this->set(JOscParameters(false)); // Initialize buffer with NuFIT NO best fit parameters      
    }


    /**
     * Constructor.
     *
     * \param  fileName             oscillation probability table filename
     * \param  parameters           oscillation parameters
     */
    JOscProbInterpolator(const char*           fileName,
			 const JOscParameters& parameters) :
      multifunction_type(),
      parameters(),
      getBaseline()
    {
      this->load(fileName);
      this->set(parameters);
    }


    /**
     * Load oscillation probability table from file.
     *
     * \param  fileName             oscillation probability table filename
     */
    void load(const char* fileName) override
    {
      using namespace std;
      using namespace JPP;
      
      try {

	NOTICE("loading oscillation probability table from file " << fileName << "... " << flush);
	
	JLANG::load<JIO::JFileStreamReader>(fileName, *this);

	NOTICE("OK" << endl);
      }
      catch(const JException& error) {
	THROW(JFileReadException, "JOscProbInterpolator::load(): Error reading file " << fileName);
      }
    }    


    /**
     * Get fixed oscillation parameters associated with this interpolation table.
     *
     * \return                    oscillation parameters
     */
    const JOscParameters& getTableParameters() const override
    {
      return parameters;
    }


    /**
     * Get baseline calculator associated with this interpolation table.
     *
     * \return                    baseline calculator
     */
    const JBaselineCalculator& getBaselineCalculator() const override
    {
      return getBaseline;
    }


    /**
     * Set oscillation parameters for interpolation.
     *
     * \param  parameters           oscillation parameters
     */
    void set(const JOscParameters& parameters) const override
    {
      using namespace JPP;

      JOscParameters pars(parameters);
      pars.join(this->parameters);
      
      const JProperties properties = pars.getProperties();

      for (JProperties::const_iterator i = properties.cbegin(); i != properties.cend(); ++i) {
	
	const JOscParameters::JParameter_t& parameter = i->second.getValue<JOscParameters::JParameter_t>();
	
	if (parameter.isDefined()) {
	    
	  const int index = std::distance(properties.cbegin(), i);

	  this->buffer[index] = parameter.getValue();
	  
	} else {

	  THROW(JNoValue,
		"JOscProbInterpolator<...>::set(JOscParameters): " <<
		"No value for parameter "                          << i->first);
	}
      }
    }


    /**
     * Get oscillation probability for a given oscillation channel.
     *
     * \param  channel            oscillation channel
     * \param  E                  neutrino energy [GeV]
     * \param  costh              cosine zenith angle
     * \return                    oscillation probability
     */
    double operator()(const JOscChannel&    channel,
		      const double          E,
		      const double          costh) const override
    {
      using namespace std;
      using namespace JPP;

      const JOscChannel* p = find(getOscChannel, getOscChannel + NUMBER_OF_OSCCHANNELS, channel);
	
      if (p != end(getOscChannel)) {
	
	const double L = getBaseline(costh);

	this->buffer[NUMBER_OF_DIMENSIONS-2] = L/E;
	this->buffer[NUMBER_OF_DIMENSIONS-1] = costh;

	const argument_type* arguments = this->buffer.data();
	
	const size_t       index         = distance(getOscChannel, p);
	const result_type& probabilities = this->evaluate(arguments);

	const double& P = probabilities[index];

	return (P > 1.0 ? 1.0 : (P < 0.0 ? 0.0 : P));
	
      } else {
	
	THROW(JValueOutOfRange, "JOscProbInterpolator<...>::operator(): Invalid oscillation channel " << channel << endl);
      }
    }


    /**
     * Get oscillation probability for a given set of oscillation parameters\n
     * and a given oscillation channel.
     *
     * \param  channel            oscillation channel
     * \param  parameters         oscillation parameters
     * \param  E                  neutrino energy [GeV]
     * \param  costh              cosine zenith angle
     * \return                    oscillation probability
     */
    double operator()(const JOscParameters& parameters,
		      const JOscChannel&    channel,
		      const double          E,
		      const double          costh) const override
    {
      set(parameters);

      return (*this)(channel, E, costh);      
    }
    
    
    /**
     * Read from input.
     *
     * \param  in                   reader
     * \return                      reader
     */
    JReader& read(JReader& in) override 
    {
      parameters.read(in);
      getBaseline.read(in);
      
      in >> static_cast<multifunction_type&>(*this);
      
      this->compile();
      
      return in;
    }
    
    
    /**
     * Write from input.
     *
     * \param  out                  writer
     * \return                      writer
     */
    JWriter& write(JWriter& out) const override 
    {
      parameters.write(out);
      getBaseline.write(out);
      
      out << static_cast<const multifunction_type&>(*this);

      return out;
    }

    
  private:

    JOscParameters      parameters;  //!< Fixed oscillation parameters corresponding to the oscillation probability table
    JBaselineCalculator getBaseline; //!< Baseline functor
  };
}


namespace JEEP {

  /**
   * JMessage template specialization for oscillation probability interpolators.
   */
  template<template<class, class> class JCollection_t, class JFunction1D_t, class JFunctionalMaplist_t>
  struct JMessage<JOSCPROB::JOscProbInterpolator<JCollection_t, JFunction1D_t, JFunctionalMaplist_t> >
  {
    static int debug;
  };


  /**
   * Default verbosity for oscillation probability interpolators.
   */
  template<template<class, class> class JCollection_t, class JFunction1D_t, class JFunctionalMaplist_t>
  int JMessage<JOSCPROB::JOscProbInterpolator<JCollection_t, JFunction1D_t, JFunctionalMaplist_t> >::debug = (int) notice_t;
  
}

#endif