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,
134 class data_t,
class vec_t=ubvector>
class mcmc_para_base {
151 std::vector<rng_gsl>
rg;
154 std::vector<o2scl::prob_cond_mdim<vec_t> *>
prop_dist;
197 O2SCL_ERR2(
"Using a proposal distribution with affine-invariant ",
198 "sampling not implemented in mcmc_para::mcmc_init().",
203 std::cout <<
"Prefix is: " <<
prefix << std::endl;
210 scr_out.setf(std::ios::scientific);
224 <<
" reject=" <<
n_reject[ic] << std::endl;
233 virtual void best_point(vec_t &best,
double w_best, data_t &dat) {
405 MPI_Comm_rank(MPI_COMM_WORLD,&this->mpi_rank);
406 MPI_Comm_size(MPI_COMM_WORLD,&this->mpi_size);
434 virtual int mcmc(
size_t n_params, vec_t &low, vec_t &high,
435 std::vector<func_t> &func, std::vector<measure_t> &meas) {
441 if (func.size()==0 || meas.size()==0) {
442 O2SCL_ERR2(
"Size of 'func' or 'meas' array is zero in ",
447 std::cout <<
"mcmc_para::mcmc(): Not enough functions for "
448 <<
n_threads <<
" threads. Setting n_threads to "
449 << func.size() <<
"." << std::endl;
455 std::cout <<
"mcmc_para::mcmc(): Not enough measurement objects for "
456 <<
n_threads <<
" threads. Setting n_threads to "
457 << meas.size() <<
"." << std::endl;
476 for(
size_t k=0;k<n_params;k++) {
485 for(
size_t ipar=0;ipar<n_params;ipar++) {
493 ") in mcmc_base::mcmc().").c_str(),
519 std::cout.setf(std::ios::scientific);
520 std::cout << i <<
" ";
522 std::cout << j <<
" ";
538 std::cout <<
"mcmc_para::mcmc(): "
539 <<
n_threads <<
" threads were requested but the "
540 <<
"-DO2SCL_OPENMP flag was not used during "
541 <<
"compilation. Setting n_threads to 1."
553 std::cout <<
"mcmc_para::mcmc(): Requested only 1 walker, "
554 <<
"forcing 2 walkers for each parameter." << std::endl;
559 std::cout <<
"mcmc_para::mcmc(): Automatically setting "
560 <<
"n_walk_per_thread to " <<
n_walk <<
"." << std::endl;
564 O2SCL_ERR2(
"More walkers per thread than total walkers ",
593 std::cout <<
"mcmc_para::mcmc(): Could not evenly "
594 <<
"organize threads and walkers." << std::endl;
595 std::cout <<
"n_threads: " <<
n_threads << std::endl;
596 std::cout <<
"n_walk: " <<
n_walk << std::endl;
597 std::cout <<
"n_walk_per_thread: "
599 std::cout <<
"n_chains_per_rank: "
601 O2SCL_ERR2(
"Values of n_walk, n_threads, and n_walk_per_thread ",
602 "not commensurate in mcmc_para::mcmc().",
617 std::cout <<
"mcmc_para::mcmc(): Requested negative or zero "
618 <<
"step_fac with aff_inv=true.\nSetting to 2.0."
622 std::cout <<
"mcmc_para::mcmc(): Requested negative or zero "
623 <<
"step_fac. Setting to 10.0." << std::endl;
630 unsigned long int seed=time(0);
631 if (this->user_seed!=0) {
636 rg[it].set_seed(seed);
658 std::vector<double> w_current(ssize);
659 for(
size_t i=0;i<ssize;i++) {
679 next[it].resize(n_params);
684 vec_t best(n_params);
689 std::vector<bool> mcmc_done_flag(
n_threads);
691 mcmc_done_flag[it]=
false;
702 O2SCL_ERR(
"Function mcmc_init() failed in mcmc_base::mcmc().",
713 scr_out <<
"mcmc: Affine-invariant step, n_params="
714 << n_params <<
", n_walk=" <<
n_walk
717 <<
", n_threads=" <<
n_threads <<
", rank="
721 scr_out <<
"mcmc: With proposal distribution, n_params="
722 << n_params <<
", n_threads=" <<
n_threads <<
", rank="
726 scr_out <<
"mcmc: Random-walk w/uniform dist., n_params="
727 << n_params <<
", n_threads=" <<
n_threads <<
", rank="
740 #pragma omp parallel default(shared)
767 for(
size_t ipar=0;ipar<n_params;ipar++) {
772 func_ret[it]=func[it](n_params,
current[sindex],
773 w_current[sindex],
data_arr[sindex]);
776 mcmc_done_flag[it]=
true;
786 meas_ret[it]=meas[it](
current[sindex],w_current[sindex],
793 mcmc_done_flag[it]=
true;
802 while (!done && !mcmc_done_flag[it]) {
805 for(
size_t ipar=0;ipar<n_params;ipar++) {
809 (
rg[it].random()*2.0-1.0)*
811 }
while (
current[sindex][ipar]>high[ipar] ||
812 current[sindex][ipar]<low[ipar]);
816 func_ret[it]=func[it](n_params,
current[sindex],
817 w_current[sindex],
data_arr[sindex]);
825 mcmc_done_flag[it]=
true;
836 meas_ret[it]=meas[it](
current[sindex],w_current[sindex],
843 mcmc_done_flag[it]=
true;
847 std::string err=((std::string)
"In loop with thread ")+
850 " in mcmc_para_base::mcmc().";
861 bool stop_early=
false;
863 if (mcmc_done_flag[it]==
true) {
866 <<
"): Returned mcmc_done "
867 <<
"(initial; ai)." << std::endl;
884 if (w_current[sindex]>w_current[0]) {
886 w_best=w_current[sindex];
903 << w_current[sindex] <<
" (initial; ai)" << std::endl;
917 #pragma omp parallel default(shared)
932 for(
size_t ipar=0;ipar<n_params;ipar++) {
939 func_ret[it]=func[it](n_params,
current[it],w_current[it],
943 func_ret[it]=func_ret[it % ip_size];
944 w_current[it]=w_current[it % ip_size];
947 for(
size_t j=0;j<
data_arr.size();j++) {
962 <<
"): Initial point returned mcmc_done."
970 O2SCL_ERR((((std::string)
"Initial weight from thread ")+
972 " vanished in mcmc_para_base::mcmc().").c_str(),
983 #pragma omp parallel default(shared)
995 func_ret[it]=func_ret[it % ip_size];
997 w_current[it]=w_current[it % ip_size];
1008 meas_ret[it]=meas[it](
current[it],w_current[it],0,
1014 mcmc_done_flag[it]=
true;
1025 bool stop_early=
false;
1027 if (mcmc_done_flag[it]==
true) {
1030 <<
"): Returned mcmc_done "
1031 <<
"(initial)." << std::endl;
1043 w_best=w_current[0];
1046 if (w_current[it]<w_best) {
1048 w_best=w_current[it];
1057 scr_out << w_current[it] <<
" ";
1059 scr_out <<
" (initial)" << std::endl;
1074 std::cout <<
"Press a key and type enter to continue. ";
1091 #pragma omp parallel default(shared)
1099 bool main_done=
false;
1100 size_t mcmc_iters=0;
1102 while (!main_done) {
1112 if (!std::isfinite(q_prop[it])) {
1113 O2SCL_ERR2(
"Proposal distribution not finite in ",
1120 for(
size_t k=0;k<n_params;k++) {
1122 next[it][k]=
current[it][k]+(
rg[it].random()*2.0-1.0)*
1125 next[it][k]=
current[it][k]+(
rg[it].random()*2.0-1.0)*
1139 for(
size_t k=0;k<n_params;k++) {
1140 if (next[it][k]<low[k] || next[it][k]>high[k]) {
1143 if (next[it][k]<low[k]) {
1144 std::cout <<
"mcmc (" << it <<
","
1145 <<
mpi_rank <<
"): Parameter with index " << k
1146 <<
" and value " << next[it][k]
1147 <<
" smaller than limit " << low[k] << std::endl;
1148 scr_out <<
"mcmc (" << it <<
","
1149 <<
mpi_rank <<
"): Parameter with index " << k
1150 <<
" and value " << next[it][k]
1151 <<
" smaller than limit " << low[k] << std::endl;
1153 std::cout <<
"mcmc (" << it <<
"," <<
mpi_rank
1154 <<
"): Parameter with index " << k
1155 <<
" and value " << next[it][k]
1156 <<
" larger than limit " << high[k] << std::endl;
1158 <<
"): Parameter with index " << k
1159 <<
" and value " << next[it][k]
1160 <<
" larger than limit " << high[k] << std::endl;
1170 func_ret[it]=func[it](n_params,next[it],w_next[it],
1174 func_ret[it]=func[it](n_params,next[it],w_next[it],
1178 mcmc_done_flag[it]=
true;
1198 double r=
rg[it].random();
1218 if (r<exp(w_next[it]-w_current[sindex]+q_prop[it])) {
1224 if (r<exp(w_next[it]-w_current[sindex])) {
1239 meas_ret[it]=meas[it](next[it],w_next[it],
1243 meas_ret[it]=meas[it](next[it],w_next[it],
1251 w_current[sindex]=w_next[it];
1263 meas_ret[it]=meas[it](next[it],w_next[it],
1267 meas_ret[it]=meas[it](next[it],w_next[it],
1282 #pragma omp critical (o2scl_mcmc_para_best_point)
1303 O2SCL_ERR((((std::string)
"Measurement function returned ")+
1305 " in mcmc_para_base::mcmc().").c_str(),
1313 if (main_done==
false) {
1323 scr_out <<
"o2scl::mcmc_para: Thread " << it
1324 <<
" finished warmup." << std::endl;
1334 scr_out <<
"o2scl::mcmc_para: Thread " << it
1335 <<
" stopping because number of iterations ("
1336 << mcmc_iters <<
") equal to max_iters (" <<
max_iters
1337 <<
")." << std::endl;
1342 if (main_done==
false) {
1351 scr_out <<
"o2scl::mcmc_para: Thread " << it
1352 <<
" stopping because elapsed (" << elapsed
1353 <<
") > max_time (" <<
max_time <<
")."
1374 bool main_done=
false;
1375 size_t mcmc_iters=0;
1377 while (!main_done) {
1383 #pragma omp parallel default(shared)
1400 ij=((size_t)(
rg[it].random()*((double)
n_walk)));
1404 double p=
rg[it].random();
1406 smove_z[it]=(1.0-2.0*p+2.0*a*p+p*p-2.0*a*p*p+a*a*p*p)/a;
1415 for(
size_t i=0;i<n_params;i++) {
1427 for(
size_t k=0;k<n_params;k++) {
1428 if (next[it][k]<low[k] || next[it][k]>high[k]) {
1431 if (next[it][k]<low[k]) {
1432 std::cout <<
"mcmc (" << it <<
","
1433 <<
mpi_rank <<
"): Parameter with index " << k
1434 <<
" and value " << next[it][k]
1435 <<
" smaller than limit " << low[k] << std::endl;
1436 scr_out <<
"mcmc (" << it <<
","
1437 <<
mpi_rank <<
"): Parameter with index " << k
1438 <<
" and value " << next[it][k]
1439 <<
" smaller than limit " << low[k] << std::endl;
1441 std::cout <<
"mcmc (" << it <<
"," <<
mpi_rank
1442 <<
"): Parameter with index " << k
1443 <<
" and value " << next[it][k]
1444 <<
" larger than limit " << high[k] << std::endl;
1446 <<
"): Parameter with index " << k
1447 <<
" and value " << next[it][k]
1448 <<
" larger than limit " << high[k] << std::endl;
1458 func_ret[it]=func[it]
1459 (n_params,next[it],w_next[it],
1463 func_ret[it]=func[it]
1464 (n_params,next[it],w_next[it],
1468 mcmc_done_flag[it]=
true;
1490 scr_out <<
"it: " << it <<
" q_prop[it]: "
1491 << q_prop[it] << std::endl;
1495 <<
"): Returned mcmc_done."
1499 <<
"): Parameter(s) out of range: " << std::endl;
1500 scr_out.setf(std::ios::showpos);
1501 for(
size_t k=0;k<n_params;k++) {
1502 scr_out << k <<
" " << low[k] <<
" "
1503 << next[it][k] <<
" " << high[k];
1504 if (next[it][k]<low[k] || next[it][k]>high[k]) {
1509 scr_out.unsetf(std::ios::showpos);
1514 <<
"): Function returned failure "
1515 << func_ret[it] <<
" at point ";
1516 for(
size_t k=0;k<n_params;k++) {
1517 scr_out << next[it][k] <<
" ";
1530 #pragma omp parallel default(shared)
1548 double r=
rg[it].random();
1551 double ai_ratio=pow(smove_z[it],((
double)n_params)-1.0)*
1552 exp(w_next[it]-w_current[sindex]);
1557 if (r<exp(w_next[it]-w_current[sindex]+q_prop[it])) {
1562 if (r<exp(w_next[it]-w_current[sindex])) {
1577 meas_ret[it]=meas[it]
1578 (next[it],w_next[it],
curr_walker[it],func_ret[it],
true,
1581 meas_ret[it]=meas[it](next[it],w_next[it],
1589 w_current[sindex]=w_next[it];
1600 meas_ret[it]=meas[it]
1601 (next[it],w_next[it],
curr_walker[it],func_ret[it],
false,
1604 meas_ret[it]=meas[it](next[it],w_next[it],
1625 scr_out.width((
int)(1.0+log10((
double)(n_params-1))));
1626 scr_out << mcmc_iters <<
" i_walk: "
1628 << w_current[sindex] << std::endl;
1657 O2SCL_ERR((((std::string)
"Measurement function returned ")+
1659 " in mcmc_para_base::mcmc().").c_str(),
1668 if (main_done==
false) {
1680 scr_out <<
"mcmc: Finished warmup." << std::endl;
1687 std::cout <<
"Press a key and type enter to continue. ";
1696 scr_out <<
"mcmc: Stopping because number of iterations ("
1697 << mcmc_iters <<
") equal to max_iters (" <<
max_iters
1698 <<
")." << std::endl;
1703 if (main_done==
false) {
1712 scr_out <<
"mcmc: Stopping because elapsed (" << elapsed
1713 <<
") > max_time (" <<
max_time <<
")."
1740 virtual int mcmc(
size_t n_params, vec_t &low, vec_t &high,
1741 func_t &func, measure_t &meas) {
1754 return mcmc(n_params,low,high,vf,vm);
1767 for(
size_t i=0;i<pv.size();i++) {
1783 for(
size_t i=0;i<pv.size();i++) {
1847 template<
class func_t,
class fill_t,
class data_t,
class vec_t=ubvector>
1848 class mcmc_para_table :
public mcmc_para_base<func_t,
1849 std::function<int(const vec_t &,double,size_t,int,bool,data_t &)>,
1855 typedef std::function<int(
const vec_t &,
double,
size_t,
int,
bool,data_t &)>
1871 std::shared_ptr<o2scl::table_units<> >
table;
1890 table=std::shared_ptr<o2scl::table_units<> >
1903 for(
size_t i=0;i<
col_names.size();i++) {
1920 std::cout <<
"mcmc: Table column names and units: " << std::endl;
1930 O2SCL_ERR2(
"Flag 'prev_read' is true but table pointer is 0 ",
1938 std::string str=((std::string)
"Table does not have correct ")+
1939 "number of columns in mcmc_para_table::mcmc_init()."+
1949 O2SCL_ERR2(
"Table does not have the correct internal columns ",
1953 O2SCL_ERR2(
"Array walker_accept_rows does not have correct size ",
1957 O2SCL_ERR2(
"Array walker_reject_rows does not have correct size ",
1980 std::vector<double> &line, data_t &dat,
1981 size_t i_walker, fill_t &fill) {
1984 size_t i_thread=omp_get_thread_num();
1992 line.push_back(i_thread);
1994 line.push_back(i_walker);
1996 line.push_back(1.0);
1997 line.push_back(log_weight);
1998 for(
size_t i=0;i<pars.size();i++) {
1999 line.push_back(pars[i]);
2001 int tempi=fill(pars,log_weight,line,dat);
2082 this->
scr_out <<
"mcmc: Start write_files(). mpi_rank: "
2084 << this->
mpi_size <<
" table_io_chunk: "
2085 << table_io_chunk << std::endl;
2088 std::vector<o2scl::table_units<> > tab_arr;
2089 bool rank_sent=
false;
2093 if (this->
mpi_rank%table_io_chunk==0) {
2099 tab_arr.push_back(t);
2100 o2scl_table_mpi_recv(child,tab_arr[tab_arr.size()-1]);
2106 o2scl_table_mpi_send(*
table,parent);
2115 int tag=0, buffer=0;
2116 if (sync_write && this->
mpi_size>1 &&
2118 MPI_Recv(&buffer,1,MPI_INT,this->
mpi_rank-table_io_chunk,
2119 tag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
2136 hf.
set_szt(
"file_update_iters",this->file_update_iters);
2137 hf.
set_szt(
"file_update_time",this->file_update_time);
2141 hf.
set_szt(
"n_params",this->n_params);
2149 hf.
seti(
"store_rejects",this->store_rejects);
2150 hf.
seti(
"table_sequence",this->table_sequence);
2161 this->ret_value_counts[0].size(),
2162 this->ret_value_counts);
2165 this->initial_points[0].size(),
2166 this->initial_points);
2168 hf.
seti(
"n_tables",tab_arr.size()+1);
2169 if (rank_sent==
false) {
2170 hdf_output(hf,*
table,
"markov_chain_0");
2172 for(
size_t i=0;i<tab_arr.size();i++) {
2173 std::string name=((std::string)
"markov_chain_")+
szttos(i+1);
2174 hdf_output(hf,tab_arr[i],name);
2180 if (sync_write && this->
mpi_size>1 &&
2182 MPI_Send(&buffer,1,MPI_INT,this->
mpi_rank+table_io_chunk,
2183 tag,MPI_COMM_WORLD);
2188 this->
scr_out <<
"mcmc: Done write_files()." << std::endl;
2210 std::vector<std::string> units) {
2211 if (names.size()!=units.size()) {
2212 O2SCL_ERR2(
"Size of names and units arrays don't match in ",
2239 int tag=0, buffer=0;
2241 MPI_Recv(&buffer,1,MPI_INT,this->
mpi_rank-1,
2242 tag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
2254 MPI_Send(&buffer,1,MPI_INT,this->
mpi_rank+1,
2255 tag,MPI_COMM_WORLD);
2263 std::cout <<
"Initial points: Finding last " << n_points
2264 <<
" points from file named "
2265 << fname <<
" ." << std::endl;
2270 for(
size_t it=0;it<this->
n_threads;it++) {
2271 for(
size_t iw=0;iw<this->
n_walk;iw++) {
2274 size_t windex=it*this->n_walk+iw;
2277 for(
int row=tip.get_nlines()-1;row>=0 && found==
false;row--) {
2278 if (tip.get(
"walker",row)==iw &&
2279 tip.get(
"thread",row)==it &&
2280 tip.get(
"mult",row)>0.5) {
2284 std::cout <<
"Function initial_point_file_last():\n\tit: "
2285 << it <<
" rank: " << this->
mpi_rank
2286 <<
" iw: " << iw <<
" row: "
2287 << row <<
" log_wgt: " << tip.get(
"log_wgt",row)
2292 for(
size_t ip=0;ip<n_param_loc;ip++) {
2301 if (found==
false && tip.get_nlines()>this->n_walk*this->n_threads) {
2302 int row=tip.get_nlines()-this->n_walk*this->n_threads+windex;
2304 for(
size_t ip=0;ip<n_param_loc;ip++) {
2311 std::cout <<
"No initial guess found for rank " << this->
mpi_rank
2312 <<
" thread " << it <<
" and walker " << iw << std::endl;
2313 O2SCL_ERR(
"Function initial_points_file_last() failed.",
2340 int tag=0, buffer=0;
2342 MPI_Recv(&buffer,1,MPI_INT,this->
mpi_rank-1,
2343 tag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
2355 MPI_Send(&buffer,1,MPI_INT,this->
mpi_rank+1,
2356 tag,MPI_COMM_WORLD);
2364 std::cout <<
"Initial points: Finding last " << n_points
2365 <<
" points from file named "
2366 << fname <<
" ." << std::endl;
2371 size_t nlines=tip.get_nlines();
2372 size_t decrement=nlines/n_points;
2373 if (decrement<1) decrement=1;
2376 for(
size_t k=0;k<n_points;k++) {
2380 std::cout <<
"Function initial_point_file_dist():\n\trow: "
2381 << row <<
" log_wgt: " << tip.get(
"log_wgt",row)
2386 for(
size_t ip=0;ip<n_param_loc;ip++) {
2406 double thresh=1.0e-6,
2414 int tag=0, buffer=0;
2416 MPI_Recv(&buffer,1,MPI_INT,this->
mpi_rank-1,
2417 tag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
2429 MPI_Send(&buffer,1,MPI_INT,this->
mpi_rank+1,
2430 tag,MPI_COMM_WORLD);
2438 std::cout <<
"Initial points: Finding best " << n_points
2439 <<
" unique points from file named "
2440 << fname <<
" ." << std::endl;
2443 typedef std::map<double,int,std::greater<double> > map_t;
2447 for(
size_t k=0;k<tip.get_nlines();k++) {
2448 m.insert(std::make_pair(tip.get(
"log_wgt",k),k));
2458 for(map_t::iterator mit=m.begin();mit!=m.end();mit++) {
2459 map_t::iterator mit2=mit;
2461 if (mit2!=m.end()) {
2462 if (fabs(mit->first-mit2->first)<thresh) {
2464 std::cout <<
"Removing duplicate weights: "
2465 << mit->first <<
" " << mit2->first << std::endl;
2474 }
while (found==
true);
2477 if (m.size()<n_points) {
2478 O2SCL_ERR2(
"Could not find enough points in file in ",
2479 "mcmc_para::initial_points_file_best().",
2485 map_t::iterator mit=m.begin();
2486 for(
size_t k=0;k<n_points;k++) {
2487 int row=mit->second;
2489 std::cout <<
"Initial point " << k <<
" at row "
2490 << row <<
" has log_weight= "
2491 << tip.get(
"log_wgt",row) << std::endl;
2494 for(
size_t ip=0;ip<n_param_loc;ip++) {
2510 virtual int mcmc(
size_t n_params_local,
2511 vec_t &low, vec_t &high, std::vector<func_t> &func,
2512 std::vector<fill_t> &fill) {
2530 std::vector<internal_measure_t> meas(this->
n_threads);
2531 for(
size_t it=0;it<this->
n_threads;it++) {
2533 (std::mem_fn<
int(
const vec_t &,
double,
size_t,
int,
bool,
2534 data_t &,
size_t, fill_t &)>
2536 std::placeholders::_2,std::placeholders::_3,
2537 std::placeholders::_4,std::placeholders::_5,
2538 std::placeholders::_6,it,std::ref(fill[it]));
2566 chain_sizes.resize(ntot);
2568 for(
size_t it=0;it<this->
n_threads;it++) {
2569 for(
size_t iw=0;iw<this->
n_walk;iw++) {
2570 size_t ix=it*this->n_walk+iw;
2574 if (
table->
get(
"mult",j)>0.5) chain_sizes[ix]++;
2594 std::string name=
"") {
2607 O2SCL_ERR2(
"Table does not have the correct internal columns ",
2608 "in mcmc_para_table::read_prev_results().",
2620 walker_accept_rows.resize(ntot);
2623 for(
size_t j=0;j<ntot;j++) {
2630 size_t i_thread=((size_t)(
table->
get(
"thread",j)+1.0e-12));
2631 size_t i_walker=((size_t)(
table->
get(
"walker",j)+1.0e-12));
2634 size_t windex=i_thread*this->n_walk+i_walker;
2638 }
else if (
table->
get(
"mult",j)<-0.5) {
2647 for(
size_t j=0;j<ntot;j++) {
2655 for(
size_t j=0;j<ntot;j++) {
2657 for(
size_t k=0;k<n_param_loc;k++) {
2662 std::cout <<
"Previous table was read, but initial points not set."
2667 std::cout <<
"mcmc_para_table::read_prev_results():" << std::endl;
2668 std::cout <<
" index walker_accept_rows walker_reject_rows"
2670 for(
size_t j=0;j<ntot;j++) {
2673 std::cout << j <<
" ";
2697 size_t total_iters=0;
2703 this->
scr_out <<
"mcmc: Writing to file. total_iters: "
2704 << total_iters <<
" file_update_iters: "
2721 this->
scr_out <<
"mcmc: Writing to file. elapsed: "
2722 << elapsed <<
" file_update_time: "
2742 virtual int add_line(
const vec_t &pars,
double log_weight,
2743 size_t walker_ix,
int func_ret,
2744 bool mcmc_accept, data_t &dat,
2745 size_t i_thread, fill_t &fill) {
2748 size_t windex=i_thread*this->
n_walk+walker_ix;
2776 #pragma omp critical (o2scl_mcmc_para_table_add_line)
2781 fabs(
table->
get(
"mult",next_row))>0.1) {
2796 for(
size_t i=0;i<this->
n_walk;i++) {
2798 table->set(
"thread",istart+j*this->n_walk+i,j);
2799 table->
set(
"walker",istart+j*this->n_walk+i,i);
2800 table->
set(
"mult",istart+j*this->n_walk+i,0.0);
2801 table->
set(
"log_wgt",istart+j*this->n_walk+i,0.0);
2813 std::vector<double> line;
2814 int fret=
fill_line(pars,log_weight,line,dat,walker_ix,fill);
2835 std::cout <<
"line: " << line.size() <<
" columns: "
2838 std::cout << k <<
". ";
2839 if (k<table->get_ncolumns()) {
2843 if (k<line.size()) std::cout << line[k] <<
" ";
2844 std::cout << std::endl;
2846 O2SCL_ERR(
"Table misalignment in mcmc_para_table::add_line().",
2855 this->
scr_out <<
"mcmc: Setting data at row " << next_row
2857 for(
size_t k=0;k<line.size();k++) {
2859 this->
scr_out << table->get_column_name(k) <<
" ";
2861 this->
scr_out <<
" " << line[k] << std::endl;
2879 O2SCL_ERR2(
"Invalid row for incrementing multiplier in ",
2889 this->
scr_out <<
"mcmc: Updating mult of row "
2890 << walker_accept_rows[windex]
2891 <<
" from " << mult_old <<
" to "
2892 << mult_old+1.0 << std::endl;
2921 if (fabs(
table->
get(
"mult",i))<0.1) {
2939 std::vector<double> &ac_coeff_avg,
2940 int loc_verbose=0) {
2942 ac_coeff_avg.resize(0);
2946 std::vector<std::vector<double> > ac_coeff_temp(n_tot);
2949 for(
size_t k=0;k<this->
n_walk;k++) {
2950 size_t tindex=j*this->n_walk+k;
2951 std::vector<double> vec, mult;
2953 if (fabs(
table->
get(
"walker",ell)-k)<0.1 &&
2954 fabs(
table->
get(
"thread",ell)-j)<0.1 &&
2956 mult.push_back(
table->
get(
"mult",ell));
2957 vec.push_back(
table->
get(icol,ell));
2960 if (mult.size()>1) {
2962 (vec.size(),vec,mult,ac_coeff_temp[tindex]);
2963 if (ac_coeff_temp[tindex].size()>max_size) {
2964 max_size=ac_coeff_temp[tindex].size();
2967 if (loc_verbose>0) {
2968 std::cout <<
"column index, thread, walker, length, mult. sum: "
2969 << icol <<
" " << j <<
" " << k <<
" "
2970 << mult.size() <<
" "
2971 << o2scl::vector_sum<std::vector<double>,
double>
2972 (mult.size(),mult) <<
" " << max_size << std::endl;
2976 for(
size_t j=0;j<max_size;j++) {
2979 for(
size_t k=0;k<n_tot;k++) {
2980 if (j<ac_coeff_temp[k].size()) {
2981 ac_sum+=ac_coeff_temp[k][j];
2986 ac_coeff_avg.push_back(ac_sum/((
double)ac_count));
2998 std::shared_ptr<o2scl::table_units<> > table2=
3002 for(
size_t j=0;j<this->
n_walk;j++) {
3003 std::string func=std::string(
"abs(walker-")+
o2scl::szttos(j)+
3027 for(
size_t it=0;it<this->
n_threads;it++) {
3031 O2SCL_ERR2(
"Cannot reblock. Not enough data in ",
3034 size_t n_block=n/n_blocks;
3036 for(
size_t j=0;j<n_blocks;j++) {
3039 for(
size_t i=0;i<m;i++) {
3042 for(
size_t k=j*n_block;k<(j+1)*n_block;k++) {
3043 mult+=(*table)[
"mult"][k];
3044 for(
size_t i=1;i<m;i++) {
3045 dat[i]+=(*table)[i][k]*(*table)[
"mult"][k];
3049 for(
size_t i=1;i<m;i++) {
3070 template<
class func_t,
class fill_t,
class data_t,
class vec_t=ubvector>
3071 class mcmc_para_cli :
public mcmc_para_table<func_t,fill_t,
3111 std::vector<std::string>
cl_args;
3152 p_file_update_iters.
help=((std::string)
"Number of MCMC successes ")+
3153 "between file updates (default 0 for no file updates).";
3154 cl.
par_list.insert(std::make_pair(
"file_update_iters",
3155 &p_file_update_iters));
3158 p_file_update_time.
help=((std::string)
"Time ")+
3159 "between file updates (default 0.0 for no file updates).";
3160 cl.
par_list.insert(std::make_pair(
"file_update_time",
3161 &p_file_update_time));
3172 p_max_time.
help=((std::string)
"Maximum run time in seconds ")+
3173 "(default 86400 sec or 1 day).";
3174 cl.
par_list.insert(std::make_pair(
"max_time",&p_max_time));
3177 p_max_iters.
help=((std::string)
"If non-zero, limit the number of ")+
3178 "iterations to be less than the specified number (default zero).";
3179 cl.
par_list.insert(std::make_pair(
"max_iters",&p_max_iters));
3182 p_prefix.
help=
"Output file prefix (default 'mcmc\').";
3183 cl.
par_list.insert(std::make_pair(
"prefix",&p_prefix));
3193 p_step_fac.
help=((std::string)
"MCMC step factor. The step size for ")+
3194 "each variable is taken as the difference between the high and low "+
3195 "limits divided by this factor (default 10.0). This factor can "+
3196 "be increased if the acceptance rate is too small, but care must "+
3197 "be taken, e.g. if the conditional probability is multimodal. If "+
3198 "this step size is smaller than 1.0, it is reset to 1.0 .";
3199 cl.
par_list.insert(std::make_pair(
"step_fac",&p_step_fac));
3202 p_n_warm_up.
help=((std::string)
"Minimum number of warm up iterations ")+
3204 cl.
par_list.insert(std::make_pair(
"n_warm_up",&p_n_warm_up));
3207 p_verbose.
help=((std::string)
"MCMC verbosity parameter ")+
3209 cl.
par_list.insert(std::make_pair(
"mcmc_verbose",&p_verbose));
3212 p_max_bad_steps.
help=((std::string)
"Maximum number of bad steps ")+
3214 cl.
par_list.insert(std::make_pair(
"max_bad_steps",&p_max_bad_steps));
3217 p_n_walk.
help=((std::string)
"Number of walkers ")+
3219 cl.
par_list.insert(std::make_pair(
"n_walk",&p_n_walk));
3222 p_user_seed.
help=((std::string)
"Seed for multiplier for random ")+
3223 "number generator. If zero is given (the default), then mcmc() "+
3224 "uses time(0) to generate a random seed.";
3225 cl.
par_list.insert(std::make_pair(
"user_seed",&p_user_seed));
3228 p_aff_inv.
help=((std::string)
"If true, then use affine-invariant ")+
3229 "sampling (default false).";
3230 cl.
par_list.insert(std::make_pair(
"aff_inv",&p_aff_inv));
3233 p_table_sequence.
help=((std::string)
"If true, then ensure equal spacing ")+
3234 "between walkers (default true).";
3235 cl.
par_list.insert(std::make_pair(
"table_sequence",&p_table_sequence));
3238 p_store_rejects.
help=((std::string)
"If true, then store MCMC rejections ")+
3240 cl.
par_list.insert(std::make_pair(
"store_rejects",&p_store_rejects));
std::vector< ubvector > initial_points
Initial points in parameter space.
std::vector< data_t > data_arr
Data array.
void set_proposal_ptrs(prob_vec_t &pv)
Set pointers to proposal distributions.
std::vector< size_t > n_accept
The number of Metropolis steps which were accepted in each independent chain (summed over all walkers...
Integer parameter for o2scl::cli.
size_t file_update_iters
Iterations between file updates (default 0 for no file updates)
size_t get_nlines() const
Return the number of lines.
void set_row(size_t row, size_vec_t &v)
Set an entire row of data.
static const int mcmc_skip
Integer to indicate rejection.
void set_nlines(size_t il)
Set the number of lines.
virtual void clear()
Clear everything.
void copy_rows(std::string func, table< vec2_t > &dest)
Copy all rows matching a particular condition to a new table.
std::string dtos(const fp_t &x, int prec=6, bool auto_prec=false)
Convert a double to a string.
void sets(std::string name, std::string s)
Set a string named name to value s.
virtual void mcmc_cleanup()
Perform cleanup after an MCMC simulation.
std::vector< std::string > col_units
Column units.
virtual void setup_cli(cli &cl)
Set up the 'cli' object.
void new_column(std::string head)
Add a new column owned by the table table .
@ exc_efailed
generic failure
std::vector< vec_t > current
Current points in parameter space for each walker and each OpenMP thread.
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.
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.
size_t n_params
Number of parameters.
bool first_write
If true, the HDF5 I/O initial info has been written to the file (set by mcmc() )
double step_fac
Stepsize factor (default 10.0)
int mpi_size
The MPI number of processors.
int verbose
Output control (default 0)
std::vector< std::string > col_names
Column names.
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
std::string prefix
Prefix for output filenames (default "mcmc")
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.
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
virtual void set_names_units(std::vector< std::string > names, std::vector< std::string > units)
Set the table names and units.
virtual int mcmc(size_t n_params_local, vec_t &low, vec_t &high, std::vector< func_t > &func, std::vector< fill_t > &fill)
Perform an MCMC simulation.
std::vector< o2scl::prob_cond_mdim< vec_t > * > prop_dist
Pointer to proposal distribution for each thread.
virtual void write_files(bool sync_write=false)
Write MCMC tables to files.
void get_chain_sizes(std::vector< size_t > &chain_sizes)
Determine the chain sizes.
size_t max_iters
If non-zero, the maximum number of MCMC iterations (default 0)
int sets_vec(std::string name, const std::vector< std::string > &s)
Set a vector of strings named name.
size_t table_prealloc
Number of rows to allocate for the table before the MCMC run.
@ exc_eunimpl
requested feature not (yet) implemented
bool always_accept
If true, accept all steps.
void set(std::string scol, size_t row, double val)
Set row row of column named col to value val . .
bool warm_up
If true, we are in the warm up phase.
size_t n_chains_per_rank
Number of fully independent chains in each MPI rank.
void close()
Close the file.
std::vector< int > walker_reject_rows
For each walker, record the last row in the table which corresponds to an reject.
bool aff_inv
If true, use affine-invariant Monte Carlo.
A generic MCMC simulation class.
void hdf_input(hdf_file &hf, o2scl::table< vec_t > &t, std::string name)
Input a o2scl::table object from a hdf_file.
virtual void reorder_table()
Reorder the table by thread and walker index.
String parameter for o2scl::cli.
String parameter for o2scl::cli.
void set_maxlines(size_t llines)
Manually set the maximum number of lines.
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.
Data table table class with units.
bool meas_for_initial
If true, call the measurement function for the initial point.
virtual void ac_coeffs(size_t icol, std::vector< double > &ac_coeff_avg, int loc_verbose=0)
Compute autocorrelation coefficient for column with index icol averaging over all walkers and all thr...
int open(std::string fname, bool write_access=false, bool err_on_fail=true)
Open a file named fname.
std::map< std::string, parameter *, std::less< std::string > > par_list
Parameter list.
int mpi_rank
The MPI processor rank.
int setd_vec_copy(std::string name, const vec_t &v)
Set vector dataset named name with v.
virtual void file_header(o2scl_hdf::hdf_file &hf)
Initial write to HDF5 file.
virtual void best_point(vec_t &best, double w_best, data_t &dat)
Function to run for the best point.
virtual void unset_proposal()
Go back to random-walk Metropolis with a uniform distribution.
void set_proposal(prob_vec_t &pv)
Set the proposal distribution.
@ exc_einval
invalid argument supplied by user
mcmc_para_base< func_t, internal_measure_t, data_t, vec_t > parent_t
Type of parent class.
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.
bool prev_read
If true, previous results have been read.
Integer parameter for o2scl::cli.
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.
size_t n_threads
Number of OpenMP threads.
size_t last_write_iters
Total number of MCMC acceptances over all threads at last file write() (default 0)
virtual void mcmc_cleanup()
Cleanup after the MCMC.
size_t n_warm_up
Number of warm up steps (successful steps not iterations) (default 0)
Configurable command-line interface.
void set_szt(std::string name, size_t u)
Set an unsigned integer named name to value u.
void reblock(size_t n_blocks)
Reaverage the data into blocks of a fixed size in order to avoid autocorrelations.
vec_t low_copy
A copy of the lower limits for HDF5 output.
Double parameter for o2scl::cli.
std::shared_ptr< o2scl::table_units<> > table
Main data table for Markov chain.
virtual void file_header(o2scl_hdf::hdf_file &hf)
Initial write to HDF5 file.
size_t get_ncolumns() const
Return the number of columns.
bool err_nonconv
If true, call the error handler if msolve() or msolve_de() does not converge (default true)
double last_write_time
Time at last file write() (default 0.0)
virtual int fill_line(const vec_t &pars, double log_weight, std::vector< double > &line, data_t &dat, size_t i_walker, fill_t &fill)
Fill line with data for insertion into the table.
std::shared_ptr< o2scl::table_units<> > get_table()
Get the output table.
std::string itos(int x)
Convert an integer to a string.
size_t max_bad_steps
Maximum number of failed steps when generating initial points with affine-invariant sampling (default...
std::vector< int > walker_accept_rows
For each walker, record the last row in the table which corresponds to an accept.
virtual int mcmc_init()
Initializations before the MCMC.
int table_io_chunk
The number of tables to combine before I/O (default 1)
int user_seed
If non-zero, use as the seed for the random number generator (default 0)
void seti(std::string name, int i)
Set an integer named name to value i.
void vector_autocorr_vector_mult(size_t n2, const vec_t &data, const vec2_t &mult, resize_vec_t &ac_vec)
Construct an autocorrelation vector using a multiplier using the first n2 elements of vectors data an...
std::ofstream scr_out
The screen output file.
bool table_sequence
If true, ensure sure walkers and OpenMP threads are written to the table with equal spacing between r...
Store data in an O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$sc...
void open_or_create(std::string fname)
Open a file named fname or create if it doesn't already exist.
void set_table(std::shared_ptr< o2scl::table_units<> > &t)
Set the output table.
double max_time
Time in seconds (default is 0)
void setd(std::string name, double d)
Set a double named name to value d.
std::vector< std::string > cl_args
The arguments sent to the command-line.
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)
@ exc_esanity
sanity check failed - shouldn't happen
o2scl::mcmc_para_table< func_t, fill_t, data_t, vec_t > parent_t
The parent typedef.
std::vector< rng_gsl > rg
Random number generators.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
std::string * str
Parameter.
size_t n_walk_per_thread
Number of walkers per thread (default 0 to set equal to n_walk)
bool is_column(std::string scol) const
Return true if scol is a column in the current table table .
virtual void critical_extra(size_t i_thread)
Additional code to execute inside the OpenMP critical section.
size_t n_walk
Number of walkers (per openMP thread) for affine-invariant MC or 1 otherwise (default 1)
bool store_rejects
If true, store MCMC rejections in the table.
std::vector< std::vector< size_t > > ret_value_counts
Return value counters, one vector independent chain.
std::vector< size_t > n_reject
The number of Metropolis steps which were rejected in each independent chain (summed over all walkers...
std::vector< bool > switch_arr
Data switch array for each walker and each OpenMP thread.
double mpi_start_time
The MPI starting time (defaults to 0.0)
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.
virtual int mcmc_init()
MCMC initialization function.
std::vector< double > step_vec
Optionally specify step sizes for each parameter.
double get(std::string scol, size_t row) const
Get value from row row of column named col. .
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.
A generic MCMC simulation class writing data to a o2scl::table_units object.
int set_szt_vec(std::string name, const std::vector< size_t > &v)
Set vector dataset named name with v.
bool pd_mode
If true, then use the user-specified proposal distribution.
std::string szttos(size_t x)
Convert a size_t to a string.
double file_update_time
Time between file updates (default 0.0 for no file updates)
double ai_initial_step
Initial step fraction for affine-invariance sampling walkers (default 0.1)
std::function< int(const vec_t &, double, size_t, int, bool, data_t &)> internal_measure_t
Measurement functor type for the parent.
vec_t high_copy
A copy of the upper limits for HDF5 output.
std::string help
Help description.
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.
static const int mcmc_done
Integer to indicate completion.
std::vector< size_t > curr_walker
Index of the current walker.
std::string get_column_name(size_t icol) const
Returns the name of column col .
Documentation generated with Doxygen. Provided under the
GNU Free Documentation License (see License Information).