Skip to content

Areustle/cubepy

Folders and files

NameName
Last commit message
Last commit date

Latest commit

9069d27 Β· May 2, 2023

History

49 Commits
May 2, 2023
May 2, 2023
May 2, 2023
May 2, 2023
Dec 14, 2021
Dec 28, 2021
May 2, 2023
May 2, 2023
May 2, 2023
Dec 30, 2021
Dec 28, 2021
May 2, 2023

Repository files navigation

🧊 CubePy 🐍

Batched Numerical Cubature: Adaptive Multi-dimensional Integration in Python.

Finally, a fully vectorized multi-dimensional numerical integrator! CubePy is the High-Performance, deterministic numerical integrator you've been waiting for.

Integrate functions with:

  • High Dimensional Domains
  • Multiple limits of integration
  • Independent batches computed concurrently.
  • Deterministic Computation (No Monte-Carlo!)

CubePy is a fully-vectorized python package for performing numerical integration on multi-dimensional vector functions. CubePy performs these operations efficiently using numpy and the Genz-Malik Adaptive cubature algorithm.

Functions with 2+D domain use an adaptive Genz Malik Cubature scheme. Functions with 1D domain use an adaptive Gauss Kronrod Quadrature scheme.

The adaptive regional subdivision is performed independently for un-converged regions. Local convergence is obtained when the global tolerance values exceed a local region's absolute and relative error estimate Using Bernsten's Error Estimation formula. Global convergence is obtained when all regions are locally converged.

Example Usage

The volume of a sphere, a 3 dimensional integral.

## The analytical volume of a sphere for validation.
def volume_of_sphere(radius):
    return (4.0 / 3.0) * np.pi * radius**3

## The integrand function.
def sphere_integrand(r, theta, phi):
    return np.sin(phi) * r**2

## The upper and lower bounds of integration.
low = [0.0, 0.0, 0.0]
high = [1.0, 2 * np.pi, np.pi]

# The Integration!
value, error = cp.integrate(sphere_integrand, low, high)

# Confirm the results
assert np.allclose(value, volume_of_sphere(1.0))

# Try a larger sphere
high[0] = 42
value, error = cp.integrate(sphere_integrand, low, high)
assert np.allclose(value, volume_of_sphere(42))

The Volume of multiple spheres in a single call

Frequently we encounter the need to compute an integral on an array of bounds, representing different events in our dataset. CubePy was designed to handle these exact cases. Simply define one or more of your integration bounds as an array of entries and call cp.integrate

## 1,000,000 radii between 1 and 42
radii = np.linspace(1, 42, int(1e6))

## The upper and lower bounds of integration.
low = [0.0, 0.0, 0.0]
high = [radii, 2 * np.pi, np.pi]

# The Integration!
value, error = cp.integrate(sphere_integrand, low, high)
# Confirm the results
assert np.allclose(value, volume_of_sphere(radii))

Single Dimensional integral support

The routine supports single dimensional integrals as well, with an adaptive Gauss Kronrod scheme.

def degree_20_polynomial(x):
    return 2.0 * np.pi * x**20 - np.e * x**3 + 3.0 * x**2 - 4.0 * x + 8.0

def exact_degree_20_polynomial(a, b):
    def exact(x):
        return (
            (2.0 * np.pi / 21.0) * x**21
            - np.e / 4.0 * x**4
            + x**3
            - 2.0 * x**2
            + 8.0 * x
        )

    return exact(b) - exact(a)

v, _ = cp.integrate(degree_20_polynomial, -1.0, 1.0)
np.allclose(v, exact_degree_20_polynomial(-1, 1))

v, _ = cp.integrate(degree_20_polynomial, -np.pi, np.pi)
np.allclose(v, exact_degree_20_polynomial(-np.pi, np.pi))

References

  • Berntsen, Espelid, Genz. An Adaptive Algorithm for the Approximate Calculation of Multiple Integrals.
  • Berntsen. Practical error estimation in adaptive multidimensional quadrature routines.
  • Genz, Malik. An adaptive algorithm for numerical integration over an N-dimensional rectangular region.
  • Genz, Malik. An imbedded family of fully symmetric numerical integration rules.
  • Van Dooren, De Ridder. An adaptive algorithm for numerical in- tegration over an N-dimensional cube.