embodied-computation-group / systole Goto Github PK

In the References section, lines 74-75, make sure the title appears as "Numba: a LLVM-based Python JIT compiler".

openjournals/joss-reviews#3832

Hi,

I have been using the Systole package recently and it is a great tool!

Two minor suggestion regarding the oxi_peaks function:
peak_enhancement -

    if peak_enhancement is True:
        # Square signal (peak enhancement)
        x = np.asarray(x) ** 2

By raising the signal to the power of 2, negative values turn to positive values causing negative peaks in the data to be detected. I recommend multiplying by the sign of the data:

    if peak_enhancement is True:
        # Square signal (peak enhancement)
        x = (np.asarray(x) ** 2) * np.sign(x)

rollingNoise - A second issue the window size of the rolling mean. When the sampling rate is too low (10Hz for example) int(sfreq * 0.05) is zero, causing the function to return all nans. I recommended changing to:

rollingNoise = max(int(sfreq * 0.05), 1)  # 0.05 second window

This would make sure that the window size is at least 1.

I can make a pull request with both changes if you want.

Hello! I am now working on adjusting my own device (Emotibit) into systole package.

for this, I am trying to modify 'read' function and 'readInwaiting' function with my own code.

I have a question about my code didn't stop after 10 seconds with below sample code.

import time
tstart=time.time()
emotibit.board.start_stream()

while time.time() - tstart < 10:
    emotibit.readInWaiting()

I am wandering what's wrong with my code...

Here is full class object with my device.

Actually, I want only fix few lines with 'readInwaiting' object... but it is not that easy!

It is now completely work, but I am trying to go on!

Thanks a lot with your beautiful work. Have a good day!

removed

Hey @LegrandNico hope you're doing great,

just saw that you were in the process of implementing recurrence matrix stuff for HRV, and thought you might be interested in this (much) faster implementation below:

https://github.com/neuropsychology/NeuroKit/blob/dev/neurokit2/complexity/complexity_recurrence.py

The neurokit one is tested against PyRQA and gives the same results, which are different from your current implementation. In order for it to match it, the np.flip() must be removed, and the subsamples must be cut short by one (I am not 100% certain that this is incorrect though, I would tend to trust pyRQA but 🤷).

Anyway, feel free to steal that code :)

Here are some benchmarks (not taking into account numba gains):

import neurokit2 as nk
import numpy as np


def rc_systole(rr: np.ndarray, m: int = 10, tau: int = 1):
    r = np.sqrt(m) * np.std(rr)  # Threshold for the Euclidean distance
    lag = (m - 1) * tau  # Lag
    j = rr.shape[0] - lag  # Dimension of the recurrence matrix

    # Initialize a (j-l) by (j-l) matrix and fill with zeros
    rc = np.zeros((j, j))

    # Iterate over the lower triangle only
    for i in range(j):
        u_i = rr[i : i + lag : tau]  # First sub-sample of RR intervals
        for ii in range(i + 1):
            u_ii = rr[ii : ii + lag : tau]  # Second sub-sample of RR intervals

            # Compare the Euclidean distance to threshold
            if np.sqrt(np.sum(np.square(u_i - u_ii))) <= r:
                rc[i, ii] = 1

    rc = rc + rc.T - np.diag(np.diag(rc))  # Make the matrix symmetric

    return rc


signal = nk.signal_simulate(duration=5, sampling_rate=500, frequency=[5, 6], noise=0.5)

# Systole (pure python)
%timeit rc_systole(signal)
> 14.3 s ± 1.57 s per loop (mean ± std. dev. of 7 runs, 1 loop each)

# NeuroKit
%timeit nk.complexity_recurrence(signal, dimension=10, tolerance = np.sqrt(10) * np.std(signal))
> 81.4 ms ± 1.87 ms per loop (mean ± std. dev. of 7 runs, 10 loops each

# PyRQA (but also incldues the computation of the indices so not really the best comparison)
%timeit nk.complexity_rqa(signal, delay=1, dimension=10, tolerance=np.sqrt(10) * np.std(signal))
> 181 ms ± 1.79 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

# Compare outputs (NeuroKit vs. PyRQA)
rc2, _ = nk.complexity_recurrence(signal, dimension=10, tolerance = np.sqrt(10) * np.std(signal))
_, rc3 = nk.complexity_rqa(signal, delay=1, dimension=10, tolerance=np.sqrt(10) * np.std(signal))

np.all(rc3['Recurrence_Matrix'] == rc2)
> True

looking forward to chatting with you sometime ☺️

Hi,

when going through the readme, I encounter the following error:

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
/var/folders/v1/lcrp_j352r16_yd3fbg7y5j40000gn/T/ipykernel_61076/2755290769.py in <module>
----> 1 from systole import import_dataset1
      2 
      3 # Import ECg recording
      4 signal = import_dataset1(modalities=['ECG']).ecg.to_numpy()

~/github/joss-reviews/systole/.venv/lib/python3.9/site-packages/systole/__init__.py in <module>
      3 from .detection import *
      4 from .hrv import *
----> 5 from .plots import *
      6 from .utils import *
      7 

~/github/joss-reviews/systole/.venv/lib/python3.9/site-packages/systole/plots/__init__.py in <module>
      1 """Plotting functions."""
      2 from .plot_circular import plot_circular
----> 3 from .plot_ectopic import plot_ectopic
      4 from .plot_events import plot_events
      5 from .plot_evoked import plot_evoked

~/github/joss-reviews/systole/.venv/lib/python3.9/site-packages/systole/plots/plot_ectopic.py in <module>
      4 
      5 import numpy as np
----> 6 from bokeh.plotting.figure import Figure
      7 from matplotlib.axes import Axes
      8 

ModuleNotFoundError: No module named 'bokeh'

Given that poth matplotlib and bokeh are imported when importing any function from this library, I recommend adding them both as actual requirements; otherwise, it will cause issues for users who install these in fresh environments.

Alternatively, the imports would have to be done conditionally depending on which backend is specified.

After updating to v0.2.3 (September 2022)
I get the following error when calling correct_rr(rr):

Cleaning the RR interval time series.
... correcting 3 missed interval(s).
... correcting 1 extra interval(s).
... correcting 89 ectopic interval(s).
... correcting 5 short interval(s).
... correcting 5 long interval(s).
Traceback (most recent call last):
  File "<string>", line 79, in <module>
TypeError: tuple indices must be integers or slices, not str

While trying to disable things one by one to debug (eg. correct_rr(rr, missed_correction=False) I got:

Cleaning the RR interval time series.
... correcting 1 extra interval(s).
... correcting 89 ectopic interval(s).
... correcting 5 short interval(s).
... correcting 5 long interval(s).
Traceback (most recent call last):
  File "<string>", line 77, in <module>
  File "/usr/lib/python3.10/site-packages/systole/correction.py", line 363, in correct_rr
    return _correct_rr(
  File "/usr/lib/python3.10/site-packages/systole/correction.py", line 263, in _correct_rr
    return clean_rr, (nMissed, nExtra, nEctopic, nShort, nLong)
UnboundLocalError: local variable 'nMissed' referenced before assignment

correct_rr(rr, n_iterations=2) doesn't work either:

Traceback (most recent call last):
  File "<string>", line 77, in <module>
TypeError: correct_rr() got an unexpected keyword argument 'n_iterations'

With version 0.2.2 correct_rr(rr, n_iterations=2, missed_correction=False) works as expected without errors.

On line 28 you write "python", while in the rest of the paper you write "Python". For consistency I suggest using "Python" throughout.

openjournals/joss-reviews#3832

On line 46, a reference shows up as "@lipponen:2019". Please make sure this is corrected.

openjournals/joss-reviews#3832

Thanks for creating this project. I found a minor error in the example/readme file
In the HRV section

from systole.hrv import plot_psd

should be

from systole.plotting import plot_psd

I am not sure if this is the only place the error occurs.

Thanks!

Hello @LegrandNico. I'm the editor of your JOSS submission, and will create some issues for editorial comments to your paper.

On lines 22-24, you write

Similarly, rising interest in understanding how heart-rate variability relates to mental illness, cognitive function, and physical wellbeing.

This is an incomplete sentence. Suggested rewrite: "Similarly, there is rising interest in understanding how heart-rate variability relates to mental illness, cognitive function, and physical wellbeing."

Hi,

I tried installing systole in 2 machines (macOS and CentOS). However, pip could not resolve the dependencies and keeps bringing up this error:

numba 0.54.1 requires numpy<1.21,>=1.17, but you'll have numpy 1.21.2 which is incompatible.

It does seem to be more of a numba than a systole issue but I thought you might have an idea for a solution or a workaround to get systole up and running.

I’m currently using the Systole to do artefact correction & rejection of cardiac ECG signal and it’s generally working well. I followed the guidance from the Systole manual, but labelling of bad segments functionality doesn’t seem to work at the moment.

After I open the display window with display(view.io_box, view.editor.commands_box, view.output) and choose the Rejection mode, the labelling of bad segments doesn’t work as intended.

I right click the mouse, but the transparent selection window is green instead of red. When I let go, it doesn’t turn into a grey transparent window that remains, but it disappears and a peak is placed on the local maxima. Then after I save the corrected peaks into the json file, "bad_segments": null.

Do you happen to know what might be the issue and when this functionality will become viable again?

Thanks.

• Integrate stim channels to the .save() method.

• Add parameters to control which time series to include.

Hello,

As an unprofessional first-timer, I have kept a record of few notes while trying Systole:

• An error message X must be a 1darray is more readable as a X must be a 1d array for example.

• * Document somewhere that in order to use Bokeh, output_notebook () should be called. I know it is a 3rd part add-on, but it took me hours until I noticed that call in one (of many) places in the documentation. After that I read their API ofc. But that time could have been saved?

• The current usage of CDSView is deprecated in Bokeh 3.0+ (e.g. plot_rr (show_artefacts = True) fails)

• Not all function can actually work with integer arrays (non numpy.float64). For the rr_ms type I had to explicit convert beforehand.

• A mixed time-frequency report, like wavelet transforms, would be great.

Just my 2c.
Thank you for the great framework!
I will make good use if it for my chest-band monitor -> Garmin watch -> hrv report notebooks.

Hello, thank you for make this wonderful package for cardioception researches.

In our laboratory, we have emotibit device and want to adjust this package with emotibit device.

For this, what procedure or classes need to be update?

Thank you for reading!

On line 8 you write "Systole is a library for cardiac signal analysis in Python.", whereas later, e.g., on line 33 you write "package". I think "package" is the most appropriate term to use here, so I suggest rewriting from "library" to "package" on line 8.

openjournals/joss-reviews#3832

Working on a Julia implementation of the reference paper used and noticed variation in results against your library. I narrowed it down to what I think is an error in correction.py. Current source reads:

    # Correct extra beats
    if extra_correction:
        if np.any(artefacts["extra"]):
            for this_id in np.where(artefacts["extra"])[0]:
                this_id -= nExtra
                clean_rr = correct_missed(clean_rr, this_id) # This is the incorrect function call.
                nExtra += 1
        artefacts = rr_artefacts(clean_rr)

Should read:

    # Correct extra beats
    if extra_correction:
        if np.any(artefacts["extra"]):
            for this_id in np.where(artefacts["extra"])[0]:
                this_id -= nExtra
                clean_rr = correct_extra(clean_rr, this_id) # Should be a call to correct_extra instead.
                nExtra += 1
        artefacts = rr_artefacts(clean_rr)

I can probably generate a PR to fix this if you allow contributions?

Hi Nico

Hope you're well :)

I use the correct_rr function to correct for ectopic beats as I get very high RMSSD values. Our range of RMSSD is 10-180 (literature suggested this should be between 19-75).

I apply the systole correct_rr function to correct for ectopic beats and we get a new RMSSD range of 6.5-167, which is still very high. My understanding is the correct_rr function corrects only peaks labelled as ectopic. Is it possible to correct for very short and very long intervals also?

Thanks in advance and best wishes

Leah

I encountered the following problem while analyzing PPG data using the cardioception package and systole for a Heartbeat counting task. While analyzing the data, I came across the following error and am posting this to inquire about a solution.

signal, peaks = ppg_peaks(ppg['2'][0], clean_extra=True, sfreq=25, new_sfreq=1000)

When I run this line, the following error occurs:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[43], line 6
      2 for nTrial in range(6):
      4     print(f'Analyzing trial number {nTrial+1}')
----> 6     signal, peaks = ppg_peaks(ppg[str(nTrial)][0], clean_extra=True, sfreq=25, new_sfreq=1000)
      7     axs = plot_raw(
      8         signal=signal, sfreq=1000, figsize=(18, 5), clean_extra=True,
      9         show_heart_rate=True
     10         );
     12     # Show the windows of interest
     13     # We need to convert sample vector into Matplotlib internal representation
     14     # so we can index it easily

File [c:\Users\user\anaconda3\envs\cardioception\lib\site-packages\systole\detection.py:133](file:///C:/Users/user/anaconda3/envs/cardioception/lib/site-packages/systole/detection.py:133), in ppg_peaks(signal, sfreq, win, new_sfreq, clipping, clipping_thresholds, moving_average, moving_average_length, peak_enhancement, distance, clean_extra, clean_nan, verbose)
    131     time = np.arange(0, len(x) / sfreq, 1 / sfreq)
    132     new_time = np.arange(0, len(x) / sfreq, 1 / new_sfreq)
--> 133     x = np.interp(new_time, time, x)
    135 # Copy resampled signal for output
    136 resampled_signal = np.copy(x)

File <__array_function__ internals>:180, in interp(*args, **kwargs)

File [c:\Users\user\anaconda3\envs\cardioception\lib\site-packages\numpy\lib\function_base.py:1570](file:///C:/Users/user/anaconda3/envs/cardioception/lib/site-packages/numpy/lib/function_base.py:1570), in interp(x, xp, fp, left, right, period)
   1567     xp = np.concatenate((xp[-1:]-period, xp, xp[0:1]+period))
   1568     fp = np.concatenate((fp[-1:], fp, fp[0:1]))
-> 1570 return interp_func(x, xp, fp, left, right)

ValueError: fp and xp are not of the same length.

To debug the issue, I used np.arange and found that the length of time and x do not match.

x = ppg['2'][0]
time = np.arange(0, len(x)/25, 0.04)
new_time = np.arange(0, len(x) / 25, 1/1000)
# x = np.interp(new_time, time, x)
print(len(time), len(new_time))
print(len(x))

947 37840
946

It doesn't match. maybe /25 and * 25 makes some error due to 946/25*25=946.0000000000001

Is there a good way to solve this issue? I appreciate any suggestions.

Hi Nico

If possible, it would be great to be able to double check the corrected peaks (saved in the json) after correcting peaks using the systole viewer.

Thanks in advance

Leah :)

Hi Nico

As requested, I am attaching the script and subject in which I get an interpolation error (numpy function interp) when using the ppg_peaks systole function. I have plotted the data and there are no obvious problems (also attached).

systole_sub-0550.zip

Thanks in advance,

Leah :)

Hi,

I thought it could be useful if the "to_rr" function could return a time series instead of a list of intervals.
I have made the following changes to the function:

def to_rr(peaks: Union[List[float], np.ndarray], sfreq: int = 1000, time_s: bool = False) -> np.ndarray:
    """Convert peaks index to intervals time series (RR, beat-to-beat...).
    Parameters
    ----------
    peaks : np.ndarray or list
        Either a boolean array or sample index. Default is *boolean*. If the
        input array does not only contain 0 or 1, will automatically try sample
        index.
    sfreq : int
        The sampling frequency (default is 1000 Hz).
    time_s: bool
        Whether to return an RR time series in that contains the distance to
        the next beat in each time point.
    Returns
    -------
    rr : np.ndarray
        Interval time series (in miliseconds).
    """
    if isinstance(peaks, list):
        peaks = np.asarray(peaks)
    if ((peaks == 1) | (peaks == 0)).all():
        rr = (np.diff(np.where(peaks)[0]) / sfreq) * 1000
    else:
        rr = (np.diff(peaks) / sfreq) * 1000
    if time_s:
        if ((peaks == 1) | (peaks == 0)).all():
            rr_time_series = np.zeros_like(peaks, dtype=float)
            non_zero_inds = np.where(peaks)[0]
            rr_time_series[non_zero_inds[:-1]] = rr
            # now replace zeros with last non-zero value
            last_non_zero = np.arange(len(rr_time_series))
            last_non_zero[rr_time_series == 0] = 0
            last_non_zero = np.maximum.accumulate(last_non_zero)
            rr_time_series = rr_time_series[last_non_zero]
            rr_time_series[:non_zero_inds[0]] = rr[0] # fill the first zero elements
            return rr_time_series
        else:
            raise ValueError('To get a RR time peaks must be a boolean array.')
    else:
        return rr

let me know if you want me to prepare a pulll request with this change.

Best,
Gidon

The function correction.correct_rr() first corrects artefacts labelled as missed or extra which change the number of RR intervals. This can cause the artefact corrections for the remaining types to fail in interpolate_bads() when idx can exceed len(rr), causing interp1d() to raise an exception. The reason for this is that interpolate_bads() is fed with cleaned_rr which may or may not correspond in length to the original rr for which the artefact detection was done.

I might make a pull request later. A quick solution is to change the order of the corrections, such that the corrections using interpolation come first.

Another question that might not be worth of an own issue, why is rr_artefacts(clean_rr) called iteratively after each correction bracket? I don't think that is the point in the original paper, as the committed corrections to clean_rr will affect the remaining detections, as the rolling thresholds are changed. The correction type annotations and the number of corrections will obviously differ running the artefact detection and the correction routine. Removing these and relying on the single pass artefact detection could cause inconsistent indexing for the short and long corrections, as they are done sequentially.