Quadrature.jl is an instantiation of the DiffEqBase.jl common QuadratureProblem
interface for the common quadrature packages of Julia. By using Quadrature.jl,
you get a single predictable interface where many of the arguments are
standardized throughout the various integrator libraries. This can be useful
for benchmarking or for library implementations, since libraries which internally
use a quadrature can easily accept a quadrature method as an argument.
For basic multidimensional quadrature we can construct and solve a QuadratureProblem
:
using Quadrature
f(x,p) = sum(sin.(x))
prob = QuadratureProblem(f,ones(2),3ones(2))
sol = solve(prob,HCubatureJL(),reltol=1e-3,abstol=1e-3)
If we would like to parallelize the computation, we can use the batch interface to compute multiple points at once. For example, here we do allocation-free multithreading with Cubature.jl:
using Quadrature, Cubature, Base.Threads
function f(dx,x,p)
Threads.@threads for i in 1:size(x,2)
dx[i] = sum(sin.(@view(x[:,i])))
end
end
prob = QuadratureProblem(f,ones(2),3ones(2),batch=2)
sol = solve(prob,CubatureJLh(),reltol=1e-3,abstol=1e-3)
If we would like to compare the results against Cuba.jl's Cuhre
method, then
the change is one argument change:
using Cuba
sol = solve(prob,CubaCuhre(),reltol=1e-3,abstol=1e-3)
Quadrature.jl is a fully differentiable quadrature library. Thus it adds the ability to perform automatic differentiation over any of the libraries that it calls. It integrates with ForwardDiff.jl for forward-mode automatic differentiation and Zygote.jl for reverse-mode automatic differentiation. For example:
using Quadrature, ForwardDiff, FiniteDiff, Zygote, Cuba
f(x,p) = sum(sin.(x .* p))
lb = ones(2)
ub = 3ones(2)
p = [1.5,2.0]
function testf(p)
prob = QuadratureProblem(f,lb,ub,p)
sin(solve(prob,CubaCuhre(),reltol=1e-6,abstol=1e-6)[1])
end
dp1 = Zygote.gradient(testf,p)
dp2 = FiniteDiff.finite_difference_gradient(testf,p)
dp3 = ForwardDiff.gradient(testf,p)
dp1[1] ≈ dp2 ≈ dp3
To use this package, you always construct a QuadratureProblem
. This has a
constructor:
QuadratureProblem(f,lb,ub,p=NullParameters();
nout=1, batch = 0, kwargs...)
f
: Either a functionf(x,p)
for out-of-place orf(dx,x,p)
for in place.lb
: Either a number or vector of lower boundsub
: Either a number or vector of upper boundsp
: The parameters associated with the problemnout
: The output size of the functionf
. Defaults to1
, i.e. a scalar integral output.batch
: The preferred number of points to batch. This allows user-side parallelization of the integrand. Ifbatch != 0
, then eachx[:,i]
is a different point of the integral to calculate, and the output should benout x batchsize
. Note thatbatch
is a suggestion for the number of points, and it is not necessarily true thatbatch
is the same asbatchsize
in all algorithms.
Additionally, we can supply iip
like QuadratureProblem{iip}(...)
as true or
false to declare at compile time whether the integrator function is in-place.
The following algorithms are available:
QuadGKJL
: Uses QuadGK.jl. Requiresnout=1
andbatch=0
.HCubatureJL
: Uses HCubature.jl. Requiresbatch=0
.VEGAS
: Uses MonteCarloIntegration.jl. Requiresnout=1
.CubatureJLh
: h-Cubature from Cubature.jl. Requiresusing Cubature
.CubatureJLp
: p-Cubature from Cubature.jl. Requiresusing Cubature
.CubaVegas
: Vegas from Cuba.jl. Requiresusing Cuba
.CubaSUAVE
: SUAVE from Cuba.jl. Requiresusing Cuba
.CubaDivonne
: Divonne from Cuba.jl. Requiresusing Cuba
.CubaCuhre
: Cuhre from Cuba.jl. Requiresusing Cuba
.
reltol
: Relative toleranceabstol
: Absolute tolerancemaxiters
: The maximum number of iterations
Additionally, the extra keyword arguments are splatted to the library calls, so see the documentation of the integrator library for all of the extra details. These extra keyword arguments are not guaranteed to act uniformly.