Full-TL-Example.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"\"\"\n",
    "This example compares the classification performance of \n",
    "linear support vector machine (LinearSVC) on the\n",
    "Riemannian Transfer Learning (RPA, Rodrigues et al., 2018) method\n",
    "and the golden-standard subject-wise train-test cross-validation method\n",
    "using real P300 BCI data.\n",
    "\n",
    "Copyright © 2023, Fahim Doumi <fahim.doumi@outlook.fr> and Fatih Altindis <fthaltindis@gmail.com>\n",
    "Team ViBS (head: Marco Congedo <marco.congedo@gipsa-lab.grenoble-inp.fr>)\n",
    "GIPSA-lab, CNRS, Université Grenoble Alpes\n",
    "License: BSD 3-Clause \n",
    "\n",
    "References:\n",
    "P.L.C. Rodrigues, C. Jutten, M. Congedo (2018)\n",
    "Riemannian procrustes analysis: transfer learning for brain–computer interfaces\n",
    "IEEE Transactions on Biomedical Engineering, 66, 8, 2390-2401.\n",
    "pdf: https://hal.science/hal-01971856/document\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "is_executing": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import warnings \n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from datetime import datetime\n",
    "from joblib import dump, load\n",
    "from tqdm import tqdm\n",
    "\n",
    "from moabb.datasets import BNCI2014008, bi2013a, bi2014a, bi2014b, bi2015a\n",
    "from moabb.paradigms import P300\n",
    "\n",
    "from pyriemann.estimation import ERPCovariances\n",
    "from pyriemann.tangentspace import TangentSpace\n",
    "from pyriemann.transfer import TLSplitter, TLCenter, TLStretch, TLRotate, TLClassifier\n",
    "from pyriemann.transfer import encode_domains, decode_domains\n",
    "\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.model_selection import StratifiedShuffleSplit\n",
    "from sklearn.svm import LinearSVC\n",
    "from sklearn.metrics import balanced_accuracy_score\n",
    "from sklearn.metrics import roc_auc_score\n",
    "from sklearn.exceptions import ConvergenceWarning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# In this program we want to consider 1 session of 1 subject as 1 target\n",
    "def get_subject_data(subject, session, X, y, metadata):\n",
    "    # Create session string based on input session number\n",
    "    session_str = f'session_{session}'\n",
    "\n",
    "    # Select data for a specific subject and session\n",
    "    X_subject_session = X[(metadata[\"subject\"] == subject) & (metadata[\"session\"] == session_str)]\n",
    "    y_subject_session = y[(metadata[\"subject\"] == subject) & (metadata[\"session\"] == session_str)]\n",
    "\n",
    "    return X_subject_session, y_subject_session\n",
    "\n",
    "warnings.filterwarnings('ignore', category=UserWarning, message='Convergence not reached.')\n",
    "warnings.filterwarnings(\"ignore\", category=ConvergenceWarning)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the database, choose the one you want to test\n",
    "dataset = BNCI2014008()\n",
    "paradigm = P300()\n",
    "X, y, metadata = paradigm.get_data(dataset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Selection of the source subject/session :\n",
    "# need to be changed for each different database ;\n",
    "# corresponds to the subject with the best score obtained \n",
    "# with the WithinSessionEvaluation function (from PyRiemann)\n",
    "# using the same pipeline you will use here. \n",
    "# Example in the repository : FindSource.ipynb \n",
    "source = 8\n",
    "session_source = 0\n",
    "\n",
    "subject_list = np.unique(metadata[\"subject\"]) # all subjects from database you want to select\n",
    "session_list = []\n",
    "target_list = []\n",
    "\n",
    "for subject in subject_list:\n",
    "    # get sessions for each subject\n",
    "    sessions = metadata[metadata[\"subject\"]==subject][\"session\"].unique()\n",
    "    # Convert each session to an integer\n",
    "    sessions = [int(s.split('_')[-1]) for s in sessions]\n",
    "    # Add the pairs (subject, session) to target_list\n",
    "    target_list.extend([(subject, session) for session in sessions if not (subject == source and session == session_source)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# number of trials of the target domain for training\n",
    "n_trials = [6, 12, 32, 48]\n",
    "\n",
    "# defining the source and target domain\n",
    "source_domain = f'subject_{source:02}_session_{session_source}'\n",
    "target_domain = ''\n",
    "\n",
    "# Object for splitting the datasets into training and validation partitions\n",
    "# the training set is composed of all data points from the source domain for the RPA\n",
    "# and only partition of the target domain will be training part for the calibration pipeline\n",
    "n_splits = 5 # how many times to split the target domain into train/test for cross-validation\n",
    "seed = 50 # set seed for reproducible results\n",
    "tl_cv = TLSplitter(\n",
    "    target_domain=target_domain,\n",
    "    cv=StratifiedShuffleSplit(n_splits=n_splits, random_state=seed),\n",
    ")\n",
    "\n",
    "# setting up base classifier for TL\n",
    "clf_base = LinearSVC(tol=1e-6, class_weight=\"balanced\")\n",
    "\n",
    "#scores (here 2 types, as many as you want)\n",
    "cumulative_scores_bac = {target: {meth: {trials: [] for trials in n_trials} for meth in ['rpa', 'calibration']} for target in target_list}\n",
    "cumulative_scores_roc = {target: {meth: {trials: [] for trials in n_trials} for meth in ['rpa', 'calibration']} for target in target_list}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load datas for the source subject by using `get_subject_data`\n",
    "d_list_source = []\n",
    "X_source, y_source = get_subject_data(source, session_source, X, y, metadata)\n",
    "erpcov_source = ERPCovariances(classes=[\"Target\"], estimator='lwf') #use SVD if N >= 32\n",
    "cov_source = erpcov_source.fit_transform(X_source, y_source)\n",
    "d_list_source = d_list_source + [f'subject_{source:02}_session_{session_source}'] * len(X_source)\n",
    "d_list_source = np.array(d_list_source)\n",
    "\n",
    "#Encoding the source for TL\n",
    "cov_source_enc, y_source_enc = encode_domains(cov_source, y_source, d_list_source)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#TL LOOP\n",
    "for target in target_list:\n",
    "    subject, session = target\n",
    "    print(\"Subject:\", subject, \"Session:\", session)    \n",
    "    \n",
    "    # Load data for subject using `get_subject_data`\n",
    "    d_list = []\n",
    "    X_target, y_target = get_subject_data(subject, session, X, y, metadata)\n",
    "    d_list = d_list + [f'subject_{subject:02}_session_{session}'] * len(X_target)\n",
    "    domains = np.array(d_list)\n",
    "\n",
    "    # Encoding datas for transfer learning\n",
    "    X_enc, y_enc = encode_domains(X_target, y_target, domains)\n",
    "    \n",
    "    for trials in tqdm(n_trials):\n",
    "        # Create dict for storing results of this particular CV split by scorer\n",
    "        scores_cv_bac = {meth: [] for meth in ['rpa', 'calibration']}\n",
    "        scores_cv_roc = {meth: [] for meth in ['rpa', 'calibration']}\n",
    "\n",
    "        # Change the target domain\n",
    "        target_domain = f'subject_{subject:02}_session_{session}'\n",
    "        tl_cv.target_domain = target_domain\n",
    "\n",
    "        # Change fraction of the target training partition\n",
    "        tl_cv.cv.train_size = trials\n",
    "        print(f\"Number of trials from target domain (Subject: {subject}, Session: {session}) for training is {trials}\")\n",
    "\n",
    "        for train_idx, test_idx in tl_cv.split(X_enc, y_enc):\n",
    "            # Split the target domain into training and testing\n",
    "            X_enc_train, X_enc_test = X_enc[train_idx], X_enc[test_idx]\n",
    "            y_enc_train, y_enc_test = y_enc[train_idx], y_enc[test_idx]\n",
    "            \n",
    "            # Fit ERPCovariances with training trials and obtain covariances matrices\n",
    "            # A different prototype of super trial is obtained with training partition\n",
    "            # of target for each different split and n_trials\n",
    "            # and X_test is transform with this prototype\n",
    "            erpcov_target = ERPCovariances(classes=[\"Target\"], estimator='lwf') #use SVD if N >= 32\n",
    "            cov_train_enc = erpcov_target.fit_transform(X_enc_train, y_target[train_idx])\n",
    "            X_test = erpcov_target.transform(X_enc_test)\n",
    "\n",
    "            # Concatenate training from source and target\n",
    "            X_train = np.concatenate((cov_source_enc, cov_train_enc))\n",
    "            y_train = np.concatenate((y_source_enc, y_enc_train))\n",
    "            \n",
    "            # (1) RPA pipeline: recenter, stretch, and rotate\n",
    "            # Classifier is trained with points from source only\n",
    "            pipeline_rpa = make_pipeline(\n",
    "                TLCenter(target_domain=target_domain),\n",
    "                TLStretch(\n",
    "                    target_domain=target_domain,\n",
    "                    final_dispersion=1,\n",
    "                    centered_data=True,\n",
    "                ),\n",
    "                TLRotate(target_domain=target_domain, metric='riemann'),\n",
    "                TangentSpace(metric=\"riemann\"),\n",
    "                TLClassifier(\n",
    "                    target_domain=target_domain,\n",
    "                    estimator=clf_base,\n",
    "                    domain_weight={source_domain: 1.0, target_domain: 0.0},\n",
    "                ),\n",
    "            )\n",
    "\n",
    "            pipeline_rpa.fit(X_train, y_train)\n",
    "            _, y_true, _ = decode_domains(X_enc_test, y_enc_test)\n",
    "            y_pred_bac_rpa = pipeline_rpa.predict(X_test)\n",
    "            y_test = np.array([y_true == i for i in np.unique(y_true)]).T\n",
    "            y_pred = np.array([y_pred_bac_rpa == i for i in np.unique(y_pred_bac_rpa)]).T\n",
    "            scores_cv_bac['rpa'].append(balanced_accuracy_score(y_true, y_pred_bac_rpa))\n",
    "            scores_cv_roc['rpa'].append(roc_auc_score(y_test, y_pred))\n",
    "\n",
    "            # (2) Calibration: use only data from target-train partition.\n",
    "            # Classifier is trained only with points from the target domain.\n",
    "            pipeline_cal = make_pipeline(\n",
    "                TangentSpace(metric=\"riemann\"),\n",
    "                TLClassifier(\n",
    "                    target_domain=target_domain,\n",
    "                    estimator=clf_base,\n",
    "                    domain_weight={source_domain: 0.0, target_domain: 1.0},\n",
    "                ),\n",
    "            )\n",
    "\n",
    "            pipeline_cal.fit(cov_train_enc, y_enc_train)\n",
    "            _, y_true, _ = decode_domains(X_enc_test, y_enc_test)\n",
    "            y_pred_bac_cal = pipeline_cal.predict(X_test)\n",
    "            y_test = np.array([y_true == i for i in np.unique(y_true)]).T\n",
    "            y_pred = np.array([y_pred_bac_cal == i for i in np.unique(y_pred_bac_cal)]).T\n",
    "            scores_cv_bac['calibration'].append(balanced_accuracy_score(y_true, y_pred_bac_cal))\n",
    "            scores_cv_roc['calibration'].append(roc_auc_score(y_test, y_pred))\n",
    "\n",
    "            # Get the average score of each pipeline\n",
    "        for meth in ['rpa', 'calibration']:\n",
    "            cumulative_scores_bac[target][meth][trials].append(np.mean(scores_cv_bac[meth]))\n",
    "            cumulative_scores_roc[target][meth][trials].append(np.mean(scores_cv_roc[meth]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate average scores for each method and each n_trials per target\n",
    "average_scores_roc = {target: {meth: {trials: np.mean(cumulative_scores_roc[target][meth][trials]) for trials in n_trials} for meth in ['rpa', 'calibration']} for target in target_list}\n",
    "average_scores_bac = {target: {meth: {trials: np.mean(cumulative_scores_bac[target][meth][trials]) for trials in n_trials} for meth in ['rpa', 'calibration']} for target in target_list}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# averages score for all targets\n",
    "average_scores_all_targets_bac = {meth: {trials: 0 for trials in n_trials} for meth in ['rpa', 'calibration']}\n",
    "for meth in ['rpa', 'calibration']:\n",
    "    for trials in n_trials:\n",
    "        all_targets_scores_bac = []\n",
    "        for target in target_list:\n",
    "            all_targets_scores_bac.append(average_scores_bac[target][meth][trials])\n",
    "        average_scores_all_targets_bac[meth][trials] = np.mean(all_targets_scores_bac)\n",
    "\n",
    "average_scores_all_targets_roc = {meth: {trials: 0 for trials in n_trials} for meth in ['rpa', 'calibration']}\n",
    "for meth in ['rpa', 'calibration']:\n",
    "    for trials in n_trials:\n",
    "        all_targets_scores_roc = []\n",
    "        for target in target_list:\n",
    "            all_targets_scores_roc.append(average_scores_roc[target][meth][trials])\n",
    "        average_scores_all_targets_roc[meth][trials] = np.mean(all_targets_scores_roc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Saving the scores \n",
    "path = '/../../'  #your file path\n",
    "database_name = dataset.__class__.__name__\n",
    "timestamp = datetime.now().strftime(\"%Y%m%d_%H%M%S\")\n",
    "dump(average_scores_roc, f'{path}{database_name}_average_scores_roc_{timestamp}.joblib')\n",
    "dump(average_scores_bac, f'{path}{database_name}_average_scores_bac_{timestamp}.joblib')\n",
    "dump(average_scores_all_targets_bac, f'{path}{database_name}_average_scores_all_targets_bac_{timestamp}.joblib')\n",
    "dump(average_scores_all_targets_roc, f'{path}{database_name}_average_scores_all_targets_roc_{timestamp}.joblib')\n",
    "\n",
    "#if i want to load scores\n",
    "#average_scores_bac = load('/../../\"database_name\"_average_scores_all_targets_bac_\"timestamp\".joblib')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plots and tables : changes scores variables if you are using more than 1 scorer "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot for all targets mean score per n_trials\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 5))\n",
    "# for each method\n",
    "for meth in ['rpa', 'calibration']:\n",
    "    # get scores for all trial numbers\n",
    "    scores = [average_scores_all_targets_bac[meth][trials] for trials in n_trials]\n",
    "    \n",
    "    # plot scores\n",
    "    ax.plot(range(len(n_trials)), scores, label=meth, lw=3.0)\n",
    "\n",
    "# set title, labels, etc.\n",
    "ax.set_title(f\"Results for {database_name}\")\n",
    "ax.set_xlabel('Number of training trials in target domain')\n",
    "ax.set_ylabel('Classification score')\n",
    "ax.set_ylim(0.48, 0.6) # change values if needed\n",
    "ax.legend(loc='lower right')\n",
    "\n",
    "# set x-axis ticks and labels\n",
    "ax.set_xticks(range(len(n_trials)))\n",
    "ax.set_xticklabels(n_trials)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot for each n_trials : highly recommanded if you are using all subjects\n",
    "\n",
    "# get a sorted list of subjects\n",
    "subjects_sorted = sorted(target_list)\n",
    "\n",
    "# one plot for each n_trials\n",
    "fig, axs = plt.subplots(len(n_trials), 1, figsize=(10, 5 * len(n_trials)))\n",
    "fig.suptitle(f\"Results for {database_name}\", fontsize=16, y=1.0)\n",
    "\n",
    "# for each ntrials\n",
    "for i, trials in enumerate(n_trials):\n",
    "    ax = axs[i]\n",
    "\n",
    "    # Collect scores and subjects in pairs for each method\n",
    "    scores_subjects = {}\n",
    "    for meth in ['rpa', 'calibration']:\n",
    "        # Pair scores and subjects together\n",
    "        scores_subjects[meth] = [(average_scores_bac[subj][meth][trials], subj) for subj in subjects_sorted]\n",
    "\n",
    "        # Sort pairs by scores in ascending order\n",
    "        scores_subjects[meth].sort()\n",
    "\n",
    "    # Unpack scores and subjects from sorted pairs for each method\n",
    "    for meth in ['rpa', 'calibration']:\n",
    "        scores, sorted_subjects = zip(*scores_subjects[meth])\n",
    "\n",
    "        # plot scores\n",
    "        ax.plot(range(len(scores)), scores, label=meth, lw=3.0)\n",
    "\n",
    "    # set title, labels, etc.\n",
    "    ax.set_title(f\"Training trials in target domain: {int(trials)}\")\n",
    "    ax.set_xlabel('Subject, Session')\n",
    "    ax.set_ylabel('Classification score')\n",
    "    ax.legend(loc='lower right')\n",
    "    ax.set_ylim(0.48, 0.6) # change values if needed\n",
    "    ax.set_xticks(range(len(sorted_subjects)))\n",
    "    ax.set_xticklabels(sorted_subjects)\n",
    "\n",
    "# adjust layout\n",
    "plt.tight_layout()\n",
    "plt.subplots_adjust(right=1.5) # change values if needed\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot for each subject/session, not recommanded if you use all subjects from database\n",
    "# create a figure with one subplot for each subject\n",
    "fig, axs = plt.subplots(len(subjects_sorted), 1, figsize=(10, 5 * len(subjects_sorted)))\n",
    "fig.suptitle(f\"Results for {database_name}\", fontsize=16, y=1.0)\n",
    "\n",
    "# for each subject\n",
    "for i, subj in enumerate(subjects_sorted):\n",
    "    ax = axs[i]\n",
    "    \n",
    "    # for each method\n",
    "    for meth in ['rpa', 'calibration']:\n",
    "        # get scores for all trial numbers\n",
    "        scores = [average_scores_bac[subj][meth][trials] for trials in n_trials]\n",
    "        \n",
    "        # plot scores\n",
    "        ax.plot(range(len(n_trials)), scores, label=meth, lw=3.0)\n",
    "\n",
    "    # set title, labels, etc.\n",
    "    ax.set_title(f\"Subject {subj}\")\n",
    "    ax.set_xlabel('Number of training trials in target domain')\n",
    "    ax.set_ylabel('Classification score')\n",
    "    ax.legend(loc='lower right')\n",
    "    ax.set_ylim(0.48, 0.6) # change values if needed\n",
    "\n",
    "    # set x-axis ticks and labels\n",
    "    ax.set_xticks(range(len(n_trials)))\n",
    "    ax.set_xticklabels(n_trials)\n",
    "\n",
    "# adjust layout\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Table for each n_trials with all targets\n",
    "for trials in n_trials:\n",
    "    data = {}\n",
    "    for meth in ['rpa', 'calibration']:\n",
    "        scores = [average_scores_bac[targets][meth][trials] for subj in subjects_sorted]\n",
    "        data[meth] = scores\n",
    "    df = pd.DataFrame(data, index=subjects_sorted)\n",
    "    print(f\"Number of trials: {trials}\")\n",
    "    print(df)\n",
    "    print(\"\\n---\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Table for all targets mean score per n_trials\n",
    "# Creat dict to store scores for each method \n",
    "data = {}\n",
    "\n",
    "# for each method\n",
    "for meth in ['rpa', 'calibration']:\n",
    "    # Obtenir les scores moyens pour tous les nombres d'essais\n",
    "    scores = [average_scores_all_targets_bac[meth][trials] for trials in n_trials]\n",
    "    # Ajouter les scores à la méthode correspondante dans le dictionnaire\n",
    "    data[meth] = scores\n",
    "\n",
    "# Créer un DataFrame à partir du dictionnaire\n",
    "df_all_subjects = pd.DataFrame(data, index=n_trials)\n",
    "\n",
    "# Renommer l'index\n",
    "df_all_subjects.index.name = 'Number of trials'\n",
    "\n",
    "# Afficher le DataFrame\n",
    "print(df_all_subjects)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.3"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}