From c85e01625ed34e7172dcb7d0963d60b5565327d8 Mon Sep 17 00:00:00 2001 From: Shuo Zhou Date: Wed, 3 Jun 2026 21:30:34 +0100 Subject: [PATCH 1/6] init examples --- examples/cmri_mpca/README.md | 31 ++ examples/cmri_mpca/__init__.py | 0 examples/cmri_mpca/config.py | 71 +++ examples/cmri_mpca/configs/README.md | 3 + examples/cmri_mpca/configs/tutorial_lr.yaml | 6 + examples/cmri_mpca/configs/tutorial_svc.yaml | 6 + examples/cmri_mpca/main.py | 157 +++++++ examples/cmri_mpca/tutorial.ipynb | 403 ++++++++++++++++++ examples/multisite_neuroimg_adapt/README.md | 42 ++ examples/multisite_neuroimg_adapt/__init__.py | 0 examples/multisite_neuroimg_adapt/config.py | 43 ++ .../configs/README.md | 3 + .../configs/tutorial.yaml | 7 + examples/multisite_neuroimg_adapt/main.py | 97 +++++ .../multisite_neuroimg_adapt/tutorial.ipynb | 293 +++++++++++++ examples/toy_domain_adaptation/README.md | 25 ++ examples/toy_domain_adaptation/__init__.py | 0 examples/toy_domain_adaptation/main.py | 98 +++++ examples/toy_domain_adaptation/tutorial.ipynb | 218 ++++++++++ setup.py | 1 + 20 files changed, 1504 insertions(+) create mode 100644 examples/cmri_mpca/README.md create mode 100644 examples/cmri_mpca/__init__.py create mode 100644 examples/cmri_mpca/config.py create mode 100644 examples/cmri_mpca/configs/README.md create mode 100644 examples/cmri_mpca/configs/tutorial_lr.yaml create mode 100644 examples/cmri_mpca/configs/tutorial_svc.yaml create mode 100644 examples/cmri_mpca/main.py create mode 100644 examples/cmri_mpca/tutorial.ipynb create mode 100644 examples/multisite_neuroimg_adapt/README.md create mode 100644 examples/multisite_neuroimg_adapt/__init__.py create mode 100644 examples/multisite_neuroimg_adapt/config.py create mode 100644 examples/multisite_neuroimg_adapt/configs/README.md create mode 100644 examples/multisite_neuroimg_adapt/configs/tutorial.yaml create mode 100644 examples/multisite_neuroimg_adapt/main.py create mode 100644 examples/multisite_neuroimg_adapt/tutorial.ipynb create mode 100644 examples/toy_domain_adaptation/README.md create mode 100644 examples/toy_domain_adaptation/__init__.py create mode 100644 examples/toy_domain_adaptation/main.py create mode 100644 examples/toy_domain_adaptation/tutorial.ipynb diff --git a/examples/cmri_mpca/README.md b/examples/cmri_mpca/README.md new file mode 100644 index 0000000..ecb1117 --- /dev/null +++ b/examples/cmri_mpca/README.md @@ -0,0 +1,31 @@ +# PAH Diagnosis from Cardiac MRI via Multilinear PCA + +## 1. Description + +This example demonstrates the multilinear PCA-based machine learning pipeline for cardiac MRI analysis [1], with application in pulmonary arterial hypertension (PAH) diagnosis. + +**Reference:** + +[1] Swift, A. J., Lu, H., Uthoff, J., Garg, P., Cogliano, M., Taylor, J., ... & Kiely, D. G. (2021). [A machine learning cardiac magnetic resonance approach to extract disease features and automate pulmonary arterial hypertension diagnosis](https://academic.oup.com/ehjcimaging/article/22/2/236/5717931). European Heart Journal-Cardiovascular Imaging. + +## 2. Usage + +* Datasets: [ShefPAH-179 v2.0 (short-axis) cardiac MRI dataset](https://github.com/pykale/data/tree/main/images/ShefPAH-179) +* Algorithms: MPCA + Linear SVM / Kernel SVM / Logistic Regression,... +* Example: Classification using SVM + +`python main.py --cfg configs/tutorial_svc.yaml` + +## 3. Related `kale` API + +`kale.interpret.model_weights`: Get model weights for interpretation. + +`kale.interpret.visualize`: Plot model weights or images. + +`kale.loaddata.image_access`: Load DICOM images as ndarray data. + +`kale.pipeline.mpca_trainer`: Pipeline of MPCA + feature selection + classification. + +`kale.prepdata.image_transform`: CMR images pre-processing. + +`kale.utils.download`: Download data. diff --git a/examples/cmri_mpca/__init__.py b/examples/cmri_mpca/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/examples/cmri_mpca/config.py b/examples/cmri_mpca/config.py new file mode 100644 index 0000000..c2e41e8 --- /dev/null +++ b/examples/cmri_mpca/config.py @@ -0,0 +1,71 @@ +""" +Default configurations for cardiac MRI data (ShefPAH) processing and classification +""" + +from yacs.config import CfgNode + +# ----------------------------------------------------------------------------- +# Config definition +# ----------------------------------------------------------------------------- + +_C = CfgNode() + +# ----------------------------------------------------------------------------- +# Dataset +# ----------------------------------------------------------------------------- +_C.DATASET = CfgNode() +_C.DATASET.SOURCE = "https://github.com/pykale/data/raw/main/images/ShefPAH-179/SA_64x64_v2.0.zip" +_C.DATASET.ROOT = "../data" +_C.DATASET.IMG_DIR = "DICOM" +_C.DATASET.BASE_DIR = "SA_64x64_v2.0" +_C.DATASET.FILE_FORAMT = "zip" +_C.DATASET.LANDMARK_FILE = "landmarks.csv" +_C.DATASET.MASK_DIR = "Mask" +# ---------------------------------------------------------------------------- # +# Image processing +# ---------------------------------------------------------------------------- # +_C.PROC = CfgNode() +_C.PROC.SCALE = 2 + +# ---------------------------------------------------------------------------- # +# Visualization +# ---------------------------------------------------------------------------- # +_C.PLT_KWS = CfgNode() +_C.PLT_KWS.PLT = CfgNode() +_C.PLT_KWS.PLT.n_cols = 10 + +_C.PLT_KWS.IM = CfgNode() +_C.PLT_KWS.IM.cmap = "gray" + +_C.PLT_KWS.MARKER = CfgNode() +_C.PLT_KWS.MARKER.marker = "+" +_C.PLT_KWS.MARKER.color = "r" +_C.PLT_KWS.MARKER.s = 100 +_C.PLT_KWS.MARKER.linewidths = 1.5 +_C.PLT_KWS.MARKER.edgecolors = "face" +# see https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.scatter.html for more options + +_C.PLT_KWS.WEIGHT = CfgNode() +_C.PLT_KWS.WEIGHT.markersize = 6 +_C.PLT_KWS.WEIGHT.alpha = 0.7 + +# ---------------------------------------------------------------------------- # +# Machine learning pipeline +# ---------------------------------------------------------------------------- # +_C.PIPELINE = CfgNode() +_C.PIPELINE.CLASSIFIER = "linear_svc" # ["svc", "linear_svc", "lr"] + +# ---------------------------------------------------------------------------- # +# Misc options +# ---------------------------------------------------------------------------- # +_C.OUTPUT = CfgNode() +_C.OUTPUT.OUT_DIR = "./outputs" # output_dir +_C.OUTPUT.SAVE_FIG = True + +_C.SAVE_FIG_KWARGS = CfgNode() +_C.SAVE_FIG_KWARGS.format = "pdf" +_C.SAVE_FIG_KWARGS.bbox_inches = "tight" + + +def get_cfg_defaults(): + return _C.clone() diff --git a/examples/cmri_mpca/configs/README.md b/examples/cmri_mpca/configs/README.md new file mode 100644 index 0000000..6a1ac03 --- /dev/null +++ b/examples/cmri_mpca/configs/README.md @@ -0,0 +1,3 @@ +# Configurations + +Configurations for experiments using [YAML](https://en.wikipedia.org/wiki/YAML). diff --git a/examples/cmri_mpca/configs/tutorial_lr.yaml b/examples/cmri_mpca/configs/tutorial_lr.yaml new file mode 100644 index 0000000..0b0ca4a --- /dev/null +++ b/examples/cmri_mpca/configs/tutorial_lr.yaml @@ -0,0 +1,6 @@ +DATASET: + LANDMARK_FILE: "landmarks_64x64.csv" +PIPELINE: + CLASSIFIER: "lr" +OUTPUT: + SAVE_FIG: True diff --git a/examples/cmri_mpca/configs/tutorial_svc.yaml b/examples/cmri_mpca/configs/tutorial_svc.yaml new file mode 100644 index 0000000..e90f2bd --- /dev/null +++ b/examples/cmri_mpca/configs/tutorial_svc.yaml @@ -0,0 +1,6 @@ +DATASET: + LANDMARK_FILE: "landmarks_64x64.csv" +PIPELINE: + CLASSIFIER: "svc" +OUTPUT: + SAVE_FIG: True diff --git a/examples/cmri_mpca/main.py b/examples/cmri_mpca/main.py new file mode 100644 index 0000000..f6014d3 --- /dev/null +++ b/examples/cmri_mpca/main.py @@ -0,0 +1,157 @@ +""" +PAH Diagnosis from Cardiac MRI via a Multilinear PCA-based Pipeline + +Reference: +Swift, A. J., Lu, H., Uthoff, J., Garg, P., Cogliano, M., Taylor, J., ... & Kiely, D. G. (2021). A machine learning +cardiac magnetic resonance approach to extract disease features and automate pulmonary arterial hypertension diagnosis. +European Heart Journal-Cardiovascular Imaging. https://academic.oup.com/ehjcimaging/article/22/2/236/5717931 +""" +import argparse +import os + +import numpy as np +import pandas as pd +from config import get_cfg_defaults +from kale.interpret import model_weights, visualize +from kale.loaddata.image_access import dicom2arraylist, read_dicom_dir +from kale.pipeline.mpca_trainer import MPCATrainer +from kale.prepdata.image_transform import mask_img_stack, normalize_img_stack, reg_img_stack, rescale_img_stack +from kale.utils.download import download_file_by_url +from sklearn.model_selection import cross_validate + + +def arg_parse(): + """Parsing arguments""" + parser = argparse.ArgumentParser(description="Machine learning pipeline for PAH diagnosis") + parser.add_argument("--cfg", required=True, help="path to config file", type=str) + + args = parser.parse_args() + return args + + +def main(): + args = arg_parse() + + # ---- setup configs ---- + cfg = get_cfg_defaults() + cfg.merge_from_file(args.cfg) + cfg.freeze() + print(cfg) + + save_figs = cfg.OUTPUT.SAVE_FIG + fig_format = cfg.SAVE_FIG_KWARGS.format + print(f"Save Figures: {save_figs}") + + # ---- initialize folder to store images ---- + save_figures_location = cfg.OUTPUT.OUT_DIR + print(f"Save Figures: {save_figures_location}") + + if not os.path.exists(save_figures_location): + os.makedirs(save_figures_location) + + # ---- setup dataset ---- + base_dir = cfg.DATASET.BASE_DIR + file_format = cfg.DATASET.FILE_FORAMT + download_file_by_url(cfg.DATASET.SOURCE, cfg.DATASET.ROOT, "%s.%s" % (base_dir, file_format), file_format) + + img_path = os.path.join(cfg.DATASET.ROOT, base_dir, cfg.DATASET.IMG_DIR) + patient_dcm_list = read_dicom_dir(img_path, sort_instance=True, sort_patient=True) + images, patient_ids = dicom2arraylist(patient_dcm_list, return_patient_id=True) + patient_ids = np.array(patient_ids, dtype=int) + n_samples = len(images) + + mask_path = os.path.join(cfg.DATASET.ROOT, base_dir, cfg.DATASET.MASK_DIR) + mask_dcm = read_dicom_dir(mask_path, sort_instance=True) + mask = dicom2arraylist(mask_dcm, return_patient_id=False)[0][0, ...] + + landmark_path = os.path.join(cfg.DATASET.ROOT, base_dir, cfg.DATASET.LANDMARK_FILE) + landmark_df = pd.read_csv(landmark_path, index_col="Subject").loc[patient_ids] # read .csv file as dataframe + landmarks = landmark_df.iloc[:, :-1].values + y = landmark_df["Group"].values + y[np.where(y != 0)] = 1 # convert to binary classification problem, i.e. no PH vs PAH + + # plot the first phase of images with landmarks + marker_names = list(landmark_df.columns[1::2]) + markers = [] + for marker in marker_names: + marker_name = marker.split(" ") + marker_name.pop(-1) + marker_name = " ".join(marker_name) + markers.append(marker_name) + + if save_figs: + n_img_per_fig = 45 + n_figures = int(n_samples / n_img_per_fig) + 1 + for k in range(n_figures): + visualize.plot_multi_images( + [images[i][0, ...] for i in range(k * n_img_per_fig, min((k + 1) * n_img_per_fig, n_samples))], + marker_locs=landmarks[k * n_img_per_fig : min((k + 1) * n_img_per_fig, n_samples), :], + im_kwargs=dict(cfg.PLT_KWS.IM), + marker_cmap="Set1", + marker_kwargs=dict(cfg.PLT_KWS.MARKER), + marker_titles=markers, + image_titles=list(patient_ids[k * n_img_per_fig : min((k + 1) * n_img_per_fig, n_samples)]), + n_cols=5, + ).savefig( + str(save_figures_location) + "/0)landmark_visualization_%s_of_%s.%s" % (k + 1, n_figures, fig_format), + **dict(cfg.SAVE_FIG_KWARGS), + ) + + # ---- data pre-processing ---- + # ----- image registration ----- + img_reg, max_dist = reg_img_stack(images.copy(), landmarks, landmarks[0]) + plt_kawargs = {**{"im_kwargs": dict(cfg.PLT_KWS.IM), "image_titles": list(patient_ids)}, **dict(cfg.PLT_KWS.PLT)} + if save_figs: + visualize.plot_multi_images([img_reg[i][0, ...] for i in range(n_samples)], **plt_kawargs).savefig( + str(save_figures_location) + "/1)image_registration.%s" % fig_format, **dict(cfg.SAVE_FIG_KWARGS) + ) + + # ----- masking ----- + img_masked = mask_img_stack(img_reg.copy(), mask) + if save_figs: + visualize.plot_multi_images([img_masked[i][0, ...] for i in range(n_samples)], **plt_kawargs).savefig( + str(save_figures_location) + "/2)masking.%s" % fig_format, **dict(cfg.SAVE_FIG_KWARGS) + ) + + # ----- resize ----- + img_rescaled = rescale_img_stack(img_masked.copy(), scale=1 / cfg.PROC.SCALE) + if save_figs: + visualize.plot_multi_images([img_rescaled[i][0, ...] for i in range(n_samples)], **plt_kawargs).savefig( + str(save_figures_location) + "/3)resize.%s" % fig_format, **dict(cfg.SAVE_FIG_KWARGS) + ) + + # ----- normalization ----- + img_norm = normalize_img_stack(img_rescaled.copy()) + if save_figs: + visualize.plot_multi_images([img_norm[i][0, ...] for i in range(n_samples)], **plt_kawargs).savefig( + str(save_figures_location) + "/4)normalize.%s" % fig_format, **dict(cfg.SAVE_FIG_KWARGS) + ) + + # ---- evaluating machine learning pipeline ---- + x = np.concatenate([img_norm[i].reshape((1,) + img_norm[i].shape) for i in range(n_samples)], axis=0) + trainer = MPCATrainer(classifier=cfg.PIPELINE.CLASSIFIER, n_features=200) + cv_results = cross_validate(trainer, x, y, cv=10, scoring=["accuracy", "roc_auc"], n_jobs=1) + + print("Averaged training time: {:.4f} seconds".format(np.mean(cv_results["fit_time"]))) + print("Averaged testing time: {:.4f} seconds".format(np.mean(cv_results["score_time"]))) + print("Averaged Accuracy: {:.4f}".format(np.mean(cv_results["test_accuracy"]))) + print("Averaged AUC: {:.4f}".format(np.mean(cv_results["test_roc_auc"]))) + + # ---- model weights interpretation ---- + trainer.fit(x, y) + + weights = trainer.mpca.inverse_transform(trainer.clf.coef_) - trainer.mpca.mean_ + weights = rescale_img_stack(weights, cfg.PROC.SCALE) # rescale weights to original shape + weights = mask_img_stack(weights, mask) # masking weights + top_weights = model_weights.select_top_weight(weights, select_ratio=0.02) # select top 2% weights + if save_figs: + visualize.plot_weights( + top_weights[0][0], + background_img=images[0][0], + im_kwargs=dict(cfg.PLT_KWS.IM), + marker_kwargs=dict(cfg.PLT_KWS.WEIGHT), + ).savefig(str(save_figures_location) + "/5)weights.%s" % fig_format, **dict(cfg.SAVE_FIG_KWARGS)) + + +if __name__ == "__main__": + main() diff --git a/examples/cmri_mpca/tutorial.ipynb b/examples/cmri_mpca/tutorial.ipynb new file mode 100644 index 0000000..059ab2a --- /dev/null +++ b/examples/cmri_mpca/tutorial.ipynb @@ -0,0 +1,403 @@ +{ + "nbformat": 4, + "nbformat_minor": 4, + "metadata": {}, + "cells": [ + { + "metadata": {}, + "source": [ + "# PyKale Tutorial: PAH Diagnosis from Cardiac MRI (CMR) via a Multilinear PCA-based Pipeline\n", + "| [Open in Colab](https://colab.research.google.com/github/pykale/pykale/blob/main/examples/cmri_mpca/tutorial.ipynb) (click `Runtime`\u2006\u2192\u2006`Run all (Ctrl+F9)` | [Launch Binder](https://mybinder.org/v2/gh/pykale/pykale/HEAD?filepath=examples%2Fcmri_mpca%2Ftutorial.ipynb) (click `Run`\u2006\u2192\u2006`Run All Cells`) |" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "## Overview\n", + "\n", + "- Pre-processing:\n", + " - Registration\n", + " - Masking\n", + " - Rescaling\n", + " - Normalization\n", + "- Machine learning pipeline:\n", + " - Multilinear principal component analysis\n", + " - Discriminative feature selection\n", + " - Linear classification model training \n", + "\n", + "**Reference:**\n", + "\n", + "Swift, A. J., Lu, H., Uthoff, J., Garg, P., Cogliano, M., Taylor, J., ... & Kiely, D. G. (2021). [A machine learning cardiac magnetic resonance approach to extract disease features and automate pulmonary arterial hypertension diagnosis](https://academic.oup.com/ehjcimaging/article/22/2/236/5717931). European Heart Journal-Cardiovascular Imaging." + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "## Setup" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "if 'google.colab' in str(get_ipython()):\n", + " print('Running on CoLab')\n", + " !pip uninstall --yes imgaug && pip uninstall --yes albumentations && pip install git+https://github.com/aleju/imgaug.git\n", + " !git clone https://github.com/pykale/pykale.git\n", + " %cd pykale\n", + " !pip install .[image,example] \n", + " %cd examples/cmri_mpca\n", + "else:\n", + " print('Not running on CoLab')" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "This imports required modules." + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "import os\n", + "\n", + "%matplotlib inline\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import pandas as pd\n", + "from config import get_cfg_defaults\n", + "\n", + "from kale.utils.download import download_file_by_url\n", + "from kale.loaddata.image_access import dicom2arraylist, read_dicom_dir\n", + "from kale.interpret import visualize" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "## Get CMR Images, Landmark Locations, and Labels" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "cfg_path = \"configs/tutorial_svc.yaml\" # Path to `.yaml` config file\n", + "\n", + "cfg = get_cfg_defaults()\n", + "cfg.merge_from_file(cfg_path)\n", + "cfg.freeze()\n", + "print(cfg)" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "### Download Data" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "base_dir = cfg.DATASET.BASE_DIR\n", + "file_format = cfg.DATASET.FILE_FORAMT\n", + "download_file_by_url(cfg.DATASET.SOURCE, cfg.DATASET.ROOT, \"%s.%s\" % (base_dir, file_format), file_format)" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "### Read DICOM Images" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "img_path = os.path.join(cfg.DATASET.ROOT, base_dir, cfg.DATASET.IMG_DIR)\n", + "patient_dcm_list = read_dicom_dir(img_path, sort_instance=True, sort_patient=True)\n", + "images, patient_ids = dicom2arraylist(patient_dcm_list, return_patient_id=True)\n", + "patient_ids = np.array(patient_ids, dtype=int)\n", + "n_samples = len(images)\n", + "\n", + "mask_path = os.path.join(cfg.DATASET.ROOT, base_dir, cfg.DATASET.MASK_DIR)\n", + "mask_dcm = read_dicom_dir(mask_path, sort_instance=True)\n", + "mask = dicom2arraylist(mask_dcm, return_patient_id=False)[0][0, ...]" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "### Read Landmarks and Get Labels" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "landmark_path = os.path.join(cfg.DATASET.ROOT, base_dir, cfg.DATASET.LANDMARK_FILE)\n", + "landmark_df = pd.read_csv(landmark_path, index_col=\"Subject\").loc[patient_ids] # read .csv file as dataframe\n", + "landmarks = landmark_df.iloc[:, :-1].values\n", + "y = landmark_df[\"Group\"].values\n", + "y[np.where(y != 0)] = 1 # convert to binary classification problem, i.e. no PH vs PAH" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "### Visualizing Data and Landmarks" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "# Get landmark names from column names\n", + "marker_names = list(landmark_df.columns[1::2])\n", + "markers = []\n", + "for marker in marker_names:\n", + " marker_name = marker.split(\" \")\n", + " marker_name.pop(-1)\n", + " marker_name = \" \".join(marker_name)\n", + " markers.append(marker_name)" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "# plot the first phase of images with landmarks\n", + "\n", + "n_img_per_fig = 35\n", + "n_figures = int(n_samples / n_img_per_fig) + 1\n", + "for k in range(n_figures):\n", + " visualize.plot_multi_images(\n", + " [images[i][0, ...] for i in range(k * n_img_per_fig, min((k + 1) * n_img_per_fig, n_samples))],\n", + " marker_locs=landmarks[k * n_img_per_fig: min((k + 1) * n_img_per_fig, n_samples), :],\n", + " im_kwargs=dict(cfg.PLT_KWS.IM),\n", + " marker_cmap=\"Set1\",\n", + " marker_kwargs=dict(cfg.PLT_KWS.MARKER),\n", + " marker_titles=markers,\n", + " image_titles=list(patient_ids[k * n_img_per_fig: min((k + 1) * n_img_per_fig, n_samples)]),\n", + " n_cols=5,\n", + " ).show()" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "## CMR Pre-processing" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "from kale.prepdata.image_transform import mask_img_stack, normalize_img_stack, reg_img_stack, rescale_img_stack" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "### Image Registration" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "img_reg, max_dist = reg_img_stack(images.copy(), landmarks, landmarks[0])" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "plt_kawargs = {**{\"im_kwargs\": dict(cfg.PLT_KWS.IM), \"image_titles\": list(patient_ids)}, **dict(cfg.PLT_KWS.PLT)}\n", + "visualize.plot_multi_images([img_reg[i][0, ...] for i in range(n_samples)], **plt_kawargs).show()" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "### Masking" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "img_masked = mask_img_stack(img_reg.copy(), mask)" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "visualize.plot_multi_images([img_masked[i][0, ...] for i in range(n_samples)], **plt_kawargs).show()" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "### Rescaling" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "img_rescaled = rescale_img_stack(img_masked.copy(), scale=1 / cfg.PROC.SCALE)" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "visualize.plot_multi_images([img_rescaled[i][0, ...] for i in range(n_samples)], **plt_kawargs).show()" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "### Normalization" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "img_norm = normalize_img_stack(img_rescaled.copy())" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "visualize.plot_multi_images([img_norm[i][0, ...] for i in range(n_samples)], **plt_kawargs).show()" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "## PAH Classification" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "from sklearn.model_selection import cross_validate\n", + "from kale.pipeline.mpca_trainer import MPCATrainer\n", + "\n", + "x = np.concatenate([img_norm[i].reshape((1,) + img_norm[i].shape) for i in range(n_samples)], axis=0)\n", + "trainer = MPCATrainer(classifier=cfg.PIPELINE.CLASSIFIER, n_features=200)\n", + "cv_results = cross_validate(trainer, x, y, cv=10, scoring=[\"accuracy\", \"roc_auc\"], n_jobs=1)" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "cv_results" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "print(\"Averaged training time: {:.4f} seconds\" .format(np.mean(cv_results['fit_time'])))\n", + "print(\"Averaged testing time: {:.4f} seconds\" .format(np.mean(cv_results['score_time'])))\n", + "print(\"Averaged Accuracy: {:.4f}\" .format(np.mean(cv_results[\"test_accuracy\"])))\n", + "print(\"Averaged AUC: {:.4f}\" .format(np.mean(cv_results[\"test_roc_auc\"])))" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "## Model Interpretation" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "from kale.interpret import model_weights\n", + "\n", + "trainer.fit(x, y)\n", + "\n", + "weights = trainer.mpca.inverse_transform(trainer.clf.coef_) - trainer.mpca.mean_\n", + "weights = rescale_img_stack(weights, cfg.PROC.SCALE) # rescale weights to original shape\n", + "weights = mask_img_stack(weights, mask) # masking weights\n", + "top_weights = model_weights.select_top_weight(weights, select_ratio=0.02) # select top 2% weights\n", + "visualize.plot_weights(\n", + " top_weights[0][0],\n", + " background_img=images[0][0],\n", + " im_kwargs=dict(cfg.PLT_KWS.IM),\n", + " marker_kwargs=dict(cfg.PLT_KWS.WEIGHT),\n", + ").show()" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + } + ] +} diff --git a/examples/multisite_neuroimg_adapt/README.md b/examples/multisite_neuroimg_adapt/README.md new file mode 100644 index 0000000..a693a03 --- /dev/null +++ b/examples/multisite_neuroimg_adapt/README.md @@ -0,0 +1,42 @@ +# Autism Detection: Domain Adaptation for Multi-Site Neuroimaging Data Analysis + +### 1. Description + +This example demonstrates multi-source domain adaptation method with application in neuroimaging data analysis for +autism detection. + +### 2. Materials and Methods + +- Data: Four largest subsets of ABIDE I + +| Site | Number of Samples | +|--------|-------------------| +| NYU | 175 | +| UM_1 | 106 | +| UCLA_1 | 72 | +| USM | 71 | +- Atlas: CC200 +- Pre-processing pipeline: cpac +- Classification problem: Autism vs Control +- Pipeline: + 1. Constructing brain networks from resting-state fMRI data + 2. Classification with Ridge Classifier or Covariate Independence Regularized Least Squares (CoIRLS) classifier + + +### 3. Related `kale` API + +`kale.interpret.visualize`: Visualize the results of a model. + +`kale.pipeline.multi_domain_adapter.CoIRLS`: Covariate Independence Regularized Least Squares (CoIRLS) classifier. + +`kale.utils.download.download_file_by_url`: Download a file from a URL. + +### References + +[1] Craddock C., Benhajali Y., Chu C., Chouinard F., Evans A., Jakab A., Khundrakpam BS., Lewis JD., Li Q., Milham M., Yan C. and Bellec P. (2013). [The Neuro Bureau Preprocessing Initiative: Open Sharing of Preprocessed Neuroimaging Data and Derivatives](https://doi.org/10.3389/conf.fninf.2013.09.00041). Frontiers in Neuroinformatics, 7. + +[2] Abraham A., Pedregosa F., Eickenberg M., Gervais P., Mueller A., Kossaifi J., Gramfort A., Thirion B. and Varoquaux G. (2014). [Machine Learning for Neuroimaging with scikit-learn](https://doi.org/10.3389/fninf.2014.00014). Frontiers in Neuroinformatics, 8. + +[3] Zhou S., Li W., Cox C. and Lu H. (2020). [Side Information Dependence as a Regularizer for Analyzing Human Brain Conditions across Cognitive Experiments](https://doi.org/10.1609/aaai.v34i04.6179). Proceedings of the AAAI Conference on Artificial Intelligence, 34(04), 6957-6964. + +[4] Zhou S. (2022). [Interpretable Domain-Aware Learning for Neuroimage Classification](https://etheses.whiterose.ac.uk/31044/1/PhD_thesis_ShuoZhou_170272834.pdf) (Doctoral Dissertation, University of Sheffield). diff --git a/examples/multisite_neuroimg_adapt/__init__.py b/examples/multisite_neuroimg_adapt/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/examples/multisite_neuroimg_adapt/config.py b/examples/multisite_neuroimg_adapt/config.py new file mode 100644 index 0000000..6bcb2a6 --- /dev/null +++ b/examples/multisite_neuroimg_adapt/config.py @@ -0,0 +1,43 @@ +""" +Default configurations for classification on resting-state fMRI of ABIDE +""" + +from yacs.config import CfgNode + +# ----------------------------------------------------------------------------- +# Config definition +# ----------------------------------------------------------------------------- + +_C = CfgNode() + +# ----------------------------------------------------------------------------- +# Dataset +# ----------------------------------------------------------------------------- +_C.DATASET = CfgNode() +_C.DATASET.ROOT = "../data" +_C.DATASET.PIPELINE = "cpac" # options: {‘cpac’, ‘css’, ‘dparsf’, ‘niak’} +_C.DATASET.ATLAS = "rois_cc200" +# options: {rois_aal, rois_cc200, rois_cc400, rois_dosenbach160, rois_ez, rois_ho, rois_tt} +_C.DATASET.SITE_IDS = None # list of site ids to use, if None, use all sites +_C.DATASET.TARGET = "NYU" # target site ids, e.g. "UM_1", "UCLA_1", "USM" +# --------------------------------------------------------- +# Solver +# --------------------------------------------------------- +_C.SOLVER = CfgNode() +_C.SOLVER.SEED = 2023 +# ---------------------------------------------------------------------------- # +# Machine learning pipeline +# ---------------------------------------------------------------------------- # +_C.MODEL = CfgNode() +_C.MODEL.KERNEL = "rbf" +_C.MODEL.ALPHA = 0.01 +_C.MODEL.LAMBDA_ = 1.0 +# ---------------------------------------------------------------------------- # +# Misc options +# ---------------------------------------------------------------------------- # +_C.OUTPUT = CfgNode() +_C.OUTPUT.OUT_DIR = "./outputs" # output_dir + + +def get_cfg_defaults(): + return _C.clone() diff --git a/examples/multisite_neuroimg_adapt/configs/README.md b/examples/multisite_neuroimg_adapt/configs/README.md new file mode 100644 index 0000000..6a1ac03 --- /dev/null +++ b/examples/multisite_neuroimg_adapt/configs/README.md @@ -0,0 +1,3 @@ +# Configurations + +Configurations for experiments using [YAML](https://en.wikipedia.org/wiki/YAML). diff --git a/examples/multisite_neuroimg_adapt/configs/tutorial.yaml b/examples/multisite_neuroimg_adapt/configs/tutorial.yaml new file mode 100644 index 0000000..0a220b5 --- /dev/null +++ b/examples/multisite_neuroimg_adapt/configs/tutorial.yaml @@ -0,0 +1,7 @@ +MODEL: + KERNEL: "linear" + ALPHA: 1.0 + LAMBDA_: 1.0 +DATASET: + ATLAS: "rois_cc200" + SITE_IDS: ['NYU', "UM_1", "UCLA_1", "USM"] diff --git a/examples/multisite_neuroimg_adapt/main.py b/examples/multisite_neuroimg_adapt/main.py new file mode 100644 index 0000000..e97b738 --- /dev/null +++ b/examples/multisite_neuroimg_adapt/main.py @@ -0,0 +1,97 @@ +""" +Autism Detection: Domain Adaptation for Multi-Site Neuroimaging Data Analysis + +Reference: +[1] Craddock C., Benhajali Y., Chu C., Chouinard F., Evans A., Jakab A., Khundrakpam BS., Lewis JD., Li Q., Milham M., Yan C. and Bellec P. (2013). The Neuro Bureau Preprocessing Initiative: Open Sharing of Preprocessed Neuroimaging Data and Derivatives. Frontiers in Neuroinformatics, 7. https://doi.org/10.3389/conf.fninf.2013.09.00041 + +[2] Abraham A., Pedregosa F., Eickenberg M., Gervais P., Mueller A., Kossaifi J., Gramfort A., Thirion B. and Varoquaux G. (2014). Machine Learning for Neuroimaging with scikit-learn. Frontiers in Neuroinformatics, 8. https://doi.org/10.3389/fninf.2014.00014 + +[3] Zhou S., Li W., Cox C. and Lu H. (2020). Side Information Dependence as a Regularizer for Analyzing Human Brain Conditions across Cognitive Experiments. Proceedings of the AAAI Conference on Artificial Intelligence, 34(04), 6957-6964. https://doi.org/10.1609/aaai.v34i04.6179 + +[4] Zhou S. (2022). Interpretable Domain-Aware Learning for Neuroimage Classification (Doctoral Dissertation, University of Sheffield). https://etheses.whiterose.ac.uk/31044/1/PhD_thesis_ShuoZhou_170272834.pdf +""" +import argparse +import os + +import kale.utils.seed as seed +import numpy as np +import pandas as pd +from config import get_cfg_defaults +from kale.evaluate import cross_validation +from nilearn.connectome import ConnectivityMeasure +from nilearn.datasets import fetch_abide_pcp +from sklearn.linear_model import RidgeClassifier + +from kalelinear.estimator import CoIRLS + + +def arg_parse(): + parser = argparse.ArgumentParser( + description="Autism Detection: Domain Adaptation for Multi-Site Neuroimaging Data Analysis" + ) + parser.add_argument("--cfg", required=True, help="path to config file", type=str) + args = parser.parse_args() + return args + + +def main(): + args = arg_parse() + + # ---- Set up configs ---- + cfg = get_cfg_defaults() + cfg.merge_from_file(args.cfg) + cfg.freeze() + seed.set_seed(cfg.SOLVER.SEED) + + # ---- Fetch ABIDE fMRI timeseries ---- + fetch_abide_pcp( + data_dir=cfg.DATASET.ROOT, + pipeline=cfg.DATASET.PIPELINE, + band_pass_filtering=True, + global_signal_regression=False, + derivatives=cfg.DATASET.ATLAS, + quality_checked=False, + SITE_ID=cfg.DATASET.SITE_IDS, + verbose=1, + ) + + # ---- Read Phenotypic data ---- + pheno_file = os.path.join(cfg.DATASET.ROOT, "ABIDE_pcp/Phenotypic_V1_0b_preprocessed1.csv") + pheno_info = pd.read_csv(pheno_file, index_col=0) + + # ---- Read timeseries from files ---- + data_dir = os.path.join(cfg.DATASET.ROOT, "ABIDE_pcp/%s/filt_noglobal" % cfg.DATASET.PIPELINE) + use_idx = [] + time_series = [] + for i in pheno_info.index: + data_file_name = "%s_%s.1D" % (pheno_info.loc[i, "FILE_ID"], cfg.DATASET.ATLAS) + data_path = os.path.join(data_dir, data_file_name) + if os.path.exists(data_path): + time_series.append(np.loadtxt(data_path, skiprows=0)) + use_idx.append(i) + + # ---- Use "DX_GROUP" (autism vs control) as labels, and "SITE_ID" as covariates ---- + pheno = pheno_info.loc[use_idx, ["SITE_ID", "DX_GROUP"]].reset_index(drop=True) + + # ---- Extracting Brain Networks Features ---- + correlation_measure = ConnectivityMeasure(kind="correlation", vectorize=True) + brain_networks = correlation_measure.fit_transform(time_series) + + # ---- Machine Learning for Multi-site Data ---- + print("Baseline") + estimator = RidgeClassifier() + results = cross_validation.leave_one_group_out( + brain_networks, pheno["DX_GROUP"].values, pheno["SITE_ID"].values, estimator + ) + print(pd.DataFrame.from_dict(results)) + + print("Domain Adaptation") + estimator = CoIRLS(kernel=cfg.MODEL.KERNEL, lambda_=cfg.MODEL.LAMBDA_, alpha=cfg.MODEL.ALPHA) + results = cross_validation.leave_one_group_out( + brain_networks, pheno["DX_GROUP"].values, pheno["SITE_ID"].values, estimator, use_domain_adaptation=True + ) + print(pd.DataFrame.from_dict(results)) + + +if __name__ == "__main__": + main() diff --git a/examples/multisite_neuroimg_adapt/tutorial.ipynb b/examples/multisite_neuroimg_adapt/tutorial.ipynb new file mode 100644 index 0000000..3ced871 --- /dev/null +++ b/examples/multisite_neuroimg_adapt/tutorial.ipynb @@ -0,0 +1,293 @@ +{ + "nbformat": 4, + "nbformat_minor": 5, + "metadata": {}, + "cells": [ + { + "metadata": {}, + "source": [ + "# PyKale Tutorial: Domain Adaptation for Autism Detection with Multi-site Brain Imaging Data\n", + "| [Open in Colab](https://colab.research.google.com/github/pykale/pykale/blob/main/examples/multisite_neuroimg_adapt/tutorial.ipynb) (click `Runtime`\u2006\u2192\u2006`Run all (Ctrl+F9)` | [Launch Binder](https://mybinder.org/v2/gh/pykale/pykale/HEAD?filepath=examples%2Fmultisite_neuroimg_adapt%2Ftutorial.ipynb) (click `Run`\u2006\u2192\u2006`Run All Cells`) |" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "## Overview\n", + "\n", + "- Pre-processing:\n", + " - [Data loading](#Data-Preparation)\n", + " - [Construct brain networks](#Extracting-Brain-Networks-Features)\n", + "- Machine learning pipeline:\n", + " - [Baseline: Ridge classifier](#Baseline-Model)\n", + " - [Domain adaptation](#Domain-Adaptation)\n", + "\n", + "**Reference:**\n", + "\n", + "[1] Craddock C., Benhajali Y., Chu C., Chouinard F., Evans A., Jakab A., Khundrakpam BS., Lewis JD., Li Q., Milham M., Yan C. and Bellec P. (2013). [The Neuro Bureau Preprocessing Initiative: Open Sharing of Preprocessed Neuroimaging Data and Derivatives](https://doi.org/10.3389/conf.fninf.2013.09.00041). Frontiers in Neuroinformatics, 7.\n", + "\n", + "[2] Abraham A., Pedregosa F., Eickenberg M., Gervais P., Mueller A., Kossaifi J., Gramfort A., Thirion B. and Varoquaux G. (2014). [Machine Learning for Neuroimaging with scikit-learn](https://doi.org/10.3389/fninf.2014.00014). Frontiers in Neuroinformatics, 8.\n", + "\n", + "[3] Zhou S., Li W., Cox C. and Lu H. (2020). [Side Information Dependence as a Regularizer for Analyzing Human Brain Conditions across Cognitive Experiments](https://doi.org/10.1609/aaai.v34i04.6179). Proceedings of the AAAI Conference on Artificial Intelligence, 34(04), 6957-6964.\n", + "\n", + "[4] Zhou S. (2022). [Interpretable Domain-Aware Learning for Neuroimage Classification](https://etheses.whiterose.ac.uk/31044/1/PhD_thesis_ShuoZhou_170272834.pdf) (Doctoral Dissertation, University of Sheffield)." + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "## Setup" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "if 'google.colab' in str(get_ipython()):\n", + " print('Running on CoLab')\n", + " !pip uninstall --yes imgaug && pip uninstall --yes albumentations && pip install git+https://github.com/aleju/imgaug.git\n", + " !git clone https://github.com/pykale/pykale.git\n", + " %cd pykale\n", + " !pip install .[image,example]\n", + " %cd examples/multisite_neuroimg_adapt\n", + "else:\n", + " print('Not running on CoLab')" + ], + "cell_type": "code", + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "Not running on CoLab\n" + ] + } + ], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "This imports required modules." + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "import os\n", + "import numpy as np\n", + "import pandas as pd\n", + "from config import get_cfg_defaults\n", + "from nilearn.connectome import ConnectivityMeasure\n", + "from nilearn.datasets import fetch_abide_pcp\n", + "from sklearn.linear_model import RidgeClassifier\n", + "\n", + "import kale.utils.seed as seed\n", + "from kale.evaluate import cross_validation\n", + "from kalelinear.estimator import CoIRLS" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "# Path to `.yaml` config file\n", + "cfg_path = \"configs/tutorial.yaml\" \n", + "cfg = get_cfg_defaults()\n", + "cfg.merge_from_file(cfg_path)\n", + "cfg.freeze()\n", + "seed.set_seed(cfg.SOLVER.SEED)\n", + "print(cfg)" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "## Data Preparation\n", + "\n", + "### Fetch ABIDE fMRI timeseries" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "root_dir = cfg.DATASET.ROOT\n", + "pipeline = cfg.DATASET.PIPELINE # fmri pre-processing pipeline\n", + "atlas = cfg.DATASET.ATLAS\n", + "site_ids = cfg.DATASET.SITE_IDS\n", + "abide = fetch_abide_pcp(data_dir=root_dir, pipeline=pipeline,\n", + " band_pass_filtering=True, global_signal_regression=False,\n", + " derivatives=atlas, quality_checked=False,\n", + " SITE_ID=site_ids,\n", + " verbose=0)" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "### Read Phenotypic data" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "pheno_file = os.path.join(cfg.DATASET.ROOT, \"ABIDE_pcp/Phenotypic_V1_0b_preprocessed1.csv\")\n", + "pheno_info = pd.read_csv(pheno_file, index_col=0)" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "View Phenotypic data" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "pheno_info.head()" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "### Read timeseries from files" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "data_dir = os.path.join(root_dir, \"ABIDE_pcp/%s/filt_noglobal\" % pipeline)\n", + "use_idx = []\n", + "time_series = []\n", + "for i in pheno_info.index:\n", + " data_file_name = \"%s_%s.1D\" % (pheno_info.loc[i, \"FILE_ID\"], atlas)\n", + " data_path = os.path.join(data_dir, data_file_name)\n", + " if os.path.exists(data_path):\n", + " time_series.append(np.loadtxt(data_path, skiprows=0))\n", + " use_idx.append(i)" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "Use \"DX_GROUP\" (autism vs control) as labels, and \"SITE_ID\" as covariates" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "pheno = pheno_info.loc[use_idx, [\"SITE_ID\", \"DX_GROUP\"]].reset_index(drop=True)" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "## Extracting Brain Networks Features" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "correlation_measure = ConnectivityMeasure(kind='correlation', vectorize=True)\n", + "brain_networks = correlation_measure.fit_transform(time_series)" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "## Machine Learning for Multi-site Data\n", + "\n", + "### Cross validation Pipeline" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "### Baseline" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "estimator = RidgeClassifier()\n", + "results = cross_validation.leave_one_group_out(\n", + " brain_networks, pheno[\"DX_GROUP\"].values, pheno[\"SITE_ID\"].values, estimator\n", + ")" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "pd.DataFrame.from_dict(results)" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "### Domain Adaptation" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "estimator = CoIRLS(kernel=cfg.MODEL.KERNEL, lambda_=cfg.MODEL.LAMBDA_, alpha=cfg.MODEL.ALPHA)\n", + "results = cross_validation.leave_one_group_out(\n", + " brain_networks, pheno[\"DX_GROUP\"].values, pheno[\"SITE_ID\"].values, estimator, use_domain_adaptation=True\n", + ")" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "pd.DataFrame.from_dict(results)" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + } + ] +} diff --git a/examples/toy_domain_adaptation/README.md b/examples/toy_domain_adaptation/README.md new file mode 100644 index 0000000..ba5482b --- /dev/null +++ b/examples/toy_domain_adaptation/README.md @@ -0,0 +1,25 @@ +# Domain Adaptation on Toy Data + +### 1. Description + +This example demonstrates domain adaptation methods on the generated toy data. + +### 2. Usage + +* Datasets: Toy data (2 domains and 2 classes) generated by `sklearn.datasets.make_blobs`. +* Algorithms: Covariate Independence Regularized Least Squares (CoIRLS) classifier. + +`python main.py` + +### 3. Related `kale` API + +`kale.interpret.visualize.distplot_1d`: Visualize 1D distributions. + +`kale.pipeline.multi_domain_adapter.CoIRLS`: Covariate Independence Regularized Least Squares (CoIRLS) classifier. + + +### References + +[1] Zhou, S., Li, W., Cox, C.R., & Lu, H. (2020). [Side Information Dependence as a Regularizer for Analyzing Human Brain Conditions across Cognitive Experiments](https://ojs.aaai.org//index.php/AAAI/article/view/6179). in *AAAI 2020*, New York, USA. + +[2] Zhou, S. (2022). [Interpretable Domain-Aware Learning for Neuroimage Classification](https://etheses.whiterose.ac.uk/31044/1/PhD_thesis_ShuoZhou_170272834.pdf) (Doctoral dissertation, University of Sheffield). diff --git a/examples/toy_domain_adaptation/__init__.py b/examples/toy_domain_adaptation/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/examples/toy_domain_adaptation/main.py b/examples/toy_domain_adaptation/main.py new file mode 100644 index 0000000..ba21764 --- /dev/null +++ b/examples/toy_domain_adaptation/main.py @@ -0,0 +1,98 @@ +import matplotlib.pyplot as plt +import numpy as np +from kale.interpret.visualize import distplot_1d +from sklearn.datasets import make_blobs +from sklearn.linear_model import RidgeClassifier +from sklearn.metrics import accuracy_score +from sklearn.preprocessing import OneHotEncoder + +from kalelinear.estimator import CoIRLS + + +def main(): + np.random.seed(29118) + # Generate toy data + n_samples = 200 + + xs, ys = make_blobs(n_samples, centers=[[0, 0], [0, 2]], cluster_std=[0.3, 0.35]) + xt, yt = make_blobs(n_samples, centers=[[2, -2], [2, 0.2]], cluster_std=[0.35, 0.4]) + + # visualize toy data + colors = ["c", "m"] + x_all = [xs, xt] + y_all = [ys, yt] + labels = ["source", "Target"] + plt.figure(figsize=(8, 5)) + for i in range(2): + idx_pos = np.where(y_all[i] == 1) + idx_neg = np.where(y_all[i] == 0) + plt.scatter( + x_all[i][idx_pos, 0], + x_all[i][idx_pos, 1], + c=colors[i], + marker="o", + alpha=0.4, + label=labels[i] + " positive", + ) + plt.scatter( + x_all[i][idx_neg, 0], + x_all[i][idx_neg, 1], + c=colors[i], + marker="x", + alpha=0.4, + label=labels[i] + " negative", + ) + plt.legend() + plt.title("Source domain and target domain blobs data", fontsize=14, fontweight="bold") + plt.show() + + clf = RidgeClassifier(alpha=1.0) + clf.fit(xs, ys) + + yt_pred = clf.predict(xt) + print("Accuracy on target domain: {:.2f}".format(accuracy_score(yt, yt_pred))) + + # visualize decision scores of non-adaptation classifier + ys_score = clf.decision_function(xs) + yt_score = clf.decision_function(xt) + title = "Ridge classifier decision score distribution" + title_kwargs = {"fontsize": 14, "fontweight": "bold"} + hist_kwargs = {"kde": True, "alpha": 0.7} + plt_labels = ["Source", "Target"] + distplot_1d( + [ys_score, yt_score], + labels=plt_labels, + xlabel="Decision Scores", + title=title, + title_kwargs=title_kwargs, + hist_kwargs=hist_kwargs, + ).show() + + # domain adaptation + clf_ = CoIRLS(lambda_=1) + # encoding one-hot domain covariate matrix + covariates = np.zeros(n_samples * 2) + covariates[:n_samples] = 1 + enc = OneHotEncoder(handle_unknown="ignore") + covariates_mat = enc.fit_transform(covariates.reshape(-1, 1)).toarray() + + x = np.concatenate((xs, xt)) + clf_.fit(x, ys, covariates_mat) + yt_pred_ = clf_.predict(xt) + print("Accuracy on target domain: {:.2f}".format(accuracy_score(yt, yt_pred_))) + + ys_score_ = clf_.decision_function(xs).detach().numpy().reshape(-1) + yt_score_ = clf_.decision_function(xt).detach().numpy().reshape(-1) + title = "Domain adaptation classifier decision score distribution" + distplot_1d( + [ys_score_, yt_score_], + labels=plt_labels, + xlabel="Decision Scores", + title=title, + title_kwargs=title_kwargs, + hist_kwargs=hist_kwargs, + ).show() + + +if __name__ == "__main__": + main() diff --git a/examples/toy_domain_adaptation/tutorial.ipynb b/examples/toy_domain_adaptation/tutorial.ipynb new file mode 100644 index 0000000..dedad53 --- /dev/null +++ b/examples/toy_domain_adaptation/tutorial.ipynb @@ -0,0 +1,218 @@ +{ + "nbformat": 4, + "nbformat_minor": 5, + "metadata": {}, + "cells": [ + { + "metadata": {}, + "source": [ + "# PyKale Tutorial: Domain Adaptation on Toy Data\n", + "| [Open in Colab](https://colab.research.google.com/github/pykale/pykale/blob/main/examples/toy_domain_adaptation/tutorial.ipynb) (click `Runtime`\u2006\u2192\u2006`Run all (Ctrl+F9)` | [Launch Binder](https://mybinder.org/v2/gh/pykale/pykale/HEAD?filepath=examples%2Ftoy_domain_adaptation%2Ftutorial.ipynb) (click `Run`\u2006\u2192\u2006`Run All Cells`) |" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "### Setup" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "# import seaborn first to avoid seaborn import error caused by newer scipy version, to be solved later\n", + "import seaborn as sns" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "if 'google.colab' in str(get_ipython()):\n", + " print('Running on CoLab')\n", + " !pip uninstall --yes imgaug && pip uninstall --yes albumentations && pip install git+https://github.com/aleju/imgaug.git\n", + " !pip install numpy>=2.0.0\n", + " !pip install git+https://github.com/pykale/pykale.git\n", + " !git clone https://github.com/pykale/pykale.git\n", + " %cd pykale/examples/toy_domain_adaptation\n", + "else:\n", + " print('Not running on CoLab')" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "### Generate toy data" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "%matplotlib inline\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "\n", + "from sklearn.datasets import make_moons, make_blobs" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "n_samples = 1000\n", + "\n", + "xs, ys = make_blobs(n_samples, centers=[[0, 0], [0, 2]], cluster_std=[0.3, 0.35])\n", + "xt, yt = make_blobs(n_samples, centers=[[2, -2], [2, 0.2]], cluster_std=[0.35, 0.4])" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "colors = [\"c\", \"m\"]\n", + "x_all = [xs, xt]\n", + "y_all = [ys, yt]\n", + "labels = [\"source\", \"Target\"]\n", + "for i in range(2):\n", + " idx_pos = np.where(y_all[i] == 1)\n", + " idx_neg = np.where(y_all[i] == 0)\n", + " plt.scatter(x_all[i][idx_pos, 0], x_all[i][idx_pos, 1], c=colors[i], marker=\"o\", alpha=0.4, \n", + " label=labels[i] + \" positive\")\n", + " plt.scatter(x_all[i][idx_neg, 0], x_all[i][idx_neg, 1], c=colors[i], marker=\"x\", alpha=0.4, \n", + " label=labels[i] + \" negative\")\n", + "plt.legend()\n", + "plt.title('Source domain and target domain blobs data',fontsize=14,fontweight='bold')" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "### Classification" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "from sklearn.preprocessing import OneHotEncoder\n", + "from sklearn.linear_model import RidgeClassifier\n", + "from sklearn.metrics import accuracy_score\n", + "from kale.interpret.visualize import distplot_1d\n", + "from kalelinear.estimator import CoIRLS" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "#### Training a standard Ridge classifier" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "clf = RidgeClassifier(alpha=1.0)\n", + "clf.fit(xs, ys)\n", + "\n", + "yt_pred = clf.predict(xt)\n", + "print('Accuracy on target domain: {:.2f}'.format(accuracy_score(yt, yt_pred)))" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "ys_score = clf.decision_function(xs)\n", + "yt_score = clf.decision_function(xt)\n", + "\n", + "title = \"Ridge classifier decision score distribution\"\n", + "title_kwargs = {\"fontsize\": 14, \"fontweight\": \"bold\"}\n", + "hist_kwargs = {\"kde\": True, \"alpha\": 0.7}\n", + "plt_labels = [\"Source\", \"Target\"]\n", + "distplot_1d(\n", + " [ys_score, yt_score],\n", + " title=title,\n", + " xlabel=\"Decision Scores\",\n", + " labels=plt_labels,\n", + " hist_kwargs=hist_kwargs,\n", + " title_kwargs=title_kwargs,\n", + ").show()" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "#### Training a domain adaptation classifier" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "clf_ = CoIRLS()\n", + "# encoding one-hot domain covariate matrix\n", + "covariates = np.zeros(n_samples * 2)\n", + "covariates[:n_samples] = 1\n", + "enc = OneHotEncoder(handle_unknown=\"ignore\")\n", + "covariates_mat = enc.fit_transform(covariates.reshape(-1, 1)).toarray()" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "x = np.concatenate((xs, xt))\n", + "clf_.fit(x, ys, covariates_mat)\n", + "yt_pred_ = clf_.predict(xt)\n", + "print(\"Accuracy on target domain: {:.2f}\".format(accuracy_score(yt, yt_pred_)))" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "ys_score_ = clf_.decision_function(xs).reshape(-1)\n", + "yt_score_ = clf_.decision_function(xt).reshape(-1)\n", + "plt.figure(figsize=(10, 5))\n", + "title = \"Domain adaptation classifier decision score distribution\"\n", + "distplot_1d(\n", + " [ys_score_, yt_score_],\n", + " title=title,\n", + " xlabel=\"Decision Scores\",\n", + " labels=plt_labels,\n", + " hist_kwargs=hist_kwargs,\n", + " title_kwargs=title_kwargs,\n", + ").show()" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + } + ] +} diff --git a/setup.py b/setup.py index ea3008b..a318ede 100644 --- a/setup.py +++ b/setup.py @@ -27,6 +27,7 @@ "ipython", "matplotlib", "nilearn", + "pykale", "seaborn", "yacs>=0.1.7", ] From dbde43d27e78c99787d582b63b6fde4c2eb58646 Mon Sep 17 00:00:00 2001 From: Shuo Zhou Date: Fri, 5 Jun 2026 00:00:17 +0100 Subject: [PATCH 2/6] init gsda exmaple --- .../brain_lateralization/configs/__init__.py | 0 .../configs/default_cfg.py | 49 ++ .../configs/demo-gsp.yaml | 13 + .../configs/demo-hcp.yaml | 13 + .../brain_lateralization/configs/test.yaml | 7 + examples/brain_lateralization/tutorial.ipynb | 158 ++++ .../brain_lateralization/utils/__init__.py | 0 .../brain_lateralization/utils/experiment.py | 713 ++++++++++++++++++ examples/brain_lateralization/utils/io_.py | 527 +++++++++++++ examples/brain_lateralization/utils/plot.py | 201 +++++ 10 files changed, 1681 insertions(+) create mode 100644 examples/brain_lateralization/configs/__init__.py create mode 100644 examples/brain_lateralization/configs/default_cfg.py create mode 100644 examples/brain_lateralization/configs/demo-gsp.yaml create mode 100644 examples/brain_lateralization/configs/demo-hcp.yaml create mode 100644 examples/brain_lateralization/configs/test.yaml create mode 100644 examples/brain_lateralization/tutorial.ipynb create mode 100644 examples/brain_lateralization/utils/__init__.py create mode 100644 examples/brain_lateralization/utils/experiment.py create mode 100644 examples/brain_lateralization/utils/io_.py create mode 100644 examples/brain_lateralization/utils/plot.py diff --git a/examples/brain_lateralization/configs/__init__.py b/examples/brain_lateralization/configs/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/examples/brain_lateralization/configs/default_cfg.py b/examples/brain_lateralization/configs/default_cfg.py new file mode 100644 index 0000000..213f426 --- /dev/null +++ b/examples/brain_lateralization/configs/default_cfg.py @@ -0,0 +1,49 @@ +""" +Default configurations +""" + +from yacs.config import CfgNode as CN + +# ----------------------------------------------------------------------------- +# Config definition +# ----------------------------------------------------------------------------- + +_C = CN() + +# ----------------------------------------------------------------------------- +# Dataset +# ----------------------------------------------------------------------------- +_C.DATASET = CN() +_C.DATASET.DATASET = "HCP" +_C.DATASET.DOWNLOAD = True +_C.DATASET.ROOT = "data" +_C.DATASET.ATLAS = "BNA" +_C.DATASET.FEATURE = "correlation" +_C.DATASET.TEST_RATIO = 0.0 +_C.DATASET.TYPE = "functional" +_C.DATASET.SESSIONS = [None] +_C.DATASET.RUN = "Fisherz" +_C.DATASET.CONNECTION = "intra" +_C.DATASET.NUM_REPEAT = 1 +_C.DATASET.MIX_GROUP = False +_C.DATASET.REST1_ONLY = False # only use rest1 data for training, HCP dataset only +# ---------------------------------------------------------------------------- # +# Solver +# ---------------------------------------------------------------------------- # +_C.SOLVER = CN() +_C.SOLVER.SEED = 2023 +_C.SOLVER.L2_HPARAM = 0.1 # hyperparameter for the l2 regularization +_C.SOLVER.LR = 0.04 # Initial learning rate +_C.SOLVER.OPTIMIZER = "lbfgs" +_C.SOLVER.MAX_ITER = 100 +_C.SOLVER.LAMBDA_ = [0.0, 1.0, 2.0, 5.0, 8.0, 10.0] + +# ---------------------------------------------------------------------------- # +# Misc options +# ---------------------------------------------------------------------------- # +_C.OUTPUT = CN() +_C.OUTPUT.ROOT = "output" + + +def get_cfg_defaults(): + return _C.clone() diff --git a/examples/brain_lateralization/configs/demo-gsp.yaml b/examples/brain_lateralization/configs/demo-gsp.yaml new file mode 100644 index 0000000..34aefd7 --- /dev/null +++ b/examples/brain_lateralization/configs/demo-gsp.yaml @@ -0,0 +1,13 @@ +DATASET: + DATASET: "GSP" + ROOT: "data" + MIX_GROUP: False + TEST_RATIO: 0.2 +SOLVER: + LAMBDA_: [0.0, 2.0, 5.0, 10] + LR: 1.0 + SEED: 2023 + OPTIMIZER: "lbfgs" + MAX_ITER: 100 +OUTPUT: + ROOT: "output" diff --git a/examples/brain_lateralization/configs/demo-hcp.yaml b/examples/brain_lateralization/configs/demo-hcp.yaml new file mode 100644 index 0000000..efa7c92 --- /dev/null +++ b/examples/brain_lateralization/configs/demo-hcp.yaml @@ -0,0 +1,13 @@ +DATASET: + DATASET: "HCP" + ROOT: "data" + MIX_GROUP: False + TEST_RATIO: 0.2 +SOLVER: + LAMBDA_: [0.0, 2.0, 5.0, 10] + LR: 1.0 + SEED: 2023 + OPTIMIZER: "lbfgs" + MAX_ITER: 100 +OUTPUT: + ROOT: "output" diff --git a/examples/brain_lateralization/configs/test.yaml b/examples/brain_lateralization/configs/test.yaml new file mode 100644 index 0000000..c8c5e71 --- /dev/null +++ b/examples/brain_lateralization/configs/test.yaml @@ -0,0 +1,7 @@ +DATASET: + ROOT: "D:/ShareFolder/BNA/Proc" + NUM_REPEAT: 5 +SOLVER: + SEED: 2022 +OUTPUT: + ROOT: 'D:/ShareFolder/BNA/Result' diff --git a/examples/brain_lateralization/tutorial.ipynb b/examples/brain_lateralization/tutorial.ipynb new file mode 100644 index 0000000..2ac58fc --- /dev/null +++ b/examples/brain_lateralization/tutorial.ipynb @@ -0,0 +1,158 @@ +{ + "nbformat": 4, + "nbformat_minor": 2, + "metadata": {}, + "cells": [ + { + "metadata": {}, + "source": [ + "# Group-Specific Discriminant Analysis for sex-specific lateralization Running Demo\n", + "\n", + "[Open in Colab](https://colab.research.google.com/github/shuo-zhou/GSDA-Lateralization/blob/main/gsda_demo.ipynb) (click `Runtime`\u2006\u2192\u2006`Run all (Ctrl+F9)`)" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "## Setup\n", + "\n", + "The first few blocks of code are necessary to set up the notebook execution environment. This checks if the notebook is running on Google Colab and installs required packages." + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "#@title Setup environment and fetch code from GitHub\n", + "\n", + "if 'google.colab' in str(get_ipython()):\n", + " print('Running on CoLab')\n", + " !pip install yacs\n", + " !git clone https://github.com/pykale/kale-linear\n", + " %cd kale-linear/examples/brain_lateralization\n", + "else:\n", + " print('Not running on CoLab')\n", + " " + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "#@title Import required modules\n", + "import os\n", + "from configs.default_cfg import get_cfg_defaults\n", + "from utils.experiment import run_experiment\n", + "\n", + "from utils.io_ import load_result, reformat_results\n", + "from utils import plot" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "### Configurations\n", + "\n", + "The customized configuration used in this demo is stored in `configs/demoh-cp.yaml`, this file overwrites defaults in `default_cfg.py` where a value is specified. Change the configuration file path to `cfg_path = \"configs/demo-gsp.yaml\"` for running the demo with GSP data." + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "#@title Setup configs using .yaml file\n", + "cfg_path = \"configs/demo-hcp.yaml\" # Path to `.yaml` config file\n", + "\n", + "cfg = get_cfg_defaults()\n", + "cfg.merge_from_file(cfg_path)\n", + "cfg.freeze()\n", + "print(cfg)" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "## Group-specific model training\n", + "\n", + "It could take a while (15 to 25 mins) to run the experiments. " + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "_ = [run_experiment(cfg, lambda_) for lambda_ in cfg.SOLVER.LAMBDA_]" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "## Classification result visualization" + ], + "cell_type": "markdown" + }, + { + "metadata": {}, + "source": [ + "#@title Load results and reformat as a dataframe\n", + "\n", + "dataset = cfg.DATASET.DATASET\n", + "model_root_dir = cfg.OUTPUT.ROOT\n", + "lambdas = cfg.SOLVER.LAMBDA_\n", + "seed_start = cfg.SOLVER.SEED\n", + "test_size = cfg.DATASET.TEST_RATIO\n", + "\n", + "res_df = load_result(dataset=dataset, root_dir=model_root_dir, \n", + " lambdas=lambdas, seed_start=seed_start, test_size=test_size)\n", + "\n", + "res_df[\"GSI_train_session\"] = 2 * (res_df[\"acc_tgt_train_session\"] * \n", + " (res_df[\"acc_tgt_train_session\"] - \n", + " res_df[\"acc_nt_train_session\"]))\n", + "res_df[\"GSI_test_session\"] = 2 * (res_df[\"acc_tgt_test_session\"] * \n", + " (res_df[\"acc_tgt_test_session\"] - \n", + " res_df[\"acc_nt_test_session\"]))\n", + "\n", + "res_df_train_session = reformat_results(res_df, [\"acc_tgt_train_session\", \"acc_nt_train_session\"])\n", + "res_df_test_session = reformat_results(res_df, [\"acc_tgt_test_session\", \"acc_nt_test_session\"])" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "#@title Plot accuracy\n", + "\n", + "if not os.path.exists(\"figures\"):\n", + " os.mkdir(\"figures\")\n", + "plot.plot_accuracy(res_df_train_session)" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "source": [ + "#@title Plot Group Specificity Index (GSI)\n", + "plot.plot_gsi(res_df, x=\"lambda\", y=\"GSI_train_session\", hue=\"target_group\")" + ], + "cell_type": "code", + "outputs": [], + "execution_count": null + } + ] +} diff --git a/examples/brain_lateralization/utils/__init__.py b/examples/brain_lateralization/utils/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/examples/brain_lateralization/utils/experiment.py b/examples/brain_lateralization/utils/experiment.py new file mode 100644 index 0000000..ab15c6b --- /dev/null +++ b/examples/brain_lateralization/utils/experiment.py @@ -0,0 +1,713 @@ +import copy +import os +import sys + +# import pickle +import numpy as np +import pandas as pd +import torch +from sklearn.metrics import accuracy_score # , roc_auc_score +from sklearn.model_selection import StratifiedShuffleSplit +from torch.hub import download_url_to_file + +sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), "../.."))) +from utils.io_ import load_half_brain, pick_half, read_tabular + +from kalelinear.estimator import GSDA # , GSLRTorch + +BASE_RESULT_DICT: dict = { + "pred_loss": [], + "hsic_loss": [], + "lambda": [], + "train_session": [], + "split": [], + "fold": [], + "target_group": [], + "time_used": [], +} + +LABEL_FILE_LINK = { + "HCP": "https://zenodo.org/records/10050233/files/HCP_half_brain.csv", + "GSP": "https://zenodo.org/records/10050234/files/gsp_half_brain.csv", +} + +GSDA_INIT_ARGS = ["lr", "max_iter", "l2_hparam", "lambda_", "optimizer", "max_iter"] +GSDA_FIT_ARGS = ["y", "groups", "target_idx"] + +GROUP_DICT = {0: "Male", 1: "Female"} + + +def run_experiment(cfg, lambda_): + atlas = cfg.DATASET.ATLAS + download = cfg.DATASET.DOWNLOAD + connection_type = cfg.DATASET.CONNECTION + data_dir = cfg.DATASET.ROOT + dataset = cfg.DATASET.DATASET.upper() + # lambda_list = cfg.SOLVER.LAMBDA_ + run_ = cfg.DATASET.RUN + random_state = cfg.SOLVER.SEED + test_size = cfg.DATASET.TEST_RATIO + + l2_hparam = cfg.SOLVER.L2_HPARAM + lr = cfg.SOLVER.LR + optimizer = cfg.SOLVER.OPTIMIZER + num_repeat = cfg.DATASET.NUM_REPEAT + max_iter = cfg.SOLVER.MAX_ITER + mix_group = cfg.DATASET.MIX_GROUP + rest1_only = cfg.DATASET.REST1_ONLY + + if mix_group: + # lambda_list = [0.0] + lambda_ = 0.0 + + label_file = "%s_%s_half_brain.csv" % (cfg.DATASET.DATASET, atlas) + label_fpath = os.path.join(data_dir, label_file) + if not os.path.exists(label_fpath): + if download: + os.makedirs(data_dir, exist_ok=True) + print("Downloading label file for %s dataset. \n" % dataset) + download_url_to_file(LABEL_FILE_LINK[dataset], label_fpath) + else: + raise ValueError("File %s does not exist" % label_file) + labels = read_tabular(label_fpath, index_col="ID") + group_label = labels["gender"].values + + # for lambda_ in lambda_list: + if mix_group: + out_folder = "lambda0_group_mix" + else: + out_folder = "lambda%s" % int(lambda_) + out_dir = os.path.join(cfg.OUTPUT.ROOT, out_folder) + if not os.path.exists(out_dir): + os.makedirs(out_dir) + + kwargs = { + "groups": group_label, + "lambda_": lambda_, + "l2_hparam": l2_hparam, + "lr": lr, + "mix_group": mix_group, + "out_dir": out_dir, + "num_repeat": num_repeat, + "test_size": test_size, + "random_state": random_state, + "rest1_only": rest1_only, + "optimizer": optimizer, + "max_iter": max_iter, + } + + if dataset == "HCP": + data = dict() + sessions = ["REST1", "REST2"] + for session in sessions: + data[session] = load_half_brain(data_dir, atlas, session, run_, connection_type, download=download) + kwargs = {**{"data": data}, **kwargs} + print("Training model with lambda = %s on %s dataset. \n" % (lambda_, dataset)) + if test_size == 0: + results = run_no_sub_hold_hcp(**kwargs) + elif 0 < test_size < 1: + results = run_sub_hold_hcp(**kwargs) + else: + raise ValueError("Invalid test_size %s" % test_size) + + elif dataset == "GSP": + data = load_half_brain( + data_dir, + atlas, + session=None, + run=run_, + connection_type=connection_type, + dataset=dataset, + download=download, + ) + kwargs = {**{"data": data}, **kwargs} + print("Training model with lambda = %s on %s dataset. \n" % (lambda_, dataset)) + if 0 < test_size < 1: + results = run_sub_hold_gsp(**kwargs) + elif test_size == 0: + results = run_sub_hold_gsp(**kwargs) + else: + raise ValueError("Invalid test_size %s." % test_size) + else: + raise ValueError("Invalid dataset %s." % dataset) + + res_df = pd.DataFrame.from_dict(results) + out_filename = "results_%s_L%s_test_size0%s_%s_%s" % ( + dataset, + int(lambda_), + str(int(test_size * 10)), + run_, + random_state, + ) + + if mix_group: + out_filename = out_filename + "_group_mix" + out_file = os.path.join(out_dir, "%s.csv" % out_filename) + print("Saving classification results to %s. \n" % out_file) + res_df.to_csv(out_file, index=False) + + # return res_df, out_file + + +def train_modal( + x_train, + init_kws, + fit_kws, + out_dir, + model_filename, +): + model_path = os.path.join(out_dir, "%s.pt" % model_filename) + + if os.path.exists(model_path): + model = torch.load(model_path) + else: + model = GSDA(**init_kws) + model.fit(x_train, **fit_kws) + torch.save(model, model_path) + print("Saving trained model to %s. \n" % model_path) + + return model + + +def save_loop_results( + model, + res_dict, + xy_test, + target_group, + train_session=None, + lambda_=0.0, + i_split=0, + train_fold=0, +): + for acc_key in xy_test: + test_x, test_y = xy_test[acc_key] + y_pred_ = model.predict(test_x) + acc_ = accuracy_score(test_y, y_pred_) + res_dict[acc_key].append(acc_) + + res_dict["pred_loss"].append(model.losses["pred"][-1]) + if "hsic" in model.losses: + res_dict["hsic_loss"].append(model.losses["hsic"][-1]) + else: + res_dict["hsic_loss"].append(model.losses["code"][-1]) + + # res['n_iter'].append(n_iter) + # n_iter += 1 + res_dict["target_group"].append(GROUP_DICT[target_group]) + if train_session is not None: + res_dict["train_session"].append(train_session) + else: + res_dict["train_session"].append(0) + res_dict["split"].append(i_split) + res_dict["fold"].append(train_fold) + + res_dict["lambda"].append(lambda_) + # res['test_size'].append(test_size) + res_dict["time_used"].append(model.losses["time"][-1]) + + return res_dict + + +def split_by_group(groups, train_subject, test_subject, target_group=0): + """Split the data by group. + + Args: + groups (np.ndarray): Group labels. + train_subject (np.ndarray): Training subjects indices. + test_subject (np.ndarray): Testing subjects indices. + target_group (int, optional): Target group. Defaults to 0. + """ + train_sub_tgt_idx = np.where(groups[train_subject] == target_group)[0] + train_sub_nt_idx = np.where(groups[train_subject] == 1 - target_group)[0] + test_sub_tgt_idx = np.where(groups[test_subject] == target_group)[0] + test_sub_nt_idx = np.where(groups[test_subject] == 1 - target_group)[0] + + return train_sub_tgt_idx, train_sub_nt_idx, test_sub_tgt_idx, test_sub_nt_idx + + +def run_no_sub_hold_hcp( + data, + groups, + lambda_, + l2_hparam, + mix_group, + out_dir, + num_repeat, + random_state, + **kwargs, +): + # main loop + res = { + **{ + "acc_tgt_train_session": [], + "acc_tgt_test_session": [], + "acc_nt_train_session": [], + "acc_nt_test_session": [], + }, + **copy.deepcopy(BASE_RESULT_DICT), + } + + for train_session, test_session in [("REST1", "REST2"), ("REST2", "REST1")]: + if train_session == "REST2" and kwargs["rest1_only"]: + continue + for i_split in range(num_repeat): + x_all = dict() + y_all = dict() + for session in [train_session, test_session]: + x, y, x1, y1 = pick_half(data[session], random_state=random_state * (i_split + 1)) + # x, y = _pick_half_subs(data[run_]) + + x_all[session] = [x, x1] + y_all[session] = [y, y1] + + for i_fold in [0, 1]: + x_train_fold = x_all[train_session][i_fold] + y_train_fold = y_all[train_session][i_fold] + + # scaler = StandardScaler() + # scaler.fit(x_train_fold) + for tgt_group in [0, 1]: + tgt_idx = np.where(groups == tgt_group)[0] + nt_idx = np.where(groups == 1 - tgt_group)[0] + + if mix_group: + if tgt_group == 1: + continue + tgt_group = "mix" + + xy_test = { + "acc_tgt_train_session": [ + x_all[train_session][1 - i_fold][tgt_idx], + y_all[train_session][1 - i_fold][tgt_idx], + ], + "acc_tgt_test_session": [ + np.concatenate([x_all[test_session][i][tgt_idx] for i in range(2)]), + np.concatenate([y_all[test_session][i][tgt_idx] for i in range(2)]), + ], + "acc_nt_train_session": [ + x_all[train_session][1 - i_fold][nt_idx], + y_all[train_session][1 - i_fold][nt_idx], + ], + "acc_nt_test_session": [ + np.concatenate([x_all[test_session][i][nt_idx] for i in range(2)]), + np.concatenate([y_all[test_session][i][nt_idx] for i in range(2)]), + ], + } + + fit_kws = { + "y": y_train_fold[tgt_idx], + "groups": groups, + "target_idx": tgt_idx, + } + model_filename = "HCP_L%s_test_size00_%s_%s_%s_group_%s_%s" % ( + int(lambda_), + train_session, + i_split, + i_fold, + tgt_group, + random_state, + ) + if mix_group: + model_filename = model_filename + "_group_mix" + fit_kws = { + "y": y_train_fold, + "groups": groups, + "target_idx": None, + } + + init_kws = {"lambda_": lambda_, "l2_hparam": l2_hparam} + for arg in GSDA_INIT_ARGS: + if arg in kwargs: + init_kws[arg] = kwargs[arg] + + model = train_modal( + x_train_fold, + init_kws, + fit_kws, + out_dir, + model_filename, + ) + res = save_loop_results( + model, + res, + xy_test, + tgt_group, + train_session, + lambda_, + i_split, + i_fold, + ) + + return res + + +def run_sub_hold_hcp( + data, + groups, + lambda_, + l2_hparam, + mix_group, + out_dir, + num_repeat, + test_size, + random_state, + **kwargs, +): + res = { + **{ + "acc_tgt_train_session": [], + "acc_tgt_test_session": [], + "acc_nt_train_session": [], + "acc_nt_test_session": [], + "acc_tgt_test_sub": [], + "acc_nt_test_sub": [], + }, + **copy.deepcopy(BASE_RESULT_DICT), + } + x_all = dict() + y_all = dict() + + for session in ["REST1", "REST2"]: + x, y, x1, y1 = pick_half(data[session], random_state=random_state) + + x_all[session] = [x, x1] + y_all[session] = [y, y1] + + for train_session, test_session in [("REST1", "REST2"), ("REST2", "REST1")]: + if train_session == "REST2" and kwargs["rest1_only"]: + continue + sss = StratifiedShuffleSplit(n_splits=num_repeat, test_size=test_size, random_state=random_state) + + for i_split, (train_sub, test_sub) in enumerate(sss.split(x, groups)): + for train_fold in [0, 1]: + x_train_fold = x_all[train_session][train_fold] + y_train_fold = y_all[train_session][train_fold] + + for tgt_group in [0, 1]: + _idx = split_by_group(groups, train_sub, test_sub, tgt_group) + ( + train_sub_tgt_idx, + train_sub_nt_idx, + test_sub_tgt_idx, + test_sub_nt_idx, + ) = _idx + + if mix_group: + if tgt_group == 1: + continue + tgt_group = "mix" + + xy_test = { + "acc_tgt_train_session": [ + x_all[train_session][1 - train_fold][train_sub][train_sub_tgt_idx], + y_all[train_session][1 - train_fold][train_sub][train_sub_tgt_idx], + ], + "acc_tgt_test_session": [ + np.concatenate([x_all[test_session][i][train_sub][train_sub_tgt_idx] for i in range(2)]), + np.concatenate([y_all[test_session][i][train_sub][train_sub_tgt_idx] for i in range(2)]), + ], + "acc_nt_train_session": [ + x_all[train_session][1 - train_fold][train_sub][train_sub_nt_idx], + y_all[train_session][1 - train_fold][train_sub][train_sub_nt_idx], + ], + "acc_nt_test_session": [ + np.concatenate([x_all[test_session][i][train_sub][train_sub_nt_idx] for i in range(2)]), + np.concatenate([y_all[test_session][i][train_sub][train_sub_nt_idx] for i in range(2)]), + ], + "acc_tgt_test_sub": [[], []], + "acc_nt_test_sub": [[], []], + } + + for _s in [train_session, test_session]: + for _fold in (0, 1): + xy_test["acc_tgt_test_sub"][0].append(x_all[_s][_fold][test_sub][test_sub_tgt_idx]) + xy_test["acc_tgt_test_sub"][1].append(y_all[_s][_fold][test_sub][test_sub_tgt_idx]) + + xy_test["acc_nt_test_sub"][0].append(x_all[_s][_fold][test_sub][test_sub_nt_idx]) + xy_test["acc_nt_test_sub"][1].append(y_all[_s][_fold][test_sub][test_sub_nt_idx]) + + for _key in ["acc_tgt_test_sub", "acc_nt_test_sub"]: + for _idx in range(2): + xy_test[_key][_idx] = np.concatenate(xy_test[_key][_idx]) + + fit_kws = { + "y": y_train_fold[train_sub][train_sub_tgt_idx], + "groups": groups[train_sub], + "target_idx": train_sub_tgt_idx, + } + model_filename = "HCP_L%s_test_size0%s_%s_%s_%s_group_%s_%s" % ( + int(lambda_), + str(int(test_size * 10)), + train_session, + i_split, + train_fold, + tgt_group, + random_state, + ) + x_train = x_train_fold[train_sub] + + if mix_group: + model_filename = model_filename + "_group_mix" + fit_kws = { + "y": y_train_fold[train_sub], + "groups": groups[train_sub], + "target_idx": None, + } + + init_kws = {"lambda_": lambda_, "l2_hparam": l2_hparam} + for arg in GSDA_INIT_ARGS: + if arg in kwargs: + init_kws[arg] = kwargs[arg] + + model = train_modal( + x_train, + init_kws, + fit_kws, + out_dir, + model_filename, + ) + res = save_loop_results( + model, + res, + xy_test, + tgt_group, + train_session, + lambda_, + i_split, + train_fold, + ) + + return res + + +def run_sub_hold_gsp( + data, + groups, + lambda_, + l2_hparam, + mix_group, + out_dir, + num_repeat, + test_size, + random_state, + **kwargs, +): + res = { + **{"acc_tgt": [], "acc_nt": [], "acc_tgt_test_sub": [], "acc_nt_test_sub": []}, + **copy.deepcopy(BASE_RESULT_DICT), + } + + x, y, x1, y1 = pick_half(data, random_state=random_state) + + x_all = [x, x1] + y_all = [y, y1] + + sss = StratifiedShuffleSplit(n_splits=num_repeat, test_size=test_size, random_state=random_state) + for i_split, (train_sub, test_sub) in enumerate(sss.split(x, groups)): + for train_fold in [0, 1]: + x_train_fold = x_all[train_fold] + y_train_fold = y_all[train_fold] + x_test_fold = x_all[1 - train_fold] + y_test_fold = y_all[1 - train_fold] + + for tgt_group in [0, 1]: + _idx = split_by_group(groups, train_sub, test_sub, tgt_group) + ( + train_sub_tgt_idx, + train_sub_nt_idx, + test_sub_tgt_idx, + test_sub_nt_idx, + ) = _idx + + if mix_group: + if tgt_group == 1: + continue + tgt_group = "mix" + + x_train_fold_test_tgt = x_test_fold[train_sub][train_sub_tgt_idx] + y_train_fold_test_tgt = y_test_fold[train_sub][train_sub_tgt_idx] + x_train_fold_test_nt = x_test_fold[train_sub][train_sub_nt_idx] + y_train_fold_test_nt = y_test_fold[train_sub][train_sub_nt_idx] + + fit_kws = { + "y": y_train_fold[train_sub_tgt_idx], + "groups": groups, + "target_idx": train_sub_tgt_idx, + } + model_filename = "GSP_L%s_test_size0%s_%s_%s_group_%s_%s" % ( + int(lambda_), + str(int(test_size * 10)), + i_split, + train_fold, + tgt_group, + random_state, + ) + if 0 < test_size < 1: + x_train = x_train_fold[train_sub] + xy_test = { + "acc_tgt": [x_train_fold_test_tgt, y_train_fold_test_tgt], + "acc_nt": [x_train_fold_test_nt, y_train_fold_test_nt], + "acc_tgt_test_sub": [ + np.concatenate( + ( + x_train_fold[test_sub][test_sub_tgt_idx], + x_test_fold[test_sub][test_sub_tgt_idx], + ) + ), + np.concatenate( + ( + y_train_fold[test_sub][test_sub_tgt_idx], + y_test_fold[test_sub][test_sub_tgt_idx], + ) + ), + ], + "acc_nt_test_sub": [ + np.concatenate( + ( + x_train_fold[test_sub][test_sub_nt_idx], + x_test_fold[test_sub][test_sub_nt_idx], + ) + ), + np.concatenate( + ( + y_train_fold[test_sub][test_sub_nt_idx], + y_test_fold[test_sub][test_sub_nt_idx], + ) + ), + ], + } + model_filename = model_filename + "_test_sub_0%s" % str(int(test_size * 10)) + fit_kws = { + "y": y_train_fold[train_sub][train_sub_tgt_idx], + "groups": groups[train_sub], + "target_idx": train_sub_tgt_idx, + } + else: + x_train = x_train_fold + xy_test = { + "acc_tgt": [x_train_fold_test_tgt, y_train_fold_test_tgt], + "acc_nt": [x_train_fold_test_nt, y_train_fold_test_nt], + } + + if mix_group: + model_filename = model_filename + "_group_mix" + if 0 < test_size < 1: + fit_kws = { + "y": y_train_fold[train_sub], + "group": groups[train_sub], + "target_idx": None, + } + else: + fit_kws = { + "y": y_train_fold, + "group": groups, + "target_idx": None, + } + init_kws = {"lambda_": lambda_, "l2_hparam": l2_hparam} + for arg in GSDA_INIT_ARGS: + if arg in kwargs: + init_kws[arg] = kwargs[arg] + + model = train_modal( + x_train, + init_kws, + fit_kws, + out_dir, + model_filename, + ) + res = save_loop_results( + model, + res, + xy_test, + tgt_group, + lambda_=lambda_, + i_split=i_split, + train_fold=train_fold, + ) + + return res + + +def run_no_sub_hold_gsp( + data, + groups, + lambda_, + l2_hparam, + mix_group, + out_dir, + num_repeat, + random_state, + **kwargs, +): + res = { + **{"acc_tgt": [], "acc_nt": []}, + **copy.deepcopy(BASE_RESULT_DICT), + } + for i_split in range(num_repeat): + x, y, x1, y1 = pick_half(data, random_state=random_state * (i_split + 1)) + x_all = [x, x1] + y_all = [y, y1] + + for i_fold in [0, 1]: + x_train = x_all[i_fold] + y_train = y_all[i_fold] + + # scaler = StandardScaler() + # scaler.fit(x_train) + for tgt_group in [0, 1]: + tgt_idx = np.where(groups == tgt_group)[0] + nt_idx = np.where(groups == 1 - tgt_group)[0] + + if mix_group: + if tgt_group == 1: + continue + tgt_group = "mix" + + xy_test = { + "acc_tgt": [x_all[1 - i_fold][tgt_idx], x_all[1 - i_fold][tgt_idx]], + "acc_nt": [x_all[1 - i_fold][nt_idx], x_all[1 - i_fold][nt_idx]], + } + + fit_kws = { + "y": y_train[tgt_idx], + "groups": groups, + "target_idx": tgt_idx, + } + model_filename = "GSP_L%s_test_size00_%s_%s_group_%s_%s" % ( + int(lambda_), + i_split, + i_fold, + tgt_group, + random_state, + ) + + if mix_group: + model_filename = model_filename + "_group_mix" + fit_kws = { + "y": y_train, + "groups": groups, + "target_idx": None, + } + init_kws = {"lambda_": lambda_, "l2_hparam": l2_hparam} + for arg in GSDA_INIT_ARGS: + if arg in kwargs: + init_kws[arg] = kwargs[arg] + + model = train_modal( + x_train, + init_kws, + fit_kws, + out_dir, + model_filename, + ) + res = save_loop_results( + model, + res, + xy_test, + tgt_group, + lambda_=lambda_, + i_split=i_split, + train_fold=i_fold, + ) + + return res diff --git a/examples/brain_lateralization/utils/io_.py b/examples/brain_lateralization/utils/io_.py new file mode 100644 index 0000000..47eb47d --- /dev/null +++ b/examples/brain_lateralization/utils/io_.py @@ -0,0 +1,527 @@ +import os + +import h5py +import matplotlib.pyplot as plt +import nibabel as nib +import numpy as np +import pandas as pd +import seaborn as sns +import torch +from scipy.io import loadmat, savemat +from scipy.stats import pearsonr +from sklearn.model_selection import train_test_split +from sklearn.preprocessing import label_binarize +from torch.hub import download_url_to_file + +HCP_LINK = { + "REST1": "https://zenodo.org/records/10050233/files/HCP_BNA_intra_half_brain_REST1_Fisherz.hdf5", + "REST2": "https://zenodo.org/records/10050233/files/HCP_BNA_intra_half_brain_REST2_Fisherz.hdf5", +} +GSP_LINK = "https://zenodo.org/records/10050234/files/gsp_BNA_intra_half_brain_Fisherz.mat" + + +def load_txt(fpaths, connection_type="intra"): + """Load data from a list of txt files。 + + Parameters + ---------- + fpaths (list): A list of .txt file paths. + + Returns + ------- + Two ndarrays + """ + left_ = [] + right_ = [] + for fpath in fpaths: + data_matrix = np.genfromtxt(fpath) + left_vec, right_vec = split_functional_brain(data_matrix, connection_type=connection_type) + left_.append(left_vec.reshape((1, -1))) + right_.append(right_vec.reshape((1, -1))) + left_ = np.concatenate(left_, axis=0) + right_ = np.concatenate(right_, axis=0) + + return left_, right_ + + +def load_hdf5(fpath): + """ + + Parameters + ---------- + fpath + + Returns + ------- + + """ + f = h5py.File(fpath, "r") + data = {"Left": f["Left"][()], "Right": f["Right"][()]} + return data + + +def read_tabular(fname, **kwargs): + """Read a table from a .xlsx or .csv file + + Parameters + ---------- + fname (string): + + Returns + ------- + + """ + file_format = fname.split(".")[-1] + if file_format == "xlsx": + df = pd.read_excel(fname, engine="openpyxl", **kwargs) + elif file_format == "csv": + df = pd.read_csv(fname, **kwargs) + else: + raise ValueError("Unsupported file type %s" % file_format) + + return df + + +def get_fpaths(fdir, idx_list, file_format="txt"): + """Check the existence of files named by subject indices under a given directory + + Parameters + ---------- + fdir (string): file directory for data + idx_list (list): a list of subject indices + file_format (string): file format, defaults to txt + + Returns + ------- + Dataframe + """ + fpaths = dict() + for idx in idx_list: + fname = "%s.%s" % (idx, file_format) + fpath = os.path.join(fdir, fname) + if os.path.exists(fpath): + fpaths[idx] = fpath + + fpath_df = pd.DataFrame(data={"File path": fpaths.values()}, index=list(fpaths.keys())) + + return fpath_df + + +def split_functional_brain(matrix, connection_type="intra"): + """ + + Parameters + ---------- + matrix + connection_type + + Returns + ------- + + """ + n_ = matrix.shape[1] / 2 + if connection_type == "intra": + left = matrix[0::2, 0::2] # even indices of rows, left brain + right = matrix[1::2, 1::2] # odd indices of rows, right brain + + idx = np.triu_indices(n_, k=1) + n_feat = int(n_ * (n_ - 1) / 2) + left_vec = np.zeros((1, n_feat)) + right_vec = np.zeros((1, n_feat)) + left_vec[0, :] = left[idx] + right_vec[0, :] = right[idx] + + elif connection_type == "inter": + left = matrix[0::2, 1::2] + right = matrix[1::2, 0::2] + left_vec = left.reshape((1, -1)) + right_vec = right.reshape((1, -1)) + + else: + raise ValueError("Invalid connection type %s" % connection_type) + + return left_vec, right_vec + + +def save_half_brain(out_dir, out_fname, data_left, data_right): + """Save two hemisphere nadarrays into a hdf5 file + + Parameters + ---------- + out_dir (string): + out_fname (string): + data_left (ndarray): + data_right (ndarray): + + Returns + ------- + + """ + f = h5py.File(os.path.join(out_dir, out_fname), "w") + f.create_dataset("Left", data=data_left) + f.create_dataset("Right", data=data_right) + f.close() + + +def save_half_brain_mat(out_dir, out_fname, data_left, data_right): + """Save two hemisphere nadarrays into a hdf5 file + + Parameters + ---------- + out_dir (string): + out_fname (string): + data_left (ndarray): + data_right (ndarray): + + Returns + ------- + + """ + savemat(os.path.join(out_dir, out_fname), {"Left": data_left, "Right": data_right}) + + +def load_half_brain( + data_dir, + atlas, + session=None, + run=None, + connection_type="intra", + data_type="functional", + dataset="HCP", + download=True, +): + """ + + Parameters + ---------- + data_dir + atlas + session + run + connection_type + data_type + dataset + download: bool, optional, download the data if the data files not exist (default=True) + + Returns + ------- + + """ + if dataset == "HCP": + if data_type == "functional": + fname = "HCP_%s_%s_half_brain_%s_%s.hdf5" % ( + atlas, + connection_type, + session, + run, + ) + fpath = os.path.join(data_dir, fname) + if not os.path.exists(fpath): + if download: + os.makedirs(data_dir, exist_ok=True) + print("Downloading %s session %s data, it may take 1-2 mins" % (dataset, session)) + download_url_to_file( + HCP_LINK[session], + fpath, + ) + else: + raise ValueError("File %s does not exist" % fpath) + data = load_hdf5(fpath) + elif data_type == "structural": + data = {"Left": [], "Right": []} + fname = "%s_Volume.mat" % atlas + data_in = loadmat(os.path.join(data_dir, fname))["%s_Volume" % atlas][0][0][0] + data["Left"] = data_in[:, 0::2] + data["Right"] = data_in[:, 1::2] + else: + raise ValueError("Invalid data type %s" % data_type) + elif dataset in ["ABIDE", "ukb", "GSP"]: + fpath = os.path.join( + data_dir, + "%s_%s_%s_half_brain_%s.mat" % (dataset, atlas, connection_type, run), + ) + if not os.path.exists(fpath): + if download: + if dataset == "GSP": + os.makedirs(data_dir, exist_ok=True) + print("Downloading %s data, it may take 1-2 mins." % dataset) + download_url_to_file(GSP_LINK, fpath) + else: + raise ValueError("File %s does not exist" % fpath) + else: + raise ValueError("File %s does not exist" % fpath) + data_file = loadmat(fpath) + data = {"Left": data_file["Left"], "Right": data_file["Right"]} + else: + raise ValueError("Invalid dataset %s" % dataset) + + return data + + +def fetch_weights(base_dir, group, lambda_, dataset, sessions, test_size="00", seed_=2023): + """ + + Args: + base_dir: + group: int, 0 or 1 + lambda_: str, "0_group_mix" or "0", "1", "2", "5", "8", "10" + dataset: str, "HCP" or "GSP" + sessions: list of session names, ["REST1_", "REST2_"] for HCP and ["_"] for GSP + test_size: str, "00" or "02", optional, default="00" + seed_: int, optional, default=2023 + + Returns: + a matrix of weights, shape (n_models, n_features) + """ + + sub_dir = os.path.join(base_dir, "lambda%s" % lambda_) + if not os.path.exists(sub_dir): + return None + if lambda_ == "0_group_mix": + lambda_ = 0 + group = "mix" + weight = [] + num_repeat = 5 + halfs = [0, 1] + for session_i in sessions: + for half_i in halfs: + for i_split in range(num_repeat): + for seed in range(52): + model_file_name = "%s_L%s_test_size%s_%s%s_%s_group_%s_%s.pt" % ( + dataset, + lambda_, + test_size, + session_i, + i_split, + half_i, + group, + seed_ - seed, + ) + if os.path.exists(os.path.join(sub_dir, model_file_name)): + weight.append(get_coef(model_file_name, sub_dir).reshape((1, -1))) + + return np.concatenate(weight, axis=0) + + +def get_coef(file_name, file_dir): + file_path = os.path.join(file_dir, file_name) + model = torch.load(file_path) + return model.theta + + +def save_results(res_dict, out_filename, output_dir, mix_group=False): + """save a dictionary to a csv file + + Args: + res_dict (_type_): _description_ + out_filename (str): _description_ + output_dir (str): _description_ + mix_group (bool, optional): _description_. Defaults to False. + """ + res_df = pd.DataFrame.from_dict(res_dict) + + if mix_group: + out_filename = out_filename + "_mix_group" + out_file = os.path.join(output_dir, "%s.csv" % out_filename) + res_df.to_csv(out_file, index=False) + + +# read: nib.load() +# write: nib.save() +# ********************************** file split ************************************* # +def file_split(file_path): + filepath, tempfilename = os.path.split(file_path) + filename, extension = os.path.splitext(tempfilename) + return filepath, filename, extension + + +# ********************************** read h5 or pickle file as a df ************************************* # + + +def read_file(file): + _, _, ext = file_split(file) + if ext == ".h5": + df = pd.read_hdf(file) + elif ext == ".pkl": + df = pd.read_pickle(file) + else: + raise ValueError("Unsupported file type %s, supported type: xxx.h5 or xxx.pkl..." % ext) + + return df + + +def select_data_subj(df, subj_id, df_column_name): + data = df[df_column_name][df["subject"] == subj_id] + return data # series datatype with dim of (1,),for the ndarray contents, data = data[0] + + +def pick_half(data, random_state=144): + x = np.zeros(data["Left"].shape) + left_idx, right_idx = train_test_split(range(x.shape[0]), test_size=0.5, random_state=random_state) + x[left_idx] = data["Left"][left_idx] + x[right_idx] = data["Right"][right_idx] + + n_sub = x.shape[0] + y = np.zeros(n_sub) + y[left_idx] = 1 + y[right_idx] = -1 + + x1 = np.zeros(data["Left"].shape) + x1[left_idx] = data["Right"][left_idx] + x1[right_idx] = data["Left"][right_idx] + + y1 = np.zeros(n_sub) + y1[left_idx] = -1 + y1[right_idx] = 1 + + y = label_binarize(y, classes=[-1, 1]).reshape(-1) + y1 = label_binarize(y1, classes=[-1, 1]).reshape(-1) + + return x, y, x1, y1 + + +def _pick_half_subs(data, random_state=144): + n_ = data["Left"].shape[0] + train_idx, hold_idx = train_test_split(range(n_), test_size=0.5, random_state=random_state) + x = np.concatenate([data["Left"][train_idx], data["Right"][train_idx]], axis=0) + y = np.ones(n_) + y[int(n_ / 2) :] = -1 + + return x, y + + +def select_data_multi_subj(df, subj_ids, df_column_name): + subj_df = df["subject"].to_numpy() + sub, index1, index2 = np.intersect1d(subj_df, subj_ids, assume_unique=False, return_indices=True) + data = df.loc[index1, df_column_name] # loc,select column name. iloc,select data as index + return data.to_numpy() # convert pd series to numpy ndarray + + +# ********************************** read MRI files ********************************* # +def read_surf(surface): + surf = nib.load(surface) + arr = surf.darrays + coord = arr[0].data # coordinate of vertices + vert = arr[1].data # index of vertices + return coord, vert + + +def read_nii(T1w): + T1 = nib.load(T1w) + data = T1.get_fdata() # 3-D ‘ray + header = T1.header # header + return data, header + + +def read_shape_gii(shape_gii): + gii = nib.load(shape_gii) + data = gii.darrays[0].data + return gii, data + + +# ********************************** Creat a .shape.gii file ********************************* # + + +def creat_shape_gii(shape_gii_base, array_write, savepath): + gii, _ = read_shape_gii(shape_gii_base) + gii.darrays[0].data = array_write + nib.save(gii, savepath) + + +# Creat a .shape gii file with atlas, e.g.,Kong400,Schaefer400 +def creat_fs_lr32k_atlas(shape_gii_base, array_write, atlas_file, hemi, savepath): + # cdata = np.zeros((32492,)) + f_atlas = nib.load(atlas_file) + f_data = f_atlas.get_fdata() # 0 mean medial wall, label starts from 1 + f_data = f_data[0, :] + + if hemi == "L": + cdata = f_data[0 : int(len(f_data) / 2)].copy() # L + print("cdata.shape: %d" % (len(cdata))) + for i, data in enumerate(array_write): + cdata[np.where(cdata == i + 1)] = data + else: + cdata = f_data[int(len(f_data) / 2) : int(len(f_data))].copy() # R + for i, data in enumerate(array_write): + cdata[np.where(cdata == i + 1 + int(np.max(f_data) / 2))] = data + + creat_shape_gii(shape_gii_base, cdata, savepath) + + +# *********************************************** Corr with pvalue ************************************************ # +def Corr(x, y): + r, p = pearsonr(x, y) + rr = round(r, 3) + pp = round(p, 3) + return rr, pp + + +# ************************************************ Join plot ***************************************************** # +# x,y should not be the object data type, if so, convert it to np.double data type. +def jointplot_fitlinear(x, y): + plt.figure(figsize=(20, 16)) + sns.jointplot(x=x, y=y, kind="reg", scatter_kws={"s": 8}) + + +# sns.regplot(x='LL_LPAC_Math_Corr',y='ReadEng_AgeAdj',data=df, color='red', scatter_kws={"s": 8}),不显示分布信息 + + +# *************************************** Extract top N in a array ************************************ # +# Extract top N element in an array and set others to 0. +def topN_array(arr, topN): + idx_all = set(np.arange(len(arr))) + idx_topN = arr.argsort()[::-1][:topN] + idx_res = idx_all - set(idx_topN) + # keep topN, set others to zero + idx_res = np.array(list(idx_res)) + arr_topN = arr.copy() + arr_topN[idx_res] = 0 + return arr_topN + + +# def main_creat_shape_gii(): +# +# tasks = ['LANGUAGE'] +# subjs = sio.loadmat(project_path + '/Subject_taskT.mat') +# subjects = subjs['Subject_taskT'] +# subjects = np.reshape(subjects,(len(subjects),)) +# for creat shape gii files +# shape_gii_base_L = project_path + '/Base_shape_gii/100206.L.thickness.32k_fs_LR.shape.gii' +# for creat shape gii files +# shape_gii_base_R = project_path + '/Base_shape_gii/100206.R.thickness.32k_fs_LR.shape.gii' +# +# for subj_id in subjects: +# df_column_name = 'task_t' +# for task in tasks: +# files = sorted(glob.glob(project_path + '/DataSort/Sort_Sub/' + task + '/*.h5')) +# +# for file in files: +# +# Dirs = os.path.split(os.path.splitext(file)[0]) +# kName = Dirs[1] +# df = rvp.read_orig_data(file, kName) +# ResultantFolder = (project_path + '/Result/Ridge_Vert_SubBased/Statistics/task_t_gii/' +# + str(subj_id) +'/' + task) +# if not os.path.exists(ResultantFolder): +# os.makedirs(ResultantFolder) +# task_t = rvp.select_data_subj(df, subj_id, 'task_t').to_numpy()[0] +# task_t = np.reshape(task_t, (len(task_t),)) +# vert_index = rvp.select_data_subj(df, subj_id, 'vert_index').to_numpy()[0] +# cdata = np.zeros((32492,)) +# vert_index = np.reshape(vert_index, (len(vert_index),)) +# cdata[vert_index - 1] = task_t # compared to matlab, the vert index should minus 1 +# +# save_name = ResultantFolder + '/' + str(subj_id) + '_' + kName + '.shape.gii' +# if 'LW' in file: +# print('LEFT SEED') +# shape_gii_base= shape_gii_base_L +# elif 'RW' in file: +# print('RIGHT SEED') +# shape_gii_base= shape_gii_base_R +# else: +# print('error!') +# +# creat_shape_gii(shape_gii_base, cdata, save_name) +# print(kName) +# +# print(str(subj_id)) diff --git a/examples/brain_lateralization/utils/plot.py b/examples/brain_lateralization/utils/plot.py new file mode 100644 index 0000000..c980aa0 --- /dev/null +++ b/examples/brain_lateralization/utils/plot.py @@ -0,0 +1,201 @@ +import matplotlib.pylab as plt +import numpy as np +import pandas as pd +import seaborn as sns +from scipy.io import loadmat +from utils.io_ import fetch_weights + + +def savefig(fig, outfile, outfig_format): + if isinstance(outfig_format, list): + for fmt in outfig_format: + fig.savefig("%s.%s" % (outfile, fmt), format=fmt, bbox_inches="tight") + else: + fig.savefig( + "%s.%s" % (outfile, outfig_format), + format=outfig_format, + bbox_inches="tight", + ) + + +# def plot_gsi( +# data, +# x="lambda", +# y="GSI", +# hue="train_gender", +# fontsize=14, +# outfile=None, +# outfig_format=None, +# ): +# fig = plt.figure() +# plt.rcParams.update({"font.size": fontsize}) +# sns.boxplot(data=data, x=x, y=y, hue=hue, showmeans=True) +# plt.legend(title="Target group") +# plt.rcParams["text.usetex"] = True +# plt.xlabel(r"$\lambda$") +# plt.rcParams["text.usetex"] = False +# plt.ylabel("Group Specificity Index (GSI)") +# if outfile is not None: +# savefig(fig, outfile, outfig_format) +# plt.show() +# +# +# def plot_accuracy( +# data, +# x="Lambda", +# y="Accuracy", +# col="Target group", +# hue="Test set", +# style="Test set", +# kind="line", +# height=4, +# fontsize=14, +# outfile=None, +# outfig_format=None, +# ): +# fig = plt.figure() +# plt.rcParams.update({"font.size": fontsize}) +# g = sns.relplot( +# data=data, +# x=x, +# y=y, +# col=col, +# hue=hue, +# style=style, +# kind=kind, +# errorbar=("sd", 1), +# height=height, +# ) +# ( +# g.map(plt.axhline, y=0.9, color=".7", dashes=(2, 1), zorder=0) +# .set_axis_labels(r"$\lambda$", "Test Accuracy") +# .set_titles("Target group: {col_name}") +# .tight_layout(w_pad=0) +# ) +# if outfile is not None: +# savefig(fig, outfile, outfig_format) +# +# plt.show() + + +def load_weight_plot_corr(dataset, base_dir, sessions, seed_start, fontsize=14): + control_weights = fetch_weights(base_dir, "mix", "0_group_mix", dataset, sessions=sessions, seed_=seed_start) + n_control_weights = control_weights.shape[0] + + lambdas = [0.0, 1.0, 2.0, 5.0, 8.0, 10.0] + corrs = {"mean": [], "sd": []} + + for lambda_ in lambdas: + for group in [0, 1]: + weights = fetch_weights( + base_dir, + group, + int(lambda_), + dataset, + sessions=sessions, + seed_=seed_start, + ) + # n_weights = weights.shape[0] + corr_matrix = np.corrcoef(control_weights, weights)[n_control_weights:, :n_control_weights] + corrs["mean"].append(np.mean(corr_matrix)) + corrs["sd"].append(np.std(corr_matrix)) + + plt.rcParams.update({"font.size": fontsize}) + fig, ax = plt.subplots(2, 1, sharex=True) + fig.set_size_inches(4, 8.5) + # .plot(x, y1) + # .plot(x, y2) + + # plt.figure(figsize=(4, 4)) + ax[0].plot(lambdas, corrs["mean"][::2], "-", c="steelblue", label="Male-specific") + ax[0].fill_between( + lambdas, + np.asarray(corrs["mean"][::2]) - np.asarray(corrs["sd"][::2]), + np.asarray(corrs["mean"][::2]) + np.asarray(corrs["sd"][::2]), + color="steelblue", + alpha=0.2, + ) + ax[0].set_ylabel("Correlation") + ax[0].legend(loc="upper right") + + ax[1].plot(lambdas, corrs["mean"][1::2], "--", color="firebrick", label="Female-specific") + ax[1].fill_between( + lambdas, + np.asarray(corrs["mean"][1::2]) - np.asarray(corrs["sd"][1::2]), + np.asarray(corrs["mean"][1::2]) + np.asarray(corrs["sd"][1::2]), + color="firebrick", + alpha=0.2, + ) + + ax[1].set_ylabel("Correlation") + ax[1].legend(loc="upper right") + plt.rcParams["text.usetex"] = True + plt.xlabel(r"$\lambda$", fontsize=fontsize) + plt.rcParams["text.usetex"] = False + + plt.savefig("figures/%s_corr.svg" % dataset, format="svg", bbox_inches="tight") + # plt.savefig('figures/%s_corr.pdf' % dataset, format='pdf', bbox_inches='tight') + # plt.savefig('figures/%s_corr.png' % dataset, format='png', bbox_inches='tight') + plt.show() + + +def load_coef_plot_corr(dataset, group_label, fontsize=14): + weights_dirs = dict() + lambdas = [0.0, 1.0, 2.0, 5.0, 8.0, 10.0] + group_dict = {0: "Male", 1: "Female"} + for lambda_ in lambdas: + weights_dirs[r"$\lambda=%s$" % (int(lambda_))] = "%s/%s_L%sG%s.mat" % ( + dataset, + dataset, + int(lambda_), + group_label, + ) + + weights = dict() + + for key in weights_dirs: + # print(key) + weights[key] = loadmat(weights_dirs[key])["mean"][0][1:] + + weight_df = pd.DataFrame(weights) + + corr = weight_df.corr() + + # Generate a mask for the upper triangle + mask = np.triu(np.ones_like(corr, dtype=bool), 1) + + plt.rcParams.update({"font.size": fontsize}) + # Set up the matplotlib figure + f, ax = plt.subplots(figsize=(11, 9)) + + # Generate a custom diverging colormap + cmap = sns.diverging_palette(230, 20, as_cmap=True) + plt.rcParams["text.usetex"] = True + # Draw the heatmap with the mask and correct aspect ratio + sns.heatmap( + corr, + mask=mask, + cmap=cmap, + vmin=0, + vmax=1, + center=0.5, + annot=True, + annot_kws={"fontsize": "xx-large"}, + square=True, + linewidths=0.5, + cbar_kws={"shrink": 0.5}, + ) + # cbar_kws={"shrink": .5, "use_gridspec": False, "location": "top"}) #, annot_kws={"rotation": 45}) + plt.rcParams["text.usetex"] = False + ax.set_xticklabels(weight_df.columns.to_list(), fontsize=fontsize + 2) + ax.set_yticklabels(weight_df.columns.to_list(), rotation=45, ha="right", fontsize=fontsize + 2) + # plt.savefig('corr.pdf', format='pdf', bbox_inches='tight') + # plt.savefig('corr.png', format='png', bbox_inches='tight') + plt.savefig( + "figures/corr_annot_%s_%s.svg" % (dataset, group_dict[group_label]), + format="svg", + bbox_inches="tight", + ) + # plt.savefig('figures/corr_annot_%s_%s.pdf' % (dataset, group_dict[sex_label]), format='pdf', bbox_inches='tight') + # plt.savefig('figures/corr_annot_%s_%s.png' % (dataset, group_dict[sex_label]), format='png', bbox_inches='tight') + plt.show() From 55127b7d4b9da501e381a0ef4d1da4099f0bb834 Mon Sep 17 00:00:00 2001 From: Shuo Zhou Date: Sat, 6 Jun 2026 00:02:52 +0100 Subject: [PATCH 3/6] update gsda exmaple --- .../brain_lateralization/utils/experiment.py | 183 +++++++----------- 1 file changed, 71 insertions(+), 112 deletions(-) diff --git a/examples/brain_lateralization/utils/experiment.py b/examples/brain_lateralization/utils/experiment.py index ab15c6b..babff69 100644 --- a/examples/brain_lateralization/utils/experiment.py +++ b/examples/brain_lateralization/utils/experiment.py @@ -31,10 +31,33 @@ "GSP": "https://zenodo.org/records/10050234/files/gsp_half_brain.csv", } -GSDA_INIT_ARGS = ["lr", "max_iter", "l2_hparam", "lambda_", "optimizer", "max_iter"] -GSDA_FIT_ARGS = ["y", "groups", "target_idx"] +GSDA_INIT_ARGS = ["lr", "max_iter", "l2_hparam", "lambda_", "optimizer"] -GROUP_DICT = {0: "Male", 1: "Female"} +GROUP_DICT = {0: "Male", 1: "Female", "mix": "Mix"} + + +def get_gsda_init_kws(lambda_, l2_hparam, extra_kws): + init_kws = {"lambda_": lambda_, "l2_hparam": l2_hparam} + for arg in GSDA_INIT_ARGS: + if arg in extra_kws: + init_kws[arg] = extra_kws[arg] + + return init_kws + + +def get_target_groups(mix_group): + if mix_group: + return ["mix"] + + return [0, 1] + + +def get_fit_kws(y, groups, target_idx=None): + return { + "y": y, + "groups": groups, + "target_idx": target_idx, + } def run_experiment(cfg, lambda_): @@ -125,7 +148,7 @@ def run_experiment(cfg, lambda_): if 0 < test_size < 1: results = run_sub_hold_gsp(**kwargs) elif test_size == 0: - results = run_sub_hold_gsp(**kwargs) + results = run_no_sub_hold_gsp(**kwargs) else: raise ValueError("Invalid test_size %s." % test_size) else: @@ -149,7 +172,7 @@ def run_experiment(cfg, lambda_): # return res_df, out_file -def train_modal( +def train_model( x_train, init_kws, fit_kws, @@ -266,14 +289,10 @@ def run_no_sub_hold_hcp( # scaler = StandardScaler() # scaler.fit(x_train_fold) - for tgt_group in [0, 1]: - tgt_idx = np.where(groups == tgt_group)[0] - nt_idx = np.where(groups == 1 - tgt_group)[0] - - if mix_group: - if tgt_group == 1: - continue - tgt_group = "mix" + for tgt_group in get_target_groups(mix_group): + group_idx = 0 if mix_group else tgt_group + tgt_idx = np.where(groups == group_idx)[0] + nt_idx = np.where(groups == 1 - group_idx)[0] xy_test = { "acc_tgt_train_session": [ @@ -294,11 +313,7 @@ def run_no_sub_hold_hcp( ], } - fit_kws = { - "y": y_train_fold[tgt_idx], - "groups": groups, - "target_idx": tgt_idx, - } + fit_kws = get_fit_kws(y_train_fold[tgt_idx], groups, tgt_idx) model_filename = "HCP_L%s_test_size00_%s_%s_%s_group_%s_%s" % ( int(lambda_), train_session, @@ -309,18 +324,11 @@ def run_no_sub_hold_hcp( ) if mix_group: model_filename = model_filename + "_group_mix" - fit_kws = { - "y": y_train_fold, - "groups": groups, - "target_idx": None, - } - - init_kws = {"lambda_": lambda_, "l2_hparam": l2_hparam} - for arg in GSDA_INIT_ARGS: - if arg in kwargs: - init_kws[arg] = kwargs[arg] - - model = train_modal( + fit_kws = get_fit_kws(y_train_fold, groups) + + init_kws = get_gsda_init_kws(lambda_, l2_hparam, kwargs) + + model = train_model( x_train_fold, init_kws, fit_kws, @@ -383,8 +391,9 @@ def run_sub_hold_hcp( x_train_fold = x_all[train_session][train_fold] y_train_fold = y_all[train_session][train_fold] - for tgt_group in [0, 1]: - _idx = split_by_group(groups, train_sub, test_sub, tgt_group) + for tgt_group in get_target_groups(mix_group): + group_idx = 0 if mix_group else tgt_group + _idx = split_by_group(groups, train_sub, test_sub, group_idx) ( train_sub_tgt_idx, train_sub_nt_idx, @@ -392,11 +401,6 @@ def run_sub_hold_hcp( test_sub_nt_idx, ) = _idx - if mix_group: - if tgt_group == 1: - continue - tgt_group = "mix" - xy_test = { "acc_tgt_train_session": [ x_all[train_session][1 - train_fold][train_sub][train_sub_tgt_idx], @@ -430,11 +434,9 @@ def run_sub_hold_hcp( for _idx in range(2): xy_test[_key][_idx] = np.concatenate(xy_test[_key][_idx]) - fit_kws = { - "y": y_train_fold[train_sub][train_sub_tgt_idx], - "groups": groups[train_sub], - "target_idx": train_sub_tgt_idx, - } + fit_kws = get_fit_kws( + y_train_fold[train_sub][train_sub_tgt_idx], groups[train_sub], train_sub_tgt_idx + ) model_filename = "HCP_L%s_test_size0%s_%s_%s_%s_group_%s_%s" % ( int(lambda_), str(int(test_size * 10)), @@ -448,18 +450,11 @@ def run_sub_hold_hcp( if mix_group: model_filename = model_filename + "_group_mix" - fit_kws = { - "y": y_train_fold[train_sub], - "groups": groups[train_sub], - "target_idx": None, - } - - init_kws = {"lambda_": lambda_, "l2_hparam": l2_hparam} - for arg in GSDA_INIT_ARGS: - if arg in kwargs: - init_kws[arg] = kwargs[arg] - - model = train_modal( + fit_kws = get_fit_kws(y_train_fold[train_sub], groups[train_sub]) + + init_kws = get_gsda_init_kws(lambda_, l2_hparam, kwargs) + + model = train_model( x_train, init_kws, fit_kws, @@ -510,8 +505,9 @@ def run_sub_hold_gsp( x_test_fold = x_all[1 - train_fold] y_test_fold = y_all[1 - train_fold] - for tgt_group in [0, 1]: - _idx = split_by_group(groups, train_sub, test_sub, tgt_group) + for tgt_group in get_target_groups(mix_group): + group_idx = 0 if mix_group else tgt_group + _idx = split_by_group(groups, train_sub, test_sub, group_idx) ( train_sub_tgt_idx, train_sub_nt_idx, @@ -519,21 +515,12 @@ def run_sub_hold_gsp( test_sub_nt_idx, ) = _idx - if mix_group: - if tgt_group == 1: - continue - tgt_group = "mix" - x_train_fold_test_tgt = x_test_fold[train_sub][train_sub_tgt_idx] y_train_fold_test_tgt = y_test_fold[train_sub][train_sub_tgt_idx] x_train_fold_test_nt = x_test_fold[train_sub][train_sub_nt_idx] y_train_fold_test_nt = y_test_fold[train_sub][train_sub_nt_idx] - fit_kws = { - "y": y_train_fold[train_sub_tgt_idx], - "groups": groups, - "target_idx": train_sub_tgt_idx, - } + fit_kws = get_fit_kws(y_train_fold[train_sub_tgt_idx], groups, train_sub_tgt_idx) model_filename = "GSP_L%s_test_size0%s_%s_%s_group_%s_%s" % ( int(lambda_), str(int(test_size * 10)), @@ -577,11 +564,9 @@ def run_sub_hold_gsp( ], } model_filename = model_filename + "_test_sub_0%s" % str(int(test_size * 10)) - fit_kws = { - "y": y_train_fold[train_sub][train_sub_tgt_idx], - "groups": groups[train_sub], - "target_idx": train_sub_tgt_idx, - } + fit_kws = get_fit_kws( + y_train_fold[train_sub][train_sub_tgt_idx], groups[train_sub], train_sub_tgt_idx + ) else: x_train = x_train_fold xy_test = { @@ -592,23 +577,12 @@ def run_sub_hold_gsp( if mix_group: model_filename = model_filename + "_group_mix" if 0 < test_size < 1: - fit_kws = { - "y": y_train_fold[train_sub], - "group": groups[train_sub], - "target_idx": None, - } + fit_kws = get_fit_kws(y_train_fold[train_sub], groups[train_sub]) else: - fit_kws = { - "y": y_train_fold, - "group": groups, - "target_idx": None, - } - init_kws = {"lambda_": lambda_, "l2_hparam": l2_hparam} - for arg in GSDA_INIT_ARGS: - if arg in kwargs: - init_kws[arg] = kwargs[arg] - - model = train_modal( + fit_kws = get_fit_kws(y_train_fold, groups) + init_kws = get_gsda_init_kws(lambda_, l2_hparam, kwargs) + + model = train_model( x_train, init_kws, fit_kws, @@ -654,25 +628,17 @@ def run_no_sub_hold_gsp( # scaler = StandardScaler() # scaler.fit(x_train) - for tgt_group in [0, 1]: - tgt_idx = np.where(groups == tgt_group)[0] - nt_idx = np.where(groups == 1 - tgt_group)[0] - - if mix_group: - if tgt_group == 1: - continue - tgt_group = "mix" + for tgt_group in get_target_groups(mix_group): + group_idx = 0 if mix_group else tgt_group + tgt_idx = np.where(groups == group_idx)[0] + nt_idx = np.where(groups == 1 - group_idx)[0] xy_test = { - "acc_tgt": [x_all[1 - i_fold][tgt_idx], x_all[1 - i_fold][tgt_idx]], - "acc_nt": [x_all[1 - i_fold][nt_idx], x_all[1 - i_fold][nt_idx]], + "acc_tgt": [x_all[1 - i_fold][tgt_idx], y_all[1 - i_fold][tgt_idx]], + "acc_nt": [x_all[1 - i_fold][nt_idx], y_all[1 - i_fold][nt_idx]], } - fit_kws = { - "y": y_train[tgt_idx], - "groups": groups, - "target_idx": tgt_idx, - } + fit_kws = get_fit_kws(y_train[tgt_idx], groups, tgt_idx) model_filename = "GSP_L%s_test_size00_%s_%s_group_%s_%s" % ( int(lambda_), i_split, @@ -683,17 +649,10 @@ def run_no_sub_hold_gsp( if mix_group: model_filename = model_filename + "_group_mix" - fit_kws = { - "y": y_train, - "groups": groups, - "target_idx": None, - } - init_kws = {"lambda_": lambda_, "l2_hparam": l2_hparam} - for arg in GSDA_INIT_ARGS: - if arg in kwargs: - init_kws[arg] = kwargs[arg] + fit_kws = get_fit_kws(y_train, groups) + init_kws = get_gsda_init_kws(lambda_, l2_hparam, kwargs) - model = train_modal( + model = train_model( x_train, init_kws, fit_kws, From c3a306ea427a29e50ade3af804128568d0b46b62 Mon Sep 17 00:00:00 2001 From: Shuo Zhou Date: Sat, 6 Jun 2026 23:52:54 +0100 Subject: [PATCH 4/6] refactor gsda example --- .../brain_lateralization/utils/data_io.py | 74 +++ .../brain_lateralization/utils/experiment.py | 5 +- .../brain_lateralization/utils/half_brain.py | 138 ++++ examples/brain_lateralization/utils/io_.py | 598 +++--------------- .../brain_lateralization/utils/neuro_io.py | 63 ++ examples/brain_lateralization/utils/plot.py | 203 +----- .../brain_lateralization/utils/plotting.py | 123 ++++ .../brain_lateralization/utils/results.py | 52 ++ examples/brain_lateralization/utils/stats.py | 32 + 9 files changed, 561 insertions(+), 727 deletions(-) create mode 100644 examples/brain_lateralization/utils/data_io.py create mode 100644 examples/brain_lateralization/utils/half_brain.py create mode 100644 examples/brain_lateralization/utils/neuro_io.py create mode 100644 examples/brain_lateralization/utils/plotting.py create mode 100644 examples/brain_lateralization/utils/results.py create mode 100644 examples/brain_lateralization/utils/stats.py diff --git a/examples/brain_lateralization/utils/data_io.py b/examples/brain_lateralization/utils/data_io.py new file mode 100644 index 0000000..05118fd --- /dev/null +++ b/examples/brain_lateralization/utils/data_io.py @@ -0,0 +1,74 @@ +import os + +import h5py +import numpy as np +import pandas as pd + + +def load_txt(fpaths, connection_type="intra"): + """Load left/right half-brain vectors from text connectivity matrices.""" + from .half_brain import split_functional_brain + + left_ = [] + right_ = [] + for fpath in fpaths: + data_matrix = np.genfromtxt(fpath) + left_vec, right_vec = split_functional_brain(data_matrix, connection_type=connection_type) + left_.append(left_vec.reshape((1, -1))) + right_.append(right_vec.reshape((1, -1))) + + return np.concatenate(left_, axis=0), np.concatenate(right_, axis=0) + + +def load_hdf5(fpath): + with h5py.File(fpath, "r") as f: + return {"Left": f["Left"][()], "Right": f["Right"][()]} + + +def read_tabular(fname, **kwargs): + """Read a table from a .xlsx or .csv file.""" + file_format = fname.split(".")[-1] + if file_format == "xlsx": + return pd.read_excel(fname, engine="openpyxl", **kwargs) + if file_format == "csv": + return pd.read_csv(fname, **kwargs) + + raise ValueError("Unsupported file type %s" % file_format) + + +def get_fpaths(fdir, idx_list, file_format="txt"): + """Return existing subject-indexed files under a directory.""" + fpaths = {} + for idx in idx_list: + fname = "%s.%s" % (idx, file_format) + fpath = os.path.join(fdir, fname) + if os.path.exists(fpath): + fpaths[idx] = fpath + + return pd.DataFrame(data={"File path": fpaths.values()}, index=list(fpaths.keys())) + + +def file_split(file_path): + filepath, tempfilename = os.path.split(file_path) + filename, extension = os.path.splitext(tempfilename) + return filepath, filename, extension + + +def read_file(file): + _, _, ext = file_split(file) + if ext == ".h5": + return pd.read_hdf(file) + if ext == ".pkl": + return pd.read_pickle(file) + + raise ValueError("Unsupported file type %s, supported type: xxx.h5 or xxx.pkl..." % ext) + + +def select_data_subj(df, subj_id, df_column_name): + return df[df_column_name][df["subject"] == subj_id] + + +def select_data_multi_subj(df, subj_ids, df_column_name): + subj_df = df["subject"].to_numpy() + _, index1, _ = np.intersect1d(subj_df, subj_ids, assume_unique=False, return_indices=True) + return df.loc[index1, df_column_name].to_numpy() diff --git a/examples/brain_lateralization/utils/experiment.py b/examples/brain_lateralization/utils/experiment.py index babff69..67ce933 100644 --- a/examples/brain_lateralization/utils/experiment.py +++ b/examples/brain_lateralization/utils/experiment.py @@ -11,10 +11,11 @@ from torch.hub import download_url_to_file sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), "../.."))) -from utils.io_ import load_half_brain, pick_half, read_tabular - from kalelinear.estimator import GSDA # , GSLRTorch +from .data_io import read_tabular +from .half_brain import load_half_brain, pick_half + BASE_RESULT_DICT: dict = { "pred_loss": [], "hsic_loss": [], diff --git a/examples/brain_lateralization/utils/half_brain.py b/examples/brain_lateralization/utils/half_brain.py new file mode 100644 index 0000000..1246a68 --- /dev/null +++ b/examples/brain_lateralization/utils/half_brain.py @@ -0,0 +1,138 @@ +import os +from urllib.request import urlretrieve + +import h5py +import numpy as np +from scipy.io import loadmat, savemat +from sklearn.model_selection import train_test_split +from sklearn.preprocessing import label_binarize + +from .data_io import load_hdf5 + +HCP_LINK = { + "REST1": "https://zenodo.org/records/10050233/files/HCP_BNA_intra_half_brain_REST1_Fisherz.hdf5", + "REST2": "https://zenodo.org/records/10050233/files/HCP_BNA_intra_half_brain_REST2_Fisherz.hdf5", +} +GSP_LINK = "https://zenodo.org/records/10050234/files/gsp_BNA_intra_half_brain_Fisherz.mat" + + +def download_url_to_file(url, fpath): + urlretrieve(url, fpath) + + +def split_functional_brain(matrix, connection_type="intra"): + n_ = matrix.shape[1] / 2 + if connection_type == "intra": + left = matrix[0::2, 0::2] + right = matrix[1::2, 1::2] + + idx = np.triu_indices(n_, k=1) + n_feat = int(n_ * (n_ - 1) / 2) + left_vec = np.zeros((1, n_feat)) + right_vec = np.zeros((1, n_feat)) + left_vec[0, :] = left[idx] + right_vec[0, :] = right[idx] + elif connection_type == "inter": + left = matrix[0::2, 1::2] + right = matrix[1::2, 0::2] + left_vec = left.reshape((1, -1)) + right_vec = right.reshape((1, -1)) + else: + raise ValueError("Invalid connection type %s" % connection_type) + + return left_vec, right_vec + + +def save_half_brain(out_dir, out_fname, data_left, data_right): + with h5py.File(os.path.join(out_dir, out_fname), "w") as f: + f.create_dataset("Left", data=data_left) + f.create_dataset("Right", data=data_right) + + +def save_half_brain_mat(out_dir, out_fname, data_left, data_right): + savemat(os.path.join(out_dir, out_fname), {"Left": data_left, "Right": data_right}) + + +def load_half_brain( + data_dir, + atlas, + session=None, + run=None, + connection_type="intra", + data_type="functional", + dataset="HCP", + download=True, +): + if dataset == "HCP": + if data_type == "functional": + fname = "HCP_%s_%s_half_brain_%s_%s.hdf5" % (atlas, connection_type, session, run) + fpath = os.path.join(data_dir, fname) + if not os.path.exists(fpath): + if download: + os.makedirs(data_dir, exist_ok=True) + print("Downloading %s session %s data, it may take 1-2 mins" % (dataset, session)) + download_url_to_file(HCP_LINK[session], fpath) + else: + raise ValueError("File %s does not exist" % fpath) + data = load_hdf5(fpath) + elif data_type == "structural": + data = {"Left": [], "Right": []} + fname = "%s_Volume.mat" % atlas + data_in = loadmat(os.path.join(data_dir, fname))["%s_Volume" % atlas][0][0][0] + data["Left"] = data_in[:, 0::2] + data["Right"] = data_in[:, 1::2] + else: + raise ValueError("Invalid data type %s" % data_type) + elif dataset in ["ABIDE", "ukb", "GSP"]: + fpath = os.path.join(data_dir, "%s_%s_%s_half_brain_%s.mat" % (dataset, atlas, connection_type, run)) + if not os.path.exists(fpath): + if download: + if dataset == "GSP": + os.makedirs(data_dir, exist_ok=True) + print("Downloading %s data, it may take 1-2 mins." % dataset) + download_url_to_file(GSP_LINK, fpath) + else: + raise ValueError("File %s does not exist" % fpath) + else: + raise ValueError("File %s does not exist" % fpath) + data_file = loadmat(fpath) + data = {"Left": data_file["Left"], "Right": data_file["Right"]} + else: + raise ValueError("Invalid dataset %s" % dataset) + + return data + + +def pick_half(data, random_state=144): + x = np.zeros(data["Left"].shape) + left_idx, right_idx = train_test_split(range(x.shape[0]), test_size=0.5, random_state=random_state) + x[left_idx] = data["Left"][left_idx] + x[right_idx] = data["Right"][right_idx] + + n_sub = x.shape[0] + y = np.zeros(n_sub) + y[left_idx] = 1 + y[right_idx] = -1 + + x1 = np.zeros(data["Left"].shape) + x1[left_idx] = data["Right"][left_idx] + x1[right_idx] = data["Left"][right_idx] + + y1 = np.zeros(n_sub) + y1[left_idx] = -1 + y1[right_idx] = 1 + + y = label_binarize(y, classes=[-1, 1]).reshape(-1) + y1 = label_binarize(y1, classes=[-1, 1]).reshape(-1) + + return x, y, x1, y1 + + +def _pick_half_subs(data, random_state=144): + n_ = data["Left"].shape[0] + train_idx, _ = train_test_split(range(n_), test_size=0.5, random_state=random_state) + x = np.concatenate([data["Left"][train_idx], data["Right"][train_idx]], axis=0) + y = np.ones(n_) + y[int(n_ / 2) :] = -1 + + return x, y diff --git a/examples/brain_lateralization/utils/io_.py b/examples/brain_lateralization/utils/io_.py index 47eb47d..e4abd7f 100644 --- a/examples/brain_lateralization/utils/io_.py +++ b/examples/brain_lateralization/utils/io_.py @@ -1,527 +1,71 @@ -import os - -import h5py -import matplotlib.pyplot as plt -import nibabel as nib -import numpy as np -import pandas as pd -import seaborn as sns -import torch -from scipy.io import loadmat, savemat -from scipy.stats import pearsonr -from sklearn.model_selection import train_test_split -from sklearn.preprocessing import label_binarize -from torch.hub import download_url_to_file - -HCP_LINK = { - "REST1": "https://zenodo.org/records/10050233/files/HCP_BNA_intra_half_brain_REST1_Fisherz.hdf5", - "REST2": "https://zenodo.org/records/10050233/files/HCP_BNA_intra_half_brain_REST2_Fisherz.hdf5", -} -GSP_LINK = "https://zenodo.org/records/10050234/files/gsp_BNA_intra_half_brain_Fisherz.mat" - - -def load_txt(fpaths, connection_type="intra"): - """Load data from a list of txt files。 - - Parameters - ---------- - fpaths (list): A list of .txt file paths. - - Returns - ------- - Two ndarrays - """ - left_ = [] - right_ = [] - for fpath in fpaths: - data_matrix = np.genfromtxt(fpath) - left_vec, right_vec = split_functional_brain(data_matrix, connection_type=connection_type) - left_.append(left_vec.reshape((1, -1))) - right_.append(right_vec.reshape((1, -1))) - left_ = np.concatenate(left_, axis=0) - right_ = np.concatenate(right_, axis=0) - - return left_, right_ - - -def load_hdf5(fpath): - """ - - Parameters - ---------- - fpath - - Returns - ------- - - """ - f = h5py.File(fpath, "r") - data = {"Left": f["Left"][()], "Right": f["Right"][()]} - return data - - -def read_tabular(fname, **kwargs): - """Read a table from a .xlsx or .csv file - - Parameters - ---------- - fname (string): - - Returns - ------- - - """ - file_format = fname.split(".")[-1] - if file_format == "xlsx": - df = pd.read_excel(fname, engine="openpyxl", **kwargs) - elif file_format == "csv": - df = pd.read_csv(fname, **kwargs) - else: - raise ValueError("Unsupported file type %s" % file_format) - - return df - - -def get_fpaths(fdir, idx_list, file_format="txt"): - """Check the existence of files named by subject indices under a given directory - - Parameters - ---------- - fdir (string): file directory for data - idx_list (list): a list of subject indices - file_format (string): file format, defaults to txt - - Returns - ------- - Dataframe - """ - fpaths = dict() - for idx in idx_list: - fname = "%s.%s" % (idx, file_format) - fpath = os.path.join(fdir, fname) - if os.path.exists(fpath): - fpaths[idx] = fpath - - fpath_df = pd.DataFrame(data={"File path": fpaths.values()}, index=list(fpaths.keys())) - - return fpath_df - - -def split_functional_brain(matrix, connection_type="intra"): - """ - - Parameters - ---------- - matrix - connection_type - - Returns - ------- - - """ - n_ = matrix.shape[1] / 2 - if connection_type == "intra": - left = matrix[0::2, 0::2] # even indices of rows, left brain - right = matrix[1::2, 1::2] # odd indices of rows, right brain - - idx = np.triu_indices(n_, k=1) - n_feat = int(n_ * (n_ - 1) / 2) - left_vec = np.zeros((1, n_feat)) - right_vec = np.zeros((1, n_feat)) - left_vec[0, :] = left[idx] - right_vec[0, :] = right[idx] - - elif connection_type == "inter": - left = matrix[0::2, 1::2] - right = matrix[1::2, 0::2] - left_vec = left.reshape((1, -1)) - right_vec = right.reshape((1, -1)) - - else: - raise ValueError("Invalid connection type %s" % connection_type) - - return left_vec, right_vec - - -def save_half_brain(out_dir, out_fname, data_left, data_right): - """Save two hemisphere nadarrays into a hdf5 file - - Parameters - ---------- - out_dir (string): - out_fname (string): - data_left (ndarray): - data_right (ndarray): - - Returns - ------- - - """ - f = h5py.File(os.path.join(out_dir, out_fname), "w") - f.create_dataset("Left", data=data_left) - f.create_dataset("Right", data=data_right) - f.close() - - -def save_half_brain_mat(out_dir, out_fname, data_left, data_right): - """Save two hemisphere nadarrays into a hdf5 file - - Parameters - ---------- - out_dir (string): - out_fname (string): - data_left (ndarray): - data_right (ndarray): - - Returns - ------- - - """ - savemat(os.path.join(out_dir, out_fname), {"Left": data_left, "Right": data_right}) - - -def load_half_brain( - data_dir, - atlas, - session=None, - run=None, - connection_type="intra", - data_type="functional", - dataset="HCP", - download=True, -): - """ - - Parameters - ---------- - data_dir - atlas - session - run - connection_type - data_type - dataset - download: bool, optional, download the data if the data files not exist (default=True) - - Returns - ------- - - """ - if dataset == "HCP": - if data_type == "functional": - fname = "HCP_%s_%s_half_brain_%s_%s.hdf5" % ( - atlas, - connection_type, - session, - run, - ) - fpath = os.path.join(data_dir, fname) - if not os.path.exists(fpath): - if download: - os.makedirs(data_dir, exist_ok=True) - print("Downloading %s session %s data, it may take 1-2 mins" % (dataset, session)) - download_url_to_file( - HCP_LINK[session], - fpath, - ) - else: - raise ValueError("File %s does not exist" % fpath) - data = load_hdf5(fpath) - elif data_type == "structural": - data = {"Left": [], "Right": []} - fname = "%s_Volume.mat" % atlas - data_in = loadmat(os.path.join(data_dir, fname))["%s_Volume" % atlas][0][0][0] - data["Left"] = data_in[:, 0::2] - data["Right"] = data_in[:, 1::2] - else: - raise ValueError("Invalid data type %s" % data_type) - elif dataset in ["ABIDE", "ukb", "GSP"]: - fpath = os.path.join( - data_dir, - "%s_%s_%s_half_brain_%s.mat" % (dataset, atlas, connection_type, run), - ) - if not os.path.exists(fpath): - if download: - if dataset == "GSP": - os.makedirs(data_dir, exist_ok=True) - print("Downloading %s data, it may take 1-2 mins." % dataset) - download_url_to_file(GSP_LINK, fpath) - else: - raise ValueError("File %s does not exist" % fpath) - else: - raise ValueError("File %s does not exist" % fpath) - data_file = loadmat(fpath) - data = {"Left": data_file["Left"], "Right": data_file["Right"]} - else: - raise ValueError("Invalid dataset %s" % dataset) - - return data - - -def fetch_weights(base_dir, group, lambda_, dataset, sessions, test_size="00", seed_=2023): - """ - - Args: - base_dir: - group: int, 0 or 1 - lambda_: str, "0_group_mix" or "0", "1", "2", "5", "8", "10" - dataset: str, "HCP" or "GSP" - sessions: list of session names, ["REST1_", "REST2_"] for HCP and ["_"] for GSP - test_size: str, "00" or "02", optional, default="00" - seed_: int, optional, default=2023 - - Returns: - a matrix of weights, shape (n_models, n_features) - """ - - sub_dir = os.path.join(base_dir, "lambda%s" % lambda_) - if not os.path.exists(sub_dir): - return None - if lambda_ == "0_group_mix": - lambda_ = 0 - group = "mix" - weight = [] - num_repeat = 5 - halfs = [0, 1] - for session_i in sessions: - for half_i in halfs: - for i_split in range(num_repeat): - for seed in range(52): - model_file_name = "%s_L%s_test_size%s_%s%s_%s_group_%s_%s.pt" % ( - dataset, - lambda_, - test_size, - session_i, - i_split, - half_i, - group, - seed_ - seed, - ) - if os.path.exists(os.path.join(sub_dir, model_file_name)): - weight.append(get_coef(model_file_name, sub_dir).reshape((1, -1))) - - return np.concatenate(weight, axis=0) - - -def get_coef(file_name, file_dir): - file_path = os.path.join(file_dir, file_name) - model = torch.load(file_path) - return model.theta - - -def save_results(res_dict, out_filename, output_dir, mix_group=False): - """save a dictionary to a csv file - - Args: - res_dict (_type_): _description_ - out_filename (str): _description_ - output_dir (str): _description_ - mix_group (bool, optional): _description_. Defaults to False. - """ - res_df = pd.DataFrame.from_dict(res_dict) - - if mix_group: - out_filename = out_filename + "_mix_group" - out_file = os.path.join(output_dir, "%s.csv" % out_filename) - res_df.to_csv(out_file, index=False) - - -# read: nib.load() -# write: nib.save() -# ********************************** file split ************************************* # -def file_split(file_path): - filepath, tempfilename = os.path.split(file_path) - filename, extension = os.path.splitext(tempfilename) - return filepath, filename, extension - - -# ********************************** read h5 or pickle file as a df ************************************* # - - -def read_file(file): - _, _, ext = file_split(file) - if ext == ".h5": - df = pd.read_hdf(file) - elif ext == ".pkl": - df = pd.read_pickle(file) - else: - raise ValueError("Unsupported file type %s, supported type: xxx.h5 or xxx.pkl..." % ext) - - return df - - -def select_data_subj(df, subj_id, df_column_name): - data = df[df_column_name][df["subject"] == subj_id] - return data # series datatype with dim of (1,),for the ndarray contents, data = data[0] - - -def pick_half(data, random_state=144): - x = np.zeros(data["Left"].shape) - left_idx, right_idx = train_test_split(range(x.shape[0]), test_size=0.5, random_state=random_state) - x[left_idx] = data["Left"][left_idx] - x[right_idx] = data["Right"][right_idx] - - n_sub = x.shape[0] - y = np.zeros(n_sub) - y[left_idx] = 1 - y[right_idx] = -1 - - x1 = np.zeros(data["Left"].shape) - x1[left_idx] = data["Right"][left_idx] - x1[right_idx] = data["Left"][right_idx] - - y1 = np.zeros(n_sub) - y1[left_idx] = -1 - y1[right_idx] = 1 - - y = label_binarize(y, classes=[-1, 1]).reshape(-1) - y1 = label_binarize(y1, classes=[-1, 1]).reshape(-1) - - return x, y, x1, y1 - - -def _pick_half_subs(data, random_state=144): - n_ = data["Left"].shape[0] - train_idx, hold_idx = train_test_split(range(n_), test_size=0.5, random_state=random_state) - x = np.concatenate([data["Left"][train_idx], data["Right"][train_idx]], axis=0) - y = np.ones(n_) - y[int(n_ / 2) :] = -1 - - return x, y - - -def select_data_multi_subj(df, subj_ids, df_column_name): - subj_df = df["subject"].to_numpy() - sub, index1, index2 = np.intersect1d(subj_df, subj_ids, assume_unique=False, return_indices=True) - data = df.loc[index1, df_column_name] # loc,select column name. iloc,select data as index - return data.to_numpy() # convert pd series to numpy ndarray - - -# ********************************** read MRI files ********************************* # -def read_surf(surface): - surf = nib.load(surface) - arr = surf.darrays - coord = arr[0].data # coordinate of vertices - vert = arr[1].data # index of vertices - return coord, vert - - -def read_nii(T1w): - T1 = nib.load(T1w) - data = T1.get_fdata() # 3-D ‘ray - header = T1.header # header - return data, header - - -def read_shape_gii(shape_gii): - gii = nib.load(shape_gii) - data = gii.darrays[0].data - return gii, data - - -# ********************************** Creat a .shape.gii file ********************************* # - - -def creat_shape_gii(shape_gii_base, array_write, savepath): - gii, _ = read_shape_gii(shape_gii_base) - gii.darrays[0].data = array_write - nib.save(gii, savepath) - - -# Creat a .shape gii file with atlas, e.g.,Kong400,Schaefer400 -def creat_fs_lr32k_atlas(shape_gii_base, array_write, atlas_file, hemi, savepath): - # cdata = np.zeros((32492,)) - f_atlas = nib.load(atlas_file) - f_data = f_atlas.get_fdata() # 0 mean medial wall, label starts from 1 - f_data = f_data[0, :] - - if hemi == "L": - cdata = f_data[0 : int(len(f_data) / 2)].copy() # L - print("cdata.shape: %d" % (len(cdata))) - for i, data in enumerate(array_write): - cdata[np.where(cdata == i + 1)] = data - else: - cdata = f_data[int(len(f_data) / 2) : int(len(f_data))].copy() # R - for i, data in enumerate(array_write): - cdata[np.where(cdata == i + 1 + int(np.max(f_data) / 2))] = data - - creat_shape_gii(shape_gii_base, cdata, savepath) - - -# *********************************************** Corr with pvalue ************************************************ # -def Corr(x, y): - r, p = pearsonr(x, y) - rr = round(r, 3) - pp = round(p, 3) - return rr, pp - - -# ************************************************ Join plot ***************************************************** # -# x,y should not be the object data type, if so, convert it to np.double data type. -def jointplot_fitlinear(x, y): - plt.figure(figsize=(20, 16)) - sns.jointplot(x=x, y=y, kind="reg", scatter_kws={"s": 8}) - - -# sns.regplot(x='LL_LPAC_Math_Corr',y='ReadEng_AgeAdj',data=df, color='red', scatter_kws={"s": 8}),不显示分布信息 - - -# *************************************** Extract top N in a array ************************************ # -# Extract top N element in an array and set others to 0. -def topN_array(arr, topN): - idx_all = set(np.arange(len(arr))) - idx_topN = arr.argsort()[::-1][:topN] - idx_res = idx_all - set(idx_topN) - # keep topN, set others to zero - idx_res = np.array(list(idx_res)) - arr_topN = arr.copy() - arr_topN[idx_res] = 0 - return arr_topN - - -# def main_creat_shape_gii(): -# -# tasks = ['LANGUAGE'] -# subjs = sio.loadmat(project_path + '/Subject_taskT.mat') -# subjects = subjs['Subject_taskT'] -# subjects = np.reshape(subjects,(len(subjects),)) -# for creat shape gii files -# shape_gii_base_L = project_path + '/Base_shape_gii/100206.L.thickness.32k_fs_LR.shape.gii' -# for creat shape gii files -# shape_gii_base_R = project_path + '/Base_shape_gii/100206.R.thickness.32k_fs_LR.shape.gii' -# -# for subj_id in subjects: -# df_column_name = 'task_t' -# for task in tasks: -# files = sorted(glob.glob(project_path + '/DataSort/Sort_Sub/' + task + '/*.h5')) -# -# for file in files: -# -# Dirs = os.path.split(os.path.splitext(file)[0]) -# kName = Dirs[1] -# df = rvp.read_orig_data(file, kName) -# ResultantFolder = (project_path + '/Result/Ridge_Vert_SubBased/Statistics/task_t_gii/' -# + str(subj_id) +'/' + task) -# if not os.path.exists(ResultantFolder): -# os.makedirs(ResultantFolder) -# task_t = rvp.select_data_subj(df, subj_id, 'task_t').to_numpy()[0] -# task_t = np.reshape(task_t, (len(task_t),)) -# vert_index = rvp.select_data_subj(df, subj_id, 'vert_index').to_numpy()[0] -# cdata = np.zeros((32492,)) -# vert_index = np.reshape(vert_index, (len(vert_index),)) -# cdata[vert_index - 1] = task_t # compared to matlab, the vert index should minus 1 -# -# save_name = ResultantFolder + '/' + str(subj_id) + '_' + kName + '.shape.gii' -# if 'LW' in file: -# print('LEFT SEED') -# shape_gii_base= shape_gii_base_L -# elif 'RW' in file: -# print('RIGHT SEED') -# shape_gii_base= shape_gii_base_R -# else: -# print('error!') -# -# creat_shape_gii(shape_gii_base, cdata, save_name) -# print(kName) -# -# print(str(subj_id)) +"""Compatibility imports for the brain lateralization utility modules. + +Prefer importing from the focused modules directly: +``data_io``, ``half_brain``, ``results``, ``neuro_io``, and ``stats``. +""" + +from .data_io import ( + file_split, + get_fpaths, + load_hdf5, + load_txt, + read_file, + read_tabular, + select_data_multi_subj, + select_data_subj, +) +from .half_brain import ( + _pick_half_subs, + GSP_LINK, + HCP_LINK, + load_half_brain, + pick_half, + save_half_brain, + save_half_brain_mat, + split_functional_brain, +) +from .neuro_io import ( + creat_fs_lr32k_atlas, + creat_shape_gii, + create_fs_lr32k_atlas, + create_shape_gii, + read_nii, + read_shape_gii, + read_surf, +) +from .results import fetch_weights, get_coef, save_results +from .stats import Corr, corr, jointplot_fitlinear, top_n_array, topN_array + +__all__ = [ + "Corr", + "GSP_LINK", + "HCP_LINK", + "_pick_half_subs", + "corr", + "creat_fs_lr32k_atlas", + "creat_shape_gii", + "create_fs_lr32k_atlas", + "create_shape_gii", + "fetch_weights", + "file_split", + "get_coef", + "get_fpaths", + "jointplot_fitlinear", + "load_half_brain", + "load_hdf5", + "load_txt", + "pick_half", + "read_file", + "read_nii", + "read_shape_gii", + "read_surf", + "read_tabular", + "save_half_brain", + "save_half_brain_mat", + "save_results", + "select_data_multi_subj", + "select_data_subj", + "split_functional_brain", + "topN_array", + "top_n_array", +] diff --git a/examples/brain_lateralization/utils/neuro_io.py b/examples/brain_lateralization/utils/neuro_io.py new file mode 100644 index 0000000..5d57edd --- /dev/null +++ b/examples/brain_lateralization/utils/neuro_io.py @@ -0,0 +1,63 @@ +import numpy as np + + +def read_surf(surface): + import nibabel as nib + + surf = nib.load(surface) + arr = surf.darrays + coord = arr[0].data + vert = arr[1].data + return coord, vert + + +def read_nii(T1w): + import nibabel as nib + + T1 = nib.load(T1w) + data = T1.get_fdata() + header = T1.header + return data, header + + +def read_shape_gii(shape_gii): + import nibabel as nib + + gii = nib.load(shape_gii) + data = gii.darrays[0].data + return gii, data + + +def create_shape_gii(shape_gii_base, array_write, savepath): + import nibabel as nib + + gii, _ = read_shape_gii(shape_gii_base) + gii.darrays[0].data = array_write + nib.save(gii, savepath) + + +def create_fs_lr32k_atlas(shape_gii_base, array_write, atlas_file, hemi, savepath): + import nibabel as nib + + f_atlas = nib.load(atlas_file) + f_data = f_atlas.get_fdata()[0, :] + + if hemi == "L": + cdata = f_data[0 : int(len(f_data) / 2)].copy() + print("cdata.shape: %d" % (len(cdata))) + for i, data in enumerate(array_write): + cdata[np.where(cdata == i + 1)] = data + else: + cdata = f_data[int(len(f_data) / 2) : int(len(f_data))].copy() + for i, data in enumerate(array_write): + cdata[np.where(cdata == i + 1 + int(np.max(f_data) / 2))] = data + + create_shape_gii(shape_gii_base, cdata, savepath) + + +def creat_shape_gii(shape_gii_base, array_write, savepath): + return create_shape_gii(shape_gii_base, array_write, savepath) + + +def creat_fs_lr32k_atlas(shape_gii_base, array_write, atlas_file, hemi, savepath): + return create_fs_lr32k_atlas(shape_gii_base, array_write, atlas_file, hemi, savepath) diff --git a/examples/brain_lateralization/utils/plot.py b/examples/brain_lateralization/utils/plot.py index c980aa0..dd56220 100644 --- a/examples/brain_lateralization/utils/plot.py +++ b/examples/brain_lateralization/utils/plot.py @@ -1,201 +1,8 @@ -import matplotlib.pylab as plt -import numpy as np -import pandas as pd -import seaborn as sns -from scipy.io import loadmat -from utils.io_ import fetch_weights +"""Compatibility wrapper for plotting utilities. +Prefer importing from ``utils.plotting`` in new code. +""" -def savefig(fig, outfile, outfig_format): - if isinstance(outfig_format, list): - for fmt in outfig_format: - fig.savefig("%s.%s" % (outfile, fmt), format=fmt, bbox_inches="tight") - else: - fig.savefig( - "%s.%s" % (outfile, outfig_format), - format=outfig_format, - bbox_inches="tight", - ) +from .plotting import load_coef_plot_corr, load_weight_plot_corr, savefig - -# def plot_gsi( -# data, -# x="lambda", -# y="GSI", -# hue="train_gender", -# fontsize=14, -# outfile=None, -# outfig_format=None, -# ): -# fig = plt.figure() -# plt.rcParams.update({"font.size": fontsize}) -# sns.boxplot(data=data, x=x, y=y, hue=hue, showmeans=True) -# plt.legend(title="Target group") -# plt.rcParams["text.usetex"] = True -# plt.xlabel(r"$\lambda$") -# plt.rcParams["text.usetex"] = False -# plt.ylabel("Group Specificity Index (GSI)") -# if outfile is not None: -# savefig(fig, outfile, outfig_format) -# plt.show() -# -# -# def plot_accuracy( -# data, -# x="Lambda", -# y="Accuracy", -# col="Target group", -# hue="Test set", -# style="Test set", -# kind="line", -# height=4, -# fontsize=14, -# outfile=None, -# outfig_format=None, -# ): -# fig = plt.figure() -# plt.rcParams.update({"font.size": fontsize}) -# g = sns.relplot( -# data=data, -# x=x, -# y=y, -# col=col, -# hue=hue, -# style=style, -# kind=kind, -# errorbar=("sd", 1), -# height=height, -# ) -# ( -# g.map(plt.axhline, y=0.9, color=".7", dashes=(2, 1), zorder=0) -# .set_axis_labels(r"$\lambda$", "Test Accuracy") -# .set_titles("Target group: {col_name}") -# .tight_layout(w_pad=0) -# ) -# if outfile is not None: -# savefig(fig, outfile, outfig_format) -# -# plt.show() - - -def load_weight_plot_corr(dataset, base_dir, sessions, seed_start, fontsize=14): - control_weights = fetch_weights(base_dir, "mix", "0_group_mix", dataset, sessions=sessions, seed_=seed_start) - n_control_weights = control_weights.shape[0] - - lambdas = [0.0, 1.0, 2.0, 5.0, 8.0, 10.0] - corrs = {"mean": [], "sd": []} - - for lambda_ in lambdas: - for group in [0, 1]: - weights = fetch_weights( - base_dir, - group, - int(lambda_), - dataset, - sessions=sessions, - seed_=seed_start, - ) - # n_weights = weights.shape[0] - corr_matrix = np.corrcoef(control_weights, weights)[n_control_weights:, :n_control_weights] - corrs["mean"].append(np.mean(corr_matrix)) - corrs["sd"].append(np.std(corr_matrix)) - - plt.rcParams.update({"font.size": fontsize}) - fig, ax = plt.subplots(2, 1, sharex=True) - fig.set_size_inches(4, 8.5) - # .plot(x, y1) - # .plot(x, y2) - - # plt.figure(figsize=(4, 4)) - ax[0].plot(lambdas, corrs["mean"][::2], "-", c="steelblue", label="Male-specific") - ax[0].fill_between( - lambdas, - np.asarray(corrs["mean"][::2]) - np.asarray(corrs["sd"][::2]), - np.asarray(corrs["mean"][::2]) + np.asarray(corrs["sd"][::2]), - color="steelblue", - alpha=0.2, - ) - ax[0].set_ylabel("Correlation") - ax[0].legend(loc="upper right") - - ax[1].plot(lambdas, corrs["mean"][1::2], "--", color="firebrick", label="Female-specific") - ax[1].fill_between( - lambdas, - np.asarray(corrs["mean"][1::2]) - np.asarray(corrs["sd"][1::2]), - np.asarray(corrs["mean"][1::2]) + np.asarray(corrs["sd"][1::2]), - color="firebrick", - alpha=0.2, - ) - - ax[1].set_ylabel("Correlation") - ax[1].legend(loc="upper right") - plt.rcParams["text.usetex"] = True - plt.xlabel(r"$\lambda$", fontsize=fontsize) - plt.rcParams["text.usetex"] = False - - plt.savefig("figures/%s_corr.svg" % dataset, format="svg", bbox_inches="tight") - # plt.savefig('figures/%s_corr.pdf' % dataset, format='pdf', bbox_inches='tight') - # plt.savefig('figures/%s_corr.png' % dataset, format='png', bbox_inches='tight') - plt.show() - - -def load_coef_plot_corr(dataset, group_label, fontsize=14): - weights_dirs = dict() - lambdas = [0.0, 1.0, 2.0, 5.0, 8.0, 10.0] - group_dict = {0: "Male", 1: "Female"} - for lambda_ in lambdas: - weights_dirs[r"$\lambda=%s$" % (int(lambda_))] = "%s/%s_L%sG%s.mat" % ( - dataset, - dataset, - int(lambda_), - group_label, - ) - - weights = dict() - - for key in weights_dirs: - # print(key) - weights[key] = loadmat(weights_dirs[key])["mean"][0][1:] - - weight_df = pd.DataFrame(weights) - - corr = weight_df.corr() - - # Generate a mask for the upper triangle - mask = np.triu(np.ones_like(corr, dtype=bool), 1) - - plt.rcParams.update({"font.size": fontsize}) - # Set up the matplotlib figure - f, ax = plt.subplots(figsize=(11, 9)) - - # Generate a custom diverging colormap - cmap = sns.diverging_palette(230, 20, as_cmap=True) - plt.rcParams["text.usetex"] = True - # Draw the heatmap with the mask and correct aspect ratio - sns.heatmap( - corr, - mask=mask, - cmap=cmap, - vmin=0, - vmax=1, - center=0.5, - annot=True, - annot_kws={"fontsize": "xx-large"}, - square=True, - linewidths=0.5, - cbar_kws={"shrink": 0.5}, - ) - # cbar_kws={"shrink": .5, "use_gridspec": False, "location": "top"}) #, annot_kws={"rotation": 45}) - plt.rcParams["text.usetex"] = False - ax.set_xticklabels(weight_df.columns.to_list(), fontsize=fontsize + 2) - ax.set_yticklabels(weight_df.columns.to_list(), rotation=45, ha="right", fontsize=fontsize + 2) - # plt.savefig('corr.pdf', format='pdf', bbox_inches='tight') - # plt.savefig('corr.png', format='png', bbox_inches='tight') - plt.savefig( - "figures/corr_annot_%s_%s.svg" % (dataset, group_dict[group_label]), - format="svg", - bbox_inches="tight", - ) - # plt.savefig('figures/corr_annot_%s_%s.pdf' % (dataset, group_dict[sex_label]), format='pdf', bbox_inches='tight') - # plt.savefig('figures/corr_annot_%s_%s.png' % (dataset, group_dict[sex_label]), format='png', bbox_inches='tight') - plt.show() +__all__ = ["load_coef_plot_corr", "load_weight_plot_corr", "savefig"] diff --git a/examples/brain_lateralization/utils/plotting.py b/examples/brain_lateralization/utils/plotting.py new file mode 100644 index 0000000..6d6ab54 --- /dev/null +++ b/examples/brain_lateralization/utils/plotting.py @@ -0,0 +1,123 @@ +import matplotlib.pylab as plt +import numpy as np +import pandas as pd +import seaborn as sns +from scipy.io import loadmat + +from .results import fetch_weights + + +def savefig(fig, outfile, outfig_format): + if isinstance(outfig_format, list): + for fmt in outfig_format: + fig.savefig("%s.%s" % (outfile, fmt), format=fmt, bbox_inches="tight") + else: + fig.savefig( + "%s.%s" % (outfile, outfig_format), + format=outfig_format, + bbox_inches="tight", + ) + + +def load_weight_plot_corr(dataset, base_dir, sessions, seed_start, fontsize=14): + control_weights = fetch_weights(base_dir, "mix", "0_group_mix", dataset, sessions=sessions, seed_=seed_start) + n_control_weights = control_weights.shape[0] + + lambdas = [0.0, 1.0, 2.0, 5.0, 8.0, 10.0] + corrs = {"mean": [], "sd": []} + + for lambda_ in lambdas: + for group in [0, 1]: + weights = fetch_weights( + base_dir, + group, + int(lambda_), + dataset, + sessions=sessions, + seed_=seed_start, + ) + corr_matrix = np.corrcoef(control_weights, weights)[n_control_weights:, :n_control_weights] + corrs["mean"].append(np.mean(corr_matrix)) + corrs["sd"].append(np.std(corr_matrix)) + + plt.rcParams.update({"font.size": fontsize}) + fig, ax = plt.subplots(2, 1, sharex=True) + fig.set_size_inches(4, 8.5) + + ax[0].plot(lambdas, corrs["mean"][::2], "-", c="steelblue", label="Male-specific") + ax[0].fill_between( + lambdas, + np.asarray(corrs["mean"][::2]) - np.asarray(corrs["sd"][::2]), + np.asarray(corrs["mean"][::2]) + np.asarray(corrs["sd"][::2]), + color="steelblue", + alpha=0.2, + ) + ax[0].set_ylabel("Correlation") + ax[0].legend(loc="upper right") + + ax[1].plot(lambdas, corrs["mean"][1::2], "--", color="firebrick", label="Female-specific") + ax[1].fill_between( + lambdas, + np.asarray(corrs["mean"][1::2]) - np.asarray(corrs["sd"][1::2]), + np.asarray(corrs["mean"][1::2]) + np.asarray(corrs["sd"][1::2]), + color="firebrick", + alpha=0.2, + ) + + ax[1].set_ylabel("Correlation") + ax[1].legend(loc="upper right") + plt.rcParams["text.usetex"] = True + plt.xlabel(r"$\lambda$", fontsize=fontsize) + plt.rcParams["text.usetex"] = False + + plt.savefig("figures/%s_corr.svg" % dataset, format="svg", bbox_inches="tight") + plt.show() + + +def load_coef_plot_corr(dataset, group_label, fontsize=14): + weights_dirs = {} + lambdas = [0.0, 1.0, 2.0, 5.0, 8.0, 10.0] + group_dict = {0: "Male", 1: "Female"} + for lambda_ in lambdas: + weights_dirs[r"$\lambda=%s$" % (int(lambda_))] = "%s/%s_L%sG%s.mat" % ( + dataset, + dataset, + int(lambda_), + group_label, + ) + + weights = {} + for key in weights_dirs: + weights[key] = loadmat(weights_dirs[key])["mean"][0][1:] + + weight_df = pd.DataFrame(weights) + corr = weight_df.corr() + mask = np.triu(np.ones_like(corr, dtype=bool), 1) + + plt.rcParams.update({"font.size": fontsize}) + _, ax = plt.subplots(figsize=(11, 9)) + + cmap = sns.diverging_palette(230, 20, as_cmap=True) + plt.rcParams["text.usetex"] = True + sns.heatmap( + corr, + mask=mask, + cmap=cmap, + vmin=0, + vmax=1, + center=0.5, + annot=True, + annot_kws={"fontsize": "xx-large"}, + square=True, + linewidths=0.5, + cbar_kws={"shrink": 0.5}, + ) + plt.rcParams["text.usetex"] = False + ax.set_xticklabels(weight_df.columns.to_list(), fontsize=fontsize + 2) + ax.set_yticklabels(weight_df.columns.to_list(), rotation=45, ha="right", fontsize=fontsize + 2) + plt.savefig( + "figures/corr_annot_%s_%s.svg" % (dataset, group_dict[group_label]), + format="svg", + bbox_inches="tight", + ) + plt.show() diff --git a/examples/brain_lateralization/utils/results.py b/examples/brain_lateralization/utils/results.py new file mode 100644 index 0000000..d7bf520 --- /dev/null +++ b/examples/brain_lateralization/utils/results.py @@ -0,0 +1,52 @@ +import os + +import numpy as np +import pandas as pd + + +def fetch_weights(base_dir, group, lambda_, dataset, sessions, test_size="00", seed_=2023): + sub_dir = os.path.join(base_dir, "lambda%s" % lambda_) + if not os.path.exists(sub_dir): + return None + if lambda_ == "0_group_mix": + lambda_ = 0 + group = "mix" + + weight = [] + num_repeat = 5 + halfs = [0, 1] + for session_i in sessions: + for half_i in halfs: + for i_split in range(num_repeat): + for seed in range(52): + model_file_name = "%s_L%s_test_size%s_%s%s_%s_group_%s_%s.pt" % ( + dataset, + lambda_, + test_size, + session_i, + i_split, + half_i, + group, + seed_ - seed, + ) + if os.path.exists(os.path.join(sub_dir, model_file_name)): + weight.append(get_coef(model_file_name, sub_dir).reshape((1, -1))) + + return np.concatenate(weight, axis=0) + + +def get_coef(file_name, file_dir): + import torch + + file_path = os.path.join(file_dir, file_name) + model = torch.load(file_path) + return model.theta + + +def save_results(res_dict, out_filename, output_dir, mix_group=False): + res_df = pd.DataFrame.from_dict(res_dict) + + if mix_group: + out_filename = out_filename + "_mix_group" + out_file = os.path.join(output_dir, "%s.csv" % out_filename) + res_df.to_csv(out_file, index=False) diff --git a/examples/brain_lateralization/utils/stats.py b/examples/brain_lateralization/utils/stats.py new file mode 100644 index 0000000..3cdc633 --- /dev/null +++ b/examples/brain_lateralization/utils/stats.py @@ -0,0 +1,32 @@ +import matplotlib.pyplot as plt +import numpy as np +import seaborn as sns +from scipy.stats import pearsonr + + +def corr(x, y): + r, p = pearsonr(x, y) + return round(r, 3), round(p, 3) + + +def Corr(x, y): + return corr(x, y) + + +def jointplot_fitlinear(x, y): + plt.figure(figsize=(20, 16)) + sns.jointplot(x=x, y=y, kind="reg", scatter_kws={"s": 8}) + + +def top_n_array(arr, top_n): + idx_all = set(np.arange(len(arr))) + idx_top_n = arr.argsort()[::-1][:top_n] + idx_res = idx_all - set(idx_top_n) + idx_res = np.array(list(idx_res)) + arr_top_n = arr.copy() + arr_top_n[idx_res] = 0 + return arr_top_n + + +def topN_array(arr, topN): + return top_n_array(arr, topN) From 61f838e4b7747db9faea414e10a39566d9c57532 Mon Sep 17 00:00:00 2001 From: Shuo Zhou Date: Mon, 8 Jun 2026 23:33:44 +0100 Subject: [PATCH 5/6] Fix import issues and typos in exmaples --- .../brain_lateralization/configs/test.yaml | 4 +- examples/brain_lateralization/tutorial.ipynb | 2 +- .../brain_lateralization/utils/data_io.py | 100 ++++++++++++++++++ .../brain_lateralization/utils/experiment.py | 2 - .../brain_lateralization/utils/half_brain.py | 4 +- examples/brain_lateralization/utils/io_.py | 4 + examples/brain_lateralization/utils/plot.py | 4 +- .../brain_lateralization/utils/plotting.py | 84 +++++++++++++++ .../brain_lateralization/utils/results.py | 7 +- examples/cmri_mpca/config.py | 2 +- examples/cmri_mpca/main.py | 2 +- examples/cmri_mpca/tutorial.ipynb | 2 +- examples/multisite_neuroimg_adapt/README.md | 14 ++- examples/multisite_neuroimg_adapt/config.py | 4 +- .../configs/tutorial.yaml | 4 +- examples/multisite_neuroimg_adapt/main.py | 2 +- .../multisite_neuroimg_adapt/tutorial.ipynb | 2 +- examples/toy_domain_adaptation/README.md | 13 +-- examples/toy_domain_adaptation/main.py | 4 +- setup.py | 1 + 20 files changed, 222 insertions(+), 39 deletions(-) diff --git a/examples/brain_lateralization/configs/test.yaml b/examples/brain_lateralization/configs/test.yaml index c8c5e71..45d00f5 100644 --- a/examples/brain_lateralization/configs/test.yaml +++ b/examples/brain_lateralization/configs/test.yaml @@ -1,7 +1,7 @@ DATASET: - ROOT: "D:/ShareFolder/BNA/Proc" + ROOT: "data" NUM_REPEAT: 5 SOLVER: SEED: 2022 OUTPUT: - ROOT: 'D:/ShareFolder/BNA/Result' + ROOT: 'output' diff --git a/examples/brain_lateralization/tutorial.ipynb b/examples/brain_lateralization/tutorial.ipynb index 2ac58fc..d6c2330 100644 --- a/examples/brain_lateralization/tutorial.ipynb +++ b/examples/brain_lateralization/tutorial.ipynb @@ -59,7 +59,7 @@ "source": [ "### Configurations\n", "\n", - "The customized configuration used in this demo is stored in `configs/demoh-cp.yaml`, this file overwrites defaults in `default_cfg.py` where a value is specified. Change the configuration file path to `cfg_path = \"configs/demo-gsp.yaml\"` for running the demo with GSP data." + "The customized configuration used in this demo is stored in `configs/demo-cp.yaml`, this file overwrites defaults in `default_cfg.py` where a value is specified. Change the configuration file path to `cfg_path = \"configs/demo-gsp.yaml\"` for running the demo with GSP data." ], "cell_type": "markdown" }, diff --git a/examples/brain_lateralization/utils/data_io.py b/examples/brain_lateralization/utils/data_io.py index 05118fd..9bba2ca 100644 --- a/examples/brain_lateralization/utils/data_io.py +++ b/examples/brain_lateralization/utils/data_io.py @@ -72,3 +72,103 @@ def select_data_multi_subj(df, subj_ids, df_column_name): subj_df = df["subject"].to_numpy() _, index1, _ = np.intersect1d(subj_df, subj_ids, assume_unique=False, return_indices=True) return df.loc[index1, df_column_name].to_numpy() + + +def load_result(dataset, root_dir, lambdas, seed_start, test_size=0.0): + """load brain left/right classification results for a dataset + + Args: + dataset (string): _description_ + root_dir (string): _description_ + lambdas (list): _description_ + seed_start (_type_): _description_ + test_size (float, optional): _description_. Defaults to 0.0. + """ + res_dict = dict() + res_list = [] + test_size_str = str(int(test_size * 10)) + for lambda_ in lambdas: + res_dict[lambda_] = [] + + for lambda_ in lambdas: + if not isinstance(lambda_, str): + lambda_str = str(int(lambda_)) + else: + lambda_str = lambda_ + model_dir = os.path.join(root_dir, "lambda%s" % lambda_str) + for seed_iter in range(51): + random_state = seed_start - seed_iter + res_fname = "results_%s_L%s_test_size0%s_Fisherz_%s.csv" % ( + dataset, + lambda_str, + test_size_str, + random_state, + ) + res_fpath = os.path.join(model_dir, res_fname) + if os.path.exists(res_fpath): + res_df = pd.read_csv(os.path.join(model_dir, res_fname)) + res_df["seed"] = random_state + res_dict[lambda_].append(res_df) + res_list.append(res_df) + + for lambda_ in lambdas: + res_dict[lambda_] = pd.concat(res_dict[lambda_]) + + res_df_all = pd.concat(res_list) + res_df_all = res_df_all.reset_index(drop=True) + + return res_df_all + + +def reformat_results(res_df, test_sets, male_label=0): + """reformat results dataframe to one accuracy per row + + Args: + res_df (_type_): _description_ + test_sets (_type_): _description_ + male_label (int, optional): _description_. Defaults to 0. + + Returns: + _type_: _description_ + """ + res_reformat = { + "Accuracy": [], + "Test set": [], + "Lambda": [], + "Target group": [], + "seed": [], + "split": [], + "fold": [], + "Train session": [], + } + for idx_ in res_df.index: + # print(idx_, res_df.iloc[idx_, 12]) + subset_ = res_df.loc[idx_, :] + for test_set in test_sets: + res_reformat["Accuracy"].append(subset_[test_set]) + res_reformat["Lambda"].append(subset_["lambda"]) + if "train_session" in subset_: + res_reformat["Train session"].append(subset_["train_session"]) + else: + res_reformat["Train session"].append(None) + if "target_group" not in subset_: # for GSP dataset + _group = subset_["train_gender"] + else: + _group = subset_["target_group"] + test_set_list = test_set.split("_") + if _group == "Male" or _group == male_label: + res_reformat["Target group"].append("Male") + if "oc" in test_set_list or "tgt" in test_set_list: + res_reformat["Test set"].append("Female") + else: + res_reformat["Test set"].append("Male") + else: + res_reformat["Target group"].append("Female") + if "oc" in test_set_list or "tgt" in test_set_list: + res_reformat["Test set"].append("Male") + else: + res_reformat["Test set"].append("Female") + + for key in ["seed", "split", "fold"]: + res_reformat[key].append(subset_[key]) + return pd.DataFrame(res_reformat) diff --git a/examples/brain_lateralization/utils/experiment.py b/examples/brain_lateralization/utils/experiment.py index 67ce933..da394d5 100644 --- a/examples/brain_lateralization/utils/experiment.py +++ b/examples/brain_lateralization/utils/experiment.py @@ -1,6 +1,5 @@ import copy import os -import sys # import pickle import numpy as np @@ -10,7 +9,6 @@ from sklearn.model_selection import StratifiedShuffleSplit from torch.hub import download_url_to_file -sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), "../.."))) from kalelinear.estimator import GSDA # , GSLRTorch from .data_io import read_tabular diff --git a/examples/brain_lateralization/utils/half_brain.py b/examples/brain_lateralization/utils/half_brain.py index 1246a68..2742662 100644 --- a/examples/brain_lateralization/utils/half_brain.py +++ b/examples/brain_lateralization/utils/half_brain.py @@ -21,13 +21,13 @@ def download_url_to_file(url, fpath): def split_functional_brain(matrix, connection_type="intra"): - n_ = matrix.shape[1] / 2 + n_ = matrix.shape[1] // 2 if connection_type == "intra": left = matrix[0::2, 0::2] right = matrix[1::2, 1::2] idx = np.triu_indices(n_, k=1) - n_feat = int(n_ * (n_ - 1) / 2) + n_feat = int(n_ * (n_ - 1) // 2) left_vec = np.zeros((1, n_feat)) right_vec = np.zeros((1, n_feat)) left_vec[0, :] = left[idx] diff --git a/examples/brain_lateralization/utils/io_.py b/examples/brain_lateralization/utils/io_.py index e4abd7f..87a100e 100644 --- a/examples/brain_lateralization/utils/io_.py +++ b/examples/brain_lateralization/utils/io_.py @@ -8,9 +8,11 @@ file_split, get_fpaths, load_hdf5, + load_result, load_txt, read_file, read_tabular, + reformat_results, select_data_multi_subj, select_data_subj, ) @@ -53,6 +55,7 @@ "jointplot_fitlinear", "load_half_brain", "load_hdf5", + "load_result", "load_txt", "pick_half", "read_file", @@ -60,6 +63,7 @@ "read_shape_gii", "read_surf", "read_tabular", + "reformat_results", "save_half_brain", "save_half_brain_mat", "save_results", diff --git a/examples/brain_lateralization/utils/plot.py b/examples/brain_lateralization/utils/plot.py index dd56220..d3c5931 100644 --- a/examples/brain_lateralization/utils/plot.py +++ b/examples/brain_lateralization/utils/plot.py @@ -3,6 +3,6 @@ Prefer importing from ``utils.plotting`` in new code. """ -from .plotting import load_coef_plot_corr, load_weight_plot_corr, savefig +from .plotting import load_coef_plot_corr, load_weight_plot_corr, plot_accuracy, plot_gsi, savefig -__all__ = ["load_coef_plot_corr", "load_weight_plot_corr", "savefig"] +__all__ = ["load_coef_plot_corr", "load_weight_plot_corr", "savefig", "plot_accuracy", "plot_gsi"] diff --git a/examples/brain_lateralization/utils/plotting.py b/examples/brain_lateralization/utils/plotting.py index 6d6ab54..8a5001e 100644 --- a/examples/brain_lateralization/utils/plotting.py +++ b/examples/brain_lateralization/utils/plotting.py @@ -121,3 +121,87 @@ def load_coef_plot_corr(dataset, group_label, fontsize=14): bbox_inches="tight", ) plt.show() + + +def plot_gsi( + data, + x="lambda", + y="GSI", + hue="train_gender", + fontsize=14, + outfile=None, + outfig_format=None, +): + """plot GSI boxplot + Args: + data (pd.DataFrame): data to plot + x (str): x-axis variable + y (str): y-axis variable + hue (str): hue variable + fontsize (int): font size + outfile (str): output file name without extension + outfig_format (str or list): output file format(s), e.g., 'pdf', 'svg', or ['pdf', 'svg'] + """ + fig = plt.figure() + plt.rcParams.update({"font.size": fontsize}) + sns.boxplot(data=data, x=x, y=y, hue=hue, showmeans=True) + plt.legend(title="Target group") + plt.rcParams["text.usetex"] = True + plt.xlabel(r"$\lambda$") + plt.rcParams["text.usetex"] = False + plt.ylabel("Group Specificity Index (GSI)") + if outfile is not None: + savefig(fig, outfile, outfig_format) + plt.show() + + +def plot_accuracy( + data, + x="Lambda", + y="Accuracy", + col="Target group", + hue="Test set", + style="Test set", + kind="line", + height=4, + fontsize=14, + outfile=None, + outfig_format=None, +): + """plot accuracy line plot + Args: + data (pd.DataFrame): data to plot + x (str): x-axis variable + y (str): y-axis variable + col (str): column variable + hue (str): hue variable + style (str): style variable + kind (str): kind of plot + height (int): height of the plot + fontsize (int): font size + outfile (str): output file name without extension + outfig_format (str or list): output file format(s), e.g., 'pdf', 'svg', or ['pdf', 'svg'] + """ + fig = plt.figure() + plt.rcParams.update({"font.size": fontsize}) + g = sns.relplot( + data=data, + x=x, + y=y, + col=col, + hue=hue, + style=style, + kind=kind, + errorbar=("sd", 1), + height=height, + ) + ( + g.map(plt.axhline, y=0.9, color=".7", dashes=(2, 1), zorder=0) + .set_axis_labels(r"$\lambda$", "Test Accuracy") + .set_titles("Target group: {col_name}") + .tight_layout(w_pad=0) + ) + if outfile is not None: + savefig(fig, outfile, outfig_format) + + plt.show() diff --git a/examples/brain_lateralization/utils/results.py b/examples/brain_lateralization/utils/results.py index d7bf520..0cc2d68 100644 --- a/examples/brain_lateralization/utils/results.py +++ b/examples/brain_lateralization/utils/results.py @@ -29,9 +29,10 @@ def fetch_weights(base_dir, group, lambda_, dataset, sessions, test_size="00", s group, seed_ - seed, ) - if os.path.exists(os.path.join(sub_dir, model_file_name)): - weight.append(get_coef(model_file_name, sub_dir).reshape((1, -1))) - + # if os.path.exists(os.path.join(sub_dir, model_file_name)): + # weight.append(get_coef(model_file_name, sub_dir).reshape((1, -1))) + if not weight: + return None return np.concatenate(weight, axis=0) diff --git a/examples/cmri_mpca/config.py b/examples/cmri_mpca/config.py index c2e41e8..e44f60c 100644 --- a/examples/cmri_mpca/config.py +++ b/examples/cmri_mpca/config.py @@ -18,7 +18,7 @@ _C.DATASET.ROOT = "../data" _C.DATASET.IMG_DIR = "DICOM" _C.DATASET.BASE_DIR = "SA_64x64_v2.0" -_C.DATASET.FILE_FORAMT = "zip" +_C.DATASET.FILE_FORMAT = "zip" _C.DATASET.LANDMARK_FILE = "landmarks.csv" _C.DATASET.MASK_DIR = "Mask" # ---------------------------------------------------------------------------- # diff --git a/examples/cmri_mpca/main.py b/examples/cmri_mpca/main.py index f6014d3..ccfb3be 100644 --- a/examples/cmri_mpca/main.py +++ b/examples/cmri_mpca/main.py @@ -51,7 +51,7 @@ def main(): # ---- setup dataset ---- base_dir = cfg.DATASET.BASE_DIR - file_format = cfg.DATASET.FILE_FORAMT + file_format = cfg.DATASET.FILE_FORMAT download_file_by_url(cfg.DATASET.SOURCE, cfg.DATASET.ROOT, "%s.%s" % (base_dir, file_format), file_format) img_path = os.path.join(cfg.DATASET.ROOT, base_dir, cfg.DATASET.IMG_DIR) diff --git a/examples/cmri_mpca/tutorial.ipynb b/examples/cmri_mpca/tutorial.ipynb index 059ab2a..bd08fa5 100644 --- a/examples/cmri_mpca/tutorial.ipynb +++ b/examples/cmri_mpca/tutorial.ipynb @@ -114,7 +114,7 @@ "metadata": {}, "source": [ "base_dir = cfg.DATASET.BASE_DIR\n", - "file_format = cfg.DATASET.FILE_FORAMT\n", + "file_format = cfg.DATASET.FILE_FORMAT\n", "download_file_by_url(cfg.DATASET.SOURCE, cfg.DATASET.ROOT, \"%s.%s\" % (base_dir, file_format), file_format)" ], "cell_type": "code", diff --git a/examples/multisite_neuroimg_adapt/README.md b/examples/multisite_neuroimg_adapt/README.md index a693a03..37e14c9 100644 --- a/examples/multisite_neuroimg_adapt/README.md +++ b/examples/multisite_neuroimg_adapt/README.md @@ -1,11 +1,11 @@ # Autism Detection: Domain Adaptation for Multi-Site Neuroimaging Data Analysis -### 1. Description +## 1. Description This example demonstrates multi-source domain adaptation method with application in neuroimaging data analysis for autism detection. -### 2. Materials and Methods +## 2. Materials and Methods - Data: Four largest subsets of ABIDE I @@ -15,6 +15,7 @@ autism detection. | UM_1 | 106 | | UCLA_1 | 72 | | USM | 71 | + - Atlas: CC200 - Pre-processing pipeline: cpac - Classification problem: Autism vs Control @@ -22,16 +23,13 @@ autism detection. 1. Constructing brain networks from resting-state fMRI data 2. Classification with Ridge Classifier or Covariate Independence Regularized Least Squares (CoIRLS) classifier +## 3. Related `kalelinear` API -### 3. Related `kale` API - -`kale.interpret.visualize`: Visualize the results of a model. - -`kale.pipeline.multi_domain_adapter.CoIRLS`: Covariate Independence Regularized Least Squares (CoIRLS) classifier. +`kalelinear.estimator.CoIRLS`: Covariate Independence Regularized Least Squares (CoIRLS) classifier. `kale.utils.download.download_file_by_url`: Download a file from a URL. -### References +## References [1] Craddock C., Benhajali Y., Chu C., Chouinard F., Evans A., Jakab A., Khundrakpam BS., Lewis JD., Li Q., Milham M., Yan C. and Bellec P. (2013). [The Neuro Bureau Preprocessing Initiative: Open Sharing of Preprocessed Neuroimaging Data and Derivatives](https://doi.org/10.3389/conf.fninf.2013.09.00041). Frontiers in Neuroinformatics, 7. diff --git a/examples/multisite_neuroimg_adapt/config.py b/examples/multisite_neuroimg_adapt/config.py index 6bcb2a6..3c773a1 100644 --- a/examples/multisite_neuroimg_adapt/config.py +++ b/examples/multisite_neuroimg_adapt/config.py @@ -30,8 +30,8 @@ # ---------------------------------------------------------------------------- # _C.MODEL = CfgNode() _C.MODEL.KERNEL = "rbf" -_C.MODEL.ALPHA = 0.01 -_C.MODEL.LAMBDA_ = 1.0 +_C.MODEL.L2_REG = 0.01 +_C.MODEL.ADAPT_REG = 1.0 # ---------------------------------------------------------------------------- # # Misc options # ---------------------------------------------------------------------------- # diff --git a/examples/multisite_neuroimg_adapt/configs/tutorial.yaml b/examples/multisite_neuroimg_adapt/configs/tutorial.yaml index 0a220b5..1bd397c 100644 --- a/examples/multisite_neuroimg_adapt/configs/tutorial.yaml +++ b/examples/multisite_neuroimg_adapt/configs/tutorial.yaml @@ -1,7 +1,7 @@ MODEL: KERNEL: "linear" - ALPHA: 1.0 - LAMBDA_: 1.0 + L2_REG: 1.0 + ADAPT_REG: 1.0 DATASET: ATLAS: "rois_cc200" SITE_IDS: ['NYU', "UM_1", "UCLA_1", "USM"] diff --git a/examples/multisite_neuroimg_adapt/main.py b/examples/multisite_neuroimg_adapt/main.py index e97b738..9c435de 100644 --- a/examples/multisite_neuroimg_adapt/main.py +++ b/examples/multisite_neuroimg_adapt/main.py @@ -86,7 +86,7 @@ def main(): print(pd.DataFrame.from_dict(results)) print("Domain Adaptation") - estimator = CoIRLS(kernel=cfg.MODEL.KERNEL, lambda_=cfg.MODEL.LAMBDA_, alpha=cfg.MODEL.ALPHA) + estimator = CoIRLS(kernel=cfg.MODEL.KERNEL, lambda_=cfg.MODEL.ADAPT_REG, sigma_=cfg.MODEL.L2_REG) results = cross_validation.leave_one_group_out( brain_networks, pheno["DX_GROUP"].values, pheno["SITE_ID"].values, estimator, use_domain_adaptation=True ) diff --git a/examples/multisite_neuroimg_adapt/tutorial.ipynb b/examples/multisite_neuroimg_adapt/tutorial.ipynb index 3ced871..d0363b6 100644 --- a/examples/multisite_neuroimg_adapt/tutorial.ipynb +++ b/examples/multisite_neuroimg_adapt/tutorial.ipynb @@ -271,7 +271,7 @@ { "metadata": {}, "source": [ - "estimator = CoIRLS(kernel=cfg.MODEL.KERNEL, lambda_=cfg.MODEL.LAMBDA_, alpha=cfg.MODEL.ALPHA)\n", + "estimator = CoIRLS(kernel=cfg.MODEL.KERNEL, lambda_=cfg.MODEL.ADAPT_REG, sigma_=cfg.MODEL.L2_REG)\n", "results = cross_validation.leave_one_group_out(\n", " brain_networks, pheno[\"DX_GROUP\"].values, pheno[\"SITE_ID\"].values, estimator, use_domain_adaptation=True\n", ")" diff --git a/examples/toy_domain_adaptation/README.md b/examples/toy_domain_adaptation/README.md index ba5482b..b7ed09a 100644 --- a/examples/toy_domain_adaptation/README.md +++ b/examples/toy_domain_adaptation/README.md @@ -1,24 +1,21 @@ # Domain Adaptation on Toy Data -### 1. Description +## 1. Description This example demonstrates domain adaptation methods on the generated toy data. -### 2. Usage +## 2. Usage * Datasets: Toy data (2 domains and 2 classes) generated by `sklearn.datasets.make_blobs`. * Algorithms: Covariate Independence Regularized Least Squares (CoIRLS) classifier. `python main.py` -### 3. Related `kale` API +## 3. Related `kalelinear` API -`kale.interpret.visualize.distplot_1d`: Visualize 1D distributions. +`kalelinear.estimator.CoIRLS`: Covariate Independence Regularized Least Squares (CoIRLS) classifier. -`kale.pipeline.multi_domain_adapter.CoIRLS`: Covariate Independence Regularized Least Squares (CoIRLS) classifier. - - -### References +## References [1] Zhou, S., Li, W., Cox, C.R., & Lu, H. (2020). [Side Information Dependence as a Regularizer for Analyzing Human Brain Conditions across Cognitive Experiments](https://ojs.aaai.org//index.php/AAAI/article/view/6179). in *AAAI 2020*, New York, USA. diff --git a/examples/toy_domain_adaptation/main.py b/examples/toy_domain_adaptation/main.py index ba21764..9d7df34 100644 --- a/examples/toy_domain_adaptation/main.py +++ b/examples/toy_domain_adaptation/main.py @@ -81,8 +81,8 @@ def main(): yt_pred_ = clf_.predict(xt) print("Accuracy on target domain: {:.2f}".format(accuracy_score(yt, yt_pred_))) - ys_score_ = clf_.decision_function(xs).detach().numpy().reshape(-1) - yt_score_ = clf_.decision_function(xt).detach().numpy().reshape(-1) + ys_score_ = clf_.decision_function(xs).reshape(-1) + yt_score_ = clf_.decision_function(xt).reshape(-1) title = "Domain adaptation classifier decision score distribution" distplot_1d( [ys_score_, yt_score_], diff --git a/setup.py b/setup.py index a318ede..6a70034 100644 --- a/setup.py +++ b/setup.py @@ -23,6 +23,7 @@ # Dependencies for all examples and tutorials example_requires = [ + "h5py", "ipykernel", "ipython", "matplotlib", From ff413af2d96e85fa5d4ec89fd171292a4dc7277e Mon Sep 17 00:00:00 2001 From: Shuo Zhou Date: Tue, 9 Jun 2026 11:44:53 +0100 Subject: [PATCH 6/6] Fix colab links and install in example notebooks --- examples/brain_lateralization/tutorial.ipynb | 2 +- examples/brain_lateralization/utils/experiment.py | 15 ++++++++------- examples/brain_lateralization/utils/results.py | 12 ++++++------ examples/cmri_mpca/tutorial.ipynb | 9 +++++---- examples/multisite_neuroimg_adapt/README.md | 2 ++ examples/multisite_neuroimg_adapt/tutorial.ipynb | 9 +++++---- examples/toy_domain_adaptation/tutorial.ipynb | 7 ++++--- setup.py | 1 + 8 files changed, 32 insertions(+), 25 deletions(-) diff --git a/examples/brain_lateralization/tutorial.ipynb b/examples/brain_lateralization/tutorial.ipynb index d6c2330..2ab3a26 100644 --- a/examples/brain_lateralization/tutorial.ipynb +++ b/examples/brain_lateralization/tutorial.ipynb @@ -8,7 +8,7 @@ "source": [ "# Group-Specific Discriminant Analysis for sex-specific lateralization Running Demo\n", "\n", - "[Open in Colab](https://colab.research.google.com/github/shuo-zhou/GSDA-Lateralization/blob/main/gsda_demo.ipynb) (click `Runtime`\u2006\u2192\u2006`Run all (Ctrl+F9)`)" + "[Open in Colab](https://colab.research.google.com/github/pykale/kale-linear/blob/main/examples/brain_lateralization/gsda_demo.ipynb) (click `Runtime`\u2006\u2192\u2006`Run all (Ctrl+F9)`)" ], "cell_type": "markdown" }, diff --git a/examples/brain_lateralization/utils/experiment.py b/examples/brain_lateralization/utils/experiment.py index da394d5..356370a 100644 --- a/examples/brain_lateralization/utils/experiment.py +++ b/examples/brain_lateralization/utils/experiment.py @@ -1,13 +1,12 @@ import copy import os +import pickle +from urllib.request import urlretrieve -# import pickle import numpy as np import pandas as pd -import torch from sklearn.metrics import accuracy_score # , roc_auc_score from sklearn.model_selection import StratifiedShuffleSplit -from torch.hub import download_url_to_file from kalelinear.estimator import GSDA # , GSLRTorch @@ -88,7 +87,7 @@ def run_experiment(cfg, lambda_): if download: os.makedirs(data_dir, exist_ok=True) print("Downloading label file for %s dataset. \n" % dataset) - download_url_to_file(LABEL_FILE_LINK[dataset], label_fpath) + urlretrieve(LABEL_FILE_LINK[dataset], label_fpath) else: raise ValueError("File %s does not exist" % label_file) labels = read_tabular(label_fpath, index_col="ID") @@ -178,14 +177,16 @@ def train_model( out_dir, model_filename, ): - model_path = os.path.join(out_dir, "%s.pt" % model_filename) + model_path = os.path.join(out_dir, "%s.pkl" % model_filename) if os.path.exists(model_path): - model = torch.load(model_path) + with open(model_path, "rb") as f: + model = pickle.load(f) else: model = GSDA(**init_kws) model.fit(x_train, **fit_kws) - torch.save(model, model_path) + with open(model_path, "wb") as f: + pickle.dump(model, f) print("Saving trained model to %s. \n" % model_path) return model diff --git a/examples/brain_lateralization/utils/results.py b/examples/brain_lateralization/utils/results.py index 0cc2d68..79f6041 100644 --- a/examples/brain_lateralization/utils/results.py +++ b/examples/brain_lateralization/utils/results.py @@ -1,4 +1,5 @@ import os +import pickle import numpy as np import pandas as pd @@ -19,7 +20,7 @@ def fetch_weights(base_dir, group, lambda_, dataset, sessions, test_size="00", s for half_i in halfs: for i_split in range(num_repeat): for seed in range(52): - model_file_name = "%s_L%s_test_size%s_%s%s_%s_group_%s_%s.pt" % ( + model_file_name = "%s_L%s_test_size%s_%s%s_%s_group_%s_%s.pkl" % ( dataset, lambda_, test_size, @@ -29,18 +30,17 @@ def fetch_weights(base_dir, group, lambda_, dataset, sessions, test_size="00", s group, seed_ - seed, ) - # if os.path.exists(os.path.join(sub_dir, model_file_name)): - # weight.append(get_coef(model_file_name, sub_dir).reshape((1, -1))) + if os.path.exists(os.path.join(sub_dir, model_file_name)): + weight.append(get_coef(model_file_name, sub_dir).reshape((1, -1))) if not weight: return None return np.concatenate(weight, axis=0) def get_coef(file_name, file_dir): - import torch - file_path = os.path.join(file_dir, file_name) - model = torch.load(file_path) + with open(file_path, "rb") as f: + model = pickle.load(f) return model.theta diff --git a/examples/cmri_mpca/tutorial.ipynb b/examples/cmri_mpca/tutorial.ipynb index bd08fa5..91f280d 100644 --- a/examples/cmri_mpca/tutorial.ipynb +++ b/examples/cmri_mpca/tutorial.ipynb @@ -7,7 +7,7 @@ "metadata": {}, "source": [ "# PyKale Tutorial: PAH Diagnosis from Cardiac MRI (CMR) via a Multilinear PCA-based Pipeline\n", - "| [Open in Colab](https://colab.research.google.com/github/pykale/pykale/blob/main/examples/cmri_mpca/tutorial.ipynb) (click `Runtime`\u2006\u2192\u2006`Run all (Ctrl+F9)` | [Launch Binder](https://mybinder.org/v2/gh/pykale/pykale/HEAD?filepath=examples%2Fcmri_mpca%2Ftutorial.ipynb) (click `Run`\u2006\u2192\u2006`Run All Cells`) |" + "| [Open in Colab](https://colab.research.google.com/github/pykale/kale-linear/blob/main/examples/cmri_mpca/tutorial.ipynb) (click `Runtime`\u2006\u2192\u2006`Run all (Ctrl+F9)` | [Launch Binder](https://mybinder.org/v2/gh/pykale/kale-linear/HEAD?filepath=examples%2Fcmri_mpca%2Ftutorial.ipynb) (click `Run`\u2006\u2192\u2006`Run All Cells`) |" ], "cell_type": "markdown" }, @@ -45,9 +45,10 @@ "if 'google.colab' in str(get_ipython()):\n", " print('Running on CoLab')\n", " !pip uninstall --yes imgaug && pip uninstall --yes albumentations && pip install git+https://github.com/aleju/imgaug.git\n", - " !git clone https://github.com/pykale/pykale.git\n", - " %cd pykale\n", - " !pip install .[image,example] \n", + " !pip install pykale\n", + " !git clone https://github.com/pykale/kale-linear.git\n", + " %cd kale-linear\n", + " !pip install .[full]\n", " %cd examples/cmri_mpca\n", "else:\n", " print('Not running on CoLab')" diff --git a/examples/multisite_neuroimg_adapt/README.md b/examples/multisite_neuroimg_adapt/README.md index 37e14c9..365ad28 100644 --- a/examples/multisite_neuroimg_adapt/README.md +++ b/examples/multisite_neuroimg_adapt/README.md @@ -22,6 +22,8 @@ autism detection. - Pipeline: 1. Constructing brain networks from resting-state fMRI data 2. Classification with Ridge Classifier or Covariate Independence Regularized Least Squares (CoIRLS) classifier +- Run: + `python main.py --cfg configs/tutorial.yaml` ## 3. Related `kalelinear` API diff --git a/examples/multisite_neuroimg_adapt/tutorial.ipynb b/examples/multisite_neuroimg_adapt/tutorial.ipynb index d0363b6..fb9a51a 100644 --- a/examples/multisite_neuroimg_adapt/tutorial.ipynb +++ b/examples/multisite_neuroimg_adapt/tutorial.ipynb @@ -7,7 +7,7 @@ "metadata": {}, "source": [ "# PyKale Tutorial: Domain Adaptation for Autism Detection with Multi-site Brain Imaging Data\n", - "| [Open in Colab](https://colab.research.google.com/github/pykale/pykale/blob/main/examples/multisite_neuroimg_adapt/tutorial.ipynb) (click `Runtime`\u2006\u2192\u2006`Run all (Ctrl+F9)` | [Launch Binder](https://mybinder.org/v2/gh/pykale/pykale/HEAD?filepath=examples%2Fmultisite_neuroimg_adapt%2Ftutorial.ipynb) (click `Run`\u2006\u2192\u2006`Run All Cells`) |" + "| [Open in Colab](https://colab.research.google.com/github/pykale/kale-linear/blob/main/examples/multisite_neuroimg_adapt/tutorial.ipynb) (click `Runtime`\u2006\u2192\u2006`Run all (Ctrl+F9)` | [Launch Binder](https://mybinder.org/v2/gh/pykale/kale-linear/HEAD?filepath=examples%2Fmultisite_neuroimg_adapt%2Ftutorial.ipynb) (click `Run`\u2006\u2192\u2006`Run All Cells`) |" ], "cell_type": "markdown" }, @@ -48,9 +48,10 @@ "if 'google.colab' in str(get_ipython()):\n", " print('Running on CoLab')\n", " !pip uninstall --yes imgaug && pip uninstall --yes albumentations && pip install git+https://github.com/aleju/imgaug.git\n", - " !git clone https://github.com/pykale/pykale.git\n", - " %cd pykale\n", - " !pip install .[image,example]\n", + " !pip install pykale\n", + " !git clone https://github.com/pykale/kale-linear.git\n", + " %cd kale-linear\n", + " !pip install .[full]\n", " %cd examples/multisite_neuroimg_adapt\n", "else:\n", " print('Not running on CoLab')" diff --git a/examples/toy_domain_adaptation/tutorial.ipynb b/examples/toy_domain_adaptation/tutorial.ipynb index dedad53..e289eb0 100644 --- a/examples/toy_domain_adaptation/tutorial.ipynb +++ b/examples/toy_domain_adaptation/tutorial.ipynb @@ -7,7 +7,7 @@ "metadata": {}, "source": [ "# PyKale Tutorial: Domain Adaptation on Toy Data\n", - "| [Open in Colab](https://colab.research.google.com/github/pykale/pykale/blob/main/examples/toy_domain_adaptation/tutorial.ipynb) (click `Runtime`\u2006\u2192\u2006`Run all (Ctrl+F9)` | [Launch Binder](https://mybinder.org/v2/gh/pykale/pykale/HEAD?filepath=examples%2Ftoy_domain_adaptation%2Ftutorial.ipynb) (click `Run`\u2006\u2192\u2006`Run All Cells`) |" + "| [Open in Colab](https://colab.research.google.com/github/pykale/kale-linear/blob/main/examples/toy_domain_adaptation/tutorial.ipynb) (click `Runtime`\u2006\u2192\u2006`Run all (Ctrl+F9)` | [Launch Binder](https://mybinder.org/v2/gh/pykale/kale-linear/HEAD?filepath=examples%2Ftoy_domain_adaptation%2Ftutorial.ipynb) (click `Run`\u2006\u2192\u2006`Run All Cells`) |" ], "cell_type": "markdown" }, @@ -36,8 +36,9 @@ " !pip uninstall --yes imgaug && pip uninstall --yes albumentations && pip install git+https://github.com/aleju/imgaug.git\n", " !pip install numpy>=2.0.0\n", " !pip install git+https://github.com/pykale/pykale.git\n", - " !git clone https://github.com/pykale/pykale.git\n", - " %cd pykale/examples/toy_domain_adaptation\n", + " !pip install git+https://github.com/pykale/kale-linear.git\n", + " !git clone https://github.com/pykale/kale-linear.git\n", + " %cd kalelinear/examples/toy_domain_adaptation\n", "else:\n", " print('Not running on CoLab')" ], diff --git a/setup.py b/setup.py index 6a70034..fe51947 100644 --- a/setup.py +++ b/setup.py @@ -27,6 +27,7 @@ "ipykernel", "ipython", "matplotlib", + "nibabel", "nilearn", "pykale", "seaborn",