Skip to content

Commit

Permalink
Add simplemaxima and friends (#42)
Browse files Browse the repository at this point in the history 
* Add fast but restrictive `simplemaxima` and `simpleminima` functions

* Include simplemaxima/minima docstrings

* Limit methoderror hint to when called with eltype <: Missing

* add tests

* Dedup simplemaxima/simpleminima function decls
  • Loading branch information
@halleysfifthinc
halleysfifthinc authored Jul 31, 2024
1 parent d8e52ce commit 76debfc
Show file tree
Hide file tree
Showing 4 changed files with 171 additions and 4 deletions.
2 changes: 2 additions & 0 deletions docs/src/reference.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,8 @@ maxima
minima
findmaxima
findminima
simplemaxima
simpleminima
```

## Peak characteristics & filtering
Expand Down
26 changes: 23 additions & 3 deletions src/Peaks.jl
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
module Peaks

export argmaxima, argminima, maxima, minima, findmaxima, findminima, findnextmaxima,
findnextminima, peakproms, peakproms!, peakwidths, peakwidths!, peakheights,
peakheights!, ismaxima, isminima, isplateau, filterpeaks!
export argmaxima, simplemaxima, argminima, simpleminima, maxima, minima, findmaxima,
findminima, findnextmaxima, findnextminima, peakproms, peakproms!, peakwidths,
peakwidths!, peakheights, peakheights!, ismaxima, isminima, isplateau, filterpeaks!

include("minmax.jl")
include("utils.jl")
Expand All @@ -11,4 +11,24 @@ include("peakwidth.jl")
include("peakheight.jl")
include("plot.jl")

function __init__()
Base.Experimental.register_error_hint(MethodError) do io, exc, argtypes, kwargs
if (exc.f == simplemaxima || exc.f == simpleminima) && eltype(argtypes[1]) <: Missing
if exc.f == simplemaxima
cmp = "max"
else
cmp = "min"
end
printstyled(io, "\nsimple$(cmp)ima"; color=:cyan)
print(io, " does not support vectors with ")
printstyled(io, "missing"; color=:cyan)
print(io, "s. See the ")
printstyled(io, "arg$(cmp)ima"; color=:cyan)
print(io, " function to find $(cmp)ima when ")
printstyled(io, "missing"; color=:cyan)
print(io, "s are present.")
end
end
end

end # module Peaks
112 changes: 112 additions & 0 deletions src/minmax.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -85,6 +85,34 @@ function findnextextrema(cmp, x::AbstractVector, i::Int, w::Int, strict::Bool)
return lastindex(x) + 1
end

function _simpleextrema(f, cmp::F, x::AbstractVector{T}) where {F,T}
T >: Missing && throw(MethodError(simplemaxima, Tuple{typeof(x)}))

pks = Int[]
i = firstindex(x) + 1
@inbounds while i < lastindex(x)
xi = x[i]
pre = cmp(x[i-1], xi)
post = cmp(x[i+1], xi)
plat = x[i+1] === xi

if plat
j = something(findnext(Base.Fix2(!==, xi), x, i+2), lastindex(x)+1)
post = j > lastindex(x) ? false : cmp(x[j], xi)
end

if pre && post
push!(pks, i)
i = plat ? j+1 : i+2
else
i += 1
end
end

return pks
end


"""
findnextmaxima(x, i[, w=1; strict=true]) -> Int
Expand Down Expand Up @@ -177,6 +205,48 @@ function argmaxima(
return pks
end

"""
simplemaxima(x) -> Vector{Int}
Find the indices of local maxima in `x`, where each maxima `i` is greater than both adjacent
elements or is the first index of a plateau.
A plateau is defined as a maxima with consecutive equal (`===`/egal) maximal values which
are bounded by lesser values immediately before and after the plateau.
This function is semantically equivalent to `argmaxima(x, w=1; strict=true)`, but is
faster because of its simplified set of features. (The difference in speed scales with
`length(x)`, up to ~2-3x faster.)
Vectors with `missing`s are not supported by `simplemaxima`, use `argmaxima` if this is
needed.
See also: [`argmaxima`](@ref)
# Examples
```julia-repl
julia> simplemaxima([0,2,0,1,1,0])
2-element Vector{Int64}:
2
4
julia> argmaxima([0,2,0,1,1,0])
2-element Vector{Int64}:
2
4
julia> simplemaxima([2,0,1,1])
Int64[]
julia> @btime simplemaxima(x) setup=(x = repeat([0,1]; outer=100));
620.759 ns (4 allocations: 1.92 KiB)
julia> @btime argmaxima(x) setup=(x = repeat([0,1]; outer=100));
1.079 μs (4 allocations: 1.92 KiB)
```
"""
simplemaxima(x::AbstractVector) = _simpleextrema(simplemaxima, <, x)

"""
maxima(x[, w=1; strict=true]) -> Vector{eltype(x)}
Expand Down Expand Up @@ -316,6 +386,48 @@ function argminima(
return pks
end

"""
simpleminima(x) -> Vector{Int}
Find the indices of local minima in `x`, where each minima `i` is less than both adjacent
elements or is the first index of a plateau.
A plateau is defined as a minima with consecutive equal (`===`/egal) minimal values which
are bounded by greater values immediately before and after the plateau.
This function is semantically equivalent to `argminima(x, w=1; strict=true)`, but is faster
because of its simplified set of features. (The difference in speed scales with `length(x)`,
up to ~2-3x faster.)
Vectors with `missing`s are not supported by `simpleminima`, use `argminima` if this is
needed.
See also: [`argminima`](@ref)
# Examples
```julia-repl
julia> simpleminima([3,2,3,1,1,3])
2-element Vector{Int64}:
2
4
julia> argminima([3,2,3,1,1,3])
2-element Vector{Int64}:
2
4
julia> simpleminima([2,3,1,1])
Int64[]
julia> @btime simpleminima(x) setup=(x = repeat([0,1]; outer=100));
671.388 ns (4 allocations: 1.92 KiB)
julia> @btime argminima(x) setup=(x = repeat([0,1]; outer=100));
1.175 μs (4 allocations: 1.92 KiB)
```
"""
simpleminima(x::AbstractVector) = _simpleextrema(simpleminima, >, x)

"""
minima(x[, w=1; strict=true]) -> Vector{eltype(x)}
Expand Down
35 changes: 34 additions & 1 deletion test/minmax.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@ t = T:T:fs

x1 = a*sin.(2*pi*f1*T*t)+b*sin.(2*pi*f2*T*t)+c*sin.(2*pi*f3*T*t);

@testset "argminima/argmaxima" begin
@testset "maxima/minima" begin
@test_throws DomainError Peaks.findnextextrema(<, x1, 1, 0, false)
@test length(argmaxima(x1)) == 30
@test length(argminima(x1)) == 30
Expand All @@ -29,6 +29,9 @@ x1 = a*sin.(2*pi*f1*T*t)+b*sin.(2*pi*f2*T*t)+c*sin.(2*pi*f3*T*t);
@test argmaxima([0,1,0,1,0,1,0,1,0,1,0,1,0,1,0], 2) == [4,8,12]
@test argmaxima([0,1,0,1,0,1,0,1,0,1,0,1,0,1,0], 2; strict=false) == [2,6,10,14]

@test simplemaxima(x1) == argmaxima(x1)
@test simpleminima(x1) == argminima(x1)

# OffsetArrays
x1_off = OffsetArray(x1, -200:length(x1)-201)
@test length(argmaxima(x1_off)) == 30
Expand All @@ -39,6 +42,8 @@ x1 = a*sin.(2*pi*f1*T*t)+b*sin.(2*pi*f2*T*t)+c*sin.(2*pi*f3*T*t);
@test length(argminima(x1_off, 1; strict=false)) == 31
@test argmaxima(x1_off, 1; strict=false) == argmaxima(x1, 1; strict=false) .- 201
@test argminima(x1_off, 1; strict=false) == argminima(x1, 1; strict=false) .- 201
@test simplemaxima(x1_off) == argmaxima(x1_off)
@test simpleminima(x1_off) == argminima(x1_off)

@testset "Plateaus" begin
p1 = [0,1,1,1,1,1]
Expand All @@ -59,6 +64,9 @@ x1 = a*sin.(2*pi*f1*T*t)+b*sin.(2*pi*f2*T*t)+c*sin.(2*pi*f3*T*t);
@test isempty(argmaxima( p3))
@test isempty(argminima(-p3))

@test isempty(simplemaxima( p1))
@test isempty(simpleminima(-p1))

@test argmaxima(reverse( p1), 2; strict=false) == [1]
@test argminima(reverse(-p1), 2; strict=false) == [1]
@test argmaxima(reverse( p2), 2; strict=false) == [1]
Expand All @@ -73,13 +81,21 @@ x1 = a*sin.(2*pi*f1*T*t)+b*sin.(2*pi*f2*T*t)+c*sin.(2*pi*f3*T*t);
@test isempty(argmaxima(reverse( p3)))
@test isempty(argminima(reverse(-p3)))

@test isempty(simplemaxima(reverse( p1)))
@test isempty(simpleminima(reverse(-p1)))

@test argmaxima( [0,1,1,1,1,0]) == [2]
@test argminima(-[0,1,1,1,1,0]) == [2]
@test argmaxima( [0,1,1,1,2,0]) == [5]
@test argminima(-[0,1,1,1,2,0]) == [5]
@test argmaxima( [0,1,1,0,2,1], 3; strict=false) == [5]
@test argminima(-[0,1,1,0,2,1], 3; strict=false) == [5]

@test simplemaxima( [0,1,1,1,1,0]) == argmaxima( [0,1,1,1,1,0])
@test simpleminima(-[0,1,1,1,1,0]) == argminima(-[0,1,1,1,1,0])
@test simplemaxima( [0,1,1,1,2,0]) == argmaxima( [0,1,1,1,2,0])
@test simpleminima(-[0,1,1,1,2,0]) == argminima(-[0,1,1,1,2,0])

# issue #4
@test isempty(argmaxima(zeros(10)))
@test argmaxima(zeros(10), 1; strict=false) == [1]
Expand All @@ -100,6 +116,7 @@ x1 = a*sin.(2*pi*f1*T*t)+b*sin.(2*pi*f2*T*t)+c*sin.(2*pi*f3*T*t);
@testset "Missings and NaNs" begin
m1 = [0,0,1,1,1,missing,missing,0,0]
n1 = [0,0,1,1,1,NaN,NaN,0,0]

@test isempty(argmaxima( m1, 1; strict=true))
@test isempty(argmaxima( n1, 1; strict=true))
@test isempty(argminima(-m1, 1; strict=true))
Expand All @@ -109,6 +126,11 @@ x1 = a*sin.(2*pi*f1*T*t)+b*sin.(2*pi*f2*T*t)+c*sin.(2*pi*f3*T*t);
@test argminima(-m1, 1; strict=false) == [3,8]
@test argminima(-n1, 1; strict=false) == [3,8]

@test_throws MethodError simplemaxima(m1)
@test_throws MethodError simpleminima(m1)
@test isempty(simplemaxima( n1))
@test isempty(simpleminima(-n1))

@test isempty(argmaxima(reverse( m1), 1; strict=true))
@test isempty(argmaxima(reverse( n1), 1; strict=true))
@test isempty(argminima(reverse(-m1), 1; strict=true))
Expand All @@ -118,6 +140,9 @@ x1 = a*sin.(2*pi*f1*T*t)+b*sin.(2*pi*f2*T*t)+c*sin.(2*pi*f3*T*t);
@test argminima(reverse(-m1), 1; strict=false) == [1,5]
@test argminima(reverse(-n1), 1; strict=false) == [1,5]

@test isempty(simplemaxima(reverse( n1)))
@test isempty(simpleminima(reverse(-n1)))

m2 = [0,1,2,1,missing,missing,missing]
n2 = [0,1,2,1,NaN,NaN,NaN]
@test argmaxima(m2, 1) == [3]
Expand All @@ -137,13 +162,21 @@ x1 = a*sin.(2*pi*f1*T*t)+b*sin.(2*pi*f2*T*t)+c*sin.(2*pi*f3*T*t);
@test isempty(argminima(reverse(-m2), 2))
@test isempty(argminima(reverse(-n2), 2))

@test simplemaxima(n2) == argmaxima(n2, 1)
@test simpleminima(-n2) == argminima(-n2, 1)
@test simplemaxima(reverse(n2)) == argmaxima(reverse(n2), 1)
@test simpleminima(reverse(-n2)) == argminima(reverse(-n2), 1)

m3 = [0,1,0,1,missing,1,0,1,0]
n3 = [0,1,0,1,NaN,1,0,1,0]
@test argmaxima(m3) == [2,8]
@test argminima(-m3) == [2,8]
@test argmaxima(n3) == [2,8]
@test argminima(-n3) == [2,8]

@test simpleminima(n3) == argminima(n3)
@test simpleminima(-n3) == argminima(-n3)

mn = [1,NaN,missing,1]
@test isempty(argmaxima(mn))
@test isempty(argminima(mn))
Expand Down

0 comments on commit 76debfc

Please sign in to comment.