From 23e679879851d2fdcb73b58b7615a399774c2603 Mon Sep 17 00:00:00 2001 From: lexverheem Date: Mon, 26 May 2025 13:56:43 +0200 Subject: [PATCH 1/8] Add optional extrapolation methods --- src/integral_inverses.jl | 26 ++++++++++++++++++-------- 1 file changed, 18 insertions(+), 8 deletions(-) diff --git a/src/integral_inverses.jl b/src/integral_inverses.jl index 075da29f..30cad22f 100644 --- a/src/integral_inverses.jl +++ b/src/integral_inverses.jl @@ -41,9 +41,10 @@ struct LinearInterpolationIntInv{uType, tType, itpType, T} <: extrapolation_right::ExtrapolationType.T iguesser::Guesser{tType} itp::itpType - function LinearInterpolationIntInv(u, t, A) + function LinearInterpolationIntInv(u, t, A, extrapolation_left = A.extrapolation_left, + extrapolation_right = A.extrapolation_right) new{typeof(u), typeof(t), typeof(A), eltype(u)}( - u, t, A.extrapolation_left, A.extrapolation_right, Guesser(t), A) + u, t, extrapolation_left, extrapolation_right, Guesser(t), A) end end @@ -57,9 +58,13 @@ function get_I(A::AbstractInterpolation) I end -function invert_integral(A::LinearInterpolation{<:AbstractVector{<:Number}}) +function invert_integral(A::LinearInterpolation{<:AbstractVector{<:Number}}, + extrapolation_left::ExtrapolationType.T = A.extrapolation_left, + extrapolation_right::ExtrapolationType.T = A.extrapolation_right) !invertible_integral(A) && throw(IntegralNotInvertibleError()) - return LinearInterpolationIntInv(A.t, get_I(A), A) + + return LinearInterpolationIntInv( + A.t, get_I(A), A, extrapolation_left, extrapolation_right) end function _interpolate( @@ -92,9 +97,11 @@ struct ConstantInterpolationIntInv{uType, tType, itpType, T} <: extrapolation_right::ExtrapolationType.T iguesser::Guesser{tType} itp::itpType - function ConstantInterpolationIntInv(u, t, A) + function ConstantInterpolationIntInv( + u, t, A, extrapolation_left = A.extrapolation_left, + extrapolation_right = A.extrapolation_right) new{typeof(u), typeof(t), typeof(A), eltype(u)}( - u, t, A.extrapolation_left, A.extrapolation_right, Guesser(t), A + u, t, extrapolation_left, extrapolation_right, Guesser(t), A ) end end @@ -103,9 +110,12 @@ function invertible_integral(A::ConstantInterpolation{<:AbstractVector{<:Number} return all(A.u .> 0) end -function invert_integral(A::ConstantInterpolation{<:AbstractVector{<:Number}}) +function invert_integral(A::ConstantInterpolation{<:AbstractVector{<:Number}}, + extrapolation_left::ExtrapolationType.T = A.extrapolation_left, + extrapolation_right::ExtrapolationType.T = A.extrapolation_right) !invertible_integral(A) && throw(IntegralNotInvertibleError()) - return ConstantInterpolationIntInv(A.t, get_I(A), A) + return ConstantInterpolationIntInv( + A.t, get_I(A), A, extrapolation_left, extrapolation_right) end function _interpolate( From ab07ff75bb18a7277082d3a940fce90f713c467c Mon Sep 17 00:00:00 2001 From: lexverheem Date: Mon, 26 May 2025 14:37:51 +0200 Subject: [PATCH 2/8] add test for configurable extrapolation method --- test/integral_inverse_tests.jl | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/test/integral_inverse_tests.jl b/test/integral_inverse_tests.jl index 6b09f9cd..b43b2ff7 100644 --- a/test/integral_inverse_tests.jl +++ b/test/integral_inverse_tests.jl @@ -1,6 +1,5 @@ using DataInterpolations using DataInterpolations: integral, derivative, invert_integral -using FiniteDifferences function test_integral_inverses(method; args = [], kwargs = []) A = method(args...; kwargs..., extrapolation = ExtrapolationType.Extension) @@ -21,6 +20,21 @@ function test_integral_inverses(method; args = [], kwargs = []) @test @inferred(A(ts[37])) == A(ts[37]) end +function test_integral_inverse_extrapolation() + # Linear function with constant extrapolation + t = collect(0:4) + u = [0.0, 2.0, 3.0, 4.0] + A = LinearInterpolation(u, t, extrapolation = ExtrapolationType.Constant) + + A_intinv = invert_integral(A, extrapolation_left = ExtrapolationType.Extension, + extrapolation_right = ExtrapolationType.Extension) + + # for a linear function, the integral is quadratic + # but the constant extrapolation part is linear. + A_5 = 1 / 2 * 4.0^2 + 5.0 + @test A_int_inv(A_5) ≈ 5.0 +end + @testset "Linear Interpolation" begin t = collect(1:5) u = [1.0, 1.0, 2.0, 4.0, 3.0] @@ -29,6 +43,8 @@ end u = [1.0, -1.0, 2.0, 4.0, 3.0] A = LinearInterpolation(u, t) @test_throws DataInterpolations.IntegralNotInvertibleError invert_integral(A) + + test_integral_inverse_extrapolation() end @testset "Constant Interpolation" begin From 778b5f5d1180c133721b7a68d16920cfc9f4582c Mon Sep 17 00:00:00 2001 From: lexverheem Date: Mon, 26 May 2025 14:56:32 +0200 Subject: [PATCH 3/8] remove optional from inner contructor --- src/integral_inverses.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/integral_inverses.jl b/src/integral_inverses.jl index 30cad22f..3820dfa5 100644 --- a/src/integral_inverses.jl +++ b/src/integral_inverses.jl @@ -41,8 +41,8 @@ struct LinearInterpolationIntInv{uType, tType, itpType, T} <: extrapolation_right::ExtrapolationType.T iguesser::Guesser{tType} itp::itpType - function LinearInterpolationIntInv(u, t, A, extrapolation_left = A.extrapolation_left, - extrapolation_right = A.extrapolation_right) + function LinearInterpolationIntInv(u, t, A, extrapolation_left, + extrapolation_right) new{typeof(u), typeof(t), typeof(A), eltype(u)}( u, t, extrapolation_left, extrapolation_right, Guesser(t), A) end @@ -98,8 +98,7 @@ struct ConstantInterpolationIntInv{uType, tType, itpType, T} <: iguesser::Guesser{tType} itp::itpType function ConstantInterpolationIntInv( - u, t, A, extrapolation_left = A.extrapolation_left, - extrapolation_right = A.extrapolation_right) + u, t, A, extrapolation_left, extrapolation_right) new{typeof(u), typeof(t), typeof(A), eltype(u)}( u, t, extrapolation_left, extrapolation_right, Guesser(t), A ) From 46298c94a8d277dd6351a36eea8523f23ab35d49 Mon Sep 17 00:00:00 2001 From: lexverheem Date: Mon, 26 May 2025 16:35:31 +0200 Subject: [PATCH 4/8] added unit test --- src/integral_inverses.jl | 3 +-- test/integral_inverse_tests.jl | 15 +++++++++------ 2 files changed, 10 insertions(+), 8 deletions(-) diff --git a/src/integral_inverses.jl b/src/integral_inverses.jl index 3820dfa5..66e66355 100644 --- a/src/integral_inverses.jl +++ b/src/integral_inverses.jl @@ -41,8 +41,7 @@ struct LinearInterpolationIntInv{uType, tType, itpType, T} <: extrapolation_right::ExtrapolationType.T iguesser::Guesser{tType} itp::itpType - function LinearInterpolationIntInv(u, t, A, extrapolation_left, - extrapolation_right) + function LinearInterpolationIntInv(u, t, A, extrapolation_left, extrapolation_right) new{typeof(u), typeof(t), typeof(A), eltype(u)}( u, t, extrapolation_left, extrapolation_right, Guesser(t), A) end diff --git a/test/integral_inverse_tests.jl b/test/integral_inverse_tests.jl index b43b2ff7..b72a57a6 100644 --- a/test/integral_inverse_tests.jl +++ b/test/integral_inverse_tests.jl @@ -1,5 +1,6 @@ using DataInterpolations using DataInterpolations: integral, derivative, invert_integral +using FiniteDifferences function test_integral_inverses(method; args = [], kwargs = []) A = method(args...; kwargs..., extrapolation = ExtrapolationType.Extension) @@ -22,17 +23,19 @@ end function test_integral_inverse_extrapolation() # Linear function with constant extrapolation - t = collect(0:4) - u = [0.0, 2.0, 3.0, 4.0] + t = collect(1:4) + u = [1.0, 2.0, 3.0, 4.0] A = LinearInterpolation(u, t, extrapolation = ExtrapolationType.Constant) - A_intinv = invert_integral(A, extrapolation_left = ExtrapolationType.Extension, - extrapolation_right = ExtrapolationType.Extension) + A_intinv = invert_integral(A, ExtrapolationType.Extension, ExtrapolationType.Extension) # for a linear function, the integral is quadratic # but the constant extrapolation part is linear. - A_5 = 1 / 2 * 4.0^2 + 5.0 - @test A_int_inv(A_5) ≈ 5.0 + area_0_to_4 = 0.5 * 4.0^2 + area_4_to_5 = 4.0 + area = area_0_to_4 + area_4_to_5 + + @test A_intinv(area) ≈ 5.0 end @testset "Linear Interpolation" begin From c709b5e1f3989bcbfbd8cb9f209933088ca9114a Mon Sep 17 00:00:00 2001 From: lexverheem Date: Mon, 26 May 2025 16:46:45 +0200 Subject: [PATCH 5/8] Added proper test cases --- src/integral_inverses.jl | 3 ++- test/integral_inverse_tests.jl | 17 +++++++++++++++++ 2 files changed, 19 insertions(+), 1 deletion(-) diff --git a/src/integral_inverses.jl b/src/integral_inverses.jl index 66e66355..84536a68 100644 --- a/src/integral_inverses.jl +++ b/src/integral_inverses.jl @@ -57,7 +57,8 @@ function get_I(A::AbstractInterpolation) I end -function invert_integral(A::LinearInterpolation{<:AbstractVector{<:Number}}, +function invert_integral( + A::LinearInterpolation{<:AbstractVector{<:Number}}, extrapolation_left::ExtrapolationType.T = A.extrapolation_left, extrapolation_right::ExtrapolationType.T = A.extrapolation_right) !invertible_integral(A) && throw(IntegralNotInvertibleError()) diff --git a/test/integral_inverse_tests.jl b/test/integral_inverse_tests.jl index b72a57a6..bfe58ed8 100644 --- a/test/integral_inverse_tests.jl +++ b/test/integral_inverse_tests.jl @@ -38,6 +38,21 @@ function test_integral_inverse_extrapolation() @test A_intinv(area) ≈ 5.0 end +function test_integral_inverse_const_extrapolation() + # Constant function with constant extrapolation + t = collect(1:4) + u = [1.0, 1.0, 1.0, 1.0] + A = ConstantInterpolation(u, t, extrapolation = ExtrapolationType.Extension) + + A_intinv = invert_integral(A, ExtrapolationType.Extension, ExtrapolationType.Extension) + + area_0_to_4 = 1.0 * (4.0 - 1.0) + area_4_to_5 = 1.0 + area = area_0_to_4 + area_4_to_5 + + @test A_intinv(area) ≈ 5.0 +end + @testset "Linear Interpolation" begin t = collect(1:5) u = [1.0, 1.0, 2.0, 4.0, 3.0] @@ -59,6 +74,8 @@ end u = [1.0, -1.0, 2.0, 4.0, 3.0] A = ConstantInterpolation(u, t) @test_throws DataInterpolations.IntegralNotInvertibleError invert_integral(A) + + test_integral_inverse_const_extrapolation() end t = collect(1:5) From d90876a503b4fc430f0417826d0b5ddbfff0ca0a Mon Sep 17 00:00:00 2001 From: lexverheem Date: Tue, 27 May 2025 09:56:46 +0200 Subject: [PATCH 6/8] add semi colon --- src/integral_inverses.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/integral_inverses.jl b/src/integral_inverses.jl index 84536a68..fc4177ca 100644 --- a/src/integral_inverses.jl +++ b/src/integral_inverses.jl @@ -58,7 +58,7 @@ function get_I(A::AbstractInterpolation) end function invert_integral( - A::LinearInterpolation{<:AbstractVector{<:Number}}, + A::LinearInterpolation{<:AbstractVector{<:Number}}; extrapolation_left::ExtrapolationType.T = A.extrapolation_left, extrapolation_right::ExtrapolationType.T = A.extrapolation_right) !invertible_integral(A) && throw(IntegralNotInvertibleError()) @@ -109,7 +109,7 @@ function invertible_integral(A::ConstantInterpolation{<:AbstractVector{<:Number} return all(A.u .> 0) end -function invert_integral(A::ConstantInterpolation{<:AbstractVector{<:Number}}, +function invert_integral(A::ConstantInterpolation{<:AbstractVector{<:Number}}; extrapolation_left::ExtrapolationType.T = A.extrapolation_left, extrapolation_right::ExtrapolationType.T = A.extrapolation_right) !invertible_integral(A) && throw(IntegralNotInvertibleError()) From 222cf637f4c18556986040e6d30fefb57c0dc3c5 Mon Sep 17 00:00:00 2001 From: lexverheem Date: Tue, 27 May 2025 09:58:29 +0200 Subject: [PATCH 7/8] run formatter --- src/interpolation_caches.jl | 4 ++-- src/interpolation_utils.jl | 3 ++- test/interpolation_tests.jl | 2 ++ 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index 2f2f1d37..47c0773b 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -1220,8 +1220,8 @@ function BSplineApprox( end for k in 2:(n - 1) q[ax_u..., - k] = u[ax_u..., k] - sc[k, 1] * u[ax_u..., 1] - - sc[k, h] * u[ax_u..., end] + k] = u[ax_u..., k] - sc[k, 1] * u[ax_u..., 1] - + sc[k, h] * u[ax_u..., end] end Q = Array{T, N}(undef, size(u)[1:(end - 1)]..., h - 2) for i in 2:(h - 1) diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index 73dd74ce..ff16d56f 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -191,7 +191,8 @@ function cumulative_integral(A::AbstractInterpolation{<:Number}, cache_parameter Base.require_one_based_indexing(A.u) idxs = cache_parameters ? (1:(length(A.t) - 1)) : (1:0) return cumsum(_integral(A, idx, t1, t2) - for (idx, t1, t2) in zip(idxs, @view(A.t[begin:(end - 1)]), @view(A.t[(begin + 1):end]))) + for (idx, t1, t2) in + zip(idxs, @view(A.t[begin:(end - 1)]), @view(A.t[(begin + 1):end]))) end function get_parameters(A::LinearInterpolation, idx) diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl index 7c4e4e44..3677da4f 100644 --- a/test/interpolation_tests.jl +++ b/test/interpolation_tests.jl @@ -1027,6 +1027,7 @@ end ut1 = Float32[0.1, 0.2, 0.3, 0.4, 0.5] ut2 = Float64[0.1, 0.2, 0.3, 0.4, 0.5] for u in (ut1, ut2), t in (ut1, ut2) + interp = @inferred(LinearInterpolation(ut1, ut2)) for xs in (u, t) ys = @inferred(interp(xs)) @@ -1109,6 +1110,7 @@ f_cubic_spline = c -> square(CubicSpline, c) iszero_allocations(u, t) = iszero(@allocated(DataInterpolations.munge_data(u, t))) for T in (String, Union{String, Missing}), dims in 1:3 + _u0 = convert(Array{T}, reshape(u0, ntuple(i -> i == dims ? 3 : 1, dims))) u, t = @inferred(DataInterpolations.munge_data(_u0, t0)) From e4c33066bbd27a19ce546cc95fdd257154a2656b Mon Sep 17 00:00:00 2001 From: lexverheem Date: Tue, 27 May 2025 11:31:01 +0200 Subject: [PATCH 8/8] fix test with keyword args --- test/integral_inverse_tests.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/test/integral_inverse_tests.jl b/test/integral_inverse_tests.jl index bfe58ed8..4974c60a 100644 --- a/test/integral_inverse_tests.jl +++ b/test/integral_inverse_tests.jl @@ -27,7 +27,8 @@ function test_integral_inverse_extrapolation() u = [1.0, 2.0, 3.0, 4.0] A = LinearInterpolation(u, t, extrapolation = ExtrapolationType.Constant) - A_intinv = invert_integral(A, ExtrapolationType.Extension, ExtrapolationType.Extension) + A_intinv = invert_integral(A, extrapolation_left = ExtrapolationType.Extension, + extrapolation_right = ExtrapolationType.Extension) # for a linear function, the integral is quadratic # but the constant extrapolation part is linear. @@ -44,7 +45,8 @@ function test_integral_inverse_const_extrapolation() u = [1.0, 1.0, 1.0, 1.0] A = ConstantInterpolation(u, t, extrapolation = ExtrapolationType.Extension) - A_intinv = invert_integral(A, ExtrapolationType.Extension, ExtrapolationType.Extension) + A_intinv = invert_integral(A, extrapolation_left = ExtrapolationType.Extension, + extrapolation_right = ExtrapolationType.Extension) area_0_to_4 = 1.0 * (4.0 - 1.0) area_4_to_5 = 1.0