mcmc_para.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2012-2018, 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_para.h
24  \brief File for definition of \ref o2scl::mcmc_para_base,
25  \ref o2scl::mcmc_para_table and \ref o2scl::mcmc_para_cli
26 */
27 #ifndef O2SCL_MCMC_PARA_H
28 #define O2SCL_MCMC_PARA_H
29 
30 #include <iostream>
31 #include <random>
32 
33 #ifdef O2SCL_OPENMP
34 #include <omp.h>
35 #endif
36 #ifdef O2SCL_MPI
37 #include <mpi.h>
38 #endif
39 
40 #include <boost/numeric/ublas/vector.hpp>
41 
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/interpm_idw.h>
50 #include <o2scl/vec_stats.h>
51 #include <o2scl/cli.h>
52 
53 namespace o2scl {
54 
57 
58  /** \brief A generic MCMC simulation class
59 
60  This class performs a Markov chain Monte Carlo simulation of a
61  user-specified function using OpenMP and/or MPI. Either the
62  Metropolis-Hastings algorithm with a user-specified proposal
63  distribution or the affine-invariant sampling method of Goodman
64  and Weare can be used.
65 
66  By default, the Metropolis-Hastings algorithm is executed with a
67  simple walk, with steps in each dimension of size \f$
68  (\mathrm{high} - \mathrm{low})/\mathrm{step\_fac} \f$ with the
69  denominator specified in \ref step_fac.
70 
71  The function type is a template type, \c func_t, which should
72  be of the form
73  \code
74  int f(size_t num_of_parameters, const vec_t &parameters,
75  double &log_pdf, data_t &dat)
76  \endcode
77  which computes \c log_pdf, the natural logarithm of the function
78  value, for any point in parameter space (any point between \c
79  low and \c high ).
80 
81  If the function being simulated returns \ref mcmc_skip then the
82  point is automatically rejected. After each acceptance or
83  rejection, a user-specified "measurement" function (of type \c
84  measure_t ) is called, which can be used to store the results.
85  In order to stop the simulation, either this function or the
86  probability distribution being simulated should return the value
87  \ref mcmc_done .
88 
89  A generic proposal distribution can be specified in \ref
90  set_proposal(). To go back to the default random walk method,
91  one can call the function \ref unset_proposal().
92 
93  If \ref aff_inv is set to true and the number of walkers, \ref
94  n_walk is set to a number larger than 1, then affine-invariant
95  sampling is used. For affine-invariant sampling, the variable
96  \ref step_fac represents the value of \f$ a \f$, the
97  limits of the distribution for \f$ z \f$.
98 
99  In order to store data at each point, the user can store this
100  data in any object of type \c data_t . If affine-invariant
101  sampling is used, then each chain has it's own data object. The
102  class keeps twice as many copies of these data object as would
103  otherwise be required, in order to avoid copying of data objects
104  in the case that the steps are accepted or rejected.
105 
106  <b>Verbose output:</b> If verbose is 0, no output is generated
107  (the default). If verbose is 1, then output to <tt>cout</tt>
108  occurs only if the settings are somehow misconfigured and the
109  class attempts to recover from them, for example if not enough
110  functions are specified for the requested number of OpenMP
111  threads, or if more than one thread was requested but
112  O2SCL_OPENMP was not defined, or if a negative value for \ref
113  step_fac was requested. When verbose is 1, a couple things are
114  also output to \ref scr_out, including a summary of the number
115  of walkers, chains, and threads at the beginning of the MCMC
116  simulation, a message indicating why the MCMC simulation
117  stopped, a message when the warm up iterations are completed, a
118  message every time files are written to disk, and a message at
119  the end counting the number of acceptances and rejections.
120  If verbose is 2, then the file prefix is output to <tt>cout</tt>
121  during initialization.
122 
123  \note This class is experimental.
124  */
125  template<class func_t, class measure_t,
126  class data_t, class vec_t=ubvector> class mcmc_para_base {
127 
128  protected:
129 
130  /// \name MPI properties
131  //@{
132  /// The MPI processor rank
133  int mpi_rank;
134 
135  /// The MPI number of processors
136  int mpi_size;
137  //@}
138 
139  /// The screen output file
140  std::ofstream scr_out;
141 
142  /// Random number generators
143  std::vector<rng_gsl> rg;
144 
145  /// Pointer to proposal distribution for each thread
146  std::vector<o2scl::prob_cond_mdim<vec_t> *> prop_dist;
147 
148  /// If true, then use the user-specified proposal distribution
149  bool pd_mode;
150 
151  /// If true, we are in the warm up phase
152  bool warm_up;
153 
154  /** \brief Current points in parameter space for each walker and
155  each OpenMP thread
156 
157  This is an array of size \ref n_threads times \ref n_walk initial
158  guesses, indexed by <tt>thread_index*n_walk+walker_index</tt> .
159  */
160  std::vector<vec_t> current;
161 
162  /** \brief Data array
163 
164  This is an array of size 2 times \ref n_threads times \ref
165  n_walk . The two copies of data objects are indexed by
166  <tt>i_copy*n_walk*n_threads+thread_index*n_walk+walker_index</tt>
167  */
168  std::vector<data_t> data_arr;
169 
170  /** \brief Data switch array for each walker and each OpenMP thread
171 
172  This is an array of size \ref n_threads times \ref n_walk initial
173  guesses, indexed by <tt>thread_index*n_walk+walker_index</tt> .
174  */
175  std::vector<bool> switch_arr;
176 
177  /** \brief Return value counters, one vector independent chain
178  */
179  std::vector<std::vector<size_t> > ret_value_counts;
180 
181  /// \name Interface customization
182  //@{
183  /** \brief Initializations before the MCMC
184  */
185  virtual int mcmc_init() {
186  if (verbose>1) {
187  std::cout << "Prefix is: " << prefix << std::endl;
188  }
189  if (verbose>0) {
190  // Open main output file for this rank
191  scr_out.open((prefix+"_"+
192  o2scl::itos(mpi_rank)+"_scr").c_str());
193  scr_out.setf(std::ios::scientific);
194  }
195  return 0;
196  }
197 
198  /** \brief Cleanup after the MCMC
199  */
200  virtual void mcmc_cleanup() {
201  if (verbose>0) {
202  for(size_t it=0;it<n_chains_per_rank;it++) {
203  scr_out << "mcmc (" << it << "," << mpi_rank
204  << "): accept=" << n_accept[it]
205  << " reject=" << n_reject[it] << std::endl;
206  }
207  scr_out.close();
208  }
209  return;
210  }
211 
212  /** \brief Function to run for the best point
213  */
214  virtual void best_point(vec_t &best, double w_best, data_t &dat) {
215  return;
216  }
217 
218  /** \brief Function to run after point evaluation and measurement steps
219  */
220  virtual void post_pointmeas() {
221  return;
222  }
223  //@}
224 
225  /** \brief Index of the current walker
226 
227  This quantity has to be a vector because different threads
228  may have different values for the current walker during
229  the initialization phase for the affine sampling algorithm.
230  */
231  std::vector<size_t> curr_walker;
232 
233  /** \brief Number of fully independent chains in each MPI rank
234  */
236 
237  public:
238 
239  /** \brief The MPI starting time (defaults to 0.0)
240 
241  This can be set by the user before mcmc() is called, so
242  that the time required for initializations before
243  the MCMC starts can be counted.
244  */
246 
247  /** \brief If non-zero, the maximum number of MCMC iterations
248  (default 0)
249 
250  If both \ref max_iters and \ref max_time are nonzero, the
251  MCMC will stop when either the number of iterations
252  exceeds \ref max_iters or the time exceeds \ref max_time,
253  whichever occurs first.
254  */
255  size_t max_iters;
256 
257  /** \brief Time in seconds (default is 0)
258 
259  If both \ref max_iters and \ref max_time are nonzero, the
260  MCMC will stop when either the number of iterations
261  exceeds \ref max_iters or the time exceeds \ref max_time,
262  whichever occurs first.
263  */
264  double max_time;
265 
266  /** \brief Prefix for output filenames (default "mcmc")
267  */
268  std::string prefix;
269 
270  /// Integer to indicate completion
271  static const int mcmc_done=-10;
272 
273  /// Integer to indicate rejection
274  static const int mcmc_skip=-20;
275 
276  /// \name Output quantities
277  //@{
278  /** \brief The number of Metropolis steps which were accepted in
279  each independent chain (summed over all walkers)
280 
281  This vector has a size equal to \ref n_chains_per_rank .
282  */
283  std::vector<size_t> n_accept;
284 
285  /** \brief The number of Metropolis steps which were rejected in
286  each independent chain (summed over all walkers)
287 
288  This vector has a size equal to \ref n_chains_per_rank .
289  */
290  std::vector<size_t> n_reject;
291  //@}
292 
293  /// \name Settings
294  //@{
295  /// If true, use affine-invariant Monte Carlo
296  bool aff_inv;
297 
298  /// Stepsize factor (default 10.0)
299  double step_fac;
300 
301  /** \brief Number of warm up steps (successful steps not
302  iterations) (default 0)
303 
304  \note Not to be confused with <tt>warm_up</tt>, which is
305  a protected boolean local variable in some functions which
306  indicates whether we're in warm up mode or not.
307  */
308  size_t n_warm_up;
309 
310  /** \brief If non-zero, use as the seed for the random number
311  generator (default 0)
312 
313  The random number generator is modified so that each thread and
314  each rank has a different set of random numbers.
315  */
317 
318  /// Output control (default 0)
319  int verbose;
320 
321  /** \brief Maximum number of failed steps when generating initial points
322  with affine-invariant sampling (default 1000)
323  */
325 
326  /** \brief Number of walkers for affine-invariant MC or 1
327  otherwise (default 1)
328  */
329  size_t n_walk;
330 
331  /** \brief Number of walkers per thread (default 1)
332  */
334 
335  /** \brief If true, call the error handler if msolve() or
336  msolve_de() does not converge (default true)
337  */
339 
340  /** \brief If true, accept all steps
341  */
343 
344  /** \brief Initial step fraction for affine-invariance sampling walkers
345  (default 0.1)
346  */
348  //@}
349 
350  mcmc_para_base() {
351  user_seed=0;
352  n_warm_up=0;
353 
354  // MC step parameters
355  aff_inv=false;
356  pd_mode=false;
357  step_fac=10.0;
358  n_walk=1;
359  err_nonconv=true;
360  verbose=1;
361  warm_up=false;
362  max_bad_steps=1000;
363 
364  always_accept=false;
365  ai_initial_step=0.1;
366 
367  n_threads=1;
368  n_walk=1;
369 
370  // Initial values for MPI paramers
371  mpi_size=1;
372  mpi_rank=0;
373  mpi_start_time=0.0;
374 
375 #ifdef O2SCL_MPI
376  // Get MPI rank, etc.
377  MPI_Comm_rank(MPI_COMM_WORLD,&this->mpi_rank);
378  MPI_Comm_size(MPI_COMM_WORLD,&this->mpi_size);
379 #endif
380 
381  prefix="mcmc";
382  max_time=0.0;
383  max_iters=0;
384  }
385 
386  /// Number of OpenMP threads
387  size_t n_threads;
388 
389  /** \brief Initial points in parameter space
390 
391  To fully specify all of the initial points, this should be
392  a vector of size \ref n_walk times \ref n_threads .
393  */
394  std::vector<ubvector> initial_points;
395 
396  /// \name Basic usage
397  //@{
398  /** \brief Perform a MCMC simulation
399 
400  Perform a MCMC simulation over \c n_params parameters starting
401  at initial point \c init, limiting the parameters to be between
402  \c low and \c high, using \c func as the objective function and
403  calling the measurement function \c meas at each MC point.
404  */
405  virtual int mcmc(size_t n_params, vec_t &low, vec_t &high,
406  std::vector<func_t> &func, std::vector<measure_t> &meas) {
407 
408  // Doxygen seems to have trouble reading the code, so we
409  // ensure it doesn't see it.
410 #ifndef DOXYGEN
411 
412  if (func.size()<n_threads) {
413  if (verbose>0) {
414  std::cout << "mcmc_para::mcmc(): Not enough functions for "
415  << n_threads << " threads. Setting n_threads to "
416  << func.size() << "." << std::endl;
417  }
418  n_threads=func.size();
419  }
420  if (meas.size()<n_threads) {
421  if (verbose>0) {
422  std::cout << "mcmc_para::mcmc(): Not enough measurement objects for "
423  << n_threads << " threads. Setting n_threads to "
424  << meas.size() << "." << std::endl;
425  }
426  n_threads=meas.size();
427  }
428 
429  // Set start time if necessary
430  if (mpi_start_time==0.0) {
431 #ifdef O2SCL_MPI
432  mpi_start_time=MPI_Wtime();
433 #else
434  mpi_start_time=time(0);
435 #endif
436  std::cout << "Set start time to: " << mpi_start_time << std::endl;
437  }
438 
439  if (initial_points.size()==0) {
440  // Setup initial guess if not specified
441  initial_points.resize(1);
442  initial_points[0].resize(n_params);
443  for(size_t k=0;k<n_params;k++) {
444  initial_points[0][k]=(low[k]+high[k])/2.0;
445  }
446  } else {
447  // If initial points are specified, make sure they're within
448  // the user-specified limits
449  for(size_t iip=0;iip<initial_points.size();iip++) {
450  for(size_t ipar=0;ipar<n_params;ipar++) {
451  if (initial_points[iip][ipar]<low[ipar] ||
452  initial_points[iip][ipar]>high[ipar]) {
453  O2SCL_ERR((((std::string)"Parameter ")+o2scl::szttos(ipar)+
454  " of "+o2scl::szttos(n_params)+" out of range (value="+
455  o2scl::dtos(initial_points[iip][ipar])+
456  ") in mcmc_base::mcmc().").c_str(),
458  }
459  }
460  }
461  }
462 
463  // Set number of threads
464 #ifdef O2SCL_OPENMP
465  omp_set_num_threads(n_threads);
466 #else
467  if (n_threads>1) {
468  std::cout << "mcmc_para::mcmc(): "
469  << n_threads << " threads were requested but the "
470  << "-DO2SCL_OPENMP flag was not used during "
471  << "compilation. Setting n_threads to 1."
472  << std::endl;
473  n_threads=1;
474  }
475 #endif
476 
477  // Storage for return values from each thread
478  std::vector<int> func_ret(n_threads), meas_ret(n_threads);
479 
480  // Fix the number of walkers if it is too small
481  if (aff_inv) {
482  if (n_walk<=1) {
483  std::cout << "mcmc_para::mcmc(): Requested only 1 walker, "
484  << "forcing 2 walkers." << std::endl;
485  n_walk=2;
486  }
487 #ifdef O2SCL_NEVER_DEFINED
488  /*
489  n_walk is always greater than or equal to n_walk_per_thread,
490  and n_chains_per_rank * n_walk = n_walk_per_thread * n_threads
491  thus n_threads is always larger than n_chains_per_rank.
492  */
493  if (n_walk_per_thread>n_walk) {
494  O2SCL_ERR2("More walkers per thread than total walkers ",
495  "in mcmc_para::mcmc().",o2scl::exc_einval);
496  }
497  n_chains_per_rank=n_walk_per_thread*n_threads/n_walk;
498  if (n_chains_per_walk*n_walk!=n_walk_per_thread*n_threads) {
499  std::cout << "mcmc_para::mcmc(): Could not evenly "
500  << "organize threads and walkers." << std::endl;
501  std::cout << "n_threads: " << n_threads << std::endl;
502  std::cout << "n_walk: " << n_walk << std::endl;
503  std::cout << "n_walk_per_thread: "
504  << n_walk_per_thread << << std::endl;
505  std::cout << "n_chains_per_rank: "
506  << n_chains_per_rank << << std::endl;
507  }
508 #endif
509  }
510 
511  n_walk_per_thread=n_walk;
512  n_chains_per_rank=n_threads;
513 
514  // Fix 'step_fac' if it's less than or equal to zero
515  if (step_fac<=0.0) {
516  if (aff_inv) {
517  std::cout << "mcmc_para::mcmc(): Requested negative or zero "
518  << "step_fac with aff_inv=true. Setting to 2.0."
519  << std::endl;
520  step_fac=2.0;
521  } else {
522  std::cout << "mcmc_para::mcmc(): Requested negative or zero "
523  << "step_fac. Setting to 10.0." << std::endl;
524  step_fac=10.0;
525  }
526  }
527 
528  // Set RNGs with a different seed for each thread and rank
529  rg.resize(n_threads);
530  unsigned long int seed=time(0);
531  if (this->user_seed!=0) {
532  seed=this->user_seed;
533  }
534  for(size_t it=0;it<n_threads;it++) {
535  seed*=(mpi_rank*n_threads+it+1);
536  rg[it].set_seed(seed);
537  }
538 
539  // Keep track of successful and failed MH moves in each
540  // independent chain
541  n_accept.resize(n_chains_per_rank);
542  n_reject.resize(n_chains_per_rank);
543  for(size_t it=0;it<n_chains_per_rank;it++) {
544  n_accept[it]=0;
545  n_reject[it]=0;
546  }
547 
548  // Warm-up flag, not to be confused with 'n_warm_up', the
549  // number of warm_up iterations.
550  warm_up=true;
551  if (n_warm_up==0) warm_up=false;
552 
553  // Storage size required
554  size_t ssize=n_walk*n_chains_per_rank;
555 
556  // Allocate current point and current weight
557  current.resize(ssize);
558  std::vector<double> w_current(ssize);
559  for(size_t i=0;i<ssize;i++) {
560  current[i].resize(n_params);
561  w_current[i]=0.0;
562  }
563 
564  // Allocate ret_value_counts and curr_walker
565  ret_value_counts.resize(n_chains_per_rank);
566  curr_walker.resize(n_threads);
567 
568  // Initialize data and switch arrays
569  data_arr.resize(2*ssize);
570  switch_arr.resize(ssize);
571  for(size_t i=0;i<switch_arr.size();i++) switch_arr[i]=false;
572 
573  // Next point and next weight for each thread
574  std::vector<vec_t> next(n_threads);
575  for(size_t it=0;it<n_threads;it++) {
576  next[it].resize(n_params);
577  }
578  std::vector<double> w_next(n_threads);
579 
580  // Best point over all threads
581  vec_t best(n_params);
582  double w_best;
583 
584  // Generally, these flags are are true for any thread if func_ret
585  // or meas_ret is equal to mcmc_done.
586  std::vector<bool> mcmc_done_flag(n_threads);
587  for(size_t it=0;it<n_threads;it++) {
588  mcmc_done_flag[it]=false;
589  }
590 
591  // Proposal weight
592  std::vector<double> q_prop(n_threads);
593 
594  // Run mcmc_init() function. The initial point, stored in
595  // current[0] can be modified by this function and the local
596  // variable 'init' is not accessible to the mcmc_init() function.
597  int init_ret=mcmc_init();
598  if (init_ret!=0) {
599  O2SCL_ERR("Function mcmc_init() failed in mcmc_base::mcmc().",
601  return init_ret;
602  }
603 
604  // ---------------------------------------------------
605  // Initial verbose output
606 
607  if (verbose>=1) {
608  if (aff_inv) {
609  scr_out << "mcmc: Affine-invariant step, n_params="
610  << n_params << ", n_walk=" << n_walk
611  << ", n_chains_per_rank=" << n_chains_per_rank
612  << ",\n\tn_walk_per_thread=" << n_walk_per_thread
613  << ", n_threads=" << n_threads << ", rank="
614  << mpi_rank << ", n_ranks="
615  << mpi_size << std::endl;
616  } else if (pd_mode==true) {
617  scr_out << "mcmc: With proposal distribution, n_params="
618  << n_params << ", n_threads=" << n_threads << ", rank="
619  << mpi_rank << ", n_ranks="
620  << mpi_size << std::endl;
621  } else {
622  scr_out << "mcmc: Random-walk w/uniform dist., n_params="
623  << n_params << ", n_threads=" << n_threads << ", rank="
624  << mpi_rank << ", n_ranks="
625  << mpi_size << std::endl;
626  }
627  }
628 
629  // --------------------------------------------------------
630  // Compute initial point and initial weights
631 
632  // --------------------------------------------------------
633  // Initial point and weights for affine-invariant sampling
634 
635  if (aff_inv) {
636 
637 #ifdef O2SCL_OPENMP
638 #pragma omp parallel default(shared)
639 #endif
640  {
641 #ifdef O2SCL_OPENMP
642 #pragma omp for
643 #endif
644  for(size_t it=0;it<n_threads;it++) {
645 
646  // Initialize each walker in turn
647  for(curr_walker[it]=0;curr_walker[it]<n_walk_per_thread &&
648  mcmc_done_flag[it]==false;curr_walker[it]++) {
649 
650  // Index in storage
651  size_t sindex=n_walk*it+curr_walker[it];
652 
653  // Index in initial_point specification
654  size_t ip_index=sindex % initial_points.size();
655 
656  size_t init_iters=0;
657  bool done=false;
658 
659  // If we already have a guess, try to use that
660  if (initial_points.size()>0) {
661 
662  // Copy from the initial points array
663  for(size_t ipar=0;ipar<n_params;ipar++) {
664  current[sindex][ipar]=initial_points[ip_index][ipar];
665  }
666 
667  // Compute the weight
668  func_ret[it]=func[it](n_params,current[sindex],
669  w_current[sindex],data_arr[sindex]);
670 
671  if (func_ret[it]==mcmc_done) {
672  mcmc_done_flag[it]=true;
673  } else if (func_ret[it]==o2scl::success) {
674 
675  // If we have a good point, update ret_value_counts
676  // and call the measurement function
677  if (func_ret[it]>=0 &&
678  func_ret[it]<((int)ret_value_counts[it].size())) {
679  ret_value_counts[it][func_ret[it]]++;
680  }
681  meas_ret[it]=meas[it](current[sindex],w_current[sindex],
682  curr_walker[it],func_ret[it],
683  true,data_arr[sindex]);
684  if (meas_ret[it]==mcmc_done) {
685  mcmc_done_flag[it]=true;
686  }
687  done=true;
688  }
689  }
690 
691  while (!done && !mcmc_done_flag[it]) {
692 
693  // Make a perturbation from the initial point
694  for(size_t ipar=0;ipar<n_params;ipar++) {
695  do {
696  current[sindex][ipar]=
697  initial_points[ip_index][ipar]+
698  (rg[it].random()*2.0-1.0)*
699  (high[ipar]-low[ipar])*ai_initial_step;
700  } while (current[sindex][ipar]>high[ipar] ||
701  current[sindex][ipar]<low[ipar]);
702  }
703 
704  // Compute the weight
705  func_ret[it]=func[it](n_params,current[sindex],
706  w_current[sindex],data_arr[sindex]);
707 
708  // ------------------------------------------------
709 
710  // Increment iteration count
711  init_iters++;
712 
713  if (func_ret[it]==mcmc_done) {
714  mcmc_done_flag[it]=true;
715  } else {
716  // If we have a good point, update ret_value_counts,
717  // call the measurement function and stop the loop
718  if (func_ret[it]==o2scl::success) {
719  if (func_ret[it]>=0 &&
720  func_ret[it]<((int)ret_value_counts[it].size())) {
721  ret_value_counts[it][func_ret[it]]++;
722  }
723  if (meas_ret[it]!=mcmc_done) {
724  meas_ret[it]=meas[it](current[sindex],w_current[sindex],
725  curr_walker[it],func_ret[it],true,
726  data_arr[sindex]);
727  } else {
728  mcmc_done_flag[it]=true;
729  }
730  done=true;
731  } else if (init_iters>max_bad_steps) {
732  O2SCL_ERR2("Initial walkers failed in ",
733  "mcmc_para_base::mcmc().",o2scl::exc_einval);
734  }
735  }
736  }
737  }
738  }
739  }
740  // End of parallel region
741 
742  // Stop early if mcmc_done was returned
743  bool stop_early=false;
744  for(size_t it=0;it<n_threads;it++) {
745  if (mcmc_done_flag[it]==true) {
746  if (verbose>=1) {
747  scr_out << "mcmc (" << it << "," << mpi_rank
748  << "): Returned mcmc_done "
749  << "(initial; ai)." << std::endl;
750  }
751  stop_early=true;
752  }
753  }
754  if (stop_early) {
755  mcmc_cleanup();
756  return 0;
757  }
758 
759  // Set initial values for best point
760  w_best=w_current[0];
761  size_t best_index=0;
762  for(size_t it=0;it<n_threads;it++) {
763  for(curr_walker[it]=0;curr_walker[it]<n_walk;
764  curr_walker[it]++) {
765  size_t sindex=n_walk*it+curr_walker[it];
766  if (w_current[sindex]>w_current[0]) {
767  best_index=sindex;
768  w_best=w_current[sindex];
769  }
770  }
771  }
772  best=current[best_index];
773  best_point(best,w_best,data_arr[best_index]);
774 
775  // Verbose output
776  if (verbose>=2) {
777  for(size_t it=0;it<n_threads;it++) {
778  for(curr_walker[it]=0;curr_walker[it]<n_walk;curr_walker[it]++) {
779  size_t sindex=n_walk*it+curr_walker[it];
780  scr_out.precision(4);
781  scr_out << "mcmc (" << it << "," << mpi_rank << "): i_walk: ";
782  scr_out.width((int)(1.0+log10((double)(n_walk-1))));
783  scr_out << curr_walker[it] << " log_wgt: "
784  << w_current[sindex] << " (initial; ai)" << std::endl;
785  scr_out.precision(6);
786  }
787  }
788  }
789 
790  // End of 'if (aff_inv)'
791 
792  } else {
793 
794  // --------------------------------------------------------
795  // Initial point and weights when aff_inv is false.
796 
797 
798 #ifdef O2SCL_OPENMP
799 #pragma omp parallel default(shared)
800 #endif
801  {
802 #ifdef O2SCL_OPENMP
803 #pragma omp for
804 #endif
805  for(size_t it=0;it<n_threads;it++) {
806 
807  // Note that this value is used (e.g. in
808  // mcmc_para_table::add_line() ) even if aff_inv is false, so we
809  // set it to zero here.
810  curr_walker[it]=0;
811 
812  // Copy from the initial points array into current point
813  size_t ip_size=initial_points.size();
814  for(size_t ipar=0;ipar<n_params;ipar++) {
815  current[it][ipar]=initial_points[it % ip_size][ipar];
816  }
817 
818  if (it<ip_size) {
819  // If we have a new unique initial point, then
820  // perform a function evaluation
821  func_ret[it]=func[it](n_params,current[it],w_current[it],
822  data_arr[it]);
823  } else {
824  func_ret[it]=0;
825  }
826 
827  }
828 
829  }
830  // End of parallel region
831 
832  // Check return values from initial point function evaluations
833  for(size_t it=0;it<n_threads;it++) {
834  if (func_ret[it]==mcmc_done) {
835  if (verbose>=1) {
836  scr_out << "mcmc (" << it
837  << "): Initial point returned mcmc_done."
838  << std::endl;
839  }
840  mcmc_cleanup();
841  return 0;
842  }
843  if (func_ret[it]!=o2scl::success) {
844  if (err_nonconv) {
845  O2SCL_ERR((((std::string)"Initial weight from thread ")+
846  o2scl::szttos(it)+
847  " vanished in mcmc_para_base::mcmc().").c_str(),
849  }
850  return 2;
851  }
852  }
853 
854 #ifdef O2SCL_OPENMP
855 #pragma omp parallel default(shared)
856 #endif
857  {
858 #ifdef O2SCL_OPENMP
859 #pragma omp for
860 #endif
861  for(size_t it=0;it<n_threads;it++) {
862 
863  size_t ip_size=initial_points.size();
864  if (it>=ip_size) {
865  // If no initial point was specified, copy one of
866  // the other initial points
867  func_ret[it]=func_ret[it % ip_size];
868  current[it]=current[it % ip_size];
869  w_current[it]=w_current[it % ip_size];
870  data_arr[it]=data_arr[it % ip_size];
871  }
872 
873  // Update the return value count
874  if (func_ret[it]>=0 &&
875  func_ret[it]<((int)ret_value_counts[it].size())) {
876  ret_value_counts[it][func_ret[it]]++;
877  }
878  // Call the measurement function
879  meas_ret[it]=meas[it](current[it],w_current[it],0,
880  func_ret[it],true,data_arr[it]);
881  if (meas_ret[it]==mcmc_done) {
882  mcmc_done_flag[it]=true;
883  }
884 
885  // End of loop over threads
886  }
887 
888 
889  }
890  // End of parallel region
891 
892  // Stop early if mcmc_done was returned from one of the
893  // measurement function calls
894  bool stop_early=false;
895  for(size_t it=0;it<n_threads;it++) {
896  if (mcmc_done_flag[it]==true) {
897  if (verbose>=1) {
898  scr_out << "mcmc (" << it << "," << mpi_rank
899  << "): Returned mcmc_done "
900  << "(initial)." << std::endl;
901  }
902  stop_early=true;
903  }
904  }
905  if (stop_early) {
906  mcmc_cleanup();
907  return 0;
908  }
909 
910  // Set initial values for best point
911  best=current[0];
912  w_best=w_current[0];
913  best_point(best,w_best,data_arr[0]);
914  for(size_t it=1;it<n_threads;it++) {
915  if (w_current[it]<w_best) {
916  best=current[it];
917  w_best=w_current[it];
918  best_point(best,w_best,data_arr[it]);
919  }
920  }
921 
922  if (verbose>=2) {
923  scr_out.precision(4);
924  scr_out << "mcmc: ";
925  for(size_t it=0;it<n_threads;it++) {
926  scr_out << w_current[it] << " ";
927  }
928  scr_out << " (initial)" << std::endl;
929  scr_out.precision(6);
930  }
931 
932  // End of initial point region for 'aff_inv=false'
933  }
934 
935  // --------------------------------------------------------
936  // Require keypress after initial point if verbose is
937  // sufficiently large.
938 
939  if (verbose>=3) {
940  std::cout << "Press a key and type enter to continue. ";
941  char ch;
942  std::cin >> ch;
943  }
944 
945  // End of initial point and weight section
946  // ---------------------------------------------------
947  // Start of main loop
948 
949  bool main_done=false;
950  size_t mcmc_iters=0;
951 
952  while (!main_done) {
953 
954  // Choose walker to move (or zero when aff_inv is false)
955  std::vector<double> smove_z(n_threads);
956  for(size_t it=0;it<n_threads;it++) {
957  curr_walker[it]=0;
958  smove_z[it]=0.0;
959  // Choose walker to move (same for all threads)
960  if (aff_inv) {
961  curr_walker[it]=mcmc_iters % n_walk;
962  }
963  }
964 
965 #ifdef O2SCL_OPENMP
966 #pragma omp parallel default(shared)
967 #endif
968  {
969 #ifdef O2SCL_OPENMP
970 #pragma omp for
971 #endif
972  for(size_t it=0;it<n_threads;it++) {
973 
974  // ---------------------------------------------------
975  // Select next point
976 
977  if (aff_inv) {
978  // Choose jth walker
979  size_t ij;
980  do {
981  ij=((size_t)(rg[it].random()*((double)n_walk)));
982  } while (ij==curr_walker[it] || ij>=n_walk);
983 
984  // Select z
985  double p=rg[it].random();
986  double a=step_fac;
987  smove_z[it]=(1.0-2.0*p+2.0*a*p+p*p-2.0*a*p*p+a*a*p*p)/a;
988 
989  // Create new trial point
990  for(size_t i=0;i<n_params;i++) {
991  next[it][i]=current[n_walk*it+ij][i]+
992  smove_z[it]*(current[n_walk*it+curr_walker[it]][i]-
993  current[n_walk*it+ij][i]);
994  }
995 
996  } else if (pd_mode) {
997 
998  // Use proposal distribution and compute associated weight
999  q_prop[it]=prop_dist[it]->log_metrop_hast(current[it],next[it]);
1000 
1001  if (!std::isfinite(q_prop[it])) {
1002  O2SCL_ERR2("Proposal distribution not finite in ",
1003  "mcmc_para_base::mcmc().",o2scl::exc_efailed);
1004  }
1005 
1006  } else {
1007 
1008  // Uniform random-walk step
1009  for(size_t k=0;k<n_params;k++) {
1010  next[it][k]=current[it][k]+(rg[it].random()*2.0-1.0)*
1011  (high[k]-low[k])/step_fac;
1012  }
1013 
1014  }
1015 
1016  // ---------------------------------------------------
1017  // Compute next weight
1018 
1019  func_ret[it]=o2scl::success;
1020  // If the next point out of bounds, ensure that the point is
1021  // rejected without attempting to evaluate the function
1022  for(size_t k=0;k<n_params;k++) {
1023  if (next[it][k]<low[k] || next[it][k]>high[k]) {
1024  func_ret[it]=mcmc_skip;
1025  if (verbose>=3) {
1026  if (next[it][k]<low[k]) {
1027  std::cout << "mcmc (" << it << ","
1028  << mpi_rank << "): Parameter with index " << k
1029  << " and value " << next[it][k]
1030  << " smaller than limit " << low[k] << std::endl;
1031  scr_out << "mcmc (" << it << ","
1032  << mpi_rank << "): Parameter with index " << k
1033  << " and value " << next[it][k]
1034  << " smaller than limit " << low[k] << std::endl;
1035  } else {
1036  std::cout << "mcmc (" << it << "," << mpi_rank
1037  << "): Parameter with index " << k
1038  << " and value " << next[it][k]
1039  << " larger than limit " << high[k] << std::endl;
1040  scr_out << "mcmc (" << it << "," << mpi_rank
1041  << "): Parameter with index " << k
1042  << " and value " << next[it][k]
1043  << " larger than limit " << high[k] << std::endl;
1044  }
1045  }
1046  }
1047  }
1048 
1049  // Evaluate the function, set the 'done' flag if
1050  // necessary, and update the return value array
1051  if (func_ret[it]!=mcmc_skip) {
1052  if (switch_arr[n_walk*it+curr_walker[it]]==false) {
1053  func_ret[it]=func[it](n_params,next[it],w_next[it],
1054  data_arr[it*n_walk+curr_walker[it]+
1055  n_walk*n_threads]);
1056  } else {
1057  func_ret[it]=func[it](n_params,next[it],w_next[it],
1058  data_arr[it*n_walk+curr_walker[it]]);
1059  }
1060  if (func_ret[it]==mcmc_done) {
1061  mcmc_done_flag[it]=true;
1062  } else {
1063  if (func_ret[it]>=0 &&
1064  func_ret[it]<((int)ret_value_counts[it].size())) {
1065  ret_value_counts[it][func_ret[it]]++;
1066  }
1067  }
1068 
1069  }
1070  }
1071  }
1072  // End of parallel region
1073 
1074  // ---------------------------------------------------------
1075  // Post-function verbose output in case parameter was out of
1076  // range, function returned "done" or a failure. More
1077  // verbose output is performed below after the possible call
1078  // to the measurement function.
1079 
1080  if (verbose>=1) {
1081  for(size_t it=0;it<n_threads;it++) {
1082  if (func_ret[it]==mcmc_done) {
1083  scr_out << "mcmc (" << it << "," << mpi_rank
1084  << "): Returned mcmc_done."
1085  << std::endl;
1086  } else if (func_ret[it]==mcmc_skip && verbose>=3) {
1087  scr_out << "mcmc (" << it
1088  << "): Parameter(s) out of range: " << std::endl;
1089  scr_out.setf(std::ios::showpos);
1090  for(size_t k=0;k<n_params;k++) {
1091  scr_out << k << " " << low[k] << " "
1092  << next[it][k] << " " << high[k];
1093  if (next[it][k]<low[k] || next[it][k]>high[k]) {
1094  scr_out << " <-";
1095  }
1096  scr_out << std::endl;
1097  }
1098  scr_out.unsetf(std::ios::showpos);
1099  } else if (func_ret[it]!=o2scl::success &&
1100  func_ret[it]!=mcmc_skip) {
1101  if (verbose>=2) {
1102  scr_out << "mcmc (" << it << "," << mpi_rank
1103  << "): Function returned failure "
1104  << func_ret[it] << " at point ";
1105  for(size_t k=0;k<n_params;k++) {
1106  scr_out << next[it][k] << " ";
1107  }
1108  scr_out << std::endl;
1109  }
1110  }
1111  }
1112  }
1113 
1114  // ----------------------------------------------------------
1115  // Parallel region to accept or reject, and call measurement
1116  // function
1117 
1118 #ifdef O2SCL_OPENMP
1119 #pragma omp parallel default(shared)
1120 #endif
1121  {
1122 #ifdef O2SCL_OPENMP
1123 #pragma omp for
1124 #endif
1125  for(size_t it=0;it<n_threads;it++) {
1126 
1127  // Index in storage
1128  size_t sindex=n_walk*it+curr_walker[it];
1129 
1130  // ---------------------------------------------------
1131  // Accept or reject
1132 
1133  bool accept=false;
1134  if (always_accept && func_ret[it]==success) accept=true;
1135 
1136  if (func_ret[it]==o2scl::success) {
1137  double r=rg[it].random();
1138 
1139  if (aff_inv) {
1140  double ai_ratio=pow(smove_z[it],((double)n_params)-1.0)*
1141  exp(w_next[it]-w_current[sindex]);
1142  if (r<ai_ratio) {
1143  accept=true;
1144  }
1145  } else if (pd_mode) {
1146  if (r<exp(w_next[it]-w_current[sindex]+q_prop[it])) {
1147  accept=true;
1148  }
1149 
1150  } else {
1151  // Metropolis algorithm
1152  if (r<exp(w_next[it]-w_current[sindex])) {
1153  accept=true;
1154  }
1155  }
1156 
1157  // End of 'if (func_ret[it]==o2scl::success)'
1158  }
1159 
1160  if (accept) {
1161 
1162  n_accept[it]++;
1163 
1164  // Store results from new point
1165  if (!warm_up) {
1166  if (switch_arr[sindex]==false) {
1167  meas_ret[it]=meas[it](next[it],w_next[it],
1168  curr_walker[it],func_ret[it],true,
1169  data_arr[sindex+n_threads*n_walk]);
1170  } else {
1171  meas_ret[it]=meas[it](next[it],w_next[it],
1172  curr_walker[it],func_ret[it],true,
1173  data_arr[sindex]);
1174  }
1175  }
1176 
1177  // Prepare for next point
1178  current[sindex]=next[it];
1179  w_current[sindex]=w_next[it];
1180  switch_arr[sindex]=!(switch_arr[sindex]);
1181 
1182  } else {
1183 
1184  // Point was rejected
1185  n_reject[it]++;
1186 
1187  // Repeat measurement of old point
1188  if (!warm_up) {
1189  if (switch_arr[sindex]==false) {
1190  meas_ret[it]=meas[it](next[it],w_next[it],
1191  curr_walker[it],func_ret[it],false,
1192  data_arr[sindex+n_threads*n_walk]);
1193  } else {
1194  meas_ret[it]=meas[it](next[it],w_next[it],
1195  curr_walker[it],func_ret[it],false,
1196  data_arr[sindex]);
1197  }
1198  }
1199 
1200  }
1201 
1202  }
1203  }
1204  // End of parallel region
1205 
1206  post_pointmeas();
1207 
1208  // -----------------------------------------------------------
1209  // Post-measurement verbose output of iteration count, weight,
1210  // and walker index for each thread
1211 
1212  if (verbose>=2) {
1213  for(size_t it=0;it<n_threads;it++) {
1214  size_t sindex=n_walk*it+curr_walker[it];
1215  scr_out.precision(4);
1216  scr_out << "mcmc (" << it << "," << mpi_rank << "): iter: ";
1217  scr_out.width((int)(1.0+log10((double)(n_params-1))));
1218  scr_out << mcmc_iters << " i_walk: "
1219  << curr_walker[it] << " log_wgt: "
1220  << w_current[sindex] << std::endl;
1221  scr_out.precision(6);
1222  }
1223  }
1224 
1225  // Collect best point
1226  for(size_t it=0;it<n_threads;it++) {
1227  if (func_ret[it]==o2scl::success && w_best>w_next[it]) {
1228  best=next[it];
1229  w_best=w_next[it];
1230  if (switch_arr[n_walk*it+curr_walker[it]]==false) {
1231  best_point(best,w_best,data_arr[curr_walker[it]+n_walk*it+
1232  n_threads*n_walk]);
1233  } else {
1234  best_point(best,w_best,data_arr[curr_walker[it]+n_walk*it]);
1235  }
1236  }
1237  }
1238 
1239  // Check to see if mcmc_done was returned or if meas_ret
1240  // returned an error
1241  for(size_t it=0;it<n_threads;it++) {
1242  if (meas_ret[it]==mcmc_done || func_ret[it]==mcmc_done) {
1243  main_done=true;
1244  }
1245  if (meas_ret[it]!=mcmc_done && meas_ret[it]!=o2scl::success) {
1246  if (err_nonconv) {
1247  O2SCL_ERR((((std::string)"Measurement function returned ")+
1248  o2scl::dtos(meas_ret[it])+
1249  " in mcmc_para_base::mcmc().").c_str(),
1251  }
1252  main_done=true;
1253  }
1254  }
1255 
1256  // Update iteration count and reset counters for
1257  // warm up iterations if necessary
1258  if (main_done==false) {
1259 
1260  mcmc_iters++;
1261 
1262  if (warm_up && mcmc_iters==n_warm_up) {
1263  warm_up=false;
1264  mcmc_iters=0;
1265  for(size_t it=0;it<n_threads;it++) {
1266  n_accept[it]=0;
1267  n_reject[it]=0;
1268  }
1269  if (verbose>=1) {
1270  scr_out << "mcmc: Finished warmup." << std::endl;
1271  }
1272 
1273  }
1274  }
1275 
1276  if (verbose>=3) {
1277  std::cout << "Press a key and type enter to continue. ";
1278  char ch;
1279  std::cin >> ch;
1280  }
1281 
1282  // Stop if iterations greater than max
1283  if (main_done==false && warm_up==false && max_iters>0 &&
1284  mcmc_iters==max_iters) {
1285  if (verbose>=1) {
1286  scr_out << "mcmc: Stopping because number of iterations ("
1287  << mcmc_iters << ") equal to max_iters (" << max_iters
1288  << ")." << std::endl;
1289  }
1290  main_done=true;
1291  }
1292 
1293  if (main_done==false) {
1294  // Check to see if we're out of time
1295 #ifdef O2SCL_MPI
1296  double elapsed=MPI_Wtime()-mpi_start_time;
1297 #else
1298  double elapsed=time(0)-mpi_start_time;
1299 #endif
1300  if (max_time>0.0 && elapsed>max_time) {
1301  if (verbose>=1) {
1302  scr_out << "mcmc: Stopping because elapsed (" << elapsed
1303  << ") > max_time (" << max_time << ")."
1304  << std::endl;
1305  }
1306  main_done=true;
1307  }
1308  }
1309 
1310  // --------------------------------------------------------------
1311  // End of main loop
1312  }
1313 
1314  // --------------------------------------------------------------
1315 
1316  mcmc_cleanup();
1317 
1318 #endif
1319 
1320  return 0;
1321  }
1322 
1323  /** \brief Perform a MCMC simulation with a thread-safe function
1324  */
1325  virtual int mcmc(size_t n_params, vec_t &low, vec_t &high,
1326  func_t &func, measure_t &meas) {
1327 
1328 #ifdef O2SCL_OPENMP
1329  omp_set_num_threads(n_threads);
1330 #else
1331  n_threads=1;
1332 #endif
1333  std::vector<func_t> vf(n_threads);
1334  std::vector<measure_t> vm(n_threads);
1335  for(size_t i=0;i<n_threads;i++) {
1336  vf[i]=func;
1337  vm[i]=meas;
1338  }
1339  return mcmc(n_params,low,high,func,meas);
1340  }
1341  //@}
1342 
1343  /// \name Proposal distribution
1344  //@{
1345  /** \brief Set the proposal distribution
1346  */
1347  template<class prob_vec_t> void set_proposal(prob_vec_t &pv) {
1348  prop_dist.resize(pv.size());
1349  for(size_t i=0;i<pv.size();i++) {
1350  prop_dist[i]=&pv[i];
1351  }
1352  pd_mode=true;
1353  aff_inv=false;
1354  n_walk=1;
1355  return;
1356  }
1357 
1358  /** \brief Go back to random-walk Metropolis with a uniform distribution
1359  */
1360  virtual void unset_proposal() {
1361  if (pd_mode) {
1362  prop_dist.clear();
1363  pd_mode=false;
1364  }
1365  aff_inv=false;
1366  n_walk=1;
1367  return;
1368  }
1369  //@}
1370 
1371  };
1372 
1373  /** \brief A generic MCMC simulation class writing data to a
1374  \ref o2scl::table_units object
1375 
1376  This class performs a MCMC simulation and stores the
1377  results in a \ref o2scl::table_units object. The
1378  user must specify the column names and units in
1379  \ref set_names_units() before \ref mcmc() is called.
1380 
1381  The function \ref add_line is the measurement function of type
1382  \c measure_t in the parent. The overloaded function \ref mcmc()
1383  in this class works a bit differently in that it takes a
1384  function object (type \c fill_t) of the form
1385  \code
1386  int fill_func(const vec_t &pars, double log_weight,
1387  std::vector<double> &line, data_t &dat);
1388  \endcode
1389  which should store any auxillary values stored in the data
1390  object to \c line, in order to be added to the table.
1391 
1392  The output table will contain the parameters, the logarithm of
1393  the function (called "log_wgt") and a multiplying factor called
1394  "mult". This "fill" function is called only when a step is
1395  accepted and the multiplier for that row is set to 1. If a
1396  future step is rejected, then the multiplier is increased by
1397  one, rather than adding the same row to the table again.
1398 
1399  There is some output which occurs in addition to the output
1400  from \ref o2scl::mcmc_para_base depending on the value
1401  of \ref o2scl::mcmc_para_base::verbose . If there is
1402  a misalignment between the number of columns in the
1403  table and the number of data points in any line, then
1404  some debugging information is sent to <tt>cout</tt>.
1405  When verbose is 2 or larger,
1406 
1407  \note This class is experimental.
1408 
1409  \future Verbose output may need improvement
1410  \future Use reorder_table() and possibly reblock()
1411  to create a full post-processing function.
1412  */
1413  template<class func_t, class fill_t, class data_t, class vec_t=ubvector>
1414  class mcmc_para_table : public mcmc_para_base<func_t,
1415  std::function<int(const vec_t &,double,size_t,int,bool,data_t &)>,
1416  data_t,vec_t> {
1417 
1418  protected:
1419 
1420  /// Measurement functor type for the parent
1421  typedef std::function<int(const vec_t &,double,size_t,int,bool,data_t &)>
1423 
1424  /// Type of parent class
1426 
1427  /// Column names
1428  std::vector<std::string> col_names;
1429 
1430  /// Column units
1431  std::vector<std::string> col_units;
1432 
1433  /// Number of parameters
1434  size_t n_params;
1435 
1436  /// Main data table for Markov chain
1437  std::shared_ptr<o2scl::table_units<> > table;
1438 
1439  /** \brief If true, the HDF5 I/O initial info has been written
1440  to the file (set by \ref mcmc() )
1441  */
1443 
1444  /** \brief MCMC initialization function
1445 
1446  This function sets the column names and units.
1447  */
1448  virtual int mcmc_init() {
1449 
1450  // -----------------------------------------------------------
1451  // Init tables
1452 
1453  table=std::shared_ptr<o2scl::table_units<> >(new o2scl::table_units<>);
1454  table->new_column("rank");
1455  table->new_column("thread");
1456  table->new_column("walker");
1457  table->new_column("mult");
1458  table->new_column("log_wgt");
1459  for(size_t i=0;i<col_names.size();i++) {
1460  table->new_column(col_names[i]);
1461  if (col_units[i].length()>0) {
1462  table->set_unit(col_names[i],col_units[i]);
1463  }
1464  }
1465 
1466  walker_accept_rows.resize(this->n_walk*this->n_threads);
1467  for(size_t i=0;i<this->n_walk*this->n_threads;i++) {
1468  walker_accept_rows[i]=-1;
1469  }
1470  walker_reject_rows.resize(this->n_walk*this->n_threads);
1471  for(size_t i=0;i<this->n_walk*this->n_threads;i++) {
1472  walker_reject_rows[i]=-1;
1473  }
1474 
1475  if (this->verbose>=2) {
1476  std::cout << "mcmc: Table column names and units: " << std::endl;
1477  for(size_t i=0;i<table->get_ncolumns();i++) {
1478  std::cout << table->get_column_name(i) << " "
1479  << table->get_unit(table->get_column_name(i)) << std::endl;
1480  }
1481  }
1482 
1483  last_write_iters=0;
1484 #ifdef O2SCL_MPI
1485  last_write_time=MPI_Wtime();
1486 #else
1487  last_write_time=time(0);
1488 #endif
1489 
1490  return parent_t::mcmc_init();
1491  }
1492 
1493  /** \brief Fill \c line with data for insertion into the table
1494  */
1495  virtual int fill_line(const vec_t &pars, double log_weight,
1496  std::vector<double> &line, data_t &dat,
1497  size_t i_walker, fill_t &fill) {
1498 
1499 #ifdef O2SCL_OPENMP
1500  size_t i_thread=omp_get_thread_num();
1501 #else
1502  size_t i_thread=0;
1503 #endif
1504 
1505  // Rank
1506  line.push_back(this->mpi_rank);
1507  // Thread
1508  line.push_back(i_thread);
1509  // Walker (set later)
1510  line.push_back(i_walker);
1511  // Initial multiplier
1512  line.push_back(1.0);
1513  line.push_back(log_weight);
1514  for(size_t i=0;i<pars.size();i++) {
1515  line.push_back(pars[i]);
1516  }
1517  int tempi=fill(pars,log_weight,line,dat);
1518  return tempi;
1519  }
1520 
1521  /** \brief For each walker, record the last row in the table which
1522  corresponds to an accept
1523  */
1524  std::vector<int> walker_accept_rows;
1525 
1526  /** \brief For each walker, record the last row in the table which
1527  corresponds to an reject
1528  */
1529  std::vector<int> walker_reject_rows;
1530 
1531  /// Likelihood estimator
1533 
1534  /** \brief Initial write to HDF5 file
1535  */
1536  virtual void file_header(o2scl_hdf::hdf_file &hf) {
1537  return;
1538  }
1539 
1540  /** \brief A copy of the lower limits for HDF5 output
1541  */
1542  vec_t low_copy;
1543 
1544  /** \brief A copy of the upper limits for HDF5 output
1545  */
1546  vec_t high_copy;
1547 
1548  /** \brief Total number of MCMC acceptances over all threads at last
1549  file write() (default 0)
1550  */
1552 
1553  /** \brief Time at last
1554  file write() (default 0.0)
1555  */
1557 
1558  public:
1559 
1560  /** \brief If true, ensure sure walkers and OpenMP threads are
1561  written to the table with equal spacing between rows (default
1562  true)
1563  */
1565 
1566  /** \brief Iterations between file updates (default 0 for no file updates)
1567  */
1569 
1570  /** \brief Time between file updates (default 0.0 for no file updates)
1571  */
1573 
1574  /** \brief The number of tables to combine before I/O (default 1)
1575  */
1577 
1578  /** \brief Write MCMC tables to files
1579  */
1580  virtual void write_files(bool sync_write=false) {
1581 
1582  if (this->verbose>=2) {
1583  this->scr_out << "mcmc: Start write_files(). mpi_rank: "
1584  << this->mpi_rank << " mpi_size: "
1585  << this->mpi_size << " table_io_chunk: "
1586  << table_io_chunk << std::endl;
1587  }
1588 
1589  std::vector<o2scl::table_units<> > tab_arr;
1590  bool rank_sent=false;
1591 
1592 #ifdef O2SCL_MPI
1593  if (table_io_chunk>1) {
1594  if (this->mpi_rank%table_io_chunk==0) {
1595  // Parent ranks
1596  for(int i=0;i<table_io_chunk-1;i++) {
1597  int child=this->mpi_rank+i+1;
1598  if (child<this->mpi_size) {
1599  table_units<> t;
1600  tab_arr.push_back(t);
1601  o2scl_table_mpi_recv(child,tab_arr[tab_arr.size()-1]);
1602  }
1603  }
1604  } else {
1605  // Child ranks
1606  size_t parent=this->mpi_rank-(this->mpi_rank%table_io_chunk);
1607  o2scl_table_mpi_send(*table,parent);
1608  rank_sent=true;
1609  }
1610  }
1611 #endif
1612 
1613 #ifdef O2SCL_MPI
1614  // Ensure that multiple threads aren't writing to the
1615  // filesystem at the same time
1616  int tag=0, buffer=0;
1617  if (sync_write && this->mpi_size>1 &&
1618  this->mpi_rank>=table_io_chunk) {
1619  MPI_Recv(&buffer,1,MPI_INT,this->mpi_rank-table_io_chunk,
1620  tag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
1621  }
1622 #endif
1623 
1625  std::string fname=this->prefix+"_"+o2scl::itos(this->mpi_rank)+"_out";
1626  hf.open_or_create(fname);
1627 
1628  if (first_write==false) {
1629  hf.set_szt("max_iters",this->max_iters);
1630  hf.sets("prefix",this->prefix);
1631  hf.seti("aff_inv",this->aff_inv);
1632  hf.setd("step_fac",this->step_fac);
1633  hf.setd("ai_initial_step",this->ai_initial_step);
1634  hf.set_szt("n_warm_up",this->n_warm_up);
1635  hf.seti("user_seed",this->user_seed);
1636  hf.seti("mpi_rank",this->mpi_rank);
1637  hf.seti("mpi_size",this->mpi_size);
1638  hf.seti("verbose",this->verbose);
1639  hf.set_szt("max_bad_steps",this->max_bad_steps);
1640  hf.set_szt("n_walk",this->n_walk);
1641  hf.set_szt("n_walk_per_thread",this->n_walk_per_thread);
1642  hf.set_szt("n_chains_per_rank",this->n_chains_per_rank);
1643  hf.set_szt("n_threads",this->n_threads);
1644  hf.set_szt("n_params",this->n_params);
1645  hf.seti("always_accept",this->always_accept);
1646  hf.seti("allow_estimates",allow_estimates);
1647  hf.setd_vec_copy("low",this->low_copy);
1648  hf.setd_vec_copy("high",this->high_copy);
1649  file_header(hf);
1650  first_write=true;
1651  }
1652 
1653  hf.set_szt_vec("n_accept",this->n_accept);
1654  hf.set_szt_vec("n_reject",this->n_reject);
1655  hf.set_szt_arr2d_copy("ret_value_counts",this->ret_value_counts.size(),
1656  this->ret_value_counts[0].size(),
1657  this->ret_value_counts);
1658  hf.setd_arr2d_copy("initial_points",this->initial_points.size(),
1659  this->initial_points[0].size(),
1660  this->initial_points);
1661 
1662  hf.seti("n_tables",tab_arr.size()+1);
1663  if (rank_sent==false) {
1664  hdf_output(hf,*table,"markov_chain_0");
1665  }
1666  for(size_t i=0;i<tab_arr.size();i++) {
1667  std::string name=((std::string)"markov_chain_")+szttos(i+1);
1668  hdf_output(hf,tab_arr[i],name);
1669  }
1670 
1671  hf.close();
1672 
1673 #ifdef O2SCL_MPI
1674  if (sync_write && this->mpi_size>1 &&
1675  this->mpi_rank<this->mpi_size-1) {
1676  MPI_Send(&buffer,1,MPI_INT,this->mpi_rank+table_io_chunk,
1677  tag,MPI_COMM_WORLD);
1678  }
1679 #endif
1680 
1681  if (this->verbose>=2) {
1682  this->scr_out << "mcmc: Done write_files()." << std::endl;
1683  }
1684 
1685  return;
1686  }
1687 
1688  /// If true, allow estimates of the weight (default false)
1690 
1691  /** \brief If true, store MCMC rejections in the table
1692  */
1694 
1695  mcmc_para_table() {
1696  allow_estimates=false;
1697  table_io_chunk=1;
1698  file_update_iters=0;
1699  file_update_time=0.0;
1700  last_write_iters=0;
1701  store_rejects=false;
1702  table_sequence=true;
1703  }
1704 
1705  /// \name Basic usage
1706  //@{
1707  /** \brief Set the table names and units
1708  */
1709  virtual void set_names_units(std::vector<std::string> names,
1710  std::vector<std::string> units) {
1711  col_names=names;
1712  col_units=units;
1713  return;
1714  }
1715 
1716  /** \brief Read initial points from the last points recorded in file
1717  named \c fname
1718 
1719  The values of \ref o2scl::mcmc_para_base::n_walk and \ref
1720  o2scl::mcmc_para_base::n_threads, must be set to their correct values before
1721  calling this function. This function requires that a table is
1722  present in \c fname which stores parameters in a block of
1723  columns and has columns named \c mult, \c thread,
1724  \c walker, and \c log_wgt.
1725  */
1726  virtual void initial_points_file_last(std::string fname,
1727  size_t n_param_loc,
1728  size_t offset=5) {
1729 
1730  table=std::shared_ptr<o2scl::table_units<> >(new o2scl::table_units<>);
1731 
1732 #ifdef O2SCL_MPI
1733  // Ensure that multiple threads aren't reading from the
1734  // filesystem at the same time
1735  int tag=0, buffer=0;
1736  if (this->mpi_size>1 && this->mpi_rank>0) {
1737  MPI_Recv(&buffer,1,MPI_INT,this->mpi_rank-1,
1738  tag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
1739  }
1740 #endif
1741 
1743  hf.open(fname);
1744  std::string tname;
1745  hdf_input(hf,*table,tname);
1746  hf.close();
1747 
1748 #ifdef O2SCL_MPI
1749  if (this->mpi_size>1 && this->mpi_rank<this->mpi_size-1) {
1750  MPI_Send(&buffer,1,MPI_INT,this->mpi_rank+1,
1751  tag,MPI_COMM_WORLD);
1752  }
1753 #endif
1754 
1755  // Determine number of points
1756  size_t n_points=this->n_walk*this->n_threads;
1757 
1758  if (this->verbose>0) {
1759  std::cout << "Initial points: Finding last " << n_points
1760  << " points from file named "
1761  << fname << " ." << std::endl;
1762  }
1763 
1764  this->initial_points.resize(n_points);
1765 
1766  for(size_t it=0;it<this->n_threads;it++) {
1767  for(size_t iw=0;iw<this->n_walk;iw++) {
1768 
1769  // The combined walker/thread index
1770  size_t windex=it*this->n_walk+iw;
1771 
1772  bool found=false;
1773  for(int row=table->get_nlines()-1;row>=0 && found==false;row--) {
1774  if (table->get("walker",row)==iw &&
1775  table->get("thread",row)==it &&
1776  table->get("mult",row)>0.5) {
1777 
1778  found=true;
1779 
1780  std::cout << "Function initial_point_file_last():\n\tit: "
1781  << it << "," << this->mpi_rank
1782  << " iw: " << iw << " row: "
1783  << row << " log_wgt: " << table->get("log_wgt",row)
1784  << std::endl;
1785 
1786  // Copy the entries from this row into the initial_points object
1787  this->initial_points[windex].resize(n_param_loc);
1788  for(size_t ip=0;ip<n_param_loc;ip++) {
1789  this->initial_points[windex][ip]=table->get(ip+offset,row);
1790  }
1791  }
1792  }
1793  if (found==false) {
1794  O2SCL_ERR("Function initial_points_file_last() failed.",
1796  }
1797  }
1798  }
1799 
1800  return;
1801  }
1802 
1803  /** \brief Read initial points from file
1804  named \c fname, distributing across the chain if necessary
1805 
1806  The values of \ref o2scl::mcmc_para_base::n_walk and \ref
1807  o2scl::mcmc_para_base::n_threads, must be set to their correct values
1808  before calling this function. This function requires that a
1809  table is present in \c fname which stores parameters in a block
1810  of columns.
1811  */
1812  virtual void initial_points_file_dist(std::string fname,
1813  size_t n_param_loc,
1814  size_t offset=5) {
1815 
1816  table=std::shared_ptr<o2scl::table_units<> >(new o2scl::table_units<>);
1817 
1818 #ifdef O2SCL_MPI
1819  // Ensure that multiple threads aren't reading from the
1820  // filesystem at the same time
1821  int tag=0, buffer=0;
1822  if (this->mpi_size>1 && this->mpi_rank>0) {
1823  MPI_Recv(&buffer,1,MPI_INT,this->mpi_rank-1,
1824  tag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
1825  }
1826 #endif
1827 
1829  hf.open(fname);
1830  std::string tname;
1831  hdf_input(hf,*table,tname);
1832  hf.close();
1833 
1834 #ifdef O2SCL_MPI
1835  if (this->mpi_size>1 && this->mpi_rank<this->mpi_size-1) {
1836  MPI_Send(&buffer,1,MPI_INT,this->mpi_rank+1,
1837  tag,MPI_COMM_WORLD);
1838  }
1839 #endif
1840 
1841  // Determine number of points
1842  size_t n_points=this->n_walk*this->n_threads;
1843 
1844  if (this->verbose>0) {
1845  std::cout << "Initial points: Finding last " << n_points
1846  << " points from file named "
1847  << fname << " ." << std::endl;
1848  }
1849 
1850  this->initial_points.resize(n_points);
1851 
1852  size_t nlines=table->get_nlines();
1853  size_t decrement=nlines/n_points;
1854  if (decrement<1) decrement=1;
1855 
1856  int row=nlines-1;
1857  for(size_t k=0;k<n_points;k++) {
1858  row-=decrement;
1859  if (row<0) row=0;
1860 
1861  std::cout << "Function initial_point_file_dist():\n\trow: "
1862  << row << " log_wgt: " << table->get("log_wgt",row)
1863  << std::endl;
1864 
1865  // Copy the entries from this row into the initial_points object
1866  this->initial_points[k].resize(n_param_loc);
1867  for(size_t ip=0;ip<n_param_loc;ip++) {
1868  this->initial_points[k][ip]=table->get(ip+offset,row);
1869  }
1870 
1871  }
1872 
1873  return;
1874  }
1875 
1876  /** \brief Read initial points from the best points recorded in file
1877  named \c fname
1878 
1879  The values of \ref o2scl::mcmc_para_base::n_walk and \ref
1880  o2scl::mcmc_para_base::n_threads, must be set to their correct values
1881  before calling this function. This function requires that a
1882  table is present in \c fname which stores parameters in a block
1883  of columns and contains a separate column named \c log_wgt .
1884  */
1885  virtual void initial_points_file_best(std::string fname,
1886  size_t n_param_loc,
1887  double thresh=1.0e-6,
1888  size_t offset=5) {
1889 
1890  table=std::shared_ptr<o2scl::table_units<> >(new o2scl::table_units<>);
1891 
1892 #ifdef O2SCL_MPI
1893  // Ensure that multiple threads aren't reading from the
1894  // filesystem at the same time
1895  int tag=0, buffer=0;
1896  if (this->mpi_size>1 && this->mpi_rank>0) {
1897  MPI_Recv(&buffer,1,MPI_INT,this->mpi_rank-1,
1898  tag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
1899  }
1900 #endif
1901 
1903  hf.open(fname);
1904  std::string tname;
1905  hdf_input(hf,*table,tname);
1906  hf.close();
1907 
1908 #ifdef O2SCL_MPI
1909  if (this->mpi_size>1 && this->mpi_rank<this->mpi_size-1) {
1910  MPI_Send(&buffer,1,MPI_INT,this->mpi_rank+1,
1911  tag,MPI_COMM_WORLD);
1912  }
1913 #endif
1914 
1915  // Determine number of points
1916  size_t n_points=this->n_walk*this->n_threads;
1917 
1918  if (this->verbose>0) {
1919  std::cout << "Initial points: Finding best " << n_points
1920  << " unique points from file named "
1921  << fname << " ." << std::endl;
1922  }
1923 
1924  typedef std::map<double,int,std::greater<double> > map_t;
1925  map_t m;
1926 
1927  // Sort by inserting into a map
1928  for(size_t k=0;k<table->get_nlines();k++) {
1929  m.insert(std::make_pair(table->get("log_wgt",k),k));
1930  }
1931 
1932  // Remove near duplicates. The map insert function will
1933  // just overwrite duplicate key entries, but we also
1934  // want to avoid near duplicates, so we have to search
1935  // manually for those.
1936  bool found;
1937  do {
1938  found=false;
1939  for(map_t::iterator mit=m.begin();mit!=m.end();mit++) {
1940  map_t::iterator mit2=mit;
1941  mit2++;
1942  if (mit2!=m.end()) {
1943  if (fabs(mit->first-mit2->first)<thresh) {
1944  if (this->verbose>0) {
1945  std::cout << "Removing duplicate weights: "
1946  << mit->first << " " << mit2->first << std::endl;
1947 
1948  }
1949  m.erase(mit2);
1950  mit=m.begin();
1951  found=true;
1952  }
1953  }
1954  }
1955  } while (found==true);
1956 
1957  // Check to see if we have enough
1958  if (m.size()<n_points) {
1959  O2SCL_ERR2("Could not find enough points in file in ",
1960  "mcmc_para::initial_points_file_best().",
1962  }
1963 
1964  // Copy the entries from this row into the initial_points object
1965  this->initial_points.resize(n_points);
1966  map_t::iterator mit=m.begin();
1967  for(size_t k=0;k<n_points;k++) {
1968  int row=mit->second;
1969  if (this->verbose>0) {
1970  std::cout << "Initial point " << k << " at row "
1971  << row << " has log_weight= "
1972  << table->get("log_wgt",row) << std::endl;
1973  }
1974  this->initial_points[k].resize(n_param_loc);
1975  for(size_t ip=0;ip<n_param_loc;ip++) {
1976  this->initial_points[k][ip]=table->get(ip+offset,row);
1977  }
1978  mit++;
1979  }
1980 
1981  return;
1982  }
1983 
1984  /** \brief Perform an MCMC simulation
1985 
1986  Perform an MCMC simulation over \c n_params parameters starting
1987  at initial point \c init, limiting the parameters to be between
1988  \c low and \c high, using \c func as the objective function and
1989  calling the measurement function \c meas at each MC point.
1990  */
1991  virtual int mcmc(size_t n_params_local,
1992  vec_t &low, vec_t &high, std::vector<func_t> &func,
1993  std::vector<fill_t> &fill) {
1994 
1995  n_params=n_params_local;
1996  low_copy=low;
1997  high_copy=high;
1998 
1999  first_write=false;
2000 
2001  // Set number of threads (this is done in the child as well, but
2002  // we need this number to set up the vector of measure functions
2003  // below).
2004 #ifdef O2SCL_OPENMP
2005  omp_set_num_threads(this->n_threads);
2006 #else
2007  this->n_threads=1;
2008 #endif
2009 
2010  // Setup the vector of measure functions
2011  std::vector<internal_measure_t> meas(this->n_threads);
2012  for(size_t it=0;it<this->n_threads;it++) {
2013  meas[it]=std::bind
2014  (std::mem_fn<int(const vec_t &,double,size_t,int,bool,
2015  data_t &, size_t, fill_t &)>
2016  (&mcmc_para_table::add_line),this,std::placeholders::_1,
2017  std::placeholders::_2,std::placeholders::_3,
2018  std::placeholders::_4,std::placeholders::_5,
2019  std::placeholders::_6,it,std::ref(fill[it]));
2020  }
2021 
2022  return parent_t::mcmc(n_params,low,high,func,meas);
2023  }
2024 
2025  /** \brief Get the output table
2026  */
2027  std::shared_ptr<o2scl::table_units<> > get_table() {
2028  return table;
2029  }
2030 
2031  /** \brief Set the output table
2032  */
2033  void set_table(std::shared_ptr<o2scl::table_units<> > &t) {
2034  table=t;
2035  return;
2036  }
2037 
2038  /** \brief Determine the chain sizes
2039 
2040  \future This algorithm could be improved by started from the end
2041  of the table and going backwards instead of starting from the
2042  front of the table and going forwards.
2043  */
2044  void get_chain_sizes(std::vector<size_t> &chain_sizes) {
2045 
2046  size_t ntot=this->n_threads*this->n_walk;
2047  chain_sizes.resize(ntot);
2048 
2049  for(size_t it=0;it<this->n_threads;it++) {
2050  for(size_t iw=0;iw<this->n_walk;iw++) {
2051  size_t ix=it*this->n_walk+iw;
2052  size_t istart=ix;
2053  chain_sizes[ix]=0;
2054  for(size_t j=istart;j<table->get_nlines();j+=ntot) {
2055  if (table->get("mult",j)>0.5) chain_sizes[ix]++;
2056  }
2057  }
2058  }
2059 
2060  return;
2061  }
2062 
2063  /** \brief A measurement function which adds the point to the
2064  table
2065  */
2066  virtual int add_line(const vec_t &pars, double log_weight,
2067  size_t walker_ix, int func_ret,
2068  bool mcmc_accept, data_t &dat,
2069  size_t i_thread, fill_t &fill) {
2070 
2071  // The combined walker/thread index
2072  size_t windex=i_thread*this->n_walk+walker_ix;
2073 
2074  // The total number of walkers * threads
2075  size_t ntot=this->n_threads*this->n_walk;
2076 
2077  int ret_value=o2scl::success;
2078 
2079  // Determine the next row
2080  int next_row;
2081  if ((mcmc_accept || store_rejects) && walker_accept_rows[windex]<0) {
2082  next_row=windex;
2083  } else {
2084  if (table_sequence) {
2085  if (walker_accept_rows[windex]>walker_reject_rows[windex]) {
2086  next_row=walker_accept_rows[windex]+ntot;
2087  } else {
2088  next_row=walker_reject_rows[windex]+ntot;
2089  }
2090  } else {
2091  if (walker_accept_rows[windex]>walker_reject_rows[windex]) {
2092  next_row=walker_accept_rows[windex]+1;
2093  } else {
2094  next_row=walker_reject_rows[windex]+1;
2095  }
2096  while (next_row<((int)table->get_nlines()) &&
2097  table->get("mult",next_row)>0.0) {
2098  next_row++;
2099  }
2100  }
2101  }
2102 
2103  // Make sure only one thread writes to the table at a time by
2104  // making this next region 'critical'. This is required because
2105  // writes to the table may require resizing the table structure.
2106 #ifdef O2SCL_OPENMP
2107 #pragma omp critical (o2scl_mcmc_para_table_add_line)
2108 #endif
2109  {
2110 
2111  // If there's not enough space in the table for this iteration,
2112  // then create it. There is not enough space if any of the
2113  // walker_accept_rows array entries is -1, if we have an
2114  // acceptance but there isn't room to store it, or if
2115  // we have a rejection and there isn't room to store it.
2116  if (next_row>=((int)table->get_nlines())) {
2117  size_t istart=table->get_nlines();
2118  // Create enough space
2119  table->set_nlines(table->get_nlines()+ntot);
2120  // Now additionally initialize the first four colums
2121  for(size_t j=0;j<this->n_threads;j++) {
2122  for(size_t i=0;i<this->n_walk;i++) {
2123  table->set("rank",istart+j*this->n_walk+i,this->mpi_rank);
2124  table->set("thread",istart+j*this->n_walk+i,j);
2125  table->set("walker",istart+j*this->n_walk+i,i);
2126  table->set("mult",istart+j*this->n_walk+i,0.0);
2127  table->set("log_wgt",istart+j*this->n_walk+i,0.0);
2128  }
2129  }
2130  }
2131 
2132  // If needed, add the line to the next row
2133  if (func_ret==0 && (mcmc_accept || store_rejects)) {
2134 
2135  if (next_row>=((int)(table->get_nlines()))) {
2136  O2SCL_ERR("Not enough space in table.",o2scl::exc_esanity);
2137  }
2138 
2139  std::vector<double> line;
2140  int fret=fill_line(pars,log_weight,line,dat,walker_ix,fill);
2141 
2142  // For rejections, set the multiplier to -1.0 (it was set to
2143  // 1.0 in the fill_line() call above)
2144  if (store_rejects && mcmc_accept==false) {
2145  line[3]=-1.0;
2146  }
2147 
2148  if (fret!=o2scl::success) {
2149 
2150  // If we're done, we stop before adding the last point to the
2151  // table. This is important because otherwise the last line in
2152  // the table will always only have unit multiplicity, which
2153  // may or may not be correct.
2154  ret_value=this->mcmc_done;
2155 
2156  } else {
2157 
2158  // First, double check that the table has the right
2159  // number of columns
2160  if (line.size()!=table->get_ncolumns()) {
2161  std::cout << "line: " << line.size() << " columns: "
2162  << table->get_ncolumns() << std::endl;
2163  for(size_t k=0;k<table->get_ncolumns() || k<line.size();k++) {
2164  std::cout << k << ". ";
2165  if (k<table->get_ncolumns()) {
2166  std::cout << table->get_column_name(k) << " ";
2167  std::cout << table->get_unit(table->get_column_name(k)) << " ";
2168  }
2169  if (k<line.size()) std::cout << line[k] << " ";
2170  std::cout << std::endl;
2171  }
2172  O2SCL_ERR("Table misalignment in mcmc_para_table::add_line().",
2173  exc_einval);
2174  }
2175 
2176  // Set the row
2177  table->set_row(((size_t)next_row),line);
2178 
2179  // Verbose output
2180  if (this->verbose>=2) {
2181  this->scr_out << "mcmc: Setting data at row " << next_row
2182  << std::endl;
2183  for(size_t k=0;k<line.size();k++) {
2184  this->scr_out << k << ". ";
2185  this->scr_out << table->get_column_name(k) << " ";
2186  this->scr_out << table->get_unit(table->get_column_name(k));
2187  this->scr_out << " " << line[k] << std::endl;
2188  }
2189  }
2190 
2191  }
2192 
2193  // End of 'if (mcmc_accept || store_rejects)'
2194  }
2195 
2196  // If necessary, increment the multiplier on the previous point
2197  if (ret_value==o2scl::success && mcmc_accept==false) {
2198  double mult_old=table->get("mult",walker_accept_rows[windex]);
2199  table->set("mult",walker_accept_rows[windex],mult_old+1.0);
2200  if (this->verbose>=2) {
2201  this->scr_out << "mcmc: Updating mult of row "
2202  << walker_accept_rows[windex]
2203  << " from " << mult_old << " to "
2204  << mult_old+1.0 << std::endl;
2205  }
2206 
2207  }
2208 
2209  }
2210  // End of critical region
2211 
2212  // Increment row counters if necessary
2213  if (ret_value==o2scl::success) {
2214  if (mcmc_accept) {
2215  walker_accept_rows[windex]=next_row;
2216  } else if (store_rejects && func_ret==0) {
2217  walker_reject_rows[windex]=next_row;
2218  }
2219  }
2220 
2221  return ret_value;
2222  }
2223  //@}
2224 
2225  /** \brief Function to run after point evaluation and measurement steps
2226  */
2227  virtual void post_pointmeas() {
2228 
2229  // If necessary, output to files
2230  bool updated=false;
2231  if (file_update_iters>0) {
2232  size_t total_iters=0;
2233  for(size_t it=0;it<this->n_chains_per_rank;it++) {
2234  total_iters+=this->n_accept[it]+this->n_reject[it];
2235  }
2236  if (total_iters>=last_write_iters+file_update_iters) {
2237  if (this->verbose>=1) {
2238  this->scr_out << "mcmc: Writing to file. total_iters: "
2239  << total_iters << " file_update_iters: "
2240  << file_update_iters << " last_write_iters: "
2241  << last_write_iters << std::endl;
2242  }
2243  write_files(false);
2244  last_write_iters=total_iters;
2245  updated=true;
2246  }
2247  }
2248  if (updated==false && file_update_time>0.0) {
2249 #ifdef O2SCL_MPI
2250  double elapsed=MPI_Wtime()-last_write_time;
2251 #else
2252  double elapsed=time(0)-last_write_time;
2253 #endif
2254  if (elapsed>file_update_time) {
2255  if (this->verbose>=1) {
2256  this->scr_out << "mcmc: Writing to file. elapsed: "
2257  << elapsed << " file_update_time: "
2258  << file_update_time << " last_write_time: "
2259  << last_write_time << std::endl;
2260  }
2261  write_files(false);
2262 #ifdef O2SCL_MPI
2263  last_write_time=MPI_Wtime();
2264 #else
2265  last_write_time=time(0);
2266 #endif
2267  }
2268  }
2269 
2270  return;
2271  }
2272 
2273  /** \brief Perform cleanup after an MCMC simulation
2274  */
2275  virtual void mcmc_cleanup() {
2276 
2277  // This section removes empty rows at the end of the
2278  // table that were allocated but not used
2279  int i;
2280  bool done=false;
2281  for(i=table->get_nlines()-1;i>=0 && done==false;i--) {
2282  done=true;
2283  if (table->get("mult",i)<1.0e-10) {
2284  done=false;
2285  }
2286  }
2287  if (i+2<((int)table->get_nlines())) {
2288  table->set_nlines(i+2);
2289  }
2290 
2291  write_files(true);
2292 
2293  return parent_t::mcmc_cleanup();
2294  }
2295 
2296  /** \brief Compute autocorrelation coefficients
2297  */
2298  virtual void ac_coeffs(size_t ncols, ubmatrix &ac_coeffs) {
2299  std::vector<size_t> chain_sizes;
2300  get_chain_sizes(chain_sizes);
2301  size_t min_size=chain_sizes[0];
2302  for(size_t i=1;i<chain_sizes.size();i++) {
2303  if (chain_sizes[i]<min_size) min_size=chain_sizes[i];
2304  }
2305  size_t N_max=min_size/2;
2306  ac_coeffs.resize(ncols,N_max-1);
2307  for(size_t i=0;i<ncols;i++) {
2308  for(size_t ell=1;ell<N_max;ell++) {
2309  ac_coeffs(i,ell-1)=0.0;
2310  }
2311  }
2312  size_t n_tot=this->n_threads*this->n_walk;
2313  size_t table_row=0;
2314  size_t cstart=table->lookup_column("log_wgt")+1;
2315  for(size_t i=0;i<ncols;i++) {
2316  for(size_t j=0;j<this->n_threads;j++) {
2317  for(size_t k=0;k<this->n_walk;k++) {
2318  size_t tindex=j*this->n_walk+k;
2319  for(size_t ell=1;ell<N_max;ell++) {
2320  const double &x=(*table)[cstart+i][table_row];
2321  double mean=o2scl::vector_mean<const double *>
2322  (chain_sizes[tindex]+1,&x);
2323  ac_coeffs(i,ell-1)+=o2scl::vector_lagk_autocorr
2324  <const double *>(chain_sizes[tindex]+1,&x,ell,mean);
2325  }
2326  table_row+=chain_sizes[tindex]+1;
2327  }
2328  }
2329  for(size_t ell=1;ell<N_max;ell++) {
2330  ac_coeffs(i,ell-1)/=((double)n_tot);
2331  }
2332  }
2333  return;
2334  }
2335 
2336  /** \brief Compute autocorrelation lengths
2337  */
2338  virtual void ac_lengths(size_t ncols, ubmatrix &ac_coeffs_cols,
2339  ubvector &ac_lengths) {
2340  size_t N_max=ac_coeffs_cols.size2();
2341  ac_lengths.resize(ncols);
2342  for(size_t icol=0;icol<ncols;icol++) {
2343  std::vector<double> tau(N_max);
2344  for(size_t i=5;i<N_max;i++) {
2345  double sum=0.0;
2346  for(size_t j=0;j<i;j++) {
2347  sum+=ac_coeffs_cols(icol,j);
2348  }
2349  tau[i]=1.0+2.0*sum;
2350  std::cout << tau[i] << " " << ((double)i)/5.0 << std::endl;
2351  }
2352  std::cout << std::endl;
2353  }
2354  return;
2355  }
2356 
2357  /** \brief Reorder the table by thread and walker index
2358  */
2359  virtual void reorder_table() {
2360 
2361  // Create a new table
2362  std::shared_ptr<o2scl::table_units<> > table2=
2363  std::shared_ptr<o2scl::table_units<> >(new o2scl::table_units<>);
2364 
2365  for(size_t i=0;i<this->n_threads;i++) {
2366  for(size_t j=0;j<this->n_walk;j++) {
2367  std::string func=std::string("abs(walker-")+o2scl::szttos(j)+
2368  ")<0.1 && abs(thread-"+o2scl::szttos(i)+")<0.1";
2369  table->copy_rows(func,*table2);
2370  }
2371  }
2372 
2373  return;
2374  }
2375 
2376  /** \brief Reaverage the data into blocks of a fixed
2377  size in order to avoid autocorrelations
2378 
2379  \note The number of blocks \c n_blocks must be larger than the
2380  current table size. This function expects to find a column named
2381  "mult" which contains the multiplicity of each column, as is the
2382  case after a call to \ref mcmc_para_base::mcmc().
2383 
2384  This function is useful to remove autocorrelations to the table
2385  so long as the autocorrelation length is shorter than the block
2386  size. This function does not compute the autocorrelation length
2387  to check that this is the case.
2388  */
2389  void reblock(size_t n_blocks) {
2390 
2391  for(size_t it=0;it<this->n_threads;it++) {
2392 
2393  size_t n=table->get_nlines();
2394  if (n_blocks>n) {
2395  O2SCL_ERR2("Cannot reblock. Not enough data in ",
2396  "mcmc_para_table::reblock().",o2scl::exc_einval);
2397  }
2398  size_t n_block=n/n_blocks;
2399  size_t m=table->get_ncolumns();
2400  for(size_t j=0;j<n_blocks;j++) {
2401  double mult=0.0;
2402  ubvector dat(m);
2403  for(size_t i=0;i<m;i++) {
2404  dat[i]=0.0;
2405  }
2406  for(size_t k=j*n_block;k<(j+1)*n_block;k++) {
2407  mult+=(*table)["mult"][k];
2408  for(size_t i=1;i<m;i++) {
2409  dat[i]+=(*table)[i][k]*(*table)["mult"][k];
2410  }
2411  }
2412  table->set("mult",j,mult);
2413  for(size_t i=1;i<m;i++) {
2414  dat[i]/=mult;
2415  table->set(i,j,dat[i]);
2416  }
2417  }
2418  table->set_nlines(n_blocks);
2419 
2420  }
2421 
2422  return;
2423  }
2424 
2425  };
2426 
2427  /** \brief MCMC class with a command-line interface
2428 
2429  This class forms the basis of the MCMC used in the Bayesian
2430  analysis of neutron star mass and radius in
2431  http://github.com/awsteiner/bamr .
2432 
2433  */
2434  template<class func_t, class fill_t, class data_t, class vec_t=ubvector>
2435  class mcmc_para_cli : public mcmc_para_table<func_t,fill_t,
2436  data_t,vec_t> {
2437 
2438  protected:
2439 
2441 
2442  /// \name Parameter objects for the 'set' command
2443  //@{
2444  o2scl::cli::parameter_double p_step_fac;
2445  o2scl::cli::parameter_size_t p_n_warm_up;
2446  o2scl::cli::parameter_int p_user_seed;
2447  o2scl::cli::parameter_size_t p_max_bad_steps;
2449  o2scl::cli::parameter_bool p_aff_inv;
2450  o2scl::cli::parameter_bool p_table_sequence;
2451  o2scl::cli::parameter_bool p_store_rejects;
2452  o2scl::cli::parameter_double p_max_time;
2453  o2scl::cli::parameter_size_t p_max_iters;
2454  //o2scl::cli::parameter_int p_max_chain_size;
2455  o2scl::cli::parameter_size_t p_file_update_iters;
2456  o2scl::cli::parameter_double p_file_update_time;
2457  //o2scl::cli::parameter_bool p_output_meas;
2459  o2scl::cli::parameter_int p_verbose;
2460  //@}
2461 
2462  /** \brief Initial write to HDF5 file
2463  */
2464  virtual void file_header(o2scl_hdf::hdf_file &hf) {
2465  hf.sets_vec("cl_args",this->cl_args);
2466  return;
2467  }
2468 
2469  public:
2470 
2471  /** \brief The arguments sent to the command-line
2472  */
2473  std::vector<std::string> cl_args;
2474 
2475  /// \name Customization functions
2476  //@{
2477  /** \brief Set up the 'cli' object
2478 
2479  This function just adds the four commands and the 'set' parameters
2480  */
2481  virtual void setup_cli(cli &cl) {
2482 
2483  // ---------------------------------------
2484  // Set commands/options
2485 
2486  /*
2487  static const size_t nopt=1;
2488  o2scl::comm_option_s options[nopt]={
2489  {'i',"initial-point","Set the starting point in the parameter space",
2490  1,-1,"<mode> [...]",
2491  ((std::string)"Mode can be one of 'best', 'last', 'N', or ")+
2492  "'values'. If mode is 'best', then it uses the point with the "+
2493  "largest weight and the second argument specifies the file. If "+
2494  "mode is 'last' then it uses the last point and the second "+
2495  "argument specifies the file. If mode is 'N' then it uses the Nth "+
2496  "point, the second argument specifies the value of N and the third "+
2497  "argument specifies the file. If mode is 'values', then the "+
2498  "remaining arguments specify all the parameter values. On the "+
2499  "command-line, enclose negative values in quotes and parentheses, "+
2500  "i.e. \"(-1.00)\" to ensure they do not get confused with other "+
2501  "options.",new o2scl::comm_option_mfptr<mcmc_para_cli>
2502  (this,&mcmc_para_cli::set_initial_point),
2503  o2scl::cli::comm_option_both}
2504  {'s',"hastings","Specify distribution for M-H step",
2505  1,1,"<filename>",
2506  ((string)"Desc. ")+"Desc2.",
2507  new comm_option_mfptr<mcmc_mpi>(this,&mcmc_mpi::hastings),
2508  cli::comm_option_both}
2509  };
2510  this->cl.set_comm_option_vec(nopt,options);
2511  */
2512 
2513  p_file_update_iters.s=&this->file_update_iters;
2514  p_file_update_iters.help=((std::string)"Number of MCMC successes ")+
2515  "between file updates (default 0 for no file updates).";
2516  cl.par_list.insert(std::make_pair("file_update_iters",
2517  &p_file_update_iters));
2518 
2519  p_file_update_time.d=&this->file_update_time;
2520  p_file_update_time.help=((std::string)"Time ")+
2521  "between file updates (default 0.0 for no file updates).";
2522  cl.par_list.insert(std::make_pair("file_update_time",
2523  &p_file_update_time));
2524 
2525  /*
2526  p_max_chain_size.i=&this->max_chain_size;
2527  p_max_chain_size.help=((std::string)"Maximum Markov chain size ")+
2528  "(default 10000).";
2529  cl.par_list.insert(std::make_pair("max_chain_size",
2530  &p_max_chain_size));
2531  */
2532 
2533  p_max_time.d=&this->max_time;
2534  p_max_time.help=((std::string)"Maximum run time in seconds ")+
2535  "(default 86400 sec or 1 day).";
2536  cl.par_list.insert(std::make_pair("max_time",&p_max_time));
2537 
2538  p_max_iters.s=&this->max_iters;
2539  p_max_iters.help=((std::string)"If non-zero, limit the number of ")+
2540  "iterations to be less than the specified number (default zero).";
2541  cl.par_list.insert(std::make_pair("max_iters",&p_max_iters));
2542 
2543  p_prefix.str=&this->prefix;
2544  p_prefix.help="Output file prefix (default 'mcmc\').";
2545  cl.par_list.insert(std::make_pair("prefix",&p_prefix));
2546 
2547  /*
2548  p_output_meas.b=&this->output_meas;
2549  p_output_meas.help=((std::string)"If true, output next point ")+
2550  "to the '_scr' file before calling TOV solver (default true).";
2551  cl.par_list.insert(std::make_pair("output_meas",&p_output_meas));
2552  */
2553 
2554  p_step_fac.d=&this->step_fac;
2555  p_step_fac.help=((std::string)"MCMC step factor. The step size for ")+
2556  "each variable is taken as the difference between the high and low "+
2557  "limits divided by this factor (default 10.0). This factor can "+
2558  "be increased if the acceptance rate is too small, but care must "+
2559  "be taken, e.g. if the conditional probability is multimodal. If "+
2560  "this step size is smaller than 1.0, it is reset to 1.0 .";
2561  cl.par_list.insert(std::make_pair("step_fac",&p_step_fac));
2562 
2563  p_n_warm_up.s=&this->n_warm_up;
2564  p_n_warm_up.help=((std::string)"Minimum number of warm up iterations ")+
2565  "(default 0).";
2566  cl.par_list.insert(std::make_pair("n_warm_up",&p_n_warm_up));
2567 
2568  p_verbose.i=&this->verbose;
2569  p_verbose.help=((std::string)"MCMC verbosity parameter ")+
2570  "(default 0).";
2571  cl.par_list.insert(std::make_pair("mcmc_verbose",&p_verbose));
2572 
2573  p_max_bad_steps.s=&this->max_bad_steps;
2574  p_max_bad_steps.help=((std::string)"Maximum number of bad steps ")+
2575  "(default 1000).";
2576  cl.par_list.insert(std::make_pair("max_bad_steps",&p_max_bad_steps));
2577 
2578  p_n_walk.s=&this->n_walk;
2579  p_n_walk.help=((std::string)"Number of walkers ")+
2580  "(default 1).";
2581  cl.par_list.insert(std::make_pair("n_walk",&p_n_walk));
2582 
2583  p_user_seed.i=&this->user_seed;
2584  p_user_seed.help=((std::string)"Seed for multiplier for random ")+
2585  "number generator. If zero is given (the default), then mcmc() "+
2586  "uses time(0) to generate a random seed.";
2587  cl.par_list.insert(std::make_pair("user_seed",&p_user_seed));
2588 
2589  p_aff_inv.b=&this->aff_inv;
2590  p_aff_inv.help=((std::string)"If true, then use affine-invariant ")+
2591  "sampling (default false).";
2592  cl.par_list.insert(std::make_pair("aff_inv",&p_aff_inv));
2593 
2594  p_table_sequence.b=&this->table_sequence;
2595  p_table_sequence.help=((std::string)"If true, then ensure equal spacing ")+
2596  "between walkers (default true).";
2597  cl.par_list.insert(std::make_pair("table_sequence",&p_table_sequence));
2598 
2599  p_store_rejects.b=&this->store_rejects;
2600  p_store_rejects.help=((std::string)"If true, then store MCMC rejections ")+
2601  "(default false).";
2602  cl.par_list.insert(std::make_pair("store_rejects",&p_store_rejects));
2603 
2604  return;
2605  }
2606 
2607  };
2608 
2609  // End of namespace
2610 }
2611 
2612 #endif
void get_chain_sizes(std::vector< size_t > &chain_sizes)
Determine the chain sizes.
Definition: mcmc_para.h:2044
double vector_lagk_autocorr(size_t n, const vec_t &data, size_t k, double mean)
Lag-k autocorrelation.
Definition: vec_stats.h:705
std::shared_ptr< o2scl::table_units<> > table
Main data table for Markov chain.
Definition: mcmc_para.h:1437
virtual void post_pointmeas()
Function to run after point evaluation and measurement steps.
Definition: mcmc_para.h:220
std::vector< data_t > data_arr
Data array.
Definition: mcmc_para.h:168
int setd_vec_copy(std::string name, const vec_t &v)
Set vector dataset named name with v.
Definition: hdf_file.h:377
std::string prefix
Prefix for output filenames (default "mcmc")
Definition: mcmc_para.h:268
static const int mcmc_done
Integer to indicate completion.
Definition: mcmc_para.h:271
std::vector< std::string > cl_args
The arguments sent to the command-line.
Definition: mcmc_para.h:2473
Double parameter for o2scl::cli.
Definition: cli.h:299
double ai_initial_step
Initial step fraction for affine-invariance sampling walkers (default 0.1)
Definition: mcmc_para.h:347
bool store_rejects
If true, store MCMC rejections in the table.
Definition: mcmc_para.h:1693
Multi-dimensional interpolation by inverse distance weighting.
Definition: interpm_idw.h:92
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
int set_szt_vec(std::string name, const std::vector< size_t > &v)
Set vector dataset named name with v.
size_t n_warm_up
Number of warm up steps (successful steps not iterations) (default 0)
Definition: mcmc_para.h:308
int * i
Parameter.
Definition: cli.h:340
A generic MCMC simulation class writing data to a o2scl::table_units object.
Definition: mcmc_para.h:1414
bool allow_estimates
If true, allow estimates of the weight (default false)
Definition: mcmc_para.h:1689
bool err_nonconv
If true, call the error handler if msolve() or msolve_de() does not converge (default true) ...
Definition: mcmc_para.h:338
int verbose
Output control (default 0)
Definition: mcmc_para.h:319
std::vector< vec_t > current
Current points in parameter space for each walker and each OpenMP thread.
Definition: mcmc_para.h:160
double last_write_time
Time at last file write() (default 0.0)
Definition: mcmc_para.h:1556
size_t file_update_iters
Iterations between file updates (default 0 for no file updates)
Definition: mcmc_para.h:1568
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.
String parameter for o2scl::cli.
Definition: cli.h:253
bool aff_inv
If true, use affine-invariant Monte Carlo.
Definition: mcmc_para.h:296
size_t max_iters
If non-zero, the maximum number of MCMC iterations (default 0)
Definition: mcmc_para.h:255
sanity check failed - shouldn&#39;t happen
Definition: err_hnd.h:65
std::vector< bool > switch_arr
Data switch array for each walker and each OpenMP thread.
Definition: mcmc_para.h:175
std::vector< std::string > col_names
Column names.
Definition: mcmc_para.h:1428
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.
Definition: mcmc_para.h:405
invalid argument supplied by user
Definition: err_hnd.h:59
virtual void ac_lengths(size_t ncols, ubmatrix &ac_coeffs_cols, ubvector &ac_lengths)
Compute autocorrelation lengths.
Definition: mcmc_para.h:2338
vec_t high_copy
A copy of the upper limits for HDF5 output.
Definition: mcmc_para.h:1546
std::vector< std::vector< size_t > > ret_value_counts
Return value counters, one vector independent chain.
Definition: mcmc_para.h:179
size_t n_chains_per_rank
Number of fully independent chains in each MPI rank.
Definition: mcmc_para.h:235
void seti(std::string name, int i)
Set an integer named name to value i.
std::vector< size_t > curr_walker
Index of the current walker.
Definition: mcmc_para.h:231
MCMC class with a command-line interface.
Definition: mcmc_para.h:2435
double step_fac
Stepsize factor (default 10.0)
Definition: mcmc_para.h:299
int sets_vec(std::string name, const std::vector< std::string > &s)
Set a vector of strings named name.
double file_update_time
Time between file updates (default 0.0 for no file updates)
Definition: mcmc_para.h:1572
size_t last_write_iters
Total number of MCMC acceptances over all threads at last file write() (default 0) ...
Definition: mcmc_para.h:1551
virtual void unset_proposal()
Go back to random-walk Metropolis with a uniform distribution.
Definition: mcmc_para.h:1360
std::function< int(const vec_t &, double, size_t, int, bool, data_t &)> internal_measure_t
Measurement functor type for the parent.
Definition: mcmc_para.h:1422
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.
Definition: mcmc_para.h:1726
bool always_accept
If true, accept all steps.
Definition: mcmc_para.h:342
std::vector< ubvector > initial_points
Initial points in parameter space.
Definition: mcmc_para.h:394
bool first_write
If true, the HDF5 I/O initial info has been written to the file (set by mcmc() )
Definition: mcmc_para.h:1442
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.
Definition: mcmc_para.h:1325
generic failure
Definition: err_hnd.h:61
int table_io_chunk
The number of tables to combine before I/O (default 1)
Definition: mcmc_para.h:1576
Configurable command-line interface.
Definition: cli.h:230
std::string help
Help description.
Definition: cli.h:242
size_t n_walk
Number of walkers for affine-invariant MC or 1 otherwise (default 1)
Definition: mcmc_para.h:329
virtual void mcmc_cleanup()
Perform cleanup after an MCMC simulation.
Definition: mcmc_para.h:2275
void sets(std::string name, std::string s)
Set a string named name to value s.
virtual int add_line(const vec_t &pars, double log_weight, size_t walker_ix, int func_ret, bool mcmc_accept, data_t &dat, size_t i_thread, fill_t &fill)
A measurement function which adds the point to the table.
Definition: mcmc_para.h:2066
std::vector< std::string > col_units
Column units.
Definition: mcmc_para.h:1431
void set_proposal(prob_vec_t &pv)
Set the proposal distribution.
Definition: mcmc_para.h:1347
double mpi_start_time
The MPI starting time (defaults to 0.0)
Definition: mcmc_para.h:245
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.
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.
Definition: mcmc_para.h:1495
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.
Definition: mcmc_para.h:1885
int mpi_size
The MPI number of processors.
Definition: mcmc_para.h:136
bool table_sequence
If true, ensure sure walkers and OpenMP threads are written to the table with equal spacing between r...
Definition: mcmc_para.h:1564
double max_time
Time in seconds (default is 0)
Definition: mcmc_para.h:264
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
std::string dtos(double x, int prec=6, bool auto_prec=false)
Convert a double to a string.
int user_seed
If non-zero, use as the seed for the random number generator (default 0)
Definition: mcmc_para.h:316
String parameter for o2scl::cli.
Definition: cli.h:276
bool * b
Parameter.
Definition: cli.h:283
void reblock(size_t n_blocks)
Reaverage the data into blocks of a fixed size in order to avoid autocorrelations.
Definition: mcmc_para.h:2389
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.
Definition: hdf_file.h:674
std::vector< size_t > n_reject
The number of Metropolis steps which were rejected in each independent chain (summed over all walkers...
Definition: mcmc_para.h:290
virtual void set_names_units(std::vector< std::string > names, std::vector< std::string > units)
Set the table names and units.
Definition: mcmc_para.h:1709
std::string * str
Parameter.
Definition: cli.h:260
std::map< std::string, parameter *, std::less< std::string > > par_list
Parameter list.
Definition: cli.h:380
bool warm_up
If true, we are in the warm up phase.
Definition: mcmc_para.h:152
std::vector< size_t > n_accept
The number of Metropolis steps which were accepted in each independent chain (summed over all walkers...
Definition: mcmc_para.h:283
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
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.
Definition: mcmc_para.h:1991
bool pd_mode
If true, then use the user-specified proposal distribution.
Definition: mcmc_para.h:149
virtual void reorder_table()
Reorder the table by thread and walker index.
Definition: mcmc_para.h:2359
double * d
Parameter.
Definition: cli.h:310
size_t n_params
Number of parameters.
Definition: mcmc_para.h:1434
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.
Definition: hdf_file.h:470
size_t n_walk_per_thread
Number of walkers per thread (default 1)
Definition: mcmc_para.h:333
virtual void mcmc_cleanup()
Cleanup after the MCMC.
Definition: mcmc_para.h:200
std::vector< int > walker_reject_rows
For each walker, record the last row in the table which corresponds to an reject. ...
Definition: mcmc_para.h:1529
std::vector< o2scl::prob_cond_mdim< vec_t > * > prop_dist
Pointer to proposal distribution for each thread.
Definition: mcmc_para.h:146
Integer parameter for o2scl::cli.
Definition: cli.h:333
virtual int mcmc_init()
MCMC initialization function.
Definition: mcmc_para.h:1448
Data table table class with units.
Definition: table_units.h:42
virtual int mcmc_init()
Initializations before the MCMC.
Definition: mcmc_para.h:185
int mpi_rank
The MPI processor rank.
Definition: mcmc_para.h:133
virtual void file_header(o2scl_hdf::hdf_file &hf)
Initial write to HDF5 file.
Definition: mcmc_para.h:2464
static const int mcmc_skip
Integer to indicate rejection.
Definition: mcmc_para.h:274
virtual void setup_cli(cli &cl)
Set up the &#39;cli&#39; object.
Definition: mcmc_para.h:2481
size_t max_bad_steps
Maximum number of failed steps when generating initial points with affine-invariant sampling (default...
Definition: mcmc_para.h:324
vec_t low_copy
A copy of the lower limits for HDF5 output.
Definition: mcmc_para.h:1542
std::vector< int > walker_accept_rows
For each walker, record the last row in the table which corresponds to an accept. ...
Definition: mcmc_para.h:1524
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:100
virtual void best_point(vec_t &best, double w_best, data_t &dat)
Function to run for the best point.
Definition: mcmc_para.h:214
interpm_idw< double * > esti
Likelihood estimator.
Definition: mcmc_para.h:1532
int open(std::string fname, bool write_access=false, bool err_on_fail=true)
Open a file named fname.
size_t n_threads
Number of OpenMP threads.
Definition: mcmc_para.h:387
mcmc_para_base< func_t, internal_measure_t, data_t, vec_t > parent_t
Type of parent class.
Definition: mcmc_para.h:1425
std::string itos(int x)
Convert an integer to a string.
A generic MCMC simulation class.
Definition: mcmc_para.h:126
std::vector< rng_gsl > rg
Random number generators.
Definition: mcmc_para.h:143
virtual void file_header(o2scl_hdf::hdf_file &hf)
Initial write to HDF5 file.
Definition: mcmc_para.h:1536
void set_table(std::shared_ptr< o2scl::table_units<> > &t)
Set the output table.
Definition: mcmc_para.h:2033
std::string szttos(size_t x)
Convert a size_t to a string.
virtual void write_files(bool sync_write=false)
Write MCMC tables to files.
Definition: mcmc_para.h:1580
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.
Definition: mcmc_para.h:1812
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
std::shared_ptr< o2scl::table_units<> > get_table()
Get the output table.
Definition: mcmc_para.h:2027
std::ofstream scr_out
The screen output file.
Definition: mcmc_para.h:140
virtual void ac_coeffs(size_t ncols, ubmatrix &ac_coeffs)
Compute autocorrelation coefficients.
Definition: mcmc_para.h:2298
Success.
Definition: err_hnd.h:47
virtual void post_pointmeas()
Function to run after point evaluation and measurement steps.
Definition: mcmc_para.h:2227

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