Add Lavery splines (#23)

[vincmarks]

• 1D Lavery splines
• 2D Lavery splines

[andrewwinters5000]

@vincmarks Would it be okay for me to push directly to your fork? I have done some optimization of the 1D implementation and prototyped the 2D version.

[vincmarks]

Sure, thanks.

[andrewwinters5000]

The 2D Lavery spline is quite slow for larger data sets. There are likely some easy optimizations that can be done when constructing the JuMP model. In 1D, I was able to speed things up by approximately 50% just by restructuring and introducing a new struct.

Below are some benchmark numbers for the one dimensional Lavery spline. Note, I changed the default value of intergral_steps to be much smaller as the original default value of 100 was overkill.

Original implementation:

julia> @benchmark LaverySpline1D($xDataMedium, $yDataMedium)
BenchmarkTools.Trial: 6 samples with 1 evaluation per sample.
 Range (min … max):  901.768 ms …   1.043 s  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     928.236 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   952.323 ms ± 59.910 ms  ┊ GC (mean ± σ):  0.67% ± 1.08%

  █     ▁         ▁                            ▁             ▁
  █▁▁▁▁▁█▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
  902 ms          Histogram: frequency by time          1.04 s <

 Memory estimate: 35.22 MiB, allocs estimate: 546856.

New implementation

julia> @benchmark LaverySpline1D($xDataMedium, $yDataMedium; integral_steps=5)
BenchmarkTools.Trial: 227 samples with 1 evaluation per sample.
Range (min … max):  19.462 ms … 55.591 ms  ┊ GC (min … max): 0.00% … 42.65%
Time  (median):     21.409 ms              ┊ GC (median):    0.00%
Time  (mean ± σ):   22.023 ms ±  2.971 ms  ┊ GC (mean ± σ):  0.47% ±  2.83%

  ▅█▂ ▁   ▃█▆▂ ▁▂    ▁
 ▇███▆█▄▆▆███████▅▆▇▆█▅▆▆▄▄▃▄▄▃▃▁▁▅▄▃▄▁▃▃▃▁▁▁▁▄▁▁▁▃▁▃▄▃▁▁▁▁▃ ▄
 19.5 ms         Histogram: frequency by time        29.1 ms <

Memory estimate: 944.72 KiB, allocs estimate: 16703.

[andrewwinters5000]

Also, here is a figure of the new 2D Lavery example.

[coveralls]

Pull Request Test Coverage Report for Build 21151938608

Details

• 217 of 220 (98.64%) changed or added relevant lines in 4 files are covered.
• No unchanged relevant lines lost coverage.
• Overall coverage increased (+0.7%) to 97.159%

---

Changes Missing Coverage	Covered Lines	Changed/Added Lines	%
src/2D/spline_cache_2D.jl	79	80	98.75%
src/1D/spline_cache_1D.jl	41	43	95.35%

Totals
Change from base Build 20641667137:	0.7%
Covered Lines:	684
Relevant Lines:	704

---

💛 - Coveralls

[codecov-commenter]

Codecov Report

❌ Patch coverage is 50.00000% with 2 lines in your changes missing coverage. Please review.
✅ Project coverage is 96.10%. Comparing base (1c361a6) to head (4236bdf).
⚠️ Report is 1 commits behind head on main.

Files with missing lines	Patch %	Lines
src/TrixiBottomTopography.jl	50.00%	2 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main      #83      +/-   ##
==========================================
- Coverage   96.48%   96.10%   -0.39%     
==========================================
  Files           9        9              
  Lines         484      488       +4     
==========================================
+ Hits          467      469       +2     
- Misses         17       19       +2     

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

Nice! Does it make sense to put the code in an extension for JuMP and HiGHS because they are not so lightweight and not necessary for every user of TrixiBottomTopography.jl?

[andrewwinters5000]

Nice! Does it make sense to put the code in an extension for JuMP and HiGHS because they are not so lightweight and not necessary for every user of TrixiBottomTopography.jl?

That is probably a better way to add this feature that compartmentalizes the implementation. I am currently looking into improved ways to construct the shape preserving splines. Right now the computational cost of the 2D version makes it not very useful for practitioners. However, I think we can remove some of the optimzation constraints...but I need to read more.

[andrewwinters5000]

I dug into the literature a bit and implemented a relaxed shape preserving spline in 2D. The relaxation of shape preservation is global versus point-wise. The best comparison I can think of is switching from TVB to TVD for shock capturing.

It is much more computationally efficient. For the final example in the thesis from Logan and Oscar (the Seaside, Oregon topography with 351 x 351 nodes) their implementation took approximately 8 hours on their laptop to produce a solution. The new implementation I put in run on the same data set runs in approximately 16 seconds. There are still some allocations that could be hunted down, but the Lavery 2D splines are much more useful for practical data.

I also think it is a good idea to have this as a package extension as suggested by @JoshuaLampert . But I want to wait and do that last, once the implementation is cleaned up.

[JoshuaLampert]

Nice! Maybe also add a bullet point in README.md and add it to docs/make.jl and docs/src/reference.md, such that the docstrings are rendered in the documentation?

[andrewwinters5000]

Maybe also add a bullet point in README.md and add it to docs/make.jl and docs/src/reference.md, such that the docstrings are rendered in the documentation?

Good point! I thought the docstrings were rendering but need to check this. I opened an issue already that someone (i.e. me) should add some examples that demonstrate how to use your RBF interpolation and the new Lavery splines.

[JoshuaLampert]

Thanks! This looks mostly good to me. I only have one question below. I cannot really judge the correctness of the implementation because I am not familiar with Lavery splines, but I assume this is fine.

[JoshuaLampert]

Thanks!

[andrewwinters5000]

Thanks @vincmarks for initiating the implementation! This is a new feature that many will find useful.

[ranocha]

I have removed what allocations I can, but most of them are due to internal allocations within JuMP. So, the "proper" way to remove them would be to hand-roll the setup of the LP problem with sparse matrices. This could be interesting for a student project or something, but for the time being I am satisfied with the current strategy.

Nice 👍 Did you already open an issue for this?

[andrewwinters5000]

I have removed what allocations I can, but most of them are due to internal allocations within JuMP. So, the "proper" way to remove them would be to hand-roll the setup of the LP problem with sparse matrices. This could be interesting for a student project or something, but for the time being I am satisfied with the current strategy.

Nice 👍 Did you already open an issue for this?

I forgot, but can do so.