From aa677c4df3f38045c5fea3a4bb2ff6517aa0b531 Mon Sep 17 00:00:00 2001
From: IlianPihlajamaa <73794090+IlianPihlajamaa@users.noreply.github.com>
Date: Tue, 19 Sep 2023 20:57:15 +0200
Subject: [PATCH] add trapezoidal rule for sampled data, attempt 2

---
 docs/pages.jl                             |  1 +
 docs/src/basics/SampledIntegralProblem.md | 73 +++++++++++++++++++++++
 src/Integrals.jl                          |  1 +
 src/algorithms.jl                         | 15 ++++-
 src/common.jl                             |  3 -
 src/sampled.jl                            | 57 ++++++++++++++++++
 src/trapezoidal.jl                        | 70 ++++++----------------
 test/sampled_tests.jl                     |  2 +
 8 files changed, 167 insertions(+), 55 deletions(-)
 create mode 100644 docs/src/basics/SampledIntegralProblem.md
 create mode 100644 src/sampled.jl

diff --git a/docs/pages.jl b/docs/pages.jl
index 76c058fc..99ef1719 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -2,6 +2,7 @@ pages = ["index.md",
     "Tutorials" => Any["tutorials/numerical_integrals.md",
         "tutorials/differentiating_integrals.md"],
     "Basics" => Any["basics/IntegralProblem.md",
+        "basics/SampledIntegralProblem.md",
         "basics/solve.md",
         "basics/FAQ.md"],
     "Solvers" => Any["solvers/IntegralSolvers.md"],
diff --git a/docs/src/basics/SampledIntegralProblem.md b/docs/src/basics/SampledIntegralProblem.md
new file mode 100644
index 00000000..02317cd4
--- /dev/null
+++ b/docs/src/basics/SampledIntegralProblem.md
@@ -0,0 +1,73 @@
+# Integrating pre-sampled data
+
+In some cases, instead of a function that acts as integrand, 
+one only possesses a list of data points `y` at a set of sampling 
+locations `x`, that must be integrated. This package contains functionality
+for doing that. 
+
+## Example
+
+Say, by some means we have generated a dataset `x` and `y`:
+```example 1
+using Integrals # hide
+f = x -> x^2
+x = range(0, 1, length=20)
+y = f.(x)
+```
+
+Now, we can integrate this data set as follows:
+
+```example 1
+problem = SampledIntegralProblem(y, x)
+method = TrapezoidalRule()
+solve(problem, method)
+```
+
+The exact aswer is of course \$ 1/3 \$.
+
+## Details
+
+### Non-equidistant grids
+
+If the sampling points `x` are provided as an `AbstractRange` 
+(constructed with the `range` function for example), faster methods are used that take advantage of
+the fact that the points are equidistantly spaced. Otherwise, general methods are used for 
+non-uniform grids.
+
+Example:
+
+```example 2
+using Integrals # hide
+f = x -> x^7
+x = [0.0; sort(rand(1000)); 1.0]
+y = f.(x)
+problem = SampledIntegralProblem(y, x)
+method = TrapezoidalRule()
+solve(problem, method)
+```
+
+### Evaluating multiple integrals at once
+
+If the provided data set `y` is a multidimensional array, the integrals are evaluated across only one
+of its axes. For performance reasons, the last axis of the array `y` is chosen by default, but this can be modified with the `dim`
+keyword argument to the problem definition.
+
+```example 3
+using Integrals # hide
+f1 = x -> x^2
+f2 = x -> x^3
+f3 = x -> x^4
+x = range(0, 1, length=20)
+y = [f1.(x) f2.(x) f3.(x)]
+problem = SampledIntegralProblem(y, x; dim=1)
+method = TrapezoidalRule()
+solve(problem, method)
+```
+
+### Supported methods
+
+Right now, only the `TrapezoidalRule` is supported, [see wikipedia](https://en.wikipedia.org/wiki/Trapezoidal_rule).
+
+```@docs
+TrapezoidalRule
+```
\ No newline at end of file
diff --git a/src/Integrals.jl b/src/Integrals.jl
index f407ae53..02b5bc5d 100644
--- a/src/Integrals.jl
+++ b/src/Integrals.jl
@@ -13,6 +13,7 @@ include("init.jl")
 include("algorithms.jl")
 include("infinity_handling.jl")
 include("quadrules.jl")
+include("sampled.jl")
 include("trapezoidal.jl")
 
 abstract type QuadSensitivityAlg end
diff --git a/src/algorithms.jl b/src/algorithms.jl
index 37f474d0..0adefdbd 100644
--- a/src/algorithms.jl
+++ b/src/algorithms.jl
@@ -127,10 +127,23 @@ end
     TrapezoidalRule
 
 Struct for evaluating an integral via the trapezoidal rule.
+
+
+Example with sampled data:
+
+```
+using Integrals 
+f = x -> x^2
+x = range(0, 1, length=20)
+y = f.(x)
+problem = SampledIntegralProblem(y, x)
+method = TrapezoidalRul()
+solve(problem, method)
+```
 """
 struct TrapezoidalRule <: SciMLBase.AbstractIntegralAlgorithm
 end
-  
+
 """
   QuadratureRule(q; n=250)
 
diff --git a/src/common.jl b/src/common.jl
index ea41f152..39876049 100644
--- a/src/common.jl
+++ b/src/common.jl
@@ -101,6 +101,3 @@ function __solvebp_call(cache::IntegralCache, args...; kwargs...)
     __solvebp_call(build_problem(cache), args...; kwargs...)
 end
 
-@inline _selectdim(y::AbstractArray{T, dims}, d, i) where {T, dims} = selectdim(y, d, i)
-@inline _selectdim(y::AbstractArray{T, 1}, _, i) where {T} = @inbounds y[i]
-@inline dimension(::Val{D}) where {D} = D
diff --git a/src/sampled.jl b/src/sampled.jl
new file mode 100644
index 00000000..9cdcbe21
--- /dev/null
+++ b/src/sampled.jl
@@ -0,0 +1,57 @@
+abstract type AbstractWeights end
+
+# must have field `n` for length, and a field `h` for stepsize
+abstract type UniformWeights <: AbstractWeights end
+@inline Base.iterate(w::UniformWeights) = (0 == w.n) ? nothing : (w[1], 1)
+@inline Base.iterate(w::UniformWeights, i) = (i == w.n) ? nothing : (w[i+1], i+1)
+Base.length(w::UniformWeights) = w.n
+Base.eltype(w::UniformWeights) = typeof(w.h)
+Base.size(w::UniformWeights) = (length(w), )
+
+# must contain field `x` which are the sampling points
+abstract type NonuniformWeights <: AbstractWeights end 
+@inline Base.iterate(w::NonuniformWeights) = (0 == length(w.x)) ? nothing : (w[firstindex(w.x)], firstindex(w.x))
+@inline Base.iterate(w::NonuniformWeights, i) = (i == lastindex(w.x)) ? nothing : (w[i+1], i+1)
+Base.length(w::NonuniformWeights) = length(w.x)
+Base.eltype(w::NonuniformWeights) = eltype(w.x)
+Base.size(w::NonuniformWeights) = (length(w), )
+
+_eachslice(data::AbstractArray; dims=ndims(data)) = eachslice(data; dims=dims)
+_eachslice(data::AbstractArray{T, 1}; dims=ndims(data)) where T = data
+
+
+# these can be removed when the Val(dim) is removed from SciMLBase
+dimension(::Val{D}) where {D} = D
+dimension(D::Int) = D 
+
+
+function evalrule(data::AbstractArray, weights, dim)
+    f = _eachslice(data, dims=dim)
+    firstidx, lastidx = firstindex(f), lastindex(f)
+    out = f[firstidx]*weights[firstidx]
+    if isbits(out)
+        for i in firstidx+1:lastidx
+            @inbounds out += f[i]*weights[i]
+        end
+    else
+        for i in firstidx+1:lastidx
+            @inbounds out .+= f[i] .* weights[i]
+        end
+    end
+    return out
+
+end
+
+
+# can be reused for other sampled rules
+function __solvebp_call(prob::SampledIntegralProblem, alg::TrapezoidalRule; kwargs...)
+    dim = dimension(prob.dim)
+    err = nothing
+    data = prob.y
+    grid = prob.x
+    weights = find_weights(grid, alg)
+    I = evalrule(data, weights, dim)
+    return SciMLBase.build_solution(prob, alg, I, err, retcode = ReturnCode.Success)
+end
+
+
diff --git a/src/trapezoidal.jl b/src/trapezoidal.jl
index a58a4fb6..1056fca3 100644
--- a/src/trapezoidal.jl
+++ b/src/trapezoidal.jl
@@ -1,55 +1,23 @@
-function __solvebp_call(prob::SampledIntegralProblem, alg::TrapezoidalRule; kwargs...)
-    dim = dimension(prob.dim)
-    err = Inf64
-    data = prob.y
-    grid = prob.x
-    # inlining is required in order to not allocate
-    integrand = @inline function (i)
-        # integrate along dimension `dim`, returning a n-1 dimensional array, or scalar if n=1
-        _selectdim(data, dim, i)
-    end
+struct TrapezoidalUniformWeights <: UniformWeights
+    n::Int
+    h::Float64
+end
+
+@inline Base.getindex(w::TrapezoidalUniformWeights, i) = ifelse((i == 1) || (i == w.n), w.h*0.5 , w.h)
 
-    firstidx, lastidx = firstindex(grid), lastindex(grid)
 
-    out = integrand(firstidx)
+struct TrapezoidalNonuniformWeights{X<:AbstractArray} <: NonuniformWeights
+    x::X
+end
 
-    if isbits(out)
-        # fast path for equidistant grids
-        if grid isa AbstractRange
-            dx = step(grid)
-            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)
+@inline function Base.getindex(w::TrapezoidalNonuniformWeights, i) 
+    x = w.x
+    (i == firstindex(x)) && return (x[i + 1] - x[i])*0.5
+    (i == lastindex(x)) && return (x[i] - x[i - 1])*0.5
+    return (x[i + 1] - x[i - 1])*0.5
 end
+
+function find_weights(x::AbstractVector, ::TrapezoidalRule)
+    x isa AbstractRange && return TrapezoidalUniformWeights(length(x), step(x))
+    return TrapezoidalNonuniformWeights(x)
+end
\ No newline at end of file
diff --git a/test/sampled_tests.jl b/test/sampled_tests.jl
index 03d838cf..c394f6a8 100644
--- a/test/sampled_tests.jl
+++ b/test/sampled_tests.jl
@@ -16,12 +16,14 @@ using Integrals, Test
                 # single dimensional y
                 y = f.(grid)
                 prob = SampledIntegralProblem(y, grid)
+                @show solve(prob, method()).u, exact
                 error = solve(prob, method()).u .- exact
                 @test all(error .< 10^-4)
 
                 # along dim=2
                 y = f.([grid grid]')
                 prob = SampledIntegralProblem(y, grid; dim=2)
+                @show solve(prob, method()).u, exact
                 error = solve(prob, method()).u .- exact
                 @test all(error .< 10^-4)
             end