Comparison with Interpolations.jl

[henry2004y]

Hi,

I compared the usage and performance of DataInterpolationsND.jl v0.1.2 and Interpolations.jl v0.16.2, to see the benefits of this new package.

Comparison script

using Test
using StaticArrays
using DataInterpolationsND
import DataInterpolationsND as dnd
using Interpolations
using Chairmarks

# Function to setup fields and return raw data for comparison
function setup_cartesian_data()
    x = range(-10.0, 10.0, length = 4)
    y = range(-10.0, 10.0, length = 4)
    z = range(-10.0, 10.0, length = 4)

    B = zeros(3, length(x), length(y), length(z)) # vector
    B[3, :, :, :] .= 1.0e-8 # 10 nT

    A = zeros(length(x), length(y), length(z)) # scalar
    A[:, :, :] .= 1.0e-8

    return B, A, x, y, z
end

struct FieldInterpolator{T} <: Function
    itp::T
end

function (fi::FieldInterpolator)(xu)
    return fi.itp(xu[1], xu[2], xu[3])
end


# 1. Get raw data
B, A, x, y, z = setup_cartesian_data()

# --- Interpolations.jl ---
# Scalar Field A
# Interpolations.jl expects (x,y,z) order if data is (x,y,z)
itp_A_raw = extrapolate(
    scale(
        interpolate(A, BSpline(Linear())),
        x, y, z
    ),
    NaN
)

itp_A = FieldInterpolator(itp_A_raw)

# Vector Field B with Interpolations.jl
# Usually easier to interpolate array of SVectors, or simpler 3xNxMxP
B_sv = reinterpret(reshape, SVector{3, Float64}, B)
itp_B_raw = extrapolate(
    scale(
        interpolate(B_sv, BSpline(Linear())),
        x, y, z
    ),
    SVector{3, Float64}(NaN, NaN, NaN)
)

itp_B = FieldInterpolator(itp_B_raw)

# --- DataInterpolationsND.jl ---
# Scalar Field A
# DataInterpolationsND expects spatial dims first in data array: (Nx, Ny, Nz)
args_grid = (
    dnd.LinearInterpolationDimension(collect(x)),
    dnd.LinearInterpolationDimension(collect(y)),
    dnd.LinearInterpolationDimension(collect(z)),
)
itp_A_dnd_raw = dnd.NDInterpolation(A, args_grid)

itp_A_dnd = FieldInterpolator(itp_A_dnd_raw)

# Vector Field B
# DataInterpolationsND expects spatial dims first, then extra components at end.
# Current B is (3, Nx, Ny, Nz). We need (Nx, Ny, Nz, 3).
B_perm = permutedims(B, (2, 3, 4, 1))
itp_B_dnd_raw = dnd.NDInterpolation(B_perm, args_grid)
# Note: This is not working!
# itp_B_dnd_raw = dnd.NDInterpolation(B_sv, args_grid)

itp_B_dnd = FieldInterpolator(itp_B_dnd_raw)

# --- Comparison ---
loc = SA[1.0, 1.0, 1.0]

# Evaluate Interpolations.jl
val_A_itp_raw = itp_A_raw(loc...)
val_B_itp_raw = itp_B_raw(loc...)

val_A_itp = itp_A(loc)
val_B_itp = itp_B(loc)

# Evaluate DataInterpolationsND.jl
val_A_dnd_raw = itp_A_dnd_raw(loc...)
val_A_dnd = itp_A_dnd(loc)
# DataInterpolationsND mutates output buffer for >1 dim
out_B = zeros(3)
itp_B_dnd_raw(out_B, loc...)
# val_B_dnd = SVector{3}(out_B)

val_B_dnd = itp_B_dnd(loc)

println("Location: $loc")
println("Vector Field (B):")
println("  Interpolations.jl:      $val_B_itp")
println("  DataInterpolationsND.jl: $val_B_dnd")
println("Scalar Field (A):")
println("  Interpolations.jl:      $val_A_itp")
println("  DataInterpolationsND.jl: $val_A_dnd")

# Tests
@test val_B_itp ≈ val_B_dnd
@test val_A_itp ≈ val_A_dnd

# Verify expected values (from setup)
@test val_B_dnd[3] ≈ 1.0e-8
@test val_A_dnd ≈ 1.0e-8

# Benchmark
println("\nBenchmarking...")

# Benchmark Interpolations.jl
println("Interpolations.jl:")
println("Scalar Field (A):")
# display(@be $itp_A_raw($loc[1], $loc[2], $loc[3]))
display(@be $itp_A($loc))
println("Vector Field (B):")
# display(@be $itp_B_raw($loc[1], $loc[2], $loc[3]))
display(@be $itp_B($loc))

# Benchmark DataInterpolationsND.jl
println("DataInterpolationsND.jl:")
println("Scalar Field (A):")
# display(@be $itp_A_dnd_raw($loc[1], $loc[2], $loc[3]))
display(@be $itp_A_dnd($loc))
println("Vector Field (B):")
# display(@be $itp_B_dnd_raw($loc[1], $loc[2], $loc[3]))
display(@be $itp_B_dnd($loc))
# display(@be $itp_B_dnd_raw($out_B, $loc[1], $loc[2], $loc[3]))

The test code linear interpolation outputs are identical:

Location: [1.0, 1.0, 1.0]
Vector Field (B):
  Interpolations.jl:      [0.0, 0.0, 1.0e-8]
  DataInterpolationsND.jl: [0.0, 0.0, 1.0e-8]
Scalar Field (A):
  Interpolations.jl:      1.0e-8
  DataInterpolationsND.jl: 1.0e-8

For benchmarking with Chairmarks.jl for the scalar and vector fields with a single point evaluation, I found that my vector field interpolation version of DataInterpolationsND.jl allocates, which I don't know how to avoid:

Benchmarking...
Interpolations.jl:
Scalar Field (A):
Benchmark: 3382 samples with 933 evaluations
 min    27.760 ns
 median 29.153 ns
 mean   30.143 ns
 max    90.890 ns
Vector Field (B):
Benchmark: 3058 samples with 778 evaluations
 min    35.733 ns
 median 37.532 ns
 mean   39.085 ns
 max    148.972 ns
DataInterpolationsND.jl:
Scalar Field (A):
Benchmark: 3333 samples with 923 evaluations
 min    28.277 ns
 median 29.794 ns
 mean   31.202 ns
 max    456.771 ns
Vector Field (B):
Benchmark: 2313 samples with 221 evaluations
 min    107.240 ns (18 allocs: 720 bytes)
 median 119.457 ns (18 allocs: 720 bytes)
 mean   186.193 ns (18 allocs: 720 bytes, 0.43% gc time)
 max    28.042 μs (18 allocs: 720 bytes, 99.03% gc time)

So my questions are:

1. How to optimize the vector field interpolation in the example above? I also notice that DataInterpolationsND expects spatial dims first, then extra components at end.
2. In what kind of benchmarks would I see significant performance benefit over Interpolations.jl?

Thanks!

[SouthEndMusic]

Hi @henry2004y,

Thanks for taking the time to make this comparison. To answer your questions:

1. The reason why DataInterpolationsND.jl allocates in this case is due to the logic where the output is generated; which at the moment for non-scalar outputs is an array similar to u of the required size. The package hasn't been optimized (yet) to generate StaticArray outputs, because the emphasis has been on vectorized interpolation evaluation. Regarding your note: the following PR generalizes things, so that the interpolating dimensions no longer have to be leading: Generalise _interpolate! to allow mixed dimensions #33.
2. I don't know, I'm not very familiar with Interpolations.jl. To get a sense of the use cases of this package see https://discourse.julialang.org/t/ann-datainterpolationsnd-jl/130883 and Interpolation of 2d data DataInterpolations.jl#402.