Spatially-distributed channel subsets for the sweeps (#91)

benchmarks/benchmark_decompose.py

@@ -318,6 +318,9 @@ def _compare(dirs, figure=None, montage=None, topo_figure=None, n_topo=6, data=N
                 {
                     "label": str(z["label"]),
                     "channels": int(z["channels"]),
+                    "channel_indices": z["channel_indices"]
+                    if "channel_indices" in z
+                    else np.arange(int(z["channels"])),
                     "frames": int(z["frames"]) if "frames" in z else -1,
                     "time": float(z["time"]),
                     "final_ll": float(z["final_ll"]),
@@ -395,16 +398,25 @@ def _compare(dirs, figure=None, montage=None, topo_figure=None, n_topo=6, data=N
     if matrix_figure:
         _plot(labels, M, target, matrix_figure)
     if topo_figure and montage:
-        # de-sphere must recompute the sphere from the SAME frames the group was fit
-        # on (#90 --frames can truncate); pass the matching slice, not the full array.
+        # de-sphere must recompute the sphere from the SAME frames AND channels the
+        # group was fit on (#90 --frames truncates frames; #91 selects distributed
+        # channels), so slice full to both, not the whole array.
+        sel = group[0]["channel_indices"]
         nf = target_cf[1]
-        topo_data = full[:, :nf] if (full is not None and nf > 0) else full
-        _plot_topomaps(group, montage, target, topo_figure, n_topo, topo_data)
+        if full is None:
+            topo_data = None
+        else:
+            topo_data = full[sel][:, :nf] if nf > 0 else full[sel]
+        _plot_topomaps(group, montage, sel, topo_figure, n_topo, topo_data)
 
 
-def _load_info(montage_tsv, n_channels):
+def _load_info(montage_tsv, channel_indices):
     """Build an MNE Info with a scalp montage from a BIDS electrodes.tsv.
 
+    ``channel_indices`` are the 0-based data-channel indices actually used (issue
+    #91 selects a distributed subset, so they need not be contiguous); the
+    electrode name for data channel ``i`` is ``EEG{i+1:03d}``.
+
     NOTE: the tsv coordinates may need a rotation to MNE's head frame (nose +y)
     for the absolute orientation to be correct -- to be checked/fixed later. A
     global rotation does NOT affect the cross-backend equivalence (every backend
@@ -423,7 +435,7 @@ def _load_info(montage_tsv, n_channels):
                 )
             except ValueError:
                 continue
-    ch_names = [f"EEG{i:03d}" for i in range(1, n_channels + 1)]
+    ch_names = [f"EEG{i + 1:03d}" for i in channel_indices]
     # some channels can be unlocalized ('n/a'); topomap only the located subset
     located = np.array([i for i, n in enumerate(ch_names) if n in pos])
     used = [ch_names[i] for i in located]
@@ -435,20 +447,23 @@ def _load_info(montage_tsv, n_channels):
     return info, located
 
 
-def _plot_topomaps(group, montage_tsv, channels, path, n_comps, data=None):
+def _plot_topomaps(group, montage_tsv, channel_indices, path, n_comps, data=None):
     """Grid of IC scalp maps: rows = backends, cols = components. Columns are
     ordered by the reference's back-projected variance (EEGLAB convention, IC1 =
     highest variance) when ``data`` is given, else by cross-backend match. Each
     cell is that backend's Hungarian-matched, sign-aligned map; identical columns
-    down the rows mean every backend recovered the same IC."""
+    down the rows mean every backend recovered the same IC.
+
+    ``channel_indices`` are the 0-based data channels used (issue #91)."""
     import matplotlib
 
     matplotlib.use("Agg")
     import matplotlib.pyplot as plt
     import mne
     from scipy.optimize import linear_sum_assignment
 
-    info, located = _load_info(montage_tsv, channels)
+    info, located = _load_info(montage_tsv, channel_indices)
+    n_ch = len(channel_indices)
     ref = group[0]  # reference backend: the first result loaded (glob order)
     Wref = ref["W"] / (np.linalg.norm(ref["W"], axis=1, keepdims=True) + 1e-12)
     n_ref = Wref.shape[0]
@@ -467,7 +482,7 @@ def _plot_topomaps(group, montage_tsv, channels, path, n_comps, data=None):
 
     # column order: EEGLAB back-projected variance if data available, else match
     if data is not None:
-        order = _variance_order(ref["W"], data, channels)[:n_comps]
+        order = _variance_order(ref["W"], data, n_ch)[:n_comps]
     else:
         order = np.argsort(-score)[: min(n_comps, n_ref)]
 
@@ -481,7 +496,7 @@ def _plot_topomaps(group, montage_tsv, channels, path, n_comps, data=None):
         col, corr = matches[bi]
         # de-sphered mixing = true sensor-space IC scalp maps (EEGLAB); the saved
         # A is whitened-space. Fall back to A only if no data to recompute the sphere.
-        ades = _desphere(r["W"], data, channels)[0] if data is not None else r["A"]
+        ades = _desphere(r["W"], data, n_ch)[0] if data is not None else r["A"]
         for ci, i in enumerate(order):
             j = col[i]
             sign = 1.0 if corr[i, j] >= 0 else -1.0
@@ -502,7 +517,7 @@ def _plot_topomaps(group, montage_tsv, channels, path, n_comps, data=None):
             transform=axes[bi][0].transAxes,
         )
     fig.suptitle(
-        f"AMICA IC scalp maps @ {channels} channels -- matched across backends "
+        f"AMICA IC scalp maps @ {n_ch} channels -- matched across backends "
         "(rows), same IC (columns)",
         fontsize=10,
     )
@@ -615,10 +630,24 @@ def main() -> int:
         else [full.shape[1]]
     )
     for nc in channels:
+        # #91: with a montage, pick spatially-distributed (whole-head) channels
+        # rather than the first nc electrodes in file order (a spatial cluster);
+        # the selected data-channel indices are saved so the topomaps use them.
+        if args.montage and nc < full.shape[0]:
+            from channel_selection import (
+                positions_for_channels,
+                select_distributed_channels,
+            )
+
+            sel = select_distributed_channels(
+                positions_for_channels(args.montage, full.shape[0]), nc
+            )
+        else:
+            sel = np.arange(nc)
         for nf in frame_list:
-            data = np.ascontiguousarray(full[:nc, :nf])
-            k = nf / nc**2
-            print(f"\n{nc}ch, {nf} frames (k={k:.1f})")
+            data = np.ascontiguousarray(full[sel][:, :nf])
+            k = nf / len(sel) ** 2
+            print(f"\n{len(sel)}ch, {nf} frames (k={k:.1f})")
             for b in backends:
                 label = f"{b}@{tag}"
                 try:
@@ -629,7 +658,8 @@ def main() -> int:
                         fname,
                         label=label,
                         backend=b,
-                        channels=nc,
+                        channels=len(sel),
+                        channel_indices=sel,
                         frames=nf,
                         iters=args.iters,
                         time=res["time"],

benchmarks/benchmark_dimsweep.py

@@ -508,6 +508,11 @@ def main() -> int:
     ap = argparse.ArgumentParser(description=__doc__)
     ap.add_argument("--data", help="path to the (70, n_samples) real-EEG .npy")
     ap.add_argument("--channels", default="16,32,48,70")
+    ap.add_argument(
+        "--montage",
+        help="BIDS electrodes.tsv; when given, reduced channel counts use "
+        "spatially-distributed (whole-head) subsets instead of the first N (#91)",
+    )
     ap.add_argument("--samples", type=int, default=30000)
     ap.add_argument("--iters", type=int, default=30)
     ap.add_argument("--repeats", type=int, default=3)
@@ -608,7 +613,24 @@ def _thread_points(b):
 
     rows = []
     for nc in channels:
-        data = np.ascontiguousarray(full[:nc, : args.samples])
+        # #91: with a montage, time whole-head distributed subsets rather than the
+        # first nc electrodes (a spatial cluster). Channel count is unchanged, so the
+        # timing is unaffected; this only makes the reduced montages physical.
+        if args.montage and nc < full.shape[0]:
+            from channel_selection import (
+                positions_for_channels,
+                select_distributed_channels,
+            )
+
+            sel = select_distributed_channels(
+                positions_for_channels(args.montage, full.shape[0]), nc
+            )
+        else:
+            sel = np.arange(nc)
+        # a montage with fewer localized electrodes than requested yields fewer
+        # channels than nc; record/print the actual count so the sweep is honest.
+        n_sel = len(sel)
+        data = np.ascontiguousarray(full[sel][:, : args.samples])
         for b in backends:
             for t in _thread_points(b):
                 tag = f"@{t}t" if t is not None else ""
@@ -626,26 +648,26 @@ def _thread_points(b):
                     rows.append(
                         {
                             "config": config,
-                            "channels": nc,
+                            "channels": n_sel,
                             "backend": b,
                             "threads": t,
                             "ms_per_iter": ms,
                             "final_ll": ll,
                         }
                     )
                     print(
-                        f"  [{config}] {nc:3d}ch {b:18s}{tag:5s} "
+                        f"  [{config}] {n_sel:3d}ch {b:18s}{tag:5s} "
                         f"{ms:9.2f} ms/it   LL={ll:+.5f}"
                     )
                 except Exception as exc:  # noqa: BLE001 - report and continue
                     print(
-                        f"  [{config}] {nc:3d}ch {b:18s}{tag:5s} "
+                        f"  [{config}] {n_sel:3d}ch {b:18s}{tag:5s} "
                         f"FAILED: {type(exc).__name__}: {str(exc)[:60]}"
                     )
                     rows.append(
                         {
                             "config": config,
-                            "channels": nc,
+                            "channels": n_sel,
                             "backend": b,
                             "threads": t,
                             "error": str(exc)[:120],

benchmarks/channel_selection.py

@@ -0,0 +1,116 @@
+"""Spatially-distributed channel subset selection (issue #91).
+
+The dimension sweeps reduce channel count by slicing ``full[:nc]`` -- the *first*
+nc electrodes in file order. For a montage like ds002718 (EEG001..EEG070) the
+first 16/32/48 electrodes are a spatial cluster (one scalp region), not a
+whole-head cap, so the reduced-channel decompositions and their IC scalp maps are
+not physically meaningful reduced montages.
+
+This module picks nc electrodes **evenly distributed across the scalp** via greedy
+farthest-point (k-center) sampling over the real 3D electrode positions from a
+BIDS ``electrodes.tsv``, so every subset is a proper whole-head montage. Selection
+is deterministic (seeded from the electrode nearest the montage centroid).
+"""
+
+from __future__ import annotations
+
+import csv
+
+import numpy as np
+
+
+def read_electrode_positions(tsv_path):
+    """Parse a BIDS ``electrodes.tsv`` into ordered names and 3D positions.
+
+    Rows whose ``x``/``y``/``z`` are non-numeric (e.g. ``n/a`` for unlocalized or
+    non-EEG channels) are skipped.
+
+    Parameters
+    ----------
+    tsv_path : str or path-like
+        Path to a BIDS ``electrodes.tsv`` (tab-separated, columns ``name x y z``).
+
+    Returns
+    -------
+    names : list of str
+        Electrode names, in file order, for the localized channels.
+    positions : ndarray of shape (n_localized, 3)
+        Corresponding 3D coordinates.
+    """
+    names, pos = [], []
+    with open(tsv_path) as f:
+        for row in csv.DictReader(f, delimiter="\t"):
+            try:
+                xyz = [float(row["x"]), float(row["y"]), float(row["z"])]
+            except (ValueError, KeyError, TypeError):
+                continue
+            names.append(row["name"])
+            pos.append(xyz)
+    return names, np.asarray(pos, dtype=float)
+
+
+def positions_for_channels(tsv_path, n_channels, name_fmt="EEG{:03d}"):
+    """Positions aligned to data-channel index for the first ``n_channels``.
+
+    Data channel ``i`` (0-based) is assumed to be electrode ``name_fmt.format(i+1)``
+    (``EEG001``, ``EEG002``, ...), matching the benchmark data layout. Channels
+    without a localized position are filled with NaN.
+
+    Returns
+    -------
+    ndarray of shape (n_channels, 3)
+        Row ``i`` is the position of data channel ``i`` (NaN if unlocalized).
+    """
+    names, pos = read_electrode_positions(tsv_path)
+    lookup = dict(zip(names, pos))
+    out = np.full((n_channels, 3), np.nan, dtype=float)
+    for i in range(n_channels):
+        p = lookup.get(name_fmt.format(i + 1))
+        if p is not None:
+            out[i] = p
+    return out
+
+
+def select_distributed_channels(positions, n):
+    """Greedy farthest-point (k-center) selection of ``n`` spread-out channels.
+
+    Starting from the channel nearest the centroid, repeatedly add the channel
+    that is farthest (in the max-min sense) from the already-selected set. This
+    yields a spatially distributed, whole-head subset rather than a cluster.
+
+    Parameters
+    ----------
+    positions : ndarray of shape (n_channels, 3)
+        Positions aligned to data-channel index; rows may be NaN (unlocalized),
+        which are excluded from selection.
+    n : int
+        Number of channels to select.
+
+    Returns
+    -------
+    ndarray of int
+        Sorted data-channel indices of the selected channels. If ``n`` is at least
+        the number of localized channels, all localized indices are returned.
+    """
+    positions = np.asarray(positions, dtype=float)
+    localized = np.where(np.isfinite(positions).all(axis=1))[0]
+    pts = positions[localized]
+    if n <= 0:
+        return localized[:0]
+    if n >= len(localized):
+        return localized
+    centroid = pts.mean(axis=0)
+    first = int(np.argmin(((pts - centroid) ** 2).sum(axis=1)))
+    selected = [first]
+    # running min distance from every point to the selected set
+    dist = np.sqrt(((pts - pts[first]) ** 2).sum(axis=1))
+    # exclude already-selected points from argmax so coincident coordinates
+    # (dist ties at 0) can never be picked twice -- keeps the subset unique.
+    avail = np.ones(len(pts), dtype=bool)
+    avail[first] = False
+    while len(selected) < n:
+        nxt = int(np.argmax(np.where(avail, dist, -np.inf)))
+        selected.append(nxt)
+        avail[nxt] = False
+        dist = np.minimum(dist, np.sqrt(((pts - pts[nxt]) ** 2).sum(axis=1)))
+    return np.sort(localized[np.asarray(selected, dtype=int)])