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

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