Drop `AbstractArray` inheritance?

[tomasaschan]

Now that #206 is implemented and the main way to use Interpolations.jl is with function call semantics, I think we should evaluate what value AbstractArray semantics provide, and whether maybe we should drop it completely.

My only strong opinion here, is that we should try to minimize the complexity of the package (it's already quite complex...), so if there are features or functionality we can drop without disappointing our users, I think we should.

There are already a couple of places where dropping or changing the semantics of itp[x] have been mentioned, e.g. https://github.com/JuliaMath/Interpolations.jl/pull/238#issuecomment-425417528 and #227.

[timholy]

I have a ton of lab-private code which uses them like arrays, but I will take a peek sometime and evaluate how important that really is.

[tomasaschan]

This would obviously need a (quite long) deprecation period; completely dropping support would be a very long term goal.

[lstagner]

As I've mentioned in other threads I would like the getindex notation to index into the parent array e.g. itp[i] == parent(itp)[i]. Being able to store an object that can act as a function or as an array depending on the context is useful for my work. For instance, I calculate particle trajectories through a magnetic field. You can do this by solving an ODE in which case the function interface is the most useful. Alternatively, you can use a contouring algorithm (e.g. marching squares) to find the level set of a constant of motion which corresponds to a particle trajectory. In that case the array/getindex interface is the most convenient.

The contouring method is a quick and dirty way of calculating the trajectories, but there are cases when accuracy/precision is required and solving the ODE is preferred. I would like to support both methods using the same object.

[tomasaschan]

@lstagner Slightly off-topic, but what is your field of work/research? Your descriptions sound a lot like the fusion plasma physics work I did in my master's thesis... :)

[lstagner]

@tlycken Yes, I am a Postdoc at General Atomics doing fusion research, specifically in energetic particles hence the need for calculating particle trajectories.

[tomasaschan]

@lstagner Fun! The reason I got involved in the work that eventually led to creating this package was that I needed cubic b-splines for calculating alpha particle trajectories in tokamaks (specifically JET and ITER) :) You're probably past the point of learning much from it, but my thesis is available online if you'd be interested.

I'd love to see the results of your work at some point :)

[tomasaschan]

Re-reading this issue after a while, and I still think it seems like there are more cases when AbstractArray inheritance yields grievance than there are cases where it is absolutely necessary.

If your purpose, @lstagner, is to get at the underlying data, why not use the data from before you created the interpolation object? Reasonably, you have something like this today:

A = get_some_data()
itp = interpolate(A, ...)

Why can't you just do marching squares on A instead of on itp?

If that's for some reason impossible, it should be quite easy to implement a wrapping type that can do the translation from getindex to function calls, and define it only for integer indices to avoid all the ambiguity we have today.

I imagine something like this:

# This is shorthand for the actual AbstractInterpolation type
# I just didn't want to deal with all the other type parameters
# I also assume 1D here, but that's not a limiting factor
struct Itp{T}
  xs::Array{T,1}
end

# This type could be defined inside the package, so it gets access to all the non-exported
# introspection methods on interpolation objects
struct Wrapper{TEl, TItp} <: AbstractArray{TEl, 1}
  itp::TItp
end
Wrapper(itp::TItp) where TItp <: Itp{TEl} where TEl = Wrapper{TEl,TItp}(itp)

import Base: size, eltype, ndims, getindex

# These helper methods are merely for display purposes
size(itp::Itp) = size(itp.xs)
size(w::Wrapper) = size(w.itp)
eltype(::Type{TItp}) where TItp <: Itp{T} where T = T
eltype(::Type{TWrapper}) where TWrapper <: Wrapper{TEl,TItp} where TItp <: Itp{TEl} where TEl = TEl
ndims(itp::Itp) = 1
ndims(w::Wrapper) = 1

# This should be enough to get @lstagner's desired behavior
# AbstractArray overloads will handle all other cases.
getindex(w::Wrapper, xs...) = getindex(w.itp.xs, xs...)

# Of course, if specific interpolation types need specific handling
# (e.g. higher-order b-splines need to actually evaluate, rather than
# just pass on to the underlying array, because the array is not with
# the original data) they will have to implement overloads for the
# corresponding wrapper type.
# That could look something like
getindex(w::Wrapper{Interpolation{BSpline (etc)...}}, xs...) = w.itp(xs...)
# with suitable type constraints to ensure that we only accept single-
# point evaluation, and punt everything else to the AbstractArray
# machinery.

# Example usage:
itp = Itp(zeros(Float64, 3))
w = Wrapper(itp)

display(w)
# output:
# 3-element Wrapper{Float64,Itp{Float64}}:
#  0.0
#  0.0
#  0.0

With this, we can drop the AbstractArray inheritance of AbstractInterpolation entirely, without dropping support for getting at an AbstractArray with the underlying data.

[lstagner]

This doesn't actually solve my problem. It's not like I am going to code up the marching squares algorithm from scratch. I want to throw the interpolation object into the Contour.jl's contours function and just have it work. Similarly for plotting libraries. (plot(x,itp) should just work) Just supporting the getindex api is all well and good if you are not using any libraries, but the moment you want to use a function whose argument is restricted to be an AbstractArray then it wont work.

The advantage of subtyping AbstractArray is that you can use all the existing libraries that act on them without issue. In my opinion, the benefits of subtyping AbstractArray (but restricted to integer indices) far outweigh any down sides. In fact, if getindex calls into the parent array are restricted to only accept integers or other integer like indices I don't see any downside.

P.S. I put my thesis online if you would like to check it out. Basically, I came up with a way to infer the entire fast-ion distribution function from experimental measurements.

[tomasaschan]

@lstagner That still doesn't answer my question, though: if you have an itp of the data, don't you always have an array of the data as well? Why must itp be the object you index into/pass to Contours etc?

@timholy Your private lab-code probably uses getindex with non-integers, which means it'll break even if we keep inheritance on AbstractArray but restrict the use of getindex to integers (which is, I guess, what we're considering as the alternative here). Now that pkg supports workspaces, it shouldn't be difficult to pin an older version of Interpolations if you need to run one of those scripts without porting it.

[lstagner]

I'm designing a library that will (hopefully) be used by other people who won't know that they would have to call parent(itp) to do work with 95% of packages. It's just cleaner to pass the interpolation object.

Let me ask you this. What do we gain by dropping the AbstractArray subtyping? Slightly cleaner code? Is that worth losing the ability to automatically work with the majority of packages? I am fully in favor of dropping the getindex interface for interpolation, the function syntax is much more natural, but for accessing the parent array I think the getindex interface has a lot of benefits with practically no downside.

[timholy]

I'm in favor of dropping this interface. For my lab I suspect it won't be entirely trivial because I think I mix Interpolations together with some of the other AbstractArray wrapper types, but I don't really remember for sure, and in any event we should do what makes most sense. I just haven't put any time in this package for a while. I'll mark this issue up-for-grabs.

[stevengj]

I don't see what you gain by having this inheritance. Hardly any function that requires an AbstractArray argument will work with interpolation objects, because they don't support getindex and other standard features of AbstractArray. (Even when you had getindex in the past, it worked in such a non-arraylike way that it still would have failed with most array code.)

What you lose is that it confuses dispatch (in cases where functions do different things on array and non-array arguments), and makes error messages more confusing (since you get a MethodError deep inside a function call rather than being told directly that function foo requires an array).

[timholy]

At one point there was a drive to generalize AbstractArray indexing to well beyond integers. But I think that drive has gone by the wayside.