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

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