diff --git a/src/integrals.jl b/src/integrals.jl index f0e4afed..f92ecacb 100644 --- a/src/integrals.jl +++ b/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 diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index f88ca382..a7801424 100644 --- a/src/interpolation_caches.jl +++ b/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]) diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index 847c3a8c..e0ee2d9d 100644 --- a/src/interpolation_methods.jl +++ b/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) diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl index 4ea6cb26..0b4dd12b 100644 --- a/test/interpolation_tests.jl +++ b/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]