diff --git a/src/phosphor/channel_grid.py b/src/phosphor/channel_grid.py new file mode 100644 index 0000000..ae0fce8 --- /dev/null +++ b/src/phosphor/channel_grid.py @@ -0,0 +1,354 @@ +"""GPU-accelerated 2D grid of one value per channel. + +Draws one square per channel at its real ``(x, y)`` position (the ``(x, y)`` is +the square's **lower-left corner**), sized per channel, and colored by value. +The squares are a single quad mesh, so they live in real world-units: the +camera's starting bounds are the grid's own bounds with aspect preserved, so +squares render square, and a channel reported as smaller than the array pitch +draws as a smaller square with a gap rather than filling its cell. + +Colour comes from one of two places. By default a value is normalized over +``vmin``/``vmax`` and run through ``cmap``, which suits a quantity that varies +continuously. A caller that classifies its values instead -- into states with +names and fixed colours, as an impedance reading or a channel-quality flag +would be -- passes the colours it computed to :meth:`ChannelGridWidget.set_colors` +and the words to :meth:`ChannelGridWidget.set_annotations`. What counts as good +is domain vocabulary, so it stays with the domain; this module renders whatever +it is handed. + +Falls back to a unit-spaced sequential layout when every channel shares one +position, which is what an unmapped source looks like. +""" + +from __future__ import annotations + +import logging +from dataclasses import dataclass + +import cmap as cmap_lib +import fastplotlib as fpl +import numpy as np +from PySide6 import QtCore, QtWidgets + +from .grid_layout import build_quad_mesh_arrays, resolve_cell_geometry + +logger = logging.getLogger(__name__) + +__all__ = [ + "ChannelGridConfig", + "ChannelGridWidget", +] + +RGBA = tuple[float, float, float, float] + + +@dataclass +class ChannelGridConfig: + positions: np.ndarray + """``(n_ch, 2)`` array of electrode ``(x, y)`` lower-left corners, in the + same units as ``sizes`` (micrometers for a CMP/device-mapped source).""" + + sizes: np.ndarray | None = None + """``(n_ch,)`` per-channel square side length. ``None`` (or a non-positive + entry) falls back to the inferred electrode pitch so squares tile their + cells, matching the old full-cell heatmap look.""" + + channel_labels: list[str] | None = None + cmap: str = "viridis_r" + vmin: float = 0.0 + vmax: float = 1000.0 + nan_color: RGBA = (0.15, 0.15, 0.18, 1.0) + """Colour for a channel with no value. Applies to the colormap path only -- + explicit colours from ``set_colors`` are used verbatim, since a caller that + computes its own colours is the one that knows what missing should look + like.""" + + show_values: bool = False + value_format: str = "{:.0f}" + text_face_color: str = "white" + text_outline_color: str = "black" + text_outline_thickness: float = 0.4 + text_font_size: float = 13.0 + value_unit: str = "" + invert_y: bool = True + """Flip the layout so low-``y`` positions render at the TOP of the screen. + Matches the "row 0 is the first row I read" expectation most people have of + a grid, and the row-0-at-bottom convention of electrode maps. ``False`` + keeps low-``y`` at the bottom, the native y-up world.""" + + +class ChannelGridWidget(QtWidgets.QWidget): + """Embeddable grid: one square per channel at its electrode position, + color-mapped by value, with optional value text and a channel tooltip.""" + + def __init__(self, config: ChannelGridConfig, parent: QtWidgets.QWidget | None = None) -> None: + super().__init__(parent) + self._config = config + n_ch = config.positions.shape[0] + + # Per-channel display rectangles: (x0, y0, side) with (x0, y0) the + # lower-left corner in display space (already y-flipped if requested), + # plus square centers for value text and tooltip hit-testing. + self._rects, self._centers = self._resolve_geometry(config, n_ch) + + self._values_per_ch = np.full(n_ch, np.nan, dtype=np.float32) + self._dirty = True + self._needs_autoscale = True + self._show_values = bool(config.show_values) + self._cmap = cmap_lib.Colormap(config.cmap) + self._nan_rgba = self._rgba_f(config.nan_color) + # Caller-supplied colours and per-channel words. None means "use the + # colormap" and "say nothing beyond the number". + self._explicit_rgba: np.ndarray | None = None + self._annotations: list[str] | None = None + + qt_layout = QtWidgets.QVBoxLayout(self) + qt_layout.setContentsMargins(0, 0, 0, 0) + + self._figure = fpl.Figure() + self._subplot = self._figure[0, 0] + self._fpl_widget = self._figure.show() + qt_layout.addWidget(self._fpl_widget) + + positions, indices = self._build_mesh_arrays() + init_colors = np.tile(self._nan_rgba, (positions.shape[0], 1)) + self._mesh = self._subplot.add_mesh(positions, indices, mode="basic", colors=init_colors) + + self._text_graphics: list = [] + if self._show_values: + self._create_text_graphics() + + self._tooltip = QtWidgets.QLabel(self._fpl_widget) + self._tooltip.setStyleSheet( + "background: rgba(25,25,30,220); color: #e8e8e8;" + " padding: 12px 24px; font-size: 36pt;" + " font-family: 'Menlo','Consolas','DejaVu Sans Mono',monospace;" + " border: 1px solid rgba(120,120,120,160);" + ) + self._tooltip.setAttribute(QtCore.Qt.WidgetAttribute.WA_TransparentForMouseEvents) + self._tooltip.hide() + + self._subplot.controller = None + self._subplot.axes.visible = False + self._subplot.title.visible = False + # Disable fastplotlib's built-in hover tooltip — it shows raw geometry + # data, meaningless for our value-coded squares; we provide our own + # channel-aware tooltip via _on_pointer_move. + try: + self._subplot.tooltip.enabled = False + except AttributeError: + pass + self._figure.add_animations(lambda: self._animation_callback()) + + renderer = self._subplot.renderer + renderer.add_event_handler(self._on_pointer_move, "pointer_move") + renderer.add_event_handler(self._on_pointer_leave, "pointer_leave") + + # ---- Public API ---------------------------------------------------- + + def push_data(self, values: np.ndarray) -> None: + v = np.asarray(values, dtype=np.float32).reshape(-1) + if v.size != self._values_per_ch.size: + logger.warning("push_data size %d != n_channels %d; ignoring.", v.size, self._values_per_ch.size) + return + self._values_per_ch = v + self._dirty = True + + def set_show_values(self, show: bool) -> None: + if show == self._show_values: + return + self._show_values = show + if show and not self._text_graphics: + self._create_text_graphics() + for _, tg in self._text_graphics: + tg.visible = show + self._dirty = True + + def set_colors(self, rgba: np.ndarray | None) -> None: + """Colour each channel explicitly, overriding the colormap. + + For values that mean something categorical rather than continuous -- a + pass/fail, a quality state -- where the mapping is the caller's to + decide. ``None`` hands colouring back to the colormap. + """ + if rgba is None: + self._explicit_rgba = None + self._dirty = True + return + arr = np.asarray(rgba, dtype=np.float32).reshape(-1, 4) + if arr.shape[0] != self._values_per_ch.size: + logger.warning("set_colors got %d rows for %d channels; ignoring.", arr.shape[0], self._values_per_ch.size) + return + self._explicit_rgba = arr + self._dirty = True + + def set_annotations(self, annotations: list[str] | None) -> None: + """Per-channel words for the tooltip, after the value. + + The counterpart to :meth:`set_colors`: a caller that colours a channel + by a state it named should be able to show that name on hover, without + this module having to know the vocabulary. + """ + if annotations is not None and len(annotations) != self._values_per_ch.size: + logger.warning( + "set_annotations got %d entries for %d channels; ignoring.", + len(annotations), + self._values_per_ch.size, + ) + return + self._annotations = list(annotations) if annotations is not None else None + self._dirty = True + + def set_color_range(self, vmin: float, vmax: float) -> None: + if vmin > vmax: + vmin, vmax = vmax, vmin + if vmin == self._config.vmin and vmax == self._config.vmax: + return + self._config.vmin = float(vmin) + self._config.vmax = float(vmax) + self._dirty = True + + def clear_values(self) -> None: + """Reset all squares to NaN. Use before changing acquisition source so + stale values don't linger until the new source overwrites them.""" + self._values_per_ch[:] = np.nan + self._dirty = True + + def close_figure(self) -> None: + """Tear down the fastplotlib figure deterministically. + + Without this, rendercanvas's deferred Qt callbacks can fire after + the C++ widget has been deleted, raising a shiboken RuntimeError + at app exit (or whenever this widget is replaced). + """ + figure = getattr(self, "_figure", None) + if figure is None: + return + try: + figure.close() + except Exception: + logger.exception("fastplotlib figure close raised; ignoring") + self._figure = None + + @property + def show_values(self) -> bool: + return self._show_values + + @property + def color_range(self) -> tuple[float, float]: + return self._config.vmin, self._config.vmax + + # ---- Geometry ------------------------------------------------------ + + @staticmethod + def _resolve_geometry(config: ChannelGridConfig, n_ch: int) -> tuple[np.ndarray, np.ndarray]: + """Per-channel display rects ``(x0, y0, side)`` and centers ``(cx, cy)``. + + Thin wrapper over the shared :func:`resolve_cell_geometry`. + """ + return resolve_cell_geometry(config.positions, config.sizes, n_ch, config.invert_y) + + def _build_mesh_arrays(self) -> tuple[np.ndarray, np.ndarray]: + """Vertex positions ``(n*4, 3)`` and triangle indices ``(n*2, 3)`` for + the per-channel quads.""" + return build_quad_mesh_arrays(self._rects) + + # ---- Internals ----------------------------------------------------- + + @staticmethod + def _rgba_f(rgba: RGBA) -> np.ndarray: + return np.asarray(rgba, dtype=np.float32) + + def _create_text_graphics(self) -> None: + for ch_idx in range(self._values_per_ch.size): + cx, cy = self._centers[ch_idx] + tg = self._subplot.add_text( + "", + font_size=self._config.text_font_size, + face_color=self._config.text_face_color, + outline_color=self._config.text_outline_color, + outline_thickness=self._config.text_outline_thickness, + # z above the mesh (z=0) so the text isn't occluded. + offset=(float(cx), float(cy), 1.0), + anchor="middle-center", + ) + tg.visible = self._show_values + self._text_graphics.append((ch_idx, tg)) + + def _animation_callback(self) -> None: + if not self._dirty: + return + self._dirty = False + + rgba = self._channel_rgba() # (n_ch, 4) float + # One color per quad → repeat across the quad's 4 vertices. + self._mesh.colors[:] = np.repeat(rgba, 4, axis=0) + + if self._show_values: + fmt = self._config.value_format + for ch_idx, tg in self._text_graphics: + v = self._values_per_ch[ch_idx] + tg.text = "" if np.isnan(v) else fmt.format(float(v)) + + if self._needs_autoscale: + self._needs_autoscale = False + # Fit the camera to the electrode grid's own bounds, equal scale on + # both axes so the squares render square (rectilinear aspect). + self._subplot.auto_scale(maintain_aspect=True, zoom=0.95) + + def _channel_rgba(self) -> np.ndarray: + if self._explicit_rgba is not None: + return self._explicit_rgba + return self._render_continuous() + + def _render_continuous(self) -> np.ndarray: + v = self._values_per_ch + cmap_span = max(self._config.vmax - self._config.vmin, 1e-12) + norm = (v - self._config.vmin) / cmap_span + norm_safe = np.nan_to_num(np.clip(norm, 0.0, 1.0), nan=0.0) + rgba = self._cmap(norm_safe).astype(np.float32) + rgba[np.isnan(v)] = self._nan_rgba + return rgba + + def _channel_at(self, wx: float, wy: float) -> int: + """Index of the square containing world point ``(wx, wy)``, or -1.""" + x0 = self._rects[:, 0] + y0 = self._rects[:, 1] + side = self._rects[:, 2] + inside = (wx >= x0) & (wx <= x0 + side) & (wy >= y0) & (wy <= y0 + side) + hits = np.flatnonzero(inside) + return int(hits[0]) if hits.size else -1 + + def _on_pointer_move(self, event) -> None: + world = self._subplot.map_screen_to_world(event) + if world is None: + self._tooltip.hide() + return + ch_idx = self._channel_at(float(world[0]), float(world[1])) + if ch_idx < 0: + self._tooltip.hide() + return + + self._tooltip.setText(self._tooltip_text(ch_idx)) + self._tooltip.adjustSize() + x = int(getattr(event, "x", 0)) + 14 + y = int(getattr(event, "y", 0)) + 14 + max_x = self._fpl_widget.width() - self._tooltip.width() - 4 + max_y = self._fpl_widget.height() - self._tooltip.height() - 4 + self._tooltip.move(min(x, max_x), min(y, max_y)) + self._tooltip.show() + + def _tooltip_text(self, ch_idx: int) -> str: + """What hovering a channel says: its name, its value, and whatever the + caller called it.""" + labels = self._config.channel_labels + label = labels[ch_idx] if labels and ch_idx < len(labels) else f"ch{ch_idx}" + note = self._annotations[ch_idx] if self._annotations else "" + value = float(self._values_per_ch[ch_idx]) + if np.isnan(value): + return f"{label}\n{note}" if note else f"{label}\nno value" + unit = f" {self._config.value_unit}" if self._config.value_unit else "" + return f"{label}\n{value:.1f}{unit} ({note})" if note else f"{label}\n{value:.1f}{unit}" + + def _on_pointer_leave(self, _event) -> None: + self._tooltip.hide() diff --git a/src/phosphor/decimate.py b/src/phosphor/decimate.py new file mode 100644 index 0000000..5cacde4 --- /dev/null +++ b/src/phosphor/decimate.py @@ -0,0 +1,112 @@ +"""Display-only min/max (envelope) decimation for waveform plotting. + +A pure, stateless companion to the plotting widgets. When a waveform has far +more samples than the screen has pixels to draw it on, uploading every sample +to the GPU is wasteful — a Utah-array evoked grid at 128 ch × 11 history × +15 000 samples is ~253 MB of vertex data, most of it invisible. + +The fix is to draw an **envelope**: split the sample axis into ``n_out // 2`` +equal buckets and keep each bucket's *min* and *max*. Unlike plain stride +decimation (``ys[..., ::k]``), which lands between samples and clips the P/N +peaks that *are* the evoked signal, min/max preserves peak amplitude at ~2 +points per bucket. + +This is the stateless form of what :class:`~phosphor.sweep_buffer.SweepBuffer` +does on a ring (``_recompute_columns`` + ``_build_multiline_array``). The sweep +reduces samples as they stream past; a grid of retained waveforms has the whole +array in hand and reduces it on demand, so it needs the same arithmetic without +the circular-buffer and thread-lock machinery. + +Decimate raw amplitudes *before* any affine cell-mapping: min/max is +order-preserving under a monotone map, so the drawn envelope is identical, and +keeping it upstream avoids re-clipping surprises. +""" + +from __future__ import annotations + +import warnings + +import numpy as np + +__all__ = ["minmax_decimate", "minmax_decimate_x", "plan_minmax_decimation"] + + +class MinMaxDecimationPlan: + """Precomputed bucket geometry for repeated min/max decimation. + + Built once from ``(n_samples, max_points)``; reused every frame. ``active`` + is False when the waveform already fits the point budget (no decimation) — + callers should pass arrays through unchanged in that case. + """ + + __slots__ = ("active", "bucket_size", "n_buckets", "n_out", "n_samples", "trim") + + def __init__(self, n_samples: int, max_points: int) -> None: + self.n_samples = int(n_samples) + n_buckets = max(1, int(max_points) // 2) + bucket_size = self.n_samples // n_buckets if n_buckets else 0 + # Only decimate when there is something to gain: each bucket must + # collapse >1 raw sample, else the envelope is the original signal. + if self.n_samples <= int(max_points) or bucket_size <= 1: + self.active = False + self.n_buckets = 0 + self.bucket_size = 0 + self.trim = 0 + self.n_out = self.n_samples + return + self.active = True + self.n_buckets = n_buckets + self.bucket_size = bucket_size + # Drop the ragged tail so every bucket is full and reshape is exact. + self.trim = n_buckets * bucket_size + self.n_out = 2 * n_buckets + + +def plan_minmax_decimation(n_samples: int, max_points: int) -> MinMaxDecimationPlan: + """Build a reusable :class:`MinMaxDecimationPlan` for ``n_samples`` waveforms + drawn with at most ``max_points`` vertices each.""" + return MinMaxDecimationPlan(n_samples, max_points) + + +def minmax_decimate(values: np.ndarray, plan: MinMaxDecimationPlan) -> np.ndarray: + """Min/max envelope decimation along the last (sample) axis. + + *values* shape ``(..., n_samples)`` → ``(..., 2 * n_buckets)``: each bucket + contributes its ``min`` then ``max``, in that fixed order. All-NaN buckets + (e.g. unfilled history slots) yield NaN, which callers treat as "empty". + + Returns *values* unchanged when ``plan.active`` is False. + """ + if not plan.active: + return values + if values.shape[-1] != plan.n_samples: + raise ValueError(f"minmax_decimate: last axis {values.shape[-1]} != plan.n_samples {plan.n_samples}") + lead = values.shape[:-1] + vb = values[..., : plan.trim].reshape(*lead, plan.n_buckets, plan.bucket_size) + with warnings.catch_warnings(): + # All-NaN buckets legitimately warn "All-NaN slice"; NaN is the intended + # result, so silence it rather than mask (matches phosphor's handling). + warnings.simplefilter("ignore", RuntimeWarning) + lo = np.nanmin(vb, axis=-1) + hi = np.nanmax(vb, axis=-1) + dec = np.stack([lo, hi], axis=-1).reshape(*lead, plan.n_out) + return dec.astype(values.dtype, copy=False) + + +def minmax_decimate_x(x: np.ndarray, plan: MinMaxDecimationPlan) -> np.ndarray: + """Decimate a monotone x coordinate to match :func:`minmax_decimate`. + + Each bucket is drawn as a (min, max) pair sharing one x — the bucket's + center — so the polyline stays monotonic in time. *x* shape + ``(..., n_samples)`` → ``(..., 2 * n_buckets)`` with each center repeated + for its pair. Returns *x* unchanged when ``plan.active`` is False. + """ + if not plan.active: + return x + if x.shape[-1] != plan.n_samples: + raise ValueError(f"minmax_decimate_x: last axis {x.shape[-1]} != plan.n_samples {plan.n_samples}") + lead = x.shape[:-1] + xb = x[..., : plan.trim].reshape(*lead, plan.n_buckets, plan.bucket_size) + centers = xb.mean(axis=-1) + out = np.repeat(centers, 2, axis=-1) + return out.astype(x.dtype, copy=False) diff --git a/src/phosphor/grid_layout.py b/src/phosphor/grid_layout.py new file mode 100644 index 0000000..e7c8c87 --- /dev/null +++ b/src/phosphor/grid_layout.py @@ -0,0 +1,172 @@ +"""Geometry for laying one cell per channel out at its real position. + +A per-channel value heatmap and a per-channel trace grid draw completely +different things into their cells, but agree entirely on where the cells go: +one square per channel at its ``(x, y)``, sized per channel, with a sensible +layout when the positions are missing or degenerate. That agreement is what +lives here. + +* :func:`tiled_grid_positions` — a square-ish fallback layout for sources that + carry no positions at all. +* :func:`tile_by_group` — fan out groups of channels that share a coordinate + range, so one group does not render on top of another. +* :func:`infer_pitch` — the spacing between adjacent channels, which is what a + cell defaults to when no size is given. +* :func:`resolve_cell_geometry` — ``(positions, sizes)`` into per-channel + rectangles ``(x0, y0, side)`` and centers, honouring ``invert_y`` and the + degenerate "every channel at one point" case. +* :func:`build_quad_mesh_arrays` — tessellate those rects into a quad mesh, for + a renderer that fills cells rather than drawing into them. + +Positions are plain arrays. Deriving them from a particular data model's +channel metadata is that model's business, not this module's. +""" + +from __future__ import annotations + +import logging + +import numpy as np + +logger = logging.getLogger(__name__) + +__all__ = [ + "build_quad_mesh_arrays", + "infer_pitch", + "resolve_cell_geometry", + "tile_by_group", + "tiled_grid_positions", +] + +# Gap between adjacent tiled blocks, in channel pitches. +_GROUP_GUTTER_PITCHES = 2.0 + +# Quad corner offsets (unit square, lower-left origin) and the two triangles +# that tessellate it. Multiplied by each channel's side length and added to its +# lower-left corner to build a mesh. +_CORNERS = np.array([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]], dtype=np.float64) +_QUAD_TRIS = np.array([[0, 1, 2], [0, 2, 3]], dtype=np.uint32) + + +def tiled_grid_positions(n_ch: int) -> np.ndarray: + """Fallback layout for sources without electrode positions (e.g. NWB + playback of plain ``TimeSeries`` streams). + + Packs channels into a roughly-square ``ceil(sqrt(n_ch)) × ceil(n_ch/side)`` + grid so they're all visible as individual tiles rather than stacking on top + of each other at the origin. + """ + if n_ch <= 0: + return np.zeros((0, 2), dtype=np.float32) + side = int(np.ceil(np.sqrt(n_ch))) + idx = np.arange(n_ch) + positions = np.column_stack([idx % side, idx // side]).astype(np.float32) + return positions + + +def tile_by_group(positions: np.ndarray, group_ids: np.ndarray) -> np.ndarray: + """Shift each group's block along x so groups render side by side. + + Devices that are physically separate often report the same coordinate + range -- two arrays each numbering their electrodes from the same origin -- + so their channels land on identical ``(x, y)`` and one block hides another. + Each group keeps its internal layout and is placed to the right of the + previous one, with a :data:`_GROUP_GUTTER_PITCHES`-pitch gutter between + blocks. + + A single distinct group is returned unchanged, so the layout is only + rewritten when there is genuinely more than one block to separate. + """ + positions = np.asarray(positions, dtype=np.float32) + hs = np.asarray(group_ids).reshape(-1) + if hs.shape[0] != positions.shape[0]: + return positions + groups = np.unique(hs) + if groups.size <= 1: + return positions + + xs = positions[:, 0].astype(np.float64) + ys = positions[:, 1].astype(np.float64) + pitch = infer_pitch(xs, ys) + gutter = _GROUP_GUTTER_PITCHES * pitch + + tiled = positions.copy() + cursor = 0.0 + for g in groups: + mask = hs == g + gx = xs[mask] + tiled[mask, 0] = (gx - gx.min() + cursor).astype(np.float32) + # Block width includes the last column's own cell footprint (+pitch). + cursor += (gx.max() - gx.min() + pitch) + gutter + logger.info("Tiled %d channel-group blocks side by side to avoid overlap.", groups.size) + return tiled + + +def infer_pitch(xs: np.ndarray, ys: np.ndarray) -> float: + """Smallest positive spacing between distinct x or y coordinates. + + The channel pitch -- for an electrode array, its inter-electrode spacing. + Defaults to ``1`` when the geometry is degenerate, giving unit-square + tiling. + """ + steps: list[float] = [] + for vals in (np.unique(xs), np.unique(ys)): + if vals.size > 1: + steps.append(float(np.min(np.diff(vals)))) + positive = [s for s in steps if s > 0] + return min(positive) if positive else 1.0 + + +def resolve_cell_geometry( + positions: np.ndarray, + sizes: np.ndarray | None, + n_ch: int, + invert_y: bool, +) -> tuple[np.ndarray, np.ndarray]: + """Per-channel display rects ``(x0, y0, side)`` and centers ``(cx, cy)``. + + ``(x0, y0)`` is each cell's lower-left corner. When ``invert_y`` is set the + layout is flipped about the x-axis so low-``y`` electrodes render at the top + while each cell stays axis-aligned with its corner anchoring intact. + + Falls back to a unit-spaced sequential layout when every channel shares one + position, which is what an unmapped source looks like: no geometry loaded, + so every channel defaulted to the same point. + """ + xs = positions[:, 0].astype(np.float64) + ys = positions[:, 1].astype(np.float64) + + if n_ch and np.allclose(xs, xs[0]) and np.allclose(ys, ys[0]): + # Degenerate (all channels at one point) — lay out a unit-spaced + # square so the cells don't stack on top of each other. + logger.info("All channel positions identical; using sequential unit layout for %d channels.", n_ch) + side = int(np.ceil(np.sqrt(n_ch))) + xs = (np.arange(n_ch) % side).astype(np.float64) + ys = (np.arange(n_ch) // side).astype(np.float64) + cell_sizes = np.ones(n_ch, dtype=np.float64) + else: + pitch = infer_pitch(xs, ys) + if sizes is None: + cell_sizes = np.full(n_ch, pitch, dtype=np.float64) + else: + cell_sizes = np.asarray(sizes, dtype=np.float64).reshape(-1).copy() + # Zero / missing electrode size → fall back to the pitch. + cell_sizes[~(cell_sizes > 0)] = pitch + + y0 = -(ys + cell_sizes) if invert_y else ys + rects = np.column_stack([xs, y0, cell_sizes]) + centers = np.column_stack([xs + cell_sizes / 2.0, y0 + cell_sizes / 2.0]) + return rects, centers + + +def build_quad_mesh_arrays(rects: np.ndarray) -> tuple[np.ndarray, np.ndarray]: + """Vertex positions ``(n*4, 3)`` and triangle indices ``(n*2, 3)`` for the + per-channel quads described by ``rects`` (``(n, 3)`` of ``(x0, y0, side)``).""" + n = rects.shape[0] + # corner_xy[i, k] = lower-left[i] + corner[k] * side[i] + corner_xy = rects[:, None, :2] + _CORNERS[None, :, :] * rects[:, None, 2:3] + positions = np.zeros((n * 4, 3), dtype=np.float32) + positions[:, :2] = corner_xy.reshape(-1, 2) + base = (np.arange(n) * 4)[:, None, None] + indices = (base + _QUAD_TRIS[None, :, :]).reshape(-1, 3).astype(np.uint32) + return positions, indices diff --git a/src/phosphor/trace_grid.py b/src/phosphor/trace_grid.py new file mode 100644 index 0000000..654ef0c --- /dev/null +++ b/src/phosphor/trace_grid.py @@ -0,0 +1,713 @@ +"""GPU-accelerated 2D grid of per-channel waveform plots. + +A reusable companion to :mod:`channel_grid`. Where the channel grid colors one +square per channel by a scalar value, this widget draws a *mini line-plot* per +channel — a rolling history of fixed-length waveforms plus their running mean — +at each channel's real ``(x, y)`` electrode position. + +It is deliberately app-agnostic so two viewers can share it: + +* **evoked potentials** — each marker-locked epoch window is one waveform per + channel (``push_epoch``). +* **spike waveforms** (future) — each detected spike snippet is a waveform; the + synchronized ``push_epoch`` model can be extended to per-channel pushes later. + +Layout, cell geometry, and the camera-anchoring background mesh are shared with +the channel grid via :mod:`grid_layout`. +""" + +from __future__ import annotations + +import logging +from dataclasses import dataclass + +import fastplotlib as fpl +import numpy as np +from PySide6 import QtCore, QtGui, QtWidgets + +from .decimate import minmax_decimate, minmax_decimate_x, plan_minmax_decimation +from .grid_layout import build_quad_mesh_arrays, resolve_cell_geometry +from .trace_grid_buffer import TraceGridBuffer + +logger = logging.getLogger(__name__) + +__all__ = ["TraceGridConfig", "TraceGridWidget"] + +RGBA = tuple[float, float, float, float] + +# World-space viewport rect shared for pan/zoom sync: (xmin, xmax, ymin, ymax). +Viewport = tuple[float, float, float, float] + + +def _viewport_differs(a: Viewport, b: Viewport, *, rtol: float = 1e-4) -> bool: + """True if two world rects differ by more than a small relative tolerance. + + Guards the per-frame camera poll against float jitter that would otherwise + emit a viewport_changed every frame even when the user isn't interacting. + """ + span = max(abs(b[1] - b[0]), abs(b[3] - b[2]), 1e-9) + return any(abs(x - y) > rtol * span for x, y in zip(a, b)) + + +@dataclass +class TraceGridConfig: + positions: np.ndarray + """``(n_ch, 2)`` array of electrode ``(x, y)`` lower-left corners.""" + + n_samples: int + """Number of samples per waveform (epoch window length).""" + + sizes: np.ndarray | None = None + """``(n_ch,)`` per-channel cell side length; ``None`` → inferred pitch.""" + + channel_labels: list[str] | None = None + history: int = 10 + """How many recent waveforms to retain and overlay per channel.""" + + show_individual: bool = True + show_mean: bool = True + + show_error: bool = False + """Draw a +/- one standard deviation band around the mean. + + Needs ``track_statistics``; meaningless without a mean to bound.""" + + track_statistics: bool = True + """Accumulate the running mean and standard deviation over every waveform. + + Off for a stack with no meaningful average -- action potentials from + possibly-different units -- where it costs two float64 arrays and a pass + over every arrival to compute a number nobody draws.""" + + age_fade: bool = True + """Fade older waveforms. + + Costs a rewrite of the whole colour buffer on every arrival, because every + retained waveform's age changes when one arrives. Worth it at evoked rates + and not at spike rates, where a flat colour keeps a push to one slot's + worth of vertices.""" + + autoscale: bool = True + """Auto-fit the shared amplitude range to the retained waveforms.""" + + y_min: float = -100.0 + y_max: float = 100.0 + """Fixed amplitude range used when ``autoscale`` is False.""" + + cell_pad_frac: float = 0.08 + """Fraction of each cell reserved as inner padding around the waveform.""" + + individual_color: RGBA = (0.45, 0.62, 0.95, 1.0) + mean_color: RGBA = (1.0, 0.82, 0.25, 1.0) + error_color: RGBA = (1.0, 0.82, 0.25, 0.45) + cell_color: RGBA = (0.12, 0.12, 0.15, 1.0) + individual_thickness: float = 1.0 + mean_thickness: float = 2.5 + error_thickness: float = 1.0 + + invert_y: bool = False + """Flip so low-``y`` electrodes render at the top (Blackrock ``.cmp`` + row-0-at-bottom convention). Matches :class:`ChannelGridWidget`'s default + of ``False`` for impedance; callers map as they prefer.""" + + x_unit: str = "s" + x_extent: tuple[float, float] | None = None + """Optional ``(t0, t1)`` of the waveform window, used only for the tooltip.""" + + display_max_points: int | None = None + """Max vertices drawn per waveform. ``None`` → derive from the monitor + (~2 samples per horizontal pixel of a single cell, which is all a min/max + envelope can show). Waveforms are min/max-decimated to this budget for + rendering only; the stored history, running mean, and autoscale stay at full + resolution. See :mod:`phosphor.decimate`.""" + + +class TraceGridWidget(QtWidgets.QWidget): + """Embeddable grid: one mini waveform-plot per channel at its electrode + position, overlaying a rolling history of waveforms and their mean.""" + + # Emitted when the user pans/zooms the camera, carrying the new visible + # world rect (xmin, xmax, ymin, ymax). Lets a companion view (e.g. the SNR + # heatmap) follow the same electrode region. Not emitted for programmatic + # :meth:`set_viewport` changes, so linked views don't echo each other. + viewport_changed = QtCore.Signal(object) + + def __init__(self, config: TraceGridConfig, parent: QtWidgets.QWidget | None = None) -> None: + super().__init__(parent) + self._config = config + self._n_ch = config.positions.shape[0] + self._n_samples = int(config.n_samples) + self._show_individual = bool(config.show_individual) + self._show_mean = bool(config.show_mean) + self._show_error = bool(config.show_error) and bool(config.track_statistics) + self._autoscale = bool(config.autoscale) + self._y_min = float(config.y_min) + self._y_max = float(config.y_max) + + # Per-channel display rects (x0, y0, side) and centers, shared layout. + self._rects, self._centers = resolve_cell_geometry(config.positions, config.sizes, self._n_ch, config.invert_y) + # Precompute the per-channel x coordinate of each sample inside its cell. + self._x_line = self._compute_x_line() # (n_ch, n_samples) + + # Display-only min/max envelope decimation: shrink the per-frame vertex + # buffer (and its GPU mirror) without touching the stored/analyzed data. + budget = self._resolve_display_max_points() + self._dec_plan = plan_minmax_decimation(self._n_samples, budget) + # Precompute the decimated x line once; it never changes. + self._x_line_dec = minmax_decimate_x(self._x_line, self._dec_plan) + # Always log the decision (even when inactive) so it's clear whether + # decimation ran — a silent no-op is indistinguishable from "not wired". + if self._dec_plan.active: + logger.info( + "TraceGrid: min/max-decimating %d samples → %d display points per waveform (%.1fx, budget=%d).", + self._n_samples, + self._dec_plan.n_out, + self._n_samples / self._dec_plan.n_out, + budget, + ) + else: + logger.info( + "TraceGrid: no decimation — %d samples already within display budget=%d.", + self._n_samples, + budget, + ) + + # Retention and running statistics live in the buffer; this widget only + # decides how to draw them. + self._buffer = TraceGridBuffer( + self._n_ch, + self._n_samples, + config.history, + track_statistics=config.track_statistics, + ) + self._graphics_version = -1 # buffer version the current graphics were built for + self._dirty = True # flag "the data changed, the screen hasn't caught up yet" + self._needs_autoscale_camera = True + + # Pan/zoom-sync state. ``_last_viewport`` is the camera rect at the last + # poll; ``_suppress_emit`` absorbs the frames right after a programmatic + # set_viewport (the aspect-preserving camera settles over a frame or two) + # so a linked view's change isn't echoed back as a user gesture. + self._last_viewport: Viewport | None = None + self._suppress_emit = False + + qt_layout = QtWidgets.QVBoxLayout(self) + qt_layout.setContentsMargins(0, 0, 0, 0) + + self._figure = fpl.Figure() + self._subplot = self._figure[0, 0] + self._fpl_widget = self._figure.show() + qt_layout.addWidget(self._fpl_widget) + + # Faint background quad per cell — frames each channel and anchors the + # camera (the waveform lines may be all-NaN at first, giving auto_scale + # nothing to fit). Reuses the channel-grid quad tessellation. + bg_positions, bg_indices = build_quad_mesh_arrays(self._rects) + bg_colors = np.tile(np.asarray(config.cell_color, dtype=np.float32), (bg_positions.shape[0], 1)) + self._bg_mesh = self._subplot.add_mesh(bg_positions, bg_indices, mode="basic", colors=bg_colors) + + self._indiv_ml = None + self._mean_ml = None + self._error_ml = None + + # NOTE: assigning None does *not* disable interaction — fastplotlib's + # setter substitutes a default PanZoomController, so the grid keeps + # mouse/trackpad pan + zoom (with maintain_aspect, cells stay square). + # We observe that camera below to broadcast viewport_changed. + self._subplot.controller = None + self._subplot.axes.visible = False + self._subplot.title.visible = False + try: + self._subplot.tooltip.enabled = False + except AttributeError: + pass + + self._tooltip = QtWidgets.QLabel(self._fpl_widget) + self._tooltip.setStyleSheet( + "background: rgba(25,25,30,220); color: #e8e8e8;" + " padding: 8px 14px; font-size: 16pt;" + " font-family: 'Menlo','Consolas','DejaVu Sans Mono',monospace;" + " border: 1px solid rgba(120,120,120,160);" + ) + self._tooltip.setAttribute(QtCore.Qt.WidgetAttribute.WA_TransparentForMouseEvents) + self._tooltip.hide() + + self._figure.add_animations(lambda: self._animation_callback()) + renderer = self._subplot.renderer + renderer.add_event_handler(self._on_pointer_move, "pointer_move") + renderer.add_event_handler(self._on_pointer_leave, "pointer_leave") + renderer.add_event_handler(self._on_double_click, "double_click") + + # ---- Public API ---------------------------------------------------- + + def push_epoch(self, data: np.ndarray) -> None: + """Append one waveform per channel. *data* shape ``(n_samples, n_ch)``. + + Rolls the per-channel history so the new waveform becomes the most + recent; the oldest is dropped once ``history`` is reached. + """ + arr = np.asarray(data, dtype=np.float32) + if arr.ndim != 2 or arr.shape[0] != self._n_samples or arr.shape[1] != self._n_ch: + logger.warning( + "push_epoch shape %s != (n_samples=%d, n_ch=%d); ignoring.", + arr.shape, + self._n_samples, + self._n_ch, + ) + return + self._buffer.push(arr.T) # (n_ch, n_samples) + self._dirty = True + + def set_history(self, history: int) -> None: + """Change how many waveforms are overlaid. + + Retained waveforms survive where they still fit, and the running mean + is untouched: how many are drawn has nothing to do with how many are + averaged. + """ + self._buffer.set_history(history) + self._dirty = True + + def set_show_error(self, show: bool) -> None: + self._show_error = bool(show) and self._buffer.track_statistics + self._dirty = True + + def set_show_individual(self, show: bool) -> None: + self._show_individual = bool(show) + self._dirty = True + + def set_show_mean(self, show: bool) -> None: + self._show_mean = bool(show) + self._dirty = True + + def set_autoscale(self, enabled: bool) -> None: + self._autoscale = bool(enabled) + self._dirty = True + + def set_y_range(self, y_min: float, y_max: float) -> None: + if y_min > y_max: + y_min, y_max = y_max, y_min + self._autoscale = False + self._y_min = float(y_min) + self._y_max = float(y_max) + self._dirty = True + + def clear(self) -> None: + """Drop all retained waveforms (e.g. before swapping acquisition source). + + Repaints immediately (do what _animation_callback does): + fastplotlib renders on-demand, so just flagging + ``_dirty`` would leave the old traces on screen until the next epoch adds + graphics and triggers a draw. We remove the line graphics now and request + a frame so the cleared grid is visible the instant the event is switched. + """ + self._buffer.clear() + self._dirty = False + self._rebuild_lines() # nothing retained → deletes the lines, recreates none + self._request_draw() + + def _request_draw(self) -> None: + """Ask the canvas to present a fresh frame now (on-demand renderer).""" + figure = getattr(self, "_figure", None) + canvas = getattr(figure, "canvas", None) if figure is not None else None + request = getattr(canvas, "request_draw", None) + if callable(request): + request() + + def close_figure(self) -> None: + """Tear down the fastplotlib figure deterministically so rendercanvas + releases its Qt callbacks before the C++ widget is destroyed.""" + figure = getattr(self, "_figure", None) + if figure is None: + return + try: + figure.close() + except Exception: + logger.exception("fastplotlib figure close raised; ignoring") + self._figure = None + + @property + def history(self) -> int: + return self._buffer.history + + @property + def show_individual(self) -> bool: + return self._show_individual + + @property + def show_mean(self) -> bool: + return self._show_mean + + @property + def autoscale(self) -> bool: + return self._autoscale + + @property + def y_range(self) -> tuple[float, float]: + return self._y_min, self._y_max + + # ---- Geometry ------------------------------------------------------ + + def _resolve_display_max_points(self) -> int: + """Vertex budget per waveform for min/max-envelope rendering. + + Honours an explicit ``config.display_max_points``; otherwise derives it + from the monitor — ~2 samples per horizontal pixel of one cell is all an + envelope can resolve. Assumes the grid may span the full screen width + (over-provisions slightly, ignoring side docks — acceptable, and a + future resize-aware variant can recompute from the live canvas). Falls + back to a fixed budget when no screen can be queried. + + **Multi-monitor:** we size for the widest to be safe. + """ + configured = self._config.display_max_points + if configured is not None: + return max(2, int(configured)) + + fallback = 1500 + n_columns = max(1, int(np.unique(self._rects[:, 0]).size)) + try: + screens = QtGui.QGuiApplication.screens() + phys_w = max( + (float(s.geometry().width()) * float(s.devicePixelRatio()) for s in screens), + default=0.0, + ) + except Exception: + logger.debug("TraceGrid: could not query screens for decimation target.", exc_info=True) + return fallback + if phys_w <= 0: + return fallback + per_cell_px = phys_w / n_columns + return max(2, int(2.0 * per_cell_px)) + + def _compute_x_line(self) -> np.ndarray: + x0 = self._rects[:, 0] + side = self._rects[:, 2] + pad = self._config.cell_pad_frac + inner = side * (1.0 - 2.0 * pad) + if self._n_samples > 1: + frac = np.arange(self._n_samples, dtype=np.float64) / (self._n_samples - 1) + else: + frac = np.zeros(self._n_samples, dtype=np.float64) + # (n_ch, n_samples) + return (x0 + pad * side)[:, None] + frac[None, :] * inner[:, None] + + def _map_y(self, values: np.ndarray) -> np.ndarray: + """Map amplitude *values* (..., n_ch, n_samples) into per-cell y.""" + y0 = self._rects[:, 1] + side = self._rects[:, 2] + pad = self._config.cell_pad_frac + inner = side * (1.0 - 2.0 * pad) + span = max(self._y_max - self._y_min, 1e-12) + frac = (values - self._y_min) / span + # Allow a little overflow past the cell, but not into neighbours. + frac = np.clip(frac, -0.1, 1.1) + return (y0 + pad * side)[..., :, None] + frac * inner[..., :, None] + + # ---- Rendering ----------------------------------------------------- + + def _animation_callback(self) -> None: + if self._dirty: + self._dirty = False + if self._autoscale: + self._recompute_autoscale() + self._refresh_lines() + + if self._needs_autoscale_camera: + self._needs_autoscale_camera = False + self._subplot.auto_scale(maintain_aspect=True, zoom=0.95) + + self._poll_viewport() + + # ---- Pan/zoom viewport sync ---------------------------------------- + + def _read_viewport(self) -> Viewport: + """Current visible world rect from the camera state.""" + state = self._subplot.camera.get_state() + zoom = state["zoom"] or 1.0 + w = state["width"] / zoom + h = state["height"] / zoom + cx, cy = float(state["position"][0]), float(state["position"][1]) + return (cx - w / 2.0, cx + w / 2.0, cy - h / 2.0, cy + h / 2.0) + + def _poll_viewport(self) -> None: + """Emit viewport_changed when the camera moved from user interaction. + + Runs every animation frame. A change originating from set_viewport is + absorbed (``_suppress_emit``) so linked views don't ping-pong; once the + rect stops changing we clear the flag and resume emitting user gestures. + """ + rect = self._read_viewport() + if self._last_viewport is not None and not _viewport_differs(rect, self._last_viewport): + self._suppress_emit = False + return + self._last_viewport = rect + if self._suppress_emit: + return + self.viewport_changed.emit(rect) + + def current_viewport(self) -> Viewport: + """The visible world rect (xmin, xmax, ymin, ymax); for initial sync.""" + return self._read_viewport() + + def set_viewport(self, rect: Viewport) -> None: + """Frame *rect* (world coords), matching a linked view's pan/zoom. + + maintain_aspect is preserved, so the region is centered and fit to this + widget's shape (it may reveal a little extra along the wider axis rather + than distorting cells). Does not re-emit viewport_changed. + """ + xmin, xmax, ymin, ymax = rect + cx = (xmin + xmax) / 2.0 + cy = (ymin + ymax) / 2.0 + cam = self._subplot.camera + state = cam.get_state() + state["position"] = (cx, cy, float(state["position"][2])) + state["width"] = max(xmax - xmin, 1e-9) + state["height"] = max(ymax - ymin, 1e-9) + state["zoom"] = 1.0 + cam.set_state(state) + self._suppress_emit = True + self._needs_autoscale_camera = False # user/linked view now owns the camera + self._request_draw() + + def reset_view(self) -> None: + """Refit the camera to all electrodes (double-click). The resulting + viewport is broadcast via the normal poll, so a linked view follows.""" + self._needs_autoscale_camera = True + self._suppress_emit = False # a reset is a user action; always broadcast + self._request_draw() + + def _on_double_click(self, _event) -> None: + self.reset_view() + + def _recompute_autoscale(self) -> None: + valid = self._buffer.recent() + if valid.size == 0 or not np.any(np.isfinite(valid)): + return + lo = float(np.nanmin(valid)) + hi = float(np.nanmax(valid)) + if not np.isfinite(lo) or not np.isfinite(hi): + return + if hi <= lo: + hi = lo + 1.0 + margin = 0.05 * (hi - lo) + self._y_min = lo - margin + self._y_max = hi + margin + + # ---- Drawing ------------------------------------------------------- + # + # Two paths. A waveform arriving changes numbers, not shapes, so its slot is + # written into the graphics that already exist. Only a change of shape -- + # history, or a toggle that adds or removes lines -- recreates them. Doing + # the recreate on every arrival is what made this unusable at spike rates: + # it tears down and re-uploads every retained waveform to draw one new one. + + def _slot_rows(self, slot: int) -> slice: + """Rows of the individual-waveform graphic belonging to one ring slot. + + Lines are laid out slot-major, so a slot's channels are contiguous and + one arrival writes one block. + """ + return slice(slot * self._n_ch, (slot + 1) * self._n_ch) + + def _slot_alpha(self) -> np.ndarray: + """Per-slot opacity, brightest for the newest, zero for never-written. + + Keyed on age rather than position because nothing is reordered when a + waveform arrives -- see :attr:`TraceGridBuffer.ages`. + """ + buf = self._buffer + ages = buf.ages + if not self._config.age_fade: + alpha = np.full(buf.history, 0.9, dtype=np.float32) + else: + ramp = np.linspace(0.9, 0.15, buf.history, dtype=np.float32) + alpha = ramp[ages] + return np.where(ages < buf.n_retained, alpha, 0.0).astype(np.float32) + + def _individual_positions(self) -> tuple[np.ndarray, int]: + """Vertex positions for every slot, ``(history * n_ch, m, 3)``.""" + buf = self._buffer + # Decimate raw amplitudes (peak-preserving) before the cell mapping; + # x is the matching precomputed envelope. Full-res data is untouched. + traces = minmax_decimate(buf.traces, self._dec_plan) + m = traces.shape[-1] + pos = np.empty((buf.history, self._n_ch, m, 3), dtype=np.float32) + pos[..., 0] = np.broadcast_to(self._x_line_dec, (buf.history, self._n_ch, m)) + pos[..., 1] = self._map_y(traces) + pos[..., 2] = 0.0 + return pos.reshape(buf.history * self._n_ch, m, 3), m + + def _individual_colors(self) -> np.ndarray: + """One RGBA per line, ``(history * n_ch, 4)``, faded by slot age. + + What ``add_multi_line`` takes. Writing into an existing graphic needs + the per-vertex form instead -- see :meth:`_expand_colors`. + """ + buf = self._buffer + rgba = np.tile(np.asarray(self._config.individual_color, dtype=np.float32), (buf.history, self._n_ch, 1)) + rgba[..., 3] = self._slot_alpha()[:, None] + return rgba.reshape(-1, 4) + + @staticmethod + def _expand_colors(graphic, per_line: np.ndarray) -> np.ndarray: + """Per-line colours as the flat per-vertex buffer a graphic holds. + + fastplotlib takes one colour per line when a graphic is built and then + stores one per vertex, so an in-place recolour has to repeat each line's + colour across its own run. The stride comes from the buffer rather than + the point count because it also covers whatever separator vertices the + renderer inserted between lines. + """ + total = graphic.colors.value.shape[0] + stride = total // max(per_line.shape[0], 1) + return np.repeat(per_line, stride, axis=0)[:total] + + def _summary_positions(self) -> tuple[np.ndarray | None, np.ndarray | None]: + """Mean and +/- one standard deviation, decimated for display. + + Both are computed at full resolution and reduced only for drawing, so + what is averaged is never what happened to fit on screen. + """ + stats = self._buffer.statistics() + if stats is None: + return None, None + mean, std = stats + mean_pos = self._curve_positions(minmax_decimate(mean, self._dec_plan)) + if not self._show_error or std is None or not np.any(np.isfinite(std)): + return mean_pos, None + band = np.concatenate([mean - std, mean + std], axis=0) + return mean_pos, self._curve_positions(minmax_decimate(band, self._dec_plan)) + + def _curve_positions(self, curve: np.ndarray) -> np.ndarray: + """One line per row of *curve*, laid into the cells.""" + n_lines, m = curve.shape[0], curve.shape[-1] + pos = np.empty((n_lines, m, 3), dtype=np.float32) + # The band is two curves per channel, so x tiles rather than broadcasts. + reps = n_lines // self._n_ch + pos[..., 0] = np.tile(self._x_line_dec, (reps, 1)) + pos[..., 1] = self._map_y(curve) + pos[..., 2] = 0.0 + return pos + + def _graphics_shape(self) -> tuple: + """Which graphics should exist and how big. A change means recreate. + + The flags say *should this be drawn at all*, not just whether the user + asked for it, because the widget is built before its first waveform + arrives. A key that described only sizes and toggles would settle on + "nothing to draw" against the empty buffer and never change when data + turned up -- every arrival would take the in-place path, find no + graphics to write into, and draw nothing. + """ + buf = self._buffer + has_data = buf.n_retained > 0 + return ( + buf.history, + self._n_ch, + self._dec_plan.n_out, + self._show_individual and has_data, + self._show_mean and has_data and buf.track_statistics, + # A spread needs a second waveform, so the band appears later than + # the mean does. + self._show_error and buf.track_statistics and buf.n_seen > 1, + ) + + def _refresh_lines(self) -> None: + """Draw the current buffer, recreating graphics only if the shape moved.""" + if self._graphics_version != self._graphics_shape(): + self._rebuild_lines() + return + self._update_lines() + + def _update_lines(self) -> None: + """Write new numbers into the graphics that already exist.""" + buf = self._buffer + if self._indiv_ml is not None and buf.n_retained: + pos, _ = self._individual_positions() + self._indiv_ml.data[:] = pos + if self._config.age_fade: + self._indiv_ml.colors[:] = self._expand_colors(self._indiv_ml, self._individual_colors()) + mean_pos, band_pos = self._summary_positions() + if self._mean_ml is not None and mean_pos is not None: + self._mean_ml.data[:] = mean_pos + if self._error_ml is not None and band_pos is not None: + self._error_ml.data[:] = band_pos + + def _rebuild_lines(self) -> None: + """Recreate the waveform, mean, and error graphics from scratch. + + Only for a change of shape: the number of lines, their length, or which + of them exist at all. + """ + for name in ("_indiv_ml", "_mean_ml", "_error_ml"): + graphic = getattr(self, name, None) + if graphic is not None: + self._subplot.delete_graphic(graphic) + setattr(self, name, None) + + buf = self._buffer + self._graphics_version = self._graphics_shape() + if buf.n_retained == 0: + return + + if self._show_individual: + pos, _ = self._individual_positions() + self._indiv_ml = self._subplot.add_multi_line( + pos, + colors=self._individual_colors(), + thickness=self._config.individual_thickness, + ) + self._indiv_ml.visible = self._show_individual + + mean_pos, band_pos = self._summary_positions() + # Drawn before the mean so the band sits behind the line it bounds. + if band_pos is not None: + self._error_ml = self._subplot.add_multi_line( + band_pos, + colors=np.tile(np.asarray(self._config.error_color, dtype=np.float32), (band_pos.shape[0], 1)), + thickness=self._config.error_thickness, + ) + self._error_ml.visible = self._show_error + if self._show_mean and mean_pos is not None: + self._mean_ml = self._subplot.add_multi_line( + mean_pos, + colors=np.tile(np.asarray(self._config.mean_color, dtype=np.float32), (self._n_ch, 1)), + thickness=self._config.mean_thickness, + ) + self._mean_ml.visible = self._show_mean + + # ---- Tooltip / hit-testing ----------------------------------------- + + def _channel_at(self, wx: float, wy: float) -> int: + x0 = self._rects[:, 0] + y0 = self._rects[:, 1] + side = self._rects[:, 2] + inside = (wx >= x0) & (wx <= x0 + side) & (wy >= y0) & (wy <= y0 + side) + hits = np.flatnonzero(inside) + return int(hits[0]) if hits.size else -1 + + def _on_pointer_move(self, event) -> None: + world = self._subplot.map_screen_to_world(event) + if world is None: + self._tooltip.hide() + return + ch_idx = self._channel_at(float(world[0]), float(world[1])) + if ch_idx < 0: + self._tooltip.hide() + return + labels = self._config.channel_labels + label = labels[ch_idx] if labels and ch_idx < len(labels) else f"ch{ch_idx}" + text = f"{label}\n{self._buffer.n_retained} shown / {self._buffer.n_seen} total" + self._tooltip.setText(text) + self._tooltip.adjustSize() + x = int(getattr(event, "x", 0)) + 14 + y = int(getattr(event, "y", 0)) + 14 + max_x = self._fpl_widget.width() - self._tooltip.width() - 4 + max_y = self._fpl_widget.height() - self._tooltip.height() - 4 + self._tooltip.move(min(x, max_x), min(y, max_y)) + self._tooltip.show() + + def _on_pointer_leave(self, _event) -> None: + self._tooltip.hide() diff --git a/src/phosphor/trace_grid_buffer.py b/src/phosphor/trace_grid_buffer.py new file mode 100644 index 0000000..fe86d61 --- /dev/null +++ b/src/phosphor/trace_grid_buffer.py @@ -0,0 +1,202 @@ +"""Retention and statistics for a grid of per-channel waveforms. + +The CPU side of :class:`~phosphor.trace_grid.TraceGridWidget`, split out for the +same reason :class:`~phosphor.sweep_buffer.SweepBuffer` is: it is pure numpy, so +it can be exercised without a GPU, and that is where the behaviour that is easy +to get wrong actually lives. + +Two things it does that a naive buffer does not: + +**Writes into a ring.** A new waveform overwrites one slot. The alternative -- +rolling the whole buffer so the newest is always at index 0 -- copies every +retained sample on every arrival, which is affordable for an evoked potential +arriving once a second and not for spikes arriving hundreds of times a second. +Nothing is reordered, so a caller draws the slots where they lie and asks +:attr:`ages` which is newest. + +**Counts every waveform, not just the retained ones.** ``history`` is how many +to *draw*; the running mean and standard deviation come from every waveform ever +pushed. An evoked response is normally the average of hundreds of sweeps while +only a handful are worth overlaying, and tying the two together would force a +caller to retain hundreds of traces to average hundreds of traces. +""" + +from __future__ import annotations + +import logging + +import numpy as np + +logger = logging.getLogger(__name__) + +__all__ = ["TraceGridBuffer"] + + +class TraceGridBuffer: + """Ring of the last ``history`` waveforms per channel, plus running stats. + + :param n_channels: Channels drawn, one cell each. + :param n_samples: Samples per waveform. + :param history: How many recent waveforms to retain for drawing. + :param track_statistics: Accumulate the running mean and standard + deviation. Off for data with no meaningful average -- a stack of + action potentials from possibly-different units -- where it would cost + two float64 arrays and a pass over every arrival to compute a number + nobody displays. + """ + + def __init__( + self, + n_channels: int, + n_samples: int, + history: int = 10, + *, + track_statistics: bool = True, + ) -> None: + self.n_channels = int(n_channels) + self.n_samples = int(n_samples) + self._history = max(1, int(history)) + self.track_statistics = bool(track_statistics) + # Bumped when the drawn shape changes, so a renderer knows its graphics + # are the wrong size. Data changes bump `updates` instead, which is the + # cheap path: same shape, new numbers. + self.version = 0 + self.updates = 0 + self._allocate() + + # ------------------------------------------------------------------ + # Allocation + # ------------------------------------------------------------------ + + def _allocate(self) -> None: + self.traces = np.full((self._history, self.n_channels, self.n_samples), np.nan, dtype=np.float32) + self._cursor = 0 + self._n_retained = 0 + self._reset_statistics() + + def _reset_statistics(self) -> None: + # float64 because these accumulate without bound: a session's worth of + # sums in float32 loses the low bits of every later arrival. + shape = (self.n_channels, self.n_samples) + self._count = np.zeros(shape, dtype=np.float64) if self.track_statistics else None + self._sum = np.zeros(shape, dtype=np.float64) if self.track_statistics else None + self._sumsq = np.zeros(shape, dtype=np.float64) if self.track_statistics else None + self.n_seen = 0 + + # ------------------------------------------------------------------ + # Writing + # ------------------------------------------------------------------ + + def push(self, trace: np.ndarray) -> None: + """Retain one waveform per channel, shape ``(n_channels, n_samples)``. + + Channels with nothing to contribute pass NaN, which is retained as a gap + and left out of the statistics rather than counted as a zero. + """ + arr = np.asarray(trace, dtype=np.float32) + if arr.shape != (self.n_channels, self.n_samples): + raise ValueError(f"push expects {(self.n_channels, self.n_samples)}, got {arr.shape}") + + self.traces[self._cursor] = arr + self._cursor = (self._cursor + 1) % self._history + self._n_retained = min(self._n_retained + 1, self._history) + + if self.track_statistics: + finite = np.isfinite(arr) + contribution = np.where(finite, arr, 0.0).astype(np.float64) + self._count += finite + self._sum += contribution + self._sumsq += contribution * contribution + self.n_seen += 1 + self.updates += 1 + + def clear(self) -> None: + """Drop every retained waveform and forget the running statistics. + + Both together: a caller clears because what came before no longer + describes what comes next -- a new source, changed conditioning -- and + a mean carried across that boundary would average two different things. + """ + self.traces[:] = np.nan + self._cursor = 0 + self._n_retained = 0 + self._reset_statistics() + self.updates += 1 + + def set_history(self, history: int) -> None: + """Change how many waveforms are retained for drawing. + + Keeps as many of the newest as still fit, so growing the history does + not blank what is on screen and shrinking it drops the oldest. The + running statistics are untouched: how many are drawn has nothing to do + with how many are averaged. + """ + history = max(1, int(history)) + if history == self._history: + return + keep = min(self._n_retained, history) + newest = self.recent(keep) + self._history = history + self.traces = np.full((history, self.n_channels, self.n_samples), np.nan, dtype=np.float32) + if keep: + # Refill oldest-first so the cursor lands past the newest. + self.traces[:keep] = newest[::-1] + self._cursor = keep % history + self._n_retained = keep + self.version += 1 + + # ------------------------------------------------------------------ + # Reading + # ------------------------------------------------------------------ + + @property + def history(self) -> int: + return self._history + + @property + def n_retained(self) -> int: + """How many slots hold a waveform. Below ``history`` until it fills.""" + return self._n_retained + + @property + def ages(self) -> np.ndarray: + """Age of each slot, ``0`` for the newest, one entry per slot. + + Nothing is reordered on write, so this is how a renderer knows which + slot is which -- to fade older waveforms, or to skip empty ones. An age + below :attr:`n_retained` means the slot holds a waveform; at or above it + the slot has never been written. + """ + return (self._cursor - 1 - np.arange(self._history)) % self._history + + def recent(self, n: int | None = None) -> np.ndarray: + """The ``n`` newest waveforms, newest first, as a copy. + + Convenience for callers that want them in order and can afford the + copy. A renderer should prefer :attr:`traces` with :attr:`ages`, which + costs nothing. + """ + n = self._n_retained if n is None else min(int(n), self._n_retained) + if n <= 0: + return np.empty((0, self.n_channels, self.n_samples), dtype=np.float32) + idx = (self._cursor - 1 - np.arange(n)) % self._history + return self.traces[idx] + + def statistics(self) -> tuple[np.ndarray, np.ndarray] | None: + """Running ``(mean, std)`` over every waveform pushed, or ``None``. + + ``None`` when statistics are switched off or nothing has been pushed. + Elements no waveform contributed to are NaN rather than zero: no data is + not the same as a flat line, and drawing it as one invents a signal. + """ + if not self.track_statistics or self.n_seen == 0: + return None + with np.errstate(invalid="ignore", divide="ignore"): + n = self._count + mean = np.where(n > 0, self._sum / np.maximum(n, 1), np.nan) + # Var = E[x^2] - E[x]^2, clamped: with every sample equal the two + # terms cancel to a small negative, and a negative variance would + # make the square root NaN and take the band off screen. + var = np.maximum(self._sumsq / np.maximum(n, 1) - mean * mean, 0.0) + std = np.where(n > 1, np.sqrt(var), np.nan) + return mean.astype(np.float32), std.astype(np.float32) diff --git a/tests/test_channel_grid.py b/tests/test_channel_grid.py new file mode 100644 index 0000000..691d512 --- /dev/null +++ b/tests/test_channel_grid.py @@ -0,0 +1,151 @@ +"""Colour and tooltip policy for the per-channel value grid. + +The widget itself needs a GPU canvas, but the two decisions worth pinning do +not: which colour a channel gets, and what hovering it says. Both exist to keep +domain vocabulary out of this module -- a caller that decides what "good" means +supplies the colours and the word, and the grid renders them without an opinion. +""" + +import numpy as np +import pytest + +from phosphor.channel_grid import ChannelGridConfig, ChannelGridWidget + + +def make_grid(n_ch: int = 4, **config_kwargs) -> ChannelGridWidget: + """A grid with only the state the colour/tooltip logic reads. + + Built without ``__init__`` because everything else it constructs needs a + canvas; these paths are numpy and string formatting. + """ + import cmap as cmap_lib + + config = ChannelGridConfig(positions=np.zeros((n_ch, 2)), **config_kwargs) + w = ChannelGridWidget.__new__(ChannelGridWidget) + w._config = config + w._values_per_ch = np.full(n_ch, np.nan, dtype=np.float32) + w._cmap = cmap_lib.Colormap(config.cmap) + w._nan_rgba = np.asarray(config.nan_color, dtype=np.float32) + w._explicit_rgba = None + w._annotations = None + w._dirty = False + return w + + +# ---- colour ---------------------------------------------------------------- + + +def test_values_are_normalized_over_the_configured_range(): + w = make_grid(3, vmin=0.0, vmax=100.0, cmap="viridis") + w.push_data(np.array([0.0, 50.0, 100.0])) + + rgba = w._channel_rgba() + assert rgba.shape == (3, 4) + # Monotone through the colormap: the ends differ and the middle is between. + assert not np.allclose(rgba[0], rgba[2]) + assert not np.allclose(rgba[1], rgba[0]) + + +def test_a_channel_with_no_value_gets_the_nan_colour(): + w = make_grid(2, nan_color=(1.0, 0.0, 1.0, 1.0)) + w.push_data(np.array([5.0, np.nan])) + np.testing.assert_allclose(w._channel_rgba()[1], [1.0, 0.0, 1.0, 1.0]) + + +def test_values_outside_the_range_clamp_rather_than_wrap(): + """An out-of-range reading should saturate at the end of the scale, not + reappear at the other end looking healthy.""" + w = make_grid(4, vmin=0.0, vmax=10.0) + w.push_data(np.array([-50.0, 0.0, 10.0, 999.0])) + rgba = w._channel_rgba() + np.testing.assert_allclose(rgba[0], rgba[1]) + np.testing.assert_allclose(rgba[3], rgba[2]) + + +def test_a_degenerate_range_does_not_divide_by_zero(): + w = make_grid(2, vmin=5.0, vmax=5.0) + w.push_data(np.array([5.0, 5.0])) + assert np.all(np.isfinite(w._channel_rgba())) + + +def test_explicit_colours_override_the_colormap(): + """The categorical path: the caller classified the values and owns the + mapping, so the colormap must not get a second say.""" + w = make_grid(2) + w.push_data(np.array([1.0, 2.0])) + w.set_colors(np.array([[1.0, 0, 0, 1.0], [0, 1.0, 0, 1.0]])) + + np.testing.assert_allclose(w._channel_rgba(), [[1, 0, 0, 1], [0, 1, 0, 1]]) + # Including for a channel with no value: a caller that computes its own + # colours is the one that decides what missing looks like. + w.push_data(np.array([np.nan, 2.0])) + np.testing.assert_allclose(w._channel_rgba()[0], [1, 0, 0, 1]) + + +def test_clearing_explicit_colours_restores_the_colormap(): + w = make_grid(2, vmin=0.0, vmax=10.0) + w.push_data(np.array([0.0, 10.0])) + w.set_colors(np.zeros((2, 4))) + w.set_colors(None) + assert not np.allclose(w._channel_rgba()[0], w._channel_rgba()[1]) + + +def test_wrongly_sized_colours_are_refused_not_applied(): + """Silently colouring the wrong channels is worse than not colouring.""" + w = make_grid(4) + w.set_colors(np.ones((3, 4))) + assert w._explicit_rgba is None + + +# ---- tooltip --------------------------------------------------------------- + + +def test_tooltip_names_the_channel_and_its_value(): + w = make_grid(2, channel_labels=["A-1", "A-2"], value_unit="kOhm") + w.push_data(np.array([12.34, np.nan])) + assert w._tooltip_text(0) == "A-1\n12.3 kOhm" + + +def test_tooltip_falls_back_to_an_index_when_unlabelled(): + w = make_grid(2) + w.push_data(np.array([1.0, 2.0])) + assert w._tooltip_text(1).startswith("ch1\n") + + +def test_tooltip_omits_a_unit_that_was_never_set(): + w = make_grid(1) + w.push_data(np.array([7.0])) + assert w._tooltip_text(0) == "ch0\n7.0" + + +def test_annotations_appear_beside_the_value(): + w = make_grid(2, channel_labels=["A-1", "A-2"], value_unit="kOhm") + w.push_data(np.array([12.0, 900.0])) + w.set_annotations(["good", "open"]) + assert w._tooltip_text(1) == "A-2\n900.0 kOhm (open)" + + +def test_an_annotation_stands_in_for_a_missing_value(): + """A channel with no reading still has something worth saying about it.""" + w = make_grid(1, channel_labels=["A-1"]) + w.push_data(np.array([np.nan])) + w.set_annotations(["unmeasured"]) + assert w._tooltip_text(0) == "A-1\nunmeasured" + + +def test_a_missing_value_says_so_without_an_annotation(): + w = make_grid(1) + assert w._tooltip_text(0) == "ch0\nno value" + + +def test_wrongly_sized_annotations_are_refused(): + w = make_grid(3) + w.set_annotations(["a", "b"]) + assert w._annotations is None + + +@pytest.mark.parametrize("value", [0.0, -1.5, 1e6]) +def test_tooltip_formats_any_finite_value(value): + w = make_grid(1, value_unit="uV") + w.push_data(np.array([value])) + assert w._tooltip_text(0).startswith("ch0\n") diff --git a/tests/test_decimate.py b/tests/test_decimate.py new file mode 100644 index 0000000..8794ffe --- /dev/null +++ b/tests/test_decimate.py @@ -0,0 +1,82 @@ +"""Tests for the min/max envelope decimation primitive. + +Pure-array display decimation shared by the plotting widgets (evoked trace grid +today). These assert the properties the widgets rely on: peak amplitude is +preserved (the reason min/max beats stride for evoked P/N peaks), output shapes +match across leading dims, empty (all-NaN) history slots survive as NaN, the +x-envelope stays monotone, and small waveforms pass through untouched. +""" + +from __future__ import annotations + +import numpy as np + +from phosphor.decimate import ( + minmax_decimate, + minmax_decimate_x, + plan_minmax_decimation, +) + + +def test_preserves_peaks_where_stride_would_clip(): + n = 15000 + sig = np.zeros(n, dtype=np.float32) + # Sharp P/N peaks on adjacent samples inside one bucket. + sig[7001] = 120.0 + sig[7002] = -90.0 + plan = plan_minmax_decimation(n, 1500) + assert plan.active and plan.n_out == 1500 + + dec = minmax_decimate(sig, plan) + assert dec.shape == (1500,) + # Both extremes retained. + assert dec.max() == 120.0 + assert dec.min() == -90.0 + # Plain stride lands between the peaks and misses at least one. + stride = sig[:: plan.bucket_size] + assert not (stride.max() == 120.0 and stride.min() == -90.0) + + +def test_multidim_shape_and_x_alignment(): + n = 15000 + plan = plan_minmax_decimation(n, 1500) + arr = np.random.default_rng(0).standard_normal((11, 128, n)).astype(np.float32) + dec = minmax_decimate(arr, plan) + assert dec.shape == (11, 128, plan.n_out) + + x = np.broadcast_to(np.linspace(0.0, 0.5, n, dtype=np.float32), (128, n)) + xd = minmax_decimate_x(x, plan) + assert xd.shape == (128, plan.n_out) + # Time axis stays monotone non-decreasing (each bucket center repeated for + # its min/max pair), so the drawn polyline never doubles back. + assert np.all(np.diff(xd, axis=-1) >= 0) + + +def test_all_nan_bucket_stays_nan(): + n = 1200 + plan = plan_minmax_decimation(n, 200) + assert plan.active + row = np.full(n, np.nan, dtype=np.float32) + dec = minmax_decimate(row, plan) + assert dec.shape == (plan.n_out,) + assert np.all(np.isnan(dec)) + + +def test_passthrough_when_within_budget(): + n = 800 + plan = plan_minmax_decimation(n, 1500) + assert not plan.active + sig = np.arange(n, dtype=np.float32) + out = minmax_decimate(sig, plan) + # Untouched (same object, same shape) so full resolution reaches the GPU. + assert out.shape == (n,) + assert np.array_equal(out, sig) + x = np.broadcast_to(sig, (4, n)) + assert np.array_equal(minmax_decimate_x(x, plan), x) + + +def test_dtype_preserved(): + n = 4000 + plan = plan_minmax_decimation(n, 500) + sig = np.random.default_rng(1).standard_normal(n).astype(np.float32) + assert minmax_decimate(sig, plan).dtype == np.float32 diff --git a/tests/test_grid_layout.py b/tests/test_grid_layout.py new file mode 100644 index 0000000..1a8a3a5 --- /dev/null +++ b/tests/test_grid_layout.py @@ -0,0 +1,162 @@ +"""Cell geometry for per-channel grids. + +These are the calculations that decide where a channel's cell lands, and they +are all in the business of degrading gracefully: real electrode geometry is +often missing, partial, or reported by two devices that each numbered their +electrodes from the same origin. Getting that wrong does not raise -- it draws +every channel on top of every other, or silently off screen, which is why the +fallbacks are pinned here rather than left to the widgets. +""" + +import numpy as np +import pytest + +from phosphor.grid_layout import ( + build_quad_mesh_arrays, + infer_pitch, + resolve_cell_geometry, + tile_by_group, + tiled_grid_positions, +) + +# ---- fallback layout -------------------------------------------------------- + + +def test_tiled_positions_are_square_ish_and_unique(): + pos = tiled_grid_positions(10) + assert pos.shape == (10, 2) + assert len({tuple(p) for p in pos}) == 10, "channels must not share a tile" + # ceil(sqrt(10)) == 4 columns, so x stays under 4 and y under 3. + assert pos[:, 0].max() < 4 and pos[:, 1].max() < 3 + + +def test_tiled_positions_handles_no_channels(): + assert tiled_grid_positions(0).shape == (0, 2) + + +# ---- pitch ------------------------------------------------------------------ + + +def test_pitch_is_the_smallest_positive_step(): + xs = np.array([0.0, 400.0, 800.0]) + ys = np.array([0.0, 400.0, 800.0]) + assert infer_pitch(xs, ys) == pytest.approx(400.0) + + +def test_pitch_of_degenerate_geometry_is_one(): + """One column of channels has no x spacing to measure; a zero pitch would + collapse every cell to nothing.""" + xs = np.zeros(4) + ys = np.array([0.0, 2.0, 4.0, 6.0]) + assert infer_pitch(xs, ys) == pytest.approx(2.0) + assert infer_pitch(np.zeros(4), np.zeros(4)) == 1.0 + + +# ---- grouping --------------------------------------------------------------- + + +def test_groups_sharing_a_coordinate_range_are_fanned_out(): + """Two devices that each number electrodes from their own origin would + otherwise render exactly on top of one another.""" + block = np.array([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]], dtype=np.float32) + positions = np.vstack([block, block]) + groups = np.array([0, 0, 0, 0, 1, 1, 1, 1]) + + tiled = tile_by_group(positions, groups) + + first, second = tiled[:4, 0], tiled[4:, 0] + assert second.min() > first.max(), "the second block must start clear of the first" + # y is untouched: only the blocks' horizontal placement changes. + np.testing.assert_array_equal(tiled[:, 1], positions[:, 1]) + # Each block keeps its own internal spacing. + assert np.ptp(first) == pytest.approx(np.ptp(block[:, 0])) + assert np.ptp(second) == pytest.approx(np.ptp(block[:, 0])) + + +def test_a_single_group_is_left_alone(): + """The common case -- one device, or an all-zero group field because no + geometry was loaded -- must not have its layout rewritten.""" + positions = np.array([[0.0, 0.0], [1.0, 0.0]], dtype=np.float32) + np.testing.assert_array_equal(tile_by_group(positions, np.zeros(2)), positions) + + +def test_mismatched_group_ids_are_ignored_rather_than_raising(): + """A stale group array is a bad layout, not a crash: the plot is a + diagnostic tool and is most wanted when the metadata is wrong.""" + positions = np.array([[0.0, 0.0], [1.0, 0.0]], dtype=np.float32) + np.testing.assert_array_equal(tile_by_group(positions, np.array([0, 1, 2])), positions) + + +# ---- cell geometry ---------------------------------------------------------- + + +def test_cells_default_to_the_pitch_so_they_tile(): + positions = np.array([[0.0, 0.0], [400.0, 0.0], [0.0, 400.0], [400.0, 400.0]]) + rects, centers = resolve_cell_geometry(positions, None, 4, invert_y=False) + + np.testing.assert_allclose(rects[:, 2], 400.0) + np.testing.assert_allclose(centers[0], [200.0, 200.0]) + + +def test_a_channel_smaller_than_the_pitch_keeps_its_size(): + """An electrode reported as smaller than the array pitch should draw + smaller, with a gap -- not fill its cell.""" + positions = np.array([[0.0, 0.0], [400.0, 0.0]]) + rects, _ = resolve_cell_geometry(positions, np.array([100.0, 400.0]), 2, invert_y=False) + assert rects[0, 2] == pytest.approx(100.0) + assert rects[1, 2] == pytest.approx(400.0) + + +def test_a_missing_size_falls_back_to_the_pitch(): + """Zero is what an unpopulated size field holds, and a zero-side cell is + invisible.""" + positions = np.array([[0.0, 0.0], [400.0, 0.0]]) + rects, _ = resolve_cell_geometry(positions, np.array([0.0, 400.0]), 2, invert_y=False) + assert rects[0, 2] == pytest.approx(400.0) + + +def test_invert_y_flips_the_layout_without_reshaping_cells(): + positions = np.array([[0.0, 0.0], [0.0, 400.0]]) + upright, up_centers = resolve_cell_geometry(positions, None, 2, invert_y=False) + flipped, flip_centers = resolve_cell_geometry(positions, None, 2, invert_y=True) + + # Row order reverses on screen... + assert up_centers[0, 1] < up_centers[1, 1] + assert flip_centers[0, 1] > flip_centers[1, 1] + # ...while cells stay the same size and keep their x. + np.testing.assert_allclose(flipped[:, 2], upright[:, 2]) + np.testing.assert_allclose(flipped[:, 0], upright[:, 0]) + + +def test_every_channel_at_one_point_lays_out_sequentially(): + """What an unmapped source looks like: no geometry, so every channel + defaulted to the same coordinate. Drawn verbatim it is one cell deep.""" + positions = np.zeros((9, 2)) + rects, centers = resolve_cell_geometry(positions, None, 9, invert_y=False) + + assert len({tuple(c) for c in centers}) == 9, "cells must not stack" + np.testing.assert_allclose(rects[:, 2], 1.0) + + +# ---- mesh ------------------------------------------------------------------- + + +def test_quad_mesh_covers_each_cell_with_two_triangles(): + rects = np.array([[0.0, 0.0, 2.0], [10.0, 0.0, 4.0]]) + positions, indices = build_quad_mesh_arrays(rects) + + assert positions.shape == (8, 3) + assert indices.shape == (4, 3) + # First quad spans its own rect, corner to corner. + np.testing.assert_allclose(positions[:4, 0].min(), 0.0) + np.testing.assert_allclose(positions[:4, 0].max(), 2.0) + # Second quad is offset and larger, and indexes its own vertices only. + np.testing.assert_allclose(positions[4:, 0].min(), 10.0) + np.testing.assert_allclose(positions[4:, 0].max(), 14.0) + assert indices[2:].min() >= 4 + + +def test_quad_mesh_of_no_cells_is_empty_not_malformed(): + positions, indices = build_quad_mesh_arrays(np.zeros((0, 3))) + assert positions.shape == (0, 3) + assert indices.shape == (0, 3) diff --git a/tests/test_trace_grid.py b/tests/test_trace_grid.py new file mode 100644 index 0000000..91a5324 --- /dev/null +++ b/tests/test_trace_grid.py @@ -0,0 +1,318 @@ +"""The trace grid's drawing decisions, minus the canvas. + +Three of them matter enough to pin: which slot is drawn how brightly (the ring +does not reorder, so brightness has to follow age), when the graphics have to be +recreated rather than written into (getting this wrong is what made the widget +tear down every retained waveform to draw one new one), and the shape of the +error band. +""" + +import numpy as np +import pytest + +from phosphor.decimate import plan_minmax_decimation +from phosphor.trace_grid import TraceGridConfig, TraceGridWidget +from phosphor.trace_grid_buffer import TraceGridBuffer + + +def make_widget(n_ch=2, n_samples=4, history=3, **config_kwargs) -> TraceGridWidget: + """A widget holding only what the drawing helpers read. + + Built without ``__init__`` because the rest of it needs a GPU canvas. + """ + config = TraceGridConfig(positions=np.zeros((n_ch, 2)), n_samples=n_samples, history=history, **config_kwargs) + w = TraceGridWidget.__new__(TraceGridWidget) + w._config = config + w._n_ch = n_ch + w._n_samples = n_samples + w._show_individual = config.show_individual + w._show_mean = config.show_mean + w._show_error = config.show_error and config.track_statistics + w._buffer = TraceGridBuffer(n_ch, n_samples, history, track_statistics=config.track_statistics) + w._dec_plan = plan_minmax_decimation(n_samples, 10_000) # inactive at these sizes + w._x_line_dec = np.tile(np.arange(n_samples, dtype=np.float32), (n_ch, 1)) + w._indiv_ml = w._mean_ml = w._error_ml = None + w._graphics_version = -1 + # _map_y is affine per channel; identity keeps these tests about layout. + w._map_y = lambda a: np.asarray(a, dtype=np.float32) + return w + + +def wave(value, n_ch=2, n_samples=4): + return np.full((n_ch, n_samples), value, dtype=np.float32) + + +# ---- age-graded brightness -------------------------------------------------- + + +def test_the_newest_waveform_is_the_brightest(): + w = make_widget(history=3) + for v in (1.0, 2.0, 3.0): + w._buffer.push(wave(v)) + + alpha = w._slot_alpha() + ages = w._buffer.ages + assert alpha[ages == 0] > alpha[ages == 1] > alpha[ages == 2] + + +def test_brightness_follows_the_slot_that_holds_the_newest(): + """Nothing is reordered on write, so after a wrap the brightest slot is not + slot 0. Keying on position instead of age would light up a stale waveform.""" + w = make_widget(history=3) + for v in (1.0, 2.0, 3.0, 4.0): # wraps: newest lands back in slot 0 + w._buffer.push(wave(v)) + + brightest = int(np.argmax(w._slot_alpha())) + assert w._buffer.traces[brightest][0, 0] == pytest.approx(4.0) + + +def test_slots_never_written_are_fully_transparent(): + w = make_widget(history=4) + w._buffer.push(wave(1.0)) + + alpha = w._slot_alpha() + assert (alpha > 0).sum() == 1 + np.testing.assert_allclose(alpha[w._buffer.ages >= w._buffer.n_retained], 0.0) + + +def test_fading_can_be_switched_off(): + """Every retained waveform's age changes when one arrives, so fading costs + a colour rewrite per arrival -- worth it at evoked rates, not spike rates.""" + w = make_widget(history=3, age_fade=False) + for v in (1.0, 2.0, 3.0): + w._buffer.push(wave(v)) + + alpha = w._slot_alpha() + assert len(set(alpha.tolist())) == 1, "all retained waveforms should look alike" + + +def test_slot_rows_are_contiguous_so_one_arrival_writes_one_block(): + w = make_widget(n_ch=8, history=4) + rows = w._slot_rows(2) + assert (rows.start, rows.stop) == (16, 24) + + +# ---- when to recreate the graphics ----------------------------------------- + + +def test_a_new_waveform_does_not_change_the_graphics_shape(): + """The whole point of the split: numbers change, shapes do not.""" + w = make_widget(history=3) + w._buffer.push(wave(1.0)) + shape = w._graphics_shape() + + w._buffer.push(wave(2.0)) + assert w._graphics_shape() == shape + + +@pytest.mark.parametrize( + "change", + [ + lambda w: w._buffer.set_history(9), + lambda w: setattr(w, "_show_mean", not w._show_mean), + lambda w: setattr(w, "_show_error", not w._show_error), + lambda w: setattr(w, "_show_individual", not w._show_individual), + ], +) +def test_anything_that_adds_or_removes_lines_does_change_it(change): + w = make_widget(history=3) + # Two, so there is a standard deviation and the error toggle has an effect. + w._buffer.push(wave(1.0)) + w._buffer.push(wave(3.0)) + shape = w._graphics_shape() + + change(w) + assert w._graphics_shape() != shape + + +# ---- mean and error band ---------------------------------------------------- + + +def test_the_band_is_two_curves_per_channel(): + w = make_widget(n_ch=2, history=4, show_error=True) + for v in (1.0, 3.0): + w._buffer.push(wave(v)) + + mean_pos, band_pos = w._summary_positions() + assert mean_pos.shape[0] == 2, "one mean line per channel" + assert band_pos.shape[0] == 4, "a lower and an upper line per channel" + + +def test_the_band_brackets_the_mean(): + w = make_widget(n_ch=1, history=4, show_error=True) + for v in (1.0, 3.0): + w._buffer.push(wave(v, n_ch=1)) + + mean_pos, band_pos = w._summary_positions() + mean_y = mean_pos[0, :, 1] + lower, upper = band_pos[0, :, 1], band_pos[1, :, 1] + assert np.all(lower < mean_y) and np.all(upper > mean_y) + + +def test_the_band_keeps_each_curve_on_its_own_channel_x(): + """Two curves per channel means x tiles rather than broadcasts; getting it + wrong draws the band across the wrong cells.""" + w = make_widget(n_ch=2, history=4, show_error=True) + for v in (1.0, 3.0): + w._buffer.push(wave(v)) + + _, band_pos = w._summary_positions() + np.testing.assert_allclose(band_pos[0, :, 0], w._x_line_dec[0]) + np.testing.assert_allclose(band_pos[2, :, 0], w._x_line_dec[0]) + + +def test_no_band_until_there_is_a_spread_to_draw(): + """One waveform has a mean but no standard deviation, and a band of NaN + would take the view off screen.""" + w = make_widget(history=4, show_error=True) + w._buffer.push(wave(1.0)) + + mean_pos, band_pos = w._summary_positions() + assert mean_pos is not None + assert band_pos is None + + +def test_no_summary_at_all_before_anything_arrives(): + w = make_widget(show_error=True) + assert w._summary_positions() == (None, None) + + +def test_the_band_is_refused_when_statistics_are_off(): + """Asking for a band around a mean that is never computed should resolve to + 'no band', not to a crash at draw time.""" + w = make_widget(history=3, show_error=True, track_statistics=False) + for v in (1.0, 3.0): + w._buffer.push(wave(v)) + + assert w._show_error is False + assert w._summary_positions() == (None, None) + + +def test_the_mean_spans_more_waveforms_than_are_drawn(): + """What the running accumulator buys: an evoked average of everything seen, + overlaid with the handful worth looking at.""" + w = make_widget(n_ch=1, history=2) + for v in (1.0, 2.0, 3.0, 4.0, 5.0): + w._buffer.push(wave(v, n_ch=1)) + + mean_pos, _ = w._summary_positions() + np.testing.assert_allclose(mean_pos[0, :, 1], 3.0) # mean of 1..5, not of 4..5 + + +# ---- graphics actually get created ------------------------------------------ + + +class FakeColors: + """fastplotlib takes one colour per line, then holds one per vertex.""" + + def __init__(self, n_lines, n_points): + self.value = np.zeros((n_lines * n_points, 4), dtype=np.float32) + + def __setitem__(self, key, value): + self.value[key] = value + + +class FakeGraphic: + def __init__(self, data): + self.data = np.asarray(data).copy() + self.colors = FakeColors(self.data.shape[0], self.data.shape[1]) + self.visible = True + + +class FakeSubplot: + """Records what a renderer would have been asked to draw.""" + + def __init__(self): + self.created: list[FakeGraphic] = [] + self.deleted: list[FakeGraphic] = [] + + def add_multi_line(self, data, colors=None, thickness=None, **_): + g = FakeGraphic(data) + self.created.append(g) + return g + + def delete_graphic(self, graphic): + self.deleted.append(graphic) + + +def drivable(**kwargs) -> TraceGridWidget: + w = make_widget(**kwargs) + w._subplot = FakeSubplot() + return w + + +def test_the_first_waveform_creates_the_graphics(): + """The regression that showed up as a grid of empty cells with epochs + arriving: the widget is built before any data, so the first frame rebuilds + against an empty buffer. If that state is indistinguishable from a + populated one, every later arrival takes the in-place path, finds nothing + to write into, and draws nothing -- for as long as the app runs. + """ + w = drivable(history=3) + + w._refresh_lines() # a frame before any data, as happens on startup + assert w._subplot.created == [], "nothing to draw yet" + + w._buffer.push(wave(1.0)) + w._refresh_lines() + + assert w._subplot.created, "the first waveform must bring its graphics with it" + assert w._indiv_ml is not None + + +def test_later_waveforms_write_in_place_instead_of_recreating(): + """The other half: once the graphics exist, arrivals must not rebuild.""" + w = drivable(history=3) + w._buffer.push(wave(1.0)) + w._refresh_lines() + created = len(w._subplot.created) + + for v in (2.0, 3.0, 4.0): + w._buffer.push(wave(v)) + w._refresh_lines() + + assert len(w._subplot.created) == created, "arrivals should not recreate graphics" + assert w._subplot.deleted == [] + + +def test_the_newest_waveform_reaches_the_graphic(): + """An in-place path that silently wrote nowhere would look identical to a + working one until you looked at the screen.""" + w = drivable(history=3, show_mean=False) + w._buffer.push(wave(1.0)) + w._refresh_lines() + + w._buffer.push(wave(7.0)) + w._refresh_lines() + + ys = w._indiv_ml.data[..., 1] + assert np.isclose(ys, 7.0).any(), "the value just pushed should be in the graphic" + + +def test_the_error_band_appears_once_there_is_a_spread(): + """It cannot be drawn from one waveform, so it arrives a frame later than + the mean and needs its own place in the shape key.""" + w = drivable(history=4, show_error=True) + + w._buffer.push(wave(1.0)) + w._refresh_lines() + assert w._error_ml is None + + w._buffer.push(wave(3.0)) + w._refresh_lines() + assert w._error_ml is not None + + +def test_clearing_and_refilling_brings_the_graphics_back(): + """clear() drops to the empty state; refilling has to leave it again.""" + w = drivable(history=3) + w._buffer.push(wave(1.0)) + w._refresh_lines() + + w._buffer.clear() + w._refresh_lines() + assert w._indiv_ml is None + + w._buffer.push(wave(2.0)) + w._refresh_lines() + assert w._indiv_ml is not None diff --git a/tests/test_trace_grid_buffer.py b/tests/test_trace_grid_buffer.py new file mode 100644 index 0000000..2765727 --- /dev/null +++ b/tests/test_trace_grid_buffer.py @@ -0,0 +1,271 @@ +"""Retention and running statistics for the trace grid. + +Pure numpy, so all of it runs headlessly. Two properties carry the weight: +a waveform arriving must not touch the ones already retained (the ring), and +what gets averaged must not be limited by what gets drawn (the accumulator). +Both are the difference between a widget that works for one waveform a second +and one that works for hundreds. +""" + +import numpy as np +import pytest + +from phosphor.trace_grid_buffer import TraceGridBuffer + + +def make_buffer(n_channels=2, n_samples=4, history=3, **kwargs) -> TraceGridBuffer: + return TraceGridBuffer(n_channels, n_samples, history, **kwargs) + + +def wave(value: float, n_channels=2, n_samples=4) -> np.ndarray: + """A waveform whose every sample is ``value``, so a slot's identity is + readable straight off its contents.""" + return np.full((n_channels, n_samples), value, dtype=np.float32) + + +# ---- the ring --------------------------------------------------------------- + + +def test_a_new_waveform_writes_one_slot_and_leaves_the_rest(): + """The whole point of the ring: arrivals cost one slot, not a copy of + everything retained.""" + buf = make_buffer(history=3) + buf.push(wave(1.0)) + before = buf.traces.copy() + + buf.push(wave(2.0)) + + # equal_nan, or the never-written slots read as "changed": NaN != NaN. + changed = [i for i in range(3) if not np.array_equal(buf.traces[i], before[i], equal_nan=True)] + assert changed == [1], "exactly one slot should differ" + + +def test_ages_identify_the_newest_without_reordering(): + buf = make_buffer(history=3) + for v in (1.0, 2.0, 3.0): + buf.push(wave(v)) + + ages = buf.ages + newest_slot = int(np.flatnonzero(ages == 0)[0]) + assert buf.traces[newest_slot][0, 0] == pytest.approx(3.0) + oldest_slot = int(np.flatnonzero(ages == 2)[0]) + assert buf.traces[oldest_slot][0, 0] == pytest.approx(1.0) + + +def test_an_age_at_or_above_the_retained_count_marks_an_empty_slot(): + """The invariant a renderer uses to skip slots that were never written, + without needing a separate mask.""" + buf = make_buffer(history=4) + buf.push(wave(1.0)) + buf.push(wave(2.0)) + + filled = buf.ages < buf.n_retained + assert filled.sum() == 2 + for slot in np.flatnonzero(~filled): + assert np.all(np.isnan(buf.traces[slot])), "an unfilled slot must stay NaN" + + +def test_the_ring_wraps_and_keeps_the_newest(): + buf = make_buffer(history=3) + for v in (1.0, 2.0, 3.0, 4.0, 5.0): + buf.push(wave(v)) + + assert buf.n_retained == 3 + retained = {float(buf.traces[i][0, 0]) for i in range(3)} + assert retained == {3.0, 4.0, 5.0} + np.testing.assert_allclose(buf.recent(1)[0], wave(5.0)) + + +def test_recent_returns_newest_first(): + buf = make_buffer(history=3) + for v in (1.0, 2.0, 3.0): + buf.push(wave(v)) + np.testing.assert_allclose([r[0, 0] for r in buf.recent()], [3.0, 2.0, 1.0]) + + +def test_recent_of_an_empty_buffer_is_empty_not_an_error(): + assert make_buffer().recent().shape == (0, 2, 4) + + +def test_a_wrongly_shaped_waveform_is_refused(): + """Retaining a mis-shaped trace would silently misalign every channel.""" + buf = make_buffer(n_channels=2, n_samples=4) + with pytest.raises(ValueError, match=r"expects \(2, 4\)"): + buf.push(np.zeros((4, 2), dtype=np.float32)) + + +# ---- history changes -------------------------------------------------------- + + +def test_growing_the_history_keeps_what_was_on_screen(): + buf = make_buffer(history=2) + buf.push(wave(1.0)) + buf.push(wave(2.0)) + + buf.set_history(5) + + assert buf.n_retained == 2 + np.testing.assert_allclose([r[0, 0] for r in buf.recent()], [2.0, 1.0]) + + +def test_shrinking_the_history_drops_the_oldest(): + buf = make_buffer(history=4) + for v in (1.0, 2.0, 3.0, 4.0): + buf.push(wave(v)) + + buf.set_history(2) + + assert buf.n_retained == 2 + np.testing.assert_allclose([r[0, 0] for r in buf.recent()], [4.0, 3.0]) + + +def test_a_resized_ring_keeps_accepting_in_order(): + """The cursor has to land past the newest, or the next push overwrites it.""" + buf = make_buffer(history=4) + for v in (1.0, 2.0, 3.0): + buf.push(wave(v)) + buf.set_history(3) + buf.push(wave(4.0)) + + np.testing.assert_allclose([r[0, 0] for r in buf.recent()], [4.0, 3.0, 2.0]) + + +def test_resizing_the_history_does_not_disturb_the_statistics(): + """How many are drawn is not how many are averaged.""" + buf = make_buffer(history=2) + for v in (1.0, 2.0, 3.0, 4.0): + buf.push(wave(v)) + mean_before, _ = buf.statistics() + + buf.set_history(8) + + mean_after, _ = buf.statistics() + np.testing.assert_allclose(mean_after, mean_before) + assert buf.n_seen == 4 + + +def test_changing_the_history_bumps_the_version_but_a_push_does_not(): + """A renderer rebuilds its graphics on a shape change and writes in place + otherwise; conflating the two is what makes a plot blank on every arrival.""" + buf = make_buffer(history=2) + version, updates = buf.version, buf.updates + + buf.push(wave(1.0)) + assert buf.version == version, "a new waveform must not force a rebuild" + assert buf.updates > updates + + buf.set_history(4) + assert buf.version > version + + +# ---- running statistics ----------------------------------------------------- + + +def test_the_mean_spans_every_waveform_not_just_the_retained_ones(): + """The reason the accumulator exists: an evoked response averages hundreds + of sweeps while only a few are worth overlaying.""" + buf = make_buffer(history=2) + for v in (1.0, 2.0, 3.0, 4.0, 5.0): + buf.push(wave(v)) + + mean, _ = buf.statistics() + assert buf.n_retained == 2 + np.testing.assert_allclose(mean, wave(3.0)) # mean of 1..5, not of 4..5 + + +def test_the_standard_deviation_is_the_population_spread(): + buf = make_buffer(history=10) + for v in (2.0, 4.0, 4.0, 4.0, 5.0, 5.0, 7.0, 9.0): + buf.push(wave(v)) + + _, std = buf.statistics() + np.testing.assert_allclose(std, wave(2.0), rtol=1e-5) + + +def test_identical_waveforms_give_a_zero_spread_not_a_nan(): + """Var = E[x^2] - E[x]^2 cancels to a small negative here, and an unclamped + square root would take the band off screen.""" + buf = make_buffer(history=4) + for _ in range(4): + buf.push(wave(1000.0)) + + _, std = buf.statistics() + np.testing.assert_allclose(std, 0.0, atol=1e-3) + + +def test_a_single_waveform_has_a_mean_but_no_spread(): + buf = make_buffer() + buf.push(wave(3.0)) + mean, std = buf.statistics() + np.testing.assert_allclose(mean, wave(3.0)) + assert np.all(np.isnan(std)), "one sample has no spread to report" + + +def test_channels_that_contributed_nothing_stay_nan(): + """A channel with no data is not a flat line at zero; drawing it as one + invents a signal that was never recorded.""" + buf = make_buffer(n_channels=2, history=4) + for v in (1.0, 3.0): + trace = wave(v) + trace[1] = np.nan # channel 1 never contributes + buf.push(trace) + + mean, std = buf.statistics() + np.testing.assert_allclose(mean[0], 2.0) + assert np.all(np.isnan(mean[1])) + assert np.all(np.isnan(std[1])) + + +def test_a_missing_channel_does_not_drag_the_others_mean_down(): + """NaN counted as zero is the easy bug, and it looks plausible.""" + buf = make_buffer(n_channels=1, history=4) + for v in (10.0, np.nan, 20.0): + buf.push(wave(v, n_channels=1)) + + mean, _ = buf.statistics() + np.testing.assert_allclose(mean, 15.0) # not 10.0, which averaging in a 0 would give + + +def test_statistics_are_absent_before_anything_arrives(): + assert make_buffer().statistics() is None + + +def test_statistics_can_be_switched_off(): + """A stack of action potentials from possibly-different units has no + meaningful average, so nothing should be paid to compute one.""" + buf = make_buffer(track_statistics=False) + buf.push(wave(1.0)) + assert buf.statistics() is None + assert buf.n_retained == 1 + + +def test_clearing_forgets_the_average_as_well_as_the_traces(): + """A caller clears at a boundary -- new source, changed conditioning -- and + a mean carried across it would average two different things.""" + buf = make_buffer(history=3) + for v in (1.0, 2.0): + buf.push(wave(v)) + + buf.clear() + + assert buf.n_retained == 0 and buf.n_seen == 0 + assert buf.statistics() is None + assert np.all(np.isnan(buf.traces)) + + buf.push(wave(9.0)) + mean, _ = buf.statistics() + np.testing.assert_allclose(mean, wave(9.0)) + + +def test_the_accumulator_holds_up_over_many_arrivals(): + """float32 sums lose the low bits of later arrivals; the accumulators are + float64 for exactly this.""" + buf = make_buffer(n_channels=1, n_samples=1, history=2) + rng = np.random.default_rng(0) + values = rng.standard_normal(5000).astype(np.float32) * 100.0 + 1000.0 + for v in values: + buf.push(np.full((1, 1), v, dtype=np.float32)) + + mean, std = buf.statistics() + np.testing.assert_allclose(mean[0, 0], values.mean(), rtol=1e-4) + np.testing.assert_allclose(std[0, 0], values.std(), rtol=1e-3)