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);
223 <<
" reject=" <<
n_reject[it] << std::endl;
236 virtual void best_point(vec_t &best,
double w_best, data_t &dat) {
403 MPI_Comm_rank(MPI_COMM_WORLD,&this->mpi_rank);
404 MPI_Comm_size(MPI_COMM_WORLD,&this->mpi_size);
435 virtual int mcmc(
size_t n_params, vec_t &low, vec_t &high,
436 std::vector<func_t> &func, std::vector<measure_t> &meas) {
442 if (func.size()==0 || meas.size()==0) {
443 O2SCL_ERR2(
"Size of 'func' or 'meas' array is zero in ",
448 std::cout <<
"mcmc_para::mcmc(): Not enough functions for "
449 <<
n_threads <<
" threads. Setting n_threads to "
450 << func.size() <<
"." << std::endl;
456 std::cout <<
"mcmc_para::mcmc(): Not enough measurement objects for "
457 <<
n_threads <<
" threads. Setting n_threads to "
458 << meas.size() <<
"." << std::endl;
477 for(
size_t k=0;k<n_params;k++) {
486 for(
size_t ipar=0;ipar<n_params;ipar++) {
488 std::cout <<
"iip,ipar: " << iip <<
" " << ipar << std::endl;
492 O2SCL_ERR2(
"Initial point vector not correctly sized ",
502 ") in mcmc_base::mcmc().").c_str(),
528 std::cout.setf(std::ios::scientific);
529 std::cout << i <<
" ";
531 std::cout << j <<
" ";
547 std::cout <<
"mcmc_para::mcmc(): "
548 <<
n_threads <<
" threads were requested but the "
549 <<
"-DO2SCL_OPENMP flag was not used during "
550 <<
"compilation. Setting n_threads to 1."
562 std::cout <<
"mcmc_para::mcmc(): Requested negative or zero "
563 <<
"step_fac with aff_inv=true.\nSetting to 2.0."
567 std::cout <<
"mcmc_para::mcmc(): Requested negative or zero "
568 <<
"step_fac. Setting to 10.0." << std::endl;
575 unsigned long int seed=time(0);
576 if (this->user_seed!=0) {
581 rg[it].set_seed(seed);
603 std::vector<double> w_current(ssize);
604 for(
size_t i=0;i<ssize;i++) {
624 next[it].resize(n_params);
629 vec_t best(n_params);
634 std::vector<bool> mcmc_done_flag(
n_threads);
636 mcmc_done_flag[it]=
false;
647 O2SCL_ERR(
"Function mcmc_init() failed in mcmc_base::mcmc().",
658 scr_out <<
"mcmc: Affine-invariant step, n_params="
659 << n_params <<
", n_walk=" <<
n_walk
660 <<
", n_threads=" <<
n_threads <<
", rank="
664 scr_out <<
"mcmc: With proposal distribution, n_params="
665 << n_params <<
", n_threads=" <<
n_threads <<
", rank="
669 scr_out <<
"mcmc: Random-walk w/uniform dist., n_params="
670 << n_params <<
", n_threads=" <<
n_threads <<
", rank="
683 #pragma omp parallel default(shared)
710 for(
size_t ipar=0;ipar<n_params;ipar++) {
715 func_ret[it]=func[it](n_params,
current[sindex],
716 w_current[sindex],
data_arr[sindex]);
719 mcmc_done_flag[it]=
true;
729 meas_ret[it]=meas[it](
current[sindex],w_current[sindex],
736 mcmc_done_flag[it]=
true;
745 while (!done && !mcmc_done_flag[it]) {
748 for(
size_t ipar=0;ipar<n_params;ipar++) {
752 (
rg[it].random()*2.0-1.0)*
754 }
while (
current[sindex][ipar]>high[ipar] ||
755 current[sindex][ipar]<low[ipar]);
759 func_ret[it]=func[it](n_params,
current[sindex],
760 w_current[sindex],
data_arr[sindex]);
768 mcmc_done_flag[it]=
true;
774 scr_out <<
"Found initial guess for thread "
775 << it <<
". func_ret,weight,params=\n "
776 << func_ret[it] <<
" "
777 << w_current[sindex] <<
" ";
778 for(
size_t iji=0;iji<n_params;iji++) {
789 meas_ret[it]=meas[it](
current[sindex],w_current[sindex],
796 mcmc_done_flag[it]=
true;
800 std::string err=((std::string)
"In loop with thread ")+
803 " in mcmc_para_base::mcmc().";
814 bool stop_early=
false;
816 if (mcmc_done_flag[it]==
true) {
819 <<
"): Returned mcmc_done "
820 <<
"(initial; ai)." << std::endl;
837 if (w_current[sindex]>w_current[0]) {
839 w_best=w_current[sindex];
855 << w_current[sindex] <<
" (initial; ai)" << std::endl;
869 #pragma omp parallel default(shared)
884 for(
size_t ipar=0;ipar<n_params;ipar++) {
891 func_ret[it]=func[it](n_params,
current[it],w_current[it],
895 func_ret[it]=func_ret[it % ip_size];
896 w_current[it]=w_current[it % ip_size];
899 for(
size_t j=0;j<
data_arr.size();j++) {
914 <<
"): Initial point returned mcmc_done."
922 O2SCL_ERR((((std::string)
"Initial weight from thread ")+
924 " vanished in mcmc_para_base::mcmc().").c_str(),
935 #pragma omp parallel default(shared)
947 func_ret[it]=func_ret[it % ip_size];
949 w_current[it]=w_current[it % ip_size];
960 meas_ret[it]=meas[it](
current[it],w_current[it],0,
966 mcmc_done_flag[it]=
true;
977 bool stop_early=
false;
979 if (mcmc_done_flag[it]==
true) {
982 <<
"): Returned mcmc_done "
983 <<
"(initial)." << std::endl;
998 if (w_current[it]<w_best) {
1000 w_best=w_current[it];
1009 scr_out << w_current[it] <<
" ";
1011 scr_out <<
" (initial)" << std::endl;
1026 std::cout <<
"Press a key and type enter to continue. ";
1043 #pragma omp parallel default(shared)
1051 bool main_done=
false;
1052 size_t mcmc_iters=0;
1054 while (!main_done) {
1064 if (!std::isfinite(q_prop[it])) {
1065 O2SCL_ERR2(
"Proposal distribution not finite in ",
1072 for(
size_t k=0;k<n_params;k++) {
1074 next[it][k]=
current[it][k]+(
rg[it].random()*2.0-1.0)*
1077 next[it][k]=
current[it][k]+(
rg[it].random()*2.0-1.0)*
1091 for(
size_t k=0;k<n_params;k++) {
1092 if (next[it][k]<low[k] || next[it][k]>high[k]) {
1095 if (next[it][k]<low[k]) {
1096 std::cout <<
"mcmc (" << it <<
","
1097 <<
mpi_rank <<
"): Parameter with index " << k
1098 <<
" and value " << next[it][k]
1099 <<
" smaller than limit " << low[k] << std::endl;
1100 scr_out <<
"mcmc (" << it <<
","
1101 <<
mpi_rank <<
"): Parameter with index " << k
1102 <<
" and value " << next[it][k]
1103 <<
" smaller than limit " << low[k] << std::endl;
1105 std::cout <<
"mcmc (" << it <<
"," <<
mpi_rank
1106 <<
"): Parameter with index " << k
1107 <<
" and value " << next[it][k]
1108 <<
" larger than limit " << high[k] << std::endl;
1110 <<
"): Parameter with index " << k
1111 <<
" and value " << next[it][k]
1112 <<
" larger than limit " << high[k] << std::endl;
1122 func_ret[it]=func[it](n_params,next[it],w_next[it],
1126 func_ret[it]=func[it](n_params,next[it],w_next[it],
1130 mcmc_done_flag[it]=
true;
1150 double r=
rg[it].random();
1170 if (r<exp(w_next[it]-w_current[sindex]+q_prop[it])) {
1176 if (r<exp(w_next[it]-w_current[sindex])) {
1191 meas_ret[it]=meas[it](next[it],w_next[it],
1195 meas_ret[it]=meas[it](next[it],w_next[it],
1203 w_current[sindex]=w_next[it];
1215 meas_ret[it]=meas[it](next[it],w_next[it],
1219 meas_ret[it]=meas[it](next[it],w_next[it],
1234 #pragma omp critical (o2scl_mcmc_para_best_point)
1255 O2SCL_ERR((((std::string)
"Measurement function returned ")+
1257 " in mcmc_para_base::mcmc().").c_str(),
1265 if (main_done==
false) {
1275 scr_out <<
"o2scl::mcmc_para: Thread " << it
1276 <<
" finished warmup." << std::endl;
1286 scr_out <<
"o2scl::mcmc_para: Thread " << it
1287 <<
" stopping because number of iterations ("
1288 << mcmc_iters <<
") equal to max_iters (" <<
max_iters
1289 <<
")." << std::endl;
1294 if (main_done==
false) {
1303 scr_out <<
"o2scl::mcmc_para: Thread " << it
1304 <<
" stopping because elapsed (" << elapsed
1305 <<
") > max_time (" <<
max_time <<
")."
1326 bool main_done=
false;
1327 size_t mcmc_iters=0;
1329 while (!main_done) {
1334 #pragma omp parallel default(shared)
1359 ij=((size_t)(
rg[it].random()*((double)n_tot)));
1363 double p=
rg[it].random();
1365 smove_z[it]=(1.0-2.0*p+2.0*a*p+p*p-2.0*a*p*p+a*a*p*p)/a;
1375 std::cout <<
"Problem 1." << std::endl;
1376 std::cout << jt <<
" " <<
n_threads <<
" " << ij <<
" "
1383 for(
size_t i=0;i<n_params;i++) {
1386 std::cout <<
"Problem 2." << std::endl;
1387 std::cout << jt <<
" " << ij <<
" " << it <<
" "
1402 for(
size_t k=0;k<n_params;k++) {
1403 if (next[it][k]<low[k] || next[it][k]>high[k]) {
1406 if (next[it][k]<low[k]) {
1407 std::cout <<
"mcmc (" << it <<
","
1408 <<
mpi_rank <<
"): Parameter with index " << k
1409 <<
" and value " << next[it][k]
1410 <<
" smaller than limit " << low[k] << std::endl;
1411 scr_out <<
"mcmc (" << it <<
","
1412 <<
mpi_rank <<
"): Parameter with index " << k
1413 <<
" and value " << next[it][k]
1414 <<
" smaller than limit " << low[k] << std::endl;
1416 std::cout <<
"mcmc (" << it <<
"," <<
mpi_rank
1417 <<
"): Parameter with index " << k
1418 <<
" and value " << next[it][k]
1419 <<
" larger than limit " << high[k] << std::endl;
1421 <<
"): Parameter with index " << k
1422 <<
" and value " << next[it][k]
1423 <<
" larger than limit " << high[k] << std::endl;
1433 func_ret[it]=func[it](n_params,next[it],w_next[it],
1437 func_ret[it]=func[it](n_params,next[it],w_next[it],
1441 mcmc_done_flag[it]=
true;
1463 scr_out <<
"it: " << it <<
" q_prop[it]: "
1464 << q_prop[it] << std::endl;
1468 <<
"): Returned mcmc_done."
1472 <<
"): Parameter(s) out of range: " << std::endl;
1473 scr_out.setf(std::ios::showpos);
1474 for(
size_t k=0;k<n_params;k++) {
1475 scr_out << k <<
" " << low[k] <<
" "
1476 << next[it][k] <<
" " << high[k];
1477 if (next[it][k]<low[k] || next[it][k]>high[k]) {
1482 scr_out.unsetf(std::ios::showpos);
1487 <<
"): Function returned failure "
1488 << func_ret[it] <<
" at point ";
1489 for(
size_t k=0;k<n_params;k++) {
1490 scr_out << next[it][k] <<
" ";
1503 #pragma omp parallel default(shared)
1521 double r=
rg[it].random();
1523 double ai_ratio=pow(smove_z[it],((
double)n_params)-1.0)*
1524 exp(w_next[it]-w_current[sindex]);
1539 meas_ret[it]=meas[it](next[it],w_next[it],
1543 meas_ret[it]=meas[it](next[it],w_next[it],
1551 w_current[sindex]=w_next[it];
1562 meas_ret[it]=meas[it](next[it],w_next[it],
1566 meas_ret[it]=meas[it](next[it],w_next[it],
1587 scr_out.width((
int)(1.0+log10((
double)(n_params-1))));
1588 scr_out << mcmc_iters <<
" i_walk: "
1590 << w_current[sindex] << std::endl;
1617 O2SCL_ERR((((std::string)
"Measurement function returned ")+
1619 " in mcmc_para_base::mcmc().").c_str(),
1628 if (main_done==
false) {
1640 scr_out <<
"mcmc: Finished warmup." << std::endl;
1647 std::cout <<
"Press a key and type enter to continue. ";
1656 scr_out <<
"mcmc: Stopping because number of iterations ("
1657 << mcmc_iters <<
") equal to max_iters (" <<
max_iters
1658 <<
")." << std::endl;
1663 if (main_done==
false) {
1672 scr_out <<
"mcmc: Stopping because elapsed (" << elapsed
1673 <<
") > max_time (" <<
max_time <<
")."
1700 virtual int mcmc(
size_t n_params, vec_t &low, vec_t &high,
1701 func_t &func, measure_t &meas) {
1714 return mcmc(n_params,low,high,vf,vm);
1727 for(
size_t i=0;i<pv.size();i++) {
1743 for(
size_t i=0;i<pv.size();i++) {
1807 template<
class func_t,
class fill_t,
class data_t,
class vec_t=ubvector>
1809 std::function<int(const vec_t &,double,size_t,int,bool,data_t &)>,
1815 typedef std::function<int(
const vec_t &,
double,
size_t,
int,
bool,data_t &)>
1831 std::shared_ptr<o2scl::table_units<> >
table;
1850 table=std::shared_ptr<o2scl::table_units<> >
1863 for(
size_t i=0;i<
col_names.size();i++) {
1880 std::cout <<
"mcmc: Table column names and units: " << std::endl;
1890 O2SCL_ERR2(
"Flag 'prev_read' is true but table pointer is 0 ",
1898 std::string str=((std::string)
"Table does not have correct ")+
1899 "number of columns in mcmc_para_table::mcmc_init()."+
1909 O2SCL_ERR2(
"Table does not have the correct internal columns ",
1913 O2SCL_ERR2(
"Array walker_accept_rows does not have correct size ",
1917 O2SCL_ERR2(
"Array walker_reject_rows does not have correct size ",
1940 std::vector<double> &line, data_t &dat,
1941 size_t i_walker, fill_t &fill) {
1944 size_t i_thread=omp_get_thread_num();
1952 line.push_back(i_thread);
1954 line.push_back(i_walker);
1956 line.push_back(1.0);
1957 line.push_back(log_weight);
1958 for(
size_t i=0;i<pars.size();i++) {
1959 line.push_back(pars[i]);
1961 int tempi=fill(pars,log_weight,line,dat);
2042 this->
scr_out <<
"mcmc: Start write_files(). mpi_rank: "
2044 << this->
mpi_size <<
" table_io_chunk: "
2045 << table_io_chunk << std::endl;
2048 std::vector<o2scl::table_units<> > tab_arr;
2049 bool rank_sent=
false;
2053 if (this->
mpi_rank%table_io_chunk==0) {
2059 tab_arr.push_back(t);
2060 o2scl_table_mpi_recv(child,tab_arr[tab_arr.size()-1]);
2066 o2scl_table_mpi_send(*
table,parent);
2075 int tag=0, buffer=0;
2076 if (sync_write && this->
mpi_size>1 &&
2078 MPI_Recv(&buffer,1,MPI_INT,this->
mpi_rank-table_io_chunk,
2079 tag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
2096 hf.
set_szt(
"file_update_iters",this->file_update_iters);
2097 hf.
set_szt(
"file_update_time",this->file_update_time);
2100 hf.
set_szt(
"n_params",this->n_params);
2107 hf.
seti(
"store_rejects",this->store_rejects);
2108 hf.
seti(
"table_sequence",this->table_sequence);
2119 this->ret_value_counts[0].size(),
2120 this->ret_value_counts);
2123 this->initial_points[0].size(),
2124 this->initial_points);
2126 hf.
seti(
"n_tables",tab_arr.size()+1);
2127 if (rank_sent==
false) {
2128 hdf_output(hf,*
table,
"markov_chain_0");
2130 for(
size_t i=0;i<tab_arr.size();i++) {
2131 std::string name=((std::string)
"markov_chain_")+
szttos(i+1);
2132 hdf_output(hf,tab_arr[i],name);
2138 if (sync_write && this->
mpi_size>1 &&
2140 MPI_Send(&buffer,1,MPI_INT,this->
mpi_rank+table_io_chunk,
2141 tag,MPI_COMM_WORLD);
2146 this->
scr_out <<
"mcmc: Done write_files()." << std::endl;
2168 std::vector<std::string> units) {
2169 if (names.size()!=units.size()) {
2170 O2SCL_ERR2(
"Size of names and units arrays don't match in ",
2197 int tag=0, buffer=0;
2199 MPI_Recv(&buffer,1,MPI_INT,this->
mpi_rank-1,
2200 tag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
2212 MPI_Send(&buffer,1,MPI_INT,this->
mpi_rank+1,
2213 tag,MPI_COMM_WORLD);
2221 std::cout <<
"Initial points: Finding last " << n_points
2222 <<
" points from file named "
2223 << fname <<
" ." << std::endl;
2228 for(
size_t it=0;it<this->
n_threads;it++) {
2229 for(
size_t iw=0;iw<this->
n_walk;iw++) {
2232 size_t windex=it*this->n_walk+iw;
2235 for(
int row=tip.get_nlines()-1;row>=0 && found==
false;row--) {
2236 if (tip.get(
"walker",row)==iw &&
2237 tip.get(
"thread",row)==it &&
2238 tip.get(
"mult",row)>0.5) {
2242 std::cout <<
"Function initial_point_file_last():\n\tit: "
2243 << it <<
" rank: " << this->
mpi_rank
2244 <<
" iw: " << iw <<
" row: "
2245 << row <<
" log_wgt: " << tip.get(
"log_wgt",row)
2250 for(
size_t ip=0;ip<n_param_loc;ip++) {
2259 if (found==
false && tip.get_nlines()>this->n_walk*this->n_threads) {
2260 int row=tip.get_nlines()-this->n_walk*this->n_threads+windex;
2262 for(
size_t ip=0;ip<n_param_loc;ip++) {
2269 std::cout <<
"No initial guess found for rank " << this->
mpi_rank
2270 <<
" thread " << it <<
" and walker " << iw << std::endl;
2271 O2SCL_ERR(
"Function initial_points_file_last() failed.",
2298 int tag=0, buffer=0;
2300 MPI_Recv(&buffer,1,MPI_INT,this->
mpi_rank-1,
2301 tag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
2313 MPI_Send(&buffer,1,MPI_INT,this->
mpi_rank+1,
2314 tag,MPI_COMM_WORLD);
2322 std::cout <<
"Initial points: Finding last " << n_points
2323 <<
" points from file named "
2324 << fname <<
" ." << std::endl;
2329 size_t nlines=tip.get_nlines();
2330 size_t decrement=nlines/n_points;
2331 if (decrement<1) decrement=1;
2334 for(
size_t k=0;k<n_points;k++) {
2338 std::cout <<
"Function initial_point_file_dist():\n\trow: "
2339 << row <<
" log_wgt: " << tip.get(
"log_wgt",row)
2344 for(
size_t ip=0;ip<n_param_loc;ip++) {
2364 double thresh=1.0e-6,
2372 int tag=0, buffer=0;
2374 MPI_Recv(&buffer,1,MPI_INT,this->
mpi_rank-1,
2375 tag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
2387 MPI_Send(&buffer,1,MPI_INT,this->
mpi_rank+1,
2388 tag,MPI_COMM_WORLD);
2396 std::cout <<
"Initial points: Finding best " << n_points
2397 <<
" unique points from file named "
2398 << fname <<
" ." << std::endl;
2401 typedef std::map<double,int,std::greater<double> > map_t;
2405 for(
size_t k=0;k<tip.get_nlines();k++) {
2406 m.insert(std::make_pair(tip.get(
"log_wgt",k),k));
2416 for(map_t::iterator mit=m.begin();mit!=m.end();mit++) {
2417 map_t::iterator mit2=mit;
2419 if (mit2!=m.end()) {
2420 if (fabs(mit->first-mit2->first)<thresh) {
2422 std::cout <<
"Removing duplicate weights: "
2423 << mit->first <<
" " << mit2->first << std::endl;
2432 }
while (found==
true);
2435 if (m.size()<n_points) {
2436 O2SCL_ERR2(
"Could not find enough points in file in ",
2437 "mcmc_para::initial_points_file_best().",
2443 map_t::iterator mit=m.begin();
2444 for(
size_t k=0;k<n_points;k++) {
2445 int row=mit->second;
2447 std::cout <<
"Initial point " << k <<
" at row "
2448 << row <<
" has log_weight= "
2449 << tip.get(
"log_wgt",row) << std::endl;
2452 for(
size_t ip=0;ip<n_param_loc;ip++) {
2468 virtual int mcmc(
size_t n_params_local,
2469 vec_t &low, vec_t &high, std::vector<func_t> &func,
2470 std::vector<fill_t> &fill) {
2488 std::vector<internal_measure_t> meas(this->
n_threads);
2489 for(
size_t it=0;it<this->
n_threads;it++) {
2491 (std::mem_fn<
int(
const vec_t &,
double,
size_t,
int,
bool,
2492 data_t &,
size_t, fill_t &)>
2494 std::placeholders::_2,std::placeholders::_3,
2495 std::placeholders::_4,std::placeholders::_5,
2496 std::placeholders::_6,it,std::ref(fill[it]));
2524 chain_sizes.resize(ntot);
2526 for(
size_t it=0;it<this->
n_threads;it++) {
2527 for(
size_t iw=0;iw<this->
n_walk;iw++) {
2528 size_t ix=it*this->n_walk+iw;
2532 if (
table->
get(
"mult",j)>0.5) chain_sizes[ix]++;
2552 std::string name=
"") {
2565 O2SCL_ERR2(
"Table does not have the correct internal columns ",
2566 "in mcmc_para_table::read_prev_results().",
2578 walker_accept_rows.resize(ntot);
2581 for(
size_t j=0;j<ntot;j++) {
2588 size_t i_thread=((size_t)(
table->
get(
"thread",j)+1.0e-12));
2589 size_t i_walker=((size_t)(
table->
get(
"walker",j)+1.0e-12));
2592 size_t windex=i_thread*this->n_walk+i_walker;
2596 }
else if (
table->
get(
"mult",j)<-0.5) {
2605 for(
size_t j=0;j<ntot;j++) {
2613 for(
size_t j=0;j<ntot;j++) {
2615 for(
size_t k=0;k<n_param_loc;k++) {
2620 std::cout <<
"Previous table was read, but initial points not set."
2625 std::cout <<
"mcmc_para_table::read_prev_results():" << std::endl;
2626 std::cout <<
" index walker_accept_rows walker_reject_rows"
2628 for(
size_t j=0;j<ntot;j++) {
2631 std::cout << j <<
" ";
2655 size_t total_iters=0;
2656 for(
size_t it=0;it<this->
n_threads;it++) {
2661 this->
scr_out <<
"mcmc: Writing to file. total_iters: "
2662 << total_iters <<
" file_update_iters: "
2679 this->
scr_out <<
"mcmc: Writing to file. elapsed: "
2680 << elapsed <<
" file_update_time: "
2700 virtual int add_line(
const vec_t &pars,
double log_weight,
2701 size_t walker_ix,
int func_ret,
2702 bool mcmc_accept, data_t &dat,
2703 size_t i_thread, fill_t &fill) {
2706 size_t windex=i_thread*this->
n_walk+walker_ix;
2734 #pragma omp critical (o2scl_mcmc_para_table_add_line)
2739 fabs(
table->
get(
"mult",next_row))>0.1) {
2754 for(
size_t i=0;i<this->
n_walk;i++) {
2756 table->set(
"thread",istart+j*this->n_walk+i,j);
2757 table->
set(
"walker",istart+j*this->n_walk+i,i);
2758 table->
set(
"mult",istart+j*this->n_walk+i,0.0);
2759 table->
set(
"log_wgt",istart+j*this->n_walk+i,0.0);
2771 std::vector<double> line;
2772 int fret=
fill_line(pars,log_weight,line,dat,walker_ix,fill);
2793 std::cout <<
"line: " << line.size() <<
" columns: "
2796 std::cout << k <<
". ";
2797 if (k<table->get_ncolumns()) {
2801 if (k<line.size()) std::cout << line[k] <<
" ";
2802 std::cout << std::endl;
2804 O2SCL_ERR(
"Table misalignment in mcmc_para_table::add_line().",
2813 this->
scr_out <<
"mcmc: Setting data at row " << next_row
2815 for(
size_t k=0;k<line.size();k++) {
2817 this->
scr_out << table->get_column_name(k) <<
" ";
2819 this->
scr_out <<
" " << line[k] << std::endl;
2837 O2SCL_ERR2(
"Invalid row for incrementing multiplier in ",
2847 this->
scr_out <<
"mcmc: Updating mult of row "
2848 << walker_accept_rows[windex]
2849 <<
" from " << mult_old <<
" to "
2850 << mult_old+1.0 << std::endl;
2879 if (fabs(
table->
get(
"mult",i))<0.1) {
2897 std::vector<double> &ac_coeff_avg,
2898 int loc_verbose=0) {
2900 ac_coeff_avg.resize(0);
2904 std::vector<std::vector<double> > ac_coeff_temp(n_tot);
2907 for(
size_t k=0;k<this->
n_walk;k++) {
2908 size_t tindex=j*this->n_walk+k;
2909 std::vector<double> vec, mult;
2911 if (fabs(
table->
get(
"walker",ell)-k)<0.1 &&
2912 fabs(
table->
get(
"thread",ell)-j)<0.1 &&
2914 mult.push_back(
table->
get(
"mult",ell));
2915 vec.push_back(
table->
get(icol,ell));
2918 if (mult.size()>1) {
2920 (vec.size(),vec,mult,ac_coeff_temp[tindex]);
2921 if (ac_coeff_temp[tindex].size()>max_size) {
2922 max_size=ac_coeff_temp[tindex].size();
2925 if (loc_verbose>0) {
2926 std::cout <<
"column index, thread, walker, length, mult. sum: "
2927 << icol <<
" " << j <<
" " << k <<
" "
2928 << mult.size() <<
" "
2929 << o2scl::vector_sum<std::vector<double>,
double>
2930 (mult.size(),mult) <<
" " << max_size << std::endl;
2934 for(
size_t j=0;j<max_size;j++) {
2937 for(
size_t k=0;k<n_tot;k++) {
2938 if (j<ac_coeff_temp[k].size()) {
2939 ac_sum+=ac_coeff_temp[k][j];
2944 ac_coeff_avg.push_back(ac_sum/((
double)ac_count));
2956 std::shared_ptr<o2scl::table_units<> > table2=
2960 for(
size_t j=0;j<this->
n_walk;j++) {
2961 std::string func=std::string(
"abs(walker-")+
o2scl::szttos(j)+
2985 for(
size_t it=0;it<this->
n_threads;it++) {
2989 O2SCL_ERR2(
"Cannot reblock. Not enough data in ",
2992 size_t n_block=n/n_blocks;
2994 for(
size_t j=0;j<n_blocks;j++) {
2997 for(
size_t i=0;i<m;i++) {
3000 for(
size_t k=j*n_block;k<(j+1)*n_block;k++) {
3001 mult+=(*table)[
"mult"][k];
3002 for(
size_t i=1;i<m;i++) {
3003 dat[i]+=(*table)[i][k]*(*table)[
"mult"][k];
3007 for(
size_t i=1;i<m;i++) {
3028 template<
class func_t,
class fill_t,
class data_t,
class vec_t=ubvector>
3110 p_file_update_iters.
help=((std::string)
"Number of MCMC successes ")+
3111 "between file updates (default 0 for no file updates).";
3112 cl.
par_list.insert(std::make_pair(
"file_update_iters",
3113 &p_file_update_iters));
3116 p_file_update_time.
help=((std::string)
"Time ")+
3117 "between file updates (default 0.0 for no file updates).";
3118 cl.
par_list.insert(std::make_pair(
"file_update_time",
3119 &p_file_update_time));
3130 p_max_time.
help=((std::string)
"Maximum run time in seconds ")+
3131 "(default 86400 sec or 1 day).";
3132 cl.
par_list.insert(std::make_pair(
"max_time",&p_max_time));
3135 p_max_iters.
help=((std::string)
"If non-zero, limit the number of ")+
3136 "iterations to be less than the specified number (default zero).";
3137 cl.
par_list.insert(std::make_pair(
"max_iters",&p_max_iters));
3140 p_prefix.
help=
"Output file prefix (default 'mcmc\').";
3141 cl.
par_list.insert(std::make_pair(
"prefix",&p_prefix));
3151 p_step_fac.
help=((std::string)
"MCMC step factor. The step size for ")+
3152 "each variable is taken as the difference between the high and low "+
3153 "limits divided by this factor (default 10.0). This factor can "+
3154 "be increased if the acceptance rate is too small, but care must "+
3155 "be taken, e.g. if the conditional probability is multimodal. If "+
3156 "this step size is smaller than 1.0, it is reset to 1.0 .";
3157 cl.
par_list.insert(std::make_pair(
"step_fac",&p_step_fac));
3160 p_n_warm_up.
help=((std::string)
"Minimum number of warm up iterations ")+
3162 cl.
par_list.insert(std::make_pair(
"n_warm_up",&p_n_warm_up));
3165 p_verbose.
help=((std::string)
"MCMC verbosity parameter ")+
3167 cl.
par_list.insert(std::make_pair(
"mcmc_verbose",&p_verbose));
3170 p_max_bad_steps.
help=((std::string)
"Maximum number of bad steps ")+
3172 cl.
par_list.insert(std::make_pair(
"max_bad_steps",&p_max_bad_steps));
3175 p_n_walk.
help=((std::string)
"Number of walkers ")+
3177 cl.
par_list.insert(std::make_pair(
"n_walk",&p_n_walk));
3180 p_user_seed.
help=((std::string)
"Seed for multiplier for random ")+
3181 "number generator. If zero is given (the default), then mcmc() "+
3182 "uses time(0) to generate a random seed.";
3183 cl.
par_list.insert(std::make_pair(
"user_seed",&p_user_seed));
3186 p_aff_inv.
help=((std::string)
"If true, then use affine-invariant ")+
3187 "sampling (default false).";
3188 cl.
par_list.insert(std::make_pair(
"aff_inv",&p_aff_inv));
3191 p_table_sequence.
help=((std::string)
"If true, then ensure equal spacing ")+
3192 "between walkers (default true).";
3193 cl.
par_list.insert(std::make_pair(
"table_sequence",&p_table_sequence));
3196 p_store_rejects.
help=((std::string)
"If true, then store MCMC rejections ")+
3198 cl.
par_list.insert(std::make_pair(
"store_rejects",&p_store_rejects));