mcmc_mpi.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2012-2017, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /** \file mcmc_mpi.h
24  \brief Definition of \ref o2scl::mcmc_mpi and \ref o2scl::mcmc_cli
25  classes
26 */
27 #ifndef MCMC_MPI_H
28 #define MCMC_MPI_H
29 
30 #include <iostream>
31 
32 #include <boost/numeric/ublas/vector.hpp>
33 
34 #ifdef O2SCL_MPI
35 #include <mpi.h>
36 #endif
37 
38 #include <o2scl/rng_gsl.h>
39 #include <o2scl/uniform_grid.h>
40 #include <o2scl/table3d.h>
41 #include <o2scl/hdf_file.h>
42 #include <o2scl/hdf_io.h>
43 #include <o2scl/exception.h>
44 #include <o2scl/prob_dens_func.h>
45 #include <o2scl/cholesky.h>
46 #include <o2scl/vector.h>
47 #include <o2scl/multi_funct.h>
48 #include <o2scl/mcmc.h>
49 
50 #ifdef O2SCL_READLINE
51 #include <o2scl/cli_readline.h>
52 #else
53 #include <o2scl/cli.h>
54 #endif
55 
56 namespace o2scl {
57 
58  /** \brief MCMC with MPI and HDF5 table I/O
59 
60  The parent, \ref o2scl::mcmc_table assembles the MCMC data into
61  objects of type \ref o2scl::table_units . This class
62  additionally writes the MCMC data to a HDF5 table and controls
63  run time on several processors via MPI.
64  */
65  template<class func_t, class fill_t, class data_t,
67  public o2scl::mcmc_table<func_t,fill_t,data_t,vec_t> {
68 
69  protected:
70 
74 
75  /// \name MPI properties
76  //@{
77  /// The MPI processor rank
78  int mpi_rank;
79 
80  /// The MPI number of processors
82 
83  /// The MPI starting time
85  //@}
86 
87  /** \brief Error handler for each thread
88  */
90 
91  /** \brief The number of parameters
92  */
93  size_t n_params;
94 
95  /** \brief A copy of the lower limits for HDF5 output
96  */
97  vec_t low_copy;
98 
99  /** \brief A copy of the upper limits for HDF5 output
100  */
101  vec_t high_copy;
102 
103  public:
104 
105  /** \brief Default method for setting the random seed
106  */
107  virtual void set_seed() {
108  // Set RNG seed
109  unsigned long int seed=time(0);
110  if (this->user_seed!=0) {
111  seed=this->user_seed;
112  }
113  seed*=(mpi_rank+1);
114  this->rg.set_seed(seed);
115 
116  return;
117  }
118 
119  /** \brief Perform an MCMC simulation
120  */
121  virtual int mcmc(size_t np, vec_t &init,
122  vec_t &low, vec_t &high, func_t &func,
123  fill_t &fill) {
124 
125  low_copy=low;
126  high_copy=high;
127 
128  // The function mcmc_init() needs to know the number of
129  // parameters, so we store it here for later use
130  n_params=np;
131  return parent_t::mcmc(np,init,low,high,func,fill);
132  }
133 
134  public:
135 
136  /// The screen output file
137  std::ofstream scr_out;
138 
139  /// If true, scr_out has been opened
141 
142  /// \name Command-line settings
143  //@{
144  /** \brief Maximum number of iterations (default 0)
145  */
146  size_t max_iters;
147 
148  /** \brief Time in seconds (default is 86400 seconds or 1 day)
149  */
150  double max_time;
151 
152  /** \brief Output each measurement
153  */
155 
156  /** \brief Prefix for output filenames
157  */
158  std::string prefix;
159 
160  /// If true, output MC accepts and rejects (default true)
161  //@}
162 
163  /// The first point in the parameter space
164  ubvector initial_point;
165 
166  /// The file containing the initial point
167  std::string initial_point_file;
168 
169  /// \name Integer designating how to set the initial point
170  //@{
171  int initial_point_type;
172  static const int fp_unspecified=-1;
173  static const int fp_last=-2;
174  static const int fp_best=-3;
175  //@}
176 
177  /// If true, then \ref first_update() has been called
179 
180  /// \name Command-line settings
181  //@{
182  /** \brief The number of MCMC successes between file updates
183  (default 40)
184  */
186 
187  /** \brief Maximum size of Markov chain (default 10000)
188  */
190  //@}
191 
192  /// Number of complete Markov chain segments
193  size_t chain_index;
194 
195  /** \brief Update files with current table
196  */
197  virtual void update_files() {
198 
200 
201  // Open main update file
202  hf.open_or_create(prefix+"_"+o2scl::itos(mpi_rank)+"_out");
203 
204  // First time, output some initial quantities
205  if (first_file_update==false) {
206  first_update(hf);
207  first_file_update=true;
208  }
209 
210  hf.set_szt("n_accept",this->n_accept);
211  hf.set_szt("n_reject",this->n_reject);
212  hf.set_szt("n_chains",chain_index+1);
213  hf.set_szt_vec("ret_value_counts",this->ret_value_counts);
214 
215  std::string ch_name="markov_chain"+o2scl::szttos(chain_index);
216  hdf_output(hf,*this->tab,ch_name);
217 
218  hf.close();
219 
220  return;
221  }
222 
223  /** \brief Add a measurement to the table
224  */
225  virtual int add_line(const ubvector &pars, double weight,
226  size_t ix, bool new_meas, data_t &dat,
227  fill_t &fill) {
228 
229  if (output_meas) {
230  scr_out << "Line: ";
231  o2scl::vector_out(scr_out,pars);
232  scr_out << " " << weight << " " << ix << " " << new_meas << std::endl;
233  }
234 
236  (pars,weight,ix,new_meas,dat,fill);
237 
238  bool files_updated=false;
239  if (((int)this->tab->get_nlines())==max_chain_size) {
240  scr_out << "Writing files and starting new chain." << std::endl;
241  update_files();
242  chain_index++;
243  this->tab->clear_data();
244  files_updated=true;
245  } else if (this->n_accept>0 && this->n_reject>0 &&
246  (this->n_accept+this->n_reject)%file_update_iters==0) {
247  update_files();
248  files_updated=true;
249  }
250 
251  if (this->max_iters>0 &&
252  (this->n_accept+this->n_reject)==this->max_iters) {
253  if (files_updated==false) {
254  update_files();
255  }
256  if (this->verbose>=1) {
257  std::cout << "mcmc: Stopping because n_accept+n_reject, "
258  << this->n_accept+this->n_reject
259  << ", is equal to max_iters." << std::endl;
260  }
261  return this->mcmc_done;
262  } else {
263  // Determine time elapsed
264 #ifdef O2SCL_MPI
265  double elapsed=MPI_Wtime()-mpi_start_time;
266 #else
267  double elapsed=time(0)-mpi_start_time;
268 #endif
269  if (elapsed>max_time) {
270  if (files_updated==false) {
271  update_files();
272  }
273  if (this->verbose>=1) {
274  std::cout << "mcmc: Stopping because elapsed > max_time."
275  << std::endl;
276  }
277  return this->mcmc_done;
278  }
279  }
280 
281  return 0;
282  }
283 
284  /// \name Customization functions
285  //@{
286  /** \brief User-defined initialization function
287  */
288  virtual int mcmc_init() {
289 
290  if (this->verbose>=2) {
291  std::cout << "Start mcmc_mpi::mcmc_init()." << std::endl;
292  }
293 
295 
296 #ifdef O2SCL_MPI
297  mpi_start_time=MPI_Wtime();
298 #else
299  mpi_start_time=time(0);
300 #endif
301 
302  chain_index=0;
303 
304  if (this->file_opened==false) {
305  // Open main output file
306  this->scr_out.open((this->prefix+"_"+
307  o2scl::itos(this->mpi_rank)+"_scr").c_str());
308  this->scr_out.setf(std::ios::scientific);
309  this->file_opened=true;
310  this->scr_out << "Opened main file in command 'mcmc'." << std::endl;
311  }
312 
313  // Fix file_update_iters if necessary
314  if (file_update_iters<1) {
315  this->scr_out << "Parameter 'file_update_iters' less "
316  << "than 1. Set equal to 1." << std::endl;
317  file_update_iters=1;
318  }
319 
320  if (max_chain_size<1) {
321  O2SCL_ERR("Parameter 'max_chain_size' must be larger than 1.",
323  }
324 
325 #ifdef O2SCL_MPI
326  int buffer=0, tag=0;
327  if (this->mpi_nprocs>1 && this->mpi_rank>0) {
328  MPI_Recv(&buffer,1,MPI_INT,this->mpi_rank-1,tag,MPI_COMM_WORLD,
329  MPI_STATUS_IGNORE);
330  }
331 #endif
332 
333  if (initial_point_file.length()>0) {
334 
335  if (initial_point_type==fp_last) {
336 
337  // Read file
338  this->scr_out << "Reading last point from file '" << initial_point_file
339  << "'." << std::endl;
341  hf.open(initial_point_file);
342 
343  // Read table
344  size_t file_n_chains;
345  hf.get_szt("n_chains",file_n_chains);
346  std::string chain_name=std::string("markov_chain")+
347  o2scl::szttos(file_n_chains-1);
348  o2scl::table_units<> file_tab;
349  hdf_input(hf,file_tab,chain_name);
350  size_t last_line=file_tab.get_nlines()-1;
351 
352  // Get parameters
353  for(size_t i=0;i<n_params;i++) {
354  std::string pname=((std::string)"param_")+this->col_names[i];
355  this->current[0][i]=file_tab.get(pname,last_line);
356  this->scr_out << "Parameter named "
357  << this->col_names[i] << " "
358  << this->current[0][i] << std::endl;
359  }
360 
361  // Finish up
362  this->scr_out << std::endl;
363  hf.close();
364 
365  } else if (initial_point_type==fp_best) {
366 
367  std::vector<double> best_point;
369  hf.open(initial_point_file);
370  hf.getd_vec("best_point",best_point);
371  hf.close();
372  this->scr_out << "Reading best point from file '" << initial_point_file
373  << "'." << std::endl;
374  for(size_t i=0;i<n_params;i++) {
375  this->current[0][i]=best_point[i];
376  this->scr_out << "Parameter " << i << " : "
377  << this->current[0][i] << std::endl;
378  }
379  this->scr_out << "Best weight: "
380  << best_point[n_params] << std::endl;
381  this->scr_out << std::endl;
382 
383  } else {
384 
385  // Read file
386  this->scr_out << "Reading "
387  << initial_point_type << "th point from file '"
388  << initial_point_file
389  << "'." << std::endl;
391  hf.open(initial_point_file);
392 
393  // Read table
394  size_t file_n_chains, row=initial_point_type;
395  hf.get_szt("n_chains",file_n_chains);
396 
397  o2scl::table_units<> file_tab;
398  for(size_t k=0;k<file_n_chains;k++) {
399  std::string chain_name=std::string("markov_chain")+o2scl::szttos(k);
400  hdf_input(hf,file_tab,chain_name);
401  if (file_tab.get_nlines()>row) {
402  k=file_n_chains;
403  } else {
404  row-=file_tab.get_nlines();
405  }
406  }
407  if (row>=file_tab.get_nlines()) {
408  this->scr_out << "Couldn't find point " << initial_point_type
409  << " in file. Using last point." << std::endl;
410  row=file_tab.get_nlines()-1;
411  }
412 
413  // Get parameters
414  for(size_t i=0;i<n_params;i++) {
415  std::string pname=((std::string)"param_")+this->col_names[i];
416  this->current[0][i]=file_tab.get(pname,row);
417  this->scr_out << "Parameter named "
418  << this->col_names[i] << " "
419  << this->current[0][i] << std::endl;
420  }
421 
422  // Finish up
423  this->scr_out << std::endl;
424  hf.close();
425  }
426 
427  } else if (initial_point.size()>0) {
428 
429  this->scr_out << "First point from command-line." << std::endl;
430  for(size_t i=0;i<n_params;i++) {
431  this->current[0][i]=initial_point[i];
432  this->scr_out << this->current[0][i] << std::endl;
433  }
434  this->scr_out << std::endl;
435 
436  } else {
437 
438  this->scr_out << "First point from default." << std::endl;
439 
440  }
441 
442 #ifdef O2SCL_MPI
443  if (this->mpi_nprocs>1 && this->mpi_rank<this->mpi_nprocs-1) {
444  MPI_Send(&buffer,1,MPI_INT,this->mpi_rank+1,tag,MPI_COMM_WORLD);
445  }
446 #endif
447 
448  this->scr_out << "First point: ";
449  o2scl::vector_out(this->scr_out,this->current[0],true);
450 
451  if (this->verbose>=2) {
452  std::cout << "End mcmc_mpi::mcmc_init()." << std::endl;
453  }
454 
455  return 0;
456  };
457 
458  /** \brief Output the best point so far
459  */
460  virtual void best_point(ubvector &best, double w_best) {
461  if (this->file_opened==false) {
462  // Open main output file
463  this->scr_out.open((this->prefix+"_"+
464  o2scl::itos(this->mpi_rank)+"_scr").c_str());
465  this->scr_out.setf(std::ios::scientific);
466  this->file_opened=true;
467  this->scr_out << "Opened main file in function 'best_point()'."
468  << std::endl;
469  }
470  this->scr_out << "Best: ";
471  o2scl::vector_out(this->scr_out,best);
472  this->scr_out << " " << w_best << std::endl;
473  return;
474  }
475 
476  /** \brief Initial write to HDF5 file
477  */
478  virtual void first_update(o2scl_hdf::hdf_file &hf) {
479 
480  hf.sets_vec("param_names",this->col_names);
481 
482  hf.set_szt("n_params",n_params);
483  hf.setd("max_time",this->max_time);
484  hf.seti("user_seed",this->user_seed);
485  hf.seti("n_warm_up",this->n_warm_up);
486  hf.setd("step_fac",this->step_fac);
487  hf.seti("max_iters",this->max_iters);
488  hf.seti("file_update_iters",file_update_iters);
489  hf.seti("output_meas",this->output_meas);
490  hf.seti("initial_point_type",initial_point_type);
491  hf.sets("initial_point_file",initial_point_file);
492  hf.setd_vec_copy("initial_point",initial_point);
493  hf.set_szt("n_chains",chain_index+1);
494  hf.set_szt_vec("ret_value_counts",this->ret_value_counts);
495  hf.setd_vec_copy("low",this->low_copy);
496  hf.setd_vec_copy("high",this->high_copy);
497 
498  return;
499  }
500 
501 #ifdef O2SCL_NEVER_DEFINED
502  /** \brief Choose a Metropolis-Hastings step
503  */
504  virtual int hastings(std::vector<std::string> &sv,
505  bool itive_com) {
506 
507  bool debug=true;
508 
509  if (file_opened==false) {
510  // Open main output file
511  scr_out.open((prefix+"_"+o2scl::itos(mpi_rank)+"_scr").c_str());
512  scr_out.setf(ios::scientific);
513  file_opened=true;
514  scr_out << "Opened main file in command 'hastings'." << endl;
515  }
516 
517  if (sv.size()<2) {
518  cout << "No arguments given to 'hastings'." << endl;
519  return exc_efailed;
520  }
521 
522  if (model_type.length()==0) {
523  cout << "No model selected in 'hastings'." << endl;
524  return exc_efailed;
525  }
526 
527 #ifdef O2SCL_MPI
528  int buffer=0, tag=0;
529  if (mpi_nprocs>1 && mpi_rank>0) {
530  MPI_Recv(&buffer,1,MPI_INT,mpi_rank-1,tag,MPI_COMM_WORLD,
531  MPI_STATUS_IGNORE);
532  }
533 #endif
534 
535  // Read the data file
536  std::string fname=sv[1];
537  scr_out << "Opening file " << fname << " for hastings." << endl;
538  hdf_file hf;
539  hf.open(fname);
540  table_units<> file_tab;
541  hdf_input(hf,file_tab,"markov_chain0");
542  hf.close();
543  scr_out << "Done opening file " << fname << " for hastings." << endl;
544 
545 #ifdef O2SCL_MPI
546  if (mpi_nprocs>1 && mpi_rank<mpi_nprocs-1) {
547  MPI_Send(&buffer,1,MPI_INT,mpi_rank+1,tag,MPI_COMM_WORLD);
548  }
549 #endif
550 
551  // Create a new column equal to mult times weight
552  file_tab.function_column("mult*weight","mwgt");
553 
554  // Remove
555  double max_mwgt=file_tab.max("mwgt");
556  if (debug) scr_out << "lines: " << file_tab.get_nlines() << endl;
557  file_tab.add_constant("max_mwgt",max_mwgt);
558  file_tab.delete_rows("mwgt<0.1*max_mwgt");
559  if (debug) scr_out << "lines: " << file_tab.get_nlines() << endl;
560 
561  // The total number of variables
562  size_t nv=n_params+n_sources;
563  if (debug) {
564  scr_out << n_params << " parameters and " << n_sources << " sources."
565  << endl;
566  }
567  hg_best.resize(nv);
568 
569  // Find the average values
570  for(size_t i=0;i<n_params;i++) {
571  string str_i=((string)"param_")+data_arr[0].modp->param_name(i);
572  hg_best[i]=wvector_mean(file_tab.get_nlines(),file_tab[str_i],
573  file_tab["mwgt"]);
574  }
575  for(size_t i=0;i<n_sources;i++) {
576  string str_i=((string)"Mns_")+source_names[i];
577  hg_best[i+n_params]=wvector_mean(file_tab.get_nlines(),file_tab[str_i],
578  file_tab["mwgt"]);
579  }
580 
581  // Construct the covariance matrix
582  ubmatrix covar(nv,nv);
583  for(size_t i=0;i<n_params;i++) {
584  string str_i=((string)"param_")+data_arr[0].modp->param_name(i);
585  for(size_t j=i;j<n_params;j++) {
586  string str_j=((string)"param_")+data_arr[0].modp->param_name(j);
587  covar(i,j)=wvector_covariance(file_tab.get_nlines(),
588  file_tab[str_i],file_tab[str_j],
589  file_tab["mult"]);
590  if (debug) {
591  scr_out << "Covar: " << i << " " << j << " "
592  << covar(i,j) << endl;
593  }
594  covar(j,i)=covar(i,j);
595  }
596  for(size_t j=0;j<n_sources;j++) {
597  string str_j=((string)"Mns_")+source_names[j];
598  covar(i,j+n_params)=wvector_covariance(file_tab.get_nlines(),
599  file_tab[str_i],file_tab[str_j],
600  file_tab["mult"]);
601  if (debug) {
602  scr_out << "Covar: " << i << " " << j+n_params << " "
603  << covar(i,j+n_params) << endl;
604  }
605  covar(j+n_params,i)=covar(i,j+n_params);
606  }
607  }
608  for(size_t i=0;i<n_sources;i++) {
609  string str_i=((string)"Mns_")+source_names[i];
610  for(size_t j=i;j<n_sources;j++) {
611  string str_j=((string)"Mns_")+source_names[j];
612  covar(i+n_params,j+n_params)=
613  wvector_covariance(file_tab.get_nlines(),
614  file_tab[str_i],file_tab[str_j],
615  file_tab["mult"]);
616  if (debug) {
617  scr_out << "Covar: " << i+n_params << " " << j+n_params << " "
618  << covar(i+n_params,j+n_params) << endl;
619  }
620  covar(j+n_params,i+n_params)=covar(i+n_params,j+n_params);
621  }
622  }
623 
624  // Perform the Cholesky decomposition
625  hg_chol=covar;
626  o2scl_linalg::cholesky_decomp(nv,hg_chol);
627 
628  // Find the inverse
629  hg_covar_inv=hg_chol;
630  o2scl_linalg::cholesky_invert<ubmatrix>(nv,hg_covar_inv);
631 
632  // Force hg_chol to be lower triangular
633  for(size_t i=0;i<nv;i++) {
634  for(size_t j=0;j<nv;j++) {
635  if (i<j) hg_chol(i,j)=0.0;
636  }
637  }
638 
639  // Compute the normalization, weighted by the likelihood function
640  hg_norm=1.0;
641  size_t step=file_tab.get_nlines()/20;
642  if (step<1) step=1;
643  double renorm=0.0;
644  double wgt_sum=0.0;
645  for(size_t i=0;i<file_tab.get_nlines();i+=step) {
646  ubvector e(n_params,n_sources);
647  for(size_t j=0;j<n_params;j++) {
648  string str_j=((string)"param_")+data_arr[0].modp->param_name(j);
649  e.params[j]=file_tab.get(str_j,i);
650  }
651  for(size_t j=0;j<n_sources;j++) {
652  string str_j=((string)"Mns_")+source_names[j];
653  e.mass[j]=file_tab.get(str_j,i);
654  }
655  double wgt=file_tab.get("mult",i)*file_tab.get("weight",i);
656  double rat=wgt/approx_like(e);
657  renorm+=wgt*wgt/approx_like(e);
658  if (debug) {
659  scr_out << wgt << " " << approx_like(e) << " " << rat << endl;
660  }
661  wgt_sum+=wgt;
662  }
663  renorm/=((double)wgt_sum);
664  hg_norm*=renorm;
665  if (debug) {
666  scr_out << "New normalization: " << hg_norm << endl;
667  }
668 
669  step=file_tab.get_nlines()/20;
670  if (step<1) step=1;
671  for(size_t i=0;i<file_tab.get_nlines();i+=step) {
672  ubvector e(n_params,n_sources);
673  for(size_t j=0;j<n_params;j++) {
674  string str_j=((string)"param_")+data_arr[0].modp->param_name(j);
675  e.params[j]=file_tab.get(str_j,i);
676  }
677  for(size_t j=0;j<n_sources;j++) {
678  string str_j=((string)"Mns_")+source_names[j];
679  e.mass[j]=file_tab.get(str_j,i);
680  }
681  double wgt=file_tab.get("mult",i)*file_tab.get("weight",i);
682  double rat=wgt/approx_like(e);
683  if (debug) {
684  scr_out << wgt << " " << approx_like(e) << " " << rat << endl;
685  }
686  }
687  hg_mode=1;
688 
689  return 0;
690  }
691 
692 #endif
693 
694  /** \brief Set the first point
695  */
696  int set_initial_point(std::vector<std::string> &sv,
697  bool itive_com) {
698 
699  if (sv.size()<2) {
700  std::cout << "No arguments given to 'initial-point'." << std::endl;
701  return o2scl::exc_efailed;
702  }
703 
704  if (sv[1]==((std::string)"values")) {
705 
706  initial_point.resize(sv.size()-1);
707  for(size_t i=2;i<sv.size();i++) {
708  initial_point[i-2]=o2scl::function_to_double(sv[i]);
709  }
710  initial_point_type=fp_unspecified;
711 
712  } else if (sv[1]==((std::string)"prefix")) {
713 
714  initial_point_type=fp_last;
715  initial_point_file=sv[2]+((std::string)"_")+
716  o2scl::itos(this->mpi_rank)+"_out";
717 
718  } else if (sv[1]==((std::string)"last")) {
719  initial_point_type=fp_last;
720  initial_point_file=sv[2];
721  } else if (sv[1]==((std::string)"best")) {
722  initial_point_type=fp_best;
723  initial_point_file=sv[2];
724  } else if (sv[1]==((std::string)"N")) {
725  initial_point_type=o2scl::stoi(sv[2]);
726  initial_point_file=sv[3];
727  }
728 
729  return 0;
730  }
731 
732  /** \brief Create an MCMC object with model \c m
733  */
735 
736  // Initial values for MPI paramers
737  mpi_nprocs=1;
738  mpi_rank=0;
739  mpi_start_time=0.0;
740 
741  // Parameters
742  prefix="mcmc";
743  output_meas=true;
744 
745  // Default to 24 hours
746  max_time=3.6e3*24;
747  max_iters=0;
748 
749  // True if scr_out has been opened
750  file_opened=false;
751  first_file_update=false;
752 
753  initial_point_file="";
754  initial_point_type=fp_unspecified;
755 
756  file_update_iters=40;
757  max_chain_size=10000;
758 
759  chain_index=0;
760 
761  // ---------------------------------------
762  // Set error handler for this thread
763 
765 
766  }
767 
768  };
769 
770  /** \brief MCMC class with a command-line interface
771  */
772  template<class func_t, class fill_t, class data_t,
774  public o2scl::mcmc_mpi<func_t,fill_t,data_t,vec_t> {
775 
776  protected:
777 
781 
782 #ifdef O2SCL_READLINE
783  /// Command-line interface
785 #else
786  /// Command-line interface
788 #endif
789 
790  /** \brief The arguments sent to the command-line
791  */
792  std::vector<std::string> cl_args;
793 
794  /// \name Parameter objects for the 'set' command
795  //@{
796  o2scl::cli::parameter_double p_step_fac;
797  o2scl::cli::parameter_size_t p_n_warm_up;
798  o2scl::cli::parameter_int p_user_seed;
799  o2scl::cli::parameter_size_t p_max_bad_steps;
801  o2scl::cli::parameter_bool p_aff_inv;
802  o2scl::cli::parameter_double p_max_time;
803  o2scl::cli::parameter_size_t p_max_iters;
804  o2scl::cli::parameter_int p_max_chain_size;
805  o2scl::cli::parameter_int p_file_update_iters;
806  o2scl::cli::parameter_bool p_output_meas;
808  o2scl::cli::parameter_int p_verbose;
809  //@}
810 
811  /** \brief Initial write to HDF5 file
812  */
813  virtual void first_update(o2scl_hdf::hdf_file &hf) {
814  parent_t::first_update(hf);
815  hf.sets_vec("cl_args",this->cl_args);
816  return;
817  }
818 
819  public:
820 
821  /// Main wrapper for parsing command-line arguments
822  virtual void run(int argc, char *argv[]) {
823 
824  // ---------------------------------------
825  // Process command-line arguments and run
826 
827  setup_cli();
828 
829 #ifdef O2SCL_MPI
830  // Get MPI rank, etc.
831  MPI_Comm_rank(MPI_COMM_WORLD,&this->mpi_rank);
832  MPI_Comm_size(MPI_COMM_WORLD,&this->mpi_nprocs);
833 #endif
834 
835  // Process arguments
836  for(int i=0;i<argc;i++) {
837  this->cl_args.push_back(argv[i]);
838  }
839 
840  this->cl.prompt="mcmc> ";
841  this->cl.run_auto(argc,argv);
842 
843  if (this->file_opened) {
844  //Close main output file
845  this->scr_out.close();
846  }
847 
848  return;
849  }
850 
851  /// \name Customization functions
852  //@{
853  /** \brief Set up the 'cli' object
854 
855  This function just adds the four commands and the 'set' parameters
856  */
857  virtual void setup_cli() {
858 
859  // ---------------------------------------
860  // Set commands/options
861 
862  static const size_t nopt=1;
863  o2scl::comm_option_s options[nopt]={
864  {'i',"initial-point","Set the starting point in the parameter space",
865  1,-1,"<mode> [...]",
866  ((std::string)"Mode can be one of 'best', 'last', 'N', or ")+
867  "'values'. If mode is 'best', then it uses the point with the "+
868  "largest weight and the second argument specifies the file. If "+
869  "mode is 'last' then it uses the last point and the second "+
870  "argument specifies the file. If mode is 'N' then it uses the Nth "+
871  "point, the second argument specifies the value of N and the third "+
872  "argument specifies the file. If mode is 'values', then the "+
873  "remaining arguments specify all the parameter values. On the "+
874  "command-line, enclose negative values in quotes and parentheses, "+
875  "i.e. \"(-1.00)\" to ensure they do not get confused with other "+
878  o2scl::cli::comm_option_both}
879  /*
880  {'s',"hastings","Specify distribution for M-H step",
881  1,1,"<filename>",
882  ((string)"Desc. ")+"Desc2.",
883  new comm_option_mfptr<mcmc_mpi>(this,&mcmc_mpi::hastings),
884  cli::comm_option_both}
885  */
886  };
887  this->cl.set_comm_option_vec(nopt,options);
888 
889  p_file_update_iters.i=&this->file_update_iters;
890  p_file_update_iters.help=((std::string)"Number of MCMC successes ")+
891  "between file upates (default 40, minimum value 1).";
892  this->cl.par_list.insert(std::make_pair("file_update_iters",
893  &p_file_update_iters));
894 
895  p_max_chain_size.i=&this->max_chain_size;
896  p_max_chain_size.help=((std::string)"Maximum Markov chain size ")+
897  "(default 10000).";
898  this->cl.par_list.insert(std::make_pair("max_chain_size",
899  &p_max_chain_size));
900 
901  p_max_time.d=&this->max_time;
902  p_max_time.help=((std::string)"Maximum run time in seconds ")+
903  "(default 86400 sec or 1 day).";
904  this->cl.par_list.insert(std::make_pair("max_time",&p_max_time));
905 
906  p_max_iters.s=&this->max_iters;
907  p_max_iters.help=((std::string)"If non-zero, limit the number of ")+
908  "iterations to be less than the specified number (default zero).";
909  this->cl.par_list.insert(std::make_pair("max_iters",&p_max_iters));
910 
911  p_prefix.str=&this->prefix;
912  p_prefix.help="Output file prefix (default 'mcmc\').";
913  this->cl.par_list.insert(std::make_pair("prefix",&p_prefix));
914 
915  p_output_meas.b=&this->output_meas;
916  p_output_meas.help=((std::string)"If true, output next point ")+
917  "to the '_scr' file before calling TOV solver (default true).";
918  this->cl.par_list.insert(std::make_pair("output_meas",&p_output_meas));
919 
920  p_step_fac.d=&this->step_fac;
921  p_step_fac.help=((std::string)"MCMC step factor. The step size for ")+
922  "each variable is taken as the difference between the high and low "+
923  "limits divided by this factor (default 10.0). This factor can "+
924  "be increased if the acceptance rate is too small, but care must "+
925  "be taken, e.g. if the conditional probability is multimodal. If "+
926  "this step size is smaller than 1.0, it is reset to 1.0 .";
927  this->cl.par_list.insert(std::make_pair("step_fac",&p_step_fac));
928 
929  p_n_warm_up.s=&this->n_warm_up;
930  p_n_warm_up.help=((std::string)"Minimum number of warm up iterations ")+
931  "(default 0).";
932  this->cl.par_list.insert(std::make_pair("n_warm_up",&p_n_warm_up));
933 
934  p_verbose.i=&this->verbose;
935  p_verbose.help=((std::string)"Verbosity parameter ")+
936  "(default 0).";
937  this->cl.par_list.insert(std::make_pair("verbose",&p_verbose));
938 
939  p_max_bad_steps.s=&this->max_bad_steps;
940  p_max_bad_steps.help=((std::string)"Maximum number of bad steps ")+
941  "(default 1000).";
942  this->cl.par_list.insert(std::make_pair("max_bad_steps",&p_max_bad_steps));
943 
944  p_n_walk.s=&this->n_walk;
945  p_n_walk.help=((std::string)"Number of walkers ")+
946  "(default 1).";
947  this->cl.par_list.insert(std::make_pair("n_walk",&p_n_walk));
948 
949  p_user_seed.i=&this->user_seed;
950  p_user_seed.help=((std::string)"Seed for multiplier for random ")+
951  "number generator. If zero is given (the default), then mcmc() "+
952  "uses time(0) to generate a random seed.";
953  this->cl.par_list.insert(std::make_pair("user_seed",&p_user_seed));
954 
955  p_aff_inv.b=&this->aff_inv;
956  p_aff_inv.help=((std::string)"If true, then use affine-invariant ")+
957  "sampling (default false).";
958  this->cl.par_list.insert(std::make_pair("aff_inv",&p_aff_inv));
959 
960  return;
961  }
962 
963  };
964 
965  // End of bamr namespace
966 }
967 
968 #endif
std::shared_ptr< o2scl::table_units<> > tab
Main data table for Markov chain.
Definition: mcmc.h:825
virtual int mcmc(size_t np, vec_t &init, vec_t &low, vec_t &high, func_t &func, fill_t &fill)
Perform an MCMC simulation.
Definition: mcmc_mpi.h:121
virtual int add_line(const vec_t &pars, double log_weight, size_t walker_ix, bool new_meas, data_t &dat, fill_t &fill)
A measurement function which adds the point to the table.
Definition: mcmc.h:943
An extension to o2scl::cli which uses readline.
Definition: cli_readline.h:45
Error handler to throw C++ exceptions.
Definition: exception.h:258
std::vector< std::string > col_names
Column names.
Definition: mcmc.h:819
double max(std::string scol) const
Return column maximum. Makes no assumptions about ordering .
Definition: table.h:1834
std::ofstream scr_out
The screen output file.
Definition: mcmc_mpi.h:137
double get(std::string scol, size_t row) const
Get value from row row of column named col. .
Definition: table.h:414
int getd_vec(std::string name, std::vector< double > &v)
Get vector dataset and place data in v.
int setd_vec_copy(std::string name, const vec_t &v)
Set vector dataset named name with v.
Definition: hdf_file.h:366
Double parameter for o2scl::cli.
Definition: cli.h:299
virtual void update_files()
Update files with current table.
Definition: mcmc_mpi.h:197
double wvector_mean(size_t n, const vec_t &data, const vec2_t &weights)
Compute the mean of weighted data.
Definition: vec_stats.h:1186
The main O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
std::string prefix
Prefix for output filenames.
Definition: mcmc_mpi.h:158
int set_szt_vec(std::string name, const std::vector< size_t > &v)
Set vector dataset named name with v.
size_t get_nlines() const
Return the number of lines.
Definition: table.h:457
virtual void setup_cli()
Set up the &#39;cli&#39; object.
Definition: mcmc_mpi.h:857
int * i
Parameter.
Definition: cli.h:340
err_hnd_type * err_hnd
The global error handler pointer.
void close()
Close the file.
Integer parameter for o2scl::cli.
Definition: cli.h:356
void set_szt(std::string name, size_t u)
Set an unsigned integer named name to value u.
void setd(std::string name, double d)
Set a double named name to value d.
size_t n_walk
Number of walkers for affine-invariant MC or 1 otherwise (default 1)
Definition: mcmc.h:207
String parameter for o2scl::cli.
Definition: cli.h:253
void set_seed(unsigned long int s)
Set the seed.
Definition: rng_gsl.h:106
virtual int mcmc(size_t nparams, vec_t &init, vec_t &low, vec_t &high, func_t &func, fill_t &fill)
Perform an MCMC simulation.
Definition: mcmc.h:914
virtual int mcmc_init()
User-defined initialization function.
Definition: mcmc_mpi.h:288
A generic MCMC simulation class.
Definition: mcmc.h:103
invalid argument supplied by user
Definition: err_hnd.h:59
o2scl::err_hnd_cpp error_handler
Error handler for each thread.
Definition: mcmc_mpi.h:89
bool first_file_update
If true, then first_update() has been called.
Definition: mcmc_mpi.h:178
virtual void best_point(ubvector &best, double w_best)
Output the best point so far.
Definition: mcmc_mpi.h:460
vec_t high_copy
A copy of the upper limits for HDF5 output.
Definition: mcmc_mpi.h:101
int mpi_rank
The MPI processor rank.
Definition: mcmc_mpi.h:78
void seti(std::string name, int i)
Set an integer named name to value i.
int run_auto(int argc, char *argv[], int debug=0)
Automatically parse arguments to main and call interactive mode if required.
int get_szt(std::string name, size_t &u)
Get an unsigned integer named name.
std::string prompt
The prompt (default "> ")
Definition: cli.h:500
size_t max_iters
Maximum number of iterations (default 0)
Definition: mcmc_mpi.h:146
generic failure
Definition: err_hnd.h:61
ubvector initial_point
If true, output MC accepts and rejects (default true)
Definition: mcmc_mpi.h:164
size_t chain_index
Number of complete Markov chain segments.
Definition: mcmc_mpi.h:193
Configurable command-line interface.
Definition: cli.h:230
int mpi_nprocs
The MPI number of processors.
Definition: mcmc_mpi.h:81
int max_chain_size
Maximum size of Markov chain (default 10000)
Definition: mcmc_mpi.h:189
std::string help
Help description.
Definition: cli.h:242
void sets(std::string name, std::string s)
Set a string named name to value s.
int stoi(std::string s, bool err_on_fail=true)
Convert a string to an integer.
MCMC with MPI and HDF5 table I/O.
Definition: mcmc_mpi.h:66
virtual void first_update(o2scl_hdf::hdf_file &hf)
Initial write to HDF5 file.
Definition: mcmc_mpi.h:478
double wvector_covariance(size_t n, const vec_t &data1, const vec2_t &data2, const vec3_t &weights)
The weighted covariance of two vectors.
Definition: vec_stats.h:1363
size_t * s
Parameter.
Definition: cli.h:363
void open_or_create(std::string fname)
Open a file named fname or create if it doesn&#39;t already exist.
size_t n_warm_up
Number of warm up steps (successful steps not iterations) (default 0)
Definition: mcmc.h:189
A generic MCMC simulation class writing data to a o2scl::table_units object.
Definition: mcmc.h:805
double max_time
Time in seconds (default is 86400 seconds or 1 day)
Definition: mcmc_mpi.h:150
virtual int add_line(const ubvector &pars, double weight, size_t ix, bool new_meas, data_t &dat, fill_t &fill)
Add a measurement to the table.
Definition: mcmc_mpi.h:225
vec_t low_copy
A copy of the lower limits for HDF5 output.
Definition: mcmc_mpi.h:97
Member function pointer for o2scl::cli command function.
Definition: cli.h:100
int set_initial_point(std::vector< std::string > &sv, bool itive_com)
Set the first point.
Definition: mcmc_mpi.h:696
void cholesky_decomp(const size_t M, mat_t &A)
Compute the in-place Cholesky decomposition of a symmetric positive-definite square matrix...
Definition: cholesky_base.h:61
int user_seed
If non-zero, use as the seed for the random number generator (default 0)
Definition: mcmc.h:194
String parameter for o2scl::cli.
Definition: cli.h:276
bool * b
Parameter.
Definition: cli.h:283
size_t max_bad_steps
Maximum number of failed steps when generating initial points with affine-invariant sampling (default...
Definition: mcmc.h:202
std::string * str
Parameter.
Definition: cli.h:260
std::map< std::string, parameter *, std::less< std::string > > par_list
Parameter list.
Definition: cli.h:381
double mpi_start_time
The MPI starting time.
Definition: mcmc_mpi.h:84
bool output_meas
Output each measurement.
Definition: mcmc_mpi.h:154
MCMC class with a command-line interface.
Definition: mcmc_mpi.h:773
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
virtual void add_constant(std::string name, double val)
Add a constant, or if the constant already exists, change its value.
Definition: table.h:2151
int set_comm_option_vec(size_t list_size, vec_t &option_list)
Add a vector containing new commands/options.
Definition: cli.h:530
double * d
Parameter.
Definition: cli.h:310
Integer parameter for o2scl::cli.
Definition: cli.h:333
std::vector< std::string > cl_args
The arguments sent to the command-line.
Definition: mcmc_mpi.h:792
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.
Definition: vector.h:2319
void delete_rows(std::string func)
Delete all rows where func evaluates to a number greater than 0.5 .
Definition: table.h:1118
virtual void set_seed()
Default method for setting the random seed.
Definition: mcmc_mpi.h:107
Data table table class with units.
Definition: table_units.h:37
Command for interactive mode in o2scl::cli.
Definition: cli.h:155
size_t n_params
The number of parameters.
Definition: mcmc_mpi.h:93
int file_update_iters
The number of MCMC successes between file updates (default 40)
Definition: mcmc_mpi.h:185
virtual void run(int argc, char *argv[])
Main wrapper for parsing command-line arguments.
Definition: mcmc_mpi.h:822
void function_column(std::string function, std::string scol)
Make a column from the function specified in function and add it to the table.
Definition: table.h:2445
std::string initial_point_file
The file containing the initial point.
Definition: mcmc_mpi.h:167
Store data in an O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$sc...
Definition: hdf_file.h:96
int sets_vec(std::string name, std::vector< std::string > &s)
Set a vector of strings named name.
o2scl::cli cl
Command-line interface.
Definition: mcmc_mpi.h:787
bool file_opened
If true, scr_out has been opened.
Definition: mcmc_mpi.h:140
int open(std::string fname, bool write_access=false, bool err_on_fail=true)
Open a file named fname.
mcmc_mpi()
Create an MCMC object with model m.
Definition: mcmc_mpi.h:734
std::string itos(int x)
Convert an integer to a string.
double function_to_double(std::string s, bool err_on_fail=true)
Convert a formula to a double.
std::string szttos(size_t x)
Convert a size_t to a string.
void hdf_input(hdf_file &hf, o2scl::table< vec_t > &t, std::string name)
Input a o2scl::table object from a hdf_file.
Definition: hdf_io.h:59
virtual void first_update(o2scl_hdf::hdf_file &hf)
Initial write to HDF5 file.
Definition: mcmc_mpi.h:813
virtual int mcmc_init()
MCMC initialization function.
Definition: mcmc.h:831

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).