Skip to content

Commit

Permalink
format
Browse files Browse the repository at this point in the history
  • Loading branch information
ArnoStrouwen committed Sep 23, 2023
1 parent 8da0a37 commit f4dbcc9
Show file tree
Hide file tree
Showing 6 changed files with 44 additions and 36 deletions.
13 changes: 7 additions & 6 deletions docs/src/basics/SampledIntegralProblem.md
Original file line number Diff line number Diff line change
@@ -1,13 +1,14 @@
# 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
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.
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
Expand All @@ -29,9 +30,9 @@ The exact aswer is of course \$ 1/3 \$.

### Non-equidistant grids

If the sampling points `x` are provided as an `AbstractRange`
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
the fact that the points are equidistantly spaced. Otherwise, general methods are used for
non-uniform grids.

Example:
Expand Down Expand Up @@ -70,4 +71,4 @@ Right now, only the `TrapezoidalRule` is supported, [see wikipedia](https://en.w

```@docs
TrapezoidalRule
```
```
6 changes: 5 additions & 1 deletion docs/src/tutorials/caching_interface.md
Original file line number Diff line number Diff line change
Expand Up @@ -55,6 +55,7 @@ If it is necessary to change the integrand `f` instead of defining a new

For sampled integral problems, it is possible to cache the weights and reuse
them for multiple data sets.

```@example cache2
using Integrals
Expand All @@ -72,14 +73,17 @@ sol1 = solve!(cache)
cache.y = cos.(x) # use .= to update in-place
sol2 = solve!(cache)
```

If the grid is modified, the weights are recomputed.

```@example cache2
cache.x = 0.0:0.2:2.0
cache.y = sin.(cache.x)
sol3 = solve!(cache)
```

For multi-dimensional datasets, the integration dimension can also be changed

```@example cache3
using Integrals
Expand All @@ -96,4 +100,4 @@ sol1 = solve!(cache)
```@example cache3
cache.dim = 1
sol2 = solve!(cache)
```
```
13 changes: 7 additions & 6 deletions src/common.jl
Original file line number Diff line number Diff line change
Expand Up @@ -95,7 +95,6 @@ function __solvebp_call(cache::IntegralCache, args...; kwargs...)
__solvebp_call(build_problem(cache), args...; kwargs...)
end


mutable struct SampledIntegralCache{Y, X, D, PK, A, K, Tc}
y::Y
x::X
Expand All @@ -117,13 +116,13 @@ end
function SciMLBase.init(prob::SampledIntegralProblem,
alg::SciMLBase.AbstractIntegralAlgorithm;
kwargs...)
NamedTuple(kwargs) == NamedTuple() || throw(ArgumentError("There are no keyword arguments allowed to `solve`"))
NamedTuple(kwargs) == NamedTuple() ||
throw(ArgumentError("There are no keyword arguments allowed to `solve`"))

cacheval = init_cacheval(alg, prob)
isfresh = true

SampledIntegralCache(
prob.y,
SampledIntegralCache(prob.y,
prob.x,
prob.dim,
prob.kwargs,
Expand All @@ -133,7 +132,6 @@ function SciMLBase.init(prob::SampledIntegralProblem,
cacheval)
end


"""
```julia
solve(prob::SampledIntegralProblem, alg::SciMLBase.AbstractIntegralAlgorithm; kwargs...)
Expand All @@ -154,5 +152,8 @@ function SciMLBase.solve!(cache::SampledIntegralCache)
end

function build_problem(cache::SampledIntegralCache)
SampledIntegralProblem(cache.y, cache.x; dim = dimension(cache.dim), cache.prob_kwargs...)
SampledIntegralProblem(cache.y,
cache.x;
dim = dimension(cache.dim),
cache.prob_kwargs...)
end
28 changes: 15 additions & 13 deletions src/sampled.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,30 +3,30 @@ 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)
@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), )
Base.size(w::UniformWeights) = (length(w),)

Check warning on line 9 in src/sampled.jl

View check run for this annotation

Codecov / codecov/patch

src/sampled.jl#L9

Added line #L9 was not covered by tests

# 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)
@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
Base.size(w::NonuniformWeights) = (length(w),)

Check warning on line 19 in src/sampled.jl

View check run for this annotation

Codecov / codecov/patch

src/sampled.jl#L19

Added line #L19 was not covered by tests

_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)
fw = zip(_eachslice(data, dims=dim), weights)
fw = zip(_eachslice(data, dims = dim), weights)
next = iterate(fw)
next === nothing && throw(ArgumentError("No points to integrate"))
(f1, w1), state = next
Expand All @@ -48,14 +48,16 @@ function evalrule(data::AbstractArray, weights, dim)
return out
end


# can be reused for other sampled rules, which should implement find_weights(x, alg)

function init_cacheval(alg::SciMLBase.AbstractIntegralAlgorithm, prob::SampledIntegralProblem)
function init_cacheval(alg::SciMLBase.AbstractIntegralAlgorithm,
prob::SampledIntegralProblem)
find_weights(prob.x, alg)
end

function __solvebp_call(cache::SampledIntegralCache, alg::SciMLBase.AbstractIntegralAlgorithm; kwargs...)
function __solvebp_call(cache::SampledIntegralCache,
alg::SciMLBase.AbstractIntegralAlgorithm;
kwargs...)
dim = dimension(cache.dim)
err = nothing
data = cache.y
Expand Down
13 changes: 7 additions & 6 deletions src/trapezoidal.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,18 +3,19 @@ struct TrapezoidalUniformWeights{T} <: UniformWeights
h::T
end

@inline Base.getindex(w::TrapezoidalUniformWeights, i) = ifelse((i == 1) || (i == w.n), w.h/2 , w.h)
@inline Base.getindex(w::TrapezoidalUniformWeights, i) = ifelse((i == 1) || (i == w.n),
w.h / 2,
w.h)


struct TrapezoidalNonuniformWeights{X<:AbstractArray} <: NonuniformWeights
struct TrapezoidalNonuniformWeights{X <: AbstractArray} <: NonuniformWeights
x::X
end

@inline function Base.getindex(w::TrapezoidalNonuniformWeights, i)
x = w.x
(i == firstindex(x)) && return (x[i + 1] - x[i])/2
(i == lastindex(x)) && return (x[i] - x[i - 1])/2
return (x[i + 1] - x[i - 1])/2
(i == firstindex(x)) && return (x[i + 1] - x[i]) / 2
(i == lastindex(x)) && return (x[i] - x[i - 1]) / 2
return (x[i + 1] - x[i - 1]) / 2
end

function find_weights(x::AbstractVector, ::TrapezoidalRule)
Expand Down
7 changes: 3 additions & 4 deletions test/sampled_tests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@ using Integrals, Test
npoints = 1000

grid1 = range(lb, ub, length = npoints)
grid2 = rand(npoints).*(ub-lb) .+ lb
grid2 = rand(npoints) .* (ub - lb) .+ lb
grid2 = [lb; sort(grid2); ub]

exact_sols = [1 / 6 * (ub^6 - lb^6), sin(ub) - sin(lb)]
Expand All @@ -21,7 +21,7 @@ using Integrals, Test

# along dim=2
y = f.([grid grid]')
prob = SampledIntegralProblem(y, grid; dim=2)
prob = SampledIntegralProblem(y, grid; dim = 2)
error = solve(prob, method()).u .- exact
@test all(error .< 10^-4)
end
Expand All @@ -30,7 +30,6 @@ using Integrals, Test
end

@testset "Caching interface" begin

x = 0.0:0.1:1.0
y = sin.(x)

Expand Down Expand Up @@ -67,5 +66,5 @@ end
cache.dim = 1
sol2 = solve!(cache)

@test sol2 == solve(SampledIntegralProblem(y, x, dim=1), alg)
@test sol2 == solve(SampledIntegralProblem(y, x, dim = 1), alg)
end

0 comments on commit f4dbcc9

Please sign in to comment.