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; }; }
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.
f
integrandstypedef void Integrand(const double x[], const int k[], const double& weight, const double aux[], double f[]);
cDim
random floats in (0, 1)
and thus covers
the integration volume for the adapted continuous dimensions
dDimSizes.size()
random integers in [0, dDimSizes[j]
) with j = 0, ..., dDimSizes.size() - 1
aDim
random numbers in (0, 1)
that have no effect on the adaptation
f
floats obtained by evaluating each integrand
at the point specified by the input parameters f[0]
) drives the adaptation. 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.
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.)
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 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).
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:
A comprehensive set of examples can be found in the file dvegas_demo.cpp
.