Interpolate along single dimension of multidimensional array

General Usage
1

I’m trying to convert a Python code to Julia that makes use of scipy.interpolate.interp1d, and I’m having difficulty finding the best way to do this in Julia. The interp1d function does 1D interpolation along a given dimension of multidimensional array.

I’ve been using Interpolations.jl but my approach does not seem the most elegant (see below). There is a very similar question asked three years ago, but the accepted answer does not provide a straightforward answer.

Here’s a simple example of what I’m looking for:

# Test data
x = [0.,  2., 5., 10., 11.]
y = repeat(x .^ 2, outer=(1, 3, 3))

# We want to interpolate y (a 3D array) along the first dimension,
# so that the output is another 3D array
itp = interp1d(x, y; dims=1, k=1)

# Expected results
itp(2.0)
3×3 Matrix{Float64}:
 4.0  4.0  4.0
 4.0  4.0  4.0
 4.0  4.0  4.0

itp(3)
3×3 Matrix{Float64}:
 11.0  11.0  11.0
 11.0  11.0  11.0
 11.0  11.0  11.0

My data are all in irregular grids.

This seems to work, but I wonder if it is the best approach:

itp = interpolate((x, 1:3, 1:3), y, Gridded(Linear()))

itp(3, 1:3, 1:3)
3×3 Matrix{Float64}:
 11.0  11.0  11.0
 11.0  11.0  11.0
 11.0  11.0  11.0

And it still seems that with Gridded the highest order interpolation supported is linear.
Anyone else has a better approach for this type of interpolation?

2

Just in case, there is an Interp1d.jl package from one of the Python contributors, with an accompanying notebook with examples.

3

Thanks for the suggestion. I had looked at Interp1d.jl before, but it seems very alpha. Not registered package, no documentation, does not work in 3D or higher (only 1D or 2D) and does not even support linear interpolation (only nearest neighbour or previous).

4

For reference, the solution above with Interpolations.jl seems to be about 50% slower than the Python implementation:

import numpy as np
from scipy.interpolate import interp1d

tmp = np.random.rand(100, 100, 100).astype('f')
x = np.sort(np.random.rand(100).astype('f'))
new_x = np.linspace(x.min(), x.max(), 100).astype('f')
itp = interp1d(x, tmp, axis=-1, kind='linear')

%timeit itp(new_x)
3.08 ms ± 10.8 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Julia:

using Interpolations
using BenchmarkTools

tmp = rand(Float32, (100, 100, 100))
x = sort(rand(Float32, 100))
new_x = LinRange(minimum(x), maximum(x), 100)
itp = interpolate((x, 1:100, 1:100), tmp, Gridded(Linear()))

@benchmark itp(new_x, 1:100, 1:100)
BenchmarkTools.Trial: 1059 samples with 1 evaluation.
 Range (min … max):  4.434 ms …   6.687 ms  ┊ GC (min … max): 0.00% … 31.21%
 Time  (median):     4.619 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   4.724 ms ± 343.316 μs  ┊ GC (mean ± σ):  2.09% ±  5.69%

      ▃██▅▂▂▁
  █▅▄▁████████▅▅▁▄▁▄▅▄▁▁▁▁▁▁▁▁▄▁▁▁▁▁▅▅▅▅▇▇▇▇▅▅▆▇▆▅▅▅▅▄▅▅▄▄▄▅▄ █
  4.43 ms      Histogram: log(frequency) by time      6.19 ms <

 Memory estimate: 3.83 MiB, allocs estimate: 7.
5

It’s common in Julia to make your function work on 1-d arrays (vectors) and then just map it over slices of n-d arrays.
Suppose you already have interpolation working exactly as you need for vectors, implemented in func(vector). Then mapslices(func, A; dims=1) should work, or map(func, eachslice(...)) - depending on how you prefer the result to be shaped.

6

Telling interpolations to not perform interpolation in the second two axes produces better performance on my machine:

julia> @benchmark itp(new_x, 1:100, 1:100)
BenchmarkTools.Trial: 376 samples with 1 evaluation.
 Range (min … max):  12.599 ms … 19.921 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     12.970 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   13.298 ms ±  1.044 ms  ┊ GC (mean ± σ):  0.31% ± 1.15%

  ▄▇█▇▆▄▃▁▁▃ ▁
  ██████████▇██▆▇▄▆▇▆▁▄▁▁▁▁▆▇▄▁▆▁▄▁▁▁▁▁▄▁▁▁▁▁▄▄▁▆▁▁▁▁▁▁▁▁▄▄▁▄ ▆
  12.6 ms      Histogram: log(frequency) by time      18.6 ms <

 Memory estimate: 3.83 MiB, allocs estimate: 7.

julia> itp2 = interpolate((x, 1:100, 1:100), tmp, (Gridded(Linear()), NoInterp(), NoInterp()));

julia> @benchmark itp2(new_x, 1:100, 1:100)
BenchmarkTools.Trial: 2001 samples with 1 evaluation.
 Range (min … max):  2.210 ms …   6.351 ms  ┊ GC (min … max): 0.00% … 18.00%
 Time  (median):     2.311 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.485 ms ± 484.836 μs  ┊ GC (mean ± σ):  1.90% ±  5.77%

  ██▇▆▅▃▃▂▁       ▂▂ ▁▁                                       ▁
  █████████████▇▇▇██████▆▇▆▇▄▅▆▁▆▅▅▅▅▆▅▆▅▅▅▄▄▄▅▅▅▅▆▅▄▅▅▅▁▄▁▁▅ █
  2.21 ms      Histogram: log(frequency) by time      4.88 ms <

 Memory estimate: 3.82 MiB, allocs estimate: 7.
2 Likes
7

Thank you! That is a great suggestion. Forgot that by default it was also interpolating in the two other axes.