Add missing matrix support for AkimaInterpolation and QuadraticSpline

src/integrals.jl

@@ -252,6 +252,33 @@ function _integral(
     Δt * (α * (t2_rel^2 + t1_rel * t2_rel + t1_rel^2) / 3 + β * (t2_rel + t1_rel) / 2 + uᵢ)
 end
 
+function _integral(
+        A::QuadraticSpline{<:AbstractMatrix}, idx::Number, t1::Number, t2::Number)
+    tᵢ = A.t[idx]
+    t1_rel = t1 - tᵢ
+    t2_rel = t2 - tᵢ
+    Δt = t2 - t1
+
+    # Compute integral for each row
+    result = similar(A.u, size(A.u, 1))
+    for row in 1:size(A.u, 1)
+        # Compute parameters on the fly for this row (same logic as in vector _integral)
+        tᵢ₊ = (A.t[idx] + A.t[idx + 1]) / 2
+        sc = zeros(eltype(A.t), length(A.t))
+        nonzero_coefficient_idxs = spline_coefficients!(sc, 2, A.k, tᵢ₊)
+        uᵢ₊ = zero(eltype(A.u))
+        for j in nonzero_coefficient_idxs
+            uᵢ₊ += sc[j] * A.c[row, j]
+        end
+        α_row = 2 * (A.u[row, idx + 1] + A.u[row, idx]) - 4uᵢ₊
+        β_row = 4 * (uᵢ₊ - A.u[row, idx]) - (A.u[row, idx + 1] - A.u[row, idx])
+        uᵢ_row = A.u[row, idx]
+        result[row] = Δt * (α_row * (t2_rel^2 + t1_rel * t2_rel + t1_rel^2) / 3 + β_row * (t2_rel + t1_rel) / 2 + uᵢ_row)
+    end
+
+    result
+end
+
 function _integral(
         A::CubicSpline{<:AbstractVector{<:Number}}, idx::Number, t1::Number, t2::Number)
     tᵢ = A.t[idx]
@@ -266,6 +293,16 @@ function _integral(A::AkimaInterpolation{<:AbstractVector{<:Number}},
     integrate_cubic_polynomial(t1, t2, A.t[idx], A.u[idx], A.b[idx], A.c[idx], A.d[idx])
 end
 
+function _integral(A::AkimaInterpolation{<:AbstractMatrix},
+        idx::Number, t1::Number, t2::Number)
+    # Compute integral for each row
+    result = similar(A.u, size(A.u, 1))
+    for row in 1:size(A.u, 1)
+        result[row] = integrate_cubic_polynomial(t1, t2, A.t[idx], A.u[row, idx], A.b[row, idx], A.c[row, idx], A.d[row, idx])
+    end
+    result
+end
+
 function _integral(A::LagrangeInterpolation, idx::Number, t1::Number, t2::Number)
     throw(IntegralNotFoundError())
 end

src/interpolation_caches.jl

@@ -281,7 +281,13 @@ function AkimaInterpolation(
     m[end] = 2m[end - 1] - m[end - 2]
 
     b = (m[4:end] .+ m[1:(end - 3)]) ./ 2
-    dm = abs.(diff(m))
+    dm = if eltype(u) <: AbstractVector
+        # For Vector of Vectors, use norm instead of abs
+        norm.(diff(m))
+    else
+        # For scalar values, use abs as before
+        abs.(diff(m))
+    end
     f1 = dm[3:(n + 2)]
     f2 = dm[1:n]
     f12 = f1 + f2
@@ -299,6 +305,57 @@ function AkimaInterpolation(
         extrapolation_right, cache_parameters, linear_lookup)
 end
 
+function AkimaInterpolation(
+        u::AbstractMatrix, t; extrapolation::ExtrapolationType.T = ExtrapolationType.None,
+        extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
+        extrapolation_right::ExtrapolationType.T = ExtrapolationType.None, cache_parameters = false, assume_linear_t = 1e-2)
+    extrapolation_left,
+    extrapolation_right = munge_extrapolation(
+        extrapolation, extrapolation_left, extrapolation_right)
+    u, t = munge_data(u, t)
+    linear_lookup = seems_linear(assume_linear_t, t)
+    n = length(t)
+    dt = diff(t)
+
+    # Handle each row of the matrix independently
+    nrows = size(u, 1)
+    b = similar(u, nrows, n)
+    c = similar(u, nrows, n - 1) 
+    d = similar(u, nrows, n - 1)
+
+    for row in 1:nrows
+        u_row = u[row, :]
+        m = Array{eltype(u_row)}(undef, n + 3)
+        m[3:(end - 2)] = diff(u_row) ./ dt
+        m[2] = 2m[3] - m[4]
+        m[1] = 2m[2] - m[3]
+        m[end - 1] = 2m[end - 2] - m[end - 3]
+        m[end] = 2m[end - 1] - m[end - 2]
+
+        b_row = (m[4:end] .+ m[1:(end - 3)]) ./ 2
+        dm = abs.(diff(m))
+        f1 = dm[3:(n + 2)]
+        f2 = dm[1:n]
+        f12 = f1 + f2
+        ind = findall(f12 .> 1e-9 * maximum(f12))
+        b_row[ind] = (f1[ind] .* m[ind .+ 1] .+
+                      f2[ind] .* m[ind .+ 2]) ./ f12[ind]
+        c_row = (3 .* m[3:(end - 2)] .- 2 .* b_row[1:(end - 1)] .- b_row[2:end]) ./ dt
+        d_row = (b_row[1:(end - 1)] .+ b_row[2:end] .- 2 .* m[3:(end - 2)]) ./ dt .^ 2
+
+        b[row, :] = b_row
+        c[row, :] = c_row  
+        d[row, :] = d_row
+    end
+
+    A = AkimaInterpolation(
+        u, t, nothing, b, c, d, extrapolation_left,
+        extrapolation_right, cache_parameters, linear_lookup)
+    I = cumulative_integral(A, cache_parameters)
+    AkimaInterpolation(u, t, I, b, c, d, extrapolation_left,
+        extrapolation_right, cache_parameters, linear_lookup)
+end
+
 """
     ConstantInterpolation(u, t; dir = :left, extrapolation::ExtrapolationType.T = ExtrapolationType.None, extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
         extrapolation_right::ExtrapolationType.T = ExtrapolationType.None, cache_parameters = false)
@@ -562,6 +619,36 @@ function QuadraticSpline(
         extrapolation_right, cache_parameters, assume_linear_t)
 end
 
+function QuadraticSpline(
+        u::AbstractMatrix, t; extrapolation::ExtrapolationType.T = ExtrapolationType.None,
+        extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
+        extrapolation_right::ExtrapolationType.T = ExtrapolationType.None, cache_parameters = false,
+        assume_linear_t = 1e-2)
+    extrapolation_left,
+    extrapolation_right = munge_extrapolation(
+        extrapolation, extrapolation_left, extrapolation_right)
+    u, t = munge_data(u, t)
+
+    n = length(t)
+    dtype_sc = typeof(t[1] / t[1])
+    sc = zeros(dtype_sc, n)
+    k, A = quadratic_spline_params(t, sc)
+
+    # For matrices, solve the system for each row/column independently
+    c = similar(u)  # Same shape as u
+    for i in axes(u, 1)
+        c[i, :] = A \ u[i, :]
+    end
+
+    p = QuadraticSplineParameterCache(u, t, k, c, sc, cache_parameters)
+    A_spline = QuadraticSpline(
+        u, t, nothing, p, k, c, sc, extrapolation_left,
+        extrapolation_right, cache_parameters, assume_linear_t)
+    I = cumulative_integral(A_spline, cache_parameters)
+    QuadraticSpline(u, t, I, p, k, c, sc, extrapolation_left,
+        extrapolation_right, cache_parameters, assume_linear_t)
+end
+
 """
     CubicSpline(u, t; extrapolation::ExtrapolationType.T = ExtrapolationType.None, extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
         extrapolation_right::ExtrapolationType.T = ExtrapolationType.None, cache_parameters = false)
@@ -814,15 +901,21 @@ function BSplineInterpolation(
     u, t = munge_data(u, t)
     n = length(t)
     n < d + 1 && error("BSplineInterpolation needs at least d + 1, i.e. $(d+1) points.")
-    s = zero(eltype(u))
+    s = zero(eltype(t))
     p = zero(t)
     k = zeros(eltype(t), n + d + 1)
-    l = zeros(eltype(u), n - 1)
+    l = zeros(eltype(t), n - 1)
     p[1] = zero(eltype(t))
     p[end] = one(eltype(t))
 
     for i in 2:n
-        s += hypot(t[i] - t[i - 1], u[i] - u[i - 1])
+        if eltype(u) <: AbstractVector
+            # For Vector of Vectors, use norm for the vector difference
+            s += sqrt((t[i] - t[i - 1])^2 + norm(u[i] - u[i - 1])^2)
+        else
+            # For scalar values, use hypot as before
+            s += hypot(t[i] - t[i - 1], u[i] - u[i - 1])
+        end
         l[i - 1] = s
     end
     if pVecType == :Uniform
@@ -1051,15 +1144,21 @@ function BSplineApprox(
     u, t = munge_data(u, t)
     n = length(t)
     h < d + 1 && error("BSplineApprox needs at least d + 1, i.e. $(d+1) control points.")
-    s = zero(eltype(u))
+    s = zero(eltype(t))
     p = zero(t)
     k = zeros(eltype(t), h + d + 1)
-    l = zeros(eltype(u), n - 1)
+    l = zeros(eltype(t), n - 1)
     p[1] = zero(eltype(t))
     p[end] = one(eltype(t))
 
     for i in 2:n
-        s += hypot(t[i] - t[i - 1], u[i] - u[i - 1])
+        if eltype(u) <: AbstractVector
+            # For Vector of Vectors, use norm for the vector difference
+            s += sqrt((t[i] - t[i - 1])^2 + norm(u[i] - u[i - 1])^2)
+        else
+            # For scalar values, use hypot as before
+            s += hypot(t[i] - t[i - 1], u[i] - u[i - 1])
+        end
         l[i - 1] = s
     end
     if pVecType == :Uniform
@@ -1108,10 +1207,10 @@ function BSplineApprox(
         end
     end
     # control points
-    c = zeros(eltype(u), h)
+    c = [zero(first(u)) for _ in 1:h]
     c[1] = u[1]
     c[end] = u[end]
-    q = zeros(eltype(u), n)
+    q = [zero(first(u)) for _ in 1:n]
     sc = zeros(eltype(t), n, h)
     for i in 1:n
         spline_coefficients!(view(sc, i, :), d, k, p[i])

src/interpolation_methods.jl

@@ -204,6 +204,12 @@ function _interpolate(A::AkimaInterpolation{<:AbstractVector}, t::Number, iguess
     @evalpoly wj A.u[idx] A.b[idx] A.c[idx] A.d[idx]
 end
 
+function _interpolate(A::AkimaInterpolation{<:AbstractMatrix}, t::Number, iguess)
+    idx = get_idx(A, t, iguess)
+    wj = t - A.t[idx]
+    @evalpoly wj A.u[:, idx] A.b[:, idx] A.c[:, idx] A.d[:, idx]
+end
+
 # Constant Interpolation
 function _interpolate(A::ConstantInterpolation{<:AbstractVector}, t::Number, iguess)
     if A.dir === :left
@@ -252,6 +258,31 @@ function _interpolate(A::QuadraticSpline{<:AbstractVector}, t::Number, iguess)
     Δt_scaled * (α * Δt_scaled + β) + uᵢ
 end
 
+function _interpolate(A::QuadraticSpline{<:AbstractMatrix}, t::Number, iguess)
+    idx = get_idx(A, t, iguess)
+    uᵢ = A.u[:, idx]
+    Δt_scaled = (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx])
+
+    # Compute parameters for each row (always compute on the fly for matrices)
+    result = similar(uᵢ)
+    for row in 1:size(A.u, 1)
+        # Compute parameters on the fly for this row
+        tᵢ₊ = (A.t[idx] + A.t[idx + 1]) / 2
+        sc = zeros(eltype(A.t), length(A.t))
+        nonzero_coefficient_idxs = spline_coefficients!(sc, 2, A.k, tᵢ₊)
+        uᵢ₊ = zero(eltype(A.u))
+        for j in nonzero_coefficient_idxs
+            uᵢ₊ += sc[j] * A.c[row, j]
+        end
+        α_row = 2 * (A.u[row, idx + 1] + A.u[row, idx]) - 4uᵢ₊
+        β_row = 4 * (uᵢ₊ - A.u[row, idx]) - (A.u[row, idx + 1] - A.u[row, idx])
+
+        result[row] = Δt_scaled * (α_row * Δt_scaled + β_row) + A.u[row, idx]
+    end
+
+    result
+end
+
 # CubicSpline Interpolation
 function _interpolate(A::CubicSpline{<:AbstractVector}, t::Number, iguess)
     idx = get_idx(A, t, iguess)

test/interpolation_tests.jl

@@ -478,6 +478,78 @@ end
     A = @inferred(AkimaInterpolation(u, t))
     @test_throws DataInterpolations.LeftExtrapolationError A(-1.0)
     @test_throws DataInterpolations.RightExtrapolationError A(11.0)
+
+    # Test AbstractMatrix interpolation
+    @testset "AbstractMatrix" begin
+        t = 0.0:0.2:1.0
+        u_matrix = [sin.(t) cos.(t)]'  # 2x6 matrix
+        A = AkimaInterpolation(u_matrix, t; extrapolation = ExtrapolationType.Extension)
+
+        # Test interpolation at original points
+        for (i, ti) in enumerate(t)
+            expected = u_matrix[:, i]
+            @test A(ti) ≈ expected
+        end
+
+        # Test interpolation at intermediate points
+        result_half = A(0.5)
+        @test length(result_half) == 2
+        @test result_half[1] ≈ sin(0.5) atol=0.1  # Akima may not be exact
+        @test result_half[2] ≈ cos(0.5) atol=0.1
+
+        # Test extrapolation
+        result_extrap = A(-0.1)
+        @test length(result_extrap) == 2
+
+        # Test output dimensions
+        @test @inferred(output_dim(A)) == 1
+        @test @inferred(output_size(A)) == (2,)
+
+        # Test edge cases
+        # Single row matrix
+        u_single_row = reshape([0.0, 1.0, 2.0, 3.0, 4.0], 1, :)
+        A_single = AkimaInterpolation(u_single_row, collect(0:4))
+        result_single = A_single(2.5)
+        @test length(result_single) == 1
+        @test @inferred(output_dim(A_single)) == 1
+        @test @inferred(output_size(A_single)) == (1,)
+
+        # Larger matrix (3x5)
+        t_large = 0.0:0.25:1.0
+        u_large = [sin.(t_large) cos.(t_large) (t_large.^2)]'  # 3x5 matrix  
+        A_large = AkimaInterpolation(u_large, t_large)
+        result_large = A_large(0.6)
+        @test length(result_large) == 3
+        @test @inferred(output_dim(A_large)) == 1
+        @test @inferred(output_size(A_large)) == (3,)
+    end
+
+    # Test Vector of Vectors support  
+    @testset "Vector of Vectors" begin
+        t = 0.0:0.2:1.0
+        u_vec_of_vec = [[sin(ti), cos(ti)] for ti in t]
+        A = AkimaInterpolation(u_vec_of_vec, t; extrapolation = ExtrapolationType.Extension)
+
+        # Test interpolation at original points
+        for (i, ti) in enumerate(t)
+            expected = u_vec_of_vec[i]
+            @test A(ti) ≈ expected
+        end
+
+        # Test interpolation at intermediate points
+        result_half = A(0.5)
+        @test length(result_half) == 2
+        @test result_half[1] ≈ sin(0.5) atol=0.1
+        @test result_half[2] ≈ cos(0.5) atol=0.1
+
+        # Test extrapolation
+        result_extrap = A(-0.1)
+        @test length(result_extrap) == 2
+
+        # Test output dimensions
+        @test @inferred(output_dim(A)) == 1
+        @test @inferred(output_size(A)) == (2,)
+    end
 end
 
 @testset "ConstantInterpolation" begin
@@ -746,6 +818,51 @@ end
     @test @inferred(output_dim(A)) == 2
     @test @inferred(output_size(A)) == (4, 3)
 
+    # Test AbstractMatrix interpolation
+    @testset "AbstractMatrix" begin
+        t = [-1.0, 0.0, 1.0]
+        u_matrix = [0.0 1.0 3.0; 2.0 3.0 5.0]  # 2x3 matrix
+        A = QuadraticSpline(u_matrix, t; extrapolation = ExtrapolationType.Extension)
+
+        # Test interpolation at original points
+        for (i, ti) in enumerate(t)
+            expected = u_matrix[:, i]
+            @test A(ti) ≈ expected
+        end
+
+        # Test interpolation at intermediate points
+        result_half = A(0.5)
+        @test length(result_half) == 2
+
+        # For the first row: uses P₁(x) = 0.5 * (x + 1) * (x + 2) with points [0.0, 1.0, 3.0]
+        # For the second row: same pattern but with points [2.0, 3.0, 5.0]
+        # At x = 0.5: P₁(0.5) = 0.5 * (1.5) * (2.5) = 1.875
+        @test result_half[1] ≈ P₁(0.5)
+
+        # Test extrapolation
+        result_extrap = A(-2.0)
+        @test length(result_extrap) == 2
+        @test result_extrap[1] ≈ P₁(-2.0)
+
+        # Test output dimensions
+        @test @inferred(output_dim(A)) == 1
+        @test @inferred(output_size(A)) == (2,)
+
+        # Test edge cases  
+        # Single row matrix
+        u_single_row = reshape([2.0, 3.0, 5.0], 1, :)
+        A_single = QuadraticSpline(u_single_row, t; extrapolation = ExtrapolationType.Extension)
+        result_single = A_single(0.5)
+        @test length(result_single) == 1
+        @test @inferred(output_dim(A_single)) == 1
+        @test @inferred(output_size(A_single)) == (1,)
+
+        # Test with cache_parameters=true
+        A_cached = QuadraticSpline(u_matrix, t; cache_parameters=true, extrapolation = ExtrapolationType.Extension)
+        result_cached = A_cached(0.5)
+        @test result_cached ≈ A(0.5)
+    end
+
     # Test extrapolation
     u = [0.0, 1.0, 3.0]
     t = [-1.0, 0.0, 1.0]