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

ext/DataInterpolationsSparseConnectivityTracerExt.jl

@@ -3,26 +3,24 @@ module DataInterpolationsSparseConnectivityTracerExt
 using SparseConnectivityTracer: AbstractTracer, Dual, primal, tracer
 using SparseConnectivityTracer: GradientTracer, gradient_tracer_1_to_1
 using SparseConnectivityTracer: HessianTracer, hessian_tracer_1_to_1
-using FillArrays: Fill # from FillArrays.jl
+using FillArrays: Fill
 using DataInterpolations:
-    AbstractInterpolation,
-    LinearInterpolation,
-    QuadraticInterpolation,
-    LagrangeInterpolation,
-    AkimaInterpolation,
-    ConstantInterpolation,
-    QuadraticSpline,
-    CubicSpline,
-    BSplineInterpolation,
-    BSplineApprox,
-    CubicHermiteSpline,
-    # PCHIPInterpolation,
-    QuinticHermiteSpline,
-    output_size
+                          AbstractInterpolation,
+                          AkimaInterpolation,
+                          BSplineApprox,
+                          BSplineInterpolation,
+                          ConstantInterpolation,
+                          CubicHermiteSpline,
+                          CubicSpline,
+                          LagrangeInterpolation,
+                          LinearInterpolation,
+                          QuadraticInterpolation,
+                          QuadraticSpline,
+                          QuinticHermiteSpline,
+                          SmoothedLinearInterpolation,
+                          output_size
 
-#===========#
 # Utilities #
-#===========#
 
 # Limit support to `u` begin an AbstractVector{<:Number} or AbstractMatrix{<:Number},
 # to avoid any cases where the output size is dependent on the input value.
@@ -33,26 +31,26 @@ function _sct_interpolate(
         uType::Type{<:AbstractVector{<:Number}},
         t::GradientTracer,
         is_der_1_zero,
-        is_der_2_zero,
-    )
+        is_der_2_zero
+)
     return gradient_tracer_1_to_1(t, is_der_1_zero)
 end
 function _sct_interpolate(
         ::AbstractInterpolation,
         uType::Type{<:AbstractVector{<:Number}},
         t::HessianTracer,
         is_der_1_zero,
-        is_der_2_zero,
-    )
+        is_der_2_zero
+)
     return hessian_tracer_1_to_1(t, is_der_1_zero, is_der_2_zero)
 end
 function _sct_interpolate(
         interp::AbstractInterpolation,
         uType::Type{<:AbstractMatrix{<:Number}},
         t::GradientTracer,
         is_der_1_zero,
-        is_der_2_zero,
-    )
+        is_der_2_zero
+)
     t = gradient_tracer_1_to_1(t, is_der_1_zero)
     N = only(output_size(interp))
     return Fill(t, N)
@@ -62,48 +60,47 @@ function _sct_interpolate(
         uType::Type{<:AbstractMatrix{<:Number}},
         t::HessianTracer,
         is_der_1_zero,
-        is_der_2_zero,
-    )
+        is_der_2_zero
+)
     t = hessian_tracer_1_to_1(t, is_der_1_zero, is_der_2_zero)
     N = only(output_size(interp))
     return Fill(t, N)
 end
 
-#===========#
 # Overloads #
-#===========#
 
 # We assume that with the exception of ConstantInterpolation and LinearInterpolation,
 # all interpolations have a non-zero second derivative at some point in the input domain.
 
 for (I, is_der1_zero, is_der2_zero) in (
-        (:ConstantInterpolation, true, true),
-        (:LinearInterpolation, false, true),
-        (:QuadraticInterpolation, false, false),
-        (:LagrangeInterpolation, false, false),
-        (:AkimaInterpolation, false, false),
-        (:QuadraticSpline, false, false),
-        (:CubicSpline, false, false),
-        (:BSplineInterpolation, false, false),
-        (:BSplineApprox, false, false),
-        (:CubicHermiteSpline, false, false),
-        (:QuinticHermiteSpline, false, false),
-    )
+    (:AkimaInterpolation, false, false),
+    (:BSplineApprox, false, false),
+    (:BSplineInterpolation, false, false),
+    (:ConstantInterpolation, true, true),
+    (:CubicHermiteSpline, false, false),
+    (:CubicSpline, false, false),
+    (:LagrangeInterpolation, false, false),
+    (:LinearInterpolation, false, true),
+    (:QuadraticInterpolation, false, false),
+    (:QuadraticSpline, false, false),
+    (:QuinticHermiteSpline, false, false),
+    (:SmoothedLinearInterpolation, false, false)
+)
     @eval function (interp::$(I){uType})(
             t::AbstractTracer
-        ) where {uType <: AbstractArray{<:Number}}
+    ) where {uType <: AbstractArray{<:Number}}
         return _sct_interpolate(interp, uType, t, $is_der1_zero, $is_der2_zero)
     end
 end
 
 # Some Interpolations require custom overloads on `Dual` due to mutation of caches.
 for I in (
-        :LagrangeInterpolation,
-        :BSplineInterpolation,
-        :BSplineApprox,
-        :CubicHermiteSpline,
-        :QuinticHermiteSpline,
-    )
+    :LagrangeInterpolation,
+    :BSplineInterpolation,
+    :BSplineApprox,
+    :CubicHermiteSpline,
+    :QuinticHermiteSpline
+)
     @eval function (interp::$(I){uType})(d::Dual) where {uType <: AbstractVector}
         p = interp(primal(d))
         t = interp(tracer(d))

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)