Public Attributes | Static Public Attributes | Protected Attributes | List of all members
o2scl::mcmc_para_base< func_t, measure_t, data_t, vec_t > Class Template Reference

A generic MCMC simulation class. More...

#include <mcmc_para.h>

Detailed Description

template<class func_t, class measure_t, class data_t, class vec_t = ubvector>
class o2scl::mcmc_para_base< func_t, measure_t, data_t, vec_t >

This class performs a Markov chain Monte Carlo simulation of a user-specified function using OpenMP and/or MPI. Either the Metropolis-Hastings algorithm with a user-specified proposal distribution or the affine-invariant sampling method of Goodman and Weare can be used.

By default, the Metropolis-Hastings algorithm is executed with a simple walk, with steps in each dimension of size $ (\mathrm{high} - \mathrm{low})/\mathrm{step\_fac} $ with the denominator specified in step_fac.

The function type is a template type, func_t, which should be of the form

int f(size_t num_of_parameters, const vec_t &parameters,
double &log_pdf, data_t &dat)

which computes log_pdf, the natural logarithm of the function value, for any point in parameter space (any point between low and high ).

If the function being simulated returns mcmc_skip then the point is automatically rejected. After each acceptance or rejection, a user-specified "measurement" function (of type measure_t ) is called, which can be used to store the results. In order to stop the simulation, either this function or the probability distribution being simulated should return the value mcmc_done .

A generic proposal distribution can be specified in set_proposal(). To go back to the default random walk method, one can call the function unset_proposal().

If aff_inv is set to true and the number of walkers, n_walk is set to a number larger than 1, then affine-invariant sampling is used. For affine-invariant sampling, the variable step_fac represents the value of $ a $, the limits of the distribution for $ z $.

In order to store data at each point, the user can store this data in any object of type data_t . If affine-invariant sampling is used, then each chain has it's own data object. The class keeps twice as many copies of these data object as would otherwise be required, in order to avoid copying of data objects in the case that the steps are accepted or rejected.

Verbose output: If verbose is 0, no output is generated (the default). If verbose is 1, then output to cout occurs only if the settings are somehow misconfigured and the class attempts to recover from them, for example if not enough functions are specified for the requested number of OpenMP threads, or if more than one thread was requested but O2SCL_OPENMP was not defined, or if a negative value for step_fac was requested. When verbose is 1, a couple messages are written to scr_out : a summary of the number of walkers, chains, and threads at the beginning of the MCMC simulation, a message indicating why the MCMC simulation stopped, a message when the warm up iterations are completed, a message every time files are written to disk, and a message at the end counting the number of acceptances and rejections. If verbose is 2, then the file prefix is output to cout during initialization.

Note
This class is experimental.
Currently, this class requires that the data_t has a good copy constructor.
Idea for Future:
The copy constructor for the data_t type is used when the user doesn't specify enough initial points for the corresponding number of threads and walkers. This requirement for a copy constructor could be removed by allowing threads and walkers to share the data_t for the initial point as needed.

Definition at line 134 of file mcmc_para.h.

Public Member Functions

Basic usage
virtual int mcmc (size_t n_params, vec_t &low, vec_t &high, std::vector< func_t > &func, std::vector< measure_t > &meas)
 Perform a MCMC simulation. More...
 
virtual int mcmc (size_t n_params, vec_t &low, vec_t &high, func_t &func, measure_t &meas)
 Perform a MCMC simulation with a thread-safe function or with only one OpenMP thread.
 
Proposal distribution
template<class prob_vec_t >
void set_proposal (prob_vec_t &pv)
 Set the proposal distribution. More...
 
template<class prob_vec_t >
void set_proposal_ptrs (prob_vec_t &pv)
 Set pointers to proposal distributions. More...
 
virtual void unset_proposal ()
 Go back to random-walk Metropolis with a uniform distribution.
 

Public Attributes

bool meas_for_initial
 If true, call the measurement function for the initial point.
 
size_t n_threads
 Number of OpenMP threads.
 
std::vector< ubvectorinitial_points
 Initial points in parameter space. More...
 
Output quantities
std::vector< size_t > n_accept
 The number of Metropolis steps which were accepted in each independent chain (summed over all walkers) More...
 
std::vector< size_t > n_reject
 The number of Metropolis steps which were rejected in each independent chain (summed over all walkers) More...
 
Settings
double mpi_start_time
 The MPI starting time (defaults to 0.0) More...
 
size_t max_iters
 If non-zero, the maximum number of MCMC iterations (default 0) More...
 
double max_time
 Time in seconds (default is 0) More...
 
std::string prefix
 Prefix for output filenames (default "mcmc")
 
bool aff_inv
 If true, use affine-invariant Monte Carlo.
 
double step_fac
 Stepsize factor (default 10.0)
 
size_t n_warm_up
 Number of warm up steps (successful steps not iterations) (default 0) More...
 
int user_seed
 If non-zero, use as the seed for the random number generator (default 0) More...
 
int verbose
 Output control (default 0)
 
size_t max_bad_steps
 Maximum number of failed steps when generating initial points with affine-invariant sampling (default 1000)
 
size_t n_walk
 Number of walkers for affine-invariant MC or 1 otherwise (default 1)
 
size_t n_walk_per_thread
 Number of walkers per thread (default 1)
 
bool err_nonconv
 If true, call the error handler if msolve() or msolve_de() does not converge (default true)
 
bool always_accept
 If true, accept all steps.
 
double ai_initial_step
 Initial step fraction for affine-invariance sampling walkers (default 0.1)
 

Static Public Attributes

static const int mcmc_done =-10
 Integer to indicate completion.
 
static const int mcmc_skip =-20
 Integer to indicate rejection.
 

Protected Member Functions

Interface customization
virtual int mcmc_init ()
 Initializations before the MCMC.
 
virtual void mcmc_cleanup ()
 Cleanup after the MCMC.
 
virtual void best_point (vec_t &best, double w_best, data_t &dat)
 Function to run for the best point.
 

Protected Attributes

std::ofstream scr_out
 The screen output file.
 
std::vector< rng_gslrg
 Random number generators.
 
std::vector< o2scl::prob_cond_mdim< vec_t > * > prop_dist
 Pointer to proposal distribution for each thread.
 
bool pd_mode
 If true, then use the user-specified proposal distribution.
 
bool warm_up
 If true, we are in the warm up phase.
 
std::vector< vec_t > current
 Current points in parameter space for each walker and each OpenMP thread. More...
 
std::vector< data_t > data_arr
 Data array. More...
 
std::vector< bool > switch_arr
 Data switch array for each walker and each OpenMP thread. More...
 
std::vector< std::vector< size_t > > ret_value_counts
 Return value counters, one vector independent chain.
 
std::vector< size_t > curr_walker
 Index of the current walker. More...
 
size_t n_chains_per_rank
 Number of fully independent chains in each MPI rank.
 
MPI properties
int mpi_rank
 The MPI processor rank.
 
int mpi_size
 The MPI number of processors.
 

Member Function Documentation

◆ mcmc()

template<class func_t, class measure_t, class data_t, class vec_t = ubvector>
virtual int o2scl::mcmc_para_base< func_t, measure_t, data_t, vec_t >::mcmc ( size_t  n_params,
vec_t &  low,
vec_t &  high,
std::vector< func_t > &  func,
std::vector< measure_t > &  meas 
)
inlinevirtual

Perform a MCMC simulation over n_params parameters starting at initial point init, limiting the parameters to be between low and high, using func as the objective function and calling the measurement function meas at each MC point.

Definition at line 423 of file mcmc_para.h.

◆ set_proposal()

template<class func_t, class measure_t, class data_t, class vec_t = ubvector>
template<class prob_vec_t >
void o2scl::mcmc_para_base< func_t, measure_t, data_t, vec_t >::set_proposal ( prob_vec_t &  pv)
inline
Note
This function automatically sets aff_inv to false and n_walk to 1.

Definition at line 1709 of file mcmc_para.h.

◆ set_proposal_ptrs()

template<class func_t, class measure_t, class data_t, class vec_t = ubvector>
template<class prob_vec_t >
void o2scl::mcmc_para_base< func_t, measure_t, data_t, vec_t >::set_proposal_ptrs ( prob_vec_t &  pv)
inline
Note
This function automatically sets aff_inv to false and n_walk to 1.

Definition at line 1725 of file mcmc_para.h.

Member Data Documentation

◆ curr_walker

template<class func_t, class measure_t, class data_t, class vec_t = ubvector>
std::vector<size_t> o2scl::mcmc_para_base< func_t, measure_t, data_t, vec_t >::curr_walker
protected

This quantity has to be a vector because different threads may have different values for the current walker during the initialization phase for the affine sampling algorithm.

Definition at line 243 of file mcmc_para.h.

◆ current

template<class func_t, class measure_t, class data_t, class vec_t = ubvector>
std::vector<vec_t> o2scl::mcmc_para_base< func_t, measure_t, data_t, vec_t >::current
protected

This is an array of size n_threads times n_walk initial guesses, indexed by thread_index*n_walk+walker_index .

Definition at line 168 of file mcmc_para.h.

◆ data_arr

template<class func_t, class measure_t, class data_t, class vec_t = ubvector>
std::vector<data_t> o2scl::mcmc_para_base< func_t, measure_t, data_t, vec_t >::data_arr
protected

This is an array of size 2 times n_threads times n_walk . The two copies of data objects are indexed by i_copy*n_walk*n_threads+thread_index*n_walk+walker_index

Definition at line 176 of file mcmc_para.h.

◆ initial_points

template<class func_t, class measure_t, class data_t, class vec_t = ubvector>
std::vector<ubvector> o2scl::mcmc_para_base< func_t, measure_t, data_t, vec_t >::initial_points

To fully specify all of the initial points, this should be a vector of size n_walk times n_threads .

Definition at line 412 of file mcmc_para.h.

◆ max_iters

template<class func_t, class measure_t, class data_t, class vec_t = ubvector>
size_t o2scl::mcmc_para_base< func_t, measure_t, data_t, vec_t >::max_iters

If both max_iters and max_time are nonzero, the MCMC will stop when either the number of iterations exceeds max_iters or the time exceeds max_time, whichever occurs first.

Definition at line 297 of file mcmc_para.h.

◆ max_time

template<class func_t, class measure_t, class data_t, class vec_t = ubvector>
double o2scl::mcmc_para_base< func_t, measure_t, data_t, vec_t >::max_time

If both max_iters and max_time are nonzero, the MCMC will stop when either the number of iterations exceeds max_iters or the time exceeds max_time, whichever occurs first.

Definition at line 306 of file mcmc_para.h.

◆ mpi_start_time

template<class func_t, class measure_t, class data_t, class vec_t = ubvector>
double o2scl::mcmc_para_base< func_t, measure_t, data_t, vec_t >::mpi_start_time

This can be set by the user before mcmc() is called, so that the time required for initializations before the MCMC starts can be counted.

Definition at line 287 of file mcmc_para.h.

◆ n_accept

template<class func_t, class measure_t, class data_t, class vec_t = ubvector>
std::vector<size_t> o2scl::mcmc_para_base< func_t, measure_t, data_t, vec_t >::n_accept

This vector has a size equal to n_chains_per_rank .

Definition at line 269 of file mcmc_para.h.

◆ n_reject

template<class func_t, class measure_t, class data_t, class vec_t = ubvector>
std::vector<size_t> o2scl::mcmc_para_base< func_t, measure_t, data_t, vec_t >::n_reject

This vector has a size equal to n_chains_per_rank .

Definition at line 276 of file mcmc_para.h.

◆ n_warm_up

template<class func_t, class measure_t, class data_t, class vec_t = ubvector>
size_t o2scl::mcmc_para_base< func_t, measure_t, data_t, vec_t >::n_warm_up
Note
Not to be confused with warm_up, which is a protected boolean local variable in some functions which indicates whether we're in warm up mode or not.

Definition at line 325 of file mcmc_para.h.

◆ switch_arr

template<class func_t, class measure_t, class data_t, class vec_t = ubvector>
std::vector<bool> o2scl::mcmc_para_base< func_t, measure_t, data_t, vec_t >::switch_arr
protected

This is an array of size n_threads times n_walk initial guesses, indexed by thread_index*n_walk+walker_index .

Definition at line 183 of file mcmc_para.h.

◆ user_seed

template<class func_t, class measure_t, class data_t, class vec_t = ubvector>
int o2scl::mcmc_para_base< func_t, measure_t, data_t, vec_t >::user_seed

The random number generator is modified so that each thread and each rank has a different set of random numbers.

Definition at line 333 of file mcmc_para.h.


The documentation for this class was generated from the following file:

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).