From d68068e5e4b6c3d10acfabb25324c62503fd674a Mon Sep 17 00:00:00 2001 From: Sathvik Bhagavan Date: Fri, 7 Jun 2024 15:03:34 +0000 Subject: [PATCH 1/2] refactor: integral function to correctly compute second index and remove `samples` --- src/integrals.jl | 11 ++--------- 1 file changed, 2 insertions(+), 9 deletions(-) diff --git a/src/integrals.jl b/src/integrals.jl index e932bc5b..2751eb56 100644 --- 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) From 031a4b271b64c71690df83a410286b34aa67cde0 Mon Sep 17 00:00:00 2001 From: Sathvik Bhagavan Date: Fri, 7 Jun 2024 15:04:04 +0000 Subject: [PATCH 2/2] test: add more integral tests --- test/integral_tests.jl | 45 +++++++++++++++++++++++++++++++++++++----- 1 file changed, 40 insertions(+), 5 deletions(-) diff --git a/test/integral_tests.jl b/test/integral_tests.jl index cb973570..86a5c17f 100644 --- 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