findpeaks from pracma

[laborg]

hi,

I've implemented findpeaks as found in the pracma R package (https://github.com/cran/pracma). The original author gave his OK for a transcription into Julia, now the question is, if you are open to accept this function as a PR. (I don't want to create/register/maintain a package just for a single function, and your package Peaks.jl looks like a good candidate to include such a function).

[halleysfifthinc]

That is generally the intent of this package, to serve as a single source of peak finding algorithms. However, I don't think that a straight port/replication of the MATLAB findpeaks function (which as I understand it, is the implementation of the findpeaks function in the pracma package) matches the typical Julian API style.

Also, the existing maxima/minima and peakprom functions in Peaks.jl already cover much of the basic functionality*. Is there a particular feature from the MATLAB/R functions that you need that isn't already implemented here?

*Peak width is the only non-trivial functionality currently missing

[laborg]

After a couple of hours playing around I still couldn't replicate the exact pracma findpeaks call with Peaks.jl (nor with FindPeaks.jl), so thats my motivation. I agree that the signature of pracmas findpeaks is not at all what I would expect from a Julia function, but for the use case of helping R users migrate to Julia I personally would accept the wart.

Here is the code without docs:

function findpeaks(x::AbstractVector{T};
    nups::Int=1,
    ndowns::Int=nups,
    zerostr::Char='0',
    peakpat=nothing, 
    minpeakheight=typemin(T), 
    minpeakdistance::Int=1,
    threshold=zero(T),
    npeaks::Int=0,
    sortstr=false) where T
    
    zerostr ∉ ('0', '+', '-') && error("zero must be one of `0`, `-` or `+`")

    # generate the peak pattern with no of ups and downs or use provided one
    peakpat = Regex(peakpat === nothing ? "[+]{$nups,}[-]{$ndowns,}" : peakpat)

    # transform x into a "+-+...-+-" character string
    xs = String(map(diff(x)) do e
        e < 0 && return '-'
        e > 0 && return '+'
        return zerostr
    end)

    # find index positions and maximum values
    peaks = map(findall(peakpat, xs)) do m
        v, i = findmax(@view x[m])
        (;value=v, idx=first(m) + i - 1, start=first(m), stop=last(m) + 1)
    end

    # eliminate peaks that are too low
    filter!(peaks) do p
        p.value >= minpeakheight && p.value - max(x[p.start], x[p.stop]) >= threshold
    end

    # sort according to peak height
    if sortstr || minpeakdistance > 1
        sort!(peaks, by=x -> x.value; rev=true)
    end

    # find peaks sufficiently distant
    if minpeakdistance > 1
        removal = falses(length(peaks))
        for i in 1:length(peaks)
            removal[i] && continue
            for j in 1:length(peaks)
                removal[j] && continue
                dist = abs(peaks[i].idx - peaks[j].idx)
                removal[j] = 0 < dist < minpeakdistance 
            end
        end
        deleteat!(peaks, removal)
    end

    # Return only the first 'npeaks' peaks
    npeaks > 0 && resize!(peaks, min(length(peaks), npeaks))

    return peaks
end

[halleysfifthinc]

I agree; it doesn't look like our current functionality would be able to replicate that. Out of curiosity, could you share what your typical function call for findpeaks (and desired output if convenient) would be?

The original author gave his OK

Was that OK also acknowledging the change of license? The Pracma package is GPL3 while Peaks.jl is MIT.

Assuming the change of license was approved, please open a PR! This definitely has some great new/different functionality.

[kongdd]

@laborg bravo. This is just what I need to translate phenofit into phenofit.jl (https://github.com/eco-hydro/phenofit/blob/master/R/findpeaks.R).

[halleysfifthinc]

Closing this due to the unresolved question about the license. Having seen this GPLv3 derivative code; I cannot personally implement this function without an OK of the license change from the original pracma author.