Merge pull request #393 from SciML/dw/cumulative_integral

Improve `cumulative_integral`

src/interpolation_utils.jl

@@ -189,14 +189,12 @@ function get_idx(A::AbstractInterpolation, t, iguess::Union{<:Integer, Guesser};
     end
 end
 
-function cumulative_integral(A, cache_parameters)
-    if cache_parameters && hasmethod(_integral, Tuple{typeof(A), Number, Number, Number})
-        integral_values = _integral.(
-            Ref(A), 1:(length(A.t) - 1), A.t[1:(end - 1)], A.t[2:end])
-        cumsum(integral_values)
-    else
-        promote_type(eltype(A.u), eltype(A.t))[]
-    end
+cumulative_integral(::AbstractInterpolation, ::Bool) = nothing
+function cumulative_integral(A::AbstractInterpolation{<:Number}, cache_parameters::Bool)
+    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])))
 end
 
 function get_parameters(A::LinearInterpolation, idx)

test/integral_tests.jl

@@ -4,6 +4,7 @@ using DataInterpolations: integral
 using Optim, ForwardDiff
 using RegularizationTools
 using StableRNGs
+using Unitful
 
 function test_integral(method; args = [], kwargs = [], name::String)
     func = method(args...; kwargs..., extrapolation_left = ExtrapolationType.Extension,
@@ -213,3 +214,22 @@ end
     @test_throws DataInterpolations.IntegralNotFoundError integral(A, 1.0, 100.0)
     @test_throws DataInterpolations.IntegralNotFoundError integral(A, 50.0)
 end
+
+# issue #385
+@testset "Integrals with unitful numbers" begin
+    u = rand(5)u"m"
+    A = ConstantInterpolation(u, (1:5)u"s")
+    @test @inferred(integral(A, 4u"s")) ≈ sum(u[1:3]) * u"s"
+end
+
+@testset "cumulative_integral" begin
+    A = ConstantInterpolation(["A", "B", "C"], [0.0, 0.25, 0.75])
+    for cache_parameter in (true, false)
+        @test @inferred(DataInterpolations.cumulative_integral(A, cache_parameter)) ===
+              nothing
+    end
+
+    A = ConstantInterpolation([3.1, 2.5, 4.7], [0.0, 0.25, 0.75])
+    @test @inferred(DataInterpolations.cumulative_integral(A, false)) == Float64[]
+    @test @inferred(DataInterpolations.cumulative_integral(A, true)) == [0.775, 2.025]
+end