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/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::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) {
394 MPI_Comm_rank(MPI_COMM_WORLD,&this->mpi_rank);
395 MPI_Comm_size(MPI_COMM_WORLD,&this->mpi_size);
401 meas_for_initial=
true;
423 virtual int mcmc(
size_t n_params, vec_t &low, vec_t &high,
424 std::vector<func_t> &func, std::vector<measure_t> &meas) {
430 if (func.size()==0 || meas.size()==0) {
431 O2SCL_ERR2(
"Size of 'func' or 'meas' array is zero in ",
436 std::cout <<
"mcmc_para::mcmc(): Not enough functions for " 437 << n_threads <<
" threads. Setting n_threads to " 438 << func.size() <<
"." << std::endl;
440 n_threads=func.size();
444 std::cout <<
"mcmc_para::mcmc(): Not enough measurement objects for " 445 << n_threads <<
" threads. Setting n_threads to " 446 << meas.size() <<
"." << std::endl;
448 n_threads=meas.size();
452 if (mpi_start_time==0.0) {
454 mpi_start_time=MPI_Wtime();
456 mpi_start_time=time(0);
460 if (initial_points.size()==0) {
463 initial_points.resize(1);
464 initial_points[0].resize(n_params);
465 for(
size_t k=0;k<n_params;k++) {
466 initial_points[0][k]=(low[k]+high[k])/2.0;
473 for(
size_t iip=0;iip<initial_points.size();iip++) {
474 for(
size_t ipar=0;ipar<n_params;ipar++) {
475 if (initial_points[iip][ipar]<low[ipar] ||
476 initial_points[iip][ipar]>high[ipar]) {
482 ") in mcmc_base::mcmc().").c_str(),
489 for(
size_t i=0;i<initial_points.size();i++) {
490 for(
size_t k=0;k<initial_points[i].size();k++) {
491 if (!std::isfinite(initial_points[i][k])) {
496 for(
size_t j=i+1;j<initial_points.size();j++) {
498 for(
size_t k=0;k<initial_points[i].size();k++) {
499 if (initial_points[i][k]!=initial_points[j][k]) {
504 std::cout.setf(std::ios::scientific);
505 std::cout << i <<
" ";
507 std::cout << j <<
" ";
519 omp_set_num_threads(n_threads);
522 std::cout <<
"mcmc_para::mcmc(): " 523 << n_threads <<
" threads were requested but the " 524 <<
"-DO2SCL_OPENMP flag was not used during " 525 <<
"compilation. Setting n_threads to 1." 532 std::vector<int> func_ret(n_threads), meas_ret(n_threads);
537 std::cout <<
"mcmc_para::mcmc(): Requested only 1 walker, " 538 <<
"forcing 2 walkers." << std::endl;
541 #ifdef O2SCL_NEVER_DEFINED 547 if (n_walk_per_thread>n_walk) {
548 O2SCL_ERR2(
"More walkers per thread than total walkers ",
551 n_chains_per_rank=n_walk_per_thread*n_threads/
n_walk;
552 if (n_chains_per_walk*n_walk!=n_walk_per_thread*n_threads) {
553 std::cout <<
"mcmc_para::mcmc(): Could not evenly " 554 <<
"organize threads and walkers." << std::endl;
555 std::cout <<
"n_threads: " << n_threads << std::endl;
556 std::cout <<
"n_walk: " << n_walk << std::endl;
557 std::cout <<
"n_walk_per_thread: " 558 << n_walk_per_thread << << std::endl;
559 std::cout <<
"n_chains_per_rank: " 560 << n_chains_per_rank << << std::endl;
571 std::cout <<
"mcmc_para::mcmc(): Requested negative or zero " 572 <<
"step_fac with aff_inv=true.\nSetting to 2.0." 576 std::cout <<
"mcmc_para::mcmc(): Requested negative or zero " 577 <<
"step_fac. Setting to 10.0." << std::endl;
583 rg.resize(n_threads);
584 unsigned long int seed=time(0);
585 if (this->user_seed!=0) {
589 seed*=(mpi_rank*n_threads+it+1);
590 rg[it].set_seed(seed);
595 n_accept.resize(n_chains_per_rank);
596 n_reject.resize(n_chains_per_rank);
605 if (n_warm_up==0) warm_up=
false;
611 current.resize(ssize);
612 std::vector<double> w_current(ssize);
613 for(
size_t i=0;i<ssize;i++) {
614 current[i].resize(n_params);
619 curr_walker.resize(n_threads);
626 data_arr.resize(2*ssize);
627 switch_arr.resize(ssize);
628 for(
size_t i=0;i<switch_arr.size();i++) switch_arr[i]=
false;
631 std::vector<vec_t> next(n_threads);
633 next[it].resize(n_params);
635 std::vector<double> w_next(n_threads);
638 vec_t best(n_params);
643 std::vector<bool> mcmc_done_flag(n_threads);
645 mcmc_done_flag[it]=
false;
649 std::vector<double> q_prop(n_threads);
656 O2SCL_ERR(
"Function mcmc_init() failed in mcmc_base::mcmc().",
667 scr_out <<
"mcmc: Affine-invariant step, n_params=" 668 << n_params <<
", n_walk=" << n_walk
669 <<
", n_chains_per_rank=" << n_chains_per_rank
670 <<
",\n\tn_walk_per_thread=" << n_walk_per_thread
671 <<
", n_threads=" << n_threads <<
", rank=" 672 << mpi_rank <<
", n_ranks=" 673 << mpi_size << std::endl;
674 }
else if (pd_mode==
true) {
675 scr_out <<
"mcmc: With proposal distribution, n_params=" 676 << n_params <<
", n_threads=" << n_threads <<
", rank=" 677 << mpi_rank <<
", n_ranks=" 678 << mpi_size << std::endl;
680 scr_out <<
"mcmc: Random-walk w/uniform dist., n_params=" 681 << n_params <<
", n_threads=" << n_threads <<
", rank=" 682 << mpi_rank <<
", n_ranks=" 683 << mpi_size << std::endl;
685 scr_out <<
"Set start time to: " << mpi_start_time << std::endl;
694 #pragma omp parallel default(shared) 703 for(curr_walker[it]=0;curr_walker[it]<n_walk_per_thread &&
704 mcmc_done_flag[it]==
false;curr_walker[it]++) {
707 size_t sindex=n_walk*it+curr_walker[it];
710 size_t ip_index=sindex % initial_points.size();
718 if (sindex<initial_points.size()) {
721 for(
size_t ipar=0;ipar<n_params;ipar++) {
722 current[sindex][ipar]=initial_points[ip_index][ipar];
726 func_ret[it]=func[it](n_params,current[sindex],
727 w_current[sindex],data_arr[sindex]);
729 if (func_ret[it]==mcmc_done) {
730 mcmc_done_flag[it]=
true;
735 if (func_ret[it]>=0 && ret_value_counts.size()>it &&
736 func_ret[it]<((int)ret_value_counts[it].size())) {
737 ret_value_counts[it][func_ret[it]]++;
739 if (meas_for_initial) {
740 meas_ret[it]=meas[it](current[sindex],w_current[sindex],
741 curr_walker[it],func_ret[it],
742 true,data_arr[sindex]);
746 if (meas_ret[it]==mcmc_done) {
747 mcmc_done_flag[it]=
true;
756 while (!done && !mcmc_done_flag[it]) {
759 for(
size_t ipar=0;ipar<n_params;ipar++) {
761 current[sindex][ipar]=
762 initial_points[ip_index][ipar]+
763 (rg[it].random()*2.0-1.0)*
764 (high[ipar]-low[ipar])*ai_initial_step;
765 }
while (current[sindex][ipar]>high[ipar] ||
766 current[sindex][ipar]<low[ipar]);
770 func_ret[it]=func[it](n_params,current[sindex],
771 w_current[sindex],data_arr[sindex]);
778 if (func_ret[it]==mcmc_done) {
779 mcmc_done_flag[it]=
true;
784 if (func_ret[it]>=0 && ret_value_counts.size()>it &&
785 func_ret[it]<((int)ret_value_counts[it].size())) {
786 ret_value_counts[it][func_ret[it]]++;
788 if (meas_ret[it]!=mcmc_done) {
789 if (meas_for_initial) {
790 meas_ret[it]=meas[it](current[sindex],w_current[sindex],
791 curr_walker[it],func_ret[it],
true,
797 mcmc_done_flag[it]=
true;
800 }
else if (init_iters>max_bad_steps) {
801 std::string err=((std::string)
"In loop with thread ")+
804 " in mcmc_para_base::mcmc().";
815 bool stop_early=
false;
817 if (mcmc_done_flag[it]==
true) {
819 scr_out <<
"mcmc (" << it <<
"," << mpi_rank
820 <<
"): Returned mcmc_done " 821 <<
"(initial; ai)." << std::endl;
835 for(curr_walker[it]=0;curr_walker[it]<
n_walk;
837 size_t sindex=n_walk*it+curr_walker[it];
838 if (w_current[sindex]>w_current[0]) {
840 w_best=w_current[sindex];
844 best=current[best_index];
850 for(curr_walker[it]=0;curr_walker[it]<
n_walk;curr_walker[it]++) {
851 size_t sindex=n_walk*it+curr_walker[it];
852 scr_out.precision(4);
853 scr_out <<
"mcmc (" << it <<
"," << mpi_rank <<
"): i_walk: ";
854 scr_out.width((
int)(1.0+log10((
double)(n_walk-1))));
855 scr_out << curr_walker[it] <<
" log_wgt: " 856 << w_current[sindex] <<
" (initial; ai)" << std::endl;
857 scr_out.precision(6);
870 #pragma omp parallel default(shared) 884 size_t ip_size=initial_points.size();
885 for(
size_t ipar=0;ipar<n_params;ipar++) {
886 current[it][ipar]=initial_points[it % ip_size][ipar];
892 func_ret[it]=func[it](n_params,current[it],w_current[it],
896 func_ret[it]=func_ret[it % ip_size];
897 w_current[it]=w_current[it % ip_size];
900 for(
size_t j=0;j<data_arr.size();j++) {
901 data_arr[it]=data_arr[it % ip_size];
912 if (func_ret[it]==mcmc_done) {
914 scr_out <<
"mcmc (" << it
915 <<
"): Initial point returned mcmc_done." 923 O2SCL_ERR((((std::string)
"Initial weight from thread ")+
925 " vanished in mcmc_para_base::mcmc().").c_str(),
936 #pragma omp parallel default(shared) 944 size_t ip_size=initial_points.size();
948 func_ret[it]=func_ret[it % ip_size];
949 current[it]=current[it % ip_size];
950 w_current[it]=w_current[it % ip_size];
951 data_arr[it]=data_arr[it % ip_size];
955 if (func_ret[it]>=0 && ret_value_counts.size()>it &&
956 func_ret[it]<((int)ret_value_counts[it].size())) {
957 ret_value_counts[it][func_ret[it]]++;
959 if (meas_for_initial) {
961 meas_ret[it]=meas[it](current[it],w_current[it],0,
962 func_ret[it],
true,data_arr[it]);
966 if (meas_ret[it]==mcmc_done) {
967 mcmc_done_flag[it]=
true;
978 bool stop_early=
false;
980 if (mcmc_done_flag[it]==
true) {
982 scr_out <<
"mcmc (" << it <<
"," << mpi_rank
983 <<
"): Returned mcmc_done " 984 <<
"(initial)." << std::endl;
999 if (w_current[it]<w_best) {
1001 w_best=w_current[it];
1007 scr_out.precision(4);
1008 scr_out <<
"mcmc: ";
1010 scr_out << w_current[it] <<
" ";
1012 scr_out <<
" (initial)" << std::endl;
1013 scr_out.precision(6);
1020 meas_for_initial=
true;
1027 std::cout <<
"Press a key and type enter to continue. ";
1038 if (aff_inv==
false) {
1044 #pragma omp parallel default(shared) 1052 bool main_done=
false;
1053 size_t mcmc_iters=0;
1055 while (!main_done) {
1063 q_prop[it]=prop_dist[it]->log_metrop_hast(current[it],next[it]);
1065 if (!std::isfinite(q_prop[it])) {
1066 O2SCL_ERR2(
"Proposal distribution not finite in ",
1073 for(
size_t k=0;k<n_params;k++) {
1074 next[it][k]=current[it][k]+(rg[it].random()*2.0-1.0)*
1075 (high[k]-low[k])/step_fac;
1087 for(
size_t k=0;k<n_params;k++) {
1088 if (next[it][k]<low[k] || next[it][k]>high[k]) {
1091 if (next[it][k]<low[k]) {
1092 std::cout <<
"mcmc (" << it <<
"," 1093 << mpi_rank <<
"): Parameter with index " << k
1094 <<
" and value " << next[it][k]
1095 <<
" smaller than limit " << low[k] << std::endl;
1096 scr_out <<
"mcmc (" << it <<
"," 1097 << mpi_rank <<
"): Parameter with index " << k
1098 <<
" and value " << next[it][k]
1099 <<
" smaller than limit " << low[k] << std::endl;
1101 std::cout <<
"mcmc (" << it <<
"," << mpi_rank
1102 <<
"): Parameter with index " << k
1103 <<
" and value " << next[it][k]
1104 <<
" larger than limit " << high[k] << std::endl;
1105 scr_out <<
"mcmc (" << it <<
"," << mpi_rank
1106 <<
"): Parameter with index " << k
1107 <<
" and value " << next[it][k]
1108 <<
" larger than limit " << high[k] << std::endl;
1116 if (func_ret[it]!=mcmc_skip) {
1117 if (switch_arr[n_walk*it+curr_walker[it]]==
false) {
1118 func_ret[it]=func[it](n_params,next[it],w_next[it],
1119 data_arr[it*n_walk+curr_walker[it]+
1122 func_ret[it]=func[it](n_params,next[it],w_next[it],
1123 data_arr[it*n_walk+curr_walker[it]]);
1125 if (func_ret[it]==mcmc_done) {
1126 mcmc_done_flag[it]=
true;
1128 if (func_ret[it]>=0 && ret_value_counts.size()>it &&
1129 func_ret[it]<((int)ret_value_counts[it].size())) {
1130 ret_value_counts[it][func_ret[it]]++;
1140 size_t sindex=n_walk*it+curr_walker[it];
1143 if (always_accept && func_ret[it]==
success) accept=
true;
1146 double r=rg[it].random();
1149 if (mcmc_iters%100==0) {
1150 std::cout.setf(std::ios::showpos);
1151 std::cout.precision(4);
1152 double v1=prop_dist[it]->log_pdf(next[it],current[it]);
1153 double v2=prop_dist[it]->log_pdf(current[it],next[it]);
1154 std::cout <<
"PD: " << w_current[it] <<
" " 1155 << w_next[it] <<
" " 1156 << v1 <<
" " << v2 <<
" " << q_prop[it] <<
" " 1157 << w_next[it]-w_current[sindex]+q_prop[it]
1161 std::cout.precision(6);
1162 std::cout.unsetf(std::ios::showpos);
1164 if (r<exp(w_next[it]-w_current[sindex]+q_prop[it])) {
1170 if (r<exp(w_next[it]-w_current[sindex])) {
1184 if (switch_arr[sindex]==
false) {
1185 meas_ret[it]=meas[it](next[it],w_next[it],
1186 curr_walker[it],func_ret[it],
true,
1187 data_arr[sindex+n_threads*
n_walk]);
1189 meas_ret[it]=meas[it](next[it],w_next[it],
1190 curr_walker[it],func_ret[it],
true,
1196 current[sindex]=next[it];
1197 w_current[sindex]=w_next[it];
1198 switch_arr[sindex]=!(switch_arr[sindex]);
1208 if (switch_arr[sindex]==
false) {
1209 meas_ret[it]=meas[it](next[it],w_next[it],
1210 curr_walker[it],func_ret[it],
false,
1211 data_arr[sindex+n_threads*
n_walk]);
1213 meas_ret[it]=meas[it](next[it],w_next[it],
1214 curr_walker[it],func_ret[it],
false,
1228 #pragma omp critical (o2scl_mcmc_para_best_point) 1233 if (switch_arr[n_walk*it+curr_walker[it]]==
false) {
1234 best_point(best,w_best,data_arr[curr_walker[it]+n_walk*it+
1237 best_point(best,w_best,data_arr[curr_walker[it]+n_walk*it]);
1244 if (meas_ret[it]==mcmc_done || func_ret[it]==mcmc_done) {
1249 O2SCL_ERR((((std::string)
"Measurement function returned ")+
1251 " in mcmc_para_base::mcmc().").c_str(),
1259 if (main_done==
false) {
1263 if (warm_up && mcmc_iters==n_warm_up) {
1269 scr_out <<
"o2scl::mcmc_para: Thread " << it
1270 <<
" finished warmup." << std::endl;
1277 if (main_done==
false && warm_up==
false && max_iters>0 &&
1278 mcmc_iters==max_iters) {
1280 scr_out <<
"o2scl::mcmc_para: Thread " << it
1281 <<
" stopping because number of iterations (" 1282 << mcmc_iters <<
") equal to max_iters (" << max_iters
1283 <<
")." << std::endl;
1288 if (main_done==
false) {
1295 if (max_time>0.0 && elapsed>max_time) {
1297 scr_out <<
"o2scl::mcmc_para: Thread " << it
1298 <<
" stopping because elapsed (" << elapsed
1299 <<
") > max_time (" << max_time <<
")." 1320 bool main_done=
false;
1321 size_t mcmc_iters=0;
1323 while (!main_done) {
1326 std::vector<double> smove_z(n_threads);
1332 curr_walker[it]=mcmc_iters %
n_walk;
1337 #pragma omp parallel default(shared) 1352 ij=((size_t)(rg[it].random()*((double)n_walk)));
1353 }
while (ij==curr_walker[it] || ij>=n_walk);
1356 double p=rg[it].random();
1358 smove_z[it]=(1.0-2.0*p+2.0*a*p+p*p-2.0*a*p*p+a*a*p*p)/a;
1361 for(
size_t i=0;i<n_params;i++) {
1362 next[it][i]=current[n_walk*it+ij][i]+
1363 smove_z[it]*(current[n_walk*it+curr_walker[it]][i]-
1364 current[n_walk*it+ij][i]);
1375 for(
size_t k=0;k<n_params;k++) {
1376 if (next[it][k]<low[k] || next[it][k]>high[k]) {
1379 if (next[it][k]<low[k]) {
1380 std::cout <<
"mcmc (" << it <<
"," 1381 << mpi_rank <<
"): Parameter with index " << k
1382 <<
" and value " << next[it][k]
1383 <<
" smaller than limit " << low[k] << std::endl;
1384 scr_out <<
"mcmc (" << it <<
"," 1385 << mpi_rank <<
"): Parameter with index " << k
1386 <<
" and value " << next[it][k]
1387 <<
" smaller than limit " << low[k] << std::endl;
1389 std::cout <<
"mcmc (" << it <<
"," << mpi_rank
1390 <<
"): Parameter with index " << k
1391 <<
" and value " << next[it][k]
1392 <<
" larger than limit " << high[k] << std::endl;
1393 scr_out <<
"mcmc (" << it <<
"," << mpi_rank
1394 <<
"): Parameter with index " << k
1395 <<
" and value " << next[it][k]
1396 <<
" larger than limit " << high[k] << std::endl;
1404 if (func_ret[it]!=mcmc_skip) {
1405 if (switch_arr[n_walk*it+curr_walker[it]]==
false) {
1406 func_ret[it]=func[it](n_params,next[it],w_next[it],
1407 data_arr[it*n_walk+curr_walker[it]+
1410 func_ret[it]=func[it](n_params,next[it],w_next[it],
1411 data_arr[it*n_walk+curr_walker[it]]);
1413 if (func_ret[it]==mcmc_done) {
1414 mcmc_done_flag[it]=
true;
1416 if (func_ret[it]>=0 && ret_value_counts.size()>it &&
1417 func_ret[it]<((int)ret_value_counts[it].size())) {
1418 ret_value_counts[it][func_ret[it]]++;
1436 scr_out <<
"it: " << it <<
" q_prop[it]: " 1437 << q_prop[it] << std::endl;
1439 if (func_ret[it]==mcmc_done) {
1440 scr_out <<
"mcmc (" << it <<
"," << mpi_rank
1441 <<
"): Returned mcmc_done." 1443 }
else if (func_ret[it]==mcmc_skip && verbose>=3) {
1444 scr_out <<
"mcmc (" << it
1445 <<
"): Parameter(s) out of range: " << std::endl;
1446 scr_out.setf(std::ios::showpos);
1447 for(
size_t k=0;k<n_params;k++) {
1448 scr_out << k <<
" " << low[k] <<
" " 1449 << next[it][k] <<
" " << high[k];
1450 if (next[it][k]<low[k] || next[it][k]>high[k]) {
1453 scr_out << std::endl;
1455 scr_out.unsetf(std::ios::showpos);
1457 func_ret[it]!=mcmc_skip) {
1459 scr_out <<
"mcmc (" << it <<
"," << mpi_rank
1460 <<
"): Function returned failure " 1461 << func_ret[it] <<
" at point ";
1462 for(
size_t k=0;k<n_params;k++) {
1463 scr_out << next[it][k] <<
" ";
1465 scr_out << std::endl;
1476 #pragma omp parallel default(shared) 1485 size_t sindex=n_walk*it+curr_walker[it];
1491 if (always_accept && func_ret[it]==
success) accept=
true;
1494 double r=rg[it].random();
1497 double ai_ratio=pow(smove_z[it],((
double)n_params)-1.0)*
1498 exp(w_next[it]-w_current[sindex]);
1502 }
else if (pd_mode) {
1503 if (r<exp(w_next[it]-w_current[sindex]+q_prop[it])) {
1508 if (r<exp(w_next[it]-w_current[sindex])) {
1522 if (switch_arr[sindex]==
false) {
1523 meas_ret[it]=meas[it](next[it],w_next[it],
1524 curr_walker[it],func_ret[it],
true,
1525 data_arr[sindex+n_threads*
n_walk]);
1527 meas_ret[it]=meas[it](next[it],w_next[it],
1528 curr_walker[it],func_ret[it],
true,
1534 current[sindex]=next[it];
1535 w_current[sindex]=w_next[it];
1536 switch_arr[sindex]=!(switch_arr[sindex]);
1545 if (switch_arr[sindex]==
false) {
1546 meas_ret[it]=meas[it](next[it],w_next[it],
1547 curr_walker[it],func_ret[it],
false,
1548 data_arr[sindex+n_threads*
n_walk]);
1550 meas_ret[it]=meas[it](next[it],w_next[it],
1551 curr_walker[it],func_ret[it],
false,
1568 size_t sindex=n_walk*it+curr_walker[it];
1569 scr_out.precision(4);
1570 scr_out <<
"mcmc (" << it <<
"," << mpi_rank <<
"): iter: ";
1571 scr_out.width((
int)(1.0+log10((
double)(n_params-1))));
1572 scr_out << mcmc_iters <<
" i_walk: " 1573 << curr_walker[it] <<
" log_wgt: " 1574 << w_current[sindex] << std::endl;
1575 scr_out.precision(6);
1584 if (switch_arr[n_walk*it+curr_walker[it]]==
false) {
1585 best_point(best,w_best,data_arr[curr_walker[it]+n_walk*it+
1588 best_point(best,w_best,data_arr[curr_walker[it]+n_walk*it]);
1596 if (meas_ret[it]==mcmc_done || func_ret[it]==mcmc_done) {
1601 O2SCL_ERR((((std::string)
"Measurement function returned ")+
1603 " in mcmc_para_base::mcmc().").c_str(),
1612 if (main_done==
false) {
1616 if (warm_up && mcmc_iters==n_warm_up) {
1624 scr_out <<
"mcmc: Finished warmup." << std::endl;
1631 std::cout <<
"Press a key and type enter to continue. ";
1637 if (main_done==
false && warm_up==
false && max_iters>0 &&
1638 mcmc_iters==max_iters) {
1640 scr_out <<
"mcmc: Stopping because number of iterations (" 1641 << mcmc_iters <<
") equal to max_iters (" << max_iters
1642 <<
")." << std::endl;
1647 if (main_done==
false) {
1654 if (max_time>0.0 && elapsed>max_time) {
1656 scr_out <<
"mcmc: Stopping because elapsed (" << elapsed
1657 <<
") > max_time (" << max_time <<
")." 1684 virtual int mcmc(
size_t n_params, vec_t &low, vec_t &high,
1685 func_t &func, measure_t &meas) {
1688 omp_set_num_threads(n_threads);
1692 std::vector<func_t> vf(n_threads);
1693 std::vector<measure_t> vm(n_threads);
1698 return mcmc(n_params,low,high,vf,vm);
1710 prop_dist.resize(pv.size());
1711 for(
size_t i=0;i<pv.size();i++) {
1712 prop_dist[i]=&pv[i];
1726 prop_dist.resize(pv.size());
1727 for(
size_t i=0;i<pv.size();i++) {
1791 template<
class func_t,
class fill_t,
class data_t,
class vec_t=ubvector>
1793 std::function<int(const vec_t &,double,size_t,int,bool,data_t &)>,
1799 typedef std::function<int(const vec_t &,double,size_t,int,bool,data_t &)>
1815 std::shared_ptr<o2scl::table_units<> >
table;
1834 table=std::shared_ptr<o2scl::table_units<> >
1839 if (table_prealloc>0) {
1840 table->set_maxlines(table_prealloc);
1842 table->new_column(
"rank");
1843 table->new_column(
"thread");
1844 table->new_column(
"walker");
1845 table->new_column(
"mult");
1846 table->new_column(
"log_wgt");
1847 for(
size_t i=0;i<col_names.size();i++) {
1848 table->new_column(col_names[i]);
1849 if (col_units[i].length()>0) {
1850 table->set_unit(col_names[i],col_units[i]);
1856 walker_accept_rows[i]=-1;
1858 walker_reject_rows.resize(this->
n_walk*this->n_threads);
1860 walker_reject_rows[i]=-1;
1864 std::cout <<
"mcmc: Table column names and units: " << std::endl;
1865 for(
size_t i=0;i<table->get_ncolumns();i++) {
1866 std::cout << table->get_column_name(i) <<
" " 1867 << table->get_unit(table->get_column_name(i)) << std::endl;
1874 O2SCL_ERR2(
"Flag 'prev_read' is true but table pointer is 0 ",
1881 if (table->get_ncolumns()!=5+col_names.size()) {
1882 std::string str=((std::string)
"Table does not have correct ")+
1883 "number of columns in mcmc_para_table::mcmc_init()."+
1888 if (!table->is_column(
"rank") ||
1889 !table->is_column(
"thread") ||
1890 !table->is_column(
"walker") ||
1891 !table->is_column(
"mult") ||
1892 !table->is_column(
"log_wgt")) {
1893 O2SCL_ERR2(
"Table does not have the correct internal columns ",
1897 O2SCL_ERR2(
"Array walker_accept_rows does not have correct size ",
1901 O2SCL_ERR2(
"Array walker_reject_rows does not have correct size ",
1913 last_write_time=MPI_Wtime();
1915 last_write_time=time(0);
1918 return parent_t::mcmc_init();
1924 std::vector<double> &line, data_t &dat,
1925 size_t i_walker, fill_t &fill) {
1928 size_t i_thread=omp_get_thread_num();
1936 line.push_back(i_thread);
1938 line.push_back(i_walker);
1940 line.push_back(1.0);
1941 line.push_back(log_weight);
1942 for(
size_t i=0;i<pars.size();i++) {
1943 line.push_back(pars[i]);
1945 int tempi=fill(pars,log_weight,line,dat);
2026 this->
scr_out <<
"mcmc: Start write_files(). mpi_rank: " 2028 << this->
mpi_size <<
" table_io_chunk: " 2029 << table_io_chunk << std::endl;
2032 std::vector<o2scl::table_units<> > tab_arr;
2033 bool rank_sent=
false;
2036 if (table_io_chunk>1) {
2037 if (this->
mpi_rank%table_io_chunk==0) {
2039 for(
int i=0;i<table_io_chunk-1;i++) {
2043 tab_arr.push_back(t);
2044 o2scl_table_mpi_recv(child,tab_arr[tab_arr.size()-1]);
2050 o2scl_table_mpi_send(*table,parent);
2059 int tag=0, buffer=0;
2060 if (sync_write && this->
mpi_size>1 &&
2062 MPI_Recv(&buffer,1,MPI_INT,this->
mpi_rank-table_io_chunk,
2063 tag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
2071 if (first_write==
false) {
2080 hf.
set_szt(
"file_update_iters",this->file_update_iters);
2081 hf.
set_szt(
"file_update_time",this->file_update_time);
2085 hf.
set_szt(
"n_params",this->n_params);
2093 hf.
seti(
"store_rejects",this->store_rejects);
2094 hf.
seti(
"table_sequence",this->table_sequence);
2112 hf.
seti(
"n_tables",tab_arr.size()+1);
2113 if (rank_sent==
false) {
2114 hdf_output(hf,*table,
"markov_chain_0");
2116 for(
size_t i=0;i<tab_arr.size();i++) {
2117 std::string name=((std::string)
"markov_chain_")+
szttos(i+1);
2118 hdf_output(hf,tab_arr[i],name);
2124 if (sync_write && this->
mpi_size>1 &&
2126 MPI_Send(&buffer,1,MPI_INT,this->
mpi_rank+table_io_chunk,
2127 tag,MPI_COMM_WORLD);
2132 this->
scr_out <<
"mcmc: Done write_files()." << std::endl;
2140 file_update_iters=0;
2141 file_update_time=0.0;
2143 store_rejects=
false;
2144 table_sequence=
true;
2154 std::vector<std::string> units) {
2155 if (names.size()!=units.size()) {
2156 O2SCL_ERR2(
"Size of names and units arrays don't match in ",
2183 int tag=0, buffer=0;
2185 MPI_Recv(&buffer,1,MPI_INT,this->
mpi_rank-1,
2186 tag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
2198 MPI_Send(&buffer,1,MPI_INT,this->
mpi_rank+1,
2199 tag,MPI_COMM_WORLD);
2207 std::cout <<
"Initial points: Finding last " << n_points
2208 <<
" points from file named " 2209 << fname <<
" ." << std::endl;
2214 for(
size_t it=0;it<this->
n_threads;it++) {
2215 for(
size_t iw=0;iw<this->
n_walk;iw++) {
2218 size_t windex=it*this->n_walk+iw;
2221 for(
int row=tip.
get_nlines()-1;row>=0 && found==
false;row--) {
2222 if (tip.
get(
"walker",row)==iw &&
2223 tip.
get(
"thread",row)==it &&
2224 tip.
get(
"mult",row)>0.5) {
2228 std::cout <<
"Function initial_point_file_last():\n\tit: " 2229 << it <<
" rank: " << this->
mpi_rank 2230 <<
" iw: " << iw <<
" row: " 2231 << row <<
" log_wgt: " << tip.
get(
"log_wgt",row)
2236 for(
size_t ip=0;ip<n_param_loc;ip++) {
2242 O2SCL_ERR(
"Function initial_points_file_last() failed.",
2269 int tag=0, buffer=0;
2271 MPI_Recv(&buffer,1,MPI_INT,this->
mpi_rank-1,
2272 tag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
2284 MPI_Send(&buffer,1,MPI_INT,this->
mpi_rank+1,
2285 tag,MPI_COMM_WORLD);
2293 std::cout <<
"Initial points: Finding last " << n_points
2294 <<
" points from file named " 2295 << fname <<
" ." << std::endl;
2301 size_t decrement=nlines/n_points;
2302 if (decrement<1) decrement=1;
2305 for(
size_t k=0;k<n_points;k++) {
2309 std::cout <<
"Function initial_point_file_dist():\n\trow: " 2310 << row <<
" log_wgt: " << tip.
get(
"log_wgt",row)
2315 for(
size_t ip=0;ip<n_param_loc;ip++) {
2335 double thresh=1.0e-6,
2343 int tag=0, buffer=0;
2345 MPI_Recv(&buffer,1,MPI_INT,this->
mpi_rank-1,
2346 tag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
2358 MPI_Send(&buffer,1,MPI_INT,this->
mpi_rank+1,
2359 tag,MPI_COMM_WORLD);
2367 std::cout <<
"Initial points: Finding best " << n_points
2368 <<
" unique points from file named " 2369 << fname <<
" ." << std::endl;
2372 typedef std::map<double,int,std::greater<double> > map_t;
2377 m.insert(std::make_pair(tip.
get(
"log_wgt",k),k));
2387 for(map_t::iterator mit=m.begin();mit!=m.end();mit++) {
2388 map_t::iterator mit2=mit;
2390 if (mit2!=m.end()) {
2391 if (fabs(mit->first-mit2->first)<thresh) {
2393 std::cout <<
"Removing duplicate weights: " 2394 << mit->first <<
" " << mit2->first << std::endl;
2403 }
while (found==
true);
2406 if (m.size()<n_points) {
2407 O2SCL_ERR2(
"Could not find enough points in file in ",
2408 "mcmc_para::initial_points_file_best().",
2414 map_t::iterator mit=m.begin();
2415 for(
size_t k=0;k<n_points;k++) {
2416 int row=mit->second;
2418 std::cout <<
"Initial point " << k <<
" at row " 2419 << row <<
" has log_weight= " 2420 << tip.
get(
"log_wgt",row) << std::endl;
2423 for(
size_t ip=0;ip<n_param_loc;ip++) {
2439 virtual int mcmc(
size_t n_params_local,
2440 vec_t &low, vec_t &high, std::vector<func_t> &func,
2441 std::vector<fill_t> &fill) {
2443 n_params=n_params_local;
2459 std::vector<internal_measure_t> meas(this->
n_threads);
2460 for(
size_t it=0;it<this->
n_threads;it++) {
2462 (std::mem_fn<
int(
const vec_t &,
double,
size_t,
int,
bool,
2463 data_t &,
size_t, fill_t &)>
2465 std::placeholders::_2,std::placeholders::_3,
2466 std::placeholders::_4,std::placeholders::_5,
2467 std::placeholders::_6,it,std::ref(fill[it]));
2470 return parent_t::mcmc(n_params,low,high,func,meas);
2495 chain_sizes.resize(ntot);
2497 for(
size_t it=0;it<this->
n_threads;it++) {
2498 for(
size_t iw=0;iw<this->
n_walk;iw++) {
2499 size_t ix=it*this->n_walk+iw;
2502 for(
size_t j=istart;j<table->get_nlines();j+=ntot) {
2503 if (table->get(
"mult",j)>0.5) chain_sizes[ix]++;
2523 std::string name=
"") {
2531 if (!table->is_column(
"rank") ||
2532 !table->is_column(
"thread") ||
2533 !table->is_column(
"walker") ||
2534 !table->is_column(
"mult") ||
2535 !table->is_column(
"log_wgt")) {
2536 O2SCL_ERR2(
"Table does not have the correct internal columns ",
2537 "in mcmc_para_table::read_prev_results().",
2549 walker_accept_rows.resize(ntot);
2550 walker_reject_rows.resize(ntot);
2552 for(
size_t j=0;j<ntot;j++) {
2553 walker_accept_rows[j]=-1;
2554 walker_reject_rows[j]=-1;
2557 for(
size_t j=0;j<table->get_nlines();j++) {
2559 size_t i_thread=((size_t)(table->get(
"thread",j)+1.0e-12));
2560 size_t i_walker=((size_t)(table->get(
"walker",j)+1.0e-12));
2563 size_t windex=i_thread*this->n_walk+i_walker;
2565 if (table->get(
"mult",j)>0.5) {
2566 walker_accept_rows[windex]=j;
2567 }
else if (table->get(
"mult",j)<-0.5) {
2568 walker_reject_rows[windex]=j;
2576 for(
size_t j=0;j<ntot;j++) {
2577 if (walker_accept_rows[j]<0) found=
false;
2584 for(
size_t j=0;j<ntot;j++) {
2586 for(
size_t k=0;k<n_param_loc;k++) {
2587 this->
initial_points[j][k]=table->get(k+5,walker_accept_rows[j]);
2591 std::cout <<
"Previous table was read, but initial points not set." 2596 std::cout <<
"mcmc_para_table::read_prev_results():" << std::endl;
2597 std::cout <<
" index walker_accept_rows walker_reject_rows" 2599 for(
size_t j=0;j<ntot;j++) {
2602 std::cout << j <<
" ";
2604 std::cout << walker_accept_rows[j] <<
" ";
2606 std::cout << walker_reject_rows[j] << std::endl;
2625 if (file_update_iters>0) {
2626 size_t total_iters=0;
2630 if (total_iters>=last_write_iters+file_update_iters) {
2632 this->
scr_out <<
"mcmc: Writing to file. total_iters: " 2633 << total_iters <<
" file_update_iters: " 2634 << file_update_iters <<
" last_write_iters: " 2635 << last_write_iters << std::endl;
2638 last_write_iters=total_iters;
2642 if (updated==
false && file_update_time>0.0) {
2644 double elapsed=MPI_Wtime()-last_write_time;
2646 double elapsed=time(0)-last_write_time;
2648 if (elapsed>file_update_time) {
2650 this->
scr_out <<
"mcmc: Writing to file. elapsed: " 2651 << elapsed <<
" file_update_time: " 2652 << file_update_time <<
" last_write_time: " 2653 << last_write_time << std::endl;
2657 last_write_time=MPI_Wtime();
2659 last_write_time=time(0);
2671 virtual int add_line(
const vec_t &pars,
double log_weight,
2672 size_t walker_ix,
int func_ret,
2673 bool mcmc_accept, data_t &dat,
2674 size_t i_thread, fill_t &fill) {
2677 size_t windex=i_thread*this->
n_walk+walker_ix;
2686 if ((mcmc_accept || store_rejects) && walker_accept_rows[windex]<0) {
2689 if (table_sequence) {
2690 if (walker_accept_rows[windex]>walker_reject_rows[windex]) {
2691 next_row=walker_accept_rows[windex]+ntot;
2693 next_row=walker_reject_rows[windex]+ntot;
2696 if (walker_accept_rows[windex]>walker_reject_rows[windex]) {
2697 next_row=walker_accept_rows[windex]+1;
2699 next_row=walker_reject_rows[windex]+1;
2705 #pragma omp critical (o2scl_mcmc_para_table_add_line) 2709 while (next_row<((
int)table->get_nlines()) &&
2710 fabs(table->get(
"mult",next_row))>0.1) {
2719 if (next_row>=((
int)table->get_nlines())) {
2720 size_t istart=table->get_nlines();
2722 table->set_nlines(table->get_nlines()+ntot);
2725 for(
size_t i=0;i<this->
n_walk;i++) {
2726 table->set(
"rank",istart+j*this->n_walk+i,this->
mpi_rank);
2727 table->set(
"thread",istart+j*this->n_walk+i,j);
2728 table->set(
"walker",istart+j*this->n_walk+i,i);
2729 table->set(
"mult",istart+j*this->n_walk+i,0.0);
2730 table->set(
"log_wgt",istart+j*this->n_walk+i,0.0);
2736 if (func_ret==0 && (mcmc_accept || store_rejects)) {
2738 if (next_row>=((
int)(table->get_nlines()))) {
2742 std::vector<double> line;
2743 int fret=fill_line(pars,log_weight,line,dat,walker_ix,fill);
2747 if (store_rejects && mcmc_accept==
false) {
2763 if (line.size()!=table->get_ncolumns()) {
2764 std::cout <<
"line: " << line.size() <<
" columns: " 2765 << table->get_ncolumns() << std::endl;
2766 for(
size_t k=0;k<table->get_ncolumns() || k<line.size();k++) {
2767 std::cout << k <<
". ";
2768 if (k<table->get_ncolumns()) {
2769 std::cout << table->get_column_name(k) <<
" ";
2770 std::cout << table->get_unit(table->get_column_name(k)) <<
" ";
2772 if (k<line.size()) std::cout << line[k] <<
" ";
2773 std::cout << std::endl;
2775 O2SCL_ERR(
"Table misalignment in mcmc_para_table::add_line().",
2780 table->set_row(((
size_t)next_row),line);
2784 this->
scr_out <<
"mcmc: Setting data at row " << next_row
2786 for(
size_t k=0;k<line.size();k++) {
2788 this->
scr_out << table->get_column_name(k) <<
" ";
2789 this->
scr_out << table->get_unit(table->get_column_name(k));
2790 this->
scr_out <<
" " << line[k] << std::endl;
2799 critical_extra(i_thread);
2806 if (walker_accept_rows[windex]<0 ||
2807 walker_accept_rows[windex]>=((
int)table->get_nlines())) {
2808 O2SCL_ERR2(
"Invalid row for incrementing multiplier in ",
2811 double mult_old=table->get(
"mult",walker_accept_rows[windex]);
2816 table->set(
"mult",walker_accept_rows[windex],mult_old+1.0);
2818 this->
scr_out <<
"mcmc: Updating mult of row " 2819 << walker_accept_rows[windex]
2820 <<
" from " << mult_old <<
" to " 2821 << mult_old+1.0 << std::endl;
2829 walker_accept_rows[windex]=next_row;
2830 }
else if (store_rejects && func_ret==0) {
2831 walker_reject_rows[windex]=next_row;
2848 for(i=table->get_nlines()-1;i>=0 && done==
false;i--) {
2850 if (fabs(table->get(
"mult",i))<0.1) {
2854 if (i+2<((
int)table->get_nlines())) {
2855 table->set_nlines(i+2);
2860 return parent_t::mcmc_cleanup();
2866 std::vector<size_t> chain_sizes;
2867 get_chain_sizes(chain_sizes);
2868 size_t min_size=chain_sizes[0];
2869 for(
size_t i=1;i<chain_sizes.size();i++) {
2870 if (chain_sizes[i]<min_size) min_size=chain_sizes[i];
2872 size_t N_max=min_size/2;
2873 ac_coeffs.resize(ncols,N_max-1);
2874 for(
size_t i=0;i<ncols;i++) {
2875 for(
size_t ell=1;ell<N_max;ell++) {
2876 ac_coeffs(i,ell-1)=0.0;
2881 size_t cstart=table->lookup_column(
"log_wgt")+1;
2882 for(
size_t i=0;i<ncols;i++) {
2884 for(
size_t k=0;k<this->
n_walk;k++) {
2885 size_t tindex=j*this->n_walk+k;
2886 for(
size_t ell=1;ell<N_max;ell++) {
2887 const double &x=(*table)[cstart+i][table_row];
2888 double mean=o2scl::vector_mean<const double *>
2889 (chain_sizes[tindex]+1,&x);
2891 <
const double *>(chain_sizes[tindex]+1,&x,ell,mean);
2893 table_row+=chain_sizes[tindex]+1;
2896 for(
size_t ell=1;ell<N_max;ell++) {
2897 ac_coeffs(i,ell-1)/=((double)n_tot);
2906 ubvector &ac_lengths) {
2907 size_t N_max=ac_coeffs_cols.size2();
2908 ac_lengths.resize(ncols);
2909 for(
size_t icol=0;icol<ncols;icol++) {
2910 std::vector<double> tau(N_max);
2911 for(
size_t i=5;i<N_max;i++) {
2913 for(
size_t j=0;j<i;j++) {
2914 sum+=ac_coeffs_cols(icol,j);
2917 std::cout << tau[i] <<
" " << ((double)i)/5.0 << std::endl;
2919 std::cout << std::endl;
2929 std::shared_ptr<o2scl::table_units<> > table2=
2933 for(
size_t j=0;j<this->
n_walk;j++) {
2934 std::string func=std::string(
"abs(walker-")+
o2scl::szttos(j)+
2936 table->copy_rows(func,*table2);
2958 for(
size_t it=0;it<this->
n_threads;it++) {
2960 size_t n=table->get_nlines();
2962 O2SCL_ERR2(
"Cannot reblock. Not enough data in ",
2965 size_t n_block=n/n_blocks;
2966 size_t m=table->get_ncolumns();
2967 for(
size_t j=0;j<n_blocks;j++) {
2970 for(
size_t i=0;i<m;i++) {
2973 for(
size_t k=j*n_block;k<(j+1)*n_block;k++) {
2974 mult+=(*table)[
"mult"][k];
2975 for(
size_t i=1;i<m;i++) {
2976 dat[i]+=(*table)[i][k]*(*table)[
"mult"][k];
2979 table->set(
"mult",j,mult);
2980 for(
size_t i=1;i<m;i++) {
2982 table->set(i,j,dat[i]);
2985 table->set_nlines(n_blocks);
3001 template<
class func_t,
class fill_t,
class data_t,
class vec_t=ubvector>
3034 hf.
sets_vec(
"cl_args",this->cl_args);
3082 p_file_update_iters.
s=&this->file_update_iters;
3083 p_file_update_iters.
help=((std::string)
"Number of MCMC successes ")+
3084 "between file updates (default 0 for no file updates).";
3085 cl.
par_list.insert(std::make_pair(
"file_update_iters",
3086 &p_file_update_iters));
3088 p_file_update_time.
d=&this->file_update_time;
3089 p_file_update_time.
help=((std::string)
"Time ")+
3090 "between file updates (default 0.0 for no file updates).";
3091 cl.
par_list.insert(std::make_pair(
"file_update_time",
3092 &p_file_update_time));
3103 p_max_time.
help=((std::string)
"Maximum run time in seconds ")+
3104 "(default 86400 sec or 1 day).";
3105 cl.
par_list.insert(std::make_pair(
"max_time",&p_max_time));
3108 p_max_iters.
help=((std::string)
"If non-zero, limit the number of ")+
3109 "iterations to be less than the specified number (default zero).";
3110 cl.
par_list.insert(std::make_pair(
"max_iters",&p_max_iters));
3113 p_prefix.
help=
"Output file prefix (default 'mcmc\').";
3114 cl.
par_list.insert(std::make_pair(
"prefix",&p_prefix));
3124 p_step_fac.
help=((std::string)
"MCMC step factor. The step size for ")+
3125 "each variable is taken as the difference between the high and low "+
3126 "limits divided by this factor (default 10.0). This factor can "+
3127 "be increased if the acceptance rate is too small, but care must "+
3128 "be taken, e.g. if the conditional probability is multimodal. If "+
3129 "this step size is smaller than 1.0, it is reset to 1.0 .";
3130 cl.
par_list.insert(std::make_pair(
"step_fac",&p_step_fac));
3133 p_n_warm_up.
help=((std::string)
"Minimum number of warm up iterations ")+
3135 cl.
par_list.insert(std::make_pair(
"n_warm_up",&p_n_warm_up));
3138 p_verbose.
help=((std::string)
"MCMC verbosity parameter ")+
3140 cl.
par_list.insert(std::make_pair(
"mcmc_verbose",&p_verbose));
3143 p_max_bad_steps.
help=((std::string)
"Maximum number of bad steps ")+
3145 cl.
par_list.insert(std::make_pair(
"max_bad_steps",&p_max_bad_steps));
3148 p_n_walk.
help=((std::string)
"Number of walkers ")+
3150 cl.
par_list.insert(std::make_pair(
"n_walk",&p_n_walk));
3153 p_user_seed.
help=((std::string)
"Seed for multiplier for random ")+
3154 "number generator. If zero is given (the default), then mcmc() "+
3155 "uses time(0) to generate a random seed.";
3156 cl.
par_list.insert(std::make_pair(
"user_seed",&p_user_seed));
3159 p_aff_inv.
help=((std::string)
"If true, then use affine-invariant ")+
3160 "sampling (default false).";
3161 cl.
par_list.insert(std::make_pair(
"aff_inv",&p_aff_inv));
3163 p_table_sequence.
b=&this->table_sequence;
3164 p_table_sequence.
help=((std::string)
"If true, then ensure equal spacing ")+
3165 "between walkers (default true).";
3166 cl.
par_list.insert(std::make_pair(
"table_sequence",&p_table_sequence));
3168 p_store_rejects.
b=&this->store_rejects;
3169 p_store_rejects.
help=((std::string)
"If true, then store MCMC rejections ")+
3171 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.
double get(std::string scol, size_t row) const
Get value from row row of column named col. .
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")
virtual void critical_extra(size_t i_thread)
Additional code to execute inside the OpenMP critical section.
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.
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
bool prev_read
If true, previous results have been read.
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.
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 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)
bool meas_for_initial
If true, call the measurement function for the initial point.
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.
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)
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() )
void set_proposal_ptrs(prob_vec_t &pv)
Set pointers to proposal distributions.
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.
requested feature not (yet) implemented
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.
o2scl::mcmc_para_table< func_t, fill_t, data_t, vec_t > parent_t
The parent typedef.
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.
size_t table_prealloc
Desc.
Integer parameter for o2scl::cli.
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.
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.
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.