Skip to content

Commit

Permalink
integrate via iterator
Browse files Browse the repository at this point in the history
  • Loading branch information
@nluetts
nluetts committed Jan 9, 2024
1 parent 596eb05 commit e119469
Show file tree
Hide file tree
Showing 4 changed files with 112 additions and 316 deletions.
121 changes: 10 additions & 111 deletions src/common.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,13 +24,13 @@ Find peak in interval `p ± w/2` and return (speak position - `w`/2, peak positi
"""
function left_right_from_peak(x, y, p, w)
# find indices of interval [left, right]
left = p - w/2
right = p + w/2
left = p - w / 2
right = p + w / 2
l = searchsortedfirst(x, left) # left index: x[l] <= p - w
r = searchsortedlast(x, right) # right index: x[r] >= p + w

m = l + argmax(view(y, l:r)) - 1 # index of local maximum at x[m]
return [x[m] - w/2, x[m] + w/2]
return [x[m] - w / 2, x[m] + w / 2]
end


Expand All @@ -41,118 +41,17 @@ end
"""Linearly interpolate y-value at position x between two points (x₀, y₀) and (x₁, y₁)."""
@inline function lininterp(x, x₀, x₁, y₀, y₁)
δx = x₁ - x₀
y = (y₁*(x - x₀) + y₀*(x₁ - x)) / δx
y = (y₁ * (x - x₀) + y₀ * (x₁ - x)) / δx
return y
end

"""Linearly interpolate y-value at position `x` for x,y data; `xs` has to be sorted."""
@inline function lininterp(x::T, xs::Vector{T}, ys::Vector{S}) where {T <: Number, S <: Number}
@inline function lininterp(x::T, xs::Vector{T}, ys::Vector{S}) where {T<:Number,S<:Number}
i = searchsortedfirst(xs, x)
(i < 2 || i > length(xs)) && throw(error("x is outside the support: $x ∉ ($(minimum(xs)), $(maximum(xs))]."))
return lininterp(x, xs[i-1], xs[i], ys[i-1], ys[i])
end

"""
trapz(x::AbstractArray{T}, y::AbstractArray{T}, left, right) where {T<:AbstractFloat}
Integrate vector `y` in interval [`left`, `right`] using trapezoidal integration.
# Notes
`left` and `right` must support conversion to the datatype `T`.
If `left` and `right` do not fall on the `x`-grid, additional data points will be interpolated linearly.
(i.e. the width of the first and last trapezoid will be somewhat smaller).
If `left` and/or `right` falls outside the `x`-range, the integration window will be cropped
to the available range.
# Examples
```jldoctest
julia> x = collect(Float64, 0:10);
julia> y = ones(Float64, 11);
julia> trapz(x, y, 1, 3)
2.0
julia> trapz(x, y, -1, 11) # at most, we can integrate the available x-range, 0 to 10
10.0
julia> trapz(x, y, -10, 20)
10.0
julia> trapz(x, y, 1.1, 1.3) ≈ 0.2 # if we integrate "between" the grid, data points are interpolated
true
```
"""
function trapz(x::AbstractArray{T}, y::AbstractArray{T}, left::T, right::T) where {T<:AbstractFloat}

left, right = left < right ? (left, right) : (right, left)

A = zero(eltype(y)) # Area to be returned

N = length(x)
N != length(y) && throw(ArgumentError("Data arrays `x` and `y` must have the same length."))
N < 2 && return A

yl = yr = nothing # the y-values of the left and right integration bound, to be interpolated
lastiter = false
j = 2

while j <= N

x₀ = x[j-1]
x₁ = x[j]
y₀ = y[j-1]
y₁ = y[j]

if x₁ <= left
j += 1
continue
elseif yl === nothing
# this will only run once, when we enter the integration window
# test whether x₀ should be replaced by `left`
if x₀ < left
y₀ = lininterp(left, x₀, x₁, y₀, y₁)
x₀ = left
else
# this case means that `left` <= x[1]
left = x₀
end
yl = y₀
end

# test whether x₁ should be replaced by `right`
if x₁ >= right
# we move out of the integration window
yr = x₁ == right ? y₁ : lininterp(right, x₀, x₁, y₀, y₁)
x₁ = right
y₁ = yr
lastiter = true # we shall break the loop after this iteration
end

A += singletrapz(x₀, x₁, y₀, y₁)

lastiter && break

j += 1
end

return A
end

function trapz(x::AbstractArray{T}, y::AbstractArray{T}, left, right) where {T<:AbstractFloat}
left = T(left)
right = T(right)
return trapz(x, y, left, right)
end


"""
@samples n::Integer e::Expression
Expand All @@ -177,9 +76,9 @@ MonteCarloMeasurements.Particles{Float64, 9999}
macro samples(n::Integer, e::Expr)
err = ArgumentError("Expression not understood.")
if length(e.args) != 3
return :( throw($err) )
return :(throw($err))
end

op, x, y = e.args

if op == :±
Expand All @@ -195,6 +94,6 @@ macro samples(n::Integer, e::Expr)
Particles($n, Uniform(a, b))
end
else
return :( throw($err) )
return :(throw($err))
end
end
end
1 change: 0 additions & 1 deletion src/curves.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -370,7 +370,6 @@ end
function Base.iterate(
iter::Tuple{Curve{Float64},Float64,Float64},
)::Union{Nothing,Tuple{Tuple{Float64,Float64},Tuple{Int64,Bool}}}

mv, left, right = iter
left, right = min(left, right), max(left, right)
if left == right || mv.x[end] <= left || mv.x[1] >= right
Expand Down
66 changes: 62 additions & 4 deletions src/integration.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,64 @@ end

_endpoint_to_endpoint_baseline(xs, ys, xₗ, xᵣ) = (xₗ, xᵣ, lininterp(xₗ, xs, ys), lininterp(xᵣ, xs, ys))

function trapz(crv::Curve{T}, left::T, right::T) where {T<:AbstractFloat}
left, right = min(left, right), max(left, right)
area = zero(T)
for ((xi, yi), (xj, yj)) in zip((crv, left, right), Iterators.drop((crv, left, right), 1))
area += singletrapz(xi, xj, yi, yj)
end
return area
end

"""
trapz(x::AbstractArray{T}, y::AbstractArray{T}, left, right) where {T<:AbstractFloat}
Integrate vector `y` in interval [`left`, `right`] using trapezoidal integration.
# Notes
`left` and `right` must support conversion to the datatype `T`.
If `left` and `right` do not fall on the `x`-grid, additional data points will be interpolated linearly.
(i.e. the width of the first and last trapezoid will be somewhat smaller).
If `left` and/or `right` falls outside the `x`-range, the integration window will be cropped
to the available range.
# Examples
```jldoctest
julia> x = collect(Float64, 0:10);
julia> y = ones(Float64, 11);
julia> trapz(x, y, 1, 3)
2.0
julia> trapz(x, y, -1, 11) # at most, we can integrate the available x-range, 0 to 10
10.0
julia> trapz(x, y, -10, 20)
10.0
julia> trapz(x, y, 1.1, 1.3) ≈ 0.2 # if we integrate "between" the grid, data points are interpolated
true
```
"""
function trapz(xs::AbstractArray{T}, ys::AbstractArray{T}, left::T, right::T) where {T<:AbstractFloat}
trapz(Curve(xs, ys), left, right)
end

function trapz(x::AbstractArray{T}, y::AbstractArray{T}, left, right) where {T<:AbstractFloat}
left = T(left)
right = T(right)
return trapz(x, y, left, right)
end


"""
mc_integrate(uc::UncertainCurve{T, N}, bnds::Vector{UncertainBound{T, M}}; intfun=trapz)
Expand All @@ -36,7 +94,7 @@ If true, for each draw a local linear baseline defined by the integration window
The y-values of the start and end point are derived from a weighted average over the start and end point distributions, see
[the documentation](https://nluetts.github.io/NoisySignalIntegration.jl/dev/baseline/#Build-in) for further information.
"""
function mc_integrate(uc::UncertainCurve{T, N}, bnds::Vector{UncertainBound{T, M}}; intfun=trapz, subtract_baseline=false, local_baseline=false) where {T, M, N}
function mc_integrate(uc::UncertainCurve{T,N}, bnds::Vector{UncertainBound{T,M}}; intfun=trapz, subtract_baseline=false, local_baseline=false) where {T,M,N}

M != N && error("Samples sizes incompatible")
subtract_baseline && @warn("subtract_baseline keyword argument is deprecated, use local_baseline instead.")
Expand All @@ -58,9 +116,9 @@ function mc_integrate(uc::UncertainCurve{T, N}, bnds::Vector{UncertainBound{T, M
end
end
end
return [Particles(areas[:,i]) for i in 1:size(areas)[2]]
return [Particles(areas[:, i]) for i in 1:size(areas)[2]]
end

function mc_integrate(uc::S, bnd::T; intfun=trapz, subtract_baseline=false, local_baseline=false) where {S <: UncertainCurve, T <: UncertainBound}
function mc_integrate(uc::S, bnd::T; intfun=trapz, subtract_baseline=false, local_baseline=false) where {S<:UncertainCurve,T<:UncertainBound}
return mc_integrate(uc, [bnd]; intfun=intfun, subtract_baseline=subtract_baseline, local_baseline=local_baseline)[1]
end
end
Loading

0 comments on commit e119469

Please sign in to comment.