# Integral algorithms

`AutoBZCore.IntegralAlgorithm`

— Type`IntegralAlgorithm`

Abstract supertype for integration algorithms.

## Quadrature

`AutoBZCore.QuadratureFunction`

— Type`QuadratureFunction(; fun=trapz, npt=50, nthreads=1)`

Quadrature rule for the standard interval [-1,1] computed from a function `x, w = fun(npt)`

. The nodes and weights should be set so the integral of `f`

on [-1,1] is `sum(w .* f.(x))`

. The default quadrature rule is `trapz`

, although other packages provide rules, e.g.

```
using FastGaussQuadrature
alg = QuadratureFunction(fun=gausslegendre, npt=100)
```

`nthreads`

sets the numbers of threads used to parallelize the quadrature only when the integrand is a `BatchIntegrand`

, in which case the user must parallelize the integrand evaluations. For no threading set `nthreads=1`

.

`AutoBZCore.QuadGKJL`

— Type`QuadGKJL(; order = 7, norm = norm)`

Duplicate of the QuadGKJL provided by Integrals.jl.

`AutoBZCore.AuxQuadGKJL`

— Type`AuxQuadGKJL(; order = 7, norm = norm)`

Generalization of the QuadGKJL provided by Integrals.jl that allows for `AuxValue`

d integrands for auxiliary integration and multi-threaded evaluation with the `batch`

argument to `IntegralProblem`

`AutoBZCore.ContQuadGKJL`

— Type`ContQuadGKJL(; order = 7, norm = norm, rho = 1.0, rootmeth = IteratedIntegration.ContQuadGK.NewtonDeflation())`

A 1d contour deformation quadrature scheme for scalar, complex-valued integrands. It defaults to regular `quadgk`

behavior on the real axis, but if it finds a root of 1/f nearby, in the sense of Bernstein ellipse for the standard segment `[-1,1]`

with semiaxes `cosh(rho)`

and `sinh(rho)`

, on either the upper/lower half planes, then it dents the contour away from the presumable pole.

`AutoBZCore.MeroQuadGKJL`

— Type`MeroQuadGKJL(; order = 7, norm = norm, rho = 1.0, rootmeth = IteratedIntegration.MeroQuadGK.NewtonDeflation())`

A 1d pole subtraction quadrature scheme for scalar, complex-valued integrands that are meromorphic. It defaults to regular `quadgk`

behavior on the real axis, but if it finds nearby roots of 1/f, in the sense of Bernstein ellipse for the standard segment `[-1,1]`

with semiaxes `cosh(rho)`

and `sinh(rho)`

, it attempts pole subtraction on that segment.

## Cubature

`AutoBZCore.HCubatureJL`

— Type`HCubatureJL(; norm=norm, initdiv=1)`

A copy of `HCubatureJL`

from Integrals.jl.

`AutoBZCore.MonkhorstPack`

— Type`MonkhorstPack(; npt=50, syms=nothing, nthreads=1)`

Periodic trapezoidal rule with a fixed number of k-points per dimension, `npt`

, using the `PTR`

rule from AutoSymPTR.jl. `nthreads`

sets the numbers of threads used to parallelize the quadrature only when the integrand is a `BatchIntegrand`

, in which case the user must parallelize the integrand evaluations. For no threading set `nthreads=1`

. **The caller should check that the integral is converged w.r.t. npt**.

`AutoBZCore.AutoSymPTRJL`

— Type`AutoSymPTRJL(; norm=norm, a=1.0, nmin=50, nmax=1000, n₀=6, Δn=log(10), keepmost=2, nthreads=1)`

Periodic trapezoidal rule with automatic convergence to tolerances passed to the solver with respect to `norm`

using the routine `autosymptr`

from AutoSymPTR.jl. `nthreads`

sets the numbers of threads used to parallelize the quadrature only when the integrand is a `BatchIntegrand`

, in which case the user must parallelize the integrand evaluations. For no threading set `nthreads=1`

. **This algorithm is the most efficient for smooth integrands**.

## Meta-algorithms

`AutoBZCore.NestedQuad`

— Type```
NestedQuad(alg::IntegralAlgorithm)
NestedQuad(algs::IntegralAlgorithm...)
```

Nested integration by repeating one quadrature algorithm or composing a list of algorithms. The domain of integration must be an `AbstractIteratedLimits`

from the IteratedIntegration.jl package. Analogous to `nested_quad`

from IteratedIntegration.jl. The integrand should expect `SVector`

inputs. Do not use this for very high-dimensional integrals, since the compilation time scales very poorly with respect to dimensionality. In order to improve the compilation time, FunctionWrappers.jl is used to enforce type stability of the integrand, so you should always pick the widest integration limit type so that inference works properly. For example, if `ContQuadGKJL`

is used as an algorithm in the nested scheme, then the limits of integration should be made complex.

`AutoBZCore.EvalCounter`

— Type`EvalCounter(::IntegralAlgorithm)`

An algorithm which counts the evaluations used by another algorithm. The count is stored in the `sol.numevals`

field.

`AutoBZCore.AbsoluteEstimate`

— Type`AbsoluteEstimate(est_alg, abs_alg; kws...)`

Most algorithms are efficient when using absolute error tolerances, but how do you know the size of the integral? One option is to estimate it using second algorithm.

A multi-algorithm to estimate an integral using an `est_alg`

to generate a rough estimate of the integral that is combined with a user's relative tolerance to re-calculate the integral to higher accuracy using the `abs_alg`

. The keywords passed to the algorithm may include `reltol`

, `abstol`

and `maxiters`

and are given to the `est_alg`

solver. They should limit the amount of work of `est_alg`

so as to only generate an order-of-magnitude estimate of the integral. The tolerances passed to `abs_alg`

are `abstol=max(abstol,reltol*norm(I))`

and `reltol=0`

.

# BZ-specific integral algorithms

In order to make algorithms domain-agnostic, the BZ loaded from `load_bz`

can be called with the algorithms below, which are wrappers for algorithms above with the additional capability of mapping integrals over the IBZ to the FBZ.

`AutoBZCore.AutoBZAlgorithm`

— Type`AutoBZAlgorithm`

Abstract supertype for Brillouin zone integration algorithms. All integration problems on the BZ get rescaled to fractional coordinates so that the Brillouin zone becomes `[0,1]^d`

, and integrands should have this periodicity. If the integrand depends on the Brillouin zone basis, then it may have to be transformed to the Cartesian coordinates as a post-processing step. These algorithms also use the symmetries of the Brillouin zone and the integrand.

`AutoBZCore.IAI`

— Type```
IAI(alg::IntegralAlgorithm=AuxQuadGKJL())
IAI(algs::IntegralAlgorithm...)
```

Iterated-adaptive integration using `nested_quad`

from IteratedIntegration.jl. **This algorithm is the most efficient for localized integrands**.

`AutoBZCore.TAI`

— Type`TAI(; norm=norm, initdivs=1)`

Tree-adaptive integration using `hcubature`

from HCubature.jl. This routine is limited to integration over hypercube domains and may not use all symmetries.

`AutoBZCore.PTR`

— Type`PTR(; npt=50, nthreads=1)`

Periodic trapezoidal rule with a fixed number of k-points per dimension, `npt`

, using the routine `ptr`

from AutoSymPTR.jl. **The caller should check that the integral is converged w.r.t. npt**.

`AutoBZCore.AutoPTR`

— Type`AutoPTR(; norm=norm, a=1.0, nmin=50, nmax=1000, n₀=6, Δn=log(10), keepmost=2, nthreads=1)`

Periodic trapezoidal rule with automatic convergence to tolerances passed to the solver with respect to `norm`

using the routine `autosymptr`

from AutoSymPTR.jl. **This algorithm is the most efficient for smooth integrands**.

`AutoBZCore.PTR_IAI`

— Function`PTR_IAI(; ptr=PTR(), iai=IAI())`

Multi-algorithm that returns an `IAI`

calculation with an `abstol`

determined from the given `reltol`

and a `PTR`

estimate, `I`

, as `reltol*norm(I)`

. This addresses the issue that `IAI`

does not currently use a globally-adaptive algorithm and may not have the expected scaling with localization length unless an `abstol`

is used since computational effort may be wasted via a `reltol`

with the naive `nested_quadgk`

.

`AutoBZCore.AutoPTR_IAI`

— Function`AutoPTR_IAI(; reltol=1.0, ptr=AutoPTR(), iai=IAI())`

Multi-algorithm that returns an `IAI`

calculation with an `abstol`

determined from an `AutoPTR`

estimate, `I`

, computed to `reltol`

precision, and the `rtol`

given to the solver as `abstol=rtol*norm(I)`

. This addresses the issue that `IAI`

does not currently use a globally-adaptive algorithm and may not have the expected scaling with localization length unless an `abstol`

is used since computational effort may be wasted via a `reltol`

with the naive `nested_quadgk`

.