OSCR

GROND: automated characterization of periodic and rhythmic EEG patterns.

Code ↔ Paper

40 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 40 matches
  1. [1] § Methods › RDA characterization pipeline (RDA-Profiler) › NB-Hilbert: iterative NB-Hilbert refinement ↔ code/evaluation/verify_v12_holdout.py, lines 1 to 47 · score 0.96 · canonical cohort LRDA, freq_cap, top_k, p2 bw, p1 hi, NB Hilbert · Report a mismatch
  2. [2] § Results › Frequency estimation ↔ code/evaluation/lrda_laterality_v13_eval.py, lines 1 to 66 · score 0.96 · asymmetric volume conduction, narrowband envelope amplitude, amplitude normalized rhythmicity, spectral peak prominence, slow drift, rhythm frequency · Report a mismatch
  3. [3] § Methods › PD characterization pipeline (PD-Profiler) › Spatial localization ↔ paper_materials/generate_discharge_topo_viewer.py, lines 136 to 247 · score 0.95 · strongly suppresses phantom, GFP aligned peak, GFP profile, low GFP, Discharge locked topographic, cross correlates · Report a mismatch
  4. [4] § Methods › Frequency ranges and filter bands ↔ code/pd_profiler.py, lines 1 to 42 · score 0.95 · expectation maximization, template refinement, channel PD probabilities, characterization pipeline, dynamic programming, PD Profiler · Report a mismatch
  5. [5] § Methods › Dataset and annotations › Annotation framework ↔ paper_materials/build_fig_irr_bars.py, lines 1 to 51 · score 0.92 · interactive narrowband overlay, Open circles, error bars, rater independent expert, frequency Spearman, paired segment · Report a mismatch
  6. [6] § Methods › PD characterization pipeline (PD-Profiler) › channel CNN (ChannelPD-Net) ↔ code/pd_channel_detector/channel_cnn.py, lines 92 to 165 · score 0.91 · temporal attention pooling, PD detection, attention weighted, log frequency, PD probability, softmax · Report a mismatch
  7. [7] § Methods › Frequency ranges and filter bands ↔ paper_materials/draw_panel_b_rda.py, lines 58 to 187 · score 0.86 · iterative narrowband refinement, Pipeline architecture, refined lateralization, envelope amplitude, NB Hilbert, dominant channels · Report a mismatch
  8. [8] § Results › Model architecture comparison ↔ paper_materials/tables/generate_table7.py, lines 40 to 140 · score 0.85 · end neural models, E2E CNN, DP encoded, neural evidence generation, winning strategy, structured inference · Report a mismatch
  9. [9] § Methods › Dataset and annotations › Annotation framework ↔ code/evaluation/lrda_laterality_v13_eval.py, lines 1 to 66 · score 0.84 · amplitude normalized rhythmicity, spectral peak prominence, peak amplitude, laterality rule, hemisphere channels, independent expert · Report a mismatch
  10. [10] § Methods › PD characterization pipeline (PD-Profiler) › channel CNN (ChannelPD-Net) ↔ code/unified_model/dataset.py, lines 215 to 316 · score 0.84 · 20–40 dB, Gaussian noise, log frequency, amplitude scaling, circular, SNR · Report a mismatch
  11. [11] § Methods › Dataset and annotations › Evaluation framework ↔ paper_materials/build_fig_irr_bars.py, lines 1 to 51 · score 0.81 · interactive narrowband overlay, IRR comparisons, rater independent expert, EA, LB, tools · Report a mismatch
  12. [12] § Methods › RDA characterization pipeline (RDA-Profiler) › NB-Hilbert: iterative NB-Hilbert refinement ↔ code/generators/labeling/generate_rda_freq_labeler.py, lines 148 to 202 · score 0.80 · 0.3–5 Hz, Coarse lateralization, 0.5–3.5 Hz, dominant hemisphere, variance, W05 · Report a mismatch
  13. [13] § Methods › RDA characterization pipeline (RDA-Profiler) › NB-Hilbert: iterative NB-Hilbert refinement ↔ code/generators/labeling/generate_rda_freq_spatial_labeler.py, lines 148 to 202 · score 0.80 · 0.3–5 Hz, Coarse lateralization, 0.5–3.5 Hz, dominant hemisphere, variance, W05 · Report a mismatch
  14. [14] § Methods › PD characterization pipeline (PD-Profiler) › channel CNN (ChannelPD-Net) ↔ code/cet_model/lrda_crnn_dataset.py, lines 105 to 167 · score 0.80 · additive Gaussian noise, 20–40 dB, amplitude scaling, SNR, augmentation, shifts · Report a mismatch
  15. [15] § Methods › Dataset and annotations › Annotation framework ↔ code/evaluation/analyze_independent_expert_v1.py, lines 1 to 65 · score 0.78 · confidence interval, independent_expert_tasks, paper_materials, expert algorithm, v1, bootstrap · Report a mismatch
  16. [16] § Methods › Frequency ranges and filter bands ↔ paper_materials/draw_panel_b.py, lines 45 to 122 · score 0.77 · channel PD probabilities, UNet, PD Profiler, bipolar channels, alignment, template · Report a mismatch
  17. [17] § Results › Spatial localization › PD spatial localization ↔ paper_materials/generate_discharge_topo_viewer.py, lines 136 to 247 · score 0.77 · Laplacian GFP alignment, Laplacian transform, Discharge locked topographic, phantom, focal, squared · Report a mismatch
  18. [18] § Methods › Dataset and annotations › Evaluation framework ↔ code/evaluation/analyze_phlbsz_cohort.py, lines 1 to 37 · score 0.76 · narrowband overlay tools, validation cohort, tools developed, EA, LB, PH · Report a mismatch
  19. [19] § Methods › RDA characterization pipeline (RDA-Profiler) › RDA-PLV: spatial extent ↔ code/rda_spatial_extent.py, lines 63 to 142 · score 0.76 · PLV amplitude, phase locking, dominant hemisphere, narrowband filtered, envelope, NB · Report a mismatch
  20. [20] § Methods › Dataset and annotations › Contest-of-methods development strategy ↔ code/evaluation/lrda_freq_hyperparam_sweep.py, lines 1 to 40 · score 0.74 · channel selection, NB Hilbert, RDA frequency, independent expert, variants, contest · Report a mismatch
  21. [21] § Methods › PD characterization pipeline (PD-Profiler) › Discharge timing: HemiCET-UNet and dynamic programming ↔ code/cet_model/train_cet_v3.py, lines 270 to 395 · score 0.74 · sharpness penalty, fold patient stratified, Training targets, augmentation, loss, Gaussian · Report a mismatch
  22. [22] § Methods › Dataset and annotations › Evaluation framework ↔ code/data_management/build_grond_h5_bank.py, lines 1 to 83 · score 0.74 · cohort flags, paper_materials, algorithm predictions, grond, bank, EEG segments · Report a mismatch
  23. [23] § Methods › PD characterization pipeline (PD-Profiler) › channel CNN (ChannelPD-Net) ↔ code/unified_model/model.py, lines 19 to 158 · score 0.74 · temporal attention pooling, softmax, kernel, sigmoid, logit, stride · Report a mismatch
  24. [24] § Methods › PD characterization pipeline (PD-Profiler) › Discharge timing: HemiCET-UNet and dynamic programming ↔ code/discharge_detector.py, lines 1 to 65 · score 0.72 · approximately periodic prior, active interval, Candidate peaks, skip, DP, trace · Report a mismatch
  25. [25] § Methods › PD characterization pipeline (PD-Profiler) › Discharge timing: HemiCET-UNet and dynamic programming ↔ code/discharge_detector.py, lines 764 to 845 · score 0.72 · post hoc confidence, template refinement, refined, median, candidates, DP · Report a mismatch
  26. [26] § Methods › PD characterization pipeline (PD-Profiler) › Discharge timing: HemiCET-UNet and dynamic programming ↔ code/cet_model/cet.py, lines 95 to 235 · score 0.72 · skip connections, discharge evidence, upsamples, activated, resolution, encoder · Report a mismatch
  27. [27] § Methods › PD characterization pipeline (PD-Profiler) ↔ code/pd_profiler.py, lines 1 to 42 · score 0.72 · expectation maximization, PD Profiler, dynamic programming, periodic discharge, EM, TKEO · Report a mismatch
  28. [28] § Methods › BIPD detection ↔ code/bipd_detector.py, lines 896 to 964 · score 0.71 · Synthetic BIPD, phase shifting, single LPD, GBT, jitter, candidates · Report a mismatch
  29. [29] § Methods › Frequency ranges and filter bands ↔ code/generators/labeling/generate_rda_freq_spatial_labeler.py, lines 148 to 202 · score 0.71 · 0.5–3.5 Hz, refined lateralization, Hilbert frequency, dominant hemisphere, coarse, variance · Report a mismatch
  30. [30] § Methods › PD characterization pipeline (PD-Profiler) › Discharge timing: HemiCET-UNet and dynamic programming ↔ code/discharge_detector.py, lines 160 to 245 · score 0.70 · skip connections, discharge evidence, transposed, resolution, encoder, sigmoid · Report a mismatch
  31. [31] § Methods › PD characterization pipeline (PD-Profiler) › Discharge timing: HemiCET-UNet and dynamic programming ↔ code/hemi_detector/train_hemi_cet.py, lines 356 to 502 · score 0.69 · sharpness penalty, fold patient stratified, floor, augmentation, loss, Gaussian · Report a mismatch
  32. [32] § Methods › Statistical analysis ↔ code/evaluation/eval_3way_classification.py, lines 212 to 253 · score 0.68 · random forest, fold patient stratified, fold prediction, RF, classifier, training · Report a mismatch
  33. [33] § Results › Frequency estimation ↔ code/evaluation/lrda_laterality_features.py, lines 88 to 172 · score 0.68 · spectral peak prominence, envelope amplitude, LRDA laterality, artifact, drift, rhythm · Report a mismatch
  34. [34] § Methods › Dataset and annotations › EEG data ↔ paper_materials/tables/generate_table1.py, lines 141 to 232 · score 0.68 · expert review, expert votes, reference montage, single rater, channel monopolar, EEG segments · Report a mismatch
  35. [35] § Methods › Dataset and annotations › Contest-of-methods development strategy ↔ paper_materials/draw_panel_b_rda.py, lines 58 to 187 · score 0.67 · iterative narrowband refinement, NB Hilbert, coherence, rounds, W05, refined · Report a mismatch
  36. [36] § Methods › Frequency ranges and filter bands ↔ code/evaluation/lrda_freq_hyperparam_sweep.py, lines 56 to 97 · score 0.66 · 0.3–5 Hz, pre filter, half width, Hilbert, clipped, 0.25 Hz · Report a mismatch
  37. [37] § Results › Spatial localization › PD spatial localization ↔ paper_materials/tables/generate_table4.py, lines 33 to 111 · score 0.66 · model expert, inter rater agreement, EE pairs, Jaccard, S2, LB · Report a mismatch
  38. [38] § Methods › PD characterization pipeline (PD-Profiler) › Discharge timing: HemiCET-UNet and dynamic programming ↔ code/evaluation/verify_c1_repro.py, lines 87 to 204 · score 0.66 · detected active interval, skip penalties, deviations, Candidate, DP, HemiCET · Report a mismatch
  39. [39] § Results › Spatial localization › PD spatial localization ↔ paper_materials/generate_fig0_examples.py, lines 238 to 365 · score 0.64 · Bottom row, rhythmic delta activity, lateralized rhythmic delta, generalized rhythmic delta, lateralized periodic discharges, Inter rater · Report a mismatch
  40. [40] § Methods › PD characterization pipeline (PD-Profiler) › Discharge timing: HemiCET-UNet and dynamic programming ↔ code/discharge_detector.py, lines 1 to 65 · score 0.63 · Product boosted, discharge detection, floor, aggregated, CET, HPP · Report a mismatch

Paper

Loaded from Europe PMC by your browser; OSCR does not store it: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 879 lines · 35 KB · CC-BY-NC-4.0 · 4 matches

  1. """
  2. Automated Periodic Discharge Detection from Raw EEG
  4. Self-contained implementation of the best-performing method for detecting
  5. periodic discharges (LPD/GPD) in 10-second EEG segments. No gold standard
  6. labels are used as input — the system operates from raw EEG only.
  8. Method: product-boosted max(HPP, CET-UNet) + CNN+ACF frequency
  9. + optimized DP + post-hoc filtering
  10. Performance: F1≈0.74, evaluated on 593 patients with expert-reviewed labels
  12. Architecture overview:
  13. EEG (18ch × 2000 samples @ 200 Hz)
  14. │
  15. ├─ CNN+Attention (5-fold ensemble)
  16. │ → per-channel PD probability + log-frequency
  17. │ → PD-weighted CNN frequency estimate
  18. │
  19. ├─ ACF on pointiness traces (per channel)
  20. │ → ACF frequency estimate
  21. │ → Ensemble: 0.8×CNN + 0.2×ACF
  22. │
  23. ├─ Handcrafted evidence (per channel):
  24. │ pointiness trace (prominence²/width at peaks)
  25. │ + TKEO on 20 Hz lowpassed signal
  26. │ → weighted combination → Gaussian smooth → E_hpp(t)
  27. │
  28. └─ CET-UNet (5-fold ensemble, per channel):
  29. → frame-level discharge evidence E_cet(t)
  30. │
  31. ├─ Aggregate channels (median, laterality-aware)
  32. └─ Combine: threshold CET (80th pct), floor (0.3),
  33. product-boost: max(HPP,CET) + 3×HPP×CET
  34. │
  35. └─ Active interval detection
  36. → candidate peak extraction (min height 5% of max)
  37. → forward DP with approximately-periodic prior
  38. (α=1.275, β=0.3, λ=0.05, max_skip=3)
  39. → EM template refinement
  40. → post-hoc filter (drop peaks < 0.3× median)
  41. → per-channel timing (±50ms)
  42. │
  43. Output: discharge times, IPI-derived frequency
  45. Dependencies: numpy, scipy, torch
  46. Model weights: data/cet_cache/cet_unet_fold{0-4}.pt
  47. data/pd_channel_cache/cnn_attn_fold{0-4}.pt
  49. Usage:
  50. from discharge_detector import DischargeDetector
  51. detector = DischargeDetector()
  52. result = detector.detect(eeg_18ch, subtype='lpd', laterality='left')
  53. # result['global_times'] = [0.42, 1.18, 1.95, 2.71, ...]
  54. # result['frequency'] = 1.31 (Hz, from IPI)
  55. # result['freq_estimate_input'] = 1.28 (Hz, CNN estimate)
  57. Authors: Westover Lab, Massachusetts General Hospital
  58. """
  60. import numpy as np
  61. import torch
  62. import torch.nn as nn
  63. from pathlib import Path
  64. from scipy.signal import find_peaks, butter, filtfilt
  65. from scipy.ndimage import gaussian_filter1d
  67. # ── Path setup ────────────────────────────────────────────────────────
  68. CODE_DIR = Path(__file__).resolve().parent
  69. PROJECT_DIR = CODE_DIR.parent
  71. # ── Constants ─────────────────────────────────────────────────────────
  72. FS = 200 # sampling rate (Hz)
  74. # Channel indices in 18-channel bipolar montage (double banana)
  75. # 0-3: Fp1-F7, F7-T3, T3-T5, T5-O1 (left temporal chain)
  76. # 4-7: Fp2-F8, F8-T4, T4-T6, T6-O2 (right temporal chain)
  77. # 8-11: Fp1-F3, F3-C3, C3-P3, P3-O1 (left parasagittal chain)
  78. # 12-15: Fp2-F4, F4-C4, C4-P4, P4-O2 (right parasagittal chain)
  79. # 16-17: Fz-Cz, Cz-Pz (midline)
  80. LEFT_INDICES = np.array([0, 1, 2, 3, 8, 9, 10, 11])
  81. RIGHT_INDICES = np.array([4, 5, 6, 7, 12, 13, 14, 15])
  83. # Handcrafted evidence parameters
  84. LOWPASS_HZ = 20.0
  85. POINTINESS_WEIGHT = 0.6
  86. TKEO_WEIGHT = 0.4
  87. SMOOTH_SIGMA_SAMPLES = 3 # ~15 ms at 200 Hz
  89. # Active interval detection
  90. ROLLING_WINDOW_S = 1.0
  91. ACTIVE_THRESHOLD_FRAC = 0.5
  92. MIN_ACTIVE_SECONDS = 3.0
  93. ACTIVE_EXPAND_S = 0.5
  95. # DP parameters (optimized via parameter sweep)
  96. DP_ALPHA = 1.275 # timing deviation penalty (loosened for noisy freq)
  97. DP_BETA = 0.3 # skip penalty
  98. DP_LAMBDA = 0.05 # new-sequence cost
  99. PEAK_HEIGHT_FRAC = 0.05
  100. MAX_SKIP = 3
  102. # EM refinement
  103. TEMPLATE_HALF_MS = 150 # ±150 ms for template extraction
  104. CHANNEL_REFINE_MS = 50 # ±50 ms for per-channel refinement
  106. # Evaluation
  107. TOLERANCE_S = 0.1 # ±100ms for discharge matching
  110. # ═══════════════════════════════════════════════════════════════════════
  111. # Neural network architectures
  112. # ═══════════════════════════════════════════════════════════════════════
  114. class ChannelPDNetAttention(nn.Module):
  115. """1D CNN with temporal attention for PD detection + frequency estimation.
  117. Input: (batch, 1, 2000) — one z-scored EEG channel
  118. Output: pd_prob (batch, 1) — PD probability [0, 1]
  119. freq_pred (batch, 1) — log frequency (Hz)
  120. attn_wts (batch, 1, T) — attention weights over time
  122. Architecture:
  123. 4 conv blocks (1→16→32→64→64), each: Conv1d(stride=2) → BN → GELU → Dropout
  124. Attention: Conv1d(64→1) → softmax → weighted pooling → (64,)
  125. PD head: Linear(64→1) → sigmoid
  126. Freq head: Linear(64→1)
  127. """
  129. def __init__(self):
  130. super().__init__()
  131. self.block1 = nn.Sequential(
  132. nn.Conv1d(1, 16, kernel_size=51, stride=2, padding=25),
  133. nn.BatchNorm1d(16), nn.GELU(), nn.Dropout(0.1))
  134. self.block2 = nn.Sequential(
  135. nn.Conv1d(16, 32, kernel_size=25, stride=2, padding=12),
  136. nn.BatchNorm1d(32), nn.GELU(), nn.Dropout(0.1))
  137. self.block3 = nn.Sequential(
  138. nn.Conv1d(32, 64, kernel_size=13, stride=2, padding=6),
  139. nn.BatchNorm1d(64), nn.GELU(), nn.Dropout(0.1))
  140. self.block4 = nn.Sequential(
  141. nn.Conv1d(64, 64, kernel_size=7, stride=2, padding=3),
  142. nn.BatchNorm1d(64), nn.GELU(), nn.Dropout(0.2))
  143. self.attn_conv = nn.Conv1d(64, 1, kernel_size=1)
  144. self.pd_head = nn.Linear(64, 1)
  145. self.freq_head = nn.Linear(64, 1)
  147. def forward(self, x):
  148. x = self.block1(x)
  149. x = self.block2(x)
  150. x = self.block3(x)
  151. x = self.block4(x)
  152. attn_logits = self.attn_conv(x)
  153. attn_wts = torch.softmax(attn_logits, dim=-1)
  154. pooled = (x * attn_wts).sum(dim=-1)
  155. pd_prob = torch.sigmoid(self.pd_head(pooled))
  156. freq_pred = self.freq_head(pooled)
  157. return pd_prob, freq_pred, attn_wts
  160. class CETUNet(nn.Module):
  161. """U-Net for frame-level discharge evidence at full temporal resolution.
  163. Input: (batch, 1, 2000) — one z-scored EEG channel
  164. Output: (batch, 1, 2000) — discharge evidence in [0, 1]
  166. Architecture:
  167. Encoder: 4 conv blocks with stride=2 downsampling
  168. e1: 1→16 (2000→1000) k=51
  169. e2: 16→32 (1000→500) k=25
  170. e3: 32→64 (500→250) k=13
  171. e4: 64→64 (250→125) k=7
  173. Decoder: 4 ConvTranspose blocks with skip connections from encoder
  174. d4: 64→64 (125→250), cat(d4, e3) → Conv(128→64)
  175. d3: 64→32 (250→500), cat(d3, e2) → Conv(64→32)
  176. d2: 32→16 (500→1000), cat(d2, e1) → Conv(32→16)
  177. d1: 16→8 (1000→2000), Conv(8→1) → Sigmoid
  178. """
  180. def __init__(self):
  181. super().__init__()
  182. self.enc1 = nn.Sequential(
  183. nn.Conv1d(1, 16, kernel_size=51, stride=2, padding=25),
  184. nn.BatchNorm1d(16), nn.GELU())
  185. self.enc2 = nn.Sequential(
  186. nn.Conv1d(16, 32, kernel_size=25, stride=2, padding=12),
  187. nn.BatchNorm1d(32), nn.GELU())
  188. self.enc3 = nn.Sequential(
  189. nn.Conv1d(32, 64, kernel_size=13, stride=2, padding=6),
  190. nn.BatchNorm1d(64), nn.GELU())
  191. self.enc4 = nn.Sequential(
  192. nn.Conv1d(64, 64, kernel_size=7, stride=2, padding=3),
  193. nn.BatchNorm1d(64), nn.GELU())
  195. self.up4 = nn.Sequential(
  196. nn.ConvTranspose1d(64, 64, kernel_size=4, stride=2, padding=1),
  197. nn.BatchNorm1d(64), nn.GELU())
  198. self.skip4 = nn.Sequential(
  199. nn.Conv1d(128, 64, kernel_size=3, padding=1),
  200. nn.BatchNorm1d(64), nn.GELU())
  202. self.up3 = nn.Sequential(
  203. nn.ConvTranspose1d(64, 32, kernel_size=4, stride=2, padding=1),
  204. nn.BatchNorm1d(32), nn.GELU())
  205. self.skip3 = nn.Sequential(
  206. nn.Conv1d(64, 32, kernel_size=3, padding=1),
  207. nn.BatchNorm1d(32), nn.GELU())
  209. self.up2 = nn.Sequential(
  210. nn.ConvTranspose1d(32, 16, kernel_size=4, stride=2, padding=1),
  211. nn.BatchNorm1d(16), nn.GELU())
  212. self.skip2 = nn.Sequential(
  213. nn.Conv1d(32, 16, kernel_size=3, padding=1),
  214. nn.BatchNorm1d(16), nn.GELU())
  216. self.up1 = nn.Sequential(
  217. nn.ConvTranspose1d(16, 8, kernel_size=4, stride=2, padding=1),
  218. nn.BatchNorm1d(8), nn.GELU())
  219. self.head = nn.Sequential(nn.Conv1d(8, 1, kernel_size=1), nn.Sigmoid())
  221. @staticmethod
  222. def _match_size(dec, enc):
  223. min_len = min(dec.shape[2], enc.shape[2])
  224. return dec[:, :, :min_len], enc[:, :, :min_len]
  226. def forward(self, x):
  227. e1 = self.enc1(x)
  228. e2 = self.enc2(e1)
  229. e3 = self.enc3(e2)
  230. e4 = self.enc4(e3)
  232. d4 = self.up4(e4)
  233. d4, e3m = self._match_size(d4, e3)
  234. d4 = self.skip4(torch.cat([d4, e3m], dim=1))
  236. d3 = self.up3(d4)
  237. d3, e2m = self._match_size(d3, e2)
  238. d3 = self.skip3(torch.cat([d3, e2m], dim=1))
  240. d2 = self.up2(d3)
  241. d2, e1m = self._match_size(d2, e1)
  242. d2 = self.skip2(torch.cat([d2, e1m], dim=1))
  244. d1 = self.up1(d2)
  245. return self.head(d1)
  248. # ═══════════════════════════════════════════════════════════════════════
  249. # Signal processing: handcrafted evidence
  250. # ═══════════════════════════════════════════════════════════════════════
  252. def compute_pointiness_trace(signal_1d, half_win=8):
  253. """Compute pointiness = prominence² / width at each local maximum.
  255. For each local max in the signal, computes:
  256. prominence = peak_value - max(left_valley, right_valley)
  257. width = number of samples above half-prominence level
  258. pointiness = prominence² / width
  260. Returns a trace that is zero everywhere except at local maxima.
  261. """
  262. n = len(signal_1d)
  263. trace = np.zeros(n)
  264. peaks, _ = find_peaks(signal_1d)
  265. for loc in peaks:
  266. if loc < half_win or loc >= n - half_win:
  267. continue
  268. peak_val = signal_1d[loc]
  269. left_valley = np.min(signal_1d[loc - half_win:loc])
  270. right_valley = np.min(signal_1d[loc + 1:loc + half_win + 1])
  271. prom = peak_val - max(left_valley, right_valley)
  272. if prom <= 0:
  273. continue
  274. half_prom_level = peak_val - 0.5 * prom
  275. width = 0
  276. for j in range(1, half_win + 1):
  277. if signal_1d[loc - j] > half_prom_level:
  278. width += 1
  279. else:
  280. break
  281. for j in range(1, half_win + 1):
  282. if loc + j < n and signal_1d[loc + j] > half_prom_level:
  283. width += 1
  284. else:
  285. break
  286. if width > 0:
  287. trace[loc] = prom ** 2 / width
  288. return trace
  291. def compute_channel_evidence(signal_1d, fs=FS):
  292. """Compute handcrafted discharge evidence E_hpp(t) for one channel.
  294. Steps:
  295. 1. Pointiness trace on raw signal
  296. 2. TKEO (Teager-Kaiser Energy Operator) on 20 Hz lowpassed signal
  297. 3. Z-score normalize both within the window
  298. 4. Combine: 0.6 × pointiness_z + 0.4 × tkeo_z
  299. 5. Gaussian smooth (σ=3 samples ≈ 15ms)
  300. 6. Clip negatives to zero
  301. """
  302. n = len(signal_1d)
  303. if n < 10:
  304. return np.zeros(n)
  306. pt = compute_pointiness_trace(signal_1d)
  308. b_lp, a_lp = butter(4, LOWPASS_HZ / (fs / 2), btype='low')
  309. try:
  310. sig_lp = filtfilt(b_lp, a_lp, signal_1d)
  311. except ValueError:
  312. sig_lp = signal_1d.copy()
  314. tkeo = np.zeros(n)
  315. if n >= 3:
  316. tkeo[1:-1] = np.abs(sig_lp[1:-1] ** 2 - sig_lp[:-2] * sig_lp[2:])
  318. def zscore(x):
  319. m, s = np.mean(x), np.std(x)
  320. return (x - m) / s if s > 1e-10 else np.zeros_like(x)
  322. evidence = POINTINESS_WEIGHT * zscore(pt) + TKEO_WEIGHT * zscore(tkeo)
  323. evidence = gaussian_filter1d(evidence, sigma=SMOOTH_SIGMA_SAMPLES)
  324. return np.clip(evidence, 0, None)
  327. # ═══════════════════════════════════════════════════════════════════════
  328. # Evidence aggregation
  329. # ═══════════════════════════════════════════════════════════════════════
  331. def aggregate_evidence(evidence_all, subtype, laterality=None):
  332. """Aggregate per-channel evidence into a single trace E(t).
  334. GPD: median across all 18 channels
  335. LPD with known laterality: median across ipsilateral 8 channels
  336. LPD without laterality: max(left_median, right_median) at each time point
  337. """
  338. if subtype == 'gpd':
  339. return np.median(evidence_all, axis=0)
  340. if laterality == 'left':
  341. return np.median(evidence_all[LEFT_INDICES], axis=0)
  342. elif laterality == 'right':
  343. return np.median(evidence_all[RIGHT_INDICES], axis=0)
  344. else:
  345. left_med = np.median(evidence_all[LEFT_INDICES], axis=0)
  346. right_med = np.median(evidence_all[RIGHT_INDICES], axis=0)
  347. return np.maximum(left_med, right_med)
  350. def combine_evidence(hpp_evidence, cet_evidence, cet_threshold_pct=80,
  351. boost_weight=3.0, cet_floor=0.3):
  352. """Combine handcrafted and CNN evidence with product-boost.
  354. Steps:
  355. 1. Normalize both traces to [0, 1]
  356. 2. Threshold CET: zero out values below the 80th percentile
  357. (suppresses CET noise floor that causes false positives)
  358. 3. Product-boost: max(HPP, CET) + boost_weight × HPP × CET
  359. (amplifies regions where both methods agree)
  360. """
  361. hpp_max = np.max(hpp_evidence) if np.max(hpp_evidence) > 0 else 1.0
  362. cet_max = np.max(cet_evidence) if np.max(cet_evidence) > 0 else 1.0
  363. hpp_norm = hpp_evidence / hpp_max
  364. cet_norm = cet_evidence / cet_max
  366. # Threshold CET to suppress noise floor
  367. if cet_threshold_pct > 0 and np.any(cet_norm > 0):
  368. thr = np.percentile(cet_norm[cet_norm > 0], cet_threshold_pct)
  369. cet_norm = np.where(cet_norm > thr, cet_norm, 0)
  371. # Floor: zero out very low CET values
  372. cet_clean = np.where(cet_norm > cet_floor, cet_norm, 0)
  374. # Product-boost: max + agreement bonus
  375. base = np.maximum(hpp_norm, cet_clean)
  376. agreement = hpp_norm * cet_clean
  377. return base + boost_weight * agreement
  380. # ═══════════════════════════════════════════════════════════════════════
  381. # Dynamic programming for discharge timing
  382. # ═══════════════════════════════════════════════════════════════════════
  384. def detect_active_interval(evidence, fs=FS):
  385. """Find the longest contiguous interval where the evidence is active.
  387. Active = rolling mean (1s window) exceeds 50% of the global rolling max.
  388. Interval is expanded by 0.5s on each side.
  389. """
  390. n = len(evidence)
  391. win = int(ROLLING_WINDOW_S * fs)
  392. cs = np.cumsum(evidence)
  393. rolling = np.zeros(n)
  394. for i in range(n):
  395. lo = max(0, i - win // 2)
  396. hi = min(n, i + win // 2 + 1)
  397. rolling[i] = (cs[hi - 1] - (cs[lo - 1] if lo > 0 else 0)) / (hi - lo)
  398. threshold = ACTIVE_THRESHOLD_FRAC * np.max(rolling)
  399. above = rolling > threshold
  400. best_start, best_len = 0, 0
  401. cur_start, cur_len = 0, 0
  402. for i in range(n):
  403. if above[i]:
  404. if cur_len == 0:
  405. cur_start = i
  406. cur_len += 1
  407. else:
  408. if cur_len > best_len:
  409. best_start, best_len = cur_start, cur_len
  410. cur_len = 0
  411. if cur_len > best_len:
  412. best_start, best_len = cur_start, cur_len
  413. if best_len < int(MIN_ACTIVE_SECONDS * fs):
  414. return 0, n - 1
  415. expand = int(ACTIVE_EXPAND_S * fs)
  416. return max(0, best_start - expand), min(n - 1, best_start + best_len - 1 + expand)
  419. def extract_candidates(evidence, fs, freq_estimate, active_start, active_end):
  420. """Find candidate discharge peaks within the active interval.
  422. Two rounds of peak detection:
  423. 1. Standard: height > 5% of max, min distance = 0.2 × T
  424. 2. Strong: height > 50% of max, min distance = 0.1 × T (captures obvious peaks)
  425. Both sets are merged and deduplicated.
  426. """
  427. segment = evidence[active_start:active_end + 1]
  428. if len(segment) < 3:
  429. return np.array([], dtype=int)
  430. T = 1.0 / freq_estimate if freq_estimate > 0 else 1.0
  431. min_dist = max(20, int(0.2 * T * fs))
  432. min_height = PEAK_HEIGHT_FRAC * np.max(segment)
  433. peaks, _ = find_peaks(segment, height=min_height, distance=min_dist)
  434. strong_height = 0.5 * np.max(segment)
  435. strong_peaks, _ = find_peaks(segment, height=strong_height,
  436. distance=max(10, int(0.1 * T * fs)))
  437. peaks = np.unique(np.concatenate([peaks, strong_peaks]))
  438. return peaks + active_start
  441. def dp_best_sequence(candidates, evidence, fs, freq_estimate):
  442. """Find optimal discharge sequence via forward dynamic programming.
  444. Scoring:
  445. Node score = evidence[candidate]^1.5 (superlinear: favor strong peaks)
  446. Edge score = -α × (deviation_from_periodic)² - β × (skips - 1)
  447. where deviation = (dt - m×T) / (m×T) for skip count m
  448. New-sequence cost = -λ per candidate
  450. The DP finds the path through candidates that maximizes:
  451. Σ(node_scores) + Σ(edge_scores) - λ × path_length
  452. """
  453. if len(candidates) == 0:
  454. return np.array([], dtype=int)
  455. if len(candidates) == 1:
  456. return candidates.copy()
  457. T = 1.0 / freq_estimate
  458. n = len(candidates)
  459. raw_scores = np.array([evidence[c] for c in candidates])
  460. node_scores = raw_scores ** 1.5
  462. best_score = np.full(n, -np.inf)
  463. best_prev = np.full(n, -1, dtype=int)
  464. for i in range(n):
  465. best_score[i] = node_scores[i] - DP_LAMBDA
  467. for j in range(1, n):
  468. for i in range(j):
  469. dt = (candidates[j] - candidates[i]) / fs
  470. if dt <= 0 or dt > 4 * T:
  471. continue
  472. best_edge = -np.inf
  473. for m in range(1, MAX_SKIP + 1):
  474. deviation = (dt - m * T) / (m * T)
  475. interval_score = -DP_ALPHA * deviation ** 2
  476. skip_penalty = -DP_BETA * (m - 1)
  477. edge = interval_score + skip_penalty
  478. if edge > best_edge:
  479. best_edge = edge
  480. total = best_score[i] + best_edge + node_scores[j] - DP_LAMBDA
  481. if total > best_score[j]:
  482. best_score[j] = total
  483. best_prev[j] = i
  485. path = []
  486. idx = int(np.argmax(best_score))
  487. while idx >= 0:
  488. path.append(idx)
  489. idx = best_prev[idx]
  490. path.reverse()
  491. return candidates[np.array(path)]
  494. def em_refine(evidence, discharge_samples, fs, freq_estimate):
  495. """Refine discharge times using template cross-correlation.
  497. Steps:
  498. 1. Extract evidence snippets (±150ms) around each detected discharge
  499. 2. Average to create a discharge template
  500. 3. Normalized cross-correlation of template with full evidence trace
  501. 4. Re-run peak detection + DP on the correlation peaks
  502. """
  503. n = len(evidence)
  504. half_win = int(TEMPLATE_HALF_MS / 1000.0 * fs)
  505. if len(discharge_samples) < 2:
  506. return discharge_samples
  508. snippets = []
  509. for s in discharge_samples:
  510. lo, hi = s - half_win, s + half_win + 1
  511. if lo >= 0 and hi <= n:
  512. snippets.append(evidence[lo:hi])
  513. if len(snippets) < 2:
  514. return discharge_samples
  516. template = np.mean(snippets, axis=0) - np.mean(np.mean(snippets, axis=0))
  517. t_norm = np.sqrt(np.sum(template ** 2))
  518. if t_norm < 1e-10:
  519. return discharge_samples
  521. corr = np.zeros(n)
  522. t_len = len(template)
  523. for i in range(half_win, n - half_win):
  524. seg = evidence[i - half_win:i + half_win + 1]
  525. if len(seg) != t_len:
  526. continue
  527. seg_c = seg - np.mean(seg)
  528. s_norm = np.sqrt(np.sum(seg_c ** 2))
  529. if s_norm < 1e-10:
  530. continue
  531. corr[i] = np.dot(template, seg_c) / (t_norm * s_norm)
  533. T = 1.0 / freq_estimate if freq_estimate > 0 else 1.0
  534. min_dist = max(30, int(0.3 * T * fs))
  535. min_height = 0.15 * np.max(corr) if np.max(corr) > 0 else 0
  536. refined_candidates, _ = find_peaks(corr, height=min_height, distance=min_dist)
  537. if len(refined_candidates) < 2:
  538. return discharge_samples
  539. return dp_best_sequence(refined_candidates, evidence, fs, freq_estimate)
  542. def per_channel_times(evidence_all, global_samples, fs=FS):
  543. """Refine global discharge times per channel within ±50ms.
  545. For each global discharge, finds the local evidence peak in each
  546. channel within a ±50ms window.
  547. """
  548. refine_win = int(CHANNEL_REFINE_MS / 1000.0 * fs)
  549. n_channels, n_samples = evidence_all.shape
  550. channel_times = {}
  551. for ch in range(n_channels):
  552. ch_times = []
  553. for gs in global_samples:
  554. lo = max(0, gs - refine_win)
  555. hi = min(n_samples, gs + refine_win + 1)
  556. window = evidence_all[ch, lo:hi]
  557. if len(window) == 0:
  558. ch_times.append(gs / fs)
  559. else:
  560. ch_times.append((lo + int(np.argmax(window))) / fs)
  561. channel_times[ch] = ch_times
  562. return channel_times
  565. # ═══════════════════════════════════════════════════════════════════════
  566. # Post-hoc confidence filtering
  567. # ═══════════════════════════════════════════════════════════════════════
  569. def posthoc_filter(discharge_samples, evidence, min_evidence_ratio=0.3):
  570. """Drop detected discharges with evidence below a fraction of the median.
  572. Removes low-confidence detections where the evidence peak at the
  573. detected discharge time is less than min_evidence_ratio × median
  574. peak value across all detections.
  575. """
  576. if len(discharge_samples) < 2 or min_evidence_ratio <= 0:
  577. return discharge_samples
  578. peak_values = np.array([evidence[int(s)] for s in discharge_samples])
  579. threshold = min_evidence_ratio * np.median(peak_values)
  580. mask = peak_values >= threshold
  581. return discharge_samples[mask]
  584. # ═══════════════════════════════════════════════════════════════════════
  585. # ACF frequency estimation (handcrafted, no neural network)
  586. # ═══════════════════════════════════════════════════════════════════════
  588. def estimate_frequency_acf(signal_1d, fs=FS):
  589. """Estimate dominant frequency from autocorrelation of pointiness trace.
  591. Steps:
  592. 1. Compute pointiness trace
  593. 2. Autocorrelation of the pointiness trace
  594. 3. Find first peak in ACF after min lag (0.4s = max 2.5 Hz)
  595. 4. Frequency = 1 / lag_at_peak
  596. """
  597. pt = compute_pointiness_trace(signal_1d)
  598. if np.max(pt) < 1e-10:
  599. return np.nan
  601. # Autocorrelation
  602. pt_centered = pt - np.mean(pt)
  603. n = len(pt_centered)
  604. acf = np.correlate(pt_centered, pt_centered, mode='full')[n-1:]
  605. if acf[0] > 0:
  606. acf = acf / acf[0]
  608. # Find first peak after minimum lag
  609. min_lag = int(0.4 * fs) # 0.4s = max 2.5 Hz
  610. max_lag = min(int(3.0 * fs), len(acf) - 1) # 3.0s = min 0.33 Hz
  611. if max_lag <= min_lag:
  612. return np.nan
  614. segment = acf[min_lag:max_lag + 1]
  615. peaks, props = find_peaks(segment, height=0.1)
  616. if len(peaks) == 0:
  617. return np.nan
  619. # Take the first (shortest lag) peak
  620. best_lag = peaks[0] + min_lag
  621. return fs / best_lag
  624. # ═══════════════════════════════════════════════════════════════════════
  625. # Main detector class
  626. # ═══════════════════════════════════════════════════════════════════════
  628. class DischargeDetector:
  629. """Automated periodic discharge detector from raw EEG.
  631. Loads CNN+Attention (frequency estimation) and CET-UNet (evidence)
  632. model ensembles on initialization, then detects discharges from raw
  633. 18-channel bipolar EEG segments.
  635. Args:
  636. cet_model_dir: directory containing cet_unet_fold{0-4}.pt
  637. cnn_model_dir: directory containing cnn_attn_fold{0-4}.pt
  638. device: 'cpu', 'cuda', or 'mps'
  639. n_folds: number of ensemble models to load
  640. """
  642. def __init__(self, cet_model_dir=None, cnn_model_dir=None,
  643. device=None, n_folds=5):
  644. if cet_model_dir is None:
  645. cet_model_dir = PROJECT_DIR / 'data' / 'cet_cache'
  646. if cnn_model_dir is None:
  647. cnn_model_dir = PROJECT_DIR / 'data' / 'pd_channel_cache'
  648. if device is None:
  649. if torch.backends.mps.is_available():
  650. device = torch.device('mps')
  651. elif torch.cuda.is_available():
  652. device = torch.device('cuda')
  653. else:
  654. device = torch.device('cpu')
  655. self.device = device
  657. self.cet_models = self._load_models(
  658. CETUNet, cet_model_dir, 'cet_unet_fold', n_folds)
  659. self.cnn_models = self._load_models(
  660. ChannelPDNetAttention, cnn_model_dir, 'cnn_attn_fold', n_folds)
  662. def _load_models(self, model_class, model_dir, prefix, n_folds):
  663. models = []
  664. model_dir = Path(model_dir)
  665. for fold in range(n_folds):
  666. path = model_dir / f'{prefix}{fold}.pt'
  667. if not path.exists():
  668. raise FileNotFoundError(f"Model not found: {path}")
  669. model = model_class()
  670. model.load_state_dict(torch.load(
  671. str(path), map_location=self.device, weights_only=True))
  672. model.to(self.device)
  673. model.eval()
  674. models.append(model)
  675. return models
  677. @torch.no_grad()
  678. def estimate_frequency(self, segment_18ch):
  679. """Estimate discharge frequency from raw EEG using CNN+Attention.
  681. Runs each channel through the 5-fold CNN+Attention ensemble.
  682. PD probability is used as weight for frequency averaging
  683. (PD-weighted mean of log-frequency predictions).
  685. Returns: float, estimated frequency in Hz, clipped to [0.3, 3.5]
  686. """
  687. n_channels = min(segment_18ch.shape[0], 18)
  688. all_pd_probs = []
  689. all_log_freqs = []
  691. for ch in range(n_channels):
  692. ch_data = segment_18ch[ch].astype(np.float32).copy()
  693. if not np.all(np.isfinite(ch_data)):
  694. all_pd_probs.append(0.0)
  695. all_log_freqs.append(0.0)
  696. continue
  698. mu, std = np.mean(ch_data), np.std(ch_data)
  699. ch_data = (ch_data - mu) / std if std > 1e-8 else ch_data - mu
  700. x = torch.from_numpy(ch_data[np.newaxis, np.newaxis, :]).to(self.device)
  702. pd_probs, log_freqs = [], []
  703. for model in self.cnn_models:
  704. pd_prob, freq_pred, _ = model(x)
  705. pd_probs.append(pd_prob.item())
  706. log_freqs.append(freq_pred.item())
  708. all_pd_probs.append(np.mean(pd_probs))
  709. all_log_freqs.append(np.mean(log_freqs))
  711. pd_weights = np.array(all_pd_probs)
  712. log_freqs = np.array(all_log_freqs)
  714. weight_sum = pd_weights.sum()
  715. if weight_sum > 1e-6:
  716. weighted_log_freq = np.sum(pd_weights * log_freqs) / weight_sum
  717. else:
  718. weighted_log_freq = np.mean(log_freqs)
  720. return float(np.clip(np.exp(weighted_log_freq), 0.3, 3.5))
  722. @torch.no_grad()
  723. def compute_cet_evidence_channel(self, channel_data):
  724. """Run CET-UNet ensemble on a single channel.
  726. Returns: numpy array (n_samples,) — evidence in [0, 1]
  727. """
  728. x = channel_data.astype(np.float32).copy()
  729. mu, std = np.mean(x), np.std(x)
  730. x = (x - mu) / std if std > 1e-8 else x - mu
  731. x_tensor = torch.from_numpy(x[np.newaxis, np.newaxis, :]).to(self.device)
  732. predictions = []
  733. for model in self.cet_models:
  734. pred = model(x_tensor).squeeze().cpu().numpy()
  735. predictions.append(pred)
  736. return np.mean(predictions, axis=0).astype(np.float32)
  738. def estimate_frequency_acf_multichannel(self, segment_18ch, subtype,
  739. laterality=None):
  740. """Estimate frequency from ACF across multiple channels.
  742. Computes ACF frequency per channel and takes the median across
  743. relevant channels (all for GPD, lateralized for LPD).
  745. Returns: float, estimated frequency in Hz (or NaN if no valid estimate)
  746. """
  747. n_channels = min(segment_18ch.shape[0], 18)
  748. b_lp, a_lp = butter(4, LOWPASS_HZ / (FS / 2), btype='low')
  750. acf_freqs = []
  751. for ch in range(n_channels):
  752. try:
  753. sig = filtfilt(b_lp, a_lp, segment_18ch[ch])
  754. except ValueError:
  755. sig = segment_18ch[ch]
  756. freq = estimate_frequency_acf(sig, FS)
  757. if np.isfinite(freq):
  758. acf_freqs.append(freq)
  760. if not acf_freqs:
  761. return np.nan
  762. return float(np.clip(np.median(acf_freqs), 0.3, 3.5))
  764. def detect(self, segment_18ch, subtype='lpd', laterality=None, refine=True):
  765. """Detect periodic discharges in a 10-second EEG segment.
  767. Args:
  768. segment_18ch: (18, 2000) bipolar EEG at 200 Hz
  769. subtype: 'lpd' or 'gpd'
  770. laterality: 'left', 'right', or None
  771. refine: whether to apply EM template refinement
  773. Returns:
  774. dict with:
  775. global_times: list of discharge times in seconds
  776. channel_times: dict ch -> list of times
  777. frequency: float, IPI-derived frequency (Hz)
  778. freq_estimate_input: float, CNN+ACF ensemble frequency (Hz)
  779. n_discharges: int
  780. active_interval: (start_s, end_s)
  781. """
  782. n_channels = min(segment_18ch.shape[0], 18)
  783. n_samples = segment_18ch.shape[1]
  784. fs = FS
  786. # 1. Estimate frequency: CNN+ACF ensemble (0.8 CNN + 0.2 ACF)
  787. cnn_freq = self.estimate_frequency(segment_18ch)
  788. acf_freq = self.estimate_frequency_acf_multichannel(
  789. segment_18ch, subtype, laterality)
  790. if np.isfinite(acf_freq):
  791. freq_estimate = 0.8 * cnn_freq + 0.2 * acf_freq
  792. else:
  793. freq_estimate = cnn_freq
  794. freq_estimate = float(np.clip(freq_estimate, 0.3, 3.5))
  796. # 2. Compute per-channel evidence (both handcrafted and CNN)
  797. hpp_all = np.zeros((n_channels, n_samples))
  798. cet_all = np.zeros((n_channels, n_samples), dtype=np.float32)
  799. for ch in range(n_channels):
  800. hpp_all[ch] = compute_channel_evidence(segment_18ch[ch], fs)
  801. if np.all(np.isfinite(segment_18ch[ch])):
  802. cet_all[ch] = self.compute_cet_evidence_channel(segment_18ch[ch])
  804. # 3. Aggregate across channels
  805. hpp_agg = aggregate_evidence(hpp_all, subtype, laterality)
  806. cet_agg = aggregate_evidence(cet_all, subtype, laterality)
  808. # 4. Combine: product-boosted max(HPP, CET) with CET thresholding
  809. evidence = combine_evidence(hpp_agg, cet_agg)
  811. # 5. DP inference
  812. active_start, active_end = detect_active_interval(evidence, fs)
  813. candidates = extract_candidates(evidence, fs, freq_estimate,
  814. active_start, active_end)
  815. discharge_samples = dp_best_sequence(candidates, evidence, fs,
  816. freq_estimate)
  818. # 6. EM refinement
  819. if refine and len(discharge_samples) >= 3:
  820. discharge_samples = em_refine(evidence, discharge_samples, fs,
  821. freq_estimate)
  823. # 7. Post-hoc confidence filter: drop weak detections
  824. discharge_samples = posthoc_filter(discharge_samples, evidence)
  826. # 8. Per-channel timing
  827. channel_times = per_channel_times(hpp_all, discharge_samples, fs)
  829. # 8. Output
  830. global_times = (discharge_samples / fs).tolist() if len(discharge_samples) > 0 else []
  832. if len(global_times) >= 2:
  833. ipis = np.diff(global_times)
  834. ipi_freq = 1.0 / float(np.median(ipis))
  835. else:
  836. ipi_freq = np.nan
  838. return {
  839. 'global_times': global_times,
  840. 'channel_times': {int(k): v for k, v in channel_times.items()},
  841. 'frequency': ipi_freq,
  842. 'freq_estimate_input': freq_estimate,
  843. 'n_discharges': len(discharge_samples),
  844. 'active_interval': (active_start / fs, active_end / fs),
  845. }
  848. # ═══════════════════════════════════════════════════════════════════════
  849. # CLI entry point
  850. # ═══════════════════════════════════════════════════════════════════════
  852. if __name__ == '__main__':
  853. import scipy.io as sio
  855. print("Discharge Detector — Self-contained implementation")
  856. print("=" * 60)
  858. detector = DischargeDetector()
  859. print(f"Models loaded on {detector.device}")
  861. # Demo: load first available EEG file
  862. eeg_dir = PROJECT_DIR / 'data' / 'eeg'
  863. mat_files = sorted(eeg_dir.glob('*.mat'))
  864. if mat_files:
  865. mat = sio.loadmat(str(mat_files[0]))
  866. data_key = [k for k in mat.keys() if not k.startswith('_')][0]
  867. seg = mat[data_key]
  868. if seg.shape[0] > seg.shape[1]:
  869. seg = seg.T
  870. seg = seg[:18, :2000]
  871. print(f"\nProcessing {mat_files[0].name}: {seg.shape}")
  873. result = detector.detect(seg, subtype='lpd')
  874. print(f" Discharges found: {result['n_discharges']}")
  875. print(f" Times: {[f'{t:.3f}' for t in result['global_times'][:10]]}")
  876. print(f" CNN freq estimate: {result['freq_estimate_input']:.2f} Hz")
  877. print(f" IPI frequency: {result['frequency']:.2f} Hz")
  878. else:
  879. print("\nNo EEG files found in data/eeg/. Provide a .mat file to test.")

discharge_detector.py at commit 58cd045, under CC-BY-NC-4.0 · at the source

Overview

Authors: Jin Jing1, Chenxi Sun1,2, Tianyu Zhang1, Matthew Byrd1, Alexandra–Maria Tăuţan3, Lara Basovic3, Peter N Hadar3, Marta P Fernandes3, Daniel Goldenholz1, Jennifer Kim4, Aaron F Struck5, Sahar F Zafar3, M Brandon Westover1,2
ORCID iDs: Chenxi Sun, Alexandra–Maria Tăuţan, Peter N Hadar, Marta P Fernandes, Daniel Goldenholz, Aaron F Struck, M Brandon Westover
  1. Department of Neurology, Clinical Data Animation Center, Beth Israel Deaconess Medical Center/Harvard Medical School, Boston, MA, United States of America
  2. Stanford University, Palo Alto, CA, United States of America
  3. Department of Neurology, Massachusetts General Hospital/Harvard Medical School, Boston, MA, United States of America
  4. Department of Neurology, Yale University School of Medicine, New Haven, CT, United States of America
  5. Department of Neurology, Washington University in St. Louis, St. Louis, MO, United States of America
Institutions: Beth Israel Deaconess Medical Center (United States); Harvard University (United States); Massachusetts General Hospital (United States); Yale University (United States); Washington University in St. Louis (United States)
Journal: Journal of neural engineering, volume 23, issue 6, article 066002
Dates: received 17 June 2026; accepted 12 October 2026; published online 3 November 2026; in print 1 December 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1088/1741-2552/aeb331 · PMID 42863487 · PMCID PMC13729739 · OpenAlex W7203690944
Open access: hybrid, a free copy (OpenAlex)
Status: code verified (how the code is found and verified)
Categories: EEG (modality), human (organism), methods / tools (subfield)
Methods: Spectral & time-frequency, Connectivity, Machine learning, Statistics, Preprocessing, Evoked potentials, Physiology & signal measures, fMRI & imaging
Keywords: EEG, periodic discharges, rhythmic delta activity, deep learning, signal processing
MeSH: Delta Rhythm*, Electroencephalography*, Periodicity*, Algorithms, Convolutional Neural Networks, Humans, Pattern Recognition, Automated (* major topic)
Topic: EEG and Brain-Computer Interfaces (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: Amazon Web Services; National Institutes of Health (1r01ns131347-01a1, 1r01ns130119-01a1); National Institute on Aging (R01AG073598, R01AG073410); National Heart, Lung, and Blood Institute (R01HL161253); National Institute of Neurological Disorders and Stroke (R01NS126282)
Citations: not cited yet (Europe PMC); 27 references in the paper

Abstract

Objective. Periodic discharges (PDs) and rhythmic delta activity (RDA) are common electroencephalographic patterns in critically ill patients that require detailed characterization—including lateralization, spatial localization, and frequency estimation—according to the American Clinical Neurophysiology Society 2021 standardized terminology. Manual characterization is subjective, time-consuming, and exhibits substantial inter-rater variability. This study presents a comprehensive automated system for characterizing all four major subtypes: lateralized PDs (LPDs), generalized PDs (GPD), lateralized RDA (LRDA), and generalized RDA (GRDA). Approach. Two complementary pipelines were developed: the PD-Profiler and the RDA-Profiler. The PD-Profiler combines a per-channel convolutional neural network (Channelpd-Net) with hemisphere-specific learned evidence traces (HemiCET-UNet), dynamic programming with an approximately-periodic prior, and discharge-locked topographic localization. The RDA-Profiler uses iterative narrowband Hilbert refinement for frequency and lateralization, with phase-locking value analysis for spatial extent. Both pipelines were trained and evaluated on 12 425 EEG segments from 11 729 unique patients. Main results. Four electroencephalographers—three (SZ, TZ, AS) not involved in algorithm development, plus the primary annotator (MW) — independently scored 200 stratified segments per pattern subtype, with analysis restricted to segments accepted by a majority of the four raters. The algorithm exceeded expert–expert (EE) inter-rater agreement on LPD frequency (mean expert–algorithm (EA) intraclass correlation coefficient (ICC) 0.931 vs. EE 0.897; paired segment-bootstrap Δ=+0.034, 95% confidence interval [+0.016,+0.055], p<0.001), and was statistically indistinguishable from EE on every other attribute tested: GPD frequency (EA ICC 0.980 vs. EE 0.973; Δ=+0.007, [−0.000,+0.017], p=0.061), GRDA frequency (Δ=−0.020, p=0.056), LRDA frequency (Δ=−0.055, p=0.435), LPD laterality (EA κ 0.945 vs. EE 0.974; Δ=−0.030, p=0.079) and LRDA laterality (EA κ 0.961 vs. EE 0.993; Δ=−0.032, [−0.087,+0.004], p=0.269). No attribute fell significantly below the EE ceiling. PD discharge timing achieved an F1 of 0.817 at single-sample precision (mean absolute timing error 3.1 ms at 200 Hz; median error 0 samples). Pattern lateralization area under the receiver operating characteristic curve was 0.989 for PD hemisphere (n=1276), 0.909 for LPD vs. GPD (n=7037), and 0.837 for RDA (n=4253). PD spatial localization reached 97.3% of expert inter-rater Jaccard agreement; RDA spatial extent matched EE ICC (0.371 vs. 0.373), at a ceiling set by the inter-rater noise of verbal region-counting labels. Significance. This work presents the first system to reach expert-level inter-rater reliability across lateralization, spatial localization, discharge timing and frequency for both periodic and rhythmic EEG patterns: it exceeded the EE ceiling on LPD frequency and was statistically indistinguishable from that ceiling on every other attribute tested, with none falling significantly below it. Agreement remains lowest, though not significantly below the ceiling, on LRDA laterality and LRDA frequency, which we identify as the priority targets for further work. In post-hoc review of discordant cases, experts judged the algorithm’s frequency estimates more accurate than the original expert labels in 94% of cases. Automated characterization can now both substitute for and improve manual annotation in critical-care EEG.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 40 matches between paragraphs and lines of code.

bdsp-core/grond

License: CC-BY-NC-4.0
State: the link answers, verified on 4 November 2026
Evidence: files inventoried
Commit: 58cd04551393b2b2282baf051ed87cfb2f5ba366, 28 September 2026
Languages: Python (429), Shell (2), Jupyter (2), JavaScript (1)
Size: 837 files, 434 scripts
Software Heritage: not archived
Found in: “Data availability statement”
Holds: README, license file, environment (code/environment.yml), documentation, 2 notebooks
Not found: CITATION.cff, tests, continuous integration
Tools: NumPy (372 files), SciPy (307 files), pandas (150 files), MNE-Python (96 files), Matplotlib (88 files), PyTorch (79 files), scikit-learn (66 files), h5py (23 files), Pillow (9 files), specparam (formerly FOOOF) (6 files), Pingouin (3 files), scikit-image (3 files), seaborn (2 files), LightGBM (1 file), Numba (1 file)
Availability: 1 check, the latest on 4 November 2026: the link answers
  • 4 November 2026: the link answers
436 files
At the source: github.com/bdsp-core/grond

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 434 scripts, each with its path and the digest of its content;
  • 40 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability statement

The data that support the findings of this study are openly available at the following URL/DOI: https://github.com/bdsp-core/grond [27].

Algorithm methods and contests available at: https://doi.org/10.1088/1741-2552/aeb331/data1.

Code. The source code for the GROND system—both pipelines, all algorithm contest entries, the montage-ablation scripts, and the figure- and table-generation scripts used to produce this manuscript—is openly available at: https://github.com/bdsp-core/grond. The repository includes REPRODUCIBILITY.md, which maps every figure and table to the script that regenerates it, and paper_materials/generate_all_figures.py, which regenerates the complete figure set from a single entry point. Trained model weights for ChannelPD-Net, HemiCET-UNet and the BIPD classifier are distributed with the data (below). Data. The EEG segments, expert-rater label files, algorithm predictions, cohort flags and trained model weights are released through the Brain Data Science Platform (https://bdsp.io) as the project ‘GROND: Automated Characterization of Periodic and Rhythmic EEG Patterns’. The data are hosted in the credentialed BDSP Amazon S3 bucket under the prefix s3://bdsp-opendata-credentialed/grond/ (2.0 GB) and comprise: a single self-contained HDF5 bank, grond_data.h5, holding all 14 307 ten-second EEG segments as 19-channel monopolar referential recordings at 200 Hz together with every label, algorithm prediction and cohort flag keyed by segment identifier; the trained model weights under data/; and the four self-contained HTML review tools under independent_expert_tasks/ that the independent expert raters used, so that the exact rating interface behind the IRR results can be inspected and re-run. After credentialing, the data can be retrieved with aws s3 sync s3://bdsp-opendata-credentialed/grond/./grond-data/ Access conditions. The EEG data are de-identified (surrogate patient identifiers and shifted dates) but remain governed by the BDSP Data Use Agreement, so access is credentialed rather than open: investigators register at https://bdsp.io, complete human-subjects research training, and sign the Data Use Agreement. The code carries no such restriction. Questions about data access may be directed to.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 4 November 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 13 authors, 5 keywords, 7 MeSH terms, 5 funders, 26 references.

Cite

This paper

Jing, J., Sun, C., Zhang, T., Byrd, M., Tăuţan, A., Basovic, L., Hadar, P. N., Fernandes, M. P., Goldenholz, D., Kim, J., Struck, A. F., Zafar, S. F., & Westover, M. B. (2026). GROND: automated characterization of periodic and rhythmic EEG patterns. Journal of neural engineering, 23(6), 066002. https://doi.org/10.1088/1741-2552/aeb331

BibTeX

@article{jing2026grond,
author = {Jing, Jin and Sun, Chenxi and Zhang, Tianyu and Byrd, Matthew and Tăuţan, Alexandra–Maria and Basovic, Lara and Hadar, Peter N and Fernandes, Marta P and Goldenholz, Daniel and Kim, Jennifer and Struck, Aaron F and Zafar, Sahar F and Westover, M Brandon},
title = {{GROND: automated characterization of periodic and rhythmic EEG patterns}},
journal = {Journal of neural engineering},
year = {2026},
month = nov,
volume = {23},
number = {6},
pages = {066002},
publisher = {IOP Publishing},
issn = {1741-2560},
doi = {10.1088/1741-2552/aeb331},
url = {https://doi.org/10.1088/1741-2552/aeb331},
pmid = {42863487},
pmcid = {PMC13729739}
}

RIS

TY - JOUR
AU - Jing, Jin
AU - Sun, Chenxi
AU - Zhang, Tianyu
AU - Byrd, Matthew
AU - Tăuţan, Alexandra–Maria
AU - Basovic, Lara
AU - Hadar, Peter N
AU - Fernandes, Marta P
AU - Goldenholz, Daniel
AU - Kim, Jennifer
AU - Struck, Aaron F
AU - Zafar, Sahar F
AU - Westover, M Brandon
TI - GROND: automated characterization of periodic and rhythmic EEG patterns
T2 - Journal of neural engineering
J2 - J Neural Eng
PY - 2026
DA - 2026/11/03
VL - 23
IS - 6
SP - 066002
SN - 1741-2560
PB - IOP Publishing
DO - 10.1088/1741-2552/aeb331
UR - https://doi.org/10.1088/1741-2552/aeb331
LA - en
ER -

CSL-JSON

{
"id": "10.1088/1741-2552/aeb331",
"type": "article-journal",
"title": "GROND: automated characterization of periodic and rhythmic EEG patterns",
"container-title": "Journal of neural engineering",
"author": [
{
"family": "Jing",
"given": "Jin"
},
{
"family": "Sun",
"given": "Chenxi"
},
{
"family": "Zhang",
"given": "Tianyu"
},
{
"family": "Byrd",
"given": "Matthew"
},
{
"family": "Tăuţan",
"given": "Alexandra–Maria"
},
{
"family": "Basovic",
"given": "Lara"
},
{
"family": "Hadar",
"given": "Peter N"
},
{
"family": "Fernandes",
"given": "Marta P"
},
{
"family": "Goldenholz",
"given": "Daniel"
},
{
"family": "Kim",
"given": "Jennifer"
},
{
"family": "Struck",
"given": "Aaron F"
},
{
"family": "Zafar",
"given": "Sahar F"
},
{
"family": "Westover",
"given": "M Brandon"
}
],
"container-title-short": "J Neural Eng",
"volume": "23",
"issue": "6",
"page": "066002",
"DOI": "10.1088/1741-2552/aeb331",
"PMID": "42863487",
"PMCID": "PMC13729739",
"ISSN": "1741-2560",
"publisher": "IOP Publishing",
"URL": "https://doi.org/10.1088/1741-2552/aeb331",
"language": "en",
"issued": {
"date-parts": [
[
2026,
11,
3
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1111/epi.18487 [code]
Harvard Electroencephalography Database: A comprehensive clinical electroencephalographic resource from four Boston hospitals.
Journal: Epilepsia
In common: MNE-Python, PyTorch, scikit-learn, 4 other tools, EEG, 3 references, 4 authors
[2] doi:10.1111/epi.18431 [code]
Automated detection of interictal epileptiform discharges with few electroencephalographic channels.
Journal: Epilepsia
In common: PyTorch, scikit-learn, pandas, 3 other tools, EEG, 1 reference, 2 authors
[3] doi:10.1186/s40708-025-00270-1 [code]
Biotuner: A python toolbox integrating music theory and signal processing for harmonic analysis of physiological and natural time series.
Journal: Brain informatics
In common: specparam (formerly FOOOF), Pingouin, Numba, 9 other tools, EEG, 1 reference
[4] doi:10.3390/biomedicines11030685 [code]
Sample Size Analysis for Machine Learning Clinical Validation Studies.
Journal: Biomedicines
In common: seaborn, scikit-learn, pandas, 3 other tools, 2 authors
[5] doi:10.7554/elife.101071 [code]
Modulation of aperiodic EEG activity provides sensitive index of cognitive state changes during working memory task.
Journal: eLife
In common: specparam (formerly FOOOF), Pingouin, Numba, 9 other tools, EEG
[6] doi:10.7554/elife.110588 [code]
Opening the black box toward a modular approach to spike sorting.
Journal: eLife
In common: LightGBM, Numba, scikit-image, 9 other tools, methods / tools
[7] doi:10.64898/2026.01.23.701239 [code]
Opening the black box: a modular approach to spike sorting
Journal: bioRxiv (preprint)
In common: LightGBM, Numba, scikit-image, 9 other tools, methods / tools
[8] doi:10.1111/jsr.70309 [code]
An Electroencephalographic Study of Sleep Spindle and Infraslow Oscillation in Children With Autism Spectrum Disorder.
Journal: Journal of sleep research
In common: Pingouin, MNE-Python, seaborn, 5 other tools, EEG, author M Brandon Westover
[9] doi:10.1590/0004-282x-anp-2022-s113
Status epilepticus: review on diagnosis, monitoring and treatment.
Journal: Arquivos de neuro-psiquiatria
In common: EEG, 7 references
[10] doi:10.1016/j.isci.2023.107522 [code]
Deep learning for risk-based stratification of cognitively impaired individuals.
Journal: iScience
In common: LightGBM, Pingouin, scikit-image, 8 other tools

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Add the badge to the code's README

OSCR: the badge of this paper

If you are an author or a maintainer of this paper's code, you can add this badge. It says what OSCR knows of the paper (its code, its validated tracing map) and links the code to this page. Paste one of these snippets under the README's title. Keep the two marks around it, so that the badge can be replaced in place and is never added twice. OSCR changes nothing in your repository.

Markdown
<!-- oscr-badge:start -->
[![OSCR](https://openscicode.org/badge/doi/10.1088/1741-2552/aeb331.svg)](https://openscicode.org/badge/doi/10.1088/1741-2552/aeb331)
<!-- oscr-badge:end -->
reStructuredText
.. oscr-badge:start

.. image:: https://openscicode.org/badge/doi/10.1088/1741-2552/aeb331.svg
   :target: https://openscicode.org/badge/doi/10.1088/1741-2552/aeb331
   :alt: OSCR

.. oscr-badge:end
HTML
<!-- oscr-badge:start -->
<a href="https://openscicode.org/badge/doi/10.1088/1741-2552/aeb331"><img src="https://openscicode.org/badge/doi/10.1088/1741-2552/aeb331.svg" alt="OSCR"></a>
<!-- oscr-badge:end -->

On GitHub, open the README, choose Edit (the pencil), paste the snippet under the title, then "Commit changes".

If you choose to, you can have the badge kept up to date for you, with oscr badge --write from the command line or with the GitHub Action in your repository's workflows. Every other form and option (a card, AsciiDoc, LaTeX, shields.io) is on the badge page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: code errors and mismatches between the code and the paper about this paper and its code, from readers and its authors, beside the paper and its tracing map.

Report a mismatch between the code and the paper · Report a code error (signed in)

Reproductions: reports from readers who ran the authors' code and did not get the paper's result, what they ran, with which environment, commit and data.

Report a reproduction failure (signed in)

Versions of its code: the releases of its code tied to a version of this paper (the preprint, the accepted manuscript, the version of record…), each with the tracing map versioned with it.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.