Add smoothed linear interpolation

docs/src/index.md

@@ -20,6 +20,7 @@ corresponding to `(u,t)` pairs.
 
   - `SmoothedConstantInterpolation(u,t)` - An integral preserving continuously differentiable approximation of constant interpolation.
   - `LinearInterpolation(u,t)` - A linear interpolation.
+  - `SmoothedLinearInterpolation(u,t;  λ = 0.25)` - A continuously differentiable linear interpolation with an interval around each data point replaced by spline section.
   - `QuadraticInterpolation(u,t)` - A quadratic interpolation.
   - `LagrangeInterpolation(u,t,n)` - A Lagrange interpolation of order `n`.
   - `QuadraticSpline(u,t)` - A quadratic spline interpolation.

docs/src/manual.md

@@ -2,6 +2,7 @@
 
 ```@docs
 LinearInterpolation
+SmoothedLinearInterpolation
 QuadraticInterpolation
 LagrangeInterpolation
 AkimaInterpolation

docs/src/methods.md

@@ -25,6 +25,17 @@ A = LinearInterpolation(u, t)
 plot(A)
 ```
 
+## Smoothed Linear Interpolation
+
+This is a continuously differentiable linear interpolation with an interval around each data point replaced by spline section.
+
+```@example tutorial
+A = LinearInterpolation(u, t)
+plot(A)
+A = SmoothedLinearInterpolation(u, t)
+plot!(A)
+```
+
 ## Quadratic Interpolation
 
 This function fits a parabola passing through the two nearest points from the input

src/DataInterpolations.jl

@@ -118,10 +118,10 @@ _output_size(::AbstractVector{<:Number}) = ()
 _output_size(u::AbstractVector) = size(first(u))
 _output_size(u::AbstractArray) = Base.front(size(u))
 
-export LinearInterpolation, QuadraticInterpolation, LagrangeInterpolation,
-       AkimaInterpolation, ConstantInterpolation, SmoothedConstantInterpolation,
-       QuadraticSpline, CubicSpline, BSplineInterpolation, BSplineApprox,
-       CubicHermiteSpline, PCHIPInterpolation, QuinticHermiteSpline,
+export LinearInterpolation, SmoothedLinearInterpolation, QuadraticInterpolation,
+       LagrangeInterpolation, AkimaInterpolation, ConstantInterpolation,
+       SmoothedConstantInterpolation, QuadraticSpline, CubicSpline, BSplineInterpolation,
+       BSplineApprox, CubicHermiteSpline, PCHIPInterpolation, QuinticHermiteSpline,
        SmoothArcLengthInterpolation, LinearInterpolationIntInv,
        ConstantInterpolationIntInv, ExtrapolationType
 export output_dim, output_size

src/derivatives.jl

@@ -100,6 +100,26 @@ function _derivative(A::LinearInterpolation, t::Number, iguess)
     slope
 end
 
+function _derivative(A::SmoothedLinearInterpolation, t::Number, iguess)
+    idx = get_idx(A, t, iguess)
+
+    # Check against boundary points between linear and spline interpolation
+    left = (t < A.p.t_tilde[2idx])
+    right = (t > A.p.t_tilde[2idx + 1])
+
+    # Spline interpolation
+    if left && idx != 1
+        return U_deriv(A, t, idx)
+    end
+
+    if right && idx != length(A.u) - 1
+        return U_deriv(A, t, idx + 1)
+    end
+
+    # Linear interpolation
+    return A.p.slope[idx + 1]
+end
+
 function _derivative(A::SmoothedConstantInterpolation{<:AbstractVector}, t::Number, iguess)
     idx = get_idx(A, t, iguess)
     d_lower, d_upper, c_lower, c_upper = get_parameters(A, idx)

src/integrals.jl

@@ -296,3 +296,7 @@ function _integral(
     integrate_quintic_polynomial(t1, t2, tᵢ, A.u[idx], A.du[idx], A.ddu[idx] / 2,
         c₁ + Δt * (-c₂ + c₃ * Δt), c₂ - 2c₃ * Δt, c₃)
 end
+
+function _integral(A::SmoothedLinearInterpolation, idx::Number, t1::Number, t2::Number)
+    throw(IntegralNotFoundError())
+end

src/interpolation_caches.jl

@@ -64,6 +64,76 @@ function LinearInterpolation(
         cache_parameters, assume_linear_t)
 end
 
+"""
+    SmoothedLinearInterpolation(
+    u, t; extrapolation::ExtrapolationType.T = ExtrapolationType.None,
+    extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
+    extrapolation_right::ExtrapolationType.T = ExtrapolationType.None,
+    cache_parameters = true, assume_linear_t = 1e-2, λ = 0.25)
+
+    A linear interpolation method where a section around each data point is replaced by spline
+    section to make the transition at the datapoint differentiable. The spline section is guaranteed
+    to be in the convex hull of the boundary points of the spline and the original data point.
+
+    ## Arguments
+
+  - `u`: data points.
+  - `t`: time points.
+
+## Keyword Arguments
+
+  - `extrapolation`: The extrapolation type applied left and right of the data. Possible options
+    are `ExtrapolationType.None` (default), `ExtrapolationType.Constant`, `ExtrapolationType.Linear`
+    `ExtrapolationType.Extension`, `ExtrapolationType.Periodic` and `ExtrapolationType.Reflective`.
+  - `extrapolation_left`: The extrapolation type applied left of the data. See `extrapolation` for
+    the possible options. This keyword is ignored if `extrapolation != Extrapolation.none`.
+  - `extrapolation_right`: The extrapolation type applied right of the data. See `extrapolation` for
+    the possible options. This keyword is ignored if `extrapolation != Extrapolation.none`.
+  - `cache_parameters`: precompute parameters at initialization for faster interpolation
+    computations. Note: if activated, `u` and `t` should not be modified. Defaults to `false`.
+  - `assume_linear_t`: boolean value to specify a faster index lookup behavior for
+    evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
+    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.
+  - `λ`: The fraction of the interval on either side of each data point that is described by a spline section. Choose
+    in [0, 1.0], defaults to 0.25
+"""
+struct SmoothedLinearInterpolation{uType, tType, IType, pType, T} <:
+       AbstractInterpolation{T}
+    u::uType
+    t::tType
+    I::IType
+    p::SmoothedLinearParameterCache{pType}
+    extrapolation_left::ExtrapolationType.T
+    extrapolation_right::ExtrapolationType.T
+    iguesser::Guesser{tType}
+    cache_parameters::Bool
+    linear_lookup::Bool
+    function SmoothedLinearInterpolation(
+            u, t, I, p, extrapolation_left, extrapolation_right,
+            cache_parameters, assume_linear_t)
+        linear_lookup = seems_linear(assume_linear_t, t)
+        new{typeof(u), typeof(t), typeof(I), typeof(p.slope), eltype(u)}(
+            u, t, I, p, extrapolation_left, extrapolation_right,
+            Guesser(t), cache_parameters, linear_lookup)
+    end
+end
+
+function SmoothedLinearInterpolation(
+        u, t; extrapolation::ExtrapolationType.T = ExtrapolationType.None,
+        extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
+        extrapolation_right::ExtrapolationType.T = ExtrapolationType.None,
+        cache_parameters = true, assume_linear_t = 1e-2, λ = 0.25)
+    extrapolation_left,
+    extrapolation_right = munge_extrapolation(
+        extrapolation, extrapolation_left, extrapolation_right)
+    u, t = munge_data(u, t)
+    p = SmoothedLinearParameterCache(u, t, λ, cache_parameters)
+    return SmoothedLinearInterpolation(
+        u, t, nothing, p, extrapolation_left,
+        extrapolation_right, cache_parameters, assume_linear_t)
+end
+
 """
     QuadraticInterpolation(u, t, mode = :Forward; extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
         extrapolation::ExtrapolationType.T = ExtrapolationType.None, extrapolation_right::ExtrapolationType.T = ExtrapolationType.None, 

src/interpolation_methods.jl

@@ -129,6 +129,27 @@ function _interpolate(A::LinearInterpolation{<:AbstractArray}, t::Number, iguess
     return A.u[ax..., idx] + slope * Δt
 end
 
+# Smoothed linear interpolation
+function _interpolate(A::SmoothedLinearInterpolation, t::Number, iguess)
+    idx = get_idx(A, t, iguess)
+
+    # Check against boundary points between linear and spline interpolation
+    left = (t < A.p.t_tilde[2idx])
+    right = (t > A.p.t_tilde[2idx + 1])
+
+    # Spline interpolation
+    if left && idx != 1
+        return U(A, t, idx)
+    end
+
+    if right && idx != length(A.u) - 1
+        return U(A, t, idx + 1)
+    end
+
+    # Linear interpolation
+    return A.p.u_tilde[2 * idx] + A.p.slope[idx + 1] * (t - A.p.t_tilde[2 * idx])
+end
+
 # Quadratic Interpolation
 function _interpolate(A::QuadraticInterpolation, t::Number, iguess)
     idx = get_idx(A, t, iguess)

src/interpolation_utils.jl

@@ -191,8 +191,7 @@ function cumulative_integral(A::AbstractInterpolation{<:Number}, cache_parameter
     Base.require_one_based_indexing(A.u)
     idxs = cache_parameters ? (1:(length(A.t) - 1)) : (1:0)
     return cumsum(_integral(A, idx, t1, t2)
-    for (idx, t1, t2) in
-        zip(idxs, @view(A.t[begin:(end - 1)]), @view(A.t[(begin + 1):end])))
+    for (idx, t1, t2) in zip(idxs, @view(A.t[begin:(end - 1)]), @view(A.t[(begin + 1):end])))
 end
 
 function get_parameters(A::LinearInterpolation, idx)
@@ -415,3 +414,73 @@ function get_transition_ts(A::SmoothedConstantInterpolation)
 end
 
 get_transition_ts(A::AbstractInterpolation) = A.t
+
+# Compute the spline parametrization value s ∈ [0,1]
+function S(A::SmoothedLinearInterpolation, t, idx)
+    (; Δt, ΔΔt, degenerate_ΔΔt, t_tilde, λ) = A.p
+    Δtᵢ = Δt[idx]
+    ΔΔtᵢ = ΔΔt[idx]
+    tdiff = t - t_tilde[2 * idx - 1]
+    @assert tdiff >= 0
+
+    s = if degenerate_ΔΔt[idx]
+        # Degenerate case Δtᵢ₊₁ ≈ Δtᵢ
+        tdiff / (λ * Δtᵢ)
+    else
+        (-Δtᵢ + sqrt(Δtᵢ^2 + 2 * ΔΔtᵢ * tdiff / λ)) / ΔΔtᵢ
+    end
+    ε = 1e-5
+    @assert -ε<=s<=1+ε "s should be in [0,1], got $s."
+    return s
+end
+
+# Compute the value of a spline section of SmoothedLinearInterpolation
+# as a function of t
+function U(A::SmoothedLinearInterpolation, t, idx)
+    s = S(A, t, idx)
+    return U_s(A, s, idx)
+end
+
+# Compute the value of a spline section of SmoothedLinearInterpolation
+# as a function of the spline parametrization value s
+function U_s(A::SmoothedLinearInterpolation, s, idx)
+    (; Δu, ΔΔu, u_tilde, λ) = A.p
+    Δuᵢ = Δu[idx]
+    ΔΔuᵢ = ΔΔu[idx]
+    return λ / 2 * ΔΔuᵢ * s^2 + λ * Δuᵢ * s + u_tilde[2 * idx - 1]
+end
+
+# Compute the derivative of a spline section of SmoothedLinearInterpolation
+# as a function of t
+function U_deriv(A::SmoothedLinearInterpolation, t, idx)
+    s = S(A, t, idx)
+    s_deriv = S_deriv(A, t, idx)
+    return U_s_deriv(A, s, idx) * s_deriv
+end
+
+# Compute the derivative of the spline parametrization value s
+# with respect to t
+function S_deriv(A::SmoothedLinearInterpolation, t, idx)
+    (; Δt, ΔΔt, degenerate_ΔΔt, t_tilde, λ) = A.p
+    Δtᵢ = Δt[idx]
+    ΔΔtᵢ = ΔΔt[idx]
+    tdiff = t - t_tilde[2 * idx - 1]
+    @assert tdiff >= 0
+
+    if degenerate_ΔΔt[idx]
+        # Degenerate case Δtᵢ₊₁ ≈ Δtᵢ
+        s_deriv = 1 / (λ * Δtᵢ)
+    else
+        s_deriv = 1 / (λ * sqrt(Δtᵢ^2 + 2 * ΔΔtᵢ * tdiff / λ))
+    end
+    return s_deriv
+end
+
+# Compute the derivative of a spline section of SmoothedLinearInterpolation
+# as a function of the spline parametrization value s
+function U_s_deriv(A::AbstractInterpolation, s, idx)
+    (; Δu, ΔΔu, λ) = A.p
+    Δuᵢ = Δu[idx]
+    ΔΔuᵢ = ΔΔu[idx]
+    return λ * ΔΔuᵢ * s + λ * Δuᵢ
+end

src/parameter_caches.jl

@@ -32,6 +32,48 @@ function linear_interpolation_parameters(u::AbstractArray{T, N}, t, idx) where {
     return slope
 end
 
+struct SmoothedLinearParameterCache{uType, tType, λType}
+    Δu::uType
+    Δt::tType
+    ΔΔu::uType
+    ΔΔt::tType
+    u_tilde::uType
+    t_tilde::tType
+    slope::uType
+    # Whether ΔΔt is sufficiently close to 0
+    degenerate_ΔΔt::Vector{Bool}
+    λ::λType
+end
+
+function get_spline_ends(u, Δu, λ)
+    u_tilde = zeros(2 * length(u))
+    u_tilde[1] = u[1]
+    u_tilde[2:2:(end - 1)] = u[1:(end - 1)] .+ (λ / 2) .* Δu[2:(end - 1)]
+    u_tilde[3:2:end] = u[2:end] .- (λ / 2) .* Δu[2:(end - 1)]
+    u_tilde[end] = u[end]
+    return u_tilde
+end
+
+function SmoothedLinearParameterCache(u::AbstractVector, t, λ, cache_parameters)
+    @assert cache_parameters "Parameter caching is mandatory for SmoothedLinearInterpolation"
+
+    Δu = diff(u)
+    Δt = diff(t)
+    pushfirst!(Δt, last(Δt))
+    push!(Δt, first(Δt))
+    pushfirst!(Δu, last(Δu))
+    push!(Δu, first(Δu))
+    ΔΔu = diff(Δu)
+    ΔΔt = diff(Δt)
+    u_tilde = get_spline_ends(u, Δu, λ)
+    t_tilde = get_spline_ends(t, Δt, λ)
+    slope = Δu ./ Δt
+    degenerate_ΔΔt = collect(isapprox.(ΔΔt, 0, atol = 1e-5))
+
+    return SmoothedLinearParameterCache(
+        Δu, Δt, ΔΔu, ΔΔt, u_tilde, t_tilde, slope, degenerate_ΔΔt, λ)
+end
+
 struct SmoothedConstantParameterCache{dType, cType}
     d::dType
     c::cType

test/derivative_tests.jl

@@ -106,6 +106,16 @@ end
         LinearInterpolation; args = [u, t], name = "Linear Interpolation with two points")
 end
 
+@testset "Smoothed Linear Interpolation" begin
+    u = [3.8, 6.7, 1.8, 2.3]
+    t = [1.05, 2.6, 6.7, 8.9]
+
+    for λ in [0.0, 0.25, 0.5]
+        test_derivatives(
+            SmoothedLinearInterpolation, args = [u, t], name = "Smoothed Linear Interpolation (λ = $λ)")
+    end
+end
+
 @testset "Quadratic Interpolation" begin
     u = [1.0, 4.0, 9.0, 16.0]
     t = [1.0, 2.0, 3.0, 4.0]

test/interpolation_tests.jl

@@ -239,6 +239,20 @@ end
     @test_throws DataInterpolations.LeftExtrapolationError A([-1.0, 11.0])
 end
 
+@testset "Smoothed Linear Interpolation" begin
+    test_interpolation_type(SmoothedLinearInterpolation)
+
+    u = [3.8, 6.7, 1.8, 2.3]
+    t = [1.0, 2.6, 6.7, 8.9]
+    A = SmoothedLinearInterpolation(u, t)
+
+    @test A.p.t_tilde ≈ [1.0, 1.2, 2.4, 3.1125, 6.1875, 6.975, 8.625, 8.9]
+    @test A.p.u_tilde ≈ [3.8, 4.1625, 6.3375, 6.0875, 2.4125, 1.8625, 2.2375, 2.3]
+
+    t_eval = A.p.t_tilde[1:(end - 1)] + diff(A.p.t_tilde) / 2
+    @test A(t_eval)≈[3.98125, 5.25, 6.41932, 4.25, 2.01298, 2.05, 2.26875] rtol=1e-4
+end
+
 @testset "Quadratic Interpolation" begin
     test_interpolation_type(QuadraticInterpolation)