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

Package Announcements
1

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
17 Likes
2

How fast is this relative to quadgk?

3 Likes
3

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.

8 Likes
4

Nice! I like the new multidimensional integration syntax.

1 Like
5

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.

8 Likes
6

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:

7

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.