DRKernel.h

Go to the documentation of this file.

#ifndef DRKERNEL_H_
#define DRKERNEL_H_
 
#include "MUQ/SamplingAlgorithms/TransitionKernel.h"
#include "MUQ/SamplingAlgorithms/MCMCProposal.h"
#include "MUQ/Utilities/VectorSlice.h"
 
#include <set>
 
namespace muq {
  namespace SamplingAlgorithms {
 
    class DRKernel : public TransitionKernel {
    public:
 
      DRKernel(boost::property_tree::ptree const& pt,
               std::shared_ptr<AbstractSamplingProblem> problem);
 
      DRKernel(boost::property_tree::ptree         const& pt,
               std::shared_ptr<AbstractSamplingProblem>   problem,
               std::vector<std::shared_ptr<MCMCProposal>> proposalsIn,
               std::vector<double>                        scales);
      
      DRKernel(boost::property_tree::ptree         const& pt,
               std::shared_ptr<AbstractSamplingProblem>   problem,
               std::vector<std::shared_ptr<MCMCProposal>> proposalsIn);
 
      virtual ~DRKernel() = default;
 
      virtual inline std::vector<std::shared_ptr<MCMCProposal>> Proposals() {return proposals;};
 
      virtual void PostStep(unsigned int const t,
                            std::vector<std::shared_ptr<SamplingState>> const& state) override;
 
      virtual std::vector<std::shared_ptr<SamplingState>> Step(unsigned int const t,
                                                               std::shared_ptr<SamplingState> prevState) override;
 
      virtual void PrintStatus(std::string prefix) const override;
 
      virtual inline Eigen::VectorXd AcceptanceRates() const{return numProposalAccepts.cast<double>().array()/numProposalCalls.cast<double>().array();};
 
 
      std::shared_ptr<SamplingState> SampleProposal(unsigned int                          stage,
                                                    std::shared_ptr<SamplingState> const& state) const;
 
      double EvaluateProposal(unsigned int                          stage,
                              std::shared_ptr<SamplingState> const& x,
                              std::shared_ptr<SamplingState> const& y) const;
 
 
      std::vector<double> GetScales() const{return propScales;};
 
 
      virtual void SetBlockInd(int newBlockInd) override;
 
 
    protected:
 
      static std::vector<std::shared_ptr<MCMCProposal>> CreateProposals(boost::property_tree::ptree const& pt,
                                                                        std::shared_ptr<AbstractSamplingProblem> const& problem);
 
      static std::vector<double> CreateScales(boost::property_tree::ptree const& pt);
 
      // A vector containing one proposal for each stage
      std::vector<std::shared_ptr<MCMCProposal>> proposals;
 
      /* A set containing unique proposals (the proposals vector might have
        duplicates but with different scales, e.g., for DRAM).  This unique set
        is used to make sure that adaptive proposals that used in multiple
        stages are updated correctly.
      */
      std::set<std::shared_ptr<MCMCProposal>> uniqueProps;
 
      // Scales for each stage
      std::vector<double> propScales;
      bool isScaled = false;
 
      // store how many times each proposal has been called
      Eigen::VectorXi numProposalCalls;
 
      // store how many times each proposal has been accepted
      Eigen::VectorXi numProposalAccepts;
 
 
 
      template<typename VecType1, typename VecType2>
      double Alpha(VecType1& likelies, VecType2& proposed_points) const
      {
        // create subcontainers of input variables to use in recursion
 
        int stage = likelies.size() - 1;
 
        double a1 =  1.0;
        double a2 = 1.0;
 
        for (int k = 1; k < stage; ++k) {
          auto slice1 = muq::Utilities::GetSlice(likelies, stage, stage - k-1,-1);
          auto slice2 = muq::Utilities::GetSlice(proposed_points, stage, stage - k-1, -1);
          a2 *= (1.0 - Alpha(slice1, slice2));
 
          slice1 = muq::Utilities::GetSlice(likelies, 0, k+1);
          slice2 = muq::Utilities::GetSlice(proposed_points, 0, k+1);
          a1 *= (1.0 - Alpha(slice1, slice2)); // forward
 
          // if a2==0, there is no chance of getting accepted, so just give up and try again
          if (a2 == 0)
            return 0.0;
        }
 
        double q = 0.0;
 
        //now put in qs
 
        for (int k = 1; k <= stage; ++k) {
          q +=  QFun(muq::Utilities::GetSlice(proposed_points, stage, stage - k-1, -1));
          q -= QFun(muq::Utilities::GetSlice(proposed_points, 0, k+1)); // forward probabilities
        }
 
        return std::min<double>(1.0, exp(likelies[stage] - likelies[0] + q) * a2 / a1);
      }
 
      template<typename VecType>
      double QFun(VecType const& proposed_points) const
      {
        //use the proposal that generated this distance move
        int stage = proposed_points.size() - 1;
        return EvaluateProposal(stage - 1, proposed_points[0], proposed_points[stage]);
      }
 
 
    };
  } // namespace SamplingAlgorithms
} // namespace muq
 
#endif