experiments: data extraction

mimic3_experiments/1_data_extraction/InclusionExclusion.ipynb

+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "import os, sys, time\n",
+    "from datetime import datetime, timedelta\n",
+    "import pickle\n",
+    "\n",
+    "from collections import Counter"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import yaml\n",
+    "config = yaml.safe_load(open('../config.yaml'))\n",
+    "data_path = config['data_path']\n",
+    "mimic3_path = config['mimic3_path']\n",
+    "\n",
+    "import pathlib\n",
+    "pathlib.Path(data_path, 'population').mkdir(parents=True, exist_ok=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "patients = pd.read_csv(mimic3_path + 'PATIENTS.csv', parse_dates=['DOB', 'DOD'], usecols=['SUBJECT_ID', 'DOB', 'DOD'])\n",
+    "admissions = pd.read_csv(mimic3_path + 'ADMISSIONS.csv', parse_dates=['DEATHTIME'], usecols=['SUBJECT_ID', 'HADM_ID', 'DEATHTIME', 'HOSPITAL_EXPIRE_FLAG'])\n",
+    "examples = pd.read_csv(data_path + 'prep/icustays_MV.csv', parse_dates=['INTIME', 'OUTTIME']).sort_values(by='ICUSTAY_ID') # Only Metavision\n",
+    "\n",
+    "examples = pd.merge(examples, patients, on='SUBJECT_ID', how='left')\n",
+    "examples = pd.merge(examples, admissions, on=['SUBJECT_ID', 'HADM_ID'], how='left')\n",
+    "examples['AGE'] = examples.apply(lambda x: (x['INTIME'] - x['DOB']).total_seconds(), axis=1) / 3600 / 24 / 365.25\n",
+    "\n",
+    "examples['LOS'] = examples['LOS'] * 24 # Convert to hours"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "tasks = ['ARF', 'Shock']\n",
+    "label_defs = { task: pd.read_csv(data_path + 'labels/{}.csv'.format(task)) for task in tasks }"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Source population 23620\n"
+     ]
+    }
+   ],
+   "source": [
+    "# Start\n",
+    "N = len(examples['ICUSTAY_ID'].unique())\n",
+    "print('Source population', N)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "assert (examples['INTIME'] <= examples['OUTTIME']).all()\n",
+    "assert (examples['DBSOURCE'] == 'metavision').all()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Exclude non-adults 23593\n"
+     ]
+    }
+   ],
+   "source": [
+    "# Remove non-adults\n",
+    "min_age = 18\n",
+    "max_age = np.inf # no max age\n",
+    "examples = examples[(examples.AGE >= min_age) & (examples.AGE <= max_age)]\n",
+    "print('Exclude non-adults', examples['ICUSTAY_ID'].nunique())\n",
+    "examples_ = examples"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "======\n",
+      "prediction time 4 hour\n",
+      "Exclude deaths 23499\n",
+      "Exclude discharges 23401\n",
+      "---\n",
+      "Outcome ARF\n",
+      "Exclude onset 15873\n",
+      "---\n",
+      "Outcome Shock\n",
+      "Exclude onset 19342\n",
+      "======\n",
+      "prediction time 12 hour\n",
+      "Exclude deaths 23319\n",
+      "Exclude discharges 23060\n",
+      "---\n",
+      "Outcome ARF\n",
+      "Exclude onset 14174\n",
+      "---\n",
+      "Outcome Shock\n",
+      "Exclude onset 17588\n"
+     ]
+    }
+   ],
+   "source": [
+    "for T in [4, 12]:\n",
+    "    print('======')\n",
+    "    print('prediction time', T, 'hour')\n",
+    "\n",
+    "    # Remove died before cutoff hour\n",

mimic3_experiments/1_data_extraction/LabelDistributions.ipynb

+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import pandas as pd\n",
+    "import numpy as np\n",
+    "import yaml\n",
+    "\n",
+    "data_path = yaml.full_load(open('../config.yaml'))['data_path']\n",
+    "\n",
+    "from matplotlib import pyplot as plt\n",
+    "import matplotlib"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "matplotlib.rcParams['figure.figsize'] = [8, 8]\n",
+    "matplotlib.rcParams['font.size'] = 15"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def visualize_labels(df, task):\n",
+    "    df['{}_ONSET_HOUR'.format(task)].plot.hist(bins=np.arange(-5, 75, 0.5), alpha=0.9)\n",
+    "    plt.xlim(-4,100)\n",
+    "    plt.xlabel('{} onset hour'.format(task))\n",
+    "    plt.savefig('Onset_{}-histogram.png'.format(task), dpi=300)\n",
+    "    plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "df_ARF   = pd.read_csv(data_path + 'labels/ARF.csv',   index_col='ICUSTAY_ID')\n",
+    "df_Shock = pd.read_csv(data_path + 'labels/Shock.csv', index_col='ICUSTAY_ID')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhsAAAHtCAYAAACwBLwkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xu4ZFV95//3h6u22NgKUVEUCSqijibpYPSXSBAvwIgYFK8ZJc5IzMSYnxgjoiSAtoI3YsCE4OigcRDBwQsqQRqDxHvaawQaUYOKEG31EIINQtrv/LH3geqiuvtUn1rn1u/X8+ynutZee9eqZdv1Ye21105VIUmS1Mp2890ASZK0tBk2JElSU4YNSZLUlGFDkiQ1ZdiQJElNGTYkSVJThg1JktSUYUOSJDVl2JAkSU3tMN8NWMh222232muvvea7GZIkzZkvf/nLP6mq3Sd5TsPGZuy1116sWbNmvpshSdKcSfK9SZ/TyyiSJKkpw4YkSWrKsCFJkpoybEiSpKYMG5IkqSnDhiRJasqwIUmSmjJsSJKkpgwbkiSpKcOGJElqyrAhSZKaMmxIkqSmDBuSJKkpw4YkSWrKsCFJkpqa87CRZJ8kf5fkG0k2JLl0C/VPTVJJ3jJi335JLkmyPsl1SU5Ksv1QnSQ5LskPktyc5LIkj57w15IkSZswHyMbDwcOBa4CvrW5ikn2A/47cOOIfSuA1UABhwMnAa8AThyqeixwPHAKcBhwE7A6yX1m9S0kSdKMzEfYuKCq9qyqI4HLt1D3NODtwNSIfS8B7gocUVUXV9UZdEHjmCTLAZLchS5svLGqTq+q1cCRdAHlpZP5OpIkaXPmPGxU1S9nUi/JM4F9gZM3UeUQ4KKqGhz1OIcugBzQv38csBw4d+Dzfw5c0B8vSZIa22G+GzBKkrsCbwWOraqfJxlVbV/gU4MFVfX9JOv7fRf0rxuAq4eOvRJ49qTbPcr+q1bf/ucvveaJc/GRkiQtKAv1bpRXA9cD79tMnRXADSPKp/p903VuqqoNI+osS7LTbBsqSZI2b8GNbCR5EPBnwIFVVfPw+UcDRwM84AEPmOuPlyRpyVmIIxsnAxcCVyW5R5J70LVz5/799DWVKWDXEcev4I4JpVPALsO3w/Z11lfVrcMHV9WZVbWyqlbuvvvuk/g+kiRt0xZi2HgocARdUJje9qS7e2QKuF9fby3dnIzbJdkTWNbvm66zPbDP0GfsO1BHkiQ1tBDDxv8ADhzafkR3R8mBwLq+3oXAU5LcfeDYZwM3A5/u33+Obo2OI6crJFlGt97Ghe2+giRJmjbnczb6H/tD+7f3A5b3t7kCfKKq1ow45hbgB1V16UDxGcDLgPOTnALsDZwAvG36dtiquiXJycDxSaboRjOOoQtZp036u0mSpDubjwmivwKcN1Q2/f5BwDUzOUlVTSU5CDid7jbXG4BT6QLHoJPpwsWrgXsBa4AnVdWPtqLtkiRpTHMeNqrqGmDkwhmbOWavTZRfATxhC8cWsKrfJEnSHFuIczYkSdISYtiQJElNGTYkSVJThg1JktSUYUOSJDVl2JAkSU0ZNiRJUlOGDUmS1JRhQ5IkNWXYkCRJTRk2JElSU4YNSZLUlGFDkiQ1ZdiQJElNGTYkSVJThg1JktSUYUOSJDVl2JAkSU0ZNiRJUlOGDUmS1JRhQ5IkNWXYkCRJTRk2JElSU4YNSZLUlGFDkiQ1ZdiQJElNGTYkSVJThg1JktSUYUOSJDVl2JAkSU0ZNiRJUlOGDUmS1JRhQ5IkNWXYkCRJTRk2JElSU4YNSZLUlGFDkiQ1ZdiQJElNGTYkSVJThg1JktSUYUOSJDVl2JAkSU0ZNiRJUlOGDUmS1JRhQ5IkNWXYkCRJTc152EiyT5K/S/KNJBuSXDq0/75J3pzk60luSvKDJO9JsseIc90vyYeS/EeSnyQ5PcmyEfVenOTqJLck+XKSgxp+RUmSNGA+RjYeDhwKXAV8a8T+3wB+D3g/cBjwSuAxwOeS7DJdKcmOwEXAA4HnAH8KHAmcOXiyJM8FzgDeCxwCXA58LMkjJvqtJEnSSDvMw2deUFUfAUjyQWC3of2fAfatqv+cLkjyFbpw8gzgPX3xM4GHAftU1b/29W4DzklyYlVd3dc7AXhPVb2ur/Np4NeAY4Hfn/zXkyRJg+Z8ZKOqfrmF/TcMBo2+7FvAemDwUsohwD9PB43eh4FbgYMBkuwNPAQ4d+jzz+uPlyRJjS2KCaJJ/guwjI0vu+wLrB2sV1W3At/p9zHwulE94Ergnkl2n3xrJUnSoAUfNpJsB7wduBr46MCuFcANIw6Z6vcx8Dpcb2po/5zaf9Vq9l+1ej4+WpKkOTcfczbG9UbgscABVXVb6w9LcjRwNMADHvCA1h8nSdKSt6BHNpL8T7q7UV5YVV8c2j0F7DrisBXcMXIx/Tpcb8XQ/ttV1ZlVtbKqVu6+u1dZJEmarQUbNpI8AzgN+POq+sCIKmu5Y07G9DE7AXtzxxyN6deN6vXvf1ZV6ybXYkmSNMqCDBtJfhf4P8BpVfWWTVS7EPjNJA8cKHsasDPwDwBV9V26SaVHDpx7u/79hZNvuSRJGjbnczb6FT4P7d/eD1ie5Jn9+0/QLdL1YbpRiQ8k+a2Bw9dV1Xf6P38QeA1wfpLj6S6VnAqcPbDGBnTrbLwvyTXAZ4EXAg8GnjfhryZJkkaYjwmiv0K3zsWg6fcPolstdFfgUcDnhuq9BzgKoKpuS3IwcDrdOhq/AM6hm+Nxu6p6f7/y6KuA4+lWEH1qVX1zQt9HkiRtxpyHjaq6BshmqpzVbzM517XA02dQ753AO2dyTkmSNFkLcs6GJElaOgwbkiSpKcOGJElqyrAhSZKaMmxIkqSmDBuSJKkpw4YkSWrKsCFJkpoybEiSpKYMG5IkqSnDhiRJasqwIUmSmjJsSJKkpgwbkiSpKcOGJElqyrAhSZKaMmxIkqSmDBuSJKkpw4YkSWrKsCFJkpoybEiSpKYMG5IkqSnDhiRJasqwIUmSmjJsSJKkpgwbkiSpKcOGJElqyrAhSZKaMmxIkqSmDBuSJKkpw4YkSWrKsCFJkpoybEiSpKYMG5IkqSnDhiRJasqwIUmSmjJsSJKkpgwbkiSpKcOGJElqyrAhSZKaMmxIkqSmDBuSJKkpw4YkSWrKsCFJkpoybMyj/VetZv9Vq+e7GZIkNWXYkCRJTc152EiyT5K/S/KNJBuSXDqiTpIcl+QHSW5OclmSR4+ot1+SS5KsT3JdkpOSbL8155IkSW3Mx8jGw4FDgauAb22izrHA8cApwGHATcDqJPeZrpBkBbAaKOBw4CTgFcCJ455LkiS1Mx9h44Kq2rOqjgQuH96Z5C50AeGNVXV6Va0GjqQLFS8dqPoS4K7AEVV1cVWdQRc0jkmyfMxzSZKkRuY8bFTVL7dQ5XHAcuDcgWN+DlwAHDJQ7xDgoqq6caDsHLoAcsCY55IkSY0sxAmi+wIbgKuHyq/s9w3WWztYoaq+D6wfqDfTc0mSpEYWYthYAdxUVRuGyqeAZUl2Gqh3w4jjp/p945zrdkmOTrImyZp169Zt9ZeQJEmdhRg25lVVnVlVK6tq5e677z7fzZEkadFbiGFjCthl+BZWulGK9VV160C9XUccv6LfN865JElSIwsxbKwFtgf2GSofnqOxlqF5F0n2BJYN1JvpuSRJUiMLMWx8DriR7hZVAJIso1sj48KBehcCT0ly94GyZwM3A58e81ySJKmRHeb6A/sf+0P7t/cDlid5Zv/+E1W1PsnJwPFJpuhGII6hC0anDZzqDOBlwPlJTgH2Bk4A3jZ9O2xV3TLDc0mSpEbmPGwAvwKcN1Q2/f5BwDXAyXSB4NXAvYA1wJOq6kfTB1TVVJKDgNPp1s24ATiVLnAM2uK5JElSO3MeNqrqGiBbqFPAqn7bXL0rgCdM4lySJKmNhThnQ5IkLSHzcRllm7X/qtXz3QRJkuacIxuSJKkpw4YkSWrKsCFJkpoybEiSpKYMG5IkqSnDhiRJasqwIUmSmjJsSJKkpgwbkiSpKcOGJElqyrAhSZKaMmxIkqSmDBuSJKkpw4YkSWrKsCFJkpoybEiSpKYMG5IkqSnDhiRJasqwIUmSmjJsSJKkpgwbkiSpKcOGJElqyrAhSZKaMmxIkqSmDBuSJKkpw4YkSWrKsCFJkpoybEiSpKYMG5IkqSnDhiRJasqwIUmSmjJsSJKkpsYKG0k+kOTJSdKqQZIkaWkZd2TjfsA/AN9P8vok+zRokyRJWkLGChtV9dvAQ4G/B14AXJXksiRHJblbiwZKkqTFbew5G1V1dVUdBzwQOBS4FngHcH2SdyX57Qm3UZIkLWJbPUG0qgr4NHAhcDmwC134uCzJl5M8ajJNlCRJi9lWhY0k/1+SdwL/BpwGfA14bFXdF3g0cCPw3om1UpIkLVo7jFM5yXHAC4F9gM8DLwc+UFXrp+tU1TeSvBa4bJINlSRJi9NYYQN4Gd2Ixbuq6qrN1FsLHL3VrZIkSUvGuGHj/lX1n1uqVFU/Bd61dU2SJElLybhzNn47yQtG7Ujy35IcMIE2SZKkJWTcsPEGYI9N7LtPv1+SJOl244aNRwBrNrHvK8DDZ9ccSZK01IwbNn4JrNjEvnttxfk2KclzknwlyU1JfpjkvUn2GKqTJMcl+UGSm/vVTB894lz7Jbkkyfok1yU5Kcn2k2qrJEnatHHDwWeBVyTZcbCwf/9y4DOTaFSSpwHvBz4HHA68Cng88PEkg20+FjgeOAU4DLgJWJ3kPgPnWgGsBqo/10nAK4ATJ9FWSZK0eePejXIcXaD4VpJzgOuB+wLPAe4J/M6E2vU84CtV9dLpgiQ3Ah+hezbLlUnuQhc23lhVp/d1Pg9cA7wUeG1/6EuAuwJHVNWNwMVJlgMnJHlTXyZJkhoZ90FsXwd+i27exouBv+pfvwQ8pqq+MaF27Qj8+1DZDf3r9OPtHwcsB84daN/PgQuAQwaOOwS4aChUnEMXQLx7RpKkxrbmQWyXV9WRVbVbVW3Xvz67qtZOsF3vBn4nyQuSLE/yEOD1wKeq6oq+zr7ABuDqoWOv7PcxUG+jtlXV94H1Q/UkSVIDE5vQOUlV9XHgKOBMuhGOq4DtgWcMVFsB3FRVG4YOnwKWJdlpoN4N3NkUIya7Jjk6yZoka9atWzer7yFJksafs0GSpwNHAPcH7jK8v6oeN9tGJTkQOAN4O91TZe8NnAB8KMkTRwSMiamqM+lCDitXrqxWnyNJ0rZi3AexHU93F8flwBXArS0aBbwV+GhVvWrgs79GdznkcOB8upGJXZJsPxQ+VgDrq2q6bVPAriM+Y0W/T5IkNTTuyMbRwJsHQ0Aj+9Ld+nq7qroqyc3Ar/ZFa+kurexDd5ll8NjBORprGZqbkWRPYNlQPUmS1MC4czbuDnyyRUOGfA/49cGCJA+ju4Pkmr7oc8CNwJEDdZbRrbdx4cChFwJPSXL3gbJnAzcDn550wyVJ0sbGHdk4F3gycEmDtgw6Azg1yXXcMWfjL+iCxicAquqWJCcDxyeZohulOIYuQJ02dK6XAecnOQXYm27+x9tcY0OSpPbGDRv/ALwlyT2Bixlxl0dVTWLk46/p5oP8Ed2iXDfQLSb26n4tjWkn04WLV9Mtl74GeFJV/WigPVNJDgJOp1uD4wbgVLrAIUmSGhs3bHywf/3v/Tas6OZRzEpVFfC3/baleqv6bXP1rgCeMNt2SZKk8Y0bNh7cpBVL0P6rVs93EyRJWhDGChtV9Z1WDZEkSUvT2CuIJtkxyYuT/F2STyTZpy9/ZpKHTr6JkiRpMRt3Ua996G593Q34Ct1TXpf3uw+ku+30hZNsoCRJWtzGHdn4a+DfgL2AJ3LHE1ihW7NiUo+YlyRJS8S4E0QPAJ5VVT9LMnzXyb8B951MsyRJ0lIx7sjGL4CdN7FvD0Y/XVWSJG3Dxg0bFwOvHlr6u5LsCLyUbtEvjWn/Vau9VVaStGSNexnllXTPJPk2cBHdIl6vAR4O3A141kRbJ0mSFr2xRjaq6vvAo4B30z1J9Xt0k0U/CvxGVV036QZKkqTFbdyRDarqp3TPIpEkSdqisRf1kiRJGse4i3pdTzdPY5Oqao9ZtUiSJC0p415GeRd3DhsrgIOAZcB7JtEoSZK0dIz7ILbXjipPsh1wHrB+Eo2SJElLx0TmbFTVL4F3Ai+bxPkkSdLSMckJog8Edprg+SRJ0hIw7gTRo0cU7wQ8DHgBcP4kGiVJkpaOcSeInjGi7D+BH9JdRvmLWbdIkiQtKeOGjR2HC6pqw4TaIkmSlqBx70YxWEiSpLGMO2fjeePUr6qzx2uOJElaasa9jPI+7ljUKwPlmyozbEiStI0b99bXx9A96fVE4L8A9+lfT+rLH0O3ougK4J6Ta6YkSVqsxh3ZOAX426p680DZj4FvJlkPvKmqDpxY6yRJ0qI37sjGbwFf38S+b9CNbEiSJN1u3LBxLXDUJvYdRbfehiRJ0u3GvYzyWuDsJPsBH6W7hPIrwNOARwLPnWzzJEnSYjfuOhvnJrkGOBb4A+DewI+Afwb+sKq+OPEWSpKkRW3ckQ2q6kvAEQ3aIkmSlqCteuprkl2TPDbJs5Lcoy+701LmkiRJY4WNJNsleQPdRNDPAu8H9u53fzTJX064fZIkaZEbd2RjFfDHwMuBh7DxiqEfppsoKkmSdLtx52y8EDi2qt6ZZPuhfd8BfnUyzZIkSUvFuCMbK4CrN7FvR2A4gEiSpG3cuGHjcuCwTex7CvDV2TVHkiQtNeNeRnkDcG6SnYHz6J7s+ogkhwF/BDx9wu2TJEmL3FgjG1V1PvAC4L8CF9NNED0L+EPgD6rqwkk3UJIkLW5bs6jX2UneDzwM2A34GXBFVf1y0o2TJEmL34zDRpK7AF8BXl5VFwFXNGuVJElaMmZ8GaWqbqEbyah2zZEkSUvNuHejvJ9uzoYkSdKMjDtn4zvAM5N8AfgE3RNfB0c6qqreOanGSZKkxW/csPFX/et9gf1H7C/AsCFJkm43btjwya6SJGksW5yzkeSTSR4KUFUbqmoDcABwl+n3g1vrBkuSpMVlJhNEnwjsOv2mfwDbxcBDWzWq/5wdkhyb5Ookv0hybZJTh+okyXFJfpDk5iSXJXn0iHPtl+SSJOuTXJfkpBEPkpMkSQ2MvahXL1uuMmtnAU8ATgTWAnsC+w3VORY4HnhlX+cYYHWSR1TVvwEkWQGsplsX5HC6J9O+lS5ovbb5t5AkaRu3tWGjqSQHA88GHlVVIxcP6xcZOxZ4Y1Wd3pd9HrgGeCl3BImXAHcFjqiqG4GLkywHTkjypr5MkiQ1MtN1NkYt5NVyca8XAZ/aVNDoPQ5YDpx7e4Oqfg5cABwyUO8Q4KKhUHEOXQA5YGItliRJI800bFyU5MdJfgxc35ddMl02uE2oXY8BvpXk9CQ39nMtzk+yx0CdfYENwNVDx17Z7xust3awQlV9H1g/VE+SJDUwk8soJzZvxZ3dBzgK+DrwHODuwJuADyX5raoqYAVw04g7YKaAZUl2qqpb+3o3jPiMqX7fRpIcDRwN8IAHPGAy30aSpG3YFsNGVc1H2Ei/HV5VPwVIcj3wabpJo5e0+uCqOhM4E2DlypU+B0aSpFka99koc2UK+JfpoNH7DHArd9yRMgXsMuIW1hXA+n5UY7rertzZin6fJElqaKGGjSsZfXttgF/2f14LbA/sM1RneI7GWobmZiTZE1g2VE+SJDWwUMPGx4BHJtltoOzxdMulf71//zngRuDI6QpJlgGHARcOHHch8JQkdx8oezZwM91lGUmS1NBCDRtnAj8FLkhyWJLnAX8PrK6qzwBU1S3AycBxSf44yUHAeXTf6bSBc50B/AI4P8kT+wmgJwBvW2hrbOy/ajX7r1o9382QJGmiFuSiXlV1Y5InAH9NtybGrcBHgJcPVT2ZLly8GrgXsAZ4UlX9aOBcU30QOZ1uDY4bgFPpAockSWpsQYYNgKr6NnDoFuoUsKrfNlfvCrq7WCRJ0hxbqJdRJEnSEmHYkCRJTRk2JElSU4YNSZLUlGFDkiQ1ZdiQJElNGTYkSVJThg1JktSUYUOSJDVl2JAkSU0ZNiRJUlOGDUmS1JRhQ5IkNWXYkCRJTRk2JElSU4YNSZLUlGFDkiQ1ZdiQJElNGTYkSVJThg1JktSUYUOSJDVl2JAkSU0ZNiRJUlOGDUmS1JRhQ5IkNWXYkCRJTRk2JElSU4YNSZLUlGFDkiQ1ZdiQJElNGTYkSVJTO8x3A3Rn+69affufv/SaJ85jSyRJmj1HNiRJUlOGDUmS1JRhQ5IkNWXYkCRJTRk2JElSU4YNSZLUlGFDkiQ1ZdiQJElNGTYkSVJThg1JktSUYUOSJDVl2JAkSU0ZNiRJUlOGDUmS1NSiCBtJ7pfkpiSVZJeB8iQ5LskPktyc5LIkjx5x/H5JLkmyPsl1SU5Ksv3cfgtJkrZNiyJsAG8GbhpRfixwPHAKcFhfZ3WS+0xXSLICWA0UcDhwEvAK4MTGbZYkSSyCsJHk8cDBwFuGyu9CFzbeWFWnV9Vq4Ei6UPHSgaovAe4KHFFVF1fVGXRB45gky+fiO0iStC1b0GGjv9RxGt1oxE+Gdj8OWA6cO11QVT8HLgAOGah3CHBRVd04UHYOXQA5oEGzJUnSgAUdNuhGJXYG3jFi377ABuDqofIr+32D9dYOVqiq7wPrh+pJkqQGFmzYSHIv4HXAMVV124gqK4CbqmrDUPkUsCzJTgP1bhhx/FS/b/hzj06yJsmadevWbf0XkCRJwAIOG8Aq4AtV9Ym5/NCqOrOqVlbVyt13330uP1qSpCVph/luwChJHg68CHh8knv0xcv6112TbKAbmdglyfZDoxsrgPVVdWv/fgrYdcTHrOj3SZKkhhZk2AAeDOwIfH7EvmuBdwFnA9sD+wBXDewfnqOxlqG5GUn2pAsvG83lkCRJk7dQw8ZngAOHyg4GXgUcCnwX+B5wI93trq8HSLKMbr2NMweOuxB4ZZK7V9V/9GXPBm4GPt3qC0iSpM6CDBtV9RPg0sGyJHv1f/ynqrqpLzsZOD7JFN0oxTF081BOGzj0DOBlwPlJTgH2Bk4A3jZ0O6wkSWpgQYaNMZxMFy5eDdwLWAM8qap+NF2hqqaSHAScTrcGxw3AqXSBQ5IkNbZowkZVnQWcNVRWdHetrNrCsVcAT2jVNkmStGkL+dZXSZK0BBg2JElSU4YNSZLUlGFDkiQ1ZdhY4PZftZr9V62e72ZIkrTVDBuSJKkpw4YkSWrKsCFJkpoybEiSpKYMG5IkqSnDhiRJasqwIUmSmjJsSJKkpgwbkiSpKcOGJElqyrAhSZKaMmxIkqSmDBuSJKkpw4YkSWrKsCFJkpoybEiSpKYMG5IkqSnDhiRJasqwIUmSmjJsSJKkpgwbkiSpKcOGJElqyrCxSOy/ajX7r1o9382QJGlshg1JktSUYUOSJDVl2JAkSU0ZNiRJUlOGDUmS1JRhQ5IkNWXYkCRJTRk2JElSU4YNSZLUlGFDkiQ1ZdiQJElNGTYkSVJThg1JktSUYUOSJDVl2FhkfNS8JGmxMWxIkqSmDBuSJKkpw4YkSWpqQYaNJEcm+WiSHya5KcmXkzx3RL0XJ7k6yS19nYNG1Llfkg8l+Y8kP0lyepJlc/NNJEnSggwbwDHATcDLgacB/wicneRPpiv04eMM4L3AIcDlwMeSPGKgzo7ARcADgecAfwocCZw5N19DkiTtMN8N2ITDquonA+8/lWQPuhByWl92AvCeqnodQJJPA78GHAv8fl/nmcDDgH2q6l/7ercB5yQ5saqubv5NJEnaxi3IkY2hoDHtq8AeAEn2Bh4CnDtwzC+B8+hGOaYdAvzzdNDofRi4FTh4ws2WJEkjLMiwsQmPBb7V/3nf/nXtUJ0rgXsm2X2g3kZ1qupW4DsD55AkSQ0tirDRT/x8OvDWvmhF/3rDUNWpof0rRtSZrrdiRDlJjk6yJsmadevWbX2jJUkSsAjCRpK9gLOBj1TVWa0/r6rOrKqVVbVy99133/IBkiRpsxZ02EhyT+BC4HvA8wd2TY9g7Dp0yIqh/VMj6kzXmxpRLkmSJmzBho1+LYyPATsBT62q9QO7p+dhDM+72Bf4WVWtG6i3UZ0kOwF7c+f5HpIkqYEFGTaS7EB3Z8mDgYOr6seD+6vqu3STRY8cOGa7/v2FA1UvBH4zyQMHyp4G7Az8Q5vWS5KkQQt1nY2/AQ6lW4TrXknuNbDvq1X1C7p1Nt6X5Brgs8AL6cLJ8wbqfhB4DXB+kuPpLqmcCpztGhuSJM2NhRo2nty/vn3EvgcB11TV+5PsArwKOJ5uBdGnVtU3pytW1W1JDgZOp1uT4xfAOcArWzZekiTdYUGGjaraa4b13gm8cwt1rqW7bVaSJM2DBTlnQ1u2/6rV7L9q9Xw3Q5KkLTJsSJKkpgwbkiSpKcPGIuflFEnSQmfYkCRJTRk2JElSU4YNSZLUlGFDkiQ1ZdiQJElNGTYkSVJThg1JktSUYWOJcL0NSdJCZdiQJElNGTYkSVJThg1JktSUYUOSJDVl2JAkSU0ZNiRJUlOGDUmS1JRhQ5IkNWXYkCRJTe0w3w3QZA2uIvql1zxxHlsiSVLHkQ1JktSUYUOSJDVl2JAkSU0ZNiRJUlOGDUmS1JRhYwnbf9Xqje5OkSRpPhg2JElSU4YNSZLUlGFjG+DlFEnSfDJsSJKkpgwbkiSpKcOGJElqyrAhSZKaMmxsg5wwKkmaS4YNSZLUlGFDkiQ1tcN8N0Bzx0snkqT5YNiYsMX0gz7c1i+95onz1BJJ0lLmZRRJktSUIxu63ahRmZmOdszm2HFMf46jMJK0eDiyobF426wkaVyObGhGhgPG1owwDB8zzjk29fmDHO2QpIXJsKHNmsQoxkyCwvC+4UAiSVq8DBtakGYTMjY1YrK50RDvzJGkdraJsJFkP+A04LHADcD/Ak6sqg3z2rAlYCajFLM5x2zM5DKNIyeS1N6SDxtJVgCrgSuAw4FfBd5KNzn2tfPYNM2RSY6SzGR0xFGRmbPPpG3Dkg8bwEuAuwJHVNWNwMVJlgMnJHlTX6YlYJKjFOPMM5nbqORFAAAMw0lEQVRpO7bmB3VzP8ZbuvSzNcfO5Htuqe5sgsPgOcf5Pls6VtL8SlXNdxuaSnIZcF1VPWeg7AHA94CnVdUFmzp25cqVtWbNmrE+z2F56c5mGlC2Zl7N5v4/t6XAMpO7miY5+rK5z5/N587mTq9x2riptmppSfLlqlo50XNuA2Hjx8DfVNUJQ+U/B06oqjdv6thRYWOc/7qStDBtzd1OM5lMPJsRo3GP2dznzjfDyOJm2NgKSW4DXllVfzVUfi3w3qo6bqj8aODo/u1Dgatm2YTdgJ/M8hy6M/u1Dfu1Dfu1Dfu1jYdW1d0necJtYc7GWKrqTODMSZ0vyZpJJ0TZr63Yr23Yr23Yr20kGW/+wAxsC8uVTwG7jihf0e+TJEkNbQthYy2w72BBkj2BZf0+SZLU0LYQNi4EnpJk8PrTs4GbgU/PwedP7JKMNmK/tmG/tmG/tmG/tjHxft0WJoiuoFvQ65vAKcDewNuAv6oqF/WSJKmxJR824Pblyk9n4+XKT3C5ckmS2tsmwoYkSZo/28KcjXmRZL8klyRZn+S6JCcl2X6+27VYJDkyyUeT/DDJTUm+nOS5I+q9OMnVSW7p6xw0H+1djJLcr+/bSrLLQHmSHJfkB0luTnJZkkfPZ1sXgyQ7JDm2//v4iyTXJjl1qI59O4Ykz0nylf7v6Q+TvDfJHkN17NPNSLJPkr9L8o0kG5JcOqLOjPpwNr9rho0GBh7+VnQPfzsJeAVw4ny2a5E5BrgJeDnwNOAfgbOT/Ml0hT58nAG8FzgEuBz4WJJHzH1zF6U30/XxsGOB4+nmOB3W11md5D5z2LbF6CzgZcBbgCfT9ePNQ3Xs2xlK8jTg/cDn6P4dfRXweODjSQZ/u+zTzXs4cCjdApXf2kSdLfbhrH/Xqsptwhvwaro1PJYPlP05sH6wzG2zfbjbiLKzgX8deH8V8O6B99sB/wK8b77bv9A3un+0fwb8Wf+Pxy59+V2Afwf+YqDu3YB1wOvnu90LdQMOBm4D9ttMHft2vD49B/jyUNnT+r+vD7NPZ9yP2w38+YPApUP7Z9SHs/1dc2SjjUOAi2rjJ8qeQ/f02QPmp0mLS1WNWoL4q8AeAEn2Bh4CnDtwzC+B8+j6X5vQD3ueRvdfJsP9/DhgORv368+BC7BfN+dFwKeq6orN1LFvx7Mj3Y/goBv61/Sv9ukW9P8ubs5M+3BWv2uGjTb2ZWjBsKr6Pl0C3HfkEZqJx3LHMOB0Pw4vzHYlcM8ku89ZqxaflwA7A+8YsW9fYANw9VD5lfh3d3MeA3wryelJbuyvaZ8/NL/Avh3Pu4HfSfKCJMuTPAR4PRuHOvt09mbah7P6XTNstLGCOxL4oKl+n8bUT/x8OvDWvmi6H4f7eWpovwYkuRfwOuCYqrptRJUVwE1159vCp4BlSXZq3cZF6j7AUcCjgecAfwD8BvChJNP/FW7fjqGqPk7Xp2fSjXBcBWwPPGOgmn06ezPtw1n9rvkgNi14Sfaim6/xkao6a14bs/itAr5QVZ+Y74YsMem3w6vqpwBJrqdbpfgJwCXz2LZFKcmBdBPA3063EvS9gRPoAtwTR/w4agEzbLThw98mJMk96f6h+R7w/IFd0/24Kxun7RVD+9VL8nC6uQWPT3KPvnhZ/7prkg10/bZLku2H/jFfAayvqlvnrsWLyhTw3emg0fsMcCuwH13YsG/H81bgo1X1qumCJF+jG8o/HDgf+3QSZtqHs/pd8zJKGz78bQKSLAM+BuwEPLWq1g/snu7H4WuF+wI/q6p1c9DExebBdJPuPk/3j8MUd8zbuJZu0uhauqHqfYaOvdP1Wm3kSu6YtDgowPQEPft2PPsCXxssqKqr6G4n/tW+yD6dvZn24ax+1wwbbcz3w98WvSQ70N1Z8mDg4Kr68eD+qvou3WTRIweO2a5/f+EcNnUx+Qxw4NB2Sr/vULp1Nz4H3MjG/bqM7t57+3XTPgY8MsluA2WPpwt3X+/f27fj+R7w64MFSR5Gd/fDNX2RfTp7M+3D2f2uzfc9wEtxoxtWuh64GHgicDTdIine9z3zPjyT7n76lwG/NbTt3Nd5Lt0s6tfS/XCe1f/Ff8R8t3+xbHQT8G5fZ6MvezXdDPM/Bg4CPk53i+y957u9C3Wju3Xw+3SjRocBzwN+AFw8VM++nXmf/indqNBb+39Hn083SfRfgbvZpzPux2XAM/vt83SLH06/XzbTPpzt79q8d8RS3eiu036q//G7nu4OgO3nu12LZaP7L5faxLbXQL0XA98GfgF8BThovtu+mLZNhI0Ar6G7tHIz8E/Ar813Wxf6RjcM/Qng53SXqM4CVgzVsW9n3p8B/gj4Rt+nPwQ+AOxtn47Vj3tt6d/SmfbhbH7XfBCbJElqyjkbkiSpKcOGJElqyrAhSZKaMmxIkqSmDBuSJKkpw4YkSWrKsCEtckn+NUklGV5umCS/2++b3v49yReTPH1E3WuG6k5vvz8332T2kjwryVEzrFtJXtq4SZIwbEiLWpLH0i3aA92KqpvyfOCxfZ2fAucnefyIemf39Qa3f5hUe+fAs+gWKpO0gPjUV2lxey7d6orf7P/8uk3U+0ZVfRMgyaV0S2n/PnDZUL3rq+oLbZqqQUm2p1t90SeTaslzZENapPofq2cBHwXeDTwsyaO2dFx1T8/9NrDnhNrxrCT/kuQXSX6QZFX/IL3p/Uf1lywemeTiJD9PsjbJEUPn+e0k/5Tkxn77WpIjh+r8jySX95/1vSR/PrDvLOAZwAEDl4BO2ELzt0/yhiTrkvw4yTuS7Dz0mY9OckmS9UmmkvyfJPce2D99qeoRQ8ddmuSDg+1LsibJ05NcDtwCPGYL7ZOWBMOGtHgdCNwbOAf4IHAbm7+UAtz+dNz70z3QasTu7DCwbb+Fcz2Z7nkVXwEOp3tM/Z8Bp4+ofjZdMPo94GrgnCT378+znO7Jqd+lCwzPBP4euMfAZ70S+Fvgw8BT+z+/bmDexeuAfwS+yh2XgP7XZjsDXgHsQTfK82bgD+keADb9mbsDl9I9zOp5wJ8ABwAXJ9lpC+ceZS/gTcAbgUMY/b+BtPTM90Ni3Nzctm4D3kX3wK+d+vcfo3uAXQbq/C7dA5ceRXfZdHe6H9UbgIcOne8a7vygpmu30IYvAP84VPbndE/jvX///qj+XC8aqHMv4D+Bl/TvV/Z17r6Jz1lO94TJvxwqPwn4N/qHQdGFrktn2H8FXDZU9mHgCwPvT+77avlA2WP6Y5871MePGDrXpcAHB96f1dd79Hz/3XFzm+vNkQ1pEer/q/oI4EN1xzX/c4AH0v0X/bCv0Y18/Bg4Bjiqqq4aUe99wG8ObIdupg3bA78OnDe06wN0o6bD7fjk9B+q6qd9W+7fF32HLkycneTwJPcYOvaxwN2A8wZHXuieQHnvgfOM65ND768YOtf+wCer6saBtn+RLpj99lZ83g+r6mtbcZy0qBk2pMXpELpLDJ9Ico/+x/lS4BeMvpTyHLrw8AzgKuB/J9ljRL0fVdWage0bm2nDbsCOwI+Gz9G/3nOo/Iah97cCdwGoqingSf35zgXWJfl4kr0HPgvgcrrQNL39Y1++tfNPNtmm3n258/ejLxv+fjMx6lzSkmfYkBan6UBxHt2llCm6O0x2Bo4cMdfi8j48nA8cRjcH4fhZtuEndD/4vzJUPj158mfjnKyqvlBVB9OFqCOAh9DN8xg811PZeORlevv6uI2foeu58/eD7jtOt+mW/nV4DseKEcfVhNolLSqGDWmRSXI3usDwfrpJooPbMXQ/hE/Y1PFV9R26iZNHJRn1QzojVbUB+DJw5NCuZwG/BD6/lee9uaouoLvDZr+++PPAzcAeQyMv09t/9PWGRyZm64vAU5LcfbogyW/STfT8TF90bf/6sIE6ewL7TrAd0qLmOhvS4nM43cjE2/v5A7dL8lngNXQjHxdv5hxvAl5Md3fFbEY4/hK4KMn/ppsz8ki6u0LeWVXXbvbIAUn+K/Aiugma3wfuR3dnyKcAquqG/jbWtyd5IN36INvRjX4cWFW/159qLXB4v0LqtcB1VXXdLL7f24A/6r/jKcAudJNG/wX4v33brk2yhu7OmPV9u45jzJEdaSlzZENafJ4LXD0cNACq6ja6OQ9HDK8XMVTve3STQf9nP1KyVarqk3TzQVYCFwD/P/BWYNxlwL9Nd4nhDXSTNt9Et3LpiwY+603A0XTzVT5CN7LzfOCfBs7zN/3x7wb+ua+/1apqHd2I0S39572j/7wn1caLcT2XLiS9r/8OJ9HNjZFEf4ucJElSK45sSJKkpgwbkiSpKcOGJElqyrAhSZKaMmxIkqSmDBuSJKkpw4YkSWrKsCFJkpr6f5wT5Qho37yvAAAAAElFTkSuQmCC\n",
+      "text/plain": [
+       "<Figure size 576x576 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhEAAAHtCAYAAACnJiztAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xu4ZFV95//3xwbUDkIOF+8oIioiJmRC2miMRMAL/kQcIkGT+Rnj/CTMSJgRxxEQEiB2BK9JJAZxzA81j0FMSBQjIo2KMZpoe02ERiKioqCgTQg0CDbf+WPvA0VR3X1q9ak65/R5v56nntO19qq9Vy2asz+99tprp6qQJEka1/0WugGSJGlpMkRIkqQmhghJktTEECFJkpoYIiRJUhNDhCRJamKIkCRJTQwRkiSpiSFCkiQ12W6hG7BQdtttt9pzzz0XuhmSJE3FF7/4xRuravf53OeyDRF77rkna9euXehmSJI0FUm+Pd/79HKGJElqYoiQJElNDBGSJKmJIUKSJDUxREiSpCaGCEmS1MQQIUmSmhgiJElSE0OEJElqYoiQJElNDBGSJKmJIUKSJDUxREiSpCaGCEmS1MQQIUmSmhgiJElSE0OEJElqYoiQJElNDBGSJKnJdgvdgG3RqtVr7v7z5193yAK2RJKkyXEkQpIkNTFESJKkJoYISZLUxBCxFVatXnOv+Q+SJC0nhghJktTEECFJkpoYIiRJUhNDhCRJamKIkCRJTQwRkiSpiSFCkiQ1MURIkqQmhghJktTEECFJkpoYIiRJUhNDhCRJamKIkCRJTQwRkiSpiSFCkiQ1MURIkqQmUw8RSbZLckKSq5L8JMm1Sd42VCdJTkry3SS3Jfl0kv1H7GvfJJcm2ZDk+0lOT7Jiet9GkqTla7sFOOa5wEHAacA6YA9g36E6JwCnAK/p6xwPrEmyX1VdD5BkBlgDXA4cDjwWeAtdMDp54t9CkqRlbqohIslzgaOAn6+qyzdR5wF0IeINVXVWX/Y54BrgWO4JCMcADwSOqKqbgUuS7AScmuSNfZkkSZqQaV/OeDnwiU0FiN7TgJ2A82cLqupW4ELg0IF6hwIXD4WF8+iCxYHz1mJJkjTStEPEU4BvJDkryc39XIYLkjx8oM4+wEbgqqHPXtFvG6y3brBCVX0H2DBUT5IkTcC0Q8RDgZcB+wMvBn4H+EXgb5OkrzMD3FJVG4c+ux5YmWSHgXo3jTjG+n7bfSQ5OsnaJGtvuOGGrfoikiQtd9OeWJn+dXhV/QggyXXAZXSTLS+d5MGr6hzgHIADDjigJnksSZK2ddMeiVgP/MtsgOh9BriDe+7QWA/sOOJWzRlgQ1XdMVBv5xHHmOm3SZKkCZp2iLiCbiRiWIC7+j+vA1YAew/VGZ4DsY6huQ9J9gBWDtWTJEkTMO0Q8RHgyUl2Gyh7BrA98NX+/WeBm4EjZyskWQkcBlw08LmLgOckedBA2VHAbXSXRyRJ0gRNO0ScA/wIuDDJYUl+E3gfsKaqPgNQVbcDZwAnJXllkoOBD/ZtffvAvs4GfgJckOSQJEcDpwJvdY0ISZImb6oTK6vq5iQHAX9Kt6bDHcCHgFcNVT2DLjScCOwKrAWeVVU/GNjX+j5gnEW3hsRNwNvogoQkSZqwqS97XVX/BjxvC3UKWN2/Nlfvcrq7OiRJ0pT5FE9JktTEECFJkpoYIiRJUhNDhCRJamKIkCRJTQwRkiSpiSFCkiQ1MURIkqQmhghJktTEECFJkpoYIiZs1eo1rFq9ZqGbIUnSvDNESJKkJlN/ANe2zBEHSdJy4kiEJElqYoiQJElNDBGSJKmJIUKSJDUxREiSpCaGCEmS1MQQIUmSmhgiJElSE0OEJElqYoiQJElNDBGSJKmJIUKSJDUxREiSpCaGCEmS1MQQIUmSmhgiJElSE0OEJElqYoiQJElNDBGSJKmJIUKSJDUxREiSpCaGCEmS1MQQIUmSmhgiJElSE0OEJElqYoiQJElNDBGSJKmJIUKSJDUxREiSpCaGCEmS1MQQIUmSmhgiJElSE0OEJElqYoiQJElNDBGSJKmJIUKSJDUxREiSpCaGCEmS1MQQIUmSmkw9RCR5WZIa8TpmoE6SnJTku0luS/LpJPuP2Ne+SS5NsiHJ95OcnmTFdL+RJEnL03YLeOyDgNsG3l898OcTgFOA1wDrgOOBNUn2q6rrAZLMAGuAy4HDgccCb6ELRidPvPWSJC1zCxkivlBVtwwXJnkAXYh4Q1Wd1Zd9DrgGOJZ7AsIxwAOBI6rqZuCSJDsBpyZ5Y18mSZImZDHOiXgasBNw/mxBVd0KXAgcOlDvUODiobBwHl2wOHAK7ZQkaVlbyBDxzSQ/TXJlkt8dKN8H2AhcNVT/in7bYL11gxWq6jvAhqF6kiRpAhbicsZ1dPMdPg+sAF4MnJ1kZVW9DZgBbqmqjUOfWw+sTLJDVd3R17tpxP7X99vuI8nRwNEAj3rUo+bju0iStGxNPURU1cXAxQNFF/XzIE5O8icTPvY5wDkABxxwQE3yWJIkbesWy5yIvwZ2AfakG0nYccStmjPAhn4Ugr7eziP2NdNvkyRJE7RYQkQN/FxHd5lj76E6w3Mg1jE09yHJHsDKoXqSJGkCFkuIeBFwI/Bt4LPAzcCRsxuTrAQOAy4a+MxFwHOSPGig7Ci6tScum3SDJUla7qY+JyLJ39BNqvwa3YjDUf3ruKq6C7g9yRnAKUnWc89iU/cD3j6wq7OB44ALkpwJ7AWcCrzVNSIkSZq8hbg740rg5cAeQOhWnHxpVb1voM4ZdKHhRGBXYC3wrKr6wWyFqlqf5GDgLLo1JG4C3kYXJCRJ0oQtxN0ZJwEnbaFOAav71+bqXU63fLYkSZqyxTInQpIkLTGGCEmS1MQQIUmSmhgiJElSE0OEJElqYoiQJElNDBGSJKmJIUKSJDUxREiSpCaGCEmS1MQQIUmSmhgiJElSE0OEJElqYoiQJElNDBFTsmr1GlatXrPQzZAkad4YIiRJUhNDhCRJamKIkCRJTQwRkiSpiSFCkiQ1MURIkqQmhghJktTEECFJkpoYIiRJUhNDhCRJamKIkCRJTQwRkiSpiSFCkiQ1MURIkqQmhghJktTEECFJkpoYIiRJUhNDhCRJamKIkCRJTQwRkiSpiSFCkiQ1MURIkqQmhghJktTEECFJkpoYIiRJUhNDhCRJarLdQjdgW7Bq9ZqFboIkSVPnSIQkSWpiiJAkSU0MEZIkqYkhYspWrV7jHApJ0jbBECFJkpoYIiRJUhNDhCRJamKIkCRJTQwRkiSpyYKGiCSPSHJLkkqy40B5kpyU5LtJbkvy6ST7j/j8vkkuTbIhyfeTnJ5kxXS/hSRJy9NCj0S8CbhlRPkJwCnAmcBhfZ01SR46WyHJDLAGKOBw4HTg1cBpE26zJEliAUNEkmcAzwXePFT+ALoQ8YaqOquq1gBH0oWFYweqHgM8EDiiqi6pqrPpAsTxSXaaxneQJGk5W5AQ0V9yeDvd6MGNQ5ufBuwEnD9bUFW3AhcChw7UOxS4uKpuHig7jy5YHDiBZkuSpAELNRJxDHB/4M9GbNsH2AhcNVR+Rb9tsN66wQpV9R1gw1A9SZI0AVMPEUl2Bf4QOL6q7hxRZQa4pao2DpWvB1Ym2WGg3k0jPr++3zbq2EcnWZtk7Q033ND2BSRJErAwIxGrgX+qqo9O+8BVdU5VHVBVB+y+++7TPrwkSduU7aZ5sCRPAl4OPCPJz/bFK/ufOyfZSDeSsGOSFUOjETPAhqq6o3+/Hth5xGFm+m2SJGmCphoigMcB2wOfG7HtWuDdwPuBFcDewJUD24fnQKxjaO5Dkj3oQsm95kpIkqT5N+0Q8RngmUNlzwVeCzwPuBr4NnAz3W2drwdIspJuvYhzBj53EfCaJA+qqv/oy44CbgMum9QXkCRJnamGiKq6EfjUYFmSPfs//kNV3dKXnQGckmQ93ajC8XTzN94+8NGzgeOAC5KcCewFnAq8dei2T0mSNAFjhYgkH6C75HBJVdVkmgTAGXSh4URgV2At8Kyq+sFshapan+Rg4Cy6NSRuAt5GFyQkSdKEjTsS8QjgY8D3krwHOLeq/m1rGlBV5wLnDpUV3V0cq7fw2cuBg7bm+JIkqc1Yt3hW1dOBJwDvA14KXNk/HOtlSX5mEg2UJEmL09jrRFTVVVV1EvBousmQ19KtPHldkncnefo8t1GSJC1CzYtN9ZccLqO7S+LrwI50oeLTSb6Y5Ofnp4mSJGkxagoRSX4lybuA6+numPgK8NSqehiwP90tmu+dt1ZKkqRFZ9y7M04CfptuIajPAa8CPlBVG2brVNXXkpwMfHo+GypJkhaXce/OOI5uhOHdVXXlZuqtA45ubpUkSVr0xg0Rj6yqn26pUlX9iG49CUmStI0ad07E05O8dNSGJP9vkgPnoU2SJGkJGDdE/BHw8E1se2i/XZIkLQPjhoj96JagHuVLwJO2rjmSJGmpGDdE3AXMbGLbrg37kyRJS9S4J/1/BF6dZPvBwv79q+ge9S1JkpaBce/OOIkuKHwjyXnAdcDDgBcDuwC/Or/NkyRJi9VYIaKqvprkl+ket/0KuuDwY+BS4A+qat28t1CSJC1K445EUFVfB46cQFskSdIS4kRISZLUZOyRiCQvBI4AHgk8YHh7VT1tHtolSZIWuXEfwHUKcBrdo78vB+6YRKMkSdLiN+5IxNHAm6rqtZNojCRJWjrGnRPxIODjk2iIJElaWsYNEecDz55EQyRJ0tIy7uWMjwFvTrILcAlw03CFqnKkQpKkZWDcEPHX/c//2r+GFbBiq1okSZKWhHFDxOMm0gpJkrTkjLvs9Tcn1RBJkrS0jL1iZZLtk7wiyTuTfDTJ3n35i5I8Yf6bKEmSFqNxF5vam+4Wz92AL9E9tXOnfvMzgcOA357PBkqSpMVp3JGIPwWuB/YEDgEysO0yfBS4JEnLxrgTKw8EfqOqfpxk+C6M64GHzU+zJEnSYjfuSMRPgPtvYtvDGbFuhCRJ2jaNGyIuAU5M8qCBskqyPXAs3WJUkiRpGRj3csZrgM8C/wZcTLe41OuAJwE/A/zGvLZOkiQtWmONRFTVd4CfB/4C2Af4Nt0kyw8Dv1hV35/vBkqSpMVp3JEIqupHwIkTaIskSVpCxl5sSpIkCcZfbOo6unkQm1RVD9+qFkmSpCVh3MsZ7+a+IWIGOBhYCbxnPholSZIWv3EfwHXyqPIk9wM+CGyYj0ZJkqTFb17mRFTVXcC7gOPmY3+SJGnxm8+JlY8GdpjH/UmSpEVs3ImVR48o3gF4IvBS4IL5aJQkSVr8xp1YefaIsp8C36O7nPH7W90iSZK0JIwbIrYfLqiqjfPUFkmStISMe3eGgUGSJAHjz4n4zXHqV9X7x2uOJElaKsa9nPGX3LPYVAbKN1VmiJAkaRs17i2eT6F7cudpwM8BD+1/nt6XP4VuBcsZYJf5a6YkSVpsxh2JOBP486p600DZD4F/TbIBeGNVPXPeWidJkhatcUcifhn46ia2fY1uJEKSJC0D44aIa4GXbWLby+jWi5AkScvAuJczTgben2Rf4MN0lzIeDLwAeDLwkvlt3rZr1eo1AHz+dYcscEskSWoz7joR5ye5BjgB+B3gIcAPgC8Av1tV/zzvLZQkSYvS2A/gqqrPV9URVbVHVe3Q/zxiLgEiyYuSfDbJj5LcnuTKJCcn2WGgTpKclOS7SW5L8ukk+4/Y175JLk2yIcn3k5yeZMW430eSJLUZ93IGAEl2BvYF9gA+XlU3Jdm+qu7cwkd3BT4BvAm4CVgFnEp3q+ixfZ0TgFOA1wDrgOOBNUn2q6rr++PPAGuAy4HDgccCb6ELRSe3fCdJkjSecVesvB/weuA4YCXdglK/BHwJ+HCSf6qq0zb1+ap651DRJ5PsBLwyye8B96cLEW+oqrP6Y34OuIYuZMwGhGOABwJHVNXNwCX9fk5N8sa+TJIkTdC4lzNWA68EXgU8nnuvUPl3dBMsx/UjuseJAzwN2Ak4f3ZjVd0KXAgcOvCZQ4GLh8LCeXTB4sCGNkiSpDGNGyJ+Gzihqt4FfGto2zfpLitsUZIVSVYmeTrdqMafV1UB+wAbgauGPnJFv23WPnSXOu5WVd8BNgzVkyRJEzLunIgZ7nuCn7U9MNeJjbfSXboAeC/d/IfZ/d8y4mmh64GVSXaoqjv6ejeN2O/6fttISY4GjgZ41KMeNcemSpKkUcYdifg6cNgmtj0H+PIc9/M04FeBV9NNjDxrzHY0qapzquqAqjpg9913n8YhJUnaZo07EvFHwPlJ7g98kG5i5X5JDgP+G/DCueykqr7U//EzSW4E3pPkLXQjCTsmWTE0GjEDbOhHIejr7Txi1zP9NkmSNGFjjURU1QXAS4H/B7iEbmLlucDvAr9TVRc1tGE2UDyGbp7DCmDvoTrDcyDWMTT3IckedHeM3GuuhCRJmoyWxabeDzwK2A/4NbpHgT+yL2/xK/3PbwGfBW4GjpzdmGQl3SWUwYByEfCcJA8aKDsKuA24rLEdkiRpDHO+nJHkAXSjBq+qqovpFnoaS5KP0S0S9XW6uzB+hW5exAeq6pt9nTOAU5Ks557Fpu4HvH1gV2fT3dVxQZIzgb3oFq16q2tESJI0HXMOEVV1e5Ld6OZBtPoC3dM+9wR+ClwNnEgXCmadQRcaTqRb4XIt8Kyq+sFAW9YnOZhuQuaFdHdqvI0uSEiSpCkYd2LlX9HNifh4y8Gq6hS6Ja03V6foFrVavYV6lwMHtbRDkiRtvXFDxDeBFyX5J+CjdE/wHByZqH4hKkmStI0bN0T8cf/zYXQPzxpWgCFCkqRlYNwQsf1EWiFJkpacLd7imeTjSZ4AUFUb+0WgDgQeMPt+8DXpBkuSpMVhLutEHMLA6pBJVtAtNPWESTVKkiQtfmMvNtXLlqtIkqRtWWuIkCRJy9xcQ8SoBaa2ZtEpSZK0xM317oyLk/x0qOzSEWVU1YO3vlmSJGmxm0uIOG3irZAkSUvOFkNEVRkiJEnSfTixUpIkNTFESJKkJoYISZLUxBAhSZKaGCIkSVITQ4QkSWpiiJAkSU0MEZIkqYkhQpIkNTFESJKkJoYISZLUxBAhSZKaGCIkSVITQ4QkSWpiiJAkSU0MEZIkqYkhYoGtWr2GVavXLHQzJEkamyFCkiQ1MURIkqQmhghJktTEECFJkpoYIiRJUhNDhCRJamKIkCRJTQwRkiSpiSFCkiQ1MURIkqQmhghJktTEELFI+AwNSdJSY4iQJElNDBGSJKmJIUKSJDUxREiSpCaGCEmS1MQQIUmSmhgiJElSE0OEJElqYoiQJElNDBGSJKmJIUKSJDWZaohIcmSSDyf5XpJbknwxyUtG1HtFkquS3N7XOXhEnUck+dsk/5HkxiRnJVk5nW8iSZKmPRJxPHAL8CrgBcAngfcn+b3ZCn2oOBt4L3Ao8HXgI0n2G6izPXAx8GjgxcD/AI4EzpnO15AkSdtN+XiHVdWNA+8/keThdOHi7X3ZqcB7quoPAZJcBvwCcALwX/o6LwKeCOxdVd/q690JnJfktKq6auLfRJKkZW6qIxFDAWLWl4GHAyTZC3g8cP7AZ+4CPkg3KjHrUOALswGi93fAHcBz57nZkiRphMUwsfKpwDf6P+/T/1w3VOcKYJckuw/Uu1edqroD+ObAPiRJ0gQtaIjoJ0y+EHhLXzTT/7xpqOr6oe0zI+rM1psZUS5JkubZgoWIJHsC7wc+VFXnTumYRydZm2TtDTfcMI1DSpK0zVqQEJFkF+Ai4NvAbw1smh1x2HnoIzND29ePqDNbb/2IcgCq6pyqOqCqDth99903VU2SJM3B1ENEv5bDR4AdgOdX1YaBzbPzHIbnNewD/Liqbhiod686SXYA9uK+8ykkSdIETHuxqe3o7rR4HPDcqvrh4PaquppukuWRA5+5X//+ooGqFwG/lOTRA2UvAO4PfGwyrZckSYOmvU7EO4Dn0S0OtWuSXQe2fbmqfkK3TsRfJrkG+Efgt+lCx28O1P1r4HXABUlOobu08Tbg/a4RIUnSdEw7RDy7//knI7Y9Brimqv4qyY7Aa4FT6FasfH5V/etsxaq6M8lzgbPo1pT4CXAe8JpJNl6SJN1jqiGiqvacY713Ae/aQp1r6W4PlSRJC2AxLDYlSZKWIEOEJElqYoiQJElNDBGSJKmJIUKSJDUxREiSpCaGCEmS1MQQIUmSmhgiJElSE0OEJElqYoiQJElNDBGSJKmJIUKSJDUxREiSpCaGCEmS1MQQIUmSmmy30A3Qva1avebuP3/+dYcsYEskSdo8RyIkSVITQ4QkSWpiiJAkSU0MEZIkqYkhQpIkNTFESJKkJoYISZLUxBAhSZKaGCIkSVITQ4QkSWpiiFjEVq1ec69lsCVJWkwMEZIkqYkhQpIkNTFESJKkJoYISZLUxBAhSZKaGCIkSVITQ4QkSWpiiJAkSU0MEZIkqYkhQpIkNTFESJKkJoYISZLUxBAhSZKaGCIkSVITQ4QkSWpiiJAkSU0MEZIkqYkhQpIkNTFELAGrVq9h1eo1C90MSZLuxRAhSZKaGCIkSVITQ4QkSWpiiJAkSU0MEZIkqcnUQ0SSvZO8M8nXkmxM8qkRdZLkpCTfTXJbkk8n2X9EvX2TXJpkQ5LvJzk9yYqpfJEF4F0akqTFZCFGIp4EPA+4EvjGJuqcAJwCnAkcBtwCrEny0NkKSWaANUABhwOnA68GTptYyyVJ0t0WIkRcWFV7VNWRwNeHNyZ5AF2IeENVnVVVa4Aj6cLCsQNVjwEeCBxRVZdU1dl0AeL4JDtN/FtIkrTMTT1EVNVdW6jyNGAn4PyBz9wKXAgcOlDvUODiqrp5oOw8umBx4Py0VpIkbcpinFi5D7ARuGqo/Ip+22C9dYMVquo7wIahepIkaQIWY4iYAW6pqo1D5euBlUl2GKh304jPr++33UeSo5OsTbL2hhtumLcGS5K0HC3GEDExVXVOVR1QVQfsvvvuC90cSZKWtMUYItYDO464VXMG2FBVdwzU23nE52f6bdssb/WUJC0GizFErANWAHsPlQ/PgVjH0NyHJHsAK4fqSZKkCViMIeKzwM10t3UCkGQl3XoRFw3Uuwh4TpIHDZQdBdwGXDaFdkqStKxtN+0D9oHgef3bRwA7JXlR//6jVbUhyRnAKUnW040qHE8XeN4+sKuzgeOAC5KcCewFnAq8dei2T0mSNAFTDxHAg4EPDpXNvn8McA1wBl1oOBHYFVgLPKuqfjD7gapan+Rg4Cy6NSRuAt5GFyQkSdKETT1EVNU1QLZQp4DV/Wtz9S4HDpq3xkmSpDlbjHMiJEnSEmCIkCRJTQwRkiSpiSFCkiQ1MURIkqQmhghJktTEECFJkpoYIiRJUhNDhCRJamKIkCRJTQwRkiSpiSFCkiQ1WYineGqerFq95l7vP/+6QxaoJZKk5ciRCEmS1MQQIUmSmhgiJElSE0OEJElq4sTKBsMTGiVJWo4ciZAkSU0MEZIkqYkhQpIkNTFESJKkJoYISZLUxBCxjVu1eo13k0iSJsIQIUmSmrhOxDZkcMTBh3FJkibNkQhJktTEECFJkpoYIiRJUhNDhCRJauLEym3U8G2dw++deClJ2lqOREiSpCaGCEmS1MQQIUmSmhgiJElSE0OEJElqYoiQJElNvMVzmRv1vI3ZMm8DlSRtjiFimRr1ePBNPTLcUCFJGsXLGRrbqtVrNhk4JEnLhyFCkiQ18XKGNmkSow1zuTTi5RNJWhocidC88lKHJC0fjkRozgwHkqRBjkRIkqQmhggtWgt1acRLMpI0N17OULNRC1UNbxtewGqu+xu3DVuahOlkTUmaf4YIzYstLVS1ufK5ntgNApK0uBgitOiNM4qxpaW7R+1rEqHEwCNpOTBEaJviXAZJmh5DhBbcuCf+ceq3jGJMwjiXb6Y9iuGoiaRWhgipt5RGMbZ0uWYugWA+v69BRFqelnSISLIv8HbgqcBNwP8BTquqjQvaMG0T5nKSHefuk605zpY+O05o2FTdcZYkn++2SVqalmyISDIDrAEuBw4HHgu8hW7ti5MXsGnSZk3i5Dp8cp/L7bfj7Hc+2jofgWe47aP2Ncnwsrl9txx3ofpXmi9LNkQAxwAPBI6oqpuBS5LsBJya5I19mTRR07oEMs2RjlH157qP+V4PZGtOjFtaq2RrgsDW3KI8TltbRoXmEsDmOhplMNGWLOUQcShw8VBYOA84EzgQuHBBWiXN0ZZOpktpjsbmLOXvsdhOovPVl4vte2npSlUtdBuaJPkh8I6qOnWo/Fbg1Kp60+Y+f8ABB9TatWvHOuZS/mUoLWXzMfdkoSzlts8aDBsGkKUryRer6oB53ecSDhF3Aq+pqj8eKr8WeG9VnTTiM0cDR/dvnwBcuZXN2A24cSv3ofuyXyfDfp0M+3X+2aeT8YSqetB87nApX84YW1WdA5wzX/tLsna+U53s10mxXyfDfp1/9ulkJBlv+H0OlvJTPNcDO48on+m3SZKkCVrKIWIdsM9gQZI9gJX9NkmSNEFLOURcBDwnyeD1naOA24DLptSGebs0onuxXyfDfp0M+3X+2aeTMe/9upQnVs7QLTT1r3S3de4FvBX446pysSlJkiZsyYYIuHvZ67O497LXp7rstSRJk7ekQ4QkSVo4S3lOxIJIsm+SS5NsSPL9JKcnWbHQ7VoqkhyZ5MNJvpfkliRfTPKSEfVekeSqJLf3dQ5eiPYuVUke0fdvJdlxoDxJTkry3SS3Jfl0kv0Xsq2LXZLtkpzQ/338SZJrk7xtqI79OqYkL07ypf7v6feSvDfJw4fq2K+bkWTvJO9M8rUkG5N8akSdOfVh67nNEDGGgYd+Fd1Dv04HXg2ctpDtWmKOB24BXgW8APgk8P4kvzdboQ8VZwPvpVve/OvAR5LsN/3mLllvouvnYScAp9DNIzqsr7MmyUOn2Lal5lzgOODNwLPp+vC2oTr26xiSvAD4K+CzdL9LXws8A/j7JIPnJft1854EPI9u4cRvbKLOFvtwq85tVeVrji/gRLo1KHYaKPvfwIbBMl+b7cPdRpS9H/jWwPsrgb8YeH8/4F+Av1zo9i+FF9388xG1AAAKf0lEQVQv4x8D/6v/pbBjX/4A4N+B3x+o+zPADcDrF7rdi/EFPBe4E9h3M3Xs1/H79Tzgi0NlL+j/vj7Rfp1zP95v4M9/DXxqaPuc+nBrzm2ORIxnUw/9eiDdQ7+0BVU1ainbLwMPB0iyF/B44PyBz9wFfJCu/7UZ/fDj2+n+JTHc108DduLefXsr3cPq7NvRXg58oqou30wd+3V829Od3Abd1P9M/9N+3YL+d+PmzLUPm89thojx7MPQQlZV9R26tLbPyE9oLp7KPUNxs/04vGDYFcAuSXafWquWpmOA+wN/NmLbPsBG4Kqh8ivw7++mPAX4RpKzktzcXy++YOjavf06vr8AfjXJS5PslOTxwOu5d2CzX7feXPuw+dxmiBjPDPek5UHr+20aUz9h8oXAW/qi2X4c7uf1Q9s1JMmuwB8Cx1fVnSOqzAC31H1vgV4PrEyyw6TbuAQ9FHgZsD/wYuB3gF8E/jbJ7L+Y7dcxVdXf0/XrOXQjElcCK4BfH6hmv269ufZh87ltWT2AS4tLkj3p5kN8qKrOXdDGbBtWA/9UVR9d6IZsQ9K/Dq+qHwEkuY5uVdyDgEsXsG1LVpJn0k2e/hO61YcfApxKF84OGXHS0yJliBiPD/2aJ0l2ofvl8W3gtwY2zfbjztw7Gc8MbdeAJE+iu37/jCQ/2xev7H/unGQjXd/tmGTF0C/pGWBDVd0xvRYvGeuBq2cDRO8zwB3AvnQhwn4d31uAD1fVa2cLknyFbkj9cOAC7Nf5MNc+bD63eTljPD70ax4kWQl8BNgBeH5VbRjYPNuPw9fh9gF+XFU3TKGJS9Hj6CarfY7uf/r13DMv4lq6yZbr6IaM9x767H2uh+puV3DPRL9BAWYntdmv49sH+MpgQVVdSXfr7GP7Ivt16821D5vPbYaI8SyGh34taUm2o7vT4nHAc6vqh4Pbq+pqukmWRw585n79+4um2NSl5jPAM4deZ/bbnke3bsRngZu5d9+upLt33L4d7SPAk5PsNlD2DLrA9tX+vf06vm8D/2mwIMkT6e4GuKYvsl+33lz7sP3cttD3uS6lF93QznXAJcAhwNF0C3d4z/Lc+/AcunvBjwN+eeh1/77OS+hmFJ9MdzI8t//LvN9Ct38pvegmrt29TkRfdiLdjOtXAgcDf093K+hDFrq9i/FFd3vcd+hGeA4DfhP4LnDJUD37dbx+/R90Izlv6X+X/hbd5MpvAT9jv865H1cCL+pfn6NbmG/2/cq59uHWnNsWvBOW2ovuOugn+pPadXSz4VcsdLuWyovuXxm1ideeA/VeAfwb8BPgS8DBC932pfbaRIgI8Dq6Sxy3Af8A/MJCt3Uxv+iGgj8K3Ep3mehcYGaojv06Xp8G+G/A1/p+/R7wAWAv+3WsftxzS79P59qHrec2H8AlSZKaOCdCkiQ1MURIkqQmhghJktTEECFJkpoYIiRJUhNDhCRJamKIkKYsycuSfDHJfyRZn+TLSd46sH3PJJXk+VNqTyU5dhrHmpQkz07yP+dY95okb550m6TlwBAhTVGSE4H/A1wMHAG8FPgQ8IKFbNc24NnAnEKEpPnjUzyl6ToWeGdVnTRQdmGS0xaqQZquJKFb4v32hW6LtLUciZCm62eB64cLa/TSsSuTvDPJvye5Nslp/cPI7pbkoCT/nOT2JD9I8o4kOw7V2bXfz3V9vSs3N/SfZL8k1yd5X5IVm6m32WMn+bX+UsmvJflgkluSXJ3kvw/t50lJPpbkx0luTXJFklcO1Tk8ydr+WNcneWOS7fttpwKvBh7dH6+SnLupdg/s81V9v65Pct7AI9Rntz8myd8lubm/9HRhkr0Hto+87JTk3CRrB96fmuTGJE9P8gXgdgYeiCQtZY5ESNP1JeD3knwH+EhV/Wgzdd8I/A3dw3QOBn6f7gE750N38gU+RvfQnF8H9gDOAPYCntvXeSDwKeDBwGl0j/Xdm/s+Gpi+/i/0+7sAOKaq7tpEvS0ee8C7gPfQPXztJcCfJVlbVZ/vt19I98jt/0L3rJQn0D34avZYvwH8FfBO4CS6R0W/ge4fQf+L7vLQ44CDgP/cf2xLj4z/DbrnNhwNPBJ4K/BHwH/vj3l/4FLgTrrnuPyUrv8uS/LkqvrxFvY/bGXfB2+ke0rt98f8vLQ4LfQDRHz5Wk4v4OeAq+kekHMXXSg4HdhpoM6e/fb3Dn32K8B5A+/PA65i4CE5dCfHAp7av//d/jj7b6ZNRXeZ5Sl0D5j6U+ieq7OZz8zl2L/Wvz99oM72dCf4M/r3u/V1nryJ44TusdH//1D5y+keFLRr//7NwDVz/G9wDfBNYLuBsj8Grh94fwxdcNhroOyRwB3AiUP/nZ4/tP9zgbUD70/t6x2+0H//fPma75eXM6QpqqqvAU+km0j5DrqT5CnA2uHLEMDHh95fTncim7UK+Nuq2jhQ9jd0J7+n9+8PAr5cVV/ZQtN+hW5U4ZyqOq6qtvRkvrkc+z7fo6rupAsfs9/jx3SP1j47yVFJHjz02ccDjwLOT7Ld7IvuaYMPAPbbQjs35ZNV9dOB95cDD569RNJ/vy9V1dUDbb8W+McR328uCriosa3SomWIkKasqn5SVRdW1bFVtS/w/9ENx//Xoao3Db2/g+7EOethwA+G9r0R+BGwS1+0K91jfbfk2XSXN987py8xt2PP2uT3qO5yybPp5on8BXB9kn/oL6tAN1IB3aO47xx4fasv32OO7R02qk0B7t+/v8/36/2A+36/uVhfVXc0fE5a1AwR0gKrqnfT/Yt8nzE/eh3dXIe79RMhd+33B91J/WFz2NfrgU8CH0+y1zwde06qal1V/TrdpNND6ALG3/eTSGf3dTTwSyNek/rX/X2+X+8hA22avbtih6E6MyM+t6WRHWlJMkRIUzRiuJ4kuwM7M/pfvpvzz8B/HrqD4gi6EYXP9O8vBX4hyc9tYV930k3g/AZwaZJHzMOxx1JVd1bVJ+gmOT6MLlRcCXwP2LOq1o54zU5MHR6l2Vr/DPxiksfMFvR98jTu+X4/pOu3Jw7U2bGvIy0L3p0hTde/JPkQ3TyBHwKPprvDYAPd7P1xvB74MvB3Sf6cbp7BmcDFVfW5vs57gVfSjTCcSndSfgzw+Ko6YXBnVXVbksOANcCaJM+oqk3d5TCXY29RH27eDHyAbsLpDPBa4KvV3wGR5NXA+5LsRDfycAfdXSAvBF5UVRvo7jp5SJKXAf8K3FhV18y1HSOc27fjoiS/D2wE/gC4ke4uEarqrv6/5auSfJvuEsmr6SZ8SsuCIUKartOBw+nugNiFbi7AZ4Gjqupbm/vgsKr6epJD6W5NvAC4me5WyP89UOf2JAfR3X55Ot2tk9fQTeoctc9b+n1+Erg4yTOr6t9bjj1H19ONwLwOeDjdifiTdCfw2WN9IMnNdLd3vpzuhH418BG6QAHdba/PpLuFcne6QPayMdtyt6r6SZJD6EZF3k03X+JTwK/XvW/vPJbu1tV30N3ZsppuJKJ1wqe0pGTLk7AlSZLuyzkRkiSpiSFCkiQ1MURIkqQmhghJktTEECFJkpoYIiRJUhNDhCRJamKIkCRJTf4vT96t7Ns0Xy8AAAAASUVORK5CYII=\n",
+      "text/plain": [
+       "<Figure size 576x576 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "visualize_labels(df_ARF, 'ARF')\n",
+    "visualize_labels(df_Shock, 'Shock')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "---"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from datetime import timedelta\n",
+    "cutoff_h = 4\n",
+    "mimic3_path = yaml.full_load(open('config.yaml'))['mimic3_path']"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "patients = pd.read_csv(mimic3_path + 'PATIENTS.csv', parse_dates=['DOB', 'DOD'], usecols=['SUBJECT_ID', 'DOB', 'DOD'])\n",
+    "admissions = pd.read_csv(mimic3_path + 'ADMISSIONS.csv', parse_dates=['DEATHTIME'], usecols=['SUBJECT_ID', 'HADM_ID', 'DEATHTIME'])\n",
+    "examples = pd.read_csv(data_path + 'prep/icustays_MV.csv', parse_dates=['INTIME', 'OUTTIME']).sort_values(by='ICUSTAY_ID') # Only Metavision\n",
+    "\n",
+    "examples = pd.merge(examples, patients, on='SUBJECT_ID', how='left')\n",
+    "examples = pd.merge(examples, admissions, on=['SUBJECT_ID', 'HADM_ID'], how='left')\n",
+    "examples['AGE'] = examples.apply(lambda x: (x['INTIME'] - x['DOB']).total_seconds(), axis=1) / 3600 / 24 / 365.25\n",
+    "\n",
+    "examples['LOS'] = examples['LOS'] * 24 # Convert to hours"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/html": [
+       "<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>LOS</th>\n",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th>count</th>\n",
+       "      <td>23620.000000</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>mean</th>\n",
+       "      <td>3.593499</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>std</th>\n",
+       "      <td>4.971162</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>min</th>\n",
+       "      <td>0.000400</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>25%</th>\n",
+       "      <td>1.151600</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>50%</th>\n",
+       "      <td>1.996050</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>75%</th>\n",
+       "      <td>3.835000</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>max</th>\n",
+       "      <td>101.739000</td>\n",
+       "    </tr>\n",
+       "  </tbody>\n",
+       "</table>\n",
+       "</div>"
+      ],
+      "text/plain": [
+       "                LOS\n",
+       "count  23620.000000\n",
+       "mean       3.593499\n",
+       "std        4.971162\n",
+       "min        0.000400\n",
+       "25%        1.151600\n",
+       "50%        1.996050\n",
+       "75%        3.835000\n",
+       "max      101.739000"
+      ]
+     },
+     "execution_count": 9,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "(examples[['LOS']] / 24.).describe()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Remove non-adults\n",
+    "min_age = 18\n",
+    "max_age = np.inf # no max age\n",
+    "examples = examples[(examples.AGE >= min_age) & (examples.AGE <= max_age)]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "23593"
+      ]
+     },
+     "execution_count": 11,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "examples['ICUSTAY_ID'].nunique()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/html": [
+       "<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>ICUSTAY_ID</th>\n",
+       "      <th>DEATHTIME</th>\n",
+       "      <th>INTIME</th>\n",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th>10</th>\n",
+       "      <td>200038</td>\n",
+       "      <td>NaT</td>\n",
+       "      <td>2143-10-24 20:35:24</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>11</th>\n",
+       "      <td>200040</td>\n",
+       "      <td>NaT</td>\n",
+       "      <td>2153-10-24 16:01:41</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>12</th>\n",
+       "      <td>200049</td>\n",
+       "      <td>NaT</td>\n",
+       "      <td>2118-08-28 08:56:44</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>13</th>\n",

mimic3_experiments/1_data_extraction/PopulationSummary.ipynb

+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import yaml\n",
+    "with open('../config.yaml') as f:\n",
+    "    config = yaml.full_load(f)\n",
+    "\n",
+    "data_path = config['data_path']\n",
+    "mimic3_path = config['mimic3_path']\n",
+    "\n",
+    "import pandas as pd\n",
+    "import itertools\n",
+    "from collections import Counter"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "icustays = pd.read_csv(data_path + 'prep/icustays_MV.csv')\n",
+    "partition = icustays.set_index('ICUSTAY_ID')[['partition']]\n",
+    "tasks = ['ARF', 'Shock']\n",
+    "Ts = [4, 12]\n",
+    "\n",
+    "populations = {}\n",
+    "for task, T in itertools.product(tasks, Ts):\n",
+    "    pop = pd.read_csv(data_path + 'population/{}_{}h.csv'.format(task, T))\n",
+    "    populations[task, T] = pop.set_index('ICUSTAY_ID')[['{}_LABEL'.format(task)]]\n",
+    "\n",
+    "populations['mortality', 48] = pd.read_csv(data_path + 'population/pop.mortality_benchmark.csv'.format('mortality', 48)) \\\n",
+    "                                 .set_index('ID')[['{}_LABEL'.format('mortality')]]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "df_out = []\n",
+    "for (task, T), labels in populations.items():\n",
+    "    df = labels.join(partition)\n",
+    "    c = Counter(df['partition'])\n",
+    "    frac = df.groupby('partition').mean()['{}_LABEL'.format(task)]\n",
+    "    df_out.append([task, T, \n",
+    "                   len(df),    df['{}_LABEL'.format(task)].mean(),\n",
+    "                   c['train'], frac['train'], \n",
+    "                   c['val'],   frac['val'], \n",
+    "                   c['test'],  frac['test']])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "df_out = pd.DataFrame(df_out, columns=['task', 'T', 'TOTAL_N', 'TOTAL_%', 'train_N', 'train_%', 'val_N', 'val_%', 'test_N', 'test_%'])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/html": [
+       "<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>task</th>\n",
+       "      <th>T</th>\n",
+       "      <th>TOTAL_N</th>\n",
+       "      <th>TOTAL_%</th>\n",
+       "      <th>train_N</th>\n",
+       "      <th>train_%</th>\n",
+       "      <th>val_N</th>\n",
+       "      <th>val_%</th>\n",
+       "      <th>test_N</th>\n",
+       "      <th>test_%</th>\n",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th>0</th>\n",
+       "      <td>ARF</td>\n",
+       "      <td>4</td>\n",
+       "      <td>15873</td>\n",
+       "      <td>0.182700</td>\n",
+       "      <td>11147</td>\n",
+       "      <td>0.182291</td>\n",
+       "      <td>2368</td>\n",
+       "      <td>0.180743</td>\n",
+       "      <td>2358</td>\n",
+       "      <td>0.186599</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>1</th>\n",
+       "      <td>ARF</td>\n",
+       "      <td>12</td>\n",
+       "      <td>14174</td>\n",
+       "      <td>0.096515</td>\n",
+       "      <td>9971</td>\n",
+       "      <td>0.097282</td>\n",
+       "      <td>2110</td>\n",
+       "      <td>0.093365</td>\n",
+       "      <td>2093</td>\n",
+       "      <td>0.096034</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>2</th>\n",
+       "      <td>Shock</td>\n",
+       "      <td>4</td>\n",
+       "      <td>19342</td>\n",
+       "      <td>0.149416</td>\n",
+       "      <td>13613</td>\n",
+       "      <td>0.148241</td>\n",
+       "      <td>2862</td>\n",
+       "      <td>0.157582</td>\n",
+       "      <td>2867</td>\n",
+       "      <td>0.146843</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>3</th>\n",
+       "      <td>Shock</td>\n",
+       "      <td>12</td>\n",
+       "      <td>17588</td>\n",
+       "      <td>0.077553</td>\n",
+       "      <td>12381</td>\n",
+       "      <td>0.075923</td>\n",
+       "      <td>2595</td>\n",
+       "      <td>0.084008</td>\n",