histogrambuildertemplate.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_HISTOGRAMBUILDERTEMPLATE_H

#define MEANWALKER_HISTOGRAMBUILDERTEMPLATE_H

#include "meanwalkerconfig.h"
#include "randomnumbergenerator.h"
#include "histogram.h"
#include <vector>

namespace meanwalker
{

template <class WalkerBase>
class MEANWALKER_IMPORTEXPORT HistogramBuilderTemplate : public WalkerBase
{
public:
        HistogramBuilderTemplate(RandomNumberGenerator *pRng) : WalkerBase(pRng)
        {
                m_init = false;
        }

        ~HistogramBuilderTemplate()
        {
        }

        bool initNormal(const double *pMinCoords, const double *pMaxCoords, size_t dim, size_t numBins)
        {
                if (m_init)
                {
                        errut::ErrorBase::setErrorString("Already initialized");
                        return false;
                }

                m_delayed = false;
                m_init = true;
                m_numBins = numBins;

                m_histograms.resize(dim);

                for (size_t i = 0 ; i < m_histograms.size() ; i++)
                {
                        m_histograms[i].setIndex(i);
                        m_histograms[i].setInstance(this);
                }

                for (size_t i = 0 ; i < m_histograms.size() ; i++)
                        m_histograms[i].init(pMinCoords[i], pMaxCoords[i], m_numBins);

                return true;
        }

        bool initDelayed(size_t dim, size_t numBins)
        {
                if (m_init)
                {
                        errut::ErrorBase::setErrorString("Already initialized");
                        return false;
                }
                m_delayed = true;
                m_init = true;
                m_numBins = numBins;

                m_histograms.resize(dim);

                for (size_t i = 0 ; i < m_histograms.size() ; i++)
                {
                        m_histograms[i].setIndex(i);
                        m_histograms[i].setInstance(this);
                }

                return true;
        }

        void destroy()
        {
                m_histograms.resize(0);
                m_init = false;
        }

        bool run(size_t initSteps, size_t finalSteps)
        {
                if (!m_init)
                {
                        errut::ErrorBase::setErrorString("Not initialized");
                        return false;
                }

                m_record = false; 
                WalkerBase::walk(initSteps); 

                if (finalSteps > 0)
                {
                        if (m_delayed)
                        {
                                for (size_t i = 0 ; i < m_histograms.size() ; i++)
                                        m_histograms[i].reserveMemory(finalSteps);
                        }

                        m_record = true; 
                        WalkerBase::walk(finalSteps); 

                        if (m_delayed)
                        {
                                m_delayed = false;
                                for (size_t i = 0 ; i < m_histograms.size() ; i++)
                                        m_histograms[i].processRecordedValues(m_numBins);

                                std::vector<const double *> recordedValues(m_histograms.size());

                                for (size_t i = 0 ; i < m_histograms.size() ; i++)
                                        recordedValues[i] = m_histograms[i].getRecordedValues();

                                onRecordedValuesProcessed(&(recordedValues[0]), m_histograms.size(), m_histograms[0].recordSize());

                                for (size_t i = 0 ; i < m_histograms.size() ; i++)
                                        m_histograms[i].clearRecordedValues();
                        }
                }

                return true;
        }

        void clearHistograms()                                                                                  
        { 
                for (size_t i = 0 ; i < m_histograms.size() ; i++) 
                        m_histograms[i].clear(); 
        }

        const Histogram &getHistogram(size_t i) const                                                           
        { 
                return m_histograms[i]; 
        }

        const std::vector<const Histogram *> getHistograms() const
        {
                std::vector<const Histogram *> h(m_histograms.size());

                for (size_t i = 0 ; i < h.size() ; i++)
                        h[i] = &(m_histograms[i]);

                return h;
        }

        void printHistograms() const
        {
                for (size_t i = 0 ; i < m_histograms.size() ; i++)
                {
                        std::cout << "# Histogram " << i << std::endl;
                        m_histograms[i].print();
                }
        }
protected:
        void onNewSample(double functionValue, const double *pCoords)
        {
                if (m_record)
                {
                        if (!m_delayed)
                        {
                                for (size_t i = 0 ; i < m_histograms.size() ; i++)
                                        m_histograms[i].process(pCoords[i]);
                        }
                        else
                        {
                                for (size_t i = 0 ; i < m_histograms.size() ; i++)
                                        m_histograms[i].record(pCoords[i]);
                        }
                }
        }

        virtual void getHistogramMinMax(size_t idx, const double *pValues, size_t numValues, double &minValue, double &maxValue)
        {
                DelayedHistogram::getDefaultHistogramMinMax(pValues, numValues, minValue, maxValue);
        }

        virtual void onRecordedValuesProcessed(const double **pValues, size_t numHistograms, size_t numRecords)
        { 
        }
private:
        class InternalDelayedHistogram : public DelayedHistogram
        {
        public:
                InternalDelayedHistogram()                                                                      { m_pInstance = 0; m_idx = 0; }

                void setInstance(HistogramBuilderTemplate *pInst)                                               { m_pInstance = pInst; }
                void setIndex(size_t idx)                                                                       { m_idx = idx; }
        protected:
                void getHistogramMinMax(const double *pValues, size_t numValues, double &minValue, double &maxValue) { m_pInstance->getHistogramMinMax(m_idx, pValues, numValues, minValue, maxValue); }
        private:
                size_t m_idx;
                HistogramBuilderTemplate *m_pInstance;
        };

        std::vector<InternalDelayedHistogram> m_histograms;
        bool m_record, m_delayed, m_init;
        size_t m_numBins;
};

} // end namespace   

#endif // MEANWALKER_HISTOGRAMBUILDERTEMPLATE_H