mcarlo_vegas.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-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 /* monte/vegas.c
24  *
25  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Michael Booth
26  *
27  * This program is free software; you can redistribute it and/or modify
28  * it under the terms of the GNU General Public License as published by
29  * the Free Software Foundation; either version 3 of the License, or (at
30  * your option) any later version.
31  *
32  * This program is distributed in the hope that it will be useful, but
33  * WITHOUT ANY WARRANTY; without even the implied warranty of
34  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
35  * General Public License for more details.
36  *
37  * You should have received a copy of the GNU General Public License
38  * along with this program; if not, write to the Free Software
39  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
40  * 02110-1301, USA.
41  */
42 #ifndef O2SCL_MCARLO_VEGAS_H
43 #define O2SCL_MCARLO_VEGAS_H
44 
45 /** \file mcarlo_vegas.h
46  \brief File defining \ref o2scl::mcarlo_vegas
47 */
48 
49 #include <iostream>
50 
51 #include <gsl/gsl_math.h>
52 #include <gsl/gsl_monte.h>
53 #include <gsl/gsl_machine.h>
54 #include <gsl/gsl_monte_vegas.h>
55 
56 #include <o2scl/mcarlo.h>
57 
58 #ifndef DOXYGEN_NO_O2NS
59 namespace o2scl {
60 #endif
61 
62  /** \brief Multidimensional integration using Vegas Monte Carlo (GSL)
63 
64  The output options are a little different than the original GSL
65  routine. The default setting of \ref mcarlo::verbose is 0,
66  which turns off all output. A verbose value of 1 prints summary
67  information about the weighted average and final result, while a
68  value of 2 also displays the grid coordinates. A value of 3
69  prints information from the rebinning procedure for each
70  iteration.
71 
72  Some original documentation from GSL:
73 
74  \verbatim
75  The input coordinates are x[j], with upper and lower limits
76  xu[j] and xl[j]. The integration length in the j-th direction is
77  delx[j]. Each coordinate x[j] is rescaled to a variable y[j] in
78  the range 0 to 1. The range is divided into bins with boundaries
79  xi[i][j], where i=0 corresponds to y=0 and i=bins to y=1. The
80  grid is refined (ie, bins are adjusted) using d[i][j] which is
81  some variation on the squared sum. A third parameter used in
82  defining the real coordinate using random numbers is called z.
83  It ranges from 0 to bins. Its integer part gives the lower index
84  of the bin into which a call is to be placed, and the remainder
85  gives the location inside the bin.
86 
87  When stratified sampling is used the bins are grouped into
88  boxes, and the algorithm allocates an equal number of function
89  calls to each box.
90 
91  The variable alpha controls how "stiff" the rebinning algorithm
92  is. alpha = 0 means never change the grid. Alpha is typically
93  set between 1 and 2.
94  \endverbatim
95 
96  \todo Mode = importance only doesn't give the same answer
97  as GSL yet.
98 
99  \future Prettify the verbose output
100 
101  \future Allow the user to get information about the how
102  the sampling was done, possibly by converting the bins
103  and boxes into a structure or class.
104 
105  \future Allow the user to change the maximum number of bins.
106 
107  Based on \ref Lepage78 . The current version of the algorithm
108  was described in the Cornell preprint CLNS-80/447 of March,
109  1980. The GSL code follows most closely the C version by D. R.
110  Yennie, coded in 1984.
111  */
112  template<class func_t=multi_funct11,
114  class rng_t=rng_gsl>
115  class mcarlo_vegas : public mcarlo<func_t,vec_t,rng_t> {
116 
117  public:
118 
122 
123  /// \name Integration mode (default is mode_importance)
124  //@{
125  int mode;
126  static const int mode_importance=1;
127  static const int mode_importance_only=0;
128  static const int mode_stratified=-1;
129  //@}
130 
131  /// Result from last iteration
132  double result;
133 
134  /// Uncertainty from last iteration
135  double sigma;
136 
137  /** \brief The stiffness of the rebinning algorithm (default 1.5)
138 
139  This usual range is between 1 and 2.
140  */
141  double alpha;
142 
143  /// Set the number of iterations (default 5)
144  unsigned int iterations;
145 
146  /** \brief The chi-squared per degree of freedom for the weighted
147  estimate of the integral
148 
149  After an integration, this should be close to 1. If it is not,
150  then this indicates that the values of the integral from
151  different iterations are inconsistent, and the error may be
152  underestimated. Further iterations of the algorithm may enable
153  one to obtain more reliable results.
154  */
155  double chisq;
156 
157  /// The output stream to send output information (default \c std::cout)
158  std::ostream *outs;
159 
160 #ifndef DOXYGEN_INTERNAL
161 
162  protected:
163 
164  /// Maximum number of bins
165  static const size_t bins_max=50;
166 
167  /// Number of dimensions
168  size_t dim;
169  /// Number of bins
170  unsigned int bins;
171  /// Number of boxes
172  unsigned int boxes;
173  /// Boundaries for each bin
174  ubvector xi;
175  /// Storage for grid refinement
176  ubvector xin;
177  /// The iteration length in each direction
178  ubvector delx;
179  /// The weight for each bin
180  ubvector weight;
181  /// The volume of the current bin
182  double vol;
183 
184  /// The bins for each direction
185  ubvector_int bin;
186  /// The boxes for each direction
187  ubvector_int box;
188 
189  /// Distribution
190  ubvector d;
191 
192  /** \name Scratch variables preserved between calls to
193  vegas_minteg_err()
194  */
195  //@{
196  double jac;
197  double wtd_int_sum;
198  double sum_wgts;
199  double chi_sum;
200  //@}
201 
202  /// The starting iteration number
203  unsigned int it_start;
204  /// The total number of iterations
205  unsigned int it_num;
206  /// Number of samples for computing chi squared
207  unsigned int samples;
208  /// Number of function calls per box
209  unsigned int calls_per_box;
210 
211  /// Initialize box coordinates
212  virtual void init_box_coord(ubvector_int &boxt) {
213  size_t i;
214  for (i=0;i<dim;i++) {
215  boxt[i]=0;
216  }
217  return;
218  }
219 
220  /** \brief Change box coordinates
221 
222  Steps through the box coordinates, e.g.
223  \verbatim
224  {0,0}, {0,1}, {0,2}, {0,3}, {1,0}, {1,1}, {1,2}, ...
225  \endverbatim
226 
227  This is among the member functions that is not virtual
228  because it is part of the innermost loop.
229  */
230  int change_box_coord(ubvector_int &boxt) {
231  int j=dim-1;
232  int ng=boxes;
233 
234  while (j >= 0) {
235 
236  boxt[j]=(boxt[j]+1) % ng;
237  if (boxt[j] != 0) {
238  return 1;
239  }
240  j--;
241 
242  }
243 
244  return 0;
245  }
246 
247  /// Initialize grid
248  virtual void init_grid(const vec_t &xl, const vec_t &xu, size_t ldim) {
249  size_t j;
250  vol=1.0;
251  bins=1;
252 
253  for (j=0;j<ldim;j++) {
254  double dx=xu[j]-xl[j];
255  delx[j]=dx;
256  vol *= dx;
257 
258  xi[j]=0.0;
259  xi[dim+j]=1.0;
260  }
261 
262  return;
263  }
264 
265  /// Reset grid values
266  virtual void reset_grid_values() {
267  size_t i, j;
268 
269  for (i=0;i<bins;i++) {
270  for (j=0;j<dim;j++) {
271  d[i*dim+j]=0.0;
272  }
273  }
274  return;
275  }
276 
277  /** \brief Add the most recently generated result to the distribution
278 
279  This is among the member functions that is not virtual
280  because it is part of the innermost loop.
281  */
282  void accumulate_distribution(ubvector_int &lbin, double y) {
283  size_t j;
284 
285  for (j=0;j<dim;j++) {
286  int i=lbin[j];
287  d[i*dim+j]+=y;
288  }
289  return;
290  }
291 
292  /** \brief Generate a random position in a given box
293 
294  Use the random number generator to return a random position x
295  in a given box. The value of bin gives the bin location of the
296  random position (there may be several bins within a given box)
297 
298  This is among the member functions that is not virtual
299  because it is part of the innermost loop.
300  */
301  void random_point(vec_t &lx, ubvector_int &lbin, double &bin_vol,
302  const ubvector_int &lbox, const vec_t &xl,
303  const vec_t &xu) {
304 
305  double lvol=1.0;
306 
307  size_t j;
308 
309  for (j=0;j<dim;++j) {
310 
311  // The equivalent of gsl_rng_uniform_pos()
312  double rdn;
313  do {
314  rdn=this->rng_dist(this->rng);
315  } while (rdn==0);
316 
317  /* lbox[j] + ran gives the position in the box units, while z
318  is the position in bin units. */
319  double z=((lbox[j]+rdn)/boxes)*bins;
320 
321  int k=(int)z;
322 
323  double y, bin_width;
324 
325  lbin[j]=k;
326 
327  if (k == 0) {
328  bin_width=xi[dim+j];
329  y=z*bin_width;
330  } else {
331  bin_width=xi[(k+1)*dim+j]-xi[k*dim+j];
332  y=xi[k*dim+j]+(z-k)*bin_width;
333  }
334 
335  lx[j]=xl[j]+y*delx[j];
336 
337  lvol *= bin_width;
338  }
339 
340  bin_vol=lvol;
341 
342  return;
343  }
344 
345  /// Resize the grid
346  virtual void resize_grid(unsigned int lbins) {
347  size_t j, k;
348 
349  /* weight is ratio of bin sizes */
350 
351  double pts_per_bin=(double) bins/(double) lbins;
352 
353  for (j=0;j<dim;j++) {
354  double xold;
355  double xnew=0;
356  double dw=0;
357  int i=1;
358 
359  for (k=1;k <= bins;k++) {
360  dw+=1.0;
361  xold=xnew;
362  xnew=xi[k*dim+j];
363 
364  for (;dw > pts_per_bin;i++) {
365  dw -= pts_per_bin;
366  (xin[i])=xnew-(xnew-xold)*dw;
367  }
368  }
369 
370  for (k=1 ;k<lbins;k++) {
371  xi[k*dim+j]=xin[k];
372  }
373 
374  xi[lbins*dim+j]=1;
375  }
376 
377  bins=lbins;
378 
379  return;
380  }
381 
382  /// Refine the grid
383  virtual void refine_grid() {
384 
385  size_t i, j, k;
386 
387  for (j=0;j<dim;j++) {
388  double grid_tot_j, tot_weight;
389 
390  double oldg=d[j];
391  double newg=d[dim+j];
392 
393  d[j]=(oldg+newg)/2;
394  grid_tot_j=d[j];
395 
396  /* This implements gs[i][j]=(gs[i-1][j]+gs[i][j]+gs[i+1][j])/3 */
397 
398  for (i=1;i<bins-1;i++) {
399  double rc=oldg+newg;
400  oldg=newg;
401  newg=d[(i+1)*dim+j];
402  d[i*dim+j]=(rc+newg)/3;
403  grid_tot_j+=d[i*dim+j];
404  }
405  d[(bins-1)*dim+j]=(newg+oldg)/2;
406 
407  grid_tot_j+=d[(bins-1)*dim+j];
408 
409  tot_weight=0;
410 
411  for (i=0;i<bins;i++) {
412  weight[i]=0;
413 
414  if (d[i*dim+j] > 0) {
415  oldg=grid_tot_j/d[i*dim+j];
416  /* damped change */
417  weight[i]=pow (((oldg-1)/oldg/log (oldg)), alpha);
418  }
419 
420  tot_weight+=weight[i];
421  }
422 
423  {
424  double pts_per_bin=tot_weight/bins;
425 
426  double xold;
427  double xnew=0;
428  double dw=0;
429  i=1;
430 
431  for (k=0;k<bins;k++) {
432  dw+=weight[k];
433  xold=xnew;
434  xnew=xi[(k+1)*dim+j];
435 
436  for (;dw > pts_per_bin;i++) {
437  dw -= pts_per_bin;
438  (xin[i])=xnew-(xnew-xold)*dw/weight[k];
439  }
440  }
441 
442  for (k=1 ;k<bins ;k++) {
443  xi[k*dim+j]=xin[k];
444  }
445 
446  xi[bins*dim+j]=1;
447  }
448  }
449  return;
450  }
451 
452  /// Print limits of integration
453  virtual void print_lim(const vec_t &xl, const vec_t &xu,
454  unsigned long ldim) {
455 
456  unsigned long j;
457 
458  (*outs) << "The limits of integration are:\n" << std::endl;
459  for (j=0;j<ldim;++j) {
460  (*outs) << "xl[" << j << "]=" << xl[j] << " xu[" << j
461  << "]=" << xu[j] << std::endl;
462  }
463  (*outs) << std::endl;
464 
465  return;
466  }
467 
468  /// Print header
469  virtual void print_head(unsigned long num_dim, unsigned long calls,
470  unsigned int lit_num, unsigned int lbins,
471  unsigned int lboxes) {
472 
473  (*outs) << "num_dim=" << num_dim << ", calls=" << calls
474  << ", it_num=" << lit_num << ", max_it_num="
475  << iterations << ", verb=" << this->verbose << ", alph=" << alpha
476  << ",\n mode=" << mode << ", boxes=" << lboxes
477  << "\n\n single.......iteration "
478  << "accumulated......results\n"
479  << "iter integral sigma integral "
480  << " sigma chi-sq/it" << std::endl;
481 
482  return;
483  }
484 
485  /// Print results
486  virtual void print_res(unsigned int itr, double res,
487  double err, double cum_res, double cum_err,
488  double chi_sq) {
489  (*outs).precision(5);
490  (*outs) << itr << " ";
491  outs->setf(std::ios::showpos);
492  (*outs) << res << " ";
493  outs->unsetf(std::ios::showpos);
494  (*outs) << err << " ";
495  outs->setf(std::ios::showpos);
496  (*outs) << cum_res << " ";
497  outs->unsetf(std::ios::showpos);
498  (*outs) << cum_err << " " << chi_sq << std::endl;
499  (*outs).precision(6);
500  return;
501  }
502 
503  /// Print distribution
504  virtual void print_dist(unsigned long ldim) {
505  unsigned long i, j;
506 
507  if (this->verbose<2) {
508  return;
509  }
510 
511  for (j=0;j<ldim;++j) {
512  (*outs) << "\n Axis " << j << std::endl;
513  (*outs) << " x g" << std::endl;
514  outs->setf(std::ios::showpos);
515  for (i=0;i<bins;i++) {
516  (*outs) << "weight [ " << (xi[(i)*dim+(j)]) << " , "
517  << xi[(i+1)*dim+j] << " ] = ";
518  (*outs) << " " << (d[(i)*dim+(j)]) << std::endl;
519  }
520  outs->unsetf(std::ios::showpos);
521  (*outs) << std::endl;
522  }
523  (*outs) << std::endl;
524  return;
525  }
526 
527  /// Print grid
528  virtual void print_grid(unsigned long ldim) {
529 
530  if (this->verbose<2) {
531  return;
532  }
533 
534  unsigned long i, j;
535  for (j=0;j<ldim;++j) {
536  (*outs) << "\n Axis " << j << std::endl;
537  (*outs) << " x " << std::endl;
538  outs->setf(std::ios::showpos);
539  for (i=0;i <= bins;i++) {
540  (*outs) << (xi[(i)*dim+(j)]) << " ";
541  if (i % 5 == 4) {
542  (*outs) << std::endl;
543  }
544  }
545  (*outs) << std::endl;
546  outs->unsetf(std::ios::showpos);
547  }
548  (*outs) << std::endl;
549  return;
550  }
551 
552  /// Point for function evaluation
553  vec_t x;
554 
555 #endif
556 
557  public:
558 
559  mcarlo_vegas() {
560  this->verbose=0;
561  outs=&std::cout;
562  alpha=1.5;
563  iterations=5;
564  mode=mode_importance;
565  chisq=0;
566  bins=bins_max;
567  dim=0;
568  }
569 
570  /// Allocate memory
571  virtual int allocate(size_t ldim) {
572 
573  delx.resize(ldim);
574  d.resize(bins_max*ldim);
575  xi.resize((bins_max+1)*ldim);
576  xin.resize(bins_max+1);
577  weight.resize(bins_max);
578  box.resize(ldim);
579  bin.resize(ldim);
580  x.resize(ldim);
581 
582  dim=ldim;
583 
584  return 0;
585  }
586 
587  /** \brief Integrate function \c func from x=a to x=b.
588 
589  Original documentation from GSL:
590 
591  Normally, <tt>stage = 0</tt> which begins with a new uniform
592  grid and empty weighted average. Calling vegas with <tt>stage
593  = 1</tt> retains the grid from the previous run but discards
594  the weighted average, so that one can "tune" the grid using a
595  relatively small number of points and then do a large run with
596  <tt>stage = 1</tt> on the optimized grid. Setting <tt>stage =
597  2</tt> keeps the grid and the weighted average from the
598  previous run, but may increase (or decrease) the number of
599  histogram bins in the grid depending on the number of calls
600  available. Choosing <tt>stage = 3</tt> enters at the main
601  loop, so that nothing is changed, and is equivalent to
602  performing additional iterations in a previous call.
603 
604  \todo Should stage be passed by reference?
605  \todo There was an update between gsl-1.12 and 1.15 which
606  has not been implemented here yet.
607  */
608  virtual int vegas_minteg_err(int stage, func_t &func, size_t ndim,
609  const vec_t &xl, const vec_t &xu,
610  double &res, double &err) {
611 
612  size_t calls=this->n_points;
613 
614  double cum_int, cum_sig;
615  size_t i, k, it;
616 
617  for (i=0;i<dim;i++) {
618  if (xu[i] <= xl[i]) {
619  std::string serr="Upper limit, "+dtos(xu[i])+", must be greater "+
620  "than lower limit, "+dtos(xl[i])+", in mcarlo_vegas::"+
621  "vegas_minteg_err().";
622  O2SCL_ERR(serr.c_str(),exc_einval);
623  }
624 
625  if (xu[i]-xl[i] > GSL_DBL_MAX) {
626  O2SCL_ERR2("Range of integration is too large, please rescale ",
627  "in mcarlo_vegas::vegas_minteg_err().",exc_einval);
628  }
629  }
630 
631  if (stage == 0) {
632  init_grid(xl,xu,dim);
633  if (this->verbose>=1) {
634  print_lim(xl,xu,dim);
635  }
636  }
637 
638  if (stage<=1) {
639  wtd_int_sum=0;
640  sum_wgts=0;
641  chi_sum=0;
642  it_num=1;
643  samples=0;
644  chisq=0;
645  }
646 
647  if (stage <= 2) {
648 
649  unsigned int lbins=bins_max;
650  unsigned int lboxes=1;
651 
652  if (mode != mode_importance_only) {
653 
654  /* shooting for 2 calls/box */
655 
656  // The original GSL code was:
657  // boxes=floor (pow (calls/2.0, 1.0/dim));
658  // but floor returns double on my machine, so
659  // we explicitly typecast here
660 
661  lboxes=((unsigned int)(floor(pow(calls/2.0,1.0/dim))));
662  mode=mode_importance;
663 
664  if (2*lboxes >= bins_max) {
665  /* if bins/box < 2 */
666  int box_per_bin=GSL_MAX(lboxes/bins_max,1);
667 
668  if (lboxes/box_per_bin<bins_max) lbins=lboxes/box_per_bin;
669  else lbins=bins_max;
670  lboxes=box_per_bin*lbins;
671 
672  mode=mode_stratified;
673  }
674 
675  }
676 
677  //double tot_boxes=gsl_pow_int((double)boxes,dim);
678  double tot_boxes=pow((double)lboxes,(double)dim);
679  calls_per_box=((unsigned int)(GSL_MAX(calls/tot_boxes,2)));
680  calls=((size_t)( calls_per_box*tot_boxes));
681 
682  /* total volume of x-space/(avg num of calls/bin) */
683  jac=vol*pow((double) lbins, (double) dim)/calls;
684 
685  boxes=lboxes;
686 
687  /* If the number of bins changes from the previous invocation, bins
688  are expanded or contracted accordingly, while preserving bin
689  density */
690 
691  if (lbins!=bins) {
692  resize_grid(lbins);
693  if (this->verbose > 2) print_grid(dim);
694  }
695  if (this->verbose >= 1) {
696  print_head(dim,calls,it_num,bins,boxes);
697  }
698  }
699 
700  it_start=it_num;
701 
702  cum_int=0.0;
703  cum_sig=0.0;
704 
705  for (it=0;it<iterations;it++) {
706 
707  double intgrl=0.0, intgrl_sq=0.0;
708  double tss=0.0;
709  double wgt, var, sig;
710  size_t lcalls_per_box=calls_per_box;
711  double jacbin=jac;
712 
713  it_num=it_start+it;
714 
716  init_box_coord(box);
717 
718  do {
719  volatile double m=0, q=0;
720  double f_sq_sum=0.0;
721 
722  for (k=0;k<lcalls_per_box;k++) {
723  double fval, bin_vol;
724 
725  random_point(x,bin,bin_vol,box,xl,xu);
726 
727  fval=func(dim,x);
728  fval*=jacbin*bin_vol;
729 
730  /* recurrence for mean and variance (sum of squares) */
731 
732  {
733  double dt=fval-m;
734  m+=dt/(k+1.0);
735  q+=dt*dt*(k/(k+1.0));
736  }
737 
738  if (mode != mode_stratified) {
739  double f_sq=fval*fval;
740  accumulate_distribution(bin,f_sq);
741  }
742  }
743 
744  intgrl+=m*lcalls_per_box;
745 
746  f_sq_sum=q*lcalls_per_box;
747 
748  tss+=f_sq_sum;
749 
750  if (mode == mode_stratified) {
751  accumulate_distribution (bin, f_sq_sum);
752  }
753 
754  } while (change_box_coord(box));
755 
756  /* Compute final results for this iteration */
757 
758  var=tss/(lcalls_per_box-1.0);
759 
760  if (var>0) {
761  wgt=1.0/var;
762  } else if (sum_wgts>0) {
763  wgt=sum_wgts/samples;
764  } else {
765  wgt=0.0;
766  }
767 
768  intgrl_sq=intgrl*intgrl;
769 
770  sig=sqrt(var);
771 
772  result=intgrl;
773  sigma=sig;
774 
775  if (wgt > 0.0) {
776  double lsum_wgts=sum_wgts;
777  double m=(sum_wgts > 0) ? (wtd_int_sum/sum_wgts) : 0;
778  double q=intgrl-m;
779 
780  samples++;
781  sum_wgts+=wgt;
782  wtd_int_sum+=intgrl*wgt;
783  chi_sum+=intgrl_sq*wgt;
784 
785  cum_int= wtd_int_sum/sum_wgts;
786  cum_sig=sqrt(1/sum_wgts);
787 
788  /* The original chisq formula from the Lepage paper is
789 
790  if ( samples > 1) {
791  chisq=(chi_sum-wtd_int_sum*cum_int)/(samples-1.0);
792  }
793 
794  This can suffer from cancellations and return a negative
795  value of chi squared. We use the new formula from
796  GSL-1.12 instead
797  */
798  if (samples==1) {
799  chisq=0;
800  } else {
801  chisq*=(samples-2.0);
802  chisq+=(wgt/(1+(wgt/lsum_wgts)))*q*q;
803  chisq/=(samples-1.0);
804  }
805 
806  } else {
807  cum_int+=(intgrl-cum_int)/(it+1.0);
808  cum_sig=0.0;
809  }
810 
811  if (this->verbose >= 1) {
812  print_res(it_num,intgrl,sig,cum_int,cum_sig,chisq);
813  if (it+1 == iterations && this->verbose > 1) {
814  print_grid(dim);
815  }
816  }
817 
818  if (this->verbose > 2) {
819  print_dist(dim);
820  }
821 
822  refine_grid ();
823 
824  if (this->verbose > 2) {
825  print_grid(dim);
826  }
827 
828  }
829 
830  /*
831  By setting stage to 1 further calls will generate independent
832  estimates based on the same grid, although it may be rebinned.
833  */
834  stage=1;
835 
836  res=cum_int;
837  err=cum_sig;
838 
839  return GSL_SUCCESS;
840  }
841 
842  virtual ~mcarlo_vegas() {}
843 
844  /// Integrate function \c func from x=a to x=b.
845  virtual int minteg_err(func_t &func, size_t ndim, const vec_t &a,
846  const vec_t &b, double &res, double &err) {
847  allocate(ndim);
848  chisq=0;
849  bins=bins_max;
850  int ret=vegas_minteg_err(0,func,ndim,a,b,res,err);
851  return ret;
852  }
853 
854  /** \brief Integrate function \c func over the hypercube from
855  \f$ x_i=a_i \f$ to \f$ x_i=b_i \f$ for
856  \f$ 0<i< \f$ ndim-1
857  */
858  virtual double minteg(func_t &func, size_t ndim, const vec_t &a,
859  const vec_t &b) {
860  double res;
861  minteg_err(func,ndim,a,b,res,this->interror);
862  return res;
863  }
864 
865  /// Return string denoting type ("mcarlo_vegas")
866  virtual const char *type() { return "mcarlo_vegas"; }
867 
868  };
869 
870 #ifndef DOXYGEN_NO_O2NS
871 }
872 #endif
873 
874 #endif
875 
ubvector_int box
The boxes for each direction.
Definition: mcarlo_vegas.h:187
rng_gsl_uniform_real rng_dist
The random number distribution.
Definition: mcarlo.h:67
virtual void print_head(unsigned long num_dim, unsigned long calls, unsigned int lit_num, unsigned int lbins, unsigned int lboxes)
Print header.
Definition: mcarlo_vegas.h:469
ubvector delx
The iteration length in each direction.
Definition: mcarlo_vegas.h:178
double chisq
The chi-squared per degree of freedom for the weighted estimate of the integral.
Definition: mcarlo_vegas.h:155
ubvector_int bin
The bins for each direction.
Definition: mcarlo_vegas.h:185
int change_box_coord(ubvector_int &boxt)
Change box coordinates.
Definition: mcarlo_vegas.h:230
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
ubvector xin
Storage for grid refinement.
Definition: mcarlo_vegas.h:176
unsigned int it_start
The starting iteration number.
Definition: mcarlo_vegas.h:203
virtual int minteg_err(func_t &func, size_t ndim, const vec_t &a, const vec_t &b, double &res, double &err)
Integrate function func from x=a to x=b.
Definition: mcarlo_vegas.h:845
void random_point(vec_t &lx, ubvector_int &lbin, double &bin_vol, const ubvector_int &lbox, const vec_t &xl, const vec_t &xu)
Generate a random position in a given box.
Definition: mcarlo_vegas.h:301
virtual void print_grid(unsigned long ldim)
Print grid.
Definition: mcarlo_vegas.h:528
virtual void refine_grid()
Refine the grid.
Definition: mcarlo_vegas.h:383
invalid argument supplied by user
Definition: err_hnd.h:59
virtual double minteg(func_t &func, size_t ndim, const vec_t &a, const vec_t &b)
Integrate function func over the hypercube from to for ndim-1.
Definition: mcarlo_vegas.h:858
virtual int vegas_minteg_err(int stage, func_t &func, size_t ndim, const vec_t &xl, const vec_t &xu, double &res, double &err)
Integrate function func from x=a to x=b.
Definition: mcarlo_vegas.h:608
void accumulate_distribution(ubvector_int &lbin, double y)
Add the most recently generated result to the distribution.
Definition: mcarlo_vegas.h:282
virtual void print_dist(unsigned long ldim)
Print distribution.
Definition: mcarlo_vegas.h:504
unsigned int iterations
Set the number of iterations (default 5)
Definition: mcarlo_vegas.h:144
ubvector d
Distribution.
Definition: mcarlo_vegas.h:190
size_t dim
Number of dimensions.
Definition: mcarlo_vegas.h:168
unsigned int it_num
The total number of iterations.
Definition: mcarlo_vegas.h:205
virtual void reset_grid_values()
Reset grid values.
Definition: mcarlo_vegas.h:266
virtual void init_grid(const vec_t &xl, const vec_t &xu, size_t ldim)
Initialize grid.
Definition: mcarlo_vegas.h:248
int verbose
Verbosity.
Definition: inte_multi.h:64
ubvector weight
The weight for each bin.
Definition: mcarlo_vegas.h:180
#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.
virtual const char * type()
Return string denoting type ("mcarlo_vegas")
Definition: mcarlo_vegas.h:866
virtual void resize_grid(unsigned int lbins)
Resize the grid.
Definition: mcarlo_vegas.h:346
unsigned long n_points
Number of integration points (default 1000)
Definition: mcarlo.h:63
virtual int allocate(size_t ldim)
Allocate memory.
Definition: mcarlo_vegas.h:571
double result
Result from last iteration.
Definition: mcarlo_vegas.h:132
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
unsigned int calls_per_box
Number of function calls per box.
Definition: mcarlo_vegas.h:209
Monte-Carlo integration [abstract base].
Definition: mcarlo.h:51
double vol
The volume of the current bin.
Definition: mcarlo_vegas.h:182
virtual void init_box_coord(ubvector_int &boxt)
Initialize box coordinates.
Definition: mcarlo_vegas.h:212
static const size_t bins_max
Maximum number of bins.
Definition: mcarlo_vegas.h:165
virtual void print_res(unsigned int itr, double res, double err, double cum_res, double cum_err, double chi_sq)
Print results.
Definition: mcarlo_vegas.h:486
vec_t x
Point for function evaluation.
Definition: mcarlo_vegas.h:553
double interror
The uncertainty for the last integration computation.
Definition: inte_multi.h:110
Multidimensional integration using Vegas Monte Carlo (GSL)
Definition: mcarlo_vegas.h:115
virtual void print_lim(const vec_t &xl, const vec_t &xu, unsigned long ldim)
Print limits of integration.
Definition: mcarlo_vegas.h:453
double sigma
Uncertainty from last iteration.
Definition: mcarlo_vegas.h:135
unsigned int boxes
Number of boxes.
Definition: mcarlo_vegas.h:172
std::ostream * outs
The output stream to send output information (default std::cout)
Definition: mcarlo_vegas.h:158
unsigned int bins
Number of bins.
Definition: mcarlo_vegas.h:170
std::function< double(size_t, const boost::numeric::ublas::vector< double > &)> multi_funct11
Multi-dimensional function typedef.
Definition: multi_funct.h:45
unsigned int samples
Number of samples for computing chi squared.
Definition: mcarlo_vegas.h:207
ubvector xi
Boundaries for each bin.
Definition: mcarlo_vegas.h:174
double alpha
The stiffness of the rebinning algorithm (default 1.5)
Definition: mcarlo_vegas.h:141
rng_t rng
The random number generator.
Definition: mcarlo.h:70

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