Histogram.h

Go to the documentation of this file.

#ifndef _utl_Histogram_h_
#define _utl_Histogram_h_

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

#include <boost/type_traits.hpp>

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


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;
  };


  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:
    const size_t fNBins;
    const double fStart;
    const double fStop;
    const double fBinSize;
  };


  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:
    const BinLowerEdgesType fBinLowerEdges;
  };


  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; }

    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;
    }

    BinType
    GetSumEntries()
      const
    {
      BinType sum = 0;
      for (const auto& entry : fBinContainer)
        sum += entry;
      return sum;
    }

    ConstBinContainerRefType GetBins() const { return fBinContainer; }

    double
    GetEntropy()
      const
    {
      double n = 0;
      double clnc = 0;
      for (size_t bin = 0; bin < BinningPolicy::GetNBins(); ++bin) {
        if (const double c = GetBinAverage(bin)) {
          n += c;
          if (c > 0)
            clnc -= c * std::log(c);
        }
      }
      return (n > 0) ? clnc / n + std::log(n) : 0;
    }

    std::vector<double>
    GetDensityVector()
      const
    {
      const auto n = BinningPolicy::GetNBins();
      std::vector<double> density;
      density.reserve(n + 1);
      density.push_back(GetBinAverage(0));
      for (size_t i = 1; i < n - 2; ++i) {
        const double value1 = GetBinAverage(i - 1);
        const double value2 = GetBinAverage(i);
        const double value3 = GetBinAverage(i + 1);
        density.push_back((1/3.)*(value1 + value2 + value3));
      }
      density.push_back(0);
      density.push_back(0);
      return density;
    }

    std::vector<double>
    GetDensityErrorVector()
      const
    {
      const auto n = BinningPolicy::GetNBins();
      std::vector<double> densityError;
      densityError.reserve(n + 1);
      densityError.push_back(sqrt(GetBinAverage(0) / BinningPolicy::GetBinSize(0)));
      for (size_t i = 1; i < n - 2; ++i) {
        const double value1 = GetBinAverage(i - 1) / BinningPolicy::GetBinSize(i - 1);
        const double value2 = GetBinAverage(i) / BinningPolicy::GetBinSize(i);
        const double value3 = GetBinAverage(i + 1) / BinningPolicy::GetBinSize(i + 1);
        densityError.push_back((1/3.)*sqrt(value1 + value2 + value3));
      }
      densityError.push_back(0);
      densityError.push_back(0);
      return densityError;
    }

  private:
    BinType&
    GetDataBin(const double x)
    {
      if (BinningPolicy::IsUnderflow(x))
        return fUnderflow;
      if (BinningPolicy::IsOverflow(x))
        return fOverflow;
      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