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

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