*For reference, this work is being completed as part of my internship at [Quansi…ght Labs](https://labs.quansight.org/), which runs until the end of September.*
#### Reference issue
Addresses #20252.
#### What does this implement/fix?
This PR would enhance support for multidimensional integration ("cubature") through a new function `scipy.integrate.cub`.
One of the current gaps in the capabilities of `scipy.integrate` is support for array-valued multidimensional integrals. At the moment, there is support for multidimensional *scalar-valued* integrands through `dblquad`, `tplquad` and `nquad`, and separate support for *1D* array-valued integrands in `quad_vec`, but not for simultaneously array-valued and multidimensional integrands.
Furthermore, `dblquad`, `tplquad` and `nquad` are implemented as iterated 1D integrals using `quad`, which means integration in a moderate number of dimensions can be slow. `quad` is implemented as a wrapper around the Fortran library QUADPACK, which makes it difficult to support other array APIs like `_tanhsinh` does.
These functions and `quad_vec` also do not support custom integration rules, which means it is hard to compute tricky multidimensional integrals with singularities.
This new `cub` function has an interface similar to `quad_vec` but works towards resolving these issues by:
- Supporting array-valued multidimensional integrands (and scalar/1D integrals as a special case).
- Using an adaptive algorithm implemented in Python similar to `quad_vec`, with the eventual goal of supporting the array API. This works by subdividing the region being integrated until a tolerance is reached, rather than computing iterated 1D integrals.
- Allowing custom integration rules, which are classes with an `estimate(f, a, b)` method and `estimate_error(f, a, b)` method. There is built-in support for Cartesian products of 1D rules like Newton-Cotes, Gauss-Legendre, Gauss-Kronrod, and support for "true cubature" rules like Genz-Malik which only work in $N>1$ dimensions, but also methods available to let users write their own rules.
#### Additional information
##### Interface details
<details>
<summary>Interface details (dropdown)</summary>
`cub` is implemented in `integrate/_cubature.py` and has the signature:
```python
cub(f, a, b,
rule="gk21",
rtol=1e-05, atol=1e-08,
max_subdivisions=10000,
args=(), kwargs=None)
```
`f` can be any callable which accepts arrays of shape `(npoints, ndim)` and outputs arrays of shape `(npoints, output_dim_1, ..., output_dim_2)`. `cub` will then return an array of shape `(output_dim_1, ..., output_dim_n)` for the integral over the region specified by corners `a` and `b`, given as rank-1 arrays of length `ndim`.
`rule` is a string specifying the underlying cubature rule to use, and supports the same options as the `quadrature` option in `quad_vec`: `gk21`, `gk15` and `trapezoid`.
`rule` can also be an instance of a new class, `Rule`, which lives in `integrate/_rules/_base.py`. This is a general interface for numerical integration algorithms:
```python
class Rule:
def estimate(self, f, a, b, args=(), kwargs=None): ...
def estimate_error(self, f, a, b, args=(), kwargs=None): ...
```
There are then several subclasses. If a subclass does not implement `error_estimate`, then a default error estimator is used. The subclasses are:
- `FixedRule`, which is a rule implemented as the weighted sum of function evaluations at fixed nodes. For example, Gauss-Legendre quadrature.
- `NestedRule`, which is a rule where the error is estimated as the absolute difference between two underlying rules.
- `NestedFixedRule`, which is a rule where the error is estimated as the absolute difference between two underlying fixed rules.
- `ProductFixed`, which finds the `n`-dim cubature rule from the product of `n` 1D fixed rules
- `ProductNestedFixed`, which finds the `n`-dim cubature rule from the product of `n` 1D nested fixed rules
Implementations are then given in `integrate/_rules` for:
- `GaussKronrodQuad`, a `NestedFixedRule` implementing Gauss-Kronrod quadrature
- `GaussLegendreQuad`, a `FixedRule` implementing Gauss-Legendre quadrature
- `NewtonCotesQuad`, a `FixedRule` implementing Newton-Cotes (trapezoid) quadrature
- `GenzMalikCub`, a `FixedRule` implementing Genz-Malik cubature
This `Rule` class and its subclasses are currently not exposed as part of the public API.
</details>
##### Examples
(Examples taken from the docstring of `cub`).
<details>
<summary>Vector-valued 1D integral (dropdown)</summary>
A simple 1D integral with vector output. Here $f(x) = x^n$ is integrated over the interval $[0, 1]$ for ten values of $n$. $f$ returns an array, so all ten integrals are done simultaneously. Since no rule is specified, the default `"gk21"` is used, which corresponds to `GaussKronrodQuad(21)`. Since the spatial dimension of the integrand is 1D, this could also be done using `quad_vec`.
```python
>>> import numpy as np
>>> from scipy.integrate import cub
>>> from scipy.integrate._rules import GaussKronrodQuad
>>> def f(x, n):
... return x.reshape(-1, 1)**n # Make sure x and n are broadcastable
>>> res = cub(
... f,
... a=[0],
... b=[1],
... args=(
... # Since f accepts arrays of shape (npoints, ndim) we need to
... # make sure n is the right shape
... np.arange(10).reshape(1, -1),
... )
... )
>>> res.estimate
array([1. , 0.5 , 0.33333333, 0.25 , 0.2 ,
0.16666667, 0.14285714, 0.125 , 0.11111111, 0.1 ])
```
</details>
<details>
<summary>Tensor-valued 7D integral (dropdown)</summary>
A 7D integral with arbitrary-shaped array output. Here $f(\pmb x) = \cos(2\pi r + \pmb \alpha^\top \pmb x)$ for some $r$ and $\pmb \alpha$, and the integral is performed over the unit hybercube, ``[0, 1]^7``.
```python
>>> import numpy as np
>>> from scipy.integrate import cub
>>> from scipy.integrate._rules import GenzMalikCub
>>> def f(x, r, alphas):
... # f(x) = cos(2*pi*r + alpha @ x)
... # Allow r and alphas to be arbitrary shape
... npoints, ndim = x.shape[0], x.shape[-1]
... alphas_reshaped = alphas[np.newaxis, :]
... x_reshaped = x.reshape(npoints, *([1]*(len(alphas.shape) - 1)), ndim)
... return np.cos(2*np.pi*r + np.sum(alphas_reshaped * x_reshaped, axis=-1))
>>> np.random.seed(1)
>>> res = cub(
... f=f,
... a=np.array([0, 0, 0, 0, 0, 0, 0]),
... b=np.array([1, 1, 1, 1, 1, 1, 1]),
... rule=GenzMalikCub(7),
... kwargs={
... "r": np.random.rand(2, 3),
... "alphas": np.random.rand(2, 3, 7),
... }
... )
>>> res.estimate
array([[-0.61336595, 0.88388877, -0.57997549],
[-0.86968418, -0.86877137, -0.64602074]])
```
</details>
##### Difficulties and next steps
It would be great to get some suggestions on the interface, especially around having `NestedRule` and `NestedFixedRule`.
Currently the following are not supported:
- Infinite endpoints
- Accepting a `points` parameter like `quad_vec`
- Broadcastable endpoints
- Non-rectangular regions
It is our hope that this will be added before the end of the internship (end of September), depending on how much further work is required for this PR.
A future PR can also update `quad_vec` to accept the 1D rule interface proposed here, perhaps just by using its existing "quadrature" argument. This would allow for plug-and-play rules in adaptive quadrature and would work nicely with integration methods such as tanh-sinh.
Many thanks to @izaid and @ev-br, who have given lots of valuable support for this PR.