global-model-1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Document dependencies\n",
    "Skip this section. To document the dependencies, execute the following cells after executing all relevant import statements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install --user watermark"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The watermark extension is already loaded. To reload it, use:\n",
      "  %reload_ext watermark\n",
      "\n",
      "CPython 3.7.4\n",
      "IPython 7.8.0\n",
      "\n",
      "pandas 0.25.1\n",
      "numpy  1.18.5\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%load_ext watermark\n",
    "print()\n",
    "%watermark -v\n",
    "print()\n",
    "%watermark -iv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Model 1: Predicting EVD(+) vs EVD(-)\n",
    "Replicate the model for predicting Ebola infection (EVD+ vs EVD-) from Annie's 2017 paper.\n",
    "It will serve along with other models, to evaluate the developed methods for model personalisation through data selection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn import model_selection as sk_ms, linear_model as sk_lm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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>Unnamed: 0</th>\n",
       "      <th>study_id</th>\n",
       "      <th>evd</th>\n",
       "      <th>evd_ct</th>\n",
       "      <th>malaria</th>\n",
       "      <th>death</th>\n",
       "      <th>sex</th>\n",
       "      <th>age</th>\n",
       "      <th>quarantine</th>\n",
       "      <th>days_admitted</th>\n",
       "      <th>...</th>\n",
       "      <th>hematuria_after</th>\n",
       "      <th>gyne_hemorrhage_after</th>\n",
       "      <th>num_hemorrhage_sites_after</th>\n",
       "      <th>dehydration_after</th>\n",
       "      <th>hepatomegaly_after</th>\n",
       "      <th>jaundice_after</th>\n",
       "      <th>edema_after</th>\n",
       "      <th>splenomegaly_after</th>\n",
       "      <th>referral(4-9)</th>\n",
       "      <th>quarantine_missing</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>0 rows × 70 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "Empty DataFrame\n",
       "Columns: [Unnamed: 0, study_id, evd, evd_ct, malaria, death, sex, age, quarantine, days_admitted, referral_time, date_admis, temp_triage, pyrexia_triage, anuria_triage, dyspnea_triage, disorientation_triage, headache_triage, hiccups_triage, myalgia_triage, rash_triage, dysphagia_triage, asthenia_triage, abdomnial_triage, anorexia_triage, diarrhea_triage, diarrhea_score_triage, vomiting_triage, hemorrhage_triage, orl_hemorrhage_triage, conjonctivitis_triage, hemoptysis_triage, hematochezia_triage, hematuria_triage, gynec_hemorrhage_triage, num_hemorrhages_triage, mult_hemorrhage_triage, dehydration_triage, hepatomegaly_triage, jaundice_triage, edema_triage, splenomegaly_triage, anuria_after, dyspnea_after, disorientation_after, headache_after, hiccups_after, myalgia_after, rash_after, dysphagia_after, asthenia_after, abdominal_after, anorexia_after, diarrhea_score_after, vomiting_after, hemorrhage_after, orl_hemorrhage_after, conjunctivitis_after, hemoptysis_after, hematochezia_after, hematuria_after, gyne_hemorrhage_after, num_hemorrhage_sites_after, dehydration_after, hepatomegaly_after, jaundice_after, edema_after, splenomegaly_after, referral(4-9), quarantine_missing]\n",
       "Index: []\n",
       "\n",
       "[0 rows x 70 columns]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Load the pre-processed dataset:\n",
    "df = pd.read_csv('../data/private/1_TreatMissing_EbolaDatabase.csv')\n",
    "df.head(0) # DO NOT PUSH A COMMIT WHERE ANY LINE IS VISIBLE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Build the feature matrix $\\mathbf{X}$ and the target vector $\\mathbf{y}$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['Unnamed: 0', 'study_id', 'evd', 'evd_ct', 'malaria', 'death',\n",
       "       'sex', 'age', 'quarantine', 'days_admitted', 'referral_time',\n",
       "       'date_admis', 'temp_triage', 'pyrexia_triage', 'anuria_triage',\n",
       "       'dyspnea_triage', 'disorientation_triage', 'headache_triage',\n",
       "       'hiccups_triage', 'myalgia_triage', 'rash_triage',\n",
       "       'dysphagia_triage', 'asthenia_triage', 'abdomnial_triage',\n",
       "       'anorexia_triage', 'diarrhea_triage', 'diarrhea_score_triage',\n",
       "       'vomiting_triage', 'hemorrhage_triage', 'orl_hemorrhage_triage',\n",
       "       'conjonctivitis_triage', 'hemoptysis_triage',\n",
       "       'hematochezia_triage', 'hematuria_triage',\n",
       "       'gynec_hemorrhage_triage', 'num_hemorrhages_triage',\n",
       "       'mult_hemorrhage_triage', 'dehydration_triage',\n",
       "       'hepatomegaly_triage', 'jaundice_triage', 'edema_triage',\n",
       "       'splenomegaly_triage', 'referral(4-9)', 'quarantine_missing'],\n",
       "      dtype='<U23')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array([c for c in df.columns if not 'after' in c])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, the prediction model is not useful if it relies on features collected upon release of the patients (by that time, the outcome is known).\n",
    "We therefore focus our efforts on the features collected during patient triage.\n",
    "\n",
    "Diarrhea is used in the model proposed in Annie's paper, and it seems to be a binary feature.\n",
    "Here, however, we have two features at our disposal:\n",
    "- `diarrhea_triage` is simply binary.\n",
    "- `diarrhea_score_triage` is an ordered variable, going from 0 to 3 in integer steps depending on the number of diarrhea episodes experienced by the patient.\n",
    "\n",
    "We plot `evd` against `diarrhea_score_triage` to identify whether a linear correlation exists:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7f1be382a110>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEICAYAAACktLTqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAW2klEQVR4nO3df5BV5X3H8fdHQDciYsSFUBYFJ1SkRDZ0paZUE8WxYpOoiYaQtkGHlkz90di0JvTHpKTND5JxQnSqJlRSiYm6ihqJsUkp0fojxQACRkUE7So7mLCiEiySKn77x3lWr+vd3bu79+5mHz6vmTv3nOc899zvuQOfffbZc85VRGBmZnk5aKALMDOz6nO4m5llyOFuZpYhh7uZWYYc7mZmGXK4m5llqKJwl3SEpBWSnpC0WdL7JB0paZWkren5namvJF0laZukRyRNr+0hmJlZR5WO3K8EfhQRk4FpwGZgIbA6IiYBq9M6wGxgUnosAK6tasVmZtYtdXcRk6TDgU3AsVHSWdIW4AMR8ZykscC9EXGcpG+l5Zs69qvZUZiZ2VtUMnI/FmgD/k3SBknXSRoOjGkP7PQ8OvUfB2wveX1rajMzs34ytMI+04FLI+IhSVfy5hRMOSrT9rZfDyQtoJi2Yfjw4b87efLkCkoxM7N269evfz4i6sttqyTcW4HWiHgora+gCPdfShpbMi2zs6T/+JLXNwA7Ou40IpYCSwGamppi3bp1FR2MmZkVJD3T2bZup2Ui4hfAdknHpaZZwOPASmBeapsH3JmWVwKfTGfNnATs9ny7mVn/qmTkDnAp8D1JBwNPAxdS/GC4RdJ84Fng/NT3buAsYBuwN/U1M7N+VFG4R8RGoKnMplll+gZwcR/rMjOzPqh05G5m9hvn1VdfpbW1lX379g10KTVVV1dHQ0MDw4YNq/g1DnczG7RaW1sZMWIEEyZMQCp3ot7gFxHs2rWL1tZWJk6cWPHrfG8ZMxu09u3bx6hRo7INdgBJjBo1qse/nTjczWxQyznY2/XmGB3uZmYDZNGiRVxxxRU12bfn3M0sGxMW/rCq+2tZ/EdV3V9/8sjdzKyPvvvd7zJjxgwaGxv51Kc+xdVXX81nP/vZN7Zff/31XHrppQB86Utf4rjjjuP0009ny5YtNavJ4W5m1gebN2+mubmZBx98kI0bNzJkyBAOO+wwbr/99jf6NDc3M2fOHNavX8/NN9/Mhg0buP3221m7dm3N6vK0jJlZH6xevZr169dz4oknAvDKK68wevRojj32WNasWcOkSZPYsmULM2fO5Morr+Tcc8/l0EMPBeDDH/5wzepyuJuZ9UFEMG/ePL7yla+8pX3ZsmXccsstTJ48mXPPPfeNM1766+weT8uYmfXBrFmzWLFiBTt3FjfGfeGFF3jmmWf4yEc+wve//31uuukm5syZA8App5zCHXfcwSuvvMKePXv4wQ9+ULO6PHI3M+uDKVOm8MUvfpEzzjiD119/nWHDhnH11VdzzDHHMGXKFB5//HFmzJgBwPTp05kzZw6NjY0cc8wxnHzyyTWrq9uv2esPvp+7mfXG5s2bOf744we6jH5R7lglrY+Icjd19LSMmVmOHO5mZhlyuJuZZcjhbmaWIYe7mVmGHO5mZhlyuJuZ9dJLL73ENddc02WflpYWbrzxxm731dLSwtSpU6tVmi9iMrOMLBpZ5f3t7nJze7hfdNFFnfZpD/dPfOIT1a2tGx65m5n10sKFC3nqqadobGzk8ssv5/LLL2fq1Km85z3vobm5+Y0+999/P42NjSxZsoSWlhZOPvlkpk+fzvTp0/npT39ak9o8cjcz66XFixfz6KOPsnHjRm677Ta++c1vsmnTJp5//nlOPPFETjnlFBYvXswVV1zBXXfdBcDevXtZtWoVdXV1bN26lblz51KLK/Qd7mZmVfDAAw8wd+5chgwZwpgxY3j/+9/P2rVrOfzww9/S79VXX+WSSy55497vTz75ZE3qcbibmVVBpffpWrJkCWPGjGHTpk28/vrr1NXV1aQez7mbmfXSiBEj2LNnD1Dczre5uZn9+/fT1tbGfffdx4wZM97SB2D37t2MHTuWgw46iBtuuIH9+/fXpDaP3M3MemnUqFHMnDmTqVOnMnv2bE444QSmTZuGJL72ta/xrne9i1GjRjF06FCmTZvGBRdcwEUXXcRHP/pRbr31Vk499VSGDx9ek9oquuWvpBZgD7AfeC0imiQdCTQDE4AW4GMR8aKKrxm5EjgL2AtcEBEPd7V/3/LXzHrDt/ytzi1/T42IxpIdLQRWR8QkYHVaB5gNTEqPBcC1PXgPMzOrgr7MuZ8NLE/Ly4FzStq/E4U1wBGSxvbhfczMrIcqDfcA/kPSekkLUtuYiHgOID2PTu3jgO0lr21NbWZm1k8q/YPqzIjYIWk0sErSE130LffV3m+b2E8/JBYAHH300RWWYWb2VhFB8ae+fPXm61ArGrlHxI70vBO4A5gB/LJ9uiU970zdW4HxJS9vAHaU2efSiGiKiKb6+voeF25mVldXx65du3oVfoNFRLBr164enw/f7chd0nDgoIjYk5bPAP4JWAnMAxan5zvTS1YCl0i6Gfg9YHf79I2ZWTU1NDTQ2tpKW1vbQJdSU3V1dTQ0NPToNZVMy4wB7ki/9gwFboyIH0laC9wiaT7wLHB+6n83xWmQ2yhOhbywRxWZmVVo2LBhTJw4caDL+I3UbbhHxNPAtDLtu4BZZdoDuLgq1ZmZWa/49gNmZhlyuJuZZcjhbmaWIYe7mVmGHO5mZhlyuJuZZcjhbmaWIYe7mVmGHO5mZhlyuJuZZcjhbmaWIYe7mVmGHO5mZhlyuJuZZcjhbmaWIYe7mVmGHO5mZhlyuJuZZcjhbmaWIYe7mVmGHO5mZhlyuJuZZcjhbmaWIYe7mVmGHO5mZhlyuJuZZcjhbmaWIYe7mVmGKg53SUMkbZB0V1qfKOkhSVslNUs6OLUfkta3pe0TalO6mZl1picj908Dm0vWvwosiYhJwIvA/NQ+H3gxIt4NLEn9zMysH1UU7pIagD8CrkvrAk4DVqQuy4Fz0vLZaZ20fVbqb2Zm/aTSkfs3gM8Cr6f1UcBLEfFaWm8FxqXlccB2gLR9d+pvZmb9pNtwl/RBYGdErC9tLtM1KthWut8FktZJWtfW1lZRsWZmVplKRu4zgQ9LagFuppiO+QZwhKShqU8DsCMttwLjAdL2kcALHXcaEUsjoikimurr6/t0EGZm9lbdhntE/G1ENETEBODjwE8i4o+Be4DzUrd5wJ1peWVaJ23/SUS8beRuZma105fz3D8HfEbSNoo59WWpfRkwKrV/BljYtxLNzKynhnbf5U0RcS9wb1p+GphRps8+4Pwq1GZmZr3kK1TNzDLkcDczy5DD3cwsQw53M7MM9egPqge8RSNrsM/d1d+nmR3wPHI3M8uQw93MLEMOdzOzDDnczcwy5HA3M8uQw93MLEMOdzOzDDnczcwy5HA3M8uQw93MLEMOdzOzDDnczcwy5HA3M8uQw93MLEMOdzOzDDnczcwy5HA3M8uQw93MLEMOdzOzDDnczcwy5HA3M8uQw93MLEPdhrukOkk/k7RJ0mOSvpDaJ0p6SNJWSc2SDk7th6T1bWn7hNoegpmZdVTJyP3XwGkRMQ1oBM6UdBLwVWBJREwCXgTmp/7zgRcj4t3AktTPzMz6UbfhHoWX0+qw9AjgNGBFal8OnJOWz07rpO2zJKlqFZuZWbcqmnOXNETSRmAnsAp4CngpIl5LXVqBcWl5HLAdIG3fDYyqZtFmZta1isI9IvZHRCPQAMwAji/XLT2XG6VHxwZJCyStk7Sura2t0nrNzKwCPTpbJiJeAu4FTgKOkDQ0bWoAdqTlVmA8QNo+EnihzL6WRkRTRDTV19f3rnozMyurkrNl6iUdkZbfAZwObAbuAc5L3eYBd6bllWmdtP0nEfG2kbuZmdXO0O67MBZYLmkIxQ+DWyLiLkmPAzdL+iKwAViW+i8DbpC0jWLE/vEa1G1mZl3oNtwj4hHgvWXan6aYf+/Yvg84vyrVmZlZr/gKVTOzDDnczcwy5HA3M8uQw93MLEMOdzOzDDnczcwy5HA3M8uQw93MLEMOdzOzDDnczcwy5HA3M8uQw93MLEMOdzOzDDnczcwy5HA3M8uQw93MLEMOdzOzDDnczcwy5HA3M8uQw93MLEMOdzOzDDnczcwy5HA3M8uQw93MLEMOdzOzDDnczcwy5HA3M8tQt+EuabykeyRtlvSYpE+n9iMlrZK0NT2/M7VL0lWStkl6RNL0Wh+EmZm9VSUj99eAv46I44GTgIslTQEWAqsjYhKwOq0DzAYmpccC4NqqV21mZl3qNtwj4rmIeDgt7wE2A+OAs4Hlqdty4Jy0fDbwnSisAY6QNLbqlZuZWad6NOcuaQLwXuAhYExEPAfFDwBgdOo2Dthe8rLW1GZmZv2k4nCXdBhwG3BZRPyqq65l2qLM/hZIWidpXVtbW6VlmJlZBSoKd0nDKIL9exFxe2r+Zft0S3remdpbgfElL28AdnTcZ0QsjYimiGiqr6/vbf1mZlZGJWfLCFgGbI6Ir5dsWgnMS8vzgDtL2j+Zzpo5CdjdPn1jZmb9Y2gFfWYCfwr8XNLG1PZ3wGLgFknzgWeB89O2u4GzgG3AXuDCqlZsZmbd6jbcI+IBys+jA8wq0z+Ai/tYlx0IFo2swT53V3+fZoOQr1A1M8uQw93MLEMOdzOzDDnczcwy5HA3M8uQw93MLEMOdzOzDDnczcwy5HA3M8uQw93MLEMOdzOzDDnczcwy5HA3M8uQw93MLEMOdzOzDDnczcwy5HA3M8uQw93MLEMOdzOzDDnczcwy5HA3M8uQw93MLEMOdzOzDDnczcwy5HA3M8uQw93MLEMOdzOzDHUb7pK+LWmnpEdL2o6UtErS1vT8ztQuSVdJ2ibpEUnTa1m8mZmVV8nI/XrgzA5tC4HVETEJWJ3WAWYDk9JjAXBtdco0M7Oe6DbcI+I+4IUOzWcDy9PycuCckvbvRGENcISksdUq1szMKtPbOfcxEfEcQHoendrHAdtL+rWmNjMz60fV/oOqyrRF2Y7SAknrJK1ra2urchlmZge23ob7L9unW9LzztTeCowv6dcA7Ci3g4hYGhFNEdFUX1/fyzLMzKyc3ob7SmBeWp4H3FnS/sl01sxJwO726RszM+s/Q7vrIOkm4APAUZJagX8EFgO3SJoPPAucn7rfDZwFbAP2AhfWoGYzM+tGt+EeEXM72TSrTN8ALu5rUWZm1je+QtXMLEMOdzOzDHU7LWMGMGHhD6u+z5a6qu/SzBKP3M3MMuRwNzPLkMPdzCxD2c65e47YzA5kHrmbmWXI4W5mliGHu5lZhhzuZmYZcribmWXI4W5mliGHu5lZhhzuZmYZcribmWXI4W5mliGHu5lZhhzuZmYZcribmWXI4W5mliGHu5lZhhzuZmYZcribmWXI4W5mliGHu5lZhhzuZmYZcribmWVoaC12KulM4EpgCHBdRCyuxfuYWbJoZA32ubv6+7R+U/WRu6QhwNXAbGAKMFfSlGq/j5mZda4WI/cZwLaIeBpA0s3A2cDjNXgvs0FnwsIfVn2fLXVV36UNcrUI93HA9pL1VuD3avA+ZmY1+mH5iarvs7+nuWoR7irTFm/rJC0AFqTVlyVtqUEtVSU4Cni+qjv9QrmP68Dgz7N6/FlW1yD6PI/pbEMtwr0VGF+y3gDs6NgpIpYCS2vw/jUjaV1ENA10Hbnw51k9/iyrK4fPsxanQq4FJkmaKOlg4OPAyhq8j5mZdaLqI/eIeE3SJcCPKU6F/HZEPFbt9zEzs87V5Dz3iLgbuLsW+x5gg2oaaRDw51k9/iyra9B/nop42986zcxskPPtB8zMMuRwNzPLUE3m3HMhaTLF1bXjKM7V3wGsjIjNA1qYHfDSv81xwEMR8XJJ+5kR8aOBq2xwkjQDiIhYm26XcibwRPr74aDkkXsnJH0OuJnioqyfUZziKeAmSQsHsrbcSLpwoGsYTCT9JXAncCnwqKSzSzZ/eWCqGrwk/SNwFXCtpK8A/wIcBiyU9PcDWlwf+A+qnZD0JPA7EfFqh/aDgcciYtLAVJYfSc9GxNEDXcdgIennwPsi4mVJE4AVwA0RcaWkDRHx3gEtcJBJn2cjcAjwC6AhIn4l6R0UvxmdMKAF9pKnZTr3OvBbwDMd2sembdYDkh7pbBMwpj9rycCQ9qmYiGiR9AFghaRjKH/7D+vaaxGxH9gr6amI+BVARLwiadD+X3e4d+4yYLWkrbx5I7SjgXcDlwxYVYPXGOAPgRc7tAv4af+XM6j9QlJjRGwESCP4DwLfBt4zsKUNSv8n6dCI2Av8bnujpJEM4oGcp2W6IOkgilsYj6MIoVZgbfopbz0gaRnwbxHxQJltN0ZEDW7DlydJDRSjzV+U2TYzIh4cgLIGLUmHRMSvy7QfBYyNiJ8PQFl95nA3M8uQz5YxM8uQw93MLEMOdzOzDDncrWYkLZL0N5L+SdLp1dhXtWr7TSPpMkmHdrH9On/RvPWET4W0mouIz/ekv6ShEfFa6Xr1q6odSaI4WaEnp9FdBnwX2Ftmf0Mi4s+qVZ8dGDxyt6qS9PeStkj6T+C41Ha9pPPS8uclrZX0qKSlKQiRdK+kL0v6L+DT6TVfl3QP8NW0+ymp39PpEvz29/wTST+TtFHStyQNSe3XSlon6TFJX+im7sWSHpf0iKQrUtsYSXdI2pQev5/aP5Pqf1TSZaltgqTNkq4BHgbGSzpD0n9LeljSrZIO6+S9/5Ligrl70vEi6eX0G89DwPvScTd1dVySzpL0hKQHJF0l6a7UPlzSt9PnvqHD7QosVxHhhx9VeVBcAPJz4FDgcGAb8DfA9cB5qc+RJf1vAD6Ulu8FrinZdj1wF8XVmACLKC52OoTiy4t3AcOA44EfAMNSv2uAT5a+F8U3gt0LnNBJ3UcCW3jz1OAj0nMzcFnJPkaWHONwivuPPAa8F5hAccHLSan/UcB9wPC0/jng8118di3AUSXrAXysZP1eoKmz4wLqKC62m5i23QTclZa/DPxJ+7EBT7bX5Ue+D4/crZpOBu6IiL1RXMJd7rtzT5X0ULqfx2nA75Rsa+7Q99Z46wVjP4yIX0fE88BOiqteZ1EE7lpJG9P6san/xyQ9DGxI79PZnPWvgH3AdZI+wptTI6cB1wJExP6I2A38QTrG/43iFgC3p+MGeCYi1qTlk9L7PZjqmkcX31Rfxn7gtk62lTuuycDTEfE/qc9NJf3PoLgJ1kaKHwZ1FFdbW8YG1VymDQqdXhUnqY5iZN0UEdslLaIImnb/2+ElHddLryLcT/HvV8DyiPjbDu81keK3hhMj4kVJ13d4rzcLLr73dwbFD4aPU9xe4rTODqOT9o71ClgVEXO76N+VfVHmSugujqurugR8NCK29LIWG4Q8crdqug84V9I7JI0APtRhe3u4Pp/mn8+rwnuuBs6TNBpA0pHpBlqHU4TtbkljgNmd7SDVMjKKe3dfRnGHwPZ9/0XqM0TS4ekYz5F0qKThwLnA/WV2uwaYKend6fWHSvrtLo5jDzCiguPt7LieAI5VcZdIgDklr/kxcGnJ3zd818gDgEfuVjUR8bCkZmAjxd007++w/SVJ/0oxZ91CcY/8vr7n45L+AfgPFfcCehW4OCLWSNpAMSf+NNDV/VZGAHem3ywE/FVq/zSwVNJ8it8U/iIi/juNln+W+lwXERtKQrW9rjZJF1Dc//+Q1PwPFPPd5SwF/l3ScxFxahfHu6nccUVxB8OLgB9Jer6kPoB/Br4BPJICvgX4YBefh2XA95Yxy4Skw6K4Q6SAq4GtEbFkoOuygeFpGbN8/Hn6o+ljFGf2fGuA67EB5JG7HVAk3QFM7ND8uYj48YHw/nbgcLibmWXI0zJmZhlyuJuZZcjhbmaWIYe7mVmGHO5mZhn6f2mQ34EsbyAeAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "df_diarrhea = df.loc[:, ['evd', 'diarrhea_score_triage']]\n",
    "df_diarrhea['total'] = 1\n",
    "df_diarrhea.groupby('diarrhea_score_triage').sum().plot(ylim=(0,600), kind='bar')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On second thought, there really are not enough samples for `diarrhea_score_triage > 1` to conduct a meaningful analysis (10 with a score of 2, and 2 with a score of 3). Therefore, we use the binary variable `diarrhea_triage`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "targets = ['evd']\n",
    "predictors = ['quarantine', 'quarantine_missing', 'conjonctivitis_triage', 'diarrhea_triage', 'referral(4-9)', 'pyrexia_triage', 'dysphagia_triage', 'hemorrhage_triage', 'myalgia_triage' ]\n",
    "assert set(predictors).issubset(df.columns)\n",
    "assert set(targets).issubset(df.columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "575\n"
     ]
    },
    {
     "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>quarantine</th>\n",
       "      <th>quarantine_missing</th>\n",
       "      <th>conjonctivitis_triage</th>\n",
       "      <th>diarrhea_triage</th>\n",
       "      <th>referral(4-9)</th>\n",
       "      <th>pyrexia_triage</th>\n",
       "      <th>dysphagia_triage</th>\n",
       "      <th>hemorrhage_triage</th>\n",
       "      <th>myalgia_triage</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Empty DataFrame\n",
       "Columns: [quarantine, quarantine_missing, conjonctivitis_triage, diarrhea_triage, referral(4-9), pyrexia_triage, dysphagia_triage, hemorrhage_triage, myalgia_triage]\n",
       "Index: []"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = df.loc[:, predictors].copy()\n",
    "print(len(X))\n",
    "X.head(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "575\n"
     ]
    },
    {
     "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>evd</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Empty DataFrame\n",
       "Columns: [evd]\n",
       "Index: []"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y = df.loc[:, targets].copy()\n",
    "print(len(y))\n",
    "y.head(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "575\n"
     ]
    },
    {
     "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>age</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Empty DataFrame\n",
       "Columns: [age]\n",
       "Index: []"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "age = df.loc[:, ['age']].copy()\n",
    "print(len(age))\n",
    "age.head(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save the feature matrix $\\mathbf{X}$, the target vector $\\mathbf{y}$, and the age (for reference)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "save_dir = '../data/private/predict_EVD'\n",
    "if not os.path.isdir(save_dir):\n",
    "    os.mkdir(save_dir)\n",
    "X.to_csv(os.path.join(save_dir, 'X.csv'), index=False)\n",
    "y.to_csv(os.path.join(save_dir, 'y.csv'), index=False)\n",
    "age.to_csv(os.path.join(save_dir, 'age.csv'), index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Train a logistic regression model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "def define_fit_eval_model(verbose=False):\n",
    "    X_train, X_test, y_train, y_test = sk_ms.train_test_split(X,y, shuffle=True, test_size=0.5)\n",
    "    model = sk_lm.LogisticRegressionCV(Cs=np.logspace(-3,1,num=60), cv=4, scoring='accuracy')\n",
    "\n",
    "    model.fit(X_train, y_train)\n",
    "    avg_scores = model.scores_[1].mean(axis=0)\n",
    "    if verbose:\n",
    "        print(\"Inverse of regularisation coefficient:\", model.C_[0])\n",
    "        print(\"Corresponding validation score (averaged across folds):\", avg_scores.max()*100, \"%\")\n",
    "\n",
    "        plt.semilogx(model.Cs_, avg_scores*100)\n",
    "        plt.title('Validation score (averaged across folds) during cross-validation')\n",
    "        plt.xlabel('C (inverse of regularization parameter)')\n",
    "        plt.ylabel('Accuracy [%]')\n",
    "\n",
    "\n",
    "    avg_scores = model.scores_[1].mean(axis=0)\n",
    "    if verbose:\n",
    "        print(\"Corresponding validation score (averaged across folds):\", avg_scores.max())\n",
    "\n",
    "    ## Evaluate the model\n",
    "\n",
    "    # We begin by evaluating its score (accuracy) on the test set:\n",
    "\n",
    "    score = model.score(X_test, y_test)\n",
    "    baseline = 1 - np.mean(y_test)[0]\n",
    "    if baseline < 0.5:\n",
    "        baseline = 1 - baseline\n",
    "    return score, baseline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/lib/python3.7/site-packages/sklearn/utils/validation.py:724: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().\n",
      "  y = column_or_1d(y, warn=True)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7291666666666667\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/lib/python3.7/site-packages/sklearn/utils/validation.py:724: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().\n",
      "  y = column_or_1d(y, warn=True)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7118055555555556\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/lib/python3.7/site-packages/sklearn/utils/validation.py:724: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().\n",
      "  y = column_or_1d(y, warn=True)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.75\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/lib/python3.7/site-packages/sklearn/utils/validation.py:724: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().\n",
      "  y = column_or_1d(y, warn=True)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7395833333333333\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/lib/python3.7/site-packages/sklearn/utils/validation.py:724: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().\n",
      "  y = column_or_1d(y, warn=True)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7256944444444444\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/lib/python3.7/site-packages/sklearn/utils/validation.py:724: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().\n",
      "  y = column_or_1d(y, warn=True)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7569444444444444\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/lib/python3.7/site-packages/sklearn/utils/validation.py:724: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().\n",
      "  y = column_or_1d(y, warn=True)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7083333333333333\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/lib/python3.7/site-packages/sklearn/utils/validation.py:724: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().\n",
      "  y = column_or_1d(y, warn=True)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7222222222222222\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/lib/python3.7/site-packages/sklearn/utils/validation.py:724: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().\n",
      "  y = column_or_1d(y, warn=True)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.6979166666666667\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/lib/python3.7/site-packages/sklearn/utils/validation.py:724: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().\n",
      "  y = column_or_1d(y, warn=True)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.6944444444444444\n"
     ]
    }
   ],
   "source": [
    "scores, baselines = [], []\n",
    "for _ in range(10):\n",
    "    score, baseline = define_fit_eval_model()\n",
    "    scores.append(100*score)\n",
    "    baselines.append(100*baseline)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfgAAAFOCAYAAABquw8XAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3zU9f3A8dc7gySsBEhYYUPCFhBcjDhQGaKgbRWr1k5Htda2vzo61La21WqXo1rb2mpd4KyKLEUZDiQMIcywAyQQRhhhZL1/f9z39IgZF3J33+9d3s/H4x53973vfb/vHCHv+37G+yOqijHGGGNiS5zbARhjjDEm9CzBG2OMMTHIErwxxhgTgyzBG2OMMTHIErwxxhgTgyzBG2OMMTHIErwxxjRRInKeiOyI8Dm/KSKLAp4fEZFekYyhqbAEb4xpEBHZKiIXuh1HtBCR+0TkObfj8CpVbamqm92OIxZZgjfGeJKIJLgdQyyxz7PpsQRvjDllTnPrhyLyZxEpEZHNIjLS2V4gIntE5PqA/f8jIk+KyFwROSwi80Wke8DrKiK3iEg+kO9sGykiS0TkoHM/0tk+VURyq8XzIxF503mcJCIPi8h2EdntnDfFee08EdkhInc4MRaKyBQRmSgiG0Rkv4j8LOC4cSJyl4hsEpF9IjJdRNo6r/Vw4r7eOddeEfm589p44GfAVU5T9Ge1fI79ReQD5zNcLSKXOdvPFpEiEYkP2PdyEVnZgLi+IyLbgXl1/Dv+zIl7q4hcE7D9EhFZLiKHnH/P+wJeSxaR55zzljj/Nh2c11JF5F/O57pTRO4P/BmqnVtFpE/A78fjIjLD+f1YLCK9A/bt5/zu7BeR9SJyZW0/k7EEb4xpvLOAlUA74AXgJeAMoA9wLfCYiLQM2P8a4DdAOrACeL7a8aY4xxzgJKsZwCPO8f8EzBCRdsCbQF8RyQp479edGAAeBLKBoU4smcA9Aft2BJIDtv/DiXc4MAa4J6Bv+DYnrnOBzsAB4PFqcY8G+gJjnff2V9VZwO+AaU5T9JDqH56IJAJvAXOA9sAPgOdFpK+qfgKUAhfU8jMGE9e5QH9gXPVzB3wO6c7ncD3wlIj0dV4rBb4BpAGXADeLyBTnteuBVKArvn+bm4BjzmvPABX4PvdhwMXAd2s5f3VXA78C2gAbgd8CiEgLYK7zs7d39vubiAwM8rhNj6razW52s1vQN2ArcKHz+JtAfsBrgwEFOgRs2wcMdR7/B3gp4LWWQCXQ1XmuwAUBr18HfFrt/B8D33QePwfc4zzOAg4DzQHBl5x6B7zvHGCL8/g8fMko3nneyjn3WQH7LwWmOI/XAmMDXusElAMJQA/nvV0CXv8UmOo8vg94ro7PcwxQBMQFbHsRuM95fD/wdECcpUD3BsTVq45zn4cvEbcI2DYd+GUt+/8F+LPz+NvAR8Bp1fbpAJwAUgK2XQ28H/A7syjgNQX6BPx+/DPgtYnAOufxVcDCauf6O3Cv2/8nvHqzPhljTGPtDnh8DEBVq28LvIIv8D9Q1SMish/f1WdB9ded7duqnW8bvqtN8F3N/RH4Nb4r2zdU9aiItMeX6JeKiP99AgQ2E+9T1crAuGv4WfxxdwdeF5GqgNcr8SUzv6KAx0er/cx16QwUqGrgsav/jB+JyM3AFcAyVfV/JsHEFfh51uSAqpZWO3dnABE5C3gAGAQ0A5KAl539/ovv6v0lEUnD92Xr505MiUBhwGcfF0QcfrV9jt2Bs0SkJOD1BCcOUwNrojfGRFpX/wOn6b4tsCvg9cAlLnfh+8MeqBuw03k8B0gXkaH4rhL9Tdd78SXogaqa5txSVTXYpFtdATAh4FhpqpqsqjvrfefJP09NdgFdRSTw7/HnP6OqrsGXdCdwcvN8sHHVd/42TvN34Ln9/x4v4OsK6aqqqcCT+L4ooarlqvorVR0AjAQm4WvOL8B3BZ8eEFNrVW1sU3oBML/az9pSVW9u5HFjliV4Y0ykTRSR0SLSDF9f/GJVre3q7h0gW0S+LiIJInIVMAB4G0BVK4BXgIfwfVGY62yvwten/mfnah4RyRSR2vqh6/Mk8FtxBgSKSIaITA7yvbuBHtUSeKDF+Jrd7xCRRBE5D7gU31gGvxfw9bfn8MUVdGPjCvQrEWkmImPwJWr/OVoB+1X1uIicie8LBs65zheRwc7guUP4ugYqVbUQ3xevP4pIa2cgYG8ROfcU4gr0Nr7fheuczylRRM4Qkf6NPG7MsgRvjIm0F4B7gf34BrRdU9uOqroPX8L5Cb6+/DuASaq6t9rxLgRedhK+3534Bml9IiKHgHfxDYI7FX/FdyU7R0QOA5/gGwgYDH+y3Cciy6q/qKplwGX4rtD3An8DvqGq6wJ2exFff/m8aj97Y+LyK8I3OG8XvgGPNwWc+/vAr51j34Ovf96vI74vV4fwjQWYj6+ZHnxX8s2ANc6xX8E3PuCUqephfIP1pjqxFuEbSJnUmOPGMnEGKhhjTNiJyH+AHar6C7djMSbW2RW8McYYE4MswRtjjDExyJrojTHGmBhkV/DGGGNMDLIEb4wxxsSgmKpkl56erj169HA7jCZt6dKle1U1w+04ws1+14wxXlDX39yYSvA9evQgNze3/h1N2IhI9bKiMcl+14wxXlDX31xrojfGGGNikCV4Y4wxJgZZgjfGGGNikCV4Y2JY4cFj9e9kjIlJluBNkyEiT4vIHhHJC9j2kIisE5GVIvK6s651TFiYX8w5v5/H0m373Q7FGOMCS/Ae9cbynYx6YB4975rBqAfm8cbyYJadNvX4DzC+2ra5wCBVPQ3YANwd6aDC5c0VviW9567Z43Ikxhg3WIL3oDeW7+Tu11axs+QYCuwsOcbdr62yJN9IqroA3xKlgdvmBCwx+gnQJeKBhUF5ZRVz1+4GfFfyxpimxxK8Bz00ez3HyitP2nasvJKHZq93KaIm49vATLeDCIXFm/dTcrScIV3TWL3rEMWHT7gdkjEmwizBe9CukpoHRtW23TSeiPwcqACer2OfG0QkV0Ryi4u9fVU8M6+QlMR4fnFJfwAWbfR2vMaY0LME70Gd01IatN00johcD0wCrtE6lldU1adUdYSqjsjI8G413soqZfbq3ZzfL4Ph3drQrkUzFmzY63ZYxpgIswTvQT8d15eUxPiTtqUkxvPTcX1diih2ich44E7gMlU96nY8obBs+wH2HjnBuIEdiYsTRmelszC/mKoqWxramKbEErwHTRmWye+vGExmWgoCZKal8PsrBjNlWKbboUU1EXkR+BjoKyI7ROQ7wGNAK2CuiKwQkSddDTIEZq4qoll8HBf0aw9ATlYGe4+UsabwkMuRmWh15EQFD8xcx4HSMrdDMQ0QU4vNxJIpwzItoYeYql5dw+Z/RTyQMFJVZq8uYkxWOq2SEwEYk50OwIL8YgZlproZnolS05YU8OT8TRw5Uc79Uwa7HY4Jkl3BGxNDVu08yM6SY4wf1PHzbe1bJdO/U2sWbLCBdqbhVJXpSwoAePHTAjbuOexyRCZYluCNiSEz84qIjxMu7N/hpO05Weks3XaA0hMVtbzTmJp9tuMg63cf5v8uzqZ5YjwPzFzndkgmSJbgsapxJjaoKrPyijinVzvatGh20ms52RmUVyofb9rnUnQmWk1bsp2UxHiuH9mD75/fh3fX7rHfoyjR5BO8VY0zsWL97sNs2Vt6UvO834gebUhJjGeBVbUzDXC0rIK3Pitk4uBOtEpO5FujepCZlsLv3llrszKiQJNP8FY1zsSKWXlFiMDFAzt86bWkhHjO7tXW+uFNg8xYWciRExVcdUZXAJKd6bqrdh7kzc92uRydqU+TT/BWNc7Eill5RYzo3ob2rZJrfD0nO4Ot+46yfV9MTPc3ETA9t4Be6S04o0ebz7ddNqQzgzJb89Ds9RyvdnFkvKXJJ3irGmdiwZa9pawrOsz4QZ1q3Scn21d9z5rpTTA27jnCkq0HuPKMrojI59vj4oSfTezPzpJj/PvDre4FaOrV5BO8G1XjbFCfCbWZeYUANfa/+/VKb0FmWoo105ugvJxbQHyccMXpX67HMbJ3Ohf2b8/f3t/IviO2kJFXNfkEH+mqcTaoz4TD7LwihnRJJbOOlicRISc7g4827aO8siqC0ZloU15ZxavLdnBBv/a1dvncNaEfR8sreeS9/AhHZ4JlleyIbNW4ugb1WeU6cyp2lhzjsx0HuWN8/a1O52an8+Kn21m+vYQze7aNQHQmGr23dg97j5Qx1RlcV5M+7Vtx9ZldeX7xdq4f2YNeGS0jGKEJRpO/go80G9RnQm1WXhEAE+rof/cb2Sed+DixZnpTp+m5BbRvlcS52XWvmnj7hdkkJ8bz4CwrfuNFluAjzAb1mVCblVdIv46t6Jneot59WycnMqxrmg20M7UqOnicD9bv4avDu5AQX3eKSG+ZxM3n9Wb26t18umV/hCI0wbIEH2GhHNRng/XMnsPHyd12oM7BddXlZGewaudB9tvKYKYGry7bQZXClSNqb54P9O1RPenYOpnfzlhjxW88JqwJXkR+JCKrRSRPRF4UkWQReV5E1jvbnhaRxFreW+ks37lCRN4MZ5yRFKpBfTZYzwDMWb0b1bpHz1eXk52BKiy0q3hTTVWVMm1JAWf3akuPIFqEAFKaxfN/4/ry2Y6DvL2qMMwRht/+0jLeWL6TihgYiBq2BC8imcBtwAhVHQTEA1OB54F+wGAgBfhuLYc4pqpDndtl4YrTDVOGZfLhXRew5YFL+PCuC05pcJ1V4DPg63/vmd6Cvh1aBf2ewZmppDVPZMGGvWGMzESjT7bsY/v+o59XrgvW5cMyGdCpNQ/OXBfVxW/KK6u44dlcbp+2gmv+uZg9h4+7HVKjhLuJPgFIEZEEoDmwS1XfUQfwKdAlzDHEJBusZw6UlvHx5n2MH9TxpEIk9YmPE0b1SWdhfjG+/4bG+ExfUkCr5ISgBmwGio8Tfn6Jr/jNsx9vDUtskfC7d9aSu+0A157djc92lDDpkUUs2Rq9YwvCluBVdSfwMLAdKAQOquoc/+tO0/x1wKxaDpEsIrki8omITAlXnNHKBuuZd9fuprJKmdCA5nm/c7My2HP4BOuKbG1v43PwaDnv5BUxZWgmydXGCQVjVJ90zu+bwWPzNnIgCsd3vPXZLv794Va+ObIH908ZzBu3jKJ5s3imPvUJ/1y4OSq/DIezib4NMBnoCXQGWojItQG7/A1YoKoLazlEN1UdAXwd+IuI9K7lPDc4XwRyi4vD16fotQFtblTgM94yK6+IzLQUBmemNvi9Y7LTAWy6nPnc/z7bSVlFVYOb5wPdPbE/R05U8Oi8jSGMLPzydx/mzldXMrx7G342sT8A/Tq25s0fjObC/u25f8ZabnlhGYePl7scacOEs4n+QmCLqharajnwGjASQETuBTKAH9f2ZlXd5dxvBj4AhtWy31OqOkJVR2Rk1D1n81R5cUBbpCvwGW85fLychfl7GTewYc3zfp1SU8ju0NKmy5nPTVtSwIBOrRl0Cl8Y/bI7tOKqM7rx30+2snVvaQijC58jJyq46bmlNG8Wz+NfP51mCV+kxdbJiTx57XB+NrEfs1fvZvLjH7Jhd/S0eoUzwW8HzhaR5uL7CzQWWCsi3wXGAVerao3DFEWkjYgkOY/TgVHAmjDGWievDmgLxWA9E53mrdtDWWUVEwY3vHneLycrgyVbDnCsLHoHRZnQyNt5kNW7DjH1zFO/evf70UVZJMbH8YfZ3i9+o6rc+cpKtuwt5ZGrh9Ex9ctleUWEG3J68/x3z+LQsQomP/Yh/1sRHbOVwtkHvxh4BVgGrHLO9RTwJNAB+NiZAncPgIiMEJF/Om/vD+SKyGfA+8ADqupagrcBbcZrZq8uIqNVEsO7tal/51rkZGdQVlnFJ1v2hTAyE42mLSmgWUIck4c0/iKhfatkbszpzTurili6zdsD1P61aAszVhVyx/h+jOydXue+Z/dqxzu3jWZQZmt++NIK7v1fHmUV3p5KF9ZR9Kp6r6r2U9VBqnqdqp5Q1QRV7R0wBe7Xzr65qvpd5/FHqjpYVYc49/8KZ5z1sQFtxkuOlVXy/rpiLh7Qgbi4hjfP+53Zsy1JCXHWD9/EHS+v5I0VO5kwqCOpzWssS9Jg38vpSftWSdw/Y61nB6d9umU/v5+5josHdODGnF5Bvad962Re+N7ZfG9MT575eBtX/v1jT1/oWSW7INiANuMl8zcUc6y8ssFTmapLToznrF7tLME3cbPyijh8vKJRg+uqa94sgf+7uC/Lt5fwzqqikB03VPYcPs6tLyyjW9vmPHzlkAaNY0mMj+Pnlwzgb9eczsY9R5j06CIW5XuzpoQl+CDYgDbjJbPyCklrnshZvRq/GlxOVjqbikvZ6eGrEBNeLy3ZTre2zTm7Z7uQHvcrw7vQr2MrHpy1zlNN2RWVVdz6wnIOHS/niWtPp3XyqbVaTBzcif/dOor0ls247unFPDYv33Olei3BB8kGtBkvKKuo4r21e7iofwcS61kIJBj+1cLsKr5p2ravlE827+fKEV0a1d1Tk/g44WcT+7N9/1H++8m2kB67Mf4we72vef6KwfTr2LpRx+qd0ZI3bhnFZUM68/CcDXz32VwOHvXOVDpL8MZEkQ837eXwiYoG1Z6vS5/2LemUmmwJvomanltAnMBXh4eueT5QTnYGOdkZPPJevicS38xVhTy1YDPXnd2dy4eFpohq82YJ/OWqofxm8kAW5hdzyaMLydt5MCTHbixL8DHOawV6TOPMWlVEy6QERmfVPeI3WCJCTlYGizbujYnFNUzwKiqreDl3B+f1bV/j9LBQuXtCPw4dL+ex9/PDdo5gbCo+wk9fWcmQrmn8YlL/kB5bRLjunB5Mv/EcKquUK574iGlLtof0HKfCEnwM82KBHnPqKiqrmLOmiAv6tScpoeGlRGszJjudw8cr+GxHSciOabxv/oZi9hw+EfSysKeqf6fWfG14F575aBvb9x0N67lqc7SsgpufW0qzhDieuOb0kP7/CTSsWxve/sFozurZljtfXcVPX/7M1cV3LMHHMK8W6DGn5tOt+zlwtDxkzfN+o/ukEycw31aXa7C1hYeidvW0aUsKSG/ZjLH924f9XD++qC/xceJK8RtV5a5XV5G/5wiPTB0W9unN7Vom8Z9vncltF/Th5aU7uOJvH7n2xcYSfAyzAj2xZVZeEcmJcZzXN7QlmdOaN+O0LmnWD99AhQePMenRRdz92iq3Q2mwPYePM2/dHr5yepeQDNasT8fUZL6X04u3VxayfPuBsJ8v0LMfb+PNz3bxk4uyQ9a1VZ/4OOHHF/fl3988g50lx5j06ELeXbM7IucOZAk+hlmBnthRVaXMyivi3OwMmjdLCPnxc7IzWLmjhJKj0bcKmFtm5xVRWaW8vnwnH6zf43Y4DfLasp1UVClfC3PzfKAbc3qR3jKJ370TueI3S7cd4P4Zaxjbrz3fP69PRM4Z6Px+7Xn7B6Pp3q4F3302l4dmr6MyglPpLMHHMCvQEzuWFxxgz+ETjS5uU5tzs9OpUli00ZrpgzUzr4heGS3o074lP389jyMnKtwOKSiqyvQlBYzo3oY+7VtG7LwtkhL4ycXZLNl6gNmrw381u/fICW55fhmdUlP405VDQz4NMFhd2zbn5ZvO4eozu/H4+5v4xtOL2XvkRETObQk+hlmBntgxK6+IxHjh/H7h6S8d0iWNVskJLLR++KDsPXKCJVv3M2lwJx78ymB2HTzGQ7O8v7gKQO62A2zeWxrSynXB+trwLmR3aMkDM9eGtfhNRWUVt724nANHy3ji2tNDVoL3VCUnxvP7Kwbz0FdPI3frASY9soil28LfVWEJPsZZgZ7op6rMzCtiVJ90UlPC84cqIT6O0X3SWZBf7Nna4V4yd81uqhTGD+rE8O5tuf6cHjz7yTZyt3p7cRWAlz4toGVSApecFp7WoLokxMdx94T+bN13lBcWh6/4zR/nbuCjTfu4f8ogBnY+9eVvQ+1rI7ry2vdHkpQYx1V//5j/fLglrP/fLMEb43Grdx1ix4FjTAjx6PnqcrIzKDx4nI17joT1PLFgZl4R3ds1p3+nVoCvO6xzagp3vrrS06PqDx8v551VhVw6pFNYxnIE47y+GYzq046/vpfPwWOhL34zZ3URT3ywiavP7BrRMQbBGtg5lTdvHc15fdtz31truO2lFZSGqXvHErxpMkTkaRHZIyJ5AdvaishcEcl37k99/dUwmZVXRJzARQPCn+DBNz/a1O7g0XI+2riX8YM6fr5ISYukBH53xWA2FZfy+PsbXY6wdm99Vsix8sqwz32vi4ivhG3JsXL+9kFoP6ute0v5yfTPGJyZyr2XDgzpsUMpNSWRp64bzp3j+zFj5S4mP/5hWL5Yx3yCt0pu9WtCn9F/gPHVtt0FvKeqWcB7znNPmZlXyFk929G2RbOwniczLYXeGS1Y4NGVsbzi3bW7qahSxg88+QvXudkZXHF6Jk98sIm1hYdciq5u05Zsp2+HVgztmuZqHAM7p3LFsC78+8Ot7DgQmjnix8oquem5pcTHC3+75nSSE8NTzCZU4uKEm8/rzXPfPYuSo2VMfmwRM1cVhvYcIT2ax1glt/o1pc9IVRcA1TtJJwPPOI+fAaZENKh65O8+zKbiUiYMDu/Vu19OdgaLN+/zdDOz22atLqJTajJDunw5Sf7ykgGkpiRy56srPVf6d13RIT7bcZArz+jaoOVRw+X/xmUjwMMhKLylqvz8jVWs332Yv1w1lK5tmzc+wAgZ2Tudt38whv6dWhPq3viYTvBWya1+9hnRQVULAZz78Jf1aoCZeb61tMcNjFyCP1FRxadbvD9YzA2lJypYsKGYcQM71jjtqk2LZvxq8kBW7jjI0x9ucSHC2k1bUkBivHC5RwbadkpN4XtjevHGil2sbGSZ5OcXb+e1ZTv54dgszuvrqf/CQemYmsz0G89h4uDQDnyM6QRvldzqZ59R8ETkBhHJFZHc4uLI9FPPyitiePc2dGgdvsVAAp3dsx3N4uOsql0t3l+/hxMVVXUOeLxkcCcu7N+BP83dwNa9pRGMrnYnKip5fflOLh7YMexdPQ1x03m9SW/ZjN/OOPXiNysKSvj1W2s4NzuD2y7ICnGEkROOefoxneCtklv97DNit4h0AnDuay1JpqpPqeoIVR2RkRHacrE12b7vKGsKD32przecUprFc0bPNizItwRfk5l5RaS3bMaIHm1r3UdEuH/KIBLj4rj7tVWemHY4Z/VuSo6Wc5XHRpW3TErghxdms3jLft5d2/BqgPtLy/j+c0vJaJXEX65yr5iNV8V0grdKbvWzz4g3geudx9cD/3MxlpPMzPMNuAn14jL1ycnKYMPuIxQetFacQMfLK3l/3R4uGtCR+HoSScfUZO6e2J+PN+9j2pKCCEVYu+m5BWSmpTC6T2RqsTfE1DO60jujBb+fuZbyBoxbqKxSfvjScvYe8RWzaeOhlgmviOkEb5Xc6teUPiMReRH4GOgrIjtE5DvAA8BFIpIPXOQ894SZeUUMymwd8QFD/ulyVtXuZAvz93K0rDLoegRTz+jK2b3a8tt31rL70PEwR1e7gv1HWbRxL18b0cWTV7iJTvGbzcWlvPRp8Guo//XdDSzM38uvJg/ktBoGPBpwp9JBBE0ZlhmTySqUmspnpKpX1/LS2IgGEoTCg8dYUVDiSktKv46taN8qifn5xVzpQjlTr5qZV0jr5ATO6d0uqP3j4oQHrjiNcX9ZwC/eyOOp64a7Mnr95aU7ADxZ9MVvbP/2nN2rLX95N58pwzJplVx3xcZ563bzyLyNfG14F6ba72itYvoK3phoNTvCo+cDiQhjsjJYlL83oitfeVlZRRXvrtnNhQM6NGh51R7pLfjxRdnMXbObd1YVhTHCmlVWKa/kFjC6TzqZHh5XIyL8fOIA9pWW8eT8TXXuW7D/KLe/tIIBnVrzmymDPDHlz6sswZuoISKnB3Eb7HacoTAzr4is9i0jutpXoJzsdA4eK2fVzoOunN9rPtm8j0PHK05pNb/vjO7pq6z2Zl7El+NdtHEvuw4ed2VhmYYa3CWVy4dl8s+FW2qdxXO83FfMBuDJa4d7vpiN22K+iT6S3li+k4dmr2dXyTE6p6Xw03F9m0TTdwTNB5YAdX1l7wn0iEg0YeJfqezW8yO/frXfmKwMRGDBhmLXq555wcy8Ilo0i2dMVsMHqSXEx/HgV07jsscW8Zu31/LHK4eEIcKaTV9SQJvmiVw0oEPEztkY/zeuLzNWFfLwnPX86cqhX3r93v+tZvWuQ/zr+hF0axc9xWzcYlfwIdKUKsK5aImqXqCq59d2Aza7HWRjBa5U5pa2LZoxODPV5sPja+aeu6aI8/u1P+UrxgGdW3Pjub14ddmOiH2m+46cYM6aIi4f1oWkhOi40s1MS+Hbo3ry+vKd5FVrPZq2ZDvTcgu49fw+jO0fHV9Y3BbWBC8iPxKR1SKSJyIvikiyiPQUkcXO4h7TRKTGuQ0icreIbBSR9SIyLpxxhoJVhAs/Vb0gFPt43cy8Irq1/WKlMrfkZGWwvKCEQ8dDv+JXNFmydT97j5Q1erriDy7IoldGC+5+bVXYVg8L9PrynZRXalQ0zwf6/vm9SUtJPKn4zaodB/nl/1Yzuk86P7oo2+UIo0fYEryIZAK3ASNUdRAQD0wFHgT+7CzucQD4Tg3vHeDsOxDf4iB/ExFPfwW1inCRIyKjRKSF8/haEfmTiHR3O65Q8K9UNiFgpTK35GRnUFmlfLSxaU+Xm5VXRFJCHOc3sgRqcmI8D37lNHaWHOPhOeH94q+qTM8tYGjXNPp2dPeLYkO1Tk7k9guz+XjzPt5fv4eSo2Xc/PxS0ls0469Th9Zbg8B8IdxN9AlAiogkAM2BQuAC4BXn9doW95gMvKSqJ1R1C7ARODPMsTaKVYSLqCeAoyIyBLgD2AY8625IofH5SmURLm5Tk2Hd0miZlMD8JjwfvqpKmb26iJzsDFokNX7I0hk92vKNc7rzn4+2snTbgRBEWLPlBSVs2H0k6q7e/b5+Vjd6prfgd++s40fTVrD70HEev+Z02qrgSwgAACAASURBVLVMcju0qBK2BK+qO4GHge34EvtBYClQoqr+9qkdQE2j0DKBwPJPte3nSn3wmlhFuIiqUF/b3WTgr6r6VyC6LlNqUddKZZGWGB/HyN7tWLCh2BPlVt3w2Y4SCg8eD7q4TTDuGN+PTq2TuevVlZyoCM+qfdOXFJCSGM+k09wbx9EYifFx3DWhHxv3HOH99cXcM2kAw7q1cTusqBPOJvo2+P4A9wQ6Ay2ACTXsWtNfjpraYGr8CxPp+uC1aUoV4TzgsIjcDVwLzHC6b+qujBEF6lupzA052RnsLDnGZo8smhJps/KKSIgTxvYL3aCulkkJ/PbyweTvOcLj79c95/tUlJ6o4K3PdnHJaZ3qLRjjZRcP6MBlQzrzrVE9uPbsmOiBi7hwTpO7ENiiqsUAIvIaMBJIE5EE5yq+C7CrhvfuAALblmrbz1OaSkU4D7gK+DrwHVUtEpFuwEMux9Ro/pXKvNA875eT5fvSvGBDMb0z3JmT7xZVZdbqIkb2SSe1eWgT5fn92jNlaGee+GAjEwd3pF/H1iE79oyVhZSWVUZ9hTcR4ZGrh7kdRlQLZx/8duBsEWkuvtFCY4E1wPvAV519alvc401gqogkiUhPIAv4NIyxmujyI1X9k6ouBFDV7fgGZEY1/0plZ9SxUlmkdWvXnB7tmjfJ6XJrCw+zbd/RkDbPB7rn0oG0Sk7kzldXhbRi4LTcAnpltGB4d2vSburC2Qe/GN9gumXAKudcTwF3Aj8WkY1AO+BfACJymYj82nnvamA6vi8Es4BbVDU8nVUmKG8s38moB+bR864ZjHpgntvz+y+qYVtN3T9RoyErlUVaTnYGn2zeH7b+Yq+alVdInBC2IjFtWzTj3ksH8FlBCf/+cEtIjrlxz2GWbjvAVSO6uj4Lw7gvrKPoVfVeVe2nqoNU9TpnVPxmVT1TVfuo6tdU9YSz75uqek/Ae3+rqr1Vta+qzgxnnKZuXiniIyI3i8gqfKvBrQy4bcH3JTJq+Vcq81LzvF9OVgbHyivJ3Rq+Ud9eNDOviDN6tCU9jCO3LxvSmbH92vPwnPVs33e00cebnruDhDjhitO7hCA6E+2skp2pl4eK+LwAXIqvC+fSgNtwVb0m0sGE0ucrlfUKbqWySDqndzsS46VJNdNv3HOE/D1HwtY87yci3H/5IBLi4rj79ZWNmq1QVlHFq0t3MLZ/ezJa2XQyYwneBMErRXxU9aCqbnWWfd0BlOObXdHSGWgXlQJXKmuW4L3/ki2SEhjevQ0L8pvOfPjZq30rv0WiXHCn1BTumtCPDzfu4+XcHad8nHnrdrOvtCxq576b0PPeXxPjOV4r4iMitwK7gbnADOf2tivBhEBjViqLlJzsDNYWHmLP4eNuhxIRM/MKGdYtjY6pyRE539fP7MaZPdty/4w17Dl0ap/xtCUFdGyd/PnMB2MswZt6ebCIz+1AX1UdqKqDndtpbgXTWDPzimh+iiuVRYo/aSxsAlXtCvYfJW/nIcYPjNx4iLg44YErBnO8oop7/re6we8vPHiM+RuK+erwLiQ0YL16E9vsN8HUy4NFfArwVUaMeqFYqSwSBnRqTXrLZizIj/1+eH/zfKRbVHpltOT2C7OYtbqImasKG/TeV3J3UKVw5QhrnjdfsPXgTVA8VsRnM/CBiMwATvg3quqf3Avp1PhXKgv3YK7GiosTxmRlMH9DMVVV6plKe+EwM6+IAZ1au7Le+PfG9GLGykLueXM1I3sHV2CnqkqZvrSAc3q1szXSzUnsCt5Eo+34+t+b4atB779FnVl5RTRLiOO8Rq5UFgk52ensLy1j9a5DbocSNrsPHWfptgOufeFKjI/jwa+cxv7SMn77zpqg3vPJ5n0U7D/G1DPt6t2czK7gTdRR1V+5HUMofL5SWVYGLUOwUlm4jfGXrc0vZnCXVJejCY85n4+ed69FZVBmKjfk9OKJDzZx2ZBMRtczNuOlJQW0Tk5gXATHDJjoYFfwJmqIyF+c+7dE5M3qN7fja6hwrFQWTuktkxjYuTXzY3g+/My8InpntCCrg7sNQj8cm0XP9Bbc/fpKjpZV1LrfwaPlzFpdxJRhmZ4ew2HcYQneRJP/OvcPA3+s4RZV/CuVXdg/PKVQwyEnO4Nl2w5w+Hi526GE3P7SMhZv2e+J6YrJifE8cMVgCvYf449zNtS63xsrdlJWUWVz302NLMGbqKGqS537+cDHwD7n9pGzLWr4Vyo7p3e7kK9UFk5jstKpqFI+3rTP7VBCbu6aIiqr1DPlgs/q1Y5rzurGvz/cwvLtXy4TrKq8tKSAQZmtGdg5NrtMTONYgjdRR0TOA/KBx4G/ARtEJMfVoBroi5XK3L9abIgR3dvSvFl8TE6Xm5VXRJc2KQzsHLqlWxvrrgn96NA6mbteXUVZRdVJr+XtPMTawkNcZVPjTC0swZto9EfgYlU9V1VzgHHAn12OqUH8K5VdPDB6mucBmiXEcU6vdiyIsYI3h46Xs2jjXiYM6uipVdhaJSdy/5RBrN99mCc+2HTSa9Nyt5OUEMdlQz0zfdV4jCV4E40SVfXzlW5UdQMQPe3cRGalsnDJyc5g+/6jbN1b6nYoITNv7R7KKzUitecbamz/Dlw2pDOPvZ/Pht2HAThWVsn/lu9i4uBOpKZE1a++iSBL8CYa5YrIv0TkPOf2D2Cp20EFy79SmVf6ehsqJ/uL6XKxYmZeIR1aJzGsa5rbodTo3ksH0DIpgTtfXUlllTIzr5DDJyqscp2pkyV4E41uBlYDtwE/BNYAN7kaUQPM9sBc68bo0a45XdumxEwz/dGyCuZvKGbcwI6erdDXrmUS9146kOXbS3jmo61MW1JA93bNObtXW7dDMx7m/eoaxlSjqidE5DHgPaAKWK+qZS6HFbSZeYUM7ZpGp1R3VuNrLBEhJyuDN5b7pmh5cYnbhpi/vpjj5VWe/8I1eWhn3lixkz/MXsfx8ip+Oq6vp8YLGO+J7v+ZpkkSkUuATcBfgceAjSIywd2oguNfqSxaitvUJic7g9KySpbVMH0r2szMK6Jti2ac2cPbV8Miwm8vH0y8CHECXx3exe2QjMfZFbyJRn8EzlfVjQAi0hvfmvAzXY0qCNHePO83snc7EuKEBRuKObtXO7fDOWUnKiqZt24Pk07rFBXLrGampfCnq4ZSsP8oHVpHZq16E728/xttzJft8Sd3x2Zgj1vBNMTMvCL6d2pN93Yt3A6lUVolJ3J6tzZRP9Duw417OXKignFR9IVr3MCOfHdML7fDMFHAEryJRqtF5B0R+aaIXA+8BSwRkStE5Aq3g6uN2yuVhVpOdjp5Ow+x98iJ+nf2qJmrimiVnMCo3nUv6GJMNLIEb6JRMrAbOBc4DygG2gKXApNO5YAi8iMRWS0ieSLyooiEvP3T3zwfOwneN11uUX50jqYvr6xi7trdXNi/Q9QPFDSmJtYHb6KOqn4rlMcTkUx8U+4GqOoxEZkOTAX+E8rzzMoroldGC/q0bxnKw7pmUOdU2jRPZMGGYqYMi75qaos376fkaHnUj4cwpjb2tdUYnwQgRUQSgObArlAe/IuVyrxVCrUx4uKE0VkZLMjfS1WVuh1Og81aXUhKYjw5zjr3xsQaS/CmyVPVnfiWoN0OFAIHVXVOKM/hX6ks2haXqU9OVjp7j5xgbdEht0NpkKoqZfbq3ZzfL4OUZraOuolNYUvwItJXRFYE3A6JyO0iMi1g21YRWVHL+7eKyCpnv9xwxWmMiLQBJgM9gc5ACxG5tob9bhCRXBHJLS5u2OhxL65UFgqfl62Nsqp2S7cfoPjwCU/WnjcmVMKW4FV1vaoOVdWhwHDgKPC6ql4VsP1V4LU6DnO+s++IcMVpopeInC0i80TkQxGZ0ohDXQhsUdViVS3H9zs5svpOqvqUqo5Q1REZGQ1r1v3NlEH88WtDYqZ53q9D62T6dWzFgg3RNV1u5qoimsXHcUG/9m6HYkzYRKqJfiywSVW3+TeI7y/dlcCLEYrBRDkRqT4a6sfAZcB44DeNOPR24GwRae78Xo4F1jbieF/SpU1zzorigjB1ycnOIHfbfkpPVLgdSlBUldmrixiTlU7LJBtnbGJXpBL8VL6cyMcAu1U1v5b3KDBHRJaKyA1hjc5EiydF5JcBU9hKgK8DVwGn3AmsqouBV4BlwCp8/y+eamSsTUZOVgbllcriLfvcDiUoq3YeZGfJMRs9b2Je2BO8iDTDd5X1crWXrqbuq/dRqno6MAG4RURyajn+KfeLmuiiqlOAFcDbInIdcDu+xWaaA41pokdV71XVfqo6SFWvU9Xord4SYSN6tCE5MS5q+uFn5hWRECdcNKCD26EYE1aRuIKfACxT1d3+Dc5UpCuAabW9SVV3Ofd7gNeBM2vZ75T7RU30UdW3gHFAGr6+8vWq+oiq2rc7lyQnxnN2r3ZR0Q+vqszKK+Kc3u1Ia97M7XCMCatIJPiartQvBNap6o6a3iAiLUSklf8xcDGQF9YojeeJyGUisgiYh+/3YSpwuVN5rre70TVtOVkZbN5bSsH+o26HUqcNu4+wZW8p4wZa87yJfbWOMBGRHwfx/lJV/Xsdx2gOXATcWO2lL/XJi0hn4J+qOhHoALzujDhOAF5Q1VlBxGNi2/3AOUAK8I6qngn8WESygN/i+70yLvh8ulx+Mdec1d3laGo3M68QEbh4oDXPm9hX1xDSnwJPAHXN67kJqDXBq+pR4EtDh1X1mzVs2wVMdB5vBobUcV7TNB3El8RTCFg9zhmoacndRb0zWpCZlsKCDd5O8LPyijije1vat7KlVk3sqyvB/1dVf13Xm53mc2Mi5XJ8XT7l+EbPG48QEXKy03n7s0LKK6tI9ODa6lv2lrKu6DD3TBrgdijGRESt/wtV9Y763hzMPsaEiqruVdVHVfVJVY2u2qhNQE5WBodPVLCioMTtUGo0K8+3ml80rf1uTGME/TW7WtWwy8MZlDE1EZFlodjHhMfI3unECZ4dTT8rr5AhXVLJTEtxOxRjIqKuQXYdVbUoYJO/apgAH+GbumZMJPUXkZV1vC5AaqSCMSdLbZ7I0K5pLNhQzE8u7ut2OCfZWXKMz3Yc5M7x/dwOxZiIqasP/kkRWQo8pKrH+aJqWBWNqBpmTCME89e5MuxRmFrlZGfw1/fy2V9aRtsW3pln7m+et+p1pimpqw8+bFXDjDkVqrotiFuNtRVMZORkZ6AKizZ6q6rd7Lwi+nVsRc90Gxdsmo46++CtapgxpiGGdEkjrXkiT36wiV0lx9wOB4A9h4+zZNt+u3o3TU6tCd6qhhljGio+Tnjoq0PYtq+USY8uYlG++1fyc1bvRhUm2Nrvpomp6wr+fnxX718BHlTVElX9MXAPvqphxrhGRLqLyIXO4xR/aWPjvosGdODNH4ymXYtmXPf0Yh6bl09VlboWz6y8InqmtyC7Q0vXYjDGDXUleH/VsKlUqxqmqlY1zLhGRL6Hb3lXfxXFLsAb7kVkquud0ZI3bhnFpad15uE5G/jus7kcPFoe8ThKjpbx8eZ9jB/UEaf0tTFNRl0J/nJ8A+oqsKphxltuAUbhzOZwStW2dzUi8yUtkhL469Sh/HryQBbmFzPpsYXk7TwY0RjmrtlNZZUywfrfTRNU1yh6qxpmvOqEqpb5nzjLD7vXBmxqJSJ845weTLvxHCoqlSue+IjpSwoidv5ZeUVkpqUwONPKI5imp65BdlY1zHjVfBH5GZAiIhcBLwNvuRyTqcPp3drw9g9Gc2aPttzx6krufGUlx8vDW7LgyIkKFubvZdxAa543TVNdhW6sapjxqruA7wCr8C1F/A7wT1cjMvVq1zKJZ759Jn95dwOPzttI3q6DPHHNcLq1ax6W881bt4eyyiomDLbmedM01ZXgrWqY8aoU4GlV/QeAiMQ72466GpWpV3yc8JOL+zKsWxq3v7SCSY8u5M9XDWVs/9Cvzz4rr5CMVkkM79Ym5Mc2JhrU1QdvVcOMV72HL6H7pQDvuhSLOQUX9OvAjNvG0LVtc77zTC4Pz15PZQin0h0rq+T9dcWMG9iBuDhrnjdNk/cWbTamfsmqesT/xHkcnnZeEzZd2zbn1ZtHMvWMrjz2/ka+8fRi9h05EZJjL8gv5lh5JeMHWnEb03RZgveqldPhz4PgvjTf/crpbkfkJaUicrr/iYgMB7xRF9U0SHJiPA985TT+8JXTyN16gEmPLmLZ9gONPu6svCLSmidyVq+2IYjSmOhUb4IXkVtFxDqxImnldHjrNjhYAKjv/q3bLMl/4XbgZRFZKCILgWnArS7HZBrhyjO68urNI0mIF676+8c889FWVE+tyb6soop31+7mov4dSIy3axjTdAXz298RWCIi00VkvNh8k/B779dQXu2CtPyYb7tBVZfgGwR6M/B9oL+qLnU3KtNYgzJTefvWMeRkZXDvm6v54UsrKD1R0eDjfLhpL4ePV9joedPk1ZvgVfUXQBbwL+CbQL6I/M4WnAmjg7WMXaxte9N0BnAaMAy4WkS+4XI8JgRSmyfyj2+M4Kfj+vL2yl1MefxDNu45Uv8bA8zOK6JlUgKj+qSHKUpjokNQ7Vfqaysrcm4VQBvgFRH5Qxhja7pSuzRsexMjIv8FHgZG40v0ZwAjXA3KhExcnHDL+X149ttnsa+0jMmPLWLGysKg3ltRWcWcNbu5oF97khLiwxypMd5W1zx4AETkNuB6YC++YiI/VdVyEYkD8oE7whtiEzT2Hl+fe2AzfWKKb7sBXzIfoKfaSWuiwuisdGbcNprvP7+MW15YxrLtPblrQr86+9U/3bqf/aVlVnveGIK7gk8HrlDVcar6sqqWA6hqFTAprNE1VaddCZc+AqldAfHdX/qIb7sByMM3NsTEuE6pKUy74Ry+ObIH/1q0hauf+oTdh47Xuv+svCKSE+M4t29GBKM0xpvqvYLHVwZ0v/+Js+72AFVdrKprwxZZU3falZbQa5cOrBGRT4HPJ06r6mXuhWTCpVlCHPddNpBh3dK469VVXPLIQh69+nTO6d3upP2qqpTZq4s4NzuD5s2C+dNmTGwL5n/BE8DpAc9La9hmTCTd53YAJvImD82kf6fW3PTcUq755yfcMb4fN+b0+nwhmeUFJew+dIIJg6y4jTEQXBO9BPZ1Ok3zwfTd9xWRFQG3QyJyu4jcJyI7A7ZPrOX940VkvYhsFJG7gv+RTKxT1fnAViDRebwEsJUNm4DsDq1489bRTBjUiQdmruPG/y7l0PFywFd7PjFeuKB/e5ejNMYbgknwm0XkNhFJdG4/BDbX9yZVXa+qQ1V1KDAc30Igrzsv/9n/mqq+U/29zuIhjwMTgAH4pkENCPaHCgurLOcZIvI94BXg786mTOAN9yIykdQyKYHHvj6MX04awLx1e7js0UWs2XWImXlFjO6TTuvkRLdDNMYTgknwNwEjgZ3ADuAs4IYGnmcssElVtwW5/5nARlXdrKplwEvA5AaeM3SsspzX3AKMAg4BqGo+YJdtTYiI8J3RPXnxhrM5WlbJ5McXsePAMcbb6HljPhdMoZs9qjpVVduragdV/bqq7mngeaYCLwY8v1VEVorI07WUwc0ECgKe73C2ucMqy3nNCeeLHwAikgDYlLkm6IwebZlx2xiGd29Dy6QELhpgCd4Yv2D60pOB7wADgWT/dlX9djAnEJFmwGXA3c6mJ4Df4PuD/Bvgj0D1Y9VUDrfGP+AicgNOi0K3bt2CCanhrLKc18wXkZ8BKSJyEb5ytW+5HJNxSUarJF747tkcKauw5nljAgTTRP9ffHOOxwHzgS7A4QacYwKwTFV3A6jqblWtdAbr/QNfc3x1O4CuAc+7ALtqOriqPqWqI1R1REZGmOa+WmU5r7kLKAZWATfim8r5i8YcUETSROQVEVknImtF5JwQxGkiJC5OLLkbU00wCb6Pqv4SKFXVZ4BLgMENOMfVBDTPi0jgHJbL8RUtqW4JkCUiPZ0WgKnAmw04Z2iNvcdXSS6QVZZzhTMA81lV/Yeqfk1Vv+o8bmwT/V+BWaraDxgCWI0HY0xUCybBlzv3JSIyCEgFegRzcBFpDlwEvBaw+Q8iskpEVgLnAz9y9u0sIu8AqGoFvuU/Z+P7QztdVVcHc86wsMpynqGqlUCG88UvJESkNZCDb0ElVLVMVUtCdXxjjHFDMIVunnIGwv0C31V0S+CXwRxcVY8C7aptu66WfXcBEwOev4Ov6dUbrLKcl2wFPhSRN/EVXgJAVf90isfrha/J/98iMgRYCvxQVUvrfpsxxnhXnVfwzoIyh1T1gKouUNVezmj6v9f1PmPCbBfwNr7f31YBt1OVgK8y4xOqOgzfl4YvFVcSkRtEJFdEcouLixtxOmOMCb86r+BVtUpEbgVswrfxDFX9FYCItAjRVfYOYIeqLnaev0INCV5VnwKeAhgxYoRNyzPGeFowffBzReT/RKSriLT138IemTG1EJFzRGQNzkA4ERkiIn871eOpahFQICJ9nU1jgTWNj9QYY9wTTB+8f476LQHbFF+/pTFu+Au+aZtvAqjqZyKS08hj/gB43hm8txn4ViOPZ4wxrqo3watqz0gEYkxDqGqBfxUxR2Ujj7cCGNGooIwxxkOCqWT3jZq2q+qzoQ/HmKAUiMhIQJ0r7tuweevGGHOSYJrozwh4nIyvf3IZYAm+upXTffXpD+7wVbkbe49NrQuPm/AVpsnEN0BuDid3IRljTJMXTBP9DwKfi0gqvvK1JpB/xTn/ojT+FefAknyIiMiDqnoncL6qXuN2PMYY42XBjKKv7iiQFepAop6tOBcJE0UkkS8WLjLGGFOLYPrg3+KLldzigAHYvPgvsxXnImEWsBdoISKH8K06qP57VW3tZnDGGOMlwfTBPxzwuALYpqqWtapL7eJrlq9puwmVX6jqT0Xkf6o62e1gjDHGy4Jpot8OLFbV+ar6IbBPRHqENapoZCvORcLHzv0hV6MwxpgoEMwV/MvAyIDnlc62M2revYnyD6SzUfTh1ExErgdGisgV1V9U1ddqeI8xxjRJwST4BFUt8z9R1bJQLtUZU2zFuXC7CbgGSAMurfaacvKyxMYY06QFk+CLReQyVX0TQEQm4xvoZExEqeoiYJGI5Krqv9yOxxhjvCyYBH8TvhrdjznPdwA1VrczJpxE5AJVnQccsCZ6Y4ypWzCFbjYBZ4tIS0BU9XD4wzKmRucC8/hy8zxYE70xxpwkmHnwvwP+oKolzvM2wE9U9RfhDs6YQKp6r3NvK70ZY0w9gmmin6CqP/M/UdUDIjIRsARvIkpEflzX66r6p0jFYowxXhdMgo8XkSRVPQEgIilAUnjDMqZGrZz7vvimab7pPL8UWOBKRMYY41HBJPjngPdE5N/4+jm/ja0kZ1ygqr8CEJE5wOn+8SAich++2gzGGGMcwQyy+4OIrAQuxFfz+zeqOjvskYWKLeEai7oBZQHPy4Ae7oRijDHeFMwVPKo6C99CH4jIKBF5XFW9v/62LeEaq/4LfCoir+NrVboceMbdkIwxxluCWi5WRIaKyIMishW4H1gX1qhCxZZwjUmq+lvgW8ABoAT4lqr+3t2ojDHGW2q9gheRbGAqcDWwD5iGbx78+RGKrfFsCdeYparLgGVux2GMMV5V1xX8OmAscKmqjlbVR/EtNBM9aluq1ZZwNcYYE+PqSvBfAYqA90XkHyIyFt8gu6CISF8RWRFwOyQit4vIQyKyTkRWisjrIpJWy/u3isgq5725DfuxHF5dwnXldPjzILgvzXe/crq78RhjjIk5tSZ4VX1dVa8C+gEfAD8COojIEyJycX0HVtX1qjpUVYcCw4GjwOvAXGCQqp4GbADuruMw5zvHGBH0TxTotCvh0kcgtSsgvvtLH3F3gJ1/4N/BAkC/GPhnST5oInKrU1HRGGNMLYKZJlcKPI9vwZm2wNeAu4A5DTjPWGCTqm4DtgVs/wT4agOO03BeW8K1roF/XorT2zoCS0RkGfA0MFtV1eWYjDHGU4IaRe+nqvtV9e+qekEDzzMVeLGG7d8GZtZ2OmCOiCwVkRsaeD7vsoF/jeasg5AF/Av4JpAvIr8Tkd6uBmaMMR7SoAR/KkSkGXAZ1SqNicjPgQp8rQM1GaWqpwMTgFtEJKeW498gIrkikltcXBzCyMPEBv6FhHPFXuTcKoA2wCsi8gdXAzPGGI8Ie4LHl6CXqepu/wYRuR6YBFxTW9Oqqu5y7vfg67s/s5b9nlLVEao6IiMjI+TBh1ykB/7F4IA+EblNRJYCfwA+BAar6s34xnp8xdXgjDHGI4KqZNdIVxPQPC8i44E7gXNV9WhNbxCRFkCcqh52Hl8MxEZ1Gn8/eyTK58ZuJb904ApnTMfnVLVKRCa5FJMxxnhKWBO8iDQHLgJuDNj8GL7V6OaKCMAnqnqTiHQG/qmqE4EOwOvO6wnAC0653NgQqYF/sTug7x1gv/+JiLQCBqjqYlVd615YxhjjHWFN8M4Vertq2/rUsu8uYKLzeDMwJJyxNQmxO6DvCeD0gOelNWwzxpgmLRJ98MYtsTugTwLHbqhqFZHpbjLGmKhhCT6WebWSX+NtdgbaJTq3HwKb3Q7KGGO8xBJ8LPNiJb/QuAkYCewEdgBnAbFTK8EYY0LAmjVjndcq+YWAM3VyqttxGGOMl1mCN1FHRJKB7wADgWT/dlX9tmtBGWOMx1gTvYlG/8VXj34cMB/oAhxu7EFFJF5ElovI2409ljHGuM0SvIlGfVT1l0Cpqj4DXAIMDsFxfwjYPHpjTEywBG+iUblzXyIig4BUoEdjDigiXfB9Ufhn40IzxhhvsD54E42ectaD/wXwJtAS+GUjj/kX4A6gVSOPY4wxnmAJ3kQVEYkDDqnqAWAB0CsEx5wE7FHVpSJyXh373YAzHa9bt26NPa0xxoSVNdFDTK64FqucqnW3hviwo4DLRGQr8BJwgYg8V8O5dHpzKQAAFfJJREFUo2vlQmNMk2YJ3r/i2sECQL9Ycc2SvJfNFZH/E5GuItLWfzvVg6nq3araRVV74JtfP09Vrw1ZtMYY4wJroo/dFddimX+++y0B25QQNNcbY0yssAQfuyuuxSxV7RnGY38AfBCu4xtjTKRYgk/t4jTP17DdeJKIfKOm7ar6bKRjMcYYr7I++Gheca3pDg48I+A2BrgPuMzNgIwxxmvsCt7fz/7er33N8qldfMnd6/3v/sGB/vED/sGB4P3YG0lVfxD4XERS8ZWvNcYY47AED9G54poNDgx0FMhyOwhjjPESS/DRqgkPDhSRt/CNmgdfN9MAoMn0TxhjTDAswUerpj048OGAxxXANlWN/W82xhjTADbILlpF8+DAxtsOLFbV+ar6IbBPRHq4G5IxxniLJfhoddqVcOkjkNoVEN/9pY80lf73l4GqgOeVzjZjjDEOa6KPZtE4ODA0ElS1zP9EVctEpJmbARljjNfYFbyJRsUi8vm8dxGZDOx1MR5jjPEcu4I30egm4HkRecx5vgOosbqdMcY0VZbgTdRR1U3A2SLSEhBVPex2TMYY4zVha6IXkb4isiLgdkhEbneW9pwrIvnOfZta3n+9s0++iFwfrjhN9BGR34lImqoeUdXDItJGRO53Oy5jjPGSsCV4VV2vqkNVdSgwHF+1sdeBu4D3VDULeM95fhJnbe97gbOAM4F7a/siYJqkCapa4n+iqgeAiS7GY4wxnhOpQXZjgU2qug2YDDzjbH8GmFLD/uOAuaq63/njPRcYH5FITTSIF5Ek/xMRSQGS6tjfGGOanEj1wU/9//buPUiq8szj+PfHJeKNYKm46sAO7hIiwsgQxqi4rIJCjIawuybR0oi6FKxKvCRrCrc2QryFqLlXonHVSBV4oQAVMcqsgktILIFxyCABAmsQBhFxuEQQVsBn/+gz2ODMMDM9TU93/z5V1Jw+1+dMn+F5z3ve877Ak8n0SRGxESAiNkrq3sD6pwLp3bTVJvM+RdJYYCxAz5492yzgolIzPd8G25kKvCLpN6S6rL0O8FCxZmZpsp7gk/eTRwK3t2SzBuZFA/OIiIeBhwEGDRrU4DrWhDwclS4i7pNUA1xI6lq5KyLm5jgsM7N25XBU0V8MvBERm5LPmySdDJD8fK+BbWqBHmmfS4B3shplsWpqVLp2LCJeioh/j4jvADsk/TLXMZmZtSeHI8FfwSfV8wCzgfpW8aOB5xrYZi4wPGkdfRwwPJlnbS1PR6WTNEDSDyWtBe4GVuY4JDOzdiWrVfSSjgIuAsalzZ4MTJf0r6QGDflasu4g4N8iYkxEbJF0F7A42ebOiNiSzViLVh6NSifpc6Tac1wB1AFPk3oP/oKcBmZm1g5lNcFHxIfA8QfNqyPVqv7gdZcAY9I+PwY8ls34jFSDuvRn8NDwqHTtoyHeSuB3wFciYg2ApFsPdxBmZvnAfdEXu+aMSlffEG/7eiA+aYhXM/1wR/svwLvAfEn/JWkYDTfINDMreu6q1g49Kl1TDfEO4118RDwDPCPpaFL9J9wKnCTpQeCZiKg8bMGYmbVzvoO3Q2tnDfEiYmdETIuIS0m9YbGUBnpENDMrZk7wdmiNNbhrBw3xkt4Ofx0RQ3Mdi5lZe+IEb4c27I5Uw7t0DTXEMzOzdsMJ3g6tOQ3xzMysXXEjO2ueQzXEMzOzdsV38GZmZgXICd7MzKwAOcGbmZkVICd4MzOzAuQEb2ZmVoCc4M3MzAqQE7wVPUk9JM2XtELSckk35zomM7NM+T14M9gLfCci3pB0LFAl6b8j4k+5DszMrLV8B29FLyI2RsQbyfQHwArg1NxGZWaWGSd4szSSSoFy4PXcRmJmlhkneLOEpGOAmcAtEfHXBpaPlbRE0pLNmzcf/gDNzFrACd4MkNSZVHKfFhGzGlonIh6OiEERMejEE088vAGambWQE7wVPUkCHgVWRMSPcx2PmVlbcII3g8HAN4GhkpYm/76c66DMzDLh1+Ss6EXEQkC5jsPMrC35Dt7MzKwAFd0d/J49e6itrWX37t25DiVvdOnShZKSEjp37pzrUMzMrJmKLsHX1tZy7LHHUlpaSqptlTUlIqirq6O2tpZevXrlOhwzM2umrCZ4Sd2AR4B+QADXAbcAfZJVugHbImJAA9uuBT4A9gF7I2JQW8S0e/duJ/cWkMTxxx+P3/s2M8sv2b6D/xnwUkRcJukzwFER8Y36hZJ+BGxvYvsLIuL9tg7Kyb1l/PsyM8s/WWtkJ6krMITU+8VExEcRsS1tuYCvA09mK4ZiUFpayvvvN10Gas46ZmZWWLLZiv40YDPwG0nVkh6RdHTa8n8ANkXE6ka2D6BSUpWksVmM08zMrOBkM8F3AgYCD0ZEObATmJC2/AqavnsfHBEDgYuBGyUNaWilbPcP/mz1BgZPnkevCS8wePI8nq3ekPE+165dy+c//3nGjBlDv379uPLKK3n55ZcZPHgwvXv3ZtGiRWzZsoVRo0ZRVlbG2WefTU1NDQB1dXUMHz6c8vJyxo0bR0Ts3+/UqVM566yzGDBgAOPGjWPfvn0Zx2pmZvkpmwm+FqiNiPpRuWaQSvhI6gT8M/B0YxtHxDvJz/eAZ4CzGlkva/2DP1u9gdtnLWPDtl0EsGHbLm6ftaxNkvyaNWu4+eabqampYeXKlTzxxBMsXLiQBx54gHvvvZeJEydSXl5OTU0N9957L1dffTUA3//+9znvvPOorq5m5MiRrFu3DoAVK1bw9NNP8/vf/56lS5fSsWNHpk2blnGcZmaWn7LWyC4i3pW0XlKfiFgFDAP+lCy+EFgZEbUNbZtU5XeIiA+S6eHAndmKtTH3z13Frj0H3gXv2rOP++euYlR5ZsOF9+rVi/79+wNwxhlnMGzYMCTRv39/1q5dy9tvv83MmTMBGDp0KHV1dWzfvp0FCxYwa1ZqLJRLLrmE4447DoBXXnmFqqoqKioqUnHu2kX37t0zitHMzPJXtlvRfwuYlrSgfwu4Npl/OQdVz0s6BXgkIr4MnAQ8k7Te7gQ8EREvZTnWT3ln264WzW+JI444Yv90hw4d9n/u0KEDe/fupVOnT3819a3ZG2rVHhGMHj2aH/zgBxnHZmZm+S+rXdVGxNKk+rwsIkZFxNZk/jUR8dBB676TJHci4q2IODP5d0ZE3JPNOBtzSrcjWzS/LQ0ZMmR/Ffurr77KCSecQNeuXQ+Y/+KLL7J161YAhg0bxowZM3jvvfcA2LJlC2+//XbW4zQzs/bJfdE34bYRfTiyc8cD5h3ZuSO3jejTyBZtZ9KkSSxZsoSysjImTJjAlClTAJg4cSILFixg4MCBVFZW0rNnTwD69u3L3XffzfDhwykrK+Oiiy5i48aNWY/TzMzaJ6W3ws53gwYNiiVLljS5zooVKzj99NObvc9nqzdw/9xVvLNtF6d0O5LbRvTJ+Pl7Pmru701SVVv1OtieNedaMzPLtqb+zy26vuhbalT5qUWZ0M3MLL+5it7MzKwAOcGbmZkVICd4MzOzAuQEb2ZmVoCc4M3MzAqQE3wOrF27ln79+mVl36+++iqXXnopALNnz2by5MlZOY6ZmbVvfk2ugI0cOZKRI0fmOgwzM8sB38EfSs10+Ek/mNQt9bNmepvsdu/evYwePZqysjIuu+wyPvzwQ+68804qKiro168fY8eO3T8U7M9//nP69u1LWVkZl19+OQA7d+7kuuuuo6KigvLycp577rlPHePxxx9n/PjxAFxzzTXcdNNNnHvuuZx22mnMmDFj/3r3338/FRUVlJWVMXHixDY5PzMzyy0n+KbUTIfnb4Lt64FI/Xz+pjZJ8qtWrWLs2LHU1NTQtWtXfvWrXzF+/HgWL17Mm2++ya5du5gzZw4AkydPprq6mpqaGh56KNWF/z333MPQoUNZvHgx8+fP57bbbmPnzp1NHnPjxo0sXLiQOXPmMGHCBAAqKytZvXo1ixYtYunSpVRVVbFgwYKMz8/MzHLLCb4pr9wJew4aOW7PrtT8DPXo0YPBgwcDcNVVV7Fw4ULmz5/PF7/4Rfr378+8efNYvnw5AGVlZVx55ZVMnTp1/yhzlZWVTJ48mQEDBnD++eeze/fu/WPDN2bUqFF06NCBvn37smnTpv37qayspLy8nIEDB7Jy5UpWr16d8fmZmVlu+Rl8U7Y3OFx94/Nb4OAhXyVxww03sGTJEnr06MGkSZPYvXs3AC+88AILFixg9uzZ3HXXXSxfvpyIYObMmfTpc+DAN/WJuyHpQ9TWV/9HBLfffjvjxo3L+JzMzKz98B18Uz5b0rL5LbBu3Tpee+01AJ588knOO+88AE444QR27Nix/xn5xx9/zPr167ngggu477772LZtGzt27GDEiBH84he/2J+oq6urWxXHiBEjeOyxx9ixYwcAGzZs2D/krJmZ5S/fwTdl2B2pZ+7p1fSdj0zNz9Dpp5/OlClTGDduHL179+b6669n69at9O/fn9LSUioqKgDYt28fV111Fdu3byciuPXWW+nWrRvf+973uOWWWygrKyMiKC0t3f/MviWGDx/OihUrOOeccwA45phjmDp1Kt27d8/4HM3MLHc8XOyh1ExPPXPfXpu6cx92B5R9PcNI84+Hiz2Qh4s1s/bAw8VmouzrRZnQzcwsv/kZvJmZWQFygjczMytARZngC6ndweHg35eZWf4pugTfpUsX6urqnLSaKSKoq6ujS5cuuQ7FzMxaoOga2ZWUlFBbW8vmzZtzHUre6NKlCyUlmb/7b2Zmh0/RJfjOnTvTq1evXIdh7YykLwE/AzoCj0SEx9k1s7xWdFX0ZgeT1BH4JXAx0Be4QlLf3EZlZpYZJ3gzOAtYExFvRcRHwFPAV3Mck5lZRpzgzeBUYH3a59pknplZ3iqoZ/BVVVXvS3q7kcUnAO8fznjaSL7F/be5DqAV1MC8T71mIWksMDb5uEPSqhYeJ9++y0z4XAtTsZxrPp1no//nFlSCj4gTG1smaUk+9pGer3HnmVqgR9rnEuCdg1eKiIeBh1t7kGL6Ln2uhalYzrVQztNV9GawGOgtqZekzwCXA7NzHJOZWUYK6g7erDUiYq+k8cBcUq/JPRYRy3MclplZRoopwbe6ajXH8jXuvBIRvwV+m+XDFNN36XMtTMVyrgVxngU1HryZmZml+Bm8mZlZASr4BC/pS5JWSVojaUKu42kuSWslLZO0VNKSXMdjmcnX67AlJPWQNF/SCknLJd2c65iyTVJHSdWS5uQ6lmyS1E3SDEkrk+/3nFzHlC2Sbk2u3zclPSkpb0faKugq+qQL0j8DF5F6FWoxcEVE/CmngTWDpLXAoIjIl3cxrRH5fB22hKSTgZMj4g1JxwJVwKhCO890kr4NDAK6RsSluY4nWyRNAX4XEY8kb5ocFRHbch1XW5N0KrAQ6BsRuyRNB34bEY/nNrLWKfQ7eHdBau1BUVyHEbExIt5Ipj8AVlDAPQJKKgEuAR7JdSzZJKkrMAR4FCAiPirE5J6mE3CkpE7AUTTQJ0a+KPQEn89dkAZQKakq6UHN8lc+X4etIqkUKAdez20kWfVT4LvAx7kOJMtOAzYDv0keRzwi6ehcB5UNEbEBeABYB2wEtkdEZW6jar1CT/DN6oK0nRocEQNJjXB2o6QhuQ7IWi2fr8MWk3QMMBO4JSL+mut4skHSpcB7EVGV61gOg07AQODBiCgHdgKF2o7kOFK1a72AU4CjJV2V26har9ATfLO6IG2PIuKd5Od7wDOkqnktP+XtddhSkjqTSu7TImJWruPJosHAyKStzFPAUElTcxtS1tQCtRFRXxszg1TCL0QXAn+JiM0RsQeYBZyb45hardATfF52QSrp6KSREklV2HDgzdxGZRnIy+uwpSSJ1HPaFRHx41zHk00RcXtElEREKanvc15E5O2dXlMi4l1gvaQ+yaxhQKE2nFwHnC3pqOR6HkaqLUleKuie7PK4C9KTgGdS1xedgCci4qXchmStlcfXYUsNBr4JLJO0NJn3H0kvgZbfvgVMSwqobwHX5jierIiI1yXNAN4A9gLV5HGvdgX9mpyZmVmxKvQqejMzs6LkBG9mZlaAnODNzMwKkBO8mZlZAXKCNzMzK0BO8GkkHZ+M3rZU0ruSNqR9/sxB686tf1e9if3VSuqW3agPON5USaMO1/HMzKz9coJPExF1ETEgIgYADwE/qf+cDBKCUjpExIhkQI2CkAysYFYU8r0wn3bchZIGHO7jWn5wgm8GSX+fjA38EKkOEE5O/4OW9HwyKMxySWMOsa9OkrZJmizpj5Jek9Q9WXbAHbikHcnPC5NxtmdIWi3pbklXS1osqSYZ2KPeCEm/k/RnSRenHfPHkhYl649J2+/Lkp4i1aGDWVEo5sJ8PRfqC58TfPP1BR6NiPJkxKF0oyPiC0AF8O1kwIKmfBb4n4g4E3gNuK4Zxz8TuBHoD4wBSiOiApgCjE9brwfwj8BXgIclHQGMJTUwxllJjDdK6pmsfzbw3Yjo34wYzApaWxbmG9jno8l2L0rqkizbfwcu6W8krUmmx0iaJWmOpL9Iul7SbUqN5vaHg2oLrkluFJZJGpRsf4ykx5NCfbWkr6Tt9ylJc4AX2+r3Zu2TE3zz/W9ELG5k2a2S/kgqWZcAf3eIfe2KiPo/riqgtBnHfz0iNkXEblJdRc5N5i87aPvpEfFxRKwiNURpb1J92V+bdB/6OtAtmQ/wWkSsa8bxzYpFWxbm6/UBfhoRZwC7gOa0lTkD+AapQvgPga3JaG5VQHq/90dExDnAzXwyNv0dwEtJoX4o8KP6QgVwDvDNiLiombFbnnIVTfPtbGimpAuBIcDZEbFL0kKgS0PrpvkobXofn3wPe0kKXZI6cuD3839p0x+nff74oPUO7ns4SA1XekNEvNJA7A2el1kRO1RhfmQyXV+YX9KMfa6JiGXJdHML9fMiYiewM3lc93wyfxnwubT1ngSIiHmSuis1XO9w4GJJ9cO6dgHqa+0qI2JrM45vec538Jn7LLAlSe5nkCrZt9Za4AvJ9D+RGpikpb6WPDv8HKnq+tWk7vZvqH/mJqmPpCMziNOskDWnMH8mUMOhC/P10gvoDRbqG9hXpoX6UWntCnpGxJ+T5S7UFwkn+My9AByVVNHfQaoKvLV+DVwkaREwgAP/wJtrDbCAVGl/bNJg6NekEv1SSW8CD+LaG7OWasvCfL21fFKov6yV+/gGgKTzgU3JXf9c4Kb6FSSVtz5Ey1f+T74RETEpbXoNqYSbvrwk7eOIRvZR0sC8vaSegdd/fgp4KpneCJyVtvp/JvNfBl5O2+a8tOn9yxobjzoi9gETkn/pDtivmTXpBWBsUphfSWaF+Xr3A09LuhaY38p9/FXSH4Bj+WQY1+8DP5W0jNSN3Brgq5kGa/nFw8WamZkVIFfRm5mZFSBX0ZuZtYGkw6rKBhadHxHbDnc8Zq6iNzMzK0CuojczMytATvBmZmYFyAnezMysADnBm5mZFSAneDMzswL0/zg8SHij65l+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(ncols=2, figsize=(8, 5), gridspec_kw={'width_ratios': [1,3], 'wspace':0.3})\n",
    "\n",
    "ax = axs[0]\n",
    "trial_numbers = np.arange(len(scores))\n",
    "ax.scatter(trial_numbers, scores)\n",
    "ax.scatter(trial_numbers, baselines)\n",
    "# ax.set_xlim(right=2*len(scores))\n",
    "ax.legend([\"model\", \"baseline\"], loc='lower center')\n",
    "ax.set_xlabel(\"Trial number\")\n",
    "ax.set_ylabel(\"Accuracy [%]\")\n",
    "ax.set_ylim(bottom=ax.get_ylim()[0]*0.95)\n",
    "\n",
    "ax = axs[1]\n",
    "ax.plot(np.array(scores - np.array(baselines)))\n",
    "ax.set_title(\"Improvement over baseline\")\n",
    "ax.set_ylabel(\"Accuracy difference [% points]\")\n",
    "ax.set_ylim(bottom=0)\n",
    "ax.set_xlabel(\"Trial_number\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model generalizes medium well: Its test error (1-accuracy) is 34% greater than its validation error. In particular, the accuracy of 79% on the test set is lower than the accuracy of 89% obtained in Annie's paper.\n",
    "\n",
    "Finally, we display the intercept and the coefficients and compare them with the model described in Annie's paper:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Intercept: [-2.82815156]\n",
      "Coefficients:\n",
      " {'quarantine': 1.819751771499741, 'quarantine_missing': 0.353819697770931, 'conjonctivitis_triage': 1.3966816480481956, 'diarrhea_triage': 1.0488543066967129, 'referral(4-9)': 1.0816016292060706, 'pyrexia_triage': 0.3769283857926169, 'dysphagia_triage': 0.5626054205326325, 'hemorrhage_triage': 0.3764649262842281, 'myalgia_triage': -0.36070330201615547}\n"
     ]
    }
   ],
   "source": [
    "print('Intercept:', model.intercept_)\n",
    "print('Coefficients:\\n', dict(zip(predictors, model.coef_.flat)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAFVCAYAAADsRFQ0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZgkVZX+8e/L7gIqiy2y2KggIiJij4D400YFQQWHARdEBUQZxw1HHZdhRsBlxn0GcUEUUHFhcWMRFVwaREHshgZkU2Rw6AFFBQEFxYb398eNhOzqrF7IjBst8X6ep56qjMyqE5GVdSryxr3nyDYREXHft1LXOxAREXUk4UdE9EQSfkRETyThR0T0RBJ+RERPJOFHRPTE2Alf0kaSfiDpckmXSjpoxGMk6aOSrpJ0saRtxo0bERHLZ5UJ/IyFwJttXyBpTWCepDNtXzb0mF2BTZuPbYFPNp8jIqKSsRO+7euB65uvb5V0ObABMJzwnw983mWV13mSHixp/eZ7p7Xuuut65syZ4+5iRERvzJs373e21xt13yTO8O8maSbwROAnU+7aALh26PaCZtsSE/7MmTOZO3fuBPcwIuK+TdKvprtvYglf0gOBrwJvtH3L1LtHfMvImg6SDgQOBJgxYwZz5syZ1C5GRPTaRBK+pFUpyf6Ltr824iELgI2Gbm8IXDfqZ9k+CjgKYNasWZ49e/YkdjEiovcmMUtHwNHA5bY/Ms3DTgFe3szW2Q64eWnj9xERMVmTOMPfAXgZcImk+c22fwU2BrB9JHA68BzgKuA2YP8JxI2IiOUwiVk65zB6jH74MQZeO26siIi497LSNiKiJ5LwIyJ6Igk/IqInJrrwKoqZb/9m6zGued9zW48REfctOcOPiOiJJPyIiJ5Iwo+I6Ikk/IiInkjCj4joiST8iIieSMKPiOiJJPyIiJ6YSMKXdIykGyT9bJr7Z0u6WdL85uOdk4gbERHLblIrbT8LfAz4/BIe80Pbz5tQvIiIWE4TOcO3fTZw4yR+VkREtKPmGP72ki6S9C1Jj6sYNyIiqFc87QLgEbb/KOk5wDeATUc98L7QxPzNj1/Yeoy/xeclIrql0oxqAj9ImgmcZnvLZXjsNcAs279b0uNmzZrluXPnTmT/akq1zIjoiqR5tmeNuq/KkI6khzXNzpH05Cbu72vEjoiIYiJDOpK+DMwG1pW0ADgEWBXubmK+F/BPkhYCtwMv9qTeWkRExDKZSMK3vfdS7v8YZdpmRER0JB2vYmJy7SJixZbSChERPZGEHxHRE0n4ERE9kYQfEdETSfgRET2RhB8R0RNJ+BERPZGEHxHRE0n4ERE9kYQfEdETSfgRET2RhB8R0ROTKo98DPA84IZRDVCaWviHA88BbgP2s33BJGJHdC1F4+JvxaSqZX6WUv7489PcvyulpeGmwLbAJ5vPERORpBuxdBMZ0rF9NnDjEh7yfODzLs4DHixp/UnEjoiIZVOrHv4GwLVDtxc0266f+sA0MV82K+Lz0uVx9zV2xPKolfA1YtvIFoe2jwKOgtLEfPbs2S3uVjv2qzG8sM/s1mMsry6Pu6+xI5ZHrVk6C4CNhm5vCFxXKXZERFAv4Z8CvFzFdsDNthcbzomIiPZMalrml4HZwLqSFgCHAKsC2D4SOJ0yJfMqyrTM/ScRNyIilt1EEr7tvZdyv4HXTiJWRETcO1lpGxHRE0n4ERE9kYQfEdETSfgRET2RhB8R0RNJ+BERPZGEHxHRE0n4ERE9kYQfEdETSfgRET2RhB8R0RNJ+BERPTGRhC9pF0lXSrpK0ttH3L+fpN9Kmt98vHIScSMiYtmNXS1T0srAx4GdKI1OfirpFNuXTXnoCbZfN268iFgxpHH8355JnOE/GbjK9tW27wCOpzQtj4iIFcgkEv50Dcqn2lPSxZK+ImmjEfdHRESLVHqTjPEDpBcAz7b9yub2y4An23790GPWAf5o+y+SXg280PYzpvl5BwIHAsyYMeNJxx9//Fj714VL/u/m1mM8foMHtR5jeXV53H2N3aW+HveKbscdd5xne9ao+yaR8LcHDrX97Ob2OwBs/+c0j18ZuNH2Un+Ts2bN8ty5c8favy70dWyzy+Pua+wu9fW4V3SSpk34kxjS+SmwqaRNJK0GvJjStHx4B9Yfurk7cPkE4kZExHIYe5aO7YWSXgd8B1gZOMb2pZLeBcy1fQrwBkm7AwuBG4H9xo0bERHLZ1JNzE8HTp+y7Z1DX78DeMckYkVExL0zkYQfEd3o6zh6X497XEn49zH5Q4iI6aSWTkRETyThR0T0RBJ+RERPJOFHRPREEn5ERE8k4UdE9EQSfkRETyThR0T0RBJ+RERPJOFHRPRErSbmq0s6obn/J5JmTiJuREQsu7ET/lAT812BLYC9JW0x5WEHADfZfjTwX8D7x40bERHLp1YT8+cDn2u+/grwTEmaQOyIiFhGtZqY3/0Y2wuBm4F1JhA7IiKWUa0m5pc2j1nQ3P5l85jfj/h5E2lingbLEdGGFT23LKmJ+STq4S8ANhq6vSFw3TSPWSBpFeBBlFaHi7F9FHAUlCbms2fPvlc7tV+NuvD7zG49RkSsWP6Wc0uVJubN7X2br/cCvu9x31pERMRyqdXE/GjgOElXUc7sXzxu3IiIWD61mpj/GXjBJGJFRMS9k5W2ERE9kYQfEdETSfgRET2RhB8R0RNJ+BERPZGEHxHRE0n4ERE9kYQfEdETSfgRET2RhB8R0RMTKa0QEdEX17zvuV3vwr2WM/yIiJ5Iwo+I6ImxEr6ktSWdKekXzeeHTPO4OyXNbz6m1sqPiIgKxj3DfzvwPdubAt9rbo9yu+2tm4/dx4wZERH3wrgJ//nA55qvPwf8/Zg/LyIiWjLuLJ0Ztq8HsH29pIdO87g1JM0FFgLvs/2N6X7glCbmzJkz517t2Jsfv/Befd/yuLf7FhHRhaUmfEnfBR424q6DlyPOxravk/RI4PuSLrH9y1EPTBPziIh2LDXh237WdPdJ+o2k9Zuz+/WBG6b5Gdc1n6+WNAd4IjAy4UdERDvGHcM/Bdi3+Xpf4OSpD5D0EEmrN1+vC+wAXDZm3IiIWE7jJvz3ATtJ+gWwU3MbSbMkfaZ5zGOBuZIuAn5AGcNPwo+IqGysi7a2fw88c8T2ucArm69/DDx+nDgRETG+rLSNiOiJJPyIiJ5Iwo+I6Ikk/IiInkjCj4joiST8iIieSMKPiOiJJPyIiJ5Iwo+I6Ikk/IiInkjCj4joiST8iIieGLeJ+QskXSrpLkmzlvC4XSRdKekqSdP1vY2IiBaNe4b/M+AfgLOne4CklYGPA7sCWwB7S9pizLgREbGcxi2PfDmApCU97MnAVbavbh57PKX5eWriR0RUVGMMfwPg2qHbC5ptERFR0VhNzG0v1tJw1I8Ysc1LiHcgcCDAjBkzmDNnzjKEWNybH7/wXn3f8ri3+xYR0YWxmpgvowXARkO3NwSuW0K8o4CjAGbNmuXZs2ffq6D7vf2b9+r7lsc1+8xuPUZExKTUGNL5KbCppE0krQa8mNL8PCIiKhp3WuYekhYA2wPflPSdZvvDJZ0OYHsh8DrgO8DlwIm2Lx1vtyMiYnmNO0vn68DXR2y/DnjO0O3TgdPHiRUREePJStuIiJ5Iwo+I6Ikk/IiInkjCj4joiST8iIieSMKPiOiJJPyIiJ5Iwo+I6Ikk/IiInkjCj4joiST8iIieSMKPiOiJWk3Mr5F0iaT5kuaOEzMiIu6dsaplck8T808tw2N3tP27MeNFRMS9VKOJeURErABqjeEbOEPSvKZnbUREVFajiTnADravk/RQ4ExJV9g+e5p4aWIeEdGCGk3MBx2wsH2DpK8DTwZGJvw0MY+IaEfrQzqSHiBpzcHXwM6Ui70REVFR603MgRnAOZIuAs4Hvmn72+PEjYiI5dd6E3PbVwNPGCdORESMLyttIyJ6Igk/IqInkvAjInoiCT8ioieS8CMieiIJPyKiJ5LwIyJ6Igk/IqInkvAjInoiCT8ioieS8CMieiIJPyKiJ8atlvlBSVdIuljS1yU9eJrH7SLpSklXSXr7ODEjIuLeGfcM/0xgS9tbAT8H3jH1AZJWBj4O7ApsAewtaYsx40ZExHIaK+HbPsP2oJfgecCGIx72ZOAq21fbvgM4Hnj+OHEjImL5TXIM/xXAt0Zs3wC4duj2gmZbRERUNJEm5pIOBhYCXxz1I0Zs8xLipYl5REQLxm5iLmlf4HnAM22PSuQLgI2Gbm8IXLeEeGliHhHRgnFn6ewCvA3Y3fZt0zzsp8CmkjaRtBrwYuCUceJGRMTyG3cM/2PAmsCZkuZLOhIWbWLeXNR9HfAd4HLgRNuXjhk3IiKW07hNzB89zfa7m5g3t08HTh8nVkREjCcrbSMieiIJPyKiJ5LwIyJ6Igk/IqInkvAjInoiCT8ioieS8CMieiIJPyKiJ5LwIyJ6Igk/IqInkvAjInoiCT8ioieS8CMiemKsapmSPgjsBtwB/BLY3/YfRjzuGuBW4E5goe1Z48SNiIjlN+4Z/pnAlra3An4OvGMJj93R9tZJ9hER3Rgr4ds+o2lwAnAepX1hRESsgMYa0pniFcAJ09xn4AxJBj7V9K0dKU3MIyLaodF9x4ceIH0XeNiIuw62fXLzmIOBWcA/jGpkLunhtq+T9FDKMNDrbZ+9tJ2bNWuW586duwyHsbiZNZqYv++5rceIiFgekuZNN3S+1DN8289ayg/fF3ge8MxRyb75Gdc1n2+Q9HXgycBSE35EREzOWGP4knYB3gbsbvu2aR7zAElrDr4GdgZ+Nk7ciIhYfuPO0vkYsCZwpqT5ko6EMoQjadC0fAZwjqSLgPOBb9r+9phxIyJiOY110db2o6fZfh3wnObrq4EnjBMnIiLGl5W2ERE9kYQfEdETSfgRET2RhB8R0RNJ+BERPZGEHxHRE0n4ERE9kYQfEdETk6yWuUJJYbOIiEXlDD8ioieS8CMieiIJPyKiJ8ZO+JLeLeniplrmGZIePs3j9pX0i+Zj33HjRkTE8pnEGf4HbW9le2vgNOCdUx8gaW3gEGBbSvOTQyQ9ZAKxIyJiGY2d8G3fMnTzAZT+tVM9GzjT9o22b6K0Odxl3NgREbHsJjItU9J7gZcDNwM7jnjIBsC1Q7cXNNsiIqKSZUr4S2tkbvtg4GBJ7wBeRxm+WeRHjPjekf1vJR0IHAgwY8YM5syZsyy7GBERS7FMCX9pjcyHfAn4Josn/AXA7KHbGwJzpol1FHAUwKxZszx79uxRD4uIiOUke+SJ9rL/AGlT279ovn498HTbe015zNrAPGCbZtMFwJNs37iUn/1b4Fdj7eCyWxf4XaVYid3v2F3HT+z7duxH2F5v1B2TGMN/n6THAHdRkvOrASTNAl5t+5W2b5T0buCnzfe8a2nJHmC6nW6DpLm2Z9WKl9j9jd11/MTuV+xhYyd823tOs30u8Mqh28cAx4wbLyIi7p2stI2I6Ikk/HscldiJ3ZP4id2v2Hcb+6JtRET8bcgZfkRETyThR0T0RBJ+RERPJOF3SNL9mjUMvSPpAR3F7fNz3stjl/QISc9qvr6fpDW73qeu9DrhS/resmxrKfZuwHzg283trSWdUin2m0Z8HCBp6wqxnyLpMuDy5vYTJH2i7bhNrM6e8ybeQyXtIem1kl4h6cmSqvwNdvx6k6SXSnpnc3tjSU+uFPtVwFeATzWbNgS+USn2DElHS/pWc3sLSQfUiD2dXiZ8SWs05R7WlfQQSWs3HzOBkQ1cWnAopTfAHwBszwdmVoo9i7IieoPm40BKraNPS3pry7H/i1Iu+/cAti8CntZyzIFD6eA5l7SjpO9Q6kztCqwPbAH8G3CJpMMkrdXybhxKd6+3TwDbA3s3t28FPl4p9muBHYBbAJoyMA+tFPuzwHe4J6f8HHhjpdgjTaQ88t+gf6Q88Q+n1PUZuIV6L8SFtm+WRhUSbd06wDa2/wgg6RDKWdDTKDWPPtBmcNvXTjnuO9uMN6Sr5/w5wKts/+/UOyStAjwP2An4aov70OXrbVvb20i6EMD2TZJWqxT7L7bvGBx383zXmou+ru0TmyrC2F4oqdZrfaReJnzbhwOHS3q97SM62o2fSXoJsLKkTYE3AD+uFHtj4I6h23+lFFy6XdJfWo59raSnAG7+6N9AM7xTQSfPue1/WcJ9C6kzxNDl6+2vklamSbSS1qPU3qrhLEn/CtxP0k7Aa4BTK8X+k6R1uOe4t6P0DOlMrxdeSbof8E/AUym/lB8CR9r+c4XY9wcOBnam9Av4DvDuSrH/HdgDOLnZtBtwCvBh4Cjb+7QYe13gcOBZlOM+AzjI9u/bijkUu8vn/FGU53wjYCHwC+DLtqskgI6PfR/gRZRquZ8D9gL+zfZJFWKvBBzAosf9GVdIfJK2AY4AtgR+BqwH7GX74rZjT7tPPU/4J1LGE7/QbNobeIjtF3S3V3U01Ux3oPwRnNMUu4sWSDqIMmxzFmV4Zz5wE+UfwGtsz+lu7+qQtDnwTMrr7Xu2a72r61QzhPQYynFfafuvne5PzxP+RbafsLRtLcXeDHgL5cLZ3UNrtp/Rduwm/srAjCmxFxtjbiHuR0dsvhmYa/vkEfdNMvapLD5+ezMwF/hUW2e7ki4BtrZ9Z3Omfbrt2ZI2Bk62/cQ24k7Zh06OvYm99ojNt9ZIfs1zP91xv6fNd5aS/mHE5puBS2zf0FbcJenlGP6QCyVtZ/s8AEnbAj+qFPsk4EjgM9S7aAnc3ajmEOA3TWxR/ii2qhB+DWBzyvED7AlcChwgaUfbbc5iuJrytvrLze0XUZ6DzYBPAy9rMfYqlOd6dWBNKP9gJa3aYsxhXR77BZShrJsor7UHA9dLuoFyMXtei7G/RXnev9TcfnHz+RbKLJrdWox9AGV20g+a27OB84DNJL3L9nEtxh7Ndm8/KBcL7wKuaT7uoiSfS4CLW449r8PjvgpYp6PY3wdWGbq9SrNtZeCylmOfPd024NIW4x4EXEypmHgFsH+zfb1R+3RfOvbm5x8JPHvo9s7AR4DtgJ+0HPtH022jnGm3GftUYMbQ7RnA14C1gZ/V+L1P/ej7Gf4uHcY+VdJrgK8Dd8+M8TJ0ApuAa+lutsAGwAOG4j8AeLjLcEfbM4TWk7Sxm6GrZkhl3ea+O6b/tvHYPlzSd4HHAh+xfUWz/bfUW4PQybE3Ztl+9eCG7TMk/YftN0laveXYD5S0re2fADQLvh7Y3Lew5dgzbf9m6PYNwGYuHQA7GcvvdcK3/StJTwU2tX1sM4NkTdv/UyH8vs3n4Sl7Bh5ZIfbVwBxJ32TRfzYfqRD7A8B8SXMob++fBvyHSqmF77Yc+83AOZJ+2cTeBHhNE/tzbQa2fSnl3SMAkg60XbNGemfHDtwo6W3A8c3tFwE3NdeR2p6e+UrgGEkPpBz3LcArm+P+z5Zj/1DSaSw6fHl2E/sPLcceqe8XbQ+hrDp9jO3NJD0cOMn2Dh3vWqua416M7cMqxV+fsupTwPm2r6sRt4m9OuUagoArXGFa4jT7cYHtbSrH7OTYmxOpQyjTnwWcAxxGeZe3se2rKuzDgyj5rlqiVVnttSdDs+GAr7rDpNv3hD8feCJwgZuZEpIutt3axUtJz7D9/Wmu4GP7a23FXlFIegiwKeUCLgC2z64Ue0tKWYPh2J+vEXvKflzoCrNzpsRcIY69NknPBR7Hosf9ru72qDu9HtIB7rBtSYOVcDUqOD6dcpFy1OwAUy7qtELSf9t+4zRT9LC9e1uxh/bhlZSLmBtS5qNvB5wLtD4dtXlnM5uS9E6n1LU5B+gi6bU5O2QxXR57s7L2rSyedGv8zo8E7g/sSJkRtxdwfttxm9jbURZePRZYjTIx4U+2266bNL0urhSvKB+UefCfooxpv4qSeF5fKfYmy7JtwjGf1Hx++qiPSsd9CeWPfn5ze3PghIqxVwIuam7PAE6tEXvEvny/crzOjp2ymvoAyqy4pwPHAO+vFPviKZ8fCJxRKfZc4NHAhZRkvz/w3pq/96kfvT7Dt/2hpr7GLZTVcO+0fWal8F+lLDUf9hXgSW0FdDPf2fZZg23N8MpGrrfc+8+2/ywJSavbvkL1arTfbvsuSQtVqlPeQIWL5JKmPreizMW+GMAtDiEO6eTYG+vYPlrSQc1r7yxJZy31uybj9ubzbc01ut9TLlhXYfsqSSvbvhM4VlKt+kUj9TrhS9oE+OEgyas0R5hp+5oWY25OeWv7oCnj+Gsx9Ha3Tc0Mmd0pv//5wG8lnWX7TRXCL5D0YErBsDMl3QTUumg7t4n9aUpV0D9S5+39NZSTivdQEpAodZtqDut0dexQivNBWWz1XMrve8NKsU9rjvuDlAVgpgzt1HCbSoHA+ZI+AFxPmYbcmb5ftJ0LPMX2Hc3t1SiLMv6uxZjPB/6eknCHG1DcChxvu/UzgMEFw2Y8fSPbh7R9sXqa/Xg68CDg24PfQcXYM4G1ar2zkbQH8M/Ah2yfIulq27XOsKfuy0zqHvvzKP/gNqKMaa8FHGa7WvOZZj9WB9ZwvYJ1j6C8k1qV8rt/EPAJV5iVNO0+9Tzhz7e99ZRttWrpbG/73LbjTBP7Espqx88BB9v+aYXZSWvZvmWauiq4xQVnkjZvho5GToO0fcGo7S3sxwOAd1PGdbex3fpZ7opy7LVlNtxovR7SoQxl7D4402jOvn9XKfZVKnW6Z7JoAbNXVIj9LkqZ2HOaZP9ISrneNn2JUjFyHuVttaZ8bvNs902Url4fHnGfqTBDCMD2n4A3SXoCpcZKDZ0du6S32v6ApCMYPSvsDW3FptvZcCfafqFGF26rdc1mpL6f4T8K+CL3tCBbALzM9i8rxP4x5W3uPIaKp9lus+tRp5qFKBu5QlXOEbFXAra3Xas43nDsJV4Xap6XDWwvaCl+J8cuaTfbp0rad9T9tltd4dsc9162T2wzzoi469u+vhnSWYztX9Xcn2G9TfjDL4bBsmvbt1aMv9hwUsXYH+CeC4jfBp4AvNH2F5b4jZOJPc92azORlhL7XNu1zqyH455EmRJ5MuUf/G8pF+gfTZkf/kzgkDZniHV47CsD7/MSun61HP9s27XqFQ3HXRn4ju1n1Y69JL1sYg5g+y7gdc3Xf6yZ7BunSXpO5ZgDO9u+hTLEsoBSIrfWH+R5klq7KL4UZ0jaszmjrsaloc6/U6b+fpzyzu4UytqPK4FnVJgO3NWx30mLU42XwZmS3iJpI0lrDz7aDtoc920qJR1WGL09wwcGrf5uB04A/jTY3uYFxKHYt1KmaP2FMm1NJXT7q/AkXWr7cZI+Tant8e2KF6svo/yD+RXlOR8cd+vjmkPP+ULgz1R8zrvW5bFL+jCllMZJLPp31vqFU0mjCiG6xgwplY562wFnsuhxt3ntYsn71POE39mLoUuS3keZGno7pYjZg4HTbG9bIfYKN67Ztulmigzc12eMSDp2xGbXmKAgaQ1PKRI3altLsUddu7A7rF/U64TfNXVbROwhwC2+p+3eWrZ/XSHucbZftrRtLcX+nu1nLm1bC3FHJbyBWomvk2Nv4uww9YLxqG0txV6sKumobS3FPsj24UvbVlOvp2U2ie5NlBKtB0ralFIq+bQKsasXERs1N3nKkG6NM83HTdmnlWl5jFfSGpQCWus2/+gGB70W98zQao3t/duOMZ2uj71xBIuXERm1bWIkPYzSbOd+kp7Iosd9/7biTrEvMDW57zdiWzW9TvjAsZRZE09pbi+gjDO2nvApyf7vgPNs76hScqHtevRdzk1+B/CvlD/AWwabKd2W2m4E8o/AGykJbh73/PHfQrmIWo3ql+rt7NglbU/521pP0nDZjrUoxcTa9GxKct2QsgZh+Lj/tc3AkvYGXgJsIml4NfGalFo+nen1kI6kubZnaag2ecWLlz+1/XcqNfm3tf2XLqdq1iLpP22/Ywn3P86lO1QbsV9v+4gl3L9Ty1MjR5bqtX1AWzGHYlc/dpXSGbOBV1P62g7cSqnU2fZiPyTtuaS1LZL2nfR6gOY61SaUjlpvH7rrVkrVzrZbK06r7wn/x5Q50D+yvU2zEOvLtp9cIfbXKeVS30gZxrkJWNV261M1VYpJvZzFV/l2NntgoNb4ahex1ZSvGPr8QOBrtnduK+Zy7Ftrxy7pEUu6KC/pCNuvbyP20nT8equ+NqLvQzqHUhYebSTpi5RWZFXGW23vMdgHST+gKSJWIzalAcZ5lBrpbfcUXV5V54lXjj2YGdJJqd6laO3Yl2EGVpctRbt8vVWpjjus1wnf9hmS5lEumAo4yHbrtXSaVb4X296y2Y9atcEH1nCdUsj3RpdvOduOfaoWL9X76ZZjLqu+vtW/L7/eFtPrhD80Le2bI7a1xqURxUWSNnYHdWWA4yS9inJx+i9D+9X6grO+av7Jf8+lifZXJZ1GxVK9Ma0uz/Cr62XCX0Gmqq0PXCrpfBZdhdd6X1nKzJgPAgdzz1lG2xUrl1XVuvhTXNPWD27+yX+Ypkqm7b8w9M92BXBNh7G7TLrVi+kNqX7cvbxoK+kg7pmq9n8sOmXr07Y/VmEfnj5qe43hHUm/pMwMqlUKejj2DpR+tn+S9FLKXOzDa620lbQlpZH38LTIKisfJR0GXEy5UFv9D6/LY18SSfvZ/mxLP3sG8B/Aw23vKmkLSuXQo9uItzwkbWn7Z1Vj9jHhDyxtqtp9VTM3+MW2b+sg9sWU6pxbAccBRwP/YHvkP8AJxz6EMk1wC8qF610pPQH2ajt2E7/LejadHbuk9YC3sfg/m9b7EEj6FmW9zcG2nyBpFeBC24+vEHs7ygKzxwKrUdYe/KnL2k29rZYJYPsISU+R9BJJLx981IgtaTtJP5X0R0l3SLpzaEFS2+6k9Nn8lKSPDj4qxV7YnN0+n3JmfzhlQUoNe1Gm4f66Wf36BGD1GoFVljQ/zvZKtlezvZbtNSv+8Xd27JSeE5dTZiQdRhk++mml2Ou61MO/C6CZA3/nkr9lYj4G7E1pLnQ/4JWUfwCd6eUY/oCk44BHUUobDF4EBmq8zf0Y8GLKyt5ZlHnxm1aIC6WB+DcqxZrq1mbV7UuBpzWlFVatFPv2Zix9oXI72L4AABHxSURBVKS1KP1Gq1y3sO1m7UVXpYI7O3ZgHdtHN3VkzgLOklRrZtqfJK1Dc62qOeuudqHc9lWSVnYpl3xss/anM71O+JREu0UX46nQ3YthaSsLJX3V9p4thX8RZdn5AbZ/LWljygXkGuY20yI/TSkz8Efg/EqxoekFYLvW2e2wLo/9r83n65vSEtdRSh7U8CZK74FHSfoRsB7l3U4Nt0lajfJu+gPA9ZQhvc70fQz/JOANtq/vIPbZwLMoS+x/TXkx7FejrMPSDJeauK+SNJNSIfTiijEvozRBuYbKvQCm7MdMKh67pOdRmr5sRBnSWAs4zE0v6QrxV6E87wKutP3XpXzLpOI+gvJOalXgnymLKz9h+6oa8UfuU88T/g+ArSlnOsPz0VufGrkivhgG2lhuLukc209tLlwOv+hqXrgUsA/wSNvvat5dPMx2lTNdddALQNLmtq+QNPL3afuCtmKvCDS6F8HNwCW2b6i9P13re8LvbGrkiqzL+iJtkvRJysW7Z9h+bLMG4wzb1VouSnoqsKntY5vZKw+0PaoRz6TiHeVS+vsHI+52mzNlJL3V9gckHcGIVaU1ajdJ+iZl7cPg+GdTyopsBrzL9nEtxDzR9gslXcLo4676jm5Yr8fwu0zsKt22Rr0YVoTFT60tCFGHDVAoaw+2kXQhgO2bmjHWKpqpkbMowwvHUt7dfYEWa8nYPrD5vGNbMZbg8ubz3A5iD9wFPNb2b+DuefmfBLYFzqZMDZ60g5rPz2vhZ4+llwl/xLDC3XdRr8fprKGv1wBeALTeXHlA0v0ojV+uHHH321oMPbUByirUm7ny12ZW0GDGxnrULR63B/BESh0dbF8nqcqUVEnvBg5tJgjQzNQ53C02Z7F9avN5ouWHl9PMQbJv3ABsZvtGSa2M5Q+uCdZaTLg8epnwbdea972kfZjaCOG/JZ0DvLPt2JJ2Az5EWQyyiaStKW9vd2/27YwWYnbZAGXgo8DXgYdKei9ltsa/VYoNcEczPXPwD6fmjI1VgPMl7Q88jHLxtMqccEmnsvgJ1s2UM/9Pud3+sj9UqVt0UnN7T+Ds5rn/Q4txpzuxHBz3m21f3Wb8kfvU5zH8Lk25iLYS5Yz/n2rM0lGpEPoMYI7vafxycY2xRS2lAUqF+JtTFiCJUszs8qV8yyRjv4Wy1mInSnOMVwBfqrXaW9KzgFMpvReeVmuCgKTDKdMhv9xsehFlZtr9KLOFWhvOay7U70kZNhNwDvDVGlOxVUppXAd8qYn9Yso/2yspf+uz296HxfYpCb8bzUW0wZO/kDJV70O2f14h9k9sb6tFO321mvBXlNkizZDODBZt/NJqxVJJq7sUS0PSTsDOlATwHbfYYWvKPjyNMnb9BeDxlOHDV9i+rkLss20/bdQ2SZfaftx03/u3bPB3NmXbeba3U6XOelP1ckhnBXEaJeEPLpAaeJ6apuK2P9Ji7J9Jegmwskrj9jcAbS/6ehNwIKW/6FSmxebtA5JeDxwC/IayslpN7Lbf2ZwLbDN0cbpKkp/iQ8ALbF8Gd09X/D6weYXY62moFHgzHXbd5r5Wq6Oq23o2d0l6IfCV5vbwgq9OzrST8LvzJEoT85MpiWc3yqyBayvEfj2lNPJfKG+zvwO8u82AHc8WGTgIeMyI6ydtW03SvsBTRs0Lt91a8/gh2w8u2A5iVixv8GbgHJUqraLU1HlNM47e9gXdUSVMHt1yzIF9gMOBT1AS/HnAS5sJE6+rtA+LyJBORySdAexp+9bm9prASbZ36XbP2iXpIuB44ETbv6wc+wfATq7cRLqZe78P8ELKMv9htv2KCvswKBO8ge1dVLlMsKTVKe8mBFzR8oXa4bhzbc8aHrKU9GPbT6kRf0WTM/zubMyib2fvoDQVb52kzYC3sHgT89aHVYDdKRftTpR0F3ACJfm3No4uadDO8WpgTrMYZ3hldZvDZ9g+h3KGO7dWgh3hszRlgpvbP6c897X2Z1PK+oM1gK0k1arFX72ezXQLzQZqLDibThJ+d46jTJP7OuXFsQftv70dOAk4klLHp1apWODuuckfAD7QXD/4d+D9lLHVtgym4f5v87Fa81Hb8ZL+jbL+4cDm+B9j+7QKsde1fWIzPRbbCyVV+d1rmlr81KlK+zLKLLjXUUqYbESZtdOmLheaLVGGdDrUzFj5f83Ns21fWCnuPNtdlekdFO96IeVM/07gBNujLuZOMubKwPts/0ubcZayDydQKlW+3PaWzVjuuba3rhB7DiXRndmsNt4OeL/rNJ65hFJ//0KXJiQzgM/Y3q3t2LGonOF3qJmK2EXxqlMlvYayCKlqE3NJP6GUFDiJMmukyuIT23dONyW0okfZfpGkvZt9ul2DaVnt67JMcGe1+FVaah4KPIJFhy9bj68OO31NJwm/n/ZtPg+f7dZqYr6v7SsqxBllvkp7x5NYtHF8jVkyAHc0Z/WDlbaPokIjc0krURLO0+mgTDDd1uI/mjKUM4/Kw5eUTl8nAM8FXk35u/tt5X1YRIZ0ogpJL7X9haELqIto+8Jpsw/Hjg5dZZaMKOPJB1DO+M6grP7cz/acCvHPtb1923GWYT9mUrcW/2KLn2oZDJ1OmSF0Vo1htOnkDL9HJD3D9vdHzQWH1s90BzMjRtUxqrHMfWXgYtv/1XasUZoaOgdRVtluRznLPsj27yrtwhmS9gS+VqOswFSSNmBoWEXS02yf3WK8wfDdDyR9EPgaiw5f1hhK7bLT10g5w+8RSYfZPqTjM90dbP9oadtaiv2DLhd+Sfo48Fl30OKwKeT1AEoZjz9Tt/HM+ykX6C9jqHe0W2w0pNH1/wdcYxxdHXf6GrlPSfj9I2kTT2m6MWpbS7EXa64yaltLsd9L6Sx2AouO4deq43MZpfHGr+iwxWFtkq4EthrUE6oYdyVgL9sn1ow7FH/tGhMhlkeGdPrpq8DUBPsVWqxLL2l74CmUuirD4/hr0e4c/GGD1ZXvGtpWpY5PY9dKcRYj6SvAMcC3bdfsAQBlwduqVLhAPayZGfQ6oJOED/xE0nzKgrdvdTGUNlUSfo+olAZ+HPCgKeP4azE0bawlqwEPpLzmhsfxb6HS9MCO6/h03RDjSGB/4AhJJ1GGllqdLTW04vQ2ygyp77HoOHqNFadnqpSlnvqursaZ92bAsyhlsI9o1mF81hUq4k4nQzo9Iun5wN9TyhsMjyPeChxvu+2KmUh6RJeJr7l49jgWnRf9rum/475F0oOAvSklFq6lTJX8QhtTNJuCcdNyhU5YKq1ER4Su20pU0o6U0tQPAC4C3m773Jr7AEn4vSRp+y5ebE3sMykLrv7Q3H4I5Z/NsyvEPhK4P7AjpazEXsD5tg9oO/aKQNI6lKmhL6XMGPki8FTg8W65GUdTz2Zzyhn/lbZbLYvcxFyJUiCu9QkB08Rfh/Jcv5zS8OVoyonW1pRCiZtU36ck/P5pVgC+isWLp9WYpXN305UlbWsp9sW2txr6/EDKNMWd247dNUlfoyTc44Bjbf966L65tmdN+83jx34O8ClguDzyP9r+Vlsxh2J3tv5A0s8pz/cxtv9vyn1vs/3+2vuUMfx+OpkyXey71F99eJcWbYbxCOo1g7i9+XybpIcDv6cknz74DGXB1w7ALJX+yZ+0/ec2k33jI8CObloqNiuMvwm0nvDpdv3BPpQ+zntIGj6x2qqLZA9J+H11f9tv6yj2wZRSwYPmG0+jdMKq4bRmif8HKTWMTEmEfbAf5QL5R5vbe1POPl9QIfYNXrR/7tWUejo1vIkybn6npNupuP6AMmb/FuBnQO2ZUSNlSKeHJL0H+LHt0zuKvy73rDY9t+Jq0+F9WB1Yw/bNtWN3QSN6qI7a1lLsT1JW2Z5I+Sf7Akoj7x9B1VpGVUk6x/ZTu96PYUn4PTS06vIvlOXfNc96FltmTwne5jL7LktKrBAkfRY40vZ5ze1tKYXsXlMh9qiV3QOtr/CWtDvlnSTAHNfpP4CkZ1LeSU2djtrZ6y0JP6oaWmZ/Kfe8zW17mf2htg9tEs/wC37wj671i9Vdk3Q5pVLmoLPYxsDllN/BfXa1r6T3UXpHf7HZtDcwz/bbK8T+AuVC+dTXemevtyT8nmqmQ27KovPRWzvLHopbfZm9pDdTEr2GPtN8XaVSZ9eai+PTanNthEprwfdQLpp/m9IM5Y22v9BWzKHYFwNbD1YXN0X0LqzxD07SJbYf33ac5ZGLtj0k6ZXAQZTKffMp4+nnUqfEQBfL7B/YfH4M5WzvZErS3w1o/Z/ciqDjVb47236rpD2ABZQx/B9QLmrW8GBgsLL2QZViApwnaQvbl1WMuURJ+P10ECXxnWd7x6bkwmGVYldfZm/7MABJZwDb2L61uX0opRlKtGvV5vNzgC/bvlHVGn3xn8CFTfVMUcby31Ep9lOBfZvVvn9hBSiWl4TfT3+2/WdJSFrd9hWSHlMp9iksWtahpo2B4RWed1AWn0W7TpV0BWVI5zXNwr8/1whs+8sq/Xz/jpJw3za86Kxlu1SKs8yS8PtpQTMf/RuU4lI3UZbat87255pl9ps1m2q22jsOOF/S1ynj93sArddz6Tvbb28u1t/i0lv4NuD5FXdhvebzysBTJFWZKdPxMNpIuWjbc5KeThnX/Hal+iazKUn2GsoZ10aU6YFVxtKbTkj/r7l5tu0La8TtM0lzKaWZv2z7psqxjwG2YgWaKdOlJPyeaQpKXWx7y47izwNeYvvK5vZmlETQWi3+6JakR1NKM78ImEupD39GjVIHki6zvUXbcf5WrNT1DkRdzfS0iyRt3NEurDpI9s3+/Jx7LurFfZDtq2wfTBnG+xLlbP9/JR0mae2Ww58rKQm/kTH8flofuFTS+SzaFKK1xU9D5ko6mjKeDqXA1LwKcaNDkraiNALZldJxbVCa+fuUcsFt+Rwl6f+aFWSmTJcypNNDzbj9YmyfNWr7hGOvDryW8scuyjz4T9RciBV1NcN4f6AUqvva8O9a0tdsjyx5MaHYV1EKqF3CUAGzFfGCag1J+BHRiqHexasAC6feX2OFs6Tv267Vs3iFlyGdHpK0HXAE8FhKr9mVgT+1WTxN0om2XyjpEhavf2/KSsj/tn1yW/sQ1Q16Fw9WOA/WX9Rc4XyFpC8Bp7KCFDDrUs7we6iZJvdiyirTWZQWbJva/tcWY65v+/ol1HRZF/ii7c3b2ofoRrPCec+hFc5rUlr8tb4waZpKnb2dlpkz/J6yfZWklW3fCRwrqdUG5ravbz5PN3b6K0n7tLkP0ZnOVjjb3r9GnL8VSfj9dFuz2nV+U8nwekp9/NY1NenfDzyUctH27lr8tjNb576psxXOzTqPTwIzbG/ZzBba3fZ7asRf0WRIp4eaYZUbKPPf/5my0vYTU9rQtRX7KmA325e3HStWHF2tcG5aaf4L8CnbT2y2/ayrhYddyxl+Dw0Nq9xOvSqZA79Jsu8f2xdQ+gjXdn/b50+pzrnYjKG+SMLvoaZc62Jv7Ww/skL4uZJOoBRu6/2siWjd7yQ9iub1LmkvyhBmLyXh99Osoa/XoDSkaHuJ+8BalJr4Ow9tM5CEH214LXAUsLmk/wP+h7K6u5cyhh8ASDrH9lO73o+ISWpWdu9FmRW0NnALZZLAu7rcr67kDL+HmgtoAytRzvjXnObhk469IWXR1w6UM/tzgINsL6gRP3rnZEpZhwuo1PNhRZYz/B5q2r0NfvELKbXpP9RUrmw79pmUiomD4mkvBfaxvVPbsaN/+jwjZ5Qk/B6S9GZKwh9MXVjkRdBmjRNJ821vvbRtEZMg6SjgCNuXdL0vK4IM6fTTkyi1TU6mJP1BbZNrK8T+naSXAl9ubu8N/L5C3OiRoZpNqwD7S7qalEfOGX4fdVzbZGPgY8D2lD/IHwNvsP2/bceO/lhCzSagv+WRc4bfT53VNgHeTelhexNA0/HoQ5TmGBET0deEvjRJ+P3UWW0TYKvhRta2b5T0xEqxI3otCb+HbL9X0re4p7bJ/rVqmwArSXrIlDP8vA4jKsgfWk91WNvkw8CPJX2F8u7ihcB7O9iPiN7JRduoTtIWwDMoMya+Z/uyjncpoheS8CMiemKlrncgIiLqSMKPiOiJJPyIiJ5Iwo+I6Ikk/IiInvj/MwsPeDIkj6gAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.bar(['Intercept']+predictors, list(model.intercept_.flat) + list(model.coef_.flat))\n",
    "plt.xticks(rotation=90)\n",
    "plt.yticks(ticks=np.arange(-3, 2.5, 0.5))\n",
    "plt.grid(axis='y')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The coefficients are similar (not in absolute value, of course) to those found in Annie's paper:\n",
    "- sick contact is the strongest (positive) predictor\n",
    "- conjunctivitis is the next strongest (its coefficient is 2/3 of the quarantine coef in the paper)\n",
    "- diarrhea and referral(4-9) are next (their coefficient is 1/2 of the quarantine coef in the paper)\n",
    "- Then come dysphagia and hemorrhage (their coefficient is 1/3 of the quarantine coef in the paper)\n",
    "- Pyrexia is the weakest positive predictor (its coefficient is 1/6 of the quarantine coef in the paper)\n",
    "- Finally, myalgia is a mild negative predictor (its coefficient is -1/3 of the quarantine coef in the paper)\n",
    "\n",
    "Hemorrhage and myalgia are found somewhat weaker here, but the similarity is good overall."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean of the available values for `quarantine`: 0.24682395644283123\n",
      "Coef for quarantine_missing VS coef for quarantine: 0.194432945917306\n"
     ]
    }
   ],
   "source": [
    "print(\"mean of the available values for `quarantine`:\", df['quarantine'].sum() / (len(df) - df['quarantine_missing'].sum()))\n",
    "print(\"Coef for quarantine_missing VS coef for quarantine:\", model.coef_.flat[1] / model.coef_.flat[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}