Example
Parallelized MLMC using MCMC components

Shows how to implement parallelized Multilevel Monte Carlo type methods using the modular parallel MIMCMC stack. Allows easily runningMultilevel MC and Multilevel / Multiindex MCMC on the same model, using a common parallel architecture.

#include "MUQ/SamplingAlgorithms/SLMCMC.h"
#include "MUQ/SamplingAlgorithms/GreedyMLMCMC.h"
#include "MUQ/SamplingAlgorithms/MIMCMC.h"

#include "MUQ/Modeling/Distributions/Gaussian.h"
#include "MUQ/Modeling/Distributions/Density.h"

#include "MUQ/SamplingAlgorithms/DummyKernel.h"
#include "MUQ/SamplingAlgorithms/MHKernel.h"
#include "MUQ/SamplingAlgorithms/MHProposal.h"
#include "MUQ/SamplingAlgorithms/CrankNicolsonProposal.h"
#include "MUQ/SamplingAlgorithms/SamplingProblem.h"
#include "MUQ/SamplingAlgorithms/SubsamplingMIProposal.h"
#include "MUQ/SamplingAlgorithms/ParallelFixedSamplesMIMCMC.h"

#include <boost/property_tree/ptree.hpp>

namespace pt = boost::property_tree;
using namespace muq::Modeling;
using namespace muq::SamplingAlgorithms;
using namespace muq::Utilities;

#include "ParallelProblem.h"

int main(int argc, char **argv){

  MPI_Init(&argc, &argv);

  pt::ptree pt;

  pt.put("NumSamples_0", 1e4);
  pt.put("NumSamples_1", 5e3);
  pt.put("NumSamples_2", 1e3);
  pt.put("NumSamples_3", 5e2);
  pt.put("MLMCMC.Subsampling", 1);
  pt.put("MCMC.BurnIn", 10); // number of samples for single level
  pt.put("verbosity", 1); // show some output

/*{
  std::cout << std::endl << "*************** multilevel" << std::endl << std::endl;
  auto componentFactory = std::make_shared<MyMIComponentFactory>(pt);

  MIMCMC mimcmc(pt, componentFactory);
  mimcmc.Run();
  std::cout << "mean QOI: " << mimcmc.MeanQOI().transpose() << std::endl;

  auto index_zero = std::make_shared<MultiIndex>(1);
  index_zero->SetValue(0, 0);
  std::cout << "coarsest level mean QOI: " << mimcmc.GetBox(index_zero)->FinestChain()->GetQOIs()->Mean().transpose() << std::endl;
}


{
  std::cout << std::endl << "*************** single level reference" << std::endl << std::endl;
  auto index = componentFactory->FinestIndex();

  auto problem = componentFactory->SamplingProblem(index);
  auto proposal = componentFactory->Proposal(index, problem);

  std::vector<std::shared_ptr<TransitionKernel>> kernels(1);
  kernels[0] = std::make_shared<DummyKernel>(pt,problem,proposal);

  Eigen::VectorXd startingPoint = componentFactory->StartingPoint(index);

  auto mcmc = std::make_shared<SingleChainMCMC>(pt,kernels);
  mcmc->Run(startingPoint);
  std::cout << "mean QOI: " << mcmc->GetQOIs()->Mean().transpose() << std::endl;
}*/

{
  auto comm = std::make_shared<parcer::Communicator>();

  auto componentFactory = std::make_shared<MyMIComponentFactory>(pt);
  StaticLoadBalancingMIMCMC parallelMIMCMC (pt, componentFactory);

  if (comm->GetRank() == 0) {
    parallelMIMCMC.Run();
    Eigen::VectorXd meanQOI = parallelMIMCMC.MeanQOI();
    std::cout << "mean QOI: " << meanQOI.transpose() << std::endl;
    parallelMIMCMC.WriteToFile("samples.h5");
  }
  parallelMIMCMC.Finalize();

}

  MPI_Finalize();

  return 0;
}

Complete Code

#include "MUQ/SamplingAlgorithms/SLMCMC.h"
#include "MUQ/SamplingAlgorithms/GreedyMLMCMC.h"
#include "MUQ/SamplingAlgorithms/MIMCMC.h"

#include "MUQ/Modeling/Distributions/Gaussian.h"
#include "MUQ/Modeling/Distributions/Density.h"

#include "MUQ/SamplingAlgorithms/DummyKernel.h"
#include "MUQ/SamplingAlgorithms/MHKernel.h"
#include "MUQ/SamplingAlgorithms/MHProposal.h"
#include "MUQ/SamplingAlgorithms/CrankNicolsonProposal.h"
#include "MUQ/SamplingAlgorithms/SamplingProblem.h"
#include "MUQ/SamplingAlgorithms/SubsamplingMIProposal.h"
#include "MUQ/SamplingAlgorithms/ParallelFixedSamplesMIMCMC.h"

#include <boost/property_tree/ptree.hpp>

namespace pt = boost::property_tree;
using namespace muq::Modeling;
using namespace muq::SamplingAlgorithms;
using namespace muq::Utilities;

#include "ParallelProblem.h"

int main(int argc, char **argv){

  MPI_Init(&argc, &argv);

  pt::ptree pt;

  pt.put("NumSamples_0", 1e4);
  pt.put("NumSamples_1", 5e3);
  pt.put("NumSamples_2", 1e3);
  pt.put("NumSamples_3", 5e2);
  pt.put("MLMCMC.Subsampling", 1);
  pt.put("MCMC.BurnIn", 10); // number of samples for single level
  pt.put("verbosity", 1); // show some output

/*{
  std::cout << std::endl << "*************** multilevel" << std::endl << std::endl;
  auto componentFactory = std::make_shared<MyMIComponentFactory>(pt);

  MIMCMC mimcmc(pt, componentFactory);
  mimcmc.Run();
  std::cout << "mean QOI: " << mimcmc.MeanQOI().transpose() << std::endl;

  auto index_zero = std::make_shared<MultiIndex>(1);
  index_zero->SetValue(0, 0);
  std::cout << "coarsest level mean QOI: " << mimcmc.GetBox(index_zero)->FinestChain()->GetQOIs()->Mean().transpose() << std::endl;
}


{
  std::cout << std::endl << "*************** single level reference" << std::endl << std::endl;
  auto index = componentFactory->FinestIndex();

  auto problem = componentFactory->SamplingProblem(index);
  auto proposal = componentFactory->Proposal(index, problem);

  std::vector<std::shared_ptr<TransitionKernel>> kernels(1);
  kernels[0] = std::make_shared<DummyKernel>(pt,problem,proposal);

  Eigen::VectorXd startingPoint = componentFactory->StartingPoint(index);

  auto mcmc = std::make_shared<SingleChainMCMC>(pt,kernels);
  mcmc->Run(startingPoint);
  std::cout << "mean QOI: " << mcmc->GetQOIs()->Mean().transpose() << std::endl;
}*/

{
  auto comm = std::make_shared<parcer::Communicator>();

  auto componentFactory = std::make_shared<MyMIComponentFactory>(pt);
  StaticLoadBalancingMIMCMC parallelMIMCMC (pt, componentFactory);

  if (comm->GetRank() == 0) {
    parallelMIMCMC.Run();
    Eigen::VectorXd meanQOI = parallelMIMCMC.MeanQOI();
    std::cout << "mean QOI: " << meanQOI.transpose() << std::endl;
    parallelMIMCMC.WriteToFile("samples.h5");
  }
  parallelMIMCMC.Finalize();

}

  MPI_Finalize();

  return 0;
}