#ifndef _utl_Histogram_h_
#define _utl_Histogram_h_

#include <vector>
#include <sstream>
#include <algorithm>

#include <boost/type_traits.hpp>

#include <utl/ErrorLogger.h>
#include <utl/AugerException.h>


namespace utl {

  // this is missing from the boost::type_traits
  template<typename T>
  struct identity {
    typedef T type;
  };


  template<typename T>
  class add_const_reference {
    typedef typename boost::add_const<T>::type const_type;
  public:
    typedef typename boost::add_reference<const_type>::type type;
  };


  /// equal binning policy for Histogram
  class EqualBinning {
  public:
    typedef double BinEdgeType;
    typedef const double& ConstBinLowerEdgesRefType;

    EqualBinning(const size_t nBins, const double start, const double stop)
      : fNBins(nBins), fStart(start), fStop(stop), fBinSize((stop-start)/nBins) { }

    size_t GetNBins() const { return fNBins; }

    double GetStart() const { return fStart; }

    double GetStop() const { return fStop; }

    bool IsUnderflow(const double x) const { return x < fStart; }

    bool IsOverflow(const double x) const { return x >= fStop; }

    double GetBinLowerEdge(const size_t bin) const
    { return fStart + bin*fBinSize; }

    double GetBinUpperEdge(const size_t bin) const
    { return fStart + (bin+1)*fBinSize; }

    double GetBinSize(const size_t) const
    { return fBinSize; }

    size_t GetBinIndexUnchecked(const double x) const
    { return (size_t)((x-fStart) / fBinSize); }

  private:
    size_t fNBins;
    double fStart;
    double fStop;
    double fBinSize;
  };


  /// policy for the variable bin boundaries in Histogram
  template<typename T = double,
           template<typename> class ContainerPolicy = boost::add_const>
  class VariableBinning {
  public:
    typedef T BinEdgeType;
    typedef typename ContainerPolicy<std::vector<BinEdgeType> >::type BinLowerEdgesType;
    typedef typename boost::add_reference<BinLowerEdgesType>::type BinLowerEdgesRefType;
    typedef typename add_const_reference<BinLowerEdgesType>::type ConstBinLowerEdgesRefType;

    VariableBinning() : fBinLowerEdges() { }

    VariableBinning(ConstBinLowerEdgesRefType edges)
      : fBinLowerEdges(edges) { }

    size_t GetNBins() const { return fBinLowerEdges.size() - 1; }

    BinEdgeType GetStart() const { return fBinLowerEdges.front(); }

    BinEdgeType GetStop() const { return fBinLowerEdges.back(); }

    bool IsUnderflow(const double x) const { return x < GetStart(); }

    bool IsOverflow(const double x) const { return x >= GetStop(); }

    BinEdgeType GetBinLowerEdge(const size_t bin) const
    { return fBinLowerEdges[bin]; }

    BinEdgeType GetBinUpperEdge(const size_t bin) const
    { return fBinLowerEdges[bin+1]; }

    BinEdgeType GetBinSize(const size_t bin) const
    { return GetBinUpperEdge(bin) - GetBinLowerEdge(bin); }

    size_t
    GetBinIndexUnchecked(const double x)
      const
    {
      return size_t(std::upper_bound(fBinLowerEdges.begin(),
                                     fBinLowerEdges.end(), x)
                    - fBinLowerEdges.begin()) - 1;
    }

    BinLowerEdgesRefType GetBinLowerEdges() { return fBinLowerEdges; }

    ConstBinLowerEdgesRefType GetBinLowerEdges() const { return fBinLowerEdges; }

  private:
    BinLowerEdgesType fBinLowerEdges;
  };


  /**
    \class Histogram Histogram.h "utl/Histogram.h"

    \brief Histogram

    \author Darko Veberic
    \author Lukas Nellen
    \date 20 Dec 2006
    \version $Id: Histogram.h 14717 2009-09-17 20:24:36Z lukas $

    histogram bin layout
    all user-side indices refer to the base-0 (C-like) numbering of bins,
    underflow is in bin -1, overflow in bin n

    \code
                              start                             stop
                                |                                 |
    edge tendency:              |>  |>  |>  |>  |>    |>    |>    |>
        bin index:              | 0 | 1 | 2 | 3 |     | n-2 | n-1 |     
                   <<_underflow_|___|___|___|___| ... |_____|_____|_overflow_>>
                                |   |   |   |   |     |     |     |
       edge index:              0   1   2   3   4    n-2   n-1    n
    \endcode

  */


  template<typename T = double,
           class BinningPolicy = EqualBinning,
           template<typename> class ContainerPolicy = identity>
  class Histogram : public BinningPolicy {
  public:
    typedef typename BinningPolicy::BinEdgeType BinEdgeType;
    typedef typename ContainerPolicy<std::vector<T> >::type BinContainerType;
    typedef T BinType;
    typedef typename boost::add_reference<BinContainerType>::type BinContainerRefType;
    typedef typename add_const_reference<BinContainerType>::type ConstBinContainerRefType;
    typedef typename BinningPolicy::ConstBinLowerEdgesRefType ConstBinLowerEdgesRefType;

    enum OutOfBounds {
      eUnderflow = -2,
      eOverflow = -1
    };

    // equal binning
    Histogram(const size_t nBins, const BinEdgeType start, const BinEdgeType stop)
      : BinningPolicy(nBins, start, stop),
        fUnderflow(0), fOverflow(0), fBinContainer()
    { fBinContainer.resize(nBins); }

    Histogram(const BinEdgeType start, const BinEdgeType stop,
               ConstBinContainerRefType dataBins)
      : BinningPolicy(dataBins.size(), start, stop),
        fUnderflow(0), fOverflow(0), fBinContainer(dataBins)
    { }

    // variable binning
    Histogram(ConstBinLowerEdgesRefType binEdges)
      : BinningPolicy(binEdges),
        fUnderflow(0), fOverflow(0), fBinContainer()
    { fBinContainer.resize(binEdges.size()-1); }

    Histogram(ConstBinLowerEdgesRefType binEdges,
              ConstBinContainerRefType dataBins)
      : BinningPolicy(binEdges),
        fUnderflow(0), fOverflow(0), fBinContainer(dataBins)
    {
      if (binEdges.size()-1 != dataBins.size()) {
        std::ostringstream err;
        err << "Initializing Histogram with number of bin edges "
            << binEdges.size() << " but number of bins "
            << dataBins.size() << '!';
        ERROR(err);
        throw OutOfBoundException(err.str());
      }
    }

    size_t
    GetBinIndex(const double x)
      const
    {
      if (BinningPolicy::IsUnderflow(x))
        return eUnderflow;
      else if (BinningPolicy::IsOverflow(x))
        return eOverflow;
      else
        return BinningPolicy::GetBinIndexUnchecked(x);
    }

    void Fill(const double x) { ++GetDataBin(x); }

    void Fill(const double x, const BinType weight)
    { GetDataBin(x) += weight; }

    void
    Clear()
    {
      fUnderflow = fOverflow = BinType();
      std::fill(fBinContainer.begin(), fBinContainer.end(), BinType());
    }

    const BinType& GetBin(const size_t i) const { return fBinContainer[i]; }

    double GetBinAverage(const size_t i) const
    { return double(GetBin(i)) / BinningPolicy::GetBinSize(i); }

    double GetBinCenter(const size_t bin) const
    { return BinningPolicy::GetBinLowerEdge(bin) +
             0.5*BinningPolicy::GetBinSize(bin); }

    BinType& GetUnderflow() { return fUnderflow; }
    const BinType& GetUnderflow() const { return fUnderflow; }

    BinType& GetOverflow() { return fOverflow; }
    const BinType& GetOverflow() const { return fOverflow; }

    /*BinValueType
    GetTotalSum()
      const
    {
      BinValueType sum = 0;
      for ( = 0; i < n; ++i)
        sum += fData[i] * BinningPolicy::GetBinSize(i);
      sum += fUnderflow + fOverflow;
      return sum;
    }

    double
    GetBinSum(const size_t low, const size_t high)
      const
    {
      if (high > BinningPolicy::GetNBins())
        throw OutOfBoundException("Out of bin bounds!");
      double sum = 0;
      for (size_t i = low; i < high; ++i)
        sum += fData[i];
      return sum;
    }

    double
    GetBinSum(const bool outOfBounds = true)
      const
    {
      const size_t n = BinningPolicy::GetNBins();
      double sum = GetBinSum(0, n);
      if (outOfBounds)
        sum += fUnderflow + fOverflow;
      return sum;
    }*/

    void SetBin(const size_t i, const BinType value)
    { fBinContainer[i] = value; }

    BinType
    GetMaximum()
      const
    {
      BinType max = fBinContainer[0];
      const size_t n = BinningPolicy::GetNBins();
      for (size_t i = 1; i < n; ++i)
        if (max < fBinContainer[i])
          max = fBinContainer[i];
      return max;
    }

    ConstBinContainerRefType GetBins() const { return fBinContainer; }

  private:
    BinType&
    GetDataBin(const double x)
    {
      if (BinningPolicy::IsUnderflow(x))
        return fUnderflow;
      else if (BinningPolicy::IsOverflow(x))
        return fOverflow;
      else
        return fBinContainer[BinningPolicy::GetBinIndexUnchecked(x)];
    }

    BinType fUnderflow;
    BinType fOverflow;
    BinContainerType fBinContainer;
  };


  //===========================================================================


  template<typename BinEdgesT, typename BinT = double>
  class VariableBinHistogram
    : public Histogram<BinT,
                       VariableBinning<BinEdgesT, boost::add_const>,
                       identity> {
  public:
    VariableBinHistogram(const std::vector<BinEdgesT>& edges)
      : Histogram<BinT,
                  VariableBinning<BinEdgesT, boost::add_const>,
                  identity>(edges)
    { }

    VariableBinHistogram(const std::vector<BinEdgesT>& edges,
                         const std::vector<BinT>& data)
      : Histogram<BinT,
                  VariableBinning<BinEdgesT, boost::add_const>,
                  identity>(edges, data)
    { }
  };


  template<typename BinEdgesT, typename BinT = double>
  class SharedBinHistogram
    : public Histogram<BinT,
                       VariableBinning<BinEdgesT, add_const_reference>,
                       identity> {
  public:
    SharedBinHistogram(const std::vector<BinEdgesT>& edges)
      : Histogram<BinT,
                  VariableBinning<BinEdgesT, add_const_reference>,
                  identity>(edges)
    { }

    SharedBinHistogram(const std::vector<BinEdgesT>& edges,
                       const std::vector<BinT>& data)
      : Histogram<BinT,
                  VariableBinning<BinEdgesT, add_const_reference>,
                  identity>(edges, data)
    { }
  };


  template<typename BinT>
  class HistogramWrap
    : public Histogram<BinT,
                       EqualBinning,
                       add_const_reference> {
  public:
    HistogramWrap(const EqualBinning::BinEdgeType start,
                  const EqualBinning::BinEdgeType stop,
                  const std::vector<BinT>& data)
      : Histogram<BinT,
                  EqualBinning,
                  add_const_reference>(start, stop, data) { }
  };


  template<typename BinEdgesT, typename BinT>
  class VariableBinHistogramWrap
    : public Histogram<BinT,
                       VariableBinning<BinEdgesT, add_const_reference>,
                       add_const_reference> {
  public:
    VariableBinHistogramWrap(const std::vector<BinEdgesT>& edges,
                             const std::vector<BinT>& data)
      : Histogram<BinT,
                  VariableBinning<BinEdgesT, add_const_reference>,
                  add_const_reference>(edges, data)
    { }
  };
}


#endif

// Configure (x)emacs for this file ...
// Local Variables:
// mode: c++
// End: