Vegas(int cdim, int nsdim, int cnbin, int ddim, int adim,
const int nvals[], int fn, Integrand fxn = 0, int ncpu = 1);
cdim: number of continuous dimensions
nsdim: number of non-separable continuous dimensions (first nsdim dimensions)
cnbin: number of bins for each continuous dimension
ddim: number of discrete dimensions
adim: number of auxilliary dimensions (no adaptation)
nvals[]: number of values for each discrete dimension (see below)
fn: number of functions to integrate, (weights are adapted for first function)
fxn: pointer to function (see next section)
ncpu: number of CPUs, default: 1 (see section 5.1)
The default integration/summation range is
Random numbers in (0, 1) are provided to facilitate the addition of auxilliary dimensions that are not VEGAS driven.
The parameters ddim or adim can be set to zero if no discrete
or auxilliary dimensions are needed.
The parameters of a constructed Vegas object can be modified with two member functions:
void set_parameters(int cdim, int cnbin, int ddim, int adim, const int nvals[], int fn); void set_integrand(Integrand fxn);
Note that these functions are questionable from a design perspective and may not be supported in the next major version, i.e. versions 2.0.0 and higher.
If you plan to use the get and put methods (see section 2.6)
directly rather than the loop method, just pass 0 for fxn.
If a function pointer fxn of type Integrand:
typedef void (*Integrand)(const double x[], const int k[], const double& wgt,
const double aux[], double f[]);
is passed to the constructor or set_integrand then the function pl_vegas
runs through iterats iterations with
ncall shots each when called, modeling the behavior of Peter Lepage's VEGAS:
void pl_vegas(Vegas& vegas, int64 ncall, int iterats, int init = 0);
The array x[] contains the random numbers for the continuous dimensions,
the array k[] contains the random numbers for the discrete dimensions,
the array aux[] contains the auxilliary random numbers. The VEGAS weight (For
each iteration the sum of the weights is normalized to 1.)
is passed in variable wgt, so that it can be passed on to functions that
fill histograms, for example. The function *fxn is expected to fill the
array f[] with the function values that correspond to the arguments
x[0], x[1], ..., k[0], k[1], ... and aux[0], aux[1], ...
The value of the first function is assigned to f[0]. The integration is optimized
for this function according to the VEGAS algorithm. Secondary functions can simultaneously be
integrated and their values are expected to be returned as f[1], f[2], etc.
The parameter init has the same meaning as in the original VEGAS program:
init = 0: fresh start (no prior information)
init = 1: use previous adaptation info, but discard shot data
init = 2: use previous adaptation info and build on previous shots data
The generic action for a single iteration is encapsulated in the member function:
void loop(int64 ncall_);
Note that ncall's type is int64, corresponding to a 64-bit integer.
int64 might be a typedef to long int or long long int depending
on your system. You can find that out by looking at the values of the
macros SIZE_LONG and SIZE_LONG_LONG.
Both, class Vegas and the associated function pl_vegas are enclosed in
namespace MonteCarlo.
fn (since the last adaptation):
double get_integral(int fn) const ; double get_error(int fn) const ; double get_redchi2(int fn) const ; // "reduced chi^2" = chi^2/(nr of estimates-1) int get_niter() const ; // number of iterations (estimates)
or accumulated over all iterations:
double get_acc_integral(int fn) const ; double get_acc_error(int fn) const ; double get_acc_redchi2(int fn) const ; int get_acc_niter() const ;
The state of a Vegas object can be saved to file and restored later with the following methods:
void save(const char* filename = "vegas.dat") const ; void save(std::ofstream& fout) const ; void restore(Integrand fxn, const char* filename = "vegas.dat"); void restore(Integrand fxn, std::ifstream& fin); void restore(std::ifstream& fin); // integrand set with set_integrand()
Saved are the parameters that can be set with set_parameters and all
weights. This does not include the function pointer of type Integrand and the
number of CPUs. The restore method does not check if the function to be integrated
is the same as when save was called.
#ifdef NO_CXX_HEADERS_FOR_C_LIB
#include <math.h>
#else
#include <cmath>
#endif
#include "dvegas.h"
#ifndef NO_CXX_HEADERS_FOR_C_LIB
// <cmath>
using std::pow ;
using std::sqrt ;
using std::exp ;
#endif
void fxn(const double x[], const int k[], const double& wgt, const double aux[], double f[])
{
f[0] = x[0]*x[1]*x[2]*x[3]*x[4]*pow(2.0,5);
}
void fxn1(const double x[], const int k[], const double& wgt, const double aux[], double f[])
{
if (k[0] < 99)
f[0] = 1 * x[0]*x[1]*x[2]*x[3]*x[4]*pow(2.0,5);
else if (k[0] == 99)
f[0] = 100 * x[0]*x[1]*x[2]*x[3]*x[4]*pow(2.0,5);
}
int main()
{
using MonteCarlo::Vegas ;
using MonteCarlo::pl_vegas ;
int nvals[2];
Vegas vegas(5, 50, 0, 0, nvals, 1, &fxn);
pl_vegas(vegas, 100000, 2); // result: 1.0
nvals[0] = 100 ;
vegas.set_parameters(5, 50, 1, 0, nvals, 1);
vegas.set_integrand(&fxn1);
pl_vegas(vegas, 100000, 1); // result: 199.0, factor sqrt(100)=10 error reduction with discrete weights
vegas.save();
pl_vegas(vegas, 100000, 1, 1);
nvals[1] = 3 ;
vegas.set_parameters(5, 50, 2, 0, nvals, 1);
pl_vegas(vegas, 100000, 2); // result: 597.0
vegas.restore(&fxn1);
vegas.save("same.dat");
pl_vegas(vegas, 100000, 1, 1); // result: 199.0
}
#if SIZEOF_LONG == 8 typedef long int int64 ; #elif SIZEOF_LONG_LONG == 8 typedef long long int int64 ; #endif typedef double float64 ;
The type int64 is a 64-bit integer and is used to allow a maximum number of
9,223,372,036,854,755,807. (The typical 32-bit integer has a maximum of about
2 billion, which I found to be too small in some cases to achieve the desired
precision.) The type float64 is a 64-bit floating point number and is used
in some member functions for variables that accumulate a large number of
temporary results (of the order of ncall). 32-bit floating point numbers provide
only a precision of typically 7 decimal digits, which can be insufficient when adding
about 100 million numbers or more. For 64-bit floating point numbers with a
precision of typically 15 decimal digits this problem occurs roughly at 10^16
additions. On most systems the type double corresponds to 64-bit floating point
numbers. In rare cases float64 might correspond to long double.
(Note that these problems can also occur in histogram packages.)
In a future dVegas release the issues addressed above will be
double-checked in the pre_loop member function,
and a warning or error will be issued when problems could arise. This kind
of checking is facilitated in C++ by specializations of the numeric_limits
template presented in <limits>. See section 22.2, Numeric Limits, in
Bjarne Stroustrup's The C++ programming language, 3rd ed.
Sampling specifies sampling techniques:
enum Sampling { importance, stratified } ;
Currently only importance sampling is implemented, since I use dVegas for calculations with high numbers of dimensions. Additional stratified sampling can lead to a substantially improved efficiency in low dimensions (cdim = 1, ..., 4) and will be implemented if the need arises.
The enumerations Reset and Info hold sets of options. They allow
to specify elementary or composite options by name rather than a crypic
(hexadecimal) number. The state of each elementary option is stored in a bit
that is part of a bit set of related options.
The enumeration Reset specifies reset options:
enum Reset { r_data=0x1, r_weights=0x2, r_none=0x0, r_all=0x3 };
The enumeration Info holds a set of options for the info member function,
which prints out information about the most recent iteration, as well as accumulated
results and statistics. The info member function is a primary target for
customization and the options listed below are intended to give an idea of the possibilities:
enum Info { i_spec =0x0001, i_int =0x0002, i_adapt =0x0004, i_grid =0x0008,
i_spec_more=0x0010, i_int_more=0x0020, i_adapt_more=0x0040, i_grid_more=0x0080,
i_spec_all =0x0011, i_int_all =0x0022, i_adapt_all =0x0044, i_grid_all =0x0088,
i_all_short=0x000f, i_all_long=0x00ff };
Here, spec relates to information like parameters and settings, int to information
about the integration itself like results, errors and statistics, adapt to information
about the progression and efficiency of the adaptation cycles, and grid options would
print out weight data (grid data for continuous dimensions). In each category the amount
of information displayed can be regulated, e.g. by appending _more and _all.
In this release the Info argument of the info member function is just a
dummy argument and the function outputs information in a style that I found useful.
Based on feedback from other users a variety of output options will be implemented in
future dVegas releases.
class Vegas
{
public:
void loop(int64 ncall_); // accumulate data
void* loop_worker() const ; // thread function for parallel mode
void adapt(); // adapt weights
void reset(Reset flag); // discard data
void info(Info flag) const ; // print info
void get(double x[], int k[], double* wgt_, double aux[], int icpu = 0);
void put(const double f[], int icpu = 0);
void pre_loop(int64 ncall_);
void post_loop();
Vegas(int cdim, int cnbin, int ddim, int adim, const int nvals[], int fn_, Integrand fxn = 0, int ncpu_ = 1);
virtual ~Vegas() { }
// the following two functions are questionable from a design perspective and may not
// be available in the next major version (i.e. versions 2.0.0 and higher)
void set_parameters(int cdim, int cnbin, int ddim, int adim, const int nvals[], int fn_);
void set_integrand(Integrand fxn);
// results since last adaptation
double get_integral(int fn_) const ;
double get_error(int fn_) const ;
double get_redchi2(int fn_) const ; // "reduced chi^2" = chi^2/(nr of estimates-1)
int get_niter() const ; // nr of estimates/iterations
// accumulated results
double get_acc_integral(int fn_) const ;
double get_acc_error(int fn_) const ;
double get_acc_redchi2(int fn_) const ; // "reduced chi^2" = chi^2/(nr of estimates-1)
int get_acc_niter() const ; // nr of estimates/iterations
// save and read weights information and configuration (possibly multiple times)
void save(const char* filename = "vegas.dat") const ;
void save(std::ofstream& fout) const ;
void restore(Integrand fxn, const char* filename = "vegas.dat");
void restore(Integrand fxn, std::ifstream& fin);
void restore(std::ifstream& fin); // integrand set with set_integrand()
friend std::ostream& operator<<(std::ostream& s, const Vegas& v);
friend std::istream& operator>>(std::istream& s, Vegas& v);
// digital watermark with version information
void print_version() const ;
// use with caution
void set_niter(int i);
private:
// ...
};