Go to the documentation of this file.
27 #ifndef O2SCL_MCMC_PARA_OLD_H
28 #define O2SCL_MCMC_PARA_OLD_H
40 #include <boost/numeric/ublas/vector.hpp>
42 #include <o2scl/uniform_grid.h>
43 #include <o2scl/table3d.h>
44 #include <o2scl/hdf_file.h>
45 #include <o2scl/exception.h>
46 #include <o2scl/prob_dens_func.h>
47 #include <o2scl/vector.h>
48 #include <o2scl/multi_funct.h>
49 #include <o2scl/vec_stats.h>
50 #include <o2scl/cli.h>
133 template<
class func_t,
class measure_t,
151 std::vector<rng_gsl>
rg;
196 O2SCL_ERR2(
"Using a proposal distribution with affine-invariant ",
197 "sampling not implemented in mcmc_para_old::mcmc_init().",
202 std::cout <<
"Prefix is: " <<
prefix << std::endl;
209 scr_out.setf(std::ios::scientific);
223 <<
" reject=" <<
n_reject[it] << std::endl;
232 virtual void best_point(vec_t &best,
double w_best, data_t &dat) {
400 MPI_Comm_rank(MPI_COMM_WORLD,&this->mpi_rank);
401 MPI_Comm_size(MPI_COMM_WORLD,&this->mpi_size);
429 virtual int mcmc(
size_t n_params, vec_t &low, vec_t &high,
430 std::vector<func_t> &func, std::vector<measure_t> &meas) {
436 if (func.size()==0 || meas.size()==0) {
437 O2SCL_ERR2(
"Size of 'func' or 'meas' array is zero in ",
442 std::cout <<
"mcmc_para_old::mcmc(): Not enough functions for "
443 <<
n_threads <<
" threads. Setting n_threads to "
444 << func.size() <<
"." << std::endl;
450 std::cout <<
"mcmc_para_old::mcmc(): Not enough measurement objects for "
451 <<
n_threads <<
" threads. Setting n_threads to "
452 << meas.size() <<
"." << std::endl;
471 for(
size_t k=0;k<n_params;k++) {
480 for(
size_t ipar=0;ipar<n_params;ipar++) {
488 ") in mcmc_base::mcmc().").c_str(),
510 std::cerr.setf(std::ios::scientific);
511 std::cerr << i <<
" ";
513 std::cerr << j <<
" ";
528 std::cout <<
"mcmc_para_old::mcmc(): "
529 <<
n_threads <<
" threads were requested but the "
530 <<
"-DO2SCL_OPENMP flag was not used during "
531 <<
"compilation. Setting n_threads to 1."
543 std::cout <<
"mcmc_para_old::mcmc(): Requested only 1 walker, "
544 <<
"forcing 2 walkers." << std::endl;
547 #ifdef O2SCL_NEVER_DEFINED
554 O2SCL_ERR2(
"More walkers per thread than total walkers ",
559 std::cout <<
"mcmc_para_old::mcmc(): Could not evenly "
560 <<
"organize threads and walkers." << std::endl;
561 std::cout <<
"n_threads: " <<
n_threads << std::endl;
562 std::cout <<
"n_walk: " <<
n_walk << std::endl;
563 std::cout <<
"n_walk_per_thread: "
565 std::cout <<
"n_chains_per_rank: "
577 std::cout <<
"mcmc_para_old::mcmc(): Requested negative or zero "
578 <<
"step_fac with aff_inv=true.\nSetting to 2.0."
582 std::cout <<
"mcmc_para_old::mcmc(): Requested negative or zero "
583 <<
"step_fac. Setting to 10.0." << std::endl;
590 unsigned long int seed=time(0);
591 if (this->user_seed!=0) {
596 rg[it].set_seed(seed);
618 std::vector<double> w_current(ssize);
619 for(
size_t i=0;i<ssize;i++) {
636 next[it].resize(n_params);
641 vec_t best(n_params);
646 std::vector<bool> mcmc_done_flag(
n_threads);
648 mcmc_done_flag[it]=
false;
659 O2SCL_ERR(
"Function mcmc_init() failed in mcmc_base::mcmc().",
670 scr_out <<
"mcmc: Affine-invariant step, n_params="
671 << n_params <<
", n_walk=" <<
n_walk
674 <<
", n_threads=" <<
n_threads <<
", rank="
678 scr_out <<
"mcmc: With proposal distribution, n_params="
679 << n_params <<
", n_threads=" <<
n_threads <<
", rank="
683 scr_out <<
"mcmc: Random-walk w/uniform dist., n_params="
684 << n_params <<
", n_threads=" <<
n_threads <<
", rank="
697 #pragma omp parallel default(shared)
724 for(
size_t ipar=0;ipar<n_params;ipar++) {
729 func_ret[it]=func[it](n_params,
current[sindex],
730 w_current[sindex],
data_arr[sindex]);
733 mcmc_done_flag[it]=
true;
738 if (func_ret[it]>=0 &&
743 meas_ret[it]=meas[it](
current[sindex],w_current[sindex],
750 mcmc_done_flag[it]=
true;
759 while (!done && !mcmc_done_flag[it]) {
762 for(
size_t ipar=0;ipar<n_params;ipar++) {
766 (
rg[it].random()*2.0-1.0)*
768 }
while (
current[sindex][ipar]>high[ipar] ||
769 current[sindex][ipar]<low[ipar]);
773 func_ret[it]=func[it](n_params,
current[sindex],
774 w_current[sindex],
data_arr[sindex]);
782 mcmc_done_flag[it]=
true;
787 if (func_ret[it]>=0 &&
793 meas_ret[it]=meas[it](
current[sindex],w_current[sindex],
800 mcmc_done_flag[it]=
true;
804 std::string err=((std::string)
"In loop with thread ")+
807 " in mcmc_para_old_base::mcmc().";
818 bool stop_early=
false;
820 if (mcmc_done_flag[it]==
true) {
823 <<
"): Returned mcmc_done "
824 <<
"(initial; ai)." << std::endl;
841 if (w_current[sindex]>w_current[0]) {
843 w_best=w_current[sindex];
859 << w_current[sindex] <<
" (initial; ai)" << std::endl;
873 #pragma omp parallel default(shared)
888 for(
size_t ipar=0;ipar<n_params;ipar++) {
895 func_ret[it]=func[it](n_params,
current[it],w_current[it],
899 func_ret[it]=func_ret[it % ip_size];
900 w_current[it]=w_current[it % ip_size];
903 for(
size_t j=0;j<
data_arr.size();j++) {
918 <<
"): Initial point returned mcmc_done."
926 O2SCL_ERR((((std::string)
"Initial weight from thread ")+
928 " vanished in mcmc_para_old_base::mcmc().").c_str(),
939 #pragma omp parallel default(shared)
951 func_ret[it]=func_ret[it % ip_size];
953 w_current[it]=w_current[it % ip_size];
958 if (func_ret[it]>=0 &&
964 meas_ret[it]=meas[it](
current[it],w_current[it],0,
970 mcmc_done_flag[it]=
true;
981 bool stop_early=
false;
983 if (mcmc_done_flag[it]==
true) {
986 <<
"): Returned mcmc_done "
987 <<
"(initial)." << std::endl;
1002 if (w_current[it]<w_best) {
1004 w_best=w_current[it];
1013 scr_out << w_current[it] <<
" ";
1015 scr_out <<
" (initial)" << std::endl;
1030 std::cout <<
"Press a key and type enter to continue. ";
1039 bool main_done=
false;
1040 size_t mcmc_iters=0;
1042 while (!main_done) {
1056 #pragma omp parallel default(shared)
1071 ij=((size_t)(
rg[it].random()*((double)
n_walk)));
1075 double p=
rg[it].random();
1077 smove_z[it]=(1.0-2.0*p+2.0*a*p+p*p-2.0*a*p*p+a*a*p*p)/a;
1080 for(
size_t i=0;i<n_params;i++) {
1091 if (!std::isfinite(q_prop[it])) {
1092 std::cout <<
"Here: " << q_prop[it] << std::endl;
1093 O2SCL_ERR2(
"Proposal distribution not finite in ",
1100 for(
size_t k=0;k<n_params;k++) {
1101 next[it][k]=
current[it][k]+(
rg[it].random()*2.0-1.0)*
1113 for(
size_t k=0;k<n_params;k++) {
1114 if (next[it][k]<low[k] || next[it][k]>high[k]) {
1117 if (next[it][k]<low[k]) {
1118 std::cout <<
"mcmc (" << it <<
","
1119 <<
mpi_rank <<
"): Parameter with index " << k
1120 <<
" and value " << next[it][k]
1121 <<
" smaller than limit " << low[k] << std::endl;
1122 scr_out <<
"mcmc (" << it <<
","
1123 <<
mpi_rank <<
"): Parameter with index " << k
1124 <<
" and value " << next[it][k]
1125 <<
" smaller than limit " << low[k] << std::endl;
1127 std::cout <<
"mcmc (" << it <<
"," <<
mpi_rank
1128 <<
"): Parameter with index " << k
1129 <<
" and value " << next[it][k]
1130 <<
" larger than limit " << high[k] << std::endl;
1132 <<
"): Parameter with index " << k
1133 <<
" and value " << next[it][k]
1134 <<
" larger than limit " << high[k] << std::endl;
1144 func_ret[it]=func[it](n_params,next[it],w_next[it],
1148 func_ret[it]=func[it](n_params,next[it],w_next[it],
1152 mcmc_done_flag[it]=
true;
1154 if (func_ret[it]>=0 &&
1174 scr_out <<
"it: " << it <<
" q_prop[it]: "
1175 << q_prop[it] << std::endl;
1179 <<
"): Returned mcmc_done."
1183 <<
"): Parameter(s) out of range: " << std::endl;
1184 scr_out.setf(std::ios::showpos);
1185 for(
size_t k=0;k<n_params;k++) {
1186 scr_out << k <<
" " << low[k] <<
" "
1187 << next[it][k] <<
" " << high[k];
1188 if (next[it][k]<low[k] || next[it][k]>high[k]) {
1193 scr_out.unsetf(std::ios::showpos);
1198 <<
"): Function returned failure "
1199 << func_ret[it] <<
" at point ";
1200 for(
size_t k=0;k<n_params;k++) {
1201 scr_out << next[it][k] <<
" ";
1214 #pragma omp parallel default(shared)
1232 double r=
rg[it].random();
1235 double ai_ratio=pow(smove_z[it],((
double)n_params)-1.0)*
1236 exp(w_next[it]-w_current[sindex]);
1241 if (r<exp(w_next[it]-w_current[sindex]+q_prop[it])) {
1246 if (r<exp(w_next[it]-w_current[sindex])) {
1261 meas_ret[it]=meas[it](next[it],w_next[it],
1265 meas_ret[it]=meas[it](next[it],w_next[it],
1273 w_current[sindex]=w_next[it];
1284 meas_ret[it]=meas[it](next[it],w_next[it],
1288 meas_ret[it]=meas[it](next[it],w_next[it],
1311 scr_out.width((
int)(1.0+log10((
double)(n_params-1))));
1312 scr_out << mcmc_iters <<
" i_walk: "
1314 << w_current[sindex] << std::endl;
1341 O2SCL_ERR((((std::string)
"Measurement function returned ")+
1343 " in mcmc_para_old_base::mcmc().").c_str(),
1352 if (main_done==
false) {
1364 scr_out <<
"mcmc: Finished warmup." << std::endl;
1371 std::cout <<
"Press a key and type enter to continue. ";
1380 scr_out <<
"mcmc: Stopping because number of iterations ("
1381 << mcmc_iters <<
") equal to max_iters (" <<
max_iters
1382 <<
")." << std::endl;
1387 if (main_done==
false) {
1396 scr_out <<
"mcmc: Stopping because elapsed (" << elapsed
1397 <<
") > max_time (" <<
max_time <<
")."
1421 virtual int mcmc(
size_t n_params, vec_t &low, vec_t &high,
1422 func_t &func, measure_t &meas) {
1435 return mcmc(n_params,low,high,func,meas);
1448 for(
size_t i=0;i<pv.size();i++) {
1464 for(
size_t i=0;i<pv.size();i++) {
1528 template<
class func_t,
class fill_t,
class data_t,
class vec_t=ubvector>
1530 std::function<int(const vec_t &,double,size_t,int,bool,data_t &)>,
1536 typedef std::function<int(
const vec_t &,
double,
size_t,
int,
bool,data_t &)>
1552 std::shared_ptr<o2scl::table_units<> >
table;
1576 for(
size_t i=0;i<
col_names.size();i++) {
1593 std::cout <<
"mcmc: Table column names and units: " << std::endl;
1606 O2SCL_ERR2(
"Table does not have the correct number of columns ",
1614 O2SCL_ERR2(
"Table does not have the correct internal columns ",
1618 O2SCL_ERR2(
"Array walker_accept_rows does not have correct size ",
1622 O2SCL_ERR2(
"Array walker_reject_rows does not have correct size ",
1645 std::vector<double> &line, data_t &dat,
1646 size_t i_walker, fill_t &fill) {
1649 size_t i_thread=omp_get_thread_num();
1657 line.push_back(i_thread);
1659 line.push_back(i_walker);
1661 line.push_back(1.0);
1662 line.push_back(log_weight);
1663 for(
size_t i=0;i<pars.size();i++) {
1664 line.push_back(pars[i]);
1666 int tempi=fill(pars,log_weight,line,dat);
1743 this->
scr_out <<
"mcmc: Start write_files(). mpi_rank: "
1745 << this->
mpi_size <<
" table_io_chunk: "
1746 << table_io_chunk << std::endl;
1749 std::vector<o2scl::table_units<> > tab_arr;
1750 bool rank_sent=
false;
1754 if (this->
mpi_rank%table_io_chunk==0) {
1760 tab_arr.push_back(t);
1761 o2scl_table_mpi_recv(child,tab_arr[tab_arr.size()-1]);
1767 o2scl_table_mpi_send(*
table,parent);
1776 int tag=0, buffer=0;
1777 if (sync_write && this->
mpi_size>1 &&
1779 MPI_Recv(&buffer,1,MPI_INT,this->
mpi_rank-table_io_chunk,
1780 tag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
1799 hf.
set_szt(
"n_params",this->n_params);
1807 hf.
seti(
"store_rejects",this->store_rejects);
1808 hf.
seti(
"table_sequence",this->table_sequence);
1818 this->ret_value_counts[0].size(),
1819 this->ret_value_counts);
1821 this->initial_points[0].size(),
1822 this->initial_points);
1824 hf.
seti(
"n_tables",tab_arr.size()+1);
1825 if (rank_sent==
false) {
1826 hdf_output(hf,*
table,
"markov_chain_0");
1828 for(
size_t i=0;i<tab_arr.size();i++) {
1829 std::string name=((std::string)
"markov_chain_")+
szttos(i+1);
1830 hdf_output(hf,tab_arr[i],name);
1836 if (sync_write && this->
mpi_size>1 &&
1838 MPI_Send(&buffer,1,MPI_INT,this->
mpi_rank+table_io_chunk,
1839 tag,MPI_COMM_WORLD);
1844 this->
scr_out <<
"mcmc: Done write_files()." << std::endl;
1865 std::vector<std::string> units) {
1866 if (names.size()!=units.size()) {
1867 O2SCL_ERR2(
"Size of names and units arrays don't match in ",
1894 int tag=0, buffer=0;
1896 MPI_Recv(&buffer,1,MPI_INT,this->
mpi_rank-1,
1897 tag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
1909 MPI_Send(&buffer,1,MPI_INT,this->
mpi_rank+1,
1910 tag,MPI_COMM_WORLD);
1918 std::cout <<
"Initial points: Finding last " << n_points
1919 <<
" points from file named "
1920 << fname <<
" ." << std::endl;
1925 for(
size_t it=0;it<this->
n_threads;it++) {
1926 for(
size_t iw=0;iw<this->
n_walk;iw++) {
1929 size_t windex=it*this->n_walk+iw;
1932 for(
int row=tip.get_nlines()-1;row>=0 && found==
false;row--) {
1933 if (tip.get(
"walker",row)==iw &&
1934 tip.get(
"thread",row)==it &&
1935 tip.get(
"mult",row)>0.5) {
1939 std::cout <<
"Function initial_point_file_last():\n\tit: "
1941 <<
" iw: " << iw <<
" row: "
1942 << row <<
" log_wgt: " << tip.get(
"log_wgt",row)
1947 for(
size_t ip=0;ip<n_param_loc;ip++) {
1953 O2SCL_ERR(
"Function initial_points_file_last() failed.",
1980 int tag=0, buffer=0;
1982 MPI_Recv(&buffer,1,MPI_INT,this->
mpi_rank-1,
1983 tag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
1995 MPI_Send(&buffer,1,MPI_INT,this->
mpi_rank+1,
1996 tag,MPI_COMM_WORLD);
2004 std::cout <<
"Initial points: Finding last " << n_points
2005 <<
" points from file named "
2006 << fname <<
" ." << std::endl;
2011 size_t nlines=tip.get_nlines();
2012 size_t decrement=nlines/n_points;
2013 if (decrement<1) decrement=1;
2016 for(
size_t k=0;k<n_points;k++) {
2020 std::cout <<
"Function initial_point_file_dist():\n\trow: "
2021 << row <<
" log_wgt: " << tip.get(
"log_wgt",row)
2026 for(
size_t ip=0;ip<n_param_loc;ip++) {
2046 double thresh=1.0e-6,
2054 int tag=0, buffer=0;
2056 MPI_Recv(&buffer,1,MPI_INT,this->
mpi_rank-1,
2057 tag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
2069 MPI_Send(&buffer,1,MPI_INT,this->
mpi_rank+1,
2070 tag,MPI_COMM_WORLD);
2078 std::cout <<
"Initial points: Finding best " << n_points
2079 <<
" unique points from file named "
2080 << fname <<
" ." << std::endl;
2083 typedef std::map<double,int,std::greater<double> > map_t;
2087 for(
size_t k=0;k<tip.get_nlines();k++) {
2088 m.insert(std::make_pair(tip.get(
"log_wgt",k),k));
2098 for(map_t::iterator mit=m.begin();mit!=m.end();mit++) {
2099 map_t::iterator mit2=mit;
2101 if (mit2!=m.end()) {
2102 if (fabs(mit->first-mit2->first)<thresh) {
2104 std::cout <<
"Removing duplicate weights: "
2105 << mit->first <<
" " << mit2->first << std::endl;
2114 }
while (found==
true);
2117 if (m.size()<n_points) {
2118 O2SCL_ERR2(
"Could not find enough points in file in ",
2119 "mcmc_para_old::initial_points_file_best().",
2125 map_t::iterator mit=m.begin();
2126 for(
size_t k=0;k<n_points;k++) {
2127 int row=mit->second;
2129 std::cout <<
"Initial point " << k <<
" at row "
2130 << row <<
" has log_weight= "
2131 << tip.get(
"log_wgt",row) << std::endl;
2134 for(
size_t ip=0;ip<n_param_loc;ip++) {
2150 virtual int mcmc(
size_t n_params_local,
2151 vec_t &low, vec_t &high, std::vector<func_t> &func,
2152 std::vector<fill_t> &fill) {
2170 std::vector<internal_measure_t> meas(this->
n_threads);
2171 for(
size_t it=0;it<this->
n_threads;it++) {
2173 (std::mem_fn<
int(
const vec_t &,
double,
size_t,
int,
bool,
2174 data_t &,
size_t, fill_t &)>
2176 std::placeholders::_2,std::placeholders::_3,
2177 std::placeholders::_4,std::placeholders::_5,
2178 std::placeholders::_6,it,std::ref(fill[it]));
2206 chain_sizes.resize(ntot);
2208 for(
size_t it=0;it<this->
n_threads;it++) {
2209 for(
size_t iw=0;iw<this->
n_walk;iw++) {
2210 size_t ix=it*this->n_walk+iw;
2214 if (
table->
get(
"mult",j)>0.5) chain_sizes[ix]++;
2234 std::string name=
"") {
2247 O2SCL_ERR2(
"Table does not have the correct internal columns ",
2248 "in mcmc_para_old_table::read_prev_results().",
2260 walker_accept_rows.resize(ntot);
2263 for(
size_t j=0;j<ntot;j++) {
2270 size_t i_thread=((size_t)(
table->
get(
"thread",j)+1.0e-12));
2271 size_t i_walker=((size_t)(
table->
get(
"walker",j)+1.0e-12));
2274 size_t windex=i_thread*this->n_walk+i_walker;
2278 }
else if (
table->
get(
"mult",j)<-0.5) {
2287 for(
size_t j=0;j<ntot;j++) {
2295 for(
size_t j=0;j<ntot;j++) {
2297 for(
size_t k=0;k<n_param_loc;k++) {
2302 std::cout <<
"Previous table was read, but initial points not set."
2307 std::cout <<
"mcmc_para_old_table::read_prev_results():" << std::endl;
2308 std::cout <<
" index walker_accept_rows walker_reject_rows"
2310 for(
size_t j=0;j<ntot;j++) {
2313 std::cout << j <<
" ";
2330 virtual int add_line(
const vec_t &pars,
double log_weight,
2331 size_t walker_ix,
int func_ret,
2332 bool mcmc_accept, data_t &dat,
2333 size_t i_thread, fill_t &fill) {
2336 size_t windex=i_thread*this->
n_walk+walker_ix;
2367 #pragma omp critical (o2scl_mcmc_para_old_table_add_line)
2375 fabs(
table->
get(
"mult",next_row))>0.1) {
2390 for(
size_t i=0;i<this->
n_walk;i++) {
2392 table->set(
"thread",istart+j*this->n_walk+i,j);
2393 table->
set(
"walker",istart+j*this->n_walk+i,i);
2394 table->
set(
"mult",istart+j*this->n_walk+i,0.0);
2395 table->
set(
"log_wgt",istart+j*this->n_walk+i,0.0);
2407 std::vector<double> line;
2408 int fret=
fill_line(pars,log_weight,line,dat,walker_ix,fill);
2429 std::cout <<
"line: " << line.size() <<
" columns: "
2432 std::cout << k <<
". ";
2433 if (k<table->get_ncolumns()) {
2437 if (k<line.size()) std::cout << line[k] <<
" ";
2438 std::cout << std::endl;
2440 O2SCL_ERR(
"Table misalignment in mcmc_para_old_table::add_line().",
2449 this->
scr_out <<
"mcmc: Setting data at row " << next_row
2451 for(
size_t k=0;k<line.size();k++) {
2453 this->
scr_out << table->get_column_name(k) <<
" ";
2455 this->
scr_out <<
" " << line[k] << std::endl;
2468 O2SCL_ERR2(
"Invalid row for incrementing multiplier in ",
2478 this->
scr_out <<
"mcmc: Updating mult of row "
2479 << walker_accept_rows[windex]
2480 <<
" from " << mult_old <<
" to "
2481 << mult_old+1.0 << std::endl;
2509 size_t total_iters=0;
2515 this->
scr_out <<
"mcmc: Writing to file. total_iters: "
2516 << total_iters <<
" file_update_iters: "
2533 this->
scr_out <<
"mcmc: Writing to file. elapsed: "
2534 << elapsed <<
" file_update_time: "
2560 if (fabs(
table->
get(
"mult",i))<0.1) {
2576 std::vector<size_t> chain_sizes;
2578 size_t min_size=chain_sizes[0];
2579 for(
size_t i=1;i<chain_sizes.size();i++) {
2580 if (chain_sizes[i]<min_size) min_size=chain_sizes[i];
2582 size_t N_max=min_size/2;
2584 for(
size_t i=0;i<ncols;i++) {
2585 for(
size_t ell=1;ell<N_max;ell++) {
2592 for(
size_t i=0;i<ncols;i++) {
2594 for(
size_t k=0;k<this->
n_walk;k++) {
2595 size_t tindex=j*this->n_walk+k;
2596 for(
size_t ell=1;ell<N_max;ell++) {
2597 const double &x=(*table)[cstart+i][table_row];
2598 double mean=o2scl::vector_mean<const double *>
2599 (chain_sizes[tindex]+1,&x);
2601 <
const double *>(chain_sizes[tindex]+1,&x,ell,mean);
2603 table_row+=chain_sizes[tindex]+1;
2606 for(
size_t ell=1;ell<N_max;ell++) {
2617 size_t N_max=ac_coeffs_cols.size2();
2619 for(
size_t icol=0;icol<ncols;icol++) {
2620 std::vector<double> tau(N_max);
2621 for(
size_t i=5;i<N_max;i++) {
2623 for(
size_t j=0;j<i;j++) {
2624 sum+=ac_coeffs_cols(icol,j);
2627 std::cout << tau[i] <<
" " << ((double)i)/5.0 << std::endl;
2629 std::cout << std::endl;
2639 std::shared_ptr<o2scl::table_units<> > table2=
2643 for(
size_t j=0;j<this->
n_walk;j++) {
2644 std::string func=std::string(
"abs(walker-")+
o2scl::szttos(j)+
2668 for(
size_t it=0;it<this->
n_threads;it++) {
2672 O2SCL_ERR2(
"Cannot reblock. Not enough data in ",
2675 size_t n_block=n/n_blocks;
2677 for(
size_t j=0;j<n_blocks;j++) {
2680 for(
size_t i=0;i<m;i++) {
2683 for(
size_t k=j*n_block;k<(j+1)*n_block;k++) {
2684 mult+=(*table)[
"mult"][k];
2685 for(
size_t i=1;i<m;i++) {
2686 dat[i]+=(*table)[i][k]*(*table)[
"mult"][k];
2690 for(
size_t i=1;i<m;i++) {
2711 template<
class func_t,
class fill_t,
class data_t,
class vec_t=ubvector>
2793 p_file_update_iters.
help=((std::string)
"Number of MCMC successes ")+
2794 "between file updates (default 0 for no file updates).";
2795 cl.
par_list.insert(std::make_pair(
"file_update_iters",
2796 &p_file_update_iters));
2799 p_file_update_time.
help=((std::string)
"Time ")+
2800 "between file updates (default 0.0 for no file updates).";
2801 cl.
par_list.insert(std::make_pair(
"file_update_time",
2802 &p_file_update_time));
2813 p_max_time.
help=((std::string)
"Maximum run time in seconds ")+
2814 "(default 86400 sec or 1 day).";
2815 cl.
par_list.insert(std::make_pair(
"max_time",&p_max_time));
2818 p_max_iters.
help=((std::string)
"If non-zero, limit the number of ")+
2819 "iterations to be less than the specified number (default zero).";
2820 cl.
par_list.insert(std::make_pair(
"max_iters",&p_max_iters));
2823 p_prefix.
help=
"Output file prefix (default 'mcmc\').";
2824 cl.
par_list.insert(std::make_pair(
"prefix",&p_prefix));
2834 p_step_fac.
help=((std::string)
"MCMC step factor. The step size for ")+
2835 "each variable is taken as the difference between the high and low "+
2836 "limits divided by this factor (default 10.0). This factor can "+
2837 "be increased if the acceptance rate is too small, but care must "+
2838 "be taken, e.g. if the conditional probability is multimodal. If "+
2839 "this step size is smaller than 1.0, it is reset to 1.0 .";
2840 cl.
par_list.insert(std::make_pair(
"step_fac",&p_step_fac));
2843 p_n_warm_up.
help=((std::string)
"Minimum number of warm up iterations ")+
2845 cl.
par_list.insert(std::make_pair(
"n_warm_up",&p_n_warm_up));
2848 p_verbose.
help=((std::string)
"MCMC verbosity parameter ")+
2850 cl.
par_list.insert(std::make_pair(
"mcmc_verbose",&p_verbose));
2853 p_max_bad_steps.
help=((std::string)
"Maximum number of bad steps ")+
2855 cl.
par_list.insert(std::make_pair(
"max_bad_steps",&p_max_bad_steps));
2858 p_n_walk.
help=((std::string)
"Number of walkers ")+
2860 cl.
par_list.insert(std::make_pair(
"n_walk",&p_n_walk));
2863 p_user_seed.
help=((std::string)
"Seed for multiplier for random ")+
2864 "number generator. If zero is given (the default), then mcmc() "+
2865 "uses time(0) to generate a random seed.";
2866 cl.
par_list.insert(std::make_pair(
"user_seed",&p_user_seed));
2869 p_aff_inv.
help=((std::string)
"If true, then use affine-invariant ")+
2870 "sampling (default false).";
2871 cl.
par_list.insert(std::make_pair(
"aff_inv",&p_aff_inv));
2874 p_table_sequence.
help=((std::string)
"If true, then ensure equal spacing ")+
2875 "between walkers (default true).";
2876 cl.
par_list.insert(std::make_pair(
"table_sequence",&p_table_sequence));
2879 p_store_rejects.
help=((std::string)
"If true, then store MCMC rejections ")+
2881 cl.
par_list.insert(std::make_pair(
"store_rejects",&p_store_rejects));
bool first_write
If true, the HDF5 I/O initial info has been written to the file (set by mcmc() )
std::vector< std::string > col_names
Column names.
std::shared_ptr< o2scl::table_units<> > table
Main data table for Markov chain.
virtual void ac_lengths(size_t ncols, ubmatrix &ac_coeffs_cols, ubvector &ac_lengths)
Compute autocorrelation lengths.
Integer parameter for o2scl::cli.
size_t n_chains_per_rank
Number of fully independent chains in each MPI rank.
std::vector< size_t > n_accept
The number of Metropolis steps which were accepted in each independent chain (summed over all walkers...
A generic MCMC simulation class.
size_t get_nlines() const
Return the number of lines.
bool table_sequence
If true, ensure sure walkers and OpenMP threads are written to the table with equal spacing between r...
size_t lookup_column(std::string lname) const
Find the index for column named name .
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.
void set_row(size_t row, size_vec_t &v)
Set an entire row of data.
std::vector< std::string > cl_args
The arguments sent to the command-line.
void set_nlines(size_t il)
Set the number of lines.
void copy_rows(std::string func, table< vec2_t > &dest)
Copy all rows matching a particular condition to a new table.
virtual void post_pointmeas()
Function to run after point evaluation and measurement steps.
std::string dtos(const fp_t &x, int prec=6, bool auto_prec=false)
Convert a double to a string.
int mpi_rank
The MPI processor rank.
void sets(std::string name, std::string s)
Set a string named name to value s.
o2scl::mcmc_para_old_table< func_t, fill_t, data_t, vec_t > parent_t
The parent typedef.
std::vector< bool > switch_arr
Data switch array for each walker and each OpenMP thread.
void new_column(std::string head)
Add a new column owned by the table table .
@ exc_efailed
generic failure
virtual void write_files(bool sync_write=false)
Write MCMC tables to files.
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 setup_cli(cli &cl)
Set up the 'cli' object.
void reblock(size_t n_blocks)
Reaverage the data into blocks of a fixed size in order to avoid autocorrelations.
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
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)
size_t n_threads
Number of OpenMP threads.
std::vector< ubvector > initial_points
Initial points in parameter space.
virtual void file_header(o2scl_hdf::hdf_file &hf)
Initial write to HDF5 file.
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.
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
bool always_accept
If true, accept all steps.
A generic MCMC simulation class writing data to a o2scl::table_units object.
size_t file_update_iters
Iterations between file updates (default 0 for no file updates)
double mpi_start_time
The MPI starting time (defaults to 0.0)
bool prev_read
If true, previous results have been read.
static const int mcmc_skip
Integer to indicate rejection.
int sets_vec(std::string name, const std::vector< std::string > &s)
Set a vector of strings named name.
@ exc_eunimpl
requested feature not (yet) implemented
std::vector< size_t > curr_walker
Index of the current walker.
void set(std::string scol, size_t row, double val)
Set row row of column named col to value val . .
double step_fac
Stepsize factor (default 10.0)
int mpi_size
The MPI number of processors.
void close()
Close the file.
std::vector< o2scl::prob_cond_mdim< vec_t > * > prop_dist
Pointer to proposal distribution for each thread.
void hdf_input(hdf_file &hf, o2scl::table< vec_t > &t, std::string name)
Input a o2scl::table object from a hdf_file.
String parameter for o2scl::cli.
mcmc_para_old_base< func_t, internal_measure_t, data_t, vec_t > parent_t
Type of parent class.
std::string prefix
Prefix for output filenames (default "mcmc")
String parameter for o2scl::cli.
std::vector< std::string > col_units
Column units.
virtual int mcmc_init()
Initializations before the MCMC.
Data table table class with units.
void set_proposal(prob_vec_t &pv)
Set the proposal distribution.
std::vector< std::vector< size_t > > ret_value_counts
Return value counters, one vector independent chain.
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 setd_vec_copy(std::string name, const vec_t &v)
Set vector dataset named name with v.
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.
virtual void post_pointmeas()
Function to run after point evaluation and measurement steps.
virtual void mcmc_cleanup()
Perform cleanup after an MCMC simulation.
virtual void file_header(o2scl_hdf::hdf_file &hf)
Initial write to HDF5 file.
virtual void initial_points_file_dist(std::string fname, size_t n_param_loc, size_t offset=5)
Read initial points from file named fname, distributing across the chain if necessary.
std::ofstream scr_out
The screen output file.
bool store_rejects
If true, store MCMC rejections in the table.
virtual void unset_proposal()
Go back to random-walk Metropolis with a uniform distribution.
size_t n_params
Number of parameters.
@ exc_einval
invalid argument supplied by user
void get_chain_sizes(std::vector< size_t > &chain_sizes)
Determine the chain sizes.
int user_seed
If non-zero, use as the seed for the random number generator (default 0)
int table_io_chunk
The number of tables to combine before I/O (default 1)
double last_write_time
Time at last file write() (default 0.0)
vec_t low_copy
A copy of the lower limits for HDF5 output.
bool aff_inv
If true, use affine-invariant Monte Carlo.
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.
virtual void reorder_table()
Reorder the table by thread and walker index.
virtual void mcmc_cleanup()
Cleanup after the MCMC.
vec_t high_copy
A copy of the upper limits for HDF5 output.
std::vector< int > walker_reject_rows
For each walker, record the last row in the table which corresponds to an reject.
Configurable command-line interface.
void set_szt(std::string name, size_t u)
Set an unsigned integer named name to value u.
Double parameter for o2scl::cli.
std::shared_ptr< o2scl::table_units<> > get_table()
Get the output table.
double max_time
Time in seconds (default is 0)
double ai_initial_step
Initial step fraction for affine-invariance sampling walkers (default 0.1)
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)
virtual int mcmc_init()
MCMC initialization function.
size_t max_bad_steps
Maximum number of failed steps when generating initial points with affine-invariant sampling (default...
std::string itos(int x)
Convert an integer to a string.
size_t max_iters
If non-zero, the maximum number of MCMC iterations (default 0)
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.
void seti(std::string name, int i)
Set an integer named name to value i.
size_t n_walk
Number of walkers for affine-invariant MC or 1 otherwise (default 1)
std::vector< vec_t > current
Current points in parameter space for each walker and each OpenMP thread.
std::vector< data_t > data_arr
Data array.
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.
size_t n_warm_up
Number of warm up steps (successful steps not iterations) (default 0)
bool meas_for_initial
If true, call the measurement function for the initial point.
static const int mcmc_done
Integer to indicate completion.
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 setd(std::string name, double d)
Set a double named name to value d.
@ exc_esanity
sanity check failed - shouldn't happen
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
std::string * str
Parameter.
double file_update_time
Time between file updates (default 0.0 for no file updates)
virtual void best_point(vec_t &best, double w_best, data_t &dat)
Function to run for the best point.
bool is_column(std::string scol) const
Return true if scol is a column in the current table table .
double vector_lagk_autocorr(size_t n, const vec_t &data, size_t k, double mean)
Lag-k autocorrelation.
MCMC class with a command-line interface.
std::vector< int > walker_accept_rows
For each walker, record the last row in the table which corresponds to an accept.
std::function< int(const vec_t &, double, size_t, int, bool, data_t &)> internal_measure_t
Measurement functor type for the parent.
bool warm_up
If true, we are in the warm up phase.
virtual void ac_coeffs(size_t ncols, ubmatrix &ac_coeffs)
Compute autocorrelation coefficients.
void set_proposal_ptrs(prob_vec_t &pv)
Set pointers to proposal distributions.
bool pd_mode
If true, then use the user-specified proposal distribution.
std::vector< rng_gsl > rg
Random number generators.
double get(std::string scol, size_t row) const
Get value from row row of column named col. .
size_t n_walk_per_thread
Number of walkers per thread (default 1)
void vector_out(std::ostream &os, size_t n, const vec_t &v, bool endline=false)
Output the first n elements of a vector to a stream, os.
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 set_szt_vec(std::string name, const std::vector< size_t > &v)
Set vector dataset named name with v.
size_t last_write_iters
Total number of MCMC acceptances over all threads at last file write() (default 0)
std::string szttos(size_t x)
Convert a size_t to a string.
virtual void set_names_units(std::vector< std::string > names, std::vector< std::string > units)
Set the table names and units.
std::string help
Help description.
void set_table(std::shared_ptr< o2scl::table_units<> > &t)
Set the output table.
std::vector< size_t > n_reject
The number of Metropolis steps which were rejected in each independent chain (summed over all walkers...
virtual void initial_points_file_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.
std::string get_column_name(size_t icol) const
Returns the name of column col .
int verbose
Output control (default 0)
Documentation generated with Doxygen. Provided under the
GNU Free Documentation License (see License Information).