Add SmoothedConstantInterpolation (#367)

* POC

* Add docstring

* derivatives

* Add to docs

* integrals

* Add tests

* Merge branch 'master' into SmoothedConstantInterpolation

* Formatting

* Make periodic extrapolation smooth

* Fix tests

* Actually fix all tests

* Actually Actually fix all tests

* Actually actually actually fix all tests

* Address comments

* Merge cleanup

* Add get_transition_ts

* Formatting

* Update extrapolation behavior

* Formatting

* Fix extrapolation derivatives

* Nit

* Use get_transition_ts in derivative testing

* Formatting

* Fix derivatives and integrals

* Clarify when the last u is ignored

* Format

* Format with JuliaFormatter v1

Author: SouthEndMusic

README.md

@@ -48,6 +48,7 @@ corresponding to `(u,t)` pairs.
 
   - `ConstantInterpolation(u,t)` - A piecewise constant interpolation.
 
+  - `SmoothedConstantInterpolation(u,t)` - An integral preserving continuously differentiable approximation of constant interpolation.
   - `LinearInterpolation(u,t)` - A linear interpolation.
   - `QuadraticInterpolation(u,t)` - A quadratic interpolation.
   - `LagrangeInterpolation(u,t,n)` - A Lagrange interpolation of order `n`.

docs/src/index.md

@@ -18,6 +18,7 @@ corresponding to `(u,t)` pairs.
 
   - `ConstantInterpolation(u,t)` - A piecewise constant interpolation.
 
+  - `SmoothedConstantInterpolation(u,t)` - An integral preserving continuously differentiable approximation of constant interpolation.
   - `LinearInterpolation(u,t)` - A linear interpolation.
   - `QuadraticInterpolation(u,t)` - A quadratic interpolation.
   - `LagrangeInterpolation(u,t,n)` - A Lagrange interpolation of order `n`.

docs/src/manual.md

@@ -6,6 +6,7 @@ QuadraticInterpolation
 LagrangeInterpolation
 AkimaInterpolation
 ConstantInterpolation
+SmoothedConstantInterpolation
 QuadraticSpline
 CubicSpline
 BSplineInterpolation

docs/src/methods.md

@@ -77,6 +77,24 @@ A = ConstantInterpolation(u, t, dir = :right)
 plot(A)
 ```
 
+## Smoothed Constant Interpolation
+
+This function is much like the constant interpolation above, but the transition
+between consecutive values is smoothed out so that the function is continuously
+differentiable. The smoothing is done in such a way that the integral of this function
+is never much off from the same integral of constant interpolation without smoothing (because of the symmetry of the smoothing sections).
+The maximum smoothing distance in the `t` direction from the data points can be set
+with the keyword argument `d_max`.
+
+```@example tutorial
+A = ConstantInterpolation(u, t)
+plot(A)
+A = SmoothedConstantInterpolation(u, t; d_max = 10)
+plot!(A)
+```
+
+Note that `u[end]` is ignored, except when using extrapolation types `Constant` or `Extension`.
+
 ## Quadratic Spline
 
 This is the quadratic spline. It is a continuously differentiable interpolation

ext/DataInterpolationsRegularizationToolsExt.jl

@@ -74,13 +74,15 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Abstrac
         extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
         extrapolation_right::ExtrapolationType.T = ExtrapolationType.None,
         cache_parameters::Bool = false)
-    extrapolation_left, extrapolation_right = munge_extrapolation(
+    extrapolation_left,
+    extrapolation_right = munge_extrapolation(
         extrapolation, extrapolation_left, extrapolation_right)
     u, t = munge_data(u, t; check_sorted = t̂, sorted_arg_name = ("third", "t̂"))
     M = _mapping_matrix(t̂, t)
     Wls½ = LA.diagm(sqrt.(wls))
     Wr½ = LA.diagm(sqrt.(wr))
-    û, λ, Aitp = _reg_smooth_solve(
+    û, λ,
+    Aitp = _reg_smooth_solve(
         u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left,
         extrapolation_right, cache_parameters)
     RegularizationSmooth(
@@ -100,15 +102,17 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, d::Int = 2;
         extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
         extrapolation_right::ExtrapolationType.T = ExtrapolationType.None,
         cache_parameters::Bool = false)
-    extrapolation_left, extrapolation_right = munge_extrapolation(
+    extrapolation_left,
+    extrapolation_right = munge_extrapolation(
         extrapolation, extrapolation_left, extrapolation_right)
     u, t = munge_data(u, t)
     t̂ = t
     N = length(t)
     M = Array{Float64}(LA.I, N, N)
     Wls½ = Array{Float64}(LA.I, N, N)
     Wr½ = Array{Float64}(LA.I, N - d, N - d)
-    û, λ, Aitp = _reg_smooth_solve(
+    û, λ,
+    Aitp = _reg_smooth_solve(
         u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left,
         extrapolation_right, cache_parameters)
     RegularizationSmooth(u,
@@ -137,14 +141,16 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Abstrac
         extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
         extrapolation_right::ExtrapolationType.T = ExtrapolationType.None,
         cache_parameters::Bool = false)
-    extrapolation_left, extrapolation_right = munge_extrapolation(
+    extrapolation_left,
+    extrapolation_right = munge_extrapolation(
         extrapolation, extrapolation_left, extrapolation_right)
     u, t = munge_data(u, t; check_sorted = t̂, sorted_arg_name = ("third", "t̂"))
     N, N̂ = length(t), length(t̂)
     M = _mapping_matrix(t̂, t)
     Wls½ = Array{Float64}(LA.I, N, N)
     Wr½ = Array{Float64}(LA.I, N̂ - d, N̂ - d)
-    û, λ, Aitp = _reg_smooth_solve(
+    û, λ,
+    Aitp = _reg_smooth_solve(
         u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left,
         extrapolation_right, cache_parameters)
     RegularizationSmooth(u,
@@ -174,14 +180,16 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Abstrac
         extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
         extrapolation_right::ExtrapolationType.T = ExtrapolationType.None,
         cache_parameters::Bool = false)
-    extrapolation_left, extrapolation_right = munge_extrapolation(
+    extrapolation_left,
+    extrapolation_right = munge_extrapolation(
         extrapolation, extrapolation_left, extrapolation_right)
     u, t = munge_data(u, t; check_sorted = t̂, sorted_arg_name = ("third", "t̂"))
     N, N̂ = length(t), length(t̂)
     M = _mapping_matrix(t̂, t)
     Wls½ = LA.diagm(sqrt.(wls))
     Wr½ = Array{Float64}(LA.I, N̂ - d, N̂ - d)
-    û, λ, Aitp = _reg_smooth_solve(
+    û, λ,
+    Aitp = _reg_smooth_solve(
         u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left,
         extrapolation_right, cache_parameters)
     RegularizationSmooth(u,
@@ -212,15 +220,17 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing
         extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
         extrapolation_right::ExtrapolationType.T = ExtrapolationType.None,
         cache_parameters::Bool = false)
-    extrapolation_left, extrapolation_right = munge_extrapolation(
+    extrapolation_left,
+    extrapolation_right = munge_extrapolation(
         extrapolation, extrapolation_left, extrapolation_right)
     u, t = munge_data(u, t)
     t̂ = t
     N = length(t)
     M = Array{Float64}(LA.I, N, N)
     Wls½ = LA.diagm(sqrt.(wls))
     Wr½ = Array{Float64}(LA.I, N - d, N - d)
-    û, λ, Aitp = _reg_smooth_solve(
+    û, λ,
+    Aitp = _reg_smooth_solve(
         u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left,
         extrapolation_right, cache_parameters)
     RegularizationSmooth(u,
@@ -251,15 +261,17 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing
         extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
         extrapolation_right::ExtrapolationType.T = ExtrapolationType.None,
         cache_parameters::Bool = false)
-    extrapolation_left, extrapolation_right = munge_extrapolation(
+    extrapolation_left,
+    extrapolation_right = munge_extrapolation(
         extrapolation, extrapolation_left, extrapolation_right)
     u, t = munge_data(u, t)
     t̂ = t
     N = length(t)
     M = Array{Float64}(LA.I, N, N)
     Wls½ = LA.diagm(sqrt.(wls))
     Wr½ = LA.diagm(sqrt.(wr))
-    û, λ, Aitp = _reg_smooth_solve(
+    û, λ,
+    Aitp = _reg_smooth_solve(
         u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left,
         extrapolation_right, cache_parameters)
     RegularizationSmooth(u,
@@ -289,7 +301,8 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing
         extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
         extrapolation_right::ExtrapolationType.T = ExtrapolationType.None,
         cache_parameters::Bool = false)
-    extrapolation_left, extrapolation_right = munge_extrapolation(
+    extrapolation_left,
+    extrapolation_right = munge_extrapolation(
         extrapolation, extrapolation_left, extrapolation_right)
     u, t = munge_data(u, t)
     t̂ = t
@@ -298,7 +311,8 @@ function RegularizationSmooth(u::AbstractVector, t::AbstractVector, t̂::Nothing
     wls, wr = _weighting_by_kw(t, d, wls)
     Wls½ = LA.diagm(sqrt.(wls))
     Wr½ = LA.diagm(sqrt.(wr))
-    û, λ, Aitp = _reg_smooth_solve(
+    û, λ,
+    Aitp = _reg_smooth_solve(
         u, t̂, d, M, Wls½, Wr½, λ, alg, extrapolation_left,
         extrapolation_right, cache_parameters)
     RegularizationSmooth(u,

src/DataInterpolations.jl

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

src/derivatives.jl

@@ -74,12 +74,43 @@ function _extrapolate_derivative_right(A, t, order)
     end
 end
 
+function _extrapolate_derivative_right(A::SmoothedConstantInterpolation, t, order)
+    if A.extrapolation_right == ExtrapolationType.None
+        throw(RightExtrapolationError())
+    elseif A.extrapolation_right in (
+        ExtrapolationType.Constant, ExtrapolationType.Extension)
+        d = min(A.t[end] - A.t[end - 1], 2A.d_max) / 2
+        if A.t[end] + d < t
+            zero(eltype(A.u))
+        else
+            c = (A.u[end] - A.u[end - 1]) / 2
+            -c * (2((t - A.t[end]) / d) - 2) / d
+        end
+
+    else
+        _extrapolate_other(A, t, A.extrapolation_right)
+    end
+end
+
 function _derivative(A::LinearInterpolation, t::Number, iguess)
     idx = get_idx(A, t, iguess; idx_shift = -1, ub_shift = -1, side = :first)
     slope = get_parameters(A, idx)
     slope
 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)
+
+    if (t - A.t[idx]) < d_lower
+        -2c_lower * ((t - A.t[idx]) / d_lower - 1) / d_lower
+    elseif (A.t[idx + 1] - t) < d_upper
+        2c_upper * (1 - (A.t[idx + 1] - t) / d_upper) / d_upper
+    else
+        zero(c_upper / oneunit(t))
+    end
+end
+
 function _derivative(A::QuadraticInterpolation, t::Number, iguess)
     idx = get_idx(A, t, iguess)
     Δt = t - A.t[idx]

src/integrals.jl

@@ -135,6 +135,50 @@ function _extrapolate_integral_right(A, t)
     end
 end
 
+function _extrapolate_integral_right(A::SmoothedConstantInterpolation, t)
+    (; extrapolation_right) = A
+    if extrapolation_right == ExtrapolationType.None
+        throw(RightExtrapolationError())
+    elseif A.extrapolation_right in (
+        ExtrapolationType.Constant, ExtrapolationType.Extension)
+        d = min(A.t[end] - A.t[end - 1], 2A.d_max) / 2
+        c = (A.u[end] - A.u[end - 1]) / 2
+
+        Δt_transition = min(t - A.t[end], d)
+        Δt_constant = max(0, t - A.t[end] - d)
+        out = Δt_transition * A.u[end - 1] -
+              c *
+              (((Δt_transition / d)^3) / (3 / d) - ((Δt_transition^2) / d) -
+               Δt_transition) +
+              Δt_constant * A.u[end]
+
+    elseif extrapolation_right == ExtrapolationType.Linear
+        slope = derivative(A, last(A.t))
+        Δt = t - last(A.t)
+        (last(A.u) + slope * Δt / 2) * Δt
+        _extrapolate_other(A, t, A.extrapolation_right)
+    elseif extrapolation_right == ExtrapolationType.Periodic
+        t_, n = transformation_periodic(A, t)
+        out = integral(A, first(A.t), t_)
+        if !iszero(n)
+            out += n * integral(A, first(A.t), last(A.t))
+        end
+        out
+    else
+        # extrapolation_right == ExtrapolationType.Reflective
+        t_, n = transformation_reflective(A, t)
+        out = if iseven(n)
+            integral(A, t_, last(A.t))
+        else
+            integral(A, t_)
+        end
+        if !iszero(n)
+            out += n * integral(A, first(A.t), last(A.t))
+        end
+        out
+    end
+end
+
 function _integral(A::LinearInterpolation{<:AbstractVector{<:Number}},
         idx::Number, t1::Number, t2::Number)
     slope = get_parameters(A, idx)
@@ -154,6 +198,38 @@ function _integral(
     end
 end
 
+function _integral(A::SmoothedConstantInterpolation{<:AbstractVector},
+        idx::Number, t1::Number, t2::Number)
+    d_lower, d_upper, c_lower, c_upper = get_parameters(A, idx)
+
+    bound_lower = A.t[idx] + d_lower
+    bound_upper = A.t[idx + 1] - d_upper
+
+    out = A.u[idx] * (t2 - t1)
+
+    # Fix extrapolation behavior as constant for now
+    if t1 <= first(A.t)
+        t1 = first(A.t)
+    elseif t2 >= last(A.t)
+        t2 = last(A.t)
+    end
+
+    if t1 < bound_lower
+        t2_ = min(t2, bound_lower)
+        out -= c_lower * d_lower *
+               (((t2_ - A.t[idx]) / d_lower - 1)^3 - ((t1 - A.t[idx]) / d_lower - 1)^3) / 3
+    end
+
+    if t2 > bound_upper
+        t1_ = max(t1, bound_upper)
+        out += c_upper * d_upper *
+               ((1 - (A.t[idx + 1] - t2) / d_upper)^3 -
+                (1 - (A.t[idx + 1] - t1_) / d_upper)^3) / 3
+    end
+
+    out
+end
+
 function _integral(A::QuadraticInterpolation{<:AbstractVector{<:Number}},
         idx::Number, t1::Number, t2::Number)
     α, β = get_parameters(A, idx)