Sequence_RNN_Attention.ipynb

{
 "nbformat": 4,
 "nbformat_minor": 0,
 "metadata": {
  "colab": {
   "provenance": [],
   "collapsed_sections": [
    "kRb_sju1lrks",
    "Ff0VIBrcEghb",
    "j3PODMTQMF1s",
    "aTrlX5jNasqx",
    "VKj5IKANavL7",
    "C5-a_FEtgGae",
    "ICJTt5J-g6Or",
    "zZlXyxkJuigm"
   ]
  },
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "PREDICTING NEXT TIME STEP HEART FAILURE FROM SEQUENCIAL PATIENT DATA"
   ],
   "metadata": {
    "id": "GMixAbHJ3qBl"
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Imports"
   ],
   "metadata": {
    "id": "kRb_sju1lrks"
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "I45d_fOd5KgF",
    "outputId": "9820aea0-0aa2-4a77-d6c2-26dca8ed0596"
   },
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "cpu is available.\n"
     ]
    }
   ],
   "source": [
    "####imports\n",
    "import pandas as pd\n",
    "from google.colab import files\n",
    "import io\n",
    "import numpy as np\n",
    "from tqdm import tqdm\n",
    "import pickle\n",
    "import torch\n",
    "from torch.utils.data import TensorDataset, DataLoader\n",
    "from sklearn.model_selection import train_test_split\n",
    "from tensorflow.keras.preprocessing.sequence import pad_sequences\n",
    "from tqdm import tqdm\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import torch.nn.functional as F\n",
    "from datetime import datetime\n",
    "from sklearn.metrics import roc_auc_score, average_precision_score\n",
    "from sklearn.metrics import precision_recall_curve, roc_curve\n",
    "# model\n",
    "import os\n",
    "import pickle\n",
    "import numpy as np\n",
    "import tensorflow as tf\n",
    "import joblib\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.model_selection import train_test_split, GridSearchCV\n",
    "from sklearn.linear_model import LogisticRegression, LogisticRegressionCV\n",
    "from sklearn.svm import LinearSVC, SVC\n",
    "from sklearn.metrics import classification_report, roc_auc_score, average_precision_score\n",
    "\n",
    "from tensorflow.keras.models import Model, Sequential, load_model\n",
    "import tensorflow.keras.backend as K\n",
    "from tensorflow.keras.optimizers import Adam\n",
    "from tensorflow.keras.layers import LSTM, GRU, SimpleRNN, Dense, Input, Masking, concatenate, LayerNormalization, Add, Dropout\n",
    "from tensorflow.keras.preprocessing.sequence import pad_sequences\n",
    "from tensorflow.keras.callbacks import ModelCheckpoint\n",
    "import tensorflow as tf\n",
    "\n",
    "from sklearn.metrics import mean_absolute_error\n",
    "\n",
    "#from keras_multi_head import MultiHeadAttention\n",
    "import seaborn as sns\n",
    "\n",
    "\n",
    "device = \"cuda\" if torch.cuda.is_available() else \"cpu\"\n",
    "print(f\"{device}\" \" is available.\")\n"
   ]
  },
  {
   "cell_type": "code",
   "source": [
    "###### IMPORTING DATA IN GOOGLE COLAB #######\n",
    "uploaded = files.upload()\n",
    "df = pd.read_csv(io.BytesIO(uploaded['synthetic_data.csv']))\n",
    "df.head(10)\n",
    "\n",
    "#NOTE IMPORTANT THIS IS DATA OF PATIENTS OLDER THAN 40 AND FOR EACH FEATURE WE ARE GIVING THE AGE AT DIAGNOSIS"
   ],
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 512
    },
    "id": "lvBV43uz54qa",
    "outputId": "adc58fbe-7cb3-4dff-ac35-265dfce47e8b"
   },
   "execution_count": null,
   "outputs": [
    {
     "output_type": "display_data",
     "data": {
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ],
      "text/html": [
       "\n",
       "     <input type=\"file\" id=\"files-f0d1235c-2b67-4c0b-8e26-f2f3cfaca089\" name=\"files[]\" multiple disabled\n",
       "        style=\"border:none\" />\n",
       "     <output id=\"result-f0d1235c-2b67-4c0b-8e26-f2f3cfaca089\">\n",
       "      Upload widget is only available when the cell has been executed in the\n",
       "      current browser session. Please rerun this cell to enable.\n",
       "      </output>\n",
       "      <script>// Copyright 2017 Google LLC\n",
       "//\n",
       "// Licensed under the Apache License, Version 2.0 (the \"License\");\n",
       "// you may not use this file except in compliance with the License.\n",
       "// You may obtain a copy of the License at\n",
       "//\n",
       "//      http://www.apache.org/licenses/LICENSE-2.0\n",
       "//\n",
       "// Unless required by applicable law or agreed to in writing, software\n",
       "// distributed under the License is distributed on an \"AS IS\" BASIS,\n",
       "// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n",
       "// See the License for the specific language governing permissions and\n",
       "// limitations under the License.\n",
       "\n",
       "/**\n",
       " * @fileoverview Helpers for google.colab Python module.\n",
       " */\n",
       "(function(scope) {\n",
       "function span(text, styleAttributes = {}) {\n",
       "  const element = document.createElement('span');\n",
       "  element.textContent = text;\n",
       "  for (const key of Object.keys(styleAttributes)) {\n",
       "    element.style[key] = styleAttributes[key];\n",
       "  }\n",
       "  return element;\n",
       "}\n",
       "\n",
       "// Max number of bytes which will be uploaded at a time.\n",
       "const MAX_PAYLOAD_SIZE = 100 * 1024;\n",
       "\n",
       "function _uploadFiles(inputId, outputId) {\n",
       "  const steps = uploadFilesStep(inputId, outputId);\n",
       "  const outputElement = document.getElementById(outputId);\n",
       "  // Cache steps on the outputElement to make it available for the next call\n",
       "  // to uploadFilesContinue from Python.\n",
       "  outputElement.steps = steps;\n",
       "\n",
       "  return _uploadFilesContinue(outputId);\n",
       "}\n",
       "\n",
       "// This is roughly an async generator (not supported in the browser yet),\n",
       "// where there are multiple asynchronous steps and the Python side is going\n",
       "// to poll for completion of each step.\n",
       "// This uses a Promise to block the python side on completion of each step,\n",
       "// then passes the result of the previous step as the input to the next step.\n",
       "function _uploadFilesContinue(outputId) {\n",
       "  const outputElement = document.getElementById(outputId);\n",
       "  const steps = outputElement.steps;\n",
       "\n",
       "  const next = steps.next(outputElement.lastPromiseValue);\n",
       "  return Promise.resolve(next.value.promise).then((value) => {\n",
       "    // Cache the last promise value to make it available to the next\n",
       "    // step of the generator.\n",
       "    outputElement.lastPromiseValue = value;\n",
       "    return next.value.response;\n",
       "  });\n",
       "}\n",
       "\n",
       "/**\n",
       " * Generator function which is called between each async step of the upload\n",
       " * process.\n",
       " * @param {string} inputId Element ID of the input file picker element.\n",
       " * @param {string} outputId Element ID of the output display.\n",
       " * @return {!Iterable<!Object>} Iterable of next steps.\n",
       " */\n",
       "function* uploadFilesStep(inputId, outputId) {\n",
       "  const inputElement = document.getElementById(inputId);\n",
       "  inputElement.disabled = false;\n",
       "\n",
       "  const outputElement = document.getElementById(outputId);\n",
       "  outputElement.innerHTML = '';\n",
       "\n",
       "  const pickedPromise = new Promise((resolve) => {\n",
       "    inputElement.addEventListener('change', (e) => {\n",
       "      resolve(e.target.files);\n",
       "    });\n",
       "  });\n",
       "\n",
       "  const cancel = document.createElement('button');\n",
       "  inputElement.parentElement.appendChild(cancel);\n",
       "  cancel.textContent = 'Cancel upload';\n",
       "  const cancelPromise = new Promise((resolve) => {\n",
       "    cancel.onclick = () => {\n",
       "      resolve(null);\n",
       "    };\n",
       "  });\n",
       "\n",
       "  // Wait for the user to pick the files.\n",
       "  const files = yield {\n",
       "    promise: Promise.race([pickedPromise, cancelPromise]),\n",
       "    response: {\n",
       "      action: 'starting',\n",
       "    }\n",
       "  };\n",
       "\n",
       "  cancel.remove();\n",
       "\n",
       "  // Disable the input element since further picks are not allowed.\n",
       "  inputElement.disabled = true;\n",
       "\n",
       "  if (!files) {\n",
       "    return {\n",
       "      response: {\n",
       "        action: 'complete',\n",
       "      }\n",
       "    };\n",
       "  }\n",
       "\n",
       "  for (const file of files) {\n",
       "    const li = document.createElement('li');\n",
       "    li.append(span(file.name, {fontWeight: 'bold'}));\n",
       "    li.append(span(\n",
       "        `(${file.type || 'n/a'}) - ${file.size} bytes, ` +\n",
       "        `last modified: ${\n",
       "            file.lastModifiedDate ? file.lastModifiedDate.toLocaleDateString() :\n",
       "                                    'n/a'} - `));\n",
       "    const percent = span('0% done');\n",
       "    li.appendChild(percent);\n",
       "\n",
       "    outputElement.appendChild(li);\n",
       "\n",
       "    const fileDataPromise = new Promise((resolve) => {\n",
       "      const reader = new FileReader();\n",
       "      reader.onload = (e) => {\n",
       "        resolve(e.target.result);\n",
       "      };\n",
       "      reader.readAsArrayBuffer(file);\n",
       "    });\n",
       "    // Wait for the data to be ready.\n",
       "    let fileData = yield {\n",
       "      promise: fileDataPromise,\n",
       "      response: {\n",
       "        action: 'continue',\n",
       "      }\n",
       "    };\n",
       "\n",
       "    // Use a chunked sending to avoid message size limits. See b/62115660.\n",
       "    let position = 0;\n",
       "    do {\n",
       "      const length = Math.min(fileData.byteLength - position, MAX_PAYLOAD_SIZE);\n",
       "      const chunk = new Uint8Array(fileData, position, length);\n",
       "      position += length;\n",
       "\n",
       "      const base64 = btoa(String.fromCharCode.apply(null, chunk));\n",
       "      yield {\n",
       "        response: {\n",
       "          action: 'append',\n",
       "          file: file.name,\n",
       "          data: base64,\n",
       "        },\n",
       "      };\n",
       "\n",
       "      let percentDone = fileData.byteLength === 0 ?\n",
       "          100 :\n",
       "          Math.round((position / fileData.byteLength) * 100);\n",
       "      percent.textContent = `${percentDone}% done`;\n",
       "\n",
       "    } while (position < fileData.byteLength);\n",
       "  }\n",
       "\n",
       "  // All done.\n",
       "  yield {\n",
       "    response: {\n",
       "      action: 'complete',\n",
       "    }\n",
       "  };\n",
       "}\n",
       "\n",
       "scope.google = scope.google || {};\n",
       "scope.google.colab = scope.google.colab || {};\n",
       "scope.google.colab._files = {\n",
       "  _uploadFiles,\n",
       "  _uploadFilesContinue,\n",
       "};\n",
       "})(self);\n",
       "</script> "
      ]
     },
     "metadata": {}
    },
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "Saving synthetic_data.csv to synthetic_data.csv\n"
     ]
    },
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "    nam      ageT0    age_end sex  Lesion  age_1stCAD  age_1stARR  \\\n",
       "0  1884  71.661537  88.588678   M  Severe   73.465811         NaN   \n",
       "1  2153  77.250255  88.897479   M   Valve   67.239079         NaN   \n",
       "2  1468  75.216748  85.942570   M   other         NaN   85.942570   \n",
       "3  6820  71.886454  82.930905   M   Shunt   52.299682         NaN   \n",
       "4  3403  66.957054  69.821617   F   other   62.555418         NaN   \n",
       "5  6654  76.571826  82.973190   M   Valve         NaN   56.247117   \n",
       "6  2521  68.743398  76.184429   F   Valve         NaN         NaN   \n",
       "7  2943  84.097805  96.972426   F   Shunt         NaN   72.049700   \n",
       "8  1948   7.663564  79.727097   F   Valve         NaN         NaN   \n",
       "9  8583  67.100556  74.339450   F  Severe   56.657575   74.339450   \n",
       "\n",
       "   age_1stVARR  age_1stPH  age_1stDIAB  ...  age_allSX5  age_allSX6  \\\n",
       "0          NaN  77.994420          NaN  ...         NaN         NaN   \n",
       "1          NaN  71.783267          NaN  ...         NaN         NaN   \n",
       "2          NaN  55.039581          NaN  ...         NaN         NaN   \n",
       "3          NaN        NaN          NaN  ...         NaN         NaN   \n",
       "4          NaN        NaN          NaN  ...         NaN         NaN   \n",
       "5          NaN        NaN    63.289523  ...         NaN         NaN   \n",
       "6          NaN        NaN    76.184429  ...         NaN         NaN   \n",
       "7          NaN        NaN          NaN  ...         NaN         NaN   \n",
       "8    63.653177        NaN          NaN  ...         NaN         NaN   \n",
       "9          NaN  46.798190          NaN  ...         NaN         NaN   \n",
       "\n",
       "   age_allSX7  age_allSX8  age_allSX9  age_allSX10  age_allSX11  age_allSX12  \\\n",
       "0         NaN         NaN         NaN          NaN          NaN          NaN   \n",
       "1         NaN         NaN         NaN          NaN          NaN          NaN   \n",
       "2         NaN         NaN         NaN          NaN          NaN          NaN   \n",
       "3         NaN         NaN         NaN          NaN          NaN          NaN   \n",
       "4         NaN         NaN         NaN          NaN          NaN          NaN   \n",
       "5         NaN         NaN         NaN          NaN          NaN          NaN   \n",
       "6         NaN         NaN         NaN          NaN          NaN          NaN   \n",
       "7         NaN         NaN         NaN          NaN          NaN          NaN   \n",
       "8         NaN         NaN         NaN          NaN          NaN          NaN   \n",
       "9         NaN         NaN         NaN          NaN          NaN          NaN   \n",
       "\n",
       "   age_allSX13  HF  \n",
       "0          NaN   1  \n",
       "1          NaN   1  \n",
       "2          NaN   1  \n",
       "3          NaN   1  \n",
       "4          NaN   1  \n",
       "5          NaN   1  \n",
       "6          NaN   1  \n",
       "7          NaN   1  \n",
       "8          NaN   1  \n",
       "9          NaN   1  \n",
       "\n",
       "[10 rows x 30 columns]"
      ],
      "text/html": [
       "\n",
       "  <div id=\"df-6f108fab-ac09-4c6d-91b1-cc59bba2c890\">\n",
       "    <div class=\"colab-df-container\">\n",
       "      <div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>nam</th>\n",
       "      <th>ageT0</th>\n",
       "      <th>age_end</th>\n",
       "      <th>sex</th>\n",
       "      <th>Lesion</th>\n",
       "      <th>age_1stCAD</th>\n",
       "      <th>age_1stARR</th>\n",
       "      <th>age_1stVARR</th>\n",
       "      <th>age_1stPH</th>\n",
       "      <th>age_1stDIAB</th>\n",
       "      <th>...</th>\n",
       "      <th>age_allSX5</th>\n",
       "      <th>age_allSX6</th>\n",
       "      <th>age_allSX7</th>\n",
       "      <th>age_allSX8</th>\n",
       "      <th>age_allSX9</th>\n",
       "      <th>age_allSX10</th>\n",
       "      <th>age_allSX11</th>\n",
       "      <th>age_allSX12</th>\n",
       "      <th>age_allSX13</th>\n",
       "      <th>HF</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1884</td>\n",
       "      <td>71.661537</td>\n",
       "      <td>88.588678</td>\n",
       "      <td>M</td>\n",
       "      <td>Severe</td>\n",
       "      <td>73.465811</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>77.994420</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2153</td>\n",
       "      <td>77.250255</td>\n",
       "      <td>88.897479</td>\n",
       "      <td>M</td>\n",
       "      <td>Valve</td>\n",
       "      <td>67.239079</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>71.783267</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1468</td>\n",
       "      <td>75.216748</td>\n",
       "      <td>85.942570</td>\n",
       "      <td>M</td>\n",
       "      <td>other</td>\n",
       "      <td>NaN</td>\n",
       "      <td>85.942570</td>\n",
       "      <td>NaN</td>\n",
       "      <td>55.039581</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>6820</td>\n",
       "      <td>71.886454</td>\n",
       "      <td>82.930905</td>\n",
       "      <td>M</td>\n",
       "      <td>Shunt</td>\n",
       "      <td>52.299682</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>3403</td>\n",
       "      <td>66.957054</td>\n",
       "      <td>69.821617</td>\n",
       "      <td>F</td>\n",
       "      <td>other</td>\n",
       "      <td>62.555418</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>6654</td>\n",
       "      <td>76.571826</td>\n",
       "      <td>82.973190</td>\n",
       "      <td>M</td>\n",
       "      <td>Valve</td>\n",
       "      <td>NaN</td>\n",
       "      <td>56.247117</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>63.289523</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>2521</td>\n",
       "      <td>68.743398</td>\n",
       "      <td>76.184429</td>\n",
       "      <td>F</td>\n",
       "      <td>Valve</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>76.184429</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>2943</td>\n",
       "      <td>84.097805</td>\n",
       "      <td>96.972426</td>\n",
       "      <td>F</td>\n",
       "      <td>Shunt</td>\n",
       "      <td>NaN</td>\n",
       "      <td>72.049700</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>1948</td>\n",
       "      <td>7.663564</td>\n",
       "      <td>79.727097</td>\n",
       "      <td>F</td>\n",
       "      <td>Valve</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>63.653177</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>8583</td>\n",
       "      <td>67.100556</td>\n",
       "      <td>74.339450</td>\n",
       "      <td>F</td>\n",
       "      <td>Severe</td>\n",
       "      <td>56.657575</td>\n",
       "      <td>74.339450</td>\n",
       "      <td>NaN</td>\n",
       "      <td>46.798190</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>10 rows × 30 columns</p>\n",
       "</div>\n",
       "      <button class=\"colab-df-convert\" onclick=\"convertToInteractive('df-6f108fab-ac09-4c6d-91b1-cc59bba2c890')\"\n",
       "              title=\"Convert this dataframe to an interactive table.\"\n",
       "              style=\"display:none;\">\n",
       "        \n",
       "  <svg xmlns=\"http://www.w3.org/2000/svg\" height=\"24px\"viewBox=\"0 0 24 24\"\n",
       "       width=\"24px\">\n",
       "    <path d=\"M0 0h24v24H0V0z\" fill=\"none\"/>\n",
       "    <path d=\"M18.56 5.44l.94 2.06.94-2.06 2.06-.94-2.06-.94-.94-2.06-.94 2.06-2.06.94zm-11 1L8.5 8.5l.94-2.06 2.06-.94-2.06-.94L8.5 2.5l-.94 2.06-2.06.94zm10 10l.94 2.06.94-2.06 2.06-.94-2.06-.94-.94-2.06-.94 2.06-2.06.94z\"/><path d=\"M17.41 7.96l-1.37-1.37c-.4-.4-.92-.59-1.43-.59-.52 0-1.04.2-1.43.59L10.3 9.45l-7.72 7.72c-.78.78-.78 2.05 0 2.83L4 21.41c.39.39.9.59 1.41.59.51 0 1.02-.2 1.41-.59l7.78-7.78 2.81-2.81c.8-.78.8-2.07 0-2.86zM5.41 20L4 18.59l7.72-7.72 1.47 1.35L5.41 20z\"/>\n",
       "  </svg>\n",
       "      </button>\n",
       "      \n",
       "  <style>\n",
       "    .colab-df-container {\n",
       "      display:flex;\n",
       "      flex-wrap:wrap;\n",
       "      gap: 12px;\n",
       "    }\n",
       "\n",
       "    .colab-df-convert {\n",
       "      background-color: #E8F0FE;\n",
       "      border: none;\n",
       "      border-radius: 50%;\n",
       "      cursor: pointer;\n",
       "      display: none;\n",
       "      fill: #1967D2;\n",
       "      height: 32px;\n",
       "      padding: 0 0 0 0;\n",
       "      width: 32px;\n",
       "    }\n",
       "\n",
       "    .colab-df-convert:hover {\n",
       "      background-color: #E2EBFA;\n",
       "      box-shadow: 0px 1px 2px rgba(60, 64, 67, 0.3), 0px 1px 3px 1px rgba(60, 64, 67, 0.15);\n",
       "      fill: #174EA6;\n",
       "    }\n",
       "\n",
       "    [theme=dark] .colab-df-convert {\n",
       "      background-color: #3B4455;\n",
       "      fill: #D2E3FC;\n",
       "    }\n",
       "\n",
       "    [theme=dark] .colab-df-convert:hover {\n",
       "      background-color: #434B5C;\n",
       "      box-shadow: 0px 1px 3px 1px rgba(0, 0, 0, 0.15);\n",
       "      filter: drop-shadow(0px 1px 2px rgba(0, 0, 0, 0.3));\n",
       "      fill: #FFFFFF;\n",
       "    }\n",
       "  </style>\n",
       "\n",
       "      <script>\n",
       "        const buttonEl =\n",
       "          document.querySelector('#df-6f108fab-ac09-4c6d-91b1-cc59bba2c890 button.colab-df-convert');\n",
       "        buttonEl.style.display =\n",
       "          google.colab.kernel.accessAllowed ? 'block' : 'none';\n",
       "\n",
       "        async function convertToInteractive(key) {\n",
       "          const element = document.querySelector('#df-6f108fab-ac09-4c6d-91b1-cc59bba2c890');\n",
       "          const dataTable =\n",
       "            await google.colab.kernel.invokeFunction('convertToInteractive',\n",
       "                                                     [key], {});\n",
       "          if (!dataTable) return;\n",
       "\n",
       "          const docLinkHtml = 'Like what you see? Visit the ' +\n",
       "            '<a target=\"_blank\" href=https://colab.research.google.com/notebooks/data_table.ipynb>data table notebook</a>'\n",
       "            + ' to learn more about interactive tables.';\n",
       "          element.innerHTML = '';\n",
       "          dataTable['output_type'] = 'display_data';\n",
       "          await google.colab.output.renderOutput(dataTable, element);\n",
       "          const docLink = document.createElement('div');\n",
       "          docLink.innerHTML = docLinkHtml;\n",
       "          element.appendChild(docLink);\n",
       "        }\n",
       "      </script>\n",
       "    </div>\n",
       "  </div>\n",
       "  "
      ]
     },
     "metadata": {},
     "execution_count": 2
    }
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Preprocessing"
   ],
   "metadata": {
    "id": "Ff0VIBrcEghb"
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "![picture](https://drive.google.com/uc?id=1V43wCFYm5q9LNZIWltohuUirHDyx1VSC)\n",
    "\n",
    "The goal of this section is to transform the event-based nature of all patient data (left picture) into a discretized sequence which is a sequence such that each time-step represents patient data within the given 6-month interval (patient data on right). Additionally, we will produced its associated HF labels needed to be predicted (patient hf on right).  \n"
   ],
   "metadata": {
    "id": "fsllvTqs7hmJ"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "####### Define the Column which will be used in Sequence #######\n",
    "surgeries_cols = [f'age_allSX{i}' for i in range(1, 14)]\n",
    "\n",
    "birth_cols = ['sex', 'Lesion']\n",
    "\n",
    "diagnoses_cols = [\n",
    "    'CAD',\n",
    "    'ARR',\n",
    "    'VARR',\n",
    "    'PH',\n",
    "    'DIAB',\n",
    "    'Stroke',\n",
    "    'CLD',\n",
    "    'CKD',\n",
    "]\n",
    "\n",
    "temporal_cols = [\n",
    "    'AMI',\n",
    "    'IE',\n",
    "    'Sepsis'\n",
    "]\n",
    "\n",
    "ignored_cols = [\n",
    "    'deces',\n",
    "    'birth',\n",
    "    'ADMDATE',\n",
    "    'timezero',\n",
    "    'age_congSX1',\n",
    "    'age_congSX2',\n",
    "    'age_congSX3',\n",
    "    'age_congSX4',\n",
    "    'age_congSX5',\n",
    "    'age_congSX6',\n",
    "    'age_congSX7',\n",
    "    'age_congSX8',\n",
    "]\n",
    "\n",
    "diagnoses_time_col = list(map(lambda st: f'age_1st{st}', diagnoses_cols))\n",
    "temporal_time_col = list(map(lambda st: f'age_1st{st}', temporal_cols))"
   ],
   "metadata": {
    "id": "2mSErHzoEkG5"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "def preprocess(df):\n",
    "  #df.drop(columns=ignored_cols, inplace=True)\n",
    "  \n",
    "  \"\"\"  SURGERY \"\"\"\n",
    "  #create surgery record for each patient -> output is a record with two columns(nam and age of surgery)\n",
    "  #row name and all surgeries\n",
    "  surgeries = df[['nam'] + surgeries_cols].copy()\n",
    "  #remove duplicates\n",
    "  surgeries = surgeries.drop_duplicates()\n",
    "  #transform into nam+variable+value basically verticalize \n",
    "  surgery_record = surgeries.melt(id_vars=['nam'])\n",
    "  #remove NaN\n",
    "  surgery_record = surgery_record.dropna().sort_values('nam')\n",
    "  #reset_index and drop variable( age_allSX)\n",
    "  surgery_record = surgery_record.reset_index(drop=True).drop(columns=['variable'])\n",
    "  # Discretize the age to be in a 6-month bucket ->>> round the value\n",
    "  surgery_record.value = surgery_record.value.apply(lambda age: round(age * 2) / 2)\n",
    "  # result is a 2 column df with nam and value showing all surgeries ages based on nam\n",
    "\n",
    "  \"\"\" Diagnosis \"\"\"\n",
    "  #row will all nam + diagnosis (8 columns) \n",
    "  diagnoses = df[['nam'] + diagnoses_time_col].drop_duplicates()\n",
    "  #transform vertically to 3 columns (nam+vairable(8)+value(age))\n",
    "  diagnoses_record = diagnoses.melt(id_vars=['nam']).sort_values(['nam', 'variable'])\n",
    "  #change NaN to inf and reset index\n",
    "  diagnoses_record = diagnoses_record.reset_index(drop=True).fillna(np.inf)\n",
    "  #set nam as index\n",
    "  diagnoses_record = diagnoses_record.set_index('nam')\n",
    "  # result is a 3 column df with nam(index), variable (8 per person), value(age)\n",
    "  # ['CAD','ARR','VARR','PH','DIAB','Stroke','CLD','CKD']\n",
    "  \n",
    "  \"\"\" temporal features \"\"\"\n",
    "  # 4 columns name+ 3 temporal\n",
    "  temporal = df[['nam'] + temporal_time_col].drop_duplicates()\n",
    "  #transform to vertical \n",
    "  temporal_record = temporal.melt(id_vars=['nam']).sort_values(['nam', 'variable'])\n",
    "  #Nan -> inf\n",
    "  temporal_record = temporal_record.reset_index(drop=True).fillna(np.inf)\n",
    "  # set name as index\n",
    "  temporal_record = temporal_record.set_index('nam')\n",
    "  # round value to age\n",
    "  temporal_record.value = np.round(temporal_record.value * 2.) / 2.\n",
    "  # result is a 3 column df with nam(index), variable (3 per person), value(age)\n",
    "  # ['AMI','IE','Sepsis']\n",
    "\n",
    "  \"\"\" birth record \"\"\"\n",
    "  # 3 column -> nam + (sexe,lesion)\n",
    "  birth_record = df[['nam'] + birth_cols].drop_duplicates()\n",
    "  #sort by nam\n",
    "  birth_record = birth_record.reset_index(drop=True).sort_values(['nam'])\n",
    "  # change sexe(M or F)-> sex_F and sex_M binary variables and Lesion into 4 binary (severe,shunt,valve,other)\n",
    "  birth_record = pd.get_dummies(birth_record, birth_cols)\n",
    "  #let nam be the index\n",
    "  birth_record = birth_record.set_index('nam')\n",
    "  # result is a 7 columns df with nam(index), 6 variable (2 sexe + 4 lesion )\n",
    "  # ['M','F', 'Lesion1','Lesion2','Lesion3','Lesion4']\n",
    "\n",
    "  \"\"\" HF record \"\"\"\n",
    "  # 4 columns -> [nam,\tHF,\tageT0,\tage_end]\n",
    "  hf_record = df[['nam', 'HF', 'ageT0', 'age_end']].copy()\n",
    "  # fancy way of saying dropna\n",
    "  non_hf_duration_threshold = 10\n",
    "  only_select_hf_positive = True\n",
    "  if only_select_hf_positive:\n",
    "      hf_record = hf_record.dropna()\n",
    "  else:\n",
    "      # Now, if patient ageT0 is missing, we automatically set it to a fixed # of years\n",
    "      # before the end of their studies\n",
    "      hf_record.ageT0 = hf_record.ageT0.fillna(\n",
    "          np.maximum(hf_record.age_end - non_hf_duration_threshold, 0)\n",
    "      )\n",
    "  # Set index to be nam\n",
    "  hf_record = hf_record.set_index('nam')\n",
    "  # Round HF age to nearest 6-month interval (0.5 year), and round up the end age\n",
    "  hf_record.ageT0 = np.round(hf_record.ageT0 * 2.) / 2.\n",
    "  hf_record.age_end = np.ceil(hf_record.age_end * 2.) / 2.\n",
    "  #extract just unique ids\n",
    "  all_patient_ids = hf_record.index.unique()\n",
    "  remaining_patient_ids = []\n",
    "  cms = []\n",
    "  hf_labels = []\n",
    "  cm_labels = []\n",
    "\n",
    "  #For each patient create the time-series sequence\n",
    "  for patient_id in tqdm(all_patient_ids):\n",
    "      #one row based all based on the same patient_id\n",
    "      patient_birth_rec = birth_record.loc[patient_id]\n",
    "      patient_diag_rec = diagnoses_record.loc[patient_id]\n",
    "      patient_hf_rec = hf_record.loc[patient_id]\n",
    "      patient_temp_rec = temporal_record.loc[patient_id]\n",
    "      #age \n",
    "      age_start = 40.\n",
    "      age_end = patient_hf_rec.age_end.max()\n",
    "      \n",
    "      if age_end < 40.:\n",
    "          continue\n",
    "    # pass through the function below which merges each subdataframe together with respect to time intervals.\n",
    "      patient_cms, patient_hf_labels = compute_patient_cms_labels(\n",
    "          patient_id,\n",
    "          patient_hf_rec, \n",
    "          patient_diag_rec, \n",
    "          patient_birth_rec, \n",
    "          patient_temp_rec,\n",
    "          surgery_record, \n",
    "          age_start, \n",
    "          age_end\n",
    "      )\n",
    "      \n",
    "      cms.append(patient_cms)\n",
    "      hf_labels.append(patient_hf_labels)\n",
    "      remaining_patient_ids.append(patient_id)\n",
    "  \n",
    "  \"\"\" MERGE \"\"\"\n",
    "  return remaining_patient_ids, cms, hf_labels"
   ],
   "metadata": {
    "id": "Jnp7WBCKEvuZ"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "def compute_patient_cms_labels(patient_id, patient_hf_rec, patient_diag_rec, patient_birth_rec,\n",
    "                               patient_temp_rec, surgery_record, age_start, age_end, \n",
    "                               min_age=40.):\n",
    "    patient_cms = []\n",
    "    patient_hf_labels = []\n",
    "    patient_cm_labels = []\n",
    "    \n",
    "    age_start = max(min_age, age_start)\n",
    "\n",
    "    #here we tranform age at health event into a binary values depending on if the event occured at the age interval\n",
    "    for age in np.arange(age_start, age_end+1, 0.5):        \n",
    "        # Compute the birth records and permanent CMs\n",
    "        single_birth_cm = patient_birth_rec.values\n",
    "        single_diag_cm = np.array(patient_diag_rec.value < age, dtype=int)\n",
    "        \n",
    "        # Compute temporal CMs and surgery\n",
    "        single_temp_cm = np.array(patient_temp_rec.value == age, dtype=int)\n",
    "        single_surg_label = np.any((surgery_record.nam == patient_id) & (surgery_record.value == age)).astype(int)\n",
    "\n",
    "        # Compute current and next HF\n",
    "        hf_current = np.any((patient_hf_rec.ageT0 == age) & (patient_hf_rec.HF == 1)).astype(int)\n",
    "        hf_next = np.any((patient_hf_rec.ageT0 == age + 0.5) & (patient_hf_rec.HF == 1)).astype(int)\n",
    "        \n",
    "        \n",
    "        \n",
    "        # Concatenate the CMs\n",
    "        single_cms = np.concatenate([\n",
    "            single_birth_cm, single_diag_cm, single_temp_cm,\n",
    "            np.array([single_surg_label, hf_current, age])\n",
    "        ])\n",
    "        patient_cms.append(single_cms)\n",
    "        patient_hf_labels.append(hf_next)\n",
    "\n",
    "\n",
    "    patient_cms = np.stack(patient_cms)\n",
    "    patient_hf_labels = np.array(patient_hf_labels)\n",
    "\n",
    "    return patient_cms, patient_hf_labels"
   ],
   "metadata": {
    "id": "dlrDI9owFTl_"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# For each patient we turn event based data into sequencial data (cms = input, hf_labels = output)\n",
    "ids,cms,hf_labels = preprocess(df)\n",
    "\n",
    "\"\"\"\n",
    "Note this step will take roughly 25 minutes to run (you only need to do this once)  \n",
    "You can run the next cell and save the pickled files by clicking on this in google colab local files on the right and dowloading them\n",
    "Then upload them and repickle them using the cell after the next to skip the preprocess steps.\n",
    "\"\"\""
   ],
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 105
    },
    "id": "uEssYVkDFzGJ",
    "outputId": "882c467e-45aa-40b5-d4a5-aa89e67793c2"
   },
   "execution_count": null,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stderr",
     "text": [
      "100%|██████████| 8428/8428 [17:32<00:00,  8.01it/s]\n"
     ]
    },
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "'\\nNote this step will take roughly 25 minutes to run (you only need to do this once)  \\nYou can run the next cell and save the pickled files by clicking on this in google colab local files on the right and dowloading them\\nThen upload them and repickle them using the cell after the next to skip the preprocess steps.\\n'"
      ],
      "application/vnd.google.colaboratory.intrinsic+json": {
       "type": "string"
      }
     },
     "metadata": {},
     "execution_count": 6
    }
   ]
  },
  {
   "cell_type": "code",
   "source": [
    "#write files to local colab files\n",
    "output1 = open('cms.pkl', 'wb')\n",
    "pickle.dump(cms,output1)\n",
    "output2 = open('ids.pkl', 'wb')\n",
    "pickle.dump(ids,output2)\n",
    "output3 = open('hf_labels.pkl', 'wb')\n",
    "pickle.dump(hf_labels,output3)"
   ],
   "metadata": {
    "id": "EiQmTCioJldH"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "#read files from local colab files (pickled files must first be uploaded) - Link to \n",
    "infile = open('cms.pkl','rb')\n",
    "cms = pickle.load(infile)\n",
    "infile2 = open('ids.pkl','rb')\n",
    "ids = pickle.load(infile2)\n",
    "infile3 = open('hf_labels.pkl','rb')\n",
    "hf_labels = pickle.load(infile3)"
   ],
   "metadata": {
    "id": "ppbJrKnfKBWr"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# This is what are data now looks like for one patient\n",
    "# if you want to see the whole array run this :  np.set_printoptions(threshold=10000) \n",
    "print(cms[0])\n",
    "\n",
    "print(hf_labels[0])\n",
    "\n",
    "#Notice that the hf_labels are the same the second to last column values  but shift over one time step (that is as we are predicting next timestep) "
   ],
   "metadata": {
    "id": "UcZUd-q_GKZ9",
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "outputId": "296a0293-a65f-455b-b29f-c1cf7054cb4c"
   },
   "execution_count": null,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "[[ 0.   1.   1.  ...  0.   0.  40. ]\n",
      " [ 0.   1.   1.  ...  0.   0.  40.5]\n",
      " [ 0.   1.   1.  ...  0.   0.  41. ]\n",
      " ...\n",
      " [ 0.   1.   1.  ...  0.   0.  88.5]\n",
      " [ 0.   1.   1.  ...  0.   0.  89. ]\n",
      " [ 0.   1.   1.  ...  0.   0.  89.5]]\n",
      "[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0\n",
      " 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0\n",
      " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "\n",
    "# Format Data "
   ],
   "metadata": {
    "id": "j3PODMTQMF1s"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "######### FORMATTING FOR SEQUENCE TO SEQUENCE BASED MODEL #########\n",
    "\n",
    "#Addional Reading:  https://www.analyticsvidhya.com/blog/2020/08/a-simple-introduction-to-sequence-to-sequence-models/\n",
    "\n",
    "\n",
    "def dataset_split(remaining_patient_ids):\n",
    "    train_id, test_id = train_test_split(remaining_patient_ids, test_size=0.25, random_state=2019)\n",
    "    train_id = set(train_id)\n",
    "    test_id = set(test_id)\n",
    "    return train_id, test_id\n",
    "\n",
    "#SPLIT DATA INTO TEST AND TRAIN (1:3 split)\n",
    "train_id, test_id = dataset_split(ids)\n",
    "\n",
    "\n",
    "def train_test_preprocess(cms, labels, train_id, test_id, all_patient_ids):\n",
    "    train_cms = []\n",
    "    train_labels = []\n",
    "    test_cms = []\n",
    "    test_labels = []\n",
    "\n",
    "    for patient_id, cm, target in tqdm(zip(all_patient_ids, cms, labels)):\n",
    "        if patient_id in train_id:\n",
    "            train_cms.append(cm)\n",
    "            train_labels.append(target)\n",
    "        else:\n",
    "            test_cms.append(cm)\n",
    "            test_labels.append(target)\n",
    "    return train_cms, train_labels, test_cms, test_labels\n",
    "\n",
    "#Using the slip ids, we split the data into train and test data and labels\n",
    "train_cms, train_labels, test_cms, test_labels = train_test_preprocess(cms, hf_labels, train_id, test_id, ids)\n",
    "\n"
   ],
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "m2pXWQvH9Atm",
    "outputId": "06f6bafa-7469-426b-a954-d49de75b36b3"
   },
   "execution_count": null,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stderr",
     "text": [
      "8428it [00:00, 967448.32it/s]"
     ]
    },
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "X_train shape: (6321, 185, 20)\n",
      "X_test shape: (2107, 185, 20)\n",
      "y_train shape: (6321, 185, 1)\n",
      "y_test shape: (2107, 185)\n"
     ]
    },
    {
     "output_type": "stream",
     "name": "stderr",
     "text": [
      "\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "PADDING \n",
    "![picture](https://drive.google.com/uc?id=1aSpXlz9pjcqeFdobOTApmyKc1fa5xCVw)\n",
    "\n",
    "The goal of padding is to transform the different lengthed sequence into the same length. Basically, making the data uniforme into order to be able to pass efficiently into our model. We padding with the value -1 so that when the model trains it will ignore this section. "
   ],
   "metadata": {
    "id": "c1sxObv5KNWz"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "#### PADDING #####\n",
    "#https://www.tensorflow.org/api_docs/python/tf/keras/utils/pad_sequences\n",
    "\n",
    "#Padding essentially makes all inputs sequences the same length based on the value 'maxlen' and fills the sequence with 'value=-1' which helps the training ignore that part of the sequence \n",
    "def pad(train_cms, train_labels, test_cms, test_labels):\n",
    "    X_train = pad_sequences(train_cms, value=-1,dtype='float32')\n",
    "    y_train = np.expand_dims(pad_sequences(train_labels), axis=-1)\n",
    "\n",
    "    X_test = pad_sequences(test_cms, value=-1,maxlen=X_train.shape[1],dtype='float32')\n",
    "    y_test = pad_sequences(test_labels,maxlen=X_train.shape[1])\n",
    "    \n",
    "    return X_train, y_train, X_test, y_test\n",
    "X_train, y_train, X_test, y_test = pad(train_cms, train_labels, test_cms, test_labels)\n",
    "\n",
    "#shows the input dimensions\n",
    "print(\"X_train shape: {}\".format(X_train.shape))\n",
    "print(\"X_test shape: {}\".format(X_test.shape))\n",
    "print(\"y_train shape: {}\".format(y_train.shape))\n",
    "print(\"y_test shape: {}\".format(y_test.shape))\n"
   ],
   "metadata": {
    "id": "xm4FfqbRKHYd"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# this normalizes the data using the min-max method but other can be used (note that we don't always need to scale)\n",
    "#https://machinelearningmastery.com/normalize-standardize-time-series-data-python/\n",
    "\"\"\"\n",
    "def min_max(x,max):\n",
    "  return (x + 1)/(max + 1)\n",
    "\n",
    "def find_max(data):\n",
    "  max = 0\n",
    "  for i in range(0,data.shape[0]):\n",
    "    for j in range(0,data.shape[1]):\n",
    "      if(data[i][j][-1]>max):\n",
    "        max = data[i][j][-1]\n",
    "\n",
    "    return max\n",
    "\n",
    "def norm(data):\n",
    "  max = find_max(data)\n",
    "  for i in range(0,data.shape[0]):\n",
    "    for j in range(0,data.shape[1]):\n",
    "      data[i][j][-1] = min_max(data[i][j][-1],max)\n",
    "  \n",
    "  return data\n",
    "\n",
    "X_train = norm(X_train)\n",
    "X_test = norm(X_test)\n",
    "\"\"\""
   ],
   "metadata": {
    "id": "LOBhyUHEToEk",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 105
    },
    "outputId": "3fcff893-a147-48fe-cd61-4e2d7bf14938"
   },
   "execution_count": null,
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "'\\ndef min_max(x,max):\\n  return (x + 1)/(max + 1)\\n\\ndef find_max(data):\\n  max = 0\\n  for i in range(0,data.shape[0]):\\n    for j in range(0,data.shape[1]):\\n      if(data[i][j][-1]>max):\\n        max = data[i][j][-1]\\n\\n    return max\\n\\ndef norm(data):\\n  max = find_max(data)\\n  for i in range(0,data.shape[0]):\\n    for j in range(0,data.shape[1]):\\n      data[i][j][-1] = min_max(data[i][j][-1],max)\\n  \\n  return data\\n\\nX_train = norm(X_train)\\nX_test = norm(X_test)\\n'"
      ],
      "application/vnd.google.colaboratory.intrinsic+json": {
       "type": "string"
      }
     },
     "metadata": {},
     "execution_count": 6
    }
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Models"
   ],
   "metadata": {
    "id": "aTrlX5jNasqx"
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    " FIXED POSITIONAL ENCODER \n",
    "\n",
    "![picture](https://drive.google.com/uc?id=1MVlhQaQgfSWnc8EZDUrD4B3vDT1UVddW)\n",
    "\n",
    "The attention layers used throughout the model see their input as a set of vectors, with no order. \n",
    "Since the model doesn't contain any recurrent or convolutional layers. It needs some way to identify visit order (i.e encode the notion of age/time)\n",
    "\n",
    "We add \"Positional Encoding\" to the input vectors.\n",
    "It uses a set of sines and cosines at different frequencies (across the sequence). By definition nearby elements will have similar position encodings.\n",
    "\n",
    "The formula for calculating the positional encoding (implemented in Python below) is as follows:\n",
    "\n",
    "$$\\Large{PE_{(pos, 2i)} = \\sin(pos / 10000^{2i / d_{model}})} $$\n",
    "$$\\Large{PE_{(pos, 2i+1)} = \\cos(pos / 10000^{2i / d_{model}})} $$\n",
    "\n",
    "If you want to learn more:\n",
    "\n",
    "https://kazemnejad.com/blog/transformer_architecture_positional_encoding/\n",
    "https://machinelearningmastery.com/a-gentle-introduction-to-positional-encoding-in-transformer-models-part-1/#:~:text=Positional%20encoding%20describes%20the%20location,item's%20position%20in%20transformer%20models.\n",
    "\n"
   ],
   "metadata": {
    "id": "kVxBqQUbLwG7"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "\n",
    "def get_angles(pos, i, d_model):\n",
    "  angle_rates = 1 / np.power(10000, (2 * (i//2)) / np.float32(d_model))\n",
    "  return pos * angle_rates\n",
    "\n",
    "def positional_encoding(position, d_model):\n",
    "  angle_rads = get_angles(np.arange(position)[:, np.newaxis],\n",
    "                          np.arange(d_model)[np.newaxis, :],\n",
    "                          d_model)\n",
    "\n",
    "  # apply sin to even indices in the array; 2i\n",
    "  angle_rads[:, 0::2] = np.sin(angle_rads[:, 0::2])\n",
    "\n",
    "  # apply cos to odd indices in the array; 2i+1\n",
    "  angle_rads[:, 1::2] = np.cos(angle_rads[:, 1::2])\n",
    "\n",
    "  pos_encoding = angle_rads[np.newaxis, ...]\n",
    "\n",
    "  return tf.cast(pos_encoding, dtype=tf.float32)"
   ],
   "metadata": {
    "id": "f1BSKwT9WUd6"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "Attention \n",
    "\n",
    "Attention has the benifit that it is fully connected and is able to account for how the each time step in a sequence impacts the others\n",
    "\n",
    "![picture](https://drive.google.com/uc?id=1jV1OW1hOAC0Kpy7HePEKy_gt3HrayyRH)\n",
    "\n",
    "The attention formula creates 3 set of matrices called the query, key and value, which have weights associated to them. \n",
    "\n",
    "The query and the kay get passed into a softmax scaled-dot-product equation. This equation yield an attention matrix of size sequence length by sequence length where each position (att_ij) in the matrix depicts the connection between (visit_i and visit_j). Thus, when we multiply this attention matix by the value matrix which is just the weighted input, we get a context-aware representation on the initial sequence called the CONTEXT VECTOR.\n",
    "\n",
    "![picture](https://drive.google.com/uc?id=1o7au_Py37XNNaU0cmidk0bGzgECvbNQf)\n",
    "\n",
    "It is important to note that the above describes and image is a single head attention mechanism. The multi-head comes in by applying the attention above multiple time in parallel then combining the output through a weighted average.\n",
    "The intuition is that different attention head will pay attention to different patterns and relationships within the sequences thus the multi-head can pick up a deeper understanding than a single one.\n",
    "\n",
    "MASKING:\n",
    "\n",
    "Finally, because we are predicting future time points we need to make sure that we do not consider future time points during the training process. In order to do that we apply a masking matrix which is a square matrix of the length of the sequences where the upper diagonal triangle is set to negative infinity and the rest is zero. We add this matrix to the attention matrix during training to ensure we do not consider any relationships between a time steps and a future time step.  \n",
    "\n",
    "![picture](https://drive.google.com/uc?id=1TqeTe7sb3kzXo9MoeJC4fqBVrL7OguIJ)\n"
   ],
   "metadata": {
    "id": "Z6TnySrwMsId"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "##### IMPLEMENTATION OF MULTI-HEAD ATTENTION MECHANISM ##########\n",
    "\"\"\"\n",
    "To read more about how it works:\n",
    "- https://theaisummer.com/self-attention/\n",
    "- https://data-science-blog.com/blog/2021/04/07/multi-head-attention-mechanism/\n",
    "- https://towardsdatascience.com/transformers-explained-visually-part-3-multi-head-attention-deep-dive-1c1ff1024853\n",
    "\"\"\"\n",
    "\n",
    "class MultiHeadAttention(tf.keras.layers.Layer):\n",
    "    #initialize weights \n",
    "    def __init__(self, d_model, num_heads):\n",
    "        super(MultiHeadAttention, self).__init__()\n",
    "        self.num_heads = num_heads\n",
    "        self.d_model = d_model\n",
    "        assert d_model % self.num_heads == 0\n",
    "        self.depth = d_model // self.num_heads\n",
    "        self.wq = tf.keras.layers.Dense(d_model)\n",
    "        self.wk = tf.keras.layers.Dense(d_model)\n",
    "        self.wv = tf.keras.layers.Dense(d_model)\n",
    "        self.dense = tf.keras.layers.Dense(d_model)\n",
    "\n",
    "    # Masking layer : https://www.tensorflow.org/guide/keras/masking_and_padding#:~:text=Masking%20is%20a%20way%20to,skipped%20when%20processing%20the%20data.\n",
    "    #Basically masking ensures that the training process doesn't cheat. It ensures that the model only considers past data when predict future time points     \n",
    "    def masked(self,size):\n",
    "        mask = 1 - tf.linalg.band_part(tf.ones((size, size)), -1, 0)\n",
    "        mask  = mask * -1e9\n",
    "        return mask  # (seq_len, seq_len)\n",
    "\n",
    "    def scaled_dot_product_attention(self, q, k, v, mask):\n",
    "        \n",
    "        matmul_qk = tf.matmul(q, k, transpose_b=True)  # (..., seq_len_q, seq_len_k)\n",
    "        dk = tf.cast(tf.shape(k)[-1], tf.float32)\n",
    "\n",
    "        # add the mask to the scaled tensor.\n",
    "        if mask is not None:\n",
    "            matmul_qk = matmul_qk+self.masked(mask)\n",
    "        attention  = matmul_qk / tf.math.sqrt(dk)\n",
    "        # softmax is normalized on the last axis (seq_len_k) so that the scores\n",
    "        # add up to 1.\n",
    "        #attention_weights = Activation('softmax')(attention)  \n",
    "        attention_weights = tf.nn.softmax(attention, axis=-1)  # (..., seq_len_q, seq_len_k)\n",
    "        output = tf.matmul(attention_weights, v)  # (..., seq_len_q, depth_v)\n",
    "        return output, attention_weights\n",
    "    \n",
    "    \n",
    "  \n",
    "    def split_heads(self, x, batch_size):\n",
    "        \"\"\"Split the last dimension into (num_heads, depth).\n",
    "        Transpose the result such that the shape is (batch_size, num_heads, seq_len, depth)\n",
    "        \"\"\"\n",
    "        x = tf.reshape(x, (batch_size, -1, self.num_heads, self.depth))\n",
    "        return tf.transpose(x, perm=[0, 2, 1, 3])\n",
    "\n",
    "    def __call__(self, q, k, v,mask):\n",
    "        batch_size = tf.shape(q)[0]\n",
    "        \n",
    "        q = self.wq(q)  # (batch_size, seq_len, d_model)\n",
    "        k = self.wk(k)  # (batch_size, seq_len, d_model)\n",
    "        v = self.wv(v)  # (batch_size, seq_len, d_model)\n",
    "        #print(q.shape)\n",
    "        q = self.split_heads(q, batch_size)  # (batch_size, num_heads, seq_len_q, depth)\n",
    "        k = self.split_heads(k, batch_size)  # (batch_size, num_heads, seq_len_k, depth)\n",
    "        v = self.split_heads(v, batch_size)  # (batch_size, num_heads, seq_len_v, depth)\n",
    "        # scaled_attention.shape == (batch_size, num_heads, seq_len_q, depth)\n",
    "        # attention_weights.shape == (batch_size, num_heads, seq_len_q, seq_len_k)\n",
    "        scaled_attention, attention_weights = self.scaled_dot_product_attention(q, k, v, mask)\n",
    "        scaled_attention = tf.transpose(scaled_attention, perm=[0, 2, 1, 3])  # (batch_size, seq_len_q, num_heads, depth)\n",
    "        concat_attention = tf.reshape(scaled_attention,  (batch_size, -1, self.d_model))  # (batch_size, seq_len_q, d_model)\n",
    "        output = self.dense(concat_attention)  # (batch_size, seq_len_q, d_model)\n",
    "        return output, attention_weights"
   ],
   "metadata": {
    "id": "HaJ1E7nXd6lD"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "#ATTENTION MODEL\n",
    "\n",
    "def build_attention(seq,n_layers,hidden_dim,heads,mask):\n",
    "  input_dim=20  # 20 patient features\n",
    "  input1= Input(shape =(None,input_dim))  #input placer we are pass data into model \n",
    "  x = input1+positional_encoding(seq,input_dim)[:, :seq, :] #apply positional encoder\n",
    "\n",
    "  fn1, att = MultiHeadAttention(d_model=hidden_dim, num_heads=heads)(q=x,k=x,v=x,mask=mask) # apply multihead attention layer\n",
    "  fn1 = LayerNormalization()(fn1) # layernormalization\n",
    "  fn1 = Dropout(0.2)(fn1) #dropout layer\n",
    "  out = Dense(1, activation='sigmoid')(fn1) # dense layer to predict next time step\n",
    "  model = Model(inputs=input1,outputs=out) \n",
    "    \n",
    "  return model\n",
    "\n",
    "  #Try adding to the model\n",
    "\n",
    "    #batch normalization\n",
    "    # convolutional layer are sometimes used\n",
    "\n",
    "    #WE ARE ABLE TO CREATE AN ATTENTION BLOCK BY COMBINING THE ATTENTION LAYER WITH:\n",
    "\n",
    "      #fn1 = LayerNormalization()(fn1) # layernormalization\n",
    "      #fn1 = Dropout(0.2)(fn1) #dropout layer\n",
    "      #feed forward network \n",
    "      #residual connections\n",
    "      #we can create a for loop to have multiple attention blocks\n",
    "\n",
    "    # Some ideas: \n",
    "    #https://towardsdatascience.com/attention-for-time-series-classification-and-forecasting-261723e0006d\n",
    "    #https://towardsdatascience.com/time-series-forecasting-with-deep-learning-and-attention-mechanism-2d001fc871fc\n",
    "    #look up attention based models!!!\n"
   ],
   "metadata": {
    "id": "CwIQq8sTd91t"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "RNN, LSTM, GRU MODELS \n",
    "\n",
    "These models are considered the norm in modelling sequences. \n",
    "\n",
    "RNN: \n",
    "\n",
    "Process sequential data by accounting for the current input and also what has been learned from previous inputs.\n",
    "\n",
    "The advantage is that it accounts for order and previous inputs in a meaningful way. The disadvantages are Processing long sequences can yield memory and computation issues this is often called the vanishing gradient problem.\n",
    "\n",
    "This is how it works: \n",
    "\n",
    "<img src=\"https://drive.google.com/uc?id=1BXg4eSDFpMl5X8ZaKXo_YQ_9AfMge_AZ\" width=\"750\" align=\"center\">\n",
    "\n",
    "Each X_i represents a time step. W is the weights that filter inputs, a is the activation of the hidden layer (a combination of weighted input and the previous hidden state), and b is the output of the hidden layer after it has been transformed, or squashed, using a rectified linear or sigmoid unit. The b_i value represents the hidden state of what was learnt in previous time steps.\n",
    "\n",
    "The vanishing gradient problem mentioned above occurs when the sequence is long (i.e. many X_i). What happens in this case is that the gradient shrink leading to it because difficutl to update the weights and biases at the early section of the sequence (i.e at X_0, X_1). Thus, lets say that a patient had a surgery at age 1, the model would not effectively learn the correlation between that and a Heart Failure at the age of 40. To overcome this, reseachers developed the LSTM and GRU blocks.\n",
    "\n",
    "LSTM:\n",
    "\n",
    "https://colah.github.io/posts/2015-08-Understanding-LSTMs/\n",
    "\n",
    "Above link give a great walk through of how LSTM works.\n",
    "Here's a quick illustration of how the LSTM Block works:\n",
    "\n",
    "\n",
    "![picture](https://drive.google.com/uc?id=1Z9gmE_kKEAVz4LJh18_0dHHYqL4ZtcLF)\n",
    "\n",
    "In the block above, \n",
    "\n",
    "- Red Box (Forget gate):  Purpose of the red box is to learn what to forget from the previous and current states. It is trained to apply a sigmoid activation layer which outputs value near 0 when the input is deemed irrelevant to predictions and 1 if important.\n",
    "\n",
    "- Pink Box (Update gate):  Purpose of the pink box is to determine what new information should be added to the networks long-term memory (cell state), given the previous hidden state and new input data. We call it the update section.\n",
    "\n",
    "- Yellow Box (Output gate): Purpose of the pink box is to decide the new hidden state. To decide this, we will use three things; the newly updated cell state, the previous hidden state and the new input data. Note that this hidden state get trained to output important information specifically about the predictive task where as the new cell state is a representation of the entire sequence at that time step.\n",
    "\n",
    "As each gate decides the relevancy of information that is stored in the cell and hidden states, when back propagating from cell to cell, the passed error is as close to 1 as possible. This ensures that there is no vanishing of gradients. If you understand residual networks it works in a similar way as the new cell state connects the layers inside the cell with information that stabilizes the propagation of error.\n",
    "\n",
    "GRU:\n",
    "\n",
    "The GRU also tries to solve the vanishing gradient problem. GRUs solve this problem through the use of two gates, the update gate and reset gate. These gates decide what information is allowed through to the output and can be trained to retain information from farther back. \n",
    "Here a link to how it works: https://towardsdatascience.com/understanding-gru-networks-2ef37df6c9be\n",
    "\n",
    "Here's a diagram of the difference between LSTM and GRU:\n",
    "\n",
    "![picture](https://drive.google.com/uc?id=1aLIWWmRsmhlbs2aYti5EoDyNpU7UcqZu)\n",
    "\n",
    "The GRU splits the learning into twwo gates, the reset and update gates. \n",
    "- Reset gate: Used to decide how much past information to forget and is responsible for the short-term memory of the network i.e the hidden state.\n",
    "- Update gate:  Used similar to the forget and input gate of an LSTM. It decides what information to throw away and what new information to add. It serves as the long-term memory update.\n",
    "\n",
    "In all, GRUs has fewer operations; therefore, they are faster to train then LSTMs and less complex. \n",
    "There isn’t a clear winner which one is better. Researchers usually try both to determine which one works better for their use case.\n",
    "\n",
    "\n",
    "It is important to mention that as more experiments were performed with the LSTMs and GRUs there is a certain degree of complexity where this network stops being able to learn. Generally it goes to the thousand time steps before it happens which is already pretty good. However, deep understanding is needed sometimes and has motivated a shift to using attention based models such as the Transformer.\n",
    "\n",
    "\n",
    "\n",
    "To read more about the differences between lstm and gru and rnn:\n",
    "\n",
    "https://towardsdatascience.com/illustrated-guide-to-lstms-and-gru-s-a-step-by-step-explanation-44e9eb85bf21\n",
    "https://www.analyticsvidhya.com/blog/2022/01/tutorial-on-rnn-lstm-gru-with-implementation/\n",
    "https://medium.com/analytics-vidhya/rnn-vs-gru-vs-lstm-863b0b7b1573\n",
    "\n",
    "\n",
    "In the models below you can add feed forward networks before the RNN layers and also you can add additional RNN layer using for loops\n"
   ],
   "metadata": {
    "id": "uG3oQ_Ghlfoo"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "#GRU Model\n",
    "\n",
    "def build_gru_model(n_layers = 2, hidden_dim=64):\n",
    "    input1 = Input(shape=(None,20)) #input placer we are pass data into model \n",
    "    mask1 = Masking(-1)(input1) #masking to ensure the training does not look into the future\n",
    "    layers = GRU(hidden_dim, return_sequences=True)(mask1) #apply GRU layer\n",
    "    out = Dense(1, activation = 'sigmoid')(layers) # dense layer to predict next time step\n",
    "    \n",
    "    model = Model(inputs=input1, outputs = out)\n",
    "    return model"
   ],
   "metadata": {
    "id": "aLMD03fFfEni"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "#LSTM Model\n",
    "def build_lstm_model(n_layers = 2, hidden_dim=64):\n",
    "    input1 = Input(shape=(None,20))\n",
    "    mask1 = Masking(-1)(input1) #masking to ensure the training does not look into the future\n",
    "    layers = LSTM(hidden_dim, return_sequences=True)(mask1) #apply LSTM layer\n",
    "    out = Dense(1, activation = 'sigmoid')(layers) # dense layer to predict next time step\n",
    "    \n",
    "    model = Model(inputs=input1, outputs = out)\n",
    "    return model"
   ],
   "metadata": {
    "id": "lQQH8C-NW1-U"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "#RNN Model\n",
    "def build_RNN_model(n_layers = 2, hidden_dim=64):\n",
    "    input1 = Input(shape=(None,20))\n",
    "    mask1 = Masking(-1)(input1) #masking to ensure the training does not look into the future\n",
    "    layers = SimpleRNN(hidden_dim, return_sequences=True)(mask1) #apply RNN layer\n",
    "    out = Dense(1, activation = 'sigmoid')(layers) # dense layer to predict next time step\n",
    "    \n",
    "    model = Model(inputs=input1, outputs = out)\n",
    "    return model"
   ],
   "metadata": {
    "id": "3jW8TpzgWSbS"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "\n",
    "\n",
    "# Training"
   ],
   "metadata": {
    "id": "VKj5IKANavL7"
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Since each patient will only have a few recorded HF hospitalizations throughout their entire life, the dataset label is heavily skewed towards negative examples, which may mislead the model to achieve a high accuracy at the cost of precision. \n",
    "    \n",
    "FOCAL LOSS overcome this limitation by added a weighted value on top of a normal cross entropy loss\n",
    "\n",
    "References: https://arxiv.org/pdf/1708.02002.pdf\n",
    "\n",
    "Learn more: https://medium.com/visionwizard/understanding-focal-loss-a-quick-read-b914422913e7\n"
   ],
   "metadata": {
    "id": "mbsKVsVWxqDd"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "def binary_focal_loss(gamma=2., alpha=.25):\n",
    "\n",
    "    def binary_focal_loss_fixed(y_true, y_pred):\n",
    "        \"\"\"\n",
    "        :param y_true: A tensor of the same shape as `y_pred`\n",
    "        :param y_pred:  A tensor resulting from a sigmoid\n",
    "        :return: Output tensor.\n",
    "        \"\"\"\n",
    "        y_true = tf.cast(y_true, tf.float32)\n",
    "        # Define epsilon so that the back-propagation will not result in NaN for 0 divisor case\n",
    "        epsilon = K.epsilon()\n",
    "        # Add the epsilon to prediction value\n",
    "        # y_pred = y_pred + epsilon\n",
    "        # Clip the prediciton value\n",
    "        y_pred = K.clip(y_pred, epsilon, 1.0 - epsilon)\n",
    "        # Calculate p_t\n",
    "        p_t = tf.where(K.equal(y_true, 1), y_pred, 1 - y_pred)\n",
    "        # Calculate alpha_t\n",
    "        alpha_factor = K.ones_like(y_true) * alpha\n",
    "        alpha_t = tf.where(K.equal(y_true, 1), alpha_factor, 1 - alpha_factor)\n",
    "        # Calculate cross entropy\n",
    "        cross_entropy = -K.log(p_t)\n",
    "        weight = alpha_t * K.pow((1 - p_t), gamma)\n",
    "        # Calculate focal loss\n",
    "        loss = weight * cross_entropy\n",
    "        # Sum the losses in mini_batch\n",
    "        loss = K.mean(K.sum(loss, axis=1))\n",
    "        return loss\n",
    "\n",
    "    return binary_focal_loss_fixed"
   ],
   "metadata": {
    "id": "vsAk_t-ml2U4"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "epoch_num = 10 #play around\n",
    "valid = 0.15 "
   ],
   "metadata": {
    "id": "8cn98XbrmBh2"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "#train GRU model\n",
    "gru = build_gru_model()\n",
    "gru.compile(loss=[binary_focal_loss(alpha=.25, gamma=2)],optimizer=Adam(learning_rate=1e-2),metrics = [tf.keras.metrics.Precision()])\n",
    "#model.summary()\n",
    "train_hist = gru.fit(\n",
    "            x=X_train, \n",
    "            y=y_train, \n",
    "            epochs=epoch_num, \n",
    "            verbose=0,\n",
    "            validation_split=valid\n",
    "        )"
   ],
   "metadata": {
    "id": "3hDACmFefIVt"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "#train LSTM model\n",
    "lstm = build_lstm_model()\n",
    "lstm.compile(loss=[binary_focal_loss(alpha=.25, gamma=2)],optimizer=Adam(learning_rate=1e-2),metrics = [tf.keras.metrics.Precision()])\n",
    "#model.summary()\n",
    "train_hist = lstm.fit(\n",
    "            x=X_train, \n",
    "            y=y_train, \n",
    "            epochs=epoch_num, \n",
    "            verbose=0,\n",
    "            validation_split=valid\n",
    "        )"
   ],
   "metadata": {
    "id": "Za5ZKLKZVknS"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "#train attention model\n",
    "attention = build_attention(X_train.shape[1],1,64,8,X_train.shape[1])\n",
    "attention.compile(loss=[binary_focal_loss(alpha=.25, gamma=2)],optimizer=Adam(learning_rate=1e-2),metrics = [tf.keras.metrics.Precision()])\n",
    "#model.summary()\n",
    "train_hist = attention.fit(\n",
    "            x=X_train, \n",
    "            y=y_train, \n",
    "            epochs=epoch_num, \n",
    "            verbose=0,\n",
    "            validation_split=valid\n",
    "        )"
   ],
   "metadata": {
    "id": "yO2nu0hCf2g3"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "#train RNN model\n",
    "rnn = build_RNN_model()\n",
    "rnn.compile(loss=[binary_focal_loss(alpha=.25, gamma=2)],optimizer=Adam(learning_rate=1e-2),metrics = [tf.keras.metrics.Precision()])\n",
    "#model.summary()\n",
    "train_hist = rnn.fit(\n",
    "            x=X_train, \n",
    "            y=y_train, \n",
    "            epochs=epoch_num, \n",
    "            verbose=0,\n",
    "            validation_split=valid\n",
    "        )"
   ],
   "metadata": {
    "id": "K_5CKYC0WYyw"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "\n",
    "# Evaluate"
   ],
   "metadata": {
    "id": "C5-a_FEtgGae"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# Calculate prediction and return them with the true labels\n",
    "def compute_flat_pred(model, X_test, test_labels):\n",
    "    with tf.device('/cpu:0'):\n",
    "        y_pred = model.predict(X_test)\n",
    "    \n",
    "    y_pred = np.squeeze(y_pred).tolist()\n",
    "    \n",
    "    true_test = []\n",
    "    pred_test = []\n",
    "\n",
    "    for pred, true in zip(y_pred, test_labels):\n",
    "        pred_test.extend(pred[-len(true):])\n",
    "        true_test.extend(true)\n",
    "\n",
    "    true_test = np.array(true_test)\n",
    "    pred_test = np.array(pred_test)\n",
    "    #print(pred_test.dtype)\n",
    "    \n",
    "    \n",
    "    return true_test, pred_test"
   ],
   "metadata": {
    "id": "d5WDv31qyHad"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "METRICS:\n",
    "\n",
    "Compute metrics we want to evaluate (area under the precision recall curve and area under the receiver operating curve)\n",
    "\n",
    "Area Under the Precision Recall Curve (AUPRC):\n",
    "- How is it calculated:  Evaluate the predicted probabilities against the actual binary labels, to find the Precision and Recall at different thresholds. Find the area under the curve\n",
    "- Ideal value:  Much Greater than the Fraction of Positive in testing set \n",
    "  (i.e all Heart Failures/ all Time Steps = 0.013) \n",
    "- Clinical Significance: For a clinical risk prediction model, the AUPRC tells you if the real diagnoses were identified (recall) without the negative diagnoses being classified as positive (precision)\n",
    "\n",
    "\n",
    "Area Under the Receiver Operating Characteristic (AUROC):\n",
    "- How is it calculated:  Evaluate the predicted probabilities against the actual binary labels, to find the TFR and FPR at different thresholds. Find the area under the curve\n",
    "- Ideal value:  >0.8 is good performance with baseline = 0.5\n",
    "- Clinical Significance: For a clinical risk prediction model, the AUROC tells you the probability that a randomly selected patient who experienced an event will have a higher predicted risk score than a randomly selected patient who did not experience an event.\n"
   ],
   "metadata": {
    "id": "zLQOpQ08yKMn"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "def pr_roc(model, test_cms, test_labels):\n",
    "    y_test_flat, y_pred = compute_flat_pred(model, X_test, test_labels)\n",
    "    precision, recall, thresh_pr = precision_recall_curve(y_test_flat, y_pred)\n",
    "    ap = average_precision_score(y_test_flat, y_pred)\n",
    "    fpr, tpr, thresh_roc = roc_curve(y_test_flat, y_pred)\n",
    "    roc = roc_auc_score(y_test_flat, y_pred)\n",
    "    return precision, recall, thresh_pr, ap, fpr, tpr, thresh_roc, roc\n",
    "\n",
    "\n",
    "def plot_roc(fpr, tpr, roc):\n",
    "    plt.figure(figsize=(6, 5))\n",
    "    plt.plot(fpr, tpr, label=f' (AUC={roc:.4f})')\n",
    "    plt.xlabel('False Positive Rate (FPR)')\n",
    "    plt.ylabel('True Positive Rate(TPR)')\n",
    "    plt.legend()\n",
    "\n",
    "def plot_pr(precision, recall, ap):\n",
    "    plt.figure(figsize=(10, 5))\n",
    "    plt.plot(recall, precision, label=f'(AUC={ap:.4f})')\n",
    "    plt.xlabel('Recall')\n",
    "    plt.ylabel('Precision')\n",
    "    plt.legend()\n"
   ],
   "metadata": {
    "id": "908zkdvhjOsy"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "#attention\n",
    "precision, recall, thresh_pr, ap, fpr, tpr, thresh_roc, roc= pr_roc(attention,X_test,y_test)\n",
    "\n",
    "#GRU\n",
    "precision1, recall1, thresh_pr1, ap1, fpr1, tpr1, thresh_roc1, roc1= pr_roc(gru,X_test,y_test)\n",
    "\n",
    "#LSTM\n",
    "precision2, recall2, thresh_pr2, ap2, fpr2, tpr2, thresh_roc2, roc2= pr_roc(lstm,X_test,y_test)\n",
    "\n",
    "#RNN\n",
    "precision3, recall3, thresh_pr3, ap3, fpr3, tpr3, thresh_roc3, roc3= pr_roc(rnn,X_test,y_test)\n",
    "\n",
    "#Precision-Recall Curve\n",
    "plt.figure(figsize=(10, 5))\n",
    "plt.plot(recall, precision, label=f'ATTENTION : (AUC={ap:.4f})')\n",
    "plt.plot(recall1, precision1, label=f'GRU : (AUC={ap1:.4f})')\n",
    "plt.plot(recall2, precision2, label=f'LSTM : (AUC={ap2:.4f})')\n",
    "plt.plot(recall3, precision3, label=f'RNN : (AUC={ap3:.4f})')\n",
    "plt.title('Precision-Recall Curve')\n",
    "plt.xlabel('Recall')\n",
    "plt.ylabel('Precision')\n",
    "plt.legend()\n",
    "\n",
    "\n",
    "#Receiver Operating Characteristic Curve\n",
    "plt.figure(figsize=(6, 5))\n",
    "plt.plot(fpr, tpr, label=f'ATTENTION : (AUC={roc:.4f})')\n",
    "plt.plot(fpr1, tpr1, label=f'GRU : (AUC={roc1:.4f})')\n",
    "plt.plot(fpr2, tpr2, label=f'LSTM : (AUC={roc2:.4f})')\n",
    "plt.plot(fpr3, tpr3, label=f'RNN : (AUC={roc3:.4f})')\n",
    "plt.title('Receiver Operating Characteristic Curve')\n",
    "plt.xlabel('False Positive Rate (FPR)')\n",
    "plt.ylabel('True Positive Rate(TPR)')\n",
    "\n"
   ],
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 700
    },
    "id": "t_9wT_BZgH07",
    "outputId": "3f0dc34b-1ed9-49e3-84bc-154bc3c709ea"
   },
   "execution_count": null,
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'True Positive Rate(TPR)')"
      ]
     },
     "metadata": {},
     "execution_count": 60
    },
    {
     "output_type": "display_data",
     "data": {
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ],
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAAFNCAYAAABIc7ibAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeXhU5d3/8fd3ZrICYUcgYdNEgbAEiIAtImCRRQV9RMTWX22ttrbSVm2pa622VsStah+s9VErLg3iwiKloBZway0Gi8piIYUICQFJWAIJWWbm/v0xk2lCEgjIOAE/r+vKxZz73Oc+33MA+XifM+eYcw4RERER+XJ5Yl2AiIiIyFeRQpiIiIhIDCiEiYiIiMSAQpiIiIhIDCiEiYiIiMSAQpiIiIhIDCiEiciXwsy+ZWavN6Hf42b2yy+jpi+DmeWb2TfCn+80s+djXZOINA8KYSJSExQOmtkBM9tpZs+YWcvjuQ/n3AvOufOa0O9a59xvjue+a5iZM7Oy8HEWmtlDZuaNxr6OhZmlmNnDZrY1XON/wssdYl2biBx/CmEiUuNC51xLYDCQDdx+aAcz833pVR1/A8PHeQ5wGXBVjOsBwMzigb8BmcB4IAU4CygBhh7DeCfD75XISU0hTETqcM4VAn8F+kFk9ug6M9sEbAq3XWBma8xsr5n93cwG1GxvZt3M7FUz22VmJWb2v+H275jZu+HPZma/M7PPzazUzD4xs5r9PWNmd9ca7xozyzOz3Wa2yMy61lrnzOxaM9sUrmW2mVkTjzMPeA/IqjXesRzXaWa2PNxWbGYvmFmboz7x8G2gO3Cxc269cy7onPvcOfcb59ySWsebXqumyLkys1FmVmBmN5nZDuBPZrbBzC6o1d8Xrn9weHl4+Dj3mtlHZjbqGOoWkWOkECYidZhZN2Ai8K9azRcBw4C+ZjYIeBr4AdAe+COwyMwSwpf2FgOfAT2BVGBuA7s5DxgJnA60BqYSmvE5tJYxwMzw+i7hcQ8d7wLgTGBAuN+4Jh5nb+BsIC+8fKzHZeEauwJ9gG7AnU2p4RDfAJY65w4cw7Y1OgPtgB7A94Ec4PJa68cBxc65D80sFfgLcHd4m58Dr5hZxy+wfxE5CgphIlJjgZntBd4F3gLuqbVupnNut3PuIKF/3P/onPuncy7gnJsDVALDCV026wrMcM6VOecqnHPvNrCvaqAV0Bsw59wG51xRA/2+BTztnPvQOVcJ3AKcZWY9a/W51zm31zm3FVhBrZmtRnxoZmXABmAl8Fi4/ZiOyzmX55x7wzlX6ZzbBTxE6FLn0WoPNHQOjkYQ+FW4loPAn4FJZpYcXv9NQsEM4ApgiXNuSXjW7Q0gl1AAF5EvgUKYiNS4yDnXxjnXwzn3o/A/4jW21frcA/hZ+BLW3nBw60YopHQDPnPO+Q+3I+fccuB/gdnA52b2hJmlNNC1K6HZp5rtDhCaMUut1WdHrc/lQEsAM1sXvrn9gJmdXavP4HCfywjN7rX4IsdlZqeY2dzwjf6lwPPAsdxIX0Jotu+L2OWcq6hZCF9y3QBcGA5ikwgFMwgd76WHHO+I41CDiDSRQpiINIWr9Xkb8NtwYKv5SXbO5YTXdW/KTeHOuUedc0OAvoQuS85ooNt2QmEBADNrQWjGqLAJ42c651qGf945ZJ1zzs0D/gHc8QWP6x5C56e/cy6F0AxTk+5LO8SbwLjwMTamHEiutdz5kPWO+mouSU4G1oeDGYSO6blDjreFc+7eY6hdRI6BQpiIHK3/A641s2HhG+xbmNn5ZtYKWEXoktq94fZEM/v6oQOY2Znh7eOAMqCC0KW0Q+UA3zWzLDNLIBR4/umcyz9Ox3IvcI2Zdf4Cx9UKOADsC99n1VCYbIrnCAWjV8yst5l5zKy9md1qZjWXCNcA3zQzr5mNp2mXPecSugfvh/x3FgxCM3YXmtm48HiJ4Zv7046xfhE5SgphInJUnHO5wDWELifuIXRj+3fC6wLAhUA6sBUoIHTZ71AphELPHkKXG0uA+xvY15vAL4FXCIWg04Bpx/FYPgHeJnSv17Ee112ELnHuI3Sj+6vHWEsloZvzPwXeAEoJhb8OwD/D3X4armMvofvlFjRh3CJCM35fA16s1b6N0OzYrcAuQgFwBvp3QeRLY841NHstIiIiItGk/+MRERERiQGFMBEREZEYUAgTERERiQGFMBEREZEYUAgTERERiYEjPlCxuenQoYPr2bNnrMsQEREROaLVq1cXO+cafCfrCRfCevbsSW5ubqzLEBERETkiM/ussXW6HCkiIiISAwphIiIiIjGgECYiIiISAyfcPWEiIvLVUF1dTUFBARUVFbEuReSIEhMTSUtLIy4ursnbKISJiEizVFBQQKtWrejZsydmFutyRBrlnKOkpISCggJ69erV5O10OVJERJqliooK2rdvrwAmzZ6Z0b59+6OetVUIExGRZksBTE4Ux/JnNWohzMyeNrPPzWxtI+vNzB41szwz+9jMBkerFhEREZHmJpozYc8A4w+zfgKQEf75PvCHKNYiIiJyTBYsWICZ8emnnwIwbNgwsrKy6N69Ox07diQrK4v+/fvTpk0bsrKy6Ny5M6mpqWRlZZGVlUVVVRVerzeynJWVxb333gvAqFGjyM7OjuwrNzeXUaNGsWzZskjfli1bcsYZZ5CVlcW3v/1tVq5cyQUXXFCnvgEDBtCnTx/69+/PggULIuu+853vkJqaSmVlJQDFxcUcy1tnrr/+et5+++3IcnFxMXFxcTz++ON1+rVs2bLO8jPPPMP06dMjy88++yz9+vWjf//+DBo0iAceeKDJNcycOZP09HTOOOMMli1b1mCfLVu2MGzYMNLT07nsssuoqqoC4PHHH6d///5kZWUxYsQI1q9fH9nm448/5qyzziIzM5P+/ftHLimuXr2a/v37k56ezk9+8hOccwD8/Oc/Z/ny5U2u+7Ccc1H7AXoCaxtZ90fg8lrL/wa6HGnMIUOGuGj6cOWrbu7PLnRb1n8Q1f2IiMjhrV+/PtYlOOecmzp1qhsxYoS744476rT/6U9/ctddd129/r/61a/c/fffX6etRYsWDY59zjnnuG7durklS5Y455z74IMP3DnnnFOvzwcf/PffpBUrVrjzzz/fOefcmjVr3GmnneY2b97snHNu8+bN7rTTTnMfffSRc865K6+80nXr1s099thjzjnndu3a5Xr06NHEIw8pLi52w4YNq9P22GOPuREjRriRI0ce9jhrn6MlS5a4QYMGucLCQueccxUVFe6JJ55oUg3r1q1zAwYMcBUVFW7z5s3u1FNPdX6/v16/Sy+91OXk5DjnnPvBD34QOe59+/ZF+ixcuNCNGzfOOedcdXW169+/v1uzZk3kWGvGPfPMM90//vEPFwwG3fjx4yO/R/n5+W7s2LEN1tnQn1kg1zWSaWJ5T1gqsK3WckG4rR4z+76Z5ZpZ7q5du6Ja1Jb3FjNg8SbyP343qvsREZHm78CBA7z77rs89dRTzJ07Nyr7mDFjBr/97W+PadsHHniAW2+9NfKNvF69enHLLbdw//33R/pcf/31/O53v8Pv9x/TPl555RXGj697YSsnJ4cHH3yQwsJCCgoKmjTOzJkzeeCBB+jatSsACQkJXHPNNU3aduHChUybNo2EhAR69epFeno6q1atqtPHOcfy5cuZMmUKAFdeeWVkVjAlJSXSr6ysLHL/1uuvv86AAQMYOHAgAO3bt8fr9VJUVERpaSnDhw/HzPj2t78dGatHjx6UlJSwY8eOJtV+OCfEIyqcc08ATwBkZ2e7GJcjIiJfsrteW8f67aXHdcy+XVP41YWZh+2zcOFCxo8fz+mnn0779u1ZvXo1Q4YMOep9HTx4kKysrMjyLbfcwmWXXQbAWWedxfz581mxYgWtWrU6qnHXrVvHz3/+8zpt2dnZzJ49O7LcvXt3RowYwXPPPceFF17Y6FgTJ07kySefjISkGu+9914k2ABs27aNoqIihg4dytSpU3nxxRf52c9+dsRa165d2+i5u//++3nhhRfqtY8cOZJHH32UwsJChg8fHmlPS0ujsLCwTt+SkhLatGmDz+drsM/s2bN56KGHqKqqilxO3LhxI2bGuHHj2LVrF9OmTeMXv/gFhYWFpKWlNbq/wYMH895773HJJZcc8bgPJ5YhrBDoVms5LdwmIiLSLOTk5PDTn/4UgGnTppGTk3NMISwpKYk1a9Y0uv7222/n7rvvZtasWcdc6+HccsstTJ48mfPPP7/RPkuWLGmwvaioiI4dO0aWX3zxRaZOnQqEzslVV1112BDWlG8NzpgxgxkzZhyx3xdx3XXXcd111/HnP/+Zu+++mzlz5uD3+3n33Xf54IMPSE5O5txzz2XIkCG0bt36sGN16tSJ7du3f+GaYhnCFgHTzWwuMAzY55wrimE9IiLSTB1pxioadu/ezfLly/nkk08wMwKBAGbG/ffff9wfnTFmzBhuv/123n///aParm/fvqxevTpyOQ1CN5RnZtY9XxkZGWRlZTFv3ryjri0pKanO869ycnLYsWNHZOZq+/btbNq0iYyMDJKSkqiqqiI+Ph4IncMOHToAkJmZyerVqxkzZky9fRxpJiw1NZVt2/57B1NBQQGpqXXvYGrfvj179+7F7/fj8/ka7AOh4PjDH/4QCM1wjRw5MlLjxIkT+fDDD7niiivqXGY9dKyKigqSkpKOcOaOLJqPqMgB/gGcYWYFZvY9M7vWzK4Nd1kCbAbygP8DfhStWkRERI7Wyy+/zP/7f/+Pzz77jPz8fLZt20avXr145513orK/22+/nfvuu++otvn5z3/OzJkzyc/PByA/P5977rmnwZmp22677ai+jVijT58+5OXlAaHLdwcOHKCwsJD8/Hzy8/O55ZZbyMnJAeCcc87h+eefB0KXYOfNm8fo0aOB0GzcjBkzIvdSVVVV8eSTTwKhmbA1a9bU+3n00UcBmDRpEnPnzqWyspItW7awadMmhg4dWqdOM2P06NG8/PLLAMyZM4fJkycDsGnTpki/v/zlL2RkZAAwbtw4PvnkE8rLy/H7/bz11lv07duXLl26kJKSwvvvv49zjmeffTYyVs156Nev31Gfy0NFLYQ55y53znVxzsU559Kcc0855x53zj0eXu+cc9c5505zzvV3zuVGqxYREZGjlZOTw8UXX1yn7ZJLLokEjqNRc09Yzc/NN99cr8/EiRPrXPZriqysLGbNmsWFF15I7969ufDCC7nvvvvq3H9WIzMzk8GDG38k58SJExu8xHb++eezcuVK4Mjn5JFHHuHVV18lKyuL4cOHc+mllzJy5MjI+NOnT+cb3/hGpJbS0qbd55eZmcnUqVPp27cv48ePZ/bs2Xi93np1z5o1i4ceeoj09HRKSkr43ve+B8D//u//kpmZSVZWFg899BBz5swBoG3bttx4442ceeaZZGVlMXjw4Mgl28cee4yrr76a9PR0TjvtNCZMmACE3mmal5dX59Eix8qcO7Huc8/Ozna5udHLa6/e8z36PPt3dt71A0Zddn3U9iMiIoe3YcMG+vTpE+syBBgxYgSLFy+mTZs2sS4l5ubPn8+HH37Ib37zm3rrGvoza2arnXMNJja9tkhEREQO68EHH2Tr1q2xLqNZ8Pv9Tfo2aFOcEI+oEBERkdgZNmxYrEtoNi699NLjNpZmwkRERERiQCFMREREJAYUwkRERERiQCFMREREJAYUwkRERBqxc+dOvvnNb3LqqacyZMiQyHseAVauXEnr1q3Jysqid+/edd7heOedd9Z7MGrPnj0pLi4+qv0//PDDPPvss5Flv99Px44d6z1n7NCxV65cyQUXXBBZ/utf/0p2djZ9+/Zl0KBBR/Xtvjlz5pCRkUFGRkbk+VqH2r17N2PHjiUjI4OxY8eyZ88eIPTuzQEDBpCVlUV2djbvvvvuEccdP348AwcOJDMzk2uvvZZAIACEHkxb887Hk4VCmIiISAOcc1x00UWMHDmSzZs3s3r1aubOnVvndTZnn302a9as4V//+heLFy/mvffeO2779/v9PP3003zzm9+MtL3xxhucfvrpvPTSSzT1OZ9r165l+vTpPP/886xfv57c3FzS09ObtO3u3bu56667+Oc//8mqVau46667IgGrtnvvvZdzzz2XTZs2ce6553LvvfcCcO655/LRRx+xZs0ann76aa6++uojjjtv3jw++ugj1q5dy65du3jppZcA+PGPfxwZ92ShECYiItKA5cuXEx8fz7XXXhtp69GjBz/+8Y/r9U1KSiIrK4vCwsLjuv/Bgwfj8/33aVI1LxTv3r07//jHP5o0zn333cdtt91G7969AfB6vZF3Jx7JsmXLGDt2LO3ataNt27aMHTuWpUuX1uu3cOFCrrzySgCuvPJKFixYAEDLli0j79ksKyuLfD7cuCkpKUAohFZVVUW26dGjByUlJZHXHp0M9JwwERFp/v56M+z45PiO2bk/TGh8ZmXdunWHfc1PbXv27GHTpk2RV/Qcjauvvpprr7223mtw3nvvPYYMGRJZrqio4M033+SPf/wje/fuJScnh6997WtHHH/t2rWNXn584YUXuP/+++u1p6en8/LLL1NYWEi3bt0i7WlpaQ0GzZ07d9KlSxcAOnfuzM6dOyPr5s+fzy233MLnn3/OX/7yF4Ajjjtu3DhWrVrFhAkTmDJlSqR98ODBvPfee1xyySVHPO4TgWbCREREmuC6665j4MCBnHnmmZG2d955h4EDB5Kamsq4cePo3LkzQGT25lANtT/55JMNvoewqKiozrskFy9ezOjRo0lKSuKSSy5hwYIFkfulGhq3sRpq+9a3vtXgi7NrXoJ9LMyszr4vvvhiPv30UxYsWMAvf/nLJo2xbNkyioqKqKysrHMfWKdOnRp8v+WJSjNhIiLS/B1mxipaMjMzeeWVVyLLs2fPpri4uE5gOvvss1m8eDFbtmxh+PDhTJ06laysLNq3b09RUVGd8fbv339U715MSkqioqIispyTk8O7775Lz549ASgpKWH58uWMHTuW9u3bs2fPHjp06ACE7rmq+ZyZmcnq1asZOHBgvX0caSYsNTU18vJugIKCAkaNGlWv/ymnnEJRURFdunShqKiITp061etTc29dcXFxk8ZNTExk8uTJLFy4kLFjxwKh2cCkpKQGz9eJSDNhIiIiDRgzZgwVFRX84Q9/iLSVl5c32LdXr17cfPPNzJo1CwgFjkWLFrF//34AXn31VQYOHIjX623y/vv06UNeXh4ApaWlvPPOO2zdupX8/Hzy8/OZPXs2OTk5AIwaNYrnnnsOgEAgwPPPP8/o0aMBmDFjBvfccw8bN24EIBgM8vjjjwNHngkbN24cr7/+Onv27GHPnj28/vrrjBs3rl6tkyZNinzDcc6cOUyePBmAvLy8yBcIPvzwQyorK2nfvn2j4x44cCASXv1+P3/5y18i97IBbNy4kX79+jX5HDZ3CmEiIiINMDMWLFjAW2+9Ra9evRg6dChXXnllJGgd6tprr+Xtt98mPz+fAQMGMH36dEaMGEFWVhaPP/44Tz75ZIPbXX311eTm5tZrnzBhAm+//TYQuq9qzJgxJCQkRNZPnjyZ1157jcrKSn75y1+Sl5fHwIEDGTRoEOnp6VxxxRUADBgwgIcffpjLL7+cPn360K9fPzZv3tykc9CuXTt++ctfcuaZZ3LmmWdyxx130K5du3p133zzzbzxxhtkZGTw5ptvRh6h8corr9CvXz+ysrK47rrrePHFFzGzRsctKytj0qRJkcdadOrUKfLFiOrqavLy8hq8dHuisqZ+xbW5yM7Odg39YT1eXr3ne/R59u/svOsHjLrs+qjtR0REDm/Dhg306dMn1mXE1MUXX8x9991HRkZGrEuJufnz5/Phhx/ym9/8JtalNKqhP7Nmtto512By1EyYiIhIM3XvvffWu7fsq8rv9x/VQ2ZPBLoxX0REpJk644wzOOOMM2JdRrNw6aWXxrqE404zYSIiIiIxoBAmIiIiEgMKYSIiIiIxoBAmIiIiEgMKYSIiIo1o2bJlvbZ///vfjBo1iqysLPr06cP3v/99li1bRlZWFllZWbRs2ZIzzjiDrKwsvv3tb7Ny5UrMrM5zwtasWYOZ8cADDxxVPUVFRVxwwQV12q6//npSU1MJBoORtjvvvLPe2D179qS4uBiAHTt2MG3aNE477TSGDBnCxIkTIw9zPZItW7YwbNgw0tPTueyyy6iqqmqw38yZM0lPT+eMM85g2bJlQOiJ90OHDmXgwIFkZmbyq1/9KtLfOcdtt93G6aefTp8+fXj00UfrjPfBBx/g8/kiD5LdtWsX48ePb1LNzZVCmIiIyFH4yU9+wg033MCaNWvYsGEDP/7xjxk3blzkafPZ2dm88MILrFmzhmeffRaAfv36MW/evMgYOTk5Db5G6EgeeughrrnmmshyMBhk/vz5dOvWjbfeeqtJYzjnuPjiixk1ahT/+c9/WL16NTNnzqzz0u3Duemmm7jhhhvIy8ujbdu2PPXUU/X6rF+/nrlz57Ju3TqWLl3Kj370IwKBAAkJCSxfvpyPPvqINWvWsHTpUt5//30AnnnmGbZt28ann37Khg0bmDZtWmS8QCDATTfdxHnnnRdp69ixI126dOG9995rUt3NkUKYiIjIUSgqKiItLS2y3L9//yNu06NHDyoqKti5cyfOOZYuXcqECROOet+vvPJKndmflStXkpmZyQ9/+MPIK4yOZMWKFcTFxUWeRA8wcOBAzj777CNu65xj+fLlTJkyBYArr7ySBQsW1Ou3cOFCpk2bRkJCAr169SI9PZ1Vq1ZhZpHZxerqaqqrqyMv+/7DH/7AHXfcgccTiia13z/5+9//nksuuaTeOykvuugiXnjhhSYdd3Ok54SJiEizN2vVLD7d/elxHbN3u97cNPSmo97uhhtuYMyYMXzta1/jvPPO47vf/W6TXsw9ZcoUXnrpJQYNGsTgwYPrvIKotkWLFpGbm8uvf/3rOu1btmyhbdu2dbbLycnh8ssvZ/Lkydx6661UV1cTFxd32DrWrl3LkCFDGly3f//+RsPYn//8Zzp16kSbNm3w+ULxIS0tjcLCwnp9CwsLGT58eGS5dr9AIMCQIUPIy8vjuuuuY9iwYQD85z//4cUXX2T+/Pl07NiRRx99lIyMDAoLC5k/fz4rVqzggw8+qLOf7Oxsbr/99sMeb3OmmTAREZGj8N3vfpcNGzZw6aWXsnLlSoYPH05lZeURt5s6dSovvfRSJDg1ZtKkSfUCGIRm4Dp27BhZrqqqYsmSJVx00UWkpKQwbNiwyL1XNbNLh2qsvUarVq0afKH3mjVr6Nu37xGPsSm8Xi9r1qyhoKCAVatWsXbtWgAqKytJTEwkNzeXa665hquuugoI3fM2a9asyAxZbZ06dWL79u3Hpa5Y0EyYiIg0e8cyYxVNXbt25aqrruKqq66iX79+h51dqtG5c2fi4uJ44403eOSRR/j73/9+VPtMSkqioqIisrxs2TL27t0buRxaXl5OUlISF1xwAe3bt6/3uqP9+/fTpk0bMjMzIze3H+pIM2F9+vRh7969+P1+fD4fBQUFpKam1uubmprKtm3bIssN9WvTpg2jR49m6dKl9OvXj7S0NP7nf/4HCL0z87vf/S4Aubm5kfvDiouLWbJkCT6fj4suuoiKigqSkpIOe96aM82EiYiIHIWlS5dSXV0NhL5lWFJS0mAQacivf/1rZs2ahdfrPer9nn766eTn50eWc3JyePLJJ8nPzyc/P58tW7bwxhtvUF5ezsiRI1m0aBH79+8H4NVXX2XgwIF4vV7GjBlDZWUlTzzxRGSsjz/+mHfeeeeIM2FmxujRoyMhbs6cOUyePLlerZMmTWLu3LlUVlayZcsWNm3axNChQ9m1axd79+4F4ODBg7zxxhv07t0bCN3ftWLFCgDeeustTj/9dCB0GbbmGKdMmcJjjz3GRRddBMDGjRvp16/fUZ/L5kIzYSIiIo0oLy+vcxP+jTfeSEFBAT/96U9JTEwE4P7776dz585NGu9rX/vaEfs0dk9YixYtOO2008jLy6Nr164sXbqUxx9/vM76ESNG8Nprr3HZZZcxffp0RowYgZnRqVOnyCMyzIz58+dHLvMlJibSs2dPHn744SYdw6xZs5g2bRq33347gwYN4nvf+169ujMzM5k6dSp9+/bF5/Mxe/ZsvF4vRUVFXHnllQQCAYLBIFOnTo08cuPmm2/mW9/6Fr/73e9o2bJlnUd6NGbFihWcf/75Taq7OTLnXKxrOCrZ2dkuNzc3auO/es/36PPs39l51w8Yddn1UduPiIgc3oYNG+jTp0+sy2hW5s+fz+rVq7n77rtjXUqzMHLkSBYuXEjbtm1jXQrQ8J9ZM1vtnMtuqL9mwkRERE4QF198MSUlJbEuo1nYtWsXN954Y7MJYMdC94SJiIicQK6++upYl9AsdOzYMXJv2IlKIUxEREQkBhTCRERERGJAIUxEREQkBhTCRERERGJAIUxERKQRXq+XrKws+vXrx4UXXhh50Gh+fj5mxu9///tI3+nTp/PMM88A8J3vfIfU1NTI64yKi4vp2bPnUe//+uuv5+23344sFxcXExcXV+f5YEDkpdg1nnnmGaZPnx5ZfvbZZ+nXrx/9+/dn0KBBPPDAA02uYebMmaSnp3PGGWdEXot0qC1btjBs2DDS09O57LLLqKqqAuChhx6ib9++DBgwgHPPPZfPPvsMCD3fKysrK/KTmJgYeRH42WefHWnv2rVr5Ob7xYsXc8cddzS57hOBQpiIiEgjkpKSWLNmDWvXrqVdu3bMnj07sq5Tp0488sgjkcBxKK/Xy9NPP33M+y4pKeH9999n5MiRkbaXXnqJ4cOHk5OT0+Rx/vrXv/Lwww/z+uuv88knn/D+++/TunXrJm27fv165s6dy7p161i6dCk/+tGPCAQC9frddNNN3HDDDeTl5dG2bVueeuopAAYNGkRubi4ff/wxU6ZM4Re/+AUAo0ePjjyJf/ny5SQnJ3PeeecB8M4770TWnXXWWZFXGZ1//vm89tprlJeXN/nYmzuFMBERkSY466yzKBJAZdQAACAASURBVCwsjCx37NiRc889lzlz5jTY//rrr+d3v/sdfr//mPb3yiuvMH78+DptOTk5PPjggxQWFlJQUNCkcWbOnMkDDzxA165dAUhISOCaa65p0rYLFy5k2rRpJCQk0KtXL9LT01m1alWdPs45li9fzpQpUwC48sorI7Nao0ePJjk5GYDhw4c3WPPLL7/MhAkTIv1qlJaWsnz58shMmJkxatQoFi9e3KTaTwRRfVirmY0HHgG8wJPOuXsPWd8dmAO0Cfe52Tm3JJo1iYjIiWfHPfdQueHT4zpmQp/edL711ib1DQQC/O1vf4u8oqfGTTfdxIQJE7jqqqvqbdO9e3dGjBjBc889x4UXXtjo2BMnTuTJJ5+MhKQa7733XiTYAGzbto2ioiKGDh3K1KlTefHFF/nZz352xNoP93Lx+++/nxdeeKFe+8iRI3n00UcpLCxk+PDhkfa0tLQ6QRRCM3Zt2rTB5/M12gfgqaeeYsKECfXa586dy4033livfcGCBZx77rmkpKRE2rKzs3nnnXeYOnVqI0d7YolaCDMzLzAbGAsUAB+Y2SLn3Ppa3W4H5jnn/mBmfYElQM9o1SQiInI0Dh48SFZWFoWFhfTp04exY8fWWX/qqacybNgw/vznPze4/S233MLkyZMP+37DJUsannsoKiqiY8eOkeUXX3wxEj6mTZvGVVddddgQZmaNrqsxY8YMZsyYccR+X9Tzzz9Pbm4ub731Vp32oqIiPvnkE8aNG1dvm5ycnHoPpu3UqRPbt2+Paq1fpmjOhA0F8pxzmwHMbC4wGagdwhxQE3FbAyfPmRURkeOmqTNWx1vNPWHl5eWMGzeO2bNn85Of/KROn1tvvZUpU6Zwzjnn1Ns+IyODrKws5s2bd0z7rqioiCzn5OSwY8eOyMzV9u3b2bRpExkZGSQlJVFVVUV8fDwAu3fvpkOHDgBkZmayevVqxowZU28fR5oJS01NZdu2bZH2goICUlNT6/Rt3749e/fuxe/34/P56vV58803+e1vf8tbb71FQkJCnW3nzZvHxRdfTFxcXJ324uJiVq1axfz58+u0V1RUkJSU1PhJO8FE856wVGBbreWCcFttdwJXmFkBoVmwH0exHhERkWOSnJzMo48+yoMPPljvHq/evXvTt29fXnvttQa3ve22247q24g1+vTpQ15eHgAbN27kwIEDFBYWkp+fT35+PrfcckvkBv1zzjmH559/HgjN3s2bN4/Ro0cDodm4GTNmsGPHDgCqqqp48skngdBMWM1N8LV/Hn30UQAmTZrE3LlzqaysZMuWLWzatImhQ4fWqdPMGD16NC+//DIAc+bMYfLkyQD861//4gc/+AGLFi2iU6dO9Y4xJyeHyy+/vF77yy+/zAUXXEBiYmKd9o0bN9KvX7+jPpfNVaxvzL8ceMY5lwZMBJ4zs3o1mdn3zSzXzHJ37dr1pRcpIiIyaNAgBgwY0OA3E2+77bZGb5TPzMxk8ODBjY47ceLEBi+xnX/++axcuRIIhZWLL764zvpLLrkkUssjjzzCq6++SlZWFsOHD+fSSy+NfKty4sSJTJ8+nW984xuRWkpLS5t0zJmZmUydOpW+ffsyfvx4Zs+ejdfrrVf3rFmzeOihh0hPT6ekpCRy79yMGTM4cOAAl156KVlZWUyaNCkydn5+Ptu2bWtwBnHu3LkNhrMVK1Yc9tLuicacc9EZ2Ows4E7n3Ljw8i0AzrmZtfqsA8Y757aFlzcDw51znzc2bnZ2tsvNzY1KzQCv3vM9+jz7d3be9QNGXXZ91PYjIiKHt2HDBvr06RPrMmJqxIgRLF68mDZt2sS6lJjbuXMn3/zmN/nb3/4W61Ia1dCfWTNb7ZzLbqh/NGfCPgAyzKyXmcUD04BFh/TZCpwbLrIPkAhoqktERAR48MEH2bp1a6zLaBa2bt3Kgw8+GOsyjquo3ZjvnPOb2XRgGaHHTzztnFtnZr8Gcp1zi4CfAf9nZjcQukn/Oy5aU3MiIiInmGHDhsW6hGbjzDPPjHUJx11UnxMWfubXkkPa7qj1eT3w9WjWICIiJy7nXJMetSASa8cyhxTrG/NFREQalJiYSElJyTH94ybyZXLOUVJSUu/bnEcS1ZkwERGRY5WWlkZBQQH6VrycCBITE0lLSzuqbRTCRESkWYqLi6NXr16xLkMkanQ5UkRERCQGFMJEREREYkAhTERERCQGFMJEREREYkAhTERERCQGFMJEREREYkAhTERERCQGFMJEREREYkAhTERERCQGFMJEREREYkAhTERERCQGFMJEREREYkAhTERERCQGFMJEREREYkAhTERERCQGFMJEREREYkAhTERERCQGFMJEREREYkAhTERERCQGFMJEREREYkAhTERERCQGFMJEREREYkAhTERERCQGFMJEREREYkAhTERERCQGFMJEREREYkAhTERERCQGFMJEREREYkAhTERERCQGFMJEREREYkAhTERERCQGFMJEREREYkAhTERERCQGFMJEREREYkAhTERERCQGFMJEREREYkAhTERERCQGohrCzGy8mf3bzPLM7OZG+kw1s/Vmts7M/hzNekRERESaC1+0BjYzLzAbGAsUAB+Y2SLn3PpafTKAW4CvO+f2mFmnaNUjIiIi0pxEcyZsKJDnnNvsnKsC5gKTD+lzDTDbObcHwDn3eRTrEREREWk2ohnCUoFttZYLwm21nQ6cbmbvmdn7ZjY+ivWIiIiINBtRuxx5FPvPAEYBacDbZtbfObe3dicz+z7wfYDu3bt/2TWKiIiIHHfRnAkrBLrVWk4Lt9VWACxyzlU757YAGwmFsjqcc08457Kdc9kdO3aMWsEiIiIiX5ZohrAPgAwz62Vm8cA0YNEhfRYQmgXDzDoQujy5OYo1iYiIiDQLUQthzjk/MB1YBmwA5jnn1pnZr81sUrjbMqDEzNYDK4AZzrmSaNUkIiIi0lxE9Z4w59wSYMkhbXfU+uyAG8M/IiIiIl8ZTQphZvZ14E6gR3gbI5ShTo1eaSIiIiInr6bOhD0F3ACsBgLRK0dERETkq6GpIWyfc+6vUa1ERERE5CukqSFshZndD7wKVNY0Ouc+jEpVIiIiIie5poawYeFfs2u1OWDM8S1HRERE5KuhSSHMOTc62oWIiIiIfJU06TlhZtbazB4ys9zwz4Nm1jraxYmIiIicrJr6sNangf3A1PBPKfCnaBUlIiIicrJr6j1hpznnLqm1fJeZrYlGQSIiIiJfBU2dCTtoZiNqFsIPbz0YnZJERERETn5NnQn7ITAnfB+YAbuB70SrKBEREZGTXVO/HbkGGGhmKeHl0qhWJSIiInKSO2wIM7MrnHPPm9mNh7QD4Jx7KIq1iYiIiJy0jjQT1iL8a6toFyIiIiLyVXLYEOac+2P417u+nHJEREREvhqa+rDW+8wsxczizOxvZrbLzK6IdnEiIiIiJ6umPqLivPDN+BcA+UA6MCNaRYmIiIic7JoawmouW54PvOSc2xelekRERES+Epr6nLDFZvYpoQe0/tDMOgIV0StLRERE5OTWpJkw59zNwNeAbOdcNVAGTI5mYSIiIiInsyM9J2yMc265mf1PrbbaXV6NVmEiIiIiJ7MjXY48B1gOXNjAOodCmIiIiMgxOdJzwn4V/vW7X045IiIiIl8NTX1O2D1m1qbWclszuzt6ZYmIiIic3Jr6iIoJzrm9NQvOuT3AxOiUJCIiInLya2oI85pZQs2CmSUBCYfpLyIiIiKH0dTnhL0A/M3M/hRe/i4wJzoliYiIiJz8mhTCnHOzzOwj4Bvhpt8455ZFrywRERGRk1tTZ8IANgB+59ybZpZsZq2cc/ujVZiIiIjIyayp3468BngZ+GO4KRVYEK2iRERERE52Tb0x/zrg60ApgHNuE9ApWkWJiIiInOyaGsIqnXNVNQtm5iP0xHwREREROQZNDWFvmdmtQJKZjQVeAl6LXlkiIiIiJ7emhrCbgF3AJ8APgCXA7dEqSkRERORkd8RvR5qZF1jnnOsN/F/0SxIRERE5+R1xJsw5FwD+bWbdv4R6RERERL4SmvqcsLbAOjNbBZTVNDrnJkWlKhEREZGTXFND2C+jWoWIiIjIV8xhQ5iZJQLXAumEbsp/yjnn/zIKExERETmZHemesDlANqEANgF4MOoViYiIiHwFHOlyZF/nXH8AM3sKWBX9kkREREROfkeaCauu+XAslyHNbLyZ/dvM8szs5sP0u8TMnJllH+0+RERERE5ER5oJG2hmpeHPRuiJ+aXhz845l9LYhuHni80GxgIFwAdmtsg5t/6Qfq2AnwL/PMZjEBERETnhHHYmzDnndc6lhH9aOed8tT43GsDChgJ5zrnN4fdOzgUmN9DvN8AsoOKYjkBERETkBNTU1xYdi1RgW63lgnBbhJkNBro55/4SxTpEREREmp1ohrDDMjMP8BDwsyb0/b6Z5ZpZ7q5du6JfnIiIiEiURTOEFQLdai2nhdtqtAL6ASvNLB8YDixq6OZ859wTzrls51x2x44do1iyiIiIyJcjmiHsAyDDzHqZWTwwDVhUs9I5t88518E519M51xN4H5jknMuNYk0iIiIizULUQlj4kRbTgWXABmCec26dmf3azPTOSREREflKa+q7I4+Jc24JsOSQtjsa6TsqmrWIiIiINCcxuzFfRERE5KtMIUxEREQkBhTCRERERGJAIUxEREQkBhTCRERERGJAIUxEREQkBhTCRERERGJAIUxEREQkBhTCRERERGJAIUxEREQkBhTCRERERGJAIUxEREQkBhTCRERERGJAIUxEREQkBhTCRERERGJAIUxEREQkBhTCRERERGJAIUxEREQkBhTCRERERGJAIUxEREQkBhTCRERERGJAIUxEREQkBhTCRERERGJAIUxEREQkBhTCRERERGJAIUxEREQkBhTCRERERGJAIUxEREQkBhTCRERERGJAIUxEREQkBhTCRERERGJAIUxEREQkBhTCRERERGJAIUxEREQkBhTCRERERGJAIUxEREQkBhTCRERERGJAIUxEREQkBhTCRERERGJAIUxEREQkBqIawsxsvJn928zyzOzmBtbfaGbrzexjM/ubmfWIZj0iIiIizUXUQpiZeYHZwASgL3C5mfU9pNu/gGzn3ADgZeC+aNUjIiIi0pxEcyZsKJDnnNvsnKsC5gKTa3dwzq1wzpWHF98H0qJYj4iIiEizEc0Qlgpsq7VcEG5rzPeAv0axHhEREZFmwxfrAgDM7AogGzinkfXfB74P0L179y+xMhEREZHoiOZMWCHQrdZyWritDjP7BnAbMMk5V9nQQM65J5xz2c657I4dO0alWBEREZEvUzRD2AdAhpn1MrN4YBqwqHYHMxsE/JFQAPs8irWIiIiINCtRC2HOOT8wHVgGbADmOefWmdmvzWxSuNv9QEvgJTNbY2aLGhlORERE5KQS1XvCnHNLgCWHtN1R6/M3orl/ERERkeZKT8wXERERiQGFMBEREZEYUAgTERERiQGFMBEREZEYUAgTERERiQGFMBEREZEYUAgTERERiQGFMBEREZEYUAgTERERiQGFMBEREZEYUAgTERERiQGFMBEREZEYUAgTERERiQGFMBEREZEYUAgTERERiQGFMBEREZEYUAgTERERiQGFMBEREZEYUAgTERERiQGFMBEREZEYUAgTERERiQGFMBEREZEYUAgTERERiQGFMBEREZEYUAgTERERiQGFMBEREZEYUAgTERERiQGFsMN48pMn+WfRP2NdhoiIiJyEfLEuoLlaVfoRz32YC8BNZ97EtN7T8Hl0ukREROT40ExYI97e80Hk86wPZjHouUFU+CtiWJGIiIicTBTCjsKZL5zJ71b/LtZliIiIyElAIewoPb32aWb/6wmcc7EuRURERE5gusnpMAIVXSnf8pPIcnKPx/Amb+Xxj3/P4/9YQeX+00jovBCArqV3M+uicxiY1hozwzmHmQGwt7wKw9iwo5Stu8tJa5NE/7TWtEqMo8of5ECln+R4Lwk+D2ZGRXUlOw4UU3rQkZQYpDp4kG17SmmV7KdH6+50SuoEBsVlezmlRXs8HmVpERGRE41C2CECBCOfy7f8iC0zJ1K49yC3zl/L2xt/RMvT78S8FdBiLQkt1kb6bk+5nf+3vGn78Jediq/FZgBcMA7zVB/XYzhavsAp+L07I8uGB5wXZ9WY89HKm0pp8DN8Fo/hoXVcZ/ZWb8fvqvCaD/DgNS+pSb3ZV/055f69VATL6JiYhtc87K/ewylJqWwv30qPlqezs3wHWIAzO32dz8t30Tk5FVw828s3c0piTyqCe9lWtoUWvhakpXTBzNE+sRPb9udzausMSspL2FO1h+4te7G/aj/x3gQCHKRLiy7sqdxNSnwbPt29gdNap5Pg81J4YDs9UrpjZsR5vGzc+28y2mRQcnAPB/3lnN4ug51lO2mb2JbdFbtpm9CWLaVb8JiHjDYZ7K7Yjde8OBzJvmQKDhTQrVU3Nu7ZSEabDNbvXk/fdn3ZU7mH1gmtAdhVvoueKT1pEd8C5xxx3jiSfcmUVZeR4E3AOUfxwWI6JXeiIlBBvDcef9CPxzx0TOpIdbAa5xw+jw+PeTjoP0iSLwmPeTAMv/OT6E0kSBCveSM/fufHa17Kq8sjY5oZhnGg+gBxnji85qUqWEWSL4mgC9IqvhX+oD/0Z8HjozJQSZIvCa95iffG4zUvAF6PF5/58Hl8kf/BcM4RcAGCLohhxHnjiPPEfbl/gEVETlB2ol1Wy87Odrm5uVEb/w+/+R9GvbCBGVd5Wbt7Fvn3nh9Z55yjojrImHnncSDwOQBtEtqyt3LPF95vsLoN4AhWdiIhsYw4bwALplBWnkjbNvsoK8nijC5xrCufj9fflTj/qcS3+IzS4GcApASGUHLAT1zL/xB0hqtqiyepAE8wieqDXcDjx5ecX2+/VtmDgLcE8xzEPAGcP5mgvxXexJ0Eq9vgAsl44krAAgSr2hOs6oAnYRfehM8J+pPx+MpD9Vd2wDkv5vHjiS8h6G+Bx1cWOm/hoBn0t8TjO/CFz5WcOIw4HKEwb2YEqSaOFgSowDkvcZ54qtwB4l0nquxzEl1XnKcUv6si2U5hv9tGh7gMiqs30cnXn0q3h2RPB/ZX7wdvKW19PdlZ9Slpif3ZfnATqckZVPnB63GU+PPomtiHiuoAewP/ITUpk+0HN5KW3BucsbVsA71aZQJQXL6Xtkkt8JJAld9ISUwAjM/2b+S0lH6UBUrw+xPo3LIdgaCjzL+HtokdMIwdB0roltIFA0orDxLnNVLiU6gIlOMCLWiVFArO5ZVeEuKraJfQiaALsq/cT5eU1lQEy0iytiTGG/GeeCoDQZJ8HiqC5SR5W1BWGaB1chw4H60TWpDoSyDOE08w6CHB58NjhpmHg1VBWifF4/MYXo9FZtY9Bh4zav5L75zD67Hw74+Bgc9jeMyI84a2rQnZIvLFmdlq51x2Q+s0E3aIgPvvTNh9UwbUWWdmJMV7+ccVfzvm8asD1cR5v8hMwZ1fYNujFww6qgKhS6bllQGqg0Eqq4OYQcmBKqoCAdq3SCDoHDtLK3HO0Soxjkp/gK27y0nwedlRWkGvDsnsK6+mV8eWFO45SCAYJDneR1FpGad2SMQfDLDnYDUEffRo34rqgOPjgr3sLC2jXUsvu8uqMIPenduwdfd+OrZKoHViAkEXYPfBMrbu2Y/f7yO1TQvat0iiOlDF5/urSE6AeJ9RUe1nX3mQeK+PDTtKGdytAzsP7Cfe58fniaP4gJ/OKS3Yua+Csqpq9h2spEUieMxHUryHTZ/vY2Bqez7bfYCubZLxOC+lldXs3h9gYPc2rM7fw6DuKXy0fTdDurehuKyUSneAzw/sJr1jB3aWltEiLpmqYDUHqxzx8ZXs2uPB78op2R+gU5sqWsYlUOY/iMdga3EV3TpWUFqajNcCOH8V7Vsks31PGV3aOpI9Rn7pPnq0TqHKs4NdJS3p2M6oqvDRKimJ/cESUmjDQd8+qva35JRkL3tcKUnBtsR5/fxndwXd27ZkW+lOOqe05/OqHbSNS6WkehctrRX7qstJ9iVSHb+D6gNtiUsMEgj4KAvsJyUxnqoqw5NwEF+gFdV+8McVE1fVngrvNloEOlLpLcAOdiLRW43fD/FmBOO34zvYAR8Ov6cKb9CLP2k7SRVpVCQUE1/ZDQ/VVCX4iT/YgerE7bSvPIWqhDx6VrXD79tIsj8Js1IsYR/xla2pSlhF16okKuL+ToegB3/gAF6M6rgy2lcnU0ku/rgy2lUnsC/uTVo440CgCL+3mhRn7N2+A2cBzBOgNBCHAc6ClOFwniBxDrYWfhz5+7Av/KsBO2r9PdkNWBP/f7ahfo1Gnib2Pdy+vcHQLG7QE8SChteFZjYDHj/xgThc6Kjxe6tJrE7A4aHaW0VidSIA5fFltKtMooo4At4qEvxJmAWockYiSVRZJYlV8QTwYXEB4gNJePyVBD1J+DyhAFoVCNLa66eKJMCD319NUnwcQXwEAg6PGQleCFocldUBkuN9eJzhDwZJwg++ZIIO4qoDxCUnEKyuxJuQjM/jxV9eTVLLBDAvYFT7Ha08HjwBP8GEOKqcjxYJ8Xgwqqv8tGiZgMeXQKUzWjkHycng8XGgKkj7lkl4vR7AiEuIx+eNJ+AgPi4Bb2ISHl88ZuCJTyYhMQlPfBzeuHiS4xNCM8W+ODw+H16fN3Ts4QBcE4TNUMCVOhTCDvFJ5RbODX+emt3t8J39VbBrA7z9ACS3g6wroNUpsHcbbHod2qdDXBL0nQyBatjxCXEHdsLez8A8cHAPVJVB1QHwV0LPs6F4I7TqDKf0g2A1dOwNrbpC0A+4UP9AFZQVh7bZvz3UVn0QfAmwdyu4IJT8BzpkQDAA+4vA44U9+bB/J7TpBr5EKFoDFaXQpjvs2wZdB8PB3RDXAso+B38lnviWJFaWktiiA2z/V+i4W54CB3aGxkhqGxrflwSBSnCOyL8cbXuFxjUvtE4L1RishpRUqNgL+wqhaxZ8ejB0zlK6QFI7KP43dMlibMVeiEuG0kBozOKNUD4YyksgpWvoPKakhdaZJ7SvQHf4z2ehfWz7J3Q4PbR++78gdUioretgyP0UElqF6t/1KbQ7DbYfDJ3PhNahcwgQnxw61lZd4NMSiG8BpT4o2wXJ7XFlJbC/La66EoqrcN4WsK0U522Bq66AoMNtbomrOADOcBYHgWpc0EKnKvyr+8zAhXfrwm154c9BcC68vvY2ziAYanPBmvX/Hee/6w9pc4dsU2vb0PiNj1Pz21tnm3DN//Vp+NctDfylKWqgrSba1J4lrflcGP718/CvZbX67G5gu6pan/fX+lxe67O/kc+BBmo7GRx6XIc75trnqfa5rH2Oa/8efHU5QmcvQN0zVSNI3W++HYwP/TVJroR9LR1BM+KroDoOKhIgaJBUARWJRlVcaLlFmaOkreHFCALlcY5k5w39F9YTmsX0mKPcjHjz8v/bu/tYyeq7juPv73mYmft8YZ+67bJAgAZ56rYlFGhsa6kKNWHTSCmkpK2gRFMsVYPR1KipMaE1amz6oFgQS6K10tisFUtMoaVWtoVKwS6NsIFdWPqwT3fv871zHr7+8Tv37t27d9nB7Mxc7nxeyeT85nd+55zvnd+dme/5zW/ODMw446MJeQRlHKZfGCVHYhiIasRmFGVJGdeI4xizkjJKIEqIpzLKkb7wXhHFxFYwaw1q9RRLEiyKiQzmqVNv1IniBOIwEjsfNejv6yNJU5rUGOjvJ6r1kdTq5EmD4cEh0lofUb2PvoEhavV+0rQGQD2NSKKI2IwogjSOSONocbR2rdPHkcvc/LELuePLJXfcHPPV3z0654sihyfug69+tG3H7pikAcuveRalIUFK+yGbgf51IdlpjIZkbt9jsPHCkACWeUjuTj8bXvpvOP2s0P7IC+FxwkMSOfVTiNOQiBZZ2LY2GJLOA8/gjWGYm8anx/DCKEugNMpkGK+fjnuM5yU+dRjPMtzq+Pwc3rcRn5sOiY7HeOH41BhlOozPTuLxID47jRNDMohPHcT7N8D0Ybx+Gl7k0Gzi6SDMjuO14ZBo5M0q90rwPMcthryoEo0SyjK0K4qQRJUlS6YQdlc45cbiCIsjiKKqHIN5WMYJFhsWGcQx5hlW64M4Ci/xcYxFEVbMQN9w2D6bhv4RzPPQLkkxcvAc6xsK+yrmsXp/WJ/NQKOqnzmIDW2CKArryiZOjJcxUV8D4jhs6w71wTDEM3MQBjdiURy2sygk92k/1YBNuJ+EURrmx6E2vNjWLAr/Z3EtPCYYNCehPnz0sZqfCAk4hJOBxsjRdflsOKGw6u2zzMObkkXh4NlsOLZ7iGf2CPSNVseKwnJpfBhWzIcTkaQW4jHCftKBUF6IM58Lf+eCIgvrknRxX5R5iGXx/WnJG5Ud/wUdO9GXdiyqMuzy6HZm4ZhxejSmhXULj+nCCVYxX+2jekwWtvEqw29OQi08xo7j+SwljkcJJSVMHaTsGw3rcEoPbTyuhbKFj03LuTGoj1KYQ1lSzjfxKMMBj1PKLMNLp2QOjxsUOKWXcGSaIp2laIxSFmWY6zvbpMgmKOrDFEVBNDWPF3NkjTpZ4UTmlKVT5k3mvEmNGnhJfHiMYrAfd6N2aJqZgYQkNqIipzadc6TmNIhxdxqzE2TeR5EC7qRTc9SaEZNDEfNesPFIyVx/QeYJlpUMZJAloW087yQORQyRQ/+UM9mARnUO2sTpK0PZytDGqofbHDZMwEz1LxaVEFe3aBW/xecR5HFIOMsIiigsl5bdFupssV3kUMtgqt/ADI+gdNiyv+Sl16YUsVHEBnFEERmWxngSU0YRUS0klSQJliRsvvwXedt1t7X179THka9ASsJx5zXf/CQ8/KcrbzCwEbZcCnv+WpEpRAAACuZJREFUE15zCUy8FO4PbgovXFM/hQPPwLlXwZ5vhZGtzdtCAjOwAab2hxf5g8/CwPowKjPxIzjnKvjxk7B/19EXuoENgMPZb4eJfWGkycuQ5KT9YfuhzeENIKmHF0YMjxshkWg2wy3Ljl0uKZfH1Gd4Vi23vj/UF008r9rsbeLZNvz5ZmiXbQ3bL+5j45JyHm7ZgWrfp+FZBt4ANq/wwM4uKTeqG0ANmAPiarmgDtbEakNYrYalg1iaYmkK8XqsmWCN14YnXiPGhuNQjmMsjojiUF6sSxaSljgkEHFyTN3J1pPER+ui0H4xEUrikPAcU15hfZIsSZ6i4+sW9xmf+M1WRHpaUeRk+TxlnlWv3fN4kVHm85RZk6II6/LmPGXRJM/nyZtN8rkZ8jwnb85R5BnF/Ax5UVIWGWWeh/3NT5OVhpc5RVZQ5M1wAlsUeJ7hRU6R5VhZUmbh5N2KksZEk7k+o8gKEkKCbaXjRUlUOlY6FI65L973wokdogJGJkomB4y+rAznBQ6bD8D4IIyMNUMCWlS3EpIckgKSFU6an5yZhjYnYS9HSdgy53A68BN+aWoK9v4X/N01xzb41Ydgy5tP6TG9KCg3XU4xOUVZ30Y5OUlxaIoyu4IiuYhycopyapJi3yTl1DTlI09Vic2LKyZT5ULiVNWT5ycP4pWokpsoTauEZ9myKkfDfcvq08VyVKtBmhItaR+WC+2WLo/fx3Hb1WohGdF8CxGRRXGcEMcJ1LsdSRe5k+cZs7OzzM9MMjczwez0GM2Zad626STTjtpMSdgyjTKMgm3Ki2MTsNufhNPOOuF2nufkhw6T799PcfgQxfg4xfgExeTE0SRqcopyaiokWUuWPjNzwv0uSlPiwUGioSGi/v7FhCTq74N0OCQl6fGJ0HHlNA2JzJK66ERtVyqnqUZdRETk1cOMJK0xlNYYGh45efsOamsSZmZXA39F+Ozo8+5+57L1deALwJuBQ8D73H1PO2M6mcb8QSA6OqFycBP5zTspDh0mf/a7ZC+8QHPvXpp79tLcu5dibAzPMoqJCcKkpuNZo0E0NEg8EJKoeGiQZNOmUDc4RDQ4GMpDQ0SDYX00FOoXEi+r1zXKIyIisoa0LQkzsxj4DPDzwD7gMTPb4e5PL2l2CzDm7uea2Q3AJ4D3tSumVkympwNHyObO5Cf5xUw/+jzNz7/12EZpSm3LFmpnnknfGy7B0pR4dJRk40aSDRtI1q0jGh4hHhkmHgpzlERERESWaudI2GXAbnd/DsDMvghsB5YmYds5euGr+4FPm5l5F7+y6dVo1kUPHOJI+m0GrryS0euvJ3nNJuLRUWpbt5Ju3hy+XSEiIiLy/9TOTOJ1wItL7u8D3nKiNu6em9k4sA44uLSRmd0K3AqwdevWdsULwKYN5wPfZfJnt/GmOz9Nsm5dW48nIiIivelVMZzj7ncBd0G4Tlg7j3XTHffy3DU7ueziK9p5GBEREelx7fya20vA0u9+buHo5a+Pa2NmCTBCmKDfNWbGOUrAREREpM3amYQ9BpxnZmebWQ24AdixrM0O4INV+TrgoW7OBxMRERHplLZ9HFnN8boNeJBwiYp73H2XmX0ceNzddwB3A/eZ2W7CD8Hd0K54RERERFaTts4Jc/cHgAeW1f3hkvIc8N52xiAiIiKyGunS5yIiIiJdoCRMREREpAuUhImIiIh0gZIwERERkS5QEiYiIiLSBUrCRERERLpASZiIiIhIF9ir7QL1ZnYA2Nvmw6xn2Y+Iy6qgfll91Cerk/pl9VGfrE6d6Jcz3X3DSitedUlYJ5jZ4+5+abfjkGOpX1Yf9cnqpH5ZfdQnq1O3+0UfR4qIiIh0gZIwERERkS5QErayu7odgKxI/bL6qE9WJ/XL6qM+WZ262i+aEyYiIiLSBRoJExEREemCnk7CzOxqM/tfM9ttZr+3wvq6mf1Ttf47ZnZW56PsPS30y2+b2dNm9pSZfd3MzuxGnL3kZH2ypN0vm5mbmb4F1mat9ImZXV89V3aZ2T90OsZe1MLr11Yze9jMnqhew97djTh7iZndY2b7zewHJ1hvZvapqs+eMrM3dSq2nk3CzCwGPgNcA1wA3GhmFyxrdgsw5u7nAn8JfKKzUfaeFvvlCeBSd78EuB/4ZGej7C0t9glmNgTcDnynsxH2nlb6xMzOA34feKu7Xwh8tOOB9pgWnyt/AHzJ3d8I3AB8trNR9qR7gatfZv01wHnV7Vbgcx2ICejhJAy4DNjt7s+5exP4IrB9WZvtwN9X5fuBq8zMOhhjLzppv7j7w+4+U93dCWzpcIy9ppXnCsCfEE5U5joZXI9qpU9+DfiMu48BuPv+DsfYi1rpFweGq/II8KMOxteT3P0R4PDLNNkOfMGDncComW3uRGy9nIS9Dnhxyf19Vd2Kbdw9B8aBdR2Jrne10i9L3QL8e1sjkpP2STV8f4a7/1snA+thrTxPXg+83sy+bWY7zezlRgLk1GilX/4YuMnM9gEPAL/ZmdDkZbzS951TJunEQUTawcxuAi4F3t7tWHqZmUXAXwAf6nIocqyE8PHKOwijxY+Y2cXufqSrUcmNwL3u/udmdgVwn5ld5O5ltwOTzuvlkbCXgDOW3N9S1a3YxswSwtDxoY5E17ta6RfM7F3Ax4Br3X2+Q7H1qpP1yRBwEfANM9sDXA7s0OT8tmrlebIP2OHumbs/DzxDSMqkfVrpl1uALwG4+6NAg/D7hdI9Lb3vtEMvJ2GPAeeZ2dlmViNMkNyxrM0O4INV+TrgIdeF1drtpP1iZm8E/oaQgGmeS/u9bJ+4+7i7r3f3s9z9LMI8vWvd/fHuhNsTWnn9+gphFAwzW0/4ePK5TgbZg1rplxeAqwDM7GcISdiBjkYpy+0APlB9S/JyYNzdf9yJA/fsx5HunpvZbcCDQAzc4+67zOzjwOPuvgO4mzBUvJswqe+G7kXcG1rslz8DBoF/rr4n8YK7X9u1oNe4FvtEOqjFPnkQ+AUzexoogDvcXSP5bdRiv/wO8Ldm9luESfof0sl9e5nZPxJOSNZXc/H+CEgB3P2vCXPz3g3sBmaAX+lYbOp7ERERkc7r5Y8jRURERLpGSZiIiIhIFygJExEREekCJWEiIiIiXaAkTERERKQLlISJyJpiZoWZfd/MfmBm/2pmo6d4/3uq625hZlOnct8i0luUhInIWjPr7tvc/SLC9f0+3O2ARERWoiRMRNayR6l+iNfMzjGzr5nZ98zsW2Z2flW/ycz+xcyerG5XVvVfqdruMrNbu/g3iMga1bNXzBeRtc3MYsLPw9xdVd0F/Lq7P2tmbwE+C7wT+BTwTXd/T7XNYNX+Znc/bGZ9wGNm9mVdcV5ETiUlYSKy1vSZ2fcJI2A/BP7DzAaBKzn6U1cA9Wr5TuADAO5eAONV/UfM7D1V+QzCj18rCRORU0ZJmIisNbPuvs3M+gm/4fdh4F7giLtva2UHZvYO4F3AFe4+Y2bfIPzQsojIKaM5YSKyJrn7DPARwg8mzwDPm9l7ASx4Q9X068BvVPWxmY0AI8BYlYCdD1ze8T9ARNY8JWEisma5+xPAU8CNwPuBW8zsSWAXsL1qdjvwc2b2P8D3gAuArwGJmf0QuBPY2enYRWTtM3fvdgwiIiIiPUcjYSIiIiJdoCRMREREpAuUhImIiIh0gZIwERERkS5QEiYiIiLSBUrCRERERLpASZiIiIhIFygJExEREemC/wNRK7yvei+TmwAAAABJRU5ErkJggg==\n"
     },
     "metadata": {
      "needs_background": "light"
     }
    },
    {
     "output_type": "display_data",
     "data": {
      "text/plain": [
       "<Figure size 432x360 with 1 Axes>"
      ],
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAFNCAYAAAANRGjoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd3wUdf7H8dcnDQJp9N5EEFFpIkVFEAFRsKIIZz+Fn/08T+/senbl9CyHJ3h6KiggoIgCShHxBOld6b33hFQ22f38/pgBN4VkBXY3yX6ej0ceyc58d+Y9s5v97NSvqCrGGGMiV1S4AxhjjAkvKwTGGBPhrBAYY0yEs0JgjDERzgqBMcZEOCsExhgT4awQlHMi8ouIdAt3jtJCRB4Xkf+Ead4ficgL4Zj3qSYiN4rI1BN8rr0nSxkrBCEkIptFJFtEMkRkt/vBkBDMearqWar6QzDncZSIVBCRl0Vkq7uc60TkERGRUMy/iDzdRGS7/zBVfUlV7wzS/EREHhCRlSKSKSLbRWSsiJwTjPmdKBF5VkRGnsw0VPVTVe0VwLwKFb8TfU+KSJybfZ27fjeLyIci0vj3TsvkZ4Ug9K5Q1QSgDdAWeCzMeX43EYk5zqixwCXA5UAicDMwGHgrCBlERErb+/ct4E/AA0BVoDkwAehzqmdUzGsQdGGc9zjgSuAPQDLQGliE8577XcK5/kolVbWfEP0Am4Eefo9fAyb5Pe4EzAFSgWVAN79xVYH/AjuBQ8AEv3F9gaXu8+YArQrOE6gLZANV/ca1BfYDse7jPwKr3Ol/BzTya6vAvcA6YFMRy3YJkAM0KDC8I+AFTncf/wC8DMwHDgNfFchU3Dr4AXgRmO0uy+nA7W7mdGAj8H9u28puGx+Q4f7UBZ4FRrptGrvLdSuw1V0XT/jNLx742F0fq4C/AtuP89o2c5ezQzGv/0fAUGCSm3ce0NRv/FvANne9LAK6+I17FueDcKQ7/k6gA/Czu652Af8C4vyecxYwDTgI7AEeB3oDHiDXXSfL3LbJwAfudHYALwDR7rjb3HX+T+CAO+424Cd3vLjj9rrZVgBn43wJyHXnlwF8XfD/AIh2c21w18kiCryH3HY93Nez0Lhi/r+Keq3vcF/rH4EpwH0FprEMuNb9u4Xf+lsD9A/3Z0jQPpvCHSCSfgr8A9R3/2Hech/Xc//JLsfZUuvpPq7hjp8EjAGqALFAV3d4W/cfsKP7T3WrO58KRczze2CQX54hwHvu31cB64EzgRjgSWCOX1t1/ymqAvFFLNsrwKzjLPcWfvuA/sH9oDkb58N6vN8/a0nr4Af3n/gsN2MszrftpjgfRl2BLKCd274bBT64j/Ph8D7Oh35r4Ahwpv8yueu8PrC84PT8pnsXsKWE1/8jd3k6uPk/BUb7jb8JqOaO+wuwG6jolzsXuNpdN/HAuTiFM8ZdllXAg277RJwP9b8AFd3HHQuuA795fwkMc1+TmjiF+uhrdhuQB9zvziue/IXgUpwP8BT3dTgTqOO3zC8U83/wCM7/wRnuc1sD1X7P+6uo6RbzWn/iLmM8cAsw2699S5yiWsFtsw3ni0YMv31pahnuz5Fg/JS2TetIMEFE0nHeZHuBZ9zhNwGTVXWyqvpUdRqwELhcROoAlwF3qeohVc1V1Vnu8wYDw1R1nqp6VfVjnA+zTkXM+zNgIDi7VoAB7jBwPsheVtVVqpoHvAS0EZFGfs9/WVUPqmp2EdOujvPBU5Rd7vijRqjqSlXNBJ4C+otIdHHrwO+5H6nqL6qa566HSaq6QR2zgKlAl+PkOJ6/q2q2qi7D+UbY2h3eH3jJXefbgbeLmUa1Ypbf35eqOt9dx5/i7CIEQFVHquoBd9lex/lAOsPvuT+r6gR33WSr6iJVneu234zzQd7VbdsX2K2qr6tqjqqmq+q8ogKJSC2cdfygqmaq6l6cb/gD/JrtVNV33HkVfP1zcQpNC0Dc91Ag6wKcLZsnVXWN+xouU9UDRbQLdP2W5Fl3GbNxip//e/xG4AtVPYKz/jar6n/dZV6C86Xl+lOQodSxQhB6V6tqIs631Rb89gHZCLheRFKP/gAXAnWABsBBVT1UxPQaAX8p8LwGOLtBChoPdHYLy0U4u03+5zedt/ymcRDnG1o9v+dvK2a59rtZi1LHHV/UdLbgfLOvTvHroMgMInKZiMwVkYNu+8vJX3QCsdvv7yzg6AH8ugXmV9zyH+D4yx/IvBCRh0VklYikucuSTP5lKbjszUXkG/fEg8M4xfto+wY4u1sC0QjnNdjlt96H4WwZFDlvf6r6Pc5uqaHAXhEZLiJJAc470JyBrt+SHFsOVU3H2dI+WvAG4hRncNZJxwLvxRuB2qcgQ6ljhSBM3G+vHwH/cAdtw/mmnOL3U1lVX3HHVRWRlCImtQ14scDzKqnqqCLmeQjnG/MNOAfcRquq+k3n/wpMJ15V5/hPophFmo7zj9PAf6CIdMT5Z//eb7B/m4Y43yj3l7AOCmUQkQo4xe0fQC1VTQEm4xSwkvIGYhfOLqGichc0A6gvIu1PZEYi0gXnGER/oIq7LGn8tixQeHn+DawGmqlqEs6+9qPttwGnHWd2BaezDWcrsrrfek9S1bOKeU7+Caq+rarn4uxeaY6zy6fE57nzblpCG3DeXx1EpH4xbTKBSn6Pi/rQLphnFDBQRDrj7EKb6ZdrVoH3YoKq3h1A1jLHCkF4vQn0FJHWOAcBrxCRS0UkWkQquqc/1nc3s6cA74pIFRGJFZGL3Gm8D9wlIh3dM2kqi0gfEUk8zjw/w9k3eh2/7RYCeA94TETOAhCRZBEJeDNYVafjfBiOF5Gz3GXo5C7Xv1V1nV/zm0SkpYhUAp4Dxqmqt7h1cJzZxuHsPtkH5InIZYD/KY17gGoikhzochTwOc46qSIi9YD7jtfQXb53gVFu5jg3/wAReTSAeSXi7IffB8SIyNNASd+qE3EOzmaISAvA/0PqG6COiDwozmm9iW5RBme9ND561pX7/poKvC4iSSISJSJNRaQrARCR89z3XyzOh3EOztbm0XkdryAB/Ad4XkSaue/fViJSrWAj9/01DfhSRM4VkRh3me4SkT+6zZYCA9z/j/Y47/GSTMb59v8cMEZVj+b+BmguIje704t1l/PMAKZZ5lghCCNV3Ydz8OppVd2Gc8D2cZwPg20436qOvkY343xzXo1zbOFBdxoLgUE4m+aHcA743lbMbCfinOGy290nfjTLl8CrwGh3N8NKnOMSv0c/nG9U3+KcJTIS50yU+wu0G4GzNbQb51vYA26GktZBPu6m/QM4H9iHcLZyJvqNX43zjW+ju3lf1O6y4jwHbAc24XwjHYfzzfl4HuC3XSSpOLs8rgG+DmBe3+Gst7U4u8tyKH5XFMDDOMucjvOFYMzREe666QlcgbOe1wEXu6PHur8PiMhi9+9bcArrrzjrchyB74pJcud/yM1+AOdEBHBe/5bu+p9QxHPfwHn9puIUtQ9wDuQW5TqcD+4xOFtLK4H2OK8NOMebmro5/k7+LzpFco8HfIFzVtJnfsPTcb5UDMA5U283zv9HhZKmWRbJb3sGjAk+EfkB50yOsFzdezJE5G5ggKoG9E3ZmLLCtgiMOQ4RqSMiF7i7Ss7AORXzy3DnMuZUs6vrjDm+OJyzZ5rg7OoZjXMcwJhyxXYNGWNMhLNdQ8YYE+GsEBhjTIQrc8cIqlevro0bNw53DGOMKVMWLVq0X1VrFDWuzBWCxo0bs3DhwnDHMMaYMkVEthxvnO0aMsaYCGeFwBhjIpwVAmOMiXBWCIwxJsJZITDGmAhnhcAYYyKcFQJjjIlwQSsEIvKhiOwVkZXHGS8i8raIrBeR5SLSLlhZjDHGHF8wtwg+AnoXM/4ynA5SmuF0wP7vIGYxxhhzHEG7slhVfxSRxsU0uQr4xO0zd66IpIhIHbfbPGNKncxd25j677/iycoMdxQToRpcdAXnXznolE83nLeYqEf+rvi2u8MKFQIRGYyz1UDDhg1DEs6Yo3LWrWX2G49QdfZaWnjCncZEsmW+CeWuEARMVYcDwwHat29vHSiYkDiybStLX76PyjPXUVPg19MgtXsXLrygHyDhjmciUI+GzYIy3XAWgh1AA7/H9d1hxoRVxs7VzPrHA9SZsY3KHvi1iTC7Sweeve9dqiRWCnc8Y065cBaCicB9IjIa6Aik2fEBE1bZqcwY+kcSP13FadmwtpEw64Lz6HLFc7zeqgEx0Xa2tSmfglYIRGQU0A2oLiLbgWeAWABVfQ+YDFwOrAeygNuDlcWYkmRl7efrRy+hxfceDiYLn/a8lEtufITXzq5tBcCUe8E8a2hgCeMVuDdY8zcmUF9+PISs8R/Sbi2srV+R7Mf/zZsXd0TEjgOYyFAmDhYbEwwHMo7wzrtP0X381zTKhl87nk7fdz6jQlJiuKMZE1JWCEzEUVWmrdjK7G9upcPUPVRJh0pvPkW/Xn8IdzRjwsIKgYko6/dmcOe/J9Oqwgv0nZxLnUOQ+NozNO41INzRjAkbKwQmYszbeIAHPprOtXnPceHXPpKPRNH4gw9I6NQp3NGMCSsrBKbcS8vK5ZVvVzNn2WSuSP2EK6f6yK0YTeP3hlkRMAYrBKac+3HtPu75dDHNY76mQ+r3XDNNOdKwKmd9NJbYunXDHc+YUsEKgSmX8rw+3v5+PW/PWMflMo6+P86lyS6IOaspZ474nKhKdoWwMUdZITDlzsFMD/ePWsyc9Xt4jCG0n3yA2DyIvecmmt7/uF0fYEwBVghMuZKa5aHbkJmk56bRO+4lOn5xhKwEocErQ6jZrU+44xlTKlkhMOXG7PX7eXjsMtJzU7ne+zo3TjhCTo042o+ZSlytWuGOZ0ypZYXAlAvLt6dyy4fzia64mj/u+ZDrfvLirRrNueNnEFOterjjGVOqWSEwZd6iLQe5/b8LaB03k6sWTaLDGkVaVaXlu2OIsiJgTImsEJgyK8/r48kJK/l84Tauy/iAq35eRVIWVOh3Pk2eG45ER4c7ojFlghUCUyapKn8avZRJK3bxF98Iuk9bxZ5qUOedd6jeuUe44xlTplghMGVOtsfLX8cvZ9KKXfwpdiwXfbmMnXWEjp+OJ6numeGOZ0yZY4XAlCmb92fy58+XsmRrKo/VmkG7j+dxMAHOGv6hFQFjTpAVAlNmLN+eSv9hP5Pj8fKv2PFUHzEXUaj3+kvUbWb3DDLmRFkhMGWCJ8/Hg2OWkpPrY0r8tzB6LltrC62GvE2d8+yYgDEnwzpjNaWez6c8On45G/dl8s9zfiDjq+msbCgkfvimFQFjTgHbIjClWp7Xx52fLOSHNXu4vtFnZH+6jIRsodHgu+hwWq9wxzOmXLBCYEq1Fyev4oc1+7hV/sEV7+0l3gMV772NFv3uD3c0Y8oNKwSm1Bq7cBv/nb2ZB3a9x2Xz9uITqD/kFZL7XhXuaMaUK1YITKm0eX8mj4xbznnp8+m5YD0b6ysXfzSBSvVbhDuaMeWOHSw2pU6WJ4+HPl9Ki/QVPPbT5xxMhlZ3XmhFwJggsS0CU+oM+mQh29ev4L3/fUxGPNS7ozUNbng/3LGMKbesEJhSZdLyXSxZtZl/zXubeA/E3Hsep9/5SbhjGVOu2a4hU2rsSM3msdELeH35q9Q6oKz7Y1NaDf443LGMKfesEJhSIeNIHrf/dz6vrnqFhtty+LF3Ilc+8jVY/8LGBJ3tGjJhl56Ty6BPFnLZ1jdosPoQs8+O4vJnx1gRMCZErBCYsMrJ9dLrnz/SbMdkuv24hmgVLnjynzRKaRLuaMZEDCsEJmxycr1c/vb/OHvtNO5aPJnD8YK++Qjt2titI4wJJTtGYMLC51PuHrmI05fM4KGfv+JwJUh87i7adf9juKMZE3GsEJiweO27Naxd/gv3rBiPVyD9mf60u/RP4Y5lTESyQmBC7qulO/jPzDW8tOh1onOVGXc34YpLnw13LGMilhUCE1I/bzjAn0Yv5s87X6DmXi/LLkjk3vsmInaGkDFhY4XAhMwXi7dz0wdzuCbnObouymB9c+GGf/+PmCg7Z8GYcLJCYEJi7MJtPPT5Ujrl/Zs7v8sgJ0G4bNhkYmMqhDuaMRHPCoEJuikrdvHIuOV0ZhpPTtzE/hSlzeiRxNRpHO5oxhjsOgITZDm5Xh7/cgWdM3/iof9NxRultB/cnZim7cIdzRjjsi0CE1SPjFtOxfTV3DNvAgk5kPziACrf/m64Yxlj/AS1EIhIbxFZIyLrReTRIsY3FJGZIrJERJaLyOXBzGNCa/HWQ/y8cA3//P49kjNh11960OCaZ8MdyxhTQNAKgYhEA0OBy4CWwEARaVmg2ZPA56raFhgA2FfFcuLXnYe567/zGDLvVZIzlV09KtB90DvhjmWMKUIwtwg6AOtVdaOqeoDRQMFexxVIcv9OBnYGMY8JkTyvj0H/nceDP79B7f05TO0aRfdXfgh3LGPMcQTzYHE9YJvf4+1AxwJtngWmisj9QGWgRxDzmBAZNX8r93z3Fm33bmNyJ7jlhXHExqeEO5Yx5jjCfbB4IPCRqtYHLgdGiEihTCIyWEQWisjCffv2hTykCZwnz8eqYe/Qbu8Gvm0nnHNrP2rUODPcsYwxxQhmIdgBNPB7XN8d5u8O4HMAVf0ZqAhULzghVR2uqu1VtX2NGjWCFNecCq/9ZzwD509m8WlCx8sa06vb8+GOZIwpQTALwQKgmYg0EZE4nIPBEwu02QpcAiAiZ+IUAvvKX0Z9uWQ7zcYNwRMDNTod5vxrPrZexowpA4JWCFQ1D7gP+A5YhXN20C8i8pyIXOk2+wswSESWAaOA21RVg5XJBM/2Q1n88ObztNqezpZ2Hnrd+Q0k2NabMWVBUK8sVtXJwOQCw572+/tX4IJgZjDBl3kkj7+9MYa/LvyRzXWUnv2vg7ptwh3LGBOgcB8sNuXAi5/N4cGvX8cTC5X6ViOx59/DHckY8ztYITAnZcqKXTT6/HUScrz8dJmProPHQ2zFcMcyxvwOVgjMCUvLzuXLNz7iok1r+aEN3Hr/KEisHe5YxpjfyQqBOWHDRv3IA3M/YXcK1OxZh6p17Y6ixpRFVgjMCRk2awMpo98hWuGrPsK1A8eFO5Ix5gRZITC/28zVexk1egZdtq5hyrnCvTe8hlSqGu5YxpgTZIXA/C7r96Zz18hF3LxxLADezsk0a943zKmMMSfDCoH5XZ74ciUtdq/k/C3bmdlB+ct1w8MdyRhzkqyrShOwpyasZN7GA7y/YRT7k6BV/95E1T4n3LGMMSfJtghMQN6btYERc7dw4/qJ1N+VzYJOUVx4+evhjmWMOQWsEJgSLd2WyitTVpPgyeKGdbPJiYXuN9wBUdHhjmaMOQWsEJhi5Xl93PnxQgDG7RpHrMfHN/1iaNP5wTAnM8acKlYITLFGzd/K/owjPNWpBr6ly1neWLi4z/9BlL11jCkvAjpYLCI1ce4SWhfIBlYCC1XVF8RsJsx2pGbz8pTVnFkniS47Z+FV+PX8KPq3vzfc0Ywxp1CxhUBELgYeBaoCS4C9OJ3HXA00FZFxwOuqejjYQU1o5Xl99H/vZ7I8Xl7o24IDN95Nak0YeP7ViHU2Y0y5UtIWweXAIFXdWnCEiMQAfYGewPggZDNh9Ndxy9mRms2gLk2o/NMYKqV6mXa1cll3u8W0MeVNsYVAVR8pZlweMOGUJzJh9791+/hiyQ4ual6Dx3o1Y17X/hysBXdc/zc7U8iYcqjEI34iEi0i1f0ex4nIYBFZFdxoJhy2Hczi/0YsonpCBd67qR3rPnydlINe9nSuRL1zbw93PGNMEBRbCERkAHAQWC4is0SkF7ARuAy4MQT5TAjl5Hq5f9QSsjxeht3cjtjsDNLf+4Q9KdDjpsfDHc8YEyQlHSN4EjhXVdeLSDvgZ+A6Vf06+NFMqA2duZ6l21IZ1KUJ5zaqyrwhD5OUrfzaN4FuZ/cLdzxjTJCUtGvIo6rrAVR1MbDOikD5tGrXYYbN2kjNxAo80aclnpwskj6YxJp60H/Q2+GOZ4wJopK2CGqKyEN+j1P8H6vqG8GJZUJJVXno82V4vD7G330+AOOHP0wboEoLoWLDzuENaIwJqpIKwftAYjGPTTnw6bytrNp1mD/3aE6DqpX4afv/iPt6JqlJyoX9B4U7njEmyEo6ffTvACJSXVX3hyaSCSWfT/ngp00A3NXtNHZn7mbKf57ipm2QdEEU0Rf9OcwJjTHBVtJZQ31FZB/OWUPbReT8EOUyITJ4xCI27c/k3oubcujIPvp9dS3dpu8lr0oede99HOwqYmPKvZIOFr8EdFHVukA/4OXgRzKhsml/JtNX7aFlnSQe7nUGQ5cOpfGaNOrvVxp0ro20tTOEjYkEJRWCPFVdDaCq87DjA+XKa9+uRgTeu+lcJm6YyIT1E7h1QQzRFb0k3T/EtgaMiRC/96yhmnbWUPkwafkupqzcTZ9z6rAn9xee+/k5ztufTIONB6jaDqKadAp3RGNMiNhZQxEo2+PllW+dO4Q8cWUDBk6+Do/Pw0PzUpAoJeWvVt+NiSQlFYIDqvqvkCQxITNh6Q62Hczm71c15fbvbuJgzkFebHw30SvepmrbOGJa9wl3RGNMCJV0jOCPIUlhQibzSB7PTPyFZjUT2MLn7MzcyS0tb+G8T74BoMrtd9uxAWMijPU3GGE++XkLnjwf7VotYezasVzS8BLu3RJH6pyNJLeqQmzPu8Md0RgTYiXtGmolIkX1PiaAqmpSEDKZIMn2ePngp01Uqb6WyTs+pFJMJV5o8zDbu/UkJt5LzTc+DXdEY0wYlFQIVqhq25AkMUE3cu4W9memU7vxF+R54aurvyL31T/hPQLV+/cmpt5p4Y5ojAkD2zUUIdKyc/nXzPVUafQFmd5U3uj2BrXia5I+awEA1Z95M8wJjTHhUlIhGBuSFCbo/jltLZnRS8mLX8olDS+hR8MeZI0fSuaeCtTs0xyJti4ojYlUJRWCaBGperyRItJdRPqe4kzmFFu9+zCfLVhPUv1viJIoXrvoNQD2Df+E6Ipeqjz9nzAnNMaEU0nHCJYDX4tIDrAY2AdUBJoBbYDpOPcjMqWUqvLI2OXEVptGrhzizW5vEhcdR9bsmWRvy6DWxUlEJdcId0xjTBiVdBvqr4CvRKQZcAFQBzgMjAQGq2p28COak/Hlkh2s2LmPxBaz6FKvC5c0vASAQ+85WwVJ9w0JZzxjTClQ0hYBAKq6DlgnIpVUNSvImcwpoqq8MW0tyY1G4QMGthgIQN7+vaQv3kRik2hizro4vCGNMWEX0FlDItJZRH4FVruPW4vIuwE8r7eIrBGR9SLy6HHa9BeRX0XkFxH57HelN8VasPkQO7PX4YtfRY+GPehSvwsAacNeRr1C9XvvDXNCY0xpENAWAfAmcCkwEUBVl4nIRcU9QUSigaFAT2A7sEBEJqrqr35tmgGPAReo6iERqXkCy2CO4+M5m4mrNosK0RX4+wV/B0C9XlInzyQm3kvF3neEOaExpjQI+DoCVd1WYJC3hKd0ANar6kZV9QCjgasKtBkEDFXVQ+489gaaxxRvR2o2U9YuIjZpBdc3v56kOOci8PQJY/AcOEKtq86CmAphTmmMKQ0CLQTb3G4qVURiReRhYFUJz6kH+BeP7e4wf82B5iIyW0TmikjvAPOYYqgqT365ggq1JwJwRdMrjo078P67xCXmknjXC+GKZ4wpZQItBHcB9+J8kO/AOXX0nlMw/xicU1G7AQOB90UkpWAjERksIgtFZOG+fftOwWzLt3d/2MCP22cTU2kzvRr1omW1lgDkrFhGzuYDpJwdh9Q+K8wpjTGlRaCF4AxVvVFVa6lqTVW9CTizhOfsABr4Pa7vDvO3HZioqrmquglYi1MY8lHV4araXlXb16hh57wXJ8uTx5Dv1pBSbwYAT3d++ti4A2+9hET7SLnr8XDFM8aUQoEWgncCHOZvAdBMRJqISBwwAPdgs58JOFsDiEh1nF1FGwPMZIowev42JG4fnugt9G7cm+QKyQB409I4PGc5Kc29RLe/IcwpjTGlSbFnDYlIZ+B8oEaBvouTgGJvTqOqeSJyH/Cd2/ZDVf1FRJ4DFqrqRHdcL/fUVC/wiKoeOPHFiWw5uV7e/d9SKjV8H4A7z7nz2LjDI94EHyT37g7RseGKaIwphUo6fTQOSHDb+fdVfBi4rqSJq+pkYHKBYU/7/a3AQ+6POUlTf91BVpX3iY49zJCuQzij6hnHxh3+ajxxiXlUvPnVMCY0xpRGJd1iYhYwS0Q+UtUtIcpkToDPpzwz+3miK2/j8fOepnfj307AyvnuA7K25VKj34VIpSphTGmMKY0CvaAsS0SGAGfh3HQOAFXtHpRU5nf7bvU6civNp2bMWQxseX2+cYc+fA9ESbnn6eM82xgTyQI9WPwpzu0lmgB/BzbjHAw2pcSQOZ8g4mNI96fyDc/bvoG0lYeJr1uRmHpNwpTOGFOaBVoIqqnqB0Cuqs5S1T8CtjVQSizYvIe9Mo2U6CacWyf/9QGpb/4V9UZR5xk7ZdQYU7RAdw3lur93iUgfYCdw3A5rTGj95YdHkZgsBrW+Od9wPZJB2qzlVKoVRYWL+ocpnTGmtAu0ELwgIsnAX3CuH0gCHgxaKhOw/23cxEFdSgzx3HJO/mMDOd8Mx5MeQ7Vb+oUpnTGmLAi0P4Jv3D/TgIsBROSCYIUygXt69jNIVB4vXfBWvuGqyr73RyLRSuIf/hymdMaYsqCkC8qigf449xj6VlVXun0UPw7EA22DH9Ecz/wdK9nvW0I17czlp+e/K3jGd5PI3JxN9Z7NiK5WLUwJjTFlQUlbBB/g3C9oPvC2iOwE2gOPquqEYIczx+fxenhwpvNN/4Uuj+Ub5/N42P3sM8TEe6n2p7+GI54xpgwpqRC0B1qpqk9EKgK7gaZ2G4jwG7f2S9K9u5G0HlzYNP9poYe/mkheahYN+lYi6vQuYUpojCkrSjp91KOqPgBVzQE2WhEIv50ZO3ltwSv48hJ45eLC+//TxnxMTLyXypdeHYZ0xpiypqQtghYistz9W4Cm7mPBuVVQq6CmM0Uas/pzvJqH7h1A77Ma5Bt3ZP16slaup6WrKa4AACAASURBVPpZWUgXO0hsjClZSYWgpD4HTIhleDL48JcP8GbX5/7OfYmKknzj9w8dSlSMjyqXdYSKyWFKaYwpS0q66ZzdaK6UGbZ8GACeAxdxU6eG+cZ5MzJInz6NlMbZxLS3aweMMYEJuPN6E36qyoR1X+HLTaHfGX1IrJi/X4G08Z+juV6Sz64MrQeGKaUxpqyxQlCGTNk0hVTPITz7LmFAhwaFxh8e/SEVUnKJv3UIiBQxBWOMKSzgQiAi8SJyRsktTTD41McLc18EoE21rrRtmL9fgSNrV5O96QBJTXzQok84IhpjyqiACoGIXAEsBb51H7cRkYL9D5sg+nbTt6TnHubI/m48dXnhC7pTh72GRCkpd/7ZtgaMMb9LoFsEzwIdgFQAVV2K0zeBCQGf+nh1wWsAtEnsR6v6KfnGq8dD2g/zqVTLQ0zX/wtHRGNMGRZoIchV1bQCw/RUhzFFG/nrSA7mHMBz8AIGnnd6ofGpoz7Bm+mlSrczISYuDAmNMWVZoLeh/kVE/gBEi0gz4AFgTvBimaNUlXHrxgEgB/twdZt6hdoc+u97xFbKI+H2Z0IdzxhTDgS6RXA/Tn/FR4DPcG5Hbf0RhMCiPYvYlLaJnD196NmyPlJg/3/e1rUc2Z1JSrsUpGGHMKU0xpRlgW4RtFDVJ4AnghnGFHZ0ayA39VyevrtlofFp//wTAAn97wtpLmNM+RHoFsHrIrJKRJ4XkbODmsgcsyF1A5M2TiI3rS0Xnd6IGokV8o33bVvBge83ULFeIhV63hKmlMaYsi6gQqCqF+P0TLYPGCYiK0TkyaAmM3zy6ycAHNl7KU/1KXzbp8yxb+M9Ek31+x8stMvIGGMCFfAFZaq6W1XfBu7Cuabg6aClMni8HqZtnoZktqZ9/SY0q5VYoEEmqZOnEx0vVL78uvCENMaUC4FeUHamiDwrIitwOq+fA9QParIIN2H9BNJz08k62IquzWsUGp83/Z9k7KhAcvfziIqzU0aNMScu0IPFHwJjgEtVdWcQ8xicU0ZHrR4NGk1eRgvu6Vb42oHDU6aBCsmDHw9DQmNMeRJQIVDVzsEOYn7z046fWJ+6Dk9ae567qlWhPgc0Yx+pC3ZQsV4VKp5ht38yxpycYguBiHyuqv3dXUL+VxJbD2VBdLTPAd1/BTd1bFRofM6nT3AkNZbaN10Z6mjGmHKopC2CP7m/+wY7iHGs2LeCZfuWcWRfd27vfGahrQF8Xg5NmgVRUSTd/EB4QhpjypViDxar6i73z3tUdYv/D3BP8ONFnufnPg9AbmpHBl10WqHxun4mGVuV+OaNiE62riiNMScv0NNHexYx7LJTGcTAtC3TWHVwFXmZTenU8DSqJ1Qo1Cbj/cfw5kRT5bbBYUhojCmPSjpGcDfON//TRGS536hEYHYwg0WisWvGApC97Vauu66Is3OzUzm0cD8xlSuReLntrTPGnBolHSP4DJgCvAw86jc8XVUPBi1VBFp/aD0/7/qZ5NyupGscV7SuW6hN7rR3yNxVgWrXdbVrB4wxp0xJhUBVdbOI3FtwhIhUtWJw6rw470WiiGbHpgu548ImxMUU2Gvn85E6YjiQQMqgv4UlozGmfApki6AvsAjn9FH/U1gUKHw00/xuy/ctZ+GehcR7WqPeRB65tPC1Abp1Hqnr46l8VgPiGhU+pdQYY05UsYVAVfu6v61byiB6a/FbAOzffjF9WtWhYmx0oTYZY/5FXnY0tW/5Y6jjGWPKuUDvNXSBiFR2/75JRN4QkYbBjRYZdmfuZv7u+TSt0APfkdrccWHRNffwj4uJqqAkXN4vxAmNMeVdoKeP/hvIEpHWwF+ADcCIoKWKIOPXjQdg17b2NKgaT7uGVQq18W5aTPpmJbljcyQ2NtQRjTHlXKCFIE9VFbgK+JeqDsU5hdScBJ/6eG/ZezSKb8ueAyk80L1Zke0yR76KeoWkP9wV4oTGmEgQaCFIF5HHgJuBSSISBZT41VREeovIGhFZLyKPFtOun4ioiLQPME+5MGnjJAB27TqN6glxXHduEdcOZO4ndeZCYpJiie9yaYgTGmMiQaCF4Aacjuv/qKq7cfoiGFLcE0QkGhiKcwVyS2CgiBTqdFdEEnHuaTTvd+Qu83Lycnhh7gtUrVCbA7va0eecOkX2Muad+Q6Zu2JJ6t0TiS58ENkYY05WoF1V7gY+BZJFpC+Qo6qflPC0DsB6Vd2oqh5gNM6upYKeB14FcgKPXfaNWTOGrLwsEjOvAaL5c8/mhRul7yZj/HBQIenam0Oe0RgTGQI9a6g/MB+4HugPzBORkvpHrAds83u83R3mP912QANVnRRw4nJi2b5lRBHFynUN6dCkKimVirhSeOGHZOyMIzqpMhVb2R2/jTHBEWgPZU8A56nqXgARqQFMB8ad6Izd4wxvALcF0HYwMBigYcOyf9bqgewDfL/1e+rFdCUN4cWrzy7cKCcN/eEfZO6pQ8KlPZCogLuXNsaY3yXQT5eoo0XAdSCA5+4AGvg9ru8OOyoROBv4QUQ2A52AiUUdMFbV4araXlXb16hRuP/esuaNRW/gVS97drQlLjqqcMf0AKsnkb0/Gm+Oj4RuXUMf0hgTMQLdIvhWRL4DRrmPbwAml/CcBUAzEWmCUwAGAH84OlJV04DqRx+LyA/Aw6q6MMBMZVKuL5fvt35Pg8rN+PVgTe7uVsQFZEfSYcI9HNraAImPpfIFF4Q+qDEmYgTaZ/EjInItcKE7aLiqflnCc/JE5D7gOyAa+FBVfxGR54CFqjrxZIKXVaNXjyYjN4O6ub0BGNyliNs1LR6Bz6tkbI8muW9fopOSQpzSGBNJSuqPoBnwD6ApsALnG/uO4p7jT1UnU2DLQVWfPk7bboFOt6zy+ryM/HUkteLrsmhxXS4+owZVKhdxkHj2W2TlNMGXfYSE7heHPqgxJqKUtJ//Q+AboB/OHUjfCXqicmzurrnszNxJQs6liAhP9DmzcKO1UyFjNxmpDZD4eCp37hz6oMaYiFLSrqFEVX3f/XuNiCwOdqDy7F9L/gXA0tWNuKZNPU6vWeAg8ZF0GHc7GptA+upUKl9wPlEVK4YhqTEmkpRUCCqKSFt+64cg3v+xqlphCNDaQ2tZeWAlcZIIvoo8elmLwo2WjQZPBhk1B5G3bxKJPXqEPqgxJuKUVAh24Zzrf9Ruv8cKdA9GqPJoyqYpABzZcg+XtKhJraQC3/S9ufDjPwA4vBGikpNJ7tMn1DGNMRGopI5p7EjlKTJ9y3TqVjyDNZlVGNChiIvilo+BjN1ojxfJ+PMnJPbsabecNsaEhF2uGgIbUzey+fBmvBnObSIuaVGzcCN3ayAz72x86ekk9uoZyojGmAhmhSAEjnY+s3FLfa5uU5eoqAJ3Gd25FA5tgoufJH3GDKIqV7azhYwxIWOFIASmbplKrMTj89Tirm5NCzeY7fRZrK3+QPr0GSR07UpUhQohTmmMiVSB3n1U3L6Kn3YfNxSRDsGNVj7szNjJ7szdaEZbmtVMoEXtAlcJ71sLv3wBtc4hY8lavIcOkdirV3jCGmMiUqBbBO8CnYGB7uN0nE5nTAmGLR8GQOquC+l1Vq3CDaY+4fzu/zHpU6cSlZREol1NbIwJoUALQUdVvRe38xhVPQQUcW8E4y/dk84X674gKboemluV685tkL9Bngc2zIRqp6PJjciYOZOErl2ROFu1xpjQCbQQ5LpdTyoc64/AF7RU5cTEDc599aIP9iM+NprG1Srlb7D6a/DlQvenyFq4EG9aGomXXBKGpMaYSBZoIXgb+BKoKSIvAj8BLwUtVTkxdfNUoiWarTvrcMeFTQr3SbzqG+d3s16kff01UZUrk9D1otAHNcZEtEBvQ/2piCwCLsG5vcTVqroqqMnKuM1pm1m8dzHNY/uzWIQBHQrsFvJ5Yd1UOLsfPo0m/bupJPboQVR8fHgCG2MiVkCFQEQaAlnA1/7DVHVrsIKVdWPWjAFgyZq6dD6tGvWrFNgttH4GeDKgWS8yf/oJX0YGSX37hiGpMSbSBdpD2SSc4wMCVASaAGuAs4KUq0zzqY/vNn9HFLH4jtTmteuK6Hh+/nCQKGjRh8OPPUt0SgqVO3UMfVhjTMQLdNfQOf6PRaQdcE9QEpUDUzZNYV/2PnJ396ddw5TCWwN7foH106DNTfh8MaT/8APJffvavYWMMWFxQlcWu7eftq+vxzFs+TAqRaeQc6g1bRpUKdxg0cfO74seJvPnuWhWFkm9Lw1tSGOMcQV6jOAhv4dRQDtgZ1ASlXG7M3ezKW0TeYdbUz2hEk8W7IUszwPzh0FiXajahMyfRiCVKhHfvn14AhtjIl6gxwj8u9LKwzlmMP7Uxyn7JqyfAIDnYGdevbpl4RvMrRjr/O72N1SVjFk/UrljR6LsIjJjTJiUWAjcC8kSVfXhEOQp877Z+A2+3CSic0/jilZ1CjdYPcn53eZGjqxeTe6OHVS/+67QhjTGGD/FHiMQkRhV9QIXhChPmbbqwCq2HN5CXvrZPND99MIXkHmyYMtP0OhCiI7l8NSpEBVFgl1NbIwJo5K2CObjHA9YKiITgbFA5tGRqvpFELOVKarKC/NeABU8+y/m1vMbF2608EPISYPz7wcgfdo0Kp13HjFVijigbIwxIRLoMYKKwAGcPoqPXk+ggBUC14r9K1i+bzk5e/vQuEptEisWOBXU54XZb0KVJtD8Uo5s3Ihn/QaqDBxY9ASNMSZESioENd0zhlbyWwE4SoOWqgz6dNWnxEgcuWntuPPKJoUbbJoFmfvg4idAhPSp0wBI7NEjxEmNMSa/kgpBNJBA/gJwlBUCl099/LTjJ+LzWnDIW5l+7eoXbjRrCFSqBq0HAM5uofjWrYmtVUQfBcYYE0IlFYJdqvpcSJKUYYv2LOKw5zDZu1tw+Tm1qRgbnb/B2u9g6xy45BmIjcezfQc5v/xCzUceCU9gY4zxU9KVxUVtCZgCftrxEwDezOZc1aZe4Qaz33Z+H90amO7uFurVMyT5jDGmOCUVAjuvsQR5vjy+3fQtCdoM9SbQ48wCu3p2LXNOGW17MyTVBSB96jQqnHkmcQ0aFDFFY4wJrWILgaoeDFWQsmrKpinszNzJvh0d6HNOHaILXkk8/33nd9e/ApC7dy/ZS5aQ2NMOEhtjSocTuumc+c2cnXMAyMs4k8EXnZZ/ZM5h55YSZ10LKQ0ByJgxA1RJ6tUr1FGNMaZIVghOQlZuFtO3TEfTz+X0Gim0bpCSv8GsVyEvB9rdcmxQ+vQZxDVuTFzTpiFOa4wxRbNCcBJ+3PEjOd4csg+25W+9WxRusGYyVK4JTS8GwJeZSdaCBSR07Vr49hPGGBMmVghOwie/fAIaS5SnMV2aVc8/cvcKOLgRWvU/Nijt629Qj4dE63vAGFOKWCE4Qfuy9rFi/wryMk7nhvanFb52YOko5/e5tx8blPrFF1Ro3pz4Nm1CmNQYY4pnheAETdk0BYAj+y9hwHkN84/0+WDpp3B6D6h+OgBZi5eQs3w5Kddfb7uFjDGlihWCE/S/HbNRXzSnJ7fg7HrJ+UeumQQ5qdC0+7FBhydNQipUIOXaa0Kc1BhjimeF4ATszNjJ3F1zyD10Ade0LeK+QnP/DfFV4LxBAKjPR/r06VTuciFRlSuHOK0xxhTPCsEJGLvW6W6yJl24s0uBO40e3glbZkOTrhDjdD+Zs3IleXv22J1GjTGlkhWCE/Djtp/weapzTq1mxEYXWIWz33J+X/CnY4PSp02HmBgSu3ULXUhjjAlQUAuBiPQWkTUisl5EHi1i/EMi8quILBeRGSLSKJh5ToXt6dtZm7qa3LQ23NSpQNz0PbDoI2h4PtRrBzg9l6VPm0blDucRnZJSeILGGBNmQSsEbqf3Q4HLgJbAQBFpWaDZEqC9qrYCxgGvBSvPqTJ69WgAakd35vymBa4dWDLCuZK490vHBnk2bMCzeTMJtlvIGFNKBXOLoAOwXlU3qqoHGA1c5d9AVWeqapb7cC5QxJHX0sPr8/L1hinkZTRjUKcO+Ufm5sC8Yc6xgbptjw0+PNk5zTTxEisExpjSKZiFoB6wze/xdnfY8dwBTAlinpO2Yv8KDh7ZS+7h1lx2Tp38I2e/BZl74YIHjg3K27ePgyNHktCtG7G1aoY4rTHGBCbQzuuDSkRuAtoDXY8zfjAwGKBhw4ZFNQmJsWvHoRrFGUkdqZ5Q4bcRuTnw42uQVB+a/taFw8ERI/FlZFDzr9YTmTGm9ArmFsEOwL/nlfrusHxEpAfwBHClqh4pakKqOlxV26tq+xo1agQlbEnSjqTxzYZJeDPO4JLmBW43Pfdd8OVBr+fBvWpYc3NJ/eILErp1o8JppxUxRWOMKR2CWQgWAM1EpImIxAEDgIn+DUSkLTAMpwjsDWKWkzZ7x2x85FFdu/JQz+a/jfB5Yd57UONMOOu3q4bTv5+Jd/9+UvpfH4a0xhgTuKAVAlXNA+4DvgNWAZ+r6i8i8pyIXOk2GwIkAGNFZKmITDzO5MJu3KoZqDeeQe0vyz9i5XjI2AMXPnhsawAg9fPPialTh4QuXUKc1Bhjfp+gHiNQ1cnA5ALDnvb7u0ycSpN2JI3F+38gL7MFV7Xx29t1JAOmPgkJtaDl1ccG5+7dS+bs2VS/914kOrqIKRpjTOlRKg4Wl3ZfrPsSLx7qR19CcnzsbyNmvepsDdz+LcRWPDY4a/4CABK6FXns2xhjShW7xUQJPF4Pw5d9gC83iTvb/3ZGEAc3OQeJ67aDRp3zPefQqFHE1KlDxTPPDHFaY4z5/awQlGDOzjlk5KVyZE8frmnnd73bj0OcA8X9/pOvffaKlWQvWkS1229DYmyDyxhT+lkhKMF7y4ahvljOr9uVCjHu/v7tC52OZzoMhmr5O6FPHTcOqViR5Gus3wFjTNlghaAYOXk5rDu0Ad+RWlx2lnuQWBW+fczpb6Bb/vvoqcfD4SlTSOzVk+jExDAkNsaY388KQTGGLx+Ox5dNxfQruKG9WwjWz4Dt8+Hc26BS1XztM+fNx3f4MEm9e4c+rDHGnCArBMVYtGsVADecfTFRUe41AnOHOr8vKnzbiNSxY4lKSqLy+eeHKqIxxpw0KwTHkXYkjaX75+E51IHrj24NZO6HjbOgyUUQl7/LSc/27aRPn06VG/oTVbFiEVM0xpjSyQrBcYxf+wU+cmlS4WJOq5HgDJz5EqgXLn25UPtDI0ZAVBRVbrwxxEmNMebkWCE4jlGrvsCb3YAbWl3gDDiSDss/d7YGap+dr603PZ3UseNIuuwyYmvXDkNaY4w5cVYIirAtfRu7szcTldWaGzu63VEuHgGedOj+dKH2qePG48vKouqtt4Y4qTHGnDwrBEV4f4nTHWXb6ucTF+OuopXjICYeGpyXr616PBz88EMqdehA/NlnhTqqMcacNCsERfhu81R8eQm82Me9pUTaDtixCNrdXKjt4WnTyNu3j2p33hHilMYYc2pYIShgc9pWsnQ3lTzn0bBaJWfgjOec3x0GF2p/6LNRxDZsSOULLwxhSmOMOXWsEBTw3oLxAPQ7fYAzIGMv/PKlc3O56s3ytc1ZvZrsRYuoMmAAEmWr0hhTNtld0QqYuW0mXk8d7ruogzNgxnPgPQJ93yjU9sD7/yGqUiVSrrX7Chljyi77Gutn46EdZEVtomHlFiRUiIE8j3NzuQYdoW7bfG2PbNrE4SlTqPKHgUSnpIQpsTHGnDwrBH5em/0xANc1d/sZXvopqA/O7leo7f6h7yJxcVS97bYQJjTGmFPPCoHLpz5+3vstPk9Vbm7n3ito1dfO7/b5zwjKWbWKw998Q9VbbiGmevUQJzXGmFPLCoFr7C8z8EUfolXCdU6/Awc2wIbvnbuMRuc/lHJw5EikUiU7ZdQYUy5YIXB9tvJb1BfLy71ucQbMfgtQuODBfO28aWkcnjyF5D6XE52UFPqgxhhzitlZQ4CqsilzIRV8Z9C4WrJzkHjJCKh1DlRtkq/t/nffRXNyqHLTTWFKa4wxp5ZtEQCfr5iDxhykXXX32MCvE5yDxJ3uztcue8UKDo78lJTrr6fiGWeEIakxxpx6VgiA95f/F4CnLx7oDFg5HqJioNUNx9r4cnLY+fAjxNSsSc2H/xKOmMYYExQRv2soJy+H3Xnzife2oEFKdUjdCmu/hU735jtInDbhKzxbttDgP/+xYwPGmHIl4rcIXp8zChGlZ71rnQGLPnJ+t/nDsTaqysGRI6jQ8kwqX2DdUBpjypeILwTfbvoO9UXzZPfr4OAm+N8b0KxXvs5nMufMwbN+A1VvvgURCWNaY4w59SK6EBzKyuSQbw1Vo86mUkwUjHFvM9310fztRowkulo1kvpcHoaUxhgTXBFdCIYv/AaJyqNP42thyUjYswI6DIL65x5rk7NqFRmzZlHlhhuIiosLY1pjjAmOiD5YPH3jAhC4vXVnGNoOYitBz+ePjVdV9rz6GtHJyVS9zbqhNMaUTxG7RZDr9bHTs4j4vJbUnD8UcrPgmvcgtuKxNhkzfyBr7lyq33efnSlkjCm3IrYQTPzlF6Iq7KdTzbaw8L9QvTm0vOrYeM3NZe+QIcQ1aUKVG/qHMakxxgRXxBaCoXO/AuB29oAnA/q+mW98+rRpeDZtosZDf0ZiY8MR0RhjQiIiC0Gu18du31zitBZtt/wEcYnQ6LfrA7zp6ez9x+vENW5MYvfuYUxqjDHBF5GFYNKqpURV3ErP2Pqwfy10+TP4XR+w+/nnyd2zh7qvvIxER4cxqTHGBF9EFoIRK75CRHlo40RnQKd7j41LmzSJwxO/pvrddxPfpk2YEhpjTOhEXCHw+Xyszf6Oup4q1PR64fqPj50plLt3L7v//hwVW7ei+l3/F+akxhgTGhF3HcHYFQsgOpMemYcgvuqxM4VUlT0vvIjm5FD3lVeQmIhbNcaYCBVxWwRDF39AjE+489BeGDgKRFCPh51/+xvpU6dS/b77qNCkSckTMsaYciKiCkGu18tBFnK2x0MViYWGnQDYM+QfznGB+++zfoiNMREnovZ//HPOWCQql2vSD8OVbwPOweFDI0ZQ9dZbqHHvvSVMwRhjyp+gbhGISG8RWSMi60Xk0SLGVxCRMe74eSLSOJh5Jm2aSCWvcGVGJr7mV7L39TfY+bdHiW/blpoPPxzMWRtjTKkVtEIgItHAUOAyoCUwUERaFmh2B3BIVU8H/gm8Gqw8h3OyOKgruDYjDZpex7bBd3Hg/fdJvuIK6r871K4eNsZErGDuGuoArFfVjQAiMhq4CvjVr81VwLPu3+OAf4mIqKqe6jBTp39G+7U+LtwkbFi1GM3zUvcf/yC5b59TPStjjClTglkI6gHb/B5vBzoer42q5olIGlAN2O/fSEQGA4MBGjZseEJhsr7+iL/O8gHxVOrWiVqP/o24xo1PaFrGGFOelImDxao6HBgO0L59+xPaWrjmqRH8eP5Iel06iNjatU9pPmOMKcuCWQh2AA38Htd3hxXVZruIxADJwIFghEms34Q+tz4VjEkbY0yZFsyzhhYAzUSkiYjEAQOAiQXaTASOdv11HfB9MI4PGGOMOb6gbRG4+/zvA74DooEPVfUXEXkOWKiqE4EPgBEish44iFMsjDHGhFBQjxGo6mRgcoFhT/v9nQNcH8wMxhhjihdRt5gwxhhTmBUCY4yJcFYIjDEmwlkhMMaYCGeFwBhjIpwVAmOMiXBWCIwxJsJJWbuQV0T2AVtO8OnVKXBDuwhgyxwZbJkjw8kscyNVrVHUiDJXCE6GiCxU1fbhzhFKtsyRwZY5MgRrmW3XkDHGRDgrBMYYE+EirRAMD3eAMLBljgy2zJEhKMscUccIjDHGFBZpWwTGGGMKKJeFQER6i8gaEVkvIo8WMb6CiIxxx88TkcahT3lqBbDMD4nIryKyXERmiEijcOQ8lUpaZr92/URERaTMn2ESyDKLSH/3tf5FRD4LdcZTLYD3dkMRmSkiS9z39+XhyHmqiMiHIrJXRFYeZ7yIyNvu+lguIu1OeqaqWq5+cDrB2QCcBsQBy4CWBdrcA7zn/j0AGBPu3CFY5ouBSu7fd0fCMrvtEoEfgblA+3DnDsHr3AxYAlRxH9cMd+4QLPNw4G7375bA5nDnPsllvghoB6w8zvjLgSmAAJ2AeSc7z/K4RdABWK+qG1XVA4wGrirQ5irgY/fvccAlIiIhzHiqlbjMqjpTVbPch3Nx+pAuywJ5nQGeB14FckIZLkgCWeZBwFBVPQSgqntDnPFUC2SZFUhy/04GdoYw3ymnqj/i9Nh4PFcBn6hjLpAiInVOZp7lsRDUA7b5Pd7uDiuyjarmAWlAtZCkC45AltnfHTjfKMqyEpfZ3WRuoKqTQhksiAJ5nZsDzUVktojMFZHeIUsXHIEs87PATSKyHadHxPtDEy1sfu//e4mC2lWlKX1E5CagPdA13FmCSUSigDeA28IcJdRicHYPdcPZ6vtRRM5R1dSwpgqugcBHqvq6iHTG6Qf9bFX1hTtYWVEetwh2AA38Htd3hxXZRkRicDYnD4QkXXAEssyISA/gCeBKVT3y/+2da4xdVRXHf3+xWGgNmtYHmGgbAZUgDG8TUq3WjAaksaFxRBIsIcYYwdgCidEKpEF8TKvBT2pLMyQSrBCLI0RHEJqOUm3rdGhL0YIWiU/4QBoKNSn498NeV25O78zccV69965fMpl7zt177bXPndnr7LXuWWuadJsqxprz64EzgS2Snqb4UvtbPGDczOf8V6Df9hHbB4D9FMPQqjQz52uAHwPY3gbMpuTkaVea+n8fD+1oCHYAp0laKOl4SjC4v9KmH/h0vF4OPOyIwrQoY85Z0jnA9ylGoNX9xjDGnG0ftD3f9gLbCyhxkaW2d86MupNCM3/b91F2A0iaT3EV/Xk6lZxkmpnzM8ASAEnvoRiC56ZVy+mlH7gqvj30PuCg7X9MV9Of5gAABZBJREFURGDbuYZsvyzpWmCA8o2DjbYfl7QG2Gm7H7iDsn18ihKU+eTMaTxxmpxzLzAXuCfi4s/YXjpjSk+QJufcVjQ55wGgW9I+4BXgRtstu9ttcs7XA+slraQEjle08o2dpLspxnx+xD1uBmYB2P4eJQ5yCfAU8BJw9YTHbOHrlSRJkkwC7egaSpIkScZBGoIkSZIOJw1BkiRJh5OGIEmSpMNJQ5AkSdLhpCFIkiTpcNIQJFOGpFckDdf9LBil7aFJGK9P0oEYayjSDYxXxgZJZ8TrL1fee3SiOoac2nXZK+lnkt4wRvuu/ye1sqSTJd0frxdLOlj3WTwU52+R9Lc6fZY2OL9P0hV1ctdK+tB49UmOXfI5gmTKkHTI9tzJbjuKjD7gftv3SuoG1to+awLyJqzTWHIl3Qnst/21UdqvoKTQvnac4/QCv7b9U0mLgRtsf6zS5hbgkO218VTuIPBm4Ka686cBvwfm2T6iUstive3u8eiTHLvkjiCZNiTNVSmKMyRpj6Sj0kbHXezWujvURXG+W9K26HuPpLEW6K3AqdF3VcjaK+mLcW6OpAckPRbne+L8FknnS/oGcELocVe8dyh+/0jSpXU690laLuk4Sb2SdqgUDPlsE5dlG5E5UtKFMcddkh6V9K5Iq7AG6AldekL3jZK2R9tG6bcBLgd+0YQOANh+AniZSp4e209SnmB9Yxz/BZgn6a3Nyk6ObdIQJFNJbSEdlrSZUhNgme1zKYVy1klH1YH4FDBguws4GxhWyZmzGvhw9N0JrBpj7MuAPZLOozyCfxEl8dxnVPIufRT4u+2zbZ9JZcG0/SXgsO0u21dWZG8CPgEQC/US4AFK8rODti8ALoixFo6koKTjom8tHcYfgEW2z6Hckd8WOfhvohQS6rK9iZI48GHbF8Z17JU0pyJ7IfB8JbngorrP4ysN9LkI+A+VPD0q6byfrOSoGgIuHmluSWvRdrmGkmOKw7GgAyBpFnCbpPdTFpy3AW8B/lnXZwewMdreZ3tY0gcolad+E3bjeMqddCN6Ja2mLGbXUBbazbZfDB1+AiyiLPzrJH2T4k4aHMe8fg7cLul1FIOy1fbhcEedJWl5tDuJkvnzQKX/CZKGY/5PAA/Wtb8zXDEm8ss0oBtYKumGOJ4NvD1k1TiZoxOvDVZdQ8FKlfTkLwA9th3XeaWkqymJ6y6r9HkWOGUE/ZIWIw1BMp1cCbwJOC98zU9TFrH/YXtrGIpLgT5J3waeBx60fUVVYANutH1v7UDSkkaNbO+PO91LgFsl/cr2mmYmYfvfkrYAHwF6KFWzoJQOvM72wBgiDtvuknQiJZna54HvUqqpPWJ7mUpgfcsI/QVcbvuPo41B5dqOwndsrx3pfASQ75D0Ttu1Sm+zY4ykDUjXUDKdnAQ8G0bgg8A7qg0iEPkv2+uBDZTarb8FLpZU8/nPkXR6k2MOAh+XdGK4T5YBg5JOAV6y/UNKZtZGBcCPxM6kEZsoLqfa7gLKov65Wh9Jp1ddNvVE6dAvANfr1boYtbzyK+qavkCpr1BjALiu5lYLV1eV/cCCkcYeD5Hhcyevpm6HsktoWFw9aT3SECTTyV3A+ZL2AFdRfOJVFgOPSdpFudu+3fZzlIXxbkm7KW6hdzczoO0hoA/YDvwO2GB7F/BeYHu4aG4Gbm3Q/QfA7lqwuMIvKVXeHgo/PhTDtQ8YkrSXUv9h1F136LKbUmXrW8DXY+71/R4BzqgFiyk7h1mh2+NxXJX7IvCnmvGcBNYAqyS9JgzdqRTjkLQB+fXRJGlTJC2juOFWT4Hcc21/dTLlJjNHxgiSpE2xvVnSvCkQ/Vpg3RTITWaI3BEkSZJ0OBkjSJIk6XDSECRJknQ4aQiSJEk6nDQESZIkHU4agiRJkg7nv7VxeYAfsiUlAAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ]
  },
  {
   "cell_type": "code",
   "source": [
    "print(f'ATTENTION : (AUC={roc:.4f})')\n",
    "print(f'GRU : (AUC={roc1:.4f})')\n",
    "print(f'LSTM : (AUC={roc2:.4f})')\n",
    "print(f'RNN : (AUC={roc3:.4f})')"
   ],
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "PmHmQRjtZI0a",
    "outputId": "d14dd8ab-b622-4d3a-899a-7ba63a598e84"
   },
   "execution_count": null,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "ATTENTION : (AUC=0.8212)\n",
      "GRU : (AUC=0.8062)\n",
      "LSTM : (AUC=0.8202)\n",
      "RNN : (AUC=0.8003)\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "NOTE: This is synthetic data thus the predictive accuracy should not be great however, aimilar models were tested and yield benchmark accuracy on real data is this paper: https://journals.plos.org/plosone/article/authors?id=10.1371/journal.pone.0245177\n"
   ],
   "metadata": {
    "id": "VmM3abvxzEeC"
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Interpret the output of attention model"
   ],
   "metadata": {
    "id": "ICJTt5J-g6Or"
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Two key important advantages of using attention based models:\n",
    "\n",
    "1. Just like in logistic regression we can examine the coefficients in front of each variable to assess their importance in the predictive task, we can examine the attention matrix to assess the relationship and importance between time steps. This is commonly done in Natural Language Processing task to see which word are important to consider to complete a sentence. For example, below to predict the word in red, we are paying more 'attention' to word in blue.\n",
    "\n",
    "![picture](https://drive.google.com/uc?id=1aYNBKjCwyfIuxYbxK3LSDjtmpSvb-XlL)\n",
    "\n",
    "2. The context vector described above allows us to output Heart Failure as a function of time, in the form of a visualized trajectory showed below. The context vector is a deep representation of a give patient medical history with repect to the target (Heart Failure). In a lot of NLP, problems the context vector is use to help computers understand paragraphs for many problems, such as translation, sentence completion, Ambiguity, speech recognotion and many others.\n",
    "\n"
   ],
   "metadata": {
    "id": "O-fv4oaGzd7o"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "def graphs(model,patient,X_test,y_test):\n",
    "  length_seq = X_test[patient].shape[0]\n",
    "\n",
    "  #### Graph Actual HF history ######\n",
    "  x = []\n",
    "  value = 0\n",
    "  for i in range(length_seq):\n",
    "    x.append(value)\n",
    "    value = value+0.5\n",
    "  plt.plot(x,y_test[patient])\n",
    "  plt.xlabel('Patient Age')\n",
    "  plt.ylabel('HF Occurence')\n",
    "  plt.title('Actual HF Trajectory')\n",
    "  plt.show()\n",
    "\n",
    "\n",
    "  #### Estimated Acutal HF history ######\n",
    "  one = model.predict(X_test[patient].reshape((1,length_seq,20)))\n",
    "  #y = (one - np.min(one)) / (np.max(one) - np.min(one)) # if you want to scale to 0,1 probability\n",
    "  x = []\n",
    "  value = 0\n",
    "  for i in range(length_seq):\n",
    "    x.append(value)\n",
    "    value = value+0.5\n",
    "  plt.plot(x,one[0])\n",
    "  plt.xlabel('Patient Age')\n",
    "  plt.ylabel('Estimated Risk')\n",
    "  plt.title('Current HF Trajectory')\n",
    "  plt.show()\n",
    "\n",
    "  #The graph shows that the predicted risk peaks around the actual HF event\n",
    "  #Note: we only have data for each patient after the age of 40 that is why it is at zero before 40\n",
    "\n",
    "patient = 0 #pick a patient from test set\n",
    "graphs(attention,patient,X_test,y_test)"
   ],
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 573
    },
    "id": "cuF7FZX5hA55",
    "outputId": "e97c1ef1-8d64-454f-c588-cb15bb13a331"
   },
   "execution_count": null,
   "outputs": [
    {
     "output_type": "display_data",
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ],
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAckUlEQVR4nO3de5hcdZ3n8fenO9zCJQETgVwgGQlonJGLeRhGwM0IMgkj4JUhgqKD5mEUFRUdRllAxnWHmV0YFRw2OAIiAuEim3ECLCK66ABDBxRJkDEThCTcOtwXhZCu7/5xftV9UnRXV5o+VUn/Pq/nqYc65/zqnG8Vlfr273oUEZiZWb66Oh2AmZl1lhOBmVnmnAjMzDLnRGBmljknAjOzzDkRmJllzonAsiDpbEnf63Qcm0LSlyR9u9Nx2NjnRGBtIeknkp6RtE2L5T8i6WdVx5WuNVfSmkH2/0TSx0plapL+X+nxL4O85sbS8VckrS9tX7QpcUXE1yLiYyN/Z0O/N7OycZ0OwMY+STOAQ4HngKOBazoZz2vwaERMa1YgIubXn0u6FFgTEWc0lpM0LiI2jH6Io2tLidNeG9cIrB0+DNwJXAqcWD4gabqk6yX1SnpK0gWS3gRcBPxJ+kv62VS2/y/0tL1RrUHS1yWtlvS8pGWSDm3Hm2uVpJD0SUm/AX6T9g0Zc2NzlqSDJP2bpGcl/VLS3NKxXSRdIunRVPO6QdL2wI3AlFKtZIqkbST9Yyr7aHq+TTrPXElrJP21pMeBSyTdL+mo0rW2krRO0v6Vf2jWFk4E1g4fBq5Ijz+TtCuApG7gh8DDwAxgKnBVRDwAnAzcERE7RMTEFq9zN7AfsAvwfeAaSduO5hsZBe8G/hiYnbZbilnSVOBfga+msqcB10manIpcDowH3gy8Hjg/Il4E5lPUZHZIj0eBLwMHpevuCxwIlGstu6Vr7AksBL4LnFA6fiTwWETc+xo+B9uMOBFYpSQdQvGDsjgilgH/CXwwHT4QmAJ8ISJejIiXImLE/QIR8b2IeCoiNkTE/wS2AfZp8eVT0l/a/Q/gkGHKHDuCMP97RDwdEb/fxJhPAJZGxNKIqEXELUAPcKSk3Sl+8E+OiGci4pWI+GmTGI4HzomIJyOiF/gK8KHS8RpwVkS8nOL8XrrOTun4hygSj40RTgRWtROB/xMR69L29xloHpoOPDxabdCSTpP0gKTn0g/5BGBSiy9/NCImlh9AY1JqLLN4BGGuHmHMewIfGCRR7U7xOT4dEc+0GMMUilpY3cNpX11vRLxU30i1iJ8D75M0kSLpXNHitWwL4M5iq4yk7YBjge7U3gzFX7wTJe1L8aO4xxAdkoMti/siRfNH3W6lax0KfBE4DFgeETVJzwAanXczavrf1ybGvBq4PCI+3ngg1Qh2kTQxIp4d6nolj1IkluVpe4+0r9lrLgM+RvGbcUdErB3szdmWyTUCq9K7gT6K9vD90uNNwO0U/Qb/DjwG/J2k7SVtK+ng9NongGmSti6d7xfAeyWNl7QXcFLp2I7ABqAXGCfpTGAnNm+bEvP3gKMk/Zmk7vRZzZU0LSIeo+gU/paknVNn7tvT654AXidpQulcVwJnSJosaRJwZjp/MzcABwCfoegzsDHEicCqdCJwSUQ8EhGP1x/ABRTt1AKOAvYCHgHWAH+RXvtjir9YH5dUb1Y6H1hP8eN2GRs3T9wM3AT8B0VTx0s0NMNshlqOOSJWA8cAX6JIHKuBLzDwb/hDwCvAr4EngVPT635N8cO/KjUpTaHocO4B7gN+BdyT9g0p9RVcB8wErh/Ru7XNlnxjGrPNk6RzgGkR8ZedjgUg1Vj2jogThi1sWxTXCMw2Q5JE0aT2UKdjgWKeAkVT3KJOx2Kjz4nAbPN0DzANuLjTgUj6OEVT1I0R8X87HY+NPjcNmZllzjUCM7PMbXHzCCZNmhQzZszodBhmZluUZcuWrYuIyYMd2+ISwYwZM+jp6el0GGZmWxRJDw91zE1DZmaZcyIwM8ucE4GZWeacCMzMMudEYGaWucoSgaTvSHpS0v1DHJekb0haKek+SQdUFYuZmQ2tyhrBpcC8JsfnA7PSYyHwTxXGYmZmQ6gsEaQ1SZ5uUuQY4LtRuJPiZiW7VxWPmdlgnn5xPTf+6rFOh9FRnewjmMrGa6+vSfteRdJCST2Senp7e9sSnJnl4YZ71/JXV9zDiy+Pyh1Tt0hbRGdxRCyKiDkRMWfy5EFnSJuZjcgrfTUANtTyXYCzk4lgLcVNt+umpX1mZm1T//3PeSXmTiaCJcCH0+ihg4Dn0r1XzczappYSQF/GNYLKFp2TdCUwF5gkaQ1wFrAVQERcBCwFjgRWAr8DPlpVLGZmQ6mlBJBxHqguEUTEgmGOB/DJqq5vZtYKNw1tIZ3FZmZVqTcN5VwjcCIws6xFfyLINxM4EZhZ1uo1AScCM7NM1RNAxnnAicDM8tbn4aNOBGaWt3DTkBOBmeXN8wicCMwsc55H4ERgZpnzPAInAjPLnOcROBGYWeY8j8CJwMwy1+d5BE4EZpa38DwCJwIzy1utlv6bcZXAicDMsuZRQ04EZpY5zyNwIjCzzIVrBE4EZpa3mucROBGYWd76PI/AicDM8tZfI6h1OJAOciIws6x5iQknAjPLnOcROBGYWeZ8q0onAjPLnBedcyIws8x5HoETgZllrs+dxU4EZpa3/qahjKsETgRmljU3DTkRmFnmvMSEE4GZZc7zCJwIzCxznkdQcSKQNE/Sg5JWSjp9kON7SLpN0r2S7pN0ZJXxmJk1Cs8jqC4RSOoGLgTmA7OBBZJmNxQ7A1gcEfsDxwHfqioeM7PB9LmzuNIawYHAyohYFRHrgauAYxrKBLBTej4BeLTCeMzMXmVg9dF8M0GViWAqsLq0vSbtKzsbOEHSGmAp8KnBTiRpoaQeST29vb1VxGpmmfISE53vLF4AXBoR04AjgcslvSqmiFgUEXMiYs7kyZPbHqSZjV2eR1BtIlgLTC9tT0v7yk4CFgNExB3AtsCkCmMyM9uI5xFUmwjuBmZJmilpa4rO4CUNZR4BDgOQ9CaKROC2HzNrm/o8gnAiGH0RsQE4BbgZeIBidNBySedIOjoV+zzwcUm/BK4EPhI5/98ws7aruWmIcVWePCKWUnQCl/edWXq+Aji4yhjMzJpx01DnO4vNzDpqYNRQZ+PoJCcCM8ua5xE4EZhZ5rzEhBOBmWXOncVOBGaWOXcWOxGYWeY8j8CJwMwy56YhJwIzy5ybhpwIzCxz/fMIMq4SOBGYWda8+qgTgZllzvcjcCIws8y5s9iJwMwyV+8b8PBRM7NMuWnIicDMMuemIScCM8tcPRH0ZZwJnAjMLGv133/3EZiZZcrzCJwIzCxz7ix2IjCzzLmz2InAzDIWEf13KHMfQROS9pZ0q6T70/ZbJJ1RfWhmZtUq1wLcNNTcxcDfAK8ARMR9wHFVBmVm1g7lH/++WgcD6bBWEsH4iPj3hn0bqgjGzKydyonATUPNrZP0BiAAJL0feKzSqMzM2iDcNATAuBbKfBJYBLxR0lrgIeCESqMyM2uD8o9/zqOGhk0EEbEKOFzS9kBXRLxQfVhmZtVzZ3GhlVFDX5M0MSJejIgXJO0s6avtCM7MrErl9YUyzgMt9RHMj4hn6xsR8QxwZHUhmZm1R2zUNJRvJmglEXRL2qa+IWk7YJsm5c3MtghuGiq0kgiuAG6VdJKkk4BbgMtaObmkeZIelLRS0ulDlDlW0gpJyyV9v/XQzcxeG88jKLTSWXyupPuAw9Kuv42Im4d7naRu4ELgncAa4G5JSyJiRanMLIrJagdHxDOSXj+SN2FmNhKeR1BoZfgoEXEjcOMmnvtAYGUadYSkq4BjgBWlMh8HLkz9DkTEk5t4DTOzEfM8gkIro4beK+k3kp6T9LykFyQ938K5pwKrS9tr0r6yvYG9Jf1c0p2S5g0Rw0JJPZJ6ent7W7i0mdnwPI+g0Eofwd8DR0fEhIjYKSJ2jIidRun644BZwFxgAXCxpImNhSJiUUTMiYg5kydPHqVLm1nuysNHXSNo7omIeGAE514LTC9tT0v7ytYASyLilYh4CPgPisRgZla58m9/xnmgpT6CHklXAzcAL9d3RsT1w7zubmCWpJkUCeA44IMNZW6gqAlcImkSRVPRqhZjNzN7TWqeRwC0lgh2An4HHFHaF0DTRBARGySdAtwMdAPfiYjlks4BeiJiSTp2hKQVQB/whYh4agTvw8xsk5X7Bfoy7iRoZfjoR0d68ohYCixt2Hdm6XkAn0sPM7O22nj4aAcD6TDfoczMsuUlJgq+Q5mZZctLTBR8hzIzy9bGw0c7GEiH+Q5lZpYtLzFRGOkdyo6vNCozszbYeImJzsXRaU0TQVo47hMR4TuUmdmYU68RdMnDR4cUEX2SDknPX2xPSGZm7VH/7R/X1ZV1Z3ErTUP3SloCXAP0J4MWZhabmW3W6j/+3V3Keh5BK4lgW+Ap4B2lfcPOLDYz29zVO4jHdck1gmZey8xiM7PNWf2uZN3dTgRNSbqENHS0LCL+spKIzMzapNZfI+hy09Awflh6vi3wHuDRasIxM2ufeiLYyjWC5iLiuvK2pCuBn1UWkZlZm9R/+7u7lPU8glZmFjeaBfgm82a2xauVOos9j6AJSS+wcR/B48BfVxaRmVmb1Eo1gvUbap0NpoNaaRrasR2BmJm1W7mzuBZ9HY6mc1q5H8F7JE0obU+U9O5qwzIzq16tNjChLOfO4lb6CM6KiOfqGxHxLHBWdSGZmbVHvWmoGDXU2Vg6qZVEMFiZVoadmplt1vqbhrq7sl6GupVE0CPpPElvSI/zgGVVB2ZmVrUINw1Ba4ngU8B64GrgKuAlinsUmJlt0QZWH/Xw0abS8tOntyEWM7O28uqjhVZGDd0iaWJpe2dJN1cblplZ9co1AjcNNTcpjRQCICKewTOLzWwMGBg+2uVRQ8OoSdqjviFpTwZZjdTMbEvjRecKrQwD/TLwM0k/BQQcCiysNCozszbobxrq9jLUTUXETZIOAA5Ku06NiHXVhmVmVr2a71AGDJMIJG0NHA+8Oe1aDrxQdVBmZu1QnkfQl3EiGLKPQNJsYAUwF3gkPeYCy9MxM7MtWnnUUATZzi5uViP4JvBXEXFLeaekw4ELgT+tMjAzs6qV5xFAcaMaqZMRdUazUUNTG5MAQET8CNitupDMzNqjPnx0XEoEufYTNEsEXZK2adwpaVtaXHRO0jxJD0paKWnI2cmS3icpJM1p5bxmZqOhPGqovJ2bZongu8B1ad4AAJJmAIuBy4c7saRuiiak+cBsYMFgfQuSdgQ+A9y1KYGbmb1W5VFD5e3cDJkIIuKrwE3A7ZLWSVoH/BS4JSLOaeHcBwIrI2JVRKynWLDumEHK/S1wLsVidmZmbTNQIxjoI8hR05nFEXFBROwBzARmRsSeEfHNFs89FVhd2l6T9vVL8xOmR8S/NjuRpIWSeiT19Pb2tnh5M7PmBoaPFj+FuQ4hbWWJCSLihYgY1fkDkrqA84DPt3D9RRExJyLmTJ48eTTDMLOMuWmo0FIiGKG1wPTS9rS0r25H4A+Bn0j6LcXM5SXuMDazdqk3DfUPH611MJgOqjIR3A3MkjQzzVA+DlhSPxgRz0XEpIiYEREzgDuBoyOip8KYzMz69Xn4KNB8ZvHXSs/fuaknjogNwCnAzcADwOKIWC7pHElHjyRYM7PRFKV7FkO+iaDZfIB5wJfS83OBV00uG05ELAWWNuw7c4iyczf1/GZmr0V5iYnydm6qbBoyM9us9XcW9w8fzTMTNKsRvF7S5yjuQVB/3i8izqs0MjOzivV3FivvGkGzRHAxxciexudmZmNCRNAl6EqJINd5BEMmgoj4SjsDMTNrt1oEXVL/iqO1TKsEQyYCSd9o9sKI+PToh2Nm1j59NejqUn+NINMKQdOmoWWl518Bzqo4FjOztupvGkrDZjx8tEFEXFZ/LunU8raZ2VhQbxrqkieUtSLPT8fMxrRa0JAIOhxQh3gegZllqxaBSqOGPI+ggaQXGKgJjJf0fP0QEBGxU9XBmZlVKfprBMW2h482iAjPGzCzMa2WOotVbxry6qNmZnnpqwXdXQM1AncWm5llphZFbSD3eQROBGaWLc8jKDgRmFm2PI+g4ERgZtnyPIKCE4GZZatxHoFrBGZmmWmcR5Dr6qNOBGaWrfrwUblpyMwsTwNNQ8V2rktMOBGYWbb6m4Z883ozszzV+m9VObCdIycCM8uW5xEUnAjMLFuNS0w4EZiZZaZ/iQmvPmpmlqeB4aPFtmsEZmaZeXXTUIcD6hAnAjPLVq1h9VHPIzAzy0x40TnAicDMMlZr7Cx2jWD0SZon6UFJKyWdPsjxz0laIek+SbdK2rPKeMzMyoolJnyrysoSgaRu4EJgPjAbWCBpdkOxe4E5EfEW4Frg76uKx8ysUXE/AtcIqqwRHAisjIhVEbEeuAo4plwgIm6LiN+lzTuBaRXGY2a2kVr/zes9j6AqU4HVpe01ad9QTgJuHOyApIWSeiT19Pb2jmKIZpaz+hITnkewGZB0AjAH+IfBjkfEooiYExFzJk+e3N7gzGzM6p9HkDoJMs0DjKvw3GuB6aXtaWnfRiQdDnwZ+C8R8XKF8ZiZbSS8+ihQbY3gbmCWpJmStgaOA5aUC0jaH/hfwNER8WSFsZiZvUr95vXdnkdQjYjYAJwC3Aw8ACyOiOWSzpF0dCr2D8AOwDWSfiFpyRCnMzMbdfV5BMp81FCVTUNExFJgacO+M0vPD6/y+mZmzdQab16faSLYLDqLzcw6oVZruDFNpm1DTgRmlq1aBF1deK2hTgdgZtYp9SUm1DWwnSMnAjPLVuPqo5nmAScCM8tXfdRQd+ajhpwIzCxb9VFDA0tMdDaeTnEiMLNs1dca8uqjZmaZKoaPlpaYyLRK4ERgZtmq+VaVgBOBmWWsPo/Ay1CbmWWqvgy1UodxOBGYmeWlvgw1FE1EbhoyM8tMfdQQFHMJ3DRkZpaZemcxFP0EfU4EZmZ5qa8+CkVCyDQPOBGYWb5qG/UReB6BmVl2akH/jevdWWxmlqFiGeriueR5BGZm2YlSZ3FXlzyPwMwsN+U+gm43DZmZ5acW0X8vAnkegZlZXiKif4kJSKOGnAjMzPJR/80vzyOo1ToYUAc5EZhZlup//W80j8A1AjOzfNQ7huvzCOTOYjOzvNT/+q/PI+jq8jLUZmZZaewj8OqjZmaZqf/od8tLTDgRmFmW+hqahrwMtZlZZiINFd14GWonAjOzbLx6+KjnEVRC0jxJD0paKen0QY5vI+nqdPwuSTOqjMfMrK4/EXQN3KHMncWjTFI3cCEwH5gNLJA0u6HYScAzEbEXcD5wblXxmJmV1TuG5c5ixlV47gOBlRGxCkDSVcAxwIpSmWOAs9Pza4ELJCkqaKhbfPdqLr591Wif1sy2UBtqGzcNdXeJO/5zHe8876cdjKq5Tx82i6P2nTLq560yEUwFVpe21wB/PFSZiNgg6TngdcC6ciFJC4GFAHvssceIgpk4fitm7brDiF5rZmPTW6ZN4O2zJgPwkbfN4NZfP9HhiJqbsN1WlZy3ykQwaiJiEbAIYM6cOSOqLRzx5t044s27jWpcZjZ2vO+t03jfW6d1OoyOqLKzeC0wvbQ9Le0btIykccAE4KkKYzIzswZVJoK7gVmSZkraGjgOWNJQZglwYnr+fuDHVfQPmJnZ0CprGkpt/qcANwPdwHciYrmkc4CeiFgC/DNwuaSVwNMUycLMzNqo0j6CiFgKLG3Yd2bp+UvAB6qMwczMmvPMYjOzzDkRmJllzonAzCxzTgRmZpnTljZaU1Iv8PAIXz6JhlnLmfLnMMCfRcGfQ2Esfw57RsTkwQ5scYngtZDUExFzOh1Hp/lzGODPouDPoZDr5+CmITOzzDkRmJllLrdEsKjTAWwm/DkM8GdR8OdQyPJzyKqPwMzMXi23GoGZmTVwIjAzy1w2iUDSPEkPSlop6fROx9MukqZLuk3SCknLJX0m7d9F0i2SfpP+u3OnY20HSd2S7pX0w7Q9U9Jd6XtxdVoyfUyTNFHStZJ+LekBSX+S4/dB0mfTv4n7JV0padscvw+QSSKQ1A1cCMwHZgMLJM3ubFRtswH4fETMBg4CPpne++nArRExC7g1befgM8ADpe1zgfMjYi/gGeCkjkTVXl8HboqINwL7UnweWX0fJE0FPg3MiYg/pFgq/zjy/D7kkQiAA4GVEbEqItYDVwHHdDimtoiIxyLinvT8BYp/9FMp3v9lqdhlwLs7E2H7SJoG/Dnw7bQt4B3AtanImP8cJE0A3k5xLxAiYn1EPEuG3weKZfi3S3dHHA88Rmbfh7pcEsFUYHVpe03alxVJM4D9gbuAXSPisXTocWDXDoXVTv8IfBGope3XAc9GxIa0ncP3YibQC1ySmsi+LWl7Mvs+RMRa4H8Aj1AkgOeAZeT3fQDySQTZk7QDcB1wakQ8Xz6Wbg86pscRS3oX8GRELOt0LB02DjgA+KeI2B94kYZmoEy+DztT1IJmAlOA7YF5HQ2qg3JJBGuB6aXtaWlfFiRtRZEEroiI69PuJyTtno7vDjzZqfja5GDgaEm/pWgafAdFW/nE1DQAeXwv1gBrIuKutH0tRWLI7ftwOPBQRPRGxCvA9RTfkdy+D0A+ieBuYFYaEbA1RafQkg7H1BapHfyfgQci4rzSoSXAien5icD/bnds7RQRfxMR0yJiBsX//x9HxPHAbcD7U7EcPofHgdWS9km7DgNWkNn3gaJJ6CBJ49O/kfrnkNX3oS6bmcWSjqRoI+4GvhMR/63DIbWFpEOA24FfMdA2/iWKfoLFwB4Uy3ofGxFPdyTINpM0FzgtIt4l6Q8oagi7APcCJ0TEy52Mr2qS9qPoMN8aWAV8lOKPwqy+D5K+AvwFxci6e4GPUfQJZPV9gIwSgZmZDS6XpiEzMxuCE4GZWeacCMzMMudEYGaWOScCM7PMORHYmCapT9Iv0gqT10ga36TsXElvK22fLOnDI7zuDEkfHKbMqZJeSuv/mHWME4GNdb+PiP3SCpPrgZOblJ0L9CeCiLgoIr47wuvOAJomAmABxWTH947wGmajwonAcnI7sJeko9Ka8/dK+pGkXdOCfCcDn001iEMlnS3pNABJb5B0k6Rlkm6X9Ma0/1JJ35D0b5JWSarPSv074NB0rs82BiLpDcAOwBkUCaG+f7ykxen+ET9Icc5Jx46QdIeke1LtZofqPirLiROBZSGtHzOfYob1z4CD0qJrVwFfjIjfAhdRrEW/X0Tc3nCKRcCnIuKtwGnAt0rHdgcOAd5FkQCgWMjt9nSu8wcJ6bh07duBfSTVV/v8BPBMun/EfwXemuKfRJE0Do+IA4Ae4HMj+jDMGowbvojZFm07Sb9Iz2+nWHdpH+DqtLja1sBDzU6Q/vJ+G3BNsSwNANuUitwQETVgRekHfTgLgPdERE3SdcAHgAsoEsrXASLifkn3pfIHUdxU6ecphq2BO1q8lllTTgQ21v0+IvYr75D0TeC8iFiS1h06e5hzdFGsU7/fEMfLa9FoiDLl6/8RMAu4pfSj/hBFIhjyZcAtEbGgSRmzEXHTkOVoAgPLC59Y2v8CsGNj4XT/hockfQCKFV0l7TvMNQY9V7IAODsiZqTHFGCKpD2BnwPHpuvMBv4oveZO4GBJe6Vj20vae5gYzFriRGA5OpuimWcZsK60/1+A99Q7ixteczxwkqRfAssZ/lan9wF9kn45SGfxccAPGvb9IO3/FjBZ0grgq+laz0VEL/AR4MrUXHQH8MZh36lZC7z6qNlmRFI3sFVEvJRGFv0I2Cfda9usEu4jMNu8jAduS3eVE/AJJwGrmmsEZmaZcx+BmVnmnAjMzDLnRGBmljknAjOzzDkRmJll7v8DprpiSWeqKagAAAAASUVORK5CYII=\n"
     },
     "metadata": {
      "needs_background": "light"
     }
    },
    {
     "output_type": "display_data",
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ],
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd5wV1fnH8c+zhd57Wap0ERHXGkWNDWwYK1hijEo0McWSqPn5M8YUNTFG87NHYyNqrBELUaOiiKCAKEqT3qQswtJZ2Huf3x8zu1wuy+5l4e7d3fm+X6993Tsz5848dxjuM3POzDnm7oiISHRlZToAERHJLCUCEZGIUyIQEYk4JQIRkYhTIhARiTglAhGRiFMiEKkhzGyjmXXPdBxS+ygRyD5nZheY2eTwh2u5mY0xs6MyHReAmf3AzD6qoMxYM7s8ad6xZrY0qczW8DuW/B2R9JmjE5ZtMjNPKt95T2J390buPn9PPpPKdxNRIpB9ysyuBe4B/gi0BToDDwDDKrGunFTmZdDV4Y9zyd+ExIXuPq5kGbB/OLtZQvnFJWWr2fcqkwX0m1EL6R9V9hkzawrcBvzE3V92903uvt3dX3P3X4ZlnjCz3yd8JvlMe6GZ3WBm04BNZtYjPJO+zMwWA++F5X5oZjPNbK2ZvWVmXRLW4WZ2pZnNMbNCM7s//BHrCzwEHBGekRdWzZ7ZlZndamYvmtkoM1sP/MDMDjWzCWHMy83sPjOrk/S9eoTv65rZXWa22MxWmtlDZlY/oewwM/vczNab2TwzG2JmfwCOBu4Lv/99YdkjzWySma0LX49MWM9YM/uDmY0HNgPXmdmUpO9yrZm9mtYdJmmlRCD70hFAPeCVvVzPCOBUoBlQHM47BugLnGxmw4BfA2cBrYFxwLNJ6zgNOAQYAJwHnOzuM4ErgQnhGXmzvYxzbw0DXiT4nv8EYsA1QCuCfXk88OPdfPYOoBcwEOgBdARuATCzQ4GngF+G6x4MLHT3/yHYVyVXMlebWQvgDeBvQEvgbuANM2uZsK2LgZFA47BctzCpJi5/qvK7QTJNiUD2pZbAancvrrBk+f7m7kvcfUvCvFvDK4wtBD/mt7v7zHBbfwQGJl4VAHe4e2FY/fI+wQ/mHsUQnpkXhlcOr1dQ5rM9XD8ECenf7h539y3uPsXdJ7p7sbsvBB4mSIA7MTMj+GG+xt3XuPsGgn0wPCxyGfAPd38nXPcyd5+1mxhOBea4+9Phdp8FZgGnJ5R5wt2nh8uLgH8BF4Wx7A90pez9IzWEEoHsS98CrfZBffeSCuZ1Ae5N+JFeAxjBWXGJFQnvNwON9jCGn7l7s5I/giuM8soM2sP1Q9L3NLNeZva6ma0Iq4v+SHB1kKw10ACYkrAP/hPOB+gEzEsxhg7AoqR5i9h5Xyb/ezwJXBAmpIuB58MEITWUEoHsSxOAIuDMcspsIvgRK9GujDJldYmbOG8J8KPEH2p3r+/uH6cQY3Xqbjc5lgcJzsZ7unsTguovK+Nzq4EtwP4J379p2CgNwf7ZL8VtfkOQWBN1Bpbt7jPuPhHYRtDecAHw9G62JTWEEoHsM+6+jqCe+n4zO9PMGphZrpkNNbM/hcU+B04xsxZm1g74RSU29RBwU1gtgZk1NbNzU/zsSiAvsRG2GmkMrAc2mlkf4KqyCrl7HPg78FczawNgZh3N7OSwyGPApWZ2vJllhcv6hMtWAonPIrwJ9LLglt8cMzsf6EfFVT1PAfcB29293NtxpfpTIpB9yt3/AlwL3AwUEJydXg38OyzyNPAFsBB4m6C+eU+38QpwJ/BcWIXyFTA0xY+/B0wHVpjZ6j3ddppdT3CGvYHgh768fXMDMBeYGO6D/wK9Adz9U+BS4K/AOuADdpz13wucE95t9Td3/5ag2us6gqq9XwGnuXtF++ZpoD8wak+/pFQ/poFpRKq/8P79GNAl8fmDDMZTH1gFDHL3OZmOR/aOrghEaob+wFZ2bgTPpKuASUoCtUO1f5pRJOrM7GzgEeAGd99WDeJZSNCIXd5NAVKDqGpIRCTiVDUkIhJxNa5qqFWrVt61a9dMhyEiUqNMmTJltbu3LmtZjUsEXbt2ZfLkyZkOQ0SkRjGz5CfIS6lqSEQk4pQIREQiTolARCTilAhERCJOiUBEJOLSlgjM7B9mtsrMvtrNcjOzv5nZXDObZmaV6c9dRET2UjqvCJ4AhpSzfCjQM/wbSdAXu4iIVLG0JQJ3/5Bg5KjdGQY85YGJQDMza5+ueERE1mzaxqiJi9hYtLejqdYumWwj6MjOQ+AtZefh8UqZ2Ugzm2xmkwsKCqokOBGpfe4cM4ub//0V371rLG9P37Uj1y+WFDLkng956IN5xOLR6YetRjQWu/sj7p7v7vmtW5f5hLSISLlWrt/KK1OXcULfNrRuXJcrR03hxSlLWb5uCx/NWc2/Ji3mwkc/YcmazdwxZhbnPTyBt6eviERCyGQXE8sIBtkukcfO46SKiOwz/xi/gOJ4nFtO25/WjetyxVOTuf6FL3Yq06NNI0Zddhjj567mrrdnM/LpKQzZvx0PXjQIs7KGj64dMpkIRgNXm9lzwGHAOndfnsF4RKSWemv6CkZNWMQpB7Snc8sGADx6ST5PfLyQ+rnZ9GrbmOYNc+nWqiF1c7I5++A8hg3swP+9N5d7353Da9OWc8aBHfZ4u58uWMNvRk/nrnMHsH+Hpvv6a+0zaUsEZvYscCzQysyWAr8BcgHc/SGCQbNPIRh3dTPBGKsiEnEFG4pYVriFgZ2a7fW64nHn5le/4plPFtOvfRNuGNKndFm93GyuPGa/3X42JzuLnx3fk7FfF3Dr6Ol0bdmAAXnNWL91O3Vzsqibk13utgs3b+Pnz01l+bqt/PSZqbz206NoWLd69vNZ4wamyc/Pd/U+KlI7xeLOWQ9+zNyVG5h268lkZ1W+Osbd+c3o6Tw1YRE/OqY715/Um9zsPW8Wnb1iA+c8+DEbiopp3iCXtZu3k5tt9GnXhAPymjK4ZytO3r/dLlVHP3nmM976agU3Du3DH96cyfn5nbjj7AGV/j57y8ymuHt+WcuqZ3oSkVpr6/YYM5avZ+7KjZwxsAP1cnecWT81YSFfLCkEYMHqjfRo07jS23ni44U8NWERIwd358YhfSpdx9+7XWPG3/RdXp6ylBnL19OtVSPWbdnOtKWFvPb5NzzzyWIG92rNXecMoE2TegCMnb2KN6Yt5/qTenH50d1ZunYLoyYu4ldD+tCiYZ1Kf6d0USIQkbSbsmgNT09YxEdzV7N6445hlxd8u6m0umZZ4Rb+/NZserRpxNxVG5n+zfpKJ4J5BRu5Y8wsju/ThpuGVj4JlGhSL5cffKfbLvNjcWfUxEXcPmYm//vqVzx8cT7biuPc9voMurVqyBWDuwNwXn4nnvh4IW9+uZyLDu+yV7GkQ424fVREaq4lazZz/sMTeXfWKgb3as01J/TiwQsHcdqA9jz20QKWrt2Mu3PzK1/iDo9+P586OVlM/2Z9mevbWFTMK1OX8vq0b5gw71tWrNu60/LiWJzrnv+CernZ3H7WAWm92yc7y7jkyK5cfHgX3p25im83FvHkxwuZX7CJW07rV9qO0Ld9Y3q0acToz7/B3VmzaVsFa65auiIQkbT6YmkhxXHn2SsOp3/HHXfOHNipGe/MWMmvX/mKw7q14P3ZBdx8al+6tmpIn3aNmf7NOiD4YX/4w/lMXVzIgXlNeW7SEpYVbtlpG+fl53H7WQPIzjIeGTefz5cU8rcRB5VW1aTbOQd34u/jFvDoRwt4esIijuvdmuP6tCldbmYMO7ADf3nnay57cjLvzVrF9Sf14ifH9agWt6UqEYhIWs1cvp6cLKNn20Y7ze/QrD7XntiL28fM4sOvCxiQ15RLw+qX/Ts0YcxXK1i3ZTuXPzmJSQvX0r5pPf47cyX7tW7IM5cfRuvGdVm5voh3Z63k8fEL2bo9zjG9WnPPO3M45YB2nD6g6nqs6d2uMQPymvLg2HnkZhv/e1q/XcqcMTBIBB9+XcAhXZtz19tfs3J9EbcN2z/jyUCJQETSatbyDezXulGZt1v+6Jj9OOfgPOas2kjPNo1K7xLq16Epz366hOue/4Ipi9Zyz/kDOfOgjqzasJXmDeqU3v3Ts21jjurZiuYN6nD3O18z+otvaNWoDr8b1r/Kf1zPOTiPaUvX8cPvdKN760a7LO/SsiF3n3cg+7VuxIC8pvzxzZn8fdwCerVtxMVHdK3SWJMpEYhIWs1asYH8rs13u7xlo7q0bFR3p3n7d2gCwH9nruT7R3ThzIOCbsjaNC67qudnx/dk+KGdWLZ2C11bNqR5Bu7MOffgThRtj3Ph4Z13W+asQXml728a2pe5qzZy2+szWFq4hQPzmjFk/3Zk7cUts5WlxmIRSZt1W7azrHALfdo12aPP9W3XhCyDpvVzueaEXil9pk3jehzUuXlGkgBA/TrZXDG4Ow3qpHZ+nZVl3HP+QRzUuTmPjVvAj//5GSOfnsKGrdvTHOmudEUgImkza3lw50+f9nt2G2j9OtmMHLwfAzs1zdgPe1Vo2iCX5390BNtjcUZNXMTv35jJz56dyuOXHlqlcSgRiMg+cffbsyncsp1fn9K39CGxWSs2AMEZ/p66cWifigvVErnZWVz6nW4sX7eVx8cvYMu2GPXrlN+Fxb6kqiER2WtfLVvH396by1MTFnH+wxNYvbEIgFkr1tO8QS5tm9StYA0CcOR+Ldkec6YsWlul21UiEJG9dud/ZtGsQS5/OfdAZq3YwI0vfcnmbcWMn/stfdo1yfjtkTVFftcWZGcZE+d/W6XbVdWQiOyVj+etZtyc1dx8al/OPjiP1RuLuH3MLM55cAJL1m7m92f2z3SINUajujkMyGvKhCpOBLoiEJG98s6MldTLzSrtQ+eyo7oxIK8pM5av59oTejG4l0YV3BOHd2/JF0sK2VSF4yorEYjIXpm8cC0HdWpe2kCck53FAxcO4ndn9ucnx/XIcHQ1zxHdW1Icr9p2AiUCEam0jUXFTP9mHYckPTCW17wBFx/eJSMPR9V0+V2bUycni/dmraqybSoRiEilfb64kLgHjZyybzSok8OxvVoz5qvlxONVM3CYEoGIVNqkhWvIMjio894PKyk7nDqgPSvXFzFlcdVUDykRiEhKyhrWdvKiNfRt34TG9XIzEFHtdXzfttTJyeKNacurZHtKBCJSoTkrN3D47e/yl7dnl84rKo4xdXEhh6haaJ9rVDeH43q35vVp3/D8pCXMWbkhrdvTcwQisotY3Jm6eC1fLltHl5YNuOXV6Xy7cRv/995cmjeoww+P6sYT4xeyeVuMk/q1zXS4tdKIQzszdnYBv3ppGgC92jbilyf34cQ07G8lAhHZyaoNWxn+8ETmr95UOq9+bjYvXnUkD42dx22vz2DVhiJGTVzE8X3acGSPVhmMtvY6tncbZtw2hMVrNjNuTgGvT1tOnZz0VOJYWfV+1Vl+fr5Pnjw502GI1EpFxTFGPDKRmcs38Ifv9efw7i2ZV7CRvOYN6NaqIVu3x7jp5S95ZeoycrKMt64ZzH5lDMIi1Y+ZTXH3/LKW6YpARErdOWY2ny0u5MELBzH0gGCoxw7N6pcur5ebzd3nHcigLs2pm5OlJFBLKBGICAAr129l1MRFDD+kU2kSKIuZcXHYnYTUDrprSEQAePiD+cTc1S1EBCkRiAir1m/lmU8XcebAjnRq0SDT4UgVUyIQqWW+3VjEx/NWl7msqDjGA2PncvxfxvLpgjUArN+6nR8+OQl3+Mlx+1VlqFJNqI1ApJa57oUvGDu7gMuP6sZNp/QlO+z4bd3m7Qz/+0RmLl9P47o5XPbkJG4c2odnP13MrOUb+Pv38+muxt9I0hWBSC3y2eK1jJ1dQJ92jXn0owX89rXpQHAl8KNRk5m7agOPXHwwY35xNA3qZPM/r3zF6g3buO+CQRzXp02Go5dM0RWBSC1yz3/n0KJhHV666kj+/NZsnvh4Iacf2IEnP17IxPlruOf8gZy0fzsAXrzySOYVbOSoHq3IydY5YZSl9V/fzIaY2Wwzm2tmN5axvLOZvW9mU81smpmdks54RGqzWSvW8+HXBYwc3J2GdXP45cm9ad+0Hhc++gmvT1vOr0/pw5kHdSwt36lFA47t3UZJQNKXCMwsG7gfGAr0A0aYWb+kYjcDz7v7QcBw4IF0xSNS2707MxjI5OxBeQA0rJvDbcP6UxyLc92JvRg5WA3BUrZ0Vg0dCsx19/kAZvYcMAyYkVDGgSbh+6bAN2mMR6RWe3/WKg7o2JTWjeuWzjuxX1s+/81JNFE30VKOdF4TdgSWJEwvDecluhW4yMyWAm8CPy1rRWY20swmm9nkgoKCdMQqUqMVbt7GZ4vXclzvXQeKVxKQimS6cnAE8IS75wGnAE+b2S4xufsj7p7v7vmtW+96oItE3YdzVhN3OFZ3/kglpDMRLAM6JUznhfMSXQY8D+DuE4B6gPq0FdlDY2etonmDXA7M05CRsufSmQgmAT3NrJuZ1SFoDB6dVGYxcDyAmfUlSASq+xHZA/G488HXBQzu1br04TGRPZG2RODuxcDVwFvATIK7g6ab2W1mdkZY7DrgCjP7AngW+IHXtAESRDLsy2Xr+HbTNo7rrWohqZy0PlDm7m8SNAInzrsl4f0M4DvpjEGktnt/9irMYHAvtZ9J5WS6sVhE9tL7swsY2KkZLRrWyXQoUkMpEYjUYN9uLGLa0kKO7aVqIak8JQKRGuyDrwtwh+P6qFpIKk+JQKQGe/mzZXRoWo/+HZpmOhSpwZQIRGqoeQUb+Wjuai44rDNZum1U9oISgUgN9fSEReRmG+cf0jnToUgNp0QgUgMVbt7GS1OWMrR/+506mROpDCUCkRpm6/YYlz85maLiOCMHd890OFILaIQykRrE3bn2+c+Zsngt940YRP+OaiSWvacrApEa5JlPF/Pmlyu4YUgfTh3QPtPhSC2hRCBSQ8xdtYHfvT6Do3u2YuTRqhKSfUeJQKQGKCqO8dNnP6dhnRz+ct6Bul1U9im1EYhUU18uXcdTExYSizvFcWfm8vU8dkk+bRrXy3RoUssoEYhUM/G487s3ZvD4+IU0rJNNbk4WhZu38/0junB837aZDk9qISUCkTRxd75Zt5WOzeqn/Jmi4hjXPf8Fr09bziVHdOG6k3tTPzeb6d+sp3+HJmmMVqJMiUAkDdyd3742gyc+XsjZg/K49Yx+NK6Xy7biONO/WUf9Otl8MLuAx8cvpE2TuhzTqzVH7teK+96fw/i533LT0D6MHNwds6AtYGAnDUEp6WM1bUCw/Px8nzx5cqbDECnXn/4ziwfGzuOwbi2YtHANjevlcnTPVkxauIaV64tKy32nR0uKtsf5bPFa4g7ZWcafzh7A2QfnZTB6qY3MbIq755e1rMIrAjOz5OEjzayuuxft7jMiUbbo2008MHYe5x6cx5/OGcDUJYWMmrCID+cU0KddE24+tR9m0LFZfQ7q3ByAdZu3M2H+ajo2a8ABeXpITKpWKlVDjwE/LJkws0bAq4SDzovIzp75dDHZWcb1J/fGzBjUuTmDwh/83WnaIJch/fWAmGRGKs8RLDWzBwDMrDnwNjAqrVGJ1FBFxTFemLyUE/q2oW0T3eYpNUOFiSAcbH6jmT1EkAT+4u6Ppz0ykRrolc+WsWbTNi48rEumQxFJ2W6rhszsrITJT4D/BT4F3MzOcveX0x2cSE1RuHkbP312KuPmrKZ328Yc1aNVpkMSSVl5bQSnJ01PBXLD+Q4oEYiE/jF+IR/NXc0NQ/pw8RFd1AWE1Ci7TQTufmlVBiJS3fznq+W0alSX/K4tyi23rTjOM58s5rjebbjq2P2qKDqRfafCNgIz+5OZNTGzXDN718wKzOyiqghOJJPuGDOLy5+azKoNW8stN+ar5azeWMQlR3atmsBE9rFU7ho6yd3XA6cBC4EewC/TGZRIdbA95hRu3s5NL33J7h68jMWdxz5aQLdWDTla7QJSQ6WSCEqqj04FXnD3dWmMR6TaiMWdxvVyeHfWKkZNXFRmmb++8zXTlq7j6uN6qF1AaqxUHih73cxmAVuAq8ysNVD+tbJILRBz55T+7SnYWMRtr88gr0UD6uVk07llAzo0rcdzk5Zw3/tzOT+/E2cN6pjpcEUqrcJE4O43mtmfgHXuHjOzTcCw9IcmklnxuJObY/z1vIGcft9HXPr4JACyDLq0bMiC1Zs4tFsLfjts/9LO4URqovKeI/iuu7+X+DxB0sGu20elViuOO9lmNG2Qy7MjD+f9WavIa16fyQvXMn7eav74vQM4/5BOZKtKSGq48q4IjgHeY9fnCSDF5wjMbAhwL5ANPOrud5RR5jzg1nCdX7j7BRWHLZJ+8biX1vt3bFafiw4PnhY+tncbrqd3JkMT2afKe47gN+FrpZ4nMLNs4H7gRGApMMnMRrv7jIQyPYGbgO+4+1oza1OZbYmkQ8ydHJ3tSwSUe9eQmWWbWauE6TpmNtLMZqaw7kOBue4+3923Ac+xa9vCFcD97r4WwN1X7Vn4IulTnHBFIFKb7TYRmNlwYA0wzcw+MLOTgPnAUODCFNbdEViSML00nJeoF9DLzMab2cSwKqmsWEaa2WQzm1xQUJDCpkX2XjxsIxCp7cprI7gZONjd55rZIGACcI67v7aPt98TOBbIAz40swPcvTCxkLs/AjwCwQhl+3D7IrulqiGJivKqhra5+1wAd/8MmLOHSWAZ0ClhOi+cl2gpMNrdt7v7AuBrgsQgklHxuOOOqoYkEsq7ImhjZtcmTDdLnHb3uytY9ySgp5l1I0gAw4HkO4L+DYwAHg/bInoRVD+JZFQs7FJCVUMSBeUlgr8DjcuZLpe7F5vZ1cBbBLeP/sPdp5vZbcBkdx8dLjvJzGYAMeCX7v7tnn4JkX0tFg8TQbYSgdR+5d0++tu9Xbm7vwm8mTTvloT3Dlwb/olUG6WJQFcEEgGpdDonEjmlVUNqI5AIUCIQKUM8rkQg0aFEIFKGYiUCiZDyOp0rt94+hbuGRGqskiuCLLURSASUd9dQyR1CvYFDgNHh9OnAp+kMSiTTStoI9ECZREGFdw2Z2YfAIHffEE7fCrxRJdGJZEhxLLwiUCKQCEiljaAtsC1hels4T6TWiuuBMomQVIaqfAr41MxeCafPBJ5MX0gimRdTY7FESCpDVf7BzMYAR4ezLnX3qekNSySzlAgkSlK9fbQBsN7d7wWWhv0HidRaeqBMoqTCRGBmvwFuIBhJDCAXGJXOoEQyLabbRyVCUrki+B5wBrAJwN2/YQ86nxOpiUoSgW4flShIJRFsCzuHcwAza5jekEQyT20EEiWpJILnzexhgvEIrgD+Czya3rBEMqvk9lE9RyBRkMpdQ3eZ2YnAeoKnjG9x93fSHplIBpU8UKaqIYmCChOBmd3p7jcA75QxT6RWKrlrSI3FEgWpVA2dWMa8ofs6EJHqJB4PXtVGIFFQXu+jVwE/Brqb2bSERY2B8ekOTCSTisNMoEQgUVBe1dAzwBjgduDGhPkb3H1NWqMSybC4HiiTCCmv99F1wDpgBICZtQHqAY3MrJG7L66aEEWqXqykakhtBBIBqTxZfLqZzQEWAB8ACwmuFERqrZiqhiRCUmks/j1wOPC1u3cDjgcmpjUqkQyLqbFYIiSVRLDd3b8Fsswsy93fB/LTHJdIRu3odC7DgYhUgVTGIyg0s0bAh8A/zWwVYb9DIrXVjqohZQKp/VI5yocBW4BrgP8A8wjGLRaptdRYLFGSShcTmwDMrAnwWtojEqkG4iXdUOuCQCIglS4mfgT8FtgKxAEj6Im0e3pDE8mc4tJuqJUJpPZLpY3geqC/u69OdzAi1UVpX0PKAxIBqRzm84DN6Q5EpDopqRpSG4FEQSpXBDcBH5vZJ0BRyUx3/1naohLJMFUNSZSkkggeBt4DviRoIxCp9dRYLFGSSiLIdfdrK7NyMxsC3AtkA4+6+x27KXc28CJwiLtPrsy2RPalmDqdkwhJ5XxnjJmNNLP2Ztai5K+iD5lZNnA/wdgF/YARZtavjHKNgZ8Dn+xh7CJpozGLJUpSuSIYEb7elDAvldtHDwXmuvt8ADN7juDhtBlJ5X4H3An8MoVYRKpETI3FEiGpPFDWrZLr7ggsSZheChyWWMDMBgGd3P0NM9ttIjCzkcBIgM6dO1cyHJHU6YpAoqS8Ecq+6+7vmdlZZS1395f3ZsNmlgXcDfygorLu/gjwCEB+fr7vzXZFUhGLO1kGpisCiYDyrgiOIbhbqKx+hRyoKBEsAzolTOeF80o0BvoDY8P/bO2A0WZ2hhqMJdNi7roakMgob4Sy34Rvb3P3BYnLzCyV6qJJQM+w7DJgOHBBwvrXAa0S1jkWuF5JQKqDeNzJ0tWAREQqdw29VMa8Fyv6kLsXA1cDbwEzgefdfbqZ3WZmZ+xZmCJVqzju5OiKQCKivDaCPsD+QNOkdoImBGMXV8jd3wTeTJp3y27KHpvKOkWqQizuZCkRSESU10bQGzgNaMbO7QQbgCvSGZRIpsXVRiARUl4bwavAq2Z2hLtPqMKYRDJOVUMSJam0EXzPzJqYWa6ZvWtmBWZ2UdojE8kgNRZLlKSSCE5y9/UE1UQLgR7oKWCp5WJxVQ1JdKSSCHLD11OBF8LbPkVqNSUCiZJU+hp6zcxmEQxgf5WZtSYYtlKk1tIDZRIlFV4RuPuNwJFAvrtvJxitbFi6AxPJpFjc1eGcRMZuE4GZ/Sph8nh3jwG4+yZAo5NJraaqIYmS8q4Ihie8vylp2ZA0xCJSbSgRSJSUlwhsN+/LmhapVeKu20clOspLBL6b92VNi9QqxXEnJ1uJQKKhvLuGDjSz9QRn//XD94TTKfU1JFJTxfRAmURIeV1MZFdlICLVifoakihJ5YEykcgpjikRSHQoEYiUIe56jkCiQ4lApAy6fVSiRIlApAxKBBIlSgQiZVBfQxIlSgQiZYjF0e2jEhlKBCJliMXjGqFMIkOJQKQMaiOQKFEiEClD3CFLiUAiQolApAzFqhqSCFEiEClDXI3FEiFKBCJlCNoIMh2FSNXQoS5ShuK4k52l/x4SDTrSRR+SjRAAAAtTSURBVMoQ9D6a6ShEqoYOdZEyaPB6iRIlApEyxFQ1JBGiI12kDGoslihJ66FuZkPMbLaZzTWzG8tYfq2ZzTCzaWb2rpl1SWc8IqmKueuBMomMtCUCM8sG7geGAv2AEWbWL6nYVCDf3QcALwJ/Slc8InsiFnc9UCaRkc4rgkOBue4+3923Ac8BwxILuPv77r45nJwI5KUxHpGUqbFYoiSdiaAjsCRhemk4b3cuA8aUtcDMRprZZDObXFBQsA9DFNlVPO6A+hqS6KgWzWFmdhGQD/y5rOXu/oi757t7fuvWras2OImc4jARqGpIoiInjeteBnRKmM4L5+3EzE4A/gc4xt2L0hiPSEririsCiZZ0XhFMAnqaWTczqwMMB0YnFjCzg4CHgTPcfVUaYxFJWSy8IlAbgURF2hKBuxcDVwNvATOB5919upndZmZnhMX+DDQCXjCzz81s9G5WJ1JlSqqGNDCNREU6q4Zw9zeBN5Pm3ZLw/oR0bl+kMuJKBBIx1aKxWKQ6ibkSgUSLEoFIkpiuCCRilAhEkqixWKJGiUAkSUwPlEnEKBGIJInpgTKJGCUCkSRqLJaoUSIQSVLa15DaCCQilAhEkuiBMokaJQKRJLp9VKJGiUAkSUmnc7p9VKJCiUAkiaqGJGqUCESSqK8hiRolApEkaiOQqFEiEEkS0+2jEjFKBCJJSh4oy8lWIpBoUCIQSaIrAokaJQKRJGojkKhRIhBJok7nJGqUCESSlDxQpqohiQolApEkeqBMokaJQCSJ2ggkapQIRJLENR6BRIwSgUiS4pg6nZNoUSIQSVJ6RaAHyiQilAhEksTiwauuCCQqlAhEksTiQSbI0v8OiQgd6iJJdjxQpv8eEg060kWShG3FqhqSyFAiEEmiqiGJGh3qIklKGotVNSRRoSNdJElpX0P63yERkdZD3cyGmNlsM5trZjeWsbyumf0rXP6JmXVNZzwiqdADZRI1OelasZllA/cDJwJLgUlmNtrdZyQUuwxY6+49zGw4cCdwfrpi2p1Y3NkeixOLOzF3YrHwNe4Ux514+BqLx4nFoTgexx3cg7PH4A88fC2ZhxPMZ8d8SpfvKO/uOMH88GSULDPMIOjlwMgyMCt5Dd7bTuWCaWdHXCXrpCSGeMnyYLvgYdkdMXoYQEksiesiab3hqoP4E79P6TZK1llSPnH/BIGW7pdqZNLCNYC6mJDoSFsiAA4F5rr7fAAzew4YBiQmgmHAreH7F4H7zMzc9/0vw78mLebhD+ezrTjOtuI4ReHrtjABSNVLTGDV7eR7/w5NsOoWlEiapDMRdASWJEwvBQ7bXRl3LzazdUBLYHViITMbCYwE6Ny5c6WCadGwLn3bN6FudhZ1coK/uuFrnexscrKNnCwjO/wL3meRncVOrzlZRlZ4Zp5lRlZWyZl6eLZOOG+Xs/gdZ+4l5Up/CBNeS356Es+6E8/OE8/cS87qE8+6DcjKKvlx3bHOHesPX8uMFRKvPhI/x07rTfwBL2MbCfsn+TtmJewPEake0pkI9hl3fwR4BCA/P79Sp+8n9mvLif3a7tO4RERqg3Q2Fi8DOiVM54XzyixjZjlAU+DbNMYkIiJJ0pkIJgE9zaybmdUBhgOjk8qMBi4J358DvJeO9gEREdm9tFUNhXX+VwNvAdnAP9x9upndBkx299HAY8DTZjYXWEOQLEREpAqltY3A3d8E3kyad0vC+63AuemMQUREyqdnJ0VEIk6JQEQk4pQIREQiTolARCTirKbdrWlmBcCiSn68FUlPLUeU9sMO2hcB7YdAbd4PXdy9dVkLalwi2BtmNtnd8zMdR6ZpP+ygfRHQfghEdT+oakhEJOKUCEREIi5qieCRTAdQTWg/7KB9EdB+CERyP0SqjUBERHYVtSsCERFJokQgIhJxkUkEZjbEzGab2VwzuzHT8VQVM+tkZu+b2Qwzm25mPw/ntzCzd8xsTvjaPNOxVgUzyzazqWb2ejjdzcw+CY+Lf4VdptdqZtbMzF40s1lmNtPMjoji8WBm14T/J74ys2fNrF4UjweISCIws2zgfmAo0A8YYWb9MhtVlSkGrnP3fsDhwE/C734j8K679wTeDaej4OfAzITpO4G/unsPYC1wWUaiqlr3Av9x9z7AgQT7I1LHg5l1BH4G5Lt7f4Ku8ocTzeMhGokAOBSY6+7z3X0b8BwwLMMxVQl3X+7un4XvNxD8p+9I8P2fDIs9CZyZmQirjpnlAacCj4bTBnwXeDEsUuv3g5k1BQYTjAWCu29z90IieDwQdMNfPxwdsQGwnIgdDyWikgg6AksSppeG8yLFzLoCBwGfAG3dfXm4aAUQhQGd7wF+BcTD6ZZAobsXh9NROC66AQXA42EV2aNm1pCIHQ/uvgy4C1hMkADWAVOI3vEARCcRRJ6ZNQJeAn7h7usTl4XDg9bq+4jN7DRglbtPyXQsGZYDDAIedPeDgE0kVQNF5HhoTnAV1A3oADQEhmQ0qAyKSiJYBnRKmM4L50WCmeUSJIF/uvvL4eyVZtY+XN4eWJWp+KrId4AzzGwhQdXgdwnqypuFVQMQjeNiKbDU3T8Jp18kSAxROx5OABa4e4G7bwdeJjhGonY8ANFJBJOAnuEdAXUIGoVGZzimKhHWgz8GzHT3uxMWjQYuCd9fArxa1bFVJXe/yd3z3L0rwb//e+5+IfA+cE5YLAr7YQWwxMx6h7OOB2YQseOBoErocDNrEP4fKdkPkToeSkTmyWIzO4Wgjjgb+Ie7/yHDIVUJMzsKGAd8yY668V8TtBM8D3Qm6Nb7PHdfk5Egq5iZHQtc7+6nmVl3giuEFsBU4CJ3L8pkfOlmZgMJGszrAPOBSwlOCiN1PJjZb4HzCe6smwpcTtAmEKnjASKUCEREpGxRqRoSEZHdUCIQEYk4JQIRkYhTIhARiTglAhGRiFMikFrNzGJm9nnYw+QLZtagnLLHmtmRCdNXmtn3K7ndrmZ2QQVlfmFmW8P+f0QyRolAarst7j4w7GFyG3BlOWWPBUoTgbs/5O5PVXK7XYFyEwEwguBhx7MquQ2RfUKJQKJkHNDDzE4P+5yfamb/NbO2YYd8VwLXhFcQR5vZrWZ2PYCZ7Wdm/zGzKWY2zsz6hPOfMLO/mdnHZjbfzEqeSr0DODpc1zXJgZjZfkAj4GaChFAyv4GZPR+OH/FKGGd+uOwkM5tgZp+FVzeN0rerJEqUCCQSwv5jhhI8Yf0RcHjY6dpzwK/cfSHwEEFf9APdfVzSKh4BfuruBwPXAw8kLGsPHAWcRpAAIOjIbVy4rr+WEdLwcNvjgN5mVtLb54+BteH4Ef8LHBzG34ogaZzg7oOAycC1ldoZIklyKi4iUqPVN7PPw/fjCPpd6g38K+xcrQ6woLwVhGfeRwIvBN3SAFA3oci/3T0OzEj4Qa/ICOB77h43s5eAc4H7CBLKvQDu/pWZTQvLH04wqNL4MIY6wIQUtyVSLiUCqe22uPvAxBlm9n/A3e4+Oux36NYK1pFF0E/9wN0sT+yLxnZTJnH7BwA9gXcSftQXECSC3X4MeMfdR5RTRqRSVDUkUdSUHd0LX5IwfwPQOLlwOH7DAjM7F4IeXc3swAq2Uea6QiOAW929a/jXAehgZl2A8cB54Xb6AQeEn5kIfMfMeoTLGppZrwpiEEmJEoFE0a0E1TxTgNUJ818DvlfSWJz0mQuBy8zsC2A6FQ91Og2ImdkXZTQWDwdeSZr3Sjj/AaC1mc0Afh9ua527FwA/AJ4Nq4smAH0q/KYiKVDvoyLViJllA7nuvjW8s+i/QO9wrG2RtFAbgUj10gB4PxxVzoAfKwlIuumKQEQk4tRGICIScUoEIiIRp0QgIhJxSgQiIhGnRCAiEnH/D36F/YFt72vdAAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ]
  },
  {
   "cell_type": "code",
   "source": [
    "# Threshold is the value at the prediction would assign a HF event. \n",
    "#Thus on the graph above when the line reaches above the threshold value it implies that there is high risk of HF\n",
    "\n",
    "def threshold(model,x,y):\n",
    "  #compute threshold \n",
    "  with tf.device('/cpu:0'):\n",
    "      y_pred_hf = model.predict(x)\n",
    "  y_pred_hf = y_pred_hf.squeeze()\n",
    "  hf_thresh = np.extract(y, y_pred_hf).mean()\n",
    "  return hf_thresh\n",
    "\n",
    "hf_thresh = threshold(attention,X_test,y_test)\n",
    "hf_thresh"
   ],
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "e54EcUq2mzaD",
    "outputId": "0a04f834-7a87-4494-f5b4-4137ac6a063d"
   },
   "execution_count": null,
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "0.14282021"
      ]
     },
     "metadata": {},
     "execution_count": 70
    }
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Addional Info: Classification SetUp "
   ],
   "metadata": {
    "id": "zZlXyxkJuigm"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "### For classification problems\n",
    "\"\"\"\n",
    "If you wanted to turn the seq2seq model into a classification probelm we have create \n",
    "minisequences and create a label which is the next HF label, i.e input = 10 times steps from 0 to i  and output = 1 HF value at i+1\n",
    "\"\"\"\n",
    "\n",
    "def trim(x):\n",
    "  z = []\n",
    "  window_size = 60\n",
    "  for i in x:\n",
    "    if(i.shape[0]>=window_size):\n",
    "      for j in range(0,i.shape[0]-window_size-1):\n",
    "        z.append(i[j:window_size+j])   \n",
    "  return np.array(z)\n",
    "\n",
    "train_cms, train_labels, test_cms, test_labels\n",
    "\n",
    "x1 = trim(train_cms)\n",
    "y1 = trim(train_labels)\n",
    "x2 = trim(test_cms)\n",
    "y2 = trim(test_labels)\n",
    "\n",
    "\n",
    "def transform_classify(y):\n",
    "  z = []\n",
    "  for i in y:\n",
    "    z.append(i[-1])   \n",
    "  return np.expand_dims(np.array(z), axis=1)\n",
    "\n",
    "y1 = transform_classify(y1)\n",
    "y2 = transform_classify(y2)\n",
    "print(y1.shape)\n",
    "print(y2.shape)\n",
    "#Padding essentially makes all inputs sequences the same length based on the value 'maxlen' and fills the sequence with 'value=-1' which helps the training ignore that part of the sequence \n",
    "def pad(train_cms, train_labels, test_cms, test_labels):\n",
    "    X_train = pad_sequences(train_cms, value=-1,dtype='float32')\n",
    "    y_train = np.expand_dims(pad_sequences(train_labels), axis=-1)\n",
    "\n",
    "    X_test = pad_sequences(test_cms, value=-1,maxlen=X_train.shape[1],dtype='float32')\n",
    "    y_test = np.expand_dims(pad_sequences(test_labels), axis=-1)\n",
    "    \n",
    "    return X_train, y_train, X_test, y_test\n",
    "    \n",
    "x1, y1, x2, y2 = pad(x1, y1, x2, y2)\n",
    "\n",
    "#shows the input dimensions\n",
    "print(\"X_train shape: {}\".format(x1.shape))\n",
    "print(\"X_test shape: {}\".format(x2.shape))\n",
    "print(\"y_train shape: {}\".format(y1.shape))\n",
    "print(\"y_test shape: {}\".format(y2.shape))"
   ],
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "eHN_WHGDuodf",
    "outputId": "fd2f9f55-4960-467c-e1a7-0e2054f0b808"
   },
   "execution_count": null,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "(181933, 1)\n",
      "(61076, 1)\n",
      "X_train shape: (181933, 60, 20)\n",
      "X_test shape: (61076, 60, 20)\n",
      "y_train shape: (181933, 1, 1)\n",
      "y_test shape: (61076, 1, 1)\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "TRY TO IMPLEMENT A MODEL TO CLASSIFY THESE SEQUENCES:\n",
    "\n",
    "Try to find the accuracy as well, any accuracy is good but try to get a value for it to see that it works\n",
    "\n",
    "This notebook was made by Harry Moroz.\n",
    "\n",
    "If you have questions or want me to look over an answer to the above problem  send email to: \n",
    "\n",
    "harry.moroz@mail.mcgill.ca\n"
   ],
   "metadata": {
    "id": "kjw0nGDT3ITw"
   }
  }
 ]
}