Cache parameters

Project.toml

@@ -7,6 +7,7 @@ FindFirstFunctions = "64ca27bc-2ba2-4a57-88aa-44e436879224"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
+ReadOnlyArrays = "988b38a3-91fc-5605-94a2-ee2116b3bd83"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 
@@ -32,6 +33,7 @@ LinearAlgebra = "1.10"
 Optim = "1.6"
 PrettyTables = "2"
 QuadGK = "2.9.1"
+ReadOnlyArrays = "0.2.0"
 RecipesBase = "1.3"
 Reexport = "1"
 RegularizationTools = "0.6"

docs/src/interface.md

@@ -17,7 +17,7 @@ t = [0.0, 62.25, 109.66, 162.66, 205.8, 252.3]
 
 All interpolation methods return an object from which we can compute the value of the dependent variable at any time point.
 
-We will use the `CubicSpline` method for demonstration, but the API is the same for all the methods. We can also pass the `extrapolate=true` keyword if we want to allow the interpolation to go beyond the range of the timepoints. The default value is `extrapolate=false`.
+We will use the `CubicSpline` method for demonstration, but the API is the same for all the methods. We can also pass the `extrapolate = true` keyword if we want to allow the interpolation to go beyond the range of the timepoints. The default value is `extrapolate = false`.
 
 ```@example interface
 A1 = CubicSpline(u, t)
@@ -35,6 +35,23 @@ A2(300.0)
 
     The values computed beyond the range of the time points provided during interpolation will not be reliable, as these methods only perform well within the range and the first/last piece polynomial fit is extrapolated on either side which might not reflect the true nature of the data.
 
+The keyword `safetycopy = false` can be passed to make sure no copies of `u` and `t` are made when initializing the interpolation object.
+
+```@example interface
+A3 = QuadraticInterpolation(u, t; safetycopy = false)
+
+# Check for same memory
+u === A3.u.parent
+```
+
+Note that this does not prevent allocation in every interpolation constructor call, because parameter values are cached for all interpolation types except [`ConstantInterpolation`](@ref).
+
+Because of the caching of parameters which depend on `u` and `t`, this data should not be mutated. Therefore `u` and `t` are wrapped in a `ReadOnlyArray` from [ReadOnlyArrays.jl](https://github.com/JuliaArrays/ReadOnlyArrays.jl).
+
+```@repl interface
+A3.t[2] = 3.14
+```
+
 ## Derivatives
 
 Derivatives of the interpolated curves can also be computed at any point for all the methods. Derivatives upto second order is supported where first order derivative is computed analytically and second order using `ForwardDiff.jl`. Order is passed as the third argument. It is 1 by default.

ext/DataInterpolationsOptimExt.jl

@@ -18,8 +18,9 @@ function Curvefit(u,
         box = false,
         lb = nothing,
         ub = nothing;
-        extrapolate = false)
-    u, t = munge_data(u, t)
+        extrapolate = false,
+        safetycopy = false)
+    u, t = munge_data(u, t, safetycopy)
     errfun(t, u, p) = sum(abs2.(u .- model(t, p)))
     if box == false
         mfit = optimize(p -> errfun(t, u, p), p0, alg)

ext/DataInterpolationsRegularizationToolsExt.jl

@@ -69,8 +69,8 @@ A = RegularizationSmooth(u, t, t̂, wls, wr, d; λ = 1.0, alg = :gcv_svd)
 """
 function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::AbstractVector,
         wls::AbstractVector, wr::AbstractVector, d::Int = 2;
-        λ::Real = 1.0, alg::Symbol = :gcv_svd, extrapolate::Bool = false)
-    u, t = munge_data(u, t)
+        λ::Real = 1.0, alg::Symbol = :gcv_svd, extrapolate::Bool = false, safetycopy::Bool = true)
+    u, t = munge_data(u, t, safetycopy)
     M = _mapping_matrix(t̂, t)
     Wls½ = LA.diagm(sqrt.(wls))
     Wr½ = LA.diagm(sqrt.(wr))
@@ -86,8 +86,8 @@ A = RegularizationSmooth(u, t, d; λ = 1.0, alg = :gcv_svd, extrapolate = false)
 """
 function RegularizationSmooth(u::AbstractVector, t::AbstractVector, d::Int = 2;
         λ::Real = 1.0,
-        alg::Symbol = :gcv_svd, extrapolate::Bool = false)
-    u, t = munge_data(u, t)
+        alg::Symbol = :gcv_svd, extrapolate::Bool = false, safetycopy::Bool = true)
+    u, t = munge_data(u, t, safetycopy)
     t̂ = t
     N = length(t)
     M = Array{Float64}(LA.I, N, N)
@@ -114,8 +114,9 @@ A = RegularizationSmooth(u, t, t̂, d; λ = 1.0, alg = :gcv_svd, extrapolate = f
 ```
 """
 function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::AbstractVector,
-        d::Int = 2; λ::Real = 1.0, alg::Symbol = :gcv_svd, extrapolate::Bool = false)
-    u, t = munge_data(u, t)
+        d::Int = 2; λ::Real = 1.0, alg::Symbol = :gcv_svd,
+        extrapolate::Bool = false, safetycopy::Bool = true)
+    u, t = munge_data(u, t, safetycopy)
     N, N̂ = length(t), length(t̂)
     M = _mapping_matrix(t̂, t)
     Wls½ = Array{Float64}(LA.I, N, N)
@@ -142,8 +143,8 @@ A = RegularizationSmooth(u, t, t̂, wls, d; λ = 1.0, alg = :gcv_svd, extrapolat
 """
 function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::AbstractVector,
         wls::AbstractVector, d::Int = 2; λ::Real = 1.0,
-        alg::Symbol = :gcv_svd, extrapolate::Bool = false)
-    u, t = munge_data(u, t)
+        alg::Symbol = :gcv_svd, extrapolate::Bool = false, safetycopy::Bool = true)
+    u, t = munge_data(u, t, safetycopy)
     N, N̂ = length(t), length(t̂)
     M = _mapping_matrix(t̂, t)
     Wls½ = LA.diagm(sqrt.(wls))
@@ -171,8 +172,8 @@ A = RegularizationSmooth(
 """
 function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing,
         wls::AbstractVector, d::Int = 2; λ::Real = 1.0,
-        alg::Symbol = :gcv_svd, extrapolate::Bool = false)
-    u, t = munge_data(u, t)
+        alg::Symbol = :gcv_svd, extrapolate::Bool = false, safetycopy::Bool = true)
+    u, t = munge_data(u, t, safetycopy)
     t̂ = t
     N = length(t)
     M = Array{Float64}(LA.I, N, N)
@@ -201,8 +202,8 @@ A = RegularizationSmooth(
 """
 function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing,
         wls::AbstractVector, wr::AbstractVector, d::Int = 2;
-        λ::Real = 1.0, alg::Symbol = :gcv_svd, extrapolate::Bool = false)
-    u, t = munge_data(u, t)
+        λ::Real = 1.0, alg::Symbol = :gcv_svd, extrapolate::Bool = false, safetycopy::Bool = true)
+    u, t = munge_data(u, t, safetycopy)
     t̂ = t
     N = length(t)
     M = Array{Float64}(LA.I, N, N)
@@ -231,8 +232,8 @@ A = RegularizationSmooth(
 """
 function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing,
         wls::Symbol, d::Int = 2; λ::Real = 1.0, alg::Symbol = :gcv_svd,
-        extrapolate::Bool = false)
-    u, t = munge_data(u, t)
+        extrapolate::Bool = false, safetycopy::Bool = true)
+    u, t = munge_data(u, t, safetycopy)
     t̂ = t
     N = length(t)
     M = Array{Float64}(LA.I, N, N)

src/DataInterpolations.jl

@@ -7,6 +7,7 @@ abstract type AbstractInterpolation{T} end
 using LinearAlgebra, RecipesBase
 using PrettyTables
 using ForwardDiff
+using ReadOnlyArrays
 import FindFirstFunctions: searchsortedfirstcorrelated, searchsortedlastcorrelated,
                            bracketstrictlymontonic
 
@@ -74,6 +75,12 @@ function Base.showerror(io::IO, e::DerivativeNotFoundError)
     print(io, DERIVATIVE_NOT_FOUND_ERROR)
 end
 
+const MUST_COPY_ERROR = "A copy must be made of u, t to filter missing data"
+struct MustCopyError <: Exception end
+function Base.showerror(io::IO, e::MustCopyError)
+    print(io, MUST_COPY_ERROR)
+end
+
 export LinearInterpolation, QuadraticInterpolation, LagrangeInterpolation,
        AkimaInterpolation, ConstantInterpolation, QuadraticSpline, CubicSpline,
        BSplineInterpolation, BSplineApprox, CubicHermiteSpline,
@@ -105,12 +112,13 @@ struct RegularizationSmooth{uType, tType, T, T2, ITP <: AbstractInterpolation{T}
             alg,
             Aitp,
             extrapolate)
-        new{typeof(u), typeof(t), eltype(u), typeof(λ), typeof(Aitp)}(u,
-            û,
-            t,
-            t̂,
-            wls,
-            wr,
+        new{typeof(u), typeof(t), eltype(u), typeof(λ), typeof(Aitp)}(
+            readonly_wrap(u),
+            readonly_wrap(û),
+            readonly_wrap(t),
+            readonly_wrap(t̂),
+            readonly_wrap(oftype(u.parent, wls)),
+            readonly_wrap(oftype(u.parent, wr)),
             d,
             λ,
             alg,

src/derivatives.jl

@@ -17,30 +17,17 @@ function derivative(A, t, order = 1)
     end
 end
 
-function _derivative(A::LinearInterpolation{<:AbstractVector}, t::Number, iguess)
+function _derivative(A::LinearInterpolation, t::Number, iguess)
     idx = get_idx(A.t, t, iguess; idx_shift = -1, ub_shift = -2, side = :first)
-    (A.u[idx + 1] - A.u[idx]) / (A.t[idx + 1] - A.t[idx]), idx
+    A.p.slope[idx], idx
 end
 
-function _derivative(A::LinearInterpolation{<:AbstractMatrix}, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess; idx_shift = -1, ub_shift = -2, side = :first)
-    (@views @. (A.u[:, idx + 1] - A.u[:, idx]) / (A.t[idx + 1] - A.t[idx])), idx
-end
-
-function _derivative(A::QuadraticInterpolation{<:AbstractVector}, t::Number, iguess)
-    i₀, i₁, i₂ = _quad_interp_indices(A, t, iguess)
-    dl₀ = (2t - A.t[i₁] - A.t[i₂]) / ((A.t[i₀] - A.t[i₁]) * (A.t[i₀] - A.t[i₂]))
-    dl₁ = (2t - A.t[i₀] - A.t[i₂]) / ((A.t[i₁] - A.t[i₀]) * (A.t[i₁] - A.t[i₂]))
-    dl₂ = (2t - A.t[i₀] - A.t[i₁]) / ((A.t[i₂] - A.t[i₀]) * (A.t[i₂] - A.t[i₁]))
-    A.u[i₀] * dl₀ + A.u[i₁] * dl₁ + A.u[i₂] * dl₂, i₀
-end
-
-function _derivative(A::QuadraticInterpolation{<:AbstractMatrix}, t::Number, iguess)
+function _derivative(A::QuadraticInterpolation, t::Number, iguess)
     i₀, i₁, i₂ = _quad_interp_indices(A, t, iguess)
-    dl₀ = (2t - A.t[i₁] - A.t[i₂]) / ((A.t[i₀] - A.t[i₁]) * (A.t[i₀] - A.t[i₂]))
-    dl₁ = (2t - A.t[i₀] - A.t[i₂]) / ((A.t[i₁] - A.t[i₀]) * (A.t[i₁] - A.t[i₂]))
-    dl₂ = (2t - A.t[i₀] - A.t[i₁]) / ((A.t[i₂] - A.t[i₀]) * (A.t[i₂] - A.t[i₁]))
-    (@views @. A.u[:, i₀] * dl₀ + A.u[:, i₁] * dl₁ + A.u[:, i₂] * dl₂), i₀
+    du₀ = A.p.l₀[i₀] * (2t - A.t[i₁] - A.t[i₂])
+    du₁ = A.p.l₁[i₀] * (2t - A.t[i₀] - A.t[i₂])
+    du₂ = A.p.l₂[i₀] * (2t - A.t[i₀] - A.t[i₁])
+    return @views @. du₀ + du₁ + du₂, i₀
 end
 
 function _derivative(A::LagrangeInterpolation{<:AbstractVector}, t::Number)
@@ -125,6 +112,10 @@ function _derivative(A::AkimaInterpolation{<:AbstractVector}, t::Number, iguess)
     (@evalpoly wj A.b[idx] 2A.c[j] 3A.d[j]), idx
 end
 
+function _derivative(A::ConstantInterpolation, t::Number, iguess)
+    return zero(first(A.u)), iguess
+end
+
 function _derivative(A::ConstantInterpolation{<:AbstractVector}, t::Number)
     ((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)
@@ -138,17 +129,18 @@ 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
 
 # CubicSpline Interpolation
 function _derivative(A::CubicSpline{<:AbstractVector}, t::Number, iguess)
     idx = get_idx(A.t, t, iguess)
-    dI = -3A.z[idx] * (A.t[idx + 1] - t)^2 / (6A.h[idx + 1]) +
-         3A.z[idx + 1] * (t - A.t[idx])^2 / (6A.h[idx + 1])
-    dC = A.u[idx + 1] / A.h[idx + 1] - A.z[idx + 1] * A.h[idx + 1] / 6
-    dD = -(A.u[idx] / A.h[idx + 1] - A.z[idx] * A.h[idx + 1] / 6)
+    Δt₁ = t - A.t[idx]
+    Δt₂ = A.t[idx + 1] - t
+    dI = (-A.z[idx] * Δt₂^2 + A.z[idx + 1] * Δt₁^2) / (2A.h[idx + 1])
+    dC = A.p.c₁[idx]
+    dD = -A.p.c₂[idx]
     dI + dC + dD, idx
 end
 
@@ -160,7 +152,8 @@ function _derivative(A::BSplineInterpolation{<:AbstractVector{<:Number}}, t::Num
     n = length(A.t)
     scale = (A.p[idx + 1] - A.p[idx]) / (A.t[idx + 1] - A.t[idx])
     t_ = A.p[idx] + (t - A.t[idx]) * scale
-    N = DataInterpolations.spline_coefficients(n, A.d - 1, A.k, t_)
+    N = t isa ForwardDiff.Dual ? zeros(eltype(t), n) : A.N
+    spline_coefficients!(N, A.d - 1, A.k, t_)
     ducum = zero(eltype(A.u))
     if t == A.t[1]
         ducum = (A.c[2] - A.c[1]) / (A.k[A.d + 2])
@@ -180,7 +173,8 @@ function _derivative(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, ig
     idx = get_idx(A.t, t, iguess)
     scale = (A.p[idx + 1] - A.p[idx]) / (A.t[idx + 1] - A.t[idx])
     t_ = A.p[idx] + (t - A.t[idx]) * scale
-    N = spline_coefficients(A.h, A.d - 1, A.k, t_)
+    N = t isa ForwardDiff.Dual ? zeros(eltype(t), A.h) : A.N
+    spline_coefficients!(N, A.d - 1, A.k, t_)
     ducum = zero(eltype(A.u))
     if t == A.t[1]
         ducum = (A.c[2] - A.c[1]) / (A.k[A.d + 2])

src/integrals.jl

@@ -25,11 +25,8 @@ end
 function _integral(A::LinearInterpolation{<:AbstractVector{<:Number}},
         idx::Number,
         t::Number)
-    t1 = A.t[idx]
-    t2 = A.t[idx + 1]
-    u1 = A.u[idx]
-    u2 = A.u[idx + 1]
-    t^2 * (u1 - u2) / (2 * t1 - 2 * t2) + t * (t1 * u2 - t2 * u1) / (t1 - t2)
+    Δt = t - A.t[idx]
+    Δt * (A.u[idx] + A.p.slope[idx] * Δt / 2)
 end
 
 function _integral(A::ConstantInterpolation{<:AbstractVector}, idx::Number, t::Number)
@@ -47,49 +44,30 @@ function _integral(A::QuadraticInterpolation{<:AbstractVector{<:Number}},
         t::Number)
     A.mode == :Backward && idx > 1 && (idx -= 1)
     idx = min(length(A.t) - 2, idx)
-    t1 = A.t[idx]
-    t2 = A.t[idx + 1]
-    t3 = A.t[idx + 2]
-    u1 = A.u[idx]
-    u2 = A.u[idx + 1]
-    u3 = A.u[idx + 2]
-    (t^3 * (-t1 * u2 + t1 * u3 + t2 * u1 - t2 * u3 - t3 * u1 + t3 * u2) /
-     (3 * t1^2 * t2 - 3 * t1^2 * t3 - 3 * t1 * t2^2 + 3 * t1 * t3^2 + 3 * t2^2 * t3 -
-      3 * t2 * t3^2) +
-     t^2 * (t1^2 * u2 - t1^2 * u3 - t2^2 * u1 + t2^2 * u3 + t3^2 * u1 - t3^2 * u2) /
-     (2 * t1^2 * t2 - 2 * t1^2 * t3 - 2 * t1 * t2^2 + 2 * t1 * t3^2 + 2 * t2^2 * t3 -
-      2 * t2 * t3^2) +
-     t *
-     (t1^2 * t2 * u3 - t1^2 * t3 * u2 - t1 * t2^2 * u3 + t1 * t3^2 * u2 + t2^2 * t3 * u1 -
-      t2 * t3^2 * u1) /
-     (t1^2 * t2 - t1^2 * t3 - t1 * t2^2 + t1 * t3^2 + t2^2 * t3 - t2 * t3^2))
+    t₀ = A.t[idx]
+    t₁ = A.t[idx + 1]
+    t₂ = A.t[idx + 2]
+
+    t_sq = (t^2) / 3
+    Iu₀ = A.p.l₀[idx] * t * (t_sq - t * (t₁ + t₂) / 2 + t₁ * t₂)
+    Iu₁ = A.p.l₁[idx] * t * (t_sq - t * (t₀ + t₂) / 2 + t₀ * t₂)
+    Iu₂ = A.p.l₂[idx] * t * (t_sq - t * (t₀ + t₁) / 2 + t₀ * t₁)
+    return Iu₀ + Iu₁ + Iu₂
 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)
-    t1 = A.t[idx]
-    t2 = A.t[idx + 1]
-    u1 = A.u[idx]
-    u2 = A.u[idx + 1]
-    z1 = A.z[idx]
-    z2 = A.z[idx + 1]
-    h2 = A.h[idx + 1]
-    (t^4 * (-z1 + z2) / (24 * h2) + t^3 * (-t1 * z2 + t2 * z1) / (6 * h2) +
-     t^2 * (h2^2 * z1 - h2^2 * z2 + 3 * t1^2 * z2 - 3 * t2^2 * z1 - 6 * u1 + 6 * u2) /
-     (12 * h2) +
-     t *
-     (h2^2 * t1 * z2 - h2^2 * t2 * z1 - t1^3 * z2 - 6 * t1 * u2 + t2^3 * z1 + 6 * t2 * u1) /
-     (6 * h2))
+    Δt₁sq = (t - A.t[idx])^2 / 2
+    Δt₂sq = (A.t[idx + 1] - t)^2 / 2
+    II = (-A.z[idx] * Δt₂sq^2 + A.z[idx + 1] * Δt₁sq^2) / (6A.h[idx + 1])
+    IC = A.p.c₁[idx] * Δt₁sq
+    ID = -A.p.c₂[idx] * Δt₂sq
+    II + IC + ID
 end
 
 function _integral(A::AkimaInterpolation{<:AbstractVector{<:Number}},
@@ -100,24 +78,9 @@ function _integral(A::AkimaInterpolation{<:AbstractVector{<:Number}},
     A.d[idx] * ((t - t1)^4 / 4)
 end
 
-integral(A::LagrangeInterpolation, t1::Number, t2::Number) = throw(IntegralNotFoundError())
-integral(A::LagrangeInterpolation, t::Number) = throw(IntegralNotFoundError())
-
-function integral(A::BSplineInterpolation{<:AbstractVector{<:Number}},
-        t1::Number,
-        t2::Number)
-    throw(IntegralNotFoundError())
-end
-function integral(A::BSplineInterpolation{<:AbstractVector{<:Number}}, t::Number)
-    throw(IntegralNotFoundError())
-end
-
-function integral(A::BSplineApprox{<:AbstractVector{<:Number}}, t1::Number, t2::Number)
-    throw(IntegralNotFoundError())
-end
-function integral(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number)
-    throw(IntegralNotFoundError())
-end
+_integral(A::LagrangeInterpolation, idx::Number, t::Number) = throw(IntegralNotFoundError())
+_integral(A::BSplineInterpolation, idx::Number, t::Number) = throw(IntegralNotFoundError())
+_integral(A::BSplineApprox, idx::Number, t::Number) = throw(IntegralNotFoundError())
 
 # Cubic Hermite Spline
 function _integral(