Add integral and localintegral

[juliohm]

Ports integral from MeshIntegrals.jl and adds localintegral for integration of functions that are defined locally in terms of parametric coordinates.

@JoshuaLampert do you have suggestions before I start writing tests?

[github-actions]

Benchmark Results (Julia vlts)

Time benchmarks

	master	bcf0275...	master / bcf0275...
clipping/SutherlandHodgman	3.64 ± 0.36 μs	3.57 ± 0.26 μs	1.02 ± 0.12
discretization/simplexify	0.671 ± 0.054 ms	0.678 ± 0.066 ms	0.989 ± 0.13
intersection/ray-triangle	0.044 ± 0.003 μs	0.049 ± 0.002 μs	0.898 ± 0.071
sideof/ring/large	5.25 ± 0.006 μs	5.26 ± 0.005 μs	0.999 ± 0.0015
sideof/ring/small	0.047 ± 0.002 μs	0.046 ± 0.001 μs	1.02 ± 0.049
topology/half-edge	4.46 ± 0.051 ms	4.46 ± 0.051 ms	0.999 ± 0.016
winding/mesh	14.5 ± 0.21 ms	14.5 ± 0.22 ms	1 ± 0.021
time_to_load	1.14 ± 0.023 s	1.15 ± 0.011 s	0.985 ± 0.022

Memory benchmarks

	master	bcf0275...	master / bcf0275...
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]

Benchmark Results (Julia v1)

Time benchmarks

	master	bcf0275...	master / bcf0275...
clipping/SutherlandHodgman	2.71 ± 0.25 μs	2.67 ± 0.24 μs	1.02 ± 0.13
discretization/simplexify	0.388 ± 0.019 ms	0.383 ± 0.017 ms	1.01 ± 0.066
intersection/ray-triangle	0.05 ± 0.01 μs	0.05 ± 0.01 μs	1 ± 0.28
sideof/ring/large	6.7 ± 0.011 μs	6.85 ± 0.06 μs	0.978 ± 0.0087
sideof/ring/small	0.06 ± 0 μs	0.06 ± 0 μs	1 ± 0
topology/half-edge	2.72 ± 0.17 ms	2.77 ± 0.27 ms	0.981 ± 0.11
winding/mesh	15.4 ± 0.12 ms	15.3 ± 0.18 ms	1 ± 0.014
time_to_load	1.05 ± 0.0053 s	1.1 ± 0.0073 s	0.955 ± 0.008

Memory benchmarks

	master	bcf0275...	master / bcf0275...
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

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 85.62%. Comparing base (0e6a436) to head (bcf0275).
⚠️ Report is 1 commits behind head on master.

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1327      +/-   ##
==========================================
+ Coverage   85.50%   85.62%   +0.12%     
==========================================
  Files         199      200       +1     
  Lines        6376     6409      +33     
==========================================
+ Hits         5452     5488      +36     
+ Misses        924      921       -3     

☔ 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.

[juliohm]

Regarding support for multiple quadrature libraries, we need to think if it is really necessary. If it is necessary, we could start by creating an IntegrationInterface.jl package after this PR is merged, similar to the DifferentiationInterface.jl approach.

[JoshuaLampert]

If it is necessary, we could start by creating an IntegrationInterface.jl package after this PR is merged, similar to the DifferentiationInterface.jl approach.

There is Integrals.jl, which basically provides exactly that.

[JoshuaLampert]

Regarding support for multiple quadrature libraries, we need to think if it is really necessary.

I think it would be nice because Gauss-Kronrod should be more efficient to get a certain accuracy in 1D and H-adaptive cubature should be better suited in multi D.

[juliohm]

There is Integrals.jl, which basically provides exactly that.

I saw it, but it is not at the same level of DifferentiationInterface.jl. It includes a bunch of dependencies, including backends and SciML stuff.

[juliohm]

Please feel free to add tests if you have some already written in MeshIntegrals.jl.

I will try to continue the work over the weekend, and next week we should have a final version to review 👍🏽

[JoshuaLampert]

I have two questions:

• How do we deal with the specializations? Do we add them in (a) later PR(s)?
• How should we write tests? In MeshIntegrals.jl we have quite extensive tests, which take some time though.
  • Should we have tests for every (supported) geometry? I think we should because just because the integral works fine for one or a few geometries does not mean it works for all, see the need for the specializations. I think we can use the non-trivial (i.e., we do not just integrate a constant function) analytical tests we have in the test suite of MeshIntegrals.jl.
  • Should we use some test infrastructure to automate the testing and reduce code duplication and redundancy like in MeshIntegrals, see here or should we rather go with a straightforward approach and testing each geometry individually?

[juliohm]

• How do we deal with the specializations? Do we add them in (a) later PR(s)?

We should add them here as well. I will try to find the time tomorrow to finish up the additions 👍🏽

• Should we have tests for every (supported) geometry?

Yes. We should have full coverage.

• Should we use some test infrastructure to automate the testing and reduce code duplication and redundancy like in MeshIntegrals, see here or should we rather go with a straightforward approach and testing each geometry individually?

I think we can be more straightforward testing each geometry individually. It gives us more freedom to add more tests later without affecting all other tests.

[juliohm]

@JoshuaLampert I've ported some tests from MeshIntegrals.jl to see if everything was working as expected, but the results are not matching. Do you have any idea of what might be happening? Could it be related to the default backend? Should we adjust the tests?

[JoshuaLampert]

@JoshuaLampert I've ported some tests from MeshIntegrals.jl to see if everything was working as expected, but the results are not matching. Do you have any idea of what might be happening? Could it be related to the default backend? Should we adjust the tests?

They are very likely failing due to the low number of Gauss-Legendre points you are using by default, i.e., related to your comment #1327 (comment). I think for general non-linear functions, we should provide a much higher value for order than 1 or 2. In the tests of MeshIntegrals.jl we use 100. Otherwise we can have big errors, which is normal and expected. Maybe 100 is a bit too much, but 50 maybe is enough?

[JoshuaLampert]

Also I just see from the failing Float32 tests, e.g., here that integral does not preserve the floating point type, i.e., the result is in Float64 even though the input is Float32. I think we want the same datatype as output as we had for the geometry, right?
Edit: Ah, I now remember that we had a kwarg FP for this in MeshIntegrals.jl. Maybe this is not strictly necessary to add now.

[juliohm]

Maybe 100 is a bit too much, but 50 maybe is enough?

Got it! I will try to find a small enough number that gives us a good default for moderately nonlinear functions.

I think we want the same datatype as output as we had for the geometry, right?

We want to preserve the type of the input. I thought the implementation would guarantee that, but maybe the backend is promoting the arguments? I will double check.

[JoshuaLampert]

I think we want the same datatype as output as we had for the geometry, right?

We want to preserve the type of the input. I thought the implementation would guarantee that, but maybe the backend is promoting the arguments? I will double check.

I think the problem is that the Gauss-Legendre nodes and weights provided by FastGaussQuadrature.jl are Float64. We would need to convert them to the datatype we want before performing the quadrature.

[juliohm]

I've increased the default order to 100 to match what you had in MeshIntegrals.jl, but the tests are still failing. I think 100 is too much, but we need to figure out why the tests are not passing.

After we sort out the source of the issue, we could pick something like 3 by default (exact approximation of polynomials of degree 5). It is hard to imagine applications where this default wouldn't be enough.

[juliohm]

Could the change in results be related to the change of variable formula? I used t -> (t+1)/2 but MeshIntegrals.jl was using rationals here.

[JoshuaLampert]

Could the change in results be related to the change of variable formula? I used t -> (t+1)/2 but MeshIntegrals.jl was using rationals here.

Hm, currently, I don't see another semantic difference between the two implementations. So give it a try?

[JoshuaLampert]

What is different is the finite difference scheme used to compute the differential. This could of course lead to numerical differences to the previous version, but I would have hoped that the scheme we have here based on FiniteDifferences.jl is a bit more accurate than the custom one from MeshIntegrals.jl.
Oh, now that I am writing this I see, we use a custom rtol in isapprox, see here and here. Maybe try that.

[juliohm]

we use a custom rtol in isapprox

That was it. I will adjust the tests to be more explicit about the n and rtol used.

[juliohm]

@JoshuaLampert I've copied/pasted the tests with infinite geometries from MeshIntegrals.jl, but they are failing too. Perhaps I messed up the change of variables. Can you spot any issue with the current implementation? Did I use the correct change of variable?

I've marked the tests with @test_broken (Ray, Line, Plane).

[JoshuaLampert]

On a brief look, I did not spot anything wrong with the infinite geometries. The transformations look correct to me. I would need to dig deeper and debug together with MeshIntegrals.jl to see where they first differ. I don't have the right now though.

[juliohm]

@JoshuaLampert we have all methods in place. The tests for Ray, Line, Plane, Triangle and Tetrahedron are failing. They are precisely the ones that require custom change of variable. Any idea of what we are missing?

[JoshuaLampert]

Thanks! This looks mostly good to me. I left a few minor suggestions/questions below.

• A test for SimpleMesh is still missing.
• Should we also add a test that the floating point type is preserved by integral or do you think this is already implicitly tested for by the accuracy tests?
• I don't see right now why the geometries, which need a custom transformation fail. It seems like they have the same reason why they fail. So to debug this I would pick the simplest one (I guess Ray), add some debug statements at different steps in Meshes.jl and MeshIntegrals.jl and see where they differ to see what might be the problem. Did you already try something like this? I do not have the time to do that right now though (I am on vacation). Otherwise maybe AI also has an idea?

[JoshuaLampert]

Why does this specialization improve performance?

[juliohm]

We had a method for general Polygon that got erased. The method for general polygon relied on discretization, and this method for Quadrangle relies on direct sampling with the parametric function. I will review the code to see if this method is still necessary.

[juliohm]

• A test for SimpleMesh is still missing.

Will add it right after we fix Triangle and Tetrahedron

• Should we also add a test that the floating point type is preserved by integral or do you think this is already implicitly tested for by the accuracy tests?

I think we can assume the implementation preserves the number type.

• I don't see right now why the geometries, which need a custom transformation fail. It seems like they have the same reason why they fail. So to debug this I would pick the simplest one (I guess Ray), add some debug statements at different steps in Meshes.jl and MeshIntegrals.jl and see where they differ to see what might be the problem.

Will check that.