`vegas` Module¶

Introduction¶

The key Python objects supported by the vegas module are:

vegas.Integrator — an object describing a multidimensional integration operator. Such objects contain information about the integration volume, and also about optimal remappings of the integration variables based upon the last integral evaluated using the object.

vegas.AdaptiveMap — an object describing the remappings used by vegas.

vegas.RAvg — an object describing the result of a vegas integration. vegas returns the weighted average of the integral estimates from each vegas iteration as an object of class vegas.RAvg. These are Gaussian random variables — that is, they have a mean and a standard deviation — but also contain information about the iterations vegas used in generating the result.

vegas.RAvgArray — an array version of vegas.RAvg used when the integrand is array-valued.

vegas.RAvgDict — a dictionary version of vegas.RAvg used when the integrand is dictionary-valued.

vegas.PDFIntegrator — a specialized integrator for evaluating Gaussian expectation values.

These are described in detail below.

Integrator Objects¶

The central component of the vegas package is the integrator class:

class vegas.Integrator¶

Adaptive multidimensional Monte Carlo integration.

vegas.Integrator objects make Monte Carlo estimates of multidimensional functions f(x) where x[d] is a point in the integration volume:

integ = vegas.Integrator(integration_region)

result = integ(f, nitn=10, neval=10000)

The integator makes nitn estimates of the integral, each using at most neval samples of the integrand, as it adapts to the specific features of the integrand. Successive estimates (iterations) typically improve in accuracy until the integrator has fully adapted. The integrator returns the weighted average of all nitn estimates, together with an estimate of the statistical (Monte Carlo) uncertainty in that estimate of the integral. The result is an object of type RAvg (which is derived from gvar.GVar).

Integrands f(x) return numbers, arrays of numbers (any shape), or dictionaries whose values are numbers or arrays (any shape). Each number returned by an integrand corresponds to a different integrand. When arrays are returned, vegas adapts to the first number in the flattened array. When dictionaries are returned, vegas adapts to the first number in the value corresponding to the first key.

vegas can generate integration points in batches for integrands built from classes derived from vegas.BatchIntegrand, or integrand functions decorated by vegas.batchintegrand(). Batch integrands are typically much faster, especially if they are coded in Cython or C/C++ or Fortran.

vegas.Integrators have a large number of parameters but the only ones that most people will care about are: the number nitn of iterations of the vegas algorithm; the maximum number neval of integrand evaluations per iteration; and the damping parameter alpha, which is used to slow down the adaptive algorithms when they would otherwise be unstable (e.g., with very peaky integrands). Setting parameter analyzer=vegas.reporter() is sometimes useful, as well, since it causes vegas to print (on sys.stdout) intermediate results from each iteration, as they are produced. This helps when each iteration takes a long time to complete (e.g., longer than an hour) because it allows you to monitor progress as it is being made (or not).

Parameters:

map (array, vegas.AdaptiveMap or vegas.Integrator) –
The integration region as specified by an array map[d, i] where d is the direction and i=0,1 specify the lower and upper limits of integration in direction d.

map could also be the integration map from another vegas.Integrator, or that vegas.Integrator itself. In this case the grid is copied from the existing integrator.
nitn (positive int) – The maximum number of iterations used to adapt to the integrand and estimate its value. The default value is 10; typical values range from 10 to 20.
neval (positive int) – Approximate number of integrand evaluations in each iteration of the vegas algorithm. Increasing neval increases the precision: statistical errors should fall at least as fast as sqrt(1./neval) and often fall much faster. The default value is 1000; real problems often require 10–10,000 times more evaluations than this.
nstrat (int array) – nstrat[d] specifies the number of stratifications to use in direction d. By default this parameter is set automatically, based on parameter neval, with nstrat[d] approximately the same for every d. Specifying nstrat explicitly makes it possible to concentrate stratifications in directions where they are most needed. If nstrat is set but neval is not, neval is set equal to 2*prod(nstrat)/(1-neval_frac).
alpha (float) – Damping parameter controlling the remapping of the integration variables as vegas adapts to the integrand. Smaller values slow adaptation, which may be desirable for difficult integrands. Small or zero alphas are also sometimes useful after the grid has adapted, to minimize fluctuations away from the optimal grid. The default value is 0.5.
beta (float) – Damping parameter controlling the redistribution of integrand evaluations across hypercubes in the stratified sampling of the integral (over transformed variables). Smaller values limit the amount of redistribution. The theoretically optimal value is 1; setting beta=0 prevents any redistribution of evaluations. The default value is 0.75.
neval_frac (float) – Approximate fraction of function evaluations used for adaptive stratified sampling. vegas distributes (1-neval_frac)*neval integrand evaluations uniformly over all hypercubes, with at least 2 evaluations per hypercube. The remaining neval_frac*neval evaluations are concentrated in hypercubes where the errors are largest. Increasing neval_frac makes more integrand evaluations available for adaptive stratified sampling, but reduces the number of hypercubes, which limits the algorithm’s ability to adapt. Ignored when beta=0. Default is neval_frac=0.75.
adapt (bool) – Setting adapt=False prevents further adaptation by vegas. Typically this would be done after training the vegas.Integrator on an integrand, in order to stabilize further estimates of the integral. vegas uses unweighted averages to combine results from different iterations when adapt=False. The default setting is adapt=True.
nproc (int or None) –
Number of processes/processors used to evalute the integrand. If nproc>1 Python’s multiprocessing module is used to spread integration points across multiple processors, thereby potentially reducing the time required to evaluate the integral. There is a significant overhead involved in using multiple processors so this option is useful only when the integrand is expensive to evaluate. When using the multiprocessing module in its default mode for MacOS and Windows, it is important that the main module can be safely imported (i.e., without launching new processes). This can be accomplished with some version of the if __name__ == '__main__`: construct in the main module: e.g.,
```
if __name__ == '__main__':
    main()
```
This is not an issue on other Unix platforms. See the multiprocessing documentation for more information. Note that setting nproc greater than 1 disables MPI support. Set nproc=None to use all the processors on the machine (equivalent to nproc=os.cpu_count()). Default value is nproc=1. (Requires Python 3.3 or later.)
analyzer –
An object with methods

analyzer.begin(itn, integrator)
analyzer.end(itn_result, result)

where: begin(itn, integrator) is called at the start of each vegas iteration with itn equal to the iteration number and integrator equal to the integrator itself; and end(itn_result, result) is called at the end of each iteration with itn_result equal to the result for that iteration and result equal to the cummulative result of all iterations so far. Setting analyzer=vegas.reporter(), for example, causes vegas to print out a running report of its results as they are produced. The default is analyzer=None.
uses_jac (bool) – Setting uses_jac=True causes vegas to call the integrand with two arguments: fcn(x, jac=jac). The second argument is the Jacobian jac[d] = dx[d]/dy[d] of the vegas map. The integral over x[d] of 1/jac[d] equals 1 (exactly). The default setting is uses_jac=False.
nhcube_batch (positive int) – The number of hypercubes (in y space) whose integration points are combined into a single batch to be passed to the integrand, together, when using vegas in batch mode. The default value is 10,000. Larger values may be lead to faster evaluations, but at the cost of more memory for internal work arrays.
maxinc_axis (positive int) – The maximum number of increments per axis allowed for the x-space grid. The default value is 1000; there is probably little need to use other values.
max_nhcube (positive int) – Maximum number of y-space hypercubes used for stratified sampling. Setting max_nhcube=1 turns stratified sampling off, which is probably never a good idea. The default setting (1e9) was chosen to correspond to the point where internal work arrays become comparable in size to the typical amount of RAM available to a processor (in a laptop in 2014). Internal memory usage is large only when beta>0 and minimize_mem=False; therefore max_nhcube is ignored if beta=0 or minimize_mem=True.
max_neval_hcube (positive int) – Maximum number of integrand evaluations per hypercube in the stratification. The default value is 1e6. Larger values might allow for more adaptation (when beta>0), but also can result in large internal work arrays.
max_mem (positive float) – Maximum number of floats allowed in internal work arrays (approx.). A MemoryError is raised if the work arrays are too large, in which case one might want to reduce max_neval_hcube or max_nhcube (or increase max_mem if there is enough RAM). Default value is 1e9.
rtol (float) – Relative error in the integral estimate at which point the integrator can stop. The default value is 0.0 which turns off this stopping condition. This stopping condition can be quite unreliable in early iterations, before vegas has converged. Use with caution, if at all.
atol (float) – Absolute error in the integral estimate at which point the integrator can stop. The default value is 0.0 which turns off this stopping condition. This stopping condition can be quite unreliable in early iterations, before vegas has converged. Use with caution, if at all.
ran_array_generator – Function that generates numpy arrays of random numbers distributed uniformly between 0 and 1. ran_array_generator(shape) should create an array whose dimensions are specified by the integer-valued tuple shape. The default generator is numpy.random.random.
sync_ran (bool) – If True (default), the default random number generator is synchronized across all processors when using MPI. If False, vegas does no synchronization (but the random numbers should synchronized some other way).
adapt_to_errors (bool) –
adapt_to_errors=False causes vegas to remap the integration variables to emphasize regions where |f(x)| is largest. This is the default mode.

adapt_to_errors=True causes vegas to remap variables to emphasize regions where the Monte Carlo error is largest. This might be superior when the number of the number of stratifications (self.nstrat) in the y grid is large (> 100). It is typically useful only in one or two dimensions.
uniform_nstrat (bool) – If True, requires that the nstrat[d] be equal for all d. If False (default), the algorithm maximizes the number of stratifications while requiring |nstrat[d1] - nstrat[d2]| <= 1. This parameter is ignored if nstrat is specified explicitly.
mpi (bool) – Setting mpi=False disables mpi support in vegas even if mpi is available; setting mpi=True (default) allows use of mpi provided module mpi4py is installed. This flag is ignored when mpi is not being used (and so has no impact on performance).
minimize_mem – When True, vegas minimizes internal workspace at the cost of extra evaluations of the integrand. This can increase execution time by 50–100% but might be desirable when the number of evaluations is very large (e.g., neval=1e9). Normally vegas uses internal work space that grows in proportion to neval. If that work space exceeds the size of the RAM available to the processor, vegas runs much more slowly. Setting minimize_mem=True greatly reduces the internal storage used by vegas; in particular memory becomes independent of neval. The default setting (minimize_mem=False), however, is much superior unless memory becomes a problem. (The large memory is needed for adaptive stratified sampling, so memory is not an issue if beta=0.)

Methods:

__call__(fcn, save=None, saveall=None, **kargs)¶

Integrate integrand fcn.

A typical integrand has the form, for example:

def f(x):
    return x[0] ** 2 + x[1] ** 4

The argument x[d] is an integration point, where index d=0... represents direction within the integration volume.

Integrands can be array-valued, representing multiple integrands: e.g.,

def f(x):
    return [x[0] ** 2, x[0] / x[1]]

The return arrays can have any shape. Dictionary-valued integrands are also supported: e.g.,

def f(x):
    return dict(a=x[0] ** 2, b=[x[0] / x[1], x[1] / x[0]])

Integrand functions that return arrays or dictionaries are useful for multiple integrands that are closely related, and can lead to substantial reductions in the errors for ratios or differences of the results.

It is usually much faster to use vegas in batch mode, where integration points are presented to the integrand in batches. A simple batch integrand might be, for example:

@vegas.batchintegrand
def f(x):
    return x[:, 0] ** 2 + x[:, 1] ** 4

where decorator @vegas.batchintegrand tells vegas that the integrand processes integration points in batches. The array x[i, d] represents a collection of different integration points labeled by i=0.... (The number of points is controlled vegas.Integrator parameter nhcube_batch.) The batch index is always first.

Batch integrands can also be constructed from classes derived from vegas.BatchIntegrand.

Batch mode is particularly useful (and fast) when the class derived from vegas.BatchIntegrand is coded in Cython. Then loops over the integration points can be coded explicitly, avoiding the need to use numpy’s whole-array operators if they are not well suited to the integrand.

Any vegas parameter can also be reset: e.g., self(fcn, nitn=20, neval=1e6).

Parameters:

fcn (callable) – Integrand function.
save (str or file or None) –
Writes results into pickle file specified by save at the end of each iteration. For example, setting save='results.pkl' means that the results returned by the last vegas iteration can be reconstructed later using:
```
import pickle
with open('results.pkl', 'rb') as ifile:
    results = pickle.load(ifile)
```
Ignored if save=None (default).
saveall (str or file or None) –
Writes (results, integrator) into pickle file specified by saveall at the end of each iteration. For example, setting saveall='allresults.pkl' means that the results returned by the last vegas iteration, together with a clone of the (adapted) integrator, can be reconstructed later using:
```
import pickle
with open('allresults.pkl', 'rb') as ifile:
    results, integrator = pickle.load(ifile)
```
Ignored if saveall=None (default).

set(ka={}, **kargs)¶

Reset default parameters in integrator.

Usage is analogous to the constructor for vegas.Integrator: for example,

old_defaults = integ.set(neval=1e6, nitn=20)

resets the default values for neval and nitn in vegas.Integrator integ. A dictionary, here old_defaults, is returned. It can be used to restore the old defaults using, for example:

integ.set(old_defaults)

settings(ngrid=0)¶

Assemble summary of integrator settings into string.

Parameters:	ngrid (int) – Number of grid nodes in each direction to include in summary. The default is 0.
Returns:	String containing the settings.

random(yield_hcube=False, yield_y=False)¶

Iterator over integration points and weights.

This method creates an iterator that returns integration points from vegas, and their corresponding weights in an integral. Each point x[d] is accompanied by the weight assigned to that point by vegas when estimating an integral. Optionally it will also return the index of the hypercube containing the integration point and/or the y-space coordinates:

integ.random()  yields  x, wgt

integ.random(yield_hcube=True) yields x, wgt, hcube

integ.random(yield_y=True) yields x, y, wgt

integ.random(yield_hcube=True, yield_y=True) yields x, y, wgt, hcube

The number of integration points returned by the iterator corresponds to a single iteration.

random_batch(yield_hcube=False, yield_y=False)¶

Iterator over integration points and weights.

This method creates an iterator that returns integration points from vegas, and their corresponding weights in an integral. The points are provided in arrays x[i, d] where i=0... labels the integration points in a batch and d=0... labels direction. The corresponding weights assigned by vegas to each point are provided in an array wgt[i].

Optionally the integrator will also return the indices of the hypercubes containing the integration points and/or the y-space coordinates of those points:

integ.random_batch()  yields  x, wgt

integ.random_batch(yield_hcube=True) yields x, wgt, hcube

integ.random_batch(yield_y=True) yields x, y, wgt

integ.random_batch(yield_hcube=True, yield_y=True) yields x, y, wgt, hcube

The number of integration points returned by the iterator corresponds to a single iteration. The number in a batch is controlled by parameter nhcube_batch.

AdaptiveMap Objects¶

vegas’s remapping of the integration variables is handled by a vegas.AdaptiveMap object, which maps the original integration variables x into new variables y in a unit hypercube. Each direction has its own map specified by a grid in x space:

$x_0 &= a \\ x_1 &= x_0 + \Delta x_0 \\ x_2 &= x_1 + \Delta x_1 \\ \cdots \\ x_N &= x_{N-1} + \Delta x_{N-1} = b$

where $a$ and $b$ are the limits of integration. The grid specifies the transformation function at the points $y=i/N$ for $i=0,1\ldots N$ :

$x(y\!=\!i/N) = x_i$

Linear interpolation is used between those points. The Jacobian for this transformation is:

$J(y) = J_i = N \Delta x_i$

vegas adjusts the increments sizes to optimize its Monte Carlo estimates of the integral. This involves training the grid. To illustrate how this is done with vegas.AdaptiveMaps consider a simple two dimensional integral over a unit hypercube with integrand:

def f(x):
   return x[0] * x[1] ** 2

We want to create a grid that optimizes uniform Monte Carlo estimates of the integral in y space. We do this by sampling the integrand at a large number ny of random points y[j, d], where j=0...ny-1 and d=0,1, uniformly distributed throughout the integration volume in y space. These samples be used to train the grid using the following code:

import vegas
import numpy as np

def f(x):
   return x[0] * x[1] ** 2

m = vegas.AdaptiveMap([[0, 1], [0, 1]], ninc=5)

ny = 1000
y = np.random.uniform(0., 1., (ny, 2))  # 1000 random y's

x = np.empty(y.shape, float)            # work space
jac = np.empty(y.shape[0], float)
f2 = np.empty(y.shape[0], float)

print('intial grid:')
print(m.settings())

for itn in range(5):                    # 5 iterations to adapt
   m.map(y, x, jac)                     # compute x's and jac

   for j in range(ny):                  # compute training data
      f2[j] = (jac[j] * f(x[j])) ** 2

   m.add_training_data(y, f2)           # adapt
   m.adapt(alpha=1.5)

   print('iteration %d:' % itn)
   print(m.settings())

In each of the 5 iterations, the vegas.AdaptiveMap adjusts the map, making increments smaller where f2 is larger and larger where f2 is smaller. The map converges after only 2 or 3 iterations, as is clear from the output:

initial grid:
    grid[ 0] = [ 0.   0.2  0.4  0.6  0.8  1. ]
    grid[ 1] = [ 0.   0.2  0.4  0.6  0.8  1. ]

iteration 0:
    grid[ 0] = [ 0.     0.412  0.62   0.76   0.883  1.   ]
    grid[ 1] = [ 0.     0.506  0.691  0.821  0.91   1.   ]

iteration 1:
    grid[ 0] = [ 0.     0.428  0.63   0.772  0.893  1.   ]
    grid[ 1] = [ 0.     0.531  0.713  0.832  0.921  1.   ]

iteration 2:
    grid[ 0] = [ 0.     0.433  0.63   0.772  0.894  1.   ]
    grid[ 1] = [ 0.     0.533  0.714  0.831  0.922  1.   ]

iteration 3:
    grid[ 0] = [ 0.     0.435  0.631  0.772  0.894  1.   ]
    grid[ 1] = [ 0.     0.533  0.715  0.831  0.923  1.   ]

iteration 4:
    grid[ 0] = [ 0.     0.436  0.631  0.772  0.895  1.   ]
    grid[ 1] = [ 0.     0.533  0.715  0.831  0.924  1.   ]

The grid increments along direction 0 shrink at larger values x[0], varying as 1/x[0]. Along direction 1 the increments shrink more quickly varying like 1/x[1]**2.

vegas samples the integrand in order to estimate the integral. It uses those same samples to train its vegas.AdaptiveMap in this fashion, for use in subsequent iterations of the algorithm.

class vegas.AdaptiveMap¶

Adaptive map y->x(y) for multidimensional y and x.

An AdaptiveMap defines a multidimensional map y -> x(y) from the unit hypercube, with 0 <= y[d] <= 1, to an arbitrary hypercube in x space. Each direction is mapped independently with a Jacobian that is tunable (i.e., “adaptive”).

The map is specified by a grid in x-space that, by definition, maps into a uniformly spaced grid in y-space. The nodes of the grid are specified by grid[d, i] where d is the direction (d=0,1...dim-1) and i labels the grid point (i=0,1...N). The mapping for a specific point y into x space is:

y[d] -> x[d] = grid[d, i(y[d])] + inc[d, i(y[d])] * delta(y[d])

where i(y)=floor(y*N), delta(y)=y*N - i(y), and inc[d, i] = grid[d, i+1] - grid[d, i]. The Jacobian for this map,

dx[d]/dy[d] = inc[d, i(y[d])] * N,

is piece-wise constant and proportional to the x-space grid spacing. Each increment in the x-space grid maps into an increment of size 1/N in the corresponding y space. So regions in x space where inc[d, i] is small are stretched out in y space, while larger increments are compressed.

The x grid for an AdaptiveMap can be specified explicitly when the map is created: for example,

m = AdaptiveMap([[0, 0.1, 1], [-1, 0, 1]])

creates a two-dimensional map where the x[0] interval (0,0.1) and (0.1,1) map into the y[0] intervals (0,0.5) and (0.5,1) respectively, while x[1] intervals (-1,0) and (0,1) map into y[1] intervals (0,0.5) and (0.5,1).

More typically, an uniform map with ninc increments is first created: for example,

m = AdaptiveMap([[0, 1], [-1, 1]], ninc=1000)

creates a two-dimensional grid, with 1000 increments in each direction, that spans the volume 0<=x[0]<=1, -1<=x[1]<=1. This map is then trained with data f[j] corresponding to ny points y[j, d], with j=0...ny-1, (usually) uniformly distributed in y space: for example,

m.add_training_data(y, f)
m.adapt(alpha=1.5)

m.adapt(alpha=1.5) shrinks grid increments where f[j] is large, and expands them where f[j] is small. Usually one has to iterate over several sets of ys and fs before the grid has fully adapted.

The speed with which the grid adapts is determined by parameter alpha. Large (positive) values imply rapid adaptation, while small values (much less than one) imply slow adaptation. As in any iterative process that involves random numbers, it is usually a good idea to slow adaptation down in order to avoid instabilities caused by random fluctuations.

Parameters:

grid (list of arrays) – Initial x grid, where grid[d][i] is the i-th node in direction d. Different directions can have different numbers of nodes.
ninc (int or array or None) – ninc[d] (or ninc, if it is a number) is the number of increments along direction d in the new x grid. The new grid is designed to give the same Jacobian dx(y)/dy as the original grid. The default value, ninc=None, leaves the grid unchanged.

Attributes and methods:

dim¶: Number of dimensions.

ninc¶: ninc[d] is the number increments in direction d.

grid¶: The nodes of the grid defining the maps are self.grid[d, i] where d=0... specifies the direction and i=0...self.ninc the node.

inc¶

The increment widths of the grid:

self.inc[d, i] = self.grid[d, i + 1] - self.grid[d, i]

adapt(alpha=0.0, ninc=None)¶

Adapt grid to accumulated training data.

self.adapt(...) projects the training data onto each axis independently and maps it into x space. It shrinks x-grid increments in regions where the projected training data is large, and grows increments where the projected data is small. The grid along any direction is unchanged if the training data is constant along that direction.

The number of increments along a direction can be changed by setting parameter ninc (array or number).

The grid does not change if no training data has been accumulated, unless ninc is specified, in which case the number of increments is adjusted while preserving the relative density of increments at different values of x.

Parameters:

alpha (double) – Determines the speed with which the grid adapts to training data. Large (postive) values imply rapid evolution; small values (much less than one) imply slow evolution. Typical values are of order one. Choosing alpha<0 causes adaptation to the unmodified training data (usually not a good idea).
ninc (int or array or None) – The number of increments in the new grid is ninc[d] (or ninc, if it is a number) in direction d. The number is unchanged from the old grid if ninc is omitted (or equals None, which is the default).

add_training_data(y, f, ny=-1)¶

Add training data f for y-space points y.

Accumulates training data for later use by self.adapt(). Grid increments will be made smaller in regions where f is larger than average, and larger where f is smaller than average. The grid is unchanged (converged?) when f is constant across the grid.

Parameters:	y (array) – `y` values corresponding to the training data. `y` is a contiguous 2-d array, where `y[j, d]` is for points along direction `d`. f (array) – Training function values. `f[j]` corresponds to point `y[j, d]` in `y`-space. ny (int) – Number of `y` points: `y[j, d]` for `d=0...dim-1` and `j=0...ny-1`. `ny` is set to `y.shape[0]` if it is omitted (or negative).

adapt_to_samples(x, f, nitn=5, alpha=1.0, nproc=1)¶

Adapt map to data {x, f(x)}.

Replace grid with one that is optimized for integrating function f(x). New grid is found iteratively

Parameters:

x (array) – x[:, d] are the components of the sample points in direction d=0,1...self.dim-1.
f (callable or array) – Function f(x) to be adapted to. If f is an array, it is assumes to contain values f[i] corresponding to the function evaluated at points x[i].
nitn (int) – Number of iterations to use in adaptation. Default is nitn=5.
alpha (float) – Damping parameter for adaptation. Default is alpha=1.0. Smaller values slow the iterative adaptation, to improve stability of convergence.
nproc (int or None) –
Number of processes/processors to use. If nproc>1 Python’s multiprocessing module is used to spread the calculation across multiple processors. There is a significant overhead involved in using multiple processors so this option is useful mainly when very high dimenions or large numbers of samples are involved. When using the multiprocessing module in its default mode for MacOS and Windows, it is important that the main module can be safely imported (i.e., without launching new processes). This can be accomplished with some version of the if __name__ == '__main__`: construct in the main module: e.g.,
```
if __name__ == '__main__':
    main()
```
This is not an issue on other Unix platforms. See the multiprocessing documentation for more information. Set nproc=None to use all the processors on the machine (equivalent to nproc=os.cpu_count()). Default value is nproc=1. (Requires Python 3.3 or later.)

__call__(y)¶

Return x values corresponding to y.

y can be a single dim-dimensional point, or it can be an array y[i,j, ..., d] of such points (d=0..dim-1).

If y=None (default), y is set equal to a (uniform) random point in the volume.

jac(y)¶

Return the map’s Jacobian at y.

y can be a single dim-dimensional point, or it can be an array y[i,j,...,d] of such points (d=0..dim-1). Returns an array jac where jac[i,j,...] is the (multidimensional) Jacobian (dx/dy) corresponding to y[i,j,...].

jac1d(y)¶

Return the map’s Jacobian at y for each direction.

y can be a single dim-dimensional point, or it can be an array y[i,j,...,d] of such points (d=0..dim-1). Returns an array jac where jac[i,j,...,d] is the (one-dimensional) Jacobian (dx[d]/dy[d]) corresponding to y[i,j,...,d].

make_uniform(ninc=None)¶

Replace the grid with a uniform grid.

The new grid has ninc[d] (or ninc, if it is a number) increments along each direction if ninc is specified. If ninc=None (default), the new grid has the same number of increments in each direction as the old grid.

map(y, x, jac, ny=-1)¶

Map y to x, where jac is the Jacobian (dx/dy).

y[j, d] is an array of ny y-values for direction d. x[j, d] is filled with the corresponding x values, and jac[j] is filled with the corresponding Jacobian values. x and jac must be preallocated: for example,

x = numpy.empty(y.shape, float)
jac = numpy.empty(y.shape[0], float)

Parameters:

y (array) – y values to be mapped. y is a contiguous 2-d array, where y[j, d] contains values for points along direction d.
x (array) – Container for x[j, d] values corresponding to y[j, d]. Must be a contiguous 2-d array.
jac (array) – Container for Jacobian values jac[j] (= dx/dy) corresponding to y[j, d]. Must be a contiguous 1-d array.
ny (int) – Number of y points: y[j, d] for d=0...dim-1 and j=0...ny-1. ny is set to y.shape[0] if it is omitted (or negative).

invmap(x, y, jac, nx=-1)¶

Map x to y, where jac is the Jacobian (dx/dy).

y[j, d] is an array of ny y-values for direction d. x[j, d] is filled with the corresponding x values, and jac[j] is filled with the corresponding Jacobian values. x and jac must be preallocated: for example,

x = numpy.empty(y.shape, float)
jac = numpy.empty(y.shape[0], float)

Parameters:

x (array) – x values to be mapped to y-space. x is a contiguous 2-d array, where x[j, d] contains values for points along direction d.
y (array) – Container for y[j, d] values corresponding to x[j, d]. Must be a contiguous 2-d array
jac (array) – Container for Jacobian values jac[j] (= dx/dy) corresponding to y[j, d]. Must be a contiguous 1-d array
nx (int) – Number of x points: x[j, d] for d=0...dim-1 and j=0...nx-1. nx is set to x.shape[0] if it is omitted (or negative).

show_grid(ngrid=40, shrink=False)¶

Display plots showing the current grid.

Parameters:

ngrid (int) – The number of grid nodes in each direction to include in the plot. The default is 40.
axes – List of pairs of directions to use in different views of the grid. Using None in place of a direction plots the grid for only one direction. Omitting axes causes a default set of pairings to be used.
shrink – Display entire range of each axis if False; otherwise shrink range to include just the nodes being displayed. The default is False.
plotter – matplotlib plotter to use for plots; plots are not displayed if set. Ignored if None, and plots are displayed using matplotlib.pyplot.

settings(ngrid=5)¶

Create string with information about grid nodes.

Creates a string containing the locations of the nodes in the map grid for each direction. Parameter ngrid specifies the maximum number of nodes to print (spread evenly over the grid).

extract_grid()¶: Return a list of lists specifying the map’s grid.

clear()¶: Clear information accumulated by AdaptiveMap.add_training_data().

PDFIntegrator Objects¶

Expectation values using probability density functions defined by collections of Gaussian random variables (see gvar) can be evaluated using the following specialized integrator:

class vegas.PDFIntegrator(g, limit=1e15, scale=1., svdcut=1e-15)¶

vegas integrator for PDF expectation values.

Given a multi-dimensional Gaussian distribtuion g (a collection of gvar.GVars), PDFIntegrator creates a vegas integrator that evaluates expectation values of arbitrary functions f(p) where p is a point in the parameter space of the distribution. Typical usage is illustrated by:

import vegas
import gvar as gv
import numpy as np

g = gv.BufferDict()
g['a'] = gv.gvar([10., 2.], [[1, 1.4], [1.4, 2]])
g['b'] = gv.BufferDict.uniform('fb', 2.9, 3.1)

expval = vegas.PDFIntegrator(g)

def f(p):
    a = p['a']
    b = p['b']
    return a[0] + np.fabs(a[1]) ** b

result = expval(f, neval=1000, nitn=5)
print('<f(p)> =', result)

Here dictionary g specifies a three-dimensional distribution where the first two variables g['a'][0] and g['a'][1] are Gaussian with means 10 and 2, respectively, and covariance matrix [[1, 1.4], [1.4, 2.]]. The last variable g['b'] is uniformly distributed on interval [2.9, 3.1]. The result is: <f(p)> = 30.233(14).

PDFIntegrator evaluates integrals of both f(p) * pdf(p) and pdf(p), where pdf(p) is the probability density function (PDF) associated with distribution g. The expectation value of f(p) is the ratio of these two integrals (so pdf(p) need not be normalized). Integration variables are chosen to optimize the integral, and the integrator is pre-adapted to pdf(p) (so it is often unnecessary to discard early iterations). The result of a PDFIntegrator integration has an extra attribute, result.pdfnorm, which is the vegas estimate of the integral over the PDF.

The default (Gaussian) PDF associated with g can be replaced by an arbitrary PDF function pdf(p), allowing PDFIntegrator to be used for non-Gaussian distributions. In such cases, Gaussian distribution g does not affect the values of the integrals; it is used only to optimize the integration variables, and pre-adapt the integrator (and so affects the accuracy of the results). Ideally g is an approximation to the real distribution — for example, it could be the prior in a Bayesian analysis, or the result of a least-squares (or other peak-finding) analysis.

Parameters:

g – gvar.GVar, array of gvar.GVars, or dictionary whose values are gvar.GVars or arrays of gvar.GVars that specifies the multi-dimensional Gaussian distribution used to construct the probability density function. The integration variables are optimized for this function.
pdf – If specified, pdf(p) is used as the probability density function rather than the Gaussian PDF associated with g. The Gaussian PDF is used if pdf=None (default). Note that PDFs need not be normalized.
adapt_to_pdf (bool) – vegas adapts to the PDF when adapt_to_pdf=True (default). vegas adapts to pdf(p) * f(p) when calculating the expectation value of f(p) if adapt_to_pdf=False.
limit (positive float) – Limits the integrations to a finite region of size limit times the standard deviation on either side of the mean. This can be useful if the functions being integrated misbehave for large parameter values (e.g., numpy.exp overflows for a large range of arguments). Default is limit=100.
scale (positive float) – The integration variables are rescaled to emphasize parameter values of order scale times the standard deviation. The rescaling does not change the value of the integral but it can reduce uncertainties in the vegas estimate. Default is scale=1.0.
svdcut (non-negative float or None) – If not None, replace correlation matrix of g with a new matrix whose small eigenvalues are modified: eigenvalues smaller than svdcut times the maximum eigenvalue eig_max are replaced by svdcut*eig_max. This can ameliorate problems caused by roundoff errors when inverting the covariance matrix. It increases the uncertainty associated with the modified eigenvalues and so is conservative. Setting svdcut=None or svdcut=0 leaves the covariance matrix unchanged. Default is svdcut=1e-15.

All other keyword parameters are passed on to the the underlying vegas.Integrator; the uses_jac keyword is ignored.

Methods:

__call__(f, f=None, pdf=None, adapt_to_pdf=None, save=None, saveall=None, **kargs)¶

Estimate expectation value of function f(p).

Uses module vegas to estimate the integral of f(p) multiplied by the probability density function associated with g (i.e., pdf(p)). At the same time it integrates the PDF. The ratio of the two integrals is the expectation value.

Parameters:

f (function) – Function f(p) to integrate. Integral is the expectation value of the function with respect to the distribution. The function can return a number, an array of numbers, or a dictionary whose values are numbers or arrays of numbers. Setting f=None means that only the PDF is integrated. Integrals can be substantially faster if f(p) (and pdf(p) if set) are batch functions (see vegas documentation).
pdf – If specified, pdf(p) is used as the probability density function rather than the Gaussian PDF associated with g. The Gaussian PDF is used if pdf=None (default). Note that PDFs need not be normalized.
adapt_to_pdf (bool) – vegas adapts to the PDF when adapt_to_pdf=True (default). vegas adapts to pdf(p) * f(p) if adapt_to_pdf=False.
save (str or file or None) –
Writes results into pickle file specified by save at the end of each iteration. For example, setting save='results.pkl' means that the results returned by the last vegas iteration can be reconstructed later using:
```
import pickle
with open('results.pkl', 'rb') as ifile:
    results = pickle.load(ifile)
```
Ignored if save=None (default).
saveall (str or file or None) –
Writes (results, integrator) into pickle file specified by saveall at the end of each iteration. For example, setting saveall='allresults.pkl' means that the results returned by the last vegas iteration, together with a clone of the (adapted) integrator, can be reconstructed later using:
```
import pickle
with open('allresults.pkl', 'rb') as ifile:
    results, integrator = pickle.load(ifile)
```
Ignored if saveall=None (default).

All other keyword arguments are passed on to a vegas integrator; see the vegas documentation for further information.

Other Objects and Functions¶

class vegas.RAvg¶

Running average of scalar-valued Monte Carlo estimates.

This class accumulates independent Monte Carlo estimates (e.g., of an integral) and combines them into a single average. It is derived from gvar.GVar (from the gvar module if it is present) and represents a Gaussian random variable.

Different estimates are weighted by their inverse variances if parameter weight=True; otherwise straight, unweighted averages are used.

Attributes and methods:

mean¶: The mean value of the weighted average.

sdev¶: The standard deviation of the weighted average.

chi2¶: chi**2 of weighted average.

dof¶: Number of degrees of freedom in weighted average.

Q¶: Q or p-value of weighted average’s chi**2.

itn_results¶: A list of the results from each iteration.

sum_neval¶: Total number of integrand evaluations used in all iterations.

add(g)¶: Add estimate g to the running average.

summary(weighted=None)¶

Assemble summary of results, iteration-by-iteration, into a string.

Parameters:	weighted (bool) – Display weighted averages of results from different iterations if `True`; otherwise show unweighted averages. Default behavior is determined by `vegas`.

extend(ravg)¶: Merge results from RAvg object ravg after results currently in self.

class vegas.RAvgArray¶

Running average of array-valued Monte Carlo estimates.

This class accumulates independent arrays of Monte Carlo estimates (e.g., of an integral) and combines them into an array of averages. It is derived from numpy.ndarray. The array elements are gvar.GVars (from the gvar module if present) and represent Gaussian random variables.

Different estimates are weighted by their inverse covariance matrices if parameter weight=True; otherwise straight, unweighted averages are used.

Attributes and methods:

chi2¶: chi**2 of weighted average.

dof¶: Number of degrees of freedom in weighted average.

Q¶: Q or p-value of weighted average’s chi**2.

itn_results¶: A list of the results from each iteration.

sum_neval¶: Total number of integrand evaluations used in all iterations.

add(g)¶: Add estimate g to the running average.

summary(extended=False, weighted=None)¶

Assemble summary of results, iteration-by-iteration, into a string.

Parameters:	extended (bool) – Include a table of final averages for every component of the integrand if `True`. Default is `False`. weighted (bool) – Display weighted averages of results from different iterations if `True`; otherwise show unweighted averages. Default behavior is determined by `vegas`.

extend(ravg)¶: Merge results from RAvgArray object ravg after results currently in self.

class vegas.RAvgDict¶

Running average of dictionary-valued Monte Carlo estimates.

This class accumulates independent dictionaries of Monte Carlo estimates (e.g., of an integral) and combines them into a dictionary of averages. It is derived from gvar.BufferDict. The dictionary values are gvar.GVars or arrays of gvar.GVars.

Different estimates are weighted by their inverse covariance matrices if parameter weight=True; otherwise straight, unweighted averages are used.

Attributes and methods:

chi2¶: chi**2 of weighted average.

dof¶: Number of degrees of freedom in weighted average.

Q¶: Q or p-value of weighted average’s chi**2.

itn_results¶: A list of the results from each iteration.

sum_neval¶: Total number of integrand evaluations used in all iterations.

add(g)¶

summary(extended=False, weighted=None)¶

Assemble summary of results, iteration-by-iteration, into a string.

Parameters:	extended (bool) – Include a table of final averages for every component of the integrand if `True`. Default is `False`. weighted (bool) – Display weighted averages of results from different iterations if `True`; otherwise show unweighted averages. Default behavior is determined by `vegas`.

extend(ravg)¶: Merge results from RAvgDict object ravg after results currently in self.

vegas.batchintegrand()¶

Decorator for batch integrand functions.

Applying vegas.batchintegrand() to a function fcn repackages the function in a format that vegas can understand. Appropriate functions take a numpy array of integration points x[i, d] as an argument, where i=0... labels the integration point and d=0... labels direction, and return an array f[i] of integrand values (or arrays of integrand values) for the corresponding points. The meaning of fcn(x) is unchanged by the decorator.

An example is

import vegas
import numpy as np

@vegas.batchintegrand      # or @vegas.lbatchintegrand
def f(x):
    return np.exp(-x[:, 0] - x[:, 1])

for the two-dimensional integrand $\exp(-x_0 - x_1)$ .

vegas.lbatchintegrand() is the same as vegas.batchintegrand().

class vegas.BatchIntegrand¶

Wrapper for l-batch integrands.

Used by vegas.batchintegrand().

vegas.LBatchIntegrand is the same as vegas.BatchIntegrand.

vegas.rbatchintegrand()¶: Same as vegas.batchintegrand() but with batch indices on the right (not left).

class vegas.RBatchIntegrand¶: Same as vegas.BatchIntegrand but with batch indices on the right (not left).

vegas.lbatchintegrand()¶

Decorator for batch integrand functions.

Applying vegas.batchintegrand() to a function fcn repackages the function in a format that vegas can understand. Appropriate functions take a numpy array of integration points x[i, d] as an argument, where i=0... labels the integration point and d=0... labels direction, and return an array f[i] of integrand values (or arrays of integrand values) for the corresponding points. The meaning of fcn(x) is unchanged by the decorator.

An example is

import vegas
import numpy as np

@vegas.batchintegrand      # or @vegas.lbatchintegrand
def f(x):
    return np.exp(-x[:, 0] - x[:, 1])

for the two-dimensional integrand $\exp(-x_0 - x_1)$ .

vegas.lbatchintegrand() is the same as vegas.batchintegrand().

vegas.LBatchIntegrand¶: alias of vegas._vegas.BatchIntegrand

vegas.ravg(reslist, weighted=None, rescale=None)¶

Create running average from list of vegas results.

This function is used to change how the weighted average of vegas results is calculated. For example, the following code discards the first five results (where vegas is still adapting) and does an unweighted average of the last five:

import vegas

def fcn(p):
    return p[0] * p[1] * p[2] * p[3] * 16.

itg = vegas.Integrator(4 * [[0,1]])
r = itg(fcn)
print(r.summary())
ur = vegas.ravg(r.itn_results[5:], weighted=False)
print(ur.summary())

The unweighted average can be useful because it is unbiased. The output is:

itn   integral        wgt average     chi2/dof        Q
-------------------------------------------------------
 1.013(19)       1.013(19)           0.00     1.00
 0.997(14)       1.002(11)           0.45     0.50
 1.021(12)       1.0112(80)          0.91     0.40
 0.9785(97)      0.9980(62)          2.84     0.04
 1.0067(85)      1.0010(50)          2.30     0.06
 0.9996(75)      1.0006(42)          1.85     0.10
 1.0020(61)      1.0010(34)          1.54     0.16
 1.0051(52)      1.0022(29)          1.39     0.21
 1.0046(47)      1.0029(24)          1.23     0.27
 0.9976(47)      1.0018(22)          1.21     0.28

itn   integral        average         chi2/dof        Q
-------------------------------------------------------
 0.9996(75)      0.9996(75)          0.00     1.00
 1.0020(61)      1.0008(48)          0.06     0.81
 1.0051(52)      1.0022(37)          0.19     0.83
 1.0046(47)      1.0028(30)          0.18     0.91
 0.9976(47)      1.0018(26)          0.31     0.87

Parameters:

reslist (list) – List whose elements are gvar.GVars, arrays of gvar.GVars, or dictionaries whose values are gvar.GVars or arrays of gvar.GVars. Alternatively reslist can be the object returned by a call to a vegas.Integrator object (i.e, an instance of any of vegas.RAvg, vegas.RAvgArray, vegas.RAvgArray, vegas.PDFRAvg, vegas.PDFRAvgArray, vegas.PDFRAvgArray).
weighted (bool) – Running average is weighted (by the inverse covariance matrix) if True. Otherwise the average is unweighted, which makes most sense if reslist items were generated by vegas with little or no adaptation (e.g., with adapt=False). If weighted is not specified (or is None), it is set equal to getattr(reslist, 'weighted', True).
rescale – Integration results are divided by rescale before taking the weighted average if weighted=True; otherwise rescale is ignored. Setting rescale=True is equivalent to setting rescale=reslist[-1]. If rescale is not specified (or is None), it is set equal to getattr(reslist, 'rescale', True).

`vegas` Module¶

Introduction¶

Integrator Objects¶

AdaptiveMap Objects¶

PDFIntegrator Objects¶

Other Objects and Functions¶

Table of Contents

Previous topic

This Page

vegas Module¶

Introduction¶

Integrator Objects¶

AdaptiveMap Objects¶

PDFIntegrator Objects¶

Other Objects and Functions¶

`vegas` Module¶