diff --git a/CHANGELOG.md b/CHANGELOG.md index fe2ac9d..010b60c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,41 @@ found 24 defects plus a dozen false README claims. All verified and fixed here. ### Behaviour changes (read these before re-running a study) +- **Atlases resolve by name to their own files.** `pipeline.atlas` is looked up in + source-localization's `registry.yaml`, so allen26 and allen64 no longer pick up + allen32's labels, mapping and `roi_categories.yaml` from the shared `allen/` + directory. For allen26 studies this changes every region-level table + (Frontal-Anterior and Olfactory were silently dropped; Deep Subcortical was built + from 4 of its 8 parcels) and every ROI mosaic (the six merged parcels drew blank). + ROI-level results are unchanged, and allen32/antwerp studies resolve to the same + files as before. An atlas name that cannot be resolved is now an error, not a guess. +- **The R region tier uses the study's categories.** Every ROI R entry point replaced + the study's (or a profile's) `roi_categories` with the atlas-directory file whenever + one existed. Python now hands R the effective map and R prefers it + (`resolve_roi_categories` in `stats_utils.R`); the file is only a fallback for a + config that carries none. +- **The 10x voxel convention is read from the NIfTI header**, as source-localization + does, not guessed from the filename. `Atlas_3DRoisLeftRight.Labels.nii` has stored + true units since source-localization 2026-03-12 but was still shrunk 10x on the + default-affine path. **No statistic changes**: ROI extraction and cluster/NBS region + labels use the raw affine, which was always right. The only visible effect is + cosmetic: the mm axis ranges and slice labels of `plot_brain_roi_mosaic` / + `plot_brain_roi` when drawn on Antwerp, and no analysis module draws them on Antwerp. + (`load_vertex_roi_labels`, the other default-affine consumer, has never run: it reads + the mapping file's top-level keys as label ids and raises on every atlas, and + `vertex_network` swallows the error and falls back to spatial node labels. Left for + the vertex split.) +- **`electrode_signature` compares only within its own paradigm**, preferring + `roi_signature` over `vertex_signature`. It used to take the first + `vertex_signature_results.csv` anywhere under the results tree, which can be a + stale table from another run. `signature_source_vs_sensor.csv` gains a + `source_module` column. +- **Signature fits run single-threaded** (`threadpoolctl`, installed with scikit-learn). + A run fits a tiny model LOOCV x (1 + n_permutations) times, and a multithreaded BLAS + spent that time synchronising threads: logistic fits ran 60-90x slower (FORGE + electrode features: 105-128 s vs 1.4-1.7 s per 21 LOOCV passes, same accuracy). + Results are unchanged; a signature module that took days now takes about an hour. + - **Vertex `absolute` band power is now a density (dB/Hz)**, `10*log10(integral / bandwidth)`, matching the ROI/electrode definition. Previously `vertex_cluster` / `vertex_specparam` reported `10*log10(integral)`. Within-band group statistics are unaffected (a per-band @@ -35,6 +70,14 @@ found 24 defects plus a dozen false README claims. All verified and fixed here. ### Added +- **`roi_signature`**: ROI-level neural signature (decoding on per-parcel relative band + power), the source-side counterpart of `electrode_signature` with the identical + feature estimator. It needs no vertex estimate, so it runs on any ROI output, + including Monte Carlo operators. `sensor_paradigm:` compares it against an + `electrode_signature` run in another paradigm. +- **`resolve_atlas` / `AtlasSpec`**, and `atlas_files:` in the study config for atlases + that are not registered. Also `header_is_inflated` and `registered_atlases`. + - **`_subnetwork_edges.csv`** next to `roi_nbs_hypotheses.csv` (ROI edge modules): one row per supra-threshold edge of every NBS component (`hypothesis, band, dv, component_id, component_p, significant, node_i, node_j, roi_i, roi_j, stat`). The diff --git a/R/roi_aperiodic_analysis.R b/R/roi_aperiodic_analysis.R index 44b6e98..35a229c 100644 --- a/R/roi_aperiodic_analysis.R +++ b/R/roi_aperiodic_analysis.R @@ -119,12 +119,8 @@ if (!is.null(args$hypothesis) && length(diag_contrasts) > 0) { message("Study: ", config$name) message("Groups: ", paste(group_order, collapse = ", ")) -# Load roi_categories from atlas file if provided -if (!is.null(args$roi_categories) && file.exists(args$roi_categories)) { - config$roi_categories <- read_yaml(args$roi_categories) - message("Loaded roi_categories from: ", args$roi_categories, - " (", length(config$roi_categories), " regions)") -} +# Study-config categories win; the atlas file is only a fallback (stats_utils.R). +config$roi_categories <- resolve_roi_categories(config$roi_categories, args$roi_categories) # ============================================================ # Aperiodic-specific LMM functions (no band dimension) diff --git a/R/roi_connectivity_analysis.R b/R/roi_connectivity_analysis.R index 424d49e..b2765c2 100644 --- a/R/roi_connectivity_analysis.R +++ b/R/roi_connectivity_analysis.R @@ -897,12 +897,8 @@ message("Study: ", config$name) message("Groups: ", paste(group_order, collapse = ", ")) message("Bands: ", paste(names(config$bands), collapse = ", ")) -# Load roi_categories from atlas file if provided -if (!is.null(args$roi_categories) && file.exists(args$roi_categories)) { - config$roi_categories <- read_yaml(args$roi_categories) - message("Loaded roi_categories from: ", args$roi_categories, - " (", length(config$roi_categories), " regions)") -} +# Study-config categories win; the atlas file is only a fallback (stats_utils.R). +config$roi_categories <- resolve_roi_categories(config$roi_categories, args$roi_categories) # =========================================================================== # 1. Global connectivity analysis diff --git a/R/roi_cross_freq_edges_analysis.R b/R/roi_cross_freq_edges_analysis.R index d477aee..4cda9c7 100644 --- a/R/roi_cross_freq_edges_analysis.R +++ b/R/roi_cross_freq_edges_analysis.R @@ -171,13 +171,8 @@ group_labels <- unlist(config$groups) group_order <- config$group_order message("Study: ", config$name) -if (!is.null(args$roi_categories) && file.exists(args$roi_categories)) { - rc <- read_yaml(args$roi_categories) - if (length(rc) == 1 && identical(names(rc), "roi_categories")) rc <- rc[["roi_categories"]] - config$roi_categories <- rc - message("Loaded roi_categories from: ", args$roi_categories, - " (", length(config$roi_categories), " regions)") -} +# Study-config categories win; the atlas file is only a fallback (stats_utils.R). +config$roi_categories <- resolve_roi_categories(config$roi_categories, args$roi_categories) metrics <- if (!is.null(args$metric)) trimws(strsplit(args$metric, ",")[[1]]) else EDGE_METRICS metrics <- intersect(EDGE_METRICS, metrics) diff --git a/R/roi_evoked_analysis.R b/R/roi_evoked_analysis.R index 9598bad..eca99ad 100644 --- a/R/roi_evoked_analysis.R +++ b/R/roi_evoked_analysis.R @@ -109,12 +109,8 @@ message("Contrasts: ", length(contrasts)) message("Study: ", config$name) message("Groups: ", paste(group_order, collapse = ", ")) -# Load roi_categories from atlas file if provided -if (!is.null(args$roi_categories) && file.exists(args$roi_categories)) { - config$roi_categories <- read_yaml(args$roi_categories) - message("Loaded roi_categories from: ", args$roi_categories, - " (", length(config$roi_categories), " regions)") -} +# Study-config categories win; the atlas file is only a fallback (stats_utils.R). +config$roi_categories <- resolve_roi_categories(config$roi_categories, args$roi_categories) # Get unique measure names measure_names <- unique(measures_df$measure_name) diff --git a/R/roi_pac_analysis.R b/R/roi_pac_analysis.R index 41fdfa8..f73cc05 100644 --- a/R/roi_pac_analysis.R +++ b/R/roi_pac_analysis.R @@ -882,12 +882,8 @@ message("Study: ", config$name) message("Groups: ", paste(group_order, collapse = ", ")) message("Freq pairs: ", paste(unique(pac$freq_pair), collapse = ", ")) -# Load roi_categories from atlas file if provided -if (!is.null(args$roi_categories) && file.exists(args$roi_categories)) { - config$roi_categories <- read_yaml(args$roi_categories) - message("Loaded roi_categories from: ", args$roi_categories, - " (", length(config$roi_categories), " regions)") -} +# Study-config categories win; the atlas file is only a fallback (stats_utils.R). +config$roi_categories <- resolve_roi_categories(config$roi_categories, args$roi_categories) # =========================================================================== # 1. Global PAC analysis (summary always needed for figures) diff --git a/R/roi_psd_analysis.R b/R/roi_psd_analysis.R index 2d25285..77011be 100644 --- a/R/roi_psd_analysis.R +++ b/R/roi_psd_analysis.R @@ -137,12 +137,8 @@ if ("delta_ref" %in% power_types && !is.null(config$delta_reference$exclude_band message(" delta_ref: excluded reference band(s) from testing: ", paste(excl, collapse = ", ")) } -# Load roi_categories from atlas file if provided -if (!is.null(args$roi_categories) && file.exists(args$roi_categories)) { - config$roi_categories <- read_yaml(args$roi_categories) - message("Loaded roi_categories from: ", args$roi_categories, - " (", length(config$roi_categories), " regions)") -} +# Study-config categories win; the atlas file is only a fallback (stats_utils.R). +config$roi_categories <- resolve_roi_categories(config$roi_categories, args$roi_categories) message("Study: ", config$name) message("Groups: ", paste(group_order, collapse = ", ")) diff --git a/R/roi_transfer_entropy_analysis.R b/R/roi_transfer_entropy_analysis.R index a61d054..7dad4a6 100644 --- a/R/roi_transfer_entropy_analysis.R +++ b/R/roi_transfer_entropy_analysis.R @@ -795,16 +795,8 @@ message("Study: ", config$name) message("Groups: ", paste(group_order, collapse = ", ")) message("Bands: ", paste(names(config$bands), collapse = ", ")) -# Load roi_categories from atlas file if provided. The pipeline passes an -# unwrapped file (categories at top level); the documented proposed file wraps -# them under a single `roi_categories:` key — unwrap that so either form works. -if (!is.null(args$roi_categories) && file.exists(args$roi_categories)) { - rc <- read_yaml(args$roi_categories) - if (length(rc) == 1 && identical(names(rc), "roi_categories")) rc <- rc[["roi_categories"]] - config$roi_categories <- rc - message("Loaded roi_categories from: ", args$roi_categories, - " (", length(config$roi_categories), " regions)") -} +# Study-config categories win; the atlas file is only a fallback (stats_utils.R). +config$roi_categories <- resolve_roi_categories(config$roi_categories, args$roi_categories) # =========================================================================== # 1. Global TE analysis (needed for figures — always compute) diff --git a/R/stats_utils.R b/R/stats_utils.R index 04220a2..7492663 100644 --- a/R/stats_utils.R +++ b/R/stats_utils.R @@ -783,3 +783,33 @@ tost_equivalent <- function(estimate, SE, df, margin, alpha = 0.05) { } NA_real_ } + + +#' Resolve the ROI category map an R entry point should use. +#' +#' The study config written by the Python side carries the EFFECTIVE categories: +#' the study's own map, a profile's narrowing of it, or the atlas's default. It +#' always wins. The --roi-categories file is only a fallback for a config that +#' carries none (an R script run by hand). The file used to win, and because it +#' was looked up per directory -- where allen32, allen26 and allen64 all live -- +#' an allen26 study got allen32's partition: two categories matched no parcel and +#' vanished, and Deep Subcortical was built from 4 of its 8 parcels. +#' +#' @param config_categories named list from config$roi_categories (may be NULL) +#' @param path optional path to a roi_categories YAML +#' @return named list of ROI name vectors (possibly empty) +resolve_roi_categories <- function(config_categories, path = NULL) { + if (length(config_categories) > 0) { + message("Using roi_categories from the study config (", + length(config_categories), " regions)") + return(config_categories) + } + if (!is.null(path) && file.exists(path)) { + rc <- yaml::read_yaml(path) + if (length(rc) == 1 && identical(names(rc), "roi_categories")) rc <- rc[["roi_categories"]] + rc[["deprecated_aliases"]] <- NULL + message("Loaded roi_categories from: ", path, " (", length(rc), " regions)") + return(rc) + } + config_categories +} diff --git a/README.md b/README.md index 9e2232f..53d3150 100644 --- a/README.md +++ b/README.md @@ -104,7 +104,7 @@ the extra to install: | Extra | Pulls in | Needed by | |---|---|---| | `mne` | mne | `roi_evoked`, `vertex_evoked`, `electrode_evoked` (Morlet TFR) | -| `mvpa` | scikit-learn | `vertex_signature`, `electrode_signature` | +| `mvpa` | scikit-learn | `vertex_signature`, `electrode_signature`, `roi_signature` | | `network` | networkx | `roi_graph`, `vertex_graph`, `*_nbs`, `*_network` | | `atlas` | nibabel | atlas readers | | `all` | all of the above + dev tools | a full study | @@ -322,6 +322,21 @@ circos_metrics: [imag_coherence, dwpli, pli, aec, coherence] # gallery circos jobs: -1 # default worker count for --jobs (-1/0 = all but one core) +# ── Atlas (optional) ─────────────────────────────────────────────── +# The parcellation the ROI data were extracted with. Resolved BY NAME to that +# atlas's own files through source-localization's registry.yaml, so atlases that +# share a directory (allen32 / allen26 / allen64 all live in allen/) never borrow +# each other's labels or categories. An unknown name is an error, not a guess. +pipeline: + atlas: allen26 +# An atlas that is not in the registry names its files instead (relative paths +# are taken from atlas_dir, else from the source-localization atlas data): +# atlas_files: +# brain_labels: /path/to/labels.nii.gz +# roi_mapping: /path/to/roi_mapping.json +# roi_categories: /path/to/roi_categories.yaml # optional; the study's own +# # roi_categories always win + # ── Random epoch sampling (global default; per-analysis override below) ── # Code defaults when the block is absent: enabled: false, n_bootstrap: 1. epoch_sampling: @@ -378,6 +393,8 @@ paradigms: | `hypotheses[]` `{name, kind, weights/groups/predictor}` | hypothesis layer | the declarative tests, run by name via `--hypothesis` | | `hypotheses[]` `{label, role}` | figures, gallery | readable labels + grouping tag (no gating) | | `bands` | all spectral/connectivity | frequency bands analysed | +| `pipeline.atlas`, `atlas_files`, `atlas_dir` | atlas I/O, R region tier, mosaics | which parcellation the ROI data use: resolved by name through source-localization's `registry.yaml` to that atlas's own labels / mapping / categories / anatomy; `atlas_files` names the files of an unregistered atlas. The 10× voxel convention is read from each NIfTI header, never inferred from its filename | +| `roi_categories` | region tier (Python + R), mosaics | category → ROI map. The study's map (or a profile's narrowing) always wins over the atlas default, on the Python and R sides alike | | `epoch_sampling` | spectral/connectivity, all levels | random-epoch resampling (`n_bootstrap: 0` = full timeseries). Precedence: global → `vertex.epoch_sampling` → per-analysis block | | `jobs` | `run --jobs` default | worker count when `--jobs` is not given | | `.{include_analyses, include_hypotheses, bands, rois}` | `run --profile` | a narrowed study written to its own tree (see below) | @@ -535,7 +552,8 @@ directed families is tracked, equation-checked, in |---|---|---|---| | `electrode_comparison` *(suppl.)* | elec | `electrode_psd` **and** `roi_psd` (same paradigm) | source-vs-electrode band-power concordance + effect-size validation | | `fcd_comparison` *(suppl.)* | elec | `electrode_connectivity` **and** `vertex_connectivity` — normally in *different* paradigms; sibling paradigm dirs are searched, or set `fcd_comparison.{sensor_dir,source_dir}` | source-vs-sensor FCD comparison (mean + spatial CV) per band × metric | -| `electrode_signature` *(suppl. of `electrode_psd`)* | elec | `electrode_psd` | sensor-level neural signature (decoding on electrode band power) — the source-vs-sensor counterpart of `vertex_signature` | +| `electrode_signature` *(suppl. of `electrode_psd`)* | elec | `electrode_psd` | sensor-level neural signature (decoding on electrode band power) — the sensor counterpart of `roi_signature` / `vertex_signature`, compared when one ran in the same paradigm | +| `roi_signature` | roi | — | ROI-level neural signature (decoding on per-parcel band power) — the source-side counterpart of `electrode_signature`; runs on any ROI output, including Monte Carlo operators. Set `sensor_paradigm:` to compare against an `electrode_signature` run in another paradigm | `ANALYSIS_METADATA` records these as `supplements` (the primary the gallery nests them under) plus `requires` (every upstream module, for run ordering). @@ -700,6 +718,7 @@ $SA resting --analysis electrode_aperiodic $SA resting --analysis electrode_comparison # ↳ after electrode_psd AND roi_psd $SA resting --analysis electrode_connectivity # sensor FC comparator $SA resting --analysis electrode_signature # ↳ after electrode_psd +$SA resting --analysis roi_signature # source side of the decoding comparison # Vertex paradigm — whole-brain $SA vertex --analysis vertex_connectivity # PRIMARY (slow; computes matrices) diff --git a/src/source_analytics/analyses/base.py b/src/source_analytics/analyses/base.py index 5b2c4aa..241da1d 100644 --- a/src/source_analytics/analyses/base.py +++ b/src/source_analytics/analyses/base.py @@ -74,14 +74,17 @@ def __init__(self, config: StudyConfig, output_dir: Path): self.fig_dir.mkdir(parents=True, exist_ok=True) self.tbl_dir.mkdir(parents=True, exist_ok=True) - # Resolve atlas directory for on-the-fly ROI extraction - from ..atlas.atlas_utils import find_atlas_dir + # Resolve the study's atlas to its OWN file set (labels, mapping, + # categories, anatomy). The attribute keeps its old name so every + # ``atlas_dir=self._atlas_dir`` call site passes the spec straight + # through: the atlas functions accept a spec wherever they took a dir. + from ..atlas.atlas_utils import AtlasSpec, resolve_atlas - atlas_name = config.raw.get("pipeline", {}).get("atlas") - atlas_dir_cfg = config.raw.get("atlas_dir") try: - self._atlas_dir: Path | None = find_atlas_dir( - atlas_dir_cfg, atlas_name=atlas_name, + self._atlas_dir: AtlasSpec | None = resolve_atlas( + config.raw.get("atlas_dir"), + atlas_name=config.raw.get("pipeline", {}).get("atlas"), + files=config.raw.get("atlas_files"), ) except FileNotFoundError: self._atlas_dir = None @@ -461,7 +464,7 @@ def _call_r_figures_only(self, r_script_name: str, data_csv: str) -> bool: # over from a prior write when this run has no process step (figures-only). import yaml config_path = data_dir / "study_config.yaml" - config_data = dict(self.config.raw) + config_data = self._r_config_data() sfreq = getattr(self, "_sfreq", None) if sfreq is None and config_path.exists(): try: @@ -507,13 +510,31 @@ def _call_r_figures_only(self, r_script_name: str, data_csv: str) -> bool: return False def _r_roi_categories_flags(self) -> list[str]: - """Return ['--roi-categories', path] if atlas roi_categories.yaml exists.""" - if self._atlas_dir is not None: - cat_path = self._atlas_dir / "roi_categories.yaml" - if cat_path.exists(): - return ["--roi-categories", str(cat_path)] + """``['--roi-categories', path]`` for the resolved atlas's own category file. + + Only a fallback for R: ``_r_config_data`` already hands R the effective + categories, and the R side prefers them (``resolve_roi_categories``). + """ + spec = self._atlas_dir + if spec is not None and spec.roi_categories is not None and spec.roi_categories.exists(): + return ["--roi-categories", str(spec.roi_categories)] return [] + def _r_config_data(self) -> dict: + """The config the R side reads (``study_config.yaml``). + + ``raw`` plus the EFFECTIVE ``roi_categories``: the study's own map, a + profile's narrowing of it, or the atlas default. ``raw`` carries neither of + the last two, and R used to fill the gap from the category file in the + atlas *directory* -- in ``allen/``, allen32's -- so allen26 studies lost two + categories and built Deep Subcortical from 4 of its 8 parcels. + """ + data = dict(self.config.raw) + if self.config.roi_categories: + data["roi_categories"] = { + cat: list(rois) for cat, rois in self.config.roi_categories.items()} + return data + def _r_no_figures_flags(self) -> list[str]: """Return ['--no-figures'] if figure generation is disabled, else [].""" if not self._generate_figures: diff --git a/src/source_analytics/analyses/electrode_analysis.py b/src/source_analytics/analyses/electrode_analysis.py index c60517f..b1e984e 100644 --- a/src/source_analytics/analyses/electrode_analysis.py +++ b/src/source_analytics/analyses/electrode_analysis.py @@ -314,7 +314,7 @@ def summary(self) -> None: config_path = data_dir / "study_config.yaml" import yaml - config_data = dict(self.config.raw) + config_data = self._r_config_data() if self._sfreq is not None: config_data["sfreq"] = self._sfreq with open(config_path, "w") as f: diff --git a/src/source_analytics/analyses/electrode_aperiodic_analysis.py b/src/source_analytics/analyses/electrode_aperiodic_analysis.py index 2f12110..c0e2678 100644 --- a/src/source_analytics/analyses/electrode_aperiodic_analysis.py +++ b/src/source_analytics/analyses/electrode_aperiodic_analysis.py @@ -256,7 +256,7 @@ def summary(self) -> None: config_path = data_dir / "study_config.yaml" import yaml - config_data = dict(self.config.raw) + config_data = self._r_config_data() if self._sfreq is not None: config_data["sfreq"] = self._sfreq with open(config_path, "w") as f: diff --git a/src/source_analytics/analyses/electrode_comparison_analysis.py b/src/source_analytics/analyses/electrode_comparison_analysis.py index ec0bd3d..08ef7de 100644 --- a/src/source_analytics/analyses/electrode_comparison_analysis.py +++ b/src/source_analytics/analyses/electrode_comparison_analysis.py @@ -657,7 +657,7 @@ def summary(self) -> None: config_path = data_dir / "study_config.yaml" import yaml - config_data = dict(self.config.raw) + config_data = self._r_config_data() with open(config_path, "w") as f: yaml.dump(config_data, f, default_flow_style=False) diff --git a/src/source_analytics/analyses/electrode_evoked_analysis.py b/src/source_analytics/analyses/electrode_evoked_analysis.py index aa94c11..bdcbc62 100644 --- a/src/source_analytics/analyses/electrode_evoked_analysis.py +++ b/src/source_analytics/analyses/electrode_evoked_analysis.py @@ -402,7 +402,7 @@ def summary(self) -> None: # Write config YAML for R config_path = data_dir / "study_config.yaml" - config_data = dict(self.config.raw) + config_data = self._r_config_data() if self._sfreq is not None: config_data["sfreq"] = self._sfreq with open(config_path, "w") as f: diff --git a/src/source_analytics/analyses/electrode_signature_analysis.py b/src/source_analytics/analyses/electrode_signature_analysis.py index c96cfbc..c000358 100644 --- a/src/source_analytics/analyses/electrode_signature_analysis.py +++ b/src/source_analytics/analyses/electrode_signature_analysis.py @@ -1,11 +1,14 @@ """Electrode (sensor-level) neural-signature classification analysis. -The sensor-space counterpart of ``vertex_signature``: it runs the SAME -classifiers (LOOCV + permutation testing) on per-electrode band power instead of -per-vertex source power, so the two can be compared directly — does source -localization buy us predictability over the raw sensor montage? When the sibling -``vertex_signature`` results are present, a source-vs-sensor accuracy comparison -(per contrast × band × classifier) is emitted. +The sensor-space counterpart of the source signatures (``roi_signature``, +``vertex_signature``): it runs the SAME classifiers (LOOCV + permutation testing) +on per-electrode band power instead of source power, so the two can be compared +directly — does source localization buy us predictability over the raw sensor +montage? When a source signature ran in the SAME paradigm, a source-vs-sensor +accuracy comparison (per contrast × band × classifier) is emitted. Only that +paradigm is searched: a search of the whole results tree once matched a stale +vertex table from an earlier study version and would have reported it as this +run's source side. Reuses the level-agnostic ``stats/signature.py`` machinery unchanged. """ @@ -40,12 +43,99 @@ logger = logging.getLogger(__name__) +def render_signature_comparison( + source_csv: Path, sensor_csv: Path, tbl_dir: Path, fig_dir: Path, *, + source_module: str, +) -> None: + """Source vs sensor decoding accuracy per contrast × band × classifier. + + Writes ``signature_source_vs_sensor.csv`` (+ ``_by_contrast``) into ``tbl_dir`` + and one scatter per classifier into ``fig_dir``, recording which source module + was compared. Called by the sensor module when a source signature ran in its + paradigm, and by ``roi_signature`` when the sensor signature ran in another. + """ + sensor = pd.read_csv(sensor_csv) + source = pd.read_csv(source_csv) + keys = ["contrast", "band", "classifier"] + if not all(k in source.columns for k in keys): + return + merged = source.merge(sensor, on=keys, suffixes=("_source", "_sensor")) + if merged.empty: + return + merged.insert(0, "source_module", source_module) + merged["accuracy_gain"] = merged["accuracy_source"] - merged["accuracy_sensor"] + for m in ("balanced_accuracy", "auc"): + cols = (f"{m}_source", f"{m}_sensor") + if all(c in merged.columns for c in cols): + merged[f"{m}_gain"] = merged[cols[0]] - merged[cols[1]] + # A cell is only INFORMATIVE about the two modalities if at least one of + # them actually decodes; contrasts where both sit at chance contribute + # noise, and averaging them into a headline gain hides real differences. + if {"p_value_source", "p_value_sensor"} <= set(merged.columns): + merged["either_significant"] = ( + (merged["p_value_source"] < 0.05) | (merged["p_value_sensor"] < 0.05)) + merged.to_csv(tbl_dir / "signature_source_vs_sensor.csv", index=False) + + # Per-contrast breakdown — the global mean is dominated by underpowered + # treated-vs-treated contrasts, so report the split explicitly. + gain_cols = [c for c in merged.columns if c.endswith("_gain")] + by_contrast = merged.groupby("contrast")[gain_cols].mean().round(4) + by_contrast.to_csv(tbl_dir / "signature_source_vs_sensor_by_contrast.csv") + + # One panel per classifier. Plot BALANCED accuracy when available (the + # unequal-n contrasts inflate raw accuracy), and mark the cells where at + # least one modality reached significance. + metric = ("balanced_accuracy" + if "balanced_accuracy_source" in merged.columns else "accuracy") + mlabel = "Balanced accuracy" if metric == "balanced_accuracy" else "Accuracy" + for clf in sorted(merged["classifier"].unique()): + sub = merged[merged["classifier"] == clf] + fig, ax = plt.subplots(figsize=(6, 6)) + sig = sub.get("either_significant") + if sig is not None: + ax.scatter(sub.loc[~sig, f"{metric}_sensor"], sub.loc[~sig, f"{metric}_source"], + c="#BDC3C7", s=45, alpha=0.7, edgecolors="white", + label="neither significant") + ax.scatter(sub.loc[sig, f"{metric}_sensor"], sub.loc[sig, f"{metric}_source"], + c="#8E44AD", s=70, alpha=0.9, edgecolors="white", + label="≥1 significant") + ax.legend(fontsize=8, loc="lower right", frameon=False) + else: + ax.scatter(sub[f"{metric}_sensor"], sub[f"{metric}_source"], + c="#8E44AD", s=60, alpha=0.8, edgecolors="white") + lo, hi = 0.3, 1.0 + ax.plot([lo, hi], [lo, hi], "--", color="grey", linewidth=1) + ax.set_xlim(lo, hi); ax.set_ylim(lo, hi); ax.set_aspect("equal") + ax.set_xlabel(f"Sensor {mlabel.lower()} (electrode)") + ax.set_ylabel(f"Source {mlabel.lower()} ({source_module})") + ax.set_title(f"Source vs Sensor decoding — {classifier_label(clf)}\n" + "(above line = source localization gains predictability)") + fig.tight_layout() + fig.savefig(fig_dir / f"signature_source_vs_sensor_{clf}.png", dpi=150) + plt.close(fig) + logger.info("Rendered source-vs-sensor comparison (%d matched cells)", len(merged)) + + class ElectrodeSignatureAnalysis(BaseAnalysis): """Sensor-level whole-montage neural-signature (classification) analysis.""" name = "electrode_signature" SELECTABLE = {"band": "frequency band"} + # Source-signature modules this sensor signature pairs with, in preference + # order. Looked for in THIS paradigm only (see _find_source_results). + _SOURCE_SIGNATURES = ("roi_signature", "vertex_signature") + + _SUMMARY_TITLE = "Electrode Neural Signature Analysis Summary" + _SUMMARY_ANALYSIS = "Sensor-level (electrode) neural signature (classification)" + _SUMMARY_METHODS = ( + "Each classifier, with LOOCV, was trained to distinguish groups from the " + "spatial pattern of per-electrode relative band power. Significance was " + "assessed by permutation testing. This is the sensor-space counterpart of " + "the source-localized neural signatures (roi_signature, vertex_signature); " + "the `signature_source_vs_sensor` table/figure compare them when one ran in " + "the same paradigm.") + def __init__(self, config: StudyConfig, output_dir: Path): super().__init__(config, output_dir) self._sfreq: float | None = None @@ -60,7 +150,7 @@ def __init__(self, config: StudyConfig, output_dir: Path): # Classifiers: `classifiers:` (list) or `classifier:` (scalar), normalised # and deduped in config order — same contract as vertex_signature. - sig_cfg = config.raw.get("electrode_signature", {}) + sig_cfg = config.raw.get(self.name, {}) raw_clfs = sig_cfg.get("classifiers") or [sig_cfg.get("classifier", "svm_linear")] seen: set[str] = set() self._classifiers: list[str] = [] @@ -197,8 +287,8 @@ def aggregate(self) -> None: if feat_df.empty: logger.warning("No electrode signature feature data collected") return - feat_df.to_csv(data_dir / "electrode_signature_features.csv", index=False) - logger.info("Exported electrode_signature_features.csv (%d rows)", len(feat_df)) + feat_df.to_csv(data_dir / f"{self.name}_features.csv", index=False) + logger.info(f"Exported {self.name}_features.csv (%d rows)", len(feat_df)) # Persist the montage layout (channel, x, y, z) for topomap rendering — # the sensor analog of source_coords.csv. @@ -272,8 +362,8 @@ def statistics(self) -> None: if all_results: pd.DataFrame(all_results).to_csv( - tbl_dir / "electrode_signature_results.csv", index=False) - logger.info("Exported electrode_signature_results.csv (%d rows)", len(all_results)) + tbl_dir / f"{self.name}_results.csv", index=False) + logger.info(f"Exported {self.name}_results.csv (%d rows)", len(all_results)) if self._signature_results: data_dir = self.output_dir / "data" @@ -292,12 +382,12 @@ def statistics(self) -> None: "n_permutations": r.n_permutations, "classifier": r.classifier, "has_weights": r.has_weights, } - with open(data_dir / "electrode_signature_results.pkl", "wb") as f: + with open(data_dir / f"{self.name}_results.pkl", "wb") as f: pickle.dump(pkl_data, f) def _load_state_from_disk(self) -> bool: data_dir = self.output_dir / "data" - pkl_path = data_dir / "electrode_signature_results.pkl" + pkl_path = data_dir / f"{self.name}_results.pkl" if not pkl_path.exists(): logger.warning("No saved electrode signature state at %s; skipping figures", pkl_path) return False @@ -363,7 +453,7 @@ def figures(self) -> None: if getattr(result, "has_weights", True) and not np.all(np.isnan(result.feature_weights)): self._plot_importance_topomap( result.feature_weights, f"Feature Importance — {key}", - fig_dir / f"electrode_signature_importance_{safe}.png") + fig_dir / f"{self.name}_importance_{safe}.png") fig, ax = plt.subplots(figsize=(8, 5)) ax.hist(result.null_distribution, bins=30, color="#3498DB", @@ -373,7 +463,7 @@ def figures(self) -> None: ax.set_xlabel("Accuracy"); ax.set_ylabel("Count") ax.set_title(f"Signature Permutation Test — {key}") ax.legend(); fig.tight_layout() - fig.savefig(fig_dir / f"electrode_signature_null_{safe}.png", dpi=150) + fig.savefig(fig_dir / f"{self.name}_null_{safe}.png", dpi=150) plt.close(fig) fig, ax = plt.subplots(figsize=(5, 4)) @@ -388,106 +478,51 @@ def figures(self) -> None: ax.set_xticklabels(["Pred 0", "Pred 1"]); ax.set_yticklabels(["True 0", "True 1"]) ax.set_title(f"Confusion Matrix — {key}") fig.tight_layout() - fig.savefig(fig_dir / f"electrode_signature_confusion_{safe}.png", dpi=150) + fig.savefig(fig_dir / f"{self.name}_confusion_{safe}.png", dpi=150) plt.close(fig) self._render_source_vs_sensor(fig_dir) - def _find_vertex_results(self) -> Path | None: - """Locate the sibling vertex_signature results table (best-effort).""" - # tbl_dir = .../results/tables//electrode_signature - tables_root = self.tbl_dir.parent.parent - for cand in tables_root.glob("**/vertex_signature_results.csv"): - return cand + def _find_source_results(self) -> tuple[str, Path] | None: + """The source-signature results table in THIS paradigm, if one ran. + + ``tbl_dir`` is ``.../tables//``, so its siblings are the + other modules of the same paradigm. Nothing outside it is searched: the + previous results-tree-wide glob returned the first + ``vertex_signature_results.csv`` anywhere, which was a stale table from an + earlier study version. + """ + for mod in self._SOURCE_SIGNATURES: + cand = self.tbl_dir.parent / mod / f"{mod}_results.csv" + if cand.exists(): + return mod, cand return None def _render_source_vs_sensor(self, fig_dir: Path) -> None: - """Source (vertex) vs sensor (electrode) decoding accuracy per - contrast × band × classifier — the headline comparison.""" - sensor_csv = self.tbl_dir / "electrode_signature_results.csv" - vertex_csv = self._find_vertex_results() - if not sensor_csv.exists() or vertex_csv is None: - logger.info("Source-vs-sensor comparison skipped (need both " - "vertex_signature and electrode_signature results).") - return - sensor = pd.read_csv(sensor_csv) - source = pd.read_csv(vertex_csv) - keys = ["contrast", "band", "classifier"] - if not all(k in source.columns for k in keys): - return - merged = source.merge(sensor, on=keys, suffixes=("_source", "_sensor")) - if merged.empty: + """Compare against the source signature of this paradigm, if one ran.""" + sensor_csv = self.tbl_dir / f"{self.name}_results.csv" + found = self._find_source_results() + if not sensor_csv.exists() or found is None: + logger.info("Source-vs-sensor comparison skipped (no source signature in " + "this paradigm; looked for %s).", ", ".join(self._SOURCE_SIGNATURES)) return - merged["accuracy_gain"] = merged["accuracy_source"] - merged["accuracy_sensor"] - for m in ("balanced_accuracy", "auc"): - cols = (f"{m}_source", f"{m}_sensor") - if all(c in merged.columns for c in cols): - merged[f"{m}_gain"] = merged[cols[0]] - merged[cols[1]] - # A cell is only INFORMATIVE about the two modalities if at least one of - # them actually decodes; contrasts where both sit at chance contribute - # noise, and averaging them into a headline gain hides real differences. - if {"p_value_source", "p_value_sensor"} <= set(merged.columns): - merged["either_significant"] = ( - (merged["p_value_source"] < 0.05) | (merged["p_value_sensor"] < 0.05)) - merged.to_csv(self.tbl_dir / "signature_source_vs_sensor.csv", index=False) - - # Per-contrast breakdown — the global mean is dominated by underpowered - # treated-vs-treated contrasts, so report the split explicitly. - gain_cols = [c for c in merged.columns if c.endswith("_gain")] - by_contrast = merged.groupby("contrast")[gain_cols].mean().round(4) - by_contrast.to_csv(self.tbl_dir / "signature_source_vs_sensor_by_contrast.csv") - - # One panel per classifier. Plot BALANCED accuracy when available (the - # unequal-n contrasts inflate raw accuracy), and mark the cells where at - # least one modality reached significance. - metric = ("balanced_accuracy" - if "balanced_accuracy_source" in merged.columns else "accuracy") - mlabel = "Balanced accuracy" if metric == "balanced_accuracy" else "Accuracy" - for clf in sorted(merged["classifier"].unique()): - sub = merged[merged["classifier"] == clf] - fig, ax = plt.subplots(figsize=(6, 6)) - sig = sub.get("either_significant") - if sig is not None: - ax.scatter(sub.loc[~sig, f"{metric}_sensor"], sub.loc[~sig, f"{metric}_source"], - c="#BDC3C7", s=45, alpha=0.7, edgecolors="white", - label="neither significant") - ax.scatter(sub.loc[sig, f"{metric}_sensor"], sub.loc[sig, f"{metric}_source"], - c="#8E44AD", s=70, alpha=0.9, edgecolors="white", - label="≥1 significant") - ax.legend(fontsize=8, loc="lower right", frameon=False) - else: - ax.scatter(sub[f"{metric}_sensor"], sub[f"{metric}_source"], - c="#8E44AD", s=60, alpha=0.8, edgecolors="white") - lo, hi = 0.3, 1.0 - ax.plot([lo, hi], [lo, hi], "--", color="grey", linewidth=1) - ax.set_xlim(lo, hi); ax.set_ylim(lo, hi); ax.set_aspect("equal") - ax.set_xlabel(f"Sensor {mlabel.lower()} (electrode)") - ax.set_ylabel(f"Source {mlabel.lower()} (vertex)") - ax.set_title(f"Source vs Sensor decoding — {classifier_label(clf)}\n" - "(above line = source localization gains predictability)") - fig.tight_layout() - fig.savefig(fig_dir / f"signature_source_vs_sensor_{clf}.png", dpi=150) - plt.close(fig) - logger.info("Rendered source-vs-sensor comparison (%d matched cells)", len(merged)) + render_signature_comparison(found[1], sensor_csv, self.tbl_dir, fig_dir, + source_module=found[0]) def summary(self) -> None: tbl_dir = self.tbl_dir models = ", ".join(classifier_label(c) for c in self._classifiers) lines = [ - "# Electrode Neural Signature Analysis Summary", "", + f"# {self._SUMMARY_TITLE}", "", f"**Study**: {self.config.name}", - "**Analysis**: Sensor-level (electrode) neural signature (classification)", + f"**Analysis**: {self._SUMMARY_ANALYSIS}", f"**Classifiers**: {models}", f"**CV method**: {self._cv_method}", f"**Permutations**: {self._n_permutations}", "", "## Methods", "", - "Each classifier, with LOOCV, was trained to distinguish groups from the " - "spatial pattern of per-electrode relative band power. Significance was " - "assessed by permutation testing. This is the sensor-space counterpart of " - "the vertex (source-localized) neural signature; the " - "`signature_source_vs_sensor` table/figure compare the two.", "", + self._SUMMARY_METHODS, "", ] - results_csv = tbl_dir / "electrode_signature_results.csv" + results_csv = tbl_dir / f"{self.name}_results.csv" if results_csv.exists(): df = pd.read_csv(results_csv) has_model = "model" in df.columns diff --git a/src/source_analytics/analyses/roi_aperiodic_analysis.py b/src/source_analytics/analyses/roi_aperiodic_analysis.py index 2b26008..8db164e 100644 --- a/src/source_analytics/analyses/roi_aperiodic_analysis.py +++ b/src/source_analytics/analyses/roi_aperiodic_analysis.py @@ -168,7 +168,7 @@ def summary(self) -> None: config_path = data_dir / "study_config.yaml" import yaml - config_data = dict(self.config.raw) + config_data = self._r_config_data() if self._sfreq is not None: config_data["sfreq"] = self._sfreq with open(config_path, "w") as f: @@ -244,5 +244,5 @@ def _render_brain_mosaics(self) -> None: correction_label="FDR", facet_cols=["hypothesis", "dv"], colorbar_label="Hedges' g", - auto_slices=True, + auto_slices=True,atlas=self._atlas_dir ) diff --git a/src/source_analytics/analyses/roi_connectivity_analysis.py b/src/source_analytics/analyses/roi_connectivity_analysis.py index 1e2a9ef..25ed606 100644 --- a/src/source_analytics/analyses/roi_connectivity_analysis.py +++ b/src/source_analytics/analyses/roi_connectivity_analysis.py @@ -319,7 +319,7 @@ def summary(self) -> None: # Write study config YAML for R config_path = data_dir / "study_config.yaml" - config_data = dict(self.config.raw) + config_data = self._r_config_data() if self._sfreq is not None: config_data["sfreq"] = self._sfreq with open(config_path, "w") as f: diff --git a/src/source_analytics/analyses/roi_cross_freq_analysis.py b/src/source_analytics/analyses/roi_cross_freq_analysis.py index 60cb423..4f80b12 100644 --- a/src/source_analytics/analyses/roi_cross_freq_analysis.py +++ b/src/source_analytics/analyses/roi_cross_freq_analysis.py @@ -295,7 +295,7 @@ def _run_r_script(self, script_name: str, *, extra_args: list[str] | None = None return False config_path = data_dir / "study_config.yaml" - config_data = dict(self.config.raw) + config_data = self._r_config_data() if self._sfreq is not None: config_data["sfreq"] = self._sfreq with open(config_path, "w") as f: @@ -355,5 +355,5 @@ def _render_brain_mosaics(self) -> None: analysis_name="roi_cross_freq", effect_col="hedges_g", roi_col="region", facet_cols=["contrast", "freq_pair"], - colorbar_label="Hedges' g", + colorbar_label="Hedges' g",atlas=self._atlas_dir ) diff --git a/src/source_analytics/analyses/roi_directed_analysis.py b/src/source_analytics/analyses/roi_directed_analysis.py index 258095e..cd8ab84 100644 --- a/src/source_analytics/analyses/roi_directed_analysis.py +++ b/src/source_analytics/analyses/roi_directed_analysis.py @@ -210,7 +210,7 @@ def summary(self) -> None: # Write study config YAML for R config_path = data_dir / "study_config.yaml" - config_data = dict(self.config.raw) + config_data = self._r_config_data() if self._sfreq is not None: config_data["sfreq"] = self._sfreq with open(config_path, "w") as f: diff --git a/src/source_analytics/analyses/roi_evoked_analysis.py b/src/source_analytics/analyses/roi_evoked_analysis.py index 938579c..d03414a 100644 --- a/src/source_analytics/analyses/roi_evoked_analysis.py +++ b/src/source_analytics/analyses/roi_evoked_analysis.py @@ -294,7 +294,7 @@ def summary(self) -> None: # Write config YAML for R config_path = data_dir / "study_config.yaml" - config_data = dict(self.config.raw) + config_data = self._r_config_data() if self._sfreq is not None: config_data["sfreq"] = self._sfreq with open(config_path, "w") as f: diff --git a/src/source_analytics/analyses/roi_network_analysis.py b/src/source_analytics/analyses/roi_network_analysis.py index 262736c..8b489ac 100644 --- a/src/source_analytics/analyses/roi_network_analysis.py +++ b/src/source_analytics/analyses/roi_network_analysis.py @@ -499,5 +499,5 @@ def summary(self) -> None: data_dir = self.output_dir / "data" config_path = data_dir / "study_config.yaml" with open(config_path, "w") as f: - yaml.dump(dict(self.config.raw), f, default_flow_style=False) + yaml.dump(self._r_config_data(), f, default_flow_style=False) self._write_summary(graph=True, nbs=True) diff --git a/src/source_analytics/analyses/roi_psd_analysis.py b/src/source_analytics/analyses/roi_psd_analysis.py index c2b016b..9309bca 100644 --- a/src/source_analytics/analyses/roi_psd_analysis.py +++ b/src/source_analytics/analyses/roi_psd_analysis.py @@ -202,7 +202,7 @@ def summary(self) -> None: config_path = data_dir / "study_config.yaml" import yaml # Always write config so sfreq is up-to-date - config_data = dict(self.config.raw) + config_data = self._r_config_data() if self._sfreq is not None: config_data["sfreq"] = self._sfreq with open(config_path, "w") as f: @@ -280,5 +280,5 @@ def _render_brain_mosaics(self) -> None: correction_label="FDR", facet_cols=["hypothesis", "band", "dv"], colorbar_label="Hedges' g", - auto_slices=True, + auto_slices=True,atlas=self._atlas_dir ) diff --git a/src/source_analytics/analyses/roi_signature_analysis.py b/src/source_analytics/analyses/roi_signature_analysis.py new file mode 100644 index 0000000..c903e1b --- /dev/null +++ b/src/source_analytics/analyses/roi_signature_analysis.py @@ -0,0 +1,153 @@ +"""ROI-level neural-signature classification analysis. + +The source-side counterpart of ``electrode_signature``: the SAME classifiers, +cross-validation and permutation testing, on per-parcel relative band power from +the source-localized ROI time series instead of per-electrode band power. The +feature estimator matches ``electrode_signature`` exactly (Welch PSD per series, +relative band power, this module's own ``epoch_sampling``), so where both ran on +the same epochs the accuracies are directly comparable. + +Unlike ``vertex_signature`` it needs no vertex-level estimate, so it runs on any +ROI output, including Monte Carlo operators, which never build one. + +``electrode_signature`` renders the comparison when both ran in the same +paradigm. The sensor signature reads the raw recordings, so one run can serve +every source arm: when it ran in another paradigm, name that paradigm with +``sensor_paradigm:`` and this module renders the comparison itself. +""" + +from __future__ import annotations + +import logging +from pathlib import Path + +import numpy as np +import pandas as pd +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt + +from ..io.discovery import SubjectInfo +from ..io.loader import SubjectLoader +from ..spectral.band_power import extract_band_power, relative_power_kwargs +from ..spectral.psd import compute_psd +from .electrode_signature_analysis import ( + ElectrodeSignatureAnalysis, + render_signature_comparison, +) + +logger = logging.getLogger(__name__) + + +class ROISignatureAnalysis(ElectrodeSignatureAnalysis): + """Per-parcel neural-signature (classification) analysis.""" + + name = "roi_signature" + + _SUMMARY_TITLE = "ROI Neural Signature Analysis Summary" + _SUMMARY_ANALYSIS = "Source-level (ROI) neural signature (classification)" + _SUMMARY_METHODS = ( + "Each classifier, with LOOCV, was trained to distinguish groups from the " + "spatial pattern of per-ROI relative band power in the source-localized " + "time series. Significance was assessed by permutation testing. The feature " + "estimator matches electrode_signature, so the `signature_source_vs_sensor` " + "comparison is like for like.") + + def setup(self) -> None: + # No electrode roster: the features come from the ROI time series. + self._feature_rows.clear() + self._subject_data.clear() + self._subject_groups.clear() + self._subject_order.clear() + self._ch_names = None + self._ch_coords = None + self._signature_results.clear() + + def process_subject(self, subject: SubjectInfo) -> None: + loader = SubjectLoader(subject.data_dir) + roi_ts = loader.load_or_extract_roi_timeseries( + signed=True, atlas_dir=self._atlas_dir, rois=self.config.rois) + sfreq = loader.load_sfreq() + if self._sfreq is None: + self._sfreq = sfreq + rois = list(roi_ts) + if self._ch_names is None: + self._ch_names = rois # feature order; statistics() matches by name + + data = np.stack([np.asarray(roi_ts[r], dtype=float) for r in rois]) + fmax = max(hi for _, hi in self.config.bands.values()) + 10 + rel: dict[tuple[str, str], list[float]] = {} + for draw in self._get_draws(data, sfreq): + for i, roi in enumerate(rois): + x = draw[i, :] + if np.all(x == 0) or np.any(np.isnan(x)): + continue + freqs, psd = compute_psd(x, sfreq, fmax=fmax) + bp = extract_band_power( + freqs, psd, self._selected_bands(), + **relative_power_kwargs(self.config.raw.get("relative_power"))) + for band, vals in bp.items(): + rel.setdefault((roi, band), []).append(vals["relative"]) + + uid = f"{subject.group}_{subject.subject_id}" + band_power: dict[str, dict[str, float]] = {} + for (roi, band), vals in rel.items(): + value = float(np.mean(vals)) + band_power.setdefault(band, {})[roi] = value + self._feature_rows.append({"subject": uid, "group": subject.group, + "roi": roi, "band": band, "relative": value}) + self._subject_groups[uid] = subject.group + self._subject_order.append(uid) + self._subject_data[uid] = {"band_power": band_power} + + def aggregate(self) -> None: + super().aggregate() + if self._ch_names: + # The results pickle stores weights without names; figures-only runs + # need the order back. + pd.DataFrame({"roi": self._ch_names}).to_csv( + self.output_dir / "data" / "roi_order.csv", index=False) + + def _load_state_from_disk(self) -> bool: + ok = super()._load_state_from_disk() + order = self.output_dir / "data" / "roi_order.csv" + if ok and order.exists(): + self._ch_names = pd.read_csv(order)["roi"].tolist() + return ok + + def _plot_importance_topomap(self, values: np.ndarray, title: str, out_path: Path) -> None: + """|weight| per parcel, ranked (parcels have no montage positions).""" + if not self._ch_names: + return + w = np.abs(np.asarray(values, dtype=float)) + ok = ~np.isnan(w) + if not ok.any(): + return + names, vals = np.asarray(self._ch_names)[ok], w[ok] + order = np.argsort(vals) + fig, ax = plt.subplots(figsize=(6, max(3.0, 0.28 * len(vals)))) + ax.barh(names[order], vals[order], color="#C0392B") + ax.set_xlabel("|weight|") + ax.set_title(title) + fig.tight_layout() + fig.savefig(out_path, dpi=150) + plt.close(fig) + + def _render_source_vs_sensor(self, fig_dir: Path) -> None: + """Compare against a sensor signature from ANOTHER paradigm, if named. + + In the same paradigm ``electrode_signature`` renders the comparison; + doing it here too would duplicate it. + """ + para = self.config.raw.get(self.name, {}).get("sensor_paradigm") + if not para or para == self.config.paradigm_name: + return + sensor_csv = (self.tbl_dir.parent.parent / para / "electrode_signature" + / "electrode_signature_results.csv") + source_csv = self.tbl_dir / f"{self.name}_results.csv" + if not (sensor_csv.exists() and source_csv.exists()): + logger.info("Source-vs-sensor comparison skipped: %s not found", sensor_csv) + return + render_signature_comparison(source_csv, sensor_csv, self.tbl_dir, fig_dir, + source_module=self.name) diff --git a/src/source_analytics/analyses/vertex_cluster_analysis.py b/src/source_analytics/analyses/vertex_cluster_analysis.py index 73fd613..88a8bfb 100644 --- a/src/source_analytics/analyses/vertex_cluster_analysis.py +++ b/src/source_analytics/analyses/vertex_cluster_analysis.py @@ -664,7 +664,7 @@ def summary(self) -> None: # Write study config for R config_path = data_dir / "study_config.yaml" - config_data = dict(self.config.raw) + config_data = self._r_config_data() if self._sfreq is not None: config_data["sfreq"] = self._sfreq with open(config_path, "w") as f: diff --git a/src/source_analytics/analyses/vertex_connectivity_analysis.py b/src/source_analytics/analyses/vertex_connectivity_analysis.py index 1e5d8ff..7c63c18 100644 --- a/src/source_analytics/analyses/vertex_connectivity_analysis.py +++ b/src/source_analytics/analyses/vertex_connectivity_analysis.py @@ -414,7 +414,7 @@ def summary(self) -> None: data_dir = self.output_dir / "data" config_path = data_dir / "study_config.yaml" - config_data = dict(self.config.raw) + config_data = self._r_config_data() if self._sfreq is not None: config_data["sfreq"] = self._sfreq with open(config_path, "w") as f: diff --git a/src/source_analytics/analyses/vertex_cross_freq_analysis.py b/src/source_analytics/analyses/vertex_cross_freq_analysis.py index c7c2e84..941c1d9 100644 --- a/src/source_analytics/analyses/vertex_cross_freq_analysis.py +++ b/src/source_analytics/analyses/vertex_cross_freq_analysis.py @@ -329,7 +329,7 @@ def figures(self) -> None: def summary(self) -> None: data_dir = self.output_dir / "data" - cfg = dict(self.config.raw) + cfg = self._r_config_data() if self._sfreq is not None: cfg["sfreq"] = self._sfreq with open(data_dir / "study_config.yaml", "w") as f: diff --git a/src/source_analytics/analyses/vertex_network_analysis.py b/src/source_analytics/analyses/vertex_network_analysis.py index 75180e5..c4977d4 100644 --- a/src/source_analytics/analyses/vertex_network_analysis.py +++ b/src/source_analytics/analyses/vertex_network_analysis.py @@ -179,7 +179,9 @@ def _merge_matrices(self, payload: dict) -> None: len(self._vertex_indices)) try: from ..atlas import find_atlas_dir, load_vertex_roi_labels - atlas_labels = load_vertex_roi_labels(self._source_coords, find_atlas_dir()) + atlas_labels = load_vertex_roi_labels( + self._source_coords, + self._atlas_dir if self._atlas_dir is not None else find_atlas_dir()) except Exception: atlas_labels = None self._vertex_labels = _generate_vertex_labels(self._source_coords, atlas_labels) @@ -577,7 +579,7 @@ def figures(self) -> None: def summary(self) -> None: data_dir = self.output_dir / "data" config_path = data_dir / "study_config.yaml" - config_data = dict(self.config.raw) + config_data = self._r_config_data() if self._sfreq is not None: config_data["sfreq"] = self._sfreq with open(config_path, "w") as f: diff --git a/src/source_analytics/analyses/vertex_signature_analysis.py b/src/source_analytics/analyses/vertex_signature_analysis.py index b93ab24..838bfac 100644 --- a/src/source_analytics/analyses/vertex_signature_analysis.py +++ b/src/source_analytics/analyses/vertex_signature_analysis.py @@ -370,7 +370,7 @@ def summary(self) -> None: data_dir = self.output_dir / "data" config_path = data_dir / "study_config.yaml" - config_data = dict(self.config.raw) + config_data = self._r_config_data() if self._sfreq is not None: config_data["sfreq"] = self._sfreq with open(config_path, "w") as f: diff --git a/src/source_analytics/analyses/vertex_specparam_analysis.py b/src/source_analytics/analyses/vertex_specparam_analysis.py index 9eb030d..4cfc789 100644 --- a/src/source_analytics/analyses/vertex_specparam_analysis.py +++ b/src/source_analytics/analyses/vertex_specparam_analysis.py @@ -700,7 +700,7 @@ def summary(self) -> None: data_dir = self.output_dir / "data" config_path = data_dir / "study_config.yaml" - config_data = dict(self.config.raw) + config_data = self._r_config_data() if self._sfreq is not None: config_data["sfreq"] = self._sfreq with open(config_path, "w") as f: diff --git a/src/source_analytics/atlas/__init__.py b/src/source_analytics/atlas/__init__.py index f1b3562..e34ed39 100644 --- a/src/source_analytics/atlas/__init__.py +++ b/src/source_analytics/atlas/__init__.py @@ -1,7 +1,11 @@ """Atlas integration: map vertex coordinates to anatomical ROI labels.""" from .atlas_utils import ( + AtlasSpec, find_atlas_dir, + header_is_inflated, + registered_atlases, + resolve_atlas, load_atlas, load_roi_mapping, load_roi_categories, @@ -10,7 +14,11 @@ ) __all__ = [ + "AtlasSpec", "find_atlas_dir", + "header_is_inflated", + "registered_atlases", + "resolve_atlas", "load_atlas", "load_roi_mapping", "load_roi_categories", diff --git a/src/source_analytics/atlas/atlas_utils.py b/src/source_analytics/atlas/atlas_utils.py index c6221f2..166e426 100644 --- a/src/source_analytics/atlas/atlas_utils.py +++ b/src/source_analytics/atlas/atlas_utils.py @@ -1,9 +1,19 @@ """Atlas utilities: load atlas NIfTI, map vertex coordinates to ROI labels. -Replicates the 10x voxel-size correction from source_localization/utils/atlas.py. -The Atlas_3DRoisLeftRight.Labels.nii header has voxel sizes 10x larger than reality, -so we apply ATLAS_VOXEL_SCALE_FACTOR = 0.1 to both the rotation/scaling block -and the translation vector of the affine. +Some bundled NIfTI files store voxel sizes 10x larger than reality and some +store true units. Which is which is read from the header, as +source_localization.utils.atlas.header_is_inflated does, never guessed from a +filename: the guess was wrong for Atlas_3DRoisLeftRight.Labels.nii, which stores +true units and was being shrunk 10x. Inflated files get +ATLAS_VOXEL_SCALE_FACTOR = 0.1 applied to the rotation/scaling block and the +translation vector of the affine. + +An atlas is resolved by NAME to its own file set (:func:`resolve_atlas`), from +source-localization's ``registry.yaml``. The registry is where an atlas is +defined; a directory is not, because several atlases share one (allen32, allen26 +and allen64 all live in ``allen/``), so "the roi_mapping.json in the atlas +directory" names allen32's files whatever the study asked for. Atlases that are +not registered can name their files explicitly. Also provides on-the-fly ROI extraction from vertex-level source time courses, ported from source_localization/steps/roi_extraction.py. @@ -13,6 +23,8 @@ import json import logging +from dataclasses import dataclass, replace +from functools import lru_cache from pathlib import Path import yaml @@ -103,18 +115,17 @@ def _find_atlas_nifti(atlas_dir: Path) -> tuple[Path, bool]: nii_path : Path Path to the NIfTI file. needs_10x_correction : bool - True for the Antwerp atlas whose header has 10x inflated voxel sizes. - False for the Allen atlas and other atlases with correct headers. + Read from the header (:func:`header_is_inflated`), not the filename. """ # Try Allen atlas first (correct header) allen_path = atlas_dir / _ALLEN_NIFTI if allen_path.exists(): - return allen_path, False + return allen_path, header_is_inflated(allen_path) # Try Antwerp atlas (needs 10x correction) antwerp_path = atlas_dir / _ATLAS_NIFTI if antwerp_path.exists(): - return antwerp_path, True + return antwerp_path, header_is_inflated(antwerp_path) raise FileNotFoundError( f"No atlas NIfTI found in {atlas_dir}. " @@ -122,6 +133,201 @@ def _find_atlas_nifti(atlas_dir: Path) -> tuple[Path, bool]: ) + +# Source-localization's threshold for the inflated-header convention: true-unit +# mouse atlases have sub-0.3 mm voxels, inflated ones 0.8-2 mm. +_INFLATED_ZOOM_THRESHOLD_MM = 0.5 +_REGISTRY_FILE = "registry.yaml" +# Names source-localization still accepts for atlases that were renamed. +_LEGACY_ATLAS_ALIASES = {"full": "antwerp", "coarse_22roi": "coarse22"} +# Atlases the fixed per-directory file names genuinely describe, for when no +# registry is reachable. Anything else must be registered or named explicitly. +_LEGACY_ATLAS_DIRS = {"allen": "allen", "allen32": "allen", "antwerp": "."} +_ATLAS_FILE_KEYS = ("labels", "roi_mapping", "roi_categories", "brain_volume", "brain_mask") + + +def header_is_inflated(nii_path: str | Path) -> bool: + """True if this NIfTI stores the 10x-inflated voxel sizes. + + Mirrors ``source_localization.utils.atlas.header_is_inflated``. Applying the + correction to a true-unit file is silent -- the affine stays well formed and + every coordinate lands 10x away -- so the convention is read, not assumed. + """ + import nibabel as nib + + zooms = nib.load(str(nii_path)).header.get_zooms()[:3] + return bool(max(float(z) for z in zooms) >= _INFLATED_ZOOM_THRESHOLD_MM) + + +@dataclass(frozen=True) +class AtlasSpec: + """One atlas, as the set of files that define it. + + Everything that needs atlas data accepts one of these wherever it used to + accept an atlas directory. ``roi_categories``, ``brain_volume`` and + ``brain_mask`` are optional: allen64 and coarse22 ship no category file. + """ + + name: str | None + labels: Path + roi_mapping: Path + roi_categories: Path | None = None + brain_volume: Path | None = None + brain_mask: Path | None = None + + @property + def labels_inflated(self) -> bool: + return header_is_inflated(self.labels) + + +@lru_cache(maxsize=None) +def _read_registry(path: str) -> dict: + with open(path) as f: + return yaml.safe_load(f) or {} + + +def _registry_path(base_dir: Path) -> Path | None: + p = Path(base_dir) / _REGISTRY_FILE + return p if p.is_file() else None + + +def _registry_file(entry_path: str, registry: Path) -> Path: + """Registry paths are relative to the source-localization package root + (``data/atlas/...``); a hand-written registry may use paths relative to itself.""" + p = Path(entry_path) + if p.is_absolute(): + return p + pkg_root = registry.parent.parent.parent + for root in (registry.parent, pkg_root): + if (root / p).exists(): + return root / p + return pkg_root / p + + +def registered_atlases(base_dir: str | Path | None = None) -> list[str]: + """Names in the atlas registry (empty if none is reachable).""" + reg = _registry_path(Path(base_dir) if base_dir is not None else find_atlas_dir()) + return sorted(_read_registry(str(reg))) if reg else [] + + +def _normalise_files(files: dict | None, base: Path | None) -> dict: + out: dict[str, Path] = {} + for key, value in (files or {}).items(): + k = "labels" if key == "brain_labels" else key + if k not in _ATLAS_FILE_KEYS: + raise ValueError( + f"atlas_files: unknown key '{key}'. Known: brain_labels (or labels), " + "roi_mapping, roi_categories, brain_volume, brain_mask.") + if value is None: + continue + p = Path(value) + out[k] = p if (p.is_absolute() or base is None) else base / p + return out + + +def _spec_from_dir(atlas_dir: Path, name: str | None) -> AtlasSpec: + """The legacy convention: fixed file names inside one directory.""" + labels, _ = _find_atlas_nifti(atlas_dir) + cats = atlas_dir / _ROI_CATEGORIES_FILE + + def _near(fname: str) -> Path | None: + for d in (atlas_dir, atlas_dir.parent): + if (d / fname).exists(): + return d / fname + return None + + return AtlasSpec( + name=name, labels=labels, roi_mapping=atlas_dir / _ROI_MAPPING_FILE, + roi_categories=cats if cats.exists() else None, + brain_volume=_near("Atlas_3DRois.nii"), + brain_mask=_near("Atlas_3DRois_brain.nii.gz"), + ) + + +def resolve_atlas( + atlas: "AtlasSpec | str | Path | None" = None, + *, + atlas_name: str | None = None, + files: dict | None = None, +) -> AtlasSpec: + """Resolve an atlas to its own file set. + + Parameters + ---------- + atlas : AtlasSpec, directory, or None + A spec is returned as is (with any ``files`` overrides applied). A + directory is the atlas base: its ``registry.yaml`` is used when present, + otherwise the legacy fixed file names inside it. None means the + source-localization atlas data. + atlas_name : str, optional + A registered atlas (renamed aliases accepted). + files : dict, optional + Explicit files, for an atlas that is not registered or to override one + that is: ``brain_labels`` (or ``labels``), ``roi_mapping``, + ``roi_categories``, ``brain_volume``, ``brain_mask``. Relative paths are + taken relative to the atlas base directory. + + Raises + ------ + ValueError + ``atlas_name`` cannot be resolved to its own files. Falling back to + whatever files share a directory is how allen32's partition was applied + to allen26 data without a word. + FileNotFoundError + The atlas files do not exist. + """ + if isinstance(atlas, AtlasSpec): + overrides = _normalise_files(files, None) + return replace(atlas, **overrides) if overrides else atlas + + base = Path(atlas) if atlas is not None else find_atlas_dir() + if not base.is_dir(): + raise FileNotFoundError(f"Atlas directory not found: {base}") + overrides = _normalise_files(files, base) + name = _LEGACY_ATLAS_ALIASES.get(atlas_name, atlas_name) if atlas_name else None + explicit = "labels" in overrides and "roi_mapping" in overrides + + spec: AtlasSpec | None = None + registry = _registry_path(base) + if name is not None and registry is not None: + entries = _read_registry(str(registry)) + if name in entries: + ins = entries[name].get("inputs", {}) + opt = {k: _registry_file(ins[k], registry) + for k in ("roi_categories", "brain_volume", "brain_mask") if ins.get(k)} + spec = AtlasSpec( + name=atlas_name, + labels=_registry_file(ins["brain_labels"], registry), + roi_mapping=_registry_file(ins["roi_mapping"], registry), + **opt) + elif not explicit: + raise ValueError( + f"Atlas '{atlas_name}' is not in {registry}. Registered: " + f"{', '.join(sorted(entries))}. Register it there, or name its files " + "under atlas_files: (brain_labels, roi_mapping, roi_categories).") + if spec is None and explicit: + spec = AtlasSpec(name=atlas_name, labels=overrides["labels"], + roi_mapping=overrides["roi_mapping"]) + if spec is None: + if name is None: + spec = _spec_from_dir(base, None) + elif (base / name).is_dir(): + spec = _spec_from_dir(base / name, atlas_name) + elif name in _LEGACY_ATLAS_DIRS and (base / _LEGACY_ATLAS_DIRS[name]).is_dir(): + spec = _spec_from_dir(base / _LEGACY_ATLAS_DIRS[name], atlas_name) + else: + raise ValueError( + f"Cannot resolve atlas '{atlas_name}' under {base}: no registry.yaml " + "there, and its files are not identifiable by name. Name them under " + "atlas_files: (brain_labels, roi_mapping, roi_categories).") + if overrides: + spec = replace(spec, **overrides) + for f in (spec.labels, spec.roi_mapping): + if not Path(f).exists(): + raise FileNotFoundError(f"Atlas file not found: {f}") + return spec + + def load_atlas( atlas_dir: str | Path, *, @@ -131,14 +337,14 @@ def load_atlas( Parameters ---------- - atlas_dir : str or Path - Directory containing the atlas NIfTI file. + atlas_dir : AtlasSpec, str or Path + An atlas: anything :func:`resolve_atlas` accepts. raw_affine : bool, default False If True, return the raw NIfTI affine without any correction. Use this when source coordinates are in the same uncorrected frame as the atlas (e.g., coordinates from source-localization pipeline). - If False (default), apply the 10x voxel correction for the Antwerp - atlas (no-op for Allen atlas which has correct headers). + If False (default), apply the 10x voxel correction when the labels + header is inflated (read from the header; see :func:`header_is_inflated`). Returns ------- @@ -149,8 +355,9 @@ def load_atlas( """ import nibabel as nib - atlas_dir = Path(atlas_dir) - nii_path, needs_correction = _find_atlas_nifti(atlas_dir) + spec = resolve_atlas(atlas_dir) + nii_path = spec.labels + needs_correction = spec.labels_inflated nii = nib.load(str(nii_path)) label_data = np.asarray(nii.dataobj, dtype=np.int32) @@ -158,7 +365,7 @@ def load_atlas( if raw_affine or not needs_correction: return label_data, nii.affine.copy() - # Apply 10x voxel correction (Antwerp atlas only) + # Apply the 10x voxel correction (inflated header) true_affine = nii.affine.copy() true_affine[:3, :3] *= ATLAS_VOXEL_SCALE_FACTOR true_affine[:3, 3] *= ATLAS_VOXEL_SCALE_FACTOR @@ -179,8 +386,7 @@ def load_roi_mapping(atlas_dir: str | Path) -> dict: dict ROI mapping: label_id (str) -> {abbreviation, name, category, color, ...}. """ - atlas_dir = Path(atlas_dir) - mapping_path = atlas_dir / _ROI_MAPPING_FILE + mapping_path = Path(resolve_atlas(atlas_dir).roi_mapping) if not mapping_path.exists(): raise FileNotFoundError(f"ROI mapping not found: {mapping_path}") @@ -191,8 +397,8 @@ def load_roi_mapping(atlas_dir: str | Path) -> dict: def load_roi_categories(atlas_dir: str | Path) -> dict[str, list[str]]: """Load canonical ROI categories from the atlas directory. - Reads ``roi_categories.yaml`` from *atlas_dir*. Returns an empty dict if - the file does not exist (studies can define their own via config). + Reads the atlas's OWN category file (``AtlasSpec.roi_categories``). Returns + an empty dict if the atlas ships none (studies can define their own via config). Parameters ---------- @@ -204,9 +410,8 @@ def load_roi_categories(atlas_dir: str | Path) -> dict[str, list[str]]: dict[str, list[str]] Mapping of category name -> list of ROI names. """ - atlas_dir = Path(atlas_dir) - categories_path = atlas_dir / _ROI_CATEGORIES_FILE - if not categories_path.exists(): + categories_path = resolve_atlas(atlas_dir).roi_categories + if categories_path is None or not Path(categories_path).exists(): return {} with open(categories_path) as f: raw = yaml.safe_load(f) or {} @@ -430,7 +635,7 @@ def label_vertices_to_rois( list[str | None] ROI name per vertex (``None`` only if a vertex mapped to no labeled ROI). """ - atlas_dir = Path(atlas_dir) + atlas_dir = resolve_atlas(atlas_dir) label_data, affine = load_atlas(atlas_dir, raw_affine=True) roi_mapping = load_roi_mapping(atlas_dir) rois = roi_mapping.get("rois", roi_mapping) @@ -526,7 +731,7 @@ def extract_roi_timeseries( sources). Only ROIs with at least one assigned source are included. """ - atlas_dir = Path(atlas_dir) + atlas_dir = resolve_atlas(atlas_dir) # Load atlas with RAW affine — source coordinates from the pipeline # are in the same uncorrected frame as the original NIfTI headers. diff --git a/src/source_analytics/config.py b/src/source_analytics/config.py index 0692123..c368bdd 100644 --- a/src/source_analytics/config.py +++ b/src/source_analytics/config.py @@ -2,6 +2,8 @@ from __future__ import annotations +import logging + import warnings as _warnings from dataclasses import dataclass, field, replace from pathlib import Path @@ -11,16 +13,27 @@ import yaml -def _load_atlas_roi_categories(atlas_name: str | None) -> dict[str, list[str]]: - """Load canonical roi_categories from atlas package data, if available.""" - if not atlas_name: +def _load_atlas_roi_categories( + atlas_name: str | None, + atlas_dir: str | None = None, + atlas_files: dict | None = None, +) -> dict[str, list[str]]: + """The atlas's own default categories, for a study that declares none. + + Resolved by atlas NAME (``resolve_atlas``), so allen26 gets allen26's + partition, not whichever category file shares its directory. + """ + if not (atlas_name or (atlas_files or {}).get("roi_categories")): return {} try: - from source_analytics.atlas import find_atlas_dir, load_roi_categories + from source_analytics.atlas import load_roi_categories, resolve_atlas - atlas_dir = find_atlas_dir(atlas_name=atlas_name) - return load_roi_categories(atlas_dir) - except Exception: + return load_roi_categories( + resolve_atlas(atlas_dir, atlas_name=atlas_name, files=atlas_files)) + except FileNotFoundError: + return {} + except ValueError as exc: + logging.getLogger(__name__).warning("No atlas roi_categories: %s", exc) return {} @@ -420,7 +433,9 @@ def _resolve(p: str) -> Path: paradigms[pname] = pcopy roi_categories = data.get("roi_categories") or _load_atlas_roi_categories( - data.get("pipeline", {}).get("atlas") + data.get("pipeline", {}).get("atlas"), + data.get("atlas_dir"), + data.get("atlas_files"), ) return cls( diff --git a/src/source_analytics/core.py b/src/source_analytics/core.py index 4d2b0cb..dfd9101 100644 --- a/src/source_analytics/core.py +++ b/src/source_analytics/core.py @@ -23,6 +23,7 @@ from .analyses.vertex_specparam_analysis import VertexSpecparamAnalysis from .analyses.vertex_signature_analysis import VertexSignatureAnalysis from .analyses.electrode_signature_analysis import ElectrodeSignatureAnalysis +from .analyses.roi_signature_analysis import ROISignatureAnalysis from .analyses.roi_network_analysis import ( ROINetworkAnalysis, ROIGraphAnalysis, @@ -60,6 +61,7 @@ "vertex_specparam": VertexSpecparamAnalysis, "vertex_signature": VertexSignatureAnalysis, "electrode_signature": ElectrodeSignatureAnalysis, + "roi_signature": ROISignatureAnalysis, "roi_graph": ROIGraphAnalysis, "roi_nbs": ROINBSAnalysis, "vertex_graph": VertexGraphAnalysis, @@ -148,7 +150,8 @@ def resolve_analysis_name(name: str) -> str: "about": "Directional (who-drives-whom) connectivity between ROIs, in two flavors. Transfer entropy (TE, Schreiber 2000) is a model-free information-theoretic measure -- a binned lag-1 estimator of how much one ROI's past reduces uncertainty about another's future -- computed for every directed ROI pair (callosal tracts excluded); its net asymmetry (te - te-transpose) gives the dominant direction. The optional Directed Transfer Function (DTF, Kaminski & Blinowska 1991) derives directed influence from a multivariate autoregressive (MVAR) model fit with ridge regularization (order 8) -- a deviation from ordinary-least-squares DTF used to stabilize the fit against collinear channels. Groups are compared on TE three ways: a global per-pair Welch t-test (Hedges' g), a within-group one-sample test of net TE against zero (is the driving direction consistent within a group), and a region-pair linear mixed model (te ~ group * region_pair + (1|subject)). Read it as: which ROI pairs show a group difference in directed influence, and which region drives which (an up arrow means the first-listed group is higher). DTF, when selected, currently emits directed edges without the R group stats."}, "vertex_directed": {"category": "resting", "level": "vertex", "domain": "Directed", "description": "Vertex DTF outflow/inflow/netflow (ridge-MVAR, cluster-corrected)", "about": "Whole-brain directed connectivity on the dorsal source surface via the Directed Transfer Function (DTF, Kaminski & Blinowska 1991) from a multivariate autoregressive model. Because source vertices are strongly collinear (mean inter-vertex |r| ~ 0.64), the MVAR is fit with ridge regularization (order 8) rather than ordinary least squares, and the fit's stability (spectral radius) is checked. The full all-to-all directed DTF matrix is reduced to three per-vertex maps: outflow (mean directed influence a vertex sends to all others), inflow (mean it receives), and netflow (outflow minus inflow -- net source vs sink). Group differences in each map are tested with a cluster-based permutation test (per-vertex t-statistics clustered by spatial adjacency, cluster-extent FWE from a permutation null; Maris & Oostenveld 2007), with per-vertex Hedges' g. Read it as: spatially-contiguous clusters where the groups differ in how strongly a region drives (outflow), is driven by (inflow), or net-drives (netflow) the rest of the brain -- a cluster with p_corrected < 0.05 marks a region of difference, the sign of its t-values gives direction."}, - "electrode_signature": {"category": "resting", "level": "electrode", "domain": "Source vs Sensor", "display_name": "Neural signature", "description": "Sensor-level neural signature (classification/decoding on electrode band power) — the source-vs-sensor counterpart of vertex_signature"}, + "electrode_signature": {"category": "resting", "level": "electrode", "domain": "Source vs Sensor", "display_name": "Neural signature", "description": "Sensor-level neural signature (classification/decoding on electrode band power) — the sensor counterpart of roi_signature / vertex_signature"}, + "roi_signature": {"category": "resting", "level": "roi", "domain": "Source vs Sensor", "display_name": "Neural signature", "description": "ROI-level neural signature (classification/decoding on per-parcel band power) — the source-side counterpart of electrode_signature; runs on any ROI output"}, "vertex_signature": {"category": "resting", "level": "vertex", "domain": "Multivariate", "display_name": "Neural signature", "description": "Multivariate/ML neural signature (classification, decoding; PCA-reduced with back-projection)"}, "vertex_cluster": {"category": "resting", "level": "vertex", "domain": "Spectral", "description": "Vertex-level cluster permutation", "about": "Whole-brain resting spectral maps on the dorsal source surface: per vertex it computes band power (absolute as mean density in dB/Hz, the same definition as roi_psd, and relative), the 1/f spectral slope, and the peak alpha frequency, then tests where the groups differ. Inference is a cluster-based permutation test -- per-vertex t-statistics are threshold-clustered over neighbouring vertices and each cluster's extent is compared to a permutation null, giving family-wise (FWE) control (Maris & Oostenveld 2007); a threshold-free TFCE variant (Smith & Nichols 2009) is available. Effect sizes are per-vertex Hedges' g. Read it as: spatially-contiguous clusters where the groups differ in a spectral measure -- a cluster with p_corrected < 0.05 marks a region of difference, and the sign of its t-values gives the direction. This is the whole-brain, unparcellated counterpart to the ROI spectral analyses."}, diff --git a/src/source_analytics/io/loader.py b/src/source_analytics/io/loader.py index 1e6d7fa..f336d2e 100644 --- a/src/source_analytics/io/loader.py +++ b/src/source_analytics/io/loader.py @@ -241,7 +241,7 @@ def load_or_extract_roi_timeseries( If True, use signed time courses; otherwise magnitude. atlas_dir : str or Path, optional Atlas directory for on-the-fly extraction. If None, auto-detected - via :func:`~source_analytics.atlas.find_atlas_dir`. + via :func:`~source_analytics.atlas.resolve_atlas`. method : str Source-to-ROI assignment: ``"nearest"`` or ``"proximity"``. proximity_radius_mm : float @@ -273,13 +273,13 @@ def load_or_extract_roi_timeseries( "No step6 ROI files in %s — extracting on-the-fly from step5", self.data_dir, ) - from ..atlas.atlas_utils import extract_roi_timeseries, find_atlas_dir + from ..atlas.atlas_utils import extract_roi_timeseries, resolve_atlas stc_data = self.load_source_timecourses(magnitude=not signed) coords = self.load_source_coords() if atlas_dir is None: - atlas_dir = find_atlas_dir() + atlas_dir = resolve_atlas() result = extract_roi_timeseries( stc_data, diff --git a/src/source_analytics/stats/signature.py b/src/source_analytics/stats/signature.py index 999ada6..f7477bf 100644 --- a/src/source_analytics/stats/signature.py +++ b/src/source_analytics/stats/signature.py @@ -106,6 +106,31 @@ class SignatureResult: balanced_accuracy_ci: tuple[float, float] = (float("nan"), float("nan")) +def _single_threaded(fn): + """Run ``fn`` with the BLAS/OpenMP thread pools limited to one thread. + + A signature run fits a tiny model (tens of subjects by tens of features) + LOOCV x (1 + n_permutations) times. At that size a multithreaded BLAS spends + its time starting and synchronising threads rather than computing: on the + FORGE treatment electrode features the logistic fit took 105-128 s per 21 + LOOCV passes with default threads and 1.4-1.7 s single-threaded, with the + same accuracy. One thread per process is also what the subject pool wants. + """ + import functools + + @functools.wraps(fn) + def wrapper(*args, **kwargs): + try: + from threadpoolctl import threadpool_limits + except ImportError: # installed with scikit-learn; absent only without it + return fn(*args, **kwargs) + with threadpool_limits(limits=1): + return fn(*args, **kwargs) + + return wrapper + + +@_single_threaded def run_signature( features: np.ndarray, labels: np.ndarray, diff --git a/src/source_analytics/viz/brain_roi.py b/src/source_analytics/viz/brain_roi.py index b24e06a..d3b64d2 100644 --- a/src/source_analytics/viz/brain_roi.py +++ b/src/source_analytics/viz/brain_roi.py @@ -25,7 +25,7 @@ import pandas as pd from matplotlib.colors import TwoSlopeNorm -from ..atlas import find_atlas_dir, load_atlas, load_roi_mapping +from ..atlas import find_atlas_dir, load_atlas, load_roi_mapping, resolve_atlas from .palettes import get_diverging_cmap_name logger = logging.getLogger(__name__) @@ -299,12 +299,12 @@ def plot_brain_roi_mosaic( sagittal_slices = DEFAULT_SAGITTAL_SLICES # Load atlas volumes - atlas_dir_path = find_atlas_dir(atlas_dir) - label_data, affine = load_atlas(atlas_dir_path) - roi_mapping = load_roi_mapping(atlas_dir_path) + spec = resolve_atlas(atlas_dir) + label_data, affine = load_atlas(spec) + roi_mapping = load_roi_mapping(spec) - # Anatomy background (skull-stripped) - anat_path = Path(atlas_dir_path) / "Atlas_3DRois_brain.nii.gz" + # Anatomy background (skull-stripped): the atlas's own brain-mask volume + anat_path = spec.brain_mask anat_data = nib.load(str(anat_path)).get_fdata() anat_norm = (anat_data / anat_data.max()) ** gamma @@ -548,15 +548,21 @@ def _pick_informative_slices( def _load_atlas_and_anat( - atlas_name: str = "allen", + atlas="allen", gamma: float = 0.5, ): - """Load atlas label volume, affine, ROI mapping, and anatomy background.""" - atlas_dir = find_atlas_dir(atlas_name=atlas_name) - label_data, affine = load_atlas(atlas_dir) - roi_mapping = load_roi_mapping(atlas_dir) + """Load atlas label volume, affine, ROI mapping, and anatomy background. + + ``atlas`` is an AtlasSpec, an atlas directory, or a registered atlas NAME. + It must be the data's own atlas: drawing allen32's volume for allen26 data + leaves the six merged parcels with no voxels, so they render blank. + """ + spec = (resolve_atlas(atlas_name=atlas) if isinstance(atlas, str) + else resolve_atlas(atlas)) + label_data, affine = load_atlas(spec) + roi_mapping = load_roi_mapping(spec) - anat_path = atlas_dir.parent / "Atlas_3DRois_brain.nii.gz" + anat_path = spec.brain_mask anat_data = nib.load(str(anat_path)).get_fdata() anat_norm = (anat_data / anat_data.max()) ** gamma @@ -715,6 +721,7 @@ def plot_significance_mosaic( alpha: float = 0.05, cmap_name: str = "YlOrRd", atlas_name: str = "allen", + atlas=None, roi_opacity: float = 0.85, dpi: int = 300, ) -> Path: @@ -742,7 +749,7 @@ def plot_significance_mosaic( atlas_name : str Atlas to load (default ``"allen"``). """ - label_data, affine, roi_mapping, anat_norm = _load_atlas_and_anat(atlas_name) + label_data, affine, roi_mapping, anat_norm = _load_atlas_and_anat(atlas if atlas is not None else atlas_name) df = df.copy() df["_q_fdr"] = fdr_bh(df[p_col].values) @@ -791,6 +798,7 @@ def plot_effect_size_mosaic( cmap_name: str = "RdBu_r", colorbar_label: str | None = None, atlas_name: str = "allen", + atlas=None, roi_opacity: float = 0.85, auto_slices: bool = True, dpi: int = 300, @@ -845,7 +853,7 @@ def plot_effect_size_mosaic( a deep thalamic locus missing every panel). Set False to restore the legacy fixed planes. """ - label_data, affine, roi_mapping, anat_norm = _load_atlas_and_anat(atlas_name) + label_data, affine, roi_mapping, anat_norm = _load_atlas_and_anat(atlas if atlas is not None else atlas_name) df = df.copy() if q_col is not None: @@ -937,6 +945,7 @@ def render_posthoc_mosaics( colorbar_label: str = "Hedges' g", alpha: float = 0.05, auto_slices: bool = True, + atlas=None, ) -> list[Path]: """Render effect-size brain mosaics from a posthoc CSV. @@ -977,6 +986,11 @@ def render_posthoc_mosaics( survivors) so no significant ROI is clipped out (MS1 revision). When False, fixed default planes are used. + atlas : AtlasSpec, atlas directory or registered name, optional + The atlas the data were extracted with. None keeps the old default + (``atlas_name="allen"``), which draws allen32's parcels whatever atlas + the data used. + Returns ------- list[Path] @@ -1051,7 +1065,7 @@ def render_posthoc_mosaics( alpha=alpha, cmap_name=cmap_name, colorbar_label=colorbar_label, - auto_slices=auto_slices, + auto_slices=auto_slices,atlas=atlas ) saved.append(out_path) except Exception as exc: @@ -1149,7 +1163,7 @@ def plot_brain_roi( if views is None: views = VIEWS - atlas_dir_resolved = find_atlas_dir(atlas_dir) + atlas_dir_resolved = resolve_atlas(atlas_dir) label_data, affine = load_atlas(atlas_dir_resolved) roi_mapping = load_roi_mapping(atlas_dir_resolved) region_to_labels = _region_to_label_ids(roi_categories, roi_mapping) diff --git a/src/source_analytics/viz/glass_brain.py b/src/source_analytics/viz/glass_brain.py index e7501d6..7a18ba1 100644 --- a/src/source_analytics/viz/glass_brain.py +++ b/src/source_analytics/viz/glass_brain.py @@ -537,7 +537,7 @@ def _load_atlas_slices( Parameters ---------- atlas_dir : Path - Directory containing Atlas_3DRois.nii. + An atlas (anything ``resolve_atlas`` accepts); its ``brain_volume`` is drawn. slice_coords : dict, optional Override slice positions: {"axial_z": 1.5, "coronal_y": 0.0, "sagittal_x": 0.0}. Defaults to midline/dorsal center. @@ -550,16 +550,21 @@ def _load_atlas_slices( """ import nibabel as nib - atlas_dir = Path(atlas_dir) - nii_path = atlas_dir / _ATLAS_INTENSITY_NIFTI - if not nii_path.exists(): - raise FileNotFoundError(f"Atlas intensity template not found: {nii_path}") + from ..atlas.atlas_utils import header_is_inflated, resolve_atlas + + spec = resolve_atlas(atlas_dir) + nii_path = spec.brain_volume + if nii_path is None or not Path(nii_path).exists(): + raise FileNotFoundError( + f"No intensity template (brain_volume) for atlas {spec.name or spec.labels}") nii = nib.load(str(nii_path)) vol = nii.get_fdata() affine = nii.affine.copy() - affine[:3, :3] *= _ATLAS_VOXEL_SCALE - affine[:3, 3] *= _ATLAS_VOXEL_SCALE + # Scale only a header that is actually inflated (read, not assumed). + if header_is_inflated(nii_path): + affine[:3, :3] *= _ATLAS_VOXEL_SCALE + affine[:3, 3] *= _ATLAS_VOXEL_SCALE nx, ny, nz = vol.shape x_coords = affine[0, 3] + np.arange(nx) * affine[0, 0] diff --git a/tests/test_atlas_registry.py b/tests/test_atlas_registry.py new file mode 100644 index 0000000..a6f4aba --- /dev/null +++ b/tests/test_atlas_registry.py @@ -0,0 +1,132 @@ +"""An atlas is resolved by NAME to its own files, never by directory convention. + +Locks out the bug that gave allen26 studies allen32's partition: allen32, allen26 +and allen64 all live in ``allen/``, so the category file and label volume found +*in the directory* were allen32's whatever the study named. Two allen26 +categories matched no parcel and vanished from the region-level tables, Deep +Subcortical was built from 4 of its 8 parcels, and the six merged parcels drew +blank in every mosaic. + +Also locks the header-based 10x convention: the filename rule treated Antwerp's +true-unit label file as inflated and shrank its affine 10x. + +Needs the source-localization atlas data; skipped when it is absent. +""" + +from __future__ import annotations + +import numpy as np +import pytest + +from source_analytics.atlas import ( + find_atlas_dir, + header_is_inflated, + load_atlas, + load_roi_categories, + load_roi_mapping, + registered_atlases, + resolve_atlas, +) + +try: + _HAVE_REGISTRY = (find_atlas_dir() / "registry.yaml").exists() +except FileNotFoundError: + _HAVE_REGISTRY = False + +pytestmark = pytest.mark.skipif(not _HAVE_REGISTRY, + reason="source-localization atlas registry not found") + +# The six bilateral pairs allen26 merges; none of these names exists in allen32. +MERGED = {"Frontal_Anterior", "Olfactory_Bulb", "Thalamus", + "Hypothalamus", "Brainstem_Tectum", "Cerebellum"} + + +def _mapping_names(spec) -> set[str]: + return {v["name"] for v in load_roi_mapping(spec)["rois"].values()} + + +def test_allen26_resolves_to_its_own_files(): + s = resolve_atlas(atlas_name="allen26") + assert s.labels.name == "allen26_labels.nii.gz" + assert s.roi_mapping.name == "roi_mapping_allen26.json" + assert s.roi_categories.name == "roi_categories_allen26.yaml" + + +def test_allen26_categories_and_labels_name_its_merged_parcels(): + s = resolve_atlas(atlas_name="allen26") + members = {r for rois in load_roi_categories(s).values() for r in rois} + assert MERGED <= members, "category map is another atlas's partition" + assert MERGED <= _mapping_names(s), "merged parcels would draw blank" + labels, _ = load_atlas(s) + assert len(np.unique(labels)) - 1 == 26 + + +def test_a_shared_directory_no_longer_decides_the_atlas(): + a26 = resolve_atlas(atlas_name="allen26") + a32 = resolve_atlas(atlas_name="allen32") + assert a26.labels.parent == a32.labels.parent # one directory... + assert a26.labels != a32.labels # ...two atlases + assert a26.roi_categories != a32.roi_categories + + +def test_aliases_and_atlases_without_categories(): + assert resolve_atlas(atlas_name="allen").labels == resolve_atlas(atlas_name="allen32").labels + assert (resolve_atlas(atlas_name="coarse_22roi").labels + == resolve_atlas(atlas_name="coarse22").labels) + a64 = resolve_atlas(atlas_name="allen64") + assert a64.roi_categories is None + assert load_roi_categories(a64) == {} + + +def test_every_registered_atlas_resolves(): + names = registered_atlases() + assert {"allen26", "allen32", "antwerp"} <= set(names) + for name in names: + s = resolve_atlas(atlas_name=name) + assert s.labels.exists() and s.roi_mapping.exists(), name + + +def test_an_unknown_atlas_refuses_to_guess(): + with pytest.raises(ValueError, match="not in"): + resolve_atlas(atlas_name="allen99") + + +def test_explicit_files_define_an_unregistered_atlas(): + ref = resolve_atlas(atlas_name="allen26") + s = resolve_atlas(atlas_name="lab_atlas", files={ + "brain_labels": str(ref.labels), + "roi_mapping": str(ref.roi_mapping), + "roi_categories": str(ref.roi_categories), + }) + assert s.labels == ref.labels + assert load_roi_categories(s) == load_roi_categories(ref) + + +def test_unknown_atlas_files_key_is_rejected(): + with pytest.raises(ValueError, match="unknown key"): + resolve_atlas(atlas_name="allen26", files={"labelz": "x.nii.gz"}) + + +def test_a_spec_passes_through_unchanged(): + s = resolve_atlas(atlas_name="allen26") + assert resolve_atlas(s) is s + + +def test_the_header_decides_the_10x_correction(): + # Antwerp's label file stores true units; the filename rule shrank it 10x. + antwerp = resolve_atlas(atlas_name="antwerp") + assert not header_is_inflated(antwerp.labels) + _, affine = load_atlas(antwerp) + assert affine[0, 0] > 0.1 # 0.203 mm, not 0.0203 + # coarse22's label file IS inflated and must still be corrected. + coarse = resolve_atlas(atlas_name="coarse22") + assert header_is_inflated(coarse.labels) + _, affine = load_atlas(coarse) + assert affine[0, 0] < 0.5 + + +def test_study_default_categories_follow_the_named_atlas(): + from source_analytics.config import _load_atlas_roi_categories + + members = {r for rois in _load_atlas_roi_categories("allen26").values() for r in rois} + assert MERGED <= members diff --git a/tests/test_r_config_categories.py b/tests/test_r_config_categories.py new file mode 100644 index 0000000..065fc7d --- /dev/null +++ b/tests/test_r_config_categories.py @@ -0,0 +1,29 @@ +"""R is handed the EFFECTIVE ROI categories, not whatever ``raw`` happens to hold. + +A profile narrows ``config.roi_categories`` but not ``config.raw``, and the atlas +default never reaches ``raw`` at all. R used to fill that gap from the category +file in the atlas *directory*, which for allen26 was allen32's. +""" + +from __future__ import annotations + +from dataclasses import replace + +from source_analytics.analyses.roi_psd_analysis import ROIPsdAnalysis +from source_analytics.config import StudyConfig + + +def test_r_config_carries_the_narrowed_categories(sample_config_yaml, tmp_path): + cfg = StudyConfig.from_yaml(sample_config_yaml) + motor = {"Motor": list(cfg.roi_categories["Motor"])} + analysis = ROIPsdAnalysis(replace(cfg, roi_categories=motor), tmp_path / "out") + assert analysis._r_config_data()["roi_categories"] == motor + # raw alone still holds the full map -- which is exactly what R used to get. + assert "Subcortical" in cfg.raw["roi_categories"] + + +def test_r_config_leaves_raw_untouched(sample_config_yaml, tmp_path): + cfg = StudyConfig.from_yaml(sample_config_yaml) + before = dict(cfg.raw) + ROIPsdAnalysis(cfg, tmp_path / "out")._r_config_data()["roi_categories"] = {} + assert cfg.raw == before diff --git a/tests/test_r_scripts_smoke.py b/tests/test_r_scripts_smoke.py index e1778d2..df1609b 100644 --- a/tests/test_r_scripts_smoke.py +++ b/tests/test_r_scripts_smoke.py @@ -185,3 +185,64 @@ def test_directed_script_uses_canonical_prefix_and_tests_every_dv(tmp_path, desi assert (tbl / name).exists(), name edges = pd.read_csv(tbl / "roi_directed_directed_edges_hypotheses.csv") assert set(edges["dv"]) == set(cols) + + +# --------------------------------------------------------------------------- +# ROI categories: the study config wins over the --roi-categories file. +# +# The file is looked up in the atlas *directory*, which several atlases share +# (allen32/allen26/allen64 all live in allen/), so for an allen26 study it was +# allen32's partition. It used to REPLACE the study's categories: two allen26 +# categories matched no parcel and vanished, and Deep Subcortical was built +# from 4 of its 8 parcels. +# --------------------------------------------------------------------------- + +def _psd_band_power_csv(data_dir: Path) -> None: + rng = np.random.default_rng(7) + rows = [] + for subj, grp in _subjects(): + for roi in ROIS: + for band in BANDS: + shift = 0.5 if (grp == "KO_VEH" and roi.startswith("Motor")) else 0.0 + rows.append({"subject": subj, "group": grp, "roi": roi, "band": band, + "absolute": float(rng.normal(-200.0 + shift, 1.0)), + "relative": float(rng.uniform(0.1, 0.3))}) + data_dir.mkdir(parents=True, exist_ok=True) + pd.DataFrame(rows).to_csv(data_dir / "band_power.csv", index=False) + + +def _psd_region_labels(out: Path) -> set[str]: + labels: set[str] = set() + for f in (out / "tables").glob("roi_psd*region*.csv"): + d = pd.read_csv(f) + for col in ("spatial", "region"): + if col in d.columns: + labels |= set(d[col].dropna().astype(str)) + return labels + + +def test_psd_region_tier_uses_study_categories_over_the_atlas_file(tmp_path, design_config): + other = tmp_path / "another_atlas_categories.yaml" + other.write_text(yaml.safe_dump({"Bogus": ["Motor_L", "Hipp_L"], "Absent": ["Not_A_Parcel"]})) + data, out = tmp_path / "data", tmp_path / "out" + _psd_band_power_csv(data) + log = _run("roi_psd_analysis.R", data, design_config, out, "--roi-categories", str(other)) + assert "Using roi_categories from the study config" in log + labels = _psd_region_labels(out) + assert labels, "no region-level table was written" + assert labels <= {"Motor", "Hipp", "Thal"}, labels + + +def test_psd_region_tier_falls_back_to_the_file_when_config_has_none(tmp_path, design_config): + cfg = yaml.safe_load(design_config.read_text()) + cfg.pop("roi_categories") + bare = tmp_path / "bare.yaml" + bare.write_text(yaml.safe_dump(cfg)) + cats = tmp_path / "cats.yaml" + cats.write_text(yaml.safe_dump({"Front": ROIS[:2], "Back": ROIS[2:]})) + data, out = tmp_path / "data", tmp_path / "out" + _psd_band_power_csv(data) + log = _run("roi_psd_analysis.R", data, bare, out, "--roi-categories", str(cats)) + assert "Loaded roi_categories from:" in log + labels = _psd_region_labels(out) + assert labels and labels <= {"Front", "Back"}, labels diff --git a/tests/test_roi_signature.py b/tests/test_roi_signature.py new file mode 100644 index 0000000..10a041f --- /dev/null +++ b/tests/test_roi_signature.py @@ -0,0 +1,53 @@ +"""roi_signature: the source-side decoding counterpart of electrode_signature. + +Runs the whole lifecycle on the synthetic study (5 WT vs 5 KO; every subject in +a group is identical and KO carries an added component, so the groups are +separable) and checks the results table, that the features are per ROI, that a +separable contrast decodes, and that figures regenerate from disk alone. +""" + +from __future__ import annotations + +import pandas as pd + +from source_analytics.config import StudyConfig +from source_analytics.core import StudyAnalyzer + +LIFECYCLE = {"setup", "process", "aggregate", "statistics", "figures", "summary"} + + +def _config(sample_config_yaml, tmp_path) -> StudyConfig: + text = sample_config_yaml.read_text() + ( + "\nroi_signature:\n classifiers: [svm_linear]\n n_permutations: 20\n") + path = tmp_path / "signature_config.yaml" + path.write_text(text) + return StudyConfig.from_yaml(path) + + +def test_roi_signature_full_lifecycle(sample_config_yaml, tmp_path): + config = _config(sample_config_yaml, tmp_path) + StudyAnalyzer(config).run_analysis("roi_signature", steps=LIFECYCLE) + + results = pd.read_csv(next(config.results_dir.rglob("roi_signature_results.csv"))) + assert set(results["band"]) == set(config.bands) + assert set(results["classifier"]) == {"svm_linear"} + assert set(results["contrast"]) == {"disease_effect"} + assert results["accuracy"].max() >= 0.9, results[["band", "accuracy"]] + + work = config.output_dir / "roi_signature" + features = pd.read_csv(work / "data" / "roi_signature_features.csv") + assert "roi" in features.columns and "channel" not in features.columns + assert features["roi"].nunique() == 8 + assert (work / "ANALYSIS_SUMMARY.md").read_text().startswith("# ROI Neural Signature") + + +def test_roi_signature_figures_regenerate_from_disk(sample_config_yaml, tmp_path): + config = _config(sample_config_yaml, tmp_path) + analyzer = StudyAnalyzer(config) + analyzer.run_analysis("roi_signature", steps=LIFECYCLE) + pngs = list(config.results_dir.rglob("roi_signature_importance_*.png")) + assert pngs, "no ROI importance figures were drawn" + for p in pngs: + p.unlink() + analyzer.run_analysis("roi_signature", steps={"figures"}) + assert list(config.results_dir.rglob("roi_signature_importance_*.png")) diff --git a/tests/test_signature_lookup.py b/tests/test_signature_lookup.py new file mode 100644 index 0000000..b341c5a --- /dev/null +++ b/tests/test_signature_lookup.py @@ -0,0 +1,64 @@ +"""The source-vs-sensor decoding comparison pairs tables from ONE paradigm. + +Locks out: electrode_signature globbed the whole results tree for +``vertex_signature_results.csv`` and took the first hit. A stale vertex table +from an earlier study version sat in that tree and would have been merged +against the new run's sensor results, labelled as if it were its source side. +""" + +from __future__ import annotations + +from pathlib import Path + +import pandas as pd + +from source_analytics.analyses.electrode_signature_analysis import ( + ElectrodeSignatureAnalysis, + render_signature_comparison, +) + + +class _Sensor(ElectrodeSignatureAnalysis): + """Just enough of the module to exercise the lookup (no config, no IO).""" + + def __init__(self, tbl_dir: Path): + self._tbl = tbl_dir + + @property + def tbl_dir(self) -> Path: + return self._tbl + + +def _table(path: Path, accuracy: float) -> Path: + path.parent.mkdir(parents=True, exist_ok=True) + pd.DataFrame({ + "contrast": ["disease_effect"] * 2, "band": ["Theta", "Beta"], + "classifier": ["svm_linear"] * 2, "accuracy": [accuracy] * 2, + "balanced_accuracy": [accuracy] * 2, "p_value": [0.01, 0.5], "auc": [accuracy] * 2, + }).to_csv(path, index=False) + return path + + +def test_other_paradigms_are_never_searched(tmp_path): + tables = tmp_path / "results" / "tables" + _table(tables / "vertex" / "vertex_signature" / "vertex_signature_results.csv", 0.9) + assert _Sensor(tables / "cartesian_mc" / "electrode_signature")._find_source_results() is None + + +def test_roi_signature_in_the_same_paradigm_is_preferred(tmp_path): + paradigm = tmp_path / "results" / "tables" / "cartesian_mc" + roi = _table(paradigm / "roi_signature" / "roi_signature_results.csv", 0.8) + _table(paradigm / "vertex_signature" / "vertex_signature_results.csv", 0.7) + found = _Sensor(paradigm / "electrode_signature")._find_source_results() + assert found == ("roi_signature", roi) + + +def test_the_comparison_records_which_source_it_used(tmp_path): + source = _table(tmp_path / "source.csv", 0.8) + sensor = _table(tmp_path / "sensor.csv", 0.6) + out = tmp_path / "out" + out.mkdir() + render_signature_comparison(source, sensor, out, out, source_module="roi_signature") + merged = pd.read_csv(out / "signature_source_vs_sensor.csv") + assert set(merged["source_module"]) == {"roi_signature"} + assert merged["accuracy_gain"].round(6).tolist() == [0.2, 0.2] diff --git a/tests/test_signature_threads.py b/tests/test_signature_threads.py new file mode 100644 index 0000000..d0a309c --- /dev/null +++ b/tests/test_signature_threads.py @@ -0,0 +1,41 @@ +"""Signature fits run with the BLAS/OpenMP pools limited to one thread. + +A signature run fits a tiny model LOOCV x (1 + n_permutations) times. With the +default multithreaded BLAS the FORGE treatment logistic fits ran 60-90x slower +(thread start-up and synchronisation on 35 x 30 problems), which turned a module +that takes about an hour into one that takes days. The results are the same. +""" + +from __future__ import annotations + +import numpy as np +import pytest + +pytest.importorskip("sklearn") +threadpoolctl = pytest.importorskip("threadpoolctl") + +from source_analytics.stats.signature import _single_threaded, run_signature # noqa: E402 + + +def test_thread_pools_are_limited_to_one_inside_a_fit(): + @_single_threaded + def probe(): + return [pool["num_threads"] for pool in threadpoolctl.threadpool_info()] + + assert all(n == 1 for n in probe()) + + +def test_run_signature_is_wrapped(): + assert hasattr(run_signature, "__wrapped__") + + +def test_single_threading_does_not_change_the_result(): + rng = np.random.default_rng(0) + X = rng.normal(size=(24, 12)) + y = np.array([0] * 12 + [1] * 12) + X[y == 1, :3] += 1.0 + kw = dict(classifier="logistic", cv_method="loocv", n_permutations=5, seed=42) + limited = run_signature(X, y, **kw) + unlimited = run_signature.__wrapped__(X, y, **kw) + assert limited.accuracy == unlimited.accuracy + assert limited.p_value == unlimited.p_value