Add support for efficient high accuracy convolution interpolation (up to 7th order accurate C11 continuous) in any number of dimensions

[NikoBiele]

Convolution-Based Interpolation

This PR adds support for convolution-based interpolation methods to Interpolations.jl.

Features

• Support for continuous interpolation in any number of dimensions
• Fast implementation with precomputed kernels (I developed this method)
• Up to C11 continuity with up to 13th degree polynomials (I derived the kernels by generalizing the methodology of 10.1109/TASSP.1981.1163711 )
• High order of accuracy (generally 7th order accurate, except the cubic polynomial from the original paper which is 4th order)
• This method does not require solution of a linear system
• Compatible with existing extrapolation types (Line, Flat, etc.)
• Test coverage for dimensions up to 4D

Basic Usage Example

using Interpolations
using Plots
N_coarse = 5
x_coarse = range(0, 1, N_coarse)
y_coarse = range(0, 1, N_coarse)
data_coarse = rand(N_coarse, N_coarse)
itp = convolution_interpolation((x_coarse, y_coarse), data_coarse)
N_fine = 100
x_fine = range(0, 1, N_fine)
y_fine = range(0, 1, N_fine)
data_fine = [itp(x,y) for x in x_fine, y in y_fine]
p1 = contourf(x_coarse, y_coarse, data_coarse', 
             title="Original Data")
p2 = contourf(x_fine, y_fine, data_fine', 
             title="Smooth Interpolation")
plot(p1, p2, layout=(1,2), size=(800,300))

[codecov]

Codecov Report

Attention: Patch coverage is 0% with 488 lines in your changes missing coverage. Please review.

Project coverage is 69.38%. Comparing base (fab0496) to head (01fd289).
Report is 6 commits behind head on master.

Files with missing lines	Patch %	Lines
src/convolution/convolution_kernels.jl	0.00%	124 Missing ⚠️
src/convolution/convolution_coefs.jl	0.00%	98 Missing ⚠️
src/convolution/convolution_fast_interpolation.jl	0.00%	75 Missing ⚠️
src/convolution/convolution.jl	0.00%	70 Missing ⚠️
src/convolution/convolution_extrapolation.jl	0.00%	67 Missing ⚠️
...rc/convolution/convolution_kernel_interpolation.jl	0.00%	31 Missing ⚠️
src/convolution/linear_predict_coefs_derivation.jl	0.00%	15 Missing ⚠️
src/convenience-constructors.jl	0.00%	7 Missing ⚠️
src/Interpolations.jl	0.00%	1 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff             @@
##           master     #609       +/-   ##
===========================================
- Coverage   87.07%   69.38%   -17.69%     
===========================================
  Files          28       35        +7     
  Lines        1888     2375      +487     
===========================================
+ Hits         1644     1648        +4     
- Misses        244      727      +483     

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

[NikoBiele]

I submitted an early version of this previously, but ended up deleting the pull request (I wanted to create my own package).
Now I changed my mind again and want to contribute to this package instead, I hope it is ok.
Since I first submitted this pull request it has been 'supercharged' with higher accuracy and more efficient implementations. :)

[mkitti]

Some of this code looks generated or computed. Could you include the code or source for generating the values or code?

[mkitti]

Is some of this code machine generated? If so, could you elaborate how?

Also I wonder if it there might be a better method for evaluating some of the polynomials. Horners?

[NikoBiele]

Some of this code looks generated or computed. Could you include the code or source for generating the values or code?

Yes, i can submit the Python script for deriving the polynomial coefficients when I get home in a couple of days. But as I already stated I simply generalized the methodology of the paper I cited. I used SymPy for generalizing these derivations. In theory this script can derive arbitrary degree polynomials, however, the coefficients will get even longer than they already are. So I thought stopping at 13th degree was sufficient (even though higher degrees would likely perform even better).

[NikoBiele]

Is some of this code machine generated? If so, could you elaborate how?

Also I wonder if it there might be a better method for evaluating some of the polynomials. Horners?

Yes, I got assistance from Claude Sonnet 3.5, mostly for debugging code I wrote which didn't work initially.

You are likely correct that this isn't the most efficient way to evaluate the polynomials. I welcome improvements. 😊 But this was also the reason I developed the precompute method which is very efficient, since it is based on a single sorted lookup and a dot product. Due to the lookup it is highly efficient and almost independent of the polynomial degree. So maybe it is worth considering deriving even higher degree polynomials in the future, as I only saw improvements with higher degrees, thats why I kept going. 😊 For the boundary condition I use a linear prediction, where I have 'presolved' the linear systems for the coefficients. This symbolic linear solve I did using Symbolics.jl. 😊

[NikoBiele]

Is some of this code machine generated? If so, could you elaborate how?

Also I wonder if it there might be a better method for evaluating some of the polynomials. Horners?

All of these coefficients are analytically /symbolically derived. I would not trust generative AI with such a task of exact math, that would be bound to go wrong. For this task, SymPy and Symbolics.jl are the ideal choices. 😊

[NikoBiele]

I just want to demonstrate the improved frequency response when going from 3rd degree polynomial to my 13th degree polynomial kernel (this is part of my derivation script which I will upload soon). In this paper (10.1109/83.743854 ) they conclude that "... However, higher order schemes only yield marginal improvement, at an increased computational cost." What they forget is that they can include extra equations in the kernel to improve the frequency response and order of accuracy. In this paper the best cutoff slope achieved (at f=-0.5 Hz) is 2.538. My 13th degree cutoff slope at the same point is 3.8, and I can keep going higher without increasing the computational cost (if the time spent at the kernel precompute step can be neglected).

[NikoBiele]

Now I have comitted the Python script I mentioned, which I wrote to derive the kernels.
It is a generalization of the paper I cited.
It finds the kernels by setting up the system of equations and looping through all combinations of free coefficients, each time checking that the FFT satisifies the desired properties and it also finds the order of accuracy.
I should also mention that this script can be used to confirm the kernels which are already published in the literature, which confirms that this script is indeed a generalization.

[NikoBiele]

I remembered incorrectly when I said I used Symbolics.jl for deriving the linear prediction coefficients.
I actually used SymPy.jl, combined with the 'back slash' (\) native Julia operator, I found this approach to work the best.
Now I added that script too, so this pull request should be complete. :)

I should also mention the 'GaussianConvolutionKernel' which is in the 'convolution_kernels.jl'-file:
It can be used for smoothing noisy data, or data where the trust in each sample is limited, but one is more interested in the general smoothed picture. The degree of smoothing is controlled by the 'B' parameter which determines the 'width' of the number of samples used.

[mkitti]

Could you organize this as follows?

1. Move the generating code under a gen folder rather than the src folder.
2. For the Python script, create a folder under gen and include a pyproject.toml and preferably a uv or pixi lock file.
3. For the Julia SymPy code, create another folder under gen and include a Project.toml and Manifest.toml.
4. Additionally, could you include in the comments of the generated code the function calls or script to rin to generate that specific code?

The purpose of this is so that we can have a reproducible environment to generate the coefficients.

[mkitti]

Please reorganize the generating code under a top level gen folder as per my prior comment such that we have reproducible environments in which to execute the code.