Task 2: notes on MCMC in dakota from MUQ

[fnrizzi]

What Dakota currently supports

the latest release of dakota: https://dakota.sandia.gov/2024/05/15/dakota-6-20/ supports: "five MCMC algorithms from MUQ: metropolis_hastings, adaptive_metropolis, delayed_rejection, dram, and mala".

Dakota files where MUQ is being included/used

fnrizzi@192 ~/D/w/m/dakota (devel)> grep -r 'MUQ/' .
./src/NonDMUQBayesCalibration.cpp:#include "MUQ/Modeling/WorkGraphPiece.h"
./src/NonDMUQBayesCalibration.cpp:#include "MUQ/Modeling/ModGraphPiece.h"
./src/NonDMUQBayesCalibration.cpp:#include "MUQ/Utilities/RandomGenerator.h"
./src/NonDMUQBayesCalibration.cpp:#include "MUQ/Utilities/AnyHelpers.h"
./src/NonDMUQBayesCalibration.cpp:#include "MUQ/SamplingAlgorithms/MHProposal.h"
./src/NonDMUQBayesCalibration.cpp:#include "MUQ/SamplingAlgorithms/AMProposal.h"
./src/NonDMUQBayesCalibration.cpp:#include "MUQ/SamplingAlgorithms/MHKernel.h"
./src/NonDMUQBayesCalibration.cpp:#include "MUQ/SamplingAlgorithms/DRKernel.h"
./src/NonDMUQBayesCalibration.cpp:#include "MUQ/SamplingAlgorithms/MALAProposal.h"
./src/unit/dakota_muq_mcmc/muq_mcmc.cpp:#include "MUQ/Modeling/Distributions/Gaussian.h"
./src/unit/dakota_muq_mcmc/muq_mcmc.cpp:#include "MUQ/Modeling/Distributions/Density.h"
./src/unit/dakota_muq_mcmc/muq_mcmc.cpp:#include "MUQ/SamplingAlgorithms/SamplingProblem.h"
./src/unit/dakota_muq_mcmc/muq_mcmc.cpp:#include "MUQ/SamplingAlgorithms/SingleChainMCMC.h"
./src/unit/dakota_muq_mcmc/muq_mcmc.cpp:#include "MUQ/SamplingAlgorithms/MCMCFactory.h"
./src/NonDMUQBayesCalibration.hpp:#include "MUQ/Modeling/WorkGraph.h"
./src/NonDMUQBayesCalibration.hpp:#include "MUQ/SamplingAlgorithms/MarkovChain.h"
./src/NonDMUQBayesCalibration.hpp:#include "MUQ/SamplingAlgorithms/SampleCollection.h"
./src/NonDMUQBayesCalibration.hpp:#include "MUQ/SamplingAlgorithms/MCMCFactory.h"
./src/NonDMUQBayesCalibration.hpp:#include "MUQ/SamplingAlgorithms/SamplingProblem.h"
./src/NonDMUQBayesCalibration.hpp:#include "MUQ/SamplingAlgorithms/SamplingState.h"
./src/NonDMUQBayesCalibration.hpp:#include "MUQ/Modeling/LinearAlgebra/IdentityOperator.h"
./src/NonDMUQBayesCalibration.hpp:#include "MUQ/Modeling/Distributions/Gaussian.h"
./src/NonDMUQBayesCalibration.hpp:#include "MUQ/Modeling/Distributions/Density.h"
./src/NonDMUQBayesCalibration.hpp:#include "MUQ/Modeling/Distributions/DensityProduct.h"

More specifically

We can see that the type of MCMC and parameters are set in the constructor:

NonDMUQBayesCalibration::NonDMUQBayesCalibration(ProblemDescDB& problem_db, Model& model):
  NonDBayesCalibration(problem_db, model),
  numBestSamples(1),
  mcmcType(probDescDB.get_string("method.nond.mcmc_type")),
  priorPropCovMult(probDescDB.get_real("method.prior_prop_cov_mult")),
  drNumStages(probDescDB.get_int("method.nond.dr_num_stages")),
  drScaleType(probDescDB.get_string("method.nond.dr_scale_type")),
  drScale(probDescDB.get_real("method.nond.dr_scale")),
  amPeriodNumSteps(probDescDB.get_int("method.nond.am_period_num_steps")),
  amStartingStep(probDescDB.get_int("method.nond.am_starting_step")),
  amScale(probDescDB.get_real("method.nond.am_scale")),
  malaStepSize(probDescDB.get_real("method.nond.mala_step_size"))

and then we have:

void NonDMUQBayesCalibration::init_bayesian_solver()
{
   //...
   if (mcmcType == "metropolis_hastings" || mcmcType == "adaptive_metropolis" || mcmcType == "mala") {
    pt.put("Kernel1.Method","MHKernel");
  }
  else if (mcmcType == "delayed_rejection" || mcmcType == "dram") {
    pt.put("Kernel1.Method","DRKernel");
  }
  // ...

  // Use the Gaussian proposal distribution to define an MCMC proposal class
  std::shared_ptr<muq::SamplingAlgorithms::MCMCProposal> proposal;
  if (mcmcType == "metropolis_hastings" || mcmcType == "delayed_rejection") {
    proposal = std::make_shared<muq::SamplingAlgorithms::MHProposal>(kernOpts.get_child("MyProposal"), problem, propDist);
  }
  else if (mcmcType == "adaptive_metropolis" || mcmcType == "dram") {
    proposal = std::make_shared<muq::SamplingAlgorithms::AMProposal>(kernOpts.get_child("MyProposal"), problem, proposalCovMatrix);
  }
  else if (mcmcType == "mala") {
    proposal = std::make_shared<muq::SamplingAlgorithms::MALAProposal>(kernOpts.get_child("MyProposal"), problem, propDist);
  }

}

[fnrizzi]

Therefore from a superficial look , it seems like extending support for say algorithm X in MUQ, involves (1) extending the paremeter choices in dakota to communicate with X and then (2) instantiate properly the components needed by the MCMC. But this only applies if the new methods need what the existing do, so maybe even if the interface needs updating.

need to make a list of the new methods needed in the SOW and make a list of what those require from the user