SciML
Performance differences with Interpolations.jl on ranges
Using the following code
using DataInterpolations, Interpolations, BenchmarkTools
N = 400
u = rand(N)
t = range(0.0, 1.0, length = N)
x = sort(rand(100)) # sorting required for Interpolations
println("cubic spline with DataInterpolations")
interp = DataInterpolations.CubicSpline(u,t)
@btime $interp.($x)
@btime $interp(0.5)
println("cubic spline with Interpolations")
interp2 = Interpolations.CubicSplineInterpolation(t, u)
@btime $interp2.($x)
@btime $interp2(0.5)
println("linear with DataInterpolations")
interplin = DataInterpolations.LinearInterpolation(u,t)
@btime $interplin.($x)
@btime $interplin(0.5)
println("linear with Interpolations")
interplin2 = Interpolations.LinearInterpolation(t, u)
@btime $interplin2.($x)
@btime $interplin2(0.5)
I get the output
cubic spline with DataInterpolations
10.900 μs (1 allocation: 896 bytes)
113.312 ns (0 allocations: 0 bytes)
cubic spline with Interpolations
2.167 μs (1 allocation: 896 bytes)
31.770 ns (0 allocations: 0 bytes)
linear with DataInterpolations
10.499 μs (1 allocation: 896 bytes)
108.972 ns (0 allocations: 0 bytes)
linear with Interpolations
3.290 μs (1 allocation: 896 bytes)
26.851 ns (0 allocations: 0 bytes)
Lots of room for mistakes in the crude benchmarking and eliminated code, but worth taking a look at. For versions
- DataInterpolations v1.0.0
- Interpolations v0.12.2
- Julia 1.2 on Windows
I suspect most of this might be because the linear and cubic splines are still using firstfirst instead of searchsortedfirst. In https://github.com/PumasAI/DataInterpolations.jl/pull/42 that made a significant difference for ZeroSpline.
I just ran the same code on with DataInterpolations v3.3.1 and Interpolations v0.13.1 and got these timings:
Julia v1.6.0-rc2:
cubic spline with DataInterpolations
38.427 μs (1 allocation: 896 bytes)
408.335 ns (0 allocations: 0 bytes)
cubic spline with Interpolations
1.156 μs (1 allocation: 896 bytes)
36.580 ns (0 allocations: 0 bytes)
linear with DataInterpolations
5.522 μs (1 allocation: 896 bytes)
53.387 ns (0 allocations: 0 bytes)
linear with Interpolations
423.192 ns (1 allocation: 896 bytes)
24.242 ns (0 allocations: 0 bytes)
Julia v1.5.4:
cubic spline with DataInterpolations
40.274 μs (1 allocation: 896 bytes)
417.327 ns (0 allocations: 0 bytes)
cubic spline with Interpolations
1.208 μs (1 allocation: 896 bytes)
37.548 ns (0 allocations: 0 bytes)
linear with DataInterpolations
4.624 μs (1 allocation: 896 bytes)
43.708 ns (0 allocations: 0 bytes)
linear with Interpolations
507.788 ns (1 allocation: 896 bytes)
26.072 ns (0 allocations: 0 bytes)
It seems the difference is a bit larger (did Interpolation.jl improve since last post?). If I'm not bringing anything useful here I apologize for the noise but I just thought I should share it since I ran it, just in case it's useful.
Current state:
julia> println("cubic spline with DataInterpolations")
cubic spline with DataInterpolations
julia> interp = DataInterpolations.CubicSpline(u,t);
julia> @btime $interp($x);
2.767 μs (1 allocation: 896 bytes)
julia> @btime $interp(0.5);
32.797 ns (0 allocations: 0 bytes)
julia> println("cubic spline with Interpolations")
cubic spline with Interpolations
julia> interp2 = Interpolations.CubicSplineInterpolation(t, u);
julia> @btime $interp2.($x);
550.802 ns (1 allocation: 896 bytes)
julia> @btime $interp2(0.5);
21.643 ns (0 allocations: 0 bytes)
julia>
julia> println("linear with DataInterpolations")
linear with DataInterpolations
julia> interplin = DataInterpolations.LinearInterpolation(u,t);
julia> @btime $interplin($x);
2.511 μs (1 allocation: 896 bytes)
julia> @btime $interplin(0.5);
31.325 ns (0 allocations: 0 bytes)
julia> println("linear with Interpolations")
linear with Interpolations
julia> interplin2 = Interpolations.LinearInterpolation(t, u);
julia> @btime $interplin2.($x);
255.398 ns (1 allocation: 896 bytes)
julia> @btime $interplin2(0.5);
16.100 ns (0 allocations: 0 bytes)
I changed the title because it seems to be about a missing specialization. If I do:
julia> N = 400;
julia> u = rand(N);
julia> t = Array(range(0.0, 1.0, length = N));
julia> x = sort(rand(100)); # sorting required for Interpolations
julia> interp2 = Interpolations.CubicSplineInterpolation(t, u);
ERROR: MethodError: no method matching cubic_spline_interpolation(::Vector{Float64}, ::Vector{Float64})
Closest candidates are:
cubic_spline_interpolation(::Tuple{Vararg{AbstractRange, N}}, ::AbstractArray{T, N}; bc, extrapolation_bc) where {N, T}
@ Interpolations C:\Users\accou\.julia\packages\Interpolations\91PhN\src\convenience-constructors.jl:28
cubic_spline_interpolation(::AbstractRange, ::AbstractVector; bc, extrapolation_bc)
@ Interpolations C:\Users\accou\.julia\packages\Interpolations\91PhN\src\convenience-constructors.jl:11
Stacktrace:
[1] CubicSplineInterpolation(::Vector{Float64}, ::Vararg{Vector{Float64}}; kwargs::@Kwargs{})
@ Interpolations .\deprecated.jl:105
[2] CubicSplineInterpolation(::Vector{Float64}, ::Vararg{Vector{Float64}})
@ Interpolations .\deprecated.jl:103
[3] top-level scope
@ REPL[53]:1
but:
julia> interp = DataInterpolations.CubicSpline(u,t)
CubicSpline with 400 points
┌────────────┬───────────┐
│ time │ u │
├────────────┼───────────┤
│ 0.0 │ 0.905018 │
│ 0.00250627 │ 0.241019 │
│ 0.00501253 │ 0.429529 │
│ 0.0075188 │ 0.649408 │
│ 0.0100251 │ 0.645401 │
│ 0.0125313 │ 0.385454 │
│ ⋮ │ ⋮ │
│ 0.987469 │ 0.513526 │
│ 0.989975 │ 0.410148 │
│ 0.992481 │ 0.0706234 │
│ 0.994987 │ 0.754565 │
│ 0.997494 │ 0.198778 │
│ 1.0 │ 0.243639 │
└────────────┴───────────┘
So what's really going on here is that Interpolations.jl only supports ranges, and it's specializing on the fact that ranges have equally spaced values so it improves performance in that case. Meanwhile, DataInterpolations.jl supports arbitrary spacing and does not specialize here.
@sathvikbhagavan maybe it would be good to look into specializing some of the common formulas like linear interpolation based on ranges.