Add matrix support for AkimaInterpolation

src/derivatives.jl

@@ -215,6 +215,13 @@ function _derivative(A::AkimaInterpolation{<:AbstractVector}, t::Number, iguess)
     return @evalpoly wj A.b[idx] 2A.c[j] 3A.d[j]
 end
 
+function _derivative(A::AkimaInterpolation{<:AbstractMatrix}, t::Number, iguess)
+    idx = get_idx(A, t, iguess; idx_shift = -1, side = :first)
+    j = min(idx, size(A.c, 2))  # for smooth derivative at A.t[end]
+    wj = t - A.t[idx]
+    return @evalpoly wj A.b[:, idx] 2A.c[:, j] 3A.d[:, j]
+end
+
 function _derivative(A::ConstantInterpolation, t::Number, iguess)
     return zero(first(A.u))
 end

src/integrals.jl

@@ -284,6 +284,15 @@ function _integral(
     return 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{<:Number}},
+        idx::Number, t1::Number, t2::Number
+    )
+    return integrate_cubic_polynomial(
+        t1, t2, A.t[idx], A.u[:, idx], A.b[:, idx], A.c[:, idx], A.d[:, idx]
+    )
+end
+
 function _integral(A::LagrangeInterpolation, idx::Number, t1::Number, t2::Number)
     throw(IntegralNotFoundError())
 end

src/interpolation_caches.jl

@@ -284,9 +284,11 @@ struct AkimaInterpolation{uType, tType, IType, bType, cType, dType, T} <:
 end
 
 function AkimaInterpolation(
-        u, t; extrapolation::ExtrapolationType.T = ExtrapolationType.None,
+        u::AbstractVector, t;
+        extrapolation::ExtrapolationType.T = ExtrapolationType.None,
         extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
-        extrapolation_right::ExtrapolationType.T = ExtrapolationType.None, cache_parameters = false, assume_linear_t = 1.0e-2
+        extrapolation_right::ExtrapolationType.T = ExtrapolationType.None,
+        cache_parameters = false, assume_linear_t = 1.0e-2
     )
     extrapolation_left,
         extrapolation_right = munge_extrapolation(
@@ -327,6 +329,53 @@ function AkimaInterpolation(
     )
 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 = 1.0e-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)
+    dt = reshape(dt, 1, :)
+    m = Array{eltype(u)}(undef, size(u, 1), n + 3)
+    m[:, 3:(end - 2)] = diff(u; dims = 2) ./ 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 = (m[:, 4:end] .+ m[:, 1:(end - 3)]) ./ 2
+    dm = abs.(diff(m; dims = 2))
+    f1 = dm[:, 3:(n + 2)]
+    f2 = dm[:, 1:n]
+    f12 = f1 + f2
+    ind = findall(f12 .> 1.0e-9 .* maximum(f12; dims = 2))
+    for i in ind
+        row, col = Tuple(i)
+        b[i] = (f1[i] * m[row, col + 1] + f2[i] * m[row, col + 2]) / f12[i]
+    end
+    c = (3 .* m[:, 3:(end - 2)] .- 2 .* b[:, 1:(end - 1)] .- b[:, 2:end]) ./ dt
+    d = (b[:, 1:(end - 1)] .+ b[:, 2:end] .- 2 .* m[:, 3:(end - 2)]) ./ dt .^ 2
+
+    A = AkimaInterpolation(
+        u, t, nothing, b, c, d, extrapolation_left,
+        extrapolation_right, cache_parameters, linear_lookup
+    )
+    I = cumulative_integral(A, cache_parameters)
+    return 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)

src/interpolation_methods.jl

@@ -205,6 +205,12 @@ function _interpolate(A::AkimaInterpolation{<:AbstractVector}, t::Number, iguess
     return @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]
+    return @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

test/interpolation_tests.jl

@@ -504,6 +504,27 @@ end
     @test @inferred(output_dim(A)) == 0
     @test @inferred(output_size(A)) == ()
 
+    # Matrix interpolation test (against Vector version)
+    u_matrix = [permutedims(u); permutedims(2u)]
+    A_matrix = @inferred(AkimaInterpolation(u_matrix, t))
+    A_row_1 = AkimaInterpolation(vec(u_matrix[1, :]), t)
+    A_row_2 = AkimaInterpolation(vec(u_matrix[2, :]), t)
+
+    for (_t, _u) in zip(t, eachcol(u_matrix))
+        @test A_matrix(_t) == _u
+    end
+    for _t in (0.5, 1.5, 5.1, 9.9)
+        @test A_matrix(_t) ≈ [A_row_1(_t), A_row_2(_t)]
+        @test DataInterpolations.derivative(A_matrix, _t) ≈ [
+            DataInterpolations.derivative(A_row_1, _t),
+            DataInterpolations.derivative(A_row_2, _t),
+        ]
+        @test DataInterpolations.derivative(A_matrix, _t, 2) ≈ [
+            DataInterpolations.derivative(A_row_1, _t, 2),
+            DataInterpolations.derivative(A_row_2, _t, 2),
+        ]
+    end
+
     # Test extrapolation
     A = @inferred(AkimaInterpolation(u, t; extrapolation = ExtrapolationType.Extension))
     @test A(-1.0) ≈ -5.0