add Trapezoidal rule

Project.toml

@@ -13,10 +13,22 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 
+[weakdeps]
+ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
+
+[extensions]
+IntegralsFastGaussQuadratureExt = "FastGaussQuadrature"
+IntegralsForwardDiffExt = "ForwardDiff"
+IntegralsZygoteExt = ["Zygote", "ChainRulesCore"]
+
 [compat]
 ChainRulesCore = "0.10.7, 1"
 CommonSolve = "0.2"
 Distributions = "0.23, 0.24, 0.25"
+FastGaussQuadrature = "0.5"
 ForwardDiff = "0.10"
 HCubature = "1.4"
 MonteCarloIntegration = "0.0.1, 0.0.2, 0.0.3"
@@ -26,16 +38,11 @@ Requires = "1"
 SciMLBase = "1.70"
 Zygote = "0.4.22, 0.5, 0.6"
 julia = "1.6"
-FastGaussQuadrature = "0.5"
-
-[extensions]
-IntegralsForwardDiffExt = "ForwardDiff"
-IntegralsZygoteExt = ["Zygote", "ChainRulesCore"]
-IntegralsFastGaussQuadratureExt = "FastGaussQuadrature"
 
 [extras]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
@@ -44,13 +51,6 @@ SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
-FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 
 [targets]
 test = ["SciMLSensitivity", "StaticArrays", "FiniteDiff", "Pkg", "SafeTestsets", "Test", "Distributions", "ForwardDiff", "Zygote", "ChainRulesCore", "FastGaussQuadrature"]
-
-[weakdeps]
-ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
-Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
-FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"

src/Integrals.jl

@@ -147,5 +147,107 @@
     SciMLBase.build_solution(prob, alg, val, err, chi = chi, retcode = ReturnCode.Success)
 end
 
-export QuadGKJL, HCubatureJL, VEGAS, GaussLegendre
+is_sampled_problem(prob::IntegralProblem) = prob.f isa AbstractArray
+import SciMLBase.IntegralProblem # this is type piracy, and belongs in SciMLBase
+function IntegralProblem(y::AbstractArray, lb, ub, args...; kwargs...)
+    IntegralProblem{false}(y, lb, ub, args...; kwargs...)
+end
+function construct_grid(prob, alg, lb, ub, dim)
+    x = alg.spec
+    @assert length(ub) == length(lb) == 1 "Multidimensional integration is not supported with the Trapezoidal method"
+    if x isa Integer
+        grid = range(lb[1], ub[1], length=x)
+    else 
+        grid = x
+        @assert ndims(grid) == 1 "Multidimensional integration is not supported with the Trapezoidal method"
+    end
+
+    @assert lb[1] ≈ grid[begin] "Lower bound in `IntegralProblem` must coincide with that of the grid"
+    @assert ub[1] ≈ grid[end] "Upper bound in `IntegralProblem` must coincide with that of the grid"
+    if is_sampled_problem(prob)
+        @assert size(prob.f, dim) == length(grid) "Integrand and grid must be of equal length along the integrated dimension"
+        @assert axes(prob.f, dim) == axes(grid,1) "Grid and integrand array must use same indexing along integrated dimension" 
+    end
+    return grid
+end
+
+@inline myselectdim(y::AbstractArray{T,dims}, d, i) where {T,dims} = selectdim(y, d, i)
+@inline myselectdim(y::AbstractArray{T,1}, _, i) where {T} = @inbounds y[i]
+
+@inline dimension(::Val{D}) where D = D
+function __solvebp_call(prob::IntegralProblem, alg::Trapezoidal{S, D}, sensealg, lb, ub, p; kwargs...) where {S,D}
+    # since all AbstractRange types are equidistant by design, we can rely on that
+    @assert prob.batch == 0
+    # using `Val`s for dimensionality is required to make `selectdim` not allocate
+    dim = dimension(D) 
+    p = p
+    if is_sampled_problem(prob)
+        @assert alg.spec isa AbstractArray "For pre-sampled problems where the integrand is an array, the integration grid must also be specified by an array."
+    end
+
+    grid = construct_grid(prob, alg, lb, ub, dim)
+
+    err = Inf64
+    if is_sampled_problem(prob)
+        data = prob.f
+        # inlining is required in order to not allocate
+        @inline function integrand(i) 
+            # integrate along dimension `dim`, returning a n-1 dimensional array, or scalar if n=1
+            myselectdim(data, dim, i) 
+        end 
+    else
+        if isinplace(prob)
+            y = zeros(eltype(lb), prob.nout)
+            integrand = i -> @inbounds (prob.f(y, grid[i], p); y)
+        else
+            integrand = i -> @inbounds prob.f(grid[i], p)
+        end
+    end
+
+    firstidx, lastidx = firstindex(grid), lastindex(grid)
+
+    out = integrand(firstidx)
+
+    if isbits(out) 
+        # fast path for equidistant grids
+        if grid isa AbstractRange 
+            dx = grid[begin+1] - grid[begin]
+            out /= 2
+            for i in (firstidx+1):(lastidx-1)
+                out += integrand(i)
+            end
+            out += integrand(lastidx)/2
+            out *= dx
+        # irregular grids:
+        else 
+            out *= (grid[firstidx + 1] - grid[firstidx])
+            for i in (firstidx+1):(lastidx-1)
+                @inbounds out += integrand(i) * (grid[i + 1] - grid[i-1])
+            end
+            out += integrand(lastidx) * (grid[lastidx] - grid[lastidx-1])
+            out /= 2
+        end
+    else # same, but inplace, broadcasted
+        out = copy(out) # to prevent aliasing
+        if grid isa AbstractRange 
+            dx = grid[begin+1] - grid[begin]
+            out ./= 2
+            for i in (firstidx+1):(lastidx-1)
+                out .+= integrand(i)
+            end
+            out .+= integrand(lastidx) ./ 2
+            out .*= dx
+        else 
+            out .*= (grid[firstidx + 1] - grid[firstidx])
+            for i in (firstidx+1):(lastidx-1)
+                @inbounds out .+= integrand(i) .* (grid[i + 1] - grid[i-1])
+            end
+            out .+= integrand(lastidx) .* (grid[lastidx] - grid[lastidx-1])
+            out ./= 2
+        end
+    end
+    return SciMLBase.build_solution(prob, alg, out, err, retcode = ReturnCode.Success)
+end
+
+export QuadGKJL, HCubatureJL, VEGAS, GaussLegendre, Trapezoidal
 end # module

src/algorithms.jl

@@ -122,3 +122,66 @@
     end
     return GaussLegendre(nodes, weights, subintervals)
 end
+
+"""
+    Trapezoidal{S, DIM}
+
+Struct for evaluating an integral via Trapezoidal rule.
+The field `spec` contains either the number of gridpoints or an array of specified gridpoints
+
+The Trapezoidal rule supports integration of pre-sampled data, stored in an array, as well as integration of 
+functions. It does not support batching or integration over multidimensional spaces.
+
+To use the Trapezoidal rule to integrate a function on a regular grid with `n` points:
+
+```@example trapz1
+using Integrals
+f = (x, p) -> x^9
+n = 1000
+method = Trapezoidal(n)
+problem = IntegralProblem(f, 0.0, 1.0)
+solve(problem, method)
+```
+
+To use the Trapezoidal rule to integrate a function on an predefined irregular grid, see the following example.
+Note that the lower and upper bound of integration must coincide with the first and last element of the grid. 
+
+```@example trapz2
+using Integrals
+f = (x, p) -> x^9
+x = sort(rand(1000))
+x = [0.0; x; 1.0]
+method = Trapezoidal(x)
+problem = IntegralProblem(f, 0.0, 1.0)
+solve(problem, method)
+```
+
+To use the Trapezoidal rule to integrate a set of sampled data, see the following example.
+By default, the integration occurs over the first dimension of the input array.
+```@example trapz3
+using Integrals
+x = sort(rand(1000))
+x = [0.0; x; 1.0]
+y1 = x' .^ 4
+y2 = x' .^ 9
+y = [y1; y2]
+method = Trapezoidal(x; dim=2)
+problem = IntegralProblem(y, 0.0, 1.0)
+solve(problem, method)
+```
+"""
+struct Trapezoidal{S, DIM} <: SciMLBase.AbstractIntegralAlgorithm
+    spec::S
+    function Trapezoidal(npoints::I; dim=1) where I<:Integer
+        @assert npoints > 1 
+        return new{I, Val(dim)}(npoints)
+    end
+    function Trapezoidal(grid::V; dim=1) where V<:AbstractVector
+        npoints = length(grid)
+        @assert npoints > 1
+        @assert isfinite(first(grid))
+        @assert isfinite(last(grid))
+        return new{V, Val(dim)}(grid)
+    end
+end
+

test/interface_tests.jl

@@ -8,7 +8,7 @@ reltol = 1e-3
 abstol = 1e-3
 
 algs = [QuadGKJL(), HCubatureJL(), CubatureJLh(), CubatureJLp(), #VEGAS(), #CubaVegas(),
-    CubaSUAVE(), CubaDivonne(), CubaCuhre()]
+    CubaSUAVE(), CubaDivonne(), CubaCuhre(), Trapezoidal(1000)]
 
 alg_req = Dict(QuadGKJL() => (nout = 1, allows_batch = false, min_dim = 1, max_dim = 1,
         allows_iip = false),
@@ -27,7 +27,9 @@ alg_req = Dict(QuadGKJL() => (nout = 1, allows_batch = false, min_dim = 1, max_d
     CubaDivonne() => (nout = Inf, allows_batch = true, min_dim = 2,
         max_dim = Inf, allows_iip = true),
     CubaCuhre() => (nout = Inf, allows_batch = true, min_dim = 2, max_dim = Inf,
-        allows_iip = true))
+        allows_iip = true),
+    Trapezoidal(1000) => (nout = Inf, allows_batch = false, min_dim = 1, max_dim = 1,
+    allows_iip = true))
 
 integrands = [
     (x, p) -> 1.0,

test/runtests.jl

@@ -22,3 +22,7 @@ end
 @time @safetestset "Gaussian Quadrature Tests" begin
     include("gaussian_quadrature_tests.jl")
 end
+
+@time @safetestset "Sampled Integration Tests" begin
+    include("sampled_tests.jl")
+end

test/sampled_tests.jl

@@ -0,0 +1,45 @@
+using Integrals, Test
+@testset "Sampled Integration" begin
+    lb = 0.0
+    ub = 1.0
+    npoints = 1000
+    for method = [Trapezoidal] # Simpson's later
+        nouts = [1,2,1,2]
+        for (i,f) = enumerate([(x,p) -> x^5, (x,p) -> [x^5, x^5], (out, x,p) -> (out[1] = x^5; out),  (out, x, p) -> (out[1] = x^5; out[2] = x^5; out)])
+
+            exact = 1/6
+            prob = IntegralProblem(f, lb, ub, nout=nouts[i])
+
+            # AbstractRange
+            error1 = solve(prob, method(npoints)).u .- exact
+            @test all(error1 .< 10^-4)
+
+            # AbstractVector equidistant
+            error2 = solve(prob, method(collect(range(lb, ub, length=npoints)))).u .- exact
+            @test all(error2 .< 10^-4)
+
+            # AbstractVector irregular
+            grid = rand(npoints)
+            grid = [lb; sort(grid); ub]
+            error3 = solve(prob, method(grid)).u .- exact
+            @test all(error3 .< 10^-4)
+
+
+        end
+        exact = 1/6
+
+        grid = rand(npoints)
+        grid = [lb; sort(grid); ub]
+        # single dimensional y
+        y = grid .^ 5
+        prob = IntegralProblem(y, lb, ub)
+        error4 = solve(prob, method(grid, dim=1)).u .- exact
+        @test all(error4 .< 10^-4) 
+
+        # along dim=2
+        y = ([grid grid]') .^ 5
+        prob = IntegralProblem(y, lb, ub)
+        error5 = solve(prob, method(grid, dim=2)).u .- exact 
+        @test all(error5 .< 10^-4) 
+    end
+end