Skip to content

Add integrals over geometries#1333

Merged
juliohm merged 24 commits intomasterfrom
integrals
Mar 28, 2026
Merged

Add integrals over geometries#1333
juliohm merged 24 commits intomasterfrom
integrals

Conversation

@juliohm
Copy link
Copy Markdown
Member

@juliohm juliohm commented Mar 17, 2026

Cleaned version of #1327.

@codecov
Copy link
Copy Markdown

codecov bot commented Mar 18, 2026

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 85.68%. Comparing base (28ca52d) to head (aeb5ae0).
⚠️ Report is 5 commits behind head on master.

Additional details and impacted files
@@            Coverage Diff             @@
##           master    #1333      +/-   ##
==========================================
+ Coverage   85.54%   85.68%   +0.13%     
==========================================
  Files         199      200       +1     
  Lines        6384     6439      +55     
==========================================
+ Hits         5461     5517      +56     
+ Misses        923      922       -1     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

@JoshuaLampert
Copy link
Copy Markdown
Member

The tests are now run in CI on Julia v1. With the new version on main of IntegrationInterface.jl many @test_broken were fixed. The remaining ones are

The Float32 tests are failing. This looks like an issue in IntegrationInterface.jl. I'll try to create an MWE and report.

@juliohm
Copy link
Copy Markdown
Member Author

juliohm commented Mar 18, 2026 via email

@JoshuaLampert
Copy link
Copy Markdown
Member

Another thing to consider are Chains like Ring and Rope. While this currently works due to the parametric function for them defined in Meshes.jl, we still might consider to add a specialization for them decomposing them into their Segments and summing the subintegrals because the parametric function for a Chain is not continuously differentiable, so differential might have issues, see also #1155 (comment). I haven't checked the difference though. Would be interesting to see if we get more accurate results when decomposing them into Segments.

@JoshuaLampert
Copy link
Copy Markdown
Member

Also another though/question: I see you increased the number of Gauss-Legendre points to 20 by default (which I find a reasonable choice). We still use some explicit relative tolerances in the tests, which I think came from the time where we had a smaller number of points by default. Did you check if we can relieve the rtol for some tests or completely remove it and use the default tolerances (as we already do in some tests)? I would like to use as strict tolerances as possible to avoid missing some subtle bug.

@juliohm
Copy link
Copy Markdown
Member Author

juliohm commented Mar 18, 2026

we still might consider to add a specialization for them decomposing them into their Segments and summing the subintegrals

Will experiment with this approach.

I see you increased the number of Gauss-Legendre points to 20 by default (which I find a reasonable choice).

I think of this current default as something temporary while II.jl provides a backend that also supports simplices. At the end we will try to pick a single backend that handles boxes and simplices, and will set the default parameters using some concrete criteria.

@JoshuaLampert
Copy link
Copy Markdown
Member

we still might consider to add a specialization for them decomposing them into their Segments and summing the subintegrals

Will experiment with this approach.

Nice! It indeed seems to be much more accurate. The test for Rope gives 7.001791171179793 A m with the previous implementation and 7.000000000000025 A m after your commit 1b36eb3.

"""
function localintegral(fun, geom::Geometry, method=GAUSSLEGENDRE)
# integrand is equal to function times differential element
integrand(uvw...) = fun(uvw...) * differential(geom, uvw)
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it would be good to allow for a custom diff_method that gets passed to differential so you can customize both the integration method and the diff_method for differentiation.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It makes sense. Out of curiosity, do you have an example where users might need to change the default differentiation and integration methods? I am wondering if there is always a "best" default in these geometric settings. If the notion of "best" default exists, i.e., a backend that is the fastest, more accurate, native Julia etc., then we need to reassess the need for customization of backends in these functions.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you want the best in terms of accuracy you should obviously take something AD based because it is exact in contrast to FiniteDiff.jl. However, (as we already saw in the differential PR) AD is easier to fail and ForwardDiff.jl conceptionally simpler. In terms of speed, I think forward mode AD is mostly faster than FD here, but that would need to be benchmarked. So I would like to allow for customization. I think @kylebeggs is also interested in customizable differentiation backends for integral since he added Enzyme.jl support to MeshIntegrals.jl.
Regarding customizability of the integration method, this mostly depends on the geometry. For 1D geometries QuadGK.jl should normally be most efficient, more multi some H-adaptive cubature. On the other hand if you want to compute many integrals, you might want to pre-allocate Gauss-Legendre (or some other Gauss) nodes and weights, so that you can reuse them for the different integrals you want to compute. So I would also like to have the option to customize this.
So to summarize a "best" default meaning best for all use-cases is probably difficult to find for both the differentiation and integral method, but of course, having sensible defaults that are robust and are reasonably fast and accurate in most use cases is still important.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

One thing we need to do in order to support customization of both backends (differentation + integration) is replace the current positional argument by two keyword arguments:

derivative(..., ∂backend=...)
integral(..., ∂backend=..., ∫backend=...)

We could do that in a separate PR after this one here is merged.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, sounds good to me.

@JoshuaLampert
Copy link
Copy Markdown
Member

JoshuaLampert commented Mar 19, 2026

I reran one Float32 test to see if pablosanjose/IntegrationInterface.jl#30 fixed the Float32 tests, see https://github.com/JuliaGeometry/Meshes.jl/actions/runs/23269401309/job/67774013904?pr=1333. The good message is yes, the previous failing Float32 tests are fixed. The bad message is we get a new failing test for BezierCurve, which returns NaN32.

@JoshuaLampert
Copy link
Copy Markdown
Member

Interesting, can you try to use an AD backend for differential?

@juliohm
Copy link
Copy Markdown
Member Author

juliohm commented Mar 20, 2026

Opened a dedicated issue: #1336

@JoshuaLampert
Copy link
Copy Markdown
Member

The Float32 tests are now running for almost an hour, while the Float64 tests finished after ~9 minutes. Seems like they hang. Do you know why that is the case, @pablosanjose? Is it an issue with HAdaptiveIntegration.jl?

@pablosanjose
Copy link
Copy Markdown

That has to be a bug. Perhaps you need a finite atol in this case? Try II.Backend.HAdaptiveIntegration(; atol =1e-7,). Can you isolate the hanging tests for me so I can investigate?

@juliohm
Copy link
Copy Markdown
Member Author

juliohm commented Mar 22, 2026

Perhaps you need a finite atol in this case? Try II.Backend.HAdaptiveIntegration(; atol =1e-7,)

Tried it and even with Float64 the test with Box3D and Hexahedron are hanging.

@pablosanjose
Copy link
Copy Markdown

Can you write here the function and the domain bounds?

@pablosanjose
Copy link
Copy Markdown

This seems to be the problem:

julia> f(x₁, x₂, x₃, a) = ((a^2 - x₁^2) + (a^2 - x₂^2) + (a^2 - x₃^2))
f (generic function with 1 method)

julia> J = Integral(f, Domain.Box{3}(a -> ((0,0,0), (a, a, a))); backend = Backend.HAdaptiveIntegration())
Integral
  Mutating   : false
  Domain     : Functional{Box{3}}
  Backend    : HAdaptiveIntegration
  Integrand  : f

julia> J(0.2)
┌ Warning: maximum number of subdivisions reached `maxsubdiv=65536`, try increasing the keyword argument `maxsubdiv`.
└ @ HAdaptiveIntegration ~/.julia/packages/HAdaptiveIntegration/JRwB2/src/integrate.jl:120
0.0037699120880203977

@pablosanjose
Copy link
Copy Markdown

Opened an issue at HAdaptiveCubature: zmoitier/HAdaptiveIntegration.jl#57

@juliohm
Copy link
Copy Markdown
Member Author

juliohm commented Mar 24, 2026

By setting rtol=1e-3 in the HAdaptiveIntegration backend we manage to pass all tests in Julia latest release (apparently LTS doesn't support [sources]). However, some tests only pass with rtol=1e-2.

@JoshuaLampert
Copy link
Copy Markdown
Member

JoshuaLampert commented Mar 24, 2026

Regarding the failing downstream test: I am seeing the same error on another repo. So I think CairoMakie.jl is currently broken.
Edit: Fixed again, see also https://discourse.julialang.org/t/makie-precompilation-failed-in-linux/136368/5.

# ------------------------------------------------------------------

# default integration method
const HADAPTIVE = II.Backend.HAdaptiveIntegration(rtol=1e-3)
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am wondering whether it makes sense to make rtol depending on the floating point type. The default for this is sqrt(eps(T)), which evaluates to 1.4901161193847656e-8 for Float64 and to 0.00034526698f0 for Float32. So setting this always to 1e-3 means we loose quite some accuracy for the Float64 case. Maybe it would be more sensible to use something like 10 * sqrt(eps(T)) (since sqrt(eps(T)) was too small), which for Float32 is close to 1e-3, but gives us more accuracy for the Float64 case. We can also add that in another PR though because I think we wanted to touch the default anyway and make it a function depending on the geometry, right?

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I played with atol=Meshes.atol(numtype(eltype(geom))) in the II.jl backend, but forgot to play with rtol too. We could investigate this issue in a separate PR to avoid delaying the progress 👍🏽

@pablosanjose
Copy link
Copy Markdown

IntegrationInterface.jl is now registered, you can add it to your deps for CI

@JoshuaLampert
Copy link
Copy Markdown
Member

IntegrationInterface.jl is now registered, you can add it to your deps for CI

Already done :)

@github-actions
Copy link
Copy Markdown
Contributor

github-actions bot commented Mar 28, 2026

Benchmark Results (Julia vlts)

Time benchmarks
master aeb5ae0... master / aeb5ae0...
clipping/SutherlandHodgman 3.68 ± 0.37 μs 3.77 ± 1.6 μs 0.976 ± 0.42
discretization/simplexify 0.663 ± 0.063 ms 0.66 ± 0.078 ms 1 ± 0.15
intersection/ray-triangle 0.05 ± 0.001 μs 0.05 ± 0.01 μs 1 ± 0.2
sideof/ring/large 6.54 ± 0.01 μs 6.53 ± 0.01 μs 1 ± 0.0022
sideof/ring/small 0.05 ± 0.001 μs 0.05 ± 0.001 μs 1 ± 0.028
topology/half-edge 2.71 ± 0.042 ms 2.7 ± 0.052 ms 1 ± 0.025
winding/mesh 16.1 ± 0.19 ms 16 ± 0.12 ms 1.01 ± 0.014
time_to_load 1.15 ± 0.0097 s 1.16 ± 0.0096 s 0.992 ± 0.012
Memory benchmarks
master aeb5ae0... master / aeb5ae0...
clipping/SutherlandHodgman 0.053 k allocs: 4.97 kB 0.053 k allocs: 4.97 kB 1
discretization/simplexify 18.1 k allocs: 1.92 MB 18.1 k allocs: 1.92 MB 1
intersection/ray-triangle 0 allocs: 0 B 0 allocs: 0 B
sideof/ring/large 0 allocs: 0 B 0 allocs: 0 B
sideof/ring/small 0 allocs: 0 B 0 allocs: 0 B
topology/half-edge 18.1 k allocs: 2.92 MB 18.1 k allocs: 2.92 MB 1
winding/mesh 23.2 k allocs: 3.08 MB 23.2 k allocs: 3.08 MB 1
time_to_load 0.153 k allocs: 14.5 kB 0.153 k allocs: 14.5 kB 1

@github-actions
Copy link
Copy Markdown
Contributor

github-actions bot commented Mar 28, 2026

Benchmark Results (Julia v1)

Time benchmarks
master aeb5ae0... master / aeb5ae0...
clipping/SutherlandHodgman 2.69 ± 0.28 μs 2.68 ± 0.25 μs 1 ± 0.14
discretization/simplexify 0.384 ± 0.018 ms 0.384 ± 0.021 ms 0.999 ± 0.073
intersection/ray-triangle 0.06 ± 0 μs 0.051 ± 0.01 μs 1.18 ± 0.23
sideof/ring/large 6.84 ± 0.019 μs 6.84 ± 0.06 μs 1 ± 0.0092
sideof/ring/small 0.06 ± 0.001 μs 0.06 ± 0 μs 1 ± 0.017
topology/half-edge 2.73 ± 0.24 ms 2.72 ± 0.21 ms 1 ± 0.12
winding/mesh 15.4 ± 0.25 ms 15.4 ± 0.33 ms 1 ± 0.027
time_to_load 1.05 ± 0.021 s 1.04 ± 0.0086 s 1.01 ± 0.022
Memory benchmarks
master aeb5ae0... master / aeb5ae0...
clipping/SutherlandHodgman 0.064 k allocs: 5.55 kB 0.064 k allocs: 5.55 kB 1
discretization/simplexify 0.0362 M allocs: 1.93 MB 0.0362 M allocs: 1.93 MB 1
intersection/ray-triangle 0 allocs: 0 B 0 allocs: 0 B
sideof/ring/large 0 allocs: 0 B 0 allocs: 0 B
sideof/ring/small 0 allocs: 0 B 0 allocs: 0 B
topology/half-edge 25.9 k allocs: 3.17 MB 25.9 k allocs: 3.17 MB 1
winding/mesh 0.0413 M allocs: 3.08 MB 0.0413 M allocs: 3.08 MB 1
time_to_load 0.145 k allocs: 11 kB 0.145 k allocs: 11 kB 1

Copy link
Copy Markdown
Member

@JoshuaLampert JoshuaLampert left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

LGTM and CI is green 🥳
The remaining three open conversations will be handled in separate PRs.

@juliohm juliohm marked this pull request as ready for review March 28, 2026 14:01
@juliohm
Copy link
Copy Markdown
Member Author

juliohm commented Mar 28, 2026

Amazing stuff! 🎉 Thank you all for the huge help. Special thanks to @pablosanjose for the new IntegrationInterface.jl package, super nice addition!

@juliohm juliohm merged commit 39a998c into master Mar 28, 2026
19 checks passed
@juliohm juliohm deleted the integrals branch March 28, 2026 14:04
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants