mcmc.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.h
24  \brief File for definition of \ref o2scl::mcmc_base and
25  \ref o2scl::mcmc_table
26 */
27 #ifndef O2SCL_MCMC_H
28 #define O2SCL_MCMC_H
29 
30 #include <iostream>
31 #include <random>
32 
33 #include <boost/numeric/ublas/vector.hpp>
34 
35 #include <o2scl/uniform_grid.h>
36 #include <o2scl/table3d.h>
37 #include <o2scl/hdf_file.h>
38 #include <o2scl/exception.h>
39 #include <o2scl/prob_dens_func.h>
40 #include <o2scl/cholesky.h>
41 #include <o2scl/vector.h>
42 #include <o2scl/multi_funct.h>
43 
44 /** \brief Main namespace
45 
46  This file is documented in mcmc.h .
47 */
48 namespace o2scl {
49 
51 
52  /** \brief A generic MCMC simulation class
53 
54  This class generates a Markov chain of a user-specified
55  function. The chain can be generated using the
56  Metropolis-Hastings algorithm with a user-specified proposal
57  distribution or using the affine-invariant sampling method of
58  Goodman and Weare.
59 
60  By default, the Metropolis-Hastings algorithm is executed with a
61  simple walk, with steps in each dimension of size \f$
62  (\mathrm{high} - \mathrm{low})/\mathrm{step\_fac} \f$ with the
63  denominator specified in \ref step_fac.
64 
65  The function type is a template type, \c func_t, which should
66  be of the form
67  \code
68  int f(size_t num_of_parameters, const vec_t &parameters,
69  double &log_pdf, data_t &dat)
70  \endcode
71  which computes \c log_pdf, the natural logarithm of the function
72  value, for any point in parameter space (any point between \c
73  low and \c high ).
74 
75  If the function being simulated returns \ref mcmc_skip then the
76  point is automatically rejected. After each acceptance or
77  rejection, a user-specified "measurement" function (of type \c
78  measure_t ) is called, which can be used to store the results.
79  In order to stop the simulation, either this function or the
80  probability distribution being simulated should return the value
81  \ref mcmc_done .
82 
83  A generic proposal distribution can be specified in \ref
84  set_proposal(). To go back to the default random walk method,
85  one can call the function \ref unset_proposal().
86 
87  If \ref aff_inv is set to true and the number of walkers, \ref
88  n_walk is set to a number larger than 1, then affine-invariant
89  sampling is used. For affine-invariant sampling, the variable
90  \ref step_fac represents the value of \f$ a \f$, the
91  limits of the distribution for \f$ z \f$.
92 
93  In order to store data at each point, the user can store this
94  data in any object of type \c data_t . If affine-invariant
95  sampling is used, then each chain has it's own data object. The
96  class keeps twice as many copies of these data object as would
97  otherwise be required, in order to avoid copying of data objects
98  in the case that the steps are accepted or rejected.
99 
100  \note This class is experimental.
101  */
102  template<class func_t, class measure_t,
103  class data_t, class vec_t=ubvector> class mcmc_base {
104 
105  protected:
106 
107  /// Random number generator
109 
110  /// Proposal distribution
112 
113  /// If true, then use the user-specified proposal distribution
114  bool pd_mode;
115 
116  /// If true, we are in the warm up phase
117  bool warm_up;
118 
119  /// Current points in parameter space
120  std::vector<vec_t> current;
121 
122  /// Data array
123  std::vector<data_t> data_arr;
124 
125  /// Data switch array
126  std::vector<bool> switch_arr;
127 
128  /// Return value counters
129  std::vector<size_t> ret_value_counts;
130 
131  /// \name Interface customization
132  //@{
133  /** \brief Initializations before the MCMC
134  */
135  virtual int mcmc_init() {
136  return 0;
137  }
138 
139  /** \brief Cleanup after the MCMC
140  */
141  virtual void mcmc_cleanup() {
142  return;
143  }
144 
145  /** \brief Function to run for the best point
146  */
147  virtual void best_point(vec_t &best, double w_best, data_t &dat) {
148  return;
149  }
150  //@}
151 
152  /// Index of the current walker
153  size_t curr_walker;
154 
155  /// Number of initial points specified by the user;
157 
158  public:
159 
160  /// Integer to indicate completion
161  static const int mcmc_done=-10;
162 
163  /// Integer to indicate rejection
164  static const int mcmc_skip=-20;
165 
166  /// \name Output quantities
167  //@{
168  /// The number of Metropolis steps which were accepted
169  size_t n_accept;
170 
171  /// The number of Metropolis steps which were rejected
172  size_t n_reject;
173  //@}
174 
175  /// \name Settings
176  //@{
177  /// If true, use affine-invariant Monte Carlo
178  bool aff_inv;
179 
180  /// Stepsize factor (default 10.0)
181  double step_fac;
182 
183  /** \brief Number of warm up steps (successful steps not
184  iterations) (default 0)
185 
186  \note Not to be confused with <tt>warm_up</tt>, which is
187  a boolean local variable in some functions not an int.
188  */
189  size_t n_warm_up;
190 
191  /** \brief If non-zero, use as the seed for the random number
192  generator (default 0)
193  */
195 
196  /// Output control (default 0)
197  int verbose;
198 
199  /** \brief Maximum number of failed steps when generating initial points
200  with affine-invariant sampling (default 1000)
201  */
203 
204  /** \brief Number of walkers for affine-invariant MC or 1
205  otherwise (default 1)
206  */
207  size_t n_walk;
208 
209  /** \brief If true, call the error handler if msolve() or
210  msolve_de() does not converge (default true)
211  */
213 
214  /** \brief If true, accept all steps
215  */
217 
218  /** \brief Initial step fraction for affine-invariance sampling walkers
219  */
221  //@}
222 
223  mcmc_base() {
224  user_seed=0;
225  n_warm_up=0;
226 
227  // MC step parameters
228  aff_inv=false;
229  pd_mode=false;
230  step_fac=10.0;
231  n_walk=1;
232  err_nonconv=true;
233  verbose=0;
234  warm_up=false;
235  max_bad_steps=1000;
236 
237  n_accept=0;
238  n_reject=0;
239  prop_dist=0;
240  always_accept=false;
241  ai_initial_step=0.1;
242 
243  n_init_points=0;
244  }
245 
246  /** \brief Default method for setting the random seed
247  */
248  virtual void set_seed() {
249  // Set RNG seed
250  unsigned long int seed=time(0);
251  if (user_seed!=0) {
252  seed=user_seed;
253  }
254  rg.set_seed(seed);
255 
256  return;
257  }
258 
259  /// \name Basic usage
260  //@{
261  /** \brief Perform an MCMC simulation
262 
263  Perform an MCMC simulation over \c nparams parameters starting
264  at initial point \c init, limiting the parameters to be between
265  \c low and \c high, using \c func as the objective function and
266  calling the measurement function \c meas at each MC point.
267  */
268  virtual int mcmc(size_t nparams, vec_t &init,
269  vec_t &low, vec_t &high, func_t &func,
270  measure_t &meas) {
271 
272  bool valid_guess=true;
273  for(size_t k=0;k<nparams;k++) {
274  if (init[k]<low[k] || init[k]>high[k]) valid_guess=false;
275  }
276  if (!valid_guess) {
277  O2SCL_ERR2("Initial guess outside of user-specified ",
278  "lower and upper bounds in mcmc_base::mcmc().",
280  }
281 
282  // Fix the number of walkers if it is too small
283  if (aff_inv && n_walk<=1) n_walk=2;
284 
285  // Fix 'step_fac' if it's less than or equal to zero
286  if (step_fac<=0.0) {
287  step_fac=10.0;
288  }
289 
290  set_seed();
291 
292  // Keep track of successful and failed MH moves
293  n_accept=0;
294  n_reject=0;
295 
296  // Warm-up flag, not to be confused with 'n_warm_up', i.e. the
297  // number of warm_up iterations.
298  warm_up=true;
299  if (n_warm_up==0) warm_up=false;
300 
301  // Allocate current point and current weight
302  current.resize(n_walk);
303  std::vector<double> w_current(n_walk);
304  w_current[0]=0.0;
305  for(size_t i=0;i<n_walk;i++) {
306  current[i].resize(nparams);
307  w_current[i]=0.0;
308  }
309 
310  // Initialize data and switch arrays
311  data_arr.resize(2*n_walk);
312  switch_arr.resize(n_walk);
313  for(size_t i=0;i<switch_arr.size();i++) switch_arr[i]=false;
314 
315  // Next and best points in parameter space
316  vec_t next(nparams), best(nparams);
317  double w_next=0.0, w_best=0.0, q_prop=0.0;
318 
319  // Set current to initial point
320  for(size_t i=0;i<nparams;i++) {
321  current[0][i]=init[i];
322  }
323 
324  // Return value of the measurement function
325  int meas_ret=0;
326 
327  // Run mcmc_init() function. The initial point, stored in
328  // current[0] can be modified by this function. Note that the
329  // local variable 'init' is not accessible to the mcmc_init()
330  // function.
331  int iret=mcmc_init();
332  if (iret!=0) {
333  O2SCL_ERR("Function mcmc_init() failed in mcmc_base::mcmc().",
335  return iret;
336  }
337 
338  if (verbose>=1) {
339  if (aff_inv) {
340  std::cout << "mcmc: Affine-invariant step, n_parameters="
341  << nparams << " n_walk=" << n_walk << std::endl;
342  } else if (pd_mode==true) {
343  std::cout << "mcmc: With proposal distribution, n_parameters="
344  << nparams << std::endl;
345  } else {
346  std::cout << "mcmc: Random-walk with uniform dist., n_parameters="
347  << nparams << std::endl;
348  }
349  }
350 
351  // ---------------------------------------------------
352  // Compute initial point and initial weights
353 
354  // Stretch-move algorithm
355  if (aff_inv) {
356 
357  // The mcmc_init() function may have changed the
358  // initial point, so we copy back to init here for
359  // temporary storage
360  for(size_t ipar=0;ipar<nparams;ipar++) {
361  init[ipar]=current[0][ipar];
362  }
363 
364  // Initialize each walker in turn
365  for(curr_walker=0;curr_walker<n_walk;curr_walker++) {
366 
367  size_t init_iters=0;
368  bool done=false;
369 
370  // If we already have a guess, use that
371  if (curr_walker<n_init_points) {
372 
373  // Compute the weight
374  int iret=func(nparams,current[curr_walker],
375  w_current[curr_walker],data_arr[curr_walker]);
376 
377  if (verbose>=1) {
378  std::cout.precision(4);
379  std::cout << "mcmc: ";
380  std::cout.width((int)(1.0+log10((double)(n_walk-1))));
381  std::cout << curr_walker << " " << w_current[curr_walker]
382  << " " << iret << " (initial; ai)" << std::endl;
383  std::cout.precision(6);
384  }
385 
386  // ------------------------------------------------
387  // If we have a good point, call the measurement function
388 
389  if (iret==o2scl::success) {
390  if (iret>=0 && iret<((int)ret_value_counts.size())) {
391  ret_value_counts[iret]++;
392  }
393  meas_ret=meas(current[curr_walker],w_current[curr_walker],
394  curr_walker,true,data_arr[curr_walker]);
395  done=true;
396  }
397  }
398 
399  while (!done) {
400 
401  // Make a perturbation from the initial point
402  for(size_t ipar=0;ipar<nparams;ipar++) {
403  if (init[ipar]<low[ipar] || init[ipar]>high[ipar]) {
404  O2SCL_ERR((((std::string)"Parameter ")+o2scl::szttos(ipar)+
405  " of "+o2scl::szttos(nparams)+
406  " out of range (value="+o2scl::dtos(init[ipar])+
407  ") in mcmc_base::mcmc().").c_str(),
409  }
410  do {
411  current[curr_walker][ipar]=init[ipar]+(rg.random()*2.0-1.0)*
412  (high[ipar]-low[ipar])*ai_initial_step;
413  } while (current[curr_walker][ipar]>high[ipar] ||
414  current[curr_walker][ipar]<low[ipar]);
415  }
416 
417  // Compute the weight
418  int iret=func(nparams,current[curr_walker],
419  w_current[curr_walker],data_arr[curr_walker]);
420 
421  if (verbose>=1) {
422  std::cout.precision(4);
423  std::cout << "mcmc: ";
424  std::cout.width((int)(1.0+log10((double)(n_walk-1))));
425  std::cout << curr_walker << " " << w_current[curr_walker]
426  << " " << iret << " (initial; ai)" << std::endl;
427  std::cout.precision(6);
428  }
429 
430  // ------------------------------------------------
431 
432  // Increment iteration count
433  init_iters++;
434 
435  // If we have a good point, call the measurement function and
436  // stop the loop
437  if (iret==o2scl::success) {
438  if (iret>=0 && iret<((int)ret_value_counts.size())) {
439  ret_value_counts[iret]++;
440  }
441  meas_ret=meas(current[curr_walker],w_current[curr_walker],
442  curr_walker,true,data_arr[curr_walker]);
443  done=true;
444  } else if (init_iters>max_bad_steps) {
445  if (err_nonconv) {
446  O2SCL_ERR2("Initial walkers failed in ",
447  "mcmc_base::mcmc().",o2scl::exc_einval);
448  }
449  return 1;
450  }
451 
452  if (verbose>=2) {
453  std::cout << "Press a key and type enter to continue. ";
454  char ch;
455  std::cin >> ch;
456  }
457 
458  }
459 
460  }
461 
462  } else {
463 
464  curr_walker=0;
465 
466  // Uniform random-walk steps
467 
468  // Compute weight for initial point
469  int iret=func(nparams,current[0],w_current[0],data_arr[0]);
470 
471  if (verbose>=1) {
472  std::cout.precision(4);
473  std::cout << "mcmc: " << w_current[0] << " (initial)" << std::endl;
474  std::cout.precision(6);
475  }
476 
477  if (iret!=o2scl::success) {
478  if (err_nonconv) {
479  O2SCL_ERR("Initial weight vanished in mcmc_base::mcmc().",
481  }
482  return 2;
483  }
484 
485  if (iret>=0 && iret<((int)ret_value_counts.size())) {
486  ret_value_counts[iret]++;
487  }
488  meas_ret=meas(current[0],w_current[0],0,true,data_arr[0]);
489 
490  best=current[0];
491  w_best=w_current[0];
492  best_point(best,w_best,data_arr[0]);
493 
494  if (meas_ret==mcmc_done) {
495  mcmc_cleanup();
496  return 0;
497  }
498 
499  }
500 
501  // ---------------------------------------------------
502 
503  bool main_done=false;
504  size_t mcmc_iters=0;
505 
506  while (!main_done) {
507 
508  // Walker to move (or zero when aff_inv is false)
509  curr_walker=0;
510  double smove_z=0.0;
511 
512  // ---------------------------------------------------
513  // Select next point
514 
515  if (aff_inv) {
516 
517  // Choose walker to move
518  curr_walker=mcmc_iters % n_walk;
519 
520  // Choose jth walker
521  size_t ij;
522  do {
523  ij=((size_t)(rg.random()*((double)n_walk)));
524  } while (ij==curr_walker || ij>=n_walk);
525 
526  // Select z
527  double p=rg.random();
528  double a=step_fac;
529  smove_z=(1.0-2.0*p+2.0*a*p+p*p-2.0*a*p*p+a*a*p*p)/a;
530 
531  // Create new trial point
532  for(size_t i=0;i<nparams;i++) {
533  next[i]=current[ij][i]+smove_z*(current[curr_walker][i]-
534  current[ij][i]);
535  }
536 
537  } else if (pd_mode) {
538 
539  // Use proposal distribution and compute associated weight
540  (*prop_dist)(current[0],next);
541  q_prop=prop_dist->log_pdf(current[0],next)-
542  prop_dist->log_pdf(next,current[0]);
543  if (!std::isfinite(q_prop)) {
544  O2SCL_ERR2("Proposal distribution not finite in ",
545  "mcmc_base::mcmc().",o2scl::exc_efailed);
546  }
547 
548  } else {
549 
550  // Uniform random-walk step
551  for(size_t k=0;k<nparams;k++) {
552  next[k]=current[0][k]+(rg.random()*2.0-1.0)*
553  (high[k]-low[k])/step_fac;
554  }
555 
556  }
557 
558  // ---------------------------------------------------
559  // Compute next weight
560 
561  int iret=o2scl::success;
562  // If the next point out of bounds, ensure that the point is rejected
563  for(size_t k=0;k<nparams;k++) {
564  if (next[k]<low[k] || next[k]>high[k]) iret=mcmc_skip;
565  }
566  if (iret!=mcmc_skip) {
567  if (switch_arr[curr_walker]==false) {
568  iret=func(nparams,next,w_next,data_arr[curr_walker+n_walk]);
569  } else {
570  iret=func(nparams,next,w_next,data_arr[curr_walker]);
571  }
572  if (iret>=0 && iret<((int)ret_value_counts.size())) {
573  ret_value_counts[iret]++;
574  }
575  if (verbose>=1 && iret!=o2scl::success && iret!=mcmc_skip) {
576  std::cout << "Function returned failure " << iret
577  << " at point ";
578  for(size_t k=0;k<nparams;k++) {
579  std::cout << next[k] << " ";
580  }
581  std::cout << std::endl;
582  }
583  } else if (verbose>=2) {
584  std::cout << "Parameter(s) out of range: " << std::endl;
585  std::cout.setf(std::ios::showpos);
586  for(size_t k=0;k<nparams;k++) {
587  std::cout << k << " " << low[k] << " " << next[k] << " "
588  << high[k];
589  if (next[k]<low[k] || next[k]>high[k]) {
590  std::cout << " <-";
591  }
592  std::cout << std::endl;
593  }
594  std::cout.unsetf(std::ios::showpos);
595  char ch;
596  std::cin >> ch;
597  }
598  if (iret==o2scl::success && w_best>w_next) {
599  best=next;
600  w_best=w_next;
601  if (switch_arr[curr_walker]==false) {
602  best_point(best,w_best,data_arr[curr_walker+n_walk]);
603  } else {
604  best_point(best,w_best,data_arr[curr_walker]);
605  }
606  }
607 
608  // ---------------------------------------------------
609 
610  bool accept=false;
611  if (always_accept && iret==success) accept=true;
612 
613  if (iret==o2scl::success) {
614  double r=rg.random();
615 
616  if (aff_inv) {
617  if (r<pow(smove_z,((double)nparams)-1.0)*
618  exp(w_next-w_current[curr_walker])) {
619  accept=true;
620  }
621  if (verbose>=1) {
622  std::cout.precision(4);
623  std::cout << "mcmc: ";
624  std::cout.width((int)(1.0+log10((double)(nparams-1))));
625  std::cout << curr_walker << " "
626  << w_current[curr_walker] << " " << w_next << " "
627  << smove_z << " ratio: "
628  << pow(smove_z,((double)nparams)-1.0)*
629  exp(w_next-w_current[curr_walker])
630  << " accept: " << accept << std::endl;
631  std::cout.precision(6);
632  }
633  } else if (pd_mode) {
634  if (r<exp(w_next-w_current[0]+q_prop)) {
635  accept=true;
636  }
637  if (verbose>=1) {
638  std::cout.precision(4);
639  std::cout << "mcmc: " << w_current[0]
640  << " " << w_next << " " << q_prop << " ratio: "
641  << exp(w_next-w_current[0]+q_prop)
642  << " accept: " << accept << std::endl;
643  std::cout.precision(6);
644  }
645  } else {
646  // Metropolis algorithm
647  if (r<exp(w_next-w_current[0])) {
648  accept=true;
649  }
650  if (verbose>=1) {
651  std::cout.precision(4);
652  std::cout << "mcmc: " << w_current[0] << " " << w_next
653  << " ratio: "
654  << exp(w_next-w_current[0])
655  << " accept: " << accept << std::endl;
656  std::cout.precision(6);
657  }
658  }
659 
660  // End of 'if (iret==o2scl::success)'
661  }
662 
663  if (accept) {
664 
665  n_accept++;
666 
667  // Store results from new point
668  if (!warm_up) {
669  if (switch_arr[curr_walker]==false) {
670  meas_ret=meas(next,w_next,curr_walker,true,
671  data_arr[curr_walker+n_walk]);
672  } else {
673  meas_ret=meas(next,w_next,curr_walker,true,
674  data_arr[curr_walker]);
675  }
676  }
677 
678  // Prepare for next point
679  current[curr_walker]=next;
680  w_current[curr_walker]=w_next;
681  switch_arr[curr_walker]=!(switch_arr[curr_walker]);
682 
683  } else {
684 
685  // Point was rejected
686  n_reject++;
687 
688  // Repeat measurement of old point
689  if (!warm_up) {
690  if (switch_arr[curr_walker]==false) {
691  meas_ret=meas(current[curr_walker],w_current[curr_walker],
692  curr_walker,false,data_arr[curr_walker]);
693  } else {
694  meas_ret=meas(current[curr_walker],w_current[curr_walker],
695  curr_walker,false,data_arr[curr_walker+n_walk]);
696  }
697  }
698 
699  }
700 
701  if (iret==mcmc_done) {
702  main_done=true;
703  } else {
704  if (meas_ret!=0) {
705  main_done=true;
706  if (meas_ret!=mcmc_done && err_nonconv) {
707  O2SCL_ERR((((std::string)"Measurement function returned ")+
708  o2scl::dtos(meas_ret)+
709  " in mcmc_base::mcmc().").c_str(),
711  }
712  }
713 
714  mcmc_iters++;
715 
716  if (warm_up && mcmc_iters==n_warm_up) {
717  warm_up=false;
718  mcmc_iters=0;
719  n_accept=0;
720  n_reject=0;
721  if (verbose>=1) {
722  std::cout << "Finished warmup." << std::endl;
723  }
724  }
725  }
726 
727  if (verbose>=2) {
728  std::cout << "Press a key and type enter to continue. ";
729  char ch;
730  std::cin >> ch;
731  }
732 
733  // --------------------------------------------------------------
734  // End of main loop
735  }
736 
737  // --------------------------------------------------------------
738 
739  mcmc_cleanup();
740 
741  return 0;
742  }
743  //@}
744 
745  /// \name Proposal distribution
746  //@{
747  /** \brief Set the proposal distribution
748  */
750  prop_dist=&p;
751  pd_mode=true;
752  aff_inv=false;
753  n_walk=1;
754  return;
755  }
756 
757  /** \brief Go back to random-walk Metropolis with a uniform distribution
758  */
759  virtual void unset_proposal() {
760  if (pd_mode) {
761  prop_dist=0;
762  pd_mode=false;
763  }
764  aff_inv=false;
765  n_walk=1;
766  return;
767  }
768  //@}
769 
770  };
771 
772  /** \brief A generic MCMC simulation class writing data to a
773  \ref o2scl::table_units object
774 
775  This class performs a MCMC simulation and stores the
776  results in a \ref o2scl::table_units object. The
777  user must specify the column names and units in
778  \ref set_names_units() before \ref mcmc() is called.
779 
780  The function \ref add_line is the measurement function of type
781  \c measure_t in the parent. The overloaded function \ref mcmc()
782  in this class works a bit differently in that it takes a
783  function object (type \c fill_t) of the form
784  \code
785  int fill_func(const vec_t &pars, double log_weight,
786  std::vector<double> &line, data_t &dat);
787  \endcode
788  which should store any auxillary values stored in the data
789  object to \c line, in order to be added to the table.
790 
791  The output table will contain the parameters, the logarithm of
792  the function (called "log_wgt") and a multiplying factor called
793  "mult". This "fill" function is called only when a step is
794  accepted and the multiplier for that row is set to 1. If a
795  future step is rejected, then the multiplier is increased by
796  one, rather than adding the same row to the table again.
797 
798  This class forms the basis of the MCMC used in the Bayesian
799  analysis of neutron star mass and radius in
800  http://github.com/awsteiner/bamr .
801 
802  \note This class is experimental.
803  */
804  template<class func_t, class fill_t, class data_t, class vec_t=ubvector>
805  class mcmc_table : public mcmc_base<func_t,
806  std::function<int(const vec_t &,double,size_t,bool,data_t &)>,
807  data_t,vec_t> {
808 
809  protected:
810 
811  /// Measurement functor type for the parent
812  typedef std::function<int(const vec_t &,double,size_t,bool,data_t &)>
814 
815  /// Type of parent class
817 
818  /// Column names
819  std::vector<std::string> col_names;
820 
821  /// Column units
822  std::vector<std::string> col_units;
823 
824  /// Main data table for Markov chain
825  std::shared_ptr<o2scl::table_units<> > tab;
826 
827  /** \brief MCMC initialization function
828 
829  This function sets the column names and units.
830  */
831  virtual int mcmc_init() {
832 
833  if (this->verbose>=2) {
834  std::cout << "Start mcmc_table::mcmc_init()." << std::endl;
835  }
836 
837  // -----------------------------------------------------------
838  // Init table
839 
840  std::string s, u;
841  tab->clear_table();
842  tab->new_column("mult");
843  tab->new_column("log_wgt");
844  for(size_t i=0;i<col_names.size();i++) {
845  tab->new_column(col_names[i]);
846  if (col_units[i].length()>0) {
847  tab->set_unit(col_names[i],col_units[i]);
848  }
849  }
850 
851  walker_rows.resize(this->n_walk);
852  for(size_t i=0;i<this->n_walk;i++) {
853  walker_rows[i]=-1;
854  }
855 
856  if (this->verbose>=2) {
857  std::cout << "mcmc: Table column names and units: " << std::endl;
858  for(size_t i=0;i<tab->get_ncolumns();i++) {
859  std::cout << tab->get_column_name(i) << " "
860  << tab->get_unit(tab->get_column_name(i)) << std::endl;
861  }
862  }
863 
864  if (this->verbose>=2) {
865  std::cout << "End mcmc_table::mcmc_init()." << std::endl;
866  }
867 
868  return 0;
869  }
870 
871  /** \brief Fill \c line with data for insertion into the table
872  */
873  virtual int fill_line(const vec_t &pars, double log_weight,
874  std::vector<double> &line, data_t &dat,
875  fill_t &fill) {
876 
877  // Initial multiplier
878  line.push_back(1.0);
879  line.push_back(log_weight);
880  for(size_t i=0;i<pars.size();i++) {
881  line.push_back(pars[i]);
882  }
883  return fill(pars,log_weight,line,dat);
884  }
885 
886  /** \brief Record the last row in the table which corresponds
887  to each walker
888  */
889  std::vector<int> walker_rows;
890 
891  public:
892 
893  mcmc_table() : tab(new o2scl::table_units<>) {
894  }
895 
896  /// \name Basic usage
897  //@{
898  /** \brief Set the table names and units
899  */
900  virtual void set_names_units(std::vector<std::string> names,
901  std::vector<std::string> units) {
902  col_names=names;
903  col_units=units;
904  return;
905  }
906 
907  /** \brief Perform an MCMC simulation
908 
909  Perform an MCMC simulation over \c nparams parameters starting
910  at initial point \c init, limiting the parameters to be between
911  \c low and \c high, using \c func as the objective function and
912  calling the measurement function \c meas at each MC point.
913  */
914  virtual int mcmc(size_t nparams, vec_t &init,
915  vec_t &low, vec_t &high, func_t &func,
916  fill_t &fill) {
917 
918  internal_measure_t meas=std::bind
919  (std::mem_fn<int(const vec_t &,double,size_t,bool,data_t &,fill_t &)>
920  (&mcmc_table::add_line),this,std::placeholders::_1,
921  std::placeholders::_2,std::placeholders::_3,std::placeholders::_4,
922  std::placeholders::_5,std::ref(fill));
923 
924  return parent_t::mcmc(nparams,init,low,high,func,meas);
925  }
926 
927  /** \brief Get the output table
928  */
929  std::shared_ptr<o2scl::table_units<> > get_table() {
930  return tab;
931  }
932 
933  /** \brief Set the output table
934  */
935  void set_table(std::shared_ptr<o2scl::table_units<> > &t) {
936  tab=t;
937  return;
938  }
939 
940  /** \brief A measurement function which adds the point to the
941  table
942  */
943  virtual int add_line(const vec_t &pars, double log_weight,
944  size_t walker_ix, bool new_meas, data_t &dat,
945  fill_t &fill) {
946 
947  // Test to see if we need to add a new line of data or increment
948  // the weight on the previous line. If the fill function has reset
949  // the table data, then the walker_rows will refer to a row which
950  // doesn't currently exist, so we have to add a new line of data.
951  if (new_meas==true ||
952  walker_rows[this->curr_walker]<0 ||
953  walker_rows[this->curr_walker]>=((int)(tab->get_nlines()))) {
954 
955  std::vector<double> line;
956  int fret=fill_line(pars,log_weight,line,dat,fill);
957 
958  if (fret!=o2scl::success) {
959  // If we're done, we stop before adding the last point to the
960  // table. This is important because otherwise the last line in
961  // the table will always only have unit multiplicity, which
962  // may or may not be correct.
963  if (fret==this->mcmc_done) {
964  if (this->verbose>=1) {
965  std::cout << "Fill function returned mcmc_done. "
966  << "Stopping run." << std::endl;
967  }
968  }
969  else {
970  if (this->verbose>=1) {
971  std::cout << "Fill function returned " << fret
972  << ". Stopping run." << std::endl;
973  }
974  }
975  return this->mcmc_done;
976  }
977 
978  if (line.size()!=tab->get_ncolumns()) {
979  std::cout << "line: " << line.size() << " columns: "
980  << tab->get_ncolumns() << std::endl;
981  O2SCL_ERR("Table misalignment in mcmc_table::add_line().",
982  exc_einval);
983  }
984 
985  if (this->verbose>=2) {
986  std::cout << "mcmc: Adding line:" << std::endl;
987  std::vector<std::string> sc_in, sc_out;
988  for(size_t k=0;k<line.size();k++) {
989  sc_in.push_back(tab->get_column_name(k)+": "+
990  o2scl::dtos(line[k]));
991  }
992  o2scl::screenify(line.size(),sc_in,sc_out);
993  for(size_t k=0;k<sc_out.size();k++) {
994  std::cout << sc_out[k] << std::endl;
995  }
996  }
997 
998  walker_rows[this->curr_walker]=tab->get_nlines();
999  tab->line_of_data(line.size(),line);
1000 
1001  } else {
1002 
1003  // Otherwise, just increment the multiplier on the previous line
1004  //std::cout << "nlines: " << tab->get_nlines() << std::endl;
1005  //std::cout << "walker: " << this->curr_walker << std::endl;
1006  //std::cout << "row: " << walker_rows[this->curr_walker]
1007  //<< std::endl;
1008  //O2SCL_ERR2("Sanity in row counting in ",
1009  //"mcmc_table::add_line().",o2scl::exc_esanity);
1010 
1011  double mult_old=tab->get("mult",walker_rows[this->curr_walker]);
1012  tab->set("mult",walker_rows[this->curr_walker],mult_old+1.0);
1013 
1014  if (this->verbose>=2) {
1015  std::cout << "mcmc: Updating line:" << std::endl;
1016  std::vector<std::string> sc_in, sc_out;
1017  for(size_t k=0;k<tab->get_ncolumns();k++) {
1018  sc_in.push_back
1019  (tab->get_column_name(k)+": "+
1020  o2scl::dtos(tab->get(tab->get_column_name(k),
1021  walker_rows[this->curr_walker])));
1022  }
1023  o2scl::screenify(tab->get_ncolumns(),sc_in,sc_out);
1024  for(size_t k=0;k<sc_out.size();k++) {
1025  std::cout << sc_out[k] << std::endl;
1026  }
1027  }
1028 
1029  }
1030 
1031  return 0;
1032  }
1033  //@}
1034 
1035  /** \brief Reaverage the data into blocks of a fixed
1036  size in order to avoid autocorrelations
1037 
1038  \note The number of blocks \c n_blocks must be larger than the
1039  current table size. This function expects to find a column named
1040  "mult" which contains the multiplicity of each column, as is the
1041  case after a call to \ref mcmc_base::mcmc().
1042 
1043  This function is useful to remove autocorrelations to the table
1044  so long as the autocorrelation length is shorter than the block
1045  size. This function does not compute the autocorrelation length
1046  to check that this is the case.
1047  */
1048  void reblock(size_t n_blocks) {
1049  size_t n=tab->get_nlines();
1050  if (n_blocks>n) {
1051  O2SCL_ERR2("Cannot reblock. Not enough data in ",
1052  "mcmc_table::reblock().",o2scl::exc_einval);
1053  }
1054  size_t n_block=n/n_blocks;
1055  size_t m=tab->get_ncolumns();
1056  for(size_t j=0;j<n_blocks;j++) {
1057  double mult=0.0;
1058  ubvector dat(m);
1059  for(size_t i=0;i<m;i++) {
1060  dat[i]=0.0;
1061  }
1062  for(size_t k=j*n_block;k<(j+1)*n_block;k++) {
1063  mult+=(*tab)["mult"][k];
1064  for(size_t i=1;i<m;i++) {
1065  dat[i]+=(*tab)[i][k]*(*tab)["mult"][k];
1066  }
1067  }
1068  tab->set("mult",j,mult);
1069  for(size_t i=1;i<m;i++) {
1070  dat[i]/=mult;
1071  tab->set(i,j,dat[i]);
1072  }
1073  }
1074  tab->set_nlines(n_blocks);
1075  return;
1076  }
1077 
1078  };
1079 
1080  // End of namespace
1081 }
1082 
1083 #endif
std::shared_ptr< o2scl::table_units<> > get_table()
Get the output table.
Definition: mcmc.h:929
std::shared_ptr< o2scl::table_units<> > tab
Main data table for Markov chain.
Definition: mcmc.h:825
virtual void set_seed()
Default method for setting the random seed.
Definition: mcmc.h:248
virtual int add_line(const vec_t &pars, double log_weight, size_t walker_ix, bool new_meas, data_t &dat, fill_t &fill)
A measurement function which adds the point to the table.
Definition: mcmc.h:943
std::vector< std::string > col_names
Column names.
Definition: mcmc.h:819
std::vector< std::string > col_units
Column units.
Definition: mcmc.h:822
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.
std::vector< vec_t > current
Current points in parameter space.
Definition: mcmc.h:120
o2scl::prob_cond_mdim< vec_t > * prop_dist
Proposal distribution.
Definition: mcmc.h:111
double step_fac
Stepsize factor (default 10.0)
Definition: mcmc.h:181
void set_table(std::shared_ptr< o2scl::table_units<> > &t)
Set the output table.
Definition: mcmc.h:935
std::vector< bool > switch_arr
Data switch array.
Definition: mcmc.h:126
size_t n_walk
Number of walkers for affine-invariant MC or 1 otherwise (default 1)
Definition: mcmc.h:207
virtual void unset_proposal()
Go back to random-walk Metropolis with a uniform distribution.
Definition: mcmc.h:759
void set_seed(unsigned long int s)
Set the seed.
Definition: rng_gsl.h:106
size_t n_reject
The number of Metropolis steps which were rejected.
Definition: mcmc.h:172
int verbose
Output control (default 0)
Definition: mcmc.h:197
virtual int mcmc(size_t nparams, vec_t &init, vec_t &low, vec_t &high, func_t &func, fill_t &fill)
Perform an MCMC simulation.
Definition: mcmc.h:914
virtual void set_proposal(o2scl::prob_cond_mdim< vec_t > &p)
Set the proposal distribution.
Definition: mcmc.h:749
A generic MCMC simulation class.
Definition: mcmc.h:103
invalid argument supplied by user
Definition: err_hnd.h:59
bool always_accept
If true, accept all steps.
Definition: mcmc.h:216
virtual void mcmc_cleanup()
Cleanup after the MCMC.
Definition: mcmc.h:141
double random()
Return a random number in .
Definition: rng_gsl.h:82
A multi-dimensional conditional probability density function.
size_t n_init_points
Number of initial points specified by the user;.
Definition: mcmc.h:156
generic failure
Definition: err_hnd.h:61
void reblock(size_t n_blocks)
Reaverage the data into blocks of a fixed size in order to avoid autocorrelations.
Definition: mcmc.h:1048
static const int mcmc_done
Integer to indicate completion.
Definition: mcmc.h:161
size_t n_accept
The number of Metropolis steps which were accepted.
Definition: mcmc.h:169
size_t n_warm_up
Number of warm up steps (successful steps not iterations) (default 0)
Definition: mcmc.h:189
A generic MCMC simulation class writing data to a o2scl::table_units object.
Definition: mcmc.h:805
bool warm_up
If true, we are in the warm up phase.
Definition: mcmc.h:117
std::vector< data_t > data_arr
Data array.
Definition: mcmc.h:123
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
virtual void set_names_units(std::vector< std::string > names, std::vector< std::string > units)
Set the table names and units.
Definition: mcmc.h:900
int user_seed
If non-zero, use as the seed for the random number generator (default 0)
Definition: mcmc.h:194
std::string dtos(double x, int prec=6, bool auto_prec=false)
Convert a double to a string.
bool aff_inv
If true, use affine-invariant Monte Carlo.
Definition: mcmc.h:178
size_t max_bad_steps
Maximum number of failed steps when generating initial points with affine-invariant sampling (default...
Definition: mcmc.h:202
size_t curr_walker
Index of the current walker.
Definition: mcmc.h:153
void screenify(size_t nin, const string_arr_t &in_cols, std::vector< std::string > &out_cols, size_t max_size=80)
Reformat the columns for output of width size.
Definition: misc.h:124
mcmc_base< func_t, internal_measure_t, data_t, vec_t > parent_t
Type of parent class.
Definition: mcmc.h:816
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
virtual void best_point(vec_t &best, double w_best, data_t &dat)
Function to run for the best point.
Definition: mcmc.h:147
std::vector< size_t > ret_value_counts
Return value counters.
Definition: mcmc.h:129
Data table table class with units.
Definition: table_units.h:37
static const int mcmc_skip
Integer to indicate rejection.
Definition: mcmc.h:164
Random number generator (GSL)
Definition: rng_gsl.h:55
virtual int mcmc_init()
Initializations before the MCMC.
Definition: mcmc.h:135
bool pd_mode
If true, then use the user-specified proposal distribution.
Definition: mcmc.h:114
bool err_nonconv
If true, call the error handler if msolve() or msolve_de() does not converge (default true) ...
Definition: mcmc.h:212
std::function< int(const vec_t &, double, size_t, bool, data_t &)> internal_measure_t
Measurement functor type for the parent.
Definition: mcmc.h:813
double ai_initial_step
Initial step fraction for affine-invariance sampling walkers.
Definition: mcmc.h:220
std::string szttos(size_t x)
Convert a size_t to a string.
Success.
Definition: err_hnd.h:47
std::vector< int > walker_rows
Record the last row in the table which corresponds to each walker.
Definition: mcmc.h:889
virtual int fill_line(const vec_t &pars, double log_weight, std::vector< double > &line, data_t &dat, fill_t &fill)
Fill line with data for insertion into the table.
Definition: mcmc.h:873
virtual int mcmc(size_t nparams, vec_t &init, vec_t &low, vec_t &high, func_t &func, measure_t &meas)
Perform an MCMC simulation.
Definition: mcmc.h:268
rng_gsl rg
Random number generator.
Definition: mcmc.h:108
virtual int mcmc_init()
MCMC initialization function.
Definition: mcmc.h:831

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