Interpolations.jl

Interpolations.jl copied to clipboard

Hi! I would like to add some new convenience copnstructors and would create a new PR if that was supposed to be a good idea and something that should be part of this package.

What I mean here is the following: The command linear_interpolation(x, y) is for example a shortcut for extrapolate(interpolate((x, ), y, Gridded(Linear())), Throw()). In the same way I would, for example for the recently created AkimaMotonicInterpolation type, make it possible to call akima_spline(x, y) as a shortcut for extrapolate(interpolate(x, y, AkimaMontonicInterpolation()), Throw()).

What are your thoughts on this idea?

One issue is that this definitely could explode. We could create a lot more convenience constructors here.

What I'm wondering about is if we should just have simplified front end overall.

interpolation(x, y, type= :linear)
interpolation(x, y, type= :akima_monotonic)
interpolation(x, y, type= :akima_monotonic, extrapolation = :throw)

While this is not going to be very type stable, maybe something like this would be preferable to creating a bunch of different methods.

What do you think?

One issue is that this definitely could explode. We could create a lot more convenience constructors here.

I see your point here. But then again I wonder if this would do the package any kind of harm in terms of performance or if this would just mean a bunch of rather boring work.

---

What I'm wondering about is if we should just have simplified front end overall.

interpolation(x, y, type= :linear)
interpolation(x, y, type= :akima_monotonic)
interpolation(x, y, type= :akima_monotonic, extrapolation = :throw)

While this is not going to be very type stable, maybe something like this would be preferable to creating a bunch of different methods.

That looks like a really nice approach, too. Obviously, that would also make it more convenient for programmers and we would get rid of that rather ugly current notation. I guess that in the end it might just be a matter of subjecitve taste. I would personally tend to vote for creating new methods, but that's probably just because I already got quite used to linear_interpolation and cubic_spline_interpolation 😂 Moreover I could not really help implementing your alternative by crating a PR unless someone showed me at least one example because I would not have any clue how to do it "correctly".

Part of the issue is that using linear_interpolation and cubic_spline_interpolation is often very bad for performance because they include extrapolation. https://github.com/JuliaMath/Interpolations.jl/issues/462

The current interfaces emphasize type stability. Based on the input types and the method, we should be able to predict the output type.

What I proposed above is not necessarily type stable in isolation. However, constant propagation may be good enough at this point to curtail those issues. Consider the following example:

julia> struct A
           x::Int
           y::Float64
       end

julia> a = A(1,2.)
A(1, 2.0)

julia> @code_warntype getproperty(a, :x)
MethodInstance for getproperty(::A, ::Symbol)
  from getproperty(x, f::Symbol) in Base at Base.jl:42
Arguments
  #self#::Core.Const(getproperty)
  x::A
  f::Symbol
Body::Union{Float64, Int64}
1 ─      nothing
│   %2 = Base.getfield(x, f)::Union{Float64, Int64}
└──      return %2

julia> f(s) = getproperty(s, :x)
f (generic function with 1 method)

julia> @code_warntype f(a)
MethodInstance for f(::A)
  from f(s) in Main at REPL[6]:1
Arguments
  #self#::Core.Const(f)
  s::A
Body::Int64
1 ─ %1 = Main.getproperty(s, :x)::Int64
└──      return %1

In the first case, we see there is ambiguity, Union{Float64, Int64}, on what will be returned since the evaluation is for getproperty(::A, ::Symbol). In the second case, we see that constant propagation is correctly able to to determine that the return type should be Int64.

These technical details are very interesting. If performance and type stabilty matter that much than I guess I must admit that this issue of unconvenient notation is much harder to solve than I expected. I might have a deeper look into that if when I have more time. For now I guess I will discard the original plan of introducting akima_spline(x, y) and so on since that does not look like the best solution here.