Merge pull request #433 from simulutions/configurable-extrapolation-method-in-invert_integral()

Add optional extrapolation methods

Author: SouthEndMusic

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(

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)

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)

test/integral_inverse_tests.jl

@@ -21,6 +21,40 @@ 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, extrapolation_left = ExtrapolationType.Extension,
+        extrapolation_right = 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, extrapolation_left = ExtrapolationType.Extension,
+        extrapolation_right = 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 +63,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 +76,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)

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))