Compatibility with negative physio traces

[m-miedema]

Peak detection returns an error if run on physio traces which have no positive values.

Expected Behavior

Ideally a large negative offset shouldn't break the peakfind_physio() operation!

Actual Behavior

I receive an error which outputs the following:

In [79]: data = operations.peakfind_physio(data, thresh=0.1, dist=100)
C:\Users\memie\miniconda3\envs\physiopy2023\lib\site-packages\peakdet\operations.py:131: RuntimeWarning: Mean of empty slice.
  cdist = np.diff(locs).mean() // 2
C:\Users\memie\miniconda3\envs\physiopy2023\lib\site-packages\numpy\core\_methods.py:192: RuntimeWarning: invalid value encountered in scalar divide
  ret = ret.dtype.type(ret / rcount)
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[79], line 1
----> 1 data = operations.peakfind_physio(data, thresh=0.1, dist=100)

File ~\miniconda3\envs\physiopy2023\lib\site-packages\peakdet\utils.py:40, in make_operation.<locals>.get_call.<locals>.wrapper(data, *args, **kwargs)
     37 params = sig.bind(data, *args, **kwargs).arguments
     39 # actually run function on data
---> 40 data = func(data, *args, **kwargs)
     42 # it shouldn't be, but don't bother appending to history if it is
     43 if data is None:

File ~\miniconda3\envs\physiopy2023\lib\site-packages\peakdet\operations.py:132, in peakfind_physio(data, thresh, dist)
    130 # second, more thorough peak detection
    131 cdist = np.diff(locs).mean() // 2
--> 132 heights = np.percentile(heights['peak_heights'], 1)
    133 locs, heights = signal.find_peaks(data[:], distance=cdist, height=heights)
    134 data._metadata['peaks'] = locs

File <__array_function__ internals>:200, in percentile(*args, **kwargs)

File ~\miniconda3\envs\physiopy2023\lib\site-packages\numpy\lib\function_base.py:4205, in percentile(a, q, axis, out, overwrite_input, method, keepdims, interpolation)
   4203 if not _quantile_is_valid(q):
   4204     raise ValueError("Percentiles must be in the range [0, 100]")
-> 4205 return _quantile_unchecked(
   4206     a, q, axis, out, overwrite_input, method, keepdims)

File ~\miniconda3\envs\physiopy2023\lib\site-packages\numpy\lib\function_base.py:4473, in _quantile_unchecked(a, q, axis, out, overwrite_input, method, keepdims)
   4465 def _quantile_unchecked(a,
   4466                         q,
   4467                         axis=None,
   (...)
   4470                         method="linear",
   4471                         keepdims=False):
   4472     """Assumes that q is in [0, 1], and is an ndarray"""
-> 4473     return _ureduce(a,
   4474                     func=_quantile_ureduce_func,
   4475                     q=q,
   4476                     keepdims=keepdims,
   4477                     axis=axis,
   4478                     out=out,
   4479                     overwrite_input=overwrite_input,
   4480                     method=method)

File ~\miniconda3\envs\physiopy2023\lib\site-packages\numpy\lib\function_base.py:3752, in _ureduce(a, func, keepdims, **kwargs)
   3749             index_out = (0, ) * nd
   3750             kwargs['out'] = out[(Ellipsis, ) + index_out]
-> 3752 r = func(a, **kwargs)
   3754 if out is not None:
   3755     return out

File ~\miniconda3\envs\physiopy2023\lib\site-packages\numpy\lib\function_base.py:4639, in _quantile_ureduce_func(a, q, axis, out, overwrite_input, method)
   4637     else:
   4638         arr = a.copy()
-> 4639 result = _quantile(arr,
   4640                    quantiles=q,
   4641                    axis=axis,
   4642                    method=method,
   4643                    out=out)
   4644 return result

File ~\miniconda3\envs\physiopy2023\lib\site-packages\numpy\lib\function_base.py:4745, in _quantile(arr, quantiles, axis, method, out)
   4737 arr.partition(
   4738     np.unique(np.concatenate(([0, -1],
   4739                               previous_indexes.ravel(),
   4740                               next_indexes.ravel(),
   4741                               ))),
   4742     axis=DATA_AXIS)
   4743 if np.issubdtype(arr.dtype, np.inexact):
   4744     slices_having_nans = np.isnan(
-> 4745         take(arr, indices=-1, axis=DATA_AXIS)
   4746     )
   4747 else:
   4748     slices_having_nans = None

File <__array_function__ internals>:200, in take(*args, **kwargs)

File ~\miniconda3\envs\physiopy2023\lib\site-packages\numpy\core\fromnumeric.py:190, in take(a, indices, axis, out, mode)
     93 @array_function_dispatch(_take_dispatcher)
     94 def take(a, indices, axis=None, out=None, mode='raise'):
     95     """
     96     Take elements from an array along an axis.
     97
   (...)
    188            [5, 7]])
    189     """
--> 190     return _wrapfunc(a, 'take', indices, axis=axis, out=out, mode=mode)

File ~\miniconda3\envs\physiopy2023\lib\site-packages\numpy\core\fromnumeric.py:57, in _wrapfunc(obj, method, *args, **kwds)
     54     return _wrapit(obj, method, *args, **kwds)
     56 try:
---> 57     return bound(*args, **kwds)
     58 except TypeError:
     59     # A TypeError occurs if the object does have such a method in its
     60     # class, but its signature is not identical to that of NumPy's. This
   (...)
     64     # Call _wrapit from within the except clause to ensure a potential
     65     # exception has a traceback chain.
     66     return _wrapit(obj, method, *args, **kwds)

IndexError: cannot do a non-empty take from an empty axes.``

Steps to Reproduce the Problem

1. load in default RESP data or generate any kind of 'negative' data to read in as physio object
2. run operations.peakfind_physio() on this physio object (error seems replicable across different thresh, dist arguments)  

Specifications

- Python version: 3.9.16
- phys2bids version:  0.2.0rc1
- Platform: Windows (Anaconda Powershell Prompt, Ipython)

Possible solution

The easiest thing could be to check and possibly change offsets in the data on read-in, but the better thing seems like finding where the issue is coming from in the percentile call.

[jstac]

Thanks @oyamad . I agree. I'll fix this when I get some time.

[shlff]

Thanks @jstac and @oyamad .

My comment on point 1 is that It is a smart change that makes the most of the analytical solution of the stationary distribution of the Markov chain when the stochastic matrix is positive.

For example, let the stochastic matrix be

$$ P = \left( \begin{matrix} 1 - \lambda & \lambda \ \alpha & 1 - \alpha \end{matrix} \right) $$

if $\alpha \in (0, 1)$ and $\lambda \in (0, 1)$, then $P$ has a unique stationary distribution.

However it cannot handle the case when $\alpha$ or $\lambda$ takes the boundary values, that is, $\alpha, \lambda$ takes $0, 1$.

I suggest that we turn this change or the original one into an exercise.