ParallelTempering.h

Go to the documentation of this file.

#ifndef PARALLELTEMPERING_H
#define PARALLELTEMPERING_H
 
#include "MUQ/SamplingAlgorithms/InferenceProblem.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 ParallelTempering
    {
 
    public:
 
 
        ParallelTempering(boost::property_tree::ptree              opts,
                          std::shared_ptr<InferenceProblem> const& problem);
 
        ParallelTempering(boost::property_tree::ptree                    opts,
                          Eigen::VectorXd                                inverseTemps,
                          std::vector<std::shared_ptr<TransitionKernel>> kernels);
 
        ParallelTempering(boost::property_tree::ptree                                 opts,
                          Eigen::VectorXd                                             inverseTemps,
                          std::vector<std::vector<std::shared_ptr<TransitionKernel>>> kernels);
 
        
 
      void SetState(std::vector<std::shared_ptr<SamplingState>> const& x0);
      void SetState(std::vector<Eigen::VectorXd> const& x0);
      void SetState(std::vector<std::vector<Eigen::VectorXd>> const& x0);
 
      template<typename... Args>
      inline void SetState(Args const&... args) {
        std::vector<Eigen::VectorXd> vec;
        SetStateRecurse(vec, args...);
      }
    
      
      double GetInverseTemp(unsigned int chainInd) const;
 
      std::vector<std::shared_ptr<TransitionKernel>> const& Kernels(unsigned int chainInd) const;
 
 
      std::shared_ptr<MarkovChain> Run();
      
      std::shared_ptr<MarkovChain> Run(Eigen::VectorXd const& x0){return Run(std::vector<Eigen::VectorXd>(numTemps,x0));}
 
      std::shared_ptr<MarkovChain> Run(std::vector<Eigen::VectorXd> const& x0);
      
      std::shared_ptr<MarkovChain> Run(std::vector<std::vector<Eigen::VectorXd>> const& initialPoints);
 
 
      //double TotalTime() { return totalTime; }
 
      void AddNumSamps(unsigned int numNewSamps){numSamps+=numNewSamps;};
 
      unsigned int NumSamps() const{return numSamps;};
 
      std::shared_ptr<MarkovChain> GetSamples() const{return chains.at(numTemps-1);}
 
      std::shared_ptr<MarkovChain> GetQOIs() const{return QOIs.at(numTemps-1);}
 
      const unsigned int numTemps;
      
    protected:
      
      Eigen::VectorXd cumulativeSwapProb; // A running sum of the swap probabilities
      Eigen::VectorXd successfulSwaps; // The number of succesfull swaps between chains i and i+1
      Eigen::VectorXd attemptedSwaps; // The number of attempted swaps between chains i and i+1
 
      
      void PrintStatus(unsigned int currInd) const{PrintStatus("",currInd);};
      void PrintStatus(std::string prefix, unsigned int currInd) const;
 
      void AdaptTemperatures();
 
      void Sample();
 
      void SwapStates();
      
      void SaveSamples(std::vector<std::vector<std::shared_ptr<SamplingState>>> const& newStates);
 
      bool ShouldSave(unsigned int chainInd, unsigned int sampNum) const;
      
      Eigen::VectorXd CollectInverseTemps() const;
 
      // Samples and quantities of interest will be stored in the MarkovChain class
      std::vector<std::shared_ptr<InferenceProblem>> problems;
      std::vector<std::shared_ptr<MarkovChain>> chains;
      std::vector<std::shared_ptr<MarkovChain>> QOIs;
 
      std::shared_ptr<SaveSchedulerBase> scheduler;
      std::shared_ptr<SaveSchedulerBase> schedulerQOI;
 
      unsigned int numSamps;
      unsigned int burnIn;
      unsigned int printLevel;
      unsigned int swapIncr;
      bool seoSwaps = false; // True if stochastic "SEO" swaps of Syed et al. 2019 should be used.  Otherwise the deterministic DEO swaps are employed.
      unsigned int adaptIncr = 2;
      unsigned int nextAdaptInd;
 
      std::vector<std::vector<std::shared_ptr<TransitionKernel>>> kernels;
    
      std::vector<std::shared_ptr<SamplingState>> prevStates;
 
    private:
 
      bool evenSwap = true;
      
      std::vector<unsigned int> sampNums;
      double totalTime = 0.0;
 
      static std::vector<std::vector<std::shared_ptr<TransitionKernel>>> StackKernels(std::vector<std::shared_ptr<TransitionKernel>> const& kerns);
      static Eigen::VectorXd ExtractTemps(boost::property_tree::ptree opts);
      static std::vector<std::vector<std::shared_ptr<TransitionKernel>>> ExtractKernels(boost::property_tree::ptree opts, std::shared_ptr<InferenceProblem> const& prob);
 
      template<typename T>
      static std::vector<std::vector<T>> StackObjects(std::vector<T> const& kerns)
      {
            std::vector<std::vector<T>> newKernels(kerns.size(), std::vector<T>(1));
            for(unsigned int i=0; i<kerns.size(); ++i)
                newKernels.at(i).at(0) = kerns.at(i);
            return newKernels;
      }
 
      static void CheckForMeta(std::shared_ptr<SamplingState> const& state);
 
 
 
    //   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 ParallelTempering 
 
  } // namespace SamplingAlgorithms
} // namespace muq
 
#endif // #ifndef SINGLECHAINMCMC_H