Add tests

src/derivatives.jl

@@ -23,6 +23,11 @@ function _derivative(A::SmoothedConstantInterpolation{<:AbstractVector}, t::Numb
     idx = get_idx(A, t, iguess)
     d_lower, d_upper, c_lower, c_upper = get_parameters(A, idx)
 
+    # Fix extrapolation behavior as constant for now
+    if t <= first(A.t) || t >= last(A.t)
+        return zero(c_upper / oneunit(t))
+    end
+
     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

src/integrals.jl

@@ -65,6 +65,13 @@ function _integral(A::SmoothedConstantInterpolation{<:AbstractVector},
 
     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 *

src/interpolation_caches.jl

@@ -308,12 +308,12 @@ drifts far from the integral of constant interpolation. In this interpolation ty
     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 SmoothedConstantInterpolation{uType, tType, dmaxType, IType, pType, T, N} <:
+struct SmoothedConstantInterpolation{uType, tType, IType, dType, cType, dmaxType, T, N} <:
        AbstractInterpolation{T, N}
     u::uType
     t::tType
     I::IType
-    p::SmoothedConstantParameterCache{uType, tType}
+    p::SmoothedConstantParameterCache{dType, cType}
     d_max::dmaxType
     extrapolate::Bool
     iguesser::Guesser{tType}
@@ -323,7 +323,8 @@ struct SmoothedConstantInterpolation{uType, tType, dmaxType, IType, pType, T, N}
             u, t, I, p, d_max, 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(d_max), typeof(I), typeof(p.d), eltype(u), N}(
+        new{typeof(u), typeof(t), typeof(I), typeof(p.d),
+            typeof(p.c), typeof(d_max), eltype(u), N}(
             u, t, I, p, d_max, extrapolate, Guesser(t), cache_parameters, linear_lookup)
     end
 end

src/interpolation_methods.jl

@@ -143,6 +143,14 @@ end
 function _interpolate(A::SmoothedConstantInterpolation{<:AbstractVector}, t::Number, iguess)
     idx = get_idx(A, t, iguess)
     d_lower, d_upper, c_lower, c_upper = get_parameters(A, idx)
+
+    # Fix extrapolation behavior as constant for now
+    if t <= first(A.t)
+        return first(A.u)
+    elseif t >= last(A.t)
+        return A.u[end - 1]
+    end
+
     out = A.u[idx]
 
     if (t - A.t[idx]) < d_lower

src/parameter_caches.jl

@@ -32,9 +32,9 @@ function linear_interpolation_parameters(u::AbstractArray{T, N}, t, idx) where {
     return slope
 end
 
-struct SmoothedConstantParameterCache{uType, tType}
-    d::tType
-    c::uType
+struct SmoothedConstantParameterCache{dType, cType}
+    d::dType
+    c::cType
 end
 
 function SmoothedConstantParameterCache(u, t, cache_parameters, d_max)
@@ -51,7 +51,7 @@ end
 function smoothed_linear_interpolation_parameters(u, t, d_max, idx)
     # TODO: Add support for making periodic extrapolation smooth
     if isone(idx) || (idx == length(t))
-        zero(eltype(t)), zero(eltype(u))
+        zero(one(eltype(t))) / 2, zero(one(eltype(u)) / 2)
     else
         min(t[idx] - t[idx - 1], t[idx + 1] - t[idx], 2d_max) / 2, (u[idx] - u[idx - 1]) / 2
     end

test/derivative_tests.jl

@@ -27,7 +27,7 @@ function test_derivatives(method; args = [], kwargs = [], name::String)
         # Interpolation time points
         for _t in t[2:(end - 1)]
             if func isa BSplineInterpolation || func isa BSplineApprox ||
-               func isa CubicHermiteSpline
+               func isa CubicHermiteSpline || func isa SmoothedConstantInterpolation
                 fdiff = forward_fdm(5, 1; geom = true)(func, _t)
                 fdiff2 = forward_fdm(5, 1; geom = true)(t -> derivative(func, t), _t)
             else
@@ -140,6 +140,13 @@ end
     @test all(derivative.(Ref(A), t2 .+ 0.1) .== 0.0)
 end
 
+@testset "SmoothedConstantInterpolation" begin
+    u = [5.5, 2.7, 5.1, 3.0]
+    t = [2.5, 5.6, 6.3, 8.9]
+    test_derivatives(SmoothedConstantInterpolation; args = [u, t],
+        name = "Smoothed constant interpolation")
+end
+
 @testset "Quadratic Spline" begin
     u = [0.0, 1.0, 3.0]
     t = [-1.0, 0.0, 1.0]

test/integral_tests.jl

@@ -65,6 +65,20 @@ end
         name = "Linear Interpolation (Vector) with random points")
 end
 
+@testset "SmoothedConstantInterpolation" begin
+    u = [1.0, 4.0, 9.0, 16.0]
+    t = [1.0, 2.0, 3.0, 4.0]
+    test_integral(
+        SmoothedConstantInterpolation; args = [u, t], name = "Smoothed constant interpolation")
+
+    A_constant = ConstantInterpolation(u, t)
+    I_ref = DataInterpolations.integral(A_constant, first(t), last(t))
+    I_smoothed = [DataInterpolations.integral(
+                      SmoothedConstantInterpolation(u, t; d_max), first(t), last(t))
+                  for d_max in 0.0:0.1:1.0]
+    @test all(I_smoothed .≈ I_ref)
+end
+
 @testset "QuadraticInterpolation" begin
     u = [1.0, 4.0, 9.0, 16.0]
     t = [1.0, 2.0, 3.0, 4.0]

test/interpolation_tests.jl

@@ -501,6 +501,19 @@ end
     @test A(-Inf) == first(u)
 end
 
+@testset "Smoothed constant Interpolation" begin
+    test_interpolation_type(SmoothedConstantInterpolation)
+    u = [0.0, 2.0, 1.0, 3.0]
+    t = [1.2, 2.5, 5.7, 8.7]
+    d_max = 0.5
+    A = SmoothedConstantInterpolation(u, t; d_max)
+    test_cached_index(A)
+
+    @test A(1.9) == u[1]
+    @test A(3.1) == u[2]
+    @test A(2.5) ≈ (u[1] + u[2]) / 2
+end
+
 @testset "QuadraticSpline Interpolation" begin
     test_interpolation_type(QuadraticSpline)
 

test/parameter_tests.jl

@@ -7,6 +7,14 @@ using DataInterpolations
     @test A.p.slope ≈ [4.0, -2.0, 1.0, 0.0]
 end
 
+@testset "Smoothed constant Interpolation" begin
+    u = [1.0, 5.0, 3.0, 4.0, 4.0]
+    t = collect(1:5)
+    A = SmoothedConstantInterpolation(u, t; cache_parameters = true)
+    A.p.d ≈ [0.0, 0.5, 0.5, 0.5, 0.0]
+    A.p.c ≈ [0.0, 2.0, -1.0, 0.5, 0.0]
+end
+
 @testset "Quadratic Interpolation" begin
     u = [1.0, 5.0, 3.0, 4.0, 4.0]
     t = collect(1:5)