Train_valid_CV_3DCNN.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "indonesian-asbestos",
   "metadata": {},
   "outputs": [],
   "source": [
    "import tensorflow as tf\n",
    "tf.keras.backend.clear_session()\n",
    "from tensorflow.keras.models import Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "rotary-roman",
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import absolute_import\n",
    "from __future__ import division\n",
    "from __future__ import print_function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "vertical-click",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "from sklearn.model_selection import RepeatedStratifiedKFold\n",
    "from pathlib import Path \n",
    "import copy\n",
    "import pandas as pd\n",
    "pd.set_option('display.max_rows', None)\n",
    "from matplotlib import pyplot as plt\n",
    "import nibabel as nib\n",
    "from scipy import ndimage\n",
    "from tensorflow import keras\n",
    "from tensorflow.keras import layers\n",
    "from sklearn import metrics\n",
    "import h5py\n",
    "from lifelines.utils import concordance_index\n",
    "from lifelines.statistics import logrank_test\n",
    "from lifelines import KaplanMeierFitter\n",
    "import random as python_random"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "turned-dragon",
   "metadata": {},
   "outputs": [],
   "source": [
    "import DenseNet3D\n",
    "import ResNet3D\n",
    "import Vgg3D\n",
    "import util\n",
    "import Aug3D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "binding-volleyball",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(123)\n",
    "python_random.seed(123)\n",
    "tf.random.set_seed(1234)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "tracked-barbados",
   "metadata": {},
   "outputs": [],
   "source": [
    "Exp_dir = Path(r'D:\\Transcend_(E)\\_IramS\\PETra_ML_DL_Analysis\\Deep_learning\\Experiments_Cox_PET')\n",
    "Exp_dir = os.path.join(Exp_dir,'Test')  # directory to store outputs\n",
    "if not os.path.exists(Exp_dir):\n",
    "    os.mkdir(Exp_dir)\n",
    "os.chdir(Exp_dir)\n",
    "# Path to directory to store weights at the end of each cv fold \n",
    "weights= os.path.join(Exp_dir,'Weights')\n",
    "if not os.path.exists(weights):\n",
    "    os.mkdir(weights)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "wired-decision",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a data frame that contains endpoint information and patient id. We need patient id for:\n",
    "# 1- startification, 2- To keep track of included patients in train and valid folds\n",
    "data_org= pd.read_csv(r'D:\\Transcend_(E)\\_IramS\\PETra_ML_DL_Analysis\\Processed_Features\\PET_Features_Training\\PET_PETra_Clinical_Imaging_Features_Robust.csv')\n",
    "\n",
    "# Provide the path to folder with 3D patches stored in nifti format\n",
    "data_org['Path']= 'D:\\Transcend_(E)\\_IramS\\PETra_ML_DL_Analysis\\Deep_learning\\Data\\PETra_Cox_model_data'+'\\\\'+ data_org.subject_id+'_0_base_img.nii.gz'\n",
    "data_org = data_org[['subject_id','LC','LCtime','Path']] # LC = event status, LCtime = time to event \n",
    "data= data_org.rename(columns={'LC': 'status', 'LCtime': 'time'}) # rename the columns time and status column for easy interprestation in remaining script \n",
    "X=data[['subject_id']] # get pateint ids for startification \n",
    "y=data[['status']] # get labels for startification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "cognitive-advocacy",
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Image pre-processing \n",
    "'''\n",
    "\n",
    "def read_nifti_file(filepath):\n",
    "    \"\"\"Read and load volume\"\"\"\n",
    "    # Read file\n",
    "    scan = nib.load(filepath)\n",
    "    # Get raw data\n",
    "    scan = scan.get_fdata()\n",
    "    return scan\n",
    "\n",
    "def normalize(volume):\n",
    "    \"\"\"Normalize the volume for PET data\"\"\"\n",
    "    volume[volume > 10] = 10\n",
    "    volume = volume/10\n",
    "    volume = volume.astype(\"float32\")\n",
    "    return volume\n",
    "\n",
    "def process_scan(path):\n",
    "    \"\"\"Read and resize volume\"\"\"\n",
    "    # Read scan\n",
    "    volume = read_nifti_file(path)\n",
    "    # Normalize\n",
    "    volume = normalize(volume)    \n",
    "    return volume\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "hairy-representative",
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "pre-processing the batch data on fly using data augmentation and adding channel\n",
    "data augmentation using batchgenerator library is imported from Aug3D.py\n",
    "'''\n",
    "\n",
    "def train_preprocessing(volume, label):\n",
    "    \"\"\"Process training data by augmenting and adding a channel.\"\"\"\n",
    "    volume = Aug3D.augment_data(volume)\n",
    "    volume = tf.squeeze(volume)\n",
    "    volume = tf.expand_dims(volume, axis=3)\n",
    "    volume = tf.reshape(volume, shape=(60, 60,44,1))\n",
    "    return volume, label\n",
    "\n",
    "\n",
    "def validation_preprocessing(volume, label):\n",
    "    \"\"\"Process validation data by only adding a channel.\"\"\"\n",
    "    volume = tf.expand_dims(volume, axis=3)\n",
    "    return volume, label"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "august-terminology",
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "tf data_loader function peform train and validation augmentation on fly\n",
    "'''\n",
    "def data_loaders(x_train,y_train,x_val,y_val):\n",
    "    train_loader = tf.data.Dataset.from_tensor_slices((x_train, y_train))\n",
    "    validation_loader = tf.data.Dataset.from_tensor_slices((x_val, y_val))\n",
    "\n",
    "    batch_size = 16\n",
    "    # Augment the on the fly during training.\n",
    "    train_dataset = (\n",
    "        train_loader.shuffle(len(x_train))\n",
    "        .map(train_preprocessing)\n",
    "        .batch(batch_size)\n",
    "        .prefetch(2)\n",
    "    )\n",
    "\n",
    "    # Only rescale.\n",
    "    validation_dataset = (\n",
    "        validation_loader.shuffle(len(x_val))\n",
    "        .map(validation_preprocessing)\n",
    "        .batch(batch_size)\n",
    "        .prefetch(2)\n",
    "    )\n",
    "    \n",
    "    return train_dataset,validation_dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "verbal-situation",
   "metadata": {},
   "outputs": [],
   "source": [
    "def cox_log_likelihood(labels, risk):\n",
    "    \"\"\"\n",
    "    The negative of the cox partial log-likelihood (since in DL one\n",
    "    typically minimizes a loss instead of maximizing likelihood).\n",
    "    Parameters\n",
    "    ----------\n",
    "    labels: 2d numpy array (n_samples x 2)\n",
    "        first column has events status, second column time\n",
    "    This does not assume that data is ordered by decreasing\n",
    "    survival time but does it beforehand.\n",
    "    \n",
    "    Method taken from: https://github.com/oncoray/cnn-hnscc\n",
    "    \"\"\"\n",
    "  \n",
    "    e = labels[:, 0] # first column in my data is event status\n",
    "    t = labels[:, 1]\n",
    "\n",
    "    # sort data\n",
    "    t_ordered, idx = tf.nn.top_k(t, k=tf.size(t))\n",
    "    e_ordered = tf.gather(e, indices=idx, axis=0)\n",
    "    risk_ordered = tf.gather(risk, indices=idx, axis=0)\n",
    "\n",
    "    # compute likelihood\n",
    "    sum_risk = tf.reduce_sum(e_ordered * risk_ordered)\n",
    "    log_sums = tf.math.log(tf.cumsum(tf.exp(risk_ordered)))\n",
    "    log_sums = tf.reduce_sum(e_ordered * log_sums)\n",
    "\n",
    "    lcpl = sum_risk - log_sums\n",
    "    return lcpl\n",
    "\n",
    "\n",
    "def neg_cox_log_likelihood(labels, risk):\n",
    "    return -1. * cox_log_likelihood(labels, risk)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "tender-extreme",
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "function to fit model with given hyperparameters\n",
    "'''\n",
    "def model_fitting(model, initial_learning_rate,decay_steps, decay_rate, epochs,train_dataset,validation_dataset,fold,weights):# Compile model.\n",
    "# if cudnn error occurs in model compile, detach and attach environment again \n",
    "    initial_learning_rate = initial_learning_rate\n",
    "    lr_schedule = keras.optimizers.schedules.ExponentialDecay(\n",
    "        initial_learning_rate, decay_steps=decay_steps, decay_rate=decay_rate, staircase=True\n",
    "    )\n",
    "    model.compile(\n",
    "        loss=neg_cox_log_likelihood,\n",
    "        optimizer=keras.optimizers.Adam(learning_rate=lr_schedule,amsgrad=False),\n",
    "    )\n",
    "\n",
    "    # Define callbacks.\n",
    "    #https://www.pyimagesearch.com/2021/06/30/how-to-use-the-modelcheckpoint-callback-with-keras-and-tensorflow/\n",
    "    # save the trained model at the end of each fold\n",
    "    fname = \"3d_image_cox_reg_\"+ str(fold) + \".hdf5\"\n",
    "    fname=os.path.join(weights,fname)\n",
    "    checkpoint_cb = keras.callbacks.ModelCheckpoint(fname, save_best_only=True, mode=\"min\",verbose=1)\n",
    "    #checkpoint_cb = keras.callbacks.ModelCheckpoint(fname, mode=\"min\",verbose=1, period=epochs) # if want to save model at the end of all epochs\n",
    "    early_stopping_cb = keras.callbacks.EarlyStopping(monitor=\"val_loss\", patience=2)\n",
    "    \n",
    "    # Train the model, doing validation at the end of each epoch\n",
    "    model.fit(\n",
    "        train_dataset,\n",
    "        validation_data=validation_dataset,\n",
    "        epochs=epochs,\n",
    "        shuffle=True,\n",
    "        verbose=2,\n",
    "        #callbacks=[checkpoint_cb, early_stopping_cb, tensorboard_callback], #use tensorboard when performance visualization is required \n",
    "        callbacks=[checkpoint_cb, early_stopping_cb], \n",
    "    )\n",
    "    return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "regulation-people",
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Once the model is trained it is applied to train and vlaidation data to collect model predictions using evaluate function \n",
    "'''\n",
    "def evaluate(model, x_train, x_valid, fold, weights):\n",
    "    fname = \"3d_image_cox_reg_\"+ str(fold) + \".hdf5\"\n",
    "    fname=os.path.join(weights,fname)\n",
    "    model.load_weights(fname)\n",
    "\n",
    "    prediction_train = []\n",
    "    prediction_valid = []\n",
    "    for i in range(0,len(x_train)):\n",
    "        prediction = model.predict(np.expand_dims(np.expand_dims(x_train[i], axis=3), axis=0))[0]\n",
    "        prediction_train.append(prediction[0])\n",
    "\n",
    "\n",
    "    for i in range(0,len(x_valid)):\n",
    "        prediction = model.predict(np.expand_dims(np.expand_dims(x_valid[i], axis=3), axis=0))[0]\n",
    "        prediction_valid.append(prediction[0])\n",
    "    return prediction_train, prediction_valid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "formed-celebration",
   "metadata": {},
   "outputs": [],
   "source": [
    "'''Compute concordance index'''\n",
    "def compute_ci(y_train, yhat_train, y_valid, yhat_valid):\n",
    "    CI_train  = concordance_index(y_train[:,1], yhat_train,y_train[:,0])\n",
    "    CI_valid  = concordance_index(y_valid[:,1], yhat_valid,y_valid[:,0])\n",
    "    \n",
    "    return CI_train, CI_valid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "unexpected-shareware",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "########### Model fitting for fold 1 of 4 ###########\n",
      "Train for 3 steps, validate for 3 steps\n",
      "Epoch 1/2\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\MBRO\\.conda\\envs\\deep\\lib\\site-packages\\tensorflow_core\\python\\framework\\indexed_slices.py:433: UserWarning: Converting sparse IndexedSlices to a dense Tensor of unknown shape. This may consume a large amount of memory.\n",
      "  \"Converting sparse IndexedSlices to a dense Tensor of unknown shape. \"\n",
      "C:\\Users\\MBRO\\.conda\\envs\\deep\\lib\\site-packages\\tensorflow_core\\python\\framework\\indexed_slices.py:433: UserWarning: Converting sparse IndexedSlices to a dense Tensor of unknown shape. This may consume a large amount of memory.\n",
      "  \"Converting sparse IndexedSlices to a dense Tensor of unknown shape. \"\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Epoch 00001: val_loss improved from inf to 344.42479, saving model to D:\\Transcend_(E)\\_IramS\\PETra_ML_DL_Analysis\\Deep_learning\\Experiments_Cox_PET\\Test\\Weights\\3d_image_cox_reg_1.hdf5\n",
      "3/3 - 24s - loss: 368.1460 - val_loss: 344.4248\n",
      "Epoch 2/2\n",
      "\n",
      "Epoch 00002: val_loss improved from 344.42479 to 335.61105, saving model to D:\\Transcend_(E)\\_IramS\\PETra_ML_DL_Analysis\\Deep_learning\\Experiments_Cox_PET\\Test\\Weights\\3d_image_cox_reg_1.hdf5\n",
      "3/3 - 1s - loss: 371.6923 - val_loss: 335.6111\n",
      "########## CI Train at fold 1 is 0.44 ############\n",
      "########## CI Valid at fold 1 is 0.57 ############\n",
      "########### Model fitting for fold 2 of 4 ###########\n",
      "Train for 3 steps, validate for 3 steps\n",
      "Epoch 1/2\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\MBRO\\.conda\\envs\\deep\\lib\\site-packages\\tensorflow_core\\python\\framework\\indexed_slices.py:433: UserWarning: Converting sparse IndexedSlices to a dense Tensor of unknown shape. This may consume a large amount of memory.\n",
      "  \"Converting sparse IndexedSlices to a dense Tensor of unknown shape. \"\n",
      "C:\\Users\\MBRO\\.conda\\envs\\deep\\lib\\site-packages\\tensorflow_core\\python\\framework\\indexed_slices.py:433: UserWarning: Converting sparse IndexedSlices to a dense Tensor of unknown shape. This may consume a large amount of memory.\n",
      "  \"Converting sparse IndexedSlices to a dense Tensor of unknown shape. \"\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Epoch 00001: val_loss improved from inf to 337.77722, saving model to D:\\Transcend_(E)\\_IramS\\PETra_ML_DL_Analysis\\Deep_learning\\Experiments_Cox_PET\\Test\\Weights\\3d_image_cox_reg_2.hdf5\n",
      "3/3 - 21s - loss: 369.1481 - val_loss: 337.7772\n",
      "Epoch 2/2\n",
      "\n",
      "Epoch 00002: val_loss improved from 337.77722 to 327.40096, saving model to D:\\Transcend_(E)\\_IramS\\PETra_ML_DL_Analysis\\Deep_learning\\Experiments_Cox_PET\\Test\\Weights\\3d_image_cox_reg_2.hdf5\n",
      "3/3 - 1s - loss: 345.3280 - val_loss: 327.4010\n",
      "########## CI Train at fold 2 is 0.45 ############\n",
      "########## CI Valid at fold 2 is 0.43 ############\n",
      "########### Model fitting for fold 3 of 4 ###########\n",
      "Train for 3 steps, validate for 3 steps\n",
      "Epoch 1/2\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\MBRO\\.conda\\envs\\deep\\lib\\site-packages\\tensorflow_core\\python\\framework\\indexed_slices.py:433: UserWarning: Converting sparse IndexedSlices to a dense Tensor of unknown shape. This may consume a large amount of memory.\n",
      "  \"Converting sparse IndexedSlices to a dense Tensor of unknown shape. \"\n",
      "C:\\Users\\MBRO\\.conda\\envs\\deep\\lib\\site-packages\\tensorflow_core\\python\\framework\\indexed_slices.py:433: UserWarning: Converting sparse IndexedSlices to a dense Tensor of unknown shape. This may consume a large amount of memory.\n",
      "  \"Converting sparse IndexedSlices to a dense Tensor of unknown shape. \"\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Epoch 00001: val_loss improved from inf to 335.25896, saving model to D:\\Transcend_(E)\\_IramS\\PETra_ML_DL_Analysis\\Deep_learning\\Experiments_Cox_PET\\Test\\Weights\\3d_image_cox_reg_3.hdf5\n",
      "3/3 - 20s - loss: 338.3135 - val_loss: 335.2590\n",
      "Epoch 2/2\n",
      "\n",
      "Epoch 00002: val_loss did not improve from 335.25896\n",
      "3/3 - 1s - loss: 337.1817 - val_loss: 335.4157\n",
      "########## CI Train at fold 3 is 0.42 ############\n",
      "########## CI Valid at fold 3 is 0.40 ############\n",
      "########### Model fitting for fold 4 of 4 ###########\n",
      "Train for 3 steps, validate for 3 steps\n",
      "Epoch 1/2\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\MBRO\\.conda\\envs\\deep\\lib\\site-packages\\tensorflow_core\\python\\framework\\indexed_slices.py:433: UserWarning: Converting sparse IndexedSlices to a dense Tensor of unknown shape. This may consume a large amount of memory.\n",
      "  \"Converting sparse IndexedSlices to a dense Tensor of unknown shape. \"\n",
      "C:\\Users\\MBRO\\.conda\\envs\\deep\\lib\\site-packages\\tensorflow_core\\python\\framework\\indexed_slices.py:433: UserWarning: Converting sparse IndexedSlices to a dense Tensor of unknown shape. This may consume a large amount of memory.\n",
      "  \"Converting sparse IndexedSlices to a dense Tensor of unknown shape. \"\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Epoch 00001: val_loss improved from inf to 332.58365, saving model to D:\\Transcend_(E)\\_IramS\\PETra_ML_DL_Analysis\\Deep_learning\\Experiments_Cox_PET\\Test\\Weights\\3d_image_cox_reg_4.hdf5\n",
      "3/3 - 21s - loss: 335.5615 - val_loss: 332.5837\n",
      "Epoch 2/2\n",
      "\n",
      "Epoch 00002: val_loss improved from 332.58365 to 327.38768, saving model to D:\\Transcend_(E)\\_IramS\\PETra_ML_DL_Analysis\\Deep_learning\\Experiments_Cox_PET\\Test\\Weights\\3d_image_cox_reg_4.hdf5\n",
      "3/3 - 1s - loss: 334.3924 - val_loss: 327.3877\n",
      "########## CI Train at fold 4 is 0.41 ############\n",
      "########## CI Valid at fold 4 is 0.41 ############\n"
     ]
    }
   ],
   "source": [
    "i=1\n",
    "splits=2 # number of splits for Cross-validation\n",
    "repeats=2 # number of repeats for Cross-validation\n",
    "kfold  = RepeatedStratifiedKFold(n_splits=splits,n_repeats=repeats,random_state=0) # create k-folds\n",
    "CI_train_all=[] # create empty list to store all concordance index train\n",
    "CI_valid_all=[] # create empty list to store all concordance index valid\n",
    "\n",
    "# here we will store all \n",
    "risk_ensemble_train =pd.DataFrame()\n",
    "risk_ensemble_train['subject_id'] = X.subject_id # store all subject ids, \n",
    "risk_ensemble_valid =pd.DataFrame()\n",
    "risk_ensemble_valid['subject_id'] = X.subject_id\n",
    "\n",
    "for train, valid in kfold.split(X, y): #train and valid contain index position of splits created \n",
    "    X_train, X_valid = X.iloc[train.tolist()], X.iloc[valid.tolist()] #get train and valid IDs  \n",
    "    X_train_subvolumes = data.loc[data.subject_id.isin(np.array(X_train.subject_id))] \n",
    "    X_valid_subvolumes =  data.loc[data.subject_id.isin(np.array(X_valid.subject_id))]\n",
    "    x_train = np.array([process_scan(path) for path in X_train_subvolumes.Path])\n",
    "    y_train = np.array(X_train_subvolumes[['status','time']])\n",
    "    x_val = np.array([process_scan(path) for path in X_valid_subvolumes.Path])\n",
    "    y_val = np.array(X_valid_subvolumes[['status','time']])\n",
    "    \n",
    "    train_dataset,validation_dataset = data_loaders(x_train,y_train,x_val,y_val)\n",
    "    \n",
    "    ''' Three models can be used here\n",
    "    1- DenseNet3D imported from DenseNet3D.py\n",
    "    2- ResNet3D imported from ResNet3D.py\n",
    "    3- Vgg3D imported from Vgg3D.py file\n",
    "    comment models that are not required\n",
    "    '''\n",
    "    \n",
    "    #model =Vgg3D.Vgg3D_3blocks(input_shape=(60,60,44,1))\n",
    "    #model =ResNet3D.Resnet3DBuilder.resnet_v3((60,60,44,1),1,reg_factor=1e-4)\n",
    "    model = DenseNet3D.DenseNet3DImageNet121(input_shape=(60,60,44,1), include_top=False)\n",
    "    \n",
    "    model = Model(inputs=model.input, outputs=model.output)\n",
    "    print(\"########### Model fitting for fold %d of %d ###########\"% (i,(splits*repeats)))\n",
    "    \n",
    "    #Model hyperparameters\n",
    "    trained_model= model_fitting(model,initial_learning_rate= 0.0001,decay_steps=1000,decay_rate=0.96,epochs=2,\n",
    "                                 train_dataset=train_dataset,validation_dataset=validation_dataset,fold=i,weights=weights)\n",
    "    \n",
    "    yhat_train, yhat_valid = evaluate(trained_model, x_train, x_val,i,weights)\n",
    "    \n",
    "    \n",
    "    \n",
    "    ##################### collect Subject_risk in current fold for computing Ensemble risk ##############################\n",
    "    \n",
    "    '''Train'''\n",
    "    for s_id_train, yhat_t in zip(np.array(X_train_subvolumes.subject_id), np.array(yhat_train)):\n",
    "        # X_train subvolumes have subject ids selected in current fold, These id are matched with risk_ensemble_train dataframe \n",
    "        # where id is matched, the prediction is placed in front of it for each cv fold \n",
    "            index = risk_ensemble_train.index[risk_ensemble_train['subject_id'] == s_id_train].tolist() \n",
    "            risk_ensemble_train.loc[index[0], 'risk_fold_'+str(i)] = yhat_t\n",
    "    '''Valid'''\n",
    "    for s_id_valid, yhat_v in zip(np.array(X_valid_subvolumes.subject_id), np.array(yhat_valid)):\n",
    "            index = risk_ensemble_valid.index[risk_ensemble_valid['subject_id'] == s_id_valid].tolist()\n",
    "            risk_ensemble_valid.loc[index[0], 'risk_fold_'+str(i)] = yhat_v\n",
    "            \n",
    "    ci_train_all, ci_valid_all = compute_ci(y_train,yhat_train,y_val,yhat_valid)\n",
    "    CI_train_all.append(ci_train_all)\n",
    "    CI_valid_all.append(ci_valid_all)\n",
    "    print('########## CI Train at fold %d is %.2f ############' % (i,(ci_train_all)))\n",
    "    print('########## CI Valid at fold %d is %.2f ############' % (i,(ci_valid_all)))\n",
    "    i=i+1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "furnished-benefit",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Calculate the ensemble risk for train and validation \n",
    "risk_ensemble_train.fillna(0)\n",
    "risk_ensemble_train['Ensemble'] = risk_ensemble_train.loc[:,'risk_fold_1': 'risk_fold_'+str(splits*repeats)].mean(axis=1)\n",
    "risk_ensemble_train.fillna(0)\n",
    "risk_ensemble_valid['Ensemble'] = risk_ensemble_valid.loc[:,'risk_fold_1': 'risk_fold_'+str(splits*repeats)].mean(axis=1)\n",
    "#Calculate CI from ensemble risk \n",
    "risk_ensemble_train['status'],risk_ensemble_train['time'] = data['status'],data['time']\n",
    "risk_ensemble_valid['status'],risk_ensemble_valid['time'] = data['status'],data['time']\n",
    "CI_ensemble_train  = concordance_index(risk_ensemble_train.time,risk_ensemble_train.Ensemble,risk_ensemble_train.status)\n",
    "CI_ensemble_valid  = concordance_index(risk_ensemble_valid.time,risk_ensemble_valid.Ensemble,risk_ensemble_valid.status)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "clear-dynamics",
   "metadata": {},
   "outputs": [],
   "source": [
    "def LogRank_test(risk_ensemble):\n",
    "    pred= logrank_test(np.array(risk_ensemble.loc[risk_ensemble.risk_group==0].time),\n",
    "                       np.array(risk_ensemble.loc[risk_ensemble.risk_group==1].time),\n",
    "                       np.array(risk_ensemble.loc[risk_ensemble.risk_group==0].status),\n",
    "                       np.array(risk_ensemble.loc[risk_ensemble.risk_group==1].status),\n",
    "                       alpha=0.95)\n",
    "    p_value = pred.p_value\n",
    "    return p_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "pressed-orientation",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get optimal cutoff and p-value for train data using log-rank test, optimal cutoff is one where p-value is minimum  \n",
    "p_value = list()\n",
    "cut_off = list()\n",
    "\n",
    "for i in range(0,risk_ensemble_train.shape[0]):\n",
    "    cut = risk_ensemble_train.Ensemble[i]\n",
    "    risk_ensemble_train['risk_group'] = np.where(risk_ensemble_train['Ensemble'] > cut, 1,0)\n",
    "    p_val= LogRank_test(risk_ensemble_train)\n",
    "    p_value.append(p_val)\n",
    "    cut_off.append(cut)\n",
    "    \n",
    "optimal_cutoff= cut_off[p_value.index(min(p_value))]\n",
    "optimal_p_value = min(p_value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "manual-prisoner",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Valid p-values at optimal threshold using log-rank test \n",
    "risk_ensemble_train['risk_group'] = np.where(risk_ensemble_train['Ensemble'] > optimal_cutoff, 1,0)\n",
    "risk_ensemble_valid['risk_group'] = np.where(risk_ensemble_valid['Ensemble'] > optimal_cutoff, 1,0)\n",
    "p_val_valid= LogRank_test(risk_ensemble_valid)\n",
    "\n",
    "# Final sheet with train and validation predictions for 25 CV folds\n",
    "risk_ensemble_train.to_excel('Ensemble_train_risk_at_optimal_cutoff.xlsx', index=None, header=True)\n",
    "risk_ensemble_valid.to_excel('Ensemble_valid_risk_at_optimal_cutoff.xlsx', index=None, header=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "successful-hobby",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA+gAAAHvCAYAAADHOfV1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABkeElEQVR4nO3deZgcVbn48e8LgZCwBkhYghAgyBLE/CQoCJJFNhVBBREiCnohuKFcUUFlGQh6Fa9cXFCJIiAQEBAFURFDFgUBAY1IQDBA2AQCJmwBAiHn90dVh06ne6ZnMtNVM/P9PE8/naml++0zk37rrTrnVKSUkCRJkiRJxVql6AAkSZIkSZIFuiRJkiRJpWCBLkmSJElSCVigS5IkSZJUAhbokiRJkiSVgAW6JEmSJEklYIEuqdtERGriMW8l3+PI/HVGdGHfC1b2/SVJqqjKSSO7sO/7IuLzPRFXd4uImRExs53178/b4ZB2trkoIl6IiDWbfM+2iEg1y1JEtK1svO3sNy5/31Vqlo/I3/vIzr6m1FkDig5AUp+yW83PvwT+DrRVLVu8ku/xm/x9Hu/CvpOB76zk+0uS1B3eB+wFnFVwHN3hWuA/wEeAy2tXRsRawPuBX6SUFq3E++wGPLoS+3dkHHAqcAawtGr54/l739+D7y0BFuiSulFK6ZbqnyNiMfB07fKabVYFIqW0pMn3eAp4qovxmVglSX1WRAxMKa3sifBOSym9GhGXAp+IiKF5rq52ELAmcOFKvk/D44melLdpIe+t/scu7pJaKu8i9rWIODEiHgReAd4UEWtExP9FxF15F7gnIuLXEbFdzf4rdHGPiHkRcXFEHBoR90TEooi4PSL2qNl3uS7uVV3WjomI0yPi8Yh4Jn/fzWr2HRwRP4yI/0TE8xHxy4h4u13eJEnV8u7VN0bEXhHx14h4Mc9t76va5gLgCGB4vSFgEbFhnnMei4jFEfHPiJhU8z6VfLhnRFwREc8At+brms2Lu0TElRHxaES8FBH3RsTXI2JQFz76hWQX/w6rs+6jwCPAzIgYGhHnRsR9eds8EhFTI2J4R29Qr4t7/hn/mbfTnIh4f539OjzGyF/31PzHVyu/l3xd3S7uEXF4RPw9Il6OiKfzbvyb1GzT1O9CqvAKuqQiHAk8AHwBWAT8GxgIrE3WrexxYH3gU8AtEbFdSumJDl7zHcC2wMnAy2Td2a+NiBEppWc62PfLwJ+BjwPDgG8DlwBjq7aZAnyQrLv+7cA7820kSaq1NdmQqv8BngaOB67M89lcshw1FNgFOCDfZzFARKwD3AQMIss5DwL7Aj/Mr5B/r+a9LgEuBQ5m+WP7ZvLi5sBs4ALgeWAUcAqwFXBoZz5wSun2iJhD1s39u5Xl+QnvccA3UkpLI2L9PJ4vk/WI2zRvn5vy9nm52feMiL2AqWTD344na9PvAKsB91Zt2swxxk+AzYD/AvYAXuvgvScB5wI/zz/LpsDXgbdFxFtSSi9Ubb4yxyjqZyzQJRUhgH1SSi/VLD9q2QZZ1/ffA0+SnY3/vw5ecx1gdEppYb7/E8BtwLvJknd7HkopTax676HAtyJi05TSvyNiW2AicGJK6cx8sz9ExGDg2A5eW5LU/2wI7JlS+hdARPyVrDA8BPh6Sun+iHgKeKVOt+3PAVsAb6rsD0yLiPWAUyPihzXDwq5MKX2pTgwd5sWU0i8qG0dEkJ0YeA74WUR8OqX0n05+7p8B34yI7VNK9+TLDifrtfuz/D3vzT9j5X1Xzd/3YeBdZPPXNOs04J/AgSmlpfnr3UPWHX1ZgZ5SepYOjjFSSo9GRGV8+63tDb3L958MzEwpHVq1/J/An8hO+H+3apeVOUZRP2MXd0lFuK5OcU5EHBIRt+bd9JaQXV1fi+ysc0duriS+3D/y582b2Pc3NT/X7vs2spMKV9Rsd2UTry1J6n/+VVVck1KaD8ynuZy0H1lX9QcjYkDlQVZQbgDsULN9o4K2w7wYEetExDcj4n6yK/ivAheR5bxtmoi11sVkV54/UrXsI2QF77KCOSI+mXcNf4Es3z+cr2om31deY1WyHghXVopzgJTSrcC8OtuvzDFGrW3Jetwt15MupXQj8BDL98CDlTtGUT9jgS6pCCvMwB4R7yXrJnYP2dXqt5El3qeANZp4zQXVP1RNktPpfXl9pvnKvpXxZPNrtnuyideWJPU/tXkFstzSTE4aBuxJVixXPyoniTeo2b7RXU2ayYvnA58gu9q7N1ne/XSd7ZqSUvo3MA04PDJjyE4oLJscLiKOBX6Qb/cB4K3Arl14zw3JurLXy8XLLeuGY4xa6+fP9dr+iar1FStzjKJ+xi7ukoqQ6iw7FJibUjqysiAiVmPFJFeESgIeRjYWsGKjAmKRJPVt/yE7Ify5Buvvrfm5Xk7tUESsARwItKWUvlO1/E1deb0qF5J12x5Hdiu5xWTFccWhwA0ppeOr3nPLLrzP02QnLurl4o3IrmRXv2d3HmNUCu6N66zbmGyuGqlLvIIuqSwGk3U5q/YRYNUCYql1K9kB0Adrltf+LElSsxaTTQRX6zpgO+DhlNLtdR7Pd9P7DyTLsa/WLD9yJV/3V8CzZOOwDwOuTSlVX0EeXOc9P9bZN0kpvUY2jvvgiFhW00TE24ARNZs3e4xRubLd0Sz295JdpV9uIr2IeDvZ/AGzOthfasgr6JLK4jrgfRHxf8C1wM7AZ4FnigwKsgltImIqMDk/CLgDmAC8N99kacOdJUmq725g/Yj4JNkV15dTSv8gmxT1Q8Cf8px4L9k9xLcD3pFSOrA73jyl9GxE3AIcHxGPk12R/jjQ4e3OOnjdlyLiCrLZ0IMV731+HXBCRHwF+AtZPj24i293KnA98KuIOJdsFvfTyLqZ175nM8cYd+fPx0fE74DXUkorXA1PKb0WEacA50bExWRj74cDXwP+RTZ0QOoSC3RJZfFj4A1kBwfHkJ0Vfy+dm821J00iuwXNl4DVgelk4/SuJbtSIElSZ/yEbOz114H1yLpkj8gL57eT3e7sBLLC7xmyQv0XdV+p6w4DfgicA7wEXE7Wtf7alXzdC8lmTX8K+F3NutPJPu9/k43BnkV2G7kHOvsmKaVpEfFhstvRXQXMBY5jxeEBzR5jXEs2Pv5TZO0f+aPee0+JiBeBLwJXAy8AvwW+VHOLNalTIqUuDVtRJ+T3fzwBGAO8mazbzJYppXlN7LsG2W0cDif7MpsNnJBS+mMPhSupSRHxReCbZAdUD3e0vaRyMC9LksrKMeitMZLsvpcLye6N2BnnAUeTncXbn2yyqt9HxOjuDFBS+yJi/4g4MSL2i4h9ImIy2VWAyy3OpV7HvCxJKiWvoLdARKxSuT9jRBxF1s2mwzP1EfFmsjPzH08pnZ8vGwDMAe5NKR3Qk3FLel1EjCW7Wr4d2VjAx8hmpT01pfRykbFJ6hzzsiSprByD3gKVg4AuOIBslstlt6ZIKS2JiMuAEyNiYNV9FCX1oJTSLF6/T6ukXsy8LEkqK7u4l9so4MGU0os1y+eQTVI1svUhSZLUb5mXJUk9ygK93NYnGx9Xa0HVekmS1BrmZUlSj7KLe7kFUG+SgLq3ewCIiElkt4NizTXX3Hm77bbrodAkSUW74447nk4pDS06jn7EvCxJaqg78rIFerktADavs3xI1frlpJSmAFMAxowZk26//faei06SVKiIeKjoGPoZ87IkqaHuyMt2cS+3OcCWETG4ZvkOwCvA3NaHJElSv2VeliT1KAv0crsGWA34YGVBfjuXDwHXO1OsJEktZV6WJPUou7i3SEQcnP9z5/z5XRHxFPBUSmlWRGwB3A+cnlI6HSClNDsifg6cHRGrAQ8CnwS2BD7c2k8gSVLfYV6WJJWRBXrrXFHz8w/y51nAOLIJZlZlxV4NHwO+BpwBrAf8HdgvpfTXngpUkqR+wLwsSSodC/QWSSk1nOE1Xz+POrPAppReAj6fPyRJUjcwL0uSysgx6JIkSZIklYBX0CWp5J577jnmz5/Pq6++WnQoapHVVluNYcOGsc466xQdiiSphnm5/2llXrZAl6QSe+6553jyyScZPnw4gwYNIqLdXrnqA1JKvPTSSzz22GMAFumSVCLm5f6n1XnZLu6SVGLz589n+PDhDB482IOAfiIiGDx4MMOHD2f+/PlFhyNJqmJe7n9anZct0CWpxF599VUGDRpUdBgqwKBBg+w+KUklY17uv1qVly3Q+7Ln7oW5U4qOQtJK8gx9/+TvvQ8yL0t9gt/P/VOrfu8W6H3Zkhdg3tSio5AkSWBeliR1yAK9LxuwVtERSJKkCvOyJKkDFuiSpJa54IILiAjmzp27wrolS5YQEbS1ta2w/bx58zr9XkceeSSbbbbZSkRbX1tbW6e7uY0bN4499tij22ORJGllmJfLx9us9XULZ8O0cc1tO2IijJwEd7bBTm09FpIkNes973kPN998M5tssknRoSxz1FFHsd9++xUdhvqTO9vMy5JKwbzc8yzQ+7KB68OQLZvbduHs7HnkJLjrNA8EJJXC0KFDGTp0aNFhALB48WIGDhzIZptt1iNXANRPdOXEuXlZUkmYl3ueXdz7soFDYa+ZzT2GjC4oSElqrF5XuhdffJFPfvKTbLDBBqy99tq8//3v589//jMRwQUXXLDCa/ztb3/jHe94B4MHD2abbbbhRz/6UYfvO3PmTCKCq666iqOPPpqhQ4ey0UYbAfW70n3nO99h++23Z9CgQQwZMoQxY8bwy1/+st33mDx5MquvvjqXXHJJxw2hvmHg+s3n24WznVBOUumYl3ueV9AlSS332muvsWTJkhWWNWPSpElcccUVtLW1MWbMGG644QY+/OEP1932ueeeY+LEiRx33HGccsopnH/++Xzyk59k2223Zfz48R2+17HHHsu73vUuLrroIl5++eW621xyySUcf/zxnHLKKbzjHe/gpZde4s4772TBggV1t1+6dCmf+cxn+NnPfsavf/1r9t1336Y+t/qAyonzZjR7lV2SuoF5uTx52QJdknqZ446D2bOLjWH0aDj77K7vv91223Vpv3vvvZepU6fyjW98gy996UsA7L333rz44ot873vfW2H7559/nh/84AfLkv6ee+7J9ddfz6WXXtrUgcBb3/pWfvKTn7S7zc0338xOO+3EKaecsmzZu9/97rrbLl68mIkTJ/LHP/6R6dOn89a3vrXDGCRJ5WZeNi93J7u463XV4+LmTikyEkl93C9/+Utuu+225R633HJLh/vdeuutpJT44Ac/uNzygw8+uO72gwcPXi7hDxw4kG222YaHH364qTjf//73d7jNLrvswuzZszn22GOZNm0aL774Yt3tnn/+efbZZx/uuOMObrzxxlIcBKjkzMuSWsS8XJ687BX0/uzOtmzimWrzZ2XPfzkme1Tb8VQnqZFKYGXOkJfFjjvuyMiRI5dbVtu1rp7HH38cgGHDhi23vDIOrdaQIUNWWDZw4MCG3eJqNTNL7Uc/+lFefvllzjvvPH7wgx+w2mqr8e53v5uzzjqLESNGLNvu4YcfZs6cORx99NFsu+22Tb2/+pk728zLUi9kXjYvdycL9P5sp7b6iX1qwLCxzY+Tk6QWqSTm+fPns+WWr9+l4sknn+yR92vmvqoRwTHHHMMxxxzDwoULuf766zn++OP50Ic+xK233rpsu1GjRvHpT3+aj3zkIwwaNIizzjqrR2JWL2ZeltTLmJe7n13cJUm9xtve9jYigiuuuGK55bU/F2XIkCF86EMf4pBDDuGuu+5aYf1hhx3GpZdeyve+9z2OO+641gcoSVI3Mi93P6+gq775s7LxbiMnFR2JJC2z7bbbMnHiRE4++WSWLl3KzjvvzPTp0/n1r38NwCqrtP6886RJk1h77bXZbbfdGDZsGPfddx8XXXQR++yzT93tP/jBD7LKKqtw2GGHsXTpUr773e+2OGJJkrqHebn7WaCrsXlTLdAllc6UKVNYe+21OfPMM3nllVeYMGEC55xzDvvvvz/rrrtuy+PZfffdOf/887nooot49tln2XTTTTn88MM57bTTGu5z0EEHcfnll/OhD32I1157je9///tNddtTP+aJc0klZV7uXpFSKuSN1fPGjBmTbr/99s7vWBnrBo53kwp2zz33sP322xcdRul961vf4oQTTmDevHlsvvnmRYfTbTr6/UfEHSmlMS0MSSthpfIyOA5dKgHzcnPMy13nFXStaMdTYf7MoqOQpLquvfZa7rrrLkaPHs0qq6zCn/70J/73f/+XQw45pE8dBEjLmJcllZh5uXtZoGtFO7W9ft9VSSqZtddem1/96ld84xvfYNGiRQwfPpzPfvaz7XZdk3o187KkEjMvdy8LdDXmeDdJJTR27FhuueWWosOQWm/h7KxQHzHR3CypNMzL3cvbrKm+EROz53lTi41DkiRleXnI6KxINzdLUp9lga76Rk56faI4SZJUrJGTsgnihowuOhJJUg+yQJckSZIkqQQs0NW+yjh0SZIkSVKPskBXY45DlySpfDx5Lkl9lgW6GnMcuiRJ5eLJc0nq0yzQJUmSegtPnktSn2aBLklqmQsuuICIYO7cuUWH0qF58+YREVxwwQVN79ObPp8kSb0pb/WXvDyg6ADUC1TGuo2cVHQkktQym2yyCTfffDNbb7110aFIktTv9Ze87BV0tc+xbpL6mZQSr7zyCgMHDmTXXXdl6NChRYckSVK/1d/ysgW62udYN0kt9uqrr3LSSScxYsQIVl99dUaMGMFJJ53Eq6++umybHXfckaOOOmrZz88++yyrrroqm2222XKvtfvuu3PIIYe0+34jRozg8MMP56c//Snbbbcdq6++Or/5zW/qdqW77bbb2Hvvvdlggw0YPHgwW221FZ/61Kfaff077riDjTbaiA984AO8/PLLnWgJSZKKZ15uLbu4S5JK5YgjjuDyyy/nK1/5CnvssQc333wzZ5xxBg888ABTp2a9eSZMmMC11167bJ+ZM2cycOBAHnvsMe677z7e+MY3smjRIm677TYOP/zwDt9zxowZzJ49m1NPPZVhw4YxYsSIFbZ54YUX2HfffXnrW9/KBRdcwNprr828efP485//3PB1r7/+eg466CA+/OEPc84557Dqqqt2vkEkSSqQebm1LNDVnIWzYdq4rMu7Y9GlYt1xXPZ/skhDRsPOZ3f7y951111ceumlnHrqqbS1tQGwzz77sOqqq3LyySdz4oknstNOOzF+/Hi+973v8dBDD7HFFlswY8YM9tprL+655x5mzJjBG9/4Rv70pz/x6quvMn78+A7fd+HChdxxxx1svPHGy5bNmzdvuW3++c9/snDhQs4880x22mmnZcuPPPLIuq95ySWX8LGPfYwTTzyR008/vdNtIUnqJczL5uVuZBd3dWzExOw//cLZjkWX1KP++Mc/Aqxwdr3y86xZswAYO3Ysq6yyCtOnTwdg+vTpTJgwgQkTJiy3bJNNNmG77bbr8H133XXX5Q4C6tlmm21Yb731OOaYY7j44ot55JFHGm579tlnc+SRR/Kd73yn1AcB6sUqJ87nTik6Ekl9mHm59byCro6NnJQ9po0rOhJJ0CNnyMtiwYIFQDZTa7VKkq6sX3/99Xnzm9/MjBkzeO9738tdd93F+PHj2Xjjjfnc5z4HZN3jmjlLX+/96ll33XWZMWMGkydP5lOf+hTPP/88o0aN4rTTTuOggw5abtvLLruM4cOHr7Bc6haVCVwrV+zs2SYVy7xsXu5GXkGXJJXG+uuvD8ATTzyx3PLKzxtssMGyZePHj2f69OnMmDGDDTbYgJ122okJEyYwf/58brrpJv72t781fSAQEU1tN3r0aH7xi1+wYMGCZbd6OeSQQ7jrrruW2+4Xv/gFAwcOZNy4cSt8FmmljZwEe83MerdJUg8yL7eeBbokqTTGjs3uGnHZZZctt/ySSy4BYM8991y2bPz48Tz22GOce+65jBs3johg2LBhjBo1ilNPPZXXXnuNCRMm9EicAwYMYNddd2Xy5MksXbqUe+65Z7n1w4cPZ+bMmSxdupTx48fz+OOP90gckiT1JPNy69nFXZ3jZHGSusF11123wtiyddddl7333pvDDjuMtrY2lixZwtvf/nZuvvlmJk+ezGGHHbbcJDB77rknq666KjfccAPnnHPOsuXjx4/n+9//PptvvjlbbbVVt8V87bXXMmXKFN73vvex5ZZbsmjRIr773e+y9tprs9tuu62w/SabbMLMmTN55zvfybhx45gxYwabbrppt8UjAeZlSd3CvFweFuhqnmPeJHWTY489doVlo0aN4q677uLCCy9kq6224qc//SlnnHEGm266KSeccAKnnnrqctuvs8467LzzzvzlL39Z7oz8hAkT+P73v990N7pmbbPNNgwaNIjJkyfz+OOPs/baa7PLLrvwhz/8YYX7vFZsvPHGzJgxY7mDgeHDh3drXOrHzMuSuol5uTwipVR0DOohY8aMSbfffnv3v3Blsri9Znb/a0tazj333MP2229fdBgqSEe//4i4I6U0poUhaSWYl6Xez7zcv7UiLzsGXZIkSZKkErBAV9fMn+W9VyVJkiSpG1mgq/MqY97mTS02DkmSlPHEuST1CRbo6ryRk2DY2KKjkCRJ4IlzSepDLNAlqeSczLN/8veupnniXGopv5/7p1b93i3QJanEBgwYwJIlS4oOQwVYsmQJAwZ4N1RJKhPzcv/VqrxsgS5JJbbGGmvwwgsvFB2GCvD888+zxhprFB2GJKmKebn/alVetkCXpBIbOnQoTz31FC+++KJd6vqJlBIvvvgiTz/9NEOHDi06HElSFfNy/9PqvGzfOUkqsTXWWIONNtqIJ554gsWLFxcdjlpk4MCBbLTRRl5Bl6SSMS/3T63MyxboklRy6667Luuuu27RYUiSJMzL6ll2cVfXec9VSZIkSeo2FujqGu+5KklSuXjiXJJ6PQt0dY33XJUkqTw8cS5JfYIFuiRJUm/niXNJ6hMs0CVJkiRJKgEL9BaIiDdExJUR8WxEPBcRV0XE5k3uu3lEXBgRD0fEixFxX0ScERFr9nTcTVk4G6aNc8ybJKnXMC9LksrK26z1sIgYDEwHFgNHAAk4A5gRETullBa1s++awDRgNeBk4GFgF+A0YBvgQz0bfQcq490Wzs6eR04qLBRJkpphXpYklZkFes87GtgK2DalNBcgIu4E/gUcA5zVzr67kyX8fVNK1+fLZkTE+sAXImJwSunFngu9AyMnZY9p4woLQZKkTjIvS5JKyy7uPe8A4JbKQQBASulB4CbgwA72XT1/fq5m+TNkv7vophglSeovzMuSpNKyQO95o4C76iyfA+zQwb7TyM7ofzMidoiItSJiAvA54EftdcNrOe+9KknqHfpHXpYk9UoW6D1vfWBhneULgCHt7ZhSehnYg+z3NAd4HrgBuBb4TPeGuRK896okqffo+3kZPHEuSb2UBXprpDrLOuwGFxFrAD8HhgEfAcYCXySbhOacBvtMiojbI+L2p556qusRd4b3XpUk9S59Oy974lySei0niet5C8nO1tcaQv0z+NX+CxgHjEwp3Z8v+2NEPAtMiYgfpZT+Xr1DSmkKMAVgzJgx9Q5AJEnqz/p+Xh45yeJcknopr6D3vDlk491q7QDc3cG+bwIWVh0EVPwlf95+JWOTJKm/MS9LkkrLAr3nXQPsGhFbVRZExAiyW7Vc08G+TwBDImJkzfK35c+PdVeQkiT1E+ZlSVJpWaD3vB8D84CrI+LAiDgAuBp4BDi3slFEbBERSyLilKp9LyCbgOa3EXFERIyPiC8C/wvcQXZLGEmS1DzzsiSptCzQe1h+y5UJwH3ARcAlwIPAhJTSC1WbBrAqVb+TlNI8YFdgNnAG8FvgaLKxbHunlJb2/CeQJKnvMC9LksrMSeJaIKX0MHBQB9vMo84Msimlu4FDeiYySZL6H/OyJKmsvIIuSZIkSVIJWKBLkiRJklQCFuiSJEmSJJWABbokSZIkSSVggS5JkiRJUglYoEuSJEmSVAIW6JIkSZIklYAFuiRJkiRJJWCBLkmSJElSCVigq/vMnwVzpxQdhSRJkiT1Shbo6h4jJmbP86YWG4ckScp44lySeh0LdHWPkZNg2Niio5AkSeCJc0nqpSzQJUmS+hpPnEtSr2SBLkmSJElSCVigS5IkSZJUAhbokiRJkiSVgAW6JEmSJEklYIEuSZIkSVIJWKCre3nPVUmSJEnqEgt0dR/vuSpJUrl44lySehULdHUf77kqSVJ5eOJcknodC3RJkqS+yBPnktTrWKBLkiRJklQCFuiSJEmSJJWABbokSZIkSSVggS5JkiRJUglYoEuSJEmSVAIW6JIkSZIklYAFuiRJkiRJJWCBLkmSJElSCVigS5IkSZJUAhbokiRJkiSVgAW6JEmSJEklYIEuSZLUl82fBXOnFB2FJKkJFuiSJEl91YiJ2fO8qcXGIUlqigV6OyJzQET8b0ScHxFb5MvHRsSmRccnSZLUrpGTYNjYoqOQJDXJAr2BiBgC/Bn4FXAU8FFgg3z10cCJxUTWC9iVTpIkSZI6zQK9sW8BbwB2BzYEomrdNOCdRQRVenalkyRJkqQusUBv7EDgqymlm4FUs+5hsuJdtexKJ0mSJEldYoHe2FrAYw3WrcHyV9QlSZIkSVopFuiN3Qvs02DdWOAfLYxFkiRJktTHDSg6gBI7BzgnIp4FKgOq14uIjwGfASYVFpkkSZIkqc+xQG8gpfTjiNgaOA04PV/8B2ApcGZK6ZLCgpMkSZIk9TkW6O1IKZ0YET8E9gaGAf8B/pBSeqDYyCRJkiRJfY0FegdSSg8BPyk6DkmSJElS3+YkcQ1ExMcioq3BuraIOKLFIfUu82fB3ClFRyFJkiRJvYYFemOfI+vSXs984LjWhdLLjJiYPc+b2v52kiSpNTxxLkm9ggV6YyOBOQ3W3QNs3cJYepeRk2DY2KKjkCRJ4IlzSepFLNAbWwJs2GDd0FYGIkmS1GWeOJekXsNJ4hr7C/AJ4PI66z4B3NbacHqhhbNh2rjll42YmB0oSJIkSZKWY4He2NeAaRFxK9ks7o8Bw4GjgLeQ3XpNjVS601VbODt7tkCXJEmSpBVYoDeQUpoVEQcDZwPnVq2aBxyUUppZQFi9x8hJKxbitVfTJUlS69Tr2Qb2bpOkErFAb0dK6Wrg6ojYFtgAeDqldF/BYUmSJHVOvZ5tYO82SSoZC/QmpJTuLToGSZKkLqvXsw3s3SZJJWOB3o6IWAd4N7A5sEbN6pRSmtz6qHq56u51dqmTJEmSpGUs0BuIiN2BXwPrNdgkARbonVHdvc4udZIklUPl5LknziWpcBbojZ1NNiHc0cA/UkqvFBpNX1Ddvc4udZIkFa9y8twT55JUCqsUHUCJbQ+clFK6w+JckiT1SSMnwV4zYcjooiORJOEV9PY8DAwsOog+rd7tXuxeJ0mSJKmf8gp6Y6cBJ+YTxam7jZi44tn6hbNh3tQiopEkSZKkwnkFvbH9gY2AByPiZmBBzfqUUjqimReKiDcA/wfsDQQwDTgupfRwk/tvD5wOjAfWJLu6/4OU0nea2b+U6t3uxXHpkqQWMC83YM82SSqcBXpje5DN1P4cMKrO+tTMi0TEYGA6sBg4It/vDGBGROyUUlrUwf5j8v1nAkcBzwLbAGs19SkkSdIy5uUGqu+0UuHEcZLUchboDaSUtuymlzoa2ArYNqU0FyAi7gT+BRwDnNVox4hYBbgQuCGl9P6qVTO6KTZJkvob83I99myTpFJwDHrPOwC4pXIQAJBSehC4CTiwg33HATvQzsGCJEnqFPOyJKm0LNDbERFrRsRnI+LKiJgREdvkyw+NiO2afJlRwF11ls8hS/Lt2SN/XiMibomIVyNifkR8NyIGNfn+kiTpdeZlSVJpWaA3kE8gcyfwLbKxZXsCa+erxwNfaPKl1gcW1lm+ABjSwb6b5s8/B64nm8zmTLIxb3WnO4+ISRFxe0Tc/tRTTzUZoiRJ/YZ5WZJUWo5Bb+zbZBPIbAP8G3ilat0soK0Tr1VvQrloYr/KCZSLU0qn5P+eGRGrAt+IiB1SSncv90YpTQGmAIwZM6apiewkSepnzMuSpFLyCnpjewOn5rdcqU2ojwHDm3ydhWRn62sNof4Z/Gr/yZ//ULP8+vx5dJMxSJKkjHlZklRaFuiNrQ4832DdusCrTb7OHOrfpm0H4O46y2v3hRVPEFTO8i9tMgZJkpQxL0uSSssCvbE7gYMarHsXcEeTr3MNsGtEbFVZEBEjgN3zde35HVk3+/1qlu+bP9/eZAySJCljXpYklZZj0Bv7FnBlRMDrE7/sEBEHAv9FdpuWZvwY+AxwdUScRHbWfTLwCHBuZaOI2AK4Hzg9pXQ6QErpPxHxP8DJEfEcMB0YA5wCXFh9ixhJktQU87IkqbS8gt5ASukq4FPAB4Fp+eKfAccBn0kpXdfk6ywCJgD3ARcBlwAPAhNSSi9UbRrAqqz4Ozkd+BJwCPBb4JNkJw+O7vSH6g3mz4K5U4qOQpLUR5mXJUll5hX0BiJiXeB8suS9GzCMbHKYP6eUGo1NryufaK5Rd/nKNvOoM4NsSikBZ+WPvm3ExKxAnzcVRk4qOhpJUh9lXpYklZVX0OuIiAFkxfg+KaVFKaVpKaWpKaXfd7Y4VyeMnATDxhYdhSRJqrBnmyS1lAV6HSmlJcCTwGtFxyJJklSIEROz53lT299OktRtLNAbuxg4quggJEmSCmHPNklqOcegNzYPmBgRtwFXA49Tc9/TlNJPC4hLkiRJktQHWaA3dk7+PBzYuc76BFigS5IkSZK6hQV6Y1sWHYAkSVLhKhPFeYcVSepxjkGvIyJWJ7vf+bCU0kONHgWHKUmS1LOcKE6SWsoCvY6U0ivAMcCgomPpl7yliyRJ5eBEcZLUUhbojf0NeFPRQfQ7nqmXJEmS1E9ZoDd2PPCFiNg/IqLoYPoNz9RLkiRJ6qecJK6xK4B1yW6xtiQi5rP8bdZSSmmLQiKTJEmSJPU5FuiN3UDNfc8lSZIkSeopFugNpJSOLDoGSZIkSVL/4Rh0SZIkSZJKwCvoDUTERzvaJqX0s1bEIkmSJEnq+yzQG7ugwfLqcekW6JIkSZKkbmGB3tiWdZZtAOwPTAQOb204kiRJkqS+zAK9gZTSQ3UWPwT8Nb8v+ufJCnVJkiRJklaak8R1zZ+A9xQdhCRJkiSp77BA75pdgReKDkKSJEmS1HfYxb2BiDilzuLVgR3Jrp5/v7URSZIkSZL6Mgv0xtrqLFtMNg79a8D/tDSa/mb+LJg7BUZOKjoSSZIkSWoJu7g3kFJapc5jUEppu5RSW0ppcdEx9lkj8rn35k0tNg5JkiRJaiELdJXPyEkwbGzRUUiSpIpKzzZJUo+yQG8gIj4WEW0N1rVFxBEtDkmSJKn17NkmSS1jgd7Y54D/NFg3HziudaFIkiQVxJ5tktQyFuiNjQTmNFh3D7B1C2ORJEmSJPVxFuiNLQE2bLBuaCsDkSRJkiT1fd5mrbG/AJ8ALq+z7hPAba0NR5Ik9Wb33gvjxnXva06cCJO8I6kk9RkW6I19DZgWEbcCPwEeA4YDRwFvAfYuMDZJktTPzZ6dPVugS1LfYYHeQEppVkQcDJwNnFu1ah5wUEppZgFhSZKkXmrbbWHmzO57ve6+Gi9JKp4FejtSSlcDV0fEtsAGwNMppfsKDkuSJEmS1AdZoDchpXRv0TFIkiRJkvo2Z3FvR0T8v4i4KiKejoglEfGWfPnXI2K/ouOTJEmSJPUdFugNRMQewM3AdsBUlm+rpWQzuasnLZwN08bB3ClFRyJJkiRJPc4CvbFvAL8HRgGfr1n3V7KZ3NVTRkyEIaOzIn3e1KKjkSRJkqQeZ4He2FuAH6aUEpBq1j0NDG19SP3IyEmw18ysSJckScWzZ5sk9TgniWvsZWBwg3WbAM+2MBZJkqTijJiYPS+cnT2P9ObrktQTvILe2I3AcRGxatWyypX0/wKmtz4kSZKkAtizTZJawivojZ0M3AT8HbiSrDg/IiLOAnYGdikwNkmSJElSH+MV9AZSSn8H9gSeBL4KBPCZfPVY740uSZIkSepOXkFvR0rpr8A7I2INYH3gmZTSiwARsVFK6clCA5QkSZIk9RleQW9CSunllNK/U0ovRsSmEfEd4IGi45IkSZIk9R1eQa8REbsAHwM2B+YCZ6eU5kXEBkAbcBSwGvDzwoKUJEmSJPU5FuhVIuK9wC+AVYGngH2BQyPiILKJ4jYCfgmcklKaU1igkiRJkqQ+xy7uyzsBmANsnlLamGzc+Z+AP5C11YSU0kEW55IkSZKk7maBvrydgK+llB4DSCk9D3wJWAM4IaU0s8DYJEmSJEl9mF3cl7cW8HDNssrPXjUvysLZMG1c5/cbMRFGTuruaCRJkiSpR1igryg1WL6kpVEoM2Ji1/ZbODt7tkCXJEmS1EtYoK9oSkQ8X2f5eRHxQtXPKaU0tlVB9VsjJ3WtyO7KFXdJkiRJKpAF+vL+SP0r6LNaHYgkSVIpOfRMknqMBXqVlNK4omOQJEkqLYeeSVKPskBX31U5w+8Ze0mSuodDzySpR1mgq2+qnOH3jL0kSZKkXsL7oKtvGjkJ9poJQ0YXHYkkSZIkNcUCXZIkST2vMvRs7pSiI5Gk0rKLuyRJknqWQ88kqSleQZckSVLPcuiZJDXFAl2SJEmSpBKwQK8SEUsj4rUmH0s68bpviIgrI+LZiHguIq6KiM27EN+XIyJFxI2d3VeSJGXMywWbP8tx6JLUgGPQl3c6kLrzBSNiMDAdWAwckb/+GcCMiNgppbSoydfZCvgqML8745MkqT8xLxdsxMSsQJ831XHoklSHBXqVlFJbD7zs0cBWwLYppbkAEXEn8C/gGOCsJl/nh8AlwLb4e+ucyqyx1UZM9MBAkvon83KRRk7KinNJUl12ce95BwC3VA4CAFJKDwI3AQc28wIRMRF4C/DlHomwLxsxccUJaRbO9uBAkvov87IkqbQ849uOiFgdeBfZ2fE1alanlNLkJl5mFHB1neVzgA82EcMQ4P+AL6WUFkREE2+pZUZOWvFKee3VdElSf2JeliSVlgV6AxGxKXAjMIJsfFolA1ePUW+mQF8fWFhn+QJgSBP7fwu4D7igiW0lSVL7zMtlUDv8zKFnkgTYxb093wKeAjYnK87fRjZm7WvA3Pzfzao38VyHp9wj4h3AR4FPppSamrwuIiZFxO0RcftTTz3ViRD7mcqBgbPISlJ/ZF4uUu3wM4eeSdIyXkFv7B3AF4B/5z8vTSnNA06JiFWB79LcWLWFZGfraw2h/hn8aucC5wGPRsR6+bIBwKr5zy+llBZX75BSmgJMARgzZky3zkjfZ4yYmD0vnJ09e8ZekvoT83LRaoefOfRMkpbxCnpjGwD/TiktBRaxfLe36cC4Jl9nDtl4t1o7AHd3sO/2wCfIDhgqj92BXfN/f7LJGFRt5CTYa+aKk8dJkvoD87IkqbS8gt7Yo8CG+b/vB/YBpuU/vxV4ucnXuQb434jYKqX0AEBEjCBL6Cd2sO/4OsvOBlYFjiXrai9JkppnXi6jytAzx6JL6ucs0BubAYwFfkXWpe2ciBgNvArsmy9rxo+BzwBXR8RJZOPeJgOPVL9GRGxBdiLg9JTS6QAppZm1LxYRzwAD6q1TC93ZBju1FRyEJKkLzMtl49AzSVrGLu6NnQT8ECCl9EPgc8BgYBPgTOD4Zl4kpbQImEA24+tFwCXAg8CElNILVZsG2Rl4fye9wV2nFR2BJKkLzMsl5NAzSVrGK+gNpJSeBp6u+vl7wPe6+FoPAwd1sM08mphBNqU0risxSJKkjHlZklRWnhVuICKGRsQbG6x7Y0RsWG+dJEmSWuzOtqIjkKRuYYHe2A9o3I39v/P1kiRJKppDzyT1EXZxb2wP4NMN1l0PfL+FsainVGaN7Ypm93NGWkmSJElNsEBvbAjwbIN1z5HdJ129WWXW2Pa8MA9efKj+uvmzVlw2eAtYa8TrPzsjrSRJkqQmWaA39ijwNuCGOuveBjze2nDU7UZO6nrhPDVgYup4u65enZckSZLU71igN3Yl8JWI+HtK6TeVhRHxHuBE8luwSZIkqZs49ExSP2eB3tjpwJ7ANRHxBPAYMBzYGLgFcDYSSZKk7uLQM0myQG8kpfRiRIwFPgLsTTbmfC7ZBHEXp5SWFBmfepHqqwGetZckqT6HnkmSBXp7UkqvAj/NH1LnVV8N8Ky9JEmSpHZYoHcgInYExgLrA08Df0op3VVsVCrcjqc2t1311QDP2kuSJElqhwV6AxExALgAOAyIqlUpIqYCR6aUXisiNpXATm1FRyBJkrqi3kR0DkGTVBKrFB1AiZ0KHAKcAmwJDMqfTwE+lD9LkiSptxgxEYaMXn7Zwtkwb2oR0UjSCryC3tjhwOSU0teqlj0EfC0iVgU+RlbES82rnLX3TL0kSd2nK0PPKhyCJqlEvILe2KbAzQ3W/TlfLzWvctbeM/WSJHUvh55J6iMs0Bv7N7B7g3Vvz9dLzRs5CfaauWLXOkmSVKz5s2DulKKjkCQL9HZcAnw1Ik6OiK0iYlBEbBkRXwa+ClxUcHySJElaWZVbotq7TVIJOAa9sTZgK+C0/N8VAVwKnN76kDrn3nth3Ljufc2JE2GSQ6dXXuVMvePQJUkq1shJFueSSsMr6A2klJaklCYCbwI+QzZr+2eAHVNKH04pLSk0wALMng1TzV8rzzP1kiRJkurwCnoHUkpzgDnVyyJiL+CslNJOxUTVnG23hZkzu+/1uvtqfL/lmXpJkiRJdVigd826wKiig5AkSVI3qdwKtZq3RZXUYhbokiRJ6t8qw8+qLZydPVugS2ohC3SpKJ6plySpHEZOWjH/1uZoSWoBC3SpCJ6plyRJklTDAr1KRGzV5KYb92gg6vs8Uy9JkiSphgX68uYCqYntosntJEmS1FvVDkdzKJqkHmaBvryPFR2A+jkPBCRJKofa4WgORZPUAhboVVJKFxYdg/oxDwQkSSqP2uFoDkWT1AIW6FJZeCAgSZIk9WurFB2AJEmSJEnyCrokSZLUnNq5YiqcM0ZSN7FAl8rMAwFJksqhdq6YCueMkdSNLNClsvJAQJKk8qidK6bCOWMkdSMLdKmsPBCQJEmS+hUniZMkSZJWRmVI2twpRUciqZfzCrokSZLUVZUhaQ5Bk9QNvIIu9UaeqVc/dOONN/L2t7+dQYMGsfHGG/P5z3+el156qal9H3nkEQ4++GDWXXdd1llnHT7wgQ/w8MMP93DEkvqFkZNgr5kwZHTRkUgt1dN5+YYbbuDwww9n6623ZtCgQWy99dZ88pOfZP78+ctt99BDD3HggQeyxRZbMGjQIDbccEPGjRvH7373u277rK1kgS71NiMmZgcBC2fDvKlFRyO1xJ133snee+/NsGHDuPbaaznjjDM4//zzOfLIIzvc98UXX2TChAn885//5MILL+Siiy7iX//6F+PHj2fRokU9H7wkSX1MK/Lyj370I/7zn/9w0kkncd111/HlL3+Za665hl133ZUXXnhh2XYvvPACG264IWeccQa//e1vOe+881hrrbV497vfzVVXXdUTH79H2cVd6m0qk8c5WZz6kVNPPZXNNtuMK664gtVWWw2A1VdfnSOOOIITTjiBt7zlLQ33/fGPf8wDDzzAvffey8iRIwHYaaed2GabbTj33HP5/Oc/35LPIElSX9GKvPyDH/yAoUOHLttv7NixvPGNb2Ts2LFcfvnlfPzjHwdg1KhRnHfeecu9x3ve8x623HJLzj//fD7wgQ9062fvaV5BlyR1i7a2NiKCf/zjH4wfP57BgwezySabcMopp7B06dIuv+6rr77KddddxyGHHLLsIADgkEMOYfXVV+fqq69ud//K2fbKQQDAlltuye67797hvpLUKfNnOfxMpdHb83J1cV6xyy67APDYY4+1+x4DBgxg3XXXXS6+3sICXerNPBBQCb3vfe9jr7324le/+hUTJ05k8uTJnH766cvWL1mypKlHxf3338/LL7/MjjvuuNz7rLHGGmy99dbcfffd7cYzZ86cFfaF7Ix7R/tKUtMqk8U5/Ewl05fy8qxZswDYfvvtV1i3dOlSlixZwhNPPMHkyZO57777+PSnP93u65WRXdyl3mrExKxAnzfVGWNVKkcffTQnnngiAPvssw/PPfcc3/72tznuuONYb731mj6bnVICYMGCBQAMGTJkhW3WX3/9ZesbWbBgQcN9Fy5c2FQsktShkZMszlVKfSUvP//88xx33HFsv/32vO9971th/Ze+9CW+/e1vA7DWWmtx2WWX8c53vrPdWMrIAl3qrTwQUEkdcsghy/186KGH8pOf/IS77rqLPfbYg9tuu61Tr1c5IIiIhus6sjL7SpLUm/WFvLxkyRIOO+wwHnvsMW666SYGDFixjD3uuOM49NBDeeKJJ/jZz37GxIkTufLKK9l///2biqksLNAlSd1qo402qvtzZbzY6NGjO/V666+/PkDdM/ILFy5k1KhR7e4/ZMiQhvvWO4MvSVJf0tvz8tKlSzniiCOYNm0av/nNb9hpp53qvu5mm23GZpttBsD+++/PuHHj+MIXvtDrCnTHoEuSutWTTz5Z9+fhw4cDsNpqqzX1qNh6660ZOHAgc+bMWe51X375ZR544AF22GGHduMZNWrUCvsC3H333R3uK0mdtnB2dqeV2odzxqggvT0vf+ITn+DnP/95p7usjxkzhrlz5za9fVl4BV3q7SoHArVGTHRsugpx+eWXLxvrBnDZZZex1lprLZsQprNd6VZffXX2228/Lr/8ctra2pZ1a7vyyitZvHgxBxxwQLv7H3DAAXzhC1/ggQceYKuttgJg3rx53HTTTXzjG9/oVCyS1K7KRHG1Fs7Ons3LKkBvzsvHH388P/nJT7jwwgvrjjtvZOnSpdx4441svfXWnfpsZWCBLvVmHgiohH784x+zdOlSdtllF37/+9/zk5/8hLa2NtZbbz0gO6PdWW1tbey2224ccsghfPrTn2bevHl88Ytf5OCDD2bnnXdett3PfvYzPv7xj3PDDTcwduxYIJsc5/vf/z4HHnggZ5xxBhHBySefzBve8AaOOeaYbvnMkgRkebde7q13Il1qkd6al7/5zW9y1lln8fGPf5xtttmGW265Zdm6oUOHLiu+29raWLBgAbvvvjsbb7wxTzzxBOeddx5/+ctfmDq1983XZIEu9WYeCKiErr76ao499lgmT57Muuuuy0knncTJJ5+8Uq85evRofv/733PCCSfwnve8h3XXXZePfvSjfP3rX19uu6VLl/Laa68tN9HMmmuuyfTp0/nv//5vPvKRj5BS4p3vfCdnn302a6211krFJRVt9mwYN67oKHqviRNhkuey1cf11rz8u9/9DoCf/vSn/PSnP13udY844gguuOACAN7ylrdw9tlnc9lll/Hss8+y8cYb8+Y3v5k//elP7L777iv1OYsQzmLbd40ZMybdfvvt3fZ6lQOAmTO77SXVUyoF+l4zi4xC/UxbWxunnXYar776at3ZVdX9IuKOlFLnL32oEN2dl6dMgV54cag0Zs+G0aNbdFwzbVx2a9S3nmvvNrWMebn1uiMv+5uS+qr5s7IJaTwQkKQ+adIkr/6ujJb2PBgxMcvL86aalyW1y1ncpb6oMjbd+6RLklS8kZNg2Niio5DUC1igS32RBwIqQFtbGyklu9FJklQC5uXeyQJdkiRJkqQSsECXJEmSJKkELNAlSZIkSSoBC3RJkiRJkkrAAl2SJEmSpBKwQJckSZLUvjvbio5A6hecc19SOd3ZBju1FRyEJEndaOFsmDau6CheN2JidmvWZtx1mnlZagELdKkv80BAkqRyGDGx6AiWt3B29txsXpbUEhboLRARbwD+D9gbCGAacFxK6eEO9hsDTAL2BDYHngb+BJyUUnqwR4NW7+eBgCTVZV5WIUZOKlcOLNMJfEnLWKD3sIgYDEwHFgNHAAk4A5gRETullBa1s/uhwCjgu8AcYDhwMnB7RIxOKT3So8Grd/NAQJJWYF6WJJWZBXrPOxrYCtg2pTQXICLuBP4FHAOc1c6+30wpPVW9ICJuAh7MX/eUHolYkqS+y7wsVXR2KNy0cZ0briap0yzQe94BwC2VgwCAlNKDeUI/kHYOBGoPAvJlD0XEU2Rn7aXepbMHAnOneBAgqbuZlyVoPBTuhXnw4kP1182flT3+cszry3Y81TljpG5kgd7zRgFX11k+B/hgZ18sIrYHhgH3rGRcUmt15UDgL8csfxAAHghIWlnmZQk6PxRuasCwsdm/95rZIyFJskBvhfWBhXWWLwCGdOaFImIA8CPgKeC8lQ9NaqGuHgh4ECCpe5mXJUmltUrRAfQTqc6y6MLrfB94O3B4SqnewQURMSkibo+I2596aoWeeJIkybwsSSopr6D3vIVkZ+trDaH+Gfy6IuJ/yG7tckRK6fpG26WUpgBTAMaMGVPvAESSpP7MvCytjHrzyThxnNRtLNB73hyy8W61dgDubuYFIuKrwInAZ1NKF3VjbFK5NZpUzgMBSV1nXpa6qt58MgtnZ8/mZalbWKD3vGuA/42IrVJKDwBExAhgd7Lk3q6I+CzZ/Vm/mlL6Xk8GKpXOkNErLvNAQNLKMS9LXbHjqfXnk+nM3VkkdcgCvef9GPgMcHVEnEQ27m0y8AhwbmWjiNgCuB84PaV0er7sUOBs4DpgekTsWvW6z6WUmjrTL/VKjWZr90BA0soxL0td4R1UpJawQO9hKaVFETEB+D/gIrJJaG4AjkspvVC1aQCrsvzEffvly/fLH9VmAeN6KGypeB4ISOoB5mVJUplZoLdASulh4KAOtplHzQyyKaUjgSN7Ki5Jkvoj87Ikqay8zZokSZIkSSVggS5JkiRJUglYoEuSJEmSVAIW6JIkSZIklYAFuiRJkqSuWzg7uw3q3ClFRyL1ehboknqf+bM8CJAkqQxGTIQho7Mifd7UoqORej0LdEm9y4iJ2bMHAZIkFW/kJNhrZlakS1ppFuiSepeRk2DY2KKjkCRJkrqdBbokSZIkSSUwoOgAJEmSJPUBlcniIBuSNnJSkdFIvZIFuqTeqfogADwQkCSpSJU5YiDL0WBelrrAAl1S71N9EAAeCEiSVLSRk17Pw9Un0CV1igW6pN6n+iAAPBCQJElSn+AkcZIkSZIklYBX0CVJkiR1r9q5YppRxHwyd7bBTm3lfb0yurNtxc9Yb1lvdmdb/c/TaHk38gq6pL6hciDQmcfcKa2N8c623vGaZXNnW3PLeqs727q2TpLKasREGDK6c/ssnA3zpvZENO2767Ryv14Z1fuMfe1zN/o8LficXkGX1PvVThrXjCImlrvrtO4/69oTr1k29T5jX/rc7X2WvvQ5JfUftXPFNMP5ZCTAAl1SX+CBgCRJkvoAu7hLkiRJklQCXkGXJEmSVLyuTCzXHbr7PftDL716n7Gvfe5Gn6eHP6cFuqT+q4gDgZ54v76WEOvp6wcC7X2WvvQ5JamRrswn0xkvzIMXH6q/bv6sFZcN3gLWGtG61yujzn7G3vi5O/sZGy3vxs9pgS6pf+rJA4GeSNoeCDS3rMyfuysHAY3WlflzSlJXdGU+me4wNWBiKu/rlVG9z9jXPnejz9Ph54yVfmsLdEn9UxEHAj2RvPpaQqynrx8ItPdZWnAgIEmSysNJ4iRJkiRJKgELdEmSJEmSSsACXZIkSZKkErBAlyRJktQ/7XhquV+vjOp9xr72uRt9nhZ8Tgt0SWqVnvhS72sJsZ6+fiDQ3mfpS59Tkspop7Zyv14Z1fuMfe1zN/o8LficFuiS1Co98aXe1xJiPX39QKC9z9KXPqckSeqQBbokSZIkSSVggS5JkiRJUglYoEuSJEmSVAIW6JIkSZIklYAFuiRJkiRJJWCBLkmSJElSCVigS5IkSZJUAhbokiRJkiSVgAW6JEmSJEklYIEuSZIkSVIJWKBLkiRJklQCFuiSJEmSJJWABbokSZIkSSVggS5JkiRJUglYoEuSJEmSVAIW6JIkSZIklYAFuiRJkiRJJWCBLkmSJElSCVigS5IkSZJUAhbokiRJkiSVgAW6JEmSJEklYIEuSZIkSVIJWKBLkiRJklQCFuiSJEmSJJWABbokSZIkSSVggS5JkiRJUglYoEuSJEmSVAIW6JIkSZIklYAFuiRJkiRJJWCBLkmSJElSCVigS5IkSZJUAhboLRARb4iIKyPi2Yh4LiKuiojNm9x3jYj4VkQ8HhEvRcTNEbFnT8csSVJfZV6WJJWVBXoPi4jBwHRgO+AI4CPANsCMiFiziZc4DzgaOAXYH3gc+H1EjO6RgCVJ6sPMy5KkMhtQdAD9wNHAVsC2KaW5ABFxJ/Av4BjgrEY7RsSbgYnAx1NK5+fLZgFzgNOBA3o2dEmS+hzzsiSptLyC3vMOAG6pHAQApJQeBG4CDmxi31eBn1ftuwS4DNg3IgZ2f7iSJPVp5mVJUmlZoPe8UcBddZbPAXZoYt8HU0ov1tl3dWDkyocnSVK/Yl6WJJWWXdx73vrAwjrLFwBDVmLfyvrlRMQkYFL+4+KIqHcQslIiuvsVS2ND4Omig+glbKvOsb06x/Zq3rZFB9AL9bm83Ie15LugDx3X+N3ZPNuqc2yv5q10XrZAb41UZ1kz6SA6u29KaQowBSAibk8pjWkqQtlenWBbdY7t1Tm2V/Mi4vaiY+ilzMu9gO3VObZX82yrzrG9mtcdedku7j1vIXXOqJOdpa93Fr7agnb2rayXJEnNMy9LkkrLAr3nzSEbs1ZrB+DuJvbdMr8lTO2+rwBzV9xFkiS1w7wsSSotC/Sedw2wa0RsVVkQESOA3fN1He27GvDBqn0HAB8Crk8pLe5g/yldCbgfs72aZ1t1ju3VObZX82yrzjMv9x62V+fYXs2zrTrH9mreSrdVpFRvKJW6S0SsCfwdeAk4iWzs2mRgbWCnlNIL+XZbAPcDp6eUTq/a/zJgX+CLwIPAJ4H9gbenlP7awo8iSVKvZ16WJJWZV9B7WEppETABuA+4CLiELKFPqBwE5AJYlRV/Jx8DzgfOAH4DvAHYz4MASZI6z7wsSSozr6BLkiRJklQCXkHvYyLiDRFxZUQ8GxHPRcRVEbF50XEVLSI2i4jvRcTNEfFiRKR8zGH1Nu+MiIsj4v6IeCl//mFEDCso7EJExMER8YuIeChvh3sj4n8iYu129jk3b9OLWxlrGUTEvhExPSKeiIjFEfFoRFweETvU2XbXiLguIp6JiEUR8Y+IOLSIuMsib48UEWfULB+Vf3/9O2+rORFxfD7et1+JiHdHxB8j4oX8e/32iJhQtX5IRPwkIp7O22paRLypyJj1OvNyY+bm5pmbm2deXjnm5Y71dF62QO9DIptVdjqwHXAE8BFgG2BGZGPu+rORwCFkt9D5U4NtPgFsQNZtcT/gf4ADgFsiYq1WBFkSXwBeA75C1g4/JBtj+YeIWOE7IyLeDnwYeK6VQZbI+sAdwGeAfYAvk80QfUtkY1gBiIj3AH8EngAmAgcCPwbWaHXAZRERhwFvrrN8U2AmsBVwHPBe4FfAt4CvtSzAEoiIY4Cryf7G3k82OdkVwOB8fZBNXLYfcCxwENkkZjMiYrMiYtbrzMsdMjc3z9zcPPNyF5mXO9aSvJxS8tFHHsDnyL68R1Yt2xJYAny+6PgKbptVqv59FNmkQCNqthlaZ789820/XvRnaGFb1WuHj+btMKFm+WrAXWTJbx5wcdHxl+EBbJu31/H5z2sD84Gzi46tLA9gPbKDosPytjqjat2kfNkba/a5DHi86Nhb2EYjyCYyO66dbQ7M22p81bJ1ye7H/d2iP0N/f5iXO2wfc3PzbWVuXrn2My933Ebm5Y7bqCV52SvofcsBwC0ppWX3YU0pPQjcRPbH0m+llJY2sc1TdRbflj8P796IyquT7fBFskmUvt2jQfU+/8mfX82fPwgMxXaqdiYwJ6V0aZ11q+fPtVd+nqF/9fz6OLAU+FE72xwA/DulNKOyIKX0LPBr+vn3fkmYl9thbm6euXmlmZc7Zl7uWEvycn9q0P5gFNkZ01pzgBXG3agpY/PnewqNongrtENEbE12i6JPpZReKSSqEomIVSNi9YjYBjiX7Cz0ZfnqPcjOnL4pH9+2JCIeiYhTI2LVomIuSkTsQXbl51MNNrkCeBr4fkRsGRHrRMT7yboH96eDqT2AfwKH5uNul0TE3Ij4dNU27X3vb97PugCXkXm5Z5ibM+bmdpiXm2deblpL8nK/G9Tfx61PNo6r1gJgSItj6fXyiVfOJkt8vyo0mAJFxHDgdGBaSun2qlU/Aq6qPkPYz90K7Jz/ey5Zl8P5+c+bko1Nmkp2v+U7gL2Ak8m6lP13SyMtUESsRnag9L8ppXvrbZNSejIidiMb4/VAZTHQllI6szWRlsKm+eNbZONO7ye76vP9iBiQUvoO2ff+vDr7LsifhwAv1Fmv1jAvdzNzc8bc3BTzchPMy53Skrxsgd731LtvXrQ8il4un5HyUrJuY7unlJYUHFIh8rN8V5ONl/xY1fLDgV3IJj5S5iPAOmQTqHyBbOKePVJK88h6K60BfDWldFa+/cyI2AD4dES05d2f+oMTgEG0M6lMRAwFrgIWAQeTdU2cAJwUEYtTSt9sRaAlsArZOMkjU0pX5cumRzbL9Zcj4rtk3+9+75ebv59uYm7OmJubZl5ujnm5eS3JyxbofctCsrM2tYZQ/wy+6shnQ72Q7Ezqe1JKdxYcUiEiYg2yWSi3AsamlB7Nl68FnAV8E3g5ItbLd1kFWC3/eVFK6dUVXrQPSylVuhjeGhG/Izt7eiLZDMSVsW9/qNnt+nz9KODPLQizUJHdWuqrZJNBDYyIgVWrB+Z/O88DXyKbiGWLlFLlu2tm3u1wckScl1J6unWRF+Y/ZDN+1/u72Q/YhOyMfKPvffC7v2jm5W5ibs6Ym5tnXu6YebnTWpKXHYPet8wh+0KptQNwd4tj6c1+BHwIODSldEPRwRQh7+70C+CtwLtTSv+oWr0h2cQqXyf7kqk83sDrt8t5T0sDLpmU0jNk3elG5ovmVFbVbFo5m9rhREl9xFZkVywuZvm/HciubiwE3pQ/5lYdBFT8hWx24pH0D3MaLK/+u2nve//hlJLd24tlXu4+5mZzc5eZlxsyL3dOS/KyBXrfcg2wa0RsVVmQd7nYPV+nDkTEt8nOIn4spfSrgsMpRH6V4hLgncCBKaVbajZ5Ahhf5/EkMC3/940tC7iEImIjsi6G9+eLfpU/71ez6b7Ay9SfTKQvmk39vx3IDg7Gkx1APQGMjIjaMbpvy58f6/FIy+GX+fO+Ncv3BR5NKT1B9t0+PCIqk0UREeuQ3aPW7/3imZe7gbnZ3LyyzMsNzca83Bktyct2ce9bfgx8Brg6Ik4iOys4GXiEbPKHfi0iDs7/WZkw5F0R8RTwVEppVkScAHwe+Cnwr4jYtWr3p1JK99M/nEM24cXXgEU17fBo3p1uZu1OEfEy8GRKaYV1fVlE/BL4K3An2e1H3kg2ucwS8plNU0p3RcQFwOn5QdZfybppHgVM7i9XOfMrGDNrl0cEwEOVv52I+BHwYeD6iPgWWZeycWRn83+ZUnqkJQEX77fADODciNiQbGKeg4F9eH3c6TXAzcDFEfFFsqsdXyY7m9+fJu4pK/NyB8zNTTM3N8m83Dzzcqe1Ji935SbtPsr7ADYn6/70HNmYkV8BI4qOqwwPsgOjeo+Z+fqZ7WxzQdHxt7Cd5rXTDm0d7Hdx0fEX0F4nkM3++gzwInAv2YH3iJrtVgfOIDswfwW4D/hc0fGX4ZH/bZ1Rs2zXPBE+TjYpzRyyWwcNKjreFrfNOmQH5k/mfzd3AhNrtlmfrHhZkP8N3gC8uejYfSz7/ZiX228fc3Nz7WRubr6tzMsr34bm5cZt0+N5OfIXkSRJkiRJBXIMuiRJkiRJJWCBLkmSJElSCVigS5IkSZJUAhbokiRJkiSVgAW6JEmSJEklYIEuSZIkSVIJWKBL/UBEpCYe8yJiRP7vI0sQ84ia+MYVHVM9eZxtEbFVnXXzIuLibniPw2vaYsDKvqYkqTjm5Z5jXlZv5x+T1D/sVvPzL4G/A21VyxYDj+fb3t+asJpyBvAb4O6iA2lgBHAqcCPwQA+9x3Vkv5ejgP/qofeQJLWOebnnjMC8rF7MAl3qB1JKt1T/HBGLgadrl+fqLSvS/Q3i7DdSSk8DT0fEfkXHIklaeebl3s28rJ5kF3dJy9TrShcRF0TEoxExJiL+HBEvRcS9EfGefP3n8y5jz0XE1RExtOY1B0TElyPinxGxOCL+HRHfjog1ViLOmRFxY0TsFxGz85j+FhFvy9/v6xHxeEQsyONfs2b/TSLiZxHxdB7TnRFxeM02R+ZtsWtEXJJ/vn9HxHcrsefd+2bku/yhUbe/iDg0Iu6JiEURcXtE7FGzfpeI+ENE/CciXoyIByLiB11tH0lS32BeXm4b87L6Ba+gS2rGOsDPgP8F/g18FfhFRJwDvBH4NLARcDZwDnBI1b4XA+8Fvgn8GdgemEzWBe2glYhpJPAt4GvAC8CZwDX5YwBwZP5e3wLmA18CyA8KZgFDgK8AjwCHAxdFxOCU0pSa97kIuBT4AFl3tjZgIVn3ub/mn/0c4LPAbfk+1d3+3gFsC5wMvJx/9msjYkRK6ZmIWAv4PfCXPObnydrm7V1uGUlSX2deNi+rj7JAl9SMtYFPpJT+CBAR/yYbK7c/sENK6bV8+Y7AsRGxakrptYh4B/Ah4IiU0s/y15oWEQuAiyNidEppdhdj2gB4e0rpgfy9VwGuBrZMKe2Vb/P7iNgT+CD5gQDwMWAbYHxKaWa+7HcRsRFwRkScV/k8uakppVOrYn8bcBhwakrpuYioJP17GnT5WwcYnVJamMf5BNkBw7uBqcB2ZAclX0op3Vm13wWdbA9JUv9hXn49dvOy+hS7uEtqxqLKQUDun/nztJqk+U+yE3+b5D/vB7xCdlZ/QOUBXJ+v33MlYrqvchBQE9Pva7b7J7BZRETVez5WdRBQcTEwFNihZvlvan7+B7B5J+K8uXIQULU/Va/xL+AZ4NzIZoV9QydeW5LUP5mXX2deVp9igS6pGc9U/5BSeiX/58Ka7SrLK+PYhgGrk3V1e7XqMT9fv8FKxNTovestHwCsmv+8PtmsuLWeqFpfbUHNz4uBgc2Hufz+KaXF+T/XyH9+FhhP1kXxB8DDEXFXRKxMN0NJUt/2TPUP5mXzsvoOu7hL6kn/IRvf9Y4G6//dwlgqFpCNPau1cf78nxbGAkDenfCg/CrGGODLwOUR8eaU0l2tjkeS1GeZl5tgXlaRvIIuqSddR3ZGet2U0u11HkUcCMwi61q3e83yiWRXEO7p5OtVzrwPWtnAUkpL8vFyJ5N9P2+/sq8pSVIV83InmJdVBK+gS+oxKaWZEXEpcGVEnEU2I+pSstlQ3w2ckFK6r8VhXQB8DrgqIr4KPAp8GNgbOKZm7F4z7gOWAB/PJ9lZDNybUnq+mZ0jYn9gEvAr4EFgTbKZZ58Hbu5kLJIkNWRe7ph5WUWzQJfU0w4HjgU+TnYbmMXAPLJJY55sdTAppUURMZbs9i/fIJsJ917gIymli7vwev+JiM8AJ5BdBViVbOzazCZf4l/AS2Rn5zchOwC4Ddg7pfRoZ+ORJKkD5uX2mZdVqEgpFR2DJK0gIkaQnbn+L7J7vb6W+ukXVj7T7arAKWQHDKullJYUG5UkqT8xL7/OvKye5Bh0SWV3HtkMs2OLDqRAHyZrg5OLDkSS1O+Zl83L6kFeQZdUShGxOrBT1aKmx4/1NRGxPrBV5eeU0u0FhiNJ6ofMy68zL6snWaBLkiRJklQCdnGXJEmSJKkELNAlSZIkSSoBC3RJkiRJkkrAAl2SJEmSpBKwQJckSZIkqQQs0CVJkiRJKgELdEmSJEmSSsACXZIkSZKkErBAlyRJkiSpBLpUoEfEuIhIEbFXg/WPRsQFdbYf14X3asv3HdCVWNt53SPz1x3RiX0uiIhHuzOOroqIHSPi3Ii4IyJeiYhUdExlFxFviIgrI+LZiHguIq6KiM2LjquMbKvOsb06x/Zqnm3Ve0TE0RHx24h4LCIWRcRdEfHFiFi96NjKKiI2i4jvRcTNEfFiZ4/L+hO/CzrH9uoc26t5rWirVl1B/yuwW/5cFr8hi+nxogPpop2BdwMPA7cXHEvpRcRgYDqwHXAE8BFgG2BGRKxZZGxlY1t1ju3VObZX82yrXucU4Angc8D+wM+BycAlRQZVciOBQ4CFwJ8KjqW0/C7oHNurc2yv5rWqrbr1qnQjKaXngFta8V4diYjVgCUppaeAp4qOZyVclFK6ECAiziA72aDGjga2ArZNKc0FiIg7gX8BxwBnFRhb2dhWnWN7dY7t1Tzbqnd5S35sUTEjIgI4LSK2Sik9UFRgJfbHlNJGABFxFLBPwfGUld8FnWN7dY7t1byWtFVLrqDX6+IeEatGxBkR8XjerWl6RGyXb9dW52W2jIjfRMQLEfFQRJwSEe3GHxEj8tf7VEScGRH/BhYD60WdLu4RMTEi/pa/x7MR8Y+IOKaD9/hYRLwaESd2oklWWkppaSvfrw84ALil8p8JIKX0IHATcGBhUZWTbdU5tlfn2F7Ns616kZrivOK2/Hl4K2PpLTyWaZrfBZ1je3WO7dW8lrTVyhboq0TEgNpHk/ueBnwF+BnZB/o9cE072/+SrEvB+4Bf5fsf0eR7fRV4IzAJeD/wcu0GEbEHcDEwK3+PDwI/BtZr9KIR8WXgXODolNI3moxFxRgF3FVn+RxghxbHUna2VefYXp1jezXPtur9xgJLgfuKDkS9mt8FnWN7dY7t1byWtNXKdnH/fVd2ioghwHHAj1JKJ+SL/xARrwLfbrDbt1NK5+f/nhYRE4DDgPMbbF/tSeD9KaVlE6llvc6WsyvwTErpuKpl1zeIfxXgO8DH89f9TRMxqFjrk41xq7UAGNLiWMrOtuoc26tzbK/m2Va9WETsRDYe/acppSeLjke9mt8FnWN7dY7t1byWtNXKXkH/NLBLncf8DvZ7E7AmcEXN8ivb2ae2CL4LaHbGvF9VF+cN3AYMiYiLI2L/iFivwXYDgMuAicBeFue9Sr2/gRXO1AiwrTrL9uoc26t5tlUvFBGbAFcD9wOfLzgc9Q1+F3SO7dU5tlfzerytVvYK+n0ppRVmEM+vhLdnk/y5tpBv7wzzgpqfFwNrdPA+FR3O1J5SmhURHwSOJetOT0TMAj6fUrqzatN1gPeQdbf/S5Pvr+ItJDvrVWsI9c+E9We2VefYXp1jezXPtuqFImID4A9kB2z7ppSeLzgk9X5+F3SO7dU5tlfzWtJWrbrNWq1KwTysZvlGPfR+Td0jPKV0ZUppLFkjv5/sRMJ1NZPRLSAr0McDl3ZizL2KNYds3EitHYC7WxxL2dlWnWN7dY7t1TzbqpeJiHXIhv9tQNbL7rGCQ1Lf4HdB59henWN7Na8lbVVUgf4PYBHZRGzVan8uRErphZTStWQTwG1Clmir188E3pU/LrNI7xWuAXaNiK0qC/IZ/Hen/ckJ+yPbqnNsr86xvZpnW/Ui+f1xfwNsCexTPcuvtJL8Lugc26tzbK/mtaStCinQU0oLgbOBT0TENyNi74g4gez+cZDNeNpSEXF6RJwbEYdGxJ4RMRH4LDC73q1TUkp/AvYju2fnz/P7q7cy3sERcXBEHAxsly87OH+MaWUsvcSPgXnA1RFxYEQcQDY+8BGyEzF6nW3VObZX59hezbOtepdfkB2ktQFrRsSuVY+hxYZWXlXHMjvni96VLxtbZFwl43dB59henWN7Na8lbVXUFXSAU4H/IbtV2jVkV6OPzNc9W0A8twIjgP8jGzv2TbJbrr2n0Q4ppZuAfYG9gCsiYvWeD3OZYWST7F0BHJQvq/z8mRbG0SuklBYBE8hudXMRcAnwIDAhpfRCkbGVjW3VObZX59hezbOtep39yMadfxe4uebR8FhCy45dPpH//IP859MKi6hk/C7oHNurc2yv5rWqraLjyc1bJ5+k7XJgz/wKtSRJkiRJ/UJhBXpEvI3sjPKtwMtkXZtOBO4F3t7EbdEkSZIkSeozipzc7AVgT7J7qa9Ddsu1y4EvW5xLkiRJkvqbUnVxlyRJkiSpvypykjhJkiRJkpSzQJckSZIkqQS6tUCPiBERkSLiqO583Z4SEfMi4oJObN/jny8iLsjfo97jnzXbbhkRV0bEMxGxKCJm9Ld7oEfEvhExPSKeiIjFEfFoRFweETvUbDckIn4SEU/nbTUtIt5UVNxFiYjNIuJ7EXFzRLyY/12NqNmm6b9BZSJifETcGBEvRcSCiLgoIjYqOq4yioiZ7fx9XVd0fGXS7Pebek5V3m/0OLRm+yMi4o6IeC4inoqIP0TEO4qKvwj5Pcx/EREP5d+J90bE/0TE2jXbjY6I6yLihby9romIkUXFXRTzcs8wLzfPvNw5rcjNRU4SVwbvB54rOogak4Ef1SwbAVxKdr94ACJiA+BG4HngGOBF4PPAjIh4a0rpnpZEW7z1gTvI7pv6FLA52d0AbomIN6WUHoqIIGu7LYFjgYXAl8naanRK6dFiQi/ESOAQsjb7E7BPnW2a+htUJj/4vh74PXAQsAFwBnBDROycUlpcZHwl9CmyiUGr7QachX9ftTr8fisyuH7icbK/z1pnAHuQ/d8HICImAeeSfX+eCAwmy8t/iIjdUkp/6/lwS+ELwMPAV4BHgf8HtAHjI+LtKaWlEbENWQ66C/gw2fHoqcAf87w8v5DIi2Fe7mbm5U4zL3dOz+fmlFK3Pci+LBJwVHe+bnc/gIG96fMBJ+fvO6pq2UnAEmBk1bI1gSeBy4tu44J/v9vm7XV8/vOB+c/jq7ZZF1gAfLfoeFvcNqtU/fuovF1GNLHfCn+DPpa1zTRgLjCgatkueXt9quj4esMDOA9YDKxfdCxlf9R+v/ko5HcwmOzk/hU1y/8M/Llm2drAK8D/FB13C9tnaJ1lH83/bifkP/8EeAZYr2qbzchuu3tm0Z+hxe1lXu7+NjUvr3wbmpc7117dmpsLGYMeEYdHxN8j4uXIuhxfFBGbVK3/fkTMrdnnjryrxciqZV+LiPn5FdJG73VB3vVgt4j4c0S8BJyZr1uui3tEbBwRF0bEv/MuC49HxLURMayd198wIm6NiHsiYvOutUiHPgrckVKaU7VsV+BfKaVl7ZRSWkR29nX/iOjPvSP+kz+/mj8fAPw7pTSjskFK6Vng12TFe7+RUlraxV3r/Q0qsyvwh5TSksqClNJtZH+H7y8sql4iIgYBHwR+nVJaUHQ8vUDt95ta7wNkhfeFNctXZ8VeeS+S/a76zZw/KaWn6iy+LX8enj/vCtycUnqmar9Hya6o96vvTfNyjzAvrwTzcpd0a25uecLIu4BdBNxDluROBPYFZkXEWvlm04GtKwVvRAwBRgMvAROqXm4CMCPlpy7asS5wGVlXoHcBUxtsdxFZl44vAnsDnyXrnjW4wWcZAdxEdsZkj5TSwx3E0WkRsTtZ96faA4HXyM7K11oMDAK27u5YyiwiVo2I1fNuc+cCT5D9zgFGkSX9WnOAzav+7lRHO3+DyrT3f3HHFsfSGzUqdpTr4PtNrXcEMB+oHZv5A2CviPiviFgvIoYD3yc7YDuvxTGWzdj8uTL8rr3vza0jYo2WRNVLmZc7ZF5eOeblJvRkbm7pVdaIWJVsHM3MlNKhVcv/SXbl9+PAd4GZ5F2Syf44xpKdlb4qXzYlL6rG0Nwfz1rA4SmlqzvYbjfgKymlS6qWXdHgs7wZ+B0wGzg4pfRiE3F0xUfJkvulNcvvBfaOiA1SSv/JY1oFeGu+fv0eiqesbgV2zv89l6wbXWUM2/rAvDr7VM4KDgFe6NHoerdGf4PK3Et2tn6ZiNgC2ASvcjbjo2TFzu+KDqTE2vt+UwvlRfcE4DvVV+cAUko/zTv0/YCsCzdkB2x7p5Tua2mgJZK32enAtJTS7fnie4G3R8RqKaVX8+3WJjuhHmR5+fEi4u0lzMvtMy+vHPNyc3osN7f6Cvq2wDCgugAmpXQj8BD5Gda8O8WdvH61fAIwi2xMyfh82Z5kJximN/G+S4Brm9juNuCLEfG5iHhTO13n96yK54CeKs4jYiDZxCHXppSerln9I7Lf388iYut8iMB3ySZCA+hql6ne6iNkX8YTyU7m/CFenwU1yE741Go4NEKZDv4GlfkO8NaIOCMihkXEdmS9cZbS//4fdkpEbArsBVxSW+xoOe19v6m1PkKWe1e4OBARBwLnkF1J2Qt4L1nvrd9GRL+8apdfTLma7DjsY1WrvkPW3f1HETE8L57OJ7ugAn53NmRebop5uYvMy53SY7m51QV65apuvbOiT7D8Vd/pvF6Mjwdm5I+N8mnsx5ONK27mrPT8lNJrTWz3IbLZCr9EdoLgsYg4Jb8yXe3dZEnk3B7+4z0QWI86BwIppQfIZj7dmeyszb/JegD8X75JvzrznFK6J6V0a0rpUuCdZL+fE/PVC6jfo2BI/rywBSH2Vg3/BpXJe9ycARxPNknj3cBjwG/pZ/8Pu+BwGhQ7el0H329qrY8Cs1NKf69emJ/QnwJcmVL6XErphpTStcB7yO62Mrn1oRYr76Z+DbAVsG+qumNKSukm4NPAwWRDCefxeq55hdd7uGlF5uUOmJdXinm5ST2Zm1tdoFe+cDeus25jXh9gD1kx/oaI2I2sy9P0lNITZOOXJuSPGSu8Sn0djVHPNkppfkrp0yml4cB2wAXAaWS3Mat2MtkEY7/LxwH1lCOAp8m+UFaQUvoF2RnoHchmc9+Z7I/jkZ4YD99b5JPOzCUbnwXZWPNRdTbdAXg4pWT39sba/RtUJqV0MrAhsBOwSUrpMGAbslshqrGPAn+vLXbUWJ3vN7VIROwCbE/9A9eNyHoI3la9MKX0CvD3fL9+IyJWA35BNuzu3Smlf9Ruk1L6AVmb7QhsnlLaC9gUuLXS7V11mZebYF7uMvNyF3R3bm51gX4v2ZmsQ6sXRsTbgS3Iuo1X/JFskofJZF9ElUm+ppNNXjCa5rq3d0lK6d6U0lfIrq7Wdk17lax70fXAdfn9FrtVRGxEdi/Mqe0lqpTSa/kZnPvzbikfAn7Y3fH0JnnbbQfcny+6BhgeEWOrtlmHrPuh93dsoNm/QWVSSotSSv9IKT0ZEfuR/Q3W3rdWuYgYQ3bizLP0nVDn+02tcwRZV+16E80uJJuA6q3VCyNidbLjlcd6OriyyHsdXkJ2RenAlNItjbZNKS1OKc1JKT0SEW8i61rbr49h2mNe7hzzcueYl7uuu3NzT00St3NEPFNn+TXAKcC5EXExcDHZFeCvAf8iG38EZLfBioi/kn3BX1E1U/sMsm5RlX93i4hYl2xM+SXAP8mK8APJukFfX7t9SunViDg03/53EfGelNKs2u1WwofJfj91/5PkZ6fPJDup8RzZf6gvk10t/nY3xlFqEfFL4K9kQxKeA94I/DfZQVSlHa4BbgYujogvkh1IfZlsDPqZrY65aBFxcP7PysQW74qIp4Cnav6G2/0bVCYi/h/Z3SH+mi/ag+xOEGemlP5cWGDl91EaFzui6e83tUCecw8FfldvEqCU0uKI+DHwmYhYQDbvzSDgM8AI4LjWRVu4c8hu0fQ1YFFEVE/W9WhK6dGI2Az4JNm94xeT5aOvAFfl3UX7FfNy9zIvd5l5uQktyc3N3Cy92QdZEkrtPDbMtzucrMvXYrJu7ReRdT+pfb1v5vt9omrZ+mQTPMxrMqYLyBJCvXXzgAvyfw8km9hlDtmM3s+RdVWbWOfzHVW1bFWyP+RFwPhubMu/A/9oZ/0AsgOAJ/N2vJ9svM3g7vydlv0BnADcATxDdr/Ze/Pf44ia7dYHfko2zOJF4AbgzUXHX1CbNfr/ObNmu3b/Bn0sa6dRZF3mniG7FeRfgY8VHVeZH8BqwFNk91gtPJ6yPpr9fvPRkt/F+/PvyYPa2WYAWUE+m2zc+VNkd6XZp+j4W9xW89rJM235NhuRXRR5Oj+GuZtsvPCAouMvqM3My93bnublzreZebn5turx3Bz5G0mSJEmSpAK1egy6JEmSJEmqwwJdkiRJkqQSsECXJEmSJKkELNAlSZIkSSoBC3RJkiRJkkrAAl2SJEmSpBKwQJckSZIkqQQs0CVJkiRJKgELdEmSJEmSSsACXZIkSZKkErBAlyRJkiSpBCzQJUmSJEkqAQt0SZIkSZJKwAJdkiRJkqQSsECXJEmSJKkELNAlSZIkSSoBC3RJkiRJkkrAAl2SJEmSpBKwQJckSZIkqQQs0CVJkiRJKgELdEmSJEmSSsACXZIkSZKkErBAlyRJkiSpBCzQJUmSJEkqAQt0SZIkSZJKwAJdkiRJkqQS+P+wpHEJ6mU/DwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1008x504 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Display final KM plots based on train and validation ensemble predictions \n",
    "#plot KM\n",
    "size=16\n",
    "data_train= risk_ensemble_train\n",
    "data_val= risk_ensemble_valid\n",
    "ix= data_train[\"risk_group\"] ==1\n",
    "ix2= data_val[\"risk_group\"] ==1\n",
    "#plotting images#\n",
    "plt.rc('xtick', labelsize=size)    # fontsize of the tick labels\n",
    "plt.rc('ytick', labelsize=size)\n",
    "fig,ax = plt.subplots(nrows=1, ncols=2,figsize= (14,7))\n",
    "\n",
    "#training kaplan meier\n",
    "kmf_train_high = KaplanMeierFitter()\n",
    "ax[0] = kmf_train_high.fit(data_train.loc[ix][\"time\"], data_train.loc[ix]['status'], label='High risk').plot(ax=ax[0],\n",
    "                            show_censors=True,ci_show=False, xlim=(0,60), color= \"blue\",yticks=(0.0,0.2,0.4,0.6,0.8,1),\n",
    "                            ylim=(0.0,1),xticks=(0,12,24,36,48,60))\n",
    "\n",
    "kmf_train_low = KaplanMeierFitter()\n",
    "ax[0] = kmf_train_low.fit(data_train.loc[~ix][\"time\"], data_train.loc[~ix]['status'], label='Low risk').plot(ax=ax[0],\n",
    "                             show_censors=True,ci_show=False, xlim=(0,60),color=\"orange\",yticks=(0.0,0.2,0.4,0.6,0.8,1),\n",
    "                             ylim=(0.0,1),xticks=(0,12,24,36,48,60))\n",
    "\n",
    "\n",
    "#validation kaplan-meier\n",
    "kmf_val_high = KaplanMeierFitter()\n",
    "ax[1]= kmf_val_high.fit(data_val.loc[ix2][\"time\"], data_val.loc[ix2]['status'], label='High risk').plot(ax=ax[1], show_censors=True,ci_show=False, xlim=(0,60), color= \"blue\",yticks=(0.0,0.2,0.4,0.6,0.8,1), ylim=(0.0,1),xticks=(0,12,24,36,48,60))\n",
    "\n",
    "kmf_val_low = KaplanMeierFitter()\n",
    "ax[1] = kmf_val_low.fit(data_val.loc[~ix2][\"time\"], data_val.loc[~ix2]['status'], label='Low risk').plot(ax=ax[1], show_censors=True,ci_show=False, xlim=(0,60),color=\"orange\",yticks=(0.0,0.2,0.4,0.6,0.8,1), ylim=(0.0,1),xticks=(0,12,24,36,48,60))\n",
    "#settings for training km\n",
    "ax[0].set_xlabel(\"Time [months]\",size=size)\n",
    "ax[0].set_ylabel(\"Local Recurrence\",size=size)\n",
    "ax[0].text(45, 0.25, \"p=\"+str(round(optimal_p_value,3)), size=size)\n",
    "#settings for validation km\n",
    "ax[1].set_xlabel(\"Time [months]\",size=size)\n",
    "ax[1].text(45, 0.25, \"p=\"+str(round(p_val_valid,3)), size=size)\n",
    "\n",
    "\n",
    "\n",
    "ax[0].tick_params(axis='both', which='major', labelsize=size)\n",
    "ax[1].tick_params(axis='both', which='major', labelsize=size)\n",
    "ax[0].legend(fontsize=size, loc=1)\n",
    "ax[1].legend(fontsize=size, loc=1)\n",
    "ax[0].set_title(\"Training\", size=size)\n",
    "ax[1].set_title(\"Internal Validation\",size=size)\n",
    "\n",
    "from lifelines.plotting import add_at_risk_counts\n",
    "util.add_at_risk(kmf_train_high, kmf_train_low, ax=ax[0],rows_to_show=['At risk'])\n",
    "util.add_at_risk(kmf_val_high, kmf_val_low, ax=ax[1],labels=[\"\", \"\"],rows_to_show=['At risk'])\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"Train_Valid_starta.png\", dpi=300)\n",
    "#dpi stands for dots per pixel per inch, good quality is above 300"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "flexible-crazy",
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Store final results in text file including\n",
    "1- train and valid p-values at optimal cutoff\n",
    "2- Median concordance index of all cv folds and concordance index with ensemble predictions\n",
    "'''\n",
    "import sys\n",
    "sys.stdout = open(\"Train_Results.txt\", \"w\")\n",
    "print(\"Median CI Train (All Subvolumes) %.2f\" % np.median(CI_train_all))\n",
    "print(\"Median CI Valid (All Subvolumes) %.2f\" % np.median(CI_valid_all))\n",
    "print('######################################')\n",
    "print(\"CI Ensemble Train (All Subvolumes) %.2f\" % CI_ensemble_train)\n",
    "print(\"CI Ensemble Valid (All Subvolumes) %.2f\" % CI_ensemble_valid)\n",
    "print('######################################')\n",
    "print(\"p_value at optimal cutoff %.3f on train data %.5f\" % (optimal_cutoff, optimal_p_value))\n",
    "print(\"p_value at optimal cutoff %.3f on valid data %.5f\" % (optimal_cutoff, p_val_valid))\n",
    "sys.stdout.close()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "deep",
   "language": "python",
   "name": "deep"
  },
  "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}