27 #ifndef O2SCL_MCMC_PARA_H 28 #define O2SCL_MCMC_PARA_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/interpm_idw.h> 50 #include <o2scl/vec_stats.h> 51 #include <o2scl/cli.h> 125 template<
class func_t,
class measure_t,
143 std::vector<rng_gsl>
rg;
187 std::cout <<
"Prefix is: " <<
prefix << std::endl;
193 scr_out.setf(std::ios::scientific);
203 scr_out <<
"mcmc (" << it <<
"," << mpi_rank
205 <<
" reject=" <<
n_reject[it] << std::endl;
214 virtual void best_point(vec_t &best,
double w_best, data_t &dat) {
377 MPI_Comm_rank(MPI_COMM_WORLD,&this->mpi_rank);
378 MPI_Comm_size(MPI_COMM_WORLD,&this->mpi_size);
405 virtual int mcmc(
size_t n_params, vec_t &low, vec_t &high,
406 std::vector<func_t> &func, std::vector<measure_t> &meas) {
414 std::cout <<
"mcmc_para::mcmc(): Not enough functions for " 415 << n_threads <<
" threads. Setting n_threads to " 416 << func.size() <<
"." << std::endl;
418 n_threads=func.size();
422 std::cout <<
"mcmc_para::mcmc(): Not enough measurement objects for " 423 << n_threads <<
" threads. Setting n_threads to " 424 << meas.size() <<
"." << std::endl;
426 n_threads=meas.size();
430 if (mpi_start_time==0.0) {
432 mpi_start_time=MPI_Wtime();
434 mpi_start_time=time(0);
436 std::cout <<
"Set start time to: " << mpi_start_time << std::endl;
439 if (initial_points.size()==0) {
441 initial_points.resize(1);
442 initial_points[0].resize(n_params);
443 for(
size_t k=0;k<n_params;k++) {
444 initial_points[0][k]=(low[k]+high[k])/2.0;
449 for(
size_t iip=0;iip<initial_points.size();iip++) {
450 for(
size_t ipar=0;ipar<n_params;ipar++) {
451 if (initial_points[iip][ipar]<low[ipar] ||
452 initial_points[iip][ipar]>high[ipar]) {
456 ") in mcmc_base::mcmc().").c_str(),
465 omp_set_num_threads(n_threads);
468 std::cout <<
"mcmc_para::mcmc(): " 469 << n_threads <<
" threads were requested but the " 470 <<
"-DO2SCL_OPENMP flag was not used during " 471 <<
"compilation. Setting n_threads to 1." 478 std::vector<int> func_ret(n_threads), meas_ret(n_threads);
483 std::cout <<
"mcmc_para::mcmc(): Requested only 1 walker, " 484 <<
"forcing 2 walkers." << std::endl;
487 #ifdef O2SCL_NEVER_DEFINED 493 if (n_walk_per_thread>n_walk) {
494 O2SCL_ERR2(
"More walkers per thread than total walkers ",
497 n_chains_per_rank=n_walk_per_thread*n_threads/
n_walk;
498 if (n_chains_per_walk*n_walk!=n_walk_per_thread*n_threads) {
499 std::cout <<
"mcmc_para::mcmc(): Could not evenly " 500 <<
"organize threads and walkers." << std::endl;
501 std::cout <<
"n_threads: " << n_threads << std::endl;
502 std::cout <<
"n_walk: " << n_walk << std::endl;
503 std::cout <<
"n_walk_per_thread: " 504 << n_walk_per_thread << << std::endl;
505 std::cout <<
"n_chains_per_rank: " 506 << n_chains_per_rank << << std::endl;
517 std::cout <<
"mcmc_para::mcmc(): Requested negative or zero " 518 <<
"step_fac with aff_inv=true. Setting to 2.0." 522 std::cout <<
"mcmc_para::mcmc(): Requested negative or zero " 523 <<
"step_fac. Setting to 10.0." << std::endl;
529 rg.resize(n_threads);
530 unsigned long int seed=time(0);
531 if (this->user_seed!=0) {
535 seed*=(mpi_rank*n_threads+it+1);
536 rg[it].set_seed(seed);
541 n_accept.resize(n_chains_per_rank);
542 n_reject.resize(n_chains_per_rank);
551 if (n_warm_up==0) warm_up=
false;
557 current.resize(ssize);
558 std::vector<double> w_current(ssize);
559 for(
size_t i=0;i<ssize;i++) {
560 current[i].resize(n_params);
565 ret_value_counts.resize(n_chains_per_rank);
566 curr_walker.resize(n_threads);
569 data_arr.resize(2*ssize);
570 switch_arr.resize(ssize);
571 for(
size_t i=0;i<switch_arr.size();i++) switch_arr[i]=
false;
574 std::vector<vec_t> next(n_threads);
576 next[it].resize(n_params);
578 std::vector<double> w_next(n_threads);
581 vec_t best(n_params);
586 std::vector<bool> mcmc_done_flag(n_threads);
588 mcmc_done_flag[it]=
false;
592 std::vector<double> q_prop(n_threads);
599 O2SCL_ERR(
"Function mcmc_init() failed in mcmc_base::mcmc().",
609 scr_out <<
"mcmc: Affine-invariant step, n_params=" 610 << n_params <<
", n_walk=" << n_walk
611 <<
", n_chains_per_rank=" << n_chains_per_rank
612 <<
",\n\tn_walk_per_thread=" << n_walk_per_thread
613 <<
", n_threads=" << n_threads <<
", rank=" 614 << mpi_rank <<
", n_ranks=" 615 << mpi_size << std::endl;
616 }
else if (pd_mode==
true) {
617 scr_out <<
"mcmc: With proposal distribution, n_params=" 618 << n_params <<
", n_threads=" << n_threads <<
", rank=" 619 << mpi_rank <<
", n_ranks=" 620 << mpi_size << std::endl;
622 scr_out <<
"mcmc: Random-walk w/uniform dist., n_params=" 623 << n_params <<
", n_threads=" << n_threads <<
", rank=" 624 << mpi_rank <<
", n_ranks=" 625 << mpi_size << std::endl;
638 #pragma omp parallel default(shared) 647 for(curr_walker[it]=0;curr_walker[it]<n_walk_per_thread &&
648 mcmc_done_flag[it]==
false;curr_walker[it]++) {
651 size_t sindex=n_walk*it+curr_walker[it];
654 size_t ip_index=sindex % initial_points.size();
660 if (initial_points.size()>0) {
663 for(
size_t ipar=0;ipar<n_params;ipar++) {
664 current[sindex][ipar]=initial_points[ip_index][ipar];
668 func_ret[it]=func[it](n_params,current[sindex],
669 w_current[sindex],data_arr[sindex]);
671 if (func_ret[it]==mcmc_done) {
672 mcmc_done_flag[it]=
true;
677 if (func_ret[it]>=0 &&
678 func_ret[it]<((
int)ret_value_counts[it].size())) {
679 ret_value_counts[it][func_ret[it]]++;
681 meas_ret[it]=meas[it](current[sindex],w_current[sindex],
682 curr_walker[it],func_ret[it],
683 true,data_arr[sindex]);
684 if (meas_ret[it]==mcmc_done) {
685 mcmc_done_flag[it]=
true;
691 while (!done && !mcmc_done_flag[it]) {
694 for(
size_t ipar=0;ipar<n_params;ipar++) {
696 current[sindex][ipar]=
697 initial_points[ip_index][ipar]+
698 (rg[it].random()*2.0-1.0)*
699 (high[ipar]-low[ipar])*ai_initial_step;
700 }
while (current[sindex][ipar]>high[ipar] ||
701 current[sindex][ipar]<low[ipar]);
705 func_ret[it]=func[it](n_params,current[sindex],
706 w_current[sindex],data_arr[sindex]);
713 if (func_ret[it]==mcmc_done) {
714 mcmc_done_flag[it]=
true;
719 if (func_ret[it]>=0 &&
720 func_ret[it]<((
int)ret_value_counts[it].size())) {
721 ret_value_counts[it][func_ret[it]]++;
723 if (meas_ret[it]!=mcmc_done) {
724 meas_ret[it]=meas[it](current[sindex],w_current[sindex],
725 curr_walker[it],func_ret[it],
true,
728 mcmc_done_flag[it]=
true;
731 }
else if (init_iters>max_bad_steps) {
743 bool stop_early=
false;
745 if (mcmc_done_flag[it]==
true) {
747 scr_out <<
"mcmc (" << it <<
"," << mpi_rank
748 <<
"): Returned mcmc_done " 749 <<
"(initial; ai)." << std::endl;
763 for(curr_walker[it]=0;curr_walker[it]<
n_walk;
765 size_t sindex=n_walk*it+curr_walker[it];
766 if (w_current[sindex]>w_current[0]) {
768 w_best=w_current[sindex];
772 best=current[best_index];
778 for(curr_walker[it]=0;curr_walker[it]<
n_walk;curr_walker[it]++) {
779 size_t sindex=n_walk*it+curr_walker[it];
780 scr_out.precision(4);
781 scr_out <<
"mcmc (" << it <<
"," << mpi_rank <<
"): i_walk: ";
782 scr_out.width((
int)(1.0+log10((
double)(n_walk-1))));
783 scr_out << curr_walker[it] <<
" log_wgt: " 784 << w_current[sindex] <<
" (initial; ai)" << std::endl;
785 scr_out.precision(6);
799 #pragma omp parallel default(shared) 813 size_t ip_size=initial_points.size();
814 for(
size_t ipar=0;ipar<n_params;ipar++) {
815 current[it][ipar]=initial_points[it % ip_size][ipar];
821 func_ret[it]=func[it](n_params,current[it],w_current[it],
834 if (func_ret[it]==mcmc_done) {
836 scr_out <<
"mcmc (" << it
837 <<
"): Initial point returned mcmc_done." 845 O2SCL_ERR((((std::string)
"Initial weight from thread ")+
847 " vanished in mcmc_para_base::mcmc().").c_str(),
855 #pragma omp parallel default(shared) 863 size_t ip_size=initial_points.size();
867 func_ret[it]=func_ret[it % ip_size];
868 current[it]=current[it % ip_size];
869 w_current[it]=w_current[it % ip_size];
870 data_arr[it]=data_arr[it % ip_size];
874 if (func_ret[it]>=0 &&
875 func_ret[it]<((
int)ret_value_counts[it].size())) {
876 ret_value_counts[it][func_ret[it]]++;
879 meas_ret[it]=meas[it](current[it],w_current[it],0,
880 func_ret[it],
true,data_arr[it]);
881 if (meas_ret[it]==mcmc_done) {
882 mcmc_done_flag[it]=
true;
894 bool stop_early=
false;
896 if (mcmc_done_flag[it]==
true) {
898 scr_out <<
"mcmc (" << it <<
"," << mpi_rank
899 <<
"): Returned mcmc_done " 900 <<
"(initial)." << std::endl;
915 if (w_current[it]<w_best) {
917 w_best=w_current[it];
923 scr_out.precision(4);
926 scr_out << w_current[it] <<
" ";
928 scr_out <<
" (initial)" << std::endl;
929 scr_out.precision(6);
940 std::cout <<
"Press a key and type enter to continue. ";
949 bool main_done=
false;
955 std::vector<double> smove_z(n_threads);
961 curr_walker[it]=mcmc_iters %
n_walk;
966 #pragma omp parallel default(shared) 981 ij=((size_t)(rg[it].random()*((double)n_walk)));
982 }
while (ij==curr_walker[it] || ij>=n_walk);
985 double p=rg[it].random();
987 smove_z[it]=(1.0-2.0*p+2.0*a*p+p*p-2.0*a*p*p+a*a*p*p)/a;
990 for(
size_t i=0;i<n_params;i++) {
991 next[it][i]=current[n_walk*it+ij][i]+
992 smove_z[it]*(current[n_walk*it+curr_walker[it]][i]-
993 current[n_walk*it+ij][i]);
996 }
else if (pd_mode) {
999 q_prop[it]=prop_dist[it]->log_metrop_hast(current[it],next[it]);
1001 if (!std::isfinite(q_prop[it])) {
1002 O2SCL_ERR2(
"Proposal distribution not finite in ",
1009 for(
size_t k=0;k<n_params;k++) {
1010 next[it][k]=current[it][k]+(rg[it].random()*2.0-1.0)*
1011 (high[k]-low[k])/step_fac;
1022 for(
size_t k=0;k<n_params;k++) {
1023 if (next[it][k]<low[k] || next[it][k]>high[k]) {
1026 if (next[it][k]<low[k]) {
1027 std::cout <<
"mcmc (" << it <<
"," 1028 << mpi_rank <<
"): Parameter with index " << k
1029 <<
" and value " << next[it][k]
1030 <<
" smaller than limit " << low[k] << std::endl;
1031 scr_out <<
"mcmc (" << it <<
"," 1032 << mpi_rank <<
"): Parameter with index " << k
1033 <<
" and value " << next[it][k]
1034 <<
" smaller than limit " << low[k] << std::endl;
1036 std::cout <<
"mcmc (" << it <<
"," << mpi_rank
1037 <<
"): Parameter with index " << k
1038 <<
" and value " << next[it][k]
1039 <<
" larger than limit " << high[k] << std::endl;
1040 scr_out <<
"mcmc (" << it <<
"," << mpi_rank
1041 <<
"): Parameter with index " << k
1042 <<
" and value " << next[it][k]
1043 <<
" larger than limit " << high[k] << std::endl;
1051 if (func_ret[it]!=mcmc_skip) {
1052 if (switch_arr[n_walk*it+curr_walker[it]]==
false) {
1053 func_ret[it]=func[it](n_params,next[it],w_next[it],
1054 data_arr[it*n_walk+curr_walker[it]+
1057 func_ret[it]=func[it](n_params,next[it],w_next[it],
1058 data_arr[it*n_walk+curr_walker[it]]);
1060 if (func_ret[it]==mcmc_done) {
1061 mcmc_done_flag[it]=
true;
1063 if (func_ret[it]>=0 &&
1064 func_ret[it]<((
int)ret_value_counts[it].size())) {
1065 ret_value_counts[it][func_ret[it]]++;
1082 if (func_ret[it]==mcmc_done) {
1083 scr_out <<
"mcmc (" << it <<
"," << mpi_rank
1084 <<
"): Returned mcmc_done." 1086 }
else if (func_ret[it]==mcmc_skip && verbose>=3) {
1087 scr_out <<
"mcmc (" << it
1088 <<
"): Parameter(s) out of range: " << std::endl;
1089 scr_out.setf(std::ios::showpos);
1090 for(
size_t k=0;k<n_params;k++) {
1091 scr_out << k <<
" " << low[k] <<
" " 1092 << next[it][k] <<
" " << high[k];
1093 if (next[it][k]<low[k] || next[it][k]>high[k]) {
1096 scr_out << std::endl;
1098 scr_out.unsetf(std::ios::showpos);
1100 func_ret[it]!=mcmc_skip) {
1102 scr_out <<
"mcmc (" << it <<
"," << mpi_rank
1103 <<
"): Function returned failure " 1104 << func_ret[it] <<
" at point ";
1105 for(
size_t k=0;k<n_params;k++) {
1106 scr_out << next[it][k] <<
" ";
1108 scr_out << std::endl;
1119 #pragma omp parallel default(shared) 1128 size_t sindex=n_walk*it+curr_walker[it];
1134 if (always_accept && func_ret[it]==
success) accept=
true;
1137 double r=rg[it].random();
1140 double ai_ratio=pow(smove_z[it],((
double)n_params)-1.0)*
1141 exp(w_next[it]-w_current[sindex]);
1145 }
else if (pd_mode) {
1146 if (r<exp(w_next[it]-w_current[sindex]+q_prop[it])) {
1152 if (r<exp(w_next[it]-w_current[sindex])) {
1166 if (switch_arr[sindex]==
false) {
1167 meas_ret[it]=meas[it](next[it],w_next[it],
1168 curr_walker[it],func_ret[it],
true,
1169 data_arr[sindex+n_threads*
n_walk]);
1171 meas_ret[it]=meas[it](next[it],w_next[it],
1172 curr_walker[it],func_ret[it],
true,
1178 current[sindex]=next[it];
1179 w_current[sindex]=w_next[it];
1180 switch_arr[sindex]=!(switch_arr[sindex]);
1189 if (switch_arr[sindex]==
false) {
1190 meas_ret[it]=meas[it](next[it],w_next[it],
1191 curr_walker[it],func_ret[it],
false,
1192 data_arr[sindex+n_threads*
n_walk]);
1194 meas_ret[it]=meas[it](next[it],w_next[it],
1195 curr_walker[it],func_ret[it],
false,
1214 size_t sindex=n_walk*it+curr_walker[it];
1215 scr_out.precision(4);
1216 scr_out <<
"mcmc (" << it <<
"," << mpi_rank <<
"): iter: ";
1217 scr_out.width((
int)(1.0+log10((
double)(n_params-1))));
1218 scr_out << mcmc_iters <<
" i_walk: " 1219 << curr_walker[it] <<
" log_wgt: " 1220 << w_current[sindex] << std::endl;
1221 scr_out.precision(6);
1230 if (switch_arr[n_walk*it+curr_walker[it]]==
false) {
1231 best_point(best,w_best,data_arr[curr_walker[it]+n_walk*it+
1234 best_point(best,w_best,data_arr[curr_walker[it]+n_walk*it]);
1242 if (meas_ret[it]==mcmc_done || func_ret[it]==mcmc_done) {
1247 O2SCL_ERR((((std::string)
"Measurement function returned ")+
1249 " in mcmc_para_base::mcmc().").c_str(),
1258 if (main_done==
false) {
1262 if (warm_up && mcmc_iters==n_warm_up) {
1270 scr_out <<
"mcmc: Finished warmup." << std::endl;
1277 std::cout <<
"Press a key and type enter to continue. ";
1283 if (main_done==
false && warm_up==
false && max_iters>0 &&
1284 mcmc_iters==max_iters) {
1286 scr_out <<
"mcmc: Stopping because number of iterations (" 1287 << mcmc_iters <<
") equal to max_iters (" << max_iters
1288 <<
")." << std::endl;
1293 if (main_done==
false) {
1300 if (max_time>0.0 && elapsed>max_time) {
1302 scr_out <<
"mcmc: Stopping because elapsed (" << elapsed
1303 <<
") > max_time (" << max_time <<
")." 1325 virtual int mcmc(
size_t n_params, vec_t &low, vec_t &high,
1326 func_t &func, measure_t &meas) {
1329 omp_set_num_threads(n_threads);
1333 std::vector<func_t> vf(n_threads);
1334 std::vector<measure_t> vm(n_threads);
1339 return mcmc(n_params,low,high,func,meas);
1348 prop_dist.resize(pv.size());
1349 for(
size_t i=0;i<pv.size();i++) {
1350 prop_dist[i]=&pv[i];
1413 template<
class func_t,
class fill_t,
class data_t,
class vec_t=ubvector>
1415 std::function<int(const vec_t &,double,size_t,int,bool,data_t &)>,
1421 typedef std::function<int(const vec_t &,double,size_t,int,bool,data_t &)>
1437 std::shared_ptr<o2scl::table_units<> >
table;
1454 table->new_column(
"rank");
1455 table->new_column(
"thread");
1456 table->new_column(
"walker");
1457 table->new_column(
"mult");
1458 table->new_column(
"log_wgt");
1459 for(
size_t i=0;i<col_names.size();i++) {
1460 table->new_column(col_names[i]);
1461 if (col_units[i].length()>0) {
1462 table->set_unit(col_names[i],col_units[i]);
1468 walker_accept_rows[i]=-1;
1470 walker_reject_rows.resize(this->
n_walk*this->n_threads);
1472 walker_reject_rows[i]=-1;
1476 std::cout <<
"mcmc: Table column names and units: " << std::endl;
1477 for(
size_t i=0;i<table->get_ncolumns();i++) {
1478 std::cout << table->get_column_name(i) <<
" " 1479 << table->get_unit(table->get_column_name(i)) << std::endl;
1485 last_write_time=MPI_Wtime();
1487 last_write_time=time(0);
1490 return parent_t::mcmc_init();
1496 std::vector<double> &line, data_t &dat,
1497 size_t i_walker, fill_t &fill) {
1500 size_t i_thread=omp_get_thread_num();
1508 line.push_back(i_thread);
1510 line.push_back(i_walker);
1512 line.push_back(1.0);
1513 line.push_back(log_weight);
1514 for(
size_t i=0;i<pars.size();i++) {
1515 line.push_back(pars[i]);
1517 int tempi=fill(pars,log_weight,line,dat);
1583 this->
scr_out <<
"mcmc: Start write_files(). mpi_rank: " 1585 << this->
mpi_size <<
" table_io_chunk: " 1586 << table_io_chunk << std::endl;
1589 std::vector<o2scl::table_units<> > tab_arr;
1590 bool rank_sent=
false;
1593 if (table_io_chunk>1) {
1594 if (this->
mpi_rank%table_io_chunk==0) {
1596 for(
int i=0;i<table_io_chunk-1;i++) {
1600 tab_arr.push_back(t);
1601 o2scl_table_mpi_recv(child,tab_arr[tab_arr.size()-1]);
1607 o2scl_table_mpi_send(*table,parent);
1616 int tag=0, buffer=0;
1617 if (sync_write && this->
mpi_size>1 &&
1619 MPI_Recv(&buffer,1,MPI_INT,this->
mpi_rank-table_io_chunk,
1620 tag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
1628 if (first_write==
false) {
1644 hf.
set_szt(
"n_params",this->n_params);
1646 hf.
seti(
"allow_estimates",allow_estimates);
1662 hf.
seti(
"n_tables",tab_arr.size()+1);
1663 if (rank_sent==
false) {
1664 hdf_output(hf,*table,
"markov_chain_0");
1666 for(
size_t i=0;i<tab_arr.size();i++) {
1667 std::string name=((std::string)
"markov_chain_")+
szttos(i+1);
1668 hdf_output(hf,tab_arr[i],name);
1674 if (sync_write && this->
mpi_size>1 &&
1676 MPI_Send(&buffer,1,MPI_INT,this->
mpi_rank+table_io_chunk,
1677 tag,MPI_COMM_WORLD);
1682 this->
scr_out <<
"mcmc: Done write_files()." << std::endl;
1696 allow_estimates=
false;
1698 file_update_iters=0;
1699 file_update_time=0.0;
1701 store_rejects=
false;
1702 table_sequence=
true;
1710 std::vector<std::string> units) {
1735 int tag=0, buffer=0;
1737 MPI_Recv(&buffer,1,MPI_INT,this->
mpi_rank-1,
1738 tag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
1750 MPI_Send(&buffer,1,MPI_INT,this->
mpi_rank+1,
1751 tag,MPI_COMM_WORLD);
1759 std::cout <<
"Initial points: Finding last " << n_points
1760 <<
" points from file named " 1761 << fname <<
" ." << std::endl;
1766 for(
size_t it=0;it<this->
n_threads;it++) {
1767 for(
size_t iw=0;iw<this->
n_walk;iw++) {
1770 size_t windex=it*this->n_walk+iw;
1773 for(
int row=table->get_nlines()-1;row>=0 && found==
false;row--) {
1774 if (table->get(
"walker",row)==iw &&
1775 table->get(
"thread",row)==it &&
1776 table->get(
"mult",row)>0.5) {
1780 std::cout <<
"Function initial_point_file_last():\n\tit: " 1782 <<
" iw: " << iw <<
" row: " 1783 << row <<
" log_wgt: " << table->get(
"log_wgt",row)
1788 for(
size_t ip=0;ip<n_param_loc;ip++) {
1794 O2SCL_ERR(
"Function initial_points_file_last() failed.",
1821 int tag=0, buffer=0;
1823 MPI_Recv(&buffer,1,MPI_INT,this->
mpi_rank-1,
1824 tag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
1836 MPI_Send(&buffer,1,MPI_INT,this->
mpi_rank+1,
1837 tag,MPI_COMM_WORLD);
1845 std::cout <<
"Initial points: Finding last " << n_points
1846 <<
" points from file named " 1847 << fname <<
" ." << std::endl;
1852 size_t nlines=table->get_nlines();
1853 size_t decrement=nlines/n_points;
1854 if (decrement<1) decrement=1;
1857 for(
size_t k=0;k<n_points;k++) {
1861 std::cout <<
"Function initial_point_file_dist():\n\trow: " 1862 << row <<
" log_wgt: " << table->get(
"log_wgt",row)
1867 for(
size_t ip=0;ip<n_param_loc;ip++) {
1887 double thresh=1.0e-6,
1895 int tag=0, buffer=0;
1897 MPI_Recv(&buffer,1,MPI_INT,this->
mpi_rank-1,
1898 tag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
1910 MPI_Send(&buffer,1,MPI_INT,this->
mpi_rank+1,
1911 tag,MPI_COMM_WORLD);
1919 std::cout <<
"Initial points: Finding best " << n_points
1920 <<
" unique points from file named " 1921 << fname <<
" ." << std::endl;
1924 typedef std::map<double,int,std::greater<double> > map_t;
1928 for(
size_t k=0;k<table->get_nlines();k++) {
1929 m.insert(std::make_pair(table->get(
"log_wgt",k),k));
1939 for(map_t::iterator mit=m.begin();mit!=m.end();mit++) {
1940 map_t::iterator mit2=mit;
1942 if (mit2!=m.end()) {
1943 if (fabs(mit->first-mit2->first)<thresh) {
1945 std::cout <<
"Removing duplicate weights: " 1946 << mit->first <<
" " << mit2->first << std::endl;
1955 }
while (found==
true);
1958 if (m.size()<n_points) {
1959 O2SCL_ERR2(
"Could not find enough points in file in ",
1960 "mcmc_para::initial_points_file_best().",
1966 map_t::iterator mit=m.begin();
1967 for(
size_t k=0;k<n_points;k++) {
1968 int row=mit->second;
1970 std::cout <<
"Initial point " << k <<
" at row " 1971 << row <<
" has log_weight= " 1972 << table->get(
"log_wgt",row) << std::endl;
1975 for(
size_t ip=0;ip<n_param_loc;ip++) {
1991 virtual int mcmc(
size_t n_params_local,
1992 vec_t &low, vec_t &high, std::vector<func_t> &func,
1993 std::vector<fill_t> &fill) {
1995 n_params=n_params_local;
2011 std::vector<internal_measure_t> meas(this->
n_threads);
2012 for(
size_t it=0;it<this->
n_threads;it++) {
2014 (std::mem_fn<
int(
const vec_t &,
double,
size_t,
int,
bool,
2015 data_t &,
size_t, fill_t &)>
2017 std::placeholders::_2,std::placeholders::_3,
2018 std::placeholders::_4,std::placeholders::_5,
2019 std::placeholders::_6,it,std::ref(fill[it]));
2022 return parent_t::mcmc(n_params,low,high,func,meas);
2047 chain_sizes.resize(ntot);
2049 for(
size_t it=0;it<this->
n_threads;it++) {
2050 for(
size_t iw=0;iw<this->
n_walk;iw++) {
2051 size_t ix=it*this->n_walk+iw;
2054 for(
size_t j=istart;j<table->get_nlines();j+=ntot) {
2055 if (table->get(
"mult",j)>0.5) chain_sizes[ix]++;
2066 virtual int add_line(
const vec_t &pars,
double log_weight,
2067 size_t walker_ix,
int func_ret,
2068 bool mcmc_accept, data_t &dat,
2069 size_t i_thread, fill_t &fill) {
2072 size_t windex=i_thread*this->
n_walk+walker_ix;
2081 if ((mcmc_accept || store_rejects) && walker_accept_rows[windex]<0) {
2084 if (table_sequence) {
2085 if (walker_accept_rows[windex]>walker_reject_rows[windex]) {
2086 next_row=walker_accept_rows[windex]+ntot;
2088 next_row=walker_reject_rows[windex]+ntot;
2091 if (walker_accept_rows[windex]>walker_reject_rows[windex]) {
2092 next_row=walker_accept_rows[windex]+1;
2094 next_row=walker_reject_rows[windex]+1;
2096 while (next_row<((
int)table->get_nlines()) &&
2097 table->get(
"mult",next_row)>0.0) {
2107 #pragma omp critical (o2scl_mcmc_para_table_add_line) 2116 if (next_row>=((
int)table->get_nlines())) {
2117 size_t istart=table->get_nlines();
2119 table->set_nlines(table->get_nlines()+ntot);
2122 for(
size_t i=0;i<this->
n_walk;i++) {
2123 table->set(
"rank",istart+j*this->n_walk+i,this->
mpi_rank);
2124 table->set(
"thread",istart+j*this->n_walk+i,j);
2125 table->set(
"walker",istart+j*this->n_walk+i,i);
2126 table->set(
"mult",istart+j*this->n_walk+i,0.0);
2127 table->set(
"log_wgt",istart+j*this->n_walk+i,0.0);
2133 if (func_ret==0 && (mcmc_accept || store_rejects)) {
2135 if (next_row>=((
int)(table->get_nlines()))) {
2139 std::vector<double> line;
2140 int fret=fill_line(pars,log_weight,line,dat,walker_ix,fill);
2144 if (store_rejects && mcmc_accept==
false) {
2160 if (line.size()!=table->get_ncolumns()) {
2161 std::cout <<
"line: " << line.size() <<
" columns: " 2162 << table->get_ncolumns() << std::endl;
2163 for(
size_t k=0;k<table->get_ncolumns() || k<line.size();k++) {
2164 std::cout << k <<
". ";
2165 if (k<table->get_ncolumns()) {
2166 std::cout << table->get_column_name(k) <<
" ";
2167 std::cout << table->get_unit(table->get_column_name(k)) <<
" ";
2169 if (k<line.size()) std::cout << line[k] <<
" ";
2170 std::cout << std::endl;
2172 O2SCL_ERR(
"Table misalignment in mcmc_para_table::add_line().",
2177 table->set_row(((
size_t)next_row),line);
2181 this->
scr_out <<
"mcmc: Setting data at row " << next_row
2183 for(
size_t k=0;k<line.size();k++) {
2185 this->
scr_out << table->get_column_name(k) <<
" ";
2186 this->
scr_out << table->get_unit(table->get_column_name(k));
2187 this->
scr_out <<
" " << line[k] << std::endl;
2198 double mult_old=table->get(
"mult",walker_accept_rows[windex]);
2199 table->set(
"mult",walker_accept_rows[windex],mult_old+1.0);
2201 this->
scr_out <<
"mcmc: Updating mult of row " 2202 << walker_accept_rows[windex]
2203 <<
" from " << mult_old <<
" to " 2204 << mult_old+1.0 << std::endl;
2215 walker_accept_rows[windex]=next_row;
2216 }
else if (store_rejects && func_ret==0) {
2217 walker_reject_rows[windex]=next_row;
2231 if (file_update_iters>0) {
2232 size_t total_iters=0;
2236 if (total_iters>=last_write_iters+file_update_iters) {
2238 this->
scr_out <<
"mcmc: Writing to file. total_iters: " 2239 << total_iters <<
" file_update_iters: " 2240 << file_update_iters <<
" last_write_iters: " 2241 << last_write_iters << std::endl;
2244 last_write_iters=total_iters;
2248 if (updated==
false && file_update_time>0.0) {
2250 double elapsed=MPI_Wtime()-last_write_time;
2252 double elapsed=time(0)-last_write_time;
2254 if (elapsed>file_update_time) {
2256 this->
scr_out <<
"mcmc: Writing to file. elapsed: " 2257 << elapsed <<
" file_update_time: " 2258 << file_update_time <<
" last_write_time: " 2259 << last_write_time << std::endl;
2263 last_write_time=MPI_Wtime();
2265 last_write_time=time(0);
2281 for(i=table->get_nlines()-1;i>=0 && done==
false;i--) {
2283 if (table->get(
"mult",i)<1.0e-10) {
2287 if (i+2<((
int)table->get_nlines())) {
2288 table->set_nlines(i+2);
2293 return parent_t::mcmc_cleanup();
2299 std::vector<size_t> chain_sizes;
2300 get_chain_sizes(chain_sizes);
2301 size_t min_size=chain_sizes[0];
2302 for(
size_t i=1;i<chain_sizes.size();i++) {
2303 if (chain_sizes[i]<min_size) min_size=chain_sizes[i];
2305 size_t N_max=min_size/2;
2306 ac_coeffs.resize(ncols,N_max-1);
2307 for(
size_t i=0;i<ncols;i++) {
2308 for(
size_t ell=1;ell<N_max;ell++) {
2309 ac_coeffs(i,ell-1)=0.0;
2314 size_t cstart=table->lookup_column(
"log_wgt")+1;
2315 for(
size_t i=0;i<ncols;i++) {
2317 for(
size_t k=0;k<this->
n_walk;k++) {
2318 size_t tindex=j*this->n_walk+k;
2319 for(
size_t ell=1;ell<N_max;ell++) {
2320 const double &x=(*table)[cstart+i][table_row];
2321 double mean=o2scl::vector_mean<const double *>
2322 (chain_sizes[tindex]+1,&x);
2324 <
const double *>(chain_sizes[tindex]+1,&x,ell,mean);
2326 table_row+=chain_sizes[tindex]+1;
2329 for(
size_t ell=1;ell<N_max;ell++) {
2330 ac_coeffs(i,ell-1)/=((double)n_tot);
2339 ubvector &ac_lengths) {
2340 size_t N_max=ac_coeffs_cols.size2();
2341 ac_lengths.resize(ncols);
2342 for(
size_t icol=0;icol<ncols;icol++) {
2343 std::vector<double> tau(N_max);
2344 for(
size_t i=5;i<N_max;i++) {
2346 for(
size_t j=0;j<i;j++) {
2347 sum+=ac_coeffs_cols(icol,j);
2350 std::cout << tau[i] <<
" " << ((double)i)/5.0 << std::endl;
2352 std::cout << std::endl;
2362 std::shared_ptr<o2scl::table_units<> > table2=
2366 for(
size_t j=0;j<this->
n_walk;j++) {
2367 std::string func=std::string(
"abs(walker-")+
o2scl::szttos(j)+
2369 table->copy_rows(func,*table2);
2391 for(
size_t it=0;it<this->
n_threads;it++) {
2393 size_t n=table->get_nlines();
2395 O2SCL_ERR2(
"Cannot reblock. Not enough data in ",
2398 size_t n_block=n/n_blocks;
2399 size_t m=table->get_ncolumns();
2400 for(
size_t j=0;j<n_blocks;j++) {
2403 for(
size_t i=0;i<m;i++) {
2406 for(
size_t k=j*n_block;k<(j+1)*n_block;k++) {
2407 mult+=(*table)[
"mult"][k];
2408 for(
size_t i=1;i<m;i++) {
2409 dat[i]+=(*table)[i][k]*(*table)[
"mult"][k];
2412 table->set(
"mult",j,mult);
2413 for(
size_t i=1;i<m;i++) {
2415 table->set(i,j,dat[i]);
2418 table->set_nlines(n_blocks);
2434 template<
class func_t,
class fill_t,
class data_t,
class vec_t=ubvector>
2465 hf.
sets_vec(
"cl_args",this->cl_args);
2513 p_file_update_iters.
s=&this->file_update_iters;
2514 p_file_update_iters.
help=((std::string)
"Number of MCMC successes ")+
2515 "between file updates (default 0 for no file updates).";
2516 cl.
par_list.insert(std::make_pair(
"file_update_iters",
2517 &p_file_update_iters));
2519 p_file_update_time.
d=&this->file_update_time;
2520 p_file_update_time.
help=((std::string)
"Time ")+
2521 "between file updates (default 0.0 for no file updates).";
2522 cl.
par_list.insert(std::make_pair(
"file_update_time",
2523 &p_file_update_time));
2534 p_max_time.
help=((std::string)
"Maximum run time in seconds ")+
2535 "(default 86400 sec or 1 day).";
2536 cl.
par_list.insert(std::make_pair(
"max_time",&p_max_time));
2539 p_max_iters.
help=((std::string)
"If non-zero, limit the number of ")+
2540 "iterations to be less than the specified number (default zero).";
2541 cl.
par_list.insert(std::make_pair(
"max_iters",&p_max_iters));
2544 p_prefix.
help=
"Output file prefix (default 'mcmc\').";
2545 cl.
par_list.insert(std::make_pair(
"prefix",&p_prefix));
2555 p_step_fac.
help=((std::string)
"MCMC step factor. The step size for ")+
2556 "each variable is taken as the difference between the high and low "+
2557 "limits divided by this factor (default 10.0). This factor can "+
2558 "be increased if the acceptance rate is too small, but care must "+
2559 "be taken, e.g. if the conditional probability is multimodal. If "+
2560 "this step size is smaller than 1.0, it is reset to 1.0 .";
2561 cl.
par_list.insert(std::make_pair(
"step_fac",&p_step_fac));
2564 p_n_warm_up.
help=((std::string)
"Minimum number of warm up iterations ")+
2566 cl.
par_list.insert(std::make_pair(
"n_warm_up",&p_n_warm_up));
2569 p_verbose.
help=((std::string)
"MCMC verbosity parameter ")+
2571 cl.
par_list.insert(std::make_pair(
"mcmc_verbose",&p_verbose));
2574 p_max_bad_steps.
help=((std::string)
"Maximum number of bad steps ")+
2576 cl.
par_list.insert(std::make_pair(
"max_bad_steps",&p_max_bad_steps));
2579 p_n_walk.
help=((std::string)
"Number of walkers ")+
2581 cl.
par_list.insert(std::make_pair(
"n_walk",&p_n_walk));
2584 p_user_seed.
help=((std::string)
"Seed for multiplier for random ")+
2585 "number generator. If zero is given (the default), then mcmc() "+
2586 "uses time(0) to generate a random seed.";
2587 cl.
par_list.insert(std::make_pair(
"user_seed",&p_user_seed));
2590 p_aff_inv.
help=((std::string)
"If true, then use affine-invariant ")+
2591 "sampling (default false).";
2592 cl.
par_list.insert(std::make_pair(
"aff_inv",&p_aff_inv));
2594 p_table_sequence.
b=&this->table_sequence;
2595 p_table_sequence.
help=((std::string)
"If true, then ensure equal spacing ")+
2596 "between walkers (default true).";
2597 cl.
par_list.insert(std::make_pair(
"table_sequence",&p_table_sequence));
2599 p_store_rejects.
b=&this->store_rejects;
2600 p_store_rejects.
help=((std::string)
"If true, then store MCMC rejections ")+
2602 cl.
par_list.insert(std::make_pair(
"store_rejects",&p_store_rejects));
void get_chain_sizes(std::vector< size_t > &chain_sizes)
Determine the chain sizes.
double vector_lagk_autocorr(size_t n, const vec_t &data, size_t k, double mean)
Lag-k autocorrelation.
std::shared_ptr< o2scl::table_units<> > table
Main data table for Markov chain.
virtual void post_pointmeas()
Function to run after point evaluation and measurement steps.
std::vector< data_t > data_arr
Data array.
int setd_vec_copy(std::string name, const vec_t &v)
Set vector dataset named name with v.
std::string prefix
Prefix for output filenames (default "mcmc")
static const int mcmc_done
Integer to indicate completion.
std::vector< std::string > cl_args
The arguments sent to the command-line.
Double parameter for o2scl::cli.
double ai_initial_step
Initial step fraction for affine-invariance sampling walkers (default 0.1)
bool store_rejects
If true, store MCMC rejections in the table.
Multi-dimensional interpolation by inverse distance weighting.
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
int set_szt_vec(std::string name, const std::vector< size_t > &v)
Set vector dataset named name with v.
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 allow_estimates
If true, allow estimates of the weight (default false)
bool err_nonconv
If true, call the error handler if msolve() or msolve_de() does not converge (default true) ...
int verbose
Output control (default 0)
std::vector< vec_t > current
Current points in parameter space for each walker and each OpenMP thread.
double last_write_time
Time at last file write() (default 0.0)
size_t file_update_iters
Iterations between file updates (default 0 for no file updates)
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.
String parameter for o2scl::cli.
bool aff_inv
If true, use affine-invariant Monte Carlo.
size_t max_iters
If non-zero, the maximum number of MCMC iterations (default 0)
sanity check failed - shouldn't happen
std::vector< bool > switch_arr
Data switch array for each walker and each OpenMP thread.
std::vector< std::string > col_names
Column names.
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.
invalid argument supplied by user
virtual void ac_lengths(size_t ncols, ubmatrix &ac_coeffs_cols, ubvector &ac_lengths)
Compute autocorrelation lengths.
vec_t high_copy
A copy of the upper limits for HDF5 output.
std::vector< std::vector< size_t > > ret_value_counts
Return value counters, one vector independent chain.
size_t n_chains_per_rank
Number of fully independent chains in each MPI rank.
void seti(std::string name, int i)
Set an integer named name to value i.
std::vector< size_t > curr_walker
Index of the current walker.
MCMC class with a command-line interface.
double step_fac
Stepsize factor (default 10.0)
int sets_vec(std::string name, const std::vector< std::string > &s)
Set a vector of strings named name.
double file_update_time
Time between file updates (default 0.0 for no file updates)
size_t last_write_iters
Total number of MCMC acceptances over all threads at last file write() (default 0) ...
virtual void unset_proposal()
Go back to random-walk Metropolis with a uniform distribution.
std::function< int(const vec_t &, double, size_t, int, bool, data_t &)> internal_measure_t
Measurement functor type for the parent.
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.
bool always_accept
If true, accept all steps.
std::vector< ubvector > initial_points
Initial points in parameter space.
bool first_write
If true, the HDF5 I/O initial info has been written to the file (set by mcmc() )
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.
int table_io_chunk
The number of tables to combine before I/O (default 1)
Configurable command-line interface.
std::string help
Help description.
size_t n_walk
Number of walkers for affine-invariant MC or 1 otherwise (default 1)
virtual void mcmc_cleanup()
Perform cleanup after an MCMC simulation.
void sets(std::string name, std::string s)
Set a string named name to value s.
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.
std::vector< std::string > col_units
Column units.
void set_proposal(prob_vec_t &pv)
Set the proposal distribution.
double mpi_start_time
The MPI starting time (defaults to 0.0)
void open_or_create(std::string fname)
Open a file named fname or create if it doesn't already exist.
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.
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.
int mpi_size
The MPI number of processors.
bool table_sequence
If true, ensure sure walkers and OpenMP threads are written to the table with equal spacing between r...
double max_time
Time in seconds (default is 0)
#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.
int user_seed
If non-zero, use as the seed for the random number generator (default 0)
String parameter for o2scl::cli.
void reblock(size_t n_blocks)
Reaverage the data into blocks of a fixed size in order to avoid autocorrelations.
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.
std::vector< size_t > n_reject
The number of Metropolis steps which were rejected in each independent chain (summed over all walkers...
virtual void set_names_units(std::vector< std::string > names, std::vector< std::string > units)
Set the table names and units.
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.
std::vector< size_t > n_accept
The number of Metropolis steps which were accepted in each independent chain (summed over all walkers...
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
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.
bool pd_mode
If true, then use the user-specified proposal distribution.
virtual void reorder_table()
Reorder the table by thread and walker index.
size_t n_params
Number of parameters.
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.
size_t n_walk_per_thread
Number of walkers per thread (default 1)
virtual void mcmc_cleanup()
Cleanup after the MCMC.
std::vector< int > walker_reject_rows
For each walker, record the last row in the table which corresponds to an reject. ...
std::vector< o2scl::prob_cond_mdim< vec_t > * > prop_dist
Pointer to proposal distribution for each thread.
Integer parameter for o2scl::cli.
virtual int mcmc_init()
MCMC initialization function.
Data table table class with units.
virtual int mcmc_init()
Initializations before the MCMC.
int mpi_rank
The MPI processor rank.
virtual void file_header(o2scl_hdf::hdf_file &hf)
Initial write to HDF5 file.
static const int mcmc_skip
Integer to indicate rejection.
virtual void setup_cli(cli &cl)
Set up the 'cli' object.
size_t max_bad_steps
Maximum number of failed steps when generating initial points with affine-invariant sampling (default...
vec_t low_copy
A copy of the lower limits for HDF5 output.
std::vector< int > walker_accept_rows
For each walker, record the last row in the table which corresponds to an accept. ...
Store data in an O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$sc...
virtual void best_point(vec_t &best, double w_best, data_t &dat)
Function to run for the best point.
interpm_idw< double * > esti
Likelihood estimator.
int open(std::string fname, bool write_access=false, bool err_on_fail=true)
Open a file named fname.
size_t n_threads
Number of OpenMP threads.
mcmc_para_base< func_t, internal_measure_t, data_t, vec_t > parent_t
Type of parent class.
std::string itos(int x)
Convert an integer to a string.
A generic MCMC simulation class.
std::vector< rng_gsl > rg
Random number generators.
virtual void file_header(o2scl_hdf::hdf_file &hf)
Initial write to HDF5 file.
void set_table(std::shared_ptr< o2scl::table_units<> > &t)
Set the output table.
std::string szttos(size_t x)
Convert a size_t to a string.
virtual void write_files(bool sync_write=false)
Write MCMC tables to files.
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.
void hdf_input(hdf_file &hf, o2scl::table< vec_t > &t, std::string name)
Input a o2scl::table object from a hdf_file.
std::shared_ptr< o2scl::table_units<> > get_table()
Get the output table.
std::ofstream scr_out
The screen output file.
virtual void ac_coeffs(size_t ncols, ubmatrix &ac_coeffs)
Compute autocorrelation coefficients.
virtual void post_pointmeas()
Function to run after point evaluation and measurement steps.