MDT_ProstateX/Result analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This Notebook loads the results (list of detections) generated by the model, and computes some metrics ans visualizations using them\n",
    "\n",
    "This directory contains a copy of the Medical Detection Toolkit (https://github.com/MIC-DKFZ/medicaldetectiontoolkit/tree/torch1x) slightly adapted to the problem of prostate lesion detection. \n",
    "\n",
    "Before running this Notebook, make sure to:\n",
    " - Go to https://github.com/MIC-DKFZ/medicaldetectiontoolkit/tree/torch1x and read the instructions in order to be able to run the Toolkit. These instructions have also been cloned to this repository, in `README.md`\n",
    " - The current experiment can be found in `experiments\\exp0`.\n",
    " - In file ``experiments\\exp0\\config.py``, the paths pointing to the data might need to be updated:\n",
    "  - ``self.root_dir``\n",
    "  - ``self.pp_dir``\n",
    " - Train the model and test it: `%run exec.py --mode train_test --exp_source experiments/exp0 --exp_dir experiments/exp0`\n",
    " - Generate the final test results from the ensemble of trained models: `%run exec.py --mode analysis --exp_source experiments/exp0 --exp_dir experiments/exp0`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load `plot_lib` (https://github.com/OscarPellicer/plot_lib)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>output_subarea { flex-direction: row; flex-wrap: wrap; }\n",
       "         .output { flex-direction: row; flex-wrap: wrap; }\n",
       "         .widget-hslider { width: auto !important}</style>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from pathlib import Path\n",
    "import sys, os\n",
    "sys.path.append(os.path.join(Path.home(), 'plot_lib'))\n",
    "from plot_lib import plot, plot_multi_mask\n",
    "from IPython.display import display, HTML\n",
    "br= lambda: print(' '*150)\n",
    "CSS = \"\"\"output_subarea { flex-direction: row; flex-wrap: wrap; }\n",
    "         .output { flex-direction: row; flex-wrap: wrap; }\n",
    "         .widget-hslider { width: auto !important}\"\"\"\n",
    "HTML('<style>{}</style>'.format(CSS))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load other required libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "#Import required libs\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline\n",
    "\n",
    "import pandas as pd\n",
    "import pickle\n",
    "import numpy as np\n",
    "\n",
    "from result_computation_lib import (plot_patient, compact_class_detections, composite_score, \n",
    "                                    plot_auc, get_optimal_thresholds, match_lesions, gen_dsc)\n",
    "pd.set_option('display.float_format', lambda x: '%.5f' % x)\n",
    "pd.set_option('display.max_rows', 200)\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "SPACING= (0.5, 0.5, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are the main configuration parameters for the analysis of the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using model located at: experiments/exp0/ Subset: _test\n",
      "Using model ensembling\n"
     ]
    }
   ],
   "source": [
    "GGG_THRESHOLD= 2     #Mark lesions as significant if their GGG is >= GGG_THRESHOLD\n",
    "HOLD_OUT= True       #Should be True\n",
    "ENSEMBLE_FOLDS= True #Set to true load results from fold ensembling\n",
    "FOLD= 3              #Fold to use if ENSEMBLE_FOLDS == False\n",
    "USE_RAW= False       #Use raw predictions (before WBC)\n",
    "SUBSET= '_test'      #Must be renamed manually if any subset is to be specified here\n",
    "N_classes= 4         #GGG0, GGG1, GGG2, GGG3+\n",
    "CLASS_WEIGHTS= {0:-1, 1:-1, 2:1*12.70/17.28, 3:1, 4:1} #Model classes: BG, GGG0, GGG1, GGG2, GGG3+\n",
    "CLASS_CONFIDENCE_THRESHOLD= 0.2 #Lesions below the threshold are marked as being class 0 (NOT USED)\n",
    "USE_MODEL= 'experiments/exp0/'\n",
    "    \n",
    "print('Using model located at:', USE_MODEL, 'Subset:', SUBSET)\n",
    "if not ENSEMBLE_FOLDS:\n",
    "    print('Using Fold %d'%FOLD)\n",
    "    print(np.load('%s/fold_%d/epoch_ranking.npy'%(USE_MODEL, FOLD)))\n",
    "else:\n",
    "    print('Using model ensembling')\n",
    "    \n",
    "sys.path.append(USE_MODEL)\n",
    "from configs import configs\n",
    "cf= configs()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we load all necessary data:\n",
    " - `test_results_list`: Contains the predictions from the model\n",
    " - `patient_info`: Contains the GTs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "#This is all the data that we will need to read\n",
    "if not ENSEMBLE_FOLDS:\n",
    "    test_results_list= pd.read_pickle('%s/fold_%d/%s_pred_boxes_%slist.pickle'%(\n",
    "        USE_MODEL, FOLD, 'processed' if not USE_RAW else 'raw', 'hold_out_' if HOLD_OUT else ''))\n",
    "else: \n",
    "    test_results_list= pd.read_pickle('%s/test/%s_pred_boxes_overall_hold_out_list%s.pickle'%(\n",
    "        USE_MODEL, 'processed' if not USE_RAW else 'raw', SUBSET))\n",
    "\n",
    "#Load info of the lesions\n",
    "patient_info= pd.read_pickle(os.path.join(cf.pp_dir, 'info_df.pickle'))\n",
    "patient_info= patient_info.set_index('pid')  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we will create a dataframe `patient_df` with patient-wide results (combined as the maximum score of all individual lesions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "ignore_first_N_classes= 1 #Ignore first class (bening)\n",
    "data= []\n",
    "for pix, patient_results in enumerate(test_results_list):\n",
    "    #Extract patient results\n",
    "    pid, boxes= patient_results[1], patient_results[0]['boxes'][0]\n",
    "\n",
    "    #Get det boxes scores (& coords) sorted by score and divided by class\n",
    "    scores_all= [sorted( [[b['box_score'], b['box_coords']] \n",
    "                          for b in boxes if b['box_type'] == 'det' and b['box_pred_class_id'] == cl], \n",
    "                       key=lambda a: a[0], reverse=True)\n",
    "                            for cl in range (1 + ignore_first_N_classes, N_classes + 1)]\n",
    "    score_per_class= [scores[0][0] if len(scores) else 0. for scores in scores_all]\n",
    "    z_slice_per_class=  [(scores[0][1][4] + scores[0][1][5])/2 if len(scores) else None for scores in scores_all]\n",
    "    detected_class= np.argmax(score_per_class) + 1\n",
    "    detected_class_thresholded= detected_class if max(score_per_class) > CLASS_CONFIDENCE_THRESHOLD else 0\n",
    "    \n",
    "    #The score depends on the GGG_THRESHOLD\n",
    "    score= max(score_per_class[GGG_THRESHOLD - ignore_first_N_classes:])\n",
    "    \n",
    "    #Find the z-axis to which the max-scoring lesions from each category belong\n",
    "    \n",
    "    #Add to list\n",
    "    data.append( [pid, score, detected_class, detected_class_thresholded, *score_per_class, *z_slice_per_class] )\n",
    "    \n",
    "patient_df= pd.DataFrame(data=data, columns=['pid', 'pred_score', 'det_class', 'det_class_thresholded', \n",
    "                                     'GGG1_score', 'GGG2_score', 'GGG3+_score', 'z1', 'z2', 'z3+']).set_index('pid')\n",
    "#patient_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot patients"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load processed results and plot them patient by patient.\n",
    "\n",
    "Only the highest-scoring prediction from highly overlapped predictions is being plotted. Change `match_iou` to 1. in `compact_class_detections` to plot all predictions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0: ProstateX-0001 | Lesion GGGs: [1]\n",
      "Max scores: GGG1: 0.19 | GGG2: 0.24 | GGG3+: 0.16\n",
      "                                                                                                                                                      \n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7ae7bd7bc07b4379b01c1135bff68115",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=0, continuous_update=False, description='c/t', max=7, style=SliderStyle(…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ba3938170e5a414cb6393f2ac0b9354f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=0, continuous_update=False, description='c/t', max=7, style=SliderStyle(…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1: ProstateX-0003 | Lesion GGGs: [10, 10]\n",
      "Max scores: GGG1: 0.04 | GGG2: 0.02 | GGG3+: 0.00\n",
      "                                                                                                                                                      \n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bf2590b80ba840718d1f62888d2ba5ef",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=0, continuous_update=False, description='c/t', max=7, style=SliderStyle(…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b3ef388218e043fdb055414fe6fbfe99",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=0, continuous_update=False, description='c/t', max=7, style=SliderStyle(…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2: ProstateX-0010 | Lesion GGGs: [10]\n",
      "Max scores: GGG1: 0.11 | GGG2: 0.14 | GGG3+: 0.05\n",
      "                                                                                                                                                      \n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9c38643b07c5426e9ba54cf9bbc54479",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=0, continuous_update=False, description='c/t', max=7, style=SliderStyle(…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4501c739868c47eb808db56987842bb7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=0, continuous_update=False, description='c/t', max=7, style=SliderStyle(…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3: ProstateX-0018 | Lesion GGGs: [2]\n",
      "Max scores: GGG1: 0.32 | GGG2: 0.20 | GGG3+: 0.31\n",
      "                                                                                                                                                      \n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f10c5e9e001f4af484b3d4cd7279b6b4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=0, continuous_update=False, description='c/t', max=7, style=SliderStyle(…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a28088084a994c8d80694fb1c510fe8d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=0, continuous_update=False, description='c/t', max=7, style=SliderStyle(…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4: ProstateX-0022 | Lesion GGGs: [10]\n",
      "Max scores: GGG1: 0.02 | GGG2: 0.04 | GGG3+: 0.01\n",
      "                                                                                                                                                      \n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a07404bb7d4f44f1b95181c081ef2819",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=0, continuous_update=False, description='c/t', max=7, style=SliderStyle(…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "deac82f449744cd28453132a7e3bcf6a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=0, continuous_update=False, description='c/t', max=7, style=SliderStyle(…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Keep only non-ovelapping (mostly) predictions. Set match_iou=1 to disable it\n",
    "#Default: benign_class=1 (i.e. class 0)\n",
    "test_results_list_compacted= compact_class_detections(test_results_list, match_iou=0.25, \n",
    "                                                      benign_class=1, class_weights=CLASS_WEIGHTS)\n",
    "\n",
    "for pix, patient_results in enumerate(test_results_list_compacted):    \n",
    "    #Get and print patient data:\n",
    "    pid= patient_results[1]\n",
    "    \n",
    "    #Stop after the fifith\n",
    "    if pix == 5: break\n",
    "    \n",
    "    #Print some info\n",
    "    pat_df= patient_info.loc[pid, :]\n",
    "    print('%d: %s | Lesion GGGs: %s'%\n",
    "          (pix, pid, pat_df.class_target))\n",
    "    print('Max scores: GGG1: %.2f | GGG2: %.2f | GGG3+: %.2f'%\n",
    "            tuple(patient_df.loc[pid, ['GGG1_score', 'GGG2_score', 'GGG3+_score']].values.tolist()) )\n",
    "    br()\n",
    "    \n",
    "    #Show plot at highest scoring lesion\n",
    "    z= np.round(patient_df.loc[pid, ['z0', 'z1', 'z2', 'z3+'][patient_df.loc[pid, 'det_class']]])\n",
    "    z= 12 if np.isnan(z) else z\n",
    "    \n",
    "    #Plot\n",
    "    plot_patient(patient_results, cf.pp_test_data_path, N_classes= N_classes, seg_threshold=0.5, \n",
    "                 min_threshold= 0.08, plot_max_per_class=100, spacing=SPACING, scale='auto', ct=0, z=z, \n",
    "                 dpi=100, text_kwargs={'linespacing':0.9},\n",
    "             class_colors={i:plt.get_cmap('Paired')(ci) for i, ci in enumerate(np.arange(0, 1, step=1/(N_classes + 2)))})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prostatex-like lesion-level results\n",
    "\n",
    "Detections must be matched with ground truth lesions to get the lesion-level results. They can be matched either by looking at the intersection of their bounding boxes, or by looking at their distance (``USE_DISTANCE= True``). In both cases, the highest-score non-benign lesion from the matched lesions is considered as the score for that lesion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "USE_DISTANCE= True #Use distance instead of IoU for matching lesions\n",
    "DISTANCE_THRESHOLD= 15 #Consider lesions within a 15mm radius of the GT lesion\n",
    "threshold_list= [] #[0.4, 0.2, 0.1, 0.05, 0.025]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lesion-level (prostatex-like). Sig if GGG >= 2\n",
      "  ProstateX\n",
      "    N= 69    AUC  | Acc.   | Sens.  | Specif.\n",
      "    t=0.38: 0.890 |  0.870 |  0.308 |  1.000\n",
      "    t=0.09: 0.890 |  0.855 |  0.923 |  0.839\n",
      "    t=0.09: 0.890 |  0.855 |  0.923 |  0.839\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEWCAYAAACT7WsrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de7xUVd3H8c9XQFHAK5YoCqiACCiipvikwOMlTMOn1BS11Ex9HlPzMbXULqSlhpfMUoMsvJSoYRbh/YLYk1Kg4A3FUEhJkouAoCAIv+ePvc9hGOac2YfhzJwz5/t+vebF7L3XWnuddQ7zm7XW3msrIjAzMytkk0pXwMzMmi4HCTMzq5ODhJmZ1clBwszM6uQgYWZmdXKQMDOzOjlI2HokXSbptkYqe7akwxqj7Jxz3C7pR415DrOWwkGiCmzsD96IuCoivr6xymtKJJ0mabWkZZI+kPSipKPz0mwm6WpJb0taLukfki6WpLx0n5P0jKSlkuZLmihpaD3n7iHp95IWSFoi6SVJF0pq1Vg/r1mpHCSsJXouItoDWwO3APdI2jrn+O+BQ4HPAx2ArwBnAT+rSSDpuDTdnUBn4NPA94EvFDqhpN2AvwHvAH0jYivgeGC/9BwNIql1Q/OYbQgHiSon6WhJ0yQtlvSspL1yjn1b0r/Sb8IzJB2a7h8u6bc56YZKejUt42lJvXKOzZZ0UfqteImkeyW1zVi3TSR9R9KbkhZKuk/StumxRySdm5f+RUlfSt/vIelxSe+ndf9yQ9smItYAdwHtgO5puYcCRwDHRsQrEfFJREwCTgG+IWn3tEdxA3BlRNwWEUsiYk1ETIyIM+s43Q+BZyPiwoiYm55/RkScFBGLJQ2SNCfv563tIaa/k7GSfivpA+CytJezbU76fdJeSpt0+2uSXpO0SNKjkro0tI3MHCSqmKT+wG+As4HtgJHAuHQ4pSdwLrB/RHQAPgfMLlBGD2AMcAGwPfAQ8GdJm+Yk+zIwBOgG7AWclrGK5wP/BQwEdgQWATenx+4GhuXUY0+gC/CgpHbA42maT6XpbpHUO+N5a8psBZwOrAL+me4+HPhbRLyTmzYi/gbMIelh9AR2BsY24HSHNTB9IcekZWwNXAs8Bxybc/wkYGxErJL0X8BlwJdIfm9/Ifk9mjWIg0R1OxMYGRF/i4jVEXEH8DFwILAa2AzYU1KbiJgdEW8WKOME4MGIeDwiVgHXAZsDB+WkuSki3o2I94E/A/0y1u9s4PKImBMRHwPDgePSoZQHgH45335PBv6QpjsamB0Ro9Nv+i8A9wPHZTzvgZIWAyvSn+eUiJiXHusIzK0j39z0+HY521lt18D0hTwXEX9Mey3LyQmkae/mxHQfJG17dUS8FhGfAFexbnuaZeIgUd26AN9Kh4kWpx+MOwM7RsRMkt7BcGCepHsk7VigjB1Z+y27ZojmHWCnnDT/znn/EdAeQNLD6QTxMkkn11G/B3Lq9hpJ8Pp0RCwFHiT54CP993c5+Q7I+7lOBnbI2C6TImJrYBtgHHBwzrEFQKc68nVKjy/M2c5qYQPTF/JO3vZYYED6ezsECJIeAyRt9LOc9nkfEOv+3syKcpCobu8AP46IrXNeW0TEGICIuDsiPkvygRLATwqU8W56HKj9xroz8K9iJ4+IIyOiffr6XYEk7wBH5tWvbUTUlD0GGCZpAEnvZUJOvol5+dpHxP9kaZSc+i0DzgG+ImmfdPcTJAFo59y0kj6T/txPATPSOuQO9RTzRJH0HwJb5JyvFckw0TpVzqv/YuAxkuG+k4AxsXZZ53eAs/PaaPOIeLYBdTZzkKgibSS1zXm1Bn4F/LekA5RoJ+koSR0k9ZT0n5I2Ixl2WU7yLT7ffcBRkg5NJ0S/RTJktTE+bH4J/LhmCETS9pKOyTn+EEmAugK4N+3FAIwHekj6iqQ26Wv/3An1rCJiIXAbyZVJRMQTwJPA/ZJ6S2ol6UCSXsytEfGP9IP4QuB7kk6XtGU6Cf9ZSaPqONUPgIMkXStph/Tn3T2diN4aeANom/5+2gDfJRkOLOZu4KskAejunP2/BC6tmaeRtJWk4xvSNmbgIFFNHiL5oK95DY+IKSTzEr8gmRSeydpJ5c2Aa0iGT/5NMgF8WX6hETGD5Mqen6dpvwB8ISJWboQ6/4xkuOcxSUuBScABOef+GPgDyaTv3Tn7l5JcgXQiSU/n3yS9oCwfqoXcCHxea6/8Opak1/IIsAz4LfBr4LycOowlma/5WlqH94AfAX8qdIJ0vmcA0BV4VdISknmUKcDSiFhC0qu5jaSX9iHJRHkx40iuzHovIl7MOd8DJG1yT3o11CvAkRnKM1uH/NAhMzOri3sSZmZWJwcJMzOrk4OEmZnVyUHCzMzqVOoiYZ71NjNrIEmPRsSQStcjC/ckzMzKr2N9ByX9RtI8Sa/UcVySbpI0U8nimv0bp5oOEmZmTdHtJItm1uVIkvtjupMsY39rY1XEQcLMrImJiGdI1tuqyzHAnZGYBGwtqdS1wQpykDAza352Yt0FH+fQSIs3+ulWZi3N8K3Wvu8xBE66F+4+Ad54JCfNEpgyGsZfsHbfsHugUz+4YY+1+/qfCkNvgpGHwNx0VZD2O8BFM2DC1TDxmrVpz3o6+XfUoLX7Bn4HBl8K1/WEZeliwp32hrOfgXHnwwt3rE174eswdxqMOXHtvqNvhP1Ob5o/05TRSd0K6yhpSs72qIioa92vQlRgX6NcSFTqshy+usmsOan/g8vKRNLzEbFfkTRdgfER0afAsZHA0zUrOkuaAQyqeerhxuThJrOWJPdbtDVn44Cvplc5HQgsaYwAAR5uMjNrciSNAQaRDEvNIVlqvg1ARPySZNXnz5Os7PwRyWN4G4WDhJlZExMRw4ocD+Ab5aiLh5vMWpJh91S6BtbMOEiYtSSd+lW6BtbMOEiYtSS5l3qaZeAgYWZmdXKQMDOzOjlImLUk/U+tdA2smfEd12ZmZZbljuumwj0Js5Zk5CGVroE1Mw4SZi1JzYJ1Zhk5SJiZWZ0cJMxakvY7VLoG1sx47aZ813SBFYuT914/3z9Ttf1MF83ArCF8dVM+r7dvZo2sOV3d5CBhZlZmzSlIeE4iX+7wgJlZC+cgYWZmdXKQMDOzOjlI5OsxpNI1MDNrMjxxbWZWZp64bs7uPqHSNTAzazIcJPLl3uRkZtbCOUiYmVmdHCTMzKxOnrg2MyszT1w3Z1NGV7oGZmZNhoNEvtyVOM3MmjlJ+0v6maQXJM2V9JakcZLOltShWH4vFW5mVqUkjQcWAn8CrgfmAW2BHsBg4EFJIyJifF1lOEiYmVWvMyLivbx9K4C/p6+fSPpUfQV4uCnfsHsqXQMzs40iN0BI6ixpcPp+M0nt0jTz6ivDQSJfp36VroGZ2UYl6WvAOOC2dFcXkiGoohwk8uU+LtLMrDqcDxwIfAAQEW8A9Q4z1XCQMDOrfisiYmXNhqRWgLJkdJAwM6t+f5V0CdA2nZe4F6jziqZcDhL5+p9a6RqYmW1slwBLgdeBbwJPApdlyehlOczMyizLshyShgA/A1oBt0XENXnHdwHuALZO03wnIh6qo6xzI+IXxfYV4p5EvpGHVLoGZtbCpXMGNwNHAnsCwyTtmZfsu8B9EbEPcCJwSz1Ffq3AvjOy1MU30+Wb+2Kla2Bm9hlgZkS8BSDpHuAYYHpOmgC2TN9vBbybX4ikE0gCSDdJf8g51AFYnKUiDhJmZuXXUdKUnO1RETEqZ3sn4J2c7TnAAXllDAcek3Qe0A44rMB5/k6yLEdnkp5JjaXA1CwVdZDI136HStfAzKrfgiJzEoUuT82fAx4G3B4R10saANwlqU9ErKnNEDELmAU8saEV9ZxEvotmVLoGZmZzgJ1ztjuz/nDSGcB9ABHxHMnCfR0LFZauBDtJ0hJJKyR9LOmDLBVxkMg34epK18DMbDLQXVI3SZuSzCuMy0vzNnAogKReJEFifh3l3QKcCrxFMh9xLnBjloo4SOSbeE3xNGZmjSgiPiH5IH8UeI3kKqZXJV0haWia7FvAmZJeBMYAp0Xd9zRsEhEzgNYRsSoifkXhOYz1eE7CzKwJSu95eChv3/dz3k8H/iNjcR+mPZIXJV0FzAXaZ8nonoSZWfU7jeTz/lxgNdAdOC5LxpZ9x/WEq9cdXjrr6eTfHfepRG3MrIXIcsd1GepwQET8rWi6FhskruvpK5nMrCLKFSQkbQIcS3LfxaMR8Vq63MdlwDYR0bdYGS13TmLZvytdAzOzxnYbsCvJ1VK3SvoHMAi4NCLGZimg5QYJM7PqdwCwV0SslrQ5sADYPSLmZi2g5U5cd9q70jUwM2tsH0fEaoCIWA7MaEiAgJY8J2FmViFlnJP4iOQZEpAs9dEz3RYQEdG/WBktd7hp3Pkw9KZK18LMrDEVnZgupuUGiRfucJAws6oWEW+WWkbLnZMwM7OiHCTMzKxOLTdIXPh68TRmZlVC0qaSdm9ovpYbJOZOq3QNzMzKQtJRwMvA4+l2P0kPZMnbcoPEmBMrXQMzs3K5guTGusUAETENyNSraLlBwsys5VgVEYvz9mW6z63lXgJrZtZyvCbpy8AmkroB3wQmZcnYcnsSR2d6cp+ZWTU4F9gXWAP8AVgBXJAlY/NbluOaLrBicXJ10txp684tHH0j7Hc6DN9q7b4eQ+Cke+HuE+CNR9buH76kfHU2M8tR7udJSOobES9vUN5mFySGb+UPeDNr1ioQJP4CbAvcB9wbEZnvAWh+w039T610DczMmpWIOBg4AlgK3CFpqqTvZMnb/HoSZmbNXCUfXyqpF3ApMCwi2hRL3/x6EiMPqXQNzMyaFUndJX1X0jTgVyRPqtslS97mdwns3BcrXQMzs+bmbuAeYGhEvN2QjM0vSJiZWYNExP4bmrf5BYn2O1S6BmZmzYKkMRExTNJU1p1DzvxkOk9cm5mVWRkfX9o5IuZI2q3Q8SwPJWp+E9cTrq50DczMmoWImJO+PSMi3sx9AWdkKaP5BYmJ11S6BmZmzc2QAvuOypKx+c1JmJlZJpLOBv4b6CHphZxDHYDns5ThIGFmVr3uA54ErgZy77BeGhHzshTQ/Cau350KO+5T9tOamW0sZZy4bhcRH0rastDxiPigWBnuSZiZVa+xwJHAqyRf6pVzLMhw13Xz60l4FVgza+YquXZTQzW/q5vMzKxBJB0oaYv0/TBJIyR1zpK36QWJ35+e/DvykLTXsBVc1zPZN+FqaLt15epmZlYmkoZImiFpZl3Lekv6sqTpkl6VdHc9xY0ClkvaC7gMeA/4XaZ6NLnhJg8nmVmVKzbcJKkV8AZwODCHZNXWYRExPSdNd5Krl/4zIhZJ+lRdVyxJeiEi+kv6HjA3Im6r2Vesrk2vJ2FmZp8BZkbEWxGxkmQF12Py0pwJ3BwRiwCKXNL6oaSLga8AD0raBCj6LAlwkDAzq4SOkqbkvM7KO74T8E7O9px0X64eJDfJ/VXSJEmF7qqucQLJlU1nR8RcoDNwQ5aKNr1LYC/M/OhVM7PmakGRq5tUYF/+8H5roDswiORD/y+S+kTE4vUyRrwr6TfAfmkwmRIRo7NUtOn1JOZOq3QNzMwqbQ6wc852Z+DdAmn+FBGrImIWMIMkaKxH0rHACyTDTV8Fpkj6YpaKNL0gMebEStfAzKzSJgPdJXWTtClwIjAuL80fgcEAkjqSDD+9VUd53wf2j4iTI+Ik4ABgeJaKNL0gYWbWwkXEJ8C5wKPAa8B9EfGqpCskDU2TPQoslDQdmABcHBEL6yhyk4h4L2d7Phk//30JrJlZmZX7jmtJNwB7kDzrGpKeyesRcVGxvE1v4vroGytdAzOzavMt4HjgsyST4neQrOtUVNPrSZiZVbly9iQkHQ3sDrwcEU82NH/Tm5MYvlWla2BmVhUk/ZzkORI7ASMkXdbQMprecJOZmW0sg4F+EfGJpHbAROCqhhTQ9HoSZma2saxMr5QiIj5kAz7zm15Pokd9d5abmVkD7JHzbGsBPdNtAZFlgT9PXJuZlVkZH1+6W33HI+LNYmWUtycx4+F176g++kbY7/R1J6u7/Aec/lBZq2VmVo2yBIFiytuT+GAubNmplPOZmTV7fnxpXW7Yo6ynMzOz0vjqJjOzFkDSppJ2b2g+Bwkzsyon6SjgZeDxdLufpAey5C1vkOh/allPZ2ZmAFxBsjz4YoCImEayVEdR5Q0SQ28q6+nMzAyAVQWeWJfpwqPyBomRh5T1dGZmBsBrkr4MbJI+yOhGYFKWjOUNEnNfLOvpzMwMSB5gtC+wBngA+Bi4IEvGprcsh5mZbVTpuk3fTl8NUt4g0X6Hsp7OzMxA0uMUmIOIiCOK5S1vkLhoRllPZ2ZmAHw3531b4FiSIaeiyhskJlwNgy8t6ynNzFq6iPhb3q6JkiZmyVveieuJ15T1dGZmBpK2zHltLelQINNCep64NjOrfq+SzEkI+ASYBZyZJaODhJlZFZO0CXB8RGS6LyJfeYebznq6rKczM2vpImINcOOG5vcCf2Zm1e9xScdsSMbyBolRg8p6OjMzA5I7rh+QtFzS+5IWSXo/S0bPSZiZVSlJu0TE20DHDS3Dw01mZtXrjwARsbrQK0sB5Q0SA79T1tOZmbVwKrmAiExLitelpMxmZi2RpOcjYr8ynGcecE9dxyPi/GJllHdO4rqeXr/JzKx8lgPPl1JAeYPEsn+X9XRmZi3cwoi4o5QCPHFtZla9VpZaQHmDRKe9y3o6M7PmStIQSTMkzZRU51U/ko6TFJIKzXGcWOQcktS53jSeuDYzK69iE9eSWgFvAIcDc4DJwLCImJ6XrgPwILApcG5ETMk7/nuSzsCfSOYm5pM8T2J3YDBwKPCDiHi8rrqUtycxruhEupmZwWeAmRHxVkSsJLlCqdCyGlcCI4AVhQqJiOOB7wE9gZuBv5AEjK8DM4D/rC9AQLknrl+4A4beVNZTmpk1QzsB7+RszwEOyE0gaR9g54gYL+miugpKex+Xb2hFvCyHmVn5dZSUOzQ0KiJG5WwXugmudng/Xf77p8BpjVO9tRwkzMzKb0GRm+nmADvnbHcG3s3Z7gD0AZ6WBLADME7S0Px5iVKVd07iwtfLejozs2ZqMtBdUjdJm5JcpTSu5mBELImIjhHRNSK6ApOAjR4goLGCxPCtkn+njE7e17zmTmuU05mZVZOI+IRkee9HgdeA+yLiVUlXSBra0PIk3S/pqHSYqmF5G+US2OFbwfAlpZRrZla1yrV2U875DgNOBw4Efg/cHhGZhnZ8x7WZWZWLiCci4mSgPzCb5El1z0o6XVKb+vI2TpDoMaRRijUzsw0jaTuSq6G+DkwFfkYSNOq9T8J3XJuZlVkFhpv+AOwB3EUy1DQ359iU+urSOD2Ju09olGLNzGyD3BYRe0bE1TUBQtJmAMWCVeMEiTceaZRizcxsg/yowL7nsmT0zXRmZlVK0g4kS3xsni7jUXMn95bAFlnKcJAwM6tenyOZrO4M3JCzfylwWZYCPHFtZlZmFZi4PjYi7t+QvI3Tk5gyGvY7vVGKNjOzbCSdEhG/BbpKujD/eETcUCDbOhpn4nr8BY1SrJmZNUi79N/2JIsC5r+K8pyEmVmVioiR6dtbImL+hpThZTnMzKrfs5Iek3SGpG0akrFxgsSwexqlWDMza7iI6A58F+gNPC9pvKRTsuRtnCDRqV+jFGtmZhsmIv4eEReSPD/7feCOLPlKm5O4pgts0wXOfgbGnZ88w7qGlwo3M2sSJG0JfJHk4UW7AQ+QBIvieUu6T2L4VuFgYGbWMBW4T2IW8EeShxdlWo6jRmk9ifY7lJTdzMzKYtfYwB6B77g2MyuzcvUkJN0YERdI+jMFPq8jouijUEvrSUy4GgZfWlIRZmbWaO5K/71uQwso7eqmideUlN3MzBpPRDyfvu0XERNzX0Cmy1B9M52ZWfU7tcC+07Jk9LIcZmZVStIw4CSgm6RxOYc6AAuzlFFakDjr6ZKym5lZo3oWmAt0BK7P2b8UeClLAe5JmJlVqYj4J/BPYMCGllHanMSoQSVlNzOzxiPp/9J/l0r6IOe1VNIHWcpwT8LMrEpFxGfTfzM9O6IQX91kZlblJO0mabP0/SBJ50vaOkve0oLEwO+UlN3MzMrifmC1pN2BXwPdgLuzZCwtSPhuazOz5mBNRHxCshLsjRHxv0CnLBlLCxLX9Swpu5mZlcWq9J6JU4Hx6b42WTKWFiSW/buk7GZmVhank1wG++OImCWpG/DbLBl9dZOZWZWLiOnA+Tnbs4BMi++VFiQ67V1SdjMza3yS/gMYDnQh+dwXEBGxa9G8fp6EmVl5VeDJdK8D/ws8D6yu2R8RRddvKm1OYtz5xdOYmVmDSRoiaYakmZLWu99A0oWSpkt6SdKTkrrUU9ySiHg4IuZFxMKaV5Z6lBYkXrijpOxmZrY+Sa2Am4EjgT2BYZL2zEs2FdgvIvYCxgIj6ilygqRrJQ2Q1L/mlaUunrg2M2t6PgPMjIi3ACTdAxwDTK9JEBETctJPAk6pp7wD0n9zh7gC+M9iFXGQMDMrv46SpuRsj4qIUTnbOwHv5GzPYe0HfSFnAA/XdTAiBm9QLSk1SFz4eknZzcxaqAVFJq5VYF/BC4UknULSQxhYZ2HSp4GrgB0j4sh06GpARPy6WEVLm5OYO62k7GZmVtAcYOec7c7Au/mJJB0GXA4MjYiP6ynvduBRYMd0+w3ggiwVKS1IjDmxpOxmZlbQZKC7pG6SNgVOBHIfP4qkfYCRJAFiXpHyOkbEfcAagHQdp9X1Z0l4qXAzsyYm/RA/l+Tb/2vAfRHxqqQrJA1Nk10LtAd+L2la3jOs830oaTvSIStJBwJLstTFE9dmZk1QRDwEPJS37/s57w9rQHEXkvREdpP0V2B74LgsGUsLEkffWFJ2MzNrfBHxgqSBQE+SSfEZEbEqS97Shpv2O72k7GZm1ngk7S9pB6gdwtoX+DFwvaRts5RRWpAYvlVJ2c3MrFGNBFYCSDqEZOXXO0nmI0bVk6+W5yTMzKpXq4h4P31/AslNe/cD90vKdA+Dr24yM6terSTVdAYOBZ7KOZapk1BaT6LHkJKym5lZoxoDTJS0AFgO/AVA0u5kvATWz5MwMyuzcj5PIr0nohPwWER8mO7rAbSPiBeK5S9tuOnuE0rKbmZmjSsiJkXEAzUBIt33RpYAAaUGiTceKSm7mZk1bZ64NjOzOjlImJlZnTxxbWZWZuWcuC5VaT2JKaM3UjXMzKwpKi1IjM/0zAozM2umvCyHVa1Vq1YxZ84cVqxYUemqWAvVtm1bOnfuTJs2bSpdlQ3mIGFVa86cOXTo0IGuXbsiFXpksFnjiQgWLlzInDlz6NatW6Wrs8FKG24ads9GqobZxrdixQq22247BwirCElst912zb4nW1qQ6NRvI1XDrHE4QFglVcPfX2lB4oY9NlI1zMysKfLNdGaNpFWrVvTr148+ffpw/PHH89FHH5Vc5u23386777670dINHTqUu+66q3b7zDPP5Nprry2a78Ybb+TOO++s3f7kk0/o2LEjl1566TrpunbtyoIFC2q3n376aY4++uja7Ycffpj99tuPXr16sccee3DRRRcVPXcxzz//PH379mX33Xfn/PPPp9C9YEuWLOELX/gCe++9N71792b06LWX83/729+mT58+9OnTh3vvvbd2/6xZszjggAPo3r07J5xwAitXrgTgF7/4xTr5q05EbPjrB1uGWVM1ffr0ip6/Xbt2te9POumkuP7669c5vmbNmli9enWDyhw4cGBMnjx5o6WbNWtWdO3aNRYtWhR//etfo2/fvrFy5cp686xatSr69u0bq1atqt334IMPxkEHHRS77rprrFmzpnZ/ly5dYv78+bXbEyZMiKOOOioiIl5++eXYdddd47XXXqst9+abby5a52L233//ePbZZ2PNmjUxZMiQeOihh9ZL8+Mf/zguueSSiIiYN29ebLPNNvHxxx/H+PHj47DDDotVq1bFsmXLYt99940lS5ZERMTxxx8fY8aMiYiIs88+O2655ZaIiPjwww+jX79+ddan0N8hMCVK+ewt46u0nkT/UzdOpDKrcgcffDAzZ85k9uzZ9OrVi3POOYf+/fvzzjvvMGbMGPr27UufPn349re/DcDq1as57bTT6NOnD3379uWnP/0pY8eOZcqUKZx88sn069eP5cuXc8UVV7D//vvTp08fzjrrLCKiYLrnn3+egQMHsu+++/K5z32OuXPnAsk3/bPOOotLLrmEc845h1/84hdFL9d86qmn6N+/P61br704csyYMXzzm99kl112YdKkSZnaZMSIEVx++eXssUcybN26dWvOOeecDWneWnPnzuWDDz5gwIABSOKrX/0qf/zjH9dLJ4mlS5cSESxbtoxtt92W1q1bM336dAYOHEjr1q1p164de++9N4888ggRwVNPPcVxxx0HwKmnnlpb7hZbbEHXrl35+9//XlLdm6wSo4xZk7XeN7inror4wZZrX/96IXnl7nvqqiTttT3W7vvlwcm+P523btol79Z7/pqexKpVq2Lo0KFxyy23xKxZs0JSPPfccxER8a9//St23nnnmDdvXqxatSoGDx4cDzzwQEyZMiUOO+yw2rIWLVoUEev3EBYuXFj7/pRTTolx48atl27lypUxYMCAmDdvXkRE3HPPPXH66afX5lu5cmXsvPPOcdJJJ2Vq1+9///tx00031W5/9NFH0alTp/jwww9j5MiRcd5559Ueq68nsc8++8S0adOKnu+pp56Kvffee73XgAED1ks7efLkOPTQQ2u3n3nmmdrz5frggw9i0KBBscMOO0S7du1i/PjxERHx6KOPxkEHHRQffvhhzJ8/P7p16xbXXXddzJ8/P3bbbbfa/G+//Xb07t27dvtHP/pRXHfddQXr39x7EqXdJzHyEDj7mY0Trcwa2+BLk1e+4QUe0HXRjPX3Db0peWW0fPly+vVLrgA8+OCDOeOMM3j33Xfp0qULBx54IACTJ09m0KBBbL/99gCcfPLJPPPMM3zve9/jrbfe4rzzzuOoo47iiCOOKHiOCRMmMGLECD766CPef/99evfuzVjgsBEAAA+qSURBVBe+8IV10syYMYNXXnmFww8/HEh6KZ06dao9/tJLLxERvP7666xZs4ZNNql/gGHu3Ln06tWrdnv8+PEMHjyYLbbYgmOPPZYrr7ySn/70p7Rq1arg1T0NveJn8ODBTJuW6XHMBecfCp3v0UcfpV+/fjz11FO8+eabHH744Rx88MEcccQRTJ48mYMOOojtt9+eAQMG0Lp166LlfupTn+L1119vwE/VfJQ23DT3xY1UDbPqs/nmmzNt2jSmTZvGz3/+czbddFMA2rVrV5um0IcPwDbbbMOLL77IoEGDuPnmm/n617++XpoVK1ZwzjnnMHbsWF5++WXOPPPMgtfkRwS9e/eurcvLL7/MY489BsCaNWs455xzuOuuu+jevTu33nprpp8r9zxjxozhiSeeoGvXruy7774sXLiQCRMmALDddtuxaNGi2rTvv/8+HTt2BKB37948//zzRc83YcIE+vXrt97roIMOWi9t586dmTNnTu32nDlz2HHHHddLN3r0aL70pS8hid13351u3brVfshffvnlTJs2jccff5yIoHv37nTs2JHFixfzySefFCx3xYoVbL755kV/lubIVzeZVdABBxzAxIkTWbBgAatXr2bMmDEMHDiQBQsWsGbNmtpv5i+8kDxErEOHDixduhSg9oO6Y8eOLFu2jLFjx9aWm5uuZ8+ezJ8/n+eeew5Ilit59dVXARg5ciTdu3dn0KBB3HDDDYwYMYL58+fXW+devXoxc+ZMAD744AP+7//+j7fffpvZs2cze/Zsbr75ZsaMGQPAoEGDaq+eWr16Nb/97W8ZPHgwABdffDFXXXUVb7zxBpAErBtuuGG989X0JPJfzz777HppO3XqRIcOHZg0aRIRwZ133skxxxyzXrpddtmFJ598EoD33nuPGTNmsOuuu7J69WoWLlwIJD2sl156iSOOOAJJDB48uLaN77jjjnXKfeONN+jTp0+97dZslTRedW2PgmNwZk1BU7q6qcasWbPWGcuOiPjd734Xffr0id69e8fFF18cERHTpk2LffbZp3b8veYKnbFjx0aPHj1i7733jo8++iguv/zy2G233eLQQw+N0047LX7wgx8UTDd16tQ4+OCDY6+99oo999wzRo0aFe+991506dIl5s6dW1uX66+/Pk477bR6f67Zs2fHwQcn8zSjR4+OE044YZ3jCxcujI4dO8aKFSti8eLFMWzYsNhrr72ib9++cfHFF69zRdef//zn6N+/f+yxxx7Rq1evuOiiizK2bt0mT54cvXv3jl133TW+8Y1v1F5tdeutt8att94aEclc0OGHH17b7nfddVdERCxfvjx69eoVvXr1igMOOCCmTp1aW+6bb74Z+++/f+y2225x3HHHxYoVK2qP7bPPPuvMveRq7nMSfp6EVa3XXnttnbFz23i++MUvMmLECLp3717pqlTc1KlTueGGG9a53yRXob/DlvM8iQlXb6RqmFlzcs0119ReRtvSLViwgCuvvLLS1Wg0pQWJiddspGqYWXPSs2dPDjnkkEpXo0k4/PDD6dq1a6Wr0Wg8cW1VrcThVLOSVMPfn4OEVa22bduycOHCqviPas1PRPI8ibZt21a6KiUp7Wa6s57eKJUwaww118wXu6TTrLHUPJmuOfOT6axqtWnTplk/EcysKShtuGnUoI1TCzMzW4ekIZJmSJop6TsFjm8m6d70+N8kdW2MenhOwsysiZHUCrgZOBLYExgmac+8ZGcAiyJid+CnwE8aoy4OEmZmTc9ngJkR8VZErATuAfLXFzkGuCN9PxY4VI3wvNSS5iT0ww9e5Ydq3k/53jg6AguKpqpuboOE2yHhdqi/DXpKmpKzPSoiRuVs7wS8k7M9Bzggr4zaNBHxiaQlwHb1nHODlDpxvaK53FremCRNaent4DZIuB0SboeS26BQjyD/Wu4saUrm4SYzs6ZnDrBzznZnIP+h5bVpJLUGtgLe39gVcZAwM2t6JgPdJXWTtClwIjAuL804oOYZ0scBT0Uj3Dla6nDTqOJJWgS3g9ughtsh4XYooQ3SOYZzgUeBVsBvIuJVSVeQLDM+Dvg1cJekmSQ9iBM3RqXzlbpUuJmZVTEPN5mZWZ0cJMzMrE6ZgkRTuT28kjK0wYWSpkt6SdKTkrpUop6NrVg75KQ7TlJIqsrLILO0g6Qvp38Tr0q6u9x1bGwZ/k/sImmCpKnp/4vPV6KejUnSbyTNk/RKHccl6aa0jV6S1L/cdSxZseebkkyavAnsCmwKvAjsmZfmHOCX6fsTgXsr/VzWjfnK2AaDgS3S9/9TbW2QtR3SdB2AZ4BJwH6VrneF/h66A1OBbdLtT1W63hVog1HA/6Tv9wRmV7rejdAOhwD9gVfqOP554GGSexoOBP5W6To39JWlJ9Fkbg+voKJtEBETIuKjdHMSyXXN1SbL3wLAlcAIoFrvxs/SDmcCN0fEIoCImFfmOja2LG0QwJbp+61Y/zr/Zi8inqH+exOOAe6MxCRga0mdylO7jSNLkCh0e/hOdaWJiE+AmtvDq0WWNsh1Bsm3h2pTtB0k7QPsHBHjy1mxMsvy99AD6CHpr5ImSRpSttqVR5Y2GA6cImkO8BBwXnmq1qQ09LOjyclyn0STuT28gjL/fJJOAfYDBjZqjSqj3naQtAnJapSnlatCFZLl76E1yZDTIJJe5V8k9YmIxY1ct3LJ0gbDgNsj4npJA0iu6e8TEWsav3pNRrP/bMzSk2gyt4dXUJY2QNJhwOXA0Ij4uEx1K6di7dAB6AM8LWk2yRjsuCqcvM76f+JPEbEqImYBM0iCRrXI0gZnAPcBRMRzQFuSRe9akkyfHU1ZliDRZG4Pr6CibZAOs4wkCRDVNv5co952iIglEdExIrpGRFeSuZmhETGlcHHNVpb/E38kuZgBSR1Jhp/eKmstG1eWNngbOBRAUi+SINHSniU7DvhqepXTgcCSiJhb6Uo1RNHhpmhCt4dXSsY2uBZoD/w+nbN/OyKGVqzSjSBjO1S9jO3wKHCEpOnAauDiiFhYuVpvXBnb4FvAryT9L8kQy2lV9uURSWNIhhQ7pnMvPwDaAETEL0nmYj4PzAQ+Ak6vTE03nJflMDOzOvmOazMzq5ODhJmZ1clBwszM6uQgYWZmdXKQMDOzOjlINAHFVpLMkP/odKXNF9NVR8/eyPW7Ir1REEkHp6uaTpO0k6SxRfLeJmnP9P1lG3DuzSVNlNRKUldJy9Nz17w2rSfvIEklLw8iabakl9P2fUzSDhtQRm4bXiBpi5xjD0naeiPW86W0zYquRJz1dyLpCUnblFpHa358CWwTIOkQYBnJQmB9Gpi3DfBP4DMRMUfSZkDXiJjRCFVF0i9JVrIcvQF5l0VE+wbm+QbQOiJ+pmQJ+vFZ20jSIOCiiDi6oXXNK2c2yWq2CyRdBbSPiPM3Rnml1Ku+ciX9ENgxIs4skifT70TSqUDniPjxxqmtNRfuSTQBGVaSrE8HkpsiF6ZlfVwTICTdLumXkv4i6Q1JR6f7W0m6VtLk9Ftnbc9D0iU535qvySnnOElfB74MfF/S79Jv9q/klHldzjfZ89L9T0vaLy1r8/Tb/+8kXSnpmznn/bGkQh+8JwN/qq8BJH1G0rNpb+pZST0LpBmY0/uYKqlDuv/inHb4YYb2fgbYPc07LP15X5H0k5x2uD3d93J6I1luG54P7AhMkDQhPTZbUkdJP5F0Tk6dh0v61gbW8zlyFpKT9EdJz6e9wLPSfev8TtJ9p0j6e7pvpKRWaRHjSNZispam0muV+5W8gK7UsSZ9hry3AfOAMSQfqpuk+28HHiH5MtCdZB2ZtsBZwHfTNJsBU4BuwJHAs6x9Lsa2OeUcV+B9bZ1JnqFxP8m3/ty8T5M+UwJYlvfzvpC+34Tk2QTb5f1cmwL/zsuzHJiWvm5O92+Zc97DgPvT94NIeh4Afwb+I33fniSwHkHyzAOldRgPHFKgfWcDHdP3vwB+QvJB/zawfVrWU8B/AfsCj+fk3bpAu9WWl7sN7ANMzNk/HdhlA+t5I3BWzrGa38fmwCs1bZ33O+mVtlObdPsW4Ks5x/+R/zvyq/pfWVaBtSYuIr4uqS/JB+RFwOGsXYn1vkhW3fyHpLeAPUg+dPaSdFyaZiuSIHIYMDrS52JEREN6N4eRPHjqkyx5I2K2pIVK1rz6NDA11l+2oiOQv2rqmxHRL2/fVsAdkrqTLP/QpsAp/wrckH5j/kMkQ3NHkLTF1DRNe5J2eKZA/gmSVgMvAd8lWeX36YiYD5CWewjJszR2lfRz4EHgsfraIVdETJX0KUk7kgSfRRHxdtr7aEg9P03ypeG7OfvPl/TF9P3Oaf789j6UJMhNVrK0zOZpOTXmkQTHqllexIpzkGgG0i7/8+nmuIj4fn6aiHgZeFnSXcAs1gaJ/EmnIPlGel5EPJp3niEF0meu5gbkvY2knjsAvylwfDlJz6eYK4EJEfHFdN7i6fwEEXGNpAdJ1tGZpGQSWcDVETEywzkGR84cglT4oVoRsUjS3sDngG+QDM99LUP5NcaSLJK5A8mDfGhoPYEPSXouVwAXpnMzhwEDIuIjSU9TuF0F3BERl9ZRdluS34m1IJ6TaAYiYnVE9Etf6wQISe3TD4Ea/UgmsmscL2kTSbuRPGpyBsmibP+jZNIbST0ktSP51vs1pVfeSNq2AdV8DPhvJUvF15V3Vc05Uw8AQ4D90zrl/9yLgFaSigWKrYB/pe9PK5RA0m4R8XJE/IRkeG2P9Jxfk9Q+TbOTpE8VOVeNvwED07mEViTj9ROVrPi6SUTcD3yP5NGW+ZaSzCUVcg/JApnHkQQMGlrPiFgOXECy+ui2JO2zKA0Qe5As4V4j93fyJHBcTdmStlV6hVQaFHcgGdKyFsQ9iSZABVaSjIhfZ80OXCJpJMm3vA9Z94NyBjCRZEjnvyNihaTbSOcE0v/884H/iohHJPUDpkhaSbKCZdbLVm8jWQ77JUmrgF+RjN/nGpUefyEiTo6Ilenk7eKIWF1HuY8BnwWeqOfcI0iGmy4kmRso5AJJg0lWZJ0OPBwRHytZwvq5tGOwDDiFdYdYCoqIuZIuBSaQ/A4eiog/pb2I0UoewARQ6Fv5KOBhSXMjYnBeua+mk+r/inRJ6Yh4rKH1TOs3hqQ3M4IkgL9E8vcwKa8utb8TSd8FHkvrvyrN/0+SYahJNcOJ1nL4EtgqJul2konbeu9lqJT0g+gF4PiI+EcdafYBLoyIr5S1crYOST8jGep8stJ1sfLycJNVhJIb7GYCT9YVICCZzCWZjG1VVxori1ccIFom9yTMzKxO7kmYmVmdHCTMzKxODhJmZlYnBwkzM6uTg4SZmdXp/wHEOp53KIgqpQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#First, we expand the list of lesions of patient_info, setting their oreder as a lesion index\n",
    "#We do it by hand to make sure that lists are expanded in order, so that lesion_id can be added\n",
    "data= []\n",
    "for row in patient_info.itertuples():\n",
    "    lst = row[1]\n",
    "    for i,col2 in enumerate(lst):\n",
    "        data.append([row[0], i, col2, row[-1]])\n",
    "        \n",
    "patient_info_expanded = pd.DataFrame(data=data, columns=['pid', 'lesion_id', 'class_target', 'massive_class_target'])\n",
    "\n",
    "#Keep from patient info only the patients in test set\n",
    "pids= [p[1] for p in test_results_list]\n",
    "patients_test= patient_info_expanded.loc[patient_info_expanded.pid.isin(pids),].set_index(['pid', 'lesion_id'])\n",
    "patients_test['score']= 0.\n",
    "patients_test.loc[patients_test.class_target == 10.,:]= 0.\n",
    "patients_test.loc[patients_test['massive_class_target'].isna(), 'massive_class_target']= -1\n",
    "patients_test['ClinSig']= patients_test['class_target'] >= GGG_THRESHOLD\n",
    "patients_test['ClinSig2']= np.nan\n",
    "patients_test['det_class']= 1\n",
    "columns_to_keep= ['ClinSig', 'ClinSig2', 'massive_class_target', 'class_target', 'score']\n",
    "patients_test= patients_test[columns_to_keep]\n",
    "\n",
    "#Match lesions\n",
    "patients_test= match_lesions(test_results_list, patients_test, drop_nans=True, class_col='det_class', \n",
    "                             cl=GGG_THRESHOLD + 1,  match_iou= 1e-5, spacing=SPACING,\n",
    "                             use_distance=USE_DISTANCE, distance_threshold=DISTANCE_THRESHOLD )\n",
    "patients_test['det_class']-= 1 #To make it correspond with GGG\n",
    "patients_test.loc[patients_test.det_class.isna(), 'det_class']=0\n",
    "patients_test['det_class']= patients_test['det_class'].astype(int)\n",
    "patients_test['det_class_thresholded']= patients_test['det_class']\n",
    "patients_test.loc[patients_test['score'] < CLASS_CONFIDENCE_THRESHOLD, 'det_class_thresholded']= 0\n",
    "\n",
    "#Prepare more dataframes in case that we want to individualize\n",
    "patients_test_prostatex= patients_test.loc[\n",
    "    patients_test.index.isin([p for p in patients_test.index.get_level_values(0) \n",
    "                              if p.lower().startswith('prostatex') and int(p[-4:]) < 204], level=0), :]\n",
    "\n",
    "#print results at different thresholds\n",
    "print('Lesion-level (prostatex-like). Sig if GGG >= %d'%GGG_THRESHOLD)\n",
    "for plot_subset, name in zip([patients_test_prostatex], ['ProstateX']):\n",
    "    print(' ', name)\n",
    "    print('    N=%3d    AUC  | Acc.   | Sens.  | Specif.'%(len(plot_subset)))\n",
    "    for t in threshold_list + [*get_optimal_thresholds(plot_subset.ClinSig, plot_subset.score)]:\n",
    "        print('    %s: %.3f |  %.3f |  %.3f |  %.3f'%('t=%.2f'%t, *composite_score(\n",
    "            plot_subset.ClinSig, plot_subset.score, t=t)))\n",
    "\n",
    "    #Plot AUC\n",
    "    plot_auc(plot_subset.ClinSig, plot_subset.score, new_fig=name=='All', legend=name + ' ', annotate=False)\n",
    "    \n",
    "_=plt.title('Lesion-level ROC Curve')\n",
    "#plt.gcf().savefig('lesion_auc.png', dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Patient-level results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The worst lesion result for a given patient is taken as the Grount Truth, while the highest scoring non-bening lesion is taken as the prediction score."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Patient-level. Sig if GGG >= 2\n",
      "  ProstateX\n",
      "    N= 45    AUC  | Acc.   | Sens.  | Specif.\n",
      "    t=0.40: 0.814 |  0.800 |  0.308 |  1.000\n",
      "    t=0.26: 0.814 |  0.822 |  0.692 |  0.875\n",
      "    t=0.09: 0.814 |  0.689 |  0.923 |  0.594\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEWCAYAAACT7WsrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de7xUdb3/8ddb1BDB684j3kDiKqBI3jhHBVIL07BOmqKe1DziL1Mzj5hmF46mGSqZJ1I4HrEsUdM0Im+lgKVSIHgXDIWURBHEK6AIn98fa+3NsJnZs2D2npk9+/18POaxZ12+3/WZtfeez1rf71rfpYjAzMwsn80qHYCZmVUvJwkzMyvIScLMzApykjAzs4KcJMzMrCAnCTMzK8hJwgCQ9B1JN5ZxeyGpewtvY5qk/2zJbZjVOieJVkzSQkkrJb0v6Q1JEyV1zFBuiKRFufMi4oqIKPkLVVLXNAFsXmpdLUnSaEmr0333tqTHJA1qtM52kq6X9LqkFZKekXRanrpOlDQrrWuxpPskHdzEtg+QdG+63bck/S1fvWbVwEmi9ftCRHQEBgL7A9+tcDytye3pvqsDpgK/qV8gaUvgT0AXYBCwLTAKuFLS+TnrnQ9cC1wB/AuwB/Bz4Jh8G0wT0cPAdKA7sCPwdeDITfkAktptSjmzzCLCr1b6AhYCh+dMXwVMSd+fBrwAvAe8DJyZzt8aWAmsBd5PX7sAo4Ff5dR1EPAY8DbwFDAkZ9k04DLg0bT+B4G6dNkrQOTUPahA7AF0T99/Arg6LfsGcAOwVbrsBeDonHKbA0uBgRnj/M8C22/8efdKY/pkOn06sATYulG549PPtQ1J4ngfOG4jfmd/AcY1sfxU4C9N7KubgeuBe4EPSA4KXgfa5az/JeDp9P1mwEXAS8Ay4A5gh0r/7frVel4+k6gRknYHPg/MSWctAY4m+TI7DfiJpIER8QHJUetrEdExfb3WqK5dgT8APwR2AC4A7pL0yZzVTkzr3QnYMl0H4ND053Zp3Y9nCP/HQE9gAMnR9a7A99Nlk4AROet+DlgaEbMzxllUetbwVZIv0eXp7COA+9L9lesuoD3J2cWg9P3dGbfTIS1z58bEl8eJwOVAJ5Lk+gHwmUbLb03fnwt8ERhMcjCwHBhX4vatDXGSaP3ukfQ2yRHqdJJmDyLiDxHxUiSmkxztH5KxzpOBeyPi3ohYGxF/BGaRJKF6EyPixYhYSXJ0OmBTgpck4AzgWxHxVkS8l36GE9JVbgWGp1+wsP4XYJY4m/KVdN+tTGM4NiI+TpfVAYsbF0iXL02X70iSsD5uvF4B25P8z21Q70b6XUQ8mn7mVeQkUkmdSD7/pHTdM4FLImJRRHxIcgZ1bLX3GVn1cJJo/b4YEdtFRJeIOCv90kbSkZJmpB2jb5N8cdRlrLMLcFzasfp2Wv5goHPOOq/nvF8BFOwwl/Rc2qn7vqTGieqTQAfgiZxt3Z/OJyLmkzQ5fSFNFMNZlySyxNmUOyJiO5K+hGeBT+csW5qvnvTLtS5dvgyo24gv3OUkzXxZ4yvk1UbTtwL/LukTwL8DsyPiH+myLsDdOfvnBWANyWc2K8pHEzUo/bK4i6QJ5XcRsVrSPYDSVYoN/fsqcEtEnLEJm9+g7ojo28T6S0mO5PtGxD8LrFN/pLwZ8HyaOEqNMze+pZLOBGZKujUiFpN0Wl8haetGTU5fBj4EZqTTq0iac4o2IUXECkmPp3VMLbDaByRJEwBJO+erqlG9z0v6B0kzYu6ZFiT76GsR8Wix+Mzy8ZlEbdqSpDP4TeBjSUcCn81Z/gawo6RtC5T/FcmR++cktZPUPr1sdrcM236T5Gi5W5ZAI2It8L8kfSY7QdInIulzOavdlsb/ddb/AiwlzsZxzAUeAC5MZ90CLAJ+k17Wu0Ua03XA6Ih4JyLeIek7GSfpi5I6pOsdKWlMgU1dCJwqaZSkHdPPu4+k29LlTwF9JQ2Q1J6keSiLW0n6Hw4l5yotkosALpfUJd3WJyXlvfLKLB8niRqUtuufS9JXsJzk6HJyzvK5JEfnL6fNELs0Kv8qySWc3yH50n+V5PLPon8vEbGCpFP10bTugzKE/G1gPjBD0rskR/G9cupcDDwO/Ctwe3PEWcBVwEhJO6Xt94endf4VeBcYS9K+f1VODGOB80muMqqP4WzgnnwbiIjHSDqZP0Oy/98CJpBcrUREvAhcmu6Dv5P0NWUxCRgCPBwRS3Pm/5Tkd/+gpPdIzoAOzFinGYrwQ4fMzCw/n0mYmVlBThJmZlaQk4SZmRXkJGFmZgWVep+Ee73NzDaSpAciYlil48jCZxJmZuXX5OgHkm6StETSswWWS9J1kuZLelrSwJYJ00nCzKwa3Qw0daZxJNAjfY0kGRm4RThJmJlVmYh4BHiriVWOAX6ZDuA5A9hOUqljguXlJGFm1vrsyvoDPS5K5zU7D/BnZm3PlV1g1dsw4jboPADG9l63bOApMPw6GH8oLH4qmddxZ7hgHkz9EUy/ct26I6clPycMWTdv8EUw9GKYNRH2K/hU2jpJs3KmJ0TEhI34BMozr0UuJCp1WA5f3WRmrc/obWH0OxXbvKQnImK/Iut0JXnSZL88y8YD0yJiUjo9j+SpjKU+q2QDbm4ys7bn6GsrHUGpJgNfTa9yOgh4pyUSBLi5yczaosLNQFVBUv2ovnWSFgE/ALYAiIgbSEYN/jzJ6MkrSB4l3DKxuLnJzNqcVtDcVC3c3GRmZgU5SZiZWUFOEmbW9vRsFcMmVQX3SZiZlZn7JMzMqtmtx1c6glbDScLM2p4X7690BK2Gk4SZmRXkJGFmZgW549rMrMzccW1mVs1mTax0BK2Gk4SZtT1Tzqt0BK2GB/iz5lU/Tn/nfeDMR2DyuTD7F+uWnz8XFj8Jk05YN+/oa5MB10Zvu25ez2Fw4u3JpYq5V6KMfic5Csz9J2/JZwJc3Qvefz2Z589UO5+p/XZYNu6TsOZ1da/kH9zMCmpNfRJOEmZmZdaakoT7JKx5Tf1RpSMws2bkJGHNK7cN2cxaPScJMzMryEnCzMwKcpKw5jVyWqUjMLNm5CRhZmYFOUlY88q90cnMWj0nCTMzK8hJwszMCnKSsOY1+KJKR2BmzcjDcpiZlZmH5bC26+pelY7AzHJI2l/STyXNlrRY0suSJks6U1KnYuU9VLg1r/rhms2s4iRNAZYBvwOuAZYA7YGewFDgD5LGRMSUQnU4SZiZ1a7TI+KNRvNWAX9LXz+WtFNTFbi5yZpX530qHYGZpXIThKTdJA1N339C0tbpOkuaqsNJwprXmY9UOgIza0TS14DJwI3prC4kTVBFOUlY85p8bqUjMLMNnQscBLwLEBEvAk02M9VzkrDmlftMYTOrFqsi4qP6CUntAGUp6CRhZlb7HpV0IdA+7Ze4HSh4RVMuJwkzs9p3IfAeMBf4JvAQ8J0sBZ0krHmdP7fSEZjVBEnDJM2TNF/SBuPdSNpD0lRJcyQ9LenzTVT39Yi4PiK+FBFfjIjrgW9kicNJwprX4icrHYFZq5f2GYwDjgT2AkZI2qvRat8F7oiIfYETgJ83UeXX8sw7PUssvpnOmtekE2D0O5WOwqy1OwCYHxEvA0i6DTgGeD5nnQC2Sd9vC7zWuBJJx5MkkD0l/TZnUSfg7SyBOEmYmZVfnaRZOdMTImJCzvSuwKs504uAAxvVMRp4UNI5wNbA4Xm28zeSYTl2IzkzqfceMCdLoE4SZmblt7TIKLD5Lk9tPOr2CODmiLhG0iDgFkn9ImJtQ4GIBcAC4E+bGqj7JKx5HX1tpSMwqwWLgN1zpndjw+ak04E7ACLicZKB++ryVZaOBDtD0juSVkn6UNK7WQJxkrDmtd9plY7ArBbMBHpI2lPSliT9CpMbrfMKcBiApD4kSeLNAvX9HDgFeJmkP+JsINMRnZOENa/R21Y6ArNWLyI+JvkifwB4geQqpuckXSppeLrafwFnSHoKmAScGoWfIrdZRMwDNo+I1RHxv+Tvw9iA+yTMzKpQRNwL3Nto3vdz3j8P/FvG6j5Iz0ieknQFsBjomKWgzyTMzGrfqSTf92cDa4AewLFZCvpMoq27sgusejt5DsSZjySjuOYO0nf+3OQGuUknrJt39LVJ30Nu01LPYXDi7dAl64GNmZVL/f0WJA8c+h6ApMaX1Oalwk1Y2bZdSmGrAlf3ggvmVToKszZF0hNFLoFtru1sBnyZ5L6LByLiBUnDSMZt2j4i+hetw0nCzKy8ypgkbgK6kVwttT/wd2AIcHFE3JmlDjc3tXVTfwRDL650FGbWMg4E9o6INZK2ApYC3SNicdYK3HHd1k2/stIRmFnL+TAi1gBExEpg3sYkCPCZhJlZLestaXb6XkCvdFpARMTAYhU4SZiZ1a6iHdPFOEm0dSOnVToCM2shEfFSqXW4T8LMzApykmjrJgypdARmVsWcJMzM2gBJW0rqvrHlnCTMzGqcpKOAZ4A/ptMDJN2dpayTRFs3+KJKR2BmLe9Skhvr3gaIiCeBTGcVThJtne+2NmsLVkfE243mZRpWyUmirbu6V6UjMLOW94KkrwCbpU+7uxaYkaWgk0Rb9/7rlY7AzFre2cCngbXAb0mGDD8vS8Hav5lu4ufhtHvh1uPhxfvXzR/9DsyaCFNy9tOI26DzABjbe928gafA8Otg/KGw+KlkXsedk+G1p/5o/bGPRk5LfuZeVjr4oqRJ5+pe676Qm+vZDc3xmdpv18TOM7Ma0S0ivg18e2ML1v5Q4aO3Tb48zcyqRLmGCs/Z3p+BHYA7gNsjYm7Wsm5uMjOrcRFxCPBZ4D3gF5LmSMp0aaOThJlZGxAR/4yIsSTPu34GuCxLudpvbjIzqzIVaG7qARwPHAu8D9wO3Jnl2RK1fyYxa2KlIzAzq7RbgZXA8Ig4OCL+J+vDh2r/6qYp5yVXB5mZtVERsf+mlq39JGFm1kZJmhQRIyTNYf3uAT+ZzszMGJX+PHZTK6j9PokRt1U6AjOzioiIRenb0yPipdwXcHqWOmo/SXQeUOkIzMwqbVieeUdlKVj7SSJ3OAozszZE0plpf0QvSbNzXn8HXshSh/skzMxq1x3AQ8CPgNw7rN+LiCVZKnCSMDOrXR9FxHxJG/Q/SNomIt4tVkHtJ4mBp1Q6AjOzSrkTOBJ4juQSWOUsC2CPYhV4WA4zszIr97Acpaj9juvxh1Y6AjOzipJ0kKQO6fsRksZI2i1L2eprbvrNaXDcxOZ7yI8fqmNmrZCkYcBPgXbAjRFxZZ51vgKMJmnVeSoiTixQ3QRgH0l7A98BbgZ+DQwuGkfVNTf5IUFmVuOKNTdJage8CBwBLAJmAiMi4vmcdXqQXL30mYhYLmmnQlcsSZodEQMlfQ9YHBE31s8rFmvtNzeZmbU+BwDzI+LliPgIuA04ptE6ZwDjImI5QJFLWj+QNAr4D+APkjYDtsgSiJOEmVn51UmalfMa2Wj5rsCrOdOL0nm5egI9JT0qaUbaPFXI8SRXNp2ZDhG+GzA2S6DV1ydxfuZHr5qZtVZLi1zdpDzzGjfvbw70AIaQfOn/WVK/iHh7g4IRr0m6CdgvTSazIiLTw3aq70xi8ZOVjsDMrNIWAbvnTO8GvJZnnd9FxOqIWADMI0kaG5D0ZWA2SXPTV4FZkr6UJZDqSxKTTqh0BGZmlTYT6CFpT0lbAicAkxutcw8wFEBSHUnz08sF6vs+sH9EnJReAXUgyVVRRVVfkjAza+Mi4mPgbOABkoH47oiI5yRdKml4utoDwDJJzwNTgVERsaxAlZtFxBs502+S8fvfl8CamZVZue+4ljQW6E3yrGtIzkzmRsQFxcpWX8f10ddWOgIzs1rzX8BxwMEkneK/IBnXqajqO5MwM6tx5TyTkHQ00B14JiIe2tjy1dcnMXrbSkdgZlYTJP0PyXMkdgXGSPrOxtZRfc1NZmbWXIYCAyLiY0lbA9OBKzamguo7kzAzs+byUXqlFBHxAZvwnV99ZxI9m7qz3MzMNkJvSbPT9yJ91nX6PrIM8OeOazOzMitXx7WkTzW1PCJeKlZH9Z1J3Ho8nHh7paMwM2v1siSBYqqvT+LF+ysdgZmZpaovSZiZWdVwkjAzawMkbSmp+8aWq74k4XGbzMyalaSjgGeAP6bTAyTdnaVs9SWJWZmeg2FmZtldSjI8+NsAEfEkyVAdRVVfkphyXqUjMDOrNavzPLEu0y0M1XcJrJmZNbcXJH0F2EzSnsA3gRlZClbfmYSZmTW3s4FPA2uBu4EPgUzNNtV3x/W8+6DXkc1erZlZtSj3Q4dKUX3NTZ0HVDoCM7OaIumP5Dmoj4jPFitbfUlibG9fBmtm1ry+m/O+PfBlkianoqovSZiZWbOKiL82mjVd0vQsZZ0kzMxqnKRtciY3I+nE7pylbPUliYGnVDoCM7Na8xxJn4SAj4EFwBlZClbf1U1mZjWunFc3SdoMOCAiMt0X0Vj13Scx/tBKR2BmVjMiYi1w7aaWr74ksfipSkdgZlZr/ijpmE0pWH19EmZm1tzOBraV9CGwknXPuN6hWMHqSxIdd650BGZmNUHSHhHxClC3yXW449rMrLzK1XEtaXZEDCyljurrk5j6o0pHYGZWK1RqBdXX3DT9Shh6caWjMDOrBbtKuq7Qwog4t1gF1ZckzMysuawEniilAicJM7PatSwiflFKBdXXJzFyWqUjMDOrFR+VWkH1JQkzM0PSMEnzJM2XdFET6x0rKSTlu1rqhCLbkKTdmlqn+pLEhCGVjsDMrKIktQPGAUcCewEjJO2VZ71OwLlA46HA610l6S5JX5XUV9JOkvaQ9BlJlwGPAn2aisV9EmZm1ecAYH5EvAwg6TbgGOD5RutdBowBLshXSUQclyaXk4CvkQwPvgJ4AbgXuDwiVjUViJOEmVn12RV4NWd6EXBg7gqS9gV2j4gpkvImCYCIeB64ZFMDqb4kMbhg05uZWa2okzQrZ3pCREzImc53E1zDCBfp8N8/AU5tmfDWqb4k4RvpzKz2LS0yLMciYPec6d2A13KmOwH9gGmSAHYGJksaHhG5yadk1ddxfXWvSkdgZlZpM4EekvaUtCXJVUqT6xdGxDsRURcRXSOiKzADaPYEAeU6k7iyC6x6G0bcBp0HwNje65YNPAWGX5c8bGjxU9B+u7KEZGZWrSLiY0lnAw8A7YCbIuI5SZcCsyJictM1rE/SXcBNwH3pQ4iyly3LKLCjt4XR75SyHTOzmlHOx5em2zscOA04CPgNcHNEzM1StjzNTUdv8pPzzMysRBHxp4g4CRgILCR5Ut1jkk6TtEVTZf08CTOzMiv3mUS6zR2Bk4H/IOkE/zVwMNA/IoYUKleeM4nR25ZlM2ZmtiFJvwX+DHQAvhARwyPi9og4B+jYVNnquwTWzMya240RcW/uDEmfiIgPi53RVN8lsGZm1tx+mGfe41kKludMouewsmzGzMzWkbQzyRAfW6XDeNTfyb0NSdNT8TrccW1mVl7l6riWdArJ0B37Abk32r1Hchnsb4vWUZYkcevxcOLtpWzHzKxmVOA+iS9HxF2bUrY8zU0v3l+WzZiZ2TqSTo6IXwFdJZ3feHlEjC1Wh69uMjOrXVunP5u8zLUpThJmZjUqIsanb38eEW9uSh3uuDYzK7MK9En8HVgA3A78NiKWZy1bnvskZk0sy2bMzGxDEdED+C7QF3hC0hRJJ2cpW54kMeW8smzGzMzyi4i/RcT5JM/Pfgv4RZZypfVJXNkFtu8CZz4Ck8+F2TnbPH8uLH4SJp3gZ0SYmVWQpG2AL5E8vOhTwN0kyaJ42ZL6JEZvG35OhJnZxqlAn8QC4B7gjojINBxHvdLOJDruXFJxMzMri26xiWcEvrrJzKzMyjgsx7URcZ6k35Pn+zoihhero7Qziak/gqEXl1SFmZm1mFvSn1dvagWlXd00/cqSipuZWcuJiCfStwMiYnruCxiQpQ4/T8LMrPadkmfeqVkKelgOM7MaJWkEcCKwp6TJOYs6Acuy1FFakhg5raTiZmbWoh4DFgN1wDU5898Dns5Sgc8kzMxqVET8A/gHMGhT6yitT2LCkJKKm5lZy5H0l/Tne5LezXm9J+ndLHX4TMLMrEZFxMHpz06bWoevbjIzq3GSPiXpE+n7IZLOlZRpUL3SksTgi0oqbmZmZXEXsEZSd+D/gD2BW7MULC1J+G5rM7PWYG1EfEwyEuy1EfEtoHOWgqUliat7lVTczMzKYnV6z8QpwJR03hZZCpaWJN5/vaTiZmZWFqeRXAZ7eUQskLQn8KssBX11k5lZjYuI54Fzc6YXAJkG3ystSXTep6TiZmbW8iT9GzAa6ELyvS8gIqJb0bJ+noSZWXlV4Ml0c4FvAU8Aa+rnR0TR8ZtK65OYfG7xdczMbKNJGiZpnqT5kja430DS+ZKel/S0pIckdWmiunci4r6IWBIRy+pfWeIoLUnM/kVJxc3MbEOS2gHjgCOBvYARkvZqtNocYL+I2Bu4ExjTRJVTJV0laZCkgfWvLLG449rMrPocAMyPiJcBJN0GHAM8X79CREzNWX8GcHIT9R2Y/sxt4grgM8UCcZIwMyu/OkmzcqYnRMSEnOldgVdzphex7os+n9OB+wotjIihmxQlpSaJ8+eWVNzMrI1aWqTjWnnm5b1QSNLJJGcIgwtWJv0LcAWwS0QcmTZdDYqI/ysWaGl9EoufLKm4mZnltQjYPWd6N+C1xitJOhy4BBgeER82Ud/NwAPALun0i8B5WQIpLUlMOqGk4mZmltdMoIekPSVtCZwA5D5+FEn7AuNJEsSSIvXVRcQdwFqAdBynNU0XSXiocDOzKpN+iZ9NcvT/AnBHRDwn6VJJw9PVrgI6Ar+R9GSjZ1g39oGkHUmbrCQdBLyTJRZ3XJuZVaGIuBe4t9G87+e8P3wjqjuf5EzkU5IeBT4JHJulYGlJ4uhrSypuZmYtLyJmSxoM9CLpFJ8XEauzlC2tuWm/00oqbmZmLUfS/pJ2hoYmrE8DlwPXSNohSx2lJYnR25ZU3MzMWtR44CMASYeSjPz6S5L+iAlNlGvgPgkzs9rVLiLeSt8fT3LT3l3AXZIy3cPgq5vMzGpXO0n1JwOHAQ/nLMt0klDamUTPYSUVNzOzFjUJmC5pKbAS+DOApO5kvATWz5MwMyuzcj5PIr0nojPwYER8kM7rCXSMiNnFypfW3HTr8SUVNzOzlhURMyLi7voEkc57MUuCgFKTxIv3l1TczMyqmzuuzcysICcJMzMryB3XZmZlVs6O61KVdiYxa2IzhWFmZtWotCQxJdMzK8zMrJXysBxWs1avXs2iRYtYtWpVpUOxNqp9+/bstttubLHFFpUOZZM5SVjNWrRoEZ06daJr165I+R4ZbNZyIoJly5axaNEi9txzz0qHs8lKa24acVszhWHW/FatWsWOO+7oBGEVIYkdd9yx1Z/JlpYkOg9opjDMWoYThFVSLfz9lZYkxvZupjDMzKwa+WY6sxbSrl07BgwYQL9+/TjuuONYsWJFyXXefPPNvPbaa8223vDhw7nlllsaps844wyuuuqqouWuvfZafvnLXzZMf/zxx9TV1XHxxRevt17Xrl1ZunRpw/S0adM4+uijG6bvu+8+9ttvP/r06UPv3r254IILim67mCeeeIL+/fvTvXt3zj33XPLdC/bOO+/whS98gX322Ye+ffsyceK6y/mHDRvGdtttt16cuc455xw6duzYMP2zn/1svfI1JyI2/fWDbcKsWj3//PMV3f7WW2/d8P7EE0+Ma665Zr3la9eujTVr1mxUnYMHD46ZM2c223oLFiyIrl27xvLly+PRRx+N/v37x0cffdRkmdWrV0f//v1j9erVDfP+8Ic/xL/+679Gt27dYu3atQ3zu3TpEm+++WbD9NSpU+Ooo46KiIhnnnkmunXrFi+88EJDvePGjSsaczH7779/PPbYY7F27doYNmxY3HvvvRusc/nll8eFF14YERFLliyJ7bffPj788MOIiPjTn/4UkydPbogz18yZM+Pkk09e73f7wQcfxIABAwrGk+/vEJgVpXz3lvFV2pnEwFOaJ1OZ1bhDDjmE+fPns3DhQvr06cNZZ53FwIEDefXVV5k0aRL9+/enX79+fPvb3wZgzZo1nHrqqfTr14/+/fvzk5/8hDvvvJNZs2Zx0kknMWDAAFauXMmll17K/vvvT79+/Rg5ciQRkXe9J554gsGDB/PpT3+az33ucyxevBhIjvRHjhzJhRdeyFlnncXPfvazopdrPvzwwwwcOJDNN193ceSkSZP45je/yR577MGMGTMy7ZMxY8ZwySWX0Lt30my9+eabc9ZZZ23K7m2wePFi3n33XQYNGoQkvvrVr3LPPfdssJ4k3nvvPSKC999/nx122KHh8xx22GF06tRpgzJr1qxh1KhRjBkzZr35HTp0oGvXrvztb38rKfaqVWKWMataGxzBPXxFxA+2Wff65+zklTvv4SuSda/quW7eDYck8353zvrrvvNak9uvP9pcvXp1DB8+PH7+85/HggULQlI8/vjjERHxz3/+M3bfffdYsmRJrF69OoYOHRp33313zJo1Kw4//PCGupYvXx4RG54hLFu2rOH9ySefHJMnT95gvY8++igGDRoUS5YsiYiI2267LU477bSGch999FHsvvvuceKJJ2bar9///vfjuuuua5hesWJFdO7cOT744IMYP358nHPOOQ3LmjqT2HfffePJJ58sur2HH3449tlnnw1egwYN2mDdmTNnxmGHHdYw/cgjj+Q9I3j33XdjyJAhsfPOO8fWW28dU6ZMWW95bpz1rr322hg7dmxErH+WGBHxwx/+MK6++uq88bf2M4nS7pMYfyic+UjzZCuzljb04uTV2Og8D+i6YN6G84Zfl7wyWrlyJQMGJFcAHnLIIZx++um89tprdOnShYMOOgiAmTNnMmTIED75yU8CcNJJJ/HII4/wve99j5dffplzzjmHo446is9+9rN5tzF16lTGjBnDihUreOutt+jbty9f+MIX1ltn3rx5PPvssxxxxBFAckTcuXPnhuVPP/00EcHcuXNZu3Ytm23WdAPD4sWL6dOnT8P0lClTGDp0KB06dODLX7nmYLUAAA80SURBVP4yl112GT/5yU9o165d3qt7NvaKn6FDh/Lkk5kex5y3/yHf9h544AEGDBjAww8/zEsvvcQRRxzBIYccwjbbbJO33tdee43f/OY3TJs2Le/ynXbaiblz52aKsbUpLUksfqqZwjCrPVtttVXeL7ett9664X2+LzWA7bffnqeeeooHHniAcePGcccdd3DTTTett86qVas466yzmDVrFrvvvjujR4/Oe01+RNC3b18ef/zxDZatXbuWs846i1tuuYUbbriB66+/nm984xtFP1fudiZNmsSjjz5K165dAVi2bBlTp07l8MMPZ8cdd2T58uXU1dUB8NZbbzW879u3L0888QT77LNPk9ubOnUq3/rWtzaY36FDBx577LH15u22224sWrSoYXrRokXssssuG5SdOHEiF110EZLo3r07e+65J3PnzuWAAw7IG8OcOXOYP38+3bt3B2DFihV0796d+fPnA8nvYquttmryc7RWvrrJrIIOPPBApk+fztKlS1mzZg2TJk1i8ODBLF26lLVr1zYcmc+enTxErFOnTrz33nsADV/UdXV1vP/++9x5550N9eau16tXL958882GJLF69Wqee+45AMaPH0+PHj0YMmQIY8eOZcyYMbz55ptNxtynT5+GL8d3332Xv/zlL7zyyissXLiQhQsXMm7cOCZNmgTAkCFDGq6eWrNmDb/61a8YOnQoAKNGjeKKK67gxRdfBJKENXbs2A22V38m0fjVOEEAdO7cmU6dOjFjxgwigl/+8pccc8wxG6y3xx578NBDDwHwxhtvMG/ePLp161bwMx911FG8/vrrDZ+xQ4cODfsA4MUXX6Rfv35N7rdWq6T2qqt65m2DM6sG1XR1U70FCxZE375915v361//Ovr16xd9+/aNUaNGRUTEk08+Gfvuu29D+3v9FTp33nln9OzZM/bZZ59YsWJFXHLJJfGpT30qDjvssDj11FPjBz/4Qd715syZE4ccckjsvffesddee8WECRPijTfeiC5dusTixYsbYrnmmmvi1FNPbfJzLVy4MA45JOmnmThxYhx//PHrLV+2bFnU1dXFqlWr4u23344RI0bE3nvvHf37949Ro0atd0XX73//+xg4cGD07t07+vTpExdccEHGvVvYzJkzo2/fvtGtW7f4xje+0XC11fXXXx/XX399RCR9QUcccUTDfr/lllsayh988MFRV1cX7du3j1133TXuv//+DbbR+He77777rtf3kqu190n4eRJWs1544YX12s6t+XzpS19izJgx9OjRo9KhVNycOXMYO3bseveb5Mr3d9h2nicx9UfNFIaZtSZXXnllw2W0bd3SpUu57LLLKh1GiyktSUy/spnCMLPWpFevXhx66KGVDqMqHHHEEQ2d9rXIHddW00psTjUrSS38/TlJWM1q3749y5Ytq4l/VGt9IpLnSbRv377SoZSktPskRk5rliDMWkL9NfPFLuk0ayn1T6ZrzfxkOqtZW2yxRat+IphZNSituWnCkOaJwszM1iNpmKR5kuZLuijP8k9Iuj1d/ldJXVsiDvdJmJlVGUntgHHAkcBewAhJezVa7XRgeUR0B34C/LglYnGSMDOrPgcA8yPi5Yj4CLgNaDy+yDHAL9L3dwKHqQWel1pSn4T++93n+G+17qd8N486YGnRtWqb90HC+yHh/dD0PuglaVbO9ISImJAzvSvwas70IuDARnU0rBMRH0t6B9ixiW1uklI7rle1llvLW5KkWW19P3gfJLwfEt4PJe+DfGcEja/lzrJOydzcZGZWfRYBu+dM7wY0fmh5wzqSNge2Bd5q7kCcJMzMqs9MoIekPSVtCZwATG60zmSg/hnSxwIPRwvcOVpqc9OE4qu0Cd4P3gf1vB8S3g8l7IO0j+Fs4AGgHXBTRDwn6VKSYcYnA/8H3CJpPskZxAnNEXRjpQ4VbmZmNczNTWZmVpCThJmZFZQpSVTL7eGVlGEfnC/peUlPS3pIUpdKxNnSiu2HnPWOlRSSavIyyCz7QdJX0r+J5yTdWu4YW1qG/4k9JE2VNCf9v/h8JeJsSZJukrRE0rMFlkvSdek+elrSwHLHWLJizzcl6TR5CegGbAk8BezVaJ2zgBvS9ycAt1f6uazN+cq4D4YCHdL3X6+1fZB1P6TrdQIeAWYA+1U67gr9PfQA5gDbp9M7VTruCuyDCcDX0/d7AQsrHXcL7IdDgYHAswWWfx64j+SehoOAv1Y65o19ZTmTqJrbwyuo6D6IiKkRsSKdnEFyXXOtyfK3AHAZMAao1bvxs+yHM4BxEbEcICKWlDnGlpZlHwSwTfp+Wza8zr/Vi4hHaPrehGOAX0ZiBrCdpM7lia55ZEkS+W4P37XQOhHxMVB/e3ityLIPcp1OcvRQa4ruB0n7ArtHxJRyBlZmWf4eegI9JT0qaYakYWWLrjyy7IPRwMmSFgH3AueUJ7SqsrHfHVUny30SVXN7eAVl/nySTgb2Awa3aESV0eR+kLQZyWiUp5YroArJ8vewOUmT0xCSs8o/S+oXEW+3cGzlkmUfjABujohrJA0iuaa/X0Ssbfnwqkar/27MciZRNbeHV1CWfYCkw4FLgOER8WGZYiunYvuhE9APmCZpIUkb7OQa7LzO+j/xu4hYHRELgHkkSaNWZNkHpwN3AETE40B7kkHv2pJM3x3VLEuSqJrbwyuo6D5Im1nGkySIWmt/rtfkfoiIdyKiLiK6RkRXkr6Z4RExK391rVaW/4l7SC5mQFIdSfPTy2WNsmVl2QevAIcBSOpDkiTa2rNkJwNfTa9yOgh4JyIWVzqojVG0uSmq6PbwSsm4D64COgK/SfvsX4mI4RULugVk3A81L+N+eAD4rKTngTXAqIhYVrmom1fGffBfwP9K+hZJE8upNXbwiKRJJE2KdWnfyw+ALQAi4gaSvpjPA/OBFcBplYl003lYDjMzK8h3XJuZWUFOEmZmVpCThJmZFeQkYWZmBTlJmJlZQU4SVaDYSJIZyh+djrT5VDrq6JnNHN+l6Y2CSDokHdX0SUm7SrqzSNkbJe2Vvv/OJmx7K0nTJbWT1FXSynTb9a8tmyg7RFLJw4NIWijpmXT/Pihp502oI3cfniepQ86yeyVt14xxPp3us6IjEWf9nUj6k6TtS43RWh9fAlsFJB0KvE8yEFi/jSy7BfAP4ICIWCTpE0DXiJjXAqEi6QaSkSwnbkLZ9yOi40aW+QaweUT8VMkQ9FOy7iNJQ4ALIuLojY21UT0LSUazXSrpCqBjRJzbHPWVEldT9Ur6b2CXiDijSJlMvxNJpwC7RcTlzROttRY+k6gCGUaSbEonkpsil6V1fVifICTdLOkGSX+W9KKko9P57SRdJWlmetTZcOYh6cKco+Yrc+o5VtJ/Al8Bvi/p1+mR/bM5dV6dcyR7Tjp/mqT90rq2So/+fy3pMknfzNnu5ZLyffGeBPyuqR0g6QBJj6VnU49J6pVnncE5Zx9zJHVK54/K2Q//nWF/PwJ0T8uOSD/vs5J+nLMfbk7nPZPeSJa7D88FdgGmSpqaLlsoqU7SjyWdlRPzaEn/tYlxPk7OQHKS7pH0RHoWODKdt97vJJ13sqS/pfPGS2qXVjGZZCwma2sqPVa5X8kL6EqBMekzlL0RWAJMIvlS3SydfzNwP8nBQA+ScWTaAyOB76brfAKYBewJHAk8xrrnYuyQU8+xed43xEzyDI27SI76c8tOI32mBPB+o887O32/GcmzCXZs9Lm2BF5vVGYl8GT6GpfO3yZnu4cDd6Xvh5CceQD8Hvi39H1HksT6WZJnHiiNYQpwaJ79uxCoS9//DPgxyRf9K8An07oeBr4IfBr4Y07Z7fLst4b6cqeBfYHpOfOfB/bYxDivBUbmLKv/fWwFPFu/rxv9Tvqk+2mLdPrnwFdzlv+98e/Ir9p/ZRkF1qpcRPynpP4kX5AXAEewbiTWOyIZdfPvkl4GepN86ewt6dh0nW1JksjhwMRIn4sRERtzdnM4yYOnPs5SNiIWSlqmZMyrfwHmxIbDVtQBjUdNfSkiBjSaty3wC0k9SIZ/2CLPJh8FxqZHzL+NpGnusyT7Yk66TkeS/fBInvJTJa0Bnga+SzLK77SIeBMgrfdQkmdpdJP0P8AfgAeb2g+5ImKOpJ0k7UKSfJZHxCvp2cfGxPkvJAcN382Zf66kL6Xvd0/LN97fh5EkuZlKhpbZKq2n3hKS5Fgzw4tYcU4SrUB6yv9EOjk5Ir7feJ2IeAZ4RtItwALWJYnGnU5BckR6TkQ80Gg7w/KsnznMTSh7I0mcOwM35Vm+kuTMp5jLgKkR8aW032Ja4xUi4kpJfyAZR2eGkk5kAT+KiPEZtjE0cvoQpPwP1YqI5ZL2AT4HfIOkee5rGeqvdyfJIJk7kzzIh42NE/iA5MzlUuD8tG/mcGBQRKyQNI38+1XALyLi4gJ1tyf5nVgb4j6JViAi1kTEgPS1XoKQ1DH9Eqg3gKQju95xkjaT9CmSR03OIxmU7etKOr2R1FPS1iRHvV9TeuWNpB02IswHgf+nZKj4QmVX128zdTcwDNg/janx514OtJNULFFsC/wzfX9qvhUkfSoinomIH5M0r/VOt/k1SR3TdXaVtFORbdX7KzA47UtoR9JeP13JiK+bRcRdwPdIHm3Z2HskfUn53EYyQOaxJAmDjY0zIlYC55GMProDyf5ZniaI3iRDuNfL/Z08BBxbX7ekHZReIZUmxZ1JmrSsDfGZRBVQnpEkI+L/shYHLpQ0nuQo7wPW/6KcB0wnadL5fxGxStKNpH0C6T//m8AXI+J+SQOAWZI+IhnBMutlqzeSDIf9tKTVwP+StN/nmpAunx0RJ0XER2nn7dsRsaZAvQ8CBwN/amLbY0iam84n6RvI5zxJQ0lGZH0euC8iPlQyhPXj6YnB+8DJrN/EkldELJZ0MTCV5Hdwb0T8Lj2LmKjkAUwA+Y7KJwD3SVocEUMb1ftc2qn+z0iHlI6IBzc2zjS+SSRnM2NIEvjTJH8PMxrF0vA7kfRd4ME0/tVp+X+QNEPNqG9OtLbDl8DWMEk3k3TcNnkvQ6WkX0SzgeMi4u8F1tkXOD8i/qOswdl6JP2UpKnzoUrHYuXl5iarCCU32M0HHiqUICDpzCXpjG1XaB0ri2edINomn0mYmVlBPpMwM7OCnCTMzKwgJwkzMyvIScLMzApykjAzs4L+P78AzmNN7Fv5AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Prepare two more dataframes in case that we want to individualize\n",
    "df= patient_df.reset_index()\n",
    "df_prostatex= df.loc[df.pid.isin([p for p in df.pid if p.lower().startswith('prostatex') and int(p[-4:]) < 204]), :]\n",
    "\n",
    "#Process dataframes to contain all, prostatex and ivo patients & plot results\n",
    "print('Patient-level. Sig if GGG >= %d'%GGG_THRESHOLD)\n",
    "for subset, name in zip([df_prostatex], ['ProstateX']):\n",
    "    #Fix patients\n",
    "    spec_df= subset.merge(patient_info, how='left', on='pid')\n",
    "    def get_worst(l):\n",
    "        l= np.array(l + [0])\n",
    "        l= l[l != 10]\n",
    "        return np.max(l)\n",
    "    spec_df['lesions_ggg']=  spec_df['class_target'].apply(get_worst)\n",
    "    spec_df['ggg']= spec_df['lesions_ggg']\n",
    "    spec_df['ClinSig']= spec_df['ggg'] >= GGG_THRESHOLD\n",
    "    spec_df.loc[spec_df['pred_score'] < CLASS_CONFIDENCE_THRESHOLD, 'det_class_thresholded']= 0\n",
    "\n",
    "    #Print results at different thresholds\n",
    "    print('  %s'%name)\n",
    "    print('    N=%3d    AUC  | Acc.   | Sens.  | Specif.'%(len(spec_df)))\n",
    "    for t in threshold_list + [*get_optimal_thresholds(spec_df.ClinSig, spec_df.pred_score)]:\n",
    "        print('    %s: %.3f |  %.3f |  %.3f |  %.3f'%('t=%.2f'%t, *composite_score(\n",
    "            spec_df.ClinSig, spec_df.pred_score, t=t)))\n",
    "        \n",
    "    #Plot AUC\n",
    "    plot_auc(spec_df.ClinSig, spec_df.pred_score, new_fig=name=='All', legend=name + ' ', annotate=name=='All')\n",
    "\n",
    "_=plt.title('Patient-level ROC Curve')\n",
    "#plt.gcf().savefig('patient_auc.png', dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute DSC scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " ProstateX | DSC: mean: 0.2705 +/- 0.2231 (N=45) | median : 0.2841 | min : 0.0001 | max : 0.7712\n"
     ]
    }
   ],
   "source": [
    "dice_threshold= 0.25\n",
    "\n",
    "dscs_prostatex= []\n",
    "for patient_results, pid in test_results_list:\n",
    "    seg_p= patient_results['seg_preds'][0] > dice_threshold\n",
    "    seg_p= np.transpose(seg_p, axes=(2,0,1))\n",
    "    seg= np.load(os.path.join(cf.pp_test_data_path,'{}_rois.npy'.format(pid)))\n",
    "    dscs_prostatex.append(gen_dsc(seg, seg_p))\n",
    "    \n",
    "metric_name= 'DSC'\n",
    "for name, metric in zip(['ProstateX'], [dscs_prostatex]):\n",
    "    print('%10s | %s: mean: %.4f +/- %.4f (N=%d) | median : %.4f | min : %.4f | max : %.4f'%\n",
    "              (name, metric_name, np.mean(metric), np.std(metric), len(metric), np.median(metric), \n",
    "               np.min(metric), np.max(metric)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compete in ProstateX challenge"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To compete in ProstateX, go to `config.py` and set `self.test_set= 'train'`. Then run `exec.py` in `test` mode, and in `analysis` mode after.\n",
    "You may want to manually rename the previous test results to avoid overwriting them.\n",
    "\n",
    "**Alternatively**, the predictions are already provided as a zip file in `train_boxes.zip`\n",
    "\n",
    "You must also set the path to the challenge ``.csv`` test file provided alongside the challenge images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "if 'train' in SUBSET:\n",
    "    #Load dataset and results\n",
    "    GGG_THRESHOLD_PROSTATEX= 2 #In prostateX clinically significant prostate cancer is considered as GGG >= 2\n",
    "    prostatex_test= pd.read_csv('../ProstateX-Findings-Test.csv') #Path to the ProstateX test data\n",
    "\n",
    "    #Keep only prostatex test patients\n",
    "    prostatex_test_results_list= [p for p in test_results_list if p[1] in prostatex_test.ProxID.values] \n",
    "\n",
    "    #Set a new column 'lesion_id' that contains the actual ID used in the model for each lesion\n",
    "    prostatex_test['lesion_id']= prostatex_test.groupby('ProxID').cumcount()\n",
    "\n",
    "    #Set the index to ProxID and fid2\n",
    "    prostatex_test= prostatex_test.set_index(['ProxID', 'lesion_id'])\n",
    "    prostatex_test['ClinSig']= 0.\n",
    "    prostatex_test['ClinSig2']= np.nan #Will be true for all since class=20 > 0\n",
    "\n",
    "    #Match lesions\n",
    "    prostatex_test= match_lesions(prostatex_test_results_list, prostatex_test, drop_nans=False, score_col='ClinSig',\n",
    "                                  cl=GGG_THRESHOLD_PROSTATEX + 1,\n",
    "                                  use_distance=USE_DISTANCE, distance_threshold=DISTANCE_THRESHOLD, \n",
    "                                  spacing=SPACING, normalize_by_distance=False)\n",
    "\n",
    "    #Save with correct formatting\n",
    "    prostatex_submission= prostatex_test.reset_index().set_index(['ProxID', 'fid']).sort_index(\n",
    "                                       ).reset_index()[['ProxID', 'fid', 'ClinSig']]\n",
    "    prostatex_submission.to_csv('prostatex_submission.csv', index=False)\n",
    "\n",
    "    #Check that the scores are not all zeros\n",
    "    prostatex_submission"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}