Merge pull request #470 from SciML/Fix-integrating-extrapolation-of-SmoothedConstantInterpolation

Fix + test Integrating `SmoothedConstantInterpolation` without smoothing

Author: ChrisRackauckas

docs/make.jl

@@ -26,6 +26,4 @@ makedocs(modules = [DataInterpolations],
         "Inverting Integrals" => "inverting_integrals.md"
     ])
 
-
-
 deploydocs(repo = "github.com/SciML/DataInterpolations.jl"; push_preview = true)

src/integrals.jl

@@ -142,16 +142,18 @@ function _extrapolate_integral_right(A::SmoothedConstantInterpolation, t)
     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]
-
+        out = Δt_constant * A.u[end]
+
+        if !iszero(d)
+            c = (A.u[end] - A.u[end - 1]) / 2
+            Δt_transition = min(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)
+        end
+        out
     elseif extrapolation_right == ExtrapolationType.Linear
         slope = derivative(A, last(A.t))
         Δt = t - last(A.t)

test/integral_tests.jl

@@ -95,6 +95,10 @@ end
                       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)
+
+    A = SmoothedConstantInterpolation(
+        u, t; extrapolation = ExtrapolationType.Constant, d_max = 0)
+    @test DataInterpolations.integral(A, 5.0, 6.0) == 16.0
 end
 
 @testset "QuadraticInterpolation" begin

test/interpolation_tests.jl

@@ -1081,15 +1081,17 @@ end
         u2 = [[u[i], u[i] + 1] for i in eachindex(u)]
         du2 = [[du[i], du[i]] for i in eachindex(du)]
         ddu2 = [[ddu[i], ddu[i]] for i in eachindex(ddu)]
-        A2 = QuinticHermiteSpline(ddu2, du2, u2, t; extrapolation = ExtrapolationType.Extension)
+        A2 = QuinticHermiteSpline(
+            ddu2, du2, u2, t; extrapolation = ExtrapolationType.Extension)
         @test u2 ≈ A2.(t)
         @test A2(100.0) ≈ repeat([10.107996], 2) + [0, 1] rtol=1e-5
         @test A2(300.0) ≈ repeat([11.364162], 2) + [0, 1] rtol=1e-5
         # Test allocation-free interpolation with Vector{StaticArrays.SVector}
         u2_s = [convert(SVector{length(u2[1])}, i) for i in u2]
         du2_s = [convert(SVector{length(du2[1])}, i) for i in du2]
         ddu2_s = [convert(SVector{length(du2[1])}, i) for i in ddu2]
-        A2_s = @inferred(QuinticHermiteSpline(ddu2_s, du2_s, u2_s, t; extrapolation = ExtrapolationType.Extension))
+        A2_s = @inferred(QuinticHermiteSpline(
+            ddu2_s, du2_s, u2_s, t; extrapolation = ExtrapolationType.Extension))
         @test A2_s(100.0) == A2(100.0)
         @test A2_s(300.0) == A2(300.0)
         @test A2_s(0.7) isa SVector{length(u2[1])}

test/runtests.jl

@@ -38,4 +38,4 @@ const GROUP = get(ENV, "GROUP", "All")
     if GROUP == "QA"
         @safetestset "Quality Assurance" include("qa.jl")
     end
-end
+end