Chapter 2: Dvegas

2.1: Dvegas interface

namespace HepSource
{
class Dvegas
{
public:
  Dvegas(const int aDim, const int f, Integrand *const i);
  Dvegas(const int cDim, const int cBin, const int f, Integrand *const i,
         const Sampling sampling = IMPORTANCE);
  Dvegas(const int cDim, const int cBin, const vector<int>& dDimSizes, const int f,
         Integrand *const i, const Sampling sampling = IMPORTANCE);
  Dvegas(const int cDim, const int cBin, const int corrDim,
         const vector<int>& dDimSizes, const int aDim, const int f, Integrand *const i,
         const Sampling sampling = IMPORTANCE);
  Dvegas(const int cDim, const int cBin, const vector<int>& corrDim,
         const vector<int>& dDimSizes, const int aDim, const int f, Integrand *const i,
         const Sampling sampling = IMPORTANCE);
  Dvegas();
  virtual ~Dvegas();

  CellAccumulators collectData(const Int64 numberOfShots);

  void info() const;     // print info about best estimate
  const bool isAdaptingContinuousDimensions() const;
  const double chiSquarePerIteration() const;

  void adapt(CellAccumulators& cellAcc);     // adapt weights
  const double getCurvature() const;
  void setCurvature(const double curvature);
  const double getRootAmpl() const;
  void setRootAmpl(const double rootAmpl);

  void disableOptimize(const int setOfCorrelatedWeights);
  
  void reset();          // discard all data
  void resetWeights();
  void resetEstimates();

  void saveStateToFile(const string& filename = "dvegas.state") const;
  void restoreStateFromFile(Integrand *const i, const string& filename = "dvegas.state");
  void restoreWeightsFromFile(Integrand *const i, const string& filename = "dvegas.state");
  void saveStateToStream(ostream& os) const;
  void restoreStateFromStream(Integrand *const i, istream& is);
  void restoreWeightsFromStream(Integrand *const i, istream& is);
  const string saveStateToString() const;
  void restoreStateFromString(Integrand *const i, const string& s);
  void restoreWeightsFromString(Integrand *const i, const string& s);

  void activateCorrelationsDetector();

  void attachOmniComp(OmniComp *const omniComp);
  void detachOmniComp();

  void printVersion() const;
};
}

2.2: Dvegas constructors

2.2.1: Constructor arguments

cDim: int: number of adapted continuous dimensions
cBin: int: resolution of adaptation (number of bins)
corrDim: vector<int>: defines sets of continuous dimensions for which the adaptation fully takes correlations into account
Example: cDim = 8 and corrDim = {3, 5}: set1 = {0, 1, 2}, set2 = {3, 4}, correlations not sampled for 5, 6, 7
dDimSizes: vector<int>: defines discrete dimensions
Example: dDimSizes = {3, 5}: two discrete dimensions with index1 in {0, 1, 2} and index2 in {0, 1, 2, 3, 4}
aDim: int: number of "auxilliary" (i.e. non- adapted) continuous dimensions
f: int: number of integrands to be evaluated (first integrand drives adapation, see below)
i: Integrand*: pointer to function that evaluates the integrand(s)
sampling: enum Sampling {NONE, IMPORTANCE, STRATIFIED}: Importance sampling concentrates sampled points in hypercubes that contribute most to the integral.
Stratified sampling concentrates sampled points in hypercubes that contribute most to the error of the integral.

The integration range for adapted and auxilliary continuous dimensions is (0, 1).

Note that, except for i, all arguments can be zero (or empty vector), thus disabling the particular functionality.

2.2.2: Function that evaluates `f` integrands

typedef void Integrand(const double x[], const int k[], const double& weight,
                       const double aux[], double f[]);

x: double array (in): contains cDim random floats in (0, 1) and thus covers the integration volume for the adapted continuous dimensions
k: int array (in): contains dDimSizes.size() random integers in [0, dDimSizes[j]) with j = 0, ..., dDimSizes.size() - 1
covering all possible index combinations for the discrete dimensions
weight: double (in): weight of the sampled point, which may, for example, be used to fill histograms
Note that by default the weight is normalized in each iteration to yield the average integrand value--corresponding to the value of the integral (assuming unit volume)--rather than the sum.
aux: double array (in): contains aDim random numbers in (0, 1) that have no effect on the adaptation
f: double array (out): passes back f floats obtained by evaluating each integrand at the point specified by the input parameters
Note that the first integrand (i.e. values passed back through f[0]) drives the adaptation.
Depending on how similar the other integrands are to the first the adapted sampling may or may not be well suited for them.

2.2.3: Constructors

Dvegas(const int aDim, const int f, Integrand *const i);

"Crude", unadapted Monte Carlo sampling.

Dvegas(const int cDim, const int cBin, const int f, Integrand *const i,
       const Sampling sampling = IMPORTANCE);

Classic VEGAS case: the sampling is adapted for a number of continuous dimensions that are assumed to be independent.

Dvegas(const int cDim, const int cBin, const vector<int>& dDimSizes, const int f,
       Integrand *const i, const Sampling sampling = IMPORTANCE);

Full adaptation for sums of integrals.

Dvegas(const int cDim, const int cBin, const int corrDim,
       const vector<int>& dDimSizes, const int aDim, const int f, Integrand *const i,
       const Sampling sampling = IMPORTANCE);

Exclude "auxilliary" dimensions from adaptation.

Dvegas(const int cDim, const int cBin, const vector<int>& corrDim,
       const vector<int>& dDimSizes, const int aDim, const int f, Integrand *const i,
       const Sampling sampling = IMPORTANCE);

Relax independence assumption for continuous dimensions and fully sample correlations for selected dimensions.

Dvegas();

Use default constructor to subsequently initialize with restoreStateFromStream() or a similar method.

2.3: 64-bit numeric types

The following type aliases are defined:

#if (defined LONG_IS_INT64) && (!defined LONG_LONG_IS_INT64)
typedef long int Int64;
#elif (defined LONG_LONG_IS_INT64) && (!defined LONG_IS_INT64)
typedef long long int Int64;
#endif
typedef double Float64;

The type alias Int64 corresponds to 64-bit integers, which allow for a maximum of 2^63 - 1 = 9,223,372,036,854,755,807 sampled points. Note that 32-bit integers only allow up to about 2 billion sampled points (2^31 - 1 = 2,147,483,647 to be precise), which in some cases may not suffice to achieve the desired precision. On 32-bit processor hardware including most Intel and AMD processors (i386 architecture), the macro LONG_LONG_IS_INT64 should be defined. On 64-bit processor hardware, on the other hand, like DEC's Alpha, Intel's Itanium or AMD's Sledgehammer processor, the macro LONG_IS_INT64 should be defined (see Makefile).

The type alias Float64 corresponds to 64-bit floating point numbers and is used internally for variables that accumulate a large number of temporary results (of O(numberOfShots)). 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 platforms the type double corresponds to 64-bit floating point numbers. In rare cases Float64 might correspond to long double. In this case the typedef should be adjusted accordingly. (Note that this issue is also relevant for histogram packages.)

2.4: Dvegas methods

CellAccumulators collectData(const Int64 numberOfShots);
void adapt(CellAccumulators& cellAcc);

Method collectData() samples numberOfShots points in the integration volume, i.e. advances the Monte Carlo integration by one iteration and returns data accumulated for each cell in an object of type CellAccumulators. This object is subsequently passed to method adapt(), which adapts all weights, and thus improves the sampling in the next iteration.

void info() const;

At any time method info() can be used to print out information regarding the currently best estimate for the integral(s).

const double getCurvature() const;
void setCurvature(const double curvature);
const double getRootAmpl() const;
void setRootAmpl(const double rootAmpl);

The parameters curvature and rootAmpl ("root amplification") control the promotion of smaller weights, which damps the adaptation in order to avoid oscillations. See implementation of helper function promoteSmallerWeights() in dvegas.cpp for details.

void disableOptimize(const int setOfCorrelatedWeights);

If correlations are sampled, an optimization is performed before the sampling occurs. In rare cases the optimization step requires more time than is subsequently saved. In such cases the optimization can be turned off with the disableOptimize() method.

void reset();
void resetWeights();
void resetEstimates();

reset() discards all state-related internal data, reverting the object back to the state immediately after construction.

void saveStateToFile(const string& filename = "dvegas.state") const;
void restoreStateFromFile(Integrand *const i, const string& filename = "dvegas.state");
void restoreWeightsFromFile(Integrand *const i, const string& filename = "dvegas.state");
void saveStateToStream(ostream& os) const;
void restoreStateFromStream(Integrand *const i, istream& is);
void restoreWeightsFromStream(Integrand *const i, istream& is);
const string saveStateToString() const;
void restoreStateFromString(Integrand *const i, const string& s);
void restoreWeightsFromString(Integrand *const i, const string& s);

Methods are provided to save Dvegas objects to a stream, file or string and subsequently restore the saved state
(or the most recent weights without the integral estimates of previous iterations).

void activateCorrelationsDetector();

By calling activateCorrelationsDetector() a separate module is activated that, during the following iteration, probes (for the first integrand) to what degree combinations of two dimensions are independent (as assumed by the VEGAS algorithm). If strong correlations are detected, they can then be fully sampled by setting corrDim accordingly. Circular correlations should be properly combined, e.g. the correlations {0, 1} and {1, 2} translate to the set {0, 1, 2}.

void attachOmniComp(OmniComp *const omniComp);
void detachOmniComp();

Before distributed integrand evaluation can be used to accelerate computations, an OmniComp object needs to be attached to the Dvegas object with attachOmniComp (see section 3 for more on OmniComp).

const bool isAdaptingContinuousDimensions() const;
const double chiSquarePerIteration() const;

These getter methods are used in VEGAS() for diagnostic warnings (see section 2.5).

2.5: VEGAS function

namespace HepSource
{
void VEGAS(Dvegas& dvegas, const Int64 numberOfShots, const int numberOfIterations,
           const int init = 0);
}

This function mimics the interface of the original VEGAS program. The parameter init has the following meaning:

init = 0: fresh start--discard all weight and estimate information
init = 1: use current weights, but discard prior estimates
init = 2: use current weights and take into account estimates of previous iterations

2.6: Examples

A comprehensive set of examples can be found in the file dvegas_demo.cpp.