diff --git a/src/derivatives.jl b/src/derivatives.jl index 268e3ec0..d8e5b11d 100644 --- a/src/derivatives.jl +++ b/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 diff --git a/src/integrals.jl b/src/integrals.jl index 92255892..799db50d 100644 --- a/src/integrals.jl +++ b/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 diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index 5c4a9f4a..5b58bd9b 100644 --- a/src/interpolation_caches.jl +++ b/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) diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index 47e18f26..b19f89de 100644 --- a/src/interpolation_methods.jl +++ b/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 diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl index 8a1c7f6f..6b5c443f 100644 --- a/test/interpolation_tests.jl +++ b/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