refactor: revert cache for akima as it is always needed

Author: sathvikbhagavan

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)

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())

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)

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

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