add Trapezoidal rule #173

IlianPihlajamaa · 2023-09-07T11:43:32Z

!! Not yet ready to merge !!
EDIT: I think it's ready now!

To do:

implement rule itself
add tests
add documentation
put the type piracy into SciMLBase

I have added a trapezoidal rule for integrating sampled data. Let me know what you think. If/once you agree with the code, I will write docs (and add Simpson's as well in the same way). Right now, I only have a docstring for Trapezoidal.

Additionally, I have some duplicate code for in- and out-of-place operations, that I don't really know how to combine. See lines 224-241 of Integrals.jl, do you have an idea how to streamline such code?

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:

using Integrals
f = (x, p) -> x^9
n = 1000
method = Trapezoidal(n)
problem = IntegralProblem(f, 0.0, 1.0)
solve(problem, method) # u: 0.10000075150155016

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.

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) # u: 0.10000490350275731

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.

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) 
#u: 2-element Vector{Float64}:
# 0.2000015680246171
# 0.10000365747228422

In order to make integration over arrays work, I needed to add a method (because isinplace(f::AbstractArray) is not defined:

import SciMLBase.IntegralProblem
function IntegralProblem(y::AbstractArray, lb, ub, args...; kwargs...)
    IntegralProblem{false}(y, lb, ub, args...; kwargs...)
end

I think this belongs in SciMLBase. I will make a PR there too if you agree this is the right approach.

ChrisRackauckas · 2023-09-09T20:06:32Z

src/Integrals.jl

+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


Yeah move this over to SciMLBase

ChrisRackauckas · 2023-09-09T20:07:15Z

src/Integrals.jl

+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


This implementation of Trapezoid should be in a separate file trapezoid.jl that implements the method, and is then included into the top level file.

ChrisRackauckas · 2023-09-09T20:08:21Z

src/Integrals.jl

+    if is_sampled_problem(prob)
+        data = prob.f
+        # inlining is required in order to not allocate
+        @inline function integrand(i) 


Suggested change

@inline function integrand(i)

integrand = function (i)

It should be an anonymous function if in an if statement

ChrisRackauckas · 2023-09-09T20:08:59Z

Just small organizational things, but otherwise looks great!

ChrisRackauckas · 2023-09-09T20:09:19Z

In order to make integration over arrays work, I needed to add a method (because isinplace(f::AbstractArray) is not defined:

Yes, let's get that merged and bumped first, then bump the lower bound and build off of that.

ChrisRackauckas · 2023-09-09T20:10:09Z

Though I think more data is actually required, since it may need to know the x at which the value f(x) is associated. So it needs two arrays. Might as well just make that a DataIntegralProblem or something at that point to keep it clean.

sathvikbhagavan · 2023-09-11T15:17:28Z

Is it also possible to implement a method to return a vector of integrated values - i.e. from beginning to that point - like in https://github.com/dextorious/NumericalIntegration.jl/blob/master/src/NumericalIntegration.jl#L44?

IlianPihlajamaa · 2023-09-11T18:34:58Z

Is it also possible to implement a method to return a vector of integrated values - i.e. from beginning to that point - like in https://github.com/dextorious/NumericalIntegration.jl/blob/master/src/NumericalIntegration.jl#L44?

Yeah, that is possible of course... We can define a CumulativeTrapezoidal() quadrature to do this. Let me implement the standard methods first though.

IlianPihlajamaa · 2023-09-11T18:40:39Z

Though I think more data is actually required, since it may need to know the x at which the value f(x) is associated. So it needs two arrays. Might as well just make that a DataIntegralProblem or something at that point to keep it clean.

Ah so you prefer the x-data to be inside the problem, not the method. I guess that makes sense.

So the API should be:

x = 0.0:0.1:1.0
y = x.^2
method = Trapezoidal()
problem = DataIntegralProblem(x, y; dim=1) # dim optional
solve(problem, method)

f = x -> x^2 
method = Trapezoidal()
problem = IntegralProblem(f, 0.0, 1.0)
solve(problem, method; x = 0.0:0.1:1.0)

Do you agree?

ChrisRackauckas · 2023-09-11T20:23:13Z

Ah so you prefer the x-data to be inside the problem, not the method. I guess that makes sense.

Yes, it's part of the problem definition, so it should be defined in solve either. So not solve(problem, method; x = 0.0:0.1:1.0), that x should move into the problem definition somehow?

ChrisRackauckas · 2023-09-11T20:24:28Z

To be clear, all solve dispatches across all of SciML are expected to have the same keyword arguments, so adding an x one for this specific case wouldn't happen. The keyword arguments are meant to be general solver controls.

IlianPihlajamaa · 2023-09-13T07:26:30Z

Ah so you prefer the x-data to be inside the problem, not the method. I guess that makes sense.

Yes, it's part of the problem definition, so it should be defined in solve either. So not solve(problem, method; x = 0.0:0.1:1.0), that x should move into the problem definition somehow?

Alright, so that means we need to put an x=nothing field into the IntegralProblem type. This also makes the DataIntegralProblem type redundant.

I will change this in the PR to SciMLBase.

EDIT: just did. Will continue this implementation when that is merged. (I dont know how to develop multiple interdependent packages at the same time).

lxvm · 2023-09-16T17:55:58Z

It sounds like this PR is suggesting 2 things:

Adding a DataIntegralProblem for sampled data (related pr here)
Adding trapezoidal integration on regular grids for functions (e.g. the first example in the OP)

Point 1 is novel and should be addressed in the pr. Point 2 could be addressed in a separate pr because it is a generalization of what is already implemented in the FastGaussQuadrature.jl extension. An api for this could be

struct QuadratureFunction{Q} <: IntegralAlgorithm
    q::Q
    n::Int
end

which accepts a function, q, of the form x, w = q(n) which given a number of points n returns vectors x and w of length n that can be used to compute the quadrature dot(w, f.(x)). It would also be possible to cache the nodes and weights so they could be used across multiple intervals [lb, ub]. I implemented something similar to this for one of my packages here so I could open a pr for this.

IlianPihlajamaa · 2023-09-16T18:27:50Z

Regarding point 2, I agree that integrating a function with the trapezoidal (or any other) quadrature rule is a different problem to that of integrating data. However i do think we should at least try to get the api for the former to be consistent with the current one of this package, so something of the form

solve(IntegralProblem(f, lb, ub), Trapezoidal(n))

The Trapezoidal struct could be responsible for the optional caching of the points and weights. Do you agree?

ChrisRackauckas · 2023-09-16T18:43:59Z

I implemented something similar to this for one of my packages here so I could open a pr for this.

Yes that would be nice to have.

Regarding point 2, I agree that integrating a function with the trapezoidal (or any other) quadrature rule is a different problem to that of integrating data. However i do think we should at least try to get the api for the former to be consistent with the current one of this package, so something of the form

Yes

lxvm · 2023-09-16T18:52:56Z

So given the QuadratureFunction I proposed above, I would implement Trapezoidal as

"""
    trapz(n::Integer)

Return the weights and nodes on the standard interval [-1,1] of the [trapezoidal
rule](https://en.wikipedia.org/wiki/Trapezoidal_rule).
"""
function trapz(n::Integer)
    @assert n > 1
    r = range(-1, 1, length=n)
    x = collect(r)
    halfh = step(r)/2
    h = step(r)
    w = [ (i == 1) || (i == n) ? halfh : h for i in 1:n ]
    return (x, w)
end

struct Trapezoidal <: IntegralAlgorithm end # use this constructor for DataIntegralProblems
Trapezoidal(n) = QuadratureFunction(trapz, n) # use this constructor for IntegralProblems

and you could implement GaussLegendre(n) = QuadratureFunction(gausslegendre, n). Then all of these different quadrature rules could be implemented in the Integrals.jl api because they are all basically computing the same thing: sum(w .* f.(x)).

The Trapezoidal struct or function doesn't need to store the weights because the SciML init interface creates a cache that we can use. However, this only applies to point 2 and if you want to apply the Trapezoidal rule to nonuniform data, then the grid should be passed to a DataIntegralProblem.

lxvm · 2023-09-16T18:56:09Z

I'll follow up with a QuadratureFunction pr today

ChrisRackauckas · 2023-09-16T18:56:54Z

BTW, the name Trapezoidal is already used in OrdinaryDiffEq so we should probably prevent the name clash.

lxvm · 2023-09-16T18:59:09Z

How about TrapezoidalRule? I think the 'rule' suffix helps clarify that this is just a quadrature rule, not a an algorithm that converges to a requested tolerance

ChrisRackauckas · 2023-09-16T19:00:30Z

👍 I like that suggestion.

I'll follow up with a QuadratureFunction pr today

Awesome! It's great to see this finally getting some love again.

IlianPihlajamaa · 2023-09-16T19:13:31Z

So given the QuadratureFunction I proposed above, I would implement Trapezoidal as

and you could implement GaussLegendre(n) = QuadratureFunction(gausslegendre, n). Then all of these different quadrature rules could be implemented in the Integrals.jl api because they are all basically computing the same thing: sum(w .* f.(x)).

The Trapezoidal struct or function doesn't need to store the weights because the SciML init interface creates a cache that we can use. However, this only applies to point 2 and if you want to apply the Trapezoidal rule to nonuniform data, then the grid should be passed to a DataIntegralProblem.

This looks great, and i agree this would be a nice way to implement many rules. In the simple cases though, it may be possible to avoid allocating. eg, x does not need to be collected, and w could be a vector-like struct that has a getindex/dot method.

lxvm · 2023-09-16T19:50:52Z

Yes, it is not necessary for x or w to be allocated. Probably iteration and a getindex method is enough. We also should also keep this compatible with batched integrands, so the exact requirements for x and w may depend on the implementation. This will have to be documented, although the point of the cache is that you only need to allocate once, which will happen at init time and not solve! time. Starting with plain Vectors should cover most cases even though it may come with a memory penalty for n very large.

IlianPihlajamaa · 2023-09-19T18:59:19Z

Sorry for the delay, here is a rewrite of the original code. If approved I can straightforwardly implement also Simpson's rule. The trapezoidal rule for integrating functions can use a combination of this code, and the QuadratureRule implementation.

Let me know how it can be improved.

Apparently the `eachslice behaviour changed in 1.9 causing tests to fail...

lxvm · 2023-09-20T19:18:01Z

We should also implement the init interface for SampledIntegralProblem since this pr adds it, so I'll add that

lxvm · 2023-09-20T20:46:26Z

@IlianPihlajamaa I'm not sure how to commit to your pr, so see my fork for how to add the init interface. It will have a type instability until the dim field of SampledIntegralProblem is made an Int

lxvm · 2023-09-20T20:58:20Z

I realized I would have to make a PR to your branch

IlianPihlajamaa · 2023-09-21T14:02:14Z

@IlianPihlajamaa I'm not sure how to commit to your pr, so see my fork for how to add the init interface. It will have a type instability until the dim field of SampledIntegralProblem is made an Int

Looks great! If you make a PR to my branch I will merge it! I don't know an easier way to include your work. Do you?

lxvm · 2023-09-21T14:47:10Z

Looks great! If you make a PR to my branch I will merge it! I don't know an easier way to include your work. Do you?

Yes, see https://github.com/IlianPihlajamaa/Integrals.jl/pull/1

add init interface for SampledIntegralProblem

ChrisRackauckas · 2023-09-21T15:14:42Z

docs/src/tutorials/caching_interface.md

+sol3 = solve!(cache)
+```
+
+For multi-dimensional datasets, the integration dimension can also be changed


oh that's pretty cool.

ChrisRackauckas · 2023-09-21T16:16:09Z

src/common.jl

+function SciMLBase.init(prob::SampledIntegralProblem,
+    alg::SciMLBase.AbstractIntegralAlgorithm;
+    kwargs...)
+    NamedTuple(kwargs) == NamedTuple() || throw(ArgumentError("There are no keyword arguments allowed to `solve`"))


that's odd, why not? There are many keyword arguments that would go here?

@lxvm is there a specific reason why not, or did you just not expect it to ever be necessary?

The current implementation of the solver doesn't use any keyword arguments, so I didn't want an api with keywords. If future algorithms for sampled integral problems need them I would expect this to change, but there are no convergence criteria for this kind of problem, so I wasn't expecting any

Oh just the sampled data methods. Okay, yeah for now might as well throw. We always try to throw if keyword arguments that are incorrect so this is good.

src/trapezoidal.jl

ChrisRackauckas · 2023-09-21T16:27:37Z

Just a few comments, other than that I think it's pretty ready to go.

put the type piracy into SciMLBase

Where is this? I didn't see a piracy.

IlianPihlajamaa · 2023-09-21T16:37:45Z

Just a few comments, other than that I think it's pretty ready to go.

put the type piracy into SciMLBase

Where is this? I didn't see a piracy.

It's already in SciMLBase

Co-authored-by: Christopher Rackauckas <[email protected]>

ChrisRackauckas · 2023-09-21T16:39:16Z

oh it wasn't checked 😅

Co-authored-by: Christopher Rackauckas <[email protected]>

IlianPihlajamaa · 2023-09-21T16:44:20Z

oh it wasn't checked 😅

Yes, my bad :)

lxvm · 2023-09-21T17:28:06Z

@IlianPihlajamaa What does this PR do about multidimensional grids? We can certainly define a multidimensional trapezoidal rule on regular grids as a product of 1d rules, but what would the API be? Would the user pass in a x array constructed like an Iterators.product? Or would x be a vector of SVectors?

If we don't want to support this then should we check x isa AbstractVector{<:Number}?

IlianPihlajamaa · 2023-09-21T17:47:41Z

@IlianPihlajamaa What does this PR do about multidimensional grids? We can certainly define a multidimensional trapezoidal rule on regular grids as a product of 1d rules, but what would the API be? Would the user pass in a x array constructed like an Iterators.product? Or would x be a vector of SVectors?

If we don't want to support this then should we check x isa AbstractVector{<:Number}?

The constructor in SciMLBase already guarantees that x isa AbstractVector. Because multidimensional integration with a product rule is so easy to do by just repeating single-dimensional integration, I havent implemented it explicitly. If we want to add it anyway for the convenience or added performance, we should indeed think of an API (pethaps in a separate issue, because it is not trivial if we also want the option to support non-product grids or other domain shapes).

lxvm · 2023-09-21T18:21:01Z

Thanks, that sounds good. I agree that multidimensional domains should be handled separately, and the user would probably have to supply some extra information about the domain. Congrats on the very nice PR!

IlianPihlajamaa and others added 6 commits September 7, 2023 13:23

add Trapezoidal rule + tests

32e438e

fix correctness bug

881804b

small change for consistency

afc0f12

remove issorted check, and return scalar instead of 0-dim arr

ade9f6b

small efficiency improvement

7265f24

clarify error msgs

c4a9208

ChrisRackauckas reviewed Sep 9, 2023

View reviewed changes

IlianPihlajamaa mentioned this pull request Sep 11, 2023

Add DataIntegralProblem SciML/SciMLBase.jl#491

Merged

implement trapezoidal rule for sampled data

efaf2e1

add trapezoidal rule for sampled data, attempt 2

aa677c4

now also works for Julia<= 1.9

e7aefc9

Apparently the `eachslice behaviour changed in 1.9 causing tests to fail...

add init interface for SampledIntegralProblem

608af76

lxvm added 3 commits September 20, 2023 18:05

bump sciml version

2eb0090

use zip for evalrule iterator

346c0ac

add docs and tests

9ab8ab2

Merge pull request #1 from lxvm/IlianPihlajamaa/master

6ba2e3c

add init interface for SampledIntegralProblem

ChrisRackauckas reviewed Sep 21, 2023

View reviewed changes

src/trapezoidal.jl Outdated Show resolved Hide resolved

ChrisRackauckas reviewed Sep 21, 2023

View reviewed changes

src/trapezoidal.jl Outdated Show resolved Hide resolved

ChrisRackauckas reviewed Sep 21, 2023

View reviewed changes

src/trapezoidal.jl Outdated Show resolved Hide resolved

Change type parameter in TrapezoidalUniformWeights

813a460

Co-authored-by: Christopher Rackauckas <[email protected]>

IlianPihlajamaa and others added 2 commits September 21, 2023 18:40

Improve type safety

7c3d710

Co-authored-by: Christopher Rackauckas <[email protected]>

Improve type safety

65cf88d

Co-authored-by: Christopher Rackauckas <[email protected]>

ChrisRackauckas approved these changes Sep 21, 2023

View reviewed changes

ChrisRackauckas merged commit c917535 into SciML:master Sep 21, 2023
6 of 7 checks passed

add Trapezoidal rule #173

add Trapezoidal rule #173

Conversation

IlianPihlajamaa commented Sep 7, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ChrisRackauckas commented Sep 9, 2023

ChrisRackauckas commented Sep 9, 2023

ChrisRackauckas commented Sep 9, 2023

sathvikbhagavan commented Sep 11, 2023

IlianPihlajamaa commented Sep 11, 2023

IlianPihlajamaa commented Sep 11, 2023

ChrisRackauckas commented Sep 11, 2023

ChrisRackauckas commented Sep 11, 2023

IlianPihlajamaa commented Sep 13, 2023 • edited Loading

lxvm commented Sep 16, 2023 • edited Loading

IlianPihlajamaa commented Sep 16, 2023

ChrisRackauckas commented Sep 16, 2023

lxvm commented Sep 16, 2023

lxvm commented Sep 16, 2023

ChrisRackauckas commented Sep 16, 2023

lxvm commented Sep 16, 2023

ChrisRackauckas commented Sep 16, 2023

IlianPihlajamaa commented Sep 16, 2023

lxvm commented Sep 16, 2023

IlianPihlajamaa commented Sep 19, 2023 • edited Loading

lxvm commented Sep 20, 2023

lxvm commented Sep 20, 2023

lxvm commented Sep 20, 2023

IlianPihlajamaa commented Sep 21, 2023

lxvm commented Sep 21, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ChrisRackauckas commented Sep 21, 2023

IlianPihlajamaa commented Sep 21, 2023

ChrisRackauckas commented Sep 21, 2023

IlianPihlajamaa commented Sep 21, 2023

lxvm commented Sep 21, 2023

IlianPihlajamaa commented Sep 21, 2023

lxvm commented Sep 21, 2023

IlianPihlajamaa commented Sep 7, 2023 •

edited

Loading

IlianPihlajamaa commented Sep 13, 2023 •

edited

Loading

lxvm commented Sep 16, 2023 •

edited

Loading

IlianPihlajamaa commented Sep 19, 2023 •

edited

Loading