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..45d00f5 --- /dev/null +++ b/examples/brain_lateralization/configs/test.yaml @@ -0,0 +1,7 @@ +DATASET: + ROOT: "data" + NUM_REPEAT: 5 +SOLVER: + SEED: 2022 +OUTPUT: + ROOT: 'output' diff --git a/examples/brain_lateralization/tutorial.ipynb b/examples/brain_lateralization/tutorial.ipynb new file mode 100644 index 0000000..2ab3a26 --- /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/pykale/kale-linear/blob/main/examples/brain_lateralization/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/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" + }, + { + "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/data_io.py b/examples/brain_lateralization/utils/data_io.py new file mode 100644 index 0000000..9bba2ca --- /dev/null +++ b/examples/brain_lateralization/utils/data_io.py @@ -0,0 +1,174 @@ +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() + + +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 new file mode 100644 index 0000000..356370a --- /dev/null +++ b/examples/brain_lateralization/utils/experiment.py @@ -0,0 +1,672 @@ +import copy +import os +import pickle +from urllib.request import urlretrieve + +import numpy as np +import pandas as pd +from sklearn.metrics import accuracy_score # , roc_auc_score +from sklearn.model_selection import StratifiedShuffleSplit + +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": [], + "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"] + +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_): + 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) + 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") + 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_no_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_model( + x_train, + init_kws, + fit_kws, + out_dir, + model_filename, +): + model_path = os.path.join(out_dir, "%s.pkl" % model_filename) + + if os.path.exists(model_path): + with open(model_path, "rb") as f: + model = pickle.load(f) + else: + model = GSDA(**init_kws) + model.fit(x_train, **fit_kws) + with open(model_path, "wb") as f: + pickle.dump(model, f) + 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 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": [ + 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 = 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, + i_split, + i_fold, + tgt_group, + random_state, + ) + if mix_group: + model_filename = model_filename + "_group_mix" + 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, + 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 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, + test_sub_tgt_idx, + test_sub_nt_idx, + ) = _idx + + 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 = 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)), + 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 = 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, + 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 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, + test_sub_tgt_idx, + test_sub_nt_idx, + ) = _idx + + 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 = 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)), + 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 = 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 = { + "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 = get_fit_kws(y_train_fold[train_sub], groups[train_sub]) + else: + 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, + 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 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], 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 = 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, + i_fold, + tgt_group, + random_state, + ) + + if mix_group: + model_filename = model_filename + "_group_mix" + fit_kws = get_fit_kws(y_train, groups) + init_kws = get_gsda_init_kws(lambda_, l2_hparam, kwargs) + + model = train_model( + 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/half_brain.py b/examples/brain_lateralization/utils/half_brain.py new file mode 100644 index 0000000..2742662 --- /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 new file mode 100644 index 0000000..87a100e --- /dev/null +++ b/examples/brain_lateralization/utils/io_.py @@ -0,0 +1,75 @@ +"""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_result, + load_txt, + read_file, + read_tabular, + reformat_results, + 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_result", + "load_txt", + "pick_half", + "read_file", + "read_nii", + "read_shape_gii", + "read_surf", + "read_tabular", + "reformat_results", + "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 new file mode 100644 index 0000000..d3c5931 --- /dev/null +++ b/examples/brain_lateralization/utils/plot.py @@ -0,0 +1,8 @@ +"""Compatibility wrapper for plotting utilities. + +Prefer importing from ``utils.plotting`` in new code. +""" + +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", "plot_accuracy", "plot_gsi"] diff --git a/examples/brain_lateralization/utils/plotting.py b/examples/brain_lateralization/utils/plotting.py new file mode 100644 index 0000000..8a5001e --- /dev/null +++ b/examples/brain_lateralization/utils/plotting.py @@ -0,0 +1,207 @@ +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() + + +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 new file mode 100644 index 0000000..79f6041 --- /dev/null +++ b/examples/brain_lateralization/utils/results.py @@ -0,0 +1,53 @@ +import os +import pickle + +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.pkl" % ( + 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))) + if not weight: + return None + return np.concatenate(weight, axis=0) + + +def get_coef(file_name, file_dir): + file_path = os.path.join(file_dir, file_name) + with open(file_path, "rb") as f: + model = pickle.load(f) + 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) 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..e44f60c --- /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_FORMAT = "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..ccfb3be --- /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_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) + 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..91f280d --- /dev/null +++ b/examples/cmri_mpca/tutorial.ipynb @@ -0,0 +1,404 @@ +{ + "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/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" + }, + { + "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", + " !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')" + ], + "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_FORMAT\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..365ad28 --- /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 +- Run: + `python main.py --cfg configs/tutorial.yaml` + +## 3. Related `kalelinear` API + +`kalelinear.estimator.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..3c773a1 --- /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.L2_REG = 0.01 +_C.MODEL.ADAPT_REG = 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..1bd397c --- /dev/null +++ b/examples/multisite_neuroimg_adapt/configs/tutorial.yaml @@ -0,0 +1,7 @@ +MODEL: + KERNEL: "linear" + 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 new file mode 100644 index 0000000..9c435de --- /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.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 + ) + 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..fb9a51a --- /dev/null +++ b/examples/multisite_neuroimg_adapt/tutorial.ipynb @@ -0,0 +1,294 @@ +{ + "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/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" + }, + { + "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", + " !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')" + ], + "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.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", + ")" + ], + "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..b7ed09a --- /dev/null +++ b/examples/toy_domain_adaptation/README.md @@ -0,0 +1,22 @@ +# 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 `kalelinear` API + +`kalelinear.estimator.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..9d7df34 --- /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).reshape(-1) + yt_score_ = clf_.decision_function(xt).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..e289eb0 --- /dev/null +++ b/examples/toy_domain_adaptation/tutorial.ipynb @@ -0,0 +1,219 @@ +{ + "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/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" + }, + { + "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", + " !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')" + ], + "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..fe51947 100644 --- a/setup.py +++ b/setup.py @@ -23,10 +23,13 @@ # Dependencies for all examples and tutorials example_requires = [ + "h5py", "ipykernel", "ipython", "matplotlib", + "nibabel", "nilearn", + "pykale", "seaborn", "yacs>=0.1.7", ]