Convert z in CubicSpline for performance

Project.toml

@@ -30,6 +30,7 @@ DataInterpolationsSymbolicsExt = "Symbolics"
 DataInterpolationsMakieExt = "Makie"
 
 [compat]
+AllocCheck = "0.2"
 Aqua = "0.8"
 BenchmarkTools = "1"
 ChainRulesCore = "1.25"
@@ -49,13 +50,15 @@ RegularizationTools = "0.6"
 SafeTestsets = "0.1"
 SparseConnectivityTracer = "1"
 StableRNGs = "1"
+StaticArrays = "1.9"
 Symbolics = "6.46"
 Test = "1.10"
 Unitful = "1.21.1"
 Zygote = "0.6.77, 0.7"
 julia = "1.10"
 
 [extras]
+AllocCheck = "9b6a8646-10ed-4001-bbdc-1d2f46dfbb1a"
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
@@ -68,10 +71,11 @@ RegularizationTools = "29dad682-9a27-4bc3-9c72-016788665182"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["Aqua", "BenchmarkTools", "SafeTestsets", "ChainRulesCore", "Optim", "RegularizationTools", "Test", "StableRNGs", "FiniteDifferences", "QuadGK", "ForwardDiff", "Symbolics", "Unitful", "Zygote", "SparseConnectivityTracer"]
+test = ["Aqua", "AllocCheck", "BenchmarkTools", "SafeTestsets", "ChainRulesCore", "Optim", "RegularizationTools", "Test", "StableRNGs", "FiniteDifferences", "QuadGK", "ForwardDiff", "StaticArrays", "Symbolics", "Unitful", "Zygote", "SparseConnectivityTracer"]

src/interpolation_caches.jl

@@ -717,7 +717,7 @@ function CubicSpline(
         1:(n + 1))
     d = transpose(reshape(reduce(hcat, d_), :, n + 1))
     z_ = reshape(transpose(tA \ d), size(u[1])..., :)
-    z = [z_s for z_s in eachslice(z_, dims = ndims(z_))]
+    z = [convert(eltype(u), z_s) for z_s in eachslice(z_, dims = ndims(z_))]
 
     p = CubicSplineParameterCache(u, h, z, cache_parameters)
     A = CubicSpline(
@@ -1463,7 +1463,7 @@ end
         interpolation_type::Type{<:AbstractInterpolation} = QuadraticSpline,
         kwargs...) where {U}
 
-Interpolate in a C¹ smooth way trough the data with unit speed by approximating
+Interpolate in a C¹ smooth way through the data with unit speed by approximating
 an interpolation (the shape interpolation) with line segments and circle segments.
 
 ## Arguments

src/interpolation_methods.jl

@@ -349,7 +349,7 @@ end
 
 # Cubic Hermite Spline
 function _interpolate(
-        A::CubicHermiteSpline{<:AbstractVector{<:Number}}, t::Number, iguess)
+        A::CubicHermiteSpline{<:AbstractVector}, t::Number, iguess)
     idx = get_idx(A, t, iguess)
     Δt₀ = t - A.t[idx]
     Δt₁ = t - A.t[idx + 1]
@@ -361,7 +361,7 @@ end
 
 # Quintic Hermite Spline
 function _interpolate(
-        A::QuinticHermiteSpline{<:AbstractVector{<:Number}}, t::Number, iguess)
+        A::QuinticHermiteSpline{<:AbstractVector}, t::Number, iguess)
     idx = get_idx(A, t, iguess)
     Δt₀ = t - A.t[idx]
     Δt₁ = t - A.t[idx + 1]

test/interpolation_tests.jl

@@ -6,6 +6,10 @@ using BenchmarkTools
 using Unitful
 using LinearAlgebra
 using Symbolics
+using AllocCheck: @check_allocs
+using StaticArrays: SVector, @SVector
+
+@check_allocs(test_allocs(itp, x) = itp(x)) # Reuse function definition to save on compilation time
 
 function test_interpolation_type(T)
     @test T <: DataInterpolations.AbstractInterpolation
@@ -80,6 +84,17 @@ end
         @test A(11) == [22.0, 33.0, 44.0]
         @test @inferred(output_dim(A)) == 1
         @test @inferred(output_size(A)) == (3,)
+
+        # Test allocation-free interpolation with StaticArrays
+        u_s = [convert(SVector{length(y), eltype(u_)}, i) for i in u]
+        @test @inferred(LinearInterpolation(
+            u_s, t; extrapolation = ExtrapolationType.Extension)) isa LinearInterpolation
+        A_s = LinearInterpolation(u_s, t; extrapolation = ExtrapolationType.Extension)
+        for _x in (0, 5.5, 11)
+            @test A(x) == A_s(x)
+        end
+        @test A_s(0) isa SVector{length(y)}
+        @test_nowarn test_allocs(A_s, 0)
     end
 
     x = 1:10
@@ -932,6 +947,21 @@ end
     test_cached_index(A)
     push!(u, 1.0)
     @test_throws AssertionError CubicHermiteSpline(du, u, t)
+
+    deleteat!(u, lastindex(u))
+    @testset "Vector of Vectors case" begin
+        u2 = [[u[i], u[i] + 1] for i in eachindex(u)]
+        du2 = [[du[i], du[i]] for i in eachindex(du)]
+        A2 = CubicHermiteSpline(du2, u2, t)
+        @test u2 ≈ A2.(t)
+    end
+    @testset "Vector of Matrices case" begin
+        u3 = [[u[i] u[i] + 1] for i in eachindex(u)]
+        du3 = [[du[i] du[i]] for i in eachindex(du)]
+        @test length(u3) == length(du3)
+        A3 = CubicHermiteSpline(du3, u3, t)
+        @test u3 ≈ A3.(t)
+    end
 end
 
 @testset "PCHIPInterpolation" begin
@@ -967,6 +997,22 @@ end
     @test_throws AssertionError QuinticHermiteSpline(ddu, du, u, t)
     @test @inferred(output_dim(A)) == 0
     @test @inferred(output_size(A)) == ()
+
+    deleteat!(u, lastindex(u))
+    @testset "Vector of Vectors case" begin
+        u2 = [[u[i], u[i] + 1] for i in eachindex(u)]
+        du2 = [[du[i], du[i]] for i in eachindex(du)]
+        ddu2 = [[ddu[i], ddu[i]] for i in eachindex(ddu)]
+        A2 = QuinticHermiteSpline(ddu2, du2, u2, t)
+        @test u2 ≈ A2.(t)
+    end
+    @testset "Vector of Matrices case" begin
+        u3 = [[u[i] u[i] + 1] for i in eachindex(u)]
+        du3 = [[du[i] du[i]] for i in eachindex(du)]
+        ddu3 = [[ddu[i] ddu[i]] for i in eachindex(ddu)]
+        A3 = QuinticHermiteSpline(ddu3, du3, u3, t)
+        @test u3 ≈ A3.(t)
+    end
 end
 
 @testset "Smooth Arc Length Interpolation" begin