SingleChainMCMC.h

Go to the documentation of this file.

#ifndef SINGLECHAINMCMC_H
#define SINGLECHAINMCMC_H
 
#include "MUQ/SamplingAlgorithms/AbstractSamplingProblem.h"
#include "MUQ/SamplingAlgorithms/TransitionKernel.h"
#include "MUQ/SamplingAlgorithms/ThinScheduler.h"
#include "MUQ/SamplingAlgorithms/MarkovChain.h"
 
#include <vector>
 
#include <boost/property_tree/ptree.hpp>
 
namespace muq{
  namespace SamplingAlgorithms{
 
    class SingleChainMCMC
    {
 
    public:
 
#if MUQ_HAS_PARCER
      SingleChainMCMC(boost::property_tree::ptree pt,
                      std::shared_ptr<parcer::Communicator> const& comm,
                      std::vector<std::shared_ptr<TransitionKernel> > const& kernelsIn);
 
      SingleChainMCMC(boost::property_tree::ptree pt,
                      std::shared_ptr<AbstractSamplingProblem> const& problem,
                      std::shared_ptr<parcer::Communicator> const& comm);
#endif
 
      SingleChainMCMC(boost::property_tree::ptree pt,
                      std::shared_ptr<AbstractSamplingProblem> const& problem);
 
      SingleChainMCMC(boost::property_tree::ptree pt,
                      std::vector<std::shared_ptr<TransitionKernel> > const& kernels);
 
      SingleChainMCMC(boost::property_tree::ptree pt, 
                      std::shared_ptr<TransitionKernel> const& kernel) : SingleChainMCMC(pt, std::vector<std::shared_ptr<TransitionKernel>>(1,kernel)){};
                      
      virtual ~SingleChainMCMC() = default;
        
 
      virtual void SetState(std::shared_ptr<SamplingState> const& x0);
      virtual void SetState(std::vector<Eigen::VectorXd> const& x0);
 
      template<typename... Args>
      inline void SetState(Args const&... args) {
        std::vector<Eigen::VectorXd> vec;
        SetStateRecurse(vec, args...);
      }
 
      virtual std::vector<std::shared_ptr<TransitionKernel>>& Kernels(){return kernels;};
 
      template<typename... Args>
      inline std::shared_ptr<MarkovChain> Run(Args const&... args) {
        std::vector<Eigen::VectorXd> vec;
        return RunRecurse(vec, args...);
      }
 
      virtual std::shared_ptr<MarkovChain> Run(std::vector<Eigen::VectorXd> const& x0);
 
      virtual double TotalTime() { return totalTime; }
 
      virtual void AddNumSamps(unsigned int numNewSamps){numSamps+=numNewSamps;};
 
      virtual unsigned int NumSamps() const{return numSamps;};
 
      void PrintStatus(std::string prefix) const;
 
      virtual std::shared_ptr<MarkovChain> GetSamples() const;
 
      std::shared_ptr<MarkovChain> GetQOIs() const;
 
#if MUQ_HAS_PARCER
      std::shared_ptr<parcer::Communicator> GetCommunicator() const;
#endif
 
    protected:
 
      virtual void Sample();
 
      std::shared_ptr<SamplingState> SaveSamples(std::vector<std::shared_ptr<SamplingState> > const& newStates,
                                                 std::shared_ptr<SamplingState>                    & lastSavedState,
                                                 unsigned int                                      & sampNum) const;
 
      bool ShouldSave(unsigned int const sampNum) const;
 
      void PrintStatus(unsigned int currInd) const{PrintStatus("",currInd);};
      void PrintStatus(std::string prefix, unsigned int currInd) const;
 
      // Samples and quantities of interest will be stored in the MarkovChain class
      std::shared_ptr<MarkovChain> samples;
      std::shared_ptr<MarkovChain> QOIs;
 
      std::shared_ptr<SaveSchedulerBase> scheduler;
      std::shared_ptr<SaveSchedulerBase> schedulerQOI;
 
      unsigned int numSamps;
      unsigned int burnIn;
      unsigned int printLevel;
 
      // A vector of transition kernels: One for each block
      std::vector<std::shared_ptr<TransitionKernel>> kernels;
 
    private:
 
#if MUQ_HAS_PARCER
      std::shared_ptr<parcer::Communicator> comm = nullptr;
#endif
 
      template<typename... Args>
      inline std::shared_ptr<MarkovChain> RunRecurse(std::vector<Eigen::VectorXd>& vec, Eigen::VectorXd const& ith, Args const&... args) {
        vec.push_back(ith);
        return RunRecurse(vec, args...);
      }
 
      inline std::shared_ptr<MarkovChain> RunRecurse(std::vector<Eigen::VectorXd>& vec, Eigen::VectorXd const& last) {
        vec.push_back(last);
        return Run(vec);
      }
 
      inline std::shared_ptr<MarkovChain> RunRecurse(std::vector<Eigen::VectorXd>& vec) {
        return Run(vec);
      }
      
      template<typename... Args>
      inline void SetStateRecurse(std::vector<Eigen::VectorXd>& vec, Eigen::VectorXd const& it, Args const&... args) {
        vec.push_back(it);
        SetStateRecurse(vec, args...);
      }
 
      inline void SetStateRecurse(std::vector<Eigen::VectorXd>& vec, Eigen::VectorXd const& last) {
        vec.push_back(last);
        SetState(vec);
      }
 
      unsigned int sampNum = 1;
      std::shared_ptr<SamplingState> prevState = nullptr;
      std::shared_ptr<SamplingState> lastSavedState = nullptr;
      double totalTime = 0.0;
 
      void Setup(boost::property_tree::ptree pt,
                 std::vector<std::shared_ptr<TransitionKernel>> const& kernelsIn);
 
      void Setup(boost::property_tree::ptree pt, std::shared_ptr<AbstractSamplingProblem> const& problem);
 
    }; // class SingleChainMCMC
  } // namespace SamplingAlgorithms
} // namespace muq
 
#endif // #ifndef SINGLECHAINMCMC_H