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