quadax.cubegm

Global adaptive cubature over a box using a Genz-Malik rule.

Integrate fun over the box interval using an h-adaptive scheme with error estimate. Breakpoints can be specified per axis in interval where integration difficulty may occur.

The general purpose cubature to reach for first. Each subdivision cuts one region along one axis, chosen by the rule itself, so a difficulty that lies along an axis aligned plane is refined without the axes it spans being refined with it.

Where the integrand has a jump or a singularity on a known axis aligned plane, passing that coordinate as a breakpoint on its axis is worth more than any change of method, since the subdivision no longer has to find it.

Parameters:

fun (callable) – Function to integrate, should have a signature of the form fun(x, *args) -> float, Array, where x has shape (ndim,). Should be JAX transformable.
interval (array-like or sequence of array-like or callable) –
Limits of integration. Either an array of shape (ndim, 2), row k giving the two limits of axis k, or a list or tuple of one entry per axis, each holding that axis’ two limits with any breakpoints between them. The second form is what carries breakpoints, and the axes may carry different numbers of them. ndim is taken from the length of either. Use np.inf to denote unbounded axes. The dtype sets the working precision: the integrand is called with an x of this dtype, and the result follows it unless the integrand upcasts. Integer types or python floats fall back to the JAX default. Must be real; complex integrands are supported, complex limits are not.

An entry of the list or tuple form may instead be a callable lim(x_prev, *args) returning those same limits, where x_prev has shape (k,) and holds the coordinates of the axes before it, so that the region need not be a box. Axis k may depend only on axes 0 to k-1, which makes the order of the axes significant. Any breakpoints the callable returns may move with the coordinates too, which is how a feature lying along a curve is marked.
args (tuple, optional) – Extra arguments passed to fun and any callable interval limit.
full_output (bool, optional) – If True, return the full state of the integrator. See below for more information.
epsabs (float, optional) – Absolute and relative error tolerance. Default is the square root of the machine precision of the working dtype. Algorithm tries to obtain an accuracy of abs(i-result) <= max(epsabs, epsrel*abs(i)) where i = integral of fun over interval, and result is the numerical approximation.
epsrel (float, optional) – Absolute and relative error tolerance. Default is the square root of the machine precision of the working dtype. Algorithm tries to obtain an accuracy of abs(i-result) <= max(epsabs, epsrel*abs(i)) where i = integral of fun over interval, and result is the numerical approximation.
max_nregion (int, optional) – An upper bound on the number of regions used in the adaptive algorithm.
degree (int) – Degree of polynomial the local rule integrates exactly, one of 7, 9, 11, 13. Higher degrees cost sharply more per region as ndim grows; see GenzMalikRule for the counts and for what the choice buys.
norm (int, callable) – Norm to use for measuring error for vector valued integrands. No effect if the integrand is scalar valued. If an int, uses p-norm of the given order, otherwise should be callable.
adjoint (AbstractAdjoint, optional) – How to compute derivatives of the cubature. DirectAdjoint, the default, differentiates the subdivision the solve settled on and is usually the cheaper choice. LeibnizAdjoint gives the derivative its own error control, at the cost of integrating over the faces of the box, which is an adaptive solve of one dimension fewer per axis whose limits move.
extrapolate (bool, optional) – Whether to accelerate convergence by applying Wynn’s epsilon algorithm to the sequence of running totals, on by default. Not needed for smooth integrands on finite domains, but can help significantly if there are algebraic singularities or infinite intervals. Can be turned off for integrands known to be smooth which can reduce overhead and improve wall time.
batch_size (int, optional) – Maximum number of points at which to evaluate the integrand in parallel. Default is all of the local rule’s nodes at once, which is fastest but makes peak memory scale with the node count. Lower it to reduce memory on an expensive integrand.
throw (bool, optional) – Whether to raise an error if the routine does not converge. If True, a run that terminates for any reason other than reaching the requested tolerance raises with the message its status carries. If False, the default, that status is reported on the returned info and left to the caller to act on.

Returns:

y (float, Array) – The integral of fun over the box.
info (QuadratureInfo) – Named tuple with the following fields:
- err : (float) Estimate of the error in the approximation.
- neval : (int) Total number of function evaluations.
- status : (quadax.STATUS) Why the routine terminated. STATUS.normal means the requested tolerances were reached; every other member names a difficulty and prints as the message explaining it.
- info : (dict or None) Other information returned by the algorithm. Only present if full_output is True. Contains the following:
  - ’nregion’ : (int) The number, K, of regions produced in the subdivision.
  - ’a_arr’, ‘b_arr’ : (ndarray) shape (max_nregion, ndim), the first K rows of which are the opposite corners of the (remapped) regions.
  - ’r_arr’ : (ndarray) the integral approximations on those regions.
  - ’e_arr’ : (ndarray) the absolute error estimates on those regions.

Notes

An unbounded axis is mapped to a finite one by the same substitution the one dimensional routines use, applied to each axis on its own. Because that map is separable its Jacobian is a product over axes, which for an integrand decaying algebraically in |x| disagrees with the decay along different directions into a corner of the mapped box. The value is unaffected, the map being a bijection onto the original domain, but the mesh converges on it at first order there rather than at the rule’s degree, and such a problem costs far more regions than its smoothness suggests. An integrand decaying exponentially, or one over an axis unbounded on only one side, is not affected. No breakpoint helps, the difficulty lying at a point the map sends to infinity.

An axis whose limits are given as a callable is integrated over a fixed interval instead, and the map onto the limits it returns is folded into the integrand along with its Jacobian. The region is then a box again, and everything above applies to it unchanged; a limit the callable returns may be unbounded like any other. What such a region costs is smoothness: a curved boundary moves its curvature into the axes before it, so a disc costs far more subdivision than a triangle of the same area does, and no breakpoint helps there either, the difficulty being a square root rather than a kink. The limit functions are evaluated once per node per axis they govern.