fix(meg): correct M100 measurement (standards-only events + in-window peak search) + study_status

matilde_plugin/engine/meg_study.py

@@ -57,12 +57,23 @@
 DEFAULT_BOUNDS: Dict[str, Any] = {
     "run": 1,            # one run only
     "max_channels": 60,  # subset of gradiometers, not the full array
-    "crop_tmax": 30.0,   # seconds — a stimulus-locked slice, never the whole recording
+    # seconds — a stimulus-locked slice, never the whole recording. Raised from
+    # 30 to 90 s so the standards-only average has enough trials for a clean
+    # peak (a 30 s crop yielded only ~16 standard tones). Still bounded: after
+    # the channel-subset + resample the loaded array stays well within budget
+    # (~60 grad ch x 90 s x 200 Hz ~ a few MB), and per-step checkpointing means
+    # an OOM still resumes from the last completed step.
+    "crop_tmax": 90.0,
     "resample_hz": 200,  # downsample to shrink the in-memory array
     "l_freq": 1.0,       # band-pass low edge (Hz)
     "h_freq": 40.0,      # band-pass high edge (Hz)
     "tmin": -0.1,        # epoch window start (s) relative to event
     "tmax": 0.3,         # epoch window end (s) relative to event
+    # Which trigger code to epoch. None -> auto: the most-frequent event code,
+    # which for an oddball auditory paradigm is the STANDARD tone. Epoching ALL
+    # triggers (standards + deviants + button presses) muddies the average and
+    # was part of the original mis-measurement.
+    "event_id": None,
 }
 
 # How far outside the expected window a measured peak must fall before the finding
@@ -71,6 +82,32 @@
 # inconclusive (not clear-cut either way).
 _REFUTE_MARGIN_MS = 50.0
 
+# Small symmetric slack (ms) added around the expected window when searching for
+# the peak, so a peak landing right at an edge is still captured. Deliberately
+# tiny: a wide search (the old +/-100 ms) grabbed the large early stimulus
+# artifact (~15 ms) instead of the in-window auditory peak (~100-120 ms).
+_PEAK_SEARCH_MARGIN_MS = 10.0
+
+
+def _peak_search_bounds(window_ms: Tuple[float, float],
+                        times_lo: float, times_hi: float) -> Tuple[float, float]:
+    """Compute the (tmin, tmax) seconds to search for the evoked peak.
+
+    Pure and mne-free so it is unit-testable without the scientific stack. The
+    search stays WITHIN the expected ``window_ms`` plus only a small symmetric
+    margin (``_PEAK_SEARCH_MARGIN_MS``), then is clamped to the available sample
+    times ``[times_lo, times_hi]`` (seconds).
+
+    This is the regression boundary for the M100 bug: a wide (+/-100 ms) search
+    over an 80-120 ms window reached back to the ~15 ms stimulus artifact and
+    reported it as the peak. Keeping the search in-window prevents that.
+    """
+    lo_s = window_ms[0] / 1000.0 - _PEAK_SEARCH_MARGIN_MS / 1000.0
+    hi_s = window_ms[1] / 1000.0 + _PEAK_SEARCH_MARGIN_MS / 1000.0
+    lo = max(times_lo, lo_s)
+    hi = min(times_hi, hi_s)
+    return lo, hi
+
 
 class MegIO(abc.ABC):
     """Structural contract for the MEG I/O backend the steps call through.
@@ -161,18 +198,38 @@ def preprocess(self, raw_handle: dict, *, l_freq: float, h_freq: float) -> dict:
         raw.save(path, overwrite=True, verbose="ERROR")
         return {"path": path, "filtered": [l_freq, h_freq]}
 
-    def epoch(self, filtered_handle: dict, *, tmin: float, tmax: float) -> dict:
+    def epoch(self, filtered_handle: dict, *, tmin: float, tmax: float,
+              event_id: Optional[int] = None) -> dict:
         import mne  # lazy
 
         raw = mne.io.read_raw_fif(filtered_handle["path"], preload=True,
                                   verbose="ERROR")
         # Read events FROM the data rather than assuming counts.
         events = mne.find_events(raw, verbose="ERROR")
-        epochs = mne.Epochs(raw, events, tmin=tmin, tmax=tmax,
-                            baseline=(None, 0), preload=True, verbose="ERROR")
+        # Filter to a SINGLE trigger code so we average one condition, not every
+        # trigger (standards + deviants + button presses). event_id semantics:
+        #   None   -> auto: the most-frequent code (the standard tone) [default]
+        #   "all"  -> no filter: epoch every trigger (legacy/diagnostic only)
+        #   int    -> that explicit trigger code
+        if event_id == "all":
+            mne_event_id = None  # mne: None means "use all event codes"
+        elif event_id is None:
+            # Most-frequent trigger code = the standard tone. Use stdlib Counter
+            # on plain ints so this needs no extra numpy surface.
+            from collections import Counter
+            codes = [int(c) for c in events[:, 2]]
+            mne_event_id = Counter(codes).most_common(1)[0][0] if codes else None
+            event_id = mne_event_id
+        else:
+            mne_event_id = int(event_id)
+            event_id = mne_event_id
+        epochs = mne.Epochs(raw, events, event_id=mne_event_id, tmin=tmin,
+                            tmax=tmax, baseline=(None, 0), preload=True,
+                            verbose="ERROR")
         path = self._path("epochs-epo.fif")
         epochs.save(path, overwrite=True, verbose="ERROR")
-        return {"path": path, "n_epochs": len(epochs), "window": [tmin, tmax]}
+        return {"path": path, "n_epochs": len(epochs), "window": [tmin, tmax],
+                "event_id": event_id}
 
     def evoked(self, epochs_handle: dict) -> dict:
         import mne  # lazy
@@ -193,11 +250,12 @@ def measure_peak(self, evoked_handle: dict, *,
         ev = mne.read_evokeds(evoked_handle["path"], verbose="ERROR")[0]
         if n_epochs <= 0 or len(ev.times) == 0:
             return {"latency_ms": None, "amplitude": None, "n_epochs": n_epochs}
-        tmin_s, tmax_s = window_ms[0] / 1000.0, window_ms[1] / 1000.0
-        # Search a slightly widened window so an out-of-window peak is still
-        # measurable (and can be reported as refuted).
-        lo = max(ev.times[0], tmin_s - 0.1)
-        hi = min(ev.times[-1], tmax_s + 0.1)
+        # Search WITHIN the expected window (+/- a small margin), clamped to the
+        # available sample times. A wide search reaches the large early stimulus
+        # artifact (~15 ms) and misreports it as the auditory peak; the pure
+        # helper below keeps the search in-window (see _peak_search_bounds).
+        lo, hi = _peak_search_bounds(window_ms, float(ev.times[0]),
+                                     float(ev.times[-1]))
         ch, latency_s, amp = ev.get_peak(tmin=lo, tmax=hi, mode="abs",
                                          return_amplitude=True)
         return {"latency_ms": float(latency_s) * 1000.0,
@@ -255,12 +313,15 @@ def fn(ctx: StepContext) -> StepResult:
         epochs = io.epoch(
             {"path": prior.get("path")},
             tmin=bounds.get("tmin", DEFAULT_BOUNDS["tmin"]),
-            tmax=bounds.get("tmax", DEFAULT_BOUNDS["tmax"]))
+            tmax=bounds.get("tmax", DEFAULT_BOUNDS["tmax"]),
+            event_id=bounds.get("event_id", DEFAULT_BOUNDS["event_id"]))
         return StepResult(
             data={"n_epochs": epochs.get("n_epochs"),
-                  "window": epochs.get("window"), "path": epochs.get("path")},
+                  "window": epochs.get("window"),
+                  "event_id": epochs.get("event_id"), "path": epochs.get("path")},
             artifacts=[{"path": epochs.get("path", ""), "kind": "epochs",
-                        "meta": {"n_epochs": epochs.get("n_epochs")}}],
+                        "meta": {"n_epochs": epochs.get("n_epochs"),
+                                 "event_id": epochs.get("event_id")}}],
         )
     return Step(name="epoch", fn=fn)
 

matilde_plugin/tools.py

@@ -17,6 +17,7 @@
 from __future__ import annotations
 
 import json
+import math
 from typing import Any
 
 # The engine lives INSIDE this plugin package (matilde_plugin/engine/), so it is
@@ -29,17 +30,46 @@
 # Shared envelope helpers
 # ---------------------------------------------------------------------------
 
+def _json_safe(obj: Any) -> Any:
+    """Recursively make *obj* strict-JSON-serializable.
+
+    Non-finite floats (``NaN``/``Infinity``) are the trap: ``json.dumps`` emits
+    them as bare ``NaN``/``Infinity`` literals, which are INVALID JSON, so a
+    strict downstream parser rejects the whole envelope and treats the tool call
+    as an error (this is what made ``matilde_study_status`` fail on a real run
+    whose finding evidence held a non-finite measured value). We convert those
+    to ``None`` so the envelope is always valid JSON. Numpy scalars are handled
+    by ``default=str`` at dump time.
+    """
+    if isinstance(obj, float):
+        return obj if math.isfinite(obj) else None
+    if isinstance(obj, dict):
+        return {k: _json_safe(v) for k, v in obj.items()}
+    if isinstance(obj, (list, tuple)):
+        return [_json_safe(v) for v in obj]
+    return obj
+
+
+def _dumps(payload: dict) -> str:
+    """Serialize an envelope to strictly-valid JSON (no bare NaN/Infinity)."""
+    # allow_nan=False makes any non-finite that slipped past _json_safe (e.g.
+    # inside a numpy type stringified by default=str — harmless) a hard failure
+    # rather than silent invalid output; _json_safe already neutralizes plain
+    # Python floats so this should not trigger.
+    return json.dumps(_json_safe(payload), default=str, allow_nan=False)
+
+
 def _tool_result(data: Any = None, **kwargs: Any) -> str:
     if data is not None:
         payload = data if isinstance(data, dict) else {"result": data}
     else:
         payload = kwargs
     payload.setdefault("success", True)
-    return json.dumps(payload, default=str)
+    return _dumps(payload)
 
 
 def _tool_error(message: str, **extra: Any) -> str:
-    return json.dumps({"error": message, "success": False, **extra}, default=str)
+    return _dumps({"error": message, "success": False, **extra})
 
 
 # ---------------------------------------------------------------------------

tests/test_meg_study.py

@@ -27,6 +27,7 @@
 
 from matilde_plugin.engine.meg_study import (  # noqa: E402
     MegIO,
+    _peak_search_bounds,
     build_steps,
 )
 from matilde_plugin.engine.pipeline import resume, run  # noqa: E402
@@ -50,6 +51,8 @@ def __init__(self, *, peak_latency_ms=100.0, peak_amp=5.0, n_epochs=40,
         self.fail_evoked = fail_evoked
         self.calls = {"fetch": 0, "preprocess": 0, "epoch": 0, "evoked": 0,
                       "measure_peak": 0}
+        # Records the event_id the study asked epoch to filter to (None=all).
+        self.epoch_event_id = "__unset__"
 
     def fetch_sample(self, *, dataset_id, bounds):
         self.calls["fetch"] += 1
@@ -60,8 +63,9 @@ def preprocess(self, raw, *, l_freq, h_freq):
         self.calls["preprocess"] += 1
         return {**raw, "filtered": [l_freq, h_freq], "path": "/tmp/fake_filt.fif"}
 
-    def epoch(self, filtered, *, tmin, tmax):
+    def epoch(self, filtered, *, tmin, tmax, event_id=None):
         self.calls["epoch"] += 1
+        self.epoch_event_id = event_id
         return {**filtered, "n_epochs": self.n_epochs,
                 "window": [tmin, tmax], "path": "/tmp/fake_epo.fif"}
 
@@ -97,6 +101,63 @@ def test_fakeio_satisfies_contract():
     assert isinstance(io, MegIO)
 
 
+# ---------------------------------------------------------------------------
+# Regression: the peak-search bounds must stay WITHIN the expected window.
+#
+# The original bug widened the search by +/-100 ms, so an 80-120 ms M100 window
+# became a ~-20..220 ms search that grabbed the large ~15 ms stimulus artifact
+# instead of the real ~100 ms auditory peak (reported 15 ms -> "refuted"). This
+# pure, mne-free helper is the unit that would have caught it.
+# ---------------------------------------------------------------------------
+
+def test_peak_search_bounds_stay_in_window_not_widened_by_100ms():
+    # Wide available times (e.g. -0.1 .. 0.3 s) must NOT let the search escape
+    # the expected window. For an 80-120 ms window we expect ~0.08..0.12 s,
+    # never the -0.02..0.22 s that the +/-100 ms bug produced.
+    lo, hi = _peak_search_bounds((80.0, 120.0), times_lo=-0.1, times_hi=0.3)
+    # A small symmetric margin is allowed, but nothing near +/-100 ms.
+    assert lo >= 0.06 and lo <= 0.08, lo      # not -0.02
+    assert hi >= 0.12 and hi <= 0.14, hi      # not 0.22
+    # And explicitly: the early stimulus artifact at ~15 ms is OUTSIDE the search.
+    assert lo > 0.015
+
+
+def test_peak_search_bounds_clamped_to_available_times():
+    # If the recording is shorter than the window, clamp to available samples.
+    lo, hi = _peak_search_bounds((80.0, 120.0), times_lo=0.09, times_hi=0.11)
+    assert lo == 0.09
+    assert hi == 0.11
+
+
+# ---------------------------------------------------------------------------
+# Regression: epoch must filter to the standard-tone event by default, not
+# average ALL triggers (standards + deviants + button presses), which muddied
+# the evoked average in the live run.
+# ---------------------------------------------------------------------------
+
+def test_epoch_applies_standards_event_filter_by_default(store):
+    sid = store.create_study(slug="evfilter", title="EvFilter", plan=_plan())
+    io = FakeMegIO()
+    run(store, sid, build_steps(dataset_id="bst_auditory", io=io))
+    # The study must NOT ask epoch to average every trigger. The default is the
+    # standards-only auto filter (event_id=None -> most-frequent code, resolved
+    # inside the real backend), never the "all events" sentinel that produced
+    # the muddied average behind the original mis-measurement.
+    assert io.epoch_event_id != "__unset__", "epoch was never called"
+    assert io.epoch_event_id != "all", (
+        "epoch was told to average ALL triggers -> muddied evoked average")
+    assert io.epoch_event_id is None, (
+        "default should be the auto standards filter (None), resolved in backend")
+
+
+def test_epoch_event_id_overridable_via_bounds(store):
+    sid = store.create_study(slug="evfilter2", title="EvFilter2", plan=_plan())
+    io = FakeMegIO()
+    run(store, sid, build_steps(dataset_id="bst_auditory", io=io,
+                                bounds={"event_id": 7}))
+    assert io.epoch_event_id == 7
+
+
 # ---------------------------------------------------------------------------
 # Happy path — peak at 100 ms is within the default 80-120 window -> supported.
 # ---------------------------------------------------------------------------

tests/test_meg_study_tools.py

@@ -39,9 +39,9 @@ def fetch_sample(self, *, dataset_id, bounds):
     def preprocess(self, raw, *, l_freq, h_freq):
         return {"path": "/tmp/fake_filt.fif", "filtered": [l_freq, h_freq]}
 
-    def epoch(self, filtered, *, tmin, tmax):
+    def epoch(self, filtered, *, tmin, tmax, event_id=None):
         return {"path": "/tmp/fake_epo.fif", "n_epochs": 30,
-                "window": [tmin, tmax]}
+                "window": [tmin, tmax], "event_id": event_id}
 
     def evoked(self, epochs):
         return {"path": "/tmp/fake_ave.fif", "n_epochs": 30}

tests/test_study_tools.py

@@ -82,6 +82,55 @@ def test_study_list_returns_created_studies(monkeypatch, tmp_path):
     assert "one" in slugs and "two" in slugs
 
 
+def _strict_loads(s: str):
+    """json.loads that REJECTS non-finite literals (NaN/Infinity), like the
+    strict JSON parsers a tool envelope is fed into downstream."""
+    def _reject(tok):
+        raise ValueError(f"non-finite JSON literal {tok!r}")
+    return json.loads(s, parse_constant=_reject)
+
+
+def test_study_status_returns_strict_json_for_completed_study(monkeypatch, tmp_path):
+    """Regression: a completed study's status must be a structured, strictly
+    JSON-parseable result — not an error envelope. A live measure_peak can yield
+    a non-finite float (NaN/Infinity) in a finding's evidence; json.dumps emits
+    those as bare NaN/Infinity, which is invalid JSON and made the framework
+    treat study_status as an error. The envelope must serialize them safely."""
+    plugin = _load_plugin()
+    monkeypatch.setenv("MATILDE_STUDY_DB", str(tmp_path / "s.db"))
+    from matilde_plugin.engine.store import StudyStore
+
+    created = json.loads(plugin._handle_study_create(
+        {"slug": "meg-nan", "title": "MEG NaN", "plan": ["validate_finding"]}))
+    sid = created["study_id"]
+
+    # Simulate what a real (lazy-mne) run can record: a finding whose evidence
+    # holds a non-finite float (e.g. a degenerate amplitude/latency).
+    store = StudyStore(str(tmp_path / "s.db"))
+    store.set_step_status(sid, "validate_finding", "done")
+    store.add_finding(
+        sid, "validate_finding",
+        claim="auditory M100 peak falls within 80-120 ms",
+        verdict="supported",
+        evidence={"latency_ms": 120.0, "amplitude": float("nan"),
+                  "extra": float("inf")})
+    store.set_study_status(sid, "done")
+    store.close()
+
+    raw = plugin._handle_study_status({"study_id": sid})
+    # Must parse under a STRICT parser (no bare NaN/Infinity tokens).
+    out = _strict_loads(raw)
+    assert out["success"] is True, out
+    assert "error" not in out
+    assert out["status"] == "done"
+    assert out["finding_count"] == 1
+    # The non-finite values are represented safely (None), not as bare NaN.
+    ev = out["findings"][0]["evidence"]
+    assert ev["amplitude"] is None
+    assert ev["extra"] is None
+    assert ev["latency_ms"] == 120.0
+
+
 def test_study_run_advances_and_is_resumable(monkeypatch, tmp_path):
     """End-to-end through the tools using the bibliography study. Coerces a string
     id, runs to completion, and a second run is a no-op (resume of a done study)."""