cannam@0: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
cannam@0: 
cannam@0: // Histogram.h: interface for the THistogram class.
cannam@0: //
cannam@0: //////////////////////////////////////////////////////////////////////
cannam@0: 
cannam@0: 
cannam@0: #ifndef HISTOGRAM_H
cannam@0: #define HISTOGRAM_H
cannam@0: 
cannam@0: 
cannam@0: #include <valarray>
cannam@0: 
cannam@0: /*! \brief A histogram class
cannam@0: 
cannam@0:   This class computes the histogram of a vector.
cannam@0: 
cannam@0: \par Template parameters
cannam@0: 	
cannam@0: 	- T type of input data (can be any: float, double, int, UINT, etc...)
cannam@0: 	- TOut type of output data: float or double. (default is double)
cannam@0: 
cannam@0: \par Moments:
cannam@0: 
cannam@0: 	The moments (average, standard deviation, skewness, etc.) are computed using 
cannam@0: the algorithm of the Numerical recipies (see Numerical recipies in C, Chapter 14.1, pg 613).
cannam@0: 
cannam@0: \par Example:
cannam@0: 
cannam@0: 	This example shows the typical use of the class:
cannam@0: \code
cannam@0: 	// a vector containing the data
cannam@0: 	vector<float> data;
cannam@0: 	// Creating histogram using float data and with 101 containers,
cannam@0: 	THistogram<float> histo(101);
cannam@0: 	// computing the histogram
cannam@0: 	histo.compute(data);
cannam@0: \endcode
cannam@0: 
cannam@0: Once this is done, you can get a vector with the histogram or the normalized histogram (such that it's area is 1):
cannam@0: \code
cannam@0: 	// getting normalized histogram
cannam@0: 	vector<float> v=histo.getNormalizedHistogram();
cannam@0: \endcode
cannam@0: 
cannam@0: \par Reference
cannam@0: 	
cannam@0: 	Equally spaced acsissa function integration (used in #GetArea): Numerical Recipies in C, Chapter 4.1, pg 130.
cannam@0: 
cannam@0: \author Jonathan de Halleux, dehalleux@auto.ucl.ac.be, 2002
cannam@0: */
cannam@0: 
cannam@0: template<class T, class TOut = double>
cannam@0: class THistogram  
cannam@0: {
cannam@0: public:
cannam@0:     //! \name Constructors
cannam@0:     //@{
cannam@0:     /*! Default constructor
cannam@0:       \param counters the number of histogram containers (default value is 10)
cannam@0:     */
cannam@0:     THistogram(unsigned int counters = 10);
cannam@0:     virtual ~THistogram()				{ clear();};
cannam@0:     //@}
cannam@0: 
cannam@0:     //! \name Histogram computation, update
cannam@0:     //@{
cannam@0:     /*! Computes histogram of vector v
cannam@0:       \param v vector to compute histogram
cannam@0:       \param computeMinMax set to true if min/max of v have to be used to get the histogram limits	
cannam@0: 	
cannam@0:       This function computes the histogram of v and stores it internally.
cannam@0:       \sa Update, GeTHistogram
cannam@0:     */
cannam@0:     void compute( const std::vector<T>& v, bool computeMinMax = true);
cannam@0:     //! Update histogram with the vector v
cannam@0:     void update( const std::vector<T>& v);
cannam@0:     //! Update histogram with t
cannam@0:     void update( const T& t);
cannam@0:     //@}
cannam@0: 
cannam@0:     //! \name Resetting functions
cannam@0:     //@{
cannam@0:     //! Resize the histogram. Warning this function clear the histogram.
cannam@0:     void resize( unsigned int counters );
cannam@0:     //! Clears the histogram
cannam@0:     void clear()						{ m_counters.clear();};
cannam@0:     //@}
cannam@0: 
cannam@0:     //! \name Setters
cannam@0:     //@{
cannam@0:     /*! This function sets the minimum of the histogram spectrum. 
cannam@0:       The spectrum is not recomputed, use it with care
cannam@0:     */
cannam@0:     void setMinSpectrum( const T& min )					{	m_min = min; computeStep();};
cannam@0:     /*! This function sets the minimum of the histogram spectrum. 
cannam@0:       The spectrum is not recomputed, use it with care
cannam@0:     */
cannam@0:     void setMaxSpectrum( const T& max )					{	m_max = max; computeStep();};
cannam@0:     //@}
cannam@0:     //! \name Getters
cannam@0:     //@{
cannam@0:     //! return minimum of histogram spectrum
cannam@0:     const T& getMinSpectrum() const							{	return m_min;};
cannam@0:     //! return maximum of histogram spectrum
cannam@0:     const T& getMaxSpectrum() const							{	return m_max;};
cannam@0:     //! return step size of histogram containers
cannam@0:     TOut getStep() const									{	return m_step;};
cannam@0:     //! return number of points in histogram
cannam@0:     unsigned int getSum() const;
cannam@0:     /*! \brief returns area under the histogram 
cannam@0: 
cannam@0:     The Simpson rule is used to integrate the histogram.
cannam@0:     */
cannam@0:     TOut getArea() const;
cannam@0: 
cannam@0:     /*! \brief Computes the moments of the histogram
cannam@0: 
cannam@0:     \param data dataset
cannam@0:     \param ave mean
cannam@0:     \f[ \bar x = \frac{1}{N} \sum_{j=1}^N x_j\f]
cannam@0:     \param adev mean absolute deviation
cannam@0:     \f[ adev(X) = \frac{1}{N} \sum_{j=1}^N | x_j - \bar x |\f]
cannam@0:     \param var average deviation:
cannam@0:     \f[ \mbox{Var}(X) = \frac{1}{N-1} \sum_{j=1}^N (x_j - \bar x)^2\f]
cannam@0:     \param sdev standard deviation:
cannam@0:     \f[ \sigma(X) = \sqrt{var(\bar x) }\f]
cannam@0:     \param skew skewness
cannam@0:     \f[ \mbox{Skew}(X) = \frac{1}{N}\sum_{j=1}^N \left[ \frac{x_j - \bar x}{\sigma}\right]^3\f]
cannam@0:     \param kurt kurtosis
cannam@0:     \f[ \mbox{Kurt}(X) = \left\{ \frac{1}{N}\sum_{j=1}^N \left[ \frac{x_j - \bar x}{\sigma}\right]^4 \right\} - 3\f]
cannam@0: 
cannam@0:     */
cannam@0:     static void getMoments(const std::vector<T>& data, TOut& ave, TOut& adev, TOut& sdev, TOut& var, TOut& skew, TOut& kurt);
cannam@0: 
cannam@0:     //! return number of containers
cannam@0:     unsigned int getSize() const							{	return m_counters.size();};
cannam@0:     //! returns i-th counter
cannam@4:     unsigned int operator [] (unsigned int i) const {
cannam@4: //	ASSERT( i < m_counters.size() );
cannam@4:         return m_counters[i];
cannam@4:     };
cannam@0:     //! return the computed histogram
cannam@0:     const std::vector<unsigned int>& geTHistogram() const	{	return m_counters;};
cannam@0:     //! return the computed histogram, in TOuts
cannam@0:     std::vector<TOut> geTHistogramD() const;
cannam@0:     /*! return the normalized computed histogram 
cannam@0: 
cannam@0:     \return the histogram such that the area is equal to 1
cannam@0:     */
cannam@0:     std::vector<TOut> getNormalizedHistogram() const;
cannam@0:     //! returns left containers position
cannam@0:     std::vector<TOut> getLeftContainers() const;
cannam@0:     //! returns center containers position
cannam@0:     std::vector<TOut> getCenterContainers() const;
cannam@0:     //@}
cannam@0: protected:
cannam@0:     //! Computes the step
cannam@0:     void computeStep()								{	m_step = (TOut)(((TOut)(m_max-m_min)) / (m_counters.size()-1));};
cannam@0:     //! Data accumulators
cannam@0:     std::vector<unsigned int> m_counters;
cannam@0:     //! minimum of dataset
cannam@0:     T m_min;
cannam@0:     //! maximum of dataset
cannam@0:     T m_max;
cannam@0:     //! width of container
cannam@0:     TOut m_step;
cannam@0: };
cannam@0: 
cannam@0: template<class T, class TOut>
cannam@0: THistogram<T,TOut>::THistogram(unsigned int counters)
cannam@0:     : m_counters(counters,0), m_min(0), m_max(0), m_step(0)
cannam@0: {
cannam@0: 
cannam@0: }
cannam@0: 
cannam@0: template<class T, class TOut>
cannam@0: void THistogram<T,TOut>::resize( unsigned int counters )
cannam@0: {
cannam@0:     clear();
cannam@0: 
cannam@0:     m_counters.resize(counters,0);
cannam@0: 
cannam@0:     computeStep();
cannam@0: }
cannam@0: 
cannam@0: template<class T, class TOut>
cannam@0: void THistogram<T,TOut>::compute( const std::vector<T>& v, bool computeMinMax)
cannam@0: {
cannam@0:     using namespace std;
cannam@0:     unsigned int i;
cannam@0:     int index;
cannam@0: 
cannam@0:     if (m_counters.empty())
cannam@0: 	return;
cannam@0: 
cannam@0:     if (computeMinMax)
cannam@0:     {
cannam@0: 	m_max = m_min = v[0];
cannam@0: 	for (i=1;i<v.size();i++)
cannam@0: 	{
cannam@0: 	    m_max = std::max( m_max, v[i]);
cannam@0: 	    m_min = std::min( m_min, v[i]);
cannam@0: 	}
cannam@0:     }
cannam@0: 
cannam@0:     computeStep();
cannam@0: 
cannam@0:     for (i = 0;i < v.size() ; i++)
cannam@0:     {
cannam@0: 	index=(int) floor( ((TOut)(v[i]-m_min))/m_step ) ;
cannam@0: 
cannam@0: 	if (index >= m_counters.size() || index < 0)
cannam@0: 	    return;
cannam@0: 
cannam@0: 	m_counters[index]++;
cannam@0:     }
cannam@0: }
cannam@0: 
cannam@0: template<class T, class TOut>
cannam@0: void THistogram<T,TOut>::update( const std::vector<T>& v)
cannam@0: {
cannam@0:     if (m_counters.empty())
cannam@0: 	return;
cannam@0: 
cannam@0:     computeStep();
cannam@0: 
cannam@0:     TOut size = m_counters.size();
cannam@0: 
cannam@0:     int index;
cannam@0:     for (unsigned int i = 0;i < size ; i++)
cannam@0:     {
cannam@0: 	index = (int)floor(((TOut)(v[i]-m_min))/m_step);
cannam@0: 
cannam@0: 	if (index >= m_counters.size() || index < 0)
cannam@0: 	    return;
cannam@0: 
cannam@0: 	m_counters[index]++;
cannam@0:     }
cannam@0: }
cannam@0: 
cannam@0: template<class T, class TOut>
cannam@0: void THistogram<T,TOut>::update( const T& t)
cannam@0: {	
cannam@0:     int index=(int) floor( ((TOut)(t-m_min))/m_step ) ;
cannam@0: 
cannam@0:     if (index >= m_counters.size() || index < 0)
cannam@0: 	return;
cannam@0: 
cannam@0:     m_counters[index]++;
cannam@0: };
cannam@0: 
cannam@0: template<class T, class TOut>
cannam@0: std::vector<TOut> THistogram<T,TOut>::geTHistogramD() const
cannam@0: {
cannam@0:     std::vector<TOut> v(m_counters.size());
cannam@0:     for (unsigned int i = 0;i<m_counters.size(); i++)
cannam@0: 	v[i]=(TOut)m_counters[i];
cannam@0: 
cannam@0:     return v;
cannam@0: }
cannam@0: 
cannam@0: template <class T, class TOut>
cannam@0: std::vector<TOut> THistogram<T,TOut>::getLeftContainers() const
cannam@0: {
cannam@0:     std::vector<TOut> x( m_counters.size());
cannam@0: 
cannam@0:     for (unsigned int i = 0;i<m_counters.size(); i++)
cannam@0: 	x[i]= m_min + i*m_step;
cannam@0: 
cannam@0:     return x;
cannam@0: }
cannam@0: 
cannam@0: template <class T, class TOut>
cannam@0: std::vector<TOut> THistogram<T,TOut>::getCenterContainers() const
cannam@0: {
cannam@0:     std::vector<TOut> x( m_counters.size());
cannam@0: 
cannam@0:     for (unsigned int i = 0;i<m_counters.size(); i++)
cannam@0: 	x[i]= m_min + (i+0.5)*m_step;
cannam@0: 
cannam@0:     return x;
cannam@0: }
cannam@0: 
cannam@0: template <class T, class TOut>
cannam@0: unsigned int THistogram<T,TOut>::getSum() const
cannam@0: {
cannam@0:     unsigned int sum = 0;
cannam@0:     for (unsigned int i = 0;i<m_counters.size(); i++)
cannam@0: 	sum+=m_counters[i];
cannam@0: 
cannam@0:     return sum;
cannam@0: }
cannam@0: 
cannam@0: template <class T, class TOut>
cannam@0: TOut THistogram<T,TOut>::getArea() const
cannam@0: {
cannam@0:     const size_t n=m_counters.size();
cannam@0:     TOut area=0;
cannam@0: 
cannam@0:     if (n>6)
cannam@0:     {
cannam@0: 	area=3.0/8.0*(m_counters[0]+m_counters[n-1])
cannam@0: 	    +7.0/6.0*(m_counters[1]+m_counters[n-2])
cannam@0: 	    +23.0/24.0*(m_counters[2]+m_counters[n-3]);
cannam@0: 	for (unsigned int i=3;i<n-3;i++)
cannam@0: 	{
cannam@0: 	    area+=m_counters[i];
cannam@0: 	}
cannam@0:     }
cannam@0:     else if (n>4)
cannam@0:     {
cannam@0: 	area=5.0/12.0*(m_counters[0]+m_counters[n-1])
cannam@0: 	    +13.0/12.0*(m_counters[1]+m_counters[n-2]);
cannam@0: 	for (unsigned int i=2;i<n-2;i++)
cannam@0: 	{
cannam@0: 	    area+=m_counters[i];
cannam@0: 	}
cannam@0:     }
cannam@0:     else if (n>1)
cannam@0:     {
cannam@0: 	area=1/2.0*(m_counters[0]+m_counters[n-1]);
cannam@0: 	for (unsigned int i=1;i<n-1;i++)
cannam@0: 	{
cannam@0: 	    area+=m_counters[i];
cannam@0: 	}
cannam@0:     }
cannam@0:     else 
cannam@0: 	area=0;
cannam@0: 
cannam@0:     return area*m_step;
cannam@0: }
cannam@0: 
cannam@0: template <class T, class TOut>
cannam@0: std::vector<TOut> THistogram<T,TOut>::getNormalizedHistogram() const
cannam@0: {
cannam@0:     std::vector<TOut> normCounters( m_counters.size());
cannam@0:     TOut area = (TOut)getArea();
cannam@0: 
cannam@0:     for (unsigned int i = 0;i<m_counters.size(); i++)
cannam@0:     {
cannam@0: 	normCounters[i]= (TOut)m_counters[i]/area;
cannam@0:     }
cannam@0: 
cannam@0:     return normCounters;
cannam@0: };
cannam@0: 
cannam@0: template <class T, class TOut>
cannam@0: void THistogram<T,TOut>::getMoments(const std::vector<T>& data, TOut& ave, TOut& adev, TOut& sdev, TOut& var, TOut& skew, TOut& kurt)
cannam@0: {
cannam@0:     int j;
cannam@0:     double ep=0.0,s,p;
cannam@0:     const size_t n = data.size();
cannam@0: 
cannam@0:     if (n <= 1)
cannam@0: 	// nrerror("n must be at least 2 in moment");
cannam@0: 	return;
cannam@0: 
cannam@0:     s=0.0; // First pass to get the mean.
cannam@0:     for (j=0;j<n;j++)
cannam@0: 	s += data[j];
cannam@0: 	
cannam@0:     ave=s/(n);
cannam@0:     adev=var=skew=kurt=0.0; 
cannam@0:     /* Second pass to get the first (absolute), second,
cannam@0:        third, and fourth moments of the
cannam@0:        deviation from the mean. */
cannam@0: 
cannam@0:     for (j=0;j<n;j++) 
cannam@0:     {
cannam@0: 	adev += fabs(s=data[j]-(ave));
cannam@0: 	ep += s;
cannam@0: 	var += (p=s*s);
cannam@0: 	skew += (p *= s);
cannam@0: 	kurt += (p *= s);
cannam@0:     }
cannam@0: 
cannam@0: 
cannam@0:     adev /= n;
cannam@0:     var=(var-ep*ep/n)/(n-1); // Corrected two-pass formula.
cannam@0:     sdev=sqrt(var); // Put the pieces together according to the conventional definitions. 
cannam@0:     if (var) 
cannam@0:     {
cannam@0: 	skew /= (n*(var)*(sdev));
cannam@0: 	kurt=(kurt)/(n*(var)*(var))-3.0;
cannam@0:     } 
cannam@0:     else
cannam@0: 	//nrerror("No skew/kurtosis when variance = 0 (in moment)");
cannam@0: 	return;
cannam@0: }
cannam@0: 
cannam@0: #endif
cannam@0: