From 64558d05435b03d083db01ff42393c959c239369 Mon Sep 17 00:00:00 2001 From: Leah South Date: Sun, 6 Aug 2017 22:53:44 +1000 Subject: [PATCH 1/2] Parameter object and adaptation This commit is the start of making it easier to perform adaptation in SMC. The main changes here are: - Adding a template parameter for the algorithm parameters to the sampler object - Creating a base class for adaptation - Changing the MCMC function to give a boolean return. --- ChangeLog | 22 ++++- R/RcppExports.R | 16 ++-- inst/include/LinReg.h | 2 +- inst/include/LinReg_LA.h | 2 +- inst/include/algParam.h | 63 ++++++++++++++ inst/include/moveset.h | 34 ++++---- inst/include/sampler.h | 174 ++++++++++++++++++++++++++++----------- src/LinReg.cpp | 37 ++++----- src/LinReg_LA.cpp | 28 +++---- src/RcppExports.cpp | 41 +++------ 10 files changed, 278 insertions(+), 141 deletions(-) create mode 100644 inst/include/algParam.h diff --git a/ChangeLog b/ChangeLog index 5e436d0..20c15e4 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,4 +1,24 @@ -2017-08-04 +2017-08-06 Leah South + + * inst/include/algParam.h: A base class to contain + additional algorithm parameters and virtual functions to + adapt them. + * inst/include/moveset.h: Changing the MCMC step to perform + multiple iterations within the library rather than on user end. + * inst/include/sampler.h: Adding an additional template + parameter for additional parameters and allowing for + optional adaptation at multiple points. + + * src/LinReg.cpp: Boolean return for MCMC function and + setting MCMC repeats separately. + * src/LinReg.cpp: Idem. + * inst/include/LinReg.h: Idem. + * inst/include/LinReg_LA.h: Idem. + + * src/RcppExports.cpp: Regenerated. + * R/RcppExports.R: Idem. + +2017-08-04 Adam M. Johansen * R/LinRegLA.R: Suppressing redundant call to data(). * R/LinReg.R: Idem. diff --git a/R/RcppExports.R b/R/RcppExports.R index b894765..a5cc2aa 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -2,26 +2,26 @@ # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 blockpfGaussianOpt_impl <- function(data, part, lag) { - .Call(`_RcppSMC_blockpfGaussianOpt_impl`, data, part, lag) + .Call('RcppSMC_blockpfGaussianOpt_impl', PACKAGE = 'RcppSMC', data, part, lag) } -LinRegLA_impl <- function(Data, intemps, lNumber) { - .Call(`_RcppSMC_LinRegLA_impl`, Data, intemps, lNumber) +LinReg_impl <- function(Data, lNumber) { + .Call('RcppSMC_LinReg_impl', PACKAGE = 'RcppSMC', Data, lNumber) } -LinReg_impl <- function(Data, lNumber) { - .Call(`_RcppSMC_LinReg_impl`, Data, lNumber) +LinRegLA_impl <- function(Data, intemps, lNumber) { + .Call('RcppSMC_LinRegLA_impl', PACKAGE = 'RcppSMC', Data, intemps, lNumber) } nonLinPMMH_impl <- function(data, lNumber, lMCMCits) { - .Call(`_RcppSMC_nonLinPMMH_impl`, data, lNumber, lMCMCits) + .Call('RcppSMC_nonLinPMMH_impl', PACKAGE = 'RcppSMC', data, lNumber, lMCMCits) } pfLineartBS_impl <- function(data, part, usef, fun) { - .Call(`_RcppSMC_pfLineartBS_impl`, data, part, usef, fun) + .Call('RcppSMC_pfLineartBS_impl', PACKAGE = 'RcppSMC', data, part, usef, fun) } pfNonlinBS_impl <- function(data, part) { - .Call(`_RcppSMC_pfNonlinBS_impl`, data, part) + .Call('RcppSMC_pfNonlinBS_impl', PACKAGE = 'RcppSMC', data, part) } diff --git a/inst/include/LinReg.h b/inst/include/LinReg.h index b9df5ed..3de683d 100644 --- a/inst/include/LinReg.h +++ b/inst/include/LinReg.h @@ -43,5 +43,5 @@ namespace LinReg { double logPosterior(long lTime, const rad_state & value); void fInitialise(rad_state & value, double & logweight); void fMove(long lTime, rad_state & value, double & logweight); - int fMCMC(long lTime, rad_state & value, double & logweight); + bool fMCMC(long lTime, rad_state & value, double & logweight); } diff --git a/inst/include/LinReg_LA.h b/inst/include/LinReg_LA.h index 783f24f..05df115 100644 --- a/inst/include/LinReg_LA.h +++ b/inst/include/LinReg_LA.h @@ -48,7 +48,7 @@ namespace LinReg_LA { void fInitialise(rad_state & value, double & logweight); void fMove(long lTime, rad_state & value, double & logweight); - int fMCMC(long lTime, rad_state & value, double & logweight); + bool fMCMC(long lTime, rad_state & value, double & logweight); double integrand_ps(long,const rad_state &, void *); double width_ps(long, void *); diff --git a/inst/include/algParam.h b/inst/include/algParam.h new file mode 100644 index 0000000..a011953 --- /dev/null +++ b/inst/include/algParam.h @@ -0,0 +1,63 @@ +// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; indent-tabs-mode: nil; -*- +// +// algParam: A class that holds all of the algorithm parameters that can be adapted +// +// Copyright (C) 2017 Dirk Eddelbuettel, Adam Johansen and Leah South +// +// This file is part of RcppSMC. +// +// RcppSMC is free software: you can redistribute it and/or modify it +// under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 2 of the License, or +// (at your option) any later version. +// +// RcppSMC is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with RInside. If not, see . + + +//! \file +//! \brief A base class containing algorithm parameters and virtual +//! functions to adapt them. +//! + +#ifndef __SMC_ALGPARAM_H +#define __SMC_ALGPARAM_H 1.0 + +#include + + +namespace smc { + + /// A base class which contains the algorithm parameters and virtual functions to adapt them. + template class algParam { + + protected: + Params param; + + public: + + /// Free the workspace allocated for the algorithm parameters. + virtual ~algParam() { + } + + /// Holder function for updates to be done before the move step. + virtual void updateForMove(const population & pop) {} + + /// Holder function for updates to be done before the MCMC step. + virtual void updateForMCMC(const population & pop, double acceptProb, int nResampled, int & nRepeats) {} + + /// Holder function for updates to be done at the end of each iteration. + virtual void updateEnd(const population & pop) {} + + /// Returns the parameters. + const Params & GetParams(void) const {return param;} + }; +} + + +#endif diff --git a/inst/include/moveset.h b/inst/include/moveset.h index d25913d..7ba0f90 100644 --- a/inst/include/moveset.h +++ b/inst/include/moveset.h @@ -46,8 +46,8 @@ namespace smc { long (*pfMoveSelect)(long , const Space &, const double &); ///The functions which perform actual moves on a single particle. void (**pfMoves)(long, Space &, double &); - ///A Markov Chain Monte Carlo move for a single particle. - int (*pfMCMC)(long, Space &,double &); + ///One iteration of a Markov Chain Monte Carlo move for a single particle. + bool (*pfMCMC)(long, Space &,double &); public: ///Create a completely unspecified moveset @@ -55,16 +55,16 @@ namespace smc { ///Create a reduced moveset with a single move moveset(void (*pfInit)(Space &, double &), void (*pfNewMoves)(long, Space &,double &), - int (*pfNewMCMC)(long,Space &,double &)); + bool (*pfNewMCMC)(long,Space &,double &)); ///Create a fully specified moveset moveset(void (*pfInit)(Space &, double &),long (*pfMoveSelector)(long , const Space &, const double &), long nMoves, void (**pfNewMoves)(long, Space &,double &), - int (*pfNewMCMC)(long,Space &, double &)); + bool (*pfNewMCMC)(long,Space &, double &)); ///Initialise the population of particles void DoInit(population & pFrom, long N); ///Perform an MCMC move on the particles - int DoMCMC(long lTime, population & pFrom, long N); + bool DoMCMC(long lTime, population & pFrom, long N, int nRepeats, int & nAccepted); ///Select an appropriate move at time lTime and apply it to pFrom void DoMove(long lTime, population & pFrom,long N); @@ -78,7 +78,7 @@ namespace smc { /// \brief Set the MCMC function /// \param pfNewMCMC The function which performs an MCMC move - void SetMCMCFunction(int (*pfNewMCMC)(long,Space &,double &)) {pfMCMC = pfNewMCMC;} + void SetMCMCFunction(bool (*pfNewMCMC)(long,Space &,double &)) {pfMCMC = pfNewMCMC;} /// \brief Set the move selection function /// \param pfMoveSelectNew returns the index of move to perform at the specified time given the specified particle @@ -114,7 +114,7 @@ namespace smc { template moveset::moveset(void (*pfInit)(Space &, double &), void (*pfNewMoves)(long, Space &,double &), - int (*pfNewMCMC)(long,Space &,double &)) + bool (*pfNewMCMC)(long,Space &,double &)) { SetInitialisor(pfInit); SetMoveSelectionFunction(NULL); @@ -133,7 +133,7 @@ namespace smc { template moveset::moveset(void (*pfInit)(Space &, double &),long (*pfMoveSelector)(long ,const Space &, const double &), long nMoves, void (**pfNewMoves)(long, Space &,double &), - int (*pfNewMCMC)(long,Space &,double &)) + bool (*pfNewMCMC)(long,Space &,double &)) { SetInitialisor(pfInit); SetMoveSelectionFunction(pfMoveSelector); @@ -157,16 +157,20 @@ namespace smc { } template - int moveset::DoMCMC(long lTime, population & pFrom, long N) - { if(pfMCMC) { - int count = 0; - for (long i=0; i::DoMCMC(long lTime, population & pFrom, long N, int nRepeats, int & nAccepted) + { + if(pfMCMC && nRepeats>0) { + nAccepted = 0; + for (int j=0; j + template class sampler { private: @@ -91,11 +97,20 @@ namespace smc { population pPopulation; ///The set of moves available. moveset Moves; + ///An object for tracking and adapting additional algorithm parameters. + algParam* pAlgParams; + ///A flag to track whether the adaptation object needs to be included in this destructor. + bool pAlgBelong; ///The number of MCMC moves which have been accepted during this iteration int nAccepted; ///A flag which tracks whether the ensemble was resampled during this iteration int nResampled; + ///The number of MCMC repeats to be performed. The default is 1 if an MCMC step is supplied. + int nRepeats; + ///The proportion of accepted MCMC proposals in the most recent MCMC step, with + /// a default of -1 if no MCMC steps have been performed. + double acceptProb; ///An estimate of the log normalising constant ratio over the entire path. double dlogNCPath; @@ -110,18 +125,28 @@ namespace smc { public: ///Create an particle system containing lSize uninitialised particles with the specified mode. sampler(long lSize, HistoryType::Enum htHistoryMode); + ///Create a particle system containing lSize uninitialised particles with the specified history mode and adaptation object. + sampler(long lSize, HistoryType::Enum htHistoryMode, algParam* adaptSet); ///Dispose of a sampler. ~sampler(); ///Calculates and Returns the Effective Sample Size. double GetESS(void) const; /// Returns the Effective Sample Size of the specified particle generation. double GetESS(long lGeneration) const; + ///Returns the number of accepted proposals from the most recent MCMC iteration + int GetAccepted(void) const {return nAccepted;} + ///Returns a flag for whether the ensemble was resampled during the most recent iteration + int GetResampled(void) const {return nResampled;} + ///Returns the number of MCMC repeats used in the most recent iteration + int GetMcmcRepeats(void) const {return nRepeats;} ///Returns the History of the particle system const std::vector > & GetHistory(void) const { return History; } ///Returns the number of particles within the system. long GetNumber(void) const {return N;} ///Returns the number of evolution times stored in the history. long GetHistoryLength(void) const {return History.size();} + ///Returns a pointer to the additional algorithm parameter and adaptation object. + algParam * GetAlgParams(void) const { return pAlgParams; } ///Return the value of particle n const Space & GetParticleValueN(long n) const { return pPopulation.GetValueN(n); } ///Return the logarithmic unnormalized weight of particle n @@ -164,6 +189,8 @@ namespace smc { void SetMoveSet(moveset& pNewMoveset) { Moves = pNewMoveset; } ///Set Resampling Parameters void SetResampleParams(ResampleType::Enum rtMode, double dThreshold); + ///Sets the number of MCMC repeats + void SetMcmcRepeats(int reps) {nRepeats = reps;} ///Dump a specified particle to the specified output stream in a human readable form std::ostream & StreamParticle(std::ostream & os, long n) const; ///Dump the entire particle set to the specified output stream in a human readable form @@ -175,9 +202,9 @@ namespace smc { private: ///Duplication of smc::sampler is not currently permitted. - sampler(const sampler & sFrom); + sampler(const sampler & sFrom); ///Duplication of smc::sampler is not currently permitted. - sampler & operator=(const sampler & sFrom); + sampler & operator=(const sampler & sFrom); protected: ///Returns the crude normalising constant ratio estimate implied by the weights. @@ -192,11 +219,16 @@ namespace smc { /// \param lSize The number of particles present in the ensemble (at time 0 if this is a variable quantity) /// \param htHM The history mode to use: set this to HistoryType::RAM to store the whole history of the system and SMC_HISTORY_NONE to avoid doing so. /// \tparam Space The class used to represent a point in the sample space. - template - sampler::sampler(long lSize, HistoryType::Enum htHM) + /// \tparam Params (optional) The class used for any additional parameters. + template + sampler::sampler(long lSize, HistoryType::Enum htHM) { N = lSize; uRSCount = arma::zeros >(static_cast(N)); + + pAlgParams = new algParam; + pAlgBelong = 1; + nRepeats = 1; //Some workable defaults. htHistoryMode = htHM; @@ -204,20 +236,46 @@ namespace smc { dResampleThreshold = 0.5 * N; } - template - sampler::~sampler() + /// The constructor prepares a sampler for use but does not assign any moves to the moveset, initialise the particles + /// or otherwise perform any sampling related tasks. Its main function is to allocate a region of memory in which to + /// store the particle set. + /// + /// \param lSize The number of particles present in the ensemble (at time 0 if this is a variable quantity) + /// \param htHM The history mode to use: set this to HistoryType::RAM to store the whole history of the system and SMC_HISTORY_NONE to avoid doing so. + /// \param adaptSet The class derived from algParam for parameter adaptation. + /// \tparam Space The class used to represent a point in the sample space. + /// \tparam Params The class used for additional algorithm parameters. + template + sampler::sampler(long lSize, HistoryType::Enum htHM, algParam* adaptSet) { + N = lSize; + pAlgParams = adaptSet; + pAlgBelong = 0; + nRepeats = 1; + uRSCount = arma::zeros >(static_cast(N)); + + //Some workable defaults. + htHistoryMode = htHM; + rtResampleMode = ResampleType::STRATIFIED; + dResampleThreshold = 0.5 * N; } - template - double sampler::GetESS(void) const + template + sampler::~sampler() + { + if(pAlgBelong) + delete pAlgParams; + } + + template + double sampler::GetESS(void) const { return expl(2*stableLogSumWeights(pPopulation.GetLogWeight())-stableLogSumWeights(2.0*pPopulation.GetLogWeight())); } - template - double sampler::GetESS(long lGeneration) const + template + double sampler::GetESS(long lGeneration) const { typename std::vector >::const_iterator it = History.begin(); std::advance(it,lGeneration); @@ -228,11 +286,12 @@ namespace smc { /// particle in the ensemble. /// /// Note that the initialisation function must be specified before calling this function. - template - void sampler::Initialise(void) + template + void sampler::Initialise(void) { T = 0; dlogNCPath = 0.0; + acceptProb = -1; //Set the initial values and log weights of the particles std::vector InitVal(N); @@ -255,17 +314,26 @@ namespace smc { double ESS = GetESS(); if(ESS < dResampleThreshold) { nResampled = 1; + pAlgParams->updateForMCMC(pPopulation,acceptProb,nResampled,nRepeats); Resample(rtResampleMode); } - else - nResampled = 0; + else { + nResampled = 0; + pAlgParams->updateForMCMC(pPopulation,acceptProb,nResampled,nRepeats); + } //A possible MCMC step should be included here. - nAccepted += Moves.DoMCMC(0,pPopulation, N); + bool didMCMC = Moves.DoMCMC(0,pPopulation, N, nRepeats, nAccepted); + if (didMCMC){ + acceptProb = static_cast(nAccepted)/(static_cast(N)*static_cast(nRepeats)); + } //Normalise the weights pPopulation.SetLogWeight(pPopulation.GetLogWeight() - CalcLogNC()); + //Perform any final updates to the additional algorithm parameters. + pAlgParams->updateEnd(pPopulation); + //Finally, the current particle set should be appended to the historical process. if(htHistoryMode != HistoryType::NONE) { History.clear(); @@ -284,8 +352,8 @@ namespace smc { /// \param pIntegrand The function to integrate with respect to the particle set /// \param pAuxiliary A pointer to any auxiliary data which should be passed to the function - template - double sampler::Integrate(double(*pIntegrand)(const Space&,void*), void * pAuxiliary) + template + double sampler::Integrate(double(*pIntegrand)(const Space&,void*), void * pAuxiliary) { long double rValue = 0; long double wSum = expl(stableLogSumWeights(pPopulation.GetLogWeight())); @@ -320,14 +388,15 @@ namespace smc { /// \param pWidth The function which returns the width of the path sampling grid at the specified evolution time. The final argument is always pAuxiliary /// \param pAuxiliary A pointer to auxiliary data to pass to both of the above functions /// \tparam Space The class used to represent a point in the sample space. + /// \tparam Params (optional) The class used for any additional parameters. /// /// The PStype parameter should be set to one of the following: /// -# PathSamplingType::RECTANGLE to use the rectangle rule for integration. /// -# PathSamplingType::TRAPEZOID1 to use the trapezoidal rule for integration. /// -# PathSamplingType::TRAPEZOID2 to use the trapezoidal rule for integration with a second order correction. - template - double sampler::IntegratePathSampling(PathSamplingType::Enum PStype, double (*pIntegrand)(long,const Space &, void*), double (*pWidth)(long,void*), void* pAuxiliary) + template + double sampler::IntegratePathSampling(PathSamplingType::Enum PStype, double (*pIntegrand)(long,const Space &, void*), double (*pWidth)(long,void*), void* pAuxiliary) { if(htHistoryMode == HistoryType::NONE) throw SMC_EXCEPTION(SMCX_MISSING_HISTORY, "The path sampling integral cannot be computed as the history of the system was not stored."); @@ -401,15 +470,15 @@ namespace smc { /// -# performs a mcmc step if required /// -# appends the current particle set to the history if desired /// -# increments the current evolution time - template - void sampler::Iterate(void) + template + void sampler::Iterate(void) { IterateEss(); return; } - template - void sampler::IterateBack(void) + template + void sampler::IterateBack(void) { if(htHistoryMode == HistoryType::NONE) throw SMC_EXCEPTION(SMCX_MISSING_HISTORY, "An attempt to undo an iteration was made; unforunately, the system history has not been stored."); @@ -424,10 +493,10 @@ namespace smc { return; } - template - double sampler::IterateEss(void) + template + double sampler::IterateEss(void) { - nAccepted = 0; + pAlgParams->updateForMove(pPopulation); //Move the particle set. MoveParticles(); @@ -444,17 +513,26 @@ namespace smc { double ESS = GetESS(); if(ESS < dResampleThreshold) { nResampled = 1; + pAlgParams->updateForMCMC(pPopulation,acceptProb,nResampled,nRepeats); Resample(rtResampleMode); } - else - nResampled = 0; + else { + nResampled = 0; + pAlgParams->updateForMCMC(pPopulation,acceptProb,nResampled,nRepeats); + } //A possible MCMC step should be included here. - nAccepted += Moves.DoMCMC(T+1,pPopulation, N); + bool didMCMC = Moves.DoMCMC(T+1,pPopulation, N, nRepeats, nAccepted); + if (didMCMC){ + acceptProb = static_cast(nAccepted)/(static_cast(N)*static_cast(nRepeats)); + } //Normalise the weights pPopulation.SetLogWeight(pPopulation.GetLogWeight() - CalcLogNC()); + //Perform any final updates to the additional algorithm parameters. + pAlgParams->updateEnd(pPopulation); + //Finally, the current particle set should be appended to the historical process. if(htHistoryMode != HistoryType::NONE){ historyelement histel; @@ -468,21 +546,21 @@ namespace smc { return ESS; } - template - void sampler::IterateUntil(long lTerminate) + template + void sampler::IterateUntil(long lTerminate) { while(T < lTerminate) Iterate(); } - template - void sampler::MoveParticles(void) + template + void sampler::MoveParticles(void) { Moves.DoMove(T+1,pPopulation, N); } - template - void sampler::Resample(ResampleType::Enum lMode) + template + void sampler::Resample(ResampleType::Enum lMode) { //Resampling is done in place. int uMultinomialCount; @@ -587,8 +665,8 @@ namespace smc { /// The dThreshold parameter can be set to a value in the range [0,1) corresponding to a fraction of the size of /// the particle set or it may be set to an integer corresponding to an actual effective sample size. - template - void sampler::SetResampleParams(ResampleType::Enum rtMode, double dThreshold) + template + void sampler::SetResampleParams(ResampleType::Enum rtMode, double dThreshold) { rtResampleMode = rtMode; if(dThreshold < 1) @@ -601,8 +679,8 @@ namespace smc { /// /// \param os The output stream to which the display should be made. /// \param n The index of the particle of interest - template - std::ostream & sampler::StreamParticle(std::ostream & os, long n) const + template + std::ostream & sampler::StreamParticle(std::ostream & os, long n) const { os << pPopulation.GetValueN(n) << "," << pPopulation.GetWeightN(n) << std::endl; return os; @@ -611,8 +689,8 @@ namespace smc { /// Produce a human-readable display of the current particle values and log weights. /// /// \param os The output stream to which the display should be made. - template - std::ostream & sampler::StreamParticles(std::ostream & os) const + template + std::ostream & sampler::StreamParticles(std::ostream & os) const { os << pPopulation << std::endl; return os; @@ -622,8 +700,8 @@ namespace smc { /// the number of moves accepted at each time instant. /// /// \param os The output stream to send the data to. - template - void sampler:: OstreamMCMCRecordToStream(std::ostream &os) const + template + void sampler:: OstreamMCMCRecordToStream(std::ostream &os) const { os << "Accepted MCMC proposals history:" << std::endl; os << "======================" << std::endl; @@ -635,8 +713,8 @@ namespace smc { /// the value 1 for those time instances when resampling occured and 0 otherwise. /// /// \param os The output stream to send the data to. - template - void sampler:: OstreamResamplingRecordToStream(std::ostream &os) const + template + void sampler:: OstreamResamplingRecordToStream(std::ostream &os) const { os << "Resampling history:" << std::endl; os << "======================" << std::endl; @@ -658,8 +736,8 @@ namespace std { /// \param os The output stream to which the display should be made. /// \param s The sampler which is to be displayed. - template - std::ostream & operator<< (std::ostream & os, smc::sampler & s) + template + std::ostream & operator<< (std::ostream & os, smc::sampler & s) { os << "Sampler Configuration:" << std::endl; os << "======================" << std::endl; diff --git a/src/LinReg.cpp b/src/LinReg.cpp index 5ec5c7b..bbfdb8f 100644 --- a/src/LinReg.cpp +++ b/src/LinReg.cpp @@ -61,6 +61,7 @@ Rcpp::List LinReg_impl(arma::mat Data, unsigned long lNumber) { Sampler.SetResampleParams(ResampleType::MULTINOMIAL, 0.5); Sampler.SetMoveSet(Moveset); + Sampler.SetMcmcRepeats(10); Sampler.Initialise(); Sampler.IterateUntil(lIterates-1); @@ -149,30 +150,22 @@ namespace LinReg { /// \param lTime The sampler iteration. /// \param value Reference to the current particle value /// \param logweight Reference to the log weight of the particle being moved - int fMCMC(long lTime, rad_state & value, double & logweight) + bool fMCMC(long lTime, rad_state & value, double & logweight) { - double MH_ratio; - double dRand; - int count = 0; - - rad_state value_prop; double logposterior_curr = logPosterior(lTime, value); - double logposterior_prop; - - for (unsigned int j=0; j<10; j++){ - value_prop.theta = value.theta + cholCovRW*Rcpp::as(Rcpp::rnorm(3)); - - logposterior_prop = logPosterior(lTime, value_prop); - - MH_ratio = exp(logposterior_prop - logposterior_curr); - dRand = R::runif(0,1); - - if (MH_ratio>dRand){ - value = value_prop; - logposterior_curr = logposterior_prop; - count++; - } + + rad_state value_prop; + value_prop.theta = value.theta + cholCovRW*Rcpp::as(Rcpp::rnorm(3)); + + double logposterior_prop = logPosterior(lTime, value_prop); + + double MH_ratio = exp(logposterior_prop - logposterior_curr); + + if (MH_ratio>R::runif(0,1)){ + value = value_prop; + logposterior_curr = logposterior_prop; + return TRUE; } - return count; + return FALSE; } } diff --git a/src/LinReg_LA.cpp b/src/LinReg_LA.cpp index c00e228..600c665 100644 --- a/src/LinReg_LA.cpp +++ b/src/LinReg_LA.cpp @@ -62,6 +62,7 @@ Rcpp::List LinRegLA_impl(arma::mat Data, arma::vec intemps, unsigned long lNumbe Sampler.SetResampleParams(ResampleType::SYSTEMATIC, 0.5); Sampler.SetMoveSet(Moveset); + Sampler.SetMcmcRepeats(10); Sampler.Initialise(); arma::cube theta(lNumber,3,lTemps); @@ -166,26 +167,19 @@ namespace LinReg_LA { ///\param lTime The sampler iteration. ///\param value Reference to the value of the particle being moved ///\param logweight Reference to the log weight of the particle being moved - int fMCMC(long lTime, rad_state & value, double & logweight) + bool fMCMC(long lTime, rad_state & value, double & logweight) { - double MH_ratio; - double dRand; - int count = 0; rad_state value_prop; + value_prop.theta = value.theta + cholCovRW*Rcpp::as(Rcpp::rnorm(3)); + value_prop.loglike = logLikelihood(value_prop); + value_prop.logprior = logPrior(value_prop); - for (unsigned int j=0; j<10; j++){ - value_prop.theta = value.theta + cholCovRW*Rcpp::as(Rcpp::rnorm(3)); - value_prop.loglike = logLikelihood(value_prop); - value_prop.logprior = logPrior(value_prop); - - MH_ratio = exp(temps(lTime)*(value_prop.loglike - value.loglike) + value_prop.logprior - value.logprior); - dRand = R::runif(0,1); - - if (MH_ratio>dRand){ - value = value_prop; - count++; - } + double MH_ratio = exp(temps(lTime)*(value_prop.loglike - value.loglike) + value_prop.logprior - value.logprior); + + if (MH_ratio>R::runif(0,1)){ + value = value_prop; + return TRUE; } - return count; + return FALSE; } } \ No newline at end of file diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index dbde803..6821e70 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -8,7 +8,7 @@ using namespace Rcpp; // blockpfGaussianOpt_impl Rcpp::List blockpfGaussianOpt_impl(arma::vec data, long part, long lag); -RcppExport SEXP _RcppSMC_blockpfGaussianOpt_impl(SEXP dataSEXP, SEXP partSEXP, SEXP lagSEXP) { +RcppExport SEXP RcppSMC_blockpfGaussianOpt_impl(SEXP dataSEXP, SEXP partSEXP, SEXP lagSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -19,34 +19,34 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } -// LinRegLA_impl -Rcpp::List LinRegLA_impl(arma::mat Data, arma::vec intemps, unsigned long lNumber); -RcppExport SEXP _RcppSMC_LinRegLA_impl(SEXP DataSEXP, SEXP intempsSEXP, SEXP lNumberSEXP) { +// LinReg_impl +Rcpp::List LinReg_impl(arma::mat Data, unsigned long lNumber); +RcppExport SEXP RcppSMC_LinReg_impl(SEXP DataSEXP, SEXP lNumberSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< arma::mat >::type Data(DataSEXP); - Rcpp::traits::input_parameter< arma::vec >::type intemps(intempsSEXP); Rcpp::traits::input_parameter< unsigned long >::type lNumber(lNumberSEXP); - rcpp_result_gen = Rcpp::wrap(LinRegLA_impl(Data, intemps, lNumber)); + rcpp_result_gen = Rcpp::wrap(LinReg_impl(Data, lNumber)); return rcpp_result_gen; END_RCPP } -// LinReg_impl -Rcpp::List LinReg_impl(arma::mat Data, unsigned long lNumber); -RcppExport SEXP _RcppSMC_LinReg_impl(SEXP DataSEXP, SEXP lNumberSEXP) { +// LinRegLA_impl +Rcpp::List LinRegLA_impl(arma::mat Data, arma::vec intemps, unsigned long lNumber); +RcppExport SEXP RcppSMC_LinRegLA_impl(SEXP DataSEXP, SEXP intempsSEXP, SEXP lNumberSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< arma::mat >::type Data(DataSEXP); + Rcpp::traits::input_parameter< arma::vec >::type intemps(intempsSEXP); Rcpp::traits::input_parameter< unsigned long >::type lNumber(lNumberSEXP); - rcpp_result_gen = Rcpp::wrap(LinReg_impl(Data, lNumber)); + rcpp_result_gen = Rcpp::wrap(LinRegLA_impl(Data, intemps, lNumber)); return rcpp_result_gen; END_RCPP } // nonLinPMMH_impl Rcpp::DataFrame nonLinPMMH_impl(arma::vec data, unsigned long lNumber, unsigned long lMCMCits); -RcppExport SEXP _RcppSMC_nonLinPMMH_impl(SEXP dataSEXP, SEXP lNumberSEXP, SEXP lMCMCitsSEXP) { +RcppExport SEXP RcppSMC_nonLinPMMH_impl(SEXP dataSEXP, SEXP lNumberSEXP, SEXP lMCMCitsSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -59,7 +59,7 @@ END_RCPP } // pfLineartBS_impl Rcpp::DataFrame pfLineartBS_impl(arma::mat data, unsigned long part, bool usef, Rcpp::Function fun); -RcppExport SEXP _RcppSMC_pfLineartBS_impl(SEXP dataSEXP, SEXP partSEXP, SEXP usefSEXP, SEXP funSEXP) { +RcppExport SEXP RcppSMC_pfLineartBS_impl(SEXP dataSEXP, SEXP partSEXP, SEXP usefSEXP, SEXP funSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -73,7 +73,7 @@ END_RCPP } // pfNonlinBS_impl Rcpp::List pfNonlinBS_impl(arma::vec data, long part); -RcppExport SEXP _RcppSMC_pfNonlinBS_impl(SEXP dataSEXP, SEXP partSEXP) { +RcppExport SEXP RcppSMC_pfNonlinBS_impl(SEXP dataSEXP, SEXP partSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -83,18 +83,3 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } - -static const R_CallMethodDef CallEntries[] = { - {"_RcppSMC_blockpfGaussianOpt_impl", (DL_FUNC) &_RcppSMC_blockpfGaussianOpt_impl, 3}, - {"_RcppSMC_LinRegLA_impl", (DL_FUNC) &_RcppSMC_LinRegLA_impl, 3}, - {"_RcppSMC_LinReg_impl", (DL_FUNC) &_RcppSMC_LinReg_impl, 2}, - {"_RcppSMC_nonLinPMMH_impl", (DL_FUNC) &_RcppSMC_nonLinPMMH_impl, 3}, - {"_RcppSMC_pfLineartBS_impl", (DL_FUNC) &_RcppSMC_pfLineartBS_impl, 4}, - {"_RcppSMC_pfNonlinBS_impl", (DL_FUNC) &_RcppSMC_pfNonlinBS_impl, 2}, - {NULL, NULL, 0} -}; - -RcppExport void R_init_RcppSMC(DllInfo *dll) { - R_registerRoutines(dll, NULL, CallEntries, NULL, NULL); - R_useDynamicSymbols(dll, FALSE); -} From 3e9d4271898378824cde3c6b6090586106844548 Mon Sep 17 00:00:00 2001 From: Leah South Date: Sun, 6 Aug 2017 23:47:33 +1000 Subject: [PATCH 2/2] Regenerating RcppExports with Latest Version of Rcpp --- R/RcppExports.R | 12 ++++++------ src/RcppExports.cpp | 27 +++++++++++++++++++++------ 2 files changed, 27 insertions(+), 12 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index a5cc2aa..e9da87c 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -2,26 +2,26 @@ # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 blockpfGaussianOpt_impl <- function(data, part, lag) { - .Call('RcppSMC_blockpfGaussianOpt_impl', PACKAGE = 'RcppSMC', data, part, lag) + .Call(`_RcppSMC_blockpfGaussianOpt_impl`, data, part, lag) } LinReg_impl <- function(Data, lNumber) { - .Call('RcppSMC_LinReg_impl', PACKAGE = 'RcppSMC', Data, lNumber) + .Call(`_RcppSMC_LinReg_impl`, Data, lNumber) } LinRegLA_impl <- function(Data, intemps, lNumber) { - .Call('RcppSMC_LinRegLA_impl', PACKAGE = 'RcppSMC', Data, intemps, lNumber) + .Call(`_RcppSMC_LinRegLA_impl`, Data, intemps, lNumber) } nonLinPMMH_impl <- function(data, lNumber, lMCMCits) { - .Call('RcppSMC_nonLinPMMH_impl', PACKAGE = 'RcppSMC', data, lNumber, lMCMCits) + .Call(`_RcppSMC_nonLinPMMH_impl`, data, lNumber, lMCMCits) } pfLineartBS_impl <- function(data, part, usef, fun) { - .Call('RcppSMC_pfLineartBS_impl', PACKAGE = 'RcppSMC', data, part, usef, fun) + .Call(`_RcppSMC_pfLineartBS_impl`, data, part, usef, fun) } pfNonlinBS_impl <- function(data, part) { - .Call('RcppSMC_pfNonlinBS_impl', PACKAGE = 'RcppSMC', data, part) + .Call(`_RcppSMC_pfNonlinBS_impl`, data, part) } diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 6821e70..5f6c425 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -8,7 +8,7 @@ using namespace Rcpp; // blockpfGaussianOpt_impl Rcpp::List blockpfGaussianOpt_impl(arma::vec data, long part, long lag); -RcppExport SEXP RcppSMC_blockpfGaussianOpt_impl(SEXP dataSEXP, SEXP partSEXP, SEXP lagSEXP) { +RcppExport SEXP _RcppSMC_blockpfGaussianOpt_impl(SEXP dataSEXP, SEXP partSEXP, SEXP lagSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -21,7 +21,7 @@ END_RCPP } // LinReg_impl Rcpp::List LinReg_impl(arma::mat Data, unsigned long lNumber); -RcppExport SEXP RcppSMC_LinReg_impl(SEXP DataSEXP, SEXP lNumberSEXP) { +RcppExport SEXP _RcppSMC_LinReg_impl(SEXP DataSEXP, SEXP lNumberSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -33,7 +33,7 @@ END_RCPP } // LinRegLA_impl Rcpp::List LinRegLA_impl(arma::mat Data, arma::vec intemps, unsigned long lNumber); -RcppExport SEXP RcppSMC_LinRegLA_impl(SEXP DataSEXP, SEXP intempsSEXP, SEXP lNumberSEXP) { +RcppExport SEXP _RcppSMC_LinRegLA_impl(SEXP DataSEXP, SEXP intempsSEXP, SEXP lNumberSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -46,7 +46,7 @@ END_RCPP } // nonLinPMMH_impl Rcpp::DataFrame nonLinPMMH_impl(arma::vec data, unsigned long lNumber, unsigned long lMCMCits); -RcppExport SEXP RcppSMC_nonLinPMMH_impl(SEXP dataSEXP, SEXP lNumberSEXP, SEXP lMCMCitsSEXP) { +RcppExport SEXP _RcppSMC_nonLinPMMH_impl(SEXP dataSEXP, SEXP lNumberSEXP, SEXP lMCMCitsSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -59,7 +59,7 @@ END_RCPP } // pfLineartBS_impl Rcpp::DataFrame pfLineartBS_impl(arma::mat data, unsigned long part, bool usef, Rcpp::Function fun); -RcppExport SEXP RcppSMC_pfLineartBS_impl(SEXP dataSEXP, SEXP partSEXP, SEXP usefSEXP, SEXP funSEXP) { +RcppExport SEXP _RcppSMC_pfLineartBS_impl(SEXP dataSEXP, SEXP partSEXP, SEXP usefSEXP, SEXP funSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -73,7 +73,7 @@ END_RCPP } // pfNonlinBS_impl Rcpp::List pfNonlinBS_impl(arma::vec data, long part); -RcppExport SEXP RcppSMC_pfNonlinBS_impl(SEXP dataSEXP, SEXP partSEXP) { +RcppExport SEXP _RcppSMC_pfNonlinBS_impl(SEXP dataSEXP, SEXP partSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -83,3 +83,18 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } + +static const R_CallMethodDef CallEntries[] = { + {"_RcppSMC_blockpfGaussianOpt_impl", (DL_FUNC) &_RcppSMC_blockpfGaussianOpt_impl, 3}, + {"_RcppSMC_LinReg_impl", (DL_FUNC) &_RcppSMC_LinReg_impl, 2}, + {"_RcppSMC_LinRegLA_impl", (DL_FUNC) &_RcppSMC_LinRegLA_impl, 3}, + {"_RcppSMC_nonLinPMMH_impl", (DL_FUNC) &_RcppSMC_nonLinPMMH_impl, 3}, + {"_RcppSMC_pfLineartBS_impl", (DL_FUNC) &_RcppSMC_pfLineartBS_impl, 4}, + {"_RcppSMC_pfNonlinBS_impl", (DL_FUNC) &_RcppSMC_pfNonlinBS_impl, 2}, + {NULL, NULL, 0} +}; + +RcppExport void R_init_RcppSMC(DllInfo *dll) { + R_registerRoutines(dll, NULL, CallEntries, NULL, NULL); + R_useDynamicSymbols(dll, FALSE); +}