Merge pull request #17 from BabaSanfour/fix/core-updates

Bugfix: ZapLine sampling rate mismatch, returning patterns for plotting & Changelog setup

Author: BabaSanfour

.github/actions/rename_towncrier/rename_towncrier.py

@@ -0,0 +1,59 @@
+#!/usr/bin/env python3
+
+# Authors: The MNE-Python contributors.
+# License: BSD-3-Clause
+# Copyright the MNE-Python contributors.
+# Copied from mne-python:
+# https://github.com/mne-tools/mne-python/blob/main/.github/actions/rename_towncrier/rename_towncrier.py
+
+import json
+import os
+import re
+import subprocess
+import sys
+from pathlib import Path
+
+from github import Github
+from tomllib import loads
+
+event_name = os.getenv("GITHUB_EVENT_NAME", "pull_request")
+if not event_name.startswith("pull_request"):
+    print(f"No-op for {event_name}")
+    sys.exit(0)
+if "GITHUB_EVENT_PATH" in os.environ:
+    with open(os.environ["GITHUB_EVENT_PATH"], encoding="utf-8") as fin:
+        event = json.load(fin)
+    pr_num = event["number"]
+    basereponame = event["pull_request"]["base"]["repo"]["full_name"]
+    real = True
+else:  # local testing
+    pr_num = 12318  # added some towncrier files
+    basereponame = "mne-tools/mne-python"
+    real = False
+
+g = Github(os.environ.get("GITHUB_TOKEN"))
+baserepo = g.get_repo(basereponame)
+
+# Grab config from upstream's default branch
+toml_cfg = loads(Path("pyproject.toml").read_text("utf-8"))
+
+config = toml_cfg["tool"]["towncrier"]
+pr = baserepo.get_pull(pr_num)
+modified_files = [f.filename for f in pr.get_files()]
+
+# Get types from config
+types = [ent["directory"] for ent in toml_cfg["tool"]["towncrier"]["type"]]
+type_pipe = "|".join(types)
+
+# Get files that potentially match the types
+directory = toml_cfg["tool"]["towncrier"]["directory"]
+assert directory.endswith("/"), directory
+
+file_re = re.compile(rf"^{directory}({type_pipe})\.rst$")
+found_stubs = [f for f in modified_files if file_re.match(f)]
+for stub in found_stubs:
+    fro = stub
+    to = file_re.sub(rf"{directory}{pr_num}.\1.rst", fro)
+    print(f"Renaming {fro} to {to}")
+    if real:
+        subprocess.check_call(["mv", fro, to])

CONTRIBUTING.md

@@ -129,6 +129,16 @@ git fetch upstream
 git rebase upstream/main
 ```
 
+### 5. Add Changelog Entry
+
+We use [towncrier](https://towncrier.readthedocs.io/) to manage our changelog. This prevents merge conflicts and ensures standardized release notes.
+
+When you create a Pull Request, please add a changelog entry file in `docs/changes/devel/`. The file name should be the change type (e.g., `feature.rst`, `bugfix.rst`).
+
+For detailed instructions and available types, see [docs/changes/README.md](https://github.com/mne-tools/mne-denoise/blob/main/docs/changes/README.md).
+
+**Author Attribution**: We encourage contributors to include their name in the changelog entry if they wish to be highlighted. In Markdown, you can link to your GitHub profile (e.g., `... (by [@YourUser](...))`).
+
 ## Code Style
 
 We use **Ruff** for linting and formatting, configured to follow PEP 8 with NumPy docstring conventions.

docs/changes/README.md

@@ -0,0 +1,35 @@
+# Changelog Guide
+
+We use `towncrier` to manage our changelog. This ensures that changes are documented as they happen, preventing merge conflicts in the changelog file and ensuring high-quality release notes.
+
+## Adding a Changelog Entry
+
+When you make a change (feature, bugfix, documentation update), you should add a fragment file to the `docs/changes/devel/` directory.
+
+The filename should be the type of change and the extension `.rst`. The PR number will be added automatically.
+
+Format: `<TYPE>.rst`
+
+### Available types:
+
+* `feature`: New feature.
+* `bugfix`: Bug fix.
+* `doc`: Documentation improvement.
+* `removal`: Deprecation or removal of a feature.
+* `misc`: Internal changes, tooling, etc.
+
+## Example
+
+If you fixed a bug in a PR, create a file `docs/changes/devel/bugfix.rst`:
+
+```rst
+Fixed a bug where the ZapLine algorithm would crash on empty data.
+```
+
+## Building the Changelog
+
+To preview the changelog (requires `towncrier`):
+
+```bash
+towncrier build --draft
+```

docs/changes/devel/10.feature.rst

@@ -0,0 +1,5 @@
+**DSS Refactor**:
+- Consolidated DSS implementation from `dev-benchmarch` branch, including:
+  - 20+ denoiser classes in `mne_denoise.dss.denoisers`.
+  - Updated `mne_denoise.utils` for consistent MNE object handling.
+  - Visualization updates.

docs/changes/devel/17.bugfix.rst

@@ -0,0 +1,2 @@
+**ZapLine**:
+- Fixed a bug in `ZapLine` adaptive mode where sampling rate mismatch caused incorrect frequency detection and potential crashes (Issue #16).

docs/changes/template.jinja

@@ -0,0 +1,14 @@
+{% for section, _ in sections.items() %}
+{% if section %}
+### {{ section }}
+
+{% endif %}
+{% for category, val in definitions.items() if category in sections[section] %}
+#### {{ definitions[category]['name'] }}
+
+{% for text, values in sections[section][category].items() %}
+- {{ text }} ({{ values|join(', ') }})
+{% endfor %}
+
+{% endfor %}
+{% endfor %}

docs/conf.py

@@ -31,7 +31,7 @@
 ]
 
 templates_path = ["_templates"]
-exclude_patterns: list[str] = ["_build", "Thumbs.db", ".DS_Store"]
+exclude_patterns: list[str] = ["_build", "Thumbs.db", ".DS_Store", "changes"]
 
 autosummary_generate = True
 napoleon_google_docstring = False

mne_denoise/dss/denoisers/spectral.py

@@ -223,60 +223,56 @@ def _apply_fft(self, data: np.ndarray) -> np.ndarray:
             raise ValueError(f"Data must be 2D or 3D, got {data.ndim}D")
 
     def _get_target_indices(self, nfft: int) -> list:
-        """Get FFT bin indices for target frequencies."""
-        freq_bins = np.fft.fftfreq(nfft, 1 / self.sfreq)
+        """Get FFT bin indices for target frequencies.
+
+        Selects exactly one bin per harmonic (no neighbor padding).
+        Negative-frequency conjugates are included automatically for
+        real-valued IFFT reconstruction.
+        """
         target_indices = []
 
         for f in self._harmonic_freqs:
-            idx = np.argmin(np.abs(freq_bins - f))
-            if idx not in target_indices:
+            # Positive-frequency bin: round(f / sfreq * nfft)
+            idx = int(round(f / self.sfreq * nfft))
+            if 0 <= idx < nfft and idx not in target_indices:
                 target_indices.append(idx)
-            # Negative frequency
-            idx_neg = np.argmin(np.abs(freq_bins + f))
-            if idx_neg not in target_indices:
+
+            # Negative-frequency (conjugate symmetric) bin
+            idx_neg = nfft - idx
+            if 0 <= idx_neg < nfft and idx_neg not in target_indices:
                 target_indices.append(idx_neg)
 
         return target_indices
 
     def _apply_fft_2d(self, data: np.ndarray) -> np.ndarray:
-        """Apply bias to 2D data using FFT."""
+        """Apply bias to 2D data using FFT.
+
+        Process the data in non-overlapping rectangular blocks of length
+        *nfft* (no windowing, no overlap-add). Short trailing blocks are
+        zero-padded to *nfft* and the output is truncated to the true block
+        length.
+        """
         n_channels, n_times = data.shape
 
         # Use data length or nfft, whichever is smaller
         actual_nfft = min(self.nfft, n_times)
         target_indices = self._get_target_indices(actual_nfft)
 
-        # If data is shorter than nfft, process as single block
-        if n_times <= actual_nfft:
-            X = fft(data, n=actual_nfft, axis=1)
-            X_bias = np.zeros_like(X)
-            for idx in target_indices:
-                X_bias[:, idx] = X[:, idx]
-            biased = np.real(ifft(X_bias, axis=1))[:, :n_times]
-            return biased
-
-        # Welch-style block processing
-        step = int(actual_nfft * (1 - self.overlap))
-        step = max(step, 1)
-
         biased = np.zeros_like(data)
-        counts = np.zeros(n_times)
+        pos = 0
 
-        for start in range(0, n_times - actual_nfft + 1, step):
-            end = start + actual_nfft
-            segment = data[:, start:end]
+        while pos < n_times:
+            end = min(pos + actual_nfft, n_times)
+            block_len = end - pos
 
-            X = fft(segment, axis=1)
+            # FFT (zero-pads short blocks automatically)
+            X = fft(data[:, pos:end], n=actual_nfft, axis=1)
             X_bias = np.zeros_like(X)
             for idx in target_indices:
                 X_bias[:, idx] = X[:, idx]
+            y = np.real(ifft(X_bias, axis=1))
 
-            segment_biased = np.real(ifft(X_bias, axis=1))
-            biased[:, start:end] += segment_biased
-            counts[start:end] += 1
-
-        # Normalize by overlap counts
-        counts = np.maximum(counts, 1)
-        biased /= counts
+            biased[:, pos:end] = y[:, :block_len]
+            pos = end
 
         return biased

mne_denoise/dss/denoisers/temporal.py

@@ -177,21 +177,34 @@ def __init__(self, window: int = 10, iterations: int = 1) -> None:
         self.iterations = iterations
 
     def apply(self, data: np.ndarray) -> np.ndarray:
-        """Apply smoothing bias."""
-        from scipy.ndimage import uniform_filter1d
+        """Apply smoothing bias.
 
+        Uses a causal running-mean filter:
+        ``y[t] = mean(x[t-W+1 : t+1])`` for ``t >= W``, with an expanding
+        window for the first ``W`` samples.  Repeated ``iterations`` passes
+        approximate a Gaussian kernel.
+        """
         orig_shape = data.shape
         if data.ndim == 3:
             data_2d = data.reshape(data.shape[0], -1)
         else:
             data_2d = data
 
+        W = int(self.window)
         smoothed = data_2d.copy()
+
         for _ in range(self.iterations):
-            # Use axis=-1 to support 1D (n_times) and 2D (n_ch, n_times)
-            smoothed = uniform_filter1d(
-                smoothed, size=self.window, axis=-1, mode="reflect"
-            )
+            mean_head = np.mean(smoothed[..., : W + 1], axis=-1, keepdims=True)
+            centered = smoothed - mean_head
+
+            # Causal running mean via cumulative sums
+            cs = np.cumsum(centered, axis=-1)
+            out = np.empty_like(centered)
+            # First W samples: expanding window
+            out[..., :W] = cs[..., :W] / np.arange(1, W + 1)
+            # Remaining samples: fixed-width causal window
+            out[..., W:] = (cs[..., W:] - cs[..., :-W]) / W
+            smoothed = out + mean_head
 
         if data.ndim == 3:
             return smoothed.reshape(orig_shape)

mne_denoise/dss/linear.py

@@ -209,7 +209,7 @@ def compute_dss(
     # =========================================================================
     dss_filters = unmixing_matrix.T
 
-    # DSS patterns: for interpretation
+    # DSS patterns: L2-normalized for topographic visualization (Haufe et al. 2014)
     dss_patterns = covariance_baseline @ unmixing_matrix
     pattern_norms = np.sqrt(np.sum(dss_patterns**2, axis=0))
     pattern_norms = np.where(pattern_norms > 1e-15, pattern_norms, 1.0)