quadax.adaptive_cubature

quadax.adaptive_cubature(rule: AbstractCubatureRule, fun: Callable[[...], Array], interval: Array | ndarray | bool | number | bool | int | float | complex | Sequence[Array | ndarray | bool | number | bool | int | float | complex | Callable], args: tuple = (), full_output: bool = False, epsabs: Any = None, epsrel: Any = None, max_nregion: int = 1000, adjoint: AbstractAdjoint = DirectAdjoint(), extrapolate: bool = True, throw: bool = False, **kwargs)Source

Global adaptive cubature with user specified local rule.

This is a lower level routine allowing for custom local cubature rules. For most applications cubegm() is preferable; this is the way to reach a TensorProductRule, which is worth its cost mainly in two or three dimensions or where one axis is far harder than the others.

Parameters:

  • rule (AbstractCubatureRule) – Local cubature rule to use. Its ndim must match interval.

  • 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.

  • 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.

  • 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.

  • kwargs (dict, optional) – Additional keyword arguments passed to rule.integrate.

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 rule 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.