JuliaMath / Interpolations.jl

Fast, continuous interpolation of discrete datasets in Julia
http://juliamath.github.io/Interpolations.jl/
Other
533 stars 109 forks source link

Add support for WENO4 interpolation #470

Open tiagopereira opened 2 years ago

tiagopereira commented 2 years ago

This PR is not currently working. I would like to add support for interpolation with a 4th order Weighted Essentially Non-Oscillatory (WENO) scheme, following Janett et al (2019). I wrote the standalone package WENO4.jl (not yet registered) that has been tested to work as it should. This PR is an attempt at bringing the same functionality into Interpolations.jl.

My experience with Julia is limited, and I don't fully understand the abstractions and architecture of Interpolations.jl, so I apologise for the many mistakes in the current PR and would appreciate some help at getting this code working and sufficiently general. This implementation follows closely what I saw in monotonic.jl, but I couldn't get my head around all the abstractions. This is the type of usage I was hoping it would get when working:

using Interpolations

itp = interpolate(x, y, interpolate(x, y, WENO4Interpolation())
new_values = itp.(new_x)

The idea is to follow existing practice of getting a way to generate an interpolant for all input data values and nodes, and then a separate function to evaluate this interpolant at a required point. I have tested and the functions inside, in particular interpolate (to build the interpolant) and the function to evaluate the interpolant, currently

function (itp::WENO4Interpolation)(x::Number)

works if I rewrite it as the standalone function

function evaluate(itp::WENO4Interpolation, x::Number)

However, it is not working in the present formulation because I don't really understand the relation between the interpolate function and the function with the same name as the struct, here called WENO4Interpolation.

Currently the functions used perform about twice as slow as in the standalone package, which I attribute to the need of generating intermediate arrays. This could also be improved.

codecov[bot] commented 2 years ago

Codecov Report

Merging #470 (7c7e11c) into master (3421847) will decrease coverage by 5.32%. The diff coverage is 0.00%.

@@            Coverage Diff             @@
##           master     #470      +/-   ##
==========================================
- Coverage   86.75%   81.42%   -5.33%     
==========================================
  Files          27       28       +1     
  Lines        1819     1938     +119     
==========================================
  Hits         1578     1578              
- Misses        241      360     +119     
Impacted Files Coverage Δ
src/Interpolations.jl 79.52% <ø> (ø)
src/weno/weno.jl 0.00% <0.00%> (ø)

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update 3421847...7c7e11c. Read the comment docs.

mkitti commented 2 years ago

interpolate is a general method that will eventually construct an instance of the interpolation object. (itp::WENO4Interpolation)(x...) is used when you try to use an instance of the struct as a function.

julia> struct Greeting 
           salutation::String
       end # define a struct

julia> hello = Greeting("Hello") # construct a Greeting instance with salutation == "Hello"
Greeting("Hello")

julia> (g::Greeting)(name) = println(g.salutation * ", " * name) # using an instance

julia> hello("Tiago") # Use the instance as a function (e.g. a Functor)
Hello, Tiago

julia> hola = Greeting("Hola") # construct another Greeting instance with salutation == "Hola"
Greeting("Hola")

julia> hola("Marcos")
Hola, Marcos

julia> Greeting() = Greeting("Hi") # define a constructor for Greeting that takes no arguments... and defaults salutation to "Hi"
Greeting

julia> default = Greeting() # use the new constructor
Greeting("Hi")

julia> default("Mark")
Hi, Mark