diff --git a/docs/src/reference.md b/docs/src/reference.md index cc39028..3c3c4b6 100644 --- a/docs/src/reference.md +++ b/docs/src/reference.md @@ -7,6 +7,8 @@ maxima minima findmaxima findminima +simplemaxima +simpleminima ``` ## Peak characteristics & filtering diff --git a/src/Peaks.jl b/src/Peaks.jl index df8ca2a..be72746 100644 --- a/src/Peaks.jl +++ b/src/Peaks.jl @@ -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") @@ -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 diff --git a/src/minmax.jl b/src/minmax.jl index fb2abcd..3a51cdf 100644 --- a/src/minmax.jl +++ b/src/minmax.jl @@ -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 @@ -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)} @@ -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)} diff --git a/test/minmax.jl b/test/minmax.jl index dbecac2..8bf2498 100644 --- a/test/minmax.jl +++ b/test/minmax.jl @@ -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 @@ -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 @@ -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] @@ -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] @@ -73,6 +81,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(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] @@ -80,6 +91,11 @@ 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,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] @@ -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)) @@ -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)) @@ -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] @@ -137,6 +162,11 @@ 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] @@ -144,6 +174,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(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))