Unitful support (#25)

* Added Unitful into test dependencies

* Handling types with x*y multiplication.

* Avoid redundant typeof

* Added test for cumulative integration

* Precompute values outside of loop for RombergEven

Author: pengwyn

Project.toml

@@ -16,6 +16,7 @@ HCubature = "19dc6840-f33b-545b-b366-655c7e3ffd49"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [targets]
-test = ["Test", "InteractiveUtils", "StaticArrays", "HCubature"]
+test = ["Test", "InteractiveUtils", "StaticArrays", "HCubature", "Unitful"]

src/NumericalIntegration.jl

@@ -26,7 +26,6 @@ RombergEven() = RombergEven(1e-12)
 
 const HALF = 1//2
 
-
 #documentation
 
 """
@@ -147,12 +146,16 @@ end
 function integrate(x::AbstractVector, y::AbstractVector, m::RombergEven)
     @assert length(x) == length(y) "x and y vectors must be of the same length!"
     @assert ((length(x) - 1) & (length(x) - 2)) == 0 "Need length of vector to be 2^n + 1"
-    maxsteps = Int(log2(length(x)-1))
-    rombaux = zeros(eltype(y), maxsteps, 2)
+    maxsteps::Integer = Int(log2(length(x)-1))
+    @inbounds h = x[end] - x[1]
+    @inbounds v = (y[1] + y[end]) * h * HALF
+    rombaux = Matrix{typeof(v)}(undef, maxsteps, 2)
+    rombaux[1,1] = v
     prevcol = 1
     currcol = 2
-    @inbounds h = x[end] - x[1]
-    @inbounds rombaux[1, 1] = (y[1] + y[end])*h*HALF
+    # Precomputed values for norm check
+    acc_compare = m.acc * oneunit(norm(v,Inf))
+    minsteps = maxsteps ÷ 3
     @inbounds for i in 1 : (maxsteps-1)
         h *= HALF
         npoints = 1 << (i-1)
@@ -167,7 +170,9 @@ function integrate(x::AbstractVector, y::AbstractVector, m::RombergEven)
             rombaux[j, currcol] = (n_k*rombaux[j-1, currcol] - rombaux[j-1, prevcol])/(n_k - 1)
         end
 
-        if i > maxsteps//3 && norm(rombaux[i, prevcol] - rombaux[i+1, currcol], Inf) < m.acc
+        # Clunky way to address Unitful compatibility, while also allowing for vectors
+        normval = norm(rombaux[i, prevcol] - rombaux[i+1, currcol], Inf)
+        if i > minsteps && normval < acc_compare
             return rombaux[i+1, currcol]
         end
 

test/runtests.jl

@@ -77,3 +77,16 @@ end
     expwarn = "RombergEven :: final step reached, but accuracy not: 1.3333333333333335 > 1.0e-12"
     @test_logs (:warn, expwarn) integrate(xs, ys, m)
 end
+
+using Unitful
+
+@testset "Unitful compatibility" begin
+    x = LinRange(1u"s", 10u"s", 9)
+    y = rand(9)*u"m/s"
+    for M in subtypes(IntegrationMethod)
+        @test typeof(integrate(x,y,M())) == typeof(1.0u"m")
+        if hasmethod(cumul_integrate, Tuple{AbstractVector, AbstractVector, M})
+            @test typeof(cumul_integrate(x,y,M())) == Vector{typeof(1.0u"m")}
+        end
+    end
+end