From f4dbcc9830029f4b137cfc857fa4878c1a696c47 Mon Sep 17 00:00:00 2001 From: Arno Strouwen Date: Sat, 23 Sep 2023 02:37:22 +0200 Subject: [PATCH] format --- docs/src/basics/SampledIntegralProblem.md | 13 ++++++----- docs/src/tutorials/caching_interface.md | 6 ++++- src/common.jl | 13 ++++++----- src/sampled.jl | 28 ++++++++++++----------- src/trapezoidal.jl | 13 ++++++----- test/sampled_tests.jl | 7 +++--- 6 files changed, 44 insertions(+), 36 deletions(-) diff --git a/docs/src/basics/SampledIntegralProblem.md b/docs/src/basics/SampledIntegralProblem.md index 02317cd4..a51e2ec3 100644 --- a/docs/src/basics/SampledIntegralProblem.md +++ b/docs/src/basics/SampledIntegralProblem.md @@ -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 @@ -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: @@ -70,4 +71,4 @@ Right now, only the `TrapezoidalRule` is supported, [see wikipedia](https://en.w ```@docs TrapezoidalRule -``` \ No newline at end of file +``` diff --git a/docs/src/tutorials/caching_interface.md b/docs/src/tutorials/caching_interface.md index 02020524..61e5f5c4 100644 --- a/docs/src/tutorials/caching_interface.md +++ b/docs/src/tutorials/caching_interface.md @@ -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 @@ -72,7 +73,9 @@ 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) @@ -80,6 +83,7 @@ sol3 = solve!(cache) ``` For multi-dimensional datasets, the integration dimension can also be changed + ```@example cache3 using Integrals @@ -96,4 +100,4 @@ sol1 = solve!(cache) ```@example cache3 cache.dim = 1 sol2 = solve!(cache) -``` \ No newline at end of file +``` diff --git a/src/common.jl b/src/common.jl index e00ae2c9..01e6ca6a 100644 --- a/src/common.jl +++ b/src/common.jl @@ -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 @@ -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, @@ -133,7 +132,6 @@ function SciMLBase.init(prob::SampledIntegralProblem, cacheval) end - """ ```julia solve(prob::SampledIntegralProblem, alg::SciMLBase.AbstractIntegralAlgorithm; kwargs...) @@ -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 diff --git a/src/sampled.jl b/src/sampled.jl index 16453b54..947ad962 100644 --- a/src/sampled.jl +++ b/src/sampled.jl @@ -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),) # 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),) +_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 @@ -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 diff --git a/src/trapezoidal.jl b/src/trapezoidal.jl index ff6777de..496a8af8 100644 --- a/src/trapezoidal.jl +++ b/src/trapezoidal.jl @@ -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) diff --git a/test/sampled_tests.jl b/test/sampled_tests.jl index ef60ded8..8a77c63f 100644 --- a/test/sampled_tests.jl +++ b/test/sampled_tests.jl @@ -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)] @@ -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 @@ -30,7 +30,6 @@ using Integrals, Test end @testset "Caching interface" begin - x = 0.0:0.1:1.0 y = sin.(x) @@ -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