SampleCollection.h

Go to the documentation of this file.

#ifndef SAMPLECOLLECTION_H
#define SAMPLECOLLECTION_H
 
#include <memory>
#include <vector>
#include <set>
 
#include <Eigen/Core>
 
#include "MUQ/Modeling/ModPiece.h"
 
#include "MUQ/SamplingAlgorithms/SamplingState.h"
#include "MUQ/SamplingAlgorithms/SampleEstimator.h"
 
#include "MUQ/Utilities/HDF5/HDF5File.h"
 
#include <boost/property_tree/ptree.hpp>
 
namespace muq{
  namespace SamplingAlgorithms{
 
    class SamplingStateIdentity{
    public:
      SamplingStateIdentity(int blockIndIn) : blockInd(blockIndIn){};
 
      virtual ~SamplingStateIdentity() = default;
 
      Eigen::VectorXd const& operator()(SamplingState const& a);
 
      const int blockInd;
 
    private:
      Eigen::VectorXd output;
    };
 
    class ExpectedModPieceValue {
    public:
      ExpectedModPieceValue(std::shared_ptr<muq::Modeling::ModPiece> const& f, std::vector<std::string> const& metains);
 
      virtual ~ExpectedModPieceValue() = default;
 
      Eigen::VectorXd const& operator()(SamplingState const& a);
 
    private:
      std::shared_ptr<muq::Modeling::ModPiece> f;
 
      const std::vector<std::string> metains;
    };
 
    class SamplingStatePartialMoment{
    public:
      SamplingStatePartialMoment(int                    blockIndIn,
                                 int                    momentPowerIn,
                                 Eigen::VectorXd const& muIn) : blockInd(blockIndIn), momentPower(momentPowerIn), mu(muIn){};
 
      virtual ~SamplingStatePartialMoment() = default;
 
      Eigen::VectorXd const& operator()(SamplingState const& a);
 
      const int blockInd;
      const int momentPower;
      const Eigen::VectorXd& mu;
 
    private:
      Eigen::VectorXd output;
    };
 
    class SampleCollection : public SampleEstimator, public std::enable_shared_from_this<SampleCollection>{
    public:
      SampleCollection() = default;
 
      virtual ~SampleCollection() = default;
 
      virtual void Add(std::shared_ptr<SamplingState> newSamp);
 
      virtual std::shared_ptr<SamplingState> at(unsigned i);
      virtual const std::shared_ptr<SamplingState> at(unsigned i) const;
 
      virtual unsigned int size() const;
 
      virtual std::shared_ptr<SampleCollection> head(unsigned int N) const{return segment(0,N,1);}
 
      virtual std::shared_ptr<SampleCollection> tail(unsigned int N) const{return segment(size()-N,N,1);}
 
 
      virtual std::shared_ptr<SampleCollection> segment(unsigned int startInd, unsigned int length, unsigned int skipBy=1) const;
 
      virtual const std::shared_ptr<SamplingState> back() const;
 
      virtual Eigen::VectorXd CentralMoment(unsigned               order, 
                                            Eigen::VectorXd const& mean,
                                            int                    blockNum=-1) const override;
      virtual Eigen::VectorXd CentralMoment(unsigned int order, 
                                            int          blockNum=-1) const override{ return CentralMoment(order, Mean(blockNum), blockNum);};
 
      virtual Eigen::VectorXd Mean(int blockInd=-1) const override;
      virtual Eigen::MatrixXd Covariance(int blockInd=-1) const override{return SampleEstimator::Covariance(blockInd);};
      virtual Eigen::MatrixXd Covariance(Eigen::VectorXd const& mean, int blockInd=-1) const override;
 
      virtual std::vector<Eigen::MatrixXd> RunningCovariance(int blockInd=-1) const;
      virtual std::vector<Eigen::MatrixXd> RunningCovariance(Eigen::VectorXd const& mean, int blockInd=-1) const;
 
      virtual Eigen::VectorXd ESS(std::string const& method="Batch") const override{return ESS(-1,method);};
 
      virtual Eigen::VectorXd ESS(int blockDim) const override{return ESS(blockDim,"Batch");};
 
      virtual Eigen::VectorXd ESS(int blockDim, std::string const& method) const override;
 
      Eigen::VectorXd BatchESS(int blockInd=-1, int batchSize=-1, int overlap=-1) const;
 
      double MultiBatchESS(int blockInd=-1, int batchSize=-1, int overlap=-1) const;
 
      virtual Eigen::VectorXd StandardError(int                blockInd, 
                                            std::string const& method) const override;
      virtual Eigen::VectorXd StandardError(int blockInd) const override{return StandardError(blockInd,"Batch");};
      virtual Eigen::VectorXd StandardError(std::string const& method="Batch") const override{return StandardError(-1,method);};
 
      virtual Eigen::VectorXd BatchError(int blockInd=-1, int batchSize=-1, int overlap=-1) const;
 
      virtual Eigen::VectorXd MultiBatchError(int blockInd=-1, int batchSize=-1, int overlap=-1) const;
 
      virtual Eigen::MatrixXd AsMatrix(int blockDim=-1) const;
 
      static std::shared_ptr<SampleCollection> FromMatrix(Eigen::MatrixXd const& samps);
 
      virtual Eigen::VectorXd Weights() const;
 
      virtual void WriteToFile(std::string const& filename, std::string const& group = "/") const;
 
      Eigen::MatrixXd GetMeta(std::string const& name) const;
 
      std::set<std::string> ListMeta(bool requireAll=true) const;
 
      virtual Eigen::VectorXd ExpectedValue(std::shared_ptr<muq::Modeling::ModPiece> const& f, 
                                            std::vector<std::string> const& metains = std::vector<std::string>()) const override;
 
      std::vector<Eigen::VectorXd> RunningExpectedValue(std::shared_ptr<muq::Modeling::ModPiece> const& f, std::vector<std::string> const& metains = std::vector<std::string>()) const;
 
 
      virtual unsigned int BlockSize(int blockInd) const override;
 
      virtual unsigned int NumBlocks() const override;
 
      virtual Eigen::VectorXd Rhat(int                         blockDim, 
                                   unsigned int                numSegments=4, 
                                   boost::property_tree::ptree options = boost::property_tree::ptree()) const;
 
      virtual Eigen::VectorXd Rhat(unsigned int                numSegments=4, 
                                   boost::property_tree::ptree options = boost::property_tree::ptree()) const{return Rhat(-1,numSegments,options);};
                                  
    protected:
 
      std::vector<std::shared_ptr<SamplingState>> samples;
 
      static std::pair<double,double> RecursiveWeightSum(std::vector<std::shared_ptr<SamplingState>>::const_iterator start,
                                                         std::vector<std::shared_ptr<SamplingState>>::const_iterator end);
 
 
      template<typename FuncType>
      static std::pair<double,Eigen::VectorXd> RecursiveSum(std::vector<std::shared_ptr<SamplingState>>::const_iterator                         start,
                                                            std::vector<std::shared_ptr<SamplingState>>::const_iterator                         end,
                                                            FuncType& f)
      {
        int numSamps = std::distance(start,end);
        const int maxSamps = 20;
 
        // If the number of samples is small enough, we can safely add them up directly
        if(numSamps<maxSamps){
 
          Eigen::VectorXd sum = (*start)->weight * f(**start);
          double weightSum = (*start)->weight;
 
          for(auto it=start+1; it!=end; ++it){
              sum += (*it)->weight * f(**it);
              weightSum += (*it)->weight;
          }
          return std::make_pair(weightSum, sum);
 
        // Otherwise, it's more numerically stable to add things pairwise
        }else{
          int halfDist = std::floor(0.5*numSamps);
          double weight1, weight2;
          Eigen::VectorXd sum1, sum2;
          std::tie(weight1,sum1) = RecursiveSum(start, start+halfDist, f);
          std::tie(weight2,sum2) = RecursiveSum(start+halfDist, end, f);
 
          return std::make_pair(weight1+weight2, (sum1+sum2).eval());
        }
      }
 
      std::unordered_map<std::string, Eigen::MatrixXd> GetMeta() const;
 
    private:
 
      bool CreateDataset(std::shared_ptr<muq::Utilities::HDF5File> hdf5file, std::string const& dataname, int const dataSize, int const totSamps) const;
    };
  }
}
 
#endif // SAMPLECOLLECTION_H