[ANN] Trapz 2.0 - Fast integration routines for discrete domains

Hi everyone, this package is not exactly new, but i never announced it, and recently i had the time&need to make it even faster and flexible

a copy of readme.md:

Trapz.jl

A simple Julia package to perform trapezoidal integration over common Julia arrays.

Docs
Build Status
codecov

the package is now registered on Julia Registry, so it can be added as follows


import Pkg
Pkg.pkg"add Trapz"

Example Usage:

using Trapz

vx=range(0,1,length=100)
vy=range(0,2,length=200)
vz=range(0,3,length=300)
M=[x^2+y^2+z^2 for x=vx,y=vy,z=vz]
I=trapz((vx,vy,vz),M)
print("result: ",I)
result: 28.000303707970264

Benchmarks

using BenchmarkTools
@btime trapz($(vx,vy,vz),$M);
2.793 ms (4 allocations: 157.30 KiB)
@btime trapz($(:,vy, vz),$M);
3.046 ms (3 allocations: 157.20 KiB)
@btime trapz($(:,vy,:),$M);
4.022 ms (2 allocations: 234.45 KiB)

Comparison to Numpy

using PyCall
np=pyimport("numpy")

timenumpy = @belapsed np.trapz(np.trapz(np.trapz($M,$vz),$vy),$vx)
timejulia = @belapsed trapz($(vx,vy,vz),$M)

how_faster=timenumpy/timejulia

print("Trapz.jl is ~ ",how_faster," times faster than numpy's trapz")
Trapz.jl is ~ 7.2802292929009 times faster than numpy's trapz
13 Likes

How fast is this relative to quadgk?

2 Likes

They have entirely different usecase, QuadGK needs a function, this operates on arrays. comparing their speed does not make much sense. in any case, quadgk can be pretty heavy on resources if you are not careful.

5 Likes

Nice! I like the new multidimensional integration syntax.

1 Like

If it is used appropriately, Gaussian quadrature can require exponentially fewer integrand evaluations than the trapezoidal rule as the required accuracy increases, assuming smooth functions. quadgk (and, in multiple dimensions, the HCubature package) are also adaptive (they put more function evaluations where the function is more oscillatory etc.).

As @francesco.alemanno says, the main utility of the trapezoidal rule is for cases where you cannot evaluate your function at arbitrary points in the domain, and only have a predetermined grid of points to work with.

6 Likes

Your PR to Trapz is to be thanked, i made some enhancements (like avoiding PermutedDimsArray) but it is almost like you designed it, thank you again @lstagner :slight_smile:

My experience is that for non-smooth/spiky functions the advantage is also enormous. Some functions seem more or less impossible to integrate accurately on a regular grid.