From 07b968c67d9147133115860d6945456b1cf4280e Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Sat, 29 Jun 2024 13:33:14 +0200 Subject: [PATCH] Cache sigma for QuadraticSplineInterpolation --- src/derivatives.jl | 2 +- src/integrals.jl | 11 +++-------- src/interpolation_caches.jl | 22 +++++++++++++--------- src/interpolation_methods.jl | 8 ++++---- src/parameter_caches.jl | 18 ++++++++++++++++-- 5 files changed, 37 insertions(+), 24 deletions(-) diff --git a/src/derivatives.jl b/src/derivatives.jl index 14635d5b..ac3ba999 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -125,7 +125,7 @@ end # QuadraticSpline Interpolation function _derivative(A::QuadraticSpline{<:AbstractVector}, t::Number, iguess) idx = get_idx(A.t, t, iguess; lb = 2, ub_shift = 0, side = :first) - σ = 1 // 2 * (A.z[idx] - A.z[idx - 1]) / (A.t[idx] - A.t[idx - 1]) + σ = A.p.σ[idx-1] A.z[idx - 1] + 2σ * (t - A.t[idx - 1]), idx end diff --git a/src/integrals.jl b/src/integrals.jl index 92c1faeb..26af6909 100644 --- a/src/integrals.jl +++ b/src/integrals.jl @@ -56,14 +56,9 @@ function _integral(A::QuadraticInterpolation{<:AbstractVector{<:Number}}, end function _integral(A::QuadraticSpline{<:AbstractVector{<:Number}}, idx::Number, t::Number) - t1 = A.t[idx] - t2 = A.t[idx + 1] - u1 = A.u[idx] - z1 = A.z[idx] - z2 = A.z[idx + 1] - t^3 * (z1 - z2) / (6 * t1 - 6 * t2) + t^2 * (t1 * z2 - t2 * z1) / (2 * t1 - 2 * t2) + - t * (-t1^2 * z1 - t1^2 * z2 + 2 * t1 * t2 * z1 + 2 * t1 * u1 - 2 * t2 * u1) / - (2 * t1 - 2 * t2) + Cᵢ = A.u[idx] + Δt = t - A.t[idx] + return A.z[idx] * Δt^2 / 2 + A.p.σ[idx] * Δt^3 / 3 + Cᵢ * Δt end function _integral(A::CubicSpline{<:AbstractVector{<:Number}}, idx::Number, t::Number) diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index 405714db..fdfc2ad3 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -17,11 +17,11 @@ Extrapolation extends the last linear polynomial on each side. struct LinearInterpolation{uType, tType, pType, T} <: AbstractInterpolation{T} u::uType t::tType - p::pType + p::LinearParameterCache{pType} extrapolate::Bool idx_prev::Base.RefValue{Int} function LinearInterpolation(u, t, p, extrapolate) - new{typeof(u), typeof(t), typeof(p), eltype(u)}(u, t, p, extrapolate, Ref(1)) + new{typeof(u), typeof(t), typeof(p.slope), eltype(u)}(u, t, p, extrapolate, Ref(1)) end end @@ -51,14 +51,14 @@ Extrapolation extends the last quadratic polynomial on each side. struct QuadraticInterpolation{uType, tType, pType, T} <: AbstractInterpolation{T} u::uType t::tType - p::pType + p::QuadraticParameterCache{pType} mode::Symbol extrapolate::Bool idx_prev::Base.RefValue{Int} function QuadraticInterpolation(u, t, p, mode, extrapolate) mode ∈ (:Forward, :Backward) || error("mode should be :Forward or :Backward for QuadraticInterpolation") - new{typeof(u), typeof(t), typeof(p), eltype(u)}(u, t, p, mode, extrapolate, Ref(1)) + new{typeof(u), typeof(t), typeof(p.l₀), eltype(u)}(u, t, p, mode, extrapolate, Ref(1)) end end @@ -235,19 +235,21 @@ Extrapolation extends the last quadratic polynomial on each side. - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`. - `safetycopy`: boolean value to make a copy of `u` and `t`. Defaults to `true`. """ -struct QuadraticSpline{uType, tType, tAType, dType, zType, T} <: +struct QuadraticSpline{uType, tType, pType, tAType, dType, zType, T} <: AbstractInterpolation{T} u::uType t::tType + p::QuadraticSplineParameterCache{pType} tA::tAType d::dType z::zType extrapolate::Bool idx_prev::Base.RefValue{Int} - function QuadraticSpline(u, t, tA, d, z, extrapolate) - new{typeof(u), typeof(t), typeof(tA), + function QuadraticSpline(u, t, p, tA, d, z, extrapolate) + new{typeof(u), typeof(t), typeof(p.σ), typeof(tA), typeof(d), typeof(z), eltype(u)}(u, t, + p, tA, d, z, @@ -272,7 +274,8 @@ function QuadraticSpline(u::uType, d = map(i -> i == 1 ? typed_zero : 2 // 1 * (u[i] - u[i - 1]) / (t[i] - t[i - 1]), 1:s) z = tA \ d - QuadraticSpline(u, t, tA, d, z, extrapolate) + p = QuadraticSplineParameterCache(z, t) + QuadraticSpline(u, t, p, tA, d, z, extrapolate) end function QuadraticSpline( @@ -290,7 +293,8 @@ function QuadraticSpline( d = transpose(reshape(reduce(hcat, d_), :, s)) z_ = reshape(transpose(tA \ d), size(u[1])..., :) z = [z_s for z_s in eachslice(z_, dims = ndims(z_))] - QuadraticSpline(u, t, tA, d, z, extrapolate) + p = QuadraticSplineParameterCache(z, t) + QuadraticSpline(u, t, p, tA, d, z, extrapolate) end """ diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index c523dd6f..43c7b08c 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -154,10 +154,10 @@ end # QuadraticSpline Interpolation function _interpolate(A::QuadraticSpline{<:AbstractVector}, t::Number, iguess) - idx = get_idx(A.t, t, iguess; lb = 2, ub_shift = 0, side = :first) - Cᵢ = A.u[idx - 1] - σ = 1 // 2 * (A.z[idx] - A.z[idx - 1]) / (A.t[idx] - A.t[idx - 1]) - return A.z[idx - 1] * (t - A.t[idx - 1]) + σ * (t - A.t[idx - 1])^2 + Cᵢ, idx + idx = get_idx(A.t, t, iguess) + Cᵢ = A.u[idx] + Δt = t - A.t[idx] + return A.z[idx] * Δt + A.p.σ[idx] * Δt^2 + Cᵢ, idx end # CubicSpline Interpolation diff --git a/src/parameter_caches.jl b/src/parameter_caches.jl index 14690625..40c7f7f3 100644 --- a/src/parameter_caches.jl +++ b/src/parameter_caches.jl @@ -3,7 +3,7 @@ struct LinearParameterCache{pType} end function LinearParameterCache(u, t) - slope = LinearInterpolationParameters.(Ref(u), Ref(t), Base.OneTo(length(t) - 1)) + slope = LinearInterpolationParameters.(Ref(u), Ref(t), 1:(length(t) - 1)) return LinearParameterCache(slope) end @@ -23,7 +23,7 @@ end function QuadraticParameterCache(u, t) parameters = QuadraticInterpolationParameters.( - Ref(u), Ref(t), Base.OneTo(length(t) - 2)) + Ref(u), Ref(t), 1:(length(t) - 2)) l₀, l₁, l₂ = collect.(eachrow(hcat(collect.(parameters)...))) return QuadraticParameterCache(l₀, l₁, l₂) end @@ -49,3 +49,17 @@ function QuadraticInterpolationParameters(u, t, idx) l₂ = u₂ / (Δt₂ * Δt₁) return l₀, l₁, l₂ end + +struct QuadraticSplineParameterCache{pType} + σ::pType +end + +function QuadraticSplineParameterCache(z, t) + σ = QuadraticSplineInterpolationParameters.(Ref(z), Ref(t), 1:(length(t) - 1)) + return QuadraticSplineParameterCache(σ) +end + +function QuadraticSplineInterpolationParameters(z, t, idx) + σ = 1 // 2 * (z[idx+1] - z[idx]) / (t[idx+1] - t[idx]) + return σ +end \ No newline at end of file