Cache sigma for QuadraticSplineInterpolation

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
 

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)

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
 
 """

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

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