Added function to perform cumulative integrals, and array inputs (#15)

* (1) Added function cumul_integrate to perform cumulative integrals (similar to functions like cumtrapz in Matlab)
(2) Added methods to compute integrals and cumulative integrals for each column of an input Array
(3) updated documentation of functions
(4) updated readme

* Modified name of function for cumulative integrals to "cumul_integrate". New methods to integrate along specific dimensions of a matrix. Use Trapzedoial() (instead of TrapezoidalFast()) as default method, for safety.

Author: nbrantut

README.md

@@ -16,22 +16,36 @@ Do note that while the code is trivial, it has not been extensively tested and d
 ```julia
 # setup some data
 x = collect(-π : π/1000 : π)
-y = sin(x)
+y = sin.(x)
 
 # integrate using the default TrapezoidalFast method
 integrate(x, y)
 
 # integrate using a specific method
 integrate(x, y, SimpsonEven())
+
+# compute cumulative integral
+Y = cumul_integrate(x, y)
+
+# compute cumulative integral for each column of an array
+z = [sin.(x) cos.(x) exp.(x/pi)]
+Z = cumul_integrate(x, z)
+
+# compute cumulative integral for each line of an array
+zp = permutedims(z) 
+ZP = cumul_integrate(x, zp, dims=1)
+
 ```
 
 The currently available methods are:
-- Trapezoidal
+- Trapezoidal (default)
 - TrapezoidalEven
-- TrapezoidalFast (default)
+- TrapezoidalFast
 - TrapezoidalEvenFast
 - SimpsonEven
 - SimpsonEvenFast
 - RombergEven
 
+Only Trapezoidal methods are available for cumulative integrals.
+
 All methods containing "Even" in the name assume evenly spaced data. All methods containing "Fast" omit basic correctness checks and focus on performance. Consequently, the fast methods will segfault or produce incorrect results if you supply incorrect data (vectors of different lengths, etc.). RombergEven needs a power of 2 + 1 points (so 9, 17, 33, 65, 129, 257, 513, 1025...) evenly spaced for it to work. Useful when control over accuracy is needed. 

src/NumericalIntegration.jl

@@ -3,7 +3,7 @@ module NumericalIntegration
 using LinearAlgebra
 using Logging
 
-export integrate
+export integrate, cumul_integrate
 export Trapezoidal, TrapezoidalEven, TrapezoidalFast, TrapezoidalEvenFast
 export SimpsonEven, SimpsonEvenFast
 export RombergEven
@@ -17,18 +17,46 @@ struct TrapezoidalFast     <: IntegrationMethod end
 struct TrapezoidalEvenFast <: IntegrationMethod end
 struct SimpsonEven         <: IntegrationMethod end # https://en.wikipedia.org/wiki/Simpson%27s_rule#Alternative_extended_Simpson.27s_rule
 struct SimpsonEvenFast     <: IntegrationMethod end
-struct RombergEven{T<:AbstractFloat}      <: IntegrationMethod
+struct RombergEven{T<:AbstractFloat} <: IntegrationMethod
     acc::T
 end # https://en.wikipedia.org/wiki/Romberg%27s_method
 RombergEven() = RombergEven(1e-12)
 
 const HALF = 1//2
 
+
+#documentation
+
+"""
+    integrate(x,y...)
+
+Compute numerical integral of y(x) from x=x[1] to x=x[end]. Return a scalar of the same type as the input. If not method is supplied, use TrapezdoialFast.
+"""
+function integrate(x,y...) end
+
+
+"""
+    cumul_integrate(x,y...)
+
+Compute cumulative numerical integral of y(x) from x=x[1] to x=x[end]. Return a vector with elements of the same type as the input. If not method is supplied, use TrapezdoialFast.
+"""
+function cumul_integrate(x,y...) end
+
+#implementation
+
 function _zero(x,y)
     ret = zero(eltype(x)) + zero(eltype(y))
-    ret / 2
 end
 
+function _zeros(x::AbstractVector,y::AbstractVector)
+    ret = zeros(eltype(x),size(x)) + zeros(eltype(y),size(y))
+end
+
+"""
+    integrate(x::AbstractVector, y::AbstractVector, ::Trapezoidal)
+
+Use Trapezoidal rule.
+"""
 function integrate(x::AbstractVector, y::AbstractVector, ::Trapezoidal)
     @assert length(x) == length(y) "x and y vectors must be of the same length!"
     retval = _zero(x,y)
@@ -38,11 +66,21 @@ function integrate(x::AbstractVector, y::AbstractVector, ::Trapezoidal)
     return HALF * retval
 end
 
+"""
+    integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalEven)
+
+Use Trapezoidal rule, assuming evenly spaced vector x.
+"""
 function integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalEven)
     @assert length(x) == length(y) "x and y vectors must be of the same length!"
     return (x[2] - x[1]) * (HALF * (y[1] + y[end]) + sum(y[2:end-1]))
 end
 
+"""
+    integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalFast)
+
+Use Trapezoidal rule. Unsafe method: no bound checking. This is the default when no method is supplied.
+"""
 function integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalFast)
     retval = _zero(x,y)
     @fastmath @simd for i in 1 : length(y)-1
@@ -51,6 +89,11 @@ function integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalFast)
     return HALF * retval
 end
 
+"""
+    integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalEvenFast)
+
+Use Trapezoidal rule, assuming evenly spaced vector x. Unsafe method: no bound checking.
+"""
 function integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalEvenFast)
     retval = _zero(x,y)
     N = length(y) - 1
@@ -60,6 +103,11 @@ function integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalEvenFast)
     @inbounds return (x[2] - x[1]) * (retval + HALF*y[1] + HALF*y[end])
 end
 
+"""
+    integrate(x::AbstractVector, y::AbstractVector, ::SimpsonEven)
+
+Use Simpson's rule, assuming evenly spaced vector x. 
+"""
 function integrate(x::AbstractVector, y::AbstractVector, ::SimpsonEven)
     @assert length(x) == length(y) "x and y vectors must be of the same length!"
     retval = (17*y[1] + 59*y[2] + 43*y[3] + 49*y[4] + 49*y[end-3] + 43*y[end-2] + 59*y[end-1] + 17*y[end]) / 48
@@ -69,6 +117,11 @@ function integrate(x::AbstractVector, y::AbstractVector, ::SimpsonEven)
     return (x[2] - x[1]) * retval
 end
 
+"""
+    integrate(x::AbstractVector, y::AbstractVector, ::SimpsonEven)
+
+Use Simpson's rule, assuming evenly spaced vector x.  Unsafe method: no bound checking.
+"""
 function integrate(x::AbstractVector, y::AbstractVector, ::SimpsonEvenFast)
     @inbounds retval = 17*y[1] + 59*y[2] + 43*y[3] + 49*y[4]
     @inbounds retval += 49*y[end-3] + 43*y[end-2] + 59*y[end-1] + 17*y[end]
@@ -79,6 +132,16 @@ function integrate(x::AbstractVector, y::AbstractVector, ::SimpsonEvenFast)
     @inbounds return (x[2] - x[1]) * retval
 end
 
+"""
+    integrate(x::AbstractVector, y::AbstractMatrix, method; dims=2)
+
+When y is an array, compute integral along dimension specified by dims (default 2: columns).
+"""
+function integrate(x::AbstractVector, y::AbstractMatrix, M::IntegrationMethod; dims=2)
+    out = [integrate(x,selectdim(y,dims,j),M) for j=1:size(y,dims)]
+    return out
+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"
@@ -113,6 +176,80 @@ function integrate(x::AbstractVector, y::AbstractVector, m::RombergEven)
     @inbounds return rombaux[maxsteps, prevrow]
 end
 
-integrate(x::AbstractVector, y::AbstractVector) = integrate(x, y, TrapezoidalFast())
 
+# cumulative integrals
+
+"""
+    cumul_integrate(x::AbstractVector, y::AbstractVector, ::Trapezoidal)
+
+Use Trapezoidal rule.
+"""
+function cumul_integrate(x::AbstractVector, y::AbstractVector, ::Trapezoidal)
+    @assert length(x) == length(y) "x and y vectors must be of the same length!"
+    retarr = _zeros(x,y)
+    for i in 2 : length(y)
+        retarr[i] = retarr[i-1] + (x[i] - x[i-1]) * (y[i] + y[i-1])
+    end
+    return HALF * retarr
 end
+
+"""
+    cumul_integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalEven)
+
+Use Trapezoidal rule, assuming evenly spaced vector x.
+"""
+function cumul_integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalEven)
+    @assert length(x) == length(y) "x and y vectors must be of the same length!"
+    retarr = _zeros(x,y)
+    for i in 2 : length(y)
+        retarr[i] = retarr[i-1] + (y[i-1] + y[i])
+    end
+    return (x[2] - x[1]) * HALF * retarr
+end
+
+"""
+    cumul_integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalFast)
+
+Use Trapezoidal rule. Unsafe method: no bound checking. This is the default when no method is supplied.
+"""
+function cumul_integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalFast)
+    retarr = _zeros(x,y)
+    @fastmath for i in 2 : length(y) #not sure if @simd can do anything here
+        @inbounds retarr[i] = retarr[i-1] + (x[i] - x[i-1]) * (y[i] + y[i-1])
+    end
+    return HALF * retarr
+end
+
+"""
+    cumul_integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalEvenFast)
+
+Use Trapezoidal rule, assuming evenly spaced vector x. Unsafe method: no bound checking.
+"""
+function cumul_integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalEvenFast)
+    retarr = _zeros(x,y)
+    @fastmath for i in 2 : length(y)
+        @inbounds retarr[i] = retarr[i-1] + (y[i] + y[i-1])
+    end
+    @inbounds return (x[2] - x[1]) * HALF * retarr
+end
+
+"""
+    cumul_integrate(x::AbstractVector, y::AbstractMatrix, method; dims=2)
+
+When y is an array, compute integral along each dimension specified by dims (default 2: columns)
+"""
+function cumul_integrate(x::AbstractVector, y::AbstractMatrix, M::IntegrationMethod; dims=2)
+    return hcat([cumul_integrate(x,selectdim(y,dims,j),M) for j=1:size(y,dims)]...)
+end
+
+
+#default behaviour
+integrate(x::AbstractVector, y::AbstractVector) = integrate(x, y, Trapezoidal())
+
+integrate(x::AbstractVector, y::AbstractMatrix; dims=2) = integrate(x, y, Trapezoidal(); dims=dims)
+
+cumul_integrate(x::AbstractVector, y::AbstractVector) = cumul_integrate(x, y, Trapezoidal())
+
+cumul_integrate(x::AbstractVector, y::AbstractMatrix; dims=2) = cumul_integrate(x, y, Trapezoidal(); dims=dims)
+
+end # module