Add optional extrapolation methods

src/integral_inverses.jl

@@ -41,9 +41,9 @@ 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, 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 +57,14 @@ 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,10 @@ 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, 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 +109,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(

test/integral_inverse_tests.jl

@@ -21,6 +21,38 @@ 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(1:4)
+    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)
+
+    # for a linear function, the integral is quadratic
+    # but the constant extrapolation part is linear.
+    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
+
+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]
@@ -29,6 +61,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
@@ -40,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)