--- a +++ b/mlp/mlp_classifier.ipynb @@ -0,0 +1,1263 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Imminent ICU Admission and Prolonged Stay Prediction using Neural Networks" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Imports & Inits" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T14:23:54.544191Z", + "start_time": "2019-08-10T14:23:54.531466Z" + } + }, + "outputs": [], + "source": [ + "%load_ext autoreload\n", + "%autoreload 2" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T14:23:56.643514Z", + "start_time": "2019-08-10T14:23:54.545371Z" + }, + "scrolled": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "{'workdir': PosixPath('../data/workdir/mlp'),\n", + " 'dataset_csv': PosixPath('../data/proc_dataset.csv'),\n", + " 'cols': ['hadm_id',\n", + " 'imminent_adm_label',\n", + " 'prolonged_stay_label',\n", + " 'processed_note',\n", + " 'charttime',\n", + " 'intime',\n", + " 'chartinterval'],\n", + " 'imminent_adm_cols': ['hadm_id', 'processed_note', 'imminent_adm_label'],\n", + " 'prolonged_stay_cols': ['hadm_id', 'processed_note', 'prolonged_stay_label'],\n", + " 'dates': ['charttime', 'intime'],\n", + " 'device': 'cuda:2',\n", + " 'start_seed': 127,\n", + " 'min_freq': 3,\n", + " 'batch_size': 128,\n", + " 'hidden_dim': 100,\n", + " 'dropout_p': 0.1,\n", + " 'lr': 0.001,\n", + " 'wd': 0.001,\n", + " 'max_lr': 0.1,\n", + " 'max_epochs': 100,\n", + " 'ia_thresh': 0.2,\n", + " 'ps_thresh': 0.25}" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import sys\n", + "sys.path.append('../')\n", + "\n", + "import pdb\n", + "import matplotlib.pyplot as plt\n", + "import seaborn as sns\n", + "sns.set_style(\"darkgrid\")\n", + "%matplotlib inline\n", + "\n", + "import numpy as np\n", + "np.set_printoptions(precision=5)\n", + "\n", + "import pandas as pd\n", + "import pickle\n", + "\n", + "import torch\n", + "from torch import nn\n", + "from torch import optim\n", + "from torch.nn import functional as F\n", + "from torch.utils.data import DataLoader, TensorDataset\n", + "\n", + "from sklearn.feature_extraction.text import TfidfVectorizer\n", + "\n", + "from skorch import NeuralNetBinaryClassifier\n", + "from skorch.toy import MLPModule\n", + "from skorch.dataset import CVSplit\n", + "from skorch.callbacks import EarlyStopping, LRScheduler, Checkpoint\n", + "\n", + "from mlp_model import MLPModel\n", + "from utils.splits import set_group_splits\n", + "from utils.metrics import BinaryAvgMetrics, get_best_model\n", + "from utils.plots import *\n", + "\n", + "from args import args\n", + "vars(args)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## NN Dev" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T14:23:58.527239Z", + "start_time": "2019-08-10T14:23:56.645116Z" + } + }, + "outputs": [], + "source": [ + "seed = 643\n", + "full_df = pd.read_csv(args.dataset_csv, usecols=args.cols, parse_dates=args.dates)\n", + "ia_df = full_df.loc[(full_df['imminent_adm_label'] != -1)][args.imminent_adm_cols].reset_index(drop=True)\n", + "ps_df = full_df.loc[(full_df['chartinterval'] != 0)][args.prolonged_stay_cols].reset_index(drop=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Imminent ICU Admission" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T14:23:58.587967Z", + "start_time": "2019-08-10T14:23:58.528632Z" + } + }, + "outputs": [], + "source": [ + "ori_df = set_group_splits(ia_df.copy(), group_col='hadm_id', seed=seed)\n", + "df = ori_df\n", + "# df = ori_df.sample(1000).reset_index(drop=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T14:24:24.952218Z", + "start_time": "2019-08-10T14:23:58.589218Z" + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "((50809, 4), (42683, 60000), (8126, 60000), (42683,), (8126,))" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "vectorizer = TfidfVectorizer(sublinear_tf=True, ngram_range=(1,2), binary=True, max_features=60_000)\n", + "\n", + "x_train = vectorizer.fit_transform(df.loc[(df['split'] == 'train')]['processed_note']).astype(np.float32)\n", + "x_test = vectorizer.transform(df.loc[(df['split'] == 'test')]['processed_note']).astype(np.float32)\n", + "\n", + "x_train = np.asarray(x_train.todense())\n", + "x_test = np.asarray(x_test.todense())\n", + "\n", + "vocab_sz = len(vectorizer.vocabulary_)\n", + "\n", + "y_train = df.loc[(df['split'] == 'train')]['imminent_adm_label'].to_numpy()\n", + "y_test = df.loc[(df['split'] == 'test')]['imminent_adm_label'].to_numpy()\n", + "df.shape, x_train.shape, x_test.shape, y_train.shape, y_test.shape" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T14:25:10.219798Z", + "start_time": "2019-08-10T14:25:09.574881Z" + } + }, + "outputs": [], + "source": [ + "train_ds = TensorDataset(torch.tensor(x_train), torch.tensor(y_train.astype(np.float32)))\n", + "train_dl = DataLoader(train_ds, batch_size=args.batch_size, shuffle=True, drop_last=True)\n", + "itr = iter(train_dl)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T14:25:10.376434Z", + "start_time": "2019-08-10T14:25:10.222543Z" + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "tensor(0.6790, grad_fn=<BinaryCrossEntropyWithLogitsBackward>)" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "clf = MLPModule(input_units=vocab_sz, output_units=1, hidden_units=100, num_hidden=1, dropout=args.dropout_p, squeeze_output=True)\n", + "\n", + "loss_fn = nn.BCEWithLogitsLoss()\n", + "x, y = next(itr)\n", + "y_pred = clf(x)\n", + "\n", + "loss_fn(y_pred, y)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T14:25:11.452893Z", + "start_time": "2019-08-10T14:25:11.426771Z" + } + }, + "outputs": [], + "source": [ + "reduce_lr = LRScheduler(\n", + " policy='ReduceLROnPlateau',\n", + " mode='min',\n", + " factor=0.5,\n", + " patience=1,\n", + ")\n", + "\n", + "checkpoint = Checkpoint(\n", + " dirname=args.workdir/'models/dev3',\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T14:25:11.480936Z", + "start_time": "2019-08-10T14:25:11.454266Z" + } + }, + "outputs": [], + "source": [ + "net = NeuralNetBinaryClassifier(\n", + " clf,\n", + " max_epochs=args.max_epochs,\n", + " lr=args.lr,\n", + " device=args.device,\n", + " optimizer=optim.Adam,\n", + " optimizer__weight_decay=args.wd,\n", + " batch_size=args.batch_size,\n", + " verbose=1,\n", + " callbacks=[EarlyStopping, checkpoint, reduce_lr],\n", + " train_split=CVSplit(cv=0.15, stratified=True),\n", + " iterator_train__shuffle=True, \n", + " threshold=args.ia_thresh,\n", + ")\n", + "\n", + "net.set_params(callbacks__valid_acc=None);" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T14:23:12.121340Z", + "start_time": "2019-08-10T14:18:03.575272Z" + } + }, + "outputs": [], + "source": [ + "# net.fit(x_train, y_train.astype(np.float32))" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T14:25:28.804887Z", + "start_time": "2019-08-10T14:25:22.921517Z" + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "<class 'skorch.classifier.NeuralNetBinaryClassifier'>[initialized](\n", + " module_=MLPModule(\n", + " (nonlin): ReLU()\n", + " (sequential): Sequential(\n", + " (0): Linear(in_features=60000, out_features=100, bias=True)\n", + " (1): ReLU()\n", + " (2): Dropout(p=0.1)\n", + " (3): Linear(in_features=100, out_features=1, bias=True)\n", + " )\n", + " ),\n", + ")" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "net.initialize()" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T14:25:32.359374Z", + "start_time": "2019-08-10T14:25:32.263137Z" + } + }, + "outputs": [], + "source": [ + "net.load_params(checkpoint=checkpoint)" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T14:25:34.483957Z", + "start_time": "2019-08-10T14:25:34.160864Z" + }, + "scrolled": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "<matplotlib.legend.Legend at 0x7f388e84a4e0>" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlkAAAHSCAYAAADfZ97BAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd5hU9d3+8feZtrszW2YbvbdDs4AFBVEUpYiVoKKxPb/okzyJMUaNmmJi1BijURNrYkliw46KgiIqiBosFEEphyqwFGnbd2ennd8fu+ACu+ywO7uzO3u/rosL5sw5cz6zH2a4OeX7NWzbRkRERETiy5HoAkRERESSkUKWiIiISDNQyBIRERFpBgpZIiIiIs1AIUtERESkGShkiYiIiDQDV6ILOFA0GrUjkeYfVsLpNGiJ/UhiqL/JTz1Ofupx8kuGHrvdzl1Afl3PtbqQFYnYFBVVNPt+/H5vi+xHEkP9TX7qcfJTj5NfMvQ4Pz9jY33P6XShiIiISDNQyBIRERFpBgpZIiIiIs2g1V2TJSIiIvERiYQpLNxJOBxMdCl1+u47g7Yyh7LL5SE7Ox+nM/bopJAlIiKSpAoLd5Ka6sXn64RhGIku5yBOp4NIJJroMhpk2zbl5SUUFu4kL69zzNvpdKGIiEiSCoeD+HyZrTJgtSWGYeDzZR72EUGFLBERkSSmgBUfjfk5KmSJiIiINAOFLBEREWkWpaWlTJ/+ymFvd+ON11JaWnrY2/3pT7cxd+77h71dc4npwnfTNCcAfwecwJOWZd1dz3pTgFeA4yzLWmiapht4Ehhes69nLMv6c1wqFxERkVatrKyU119/hcmTL9hveSQSwel01rvdX//6YHOX1iIaDFmmaTqBR4AzgALgS9M0Z1iWteKA9TKAa4HPay2+AEixLOsI0zS9wArTNF+wLOvbeL0BERERadjM5d8x45vtcX3Nc4Z2YtKQjvU+/49/PMSWLVu48spLcLlcpKWlkZubx9q1q3nuuVe4+ebr2b59O8FgkAsumMq5504GYMqUs3nyyWeprKzgxhuv5cgjj+brr5eRn5/P3XffR0pKaoO1LVz4BY888jcikQgDBw7mxht/jcfj4bHHHuLTT+fjdDo57rgTuOaa6/jww/f5978fx+Fwkp6eziOPPBGXn08sR7KOB9ZalrUewDTNF4FzgRUHrHcHcA9wY61lNuAzTdMFpAFBoKSpRYuIiEjr95Of/Jz169fxn/9MY/Hihdx003U888xLdOnSFYDf/OYPpKdnUFUV4KqrLmfMmNPIyvLv9xoFBZu57bY/cfPNv+PWW29h3rwPGT/+zEPut6qqirvu+iN/+9uj9OjRkzvu+D1vvPEqEyZMYv78uUyb9hqGYew7Jfmf/zzB/fc/TH5+h0adpqxPLCGrK7C51uMCYETtFUzTHAZ0tyzrbdM0a4esV6kOZNsAL/BLy7L2NK1kEREROVyThnQ85FGnljBo0JB9AQvglVdeYN68uQDs2PEdmzdvPihkde7chf79TQBMcyDbtm1tcD+bNm2kc+cu9OjRE4CJE89i+vRXmDz5QjyeFO6++w5GjjyJkSNHA3DEEUfxpz/dxmmnncEpp5wal/cKsYWsuu5Z3Dc8q2maDuAB4Mo61jseiABdgGzgY9M03997VKwuTqeB3++NoaymcTodLbIfSQz1N/mpx8lPPW66774zcDoTd4+b0+nAMKprcDodpKWl7atn8eKFfPnlFzz55H9ITU3jpz+9mkgktO95p7N6O4/Hs2+Zy+UiFArW+54Mw8DhcOBwGBgG+9bb+zglxcO//vUsCxd+wZw5s5k+/WUefvhxbrnldyxf/jWffvoJ//M/P+SZZ144KOztff3D+TsZS8gqALrXetwNqB0jM4ChwDzTNAE6ATNM0zwHuAR417KsELDDNM1PgWOBekNWJGJTVFQR8xtoLL/f2yL7kcRQf5Ofepz81OOms207oSOqp6SkUl5eTiQS3VfH3t9LSkrIyMjA7U5h/fr1LF/+9QHrVdde+z1EozbRaP3vybZtotEo3br1YOvWrWzcuJFu3brzzjszOeqo4ZSWllFVFWDEiJEMGjSEiy46n0gkypYtBQwcOISBA4fwyScfsW3bNtLTM+t8/QP/TubnZ9T7/mMJWV8C/U3T7A1sAaZSHZ4AsCyrGMjb+9g0zXnAjTV3F44FTjNN8zmqTxeeAPwthn2KiIhIG5eV5eeII47isssuJCUllZycnH3PjRgxkjffnM4VV0yle/eeDB48NG77TUlJ4Te/+QO33nrzvgvfzzvvB5SUlPDrX19PMBjEtm2uvfZ6AB555O8UFGzCtm2OOeZ4+vUbEJc6jFgmZjRN80yqw5ET+JdlWX8yTfN2YKFlWTMOWHce34esdODfwGCqTzv+27Ksew+1r1AoYutIljSV+pv81OPkpx433fbtG+nUqWeiy6hXW5m7cK+6fp75+RmLqD5Ld5CYQlZLau6QFY7anPfkF9x29mCO7Vz/IT5p2/TlnPzU4+SnHjedQlZ8HW7Iimkw0mRiALvKqli2pVghS0REpA26776/8PXXS/dbdsEFU5k06ZwEVVS3dheynA6DbK+H3WWHN5O2iIiItA433HBzokuISbucuzDX52FnWVWiyxAREZEk1k5DlptdClkiIiLSjNplyMrxethZqtOFIiIi0nzaZcjK9XnYXV5FtJXdWSkiIiLJo92GrFDEpiQQTnQpIiIiUuOMM6rnEty1aye/+91Nda5zzTX/y6pVK+p9jSlTzqaoqKhZ6jtc7TNked0A7KnQKUMREZHWJi8vnzvvvCfRZTRZuxvCAaqPZAHsLg/SJ9eX4GpERESaX8qqV0ld+WJcXzMwaCpVA6fU+/yjjz5Ip06dmTz5AgCeeuqfGIbB0qVLKC0tIRyOcPXVP2H06DH7bbdt21Zuuuk6nn32ZaqqAtx11x/59tsN9OzZm6qq2G9ce/HF55g5s3pimrPPPo8LL7yEyspKfv/7W9ixYwfRaIQrr7yKsWPH8dhjD/Hpp/NxOp0cd9wJXHPNdYf/AzlA+wxZ3r0hK5TgSkRERJLX6aeP48EH798XsubOfZ+//vUhLrroEny+dEpLi7nqqis46aRTMAyjztd4/fVXSUlJ5emnX2Tt2jX86EeXxrTvVatWMmvWWzz++NPYts3//u+VHH30cLZu3UJeXj733vt3AMrKyigpKWb+/LlMm/YahmFQWloal/ffPkNWrSNZIiIi7UHVwCmHPOrUHAYMGEhh4R527dpJYWEhGRkZ5OXl8eCD97F06RIcDgc7d+5kz57d5Obm1fkaS5cuYcqUqQD069efvn37xbTvZcu+4uSTTyUtLQ2AU045laVLv2LEiBN55JG/8+ijDzJq1GiOOmoY4XAYjyeFu+++g5EjT2LkyNFxef/t8pqs9BQnHpdDIUtERKSZjRkzlrlzP+DDD+cwduw43nvvHYqKinjqqed45pkXycnJIRg89L/H9R3lOrS6RxDo0aMnTz31LH379uMf/3iYf//7CVwuF0888TRjxpzG/PnzuOGGnzdifwdrlyHLMAzy0z268F1ERKSZjR07jg8+eI+5cz9gzJjTKSsrIzs7G5fLxaJFX7J9+7ZDbn/UUcN47713AFi/fi3r1q2Nab9HHTWcjz+eRyAQoLKykvnz53LUUUeza9dOUlJSGT/+TC6++DJWr15FRUUF5eVlnHjiSfziFzewZs3qJr9vaKenCwHy0lN0TZaIiEgz69OnLxUV5eTn55OXl8e4cRO5+eZf8qMfXcaAASY9e/Y65Pbnnz+Fu+76I1dcMZV+/QYwaNCQmPZrmgOZOPEsrr76cqD6wvcBAwby+ecLePTRv2MYDlwuFzfeeAsVFRX8+tfXEwwGsW2ba6+9vqlvGwDDbmUDcoZCEbuoqKLZ9/PrmavYuLucaZcf0+z7kpbn93tpib9HkjjqcfJTj5tu+/aNdOrUM9Fl1MvpdBCJRBNdRszq+nnm52csAo6ta/12eboQIDfdo2uyREREpNm029OF+ekpFFaECEdtXI7GXFAnIiIiiXL11VcQCu1/2c+tt94e892HLaHdhqy8jBRsoKgiSF56SqLLERERkcPwxBNPJ7qEBrXb04X5NcFqd4UufhcRkeTV2q69bqsa83NsxyFLA5KKiEhyc7k8lJeXKGg1kW3blJeX4HJ5Dmu7dnu6MHfvkSyFLBERSVLZ2fkUFu6krKwo0aXUyTCMNhMAXS4P2dn5h7dNM9XS6uXpSJaIiCQ5p9NFXl7nRJdRr2QfpqPdni70elz4PE5dkyUiIiLNot2GLKieKHqPjmSJiIhIM2jfIcvrZrfmLxQREZFm0L5Dlk+jvouIiEjzaNchK8fr0STRIiIi0izadcjK9XkorQpTFW47k1OKiIhI29DOQ5YbgD26LktERETirJ2HrOqxsnSHoYiIiMSbQhawS9dliYiISJy165CV460Z9V2nC0VERCTO2nnIqr4mS8M4iIiISLy165DldjrISnUpZImIiEjcteuQBTVT62j+QhEREYkzhSyN+i4iIiLNQCFLIUtERESaQbsPWTleN7vLg9i2nehSREREJIm0+5CV5/MQCEepCEUSXYqIiIgkkXYfsvYOSKqJokVERCSeFLK8mlpHRERE4s8Vy0qmaU4A/g44gScty7q7nvWmAK8Ax1mWtbBm2ZHAP4FMIFrzXCAOtcfFviNZGvVdRERE4qjBI1mmaTqBR4CJwGDgYtM0B9exXgZwLfB5rWUu4DngJ5ZlDQHGAK3qvFyOT6O+i4iISPzFcrrweGCtZVnrLcsKAi8C59ax3h3APUDto1TjgGWWZS0FsCxrt2VZreoK86xUN05DIUtERETiK5aQ1RXYXOtxQc2yfUzTHAZ0tyzr7QO2HQDYpmnONk1zsWmaNzWp2mbgdBhkez268F1ERETiKpZrsow6lu0bVMo0TQfwAHBlPa9/EnAcUAF8YJrmIsuyPqhvZ06ngd/vjaGspnE6Hfv20yEzlZJQpEX2Ky2jdn8lOanHyU89Tn7J3uNYQlYB0L3W427A1lqPM4ChwDzTNAE6ATNM0zynZtuPLMvaBWCa5ixgOFBvyIpEbIqKKg7nPTSK3+/dtx9/qpPtRZUtsl9pGbX7K8lJPU5+6nHyS4Ye5+dn1PtcLCHrS6C/aZq9gS3AVOCSvU9allUM5O19bJrmPOBGy7IWmqa5DrjJNE0vEAROofqoV6uS6/Wwdmd5ossQERGRJNLgNVmWZYWBa4DZwErgZcuylpumeXvN0apDbVsI3E91UPsKWGxZ1symlx1fOT4PeypCRDW1joiIiMRJTONkWZY1C5h1wLLf17PumAMeP0f1MA6tVq7PQzhqUxII409zJ7ocERERSQLtfsR3gFyvxsoSERGR+FLI4vtR3/do1HcRERGJE4UsNEm0iIiIxJ9CFt9PEq3ThSIiIhIvCllAeooTj9NQyBIREZG4UcgCDMMg1+dht67JEhERkThRyKqR6/PoSJaIiIjEjUJWjVxv9YCkIiIiIvGgkFVDR7JEREQknhSyauR43RRWhAhHNbWOiIiINJ1CVo1cnwcbKNLF7yIiIhIHClk1NCCpiIiIxJNCVo19IUtHskRERCQOFLJq5Po0SbSIiIjEj0JWDU2tIyIiIvGkkFUj1e3E53GyW2NliYiISBwoZNWisbJEREQkXhSyasn1uhWyREREJC4UsmrJ9XnYo7sLRUREJA4UsmqpPl2oa7JERESk6RSyasnxeiitClMVjia6FBEREWnjFLJq2TtWlk4ZioiISFMpZNXy/dQ6ClkiIiLSNApZtWj+QhEREYkXhaxa9o36rtOFIiIi0kQKWbXkeDV/oYiIiMSHQlYtLqeDrFSXQpaIiIg0mULWATS1joiIiMSDQtYBNCCpiIiIxINC1gE0tY6IiIjEg0LWAXK91acLbdtOdCkiIiLShilkHSDX5yYQjlIRiiS6FBEREWnDFLIOoAFJRUREJB4Usg6wb0BS3WEoIiIiTaCQdYC9R7J08buIiIg0hULWAXJ9GvVdREREmk4h6wBZaW6chkKWiIiINI1C1gEchkG2VwOSioiISNMoZNUh1+dht67JEhERkSZQyKpDrs+t04UiIiLSJApZddg76ruIiIhIY7liWck0zQnA3wEn8KRlWXfXs94U4BXgOMuyFtZa3gNYAdxmWdZfm1x1M6uevzBE1LZxGEaiyxEREZE2qMEjWaZpOoFHgInAYOBi0zQH17FeBnAt8HkdL/MA8E7TSm05OT4P4ahNSSCc6FJERESkjYrldOHxwFrLstZblhUEXgTOrWO9O4B7gEDthaZpngesB5Y3sdYWk+vVWFkiIiLSNLGErK7A5lqPC2qW7WOa5jCgu2VZbx+w3AfcDPyxiXXGlf+18zBWvV3v89/PX6iQJSIiIo0TyzVZdV2UZO/9g2maDqpPB15Zx3p/BB6wLKvMNM2YCnI6Dfx+b0zrNpbTYWPM+Q3+//sCXKkHPd87FAUgQPPXIs3D6XSod0lOPU5+6nHyS/YexxKyCoDutR53A7bWepwBDAXm1QSpTsAM0zTPAUYAU0zTvAfwA1HTNAOWZT1c384iEZuioorDexeHyX3cTfjfvIiqT/5J5dFXH/S8J1odsjbvLGv2WqR5+P1e9S7JqcfJTz1OfsnQ4/z8jHqfiyVkfQn0N02zN7AFmApcsvdJy7KKgby9j03TnAfcWHN34ehay28Dyg4VsFpKqNsoor3H4F30EIHBU7E9+/+AfB4nHqeh04UiIiLSaA1ek2VZVhi4BpgNrARetixruWmat9ccrWqTomNuxRHYQ9pXjx/0nGEYGvVdREREmiSmcbIsy5oFzDpg2e/rWXdMPctvO8zampXdZRhVfSeR9tXjVA69Atubt9/zuT4NSCoiIiKN165HfC8fcRNGOIB30UMHPZerSaJFRESkCdp1yIpk9yUw6ELSvnkWR0nBfs9Vj/quI1kiIiLSOO06ZAFUHPdLMAx8X96/3/Jcn5vCihDhqF3PliIiIiL1a/chK5rehcojriTFehXnbmvf8hyvBxso0tEsERERaYR2H7IAKo65Btvtw/f5PfuWfT/qu67LEhERkcOnkAXYqdlUDvsJKRtm49q+CPg+ZO3SkSwRERFpBIWsGhVHXkU0LQ/fZ3eDbZPrq54keo+GcRAREZFGUMjay+Oj/Nhf4NmyAPfm+eR6NUm0iIiINJ5CVi2BIT8kktEd32d3k+oy8Hmc7K7QNVkiIiJy+BSyanN6KB9xA+6dX5OydqZGfRcREZFGU8g6QFX/8wnnmHg/v4cOaZokWkRERBpHIetADiflJ9yCq3gD5zBXIUtEREQaRSGrDsFepxPqdCznljxPeUV5ossRERGRNkghqy6GQfmJt5AV3sXk8EyqwtFEVyQiIiJtjEJWPUJdTmBz9kh+6ppBUdGuRJcjIiIibYxC1iGsGfwL/EY56V/9I9GliIiISBujkHUIrs5H8mZkJN3WPoNRviPR5YiIiEgbopB1CLleD/eHp+CIhPAueyrR5YiIiEgbopB1CDleNxvtTmzIPJaUtW+DbSe6JBEREWkjFLIOweV0kJXqYlHaaJwlG3HtWp7okkRERKSNUMhqQK7Pw3zH8diGE8+6mYkuR0RERNoIhawG5Po8bAx4CXUdqVOGIiIiEjOFrAbk+jzsrghS1XcSruINOPesSnRJIiIi0gYoZDUg1+thT3mQQO9x2IaDlLU6ZSgiIiINU8hqQK7PTSAcpdydQ6jLCFLWzUp0SSIiItIGKGQ1INfnAWB3eaj6lGHhapx7Vie4KhEREWntFLIakOvdG7KCBPtMwMYgRXcZioiISAMUshrw/ZGsIFFfJ8Kdj1PIEhERkQYpZDUg1+cGYE9FEICqvmfi2r0KZ+G6RJYlIiIirZxCVgOy0tx4nAYb91QC1SEL0AXwIiIickgKWQ1wGAYje+fw4ZpdRKI20fQuhDodo9HfRURE5JAUsmIwfmAHdpUHWVxQBEBV30m4d32Do/jbxBYmIiIirZZCVgxO6pOD1+1k9qqdAFT10SlDEREROTSFrBikup2M6Z/Lh6t3EQxHiWZ2I9ThKN1lKCIiIvVSyIrRuIEdKK0Ks+DbQqDmlOGOpThKChJcmYiIiLRGClkxGtHDjz/NzexVO4Badxmu1ylDEREROZhCVoxcTgdjB+Qxf91uKoIRolm9COUN1SlDERERqZNC1mGYMLADVeEoH63bBUCw7yTc2xfhKN2a4MpERESktVHIOgxHds2kY0YKs1fW3GXYbxKgU4YiIiJyMIWsw+AwDMYPzOezjYUUVYSI+PsQzh2koRxERETkIApZh2ncwA5EojYfrKk5mtV3Eq5tX+Io357gykRERKQ1iSlkmaY5wTRNyzTNtaZp3nKI9aaYpmmbpnlszeMzTNNcZJrm1zW/nxavwhNlQL6P3jne7wcm7TsJAxvP+ncTXJmIiIi0Jg2GLNM0ncAjwERgMHCxaZqD61gvA7gW+LzW4l3A2ZZlHQFcATwbj6ITyTAMxg3MZ0lBMdtLAkRy+hPOHqC7DEVERGQ/sRzJOh5Ya1nWesuygsCLwLl1rHcHcA8Q2LvAsqwllmXtvfVuOZBqmmZKE2tOuPEDOwAwx9p7NOtM3Fs/x6jYmciyREREpBWJJWR1BTbXelxQs2wf0zSHAd0ty3r7EK/zA2CJZVlVh11lK9M9O40hnTK+P2XYbxKGHSVl/ewEVyYiIiKthSuGdYw6ltl7/2CapgN4ALiyvhcwTXMI8BdgXEM7czoN/H5vDGU1jdPpaNJ+zh3WlbveWcXuUJS+fYZj5/TDt+kdUk/63zhWKY3V1P5K66ceJz/1OPkle49jCVkFQPdaj7sBtUffzACGAvNM0wToBMwwTfMcy7IWmqbZDXgduNyyrHUN7SwSsSkqqoi1/kbz+71N2s/oHlkYwKtfbOLHo3rh7T0R7+JHKd5WgJ2WE79CpVGa2l9p/dTj5KceJ79k6HF+fka9z8VyuvBLoL9pmr1N0/QAU4EZe5+0LKvYsqw8y7J6WZbVC/gM2Buw/MBM4NeWZX3alDfR2uSlp3BMDz/vWTuxbZtg30kYdoSUDbrLUERERGIIWZZlhYFrgNnASuBly7KWm6Z5u2ma5zSw+TVAP+BW0zS/qvnVoclVtxITBuazqbCSld+VEc4bQiSzp+4yFBEREQAM27YbXqsFhUIRuy2cLgQoCYQY/9hnXDisC78c0xffgrtI++pxdv/PEuzU7DhVKo2RDIeg5dDU4+SnHie/ZOhxfn7GIuDYup7TiO9NkJnqZlTvHOZYO4lE7eqBSaNhPBvmJLo0ERERSTCFrCYaNzCfnWVBvtpSTDj/SCIZ3XTKUERERBSymurkvrmkuR28u3IHGAZVfSfh2Twfo6ok0aWJiIhIAilkNVGq28kp/fL4cM0uQpEoVX3PxIiG8Hz7XqJLExERkQRSyIqDCQM7UBIIs+DbQsIdhxHJ7Ilvwd04yrY2vLGIiIgkJYWsOBjR009WqovZK3eA4aB44hMYwTKy3roMo6o40eWJiIhIAihkxYHL6eB0M5/563ZTEYwQyRtMycQncRatJ/OdqyDS5qdrFBERkcOkkBUn4wbmEwhHmb9uNwCh7idRetp9eLYsIOOD68GOJrhCERERaUkKWXFydNcsOqR7mL1qx75lVeZkyk64hdQ1b+JbcFcCqxMREZGWppAVJw7DYNzADiz4tpCiytC+5ZXDf0bl0CvwLvkHqcv+lcAKRUREpCUpZMXRhIEdiERtPlyz6/uFhkHZ6Nup6j2e9I//gGfdrMQVKCIiIi1GISuOBnTw0Ssnrfouw9ocTkrOeJhwx2Fkzvk5rm1fJqZAERERaTEKWXFk1JwyXFJQzHelB9xR6E6jeNJ/iKR3IWvm/+AsXJeYIkVERKRFKGTF2fiBHbDh4KNZgJ2WQ/HZz4HDRdZbl2KUH7yOiIiIJAeFrDjrkZ3Gsd2zeHZhASWB0EHPR7N6UnzW0zgqd5E18wqMYFkCqhQREZHmppDVDH45pi8lgRCP/3djnc+HOxxFyfh/4Nq1gszZP4bIwWFMRERE2jaFrGYwoEM6k4/szKtfbWXtzvI61wn2GkvZmLvxbPqI9Hm3gG23cJUiIiLSnBSymslPRvUiPcXFX+euxa4nQAUGX0z5cb8kbdVL+P57p0aFFxERSSIKWc0kK83N/53Ui0Wbi3l/9a5616s47noqj7gC71f/JOO9ayAcaMEqRUREpLkoZDWj847ozIB8H3+bt47KUKTulQyDstF3Unbib0ldOwP/m1MxKve0bKEiIiISdwpZzcjpMPjVaf3YURbkP19srn9Fw6By+P9RPP4fuHZ+TfarZ+MsWt9yhYqIiEjcKWQ1s6O7ZTFhUAee+3IzBUWVh1w32O8sis57GSNUhv/Vc3Bv/byFqhQREZF4U8hqAT8f3Runw+Bv8xo+OhXudAyFP5hB1JtH1psXk2JNb4EKRUREJN4UslpAh4wUfnRCTz5at5sF3zZ8vVU0qydFk98g1PkYMt+/Fu+XD2iIBxERkTZGIauFXDy8Kz2y07jvw3WEIg0P1WCn+ik++3kC5hR8X9xHxofXQyTYApWKiIhIPChktRCPy8H1Y/qysbCSFxdviW0jp4fSsQ9QfvwNpK56hay3fogRKGreQkVERCQuFLJa0Kg+OZzUJ4cnF2xiV1lVbBsZBhXH/ZKS0/+Oe9si/NPPw1Fc93Q9IiIi0nooZLWw68f0JRSN8vDHGw5ruyrzBxSfOw1HxS6yXzuHFGs6rh3LqsfU0vVaIiIirY4r0QW0N92z0/jhMd34zxebOf/IzhzVNSvmbUNdTqBoygyy3rqMzPev3bfcdnmJZHQlmtGVSEb3mj93I5LRjWhGV6K+jmAoT4uIiLQkhawE+J8RPZi14jvum7uOf18yDKfDiHnbiL8Pey7+AFfhWhylBThLC2r9voWUHUtxBAr328Z2uAl1G0nJ2L9he/Pj/XZERESkDgpZCeD1OLn25D78btYqZnyznfOP7Hx4L+BKJZw/FPKH1v18sBzn3uBVtgVn0bekLX+G7JcnUjLhccKdhjf9Te6ggcEAACAASURBVIiIiMgh6RxSgowbmM+wrpk8+sm3lARC8X1xj49Irkmw11gCQy+n/KTfU/iDGeD04H99CqkrpsV3fyIiInIQhawEMQyDG07rR0kgxOP/bf67BSN5gym8YCahrieSMfcm0ufeDJEY73AUERGRw6aQlUBmh3TOP7Izr361lbU7y5t9f3ZqNsVnPUPF8GtIW/E8/tcvwFG2rcmvawQKNaG1iIjIARSyEuwno3qRnuLid7NWUlQR59OGdXE4KT/xFoon/BPnHovsl89s3ETUto1r+yIy3r+O3P8cS/YLY3HuWhH/ekVERNoohawE86e5ueusQRQUBfjZq8sormyBoAUE+06iaMpbRD3pZL15EanL/h3beFuhClKXP4//5Ylkv3YunvXvEBh4IXaKn4wPfqmpf0RERGooZLUCx/fM5t5zB7NhTwU/f+1rSgPhFtlvJGcARRfMJNjjVDI+vrU6JIUr61zXuWc16fN/R+5/jiFj3s0YdoTSU/7MnisXUTbmz5SOuRv3ruV4Fz3UIrWLiIi0dgpZrcSJvXK495whrNlZzs9f+5qyqpYJWnZKJiVnPlU9P6L1Kv7pk3GUFFQ/GQmSsmYGWa9PIeeF00hdPo1gr9MpnPw6hRe9R2DoZdiedACCfcYTGDAZ76KHcO38ukVqFxERac0Mu5VNyRIKReyioopm34/f76Ul9nO4Plq7m5vfWsHgjhk8NGUoPk/LDWXm+fZ9Mub8HBwuAgMmk7pmBo7KnUQye1A55FICgy7CTsutd3sjUET2C2OxU/0UXjgLnCktVvuBWmt/JX7U4+SnHie/ZOhxfn7GIuDYup7TkaxW5pR+udx11iBWbC/huunfUBGMtNi+g71Op+iCmUTT8kn7+t+EOh5N8VnPsOfST6gc/tNDBiwAO9VP2an34Npj4f3yby1UtYiISOukkNUKndY/jzsmDWLZ1hKuf+MbAqGWC1oRfx8Kp77H7v+3jJJJ/ybY87TDmvcw2GsslQMvwrv4EVzffdWMlYqIiLRuMf3raZrmBNM0LdM015qmecsh1ptimqZtmuaxtZb9umY7yzTN8fEouj04w8znjxMHsqSgmOvfWN6iQQuHCzvV3+jNy0/6A1Ffx5oL6QNxLExERKTtaDBkmabpBB4BJgKDgYtN0xxcx3oZwLXA57WWDQamAkOACcCjNa8nMZgwqAO/H2+ycFMRv3pzBVXhaKJLiomdkknpqX/FVbgG3xd/TXQ5IiIiCRHLkazjgbWWZa23LCsIvAicW8d6dwD3ALUPXZwLvGhZVpVlWRuAtTWvJzGaNKQjvxs3gM82FnLzjBUE20jQCvU4hcrBPyRtyT9xbVuY6HJERERaXCy3rnUFNtd6XACMqL2CaZrDgO6WZb1tmuaNB2z72QHbdj3UzpxOA7/fG0NZTeN0OlpkP/Fw+eg+eFLd3DpjOb+fvZoHLzoaj6sNXE535l2wZT7+eTcQvuojcLfcz7st9VcaRz1Ofupx8kv2HscSsow6lu0b98E0TQfwAHDl4W5bl0jEbpHbOdvabaMT+udSOrYf93ywlp89v4g/nzUIl7O1By0n7jF/xf/mRYRm30b5Sbe12J7bWn/l8KnHyU89Tn7J0OP8/Ix6n4vlX+kCoHutx92ArbUeZwBDgXmmaX4LnADMqLn4vaFt5TBccHQXbji1L/PW7ubO91YnupyYhLqNovKIK0hb+hTurZ81vIGIiEiSiOVI1pdAf9M0ewNbqL6Q/ZK9T1qWVQzk7X1smuY84EbLshaaplkJTDNN836gC9Af+CJ+5bc/U4d3pbgyxJOfbeKEXjlMGNQh0SU1qOzE3+LZOI+MD25gz0XvgceX6JKSjqN8O1FPFrjTEl2KiIjUaPBIlmVZYeAaYDawEnjZsqzlpmnebprmOQ1suxx4GVgBvAv8zLKsFhyLIDn96MSeHNklk798sIbtJW1giAS3l9Kx9+Mo2UT6Z3clupr4iLTMRN4Nsm1Sv36anGdH4Z9+HkZVcaIrEhGRGppWp40qKKrkh88sZmDHdB694Eicjrouf2tdfJ/chnfpkxSd+xKhbqOadV/N0V9n4To8G94lZf27uHYsJTBoKmUn35Gw6YOMQBEZc39Fyvp3CHU+Htd3Swh3HEbR2c+3iyNabf0zLA1Tj5NfMvRY0+okoW7+NG48rS+LC4qZtqgg0eXEpHzEzYT9fcj48AaMYGmiy2mYbePasRTvZ38he9qp5Ew7hfQFfwY7SlX/80hbMQ3/6xfgKN/e4qW5ti8i++UJeL6dQ9nI31F0/quUnv4grm1fkvnu/0Ik2OI1iYjI/lpu9mGJu7OGdOST9Xt49JNvOb5nNmaH9ESXdGjuNErHPoB/+vn4Pr6NsjF3g9Od6Kr2Fw3j3vo5nvXvkrJhNs6yrdiGk1CXEZQOvYxg7wlEM7oAUNVnApnvX4f/5TMpmfA44c51/kcmvuwoaUsew/fZPUTTu1B0/nTCnYZX19P/bIxgCRnzbibj/esoPeMhcGjsXxGRRFHIasMMw+DXZ/Rn2dYSbp21imd+OIxUd+v+RzXc6Rgqh/0E7+JH8WyaS2DQVAKDLyaa2b3hjRvLtiFShRGpgnD179//OYARCWJU7iZl44d4vp2DI1CI7Uwh2P0Uyo+/kWDvM7BTsw962WDfMyn09yFr1o/wv3EBZSffQWDIpc32NoyKXWR+8As8mz4i0Pcsyk79C3ZK1n7rBIb8EKOqhPQFf8L2ZFQHWaP1n0pubRwlm/F+9U/Kj7+xSVNMiUj7pmuyksBn3+7h5699w0XDunDjaf0SXU7D7Ciebz8gdflzeDZ+CNSMED/kUoK9TgdH47K/o2wbnk0f4dk0F893i7GDZRjhKoxobKfOop5Mgr3GUtVnAsHuY2K+C9IIFJE55xo8m+ZROfiHlJ18e9yv03Jv/oSM96/FUVVM2Ul/JDDkh4cMT74Fd+Nd/DAVw39K+Ym/iWstrUVzfoYz3ruG1DVvUNX3TErG/1NBNUGS6Xta6pYMPT7UNVk6kpUETuiVw0XDuvDSkq2M6pPDib1yEl3SoRkOgr3PINj7DBylW0hd8QKpK18g652riHg7Ehg8lcCgi4lmdjv060SCuLcvxLNxLp5N83DtXlm92NcRu/cpBJyZ4PRgO1OwnSngSq315xRspwdqHttuL+G8IeD0HPbbsVP9FE96Gt/n9+Jd/DCuPasomfA4UV/Hxvx09hcN4/3ifryLHiKS3ZfCc54nkjuowc3KT7gZI1iCd/GjRFOyqBz+s6bX0k44SjaRsnYG4azepKybRery5wgMvSzRZYlIG6QjWUkiEIpw+fNLKAmEefHyY/B7W9m1Tg2Jhr8/urVpHgDBnqcSGHIpwZ6n7Tu65SjdUn2kauNc3AWf4giVYTtchDofR7DHqQR7jCGSOwh/ti8h/fWsfZvMD35J1JNJycTHCXc6ptGv5SjdSuaca3Bv+4LKgRdV38l4OFMT2VEy5vyc1DVvUnrK3QSGNt+pzERors9w+vzfkrp8Gnsu+5SMub/CveUzCi94O6ZwK/GVbN/TcrBk6PGhjmQpZCURa0cZVz6/hNF9c/nL2YMw2ugpDkdJAakrXyB1xYs4K74j4utEsMcpuLcvwVVYPdJ9JL1LdajqeSqhbqOwPftPa5DI/jp3ryRr1lU4yrZSdvKd1af2YmVHcZQW4N76Bemf3AbREGWn/Jkqc3LjiomEyHznKjwbP6R03MNU9a9rbve2qTl6bFTsIveZEQQGnEfZafdhVOwi+6Vx2ClZFF4ws0Xn35Tk/J6W/SVDjxWy6pAMja3Ls19u5sH5G7h1/ADOGdop0eU0TSSEZ+P7pC1/HvfWLwh1OoZgz5qjVdn9D3mdTKL7awQKyXzvGjybP6JyyKWUjb59/1OR4UqcRRtwFa7FWbgGZ+G66j8Xrau+KB8I5Q2ldPyjRPx9mlZMuJKsty7FvX0RJROfIthrbNNer5Vojh57P78X78IHKbxkLpHs6usb3Zs/JmvGJQQGT6Xs1Hvjuj85tER/jqX5JUOPFbLqkAyNrUskavOzV5excnsZz18+nG7+5B+Usi6tor/RCL7P/4J38aOEOh1LqONwnIVrcBWtw1GyGaNmrnQbg2hmD8LZfYlk9yeS3Zewvx/hjsPiNsSFUVVC1psX4dqzmuJzphHqMiIur5tI8e6xESwj55kRhLqOpGTiE/s9t/dGgpJxj1LV/5ATXUgctYrPsTSrZOixQlYdkqGx9dleEuDiZxbRO8fH41OPwtUGRoOPt9bU35Q1M0if+yuMaJiIvy/hnP5E/H2JZPcjnN2PiL83uJo/DBuVu/FPn4yjYgfF571COH9os++zOcW7x2lfPU76p7dTOOWt6oBbWySE/40pOPespvDCd4lm9YzbfqV+relzLM0jGXqskFWHZGjsocxeuYPfzVrFj0f25KoT298/CK2uv5Fg9cX7RmInWXCUbsU//XyMSICis6cRyR+S0HqaIq49jgTJeXYkEX8fis97uc5VHCWbyX55ApGs3hRNnt6oO1Hl8LS6z7HEXTL0WNPqtEPjB3Vg/MB8nlywkeXbShJdjjg9CQ9YANGMLhSf+wJgkPPyePwvTSBt8SM4ijcmurSESln9Os7y7VQM/2m960Qzu1N66j24d3yF7/N7WrA6EWmrEv+tL83m5rH9yUtP4ffvWFSGIokuR1qJiL8Pey56j7JRfwCnm/QFfyb3uVH4X5lE2pJ/4CjdkugSW5Ydrb5uLm8Ioe6nHHLVYN9JVA65DO+Sf+DeOLeFChSRtkohK4llpLr440STzYWV3Dl7NQEFLalh+zpQefTVFE15i92XLaDsxN8CkP7fO8l9ZgT+V88h7asncJRtTXClzc+zYTauonVUDv9pTCO7l530e8I5JpkfXIej/LsWqFBE2iqFrCR3THc/PxnVi/esnVz67GK+0alDOUA0szuVw/+PogtmsvvSTyg74RaIVJH+6R/Jffp4/NPPJ3XZv5LzCJdt4130CJHMnlT1nRTbNq40SsY/hhEqJ+P9X4Adbd4aRaTNct52222JrmE/0ah9WyAQavb9pKa6aYn9tAbDumVxVNdMPlyzixcWbyEYjnJ01yycSXzXYXvqbzzZqX7CXY4nMPQyqgacR9Sbh2vn16StegXv0idJtV7DuXsFRrAMO9V/0CCwLSkePXZvXYBv0UOUn3DzwXcUHoKdlkvUm4932VPYzhTCXY5vUh1SN32Ok18y9NjnS9kGPF7Xc7q7sB0pqwpz/9x1vLX8O/rl+bhtoonZIT3RZTWL9tjf5uTcswbP5o9wb1mAe+tnOKqKAYhk9iTY9URCNb+i6V1arKZ49DjrrUtx7fyG3ZcvOPxhNGybjPd+Rsq6mRSd/xrhznXeXCRNoM9x8kuGHmsIhzokQ2Mb6+N1u/nTnDUUVYa46oQeXHl8d1zO5Dpz3J772+yiEZy7V+HZuuCg0BXO6lUduLqcSLD7ydjevGYro6k9du5cTs7L4yk74RYqj7mmUa9hVJWQ/fIEiEYovGg2dqq/0fXIwfQ5Tn7J0GOFrDokQ2OboqgyxF8/XMvsVTsZ1DGd2yaa9Mn1JbqsuGnv/W1R9YQu25VKxdE/pmLYT8ET579bwXL8rlKKHI2fOirjvZ/h+fYD9lzxOXZKVqNfx/XdEvzTzyfY6wxKJjwe08XzycBZtJ5oajZ2anaz7aPZP8eRUNxmVZDGSYbvao2TJQfxp7m5c9Ig7j57ENtKqrjs2cU8++VmItHWFbqlDXA4ieQPofKoqyg58yl2/79lFF4wi6re4/Et/Du5z51E6vLnIdr0u1uNQCHeL+4n95njcf3jeLwLH2rUheeO4o2krH2LwNBLmxSwAMIdh1F+wi2krH+HzHevxggUNun1Wr1IFd7P/kL2tDFkv3IWjrJtia6oUdI/vJGcaWMgWJ7oUiSJ6cL3dq5Pro9JgzuysbCSl5Zs5YtNRQzvlkVWWtv+3536m0CGg6ivI8G+kwj2GIN75zLSvnmalPXvEMnsQdTf+/BfsnwHvoUPkDHnWlI2zyfY/RQcnYeQ8tVTuHYuJ9hzDLhSY3493+f34tq1gtLxj2B7mn5dYrjTMdieDNK+eYbU1dMJdziSaEa3Jr9ua+PctQL/25eTum4mVf3OwrV7JSnrZhLse2Zcfo4Haq7PcYr1Kulf/LX6NLfblxRzebZVyfBdrQvf65AMhyjjybZt3lm5g3s/XEs4YnPdmD5MPrIzRhs99aH+tiK2jWf9LNL/exfOko0Eu59C2ajfEckd1OCmjpICvF89RuqKFyEaoqrf2VQccw2R3EH4s9Komv8wvv/eQTS9K8UTnyCSN7jB1zQqdpH7zAgC5mTKTr03Hu9wH9eOpWS89zOcJZuoOO6XVBxzLTiccd1HQkTDeBc/ivfLB4imZlN26j0Ee52Oa/sismZcQtTXkeLzXiHq6xjX3TbH59hRtIHslycQzh+K7U7Hve1L9lz2abOe9pT6JcN39aFOF+pIlgBgGAb989M5c1BHVu8s58XFWykJhBnRMxtHGwxa6m8rYhhEcgZQOfQy7NRsUta8Qdqyp3CUbiHc4cg6j4A4C9eR/t87yZj3K1w7vyFgTqZ03CMEhlyK7c0HIDXNQ5n/CILdRpGy5g28X/+HSEaXBoOWd9HDuLcsoHTcI3H/hzXq60TVwAtxlm7Bu+xfuLd+Rqj7SQkd6qKpnIVryZp5Bamrp1PV72xKznqaSF71nJfR9C6EupyA9+tn8GyYTVWfM+N6/V3cP8eRIFkzr8ARLKH4nBcIdTmOtKVPYRgQ6j46fvuRmCXDd/WhjmQpZMl+fCkuJgzqQEUwwouLt7JmZzkn981tc3cfqr+tkMNJuNNwAoMvwYiESVsxjbSvn4ZomFCHo8DpxrlrBekf/570j35TPQr7kEspHfcYVQN/cFAg2tvjaEZXAgMm496+CO/SJ3FU7iLY/aTqCbkPYARLyZzzc0K9TiMw9PLmeZ9OD8G+E4lk9iBtxTRSV0wjkt2PSHbf5tlfc7GjpC19kszZP8EIV1B62gNUHP/Lg4a6iGZ0JdTleNK+foaUDe9VD+rq9salhHh/jn2f3U3qurcpOeMhwp2GY3vzcRZvIHXVKwQGXdQspzzl0JLhu1ohqw7J0NjmYhgGJ/bKITPVxYuLt/D5xiJO7pdLmrvtnPZQf1sxVxqhHqcQGHAeztIt1dcxrXwZT8F80v97J47SrQSOvpqScY8R7DcJO6Xuo0D79djto8qcjBEJ4F32FJ7NHxPsPuagbdOW/ZuUb+dQOvYBoumNvzMxFpG8wVT1m4R788d4lz6BESgk1HVkneGvtXEUbyTznatIW/E8wZ5jKT7rWcKdhte7fjSjG6HOx1Zfe/ftHKr6nRWXoBXPz7F783wyPvo1lUMupXL4/+1bHs4bTNqyf2GEKgn2GhuXfUnskuG7Wtdk1SEZzgO3hHlrdvG7WavI83n4++Sh9MyJz/9Qm5v623a4ti0k/b934izeSOWRV1J5xJUx3fFXX48962aS8cEN4PRQMu5RQt1Pqn4iUkXOsyOJ+PtRfN5L8X4b9YtU4VvwZ7xLnyScO5iS8Y8Sye4X06ZGVTGu3Stx7lqJs3gDoc7HE+x9BjhTmqdW2yZ1+XOkf3oHtsNJ2ejbqTKnxDwshbvgU7JmXkEkqxdF576MnZbTpHLi9Tk2KnaR/dI47BQ/hRfMBPf+R+PS5/2a1JUvsOeSj4hm9Wzy/iR2yfBdrXGy6pAMjW0p32wr4frXlxO1be47bwhHdW3aLe8tQf1NfofqsbNwHZnvXI2zaC3lI26icvhPSV35Ehlzf0XROdMIdT+5hasFz7cfkPHBLzHClZSNvoPAoIu+Dy/RCM7iDbh2rcS5ewWu3Sur/1z2/XyRtjMFI1JFNDWHgDmFwOCpRHIGxK0+R+lWMub+Cs/mjwh2G03pafcRzTj8Efzdmz+pDlr+PhSd93KTrnuLy+fYtsmceQWegk8pvODtOm+4cJRvJ+fZUVT1O4vS0//etP3JYUmG72qFrDokQ2NbUkFRJb+Y/g3bSwL8ceJATjfzE13SIam/ya/BHgfLyZh3E6lr3qSq1zicRWuxXV6KLnwnYQOGOsq3kzHnF3i2fEpVnwlEU7Jw7V6Fa/cqjEgVALbDRcTfl3DuIMJ5gwnnDiKSN4hoWj7ugo9JW/ECng3vYURDhDodS+Xgi6nqd/bhn56LBHF/txj35o/xFHyC67uvwOmhbNStBIZc1qSfkXvzfLJm/g/h7H4Un/vi4QetaBjXjmWkd+hMkaNzo+sASFv6JOmf3EbpyXcSOOLKetfz/fdO0pb8k8Kpc4jkDmzSPiV2yfBdrZBVh2RobEsrqghxw5vLWba1hF+c0ocfHtO11Q7xoP4mv5h6bNukLfsXvv/egRENUzLuMar6n90yBdYnGqkZDuF+7JRMwrmD94WpcO4gIjn9GjwdaFTsItV6jdSVL+AqXEvUnU5V/3MJDL6YcIej6g5IdhTnrpV4Cj7BU/Ax7q2fY4QrsQ0H4Q5HEex2EoFBU+N2usy9aR5Zs35EOGcAxee8cOgph2wbZ/GG6sC3eT7uLQtwBEuwDScVw39GxXG/aNQpUtfOb/C/eg7BHmMoOfOpQwZHI1BIzrMjCXUdWb2utIhk+K5WyKpDMjQ2EQKhCLe9a/HB6l1ceHQXrj+1L05H6wta6m/yO5weu7YvwrNpHhXHXtd6xq2Kx5Quto1r+0LSVrxAytoZGOEA4dzB1Ue3BpyPESzBs/lj3AWf4in4BEdgDwDh7P6Euo0i2G00oa4nNHnU+/p4Nn5I5qyrCOcOpPjcF/bbjxEoxLP5E9wF8/Fs/hhnaQEAkYxuBLuPJtRtNOnbP8ax7AXC2f0pPe2vhDsdE/vOQxVkvzwRI1RG4UVzYro+zLvw7/g+v5fCH8w45IX+Ej/J8F2tkFWHZGhsokRtm4fmb+C5hQWc3DeXP00aSGoru/NQ/U1+6vH+jKoSUtbMIHXFNNw7l2EbDoyaKYcivo6Euo0m2O0kQt1GEU1v2im4w+H59n0y37macN5gykfchGfLAtyb5+Pa+TUGNlFPBqGuIwl2P5lQ99FEsnrvO+Lk93spXzaLjLk34SjbRuVRV1E+4qaDLlyvS/rcX5G64kWKz32RULdRMdVqBMvIeW4U4ZyBLXtzRDuWDJ9jhaw6JENjE+3lJVu4b+46BnXM4P7zh5Dj9SS6pH3U3+SnHtfPuXM5qWtnEPF1ItTtpOq7GRN4at+z4T0y3/0xRjSE7XAR7jicYPfRBLufXH16s55hLfb22AiW4VvwZ9K+eZpIZk9KT73nkMHJs/Ztsmb/hPJjfk7FCTcfVq17r+EqOufF7+9MlQYZgSJ8C/5EVd9JhHqMiXm7ZPgcK2TVIRka2xp8tHY3v525kjyfh4d+cATdsxv+H2ZLUH+Tn3rctjh3rcBZtpVQlxNiHvTzwB67tywgfe6vcBV/S+WQyygf+ZuDRtN3lBSQ/dI4Itl9KTp/+uGfko1UkfPcaKLeDhRNeavR4dQIFOLe+hnBnqc135AbrYSjdAtZb12Gq3A1tsND8Vn/ifkO3mT4HB8qZLWtYbyl1TmlXy7/vPBIyqrCXPXiV1g7yhJdkoi0QpG8wQR7nd6kUdVDXU+k8KI5VBz9Y1JXPE/2C2PxbPzw+xWiYTLnXAPYlIx7pHHXvDlTqDjuetw7vsKz4d1G1en6bgnZL00g652ryXl2JGlfPQGhth0k6uPcvRL/a+fiKN9G8YTHiWT3IWvW/8O99bNEl9YqKGRJkw3pnMmTU4/G7XTw45eWsrigKNEliUiycqdRPupWiia/ge1OJ+vty8l4/zqMQCHeL/+Ge/tCyk75M9HMHo3eRWDgFMLZ/fB9di9EI7FvaNukfv00/umTwTAoOe0+Iv4+pH/6R3KfOQHvwgcxqkoaXVdr497yX/zTfwDYFJ3/GsG+Z1J0zotEMrqR+fYVuLYvTnSJCaeQJXHRK9fLk1OPokN6Cte+9g0frd2d6JJEJImFOw2n8KJ3KD/2F6SseYOcaWPwLnqQwMALqRpwXtNe3OGi/PgbcRWuJmX167FtEywnY841ZMz/LcHuJ1N44TtUDbqI4vNeoXDyG4Q6Ho3v83vIeeYEvJ/dg1G5p2k1Jphn7dtkzbiUqK8jRT+YsW9idtubVz1kR1oeWW9fhmvnNwmuNLF0TZbEVVFliOumf8Oq70r57bgBnD20eeeHq4/6m/zU4+QXa4+du1aQ8eGNGJEqCn8wAzy+pu/cjuJ/ZRKOQBF7fvgROOu/sce5ZzWZ7/4YZ9G6fTMMYBx8DMO18xu8ix7Cs24WuFKr51Ec9mOivsR8TzZW2tKn8H1yG+HOx1J85r/qHGzWUVKA//XJGOFKis57lUiuWedrJcPnWBe+1yEZGttaVQQj3DRjOZ9vLOLak3tz2XHdW7wG9Tf5qcfJ77B6bNtgR+M6Dpp70zz8b116yNHiU1a/Tsbcm7DdPkrGPRLTcBHOPWvwLn6k+iiZ4SQw6EIqhv/04FOc0QhG5W6cFTtwlH+Ho+I7HOU7cOx9HCjEdqZge3zY7jp+HbA8mpJFJMds/M/IjlbPw7nkMap6j6dk3MPgqv9mJ0fRBvyvT8GwoxRNfo2Iv89B6yTD51ghqw7J0NjWLBiO8od3LN5fvZPLju3Gz0/u3aKjw6u/yU89Tn4J77Ftk/XGBbgK17H7sk/3n7ooUkX6J7eT9s3TBDuPoHT8I4d9RMpRsgnv4sdIXfkS2BGCPceCHakJUTtwVO7cN9ZZbdHUbKLeDkTTcjAiQYxQOUawvPr3UPm+KZrqEvF2oKr/OVT1P6/+2QHq3DBIxoc3kLr6dSqHXk7Z6DtiCmvOPWvwvzEF2+mh6PzXDgqSCe9xHChk1SEZGtvaRaI29364ltf+f3t3Hl9leed9/HOfJftykpA9WfDdaQAAIABJREFUgZAQLiBsYVVBpbRS3JeKImqxtrWO2pnWsdOZTudpp53O0+nmOE+tHfddXKtYtVbrQt2QJexwsUQgKwmQjSWQ5Tx/nEApBQwmJye5832/XrzOue9zn5xfXj/Pydf7vs51ra7h4pJMvjdnJL4+mh1e/XU/9dj9+kOPfTXLSXnhMvad8c8cnHwbAJ7mCpJevxl/3WoOlN7M/jP++aRzfXWHZ18NsavuI7r8NTqjk+iMzwyFqKO3GaHbuEw649M/fUqIjjac9gN/E748+3cRXf4qUdvfwuk8THtyAYeKL+PQyMvpSCk66Y9zDu8j6Q83EVWxhP3Tv8uBybed1tQW3t0bCLw4j2B0Mo2XP0dnwl8WHu8PPe6pHocsY8xc4C7AC9xvrf3pcY/fDNwKdAD7gJustRuMMX7gfmAS4AMetdb+31O9lkKWuwSDQe77cAf3fbiTc4vS+I8+mh1e/XU/9dj9+kuPk165AX/NMvZe9z7+2hUkvvn3EAzS8vlfcbhwbqTLO23OoSait71G9JYX8Ve+j0OQtvRxocBVfMlfrQjg7K8j+fdfxrdnIy2f+zmHRl/1mV7Tt2sVyYuvoTMuncbLnycYlw70nx73RI9CljHGC2wGzgMqgWXANdbaDccck2Stbe66fwlwi7V2rjFmAXCJtXa+MSYO2ADMstZuP9nrKWS50zNlVfz8rW1Mykvml5eVkBD92f+vrzvUX/dTj92vv/TYu3sDqU/PoS19HP76tbQNKaF57v/SmVwQ6dJ6zLN/F9FbXyZ68+/w160miENb7hmhy4lDxpD0x1vxHNhN89zfhiZW7QFfzTICixfQkTSUxsueJRib2m963BM9nYx0GrDVWlturT0MLAIuPfaAIwGrSzxwJLkFgXhjjA+IBQ4D7pkkRLrtqtJc/uOCUayubuYbT69mz/7DkS5JRKRbOoaMobX4Uvz1azk4ZgGNX3rRFQELoDM+k4MTvkbjvFfYe+0SDky7Hc/+XSS+811SnrsYp20/jZc90+OABdCePZWmCx/G27Sd5JevxTnU1Au/Qf/WndMJuUDFMduVwPTjDzLG3ArcDkQBR7rxHKFAVgPEAd+21g7syUHkM/vi6AwSY3x8d/EGvvJkGV8/cxhzR2fg92q6NhHp31pm/YyD426gPXtqpEsJm45AIQemfpsDU76Fb/c6/JXvc6jwfDqTh/Xaa7TlzaDp/PtJfvVGkl++Hi78Jd79baEHHQfwdN06oWkwusZ+BY9ue8HjIej4QgPvHS94fAQdb2jb4zvh9BmR0p3LhfOAL1prv9a1fT0wzVr7zZMcv6Dr+IXGmBnALcANQArwZ+B8a235yV6vs7Mz2NER/sH4Xq+Hjo6//daGhN+qikb+7aX1bNrVQnZyDDfOKOCqyXnERfXeJUT11/3UY/dTj93Lsa/gff4GnOBpzKh/GoIeH3h8dJ79XTrP+oewvMYRfr/3pJcLu/NXrRI4dqKjPKD6FMcvAu7pur8A+IO1tg2oM8a831XISUNWR0ewT67PuuE68EBVkBjFo9dO5INPGnjk45385NVN/PqtrVw9KZerJuaQHPsZ1hs7jvrrfuqx+6nHLpb5ObxXv05Sew3797f+ZZ4zwAkGgc7QPo7sDx49xulsD+3rbA+FtM4OCLbjdHZCsB26Hnc62zmUNoX2MP83lJ6eeNLHuhOylgHFxpjhQBUwn1B4OsoYU2yt3dK1eSFw5P5OYLYx5nFClwvPAP77tKoXV3IchxmFqcwoTGV1VRMPf1zBvR/s4LFlFVw+PpsFk/PITHT3yvUiIoNZR9oogoFJHHZxkP7UkGWtbTfG3Aa8TmgKhwetteuNMT8ClltrFwO3GWO+ALQBDcDCrqffDTwErAMc4CFr7Zow/B4ygE3ITebOy5PZWr+fR5dV8PTKKp4pq+aCMRlcPzWfgtS4T/8hIiIi/YwmI5V+p7qplSeWV/LSuloOt3cyq3gIN0zLZ0zWyU/JHk/9dT/12P3UY/dzQ497OoWDSJ/KSY7hO58fweKvT+Mr0/NZvrORhU+Ucde75bRpEKyIiAwQClnSb6XGRfF3M4ez+OvTuHJCNo8vr+SrT62iouFgpEsTERH5VApZ0u8lRPv47heK+dklY6hqauW6x1by2sZdkS5LRETklBSyZMD4XPEQnrh+EiMz4vk/r1p++NomDhwOzxwrIiIiPaWQJQNKVlIM91w1ga+dMZTXNtZx/eMr2bSrJdJliYiI/A2FLBlwfB6Hb8wo4DfzxtPa1sGNT63iyRWV9LdvyoqIyOCmkCUD1uT8AE98eTJnFqRy5zvl3P7iehoOaOFpERHpHxSyZEALxPr5xaVj+M7sIpbuaGDBoytZtrMh0mWJiIgoZMnA5zgOV5Xm8vCCUuKjvNz67Fp++cZmDYoXEZGIUsgS1xiZkcBj10/ikrFZ/HZJORfe+xF3vrONykbNqyUiIn1Py+qIK33SfIj7l5Tzpy276ewMcnZRGleX5jB1aADHcSJdnvQCvYfdTz12Pzf0+FTL6nzqAtEiA1Hp0BR+ctFo/qHlEM+vruaFNbUs2baHwrQ4rp6UywWjM4jxeyNdpoiIuJjOZIkrHd/fQ+2d/HFTHYtWVrG5fj9JMT4uG5fFlRNzyE6KiWCl8lnpPex+6rH7uaHHOpMlg160z8PFY7O4qCSTVVXNPF1WxePLK3l8eSXnjhjC1aU5TMhNxufRpUQREekdClkyqDiOQ2leMqV5ydQ2t/LsqhpeWlvD21t2E+f3MiY7kfE5SYzPTmJsdiLJsf5IlywiIgOUQpYMWllJMXzznOF8/cyhLNm2h1VVzaytbuaRpTvp6LqKXpAay7jsJMbnJDEuJ4nhaXF4NHBeRES6QSFLBr0Yv5c5ozKYMyoDgINtHWyobWFNdTNrqptZsm0PL6/fBUBCtJexXaFr+rAUSrIS8eoSo4iInIBClshxYv1eJucHmJwfACAYDLKz4SBra0Kha211C/d9sIN7P9hBcoyPMwpSmFmYxhkFKQR0eVFERLooZIl8CsdxGJYax7DUOC4qyQKgubWNj7Y38MEne/ngkwZe31SPx4GSrCRmFqYyY3gqIzPiNSeXiMggppAl8hkkxfiPXmLsDAbZWNvC+5/s5b3yvdzz/nbueX87Q+KjOGt4CjMK05g2NEBCtN5uIiKDiT71RXrI4ziUZCdRkp3ETWcVsGf/YT7cvpf3yxt4a8tuFq/bhc/j8JXp+XztzGEaOC8iMkgoZIn0srT4KC4qyeKikizaOzpZU9PMC6truO/DnWyp388PzzfER+mtJyLidlogWiSMfF4Pk/IC/PiCUdz+uSKWbNvDV59apUWrRUQGAYUskT7gOA7XTMrlf64YR/2+w9zwRBnLdzZGuiwREQkjhSyRPjS9IIWHF5SSGhfFbc+t4Zmyavrb+qEiItI7FLJE+lh+SiwPLpjImcNT+flbW/nPN7bQ1tEZ6bJERKSXKWSJREBCtI9fXFrCDdPyeXFtLbc8u4a9Bw5HuiwREelFClkiEeL1ONx69nB+cuEoNu7ax8LHy7C79kW6LBER6SUKWSIRNmdUBvfPn0BnMMhXF63iDVsf6ZJERKQXKGSJ9AOjMhN59LpJmIwEvvf7jdzz3id0akC8iMiAppAl0k+kxUdxz7zxXDo2iweXVvC1p1aztro50mWJiMhnpJAl0o9E+Tz865xifjB3JNXNrdz41Cq+9/uNVDVp8lIRkYFGa3uI9DOO43BRSRazi9N5dFkFjy+v5J2tu5lfmstXpg8lMUZvWxGRgUBnskT6qbgoLzfPKOD5G6fyxVEZPL68kssf+Jhnyqpo17xaIiL9nkKWSD+XmRjND+YaHrtuEsUZCfz8rW3Mf2QF727do9niRUT6MYUskQHCZCbwmyvH8cvLSgC446X1/N2za9i0qyXClYmIyIlocIfIAOI4DucUpXFWQQq/W1vLvR/s4MuPl3HBmAw+PzKdWL+XWL+HGL+XGL+HWL+XGF/ovsdxIl2+iMigopAlMgD5vB7mTczh/NEZPLS0gkUrK3llQ90pnxPt8xDj6wpefg9T8gPcMnO4BtKLiISJPl1FBrCEaB/fPGc4107Jpaapldb2Tg62ddDa1nXb3klrW8ff7GtubeeFNTW8vXUP35ldxOziITg60yUi0qsUskRcIDUuitS4qNN6zqZdLfzkj1v455c3ck5RGt+ZXURWUkyYKhQRGXw08F1kkBqVmchD15byrXML+XhHA1c/vIJnyqro6NQ3FkVEekO3zmQZY+YCdwFe4H5r7U+Pe/xm4FagA9gH3GSt3dD12Hjgf4EkoBOYaq1t7bXfQEQ+M5/H4dopecwqTuOnb27l529t47WNdfzreSMZkR4f6fJERAa0Tz2TZYzxAncD5wNjgGuMMWOOO+xJa+04a+1E4GfAr7qe6wMeB2621pYAs4C23itfRHpDbnIs/3PFWH58wSiqGlu57vGV/Oa9T2ht64h0aSIiA1Z3zmRNA7Zaa8sBjDGLgEuBDUcOsNYeu4ptPHDkesMcYI21dnXXcXt6o2gR6X2O4zB3dAZnFKRw17vlPLS0gjdtPd87byRThgYiXZ6IyIDTnTFZuUDFMduVXfv+ijHmVmPMNkJnsv6+a/dIIGiMed0Ys9IY8089LVhEwisQ6+cHcw2/mTeOIPB3z67hR3+wNB7QSWgRkdPRnTNZJ/pe99+MjLXW3g3cbYxZAHwfWNj182cCU4EDwJ+MMSustX862Yt5vQ6BQFx3au8Rr9fTJ68jkaH+9tx5gTjOHp3F3e9s4/73PuEPm+o4tzidC8dlM3tUOnFRkf1ysnrsfuqx+7m9x935lKwE8o/ZzgOqT3H8IuCeY577rrV2N4Ax5lVgEnDSkNXREaSx8UA3yuqZQCCuT15HIkP97T1fnZrHrOEpvLS2ljc31/PmpjpifB5mFqYxZ1Q6Zw1PJdrX919UVo/dTz12Pzf0OD098aSPdSdkLQOKjTHDgSpgPrDg2AOMMcXW2i1dmxcCR+6/DvyTMSYOOAycC9x5WtWLSMQVDYnn9s8V8a1ZhayqauKPm+r50+bdvLm5nvgoL+eOSGOOyWDasAB+r2aGERGBboQsa227MeY2QoHJCzxorV1vjPkRsNxauxi4zRjzBULfHGwgdKkQa22DMeZXhIJaEHjVWvtKmH4XEQkzj+MwKS/ApLwAd8wewfKdDbxh63l7yx5e3VBHcoyPWcVDmGPSmZwfwOvRLPIiMng5wWD/mniwra0jqMuF0lPqb9863N7JRztCgevdrbs52NZJTnIMt509nC+MDM+SPeqx+6nH7ueGHqenJ64AppzoMS2rIyI9FuXzcE5RGucUpdHa1sGfy/fy0NKdfO/3G3kqO4lvzSpkfE5SpMsUEelTGjwhIr0qxu/lPJPOY9dN4vtziqlubuWrT63iX17eSFXTwUiXJyLSZ3QmS0TCwutxuHRcNueZDB5bVsFjyyt5d9tu5pfm8pXpQ0mM0cePiLibzmSJSFjFRXn5xowCXrhxKl8clcHjyyu5/IGPeaasivaOzkiXJyISNgpZItInMhKj+cFcw2PXTaI4PZ6fv7WN+Y+s4N2te+hvX8AREekNClki0qdMZgK/mTeeX15WAsAdL63nlmfXsKG2JcKViYj0Lg2KEJE+5zgO5xSlcVZBCi+sqeW+D3ew8IkyJucnM780l7OL0jTHlogMeApZIhIxPq+Hq0pzuGBMBr9bU8Ozq6r5zuIN5CTHcHVpDpeMzSIhWh9TIjIwaTJScSX1d2Bq7wyyZOtuFq2soqyqmTi/l4vHZnJVaS5DU2L/6lj12P3UY/dzQ481GamIDAg+j8PskenMHpnOpl0tLFpZxfOra3imrJoZhanMn5TLtKGBsMwgLyLS23QmS1xJ/XWP3fsP88Lqap5fXcPeA20UpsUxf1Iu888o4NCBQ5EuT8JI72P3c0OPT3UmSyFLXEn9dZ9D7Z38cVMdi1ZWsbl+P4FYP5ePz2LexBzSE6IjXZ6Egd7H7ueGHutyoYgMeNE+DxePzeKikkxWVjbx/NpaHl5awWPLKpkzKp0Fk/IwmQmRLlNE5CiFLBEZUBzHYXJ+gM+Py2Hd9j0sWlnF4nW1vLqhjkl5ySyYnMfZRal4NG5LRCJMIUtEBqy8QCx3zB7BN84q4MW1NTxdVs0dL60nPxDD/Em5XFSSRVyUN9JlisggpTFZ4krqr/udqMftnUHe2lzPUyurWFfTQmK0j8vGZXFVaQ5ZSTERqlQ+K72P3c8NPdaYLBEZFHwehzmjMpgzKoM11c08taKSJ1ZU8uSKSs4ankpOcgzJsX4CsX5Sum5D/3wkx/rxe7XSmIj0HoUsEXGl8TlJjM8ZQ3VTK0+XVfFe+V5WVTXTcqj9pM+Jj/IeDV5ZSdFcPyWPkuykPqxaRNxEIUtEXC0nOYZvzyri27OKAGjv6KSxtZ3Gg200HWyj8Zh/DQfauva3s6KiiT9t3s3c0RncOrNAlxtF5LQpZInIoOLzehgSH8WQ+KhTHrf/cDuPfFzBE8sreXvLbq6dksfCqfkaSC8i3aYBCCIiJxAf5eOWmcN57sapzBqRxoMf7eSKB5exeG0tHZ396wtDItI/KWSJiJxCdlIM/3HhaB68ZiI5STH8+I+buf7xlSzb2RDp0kSkn1PIEhHphnE5STxwzQR+cuEo9h1q55Zn13L779axfe/A/vq5iISPQpaISDc5TmiKiGe/MpVbZxawsrKJ+Y+s4BdvbaXxYFukyxORfkYD30VETlO0z8MN04dy8dgs7v1gB8+uqub363cxdWiA0rxkSvOSKU5PwOfR0j4ig5lClojIZ5QWH8W/nFfMvNIcnlxeycrKJt7ZugeAOL+X8blJlOaGQteYrESifbp4IDKYKGSJiPTQiCHx/J+5BoC6lkOsqmpiZWUTq6qauOf97QD4vQ5jsxKZmJfMxNxkxuckkRCtj2ARN9M7XESkF2UkRh9d2geg8WAbq6uaKesKXY9+XMFDwQq8Dpw3KoOFU/MZkR4f4apFJBwUskREwigQ6+fcEWmcOyINgAOHO1hb08x75Xt5aW0Nf9hYx8zCVBZOzWdiXnKEqxWR3qSQJSLSh+KivEwflsL0YSl87YyhPLuqmqfLqvn606uZkJPEwmn5zChMxeNo0LzIQKeQJSISIcmxfr525jCum5LH4nW1PL68kttfXE/RkDi+PDWfOSYdn1eD5UUGKr17RUQiLMbv5arSXF64cSr/fr4hGIQfvGa5/IFlPL2yita2jkiXKCKfgc5kiYj0Ez6vhwvGZDJ3dAbvl+/lkY8r+MXb27j/o51cVZrDrBFpFKbF49X8WyIDgkKWiEg/43Eczi5K4+yiNFZVNvHIsgru/WAH936wg1i/h9GZiZRkJTI2O5GS7CQyEqJwNIZLpN9RyBIR6ccm5iUzMS+ZqqaDrK5qZkNtC+tqWlhUVkXb8iAAQ+KjQoErK5GS7ERGZyZqDi6RfkDvQhGRASA3OZbc5FguGJMJwOH2TrbU72NdTQvra0P/jsw27wDD0+KYOjTA2YVplOYlE6XZ5kX6nEKWiMgAFOXzUJKdREl20tF9TQfb2LArdKZrbXUzL66t5emyauL8XqYNCwWuswpTGRIfFcHKRQYPhSwREZdIjvVzZkEqZxakAtDa1sHyikbeK9/Le+V7j57pGp2ZwNmFacwoTGVUZoLm5BIJE4UsERGXivF7mVmYxszCNILBIFt37+e98r38edte7vtwB/d+uIO0+ChmDk9lRmEqZw1P1SLWIr1IIUtEZBBwHIfi9ASK0xP4yvShNBw4zIfbG/jztr28ubmel9bVMiQ+ioXT8rlsXBYxfm+kSxYZ8BSyREQGoZS4KC4Yk8kFYzJp7+hkWUUjDy+t4Jdvb+PhjytYOC2fyxW2RHqkWyHLGDMXuAvwAvdba3963OM3A7cCHcA+4CZr7YZjHh8KbAB+aK39RS/VLiIivcDn9Rwdy7WiopH7PtzBr97exiMKWyI98qkX340xXuBu4HxgDHCNMWbMcYc9aa0dZ62dCPwM+NVxj98JvNYL9YqISBhNzg/w26sm8NurxlOQGsuv3t7GZQ8s48kVlVreR+Q0dedM1jRgq7W2HMAYswi4lNCZKQCstc3HHB8PBI9sGGMuA8qB/b1RsIiIhN/k/ACT8wOsqGjk/g93cOc75Ty6rJIvT83jivHZOrMl0g3dCVm5QMUx25XA9OMPMsbcCtwORAGzu/bFA98FzgPu6GmxIiLSt04Uto5cRlTYEjm17oSsE02gEjx+h7X2buBuY8wC4PvAQuDfgTuttfuMMd0qyOt1CATiunVsT3i9nj55HYkM9df91OO+9flAHJ8fl8PST/by67e3cuc75Ty4tIIrJ+WyYNpQhqb2fi/UY/dze4+7E7IqgfxjtvOA6lMcvwi4p+v+dOBKY8zPgADQaYxptdb++mRP7ugI0th4oBtl9UwgENcnryORof66n3ocGSYlhv93xVhWVTbxdFkVD3+wnQff386MwlTmTczhjIKUXpvcVD12Pzf0OD098aSPdSdkLQOKjTHDgSpgPrDg2AOMMcXW2i1dmxcCWwCstWcfc8wPgX2nClgiIjIwHFm4uq7lEC+sqeF3a2r4hxfWkR+I4cqJOVxckkVijGYJksHtU98B1tp2Y8xtwOuEpnB40Fq73hjzI2C5tXYxcJsx5gtAG9BA6FKhiIi4XEZiNDfPKOCrZwzlrc27eWZVNXe+U849723n/DEZzJuYQ3F6QqTLFIkIJxj8m+FVEdXW1hHU5ULpKfXX/dTj/svu2sezq6r5w6Y6DrV3UpqXzLyJOcwakYbf2/1le9Rj93NDj9PTE1cAU070mEKWuJL6637qcf/XdLCNxetqeW51DdVNrfi9DkVp8YzMiMdkJDAyPYHijHjio058UUU9dj839PhUIUsXzEVEJCySY/1cPzWfBZPz+GhHAysrGrF1+1iybS+L1+06elx+ICYUurqCl8mIJy0+KoKVi/QOhSwREQkrr8dhxvBUZgxPBSAYDFK/7zCb6/dh6/axuW4/G3ft483Nu48+JzXOz4wRQ5g3PovRmSf/9pZIf6aQJSIifcpxHDISo8lIjGZmYdrR/fsOtbOlfj+2bh+b6vbxtq3n5TU1TM5P5ropeZw1PLXXpocQ6QsKWSIi0i8kRPsozUumNC8ZAG+Mn0fe+4SnVlTy7d+tZ3haHNdNzmPu6AyifN0fQC8SKfqvVERE+qXEGD/XTcnjpa9N40cXGHwehx//cTOX3P8xDy3dSdPBtkiXKHJKOpMlIiL9ms/r4fzRmcwdlcHHOxt5fHklv3lvOw8t3cklY7O4ZnIuucmxkS5T5G8oZImIyIDgOA7Th6UwfVgKW+r38cSKKp5fXcOzq6qZXZzOlyZkMy4niWhdSpR+QiFLREQGnOL0BH4413DLjAKeLguFrTc31xPldSjJSqS0a9mf8TlJJ52HSyTcNBmpuJL6637qsfudTo/3H25n+c4mVlU1UVbZxKZdLXQEweOAyUgIDajPTWZibjKBOH+YK5fucsP7WJORioiIq8VH+Th3RBrnjghNCXHgcAdra5opqwwFr+dX1/DkiioAhqfFUZqbzKjMBFJi/STH+kmO9ZEc4yc5xofvNJb+ETkVhSwREXGduCjv0fFbAIfbO9m4q4WyyibKqpp4fVMdL6ypOeFz46O8oeAV4zt6G4j1k50Uw3kmnYzE6L78VWQAU8gSERHXi/J5mJCbzITcZG4AOjqD1O07RNPBNpoOttPU2kZj123TwTaaWtuP3lY0HKSptY19hzr4nyXlTBuWwsUlmZw7YogG2cspKWSJiMig4/U4ZCfFkJ0U0+3nVDYe5OX1u3hl/S7+9ZVNJMX4mGPSuXhsFqMzE3A0G70cRyFLRESkG/ICsfzdjAJuOnMYyysaeXldLS+v38Vzq2soGhLHRSVZnD86Q4tby1EKWSIiIqfB6/nLfF0tre28Yev4/fpd3PVuOb/+8yfMGJ7KxSWZzChMxa9B9IOaQpaIiMhnlBjj44oJOVwxIYdP9hzg9+treWVDHUu27SE1zs/CaflcMT6bGL830qVKBGieLHEl9df91GP3G6g9bu8MsnR7A0+uqOTjnY2kJ0Rxw7ShXDYuSwtbH2eg9vhYmidLRESkj/g8DjMKU5lRmMqKikb+9/3t/PytrTy2rIIbzxjKxSWZmotrkFCXRUREwmRyfoD/vXoCv/7SOIYkRPGfb2zhyoeW8/v1tbR39q8rSdL7FLJERETCyHEcphek8OA1E7nz8hISon38+x82M//h5by+sY7OfjZsR3qPQpaIiEgfcByHmYVpPHZdKT+7ZAw+r8P3X93EgkdX8NaW3fS3MdLScxqTJSIi0occx+FzxUM4d0Qab9p67v1gB99dvIHi9HhmjUhj2tAUxmYnatyWCyhkiYiIRIDHcZgzKoPZI9N5fWMdz6yq5v4Pd3LfhzuJ9XuYlBdg2rAAU4cGKBoSj0czyg84ClkiIiIR5PM4XFiSyYUlmTQdbGNFZRPLdjTw8c5G3v9kLwCpcX6m5IdC17RhKae1HJBEjkKWiIhIP5Ec62d28RBmFw8BoLa5lWU7G/l4ZyPLdjbyR1sPQF4ghmlDU5g+LMCUoQGSYvyRLFtOQiFLRESkn8pKiuHisVlcPDaLYDBI+Z4DodC1o4HXN9XxwpoaPA6MyUpk+rAUzhim8Vz9iUKWiIjIAOA4DkVD4ikaEs/8Sbm0d3SyvraFj7Y3sHRHAw8t3ckDH+0kPsrL5PxAKHQVpJAfiMHReK6IUMgSEREZgHxeDxNyk5mQm8w3ZhTQ0trOsorQWa4PtzewZNseALKTopk+LIWpQwMkRIf+7B/JXEeil8Nf73AILYRdkpWodRd7QCFLRETEBRJjfH81nquy8SBLdzTw0fYG3rD1vLi29rR/Zmp0rcEVAAAI5ElEQVScn6tLc7lyYrbGfX0GClkiIiIulBeIJS8Qy5cm5NDeGWTb7v0cbu/kyJSnx09+emTzyN7m1naeX13NPe9v55GPK7hsfBbXTMolS99s7DaFLBEREZfzeRxMRsJpP+/cEWlsqd/Ho8sqeXplFU+XVTN3dAbXT8mjaEh8GCp1F4UsEREROani9AR+fMEobplZwBPLK3lpbS2vrN/FzMJUFk7NZ2JecqRL7Lf0HU8RERH5VNlJMdwxewQv3zSdm84axrqaFr7+9Gq++tQq3t26Wwtdn4DOZImIiEi3BWL9fP3MYVw/JY/F63bxxPIK7nhpAwWpsVxcksU5RWkMS43VtBGA099W/W5r6wg2Nh4I++sEAnH0xetIZKi/7qceu596PDC0dwZ5a3M9T66oYn1tCwD5gRjOLkrjnKI0JuQm4/OcOHC5ocfp6YkrgCknekxnskREROQz83lCC13PGZVBbXMr75Xv5c/le3h2VTVPrqgiKcbHWcNTOacojTMLUo7O1TUYDJ7fVERERMIqKymGKyfmcOXEHPYfbmfpjkaWbNvDe9v28IeNdfg8DpPzkzm7MI2zi9IIBOIiXXJY6XKhuJL6637qsfupx+7R0RlkXU0zS7btYcm2PWzfexCA/JRYxmQmUJKdxNisRExGAlG+gfWdPF0uFBERkYjxepyjSwB985xCduw9wHvle9lQv5+ynQ28vqkeCF16HJmRwNisREqyExmbnTSg115UyBIREZE+NSw1jmGpcUfPVta1HGJdbQvra5pZV9PCy+treWZVNQDJMT7GZCUyNjuRGcNTGZOVOGBCV7dCljFmLnAX4AXut9b+9LjHbwZuBTqAfcBN1toNxpjzgJ8CUcBh4DvW2rd6sX4REREZ4DISo5mdGH103cX2ziCf7NnP2pq/BK+Ptjdw34c7KRoSxyVjs7hgdCaBuP69nuKnhixjjBe4GzgPqASWGWMWW2s3HHPYk9ba33YdfwnwK2AusBu42FpbbYwZC7wO5Pby7yAiIiIu4vM4FKcnUJyewBXjswFoaW3njc31vLS2ljvfKefXf/6Ec4vSuGRcFtOGpuA9yTQRkdSdM1nTgK3W2nIAY8wi4FLgaMiy1jYfc3w8XetLWmvLjtm/HogxxkRbaw/1tHAREREZPBJjfFwxPpsrxmeztX4/L62r5bUNu3hz826yEqO5qCSTi8dmkZPcfxaw7k7IygUqjtmuBKYff5Ax5lbgdkKXBmef4Od8CShTwBIREZGeGJEezz9+rohvnj2cd7ftYfHaWh74aCcPfLSTacMCXDI2i1kjhkT8m4rdCVknOv/2N/M+WGvvBu42xiwAvg8sPPKYMaYE+C9gzqe9mNfr9Mm8GV6vx/Xzcwxm6q/7qcfupx67X2/0eN6QBOZNH0ZV40FeWFnF82WV/OsrmwjE+vnBRWO4qOtyYyR0J2RVAvnHbOcB1ac4fhFwz5ENY0we8Dvgy9babZ/2Yh0dwT6ZF0Xzr7ib+ut+6rH7qcfu15s9jgeun5TDtaXZLNvRyKsbd9GyvzXs/w2lpyee9LHuhKxlQLExZjhQBcwHFhx7gDGm2Fq7pWvzQmBL1/4A8ArwL9ba90+/dBEREZHu8zgO0wtSmF6QEulSPj1kWWvbjTG3EfpmoBd40Fq73hjzI2C5tXYxcJsx5gtAG9DAXy4V3gaMAP7NGPNvXfvmWGvrevsXEREREelPtKyOuJL6637qsfupx+7nhh6falmdgbVAkIiIiMgAoZAlIiIiEgYKWSIiIiJhoJAlIiIiEgYKWSIiIiJhoJAlIiIiEgYKWSIiIiJhoJAlIiIiEgYKWSIiIiJhoJAlIiIiEgYKWSIiIiJhoJAlIiIiEgYKWSIiIiJhoJAlIiIiEgYKWSIiIiJhoJAlIiIiEgZOMBiMdA3Hqwd2RLoIERERkW4YBqSf6IH+GLJEREREBjxdLhQREREJA4UsERERkTBQyBIREREJA4UsERERkTBQyBIREREJA1+kC+hrxpi5wF2AF7jfWvvTCJckPWSMeRC4CKiz1o7t2pcKPA0UANuBq6y1DZGqUXrGGJMPPApkAZ3Avdbau9Rn9zDGxABLgGhCf5ues9b+wBgzHFgEpAIrgeuttYcjV6n0hDHGCywHqqy1F7m9v4PqTFZXc+8GzgfGANcYY8ZEtirpBQ8Dc4/b98/An6y1xcCfurZl4GoH/tFaOxo4A7i1672rPrvHIWC2tXYCMBGYa4w5A/gv4M6uHjcAX41gjdJz/wBsPGbb1f0dVCELmAZstdaWdyXlRcClEa5JeshauwTYe9zuS4FHuu4/AlzWp0VJr7LW1lhrV3bdbyH0IZ2L+uwa1tqgtXZf16a/618QmA0817VfPR7AjDF5wIXA/V3bDi7v72ALWblAxTHblV37xH0yrbU1EPoDDWREuB7pJcaYAqAUWIr67CrGGK8xZhVQB7wBbAMarbXtXYfoM3tg+2/gnwhd8gdIw+X9HWwhyznBPk15LzJAGGMSgOeBb1lrmyNdj/Qua22HtXYikEfoysPoExymz+wByBhzZNzsimN2u/5v8mALWZVA/jHbeUB1hGqR8NpljMkG6Lqti3A90kPGGD+hgPWEtfaFrt3qswtZaxuBdwiNvwsYY458SUuf2QPXDOASY8x2QkN1ZhM6s+Xq/g62kLUMKDbGDDfGRAHzgcURrknCYzGwsOv+QuClCNYiPdQ1duMBYKO19lfHPKQ+u4QxJt0YE+i6Hwt8gdDYu7eBK7sOU48HKGvtv1hr86y1BYT+9r5lrb0Wl/d30C0QbYy5gFB69gIPWmt/EuGSpIeMMU8Bs4AhwC7gB8CLwDPAUGAnMM9ae/zgeBkgjDEzgT8Da/nLeI7vERqXpT67gDFmPKGBz15CJwCesdb+yBhTyF++4l8GXGetPRS5SqWnjDGzgDu6pnBwdX8HXcgSERER6QuD7XKhiIiISJ9QyBIREREJA4UsERERkTBQyBIREREJA4UsERERkTBQyBIREREJA4UsERERkTBQyBIREREJg/8PgJ6OFKh/giEAAAAASUVORK5CYII=\n", + "text/plain": [ + "<Figure size 720x576 with 1 Axes>" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = plt.subplots(figsize=(10, 8))\n", + "n_epochs = len(net.history)\n", + "sns.lineplot(range(n_epochs), net.history[:, 'train_loss'])\n", + "sns.lineplot(range(n_epochs), net.history[:, 'valid_loss'])\n", + "ax.legend(['train_loss', 'valid_loss'])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T02:07:54.326248Z", + "start_time": "2019-08-10T02:07:51.779447Z" + }, + "scrolled": true + }, + "outputs": [], + "source": [ + "prob = net.predict_proba(x_test)\n", + "y_pred = net.predict(x_test)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T02:07:55.088517Z", + "start_time": "2019-08-10T02:07:54.327862Z" + } + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(20, 10))\n", + "plot_youden(ax, y_test, prob, 0.1, 0.9, 40)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T02:09:14.799918Z", + "start_time": "2019-08-10T02:09:13.571047Z" + } + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(20, 10))\n", + "plot_thresh_range(ax, y_test, prob, 0.1, 0.9, 40)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T02:09:22.235184Z", + "start_time": "2019-08-10T02:09:21.973392Z" + }, + "scrolled": false + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(10,8))\n", + "plot_roc(ax, y_test, prob)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T02:09:25.384728Z", + "start_time": "2019-08-10T02:09:25.127586Z" + } + }, + "outputs": [], + "source": [ + "threshold = 0.2\n", + "y_pred = (prob > threshold).astype(np.int64)\n", + "cm = confusion_matrix(y_test, y_pred)\n", + "tn,fp,fn,tp = cm[0][0],cm[0][1],cm[1][0],cm[1][1]\n", + "sensitivity = tp/(tp+fn)\n", + "specificity = tn/(tn+fp)\n", + "ppv = tp/(tp+fp)\n", + "npv = tn/(tn+fn)\n", + "f1 = (2*ppv*sensitivity)/(ppv+sensitivity)\n", + "auroc = roc_auc_score(y_test, prob)\n", + "\n", + "d = {\n", + " 'sensitivity': np.round(sensitivity, 3),\n", + " 'specificity': np.round(specificity, 3),\n", + " 'ppv': np.round(ppv, 3),\n", + " 'npv': np.round(npv, 3),\n", + " 'f1': np.round(f1, 3),\n", + " 'auroc': np.round(auroc, 3),\n", + "}\n", + "\n", + "metrics = pd.DataFrame(d.values(), index=d.keys(), columns=['Value'])\n", + "metrics" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T02:10:59.034655Z", + "start_time": "2019-08-10T02:10:58.573846Z" + } + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(10, 8))\n", + "plot_confusion_matrix(ax, cm, classes=['Delayed', 'Imminent'], normalize=False, title='Confusion matrix')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Prolonged ICU Stay" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T02:11:48.636198Z", + "start_time": "2019-08-10T02:11:48.514559Z" + } + }, + "outputs": [], + "source": [ + "ori_df = set_group_splits(ps_df.copy(), group_col='hadm_id', seed=seed)\n", + "df = ori_df\n", + "# df = ori_df.sample(1000).reset_index(drop=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T02:12:51.529544Z", + "start_time": "2019-08-10T02:12:18.693490Z" + } + }, + "outputs": [], + "source": [ + "vectorizer = TfidfVectorizer(sublinear_tf=True, ngram_range=(1,2), binary=True, max_features=60_000)\n", + "\n", + "x_train = vectorizer.fit_transform(df.loc[(df['split'] == 'train')]['processed_note']).astype(np.float32)\n", + "x_test = vectorizer.transform(df.loc[(df['split'] == 'test')]['processed_note']).astype(np.float32)\n", + "\n", + "x_train = np.asarray(x_train.todense())\n", + "x_test = np.asarray(x_test.todense())\n", + "\n", + "vocab_sz = len(vectorizer.vocabulary_)\n", + "\n", + "y_train = df.loc[(df['split'] == 'train')]['prolonged_stay_label'].to_numpy()\n", + "y_test = df.loc[(df['split'] == 'test')]['prolonged_stay_label'].to_numpy()\n", + "df.shape, x_train.shape, x_test.shape, y_train.shape, y_test.shape" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T02:13:13.127114Z", + "start_time": "2019-08-10T02:13:12.459632Z" + } + }, + "outputs": [], + "source": [ + "train_ds = TensorDataset(torch.tensor(x_train), torch.tensor(y_train.astype(np.float32)))\n", + "train_dl = DataLoader(train_ds, batch_size=args.batch_size, shuffle=True, drop_last=True)\n", + "itr = iter(train_dl)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T02:13:13.727112Z", + "start_time": "2019-08-10T02:13:13.663448Z" + } + }, + "outputs": [], + "source": [ + "clf = MLPModule(input_units=vocab_sz, output_units=1, hidden_units=args.hidden_dim, num_hidden=1, dropout=args.dropout_p, squeeze_output=True)\n", + "loss_fn = nn.BCEWithLogitsLoss()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T02:13:14.414141Z", + "start_time": "2019-08-10T02:13:14.370681Z" + } + }, + "outputs": [], + "source": [ + "x, y = next(itr)\n", + "y_pred = clf(x)\n", + "loss_fn(y_pred, y)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T02:13:15.663594Z", + "start_time": "2019-08-10T02:13:15.631754Z" + } + }, + "outputs": [], + "source": [ + "reduce_lr = LRScheduler(\n", + " policy='ReduceLROnPlateau',\n", + " mode='min',\n", + " factor=0.5,\n", + " patience=1,\n", + ")\n", + "\n", + "checkpoint = Checkpoint(\n", + " dirname=args.workdir/'models/ps_dev_run1',\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T02:13:16.561491Z", + "start_time": "2019-08-10T02:13:16.528937Z" + } + }, + "outputs": [], + "source": [ + "net = NeuralNetBinaryClassifier(\n", + " clf,\n", + " max_epochs=args.max_epochs,\n", + " lr=args.lr,\n", + " device=args.device,\n", + " optimizer=optim.Adam,\n", + " optimizer__weight_decay=args.wd,\n", + " batch_size=args.batch_size,\n", + " verbose=1,\n", + " callbacks=[EarlyStopping, checkpoint, reduce_lr],\n", + " train_split=CVSplit(cv=0.15, stratified=True),\n", + " iterator_train__shuffle=True, \n", + ")\n", + "\n", + "net.set_params(callbacks__valid_acc=None);" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T02:19:50.994267Z", + "start_time": "2019-08-10T02:13:18.983380Z" + }, + "scrolled": false + }, + "outputs": [], + "source": [ + "net.fit(x_train, y_train.astype(np.float32))\n", + "# net.initialize()\n", + "# net.load_params(checkpoint=checkpoint)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T02:19:51.519879Z", + "start_time": "2019-08-10T02:19:50.996110Z" + } + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(10, 8))\n", + "n_epochs = len(net.history)\n", + "sns.lineplot(range(n_epochs), net.history[:, 'train_loss'])\n", + "sns.lineplot(range(n_epochs), net.history[:, 'valid_loss'])\n", + "ax.legend(['train_loss', 'valid_loss'])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T02:19:54.687943Z", + "start_time": "2019-08-10T02:19:51.521325Z" + }, + "scrolled": true + }, + "outputs": [], + "source": [ + "prob = net.predict_proba(x_test)\n", + "y_pred = net.predict(x_test)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T02:19:55.501191Z", + "start_time": "2019-08-10T02:19:54.689503Z" + } + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(20, 10))\n", + "plot_youden(ax, y_test, prob, 0.1, 0.9, 40)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T02:29:33.089367Z", + "start_time": "2019-08-10T02:29:31.801492Z" + } + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(20, 10))\n", + "plot_thresh_range(ax, y_test, prob, 0.1, 0.9, 40)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T02:29:40.317871Z", + "start_time": "2019-08-10T02:29:40.049324Z" + }, + "scrolled": false + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(10,8))\n", + "plot_roc(ax, y_test, prob)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-08-10T02:30:01.696081Z", + "start_time": "2019-08-10T02:30:01.657269Z" + } + }, + "outputs": [], + "source": [ + "threshold = 0.25\n", + "y_pred = (prob > threshold).astype(np.int64)\n", + "cm = confusion_matrix(y_test, y_pred)\n", + "tn,fp,fn,tp = cm[0][0],cm[0][1],cm[1][0],cm[1][1]\n", + "sensitivity = tp/(tp+fn)\n", + "specificity = tn/(tn+fp)\n", + "ppv = tp/(tp+fp)\n", + "npv = tn/(tn+fn)\n", + "f1 = (2*ppv*sensitivity)/(ppv+sensitivity)\n", + "auroc = roc_auc_score(y_test, prob)\n", + "\n", + "d = {\n", + " 'sensitivity': np.round(sensitivity, 3),\n", + " 'specificity': np.round(specificity, 3),\n", + " 'ppv': np.round(ppv, 3),\n", + " 'npv': np.round(npv, 3),\n", + " 'f1': np.round(f1, 3),\n", + " 'auroc': np.round(auroc, 3),\n", + "}\n", + "\n", + "metrics = pd.DataFrame(d.values(), index=d.keys(), columns=['Value'])\n", + "metrics" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-07-16T17:27:21.639571Z", + "start_time": "2019-07-16T17:27:21.397075Z" + } + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(10, 8))\n", + "plot_confusion_matrix(ax, cm, classes=['Delayed', 'Imminent'], normalize=False, title='Confusion matrix')" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "heading_collapsed": true + }, + "source": [ + "## Metrics" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "hidden": true + }, + "source": [ + "### Imminent ICU Admission" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-07-22T11:55:29.364845Z", + "start_time": "2019-07-22T11:55:28.722095Z" + }, + "hidden": true + }, + "outputs": [], + "source": [ + "with open(args.workdir/f'imminent_adm_preds.pkl', 'rb') as f:\n", + " targs = pickle.load(f)\n", + " preds = pickle.load(f)\n", + " probs = pickle.load(f)\n", + "\n", + "fnames = [f'imminent_adm_seed_{seed}.pkl' for seed in range(args.start_seed, args.start_seed + 100)]\n", + "bam = BinaryAvgMetrics(targs, preds, probs)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-07-22T11:55:29.586114Z", + "start_time": "2019-07-22T11:55:29.366363Z" + }, + "hidden": true + }, + "outputs": [], + "source": [ + "bam.get_avg_metrics(conf=0.95)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-07-22T11:55:29.818141Z", + "start_time": "2019-07-22T11:55:29.587526Z" + }, + "hidden": true + }, + "outputs": [], + "source": [ + "get_best_model(bam, fnames)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-07-22T11:55:30.206128Z", + "start_time": "2019-07-22T11:55:29.819466Z" + }, + "hidden": true + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(10, 8))\n", + "plot_mean_roc(ax, bam.targs, bam.probs)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-07-22T11:55:30.515657Z", + "start_time": "2019-07-22T11:55:30.207420Z" + }, + "hidden": true + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(1, 2, figsize=(15, 6))\n", + "\n", + "plot_confusion_matrix(ax[0], bam.cm_avg, classes=['not imminent', 'imminent'], normalize=False,\\\n", + " title='Confusion Matrix Over Runs')\n", + "plot_confusion_matrix(ax[1], bam.cm_avg, classes=['not imminent', 'imminent'], normalize=True,\\\n", + " title='Normalized Confusion Matrix Over Runs')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "hidden": true + }, + "source": [ + "### Prolonged ICU Stay" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-07-22T11:55:55.670109Z", + "start_time": "2019-07-22T11:55:54.957696Z" + }, + "hidden": true + }, + "outputs": [], + "source": [ + "with open(args.workdir/f'prolonged_stay_preds.pkl', 'rb') as f:\n", + " targs = pickle.load(f)\n", + " preds = pickle.load(f)\n", + " probs = pickle.load(f)\n", + "\n", + "fnames = [f'prolonged_stay_seed_{seed}.pkl' for seed in range(args.start_seed, args.start_seed + 100)]\n", + "bam = BinaryAvgMetrics(targs, preds, probs)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-07-22T11:55:59.488758Z", + "start_time": "2019-07-22T11:55:59.229856Z" + }, + "hidden": true + }, + "outputs": [], + "source": [ + "bam.get_avg_metrics(conf=0.95)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-07-22T11:56:04.618880Z", + "start_time": "2019-07-22T11:56:04.366173Z" + }, + "hidden": true + }, + "outputs": [], + "source": [ + "get_best_model(bam, fnames)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-07-22T11:56:08.900406Z", + "start_time": "2019-07-22T11:56:08.487768Z" + }, + "hidden": true + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(10, 8))\n", + "plot_mean_roc(ax, bam.targs, bam.probs)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-07-22T11:56:16.786590Z", + "start_time": "2019-07-22T11:56:16.425934Z" + }, + "hidden": true + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(1, 2, figsize=(15, 6))\n", + "\n", + "plot_confusion_matrix(ax[0], bam.cm_avg, classes=['Discharge within 5 days', 'Discharge after 5 days'], normalize=False, title='Confusion Matrix Over Runs')\n", + "plot_confusion_matrix(ax[1], bam.cm_avg, classes=['Discharge within 5 days', 'Discharge after 5 days'], normalize=True, title='Normalized Confusion Matrix Over Runs')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Full Run" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-07-22T12:06:04.466706Z", + "start_time": "2019-07-22T12:06:02.766145Z" + } + }, + "outputs": [], + "source": [ + "seed = 643\n", + "ori_df = pd.read_csv(args.dataset_csv, usecols=args.cols, parse_dates=args.dates)\n", + "ori_df['relative_charttime'] = (ori_df['charttime'] - ori_df['intime'])\n", + "\n", + "ia_df = ori_df.loc[(ori_df['imminent_adm_label'] != -1)][args.imminent_adm_cols + ['relative_charttime']].reset_index(drop=True)\n", + "\n", + "ps_df = ori_df.loc[(ori_df['chartinterval'] != 0)][args.prolonged_stay_cols + ['relative_charttime']].reset_index(drop=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-07-22T12:06:04.491622Z", + "start_time": "2019-07-22T12:06:04.468225Z" + } + }, + "outputs": [], + "source": [ + "interval_hours = 12\n", + "starting_day = -20\n", + "ending_day = -1" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Imminent ICU Admission" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-07-22T12:06:04.550287Z", + "start_time": "2019-07-22T12:06:04.492938Z" + } + }, + "outputs": [], + "source": [ + "df = set_group_splits(ia_df.copy(), pct=0.25, group_col='hadm_id', seed=seed)\n", + "df['prob'] = -1" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-07-22T12:06:13.284651Z", + "start_time": "2019-07-22T12:06:04.551535Z" + } + }, + "outputs": [], + "source": [ + "vectorizer = TfidfVectorizer(min_df=args.min_freq, binary=True, analyzer=str.split, sublinear_tf=True)\n", + "\n", + "x_train = vectorizer.fit_transform(df.loc[(df['split'] == 'train')]['processed_note']).astype(np.float32)\n", + "x_test = vectorizer.transform(df.loc[(df['split'] == 'test')]['processed_note']).astype(np.float32)\n", + "\n", + "x_train = np.asarray(x_train.todense())\n", + "x_test = np.asarray(x_test.todense())\n", + "\n", + "vocab_sz = len(vectorizer.vocabulary_)\n", + "\n", + "y_train = df.loc[(df['split'] == 'train')]['imminent_adm_label'].to_numpy()\n", + "y_test = df.loc[(df['split'] == 'test')]['imminent_adm_label'].to_numpy()\n", + "\n", + "df.shape, x_train.shape, x_test.shape, y_train.shape, y_test.shape" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-07-22T12:04:38.168048Z", + "start_time": "2019-07-22T12:04:37.823108Z" + } + }, + "outputs": [], + "source": [ + "train_ds = TensorDataset(torch.tensor(x_train), torch.tensor(y_train.astype(np.float32)))\n", + "train_dl = DataLoader(train_ds, batch_size=args.batch_size, shuffle=True, drop_last=True)\n", + "itr = iter(train_dl)\n", + "\n", + "clf = MLPModule(input_units=vocab_sz, output_units=1, hidden_units=args.hidden_dim, num_hidden=1, dropout=args.dropout_p, squeeze_output=True)\n", + "loss_fn = nn.BCEWithLogitsLoss()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-07-22T12:04:38.239707Z", + "start_time": "2019-07-22T12:04:38.170461Z" + } + }, + "outputs": [], + "source": [ + "x, y = next(itr)\n", + "y_pred = clf(x)\n", + "\n", + "loss_fn(y_pred, y)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-07-22T12:06:26.700111Z", + "start_time": "2019-07-22T12:06:26.457513Z" + } + }, + "outputs": [], + "source": [ + "reduce_lr = LRScheduler(\n", + " policy='ReduceLROnPlateau',\n", + " mode='min',\n", + " factor=0.5,\n", + " patience=1,\n", + ")\n", + "\n", + "checkpoint = Checkpoint(\n", + " dirname=args.workdir/'models/ia_full_run_01',\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-07-22T12:06:51.913709Z", + "start_time": "2019-07-22T12:06:51.852150Z" + } + }, + "outputs": [], + "source": [ + "clf = MLPModule(input_units=vocab_sz, output_units=1, hidden_units=args.hidden_dim, num_hidden=1, dropout=args.dropout_p, squeeze_output=True)\n", + "\n", + "args.batch_size=64\n", + "args.device='cuda:2'" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-07-22T12:06:51.940021Z", + "start_time": "2019-07-22T12:06:51.915151Z" + } + }, + "outputs": [], + "source": [ + "net = NeuralNetBinaryClassifier(\n", + " clf,\n", + " max_epochs=args.max_epochs,\n", + " lr=args.lr,\n", + " device=args.device,\n", + " optimizer=optim.Adam,\n", + " optimizer__weight_decay=args.wd,\n", + " batch_size=args.batch_size,\n", + " verbose=1,\n", + " callbacks=[EarlyStopping, checkpoint, reduce_lr],\n", + " train_split=CVSplit(cv=0.15, stratified=True),\n", + " iterator_train__shuffle=True, \n", + "# iterator_train__num_workers=4,\n", + "# iterator_train__pin_memory=True,\n", + "# iterator_train__drop_last=True,\n", + "# iterator_valid__num_workers=4,\n", + "# iterator_valid__pin_memory=True,\n", + " threshold=args.ia_thresh,\n", + ")\n", + "\n", + "net.set_params(callbacks__valid_acc=None);" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2019-07-22T12:07:07.392335Z", + "start_time": "2019-07-22T12:06:53.182073Z" + } + }, + "outputs": [], + "source": [ + "net.fit(x_train, y_train.astype(np.float32))\n", + "# net.initialize()\n", + "# net.load_params(checkpoint=checkpoint)" + ] + } + ], + "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.4" + }, + "toc": { + "base_numbering": 1, + "nav_menu": {}, + "number_sections": true, + "sideBar": true, + "skip_h1_title": true, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": false, + "toc_position": {}, + "toc_section_display": true, + "toc_window_display": false + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}