findpeaks bug? ("invalid MinPeakDistance" error)

[William Frane]

Hello,

I've been using the findpeaks function to individually process hundreds of datasets (each containing around 400 values) and yesterday it produced an odd error when applied to one of them. The error was:
??? Error using ==> findpeaks at 85
Invalid MinPeakDistance. Set MinPeakDistance in the range (0,1).

The odd thing is that I don't set MinPeakDistance (although I do set MinPeakHeight) so presumably it's set to 1 by default (as indicated in the documentation for findpeaks). As of yet this error has not occurred for any other dataset.

Using the first 12 values of the dataset that caused the error, it can be reproduced as follows:
a = [58 -41 -8 4 -3 7 7 -6 -2 -1 2 -2];
p = findpeaks(a,'minpeakheight',2*std(a));

It seems that I'm not the only person to have encountered this problem (see http://www.mathworks.com/matlabcentral/newsreader/view_thread/236916), but I've been unable to find a solution to it or even a clear cause.

Any suggestions would be appreciated (I'm using MATLAB R2008a).

W Frane

[dpb]

I don't have findpeaks() so can't look but I'll venture a couple of
suggestions...

1. What is the code at and around line 85 that generates the error?
Should give a clue as to what's going on...

2. Does it error if you take the same data set but instead of 58 at the
first location make it under the criterion level (or raise the height to
say 3*std(s))? IOW, I'm wondering if it's the max value as the first
point not having a lesser value to its left that's causing the problem
by an end effect not being accounted for properly.

After determining the above, the results therefrom may raise other ideas...

--

[William Frane]

> I don't have findpeaks() so can't look but I'll venture a couple of
> suggestions...
>
> 1. What is the code at and around line 85 that generates the error?
> Should give a clue as to what's going on...

Since the line 85 error message is for findpeaks itself, not the function calling it, I've only glanced at the code in that area of findpeaks. However, I didn't see anything in that vicinity that would explain the error.

> 2. Does it error if you take the same data set but instead of 58 at the
> first location make it under the criterion level (or raise the height to
> say 3*std(s))? IOW, I'm wondering if it's the max value as the first
> point not having a lesser value to its left that's causing the problem
> by an end effect not being accounted for properly.

Changing the first value such that it's under the peak threshold removes the error. If I gradually raise the value specified for MinPeakHeight (when applying findpeaks to the original data), the error persists until it's replaced by the "no peaks found" warning (i.e., when MinPeakHeight exceeds the maximum value in the input data).

I've yet to find another dataset that causes findpeaks to generate this error. I've tried a few sets with maximum values in the first element and different combinations of signs for each element, but haven't been able to trigger the error.

[dpb]

On 5/29/2012 12:34 PM, William Frane wrote:
>> I don't have findpeaks() so can't look but I'll venture a couple of
>> suggestions...
>>
>> 1. What is the code at and around line 85 that generates the error?
>> Should give a clue as to what's going on...
>
> Since the line 85 error message is for findpeaks itself, not the
> function calling it, I've only glanced at the code in that area of
> findpeaks. However, I didn't see anything in that vicinity that would
> explain the error.

Well, unless one looks it's unlikely will find anything, indeed. :)

Why not post the area of code in that region w/ the specific message?
More eyes may well spot something.

>> 2. Does it error if you take the same data set but instead of 58 at
>> the first location make it under the criterion level (or raise the
>> height to say 3*std(s))? IOW, I'm wondering if it's the max value as
>> the first point not having a lesser value to its left that's causing
>> the problem by an end effect not being accounted for properly.
>
> Changing the first value such that it's under the peak threshold removes
> the error. If I gradually raise the value specified for MinPeakHeight
> (when applying findpeaks to the original data), the error persists until
> it's replaced by the "no peaks found" warning (i.e., when MinPeakHeight
> exceeds the maximum value in the input data).
> I've yet to find another dataset that causes findpeaks to generate this
> error. I've tried a few sets with maximum values in the first element
> and different combinations of signs for each element, but haven't been
> able to trigger the error.

I'd suggest since findpeaks is a toolbox function to submit the sample
dataset to official support as a bug report. Having a repeatable
dataset should enable TMW to uncover the problem. (From your
description of behavior w/ changing threshold, etc., it still sounds to
me like an end-effect that may be related to the actual values as well.
The mixture of positive and negative values is another area I wonder
if there's a logic error concerning in the implementation. What happens
if you take the same dataset but add abs(min(x))+C) where C is a
constant >=0? IOW, same shape but all positive values.)

--

[William Frane]

> Well, unless one looks it's unlikely will find anything, indeed. :)
>
> Why not post the area of code in that region w/ the specific message?
> More eyes may well spot something.

True enough! The code around line 85 in findpeaks is as follows (line 85 is the line on which the error function is called):

% validate value of Pd and set default values for Pd and Np
[Pd,Np] = setvalues(Pd,Np,L);
if(Pd >=Indx(L))
pdmsgid = generatemsgid('largeMinPeakDistance');
error(pdmsgid,strcat('Invalid MinPeakDistance. Set MinPeakDistance in the range (',num2str(0), ',', num2str(Indx(L)), ').'));
else
% Call a nested function that finds the peaks
end

Pd is the peak distance; if not specified by the user (as in my case) it is set to 1 by setvalues. Indx is a list of the indices of all the input data values greater than MinPeakHeight. L is the length of this vector.

For the given data set, Pd should be 1 (the default) and Indx(L) should also be 1 (since the first element is the only one that exceeds MinPeakHeight, Indx should be [1] and L, being effectively the length of Indx, will also be 1). Thus Pd >= Indx(L) will be true, triggering the error. The mechanics of what's causing the error are simple enough, but I'm curious about the logic behind them.

> I'd suggest since findpeaks is a toolbox function to submit the sample
> dataset to official support as a bug report. Having a repeatable
> dataset should enable TMW to uncover the problem. (From your
> description of behavior w/ changing threshold, etc., it still sounds to
> me like an end-effect that may be related to the actual values as well.
> The mixture of positive and negative values is another area I wonder
> if there's a logic error concerning in the implementation. What happens
> if you take the same dataset but add abs(min(x))+C) where C is a
> constant >=0? IOW, same shape but all positive values.)

Adding an offset to all values removes the error.

Another dataset that triggers the error is:
a = [100 1 1 1];

p = findpeaks(a,'minpeakheight',2*std(a));

Interestingly, the command sets
b = [100 1 1 1 100];
p = findpeaks(b,'minpeakheight',2*std(a)); % Using same MinPH as above
and
c= [100 1 1 1];
p = findpeaks(c);
don't trigger the error (the result for each is that no peaks are found). This appears to indicate that the error occurs only with a single peak above MinPeakHeight at the start of the input dataset. Based on datasets b and c above, consistent behaviour in the error case would be for findpeaks to simply report that no peaks were found.

I think I'll take your suggestion and file this as a bug report, but if you (or anyone else) have any additional thoughts on this I'd still welcome them.

W Frane

[dpb]

On 5/30/2012 9:05 AM, William Frane wrote:
...

> Adding an offset to all values removes the error.

...

OK, that shows that the sign is significant as well as the shape as I
suspected.

I've got to run now; I'll look at the code snippet a little but not sure
will be able to do much w/o being able to run it.

If it were me and a solution were important for immediate work, I'd use
debug and work my way through the logic to understand what test actually
triggered the message. With that it should be possible to see what is
going on that caused it and where the logic of what they're testing got
fouled up.

If it's only of academic rather than immediate practical concern then
the bug report should be satisfactory (assuming, of course, you have a
support contract so TMW will pay any attention).

--

[dpb]

On 5/30/2012 9:05 AM, William Frane wrote:
...

> % validate value of Pd and set default values for Pd and Np
> [Pd,Np] = setvalues(Pd,Np,L);
if(Pd >=Indx(L))
> pdmsgid = generatemsgid('largeMinPeakDistance');
> error(pdmsgid,strcat('Invalid MinPeakDistance. Set MinPeakDistance in
> the range (',num2str(0), ',', num2str(Indx(L)), ').'));
> else % Call a nested function that finds the peaks
> end
>
> Pd is the peak distance; if not specified by the user (as in my case) it
> is set to 1 by setvalues. Indx is a list of the indices of all the input
> data values greater than MinPeakHeight. L is the length of this vector.
>
> For the given data set, Pd should be 1 (the default) and Indx(L) should
> also be 1 (since the first element is the only one that exceeds
> MinPeakHeight, Indx should be [1] and L, being effectively the length of
> Indx, will also be 1). Thus Pd >= Indx(L) will be true, triggering the
> error. The mechanics of what's causing the error are simple enough, but
> I'm curious about the logic behind them.

...

OK, I had time to look at the logic a little--it's an end effect
triggered w/ the choice of MinPeakHeight and the dataset as figured.

Consider your original dataset and level again==

a = [58 -41 -8 4 -3 7 7 -6 -2 -1 2 -2];

p = findpeaks(a,'minpeakheight',2*std(a));

hMin=2*std(a)
hMin =
43.7441
idxa=find(a>=hMin)
idxa =
1

I asked about the magnitude shift;

C=0;
b=a+abs(min(a))+C
b =
59.25 -39.75 -6.75 5.25 -1.75 8.25 8.25 -4.75 -0.75 0.25 3.25 -0.75

idxb=find(b>hMin)
idxb =
1 4 6 7

So shifting the magnitude but not compensating for the mean shift makes
the problem go away.

Of course, if

C=mean(a);
idxb=find(b>=h+C)
idxb =
1

back to where was originally.

The flaw that I see it is that logic doesn't consider a single maximum
at the first point. The question then is whether a value greater than
the threshold w/o a point to the left which is less than that first
point should be qualified as a peak as it could possibly be a point on
the right side of the real peak. If it should be counted then there
needs to be a special case made for the index vector of length one.
OTOH, if magnitude of a peak is important as well as it's existence,
perhaps it should error (altho in this case I'd change the error message
to reflect the condition more accurately).

A similar problem arises if the last point is a maximum it would appear.

--

[Star Strider]

> a = [58 -41 -8 4 -3 7 7 -6 -2 -1 2 -2];
> p = findpeaks(a,'minpeakheight',2*std(a));

---------------------------------------------------------

In 2012a, ‘p’ is empty in spite of the first value being larger than ‘minpeakheight’. It doesn't generate your error. Reading the documentation, not detecting a peak at the origin also seems to be a bug.