NeuroKit
NeuroKit copied to clipboard
ecg_process failed while simulating simple ecg signal
See the below example where I experienced a bug, and whatever the sampling_rate is:
# Simulate data: heart beat = 75 and 100 bpm
duration = 600; sampling_rate = 128; heart_rate = 75; noise = 0.0; method='simple'
ecg75 = nk.ecg_simulate(duration=duration,
sampling_rate=sampling_rate,
noise=noise,
heart_rate=heart_rate,
method=method,
random_state=10)
_, info = nk.ecg_process(ecg_signal=ecg75, sampling_rate=sampling_rate, method="neurokit")
Describe the bug a bug is raised in ecg_process, division by zero in find_artifacts. See the stack: C:\ProgramData\Anaconda3\envs\hop-child-covid\lib\site-packages\neurokit2\ecg\ecg_process.py:97: in ecg_process instant_peaks, rpeaks, = ecg_peaks( C:\ProgramData\Anaconda3\envs\hop-child-covid\lib\site-packages\neurokit2\ecg\ecg_peaks.py:74: in ecg_peaks _, rpeaks = signal_fixpeaks(rpeaks, sampling_rate=sampling_rate, iterative=True, method="Kubios") C:\ProgramData\Anaconda3\envs\hop-child-covid\lib\site-packages\neurokit2\signal\signal_fixpeaks.py:116: in signal_fixpeaks return _signal_fixpeaks_kubios(peaks, sampling_rate=sampling_rate, iterative=iterative, show=show) C:\ProgramData\Anaconda3\envs\hop-child-covid\lib\site-packages\neurokit2\signal\signal_fixpeaks.py:154: in _signal_fixpeaks_kubios artifacts, subspaces = _find_artifacts(peaks, sampling_rate=sampling_rate)
peaks = array([ 112, 214, 317, 419, 522, 624, 726, 829, 931, 1034, 1136, 1238, 1341, 1443, 1546,...75276, 75378, 75481, 75583, 75686, 75788, 75890, 75993, 76095, 76198, 76300, 76402, 76505, 76607, 76710]) c1 = 0.13, c2 = 0.17, alpha = 5.2, window_width = 91, medfilt_order = 11 sampling_rate = 128
def _find_artifacts(peaks, c1=0.13, c2=0.17, alpha=5.2, window_width=91, medfilt_order=11, sampling_rate=1000):
# Compute period series (make sure it has same numer of elements as peaks);
# peaks are in samples, convert to seconds.
rr = np.ediff1d(peaks, to_begin=0) / sampling_rate
# For subsequent analysis it is important that the first element has
# a value in a realistic range (e.g., for median filtering).
rr[0] = np.mean(rr[1:])
# Artifact identification #################################################
###########################################################################
# Compute dRRs: time series of differences of consecutive periods (dRRs).
drrs = np.ediff1d(rr, to_begin=0)
drrs[0] = np.mean(drrs[1:])
# Normalize by threshold.
th1 = _compute_threshold(drrs, alpha, window_width)
drrs /= th1
E RuntimeWarning: divide by zero encountered in true_divide
To Reproduce Just run the code above.
System Specifications Windows 10 Pro Conda env Python 3.8.6 neurokit2 0.1.2
Thanks for your help.
Hi there @FrancisThibaultNRC! Do you mean the error occurred when you ran this line:
_, info = nk.ecg_process(ecg_signal=ecg75, sampling_rate=sampling_rate, method="neurokit")
?
Not sure why but it works fine for me so I can't replicate the error 🤔 there's a line of code in _find_artifacts
to ignore the division by zero warning: np.seterr(divide="ignore", invalid="ignore")
. Can you let us know what's your numpy version?
Thanks for your follow-up.
From the stack, I have got this message while running ecg_process. The issue happened into the _find_artifacts function. I have numpy version 1.19.2. By ignoring the division by zero (by setting np.seterr(...)), this will hide an division by zero, which is an issue in _find_artifacts.
Thanks for clarifying. Hmm I'm not too sure about this, tagging @TiagoTostas here to help assist! Is ignoring division by zero appropriate in this case?
Hi all, I am still having this issue in _find_artifacts function, so a division by zero drrs /= th1. If I understand, the line np.seterr(divide='ignore', invalid='ignore') is set just right after the division by zero, and will not be effective in that case, isn't?
Hi @FrancisThibaultNRC
Sorry for this very delayed reply. And yes, the np.seterr()
function should be placed before and not after the division by zero. I will make a quick fix for this and it should be fixed in the next release.
This issue has been automatically marked as inactive because it has not had recent activity. It will eventually be closed if no further activity occurs.
Is this still relevant? If so, what is blocking it? Is there anything you can do to help move it forward?
This issue has been automatically marked as stale because it has not had recent activity. It will be closed if no further activity occurs.
Should be fixed now