mcmc_para.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2012-2017, 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 and
25  \ref o2scl::mcmc_para_table
26 */
27 #ifndef O2SCL_MCMC_PARA_H
28 #define O2SCL_MCMC_PARA_H
29 
30 #include <iostream>
31 #include <random>
32 
33 #ifdef O2SCL_OPENMP
34 #include <omp.h>
35 #endif
36 #ifdef O2SCL_MPI
37 #include <mpi.h>
38 #endif
39 
40 #include <boost/numeric/ublas/vector.hpp>
41 
42 #include <o2scl/uniform_grid.h>
43 #include <o2scl/table3d.h>
44 #include <o2scl/hdf_file.h>
45 #include <o2scl/exception.h>
46 #include <o2scl/prob_dens_func.h>
47 #include <o2scl/vector.h>
48 #include <o2scl/multi_funct.h>
49 #include <o2scl/interpm_idw.h>
50 #include <o2scl/vec_stats.h>
51 
52 namespace o2scl {
53 
55 
56  /** \brief A generic MCMC simulation class
57 
58  This class performs a Markov chain Monte Carlo simulation of a
59  user-specified function using OpenMP and/or MPI. Either the
60  Metropolis-Hastings algorithm with a user-specified proposal
61  distribution or the affine-invariant sampling method of Goodman
62  and Weare can be used.
63 
64  By default, the Metropolis-Hastings algorithm is executed with a
65  simple walk, with steps in each dimension of size \f$
66  (\mathrm{high} - \mathrm{low})/\mathrm{step\_fac} \f$ with the
67  denominator specified in \ref step_fac.
68 
69  The function type is a template type, \c func_t, which should
70  be of the form
71  \code
72  int f(size_t num_of_parameters, const vec_t &parameters,
73  double &log_pdf, data_t &dat)
74  \endcode
75  which computes \c log_pdf, the natural logarithm of the function
76  value, for any point in parameter space (any point between \c
77  low and \c high ).
78 
79  If the function being simulated returns \ref mcmc_skip then the
80  point is automatically rejected. After each acceptance or
81  rejection, a user-specified "measurement" function (of type \c
82  measure_t ) is called, which can be used to store the results.
83  In order to stop the simulation, either this function or the
84  probability distribution being simulated should return the value
85  \ref mcmc_done .
86 
87  A generic proposal distribution can be specified in \ref
88  set_proposal(). To go back to the default random walk method,
89  one can call the function \ref unset_proposal().
90 
91  If \ref aff_inv is set to true and the number of walkers, \ref
92  n_walk is set to a number larger than 1, then affine-invariant
93  sampling is used. For affine-invariant sampling, the variable
94  \ref step_fac represents the value of \f$ a \f$, the
95  limits of the distribution for \f$ z \f$.
96 
97  In order to store data at each point, the user can store this
98  data in any object of type \c data_t . If affine-invariant
99  sampling is used, then each chain has it's own data object. The
100  class keeps twice as many copies of these data object as would
101  otherwise be required, in order to avoid copying of data objects
102  in the case that the steps are accepted or rejected.
103 
104  \note This class is experimental.
105  */
106  template<class func_t, class measure_t,
107  class data_t, class vec_t=ubvector> class mcmc_para_base {
108 
109  protected:
110 
111  /// \name MPI properties
112  //@{
113  /// The MPI processor rank
114  int mpi_rank;
115 
116  /// The MPI number of processors
118 
119  /// The MPI starting time
121  //@}
122 
123  /** \brief Time in seconds (default is 0)
124  */
125  double max_time;
126 
127  /// The screen output file
128  std::ofstream scr_out;
129 
130  /// Random number generators
131  std::vector<rng_gsl> rg;
132 
133  /// Proposal distribution
135 
136  /// If true, then use the user-specified proposal distribution
137  bool pd_mode;
138 
139  /// If true, we are in the warm up phase
140  bool warm_up;
141 
142  /** \brief Current points in parameter space for each walker and
143  each OpenMP thread
144 
145  This is an array of size \ref n_threads times \ref n_walk initial
146  guesses, indexed by <tt>thread_index*n_walk+walker_index</tt> .
147  */
148  std::vector<vec_t> current;
149 
150  /** \brief Data array
151 
152  This is an array of size 2 times \ref n_threads times \ref
153  n_walk . The two copies of data objects are indexed by
154  <tt>i_copy*n_walk*n_threads+thread_index*n_walk+walker_index</tt>
155  .
156  */
157  std::vector<data_t> data_arr;
158 
159  /** \brief Data switch array for each walker and each OpenMP thread
160 
161  This is an array of size \ref n_threads times \ref n_walk initial
162  guesses, indexed by <tt>thread_index*n_walk+walker_index</tt> .
163  */
164  std::vector<bool> switch_arr;
165 
166  /** \brief Return value counters, one vector for each OpenMP thread
167  */
168  std::vector<std::vector<size_t> > ret_value_counts;
169 
170  /// \name Interface customization
171  //@{
172  /** \brief Initializations before the MCMC
173  */
174  virtual int mcmc_init() {
175  if (verbose>0) {
176  // Open main output file for this rank
177  scr_out.open((prefix+"_"+
178  o2scl::itos(mpi_rank)+"_scr").c_str());
179  scr_out.setf(std::ios::scientific);
180  }
181  return 0;
182  }
183 
184  /** \brief Cleanup after the MCMC
185  */
186  virtual void mcmc_cleanup() {
187  if (verbose>0) {
188  for(size_t it=0;it<n_threads;it++) {
189  scr_out << "mcmc (" << it << "): accept=" << n_accept[it]
190  << " reject=" << n_reject[it] << std::endl;
191  }
192  scr_out.close();
193  }
194  return;
195  }
196 
197  /** \brief Function to run for the best point
198  */
199  virtual void best_point(vec_t &best, double w_best, data_t &dat) {
200  return;
201  }
202  //@}
203 
204  /** \brief Index of the current walker
205 
206  This quantity has to be a vector because different threads
207  may have different values for the current walker during
208  the initialization phase for the affine sampling algorithm.
209  */
210  std::vector<size_t> curr_walker;
211 
212  public:
213 
214  /// If non-zero, the maximum number of MCMC iterations (default 0)
215  size_t max_iters;
216 
217  /** \brief Prefix for output filenames
218  */
219  std::string prefix;
220 
221  /// Integer to indicate completion
222  static const int mcmc_done=-10;
223 
224  /// Integer to indicate rejection
225  static const int mcmc_skip=-20;
226 
227  /// \name Output quantities
228  //@{
229  /** \brief The number of Metropolis steps which were accepted in
230  each thread (summed over all walkers)
231  */
232  std::vector<size_t> n_accept;
233 
234  /** \brief The number of Metropolis steps which were rejected in
235  each thread (summed over all walkers)
236  */
237  std::vector<size_t> n_reject;
238  //@}
239 
240  /// \name Settings
241  //@{
242  /// If true, use affine-invariant Monte Carlo
243  bool aff_inv;
244 
245  /// Stepsize factor (default 10.0)
246  double step_fac;
247 
248  /** \brief Number of warm up steps (successful steps not
249  iterations) (default 0)
250 
251  \note Not to be confused with <tt>warm_up</tt>, which is
252  a boolean local variable in some functions not an int.
253  */
254  size_t n_warm_up;
255 
256  /** \brief If non-zero, use as the seed for the random number
257  generator (default 0)
258  */
260 
261  /// Output control (default 0)
262  int verbose;
263 
264  /** \brief Maximum number of failed steps when generating initial points
265  with affine-invariant sampling (default 1000)
266  */
268 
269  /** \brief Number of walkers for affine-invariant MC or 1
270  otherwise (default 1)
271  */
272  size_t n_walk;
273 
274  /** \brief If true, call the error handler if msolve() or
275  msolve_de() does not converge (default true)
276  */
278 
279  /** \brief If true, accept all steps
280  */
282 
283  /** \brief Initial step fraction for affine-invariance sampling walkers
284  (default 0.1)
285  */
287  //@}
288 
289  mcmc_para_base() {
290  user_seed=0;
291  n_warm_up=0;
292 
293  // MC step parameters
294  aff_inv=false;
295  pd_mode=false;
296  step_fac=10.0;
297  n_walk=1;
298  err_nonconv=true;
299  verbose=0;
300  warm_up=false;
301  max_bad_steps=1000;
302 
303  prop_dist=0;
304  always_accept=false;
305  ai_initial_step=0.1;
306 
307  n_threads=1;
308 
309  // Initial values for MPI paramers
310  mpi_nprocs=1;
311  mpi_rank=0;
312  mpi_start_time=0.0;
313 
314 #ifdef O2SCL_MPI
315  // Get MPI rank, etc.
316  MPI_Comm_rank(MPI_COMM_WORLD,&this->mpi_rank);
317  MPI_Comm_size(MPI_COMM_WORLD,&this->mpi_nprocs);
318 #endif
319 
320  prefix="mcmc";
321  max_time=0.0;
322  max_iters=0;
323  }
324 
325  /// Number of OpenMP threads
326  size_t n_threads;
327 
328  /** \brief Initial points in parameter space
329 
330  To fully specify all of the initial points, this should be
331  a vector of size \ref n_walk times \ref n_threads .
332  */
333  std::vector<ubvector> initial_points;
334 
335  /// \name Basic usage
336  //@{
337  /** \brief Perform an MCMC simulation
338 
339  Perform an MCMC simulation over \c nparams parameters starting
340  at initial point \c init, limiting the parameters to be between
341  \c low and \c high, using \c func as the objective function and
342  calling the measurement function \c meas at each MC point.
343  */
344  virtual int mcmc(size_t nparams, vec_t &low, vec_t &high,
345  std::vector<func_t> &func, std::vector<measure_t> &meas) {
346 
347  if (initial_points.size()==0) {
348  // Setup initial guess if not specified
349  initial_points.resize(1);
350  for(size_t k=0;k<nparams;k++) {
351  initial_points[0][k]=(low[k]+high[k])/2.0;
352  }
353  } else {
354  // If initial points are specified, make sure they're within
355  // the user-specified limits
356  for(size_t iip=0;iip<initial_points.size();iip++) {
357  for(size_t ipar=0;ipar<nparams;ipar++) {
358  if (initial_points[iip][ipar]<low[ipar] ||
359  initial_points[iip][ipar]>high[ipar]) {
360  O2SCL_ERR((((std::string)"Parameter ")+o2scl::szttos(ipar)+
361  " of "+o2scl::szttos(nparams)+" out of range (value="+
362  o2scl::dtos(initial_points[iip][ipar])+
363  ") in mcmc_base::mcmc().").c_str(),
365  }
366  }
367  }
368  }
369 
370  // Set number of threads
371 #ifdef O2SCL_OPENMP
372  omp_set_num_threads(n_threads);
373  n_threads=omp_get_num_threads();
374 #else
375  n_threads=1;
376 #endif
377 
378  // Set starting time
379 #ifdef O2SCL_MPI
380  mpi_start_time=MPI_Wtime();
381 #else
382  mpi_start_time=time(0);
383 #endif
384 
385  // Storage for return values from each thread
386  std::vector<int> func_ret(n_threads), meas_ret(n_threads);
387 
388  // Fix the number of walkers if it is too small
389  if (aff_inv) {
390  if (n_walk<=1) n_walk=2;
391  if (n_walk<n_threads) n_walk=n_threads;
392  }
393 
394  // Fix 'step_fac' if it's less than or equal to zero
395  if (step_fac<=0.0) {
396  if (aff_inv) step_fac=2.0;
397  else step_fac=10.0;
398  }
399 
400  // Set RNGs with a different seed for each thread and rank
401  rg.resize(n_threads);
402  unsigned long int seed=time(0);
403  if (this->user_seed!=0) {
404  seed=this->user_seed;
405  }
406  for(size_t it=0;it<n_threads;it++) {
407  seed*=(mpi_rank*n_threads+it+1);
408  rg[it].set_seed(seed);
409  }
410 
411  // Keep track of successful and failed MH moves
412  n_accept.resize(n_threads);
413  n_reject.resize(n_threads);
414  for(size_t it=0;it<n_threads;it++) {
415  n_accept[it]=0;
416  n_reject[it]=0;
417  }
418 
419  // Warm-up flag, not to be confused with 'n_warm_up', i.e. the
420  // number of warm_up iterations.
421  warm_up=true;
422  if (n_warm_up==0) warm_up=false;
423 
424  // Storage size required
425  size_t ssize=n_walk*n_threads;
426 
427  // Allocate current point and current weight
428  current.resize(ssize);
429  std::vector<double> w_current(ssize);
430  for(size_t i=0;i<ssize;i++) {
431  current[i].resize(nparams);
432  w_current[i]=0.0;
433  }
434 
435  // Allocate ret_value_counts and curr_walker
436  ret_value_counts.resize(n_threads);
437  curr_walker.resize(n_threads);
438 
439  // Initialize data and switch arrays
440  data_arr.resize(2*ssize);
441  switch_arr.resize(ssize);
442  for(size_t i=0;i<switch_arr.size();i++) switch_arr[i]=false;
443 
444  // Next point and next weight for each thread
445  std::vector<vec_t> next(n_threads);
446  for(size_t it=0;it<n_threads;it++) {
447  next[it].resize(nparams);
448  }
449  std::vector<double> w_next(n_threads);
450 
451  // Best point over all threads
452  vec_t best(nparams);
453  double w_best;
454 
455  // Generally, these flags are are true for any thread if func_ret
456  // or meas_ret is equal to mcmc_done.
457  std::vector<bool> mcmc_done_flag(n_threads);
458  for(size_t it=0;it<n_threads;it++) {
459  mcmc_done_flag[it]=false;
460  }
461 
462  // Proposal weight
463  double q_prop;
464 
465  // Run mcmc_init() function. The initial point, stored in
466  // current[0] can be modified by this function and the local
467  // variable 'init' is not accessible to the mcmc_init() function.
468  int init_ret=mcmc_init();
469  if (init_ret!=0) {
470  O2SCL_ERR("Function mcmc_init() failed in mcmc_base::mcmc().",
472  return init_ret;
473  }
474 
475  // Initial verbose output
476  if (verbose>=1) {
477  if (aff_inv) {
478  scr_out << "mcmc: Affine-invariant step, n_params="
479  << nparams << ", n_walk=" << n_walk
480  << ", n_threads=" << n_threads << ", n_ranks="
481  << mpi_nprocs << std::endl;
482  } else if (pd_mode==true) {
483  scr_out << "mcmc: With proposal distribution, n_params="
484  << nparams << ", n_threads=" << n_threads << ", n_ranks="
485  << mpi_nprocs << std::endl;
486  } else {
487  scr_out << "mcmc: Random-walk w/uniform dist., n_params="
488  << nparams << ", n_threads=" << n_threads << ", n_ranks="
489  << mpi_nprocs << std::endl;
490  }
491  }
492 
493  // ---------------------------------------------------
494  // Compute initial point and initial weights
495 
496  // Initial point and weights for stretch-move algorithm
497  if (aff_inv) {
498 
499 #ifdef O2SCL_OPENMP
500 #pragma omp parallel default(shared)
501 #endif
502  {
503 #ifdef O2SCL_OPENMP
504 #pragma omp for
505 #endif
506  for(size_t it=0;it<n_threads;it++) {
507 
508  // Initialize each walker in turn
509  for(curr_walker[it]=0;curr_walker[it]<n_walk &&
510  mcmc_done_flag[it]==false;curr_walker[it]++) {
511 
512  // Index in storage
513  size_t sindex=n_walk*it+curr_walker[it];
514 
515  size_t init_iters=0;
516  bool done=false;
517 
518  // If we already have a guess, try to use that
519  if (sindex<initial_points.size()) {
520 
521  // Copy from the initial points array
522  for(size_t ipar=0;ipar<nparams;ipar++) {
523  current[sindex][ipar]=initial_points[sindex][ipar];
524  }
525 
526  // Compute the weight
527  func_ret[it]=func[it](nparams,current[sindex],
528  w_current[sindex],data_arr[sindex]);
529 
530  if (func_ret[it]==mcmc_done) {
531  mcmc_done_flag[it]=true;
532  } else if (func_ret[it]==o2scl::success) {
533 
534  // If we have a good point, update ret_value_counts
535  // and call the measurement function
536  if (func_ret[it]>=0 &&
537  func_ret[it]<((int)ret_value_counts[it].size())) {
538  ret_value_counts[it][func_ret[it]]++;
539  }
540  meas_ret[it]=meas[it](current[sindex],w_current[sindex],
541  curr_walker[it],true,data_arr[sindex]);
542  if (meas_ret[it]==mcmc_done) {
543  mcmc_done_flag[it]=true;
544  }
545  done=true;
546  }
547  }
548 
549  while (!done && !mcmc_done_flag[it]) {
550 
551  // Make a perturbation from the initial point
552  for(size_t ipar=0;ipar<nparams;ipar++) {
553  do {
554  current[sindex][ipar]=
555  initial_points[sindex % initial_points.size()][ipar]+
556  (rg[it].random()*2.0-1.0)*
557  (high[ipar]-low[ipar])*ai_initial_step;
558  } while (current[sindex][ipar]>high[ipar] ||
559  current[sindex][ipar]<low[ipar]);
560  }
561 
562  // Compute the weight
563  func_ret[it]=func[it](nparams,current[sindex],
564  w_current[sindex],data_arr[sindex]);
565 
566  // ------------------------------------------------
567 
568  // Increment iteration count
569  init_iters++;
570 
571  if (func_ret[it]==mcmc_done) {
572  mcmc_done_flag[it]=true;
573  } else {
574  // If we have a good point, update ret_value_counts,
575  // call the measurement function and stop the loop
576  if (func_ret[it]==o2scl::success) {
577  if (func_ret[it]>=0 &&
578  func_ret[it]<((int)ret_value_counts[it].size())) {
579  ret_value_counts[it][func_ret[it]]++;
580  }
581  if (meas_ret[it]!=mcmc_done) {
582  meas_ret[it]=meas[it](current[sindex],w_current[sindex],
583  curr_walker[it],true,
584  data_arr[sindex]);
585  } else {
586  mcmc_done_flag[it]=true;
587  }
588  done=true;
589  } else if (init_iters>max_bad_steps) {
590  O2SCL_ERR2("Initial walkers failed in ",
591  "mcmc_para_base::mcmc().",o2scl::exc_einval);
592  }
593  }
594  }
595  }
596  }
597  }
598  // End of parallel region
599 
600  // Stop early if mcmc_done was returned
601  bool stop_early=false;
602  for(size_t it=0;it<n_threads;it++) {
603  if (mcmc_done_flag[it]==true) {
604  if (verbose>=1) {
605  scr_out << "mcmc (" << it << "): Returned mcmc_done "
606  << "(initial; ai)." << std::endl;
607  }
608  stop_early=true;
609  }
610  }
611  if (stop_early) {
612  mcmc_cleanup();
613  return 0;
614  }
615 
616  // Set initial values for best point
617  w_best=w_current[0];
618  size_t best_index=0;
619  for(size_t it=0;it<n_threads;it++) {
620  for(curr_walker[it]=0;curr_walker[it]<n_walk;curr_walker[it]++) {
621  size_t sindex=n_walk*it+curr_walker[it];
622  if (w_current[sindex]>w_current[0]) {
623  best_index=sindex;
624  w_best=w_current[sindex];
625  }
626  }
627  }
628  best=current[best_index];
629  best_point(best,w_best,data_arr[best_index]);
630 
631  // Verbose output
632  if (verbose>=2) {
633  for(size_t it=0;it<n_threads;it++) {
634  for(curr_walker[it]=0;curr_walker[it]<n_walk;curr_walker[it]++) {
635  size_t sindex=n_walk*it+curr_walker[it];
636  scr_out.precision(4);
637  scr_out << "mcmc (" << it << "): i_walk: ";
638  scr_out.width((int)(1.0+log10((double)(n_walk-1))));
639  scr_out << curr_walker[it] << " log_wgt: "
640  << w_current[sindex] << " (initial; ai)" << std::endl;
641  scr_out.precision(6);
642  }
643  }
644  }
645 
646  } else {
647 
648  // Note that this value is used (e.g. in
649  // mcmc_para_table::add_line() ) even if aff_inv is false, so we
650  // set it to zero here.
651  for(size_t it=0;it<n_threads;it++) {
652  curr_walker[it]=0;
653  }
654 
655  // Copy from the initial points array
656  for(size_t ipar=0;ipar<nparams;ipar++) {
657  current[0][ipar]=initial_points[0][ipar];
658  }
659 
660  // Initial point and weights without stretch-move
661 
662  func_ret[0]=func[0](nparams,current[0],w_current[0],data_arr[0]);
663  if (func_ret[0]==mcmc_done) {
664  if (verbose>=1) {
665  scr_out << "mcmc: Initial point returned mcmc_done."
666  << std::endl;
667  }
668  mcmc_cleanup();
669  return 0;
670  }
671  if (func_ret[0]!=o2scl::success) {
672  if (err_nonconv) {
673  O2SCL_ERR("Initial weight vanished in mcmc_para_base::mcmc().",
675  }
676  return 2;
677  }
678 
679 #ifdef O2SCL_OPENMP
680 #pragma omp parallel default(shared)
681 #endif
682  {
683 #ifdef O2SCL_OPENMP
684 #pragma omp for
685 #endif
686  for(size_t it=0;it<n_threads;it++) {
687  // Copy the results over from the initial point
688  if (it!=0) {
689  func_ret[it]=func_ret[0];
690  current[it]=current[0];
691  w_current[it]=w_current[0];
692  data_arr[it]=data_arr[0];
693  }
694  // Update the return value count
695  if (func_ret[it]>=0 &&
696  func_ret[it]<((int)ret_value_counts[it].size())) {
697  ret_value_counts[it][func_ret[it]]++;
698  }
699  // Call the measurement function
700  meas_ret[it]=meas[it](current[it],w_current[it],0,
701  true,data_arr[it]);
702  if (meas_ret[it]==mcmc_done) {
703  mcmc_done_flag[it]=true;
704  }
705  }
706  }
707  // End of parallel region
708 
709  // Stop early if mcmc_done was returned
710  bool stop_early=false;
711  for(size_t it=0;it<n_threads;it++) {
712  if (mcmc_done_flag[it]==true) {
713  if (verbose>=1) {
714  scr_out << "mcmc (" << it << "): Returned mcmc_done "
715  << "(initial)." << std::endl;
716  }
717  stop_early=true;
718  }
719  }
720  if (stop_early) {
721  mcmc_cleanup();
722  return 0;
723  }
724 
725  // Set initial values for best point
726  best=current[0];
727  w_best=w_current[0];
728  best_point(best,w_best,data_arr[0]);
729 
730  if (verbose>=2) {
731  scr_out.precision(4);
732  scr_out << "mcmc (0): "
733  << w_current[0] << " (initial)" << std::endl;
734  scr_out.precision(6);
735  }
736 
737  }
738 
739  if (verbose>=3) {
740  std::cout << "Press a key and type enter to continue. ";
741  char ch;
742  std::cin >> ch;
743  }
744 
745  // End of initial point and weight section
746  // ---------------------------------------------------
747 
748  bool main_done=false;
749  size_t mcmc_iters=0;
750 
751  while (!main_done) {
752 
753  // Walker to move (or zero when aff_inv is false)
754  std::vector<double> smove_z(n_threads);
755  for(size_t it=0;it<n_threads;it++) {
756  curr_walker[it]=0;
757  smove_z[it]=0.0;
758  // Choose walker to move (same for all threads)
759  if (aff_inv) {
760  curr_walker[it]=mcmc_iters % n_walk;
761  }
762  }
763 
764 #ifdef O2SCL_OPENMP
765 #pragma omp parallel default(shared)
766 #endif
767  {
768 #ifdef O2SCL_OPENMP
769 #pragma omp for
770 #endif
771  for(size_t it=0;it<n_threads;it++) {
772 
773  // ---------------------------------------------------
774  // Select next point
775 
776  if (aff_inv) {
777  // Choose jth walker
778  size_t ij;
779  do {
780  ij=((size_t)(rg[it].random()*((double)n_walk)));
781  } while (ij==curr_walker[it] || ij>=n_walk);
782 
783  // Select z
784  double p=rg[it].random();
785  double a=step_fac;
786  smove_z[it]=(1.0-2.0*p+2.0*a*p+p*p-2.0*a*p*p+a*a*p*p)/a;
787 
788  // Create new trial point
789  for(size_t i=0;i<nparams;i++) {
790  next[it][i]=current[n_walk*it+ij][i]+
791  smove_z[it]*(current[n_walk*it+curr_walker[it]][i]-
792  current[n_walk*it+ij][i]);
793  }
794 
795  } else if (pd_mode) {
796 
797  // Use proposal distribution and compute associated weight
798  (*prop_dist)(current[it],next[it]);
799  q_prop=prop_dist->log_pdf(current[it],next[it])-
800  prop_dist->log_pdf(next[it],current[it]);
801  if (!std::isfinite(q_prop)) {
802  O2SCL_ERR2("Proposal distribution not finite in ",
803  "mcmc_para_base::mcmc().",o2scl::exc_efailed);
804  }
805 
806  } else {
807 
808  // Uniform random-walk step
809  for(size_t k=0;k<nparams;k++) {
810  next[it][k]=current[it][k]+(rg[it].random()*2.0-1.0)*
811  (high[k]-low[k])/step_fac;
812  }
813 
814  }
815 
816  // ---------------------------------------------------
817  // Compute next weight
818 
819  func_ret[it]=o2scl::success;
820  // If the next point out of bounds, ensure that the point is rejected
821  for(size_t k=0;k<nparams;k++) {
822  if (next[it][k]<low[k] || next[it][k]>high[k]) {
823  func_ret[it]=mcmc_skip;
824  }
825  }
826  if (func_ret[it]!=mcmc_skip) {
827  if (switch_arr[n_walk*it+curr_walker[it]]==false) {
828  func_ret[it]=func[it](nparams,next[it],w_next[it],
829  data_arr[it*n_walk+curr_walker[it]+
830  n_walk*n_threads]);
831  } else {
832  func_ret[it]=func[it](nparams,next[it],w_next[it],
833  data_arr[it*n_walk+curr_walker[it]]);
834  }
835  if (func_ret[it]==mcmc_done) {
836  mcmc_done_flag[it]=true;
837  } else {
838  if (func_ret[it]>=0 &&
839  func_ret[it]<((int)ret_value_counts[it].size())) {
840  ret_value_counts[it][func_ret[it]]++;
841  }
842  }
843 
844  }
845  }
846  }
847  // End of parallel region
848 
849  // Handle verbose output serially
850  if (verbose>=1) {
851  for(size_t it=0;it<n_threads;it++) {
852  if (func_ret[it]==mcmc_done) {
853  scr_out << "mcmc (" << it << "): Returned mcmc_done."
854  << std::endl;
855  } else if (func_ret[it]==mcmc_skip && verbose>=3) {
856  scr_out << "mcmc (" << it
857  << "): Parameter(s) out of range: " << std::endl;
858  scr_out.setf(std::ios::showpos);
859  for(size_t k=0;k<nparams;k++) {
860  scr_out << k << " " << low[k] << " "
861  << next[it][k] << " " << high[k];
862  if (next[it][k]<low[k] || next[it][k]>high[k]) {
863  scr_out << " <-";
864  }
865  scr_out << std::endl;
866  }
867  scr_out.unsetf(std::ios::showpos);
868  } else if (func_ret[it]!=o2scl::success &&
869  func_ret[it]!=mcmc_skip) {
870  if (verbose>=2) {
871  scr_out << "mcmc (" << it << "): Function returned failure "
872  << func_ret[it] << " at point ";
873  for(size_t k=0;k<nparams;k++) {
874  scr_out << next[it][k] << " ";
875  }
876  scr_out << std::endl;
877  }
878  }
879  }
880  if (verbose>=3) {
881  std::cout << "Press a key and type enter to continue. ";
882  char ch;
883  std::cin >> ch;
884  }
885  }
886 
887 #ifdef O2SCL_OPENMP
888 #pragma omp parallel default(shared)
889 #endif
890  {
891 #ifdef O2SCL_OPENMP
892 #pragma omp for
893 #endif
894  for(size_t it=0;it<n_threads;it++) {
895 
896  // Index in storage
897  size_t sindex=n_walk*it+curr_walker[it];
898 
899  // ---------------------------------------------------
900  // Accept or reject
901 
902  bool accept=false;
903  if (always_accept && func_ret[it]==success) accept=true;
904 
905  if (func_ret[it]==o2scl::success) {
906  double r=rg[it].random();
907 
908  if (aff_inv) {
909  double ai_ratio=pow(smove_z[it],((double)nparams)-1.0)*
910  exp(w_next[it]-w_current[sindex]);
911  if (r<ai_ratio) {
912  accept=true;
913  }
914  } else if (pd_mode) {
915  if (r<exp(w_next[it]-w_current[sindex]+q_prop)) {
916  accept=true;
917  }
918  } else {
919  // Metropolis algorithm
920  if (r<exp(w_next[it]-w_current[sindex])) {
921  accept=true;
922  }
923  }
924 
925  // End of 'if (func_ret[it]==o2scl::success)'
926  }
927 
928  if (accept) {
929 
930  n_accept[it]++;
931 
932  // Store results from new point
933  if (!warm_up) {
934  if (switch_arr[sindex]==false) {
935  meas_ret[it]=meas[it](next[it],w_next[it],curr_walker[it],true,
936  data_arr[sindex+n_threads*n_walk]);
937  } else {
938  meas_ret[it]=meas[it](next[it],w_next[it],curr_walker[it],true,
939  data_arr[sindex]);
940  }
941  }
942 
943  // Prepare for next point
944  current[sindex]=next[it];
945  w_current[sindex]=w_next[it];
946  switch_arr[sindex]=!(switch_arr[sindex]);
947 
948  } else {
949 
950  // Point was rejected
951  n_reject[it]++;
952 
953  // Repeat measurement of old point
954  if (!warm_up) {
955  if (switch_arr[sindex]==false) {
956  meas_ret[it]=meas[it](current[sindex],
957  w_current[sindex],
958  curr_walker[it],false,data_arr[sindex]);
959  } else {
960  meas_ret[it]=meas[it](current[sindex],
961  w_current[sindex],
962  curr_walker[it],false,
963  data_arr[sindex+n_walk*n_threads]);
964  }
965  }
966 
967  }
968 
969  }
970  }
971  // End of parallel region
972 
973  // Verbose output
974  if (verbose>=2) {
975  for(size_t it=0;it<n_threads;it++) {
976  size_t sindex=n_walk*it+curr_walker[it];
977  scr_out.precision(4);
978  scr_out << "mcmc (" << it << "): iter: ";
979  scr_out.width((int)(1.0+log10((double)(nparams-1))));
980  scr_out << mcmc_iters << " i_walk: "
981  << curr_walker[it] << " log_wgt: "
982  << w_current[sindex] << std::endl;
983  scr_out.precision(6);
984  }
985  }
986 
987  // Collect best point
988  for(size_t it=0;it<n_threads;it++) {
989  if (func_ret[it]==o2scl::success && w_best>w_next[it]) {
990  best=next[it];
991  w_best=w_next[it];
992  if (switch_arr[n_walk*it+curr_walker[it]]==false) {
993  best_point(best,w_best,data_arr[curr_walker[it]+n_walk*it+
994  n_threads*n_walk]);
995  } else {
996  best_point(best,w_best,data_arr[curr_walker[it]+n_walk*it]);
997  }
998  }
999  }
1000 
1001  // Check to see if mcmc_done was returned or if meas_ret
1002  // returned an error
1003  for(size_t it=0;it<n_threads;it++) {
1004  if (meas_ret[it]==mcmc_done || func_ret[it]==mcmc_done) {
1005  main_done=true;
1006  }
1007  if (meas_ret[it]!=mcmc_done && meas_ret[it]!=o2scl::success) {
1008  if (err_nonconv) {
1009  O2SCL_ERR((((std::string)"Measurement function returned ")+
1010  o2scl::dtos(meas_ret[it])+
1011  " in mcmc_para_base::mcmc().").c_str(),
1013  }
1014  main_done=true;
1015  }
1016  }
1017 
1018  if (main_done==false) {
1019 
1020  mcmc_iters++;
1021 
1022  if (warm_up && mcmc_iters==n_warm_up) {
1023  warm_up=false;
1024  mcmc_iters=0;
1025  for(size_t it=0;it<n_threads;it++) {
1026  n_accept[it]=0;
1027  n_reject[it]=0;
1028  }
1029  if (verbose>=1) {
1030  scr_out << "Finished warmup." << std::endl;
1031  }
1032 
1033  }
1034  }
1035 
1036  if (verbose>=3) {
1037  std::cout << "Press a key and type enter to continue. ";
1038  char ch;
1039  std::cin >> ch;
1040  }
1041 
1042  if (main_done==false && warm_up==false && max_iters>0 &&
1043  mcmc_iters==max_iters) {
1044  scr_out << "mcmc (0): Stopping because number of iterations "
1045  << "equal to 'max_iters'." << std::endl;
1046  main_done=true;
1047  }
1048 
1049  if (main_done==false) {
1050  // Check to see if we're out of time
1051 #ifdef O2SCL_MPI
1052  double elapsed=MPI_Wtime()-mpi_start_time;
1053 #else
1054  double elapsed=time(0)-mpi_start_time;
1055 #endif
1056  if (max_time>0.0 && elapsed>max_time) {
1057  if (verbose>=0) {
1058  scr_out << "mcmc (0): Stopping because elapsed > max_time."
1059  << std::endl;
1060  }
1061  main_done=true;
1062  }
1063  }
1064 
1065  // --------------------------------------------------------------
1066  // End of main loop
1067  }
1068 
1069  // --------------------------------------------------------------
1070 
1071  mcmc_cleanup();
1072 
1073  return 0;
1074  }
1075  //@}
1076 
1077  /// \name Proposal distribution
1078  //@{
1079  /** \brief Set the proposal distribution
1080  */
1082  prop_dist=&p;
1083  pd_mode=true;
1084  aff_inv=false;
1085  n_walk=1;
1086  return;
1087  }
1088 
1089  /** \brief Go back to random-walk Metropolis with a uniform distribution
1090  */
1091  virtual void unset_proposal() {
1092  if (pd_mode) {
1093  prop_dist=0;
1094  pd_mode=false;
1095  }
1096  aff_inv=false;
1097  n_walk=1;
1098  return;
1099  }
1100  //@}
1101 
1102  };
1103 
1104  /** \brief A generic MCMC simulation class writing data to a
1105  \ref o2scl::table_units object
1106 
1107  This class performs a MCMC simulation and stores the
1108  results in a \ref o2scl::table_units object. The
1109  user must specify the column names and units in
1110  \ref set_names_units() before \ref mcmc() is called.
1111 
1112  The function \ref add_line is the measurement function of type
1113  \c measure_t in the parent. The overloaded function \ref mcmc()
1114  in this class works a bit differently in that it takes a
1115  function object (type \c fill_t) of the form
1116  \code
1117  int fill_func(const vec_t &pars, double log_weight,
1118  std::vector<double> &line, data_t &dat);
1119  \endcode
1120  which should store any auxillary values stored in the data
1121  object to \c line, in order to be added to the table.
1122 
1123  The output table will contain the parameters, the logarithm of
1124  the function (called "log_wgt") and a multiplying factor called
1125  "mult". This "fill" function is called only when a step is
1126  accepted and the multiplier for that row is set to 1. If a
1127  future step is rejected, then the multiplier is increased by
1128  one, rather than adding the same row to the table again.
1129 
1130  This class forms the basis of the MCMC used in the Bayesian
1131  analysis of neutron star mass and radius in
1132  http://github.com/awsteiner/bamr .
1133 
1134  \note This class is experimental.
1135  */
1136  template<class func_t, class fill_t, class data_t, class vec_t=ubvector>
1137  class mcmc_para_table : public mcmc_para_base<func_t,
1138  std::function<int(const vec_t &,double,size_t,bool,data_t &)>,
1139  data_t,vec_t> {
1140 
1141  protected:
1142 
1143  /// Measurement functor type for the parent
1144  typedef std::function<int(const vec_t &,double,size_t,bool,data_t &)>
1146 
1147  /// Type of parent class
1149 
1150  /// Column names
1151  std::vector<std::string> col_names;
1152 
1153  /// Column units
1154  std::vector<std::string> col_units;
1155 
1156  /// Main data table for Markov chain
1157  std::shared_ptr<o2scl::table_units<> > table;
1158 
1159  bool first_write;
1160 
1161  /** \brief MCMC initialization function
1162 
1163  This function sets the column names and units.
1164  */
1165  virtual int mcmc_init() {
1166 
1167  // -----------------------------------------------------------
1168  // Init tables
1169 
1170  table=std::shared_ptr<o2scl::table_units<> >(new o2scl::table_units<>);
1171  table->new_column("thread");
1172  table->new_column("walker");
1173  table->new_column("mult");
1174  table->new_column("log_wgt");
1175  for(size_t i=0;i<col_names.size();i++) {
1176  table->new_column(col_names[i]);
1177  if (col_units[i].length()>0) {
1178  table->set_unit(col_names[i],col_units[i]);
1179  }
1180  }
1181 
1182  walker_rows.resize(this->n_walk*this->n_threads);
1183  for(size_t i=0;i<this->n_walk*this->n_threads;i++) {
1184  walker_rows[i]=-1;
1185  }
1186 
1187  /*
1188  if (this->verbose>=2) {
1189  std::cout << "mcmc: Table column names and units: " << std::endl;
1190  for(size_t i=0;i<tab->get_ncolumns();i++) {
1191  std::cout << tab->get_column_name(i) << " "
1192  << tab->get_unit(tab->get_column_name(i)) << std::endl;
1193  }
1194  }
1195 
1196  if (this->verbose>=2) {
1197  std::cout << "End mcmc_para_table::mcmc_init()." << std::endl;
1198  }
1199  */
1200 
1201  return parent_t::mcmc_init();
1202  }
1203 
1204  /** \brief Fill \c line with data for insertion into the table
1205  */
1206  virtual int fill_line(const vec_t &pars, double log_weight,
1207  std::vector<double> &line, data_t &dat,
1208  size_t i_walker, fill_t &fill) {
1209 
1210 #ifdef O2SCL_OPENMP
1211  size_t i_thread=omp_get_thread_num();
1212 #else
1213  size_t i_thread=0;
1214 #endif
1215 
1216  // Thread
1217  line.push_back(i_thread);
1218  // Walker (set later)
1219  line.push_back(i_walker);
1220  // Initial multiplier
1221  line.push_back(1.0);
1222  line.push_back(log_weight);
1223  for(size_t i=0;i<pars.size();i++) {
1224  line.push_back(pars[i]);
1225  }
1226  int tempi=fill(pars,log_weight,line,dat);
1227  return tempi;
1228  }
1229 
1230  /** \brief Record the last row in the table which corresponds
1231  to each walker
1232  */
1233  std::vector<int> walker_rows;
1234 
1235  /// Likelihood estimator
1237 
1238  public:
1239 
1240  /** \brief Write MCMC tables to files
1241  */
1242  virtual void write_files() {
1243 
1244 #ifdef O2SCL_MPI
1245  // Ensure that multiple threads aren't writing to the
1246  // filesystem at the same time
1247  int tag=0, buffer=0;
1248  if (this->mpi_nprocs>1 && this->mpi_rank>0) {
1249  MPI_Recv(&buffer,1,MPI_INT,this->mpi_rank-1,tag,MPI_COMM_WORLD,
1250  MPI_STATUS_IGNORE);
1251  }
1252 #endif
1253 
1255  std::string fname=this->prefix+"_"+o2scl::itos(this->mpi_rank)+"_out";
1256  hf.open_or_create(fname);
1257 
1258  if (first_write==false) {
1259  hf.set_szt("max_iters",this->max_iters);
1260  hf.sets("prefix",this->prefix);
1261  hf.seti("aff_inv",this->aff_inv);
1262  hf.setd("step_fac",this->step_fac);
1263  hf.setd("ai_initial_step",this->ai_initial_step);
1264  hf.set_szt("n_warm_up",this->n_warm_up);
1265  hf.seti("user_seed",this->user_seed);
1266  hf.seti("mpi_rank",this->mpi_rank);
1267  hf.seti("mpi_nprocs",this->mpi_nprocs);
1268  hf.seti("verbose",this->verbose);
1269  hf.set_szt("max_bad_steps",this->max_bad_steps);
1270  hf.set_szt("n_walk",this->n_walk);
1271  hf.set_szt("n_threads",this->n_threads);
1272  hf.seti("always_accept",this->always_accept);
1273  hf.seti("allow_estimates",allow_estimates);
1274  first_write=true;
1275  }
1276 
1277  hf.set_szt_vec("n_accept",this->n_accept);
1278  hf.set_szt_vec("n_reject",this->n_reject);
1279  hf.set_szt_arr2d_copy("ret_value_counts",this->ret_value_counts.size(),
1280  this->ret_value_counts[0].size(),
1281  this->ret_value_counts);
1282  hf.setd_arr2d_copy("initial_points",this->initial_points.size(),
1283  this->initial_points[0].size(),
1284  this->initial_points);
1285 
1286  hdf_output(hf,*table,"markov_chain0");
1287 
1288  hf.close();
1289 
1290 #ifdef O2SCL_MPI
1291  if (this->mpi_nprocs>1 && this->mpi_rank>0) {
1292  MPI_Send(&buffer,1,MPI_INT,this->mpi_rank+1,tag,MPI_COMM_WORLD);
1293  }
1294 #endif
1295 
1296  return;
1297  }
1298 
1299  /// If true, allow estimates of the weight (default false)
1301 
1302  mcmc_para_table() {
1303  allow_estimates=false;
1304  first_write;
1305  }
1306 
1307  /// \name Basic usage
1308  //@{
1309  /** \brief Set the table names and units
1310  */
1311  virtual void set_names_units(std::vector<std::string> names,
1312  std::vector<std::string> units) {
1313  col_names=names;
1314  col_units=units;
1315  return;
1316  }
1317 
1318  /** \brief Perform an MCMC simulation
1319 
1320  Perform an MCMC simulation over \c nparams parameters starting
1321  at initial point \c init, limiting the parameters to be between
1322  \c low and \c high, using \c func as the objective function and
1323  calling the measurement function \c meas at each MC point.
1324  */
1325  virtual int mcmc(size_t nparams,
1326  vec_t &low, vec_t &high, std::vector<func_t> &func,
1327  std::vector<fill_t> &fill) {
1328 
1329  // Set number of threads (this is done in the child as well, but
1330  // we need this number to set up the vector of measure functions
1331  // below).
1332 #ifdef O2SCL_OPENMP
1333  omp_set_num_threads(this->n_threads);
1334  this->n_threads=omp_get_num_threads();
1335 #else
1336  this->n_threads=1;
1337 #endif
1338 
1339  // Setup the vector of measure functions
1340  std::vector<internal_measure_t> meas(this->n_threads);
1341  for(size_t it=0;it<this->n_threads;it++) {
1342  meas[it]=std::bind
1343  (std::mem_fn<int(const vec_t &,double,size_t,bool,
1344  data_t &, size_t, fill_t &)>
1345  (&mcmc_para_table::add_line),this,std::placeholders::_1,
1346  std::placeholders::_2,std::placeholders::_3,std::placeholders::_4,
1347  std::placeholders::_5,it,std::ref(fill[it]));
1348  }
1349 
1350  return parent_t::mcmc(nparams,low,high,func,meas);
1351  }
1352 
1353  /** \brief Get the output table
1354  */
1355  std::shared_ptr<o2scl::table_units<> > get_table() {
1356  return table;
1357  }
1358 
1359  /** \brief Set the output table
1360  */
1361  void set_table(std::shared_ptr<o2scl::table_units<> > &t) {
1362  table=t;
1363  return;
1364  }
1365 
1366  /** \brief Determine the chain sizes
1367  */
1368  void get_chain_sizes(std::vector<size_t> &csizes) {
1369 
1370  size_t ntot=this->n_threads*this->n_walk;
1371  csizes.resize(ntot);
1372 
1373  for(size_t it=0;it<this->n_threads;it++) {
1374  for(size_t iw=0;iw<this->n_walk;iw++) {
1375  size_t ix=it*this->n_walk+iw;
1376  size_t istart=ix;
1377  csizes[ix]=0;
1378  for(size_t j=istart;j<table->get_nlines();j+=ntot) {
1379  if (table->get("mult",j)>0.5) csizes[ix]++;
1380  }
1381  }
1382  }
1383 
1384  return;
1385  }
1386 
1387  /** \brief A measurement function which adds the point to the
1388  table
1389  */
1390  virtual int add_line(const vec_t &pars, double log_weight,
1391  size_t walker_ix, bool new_meas, data_t &dat,
1392  size_t i_thread, fill_t &fill) {
1393 
1394  // The combined walker/thread index
1395  size_t windex=i_thread*this->n_walk+walker_ix;
1396 
1397  // The total number of walkers * threads
1398  size_t ntot=this->n_threads*this->n_walk;
1399 
1400  int ret_value=o2scl::success;
1401 
1402  // Make sure only one thread writes to the table at a time by
1403  // making this next region 'critical'. This is required because
1404  // writes to the table may require resizing the table structure.
1405 #ifdef O2SCL_OPENMP
1406 #pragma omp critical (o2scl_mcmc_para_table_add_line)
1407 #endif
1408  {
1409 
1410  // If there's not enough space in the table for this iteration,
1411  // create it
1412  if ((walker_rows[windex]<0 && table->get_nlines()<ntot) ||
1413  table->get_nlines()<=walker_rows[windex]+ntot) {
1414  size_t istart=table->get_nlines();
1415  // Create enough space
1416  table->set_nlines(table->get_nlines()+ntot);
1417  // Now additionally initialize the first four colums
1418  for(size_t j=0;j<this->n_threads;j++) {
1419  for(size_t i=0;i<this->n_walk;i++) {
1420  table->set("thread",istart+j*this->n_walk+i,j);
1421  table->set("walker",istart+j*this->n_walk+i,i);
1422  table->set("mult",istart+j*this->n_walk+i,0.0);
1423  table->set("log_wgt",istart+j*this->n_walk+i,-1.0);
1424  }
1425  }
1426  }
1427 
1428  // Test to see if we need to add a new line of data or increment
1429  // the weight on the previous line. If the fill function has reset
1430  // the table data, then the walker_rows will refer to a row which
1431  // doesn't currently exist, so we have to add a new line of data.
1432 
1433  if (new_meas==true || walker_rows[windex]<0) {
1434 
1435  // We need to create a new measurement, so update the
1436  // row corresponding to this index
1437  if (walker_rows[windex]>=0) {
1438  walker_rows[windex]+=ntot;
1439  } else {
1440  walker_rows[windex]=windex;
1441  }
1442 
1443  if (walker_rows[windex]>=((int)(table->get_nlines()))) {
1444  O2SCL_ERR("Not enough space in table.",o2scl::exc_esanity);
1445  }
1446 
1447  std::vector<double> line;
1448  int fret=fill_line(pars,log_weight,line,dat,walker_ix,fill);
1449 
1450  // The fill_line() function doesn't set the walker index,
1451  // so we do this here
1452  //line[1]=walker_ix;
1453 
1454  if (fret!=o2scl::success) {
1455  // If we're done, we stop before adding the last point to the
1456  // table. This is important because otherwise the last line in
1457  // the table will always only have unit multiplicity, which
1458  // may or may not be correct.
1459  ret_value=this->mcmc_done;
1460  } else {
1461 
1462  if (line.size()!=table->get_ncolumns()) {
1463  std::cout << "line: " << line.size() << " columns: "
1464  << table->get_ncolumns() << std::endl;
1465  O2SCL_ERR("Table misalignment in mcmc_para_table::add_line().",
1466  exc_einval);
1467  }
1468 
1469  table->set_row(((size_t)walker_rows[windex]),line);
1470 
1471  }
1472 
1473  } else {
1474 
1475  // Otherwise, just increment the multiplier on the previous line
1476 
1477  double mult_old=table->get("mult",walker_rows[windex]);
1478  table->set("mult",walker_rows[windex],mult_old+1.0);
1479 
1480  }
1481 
1482  }
1483  // End of parallel region
1484 
1485  return ret_value;
1486  }
1487  //@}
1488 
1489  /** \brief Desc
1490  */
1491  virtual void mcmc_cleanup() {
1492 
1493  // This section removes empty rows at the end of the
1494  // table that were allocated but not used
1495  int i;
1496  bool done=false;
1497  for(i=table->get_nlines()-1;i>=0 && done==false;i--) {
1498  done=true;
1499  if (table->get("mult",i)<1.0e-10) {
1500  done=false;
1501  }
1502  }
1503  if (i+2<((int)table->get_nlines())) {
1504  table->set_nlines(i+2);
1505  }
1506  return parent_t::mcmc_cleanup();
1507  }
1508 
1509  /** \brief Desc
1510  */
1511  virtual void ac_coeffs(size_t ncols, ubmatrix &ac_coeffs) {
1512  std::vector<size_t> csizes;
1513  get_chain_sizes(csizes);
1514  size_t min_size=csizes[0];
1515  for(size_t i=1;i<csizes.size();i++) {
1516  if (csizes[i]<min_size) min_size=csizes[i];
1517  }
1518  size_t N_max=min_size/2;
1519  ac_coeffs.resize(ncols,N_max-1);
1520  for(size_t i=0;i<ncols;i++) {
1521  for(size_t ell=1;ell<N_max;ell++) {
1522  ac_coeffs(i,ell-1)=0.0;
1523  }
1524  }
1525  size_t n_tot=n_threads*n_walk;
1526  size_t table_row=0;
1527  size_t cstart=table->lookup_column("log_wgt")+1;
1528  for(size_t i=0;i<ncols;i++) {
1529  for(size_t j=0;j<n_threads;j++) {
1530  for(size_t k=0;k<n_walk;k++) {
1531  size_t tindex=j*n_walk+k;
1532  for(size_t ell=1;ell<N_max;ell++) {
1533  double mean=o2scl::vector_mean_double
1534  (csizes[tindex]+1,&((*table)[cstart+i][table_row]));
1535  ac_coeffs(i,ell-1)+=o2scl::vector_lagk_autocorr
1536  (csizes[tindex]+1,&((*table)[cstart+i][table_row]),
1537  ell,mean);
1538  }
1539  table_row+=csizes[tindex]+1;
1540  }
1541  }
1542  for(size_t ell=1;ell<N_max;ell++) {
1543  ac_coeffs(i,ell-1)/=((double)n_tot);
1544  }
1545  }
1546  return;
1547  }
1548 
1549  /** \brief Desc
1550  */
1551  virtual void ac_lengths(size_t ncols, ubmatrix &ac_coeffs_cols,
1552  ubvector &ac_lengths) {
1553  size_t N_max=ac_coeffs_cols.size2();
1554  ac_lengths.resize(ncols);
1555  for(size_t icol=0;icol<ncols;icol++) {
1556  std::vector<double> tau(N_max);
1557  for(size_t i=5;i<N_max;i++) {
1558  double sum=0.0;
1559  for(size_t j=0;j<i;j++) {
1560  sum+=ac_coeffs_cols(icol,j);
1561  }
1562  tau[i]=1.0+2.0*sum;
1563  std::cout << tau[i] << " " << ((double)i)/5.0 << std::endl;
1564  }
1565  std::cout << std::endl;
1566  }
1567  return;
1568  }
1569 
1570  /** \brief Reorder the table by thread and walker index
1571  */
1572  virtual void reorder_table() {
1573 
1574  // Create a new table
1575  std::shared_ptr<o2scl::table_units<> > table2=
1576  std::shared_ptr<o2scl::table_units<> >(new o2scl::table_units<>);
1577 
1578  for(size_t i=0;i<n_threads;i++) {
1579  for(size_t j=0;j<n_walk;j++) {
1580  std::string func=std::string("abs(walker-")+o2scl::szttos(j)+
1581  ")<0.1 && abs(thread-"+o2scl::szttos(i)+")<0.1";
1582  table->copy_rows(func,*table2);
1583  }
1584  }
1585 
1586  return;
1587  }
1588 
1589  /** \brief Reaverage the data into blocks of a fixed
1590  size in order to avoid autocorrelations
1591 
1592  \note The number of blocks \c n_blocks must be larger than the
1593  current table size. This function expects to find a column named
1594  "mult" which contains the multiplicity of each column, as is the
1595  case after a call to \ref mcmc_para_base::mcmc().
1596 
1597  This function is useful to remove autocorrelations to the table
1598  so long as the autocorrelation length is shorter than the block
1599  size. This function does not compute the autocorrelation length
1600  to check that this is the case.
1601  */
1602  void reblock(size_t n_blocks) {
1603 
1604  for(size_t it=0;it<this->n_threads;it++) {
1605 
1606  size_t n=table->get_nlines();
1607  if (n_blocks>n) {
1608  O2SCL_ERR2("Cannot reblock. Not enough data in ",
1609  "mcmc_para_table::reblock().",o2scl::exc_einval);
1610  }
1611  size_t n_block=n/n_blocks;
1612  size_t m=table->get_ncolumns();
1613  for(size_t j=0;j<n_blocks;j++) {
1614  double mult=0.0;
1615  ubvector dat(m);
1616  for(size_t i=0;i<m;i++) {
1617  dat[i]=0.0;
1618  }
1619  for(size_t k=j*n_block;k<(j+1)*n_block;k++) {
1620  mult+=(*table)["mult"][k];
1621  for(size_t i=1;i<m;i++) {
1622  dat[i]+=(*table)[i][k]*(*table)["mult"][k];
1623  }
1624  }
1625  table->set("mult",j,mult);
1626  for(size_t i=1;i<m;i++) {
1627  dat[i]/=mult;
1628  table->set(i,j,dat[i]);
1629  }
1630  }
1631  table->set_nlines(n_blocks);
1632 
1633  }
1634 
1635  return;
1636  }
1637 
1638  };
1639 
1640  // End of namespace
1641 }
1642 
1643 #endif
double vector_lagk_autocorr(size_t n, const vec_t &data, size_t k, double mean)
Lag-k autocorrelation.
Definition: vec_stats.h:705
std::shared_ptr< o2scl::table_units<> > table
Main data table for Markov chain.
Definition: mcmc_para.h:1157
std::vector< data_t > data_arr
Data array.
Definition: mcmc_para.h:157
std::string prefix
Prefix for output filenames.
Definition: mcmc_para.h:219
static const int mcmc_done
Integer to indicate completion.
Definition: mcmc_para.h:222
double ai_initial_step
Initial step fraction for affine-invariance sampling walkers (default 0.1)
Definition: mcmc_para.h:286
Multi-dimensional interpolation by inverse distance weighting.
Definition: interpm_idw.h:92
The main O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
virtual double log_pdf(const vec_t &x, const vec_t &x2) const =0
The log of the normalized density.
int set_szt_vec(std::string name, const std::vector< size_t > &v)
Set vector dataset named name with v.
size_t n_warm_up
Number of warm up steps (successful steps not iterations) (default 0)
Definition: mcmc_para.h:254
A generic MCMC simulation class writing data to a o2scl::table_units object.
Definition: mcmc_para.h:1137
bool allow_estimates
If true, allow estimates of the weight (default false)
Definition: mcmc_para.h:1300
bool err_nonconv
If true, call the error handler if msolve() or msolve_de() does not converge (default true) ...
Definition: mcmc_para.h:277
int verbose
Output control (default 0)
Definition: mcmc_para.h:262
std::vector< vec_t > current
Current points in parameter space for each walker and each OpenMP thread.
Definition: mcmc_para.h:148
virtual int mcmc(size_t nparams, vec_t &low, vec_t &high, std::vector< func_t > &func, std::vector< fill_t > &fill)
Perform an MCMC simulation.
Definition: mcmc_para.h:1325
void close()
Close the file.
void set_szt(std::string name, size_t u)
Set an unsigned integer named name to value u.
void setd(std::string name, double d)
Set a double named name to value d.
bool aff_inv
If true, use affine-invariant Monte Carlo.
Definition: mcmc_para.h:243
size_t max_iters
If non-zero, the maximum number of MCMC iterations (default 0)
Definition: mcmc_para.h:215
sanity check failed - shouldn&#39;t happen
Definition: err_hnd.h:65
std::vector< bool > switch_arr
Data switch array for each walker and each OpenMP thread.
Definition: mcmc_para.h:164
std::vector< std::string > col_names
Column names.
Definition: mcmc_para.h:1151
invalid argument supplied by user
Definition: err_hnd.h:59
virtual void ac_lengths(size_t ncols, ubmatrix &ac_coeffs_cols, ubvector &ac_lengths)
Desc.
Definition: mcmc_para.h:1551
int mpi_nprocs
The MPI number of processors.
Definition: mcmc_para.h:117
std::function< int(const vec_t &, double, size_t, bool, data_t &)> internal_measure_t
Measurement functor type for the parent.
Definition: mcmc_para.h:1145
std::vector< std::vector< size_t > > ret_value_counts
Return value counters, one vector for each OpenMP thread.
Definition: mcmc_para.h:168
void seti(std::string name, int i)
Set an integer named name to value i.
std::vector< size_t > curr_walker
Index of the current walker.
Definition: mcmc_para.h:210
double step_fac
Stepsize factor (default 10.0)
Definition: mcmc_para.h:246
A multi-dimensional conditional probability density function.
virtual void unset_proposal()
Go back to random-walk Metropolis with a uniform distribution.
Definition: mcmc_para.h:1091
virtual int add_line(const vec_t &pars, double log_weight, size_t walker_ix, bool new_meas, data_t &dat, size_t i_thread, fill_t &fill)
A measurement function which adds the point to the table.
Definition: mcmc_para.h:1390
bool always_accept
If true, accept all steps.
Definition: mcmc_para.h:281
std::vector< ubvector > initial_points
Initial points in parameter space.
Definition: mcmc_para.h:333
generic failure
Definition: err_hnd.h:61
virtual void write_files()
Write MCMC tables to files.
Definition: mcmc_para.h:1242
o2scl::prob_cond_mdim< vec_t > * prop_dist
Proposal distribution.
Definition: mcmc_para.h:134
size_t n_walk
Number of walkers for affine-invariant MC or 1 otherwise (default 1)
Definition: mcmc_para.h:272
virtual void mcmc_cleanup()
Desc.
Definition: mcmc_para.h:1491
void sets(std::string name, std::string s)
Set a string named name to value s.
virtual int mcmc(size_t nparams, vec_t &low, vec_t &high, std::vector< func_t > &func, std::vector< measure_t > &meas)
Perform an MCMC simulation.
Definition: mcmc_para.h:344
std::vector< std::string > col_units
Column units.
Definition: mcmc_para.h:1154
double mpi_start_time
The MPI starting time.
Definition: mcmc_para.h:120
void open_or_create(std::string fname)
Open a file named fname or create if it doesn&#39;t already exist.
virtual int fill_line(const vec_t &pars, double log_weight, std::vector< double > &line, data_t &dat, size_t i_walker, fill_t &fill)
Fill line with data for insertion into the table.
Definition: mcmc_para.h:1206
double max_time
Time in seconds (default is 0)
Definition: mcmc_para.h:125
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
std::string dtos(double x, int prec=6, bool auto_prec=false)
Convert a double to a string.
int user_seed
If non-zero, use as the seed for the random number generator (default 0)
Definition: mcmc_para.h:259
void reblock(size_t n_blocks)
Reaverage the data into blocks of a fixed size in order to avoid autocorrelations.
Definition: mcmc_para.h:1602
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:663
std::vector< size_t > n_reject
The number of Metropolis steps which were rejected in each thread (summed over all walkers) ...
Definition: mcmc_para.h:237
virtual void set_names_units(std::vector< std::string > names, std::vector< std::string > units)
Set the table names and units.
Definition: mcmc_para.h:1311
bool warm_up
If true, we are in the warm up phase.
Definition: mcmc_para.h:140
std::vector< size_t > n_accept
The number of Metropolis steps which were accepted in each thread (summed over all walkers) ...
Definition: mcmc_para.h:232
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
bool pd_mode
If true, then use the user-specified proposal distribution.
Definition: mcmc_para.h:137
virtual void reorder_table()
Reorder the table by thread and walker index.
Definition: mcmc_para.h:1572
std::vector< int > walker_rows
Record the last row in the table which corresponds to each walker.
Definition: mcmc_para.h:1233
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:459
virtual void mcmc_cleanup()
Cleanup after the MCMC.
Definition: mcmc_para.h:186
virtual int mcmc_init()
MCMC initialization function.
Definition: mcmc_para.h:1165
virtual void set_proposal(o2scl::prob_cond_mdim< vec_t > &p)
Set the proposal distribution.
Definition: mcmc_para.h:1081
Data table table class with units.
Definition: table_units.h:37
virtual int mcmc_init()
Initializations before the MCMC.
Definition: mcmc_para.h:174
int mpi_rank
The MPI processor rank.
Definition: mcmc_para.h:114
static const int mcmc_skip
Integer to indicate rejection.
Definition: mcmc_para.h:225
size_t max_bad_steps
Maximum number of failed steps when generating initial points with affine-invariant sampling (default...
Definition: mcmc_para.h:267
Store data in an O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$sc...
Definition: hdf_file.h:96
virtual void best_point(vec_t &best, double w_best, data_t &dat)
Function to run for the best point.
Definition: mcmc_para.h:199
interpm_idw< double * > esti
Likelihood estimator.
Definition: mcmc_para.h:1236
size_t n_threads
Number of OpenMP threads.
Definition: mcmc_para.h:326
mcmc_para_base< func_t, internal_measure_t, data_t, vec_t > parent_t
Type of parent class.
Definition: mcmc_para.h:1148
std::string itos(int x)
Convert an integer to a string.
A generic MCMC simulation class.
Definition: mcmc_para.h:107
std::vector< rng_gsl > rg
Random number generators.
Definition: mcmc_para.h:131
void set_table(std::shared_ptr< o2scl::table_units<> > &t)
Set the output table.
Definition: mcmc_para.h:1361
std::string szttos(size_t x)
Convert a size_t to a string.
std::shared_ptr< o2scl::table_units<> > get_table()
Get the output table.
Definition: mcmc_para.h:1355
std::ofstream scr_out
The screen output file.
Definition: mcmc_para.h:128
virtual void ac_coeffs(size_t ncols, ubmatrix &ac_coeffs)
Desc.
Definition: mcmc_para.h:1511
Success.
Definition: err_hnd.h:47
void get_chain_sizes(std::vector< size_t > &csizes)
Determine the chain sizes.
Definition: mcmc_para.h:1368

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