Public Types | Public Member Functions | Public Attributes | Protected Member Functions | Protected Attributes | Static Protected Attributes | List of all members
o2scl::mcarlo_vegas< func_t, vec_t, rng_t > Class Template Reference

Multidimensional integration using Vegas Monte Carlo (GSL) More...

#include <mcarlo_vegas.h>

Inheritance diagram for o2scl::mcarlo_vegas< func_t, vec_t, rng_t >:
o2scl::mcarlo< func_t, vec_t, rng_t > o2scl::inte_multi< func_t, vec_t >

Public Types

typedef boost::numeric::ublas::vector< double > ubvector
 
typedef boost::numeric::ublas::vector< size_t > ubvector_size_t
 
typedef boost::numeric::ublas::vector< int > ubvector_int
 

Public Member Functions

virtual int allocate (size_t ldim)
 Allocate memory.
 
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. More...
 
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.
 
virtual double minteg (func_t &func, size_t ndim, const vec_t &a, const vec_t &b)
 Integrate function func over the hypercube from $ x_i=a_i $ to $ x_i=b_i $ for $ 0<i< $ ndim-1.
 
virtual const char * type ()
 Return string denoting type ("mcarlo_vegas")
 
- Public Member Functions inherited from o2scl::inte_multi< func_t, vec_t >
double get_error ()
 Return the error in the result from the last call to minteg() or minteg_err() More...
 
const char * type ()
 Return string denoting type ("inte_multi")
 

Public Attributes

double result
 Result from last iteration.
 
double sigma
 Uncertainty from last iteration.
 
double alpha
 The stiffness of the rebinning algorithm (default 1.5) More...
 
unsigned int iterations
 Set the number of iterations (default 5)
 
double chisq
 The chi-squared per degree of freedom for the weighted estimate of the integral. More...
 
std::ostream * outs
 The output stream to send output information (default std::cout)
 
- Public Attributes inherited from o2scl::mcarlo< func_t, vec_t, rng_t >
unsigned long n_points
 Number of integration points (default 1000)
 
rng_gsl_uniform_real rng_dist
 The random number distribution.
 
rng_t rng
 The random number generator.
 
- Public Attributes inherited from o2scl::inte_multi< func_t, vec_t >
bool err_nonconv
 If true, call the error handler if the routine does not "converge".
 
int verbose
 Verbosity.
 
double tol_rel
 The maximum "uncertainty" in the value of the integral (default $ 10^{-8} $).
 

Protected Member Functions

virtual void init_box_coord (ubvector_int &boxt)
 Initialize box coordinates.
 
int change_box_coord (ubvector_int &boxt)
 Change box coordinates. More...
 
virtual void init_grid (const vec_t &xl, const vec_t &xu, size_t ldim)
 Initialize grid.
 
virtual void reset_grid_values ()
 Reset grid values.
 
void accumulate_distribution (ubvector_int &lbin, double y)
 Add the most recently generated result to the distribution. More...
 
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. More...
 
virtual void resize_grid (unsigned int lbins)
 Resize the grid.
 
virtual void refine_grid ()
 Refine the grid.
 
virtual void print_lim (const vec_t &xl, const vec_t &xu, unsigned long ldim)
 Print limits of integration.
 
virtual void print_head (unsigned long num_dim, unsigned long calls, unsigned int lit_num, unsigned int lbins, unsigned int lboxes)
 Print header.
 
virtual void print_res (unsigned int itr, double res, double err, double cum_res, double cum_err, double chi_sq)
 Print results.
 
virtual void print_dist (unsigned long ldim)
 Print distribution.
 
virtual void print_grid (unsigned long ldim)
 Print grid.
 

Protected Attributes

size_t dim
 Number of dimensions.
 
unsigned int bins
 Number of bins.
 
unsigned int boxes
 Number of boxes.
 
ubvector xi
 Boundaries for each bin.
 
ubvector xin
 Storage for grid refinement.
 
ubvector delx
 The iteration length in each direction.
 
ubvector weight
 The weight for each bin.
 
double vol
 The volume of the current bin.
 
ubvector_int bin
 The bins for each direction.
 
ubvector_int box
 The boxes for each direction.
 
ubvector d
 Distribution.
 
unsigned int it_start
 The starting iteration number.
 
unsigned int it_num
 The total number of iterations.
 
unsigned int samples
 Number of samples for computing chi squared.
 
unsigned int calls_per_box
 Number of function calls per box.
 
vec_t x
 Point for function evaluation.
 
Scratch variables preserved between calls to
double jac
 
double wtd_int_sum
 
double sum_wgts
 
double chi_sum
 
- Protected Attributes inherited from o2scl::inte_multi< func_t, vec_t >
double interror
 The uncertainty for the last integration computation.
 

Static Protected Attributes

static const size_t bins_max =50
 Maximum number of bins.
 

Integration mode (default is mode_importance)

int mode
 
static const int mode_importance =1
 
static const int mode_importance_only =0
 
static const int mode_stratified =-1
 

Detailed Description

template<class func_t = multi_funct11, class vec_t = boost::numeric::ublas::vector<double>, class rng_t = rng_gsl>
class o2scl::mcarlo_vegas< func_t, vec_t, rng_t >

The output options are a little different than the original GSL routine. The default setting of mcarlo::verbose is 0, which turns off all output. A verbose value of 1 prints summary information about the weighted average and final result, while a value of 2 also displays the grid coordinates. A value of 3 prints information from the rebinning procedure for each iteration.

Some original documentation from GSL:

The input coordinates are x[j], with upper and lower limits
xu[j] and xl[j]. The integration length in the j-th direction is
delx[j]. Each coordinate x[j] is rescaled to a variable y[j] in
the range 0 to 1. The range is divided into bins with boundaries
xi[i][j], where i=0 corresponds to y=0 and i=bins to y=1. The
grid is refined (ie, bins are adjusted) using d[i][j] which is
some variation on the squared sum. A third parameter used in
defining the real coordinate using random numbers is called z.
It ranges from 0 to bins. Its integer part gives the lower index
of the bin into which a call is to be placed, and the remainder
gives the location inside the bin.

When stratified sampling is used the bins are grouped into
boxes, and the algorithm allocates an equal number of function
calls to each box.

The variable alpha controls how "stiff" the rebinning algorithm
is. alpha = 0 means never change the grid. Alpha is typically
set between 1 and 2. 
Todo:
Mode = importance only doesn't give the same answer as GSL yet.
Idea for Future:
Prettify the verbose output
Idea for Future:
Allow the user to get information about the how the sampling was done, possibly by converting the bins and boxes into a structure or class.
Idea for Future:
Allow the user to change the maximum number of bins.

Based on Lepage78 . The current version of the algorithm was described in the Cornell preprint CLNS-80/447 of March,

  1. The GSL code follows most closely the C version by D. R. Yennie, coded in 1984.

Definition at line 115 of file mcarlo_vegas.h.

Member Function Documentation

◆ accumulate_distribution()

template<class func_t = multi_funct11, class vec_t = boost::numeric::ublas::vector<double>, class rng_t = rng_gsl>
void o2scl::mcarlo_vegas< func_t, vec_t, rng_t >::accumulate_distribution ( ubvector_int lbin,
double  y 
)
inlineprotected

This is among the member functions that is not virtual because it is part of the innermost loop.

Definition at line 282 of file mcarlo_vegas.h.

◆ change_box_coord()

template<class func_t = multi_funct11, class vec_t = boost::numeric::ublas::vector<double>, class rng_t = rng_gsl>
int o2scl::mcarlo_vegas< func_t, vec_t, rng_t >::change_box_coord ( ubvector_int boxt)
inlineprotected

Steps through the box coordinates, e.g.

{0,0}, {0,1}, {0,2}, {0,3}, {1,0}, {1,1}, {1,2}, ...

This is among the member functions that is not virtual because it is part of the innermost loop.

Definition at line 230 of file mcarlo_vegas.h.

◆ random_point()

template<class func_t = multi_funct11, class vec_t = boost::numeric::ublas::vector<double>, class rng_t = rng_gsl>
void o2scl::mcarlo_vegas< func_t, vec_t, rng_t >::random_point ( vec_t &  lx,
ubvector_int lbin,
double &  bin_vol,
const ubvector_int lbox,
const vec_t &  xl,
const vec_t &  xu 
)
inlineprotected

Use the random number generator to return a random position x in a given box. The value of bin gives the bin location of the random position (there may be several bins within a given box)

This is among the member functions that is not virtual because it is part of the innermost loop.

Definition at line 301 of file mcarlo_vegas.h.

◆ vegas_minteg_err()

template<class func_t = multi_funct11, class vec_t = boost::numeric::ublas::vector<double>, class rng_t = rng_gsl>
virtual int o2scl::mcarlo_vegas< func_t, vec_t, rng_t >::vegas_minteg_err ( int  stage,
func_t &  func,
size_t  ndim,
const vec_t &  xl,
const vec_t &  xu,
double &  res,
double &  err 
)
inlinevirtual

Original documentation from GSL:

Normally, stage = 0 which begins with a new uniform grid and empty weighted average. Calling vegas with stage = 1 retains the grid from the previous run but discards the weighted average, so that one can "tune" the grid using a relatively small number of points and then do a large run with stage = 1 on the optimized grid. Setting stage = 2 keeps the grid and the weighted average from the previous run, but may increase (or decrease) the number of histogram bins in the grid depending on the number of calls available. Choosing stage = 3 enters at the main loop, so that nothing is changed, and is equivalent to performing additional iterations in a previous call.

Todo:

Should stage be passed by reference?

There was an update between gsl-1.12 and 1.15 which has not been implemented here yet.

Definition at line 608 of file mcarlo_vegas.h.

Member Data Documentation

◆ alpha

template<class func_t = multi_funct11, class vec_t = boost::numeric::ublas::vector<double>, class rng_t = rng_gsl>
double o2scl::mcarlo_vegas< func_t, vec_t, rng_t >::alpha

This usual range is between 1 and 2.

Definition at line 141 of file mcarlo_vegas.h.

◆ chisq

template<class func_t = multi_funct11, class vec_t = boost::numeric::ublas::vector<double>, class rng_t = rng_gsl>
double o2scl::mcarlo_vegas< func_t, vec_t, rng_t >::chisq

After an integration, this should be close to 1. If it is not, then this indicates that the values of the integral from different iterations are inconsistent, and the error may be underestimated. Further iterations of the algorithm may enable one to obtain more reliable results.

Definition at line 155 of file mcarlo_vegas.h.


The documentation for this class was generated from the following file:

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