27 #ifndef O2SCL_MCMC_PARA_OLD_H 28 #define O2SCL_MCMC_PARA_OLD_H 40 #include <boost/numeric/ublas/vector.hpp> 42 #include <o2scl/uniform_grid.h> 43 #include <o2scl/table3d.h> 44 #include <o2scl/hdf_file.h> 45 #include <o2scl/exception.h> 46 #include <o2scl/prob_dens_func.h> 47 #include <o2scl/vector.h> 48 #include <o2scl/multi_funct.h> 49 #include <o2scl/vec_stats.h> 50 #include <o2scl/cli.h> 133 template<
class func_t,
class measure_t,
151 std::vector<rng_gsl>
rg;
196 O2SCL_ERR2(
"Using a proposal distribution with affine-invariant ",
197 "sampling not implemented in mcmc_para_old::mcmc_init().",
202 std::cout <<
"Prefix is: " <<
prefix << std::endl;
209 scr_out.setf(std::ios::scientific);
221 scr_out <<
"mcmc (" << it <<
"," << mpi_rank
223 <<
" reject=" <<
n_reject[it] << std::endl;
232 virtual void best_point(vec_t &best,
double w_best, data_t &dat) {
400 MPI_Comm_rank(MPI_COMM_WORLD,&this->mpi_rank);
401 MPI_Comm_size(MPI_COMM_WORLD,&this->mpi_size);
407 meas_for_initial=
true;
429 virtual int mcmc(
size_t n_params, vec_t &low, vec_t &high,
430 std::vector<func_t> &func, std::vector<measure_t> &meas) {
436 if (func.size()==0 || meas.size()==0) {
437 O2SCL_ERR2(
"Size of 'func' or 'meas' array is zero in ",
442 std::cout <<
"mcmc_para_old::mcmc(): Not enough functions for " 443 << n_threads <<
" threads. Setting n_threads to " 444 << func.size() <<
"." << std::endl;
446 n_threads=func.size();
450 std::cout <<
"mcmc_para_old::mcmc(): Not enough measurement objects for " 451 << n_threads <<
" threads. Setting n_threads to " 452 << meas.size() <<
"." << std::endl;
454 n_threads=meas.size();
458 if (mpi_start_time==0.0) {
460 mpi_start_time=MPI_Wtime();
462 mpi_start_time=time(0);
466 if (initial_points.size()==0) {
469 initial_points.resize(1);
470 initial_points[0].resize(n_params);
471 for(
size_t k=0;k<n_params;k++) {
472 initial_points[0][k]=(low[k]+high[k])/2.0;
479 for(
size_t iip=0;iip<initial_points.size();iip++) {
480 for(
size_t ipar=0;ipar<n_params;ipar++) {
481 if (initial_points[iip][ipar]<low[ipar] ||
482 initial_points[iip][ipar]>high[ipar]) {
488 ") in mcmc_base::mcmc().").c_str(),
495 for(
size_t i=0;i<initial_points.size();i++) {
496 for(
size_t k=0;k<initial_points[i].size();k++) {
497 if (!std::isfinite(initial_points[i][k])) {
502 for(
size_t j=i+1;j<initial_points.size();j++) {
504 for(
size_t k=0;k<initial_points[i].size();k++) {
505 if (initial_points[i][k]!=initial_points[j][k]) {
510 std::cerr.setf(std::ios::scientific);
511 std::cerr << i <<
" ";
513 std::cerr << j <<
" ";
525 omp_set_num_threads(n_threads);
528 std::cout <<
"mcmc_para_old::mcmc(): " 529 << n_threads <<
" threads were requested but the " 530 <<
"-DO2SCL_OPENMP flag was not used during " 531 <<
"compilation. Setting n_threads to 1." 538 std::vector<int> func_ret(n_threads), meas_ret(n_threads);
543 std::cout <<
"mcmc_para_old::mcmc(): Requested only 1 walker, " 544 <<
"forcing 2 walkers." << std::endl;
547 #ifdef O2SCL_NEVER_DEFINED 553 if (n_walk_per_thread>n_walk) {
554 O2SCL_ERR2(
"More walkers per thread than total walkers ",
557 n_chains_per_rank=n_walk_per_thread*n_threads/
n_walk;
558 if (n_chains_per_walk*n_walk!=n_walk_per_thread*n_threads) {
559 std::cout <<
"mcmc_para_old::mcmc(): Could not evenly " 560 <<
"organize threads and walkers." << std::endl;
561 std::cout <<
"n_threads: " << n_threads << std::endl;
562 std::cout <<
"n_walk: " << n_walk << std::endl;
563 std::cout <<
"n_walk_per_thread: " 564 << n_walk_per_thread << << std::endl;
565 std::cout <<
"n_chains_per_rank: " 566 << n_chains_per_rank << << std::endl;
577 std::cout <<
"mcmc_para_old::mcmc(): Requested negative or zero " 578 <<
"step_fac with aff_inv=true.\nSetting to 2.0." 582 std::cout <<
"mcmc_para_old::mcmc(): Requested negative or zero " 583 <<
"step_fac. Setting to 10.0." << std::endl;
589 rg.resize(n_threads);
590 unsigned long int seed=time(0);
591 if (this->user_seed!=0) {
595 seed*=(mpi_rank*n_threads+it+1);
596 rg[it].set_seed(seed);
601 n_accept.resize(n_chains_per_rank);
602 n_reject.resize(n_chains_per_rank);
611 if (n_warm_up==0) warm_up=
false;
617 current.resize(ssize);
618 std::vector<double> w_current(ssize);
619 for(
size_t i=0;i<ssize;i++) {
620 current[i].resize(n_params);
625 ret_value_counts.resize(n_chains_per_rank);
626 curr_walker.resize(n_threads);
629 data_arr.resize(2*ssize);
630 switch_arr.resize(ssize);
631 for(
size_t i=0;i<switch_arr.size();i++) switch_arr[i]=
false;
634 std::vector<vec_t> next(n_threads);
636 next[it].resize(n_params);
638 std::vector<double> w_next(n_threads);
641 vec_t best(n_params);
646 std::vector<bool> mcmc_done_flag(n_threads);
648 mcmc_done_flag[it]=
false;
652 std::vector<double> q_prop(n_threads);
659 O2SCL_ERR(
"Function mcmc_init() failed in mcmc_base::mcmc().",
670 scr_out <<
"mcmc: Affine-invariant step, n_params=" 671 << n_params <<
", n_walk=" << n_walk
672 <<
", n_chains_per_rank=" << n_chains_per_rank
673 <<
",\n\tn_walk_per_thread=" << n_walk_per_thread
674 <<
", n_threads=" << n_threads <<
", rank=" 675 << mpi_rank <<
", n_ranks=" 676 << mpi_size << std::endl;
677 }
else if (pd_mode==
true) {
678 scr_out <<
"mcmc: With proposal distribution, n_params=" 679 << n_params <<
", n_threads=" << n_threads <<
", rank=" 680 << mpi_rank <<
", n_ranks=" 681 << mpi_size << std::endl;
683 scr_out <<
"mcmc: Random-walk w/uniform dist., n_params=" 684 << n_params <<
", n_threads=" << n_threads <<
", rank=" 685 << mpi_rank <<
", n_ranks=" 686 << mpi_size << std::endl;
688 scr_out <<
"Set start time to: " << mpi_start_time << std::endl;
697 #pragma omp parallel default(shared) 706 for(curr_walker[it]=0;curr_walker[it]<n_walk_per_thread &&
707 mcmc_done_flag[it]==
false;curr_walker[it]++) {
710 size_t sindex=n_walk*it+curr_walker[it];
713 size_t ip_index=sindex % initial_points.size();
721 if (sindex<initial_points.size()) {
724 for(
size_t ipar=0;ipar<n_params;ipar++) {
725 current[sindex][ipar]=initial_points[ip_index][ipar];
729 func_ret[it]=func[it](n_params,current[sindex],
730 w_current[sindex],data_arr[sindex]);
732 if (func_ret[it]==mcmc_done) {
733 mcmc_done_flag[it]=
true;
738 if (func_ret[it]>=0 &&
739 func_ret[it]<((
int)ret_value_counts[it].size())) {
740 ret_value_counts[it][func_ret[it]]++;
742 if (meas_for_initial) {
743 meas_ret[it]=meas[it](current[sindex],w_current[sindex],
744 curr_walker[it],func_ret[it],
745 true,data_arr[sindex]);
749 if (meas_ret[it]==mcmc_done) {
750 mcmc_done_flag[it]=
true;
759 while (!done && !mcmc_done_flag[it]) {
762 for(
size_t ipar=0;ipar<n_params;ipar++) {
764 current[sindex][ipar]=
765 initial_points[ip_index][ipar]+
766 (rg[it].random()*2.0-1.0)*
767 (high[ipar]-low[ipar])*ai_initial_step;
768 }
while (current[sindex][ipar]>high[ipar] ||
769 current[sindex][ipar]<low[ipar]);
773 func_ret[it]=func[it](n_params,current[sindex],
774 w_current[sindex],data_arr[sindex]);
781 if (func_ret[it]==mcmc_done) {
782 mcmc_done_flag[it]=
true;
787 if (func_ret[it]>=0 &&
788 func_ret[it]<((
int)ret_value_counts[it].size())) {
789 ret_value_counts[it][func_ret[it]]++;
791 if (meas_ret[it]!=mcmc_done) {
792 if (meas_for_initial) {
793 meas_ret[it]=meas[it](current[sindex],w_current[sindex],
794 curr_walker[it],func_ret[it],
true,
800 mcmc_done_flag[it]=
true;
803 }
else if (init_iters>max_bad_steps) {
804 std::string err=((std::string)
"In loop with thread ")+
807 " in mcmc_para_old_base::mcmc().";
818 bool stop_early=
false;
820 if (mcmc_done_flag[it]==
true) {
822 scr_out <<
"mcmc (" << it <<
"," << mpi_rank
823 <<
"): Returned mcmc_done " 824 <<
"(initial; ai)." << std::endl;
838 for(curr_walker[it]=0;curr_walker[it]<
n_walk;
840 size_t sindex=n_walk*it+curr_walker[it];
841 if (w_current[sindex]>w_current[0]) {
843 w_best=w_current[sindex];
847 best=current[best_index];
853 for(curr_walker[it]=0;curr_walker[it]<
n_walk;curr_walker[it]++) {
854 size_t sindex=n_walk*it+curr_walker[it];
855 scr_out.precision(4);
856 scr_out <<
"mcmc (" << it <<
"," << mpi_rank <<
"): i_walk: ";
857 scr_out.width((
int)(1.0+log10((
double)(n_walk-1))));
858 scr_out << curr_walker[it] <<
" log_wgt: " 859 << w_current[sindex] <<
" (initial; ai)" << std::endl;
860 scr_out.precision(6);
873 #pragma omp parallel default(shared) 887 size_t ip_size=initial_points.size();
888 for(
size_t ipar=0;ipar<n_params;ipar++) {
889 current[it][ipar]=initial_points[it % ip_size][ipar];
895 func_ret[it]=func[it](n_params,current[it],w_current[it],
899 func_ret[it]=func_ret[it % ip_size];
900 w_current[it]=w_current[it % ip_size];
903 for(
size_t j=0;j<data_arr.size();j++) {
904 data_arr[it]=data_arr[it % ip_size];
915 if (func_ret[it]==mcmc_done) {
917 scr_out <<
"mcmc (" << it
918 <<
"): Initial point returned mcmc_done." 926 O2SCL_ERR((((std::string)
"Initial weight from thread ")+
928 " vanished in mcmc_para_old_base::mcmc().").c_str(),
939 #pragma omp parallel default(shared) 947 size_t ip_size=initial_points.size();
951 func_ret[it]=func_ret[it % ip_size];
952 current[it]=current[it % ip_size];
953 w_current[it]=w_current[it % ip_size];
954 data_arr[it]=data_arr[it % ip_size];
958 if (func_ret[it]>=0 &&
959 func_ret[it]<((
int)ret_value_counts[it].size())) {
960 ret_value_counts[it][func_ret[it]]++;
962 if (meas_for_initial) {
964 meas_ret[it]=meas[it](current[it],w_current[it],0,
965 func_ret[it],
true,data_arr[it]);
969 if (meas_ret[it]==mcmc_done) {
970 mcmc_done_flag[it]=
true;
981 bool stop_early=
false;
983 if (mcmc_done_flag[it]==
true) {
985 scr_out <<
"mcmc (" << it <<
"," << mpi_rank
986 <<
"): Returned mcmc_done " 987 <<
"(initial)." << std::endl;
1002 if (w_current[it]<w_best) {
1004 w_best=w_current[it];
1010 scr_out.precision(4);
1011 scr_out <<
"mcmc: ";
1013 scr_out << w_current[it] <<
" ";
1015 scr_out <<
" (initial)" << std::endl;
1016 scr_out.precision(6);
1023 meas_for_initial=
true;
1030 std::cout <<
"Press a key and type enter to continue. ";
1039 bool main_done=
false;
1040 size_t mcmc_iters=0;
1042 while (!main_done) {
1045 std::vector<double> smove_z(n_threads);
1051 curr_walker[it]=mcmc_iters %
n_walk;
1056 #pragma omp parallel default(shared) 1071 ij=((size_t)(rg[it].random()*((double)n_walk)));
1072 }
while (ij==curr_walker[it] || ij>=n_walk);
1075 double p=rg[it].random();
1077 smove_z[it]=(1.0-2.0*p+2.0*a*p+p*p-2.0*a*p*p+a*a*p*p)/a;
1080 for(
size_t i=0;i<n_params;i++) {
1081 next[it][i]=current[n_walk*it+ij][i]+
1082 smove_z[it]*(current[n_walk*it+curr_walker[it]][i]-
1083 current[n_walk*it+ij][i]);
1086 }
else if (pd_mode) {
1089 q_prop[it]=prop_dist[it]->log_metrop_hast(current[it],next[it]);
1091 if (!std::isfinite(q_prop[it])) {
1092 std::cout <<
"Here: " << q_prop[it] << std::endl;
1093 O2SCL_ERR2(
"Proposal distribution not finite in ",
1100 for(
size_t k=0;k<n_params;k++) {
1101 next[it][k]=current[it][k]+(rg[it].random()*2.0-1.0)*
1102 (high[k]-low[k])/step_fac;
1113 for(
size_t k=0;k<n_params;k++) {
1114 if (next[it][k]<low[k] || next[it][k]>high[k]) {
1117 if (next[it][k]<low[k]) {
1118 std::cout <<
"mcmc (" << it <<
"," 1119 << mpi_rank <<
"): Parameter with index " << k
1120 <<
" and value " << next[it][k]
1121 <<
" smaller than limit " << low[k] << std::endl;
1122 scr_out <<
"mcmc (" << it <<
"," 1123 << mpi_rank <<
"): Parameter with index " << k
1124 <<
" and value " << next[it][k]
1125 <<
" smaller than limit " << low[k] << std::endl;
1127 std::cout <<
"mcmc (" << it <<
"," << mpi_rank
1128 <<
"): Parameter with index " << k
1129 <<
" and value " << next[it][k]
1130 <<
" larger than limit " << high[k] << std::endl;
1131 scr_out <<
"mcmc (" << it <<
"," << mpi_rank
1132 <<
"): Parameter with index " << k
1133 <<
" and value " << next[it][k]
1134 <<
" larger than limit " << high[k] << std::endl;
1142 if (func_ret[it]!=mcmc_skip) {
1143 if (switch_arr[n_walk*it+curr_walker[it]]==
false) {
1144 func_ret[it]=func[it](n_params,next[it],w_next[it],
1145 data_arr[it*n_walk+curr_walker[it]+
1148 func_ret[it]=func[it](n_params,next[it],w_next[it],
1149 data_arr[it*n_walk+curr_walker[it]]);
1151 if (func_ret[it]==mcmc_done) {
1152 mcmc_done_flag[it]=
true;
1154 if (func_ret[it]>=0 &&
1155 func_ret[it]<((
int)ret_value_counts[it].size())) {
1156 ret_value_counts[it][func_ret[it]]++;
1174 scr_out <<
"it: " << it <<
" q_prop[it]: " 1175 << q_prop[it] << std::endl;
1177 if (func_ret[it]==mcmc_done) {
1178 scr_out <<
"mcmc (" << it <<
"," << mpi_rank
1179 <<
"): Returned mcmc_done." 1181 }
else if (func_ret[it]==mcmc_skip && verbose>=3) {
1182 scr_out <<
"mcmc (" << it
1183 <<
"): Parameter(s) out of range: " << std::endl;
1184 scr_out.setf(std::ios::showpos);
1185 for(
size_t k=0;k<n_params;k++) {
1186 scr_out << k <<
" " << low[k] <<
" " 1187 << next[it][k] <<
" " << high[k];
1188 if (next[it][k]<low[k] || next[it][k]>high[k]) {
1191 scr_out << std::endl;
1193 scr_out.unsetf(std::ios::showpos);
1195 func_ret[it]!=mcmc_skip) {
1197 scr_out <<
"mcmc (" << it <<
"," << mpi_rank
1198 <<
"): Function returned failure " 1199 << func_ret[it] <<
" at point ";
1200 for(
size_t k=0;k<n_params;k++) {
1201 scr_out << next[it][k] <<
" ";
1203 scr_out << std::endl;
1214 #pragma omp parallel default(shared) 1223 size_t sindex=n_walk*it+curr_walker[it];
1229 if (always_accept && func_ret[it]==
success) accept=
true;
1232 double r=rg[it].random();
1235 double ai_ratio=pow(smove_z[it],((
double)n_params)-1.0)*
1236 exp(w_next[it]-w_current[sindex]);
1240 }
else if (pd_mode) {
1241 if (r<exp(w_next[it]-w_current[sindex]+q_prop[it])) {
1246 if (r<exp(w_next[it]-w_current[sindex])) {
1260 if (switch_arr[sindex]==
false) {
1261 meas_ret[it]=meas[it](next[it],w_next[it],
1262 curr_walker[it],func_ret[it],
true,
1263 data_arr[sindex+n_threads*
n_walk]);
1265 meas_ret[it]=meas[it](next[it],w_next[it],
1266 curr_walker[it],func_ret[it],
true,
1272 current[sindex]=next[it];
1273 w_current[sindex]=w_next[it];
1274 switch_arr[sindex]=!(switch_arr[sindex]);
1283 if (switch_arr[sindex]==
false) {
1284 meas_ret[it]=meas[it](next[it],w_next[it],
1285 curr_walker[it],func_ret[it],
false,
1286 data_arr[sindex+n_threads*
n_walk]);
1288 meas_ret[it]=meas[it](next[it],w_next[it],
1289 curr_walker[it],func_ret[it],
false,
1308 size_t sindex=n_walk*it+curr_walker[it];
1309 scr_out.precision(4);
1310 scr_out <<
"mcmc (" << it <<
"," << mpi_rank <<
"): iter: ";
1311 scr_out.width((
int)(1.0+log10((
double)(n_params-1))));
1312 scr_out << mcmc_iters <<
" i_walk: " 1313 << curr_walker[it] <<
" log_wgt: " 1314 << w_current[sindex] << std::endl;
1315 scr_out.precision(6);
1324 if (switch_arr[n_walk*it+curr_walker[it]]==
false) {
1325 best_point(best,w_best,data_arr[curr_walker[it]+n_walk*it+
1328 best_point(best,w_best,data_arr[curr_walker[it]+n_walk*it]);
1336 if (meas_ret[it]==mcmc_done || func_ret[it]==mcmc_done) {
1341 O2SCL_ERR((((std::string)
"Measurement function returned ")+
1343 " in mcmc_para_old_base::mcmc().").c_str(),
1352 if (main_done==
false) {
1356 if (warm_up && mcmc_iters==n_warm_up) {
1364 scr_out <<
"mcmc: Finished warmup." << std::endl;
1371 std::cout <<
"Press a key and type enter to continue. ";
1377 if (main_done==
false && warm_up==
false && max_iters>0 &&
1378 mcmc_iters==max_iters) {
1380 scr_out <<
"mcmc: Stopping because number of iterations (" 1381 << mcmc_iters <<
") equal to max_iters (" << max_iters
1382 <<
")." << std::endl;
1387 if (main_done==
false) {
1394 if (max_time>0.0 && elapsed>max_time) {
1396 scr_out <<
"mcmc: Stopping because elapsed (" << elapsed
1397 <<
") > max_time (" << max_time <<
")." 1421 virtual int mcmc(
size_t n_params, vec_t &low, vec_t &high,
1422 func_t &func, measure_t &meas) {
1425 omp_set_num_threads(n_threads);
1429 std::vector<func_t> vf(n_threads);
1430 std::vector<measure_t> vm(n_threads);
1435 return mcmc(n_params,low,high,func,meas);
1447 prop_dist.resize(pv.size());
1448 for(
size_t i=0;i<pv.size();i++) {
1449 prop_dist[i]=&pv[i];
1463 prop_dist.resize(pv.size());
1464 for(
size_t i=0;i<pv.size();i++) {
1528 template<
class func_t,
class fill_t,
class data_t,
class vec_t=ubvector>
1530 std::function<int(const vec_t &,double,size_t,int,bool,data_t &)>,
1536 typedef std::function<int(const vec_t &,double,size_t,int,bool,data_t &)>
1552 std::shared_ptr<o2scl::table_units<> >
table;
1571 table->new_column(
"rank");
1572 table->new_column(
"thread");
1573 table->new_column(
"walker");
1574 table->new_column(
"mult");
1575 table->new_column(
"log_wgt");
1576 for(
size_t i=0;i<col_names.size();i++) {
1577 table->new_column(col_names[i]);
1578 if (col_units[i].length()>0) {
1579 table->set_unit(col_names[i],col_units[i]);
1585 walker_accept_rows[i]=-1;
1587 walker_reject_rows.resize(this->
n_walk*this->n_threads);
1589 walker_reject_rows[i]=-1;
1593 std::cout <<
"mcmc: Table column names and units: " << std::endl;
1594 for(
size_t i=0;i<table->get_ncolumns();i++) {
1595 std::cout << table->get_column_name(i) <<
" " 1596 << table->get_unit(table->get_column_name(i)) << std::endl;
1605 if (table->get_ncolumns()!=5+col_names.size()) {
1606 O2SCL_ERR2(
"Table does not have the correct number of columns ",
1609 if (!table->is_column(
"rank") ||
1610 !table->is_column(
"thread") ||
1611 !table->is_column(
"walker") ||
1612 !table->is_column(
"mult") ||
1613 !table->is_column(
"log_wgt")) {
1614 O2SCL_ERR2(
"Table does not have the correct internal columns ",
1618 O2SCL_ERR2(
"Array walker_accept_rows does not have correct size ",
1622 O2SCL_ERR2(
"Array walker_reject_rows does not have correct size ",
1634 last_write_time=MPI_Wtime();
1636 last_write_time=time(0);
1639 return parent_t::mcmc_init();
1645 std::vector<double> &line, data_t &dat,
1646 size_t i_walker, fill_t &fill) {
1649 size_t i_thread=omp_get_thread_num();
1657 line.push_back(i_thread);
1659 line.push_back(i_walker);
1661 line.push_back(1.0);
1662 line.push_back(log_weight);
1663 for(
size_t i=0;i<pars.size();i++) {
1664 line.push_back(pars[i]);
1666 int tempi=fill(pars,log_weight,line,dat);
1743 this->
scr_out <<
"mcmc: Start write_files(). mpi_rank: " 1745 << this->
mpi_size <<
" table_io_chunk: " 1746 << table_io_chunk << std::endl;
1749 std::vector<o2scl::table_units<> > tab_arr;
1750 bool rank_sent=
false;
1753 if (table_io_chunk>1) {
1754 if (this->
mpi_rank%table_io_chunk==0) {
1756 for(
int i=0;i<table_io_chunk-1;i++) {
1760 tab_arr.push_back(t);
1761 o2scl_table_mpi_recv(child,tab_arr[tab_arr.size()-1]);
1767 o2scl_table_mpi_send(*table,parent);
1776 int tag=0, buffer=0;
1777 if (sync_write && this->
mpi_size>1 &&
1779 MPI_Recv(&buffer,1,MPI_INT,this->
mpi_rank-table_io_chunk,
1780 tag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
1788 if (first_write==
false) {
1799 hf.
set_szt(
"n_params",this->n_params);
1807 hf.
seti(
"store_rejects",this->store_rejects);
1808 hf.
seti(
"table_sequence",this->table_sequence);
1824 hf.
seti(
"n_tables",tab_arr.size()+1);
1825 if (rank_sent==
false) {
1826 hdf_output(hf,*table,
"markov_chain_0");
1828 for(
size_t i=0;i<tab_arr.size();i++) {
1829 std::string name=((std::string)
"markov_chain_")+
szttos(i+1);
1830 hdf_output(hf,tab_arr[i],name);
1836 if (sync_write && this->
mpi_size>1 &&
1838 MPI_Send(&buffer,1,MPI_INT,this->
mpi_rank+table_io_chunk,
1839 tag,MPI_COMM_WORLD);
1844 this->
scr_out <<
"mcmc: Done write_files()." << std::endl;
1852 file_update_iters=0;
1853 file_update_time=0.0;
1855 store_rejects=
false;
1856 table_sequence=
true;
1865 std::vector<std::string> units) {
1866 if (names.size()!=units.size()) {
1867 O2SCL_ERR2(
"Size of names and units arrays don't match in ",
1894 int tag=0, buffer=0;
1896 MPI_Recv(&buffer,1,MPI_INT,this->
mpi_rank-1,
1897 tag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
1909 MPI_Send(&buffer,1,MPI_INT,this->
mpi_rank+1,
1910 tag,MPI_COMM_WORLD);
1918 std::cout <<
"Initial points: Finding last " << n_points
1919 <<
" points from file named " 1920 << fname <<
" ." << std::endl;
1925 for(
size_t it=0;it<this->
n_threads;it++) {
1926 for(
size_t iw=0;iw<this->
n_walk;iw++) {
1929 size_t windex=it*this->n_walk+iw;
1932 for(
int row=tip.
get_nlines()-1;row>=0 && found==
false;row--) {
1933 if (tip.
get(
"walker",row)==iw &&
1934 tip.
get(
"thread",row)==it &&
1935 tip.
get(
"mult",row)>0.5) {
1939 std::cout <<
"Function initial_point_file_last():\n\tit: " 1941 <<
" iw: " << iw <<
" row: " 1942 << row <<
" log_wgt: " << tip.
get(
"log_wgt",row)
1947 for(
size_t ip=0;ip<n_param_loc;ip++) {
1953 O2SCL_ERR(
"Function initial_points_file_last() failed.",
1980 int tag=0, buffer=0;
1982 MPI_Recv(&buffer,1,MPI_INT,this->
mpi_rank-1,
1983 tag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
1995 MPI_Send(&buffer,1,MPI_INT,this->
mpi_rank+1,
1996 tag,MPI_COMM_WORLD);
2004 std::cout <<
"Initial points: Finding last " << n_points
2005 <<
" points from file named " 2006 << fname <<
" ." << std::endl;
2012 size_t decrement=nlines/n_points;
2013 if (decrement<1) decrement=1;
2016 for(
size_t k=0;k<n_points;k++) {
2020 std::cout <<
"Function initial_point_file_dist():\n\trow: " 2021 << row <<
" log_wgt: " << tip.
get(
"log_wgt",row)
2026 for(
size_t ip=0;ip<n_param_loc;ip++) {
2046 double thresh=1.0e-6,
2054 int tag=0, buffer=0;
2056 MPI_Recv(&buffer,1,MPI_INT,this->
mpi_rank-1,
2057 tag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
2069 MPI_Send(&buffer,1,MPI_INT,this->
mpi_rank+1,
2070 tag,MPI_COMM_WORLD);
2078 std::cout <<
"Initial points: Finding best " << n_points
2079 <<
" unique points from file named " 2080 << fname <<
" ." << std::endl;
2083 typedef std::map<double,int,std::greater<double> > map_t;
2088 m.insert(std::make_pair(tip.
get(
"log_wgt",k),k));
2098 for(map_t::iterator mit=m.begin();mit!=m.end();mit++) {
2099 map_t::iterator mit2=mit;
2101 if (mit2!=m.end()) {
2102 if (fabs(mit->first-mit2->first)<thresh) {
2104 std::cout <<
"Removing duplicate weights: " 2105 << mit->first <<
" " << mit2->first << std::endl;
2114 }
while (found==
true);
2117 if (m.size()<n_points) {
2118 O2SCL_ERR2(
"Could not find enough points in file in ",
2119 "mcmc_para_old::initial_points_file_best().",
2125 map_t::iterator mit=m.begin();
2126 for(
size_t k=0;k<n_points;k++) {
2127 int row=mit->second;
2129 std::cout <<
"Initial point " << k <<
" at row " 2130 << row <<
" has log_weight= " 2131 << tip.
get(
"log_wgt",row) << std::endl;
2134 for(
size_t ip=0;ip<n_param_loc;ip++) {
2150 virtual int mcmc(
size_t n_params_local,
2151 vec_t &low, vec_t &high, std::vector<func_t> &func,
2152 std::vector<fill_t> &fill) {
2154 n_params=n_params_local;
2170 std::vector<internal_measure_t> meas(this->
n_threads);
2171 for(
size_t it=0;it<this->
n_threads;it++) {
2173 (std::mem_fn<
int(
const vec_t &,
double,
size_t,
int,
bool,
2174 data_t &,
size_t, fill_t &)>
2176 std::placeholders::_2,std::placeholders::_3,
2177 std::placeholders::_4,std::placeholders::_5,
2178 std::placeholders::_6,it,std::ref(fill[it]));
2181 return parent_t::mcmc(n_params,low,high,func,meas);
2206 chain_sizes.resize(ntot);
2208 for(
size_t it=0;it<this->
n_threads;it++) {
2209 for(
size_t iw=0;iw<this->
n_walk;iw++) {
2210 size_t ix=it*this->n_walk+iw;
2213 for(
size_t j=istart;j<table->get_nlines();j+=ntot) {
2214 if (table->get(
"mult",j)>0.5) chain_sizes[ix]++;
2234 std::string name=
"") {
2242 if (!table->is_column(
"rank") ||
2243 !table->is_column(
"thread") ||
2244 !table->is_column(
"walker") ||
2245 !table->is_column(
"mult") ||
2246 !table->is_column(
"log_wgt")) {
2247 O2SCL_ERR2(
"Table does not have the correct internal columns ",
2248 "in mcmc_para_old_table::read_prev_results().",
2260 walker_accept_rows.resize(ntot);
2261 walker_reject_rows.resize(ntot);
2263 for(
size_t j=0;j<ntot;j++) {
2264 walker_accept_rows[j]=-1;
2265 walker_reject_rows[j]=-1;
2268 for(
size_t j=0;j<table->get_nlines();j++) {
2270 size_t i_thread=((size_t)(table->get(
"thread",j)+1.0e-12));
2271 size_t i_walker=((size_t)(table->get(
"walker",j)+1.0e-12));
2274 size_t windex=i_thread*this->n_walk+i_walker;
2276 if (table->get(
"mult",j)>0.5) {
2277 walker_accept_rows[windex]=j;
2278 }
else if (table->get(
"mult",j)<-0.5) {
2279 walker_reject_rows[windex]=j;
2287 for(
size_t j=0;j<ntot;j++) {
2288 if (walker_accept_rows[j]<0) found=
false;
2295 for(
size_t j=0;j<ntot;j++) {
2297 for(
size_t k=0;k<n_param_loc;k++) {
2298 this->
initial_points[j][k]=table->get(k+5,walker_accept_rows[j]);
2302 std::cout <<
"Previous table was read, but initial points not set." 2307 std::cout <<
"mcmc_para_old_table::read_prev_results():" << std::endl;
2308 std::cout <<
" index walker_accept_rows walker_reject_rows" 2310 for(
size_t j=0;j<ntot;j++) {
2313 std::cout << j <<
" ";
2315 std::cout << walker_accept_rows[j] <<
" ";
2317 std::cout << walker_reject_rows[j] << std::endl;
2330 virtual int add_line(
const vec_t &pars,
double log_weight,
2331 size_t walker_ix,
int func_ret,
2332 bool mcmc_accept, data_t &dat,
2333 size_t i_thread, fill_t &fill) {
2336 size_t windex=i_thread*this->
n_walk+walker_ix;
2345 if ((mcmc_accept || store_rejects) && walker_accept_rows[windex]<0) {
2348 if (table_sequence) {
2349 if (walker_accept_rows[windex]>walker_reject_rows[windex]) {
2350 next_row=walker_accept_rows[windex]+ntot;
2352 next_row=walker_reject_rows[windex]+ntot;
2355 if (walker_accept_rows[windex]>walker_reject_rows[windex]) {
2356 next_row=walker_accept_rows[windex]+1;
2358 next_row=walker_reject_rows[windex]+1;
2367 #pragma omp critical (o2scl_mcmc_para_old_table_add_line) 2374 while (next_row<((
int)table->get_nlines()) &&
2375 fabs(table->get(
"mult",next_row))>0.1) {
2384 if (next_row>=((
int)table->get_nlines())) {
2385 size_t istart=table->get_nlines();
2387 table->set_nlines(table->get_nlines()+ntot);
2390 for(
size_t i=0;i<this->
n_walk;i++) {
2391 table->set(
"rank",istart+j*this->n_walk+i,this->
mpi_rank);
2392 table->set(
"thread",istart+j*this->n_walk+i,j);
2393 table->set(
"walker",istart+j*this->n_walk+i,i);
2394 table->set(
"mult",istart+j*this->n_walk+i,0.0);
2395 table->set(
"log_wgt",istart+j*this->n_walk+i,0.0);
2401 if (func_ret==0 && (mcmc_accept || store_rejects)) {
2403 if (next_row>=((
int)(table->get_nlines()))) {
2407 std::vector<double> line;
2408 int fret=fill_line(pars,log_weight,line,dat,walker_ix,fill);
2412 if (store_rejects && mcmc_accept==
false) {
2428 if (line.size()!=table->get_ncolumns()) {
2429 std::cout <<
"line: " << line.size() <<
" columns: " 2430 << table->get_ncolumns() << std::endl;
2431 for(
size_t k=0;k<table->get_ncolumns() || k<line.size();k++) {
2432 std::cout << k <<
". ";
2433 if (k<table->get_ncolumns()) {
2434 std::cout << table->get_column_name(k) <<
" ";
2435 std::cout << table->get_unit(table->get_column_name(k)) <<
" ";
2437 if (k<line.size()) std::cout << line[k] <<
" ";
2438 std::cout << std::endl;
2440 O2SCL_ERR(
"Table misalignment in mcmc_para_old_table::add_line().",
2445 table->set_row(((
size_t)next_row),line);
2449 this->
scr_out <<
"mcmc: Setting data at row " << next_row
2451 for(
size_t k=0;k<line.size();k++) {
2453 this->
scr_out << table->get_column_name(k) <<
" ";
2454 this->
scr_out << table->get_unit(table->get_column_name(k));
2455 this->
scr_out <<
" " << line[k] << std::endl;
2466 if (walker_accept_rows[windex]<0 ||
2467 walker_accept_rows[windex]>=((
int)table->get_nlines())) {
2468 O2SCL_ERR2(
"Invalid row for incrementing multiplier in ",
2471 double mult_old=table->get(
"mult",walker_accept_rows[windex]);
2476 table->set(
"mult",walker_accept_rows[windex],mult_old+1.0);
2478 this->
scr_out <<
"mcmc: Updating mult of row " 2479 << walker_accept_rows[windex]
2480 <<
" from " << mult_old <<
" to " 2481 << mult_old+1.0 << std::endl;
2492 walker_accept_rows[windex]=next_row;
2493 }
else if (store_rejects && func_ret==0) {
2494 walker_reject_rows[windex]=next_row;
2508 if (file_update_iters>0) {
2509 size_t total_iters=0;
2513 if (total_iters>=last_write_iters+file_update_iters) {
2515 this->
scr_out <<
"mcmc: Writing to file. total_iters: " 2516 << total_iters <<
" file_update_iters: " 2517 << file_update_iters <<
" last_write_iters: " 2518 << last_write_iters << std::endl;
2521 last_write_iters=total_iters;
2525 if (updated==
false && file_update_time>0.0) {
2527 double elapsed=MPI_Wtime()-last_write_time;
2529 double elapsed=time(0)-last_write_time;
2531 if (elapsed>file_update_time) {
2533 this->
scr_out <<
"mcmc: Writing to file. elapsed: " 2534 << elapsed <<
" file_update_time: " 2535 << file_update_time <<
" last_write_time: " 2536 << last_write_time << std::endl;
2540 last_write_time=MPI_Wtime();
2542 last_write_time=time(0);
2558 for(i=table->get_nlines()-1;i>=0 && done==
false;i--) {
2560 if (fabs(table->get(
"mult",i))<0.1) {
2564 if (i+2<((
int)table->get_nlines())) {
2565 table->set_nlines(i+2);
2570 return parent_t::mcmc_cleanup();
2576 std::vector<size_t> chain_sizes;
2577 get_chain_sizes(chain_sizes);
2578 size_t min_size=chain_sizes[0];
2579 for(
size_t i=1;i<chain_sizes.size();i++) {
2580 if (chain_sizes[i]<min_size) min_size=chain_sizes[i];
2582 size_t N_max=min_size/2;
2583 ac_coeffs.resize(ncols,N_max-1);
2584 for(
size_t i=0;i<ncols;i++) {
2585 for(
size_t ell=1;ell<N_max;ell++) {
2586 ac_coeffs(i,ell-1)=0.0;
2591 size_t cstart=table->lookup_column(
"log_wgt")+1;
2592 for(
size_t i=0;i<ncols;i++) {
2594 for(
size_t k=0;k<this->
n_walk;k++) {
2595 size_t tindex=j*this->n_walk+k;
2596 for(
size_t ell=1;ell<N_max;ell++) {
2597 const double &x=(*table)[cstart+i][table_row];
2598 double mean=o2scl::vector_mean<const double *>
2599 (chain_sizes[tindex]+1,&x);
2601 <
const double *>(chain_sizes[tindex]+1,&x,ell,mean);
2603 table_row+=chain_sizes[tindex]+1;
2606 for(
size_t ell=1;ell<N_max;ell++) {
2607 ac_coeffs(i,ell-1)/=((double)n_tot);
2616 ubvector &ac_lengths) {
2617 size_t N_max=ac_coeffs_cols.size2();
2618 ac_lengths.resize(ncols);
2619 for(
size_t icol=0;icol<ncols;icol++) {
2620 std::vector<double> tau(N_max);
2621 for(
size_t i=5;i<N_max;i++) {
2623 for(
size_t j=0;j<i;j++) {
2624 sum+=ac_coeffs_cols(icol,j);
2627 std::cout << tau[i] <<
" " << ((double)i)/5.0 << std::endl;
2629 std::cout << std::endl;
2639 std::shared_ptr<o2scl::table_units<> > table2=
2643 for(
size_t j=0;j<this->
n_walk;j++) {
2644 std::string func=std::string(
"abs(walker-")+
o2scl::szttos(j)+
2646 table->copy_rows(func,*table2);
2668 for(
size_t it=0;it<this->
n_threads;it++) {
2670 size_t n=table->get_nlines();
2672 O2SCL_ERR2(
"Cannot reblock. Not enough data in ",
2675 size_t n_block=n/n_blocks;
2676 size_t m=table->get_ncolumns();
2677 for(
size_t j=0;j<n_blocks;j++) {
2680 for(
size_t i=0;i<m;i++) {
2683 for(
size_t k=j*n_block;k<(j+1)*n_block;k++) {
2684 mult+=(*table)[
"mult"][k];
2685 for(
size_t i=1;i<m;i++) {
2686 dat[i]+=(*table)[i][k]*(*table)[
"mult"][k];
2689 table->set(
"mult",j,mult);
2690 for(
size_t i=1;i<m;i++) {
2692 table->set(i,j,dat[i]);
2695 table->set_nlines(n_blocks);
2711 template<
class func_t,
class fill_t,
class data_t,
class vec_t=ubvector>
2744 hf.
sets_vec(
"cl_args",this->cl_args);
2792 p_file_update_iters.
s=&this->file_update_iters;
2793 p_file_update_iters.
help=((std::string)
"Number of MCMC successes ")+
2794 "between file updates (default 0 for no file updates).";
2795 cl.
par_list.insert(std::make_pair(
"file_update_iters",
2796 &p_file_update_iters));
2798 p_file_update_time.
d=&this->file_update_time;
2799 p_file_update_time.
help=((std::string)
"Time ")+
2800 "between file updates (default 0.0 for no file updates).";
2801 cl.
par_list.insert(std::make_pair(
"file_update_time",
2802 &p_file_update_time));
2813 p_max_time.
help=((std::string)
"Maximum run time in seconds ")+
2814 "(default 86400 sec or 1 day).";
2815 cl.
par_list.insert(std::make_pair(
"max_time",&p_max_time));
2818 p_max_iters.
help=((std::string)
"If non-zero, limit the number of ")+
2819 "iterations to be less than the specified number (default zero).";
2820 cl.
par_list.insert(std::make_pair(
"max_iters",&p_max_iters));
2823 p_prefix.
help=
"Output file prefix (default 'mcmc\').";
2824 cl.
par_list.insert(std::make_pair(
"prefix",&p_prefix));
2834 p_step_fac.
help=((std::string)
"MCMC step factor. The step size for ")+
2835 "each variable is taken as the difference between the high and low "+
2836 "limits divided by this factor (default 10.0). This factor can "+
2837 "be increased if the acceptance rate is too small, but care must "+
2838 "be taken, e.g. if the conditional probability is multimodal. If "+
2839 "this step size is smaller than 1.0, it is reset to 1.0 .";
2840 cl.
par_list.insert(std::make_pair(
"step_fac",&p_step_fac));
2843 p_n_warm_up.
help=((std::string)
"Minimum number of warm up iterations ")+
2845 cl.
par_list.insert(std::make_pair(
"n_warm_up",&p_n_warm_up));
2848 p_verbose.
help=((std::string)
"MCMC verbosity parameter ")+
2850 cl.
par_list.insert(std::make_pair(
"mcmc_verbose",&p_verbose));
2853 p_max_bad_steps.
help=((std::string)
"Maximum number of bad steps ")+
2855 cl.
par_list.insert(std::make_pair(
"max_bad_steps",&p_max_bad_steps));
2858 p_n_walk.
help=((std::string)
"Number of walkers ")+
2860 cl.
par_list.insert(std::make_pair(
"n_walk",&p_n_walk));
2863 p_user_seed.
help=((std::string)
"Seed for multiplier for random ")+
2864 "number generator. If zero is given (the default), then mcmc() "+
2865 "uses time(0) to generate a random seed.";
2866 cl.
par_list.insert(std::make_pair(
"user_seed",&p_user_seed));
2869 p_aff_inv.
help=((std::string)
"If true, then use affine-invariant ")+
2870 "sampling (default false).";
2871 cl.
par_list.insert(std::make_pair(
"aff_inv",&p_aff_inv));
2873 p_table_sequence.
b=&this->table_sequence;
2874 p_table_sequence.
help=((std::string)
"If true, then ensure equal spacing ")+
2875 "between walkers (default true).";
2876 cl.
par_list.insert(std::make_pair(
"table_sequence",&p_table_sequence));
2878 p_store_rejects.
b=&this->store_rejects;
2879 p_store_rejects.
help=((std::string)
"If true, then store MCMC rejections ")+
2881 cl.
par_list.insert(std::make_pair(
"store_rejects",&p_store_rejects));
int verbose
Output control (default 0)
std::string prefix
Prefix for output filenames (default "mcmc")
double vector_lagk_autocorr(size_t n, const vec_t &data, size_t k, double mean)
Lag-k autocorrelation.
A generic MCMC simulation class.
double get(std::string scol, size_t row) const
Get value from row row of column named col. .
int setd_vec_copy(std::string name, const vec_t &v)
Set vector dataset named name with v.
double step_fac
Stepsize factor (default 10.0)
std::function< int(const vec_t &, double, size_t, int, bool, data_t &)> internal_measure_t
Measurement functor type for the parent.
size_t n_params
Number of parameters.
Double parameter for o2scl::cli.
int user_seed
If non-zero, use as the seed for the random number generator (default 0)
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
virtual void write_files(bool sync_write=false)
Write MCMC tables to files.
int set_szt_vec(std::string name, const std::vector< size_t > &v)
Set vector dataset named name with v.
size_t get_nlines() const
Return the number of lines.
bool pd_mode
If true, then use the user-specified proposal distribution.
bool aff_inv
If true, use affine-invariant Monte Carlo.
size_t n_threads
Number of OpenMP threads.
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.
virtual void setup_cli(cli &cl)
Set up the 'cli' object.
virtual void ac_coeffs(size_t ncols, ubmatrix &ac_coeffs)
Compute autocorrelation coefficients.
void close()
Close the file.
Integer parameter for o2scl::cli.
void set_szt(std::string name, size_t u)
Set an unsigned integer named name to value u.
void setd(std::string name, double d)
Set a double named name to value d.
std::vector< data_t > data_arr
Data array.
String parameter for o2scl::cli.
size_t last_write_iters
Total number of MCMC acceptances over all threads at last file write() (default 0) ...
vec_t high_copy
A copy of the upper limits for HDF5 output.
sanity check failed - shouldn't happen
size_t max_iters
If non-zero, the maximum number of MCMC iterations (default 0)
invalid argument supplied by user
std::vector< size_t > n_accept
The number of Metropolis steps which were accepted in each independent chain (summed over all walkers...
void seti(std::string name, int i)
Set an integer named name to value i.
bool always_accept
If true, accept all steps.
int sets_vec(std::string name, const std::vector< std::string > &s)
Set a vector of strings named name.
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.
static const int mcmc_skip
Integer to indicate rejection.
MCMC class with a command-line interface.
void reblock(size_t n_blocks)
Reaverage the data into blocks of a fixed size in order to avoid autocorrelations.
o2scl::mcmc_para_old_table< func_t, fill_t, data_t, vec_t > parent_t
The parent typedef.
virtual void read_prev_results(o2scl_hdf::hdf_file &hf, size_t n_param_loc, std::string name="")
Read previous results (number of threads and walkers must be set first)
double file_update_time
Time between file updates (default 0.0 for no file updates)
void set_proposal_ptrs(prob_vec_t &pv)
Set pointers to proposal distributions.
size_t n_walk_per_thread
Number of walkers per thread (default 1)
requested feature not (yet) implemented
bool prev_read
If true, previous results have been read.
Configurable command-line interface.
size_t n_chains_per_rank
Number of fully independent chains in each MPI rank.
std::string help
Help description.
virtual void post_pointmeas()
Function to run after point evaluation and measurement steps.
void sets(std::string name, std::string s)
Set a string named name to value s.
static const int mcmc_done
Integer to indicate completion.
virtual void mcmc_cleanup()
Cleanup after the MCMC.
virtual void reorder_table()
Reorder the table by thread and walker index.
bool err_nonconv
If true, call the error handler if msolve() or msolve_de() does not converge (default true) ...
std::vector< o2scl::prob_cond_mdim< vec_t > * > prop_dist
Pointer to proposal distribution for each thread.
double mpi_start_time
The MPI starting time (defaults to 0.0)
std::vector< vec_t > current
Current points in parameter space for each walker and each OpenMP thread.
bool table_sequence
If true, ensure sure walkers and OpenMP threads are written to the table with equal spacing between r...
void open_or_create(std::string fname)
Open a file named fname or create if it doesn't already exist.
std::vector< size_t > n_reject
The number of Metropolis steps which were rejected in each independent chain (summed over all walkers...
virtual void initial_points_file_best(std::string fname, size_t n_param_loc, double thresh=1.0e-6, size_t offset=5)
Read initial points from the best points recorded in file named fname.
virtual void ac_lengths(size_t ncols, ubmatrix &ac_coeffs_cols, ubvector &ac_lengths)
Compute autocorrelation lengths.
std::vector< rng_gsl > rg
Random number generators.
std::vector< int > walker_reject_rows
For each walker, record the last row in the table which corresponds to an reject. ...
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
std::string dtos(double x, int prec=6, bool auto_prec=false)
Convert a double to a string.
std::vector< bool > switch_arr
Data switch array for each walker and each OpenMP thread.
String parameter for o2scl::cli.
void set_table(std::shared_ptr< o2scl::table_units<> > &t)
Set the output table.
int set_szt_arr2d_copy(std::string name, size_t r, size_t c, const arr2d_t &a2d)
Set a two-dimensional array dataset named name with m.
void get_chain_sizes(std::vector< size_t > &chain_sizes)
Determine the chain sizes.
std::vector< std::vector< size_t > > ret_value_counts
Return value counters, one vector independent chain.
int mpi_rank
The MPI processor rank.
size_t file_update_iters
Iterations between file updates (default 0 for no file updates)
bool store_rejects
If true, store MCMC rejections in the table.
std::ofstream scr_out
The screen output file.
virtual int mcmc(size_t n_params_local, vec_t &low, vec_t &high, std::vector< func_t > &func, std::vector< fill_t > &fill)
Perform an MCMC simulation.
std::string * str
Parameter.
std::map< std::string, parameter *, std::less< std::string > > par_list
Parameter list.
bool warm_up
If true, we are in the warm up phase.
int mpi_size
The MPI number of processors.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
int setd_arr2d_copy(std::string name, size_t r, size_t c, const arr2d_t &a2d)
Set a two-dimensional array dataset named name with m.
bool meas_for_initial
If true, call the measurement function for the initial point.
std::shared_ptr< o2scl::table_units<> > table
Main data table for Markov chain.
Integer parameter for o2scl::cli.
int table_io_chunk
The number of tables to combine before I/O (default 1)
void vector_out(std::ostream &os, size_t n, const vec_t &v, bool endline=false)
Output the first n elements of a vector to a stream, os.
size_t max_bad_steps
Maximum number of failed steps when generating initial points with affine-invariant sampling (default...
mcmc_para_old_base< func_t, internal_measure_t, data_t, vec_t > parent_t
Type of parent class.
std::vector< std::string > col_names
Column names.
Data table table class with units.
double last_write_time
Time at last file write() (default 0.0)
virtual int mcmc_init()
Initializations before the MCMC.
virtual void best_point(vec_t &best, double w_best, data_t &dat)
Function to run for the best point.
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.
size_t n_walk
Number of walkers for affine-invariant MC or 1 otherwise (default 1)
Store data in an O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$sc...
virtual void set_names_units(std::vector< std::string > names, std::vector< std::string > units)
Set the table names and units.
std::vector< int > walker_accept_rows
For each walker, record the last row in the table which corresponds to an accept. ...
double max_time
Time in seconds (default is 0)
virtual void file_header(o2scl_hdf::hdf_file &hf)
Initial write to HDF5 file.
virtual void unset_proposal()
Go back to random-walk Metropolis with a uniform distribution.
virtual void initial_points_file_last(std::string fname, size_t n_param_loc, size_t offset=5)
Read initial points from the last points recorded in file named fname.
vec_t low_copy
A copy of the lower limits for HDF5 output.
virtual void file_header(o2scl_hdf::hdf_file &hf)
Initial write to HDF5 file.
std::vector< std::string > col_units
Column units.
virtual int add_line(const vec_t &pars, double log_weight, size_t walker_ix, int func_ret, bool mcmc_accept, data_t &dat, size_t i_thread, fill_t &fill)
A measurement function which adds the point to the table.
int open(std::string fname, bool write_access=false, bool err_on_fail=true)
Open a file named fname.
std::vector< size_t > curr_walker
Index of the current walker.
std::vector< std::string > cl_args
The arguments sent to the command-line.
std::string itos(int x)
Convert an integer to a string.
virtual void post_pointmeas()
Function to run after point evaluation and measurement steps.
std::string szttos(size_t x)
Convert a size_t to a string.
virtual void mcmc_cleanup()
Perform cleanup after an MCMC simulation.
void hdf_input(hdf_file &hf, o2scl::table< vec_t > &t, std::string name)
Input a o2scl::table object from a hdf_file.
virtual void initial_points_file_dist(std::string fname, size_t n_param_loc, size_t offset=5)
Read initial points from file named fname, distributing across the chain if necessary.
virtual int fill_line(const vec_t &pars, double log_weight, std::vector< double > &line, data_t &dat, size_t i_walker, fill_t &fill)
Fill line with data for insertion into the table.
std::shared_ptr< o2scl::table_units<> > get_table()
Get the output table.
bool first_write
If true, the HDF5 I/O initial info has been written to the file (set by mcmc() )
void set_proposal(prob_vec_t &pv)
Set the proposal distribution.
virtual int mcmc_init()
MCMC initialization function.
double ai_initial_step
Initial step fraction for affine-invariance sampling walkers (default 0.1)
std::vector< ubvector > initial_points
Initial points in parameter space.