#ifndef HISTOGRAM_H
#define HISTOGRAM_H

#include <cassert>
#include <iostream>
#include <vector>
using namespace std;

#include <typedefs.h>

/*!
\brief Count instances in 1D. Equally spaced intervals
  except the end intervals (-infinty,xlow) and (xhigh,infinity).

xlow is the lowest partition point.
  xhigh is the highest partition point.
  N is the number of inner divisions. 
    ie N-1 partitions between xlow and xhigh.
  The other two partitions are for x>xhigh and x<xlow.
  n+1 partitions in total.
*/
template< typename T >
class histogram
{
  T xlow;
  T xhigh;
  uint N;
  T dw;
public:

  /** Keeps interval counts. */
  vector<T> counter;

  /** Sets counts to zero. */
  void reset();

  /** Constructor defines the intervals. */
  histogram
  (
    T const xlow_,
    T const xhigh_,
    uintc N_
  );

  /** Count the instance. */
  void eval(T const x);

  /** Construct a frequency distribution from data.
      ie c[i]/sum(c[i]) */
  void frequencydist
  (
    vector<T> & vd,
    boolc exclude_low=true,  // Ignore lowest interval
    boolc exclude_high=true  // Ignore highest interval
  ) const;

  /** Interlace the middle of the range and the frequency 
      distribution. Useful for plots. */
  void frequencydistpair
  (
    vector<T> & vd,
    boolc exclude_low=true,  // Ignore lowest interval
    boolc exclude_high=true  // Ignore highest interval
  ) const;

  /** Same as frequencydistpair but prints a return after
      each pair. */
  ostream & printfrequencydistpair
  (
    ostream & os,
    boolc exclude_low=true,  // Ignore lowest interval
    boolc exclude_high=true  // Ignore highest interval
  ) const;

  /** Prints the counter including lowest and highest intervals. */
  ostream & print(ostream & os) const;
};

template< typename T >
ostream & operator << (ostream & os, histogram<T> const & hg)
{
  return hg.print(os);
}


//----------------------------------------------------------
// Implementation


template< typename T >
inline void histogram<T>::eval(T const x)
{
  assert(counter.size()==N+1);
  assert(N>0);
  if ((xlow<x)&&(x<xhigh))
  {
    uint k = (uint)((x-xlow)*dw+1);
    ++counter[k];
    return;
  }

  if ((xlow<x)==false)
  {
    ++counter[0];
    return;
  }

  ++counter[N];
}

template< typename T >
ostream & histogram<T>::printfrequencydistpair
(
  ostream & os,
  boolc exclude_low,  // Ignore lowest interval
  boolc exclude_high  // Ignore highest interval
) const
{
  vector<T> fd;
  frequencydistpair(fd,exclude_low,exclude_high);
  uint imax=fd.size()/2;
  for (uint i=0; i<imax; ++i)
  {
    os << fd[2*i] << " ";
    os << fd[2*i+1] << endl;
  }

  return os;
}

template< typename T >
void histogram<T>::frequencydistpair
(
  vector<T> & vd,
  boolc exclude_low,
  boolc exclude_high
) const
{
  if (counter.empty()||(N<3))
    return;

  vector<T> v;
  frequencydist(v,exclude_low,exclude_high);

  uintc imax = v.size();
  vd.resize(imax*2);
  for (uint i=0; i<imax; ++i)
  {
    vd[i*2] = xlow+(xhigh-xlow)*(.5+i)/(N-1);
    vd[i*2+1] = v[i];
  }

}

template< typename T >
void histogram<T>::frequencydist
(
  vector<T> & vd,
  boolc exclude_low,
  boolc exclude_high
) const
{
  if (counter.empty()||(N<3))
    return;

  uint i=0;
  uint imax=counter.size();
  if (exclude_low)
    ++i;

  if (exclude_high)
    --imax;

  T sum = 0;
  for (uint k=i;k<imax;++k)
    sum += counter[k];

  vd.resize(imax-i);

  if (sum==0)
    return;

  T suminv = T(1.0) / sum;
  for (uint k=i;k<imax;++k)
    vd[k-i] = counter[k]*suminv;
}

template< typename T >
ostream & histogram<T>::print(ostream & os) const
{
  uint imax = counter.size();
  for (uint i=0; i<imax; ++i)
    os << counter[i] << " ";

  return os;
}

template< typename T >
histogram<T>::histogram
(
  T const xlow_,
  T const xhigh_,
  uintc N_
)
  : xlow(xlow_), xhigh(xhigh_), N(N_)
{
  assert(xlow<xhigh);
  assert(N>0);

  dw = (N-1.0)/(xhigh-xlow);
  reset();
}

template< typename T >
void histogram<T>::reset()
{
  counter.resize(N+1);
  assert(counter.size()==N+1);
  uintc imax=counter.size();

  for (uint i=0; i<imax; ++i)
    counter[i]=0;
}

#endif