mcmc.hpp

Go to the documentation of this file.

#ifndef SRC_INVERSION_MCMC_HPP_
#define SRC_INVERSION_MCMC_HPP_
#include <cstddef>
 
#ifdef MPI_PARALLEL
#include <mpi.h>
#endif
 
void _choldc(double **a, int n, double *p);
 
int _acor(double *mean, double *sigma, double *tau, double *X, int L);
 
struct mcmc_opts {
  double a, b, c, d;
  int p, q, r;
 
  int print;
  char logfile[80];
 
#ifdef MPI_PARALLEL
  MPI_Comm mpi_comm;
#endif
};
 
struct mcmc_recs {
  int ndim,    
      nvalue,  
      nwalker, 
      nstep;   
  double ***par, 
      ***val,    
      **lnp;     
  double *opt_par, 
      *opt_val;    
  int **newstate; 
  int cur,        
      accept,     
      reset;      
  double opt_lnp; 
};
 
void mcmc_alloc(mcmc_recs *recs, int nstep, int nwalker, int ndim, int nvalue);
 
void mcmc_free(mcmc_recs *recs);
 
void mcmc_statistics(double *mean, double *sigma, double *tau, mcmc_recs *recs);
 
void mcmc_save_fits(char const *fname, mcmc_opts *opts, mcmc_recs *recs,
                    int include_last = false);
 
void mcmc_load_fits(char const *fname, mcmc_opts *opts, mcmc_recs *recs,
                    int alloc = true);
 
void mcmc_append_recs(mcmc_recs *dst, mcmc_recs *src);
 
// report MCMC sampler result
void mcmc_report(mcmc_opts *opts, mcmc_recs *recs, char const *mode);
 
typedef double (*ObjectiveFunction_t)(double *, double *, int, int, void *);
 
void mcmc_init(ObjectiveFunction_t lnprob, double **par, mcmc_opts *opts,
               mcmc_recs *recs, void *obj);
 
void mcmc_advance(ObjectiveFunction_t lnprob, mcmc_opts *opts, mcmc_recs *recs,
                  void *obj);
 
double mcmc_stretch_move(double *newp, double **oldp, int iwalker, int nwalker,
                         int ndim, mcmc_opts *opts);
 
void mcmc_walk_move(double *newp, double **oldp, int k, int nwalker, int np,
                    mcmc_opts *opts);
 
#endif  // SRC_INVERSION_MCMC_HPP_