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

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