From 802f7361ecbb9efb3d200d85770b406787f90f87 Mon Sep 17 00:00:00 2001 From: Clara El Khantour Date: Fri, 17 Jul 2026 16:43:49 -0400 Subject: [PATCH 01/33] add site flag and in qcfc --- .../features/quality_control_connectivity.py | 9 +++++++-- wonkyconn/run.py | 7 +++++++ wonkyconn/workflow.py | 18 ++++++++++++++++-- 3 files changed, 30 insertions(+), 4 deletions(-) diff --git a/wonkyconn/features/quality_control_connectivity.py b/wonkyconn/features/quality_control_connectivity.py index dcf9fb02..1d6cf891 100644 --- a/wonkyconn/features/quality_control_connectivity.py +++ b/wonkyconn/features/quality_control_connectivity.py @@ -20,6 +20,7 @@ def calculate_qcfc( data_frame: pd.DataFrame, connectivity_matrices: Iterable[ConnectivityMatrix], metric_key: str = "MeanFramewiseDisplacement", + site_correction: bool = False, ) -> pd.DataFrame: """ metric calculation: quality control / functional connectivity @@ -30,7 +31,8 @@ def calculate_qcfc( accounted for participant age and sex Parameters: - data_frame (pd.DataFrame): The data frame containing the covariates "age" and "gender". + data_frame (pd.DataFrame): The data frame containing the covariates "age" and "gender". + "site" is also required if site correction is applied. It needs to have one row for each connectivity matrix. connectivity_matrices (Iterable[ConnectivityMatrix]): The connectivity matrices to calculate QCFC for. metric_key (str, optional): The key of the metric to use for QCFC calculation. Defaults to "MeanFramewiseDisplacement". @@ -44,7 +46,10 @@ def calculate_qcfc( ) if np.isnan(metrics).all(): raise ValueError(f"None of the connectivity matrices have a metric with key '{metric_key}'") - covariates = np.asarray(dmatrix("age + gender", data_frame)) + if site_correction: + covariates = np.asarray(dmatrix("age + C(gender) + C(site)", data_frame)) + else: + covariates = np.asarray(dmatrix("age + C(gender)", data_frame)) connectivity_arrays = [ connectivity_matrix.load() diff --git a/wonkyconn/run.py b/wonkyconn/run.py index c8982001..154123f8 100644 --- a/wonkyconn/run.py +++ b/wonkyconn/run.py @@ -66,6 +66,13 @@ def global_parser(exit_on_error: bool = True) -> argparse.ArgumentParser: default=False, help="Disable sex and age prediction to reduce runtime.", ) + parser.add_argument( + "--site_correction", + required=False, + action="store_true", + default=False, + help="Apply site correction to the data.", + ) parser.add_argument( "--verbosity", help=""" diff --git a/wonkyconn/workflow.py b/wonkyconn/workflow.py index ef556225..3d2f67d4 100644 --- a/wonkyconn/workflow.py +++ b/wonkyconn/workflow.py @@ -59,6 +59,8 @@ def workflow(args: argparse.Namespace) -> None: # check if light mode is enabled - if so, it will not run the age and sex prediction and gradient similarity disable_prediction_gradient = getattr(args, "light_mode", False) + enable_site_correction = getattr(args, "site_correction", True) + # Check BIDS path bids_dir = args.bids_dir index = BIDSIndex() @@ -138,6 +140,7 @@ def workflow(args: argparse.Namespace) -> None: seg_key, atlases, disable_prediction_gradient, + enable_site_correction ) record.update(dict(zip(group_by, key, strict=False))) if len(group_by) == 2: @@ -171,6 +174,7 @@ def make_record( seg_key: str, atlases: dict[str, Atlas], disable_prediction_gradient: bool, + site_correction: bool = False, ) -> dict[str, Any]: """Compute all QC metrics for a single group of connectivity matrices.""" seg_subjects: list[str] = list() @@ -197,7 +201,7 @@ def make_record( # Slice phenotypes (age, gender, etc.) for just this group seg_data_frame = data_frame.loc[seg_subjects] - qcfc = calculate_qcfc(seg_data_frame, connectivity_matrices, metric_key) + qcfc = calculate_qcfc(seg_data_frame, connectivity_matrices, metric_key, site_correction) (seg,) = index.get_tag_values(seg_key, {c.path for c in connectivity_matrices}) distance_matrix = distance_matrices[seg] @@ -241,6 +245,10 @@ def make_record( try: ages = seg_data_frame["age"].to_numpy() genders = seg_data_frame["gender"].to_numpy() + if site_correction: + # Might need to use get_dummies to convert categorical + # site variable into one-hot encoding for regression + sites = seg_data_frame["site"].to_numpy() scores = age_sex_scores( connectivity_matrices, @@ -285,7 +293,9 @@ def make_record( def load_data_frame(args: argparse.Namespace) -> pd.DataFrame: - """Load a phenotype TSV with ``participant_id``, ``gender``, and ``age`` columns.""" + """Load a phenotype TSV with ``participant_id``, ``gender``, and ``age`` columns. + If site correction is enabled, the ``site`` column is also required. + """ data_frame = pd.read_csv( args.phenotypes, sep="\t", @@ -296,4 +306,8 @@ def load_data_frame(args: argparse.Namespace) -> pd.DataFrame: raise ValueError('Phenotypes file is missing the "gender" column') if "age" not in data_frame.columns: raise ValueError('Phenotypes file is missing the "age" column') + if getattr(args, "site_correction", True): + logger.info("Site correction is enabled - checking for 'site' column in phenotypes file.") + if "site" not in data_frame.columns: + raise ValueError('Phenotypes file is missing the "site" column required for site correction') return data_frame From d322b28f0544bb0a5b1908a08c43c9d6a189036c Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 17 Jul 2026 20:54:21 +0000 Subject: [PATCH 02/33] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- wonkyconn/features/quality_control_connectivity.py | 2 +- wonkyconn/workflow.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/wonkyconn/features/quality_control_connectivity.py b/wonkyconn/features/quality_control_connectivity.py index 1d6cf891..75498b10 100644 --- a/wonkyconn/features/quality_control_connectivity.py +++ b/wonkyconn/features/quality_control_connectivity.py @@ -31,7 +31,7 @@ def calculate_qcfc( accounted for participant age and sex Parameters: - data_frame (pd.DataFrame): The data frame containing the covariates "age" and "gender". + data_frame (pd.DataFrame): The data frame containing the covariates "age" and "gender". "site" is also required if site correction is applied. It needs to have one row for each connectivity matrix. connectivity_matrices (Iterable[ConnectivityMatrix]): The connectivity matrices to calculate QCFC for. diff --git a/wonkyconn/workflow.py b/wonkyconn/workflow.py index 3d2f67d4..d6e1e9f4 100644 --- a/wonkyconn/workflow.py +++ b/wonkyconn/workflow.py @@ -140,7 +140,7 @@ def workflow(args: argparse.Namespace) -> None: seg_key, atlases, disable_prediction_gradient, - enable_site_correction + enable_site_correction, ) record.update(dict(zip(group_by, key, strict=False))) if len(group_by) == 2: @@ -246,7 +246,7 @@ def make_record( ages = seg_data_frame["age"].to_numpy() genders = seg_data_frame["gender"].to_numpy() if site_correction: - # Might need to use get_dummies to convert categorical + # Might need to use get_dummies to convert categorical # site variable into one-hot encoding for regression sites = seg_data_frame["site"].to_numpy() From b63998f9feb9c8e0dd403d62657f55f5ae39c8cb Mon Sep 17 00:00:00 2001 From: Clara El Khantour Date: Mon, 20 Jul 2026 10:20:12 -0400 Subject: [PATCH 03/33] add site_correction in WonkyConnConfig --- wonkyconn/config.py | 3 +++ wonkyconn/tests/test_cli.py | 2 ++ 2 files changed, 5 insertions(+) diff --git a/wonkyconn/config.py b/wonkyconn/config.py index b3bd4b53..291375b9 100644 --- a/wonkyconn/config.py +++ b/wonkyconn/config.py @@ -27,6 +27,7 @@ class WonkyConnConfig: light_mode: bool = False theme: str | None = None # GUI-only suppress_warnings: bool = False + site_correction: bool = False @classmethod def from_cli_args(cls, args: argparse.Namespace | None) -> "WonkyConnConfig": @@ -52,6 +53,7 @@ def from_cli_args(cls, args: argparse.Namespace | None) -> "WonkyConnConfig": debug=bool(getattr(args, "debug", False)), light_mode=bool(getattr(args, "light_mode", False)), suppress_warnings=bool(getattr(args, "suppress_warnings", False)), + site_correction=bool(getattr(args, "site_correction", False)), ) def to_namespace(self) -> argparse.Namespace: @@ -75,4 +77,5 @@ def to_namespace(self) -> argparse.Namespace: verbosity=self.verbosity, debug=self.debug, light_mode=self.light_mode, + site_correction=self.site_correction, ) diff --git a/wonkyconn/tests/test_cli.py b/wonkyconn/tests/test_cli.py index 25343b29..0dd69dc1 100644 --- a/wonkyconn/tests/test_cli.py +++ b/wonkyconn/tests/test_cli.py @@ -85,7 +85,9 @@ def test_cli_and_textual_namespace_consistency(tmp_path: Path): # Get the attribute names from both namespaces cli_attrs = set(vars(cli_args).keys()) + print(f"CLI attributes: {cli_attrs}") config_attrs = set(vars(config_namespace).keys()) + print(f"Config attributes: {config_attrs}") # Attributes that are interface-specific and not passed to workflow() # These are handled separately before calling workflow() From 3ac2dcc42daaf2ca88b265f1bc8e685d60d1b180 Mon Sep 17 00:00:00 2001 From: Clara El Khantour Date: Mon, 20 Jul 2026 10:21:41 -0400 Subject: [PATCH 04/33] Remove print in test --- wonkyconn/tests/test_cli.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/wonkyconn/tests/test_cli.py b/wonkyconn/tests/test_cli.py index 0dd69dc1..25343b29 100644 --- a/wonkyconn/tests/test_cli.py +++ b/wonkyconn/tests/test_cli.py @@ -85,9 +85,7 @@ def test_cli_and_textual_namespace_consistency(tmp_path: Path): # Get the attribute names from both namespaces cli_attrs = set(vars(cli_args).keys()) - print(f"CLI attributes: {cli_attrs}") config_attrs = set(vars(config_namespace).keys()) - print(f"Config attributes: {config_attrs}") # Attributes that are interface-specific and not passed to workflow() # These are handled separately before calling workflow() From ad50da7441fc0d65fe0fb247a16f0e4eae54145c Mon Sep 17 00:00:00 2001 From: Clara El Khantour Date: Mon, 20 Jul 2026 10:24:37 -0400 Subject: [PATCH 05/33] Rename variable for pre-commit --- wonkyconn/features/network.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/wonkyconn/features/network.py b/wonkyconn/features/network.py index 5b5b5616..2a2c895d 100644 --- a/wonkyconn/features/network.py +++ b/wonkyconn/features/network.py @@ -45,8 +45,8 @@ def single_subject_within_network_connectivity( roi_index -= 1 # calculate similarity of the thresholded individual level seed based connectivity with the binary mask - subj_average_connectivity_within_network = [] - subj_variance_connectivity_within_network = [] + subj_average_connectivity_within_network_list = [] + subj_variance_connectivity_within_network_list = [] subj_corr_with_network = [] for idx_roi in roi_index: seed_based_map = connectivity_matrix.load()[int(idx_roi), :] @@ -72,13 +72,13 @@ def single_subject_within_network_connectivity( seed_based_map[mask], region_membership[f"yeo7-{yeo_network_index}"].to_numpy()[mask] ) - subj_average_connectivity_within_network.append(mean_within_network_connection) - subj_variance_connectivity_within_network.append(std_within_network_connection) + subj_average_connectivity_within_network_list.append(mean_within_network_connection) + subj_variance_connectivity_within_network_list.append(std_within_network_connection) subj_corr_with_network.append(correlation_with_given_network[1, 0]) # summarise of the given subject - subj_average_connectivity_within_network = np.asarray(subj_average_connectivity_within_network).mean(axis=0) - subj_variance_connectivity_within_network = np.asarray(subj_variance_connectivity_within_network).mean(axis=0) + subj_average_connectivity_within_network = np.asarray(subj_average_connectivity_within_network_list).mean(axis=0) + subj_variance_connectivity_within_network = np.asarray(subj_variance_connectivity_within_network_list).mean(axis=0) return ( subj_average_connectivity_within_network, subj_variance_connectivity_within_network, From b34a6b47a21c4bf79cbc712112491fc603900d34 Mon Sep 17 00:00:00 2001 From: Clara El Khantour Date: Mon, 20 Jul 2026 11:18:39 -0400 Subject: [PATCH 06/33] Remove site for prediction for now --- wonkyconn/workflow.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/wonkyconn/workflow.py b/wonkyconn/workflow.py index d6e1e9f4..527d1120 100644 --- a/wonkyconn/workflow.py +++ b/wonkyconn/workflow.py @@ -245,10 +245,6 @@ def make_record( try: ages = seg_data_frame["age"].to_numpy() genders = seg_data_frame["gender"].to_numpy() - if site_correction: - # Might need to use get_dummies to convert categorical - # site variable into one-hot encoding for regression - sites = seg_data_frame["site"].to_numpy() scores = age_sex_scores( connectivity_matrices, From 02d2423e55c0a0567bcb94b3eb9a470b7d25e546 Mon Sep 17 00:00:00 2001 From: Clara El Khantour Date: Wed, 22 Jul 2026 11:16:52 -0400 Subject: [PATCH 07/33] draft site correction for prediction --- wonkyconn/features/age_sex_prediction.py | 115 +++++++++++++++++++---- wonkyconn/workflow.py | 4 +- 2 files changed, 101 insertions(+), 18 deletions(-) diff --git a/wonkyconn/features/age_sex_prediction.py b/wonkyconn/features/age_sex_prediction.py index d2d82420..b7d6d6d7 100644 --- a/wonkyconn/features/age_sex_prediction.py +++ b/wonkyconn/features/age_sex_prediction.py @@ -13,10 +13,51 @@ from sklearn.model_selection import StratifiedShuffleSplit, cross_validate # type: ignore[import-not-found] from sklearn.pipeline import Pipeline # type: ignore[import-not-found] from sklearn.preprocessing import LabelEncoder, StandardScaler # type: ignore[import-not-found] +from sklearn.metrics import ( # type: ignore[import-not-found] + accuracy_score, + mean_absolute_error, + r2_score, + roc_auc_score, +) if TYPE_CHECKING: from ..base import ConnectivityMatrix +def regress_site( + X_train: NDArray[np.float32], + X_test: NDArray[np.float32], + site_train: NDArray[np.str_], + site_test: NDArray[np.str_], +) -> tuple[NDArray[np.float32], NDArray[np.float32]]: + """Regress out site effects from the training and test data.""" + + train_site = pd.get_dummies(site_train, drop_first=True, dtype=float) + test_site = pd.get_dummies(site_test, drop_first=True, dtype=float) + + test_site = test_site.reindex( + columns=train_site.columns, + fill_value=0.0, + ) + + design_train = np.column_stack( + [np.ones(len(train_site)), train_site.to_numpy()] + ) + + design_test = np.column_stack( + [np.ones(len(test_site)), test_site.to_numpy()] + ) + + beta, *_ = np.linalg.lstsq( + design_train, + X_train, + rcond=None, + ) + + beta_site = beta[1:, :] + X_train_corr = X_train - design_train[:, 1:] @ beta_site + X_test_corr = X_test - design_test[:, 1:] @ beta_site + + return X_train_corr, X_test_corr def training_pipeline( connectivity_data: NDArray[np.float32], @@ -26,6 +67,7 @@ def training_pipeline( n_pca: int, n_jobs: int = 4, random_state: int = 1, + sites: NDArray[np.str_] | None = None, ) -> pd.DataFrame: """Runs a cross-validation pipeline for age or sex prediction. @@ -37,6 +79,7 @@ def training_pipeline( n_pca (int): Number of principal components to extract. n_jobs (int): Number of cores for parallel calculation. random_state (int): Seed for reproducibility. + sites (NDArray[np.str_] | None): Site labels for the data. Returns: pd.DataFrame: Statistics (mean, 95% CI) of the scores obtained. @@ -44,18 +87,19 @@ def training_pipeline( connectivity_data = np.asarray(connectivity_data, dtype=np.float32, order="C") if task_type == "classification": - y_train = LabelEncoder().fit_transform(target_labels) + y = LabelEncoder().fit_transform(target_labels) estimator = LogisticRegression(max_iter=5000, solver="saga", penalty="l2", n_jobs=n_jobs, random_state=random_state) cv_strategy = StratifiedShuffleSplit(n_splits=n_splits, test_size=0.2, random_state=random_state) scoring_metrics = {"accuracy": "accuracy", "roc_auc": "roc_auc"} + splits = cv_strategy.split(connectivity_data, y) + else: - y_train = np.asarray(target_labels) + y = np.asarray(target_labels) estimator = Ridge(alpha=1.0) - bins = pd.qcut(y_train, q=5, labels=False, duplicates="drop") + bins = pd.qcut(y, q=5, labels=False, duplicates="drop") cv_strategy = StratifiedShuffleSplit(n_splits=n_splits, test_size=0.2, random_state=random_state) splits = list(cv_strategy.split(np.zeros_like(bins), bins)) - cv_strategy = splits scoring_metrics = {"mae": "neg_mean_absolute_error", "r2": "r2"} @@ -67,21 +111,54 @@ def training_pipeline( ("estimator", estimator), ] ) + + scores = {metric: [] for metric in scoring_metrics} with parallel_backend("threading", n_jobs=n_jobs): - cv_results = cross_validate( - pipe, - connectivity_data, - y_train, - cv=cv_strategy, - scoring=scoring_metrics, - n_jobs=n_jobs, - ) - - scores_df = pd.DataFrame({k.replace("test_", ""): v for k, v in cv_results.items() if k.startswith("test_")}) + for train_idx, test_idx in splits: + X_train = connectivity_data[train_idx] + X_test = connectivity_data[test_idx] + + y_train = y[train_idx] + y_test = y[test_idx] + + if sites is not None: + X_train, X_test = regress_site( + X_train, + X_test, + sites[train_idx], + sites[test_idx], + ) + + pipe.fit(X_train, y_train) + y_pred = pipe.predict(X_test) + + if task_type == "classification": + scores["accuracy"].append( + accuracy_score(y_test, y_pred) + ) + + y_prob = pipe.predict_proba(X_test)[:, 1] + + scores["roc_auc"].append( + roc_auc_score(y_test, y_prob) + ) + + else: + scores["mae"].append( + mean_absolute_error(y_test, y_pred) + ) + + scores["r2"].append( + r2_score(y_test, y_pred) + ) + + scores_df = pd.DataFrame(scores) + summary = scores_df.agg(["mean"]).T summary["ci_lower"] = scores_df.quantile(0.025) summary["ci_upper"] = scores_df.quantile(0.975) + return summary @@ -89,6 +166,7 @@ def age_sex_scores( connectivity_matrices: List[ConnectivityMatrix], ages: NDArray[np.float64], genders: NDArray[np.str_], + sites: NDArray[np.str_] | None, n_splits: int, n_pca: int, n_jobs: int = 4, @@ -100,6 +178,7 @@ def age_sex_scores( connectivity_matrices (List[ConnectivityMatrix]): List of matrix objects. ages: Vector of subject ages. genders: Vector of subject genders. + sites: Vector of subject sites. n_splits (int): Number of splits for cross-validation. n_pca (int): Number of PCA components. n_jobs (int): Number of joblib threads. @@ -119,6 +198,7 @@ def age_sex_scores( n_pca=n_pca, n_jobs=n_jobs, random_state=random_state, + sites=sites, ) age_summary = training_pipeline( @@ -129,6 +209,7 @@ def age_sex_scores( n_pca=n_pca, n_jobs=n_jobs, random_state=random_state, + sites=sites, ) return { @@ -136,8 +217,8 @@ def age_sex_scores( "sex_auc_ci_lower": float(sex_summary.loc["roc_auc", "ci_lower"]), # type: ignore[arg-type] "sex_auc_ci_upper": float(sex_summary.loc["roc_auc", "ci_upper"]), # type: ignore[arg-type] "sex_accuracy": float(sex_summary.loc["accuracy", "mean"]), # type: ignore[arg-type] - "age_mae": float(-age_summary.loc["mae", "mean"]), # type: ignore[arg-type, operator] - "age_mae_ci_lower": float(-age_summary.loc["mae", "ci_upper"]), # type: ignore[arg-type, operator] - "age_mae_ci_upper": float(-age_summary.loc["mae", "ci_lower"]), # type: ignore[arg-type, operator] + "age_mae": float(age_summary.loc["mae", "mean"]), # type: ignore[arg-type, operator] + "age_mae_ci_lower": float(age_summary.loc["mae", "ci_lower"]), # type: ignore[arg-type, operator] + "age_mae_ci_upper": float(age_summary.loc["mae", "ci_upper"]), # type: ignore[arg-type, operator] "age_r2": float(age_summary.loc["r2", "mean"]), # type: ignore[arg-type] } diff --git a/wonkyconn/workflow.py b/wonkyconn/workflow.py index 527d1120..1ff463f4 100644 --- a/wonkyconn/workflow.py +++ b/wonkyconn/workflow.py @@ -154,7 +154,7 @@ def workflow(args: argparse.Namespace) -> None: record["dmn_similarity_mean"] = dmn_similarity_avg records.append(record) - + plot(records, group_by, output_dir) for record in records: @@ -245,11 +245,13 @@ def make_record( try: ages = seg_data_frame["age"].to_numpy() genders = seg_data_frame["gender"].to_numpy() + sites = seg_data_frame["site"].to_numpy() if site_correction else None scores = age_sex_scores( connectivity_matrices, ages=ages, genders=genders, + sites=sites, n_splits=_DEFAULT_N_SPLITS, random_state=42, n_pca=_DEFAULT_N_PCA, From 80bc3e8513cfd8148d2dd7d722d8717f3ad0f716 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 22 Jul 2026 15:17:14 +0000 Subject: [PATCH 08/33] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- wonkyconn/features/age_sex_prediction.py | 38 +++++++++--------------- wonkyconn/workflow.py | 2 +- 2 files changed, 15 insertions(+), 25 deletions(-) diff --git a/wonkyconn/features/age_sex_prediction.py b/wonkyconn/features/age_sex_prediction.py index b7d6d6d7..d7b8b489 100644 --- a/wonkyconn/features/age_sex_prediction.py +++ b/wonkyconn/features/age_sex_prediction.py @@ -10,19 +10,20 @@ from sklearn.decomposition import PCA # type: ignore[import-not-found] from sklearn.impute import SimpleImputer # type: ignore[import-not-found] from sklearn.linear_model import LogisticRegression, Ridge # type: ignore[import-not-found] -from sklearn.model_selection import StratifiedShuffleSplit, cross_validate # type: ignore[import-not-found] -from sklearn.pipeline import Pipeline # type: ignore[import-not-found] -from sklearn.preprocessing import LabelEncoder, StandardScaler # type: ignore[import-not-found] -from sklearn.metrics import ( # type: ignore[import-not-found] +from sklearn.metrics import ( # type: ignore[import-not-found] accuracy_score, mean_absolute_error, r2_score, roc_auc_score, ) +from sklearn.model_selection import StratifiedShuffleSplit # type: ignore[import-not-found] +from sklearn.pipeline import Pipeline # type: ignore[import-not-found] +from sklearn.preprocessing import LabelEncoder, StandardScaler # type: ignore[import-not-found] if TYPE_CHECKING: from ..base import ConnectivityMatrix + def regress_site( X_train: NDArray[np.float32], X_test: NDArray[np.float32], @@ -39,13 +40,9 @@ def regress_site( fill_value=0.0, ) - design_train = np.column_stack( - [np.ones(len(train_site)), train_site.to_numpy()] - ) + design_train = np.column_stack([np.ones(len(train_site)), train_site.to_numpy()]) - design_test = np.column_stack( - [np.ones(len(test_site)), test_site.to_numpy()] - ) + design_test = np.column_stack([np.ones(len(test_site)), test_site.to_numpy()]) beta, *_ = np.linalg.lstsq( design_train, @@ -59,6 +56,7 @@ def regress_site( return X_train_corr, X_test_corr + def training_pipeline( connectivity_data: NDArray[np.float32], target_labels: NDArray[np.float64] | NDArray[np.str_], @@ -92,7 +90,7 @@ def training_pipeline( cv_strategy = StratifiedShuffleSplit(n_splits=n_splits, test_size=0.2, random_state=random_state) scoring_metrics = {"accuracy": "accuracy", "roc_auc": "roc_auc"} splits = cv_strategy.split(connectivity_data, y) - + else: y = np.asarray(target_labels) estimator = Ridge(alpha=1.0) @@ -111,7 +109,7 @@ def training_pipeline( ("estimator", estimator), ] ) - + scores = {metric: [] for metric in scoring_metrics} with parallel_backend("threading", n_jobs=n_jobs): @@ -134,24 +132,16 @@ def training_pipeline( y_pred = pipe.predict(X_test) if task_type == "classification": - scores["accuracy"].append( - accuracy_score(y_test, y_pred) - ) + scores["accuracy"].append(accuracy_score(y_test, y_pred)) y_prob = pipe.predict_proba(X_test)[:, 1] - scores["roc_auc"].append( - roc_auc_score(y_test, y_prob) - ) + scores["roc_auc"].append(roc_auc_score(y_test, y_prob)) else: - scores["mae"].append( - mean_absolute_error(y_test, y_pred) - ) + scores["mae"].append(mean_absolute_error(y_test, y_pred)) - scores["r2"].append( - r2_score(y_test, y_pred) - ) + scores["r2"].append(r2_score(y_test, y_pred)) scores_df = pd.DataFrame(scores) diff --git a/wonkyconn/workflow.py b/wonkyconn/workflow.py index 1ff463f4..2ce00aa5 100644 --- a/wonkyconn/workflow.py +++ b/wonkyconn/workflow.py @@ -154,7 +154,7 @@ def workflow(args: argparse.Namespace) -> None: record["dmn_similarity_mean"] = dmn_similarity_avg records.append(record) - + plot(records, group_by, output_dir) for record in records: From 59d4e53f574785fbd1a11010f856086262a6e32f Mon Sep 17 00:00:00 2001 From: Clara El Khantour Date: Thu, 23 Jul 2026 09:10:39 -0400 Subject: [PATCH 09/33] Fit both train and test --- wonkyconn/features/age_sex_prediction.py | 17 ++++++----------- 1 file changed, 6 insertions(+), 11 deletions(-) diff --git a/wonkyconn/features/age_sex_prediction.py b/wonkyconn/features/age_sex_prediction.py index d7b8b489..f962017d 100644 --- a/wonkyconn/features/age_sex_prediction.py +++ b/wonkyconn/features/age_sex_prediction.py @@ -44,15 +44,10 @@ def regress_site( design_test = np.column_stack([np.ones(len(test_site)), test_site.to_numpy()]) - beta, *_ = np.linalg.lstsq( - design_train, - X_train, - rcond=None, - ) - - beta_site = beta[1:, :] - X_train_corr = X_train - design_train[:, 1:] @ beta_site - X_test_corr = X_test - design_test[:, 1:] @ beta_site + beta_train, *_ = np.linalg.lstsq(design_train, X_train, rcond=None,) + beta_test, *_ = np.linalg.lstsq(design_test, X_test, rcond=None,) + X_train_corr = X_train - design_train[:, 1:] @ beta_train[1:, :] + X_test_corr = X_test - design_test[:, 1:] @ beta_test[1:, :] return X_train_corr, X_test_corr @@ -88,7 +83,7 @@ def training_pipeline( y = LabelEncoder().fit_transform(target_labels) estimator = LogisticRegression(max_iter=5000, solver="saga", penalty="l2", n_jobs=n_jobs, random_state=random_state) cv_strategy = StratifiedShuffleSplit(n_splits=n_splits, test_size=0.2, random_state=random_state) - scoring_metrics = {"accuracy": "accuracy", "roc_auc": "roc_auc"} + scoring_metrics = ["accuracy", "roc_auc"] splits = cv_strategy.split(connectivity_data, y) else: @@ -99,7 +94,7 @@ def training_pipeline( cv_strategy = StratifiedShuffleSplit(n_splits=n_splits, test_size=0.2, random_state=random_state) splits = list(cv_strategy.split(np.zeros_like(bins), bins)) - scoring_metrics = {"mae": "neg_mean_absolute_error", "r2": "r2"} + scoring_metrics = ["mae", "r2"] pipe = Pipeline( [ From d8edae27fd294340bde1f5ea712af46504d61c9e Mon Sep 17 00:00:00 2001 From: Clara El Khantour Date: Fri, 24 Jul 2026 16:33:15 -0400 Subject: [PATCH 10/33] Refactor site correction for prediction --- wonkyconn/features/age_sex_prediction.py | 173 +++++++++++++---------- 1 file changed, 98 insertions(+), 75 deletions(-) diff --git a/wonkyconn/features/age_sex_prediction.py b/wonkyconn/features/age_sex_prediction.py index f962017d..97c32131 100644 --- a/wonkyconn/features/age_sex_prediction.py +++ b/wonkyconn/features/age_sex_prediction.py @@ -7,16 +7,11 @@ from joblib import parallel_backend # type: ignore[import-not-found] from nilearn.connectome import sym_matrix_to_vec # type: ignore[import-not-found] from numpy.typing import NDArray +from sklearn.base import BaseEstimator, TransformerMixin # type: ignore[import-not-found] from sklearn.decomposition import PCA # type: ignore[import-not-found] from sklearn.impute import SimpleImputer # type: ignore[import-not-found] from sklearn.linear_model import LogisticRegression, Ridge # type: ignore[import-not-found] -from sklearn.metrics import ( # type: ignore[import-not-found] - accuracy_score, - mean_absolute_error, - r2_score, - roc_auc_score, -) -from sklearn.model_selection import StratifiedShuffleSplit # type: ignore[import-not-found] +from sklearn.model_selection import StratifiedShuffleSplit, cross_validate # type: ignore[import-not-found] from sklearn.pipeline import Pipeline # type: ignore[import-not-found] from sklearn.preprocessing import LabelEncoder, StandardScaler # type: ignore[import-not-found] @@ -24,32 +19,60 @@ from ..base import ConnectivityMatrix -def regress_site( - X_train: NDArray[np.float32], - X_test: NDArray[np.float32], - site_train: NDArray[np.str_], - site_test: NDArray[np.str_], -) -> tuple[NDArray[np.float32], NDArray[np.float32]]: - """Regress out site effects from the training and test data.""" +class SiteRegressor(BaseEstimator, TransformerMixin): + """ + Regress out site effects from connectivity features. + """ - train_site = pd.get_dummies(site_train, drop_first=True, dtype=float) - test_site = pd.get_dummies(site_test, drop_first=True, dtype=float) + def __init__(self: SiteRegressor, n_connectivity_features: int) -> None: + """ + Args: + n_connectivity_features: Number of connectivity features in the input data. + """ + self.n_connectivity_features = n_connectivity_features - test_site = test_site.reindex( - columns=train_site.columns, - fill_value=0.0, - ) + def fit(self: SiteRegressor, connectivity_data_site: NDArray[np.float32], y: None = None) -> SiteRegressor: + """ + Fit the site regressor to the connectivity data. + + Args: + connectivity_data_site: A 2D array where the first n_connectivity_features columns + are connectivity features and the remaining columns are site dummy variables. + y: Ignored. This parameter exists for compatibility with the scikit-learn API. + + Returns: + self: Returns the instance itself. + """ + connectivity_feat = connectivity_data_site[:, : self.n_connectivity_features] + site = connectivity_data_site[:, self.n_connectivity_features :] + + # Add intercept + design = np.column_stack([np.ones(len(site)), site]) + + # Estimate coefficients on the training fold only + self.beta_ = np.linalg.lstsq(design, connectivity_feat, rcond=None)[0] - design_train = np.column_stack([np.ones(len(train_site)), train_site.to_numpy()]) + return self - design_test = np.column_stack([np.ones(len(test_site)), test_site.to_numpy()]) + def transform(self: SiteRegressor, connectivity_data_site: NDArray[np.float32]) -> NDArray[np.float64]: + """ + Transform the connectivity data by regressing out site effects. + Args: + connectivity_data_site: A 2D array where the first n_connectivity_features columns + are connectivity features and the remaining columns are site dummy variables. - beta_train, *_ = np.linalg.lstsq(design_train, X_train, rcond=None,) - beta_test, *_ = np.linalg.lstsq(design_test, X_test, rcond=None,) - X_train_corr = X_train - design_train[:, 1:] @ beta_train[1:, :] - X_test_corr = X_test - design_test[:, 1:] @ beta_test[1:, :] + Returns: + A 2D array of connectivity features with site effects regressed out. + """ - return X_train_corr, X_test_corr + connectivity_feat = connectivity_data_site[:, : self.n_connectivity_features] + site = connectivity_data_site[:, self.n_connectivity_features :] + + # Add intercept + design = np.column_stack([np.ones(len(site)), site]) + + # Remove only the site contribution (keep the intercept) + return connectivity_feat - design[:, 1:] @ self.beta_[1:] def training_pipeline( @@ -79,71 +102,71 @@ def training_pipeline( """ connectivity_data = np.asarray(connectivity_data, dtype=np.float32, order="C") + # Transform site labels into dummy variables and concatenate with connectivity data + if sites is not None: + site = pd.get_dummies( + sites, + drop_first=True, + dtype=np.float32, + ).to_numpy() + + connectivity_data_site = np.concatenate([connectivity_data, site], axis=1) + else: + connectivity_data_site = connectivity_data + if task_type == "classification": - y = LabelEncoder().fit_transform(target_labels) + y_train = LabelEncoder().fit_transform(target_labels) estimator = LogisticRegression(max_iter=5000, solver="saga", penalty="l2", n_jobs=n_jobs, random_state=random_state) cv_strategy = StratifiedShuffleSplit(n_splits=n_splits, test_size=0.2, random_state=random_state) - scoring_metrics = ["accuracy", "roc_auc"] - splits = cv_strategy.split(connectivity_data, y) - + scoring_metrics = {"accuracy": "accuracy", "roc_auc": "roc_auc"} else: - y = np.asarray(target_labels) + y_train = np.asarray(target_labels) estimator = Ridge(alpha=1.0) - bins = pd.qcut(y, q=5, labels=False, duplicates="drop") + bins = pd.qcut(y_train, q=5, labels=False, duplicates="drop") cv_strategy = StratifiedShuffleSplit(n_splits=n_splits, test_size=0.2, random_state=random_state) splits = list(cv_strategy.split(np.zeros_like(bins), bins)) + cv_strategy = splits + + scoring_metrics = {"mae": "neg_mean_absolute_error", "r2": "r2"} - scoring_metrics = ["mae", "r2"] + steps = [ + ("imputer", SimpleImputer(strategy="median")), + ] - pipe = Pipeline( + if sites is not None: + steps.append(("site_regression", SiteRegressor(connectivity_data.shape[1]))) + + steps.extend( [ - ("imputer", SimpleImputer(strategy="median")), ("scaler", StandardScaler()), - ("pca", PCA(n_components=n_pca, svd_solver="randomized", random_state=random_state)), + ( + "pca", + PCA( + n_components=n_pca, + svd_solver="randomized", + random_state=random_state, + ), + ), ("estimator", estimator), ] ) - scores = {metric: [] for metric in scoring_metrics} - + pipe = Pipeline(steps) with parallel_backend("threading", n_jobs=n_jobs): - for train_idx, test_idx in splits: - X_train = connectivity_data[train_idx] - X_test = connectivity_data[test_idx] - - y_train = y[train_idx] - y_test = y[test_idx] - - if sites is not None: - X_train, X_test = regress_site( - X_train, - X_test, - sites[train_idx], - sites[test_idx], - ) - - pipe.fit(X_train, y_train) - y_pred = pipe.predict(X_test) - - if task_type == "classification": - scores["accuracy"].append(accuracy_score(y_test, y_pred)) - - y_prob = pipe.predict_proba(X_test)[:, 1] - - scores["roc_auc"].append(roc_auc_score(y_test, y_prob)) - - else: - scores["mae"].append(mean_absolute_error(y_test, y_pred)) - - scores["r2"].append(r2_score(y_test, y_pred)) - - scores_df = pd.DataFrame(scores) - + cv_results = cross_validate( + pipe, + connectivity_data_site, + y_train, + cv=cv_strategy, + scoring=scoring_metrics, + n_jobs=n_jobs, + ) + + scores_df = pd.DataFrame({k.replace("test_", ""): v for k, v in cv_results.items() if k.startswith("test_")}) summary = scores_df.agg(["mean"]).T summary["ci_lower"] = scores_df.quantile(0.025) summary["ci_upper"] = scores_df.quantile(0.975) - return summary @@ -202,8 +225,8 @@ def age_sex_scores( "sex_auc_ci_lower": float(sex_summary.loc["roc_auc", "ci_lower"]), # type: ignore[arg-type] "sex_auc_ci_upper": float(sex_summary.loc["roc_auc", "ci_upper"]), # type: ignore[arg-type] "sex_accuracy": float(sex_summary.loc["accuracy", "mean"]), # type: ignore[arg-type] - "age_mae": float(age_summary.loc["mae", "mean"]), # type: ignore[arg-type, operator] - "age_mae_ci_lower": float(age_summary.loc["mae", "ci_lower"]), # type: ignore[arg-type, operator] - "age_mae_ci_upper": float(age_summary.loc["mae", "ci_upper"]), # type: ignore[arg-type, operator] + "age_mae": float(-age_summary.loc["mae", "mean"]), # type: ignore[arg-type, operator] + "age_mae_ci_lower": float(-age_summary.loc["mae", "ci_upper"]), # type: ignore[arg-type, operator] + "age_mae_ci_upper": float(-age_summary.loc["mae", "ci_lower"]), # type: ignore[arg-type, operator] "age_r2": float(age_summary.loc["r2", "mean"]), # type: ignore[arg-type] } From 2ea2be86627bae3ae1d7195aae6d995f1be815a6 Mon Sep 17 00:00:00 2001 From: Lea Waller Date: Sun, 2 Aug 2026 10:52:50 +0200 Subject: [PATCH 11/33] Fix typing issues for mypy and pyright --- .pre-commit-config.yaml | 2 +- pyproject.toml | 5 +++++ wonkyconn/atlas.py | 16 ++++++++-------- wonkyconn/features/age_sex_prediction.py | 18 +++++++++--------- .../calculate_gradients_correlation.py | 8 ++++---- wonkyconn/features/network.py | 2 +- wonkyconn/run.py | 2 +- wonkyconn/textual_app.py | 14 +++++++------- 8 files changed, 36 insertions(+), 31 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 3a7f90b2..877304a1 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -20,7 +20,7 @@ repos: rev: v1.19.1 hooks: - id: mypy - additional_dependencies: [pandas-stubs, types-tqdm, types-setuptools, types-Jinja2] + additional_dependencies: [pandas-stubs, types-tqdm, types-setuptools, types-Jinja2, textual] args: [--config-file=pyproject.toml] - repo: https://github.com/codespell-project/codespell rev: v2.4.2 diff --git a/pyproject.toml b/pyproject.toml index f162fb4c..fc1dd056 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -164,13 +164,18 @@ warn_unused_ignores = true ignore_missing_imports = true module = [ "bids.*", + "brainspace.*", "datalad.*", + "joblib", "matplotlib.*", + "nilearn", + "nilearn.*", "numba.*", "patsy.*", "rich.*", "scipy.*", "seaborn.*", + "sklearn.*", "statsmodels.*", "templateflow.*", "nibabel.*", diff --git a/wonkyconn/atlas.py b/wonkyconn/atlas.py index 3c28807f..3fa8cc86 100644 --- a/wonkyconn/atlas.py +++ b/wonkyconn/atlas.py @@ -7,8 +7,8 @@ import numpy as np import pandas as pd import scipy -from nilearn.image import iter_img, load_img, math_img, resample_to_img # type: ignore[import-not-found] -from nilearn.maskers import NiftiLabelsMasker, NiftiMasker # type: ignore[import-not-found] +from nilearn.image import iter_img, load_img, math_img, resample_to_img +from nilearn.maskers import NiftiLabelsMasker, NiftiMasker from numpy import typing as npt from .logger import logger @@ -29,7 +29,7 @@ class Atlas(ABC): """ seg: str - image: nib.nifti1.Nifti1Image + image: nib.nifti1.Nifti1Image # pyright: ignore[reportAttributeAccessIssue] structure: npt.NDArray[np.bool_] = field(default_factory=lambda: np.ones((3, 3, 3), dtype=bool)) @@ -51,7 +51,7 @@ def get_centroids(self) -> npt.NDArray[np.float64]: npt.NDArray[np.float64]: An array of centroid coordinates. """ centroid_points = self.get_centroid_points() - centroid_coordinates = nib.affines.apply_affine(self.image.affine, centroid_points) + centroid_coordinates = nib.affines.apply_affine(self.image.affine, centroid_points) # pyright: ignore[reportAttributeAccessIssue] return centroid_coordinates def get_distance_matrix(self) -> npt.NDArray[np.float64]: @@ -65,7 +65,7 @@ def get_distance_matrix(self) -> npt.NDArray[np.float64]: centroids = self.get_centroids() return scipy.spatial.distance.squareform(scipy.spatial.distance.pdist(centroids)) - def load_yeo7_network(self) -> nib.nifti1.Nifti1Image: + def load_yeo7_network(self) -> nib.nifti1.Nifti1Image: # pyright: ignore[reportAttributeAccessIssue] """ Load and resample the yeo 7 networks to the atlas's space. @@ -102,10 +102,10 @@ def create(seg: str, path: Path) -> "Atlas": None """ - image = nib.nifti1.load(path) + image = nib.nifti1.load(path) # pyright: ignore[reportAttributeAccessIssue] if image.ndim <= 3 or image.shape[3] == 1: - return DsegAtlas(seg, nib.funcs.squeeze_image(image)) + return DsegAtlas(seg, nib.funcs.squeeze_image(image)) # pyright: ignore[reportAttributeAccessIssue] else: return ProbsegAtlas(seg, image) @@ -161,7 +161,7 @@ def _get_centroid_point(self, i: int, array: npt.NDArray[np.float64]) -> tuple[f def get_centroid_points(self) -> npt.NDArray[np.float64]: return np.asarray( - [self._get_centroid_point(i, image.get_fdata()) for i, image in enumerate(nib.funcs.four_to_three(self.image))] + [self._get_centroid_point(i, image.get_fdata()) for i, image in enumerate(nib.funcs.four_to_three(self.image))] # pyright: ignore[reportAttributeAccessIssue] ) def get_yeo7_membership(self) -> pd.DataFrame: diff --git a/wonkyconn/features/age_sex_prediction.py b/wonkyconn/features/age_sex_prediction.py index 97c32131..324909d6 100644 --- a/wonkyconn/features/age_sex_prediction.py +++ b/wonkyconn/features/age_sex_prediction.py @@ -4,16 +4,16 @@ import numpy as np import pandas as pd -from joblib import parallel_backend # type: ignore[import-not-found] -from nilearn.connectome import sym_matrix_to_vec # type: ignore[import-not-found] +from joblib import parallel_backend +from nilearn.connectome import sym_matrix_to_vec from numpy.typing import NDArray -from sklearn.base import BaseEstimator, TransformerMixin # type: ignore[import-not-found] -from sklearn.decomposition import PCA # type: ignore[import-not-found] -from sklearn.impute import SimpleImputer # type: ignore[import-not-found] -from sklearn.linear_model import LogisticRegression, Ridge # type: ignore[import-not-found] -from sklearn.model_selection import StratifiedShuffleSplit, cross_validate # type: ignore[import-not-found] -from sklearn.pipeline import Pipeline # type: ignore[import-not-found] -from sklearn.preprocessing import LabelEncoder, StandardScaler # type: ignore[import-not-found] +from sklearn.base import BaseEstimator, TransformerMixin +from sklearn.decomposition import PCA +from sklearn.impute import SimpleImputer +from sklearn.linear_model import LogisticRegression, Ridge +from sklearn.model_selection import StratifiedShuffleSplit, cross_validate +from sklearn.pipeline import Pipeline +from sklearn.preprocessing import LabelEncoder, StandardScaler if TYPE_CHECKING: from ..base import ConnectivityMatrix diff --git a/wonkyconn/features/calculate_gradients_correlation.py b/wonkyconn/features/calculate_gradients_correlation.py index f8bb087e..df3449f4 100644 --- a/wonkyconn/features/calculate_gradients_correlation.py +++ b/wonkyconn/features/calculate_gradients_correlation.py @@ -4,10 +4,10 @@ import nibabel as nib import numpy as np -from brainspace.gradient import GradientMaps # type: ignore[import-not-found] -from nilearn import image # type: ignore[import-not-found] -from nilearn.connectome import sym_matrix_to_vec, vec_to_sym_matrix # type: ignore[import-not-found] -from nilearn.maskers import NiftiLabelsMasker # type: ignore[import-not-found] +from brainspace.gradient import GradientMaps +from nilearn import image +from nilearn.connectome import sym_matrix_to_vec, vec_to_sym_matrix +from nilearn.maskers import NiftiLabelsMasker from scipy import stats from ..base import ConnectivityMatrix diff --git a/wonkyconn/features/network.py b/wonkyconn/features/network.py index 2a2c895d..7d92121d 100644 --- a/wonkyconn/features/network.py +++ b/wonkyconn/features/network.py @@ -82,7 +82,7 @@ def single_subject_within_network_connectivity( return ( subj_average_connectivity_within_network, subj_variance_connectivity_within_network, - np.nanmean(subj_corr_with_network), + np.nanmean(subj_corr_with_network), # pyright: ignore[reportReturnType] ) diff --git a/wonkyconn/run.py b/wonkyconn/run.py index 154123f8..5b9864d8 100644 --- a/wonkyconn/run.py +++ b/wonkyconn/run.py @@ -67,7 +67,7 @@ def global_parser(exit_on_error: bool = True) -> argparse.ArgumentParser: help="Disable sex and age prediction to reduce runtime.", ) parser.add_argument( - "--site_correction", + "--site-correction", required=False, action="store_true", default=False, diff --git a/wonkyconn/textual_app.py b/wonkyconn/textual_app.py index 0e149e50..6f742fcc 100644 --- a/wonkyconn/textual_app.py +++ b/wonkyconn/textual_app.py @@ -3,10 +3,10 @@ from pathlib import Path from typing import Iterable -from textual.app import App, ComposeResult # type: ignore[import-not-found] -from textual.containers import Center, Container, Horizontal, Vertical # type: ignore[import-not-found] -from textual.events import DescendantFocus # type: ignore[import-not-found] -from textual.widgets import ( # type: ignore[import-not-found] +from textual.app import App, ComposeResult +from textual.containers import Center, Container, Horizontal, Vertical +from textual.events import DescendantFocus +from textual.widgets import ( Button, Checkbox, DirectoryTree, @@ -17,9 +17,9 @@ Select, Static, ) -from textual.widgets._directory_tree import DirEntry # type: ignore[import-not-found] -from textual.widgets._select import NoSelection # type: ignore[import-not-found] -from textual.widgets._tree import Tree # type: ignore[import-not-found] +from textual.widgets._directory_tree import DirEntry +from textual.widgets._select import NoSelection +from textual.widgets._tree import Tree from .config import WonkyConnConfig From c5c593eaef10b8dc1f54e026a9a3a967402254a8 Mon Sep 17 00:00:00 2001 From: Lea Waller Date: Sun, 2 Aug 2026 10:55:38 +0200 Subject: [PATCH 12/33] Add stratification by site to CV sampler Also simplify SiteRegressor by using pandas DataFrames and add tests to ensure that site effects are removed and do not confound prediction after regression --- .pre-commit-config.yaml | 8 +- wonkyconn/features/age_sex_prediction.py | 90 ++++++------ wonkyconn/features/tests/test_prediction.py | 146 ++++++++++++++++++++ 3 files changed, 195 insertions(+), 49 deletions(-) create mode 100644 wonkyconn/features/tests/test_prediction.py diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 877304a1..93eead49 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -20,7 +20,13 @@ repos: rev: v1.19.1 hooks: - id: mypy - additional_dependencies: [pandas-stubs, types-tqdm, types-setuptools, types-Jinja2, textual] + additional_dependencies: + - pandas-stubs + - pytest + - types-tqdm + - types-setuptools + - types-Jinja2 + - textual args: [--config-file=pyproject.toml] - repo: https://github.com/codespell-project/codespell rev: v2.4.2 diff --git a/wonkyconn/features/age_sex_prediction.py b/wonkyconn/features/age_sex_prediction.py index 324909d6..11afcf22 100644 --- a/wonkyconn/features/age_sex_prediction.py +++ b/wonkyconn/features/age_sex_prediction.py @@ -1,5 +1,6 @@ from __future__ import annotations +from dataclasses import dataclass, field from typing import TYPE_CHECKING, Dict, List import numpy as np @@ -10,7 +11,7 @@ from sklearn.base import BaseEstimator, TransformerMixin from sklearn.decomposition import PCA from sklearn.impute import SimpleImputer -from sklearn.linear_model import LogisticRegression, Ridge +from sklearn.linear_model import LinearRegression, LogisticRegression, Ridge from sklearn.model_selection import StratifiedShuffleSplit, cross_validate from sklearn.pipeline import Pipeline from sklearn.preprocessing import LabelEncoder, StandardScaler @@ -19,19 +20,25 @@ from ..base import ConnectivityMatrix +@dataclass class SiteRegressor(BaseEstimator, TransformerMixin): - """ - Regress out site effects from connectivity features. - """ + sites: NDArray[np.str_] + + model: LinearRegression = field(default_factory=LinearRegression) - def __init__(self: SiteRegressor, n_connectivity_features: int) -> None: + def _get_dummies(self: SiteRegressor, X: pd.DataFrame) -> pd.DataFrame: # noqa: N803 """ + Convert site labels to dummy variables. + Args: - n_connectivity_features: Number of connectivity features in the input data. + X: A DataFrame containing the site labels. + + Returns: + A DataFrame with dummy variables for each site. """ - self.n_connectivity_features = n_connectivity_features + return pd.get_dummies(self.sites[X.index], drop_first=False, dtype=np.float32) - def fit(self: SiteRegressor, connectivity_data_site: NDArray[np.float32], y: None = None) -> SiteRegressor: + def fit(self: SiteRegressor, X: pd.DataFrame, y: pd.DataFrame | None = None) -> SiteRegressor: # noqa: N803 """ Fit the site regressor to the connectivity data. @@ -43,18 +50,17 @@ def fit(self: SiteRegressor, connectivity_data_site: NDArray[np.float32], y: Non Returns: self: Returns the instance itself. """ - connectivity_feat = connectivity_data_site[:, : self.n_connectivity_features] - site = connectivity_data_site[:, self.n_connectivity_features :] + fit_sites = np.asarray(self.sites[X.index]) + if np.unique(fit_sites).size < 2: + raise ValueError("SiteRegressor requires at least two sites in the training data.") - # Add intercept - design = np.column_stack([np.ones(len(site)), site]) + y = X # Estimate coefficients on the training fold only - self.beta_ = np.linalg.lstsq(design, connectivity_feat, rcond=None)[0] - + self.model.fit(self._get_dummies(X), y) return self - def transform(self: SiteRegressor, connectivity_data_site: NDArray[np.float32]) -> NDArray[np.float64]: + def transform(self: SiteRegressor, X: pd.DataFrame) -> pd.DataFrame: # noqa: N803 """ Transform the connectivity data by regressing out site effects. Args: @@ -64,15 +70,7 @@ def transform(self: SiteRegressor, connectivity_data_site: NDArray[np.float32]) Returns: A 2D array of connectivity features with site effects regressed out. """ - - connectivity_feat = connectivity_data_site[:, : self.n_connectivity_features] - site = connectivity_data_site[:, self.n_connectivity_features :] - - # Add intercept - design = np.column_stack([np.ones(len(site)), site]) - - # Remove only the site contribution (keep the intercept) - return connectivity_feat - design[:, 1:] @ self.beta_[1:] + return X - self.model.predict(self._get_dummies(X)) def training_pipeline( @@ -100,42 +98,37 @@ def training_pipeline( Returns: pd.DataFrame: Statistics (mean, 95% CI) of the scores obtained. """ - connectivity_data = np.asarray(connectivity_data, dtype=np.float32, order="C") + connectivity_data_frame = pd.DataFrame(connectivity_data, dtype=np.float32) - # Transform site labels into dummy variables and concatenate with connectivity data - if sites is not None: - site = pd.get_dummies( - sites, - drop_first=True, - dtype=np.float32, - ).to_numpy() + if task_type == "classification": + y_train = pd.Series(LabelEncoder().fit_transform(target_labels)) # pyright: ignore[reportArgumentType, reportCallIssue] + estimator = LogisticRegression(max_iter=5000, solver="saga", l1_ratio=0.0, random_state=random_state) - connectivity_data_site = np.concatenate([connectivity_data, site], axis=1) - else: - connectivity_data_site = connectivity_data + bins = y_train - if task_type == "classification": - y_train = LabelEncoder().fit_transform(target_labels) - estimator = LogisticRegression(max_iter=5000, solver="saga", penalty="l2", n_jobs=n_jobs, random_state=random_state) - cv_strategy = StratifiedShuffleSplit(n_splits=n_splits, test_size=0.2, random_state=random_state) scoring_metrics = {"accuracy": "accuracy", "roc_auc": "roc_auc"} else: - y_train = np.asarray(target_labels) + y_train = pd.Series(target_labels) estimator = Ridge(alpha=1.0) bins = pd.qcut(y_train, q=5, labels=False, duplicates="drop") - cv_strategy = StratifiedShuffleSplit(n_splits=n_splits, test_size=0.2, random_state=random_state) - splits = list(cv_strategy.split(np.zeros_like(bins), bins)) - cv_strategy = splits scoring_metrics = {"mae": "neg_mean_absolute_error", "r2": "r2"} - steps = [ - ("imputer", SimpleImputer(strategy="median")), + if sites is not None: + data_frame = pd.DataFrame({"site": sites, "bins": bins}) + # Get unique row indices as combined bins + bins = data_frame.groupby(data_frame.columns.tolist(), sort=False).ngroup() + + cv_strategy = StratifiedShuffleSplit(n_splits=n_splits, test_size=0.2, random_state=random_state) + splits = list(cv_strategy.split(np.zeros(len(bins)), bins)) + + steps: list[tuple[str, BaseEstimator]] = [ + ("imputer", SimpleImputer(strategy="median").set_output(transform="pandas")), ] if sites is not None: - steps.append(("site_regression", SiteRegressor(connectivity_data.shape[1]))) + steps.append(("site_regression", SiteRegressor(sites))) steps.extend( [ @@ -156,11 +149,12 @@ def training_pipeline( with parallel_backend("threading", n_jobs=n_jobs): cv_results = cross_validate( pipe, - connectivity_data_site, + connectivity_data_frame, y_train, - cv=cv_strategy, + cv=splits, scoring=scoring_metrics, n_jobs=n_jobs, + error_score="raise", ) scores_df = pd.DataFrame({k.replace("test_", ""): v for k, v in cv_results.items() if k.startswith("test_")}) diff --git a/wonkyconn/features/tests/test_prediction.py b/wonkyconn/features/tests/test_prediction.py new file mode 100644 index 00000000..85c9c50b --- /dev/null +++ b/wonkyconn/features/tests/test_prediction.py @@ -0,0 +1,146 @@ +from dataclasses import dataclass +from typing import Literal + +import numpy as np +import pandas as pd +import pytest +from numpy.typing import NDArray + +from wonkyconn.features.age_sex_prediction import training_pipeline + +subject_count = 128 +feature_count = 4_096 +random_state = 1 + +Task = Literal["classification", "regression"] + +# Alternating site assignment reused both for confounded data generation and for +# the ``sites`` parametrization below, so the two stay perfectly aligned. +site_labels = np.array(["site-a", "site-b"] * (subject_count // 2)) + + +@dataclass(frozen=True) +class Config: + task: Task + is_predictable: bool = False + # When True, the site drives both the features and the labels, so the target + # is only predictable while the site effect is present in the features. + is_site_confounded: bool = False + + +@dataclass(frozen=True) +class Dataset: + config: Config + connectivity_data: NDArray[np.float32] + target_labels: NDArray[np.float64] | NDArray[np.str_] + + +@pytest.fixture( + params=[ + pytest.param(Config(task="classification", is_predictable=False), id="binary-random-labels"), + pytest.param(Config(task="classification", is_predictable=True), id="binary-predictable-labels"), + pytest.param(Config(task="regression", is_predictable=False), id="continuous-random-labels"), + pytest.param(Config(task="regression", is_predictable=True), id="continuous-predictable-labels"), + pytest.param(Config(task="classification", is_site_confounded=True), id="binary-site-confounded-labels"), + pytest.param(Config(task="regression", is_site_confounded=True), id="continuous-site-confounded-labels"), + ], +) +def dataset(request: pytest.FixtureRequest) -> Dataset: + config: Config = request.param + + rng = np.random.default_rng(random_state) + + labels: NDArray[np.float64] | NDArray[np.str_] + if config.is_site_confounded: + # Site drives the features (via ``site_code``) and the labels together, + # so the target is recoverable only until the site effect is removed. + site_code = np.where(site_labels == "site-a", 1.0, -1.0) + loadings = rng.normal(size=feature_count) + connectivity_data = ( + np.outer(site_code, loadings) + rng.normal(scale=0.3, size=(subject_count, feature_count)) + ).astype(np.float32) + match config.task: + case "classification": + labels = np.where(site_labels == "site-a", "male", "female") + case "regression": + labels = np.where(site_labels == "site-a", 30.0, 60.0) + rng.normal(scale=2.0, size=subject_count) + else: + latent = rng.normal(size=subject_count) + loadings = rng.normal(size=feature_count) + connectivity_data = (np.outer(latent, loadings) + rng.normal(scale=0.3, size=(subject_count, feature_count))).astype( + np.float32 + ) + match config.task: + case "classification": + if config.is_predictable: + score = latent + rng.normal(scale=0.5, size=subject_count) + labels = np.where(score > np.median(score), "male", "female") + else: + labels = np.array(["female", "male"] * (subject_count // 2)) + rng.shuffle(labels) + case "regression": + noise = rng.normal(scale=0.5, size=subject_count) + raw = latent + noise if config.is_predictable else rng.normal(size=subject_count) + # Rescale to a plausible age range. + labels = 18.0 + (raw - raw.min()) / (raw.max() - raw.min()) * 62.0 + + return Dataset( + config=config, + connectivity_data=connectivity_data, + target_labels=labels, + ) + + +@pytest.mark.parametrize( + "sites", + [ + pytest.param(None, id="without-sites"), + pytest.param(site_labels, id="with-sites"), + ], +) +def test_training_pipeline( + dataset: Dataset, + sites: NDArray[np.str_] | None, +) -> None: + if dataset.config.task == "classification": + expected_metrics = frozenset({"accuracy", "roc_auc"}) + primary_metric = "roc_auc" + learnable_threshold = 0.8 + chance_ceiling = 0.65 + else: + expected_metrics = frozenset({"mae", "r2"}) + primary_metric = "r2" + learnable_threshold = 0.3 + chance_ceiling = 0.2 + + summary = training_pipeline( + dataset.connectivity_data, + dataset.target_labels, + task_type=dataset.config.task, + n_splits=3, + n_pca=5, + n_jobs=1, + random_state=42, + sites=sites, + ) + + assert isinstance(summary, pd.DataFrame) + assert set(summary.index) == expected_metrics + assert list(summary.columns) == ["mean", "ci_lower", "ci_upper"] + assert np.isfinite(summary.to_numpy()).all() + assert (summary["ci_lower"] <= summary["mean"]).all() + assert (summary["mean"] <= summary["ci_upper"]).all() + + score = float(summary.loc[primary_metric, "mean"]) # type: ignore[arg-type] + + if dataset.config.is_site_confounded: + if sites is None: + # Without correction the site confound is fully exploitable. + assert score > learnable_threshold + else: + # Site-effect correction must strip the confound, dropping the + # otherwise-perfect prediction back to chance level. + assert score < chance_ceiling + elif dataset.config.is_predictable: + # Labels that carry signal should be recovered above the chance baseline. + assert score > learnable_threshold From 39f306e8492a02042de4d35fd5da6e235b716cfc Mon Sep 17 00:00:00 2001 From: Lea Waller Date: Sun, 2 Aug 2026 11:17:09 +0200 Subject: [PATCH 13/33] Fix warnings that come up during unit tests - The stochastic gradient descent (saga) solver for logistic regression fails to converge on small datasets, which is a core use case for wonkyconn. We switch to lbfgs for greater performance on small data and better convergence - Replace fit_transform call --- wonkyconn/atlas.py | 3 ++- wonkyconn/features/age_sex_prediction.py | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/wonkyconn/atlas.py b/wonkyconn/atlas.py index 3fa8cc86..6131d4a1 100644 --- a/wonkyconn/atlas.py +++ b/wonkyconn/atlas.py @@ -142,7 +142,8 @@ def get_yeo7_membership(self) -> pd.DataFrame: for n in network_labels: cur_region = math_img(f"img=={n}", img=yeo7_nii) masker = NiftiMasker(cur_region) - atlas_parcel_in_network = masker.fit_transform(self.image) + masker.fit() + atlas_parcel_in_network = masker.transform(self.image) atlas_parcel_in_network = np.unique(atlas_parcel_in_network)[1:] region_membership.loc[atlas_parcel_in_network, f"yeo7-{int(n)}"] = 1 return region_membership diff --git a/wonkyconn/features/age_sex_prediction.py b/wonkyconn/features/age_sex_prediction.py index 11afcf22..98a2ae2c 100644 --- a/wonkyconn/features/age_sex_prediction.py +++ b/wonkyconn/features/age_sex_prediction.py @@ -102,7 +102,7 @@ def training_pipeline( if task_type == "classification": y_train = pd.Series(LabelEncoder().fit_transform(target_labels)) # pyright: ignore[reportArgumentType, reportCallIssue] - estimator = LogisticRegression(max_iter=5000, solver="saga", l1_ratio=0.0, random_state=random_state) + estimator = LogisticRegression(max_iter=5000, solver="lbfgs", random_state=random_state) bins = y_train From dbfdac5cf72e3da8c40087473b714ac3318461a1 Mon Sep 17 00:00:00 2001 From: Lea Waller Date: Sun, 2 Aug 2026 11:29:33 +0200 Subject: [PATCH 14/33] Update test code for readability --- wonkyconn/features/tests/test_prediction.py | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/wonkyconn/features/tests/test_prediction.py b/wonkyconn/features/tests/test_prediction.py index 85c9c50b..91bc76a3 100644 --- a/wonkyconn/features/tests/test_prediction.py +++ b/wonkyconn/features/tests/test_prediction.py @@ -70,17 +70,15 @@ def dataset(request: pytest.FixtureRequest) -> Dataset: connectivity_data = (np.outer(latent, loadings) + rng.normal(scale=0.3, size=(subject_count, feature_count))).astype( np.float32 ) + noise = rng.normal(scale=0.5, size=subject_count) match config.task: case "classification": - if config.is_predictable: - score = latent + rng.normal(scale=0.5, size=subject_count) - labels = np.where(score > np.median(score), "male", "female") - else: - labels = np.array(["female", "male"] * (subject_count // 2)) + score = latent + noise + labels = np.where(score > np.median(score), "male", "female") + if not config.is_predictable: rng.shuffle(labels) case "regression": - noise = rng.normal(scale=0.5, size=subject_count) - raw = latent + noise if config.is_predictable else rng.normal(size=subject_count) + raw = (latent + noise) if config.is_predictable else (noise + noise) # Rescale to a plausible age range. labels = 18.0 + (raw - raw.min()) / (raw.max() - raw.min()) * 62.0 @@ -144,3 +142,6 @@ def test_training_pipeline( elif dataset.config.is_predictable: # Labels that carry signal should be recovered above the chance baseline. assert score > learnable_threshold + else: + # Random labels should not be recoverable above chance level. + assert score < chance_ceiling From 468177e14af1989ce0dff09ba77b509c9dfe01ac Mon Sep 17 00:00:00 2001 From: Lea Waller Date: Sun, 2 Aug 2026 11:38:03 +0200 Subject: [PATCH 15/33] [DATALAD] Recorded changes --- data/halfpipe/participants.tsv | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/data/halfpipe/participants.tsv b/data/halfpipe/participants.tsv index 9d371793..9f8d3306 120000 --- a/data/halfpipe/participants.tsv +++ b/data/halfpipe/participants.tsv @@ -1 +1 @@ -../../.git/annex/objects/gG/5m/MD5E-s16387--b0d7908c25384a3ecd7c317238748596.tsv/MD5E-s16387--b0d7908c25384a3ecd7c317238748596.tsv \ No newline at end of file +../../.git/annex/objects/Fj/K7/MD5E-s16372--3c23f937d8662e3bbd2f588f3d992e8b.tsv/MD5E-s16372--3c23f937d8662e3bbd2f588f3d992e8b.tsv \ No newline at end of file From e1eb6f3b44aea9b3f78f484873ca927e6075ca67 Mon Sep 17 00:00:00 2001 From: Lea Waller Date: Sun, 2 Aug 2026 12:17:40 +0200 Subject: [PATCH 16/33] Ensure full code coverage --- codecov.yml | 15 ++++++---- wonkyconn/features/tests/test_prediction.py | 31 ++++++++++++++++++++- 2 files changed, 39 insertions(+), 7 deletions(-) diff --git a/codecov.yml b/codecov.yml index 7b9ef93d..e4af7144 100644 --- a/codecov.yml +++ b/codecov.yml @@ -14,13 +14,16 @@ # # noisy red coverage status on github PRs. # target: auto # threshold: 1% -comment: # this is a top-level key +comment: # this is a top-level key layout: reach, diff, flags, files behavior: default - require_changes: false # if true: only post the comment if coverage changes - require_base: no # [yes :: must have a base report to post] - require_head: yes # [yes :: must have a head report to post] + require_changes: false # if true: only post the comment if coverage changes + require_base: false # if true: must have a base report to post + require_head: false # if true: must have a head report to post ignore: -- '*/tests/' # ignore folders related to testing -- '*/data/' + - "**/tests/**" + - "**/test_*.py" + - "**/conftest.py" + - "**/data/**" + - "**/setup.py" diff --git a/wonkyconn/features/tests/test_prediction.py b/wonkyconn/features/tests/test_prediction.py index 91bc76a3..aa89718e 100644 --- a/wonkyconn/features/tests/test_prediction.py +++ b/wonkyconn/features/tests/test_prediction.py @@ -1,4 +1,5 @@ from dataclasses import dataclass +from pathlib import Path from typing import Literal import numpy as np @@ -6,7 +7,8 @@ import pytest from numpy.typing import NDArray -from wonkyconn.features.age_sex_prediction import training_pipeline +from wonkyconn.base import ConnectivityMatrix +from wonkyconn.features.age_sex_prediction import SiteRegressor, age_sex_scores, training_pipeline subject_count = 128 feature_count = 4_096 @@ -145,3 +147,30 @@ def test_training_pipeline( else: # Random labels should not be recoverable above chance level. assert score < chance_ceiling + + +def test_age_sex_scores(tmp_path: Path) -> None: + rng = np.random.default_rng(random_state) + + region_count = np.sqrt(feature_count).astype(int).item() + + matrices = [] + for i in range(subject_count): + square = rng.normal(size=(region_count, region_count)).astype(np.float32) + path = tmp_path / f"sub-{i}.tsv" + np.savetxt(path, square + square.T, delimiter="\t") + matrices.append(ConnectivityMatrix(path, metadata=dict())) + + ages = rng.uniform(18.0, 80.0, subject_count) + genders = np.array(["m", "f"] * (subject_count // 2)) + sites = np.array(["site-a", "site-b"] * (subject_count // 2)) + + scores = age_sex_scores(matrices, ages, genders, sites, n_splits=3, n_pca=5, n_jobs=1) + assert len(scores) == 8 + assert all(np.isfinite(value) for value in scores.values()) + + +def test_site_regressor_requires_two_sites() -> None: + regressor = SiteRegressor(np.array(["site-a", "site-a"])) + with pytest.raises(ValueError, match="at least two sites"): + regressor.fit(pd.DataFrame(np.zeros((2, 3), dtype=np.float32))) From bea3359a618fc36c7c5d1306e0a5142ac2d16a0a Mon Sep 17 00:00:00 2001 From: Lea Waller Date: Sun, 2 Aug 2026 12:28:10 +0200 Subject: [PATCH 17/33] Upload code coverage from unittest run too --- .github/workflows/test.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 60601758..8ec1b24a 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -65,7 +65,7 @@ jobs: environments: test frozen: true - run: pixi run unittest - - run: pixi run smoketestlight + - run: pixi run smoketestlight --cov-append - uses: codecov/codecov-action@v6 if: ${{ always() }} with: From 24298ca1a34945e41520655c24963745596069ce Mon Sep 17 00:00:00 2001 From: Lea Waller Date: Sun, 2 Aug 2026 13:14:24 +0200 Subject: [PATCH 18/33] Fix or ignore all pyright errors --- .../features/calculate_degrees_of_freedom.py | 2 +- .../calculate_gradients_correlation.py | 19 ++++++++++--------- wonkyconn/features/distance_dependence.py | 4 ++-- .../features/quality_control_connectivity.py | 4 ++-- wonkyconn/features/tests/test_network.py | 6 +++--- wonkyconn/tests/test_atlas.py | 15 +++++---------- wonkyconn/tests/test_cli.py | 10 +++++----- wonkyconn/tests/test_correlation.py | 2 +- wonkyconn/tests/test_gradients_correlation.py | 3 ++- 9 files changed, 31 insertions(+), 34 deletions(-) diff --git a/wonkyconn/features/calculate_degrees_of_freedom.py b/wonkyconn/features/calculate_degrees_of_freedom.py index cfb9772f..a023f0aa 100644 --- a/wonkyconn/features/calculate_degrees_of_freedom.py +++ b/wonkyconn/features/calculate_degrees_of_freedom.py @@ -98,5 +98,5 @@ def calculate_for_key( else: raise ValueError(f"Unexpected value for `{keys}`: {value}") - percentages = pd.Series(proportions) * 100 + percentages: pd.Series = pd.Series(proportions) * 100 return percentages.mean() diff --git a/wonkyconn/features/calculate_gradients_correlation.py b/wonkyconn/features/calculate_gradients_correlation.py index df3449f4..ec30162d 100644 --- a/wonkyconn/features/calculate_gradients_correlation.py +++ b/wonkyconn/features/calculate_gradients_correlation.py @@ -5,6 +5,7 @@ import nibabel as nib import numpy as np from brainspace.gradient import GradientMaps +from nibabel.nifti1 import Nifti1Image from nilearn import image from nilearn.connectome import sym_matrix_to_vec, vec_to_sym_matrix from nilearn.maskers import NiftiLabelsMasker @@ -38,7 +39,7 @@ def remove_nan_from_matrix(matrix: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: return conn_clean, kept_idx -def remove_nan_roi_atlas(atlas: nib.Nifti1Image, kept_idx: np.ndarray) -> nib.Nifti1Image: +def remove_nan_roi_atlas(atlas: nib.Nifti1Image, kept_idx: np.ndarray) -> Nifti1Image: """Remove ROIs from an atlas that are not present in a connectivity matrix. Args: @@ -62,10 +63,10 @@ def remove_nan_roi_atlas(atlas: nib.Nifti1Image, kept_idx: np.ndarray) -> nib.Ni keep_mask = np.isin(atlas_data, kept_labels) kept_atlas_data = atlas_data.copy() kept_atlas_data[~keep_mask] = 0 - return nib.Nifti1Image(kept_atlas_data, atlas.affine, atlas.header), kept_labels + return Nifti1Image(kept_atlas_data, atlas.affine, atlas.header) -def overlapping_atlas_with_mask(subject_atlas: nib.Nifti1Image, group_mask: nib.Nifti1Image) -> nib.Nifti1Image: +def overlapping_atlas_with_mask(subject_atlas: Nifti1Image, group_mask: Nifti1Image) -> Nifti1Image: """Create a new atlas containing only regions that overlap with the group gradient mask. Args: @@ -76,9 +77,9 @@ def overlapping_atlas_with_mask(subject_atlas: nib.Nifti1Image, group_mask: nib. nib.Nifti1Image: Atlas with non-overlapping regions zeroed out. """ - mask_gradient_resampled = image.resample_to_img( + mask_gradient_resampled: Nifti1Image = image.resample_to_img( group_mask, subject_atlas, interpolation="nearest", copy_header=True, force_resample=True - ) + ) # pyright: ignore[reportAssignmentType] # Get arrays atlas_data = subject_atlas.get_fdata() @@ -153,7 +154,7 @@ def process_single_matrix( """ matrix = np.asarray(connectivity_matrix, dtype=np.float64) conn_clean, kept_idx = remove_nan_from_matrix(matrix) - atlas_mask_without_nan, _ = remove_nan_roi_atlas(atlas, kept_idx) + atlas_mask_without_nan = remove_nan_roi_atlas(atlas, kept_idx) # filter out labels > 400 atlas_data = atlas_mask_without_nan.get_fdata() @@ -180,7 +181,7 @@ def process_single_matrix( gm = GradientMaps(approach="pca", n_components=5, alignment="procrustes", kernel="normalized_angle") ind_gradient = gm.fit(masked_matrix, reference=group_gradients_np) - return ind_gradient.aligned_, group_gradients_np + return ind_gradient.aligned_, group_gradients_np # pyright: ignore[reportReturnType] def extract_gradients( @@ -201,11 +202,11 @@ def extract_gradients( repo_root = Path(__file__).resolve().parent.parent path_gradients = repo_root / "data" / "gradients" - gradient_mask = nib.load(path_gradients / "gradientmask_cortical.nii.gz") + gradient_mask = nib.nifti1.load(path_gradients / "gradientmask_cortical.nii.gz") # pyright: ignore[reportAttributeAccessIssue] # Load all group gradient templates gradient_files = sorted(glob.glob(str(path_gradients / "templates" / "gradient*_cortical_only.nii.gz"))) - gradient_imgs = [nib.load(fname) for fname in gradient_files] + gradient_imgs = [nib.nifti1.load(fname) for fname in gradient_files] # pyright: ignore[reportAttributeAccessIssue] mean_connectome = group_mean_connectivity(connectivity_matrices) gradient_aligned, template_gradient = process_single_matrix(mean_connectome, atlas, gradient_mask, gradient_imgs) diff --git a/wonkyconn/features/distance_dependence.py b/wonkyconn/features/distance_dependence.py index 03279a99..a1698274 100644 --- a/wonkyconn/features/distance_dependence.py +++ b/wonkyconn/features/distance_dependence.py @@ -19,5 +19,5 @@ def calculate_distance_dependence(qcfc: pd.DataFrame, distance_matrix: npt.NDArr """ i, j = map(np.asarray, zip(*qcfc.index, strict=False)) distance_vector = distance_matrix[i, j] - r, _ = spearmanr(distance_vector, qcfc.correlation, nan_policy="omit") - return np.abs(r) + statistic, _ = spearmanr(distance_vector, qcfc.correlation, nan_policy="omit") + return np.abs(statistic).item() # pyright: ignore[reportArgumentType, reportCallIssue] diff --git a/wonkyconn/features/quality_control_connectivity.py b/wonkyconn/features/quality_control_connectivity.py index 75498b10..24c46bb5 100644 --- a/wonkyconn/features/quality_control_connectivity.py +++ b/wonkyconn/features/quality_control_connectivity.py @@ -70,7 +70,7 @@ def calculate_qcfc( axis=1, ) - correlation, count = partial_correlation(connectivity_array, metrics, covariates) + correlation, count = partial_correlation(connectivity_array, metrics, covariates) # pyright: ignore[reportCallIssue] p_value = correlation_p_value(correlation, count) qcfc = pd.DataFrame(dict(i=i, j=j, correlation=correlation, p_value=p_value)) @@ -100,7 +100,7 @@ def significant_level(x: "pd.Series[float]", alpha: float = 0.05, correction: st res, _, _, _ = multipletests(x, alpha=alpha, method=correction) else: res = x < alpha - return res + return np.asarray(res) def calculate_qcfc_percentage(qcfc: pd.DataFrame) -> float: diff --git a/wonkyconn/features/tests/test_network.py b/wonkyconn/features/tests/test_network.py index a6e314a8..1efa365b 100644 --- a/wonkyconn/features/tests/test_network.py +++ b/wonkyconn/features/tests/test_network.py @@ -22,9 +22,9 @@ def test_single_subject_within_network_connectivity(data_path: Path) -> None: # / "halfpipe/derivatives/halfpipe/" / "sub-10171/func/task-rest/sub-10171_task-rest_feature-cCompCor_atlas-Schaefer2018Combined_timeseries.json" ) - dl.get(str(dseg_path)) - dl.get(str(relmat_path)) - dl.get(str(metadata_path)) + dl.get(str(dseg_path)) # pyright: ignore[reportAttributeAccessIssue] + dl.get(str(relmat_path)) # pyright: ignore[reportAttributeAccessIssue] + dl.get(str(metadata_path)) # pyright: ignore[reportAttributeAccessIssue] atlas = Atlas.create("Schaefer2018Combined", dseg_path) with metadata_path.open("r") as file: metadata = json.load(file) diff --git a/wonkyconn/tests/test_atlas.py b/wonkyconn/tests/test_atlas.py index cf72ece1..e0c04e14 100644 --- a/wonkyconn/tests/test_atlas.py +++ b/wonkyconn/tests/test_atlas.py @@ -9,10 +9,7 @@ from wonkyconn.atlas import Atlas -def test_dseg_atlas(data_path: Path) -> None: - # atlas_path = data_path / YEO_NETWORK_MAP - # dl.get(str(atlas_path)) - +def test_dseg_atlas() -> None: url = ( "https://raw.githubusercontent.com/ThomasYeoLab/CBIG/master/" "stable_projects/brain_parcellation/Schaefer2018_LocalGlobal/" @@ -29,7 +26,7 @@ def test_dseg_atlas(data_path: Path) -> None: resolution=2, suffix="dseg", extension=".nii.gz", - ) + ) # pyright: ignore[reportCallIssue] assert isinstance(path, Path) atlas = Atlas.create("Schaefer2018400Parcels7Networks", path) @@ -60,10 +57,8 @@ def _get_centroids(path: Path): return centroids -def test_probseg_atlas(data_path: Path) -> None: - # "TODO: @haoting wants to revisit this test, to check if the assertion values make sense" - # atlas_path = data_path / YEO_NETWORK_MAP - # dl.get(str(atlas_path)) +def test_probseg_atlas() -> None: + # TODO: @haoting wants to revisit this test, to check if the assertion values make sense path = get_template( template="MNI152NLin2009cAsym", @@ -72,7 +67,7 @@ def test_probseg_atlas(data_path: Path) -> None: suffix="probseg", resolution=3, # matches “res-03” extension=".nii.gz", - ) + ) # pyright: ignore[reportCallIssue] assert isinstance(path, Path) _centroids = _get_centroids(path) diff --git a/wonkyconn/tests/test_cli.py b/wonkyconn/tests/test_cli.py index 25343b29..99885ceb 100644 --- a/wonkyconn/tests/test_cli.py +++ b/wonkyconn/tests/test_cli.py @@ -130,7 +130,7 @@ def _copy_file(path: Path, new_path: Path, sub: str) -> None: @pytest.mark.heavy_smoke def test_giga_connectome(data_path: Path, tmp_path: Path): data_path = data_path / "giga_connectome" / "connectome_Schaefer20187Networks_dev" - dl.get(str(data_path)) + dl.get(str(data_path)) # pyright: ignore[reportAttributeAccessIssue] bids_dir = tmp_path / "bids" bids_dir.mkdir() @@ -187,7 +187,7 @@ def test_giga_connectome(data_path: Path, tmp_path: Path): @pytest.mark.smoke def test_halfpipe(data_path: Path, tmp_path: Path): bids_dir = data_path / "halfpipe" - dl.get(str(bids_dir)) + dl.get(str(bids_dir)) # pyright: ignore[reportAttributeAccessIssue] index = BIDSIndex() index.put(bids_dir) @@ -198,7 +198,7 @@ def test_halfpipe(data_path: Path, tmp_path: Path): phenotypes_path = bids_dir / "participants.tsv" atlas_path = data_path / "atlases" - dl.get(str(atlas_path)) + dl.get(str(atlas_path)) # pyright: ignore[reportAttributeAccessIssue] atlas_args: list[str] = list() atlas_args.append("--atlas") @@ -237,7 +237,7 @@ def test_halfpipe(data_path: Path, tmp_path: Path): @pytest.mark.heavy_smoke def test_halfpipe_with_full_metrics(data_path: Path, tmp_path: Path): bids_dir = data_path / "halfpipe" - dl.get(str(bids_dir)) + dl.get(str(bids_dir)) # pyright: ignore[reportAttributeAccessIssue] index = BIDSIndex() index.put(bids_dir) @@ -248,7 +248,7 @@ def test_halfpipe_with_full_metrics(data_path: Path, tmp_path: Path): phenotypes_path = bids_dir / "participants.tsv" atlas_path = data_path / "atlases" - dl.get(str(atlas_path)) + dl.get(str(atlas_path)) # pyright: ignore[reportAttributeAccessIssue] atlas_args: list[str] = list() atlas_args.append("--atlas") diff --git a/wonkyconn/tests/test_correlation.py b/wonkyconn/tests/test_correlation.py index 0eb2cec2..0ea7320b 100644 --- a/wonkyconn/tests/test_correlation.py +++ b/wonkyconn/tests/test_correlation.py @@ -15,7 +15,7 @@ def test_correlation() -> None: y = np.random.normal(size=(m,)) cov = np.random.normal(size=(m, 2)) - correlation, count = partial_correlation(x, y, cov) + correlation, count = partial_correlation(x, y, cov) # pyright: ignore[reportCallIssue] p_value = correlation_p_value(correlation, count) assert np.all(np.isfinite(correlation)) diff --git a/wonkyconn/tests/test_gradients_correlation.py b/wonkyconn/tests/test_gradients_correlation.py index 7226d83e..69f3aea9 100644 --- a/wonkyconn/tests/test_gradients_correlation.py +++ b/wonkyconn/tests/test_gradients_correlation.py @@ -45,7 +45,8 @@ def test_gradients(): connectivity_matrices = create_fake_connectivity(n_regions=434, n_subjects=3) random_gradient, group_gradients = calculate_gradients_correlation.extract_gradients( - connectivity_matrices, atlas=nib.load(atlas_file) + connectivity_matrices, + atlas=nib.nifti1.load(atlas_file), # pyright: ignore[reportAttributeAccessIssue] ) # Calculate similarity for random gradient From dcd9fd750c03ee1406730bd60ee5633a896dc6f7 Mon Sep 17 00:00:00 2001 From: Lea Waller Date: Sun, 2 Aug 2026 13:27:58 +0200 Subject: [PATCH 19/33] Parametrize site correction for smoke tests Also add type annotations to tests --- wonkyconn/tests/test_cli.py | 18 ++++++++++++------ wonkyconn/tests/test_gradients_correlation.py | 2 +- 2 files changed, 13 insertions(+), 7 deletions(-) diff --git a/wonkyconn/tests/test_cli.py b/wonkyconn/tests/test_cli.py index 99885ceb..318b99a2 100644 --- a/wonkyconn/tests/test_cli.py +++ b/wonkyconn/tests/test_cli.py @@ -22,7 +22,7 @@ from wonkyconn.workflow import workflow -def test_version(capsys): +def test_version(capsys: pytest.CaptureFixture[str]) -> None: try: main(["-v"]) except SystemExit: @@ -31,7 +31,7 @@ def test_version(capsys): assert __version__ == captured.out.split()[0] -def test_help(capsys): +def test_help(capsys: pytest.CaptureFixture[str]) -> None: try: main(["-h"]) except SystemExit: @@ -40,7 +40,7 @@ def test_help(capsys): assert "Evaluating the residual motion in fMRI connectome and visualize reports" in captured.out -def test_cli_and_textual_namespace_consistency(tmp_path: Path): +def test_cli_and_textual_namespace_consistency(tmp_path: Path) -> None: """Ensure CLI (run.py) and Textual UI produce namespaces with the same attributes. Both interfaces should produce namespaces that workflow() can consume, @@ -128,7 +128,7 @@ def _copy_file(path: Path, new_path: Path, sub: str) -> None: @pytest.mark.heavy_smoke -def test_giga_connectome(data_path: Path, tmp_path: Path): +def test_giga_connectome(data_path: Path, tmp_path: Path) -> None: data_path = data_path / "giga_connectome" / "connectome_Schaefer20187Networks_dev" dl.get(str(data_path)) # pyright: ignore[reportAttributeAccessIssue] @@ -185,7 +185,8 @@ def test_giga_connectome(data_path: Path, tmp_path: Path): @pytest.mark.smoke -def test_halfpipe(data_path: Path, tmp_path: Path): +@pytest.mark.parametrize("site_correction", [False, True], ids=["no-site-correction", "site-correction"]) +def test_halfpipe(data_path: Path, tmp_path: Path, site_correction: bool) -> None: bids_dir = data_path / "halfpipe" dl.get(str(bids_dir)) # pyright: ignore[reportAttributeAccessIssue] @@ -216,6 +217,8 @@ def test_halfpipe(data_path: Path, tmp_path: Path): str(output_dir), "group", ] + if site_correction: + argv.insert(0, "--site-correction") args = parser.parse_args(argv) workflow(args) @@ -235,7 +238,8 @@ def test_halfpipe(data_path: Path, tmp_path: Path): @pytest.mark.heavy_smoke -def test_halfpipe_with_full_metrics(data_path: Path, tmp_path: Path): +@pytest.mark.parametrize("site_correction", [False, True], ids=["no-site-correction", "site-correction"]) +def test_halfpipe_with_full_metrics(data_path: Path, tmp_path: Path, site_correction: bool) -> None: bids_dir = data_path / "halfpipe" dl.get(str(bids_dir)) # pyright: ignore[reportAttributeAccessIssue] @@ -265,6 +269,8 @@ def test_halfpipe_with_full_metrics(data_path: Path, tmp_path: Path): str(output_dir), "group", ] + if site_correction: + argv.insert(0, "--site-correction") args = parser.parse_args(argv) workflow(args) diff --git a/wonkyconn/tests/test_gradients_correlation.py b/wonkyconn/tests/test_gradients_correlation.py index 69f3aea9..b4e840fb 100644 --- a/wonkyconn/tests/test_gradients_correlation.py +++ b/wonkyconn/tests/test_gradients_correlation.py @@ -36,7 +36,7 @@ def create_fake_connectivity( return connectivity_matrices -def test_gradients(): +def test_gradients() -> None: repo_root = Path(__file__).resolve().parent.parent path = repo_root / "data" / "gradients" From 93818986f9fe9f110bf569e6d02cfacda16622de Mon Sep 17 00:00:00 2001 From: Lea Waller Date: Sun, 2 Aug 2026 14:53:12 +0200 Subject: [PATCH 20/33] Silence warnings for missing data --- wonkyconn/features/network.py | 15 ++++++++++----- .../features/quality_control_connectivity.py | 3 ++- 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/wonkyconn/features/network.py b/wonkyconn/features/network.py index 7d92121d..ffa76e1a 100644 --- a/wonkyconn/features/network.py +++ b/wonkyconn/features/network.py @@ -4,6 +4,7 @@ annotations, ) +import warnings from typing import Tuple import numpy as np @@ -57,8 +58,10 @@ def single_subject_within_network_connectivity( ] networks = np.asarray(isolate_parcel) networks[networks == 0] = np.nan - mean_within_network_connection = np.nanmean(networks, axis=1) - std_within_network_connection = np.nanstd(networks, axis=1) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", RuntimeWarning) + mean_within_network_connection = np.nanmean(networks, axis=1) + std_within_network_connection = np.nanstd(networks, axis=1) # Remove rows with NaN values mask = np.array( @@ -68,9 +71,11 @@ def single_subject_within_network_connectivity( ] ) - correlation_with_given_network = np.corrcoef( - seed_based_map[mask], region_membership[f"yeo7-{yeo_network_index}"].to_numpy()[mask] - ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", RuntimeWarning) + correlation_with_given_network = np.corrcoef( + seed_based_map[mask], region_membership[f"yeo7-{yeo_network_index}"].to_numpy()[mask] + ) subj_average_connectivity_within_network_list.append(mean_within_network_connection) subj_variance_connectivity_within_network_list.append(std_within_network_connection) diff --git a/wonkyconn/features/quality_control_connectivity.py b/wonkyconn/features/quality_control_connectivity.py index 24c46bb5..0b875f33 100644 --- a/wonkyconn/features/quality_control_connectivity.py +++ b/wonkyconn/features/quality_control_connectivity.py @@ -70,7 +70,8 @@ def calculate_qcfc( axis=1, ) - correlation, count = partial_correlation(connectivity_array, metrics, covariates) # pyright: ignore[reportCallIssue] + with np.errstate(divide="ignore", invalid="ignore"): + correlation, count = partial_correlation(connectivity_array, metrics, covariates) # pyright: ignore[reportCallIssue] p_value = correlation_p_value(correlation, count) qcfc = pd.DataFrame(dict(i=i, j=j, correlation=correlation, p_value=p_value)) From 253a9a8c3a806e42b7a3481ed3e96eff1acafec2 Mon Sep 17 00:00:00 2001 From: Lea Waller Date: Sun, 2 Aug 2026 15:56:49 +0200 Subject: [PATCH 21/33] Improve readability even more --- wonkyconn/features/tests/test_prediction.py | 112 ++++++++++---------- 1 file changed, 57 insertions(+), 55 deletions(-) diff --git a/wonkyconn/features/tests/test_prediction.py b/wonkyconn/features/tests/test_prediction.py index aa89718e..67608f7d 100644 --- a/wonkyconn/features/tests/test_prediction.py +++ b/wonkyconn/features/tests/test_prediction.py @@ -1,6 +1,6 @@ from dataclasses import dataclass from pathlib import Path -from typing import Literal +from typing import Callable, Literal import numpy as np import pandas as pd @@ -10,92 +10,94 @@ from wonkyconn.base import ConnectivityMatrix from wonkyconn.features.age_sex_prediction import SiteRegressor, age_sex_scores, training_pipeline -subject_count = 128 +subject_count = 192 feature_count = 4_096 random_state = 1 -Task = Literal["classification", "regression"] - -# Alternating site assignment reused both for confounded data generation and for -# the ``sites`` parametrization below, so the two stay perfectly aligned. site_labels = np.array(["site-a", "site-b"] * (subject_count // 2)) @dataclass(frozen=True) class Config: - task: Task - is_predictable: bool = False - # When True, the site drives both the features and the labels, so the target - # is only predictable while the site effect is present in the features. - is_site_confounded: bool = False + task: Literal["classification", "regression"] + # How the target relates to the connectivity features: + # "random" - labels are independent noise (nothing to learn) + # "predictable" - a latent factor drives both features and labels + # "site-confounded" - the site drives both features and labels + signal: Literal["random", "predictable", "site-confounded"] @dataclass(frozen=True) class Dataset: config: Config connectivity_data: NDArray[np.float32] - target_labels: NDArray[np.float64] | NDArray[np.str_] + labels: NDArray[np.float64] | NDArray[np.str_] + + +def sex(score: NDArray[np.float64]) -> NDArray[np.str_]: + """Convert a continuous score to a binary sex label based on the mean.""" + return np.where(score > np.mean(score), "male", "female") + + +def age(score: NDArray[np.float64]) -> NDArray[np.float64]: + """Convert a continuous score to an age label based on a linear transformation.""" + return 18.0 + (score - score.min()) / (score.max() - score.min()) * 62.0 @pytest.fixture( params=[ - pytest.param(Config(task="classification", is_predictable=False), id="binary-random-labels"), - pytest.param(Config(task="classification", is_predictable=True), id="binary-predictable-labels"), - pytest.param(Config(task="regression", is_predictable=False), id="continuous-random-labels"), - pytest.param(Config(task="regression", is_predictable=True), id="continuous-predictable-labels"), - pytest.param(Config(task="classification", is_site_confounded=True), id="binary-site-confounded-labels"), - pytest.param(Config(task="regression", is_site_confounded=True), id="continuous-site-confounded-labels"), + pytest.param(Config(task="classification", signal="random"), id="classification-random"), + pytest.param(Config(task="classification", signal="predictable"), id="classification-predictable"), + pytest.param(Config(task="classification", signal="site-confounded"), id="classification-site-confounded"), + pytest.param(Config(task="regression", signal="random"), id="regression-random"), + pytest.param(Config(task="regression", signal="predictable"), id="regression-predictable"), + pytest.param(Config(task="regression", signal="site-confounded"), id="regression-site-confounded"), ], ) def dataset(request: pytest.FixtureRequest) -> Dataset: config: Config = request.param - rng = np.random.default_rng(random_state) + loadings = rng.normal(size=feature_count) + + connectivity_data: NDArray[np.float32] labels: NDArray[np.float64] | NDArray[np.str_] - if config.is_site_confounded: - # Site drives the features (via ``site_code``) and the labels together, - # so the target is recoverable only until the site effect is removed. - site_code = np.where(site_labels == "site-a", 1.0, -1.0) - loadings = rng.normal(size=feature_count) - connectivity_data = ( - np.outer(site_code, loadings) + rng.normal(scale=0.3, size=(subject_count, feature_count)) - ).astype(np.float32) - match config.task: - case "classification": - labels = np.where(site_labels == "site-a", "male", "female") - case "regression": - labels = np.where(site_labels == "site-a", 30.0, 60.0) + rng.normal(scale=2.0, size=subject_count) + + if config.signal == "site-confounded": + latent = np.where(site_labels == "site-a", 1.0, -1.0) else: latent = rng.normal(size=subject_count) - loadings = rng.normal(size=feature_count) - connectivity_data = (np.outer(latent, loadings) + rng.normal(scale=0.3, size=(subject_count, feature_count))).astype( - np.float32 - ) - noise = rng.normal(scale=0.5, size=subject_count) - match config.task: - case "classification": - score = latent + noise - labels = np.where(score > np.median(score), "male", "female") - if not config.is_predictable: - rng.shuffle(labels) - case "regression": - raw = (latent + noise) if config.is_predictable else (noise + noise) - # Rescale to a plausible age range. - labels = 18.0 + (raw - raw.min()) / (raw.max() - raw.min()) * 62.0 - + noise = rng.normal(scale=0.5, size=subject_count) + + func: Callable[[NDArray[np.float64]], NDArray[np.float64] | NDArray[np.str_]] + match config.task: + case "classification": + func = sex + case "regression": + func = age + + match config.signal: + case "random": + labels = func(rng.normal(size=subject_count)) + case "predictable": + labels = func(latent + noise) + case "site-confounded": + labels = func(latent) + + noise = rng.normal(scale=0.3, size=(subject_count, feature_count)) + connectivity_data = (np.outer(latent, loadings) + noise).astype(np.float32) return Dataset( config=config, connectivity_data=connectivity_data, - target_labels=labels, + labels=labels, ) @pytest.mark.parametrize( "sites", [ - pytest.param(None, id="without-sites"), - pytest.param(site_labels, id="with-sites"), + pytest.param(None, id="without-site-correction"), + pytest.param(site_labels, id="with-site-correction"), ], ) def test_training_pipeline( @@ -115,12 +117,12 @@ def test_training_pipeline( summary = training_pipeline( dataset.connectivity_data, - dataset.target_labels, + dataset.labels, task_type=dataset.config.task, n_splits=3, n_pca=5, n_jobs=1, - random_state=42, + random_state=random_state, sites=sites, ) @@ -133,7 +135,7 @@ def test_training_pipeline( score = float(summary.loc[primary_metric, "mean"]) # type: ignore[arg-type] - if dataset.config.is_site_confounded: + if dataset.config.signal == "site-confounded": if sites is None: # Without correction the site confound is fully exploitable. assert score > learnable_threshold @@ -141,7 +143,7 @@ def test_training_pipeline( # Site-effect correction must strip the confound, dropping the # otherwise-perfect prediction back to chance level. assert score < chance_ceiling - elif dataset.config.is_predictable: + elif dataset.config.signal == "predictable": # Labels that carry signal should be recovered above the chance baseline. assert score > learnable_threshold else: From dedee88c70e29ca0045273dd0265040ef41375ee Mon Sep 17 00:00:00 2001 From: Lea Waller Date: Sun, 2 Aug 2026 17:40:04 +0200 Subject: [PATCH 22/33] Fix or ignore remaining warnings --- wonkyconn/features/age_sex_prediction.py | 4 +++- .../calculate_gradients_correlation.py | 6 ++++- wonkyconn/tests/test_textual_app.py | 7 ++++++ wonkyconn/tests/test_workflow.py | 23 +++++++++++++++++++ 4 files changed, 38 insertions(+), 2 deletions(-) create mode 100644 wonkyconn/tests/test_textual_app.py create mode 100644 wonkyconn/tests/test_workflow.py diff --git a/wonkyconn/features/age_sex_prediction.py b/wonkyconn/features/age_sex_prediction.py index 98a2ae2c..297df741 100644 --- a/wonkyconn/features/age_sex_prediction.py +++ b/wonkyconn/features/age_sex_prediction.py @@ -124,7 +124,9 @@ def training_pipeline( splits = list(cv_strategy.split(np.zeros(len(bins)), bins)) steps: list[tuple[str, BaseEstimator]] = [ - ("imputer", SimpleImputer(strategy="median").set_output(transform="pandas")), + # keep_empty_features=True avoids sklearn's "Skipping features without any + # observed values" warning. + ("imputer", SimpleImputer(strategy="median", keep_empty_features=True).set_output(transform="pandas")), ] if sites is not None: diff --git a/wonkyconn/features/calculate_gradients_correlation.py b/wonkyconn/features/calculate_gradients_correlation.py index ec30162d..dfa20c67 100644 --- a/wonkyconn/features/calculate_gradients_correlation.py +++ b/wonkyconn/features/calculate_gradients_correlation.py @@ -1,4 +1,5 @@ import glob +import warnings from pathlib import Path from typing import Iterable, List, Tuple @@ -128,7 +129,10 @@ def group_mean_connectivity( matrices = [np.asarray(cm.load(), dtype=np.float64) for cm in connectivity_matrices] matrices_vec = [sym_matrix_to_vec(mat, discard_diagonal=False) for mat in matrices] - mean_vec = np.nanmean(matrices_vec, axis=0) + with warnings.catch_warnings(): + # Edges that are NaN across all subjects yield all-NaN slices; the NaN is intended. + warnings.simplefilter("ignore", RuntimeWarning) + mean_vec = np.nanmean(matrices_vec, axis=0) mean_matrix = vec_to_sym_matrix(mean_vec, diagonal=None) return mean_matrix diff --git a/wonkyconn/tests/test_textual_app.py b/wonkyconn/tests/test_textual_app.py new file mode 100644 index 00000000..17ad2e24 --- /dev/null +++ b/wonkyconn/tests/test_textual_app.py @@ -0,0 +1,7 @@ +import textual.app + + +def test_import_textual_app() -> None: + from wonkyconn import textual_app + + assert issubclass(textual_app.WonkyConnApp, textual.app.App) diff --git a/wonkyconn/tests/test_workflow.py b/wonkyconn/tests/test_workflow.py new file mode 100644 index 00000000..ada66323 --- /dev/null +++ b/wonkyconn/tests/test_workflow.py @@ -0,0 +1,23 @@ +import argparse +from pathlib import Path + +import pandas as pd +import pytest + +from wonkyconn.workflow import load_data_frame + + +def test_load_data_frame_requires_site_column(tmp_path: Path) -> None: + """Enabling site correction requires a 'site' column in the phenotypes file.""" + phenotypes_path = tmp_path / "participants.tsv" + pd.DataFrame( + dict( + participant_id=["sub-1", "sub-2"], + age=[30.0, 40.0], + gender=["m", "f"], + ) + ).to_csv(phenotypes_path, sep="\t", index=False) + + args = argparse.Namespace(phenotypes=phenotypes_path, site_correction=True) + with pytest.raises(ValueError, match="site"): + load_data_frame(args) From 05ae4b043bb7972c12d54d2dccfe45f5d0aa3cfb Mon Sep 17 00:00:00 2001 From: Lea Waller Date: Sun, 2 Aug 2026 17:54:13 +0200 Subject: [PATCH 23/33] Update outdated docstrings --- wonkyconn/features/age_sex_prediction.py | 14 +++++++++----- .../features/calculate_gradients_correlation.py | 4 ++-- wonkyconn/features/quality_control_connectivity.py | 4 +++- 3 files changed, 14 insertions(+), 8 deletions(-) diff --git a/wonkyconn/features/age_sex_prediction.py b/wonkyconn/features/age_sex_prediction.py index 297df741..33b583fb 100644 --- a/wonkyconn/features/age_sex_prediction.py +++ b/wonkyconn/features/age_sex_prediction.py @@ -43,12 +43,15 @@ def fit(self: SiteRegressor, X: pd.DataFrame, y: pd.DataFrame | None = None) -> Fit the site regressor to the connectivity data. Args: - connectivity_data_site: A 2D array where the first n_connectivity_features columns - are connectivity features and the remaining columns are site dummy variables. + X: Connectivity data with one row per subject. Its index is used to + look up the corresponding site labels in ``self.sites``. y: Ignored. This parameter exists for compatibility with the scikit-learn API. Returns: self: Returns the instance itself. + + Raises: + ValueError: If the training data contains fewer than two sites. """ fit_sites = np.asarray(self.sites[X.index]) if np.unique(fit_sites).size < 2: @@ -63,12 +66,13 @@ def fit(self: SiteRegressor, X: pd.DataFrame, y: pd.DataFrame | None = None) -> def transform(self: SiteRegressor, X: pd.DataFrame) -> pd.DataFrame: # noqa: N803 """ Transform the connectivity data by regressing out site effects. + Args: - connectivity_data_site: A 2D array where the first n_connectivity_features columns - are connectivity features and the remaining columns are site dummy variables. + X: Connectivity data to correct, with one row per subject. Its index + is used to look up the corresponding site labels in ``self.sites``. Returns: - A 2D array of connectivity features with site effects regressed out. + pd.DataFrame: Connectivity data with the fitted site effects removed. """ return X - self.model.predict(self._get_dummies(X)) diff --git a/wonkyconn/features/calculate_gradients_correlation.py b/wonkyconn/features/calculate_gradients_correlation.py index dfa20c67..f9214050 100644 --- a/wonkyconn/features/calculate_gradients_correlation.py +++ b/wonkyconn/features/calculate_gradients_correlation.py @@ -48,8 +48,8 @@ def remove_nan_roi_atlas(atlas: nib.Nifti1Image, kept_idx: np.ndarray) -> Nifti1 kept_idx (np.ndarray): Indices of rows/columns kept after NaN removal. Returns: - tuple[nib.Nifti1Image, list[int]]: Atlas image with only kept ROIs and - the list of kept label values. + nib.Nifti1Image: Atlas image containing only the kept ROIs, with removed + regions set to zero. """ # Load atlas diff --git a/wonkyconn/features/quality_control_connectivity.py b/wonkyconn/features/quality_control_connectivity.py index 0b875f33..f79903dc 100644 --- a/wonkyconn/features/quality_control_connectivity.py +++ b/wonkyconn/features/quality_control_connectivity.py @@ -28,7 +28,7 @@ def calculate_qcfc( For each edge, we then computed the correlation between the weight of that edge and the mean relative RMS motion. QC-FC relationships were calculated as partial correlations that - accounted for participant age and sex + accounted for participant age and sex (and site when ``site_correction`` is enabled). Parameters: data_frame (pd.DataFrame): The data frame containing the covariates "age" and "gender". @@ -36,6 +36,8 @@ def calculate_qcfc( It needs to have one row for each connectivity matrix. connectivity_matrices (Iterable[ConnectivityMatrix]): The connectivity matrices to calculate QCFC for. metric_key (str, optional): The key of the metric to use for QCFC calculation. Defaults to "MeanFramewiseDisplacement". + site_correction (bool, optional): If True, include site as an additional covariate + (requires a "site" column in ``data_frame``). Defaults to False. Returns: pd.DataFrame: The QCFC values between connectivity matrices and the metric. From 5c0d2f5e1e66b9366b03be344470c1156b7aa097 Mon Sep 17 00:00:00 2001 From: Lea Waller Date: Sun, 2 Aug 2026 17:56:39 +0200 Subject: [PATCH 24/33] Remove getattr calls for arguments Argparse already makes sure that every argument exists as an attribute, set to either the user-specified value, the default value or None if no default was given. This means that the default values defined in the getattr call will be ignored in most cases, so keeping them may lead to misunderstandings down the line. --- wonkyconn/config.py | 22 +++++++++++----------- wonkyconn/run.py | 8 ++++---- wonkyconn/tests/test_cli.py | 4 ++-- wonkyconn/textual_app.py | 10 +++++----- wonkyconn/workflow.py | 12 ++++-------- 5 files changed, 26 insertions(+), 30 deletions(-) diff --git a/wonkyconn/config.py b/wonkyconn/config.py index 291375b9..d99cf638 100644 --- a/wonkyconn/config.py +++ b/wonkyconn/config.py @@ -14,7 +14,7 @@ def _coerce_path(value: str | Path | None) -> Path | None: @dataclass -class WonkyConnConfig: +class WonkyconnConfig: """Shared configuration for CLI and GUI.""" bids_dir: Path | None = None @@ -30,7 +30,7 @@ class WonkyConnConfig: site_correction: bool = False @classmethod - def from_cli_args(cls, args: argparse.Namespace | None) -> "WonkyConnConfig": + def from_cli_args(cls, args: argparse.Namespace | None) -> "WonkyconnConfig": """Create a config from argparse args (may be partial when GUI is requested).""" if args is None: return cls() @@ -40,20 +40,20 @@ def from_cli_args(cls, args: argparse.Namespace | None) -> "WonkyConnConfig": verbosity = verbosity[0] atlas_entries: list[tuple[str, Path]] = list() - for label, atlas_path in getattr(args, "atlas", []) or []: + for label, atlas_path in args.atlas or []: atlas_entries.append((label, Path(atlas_path).expanduser().resolve())) return cls( - bids_dir=_coerce_path(getattr(args, "bids_dir", None)), - output_dir=_coerce_path(getattr(args, "output_dir", None)), - analysis_level=getattr(args, "analysis_level", "group"), - phenotypes=_coerce_path(getattr(args, "phenotypes", None)), + bids_dir=_coerce_path(args.bids_dir), + output_dir=_coerce_path(args.output_dir), + analysis_level=args.analysis_level, + phenotypes=_coerce_path(args.phenotypes), atlas=atlas_entries, verbosity=int(verbosity) if verbosity is not None else 2, - debug=bool(getattr(args, "debug", False)), - light_mode=bool(getattr(args, "light_mode", False)), - suppress_warnings=bool(getattr(args, "suppress_warnings", False)), - site_correction=bool(getattr(args, "site_correction", False)), + debug=bool(args.debug), + light_mode=bool(args.light_mode), + suppress_warnings=bool(args.suppress_warnings), + site_correction=bool(args.site_correction), ) def to_namespace(self) -> argparse.Namespace: diff --git a/wonkyconn/run.py b/wonkyconn/run.py index 5b9864d8..a6e735aa 100644 --- a/wonkyconn/run.py +++ b/wonkyconn/run.py @@ -6,7 +6,7 @@ from typing import Sequence from . import __version__ -from .config import WonkyConnConfig +from .config import WonkyconnConfig from .logger import logger from .workflow import workflow @@ -108,11 +108,11 @@ def _loosen_parser_for_gui(parser: argparse.ArgumentParser) -> None: action.required = False -def _build_initial_config(args: argparse.Namespace | None) -> WonkyConnConfig: - return WonkyConnConfig.from_cli_args(args) +def _build_initial_config(args: argparse.Namespace | None) -> WonkyconnConfig: + return WonkyconnConfig.from_cli_args(args) -def _run_textual_ui(config: WonkyConnConfig) -> WonkyConnConfig | None: +def _run_textual_ui(config: WonkyconnConfig) -> WonkyconnConfig | None: if not sys.stdout.isatty(): print("The Textual UI requires an interactive terminal; run without --textual instead.", file=sys.stderr) sys.exit(2) diff --git a/wonkyconn/tests/test_cli.py b/wonkyconn/tests/test_cli.py index 318b99a2..7832b8dd 100644 --- a/wonkyconn/tests/test_cli.py +++ b/wonkyconn/tests/test_cli.py @@ -16,7 +16,7 @@ from tqdm.auto import tqdm from wonkyconn import __version__ -from wonkyconn.config import WonkyConnConfig +from wonkyconn.config import WonkyconnConfig from wonkyconn.file_index.bids import BIDSIndex from wonkyconn.run import global_parser, main from wonkyconn.workflow import workflow @@ -72,7 +72,7 @@ def test_cli_and_textual_namespace_consistency(tmp_path: Path) -> None: ) # Get namespace from WonkyConnConfig (used by Textual UI) - config = WonkyConnConfig( + config = WonkyconnConfig( bids_dir=bids_dir, output_dir=output_dir, analysis_level="group", diff --git a/wonkyconn/textual_app.py b/wonkyconn/textual_app.py index 6f742fcc..19bb9cee 100644 --- a/wonkyconn/textual_app.py +++ b/wonkyconn/textual_app.py @@ -21,10 +21,10 @@ from textual.widgets._select import NoSelection from textual.widgets._tree import Tree -from .config import WonkyConnConfig +from .config import WonkyconnConfig -class WonkyConnApp(App[WonkyConnConfig | None]): +class WonkyConnApp(App[WonkyconnConfig | None]): """Textual UI for configuring wonkyconn.""" CSS = """ @@ -107,7 +107,7 @@ class WonkyConnApp(App[WonkyConnConfig | None]): BINDINGS = [("escape", "cancel", "Cancel"), ("ctrl+s", "run", "Run")] - def __init__(self, initial_config: WonkyConnConfig): + def __init__(self, initial_config: WonkyconnConfig): super().__init__() self.initial_config = initial_config self.selected_path: Path | None = None @@ -301,7 +301,7 @@ def _set_status(self, message: str, error: bool = False) -> None: if error: status.add_class("status-error") - def _validate_and_build_config(self) -> WonkyConnConfig | None: + def _validate_and_build_config(self) -> WonkyconnConfig | None: errors: list[str] = list() bids_str = self._path_values.get("bids_dir", "").strip() @@ -348,7 +348,7 @@ def _validate_and_build_config(self) -> WonkyConnConfig | None: return None assert atlas_path is not None # validated above - config = WonkyConnConfig( + config = WonkyconnConfig( bids_dir=bids_dir, output_dir=output_dir, analysis_level="group", diff --git a/wonkyconn/workflow.py b/wonkyconn/workflow.py index 2ce00aa5..b5162f00 100644 --- a/wonkyconn/workflow.py +++ b/wonkyconn/workflow.py @@ -56,11 +56,6 @@ def workflow(args: argparse.Namespace) -> None: set_verbosity(args.verbosity) logger.debug(vars(args)) - # check if light mode is enabled - if so, it will not run the age and sex prediction and gradient similarity - disable_prediction_gradient = getattr(args, "light_mode", False) - - enable_site_correction = getattr(args, "site_correction", True) - # Check BIDS path bids_dir = args.bids_dir index = BIDSIndex() @@ -139,8 +134,9 @@ def workflow(args: argparse.Namespace) -> None: metric_key, seg_key, atlases, - disable_prediction_gradient, - enable_site_correction, + site_correction=args.site_correction, + # check if light mode is enabled - if so, it will not run the age and sex prediction and gradient similarity + disable_prediction_gradient=args.light_mode, ) record.update(dict(zip(group_by, key, strict=False))) if len(group_by) == 2: @@ -304,7 +300,7 @@ def load_data_frame(args: argparse.Namespace) -> pd.DataFrame: raise ValueError('Phenotypes file is missing the "gender" column') if "age" not in data_frame.columns: raise ValueError('Phenotypes file is missing the "age" column') - if getattr(args, "site_correction", True): + if args.site_correction: logger.info("Site correction is enabled - checking for 'site' column in phenotypes file.") if "site" not in data_frame.columns: raise ValueError('Phenotypes file is missing the "site" column required for site correction') From a2b562b233febf78463e9dc3f8a1588e57b578b6 Mon Sep 17 00:00:00 2001 From: Lea Waller Date: Sun, 2 Aug 2026 18:03:41 +0200 Subject: [PATCH 25/33] Expand test case to check for all column errors --- wonkyconn/tests/test_workflow.py | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/wonkyconn/tests/test_workflow.py b/wonkyconn/tests/test_workflow.py index ada66323..00309dab 100644 --- a/wonkyconn/tests/test_workflow.py +++ b/wonkyconn/tests/test_workflow.py @@ -7,17 +7,21 @@ from wonkyconn.workflow import load_data_frame -def test_load_data_frame_requires_site_column(tmp_path: Path) -> None: +@pytest.mark.parametrize("column", ["participant_id", "age", "gender", "site"]) +def test_load_data_frame_missing_column(tmp_path: Path, column: str) -> None: """Enabling site correction requires a 'site' column in the phenotypes file.""" + + phenotypes = dict( + participant_id=["sub-1", "sub-2"], + age=[30.0, 40.0], + gender=["m", "f"], + site=["site-a", "site-b"], + ) + del phenotypes[column] + phenotypes_path = tmp_path / "participants.tsv" - pd.DataFrame( - dict( - participant_id=["sub-1", "sub-2"], - age=[30.0, 40.0], - gender=["m", "f"], - ) - ).to_csv(phenotypes_path, sep="\t", index=False) + pd.DataFrame(phenotypes).to_csv(phenotypes_path, sep="\t", index=False) args = argparse.Namespace(phenotypes=phenotypes_path, site_correction=True) - with pytest.raises(ValueError, match="site"): + with pytest.raises(ValueError, match=column): load_data_frame(args) From c2ce58aabe850d4642d8280d00682d6daf459618 Mon Sep 17 00:00:00 2001 From: Lea Waller Date: Fri, 7 Aug 2026 15:35:09 +0200 Subject: [PATCH 26/33] Make stratified split failure more verbose --- wonkyconn/features/age_sex_prediction.py | 10 ++++++---- wonkyconn/features/tests/test_prediction.py | 21 +++++++++++++++++++++ 2 files changed, 27 insertions(+), 4 deletions(-) diff --git a/wonkyconn/features/age_sex_prediction.py b/wonkyconn/features/age_sex_prediction.py index 33b583fb..1bfcee37 100644 --- a/wonkyconn/features/age_sex_prediction.py +++ b/wonkyconn/features/age_sex_prediction.py @@ -108,21 +108,23 @@ def training_pipeline( y_train = pd.Series(LabelEncoder().fit_transform(target_labels)) # pyright: ignore[reportArgumentType, reportCallIssue] estimator = LogisticRegression(max_iter=5000, solver="lbfgs", random_state=random_state) - bins = y_train + bins = pd.Series(target_labels) scoring_metrics = {"accuracy": "accuracy", "roc_auc": "roc_auc"} else: y_train = pd.Series(target_labels) estimator = Ridge(alpha=1.0) - bins = pd.qcut(y_train, q=5, labels=False, duplicates="drop") + bins, edges = pd.qcut(y_train, q=5, labels=False, retbins=True, duplicates="drop") + labels = pd.Series([f"age-group-{int(round(edges[i]))}-{int(round(edges[i + 1]))}" for i in range(len(edges) - 1)]) + bins = bins.map(labels) scoring_metrics = {"mae": "neg_mean_absolute_error", "r2": "r2"} if sites is not None: data_frame = pd.DataFrame({"site": sites, "bins": bins}) - # Get unique row indices as combined bins - bins = data_frame.groupby(data_frame.columns.tolist(), sort=False).ngroup() + # Combine bins and sites + bins = data_frame.agg("_".join, axis=1) cv_strategy = StratifiedShuffleSplit(n_splits=n_splits, test_size=0.2, random_state=random_state) splits = list(cv_strategy.split(np.zeros(len(bins)), bins)) diff --git a/wonkyconn/features/tests/test_prediction.py b/wonkyconn/features/tests/test_prediction.py index 67608f7d..6bcf7eeb 100644 --- a/wonkyconn/features/tests/test_prediction.py +++ b/wonkyconn/features/tests/test_prediction.py @@ -176,3 +176,24 @@ def test_site_regressor_requires_two_sites() -> None: regressor = SiteRegressor(np.array(["site-a", "site-a"])) with pytest.raises(ValueError, match="at least two sites"): regressor.fit(pd.DataFrame(np.zeros((2, 3), dtype=np.float32))) + + +def test_training_pipeline_singleton_classes() -> None: + labels = np.array(["male"] * 10 + ["female"] * 10) + sites = np.array(["site-a"] * 5 + ["site-b"] * 5 + ["site-a"] * 9 + ["site-b"] * 1) + + rng = np.random.default_rng(random_state) + connectivity_data = rng.normal(size=(len(labels), feature_count)).astype(np.float32) + + with pytest.raises(ValueError, match=r"least populated classes in y have only 1 member") as exc_info: + training_pipeline( + connectivity_data, + labels, + task_type="classification", + n_splits=3, + n_pca=5, + n_jobs=1, + random_state=random_state, + sites=sites, + ) + assert "site-b_female" in str(exc_info.value) From e221347d7fbb7d28201c6f0c3c6dcc81308b8db6 Mon Sep 17 00:00:00 2001 From: Lea Waller Date: Fri, 7 Aug 2026 18:06:16 +0200 Subject: [PATCH 27/33] Change error to warning and exclude subjects that are the only subjects in their subgroups from the training pipeline --- wonkyconn/features/age_sex_prediction.py | 15 ++++++++++++++- wonkyconn/features/tests/test_prediction.py | 12 ++++++++---- 2 files changed, 22 insertions(+), 5 deletions(-) diff --git a/wonkyconn/features/age_sex_prediction.py b/wonkyconn/features/age_sex_prediction.py index 1bfcee37..e29a90b3 100644 --- a/wonkyconn/features/age_sex_prediction.py +++ b/wonkyconn/features/age_sex_prediction.py @@ -16,6 +16,8 @@ from sklearn.pipeline import Pipeline from sklearn.preprocessing import LabelEncoder, StandardScaler +from ..logger import logger + if TYPE_CHECKING: from ..base import ConnectivityMatrix @@ -126,8 +128,19 @@ def training_pipeline( # Combine bins and sites bins = data_frame.agg("_".join, axis=1) + counts = bins.value_counts() + singletons = counts.index[counts == 1] + mask = bins.isin(singletons).to_numpy() + if mask.any(): + count = mask.sum().item() + logger.warning(f"Excluding {count} {'subject'} that are the only subjects in their subgroups: {singletons.tolist()}") + keep = np.flatnonzero(np.logical_not(mask)) + cv_strategy = StratifiedShuffleSplit(n_splits=n_splits, test_size=0.2, random_state=random_state) - splits = list(cv_strategy.split(np.zeros(len(bins)), bins)) + splits = [ + (keep[train_indices], keep[test_indices]) + for train_indices, test_indices in cv_strategy.split(np.zeros(keep.size), bins.to_numpy()[keep]) + ] steps: list[tuple[str, BaseEstimator]] = [ # keep_empty_features=True avoids sklearn's "Skipping features without any diff --git a/wonkyconn/features/tests/test_prediction.py b/wonkyconn/features/tests/test_prediction.py index 6bcf7eeb..b0568a43 100644 --- a/wonkyconn/features/tests/test_prediction.py +++ b/wonkyconn/features/tests/test_prediction.py @@ -178,15 +178,15 @@ def test_site_regressor_requires_two_sites() -> None: regressor.fit(pd.DataFrame(np.zeros((2, 3), dtype=np.float32))) -def test_training_pipeline_singleton_classes() -> None: +def test_training_pipeline_exclude_singleton_classes(caplog: pytest.LogCaptureFixture) -> None: labels = np.array(["male"] * 10 + ["female"] * 10) sites = np.array(["site-a"] * 5 + ["site-b"] * 5 + ["site-a"] * 9 + ["site-b"] * 1) rng = np.random.default_rng(random_state) connectivity_data = rng.normal(size=(len(labels), feature_count)).astype(np.float32) - with pytest.raises(ValueError, match=r"least populated classes in y have only 1 member") as exc_info: - training_pipeline( + with caplog.at_level("WARNING"): + summary = training_pipeline( connectivity_data, labels, task_type="classification", @@ -196,4 +196,8 @@ def test_training_pipeline_singleton_classes() -> None: random_state=random_state, sites=sites, ) - assert "site-b_female" in str(exc_info.value) + + assert isinstance(summary, pd.DataFrame) + assert set(summary.index) == {"accuracy", "roc_auc"} + assert np.isfinite(summary.to_numpy()).all() + assert "site-b_female" in caplog.text From 2c17820c2ad11e4536c4bdca3f017d4a52a35c7b Mon Sep 17 00:00:00 2001 From: Lea Waller Date: Fri, 7 Aug 2026 20:44:16 +0200 Subject: [PATCH 28/33] Fix stratify for numeric site column Smoke test fails for ds000030 --- wonkyconn/features/age_sex_prediction.py | 2 +- wonkyconn/features/tests/test_prediction.py | 23 +++++++++++++++++++++ 2 files changed, 24 insertions(+), 1 deletion(-) diff --git a/wonkyconn/features/age_sex_prediction.py b/wonkyconn/features/age_sex_prediction.py index e29a90b3..7447f0a9 100644 --- a/wonkyconn/features/age_sex_prediction.py +++ b/wonkyconn/features/age_sex_prediction.py @@ -124,7 +124,7 @@ def training_pipeline( scoring_metrics = {"mae": "neg_mean_absolute_error", "r2": "r2"} if sites is not None: - data_frame = pd.DataFrame({"site": sites, "bins": bins}) + data_frame = pd.DataFrame({"site": sites, "bins": bins}).astype(str) # Combine bins and sites bins = data_frame.agg("_".join, axis=1) diff --git a/wonkyconn/features/tests/test_prediction.py b/wonkyconn/features/tests/test_prediction.py index b0568a43..5a0e464f 100644 --- a/wonkyconn/features/tests/test_prediction.py +++ b/wonkyconn/features/tests/test_prediction.py @@ -201,3 +201,26 @@ def test_training_pipeline_exclude_singleton_classes(caplog: pytest.LogCaptureFi assert set(summary.index) == {"accuracy", "roc_auc"} assert np.isfinite(summary.to_numpy()).all() assert "site-b_female" in caplog.text + + +def test_training_pipeline_numeric_sites() -> None: + labels = np.array(["male", "female"] * 10) + sites = np.array([1.0, 2.0] * 10) + + rng = np.random.default_rng(random_state) + connectivity_data = rng.normal(size=(len(labels), feature_count)).astype(np.float32) + + summary = training_pipeline( + connectivity_data, + labels, + task_type="classification", + n_splits=3, + n_pca=5, + n_jobs=1, + random_state=random_state, + sites=sites, + ) + + assert isinstance(summary, pd.DataFrame) + assert set(summary.index) == {"accuracy", "roc_auc"} + assert np.isfinite(summary.to_numpy()).all() From d12c7924faa101d72fadd0d1195e86896431b2e4 Mon Sep 17 00:00:00 2001 From: Lea Waller Date: Mon, 10 Aug 2026 16:48:50 +0200 Subject: [PATCH 29/33] Allow specifying specific metrics to run in cli --- .github/workflows/docker.yml | 2 +- wonkyconn/config.py | 42 +++---- wonkyconn/run.py | 28 +++-- wonkyconn/tests/test_cli.py | 71 +----------- wonkyconn/tests/test_workflow.py | 6 +- wonkyconn/textual_app.py | 4 +- wonkyconn/visualization/plot.py | 21 +--- wonkyconn/workflow.py | 190 ++++++++++++++----------------- 8 files changed, 135 insertions(+), 229 deletions(-) diff --git a/.github/workflows/docker.yml b/.github/workflows/docker.yml index 01b0ce5f..8b9c0985 100644 --- a/.github/workflows/docker.yml +++ b/.github/workflows/docker.yml @@ -62,7 +62,7 @@ jobs: cache-to: type=gha,mode=max tags: wonky:conn load: true - # always use light mode for testing container since the full version was tested during build + # Always use light mode for testing container since the full version was tested during build - run: | datalad get data/halfpipe data/atlases docker run --rm \ diff --git a/wonkyconn/config.py b/wonkyconn/config.py index d99cf638..ac99977d 100644 --- a/wonkyconn/config.py +++ b/wonkyconn/config.py @@ -3,7 +3,7 @@ import argparse from dataclasses import dataclass, field from pathlib import Path -from typing import Iterable, Sequence +from typing import Literal, Sequence, TypeAlias def _coerce_path(value: str | Path | None) -> Path | None: @@ -13,6 +13,12 @@ def _coerce_path(value: str | Path | None) -> Path | None: return Path(value).expanduser().resolve() +Metric: TypeAlias = Literal["motion", "analytic-insights", "gradients", "prediction"] + +light_mode_metrics: set[Metric] = {"motion", "analytic-insights"} +all_metrics: set[Metric] = {"motion", "analytic-insights", "gradients", "prediction"} + + @dataclass class WonkyconnConfig: """Shared configuration for CLI and GUI.""" @@ -24,11 +30,15 @@ class WonkyconnConfig: atlas: list[tuple[str, Path]] = field(default_factory=list) verbosity: int = 2 debug: bool = False - light_mode: bool = False + metrics: set[Metric] | None = None theme: str | None = None # GUI-only suppress_warnings: bool = False site_correction: bool = False + @property + def light_mode(self) -> bool: + return self.metrics == light_mode_metrics + @classmethod def from_cli_args(cls, args: argparse.Namespace | None) -> "WonkyconnConfig": """Create a config from argparse args (may be partial when GUI is requested).""" @@ -43,6 +53,8 @@ def from_cli_args(cls, args: argparse.Namespace | None) -> "WonkyconnConfig": for label, atlas_path in args.atlas or []: atlas_entries.append((label, Path(atlas_path).expanduser().resolve())) + metrics: set[Metric] = light_mode_metrics if args.light_mode else set(args.metrics) + return cls( bids_dir=_coerce_path(args.bids_dir), output_dir=_coerce_path(args.output_dir), @@ -51,31 +63,7 @@ def from_cli_args(cls, args: argparse.Namespace | None) -> "WonkyconnConfig": atlas=atlas_entries, verbosity=int(verbosity) if verbosity is not None else 2, debug=bool(args.debug), - light_mode=bool(args.light_mode), + metrics=metrics, suppress_warnings=bool(args.suppress_warnings), site_correction=bool(args.site_correction), ) - - def to_namespace(self) -> argparse.Namespace: - """Convert to argparse.Namespace expected by workflow.""" - if self.bids_dir is None: - raise ValueError("bids_dir is required") - if self.output_dir is None: - raise ValueError("output_dir is required") - if self.phenotypes is None: - raise ValueError("phenotypes is required") - if not self.atlas: - raise ValueError("At least one atlas entry is required") - - atlas_as_str: Iterable[tuple[str, str]] = ((label, str(path)) for label, path in self.atlas) - return argparse.Namespace( - bids_dir=self.bids_dir, - output_dir=self.output_dir, - analysis_level=self.analysis_level, - phenotypes=str(self.phenotypes), - atlas=list(atlas_as_str), - verbosity=self.verbosity, - debug=self.debug, - light_mode=self.light_mode, - site_correction=self.site_correction, - ) diff --git a/wonkyconn/run.py b/wonkyconn/run.py index a6e735aa..26d00120 100644 --- a/wonkyconn/run.py +++ b/wonkyconn/run.py @@ -57,15 +57,25 @@ def global_parser(exit_on_error: bool = True) -> argparse.ArgumentParser: help="Specify the atlas label and the path to the atlas file (for example --atlas Schaefer2018 /path/to/atlas.nii.gz)", ) - parser.add_argument("-v", "--version", action="version", version=__version__) - parser.add_argument("--debug", action="store_true", default=False) - parser.add_argument( + metrics_group = parser.add_mutually_exclusive_group(required=False) + metrics_group.add_argument( "--light-mode", required=False, action="store_true", default=False, help="Disable sex and age prediction to reduce runtime.", ) + metrics = ["motion", "analytic-insights", "gradients", "prediction"] + metrics_group.add_argument( + "--metrics", + required=False, + nargs="+", + choices=metrics, + default=metrics, + ) + + parser.add_argument("-v", "--version", action="version", version=__version__) + parser.add_argument("--debug", action="store_true", default=False) parser.add_argument( "--site-correction", required=False, @@ -162,15 +172,13 @@ def main(argv: None | Sequence[str] = None) -> None: if use_gui: initial_config = _build_initial_config(parsed_args) - result_config = _run_textual_ui(initial_config) - if result_config is None: + config = _run_textual_ui(initial_config) + if config is None: return - args_for_workflow = result_config.to_namespace() - debug_enabled = result_config.debug - suppress_warnings = result_config.suppress_warnings + debug_enabled = config.debug + suppress_warnings = config.suppress_warnings else: config = _build_initial_config(parsed_args) - args_for_workflow = config.to_namespace() debug_enabled = config.debug suppress_warnings = config.suppress_warnings @@ -180,7 +188,7 @@ def main(argv: None | Sequence[str] = None) -> None: warnings.filterwarnings("ignore", category=RuntimeWarning) try: - workflow(args_for_workflow) + workflow(config) except Exception as e: logger.exception("Exception: %s", e, exc_info=True) if debug_enabled: diff --git a/wonkyconn/tests/test_cli.py b/wonkyconn/tests/test_cli.py index 7832b8dd..e6edf7a2 100644 --- a/wonkyconn/tests/test_cli.py +++ b/wonkyconn/tests/test_cli.py @@ -40,68 +40,6 @@ def test_help(capsys: pytest.CaptureFixture[str]) -> None: assert "Evaluating the residual motion in fMRI connectome and visualize reports" in captured.out -def test_cli_and_textual_namespace_consistency(tmp_path: Path) -> None: - """Ensure CLI (run.py) and Textual UI produce namespaces with the same attributes. - - Both interfaces should produce namespaces that workflow() can consume, - meaning they must have identical attribute names. - """ - # Create minimal valid paths for testing - bids_dir = tmp_path / "bids" - bids_dir.mkdir() - output_dir = tmp_path / "output" - output_dir.mkdir() - phenotypes = tmp_path / "participants.tsv" - phenotypes.touch() - atlas_path = tmp_path / "atlas.nii.gz" - atlas_path.touch() - - # Get namespace from CLI parser - parser = global_parser() - cli_args = parser.parse_args( - [ - str(bids_dir), - str(output_dir), - "group", - "--phenotypes", - str(phenotypes), - "--atlas", - "TestAtlas", - str(atlas_path), - ] - ) - - # Get namespace from WonkyConnConfig (used by Textual UI) - config = WonkyconnConfig( - bids_dir=bids_dir, - output_dir=output_dir, - analysis_level="group", - phenotypes=phenotypes, - atlas=[("TestAtlas", atlas_path)], - verbosity=2, - debug=False, - ) - config_namespace = config.to_namespace() - - # Get the attribute names from both namespaces - cli_attrs = set(vars(cli_args).keys()) - config_attrs = set(vars(config_namespace).keys()) - - # Attributes that are interface-specific and not passed to workflow() - # These are handled separately before calling workflow() - interface_specific_attrs = {"textual", "version", "suppress_warnings"} - - # The workflow-relevant CLI attributes (excluding interface-specific ones) - workflow_cli_attrs = cli_attrs - interface_specific_attrs - - # Both should have the same attributes for workflow consumption - assert workflow_cli_attrs == config_attrs, ( - f"Namespace mismatch!\n" - f"CLI-only attrs (not in config): {workflow_cli_attrs - config_attrs}\n" - f"Config-only attrs (not in CLI): {config_attrs - workflow_cli_attrs}" - ) - - def _copy_file(path: Path, new_path: Path, sub: str) -> None: new_path = Path(re.sub(r"sub-\d+", sub, str(new_path))) new_path.parent.mkdir(parents=True, exist_ok=True) @@ -178,7 +116,8 @@ def test_giga_connectome(data_path: Path, tmp_path: Path) -> None: ] args = parser.parse_args(argv) - workflow(args) + config = WonkyconnConfig.from_cli_args(args) + workflow(config) assert (output_dir / "metrics.tsv").is_file() assert (output_dir / "metrics.png").is_file() @@ -221,7 +160,8 @@ def test_halfpipe(data_path: Path, tmp_path: Path, site_correction: bool) -> Non argv.insert(0, "--site-correction") args = parser.parse_args(argv) - workflow(args) + config = WonkyconnConfig.from_cli_args(args) + workflow(config) # Add persistent storage to extract figure as artifact persistent_dir = Path("figures_artifacts") @@ -273,7 +213,8 @@ def test_halfpipe_with_full_metrics(data_path: Path, tmp_path: Path, site_correc argv.insert(0, "--site-correction") args = parser.parse_args(argv) - workflow(args) + config = WonkyconnConfig.from_cli_args(args) + workflow(config) # Add persistent storage to extract figure as artifact persistent_dir = Path("figures_artifacts") diff --git a/wonkyconn/tests/test_workflow.py b/wonkyconn/tests/test_workflow.py index 00309dab..b60787c1 100644 --- a/wonkyconn/tests/test_workflow.py +++ b/wonkyconn/tests/test_workflow.py @@ -1,9 +1,9 @@ -import argparse from pathlib import Path import pandas as pd import pytest +from wonkyconn.config import WonkyconnConfig from wonkyconn.workflow import load_data_frame @@ -22,6 +22,6 @@ def test_load_data_frame_missing_column(tmp_path: Path, column: str) -> None: phenotypes_path = tmp_path / "participants.tsv" pd.DataFrame(phenotypes).to_csv(phenotypes_path, sep="\t", index=False) - args = argparse.Namespace(phenotypes=phenotypes_path, site_correction=True) + config = WonkyconnConfig(phenotypes=phenotypes_path, site_correction=True) with pytest.raises(ValueError, match=column): - load_data_frame(args) + load_data_frame(config) diff --git a/wonkyconn/textual_app.py b/wonkyconn/textual_app.py index 19bb9cee..cfe5f966 100644 --- a/wonkyconn/textual_app.py +++ b/wonkyconn/textual_app.py @@ -21,7 +21,7 @@ from textual.widgets._select import NoSelection from textual.widgets._tree import Tree -from .config import WonkyconnConfig +from .config import WonkyconnConfig, all_metrics, light_mode_metrics class WonkyConnApp(App[WonkyconnConfig | None]): @@ -356,7 +356,7 @@ def _validate_and_build_config(self) -> WonkyconnConfig | None: atlas=[(atlas_label, atlas_path)], verbosity=verbosity, debug=debug, - light_mode=light_mode, + metrics=light_mode_metrics if light_mode else all_metrics, theme="dark" if self.dark else "light", suppress_warnings=suppress_warnings, ) diff --git a/wonkyconn/visualization/plot.py b/wonkyconn/visualization/plot.py index 6be66469..5d565d3a 100644 --- a/wonkyconn/visualization/plot.py +++ b/wonkyconn/visualization/plot.py @@ -39,23 +39,6 @@ def plot(records: list[dict[str, Any]], group_by: list[str], output_dir: Path) - group_by (list[str]): The list of columns that the results are grouped by. output_dir (Path): The directory to save the plot image into as "metrics.png". """ - # separate dmn similarity from the rest of the metrics - dmn_sim_array = [] - for record in records: - df_dmn_similarity = record.pop("dmn_similarity") - for g in group_by: - df_dmn_similarity[g] = record[g] - dmn_sim_array.append(df_dmn_similarity[["corr_with_dmn"] + group_by]) - - df_dmn_sim_array = pd.concat(dmn_sim_array, ignore_index=True) - if len(group_by) == 2: - df_dmn_sim_array["group_labels"] = df_dmn_sim_array[group_by].apply(lambda x: "-".join(x.astype(str)), axis=1) - else: - df_dmn_sim_array["group_labels"] = df_dmn_sim_array[group_by[0]] - - # summarize the info - for record, dmn_sim in zip(records, dmn_sim_array, strict=True): - record["dmn_similarity"] = dmn_sim["corr_with_dmn"].mean() result_frame = pd.DataFrame.from_records(records, index=group_by) data_frame = result_frame.reset_index() if len(group_by) == 2: # halfpipe @@ -113,7 +96,7 @@ def plot(records: list[dict[str, Any]], group_by: list[str], output_dir: Path) - gcor_axes.set_title("Global correlation (GCOR)") gcor_axes.set_xlabel("Mean correlation") - sns.barplot(data=df_dmn_sim_array, y="group_labels", x="corr_with_dmn", color=palette[4], ax=dmn_mean_axes, errorbar="sd") + sns.barplot(data=data_frame, y="group_labels", x="dmn_similarity_mean", color=palette[4], ax=dmn_mean_axes, errorbar="sd") dmn_mean_axes.set_title("Similarity with DMN") dmn_mean_axes.set_xlabel("Mean correlation") @@ -181,7 +164,7 @@ def plot_degrees_of_freedom_loss( result_frame: pd.DataFrame, degrees_of_freedom_loss_axes: Axes, legend_axes: Axes, - colors: list[str], + colors: list[Any], ) -> None: """Plot stacked bars showing degrees-of-freedom loss by source.""" sns.barplot( diff --git a/wonkyconn/workflow.py b/wonkyconn/workflow.py index b5162f00..8fcf6e0a 100644 --- a/wonkyconn/workflow.py +++ b/wonkyconn/workflow.py @@ -2,7 +2,6 @@ Process fMRIPrep outputs to timeseries based on denoising strategy. """ -import argparse import sys from collections import defaultdict, namedtuple from pathlib import Path @@ -13,6 +12,8 @@ from numpy import typing as npt from tqdm.auto import tqdm +from wonkyconn.config import Metric, WonkyconnConfig + from .atlas import Atlas from .base import ConnectivityMatrix from .features.age_sex_prediction import age_sex_scores @@ -50,14 +51,16 @@ def is_halfpipe(index: BIDSIndex) -> bool: return False -def workflow(args: argparse.Namespace) -> None: +def workflow(config: WonkyconnConfig) -> None: """Run the group-level connectivity quality-control pipeline.""" if "pytest" not in sys.modules: - set_verbosity(args.verbosity) - logger.debug(vars(args)) + set_verbosity(config.verbosity) + logger.debug(vars(config)) # Check BIDS path - bids_dir = args.bids_dir + bids_dir = config.bids_dir + if bids_dir is None: + raise ValueError("BIDS directory is not specified in the configuration") index = BIDSIndex() index.put(bids_dir) @@ -75,14 +78,16 @@ def workflow(args: argparse.Namespace) -> None: has_header = False # Check output path - output_dir = args.output_dir + output_dir = config.output_dir + if output_dir is None: + raise ValueError("Output directory is not specified in the configuration") output_dir.mkdir(parents=True, exist_ok=True) # Load data frame (participants: age, gender, etc.) - data_frame = load_data_frame(args) + data_frame = load_data_frame(config) # Load atlases - atlases: dict[str, Atlas] = {name: Atlas.create(name, Path(atlas_path_str)) for name, atlas_path_str in args.atlas} + atlases: dict[str, Atlas] = {name: Atlas.create(name, Path(atlas_path_str)) for name, atlas_path_str in config.atlas} logger.debug(f"Atlas dictionary contains: {list(atlases.keys())}") Group = namedtuple("Group", group_by) # type: ignore[misc] @@ -125,6 +130,10 @@ def workflow(args: argparse.Namespace) -> None: records: list[dict[str, Any]] = list() for key, connectivity_matrices in tqdm(grouped_connectivity_matrix.items(), unit="groups"): + if len(group_by) == 2: + dmn_similarity_path = output_dir / f"dmn_similarity_{'-'.join(group_by)}.tsv" + else: + dmn_similarity_path = output_dir / f"dmn_similarity_{group_by[0]}.tsv" record = make_record( index, data_frame, @@ -134,28 +143,16 @@ def workflow(args: argparse.Namespace) -> None: metric_key, seg_key, atlases, - site_correction=args.site_correction, - # check if light mode is enabled - if so, it will not run the age and sex prediction and gradient similarity - disable_prediction_gradient=args.light_mode, + dmn_similarity_path, + site_correction=config.site_correction, + metrics=config.metrics, ) record.update(dict(zip(group_by, key, strict=False))) - if len(group_by) == 2: - record["dmn_similarity"].to_csv(output_dir / f"dmn_similarity_{'-'.join(group_by)}.tsv", sep="\t") - else: - record["dmn_similarity"].to_csv(output_dir / f"dmn_similarity_{group_by[0]}.tsv", sep="\t") - - dmn_similarity_std = record["dmn_similarity"].loc[:, "corr_with_dmn"].std() - dmn_similarity_avg = record["dmn_similarity"].loc[:, "corr_with_dmn"].mean() - record["dmn_similarity_std"] = dmn_similarity_std - record["dmn_similarity_mean"] = dmn_similarity_avg records.append(record) plot(records, group_by, output_dir) - for record in records: - record.pop("dmn_similarity") - result_frame = pd.DataFrame.from_records(records, index=group_by) result_frame.to_csv(output_dir / "metrics.tsv", sep="\t") @@ -169,7 +166,8 @@ def make_record( metric_key: str, seg_key: str, atlases: dict[str, Atlas], - disable_prediction_gradient: bool, + dmn_similarity_path: Path, + metrics: set[Metric] | None = None, site_correction: bool = False, ) -> dict[str, Any]: """Compute all QC metrics for a single group of connectivity matrices.""" @@ -192,106 +190,94 @@ def make_record( else: logger.info(f"Skipping subject {sub}: not found in phenotype file.") - # Renaming for consistency + # Renaming for consistency connectivity_matrices[:] = filtered - - # Slice phenotypes (age, gender, etc.) for just this group seg_data_frame = data_frame.loc[seg_subjects] - qcfc = calculate_qcfc(seg_data_frame, connectivity_matrices, metric_key, site_correction) (seg,) = index.get_tag_values(seg_key, {c.path for c in connectivity_matrices}) distance_matrix = distance_matrices[seg] - - gcor = calculate_gcor(connectivity_matrices) - - dmn_similarity_summary, t_stats_dmn_vis_fpn = network_similarity(connectivity_matrices, region_memberships[seg]) atlas = atlases[seg].image - record = dict( - median_absolute_qcfc=calculate_median_absolute(qcfc.correlation), - percentage_significant_qcfc=calculate_qcfc_percentage(qcfc), - distance_dependence=calculate_distance_dependence(qcfc, distance_matrix), - gcor=gcor, - dmn_similarity=dmn_similarity_summary, - dmn_vis_distance_vs_dmn_fpn=t_stats_dmn_vis_fpn, + record: dict[str, Any] = dict( + # Motion + median_absolute_qcfc=np.nan, + percentage_significant_qcfc=np.nan, + distance_dependence=np.nan, + gcor=np.nan, + # Analytic insights + dmn_similarity_std=np.nan, + dmn_similarity_mean=np.nan, + dmn_vis_distance_vs_dmn_fpn=np.nan, + # Gradients + gradients_similarity=np.nan, + # Prediction + sex_auc=np.nan, + sex_auc_ci_lower=np.nan, + sex_auc_ci_upper=np.nan, + sex_accuracy=np.nan, + age_mae=np.nan, + age_mae_ci_lower=np.nan, + age_mae_ci_upper=np.nan, + age_r2=np.nan, + # Degrees of freedom loss **calculate_degrees_of_freedom_loss(connectivity_matrices)._asdict(), ) - if disable_prediction_gradient: - logger.info("Light mode enabled - skipping age and sex prediction, gradient similarity.") + if metrics is None: + raise ValueError("Metrics set is not specified") + + if "motion" in metrics: + qcfc = calculate_qcfc(seg_data_frame, connectivity_matrices, metric_key, site_correction) record.update( - dict( - sex_auc=np.nan, - sex_auc_ci_lower=np.nan, - sex_auc_ci_upper=np.nan, - sex_accuracy=np.nan, - age_mae=np.nan, - age_mae_ci_lower=np.nan, - age_mae_ci_upper=np.nan, - age_r2=np.nan, - gradients_similarity=np.nan, - ) - ) # place holders - return record - # Gradient similarity - gradients, gradients_group = extract_gradients(connectivity_matrices, atlas) - record["gradients_similarity"] = calculate_gradients_similarity(gradients, gradients_group) - - # age / sex predictability metrics - try: - ages = seg_data_frame["age"].to_numpy() - genders = seg_data_frame["gender"].to_numpy() - sites = seg_data_frame["site"].to_numpy() if site_correction else None - - scores = age_sex_scores( - connectivity_matrices, - ages=ages, - genders=genders, - sites=sites, - n_splits=_DEFAULT_N_SPLITS, - random_state=42, - n_pca=_DEFAULT_N_PCA, - n_jobs=_DEFAULT_N_JOBS, + median_absolute_qcfc=calculate_median_absolute(qcfc.correlation), + percentage_significant_qcfc=calculate_qcfc_percentage(qcfc), + distance_dependence=calculate_distance_dependence(qcfc, distance_matrix), ) - # scores is: - # { - # "sex_auc": float, - # "sex_auc_ci_lower": float, - # "sex_auc_ci_upper": float, - # "sex_accuracy": float, - # "age_mae": float, - # "age_mae_ci_lower": float, - # "age_mae_ci_upper": float, - # "age_r2": float, - # } - record.update(scores) - - except (ValueError, np.linalg.LinAlgError) as exc: - logger.warning(f"[age_sex_prediction] Skipping age/sex prediction for this group due to error: {exc!r}") - # If it fails, we still want consistent columns in the output. + if "analytic-insights" in metrics: + gcor = calculate_gcor(connectivity_matrices) + dmn_similarity_summary, t_stats_dmn_vis_fpn = network_similarity(connectivity_matrices, region_memberships[seg]) + dmn_similarity_summary.to_csv(dmn_similarity_path, sep="\t") record.update( - dict( - sex_auc=np.nan, - sex_auc_ci_lower=np.nan, - sex_auc_ci_upper=np.nan, - sex_accuracy=np.nan, - age_mae=np.nan, - age_mae_ci_lower=np.nan, - age_mae_ci_upper=np.nan, - age_r2=np.nan, - ) + gcor=gcor, + dmn_similarity_std=dmn_similarity_summary.loc[:, "corr_with_dmn"].std(), + dmn_similarity_mean=dmn_similarity_summary.loc[:, "corr_with_dmn"].mean(), + dmn_vis_distance_vs_dmn_fpn=t_stats_dmn_vis_fpn, ) + if "gradients" in metrics: + gradients, gradients_group = extract_gradients(connectivity_matrices, atlas) + record["gradients_similarity"] = calculate_gradients_similarity(gradients, gradients_group) + + if "prediction" in metrics: + try: + record.update( + age_sex_scores( + connectivity_matrices, + ages=seg_data_frame["age"].to_numpy(), + genders=seg_data_frame["gender"].to_numpy(), + sites=seg_data_frame["site"].to_numpy() if site_correction else None, + n_splits=_DEFAULT_N_SPLITS, + random_state=42, + n_pca=_DEFAULT_N_PCA, + n_jobs=_DEFAULT_N_JOBS, + ) + ) + except (ValueError, np.linalg.LinAlgError) as exc: + logger.warning(f"[age_sex_prediction] Skipping age/sex prediction for this group due to error: {exc!r}") + return record -def load_data_frame(args: argparse.Namespace) -> pd.DataFrame: +def load_data_frame(config: WonkyconnConfig) -> pd.DataFrame: """Load a phenotype TSV with ``participant_id``, ``gender``, and ``age`` columns. If site correction is enabled, the ``site`` column is also required. """ + path = config.phenotypes + if path is None: + raise ValueError("Phenotypes file path is not specified in the configuration") data_frame = pd.read_csv( - args.phenotypes, + path, sep="\t", index_col="participant_id", dtype={"participant_id": str}, @@ -300,7 +286,7 @@ def load_data_frame(args: argparse.Namespace) -> pd.DataFrame: raise ValueError('Phenotypes file is missing the "gender" column') if "age" not in data_frame.columns: raise ValueError('Phenotypes file is missing the "age" column') - if args.site_correction: + if config.site_correction: logger.info("Site correction is enabled - checking for 'site' column in phenotypes file.") if "site" not in data_frame.columns: raise ValueError('Phenotypes file is missing the "site" column required for site correction') From 189a424803390c148f35316e23b70dd9c3c0b3e7 Mon Sep 17 00:00:00 2001 From: Lea Waller Date: Mon, 10 Aug 2026 19:19:21 +0200 Subject: [PATCH 30/33] Improve index performance --- wonkyconn/file_index/base.py | 67 ++++++++++++++++++++++-------------- wonkyconn/file_index/bids.py | 1 + 2 files changed, 42 insertions(+), 26 deletions(-) diff --git a/wonkyconn/file_index/base.py b/wonkyconn/file_index/base.py index adf98c14..9c52c62a 100644 --- a/wonkyconn/file_index/base.py +++ b/wonkyconn/file_index/base.py @@ -5,33 +5,49 @@ from __future__ import annotations from collections import defaultdict -from hashlib import sha1 +from dataclasses import dataclass, field +from functools import cached_property from pathlib import Path -from typing import Mapping +from typing import Any, Mapping -def create_defaultdict_of_set() -> defaultdict[str, set[Path]]: - """Factory for a ``defaultdict`` whose values are sets of ``Path``.""" +def _defaultdict_of_sets() -> dict[str, set[Any]]: return defaultdict(set) -class FileIndex: - def __init__(self) -> None: - self.paths_by_tags: dict[str, dict[str, set[Path]]] = defaultdict(create_defaultdict_of_set) - self.tags_by_paths: dict[Path, dict[str, str]] = defaultdict(dict) +def _defaultdict_of_defaultdict_of_sets() -> dict[str, dict[str, set[Path]]]: + return defaultdict(_defaultdict_of_sets) + + +def _defaultdict_of_dict() -> dict[Path, dict[str, str]]: + return defaultdict(dict) - @property - def hexdigest(self) -> str: - """ - A forty character hash code of the paths in the index, obtained using the `sha1` algorithm. - """ - hash_algorithm = sha1() - for path in sorted(self.tags_by_paths.keys(), key=str): - path_bytes = str(path).encode() - hash_algorithm.update(path_bytes) +@dataclass +class FileIndex: + paths_by_tags: dict[str, dict[str, set[Path]]] = field(default_factory=_defaultdict_of_defaultdict_of_sets) + tags_by_paths: dict[Path, dict[str, str]] = field(default_factory=_defaultdict_of_dict) - return hash_algorithm.hexdigest() + @cached_property + def paths(self) -> set[Path]: + return set(self.tags_by_paths.keys()) + + def _get_phenotypes_without_key(self, key: str) -> set[Path]: + if key not in self.paths_by_tags: + return self.paths.copy() + else: + cache = self.__dict__.setdefault("_phenotypes_without_key", dict()) + phenotypes = cache.get(key) + if phenotypes is None: + phenotypes = self.paths.difference(*self.paths_by_tags[key].values()) + cache[key] = phenotypes + return phenotypes + + def _invalidate_caches(self) -> None: + if "phenotypes" in self.__dict__: + del self.__dict__["phenotypes"] + if "_phenotypes_without_key" in self.__dict__: + del self.__dict__["_phenotypes_without_key"] def get(self, **tags: str | None) -> set[Path]: """ @@ -47,21 +63,20 @@ def get(self, **tags: str | None) -> set[Path]: """ matches: set[Path] | None = None - for key, value in tags.items(): + for key, query in tags.items(): if key not in self.paths_by_tags: return set() values = self.paths_by_tags[key] - if value is not None: - if value not in values: - return set() - paths: set[Path] = values[value] + if query is None: + paths = self._get_phenotypes_without_key(key) + elif query in values: + paths = values[query] else: - paths_in_index = set(self.tags_by_paths.keys()) - paths = paths_in_index.difference(*values.values()) + return set() if matches is not None: - matches &= paths + matches.intersection_update(paths) else: matches = paths.copy() diff --git a/wonkyconn/file_index/bids.py b/wonkyconn/file_index/bids.py index 96e01a99..802a87a9 100644 --- a/wonkyconn/file_index/bids.py +++ b/wonkyconn/file_index/bids.py @@ -102,6 +102,7 @@ def parse(path: Path) -> dict[str, str] | None: class BIDSIndex(FileIndex): def put(self, root: Path) -> None: + self._invalidate_caches() for path in root.glob("**/*"): tags = parse(path) From 08019197b05cd20507c2277f5230cd0bf1b94717 Mon Sep 17 00:00:00 2001 From: Lea Waller Date: Mon, 10 Aug 2026 19:39:04 +0200 Subject: [PATCH 31/33] Add --log-level cli variable as an alternative to --verbosity --- wonkyconn/config.py | 10 +++----- wonkyconn/logger.py | 14 ----------- wonkyconn/run.py | 51 ++++++++++++++++++++++++++-------------- wonkyconn/textual_app.py | 18 +++++++------- wonkyconn/workflow.py | 6 ++--- 5 files changed, 48 insertions(+), 51 deletions(-) diff --git a/wonkyconn/config.py b/wonkyconn/config.py index ac99977d..54ee2d19 100644 --- a/wonkyconn/config.py +++ b/wonkyconn/config.py @@ -3,7 +3,7 @@ import argparse from dataclasses import dataclass, field from pathlib import Path -from typing import Literal, Sequence, TypeAlias +from typing import Literal, TypeAlias def _coerce_path(value: str | Path | None) -> Path | None: @@ -28,7 +28,7 @@ class WonkyconnConfig: analysis_level: str = "group" phenotypes: Path | None = None atlas: list[tuple[str, Path]] = field(default_factory=list) - verbosity: int = 2 + log_level: str | None = None debug: bool = False metrics: set[Metric] | None = None theme: str | None = None # GUI-only @@ -45,10 +45,6 @@ def from_cli_args(cls, args: argparse.Namespace | None) -> "WonkyconnConfig": if args is None: return cls() - verbosity = args.verbosity - if isinstance(verbosity, Sequence) and not isinstance(verbosity, (str, bytes)): - verbosity = verbosity[0] - atlas_entries: list[tuple[str, Path]] = list() for label, atlas_path in args.atlas or []: atlas_entries.append((label, Path(atlas_path).expanduser().resolve())) @@ -61,7 +57,7 @@ def from_cli_args(cls, args: argparse.Namespace | None) -> "WonkyconnConfig": analysis_level=args.analysis_level, phenotypes=_coerce_path(args.phenotypes), atlas=atlas_entries, - verbosity=int(verbosity) if verbosity is not None else 2, + log_level=args.log_level, debug=bool(args.debug), metrics=metrics, suppress_warnings=bool(args.suppress_warnings), diff --git a/wonkyconn/logger.py b/wonkyconn/logger.py index 5bcee33d..b54fd678 100644 --- a/wonkyconn/logger.py +++ b/wonkyconn/logger.py @@ -20,17 +20,3 @@ def _setup_logger(log_level: str = "INFO") -> logging.Logger: logger = _setup_logger() - - -def set_verbosity(verbosity: int | list[int]) -> None: - """Set the logger verbosity level (0=ERROR, 1=WARNING, 2=INFO, 3=DEBUG).""" - if isinstance(verbosity, list): - verbosity = verbosity[0] - if verbosity == 0: - logger.setLevel("ERROR") - elif verbosity == 1: - logger.setLevel("WARNING") - elif verbosity == 2: - logger.setLevel("INFO") - elif verbosity == 3: - logger.setLevel("DEBUG") diff --git a/wonkyconn/run.py b/wonkyconn/run.py index 26d00120..4aeb1e94 100644 --- a/wonkyconn/run.py +++ b/wonkyconn/run.py @@ -3,7 +3,7 @@ import argparse import sys from pathlib import Path -from typing import Sequence +from typing import Any, Sequence from . import __version__ from .config import WonkyconnConfig @@ -11,6 +11,17 @@ from .workflow import workflow +class VerbosityAction(argparse.Action): + def __call__( + self, + parser: argparse.ArgumentParser, + namespace: argparse.Namespace, + values: Any, + option_string: str | None = None, + ) -> None: + namespace.log_level = {0: "ERROR", 1: "WARNING", 2: "INFO", 3: "DEBUG"}[values] + + def global_parser(exit_on_error: bool = True) -> argparse.ArgumentParser: """Build the CLI argument parser for wonkyconn.""" parser = argparse.ArgumentParser( @@ -74,6 +85,21 @@ def global_parser(exit_on_error: bool = True) -> argparse.ArgumentParser: default=metrics, ) + logging_group = parser.add_mutually_exclusive_group(required=False) + logging_group.add_argument( + "--verbosity", + choices=[0, 1, 2, 3], + type=int, + action=VerbosityAction, + ) + logging_group.add_argument( + "--log-level", + required=False, + choices=["DEBUG", "INFO", "WARNING", "ERROR", "CRITICAL"], + default="INFO", + type=str, + ) + parser.add_argument("-v", "--version", action="version", version=__version__) parser.add_argument("--debug", action="store_true", default=False) parser.add_argument( @@ -83,17 +109,6 @@ def global_parser(exit_on_error: bool = True) -> argparse.ArgumentParser: default=False, help="Apply site correction to the data.", ) - parser.add_argument( - "--verbosity", - help=""" - Verbosity level. - """, - required=False, - choices=[0, 1, 2, 3], - default=2, - type=int, - nargs=1, - ) parser.add_argument( "--textual", action="store_true", @@ -157,28 +172,28 @@ def main(argv: None | Sequence[str] = None) -> None: if use_gui: _loosen_parser_for_gui(parser) - parsed_args: argparse.Namespace | None + args: argparse.Namespace | None try: - parsed_args = parser.parse_args(raw_args) + args = parser.parse_args(raw_args) except SystemExit as exc: if use_gui and exc.code != 0: - parsed_args = None + args = None else: raise except argparse.ArgumentError: if not use_gui: raise - parsed_args = None + args = None if use_gui: - initial_config = _build_initial_config(parsed_args) + initial_config = _build_initial_config(args) config = _run_textual_ui(initial_config) if config is None: return debug_enabled = config.debug suppress_warnings = config.suppress_warnings else: - config = _build_initial_config(parsed_args) + config = _build_initial_config(args) debug_enabled = config.debug suppress_warnings = config.suppress_warnings diff --git a/wonkyconn/textual_app.py b/wonkyconn/textual_app.py index cfe5f966..d61369cc 100644 --- a/wonkyconn/textual_app.py +++ b/wonkyconn/textual_app.py @@ -195,12 +195,12 @@ def compose(self) -> ComposeResult: with Horizontal(): yield Select( options=[ - ("Errors only (0)", "0"), - ("Warnings (1)", "1"), - ("Info (2)", "2"), - ("Debug (3)", "3"), + ("Errors only", "ERROR"), + ("Warnings and errors", "WARNING"), + ("Info", "INFO"), + ("Debug", "DEBUG"), ], - id="verbosity", + id="log_level", ) yield Checkbox("Debug", id="debug") yield Checkbox("Skip age/sex prediction and gradient", id="light_mode") @@ -289,7 +289,7 @@ def _load_initial_values(self) -> None: self._path_values["atlas_path"] = str(path) self.query_one("#atlas_path_display", Button).label = f"Atlas Path: {path}" - self.query_one("#verbosity", Select).value = str(self.initial_config.verbosity) + self.query_one("#log_level", Select).value = str(self.initial_config.log_level) self.query_one("#debug", Checkbox).value = bool(self.initial_config.debug) self.query_one("#light_mode", Checkbox).value = bool(self.initial_config.light_mode) self.query_one("#suppress_warnings", Checkbox).value = bool(self.initial_config.suppress_warnings) @@ -337,8 +337,8 @@ def _validate_and_build_config(self) -> WonkyconnConfig | None: elif not (atlas_path.exists() or atlas_path.is_symlink()) or (atlas_path.exists() and not atlas_path.is_file()): errors.append(f"Atlas file must exist: {atlas_path}") - raw_verbosity = self.query_one("#verbosity", Select[str]).value - verbosity = int("2" if isinstance(raw_verbosity, NoSelection) else raw_verbosity) + raw_log_level = self.query_one("#log_level", Select[str]).value + log_level = "INFO" if isinstance(raw_log_level, NoSelection) else raw_log_level debug = self.query_one("#debug", Checkbox).value light_mode = self.query_one("#light_mode", Checkbox).value suppress_warnings = self.query_one("#suppress_warnings", Checkbox).value @@ -354,7 +354,7 @@ def _validate_and_build_config(self) -> WonkyconnConfig | None: analysis_level="group", phenotypes=phenotypes, atlas=[(atlas_label, atlas_path)], - verbosity=verbosity, + log_level=log_level, debug=debug, metrics=light_mode_metrics if light_mode else all_metrics, theme="dark" if self.dark else "light", diff --git a/wonkyconn/workflow.py b/wonkyconn/workflow.py index 8fcf6e0a..9cacd260 100644 --- a/wonkyconn/workflow.py +++ b/wonkyconn/workflow.py @@ -30,7 +30,7 @@ calculate_qcfc_percentage, ) from .file_index.bids import BIDSIndex -from .logger import logger, set_verbosity +from .logger import logger from .visualization.plot import plot _DEFAULT_N_SPLITS = 20 @@ -53,8 +53,8 @@ def is_halfpipe(index: BIDSIndex) -> bool: def workflow(config: WonkyconnConfig) -> None: """Run the group-level connectivity quality-control pipeline.""" - if "pytest" not in sys.modules: - set_verbosity(config.verbosity) + if "pytest" not in sys.modules and config.log_level: + logger.setLevel(config.log_level) logger.debug(vars(config)) # Check BIDS path From b7cb04baa5ea208e1e42b04df682f26b80ac1e42 Mon Sep 17 00:00:00 2001 From: Lea Waller Date: Mon, 10 Aug 2026 19:39:52 +0200 Subject: [PATCH 32/33] Reduce the amount of log messages --- wonkyconn/workflow.py | 22 ++++++++++++++++------ 1 file changed, 16 insertions(+), 6 deletions(-) diff --git a/wonkyconn/workflow.py b/wonkyconn/workflow.py index 9cacd260..854d2cc8 100644 --- a/wonkyconn/workflow.py +++ b/wonkyconn/workflow.py @@ -95,7 +95,7 @@ def workflow(config: WonkyconnConfig) -> None: grouped_connectivity_matrix: defaultdict[tuple[str, ...], list[ConnectivityMatrix]] = defaultdict(list) segs: set[str] = set() - for timeseries_path in index.get(suffix="timeseries", extension=".tsv"): + for timeseries_path in tqdm(index.get(suffix="timeseries", extension=".tsv"), unit="connectivity matrices"): query = dict(**index.get_tags(timeseries_path)) del query["suffix"] @@ -171,7 +171,9 @@ def make_record( site_correction: bool = False, ) -> dict[str, Any]: """Compute all QC metrics for a single group of connectivity matrices.""" - seg_subjects: list[str] = list() + subjects: list[str] = list() + missing_subjects: list[str] = list() + filtered: list[ConnectivityMatrix] = list() for c in connectivity_matrices: @@ -185,14 +187,22 @@ def make_record( found = next((s for s in candidates if s in data_frame.index), None) if found: - seg_subjects.append(found) + subjects.append(found) filtered.append(c) else: - logger.info(f"Skipping subject {sub}: not found in phenotype file.") + missing_subjects.append(sub) + + if missing_subjects: + logger.info( + f"Found {len(subjects)} subjects, skipped {len(missing_subjects)} subjects not in " + f"phenotype file: {', '.join(missing_subjects)}" + ) + else: + logger.info(f"Found {len(subjects)} subjects") # Renaming for consistency connectivity_matrices[:] = filtered - seg_data_frame = data_frame.loc[seg_subjects] + seg_data_frame = data_frame.loc[subjects] (seg,) = index.get_tag_values(seg_key, {c.path for c in connectivity_matrices}) distance_matrix = distance_matrices[seg] @@ -287,7 +297,7 @@ def load_data_frame(config: WonkyconnConfig) -> pd.DataFrame: if "age" not in data_frame.columns: raise ValueError('Phenotypes file is missing the "age" column') if config.site_correction: - logger.info("Site correction is enabled - checking for 'site' column in phenotypes file.") if "site" not in data_frame.columns: raise ValueError('Phenotypes file is missing the "site" column required for site correction') + logger.info("Site correction is enabled") return data_frame From 3af2f83c5045d5bb05c61505dbf9e2aba5b8e41b Mon Sep 17 00:00:00 2001 From: Lea Waller Date: Mon, 10 Aug 2026 19:48:11 +0200 Subject: [PATCH 33/33] Store the dmn similarity results in a subdirectory Previously, each group would be assigned the same file name, leading to each subsequent group overwriting the previous group's results. Now, each group will have its own file name. --- wonkyconn/workflow.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/wonkyconn/workflow.py b/wonkyconn/workflow.py index 854d2cc8..d2a78c28 100644 --- a/wonkyconn/workflow.py +++ b/wonkyconn/workflow.py @@ -130,10 +130,9 @@ def workflow(config: WonkyconnConfig) -> None: records: list[dict[str, Any]] = list() for key, connectivity_matrices in tqdm(grouped_connectivity_matrix.items(), unit="groups"): - if len(group_by) == 2: - dmn_similarity_path = output_dir / f"dmn_similarity_{'-'.join(group_by)}.tsv" - else: - dmn_similarity_path = output_dir / f"dmn_similarity_{group_by[0]}.tsv" + suffix = "_".join(f"{key}-{value}" for key, value in zip(group_by, key, strict=False)) + dmn_similarity_path = output_dir / "dmn-similarity" / f"{suffix}.tsv" + record = make_record( index, data_frame, @@ -247,7 +246,10 @@ def make_record( if "analytic-insights" in metrics: gcor = calculate_gcor(connectivity_matrices) dmn_similarity_summary, t_stats_dmn_vis_fpn = network_similarity(connectivity_matrices, region_memberships[seg]) + + dmn_similarity_path.parent.mkdir(parents=True, exist_ok=True) dmn_similarity_summary.to_csv(dmn_similarity_path, sep="\t") + record.update( gcor=gcor, dmn_similarity_std=dmn_similarity_summary.loc[:, "corr_with_dmn"].std(),