Visualizations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c97fe7ca",
   "metadata": {},
   "source": [
    "# Notebook to analyze and visualize results\n",
    "\n",
    "Implemented are\n",
    "* plotting of \n",
    "    * receiver-operator-curves (ROC)\n",
    "    * precision-recall-curves (PRC)\n",
    "    * confusion matrices\n",
    "* sensitivity / specificity / negative predictive value analysis\n",
    "* data distributions for properties such as classes, tiles, age\n",
    "* analysis of clinicopathological fetures (via plotting of AUROCs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2305238",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fbe09984",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import argparse\n",
    "from pathlib import Path\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.pylab as pl\n",
    "import pandas as pd\n",
    "from tqdm import tqdm\n",
    "from sklearn.metrics import roc_curve, roc_auc_score, auc, precision_recall_curve, average_precision_score, confusion_matrix\n",
    "import yaml\n",
    "import h5py\n",
    "import torch\n",
    "\n",
    "from options import Options"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ff93680",
   "metadata": {},
   "outputs": [],
   "source": [
    "parser = Options()\n",
    "args = parser.parser.parse_args('')  \n",
    "\n",
    "# Load the configuration from the YAML file\n",
    "with open(args.config_file, 'r') as f:\n",
    "    config = yaml.load(f, Loader=yaml.FullLoader)\n",
    "\n",
    "# Update the configuration with the values from the argument parser\n",
    "for arg_name, arg_value in vars(args).items():\n",
    "    if arg_value is not None and arg_name != 'config_file':\n",
    "        config[arg_name] = getattr(args, arg_name)\n",
    "\n",
    "print('\\n--- load options ---')\n",
    "for name, value in sorted(config.items()):\n",
    "    print(f'{name}: {str(value)}')\n",
    "\n",
    "cfg = argparse.Namespace(**config)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a123f0ec",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "## Load predictions for runs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cdfb3873",
   "metadata": {},
   "outputs": [],
   "source": [
    "cfg.name = 'multi-all-cohorts'\n",
    "# cfg.name = 'single'\n",
    "\n",
    "# cfg.cohorts = ['DACHS', 'QUASAR', 'RAINBOW', 'TCGA', 'MCO']\n",
    "cfg.cohorts = ['CPTAC', 'DACHS', 'DUSSEL', 'Epi700', 'ERLANGEN', 'FOXTROT', 'MCO', 'MECC', 'MUNICH', 'QUASAR', 'RAINBOW', 'TCGA', 'TRANSCOT']\n",
    "# cfg.cohorts = ['FOXTROT']\n",
    "\n",
    "cfg.target = 'isMSIH'\n",
    "# cfg.target = 'BRAF'\n",
    "# cfg.target = 'KRAS'\n",
    "true_label = 1\n",
    "\n",
    "cfg.norm = 'histaugan'\n",
    "\n",
    "cfg.model = 'transformer'\n",
    "\n",
    "cfg.logging_name = f'{cfg.name}_{cfg.model}_{\"-\".join(cfg.cohorts)}_{cfg.norm}_{cfg.target}'\n",
    "# base_path = Path('/Volumes/SSD/logs/idkidc/') / cfg.logging_name\n",
    "base_path = Path(cfg.save_dir) / cfg.logging_name\n",
    "result_path = base_path / 'results'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42d0913d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# --- choose training and test cohort\n",
    "# test_cohorts = 'YCR-BCIP-biopsies'\n",
    "test_cohorts = 'YCR-BCIP-resections'\n",
    "# test_cohorts = 'CHINA'\n",
    "# test_cohorts = 'MAINZ'\n",
    "# test_cohorts = 'Epi700'\n",
    "# test_cohorts = 'FOXTROT'\n",
    "csv_paths = result_path.glob(f'fold*/outputs_{test_cohorts}.csv')\n",
    "pred_csvs = list(csv_paths)\n",
    "\n",
    "# train_cohorts = 'DACHS, QUASAR, RAINBOW, TCGA, MCO'\n",
    "train_cohorts = 'CPTAC, DACHS, DUSSEL, Epi700, ERLANGEN, FOXTROT, MCO, MECC, MUNICH, QUASAR, RAINBOW, TCGA, TRANSCOT'\n",
    "# train_cohorts = 'FOXTROT'\n",
    "csv_paths_indomain = result_path.glob(f'fold*/outputs_{train_cohorts}.csv')\n",
    "pred_csvs_indomain = list(csv_paths_indomain)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a75233c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load predictions for test cohort\n",
    "pred_dfs = [pd.read_csv(p, dtype=str) for p in pred_csvs]\n",
    "\n",
    "y_trues = [df['ground_truth'] == f'{true_label:.1f}' for df in pred_dfs]\n",
    "y_preds = [pd.to_numeric(df[f'logits']) for df in pred_dfs]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62570c47",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load predictions for indomain test set\n",
    "pred_dfs_indomain = [pd.read_csv(p, dtype=str) for p in pred_csvs_indomain]\n",
    "\n",
    "y_trues_indomain = [df['ground_truth'] == f'{true_label:.1f}' for df in pred_dfs_indomain]\n",
    "y_preds_indomain = [pd.to_numeric(df[f'logits']) for df in pred_dfs_indomain]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27ce69b4",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "## Plot ROC curves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2252a9b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# figure_path = Path('/Users/sophia.wagner/Documents/PhD/projects/2022_MSI_transformer/figures') / 'curves'\n",
    "figure_path = Path('figures') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f61d776e",
   "metadata": {},
   "outputs": [],
   "source": [
    "tprs = []\n",
    "aucs = []\n",
    "mean_fpr = np.linspace(0, 1, 100)\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(5, 5))\n",
    "\n",
    "for i in range(len(y_trues)):\n",
    "    fpr, tpr, _ = roc_curve(y_trues[i], y_preds[i])\n",
    "    roc_auc = roc_auc_score(y_trues[i], y_preds[i])\n",
    "    \n",
    "#     ax.plot(fpr, tpr, label=f'AUC = {roc_auc:0.2f}', alpha=0.3, color='black')\n",
    "\n",
    "\n",
    "    interp_tpr = np.interp(mean_fpr, fpr, tpr)\n",
    "    interp_tpr[0] = 0.0\n",
    "    tprs.append(interp_tpr)\n",
    "    aucs.append(roc_auc)\n",
    "print(aucs)\n",
    "ax.plot([0, 1], [0, 1], linestyle=\"--\", lw=1, color=\"gray\", label=\"Chance\", alpha=0.8)\n",
    "\n",
    "mean_tpr = np.mean(tprs, axis=0)\n",
    "mean_tpr[-1] = 1.0\n",
    "mean_auc = auc(mean_fpr, mean_tpr)\n",
    "std_auc = np.std(aucs)\n",
    "# plot mean curve\n",
    "ax.plot(\n",
    "    mean_fpr,\n",
    "    mean_tpr,\n",
    "    color=\"black\",\n",
    "    label=r\"Mean ROC (AUC = %0.2f $\\pm$ %0.3f)\" % (np.round(mean_auc, 2), np.round(std_auc, 3)),\n",
    "    lw=2,\n",
    "    alpha=0.8,\n",
    ")\n",
    "\n",
    "std_tpr = np.std(tprs, axis=0)\n",
    "tprs_upper = np.minimum(mean_tpr + std_tpr, 1)\n",
    "tprs_lower = np.maximum(mean_tpr - std_tpr, 0)\n",
    "ax.fill_between(\n",
    "    mean_fpr,\n",
    "    tprs_lower,\n",
    "    tprs_upper,\n",
    "    color=\"grey\",\n",
    "    alpha=0.5,\n",
    "    label=r\"$\\pm$ 1 std. dev.\",\n",
    ")\n",
    "\n",
    "ax.set(\n",
    "    xlim=[-0.05, 1.05],\n",
    "    ylim=[-0.05, 1.05],\n",
    "    # title=f\"ROC for {test_cohorts} (target: {cfg.target})\",\n",
    ")\n",
    "\n",
    "ax.set_xlabel('1 - Specificity')\n",
    "ax.set_ylabel('Sensitivity')\n",
    "\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "\n",
    "ax.legend(loc=\"lower right\")\n",
    "fig.savefig(figure_path / f'auroc_{cfg.logging_name}_{test_cohorts}.svg', format='svg', bbox_inches = 'tight', pad_inches = 0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06221e2e",
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_auc, std_auc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11183e59",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot zoom in\n",
    "tprs = []\n",
    "aucs = []\n",
    "mean_fpr = np.linspace(0, 1, 100)\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(3,3))\n",
    "\n",
    "for i in range(len(y_trues)):\n",
    "    fpr, tpr, _ = roc_curve(y_trues[i], y_preds[i])\n",
    "    roc_auc = roc_auc_score(y_trues[i], y_preds[i])\n",
    "    \n",
    "    ax.plot(fpr, tpr, lw=1, label=f'AUC = {roc_auc:0.2f}', color=pl.cm.viridis(i*(1/len(y_trues))))\n",
    "\n",
    "    interp_tpr = np.interp(mean_fpr, fpr, tpr)\n",
    "    interp_tpr[0] = 0.0\n",
    "    tprs.append(interp_tpr)\n",
    "    aucs.append(roc_auc)\n",
    "ax.plot([0, 1], [0, 1], linestyle=\"--\", lw=1, color=\"gray\", label=\"Chance\", alpha=0.8)\n",
    "\n",
    "mean_tpr = np.mean(tprs, axis=0)\n",
    "mean_tpr[-1] = 1.0\n",
    "mean_auc = auc(mean_fpr, mean_tpr)\n",
    "std_auc = np.std(aucs)\n",
    "# plot mean curve\n",
    "ax.plot(\n",
    "    mean_fpr,\n",
    "    mean_tpr,\n",
    "    color=\"black\",\n",
    "    label=r\"Mean ROC (AUC = %0.2f $\\pm$ %0.3f)\" % (mean_auc, std_auc),\n",
    "    lw=2,\n",
    "    alpha=0.8,\n",
    ")\n",
    "\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "\n",
    "ax.set(\n",
    "    xlim=[0., .25],\n",
    "    ylim=[0.75, 1.],\n",
    ")\n",
    "\n",
    "# fig.savefig(figure_path / f'auroc_{cfg.logging_name}_{test_cohorts}_zoom.svg',  format='svg', bbox_inches = 'tight', pad_inches = 0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01a11224",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "## Plot PR curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73dd5c03",
   "metadata": {},
   "outputs": [],
   "source": [
    "precs = []\n",
    "aucs = []\n",
    "mean_recall = np.linspace(0, 1, 100)\n",
    "\n",
    "y_trues = [df[\"ground_truth\"] == f'{true_label:.1f}' for df in pred_dfs]\n",
    "y_preds = [pd.to_numeric(df[f'logits']) for df in pred_dfs]\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(5, 5))\n",
    "\n",
    "for i in range(len(y_trues)):\n",
    "    prec, rec, _ = precision_recall_curve(y_trues[i], y_preds[i])\n",
    "    prec, rec = prec[::-1], rec[::-1]  # reverse order for interp\n",
    "    pr_auc = average_precision_score(y_trues[i], y_preds[i])\n",
    "    \n",
    "    ax.plot(rec, prec, label=f'AUC = {pr_auc:0.2f}', color=pl.cm.viridis(i*(1/len(y_trues))))\n",
    "    \n",
    "    baseline = y_trues[i].sum()/len(y_trues[i])\n",
    "\n",
    "    interp_prec = np.interp(mean_recall, rec, prec)\n",
    "    interp_prec[0] = 0.0\n",
    "    precs.append(interp_prec)\n",
    "    aucs.append(pr_auc)\n",
    "    \n",
    "y_trues = np.concatenate(y_trues, dtype=int)\n",
    "y_preds = np.concatenate(y_preds, dtype=float)\n",
    "prec, rec, _ = precision_recall_curve(y_trues, y_preds)\n",
    "\n",
    "pr_auc = average_precision_score(y_trues, y_preds)\n",
    "\n",
    "ax.plot([0, 1], [baseline, baseline], '--', label=\"Chance\", alpha=0.8, color=\"gray\", lw=1)\n",
    "\n",
    "ax.set(\n",
    "    xlim=[-0.05, 1.05],\n",
    "    ylim=[-0.05, 1.05],\n",
    "    # title=f\"PRC for {test_cohorts} (target: {cfg.target})\",\n",
    ")\n",
    "\n",
    "ax.set_xlabel('Recall')\n",
    "ax.set_ylabel('Precision')\n",
    "\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "\n",
    "ax.legend(loc=\"lower left\")\n",
    "fig.savefig(figure_path / f'prauc_{cfg.logging_name}_{test_cohorts}.svg', format='svg', bbox_inches = 'tight', pad_inches = 0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "166173df",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "## Plot confusion matrix, analyze sensitifity and specificity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1b1ff01",
   "metadata": {},
   "outputs": [],
   "source": [
    "desired_sensitivity = 0.95"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc11da2d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# determine classification threshold on test data\n",
    "thresholds_test = []\n",
    "sensitivity = []\n",
    "specificity = []\n",
    "# Calculate the ROC curve\n",
    "for k in range(5):\n",
    "    fpr, tpr, thresholds = roc_curve(y_trues[k], y_preds[k])\n",
    "\n",
    "    # Find the threshold for a desired sensitivity\n",
    "    thresholds_test.append(thresholds[tpr >= desired_sensitivity][0])\n",
    "    sensitivity.append(tpr[tpr >= desired_sensitivity][0])\n",
    "    specificity.append(1-fpr[tpr >= desired_sensitivity][0])\n",
    "#     print(1-fpr[tpr >= 0.95])\n",
    "\n",
    "threshold_test = np.array(thresholds_test).mean()\n",
    "# Print the threshold\n",
    "print(f\"Mean treshold for desired sensitivity of {desired_sensitivity:.2f}: {np.array(thresholds_test).mean():.4f}\")\n",
    "print(f\"The model achieves sensitivity: {np.array(sensitivity).mean():.4f} with specificty {np.array(specificity).mean():.4f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ebb79b30",
   "metadata": {},
   "outputs": [],
   "source": [
    "csv_paths = result_path.glob(f'fold*/outputs_{test_cohorts}.csv')\n",
    "pred_csvs = list(csv_paths)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4c66960",
   "metadata": {},
   "outputs": [],
   "source": [
    "csv_paths = result_path.glob(f'fold*/outputs_{test_cohorts}.csv')\n",
    "pred_csvs = list(csv_paths)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b45fec9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# determine best fold according to AUROC score\n",
    "rocs = []\n",
    "for i in range(len(y_trues)):\n",
    "    rocs.append(roc_auc_score(y_trues[i], y_preds[i])) \n",
    "    # rocs.append(roc_auc_score(y_trues_train[i], y_preds_train[i]))\n",
    "best_fold = np.argmax(rocs)\n",
    "print(f'best fold: {best_fold}, with AUROC {rocs[best_fold]:.4f}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c91f40dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "rocs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3fd61153",
   "metadata": {},
   "outputs": [],
   "source": [
    "desired_sensitivity = 0.9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f745835a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# determine classification threshold on test data\n",
    "thresholds_train = []\n",
    "sensitivity_indomain = []\n",
    "specificity_indomain =[]\n",
    "\n",
    "sensitivity = []\n",
    "specificity =[]\n",
    "# Calculate the ROC curve\n",
    "for k in range(5):\n",
    "    fpr_indomain, tpr_indomain, thresholds_indomain = roc_curve(y_trues_indomain[k], y_preds_indomain[k])\n",
    "    fpr, tpr, thresholds = roc_curve(y_trues[k], y_preds[k])\n",
    "\n",
    "    # Find the threshold for a desired sensitivity\n",
    "    thresholds_train.append(thresholds_indomain[tpr_indomain >= desired_sensitivity][0])\n",
    "    sensitivity_indomain.append(tpr_indomain[tpr_indomain >= desired_sensitivity][0])\n",
    "    specificity_indomain.append(1-fpr_indomain[tpr_indomain >= desired_sensitivity][0])\n",
    "    \n",
    "    sensitivity.append(tpr[thresholds >= thresholds_train[k]][-1])\n",
    "    specificity.append(1-fpr[thresholds >= thresholds_train[k]][-1])\n",
    "    \n",
    "threshold_train = np.array(thresholds_train).mean()\n",
    "# Print the threshold\n",
    "print(thresholds_train)\n",
    "print(f\"Mean treshold for desired sensitivity of {desired_sensitivity:.2f}: {np.array(thresholds_train).mean():.4f}\")\n",
    "print(f\"Sensitivity (with in-domain threshold): {np.array(sensitivity_indomain).mean():.4f} with specificty {np.array(specificity_indomain).mean():.4f}\")\n",
    "print(f\"Sensitivity (with external threshold): {np.array(sensitivity).mean():.4f} with specificty {np.array(specificity).mean():.4f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45f30ad7",
   "metadata": {},
   "outputs": [],
   "source": [
    "npv = []\n",
    "for k in range(5):\n",
    "    tn, fp, fn, tp = confusion_matrix(y_trues[k], y_preds[k] > thresholds_train[k]).flatten()\n",
    "    npv.append(tn / (tn + fn))\n",
    "print(f\"Negative predictive value (with in-domain threshold): {np.array(npv).mean():.4f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54f9eb95",
   "metadata": {},
   "outputs": [],
   "source": [
    "npv = []\n",
    "for k in range(5):\n",
    "    tn, fp, fn, tp = confusion_matrix(y_trues[k], y_preds[k] > thresholds_test[k]).flatten()\n",
    "    npv.append(tn / (tn + fn))\n",
    "print(f\"Negative predictive value (with external threshold): {np.array(npv).mean():.4f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a9b17ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "threshold_train, threshold_test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "deb6e1cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# determine negative predictive value\n",
    "print(f\"Negative predictive value (with in-domain threshold): {np.array(specificity_indomain).mean()/(1-np.array(sensitivity_indomain).mean()):.4f}\")\n",
    "print(f\"Negative predictive value (with external threshold): {np.array(specificity).mean()/(1-np.array(sensitivity).mean()):.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa3062af",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "## determine thresholds for each fold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a557bb21",
   "metadata": {},
   "outputs": [],
   "source": [
    "# determine classification threshold on in-domain test data\n",
    "best_thresholds = []\n",
    "# choose threshold according to desired sensitivity or geometric mean\n",
    "desired_sensitivity = 0.9\n",
    "\n",
    "for k in range(5):\n",
    "    # in-domain test set\n",
    "    fpr, tpr, thresholds = roc_curve(y_trues_indomain[k], y_preds_indomain[k])\n",
    "    \n",
    "    if desired_sensitivity:\n",
    "        best_threshold = thresholds[tpr >= desired_sensitivity][0]\n",
    "    else:\n",
    "        gm = np.sqrt(tpr * (1 - fpr))\n",
    "        best_threshold = thresholds[np.argmax(gm)]\n",
    "\n",
    "    best_thresholds.append(best_threshold)\n",
    "    \n",
    "print(best_thresholds)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3587ca60",
   "metadata": {},
   "outputs": [],
   "source": [
    "sensitivity = []\n",
    "specificity = []\n",
    "precision = []\n",
    "npv = []\n",
    "\n",
    "# threshold = [0.6940] * 5\n",
    "# threshold = [0.5] * 5\n",
    "threshold = best_thresholds\n",
    "for k in range(5):\n",
    "    tn, fp, fn, tp = confusion_matrix(y_trues_indomain[k], y_preds_indomain[k] > threshold[k]).ravel()\n",
    "    sensitivity.append(tp / (tp + fn))  # recall\n",
    "    specificity.append(tn / (tn + fp))\n",
    "    precision.append(tp / (tp + fp))\n",
    "    npv.append(tn / (tn + fn))\n",
    "\n",
    "print(\"In-domain test set\")\n",
    "print(f\"Mean treshold: {np.mean(threshold):.3f} with thresholds in each fold: {threshold}\")\n",
    "print(f\"Sensitivity {np.mean(sensitivity):.3f} ({np.std(sensitivity):.4f}) with specificty {np.mean(specificity):.3f} ({np.std(specificity):.4f})\")\n",
    "print(f\"Negative predictive value: {np.mean(npv):.3f} ({np.std(npv):.4f})\")\n",
    "print(f\"Precision: {np.mean(precision):.3f} ({np.std(precision):.4f})\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9b66791",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load evaluation cohort\n",
    "# test_cohorts = 'CHINA' \n",
    "# test_cohorts = 'Epi700' \n",
    "test_cohorts = 'MAINZ'\n",
    "# test_cohorts = 'YCR-BCIP-resections'\n",
    "# test_cohorts = 'YCR-BCIP-biopsies'\n",
    "\n",
    "csv_paths = result_path.glob(f'fold*/outputs_{test_cohorts}.csv')\n",
    "pred_csvs = list(csv_paths)\n",
    "pred_dfs = [pd.read_csv(p, dtype=str) for p in pred_csvs]\n",
    "\n",
    "y_trues = [df['ground_truth'] == f'{true_label:.1f}' for df in pred_dfs]\n",
    "y_preds = [pd.to_numeric(df[f'logits']) for df in pred_dfs]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2a4f8e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "sensitivity = []\n",
    "specificity = []\n",
    "precision = []\n",
    "npv = []\n",
    "\n",
    "# threshold = [0.6940] * 5\n",
    "# threshold = [0.5] * 5\n",
    "threshold = best_thresholds\n",
    "for k in range(5):\n",
    "    tn, fp, fn, tp = confusion_matrix(y_trues[k], y_preds[k] > threshold[k]).ravel()\n",
    "    sensitivity.append(tp / (tp + fn))  # recall\n",
    "    specificity.append(tn / (tn + fp))\n",
    "    precision.append(tp / (tp + fp))\n",
    "    npv.append(tn / (tn + fn))\n",
    "\n",
    "print(f\"{test_cohorts} test set\")\n",
    "print(f\"Mean treshold: {np.mean(threshold):.3f} with thresholds in each fold: {threshold}\")\n",
    "print(f\"Sensitivity {np.mean(sensitivity):.3f} ({np.std(sensitivity):.4f}) with specificty {np.mean(specificity):.3f} ({np.std(specificity):.4f})\")\n",
    "print(f\"Negative predictive value: {np.mean(npv):.3f} ({np.std(npv):.4f})\")\n",
    "print(f\"Precision: {np.mean(precision):.3f} ({np.std(precision):.4f})\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb1c2034",
   "metadata": {},
   "outputs": [],
   "source": [
    "specificity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89058319",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the class labels\n",
    "labels = [\"non MSIH\", \"MSIH\"]\n",
    "\n",
    "# Calculate the confusion matrix for different classification thresholds\n",
    "# thresholds_test = np.array(thresholds_test).mean().repeat(5)\n",
    "# thresholds_train = np.array(thresholds_train).mean().repeat(5)\n",
    "thresholds = [thresholds_test, [0.25] * 5, [0.5] * 5, [0.75] * 5, thresholds_train]\n",
    "\n",
    "# Plot the confusion matrices\n",
    "fig, axs = plt.subplots(1, len(thresholds), figsize=(12,10), gridspec_kw={'width_ratios': [4, 4, 4, 4, 4.4]})\n",
    "for ax, threshold in zip(axs.flatten(), thresholds):\n",
    "    cm = []\n",
    "    for k in range(5):\n",
    "        cm.append(confusion_matrix(y_trues[k], y_preds[k] > threshold[k]))\n",
    "    cm = np.stack(cm, axis=0).mean(axis=0)    \n",
    "    print(threshold[k])\n",
    "#     cm = confusion_matrix(y_trues[best_fold], y_preds[best_fold] > threshold)\n",
    "\n",
    "    # Normalize the confusion matrices by the row\n",
    "    cm_normalized = cm.astype(\"float\") / cm.sum(axis=1)[:, np.newaxis]\n",
    "    cm_perc = cm.astype(\"float\") / cm.sum(axis=(0, 1))[np.newaxis, np.newaxis]\n",
    "\n",
    "    im = ax.imshow(cm_normalized, interpolation='nearest', cmap=cmap)\n",
    "    ax.set(xticks=np.arange(cm.shape[1]),\n",
    "           yticks=np.arange(cm.shape[0]),\n",
    "           xticklabels=labels, yticklabels=labels,\n",
    "           title=f\"Threshold = {threshold[k]:.2f}\",\n",
    "#            ylabel='True label',\n",
    "           xlabel='Predicted label')\n",
    "    ax.set_yticklabels(labels, rotation='vertical', va='center', ha='center')\n",
    "    for i in range(cm.shape[0]):\n",
    "        for j in range(cm.shape[1]):\n",
    "            color = \"w\" if cm_normalized[i, j] > 0.5 else \"k\"\n",
    "            ax.text(j, i, int(cm[i, j].round()), ha=\"center\", va=\"bottom\", color=color)\n",
    "            ax.text(j, i, f'{cm_normalized[i, j]*100:.01f}%', ha=\"center\", va=\"top\", color=color)\n",
    "#             ax.text(j, i, f'{cm_perc[i, j]*100:.01f}%', ha=\"center\", va=\"top\", color=color)\n",
    "\n",
    "cbar = fig.colorbar(im, ax=axs[-1], fraction=0.046, pad=0.04, ticks=[0.02, 0.2, 0.4, 0.6, 0.8, 0.98])\n",
    "cbar.ax.set_yticklabels(np.linspace(0, 1, 6, dtype=np.float16))  # horizontal colorbar\n",
    "axs[0].set(ylabel='True label')\n",
    "axs[0].set(title='External threshold')\n",
    "axs[-1].set(title='In-domain threshold')\n",
    "# plt.tight_layout()\n",
    "fig.savefig(output_path / f'confusion_matrices_{log_name}_{target_label}_{eval_name}_mean.svg',  format='svg', bbox_inches = 'tight', pad_inches = 0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a0453d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "cm.astype(\"float\") / cm.sum(axis=(0, 1))[np.newaxis, np.newaxis]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c70fdc5",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.array(thresholds_test).mean().repeat(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee1d5fa4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the class labels\n",
    "labels = [\"non MSIH\", \"MSIH\"]\n",
    "\n",
    "# Calculate the confusion matrix for different classification thresholds\n",
    "# threshsolds = [np.array(thresholds_test).mean()]\n",
    "thresholds = [thresholds_test, [0.25] * 5, [0.5] * 5, [0.75] * 5, thresholds_train]\n",
    "\n",
    "# Plot the confusion matrices\n",
    "fig, axs = plt.subplots(5, len(thresholds), figsize=(12,15), gridspec_kw={'width_ratios': [4, 4, 4, 4, 4.4]})\n",
    "for k in range(5):\n",
    "    for t in range(5):\n",
    "        cm = confusion_matrix(y_trues[k], y_preds[k] > thresholds[t][k])\n",
    "        print(t, k, thresholds[t][k], cm.sum(axis=1))\n",
    "\n",
    "        # Normalize the confusion matrices by the row\n",
    "        cm_normalized = cm.astype(\"float\") / cm.sum(axis=1)[:, np.newaxis]\n",
    "\n",
    "        im = axs[k, t].imshow(cm_normalized, interpolation='nearest', cmap=cmap)\n",
    "        axs[k, t].set(xticks=np.arange(cm.shape[1]),\n",
    "               yticks=np.arange(cm.shape[0]),\n",
    "               xticklabels=labels, yticklabels=labels,\n",
    "               title=f\"Threshold = {thresholds[t][k]:.2f}\",\n",
    "    #            ylabel='True label',\n",
    "               xlabel='Predicted label')\n",
    "        axs[k, t].set_yticklabels(labels, rotation='vertical', va='center', ha='center')\n",
    "        for i in range(cm.shape[0]):\n",
    "            for j in range(cm.shape[1]):\n",
    "                color = \"w\" if cm_normalized[i, j] > 0.5 else \"k\"\n",
    "                axs[k, t].text(j, i, int(cm[i, j].round()), ha=\"center\", va=\"bottom\", color=color)\n",
    "                axs[k, t].text(j, i, f'{cm_normalized[i, j]*100:.01f}%', ha=\"center\", va=\"top\", color=color)\n",
    "\n",
    "    cbar = fig.colorbar(im, ax=axs[k, 4], fraction=0.046, pad=0.04, ticks=[0.02, 0.2, 0.4, 0.6, 0.8, 0.98])\n",
    "    cbar.ax.set_yticklabels(np.linspace(0, 1, 6, dtype=np.float16))  # horizontal colorbar\n",
    "    axs[k, 0].set(ylabel='True label')\n",
    "    axs[k, 0].set(title=f'Ext. thresh. ({thresholds[0][k]:.03f})')\n",
    "    axs[k, 4].set(title=f'In-dom. thresh. ({thresholds[4][k]:.03f})')\n",
    "# plt.tight_layout()\n",
    "fig.savefig(output_path / f'confusion_matrices_{log_name}_{target_label}_{eval_name}_all_folds.svg',  format='svg', bbox_inches = 'tight', pad_inches = 0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e312e1d",
   "metadata": {},
   "outputs": [],
   "source": [
    "thresholds_test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4dfeba7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10, 5))\n",
    "# Split the predictions and labels by label\n",
    "preds_0 = y_preds[best_fold][y_trues[best_fold] == 0]\n",
    "preds_1 = y_preds[best_fold][y_trues[best_fold] == 1]\n",
    "\n",
    "# Plot the histograms\n",
    "plt.hist([preds_0, preds_1], bins=50, density=True, label=[\"nonMSIH\", \"MSIH\"])\n",
    "# plt.hist([preds_0, preds_1], bins=50, density=False, label=[\"nonMSIH\", \"MSIH\"], alpha=0.2)\n",
    "plt.xlabel(\"Prediction\")\n",
    "plt.ylabel(\"Density\")\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0f15ab5",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_folds = 5\n",
    "\n",
    "# Set up the subplots\n",
    "fig, axs = plt.subplots(1, num_folds, figsize=(num_folds * 5, 5))\n",
    "\n",
    "# Plot the histograms for each model\n",
    "for k, ax, in enumerate(axs):\n",
    "    preds_0 = y_preds[k][y_trues[k] == 0]\n",
    "    preds_1 = y_preds[k][y_trues[k] == 1]\n",
    "    \n",
    "    ax.hist([preds_0, preds_1], bins=50, density=True, label=[\"nonMSIH\", \"MSIH\"])\n",
    "    ax.set_xlabel(\"Prediction\")\n",
    "    ax.set_ylabel(\"Density\")\n",
    "    ax.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00e33fde",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_trues"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45713b35",
   "metadata": {},
   "source": [
    "### compute missing statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18442165",
   "metadata": {},
   "outputs": [],
   "source": [
    "pred_dfs = [pd.read_csv(p, dtype=str) for p in pred_csvs]\n",
    "\n",
    "y_trues = [df['ground_truth'] == f'{true_label:.1f}' for df in pred_dfs]\n",
    "y_preds = [pd.to_numeric(df[f'logits']) for df in pred_dfs]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0a454da",
   "metadata": {},
   "outputs": [],
   "source": [
    "pr_auc = []\n",
    "# Compute AUPRC mean and std\n",
    "for i in range(len(y_trues)):\n",
    "    prec, rec, _ = precision_recall_curve(y_trues[i], y_preds[i])\n",
    "    prec, rec = prec[::-1], rec[::-1]  # reverse order for interp\n",
    "    pr_auc.append(average_precision_score(y_trues[i], y_preds[i]))\n",
    "    \n",
    "np.mean(pr_auc), np.std(pr_auc)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd1c9873",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "## plot num samples curves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "402fc146",
   "metadata": {},
   "outputs": [],
   "source": [
    "# dir_path = Path('/lustre/groups/peng/workspace/sophia.wagner/logs/genalt/results')\n",
    "dir_path = Path('/home/haicu/sophia.wagner/projects/genetic_alteration_prediction')\n",
    "ext_cohort = 'Yorkshire_resections'\n",
    "# ext_cohort = 'Yorkshire_biopsies'\n",
    "model = 'transformer'\n",
    "feats = 'ctranspath'\n",
    "# ext_cohort = 'Yorkshire_resections'\n",
    "filename = dir_path / f'MSI_prediction_num_samples_{ext_cohort}_{model}_{feats}.csv'\n",
    "num_df = pd.read_csv(filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec4d78fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "#             # save results in data frame \n",
    "#             for ext_cohort in ext_cohorts:\n",
    "#                 results_csv = result_path / f'results_num_samples_{ext_cohort}_{args.model}_{args.feats}{args.name}.csv'\n",
    "#                 columns = [\"num samples\"] + [f'fold {i}' for i in range(5)]\n",
    "#                 new_results = pd.DataFrame(data=np.array([[n, *ext_auc_dict[ext_cohort][n]]]), columns=columns)\n",
    "\n",
    "#                 # append results to existing results data frame\n",
    "#                 if results_csv.is_file():\n",
    "#                     results_df = pd.read_csv(results_csv, dtype=str)\n",
    "#                     if n in results_df[\"num samples\"]:\n",
    "#                         continue\n",
    "#                     else:\n",
    "#                         results_df = results_df.append(new_results)\n",
    "#                 else:\n",
    "#                     results_df = new_results\n",
    "\n",
    "#                 results_df.to_csv(results_csv, sep=',', index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29df18a9-0841-4693-b0f3-59e76e95284b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c5c37ca-effc-4a29-b212-cc9a8c30dbec",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "560cbf0c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# cohort for evaluation\n",
    "# ext_cohort = 'Yorkshire-resections'\n",
    "ext_cohort = 'Yorkshire-biopsies'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da3bca1b",
   "metadata": {},
   "outputs": [],
   "source": [
    "model = 'vit-cls'\n",
    "feats = 'ctranspath'\n",
    "name = ''\n",
    "\n",
    "label = 'Ours'\n",
    "\n",
    "result_path = Path('/lustre/groups/peng/workspace/sophia.wagner/logs/genalt/results')\n",
    "filename = result_path / f'results_num_samples_{ext_cohort}_{model}_{feats}{name}.csv'\n",
    "\n",
    "num_df = pd.read_csv(filename)\n",
    "num_df = num_df.rename({'fold 4 ': 'fold 4'}, axis=1)\n",
    "num_df = num_df.rename({'training samples': 'num samples'}, axis=1)\n",
    "num_df = num_df.sort_values(by='num samples')\n",
    "num_df = num_df.drop_duplicates(subset='num samples')\n",
    "num_samples_array = [num_df[f'fold {i}'].values for i in range(5)]\n",
    "num_samples_x = num_df['num samples'].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb0c9f54",
   "metadata": {},
   "outputs": [],
   "source": [
    "model = 'attmil'\n",
    "# model = 'vit-cls'\n",
    "feats = 'ctranspath'\n",
    "name = ''\n",
    "\n",
    "label2 = 'AttentionMIL'\n",
    "\n",
    "filename = result_path / f'results_num_samples_{ext_cohort}_{model}_{feats}{name}.csv'\n",
    "\n",
    "num_df2 = pd.read_csv(filename)\n",
    "num_df2 = num_df2.rename({'fold 4 ': 'fold 4'}, axis=1)\n",
    "num_df2 = num_df2.rename({'training samples': 'num samples'}, axis=1)\n",
    "num_df2 = num_df2.sort_values(by='num samples')\n",
    "num_df2 = num_df2.drop_duplicates(subset='num samples')\n",
    "num_samples_array2 = [num_df2[f'fold {i}'].values for i in range(5)]\n",
    "num_samples_x2 = num_df2['num samples'].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ac317f2",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2bc0992",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_df2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b09e3d6-e828-4f2a-a2e5-76b09839ad98",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2013c02",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(5, 5))\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(5, 5))\n",
    "\n",
    "ax.plot((0, 10000), (0.9, 0.9), linewidth=0.5, color='gray')\n",
    "# ax.plot((0, 10000), (0.95, 0.95), linewidth=0.5, color='gray')\n",
    "# plt.plot((0, 10000), (0.95, 0.95), linestyle='--', linewidth=1, color='gray')\n",
    "\n",
    "\n",
    "num_samples = np.stack(num_samples_array)\n",
    "mean = num_samples.mean(axis=0)\n",
    "std = num_samples.std(axis=0)\n",
    "ax.plot(num_samples_x, mean, color='tab:blue', label=label)\n",
    "# ax.plot((0, 10000), (0.9686, 0.9686), linestyle='--', linewidth=1, color='tab:blue', label='Best model')\n",
    "plt.plot((0, 10000), (0.9072, 0.9072), linestyle='--', linewidth=1, color='tab:blue', label='Best model')\n",
    "ax.plot(num_samples_x, num_samples[0], alpha=0.3, color='tab:blue')\n",
    "ax.plot(num_samples_x, num_samples[1], alpha=0.3, color='tab:blue')\n",
    "ax.plot(num_samples_x, num_samples[2], alpha=0.3, color='tab:blue')\n",
    "ax.plot(num_samples_x, num_samples[3], alpha=0.3, color='tab:blue')\n",
    "ax.plot(num_samples_x, num_samples[4], alpha=0.3, color='tab:blue')\n",
    "ax.fill_between(num_samples_x, mean - std, mean + std, alpha=0.5, color='tab:blue')\n",
    "\n",
    "\n",
    "num_samples = np.stack(num_samples_array2)\n",
    "mean = num_samples.mean(axis=0)\n",
    "std = num_samples.std(axis=0)\n",
    "ax.plot(num_samples_x2, mean, color='tab:orange', label=label2)\n",
    "# plt.plot((0, 8000), (0.9633, 0.9633), linestyle='--', linewidth=1, color='tab:orange', label='Best model')\n",
    "ax.plot(num_samples_x2, num_samples[0], alpha=0.3, color='tab:orange')\n",
    "ax.plot(num_samples_x2, num_samples[1], alpha=0.3, color='tab:orange')\n",
    "ax.plot(num_samples_x2, num_samples[2], alpha=0.3, color='tab:orange')\n",
    "ax.plot(num_samples_x2, num_samples[3], alpha=0.3, color='tab:orange')\n",
    "ax.plot(num_samples_x2, num_samples[4], alpha=0.3, color='tab:orange')\n",
    "ax.fill_between(num_samples_x2, mean - std, mean + std, alpha=0.6, color='tab:orange')\n",
    "\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "\n",
    "ax.legend(loc=\"lower right\")\n",
    "ax.set(\n",
    "    ylim=(0.5, 1), \n",
    "    yticks=np.linspace(0.5, 1, 11), \n",
    "    ylabel='AUROC',\n",
    "    xlim=(min(num_samples_x)-2, 10100), \n",
    "    xscale='log', \n",
    "    xlabel='Number of samples in the training set (log scale)'\n",
    ")\n",
    "# xticks=num_samples_x, labels=[x for x in num_samples_x]\n",
    "# plt.savefig(output_path / f'num_samples_{ext_cohort}.svg', format='svg', bbox_inches = 'tight', pad_inches = 0)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b84307e",
   "metadata": {},
   "outputs": [],
   "source": [
    "figure_path = Path('/home/ubuntu/projects/idkidc/figures')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5feaf0e-f0aa-4000-90d0-bbc6b43450b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "cfg.name = 'full'\n",
    "\n",
    "# isMSIH experiments\n",
    "ext_cohort = 'YCR-BCIP-resections'\n",
    "# ext_cohort = 'YCR-BCIP-biopsies'\n",
    "# ext_cohort = 'MAINZ'\n",
    "cfg.target = 'isMSIH'\n",
    "best_model = 0.97\n",
    "max_samples = 10095\n",
    "\n",
    "# ext_cohort = 'Epi700'\n",
    "# cfg.target = 'BRAF'\n",
    "# best_model = 0.88\n",
    "# max_samples = 6829\n",
    "\n",
    "# ext_cohort = 'Epi700'\n",
    "# cfg.target = 'KRAS'\n",
    "# best_model = 0.80\n",
    "# max_samples = 6777"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2794963-6ee0-41f5-a99e-a7ab830edb65",
   "metadata": {},
   "outputs": [],
   "source": [
    "cfg.model = 'Transformer'\n",
    "label = 'Ours'\n",
    "\n",
    "cfg.logging_name = f'num_samples_{cfg.name}_{cfg.model}_MCO_{cfg.norm}_{cfg.target}'\n",
    "# base_path = Path('/lustre/groups/peng/workspace/sophia.wagner/logs/idkidc/') / cfg.logging_name\n",
    "base_path = Path('/home/ubuntu/projects/idkidc/temp/')\n",
    "filename = base_path / f'results_num_samples_{ext_cohort}_{cfg.model}_{cfg.feats}_{cfg.target}.csv'\n",
    "\n",
    "num_df = pd.read_csv(filename)\n",
    "num_df = num_df.rename({'fold 4 ': 'fold 4'}, axis=1)\n",
    "num_df = num_df.rename({'training samples': 'num samples'}, axis=1)\n",
    "num_df = num_df.sort_values(by='num samples')\n",
    "num_df = num_df.drop_duplicates(subset='num samples')\n",
    "num_samples_array = [num_df[f'fold {i}'].values for i in range(5)]\n",
    "num_samples_x = num_df['num samples'].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe43db02-a111-4b28-97db-9229eb739094",
   "metadata": {},
   "outputs": [],
   "source": [
    "cfg.model = 'AttentionMIL'\n",
    "label2 = 'AttentionMIL'\n",
    "\n",
    "cfg.logging_name = f'num_samples_{cfg.name}_{cfg.model}_MCO_{cfg.norm}_{cfg.target}'\n",
    "# base_path = Path('/lustre/groups/peng/workspace/sophia.wagner/logs/idkidc/') / cfg.logging_name\n",
    "base_path = Path('/home/ubuntu/projects/idkidc/temp/')\n",
    "filename = base_path / f'results_num_samples_{ext_cohort}_{cfg.model}_{cfg.feats}_{cfg.target}.csv'\n",
    "\n",
    "num_df2 = pd.read_csv(filename)\n",
    "num_df2 = num_df2.rename({'fold 4 ': 'fold 4'}, axis=1)\n",
    "num_df2 = num_df2.rename({'training samples': 'num samples'}, axis=1)\n",
    "num_df2 = num_df2.sort_values(by='num samples')\n",
    "num_df2 = num_df2.drop_duplicates(subset='num samples')\n",
    "num_samples_array2 = [num_df2[f'fold {i}'].values for i in range(5)]\n",
    "num_samples_x2 = num_df2['num samples'].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d055c5b8-83ce-451d-8c03-f53fc9d46aaa",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(5, 5))\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(5, 5))\n",
    "\n",
    "# ax.plot((0, 10000), (0.9, 0.9), linewidth=0.5, color='gray')\n",
    "# ax.plot((0, 10000), (0.95, 0.95), linewidth=0.5, color='gray')\n",
    "# plt.plot((0, 10000), (0.95, 0.95), linestyle='--', linewidth=1, color='gray')\n",
    "\n",
    "\n",
    "num_samples = np.stack(num_samples_array)\n",
    "mean = num_samples.mean(axis=0)\n",
    "std = num_samples.std(axis=0)\n",
    "ax.plot(num_samples_x, mean, color='tab:blue', label=label)\n",
    "ax.plot((0, 10000), (best_model, best_model), linestyle='--', linewidth=1, color='tab:blue', label='Best model')\n",
    "# plt.plot((0, 10000), (0.9072, 0.9072), linestyle='--', linewidth=1, color='tab:blue', label='Best model')\n",
    "ax.plot(num_samples_x, num_samples[0], alpha=0.3, color='tab:blue')\n",
    "ax.plot(num_samples_x, num_samples[1], alpha=0.3, color='tab:blue')\n",
    "ax.plot(num_samples_x, num_samples[2], alpha=0.3, color='tab:blue')\n",
    "ax.plot(num_samples_x, num_samples[3], alpha=0.3, color='tab:blue')\n",
    "ax.plot(num_samples_x, num_samples[4], alpha=0.3, color='tab:blue')\n",
    "ax.fill_between(num_samples_x, mean - std, mean + std, alpha=0.5, color='tab:blue')\n",
    "\n",
    "\n",
    "num_samples = np.stack(num_samples_array2)\n",
    "mean = num_samples.mean(axis=0)\n",
    "std = num_samples.std(axis=0)\n",
    "ax.plot(num_samples_x2, mean, color='tab:orange', label=label2)\n",
    "# plt.plot((0, 8000), (0.9633, 0.9633), linestyle='--', linewidth=1, color='tab:orange', label='Best model')\n",
    "ax.plot(num_samples_x2, num_samples[0], alpha=0.3, color='tab:orange')\n",
    "ax.plot(num_samples_x2, num_samples[1], alpha=0.3, color='tab:orange')\n",
    "ax.plot(num_samples_x2, num_samples[2], alpha=0.3, color='tab:orange')\n",
    "ax.plot(num_samples_x2, num_samples[3], alpha=0.3, color='tab:orange')\n",
    "ax.plot(num_samples_x2, num_samples[4], alpha=0.3, color='tab:orange')\n",
    "ax.fill_between(num_samples_x2, mean - std, mean + std, alpha=0.6, color='tab:orange')\n",
    "\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "\n",
    "ax.legend(loc=\"lower right\")\n",
    "ax.set(\n",
    "    # ylim=(0.5, 1), \n",
    "    # yticks=np.linspace(0.5, 1, 11), \n",
    "    ylabel='AUROC',\n",
    "    xlim=(min(num_samples_x)-2, 10100), \n",
    "    xscale='log', \n",
    "    xlabel='Number of samples in the training set (log scale)'\n",
    ")\n",
    "# xticks=num_samples_x, labels=[x for x in num_samples_x]\n",
    "# plt.savefig(output_path / f'num_samples_{ext_cohort}.svg', format='svg', bbox_inches = 'tight', pad_inches = 0)\n",
    "plt.savefig(figure_path / f'num_samples_{ext_cohort}_{cfg.target}.svg', format='svg', bbox_inches = 'tight', pad_inches = 0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11869c06",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_samples = np.stack(num_samples_array)\n",
    "mean = num_samples.mean(axis=0)\n",
    "std = num_samples.std(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8627278e",
   "metadata": {},
   "outputs": [],
   "source": [
    "mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e4de35d",
   "metadata": {},
   "outputs": [],
   "source": [
    "std"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf80aaad",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_samples = np.stack(num_samples_array)\n",
    "mean = num_samples.mean(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "536ee00c",
   "metadata": {},
   "outputs": [],
   "source": [
    "mean, num_samples_x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd749df6",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_samples = np.stack(num_samples_array2)\n",
    "mean = num_samples.mean(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25af967e",
   "metadata": {},
   "outputs": [],
   "source": [
    "mean, num_samples_x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c59f6b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "t_ctp_n = np.array([50, 100, 250, 500, 1000, 1500, 2000, 2500, 3000, 3500, 4000, 4500, 5000, 5500, 6000, 6500, 7000, 7500, 8000])\n",
    "t_ctp = np.array([\n",
    "    [0.7652, 0.8778, 0.7515, 0.7836, 0.8378,],\n",
    "    [0.8041, 0.8283, 0.9115, 0.7765, 0.7524,],\n",
    "    [0.9340, 0.9329, 0.8734, 0.9070, 0.9298,],\n",
    "    [0.9313, 0.9485, 0.9064, 0.8828, 0.9169,],\n",
    "    [0.9601, 0.9477, 0.9599, 0.9544, 0.9588,],\n",
    "    [0.9452, 0.9513, 0.9575, 0.9680, 0.9475,],\n",
    "    [0.9634, 0.9727, 0.9631, 0.9518, 0.9489,],\n",
    "    [0.9676, 0.9677, 0.9551, 0.9600, 0.9635,],\n",
    "    [0.9583, 0.9637, 0.9565, 0.9655, 0.9591,],\n",
    "    [0.9619, 0.9637, 0.9707, 0.9633, 0.9654,],\n",
    "    [0.9657, 0.9573, 0.9683, 0.9726, 0.9583,],\n",
    "    [0.9647, 0.9672, 0.9729, 0.9700, 0.9651,],\n",
    "    [0.9668, 0.9768, 0.9703, 0.9613, 0.9647,],\n",
    "    [0.9669, 0.9682, 0.9604, 0.9667, 0.9594,],\n",
    "    [0.9728, 0.9702, 0.9682, 0.9746, 0.9663,],\n",
    "    [0.9636, 0.9646, 0.9696, 0.9785, 0.9730,],\n",
    "    [0.9709, 0.9702, 0.9705, 0.9781, 0.9679,],\n",
    "    [0.9724, 0.9662, 0.9677, 0.9711, 0.9741,],\n",
    "    [0.9692, 0.9684, 0.9706, 0.9590, 0.9686,],\n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89a8af9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_samples = t_ctp\n",
    "mean = num_samples.mean(axis=1)\n",
    "std = num_samples.std(axis=1)\n",
    "plt.plot(t_ctp_n, mean)\n",
    "plt.plot(t_ctp_n, num_samples[:, 0], alpha=0.2)\n",
    "plt.plot(t_ctp_n, num_samples[:, 1], alpha=0.2)\n",
    "plt.plot(t_ctp_n, num_samples[:, 2], alpha=0.2)\n",
    "plt.plot(t_ctp_n, num_samples[:, 3], alpha=0.2)\n",
    "plt.plot(t_ctp_n, num_samples[:, 4], alpha=0.2)\n",
    "plt.fill_between(t_ctp_n, mean - std, mean + std, alpha=0.5)\n",
    "\n",
    "# plt.xlim(1, num_samples.shape[1]+1)\n",
    "plt.xticks(ticks=t_ctp_n) # labels=[x for x in range(500, num_samples.shape[1]*500+1, 500)])\n",
    "plt.xscale('log')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01896018",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_samples = t_ctp\n",
    "mean = num_samples.mean(axis=1)\n",
    "std = num_samples.std(axis=1)\n",
    "plt.plot(t_ctp_n, mean)\n",
    "plt.plot(t_ctp_n, num_samples[:, 0], alpha=0.2)\n",
    "plt.plot(t_ctp_n, num_samples[:, 1], alpha=0.2)\n",
    "plt.plot(t_ctp_n, num_samples[:, 2], alpha=0.2)\n",
    "plt.plot(t_ctp_n, num_samples[:, 3], alpha=0.2)\n",
    "plt.plot(t_ctp_n, num_samples[:, 4], alpha=0.2)\n",
    "plt.fill_between(t_ctp_n, mean - std, mean + std, alpha=0.5)\n",
    "\n",
    "# plt.xlim(1, num_samples.shape[1]+1)\n",
    "plt.xticks(ticks=t_ctp_n) # labels=[x for x in range(500, num_samples.shape[1]*500+1, 500)])\n",
    "plt.xscale('log')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e440129b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "380ccae9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import module for data manipulation\n",
    "import pandas as pd\n",
    "# Import module for linear algebra\n",
    "import numpy as np\n",
    "# Import module for data simulation\n",
    "from sklearn.datasets import make_classification     # Create a synthetic dataframe\n",
    "from sklearn.linear_model import LogisticRegression  # Classification model\n",
    "from sklearn.model_selection import train_test_split # Split the dataframe\n",
    "from sklearn.metrics import roc_curve                # Calculate the ROC curve\n",
    "from sklearn.metrics import precision_recall_curve   # Calculate the Precision-Recall curve\n",
    "from sklearn.metrics import f1_score                 # Calculate the F-score\n",
    "# Import module for data visualization\n",
    "from plotnine import *\n",
    "import plotnine\n",
    "\n",
    "# Generate the dataset\n",
    "X, y = make_classification(n_samples = 10000, n_features = 2, n_redundant = 0,\n",
    "                           n_clusters_per_class = 1, weights = [0.99], flip_y = 0, random_state = 0)\n",
    "\n",
    "# Data partitioning\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.5, random_state=0, stratify=y)\n",
    "\n",
    "# Fit the model\n",
    "reglogModel = LogisticRegression(random_state = 0)\n",
    "reglogModel.fit(X_train, y_train)\n",
    "\n",
    "# Predict the probabilities\n",
    "y_pred = reglogModel.predict_proba(X_test)\n",
    "\n",
    "# Get the probabilities for positive class\n",
    "y_pred = y_pred[:, 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94a47114",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the ROC curve\n",
    "fpr, tpr, thresholds = roc_curve(y_test, y_pred)\n",
    "\n",
    "# Plot the ROC curve\n",
    "df_fpr_tpr = pd.DataFrame({'FPR':fpr, 'TPR':tpr, 'Threshold':thresholds})\n",
    "df_fpr_tpr.head()\n",
    "\n",
    "# Create the data viz\n",
    "plotnine.options.figure_size = (8, 4.8)\n",
    "(\n",
    "    ggplot(data = df_fpr_tpr)+\n",
    "    geom_point(aes(x = 'FPR',\n",
    "                   y = 'TPR'),\n",
    "               size = 0.4)+\n",
    "    geom_line(aes(x = 'FPR',\n",
    "                  y = 'TPR'))+\n",
    "    labs(title = 'ROC Curve')+\n",
    "    xlab('False Positive Rate')+\n",
    "    ylab('True Positive Rate')+\n",
    "    theme_minimal()\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61b06521",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the G-mean\n",
    "gmean = np.sqrt(tpr * (1 - fpr))\n",
    "\n",
    "# Find the optimal threshold\n",
    "index = np.argmax(gmean)\n",
    "thresholdOpt = round(thresholds[index], ndigits = 4)\n",
    "gmeanOpt = round(gmean[index], ndigits = 4)\n",
    "fprOpt = round(fpr[index], ndigits = 4)\n",
    "tprOpt = round(tpr[index], ndigits = 4)\n",
    "print('Best Threshold: {} with G-Mean: {}'.format(thresholdOpt, gmeanOpt))\n",
    "print('FPR: {}, TPR: {}'.format(fprOpt, tprOpt))\n",
    "\n",
    "# Create data viz\n",
    "plotnine.options.figure_size = (8, 4.8)\n",
    "(\n",
    "    ggplot(data = df_fpr_tpr)+\n",
    "    geom_point(aes(x = 'FPR',\n",
    "                   y = 'TPR'),\n",
    "               size = 0.4)+\n",
    "    # Best threshold\n",
    "    geom_point(aes(x = fprOpt,\n",
    "                   y = tprOpt),\n",
    "               color = '#981220',\n",
    "               size = 4)+\n",
    "    geom_line(aes(x = 'FPR',\n",
    "                  y = 'TPR'))+\n",
    "    geom_text(aes(x = fprOpt,\n",
    "                  y = tprOpt),\n",
    "              label = 'Optimal threshold \\n for class: {}'.format(thresholdOpt),\n",
    "              nudge_x = 0.14,\n",
    "              nudge_y = -0.10,\n",
    "              size = 10,\n",
    "              fontstyle = 'italic')+\n",
    "    labs(title = 'ROC Curve')+\n",
    "    xlab('False Positive Rate (FPR)')+\n",
    "    ylab('True Positive Rate (TPR)')+\n",
    "    theme_minimal()\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7fa7db0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the sensitivity threshold\n",
    "sensitivity_thresh = 0.90\n",
    "\n",
    "# Find the optimal threshold\n",
    "index = np.argwhere(tpr >= sensitivity_thresh)[0].item()\n",
    "thresholdOpt = round(thresholds[index], ndigits = 4)\n",
    "gmeanOpt = round(gmean[index], ndigits = 4)\n",
    "fprOpt = round(fpr[index], ndigits = 4)\n",
    "tprOpt = round(tpr[index], ndigits = 4)\n",
    "print('Best Threshold: {} with G-Mean: {}'.format(thresholdOpt, gmeanOpt))\n",
    "print('FPR: {}, TPR: {}'.format(fprOpt, tprOpt))\n",
    "\n",
    "# Create data viz\n",
    "plotnine.options.figure_size = (8, 4.8)\n",
    "(\n",
    "    ggplot(data = df_fpr_tpr)+\n",
    "    geom_point(aes(x = 'FPR',\n",
    "                   y = 'TPR'),\n",
    "               size = 0.4)+\n",
    "    # Best threshold\n",
    "    geom_point(aes(x = fprOpt,\n",
    "                   y = tprOpt),\n",
    "               color = '#981220',\n",
    "               size = 4)+\n",
    "    geom_line(aes(x = 'FPR',\n",
    "                  y = 'TPR'))+\n",
    "    geom_text(aes(x = fprOpt,\n",
    "                  y = tprOpt),\n",
    "              label = 'Optimal threshold \\n for class: {}'.format(thresholdOpt),\n",
    "              nudge_x = 0.14,\n",
    "              nudge_y = -0.10,\n",
    "              size = 10,\n",
    "              fontstyle = 'italic')+\n",
    "    labs(title = 'ROC Curve')+\n",
    "    xlab('False Positive Rate (FPR)')+\n",
    "    ylab('True Positive Rate (TPR)')+\n",
    "    theme_minimal()\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f49d5715",
   "metadata": {},
   "outputs": [],
   "source": [
    "thresholds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e3d6827",
   "metadata": {},
   "outputs": [],
   "source": [
    "index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad01af8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "tpr[tpr >= 0.9]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61383b25",
   "metadata": {},
   "outputs": [],
   "source": [
    "tpr[tpr >= 0.9][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7558752b",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.argwhere(tpr >= 0.9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "449ae99a",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.argwhere(tpr >= 0.9)[0].item()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6553976",
   "metadata": {},
   "outputs": [],
   "source": [
    "thresholds[np.where(tpr == tpr[tpr >= 0.9][0])[0][0]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ad2e624",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63febbde",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77697acb",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6dd560b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "d = {'key':3}, {'te': 45\n",
    "               }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03c21a32",
   "metadata": {},
   "outputs": [],
   "source": [
    "d"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "578d306d",
   "metadata": {},
   "source": [
    "# Analysis of clinicopathological features"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff40c2bc",
   "metadata": {},
   "source": [
    "## Plot AUROCs for subgroups"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b617690f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# output_path = Path('/home/haicu/sophia.wagner/projects/idkidc/figures') \n",
    "output_path = Path('/Users/sophia.wagner/Documents/PhD/projects/2022_MSI_transformer/figures/curves') \n",
    "\n",
    "# test_cohorts = 'YCR-BCIP-resections'\n",
    "# cohort = 'YCR-BCIP-resections'\n",
    "# cohort_label = 'YCR-BCIP'\n",
    "\n",
    "cfg.norm = norm = 'histaugan'\n",
    "cfg.feats = feats =  'ctranspath'\n",
    "cfg.target = 'isMSIH'\n",
    "categories = ['Not mut.', 'Mutat.', 'nonMSIH', 'MSIH', 'WT', 'MUT', 'wt', 'MT']\n",
    "target_labels = [cfg.target]\n",
    "true_label = 1.\n",
    "clini_info = {'AGE': None, 'GENDER': None, 'LEFT_RIGHT': None, 'STAGE': None}\n",
    "\n",
    "label_dict = {\n",
    "    'Not mut.': 0,\n",
    "    'Mutat.': 1,\n",
    "    'nonMSIH': 0,\n",
    "    'MSIH': 1,\n",
    "    'WT': 0,\n",
    "    'MUT': 1,\n",
    "    'wt': 0,\n",
    "    'MT': 1,\n",
    "    'left': 1,\n",
    "    'right': 0,\n",
    "    'female': 1,\n",
    "    'male': 0,\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "790b953e",
   "metadata": {},
   "outputs": [],
   "source": [
    "csv_paths = result_path.glob(f'fold*/outputs_{test_cohorts}.csv')\n",
    "pred_csvs = list(csv_paths)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "899d1060",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_config = Path('data_config.yaml')\n",
    "with open(data_config, 'r') as f:\n",
    "    data_config = yaml.safe_load(f)\n",
    "\n",
    "    clini_table = Path(data_config[cohort]['clini_table'])\n",
    "    slide_csv = Path(data_config[cohort]['slide_csv'])\n",
    "    feature_dir = Path(data_config[cohort]['feature_dir'][norm][feats])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40bf4fb0",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_excel('/Users/sophia.wagner/Downloads/clini_tables/FOXTROT-CRC-DX_CLINI_12062023.xlsx')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "391de8c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# from data import get_cohort_df\n",
    "\n",
    "# df = get_cohort_df(clini_table, slide_csv, feature_dir, [cfg.target], categories, cohort, clini_info) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19e9becc",
   "metadata": {},
   "outputs": [],
   "source": [
    "df['LEFT_RIGHT'].value_counts(dropna=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e758610",
   "metadata": {},
   "outputs": [],
   "source": [
    "df['GENDER'].value_counts(dropna=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61e82e9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(df['AGE'].astype(int))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b45def5d",
   "metadata": {},
   "outputs": [],
   "source": [
    "df['STAGE'].value_counts(dropna=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29a10165",
   "metadata": {},
   "outputs": [],
   "source": [
    "df['AGE'] = df['AGE'].astype(int)\n",
    "df['STAGE'] = pd.to_numeric(df['STAGE'], errors='coerce', downcast='integer')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45173657",
   "metadata": {},
   "outputs": [],
   "source": [
    "df['prepost'].value_counts(dropna=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "363d62b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# --------------------------\n",
    "# state the condition\n",
    "# --------------------------\n",
    "\n",
    "# age\n",
    "# lower = 0\n",
    "# upper = 60\n",
    "# df = df[df['AGE'] < upper]\n",
    "# df = df[df['AGE'] >= lower]\n",
    "\n",
    "# # gender\n",
    "# gender = 'female'\n",
    "# df = df[df['GENDER'] == gender]\n",
    "\n",
    "# # tumor site\n",
    "# site = 'left'\n",
    "# df = df[df['LEFT_RIGHT'] == site]\n",
    "\n",
    "# # stage\n",
    "# stage = 3\n",
    "# df = df[df['STAGE'] == stage]\n",
    "\n",
    "# pre/post treatment (FOXTROT)\n",
    "prepost = 'Control'\n",
    "df = df[df['prepost'] == prepost]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58437580",
   "metadata": {},
   "outputs": [],
   "source": [
    "# test_df = pd.read_csv(pred_csvs[0], dtype=str).rename({'patient': 'PATIENT'}, axis=1)\n",
    "# logits = test_df['logits'].astype(float)\n",
    "# test_df['predictions'] = torch.tensor(logits.values).sigmoid().round()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "681032b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "df['PATIENT'] = df['PATIENT'].astype(str)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c676454",
   "metadata": {},
   "outputs": [],
   "source": [
    "pred_dfs = [pd.read_csv(p, dtype=str).rename({'patient': 'PATIENT'}, axis=1) for p in pred_csvs]\n",
    "pred_dfs = [pd.merge(pred_dfs[i], df, on=\"PATIENT\", how='inner', ) for i in range(len(pred_dfs))]\n",
    "\n",
    "y_trues = [df[\"ground_truth\"] == f'{true_label:.1f}' for df in pred_dfs]\n",
    "# y_trues = [label_dict[df[cfg.target]] == 'ground_truth' for df in pred_dfs]\n",
    "y_preds = [pd.to_numeric(df[f'logits']) for df in pred_dfs]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72a9bd75",
   "metadata": {},
   "outputs": [],
   "source": [
    "tprs = []\n",
    "aucs = []\n",
    "mean_fpr = np.linspace(0, 1, 100)\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(5, 5))\n",
    "\n",
    "for i in range(len(y_trues)):\n",
    "    fpr, tpr, _ = roc_curve(y_trues[i], y_preds[i])\n",
    "    roc_auc = roc_auc_score(y_trues[i], y_preds[i])\n",
    "    \n",
    "#     ax.plot(fpr, tpr, label=f'AUC = {roc_auc:0.2f}', alpha=0.3, color='black')\n",
    "\n",
    "\n",
    "    interp_tpr = np.interp(mean_fpr, fpr, tpr)\n",
    "    interp_tpr[0] = 0.0\n",
    "    tprs.append(interp_tpr)\n",
    "    aucs.append(roc_auc)\n",
    "ax.plot([0, 1], [0, 1], linestyle=\"--\", lw=1, color=\"gray\", label=\"Chance\", alpha=0.8)\n",
    "\n",
    "mean_tpr = np.mean(tprs, axis=0)\n",
    "mean_tpr[-1] = 1.0\n",
    "mean_auc = auc(mean_fpr, mean_tpr)\n",
    "std_auc = np.std(aucs)\n",
    "# plot mean curve\n",
    "ax.plot(\n",
    "    mean_fpr,\n",
    "    mean_tpr,\n",
    "    color=\"black\",\n",
    "    label=r\"Mean ROC (AUC = %0.2f $\\pm$ %0.3f)\" % (mean_auc, std_auc),\n",
    "    lw=2,\n",
    "    alpha=0.8,\n",
    ")\n",
    "\n",
    "std_tpr = np.std(tprs, axis=0)\n",
    "tprs_upper = np.minimum(mean_tpr + std_tpr, 1)\n",
    "tprs_lower = np.maximum(mean_tpr - std_tpr, 0)\n",
    "ax.fill_between(\n",
    "    mean_fpr,\n",
    "    tprs_lower,\n",
    "    tprs_upper,\n",
    "    color=\"grey\",\n",
    "    alpha=0.5,\n",
    "    label=r\"$\\pm$ 1 std. dev.\",\n",
    ")\n",
    "\n",
    "ax.set(\n",
    "    xlim=[-0.05, 1.05],\n",
    "    ylim=[-0.05, 1.05],\n",
    "    # title=f\"ROC for patients with: {lower} ≤ age < {upper} (n={len(df)})\",\n",
    "    # title=f\"ROC for patients with: gender = {gender} (n={len(df)})\",\n",
    "    # title=f\"ROC for patients with: tumor position = {site} (n={len(df)})\",\n",
    "    title=f\"ROC for patients in control group (n={len(df)})\",\n",
    ")\n",
    "\n",
    "ax.set_xlabel('1 - Specificity')\n",
    "ax.set_ylabel('Sensitivity')\n",
    "\n",
    "ax.legend(loc=\"lower right\")\n",
    "# fig.savefig(output_path / f'auroc_{test_cohorts}_{cfg.target}_age_{lower}_{upper}.svg', format='svg', bbox_inches = 'tight', pad_inches = 0)\n",
    "# fig.savefig(output_path / f'auroc_{test_cohorts}_{cfg.target}_gender_{gender}.svg', format='svg', bbox_inches = 'tight', pad_inches = 0)\n",
    "# fig.savefig(output_path / f'auroc_{test_cohorts}_{cfg.target}_site_{site}.svg', format='svg', bbox_inches = 'tight', pad_inches = 0)\n",
    "# fig.savefig(output_path / f'auroc_{test_cohorts}_{cfg.target}_stage_{stage}.svg', format='svg', bbox_inches = 'tight', pad_inches = 0)\n",
    "fig.savefig(output_path / f'auroc_{test_cohorts}_{cfg.target}_prepost_{prepost}.svg', format='svg', bbox_inches = 'tight', pad_inches = 0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0282813d",
   "metadata": {},
   "outputs": [],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88f1af58",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "idkidc",
   "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}