From caf0b764dd4c1b4cfc961ea263785cb9aee3ede7 Mon Sep 17 00:00:00 2001 From: Sathvik Bhagavan Date: Thu, 17 Oct 2024 22:59:40 +0200 Subject: [PATCH 1/5] refactor: methods for higher dim arrays --- src/interpolation_caches.jl | 27 +++++++++++++++++++++++++++ src/interpolation_methods.jl | 27 ++++++++++++++++++++------- src/parameter_caches.jl | 19 +++++++++++++------ 3 files changed, 60 insertions(+), 13 deletions(-) diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index ca8e7b36..be7ce1ca 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -389,6 +389,33 @@ function QuadraticSpline( QuadraticSpline(u, t, I, p, tA, d, z, extrapolate, cache_parameters, linear_lookup) end +function QuadraticSpline( + u::uType, t; extrapolate = false, cache_parameters = false, + assume_linear_t = 1e-2) where {uType <: + AbstractArray{T, N}} where {T, N} + u, t = munge_data(u, t) + linear_lookup = seems_linear(assume_linear_t, t) + s = length(t) + dl = ones(eltype(t), s - 1) + d_tmp = ones(eltype(t), s) + du = zeros(eltype(t), s - 1) + tA = Tridiagonal(dl, d_tmp, du) + ax = axes(u)[1:(end - 1)] + d_ = map( + i -> i == 1 ? zeros(eltype(t), size(u[ax..., 1])) : + 2 // 1 * (u[ax..., i] - u[ax..., i - 1]) / (t[i] - t[i - 1]), + 1:s) + d = transpose(reshape(reduce(hcat, d_), :, s)) + z_ = reshape(transpose(tA \ d), size(u[ax..., 1])..., :) + z = [z_s for z_s in eachslice(z_, dims = ndims(z_))] + + p = QuadraticSplineParameterCache(z, t, cache_parameters) + A = QuadraticSpline( + u, t, nothing, p, tA, d, z, extrapolate, cache_parameters, linear_lookup) + I = cumulative_integral(A, cache_parameters) + QuadraticSpline(u, t, I, p, tA, d, z, extrapolate, cache_parameters, linear_lookup) +end + """ CubicSpline(u, t; extrapolate = false, cache_parameters = false) diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index 5223d416..240f7244 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -87,13 +87,15 @@ function _interpolate(A::LagrangeInterpolation{<:AbstractVector}, t::Number, igu N / D end -function _interpolate(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number, iguess) +function _interpolate( + A::LagrangeInterpolation{<:AbstractArray{T, N}}, t::Number, iguess) where {T, N} idx = get_idx(A, t, iguess) findRequiredIdxs!(A, t, idx) + ax = axes(A.u)[1:(end - 1)] if A.t[A.idxs[1]] == t - return A.u[:, A.idxs[1]] + return A.u[ax..., A.idxs[1]] end - N = zero(A.u[:, 1]) + N1 = zero(A.u[ax..., 1]) D = zero(A.t[1]) tmp = D for i in 1:length(A.idxs) @@ -111,9 +113,9 @@ function _interpolate(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number, igu end tmp = inv((t - A.t[A.idxs[i]]) * mult) D += tmp - @. N += (tmp * A.u[:, A.idxs[i]]) + @. N1 += (tmp * A.u[ax..., A.idxs[i]]) end - N / D + N1 / D end function _interpolate(A::AkimaInterpolation{<:AbstractVector}, t::Number, iguess) @@ -134,7 +136,8 @@ function _interpolate(A::ConstantInterpolation{<:AbstractVector}, t::Number, igu A.u[idx] end -function _interpolate(A::ConstantInterpolation{<:AbstractMatrix}, t::Number, iguess) +function _interpolate( + A::ConstantInterpolation{<:AbstractArray{T, N}}, t::Number, iguess) where {T, N} if A.dir === :left # :left means that value to the left is used for interpolation idx = get_idx(A, t, iguess; lb = 1, ub_shift = 0) @@ -142,7 +145,7 @@ function _interpolate(A::ConstantInterpolation{<:AbstractMatrix}, t::Number, igu # :right means that value to the right is used for interpolation idx = get_idx(A, t, iguess; side = :first, lb = 1, ub_shift = 0) end - A.u[:, idx] + A.u[axes(A.u)[1:(end - 1)]..., idx] end # QuadraticSpline Interpolation @@ -154,6 +157,16 @@ function _interpolate(A::QuadraticSpline{<:AbstractVector}, t::Number, iguess) return A.z[idx] * Δt + σ * Δt^2 + Cᵢ end +function _interpolate( + A::QuadraticSpline{<:AbstractArray{T, N}}, t::Number, iguess) where {T, N} + idx = get_idx(A, t, iguess) + ax = axes(A.u)[1:(end - 1)] + Cᵢ = A.u[ax..., idx] + Δt = t - A.t[idx] + σ = get_parameters(A, idx) + return A.z[idx] * Δt + σ * Δt^2 + Cᵢ +end + # CubicSpline Interpolation function _interpolate(A::CubicSpline{<:AbstractVector}, t::Number, iguess) idx = get_idx(A, t, iguess) diff --git a/src/parameter_caches.jl b/src/parameter_caches.jl index 7c2d6f0b..424a7b9b 100644 --- a/src/parameter_caches.jl +++ b/src/parameter_caches.jl @@ -52,11 +52,12 @@ function QuadraticParameterCache(u, t, cache_parameters) end end -function quadratic_interpolation_parameters(u, t, idx) - if u isa AbstractMatrix - u₀ = u[:, idx] - u₁ = u[:, idx + 1] - u₂ = u[:, idx + 2] +function quadratic_interpolation_parameters(u::AbstractArray{T, N}, t, idx) where {T, N} + if N > 1 + ax = axes(u) + u₀ = u[ax[1:(end - 1)]..., idx] + u₁ = u[ax[1:(end - 1)]..., idx + 1] + u₂ = u[ax[1:(end - 1)]..., idx + 2] else u₀ = u[idx] u₁ = u[idx + 1] @@ -89,11 +90,17 @@ function QuadraticSplineParameterCache(z, t, cache_parameters) end end -function quadratic_spline_parameters(z, t, idx) +function quadratic_spline_parameters(z::AbstractVector, t, idx) σ = 1 // 2 * (z[idx + 1] - z[idx]) / (t[idx + 1] - t[idx]) return σ end +function quadratic_spline_parameters(z::AbstractArray, t, idx) + ax = axes(z)[1:(end - 1)] + σ = 1 // 2 * (z[ax..., idx + 1] - z[ax..., idx]) / (t[idx + 1] - t[idx]) + return σ +end + struct CubicSplineParameterCache{pType} c₁::pType c₂::pType From e9e52d78cebda7b5e82e89ba721e21deda9b5f4f Mon Sep 17 00:00:00 2001 From: Sathvik Bhagavan Date: Sat, 19 Oct 2024 12:22:52 +0200 Subject: [PATCH 2/5] refactor: interpolation with higher dim arrays --- src/interpolation_caches.jl | 42 ++++------------ src/interpolation_methods.jl | 47 ++++++++++++++--- src/parameter_caches.jl | 97 +++++++++++++++++++++++++++++++++++- 3 files changed, 143 insertions(+), 43 deletions(-) diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index be7ce1ca..226b4f7a 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -173,29 +173,24 @@ Extrapolation extends the last cubic polynomial on each side. for a test based on the normalized standard deviation of the difference with respect to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2. """ -struct AkimaInterpolation{uType, tType, IType, bType, cType, dType, T, N} <: +struct AkimaInterpolation{uType, tType, IType, pType, T, N} <: AbstractInterpolation{T, N} u::uType t::tType I::IType - b::bType - c::cType - d::dType + p::pType extrapolate::Bool iguesser::Guesser{tType} cache_parameters::Bool linear_lookup::Bool function AkimaInterpolation( - u, t, I, b, c, d, extrapolate, cache_parameters, assume_linear_t) + u, t, I, p, extrapolate, cache_parameters, assume_linear_t) linear_lookup = seems_linear(assume_linear_t, t) N = get_output_dim(u) - new{typeof(u), typeof(t), typeof(I), typeof(b), typeof(c), - typeof(d), eltype(u), N}(u, + new{typeof(u), typeof(t), typeof(I), typeof(p), eltype(u), N}(u, t, I, - b, - c, - d, + p, extrapolate, Guesser(t), cache_parameters, @@ -208,30 +203,11 @@ function AkimaInterpolation( u, t; extrapolate = false, cache_parameters = false, assume_linear_t = 1e-2) u, t = munge_data(u, t) linear_lookup = seems_linear(assume_linear_t, t) - n = length(t) - dt = diff(t) - m = Array{eltype(u)}(undef, n + 3) - m[3:(end - 2)] = diff(u) ./ 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 = 0.5 .* (m[4:end] .+ m[1:(end - 3)]) - dm = abs.(diff(m)) - f1 = dm[3:(n + 2)] - f2 = dm[1:n] - f12 = f1 + f2 - ind = findall(f12 .> 1e-9 * maximum(f12)) - b[ind] = (f1[ind] .* m[ind .+ 1] .+ - f2[ind] .* m[ind .+ 2]) ./ f12[ind] - c = (3.0 .* m[3:(end - 2)] .- 2.0 .* b[1:(end - 1)] .- b[2:end]) ./ dt - d = (b[1:(end - 1)] .+ b[2:end] .- 2.0 .* m[3:(end - 2)]) ./ dt .^ 2 - + p = AkimaParameterCache(u, t) A = AkimaInterpolation( - u, t, nothing, b, c, d, extrapolate, cache_parameters, linear_lookup) + u, t, nothing, p, extrapolate, cache_parameters, linear_lookup) I = cumulative_integral(A, cache_parameters) - AkimaInterpolation(u, t, I, b, c, d, extrapolate, cache_parameters, linear_lookup) + AkimaInterpolation(u, t, I, p, extrapolate, cache_parameters, linear_lookup) end """ @@ -392,7 +368,7 @@ end function QuadraticSpline( u::uType, t; extrapolate = false, cache_parameters = false, assume_linear_t = 1e-2) where {uType <: - AbstractArray{T, N}} where {T, N} + AbstractArray} u, t = munge_data(u, t) linear_lookup = seems_linear(assume_linear_t, t) s = length(t) diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index 240f7244..f1a21741 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -88,14 +88,14 @@ function _interpolate(A::LagrangeInterpolation{<:AbstractVector}, t::Number, igu end function _interpolate( - A::LagrangeInterpolation{<:AbstractArray{T, N}}, t::Number, iguess) where {T, N} + A::LagrangeInterpolation{<:AbstractArray}, t::Number, iguess) idx = get_idx(A, t, iguess) findRequiredIdxs!(A, t, idx) ax = axes(A.u)[1:(end - 1)] if A.t[A.idxs[1]] == t return A.u[ax..., A.idxs[1]] end - N1 = zero(A.u[ax..., 1]) + N = zero(A.u[ax..., 1]) D = zero(A.t[1]) tmp = D for i in 1:length(A.idxs) @@ -113,15 +113,22 @@ function _interpolate( end tmp = inv((t - A.t[A.idxs[i]]) * mult) D += tmp - @. N1 += (tmp * A.u[ax..., A.idxs[i]]) + @. N += (tmp * A.u[ax..., A.idxs[i]]) end - N1 / D + N / D end function _interpolate(A::AkimaInterpolation{<:AbstractVector}, 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] + @evalpoly wj A.u[idx] A.p.b[idx] A.p.c[idx] A.p.d[idx] +end + +function _interpolate(A::AkimaInterpolation{<:AbstractArray}, t::Number, iguess) + idx = get_idx(A, t, iguess) + wj = t - A.t[idx] + ax = axes(A.u)[1:(end - 1)] + @. @evalpoly wj A.u[ax..., idx] A.p.b[ax..., idx] A.p.c[ax..., idx] A.p.d[ax..., idx] end # ConstantInterpolation Interpolation @@ -137,7 +144,7 @@ function _interpolate(A::ConstantInterpolation{<:AbstractVector}, t::Number, igu end function _interpolate( - A::ConstantInterpolation{<:AbstractArray{T, N}}, t::Number, iguess) where {T, N} + A::ConstantInterpolation{<:AbstractArray}, t::Number, iguess) if A.dir === :left # :left means that value to the left is used for interpolation idx = get_idx(A, t, iguess; lb = 1, ub_shift = 0) @@ -158,7 +165,7 @@ function _interpolate(A::QuadraticSpline{<:AbstractVector}, t::Number, iguess) end function _interpolate( - A::QuadraticSpline{<:AbstractArray{T, N}}, t::Number, iguess) where {T, N} + A::QuadraticSpline{<:AbstractArray}, t::Number, iguess) idx = get_idx(A, t, iguess) ax = axes(A.u)[1:(end - 1)] Cᵢ = A.u[ax..., idx] @@ -179,7 +186,7 @@ function _interpolate(A::CubicSpline{<:AbstractVector}, t::Number, iguess) I + C + D end -function _interpolate(A::CubicSpline{<:AbstractArray{T, N}}, t::Number, iguess) where {T, N} +function _interpolate(A::CubicSpline{<:AbstractArray}, t::Number, iguess) idx = get_idx(A, t, iguess) Δt₁ = t - A.t[idx] Δt₂ = A.t[idx + 1] - t @@ -238,6 +245,18 @@ function _interpolate( out end +function _interpolate( + A::CubicHermiteSpline{<:AbstractArray}, t::Number, iguess) + idx = get_idx(A, t, iguess) + Δt₀ = t - A.t[idx] + Δt₁ = t - A.t[idx + 1] + ax = axes(A.u)[1:(end - 1)] + out = A.u[ax..., idx] .+ Δt₀ .* A.du[ax..., idx] + c₁, c₂ = get_parameters(A, idx) + out .+= Δt₀^2 .* (c₁ .+ Δt₁ .* c₂) + out +end + # Quintic Hermite Spline function _interpolate( A::QuinticHermiteSpline{<:AbstractVector{<:Number}}, t::Number, iguess) @@ -249,3 +268,15 @@ function _interpolate( out += Δt₀^3 * (c₁ + Δt₁ * (c₂ + c₃ * Δt₁)) out end + +function _interpolate( + A::QuinticHermiteSpline{<:AbstractArray}, t::Number, iguess) + idx = get_idx(A, t, iguess) + Δt₀ = t - A.t[idx] + Δt₁ = t - A.t[idx + 1] + ax = axes(A.u)[1:(end - 1)] + out = A.u[ax..., idx] + Δt₀ * (A.du[ax..., idx] + A.ddu[ax..., idx] * Δt₀ / 2) + c₁, c₂, c₃ = get_parameters(A, idx) + out .+= Δt₀^3 .* (c₁ .+ Δt₁ .* (c₂ .+ c₃ .* Δt₁)) + out +end diff --git a/src/parameter_caches.jl b/src/parameter_caches.jl index 424a7b9b..be6adbdd 100644 --- a/src/parameter_caches.jl +++ b/src/parameter_caches.jl @@ -75,6 +75,72 @@ function quadratic_interpolation_parameters(u::AbstractArray{T, N}, t, idx) wher return l₀, l₁, l₂ end +struct AkimaParameterCache{pType} + b::pType + c::pType + d::pType +end + +function AkimaParameterCache(u, t) + b, c, d = akima_interpolation_parameters(u, t) + AkimaParameterCache(b, c, d) +end + +function akima_interpolation_parameters(u::AbstractVector, t) + n = length(t) + dt = diff(t) + m = Array{eltype(u)}(undef, n + 3) + m[3:(end - 2)] = diff(u) ./ 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 = 0.5 .* (m[4:end] .+ m[1:(end - 3)]) + dm = abs.(diff(m)) + f1 = dm[3:(n + 2)] + f2 = dm[1:n] + f12 = f1 + f2 + ind = findall(f12 .> 1e-9 * maximum(f12)) + b[ind] = (f1[ind] .* m[ind .+ 1] .+ + f2[ind] .* m[ind .+ 2]) ./ f12[ind] + c = (3.0 .* m[3:(end - 2)] .- 2.0 .* b[1:(end - 1)] .- b[2:end]) ./ dt + d = (b[1:(end - 1)] .+ b[2:end] .- 2.0 .* m[3:(end - 2)]) ./ dt .^ 2 + return b, c, d +end + +function akima_interpolation_parameters(u::AbstractArray, t) + n = length(t) + dt = diff(t) + ax = axes(u)[1:(end - 1)] + su = size(u) + m = zeros(eltype(u), su[1:(end - 1)]..., n + 3) + m[ax..., 3:(end - 2)] .= mapslices( + x -> x ./ dt, diff(u, dims = length(su)); dims = length(su)) + m[ax..., 2] .= 2m[ax..., 3] .- m[ax..., 4] + m[ax..., 1] .= 2m[ax..., 2] .- m[3] + m[ax..., end - 1] .= 2m[ax..., end - 2] - m[ax..., end - 3] + m[ax..., end] .= 2m[ax..., end - 1] .- m[ax..., end - 2] + b = 0.5 .* (m[ax..., 4:end] .+ m[ax..., 1:(end - 3)]) + dm = abs.(diff(m, dims = length(su))) + f1 = dm[ax..., 3:(n + 2)] + f2 = dm[ax..., 1:n] + f12 = f1 .+ f2 + ind = findall(f12 .> 1e-9 * maximum(f12)) + indi = map(i -> i.I, ind) + b[ind] .= (f1[ind] .* + m[CartesianIndex.(map(i -> (i[1:(end - 1)]..., i[end] + 1), indi))] .+ + f2[ind] .* + m[CartesianIndex.(map(i -> (i[1:(end - 1)]..., i[end] + 2), indi))]) ./ + f12[ind] + c = mapslices(x -> x ./ dt, + (3.0 .* m[ax..., 3:(end - 2)] .- 2.0 .* b[ax..., 1:(end - 1)] .- b[ax..., 2:end]); + dims = length(su)) + d = mapslices(x -> x ./ dt .^ 2, + (b[ax..., 1:(end - 1)] .+ b[ax..., 2:end] .- 2.0 .* m[ax..., 3:(end - 2)]); + dims = length(su)) + return b, c, d +end + struct QuadraticSplineParameterCache{pType} σ::pType end @@ -152,7 +218,19 @@ function CubicHermiteParameterCache(du, u, t, cache_parameters) end end -function cubic_hermite_spline_parameters(du, u, t, idx) +function cubic_hermite_spline_parameters(du::AbstractArray, u, t, idx) + ax = axes(u)[1:(end - 1)] + Δt = t[idx + 1] - t[idx] + u₀ = u[ax..., idx] + u₁ = u[ax..., idx + 1] + du₀ = du[ax..., idx] + du₁ = du[ax..., idx + 1] + c₁ = (u₁ - u₀ - du₀ * Δt) / Δt^2 + c₂ = (du₁ - du₀ - 2c₁ * Δt) / Δt^2 + return c₁, c₂ +end + +function cubic_hermite_spline_parameters(du::AbstractVector, u, t, idx) Δt = t[idx + 1] - t[idx] u₀ = u[idx] u₁ = u[idx + 1] @@ -183,7 +261,7 @@ function QuinticHermiteParameterCache(ddu, du, u, t, cache_parameters) end end -function quintic_hermite_spline_parameters(ddu, du, u, t, idx) +function quintic_hermite_spline_parameters(ddu::AbstractVector, du, u, t, idx) Δt = t[idx + 1] - t[idx] u₀ = u[idx] u₁ = u[idx + 1] @@ -196,3 +274,18 @@ function quintic_hermite_spline_parameters(ddu, du, u, t, idx) c₃ = (6u₁ - 6u₀ - 3(du₀ + du₁)Δt + (ddu₁ - ddu₀)Δt^2 / 2) / Δt^5 return c₁, c₂, c₃ end + +function quintic_hermite_spline_parameters(ddu::AbstractArray, du, u, t, idx) + ax = axes(ddu)[1:(end - 1)] + Δt = t[idx + 1] - t[idx] + u₀ = u[ax..., idx] + u₁ = u[ax..., idx + 1] + du₀ = du[ax..., idx] + du₁ = du[ax..., idx + 1] + ddu₀ = ddu[ax..., idx] + ddu₁ = ddu[ax..., idx + 1] + c₁ = (u₁ - u₀ - du₀ * Δt - ddu₀ * Δt^2 / 2) / Δt^3 + c₂ = (3u₀ - 3u₁ + 2(du₀ + du₁ / 2)Δt + ddu₀ * Δt^2 / 2) / Δt^4 + c₃ = (6u₁ - 6u₀ - 3(du₀ + du₁)Δt + (ddu₁ - ddu₀)Δt^2 / 2) / Δt^5 + return c₁, c₂, c₃ +end From 8a6f945bbb0dd111c4b745016d9603b8e9c08364 Mon Sep 17 00:00:00 2001 From: Sathvik Bhagavan Date: Sat, 19 Oct 2024 12:33:09 +0200 Subject: [PATCH 3/5] refactor: add missing derivatives --- src/derivatives.jl | 66 +++++++++++++++++++++++++++++++++++++++------- 1 file changed, 57 insertions(+), 9 deletions(-) diff --git a/src/derivatives.jl b/src/derivatives.jl index 86640360..90c84e38 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -62,9 +62,10 @@ function _derivative(A::LagrangeInterpolation{<:AbstractVector}, t::Number) der end -function _derivative(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number) +function _derivative(A::LagrangeInterpolation{<:AbstractArray}, t::Number) ((t < A.t[1] || t > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError()) - der = zero(A.u[:, 1]) + ax = axes(A.u)[1:(end - 1)] + der = zero(A.u[ax..., 1]) for j in eachindex(A.t) tmp = zero(A.t[1]) if isnan(A.bcache[j]) @@ -91,7 +92,7 @@ function _derivative(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number) tmp += k end end - der += A.u[:, j] * tmp + der += A.u[ax..., j] * tmp end der end @@ -99,15 +100,23 @@ end function _derivative(A::LagrangeInterpolation{<:AbstractVector}, t::Number, idx) _derivative(A, t) end -function _derivative(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number, idx) +function _derivative(A::LagrangeInterpolation{<:AbstractArray}, t::Number, idx) _derivative(A, t) end function _derivative(A::AkimaInterpolation{<:AbstractVector}, t::Number, iguess) idx = get_idx(A, t, iguess; idx_shift = -1, side = :first) - j = min(idx, length(A.c)) # for smooth derivative at A.t[end] + j = min(idx, length(A.p.c)) # for smooth derivative at A.t[end] wj = t - A.t[idx] - @evalpoly wj A.b[idx] 2A.c[j] 3A.d[j] + @evalpoly wj A.p.b[idx] 2A.p.c[j] 3A.p.d[j] +end + +function _derivative(A::AkimaInterpolation{<:AbstractArray}, t::Number, iguess) + idx = get_idx(A, t, iguess; idx_shift = -1, side = :first) + j = min(idx, length(A.p.c)) # for smooth derivative at A.t[end] + wj = t - A.t[idx] + ax = axes(A.u)[1:(end - 1)] + @. @evalpoly wj A.p.b[ax..., idx] 2A.p.c[ax..., j] 3A.p.d[ax..., j] end function _derivative(A::ConstantInterpolation, t::Number, iguess) @@ -119,13 +128,15 @@ function _derivative(A::ConstantInterpolation{<:AbstractVector}, t::Number, igue return isempty(searchsorted(A.t, t)) ? zero(A.u[1]) : eltype(A.u)(NaN) end -function _derivative(A::ConstantInterpolation{<:AbstractMatrix}, t::Number, iguess) +function _derivative(A::ConstantInterpolation{<:AbstractArray}, t::Number, iguess) ((t < A.t[1] || t > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError()) - return isempty(searchsorted(A.t, t)) ? zero(A.u[:, 1]) : eltype(A.u)(NaN) .* A.u[:, 1] + ax = axes(A.u)[1:(end - 1)] + return isempty(searchsorted(A.t, t)) ? zero(A.u[ax..., 1]) : + eltype(A.u)(NaN) .* A.u[ax..., 1] end # QuadraticSpline Interpolation -function _derivative(A::QuadraticSpline{<:AbstractVector}, t::Number, iguess) +function _derivative(A::QuadraticSpline, t::Number, iguess) idx = get_idx(A, t, iguess; lb = 2, ub_shift = 0, side = :first) σ = get_parameters(A, idx - 1) A.z[idx - 1] + 2σ * (t - A.t[idx - 1]) @@ -143,6 +154,18 @@ function _derivative(A::CubicSpline{<:AbstractVector}, t::Number, iguess) dI + dC + dD end +function _derivative(A::CubicSpline{<:AbstractArray}, t::Number, iguess) + idx = get_idx(A, t, iguess) + Δt₁ = t - A.t[idx] + Δt₂ = A.t[idx + 1] - t + ax = axes(A.u)[1:(end - 1)] + dI = (-A.z[ax..., idx] * Δt₂^2 + A.z[ax..., idx + 1] * Δt₁^2) / (2A.h[idx + 1]) + c₁, c₂ = get_parameters(A, idx) + dC = c₁ + dD = -c₂ + dI + dC + dD +end + function _derivative(A::BSplineInterpolation{<:AbstractVector{<:Number}}, t::Number, iguess) # change t into param [0 1] t < A.t[1] && return zero(A.u[1]) @@ -197,6 +220,18 @@ function _derivative( out end +function _derivative( + A::CubicHermiteSpline{<:AbstractArray}, t::Number, iguess) + idx = get_idx(A, t, iguess) + Δt₀ = t - A.t[idx] + Δt₁ = t - A.t[idx + 1] + ax = axes(A.u)[1:(end - 1)] + out = A.du[ax..., idx] + c₁, c₂ = get_parameters(A, idx) + out .+= Δt₀ .* (Δt₀ .* c₂ .+ 2(c₁ .+ Δt₁ .* c₂)) + out +end + # Quintic Hermite Spline function _derivative( A::QuinticHermiteSpline{<:AbstractVector{<:Number}}, t::Number, iguess) @@ -209,3 +244,16 @@ function _derivative( (3c₁ + (3Δt₁ + Δt₀) * c₂ + (3Δt₁^2 + Δt₀ * 2Δt₁) * c₃) out end + +function _derivative( + A::QuinticHermiteSpline{<:AbstractArray}, t::Number, iguess) + idx = get_idx(A, t, iguess) + Δt₀ = t - A.t[idx] + Δt₁ = t - A.t[idx + 1] + ax = axes(A.u)[1:(end - 1)] + out = A.du[ax..., idx] + A.ddu[ax..., idx] * Δt₀ + c₁, c₂, c₃ = get_parameters(A, idx) + out .+= Δt₀^2 .* + (3c₁ .+ (3Δt₁ .+ Δt₀) .* c₂ + (3Δt₁^2 .+ Δt₀ .* 2Δt₁) .* c₃) + out +end From d76d1e392e75a971c7148ed71182881a10f171ca Mon Sep 17 00:00:00 2001 From: Sathvik Bhagavan Date: Sat, 19 Oct 2024 12:33:23 +0200 Subject: [PATCH 4/5] refactor: integral for Akima --- src/integrals.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/integrals.jl b/src/integrals.jl index 3ae893f7..557760ce 100644 --- a/src/integrals.jl +++ b/src/integrals.jl @@ -61,7 +61,6 @@ function _integral(A::QuadraticInterpolation{<:AbstractVector{<:Number}}, t₀ = A.t[idx] t₁ = A.t[idx + 1] t₂ = A.t[idx + 2] - t_sq = (t^2) / 3 l₀, l₁, l₂ = get_parameters(A, idx) Iu₀ = l₀ * t * (t_sq - t * (t₁ + t₂) / 2 + t₁ * t₂) @@ -91,8 +90,8 @@ function _integral(A::AkimaInterpolation{<:AbstractVector{<:Number}}, idx::Number, t::Number) t1 = A.t[idx] - A.u[idx] * (t - t1) + A.b[idx] * ((t - t1)^2 / 2) + A.c[idx] * ((t - t1)^3 / 3) + - A.d[idx] * ((t - t1)^4 / 4) + A.u[idx] * (t - t1) + A.p.b[idx] * ((t - t1)^2 / 2) + A.p.c[idx] * ((t - t1)^3 / 3) + + A.p.d[idx] * ((t - t1)^4 / 4) end _integral(A::LagrangeInterpolation, idx::Number, t::Number) = throw(IntegralNotFoundError()) From d6a0824bbb8cc6afbba9865b181895a7d3c41070 Mon Sep 17 00:00:00 2001 From: Sathvik Bhagavan Date: Sat, 19 Oct 2024 12:58:56 +0200 Subject: [PATCH 5/5] refactor: revert cache for akima as it is always needed --- src/derivatives.jl | 8 ++-- src/integrals.jl | 4 +- src/interpolation_caches.jl | 83 +++++++++++++++++++++++++++++++----- src/interpolation_methods.jl | 4 +- src/parameter_caches.jl | 66 ---------------------------- 5 files changed, 80 insertions(+), 85 deletions(-) diff --git a/src/derivatives.jl b/src/derivatives.jl index 90c84e38..cc57f066 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -106,17 +106,17 @@ end function _derivative(A::AkimaInterpolation{<:AbstractVector}, t::Number, iguess) idx = get_idx(A, t, iguess; idx_shift = -1, side = :first) - j = min(idx, length(A.p.c)) # for smooth derivative at A.t[end] + j = min(idx, length(A.c)) # for smooth derivative at A.t[end] wj = t - A.t[idx] - @evalpoly wj A.p.b[idx] 2A.p.c[j] 3A.p.d[j] + @evalpoly wj A.b[idx] 2A.c[j] 3A.d[j] end function _derivative(A::AkimaInterpolation{<:AbstractArray}, t::Number, iguess) idx = get_idx(A, t, iguess; idx_shift = -1, side = :first) - j = min(idx, length(A.p.c)) # for smooth derivative at A.t[end] + j = min(idx, length(A.c)) # for smooth derivative at A.t[end] wj = t - A.t[idx] ax = axes(A.u)[1:(end - 1)] - @. @evalpoly wj A.p.b[ax..., idx] 2A.p.c[ax..., j] 3A.p.d[ax..., j] + @. @evalpoly wj A.b[ax..., idx] 2A.c[ax..., j] 3A.d[ax..., j] end function _derivative(A::ConstantInterpolation, t::Number, iguess) diff --git a/src/integrals.jl b/src/integrals.jl index 557760ce..8c15d094 100644 --- a/src/integrals.jl +++ b/src/integrals.jl @@ -90,8 +90,8 @@ function _integral(A::AkimaInterpolation{<:AbstractVector{<:Number}}, idx::Number, t::Number) t1 = A.t[idx] - A.u[idx] * (t - t1) + A.p.b[idx] * ((t - t1)^2 / 2) + A.p.c[idx] * ((t - t1)^3 / 3) + - A.p.d[idx] * ((t - t1)^4 / 4) + A.u[idx] * (t - t1) + A.b[idx] * ((t - t1)^2 / 2) + A.c[idx] * ((t - t1)^3 / 3) + + A.d[idx] * ((t - t1)^4 / 4) end _integral(A::LagrangeInterpolation, idx::Number, t::Number) = throw(IntegralNotFoundError()) diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index 226b4f7a..cfca16cf 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -173,24 +173,28 @@ Extrapolation extends the last cubic polynomial on each side. for a test based on the normalized standard deviation of the difference with respect to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2. """ -struct AkimaInterpolation{uType, tType, IType, pType, T, N} <: +struct AkimaInterpolation{uType, tType, IType, bType, cType, dType, T, N} <: AbstractInterpolation{T, N} u::uType t::tType I::IType - p::pType + b::bType + c::cType + d::dType extrapolate::Bool iguesser::Guesser{tType} cache_parameters::Bool linear_lookup::Bool function AkimaInterpolation( - u, t, I, p, extrapolate, cache_parameters, assume_linear_t) + u, t, I, b, c, d, extrapolate, cache_parameters, assume_linear_t) linear_lookup = seems_linear(assume_linear_t, t) N = get_output_dim(u) - new{typeof(u), typeof(t), typeof(I), typeof(p), eltype(u), N}(u, + new{typeof(u), typeof(t), typeof(I), typeof(b), typeof(c), typeof(d), eltype(u), N}(u, t, I, - p, + b, + c, + d, extrapolate, Guesser(t), cache_parameters, @@ -200,14 +204,72 @@ struct AkimaInterpolation{uType, tType, IType, pType, T, N} <: end function AkimaInterpolation( - u, t; extrapolate = false, cache_parameters = false, assume_linear_t = 1e-2) + u::uType, t; extrapolate = false, cache_parameters = false, + assume_linear_t = 1e-2) where {uType <: + AbstractVector{<:Number}} u, t = munge_data(u, t) linear_lookup = seems_linear(assume_linear_t, t) - p = AkimaParameterCache(u, t) + n = length(t) + dt = diff(t) + m = Array{eltype(u)}(undef, n + 3) + m[3:(end - 2)] = diff(u) ./ 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 = 0.5 .* (m[4:end] .+ m[1:(end - 3)]) + dm = abs.(diff(m)) + f1 = dm[3:(n + 2)] + f2 = dm[1:n] + f12 = f1 + f2 + ind = findall(f12 .> 1e-9 * maximum(f12)) + b[ind] = (f1[ind] .* m[ind .+ 1] .+ + f2[ind] .* m[ind .+ 2]) ./ f12[ind] + c = (3.0 .* m[3:(end - 2)] .- 2.0 .* b[1:(end - 1)] .- b[2:end]) ./ dt + d = (b[1:(end - 1)] .+ b[2:end] .- 2.0 .* m[3:(end - 2)]) ./ dt .^ 2 A = AkimaInterpolation( - u, t, nothing, p, extrapolate, cache_parameters, linear_lookup) + u, t, nothing, b, c, d, extrapolate, cache_parameters, linear_lookup) I = cumulative_integral(A, cache_parameters) - AkimaInterpolation(u, t, I, p, extrapolate, cache_parameters, linear_lookup) + AkimaInterpolation(u, t, I, b, c, d, extrapolate, cache_parameters, linear_lookup) +end + +function AkimaInterpolation( + u::uType, t; extrapolate = false, cache_parameters = false, + assume_linear_t = 1e-2) where {uType <: + AbstractArray} + n = length(t) + dt = diff(t) + ax = axes(u)[1:(end - 1)] + su = size(u) + m = zeros(eltype(u), su[1:(end - 1)]..., n + 3) + m[ax..., 3:(end - 2)] .= mapslices( + x -> x ./ dt, diff(u, dims = length(su)); dims = length(su)) + m[ax..., 2] .= 2m[ax..., 3] .- m[ax..., 4] + m[ax..., 1] .= 2m[ax..., 2] .- m[3] + m[ax..., end - 1] .= 2m[ax..., end - 2] - m[ax..., end - 3] + m[ax..., end] .= 2m[ax..., end - 1] .- m[ax..., end - 2] + b = 0.5 .* (m[ax..., 4:end] .+ m[ax..., 1:(end - 3)]) + dm = abs.(diff(m, dims = length(su))) + f1 = dm[ax..., 3:(n + 2)] + f2 = dm[ax..., 1:n] + f12 = f1 .+ f2 + ind = findall(f12 .> 1e-9 * maximum(f12)) + indi = map(i -> i.I, ind) + b[ind] .= (f1[ind] .* + m[CartesianIndex.(map(i -> (i[1:(end - 1)]..., i[end] + 1), indi))] .+ + f2[ind] .* + m[CartesianIndex.(map(i -> (i[1:(end - 1)]..., i[end] + 2), indi))]) ./ + f12[ind] + c = mapslices(x -> x ./ dt, + (3.0 .* m[ax..., 3:(end - 2)] .- 2.0 .* b[ax..., 1:(end - 1)] .- b[ax..., 2:end]); + dims = length(su)) + d = mapslices(x -> x ./ dt .^ 2, + (b[ax..., 1:(end - 1)] .+ b[ax..., 2:end] .- 2.0 .* m[ax..., 3:(end - 2)]); + dims = length(su)) + A = AkimaInterpolation( + u, t, nothing, b, c, d, extrapolate, cache_parameters, linear_lookup) + I = cumulative_integral(A, cache_parameters) + AkimaInterpolation(u, t, I, b, c, d, extrapolate, cache_parameters, linear_lookup) end """ @@ -367,8 +429,7 @@ end function QuadraticSpline( u::uType, t; extrapolate = false, cache_parameters = false, - assume_linear_t = 1e-2) where {uType <: - AbstractArray} + assume_linear_t = 1e-2) where {uType <: AbstractArray} u, t = munge_data(u, t) linear_lookup = seems_linear(assume_linear_t, t) s = length(t) diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index f1a21741..591dd971 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -121,14 +121,14 @@ end function _interpolate(A::AkimaInterpolation{<:AbstractVector}, t::Number, iguess) idx = get_idx(A, t, iguess) wj = t - A.t[idx] - @evalpoly wj A.u[idx] A.p.b[idx] A.p.c[idx] A.p.d[idx] + @evalpoly wj A.u[idx] A.b[idx] A.c[idx] A.d[idx] end function _interpolate(A::AkimaInterpolation{<:AbstractArray}, t::Number, iguess) idx = get_idx(A, t, iguess) wj = t - A.t[idx] ax = axes(A.u)[1:(end - 1)] - @. @evalpoly wj A.u[ax..., idx] A.p.b[ax..., idx] A.p.c[ax..., idx] A.p.d[ax..., idx] + @. @evalpoly wj A.u[ax..., idx] A.b[ax..., idx] A.c[ax..., idx] A.d[ax..., idx] end # ConstantInterpolation Interpolation diff --git a/src/parameter_caches.jl b/src/parameter_caches.jl index be6adbdd..ac09ac36 100644 --- a/src/parameter_caches.jl +++ b/src/parameter_caches.jl @@ -75,72 +75,6 @@ function quadratic_interpolation_parameters(u::AbstractArray{T, N}, t, idx) wher return l₀, l₁, l₂ end -struct AkimaParameterCache{pType} - b::pType - c::pType - d::pType -end - -function AkimaParameterCache(u, t) - b, c, d = akima_interpolation_parameters(u, t) - AkimaParameterCache(b, c, d) -end - -function akima_interpolation_parameters(u::AbstractVector, t) - n = length(t) - dt = diff(t) - m = Array{eltype(u)}(undef, n + 3) - m[3:(end - 2)] = diff(u) ./ 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 = 0.5 .* (m[4:end] .+ m[1:(end - 3)]) - dm = abs.(diff(m)) - f1 = dm[3:(n + 2)] - f2 = dm[1:n] - f12 = f1 + f2 - ind = findall(f12 .> 1e-9 * maximum(f12)) - b[ind] = (f1[ind] .* m[ind .+ 1] .+ - f2[ind] .* m[ind .+ 2]) ./ f12[ind] - c = (3.0 .* m[3:(end - 2)] .- 2.0 .* b[1:(end - 1)] .- b[2:end]) ./ dt - d = (b[1:(end - 1)] .+ b[2:end] .- 2.0 .* m[3:(end - 2)]) ./ dt .^ 2 - return b, c, d -end - -function akima_interpolation_parameters(u::AbstractArray, t) - n = length(t) - dt = diff(t) - ax = axes(u)[1:(end - 1)] - su = size(u) - m = zeros(eltype(u), su[1:(end - 1)]..., n + 3) - m[ax..., 3:(end - 2)] .= mapslices( - x -> x ./ dt, diff(u, dims = length(su)); dims = length(su)) - m[ax..., 2] .= 2m[ax..., 3] .- m[ax..., 4] - m[ax..., 1] .= 2m[ax..., 2] .- m[3] - m[ax..., end - 1] .= 2m[ax..., end - 2] - m[ax..., end - 3] - m[ax..., end] .= 2m[ax..., end - 1] .- m[ax..., end - 2] - b = 0.5 .* (m[ax..., 4:end] .+ m[ax..., 1:(end - 3)]) - dm = abs.(diff(m, dims = length(su))) - f1 = dm[ax..., 3:(n + 2)] - f2 = dm[ax..., 1:n] - f12 = f1 .+ f2 - ind = findall(f12 .> 1e-9 * maximum(f12)) - indi = map(i -> i.I, ind) - b[ind] .= (f1[ind] .* - m[CartesianIndex.(map(i -> (i[1:(end - 1)]..., i[end] + 1), indi))] .+ - f2[ind] .* - m[CartesianIndex.(map(i -> (i[1:(end - 1)]..., i[end] + 2), indi))]) ./ - f12[ind] - c = mapslices(x -> x ./ dt, - (3.0 .* m[ax..., 3:(end - 2)] .- 2.0 .* b[ax..., 1:(end - 1)] .- b[ax..., 2:end]); - dims = length(su)) - d = mapslices(x -> x ./ dt .^ 2, - (b[ax..., 1:(end - 1)] .+ b[ax..., 2:end] .- 2.0 .* m[ax..., 3:(end - 2)]); - dims = length(su)) - return b, c, d -end - struct QuadraticSplineParameterCache{pType} σ::pType end