Skip to content

ENH: change to use array view in signal/_peak_finding.py #3409

Merged
merged 3 commits into from Mar 12, 2014

5 participants

@yotam
yotam commented Feb 27, 2014

Unnecessary call to arange and copying of data. Changed to call scoreatpercentile on a view of the array row_one .

@coveralls

Coverage Status

Coverage remained the same when pulling 09c6636 on yotam:master into 4f4b037 on scipy:master.

@yotam yotam closed this Feb 28, 2014
@yotam yotam reopened this Feb 28, 2014
@coveralls

Coverage Status

Coverage remained the same when pulling 4f0808a on yotam:master into 3879d21 on scipy:master.

@coveralls

Coverage Status

Coverage remained the same when pulling 4f0808a on yotam:master into 3879d21 on scipy:master.

@WarrenWeckesser

It looked harmless at first, but there are cases where this PR introduces a change in behavior.

Suppose num_points = 100 and window_size = 5, so hf_window = 2.5. Consider the first iteration of the loop at line 394, when ind = 0. Because arange(0, 2.5) is array([0., 1., 2.]), the old code creates window = arange([0, 1, 2]). The code in this PR creates window_start = 0, window_end = 2, so the slice row_one[window_start:window_end] is not the same as the old row_one[window].

I haven't looked at what is going on closely enough to know if this is a bug in the PR, or if the PR fixes a bug in the existing code. :)

@yotam
yotam commented Mar 3, 2014

Good point. I didn't spot that.

The method implements this paper: http://basic.northwestern.edu/publications/peakdetection/

There's only mention of a sliding window in the original paper, no details about indexing.

I can try and pick through their code (there's equivalent code in the file IdentifyMajorPeaks.R) which itself may have bugs. Or I can just change this to match the old behaviour.

@rgommers
SciPy member
rgommers commented Mar 9, 2014

@yotam you can't look at R code because it's GPL. Trying to figure out if this is a bug would be quite helpful if you have the time for that. If not, please make this match the old behavior and open an issue with a note to investigate the possible bug later.

@WarrenWeckesser

It looks like both the original code and this PR have an off-by-one error. I'm assuming that when window_size is odd, the window should be centered around ind. So, for example, if ind is 10, and window_size is 5, the indices of the window should be [8, 9, 10, 11, 12]. But the original code and this PR create a window with indices [7, 8, 9, 10, 11].

I suggest removing all use of floating point values when working with the indices. So first make sure window_size is an integer by doing:

window_size = int(window_size)

before entering the loop. Use divmod to compute hf_window, and at the same time note whether window_size is odd or even:

hf_window, odd = divmod(window_size, 2)

In the loop, compute window_start and window_end as

    window_start = max(ind - hf_window, 0)
    window_end = min(ind + hf_window + odd, num_points)

When ind is away from the edges, the window is centered around ind when window_size is odd. When window_size is even, the extra sample is taken from the left.

@yotam
yotam commented Mar 10, 2014

Thank you both for the feedback and analysis. I'll update the PR this evening.

@coveralls

Coverage Status

Coverage remained the same when pulling ac545ed on yotam:master into 3879d21 on scipy:master.

@WarrenWeckesser WarrenWeckesser commented on an outdated diff Mar 11, 2014
scipy/signal/_peak_finding.py
#Filter based on SNR
row_one = cwt[0, :]
noises = np.zeros_like(row_one)
for ind, val in enumerate(row_one):
- window = np.arange(max([ind - hf_window, 0]), min([ind + hf_window, num_points]))
- window = window.astype(int)
- noises[ind] = scoreatpercentile(row_one[window], per=noise_perc)
+ window_start = max(ind - hf_window, 0)
+ window_end = min(ind + hf_window + odd, num_points)
+ noises[ind] = scoreatpercentile(row_one[window_start:window_end], per=noise_perc)

This line is now too long. PEP8 limits lines to at most 79 characters.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@WarrenWeckesser

One last tweak needed: line 399 is too long (http://legacy.python.org/dev/peps/pep-0008/). Once that is fixed, this looks ready to merge.

@coveralls

Coverage Status

Coverage remained the same when pulling ee6da29 on yotam:master into 3879d21 on scipy:master.

@WarrenWeckesser WarrenWeckesser merged commit 4aa534f into scipy:master Mar 12, 2014

1 check passed

Details default The Travis CI build passed
@WarrenWeckesser

@yotam, thanks for spotting and fixing this.

@WarrenWeckesser WarrenWeckesser added this to the 0.15.0 milestone Mar 12, 2014
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Something went wrong with that request. Please try again.