Enhancement/gating interpolation

[rspwarnaar]

Summary

Improvements in the gating interpolation algorithms and addition of a quadratic interpolation method.

Related Issues

N/A

Changes Introduced

List the major changes made in this PR:

• Interpolation with numpy rolling mean in Gating method 1 and 3 yields ~ 83% and 96% reduction in evaluation time.
• Implemented gating method 4: quadratic fit to the surrounding data

Motivation and Context

Gating was one of the most computationally heavy procedures in a sEMG processing pipeline. The updated interpolation largely resolves this issue.
Quadradatic interpolation might offer more realistic filling of the gates than linear interpolation

Testing

Describe how you tested your changes. Include:

• Unit tests added or updated
• Benchmarking was performed against the original implementation to test similarity of results and processing time. Method 3 does not yield exact replication of original results in relation to pandas rolling average, but error is negligible (< 3.552713678801e-15 for max difference, 1.426964945016e-17 mean error).

======================================================================
Method 1 comparison
======================================================================
Original ReSurfEMG v1.1.0 (mean ± std): 0.043134 ± 0.006805 s
New method (mean ± std): 0.007436 ± 0.001630 s
Average speed improvement: 82.681% (mean over 100 reps)
Average time per peak (v1.1.0): 0.00015350 s/peak
Average time per peak (new):     0.00002646 s/peak
---------------------------------------------------------------------
Runs with exact zero max-difference: 100/100
Runs with max-difference ≈ 0 (tol=1e-12): 100/100
Max of max-differences across runs: 0.000000000000e+00
Mean of max-differences across runs: 0.000000000000e+00
Mean of mean-differences across runs: 0.000000000000e+00
======================================================================
Method 2 comparison
======================================================================
Original ReSurfEMG v1.1.0 (mean ± std): 0.026278 ± 0.009667 s
New method (mean ± std): 0.025894 ± 0.006778 s
Average speed improvement: -2.686% (mean over 100 reps)
Average time per peak (v1.1.0): 0.00009351 s/peak
Average time per peak (new):     0.00009215 s/peak
---------------------------------------------------------------------
Runs with exact zero max-difference: 100/100
Runs with max-difference ≈ 0 (tol=1e-12): 100/100
Max of max-differences across runs: 0.000000000000e+00
Mean of max-differences across runs: 0.000000000000e+00
Mean of mean-differences across runs: 0.000000000000e+00
======================================================================
Method 3 comparison
======================================================================
Original ReSurfEMG v1.1.0 (mean ± std): 3.198082 ± 0.887588 s
New method (mean ± std): 0.106740 ± 0.025473 s
Average speed improvement: 96.587% (mean over 5 reps)
Average time per peak (v1.1.0): 0.01138108 s/peak
Average time per peak (new):     0.00037986 s/peak
---------------------------------------------------------------------
Runs with exact zero max-difference: 0/5
Runs with max-difference ≈ 0 (tol=1e-12): 5/5
Max of max-differences across runs: 3.552713678801e-15
Mean of max-differences across runs: 3.552713678801e-15
Mean of mean-differences across runs: 1.426964945016e-17
======================================================================

Checklist

• I have read the contributing guidelines
• My code follows the project's coding style
• I have added tests that prove my fix is effective or my feature works
• I have updated the documentation accordingly
• All existing and new tests pass

Additional Notes

Any other context or information reviewers should know?

[codecov]

Codecov Report

❌ Patch coverage is 73.21429% with 15 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
resurfemg/preprocessing/ecg_removal.py	73.21%	15 Missing ⚠️

📢 Thoughts on this report? Let us know!

[ChiaraM96]

Using a cumsum-based implementation instead of the concatenate/for loop would significantly cut execution time here (I did some benchmarks and the difference is about ~10x)

[ChiaraM96]

Suggested change

	starts = np.maximum(0, peaks - half_gate_width)
	ends = np.minimum(max_sample, peaks + half_gate_width)
	# build ranges per peak and concat, then unique to avoid duplicates
	gate_samples = np.concatenate(
	[np.arange(s, e) for s, e in zip(starts, ends)])
	gate_samples = np.unique(gate_samples).astype(int)
	starts= np.clip(peaks - half_gate_width, 0, max_sample).astype(int)
	ends= np.clip(peaks + half_gate_width, 0, max_sample).astype(int)
	delta = np.zeros(max_sample + 1, dtype=np.int64)
	np.add.at(delta, starts, 1) # +1 entering each gate
	np.add.at(delta, ends, -1) # -1 leaving each gate
	gate_samples= np.cumsum(delta)[:max_sample] > 0

[rspwarnaar]

I tested the vectorized method versus the original new version of method 0 (zeros), and it gave me a quite different result (test dataset, 860 k samples, 281 peaks): the vectorized method would be over twice as slow as compared to the original new method. Any idea where this difference might come from? I'll push the notebook with both versions for comparison

======================================================================
Method 0 comparison - New vs Vectorized
======================================================================
New method (mean ± std): 0.007217 ± 0.001819 s
New method vectorized (mean ± std): 0.018283 ± 0.003342 s
Average speed improvement: -167.662% (mean over 100 reps)
Average time per peak (new): 0.00002568 s/peak
Average time per peak (new vectorized):     0.00006506 s/peak
---------------------------------------------------------------------
Runs with exact zero max-difference: 100/100
Runs with max-difference ≈ 0 (tol=1e-12): 100/100
Max of max-differences across runs: 0.000000000000e+00
Mean of max-differences across runs: 0.000000000000e+00
Mean of mean-differences across runs: 0.000000000000e+00
======================================================================
`

[ChiaraM96]

I pushed a new version of the notebook including some changes in the gating_new_vectorized algorithm and a profiling/comparison method. Here's the changes in speed I get now:

======================================================================
Method 0 comparison - Vectorized
======================================================================
Original ReSurfEMG v1.1.0 (mean ± std): 0.006804 ± 0.001617 s
New vectorized method (mean ± std): 0.002189 ± 0.000361 s
Average speed improvement: 67.171% (mean over 100 reps)
Average time per peak (v1.1.0): 0.00002421 s/peak
Average time per peak (new):     0.00000779 s/peak
---------------------------------------------------------------------
Runs with exact zero max-difference: 100/100
Runs with max-difference ≈ 0 (tol=1e-12): 100/100
Max of max-differences across runs: 0.000000000000e+00
Mean of max-differences across runs: 0.000000000000e+00
Mean of mean-differences across runs: 0.000000000000e+00
======================================================================

Aside from the vectorization, I think the largest impact on speed in this version is due to me switching from
emg_raw_gated = copy.deepcopy(emg_raw)
to
emg_raw_gated = np.array(emg_raw, dtype=float, copy=True)
as on np arrays, deepcopy is only necessary if the array's dtype is object

[ChiaraM96]

in the if fit_available branch of the quadratic gate filling method, the polyfit and polyval functions are passed the absolute valid_idx indexes; for long recordings, these indexes become quite large, potentially resulting in an ill-conditioned matrix for the solver and in approximation errors. Using offsetted support times (running the polyfit on valid_idxs - gs, and then polyval on ks - gs) largely reduces the magnitude of the x values used in the calculations, improving the stability of the algorithm.