diff --git a/docs/auto.rst b/docs/auto.rst index 520469c..e3f684f 100644 --- a/docs/auto.rst +++ b/docs/auto.rst @@ -39,17 +39,26 @@ Pipeline target that tests nonlinear only if the Gaussian estimator's parents predict it at least as well as the plug-in's (below). Every selection stage is tested at 0.001, with enough permutations for the split level - and curtailment (which never changes a decision). + and curtailment (which never changes a decision). That check runs a Gaussian + selection with a local permutation null per nonlinear target, which can take + minutes per target on long, finely binned series; with a ``time_budget``, checks + not started within half the budget are skipped and those targets select with the + plug-in only (a recorded downgrade). 3. **Pilot.** The selection is run on three sampled targets, with the seeds the full run will use, so their skeletons are reused. The pilot measures the time per target of each estimator, the time per edge of the layers (at the resample count - the run will use), and of interpretation. + the run will use), and of interpretation. If no pilot target uses the Gaussian + estimator, its time is taken from the check above. 4. **Fit the budget.** The cost of the whole run is predicted from the pilot (targets, then edges, run ``threads`` at a time). If it exceeds ``time_budget``, options are dropped in a fixed order, KSG selection, layer estimates for shared-only candidates, hyperedges, then the source pre-screen is turned on, until it fits; the significance level is never loosened. With - ``strict_budget=True`` an unattainable budget raises instead of warning. + ``strict_budget=True`` an unattainable budget raises instead of warning. Because + selection times vary by orders of magnitude between targets (one Gaussian selection + with a local null took 47 s and another nearly an hour on the same iEEG data), the + prediction is also enforced: targets not started within 80% of the budget select + with the plug-in only, and the report lists them. 5. **Run** :func:`~infoflow.infer_multiplex` and return an :class:`~infoflow.auto.AutoResult` with the network and an :class:`~infoflow.auto.AutoReport`: diagnostics, gated targets, the final diff --git a/docs/preprocessing.rst b/docs/preprocessing.rst index b9d354d..3190da4 100644 --- a/docs/preprocessing.rst +++ b/docs/preprocessing.rst @@ -50,7 +50,9 @@ resolutions for parent selection instead (:func:`~infoflow.preprocess.scoring.network_scores`: each target predicted from its own past and its two most correlated sources, every variable in the same equal-frequency bins) and gives every node the best one, provided every node keeps -at least ``network_retention`` (default 0.7) of its own per-node score. +at least ``network_retention`` (default 0.7) of its own per-node score. The uniform +history sets the resolution of the contexts; each node keeps its own lag budget +(at least that history), so sources acting at longer delays remain candidates. The guard protects nodes whose own dynamics need fine symbols. On the MuTE network the oscillating driver keeps only 30-60% of its score at the network resolution; diff --git a/docs/references.bib b/docs/references.bib index 6b65bbd..dbe5788 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -31,6 +31,17 @@ @article{Benjamini2001 doi = {10.1214/aos/1013699998} } +@article{Besag1991, + author = {Besag, Julian and Clifford, Peter}, + title = {{Sequential Monte Carlo p-values}}, + journal = {Biometrika}, + volume = {78}, + number = {2}, + pages = {301--304}, + year = {1991}, + doi = {10.1093/biomet/78.2.301} +} + @article{Bertschinger2014, author = {Bertschinger, Nils and Rauh, Johannes and Olbrich, Eckehard and Jost, Jürgen and Ay, Nihat}, title = {Quantifying Unique Information}, diff --git a/docs/workflow.rst b/docs/workflow.rst index 400405e..f0f292b 100644 --- a/docs/workflow.rst +++ b/docs/workflow.rst @@ -10,6 +10,24 @@ conditions on the union of their edges, for trials within a subject (``dependent=True``) or independent units. Condition labels are permuted (swapped within pairs when dependent) and p-values are corrected per layer. +The permutation test is exact only if the tested edges do not depend on the labels. +Edges selected separately in each condition are partly selected *because* of +chance differences between the conditions, which the permutations then destroy: on +random halves of one condition of an iEEG dataset, 13-15% of splits gave a +significant intrinsic or synergistic difference instead of 5%. For a calibrated +test, select one network on both conditions pooled and pass it as both networks: + +.. code-block:: python + + pooled = infoflow.infer_multiplex(data.subset_trials(idx_a + idx_b), prng=0) + cmp = compare_networks(pooled, pooled, data.subset_trials(idx_a), data.subset_trials(idx_b)) + +Every permutation re-estimates every tested edge, which is slow for finely binned +contexts. ``curtail=10`` stops an edge's permutations once each layer's permuted +difference has reached the observed one ten times (sequential Monte Carlo p-values, +:cite:`Besag1991`), which leaves small p-values unchanged and stops edges with large +ones early; ``map_fn`` runs edges in parallel. + Parallelism ----------- diff --git a/examples/catfr/cohort.py b/examples/catfr/cohort.py new file mode 100644 index 0000000..630838c --- /dev/null +++ b/examples/catfr/cohort.py @@ -0,0 +1,197 @@ +""" +Cohort follow-up: region-level subsequent-memory contrast across subjects. + +1. Choose one session per subject with at least two non-overlapping bipolar pairs in + each of MTL, LTC and PFC and at least ``--min-trials`` words in the smaller of the + recalled and forgotten classes, ranked by that count. +2. For each session, download the bipolar EDF, write epochs with ``prepare.py``'s + ``prepare`` (up to ``--per-roi`` pairs per region), and delete the EDF again unless + it was already on disk. +3. Average each region's HFA (each pair z-scored over the session's epochs first) into + one node, so every subject has the same three nodes. +4. Discretize every subject with the same 4 equal-frequency bins and test, across + subjects, recalled minus forgotten flow on all six directed region pairs with + ``compare_networks(design="between", dependent=True)``. The tested variables are + fixed in advance (source at ``--source-lag``, target at lags 1-2, the third region + at ``--source-lag``), so nothing is selected on the data. +5. Within each subject, the pilot's analysis on the three nodes: one network selected + on both conditions pooled, and a within-subject comparison on its edges. + + uv run --project ../.. --with mne --with netcdf4 python cohort.py --subjects 25 +""" + +import argparse +import json +import subprocess +from collections import Counter +from types import SimpleNamespace + +import numpy as np +import xarray as xr +from common import BUCKET, OUT, ROIS, ROOT, pair_regions, sampling_rate, session_files, word_events +from prepare import pick_pairs, prepare +from run import condition_network, layer_counts, significant_edges, tested + +from infoflow.compare import compare_networks +from infoflow.layers import LAYERS +from infoflow.parallel import thread_map +from infoflow.preprocess import EqualFrequency, preprocess +from infoflow.selection import TargetSkeleton + +COHORT = OUT / "cohort" + + +def choose(min_trials, per_roi): + rows = [] + for events in sorted(ROOT.glob("sub-*/ses-*/ieeg/*_task-catFR1_events.tsv")): + subject, session = events.parts[-4][4:], events.parts[-3][4:] + files = session_files(subject, session) + if files["electrodes"] is None or not files["channels"].exists(): + continue + words = word_events(events, sampling_rate(files)) + n_rec = sum(w["recalled"] for w in words) + regions = pair_regions(files["channels"], files["electrodes"]) + counts = Counter(regions[p][0] for p in pick_pairs(regions, per_roi)) + balanced = min(n_rec, len(words) - n_rec) + if balanced >= min_trials and all(counts[r] >= 2 for r in ROIS): + rows.append({"subject": subject, "session": session, "balanced": balanced, **{r: counts[r] for r in ROIS}}) + best = {} + for r in sorted(rows, key=lambda r: -r["balanced"]): + best.setdefault(r["subject"], r) + return sorted(best.values(), key=lambda r: -r["balanced"]) + + +def fetch_and_prepare(row, per_roi): + stem = f"sub-{row['subject']}_ses-{row['session']}" + if (OUT / f"{stem}_hfa.nc").exists(): + return + files = session_files(row["subject"], row["session"]) + downloaded = not files["edf"].exists() + if downloaded: + subprocess.run( + ["aws", "s3", "cp", "--no-sign-request", "--only-show-errors", BUCKET + files["key"], str(files["edf"])], + check=True, + ) + try: + prepare(row["subject"], row["session"], per_roi) + finally: + if downloaded: + files["edf"].unlink(missing_ok=True) + + +def region_trials(path): + """Trials of shape (time, 3): each region's z-scored HFA averaged over its pairs.""" + ds = xr.open_dataset(path).load() + x = ds.signal.values + x = (x - x.mean(axis=(0, 1))) / x.std(axis=(0, 1)) + region = ds.region.values + nodes = np.stack([x[:, :, region == r].mean(axis=2) for r in ROIS], axis=2) + return [nodes[i] for i in range(len(nodes))], ds.recalled.values.astype(bool) + + +def fixed_network(source_lag): + skeleton = {} + for t in range(len(ROIS)): + sources = [(s, source_lag) for s in range(len(ROIS)) if s != t] + skeleton[t] = TargetSkeleton(target=t, target_past=[(t, 1), (t, 2)], sources=sources, significant=True) + return SimpleNamespace(names=list(ROIS), skeleton=skeleton) + + +def main(): + parser = argparse.ArgumentParser() + parser.add_argument("--subjects", type=int, default=25) + parser.add_argument("--min-trials", type=int, default=40) + parser.add_argument("--per-roi", type=int, default=10) + parser.add_argument("--source-lag", type=int, default=2) + parser.add_argument("--n-perm", type=int, default=5000) + parser.add_argument("--threads", type=int, default=4) + args = parser.parse_args() + COHORT.mkdir(parents=True, exist_ok=True) + + rows = choose(args.min_trials, args.per_roi) + print(f"{len(rows)} eligible subjects", flush=True) + data_rec, data_forg, per_subject, skipped = [], [], [], {} + for row in rows: + if len(per_subject) == args.subjects: + break + stem = f"sub-{row['subject']}_ses-{row['session']}" + try: + fetch_and_prepare(row, args.per_roi) + except ValueError as error: + # e.g. no WORD onset falls inside the recording + skipped[stem] = str(error) + print(f"{stem}: skipped ({error})", flush=True) + continue + trials, recalled = region_trials(OUT / f"{stem}_hfa.nc") + rng = np.random.default_rng(0) + rec, forg = np.flatnonzero(recalled), np.flatnonzero(~recalled) + n = min(len(rec), len(forg)) + idx_rec = np.sort(rng.choice(rec, n, replace=False)) + idx_forg = np.sort(rng.choice(forg, n, replace=False)) + pre = preprocess(trials, names=list(ROIS), max_lag=15, discretizer=EqualFrequency(4), prng=0) + data_rec.append(pre.data.subset_trials(list(idx_rec))) + data_forg.append(pre.data.subset_trials(list(idx_forg))) + + net = condition_network(pre, np.concatenate([idx_rec, idx_forg]), list(ROIS), args.threads) + within = compare_networks( + net, + net, + data_rec[-1], + data_forg[-1], + design="within", + n_perm=500, + curtail=10, + map_fn=thread_map(args.threads), + prng=0, + ) + per_subject.append( + { + **row, + "n_per_condition": int(n), + "edges": layer_counts(net), + "tested": tested(within), + "sme_edges": significant_edges(within, list(ROIS)), + } + ) + print( + f"{stem}: {per_subject[-1]['n_per_condition']} trials/condition, edges {per_subject[-1]['edges']}, " + f"SME {per_subject[-1]['sme_edges']}", + flush=True, + ) + + fixed = fixed_network(args.source_lag) + between = compare_networks( + fixed, fixed, data_rec, data_forg, design="between", dependent=True, n_perm=args.n_perm, prng=0 + ) + between.dataset.to_netcdf(COHORT / "between.nc") + from infoflow.compare import _edge_values + + per_unit = {} + for (s, t), (x_vars, w_vars) in between.edges.items(): + per_unit[f"{ROIS[s]}->{ROIS[t]}"] = { + "recalled": [[_edge_values(d, t, x_vars, w_vars)[k] for k in LAYERS] for d in data_rec], + "forgotten": [[_edge_values(d, t, x_vars, w_vars)[k] for k in LAYERS] for d in data_forg], + } + summary = { + "subjects": per_subject, + "skipped": skipped, + "source_lag": args.source_lag, + "n_perm": args.n_perm, + "between": { + f"{ROIS[s]}->{ROIS[t]}": { + k: { + v: float(between.dataset[v].sel(layer=k).values[s, t]) + for v in ("a", "b", "difference", "pvalue", "qvalue") + } + for k in LAYERS + } + for s, t in between.edges + }, + "per_unit": per_unit, + } + (COHORT / "summary.json").write_text(json.dumps(summary, indent=2, default=str)) + print(json.dumps(summary["between"], indent=1), flush=True) + + +if __name__ == "__main__": + main() diff --git a/examples/catfr/common.py b/examples/catfr/common.py new file mode 100644 index 0000000..c0ee3de --- /dev/null +++ b/examples/catfr/common.py @@ -0,0 +1,128 @@ +""" +Shared readers for the catFR1 pilot (OpenNeuro ds004809, Herrema & Kahana). + +The BIDS events carry no ``recalled`` flag: a WORD is recalled when a REC_WORD +in the same list names the same item. Regions come from the hippocampal-subfield +atlases first (``das.region``, ``stein.region``), then the whole-brain +(``wb.region``) and Desikan-Killiany (``ind.region``) labels. +""" + +import csv +import json +import os +import re +from pathlib import Path + +DATASET = "ds004809" +BUCKET = f"s3://openneuro.org/{DATASET}/" +ROOT = Path(os.environ.get("CATFR_ROOT", Path.home() / "data" / "openneuro" / DATASET)) +OUT = Path(os.environ.get("CATFR_OUT", Path.home() / "data" / "catfr_pilot")) + +ROIS = ("MTL", "LTC", "PFC") + +_MTL = re.compile( + r"hippocamp|amygdala|\bamy\b|parahippocampal|\bphg\b|\bphc\b|entorhinal|\berc\b|\bec\b|" + r"\bprc\b|\bba3[56]\b|\bca[123]\b|\bdg\b|\bsub\b", + re.I, +) +_LTC = re.compile( + r"superiortemporal|middletemporal|inferiortemporal|bankssts|temporalpole|" + r"temporal gyrus|\b[smi]tg\b|\btmp\b", + re.I, +) +_PFC = re.compile( + r"middlefrontal|superiorfrontal|pars(triangularis|opercularis|orbitalis)|orbitofrontal|" + r"frontalpole|frontal gyrus|\b[sm]fg\b|\bifg\b|dlpfc|\bofc\b", + re.I, +) +_SKIP = re.compile(r"white matter|\bwm\b|^n/a$|^$|misc|sulcus", re.I) + + +def read_tsv(path): + with open(path, newline="") as f: + return list(csv.DictReader(f, delimiter="\t")) + + +def roi(labels): + """Map an ordered list of atlas labels (most specific first) to MTL, LTC, PFC or None.""" + for label in labels: + label = (label or "").strip() + if _SKIP.search(label): + continue + for name, pattern in (("MTL", _MTL), ("LTC", _LTC), ("PFC", _PFC)): + if pattern.search(label): + return name + return None + return None + + +def contact_regions(electrodes_tsv): + """``contact -> (roi, label)`` from an electrodes table.""" + out = {} + for row in read_tsv(electrodes_tsv): + labels = [row.get(k) for k in ("das.region", "stein.region", "wb.region", "ind.region")] + named = next((l for l in labels if l and not _SKIP.search(l)), "n/a") + out[row["name"]] = (roi(labels), named) + return out + + +def pair_regions(channels_tsv, electrodes_tsv): + """``bipolar pair -> (roi, label)`` for pairs whose contacts agree on one ROI.""" + contacts = contact_regions(electrodes_tsv) + out = {} + for row in read_tsv(channels_tsv): + a, _, b = row["name"].partition("-") + ra, la = contacts.get(a, (None, "n/a")) + rb, _ = contacts.get(b, (None, "n/a")) + if row.get("status", "good") == "bad": + continue + if ra is not None and ra == rb: + out[row["name"]] = (ra, la) + return out + + +def sampling_rate(files): + return float(json.loads(files["sidecar"].read_text())["SamplingFrequency"]) + + +def word_events(events_tsv, sfreq): + """ + Non-practice WORD events as dicts with ``onset``, ``list``, ``item``, ``category``, + ``serialpos``, ``recalled``. + + ``onset`` is ``sample / sfreq`` seconds from the start of the recording. The BIDS + ``onset`` column of this dataset is not: it is offset from the recording by 1 s to + 15 min (drifting within a session), and only ``sample``-aligned epochs show the + word-locked evoked response. + """ + rows = read_tsv(events_tsv) + recalled = {(r["list"], r["item_name"].upper()) for r in rows if r["trial_type"] == "REC_WORD"} + words = [] + for r in rows: + if r["trial_type"] != "WORD" or int(r["list"]) < 1 or int(r["sample"]) < 0: + continue + words.append( + { + "onset": int(r["sample"]) / sfreq, + "list": int(r["list"]), + "item": r["item_name"], + "category": r["category"], + "serialpos": int(r["serialpos"]), + "recalled": (r["list"], r["item_name"].upper()) in recalled, + } + ) + return words + + +def session_files(subject, session, root=ROOT): + base = root / f"sub-{subject}" / f"ses-{session}" / "ieeg" + stem = f"sub-{subject}_ses-{session}_task-catFR1" + electrodes = sorted(base.glob(f"{stem}_space-*_electrodes.tsv")) + return { + "events": base / f"{stem}_events.tsv", + "channels": base / f"{stem}_acq-bipolar_channels.tsv", + "electrodes": electrodes[0] if electrodes else None, + "edf": base / f"{stem}_acq-bipolar_ieeg.edf", + "sidecar": base / f"{stem}_acq-bipolar_ieeg.json", + "key": f"sub-{subject}/ses-{session}/ieeg/{stem}_acq-bipolar_ieeg.edf", + } diff --git a/examples/catfr/fetch.py b/examples/catfr/fetch.py new file mode 100644 index 0000000..06c44c6 --- /dev/null +++ b/examples/catfr/fetch.py @@ -0,0 +1,106 @@ +""" +Rank catFR1 sessions for the pilot and download their bipolar EDFs. + +1. Sync only the small sidecars (events, bipolar channels, electrodes) from the + public OpenNeuro bucket. +2. Score each session by bipolar-pair coverage of MTL, lateral temporal and + prefrontal cortex, and by the number of recalled words. +3. Download the bipolar EDF of the best session of the top subjects. + + uv run --with mne python examples/catfr/fetch.py --top 3 +""" + +import argparse +import json +import subprocess +from collections import Counter + +from common import BUCKET, OUT, ROIS, ROOT, pair_regions, sampling_rate, session_files, word_events + + +def sync_sidecars(): + includes = [ + "participants.tsv", + "*_events.tsv", + "*acq-bipolar_channels.tsv", + "*_electrodes.tsv", + "*acq-bipolar_ieeg.json", + ] + cmd = ["aws", "s3", "sync", "--no-sign-request", "--only-show-errors", BUCKET, str(ROOT), "--exclude", "*"] + for pattern in includes: + cmd += ["--include", pattern] + subprocess.run(cmd, check=True) + + +def score_sessions(min_recalled=60, min_pairs=2): + rows = [] + for events in sorted(ROOT.glob("sub-*/ses-*/ieeg/*_task-catFR1_events.tsv")): + subject = events.parts[-4][4:] + session = events.parts[-3][4:] + files = session_files(subject, session) + if files["electrodes"] is None or not files["channels"].exists(): + continue + words = word_events(events, sampling_rate(files)) + n_recalled = sum(w["recalled"] for w in words) + counts = Counter(r for r, _ in pair_regions(files["channels"], files["electrodes"]).values()) + rows.append( + { + "subject": subject, + "session": session, + "n_words": len(words), + "n_recalled": n_recalled, + **{roi: counts.get(roi, 0) for roi in ROIS}, + } + ) + eligible = [r for r in rows if r["n_recalled"] >= min_recalled and all(r[k] >= min_pairs for k in ROIS)] + eligible.sort(key=lambda r: (min(min(r[k] for k in ROIS), 6), r["n_recalled"]), reverse=True) + return rows, eligible + + +def main(): + parser = argparse.ArgumentParser() + parser.add_argument("--top", type=int, default=3) + parser.add_argument("--skip-sync", action="store_true") + parser.add_argument("--no-download", action="store_true") + args = parser.parse_args() + + if not args.skip_sync: + sync_sidecars() + rows, eligible = score_sessions() + chosen, seen = [], set() + for r in eligible: + if r["subject"] not in seen: + chosen.append(r) + seen.add(r["subject"]) + if len(chosen) == args.top: + break + + OUT.mkdir(parents=True, exist_ok=True) + (OUT / "sessions.json").write_text( + json.dumps({"n_sessions": len(rows), "n_eligible": len(eligible), "chosen": chosen}, indent=2) + ) + print(f"{len(rows)} sessions, {len(eligible)} eligible") + for r in chosen: + print(r) + + if args.no_download: + return + for r in chosen: + files = session_files(r["subject"], r["session"]) + if not files["edf"].exists(): + subprocess.run( + [ + "aws", + "s3", + "cp", + "--no-sign-request", + "--only-show-errors", + BUCKET + files["key"], + str(files["edf"]), + ], + check=True, + ) + + +if __name__ == "__main__": + main() diff --git a/examples/catfr/null_check.py b/examples/catfr/null_check.py new file mode 100644 index 0000000..84f8f96 --- /dev/null +++ b/examples/catfr/null_check.py @@ -0,0 +1,74 @@ +""" +Per-edge p-values of same-condition null splits, to check the calibration of +``compare_networks`` under pooled selection. + +Each split draws two disjoint halves of the larger class exactly as ``run.py`` does +(seeds ``1000 + r``), selects one network on both halves together, and records every +tested edge's raw p-value per layer along with the FDR decisions. Under the null the +raw p-values should be (super-)uniform; a per-split excess with uniform p-values points +at the multiple-testing step instead. + + uv run --project ../.. --with netcdf4 python null_check.py --subject R1332M --start 0 --count 25 +""" + +import argparse +import json + +import numpy as np +import xarray as xr +from common import OUT +from run import DEFAULT_MAX_LAG, contrast + +from infoflow.layers import LAYERS +from infoflow.parallel import thread_map +from infoflow.preprocess import preprocess + + +def main(): + parser = argparse.ArgumentParser() + parser.add_argument("--subject", required=True) + parser.add_argument("--signal", default="hfa") + parser.add_argument("--start", type=int, default=0) + parser.add_argument("--count", type=int, default=25) + parser.add_argument("--n-perm", type=int, default=500) + parser.add_argument("--threads", type=int, default=4) + args = parser.parse_args() + + path = next(OUT.glob(f"sub-{args.subject}_*_{args.signal}.nc")) + ds = xr.open_dataset(path).load() + regions = [str(r) for r in ds.region.values] + names = [f"{r}:{c}" for r, c in zip(regions, ds.channel.values, strict=True)] + trials = [ds.signal.values[i].astype(float) for i in range(ds.sizes["trial"])] + recalled = ds.recalled.values.astype(bool) + rec, forg = np.flatnonzero(recalled), np.flatnonzero(~recalled) + major = rec if len(rec) >= len(forg) else forg + m = min(len(rec), len(forg), len(major) // 2) + pre = preprocess( + trials, + names=names, + max_lag=DEFAULT_MAX_LAG[args.signal], + network_resolution="auto", + map_fn=thread_map(args.threads), + prng=0, + ) + + out = OUT / f"nullcheck_{path.stem}_{args.start:03d}.jsonl" + with out.open("w") as f: + for r in range(args.start, args.start + args.count): + draw = np.random.default_rng(1000 + r).permutation(major) + _, _, cmp = contrast( + pre, np.sort(draw[:m]), np.sort(draw[m : 2 * m]), names, args.n_perm, args.threads, 1000 + r, "pooled" + ) + record = {"split": r, "pvalues": {}, "significant": {}} + for k in LAYERS: + p = cmp.dataset["pvalue"].sel(layer=k).values + s, t = np.nonzero(np.isfinite(p)) + record["pvalues"][k] = [[names[i], names[j], float(p[i, j])] for i, j in zip(s, t, strict=True)] + record["significant"][k] = int(cmp.dataset["significant"].sel(layer=k).sum()) + f.write(json.dumps(record) + "\n") + f.flush() + print(f"split {r}: {record['significant']}", flush=True) + + +if __name__ == "__main__": + main() diff --git a/examples/catfr/prepare.py b/examples/catfr/prepare.py new file mode 100644 index 0000000..5394f6e --- /dev/null +++ b/examples/catfr/prepare.py @@ -0,0 +1,133 @@ +""" +Turn a catFR1 bipolar session into encoding epochs for infoflow. + +For each chosen session (``sessions.json`` from ``fetch.py``) this picks up to +``--per-roi`` bipolar pairs in each of MTL, LTC and PFC (spread over atlas +labels), notch-filters line noise, and writes two signals of WORD epochs +(0 to 1.6 s after onset): + +* ``broadband``: 1-80 Hz, resampled to 200 Hz; +* ``hfa``: log mean Morlet power over 8 log-spaced frequencies in 70-150 Hz, + computed on epochs padded by 1 s each side and then cropped, averaged down to + 50 Hz. + +Every non-practice word is kept with its ``recalled`` label; ``run.py`` balances. + + uv run --with mne python examples/catfr/prepare.py +""" + +import argparse +import json +from collections import defaultdict +from itertools import zip_longest + +import mne +import numpy as np +import xarray as xr +from common import OUT, ROIS, pair_regions, session_files, word_events + +EPOCH = 1.6 +PAD = 1.0 +BROADBAND_SFREQ = 200.0 +HFA_SFREQ = 50.0 +HFA_FREQS = np.geomspace(70, 150, 8) + + +def pick_pairs(regions, per_roi): + """ + Round-robin over atlas labels within each ROI so one electrode shaft doesn't + dominate. Pairs sharing a contact share a signal term, so at most one is kept. + """ + picks, used = [], set() + for roi in ROIS: + by_label = defaultdict(list) + for name, (r, label) in sorted(regions.items()): + if r == roi: + by_label[label].append(name) + ordered = [n for group in zip_longest(*by_label.values()) for n in group if n is not None] + chosen = [] + for name in ordered: + contacts = set(name.split("-")) + if len(chosen) < per_roi and not contacts & used: + chosen.append(name) + used |= contacts + picks += chosen + return picks + + +def epochs_of(raw, onsets, tmin, tmax): + sfreq = raw.info["sfreq"] + starts = np.round((onsets + tmin) * sfreq).astype(int) + n = int(round((tmax - tmin) * sfreq)) + data = raw.get_data() + keep = (starts >= 0) & (starts + n <= data.shape[1]) + return np.stack([data[:, s : s + n] for s in starts[keep]]), keep + + +def to_dataset(values, sfreq, words, keep, picks, regions, attrs): + kept = [w for w, k in zip(words, keep) if k] + return xr.Dataset( + {"signal": (("trial", "time", "channel"), values.transpose(0, 2, 1).astype(np.float32))}, + coords={ + "time": np.arange(values.shape[2]) / sfreq, + "channel": picks, + "region": ("channel", [regions[p][0] for p in picks]), + "label": ("channel", [regions[p][1] for p in picks]), + "recalled": ("trial", [w["recalled"] for w in kept]), + "list": ("trial", [w["list"] for w in kept]), + "serialpos": ("trial", [w["serialpos"] for w in kept]), + "item": ("trial", [w["item"] for w in kept]), + }, + attrs={**attrs, "sfreq": sfreq}, + ) + + +def prepare(subject, session, per_roi): + files = session_files(subject, session) + regions = pair_regions(files["channels"], files["electrodes"]) + raw = mne.io.read_raw_edf(files["edf"], preload=False, verbose="error") + words = word_events(files["events"], raw.info["sfreq"]) + onsets = np.array([w["onset"] for w in words]) + + regions = {k: v for k, v in regions.items() if k in raw.ch_names} + picks = pick_pairs(regions, per_roi) + raw.pick(picks).load_data(verbose="error") + sfreq = raw.info["sfreq"] + raw.notch_filter(np.arange(60, sfreq / 2 - 1, 60), verbose="error") + attrs = {"subject": subject, "session": session, "source_sfreq": sfreq} + + bb = raw.copy().filter(1, 80, verbose="error").resample(BROADBAND_SFREQ, verbose="error") + values, keep = epochs_of(bb, onsets, 0.0, EPOCH) + broadband = to_dataset(values, BROADBAND_SFREQ, words, keep, picks, regions, {**attrs, "signal": "broadband"}) + + hf = raw.copy().resample(500.0, verbose="error") + padded, keep = epochs_of(hf, onsets, -PAD, EPOCH + PAD) + power = mne.time_frequency.tfr_array_morlet(padded, 500.0, HFA_FREQS, n_cycles=7, output="power", verbose="error") + lo = int(PAD * 500) + env = np.log10(power.mean(axis=2))[:, :, lo : lo + int(EPOCH * 500)] + step = int(500 / HFA_SFREQ) + env = env.reshape(*env.shape[:2], -1, step).mean(axis=-1) + hfa = to_dataset(env, HFA_SFREQ, words, keep, picks, regions, {**attrs, "signal": "hfa"}) + + stem = OUT / f"sub-{subject}_ses-{session}" + broadband.to_netcdf(f"{stem}_broadband.nc") + hfa.to_netcdf(f"{stem}_hfa.nc") + return broadband, hfa + + +def main(): + parser = argparse.ArgumentParser() + parser.add_argument("--per-roi", type=int, default=6) + args = parser.parse_args() + for s in json.loads((OUT / "sessions.json").read_text())["chosen"]: + bb, hfa = prepare(s["subject"], s["session"], args.per_roi) + rec = bb.recalled.values + print( + f"{s['subject']} ses-{s['session']}: {bb.sizes['channel']} pairs " + f"{dict(zip(*np.unique(bb.region.values, return_counts=True)))}, " + f"{rec.sum()} recalled / {(~rec).sum()} forgotten, broadband {bb.signal.shape}, hfa {hfa.signal.shape}" + ) + + +if __name__ == "__main__": + main() diff --git a/examples/catfr/report.py b/examples/catfr/report.py new file mode 100644 index 0000000..f41b1bd --- /dev/null +++ b/examples/catfr/report.py @@ -0,0 +1,228 @@ +""" +Tables and figures for ``reports/catfr_pilot.md`` from ``run.py``'s outputs. + + uv run --project ../.. --with netcdf4 --with matplotlib python report.py +""" + +import argparse +import json +from collections import Counter +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np +import xarray as xr +from common import OUT +from scipy.stats import binomtest + +from infoflow.layers import LAYERS + +FIGURES = Path(__file__).resolve().parents[2] / "reports" / "figures" +SIGNALS = ("hfa", "hfa-long", "hfa-pooled", "broadband", "broadband-long", "broadband-pooled") + + +def from_outputs(stem): + """A summary for a run stopped during its null splits, from its saved contrast and networks.""" + cmp = xr.open_dataset(OUT / f"{stem}_compare.nc") + names = [str(n) for n in cmp.source.values] + signal = "hfa" if "_hfa" in stem else "broadband" + epochs = xr.open_dataset(OUT / f"{stem.split(f'_{signal}')[0]}_{signal}.nc") + rec = epochs.recalled.values.astype(bool) + edges = {} + for label in ("recalled", "forgotten"): + net = xr.open_dataset(OUT / f"{stem}_net_{label}.nc") + edges[label] = {k: int(net["significant"].sel(layer=k).sum()) for k in LAYERS} + sme = [] + for k in LAYERS: + for s, t in zip(*np.nonzero(cmp["significant"].sel(layer=k).values)): + sme.append( + { + "layer": k, + "source": names[s], + "target": names[t], + "regions": f"{names[s].split(':')[0]}->{names[t].split(':')[0]}", + "recalled": float(cmp["a"].sel(layer=k).values[s, t]), + "forgotten": float(cmp["b"].sel(layer=k).values[s, t]), + "qvalue": float(cmp["qvalue"].sel(layer=k).values[s, t]), + } + ) + return { + "file": f"{stem.split(f'_{signal}')[0]}_{signal}.nc", + "signal": signal, + "n_per_condition": int(min(rec.sum(), (~rec).sum())), + "edges": edges, + "tested": {k: int(np.isfinite(cmp["pvalue"].sel(layer=k)).sum()) for k in LAYERS}, + "sme_edges": sme, + "null": {"reps": []}, + "timing": {}, + "lags": "long" if "_long" in stem else "auto", + "selection": "pooled" if "_pooled" in stem else "condition", + "partial": True, + } + + +def load(): + out = {} + stems = {p.name.removesuffix("_summary.json") for p in OUT.glob("sub-*_summary.json")} + partial = {p.name.removesuffix("_compare.nc") for p in OUT.glob("sub-*_compare.nc")} - stems + for stem in sorted(stems | partial): + s = from_outputs(stem) if stem in partial else json.loads((OUT / f"{stem}_summary.json").read_text()) + s["stem"] = stem + variant = s["signal"] + ("-long" if s.get("lags") == "long" else "") + variant += "-pooled" if s.get("selection") == "pooled" else "" + out[(s["file"].split("_ses")[0][4:], variant)] = s + return out + + +def runtime(seconds): + return f"{seconds / 60:.0f} min" if seconds < 5400 else f"{seconds / 3600:.1f} h" + + +def tables(summaries): + lines = [ + "| Subject | Signal | Trials / condition | Edges recalled (I/S/Sh) | Edges forgotten (I/S/Sh) " + "| Tested (I/S/Sh) | Significant SME (I/S/Sh) | Null splits with any significant (I/S/Sh) |", + "|---|---|---|---|---|---|---|---|", + ] + for (sub, sig), s in summaries.items(): + sme = Counter(e["layer"] for e in s["sme_edges"]) + reps = s["null"]["reps"] + any_sig = [sum(r["significant"][k] > 0 for r in reps) for k in LAYERS] + fmt = lambda d: "/".join(str(d.get(k, 0)) for k in LAYERS) # noqa: E731 + lines.append( + f"| {sub} | {sig} | {s['n_per_condition']} | {fmt(s['edges']['recalled'])} | " + f"{fmt(s['edges']['forgotten'])} | {fmt(s['tested'])} | {fmt(sme)} | " + f"{'/'.join(map(str, any_sig))} of {len(reps)} |" + ) + + lines += [ + "", + "| Subject | Signal | Layer | Edge | Regions | Recalled (bits) | Forgotten (bits) | q |", + "|---|---|---|---|---|---|---|---|", + ] + for (sub, sig), s in summaries.items(): + for e in s["sme_edges"]: + lines.append( + f"| {sub} | {sig} | {e['layer']} | {e['source']} -> {e['target']} | {e['regions']} | " + f"{e['recalled']:.4f} | {e['forgotten']:.4f} | {e['qvalue']:.3f} |" + ) + + lines += [ + "", + "| Subject | Signal | Null: significant / tested edges, pooled over splits (I/S/Sh) |", + "|---|---|---|", + ] + for (sub, sig), s in summaries.items(): + reps = s["null"]["reps"] + cells = [f"{sum(r['significant'][k] for r in reps)}/{sum(r['tested'][k] for r in reps)}" for k in LAYERS] + lines.append(f"| {sub} | {sig} | {' '.join(cells)} |") + + lines += [ + "", + "| Variant | Null splits | Any significant: intrinsic | synergistic | shared |", + "|---|---|---|---|---|", + ] + for variant in SIGNALS: + reps = [r for (_, sig), s in summaries.items() if sig == variant for r in s["null"]["reps"]] + if not reps: + continue + cells = [] + for k in LAYERS: + hits = sum(r["significant"][k] > 0 for r in reps) + p = binomtest(hits, len(reps), 0.05, alternative="greater").pvalue + cells.append(f"{hits} ({hits / len(reps):.0%}, p={p:.3f})") + lines.append(f"| {variant} | {len(reps)} | {' | '.join(cells)} |") + + lines += [ + "", + "| Subject | Signal | Preprocess | Contrast | Null (all splits) | auto.infer | auto edges (I/S/Sh) |", + "|---|---|---|---|---|---|---|", + ] + for (sub, sig), s in summaries.items(): + t = s["timing"] + if not t: + continue + auto = s.get("auto_edges", {}) + lines.append( + f"| {sub} | {sig} | {runtime(t['preprocess'])} | {runtime(t['contrast'])} | " + f"{runtime(t['null'])} | {runtime(t['auto']) if 'auto' in t else 'n/a'} | " + f"{'/'.join(str(auto.get(k, 0)) for k in LAYERS) if auto else 'n/a'} |" + ) + + lines += ["", "Broadband vs HFA overlap of significant SME edges (layer, source, target):", ""] + for suffix in ("", "-long", "-pooled"): + for sub in sorted({k[0] for k in summaries}): + keys = (f"hfa{suffix}", f"broadband{suffix}") + if all((sub, k) in summaries for k in keys): + hfa, bb = ( + {(e["layer"], e["source"], e["target"]) for e in summaries[(sub, k)]["sme_edges"]} for k in keys + ) + lines.append(f"- {sub}{suffix}: HFA {len(hfa)}, broadband {len(bb)}, shared {len(hfa & bb)}") + return "\n".join(lines) + + +def figure(summaries): + subjects = sorted({k[0] for k in summaries}) + for sig in SIGNALS: + rows = [s for s in subjects if (s, sig) in summaries] + if not rows: + continue + fig, axes = plt.subplots(len(rows), len(LAYERS), figsize=(4 * len(LAYERS), 3.8 * len(rows)), squeeze=False) + for i, sub in enumerate(rows): + ds = xr.open_dataset(OUT / f"{summaries[(sub, sig)]['stem']}_compare.nc") + names = list(ds.source.values) + for j, layer in enumerate(LAYERS): + ax = axes[i, j] + d = ds["difference"].sel(layer=layer).values + lim = float(np.max(np.nan_to_num(np.abs(ds["difference"].values)))) or 1.0 + ax.imshow(d, cmap="RdBu_r", vmin=-lim, vmax=lim) + for s, t in zip(*np.nonzero(ds["significant"].sel(layer=layer).values)): + ax.text(t, s, "*", ha="center", va="center", fontsize=14) + ax.set_xticks(range(len(names)), [n.split(":")[0] for n in names], rotation=90, fontsize=6) + ax.set_yticks(range(len(names)), names, fontsize=6) + ax.set_title(f"{sub} {sig}: {layer}", fontsize=9) + ax.set_xlabel("target", fontsize=7) + ax.set_ylabel("source", fontsize=7) + fig.suptitle("Recalled minus forgotten flow (bits); * significant (FDR 0.05); blank: not tested", fontsize=10) + fig.tight_layout() + FIGURES.mkdir(parents=True, exist_ok=True) + fig.savefig(FIGURES / f"catfr_sme_{sig}.png", dpi=130) + plt.close(fig) + + +def cohort_figure(summary_path): + """Per-subject recalled minus forgotten flow for each region pair and layer (``cohort.py``).""" + s = json.loads(Path(summary_path).read_text()) + edges = list(s["between"]) + fig, axes = plt.subplots(1, len(LAYERS), figsize=(4.5 * len(LAYERS), 3.6), sharey=False) + rng = np.random.default_rng(0) + for i, (ax, layer) in enumerate(zip(axes, LAYERS, strict=True)): + for j, e in enumerate(edges): + d = np.array(s["per_unit"][e]["recalled"])[:, i] - np.array(s["per_unit"][e]["forgotten"])[:, i] + ax.scatter(j + rng.uniform(-0.15, 0.15, len(d)), d, s=10, alpha=0.6) + ax.plot([j - 0.3, j + 0.3], [np.nanmean(d)] * 2, color="k") + ax.axhline(0, color="grey", lw=0.5) + ax.set_xticks(range(len(edges)), edges, rotation=45, fontsize=7) + q = [s["between"][e][layer]["qvalue"] for e in edges] + ax.set_title(f"{layer} (smallest q = {min(q):.2f})", fontsize=9) + ax.set_ylabel("recalled - forgotten (bits)", fontsize=7) + fig.suptitle(f"Cohort: {len(s['subjects'])} subjects, region-averaged HFA; line: mean", fontsize=10) + fig.tight_layout() + FIGURES.mkdir(parents=True, exist_ok=True) + fig.savefig(FIGURES / "catfr_cohort.png", dpi=130) + plt.close(fig) + + +def main(): + parser = argparse.ArgumentParser() + parser.add_argument("--cohort", help="cohort summary.json from cohort.py") + args = parser.parse_args() + summaries = load() + print(tables(summaries)) + figure(summaries) + if args.cohort: + cohort_figure(args.cohort) + + +if __name__ == "__main__": + main() diff --git a/examples/catfr/run.py b/examples/catfr/run.py new file mode 100644 index 0000000..b91f2b2 --- /dev/null +++ b/examples/catfr/run.py @@ -0,0 +1,248 @@ +""" +Subsequent-memory contrast of infoflow networks on catFR1 encoding epochs. + +For each prepared session and signal: + +1. Preprocess (discretize, choose embeddings) once on all of the session's + trials, so every condition and every null split shares one symbol encoding. +2. Balance recalled and forgotten words by subsampling the larger class, select + one network on both conditions together with ``infer_multiplex`` on the shared + encoding (plug-in selection at auto's level), and test the per-layer + differences on its edges with ``compare_networks(design="within")``. With + ``--selection condition`` the tested edges are instead the union of one network + per condition, which makes the test anti-conservative. +3. Calibration null: two disjoint same-condition draws from the larger class, + each the size of the main contrast's arms (capped at half the class), pushed + through the same steps ``--n-null`` times. +4. ``auto.infer`` on the balanced pooled trials, for its report. + + uv run --project ../.. --with netcdf4 python run.py --signal both +""" + +import argparse +import json +import time + +import numpy as np +import xarray as xr +from common import OUT + +from infoflow import infer_multiplex +from infoflow.auto import _skeleton_settings, infer +from infoflow.compare import compare_networks +from infoflow.embedding import Embedding +from infoflow.layers import LAYERS +from infoflow.parallel import thread_map +from infoflow.preprocess import preprocess + +DEFAULT_MAX_LAG = {"broadband": 40, "hfa": 15} + + +def condition_network(pre, idx, names, threads): + settings = _skeleton_settings(("plugin",), False, 1, None, {}) + return infer_multiplex( + pre.data.subset_trials(list(idx)), + names=names, + embeddings=pre.embeddings, + skeleton_settings=settings, + include_shared_candidates=False, + interpret=False, + map_fn=thread_map(threads), + prng=0, + ) + + +def contrast(pre, idx_a, idx_b, names, n_perm, threads, seed, selection="condition", curtail=10): + """ + ``selection='condition'`` tests the union of the per-condition skeletons, as + ``compare_networks`` is designed; ``'pooled'`` tests the skeleton selected on both + conditions together, which does not depend on the labels being permuted. + """ + data = pre.data + if selection == "pooled": + net_a = net_b = condition_network(pre, np.concatenate([idx_a, idx_b]), names, threads) + else: + net_a = condition_network(pre, idx_a, names, threads) + net_b = condition_network(pre, idx_b, names, threads) + cmp = compare_networks( + net_a, + net_b, + data.subset_trials(list(idx_a)), + data.subset_trials(list(idx_b)), + design="within", + n_perm=n_perm, + curtail=curtail, + map_fn=thread_map(threads), + prng=seed, + ) + return net_a, net_b, cmp + + +def layer_counts(net): + return {k: int(net.dataset["significant"].sel(layer=k).sum()) for k in LAYERS if k in net.dataset.layer} + + +def tested(cmp): + return {k: int(np.isfinite(cmp.dataset["pvalue"].sel(layer=k)).sum()) for k in LAYERS} + + +def significant_edges(cmp, regions): + names = list(cmp.dataset.source.values) + out = [] + for k in LAYERS: + sig = cmp.dataset["significant"].sel(layer=k).values + for s, t in zip(*np.nonzero(sig)): + out.append( + { + "layer": k, + "source": names[s], + "target": names[t], + "regions": f"{regions[s]}->{regions[t]}", + "recalled": float(cmp.dataset["a"].sel(layer=k).values[s, t]), + "forgotten": float(cmp.dataset["b"].sel(layer=k).values[s, t]), + "qvalue": float(cmp.dataset["qvalue"].sel(layer=k).values[s, t]), + } + ) + return out + + +def run_one(path, args): + ds = xr.open_dataset(path).load() + signal = ds.attrs["signal"] + stem = path.name.removesuffix(".nc") + regions = [str(r) for r in ds.region.values] + names = [f"{r}:{c}" for r, c in zip(regions, ds.channel.values)] + trials = [ds.signal.values[i].astype(float) for i in range(ds.sizes["trial"])] + recalled = ds.recalled.values.astype(bool) + rng = np.random.default_rng(0) + rec, forg = np.flatnonzero(recalled), np.flatnonzero(~recalled) + n = min(len(rec), len(forg)) + idx_rec = np.sort(rng.choice(rec, n, replace=False)) + idx_forg = np.sort(rng.choice(forg, n, replace=False)) + max_lag = args.max_lag or DEFAULT_MAX_LAG[signal] + summary = { + "file": path.name, + "signal": signal, + "sfreq": ds.attrs["sfreq"], + "n_per_condition": int(n), + "n_recalled": int(len(rec)), + "n_forgotten": int(len(forg)), + "max_lag": max_lag, + "names": names, + "timing": {}, + } + + t0 = time.perf_counter() + pre = preprocess( + trials, + names=names, + max_lag=max_lag, + network_resolution="auto", + map_fn=thread_map(args.threads), + prng=0, + ) + summary["timing"]["preprocess"] = time.perf_counter() - t0 + summary["preprocessing"] = pre.report.summary() + summary["lags"] = args.lags + if args.lags == "long": + pre.embeddings = [Embedding(max_lag=max_lag, log_spaced=True) for _ in names] + stem += "_long" + summary["selection"] = args.selection + if args.selection == "pooled": + stem += "_pooled" + summary["candidate_lags"] = sorted({lag for e in pre.embeddings for lag in e.lags()}) + print(f"[{stem}] preprocessed in {summary['timing']['preprocess']:.0f}s\n{summary['preprocessing']}", flush=True) + + t0 = time.perf_counter() + net_rec, net_forg, cmp = contrast(pre, idx_rec, idx_forg, names, args.n_perm, args.threads, 0, args.selection) + summary["timing"]["contrast"] = time.perf_counter() - t0 + summary["edges"] = {"recalled": layer_counts(net_rec), "forgotten": layer_counts(net_forg)} + summary["tested"] = tested(cmp) + summary["sme_edges"] = significant_edges(cmp, regions) + cmp.dataset.to_netcdf(OUT / f"{stem}_compare.nc") + for label, net in (("recalled", net_rec), ("forgotten", net_forg)): + net.dataset.drop_vars([v for v in net.dataset.data_vars if net.dataset[v].dtype == object]).to_netcdf( + OUT / f"{stem}_net_{label}.nc" + ) + print( + f"[{stem}] contrast in {summary['timing']['contrast']:.0f}s: edges {summary['edges']}, " + f"tested {summary['tested']}, significant SME {len(summary['sme_edges'])}", + flush=True, + ) + + major = rec if len(rec) >= len(forg) else forg + m = min(n, len(major) // 2) + null = [] + t0 = time.perf_counter() + for r in range(args.n_null): + draw = np.random.default_rng(1000 + r).permutation(major) + _, _, c = contrast( + pre, np.sort(draw[:m]), np.sort(draw[m : 2 * m]), names, args.n_perm, args.threads, 1000 + r, args.selection + ) + sig = {k: int(c.dataset["significant"].sel(layer=k).sum()) for k in LAYERS} + null.append({"significant": sig, "tested": tested(c)}) + print(f"[{stem}] null {r}: {sig}", flush=True) + summary["timing"]["null"] = time.perf_counter() - t0 + summary["null"] = {"class": "recalled" if major is rec else "forgotten", "n_per_split": int(m), "reps": null} + + path_out = OUT / f"{stem}_summary.json" + path_out.write_text(json.dumps(summary, indent=2, default=str)) + + if not args.skip_auto and args.lags == "auto": + t0 = time.perf_counter() + pooled = [trials[i] for i in np.concatenate([idx_rec, idx_forg])] + try: + result = infer( + pooled, + preset="balanced", + time_budget=args.time_budget, + names=names, + max_lag=max_lag, + threads=args.threads, + prng=0, + ) + except Exception as error: # noqa: BLE001 + summary["auto_error"] = f"{type(error).__name__}: {error}" + print(f"[{stem}] auto failed: {summary['auto_error']}", flush=True) + else: + summary["auto_report"] = str(result.report) + summary["auto_edges"] = layer_counts(result.network) + print(f"[{stem}] auto in {time.perf_counter() - t0:.0f}s\n{result.report}", flush=True) + summary["timing"]["auto"] = time.perf_counter() - t0 + path_out.write_text(json.dumps(summary, indent=2, default=str)) + return summary + + +def main(): + parser = argparse.ArgumentParser() + parser.add_argument("--signal", choices=("broadband", "hfa", "both"), default="both") + parser.add_argument("--subjects", nargs="*") + parser.add_argument("--max-lag", type=int) + parser.add_argument("--n-perm", type=int, default=500) + parser.add_argument("--n-null", type=int, default=20) + parser.add_argument("--threads", type=int) + parser.add_argument("--time-budget", default="2h") + parser.add_argument("--skip-auto", action="store_true") + parser.add_argument( + "--lags", + choices=("auto", "long"), + default="auto", + help="auto: preprocessing's lag budgets; long: log-spaced lags up to --max-lag", + ) + parser.add_argument( + "--selection", + choices=("condition", "pooled"), + default="pooled", + help="edges tested: union of per-condition skeletons, or the pooled skeleton", + ) + args = parser.parse_args() + signals = ("hfa", "broadband") if args.signal == "both" else (args.signal,) + for signal in signals: + for path in sorted(OUT.glob(f"sub-*_{signal}.nc")): + if args.subjects and not any(s in path.name for s in args.subjects): + continue + run_one(path, args) + + +if __name__ == "__main__": + main() diff --git a/infoflow/auto.py b/infoflow/auto.py index 63e0006..21c151c 100644 --- a/infoflow/auto.py +++ b/infoflow/auto.py @@ -86,6 +86,11 @@ _PILOT_TARGETS = 3 _MAX_RESAMPLES = 2000 _LAYER_ALPHA = 0.05 +# With a budget, checks and targets not started by these fractions of it are skipped or +# selected with the plug-in only: per-target selection times are heavy-tailed, so the +# pilot's prediction alone cannot keep a run within its budget. +_GATE_BUDGET_FRACTION = 0.5 +_SELECTION_BUDGET_FRACTION = 0.8 @dataclass @@ -295,16 +300,18 @@ def parent_prediction_error(data, target, variables, folds=5, k_fraction=0.02): def _gaussian_ratio(data, target, embeddings, plugin_settings, gaussian_settings, seed): """ The nearest-neighbour prediction error of the Gaussian estimator's parent set over the - plug-in's, each selected alone. + plug-in's, each selected alone, and the seconds the Gaussian selection took. """ def chosen(settings): sk = select_parents(data, target, embeddings, None, settings, seed) - return list(sk.target_past) + (list(sk.sources) if sk.significant else []) + return list(sk.target_past) + (list(sk.sources) if sk.significant else []), sk.timing - plugin = parent_prediction_error(data, target, chosen(plugin_settings)) - gaussian = parent_prediction_error(data, target, chosen(gaussian_settings)) - return gaussian / plugin if plugin > 0 else 1.0 + plugin_vars, _ = chosen(plugin_settings) + gaussian_vars, timing = chosen(gaussian_settings) + plugin = parent_prediction_error(data, target, plugin_vars) + gaussian = parent_prediction_error(data, target, gaussian_vars) + return (gaussian / plugin if plugin > 0 else 1.0), float(timing.get("gaussian", 0.0)) def linearity_pvalues(data, targets=None, max_lag=5, max_regressors=None, bins=4): @@ -623,6 +630,7 @@ def infer( # Gaussian estimator alone (local null), and keep the Gaussian estimator only where its # parents predict the target at least as well by a nonparametric referee. keep_gaussian, drop_gaussian = [], list(nonlinear) + gate_downgrades, gaussian_seconds = [], [] if nonlinear: pinned_base = {k: v for k, v in pinned_skeleton.items() if k != "estimator"} plugin_alone = _skeleton_settings(("plugin",), options["prescreen"], 1, device, pinned_base) @@ -633,11 +641,24 @@ def infer( def ratio(args): t, seed = args + if budget is not None and time.perf_counter() - started > _GATE_BUDGET_FRACTION * budget: + return t, None return t, _gaussian_ratio(discrete, t, embeddings, plugin_alone, gaussian_alone, int(seed)) - ratios = dict(thread_map(threads)(ratio, list(zip(nonlinear, choice_seeds, strict=True)))) - keep_gaussian = [t for t in nonlinear if ratios[t] < 1.0] - drop_gaussian = [t for t in nonlinear if ratios[t] >= 1.0] + ratios, skipped = {}, [] + for t, out in thread_map(threads)(ratio, list(zip(nonlinear, choice_seeds, strict=True))): + if out is None: + skipped.append(t) + else: + ratios[t] = out[0] + gaussian_seconds.append(out[1]) + if skipped: + gate_downgrades.append( + f"skipped the Gaussian check for {len(skipped)} nonlinear target(s) not started within " + f"{_GATE_BUDGET_FRACTION:.0%} of the budget; they select with the plug-in only" + ) + keep_gaussian = [t for t in nonlinear if ratios.get(t, 1.0) < 1.0] + drop_gaussian = [t for t in nonlinear if t not in keep_gaussian] diagnostics["gaussian_ratio"] = ratios diagnostics["gaussian_kept"] = keep_gaussian @@ -671,6 +692,9 @@ def pilot_run(t): name: float(np.mean([sk.timing[name] for sk in pilot_skeletons.values() if name in sk.timing] or [0.0])) for name in cheap["estimators"] } + if "gaussian" in member_seconds and member_seconds["gaussian"] == 0.0 and gaussian_seconds: + # No pilot target used the Gaussian estimator; the gate timed it (local null) on others. + member_seconds["gaussian"] = float(np.mean(gaussian_seconds)) if "ksg" in estimators and "ksg" not in cheap["estimators"]: t0 = pilot_targets[0] member_seconds["ksg"] = _ksg_seconds( @@ -700,9 +724,10 @@ def pilot_run(t): # Preprocessing, the linearity test, and the pilot are already spent. elapsed = time.perf_counter() - started remaining = None if budget is None else max(0.0, budget - elapsed) - final, report.downgrades, predicted = _fit_budget( + final, downgrades, predicted = _fit_budget( options, pilot, len(targets), map_threads(options), remaining, pinned_skeleton | pinned_infer ) + report.downgrades = gate_downgrades + downgrades if budget is not None and pinned_infer.get("machines"): report.warnings.append("the predicted time excludes the machine layer (machines=True runs outside the budget)") report.predicted_seconds = elapsed + predicted @@ -729,6 +754,31 @@ def pilot_run(t): if final_settings == pilot_settings and pinned_infer.get("holdout") is None: for t, sk in pilot_skeletons.items(): checkpoint.put(f"skeleton-{t}", sk) + if budget is not None and pinned_infer.get("holdout") is None: + # Select here, with the seeds infer_skeleton would use, so targets that start late + # can fall back to the plug-in; infer_multiplex then reads them from the checkpoint. + fallback = final_settings + if "estimator" not in overrides: + fallback = settings_for(dict(final, estimators=("plugin",)), map_threads(final)) + deadline = started + _SELECTION_BUDGET_FRACTION * budget + late = [] + + def select(t): + if f"skeleton-{t}" in checkpoint: + return t, checkpoint.get(f"skeleton-{t}") + settings = final_settings + if time.perf_counter() > deadline: + settings = fallback + late.append(t) + return t, select_parents(discrete, t, embeddings, None, settings, seed_of[t]) + + for t, sk in thread_map(map_threads(final))(select, targets): + checkpoint.put(f"skeleton-{t}", sk) + if late and fallback != final_settings: + report.downgrades.append( + f"selected {len(late)} target(s) with the plug-in only: started after " + f"{_SELECTION_BUDGET_FRACTION:.0%} of the budget ({', '.join(map(str, sorted(late)))})" + ) kwargs = { "embeddings": embeddings, "names": names, diff --git a/infoflow/compare.py b/infoflow/compare.py index f83356e..aaea8ea 100644 --- a/infoflow/compare.py +++ b/infoflow/compare.py @@ -93,6 +93,8 @@ def compare_networks( dependent=False, n_perm=200, alpha=0.05, + curtail=None, + map_fn=map, prng=None, ): """ @@ -102,6 +104,9 @@ def compare_networks( ---------- net_a, net_b : MultiplexNetwork Networks inferred under conditions A and B (their skeletons define the union). + Skeletons selected on each condition's own units make the test + anti-conservative; pass one network selected on both conditions pooled as + both arguments for an exact permutation test. data_a, data_b : DiscreteData or list of DiscreteData For ``design='within'``, one DiscreteData per condition whose trials are the units. For ``design='between'``, a list of per-subject DiscreteData. @@ -111,6 +116,15 @@ def compare_networks( n_perm : int alpha : float FDR level per layer. + curtail : int, None + Sequential Monte Carlo p-values :cite:`Besag1991`: an edge's permutations stop + once every layer's permuted difference has reached the observed one `curtail` + times, and a layer that reached it after ``L`` permutations gets p-value + ``curtail / L``. Still valid, and edges with large p-values stop after a few + dozen permutations instead of `n_perm`. None runs all `n_perm`. + map_fn : callable + ``map``-like function over edges (e.g. from :mod:`infoflow.parallel`). Each + edge's permutations use a seed drawn up front, so results do not depend on it. prng : None, int, Generator Returns @@ -148,24 +162,39 @@ def statistic(lab, e): b = _edge_values(_pool([u for u, g in zip(units, lab, strict=True) if g == 1]), e[1], *edges[e]) return np.array([a[k] for k in LAYERS]), np.array([b[k] for k in LAYERS]) - def permuted(): + def permuted(edge_rng): if dependent: n = len(units_a) - swap = rng.random(n) < 0.5 + swap = edge_rng.random(n) < 0.5 lab = labels.copy() lab[:n][swap], lab[n:][swap] = 1, 0 return lab - return rng.permutation(labels) + return edge_rng.permutation(labels) + + def test(job): + e, seed = job + edge_rng = np.random.default_rng(seed) + a, b = statistic(labels, e) + observed = np.abs(a - b) - 1e-12 + exceed = np.zeros(len(LAYERS), dtype=np.int64) + reached = np.zeros(len(LAYERS), dtype=np.int64) + for i in range(1, n_perm + 1): + exceed += np.abs(np.subtract(*statistic(permuted(edge_rng), e))) >= observed + if curtail is not None: + reached[(reached == 0) & (exceed >= curtail)] = i + if reached.all(): + break + p = (1 + exceed) / (1 + n_perm) + if curtail is not None: + p = np.where(reached > 0, curtail / np.maximum(reached, 1), p) + return e, a, b, np.where(np.isnan(observed), np.nan, p) shape = (len(LAYERS), P, P) va, vb, diff, pvalue = (np.full(shape, np.nan) for _ in range(4)) - for e in edges: - a, b = statistic(labels, e) - observed = a - b - null = np.array([np.subtract(*statistic(permuted(), e)) for _ in range(n_perm)]) - s, t = e - va[:, s, t], vb[:, s, t], diff[:, s, t] = a, b, observed - pvalue[:, s, t] = (1 + np.sum(np.abs(null) >= np.abs(observed) - 1e-12, axis=0)) / (1 + n_perm) + seeds = rng.integers(0, 2**32, size=len(edges)) + for (s, t), a, b, p in map_fn(test, list(zip(edges, (int(x) for x in seeds), strict=True))): + va[:, s, t], vb[:, s, t], diff[:, s, t] = a, b, a - b + pvalue[:, s, t] = p significant, adjusted = adjust_layer_pvalues({k: pvalue[i] for i, k in enumerate(LAYERS)}, alpha) ds = xr.Dataset( { @@ -177,6 +206,6 @@ def permuted(): "significant": (("layer", "source", "target"), np.stack([significant[k] for k in LAYERS])), }, coords={"layer": list(LAYERS), "source": names, "target": names}, - attrs={"design": design, "dependent": int(dependent), "n_perm": n_perm}, + attrs={"design": design, "dependent": int(dependent), "n_perm": n_perm, "curtail": curtail or 0}, ) return NetworkComparison(ds, edges) diff --git a/infoflow/preprocess/pipeline.py b/infoflow/preprocess/pipeline.py index 2ed1b4c..a3862ea 100644 --- a/infoflow/preprocess/pipeline.py +++ b/infoflow/preprocess/pipeline.py @@ -542,7 +542,8 @@ def _network_resolution(trials, results, discretizers, embeddings, settings, mod network = {"bins": B, "history": L, "scores": scores, "retained": retained, "adopted": bool(adopted)} if adopted: discretizers = [EqualFrequency(B) for _ in discretizers] - embeddings = [Embedding(max_lag=L) for _ in embeddings] + # The uniform history sets the context resolution, not how far back sources may act. + embeddings = [Embedding(max_lag=max(L, e.max_lag), tau=e.tau) for e in embeddings] return discretizers, embeddings, network diff --git a/infoflow/preprocess/scoring.py b/infoflow/preprocess/scoring.py index fc9059c..aa790e9 100644 --- a/infoflow/preprocess/scoring.py +++ b/infoflow/preprocess/scoring.py @@ -63,6 +63,8 @@ def _fold_scores(contexts, targets, B, folds, gap, alpha=0.5): ctx = np.concatenate([c for c, _ in contexts]) tgt = np.concatenate([targets[i][t] for i, (_, t) in enumerate(contexts)]) trial = np.concatenate([np.full(len(c), i) for i, (c, _) in enumerate(contexts)]) + if len(ctx) == 0: + return np.array([]) _, ctx = np.unique(ctx, return_inverse=True) ctx = ctx.ravel() n = len(ctx) @@ -143,7 +145,8 @@ def score_contexts(contexts, raw_series, resolutions=(2, 4, 8), folds=5, gap=1, per[B] = _fold_scores(contexts, reference_codes(raw_series, B), B, folds, gap) k = min(len(v) for v in per.values()) if k == 0: - return 0.0, np.inf, per, np.zeros(0) + # No fold has enough held-out windows (e.g. the context spans most of each trial). + return 0.0, np.inf, {B: (0.0, v) for B, v in per.items()}, np.zeros(0) fold_scores = np.mean([v[:k] for v in per.values()], axis=0) se = float(np.std(fold_scores, ddof=1) / np.sqrt(k)) if k > 1 else np.inf return float(np.mean(fold_scores)), se, {B: (float(np.mean(v)), v) for B, v in per.items()}, fold_scores diff --git a/reports/catfr_pilot.md b/reports/catfr_pilot.md new file mode 100644 index 0000000..e367d03 --- /dev/null +++ b/reports/catfr_pilot.md @@ -0,0 +1,259 @@ +# infoflow on catFR1: a subsequent-memory pilot and cohort + +*Branch `catfr-pilot`, October 2026. Scripts: `examples/catfr/{fetch,prepare,run,cohort,null_check,report}.py`; +library changes on `catfr-library-fixes`. Tables are printed by `report.py` from the +per-run summaries.* + +## Summary + +- **The dataset's event onsets are misaligned; the first pilot runs used them.** The + BIDS `onset` column of ds004809 is offset from the recording by 1 s to 15 min per + session (drifting slowly within a session); the `sample` column is aligned. Word-locked + evoked responses appear only with `sample / sfreq` (evoked power after onset 7.5x and + 4.3x the pre-stimulus level in two subjects, against 0.8x with `onset`). Every + memory-effect result from the first runs is superseded by the corrected runs below; + their calibration and engineering findings do not depend on alignment and stand. +- **No subsequent-memory effect in directed flow.** On correctly aligned epochs, the + three-subject channel-level pilot found one significant edge in six contrasts (HFA, + shared flow between two MTL contacts in R1065J, q = 0.036) and the 25-subject cohort, + with one node per region (MTL, lateral temporal, prefrontal), found no region pair and + layer whose recalled-minus-forgotten flow differs across subjects (18 tests, smallest + q = 0.48). The same cohort data show the classical amplitude effect (prefrontal HFA + higher for recalled words in 18 of 25 subjects, p = 0.001), so the data and labels are + sound and the null is specific to flow. +- **Comparisons need pooled selection.** Testing the union of edges selected in each + condition gave significant differences between random halves of one condition in + 13-15% of splits instead of 5%. Selecting the tested edges on both conditions pooled + makes the permutation test exact: intrinsic and shared flow come back to nominal. In a + 100-split check, synergistic flow's per-edge p-values are mildly low between 0.01 and + 0.05 (9.4% below 0.05, uniform below 0.01); modest evidence, from overlapping splits. +- **Library changes** (separate PR): preprocessing no longer crashes on short trials and + keeps each node's lag budget under the network resolution; `compare_networks` gains + sequential p-values (7.7x faster on the pilot's HFA contrast, same decisions), a map + over edges, and pooled-selection guidance; `auto.infer` enforces its time budget with + deadlines instead of trusting a three-target pilot. + +## 1. Data and pipeline + +**Dataset.** OpenNeuro [ds004809](https://openneuro.org/datasets/ds004809) (Herrema & +Kahana; catFR1, CC0): BIDS iEEG with bipolar montages, atlas-labelled contacts and +behavioural events, about 250 subjects. + +- *Event timing.* Epochs are cut at `sample / sfreq`. The `onset` column differs from it + by a session-specific offset (R1364C 1.2 s, R1332M 1.5 s, R1065J 56.6 s, R1108J 228 s, + R1004D 852 s; one session, R1293P, has onsets 23 h into a 68 min recording) that drifts + by tens of milliseconds over a session, as if `onset` came from the behavioural clock. + `common.word_events` documents this. +- *Recall labels.* The events carry no `recalled` column; a `WORD` is recalled when a + `REC_WORD` in the same list names the same item. Practice lists are dropped. +- *No channel exclusion.* The release has no bad-channel or seizure-onset flags. + +**Channels.** A bipolar pair keeps a region only if both contacts agree (hippocampal +subfield atlases first, then whole-brain and Desikan-Killiany labels); pairs sharing a +contact are never both used, since they share a signal term (an early smoke test's only +significant edges joined such pairs). + +**Signals.** 60 Hz notch and harmonics; WORD epochs 0-1.6 s after onset: + +- *broadband*: 1-80 Hz at 200 Hz (320 samples per epoch); +- *HFA*: log mean Morlet power (7 cycles, 8 log-spaced frequencies in 70-150 Hz), + computed on epochs padded by 1 s and cropped, averaged down to 50 Hz (80 samples). + +**Pilot** (`run.py`): three subjects (R1065J, R1364C, R1332M), six pairs per region +(18 nodes), 75-130 trials per condition after balancing by subsampling the larger class. +Preprocessing runs once on all of a session's trials so every condition and split shares +one symbol encoding; one network is selected on both conditions pooled +(`infer_multiplex`, plug-in selection at `auto`'s 0.001 level); `compare_networks` +(within subjects, 500 permutations, `curtail=10`, FDR 0.05 per layer) tests each layer's +recalled-minus-forgotten difference on its edges. Calibration null: random halves of the +larger class (recalled for all three), the size of the main contrast's arms. + +**Cohort** (`cohort.py`): one session for each of 25 subjects with at least two +non-overlapping pairs in each region, ranked by trials per condition after balancing +(the 25 chosen have 95-149). Each region's HFA, each pair z-scored first, +is averaged into one node. Every subject is discretized with the same 4 equal-frequency +bins, and `compare_networks(design="between", dependent=True, n_perm=5000)` tests the +paired recalled-minus-forgotten difference across subjects on all six directed region +pairs, with variables fixed in advance (source at lag 2 = 40 ms, target at lags 1-2, +third region at lag 2), so nothing is selected on the data. Each subject also gets the +pilot's within-subject analysis on its three nodes. EDFs are downloaded, epoched and +deleted one at a time (17.6 GB in total). + +## 2. Results + +### Pilot (channel level, correctly aligned) + +| Subject | Signal | Trials / condition | Pooled edges (I/S/Sh) | Significant SME (I/S/Sh) | Null splits with any significant (I/S/Sh) | +|---|---|---|---|---|---| +| R1065J | hfa | 75 | 1/0/1 | 0/0/1 | 0/0/0 of 20 | +| R1332M | hfa | 130 | 2/0/2 | 0/0/0 | 1/2/0 of 20 | +| R1364C | hfa | 79 | 5/0/6 | 0/0/0 | 0/1/1 of 20 | +| R1065J | broadband | 75 | 0/1/3 | 0/0/0 | 0/0/0 of 3 | +| R1332M | broadband | 130 | 3/7/18 | 0/0/0 | 0/0/0 of 3 | +| R1364C | broadband | 79 | 1/4/18 | 0/0/0 | 0/0/1 of 3 | + +The one significant edge is shared flow from MTL:LC1-LC2 to MTL:LD1-LD2 in R1065J (two +shafts): 0.023 bits for recalled words against 0.012 for forgotten, q = 0.036. Nothing +comparable appears in the other subjects or in broadband. HFA networks are small (1-6 +edges per layer from 18 nodes) and carry intrinsic and shared flow; broadband networks +are mostly shared and synergistic flow (up to 18 shared edges), the signature of common +input or residual volume conduction rather than directed transmission. + +![Pilot, HFA](figures/catfr_sme_hfa-pooled.png) + +**Label check.** Mean HFA 0.4-1.6 s after onset, recalled vs forgotten (Welch t averaged +over each region's pairs): MTL +1.6, +1.2, +2.5 (R1065J, R1332M, R1364C), with 7, 3 and 7 +of 18 pairs at p < 0.05, in line with the HFA subsequent-memory literature. With the +misaligned epochs the same check gave 0-2 pairs. + +### Cohort (region level, 25 subjects) + +| Edge | Intrinsic: rec / forg (bits), q | Synergistic: rec / forg, q | Shared: rec / forg, q | +|---|---|---|---| +| LTC -> MTL | 0.0001 / 0.0002, 0.56 | 0.0069 / 0.0075, 0.72 | 0.0036 / 0.0040, 0.86 | +| PFC -> MTL | 0.0002 / 0.0002, 0.82 | 0.0060 / 0.0074, 0.48 | 0.0007 / 0.0009, 0.67 | +| MTL -> LTC | 0.0002 / 0.0001, 0.56 | 0.0061 / 0.0072, 0.48 | 0.0035 / 0.0037, 0.86 | +| PFC -> LTC | 0.0001 / 0.0002, 0.56 | 0.0067 / 0.0070, 0.77 | 0.0034 / 0.0027, 0.67 | +| MTL -> PFC | 0.0002 / 0.0001, 0.56 | 0.0056 / 0.0059, 0.77 | 0.0010 / 0.0007, 0.67 | +| LTC -> PFC | 0.0002 / 0.0001, 0.56 | 0.0061 / 0.0068, 0.72 | 0.0033 / 0.0028, 0.67 | + +Means over subjects; q is FDR-adjusted per layer over the six edges. In every cell +between 4 and 13 of the 25 subjects have recalled above forgotten, so there is no +consistent direction either. Within subjects, pooled selection on the three nodes found +0-4 edges per layer (none in 13 subjects) and one significant difference in 25 subjects (R1201P, +shared MTL to LTC, q = 0.048), about what chance gives. + +![Cohort](figures/catfr_cohort.png) + +Over the same 25 subjects and nodes, the classical amplitude effect (mean HFA 0.4-1.6 s, +recalled minus forgotten) is present in PFC (18/25 positive, t = 3.60, p = 0.001), leans +the same way in MTL (16/25, p = 0.16), and is absent in LTC (13/25, p = 0.66). + +## 3. Calibration of the comparison + +Random halves of one condition have no true difference, however much flow of any kind +the data contain, so every significant result in these splits is a false positive. With +the edges selected label-independently, the permutation test is exact per edge. + +| Variant | Null splits | Any significant: intrinsic | synergistic | shared | +|---|---|---|---|---| +| Per-condition selection (misaligned HFA) | 60 | 8 (13%, p = 0.010) | 9 (15%, p = 0.003) | 6 (10%, p = 0.079) | +| Pooled selection (misaligned HFA) | 60 | 4 (7%, p = 0.35) | 8 (13%, p = 0.010) | 4 (7%, p = 0.35) | +| Pooled selection, R1332M, 100 splits | 100 | 6 (6%, p = 0.38) | 10 (10%, p = 0.028) | 3 (3%, p = 0.88) | +| Pooled selection, `curtail=10` (aligned HFA) | 60 | 1 (2%) | 3 (5%) | 1 (2%) | +| Pooled selection, `curtail=10` (aligned broadband) | 9 | 0 | 0 | 1 | + +(p: one-sided binomial test against 5%.) + +- **Per-condition selection inflates** because an edge strong by chance in one half is + selected for that reason, and permuting the labels destroys the chance alignment. + `compare_networks`'s docstring and `docs/workflow.rst` now say so and show the pooled + pattern. +- **The synergistic excess** (`null_check.py`, R1332M, 100 splits, 170 edge tests, + per-edge raw p-values recorded): 9.4% of synergistic p-values fall below 0.05 but 1.8% + below 0.01, and a Kolmogorov-Smirnov test does not reject uniformity (p = 0.43); with + about 1.7 edges per split, the per-split rate is the per-edge rate, so the FDR step is + not the cause. The excess is mild, confined to p between 0.01 and 0.05, and the splits + overlap heavily (all from the same 170 trials), so the binomial p overstates the + evidence. Intrinsic (6.5% below 0.05) and shared (4.7%) are nominal. Splits 0-19 + reproduce the earlier run exactly. + +## 4. Library changes + +All on `catfr-library-fixes`, rebased on `main` (after the spike-train merge); the full +test suite passes (281 passed, 2 skipped, plus the new `auto` tests). + +1. **Short-trial crash** (`preprocess/scoring.py`). A candidate whose context leaves no + fold with enough held-out windows returned its per-resolution scores in a different + shape, and one spanning the whole trial failed on an empty array; both now score 0 and + are never chosen. `auto.infer` now runs on 80-sample epochs; two regression tests. +2. **Lag budgets under the network resolution** (`preprocess/pipeline.py`). Adoption kept + only the uniform history (often one sample) as every node's lag budget; nodes now keep + their own budgets (at least that history). On this data candidate lags became 1-5 + samples for HFA and up to 10 for broadband, instead of 1. +3. **`compare_networks`** (`compare.py`). `curtail` gives sequential Monte Carlo p-values + (Besag & Clifford 1991): a layer's permutations stop once its permuted difference has + reached the observed one `curtail` times. On R1364C's HFA contrast (18 tested + edge-layers), `curtail=10` with 4 threads took 9 s against 69 s for all 500 + permutations serially, with identical decisions. `map_fn` runs edges in parallel with + seeds drawn up front, and a NaN observed difference now gives a NaN p-value instead of + the smallest attainable one. Pooled-selection guidance as above. +4. **`auto.infer` time budget** (`auto.py`). Two causes of the broadband overrun: the + check of nonlinear targets ran a Gaussian selection with a local permutation null per + target before the pilot, outside the budget; and the pilot's three targets could not + represent selection times that ranged from 47 s to nearly an hour per target on the + same data, and predicted zero Gaussian cost when none of them used the Gaussian + estimator. Now checks not started within half the budget are skipped, targets not + started within 80% of it select with the plug-in only (both recorded as downgrades), + and the pilot takes the Gaussian time from the check when it has none. A target + started just before a deadline still runs to completion, so a run can exceed its + budget by up to one target's selection. On R1065J's 150 broadband epochs, a 2 h budget + had overrun by more than an hour before (stopped unfinished after 3.1 h, then 3.2 h + with only the round-based check); with the deadlines, the targets started after 80% + of the budget took about 10 s each with the plug-in, and the last Gaussian selection + that started before it finished at 1.9 h. With a deliberately infeasible 30 min + budget, every downgrade fired (three checks skipped, pre-screen on, all 18 targets + plug-in only) and the run took 97 min: about 16 min of preprocessing, which runs + before any budget logic, and checks already in flight at the deadline. A hard stop + would need cancellation inside `select_parents`. + +## 5. Runtime + +16-core laptop shared with other jobs (load 15-25); 4 worker threads per run. + +| Run | Preprocess | Contrast (selection + comparison) | Null splits | auto.infer (2 h budget) | +|---|---|---|---|---| +| R1065J HFA | 7 min | 3 min | 44 min (20) | 45 min (predicted 15) | +| R1332M HFA | 7 min | 4 min | 31 min (20) | 49 min (predicted 37) | +| R1364C HFA | 7 min | 6 min | 60 min (20) | 43 min (predicted 15) | +| R1065J broadband | 16 min | 8 min | 23 min (3) | | +| R1332M broadband | 16 min | 91 min | 91 min (3) | | +| R1364C broadband | 17 min | 88 min | 3.9 h (3) | | +| Cohort, 25 subjects | | | | about 2 h end to end, including downloads | + +With `curtail`, the comparison is no longer the bottleneck; broadband contrast time is +now mostly the pooled selection on 160-260 trials. + +## 6. Caveats + +- One session per subject, montages that differ between subjects, and no seizure-onset + or bad-channel exclusion. +- The cohort's fixed variable set (one source lag, 40 ms) is one pre-specified choice; + other lags were not tested, to avoid selecting on the data. +- Region averaging trades spatial specificity for power and comparability; it dilutes + the MTL amplitude effect that the channel-level check shows. +- Bipolar referencing and the shared-contact rule reduce but do not remove volume + conduction; broadband shared flow is likely partly that. +- A 1.6 s encoding window is short for directed estimates, and pooling it across trials + assumes stationarity the nonstationary flags dispute. + +## 7. Next steps + +1. Report the onset misalignment to the dataset maintainers. +2. Look at synergistic flow's p-values between 0.01 and 0.05 directly: the synergy + estimate is a difference of two estimates, one from a numerical optimisation, and its + permutation distribution may be lumpy. +3. Cooperative cancellation in `select_parents` (a deadline checked between permutation + batches), so `auto.infer` can stop in-flight selections instead of waiting for them, + and a budget check that counts preprocessing. +4. If directed flow is still of interest: longer windows (encoding plus the following + item), more subjects, or a pre-specified hypothesis (e.g. MTL to PFC at several lags) + rather than all pairs. + +## Reproducing + +```bash +cd examples/catfr +uv run --project ../.. python fetch.py --top 3 # needs the aws CLI +uv run --project ../.. --with mne --with netcdf4 --with "decorator<5.3" python prepare.py +export OMP_NUM_THREADS=1 +uv run --project ../.. --with netcdf4 python run.py --signal hfa --n-null 20 +uv run --project ../.. --with netcdf4 python run.py --signal broadband --n-null 3 --skip-auto +uv run --project ../.. --with netcdf4 python null_check.py --subject R1332M --start 0 --count 100 +CATFR_OUT=~/data/catfr_cohort uv run --project ../.. --with mne --with netcdf4 --with "decorator<5.3" \ + python cohort.py --subjects 25 +uv run --project ../.. --with netcdf4 --with matplotlib python report.py --cohort ~/data/catfr_cohort/cohort/summary.json +``` + +`CATFR_ROOT` and `CATFR_OUT` override the data and output directories. `decorator` +5.3.1 is broken on PyPI (it installs without its module), hence the pin for `mne`. diff --git a/reports/figures/catfr_cohort.png b/reports/figures/catfr_cohort.png new file mode 100644 index 0000000..6dde631 Binary files /dev/null and b/reports/figures/catfr_cohort.png differ diff --git a/reports/figures/catfr_sme_broadband-pooled.png b/reports/figures/catfr_sme_broadband-pooled.png new file mode 100644 index 0000000..decf2b8 Binary files /dev/null and b/reports/figures/catfr_sme_broadband-pooled.png differ diff --git a/reports/figures/catfr_sme_hfa-pooled.png b/reports/figures/catfr_sme_hfa-pooled.png new file mode 100644 index 0000000..03537a4 Binary files /dev/null and b/reports/figures/catfr_sme_hfa-pooled.png differ diff --git a/tests/test_auto.py b/tests/test_auto.py index 8a33adb..f68d9d4 100644 --- a/tests/test_auto.py +++ b/tests/test_auto.py @@ -260,6 +260,41 @@ def test_infer_gates_gaussian_per_target(): assert all(r >= 1.0 for r in report.diagnostics["gaussian_ratio"].values()) +def test_gaussian_gate_stops_within_the_budget(monkeypatch): + calls = [] + + def ratio(data, target, *args): + calls.append(target) + return 2.0, 5.0 + + monkeypatch.setattr(auto, "_gaussian_ratio", ratio) + monkeypatch.setattr(auto, "_GATE_BUDGET_FRACTION", 0.0) + result = auto.infer( + datasets.mute_network(3000, seed=0), preset="fast", max_lag=5, threads=1, time_budget="1h", prng=0 + ) + report = result.report + assert report.diagnostics["nonlinear_targets"] == [1, 3] + assert calls == [] and report.diagnostics["gaussian_ratio"] == {} + assert report.downgrades[0].startswith("skipped the Gaussian check for 2 nonlinear target(s)") + # Unchecked targets select with the plug-in only, like targets that failed the check. + assert result.network.settings["skeleton_settings"].target_estimators == {1: ("plugin",), 3: ("plugin",)} + + +def test_late_targets_select_with_the_plugin_within_the_budget(monkeypatch): + data = datasets.mute_network(3000, seed=0) + on_time = auto.infer(data, preset="fast", max_lag=5, threads=1, time_budget="1h", prng=0) + assert not any("plug-in only: started after" in d for d in on_time.report.downgrades) + monkeypatch.setattr(auto, "_SELECTION_BUDGET_FRACTION", 0.0) + late = auto.infer(data, preset="fast", max_lag=5, threads=1, time_budget="1h", prng=0) + pilot = set(late.report.pilot["targets"]) + message = next(d for d in late.report.downgrades if "plug-in only: started after" in d) + assert f"selected {len(late.network.skeleton) - len(pilot)} target(s)" in message + assert all(name == "plugin" for (s, t), name in late.report.admitted.items() if t not in pilot) + # Without a budget the selection is unchanged. + free = auto.infer(data, preset="fast", max_lag=5, threads=1, prng=0) + assert np.array_equal(free.network.dataset["significant"].values, on_time.network.dataset["significant"].values) + + @pytest.mark.slow def test_presets_use_network_resolution_on_lorenz96(): from infoflow.benchmarks import lorenz96_graph diff --git a/tests/test_engineering.py b/tests/test_engineering.py index 7812672..87cb14a 100644 --- a/tests/test_engineering.py +++ b/tests/test_engineering.py @@ -49,6 +49,32 @@ def test_compare_within_no_difference(): assert np.nanmax(np.abs(result.dataset["difference"].values)) < 0.03 +def test_compare_curtailed_and_parallel(): + a, b = _coupled_trials(0.9, 10, 300, seed=0), _coupled_trials(0.5, 10, 300, seed=1) + net_a, net_b = _net(a), _net(b) + full = compare_networks(net_a, net_b, a, b, n_perm=99, prng=0) + curtailed = compare_networks(net_a, net_b, a, b, n_perm=99, curtail=5, prng=0) + parallel = compare_networks(net_a, net_b, a, b, n_perm=99, curtail=5, map_fn=thread_map(4), prng=0) + assert np.array_equal(curtailed.dataset["pvalue"].values, parallel.dataset["pvalue"].values, equal_nan=True) + assert np.array_equal(full.dataset["significant"].values, curtailed.dataset["significant"].values) + p = curtailed.dataset["pvalue"].values + # A curtailed p-value is curtail / L for some L <= n_perm. + tested = np.isfinite(p) & (p > 0.06) + assert np.allclose(5 / p[tested], np.round(5 / p[tested])) + + +def test_compare_within_null_is_calibrated_with_pooled_selection(): + pooled = _coupled_trials(0.8, 40, 200, seed=7) + halves = [] + for seed in range(20): + idx = np.random.default_rng(seed).permutation(40) + a, b = pooled.subset_trials(list(idx[:20])), pooled.subset_trials(list(idx[20:])) + net = _net(pooled) + result = compare_networks(net, net, a, b, n_perm=39, curtail=5, prng=seed) + halves.append(bool(result.dataset["significant"].values.any())) + assert sum(halves) <= 3 + + def test_compare_between_dependent(): subjects_a = [_coupled_trials(0.9, 2, 400, seed=s) for s in range(6)] subjects_b = [_coupled_trials(0.4, 2, 400, seed=100 + s) for s in range(6)] diff --git a/tests/test_preprocess.py b/tests/test_preprocess.py index afbdbcb..5d7c563 100644 --- a/tests/test_preprocess.py +++ b/tests/test_preprocess.py @@ -125,7 +125,8 @@ def test_network_resolution_adopted_only_when_every_node_keeps_its_score(): net = lorenz.report.network assert net["adopted"] and min(net["retained"]) >= 0.7 assert {str(d) for d in lorenz.data.discretizers} == {f"EqualFrequency(bins={net['bins']})"} - assert all(e.max_lag == net["history"] for e in lorenz.embeddings) + budgets = [n.lag_budget for n in lorenz.report.nodes] + assert [e.max_lag for e in lorenz.embeddings] == [min(2, max(net["history"], b)) for b in budgets] assert "network resolution" in lorenz.report.summary() mute = preprocess(datasets.mute_network(3000, seed=0), max_lag=5, network_resolution="auto", prng=0) @@ -133,3 +134,28 @@ def test_network_resolution_adopted_only_when_every_node_keeps_its_score(): forced = preprocess(datasets.mute_network(3000, seed=0), max_lag=5, network_resolution=True, prng=0) assert forced.report.network["adopted"] assert preprocess(datasets.mute_network(500, seed=0), max_lag=3, prng=0).report.network == {} + + +def test_network_resolution_keeps_each_nodes_lag_budget(): + forced = preprocess(datasets.mute_network(3000, seed=0), max_lag=8, network_resolution=True, prng=0) + history = forced.report.network["history"] + budgets = [min(8, max(history, n.lag_budget)) for n in forced.report.nodes] + assert [e.max_lag for e in forced.embeddings] == budgets + assert max(budgets) > history + + +def test_score_contexts_without_scorable_folds(): + from infoflow.preprocess.scoring import context_codes, score_contexts + + series = [np.random.default_rng(i).normal(size=12) for i in range(3)] + symbols = [np.searchsorted([0.0], s) for s in series] + mean, se, per, folds = score_contexts(context_codes(symbols, [1, 2], 2), series) + assert mean == 0.0 and np.isinf(se) and len(folds) == 0 + assert all(v[0] == 0.0 for v in per.values()) + + +def test_preprocess_short_trials(): + rng = np.random.default_rng(0) + trials = [np.cumsum(rng.normal(size=(80, 2)), axis=0) for _ in range(60)] + result = preprocess(trials, max_lag=5, prng=0) + assert result.data.n_trials == 60 and len(result.embeddings) == 2