echino_loader.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "97d4f49a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imports complete\n"
     ]
    }
   ],
   "source": [
    "import glob\n",
    "import os\n",
    "import pickle\n",
    "import sys\n",
    "from tqdm import tqdm\n",
    "from pathlib import Path \n",
    "\n",
    "# change these to your local paths\n",
    "sys.path.append('~/dev/pathml') \n",
    "sys.path.append('~/dev/pathml/Pytorch-UNet')\n",
    "\n",
    "from pathml.slide import Slide\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import copy\n",
    "import torch\n",
    "import torchvision\n",
    "from torchvision import datasets, transforms\n",
    "from sklearn.metrics import accuracy_score, balanced_accuracy_score, f1_score, precision_score, recall_score, roc_auc_score, roc_curve\n",
    "import time\n",
    "\n",
    "print(\"Imports complete\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "d5ea2d1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def caseRandomizer():\n",
    "    wsi_path = Path(\"/pml/wsi_data/\")\n",
    "    dirs = os.listdir(wsi_path)\n",
    "    print(f\"Loaded {len(dirs)} images from {wsi_path} \\r\\n\")\n",
    "    res = np.char.rstrip(dirs, '.tif')\n",
    "    df = pd.DataFrame(res)\n",
    "    train_cases, val_cases, test_cases = np.split(df.sample(frac=1, random_state=123), [int(.6*len(df)), int(.8*len(df))])\n",
    "\n",
    "    train_cases = train_cases[0].tolist()\n",
    "    val_cases = val_cases[0].tolist()\n",
    "    test_cases = test_cases[0].tolist()\n",
    "\n",
    "    print(\"Training dataset: \", sorted(train_cases), \"\\r\\n\")\n",
    "    print(\"Validation dataset: \", sorted(val_cases), \"\\r\\n\")\n",
    "    print(\"Ttest: \", sorted(test_cases), \"\\r\\n\")\n",
    "\n",
    "    return train_cases, val_cases, test_cases\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "7f080377",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loaded 92 images from /pml/wsi_data \r\n",
      "\n",
      "Paths are all set, dataset splitting done \r\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Since images and annotations are inputs, we load them from out of the working directory\n",
    "wsi_path = Path(\"/pml/wsi_data/\")\n",
    "annotations_dir_path = Path(\"/pml/annotations/\")\n",
    "\n",
    "# The analysis path being an output remains to be local\n",
    "analysis_dir_path = '~/dev/pml-echino/analysis/'\n",
    "\n",
    "tiles_path = \"/pml/\"\n",
    "\n",
    "# Storing pathml slides (result of tissue detection) outside, as they don`t change if images and annotations are the same\n",
    "pathml_slide_dir_path = Path(\"/pml/slides\")\n",
    "\n",
    "os.makedirs(analysis_dir_path, exist_ok=True)\n",
    "os.makedirs(os.path.join(analysis_dir_path, 'classification_results'), exist_ok=True)\n",
    "os.makedirs(os.path.join(analysis_dir_path, 'segmentation_results'), exist_ok=True)\n",
    "\n",
    "# Splitting train, validation and test dataset randomly\n",
    "train_cases, val_cases, test_cases = caseRandomizer()\n",
    "\n",
    "train_wsi_paths = [os.path.join(wsi_path, train_case+'.tif') for train_case in train_cases]\n",
    "val_wsi_paths = [os.path.join(wsi_path, val_case+'.tif') for val_case in val_cases]\n",
    "test_wsi_paths = [os.path.join(wsi_path, test_case+'.tif') for test_case in test_cases]\n",
    "\n",
    "\n",
    "print(\"Paths are all set, dataset splitting done\", \"\\r\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98173a5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def tissueDetection(tile_size):\n",
    "    start_time = time.time()\n",
    "\n",
    "    for wsi_path in train_wsi_paths+val_wsi_paths+test_wsi_paths:\n",
    "        case = Path(wsi_path).stem\n",
    "        # check for existing pml files\n",
    "        case_pml_file = pathml_slide_dir_path / Path(case+'.pml')\n",
    "        # print(f\"DEBUG: {type(case_pml_file)}, {case_pml_file}\")\n",
    "        if case_pml_file.is_file():\n",
    "            print(f\"pml file exists for {case} as {case_pml_file}\")\n",
    "            continue\n",
    "        \n",
    "        pathml_slide = Slide(wsi_path, level=0).setTileProperties(tileSize=tile_size)\n",
    "        pathml_slide.detectTissue(numWorkers=0) \n",
    "        pathml_slide.detectForeground()\n",
    "\n",
    "        if 'echino' in case:\n",
    "            annotation_path = os.path.join(annotations_dir_path, case+'.geojson')\n",
    "            \n",
    "            print(annotation_path)\n",
    "            pathml_slide.addAnnotations(annotation_path, negativeClass='negative')\n",
    "\n",
    "        pathml_slide.save(folder=pathml_slide_dir_path)\n",
    "\n",
    "    time_elapsed = time.time() - start_time\n",
    "    print('Complete in {:.0f}m {:.0f}s'.format(time_elapsed // 60, time_elapsed % 60))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "2f5b277a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def detectionchecking():\n",
    "    for case in val_cases:\n",
    "        pathml_slide = Slide(os.path.join(pathml_slide_dir_path, case+'.pml'))\n",
    "        pathml_slide.visualizeThumbnail(folder=os.path.join(analysis_dir_path, 'classification_results'))\n",
    "        pathml_slide.visualizeTissueDetection(folder=os.path.join(analysis_dir_path, 'classification_results'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "83486e04",
   "metadata": {},
   "outputs": [],
   "source": [
    "def tilextractor(tile_size):\n",
    "    \n",
    "    print(\"Tile size is set to\", tile_size, \"pixels\")\n",
    "    \n",
    "    numTilesToExtract = 500\n",
    "\n",
    "    start_time = time.time()\n",
    "    \n",
    "    global global_channel_sums\n",
    "    global global_channel_squared_sums\n",
    "    global global_tile_count\n",
    "\n",
    "    global_channel_sums = np.zeros(3)\n",
    "    global_channel_squared_sums = np.zeros(3)\n",
    "    global_tile_count = 0\n",
    "\n",
    "    for case in train_cases + val_cases:\n",
    "        \n",
    "        print(case)\n",
    "               \n",
    "        slide_file = str(pathml_slide_dir_path / Path(case+'.pml'))\n",
    "        # print(\"slide file: \", slide_file)\n",
    "        pathml_slide = Slide(slide_file)\n",
    "\n",
    "        if 'echino' in case: #####change to echino\n",
    "            channel_data = pathml_slide.extractAnnotationTiles(tiles_path, \n",
    "                                                               tileAnnotationOverlapThreshold=0.7,\n",
    "                                                               otherClassNames='non_echinococcus',\n",
    "                                                               numTilesToExtractPerClass=numTilesToExtract, \n",
    "                                                               extractSegmentationMasks=True, \n",
    "                                                               tissueLevelThreshold=0.995, \n",
    "                                                               foregroundLevelThreshold=88)\n",
    "\n",
    "            global_channel_sums = np.add(global_channel_sums, channel_data['channel_sums'])\n",
    "            global_channel_squared_sums = np.add(global_channel_squared_sums, channel_data['channel_squared_sums'])\n",
    "            global_tile_count = global_tile_count + channel_data['num_tiles']\n",
    "\n",
    "        else: \n",
    "            channel_data = pathml_slide.extractRandomUnannotatedTiles(tiles_path, \n",
    "                                                                      numTilesToExtract=numTilesToExtract, \n",
    "                                                                      otherClassNames='Tumor',\n",
    "                                                                      unannotatedClassName='non_echinococcus',\n",
    "                                                                      extractSegmentationMasks=True, \n",
    "                                                                      tissueLevelThreshold=0.995, \n",
    "                                                                      foregroundLevelThreshold=88)\n",
    "\n",
    "            global_channel_sums = np.add(global_channel_sums, channel_data['channel_sums'])\n",
    "            global_channel_squared_sums = np.add(global_channel_squared_sums, channel_data['channel_squared_sums'])\n",
    "            global_tile_count = global_tile_count + channel_data['num_tiles']\n",
    "            \n",
    "    %store global_channel_sums\n",
    "    %store global_channel_squared_sums\n",
    "    %store global_tile_count\n",
    "\n",
    "    time_elapsed = time.time() - start_time\n",
    "    print('Complete in {:.0f}m {:.0f}s'.format(time_elapsed // 60, time_elapsed % 60))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "27d982cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "def tileworks(batch_size):\n",
    "    \n",
    "    total_pixels_per_channel = global_tile_count * tile_size * tile_size\n",
    "    global_channel_means = np.divide(global_channel_sums, total_pixels_per_channel)\n",
    "    global_channel_squared_means = np.divide(global_channel_squared_sums, total_pixels_per_channel)\n",
    "    global_channel_variances = np.subtract(global_channel_squared_means, np.square(global_channel_means))\n",
    "    global_channel_stds = np.sqrt(global_channel_variances * (total_pixels_per_channel / (total_pixels_per_channel-1)))\n",
    "    \n",
    "    global means_and_stds\n",
    "    \n",
    "    means_and_stds = {'channel_means': global_channel_means.tolist(), 'channel_stds': global_channel_stds.tolist()}\n",
    "    \n",
    "    pickle.dump(means_and_stds,\n",
    "        open(os.path.join(analysis_dir_path, 'classification_results', 'trainval_channel_means_and_stds.p'), 'wb'))\n",
    "\n",
    "    print('Channel means:', global_channel_means.tolist())\n",
    "    print('Channel standard deviations:', global_channel_stds.tolist())\n",
    "\n",
    "    echninococcus_tiles = glob.glob(os.path.join(tiles_path, 'tiles', '*', 'Tumor', '*.jpg'))\n",
    "    non_echninococcus_tiles = glob.glob(os.path.join(tiles_path, 'tiles', '*', 'non_echinococcus', '*.jpg'))\n",
    "\n",
    "    print('Total number of echinococcus tiles:', len(echninococcus_tiles))\n",
    "    print('Total number of non-echninococcus tiles:', len(non_echninococcus_tiles))\n",
    "\n",
    "    data_transforms = {\n",
    "    'train': transforms.Compose([\n",
    "        transforms.Resize(224),\n",
    "        transforms.RandomVerticalFlip(),\n",
    "        transforms.RandomHorizontalFlip(),\n",
    "        transforms.ColorJitter(brightness=0, contrast=0, saturation=1, hue=.5),\n",
    "        transforms.ToTensor(),\n",
    "        transforms.Normalize(means_and_stds['channel_means'], means_and_stds['channel_stds'])\n",
    "    ]),\n",
    "    'val': transforms.Compose([\n",
    "        transforms.Resize(224),\n",
    "        transforms.ToTensor(),\n",
    "        transforms.Normalize(means_and_stds['channel_means'], means_and_stds['channel_stds'])\n",
    "    ]),\n",
    "    }\n",
    "    \n",
    "    global dataset_sizes\n",
    "    global dataloaders\n",
    "    \n",
    "    train_dataset = torch.utils.data.ConcatDataset([datasets.ImageFolder(os.path.join(\n",
    "            tiles_path, 'tiles', train_case), data_transforms['train']) for train_case in train_cases])\n",
    "    val_dataset = torch.utils.data.ConcatDataset([datasets.ImageFolder(os.path.join(\n",
    "            tiles_path, 'tiles', val_case), data_transforms['val']) for val_case in val_cases])\n",
    "\n",
    "    image_datasets = {'train': train_dataset, 'val': val_dataset}\n",
    "    dataset_sizes = {x: len(image_datasets[x]) for x in ['train', 'val']}\n",
    "    \n",
    "    dataloaders = {}\n",
    "    batch_size = 48\n",
    "    dataloaders['train'] = torch.utils.data.DataLoader(\n",
    "        image_datasets['train'], batch_size=batch_size, num_workers=0, shuffle=True)\n",
    "    dataloaders['val'] = torch.utils.data.DataLoader(\n",
    "        image_datasets['val'], batch_size=batch_size, num_workers=0, shuffle=True)\n",
    "\n",
    "    %store dataloaders\n",
    "    %store means_and_stds\n",
    "    %store dataset_sizes\n",
    "  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "a7782c0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def vizBatch():\n",
    "\n",
    "    def visualize_batch(inp, title=None):\n",
    "        inp = inp.numpy().transpose((1, 2, 0))\n",
    "        mean = means_and_stds['channel_means']\n",
    "        std = means_and_stds['channel_stds']\n",
    "        inp = std * inp + mean\n",
    "        inp = np.clip(inp, 0, 1)\n",
    "        plt.imshow(inp)\n",
    "        if title is not None:\n",
    "            plt.title(title)\n",
    "        plt.pause(0.001)\n",
    "\n",
    "    inputs, classes = next(iter(dataloaders['train']))\n",
    "    out = torchvision.utils.make_grid(inputs)\n",
    "    visualize_batch(out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "8d6aac34",
   "metadata": {},
   "outputs": [],
   "source": [
    "def train_classification_model(model, criterion, optimizer, scheduler=False, num_epochs=30):\n",
    "    best_model_state_dict = copy.deepcopy(model.state_dict())\n",
    "    best_accuracy = 0.0\n",
    "    learning_stats = {'train': [], 'val': []}\n",
    "        \n",
    "    for epoch in range(num_epochs):\n",
    "        for phase in ['train', 'val']:\n",
    "            if phase == 'train':\n",
    "                model.train()\n",
    "            else:\n",
    "                model.eval()\n",
    "\n",
    "            epoch_ground_truth = []\n",
    "            epoch_predictions = []\n",
    "            running_loss = 0.0#\n",
    "            \n",
    "           \n",
    "        \n",
    "            for inputs, labels in tqdm(dataloaders[phase]):\n",
    "                inputs = inputs.to(device)\n",
    "                labels = labels.to(device)\n",
    "                # print(labels)\n",
    "                \n",
    "                \n",
    "                optimizer.zero_grad()\n",
    "                \n",
    "                with torch.set_grad_enabled(phase == 'train'):\n",
    "                    outputs = model(inputs)\n",
    "                    _, preds = torch.max(outputs, 1)\n",
    "\n",
    "                    loss = criterion(outputs, labels)\n",
    "                    \n",
    "                    if phase == 'train':\n",
    "                        loss.backward()\n",
    "                        optimizer.step()\n",
    "                      \n",
    "                epoch_ground_truth = epoch_ground_truth + labels.data.tolist()\n",
    "                epoch_predictions = epoch_predictions + preds.tolist()\n",
    "                running_loss += loss.item() * inputs.size(0)\n",
    "                \n",
    "            if scheduler:\n",
    "                if phase == 'train':\n",
    "                    scheduler.step()\n",
    "                \n",
    "            epoch_loss = running_loss / dataset_sizes[phase]\n",
    "            epoch_acc = accuracy_score(epoch_ground_truth, epoch_predictions)\n",
    "            epoch_weighted_acc = balanced_accuracy_score(epoch_ground_truth, epoch_predictions)\n",
    "            epoch_weighted_rec = recall_score(epoch_ground_truth, epoch_predictions, average='weighted')\n",
    "            epoch_weighted_prec = precision_score(epoch_ground_truth, epoch_predictions, average='weighted')\n",
    "            epoch_weighted_f1 = f1_score(epoch_ground_truth, epoch_predictions, average='weighted')\n",
    "            \n",
    "            learning_stats[phase].append(\n",
    "                {'loss': epoch_loss, \n",
    "                 'accuracy': epoch_acc, \n",
    "                 'weighted_accuracy': epoch_weighted_acc,\n",
    "                 'weighted_precision': epoch_weighted_prec,\n",
    "                 'weighted_recall': epoch_weighted_rec,\n",
    "                 'weighted_f1': epoch_weighted_f1})\n",
    "            \n",
    "            print('Phase Loss: {:.4f} Acc: {:.4f} Weighted Acc: {:.4f} Weighted Pre: {:.4f} Weighted Rec: {:.4f} Weighted F1: {:.4f}')\n",
    "            print(phase, epoch_loss, epoch_acc, epoch_weighted_acc, epoch_weighted_prec, epoch_weighted_rec, epoch_weighted_f1)\n",
    "     \n",
    "            \n",
    "            if (phase == 'val') and (epoch_acc > best_accuracy):\n",
    "                best_accuracy = epoch_acc\n",
    "                best_model_state_dict = copy.deepcopy(model.state_dict())\n",
    "                  \n",
    "    return best_model_state_dict, learning_stats"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eba61a90",
   "metadata": {},
   "source": [
    "**Plot classification learning curve**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "d6735c7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotLearningCurve():\n",
    "    trainLoss = [epoch[\"loss\"] for epoch in learningStats['train']]\n",
    "    valLoss = [epoch[\"loss\"] for epoch in learningStats['val']]\n",
    "    trainAcc = [epoch[\"weighted_accuracy\"] for epoch in learningStats['train']]\n",
    "    valAcc = [epoch[\"weighted_accuracy\"] for epoch in learningStats['val']]\n",
    "    numEpochs = len(learningStats['train'])\n",
    "\n",
    "    fig, (ax1, ax2) = plt.subplots(2)\n",
    "    fig.suptitle('Echinococcus classification')\n",
    "    ax1.plot(np.arange(numEpochs)+1,trainAcc,'bo-.',label=\"Training\",alpha=0.6,markersize=4)\n",
    "    ax1.plot(np.arange(numEpochs)+1,valAcc,'go-',label=\"Validation\",markersize=4)\n",
    "    ax1.axhline(y=np.max(valAcc),color=\"r\",alpha=0.4)\n",
    "    ax1.set(ylabel=\"Weighted accuracy\")\n",
    "    ax1.label_outer()\n",
    "    ax2.plot(np.arange(numEpochs)+1,trainLoss,'bo-.',label=\"Training\",alpha=0.6,markersize=4)\n",
    "    ax2.plot(np.arange(numEpochs)+1,valLoss,'go-',label=\"Validation\",markersize=4)\n",
    "    ax2.set(xlabel=\"Epoch\", ylabel=\"Loss\")\n",
    "    fig.set_size_inches(7,9)\n",
    "    plt.legend()\n",
    "    plt.savefig(os.path.join(analysis_dir_path, 'classification_results', 'classification_learning_curves.png'))\n",
    "    plt.show(block=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf9b8411",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "34f2f9a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def vizInferHeatmap_Val():\n",
    "    print(\"\\r\\n Visualizing val cases\\r\\n\")\n",
    "    for case in val_cases:\n",
    "        pathml_slide = Slide(os.path.join(pathml_slide_dir_path, case+'.pml'))\n",
    "        pathml_slide.visualizeClassifierInference('Tumor', folder=os.path.join(analysis_dir_path, 'classification_results'))\n",
    "\n",
    "def vizInferHeatmap_Test():\n",
    "    print(\"\\r\\n Visualizing test cases\\r\\n\")\n",
    "    \n",
    "    for case in test_cases:\n",
    "        pathml_slide = Slide(os.path.join(pathml_slide_dir_path, case+'.pml'))\n",
    "        pathml_slide.visualizeClassifierInference('Tumor', folder=os.path.join(analysis_dir_path, 'classification_results'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "4b7e2284",
   "metadata": {},
   "outputs": [],
   "source": [
    "def probThresholds():   \n",
    "    \n",
    "    global best_classification_threshold\n",
    "    \n",
    "    probability_thresholds = np.append(np.arange(0, 1, 0.005),[0.9925, 0.995, 0.9975, 0.999, 0.99925, 0.9995, 0.99975, 0.9999, 0.999925, 0.99995, 0.999975, 0.99999, 0.9999925, 0.999995, 0.9999975, 0.999999]).tolist()\n",
    "    threshold_accuracies_all_slides = []\n",
    "\n",
    "    for val_case in val_cases:\n",
    "        val_slide_path = os.path.join(pathml_slide_dir_path, val_case+'.pml')\n",
    "        val_slide = Slide(val_slide_path)\n",
    "\n",
    "        threshold_accuracies = val_slide.classifierMetricAtThreshold('Tumor', probability_thresholds, tileAnnotationOverlapThreshold=0.3, metric='accuracy', assignZeroToTilesWithoutAnnotationOverlap=True)\n",
    "        threshold_accuracies_all_slides.append(threshold_accuracies)\n",
    "\n",
    "    threshold_avg_accuracies = np.mean(np.array(threshold_accuracies_all_slides), axis=0)\n",
    "    index_of_best_classification_threshold = np.argmax(threshold_avg_accuracies)\n",
    "    best_classification_threshold = probability_thresholds[index_of_best_classification_threshold]\n",
    "\n",
    "    print(\"Best tile-level validation accuracy:\", threshold_avg_accuracies[index_of_best_classification_threshold])\n",
    "    print(\"Threshold that gives the best tile-level validation accuracy:\", best_classification_threshold)\n",
    "    \n",
    "    plt.figure(figsize=(7.5,6))\n",
    "    plt.plot(probability_thresholds, threshold_avg_accuracies)\n",
    "    plt.xlim(-0.05, 1.05)\n",
    "    plt.ylim(0.6, 1.0)\n",
    "    plt.title('Tile-level accuracy on val vs. tile probability threshold')\n",
    "    plt.xlabel('Probability threshold for determination\\nof number of tiles showing echinococcus')\n",
    "    plt.ylabel('Tile accuracy for echinococcus detection with\\nthresholded number of tiles')\n",
    "    plt.savefig(os.path.join(analysis_dir_path, 'classification_results', 'classification_validation_tile_accuracy_vs_probability_threshold.png'))\n",
    "    plt.show(block=False)\n",
    "    \n",
    "    %store best_classification_threshold"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9864fd2b",
   "metadata": {},
   "source": [
    "**Test set accuracy by threshold**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "0962a6dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "def testaccuracy():\n",
    "\n",
    "    test_accuracies_df_rows = []\n",
    "\n",
    "    for test_case in test_cases:\n",
    "        test_slide_path = os.path.join(pathml_slide_dir_path, test_case+'.pml')\n",
    "        test_slide = Slide(test_slide_path)\n",
    "        test_accuracy = test_slide.classifierMetricAtThreshold('Tumor', best_classification_threshold, tileAnnotationOverlapThreshold=0.3, metric='accuracy', assignZeroToTilesWithoutAnnotationOverlap=True)\n",
    "        test_accuracies_df_rows.append([test_case, test_accuracy])\n",
    "\n",
    "    test_accuracies_df = pd.DataFrame(test_accuracies_df_rows, columns=['Slide_ID', 'echinococcus accuracy (>'+str(round(best_classification_threshold, 6))+')'])\n",
    "    test_accuracies_df.to_csv(os.path.join(analysis_dir_path, 'classification_results', 'classification_test_accuracies_with_'+str(round(best_classification_threshold, 6))+'_df.csv'), index=False)\n",
    "\n",
    "    print(\"Average test set accuracy at echinococcus probability threshold of \"+str(best_classification_threshold)+\":\", test_accuracies_df['echinococcus accuracy (>'+str(round(best_classification_threshold, 6))+')'].mean())\n",
    "    print(\"Case-by-case test set accuracies:\")\n",
    "    print(test_accuracies_df)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a09f7882",
   "metadata": {},
   "source": [
    "**Slide level tile counts above threshold**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "ac96ac56",
   "metadata": {},
   "outputs": [],
   "source": [
    "def tilecountslidelevel():\n",
    "    \n",
    "    global test_tile_counts_df\n",
    "    test_tile_counts_df_rows = []\n",
    "\n",
    "    for test_case in test_cases:\n",
    "        test_slide = Slide(os.path.join(pathml_slide_dir_path, test_case+'.pml'))\n",
    "        num_metastasis_tiles_above_threshold = test_slide.numTilesAboveClassPredictionThreshold('Tumor', best_classification_threshold)\n",
    "\n",
    "        if 'echino' in test_case:\n",
    "            metastasis_ground_truth = 1\n",
    "        else:\n",
    "            metastasis_ground_truth = 0\n",
    "\n",
    "        test_tile_counts_df_rows.append([test_case, metastasis_ground_truth, num_metastasis_tiles_above_threshold])\n",
    "\n",
    "    test_tile_counts_df = pd.DataFrame(test_tile_counts_df_rows, columns=['Slide_ID', 'echinococcus_ground_truth', 'echinococcus tile count (>'+str(round(best_classification_threshold, 6))+')'])\n",
    "    test_tile_counts_df.to_csv(os.path.join(analysis_dir_path, 'classification_results', 'classification_test_tile_counts_above_'+str(round(best_classification_threshold, 6))+'_df.csv'), index=False)\n",
    "\n",
    "    print(test_tile_counts_df)\n",
    "    \n",
    "    %store test_tile_counts_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a9f0e19",
   "metadata": {},
   "source": [
    "**AUC**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "22892c73",
   "metadata": {},
   "outputs": [],
   "source": [
    "def finalAUC():\n",
    "    testAuc = roc_auc_score(test_tile_counts_df['echinococcus_ground_truth'], \n",
    "        test_tile_counts_df['echinococcus tile count (>'+str(round(best_classification_threshold, 6))+')'])\n",
    "    fpr, tpr, thresholds = roc_curve(test_tile_counts_df['echinococcus_ground_truth'], \n",
    "        test_tile_counts_df['echinococcus tile count (>'+str(round(best_classification_threshold, 6))+')'])\n",
    "\n",
    "    plt.figure()\n",
    "    plt.title('ROC with echinococcus probability threshold of '+str(round(best_classification_threshold, 6)))\n",
    "    plt.plot(fpr, tpr, 'b', label = 'AUC = %0.2f' % testAuc)\n",
    "    plt.legend(loc = 'lower right')\n",
    "    plt.plot([0, 1], [0, 1],'r--')\n",
    "    plt.xlim([0, 1])\n",
    "    plt.ylim([0, 1])\n",
    "    plt.ylabel('True Positive Rate')\n",
    "    plt.xlabel('False Positive Rate')\n",
    "    plt.savefig(os.path.join(analysis_dir_path, 'classification_results', 'classification_test_auc_roc.png'))\n",
    "    plt.show(block=False)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pathml-env",
   "language": "python",
   "name": "pathml-env"
  },
  "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.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}