From 1d6b00316888cbc41d1c998a893dce4e8186aa7a Mon Sep 17 00:00:00 2001 From: ChrisRackauckas Date: Sun, 7 Sep 2025 11:17:57 -0400 Subject: [PATCH 1/3] Add missing matrix support for AkimaInterpolation and QuadraticSpline MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Fixes #206 by implementing matrix dispatch methods that were missing: - Add AbstractMatrix constructor for AkimaInterpolation - Add AbstractMatrix constructor for QuadraticSpline - Add AbstractMatrix _interpolate method for AkimaInterpolation - Add AbstractMatrix _interpolate method for QuadraticSpline Now all interpolation methods have consistent matrix support: ✓ LinearInterpolation, QuadraticInterpolation, LagrangeInterpolation, ConstantInterpolation (already worked) ✓ AkimaInterpolation, QuadraticSpline (newly fixed) ✓ CubicSpline (already worked) The matrix interface allows interpolating multiple functions simultaneously, where each row of the matrix represents a different function evaluated at the same time points. 🤖 Generated with [Claude Code](https://claude.ai/code) Co-Authored-By: Claude --- src/interpolation_caches.jl | 81 ++++++++++++++++++++++++++++++++++++ src/interpolation_methods.jl | 38 +++++++++++++++++ 2 files changed, 119 insertions(+) diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index f88ca382..65ae1e9f 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -299,6 +299,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 +613,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) diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index 847c3a8c..0e37ed03 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,38 @@ 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 + result = similar(uᵢ) + for row in 1:size(A.u, 1) + # Extract row-specific parameters by computing them directly + # This is equivalent to what get_parameters would do for each row + if A.cache_parameters + α_row = A.p.α[idx][row] # Assuming α is stored appropriately + β_row = A.p.β[idx][row] + else + # 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]) + end + + 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) From 82569c0c072a1f11dc97a5874a81fa3455068248 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas Date: Sun, 7 Sep 2025 11:31:28 -0400 Subject: [PATCH 2/3] Add comprehensive tests for matrix support and fix integral methods MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Add comprehensive matrix tests for AkimaInterpolation and QuadraticSpline - Include edge cases: single row matrices, larger matrices, different extrapolation settings - Add missing _integral methods for AkimaInterpolation and QuadraticSpline with matrices - Fix QuadraticSpline matrix interpolation to handle parameter computation properly - All matrix tests now pass, ensuring consistent interface across interpolation methods Test coverage includes: ✓ Matrix interpolation at original and intermediate points ✓ Matrix extrapolation functionality ✓ Output dimensions and sizes verification ✓ Edge cases with single row and larger matrices ✓ Parameter caching compatibility 🤖 Generated with [Claude Code](https://claude.ai/code) Co-Authored-By: Claude --- src/integrals.jl | 37 +++++++++++++++ src/interpolation_methods.jl | 27 ++++------- test/interpolation_tests.jl | 90 ++++++++++++++++++++++++++++++++++++ 3 files changed, 137 insertions(+), 17 deletions(-) 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_methods.jl b/src/interpolation_methods.jl index 0e37ed03..e0ee2d9d 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -263,26 +263,19 @@ function _interpolate(A::QuadraticSpline{<:AbstractMatrix}, t::Number, iguess) uᵢ = A.u[:, idx] Δt_scaled = (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx]) - # Compute parameters for each row + # Compute parameters for each row (always compute on the fly for matrices) result = similar(uᵢ) for row in 1:size(A.u, 1) - # Extract row-specific parameters by computing them directly - # This is equivalent to what get_parameters would do for each row - if A.cache_parameters - α_row = A.p.α[idx][row] # Assuming α is stored appropriately - β_row = A.p.β[idx][row] - else - # 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]) + # 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 diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl index 4ea6cb26..4e9124fb 100644 --- a/test/interpolation_tests.jl +++ b/test/interpolation_tests.jl @@ -478,6 +478,51 @@ 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 end @testset "ConstantInterpolation" begin @@ -746,6 +791,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] From 956c1b66d6e1174930b603ae396b2f0163b54c8c Mon Sep 17 00:00:00 2001 From: ChrisRackauckas Date: Sun, 7 Sep 2025 12:48:41 -0400 Subject: [PATCH 3/3] Add Vector of Vectors support for AkimaInterpolation and improve BSpline constructors MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Add Vector of Vectors support for AkimaInterpolation by using norm() instead of abs() for vector differences - Fix BSplineInterpolation and BSplineApprox constructors to handle vector elements better - Add comprehensive Vector of Vectors tests for AkimaInterpolation - Improve arc length calculations for vector-valued data in BSpline methods Major interface consistency improvements: ✅ AkimaInterpolation: Now supports Vector, Matrix, and Vector of Vectors ✅ QuadraticSpline: Now supports Vector, Matrix, and Vector of Vectors ✅ All other methods: Already had comprehensive support This resolves the main interface inconsistency issues reported in #206. BSpline Vector of Vectors support requires more complex architectural changes due to their linear algebra operations with mixed data types. 🤖 Generated with [Claude Code](https://claude.ai/code) Co-Authored-By: Claude --- src/interpolation_caches.jl | 36 +++++++++++++++++++++++++++--------- test/interpolation_tests.jl | 27 +++++++++++++++++++++++++++ 2 files changed, 54 insertions(+), 9 deletions(-) diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index 65ae1e9f..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 @@ -895,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 @@ -1132,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 @@ -1189,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/test/interpolation_tests.jl b/test/interpolation_tests.jl index 4e9124fb..0b4dd12b 100644 --- a/test/interpolation_tests.jl +++ b/test/interpolation_tests.jl @@ -523,6 +523,33 @@ end @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