Interpolating along a single dimension of a multi dimensional array for particular points

General Usage
1

I have some simple linear interpolation code in a package. It would be nice to support spline etc and not write that myself - but it has to be extremely fast and non-allocating at the point of use. It can take as long as it likes to define an interpolator object as that happens outside the main loop.

The problem is to interpolate along a single dimension (time) of a 3 dimensional long/lat/time array. This happens one lat/long position at a time only a small (but changing) fraction of the lat/long grid is actually calculated, for performance reasons.

Can Interpolations.jl or another package do this without creating a view along the time axis?

2

You can broadcast over GitHub - PumasAI/DataInterpolations.jl: A library of data interpolation and smoothing functions . More can be done but that should be a good start. I don’t think we implemented multiple 1D interpolations at one with this yet, but it would be a good idea.

1 Like
3

Thanks I wasn’t aware of that package, it looks great.

Unfortunately we can’t use broadcast, the key optimisation in DynamicGrids.jl is to avoid unnecessary work by not broadcasting the whole array. It would be good to do something like this where the specified dimension(s) are interpolated:

u = rand(5, 5, 5)
t = 0:4
interp = LinearInterpolation(u, t; dims=3)
interp(1, 2, 3.5)

As a side point DataInterpolation.jl could play really well with DimensionalData.jl if it was generalised like this, you could specify the interpolation dimension instead of t and dims.

4

that looks like a fairly reasonable suggestion

5

It would also be good to get an intermediate object that can apply an interpolation, like:

u = rand(5, 5, 5)
t = 0:4
interp = LinearInterpolation(u, t; dims=3)
partialinterp = interp(3.5) 
partialinterp(1, 2)

Then the calculation of indices could be done before the actual interpolation, outside of main loops.

But I’m not sure how widely this kind of selective/non-broadcasted version would be required.

6

Interpolations should be just fine doing this kind of thing. It’s not well documented, but Interpolations (now) does interpolation in the naive way, which IMO is a strength. For example, this is how the internals for linear interpolation in 1d actually work:

julia> using Interpolations

# Build a "weighted adjacent index" for indexing into any array: compute
#    0.2*a[j] + 0.8*a[j+1], with j=3

julia> widx = Interpolations.WeightedAdjIndex(3, (0.2, 0.8))
Interpolations.WeightedAdjIndex{2,Float64}(3, (0.2, 0.8))

julia> r = 1:10
1:10

julia> r[widx]    # you can index arbitrary arrays with these weird indexes
3.8000000000000003

# Now do something kinda silly, but it shows you're not constrained in what you can do

julia> widx = Interpolations.WeightedAdjIndex(3, (0.0, 0.0, 0.0, 0.7))
Interpolations.WeightedAdjIndex{4,Float64}(3, (0.0, 0.0, 0.0, 0.7))

julia> r[widx]
4.199999999999999

julia> 0.7*r[6]      # the 4 weights apply to positions 3, 4, 5, and 6, hence the r[6]
4.199999999999999

# Try a 2d example

julia> A = reshape(1:20, 4, 5)
4×5 reshape(::UnitRange{Int64}, 4, 5) with eltype Int64:
 1  5   9  13  17
 2  6  10  14  18
 3  7  11  15  19
 4  8  12  16  20

# Try indexing the first dimension with an integer and the 2nd with a weighted index
julia> i, j = 2, Interpolations.WeightedAdjIndex(2, (0.3, 0.7))
(2, Interpolations.WeightedAdjIndex{2,Float64}(2, (0.3, 0.7)))

julia> A[i,j]
8.8

julia> 0.3*A[i,2] + 0.7*A[i,3]   # see, same as linear interpolation along dimension 2
8.8

# Interpolate in both dimensions, this time trying out the other index type, WeightedArbIndex.
# This type allows the index locations to be arbitrary, rather than incremental. It is used internally to
# implement things like periodic boundary conditions, reflecting boundary conditions, etc.
julia> i = Interpolations.WeightedArbIndex((1, 3), (-0.1, 100.0))
Interpolations.WeightedArbIndex{2,Float64}((1, 3), (-0.1, 100.0))

julia> A[i, j]
979.2199999999999

# Indeed, this is the full 2d interpolation
julia> A[1,2]*(-0.1*0.3) + A[1,3]*(-0.1*0.7) + A[3,2]*(100*0.3) + A[3,3]*(100*0.7)
979.22

As you can see, Interpolations puts all the magic into the index objects, not the arrays. This makes it very flexible. Moreover, if you put stuff inside @inbounds then the generated code should be very efficient.

Interpolations does include a NoInterp mode as an exported interface for doing these kinds of things, but the real magic is in the indexes.

2 Likes
7

That’s perfect, thanks for the detailed example. I agree that kind naive composability is a real strength. Most of the interpolation I need day to day (in dynamic ecological models) are only for some subset of an array, optimizations on broadcast just cant compete with not doing 90% of the work.

It might be worth putting this in the docs.

Edit: I’m wondering if splines work in a similar way? It would be good to be able to structure a package so that users can easily switch between rounding, linear interpolation and splines - depending on their accuracy/performance tradeoff.