Integration of functions over geometries

GitHub - JuliaApproximation/FastGaussQuadrature.jl: Julia package for Gaussian quadrature?

Amazing links. Thank you all. The plan is to implement some of these for Meshes.jl geometries.

3 Likes

@Perrin_Meyer is quadpy closed-source? Can’t find their source code to understand what is going on.

@Oscar_Smith I think FastGaussQuadrature.jl doesn’t implement integration over geometries. It is mostly hypercubes as far as I know.

This is a very good C++ library for integration over convex polytopes that I used some time ago. No Julia interface as far as I know, however…

1 Like

You can take a quadrangle or a hexahedron and easily map into into a hypercube by a coordinate transformation, at which point you can apply any quadrature rule for a hypercube (with a corresponding Jacobian factor), whether it is a simple first-order bilinear rule or a tensor product of Gaussian quadrature rules or an adaptive quadrature scheme.

More generally, things like triangular domains can also be converted to rectangles with a change of variables, and general polygons could be decomposed into triangles, etcetera.

2 Likes

Weird, because last time I looked (six months ago) all the python code was right on github, and I thought it was all open source… My best guess is this is some sort of code-reorganization?

Yes, it is unfortunate that the source code is not open, it would be a great resource to understand the tradeoff of different methods and designs.

@juliohm Did you ever come around to making an integration scheme work on Meshes.jl? I can’t seem to find anything…

@RomeoV we added some basic functionality in the GeoStats.jl modules on top of Meshes.jl. I think there is a integrate utility function in GeoStatsBase.jl if I remember correctly.

There is also some ongoing work by @mikeingold on line integrals with unitful coordinates over 1D polytopes from Meshes.jl: https://github.com/mikeingold/LineIntegrals.jl

We are in the process of generalizing the coordinate infrastructure in Meshes.jl to accommodate units and non Euclidean coordinates.

1 Like

As-is anyone interested in playing around with it can Pkg.add LineIntegals directly from GitHub. The core functionality is to link up Meshes.jl with QuadGK.jl to enable line/path/contour integrals where the geometry is specified using Meshes, e.g.:

using Meshes
using QuadGK
using LineIntegrals

# Construct a path that approximates a unit circle on the xy-plane
#   embedded in a 3D coordinate system
unit_circle = BezierCurve(
    [Point(cos(t), sin(t), 0.0) for t in range(0,2pi,length=361)]
)

# Real function
f(x, y, z) = abs(x + y)
f(p::Point) = f(p.coords...)
quadgk(f, unit_circle)       # -> (5.551055333711397, 1.1102230246251565e-16)

The implementation for 1D polytopes is actually very simple since Meshes already defines a parametric function (::Segment)(t) -> ::Point. The new methods just define a parameterized 1D integral from 0 to 1, and as long as you have a method f(::Point) defined it just works.

I’ve done some testing and everything seems to also compose well with dimensionful types like those from Unitful.jl and DynamicQuantities.jl so you could also have, say, a f that outputs values in units Ohms/meter, integrate f along some Meshes geometry with dimensionful coordinates and the quadgk result will be dimensionful in units Ohms.

I’d considered making this a standalone utility package for General, but am now leaning toward just trying to upstream the useful bits into a Meshes extension. Defining QuadGK.quadgk(f, ::Meshes.*) from a third-party package is type-piracy; I could technically work around this by just defining some LineIntegrals.integral(f, ::Meshes.*) but I don’t know that I’d be adding any real value with a thin rebranded wrapper. Probably preferable to just have a Meshes extension defined for QuadGK that enables the example above without a third package dependency.

1 Like

Thanks @mike.ingold for leading this effort. Have you considered

as an alternative to QuadGK.jl? It is native Julia and is super fast from what I remember.

I’ve heard of it but haven’t yet had a chance to try it out.

I’ve tended towards QuadGK for 1D problems in the past largely because it’s native Julia that seems to compose well with arbitrary f output types (basically anything that supports (+,-,*,norm) operations) and “just works” for arbitrary atol or rtol specs. The downside I’ve noticed with it, and other adaptive quadrature/cubature implementations, is that if you treat them like a magic bullet and happen to have an f with certain behaviors (like a step change) then solution times will go through the roof.

I see that the FastGaussQuadrature.jl docs promise “16-digit accuracy and in O(n) time”, which seems pretty promising. It looks like all of the built-in quadratures operate on fixed domains, so I’d probably need to implement a change of variables internally to map t instead to [-1,1].

I’ll definitely check it out. Thanks!

1 Like

Also, thanks to @RomeoV for bumping this thread I didn’t realize existed. Longer term, I’m also very interested in having a way to compute surface integrals over arbitrary meshes and 2D geometries. I’ll definitely be reading further into the references SGJ provided above.

1 Like

I did some more research into FastGaussQuadrature.jl yesterday; the paper they reference was a nice read to help me get up to speed. It looks like the implicit caveat to the “16-digit accuracy and in O(n) time” bumper sticker is that this is only precisely true for functions that can be exactly represented by a polynomial of order 2n-1 over the specified interval. Off the top of my head, I’d guess that for functions with periodic content you’d want to configure n to a bare minimum of the number of expected periods per line/geometric feature, i.e, n>l/\lambda.

One of the features I really like about adaptive quadratures like QuadGK.jl is that a quadgk call can be made without any explicit keyword arguments and the solver will produce a result with some “guaranteed” error tolerance, and also explicitly provide an estimate of the actual resultant error tolerance. Granted, from what I’ve now learned, it sounds like these error estimates are by default predicated on the assumption that the given f is well-approximated over the specified interval by a polynomial of order 15 (order=7).

Within LineIntegrals.jl I originally had methods integral(f, geometry) which were really just wrappers over quadgk that returned only the integral value. That seemed like a weak abstraction since it was really just passing through kwargs and only picking off the result[1], so I re-implemented those as simple QuadGK.quadgk methods.

That being said, I definitely think a case could be made that implementing an interface integral(f, geometry) that internally applies tools from FastGaussQuadrature.jl is a non-trivial value-add. I only had a few minutes last night, but I made a quick-and-dirty integral(f, ::Segment) method that works for scalar functions; I’m going to take another pass to generalize it and make sure I’m minimizing allocations (maybe later implement an in-place integrate! version). The only thing I’m not sure about is how to set a default order n. I’m leaning towards picking a reasonable/conservative default like other packages do (7?) and just making sure that the decision criteria is clear in the docstring.

1 Like

That is very exciting @mike.ingold. This approach where you define the high-level interface as integral(f, g) and workout the details with different backends is great. I would suggest that backend options are passed as kwargs so that users can fine-tune the results in specific applications. Something like integral(f, g, order = 7).

If QuadGK.jl has comparable computational performance and better numerical performance, then it makes sense to use it as the default Julia native backend.

Looking forward to seeing this effort evolve into a whole set of integrals over Meshes.jl geometries :100: All our geometries are already parametric in the sense that you can g(u, v) and g(u, v, w) with unitless coordinates u, v, w in [0,1].

1 Like

No, that’s not the case. For polynomials Gaussian quadrature is exact, but for analytic functions Gaussian quadrature converges exponentially fast with the number of quadrature points n, so you can indeed get to nearly machine precision (nearly 16 digits) rapidly, and if you need to go to large n then FastGaussQuadrature.jl computes the quadrature rule in O(n) time.

A downside of FastGaussQuadrature.jl is that it just gives you an n-point rule, it doesn’t tell you what n should be for a given tolerance and a given function, and the rule is not a “nested” rule — it doesn’t give you an error estimate. So you are only your own in finding the correct n for your problem, e.g. by repeatedly doubling n until it converges. (Repeatedly increasing the order is also called p-adaptive quadrature.)

QuadGK is h-adaptive, in that it repeatedly subdivides the integration interval until a tolerance is achieved, and it uses an embedded rule to get an error estimate. This is convenient, and it is also fairly robust if your integrand is badly behaved in some localized region because it will only refine the quadrature around that region.

If your function is analytic, then using high-order quadrature is theoretically more efficient, but on the other hand you rarely have to go to very high order in practice, and the time for the computation of the rule is less important if you are doing integration over and over for many different functions and/or domains. So for smooth functions you can simply pass order=n for a larger n to QuadGK.jl to use a higher-order rule (which is computed once and cached). The computation of the rule is O(n^2), so it will be much slower than FastGaussQuadrature.jl for large n, but since in practice you rarely need n larger than a few hundred it doesn’t matter too much in my experience.

The other big advantage of FastGaussQuadrature.jl is that it has built-in rules for weighted quadrature with many different weight functions, which allow you to efficiently handle known endpoint singularities. See the discussion of weighted Gauss rules in the QuadGK manual, and integrands with singularities and discontinuities.

PS. Although QuadGK is not asymptotically exponentially convergent for a fixed order, unlike repeatedly increasing the order, it often gives the illusion of exponential convergence up to a given precision. See also this post: (Insanely-) Higher Order Derivatives - #17 by stevengj

5 Likes

It’s for exactly these reasons that I generally tend toward h-adaptive quadrature/cubature methods. These days I find myself writing a lot of tooling that will be used to run calculations/models, moreso than direct data-in-hand problems. That being the case, I find it hugely beneficial to have an interface where I can specify some acceptable atol or rtol and get a result within those specs reliably since I can’t know at time of writing what function/signal will even be integrated.

Updates on LineIntegrals:

  • I’ve re-implemented all of the existing integral methods using FGQ.
  • All QuadGK.quadgk(f, <:Meshes.Geometry) methods remain defined. Again, this is type piracy so would probably need to be PR’d into another package.
  • I’m planning to implement an internal function with in-place memory re-use for geometries with multiple segments.
  • The README has a toy example integrating f(x,y) = abs(x+y) over a Bezier curve approximation of the unit circle using both techniques. The integrate methods default to n=100 for Gauss-Legendre nodes, and quadgk uses all default kwargs. Both are solved in similar time with similar allocations, but QuadGK.jl seems to be more accurate. Increasing n with integral causes its result to converge toward quadgk’s, but at the expense of a significant increase in allocations. Not a scientific test by any means; your mileage may vary.
2 Likes

Hi,

I’m new to the thread, so I may have missed some details, and I wanted to share a package that allows integration over polytopes using nested integration as an algorithm, which is convenient for integration in low-dimensions (up to 3 or 4), and even ideal for nearly singular integrands. In two (and more) dimensions the nested integration algorithm evaluates the multivariate integral in terms of recursive 1d integrals, i.e:

\int_{a_x}^{b_x} dx\ \int_{a_y(x)}^{b_y(x)} dy\ f(x,y) = \int_{a_x}^{b_x} dx\ I_y(x), \quad I_y(x) = \int_{a_y(x)}^{b_y(x)} dy\ f(x,y)

What is nice about this algorithm is that it is compatible with any quadrature scheme and any geometry (polytopes, manifolds) as long as either a parametrization of the geometry is available, or optimization on that geometry is possible and fast compared to evaluation of the integrand.
The example below is with the package AutoBZCore.jl, which emulates the Integrals.jl interface, and it allows for integration over polytopes from Polyhedra.jl

julia> using AutoBZCore, Polyhedra

# make convex hull in 3d
julia> domain = load_limits(polyhedron(vrep(rand(10, 3))))
PolyhedraExt.PolyhedralLimits{3, Float64, DefaultPolyhedron{Float64, MixedMatHRep{Float64, Matrix{Float64}}, MixedMatVRep{Float64, Matrix{Float64}}}}(convexhull([0.900543, 0.738892, 0.00729183], [0.165838, 0.909174, 0.316936], [0.552774, 0.774557, 0.908199], [0.906346, 0.161918, 0.456898], [0.979034, 0.303522, 0.0330761], [0.146472, 0.13728, 0.183881], [0.673617, 0.694496, 0.738391], [0.380986, 0.0822024, 0.346046], [0.894882, 0.748271, 0.0540279], [0.40436, 0.03022, 0.151078]))

# setup integrand
julia> prob = AutoBZCore.IntegralProblem((x,p) -> sin(prod(x)), domain)
AutoBZCore.IntegralProblem{var"#5#6", PolyhedraExt.PolyhedralLimits{3, Float64, DefaultPolyhedron{Float64, MixedMatHRep{Float64, Matrix{Float64}}, MixedMatVRep{Float64, Matrix{Float64}}}}, AutoBZCore.NullParameters}(var"#5#6"(), PolyhedraExt.PolyhedralLimits{3, Float64, DefaultPolyhedron{Float64, MixedMatHRep{Float64, Matrix{Float64}}, MixedMatVRep{Float64, Matrix{Float64}}}}(convexhull([0.900543, 0.738892, 0.00729183], [0.165838, 0.909174, 0.316936], [0.552774, 0.774557, 0.908199], [0.906346, 0.161918, 0.456898], [0.979034, 0.303522, 0.0330761], [0.146472, 0.13728, 0.183881], [0.673617, 0.694496, 0.738391], [0.380986, 0.0822024, 0.346046], [0.894882, 0.748271, 0.0540279], [0.40436, 0.03022, 0.151078])), AutoBZCore.NullParameters())

# call quadgk recursively
julia> AutoBZCore.solve(prob, NestedQuad(QuadGKJL()))
AutoBZCore.IntegralSolution{Float64, Float64}(0.019284903620528718, 1.7385351242394514e-18, true, -1)

julia> using FastGaussQuadrature

# use a product rule of 10-point Gauss-Legendre quadratures
julia> AutoBZCore.solve(prob, NestedQuad(QuadratureFunction(fun=gausslegendre, npt=10)))
AutoBZCore.IntegralSolution{Float64, Nothing}(0.01928490362052873, nothing, true, -1)

# pick quadgk for two variables and Gauss-Legendre for the last
julia> AutoBZCore.solve(prob, NestedQuad(QuadGKJL(), QuadratureFunction(fun=gausslegendre, npt=10), QuadGKJL()))
AutoBZCore.IntegralSolution{Float64, Float64}(0.01928490362052872, 5.454892180317694e-19, true, -1)

My hope is to integrate this functionality with Integrals.jl (no pun intended) when I find the time, and the main idea of the implementation is in
IteratedIntegration.jl.

2 Likes