From 51d3dee5343f67e59a86e55bf988dbdc502a2efe Mon Sep 17 00:00:00 2001 From: Ashutosh Bharambe Date: Mon, 7 Oct 2024 08:23:57 +0000 Subject: [PATCH 1/8] feat(CubicSpline): add `parameter` dispatch for AbstractMatrix --- src/parameter_caches.jl | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/parameter_caches.jl b/src/parameter_caches.jl index 824b2e2a..6f3826ef 100644 --- a/src/parameter_caches.jl +++ b/src/parameter_caches.jl @@ -113,12 +113,19 @@ function CubicSplineParameterCache(u, h, z, cache_parameters) end end -function cubic_spline_parameters(u, h, z, idx) +function cubic_spline_parameters(u::AbstractVector, h, z, idx) c₁ = (u[idx + 1] / h[idx + 1] - z[idx + 1] * h[idx + 1] / 6) c₂ = (u[idx] / h[idx + 1] - z[idx] * h[idx + 1] / 6) return c₁, c₂ end +function cubic_spline_parameters(u::AbstractArray, h, z, idx) + ax = axes(u)[1:end-1] + c₁ = (u[ax..., idx + 1] / h[idx + 1] - z[ax..., idx + 1] * h[idx + 1] / 6) + c₂ = (u[ax..., idx] / h[idx + 1] - z[ax..., idx] * h[idx + 1] / 6) + return c₁, c₂ +end + struct CubicHermiteParameterCache{pType} c₁::pType c₂::pType From e9a6d733af2ad36cc0df591ab377a519c6ea1ca1 Mon Sep 17 00:00:00 2001 From: Ashutosh Bharambe Date: Mon, 7 Oct 2024 08:24:39 +0000 Subject: [PATCH 2/8] feat(CubicSpline): add interpolation dispatch for AbstractMatrix --- src/interpolation_methods.jl | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index 03adf558..69159b20 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -166,6 +166,17 @@ function _interpolate(A::CubicSpline{<:AbstractVector}, t::Number, iguess) I + C + D end +function _interpolate(A::CubicSpline{<:AbstractMatrix}, t::Number, iguess) + idx = get_idx(A, t, iguess) + Δt₁ = t - A.t[idx] + Δt₂ = A.t[idx + 1] - t + I = (A.z[:, idx] * Δt₂^3 + A.z[:, idx + 1] * Δt₁^3) / (6A.h[idx + 1]) + c₁, c₂ = get_parameters(A, idx) + C = c₁ * Δt₁ + D = c₂ * Δt₂ + I + C + D +end + # BSpline Curve Interpolation function _interpolate(A::BSplineInterpolation{<:AbstractVector{<:Number}}, t::Number, From e3406e51254f26562190921f772ef79a504661ea Mon Sep 17 00:00:00 2001 From: Ashutosh Bharambe Date: Mon, 7 Oct 2024 08:25:17 +0000 Subject: [PATCH 3/8] feat(CubicSpline): add dispatch for AbstractMatrix --- src/interpolation_caches.jl | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index d7dce336..4059604a 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -460,6 +460,37 @@ function CubicSpline(u::uType, CubicSpline(u, t, I, p, h[1:(n + 1)], z, extrapolate, cache_parameters, linear_lookup) end +function CubicSpline(u::uType, + t; + extrapolate = false, cache_parameters = false, + assume_linear_t = 1e-2) where {uType <: + AbstractMatrix} + u, t = munge_data(u, t) + n = length(t) - 1 + h = vcat(0, map(k -> t[k + 1] - t[k], 1:(length(t) - 1)), 0) + dl = vcat(h[2:n], zero(eltype(h))) + d_tmp = 2 .* (h[1:(n + 1)] .+ h[2:(n + 2)]) + du = vcat(zero(eltype(h)), h[3:(n + 1)]) + tA = Tridiagonal(dl, d_tmp, du) + + # zero for element type of d, which we don't know yet + ax = axes(u)[1:end-1] + typed_zero = zero(6(u[ax..., begin + 2] - u[ax..., begin + 1]) / h[begin + 2] - + 6(u[ax..., begin + 1] - u[ax..., begin]) / h[begin + 1]) + + h_ = reshape(h, 1, :) + d = 6*((u[ax..., 3:n+1] - u[ax..., 2:n]) ./ h_[:, 3:n+1]) - 6*((u[ax..., 2:n] - u[ax..., 1:n-1]) ./ h_[:, 2:n]) + d = cat(typed_zero, d, typed_zero; dims = ndims(d)) + + z = (tA \ d')' + linear_lookup = seems_linear(assume_linear_t, t) + p = CubicSplineParameterCache(u, h, z, cache_parameters) + A = CubicSpline( + u, t, nothing, p, h[1:(n + 1)], z, extrapolate, cache_parameters, linear_lookup) + I = cumulative_integral(A, cache_parameters) + CubicSpline(u, t, I, p, h[1:(n + 1)], z, extrapolate, cache_parameters, linear_lookup) +end + function CubicSpline( u::uType, t; extrapolate = false, cache_parameters = false, assume_linear_t = 1e-2) where {uType <: @@ -486,6 +517,8 @@ function CubicSpline( CubicSpline(u, t, I, p, h[1:(n + 1)], z, extrapolate, cache_parameters, assume_linear_t) end + + """ BSplineInterpolation(u, t, d, pVecType, knotVecType; extrapolate = false, safetycopy = true) From b9637bdc719db419670d53188ef60ebbedd5b90a Mon Sep 17 00:00:00 2001 From: Ashutosh Bharambe Date: Mon, 7 Oct 2024 08:30:03 +0000 Subject: [PATCH 4/8] test(CubicSpline): add tests for `AbstractMatrix` --- test/interpolation_tests.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl index 69e5c197..97e2a8d4 100644 --- a/test/interpolation_tests.jl +++ b/test/interpolation_tests.jl @@ -599,6 +599,16 @@ end A = CubicSpline(u, t) @test_throws DataInterpolations.ExtrapolationError A(-2.0) @test_throws DataInterpolations.ExtrapolationError A(2.0) + + @testset "AbstractMatrix" begin + t = 0.1:0.1:1.0 + u = [sin.(t) cos.(t)]' |> collect + c = CubicSpline(u, t) + t_test = 0.1:0.05:1.0 + u_test = reduce(hcat, c.(t_test)) + @test isapprox(u_test[1, :], sin.(t_test), atol = 1e-3) + @test isapprox(u_test[2, :], cos.(t_test), atol = 1e-3) + end end @testset "BSplines" begin From 2b9266217549eae3735757cf21618295b1276627 Mon Sep 17 00:00:00 2001 From: Ashutosh Bharambe Date: Mon, 7 Oct 2024 08:41:07 +0000 Subject: [PATCH 5/8] chore(DataInterpolations): fix formatting --- src/interpolation_caches.jl | 17 ++++++++--------- src/parameter_caches.jl | 2 +- test/interpolation_tests.jl | 6 +++--- 3 files changed, 12 insertions(+), 13 deletions(-) diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index 4059604a..0baa65ad 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -461,10 +461,10 @@ function CubicSpline(u::uType, end function CubicSpline(u::uType, - t; - extrapolate = false, cache_parameters = false, - assume_linear_t = 1e-2) where {uType <: - AbstractMatrix} + t; + extrapolate = false, cache_parameters = false, + assume_linear_t = 1e-2) where {uType <: + AbstractMatrix} u, t = munge_data(u, t) n = length(t) - 1 h = vcat(0, map(k -> t[k + 1] - t[k], 1:(length(t) - 1)), 0) @@ -474,12 +474,13 @@ function CubicSpline(u::uType, tA = Tridiagonal(dl, d_tmp, du) # zero for element type of d, which we don't know yet - ax = axes(u)[1:end-1] + ax = axes(u)[1:(end - 1)] typed_zero = zero(6(u[ax..., begin + 2] - u[ax..., begin + 1]) / h[begin + 2] - - 6(u[ax..., begin + 1] - u[ax..., begin]) / h[begin + 1]) + 6(u[ax..., begin + 1] - u[ax..., begin]) / h[begin + 1]) h_ = reshape(h, 1, :) - d = 6*((u[ax..., 3:n+1] - u[ax..., 2:n]) ./ h_[:, 3:n+1]) - 6*((u[ax..., 2:n] - u[ax..., 1:n-1]) ./ h_[:, 2:n]) + d = 6 * ((u[ax..., 3:(n + 1)] - u[ax..., 2:n]) ./ h_[:, 3:(n + 1)]) - + 6 * ((u[ax..., 2:n] - u[ax..., 1:(n - 1)]) ./ h_[:, 2:n]) d = cat(typed_zero, d, typed_zero; dims = ndims(d)) z = (tA \ d')' @@ -517,8 +518,6 @@ function CubicSpline( CubicSpline(u, t, I, p, h[1:(n + 1)], z, extrapolate, cache_parameters, assume_linear_t) end - - """ BSplineInterpolation(u, t, d, pVecType, knotVecType; extrapolate = false, safetycopy = true) diff --git a/src/parameter_caches.jl b/src/parameter_caches.jl index 6f3826ef..09362e58 100644 --- a/src/parameter_caches.jl +++ b/src/parameter_caches.jl @@ -120,7 +120,7 @@ function cubic_spline_parameters(u::AbstractVector, h, z, idx) end function cubic_spline_parameters(u::AbstractArray, h, z, idx) - ax = axes(u)[1:end-1] + ax = axes(u)[1:(end - 1)] c₁ = (u[ax..., idx + 1] / h[idx + 1] - z[ax..., idx + 1] * h[idx + 1] / 6) c₂ = (u[ax..., idx] / h[idx + 1] - z[ax..., idx] * h[idx + 1] / 6) return c₁, c₂ diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl index 97e2a8d4..861f0aa7 100644 --- a/test/interpolation_tests.jl +++ b/test/interpolation_tests.jl @@ -600,11 +600,11 @@ end @test_throws DataInterpolations.ExtrapolationError A(-2.0) @test_throws DataInterpolations.ExtrapolationError A(2.0) - @testset "AbstractMatrix" begin - t = 0.1:0.1:1.0 + @testset "AbstractMatrix" begin + t = 0.1:0.1:1.0 u = [sin.(t) cos.(t)]' |> collect c = CubicSpline(u, t) - t_test = 0.1:0.05:1.0 + t_test = 0.1:0.05:1.0 u_test = reduce(hcat, c.(t_test)) @test isapprox(u_test[1, :], sin.(t_test), atol = 1e-3) @test isapprox(u_test[2, :], cos.(t_test), atol = 1e-3) From d6a944270bdba72f452b281d41946e3c15c2401b Mon Sep 17 00:00:00 2001 From: Ashutosh Bharambe Date: Sun, 13 Oct 2024 18:11:49 +0000 Subject: [PATCH 6/8] feat(CubicSpline): update AbstractMatrix codepath for AbstractArrays --- src/interpolation_caches.jl | 14 ++++++++------ src/interpolation_methods.jl | 5 +++-- 2 files changed, 11 insertions(+), 8 deletions(-) diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index 0baa65ad..1462a46f 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -464,7 +464,7 @@ function CubicSpline(u::uType, t; extrapolate = false, cache_parameters = false, assume_linear_t = 1e-2) where {uType <: - AbstractMatrix} + AbstractArray{T, N}} where {T, N} u, t = munge_data(u, t) n = length(t) - 1 h = vcat(0, map(k -> t[k + 1] - t[k], 1:(length(t) - 1)), 0) @@ -478,12 +478,14 @@ function CubicSpline(u::uType, typed_zero = zero(6(u[ax..., begin + 2] - u[ax..., begin + 1]) / h[begin + 2] - 6(u[ax..., begin + 1] - u[ax..., begin]) / h[begin + 1]) - h_ = reshape(h, 1, :) - d = 6 * ((u[ax..., 3:(n + 1)] - u[ax..., 2:n]) ./ h_[:, 3:(n + 1)]) - - 6 * ((u[ax..., 2:n] - u[ax..., 1:(n - 1)]) ./ h_[:, 2:n]) + h_ = reshape(h, ones(Int64, N - 1)..., :) + ax_h = axes(h_)[1:end-1] + d = 6 * ((u[ax..., 3:(n + 1)] - u[ax..., 2:n]) ./ h_[ax_h..., 3:(n + 1)]) - + 6 * ((u[ax..., 2:n] - u[ax..., 1:(n - 1)]) ./ h_[ax_h..., 2:n]) d = cat(typed_zero, d, typed_zero; dims = ndims(d)) - - z = (tA \ d')' + d_reshaped = reshape(d, prod(size(d)[1:end-1]), :) + z = (tA \ d_reshaped')' + z = reshape(z, size(u)...) linear_lookup = seems_linear(assume_linear_t, t) p = CubicSplineParameterCache(u, h, z, cache_parameters) A = CubicSpline( diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index 69159b20..eb018064 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -166,11 +166,12 @@ function _interpolate(A::CubicSpline{<:AbstractVector}, t::Number, iguess) I + C + D end -function _interpolate(A::CubicSpline{<:AbstractMatrix}, t::Number, iguess) +function _interpolate(A::CubicSpline{<:AbstractArray{T, N}}, t::Number, iguess) where {T, N} idx = get_idx(A, t, iguess) Δt₁ = t - A.t[idx] Δt₂ = A.t[idx + 1] - t - I = (A.z[:, idx] * Δt₂^3 + A.z[:, idx + 1] * Δt₁^3) / (6A.h[idx + 1]) + ax = axes(A.z)[1:end-1] + I = (A.z[ax..., idx] * Δt₂^3 + A.z[ax..., idx + 1] * Δt₁^3) / (6A.h[idx + 1]) c₁, c₂ = get_parameters(A, idx) C = c₁ * Δt₁ D = c₂ * Δt₂ From 6330a49fc95e7910e3b4f0a5dcff5ae0fe31e2f2 Mon Sep 17 00:00:00 2001 From: Ashutosh Bharambe Date: Sun, 13 Oct 2024 18:17:39 +0000 Subject: [PATCH 7/8] test(CubicSpline): add tests for AbstractArray dispatch --- test/interpolation_tests.jl | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl index 861f0aa7..e0a02cbc 100644 --- a/test/interpolation_tests.jl +++ b/test/interpolation_tests.jl @@ -609,6 +609,17 @@ end @test isapprox(u_test[1, :], sin.(t_test), atol = 1e-3) @test isapprox(u_test[2, :], cos.(t_test), atol = 1e-3) end + @testset "AbstractArray{T, 3}" begin + f3d(t) = [sin(t) cos(t); + 0. cos(2t)] + t = 0.1:0.1:1.0 + u3d = f3d.(t) + c = CubicSpline(u3d, t) + t_test = 0.1:0.05:1.0 + u_test = reduce(hcat, c.(t_test)) + f_test = reduce(hcat, f3d.(t_test)) + @test isapprox(u_test, f_test, atol = 1e-2) + end end @testset "BSplines" begin From 6588747e34e6a109fc0305e9eca3327368d35330 Mon Sep 17 00:00:00 2001 From: Ashutosh Bharambe Date: Sun, 13 Oct 2024 18:20:22 +0000 Subject: [PATCH 8/8] chore(DataInterpolations): fix formatting --- src/interpolation_caches.jl | 4 ++-- src/interpolation_methods.jl | 2 +- test/interpolation_tests.jl | 4 ++-- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index 1462a46f..e5825ffc 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -479,11 +479,11 @@ function CubicSpline(u::uType, 6(u[ax..., begin + 1] - u[ax..., begin]) / h[begin + 1]) h_ = reshape(h, ones(Int64, N - 1)..., :) - ax_h = axes(h_)[1:end-1] + ax_h = axes(h_)[1:(end - 1)] d = 6 * ((u[ax..., 3:(n + 1)] - u[ax..., 2:n]) ./ h_[ax_h..., 3:(n + 1)]) - 6 * ((u[ax..., 2:n] - u[ax..., 1:(n - 1)]) ./ h_[ax_h..., 2:n]) d = cat(typed_zero, d, typed_zero; dims = ndims(d)) - d_reshaped = reshape(d, prod(size(d)[1:end-1]), :) + d_reshaped = reshape(d, prod(size(d)[1:(end - 1)]), :) z = (tA \ d_reshaped')' z = reshape(z, size(u)...) linear_lookup = seems_linear(assume_linear_t, t) diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index eb018064..4b0ef3e1 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -170,7 +170,7 @@ function _interpolate(A::CubicSpline{<:AbstractArray{T, N}}, t::Number, iguess) idx = get_idx(A, t, iguess) Δt₁ = t - A.t[idx] Δt₂ = A.t[idx + 1] - t - ax = axes(A.z)[1:end-1] + ax = axes(A.z)[1:(end - 1)] I = (A.z[ax..., idx] * Δt₂^3 + A.z[ax..., idx + 1] * Δt₁^3) / (6A.h[idx + 1]) c₁, c₂ = get_parameters(A, idx) C = c₁ * Δt₁ diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl index e0a02cbc..9ee6e5ef 100644 --- a/test/interpolation_tests.jl +++ b/test/interpolation_tests.jl @@ -609,9 +609,9 @@ end @test isapprox(u_test[1, :], sin.(t_test), atol = 1e-3) @test isapprox(u_test[2, :], cos.(t_test), atol = 1e-3) end - @testset "AbstractArray{T, 3}" begin + @testset "AbstractArray{T, 3}" begin f3d(t) = [sin(t) cos(t); - 0. cos(2t)] + 0.0 cos(2t)] t = 0.1:0.1:1.0 u3d = f3d.(t) c = CubicSpline(u3d, t)