histogram.h

Go to the documentation of this file.

/*

  This file is a part of MeanWalker, a library which provides utilities to
  sample from probability distributions using methods like the 
  Metropilis-Hastings algorithm and Goodman-Weare algorithm.

  Copyright (C) 2012 Jori Liesenborgs

  Contact: jori.liesenborgs@gmail.com
  
  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 2 of the License, or
  (at your option) any later version.
  
  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
  
  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  
*/

#ifndef MEANWALKER_HISTOGRAM_H

#define MEANWALKER_HISTOGRAM_H

#include "meanwalkerconfig.h"
#include <errut/errorbase.h>
#include <stdint.h>
#include <iostream>
#include <vector>
#include <list>

namespace meanwalker
{

class MEANWALKER_IMPORTEXPORT Histogram : public errut::ErrorBase
{
public:
        Histogram();
        ~Histogram();

        bool init(double minValue, double maxValue, size_t numBins);

        void clear();

        void destroy();

        void process(double x);

        double getMinValue() const                                                                      { return m_minValue; }

        double getMaxValue() const                                                                      { return m_maxValue; }

        size_t getNumberOfBins() const                                                                  { return m_histogram.size(); }

        size_t getBinCount(size_t bin) const                                                            { return m_histogram[bin]; }

        double getProbabilityDensity(size_t bin) const                                                  { return ((double)m_histogram[bin]/(double)m_totalCount)/m_binSize; }
        
        double getInsideProbabilityDensity(size_t bin) const                                            { return ((double)m_histogram[bin]/(double)(m_totalCount-m_outsideCount))/m_binSize; }

        double getBinCenter(size_t bin) const                                                           { return ((double)bin+0.5)*m_binSize+m_minValue; }

        double getBinWidth() const                                                                      { return m_binSize; }

        size_t getTotalCount() const                                                                    { return m_totalCount; }

        size_t getOutsideCount() const                                                                  { return m_outsideCount; }

        size_t getInsideCount() const                                                                   { return m_totalCount - m_outsideCount; }

        void print() const;
protected:
        std::vector<size_t> m_histogram;
        double m_minValue, m_maxValue, m_binSize;
        double m_doubleBins;
        size_t m_totalCount, m_outsideCount;
};

inline void Histogram::process(double x)
{
        double binDouble = (x-m_minValue)/m_binSize;

        if (binDouble < 0 || binDouble >= m_doubleBins)
                m_outsideCount++;
        else
                m_histogram[(int)binDouble]++;

        m_totalCount++;
}

class MEANWALKER_IMPORTEXPORT DelayedHistogram : public Histogram
{
public:
        DelayedHistogram();
        void record(double value);
        bool processRecordedValues(size_t numBins);
        void clearRecordedValues()                                                                      { m_recordedValues.clear(); m_recPos = 0; }

        const double *getRecordedValues() const                                                         { return &(m_recordedValues[0]); }
        size_t recordSize() const                                                                       { return m_recPos; }
        void reserveMemory(size_t entries)                                                              { m_recordedValues.resize(entries); m_recPos = 0; }

        static void getDefaultHistogramMinMax(const double *pValues, size_t numValues, double &minValue, double &maxValue);
protected:
        virtual void getHistogramMinMax(const double *pValues, size_t numValues, double &minValue, double &maxValue)
                                                                                                        { getDefaultHistogramMinMax(pValues, numValues, minValue, maxValue); }
private:
        std::vector<double> m_recordedValues;
        size_t m_recPos;
};

inline void DelayedHistogram::record(double value)
{
        if (m_histogram.size() != 0)
                return;

        if (m_recPos < m_recordedValues.size())
        {
                m_recordedValues[m_recPos] = value;
                m_recPos++;
        }
        else
                std::cerr << "WARNING: " << m_recPos << " exceeds " << m_recordedValues.size() << std::endl;
}

} // end namespace   

#endif // MEANWALKER_HISTOGRAM_H