diff --git a/matilde_plugin/engine/meg_study.py b/matilde_plugin/engine/meg_study.py index cf1fd50..8c98250 100644 --- a/matilde_plugin/engine/meg_study.py +++ b/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) diff --git a/matilde_plugin/tools.py b/matilde_plugin/tools.py index 7778412..9c5a898 100644 --- a/matilde_plugin/tools.py +++ b/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}) # --------------------------------------------------------------------------- diff --git a/tests/test_meg_study.py b/tests/test_meg_study.py index bd7fc58..13cb0bf 100644 --- a/tests/test_meg_study.py +++ b/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. # --------------------------------------------------------------------------- diff --git a/tests/test_meg_study_tools.py b/tests/test_meg_study_tools.py index 46dc1bd..62ff5f9 100644 --- a/tests/test_meg_study_tools.py +++ b/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} diff --git a/tests/test_study_tools.py b/tests/test_study_tools.py index 02cc136..ab9a6e1 100644 --- a/tests/test_study_tools.py +++ b/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)."""