#ifndef __JTOOLKIT__
#define __JTOOLKIT__

#include "JLang/JException.hh"
#include "JTools/JCollection.hh"
#include "JTools/JElement.hh"
#include "JTools/JVector.hh"
#include "JTools/JPolint.hh"
#include "JTools/JSpline.hh"
#include "JTools/JHistogram.hh"
#include "JTools/JMultiHistogram.hh"
#include "JTools/JQuadrature.hh"


namespace JTOOLS {

  namespace {
    using JLANG::JException;
    using JLANG::JEmptyCollection;
    using JLANG::JDivisionByZero;
  }


  /**
   * Garbage collection.
   * This class implements the JMappableCollection<> interface but does nothing.
   */
  template<class JKey_t, class JValue_t>
  class JGarbageCollection :
    public JMappableCollection<JKey_t, JValue_t>
  {
  public:

    typedef JKey_t     key_type;
    typedef JValue_t   mapped_type;

    /**
     * Value assignment.
     *
     * \param  key          key
     * \return              value
     */
    virtual mapped_type& operator[](const key_type& key)
    {
      throw JException("Run time error.");
    }


    /**
     * Associates the specified value with the specified key.
     *
     * \param  key          key
     * \param  value        value
     * \return              value
     */
    virtual void put(const key_type& key, const mapped_type& value) 
    {}


    /**
     * Clear memory.
     */
    void clear() 
    {}
  };


  /**
   * Conversion of data points to integral values.
   * The integration is based on the trapezoidal rule applied to the input data points.
   *
   * \param  input      JAbstractCollection<>   
   * \param  output     JMappableCollection<>
   * \return            integral
   */
  template<class JElement_t, template<class> class JContainer_t, class JKey_t, class JValue_t>
  inline double convertToIntegral(const JAbstractCollection<JElement_t, JContainer_t>& input, 
				  JMappableCollection<JKey_t, JValue_t>&               output)
  {
    if (input.empty()) 
      throw JEmptyCollection("Method convertToIntegral() no data.");
	      	  
    double V = 0.0;
	  
    output.put(input.begin()->getX(), V);
	    
    for (typename JAbstractCollection<JElement_t, JContainer_t>::const_iterator j = input.begin(), i = j++; j != input.end(); ++i, ++j) {
      
      const double dx = j->getX() - i->getX();
      const double y  = i->getY() + j->getY();
      const double v  = dx * 0.50 * y;
      
      V += v;
      
      output.put(j->getX(), V);
    }

    return V;
  }


  /**
   * Conversion of data points to integral values.
   * The integration includes the 2nd derivatives of the data points of the input spline interpolating function.
   *
   * \param  input      JSplineFunction1D<>
   * \param  output     JMappableCollection<>
   * \return            integral
   */
  template<class JElement_t, template<class> class JCollection_t, class JResult_t, class JKey_t, class JValue_t>
  inline double convertToIntegral(const JSplineFunction1D<JElement_t, JCollection_t, JResult_t>& input, 
				  JMappableCollection<JKey_t, JValue_t>&                         output)
  {
    if (input.empty()) 
      throw JEmptyCollection("Method convertToIntegral() no data.");
	    
    double V = 0.0;
    
    output.put(input.begin()->getX(), V);
    
    for (typename JSplineFunction1D<JElement_t, JCollection_t, JResult_t>::const_iterator j = input.begin(), i = j++; j != input.end(); ++i, ++j) {
      
      const double dx = j->getX() - i->getX();
      const double y  = i->getY() + j->getY();
      const double z  = i->get2ndDerivative() + j->get2ndDerivative();	      
      const double v  = dx * 0.50 * y;
      const double w  = dx * 0.25 * z*dx*dx/6;
      
      V += v - w;
      
      output.put(j->getX(), V);
    }

    return V;
  }


  /**
   * Conversion of data points to integral values.
   * The integration includes the 2nd derivatives of the data points of the input spline interpolating function.
   *
   * \param  input      JSplineFunction1D<>
   * \param  output     JMappableCollection<>
   * \return            integral
   */
  template<class JElement_t, template<class> class JCollection_t, class JResult_t, class JKey_t, class JValue_t>
  inline double convertToIntegral(const JSplineFunction1D<JElement_t, JCollection_t, JFunctionValuePDF<JResult_t> >& input,
				  JMappableCollection<JKey_t, JValue_t>&                                             output)
  {
    if (input.empty()) 
      throw JEmptyCollection("Method convertToIntegral() no data.");
	    
    for (typename JSplineFunction1D<JElement_t, JCollection_t, JResult_t>::const_iterator i = input.begin(); i != input.end(); ++i)
      output.put(i->getX(), i->getIntegral());

    return input.rbegin()->getIntegral();
  }


  /**
   * Conversion of data points to integral values.
   * The integration uses the Gauss-Legendre quadratures with the number of points set
   * to the degree of the input polynomial interpolating function.
   *
   * \param  input      JPolintFunction1D<>
   * \param  output     JMappableCollection<>
   * \return            integral
   */
  template<unsigned int N, class JElement_t, template<class> class JCollection_t, class JKey_t, class JValue_t>
  inline double convertToIntegral(const JPolintFunction1D<N, JElement_t, JCollection_t>& input, 
				  JMappableCollection<JKey_t, JValue_t>&                 output)
  {
    if (input.empty()) 
      throw JEmptyCollection("Method convertToIntegral() no data.");

    const JGaussLegendre engine(N);
    
    double V = 0.0;
    
    output.put(input.begin()->getX(), V);
    
    for (typename JPolintFunction1D<N, JElement_t, JCollection_t>::const_iterator j = input.begin(), i = j++; j != input.end(); ++i, ++j) {
      
      const double xmin = i->getX();
      const double xmax = j->getX();
      
      for (JGaussLegendre::const_iterator m = engine.begin(); m != engine.end(); ++m) {
	
	const double x = 0.5 * (xmax + xmin  +  m->getX() * (xmax - xmin));
	const double y = 0.5 * (xmax - xmin) *  m->getY() * input(x);
	
	V += y;
      }
      
      output.put(j->getX(), V);
    }

    return V;
  }


  /**
   * Conversion of data points to integral values.
   * The integration is based on the sum of y-values of the input data points.
   *
   * \param  input      JPolintFunction1D<0, >
   * \param  output     JMappableCollection<>
   * \return            integral
   */
  template<class JElement_t, template<class> class JCollection_t, class JKey_t, class JValue_t>
  inline double convertToIntegral(const JPolintFunction1D<0, JElement_t, JCollection_t>& input, 
				  JMappableCollection<JKey_t, JValue_t>&                 output)
  {
    if (input.empty()) 
      throw JEmptyCollection("Method convertToIntegral() no data.");

    double V = 0.0;
    
    output.put(input.begin()->getX(), V);
    
    for (typename JPolintFunction1D<0, JElement_t, JCollection_t>::const_iterator j = input.begin(), i = j++; j != input.end(); ++i, ++j) {
      
      const double dx = j->getX() - i->getX();
      const double y  = i->getY();
      const double v  = dx * y;
      
      V += v;
	      
      output.put(j->getX(), V);
    }

    return V;
  }


  /**
   * Conversion of data points to integral values.
   * The integration is based on the trapezoidal rule applied to the input data points.
   *
   * \param  input      JPolintFunction1D<1, >
   * \param  output     JMappableCollection<>
   */
  template<class JElement_t, template<class> class JCollection_t, class JKey_t, class JValue_t>
  inline double convertToIntegral(const JPolintFunction1D<1, JElement_t, JCollection_t>& input, 
				  JMappableCollection<JKey_t, JValue_t>&                 output)
  {
    if (input.empty()) 
      throw JEmptyCollection("Method convertToIntegral() no data.");

    double V = 0.0;
    
    output.put(input.begin()->getX(), V);
    
    for (typename JPolintFunction1D<1, JElement_t, JCollection_t>::const_iterator j = input.begin(), i = j++; j != input.end(); ++i, ++j) {
      
      const double dx = j->getX() - i->getX();
      const double y  = i->getY() + j->getY();
      const double v  = dx * 0.50 * y;
      
      V += v;
      
      output.put(j->getX(), V);
    }

    return V;
  }


  /**
   * Conversion of data points to integral values.
   * The integration is based on the sum of y-values of the input data points.
   *
   * \param  input      JHistogram1D<>   
   * \param  output     JMappableCollection<>
   * \return            integral
   */
  template<class JElement_t, template<class> class JContainer_t, class JKey_t, class JValue_t>
  inline double convertToIntegral(const JHistogram1D<JElement_t, JContainer_t>& input, 
				  JMappableCollection<JKey_t, JValue_t>&        output)
  {
    if (input.empty()) 
      throw JEmptyCollection("Method convertToIntegral() no data.");
	      	  
    double V = 0.0;
    
    output.put(input.begin()->getX(), V);
    
    for (typename JHistogram1D<JElement_t, JContainer_t>::const_iterator j = input.begin(), i = j++; j != input.end(); ++i, ++j) {
      
      V += i->getY();
      
      output.put(j->getX(), V);
    }

    return V;
  }


  /**
   * Get integral of input data points.
   *
   * \param  input      input data
   * \return            integral
   */
  template<class JCollection_t>
  inline double getIntegral(JCollection_t& input)
  { 
    typedef typename JCollection_t::key_type     key_type;
    typedef typename JCollection_t::mapped_type  mapped_type;

    JGarbageCollection<key_type, mapped_type> garbage;

    return convertToIntegral(input, garbage);
  }
  

  /**
   * Conversion of data points to probability density distribition.
   * The y-values of the data points are divided by the integral.
   *
   * \param  input      JAbstractCollection<>   
   * \param  output     JMappableCollection<>
   */
  template<class JElement_t, template<class> class JContainer_t, class JKey_t, class JValue_t>
  inline double convertToPDF(const JAbstractCollection<JElement_t, JContainer_t>& input, 
			     JMappableCollection<JKey_t, JValue_t>&               output)
  {
    if (input.empty()) 
      throw JEmptyCollection("Method convertToPDF() no data.");
    
    for (typename JAbstractCollection<JElement_t, JContainer_t>::const_iterator i = input.begin(); i != input.end(); ++i) {
      
      const JKey_t   x = i->getX();
      const JValue_t y = i->getY();
      
      output.put(x,y);
    }
  }


  /**
   * Conversion of data points to probability density distribition.
   * This histogram abscissa and contents are set at the bin center and divided by the integral and the bin width, respectively.
   *
   * \param  input      JHistogram1D<>
   * \param  output     JMappableCollection<>
   */
  template<class JElement_t, template<class> class JContainer_t, class JKey_t, class JValue_t>
  inline void convertToPDF(const JHistogram1D<JElement_t, JContainer_t>& input, 
			   JMappableCollection<JKey_t, JValue_t>&        output)
  {
    if (input.empty()) 
      throw JEmptyCollection("Method convertToPDF() no data.");

    for (typename JHistogram1D<JElement_t, JContainer_t>::const_iterator j = input.begin(), i = j++; j != input.end(); ++i, ++j) {
	    
      const JKey_t   x = 0.5 * (i->getX() + j->getX());
      const JValue_t y = i->getY() / (j->getX() - i->getX());
      
      output.put(x,y);
    }
  }


  /**
   * Conversion of data points to cumulative probability distribition.
   *
   * \param  input      see template methods convertToIntegral(..)
   * \param  output     JMappableCollection<>
   * \param  eps        minimal step size
   * \return            integral value
   */
  template<class JCollection_t, class JKey_t, class JValue_t>
  inline double convertToCDF(const JCollection_t&                   input,
			     JMappableCollection<JKey_t, JValue_t>& output, 
			     const double eps = 0.0)
  {
    if (input.empty()) 
      throw JEmptyCollection("Method convertToPDF() no data.");

    JVector< JElement2D<JKey_t, JValue_t> > buffer;

    const double V = convertToIntegral(input, buffer);
	  
    if (V == 0.0) 
      throw JDivisionByZero("Method convertToCDF() integral equals zero.");
	    
    output.clear();

    typename JVector< JElement2D<JKey_t, JValue_t> >::const_iterator i = buffer.begin();
	    
    for ( ; i != buffer.end() && i->getY() <= 0.5*eps * V; ++i) {}
    
    if (i != buffer.end()) {
      
      // force first point at (0, ymin)
      
      JKey_t   x = 0.0;
      JValue_t y = i->getX();
      
      output.put(x,y);
      
      JKey_t   xmax = x;
      JValue_t ymax = y;
      
      for (++i; i != buffer.end(); ++i) {
	
	x = i->getY() / V;
	y = i->getX();
	
	if (x > xmax) {
	  
	  ymax = y;
	  
	  if (x > xmax + eps /*&& x < 1.0 - 0.5*eps*/) {
	    
	    output.put(x,y);
	    
	    xmax = x;
	  }
	}
      }
      
      // force last point at (1, ymax)
      
      x = 1.0;
      y = ymax;

      output.put(x,y);

    } else

      throw JDivisionByZero("Method convertToCDF() integral equals zero.");

    return V;
  }
}

#endif