From aa6ce3093fc0549601269a606b85e1e77d6fe237 Mon Sep 17 00:00:00 2001 From: Alex Edmondson Date: Thu, 10 Sep 2026 10:24:27 -0400 Subject: [PATCH 1/4] Resolve atlases by name to their own files, and let R keep the study's categories An allen26 study was being analysed with allen32's partition. Several atlases share one directory (allen32, allen26 and allen64 all live in allen/), and every lookup went "the file with the fixed name in the atlas directory": roi_categories.yaml, roi_mapping.json and allen_labels.nii.gz are allen32's. Measured on the FORGE treatment run (allen26, Monte Carlo ROI operators): - the R region tier replaced the study's categories with that file, so Frontal-Anterior and Olfactory matched no parcel and vanished from every region-level table, and Deep Subcortical was built from 4 of its 8 parcels - the ROI mosaics drew allen32's label volume, so the six merged parcels had no voxels and rendered blank (6/26) Nothing warned. ROI-level results were never affected. An atlas is now resolved by NAME through source-localization's registry.yaml, the place an atlas is actually defined, into an AtlasSpec carrying its own labels, mapping, categories, brain volume and brain mask. base.py stores the spec where _atlas_dir lived, and every atlas function accepts a spec wherever it accepted a directory, so the existing call sites pass it through unchanged. An atlas that is not registered names its files under atlas_files:. A name that cannot be resolved raises instead of falling back to whatever shares a directory, which is the fallback that caused this. R now receives the EFFECTIVE categories (the study's map, a profile's narrowing, or the atlas default) through BaseAnalysis._r_config_data(), and a single resolve_roi_categories() in stats_utils.R, replacing seven copies of the overwrite block, lets them win. The --roi-categories file is only a fallback. The 10x voxel convention is now read from the NIfTI header, as source-localization does. The filename rule treated Antwerp's true-unit label file as inflated and shrank its affine 10x on the default-affine path; this moves Antwerp studies' vertex ROI labels and mosaics to their correct places. Tests: atlas resolution for every registered atlas (and the allen26 regression itself), the header rule, explicit atlas_files, the effective categories handed to R, and an R-side test that the study's categories beat a conflicting file. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01RibM8Zep2YEjjUbc2LLkgj --- R/roi_aperiodic_analysis.R | 8 +- R/roi_connectivity_analysis.R | 8 +- R/roi_cross_freq_edges_analysis.R | 9 +- R/roi_evoked_analysis.R | 8 +- R/roi_pac_analysis.R | 8 +- R/roi_psd_analysis.R | 8 +- R/roi_transfer_entropy_analysis.R | 12 +- R/stats_utils.R | 30 +++ src/source_analytics/analyses/base.py | 45 +++- .../analyses/electrode_analysis.py | 2 +- .../analyses/electrode_aperiodic_analysis.py | 2 +- .../analyses/electrode_comparison_analysis.py | 2 +- .../analyses/electrode_evoked_analysis.py | 2 +- .../analyses/roi_aperiodic_analysis.py | 4 +- .../analyses/roi_connectivity_analysis.py | 2 +- .../analyses/roi_cross_freq_analysis.py | 4 +- .../analyses/roi_directed_analysis.py | 2 +- .../analyses/roi_evoked_analysis.py | 2 +- .../analyses/roi_network_analysis.py | 2 +- .../analyses/roi_psd_analysis.py | 4 +- .../analyses/vertex_cluster_analysis.py | 2 +- .../analyses/vertex_connectivity_analysis.py | 2 +- .../analyses/vertex_cross_freq_analysis.py | 2 +- .../analyses/vertex_network_analysis.py | 6 +- .../analyses/vertex_signature_analysis.py | 2 +- .../analyses/vertex_specparam_analysis.py | 2 +- src/source_analytics/atlas/__init__.py | 8 + src/source_analytics/atlas/atlas_utils.py | 253 ++++++++++++++++-- src/source_analytics/config.py | 31 ++- src/source_analytics/io/loader.py | 6 +- src/source_analytics/viz/brain_roi.py | 46 ++-- src/source_analytics/viz/glass_brain.py | 19 +- tests/test_atlas_registry.py | 132 +++++++++ tests/test_r_config_categories.py | 29 ++ tests/test_r_scripts_smoke.py | 61 +++++ 35 files changed, 627 insertions(+), 138 deletions(-) create mode 100644 tests/test_atlas_registry.py create mode 100644 tests/test_r_config_categories.py 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/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/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/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/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/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 From b1c4c82e6ca44959598cda253767a4ef5dcafded Mon Sep 17 00:00:00 2001 From: Alex Edmondson Date: Thu, 10 Sep 2026 10:24:27 -0400 Subject: [PATCH 2/4] Add roi_signature, and pair signatures within one paradigm vertex_signature was the only source-side decoding module, so a study with ROI output alone had no source half of the source-vs-sensor comparison, and Monte Carlo ROI operators never build the vertex estimate it needs. roi_signature runs the same classifiers, cross-validation and permutation test as electrode_signature on per-parcel relative band power, with the identical feature estimator, so where both run on the same epochs their accuracies are like for like. It subclasses electrode_signature, whose table, pickle and figure names now come from the module name rather than a hard-coded prefix. electrode_signature used to find its source partner by globbing the whole results tree for vertex_signature_results.csv and taking the first hit. In the FORGE treatment study that was a stale vertex table from the previous analysis version, which would have been merged against the new run's sensor results as if it were its source side. It now looks only in its own paradigm, preferring roi_signature, and the comparison table records which source module it used. When the sensor signature ran in another paradigm (it reads the raw recordings, so one run can serve every source arm), roi_signature takes sensor_paradigm: and renders the comparison itself. Tests: a full roi_signature lifecycle including figures regenerated from disk, and the in-paradigm lookup (a stale table elsewhere is never used). Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01RibM8Zep2YEjjUbc2LLkgj --- CHANGELOG.md | 33 +++ README.md | 23 +- .../analyses/electrode_signature_analysis.py | 227 ++++++++++-------- .../analyses/roi_signature_analysis.py | 153 ++++++++++++ src/source_analytics/core.py | 5 +- tests/test_roi_signature.py | 53 ++++ tests/test_signature_lookup.py | 64 +++++ 7 files changed, 459 insertions(+), 99 deletions(-) create mode 100644 src/source_analytics/analyses/roi_signature_analysis.py create mode 100644 tests/test_roi_signature.py create mode 100644 tests/test_signature_lookup.py diff --git a/CHANGELOG.md b/CHANGELOG.md index fe2ac9d..efbd007 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,31 @@ 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` stores true + units and was being shrunk 10x on the default-affine path (`load_atlas`, + `load_vertex_roi_labels`, the ROI mosaics); raw-affine extraction was unaffected. + **Antwerp-based studies' vertex ROI labels and mosaics move to their correct + positions.** +- **`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. + - **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 +60,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/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/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_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/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/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] From 33e2a986f75ef9ed1ebc5d9ce496fe50da99eef1 Mon Sep 17 00:00:00 2001 From: Alex Edmondson Date: Thu, 10 Sep 2026 10:35:45 -0400 Subject: [PATCH 3/4] Run signature fits single-threaded A signature run fits a tiny model -- tens of subjects by tens of features -- LOOCV x (1 + n_permutations) times. With the default multithreaded BLAS each of those fits paid for starting and synchronising a thread pool, and that cost dominated. Measured on the FORGE treatment electrode features (35 subjects x 30 electrodes, 21 LOOCV passes): band default threads one thread accuracy Alpha 127.8 s 1.4 s 0.486 in both Low Gamma 105.4 s 1.7 s 0.600 in both That is why one electrode_signature run was heading for several days, and why v1's logistic fits averaged ~16 minutes with a 15-hour outlier. run_signature is now wrapped in threadpoolctl.threadpool_limits(1). threadpoolctl comes with scikit-learn, which the signature modules already need, and the wrapper falls back to a plain call if it is missing. Results do not change: the tests check that the limit is really applied inside a fit and that the limited and unlimited calls agree. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01RibM8Zep2YEjjUbc2LLkgj --- CHANGELOG.md | 5 +++ src/source_analytics/stats/signature.py | 25 +++++++++++++++ tests/test_signature_threads.py | 41 +++++++++++++++++++++++++ 3 files changed, 71 insertions(+) create mode 100644 tests/test_signature_threads.py diff --git a/CHANGELOG.md b/CHANGELOG.md index efbd007..f18218f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -31,6 +31,11 @@ found 24 defects plus a dozen false README claims. All verified and fixed here. `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` 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/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 From 9e4fc8eb7ae9059952de861981f773f0ae24dedf Mon Sep 17 00:00:00 2001 From: Alex Edmondson Date: Thu, 10 Sep 2026 10:50:57 -0400 Subject: [PATCH 4/4] CHANGELOG: correct the Antwerp affine entry The entry (and aa6ce30's message) said reading the 10x convention from the header moves Antwerp studies' vertex ROI labels and mosaics to their correct positions. Measured, it does neither: - load_vertex_roi_labels, the vertex path that used the corrected affine, has never produced a label. It reads the mapping file's top-level keys ("atlas_name", "n_rois", ...) as label ids and raises on every atlas; vertex_network catches that and falls back to spatial node labels. - the mosaics use the affine only for mm axis ranges and slice labels, and the analysis modules never draw them on Antwerp. No statistic changes. ROI extraction and cluster/NBS region labels use the raw affine, which was always right. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01RibM8Zep2YEjjUbc2LLkgj --- CHANGELOG.md | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f18218f..010b60c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,11 +21,16 @@ found 24 defects plus a dozen false README claims. All verified and fixed here. (`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` stores true - units and was being shrunk 10x on the default-affine path (`load_atlas`, - `load_vertex_roi_labels`, the ROI mosaics); raw-affine extraction was unaffected. - **Antwerp-based studies' vertex ROI labels and mosaics move to their correct - positions.** + 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