33 #include <boost/numeric/ublas/vector.hpp> 35 #include <o2scl/uniform_grid.h> 36 #include <o2scl/table3d.h> 37 #include <o2scl/hdf_file.h> 38 #include <o2scl/exception.h> 39 #include <o2scl/prob_dens_func.h> 40 #include <o2scl/cholesky.h> 41 #include <o2scl/vector.h> 42 #include <o2scl/multi_funct.h> 102 template<
class func_t,
class measure_t,
147 virtual void best_point(vec_t &best,
double w_best, data_t &dat) {
250 unsigned long int seed=time(0);
268 virtual int mcmc(
size_t nparams, vec_t &init,
269 vec_t &low, vec_t &high, func_t &func,
272 bool valid_guess=
true;
273 for(
size_t k=0;k<nparams;k++) {
274 if (init[k]<low[k] || init[k]>high[k]) valid_guess=
false;
277 O2SCL_ERR2(
"Initial guess outside of user-specified ",
278 "lower and upper bounds in mcmc_base::mcmc().",
283 if (aff_inv && n_walk<=1) n_walk=2;
299 if (n_warm_up==0) warm_up=
false;
302 current.resize(n_walk);
303 std::vector<double> w_current(n_walk);
305 for(
size_t i=0;i<
n_walk;i++) {
306 current[i].resize(nparams);
311 data_arr.resize(2*n_walk);
312 switch_arr.resize(n_walk);
313 for(
size_t i=0;i<switch_arr.size();i++) switch_arr[i]=
false;
316 vec_t next(nparams), best(nparams);
317 double w_next=0.0, w_best=0.0, q_prop=0.0;
320 for(
size_t i=0;i<nparams;i++) {
321 current[0][i]=init[i];
333 O2SCL_ERR(
"Function mcmc_init() failed in mcmc_base::mcmc().",
340 std::cout <<
"mcmc: Affine-invariant step, n_parameters=" 341 << nparams <<
" n_walk=" << n_walk << std::endl;
342 }
else if (pd_mode==
true) {
343 std::cout <<
"mcmc: With proposal distribution, n_parameters=" 344 << nparams << std::endl;
346 std::cout <<
"mcmc: Random-walk with uniform dist., n_parameters=" 347 << nparams << std::endl;
360 for(
size_t ipar=0;ipar<nparams;ipar++) {
361 init[ipar]=current[0][ipar];
365 for(curr_walker=0;curr_walker<
n_walk;curr_walker++) {
371 if (curr_walker<n_init_points) {
374 int iret=func(nparams,current[curr_walker],
375 w_current[curr_walker],data_arr[curr_walker]);
378 std::cout.precision(4);
379 std::cout <<
"mcmc: ";
380 std::cout.width((
int)(1.0+log10((
double)(n_walk-1))));
381 std::cout << curr_walker <<
" " << w_current[
curr_walker]
382 <<
" " << iret <<
" (initial; ai)" << std::endl;
383 std::cout.precision(6);
390 if (iret>=0 && iret<((
int)ret_value_counts.size())) {
391 ret_value_counts[iret]++;
393 meas_ret=meas(current[curr_walker],w_current[curr_walker],
394 curr_walker,
true,data_arr[curr_walker]);
402 for(
size_t ipar=0;ipar<nparams;ipar++) {
403 if (init[ipar]<low[ipar] || init[ipar]>high[ipar]) {
407 ") in mcmc_base::mcmc().").c_str(),
412 (high[ipar]-low[ipar])*ai_initial_step;
413 }
while (current[curr_walker][ipar]>high[ipar] ||
414 current[curr_walker][ipar]<low[ipar]);
418 int iret=func(nparams,current[curr_walker],
419 w_current[curr_walker],data_arr[curr_walker]);
422 std::cout.precision(4);
423 std::cout <<
"mcmc: ";
424 std::cout.width((
int)(1.0+log10((
double)(n_walk-1))));
425 std::cout << curr_walker <<
" " << w_current[
curr_walker]
426 <<
" " << iret <<
" (initial; ai)" << std::endl;
427 std::cout.precision(6);
438 if (iret>=0 && iret<((
int)ret_value_counts.size())) {
439 ret_value_counts[iret]++;
441 meas_ret=meas(current[curr_walker],w_current[curr_walker],
442 curr_walker,
true,data_arr[curr_walker]);
444 }
else if (init_iters>max_bad_steps) {
453 std::cout <<
"Press a key and type enter to continue. ";
469 int iret=func(nparams,current[0],w_current[0],data_arr[0]);
472 std::cout.precision(4);
473 std::cout <<
"mcmc: " << w_current[0] <<
" (initial)" << std::endl;
474 std::cout.precision(6);
479 O2SCL_ERR(
"Initial weight vanished in mcmc_base::mcmc().",
485 if (iret>=0 && iret<((
int)ret_value_counts.size())) {
486 ret_value_counts[iret]++;
488 meas_ret=meas(current[0],w_current[0],0,
true,data_arr[0]);
494 if (meas_ret==mcmc_done) {
503 bool main_done=
false;
518 curr_walker=mcmc_iters %
n_walk;
523 ij=((size_t)(rg.
random()*((double)n_walk)));
524 }
while (ij==curr_walker || ij>=n_walk);
529 smove_z=(1.0-2.0*p+2.0*a*p+p*p-2.0*a*p*p+a*a*p*p)/a;
532 for(
size_t i=0;i<nparams;i++) {
533 next[i]=current[ij][i]+smove_z*(current[
curr_walker][i]-
537 }
else if (pd_mode) {
540 (*prop_dist)(current[0],next);
541 q_prop=prop_dist->
log_pdf(current[0],next)-
542 prop_dist->
log_pdf(next,current[0]);
543 if (!std::isfinite(q_prop)) {
544 O2SCL_ERR2(
"Proposal distribution not finite in ",
551 for(
size_t k=0;k<nparams;k++) {
552 next[k]=current[0][k]+(rg.
random()*2.0-1.0)*
553 (high[k]-low[k])/step_fac;
563 for(
size_t k=0;k<nparams;k++) {
564 if (next[k]<low[k] || next[k]>high[k]) iret=
mcmc_skip;
566 if (iret!=mcmc_skip) {
567 if (switch_arr[curr_walker]==
false) {
568 iret=func(nparams,next,w_next,data_arr[curr_walker+n_walk]);
570 iret=func(nparams,next,w_next,data_arr[curr_walker]);
572 if (iret>=0 && iret<((
int)ret_value_counts.size())) {
573 ret_value_counts[iret]++;
576 std::cout <<
"Function returned failure " << iret
578 for(
size_t k=0;k<nparams;k++) {
579 std::cout << next[k] <<
" ";
581 std::cout << std::endl;
583 }
else if (verbose>=2) {
584 std::cout <<
"Parameter(s) out of range: " << std::endl;
585 std::cout.setf(std::ios::showpos);
586 for(
size_t k=0;k<nparams;k++) {
587 std::cout << k <<
" " << low[k] <<
" " << next[k] <<
" " 589 if (next[k]<low[k] || next[k]>high[k]) {
592 std::cout << std::endl;
594 std::cout.unsetf(std::ios::showpos);
601 if (switch_arr[curr_walker]==
false) {
602 best_point(best,w_best,data_arr[curr_walker+n_walk]);
604 best_point(best,w_best,data_arr[curr_walker]);
611 if (always_accept && iret==
success) accept=
true;
617 if (r<pow(smove_z,((
double)nparams)-1.0)*
618 exp(w_next-w_current[curr_walker])) {
622 std::cout.precision(4);
623 std::cout <<
"mcmc: ";
624 std::cout.width((
int)(1.0+log10((
double)(nparams-1))));
625 std::cout << curr_walker <<
" " 627 << smove_z <<
" ratio: " 628 << pow(smove_z,((
double)nparams)-1.0)*
629 exp(w_next-w_current[curr_walker])
630 <<
" accept: " << accept << std::endl;
631 std::cout.precision(6);
633 }
else if (pd_mode) {
634 if (r<exp(w_next-w_current[0]+q_prop)) {
638 std::cout.precision(4);
639 std::cout <<
"mcmc: " << w_current[0]
640 <<
" " << w_next <<
" " << q_prop <<
" ratio: " 641 << exp(w_next-w_current[0]+q_prop)
642 <<
" accept: " << accept << std::endl;
643 std::cout.precision(6);
647 if (r<exp(w_next-w_current[0])) {
651 std::cout.precision(4);
652 std::cout <<
"mcmc: " << w_current[0] <<
" " << w_next
654 << exp(w_next-w_current[0])
655 <<
" accept: " << accept << std::endl;
656 std::cout.precision(6);
669 if (switch_arr[curr_walker]==
false) {
670 meas_ret=meas(next,w_next,curr_walker,
true,
671 data_arr[curr_walker+n_walk]);
673 meas_ret=meas(next,w_next,curr_walker,
true,
674 data_arr[curr_walker]);
690 if (switch_arr[curr_walker]==
false) {
691 meas_ret=meas(current[curr_walker],w_current[curr_walker],
692 curr_walker,
false,data_arr[curr_walker]);
694 meas_ret=meas(current[curr_walker],w_current[curr_walker],
695 curr_walker,
false,data_arr[curr_walker+n_walk]);
701 if (iret==mcmc_done) {
706 if (meas_ret!=mcmc_done && err_nonconv) {
707 O2SCL_ERR((((std::string)
"Measurement function returned ")+
709 " in mcmc_base::mcmc().").c_str(),
716 if (warm_up && mcmc_iters==n_warm_up) {
722 std::cout <<
"Finished warmup." << std::endl;
728 std::cout <<
"Press a key and type enter to continue. ";
804 template<
class func_t,
class fill_t,
class data_t,
class vec_t=ubvector>
806 std::function<int(const vec_t &,double,size_t,bool,data_t &)>,
812 typedef std::function<int(const vec_t &,double,size_t,bool,data_t &)>
825 std::shared_ptr<o2scl::table_units<> >
tab;
834 std::cout <<
"Start mcmc_table::mcmc_init()." << std::endl;
842 tab->new_column(
"mult");
843 tab->new_column(
"log_wgt");
844 for(
size_t i=0;i<col_names.size();i++) {
845 tab->new_column(col_names[i]);
846 if (col_units[i].length()>0) {
847 tab->set_unit(col_names[i],col_units[i]);
851 walker_rows.resize(this->
n_walk);
852 for(
size_t i=0;i<this->
n_walk;i++) {
857 std::cout <<
"mcmc: Table column names and units: " << std::endl;
858 for(
size_t i=0;i<tab->get_ncolumns();i++) {
859 std::cout << tab->get_column_name(i) <<
" " 860 << tab->get_unit(tab->get_column_name(i)) << std::endl;
865 std::cout <<
"End mcmc_table::mcmc_init()." << std::endl;
873 virtual int fill_line(
const vec_t &pars,
double log_weight,
874 std::vector<double> &line, data_t &dat,
879 line.push_back(log_weight);
880 for(
size_t i=0;i<pars.size();i++) {
881 line.push_back(pars[i]);
883 return fill(pars,log_weight,line,dat);
901 std::vector<std::string> units) {
914 virtual int mcmc(
size_t nparams, vec_t &init,
915 vec_t &low, vec_t &high, func_t &func,
919 (std::mem_fn<
int(
const vec_t &,
double,
size_t,
bool,data_t &,fill_t &)>
921 std::placeholders::_2,std::placeholders::_3,std::placeholders::_4,
922 std::placeholders::_5,std::ref(fill));
924 return parent_t::mcmc(nparams,init,low,high,func,meas);
943 virtual int add_line(
const vec_t &pars,
double log_weight,
944 size_t walker_ix,
bool new_meas, data_t &dat,
951 if (new_meas==
true ||
953 walker_rows[this->
curr_walker]>=((
int)(tab->get_nlines()))) {
955 std::vector<double> line;
956 int fret=fill_line(pars,log_weight,line,dat,fill);
965 std::cout <<
"Fill function returned mcmc_done. " 966 <<
"Stopping run." << std::endl;
971 std::cout <<
"Fill function returned " << fret
972 <<
". Stopping run." << std::endl;
978 if (line.size()!=tab->get_ncolumns()) {
979 std::cout <<
"line: " << line.size() <<
" columns: " 980 << tab->get_ncolumns() << std::endl;
981 O2SCL_ERR(
"Table misalignment in mcmc_table::add_line().",
986 std::cout <<
"mcmc: Adding line:" << std::endl;
987 std::vector<std::string> sc_in, sc_out;
988 for(
size_t k=0;k<line.size();k++) {
989 sc_in.push_back(tab->get_column_name(k)+
": "+
993 for(
size_t k=0;k<sc_out.size();k++) {
994 std::cout << sc_out[k] << std::endl;
999 tab->line_of_data(line.size(),line);
1011 double mult_old=tab->get(
"mult",walker_rows[this->
curr_walker]);
1012 tab->set(
"mult",walker_rows[this->curr_walker],mult_old+1.0);
1015 std::cout <<
"mcmc: Updating line:" << std::endl;
1016 std::vector<std::string> sc_in, sc_out;
1017 for(
size_t k=0;k<tab->get_ncolumns();k++) {
1019 (tab->get_column_name(k)+
": "+
1024 for(
size_t k=0;k<sc_out.size();k++) {
1025 std::cout << sc_out[k] << std::endl;
1049 size_t n=tab->get_nlines();
1051 O2SCL_ERR2(
"Cannot reblock. Not enough data in ",
1054 size_t n_block=n/n_blocks;
1055 size_t m=tab->get_ncolumns();
1056 for(
size_t j=0;j<n_blocks;j++) {
1059 for(
size_t i=0;i<m;i++) {
1062 for(
size_t k=j*n_block;k<(j+1)*n_block;k++) {
1063 mult+=(*tab)[
"mult"][k];
1064 for(
size_t i=1;i<m;i++) {
1065 dat[i]+=(*tab)[i][k]*(*tab)[
"mult"][k];
1068 tab->set(
"mult",j,mult);
1069 for(
size_t i=1;i<m;i++) {
1071 tab->set(i,j,dat[i]);
1074 tab->set_nlines(n_blocks);
std::shared_ptr< o2scl::table_units<> > get_table()
Get the output table.
std::shared_ptr< o2scl::table_units<> > tab
Main data table for Markov chain.
virtual void set_seed()
Default method for setting the random seed.
virtual int add_line(const vec_t &pars, double log_weight, size_t walker_ix, bool new_meas, data_t &dat, fill_t &fill)
A measurement function which adds the point to the table.
std::vector< std::string > col_names
Column names.
std::vector< std::string > col_units
Column units.
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
virtual double log_pdf(const vec_t &x, const vec_t &x2) const =0
The log of the normalized density.
std::vector< vec_t > current
Current points in parameter space.
o2scl::prob_cond_mdim< vec_t > * prop_dist
Proposal distribution.
double step_fac
Stepsize factor (default 10.0)
void set_table(std::shared_ptr< o2scl::table_units<> > &t)
Set the output table.
std::vector< bool > switch_arr
Data switch array.
size_t n_walk
Number of walkers for affine-invariant MC or 1 otherwise (default 1)
virtual void unset_proposal()
Go back to random-walk Metropolis with a uniform distribution.
void set_seed(unsigned long int s)
Set the seed.
size_t n_reject
The number of Metropolis steps which were rejected.
int verbose
Output control (default 0)
virtual int mcmc(size_t nparams, vec_t &init, vec_t &low, vec_t &high, func_t &func, fill_t &fill)
Perform an MCMC simulation.
virtual void set_proposal(o2scl::prob_cond_mdim< vec_t > &p)
Set the proposal distribution.
A generic MCMC simulation class.
invalid argument supplied by user
bool always_accept
If true, accept all steps.
virtual void mcmc_cleanup()
Cleanup after the MCMC.
double random()
Return a random number in .
A multi-dimensional conditional probability density function.
size_t n_init_points
Number of initial points specified by the user;.
void reblock(size_t n_blocks)
Reaverage the data into blocks of a fixed size in order to avoid autocorrelations.
static const int mcmc_done
Integer to indicate completion.
size_t n_accept
The number of Metropolis steps which were accepted.
size_t n_warm_up
Number of warm up steps (successful steps not iterations) (default 0)
A generic MCMC simulation class writing data to a o2scl::table_units object.
bool warm_up
If true, we are in the warm up phase.
std::vector< data_t > data_arr
Data array.
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
virtual void set_names_units(std::vector< std::string > names, std::vector< std::string > units)
Set the table names and units.
int user_seed
If non-zero, use as the seed for the random number generator (default 0)
std::string dtos(double x, int prec=6, bool auto_prec=false)
Convert a double to a string.
bool aff_inv
If true, use affine-invariant Monte Carlo.
size_t max_bad_steps
Maximum number of failed steps when generating initial points with affine-invariant sampling (default...
size_t curr_walker
Index of the current walker.
void screenify(size_t nin, const string_arr_t &in_cols, std::vector< std::string > &out_cols, size_t max_size=80)
Reformat the columns for output of width size.
mcmc_base< func_t, internal_measure_t, data_t, vec_t > parent_t
Type of parent class.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
virtual void best_point(vec_t &best, double w_best, data_t &dat)
Function to run for the best point.
std::vector< size_t > ret_value_counts
Return value counters.
Data table table class with units.
static const int mcmc_skip
Integer to indicate rejection.
Random number generator (GSL)
virtual int mcmc_init()
Initializations before the MCMC.
bool pd_mode
If true, then use the user-specified proposal distribution.
bool err_nonconv
If true, call the error handler if msolve() or msolve_de() does not converge (default true) ...
std::function< int(const vec_t &, double, size_t, bool, data_t &)> internal_measure_t
Measurement functor type for the parent.
double ai_initial_step
Initial step fraction for affine-invariance sampling walkers.
std::string szttos(size_t x)
Convert a size_t to a string.
std::vector< int > walker_rows
Record the last row in the table which corresponds to each walker.
virtual int fill_line(const vec_t &pars, double log_weight, std::vector< double > &line, data_t &dat, fill_t &fill)
Fill line with data for insertion into the table.
virtual int mcmc(size_t nparams, vec_t &init, vec_t &low, vec_t &high, func_t &func, measure_t &meas)
Perform an MCMC simulation.
rng_gsl rg
Random number generator.
virtual int mcmc_init()
MCMC initialization function.