mcarlo_miser.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/miser.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_MISER_H
43 #define O2SCL_MCARLO_MISER_H
44 
45 /** \file mcarlo_miser.h
46  \brief File defining \ref o2scl::mcarlo_miser
47 */
48 #include <iostream>
49 
50 #include <boost/numeric/ublas/vector.hpp>
51 
52 #include <o2scl/misc.h>
53 #include <o2scl/mcarlo.h>
54 
55 #include <gsl/gsl_math.h>
56 #include <gsl/gsl_monte.h>
57 #include <gsl/gsl_machine.h>
58 #include <gsl/gsl_monte_miser.h>
59 
60 #ifndef DOXYGEN_NO_O2NS
61 namespace o2scl {
62 #endif
63 
64  /** \brief Multidimensional integration using the MISER Monte Carlo
65  algorithm (GSL)
66 
67  This class uses recursive stratified sampling to estimate the
68  value of an integral over a hypercubic region.
69 
70  By default the minimum number of calls to estimate the variance
71  is 16 times the number of dimensions. This ratio is
72  stored in \ref calls_per_dim. By default the minimum number of
73  calls per bisection is 32 times \ref calls_per_dim times
74  the number of dimensions. This ratio is stored in
75  \ref bisection_ratio. These ratios are employed by
76  \ref minteg_err().
77 
78  Alternatively, the user can directly set these minimums by \ref
79  set_min_calls() and \ref set_min_calls_per_bisection() and use
80  \ref miser_minteg_err() which ignores \ref calls_per_dim and
81  \ref bisection_ratio.
82 
83  If \ref mcarlo::verbose is greater than zero, then the
84  status of the integration is output at every level of bisection
85  less than \ref n_levels_out. If it is greater than 1, then the
86  boundaries of the current region are also output. Finally, if it
87  is greater than 2, a keypress is required after each output.
88 
89  Based on \ref Press90 .
90  */
91  template<class func_t=multi_funct11,
93  class rng_t=rng_gsl>
94  class mcarlo_miser : public mcarlo<func_t,vec_t,rng_t> {
95 
96  public:
97 
100 
101  /** \brief Number of calls per dimension (default 16)
102  */
104 
105  /** \brief Factor to set \ref min_calls_per_bisection (default 32)
106  */
108 
109  /** \brief Introduce random variation into bisection (default 0.0)
110 
111  From GSL documentation:
112  \verbatim
113  This parameter introduces a random fractional variation of size
114  DITHER into each bisection, which can be used to break the
115  symmetry of integrands which are concentrated near the exact
116  center of the hypercubic integration region. The default value of
117  dither is zero, so no variation is introduced. If needed, a
118  typical value of DITHER is 0.1.
119  \endverbatim
120  */
121  double dither;
122 
123  /** \brief Specify fraction of function calls for estimating variance
124  (default 0.1)
125 
126  From GSL documentation:
127  \verbatim
128  This parameter specifies the fraction of the currently available
129  number of function calls which are allocated to estimating the
130  variance at each recursive step. The default value is 0.1.
131  \endverbatim
132  */
134 
135  /** \brief How estimated variances for two sub-regions are combined
136  (default 2.0)
137 
138  The error handler will be called if this is less than zero.
139 
140  From GSL documentation:
141  \verbatim
142  This parameter controls how the estimated variances for the two
143  sub-regions of a bisection are combined when allocating points.
144  With recursive sampling the overall variance should scale better
145  than 1/N, since the values from the sub-regions will be obtained
146  using a procedure which explicitly minimizes their variance. To
147  accommodate this behavior the MISER algorithm allows the total
148  variance to depend on a scaling parameter \alpha,
149 
150  \Var(f) = {\sigma_a \over N_a^\alpha} + {\sigma_b \over N_b^\alpha}.
151 
152  The authors of the original paper describing MISER recommend the
153  value \alpha = 2 as a good choice, obtained from numerical
154  experiments, and this is used as the default value in this
155  implementation.
156  \endverbatim
157  */
158  double alpha;
159 
160  /** \brief Minimum number of calls to estimate the variance
161 
162  This is set by minteg() and minteg_err() to be \ref
163  calls_per_dim times the number of dimensions in the problem. The
164  default value of calls_per_dim is 16 (which is the GSL default).
165 
166  From GSL documentation:
167  \verbatim
168  This parameter specifies the minimum number of function calls
169  required for each estimate of the variance. If the number of
170  function calls allocated to the estimate using ESTIMATE_FRAC falls
171  below MIN_CALLS then MIN_CALLS are used instead. This ensures
172  that each estimate maintains a reasonable level of accuracy.
173  \endverbatim
174  */
175  int set_min_calls(size_t &mc) {
176  min_calls=mc;
177  return 0;
178  }
179 
180  /** \brief Minimum number of calls required to proceed with bisection
181 
182  This is set by minteg() and minteg_err() to be \ref
183  calls_per_dim times \ref bisection_ratio times the number of
184  dimensions in the problem. The default values give 512 times the
185  number of dimensions (also the GSL default).
186 
187  From GSL documentation:
188  \verbatim
189  This parameter specifies the minimum number of function calls
190  required to proceed with a bisection step. When a recursive step
191  has fewer calls available than MIN_CALLS_PER_BISECTION it performs
192  a plain Monte Carlo estimate of the current sub-region and
193  terminates its branch of the recursion.
194  \endverbatim
195  */
196  int set_min_calls_per_bisection(size_t &mcb) {
198  return 0;
199  }
200 
201  /** \brief The number of recursive levels to output when verbose is
202  greater than zero (default 5)
203  */
204  size_t n_levels_out;
205 
206 #ifndef DOXYGEN_INTERNAL
207 
208  protected:
209 
210  /// Minimum number of calls to estimate the variance
211  size_t min_calls;
212 
213  /// Minimum number of calls required to proceed with bisection
215 
216  /// The number of dimensions of memory allocated
217  size_t dim;
218 
219  /// \name Arrays which contain a value for each dimension
220  //@{
221  /// The current midpoint
222  ubvector xmid;
223  /// The left variance
224  ubvector sigma_l;
225  /// The right variance
226  ubvector sigma_r;
227  /// The maximum function value in the left half
228  ubvector fmax_l;
229  /// The maximum function value in the right half
230  ubvector fmax_r;
231  /// The minimum function value in the left half
232  ubvector fmin_l;
233  /// The minimum function value in the right half
234  ubvector fmin_r;
235  /// The sum in the left half
236  ubvector fsum_l;
237  /// The sum in the right half
238  ubvector fsum_r;
239  /// The sum of the squares in the left half
240  ubvector fsum2_l;
241  /// The sum of the squares in the right half
242  ubvector fsum2_r;
243  /// The number of evaluation points in the left half
244  ubvector_size_t hits_l;
245  /// The number of evaluation points in the right half
246  ubvector_size_t hits_r;
247  //@}
248 
249  /** \brief Estimate the variance
250 
251  \future Remove the reference to GSL_POSINF and replace with a
252  function parameter.
253  */
254  virtual int estimate_corrmc(func_t &func, size_t ndim,
255  const vec_t &xl, const vec_t &xu,
256  size_t calls, double &res,
257  double &err, const ubvector &lxmid,
258  ubvector &lsigma_l, ubvector &lsigma_r) {
259  size_t i, n;
260 
261  double m=0.0, q=0.0;
262  double vol=1.0;
263 
264  for (i=0;i<dim;i++) {
265  vol*=xu[i]-xl[i];
266  hits_l[i]=hits_r[i]=0;
267  fsum_l[i]=fsum_r[i]=0.0;
268  fsum2_l[i]=fsum2_r[i]=0.0;
269  lsigma_l[i]=lsigma_r[i]=-1;
270  }
271 
272  for (n=0;n<calls;n++) {
273  double fval;
274 
275  unsigned int j=(n/2) % dim;
276  unsigned int side=(n % 2);
277 
278  for (i=0;i<dim;i++) {
279 
280  // The equivalent of gsl_rng_uniform_pos()
281  double z;
282  do {
283  z=this->rng_dist(this->rng);
284  } while (z==0);
285 
286  if (i != j) {
287  x[i]=xl[i]+z*(xu[i]-xl[i]);
288  } else {
289  if (side == 0) {
290  x[i]=lxmid[i]+z*(xu[i]-lxmid[i]);
291  } else {
292  x[i]=xl[i]+z*(lxmid[i]-xl[i]);
293  }
294  }
295  }
296 
297  fval=func(ndim,x);
298 
299  /* recurrence for mean and variance */
300  {
301  double d=fval-m;
302  m+=d/(n+1.0);
303  q+=d*d*(n/(n+1.0));
304  }
305 
306  /* compute the variances on each side of the bisection */
307  for (i=0;i<dim;i++) {
308  if (x[i] <= lxmid[i]) {
309  fsum_l[i]+=fval;
310  fsum2_l[i]+=fval*fval;
311  hits_l[i]++;
312  } else {
313  fsum_r[i]+=fval;
314  fsum2_r[i]+=fval*fval;
315  hits_r[i]++;
316  }
317  }
318  }
319 
320  for (i=0;i<dim;i++) {
321  double fraction_l=(lxmid[i]-xl[i])/(xu[i]-xl[i]);
322 
323  if (hits_l[i] > 0) {
324  fsum_l[i] /= hits_l[i];
325  lsigma_l[i]=sqrt(fsum2_l[i]-fsum_l[i]*fsum_l[i]/hits_l[i]);
326  lsigma_l[i]*=fraction_l*vol/hits_l[i];
327  }
328 
329  if (hits_r[i] > 0) {
330  fsum_r[i] /= hits_r[i];
331  lsigma_r[i]=sqrt(fsum2_r[i]-fsum_r[i]*fsum_r[i]/hits_r[i]);
332  lsigma_r[i]*=(1-fraction_l)*vol/hits_r[i];
333  }
334  }
335 
336  res=vol*m;
337 
338  if (calls<2) {
339  err=GSL_POSINF;
340  } else {
341  err=vol*sqrt(q/(calls*(calls-1.0)));
342  }
343 
344  return success;
345  }
346 
347  /// The most recent integration point
348  vec_t x;
349 
350 #endif
351 
352  public:
353 
354  mcarlo_miser() {
355  estimate_frac=0.1;
356  alpha=2.0;
357  dither=0.0;
358  min_calls=100;
359  min_calls_per_bisection=3000;
360  dim=0;
361  n_levels_out=5;
362  calls_per_dim=16;
363  bisection_ratio=32;
364  }
365 
366  virtual ~mcarlo_miser() {
367  }
368 
369  /// Allocate memory
370  virtual int allocate(size_t ldim) {
371 
372  if (ldim==0) {
373  O2SCL_ERR2("Can't allocate zero memory in ",
374  "mcarlo_miser::allocate().",exc_efailed);
375  }
376 
377  if (ldim!=dim) {
378  x.resize(ldim);
379  xmid.resize(ldim);
380  sigma_l.resize(ldim);
381  sigma_r.resize(ldim);
382  fmax_l.resize(ldim);
383  fmax_r.resize(ldim);
384  fmin_l.resize(ldim);
385  fmin_r.resize(ldim);
386  fsum_l.resize(ldim);
387  fsum_r.resize(ldim);
388  fsum2_l.resize(ldim);
389  fsum2_r.resize(ldim);
390  hits_l.resize(ldim);
391  hits_r.resize(ldim);
392  }
393 
394  dim=ldim;
395 
396  return 0;
397  }
398 
399  /** \brief Integrate function \c func over the hypercube from
400  \f$ x_i=\mathrm{xl}_i \f$ to \f$ x_i=\mathrm{xu}_i \f$ for
401  \f$ 0<i< \f$ ndim-1
402 
403  \note The values of \ref min_calls and \ref
404  min_calls_per_bisection should be set before calling this
405  function. The default values if not set are 100 and 3000
406  respectively, which correspond to the GSL default setting
407  for a 6 dimensional problem.
408  */
409  virtual int miser_minteg_err(func_t &func, size_t ndim, const vec_t &xl,
410  const vec_t &xu, size_t calls, size_t level,
411  double &res, double &err) {
412 
413  if (min_calls==0 || min_calls_per_bisection==0) {
414  O2SCL_ERR2("Variables min_calls or min_calls_per_bisection ",
415  "are zero in mcarlo_miser::miser_minteg_err().",
416  exc_einval);
417  }
418 
419  size_t n, estimate_calls, calls_l, calls_r;
420  size_t i;
421  size_t i_bisect;
422  int found_best;
423 
424  double res_est=0, err_est=0;
425  double res_r=0, err_r=0, res_l=0, err_l=0;
426  double xbi_l, xbi_m, xbi_r, s;
427 
428  double vol;
429  double weight_l, weight_r;
430 
431  for (i=0;i<dim;i++) {
432  if (xu[i] <= xl[i]) {
433  std::string str="Upper limit, "+dtos(xu[i])+", must be greater "+
434  "than lower limit, "+dtos(xl[i])+", in mcarlo_miser::"+
435  "miser_minteg_err().";
436  O2SCL_ERR(str.c_str(),exc_einval);
437  }
438  if (xu[i]-xl[i]>GSL_DBL_MAX) {
439  O2SCL_ERR2("Range of integration is too large ",
440  "in mcarlo_miser::miser_minteg_err().",exc_einval);
441  }
442  }
443 
444  if (alpha<0) {
445  std::string str="Parameter alpha, "+dtos(alpha)+", must be non-"+
446  "negative in mcarlo_miser::mister_minteg_err().";
447  O2SCL_ERR(str.c_str(),exc_einval);
448  }
449 
450  /* [GSL] Compute volume */
451 
452  vol=1;
453 
454  for (i=0;i<dim;i++) {
455  vol*=xu[i]-xl[i];
456  }
457 
458  if (calls<min_calls_per_bisection) {
459  double m=0.0, q=0.0;
460 
461  if (calls<2) {
462  O2SCL_ERR2("Insufficient calls for subvolume ",
463  "in mcarlo_miser::miser_minteg_err().",exc_einval);
464  }
465 
466  for (n=0;n<calls;n++) {
467  /* Choose a random point in the integration region */
468 
469  for (i=0;i<dim;i++) {
470 
471  // The equivalent of gsl_rng_uniform_pos()
472  double rdn;
473  do {
474  rdn=this->rng_dist(this->rng);
475  } while (rdn==0);
476 
477  x[i]=xl[i]+rdn*(xu[i]-xl[i]);
478  }
479 
480  {
481  double fval;
482  fval=func(ndim,x);
483 
484  /* [GSL] recurrence for mean and variance */
485 
486  double d=fval-m;
487  m+=d/(n+1.0);
488  q+=d*d*(n/(n+1.0));
489  }
490  }
491 
492  res=vol*m;
493 
494  err=vol*sqrt(q/(calls*(calls-1.0)));
495 
496  return success;
497  }
498 
499  // The equivalent of what's in GSL but rewritten with explicit
500  // typecasts
501  size_t prod=(size_t)(((double)calls)*estimate_frac);
502  estimate_calls=(min_calls > prod ? min_calls : prod);
503 
504  if (estimate_calls<4*dim) {
505  O2SCL_ERR2("Insufficient calls to sample all halfspaces ",
506  "in mcarlo_miser::miser_minteg_err().",exc_esanity);
507  }
508 
509  // [GSL] Flip coins to bisect the integration region with some fuzz
510 
511  for (i=0;i<dim;i++) {
512  s=(this->rng_dist(this->rng)-0.5) >= 0.0 ? dither : -dither;
513  xmid[i]=(0.5+s)*xl[i]+(0.5-s)*xu[i];
514  }
515 
516  /*
517  [GSL] The idea is to chose the direction to bisect based on
518  which will give the smallest total variance. We could (and may
519  do so later) use MC to compute these variances. But the NR guys
520  simply estimate the variances by finding the min and max
521  function values for each half-region for each bisection.
522  */
523  estimate_corrmc(func,dim,xl,xu,estimate_calls,
524  res_est,err_est,xmid,sigma_l,sigma_r);
525 
526  // [GSL] We have now used up some calls for the estimation
527 
528  calls -= estimate_calls;
529 
530  // [GSL] Now find direction with the smallest total "variance"
531 
532  {
533  double best_var=GSL_DBL_MAX;
534  double beta=2.0/(1.0+alpha);
535  found_best=0;
536  i_bisect=0;
537  weight_l=weight_r=1.0;
538 
539  for (i=0;i<dim;i++) {
540  if (sigma_l[i] >= 0 && sigma_r[i] >= 0) {
541  // [GSL] Estimates are okay
542  double var=pow (sigma_l[i], beta)+pow (sigma_r[i], beta);
543 
544  if (var <= best_var) {
545  found_best=1;
546  best_var=var;
547  i_bisect=i;
548  weight_l=pow (sigma_l[i], beta);
549  weight_r=pow (sigma_r[i], beta);
550  if (weight_l==0 && weight_r==0) {
551  weight_l=1;
552  weight_r=1;
553  }
554  }
555  } else {
556  if (sigma_l[i]<0) {
557  O2SCL_ERR2("No points in left-half space ",
558  "in mcarlo_miser::miser_minteg_err().",
559  exc_esanity);
560  }
561  if (sigma_r[i]<0) {
562  O2SCL_ERR2("No points in right-half space ",
563  "in mcarlo_miser::miser_minteg_err().",
564  exc_esanity);
565  }
566  }
567  }
568  }
569 
570  if (!found_best) {
571  // [GSL] All estimates were the same, so chose a direction at
572  // random
573  i_bisect=((int)(this->rng_dist(this->rng)*(dim-1.0e-10)));
574  //std::uniform_int_distribution<int> int_dist(0,dim-1);
575  //i_bisect=int_dist(this->rng);
576  //gsl_rnga gr;
577  //i_bisect=gr.random_int(dim);
578  }
579 
580  xbi_l=xl[i_bisect];
581  xbi_m=xmid[i_bisect];
582  xbi_r=xu[i_bisect];
583 
584  // [GSL] Get the actual fractional sizes of the two "halves", and
585  // distribute the remaining calls among them
586  {
587  double fraction_l=fabs ((xbi_m-xbi_l)/(xbi_r-xbi_l));
588  double fraction_r=1-fraction_l;
589 
590  double a=fraction_l*weight_l;
591  double b=fraction_r*weight_r;
592 
593  calls_l=(size_t)(min_calls+(calls-2*min_calls)*a/(a+b));
594  calls_r=(size_t)(min_calls+(calls-2*min_calls)*b/(a+b));
595  }
596 
597  if (this->verbose>0 && level<n_levels_out) {
598  std::cout << "mcarlo_miser: level,calls_l,calls_r,calls,min_calls_pb: "
599  << level << " " << calls_l << " " << calls_r << " " << calls << " "
600  << min_calls_per_bisection << std::endl;
601  std::cout << "\tres,err: " << res_est << " " << err_est << std::endl;
602  if (this->verbose>1) {
603  std::cout << "\ti,left,mid,right: " << i_bisect << " "
604  << xbi_l << " " << xbi_m << " " << xbi_r << std::endl;
605  for(size_t j=0;j<dim;j++) {
606  std::cout << "\t\ti,low,high: " << j << " " << xl[j] << " " << xu[j]
607  << std::endl;
608  }
609  }
610  if (this->verbose>2) {
611  char ch;
612  std::cin >> ch;
613  }
614  }
615 
616  /* [GSL] Compute the integral for the left hand side of the
617  bisection. Due to the recursive nature of the algorithm we must
618  allocate some new memory for the integration limits for each
619  recursive call
620  */
621  {
622  int status;
623 
624  vec_t xu_tmp(dim);
625 
626  for (i=0;i<dim;i++) {
627  xu_tmp[i]=xu[i];
628  }
629 
630  xu_tmp[i_bisect]=xbi_m;
631 
632  status=miser_minteg_err(func,dim,xl,xu_tmp,calls_l,level+1,
633  res_l,err_l);
634 
635  if (status != success) {
636  return status;
637  }
638  }
639 
640  // [GSL] Compute the integral for the right hand side of the
641  // bisection
642  {
643  int status;
644 
645  vec_t xl_tmp(dim);
646 
647  for (i=0;i<dim;i++) {
648  xl_tmp[i]=xl[i];
649  }
650 
651  xl_tmp[i_bisect]=xbi_m;
652 
653  status=miser_minteg_err(func,dim,xl_tmp,xu,calls_r,level+1,
654  res_r,err_r);
655 
656  if (status != success) {
657  return status;
658  }
659  }
660 
661  res=res_l+res_r;
662  err=sqrt(err_l*err_l+err_r*err_r);
663 
664  return 0;
665  }
666 
667  /** \brief Integrate function \c func from x=a to x=b.
668 
669  This function is just a wrapper to miser_minteg_err() which
670  allocates the memory if necessary, sets \c min_calls and \c
671  min_calls_per_bisection, calls \ref miser_minteg_err(), and then
672  frees the previously allocated memory.
673  */
674  virtual int minteg_err(func_t &func, size_t ndim, const vec_t &a,
675  const vec_t &b, double &res, double &err) {
676  if (ndim!=dim) allocate(ndim);
677  min_calls=calls_per_dim*ndim;
678  min_calls_per_bisection=bisection_ratio*min_calls;
679  int ret=miser_minteg_err(func,ndim,a,b,this->n_points,0,res,err);
680  // Set these to back to zero to ensure that if
681  // miser_minteg_err() is called directly, the user sets
682  // min_calls and min_calls_per_bisection first.
683  min_calls=0;
684  min_calls_per_bisection=0;
685  return ret;
686  }
687 
688  /** \brief Integrate function \c func over the hypercube from
689  \f$ x_i=a_i \f$ to \f$ x_i=b_i \f$ for
690  \f$ 0<i< \f$ ndim-1
691 
692  This function is just a wrapper to minteg_err() which allocates
693  the memory, sets min_calls and min_calls_per_bisection, calls
694  miser_minteg_err(), and then frees the previously allocated
695  memory.
696  */
697  virtual double minteg(func_t &func, size_t ndim, const vec_t &a,
698  const vec_t &b) {
699  double res;
700  minteg_err(func,ndim,a,b,res,this->interror);
701  return res;
702  }
703 
704  /// Return string denoting type ("mcarlo_miser")
705  virtual const char *type() { return "mcarlo_miser"; }
706 
707  };
708 
709 #ifndef DOXYGEN_NO_O2NS
710 }
711 #endif
712 
713 #endif
714 
rng_gsl_uniform_real rng_dist
The random number distribution.
Definition: mcarlo.h:67
ubvector xmid
The current midpoint.
Definition: mcarlo_miser.h:222
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
size_t bisection_ratio
Factor to set min_calls_per_bisection (default 32)
Definition: mcarlo_miser.h:107
double dither
Introduce random variation into bisection (default 0.0)
Definition: mcarlo_miser.h:121
sanity check failed - shouldn&#39;t happen
Definition: err_hnd.h:65
ubvector fsum2_l
The sum of the squares in the left half.
Definition: mcarlo_miser.h:240
invalid argument supplied by user
Definition: err_hnd.h:59
int set_min_calls_per_bisection(size_t &mcb)
Minimum number of calls required to proceed with bisection.
Definition: mcarlo_miser.h:196
generic failure
Definition: err_hnd.h:61
int set_min_calls(size_t &mc)
Minimum number of calls to estimate the variance.
Definition: mcarlo_miser.h:175
ubvector fmin_r
The minimum function value in the right half.
Definition: mcarlo_miser.h:234
virtual const char * type()
Return string denoting type ("mcarlo_miser")
Definition: mcarlo_miser.h:705
size_t calls_per_dim
Number of calls per dimension (default 16)
Definition: mcarlo_miser.h:103
int verbose
Verbosity.
Definition: inte_multi.h:64
size_t min_calls
Minimum number of calls to estimate the variance.
Definition: mcarlo_miser.h:211
#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.
ubvector sigma_r
The right variance.
Definition: mcarlo_miser.h:226
virtual int allocate(size_t ldim)
Allocate memory.
Definition: mcarlo_miser.h:370
unsigned long n_points
Number of integration points (default 1000)
Definition: mcarlo.h:63
vec_t x
The most recent integration point.
Definition: mcarlo_miser.h:348
double alpha
How estimated variances for two sub-regions are combined (default 2.0)
Definition: mcarlo_miser.h:158
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
Monte-Carlo integration [abstract base].
Definition: mcarlo.h:51
size_t n_levels_out
The number of recursive levels to output when verbose is greater than zero (default 5) ...
Definition: mcarlo_miser.h:204
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_miser.h:697
ubvector fmin_l
The minimum function value in the left half.
Definition: mcarlo_miser.h:232
ubvector_size_t hits_l
The number of evaluation points in the left half.
Definition: mcarlo_miser.h:244
ubvector fsum2_r
The sum of the squares in the right half.
Definition: mcarlo_miser.h:242
size_t dim
The number of dimensions of memory allocated.
Definition: mcarlo_miser.h:217
double interror
The uncertainty for the last integration computation.
Definition: inte_multi.h:110
Multidimensional integration using the MISER Monte Carlo algorithm (GSL)
Definition: mcarlo_miser.h:94
ubvector sigma_l
The left variance.
Definition: mcarlo_miser.h:224
ubvector fsum_r
The sum in the right half.
Definition: mcarlo_miser.h:238
ubvector fsum_l
The sum in the left half.
Definition: mcarlo_miser.h:236
double estimate_frac
Specify fraction of function calls for estimating variance (default 0.1)
Definition: mcarlo_miser.h:133
ubvector fmax_r
The maximum function value in the right half.
Definition: mcarlo_miser.h:230
virtual int miser_minteg_err(func_t &func, size_t ndim, const vec_t &xl, const vec_t &xu, size_t calls, size_t level, double &res, double &err)
Integrate function func over the hypercube from to for ndim-1.
Definition: mcarlo_miser.h:409
size_t min_calls_per_bisection
Minimum number of calls required to proceed with bisection.
Definition: mcarlo_miser.h:214
ubvector_size_t hits_r
The number of evaluation points in the right half.
Definition: mcarlo_miser.h:246
std::function< double(size_t, const boost::numeric::ublas::vector< double > &)> multi_funct11
Multi-dimensional function typedef.
Definition: multi_funct.h:45
ubvector fmax_l
The maximum function value in the left half.
Definition: mcarlo_miser.h:228
virtual int estimate_corrmc(func_t &func, size_t ndim, const vec_t &xl, const vec_t &xu, size_t calls, double &res, double &err, const ubvector &lxmid, ubvector &lsigma_l, ubvector &lsigma_r)
Estimate the variance.
Definition: mcarlo_miser.h:254
Success.
Definition: err_hnd.h:47
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_miser.h:674
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).