view maths/Histogram.h @ 56:a0f987c06bec

* Make it possible to retrieve beat spectral difference function from
  bar detector
* Avoid crashes when bar detector is run with very short input

/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */

// Histogram.h: interface for the THistogram class.
//
//////////////////////////////////////////////////////////////////////


#ifndef HISTOGRAM_H
#define HISTOGRAM_H


#include <valarray>

/*! \brief A histogram class

  This class computes the histogram of a vector.

\par Template parameters
	
	- T type of input data (can be any: float, double, int, UINT, etc...)
	- TOut type of output data: float or double. (default is double)

\par Moments:

	The moments (average, standard deviation, skewness, etc.) are computed using 
the algorithm of the Numerical recipies (see Numerical recipies in C, Chapter 14.1, pg 613).

\par Example:

	This example shows the typical use of the class:
\code
	// a vector containing the data
	vector<float> data;
	// Creating histogram using float data and with 101 containers,
	THistogram<float> histo(101);
	// computing the histogram
	histo.compute(data);
\endcode

Once this is done, you can get a vector with the histogram or the normalized histogram (such that it's area is 1):
\code
	// getting normalized histogram
	vector<float> v=histo.getNormalizedHistogram();
\endcode

\par Reference
	
	Equally spaced acsissa function integration (used in #GetArea): Numerical Recipies in C, Chapter 4.1, pg 130.

\author Jonathan de Halleux, [e-mail], 2002
*/

template<class T, class TOut = double>
class THistogram  
{
public:
    //! \name Constructors
    //@{
    /*! Default constructor
      \param counters the number of histogram containers (default value is 10)
    */
    THistogram(unsigned int counters = 10);
    virtual ~THistogram()				{ clear();};
    //@}

    //! \name Histogram computation, update
    //@{
    /*! Computes histogram of vector v
      \param v vector to compute histogram
      \param computeMinMax set to true if min/max of v have to be used to get the histogram limits	
	
      This function computes the histogram of v and stores it internally.
      \sa Update, GeTHistogram
    */
    void compute( const std::vector<T>& v, bool computeMinMax = true);
    //! Update histogram with the vector v
    void update( const std::vector<T>& v);
    //! Update histogram with t
    void update( const T& t);
    //@}

    //! \name Resetting functions
    //@{
    //! Resize the histogram. Warning this function clear the histogram.
    void resize( unsigned int counters );
    //! Clears the histogram
    void clear()						{ m_counters.clear();};
    //@}

    //! \name Setters
    //@{
    /*! This function sets the minimum of the histogram spectrum. 
      The spectrum is not recomputed, use it with care
    */
    void setMinSpectrum( const T& min )					{	m_min = min; computeStep();};
    /*! This function sets the minimum of the histogram spectrum. 
      The spectrum is not recomputed, use it with care
    */
    void setMaxSpectrum( const T& max )					{	m_max = max; computeStep();};
    //@}
    //! \name Getters
    //@{
    //! return minimum of histogram spectrum
    const T& getMinSpectrum() const							{	return m_min;};
    //! return maximum of histogram spectrum
    const T& getMaxSpectrum() const							{	return m_max;};
    //! return step size of histogram containers
    TOut getStep() const									{	return m_step;};
    //! return number of points in histogram
    unsigned int getSum() const;
    /*! \brief returns area under the histogram 

    The Simpson rule is used to integrate the histogram.
    */
    TOut getArea() const;

    /*! \brief Computes the moments of the histogram

    \param data dataset
    \param ave mean
    \f[ \bar x = \frac{1}{N} \sum_{j=1}^N x_j\f]
    \param adev mean absolute deviation
    \f[ adev(X) = \frac{1}{N} \sum_{j=1}^N | x_j - \bar x |\f]
    \param var average deviation:
    \f[ \mbox{Var}(X) = \frac{1}{N-1} \sum_{j=1}^N (x_j - \bar x)^2\f]
    \param sdev standard deviation:
    \f[ \sigma(X) = \sqrt{var(\bar x) }\f]
    \param skew skewness
    \f[ \mbox{Skew}(X) = \frac{1}{N}\sum_{j=1}^N \left[ \frac{x_j - \bar x}{\sigma}\right]^3\f]
    \param kurt kurtosis
    \f[ \mbox{Kurt}(X) = \left\{ \frac{1}{N}\sum_{j=1}^N \left[ \frac{x_j - \bar x}{\sigma}\right]^4 \right\} - 3\f]

    */
    static void getMoments(const std::vector<T>& data, TOut& ave, TOut& adev, TOut& sdev, TOut& var, TOut& skew, TOut& kurt);

    //! return number of containers
    unsigned int getSize() const							{	return m_counters.size();};
    //! returns i-th counter
    unsigned int operator [] (unsigned int i) const {
//	ASSERT( i < m_counters.size() );
        return m_counters[i];
    };
    //! return the computed histogram
    const std::vector<unsigned int>& geTHistogram() const	{	return m_counters;};
    //! return the computed histogram, in TOuts
    std::vector<TOut> geTHistogramD() const;
    /*! return the normalized computed histogram 

    \return the histogram such that the area is equal to 1
    */
    std::vector<TOut> getNormalizedHistogram() const;
    //! returns left containers position
    std::vector<TOut> getLeftContainers() const;
    //! returns center containers position
    std::vector<TOut> getCenterContainers() const;
    //@}
protected:
    //! Computes the step
    void computeStep()								{	m_step = (TOut)(((TOut)(m_max-m_min)) / (m_counters.size()-1));};
    //! Data accumulators
    std::vector<unsigned int> m_counters;
    //! minimum of dataset
    T m_min;
    //! maximum of dataset
    T m_max;
    //! width of container
    TOut m_step;
};

template<class T, class TOut>
THistogram<T,TOut>::THistogram(unsigned int counters)
    : m_counters(counters,0), m_min(0), m_max(0), m_step(0)
{

}

template<class T, class TOut>
void THistogram<T,TOut>::resize( unsigned int counters )
{
    clear();

    m_counters.resize(counters,0);

    computeStep();
}

template<class T, class TOut>
void THistogram<T,TOut>::compute( const std::vector<T>& v, bool computeMinMax)
{
    using namespace std;
    unsigned int i;
    int index;

    if (m_counters.empty())
	return;

    if (computeMinMax)
    {
	m_max = m_min = v[0];
	for (i=1;i<v.size();i++)
	{
	    m_max = std::max( m_max, v[i]);
	    m_min = std::min( m_min, v[i]);
	}
    }

    computeStep();

    for (i = 0;i < v.size() ; i++)
    {
	index=(int) floor( ((TOut)(v[i]-m_min))/m_step ) ;

	if (index >= m_counters.size() || index < 0)
	    return;

	m_counters[index]++;
    }
}

template<class T, class TOut>
void THistogram<T,TOut>::update( const std::vector<T>& v)
{
    if (m_counters.empty())
	return;

    computeStep();

    TOut size = m_counters.size();

    int index;
    for (unsigned int i = 0;i < size ; i++)
    {
	index = (int)floor(((TOut)(v[i]-m_min))/m_step);

	if (index >= m_counters.size() || index < 0)
	    return;

	m_counters[index]++;
    }
}

template<class T, class TOut>
void THistogram<T,TOut>::update( const T& t)
{	
    int index=(int) floor( ((TOut)(t-m_min))/m_step ) ;

    if (index >= m_counters.size() || index < 0)
	return;

    m_counters[index]++;
};

template<class T, class TOut>
std::vector<TOut> THistogram<T,TOut>::geTHistogramD() const
{
    std::vector<TOut> v(m_counters.size());
    for (unsigned int i = 0;i<m_counters.size(); i++)
	v[i]=(TOut)m_counters[i];

    return v;
}

template <class T, class TOut>
std::vector<TOut> THistogram<T,TOut>::getLeftContainers() const
{
    std::vector<TOut> x( m_counters.size());

    for (unsigned int i = 0;i<m_counters.size(); i++)
	x[i]= m_min + i*m_step;

    return x;
}

template <class T, class TOut>
std::vector<TOut> THistogram<T,TOut>::getCenterContainers() const
{
    std::vector<TOut> x( m_counters.size());

    for (unsigned int i = 0;i<m_counters.size(); i++)
	x[i]= m_min + (i+0.5)*m_step;

    return x;
}

template <class T, class TOut>
unsigned int THistogram<T,TOut>::getSum() const
{
    unsigned int sum = 0;
    for (unsigned int i = 0;i<m_counters.size(); i++)
	sum+=m_counters[i];

    return sum;
}

template <class T, class TOut>
TOut THistogram<T,TOut>::getArea() const
{
    const size_t n=m_counters.size();
    TOut area=0;

    if (n>6)
    {
	area=3.0/8.0*(m_counters[0]+m_counters[n-1])
	    +7.0/6.0*(m_counters[1]+m_counters[n-2])
	    +23.0/24.0*(m_counters[2]+m_counters[n-3]);
	for (unsigned int i=3;i<n-3;i++)
	{
	    area+=m_counters[i];
	}
    }
    else if (n>4)
    {
	area=5.0/12.0*(m_counters[0]+m_counters[n-1])
	    +13.0/12.0*(m_counters[1]+m_counters[n-2]);
	for (unsigned int i=2;i<n-2;i++)
	{
	    area+=m_counters[i];
	}
    }
    else if (n>1)
    {
	area=1/2.0*(m_counters[0]+m_counters[n-1]);
	for (unsigned int i=1;i<n-1;i++)
	{
	    area+=m_counters[i];
	}
    }
    else 
	area=0;

    return area*m_step;
}

template <class T, class TOut>
std::vector<TOut> THistogram<T,TOut>::getNormalizedHistogram() const
{
    std::vector<TOut> normCounters( m_counters.size());
    TOut area = (TOut)getArea();

    for (unsigned int i = 0;i<m_counters.size(); i++)
    {
	normCounters[i]= (TOut)m_counters[i]/area;
    }

    return normCounters;
};

template <class T, class TOut>
void THistogram<T,TOut>::getMoments(const std::vector<T>& data, TOut& ave, TOut& adev, TOut& sdev, TOut& var, TOut& skew, TOut& kurt)
{
    int j;
    double ep=0.0,s,p;
    const size_t n = data.size();

    if (n <= 1)
	// nrerror("n must be at least 2 in moment");
	return;

    s=0.0; // First pass to get the mean.
    for (j=0;j<n;j++)
	s += data[j];
	
    ave=s/(n);
    adev=var=skew=kurt=0.0; 
    /* Second pass to get the first (absolute), second,
       third, and fourth moments of the
       deviation from the mean. */

    for (j=0;j<n;j++) 
    {
	adev += fabs(s=data[j]-(ave));
	ep += s;
	var += (p=s*s);
	skew += (p *= s);
	kurt += (p *= s);
    }


    adev /= n;
    var=(var-ep*ep/n)/(n-1); // Corrected two-pass formula.
    sdev=sqrt(var); // Put the pieces together according to the conventional definitions. 
    if (var) 
    {
	skew /= (n*(var)*(sdev));
	kurt=(kurt)/(n*(var)*(var))-3.0;
    } 
    else
	//nrerror("No skew/kurtosis when variance = 0 (in moment)");
	return;
}

#endif