Merge pull request #23 from rmarkello/find_peaks

[REF, FIX, ENH] General updates + relicensing

README.rst

@@ -11,8 +11,8 @@ photoplethysmography, or electrocardiogram (ECG/EKG) monitors.
    :target: https://codecov.io/gh/rmarkello/peakdet
 .. image:: https://readthedocs.org/projects/peakdet/badge/?version=latest
    :target: http://peakdet.readthedocs.io/en/latest
-.. image:: https://img.shields.io/badge/License-GPL%20v3-blue.svg
-   :target: https://www.gnu.org/licenses/gpl-3.0
+.. image:: https://img.shields.io/badge/license-Apache%202-blue.svg
+   :target: http://www.apache.org/licenses/LICENSE-2.0
 
 .. _overview:
 
@@ -55,31 +55,12 @@ anything you might like to work on. Alternatively, if you've found a bug, are
 experiencing a problem, or have a question, create a new issue with some
 information about it!
 
-.. _acknowledgments:
-
-Acknowledgments
----------------
-
-While this package was initially created in 2016, some of the behind-the-scenes
-functions in the project were adopted from the `PhysIO <https://github.com/
-translationalneuromodeling/tapas/tree/master/PhysIO>`_ MATLAB toolbox. As such,
-if you use this code it would be good to (1) provide a link back to the
-``peakdet`` repository with the version of the code used, and (2) cite `their
-paper <http://www.sciencedirect.com/science/article/pii/S016502701630259X>`_
-([1]_):
-
-.. [1] Kasper, L., Bollmann, S., Diaconescu, A. O., Hutton, C., Heinzle, J.,
-   Iglesias, S., ... & Stephan, K. E. (2017). The PhysIO toolbox for modeling
-   physiological noise in fMRI data. Journal of Neuroscience Methods, 276,
-   56-72.
-
 .. _licensing:
 
 License Information
 -------------------
 
-This codebase is licensed under the GNU General Public License version 3.
-The full license can be found in the `LICENSE <https://github.com/rmarkello/
-peakdet/blob/master/LICENSE>`_ file in the ``peakdet`` distribution.
-
-All trademarks referenced herein are property of their respective holders.
+This codebase is licensed under the Apache License, Version 2.0. The full
+license can be found in the `LICENSE <https://github.com/rmarkello/peakdet/
+blob/master/LICENSE>`_ file in the ``peakdet`` distribution. You may also
+obtain a copy of the license at: http://www.apache.org/licenses/LICENSE-2.0.

docs/conf.py

@@ -101,7 +101,6 @@
     'matplotlib': ('https://matplotlib.org', None),
     'numpy': ('https://docs.scipy.org/doc/numpy', None),
     'scipy': ('https://docs.scipy.org/doc/scipy/reference', None),
-    'sklearn': ('http://scikit-learn.org/stable', None),
 }
 
 doctest_global_setup = """

docs/user_guide/loading.rst

@@ -33,7 +33,7 @@ This way, the loading of the data is retained in the object's history:
 .. doctest::
 
     >>> ecg.history
-    [('load_physio', {'data': 'ECG.csv', 'dtype': None, 'fs': 1000.0, 'history': None})]
+    [('load_physio', {'allow_pickle': False, 'data': 'ECG.csv', 'dtype': None, 'fs': 1000.0, 'history': None})]
 
 There are also a number of functions for loading data from "standard" formats.
 For example, if your data were collected using the `rtpeaks <https://github.com

docs/user_guide/metrics.rst

@@ -7,6 +7,7 @@
     data = operations.interpolate_physio(data, target_fs=250.0)
     data = operations.filter_physio(data, cutoffs=1.0, method='lowpass')
     data = operations.peakfind_physio(data, thresh=0.1, dist=100)
+    data = operations.reject_peaks(data, [81441, 160786, 163225])
 
 Deriving physiological metrics
 ------------------------------
@@ -26,7 +27,7 @@ directly to the :py:class:`~.HRV` class:
     >>> from peakdet import HRV
     >>> metrics = HRV(data)
     >>> print(f'{metrics.rmssd:.2f} ms')
-    27.75 ms
+    26.61 ms
 
 The :py:class:`~.HRV` class contains many common heart rate variability metrics
 including the root mean square of successive differences, as shown above. It

docs/user_guide/physio.rst

@@ -68,9 +68,9 @@ store the indices of the detected :py:attr:`~.Physio.peaks` and
 
     >>> phys = operations.peakfind_physio(phys)
     >>> phys.peaks
-    array([ 477, 2120, 3253, 4128])
+    array([ 477, 2120, 3253])
     >>> phys.troughs
-    array([1413, 2611, 3756])
+    array([1413, 2611])
 
 Next, we'll move on to how you can load your data into a :py:class:`~.Physio`
 object in a more reproducible manner. Feel free to refer to the :ref:`api_ref`

docs/user_guide/saving.rst

@@ -31,13 +31,16 @@ If later on you want to reload the processed data you can do so:
 .. doctest::
 
     >>> from peakdet import load_physio
-    >>> data = load_physio('out.phys')
+    >>> data = load_physio('out.phys', allow_pickle=True)
 
-.. note::
+.. warning::
 
     :py:func:`peakdet.save_physio` uses :py:func:`numpy.savez` to save the data
     objects, meaning the generated files can actually be quite a bit larger
-    than the original data, themselves!
+    than the original data, themselves! Moreover, you NEED to pass the
+    ``allow_pickle=True`` parameter; this is turned off by default for safety
+    reasons, as you should never load pickle files not generated by a trusted
+    source.
 
 Saving history
 ^^^^^^^^^^^^^^

peakdet/editor.py

@@ -112,19 +112,18 @@ def undo(self):
         func, peaks = self.data._history.pop()
 
         if func == 'reject_peaks':
-            self.data._metadata.reject = np.setdiff1d(
-                self.data._metadata.reject, peaks['remove']
+            self.data._metadata['reject'] = np.setdiff1d(
+                self.data._metadata['reject'], peaks['remove']
             )
             self.rejected.difference_update(peaks['remove'])
         elif func == 'delete_peaks':
-            self.data._metadata.peaks = np.insert(
-                self.data._metadata.peaks,
-                np.searchsorted(self.data._metadata.peaks, peaks['remove']),
+            self.data._metadata['peaks'] = np.insert(
+                self.data._metadata['peaks'],
+                np.searchsorted(self.data._metadata['peaks'], peaks['remove']),
                 peaks['remove']
             )
             self.deleted.difference_update(peaks['remove'])
 
-        self.data._metadata.troughs = utils.check_troughs(self.data,
-                                                          self.data.peaks,
-                                                          self.data.troughs)
+        self.data._metadata['troughs'] = utils.check_troughs(self.data,
+                                                             self.data.peaks)
         self.plot_signals()

peakdet/info.py

@@ -5,7 +5,7 @@
 MAINTAINER = 'Ross Markello'
 EMAIL = 'rossmarkello@gmail.com'
 VERSION = __version__
-LICENSE = 'GPLv3'
+LICENSE = 'Apache-2.0'
 DESCRIPTION = """\
 A toolbox for reproducible physiological data analysis\
 """
@@ -18,7 +18,6 @@
 INSTALL_REQUIRES = [
     'matplotlib',
     'numpy',
-    'scikit-learn',
     'scipy',
 ]
 
@@ -49,6 +48,6 @@
 CLASSIFIERS = [
     'Development Status :: 3 - Alpha',
     'Intended Audience :: Science/Research',
-    'License :: OSI Approved :: GNU General Public License v3 (GPLv3)',
+    'License :: OSI Approved :: Apache Software License',
     'Programming Language :: Python :: 3.6',
 ]

peakdet/io.py

@@ -12,7 +12,8 @@
 EXPECTED = ['data', 'fs', 'history', 'metadata']
 
 
-def load_physio(data, *, fs=None, dtype=None, history=None):
+def load_physio(data, *, fs=None, dtype=None, history=None,
+                allow_pickle=False):
     """
     Returns `Physio` object with provided data
 
@@ -26,6 +27,9 @@ def load_physio(data, *, fs=None, dtype=None, history=None):
         Data type to convert `data` to, if conversion needed. Default: None
     history : list of tuples, optional
         Functions that have been performed on `data`. Default: None
+    allow_pickle : bool, optional
+        Whether to allow loading if `data` contains pickled objects. Default:
+        False
 
     Returns
     -------
@@ -42,7 +46,7 @@ def load_physio(data, *, fs=None, dtype=None, history=None):
     # load it as a plain text file and instantiate a history
     if isinstance(data, str):
         try:
-            inp = dict(np.load(data))
+            inp = dict(np.load(data, allow_pickle=allow_pickle))
             for attr in EXPECTED:
                 try:
                     inp[attr] = inp[attr].dtype.type(inp[attr])
@@ -53,7 +57,7 @@ def load_physio(data, *, fs=None, dtype=None, history=None):
             # fix history, which needs to be list-of-tuple
             if inp['history'] is not None:
                 inp['history'] = list(map(tuple, inp['history']))
-        except IOError:
+        except (IOError, OSError, ValueError):
             inp = dict(data=np.loadtxt(data),
                        history=[utils._get_call(exclude=[])])
         phys = physio.Physio(**inp)
@@ -86,33 +90,33 @@ def load_physio(data, *, fs=None, dtype=None, history=None):
     return phys
 
 
-def save_physio(file, data):
+def save_physio(fname, data):
     """
     Saves `data` to `fname`
 
     Parameters
     ----------
-    file : str
+    fname : str
         Path to output file; .phys will be appended if necessary
     data : Physio_like
         Data to be saved to file
 
     Returns
     -------
-    file : str
+    fname : str
         Full filepath to saved output
     """
 
     from peakdet.utils import check_physio
 
     data = check_physio(data)
-    file += '.phys' if not file.endswith('.phys') else ''
-    with open(file, 'wb') as dest:
+    fname += '.phys' if not fname.endswith('.phys') else ''
+    with open(fname, 'wb') as dest:
         hist = data.history if data.history != [] else None
         np.savez_compressed(dest, data=data.data, fs=data.fs,
                             history=hist, metadata=data._metadata)
 
-    return file
+    return fname
 
 
 def load_history(file, verbose=False):

peakdet/modalities.py

@@ -12,14 +12,14 @@ def iHR(self, step=1, start=0, end=None, TR=None):
         time = np.arange(start - mod, end + mod + 1, self.TR, dtype='int')
         HR = np.zeros(len(time) - step)
 
-        for l in range(step, time.size):
-            inds = np.logical_and(self.rrtime >= time[l - step],
-                                  self.rrtime < time[l])
+        for tpoint in range(step, time.size):
+            inds = np.logical_and(self.rrtime >= time[tpoint - step],
+                                  self.rrtime < time[tpoint])
             relevant = self.rrint[inds]
 
             if relevant.size == 0:
                 continue
-            HR[l - step] = (60 / relevant).mean()
+            HR[tpoint - step] = (60 / relevant).mean()
 
         return HR
 

peakdet/operations.py

@@ -5,8 +5,7 @@
 
 import matplotlib.pyplot as plt
 import numpy as np
-from scipy import signal
-from scipy.interpolate import InterpolatedUnivariateSpline
+from scipy import interpolate, signal
 from peakdet import editor, utils
 
 
@@ -63,7 +62,7 @@ def filter_physio(data, cutoffs, method, *, order=3):
 
 
 @utils.make_operation()
-def interpolate_physio(data, target_fs):
+def interpolate_physio(data, target_fs, *, kind='cubic'):
     """
     Interpolates `data` to desired sampling rate `target_fs`
 
@@ -73,6 +72,9 @@ def interpolate_physio(data, target_fs):
         Input physiological data to be interpolated
     target_fs : float
         Desired sampling rate for `data`
+    kind : str or int, optional
+        Type of interpolation to perform. Must be one of available kinds in
+        :func:`scipy.interpolate.interp1d`. Default: 'cubic'
 
     Returns
     -------
@@ -89,7 +91,7 @@ def interpolate_physio(data, target_fs):
     t_new = np.linspace(0, len(data) / data.fs, int(len(data) * factor))
 
     # interpolate data and generate new Physio object
-    interp = InterpolatedUnivariateSpline(t_orig, data[:])(t_new)
+    interp = interpolate.interp1d(t_orig, data, kind=kind)(t_new)
     interp = utils.new_physio_like(data, interp, fs=target_fs)
 
     return interp
@@ -122,12 +124,16 @@ def peakfind_physio(data, *, thresh=0.2, dist=None):
 
     # first pass peak detection to get approximate distance between peaks
     cdist = data.fs // 4 if dist is None else dist
-    locs = utils.find_peaks(data, dist=cdist, thresh=thresh)
-    cdist = np.diff(locs).mean() // 2
+    thresh = np.squeeze(np.diff(np.percentile(data, [5, 95]))) * thresh
+    locs, heights = signal.find_peaks(data[:], distance=cdist, height=thresh)
+
     # second, more thorough peak detection
-    data._metadata.peaks = utils.find_peaks(data, dist=cdist, thresh=thresh)
+    cdist = np.diff(locs).mean() // 2
+    heights = np.percentile(heights['peak_heights'], 1)
+    locs, heights = signal.find_peaks(data[:], distance=cdist, height=heights)
+    data._metadata['peaks'] = locs
     # perform trough detection based on detected peaks
-    data._metadata.troughs = utils.check_troughs(data, data.peaks, [])
+    data._metadata['troughs'] = utils.check_troughs(data, data.peaks)
 
     return data
 
@@ -148,9 +154,8 @@ def delete_peaks(data, remove):
     """
 
     data = utils.check_physio(data, ensure_fs=False, copy=True)
-    data._metadata.peaks = np.setdiff1d(data._metadata.peaks, remove)
-    data._metadata.troughs = utils.check_troughs(data, data.peaks,
-                                                 data.troughs)
+    data._metadata['peaks'] = np.setdiff1d(data._metadata['peaks'], remove)
+    data._metadata['troughs'] = utils.check_troughs(data, data.peaks)
 
     return data
 
@@ -171,9 +176,8 @@ def reject_peaks(data, remove):
     """
 
     data = utils.check_physio(data, ensure_fs=False, copy=True)
-    data._metadata.reject = np.append(data._metadata.reject, remove)
-    data._metadata.troughs = utils.check_troughs(data, data.peaks,
-                                                 data.troughs)
+    data._metadata['reject'] = np.append(data._metadata['reject'], remove)
+    data._metadata['troughs'] = utils.check_troughs(data, data.peaks)
 
     return data
 

peakdet/physio.py

@@ -4,7 +4,6 @@
 """
 
 import numpy as np
-from sklearn.utils import Bunch
 
 
 class Physio():
@@ -63,11 +62,11 @@ def __init__(self, data, fs=None, history=None, metadata=None):
                     except TypeError:
                         raise TypeError('Provided metadata must be dict-like'
                                         'with integer array entries.')
-            self._metadata = Bunch(**metadata)
+            self._metadata = dict(**metadata)
         else:
-            self._metadata = Bunch(peaks=np.empty(0, dtype=int),
-                                   troughs=np.empty(0, dtype=int),
-                                   reject=np.empty(0, dtype=int))
+            self._metadata = dict(peaks=np.empty(0, dtype=int),
+                                  troughs=np.empty(0, dtype=int),
+                                  reject=np.empty(0, dtype=int))
 
     def __array__(self):
         return self.data
@@ -110,10 +109,10 @@ def peaks(self):
     @property
     def troughs(self):
         """ Indices of detected troughs in `data` """
-        return self._metadata.troughs
+        return self._metadata['troughs']
 
     @property
     def _masked(self):
-        return np.ma.masked_array(self._metadata.peaks,
-                                  mask=np.isin(self._metadata.peaks,
-                                               self._metadata.reject))
+        return np.ma.masked_array(self._metadata['peaks'],
+                                  mask=np.isin(self._metadata['peaks'],
+                                               self._metadata['reject']))