Fix bugs in jacobian approximation

[mikeingold]

Completed

• Discovered and fixed a bug that reduced accuracy of some integral calculations.
• Discovered and fixed a bug causing Floating point argument FP not working for Float32 #74.
• Improved performance of jacobian function, leading to a performance improvement of about 10% in most integral calculations.
• Unit tests
  • Added combinations.jl tests for Box in (parametric) 1D, 2D, and 3D.
  • Enabled verbose for auto_tests.jl sets, allowing for direct comparison of runtimes versus combinations.jl tests.
  • Removed all Box tests from auto_tests.jl as now obsoleted.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 100.00%. Comparing base (a445f66) to head (9fcaedc).
Report is 2 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff            @@
##              main       #96   +/-   ##
=========================================
  Coverage   100.00%   100.00%           
=========================================
  Files           17        17           
  Lines          269       258   -11     
=========================================
- Hits           269       258   -11     

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

[mikeingold]

Whoa. It looks like these tests uncovered a bug that somehow the auto tests weren't finding. For some reason the Jacobian calculation is using a central finite diff near the interval boundaries. The conditional inside ∂ₙr! is supposed to switch to left or right finite diff when that close to the domain boundaries to prevent this.

[mikeingold]

It looks like I accidentally introduced this bug with commit #8d5f81a in v0.13.5, where I generalized the `jacobian` function to allow `ts::V` where `V<:Union{AbstractVector, Tuple}`. The internal finite-diff functions made use of a re-usable buffer vector that was no longer being re-zero'd between usages. This should have resulted in inaccurate `differential` results, so I'm surprised that the impact was apparently so small that tests still passed (though a few received slightly-increased `atol`/`rtol` values). TODO (maybe in a new branch/PR): - Overhaul jacobian's finite-diff code: eliminate re-usable buffer without new allocations - Re-run prior tests to determine whether explicit/relaxed tolerances are still required

…

On Tue, Oct 1, 2024 at 10:08 AM github-actions[bot] < ***@***.***> wrote: ***@***.***[bot]* commented on this pull request. ------------------------------ In test/combinations.jl <#96 (comment)> : > + + # Vector integrand + vsol = fill(sol, 3) + @test integral(fv, box, GaussLegendre(100))≈vsol rtol=1e-6 + @test integral(fv, box, GaussKronrod()) ≈ vsol + @test integral(fv, box, HAdaptiveCubature()) ≈ vsol + + # Integral aliases + @test lineintegral(f, box) ≈ sol + @test_throws "not supported" surfaceintegral(f, box) + @test_throws "not supported" volumeintegral(f, box) +end + ***@***.*** "Meshes.Box 2D" setup=[Setup] begin + a = π + box = Box(Point(0,0), Point(a,a)) *[JuliaFormatter]* reported by reviewdog <https://github.com/reviewdog/reviewdog> 🐶 ⬇️ Suggested change - box = Box(Point(0,0), Point(a,a)) + box = Box(Point(0, 0), Point(a, a)) — Reply to this email directly, view it on GitHub <#96 (review)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AA657OQXLL32INRTF46YBL3ZZKUEHAVCNFSM6AAAAABPELP43CVHI2DSMVQWIX3LMV43YUDVNRWFEZLROVSXG5CSMV3GSZLXHMZDGNBQGQ4DKMRTGY> . You are receiving this because you authored the thread.Message ID: ***@***.***>

[mikeingold]

Good news! It looks like the current rewrite of jacobian fixes Issue #74. It turns out the issue was caused by the default epsilon value being specified as a Float64. When testing with `FP=BigFloat` this value was being promoted to BigFloat anyway, but when testing with `FP=Float32` this caused the parametric coordinates to be promoted to Float64.

…

On Tue, Oct 1, 2024 at 7:29 PM github-actions[bot] ***@***.***> wrote: ***@***.***[bot]* commented on this pull request. ------------------------------ In test/combinations.jl <#96 (comment)> : > + @test_throws "not supported" volumeintegral(f, box) +end + ***@***.*** "Meshes.Box 3D" setup=[Setup] begin + a = π + box = Box(Point(0, 0, 0), Point(a, a, a)) + + function f(p::P) where {P <: Meshes.Point} + x, y, z = ustrip.((p.coords.x, p.coords.y, p.coords.z)) + (sqrt(a^2 - x^2) + sqrt(a^2 - y^2) + sqrt(a^2 - z^2)) * u"Ω/m^3" + end + fv(p) = fill(f(p), 3) + + # Scalar integrand + sol = 3a^2 * (π * a^2 / 4) * u"Ω" + @test integral(f, box, GaussLegendre(100))≈sol *[JuliaFormatter]* reported by reviewdog <https://github.com/reviewdog/reviewdog> 🐶 ⬇️ Suggested change - @test integral(f, box, GaussLegendre(100))≈sol + @test integral(f, box, GaussLegendre(100)) ≈ sol ------------------------------ In test/combinations.jl <#96 (comment)> : > + + function f(p::P) where {P <: Meshes.Point} + x, y, z = ustrip.((p.coords.x, p.coords.y, p.coords.z)) + (sqrt(a^2 - x^2) + sqrt(a^2 - y^2) + sqrt(a^2 - z^2)) * u"Ω/m^3" + end + fv(p) = fill(f(p), 3) + + # Scalar integrand + sol = 3a^2 * (π * a^2 / 4) * u"Ω" + @test integral(f, box, GaussLegendre(100))≈sol + @test_throws "not supported" integral(f, box, GaussKronrod()) + @test integral(f, box, HAdaptiveCubature()) ≈ sol + + # Vector integrand + vsol = fill(sol, 3) + @test integral(fv, box, GaussLegendre(100))≈vsol *[JuliaFormatter]* reported by reviewdog <https://github.com/reviewdog/reviewdog> 🐶 ⬇️ Suggested change - @test integral(fv, box, GaussLegendre(100))≈vsol + @test integral(fv, box, GaussLegendre(100)) ≈ vsol — Reply to this email directly, view it on GitHub <#96 (review)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AA657OSFPUFHXQHE265NJSTZZMV7PAVCNFSM6AAAAABPELP43CVHI2DSMVQWIX3LMV43YUDVNRWFEZLROVSXG5CSMV3GSZLXHMZDGNBRGYZDCOJQGU> . You are receiving this because you authored the thread.Message ID: ***@***.***>

[mikeingold]

I never cease to be amazed how drastically even just a few allocations can slow down code. I ran some @benchmarks and for a 3D Box the mean execution time for jacobian drops by an order of magnitude with this code change.

[JoshuaLampert]

That's really good news. Thanks a lot @mikeingold. Great work! You say you fixed a bug that lead to reduced accuracy. How big was the loss of accuracy approximately? Can we maybe make some tolerances more strict after this fix?

[mikeingold]

That's really good news. Thanks a lot @mikeingold. Great work! You say you fixed a bug that lead to reduced accuracy. How big was the loss of accuracy approximately? Can we maybe make some tolerances more strict after this fix?

I'm honestly a bit puzzled that this bug didn't cause test failures across the board.

Essentially, for a parametricly-3D geometry, the Jacobian should have been finding tangent-like vectors in the parametric directions [1, 0, 0], [0, 1, 0], [0, 0, 1]... but in my mind the bug would have instead been finding tangents in the directions [1, 0, 0], [1, 1, 0], [1, 1, 1]. Then the differential element gets calculated via the outer products of all of these. Somehow, it seems like the difference was so small that it didn't even fail most tests with default precision.

[mikeingold]

Either way, this is a significantly-improved version of the jacobian function. It performs exactly the same work without all the extra abstraction and without all of the extra allocations. On my machine, the old jacobian runs in a mean time of around 8 us, while the new one is around 800 ns. At the aggregate level this seems to have reduced CI time by roughly 10%.

Longer-term I'd like to re-evaluate all of the explicit rtol values and rules, e.g. GaussLegendre(100) is just a holdover from the auto-tests that seemed to be good enough to universally pass all unit tests. Also, weirdly, the current Box in 2D tests seem to take twice as long as the Box in 3D tests.

[mikeingold]

Did you see the comment I made earlier about anonymous function arguments? I think it got rolled into a resolved review.

[JoshuaLampert]

It's indeed strange that this bug didn't cause more trouble before, but anyway. I made a small suggestion regarding the anonymous function, but it's also not ideal. Otherwise this looks good to me. Thanks again!

[mikeingold]

Once this gets merged, I'll get #98 merged and release a v0.14.1.