SingleChainMCMC.cpp

Go to the documentation of this file.

#include "MUQ/SamplingAlgorithms/SingleChainMCMC.h"
 
#include <chrono>
 
#include "MUQ/Utilities/AnyHelpers.h"
#include "MUQ/Utilities/StringUtilities.h"
 
#include "MUQ/SamplingAlgorithms/MarkovChain.h"
 
namespace pt = boost::property_tree;
using namespace muq::Utilities;
using namespace muq::SamplingAlgorithms;
 
SingleChainMCMC::SingleChainMCMC(pt::ptree pt,
                                 std::shared_ptr<AbstractSamplingProblem> const& problem) :
                 samples(std::make_shared<MarkovChain>()),
                 QOIs(std::make_shared<MarkovChain>()),
                 printLevel(pt.get("PrintLevel",3))
{
  Setup(pt, problem);
}
 
SingleChainMCMC::SingleChainMCMC(boost::property_tree::ptree pt,
                                 std::vector<std::shared_ptr<TransitionKernel> > const& kernelsIn) :
                 samples(std::make_shared<MarkovChain>()),
                 QOIs(std::make_shared<MarkovChain>()),
                 printLevel(pt.get("PrintLevel",3))
{
  Setup(pt,kernelsIn);
}
 
#if MUQ_HAS_PARCER
SingleChainMCMC::SingleChainMCMC(pt::ptree pt,
                                 std::shared_ptr<AbstractSamplingProblem> const& problem,
                                 std::shared_ptr<parcer::Communicator> const& comm) :
                 samples(std::make_shared<MarkovChain>()),
                 QOIs(std::make_shared<MarkovChain>()),
                 printLevel(pt.get("PrintLevel",3)),
                 comm(comm)
{
  Setup(pt, problem);
}
 
SingleChainMCMC::SingleChainMCMC(pt::ptree pt,
                                 std::shared_ptr<parcer::Communicator> const& comm,
                                 std::vector<std::shared_ptr<TransitionKernel> > const& kernelsIn) :
                 samples(std::make_shared<MarkovChain>()),
                 QOIs(std::make_shared<MarkovChain>()),
                 printLevel(pt.get("PrintLevel",3)),
                 comm(comm)
{
  Setup(pt, kernelsIn);
}
 
#endif
 
 
 
void SingleChainMCMC::Setup(pt::ptree pt,
                            std::vector<std::shared_ptr<TransitionKernel>> const& kernelsIn) {
 
  assert(kernelsIn.size()>0);
 
  numSamps = pt.get<unsigned int>("NumSamples");
  burnIn = pt.get("BurnIn",0);
 
  kernels = kernelsIn;
 
  scheduler = std::make_shared<ThinScheduler>(pt);
  schedulerQOI = std::make_shared<ThinScheduler>(pt);
  assert(scheduler);
  assert(schedulerQOI);
 
#if MUQ_HAS_PARCER
 
  for(int i=0; i<kernels.size(); ++i)
    kernels.at(i)->SetCommunicator(comm);
 
#endif
 
}
 
void SingleChainMCMC::Setup(pt::ptree pt, std::shared_ptr<AbstractSamplingProblem> const& problem) {
 
  std::string kernelString = pt.get<std::string>("KernelList");
  std::vector<std::string> kernelNames = StringUtilities::Split(kernelString, ',');
 
  std::vector<std::shared_ptr<TransitionKernel>> kernelVec;
  unsigned int numBlocks = kernelNames.size();
  assert(numBlocks>0);
  kernelVec.resize(numBlocks);
 
  // Add the block id to a child tree and construct a kernel for each block
  for(int i=0; i<numBlocks; ++i) {
    boost::property_tree::ptree subTree = pt.get_child(kernelNames.at(i));
    subTree.put("BlockIndex",i);
 
    problem->AddOptions(subTree);
    kernelVec.at(i) = TransitionKernel::Construct(subTree, problem);
  }
 
  Setup(pt, kernelVec);
}
 
void SingleChainMCMC::PrintStatus(std::string prefix) const {
  PrintStatus(prefix, numSamps);
}
 
void SingleChainMCMC::PrintStatus(std::string prefix, unsigned int currInd) const
{
  std::cout << prefix << int(std::floor(double((currInd - 1) * 100) / double(numSamps))) << "% Complete" << std::endl;
  if(printLevel>1){
    for(int blockInd=0; blockInd<kernels.size(); ++blockInd){
      std::cout << prefix << "  Block " << blockInd << ":\n";
      kernels.at(blockInd)->PrintStatus(prefix + "    ");
    }
  }
}
 
std::shared_ptr<MarkovChain> SingleChainMCMC::Run(std::vector<Eigen::VectorXd> const& x0) {
  if( !x0.empty() ) { SetState(x0); }
 
  // What is the next iteration that we want to print at
  const unsigned int printIncr = std::floor(numSamps / double(10));
  unsigned int nextPrintInd = printIncr;
 
  // Run until we've received the desired number of samples
  if(printLevel>0)
    std::cout << "Starting single chain MCMC sampler..." << std::endl;
 
  while(sampNum < numSamps)
  {
    // Should we print
    if(sampNum > nextPrintInd){
      if(printLevel>0){
        PrintStatus("  ", sampNum);
      }
      nextPrintInd += printIncr;
    }
 
    Sample();
  }
 
 
  if(printLevel>0){
    PrintStatus("  ", numSamps+1);
    std::cout << "Completed in " << totalTime << " seconds." << std::endl;
  }
 
  return samples;
}
 
void SingleChainMCMC::Sample() {
  if(prevState==nullptr){
    std::stringstream msg;
    msg << "\nERROR in SingleChainMCMC::Sample!  Trying to sample chain but previous (or initial) state has not been set.\n";
    throw std::runtime_error(msg.str());
  }
 
  auto startTime = std::chrono::high_resolution_clock::now();
 
  std::vector<std::shared_ptr<SamplingState> > newStates;
 
  // Loop through each parameter block
  for(int blockInd=0; blockInd<kernels.size(); ++blockInd){
 
    // Set some metadata that might be needed by the expensive sampling problem
    prevState->meta["iteration"] = sampNum;
    prevState->meta["IsProposal"] = false;
 
    // kernel prestep
    kernels.at(blockInd)->PreStep(sampNum, prevState);
 
    // use the kernel to get the next state(s)
    newStates = kernels.at(blockInd)->Step(sampNum, prevState);
 
    // save when these samples where created
    const double now = std::chrono::duration<double>(std::chrono::high_resolution_clock::now()-startTime).count();
    for( auto it=newStates.begin(); it!=newStates.end(); ++it ) {
       (*it)->meta["time"] = now;
    }
 
    // kernel post-processing
    kernels.at(blockInd)->PostStep(sampNum, newStates);
 
    // add the new states to the SampleCollection (this also increments sampNum)
    prevState = SaveSamples(newStates, lastSavedState, sampNum);
  }
 
  auto endTime = std::chrono::high_resolution_clock::now();
  totalTime += std::chrono::duration<double>(endTime - startTime).count();
}
 
std::shared_ptr<SamplingState> SingleChainMCMC::SaveSamples(std::vector<std::shared_ptr<SamplingState> > const& newStates,
                                                            std::shared_ptr<SamplingState>& lastSavedState,
                                                            unsigned int& sampNum) const {
  for( auto it : newStates ) {
    // save the sample, if we want to
    if( ShouldSave(sampNum) ) {
      samples->Add(it);
      if (it->HasMeta("QOI")) {
        std::shared_ptr<SamplingState> qoi = AnyCast(it->meta["QOI"]);
        QOIs->Add(qoi);
      }
    }
 
    // increment the number of samples and break of we hit the max. number
    if( ++sampNum>=numSamps ) { return it; }
  }
 
  return newStates.back();
}
 
bool SingleChainMCMC::ShouldSave(unsigned int const sampNum) const { return ((sampNum>=burnIn) && (scheduler->ShouldSave(sampNum))); }
 
void SingleChainMCMC::SetState(std::vector<Eigen::VectorXd> const& x0) {
  SetState(std::make_shared<SamplingState>(x0));
}
 
void SingleChainMCMC::SetState(std::shared_ptr<SamplingState> const& x0) {
  prevState = x0;
  if(ShouldSave(0))
    samples->Add(prevState);
}
 
 
std::shared_ptr<MarkovChain> SingleChainMCMC::GetSamples() const { return samples; }
 
std::shared_ptr<MarkovChain> SingleChainMCMC::GetQOIs() const { return QOIs; }