Home

Awesome

Cubature

This is a simple package for adaptive multidimensional integration (cubature) of vector-valued integrands over hypercubes, originally written by Steven G. Johnson and ported to Java by Jonathan Schilling. That is, it computes integrals of the form:

n-dimensional integral

(Of course, it can handle scalar integrands as the special case where f is a one-dimensional vector: the dimensionalities of f and x are independent.) The integrand can be evaluated for an array of points at once to enable easy parallelization. The code, which is distributed as free software under the terms of the GNU General Public License (v2 or later), implements an algorithm for adaptive integration.

The h-adaptive integration algorithm (recursively partitioning the integration domain into smaller subdomains, applying the same integration rule to each, until convergence is achieved), is based on the algorithms described in:

This algorithm is best suited for a moderate number of dimensions (say, < 7), and is superseded for high-dimensional integrals by other methods (e.g. Monte Carlo variants or sparse grids).

(Note that we do not use any of the original DCUHRE code by Genz, which is not under a free/open-source license.) Our code is based in part on code borrowed from the HIntLib numeric-integration library by Rudolf Schürer and from code for Gauss-Kronrod quadrature (for 1d integrals) from the GNU Scientific Library, both of which are free software under the GNU GPL. (Another free-software multi-dimensional integration library, unrelated to our code here but also implementing the Genz–Malik algorithm among other techniques, is Cuba.)

We are also grateful to Dmitry Turbiner (dturbiner ατ alum.mit.edu), who implemented an initial prototype of the “vectorized” functionality (see below) for evaluating an array of points in a single call, which facilitates parallelization of the integrand evaluation.

Download

The current version of the code is Cubature.java on Github, and the latest "official" version is available in Maven Central:

<dependency>
	<groupId>de.labathome</groupId>
	<artifactId>Cubature</artifactId>
	<version>1.2.0</version>
</dependency>

The TestCubature.java file contains some JUnit tests.

Usage

The central subroutine you will be calling for h-adaptive cubature is:

public static double[][] integrate(
	UnaryOperator<double[][]> integrand,
	double[] xmin, double[] xmax,
	double relTol, double absTol, Error norm,
	int maxEval);

This integrates a function F(x), returning a vector of fdim integrands, where x is a dim-dimensional vector ranging from xmin to xmax (i.e. in a hypercube xminᵢ ≤ xᵢ ≤ xmaxᵢ).

maxEval specifies a maximum number of function evaluations (0 for no limit). (Note: the actual number of evaluations may somewhat exceed maxEval: maxEval is rounded up to an integer number of subregion evaluations.) Otherwise, the integration stops when the estimated |error| is less than absTol (the absolute error requested) or when the estimated |error| is less than relTol × |integral value| (the relative error requested). (Either of the error tolerances can be set to Double.NaN to ignore it.)

For vector-valued integrands (fdim > 1), norm specifies the norm that is used to measure the error and determine convergence properties. (The norm argument is irrelevant for fdim ≤ 1 and is ignored.) Given vectors v and e of estimated integrals and errors therein, respectively, the norm argument takes on one of the following enumerated constant values:

The integrate subroutine returns an array of dimensions [2][fdim], in which the integral value and the associated error estimates are stored. The estimated errors are based on an embedded cubature rule of lower order; for smooth functions, this estimate is usually conservative.

The integrand function F should have the following signature (or be a UnaryOperator<double[][]>):

double[][] eval(final double[][] x)

Here, the input is a matrix x of dimensions [dim][nPoints] (the points to be evaluated), the output is an array fval of dimensions [fdim][nPoints] (the vectors of function values at the point x).

“Vectorized” interface

This integration algorithm actually evaluates the integrand in “batches” of several points at a time. It is often useful to have access to this information so that your integrand function is not called for one point at a time, but rather for a whole “vector” of many points at once.

The value of nPoints will vary with the dimensionality of the problem; higher-dimensional problems will have (exponentially) larger nPoints, allowing for the possibility of more parallelism. Currently, for integrate, nPoints starts at 15 in 1d, 17 in 2d, and 33 in 3d, but as integrate calls your integrand more and more times the value of nPoints will grow. e.g. if you end up requiring several thousand points in total, nPoints may grow to several hundred. We utilize an algorithm from:

as described in the article “Parallel globally adaptive algorithms for multi-dimensional integration” by Bull and Freeman (1994).

Example

As a simple example, consider the Gaussian integral of the scalar function f(x) = exp(-sigma |x|²) over the hypercube [-2,2]³ in 3 dimensions. Using Cubature, you can compute this integral e.g. using a static member function as shown in ThreeDimGaussianStaticExample.java:

import java.util.Locale;
import java.util.function.UnaryOperator;

import de.labathome.cubature.Cubature;
import de.labathome.cubature.CubatureError;

public class ThreeDimGaussianExample {

    public static void ex_ThreeDimGaussian() {

        double[] xmin = { -2.0, -2.0, -2.0 };
        double[] xmax = { 2.0, 2.0, 2.0 };

        double sigma = 0.5;

        UnaryOperator<double[][]> gaussianNd = (double[][] x) -> {
            int dim = x.length;
            int nPoints = x[0].length;
            double[][] fval = new double[1][nPoints];
            for (int i = 0; i < nPoints; ++i) {
                double sum = 0.0;
                for (int d = 0; d < dim; ++d) {
                    sum += x[d][i] * x[d][i];
                }
                fval[0][i] = Math.exp(-sigma * sum);
            }
            return fval;
        };

        double[][] val_err = Cubature.integrate(gaussianNd,
                xmin, xmax,
                1.0e-4, 0.0, CubatureError.INDIVIDUAL,
                0);

        System.out.println(String.format(Locale.ENGLISH,
                "Computed integral = %.8f +/- %g", val_err[0][0], val_err[1][0]));
    }

    public static void main(String[] args) {
        ex_ThreeDimGaussian();
    }
}

or using an object's member function as shown in ThreeDimGaussianMemberExample.java:

Here, we have specified a relative error tolerance of $10^{-4}$ (and no absolute error tolerance or maximum number of function evaluations).

The output should be:

Computed integral = 13.69609043 +/- 0.00136919

Note that the estimated relative error is 0.00136919/13.69609043 = 9.9969×10⁻⁵, within our requested tolerance of 10⁻⁴. The actual error in the integral value, as can be determined e.g. by running the integration with a much lower tolerance, is much smaller: the integral is too small by about 0.00002, for an actual relative error of about 1.4×10⁻⁶. As mentioned above, for smooth integrands the estimated error is almost always conservative (which means, unfortunately, that the integrator usually does more function evaluations than it needs to).

Infinite intervals

Integrals over infinite or semi-infinite intervals is possible by a change of variables. This is best illustrated in one dimension.

To compute an integral over a semi-infinite interval, you can perform the change of variables x=a+t/(1-t):

semi-infinite integral

For an infinite interval, you can perform the change of variables x=t/(1-t²):

infinite integral

Note the Jacobian factors multiplying f(⋅⋅⋅) in both integrals, and also that the limits of the t integrals are different in the two cases.

In multiple dimensions, one simply performs this change of variables on each dimension separately, as desired, multiplying the integrand by the corresponding Jacobian factor for each dimension being transformed.

The Jacobian factors diverge as the endpoints are approached. However, if f(x) goes to zero at least as fast as 1/x², then the limit of the integrand (including the Jacobian factor) is finite at the endpoints. If your f(x) vanishes more slowly than 1/x² but still faster than 1/x, then the integrand blows up at the endpoints but the integral is still finite (it is an integrable singularity), so the code will work (although it may take many function evaluations to converge). If your f(x) vanishes only as 1/x, then it is not absolutely convergent and much more care is required even to define what you are trying to compute. In any case, the h-adaptive quadrature/cubature rules currently employed in Cubature.java do not evaluate the integrand at the endpoints, so you need not implement special handling for |t|=1.

Tests

JUnit tests should be added based on the test.c suite of the original cubature package.

The different test integrands are: