Master_thesis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Declaration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Declaration of variables and used libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import os\n",
    "from PIL import Image  ### conda install pillow  \n",
    "import matplotlib.pyplot as plt\n",
    "from itertools import chain\n",
    "\n",
    "import pymrmr\n",
    "import IPython\n",
    "from IPython.core.debugger import set_trace\n",
    "from sklearn import preprocessing\n",
    "from sklearn.discriminant_analysis import LinearDiscriminantAnalysis\n",
    "import SimpleITK as sitk\n",
    "from medpy.io import load\n",
    "from radiomics import featureextractor\n",
    "from sklearn.linear_model import Lasso\n",
    "# import PyWavelets\n",
    "\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.ensemble import ExtraTreesClassifier\n",
    "from sklearn.metrics import classification_report, confusion_matrix, accuracy_score\n",
    "from sklearn.model_selection import cross_val_score\n",
    "from sklearn.metrics import roc_auc_score\n",
    "import re\n",
    "\n",
    "\n",
    "import Dclr_of_class_variables as DCV\n",
    "import Dclr_of_functions as DF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Master Thesis\n",
    "\n",
    "https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6252243/\n",
    "\n",
    "This thesis deals with a classification of glioma grade in high and low agressive tumours. For calculating features is used special radiomics library called PyRadiomics. Goal of this work is calculate all the features from PyRadiomics from different weightening sequence of MRI, then select the best features and choose a best classifier, that will provide the glioma grade. Another important task of this thesis is survival prediction - predicting number of days or years that the patient will survive after GTR-Gross Rotal Resection."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Image import and ROI creation\n",
    "In order to properly calculate features for classification of gliomas, the new mask-files were created that represented specific regions in tumours - regions of interest(ROI)\n",
    "The first ROI was created  by merging the non-enhancing tumor and necrotic region\n",
    "The second ROI was created by adding the tumor region with enhancememnt to the first region\n",
    "The third ROI combined the second region with the area of edema\n",
    "Using various ROI should help to distinguish the largely varying parts of tumour between grades\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Succesful\n"
     ]
    }
   ],
   "source": [
    "BRATS=DF.Create_Image_repository();\n",
    "print('Succesful')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Succesful\n"
     ]
    }
   ],
   "source": [
    "DF.Create_ROI(BRATS);\n",
    "print('Succesful')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Radiomics\n",
    "The histogram-based features were computed from 128 bin histogram computed over the whole intensity range. \n",
    "The shape descriptors are not correlated features from radionomics library as Sphericity,Major Axis Length,Flatness...\n",
    "The texture based descriptors are calculated from given matrixes:GLCM,GLSZM,GLRLM,NGTDM,GLDM."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Fixed bin Count enabled! However, we recommend using a fixed bin Width. See http://pyradiomics.readthedocs.io/en/latest/faq.html#radiomics-fixed-bin-width for more details\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ROI_I finished\n",
      "ROI_II finished\n",
      "ROI_III finished\n"
     ]
    }
   ],
   "source": [
    "extractor2=featureextractor.RadiomicsFeatureExtractor('D:/data/testovanie.yaml')\n",
    "result_HGG_I=DF.Radionomics_structure(extractor2,BRATS,'ROI_I.nii.gz','HGG')\n",
    "result_LGG_I=DF.Radionomics_structure(extractor2,BRATS,'ROI_I.nii.gz','LGG')\n",
    "print('ROI_I finished')\n",
    "result_HGG_II=DF.Radionomics_structure(extractor2,BRATS,'ROI_II.nii.gz','HGG')\n",
    "result_LGG_II=DF.Radionomics_structure(extractor2,BRATS,'ROI_II.nii.gz','LGG')\n",
    "print('ROI_II finished')\n",
    "result_HGG_III=DF.Radionomics_structure(extractor2,BRATS,'ROI_III.nii.gz','HGG')\n",
    "result_LGG_III=DF.Radionomics_structure(extractor2,BRATS,'ROI_III.nii.gz','LGG')\n",
    "print('ROI_III finished')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Z-score normalization and saving/reading CSV file\n",
    "Feature values of the training cohort were normalized to z-scores for each feature across subjects as the different radiomics features have different units and range.\n",
    "Now is necessary to divide samples in two classes- LG and HG and save them. \n",
    "As the calculation of features is computationally expensive is better to save the features after calculation as csv file and than, while tuning the algorithm, read features from csv."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.exists(DCV.newpath):\n",
    "    os.makedirs(DCV.newpath)\n",
    "    \n",
    "count=0\n",
    "for ROI_name in range(0,len(DCV.List_of_ROI_names)-1,2):\n",
    "    DF.Saving_normalized_CSV(eval(DCV.List_of_ROI_names[ROI_name]),eval(DCV.List_of_ROI_names[ROI_name+1]),DCV.ROI_list[count])\n",
    "    count=count+1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "folder=os.listdir(DCV.newpath)\n",
    "    \n",
    "ROI_I= pd.read_csv(DF.Add_folder(DCV.newpath,folder[0]),header=[0,1])\n",
    "ROI_II= pd.read_csv(DF.Add_folder(DCV.newpath,folder[1]),header=[0,1])\n",
    "ROI_III= pd.read_csv(DF.Add_folder(DCV.newpath,folder[2]),header=[0,1])\n",
    "target=ROI_I.iloc[:,0]\n",
    "No_HG=target.value_counts()[1]\n",
    "\n",
    "if 'BRATS' not in globals():\n",
    "    BRATS=DF.Create_Image_repository();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Vyobrazenie parametrov\n",
    "The  random forest with gini impurity was used to select features useful for classifying glioma grades in the training cohort. In this moment for classification are used features calculated only from ROI I"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEiCAYAAADjxEWuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3de5xVVf3/8dc7ULxr6pQKKBSoQZnliD9/lTfKtG+JfUOFyktZ2sXMW980+5XZr/paFllqSaGgX1NJs6gsJE0r8zZ4QZFQvDJhhWEipiL4+f6x19T2zJqZM8PsOTPj+/l4zMNz9l57n7UOx/Pee62z11ZEYGZmVusVja6AmZn1Tw4IMzPLckCYmVmWA8LMzLIcEGZmluWAMDOzLAeEVU5SSBrT4DocJekPnax/r6SlklZJelNf1q1Kkt4i6YHUroMbXR8bWBwQVhdJb5X0R0lPSVoh6SZJu/fh698g6SMVvsTZwHERsUlE3LkuO+oPgVhyJnBuatdP13VnkmZKWp0CZ4WkeZJ2rikzQtKlkv4u6RlJt0l6d02ZDt+j9Br/f13rauvOAWFdkrQZ8Avgu8CWwHDgS8DzjaxXL9sBWNjoSgBIGtKLu+txuyQN7WDV1yNiE4rPwZ+BGaVttgT+AKwGxgNbA9OAH0ma3JN6WOM4IKweOwJExGURsTYino2IayNiQVsBSR+WtEjSk5LmStohtyNJwySdLekxSX+V9H1JG5bWT5J0l6SVkh6UdICkrwBvA85NR67nprI7pyPYFZIWSzq0tJ+tJM1J+7kNeG0n9VkFDAHulvRgWr6dpKskLZf0sKTjS9tMkHSzpH9IelzSuZLWT+t+l4rdnep6WK57q3wEnY6YvyfpGknPAPt29j5J2lrSL9Lrr5D0e0nt/l9ObXkN8PNUl2GpXXPSdkskfbRU/gxJV0r6H0krgaNy71mbiHgWmA3sWlp8IrAKODoi/pI+K5cBXwG+KUmd7dP6FweE1eN+YK2kWZIOlPTK8srUt/054D+BJuD3wGUd7OssisDZFRhDcRT6hbSfCcDFwGeALYC9gEci4vS0z7YuoOMkbQzMA34EvAqYCpwvaXx6nfOA54BtgQ+nv3Yi4vl0NAzwxoh4bfqy/Tlwd6rfROAESe9M5dZSfBFuDeyZ1n8i7W+v0r42iYgrOngfar2f4kt0U4oj8A7fJ+BkoJXivX41xXvfbs6ciHgt8BjwnlSX5yn+XVqB7YDJwFclTSxtNgm4kuL9v7SzCqd/g6nAktLidwBXRcSLNcVnA9unNtkA4YCwLkXESuCtFF9CPwCWp6PQV6cixwJfi4hFEbEG+Cqwa+1ZRDp6/ChwYkSsiIinU9kpqcjRwIURMS8iXoyIP0fEnzqo1rspwuOiiFgTEXcAVwGTUxfN+4AvRMQzEXEvMKsbTd4daIqIMyNidUQ8lNo9Jb0f8yPilvS6jwAXAHt3Y/85P4uIm9IX6/N0/j69QBF8O0TECxHx+6hjUjVJIyn+HT8bEc9FxF3AD4HDS8Vujoifpvf/2Q52dYqkfwBPp/2Vt98aeDyzzeOl9TZAOCCsLunL/6iIGAG8nuII9Ntp9Q7AOanL4x/ACkAUR71lTcBGwPxS2V+n5QAjgQfrrNIOwB5t+0n7+gCwTdrfUGBpqfyj3WjuDsB2Nfv+HMXROpJ2TF08f0ldMV9l3b/4ynXt6n36BsVR+7WSHpJ0ap2vsR3QFjhtHuWl/05L6drZEbEFMAp4FtiptO4JivCqtW1pvQ0QDgjrtnRUP5MiKKD4Ujk2IrYo/W0YEX+s2fQJii+U8aVym5e6eJbSwVgB7btQlgI31rzmJhHxcWA5sIYicNps340mLgUertn3phHxrrT+e8CfgLERsRlFeHTWt/4MxRc+AJK26aJ9nb5PEfF0RJwcEa8B3gOcVNNN1JFlwJaSNi0t255ioDlXj05FxGPApykODtrGkX4DvC8zJnIoxft6f737t8ZzQFiX0mDwyZJGpOcjKfqeb0lFvg+c1tb/L2lzSYfU7id1n/wAmCbpVans8FLf/gzgQ5ImSnpFWtf2E8q/Ugy4tvkFsKOkwyWtl/52l/S6iFgL/AQ4Q9JGksYBR3ajybcBKyV9VtKGkoZIer3+/bPeTYGVwKpUv4/XbF9b17uB8ZJ2lbQBcEZnL97V+yTp3ZLGpC67lRRjImu7alRELAX+CHxN0gaSdqHo1ut0rKGLfc6jCJ5j0qJpwGbADEnbpNeZCpwOfKaerjDrPxwQVo+ngT2AW9OvbG4B7qUYLCUirqYYVL08dbncCxzYwb4+S9E9cksq+xtSF0VE3AZ8iOJL5ingRoruHoBzKMYXnpT0ndRNsj9Fv/wy4C+pDsNS+eOATdLymcBF9TY2Bcx7KAaIH6Y4ov8hsHkqcgrFoPLTFF/ktQPRZwCzUvfQoRFxP8X1CL8BHqAYhO5Kh+8TMDY9XwXcDJwfETfU2bypFF1Dy4CrgS+mL/l18Q3gvyQNi4i/U4xLbADcB/wdOAk4vBsD9tZPyIFuZmY5PoMwM7MsB4SZmWU5IMzMLMsBYWZmWQ4IMzPL6mi2xgFn6623jlGjRjW6GmZmA8r8+fOfiIim3LpBExCjRo2ipaWl0dUwMxtQJHU4DY27mMzMLMsBYWZmWZUGhIqbvSxONyZpN+OkpL0k3SFpjUp3m0pz1twsaaGkBZIOq7KeZmbWXmUBkebkP49iTp5xwNQ0aVrZYxR3rfpRzfJ/AkdExHjgAODbkraoqq5mZtZelYPUE4Al6WYrSLqc4m5V97UVSDdbQdJL7j6VJjdre7xM0t8o5sL/R4X1NTOzkiq7mIbz0puPtNL+BjJdSrehXJ/6byRjZma9oMqAyN1ApVtTx0raFrgE+FDmHrdIOkZSi6SW5cuX97CaZmaWU2VAtPLSO3qNoJiDvi6SNgN+CXw+Im7JlYmI6RHRHBHNTU3Z6zzMzKyHqhyDuB0YK2k0xS0Np1DcZKVLktanuJnJxRHx4+qqWAF1dufJfsr3BDGzjMrOICJiDcVdveYCi4DZEbFQ0pmSDgJIt4hsBQ4BLpC0MG1+KLAXcJSku9LfrlXV1czM2hs0d5Rrbm6OfjHVhs8gzGwAkTQ/Ippz63wltZmZZTkgzMwsywFhZmZZDggzM8tyQJiZWZYDwszMshwQZmaW5YAwM7MsB4SZmWU5IMzMLMsBYWZmWQ4IMzPLckCYmVmWA8LMzLIcEGZmluWAMDOzLAeEmZllOSDMzCzLAWFmZlkOCDMzy3JAmJlZlgPCzMyyHBBmZpblgDAzsywHhJmZZTkgzMwsq9KAkHSApMWSlkg6NbN+L0l3SFojaXLNuiMlPZD+jqyynmZm1l5lASFpCHAecCAwDpgqaVxNsceAo4Af1Wy7JfBFYA9gAvBFSa+sqq5mZtZelWcQE4AlEfFQRKwGLgcmlQtExCMRsQB4sWbbdwLzImJFRDwJzAMOqLCuZmZWo8qAGA4sLT1vTcuq3tbMzHpBlQGhzLLozW0lHSOpRVLL8uXLu1U5MzPrXJUB0QqMLD0fASzrzW0jYnpENEdEc1NTU48ramZm7VUZELcDYyWNlrQ+MAWYU+e2c4H9Jb0yDU7vn5aZmVkfqSwgImINcBzFF/siYHZELJR0pqSDACTtLqkVOAS4QNLCtO0K4MsUIXM7cGZaZmZmfUQR9Q4L9G/Nzc3R0tLS6GqAcsMn/dwg+QyYWfdJmh8Rzbl1vpLazMyyHBBmZpblgDAzsywHhJmZZTkgzMwsywFhZmZZDggzM8tyQJiZWZYDwszMshwQZmaW5YAwM7MsB4SZmWU5IMzMLMsBYWZmWQ4IMzPLckCYmVmWA8LMzLIcEGZmluWAMDOzLAeEmZllOSDMzCzLAWFmZlkOCDMzy3JAmJlZlgPCzMyyHBBmZpZVaUBIOkDSYklLJJ2aWT9M0hVp/a2SRqXl60maJekeSYsknVZlPc3MrL3KAkLSEOA84EBgHDBV0riaYkcDT0bEGGAacFZafggwLCLeAOwGHNsWHmZm1jeqPIOYACyJiIciYjVwOTCppswkYFZ6fCUwUZKAADaWNBTYEFgNrKywrmZmVqPKgBgOLC09b03LsmUiYg3wFLAVRVg8AzwOPAacHREral9A0jGSWiS1LF++vPdbYGb2MlZlQCizLOosMwFYC2wHjAZOlvSadgUjpkdEc0Q0NzU1rWt9zcyspMqAaAVGlp6PAJZ1VCZ1J20OrADeD/w6Il6IiL8BNwHNFdbVzMxqVBkQtwNjJY2WtD4wBZhTU2YOcGR6PBm4PiKColtpPxU2Bv4P8KcK62pmZjUqC4g0pnAcMBdYBMyOiIWSzpR0UCo2A9hK0hLgJKDtp7DnAZsA91IEzUURsaCqupqZWXsqDtgHvubm5mhpaWl0NUC5YZV+bpB8Bsys+yTNj4hsF76vpDYzsywHhJmZZTkgzMwsywFhZmZZDggzM8tyQJiZWZYDwszMshwQZmaW5YAwM7MsB4SZmWU5IMzMLMsBYWZmWUMbXQEbYAbaZISeiNCsx3wGYWZmWXUHhKS3SvpQetwkaXR11TIzs0arKyAkfRH4LHBaWrQe8D9VVcrMzBqv3jOI9wIHAc8ARMQyYNOqKmVmZo1Xb0CsTveKDoB0n2gzMxvE6v0V02xJFwBbSPoo8GHgB9VVy6xB/Csts3+pKyAi4mxJ7wBWAjsBX4iIeZXWzMzMGqrLgJA0BJgbEW8HHApmZi8TXY5BRMRa4J+SNu+D+piZWT9R7xjEc8A9kuaRfskEEBHHV1IrMzNruHoD4pfpz8zMXibqHaSeJWl9YMe0aHFEvFBdtczMrNHqCghJ+wCzgEcAASMlHRkRv6uuamZm1kj1Xij3TWD/iNg7IvYC3glM62ojSQdIWixpiaRTM+uHSboirb9V0qjSul0k3SxpoaR7JG1QZ13NzKwX1BsQ60XE4rYnEXE/xXxMHUo/jz0POBAYB0yVNK6m2NHAkxExhiJwzkrbDqWY6+ljETEe2Adwl5aZWR+qNyBaJM2QtE/6+wEwv4ttJgBLIuKhiFgNXA5MqikziaLrCuBKYKIkAfsDCyLiboCI+Hv6ua2ZmfWRegPi48BC4Hjg08B9wMe62GY4sLT0vDUty5aJiDXAU8BWFIPhIWmupDsk/VfuBSQdI6lFUsvy5cvrbIqZmdWj3p+5DgXOiYhvwb+6j4Z1sU1uUpvaiWM6KjMUeCuwO/BP4DpJ8yPiupcUjJgOTAdobm72pDRmZr2o3jOI64ANS883BH7TxTatwMjS8xHAso7KpHGHzYEVafmNEfFERPwTuAZ4c511NTOzXlBvQGwQEavanqTHG3Wxze3AWEmj0zUUU4A5NWXmAEemx5OB69O04nOBXSRtlIJjb4puLTMz6yP1BsQzkv51BC+pGXi2sw3SmMJxFF/2i4DZEbFQ0pmSDkrFZgBbSVoCnAScmrZ9EvgWRcjcBdwREb6S28ysDynqmE9e0u4Uv0JaRjFGsB1wWER09UumPtPc3BwtLS2NrsbAu58AdO+eAgOtfd29X8Jgb59ZjTS+25xb1+kZhKTdJW0TEbcDOwNXAGuAXwMP93pNzcys3+iqi+kCYHV6vCfwOYqL354k/XrIzMwGp65+5jokIlakx4cB0yPiKuAqSXdVWzUzM2ukrs4ghqRfEQFMBK4vrav3GgozMxuAuvqSvwy4UdITFL9a+j2ApDEUVz2b2UAx0AbgwYPwDdZpQETEVyRdB2wLXBv//snTK4BPVV05MzNrnC67iSLilsyy+6upjplZD/kMqdfVe6GcmZm9zDggzMwsywFhZmZZDggzM8tyQJiZWZYDwszMshwQZmaW5YAwM7MsB4SZmWU5IMzMLMsBYWZmWQ4IMzPLckCYmVmWA8LMzLIcEGZmluWAMDOzLAeEmZllOSDMzCzLAWFmZlmVBoSkAyQtlrRE0qmZ9cMkXZHW3yppVM367SWtknRKlfU0M7P2KgsISUOA84ADgXHAVEnjaoodDTwZEWOAacBZNeunAb+qqo5mZtaxKs8gJgBLIuKhiFgNXA5MqikzCZiVHl8JTJQkAEkHAw8BCyuso5mZdaDKgBgOLC09b03LsmUiYg3wFLCVpI2BzwJf6uwFJB0jqUVSy/Lly3ut4mZmVm1AKLMs6izzJWBaRKzq7AUiYnpENEdEc1NTUw+raWZmOUMr3HcrMLL0fASwrIMyrZKGApsDK4A9gMmSvg5sAbwo6bmIOLfC+pqZWUmVAXE7MFbSaODPwBTg/TVl5gBHAjcDk4HrIyKAt7UVkHQGsMrhYGbWtyoLiIhYI+k4YC4wBLgwIhZKOhNoiYg5wAzgEklLKM4cplRVHzMz6x4VB+wDX3Nzc7S0tDS6GqDcsEo/153PwEBrX3c/34O5fQOtbeD29QFJ8yOiObfOV1KbmVmWA8LMzLIcEGZmluWAMDOzLAeEmZllOSDMzCzLAWFmZlkOCDMzy3JAmJlZlgPCzMyyHBBmZpblgDAzsywHhJmZZTkgzMwsywFhZmZZDggzM8tyQJiZWZYDwszMshwQZmaW5YAwM7MsB4SZmWU5IMzMLMsBYWZmWQ4IMzPLckCYmVmWA8LMzLIqDQhJB0haLGmJpFMz64dJuiKtv1XSqLT8HZLmS7on/Xe/KutpZmbtVRYQkoYA5wEHAuOAqZLG1RQ7GngyIsYA04Cz0vIngPdExBuAI4FLqqqnmZnlVXkGMQFYEhEPRcRq4HJgUk2ZScCs9PhKYKIkRcSdEbEsLV8IbCBpWIV1NTOzGlUGxHBgael5a1qWLRMRa4CngK1qyrwPuDMinq99AUnHSGqR1LJ8+fJeq7iZmVUbEMosi+6UkTSeotvp2NwLRMT0iGiOiOampqYeV9TMzNqrMiBagZGl5yOAZR2VkTQU2BxYkZ6PAK4GjoiIByusp5mZZVQZELcDYyWNlrQ+MAWYU1NmDsUgNMBk4PqICElbAL8ETouImyqso5mZdaCygEhjCscBc4FFwOyIWCjpTEkHpWIzgK0kLQFOAtp+CnscMAb4f5LuSn+vqqquZmbWniJqhwUGpubm5mhpaWl0NUC5YZV+rjufgYHWvu5+vgdz+wZa28Dt6wOS5kdEc26dr6Q2M7MsB4SZmWU5IMzMLMsBYWZmWQ4IMzPLckCYmVmWA8LMzLIcEGZmluWAMDOzLAeEmZllOSDMzCzLAWFmZlkOCDMzy3JAmJlZlgPCzMyyHBBmZpblgDAzsywHhJmZZTkgzMwsywFhZmZZDggzM8tyQJiZWZYDwszMshwQZmaW5YAwM7MsB4SZmWVVGhCSDpC0WNISSadm1g+TdEVaf6ukUaV1p6XliyW9s8p6mplZe5UFhKQhwHnAgcA4YKqkcTXFjgaejIgxwDTgrLTtOGAKMB44ADg/7c/MzPpIlWcQE4AlEfFQRKwGLgcm1ZSZBMxKj68EJkpSWn55RDwfEQ8DS9L+zMysjwytcN/DgaWl563AHh2ViYg1kp4CtkrLb6nZdnjtC0g6BjgmPV0laXHvVL1f2hp4opI9S5XstpuqaV//aBu4fT3TP9o32P/f26GjFVUGRK7lUWeZerYlIqYD07tftYFHUktENDe6HlVx+wa2wdy+wdy2rlTZxdQKjCw9HwEs66iMpKHA5sCKOrc1M7MKVRkQtwNjJY2WtD7FoPOcmjJzgCPT48nA9RERafmU9Cun0cBY4LYK62pmZjUq62JKYwrHAXOBIcCFEbFQ0plAS0TMAWYAl0haQnHmMCVtu1DSbOA+YA3wyYhYW1VdB4jB3pXm9g1sg7l9g7ltnVJxwG5mZvZSvpLazMyyHBBmZpblgOgnJJ0uaaGkBZLukrRHWt4k6QVJx9aUf0TSPansXZK+05ia12cwtS/Tlr1L9Wz7WympbWaARyRdVdp+sqSZ6fEZkk6p2f8jkrauuM611ySVy86UNHkdX+8aSVt0sv5jko5Ij4+StN06vt6gbl+jVHkdhNVJ0p7Au4E3R8Tz6cth/bT6EIqLBqcCF9Rsum9EVHMBTy8aTO3rqC0RsWupzBuAa4BvlzZtljQ+Ihb2bY27fP97+7VEMbb5rs7KRcT3S0+PAu6lhz9lH+ztaySfQfQP2wJPRMTzABHxRES0fZimAicDIyS1u5p8gBhM7eusLUjaALiU4pd3j5e2Oxv4XJ/W9N+ydU5nKmdJui39jSlts5ekP0p6qHy0Lekzkm5PR+pfSstGSVok6XzgDmBk+SxI0hGp/N2SLknLzpB0Stp3M3BpOvL/D0lXl17vHZJ+8jJvX8M4IPqHayk+dPdLOl/S3gCSRgLbRMRtwGzgsJrtflvq0jixj+vcHYOpfdm2lHwduCn9jLtsNvDmmi+pvtJZnVdGxATgXF56xrMt8FaKI/P/BpC0P8U1SROAXYHdJO2Vyu8EXBwRb4qIR9t2Imk8cDqwX0S8Efh0uWIRcSXQAnwgnYVdA7xOUlMq8iHgopd5+xrGAdEPRMQqYDeKeaWWA1dIOoriupDZqdjlFEfbZftGxK7pb1pf1be7BlP7OmkLkg4E3k5xRlRrLfAN4LTaXXb0Ur1RX+i8zsBlpf/uWdrspxHxYkTcB7w6Lds//d1JcSS9M8UXKsCjEVGeP63NfsCVbV2FEbGii7oGcAnwwdTHvyfwq5dz+xrJYxD9RLoQ8AbgBkn3UFxhPhx4taQPpGLbSRobEQ80qJo9Npjal2uLpF9SjKFMioh/drDpJRQBUR6H+DvF0WzZpsA/qq5z26pysdLj50uPVfrv1yLiJWNFKu7j8kwHLy26H3YXAT8HngN+HBFrutpgsLevUXwG0Q9I2knS2NKiXSnCe+OIGB4RoyJiFPA10tXmA8lgal8HbXkUuBD4bkTc2dG2EfECxX1PTigt/h1wkKRN0/7/E7i7N2cO6KTO8O9uvcOAm7vY1Vzgw5I2SfsdLulVXWxzHXCopK3SNltmyjxNEYoApDGdZcDngZld7H/Qt6+RfAbRP2wCfDedcq6huP/Fg8DdNeWuouiK+XJ6/ltJbV8kCyLiiL6obA8Mpvbl2jKN4ot++9LZEMC8iPhMzfYzKL4YAIiIBZLOBf4gKYC/AR/pgzofQ9H/PkzSrRQHi7VdfC8REddKeh1ws4ppqlcBH6ToPutom4WSvgLcmP4t76T4VU/ZTOD7kp4F9oyIZykG+ptSF9DLvX0N46k2zF6mJD0CNPe3nxIDpNC8MyJmrMM+HmEQt68v+AzCzPoVSfMp+vxzg/0D3kBqn88gzMwsy2cQ/UAa4LouPd2Gos9zeXp+B0Vf6t8i4vU1251C0V+9Jm3zzYi4uE8q3UOdtHVT4LG07EVgekSc05BKdiH1Nd9TWnQwMAo4JSLenSnfRDEoeVz5FzLlLpDSPocCDwOHR0Sv/pJpXfX0czpQDIbPZm/zGUQ/I+kMYFVEnJ2e70UxWHZx+X88SR8D3gscEhErJW0OHBwRsxpQ7R4pt1XStsC2EXFH+kXPfIr29LtBPEmrImKTmmX70HFAfIJigHRtROxTWv4I/w6If+1T0izg/oj4SnWtWDf1fk4HqoH62ext/plrPxcRv6O4mVKtzwGfiIiVqdxTAykcakXE4xFxR3r8NLCI4jqJwaC704nczABreyef0wFvkH82O+WAGIDSUcymEfFgo+tShXRh0puAWxtbkw5tWJoC5OrOCtYxnUht+SHARNrfntf6gQHw2exVDoiBqSdXbw4I6SKlq4AT2s6O+qFnS1OAvLeLsl1NJ9JmQ0l3UVxZvSUwr3eqar1lgHw2e5UDYgBKH85nJL2m0XXpTZLWo/gf8NKI6LczXHbTVOCoNN4wB3hjzVW/bZ5Nk7ntQDFV9Sf7rorWlUH62eySA2Lg+hpwnqTNACRtJumYBtepx1RcujoDWBQR32p0fXqDpJ3o5nQiEfEUcDxwSvpSsgYbjJ/Nejkg+jlJl1EMWu4kqVXS0WnV94DfArdLuhe4EehokriB4C3A4cB+pf79Tm/K0g9NTP9GrZJage8AtWMUV9H1lA93UkxD0q/npSrr5HM6GAyGz2aP+GeuZmaW5TMIMzPLckCYmVmWA8LMzLIcEGZmluWAMDOzLAeE9QuSTpe0UNKC9DPCPbooP1PS5B68zihJ7+/BdtnXk7Rzqu+dkl7bg/2eIGmj7m63LiQ9Iume9F7fKGmH0roRkn4m6QFJD0o6R9L6ad0+kn6R2V92uQ18DghrOEl7UkwV/eaI2AV4O7C0opcbBXQ7IDpxMPCziHhTD+fGOgHoVkBI6o1p+vdN7/UNpFugpgvCfgL8NCLGAjtS3M6z384qa9VyQFh/sC3wREQ8DxART6QbuyNpt3SUO1/S3DT18kt0VEbSGEm/kXS3pDvSEf5/A29LR/0nShoi6RuSbk9H1MembSXpXEn3Sfol0O7m9eliqROAj0j6bVr2QUm3pf1fkCbfQ9L3JLWks6QvpWXHA9tR3Hu7bftVpf1PljQzPZ4p6Vup3FmSNpZ0Yar3nZImpXLjS6+/oINpPcrKM8fuBzwXERelf4e1wInAh/v6LMf6BweE9QfXAiMl3S/pfEl7w7/mv/kuMDkidgMupOZotosylwLnRcQbgf8LPA6cCvw+TbQ3DTgaeCoidgd2Bz4qaTTFvTZ2At4AfDRt/xIRcQ3wfWBaROyr4ob3hwFvSfMqrQU+kIqfHhHNwC7A3pJ2iYjvUNxIaN+I2LeO92lH4O0RcTJwOnB9qve+wDckbQx8DDgnvX4z0NrFPg8Afpoej6e410G5jSspbpYzpo762SDjO8pZw0XEKkm7AW+j+LK7QtKpQAvwemBe0fvBEIov+bKdcmVUTIk+PCKuTq/xHEAqU7Y/sEtpfGFzYCywF3BZOopeJun6OpoyEdiNYvoTgA2Bv6V1h6qYK2soxRnTOGBBHfss+3GqT1u9D1JxV0GADYDtKc4ITpc0AvhJRDzQwb5+K+nVqX6fT8s6miV40M4ebJ1zQFi/kL74bgBukHQPcCTF0ezCiNizk02VK6M0iWEdBHwqIubWbP8uuv+lKGBWRJxWs6/RwCnA7hHxZOo22qCDfZRfs7bMMzWv9b6IWFxTZpGkW4H/AOZK+khE5MJt37S/mcCZwDABOLYAAAGLSURBVEnAQuB9NXXfDBgJPAhs1UGdbZByF5M1nKSdavrKdwUeBRYDTWkQG0nrSRpfs3m2TOoaaZV0cFo+LPWjP01xj+E2c4GPp64qJO2Yump+B0xJYxTbUnyhduU6YLKkV6V9bZl+IbQZxZfxU+mo/cDSNrX1+auk10l6BUU3V0fmAp9KA8tIelP672uAh1L31RyKLq2siHiWYgzlCElbpvpvJOmItK8hwDeBmRExkCeCtB5yQFh/sAkwKw0IL6DofjkjIlYDkykGZe8G7qJmLKCLMocDx6d9/pHipvMLgDVp4PpE4IfAfcAdKmbFvYDizPpq4AHgHoqZc2/sqhHpHsWfB65NrzmP4l7GdwN3UhyhXwjcVNpsOvCrtkFqijGSXwDX0747rezLwHrAglTvL6flhwH3qrj50M7AxV3U+XHgMuCTUczc+V7gEEkPAPcDz1Hc3tZehjybq5mZZfkMwszMshwQZmaW5YAwM7MsB4SZmWU5IMzMLMsBYWZmWQ4IMzPLckCYmVnW/wI2HlFcB5UaiAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "model = ExtraTreesClassifier(n_estimators=2000,criterion='gini')\n",
    "model.fit(ROI_I.iloc[:,1:21], target)\n",
    "importances = model.feature_importances_\n",
    "\n",
    "indic = np.argsort(importances)[::-1]\n",
    "ticks=DF.Create_ticks(ROI_I,indic)\n",
    "pom=importances[indic]\n",
    "indices=pom[0:5]\n",
    "plt.figure()\n",
    "plt.title(\"Selected features for ROI I \")\n",
    "plt.bar(np.arange(5), indices,color=\"r\")\n",
    "plt.xticks(np.arange(5),(ticks))\n",
    "plt.xlim([-1,5])\n",
    "plt.xlabel(\"Selected features ROI I\")\n",
    "plt.ylabel(\"Score\")\n",
    "\n",
    "plt.savefig('Selected_features_ROI_I.svg')\n",
    "plt.show()\n",
    "\n",
    "number_of_selected_features=7\n",
    "Indices=indic[0:number_of_selected_features]+1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 198,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAExCAYAAACwI+/+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOydd3gc1dX/P2dmi1ZdtuQi94bBNGOMCaF3CBAIIaEkhOSXlxrSCCSkkEAS0shLGklIoXdeCAQCDt04oRjb4N67ZbnI6tL2mfP7Y0byarUqNsKW7ft5nn00e+feO2dGO/Od284RVcVgMBgMhr7E2tMGGAwGg2Hfw4iLwWAwGPocIy4Gg8Fg6HOMuBgMBoOhzzHiYjAYDIY+x4iLwWAwGPocIy4Gg2GnEJHjRWT5nrbD0L8x4mLol4jIOhGJiUiziDSIyNsico2I9Oo3KyKjRURFJPBR27o3ISIniUjVTpZRERnf9l1V/6OqE/veOsO+hBEXQ3/mPFUtAkYBvwC+A9yzJwzpa5EyomfY1zHiYuj3qGqjqj4HXAxcISKHAIjIOSLygYg0ichGEbk1o9hM/2+DiLSIyDEiMk5EXheRWhHZLiKPiEhpV8f139i/IiIrgZV+2oEi8oqI1InIchH5bEb++0Xkbn9/s4i8KSKjPkR9nxCRJX5dm0Tkxox954rIvIxW3WEZ+9aJyI0iskBEGkXkCRHJE5ECYDpQ6V+TFhGpFJFpIvKOX9dmEblLREJ+XW3Xcb6f/+Ls1o+IHCQiM/zyi0Xkk1nX5I8i8oJ/HrNEZFxP/3PDPoCqmo/59LsPsA44LUf6BuBaf/sk4FC8l6TDgK3ABf6+0YACgYyy44HTgTBQgSdAv+3GBgVeAQYAEaAA2Ah8CQgAU4DtwMF+/vuBZuAE/xi/A/77IerbDBzvb5cBU/ztKcA24GjABq7wr1c449q9B1T6x1oKXJNxzaqyzvNI4GO+DaP9/N/Isnt8xvf2OoAgsAr4HhACTvGvwcSMa1IHTPPrfwR4fE//vszno/+Ylothb6Ma74GJqs5Q1YWq6qrqAuAx4MSuCqrqKlV9RVUTqloD3Nldfp+fq2qdqsaAc4F1qnqfqqZV9X3gaeCijPwvqOpMVU0A3weOEZERu1hfCpgkIsWqWu/vB7gS+IuqzlJVR1UfABJ4AtHG71W1WlXrgOeByd1cl7mq+q5vwzrgL724Lm18DCgEfqGqSVV9HfgXcGlGnn+o6nuqmsYTly5tMew7GHEx7G0Mw3sTRkSOFpE3RKRGRBqBa4DyrgqKyCARedzvYmoCHu4uv8/GjO1RwNF+90+DiDQAnwOG5Mqvqi2+rZW7WN+ngU8A6/0utmMyyn0rq9yIrONsydiO4glATkTkABH5l4hs8a/Lz+j5urRRCWxUVTcjbT3e/2mnbTHsOxhxMew1iMhReA+t//pJjwLPASNUtQS4GxB/Xy533z/30w9T1WLg8xn5uyKzno3Am6pamvEpVNVrM/K0t1JEpBCvlVW9K/Wp6mxVPR8YBDwLPJlR7vascvmq+lgP55J9/Db+DCwDJvjX5Xv0fF3aqAZGZM3iGwls6mV5wz6KERdDv0dEikXkXOBx4GFVXejvKgLqVDUuItOAyzKK1QAuMDYjrQhowRvkHwbctJOm/As4QEQuF5Gg/zlKRA7KyPMJETnOHxD/CTBLVTfmrq7r+kQkJCKfE5ESVU0BTYDjl/sbcI3fchMRKRBvckNRL85hKzBQREoy0or8+ltE5EDg2hxlxpKbWUAr8G3f/pOA8/D+V4b9GCMuhv7M8yLSjPem/n28MZIvZey/Dvixn+eH7HizR1WjwO3AW37X0ceA2/AGwxuBF4B/7IwxqtoMnAFcgvfGvgX4Jd7gfRuPAj/C6w47Eq+ba1fruxxY53dVXYPX0kJV5+CNu9wF1OMNqH+xl+ewDG9sao1/XSqBG/GEuRlPuJ7IKnYr8ICf/7NZ9SWBTwJn401G+BPwBf84hv0YUTXBwgyGvkBE7sebRfWDPW2LwbCnMS0Xg8FgMPQ5RlwMBoPB0OeYbjGDwWAw9Dmm5WIwGAyGPseIy16OZHms3UM2rBOR0/akDfsaIjLS9+Vl72lb9iXE8/12i7+90x6iDb3HiIthr2FfEjHf0eP/dLVfVTf4CyqdrvLsLqSPwxd82Id6LntE5Isi8t/uygGo6jWq+pNdPXaWHXv8xa4/Y8TFYDB0SV8JSn/AtAJ3M3vac6b5fLgPnjuPrwFr8Bax3QFY/r5xwOtArb/vEaA0o+x38Nx0NAPLgVP9dAu4GVjtl30SGJBR7nI8/1G1eIsb15HDg7GftwR4EG/F/HrgBxn2fRHPlcuv8RYDrgXO7qKeh/BW3MfwVtl/G3gA+Ja/f5h/La7zv4/HW8jYNmnlSrzFhnV4LmMquzjOaL+eL+Et3qzHW8B4FLAAaADuyshv+ee0Hs9T8YNAib8vD89/Wa1fbjYwGG9xpwPE/XO5qxs7Av73GcBPgbf9Ms8DA/3/aZNf9+he/i66s7ntuF/G80A90/+r/nFbgGPo+be1Dm9x5gK8RatP+NejwP8fuhn1dfpfAOcAH/jnthG4NWNfLnvi/jVtARr8fPfjubZ5Ec+LwGl+2k/9/ScBVXjubrb7Nn8u4zgzgP/J+P5FfC/X/nVRv94W4GI//Vxgnv//fhvP1VC399u++tnjBpjPh/wHej/wN/B8WI0EVrTdEHTjYh6Y6N+0lf730cA4f/sbwLvAcL/sX4DH/H2T/Jupza38nUCarsXlQeCfeC5GRvv2fdnf90U8z79X4rmOvxZvpbp0Ude6zOMA/w943t++DE8Mn8jY909/+xT/4THFt/kPwMwujjHav6Z34z0Mz8B7cD2L5+NrGN4D+cSM46zCc49SiLfq/yF/39V4IpDvn9+RQLG/r8ODqxs7MsVlFd5DvQRY4l/L0/Bc2T8I3NfL30V3Nrcd90E8IYhk29LTbyvjf9Vrt/85zv8kdi6cwhfJCG/gp92PJ2zH+vXk0Vlc0ni/4TCeJ+hWdoQL6PA/yj4GnUMRdBkKgW7ut331s8cNMJ8P+Q/0fuBnZXy/Dniti7wXAB/42+P9G+E0IJiVbykZb1XAUDwRCOC5WXk8Y18BkCR37BUbzxX8pIy0q4EZ/vYXgVUZ+/L98xnShf3r6Cgu4/DeEC08MbiaHXFGHgBu8LfvAX6VUa7QP5/ROY7R9uAalpFWi/9m6n9/Gj/eCfAafmvJ/z4x41r9P7LeXjPydXhwdWNHprh8P2P//wLTM76fB8zrze+iB5vbjju2K1t6+m1l/K8+n/H9V8Dd/vZJ9CAuOer/LfCbruyha3F5MEdatrgUZOx/Ergl1/8o+xh0Fpc/Az/JOt5yPNHq8n7bVz9mzGXfINMx4np81+vduZhX1VV4LZRbgW1+vjaX7aOAZzLcuS/F63IY7Ned6Va+Fe/hm4tyvABS67Psy+mOXT1/YNBLl+yquhqvFTUZOB7PEWS1iEzEu6Hf9LNWZtqgniv82iw7stmasR3L8b3Nxg51+9sBvGv1EPAS8LiIVIvIr0Qk2Jtz+5A2tZHzd9GDzbnKdqKX4Qt22dX+zoZT6IZuzwOo93/DbWRep51lFF2EQujhftsnMeKyb5AZjGokO1y8d+tiXlUfVdXj8G4KxXOaCN4NebZ2dOmep6qb8KIjZrqVz8fr+8/Fdrw34lFZ9u2qO3bNkfYmXnCtkG/fm8AX8CI3zvPzVGfaIF6434Efwo5MOtSNd35pYKuqplT1NlWdBHwcrz/+C92cS1/T1e+iS5sz0rSL7TZ2JXxBd/Vls7PhFLqqs6djlfm/hzYyr1MrXmu6jcy4PbnoNhRCN/fbPokRl32Dm0SkTLyIh19nh1fbLl3Mi8hEETlFRMJ4Ywoxdrh0vxu4Xfz47yJSISLn+/ueAs7NcCv/Y7r4Hak3jfZJv64iv74b8N5yd4Vcrt/fBK7H6/MHryvjq3jdF23n8yjwJRGZ7J/vz/Bc4a/bRTsyeQz4poiM8eO3/Axv3CctIieLyKH+LKUmPKFts6k7N/Z9RVe/iy5t7qKevg5fkMvtfzY7G05hKzDc/03uLLf5IQ6Ox3sB+D8/fR5woYjk+1OOv5zjPDJt6DIUQg/32z6JEZd9g38Cc/FuhhfwxhigexfzYeAXeK2LLXiD1d/z9/0O763xZd+d/bt4g5So6mLgK3gP7M14s6m6W7PwVbw3wDV4M8MeBe7dxfP8OfADv8vhRj/tTbwHUZu4/BfvbbPtO6r6GnAL3ljJZryxmkt20YZs7sXr/pqJN9stjnfO4L3pPoUnLEt9W9uE9XfARSJSLyK/7yNbsunqd9GdzZ3QPg5foLnd/mezs+EUXgcWA1tEZHtvbcH77dfjtVYewZt00BYu4Dd444lb8cbwHskqeysZoQi0+1AI3d1v+yTGt5jBsA8iIooXWXLVnrbFsH9iWi4Gg8Fg6HOMuBgMBoOhzzHdYgaDwWDoc0zLxWAwGAx9Tr8RFxE5S0SWi8gqEbk5x/6wiDzh758lIqP99JCI3CciC0VkvoictJtNNxgMBkMW/cLjqb8O4I94voqqgNki8pyqLsnI9mW81bTjReQSvAVIF+NN/UNVDxWRQcB0ETlKVd3ujlleXq6jR4/+CM7GYDAY9l3mzp27XVUresrXL8QFmIbnY2oNgIg8DpyP55yvjfPx5pWDt3bgLhERPEeKrwGo6jbf5cJUPKd5XTJ69GjmzJnTl+dgMBgM+zwisr7nXP2nW2wYHX0AVdHZ71N7Hn8lcSOeC4/5wPkiEhCRMXieZ0dgMBgMhj1Gf2m55PJJlD2Nras89wIHAXPwnM69jecnqfNBRK4CrgIYOXLkrtpqMBgMhh7oLy2XKjq2Noazw3lcpzx+dLwSPN9DaVX9pqpOVtXzgVJgZa6DqOpfVXWqqk6tqOixy9BgMBgMu0h/EZfZwATfkV4Iz+/Tc1l5nsMLvgOeF9zXVVV9p3IFACJyOpDOmghgMBgMht1Mv+gW8z3IXo8X+8IG7lXVxSLyY2COqj6H53TvIRFpC1Xb5nhwEPCSiLh4LtQv3/1nYDAYDIZM9tsV+lOnTlUzW8xgMBh2DhGZq6pTe8rXX7rFDAaDwbAPYcTFYDAYDH2OEReDwWAw9DlGXAwGg8HQ5xhxMRgMBkOfY8TFYDAYDH2OEReDwWAw9DlGXAwGg8HQ5xhxMRgMBkOfY8TFYDAYDH2OEReDwWAw9DlGXAwGg8HQ5xhxMRgMBkOfY8TFYDAYDH2OEReDwWAw9DlGXAwGg8HQ5xhxMRgMBkOfY8TFYDAYDH2OEReDwWAw9DlGXAwGg8HQ5xhxMRgMBkOfY8TFYDAYDH2OEReDwWAw9DlGXAwGg8HQ5xhxMRgMBkOf06/ERUTOEpHlIrJKRG7OsT8sIk/4+2eJyGg/PSgiD4jIQhFZKiLf3d22GwwGg2EH/UZcRMQG/gicDUwCLhWRSVnZvgzUq+p44DfAL/30zwBhVT0UOBK4uk14DAaDwbD76TfiAkwDVqnqGlVNAo8D52flOR94wN9+CjhVRARQoEBEAkAESAJNu8dsg8FgMGTTn8RlGLAx43uVn5Yzj6qmgUZgIJ7QtAKbgQ3Ar1W17qM22GAwGAy56U/iIjnStJd5pgEOUAmMAb4lImM7HUDkKhGZIyJzampqPqy9BoPBYOiC/iQuVcCIjO/Dgequ8vhdYCVAHXAZ8G9VTanqNuAtYGr2AVT1r6o6VVWnVlRUfASnYDAYDAboX+IyG5ggImNEJARcAjyXlec54Ap/+yLgdVVVvK6wU8SjAPgYsGw32W0wGAyGLPqNuPhjKNcDLwFLgSdVdbGI/FhEPulnuwcYKCKrgBuAtunKfwQKgUV4InWfqi7YrSdgMBgMhnbEe/Hf/5g6darOmTNnT5thMBgMexUiMldVOw07ZNNvWi4Gg8Fg2Hcw4mIwGAyGPseIi8FgMBj6HCMuBoPBYOhzjLgYDAaDoc8x4mIwGAyGPseIi8FgMBj6HCMuBoPBYOhzjLgYDAaDoc8x4mIwGAyGPseIi8FgMBj6nMCeNsBgMBj2FO9tX81vlr3ImpYaysOFXDn+FM4ffiRegFvDh8GIi8Fg2C95v24t35j7EAk3BcCWeCN3LP0XLak4nx973B62bu/HdIsZDIb9krtWvNwuLG3EnRR/W/06adfZQ1btOxhxMRgM+yVrm3OHOk+6Do2p2G62Zt/DiIvBYNgvGVEwMGd6QCyKg5HdbM2+hxEXg8GwX3LthNMIW8EOaXl2kMvHHEfQsveQVfsORlwMBsN+yTEVE/jp4Z+hMlKGAMXBCFeNO4Urx5+yp03bJzCzxQwGw37LKUMO5pQhB5NyHQJimSnIfYgRF4PBsN/T1g2WctO8unkJs7avZnCkhE+NOJIhkZI9bN3eiREXg8FgAKLpBFe89TeqovVEnSRBy+b+1f/lD0d9nmnlY/e0eXsdZszFYDAYgEfWvsP61lqiThKAlOsQd1J894P/w1V3D1u392HExWAwGIDpmxaScNOd0lvTSda05F4TY+gaIy4Gg8EAhO3cowSuuuRlTVk29IwRF4PBYAA+M2oaEbujiAjCiIIBDC8YsIes2nvpN+IiImeJyHIRWSUiN+fYHxaRJ/z9s0RktJ/+ORGZl/FxRWTy7rbfYDDs3Vww4ghOHnIQYStAnhWkIBCmPFzInUdetqdN2ysRVd3TNiAiNrACOB2oAmYDl6rqkow81wGHqeo1InIJ8ClVvTirnkOBf6pqj1M7pk6dqnPmzOnL0zAYDHuYjS0N3DHvTd7aspaiUB5fPvAoPjdhCtZOrF9Z07yNefUbKQ8X8vGK8QTMav0OiMhcVZ3aU77+MhV5GrBKVdcAiMjjwPnAkow85wO3+ttPAXeJiGhHdbwUeOyjN9dgMPQ3tsVaOG/6fTQnE7godYkYP3//DVY1bOe2aWf2up6xRYMYWzToI7R0/6C/dIsNAzZmfK/y03LmUdU00Ahke567GCMuBsM+TV1TlOffWcz095bRHI23p9+7bDbRdAqXHe+bMSfFE6vnsz3euidM3a/pLy2XXG3W7P66bvOIyNFAVFUXdXkQkauAqwBGjhy5C2YaDIY9yVP/mc+vn3wT27YQwHWV2//f2Zw8eTyzt20klSMOS8gOsKKhhvIhBbvf4P2Y/tJyqQJGZHwfDlR3lUdEAkAJUJex/xJ6aLWo6l9VdaqqTq2oqPjQRhsMht3H+q31/O//zSSZdoglUkQTKeKpNN+7dzqNrXHGFg/AzjG2knIdmpbW8u0L7uRzh93MLZfdxcr56/fAGexf9BdxmQ1MEJExIhLCE4rnsvI8B1zhb18EvN423iIiFvAZ4PHdZK/BYNjNTJ+9jLTbeaV8Kp7k0nNuJ//5bQSl4+B7yLIZJyXc/aUHWfjOSuq2NjL39cXcdP6dLJ29ZneZvl/SL8TFH0O5HngJWAo8qaqLReTHIvJJP9s9wEARWQXcAGROVz4BqGqbEGAwGPZ+WpNJ1jc2kHDSqCqboo3EIw6u3bHH3FVoaY0z9++zGPdIHZWRIkKWTciyOX34BEp+u55ELNmeXxUSsSR/v+3p3X1K+xX9ZcwFVX0ReDEr7YcZ23G81kmusjOAj32U9hkMht1DynG4beYbPLV0EbZloark5wWpj8VwRnnCEq4VCrZ478YVgxo45pJFFA5IsHl2GSc0nsixlx5PfiCIlYRPr38153FWL6rKmb6lrpln31rEltomph44gjOOPIBQsN88KvcazBUzGAz9il+8PZOnly0m4TjgOKilRKMpb0qP39eSGKhI0uWIsg2cc/a72LaLFYAhU1toit5BcfBYQnYI13IJ5QWJtyY6Hae0oqhT2uzlG/n6Xc/iuC6ptMsr76/k/n/P5oGbL6UgL/QRn/m+Rb/oFjMYDPsOqsryhVX89+VFbN1Uv1Nlk47Do4sWEE97DiS1bUJo9ji9BU55mrNPm0Uw7AkLQDCiSEkrN/zzW/zhtbdR4IIrTyYc6SgMoUiIC67qGHHSdZUf3DOdeDJNKu2N7cQSKapqGnn4lbk7dR4G03IxGAx9SENtC9/98r1s2ViHZQuppMPJ5x7O13/8KSyr53fZaCqJ4/bOa8jg0kYkx+L5YMBh/LB13PvyXJoTCb7z7fNIxFO8cP9MXBHSCG4gwL13TGfJvI3ccMcl5OWH2bCtnpZY5xZOMu3w0pzlXH3eMb2yy+BhWi4Gg6HXrKuv52vP/Yuj/3g3n7jvAZ5dvIRMJxm/vOkJNq7ZRjyWJNqSIJVM8+b0Bbz4xHsd6mmNJnjulQXc88RbzJq3FtcXlJJwHmV5eT0bopBsDeYUF4B4MkQ8lebJOQuJp9NcddtFfP++a7DCIUBIpx1SyTSzXl3Er294FIBQMEBXupYXNF6RdxbTcjEYDL2iqrGRTz30CK3JFK4qtdEot7zyKhsbGvnqscfQ1BBl0dx1OOmO04UTsRTPPfIO517qzblZsXYbX/3RE6Qdh3giTSQvyLiRFfz+R58hHA7ywxNO5qZX/00snUaQDuLlrZz0PrG1RTQdUkDZ0GasjBlkyZTNrGUHABCwLLY0NTOuYiDP3PMmyUTHeC3JRJr33lhCQ20LlQOLGTW4jFWbtuNmHDMSCnDRiYf17cXcDzAtF4PB0Cv+Mms2sVSqw4M3lkrzl/dm05JMEo8msazcDiJj/oC6qvLDO5+jJZog7j/oY/EUK9Zt5YkXvHGNcyZM5G/nXsC0yuEMKSzkjDHjuXjMIQTjghWHQDMUrYWBy9L8+88fp7k2n2QqSDwZJOVYzFo2kWUbhwPguMrQkmIAtm2qy2EZBII29TVNAPz6mnMpLy2gIC9IJBQkHAxw0uTxXHDsIX1wBfcvTMvFYDD0ijlVm0jn6DcKWBZr6+o4ZOhgSsoKqNnS2GG/HbA45tRJAGypaWJbbUunOpJJh+kzFvOFC73WzbEjRnHsiFFE043UJqsoCQ7mnNIx/Ph3L5BsTRLwh0Za6gp46ief4LLrogyomAGuEmgJEbBcgnaIr5ycR1PsL8SSg5h83DC2bKzDdTq2rNRVKkeXAzC8opQXfvZlZi3dQE1DK4eNG8qYISaWy65gxMVgMPSKkaUlrKyt7ZSedBwGFxYiItzws4u49boHSacdnLRLOC9IYXGEy671ZmaJdHYa2IZkuG5RdXlly1+Z1zAdW4I4mmL0sCkMCJRSm+xYxymnL+acI/6LRQrbUqYO28jFhy0kFSmnOLKK6sYEloQ57iqLZYuOZ+2CAtQXyXAkxOU3nE04Y5qxbVl8/ODRH/Zy7fcYcTEYDL3i6qOn8db6De3ThAHCts3xY0YzqLAQgMkfG8efnvkazz/2DtXrazn86HGc+empFBR5g/RDKkoYWlHMhuo6ModSwqEA55y8o+tpTt3zzG94CUdTOJoCYEXTbI67aTLv/34kW7d53VhuaQtfvOhNgtYOh5X5wRQTK7ZR69bQ7HplXY0CcOXv5vPyb77CwllrGDCoiM9eexrHnmXGUz4K+kWwsD2BCRZm2F9pTtWwJbaU/EAZlZFDOrQYeuLfy1dw66uv05xMoqqcdcAEfnrG6eSHej+bas2GGr7ywydIpR2SyTShYICDxg/hf3/w6faV8H9a+UUaU9s6lXUcizPL/kiwOsUPPnUHFbc4/PjsNygKJjvlbXVhq9vx3CyJcOCQF8gLju+1vYaO7G3BwgwGw0eMqjJz259Z2PA8lgRBlfxAGReOvIPi4OBe1XHWxAM444AJbGtpoSgcpiCUe9V69fZG5i6rojg/zDGHju7gPmXsyAqe+cvVzHh3BTV1LRxyQCWTJw3vIHLRdOdxGQBEuWfuLI54K4HTFKWhdCBWjo42VW9SWY4KdizMNHykGHExGPYTVjTPYFHDix26mppSCf5V9SMuG3N3j+VVlVTaYeO2BsqK8nMKi6ryuydm8n+vz8O2LCwRbNviTzddhISEgGUxdvBA8sJBzjrxYNY3LGP51pl8sHUwhw46jqDldZ+VBQ5gS2Ie2esuW2N5rK1rpXDuBtKJFIvrKkiKTURTZE5US6lFbVrA6igxAauMvIBptewOjLgYDPsJ8+ufIa3xDmmKS31yA43JakpClTnLvT97DX/87cusrK8nPihEIGgjIkybNJKfXP0JCiPh9rxvL1jLP2YsIJlyAG8cJB2CT//hkXYXLpFQgD9ccS7vVP2e5tJV3iB/kzC9+k98ccIvGVY0gTOGXsm9q7+OjYNtKa4LjmvxzuKDOHbYcEZMtlnyzgoqqeFX1cdwc+Xb3riLQkBc/lF3IHluExOKt5Fnp7CtCGAzpvzuneoGNOw6RlwMhr0AVaUlnSBihwj0wo1KLpJONGe6iE3SjeXct2RRFT/8zpO02C6xYWGwhJQ/lXfWkg384C8v8ttvfKo9/9MzFhBLpNq/uxaksvxDxpJpfvbGn5g6dRV22O+iCinqJnlw6S3cfNRjjCgcQ1H0GuY3PsfA4noaWgpYtG40ifgArrlgGqGJSabf+zpFAxNUp4v5xoYzODCvlnw7xfLYQBpSeazYOJgDCrZy5ogYlx54MmX552Fbxbt07Qw7jxEXg6Gf88qmZfxk3r/ZHm8lYNlcOnYKNx56KkGrC98nXTC++AQaaje1d4m1YRNgYHh0zjIP3jOTRCJNYrgnLJmk0g6zl2xge2Mr5SVeCOHmZk+k7NIkoTFRnLDibisg1RTuUHbshGrsvI5jH2JBym5la2wdQ/LH8M2Pnc0zi0dx93uzqYvFmDZ8GDcdfzwjSkugFO5841a+u+Q3JHBwsVgS3xFdVkRxHJsVTaM5xTqa8sLjdupaGT48RlwMhn7M7Jr1fOu9Z4g73vTftOPy2Jq5RNMpfnLkOTtV1xFlF7K88TVa0ttJawLBxpYApw+9CasLJ10b1m0HwA3k7koKBmzqGqOUFkaYOX81a1dvJXxoI/mHN4GlIJA/voXWdUU0LStrL2dZuYfbVYVEPA753rqXCw85mAsPOThn3nGTR/P14VfwuxUPk3B3zBZzFVrjYVzHpjAc5AsHHdGr62PoW4y4GAz9mLuWzmwXljbiTppn1y/gO4edRmEw3EXJzoTtQi4bczdLGxgtEF4AACAASURBVF9lfetsigKDOKzsPAaER+XM31DXyqBBxdRsayIQdUiGxFsFmYGrytDyYq781ZOsrKohUZqk7PBGJOPJIgGlYEwz0U0FpJu9SQBrVw2lpKSVQKijyDhpm23rQozq5aL4jw88gnXDqnmm6lWCVoC4kyKVDFG7fSifGDOem6eeSHmkoNfXyNB3GHExGPox61tyx0OxLYuaeEuP4hJPp7h74Xs8tXIRClw4fhLXHnomh5Wd12WZVCrNb3/4DDNfXoQVskGEvLo0qeIAamm7wOSFAnzl08cyfdZSVlTVEE+myRuTe+xGRIkMidLsi8vG9SMYMW4rA8paCAYd0mlB1eLt1w5m/vwX+Mdj12PbPY8tiQifG3Uun6w8mTWtVQwIlTAif0iP5QwfPUZcDIZ+zCGlQ9kcbey0MkNVGZrf/eC0qnLpv59gce02En7r5y8L3mNG1VqePffz2DkmBmyNNfKzB59m0ZZ12OJgtaSxbAuJhCheH8cZEiEwMMyYEeVc8YmjCEeVb9/3InHXb12poJ4v42xjCBcnsMY0MT48loQb5OUZR1I5pJYhg+qJxUOsXT+EeDzECNdh7sL1TD1sdJeOMLMpChZweOnEXuU17B6MuBgM/ZivHnwC/9m6mpizYxA+Yge5+sBjybO7XxX/1uYNLK+raRcWgITrsLqhljc3reWUEeM65L931UzuXvkGqUFp9FyFTwplDynhlS7SEmdQQTEPP35De/5/PD6L+/7yBq1DbCj0HiWJjREKpnRsbQlKMOASGhxFB0epk1qimyqAQqq3lFO9pTwjt7JtqHL9vc9TlB/mpotO4swjxwMOIibM8N6EcblvMPRjJpYM5tGTruCYQaMpCIQYWVDGLZPP5NoDe579NL9mc6fxGoDWdIr5NZtJpx3mL9zIgoUbWVRXxV9XziDpptEQkCcQEho+L6QLLKJD8qiKuLy1YA2uq8SiSe67+3US8RR5dSlwvJaKxgK0vDsATQuaEgKuRUBcRLyxFyugYCt5lTXkFXWO+gggrpByHBKJBtyGG3G2HIZuPQx3+4VoasmHup6G3Ue3LRcRyR5WU6BBd9IhmYjYwCWq+shO2mcw7Ne0xpP84dl3mbeqDiVM5YgKTj5+Yq8WAlYWFpFnB2hNd5x6nB8IkqpN8qmL78J1vQH1xmOixA9KdXrdVIT6U/OxNoXAEr775xc4eMwQrjl9GnbAhkSaUJNDuCBFoiwICqm1hThNpXzhigOpYStvbHkPCWZPO1aKK5uJL88YM1IlEJW2tZb8/tIXOKiyBkt8p5TpRWjd56B8OmKbcZX+Tk8tl7nAHP/vXOB9YJuIvCoio7Mzi0ixiHxXRO4SkTPE46vAGuCzfWu6wbBv47oup//qHv6zYh1p18Vxlbnrqzn9l38nnuzcIsnmrFEHELIDZMqQAAGxeO3PH9DSEicaTdIaTZLUNOTSKwvvFdQf+4glUixas5n311XjpJ32Ogs3JyldFaVgS4IjCsv49+3XccWk0xlhV+asFoECxyK/WbBSYMeUUKMQ9B0IHDB4OxOHbCcccDqW0xQafazHczfseboVF1Udo6pj/b9tnwrgT0AuZ0QPAROBhcD/AC8DFwHnq+r5fWy7wbBP8/TsRTTG4p3SE2mHP7/+To/lI4EgT597GZMGDCJk2YQsm4llFVyTPxkyGjMqEFgdgFx6ZYFs7Ti2E0+meXvZBsaMH4wd2PEIsZNKSUy4/gsnEwzY1Da1snGxm1O0NC2kVxcR2ehQttKlZK1iJ3e0bkYMaMRxc8lSEtLLezx3w55nlwb0VfUfIvKDHLvGquqhACLyd2A7MFJVmz+EjQbDfsnMFeu63PfOqg18sxd1jCsZyIsXfJGaWCuqyqD8Qh54+C1SqY5KEqgKYG8I4oxKge15FUbBWRPBDVuEUiAKmufiHhhj0cjljPn4QAa8UUzDjGaCgQCu63L1V0/nsCmjqKpp4PKfPUY8kcIePZDSo2q9CixPWFIbimhdZrXrjgCBViVdCIiwetsAAnauhZZhCJr4K3sDuzSgLyKFXZRtfx9SVQdY21thEZGzRGS5iKwSkZtz7A+LyBP+/lmZ3XIicpiIvCMii0VkoYjk7fRJGQy7GUddljZuYmnjJlzt/CAdMaCky7KVpTvnI6siUsCgfC+g15QjRhEOd2yNiII1Ox/7rUKcTWHcqjDpeSXo1ghOBOJl4IZc3NMaYUyCRCjJstbNrD2mjtKvV7B1aj4bjsjj8Y0rWbJxK79/+j+0RBMk0w6xVcVsfqOS2i1F1G0voHX9cG6eegVhq+O7bSgGoSbFdmBz7QDmrxmJ42TaaYFEkPxLdurcdwZXXebWLeDhdU/zYvVrNKW6cP1v6JFug4WJyA05ksuATwJ3qerfsvI7QGvbVyACRP1tVdWcd4Q/4L8COB2oAmYDl6rqkow81wGHqeo1InIJ8ClVvVhEAnhjQZer6nwRGYg36cDpfKQdmGBhhj3J/Pr1fOf9R4g5ntuSiB3il1M+x+FlO1bLt8QTfOy2P+WMPjL9xi8xcmDpLh1bVfnRT59l9py1xOMpFFBbiA20SRaJ52gye32Jq4QHt8KEBGR7inEg/Z8ySHvvmzaC1ey2L7ZMlju0HpT0ngIWhCybwmCYs7aP5603V/kelD1sEUKbmwitrcMOOHzhurWc/elN5BcA4eOQopuRwIhdOu+eSDpJblvyG6qi1cTdBCEJYonF9yd9nQOKxn4kx9wb6W2wsJ5aLkVZn0JgC/D5bGEBUFVbVYv9T5GqBjK2u3vVmgasUtU1qpoEHgeyx2jOBx7wt58CThVvyswZwAJVne/bUNuTsBgMe5KmVIyvzb6P2mQLUSdJ1ElSm2zha7Pvozm1Y4V7YV6Yv3zpU4QDO57mAUv4xWfO2mVhAW9V+63fv4CbvnkWU48czbQjx/CFS46hUAJokE7ComEXZ0IMZ3wOYQHUBQp33HKOuqT9JSmKEj0g6ZXznzZJ16ExGSNxiHLIAZWEQwEKIiFCQRurPkpgbS2qSjplce/vxnH5Wacxf+WTWGV//MiEBWD6ljfYEK0i7npTpJOaIu4m+O2Kv7G/Ruz9MHQ75qKqt+1MZTmmLmfXV9fFrmHAxozvVcDRXeVR1bSINAIDgQMAFZGXgArgcVX91c7YbTDsTl7dvAA3R3vERXl1yyI+NeKo9rSPTxjObScW8q+35rO6ajjppgD//fcyTj1oHPmRXV9UaFnCKSdN4pSTJlG1ZhuvPDmLsyoG8mSihhS+ixcFAi7upFZvrARBVbPdi3mikch4TxXx+tlcxc1XNIcgOaq8XbOet269njUbt7Oxup6qeRt46o4XiWddmng0ybsvL+CI4z/aFfj/qZlF0k11Sm9Jt7I5vpXKiJn+vDP0tM7lPugyJqiq6pez0ub6+XNN81Cgq7ZlV/l7kycAHAcchdcF95rfbHut00FErgKuAhg5cmQXphgMHy31ydYOq+bbSDppGpKt7d9jidms33I5o0bHuW6U55/rbw+dxXvz4fY/TOf2b3ds3K9rruWh1bPZ0FrPtPJRXDxmCsWh3MOPScchmkoy7+Ul/ObGx3HSDvVj8mk9eoDfwvBuYx2Uau/OSjs2tp1ltwvaHIBYhoKoEogrdgLcILnvXKA0FAFg7Ihyxo4oZ8bmRiSHuxc7YBEp2HEe8XiK1tYEZWUFvXYP0xssyd2Ro9Cl12hD1/Q0W+xfOdJGAt8gRwNZVcfsoh1VQGZ7dzhQ3UWeKn+cpQSo89PfVNXtACLyIjAF6CQuqvpX4K/gjbnsoq0Gw4diyoAx5NnB9vGWNsJ2kCkDvFvIdVvZtP1z2HYLkYw77eorpvO92wfz9tw1NLXEKS70Hrpvb1vLtW8/Qcp1SKvLrJp1PLBqFs+eeiXleYXt5ePpNLe98xr/WLkYVxUaUgwZbhHaptRPKwM762Fd5LTf6a5aJFM2oaDXBRYQi9J0EQ2LQyQyI9YrhBsVOw3BFiVRDfEhChlu+yN2kCsP6tg5Me20Q3K+ydoBm9MumkYymeZ3v3+Z115bDAj5+SGuv/40Tjl5Ui+ues+cMuhYHt3wTAf3/QDloTIGh8u7KGXoip7WuTzd9gE+AM4GrgV+QdetkA6IyDgR+b6ILOom22xggoiMEc+B0CXAc1l5ngOu8LcvAl73PQW8BBwmIvm+6JwIGB8Rhn7L5LLR7QLTRpgAE6whHJjvhRpuib9Erk4DS1yOmboEyxLe/mANv3t4Bn99+i1umvUsMSdF2p91FnfS1CWi3LV0Jls3N7B21VactMONb07nHyuX4JCgKL+J0KAUVZ8eQP2UIjRXKyBmkakbjmsTSwSJtwYo+bPL+HsDXHDoJCTlgipWQgk2KYkyIV4iIDDk+QSRjWkkpVgxF0kp5xSP5/zRHeO05Bfm8aN7ryRSGCa/MEykMEwoL8h1t3+G4eMHc8evX+T115eQTDokk2kaGqLccceLzJu3/sP/U4DTh5zAwcUTCVshAhIgzwpTGCjgholXm9DIu0CP61xE5CDg+8ARwB3ANara7fJgERkKXAxcBhwG/By4tKv8/hjK9XhCYQP3qupiEfkxMEdVnwPuAR4SkVV4LZZL/LL1InInnkAp8KKqvtDTeRkMe5Krh53N66GF/HfbIjauqSH/A4u6BTVcqr/kOz//DBOmNOA4KbIdF9u2S0FRnITl8ov7XiWRSEEB1J4W7dSXkFaXpxd/wNzr38e2LdwCiyUXBhk5eDNjhmzHVc97cTQeZmF0BDTmMHRbAK1IdejZUlfQpI1shnWxGvShhZw3IJ+Z1TU0HlxIqli8SQGukiyxyVubZOQ9KZIDLNIlFnlbHYpPaELO7vzAnnzcRB6b9zM++M9yUok0k4+fSFFpPk1NMWbOXE4q1XGuTiKR5pFH32Hy5NwxaXYGW2y+c9BXWNW8luXNqykLlTJ1wOGErO4dhBpy09OYy/8BU4FfA98EHKC4TcWzB+hF5Eo8ERkOPIm3Sv+fvZkYoKovAi9mpf0wYzsOfKaLsg8DD/d0DINhT7O2po5rH/wnNU2tWAKxqEPZvBDW1hQxHFTghn89RkmsmV9NcsiO1pJMBZm/YgxpdUn7serdZNc9vE7MJZFwECCebzGwoIXRg7djW4rtt4wKI3EOnlLFwjcO6FTetSxSNXkEBySQgJffiQdwaoO0jIKCaoeFlU1M+eQwojMacF1F2sYuLG8CQMPkUga9WUOoziVU5zWD6rY1dWlzOBLiY2cc2iGttraFYNDuJC4A1dUNXda1K4wvGsP4ol3t4Te00VPL5Si81sCNwLf8tLbXjVwD9H8E3gEuU9U5ACJixjYMBsBxXb7096eoaW7d0eEVEOqOKISFMbACRKcmiB3sUh8v4pWaiZxavoKIP4geSwdYtmY4i1eN6hAR0koJge0W6Qq3Y0d3CiLLbBLlQULb4wQbXQYNrSVgd7wlLQuKCuPkWSnibtYMNAV1bJLb89nRTSdgK9FRQs0FYRB4MboKORIkmiRvaRhR3z4RUqVBnBDY/lBGKBzgqJMP2qlrV1lZiuvm6Ca0hEMOHrZTdRl2Dz1NRR69k/VV4rUu7hSRwXitF9OmNBiAWas30ppMdQ78ZUHr8Dwi9S6xQ9z2u/JvG4/lg6YRnF6+DFtcXtt8IAv+eRBiC/GB4IYg1ATBZiicHabl+AROgevN4LIgb62F1AdomgDhAXkUrEuR53SeaguAI+Sl0zgtIdwwoGAnIJ0PyVL8LrcdgiYpJXqooqEdaWqD5rvERichFsRKQbBJsJNQe2Q5g97ZTjAcoHRgEeddfuxOXbtwOMjln/84Dz38NvG4dw4iXvrnP79zdRl2Dz11i/0c+I2qbsux75eq+p3MNH/G1p+BP4vIcLxxkW0ishR4RlW/13emGwx7Dkdd6pMtFAUihHsI2tVGXWvMd9qVhQhuUHHzNGsMX5jTOIo5jaNQBce1kGOTpJtAY0FPlFw4PLKJz46dy8DSFuasHc1L7x+KBITUYXHiowI0NOaDYzPsxGbG59dQlB/DQahLFVCbKgS8NS2J+jC2Km5IiQ9zcPMUuwnEsVGRHa0iV7HURXM9PWxwy1xUwQ1DukDJ2ypYZWHGTark6FMP5oL/dzyFxZEOxbYntjOz5r/UJeuYVDyJaQOmEshyD3PppccwaHAJjz7yNnX1UQ45ZBj/8+UTGTGi2+V1hj1ET+5fGoEG4LrsQXIReV9Vp/TqICITgd+q6tkfxti+xLh/Mewq/9r0Hn9a9SJxJ4kgnFc5jesPOJeA1f1aiOqGJj7xv/eTdLLGDRylYKsSanGp/2wrZPdM+beoq97TXV1wGkK4sSBnjVzEt454haDlYAkk0jYtbphfbTqFuARxXcFVYVt1CRceMpdghjNIV6E+lc/WZAktyRCx1jzqWiOk84VEIkBjaz7ppA0uBLYEPKeSQKRGqXi/leqLA37LJaO7DNCUoPVhUC9FUlC4Tpj/h1zepGBx4xJ+u/IPuOqQVoewFWZQuIJbJn2PsJ096mTY0/SV+5c1wDnAz/wYLZkrsno9N09VlwMH95jRYOjnvFWzhN8s/ydNqShJN03CTfF89Xv8fsXz3ZZraI7xo7unQ9TBDYATAtf2biI7BaFmRVwhsiCU0/W9qxldUhbYxSkG5DXz7SkvEbaddo8t4YBDcSDOKWUrEQHbVgK2y/ETlmNLR+eYlkBZMErUCZCWAMHCNBWDmgkGHPLzkwwpbyAYSntdYnkuw19VRryijFwS4GNFAyl0oCgSo6woSllRlIJIHFA03lFk1YaBAwpyXhdXXe5e/TeSbpK077Up4SbYEt/Ky1tf7faaGvo3PQ3oq6ouEpGjgV8Bc0TkMlVdsAvHMhPFDXs99619lUSWi5CEm+KF6vc4KjqJBx55l41baxh+8jYGHN7gDV6XHctzjzgs2Lid6BD1XKOEvRZJREOUbEiT9l/+85YEkRjEJydxCrxYK6pCp9vHVr497QVyzZcJWi6H5G9meqO3uFAESsOxTtOawROtoCiZywYjgRRR3zlYSXGU7XXFOMUu6bBFMA4ydw0LWpoIXjgUDVjtcwtCAQc7P05Djd/ayGjU3H7ZWTmvZ3WsmoTbOdxxSlO8W/se51Wek7Ocof/Tq3gu/jTgr4nI2cDzIvL7XTiWmTVm2OvZGqvPme66yo9+/09SDcoB/7OCvIo4cUuJp2DGtpeQI0K0RA+EfBdreAz8ab2xlDBpyFi2vVQLeK5NIssDFM92aRkepuGTcSSkWOKiKqgvMsUSZ3hRXZdvbC1OZneSEnOCRKxUJ79glkDCDXT43tbCEYFwME0wmKaoJI7zJRfHgeb1A7Gri7ACbof6RMCyXIIRh1TUUzILOHbkKI6ZmHsdStAKol08GkLWrvtOM+x5ehKXDj9FVZ0uIkcC9+ItjuyYWeR5couI4DmZNBj6PWnX4fWtC3hlyzxCVoBPDpvGtIEHICJMLB7Ou7WdIyGmk4o2KsXjmgkPTGBlxIx3SBMucRk4vJHaMptMN1UaUj4IruGtv9/I169/mOqaGJpW4mVBmsbaRFpdDh6xkZK8OApsjxawtGYw4ZhLQm1SKoTo6EzSVZjRNL79e9hKsy1ZSGkw1r62BcBRoSZZSDrDs6SrkHQzvwulJdEd9QcgODqJM9iCHG7LBLCDDikCBNXiiMpK/nDOuV1e68F5gykPlbM5vrmDyIStEKcOOqnLcob+T0/i8nURGaKqWwBE5AvAp4F1wAk58v+6m7q622cw9Atcdfn2vPtZ0LCu3ffXu7UruGDY0Xx14rkMDg3BG0LcsdREFWJNAfJVyB8axQp2DvxlBVyKBzVTm+7oKl/EK3/NQ4+xMdaMFNsEmxxahwWRYodph6whEHA8R8NAeX4rRw7axHuvHUhqaoDtgTAVdoJAxsScJbEKFseGtA2vE7ZSIEJVvJTBoWbCVhoXYXOsmHXRcqx232GemCQc77HgupBM2WTPUxAb7HwXdaCTP0cFJyqULGjkF5dcwFmnHtnjNf/6Adfz86W/JO4kUFxcVY4aMJWPlx/TY1lD/6UncfktcBqAiJyA51Psq8BkvBX7b2XlX6uqG/raSINhdzGrdkUHYQGIO0n+UfUOF444hjc2r6Y5nkckmMS2XVxXiKeCpEIB8vKUREMYN2VhhzsKjK1BEolQF3ecsrxmG+FgGAIWqXwLtWFoZS2W1bHryRLIC6UoL4jy0Gun8OUzXqbatQhbDuoKjW6YP246noSGCFguhYE4IcsTp7iGWJ8YiKqSUpvGeITq9eWUD2lEREk2BogV2LhYINBSXYBdmkICWRMBcLFESaqNut4EAwBSimx3GfGzNYwdOIAz/9qryaQMyRvMnZPvYFHjEhpTjUwoHM9Q495+r6cncbEzXLxcDPzVd2L5tIjMy5H/WTyPxIjI06r66b4z1WDwcF2XVR+sRV1l/JQx2HbfuUN/p2ZZJ2/FABbCu7XLiGk9BaEEKceiNRGmPBJleGkD6lo0Dy+ncVkJw87chBV0dzx0XQgHw5w5+BTur3m33QeYiDeWAmC1jWm3KYlAcV4M287Vy6xEChPMXzWOXz11EadMnsfAsiaWJwcxLzmchD+XOeXaDIk0dhpnERGCuDQ35FNTPYCaam+dSKDVpXidQ+skcJIBcIVIsBErlPTPRSkNx8gP+gtBXaG5LkSrHYY0hGa0UvhYI0ccPYmb/n7NTjl7tMXm8NJDe85o2GvoUVxEJOA7qjwVPxZKN2Uzf00mLqihz1ny7gpuu/AOYi1xEM8P1S1PfovDTugbt+vFwXxssXCyY9qL8udVzyC2S8SGSFA5cEAjETuNbSmqMOizMbbOHMjKe8Yz4rMbKBjkRZVsac7n3cVjueHYoUxhBB/UbMQKOliWLxwWOEdGIeJSkE4TObIRK9/BcpI4acEOaJYpQmu9twixqqace/9zGm6eEKlsoXhkI5alWCguQtDqOijr5nWeG/m2sY50BBKjFMrTWKk0bl2YxKYCAqWe2BaH4kSCqfYuOmxlYIXDDWPO4NzKk2k4rongTwIUFOd/yP+CYV+gJ3F5DHhTRLYDMeA/ACIyntw+VLWLbYPhQ9Pa2Mp3z/wp0eYdoYBjzXG+f+7PeWTtnygeWPShj3FW5ZE8un5mJ3FR4h2iMBYEkuQHUu3rS0RAbYehJ9dxhHMYj8wrRELegpV02rvNbn9tBneefyazNz1OXjgromMAwkc0UhBOtjuIjNtB1BJcV9unETtpoaU+QnNtATiKKIQaHOKDbVJ1IcoPaKEonEAARy3SrkXI7jwGlErZtKZDOKWp9uBgBcVxNOhgB7yFmvaQOOk1hUSXlRIZ1ULBkGR2BGTSpHl206ucW3kypRXdRTI37G/0FM/ldjyHlfcDx+mO5fwW3thLNoeLSJOINOPFWGlq+y4iXbtBNRh6wcyn3sV1Oz8o1XV54/Hs4b/e47rNxBJvk0ytZHj+QL5/8GfIs0MU2GHy7TDFgTxEOopBQSDR6UELELBswqfY2EGbdDrQLiwAqsqd8//beRDc20txZIeweClCdayEltY8nIRFOmaz/f1yVj82DivtkBqdJDo5TsuUNOkSl1ETt1AU8uwSgYDlBVPO9vfopIUVqytxilyvi04glJ/CDjvtr5tieYP1gTGtOM0BjpejCHTRzdWcbs2Zbti/6XGdi6q+myNtRRd5TSzQvZSUmyTmRCkMFHcZ7nVP07i9mVSi8/L1RCxJ0/bmXaqzoflP1DfdgedfNU0wMIETKx7i2BNvYX79WoJWgMpIGZ9793YyG+OuWqjSeTwDASeAk0MEXYXNyeYOERl3lCPngkhHLbanClnxyAiCpRbJmhSElJap6XZhAMUujFExuLHTQknFoq4lj5DtUpCXIJ4Isn7BUDa5/7+9Ow+TqyoTP/59762tq6q3dKe7s+8JISEEEhNAVtlhMKAwIKPiqD8FBhdmFGF0ZtwdcZ8RF2Z0BAYEBlBQERRQ1rAEAiEJJGRfO72vtd/7/v6oSqe3hCyVToe8n+fp56m6derec/p211v3nnPeUw6RXccLhnMMftqV2TWV/OD9F3Htyy+zPdU4oMT0+MTB3miOcHs1idK8c+X8HPdv+RUvtjwFQMQt4ZIxH2b+iOGXafbY02cRDAfwcn37ESLRMHNOG7zPpTndyb0bn2Np63omxKq5YsLJTCnNj0RKpJ6gteO75OcIpwDIZFewo/nvGVPzBxZWz+jZT2UwTkumsyeYdGQjlIeS+XkdeNQGO/AR2nK1PH5fC94IHXBfwBVhTDDGm+kE4ZK+s/wVBg1W+Tq5NJ6QI+g6ZD3FCXoD/nPDwRy+Ck6/u9E7U8A890ZhrRYfSrY76Chv73JmCKxMt9CRTPGJKZfzrTd+TtbPd+g7CEEnyEcmvW8vdmSONMPzK6oZMvdu/m9eanmanGbJaZauXAd3b7qV1Z17WpX64FFVcn6WwRKqHrVgKvPOOZZIbNfs80gszDGnHj1ocNmebOXKZ3/IXRufZlnbRv6w9RU+9vxPeKEpf+Hd3nkrqsl+7/LIZFeSzfVdOvcbcz6OK05PEsm0F6ArGeKC2Bt8duQzXFy+grPiK1kUe4bK8HpKmjR/qaL5TMeOD6WLG8ndsgpNQTIRyr/k7/wREqvjaDb/ie/gU+qmcH2PhvYyNAAZyWciHiwbcSobHPTKx/ehI9E3A3GgSwgkpU9sySSDAxI2q4J6QjYp3PzYUxxbfhTfmH09C0bMYXSkhndXz+M7x36BKfHxAytkjnh7zIr8TmZZkSHpJfjS61eT04FrfEyNH82npv3LkNbn5ZZn+N32u+jItlLixjir9mJOH3lhnyGtnufxxJ3P8MdfPI6qcu5HzuDMD86lo/tmupMPAUKs5BIqy7/AV5Y/zJ+3v4bf79t8baSC3556A1sbziaTXTGgHiJxRo38PyKhuX22d2S7uXPD4/x12StEKy9LjwAAIABJREFUVq7h65c/TSTm7VzRFw/YnAvQnQrzpZs/StIJkin0ccc2Jqh8ckt+P7OjNLyvCr9McIMe4Xqf6LIITsYlNKeD2aetYUZFPUp+GeJVLXU8sm4WXiGqOEFv0Cuc8RXNjKto7Rm+nE/TLyxZNZlkKp+lONwgRBod5s8cw7OxrSRzO28zKtHSNIFwbudTUCG9MYZbliMwIoM4yuyKMXxxzoXMqrAFuo5Ue5sV2W6LHcE6sm244g4aXJrTA5bwOaheb3uJezbfSlbzw14TXheP1P8fqj7vqX1vTznXdTn7w6dx9odPA0A1x7b60xBvMwFyZFE6u+8glV7Mi02nDwgsAK2ZLloyXZREziKTfQvoP69FCQVnDHhfWTDGNdPey9QXR1B5wr1EoruyETuSz0dZ5XgkRJk9Yz1Ll08nlMoSntlJ4Og0WhmGJ5XS5QlKlyeomjGSji6lozpCpib/r1gX7GB62Q4CjrKzj2d65Q5yEx0eXT8LcFAvvxJk7wCjChsaq0gmQ0wY2Uww4NHRHmXD2hr8tiBhzyfc4BBp9RnV7XDrDy9jWXM9X3vhLyxvrqcyXMLRgSk8s2Y9btxDc4LfGSQ4JoFTmh9RpsDrbVv56LO/4r7Tr2FczNZRMbtnweUINiJUPWjSQEGYGJs6yDsOnofr7+kJLDtl/DSPNfyW02v+ho7sDlZ1PIviMbX0BKrD+USIXe3fJO5vBlEUFwG6NUvW28i8snoea945PFYZG+pCgYZsOSVuiHD8E3R134vnNwH5fgQhhITOpLHrYapi5xJw4gPqevLFR9HZ3Dmg81wE4o4iopREMgRGpSg9pwEcRVzQ8cDfOHjfr8ZPB2hQh3mnT2bpS2vJ+IoXU2aduYZgvxnxQddnVvV2Ht94FDlfUL8QVXpNsFQPNOdQ31ZJfWsl5AQn4VLmBJn2aCeBgIuqRzQa4qavvZ8H7n2R7dva+NSxczn5744iGHTpTqWZ/92fkksH8v0xAR+nNDugoz/j5/jV2mf5lzkX7cspNkcYCy5HsKAT4vy69/PH+vvJ9Ep7HnRCnDfq0iGtS0umadDtaS/Nkubf8kzjbfjqA8pzjXexoOpSTqo6Fz/xP4Vv8Lv6EGK4tGuK82uiPNMaZFx4B18d/xyVgRQikPDLiOiV+M5UXs3MIe6/SqXbjotPQl1yub/iJBaztuVrHFN3B6XhvjPHIyVRul2XwRZeUfJXMW+tH0fsrCbi4SQ+Dkk/hIRBXUUuzMLDUXzgtZVbqZtUQzDis25uG9HgwOwAO4UDOciFEARHHdKZbD4I+Pnu9Z4KKEjSwRHh1KmT+ebv38OtrzzMU/5rNGozn135EwKPxelIhsmsXon+7lHix8ZokATuWEE6FL8rn4BysNlqOfVZ1V6/p9NpjAWXI917ai+iIlTFn+p/Q0e2jQmxaVw0+grqIkN7T70mPIotyfXkEy3mcPHJaICIE+Xpxtvwel3V+OrxYvP9zAw1EsQh39vRV0iCzKl6N5eMzfDB0nuIu7tu/YWdVrLNV7Iq8FkaM9uo13IqPJe6QHvhNpeHrwkAVjZcy4KxTyEipL0szzeuJ+3lmBc6AbKL+xzbV+jUEPHgZZxxSi2Tpj1IZTABKPWZch5sOo52ojgzM3gP59+TyXg0N3cx+5rprEq3sC1RzuTSpgFzaNJegEQmRNk2+OZJ59FZ53Pz43+mM50m9kYCH5eOOXEEB7JCQBxKgkGuP/Ek7m38K4/qS2R23v4c5ZO+sp3001X4rQHaZ3q05Dp64pMbE0Ixl6mlNWwNbOpZxGungDjMLB+97yfZHFEsuBiOrzyJ4ytPGrLjbexezV8aHqIl08CU+CxOr7mIvxl9Jb9a9y0qgi1I4euyAKXBGN2DrMzoaZaG1CrG7GY8rSNB4tHzuW7Cb8m2B4BdwSU/NSRHInEvnuav2CrcxKCTItNeE4nsat7oiHDt83f13EYsdyfxq2PeosTtQjWD4uNJOVWV36PWnUbXvIsQ7SaAjyKMCbXyodrF3LLtDPx03+lgiXSG33avRAPwbP0UxsdaccXDLXzYZ32HxzbOQHyh/Bmf/3z2T/z8lqv43bkf5lPnfIt0IoPvK5VLu2l99wgiEys4Y+ZU/mHhQmriUe597Uky/RY4E0cJzeukaVUAdegzbtRTJak5Xm9uJD7CIRSlT4AJOQGumjJ0fy/m8GTBxQyp11oXc/fmn/T0r+xIbWFJ65NcP/3bTIi6dOa0T7hI5BrZ3YSMTjkKeGLgC+JSUflTRMKo14AU5rD0piSJy65hyLub8uGpx5qO9Vzz/Ct053ZdPXXnwly8dBH3vfsYapxXaWl9mdZEJ5sbvo0Xb6Q20EGo17LCkcKKkpNDTbzQMo1IHEJdhWOPcPPpY3yhJR3jrrXvYuHI9YyKttOeKeGF+klsaqqm+hlwM0LO9Xjwd0v59D+czY8e+QK//Pn9rH5lI9ODVVx+zrmcdP6uUW5bEo34g4wIFYFQJEuujN1PSBDobg6Dl6O8Mkgil2Fu5ThuOuYCxsYqd/MmY/IsuJgh46vPA1t/0afj3sMj5SV5eNv/kPLaB3zI++QozF/vs90Vl6ymWJUpYUowk18ES0CIEI79HaGSMwFwQvPIFYb09qHKiAAEJEJOU7R5UcLSgdNvroinDt956wFqSkvY3F5GzncJiMfocDvdXglrW55jZMnvqYlkqSuBaQoZhQ1ZH+3XmDI3TTiboSlVBmcL4xYnCac93HQk3znv5/OItaajPLJlVr6aHpQ+H2DsGsFN53foecqOHR2s7drAf3TeSttlHQQuU4LhEsZPq6EtlSSVy1Ebi1MVKhuYhDPffDKZQE8fzaDRVQEE6Y7x0zMv49jaUYMUMmZwwya4iMh5wI/IJ7X4b1X9936vh4HbgXlAM3C5qm4QkYnAG8DO5QGfV9Wrh6reZu+1ZprI+Gkc8h92fiFoKD4bEqsodZ1BO5DjgSq6c12Aoqo44nB06Ww2dj2Kp1G2ewHGukkccSF8JnPL/23Xm4MLyGqWIE7PfBlVJYsP/kbGx69iU9dztHgQd1JEnUxPRmEFVqRGk1KXoJtlZtUOZoa2cfHI11EgJD5TQ2HcXh/MAcm3qNJ1afEH9gUl2qIEgx5Tp2xj7imbOL50AzE3Q31iJLevO55NySp8r/BL8IHGEJoR3PSuABEOB5jzrlF8443vk/R2XZVtTW7n80u/watLp4I6jIzF+PwJJzMiO5oGd0uf0W2q0NEaw8mBH6BvgNk5KKAQHbO+z4gSy3Rs9s2wCC4i4gK3AGcDW4CXROQhVV3Zq9jHgFZVnSoiVwDfJr/GDMBaVe07480MO1uTqxG6+/zReTh4uJQ4VQSliyx9Z8wHJMTcyouYVX4mqzuewy8MRX5y26fxNP/B2u6HaPdDlEiGmPcU01JLiIbnISJs7X4d8bNExSVcWEglrR4pPHxKeGbxu1mxI8roCZtom7CGUjdFuZskowHqs+XsyJTh4DMl3srYYCtnl75BsHC7q0QcBvva74pQ6jiDBpezpi+nozzMjLJtLCxbT9DJ72t8bDs3zHqE/950Mi1+KdmcQ2t9GalsAL8cKATkYNClqipObG433vaBVySe5oiVt9PSXMaWjg4+8+jDhBJKZEqc8qpuxFG8nIOIUju6FR0N3W0lNLWUgUh+lr5PzziFgOMwp6aWcWXl+3KqjRkewQVYAKxR1XUAInI3sAjoHVwWAV8uPL4P+LHsy2pE5pBKed08sPnHFEYN93DVx5Ew76ldxIhQlN9s/jK+evhkCUqEitBojh9xMSEnwryqRb3219bzWFCmBxuJO2kU2Np0OeHgTMaOvIvX2v6PWi9GXaCbVK+RXZ4KS1uruO/F5UCQlRsnUzG1khNmv8WI0i48FXaky9ieLmditJmA+BxbsqUnsOT3oYPfTVLIDtLPoQjqKO+uXs3kSFNPYIF8H0jQ8TindiUPNh7H6EiCmkmd+ONcwuTwZoTghTFcsGABf3vpAh5s/j1ZP9XTXaKF/YsDoVDfzvucI3StL6djeynB8jR141t6rmIEiFckqS2N8uDZN3L760v5zgtP4wYdcr7P7JE1/Py8RRizr4ZLcBkDbO71fAuwcHdlVDUnIu1AVeG1SSKyFOgAvqSqTx/k+pp9tKrjZUQGue0lML5kAnMrTkJE+OiU/+L1tkfozDYxIXYc08tOxpXggP3VROawNbEYUMYGWil1Uj2jvVSTpDLLaWj9Er66rMiMpNTJEnV2feh2+GGe3DalVzWEtnWV/G7HAkZOamfkuGYcVwlKjmBhmeCo02+SJ0oOJaj0TVEDbM65CDmEXU3e4ZXiOEJdqH1gHxD5+THjIy1cN/YJfARXfDanq3imYxq56T5MXEv3slHEYqewasObfW7H5VOYKZ4KXZ29cokVuqucnOCkXSrKEgMzOTuQddvpyHXz0WPnccXRx7CqpZnqkhLGlVXgaZa1HU/QkHqT8tBYppaeSciNDai/Mb0Nl+Cy2+7EvSizHRivqs0iMg/4rYjMUtUB68eIyCcorKY5frwl2xtKng4ynpj8SR0bndTz4VwWrOGk6g/RmfwjrV0/Y1PD9ygrWURl/IM4zq4PzeOrP8WOzUtx6GSk2z3IMOIMbd0P0JS5kJwGeCY1lkonRczJ0uWHaPHCvLoyv1hqzx+aL4Q64YaZl/IYv6Ex3YKQRQshYlu2nDIn2edDfXsuzSg3jFsIF77v8FJbFc3hGiLOBkokP4myyw/jFa4zfJxBc4MBBMSj95T4ceFmFpY6PNc5HQ3BM5FnOaNpAduS2/r+HgtRLJUM0dXVK7goOFlwUxAucSkrFXKDHDvkBGhMtTEiVEo0GOK4Qud92uvkNxuvoTvXRE6TBCTCi423smj8LVQWsiQYM5jhkhV5CzCu1/OxwLbdlRGRAFAOtKhqWlWbAVT1ZWAtMH2wg6jqrao6X1Xnjxw5sshNMHsyrfQ4fB3YBxF0wswu7ztnYkfbV9ja8hm608+QzCylof3fWd+wCM9P4qX+TKbtBqKp/2VGbDJBYUCa+Z59o0yUVzgp0sJ4N0m7H2FLrowWr4RN20ZTX1+N70IuBrnSwk8YKktifGfODfzt2PMZUzIJV/LfwZanxuDj9MkenFFledZneQYavTTNZJlS2cy7Sl5jSmQ0bX457X5JT2DxcsK61lpWdo8i6/f999PC6KzeAqJMjjTi4uUX8Jqm3PTyr0l5A2fyC5DoLN21j0L/iZvOb5k/ZgznjJ+DM8j3tJx6jI0O/J94sfG/6cxuJ1fIHp3TFGm/i7/Wf2vQ37kxOw2XK5eXgGkiMgnYClwBXNmvzEPAVcBi4FLgCVVVERlJPsh4IjIZmAasG7qqm71RGqzkvFFX8ej22/E0h6IEJMSc8pOZGNuVLj+b20pr1+0ou9LRKCkyuXUkm96H628CTQAu09Uj55bQ5YeJO31XhozgUBUIUhdowhUYF0ww20+xODuXWSOu5P43wvjO63hR+nyeawj+sHo1p82czPvHncv7x53L6s5XuWPDd8iqsC5TyahAJ1Eni4fQ7YfIEKDOTaB4CLnC+jAQ1VWMDU1kQ2cnqh5+TuhsCvHrZ08mNj7Jxce9yDGlWxAgrQFcPILObgKleHjqkpYgDZkscRUC/YZNh9wgp42bSVdDksbubpw0BLrz665EggH+8dxTqKh0+WvDqyS9TM+E0IgT5NJxpxILRAYcd13XXwvDwXtTmlJvkfG67faY2a1hEVwKfSjXAY+SH4r8S1VdISJfBZao6kPAL4A7RGQN0EI+AAGcCnxVRHLkb3dfraotQ98K83ZOrL6AKfFjeLX1KbKaZlb5iUyIHtWnvyKRfhEkAJru894SUoj3FrvSrXi4AkeHkjydjOf7Q1RxBDyFqkCgz+2rkChhN81FwXUEvIe46byPcP8bQbq8fhmhBX6/6k28v/isa29hdnUtHz9mPheM+Gdea/8TOVlDsx+lTX0ChcHULj4RGSwNfpJx4QCj4//Bb267nWWPNNGwrILUVWESm2Lcuvl8GJGmdnwz8bI0761+jckljQNu8WX8ACkN4vnC5rZKWrujTK5upv+dYweH6+e+j3+ZH+f5tZv46RPPs9ltZ9boWj511klMr6sG4CfzP8Ota//A623rKA/GuWL86Vwwun8XJyxvX0ZntpPgblaotPE0Zk9sPRczrHSlnmJL0//D164+20c6QeKD5GfJKryeidLkBahwkpRIlpyGOTqUJdxrNFao8D1KJN83on6Ybzw2nzs3z+mzPz+g5Co8XEfw1O9ZICycDRLE5UPv+gPHVG/sFUiUEB6Tg52Dpo9x3XGMrnsRgLaGDprr27jrzy9wZ9sm1BFwhUh5N6PnNFIV7uIDtS8QkHwqf9X8JM4nm6exMTGSjS2VbE/mZ8bHwmlm1tUTdD0cEarC5Xxq6tVMLZ0ysBL7oT3bxk3LbqAmuJ2aUCdur6skwWFU9FguGvfDohzLHF5sPRdzWIqFT8KRKL520/ubueKj6g66Zr2ngodDsx/HlQgN6TJmhd7qKeMW+jt2ftMWQJw0/3TaC9z3v0eRJtRT1o/nOHXiKhaOXUckkGVHVxmPvDWHTa1VpFMud79yIqNOa6E63NlTg4wEEIkCiX6tCVISOafnWUVNGRU1ZXw4HODZ9z1L87Q4uXiQ6KZOiHu0TohyZ/0JLChbz6hwO53JGM99upytb4ZIRVLU/1sZFBbh7E6HWbJxPPFQjjPHTuUHcz9Q1CuJF5pfQFG2ZSopDaQpcTIIiiKE3Thn1P1z0Y5l3pksuJhhRSTA+JF3sq3pQ3h+Kx4uUZL59C6DTFh0pISMO50yNlPhduPSwbiQ12uIcP520WAfvJpVppW3sLytNj+BEOXcmcs4dtQmQm7+qqeutIMrj13M/7x8Co2pEbQm43zlT5dy3oxWPnKcUhKoZFLZByH3Os2t16CaIX/rLoLjlFFW+hkANnX+ntWtPyfpNVAWn8ZlnziW331vI6lEvmM+9O0wdf9YQ/exKV7squGUyHEsGn0azzV+k3CyGZrTRLamSY0vQQM72yJ4XgmfPOqsnvZ151I8Xr+MpnQHsyvGM3/EVJz+C7LshW6vi6xmAYc3EqModfPZC7Ia4j21HyIerNnnfZojiwUXM6ykU3+iq/VTlImgrlNYYB6QfMqWoOY/KPOr8IYoqbqD47Nb2NjyOVRT+WsdfztbFSYEXRwFV3TQEVKBgE/wrnoiR8dIj4kS9VPMHbWJoNt35nvA8Thl4ioeaDwRAM93WbZtCvPO+xgAHS1dPPFrJZW4hrnnLKNyVBeRyGnE41fhOpWsa7+LlS3/0ZNRoC29nJr3ruELC27kqdua8XMe77nsBBacc8yAIPhfy77Hi39cyvrXN1E+tZp7RmzlhcbNuOIQD4b41sLzmVmZ/6Bf3bGV616+Fc/3SfsZIm6I6aWj+eG8/0fI2bd/9aPLZvHYjj+R9tOA0OmV0OmVEJIQs8rnvO37jbHgYoYFz0+i3g46Wq6BXlmMhV23s3L4ePlZIvgKW2U+x4SOZ+OOj/aM0NopJ7Aqq/hEGCHKpKD26eDPZuGt1+O0LYNRy9YCUDPHx7kyQP9ljx2BmtiuaVNBx+Gsyfm+jVUvreGGs7+O53mkExlKvh5h/MyZfPeJa3CdML7meKP1pz2BZSdfU3gTfscXfvaLPf5e3IDLiRfN58SL8re4LwRa0wm6shnGxMpxeuVL+9KyO+nO7TpO0svwZsdW7tv0LFdOPG2Px+lvenwGM8tm8UbHikKAgbATZl7luxgXtTli5u1ZcDFDQlXZntpExk8ztmQygcI36Y7Us2xo+SLp3GYqHJ9yJ7fb9PeQv2LxCoN+KyPH05Zez+4GB7r4NPkBWoAsOaYFfUBwJcSWdaV8+592LeUcjoaYMm0aocBWvH5jXHwfGtorAIgEAlREIly7YAGqytc/8CMSnbvyoSW7UqxfvokHfvQwV950CRmvDV8HX12yI7Pm7X5tg6oMR6kM900kuTXZQlN6wLxh0n6WP2x7eZ+Di4hw7ZTreKnlRZ5rfhZHHE6uPoXjK+btV53NkceCiznoGlLb+K91N9ORa0UKN6iuGHc102Jx3mr8OH5hgp5DbkBalJ3PBgYch7qyq3mt5XZy6vTJ+bWT1+tdG70AWzylOjiS4+tuZ2x5HRdf82ce+/VzuEGXC/7+NM7/yOm82BzkrY4/kOt1pRFwwpTJRZwyoYSTx0/g8tmzKYtE2La2ntb6tgHHzSSzPHbHU1x50yWE3PL86pCDiAXHvt2vbq/tKSDvbze/Iw4Lq05gYdUJ+7kHcySz4GIOKl99frL2a7RnW/tsv3PTLXy4LoDf60M8qQ5x9QYd0qsE86lYdGdZpb3rZ6S9LE25MkYF2vsMl/UVujTcty6EiEUXIYHxhANw2fUXctn1F/Yps7Dm05QERrCi9V4yfidV4eksrPkMtdOPGVAnx91dbgBwA/mA4kiQKeUfZG37HX1ujbkS4ajKa3fz7n03umQENZFyNiea+mwPO0EuHPO2o0aNKToLLuageqtrBSkvOWC7pzk60m8S7PXxnFQhrUIY3RVgJEoocgnNyUcI+S0ERPFVCeGT7vw+E9wZbNAq6nNQG+jAxSeHQ7NXlk+UiQ/4uFJC2B3JlPKP7bG+jrjMrbqKuVVXvW3b6ibWMGpiDZve3NInJUw4GuL8j53Z8/yoymtwCLCm/XY8TRF2q5g14h+pjb77bY+xt0SEr835O65bciueeqS9LGE3yMyycVw6zpYkNkPPJlGag8bLvExryw0E/bUk/ABPd49mcXIUO2/UnF/ZQaWsAnrnHFPi4lIXmYvjxIhEP0hWKtnY9EEq6OrTKZ8v7bKNY3kzmSCnyfyMeYkws+IKZla8l02d95HMbqGqZCGjYxfgOgNTnByIjSu38I9nfJlsKks2kyMQdJl98lF87cEbCAT7fndT9fE0jSuRgza7PZFL89eG5TSm2zmmfALHVU62mfSmqPZ2EqUFF3NQeNkVpJreD70W/8r4DouTdfy5awIhJ8yloy/ESX4ZX3dNPnSkhOr4B5hQ+a892zqSj7Kt+VOU0zVoNuFAcDadJZ9nfeefcSXC9PL3Mio6dLeCMqkMzz24hKZtLRx94nRmLpxmH+jmHctm6JtDKtv5I3oPKQYIOT4nRut5NjGJylAtx424hEzuGDa3fp2uzCu4TgV1pR+ntvTv+7wvGpqHahaV3XdOTyo9m0mlZw/62tbEep5sfIjmdD2T4jM5deRFlAUri9DKQrsiIU6/3G49GdObBRdzUPi5lQxckgdAmBaup81v4a4N/8CZdZ9hRu2de9xXwK2msvQf8Lp+MMg+I5SU/O1u37uy42Xu3PB9cprNpzNJbeCllr/w2ek3UxmyZReMOViGy3ou5h3GCUwbfDs+Xb6g+DSnN/DApptoSW8etGxvdRWfI1r+ryguuvPPVqIEQ8cTjF7Okzt+zX+u/iQ/Xn01zzbeR87P4qvPA5t/TlZ3pZf3NEfKS/Bo/d0AZP0MGT+1u8MaY/aTXbmYgyIY/zRe+jl63xrLqrAqU04Wt2ebp1leabmfs0Z99m33WVn6SfzY5aQSD+J5DYTCCwmETuJ/1n+BhtRGcoXJik823M2azqUsGns9Ca97wH4Un9Udr3LXxq+zpvMVQKmLTGbR2E9TG7HVFY0pBrtyMQeFGzqOyIhfIO4UQPApYUW6lsXJvgkPFY+m9Ia93q/jVBCNX0Vp+ecJR05lbfdSmlKbewILQE4zbEu9RWN6I8rAyZX51ebbWdP5SiGhjM+21Bp+ue5GErmBs9yNMfvOgos5aNzwyURrniBatwapeoaX09WF9eh3cQgwqmTmfh9ja2IVGR14W8vzszSkNjCzbF7PMsU7hSSAI+DTd9llT3MsbX18v+tijNnFgos56EQCxEPVzCg7g4D0njUvBJwQx4+4ZL/3XRasJijhAdsDEqQsWM3fjruWSbGZBCRExIkSkCBT4jNxBvnTz2mGpvSW/a6LMWYX63MxQ+bsUddTGRrLqy2/Je13MzY6h9NqP0npAawNMqv8FB6r/1W/QWSC64Q4quwEgk6YT0z5V5rTO2jLNlEXGU9bpp5frf/igH0FJcKYksEHIhhj9o0FFzNkHHFZUH0FC6qvKNo+I26MqyZ9i/s330x7tgGAEaHRvH/cDQSdXVc0VeFaqsK1AETdOGNKprElubqnr8bBJeLGmFN5etHqZsyRzGbom3eM9kwjIkJZsPpty2b9NH9tuIdXWx/D0xwzyhZyVu2HKS3i5Epj3oks/cvbsOBijDH7bm+Di3XoG2OMKToLLsYYY4rOgosxxpiis+BijDGm6IZVcBGR80RklYisEZEbB3k9LCL3FF5/QUQm9nt9vIh0icjnhqrOxhhjBho2wUVEXOAW4HzgaOADInJ0v2IfA1pVdSrwA+Db/V7/AfDHg11XY4wxezZsgguwAFijqutUNQPcDSzqV2YRcFvh8X3AmVJY8k9ELgbWASuGqL7GGGN2YzgFlzFA74U9thS2DVpGVXNAO1AlIjHgC8BXhqCexhhj3sZwCi6DrWDbf4bn7sp8BfiBqnbt8QAinxCRJSKypLGxcT+raYwx5u0Mp9xiW4BxvZ6PBbbtpswWEQkA5UALsBC4VERuBioAX0RSqvrj3m9W1VuBWyE/Q/+gtMIYY8ywCi4vAdNEZBKwFbgCuLJfmYeAq4DFwKXAE5rPX3PKzgIi8mWgq39gMcYYM3SGTXBR1ZyIXAc8CrjAL1V1hYh8FViiqg8BvwDuEJE15K9Yipde1xhjTNFY4kpjjDF7zRJXGmOMOWQsuBhjjCk6Cy7GGGOKzoKLMcaYorPgYowxpugsuBhjjCk6Cy7GGGOKzoKLMcaYorPgYowxpugsuBhjjCk6Cy7GGGOKzoKLMcaYorPgYowxpugsuBhjjCk6Cy7GGGMxcjGwAAAIwUlEQVSKzoKLMcaYorPgYowxpugsuBhjjCk6Cy7GGGOKzoKLMcaYorPgYowxpugsuBhjjCk6Cy7GGGOKzoKLMcaYorPgYowxpuiGVXARkfNEZJWIrBGRGwd5PSwi9xRef0FEJha2LxCRVws/r4nIJUNdd2OMMbsMm+AiIi5wC3A+cDTwARE5ul+xjwGtqjoV+AHw7cL25cB8VZ0LnAf8XEQCQ1NzY4wx/Q2b4AIsANao6jpVzQB3A4v6lVkE3FZ4fB9wpoiIqiZUNVfYHgF0SGpsjDFmUMMpuIwBNvd6vqWwbdAyhWDSDlQBiMhCEVkBvA5c3SvYGGOMGWLDKbjIINv6X4HstoyqvqCqs4B3ATeJSGTAAUQ+ISJLRGRJY2PjAVfYGGPM4IZTcNkCjOv1fCywbXdlCn0q5UBL7wKq+gbQDczufwBVvVVV56vq/JEjRxax6sYYY3obTsHlJWCaiEwSkRBwBfBQvzIPAVcVHl8KPKGqWnhPAEBEJgAzgA1DU21jjDH9DZsRVaqaE5HrgEcBF/ilqq4Qka8CS1T1IeAXwB0isob8FcsVhbefDNwoIlnAB65V1aahb4UxxhgAUT0yB1bNnz9flyxZcqirYYwxhxUReVlV579dueF0W8wYY8w7hAUXY4wxRWfBxRhjTNFZcDHGGFN0FlyMMcYU3RE7WkxEGoGNB/kw1cCRNiT6SGvzkdZeOPLafKS1F/bc5gmq+raz0I/Y4DIURGTJ3gzZeyc50tp8pLUXjrw2H2ntheK02W6LGWOMKToLLsYYY4rOgsvBdeuhrsAhcKS1+UhrLxx5bT7S2gtFaLP1uRhjjCk6u3IxxhhTdBZcjDHGFJ0Fl/0kIueJyCoRWSMiNw7yelhE7im8/oKITCxsXyAirxZ+XhORS4a67vtjf9vb6/XxItIlIp8bqjofqAM4xxNFJNnrPP9sqOu+Pw7kHIvIHBFZLCIrROT1wVaCHY4O4Bz/Xa/z+6qI+CIyd6jrv68OoL1BEbmtcG7fEJGb3vZgqmo/+/hDfr2ZtcBkIAS8Bhzdr8y1wM8Kj68A7ik8jgKBwuNRQMPO58P150Da2+v1+4H/Az53qNszBOd4IrD8ULdhCNsbAJYBxxaeVwHuoW7TwWxzvzLHAOsOdXsO8jm+Eri78DhKfjHGiXs6nl257J8FwBpVXaeqGeBuYFG/MouA2wqP7wPOFBFR1YSq5grbI8DhMKJiv9sLICIXA+uAFUNU32I4oDYfhg6kvecAy1T1NQBVbVZVb4jqfSCKdY4/APz6oNa0OA6kvQrECiv+lgAZoGNPB7Pgsn/GAJt7Pd9S2DZomUIwaSf/jQ4RWSgiK4DXgat7BZvhar/bKyIx4AvAV4agnsV0QOcYmCQiS0XkSRE55WBXtggOpL3TARWRR0XkFRG5YQjqWwwHeo53upzDI7gcSHvvA7qB7cAm4Luq2rKngw2bZY4PM4N9O+1/BbLbMqr6AjBLRGYCt4nIH1U1VeQ6FtOBtPcrwA9Utesw+1J/IG3eDoxX1WYRmQf8VkRmqeoev+kdYgfS3gD5pcbfBSSAxwurFT5e3CoW3QH9H0P+iyKQUNXlxazYQXIg7V0AeMBooBJ4WkQeU9V1uzuYXbnsny3AuF7PxwLbdlemcClZDvSJ9Kr6BvlvA7MPWk2L40DauxC4WUQ2AJ8F/llErjvYFS6C/W6zqqZVtRlAVV8mf597+kGv8YE5kHO8BXhSVZtUNQE8DBx/0Gt84Irxf3wFh8dVCxxYe68EHlHVrKo2AM8Ce8w9ZsFl/7wETBORSSISIv8H9lC/Mg8BVxUeXwo8oapaeE8AQEQmADPId44NZ/vdXlU9RVUnqupE4IfAN1X1x0NV8QNwIOd4pIi4ACIyGZhGvs9pONvv9gKPAnNEJFr42z4NWDlE9T4QB9JmRMQBLiPfd3E4OJD2bgLeI3kx4ATgzT0e7VCPYDhcf4ALgNXkv5V+sbDtq8B7C48j5EdHrQFeBCYXtn+IfMf2q8ArwMWHui0Hs7399vFlDpPRYgd4jt9fOMevFc7xRYe6LQf7HAMfLLR5OXDzoW7LELX5dOD5Q92GoWgvEC9sX0H+i8Pn3+5Ylv7FGGNM0dltMWOMMUVnwcUYY0zRWXAxxhhTdBZcjDHGFJ0FF2OMMUVnwcWYISYiXyxkD15WyKi7sLB9pIhkReST/cpvKGSj3ZmB9z8OTc2N2XuW/sWYISQiJwJ/AxyvqmkRqSafoRbyE/KeJ58I8ef93nqGqjYNXU2NOTB25WLM0BoFNKlqGkDzKVN2puD4APBPwFgR6Z9Q0JjDigUXY4bWn4BxIrJaRH4iIqcBiMg4oE5VXwTuJZ9pt7e/9Lotdv0Q19mYfWYz9I0ZYoW8Y6cAZwCfBG4ERgIVqvpFEZkD/EJV31UovwGYb7fFzOHEgosxh5CIXEo+UeAYoBbIFl4aDcxS1bcsuJjDkd0WM2YIicgMEZnWa9Nc8gNrYqo6RndlkP4W+ay1xhyW7MrFmCFUWDzsP4EKIEc+++w2oEtVb+xVbg75NcuPLly5dJJfrAnySwp/eEgrbsw+suBijDGm6Oy2mDHGmKKz4GKMMabobIa+MUNIRKqAxwtP68j3ozQWnr9CfvZ+g6rO7ve+zwEfJ99P4wHfU9Xbh6TSxuwHCy7GDCFVbSY/QgwR+TL5jvzvFp6fCvwY6BM0RORq4Gxggap2iEg5cPFQ1tuYfWXBxZhhQlWfEpGJg7z0z+Rzi3UUyrUDtw1h1YzZZ9bnYswwJiKlQKmqrj3UdTFmX1hwMWZ4E8DmC5jDjgUXY4axwq2wbhGZfKjrYsy+sOBizPD3LeAWESkDEJEyEfnEIa6TMXtkwcWYYUJEfg0sBmaIyBYR+VjhpZ8CfwFeEpHlwJNA4hBV05i9YulfjDHGFJ1duRhjjCk6Cy7GGGOKzoKLMcaYorPgYowxpugsuBhjjCk6Cy7GGGOKzoKLMcaYorPgYowxpuj+PywXeBFFG0x5AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(ROI_I.iloc[:,Indices[0]].values,ROI_I.iloc[:,Indices[1]].values,c=ROI_I.iloc[:,Indices[1]].axes[0])\n",
    "plt.title(\"Data representation  \\n based on two most important attributes\")\n",
    "plt.xlabel(ticks[0])\n",
    "plt.ylabel(ticks[1])\n",
    "plt.savefig('graph_ROI_I.svg')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Training\n",
    "According to information provided in multiple scientific articles dedicated to problem of automatic glioma grading, the best classifier for this task is random forest, as there are not enough data to train neural network.\n",
    "To correctly evaluate classification, the accuracy in training cohort is calculated and then, using k-fold cross validation, the accuracy of testing cohort. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "X_train, X_test, y_train, y_test,idx1,idx2 = train_test_split(ROI_I.iloc[:,Indices],target,np.arange(len(target)), test_size=0.3)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy of training cohort: 1.0\n",
      "Confusion Matrix: [[19 11]\n",
      " [ 3 68]]\n",
      "Area under curve for random forest: 0.9204225352112676\n",
      "Accuracy_score of testing cohort: 0.8613861386138614\n",
      "Cross_val_accuracy: 0.8804812834224599\n"
     ]
    }
   ],
   "source": [
    "\n",
    "clasifier= RandomForestClassifier(n_estimators=200,bootstrap = True,)\n",
    "clasifier.fit(X_train, y_train)\n",
    "print(f'Accuracy of training cohort: {clasifier.score(X_train, y_train)}')\n",
    "\n",
    "y_pred = clasifier.predict(X_test)\n",
    "rf_probs = clasifier.predict_proba(X_test)[:, 1]\n",
    "\n",
    "print(f'Confusion Matrix: {confusion_matrix(y_test,y_pred)}')\n",
    "print(f'Area under curve for random forest: {roc_auc_score(y_test, rf_probs)}')\n",
    "print(f'Accuracy_score of testing cohort: {accuracy_score(y_test, y_pred)}')\n",
    "\n",
    "scores = cross_val_score(clasifier,ROI_I.iloc[:,Indices],target,cv=10)\n",
    "print(f'Cross_val_accuracy: {np.mean(scores)}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.85294118 0.94117647 0.88235294 0.88235294 0.88235294 0.78787879\n",
      " 0.90909091 0.84848485 0.87878788 0.93939394]\n"
     ]
    }
   ],
   "source": [
    "print(scores)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Errors evaluation\n",
    "Detailed analysis of incorectly clasified samples is performed in different open-source software as ImageJ or VV4\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Missclasified: [310 314 271 324 175 106 303 301   1 333 276 296 293 284]\n"
     ]
    }
   ],
   "source": [
    "y = 0 \n",
    "Output = []  \n",
    "for x in y_pred: \n",
    "    if x != y_test.values[y]: \n",
    "        Output.append(y) \n",
    "    y = y + 1\n",
    "print(f'Missclasified: {idx2[Output]}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Missclasified LGG: [49, 53, 10, 63, 42, 40, 72, 15, 35, 32, 23]\n",
      "Missclasified HGG: [175, 106, 1]\n"
     ]
    }
   ],
   "source": [
    "misclass_LGG=[]\n",
    "misclass_HGG=[]\n",
    "for x in idx2[Output]:\n",
    "    if (No_HG-1)< x:\n",
    "         misclass_LGG.append(x-No_HG-1)\n",
    "    else :\n",
    "        misclass_HGG.append(x)\n",
    "print(f'Missclasified LGG: {misclass_LGG}')\n",
    "print(f'Missclasified HGG: {misclass_HGG}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BAD HGG IMAGE\n",
      "Brats18_2013_11_1\n",
      "Brats18_CBICA_AXO_1\n",
      "Brats18_TCIA02_151_1\n",
      "BAD LGG IMAGE\n",
      "Brats18_TCIA09_141_1\n",
      "Brats18_TCIA09_402_1\n",
      "Brats18_TCIA10_130_1\n",
      "Brats18_TCIA10_299_1\n",
      "Brats18_TCIA10_325_1\n",
      "Brats18_TCIA10_393_1\n",
      "Brats18_TCIA10_410_1\n",
      "Brats18_TCIA10_628_1\n",
      "Brats18_TCIA10_639_1\n",
      "Brats18_TCIA13_618_1\n",
      "Brats18_TCIA13_650_1\n"
     ]
    }
   ],
   "source": [
    "print('BAD HGG IMAGE')\n",
    "for bad_HGG_image in range(No_HG):\n",
    "    if bad_HGG_image in misclass_HGG:\n",
    "        print(BRATS[0].type[0].grade[bad_HGG_image].foldername)\n",
    "        \n",
    "        \n",
    "print('BAD LGG IMAGE')\n",
    "for bad_LGG_image in range(No_HG):\n",
    "    if bad_LGG_image in misclass_LGG:\n",
    "        print(BRATS[0].type[1].grade[bad_LGG_image].foldername)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Survival prediction\n",
    "All patients in this dataset died in scope of few years except 2 subjects. The goal is to estimate number of years this subjects have based on age and calculated radiomic features.\n",
    "One type of analysis is classifying subject in three groups based on numbers of years they lived after being diagnosed with GBM - Glioblastoma multiforme.\n",
    "Second type of analysis will be prediction values."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Classification according to number of years\n",
    "Separation of data into 3 classes and joining two different datasets in order to combine given and calculated information are necessary steps in preprocessing data for classification."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "survival_table = pd.read_csv(\"D:\\\\Data\\\\Nove_datasety\\\\Survival_data\\\\survival_data.csv\",sep=',',converters={\"Survival\" : lambda x: 0 if x == \"NA\" else int(x)},index_col =\"ResectionStatus\")\n",
    "valid_survivors=survival_table.loc[\"GTR\"] \n",
    "survival_class=[]\n",
    "\n",
    "for j in range(len(valid_survivors)):\n",
    "    if valid_survivors.Survival.iloc[j] <=365:\n",
    "        survival_class.append(0) # less than 1 year\n",
    "    if (valid_survivors.Survival.iloc[j] >365 and valid_survivors.Survival.iloc[j]<=365*2):\n",
    "        survival_class.append(1) # less than 2 year but more than 1 year\n",
    "    if valid_survivors.Survival.iloc[j] >365*2:\n",
    "        survival_class.append(2) # more than 3 year\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "newDF = pd.DataFrame(columns=ROI_I.columns[Indices].values, index=valid_survivors.Brats18ID.values)\n",
    "\n",
    "for j in range(len(BRATS[0].type[0].grade)):\n",
    "    if BRATS[1].type[0].grade[j].foldername in valid_survivors.Brats18ID.values:\n",
    "        newDF.loc[BRATS[1].type[0].grade[j].foldername] = ROI_I.iloc[j,Indices].values\n",
    "\n",
    "survivors_table=valid_survivors.set_index('Brats18ID').join(newDF);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy of training cohort: 1.0\n",
      "Confusion Matrix: [[16  0  0]\n",
      " [ 0 13  0]\n",
      " [ 0  1  2]]\n",
      "Accuracy_score of testing cohort: 0.96875\n",
      "Cross_val_accuracy: 0.9708589655958078\n"
     ]
    }
   ],
   "source": [
    "X_train, X_test, y_train, y_test = train_test_split(survivors_table.values,survival_class,test_size=0.3)\n",
    "clasifier= RandomForestClassifier(n_estimators=200,bootstrap = True,)\n",
    "clasifier.fit(X_train, y_train)\n",
    "print(f'Accuracy of training cohort: {clasifier.score(X_train, y_train)}')\n",
    "\n",
    "y_pred = clasifier.predict(X_test)\n",
    "rf_probs = clasifier.predict_proba(X_test)[:, 1]\n",
    "\n",
    "print(f'Confusion Matrix: {confusion_matrix(y_test,y_pred)}')\n",
    "print(f'Accuracy_score of testing cohort: {accuracy_score(y_test, y_pred)}')\n",
    "\n",
    "scores = cross_val_score(clasifier,survivors_table.values,survival_class,cv=5)\n",
    "print(f'Cross_val_accuracy: {np.mean(scores)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Regresion\n",
    "As we have too little subject for accurate regresion, we have to regularization to obtain correct results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RMS: 0.26407783680148583\n"
     ]
    }
   ],
   "source": [
    "alpha=1\n",
    "regr=Lasso(alpha)\n",
    "regr.fit(X_train,y_train)\n",
    "pred=regr.predict(X_test) #predikovany y\n",
    "w=regr.coef_\n",
    "print(f'RMS: {np.sqrt(np.mean((y_test-pred)**2))}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.3"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "oldHeight": 499.44999999999993,
   "position": {
    "height": "521.05px",
    "left": "417px",
    "right": "20px",
    "top": "83px",
    "width": "671.2px"
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "varInspector_section_display": "block",
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}