Merge pull request #386 from SouthEndMusic/fix_extrapolation_typo

Fix integrating over interval that is fully outside data range

Author: ChrisRackauckas

src/integrals.jl

@@ -24,7 +24,7 @@ function integral(A::AbstractInterpolation, t1::Number, t2::Number)
     if t1 < first(A.t)
         if t2 < first(A.t)
             # If interval is entirely below data
-            return _extrapolate_integral_left(A, t2) - extrapolate_integral_left(A.t1)
+            return _extrapolate_integral_left(A, t1) - _extrapolate_integral_left(A, t2)
         end
 
         idx1 -= 1 # Make sure lowest complete interval is included
@@ -37,7 +37,7 @@ function integral(A::AbstractInterpolation, t1::Number, t2::Number)
     if t2 > last(A.t)
         if t1 > last(A.t)
             # If interval is entirely above data
-            return _extrapolate_integral_right(A, t2) - extrapolate_integral_right(A.t, t1)
+            return _extrapolate_integral_right(A, t2) - _extrapolate_integral_right(A, t1)
         end
 
         idx2 += 1 # Make sure highest complete interval is included

test/integral_tests.jl

@@ -47,6 +47,15 @@ function test_integral(method; args = [], kwargs = [], name::String)
         qint, err = quadgk(func, (t1 + t2) / 2, t2 + 5.0; atol = 1e-12, rtol = 1e-12)
         aint = integral(func, (t1 + t2) / 2, t2 + 5.0)
         @test isapprox(qint, aint, atol = 1e-6, rtol = 1e-8)
+
+        # Integrate intervals fully outside data
+        qint, err = quadgk(func, t1 - 5.0, t1 - 2.5; atol = 1e-12, rtol = 1e-12)
+        aint = integral(func, t1 - 5.0, t1 - 2.5)
+        @test isapprox(qint, aint, atol = 1e-6, rtol = 1e-8)
+
+        qint, err = quadgk(func, t2 + 2.5, t2 + 5.0; atol = 1e-12, rtol = 1e-12)
+        aint = integral(func, t2 + 2.5, t2 + 5.0)
+        @test isapprox(qint, aint, atol = 1e-6, rtol = 1e-8)
     end
     func = method(args...; kwargs...)
     @test_throws DataInterpolations.LeftExtrapolationError integral(func, t[1] - 1.0)

test/interface.jl

@@ -57,7 +57,7 @@ end
 end
 
 @testset "Output Type" begin
-        # Test consistency between eltype(u) and type of the output
+    # Test consistency between eltype(u) and type of the output
     u = Float32[-0.676367f0, 0.8449812f0, 1.2366607f0, -0.13347931f0, 1.9928657f0,
         -0.63596356f0, 0.76009744f0, -0.30632544f0, 0.34649512f0, -0.3846099f0]
     t = 0.1f0:0.1f0:1.0f0