Sampler

Defines the overall sampler object.

This file defines the smc::sampler class which is used to implement entire particle systems.

Defines

__SMC_SAMPLER_HH

namespace ResampleType

Specifiers for various resampling algorithms:

Enums

enum Enum

Values:

enumerator MULTINOMIAL

enumerator RESIDUAL

enumerator STRATIFIED

enumerator SYSTEMATIC

namespace PathSamplingType

Specifiers for various path sampling methods:

Enums

enum Enum

Values:

enumerator TRAPEZOID2

enumerator TRAPEZOID1

enumerator RECTANGLE

namespace HistoryType

Storage types for the history of the particle system.

Enums

enum Enum

Values:

enumerator NONE

enumerator RAM

enumerator AL

namespace smc

The Sequential Monte Carlo namespace.

The classes and functions within this namespace are intended to be used for producing implemenetations of SMC samplers and related simulation techniques.

class nullParams: #include <sampler.h>

An empty class for use when additional algorithm parameters are not required.

template<class Space, class Params = nullParams> class sampler

#include <sampler.h>

A template class for an interacting particle system suitable for SMC sampling.

Public Functions

sampler(long lSize, HistoryType::Enum htHistoryMode)

Create an particle system containing lSize uninitialised particles with the specified mode.

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.

Parameters

lSize – The number of particles present in the ensemble (at time 0 if this is a variable quantity)
htHM – The history mode to use: set this to HistoryType::RAM to store the whole history of the system and HistoryType::NONE to avoid doing so.

Template Parameters

Space – The class used to represent a point in the sample space.
Params – (optional) The class used for any additional parameters.

sampler(long lSize, HistoryType::Enum htHistoryMode, moveset<Space, Params>*)

Create an particle system containing lSize uninitialised particles with the specified mode.

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.

Parameters

lSize – The number of particles present in the ensemble (at time 0 if this is a variable quantity)
htHM – The history mode to use: set this to HistoryType::RAM to store the whole history of the system and HistoryType::NONE to avoid doing so.
pNewMoves – Pointer to a moveset.

Template Parameters

Space – The class used to represent a point in the sample space.
Params – (optional) The class used for any additional parameters.

~sampler(): Dispose of a sampler.

sampler(const sampler<Space, Params> &sFrom): Copy constructor.

sampler<Space, Params> &operator=(const sampler<Space, Params> &sFrom): Assignment overloading.

double GetESS(void) const: Calculates and Returns the Effective Sample Size.

inline int GetAccepted(void) const: Returns the number of accepted proposals from the most recent MCMC iteration.

inline int GetResampled(void) const: Returns a flag for whether the ensemble was resampled during the most recent iteration.

inline int GetMcmcRepeats(void) const: Returns the number of MCMC repeats used in the most recent iteration.

inline const std::vector<historyelement<Space>> &GetHistory(void) const: Returns the History of the particle system.

inline long GetNumber(void) const: Returns the number of particles within the system.

inline long GetHistoryLength(void) const: Returns the number of evolution times stored in the history.

inline population<Space> GetHistoryPopulation(long n) const: Returns the current particle set stored in the history.

inline const population<Space> &GetHistoryPopulationRefs(long n) const: Returns a reference to the particle set stored in the history.

inline historyflags GetHistoryFlags(long n) const: Returns the history flags.

inline double GetESS(long n) const: Returns the Effective Sample Size of the specified particle generation.

std::vector<population<Space>> GetAPop(void) const: Returns a vector of population objects re-ordered to match ancestral lines of the particle system.

arma::Col<unsigned int> GetALineInd(long n) const: Returns indices for ancestral line of particle ‘n’.

std::vector<Space> GetALineSpace(long n) const: Returns ancestral line of particle ‘n’.

inline arma::Col<unsigned int> GetuRSIndices(void) const: Returns current iteration resampled indices.

inline int GetuRSIndex(long n) const: Returns current iteration resampled index number ‘n’.

inline int GetHistorymcmcRepeats(long n): Returns the history number of MCMC iterations performed during this iteration.

inline const Params &GetAlgParams(void) const: Returns the additional algorithm parameters.

inline const Space &GetParticleValueN(long n) const: Return the value of particle n.

inline double GetParticleLogWeightN(long n) const: Return the logarithmic unnormalized weight of particle n.

inline double GetParticleWeightN(long n) const: Return the unnormalized weight of particle n.

inline arma::vec GetParticleWeight(void) const: Return the unnormalized particle weights.

inline long GetTime(void) const: Returns the current evolution time of the system.

inline double GetLogNCPath(void) const: Returns the current estimate of the log normalising constant ratio over the entire path.

inline double GetLogNCStep(void) const: Returns the current estimate of the log normalising constant ratio over the last step.

inline double GetNCPath(void) const: Returns the current estimate of the normalising constant ratio over the entire path.

inline double GetNCStep(void) const: Returns the current estimate of the normalising constant ratio over the last step.

void Initialise(void)

Initialise the sampler and its constituent particles.

At present this function resets the system evolution time to 0 and calls the moveset initialisor to assign each particle in the ensemble.

Note that the initialisation function must be specified before calling this function.

double Integrate(double (*pIntegrand)(const Space&, void*), void *pAuxiliary)

Integrate the supplied function with respect to the current particle set.

This function returns the result of integrating the supplied function under the empirical measure associated with the particle set at the present time. The final argument of the integrand function is a pointer which will be supplied with pAuxiliary to allow for arbitrary additional information to be passed to the function being integrated.

Parameters

pIntegrand – The function to integrate with respect to the particle set
pAuxiliary – A pointer to any auxiliary data which should be passed to the function

double IntegratePathSampling(PathSamplingType::Enum, double (*pIntegrand)(long, const Space&, void*), double (*pWidth)(long, void*), void *pAuxiliary)

Integrate the supplied function over the path using the supplied width function and integration method.

This function is intended to be used to estimate integrals of the sort which must be evaluated to determine the normalising constant of a distribution obtained using a sequence of potential functions proportional to densities with respect to the initial distribution to define a sequence of distributions leading up to the terminal, interesting distribution.

In this context, the particle set at each time is used to make an estimate of the path sampling integrand, and numerical integration is then performed to obtain an estimate of the path sampling integral which is the natural logarithm of the ratio of normalising densities.

The integrand is integrated at every time point in the population history. The results of this integration are taken to be point-evaluations of the path sampling integrand which are spaced on a grid of intervals given by the width function. The path sampling integral is then calculated by performing a suitable numerical integration and the results of this integration is returned.

pAuxiliary is passed to both of the user specified functions to allow the user to pass additional data to either or both of these functions in a convenient manner. It is safe to use NULL if no such data is used by either function.

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.

Parameters

PStype – The numerical integration method to use
pIntegrand – The function to integrated. The first argument is evolution time, the second the particle value at which the function is to be evaluated and the final argument is always pAuxiliary.
pWidth – The function which returns the width of the path sampling grid at the specified evolution time. The final argument is always pAuxiliary
pAuxiliary – A pointer to auxiliary data to pass to both of the above functions

Template Parameters

Space – The class used to represent a point in the sample space.
Params – (optional) The class used for any additional parameters.

inline double IntegratePathSampling(double (*pIntegrand)(long, const Space&, void*), double (*pWidth)(long, void*), void *pAuxiliary): Integrate the supplied function over the path using the supplied width function and the default integration method (the corrected trapezoid rule).

void Iterate(void)

Perform one iteration of the simulation algorithm.

The iterate function:

moves the current particle set
checks the effective sample size and resamples if necessary
performs a mcmc step if required
appends the current particle set to the history if desired
increments the current evolution time

void IterateBack(void): Cancel one iteration of the simulation algorithm.

double IterateEss(void): Perform one iteration of the simulation algorithm and return the resulting ess.

void IterateUntil(long lTerminate): Perform iterations until the specified evolution time is reached.

void MoveParticles(void): Move the particle set by proposing and applying an appropriate move to each particle.

void Resample(ResampleType::Enum lMode): Resample the particle set using the specified resampling scheme.

inline void SetMoveSet(moveset<Space, Params> *pNewMoves): Sets the entire moveset to the one which is supplied.

inline void SetMoveSet(moveset<Space, Params> &NewMoves): Sets the entire moveset to the one which is supplied - backwards compatibility.

void SetResampleParams(ResampleType::Enum rtMode, double dThreshold)

Set Resampling Parameters.

This function configures the resampling parameters, allowing the specification of both the resampling mode and the threshold at which resampling is used.

The rtMode parameter should be set to one of the following:

ResampleType::MULTINOMIAL to use multinomial resampling
ResampleType::RESIDUAL to use residual resampling
ResampleType::STRATIFIED to use stratified resampling
ResampleType::SYSTEMATIC to use systematic resampling

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.

Parameters

rtMode – The resampling mode to be used.
dThreshold – The threshold at which resampling is deemed necesary.

inline void SetAlgParam(Params parameters): Set additional algorithm parameters.

inline void SetAdaptMethods(adaptMethods<Space, Params> *adaptMethod): Set the methods to adapt the additional algorithm parameters.

inline void SetMcmcRepeats(int reps): Sets the number of MCMC repeats.

std::ostream &StreamParticle(std::ostream &os, long n) const

Dump a specified particle to the specified output stream in a human readable form.

Produce a human-readable display of the current nth particle value and log weight.

Parameters

os – The output stream to which the display should be made.
n – The index of the particle of interest

std::ostream &StreamParticles(std::ostream &os) const

Dump the entire particle set to the specified output stream in a human readable form.

Produce a human-readable display of the current particle values and log weights.

Parameters: os – The output stream to which the display should be made.

void OstreamMCMCRecordToStream(std::ostream &os) const

Output a vector indicating the number of accepted MCMC moves at each time instance.

This function records the MCMC acceptance history to the specified output stream as a list of the number of moves accepted at each time instant.

Parameters: os – The output stream to send the data to.

void OstreamResamplingRecordToStream(std::ostream &os) const

Output a 0-1 value vector indicating the times at which resampling occured to an output stream.

This function records the resampling history to the specified output stream as a 0-1 valued list which takes the value 1 for those time instances when resampling occured and 0 otherwise.

Parameters: os – The output stream to send the data to.

Protected Functions

void _copy(const sampler<Space, Params> &sFrom): Helper function for copy constructor and assignment overloading.

inline double CalcLogNC(void) const: Returns the crude normalising constant ratio estimate implied by the weights.

Protected Attributes

long N: Number of particles in the system.

long T: The current evolution time of the system.

ResampleType::Enum rtResampleMode: The resampling mode which is to be employed.

double dResampleThreshold: The effective sample size at which resampling should be used.

arma::vec dRSWeights: Structure used internally for resampling.

arma::Col<int> uRSCount: Structure used internally for resampling.

arma::Col<unsigned int> uRSIndices: Structure used internally for resampling.

population<Space> pPopulation: The particles within the system.

moveset<Space, Params> *pMoves: The set of moves available.

bool movesetBelong: A flag to track whether the moveset object needs to be included in this destructor.

Params algParams: The additional algorithm parameters.

adaptMethods<Space, Params> *pAdapt: An object for adapting additional algorithm parameters.

bool adaptBelong: A flag to track whether the adaptation object needs to be included in this destructor.

int nAccepted: The number of MCMC moves which have been accepted during this iteration.

int nResampled: A flag which tracks whether the ensemble was resampled during this iteration.

int nRepeats: The number of MCMC repeats to be performed. The default is 1 if an MCMC step is supplied.

double acceptProb: 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 dlogNCPath: An estimate of the log normalising constant ratio over the entire path.

double dlogNCIt: An estimate of the log normalising constant ratio over the last step.

HistoryType::Enum htHistoryMode: A mode flag which indicates whether historical information is stored.

std::vector<historyelement<Space>> History: The historical process associated with the particle system.

namespace std

The standard namespace.

The classes provided within the standard libraries reside within this namespace and the TDSMC class library adds a number of additional operator overloads to some of the standard classes to allow them to deal with our structures.

Functions

template<class Space, class Params> std::ostream &operator<<(std::ostream &os, smc::sampler<Space, Params> &s)

Produce a human-readable display of the state of an smc::sampler class using the stream operator.

Parameters

os – The output stream to which the display should be made.
s – The sampler which is to be displayed.