timates of the integral over smaller
        regions of the domain.

    Each object in ``regions`` has the following attributes:

    a, b : ndarray
        Points describing the corners of the region. If the original integral
        contained infinite limits or was over a region described by `region`,
        then `a` and `b` are in the transformed coordinates.
    estimate : ndarray
        Estimate of the value of the integral over this region.
    error : ndarray
        Estimate of the error of the approximation over this region.

Notes
-----
The algorithm uses a similar algorithm to `quad_vec`, which itself is based on the
implementation of QUADPACK's DQAG* algorithms, implementing global error control and
adaptive subdivision.

The source of the nodes and weights used for Gauss-Kronrod quadrature can be found
in [1]_, and the algorithm for calculating the nodes and weights in Genz-Malik
cubature can be found in [2]_.

The rules currently supported via the `rule` argument are:

- ``"gauss-kronrod"``, 21-node Gauss-Kronrod
- ``"genz-malik"``, n-node Genz-Malik

If using Gauss-Kronrod for an ``n``-dim integrand where ``n > 2``, then the
corresponding Cartesian product rule will be found by taking the Cartesian product
of the nodes in the 1D case. This means that the number of nodes scales
exponentially as ``21^n`` in the Gauss-Kronrod case, which may be problematic in a
moderate number of dimensions.

Genz-Malik is typically less accurate than Gauss-Kronrod but has much fewer nodes,
so in this situation using "genz-malik" might be preferable.

Infinite limits are handled with an appropriate variable transformation. Assuming
``a = [a_1, ..., a_n]`` and ``b = [b_1, ..., b_n]``:

If :math:`a_i = -\infty` and :math:`b_i = \infty`, the i-th integration variable
will use the transformation :math:`x = \frac{1-|t|}{t}` and :math:`t \in (-1, 1)`.

If :math:`a_i \ne \pm\infty` and :math:`b_i = \infty`, the i-th integration variable
will use the transformation :math:`x = a_i + \frac{1-t}{t}` and
:math:`t \in (0, 1)`.

If :math:`a_i = -\infty` and :math:`b_i \ne \pm\infty`, the i-th integration
variable will use the transformation :math:`x = b_i - \frac{1-t}{t}` and
:math:`t \in (0, 1)`.

References
----------
.. [1] R. Piessens, E. de Doncker, Quadpack: A Subroutine Package for Automatic
    Integration, files: dqk21.f, dqk15.f (1983).

.. [2] A.C. Genz, A.A. Malik, Remarks on algorithm 006: An adaptive algorithm for
    numerical integration over an N-dimensional rectangular region, Journal of
    Computational and Applied Mathematics, Volume 6, Issue 4, 1980, Pages 295-302,
    ISSN 0377-0427
    :doi:`10.1016/0771-050X(80)90039-X`

Examples
--------
**1D integral with vector output**:

.. math::

    \int^1_0 \mathbf f(x) \text dx

Where ``f(x) = x^n`` and ``n = np.arange(10)`` is a vector. Since no rule is
specified, the default "gk21" is used, which corresponds to Gauss-Kronrod
integration with 21 nodes.

>>> import numpy as np
>>> from scipy.integrate import cubature
>>> def f(x, n):
...    # Make sure x and n are broadcastable
...    return x[:, np.newaxis]**n[np.newaxis, :]
>>> res = cubature(
...     f,
...     a=[0],
...     b=[1],
...     args=(np.arange(10),),
... )
>>> res.estimate
 array([1.        , 0.5       , 0.33333333, 0.25      , 0.2       ,
        0.16666667, 0.14285714, 0.125     , 0.11111111, 0.1       ])

**7D integral with arbitrary-shaped array output**::

    f(x) = cos(2*pi*r + alphas @ x)

for some ``r`` and ``alphas``, and the integral is performed over the unit
hybercube, :math:`[0, 1]^7`. Since the integral is in a moderate number of
dimensions, "genz-malik" is used rather than the default "gauss-kronrod" to
avoid constructing a product rule with :math:`21^7 \approx 2 \times 10^9` nodes.

>>> import numpy as np
>>> from scipy.integrate import cubature
>>> def f(x, r, alphas):
...     # f(x) = cos(2*pi*r + alphas @ x)
...     # Need to allow r and alphas to be arbitrary shape
...     npoints, ndim = x.shape[0], x.shape[-1]
...     alphas = alphas[np.newaxis, ...]
...     x = x.reshape(npoints, *([1]*(len(alphas.shape) - 1)), ndim)
...     return np.cos(2*np.pi*r + np.sum(alphas * x, axis=-1))
>>> rng = np.random.default_rng()
>>> r, alphas = rng.random((2, 3)), rng.random((2, 3, 7))
>>> res = cubature(
...     f=f,
...     a=np.array([0, 0, 0, 0, 0, 0, 0]),
...     b=np.array([1, 1, 1, 1, 1, 1, 1]),
...     rtol=1e-5,
...     rule="genz-malik",
...     args=(r, alphas),
... )
>>> res.estimate
 array([[-0.79812452,  0.35246913, -0.52273628],
        [ 0.88392779,  0.59139899,  0.41895111]])

**Parallel computation with** `workers`:

>>> from concurrent.futures import ThreadPoolExecutor
>>> with ThreadPoolExecutor() as executor:
...     res = cubature(
...         f=f,
...         a=np.array([0, 0, 0, 0, 0, 0, 0]),
...         b=np.array([1, 1, 1, 1, 1, 1, 1]),
...         rtol=1e-5,
...         rule="genz-malik",
...         args=(r, alphas),
...         workers=executor.map,
...      )
>>> res.estimate
 array([[-0.79812452,  0.35246913, -0.52273628],
        [ 0.88392779,  0.59139899,  0.41895111]])

**2D integral with infinite limits**:

.. math::

    \int^{ \infty }_{ -\infty }
    \int^{ \infty }_{ -\infty }
        e^{-x^2-y^2}
    \text dy
    \text dx

>>> def gaussian(x):
...     return np.exp(-np.sum(x**2, axis=-1))
>>> res = cubature(gaussian, [-np.inf, -np.inf], [np.inf, np.inf])
>>> res.estimate
 3.1415926

**1D integral with singularities avoided using** `points`:

.. math::

    \int^{ 1 }_{ -1 }
      \frac{\sin(x)}{x}
    \text dx

It is necessary to use the `points` parameter to avoid evaluating `f` at the origin.

>>> def sinc(x):
...     return np.sin(x)/x
>>> res = cubature(sinc, [-1], [1], points=[[0]])
>>> res.estimate
 1.8921661
NÚ