Merge pull request #263 from sathvikbhagavan/sb/integral

refactor: integral function to correctly compute second index and remove `samples`
SciML · Jun 7, 2024 · 2464e55 · 2464e55
2 parents daebb52 + 031a4b2
commit 2464e55
Show file tree

Hide file tree

Showing 2 changed files with 42 additions and 14 deletions.
diff --git a/src/integrals.jl b/src/integrals.jl
@@ -6,11 +6,10 @@ end
 function integral(A::AbstractInterpolation, t1::Number, t2::Number)
     ((t1 < A.t[1] || t1 > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError())
     ((t2 < A.t[1] || t2 > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError())
-    bw, fw = samples(A)
     # the index less than or equal to t1
-    idx1 = max(1 + bw, min(searchsortedlast(A.t, t1), length(A.t) - fw))
+    idx1 = max(1, min(searchsortedlast(A.t, t1), length(A.t) - 1))
     # the index less than t2
-    idx2 = max(2 + bw, min(searchsortedlast(A.t, t2), length(A.t) - fw))
+    idx2 = max(1, min(searchsortedlast(A.t, t2), length(A.t) - 1))
     if A.t[idx2] == t2
         idx2 -= 1
     end
@@ -23,7 +22,6 @@ function integral(A::AbstractInterpolation, t1::Number, t2::Number)
     total
 end
 
-samples(A::LinearInterpolation{<:AbstractVector}) = (0, 1)
 function _integral(A::LinearInterpolation{<:AbstractVector{<:Number}},
         idx::Number,
         t::Number)
@@ -34,7 +32,6 @@ function _integral(A::LinearInterpolation{<:AbstractVector{<:Number}},
     t^2 * (u1 - u2) / (2 * t1 - 2 * t2) + t * (t1 * u2 - t2 * u1) / (t1 - t2)
 end
 
-samples(A::ConstantInterpolation{<:AbstractVector}) = (0, 1)
 function _integral(A::ConstantInterpolation{<:AbstractVector}, idx::Number, t::Number)
     if A.dir === :left
         # :left means that value to the left is used for interpolation
@@ -45,7 +42,6 @@ function _integral(A::ConstantInterpolation{<:AbstractVector}, idx::Number, t::N
     end
 end
 
-samples(A::QuadraticInterpolation{<:AbstractVector}) = (0, 1)
 function _integral(A::QuadraticInterpolation{<:AbstractVector{<:Number}},
         idx::Number,
         t::Number)
@@ -69,7 +65,6 @@ function _integral(A::QuadraticInterpolation{<:AbstractVector{<:Number}},
      (t1^2 * t2 - t1^2 * t3 - t1 * t2^2 + t1 * t3^2 + t2^2 * t3 - t2 * t3^2))
 end
 
-samples(A::QuadraticSpline{<:AbstractVector{<:Number}}) = (0, 1)
 function _integral(A::QuadraticSpline{<:AbstractVector{<:Number}}, idx::Number, t::Number)
     t1 = A.t[idx]
     t2 = A.t[idx + 1]
@@ -81,7 +76,6 @@ function _integral(A::QuadraticSpline{<:AbstractVector{<:Number}}, idx::Number,
     (2 * t1 - 2 * t2)
 end
 
-samples(A::CubicSpline{<:AbstractVector{<:Number}}) = (0, 1)
 function _integral(A::CubicSpline{<:AbstractVector{<:Number}}, idx::Number, t::Number)
     t1 = A.t[idx]
     t2 = A.t[idx + 1]
@@ -98,7 +92,6 @@ function _integral(A::CubicSpline{<:AbstractVector{<:Number}}, idx::Number, t::N
      (6 * h2))
 end
 
-samples(A::AkimaInterpolation{<:AbstractVector{<:Number}}) = (0, 1)
 function _integral(A::AkimaInterpolation{<:AbstractVector{<:Number}},
         idx::Number,
         t::Number)

diff --git a/test/integral_tests.jl b/test/integral_tests.jl
@@ -13,21 +13,38 @@ function test_integral(method, u, t; args = [], kwargs = [], name::String)
         # integral(A, t1, t2)
         qint, err = quadgk(func, t1, t2; atol = 1e-12, rtol = 1e-12)
         aint = integral(func, t1, t2)
-        @test isapprox(qint, aint, atol = 1e-8)
+        @test isapprox(qint, aint, atol = 1e-6, rtol = 1e-8)
 
         # integral(A, t)
         qint, err = quadgk(func, t1, (t1 + t2) / 2; atol = 1e-12, rtol = 1e-12)
         aint = integral(func, (t1 + t2) / 2)
-        @test isapprox(qint, aint, atol = 1e-8)
+        @test isapprox(qint, aint, atol = 1e-6, rtol = 1e-8)
+
+        # integral(A, t1, t), integral(A, t, t2), integral(A, t, t)
+        ts = range(t1, t2; length = 100)
+        for t in ts
+            qint, err = quadgk(func, t1, t; atol = 1e-12, rtol = 1e-12)
+            aint1 = integral(func, t1, t)
+            @test isapprox(qint, aint1, atol = 1e-5, rtol = 1e-8)
+            aint2 = integral(func, t)
+            @test aint1 == aint2
+
+            qint, err = quadgk(func, t, t2; atol = 1e-12, rtol = 1e-12)
+            aint = integral(func, t, t2)
+            @test isapprox(qint, aint, atol = 1e-5, rtol = 1e-8)
+
+            aint = integral(func, t, t)
+            @test aint == 0.0
+        end
 
         # integrals with extrapolation
         qint, err = quadgk(func, t1 - 5.0, (t1 + t2) / 2; atol = 1e-12, rtol = 1e-12)
         aint = integral(func, t1 - 5.0, (t1 + t2) / 2)
-        @test isapprox(qint, aint, atol = 1e-8)
+        @test isapprox(qint, aint, atol = 1e-6, rtol = 1e-8)
 
         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-8)
+        @test isapprox(qint, aint, atol = 1e-6, rtol = 1e-8)
     end
     func = method(u, t, args...; kwargs...)
     @test_throws DataInterpolations.ExtrapolationError integral(func, t[1] - 1.0)
@@ -40,6 +57,10 @@ end
     u = 2.0collect(1:10)
     t = 1.0collect(1:10)
     test_integral(LinearInterpolation, u, t; name = "Linear Interpolation (Vector)")
+    u = round.(rand(100), digits = 5)
+    t = 1.0collect(1:100)
+    test_integral(LinearInterpolation, u, t;
+        name = "Linear Interpolation (Vector) with random points")
 end
 
 @testset "QuadraticInterpolation" begin
@@ -53,6 +74,10 @@ end
         t;
         args = [:Backward],
         name = "Quadratic Interpolation (Vector)")
+    u = round.(rand(100), digits = 5)
+    t = 1.0collect(1:10)
+    test_integral(QuadraticInterpolation, u, t;
+        name = "Quadratic Interpolation (Vector) with random points")
 end
 
 @testset "LagrangeInterpolation" begin
@@ -67,18 +92,28 @@ end
     u = [0.0, 1.0, 3.0]
     t = [-1.0, 0.0, 1.0]
     test_integral(QuadraticSpline, u, t; name = "Quadratic Spline (Vector)")
+    u = round.(rand(100), digits = 5)
+    t = 1.0collect(1:100)
+    test_integral(
+        QuadraticSpline, u, t; name = "Quadratic Spline (Vector) with random points")
 end
 
 @testset "CubicSpline" begin
     u = [0.0, 1.0, 3.0]
     t = [-1.0, 0.0, 1.0]
-    test_integral(CubicSpline, u, t; name = "Cubic Spline Interpolation (Vector)")
+    test_integral(CubicSpline, u, t; name = "Cubic Spline (Vector)")
+    u = round.(rand(100), digits = 5)
+    t = 1.0collect(1:100)
+    test_integral(CubicSpline, u, t; name = "Cubic Spline (Vector) with random points")
 end
 
 @testset "AkimaInterpolation" begin
     u = [0.0, 2.0, 1.0, 3.0, 2.0, 6.0, 5.5, 5.5, 2.7, 5.1, 3.0]
     t = collect(0.0:10.0)
     test_integral(AkimaInterpolation, u, t; name = "Akima Interpolation (Vector)")
+    u = round.(rand(100), digits = 5)
+    t = 1.0collect(1:100)
+    test_integral(AkimaInterpolation, u, t; name = "Akima Interpolation (Vector)")
 end
 
 @testset "RegularizationSmooth" begin