sujitpal/02-symptom-checker-experiments.ipynb

## 02-symptom-checker-experiments.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Symptom Checker and other Experiments\n",
    "\n",
    "The data for this notebook comes from the Medium article [Profiling presenting symptoms of patients screened for SARS-Cov-2](https://medium.com/@nigam/an-ehr-derived-summary-of-the-presenting-symptoms-of-patients-screened-for-sars-cov-2-910ceb1b22b9) by Prof. Nigam Shah.\n",
    "\n",
    "The dataset represents 895 patients admitted to Stanford Health Care Emergency with possible symptoms of COVID-19, of which 64 (7.2%) tested positive. Text of the patient notes for these patients were mined for observation terms, and their frequencies calculated across the corpus. Frequencies and conditional probabilities for the top 50 observations were provided as an Excel spreadsheet. Please see the article for details on how this was done.\n",
    "\n",
    "This notebook uses a manually converted CSV (File :: Save As :: CSV) from the Data tab of the provided Excel worksheet, followed by some edits to remove extraneous rows and columns. It then uses the data to run various experiments and analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import collections\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import os\n",
    "import pandas as pd\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "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>Clinical observation</th>\n",
       "      <th>Count (observation)</th>\n",
       "      <th>Count (observation &amp; +ve)</th>\n",
       "      <th>Count (observation &amp; -ve)</th>\n",
       "      <th>P(observation)</th>\n",
       "      <th>P(observation|+ve)</th>\n",
       "      <th>P(observation|-ve)</th>\n",
       "      <th>P(+ve|observation)</th>\n",
       "      <th>P(-ve|observation)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>cough</td>\n",
       "      <td>577</td>\n",
       "      <td>51</td>\n",
       "      <td>526</td>\n",
       "      <td>0.645</td>\n",
       "      <td>0.797</td>\n",
       "      <td>0.633</td>\n",
       "      <td>0.088</td>\n",
       "      <td>0.912</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>dyspnea</td>\n",
       "      <td>526</td>\n",
       "      <td>41</td>\n",
       "      <td>485</td>\n",
       "      <td>0.588</td>\n",
       "      <td>0.641</td>\n",
       "      <td>0.584</td>\n",
       "      <td>0.078</td>\n",
       "      <td>0.922</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>febrile</td>\n",
       "      <td>396</td>\n",
       "      <td>44</td>\n",
       "      <td>352</td>\n",
       "      <td>0.442</td>\n",
       "      <td>0.688</td>\n",
       "      <td>0.424</td>\n",
       "      <td>0.111</td>\n",
       "      <td>0.889</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>sore throat</td>\n",
       "      <td>244</td>\n",
       "      <td>13</td>\n",
       "      <td>231</td>\n",
       "      <td>0.273</td>\n",
       "      <td>0.203</td>\n",
       "      <td>0.278</td>\n",
       "      <td>0.053</td>\n",
       "      <td>0.947</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>chest pain</td>\n",
       "      <td>129</td>\n",
       "      <td>11</td>\n",
       "      <td>118</td>\n",
       "      <td>0.144</td>\n",
       "      <td>0.172</td>\n",
       "      <td>0.142</td>\n",
       "      <td>0.085</td>\n",
       "      <td>0.915</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  Clinical observation  Count (observation)  Count (observation & +ve)  \\\n",
       "0                cough                  577                         51   \n",
       "1              dyspnea                  526                         41   \n",
       "2              febrile                  396                         44   \n",
       "3          sore throat                  244                         13   \n",
       "4           chest pain                  129                         11   \n",
       "\n",
       "   Count (observation & -ve)  P(observation)  P(observation|+ve)  \\\n",
       "0                        526           0.645               0.797   \n",
       "1                        485           0.588               0.641   \n",
       "2                        352           0.442               0.688   \n",
       "3                        231           0.273               0.203   \n",
       "4                        118           0.144               0.172   \n",
       "\n",
       "   P(observation|-ve)  P(+ve|observation)  P(-ve|observation)  \n",
       "0               0.633               0.088               0.912  \n",
       "1               0.584               0.078               0.922  \n",
       "2               0.424               0.111               0.889  \n",
       "3               0.278               0.053               0.947  \n",
       "4               0.142               0.085               0.915  "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nigam_df = pd.read_csv(\"../data/nigam_data.csv\")\n",
    "nigam_df.drop(columns=[\"Unnamed: 4\", \"Unnamed: 8\"], inplace=True)\n",
    "nigam_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.02941785635414167"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def compute_prob_positive(prevalence, symptoms, nigam_df):\n",
    "    # populate the array with S(k)=False probabilities (1 - P(observation|+ve))\n",
    "    ps_s_given_d = [(1 - x) for x in nigam_df[\"P(observation|+ve)\"].values]\n",
    "    # find symptom indexes for symptoms provided\n",
    "    all_symptoms = nigam_df[\"Clinical observation\"].values.tolist()\n",
    "    sym_idxs = [all_symptoms.index(x) for x in symptoms]\n",
    "    # update these positions with S(k)=True probabilities, i.e., 1 - P(S(k)=False)\n",
    "    for i in sym_idxs:\n",
    "        ps_s_given_d[i] = 1 - ps_s_given_d[i]\n",
    "    p_d_given_s = np.prod(ps_s_given_d) * prevalence\n",
    "    # we have to normalize, so calculate p_not_d_given_s\n",
    "    ps_s_given_not_d = [(1 - x) for x in nigam_df[\"P(observation|-ve)\"].values]\n",
    "    for i in sym_idxs:\n",
    "        ps_s_given_not_d[i] = 1 - ps_s_given_not_d[i]\n",
    "    p_not_d_given_s = np.prod(ps_s_given_not_d) * (1 - prevalence)\n",
    "    norm_denom = p_d_given_s + p_not_d_given_s\n",
    "    return p_d_given_s / norm_denom\n",
    "\n",
    "\n",
    "compute_prob_positive(0.072, [\"cough\"], nigam_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### How does P(D=+|S) change as we add symptoms?\n",
    "\n",
    "The observations are listed in decreasing order of frequency, but the ones towards the beginning are more indicative of COVID-19 than the ones later in the list. That is why we see the probability of the patient having COVID-19 increase as more symptoms are added initially, but falls off as symptoms that are more indicative of other related diseases are added."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3xUZdbA8d9J74UkFElCLwJCwIAKgoCKYsO6ghV1dbHrrrvruu/a9tXdd3ctqLguuohrwbIrVlyxUESkE0poCTWhpJJCenneP2YCISSQgUxm7p3z/Xzmw8ydOzPnKsyZp5znEWMMSimlfJefpwNQSinlWZoIlFLKx2kiUEopH6eJQCmlfJwmAqWU8nEBng7gZMTHx5vu3bt7OgyllLKU1atX5xtjEpoet2Qi6N69O6tWrfJ0GEopZSkisru549o1pJRSPk4TgVJK+ThNBEop5eM0ESillI/TRKCUUj5OE4FSSvk4TQRKKeXjNBEoy/shI4/12UWeDkMpy9JEoCztu8053DprBTe9sZy9RRWeDkcpS9JEoCxr495i7p+zlr6dIqk38PD7adTW1Xs6LKUsRxOBsqR9RRXcPnslMaGB/Ov2EfzvlYNYsauQVxZkejo0pSxHE4GynNLKGm6fvZLy6jpm3TacjlEhXDm0K1cP7cpL32Wwclehp0NUylI0EShLqa2r57731pKRe4hXbxxG/85Rh597+spBJHcI48E5aykur/FglEpZiyYCZRnGGB7/LJ1F2/J45spBjOl79Gq6EcEBTJ88lNzSKh79eD3GGA9FqpS1aCJQljFz8Q7eW76Hu8f2YvKI5GbPGZIUw68v6sdXGw/w/sqsdo5QKWvSRKAsYfmOAv783y1cOrgLv57Q77jn3jm6J6P7xPPU5+lk5JS2U4RKWZcmAuX1iitq+OWH6+geF85frx2Mn58c93w/P+G564YQFhTA1DdXMndttk4rVeo4NBEor/fEpxs5UFLJC9enEBbUuk31OkaF8PotZxIZEsDDH6zj/OcX8eHKLGo0ISh1DE0Eyqt9mraXT9L28eD5fUhJinHptWd268C8B0bzj5vPJCI4gN/8Zz3j/raQd5fvpqq2zk0RK2U9mgiU19pbVMH/fLKRYckx3DO210m9h5+fcNHAznxx/7nMmppKXEQwv5+7kfOfW8T+4tYvSbG7oIzy6tqTikEpb6eJQHml+nrDrz5Mo77e8ML1KQT4n9pfVRFhfP9OfHLPSGbfNpzCsmoemLO2VWMHq3YVcsHzi/j7wu2nFINS3koTgfJKbyzZwbIdhTxx+UC6xYW32fuKCGP7deTZq85g5a6DPP/NtuOev7+4gmnvrKGmzrAjr6zN4lDKm2giUF5n074S/vr1Vi4a2InrUhPd8hlXDu3K5OFJvLpwO4u25TV7TmVNHdPeXk1FdS29EsLJPljulliU8jRNBMqrVNbU8dAHa4kJC+JPVw9G5PhTRU/Fk1cMpH/nSB7+II0DxZVHPWeM4bG5G1iXXcwL16cwokcc2Qd1mWtlT5oIlFeZvymHbTmHePaqM+gQHuTWzwoJ9OeVG4ZRWVPHA+8fPV4w68ddfLxmLw9f0JcJAzuTGBtKQVm1DhgrW9JEoLzK/PQDxEcEMb5/x3b5vN4dI3jmqkGs2FnIi99mAPBjZj7PztvMhAGduH98bwASY0MB2KutAmVDravOUaodVNXWsXBrHpcP6YL/CaqH29JVQxP5aXsBMxZmclpMKH/5egu9EsJ5/vqUw1XMibFhAGQfrKBPp8h2i02p9qAtAuU1ftpewKGqWiYM6Nzun/3UFYPo0zGCx+ZuoL7eMPPmVCKCj/xOSnK2CHTAWNmRJgLlNb5OzyE8yJ9zesW1+2eHBvnz6o3DGJIUw6s3nkn3+KOnrMZHBBMU4KcDxsqWtGtIeYX6esM3m3IY268jIYH+Homhd8dIPr13VLPP+fkJiTGhmgiULWmLQHmFtVlF5B+qYsLATp4OpUWJHcK0a0jZkiYC5RXmbzpAgJ+j6tdbJcZqi0DZkyYC5XHGGOan53BOrziiQwM9HU6LtJZA2ZUmAuVx2/MOsTO/jAkDvLdbCI5MIdVaAmU3mgiUx32dngPABV6fCBqmkGoiUPbi1kQgIheLyFYRyRSRR5t5PlpEPheRdSKSLiK3uTMe5Z3mb8phSGI0XaJDPR3KcSVqLYGyKbclAhHxB2YAE4EBwBQRGdDktHuBTcaYIcBY4DkRce8CM8qrHCiuZF1WERMGtn8RmasSIoIJ1loCZUPubBGMADKNMTuMMdXA+8CkJucYIFIcS0xGAIWAjsT5kG82O7qFvH18ABx7GXTVmUPKhtyZCLoCWY0eZzuPNfYKcDqwD9gAPGiMaXbLKBG5S0RWiciqvLzm149X1jM//QA94sPp3THC06G0SmKs1hIo+3FnImhu1TDT5PFFQBpwGpACvCIiUc29mTFmpjEm1RiTmpCQ0LaRKo8orqjhp+0FTBjQya37DrQlrSVQduTORJANJDV6nIjjl39jtwEfG4dMYCfQ340xKS+ycGsutfXGq6uJm9JaAmVH7kwEK4E+ItLDOQA8GfisyTl7gPMBRKQT0A/Y4caYlBeZn55DfEQwQ5NiPR1Kq2ktgbIjtyUCY0wtcB/wNbAZ+NAYky4i00RkmvO0PwIjRWQD8B3wW2NMvrtiUt6jsqaOhVtzuXBAp8Nr/luB1hIoO3Lr6qPGmHnAvCbHXmt0fx8wwZ0xKO/00/YCyqrrLNUtBFpLoOxJK4uVR6zNKkIEzunZ/nsPnAqtJVB2pIlAeURWYTldokI8tvfAydJaAmVHmgiUR2QVlpPUIczTYZwUrSVQdqOJQHlE1kErJwJtESh70USg2l1lTR05JVUkxVo3EWgtgbITTQSq3TX8mk6O8+7VRluitQTKbjQRqHaX5exft3KLALSWQNmHJgLV7rIKnYnAwmMEoLUEyj40Eah2l1VYTnCAHwkRwZ4O5aRoLYGyG00Eqt1lFVaQGBtqqaUlGtNaAmU3mghUu9tTWE6yRbuFGmgtgbITTQSq3Vm5hqCB1hIoO9FEoNpVcXkNpZW1lp0x1EBrCZSdaCJQ7WqPxWcMNdBaAmUnmghUuzpcQ9DBmsVkDbSWQNmJJgLVruzTItBaAmUfmghUu8oqLCcmLJCokEBPh3JKtJZA2YkmAtWusg5WWH6gGLSWQNmLJgLVrhz7EFh7fKCB1hIou9BEoNpNfb1h78EKy48PNNBaAmUXmghUu8kpraS6rt4WXUOgtQTKPjQRqHazp8AeM4YaaC2BsgtNBKrdZDVsSGObRKC1BMoeNBGodpNVWI4InBYT4ulQ2oTWEii70ESg2k1WYTmdo0IIDvD3dChtQmsJlF1oIlDtxg6rjjamtQTKLjQRqHaTVWiPYrLGTosOZV+xJgJlbZoIVLuorKnjQEmlbQaKG8RHBJF/qMrTYSh1SgJae6KIdARGAacBFcBGYJUxpt5NsSkb2Vvk+NVsl6riBvERwRQcqvZ0GEqdkhMmAhEZBzwKdADWArlACHAl0EtE/g08Z4wpcWegytqybLLqaFPxkcGUV9dRXl1LWFCrf1cp5VVa8zf3EuBOY8yepk+ISABwGXAh8J82jk3ZSEMisFvXUFx4EAD5pdUkx2kiUNZ0wr+5xphfH+e5WuCTNo1I2VLWwQqCAvxIiAj2dChtKj7ScT15h6pIjrNXklO+o9WDxSLyoIhEicM/RWSNiExwZ3DKPvYUlJMYG4qfn3g6lDbVkNgKdMBYWZgrs4Zud44DTAASgNuAP7slKmU7WQfLbdctBI7BYoB8HTBWFuZKImj4KXcJ8KYxZl2jY82/QORiEdkqIpki8mgL54wVkTQRSReRRS7Eoywkq7DcdjUEAB0axgi0RaAszJXRrdUiMh/oAfxORCKBFqeOiog/MAPHQHI2sFJEPjPGbGp0TgzwKnCxMWaPc4qqspni8hpKKmttN3UUICjAj+jQQE0EytJcSQR3ACnADmNMuYjE4egeaskIINMYswNARN4HJgGbGp1zA/Bxw4wkY0yuK8Era8g6aM8ZQw20qExZXasTgTGmXkRqgTHOaaMN1rfwkq5AVqPH2cBZTc7pCwSKyEIgEphujPlXc28mIncBdwEkJye3NmzlBRqmjibasGsIHOMEOkagrMyVyuJZwGAgnSNdQgb4uKWXNHPMNPP5ZwLnA6HATyKyzBiz7ZgXGjMTmAmQmpra9H2UF9tj02KyBvGRwWzer/WUyrpc6Ro62xgzwIXzs4GkRo8TgX3NnJNvjCkDykRkMTAEOCYRKOvKOlhOdGgg0aGBng7FLeLDg8gv1a4hZV2uzBr6SURcSQQrgT4i0kNEgoDJwGdNzvkUGC0iASIShqPraLMLn6G8xO6CMorKm+8eySqssOVAcYP4iGBKKmupqq3zdChKnRRXWgRv4UgGB4AqHF0/xhgzuLmTjTG1InIf8DXgD8wyxqSLyDTn868ZYzaLyH9xjDPUA28YYzaewvUoD5n65krq6g0f/OJsukQf/aWfVVhOv86RHorM/RqqiwvLqo+5dqWswJVEMAu4GdjAcaaNNmaMmQfMa3LstSaP/wr81YU4lJeprzdkHyynps5ww+vL+eCus+kYFdLouQouHNDJw1G6z+GislJNBMqaXOka2mOM+cwYs9MYs7vh5rbIlGUcLK+mps4wKeU0ckoqmfL6MvKcfea5pVVU19WTaNOBYoC4CC0qU9bmSiLYIiLvicgUEbm64ea2yJRl5JQ4vgAvHtiZN6cOZ19RJTe+sYyCQ1WHZwzZtYYAjqw3lKeJQFmUK4kgFMfYwATgcuftMncEpawlp7QSgI5RIZzVM45/3prK7oJybvrnCjbsLQYgKda+XSZH1hvSRKCsyZWCsuNVESsfllviTATOQdORveN5/ZZUfv7WKp6dtxkR6GrjRBAa5E94kL/uVKYsy5VlqBNFZK6I5IpIjoj8R0QS3RmcsoZcZ9dQx6gjew2M6ZvAazcPw0+gc1QIwQH+ngqvXcRFBGuLQFmWK7OG3gTeA65zPr7JeezCtg5KWUtOaSWxYYHHfNmP79+Jf91+FsUVNR6KrP3oekPKylxJBAnGmDcbPZ4tIg+1dUDKenJKqujknC7a1Dm94to5Gs+Ijwhmd0G5p8NQ6qS4MlicLyI3iYi/83YTUOCuwJR15JZUHq4b8FXxkdo1pKzLpR3KgJ8BB4D9wLUcfxlq5SNySqroFGmvvYhdFR8eRGF5NXX1uh6ish5XuoaSjDFXND4gIqOAPW0bkrKSunpD3qGqowaKfVF8ZDDGOJaZSPDxpKisx5UWwcutPKZ8SGGZ41dwS2MEvkJrCZSVnbBFICLnACOBBBH5ZaOnonAsJqd8WM7hGgJNBKCJQFlTa7qGgoAI57mNl5AswTFOoHxYrrOquJOvdw051xvSojJlRSdMBMaYRcAiEZltjNktIlGOw6bU/eEpb9ewzpCvdw3FaYtAWZgrYwQJIrIBx94BG0RknYic6aa4lEU0dA01dI34qqiQAIL8/XThOWVJru5HcI8x5gcAETkXR2VxsxvTKN+QW1pFXHgQQQGu/KawHxFxVBeXateQsh5X/vWWNiQBAGPMEkC7h3ycFpMdoUVlyqpcaRGsEJF/AHMAA1wPLBSRYQDGmDVuiE95OcfyEr7dLdQgLjxIu4aUJbmSCFKcfz7R5PhIHIlhfJtEpCwlp6SSAV2iPB2GV4iPCGbzfm0kK+txZT+Cce4MRFlPXb0h/5C2CBrERwZTUFaFMQYR8XQ4SrVaqxOBiMQAtwDdG7/OGPNA24elrKDgUBX1BhJ0jABwtAhq6gzFFTXEhAV5OhylWs2VrqF5wDJgA1DvnnCUlRyuIdC1dYAjRWX5h6o1EShLcSURhBhjfnni05SvaKgh8PVisgaNl5no3THCw9Eo1XquTB99W0TuFJEuItKh4ea2yJTXyynVRNCYrjekrMqVFkE18Ffg9zhmCeH8s2dbB6WsIaekCpEjXSK+7nDXUKkmAmUtriSCXwK9jTH57gpGWUteaSVx4cEE+Pt2VXGD2LAg/MQxRqCUlbjyLzgd0E1Z1WFaTHY0Pz+hQ7hjCqlSVuJKi6AOSBORBcDhv+k6fdR35ZRU6vhAE/ERQeTpekPKYlxJBJ84b0oBjhbB4MRoT4fhVRJ0vSFlQa4kggJgnjFGawgUNXX1FJRV+fzOZE3FRwSzM7/M02Eo5RJXxggmAxki8hcROd1dASlryD9UhTH4/Kb1TcWFB+kuZcpyWp0IjDE3AUOB7cCbIvKTiNwlIpEneKmyodzDVcXaImgsPjKYipo6yqpqPR2KUq3m0rw/Y0wJ8B/gfaALcBWwRkTud0NsyotpVXHztKhMWVGrE4GIXC4ic4HvgUBghDFmIjAEeMRN8SkvlVPasFexdg01dmS9IU0EyjpcGSy+DnjBGLO48UFjTLmI3N62YSlvl1tSiZ8c2bRdOTS0CHQKqbISV7qG7gaWAIhIXxG5QkQCAYwx3zX3AhG5WES2ikimiDza0huLyHARqRORa10JXnlOTkklCZHB+PvpuvuNNSQCLSpTVuJKIlgMhIhIV+A74DZgdksni4g/MAOYCAwApojIgBbO+z/gaxdiUR6WW6pTR5sTd3i9IW0RKOtwJRGIMaYcuBp42RhzFY4v+JaMADKNMTuMMdU4BpgnNXPe/TgGoHNdiEV5mC4v0bxAfz9iwgJ1jEBZikuJQETOAW4EvnQeO94YQ1cgq9HjbOexxm/YFcfMo9da8eF3icgqEVmVl5fnQtjKHXJLKumoM4aaFRcepIlAWYorieAh4HfAXGNMuoj0BBYc5/zmOo9Nk8cvAr81xtSd6MONMTONManGmNSEhIRWB63aXnVtPQVl1VpD0IL4CF1mQlmLK5vXLwIWNXq8AzjegnPZQFKjx4nAvibnpALvOzf6jgcuEZFaY4yuaeTF8g7p1NHjiY8MZvO+Ek+HoVSrnTARiMgtrXyvNGPM+kaPVwJ9RKQHsBfHEhU3NH6BMaZHo8+ZDXyhScD75TqLyXR5ieYlRASzWFsEykJa0yLoceJTANjV+IExplZE7sMxG8gfmOXsUprmfP6E4wLKOzVsWq+zhpoXHxFEaWUtlTV1hAT6ezocpU7ohInAGPPUyb65MWYeMK/JsWYTgDFm6sl+jmpfubpX8XHFHa4lqKZrTKiHo1HqxFzeY1BEBrojEGUdOSWV+PsJceG6V3FzDheVafeQsoiT2Wz27TaPQllKTkkVHSOD8dOq4mbpekPKak4mEei/fh+XW1qlNQTHcXgFUq0uVhbRqumjIvIEjhoAATqJyOMNzxljnnZTbMpL5ZZUktQhzNNheK3DC89pi0BZRGvrCHY1ul8D7G77UJRV5JRUkto91tNheK3QIH/Cg/y1a0hZRqsSgTHmrYb7IvJg48fKt1TV1nGwvEarik8gPjJYt6xUlqFjBMolh7eo1DGC49JlJpSVnEwiOL/No1Bepbiihufnb6WksuaY5xpqCLSq+Ph04TllJSdMBOJcCKiBMabwROcoa/tqw35e+j6TX324jvr6o9cJzNWq4laJjwzmQHElFdUnXE9RKY9rTYtggYjcLyLJjQ+KSJCIjBeRt4Bb3ROe8oS0rCJE4JtNOfx90fajnjuyab22CI7nooGdKa2q5f45a6itq/d0OEodV2sSwcVAHTBHRPaJyCYR2QFkAFNw7GM8240xqnaWllXEmD4JXD7kNP42fyuLtx3Z/yGntIpAfyE2TKuKj+e8vgk8fcVAvt2cy2NzN2BM0xXYlfIerVlrqBJ4FXjVuUdxPFBhjClyd3Cq/ZVV1bItp5SLBnbmF+f1ZNuBUh54fy2f33cuSR3CyCmppGNkiFYVt8LN53Qnr7SKl77PJCEymF9f1N/TISnVrNaMEYSIyEMi8gqOfYrzNAnY1/rsYuoNpCTHEBYUwGs3n0ldneHud1dTWVNHbkmVDhS74OEL+zJlRBIzFmxn9o87PR2OUs1qTdfQWzg2kNkAXAI859aIlEetzToIQEpiDAA94sN54foUNu4t4Q+fbCSnpFJrCFwgIvxx0iAmDOjEU19s4vN1TfdmUsrzWpMIBhhjbjLG/AO4Fhjt5piUB6XtKaJ7XBixjVYWvWBAJ+4f35uPVmeTkXtIWwQuCvD346UpQxnerQO//DCNJRn5ng5JqaO0JhEcnkxujKl1YyzKw4wxpGUVkZIUc8xzD13QlzF9HXtFazGZ60IC/Xn91lR6JUTwi7dXHa7HUMobtCYRDBGREuetFBjccF9EdGNWG9lfXEluaVWzicDfT3hpcgoXDezEmD4JHojO+qJDA3n+ZymUVdexYEuup8NR6rDWzBrSvfZ8RFqWYw5ASnLzC8rFhAXxj5tT2zMk2zm9SySdooJZvC2f64cnn/gFSrWDk1liQtlUWlYRQf5+nN4l0tOh2JaIMLpPAksy86mr19oC5R00EajD0vYUMeC0KIIDtBHoTmP6JlBcUcO6bJ2FrbyDJgIFQG1dPRv2Fjc7PqDa1uje8YhwVMW2Up6kiUABsDWnlIqaOoYmayJwt9jwIAZ3jdZEoLyGJgIFNBoo1hZBuxjTN4G0rCKKK45d6lup9qaJQAGO8YEO4UEk617E7WJM3wTqDSzN1OIy5XmaCBTgaBEMSYxGt5ZoHylJMUQGB7A4Q7uHlOdpIlCUVtaQmXeIlCTdkL69BPr7MbJ3HIu35esS1crjNBEo1mcXY5wrjqr2M6ZvAnuLKtieV+bpUJSP00SgjgwUJ2oiaE8NS3Xo7CHlaZoIFGv3FNEzIZzosEBPh+JTkjqE0TM+XMcJlMdpIvBxx1txVLnf6D7xLNtRQGWNbnKvPEcTgY/bW1RB/qEqhmoi8IgxfROorKln1a6Dng7F0owxzFy8nQueX0T2wXJPh2M5mgh83JFCMp0x5Aln94wj0F+0e+gU1NTV89jcDTw7bwuZuYd4+vNNng7JcjQR+Li0PUUEB/jRX1cc9Yjw4ABSu3XQAeOTVFxRw21vrmTOiizuG9ebX1/Uj/mbcnS/BxdpIvBxaVlFDOoaTaC//lXwlDF9E9hyoJScEt21zBVZheVc8/elLN9ZwN+uG8IjF/XjztE96ZkQzpOfp+u4iwvc+q9fRC4Wka0ikikijzbz/I0ist55WyoiQ9wZjzpaja446hXG9I0HdBqpK1bvPsiVM34kr7SKf91+FteemQhAUIAfT18xiN0F5cxcvMPDUVqH2xKBiPgDM4CJwABgiogMaHLaTuA8Y8xg4I/ATHfFo4619UApVbX1mgg87PTOUcRHBLNYN7VvlU/W7mXK68uIDAlg7j0jOadX3FHPn9snnksHd2HGgkyyCnXguDXc2SIYAWQaY3YYY6qB94FJjU8wxiw1xjRMl1gGJLoxHtXEip2FgK446ml+fsKYPvEsycijXncta1FNXT1PfpbOQx+kkZIUw8f3jKJnQkSz5/7h0gH4+wlPfZ7ezlFakzsTQVcgq9HjbOexltwBfNXSkyJyl4isEpFVeXnahD5VWYXlTP8ug8GJ0STGhno6HJ83pm8CB8trWLhNBzmbk1tSyZSZy5i9dBd3nNuDd39+Fh3Cg1o8v3N0CA9d0IdvN+fy7aacdozUmtyZCJpbxrLZnzsiMg5HIvhtS29mjJlpjEk1xqQmJCS0UYi+qaq2jnvfW0O9MbwyZZiuOOoFzuubQHxEMLfPXsWNbyxj8bY8XYzOaeWuQi59eQnp+0p4acpQ/nDZgFZNbrhtVA/6dIzgqS904PhE3JkIsoGkRo8TgX1NTxKRwcAbwCRjTIEb41FOz3y5mfXZxfztuiEkx+n+A94gNjyI7x85j99N7E9GziFumbWCS19awqdpe6mtq/d0eB5hjOHNH3cyZeYyIoID+OTeUVwx5LRWvz7Q34+nJw0iq7CCVxdud2Ok1ifu+tUhIgHANuB8YC+wErjBGJPe6Jxk4HvgFmPM0ta+d2pqqlm1alUbR+wbPlu3jwfmrOXO0T34/aVNx+6VN6iqrePTtH38Y9F2tueVkRgbSkpSjEstt8TYUH5zUT9Lt/ae/Cyd2Ut3ccHpnXj++iFEhZzcWlgPvr+WrzYe4OuHxtAjPryNo7QWEVltjEk95rg7m58icgnwIuAPzDLGPCMi0wCMMa+JyBvANcBu50tqmwuyKU0EJycz9xBXvLKEAV2imHPX2Vo74OXq6w3fb8ll1o87OVDc+hqDqtp69hZVMGtqKuP7d3JjhO6TmVvKhS8sZsqIZP530iD8/E4+oeWWVHLhC4vpHhfGR9NGEhTgu3/vPZII3EUTgevKq2u5csaP5B+q5ssHzqVLtA4Q21VNXT1j/7qQhMhg5t4z0pKtgoc/SOO/Gw+w5LfjiIsIPuX3++/G/Ux7Zw3TzuvFoxP7t0GE1tRSIvDd1OhDjDH8zycbycg9xPTJKZoEbC7Q34+7x/YiLauIJRbcE3lXfhmfpu3lprOT2yQJAFw8qAtTRiTzj8XbWaL1GsfQRGBzeaVVPP/NNj5es5cHxvdhdB+dceULrktNpHNUCC9/n+npUFz26sJMAv39uHNMzzZ938cvG0CvhAge/jCNgkNVbfreVqeJwIaKy2v4YOUebnpjOWc9+y0vf5/JRQM78cD5fTwdmmonwQH+/OK8nqzYWcjyHdaZjJdVWM7Ha/YyZUQyHSND2vS9Q4P8eWnyUIoravj1v9fr9NxGAjwdgGo7SzLyeeunXSzamkd1XT3d4sK4d1xvrhhyGn066eqivmbKiGRmLNjOy99nclbPuBO/wAu8tmg7IvCL89q2NdBgwGlRPDaxP09+vom3lu5i6qgebvkcq9FEYBOZuYeY+uYKOoQHcfM53bhiyGkMToy25EChahshgf7cNaYHz87bwpo9BxmW7N17ThworuSjVdlce2aSW8exbh3ZnR8y8nl23hZG9IhjwGlRrX5tXb2htLKGmLCWq5qtSLuGbOJ/v9xEaJA/Xz04mj9cNoAhLs47V/Z041ndiA0L5OXvMjwdygm9tmg7dcZwz9hebv0cEeEv1w4mJiyQ++esoaK6dVXHtXX1/OLt1Yx45jte+T6DGhsV+mkisIEFW3JZuKWK9dUAABOaSURBVDWPB8/v02azLJQ9hAcHcMe5PViwNY8N2cWeDqdFuaWVzFmxh6uGdiWpg/ur3eMignnh+hR25JfxwPtrqao9fjJomHn37eYczkiM5m/zt3H5y0u8+r+pKzQRWFx1bT1//HITPePDueWc7p4OR3mhW0Z2JyokgFcWeG+r4I0fdlJTV8+943q322eO6h3PU1cM5JtNOdwxexXl1bUtnvvitxm8v9KxC9p/7h7JzJvPpLCsmitf/ZE/fbXZ8msZaSKwuH/9tIsdeWX84bIBPl0xqVoWFRLI1FE9+Do9hy0HSjwdzjEKy6p5Z9luLh9yWrsvAXHLOd3523VDWLo9n5veWE5xec0x57y7fDfTv8vgZ6mJ/GpCXwAmDOzMN788j5+lJvKPRTuYOP0HlllodlZT+s1hYQWHqpj+XQbn9U1gXP+Ong5HebHbR3UnPMifV7ywrmDWkp1U1NRxXzu2Bhq79sxEXr1xGBv3ljD59WXklR6pMZiffoA/fLKRcf0SeOaqM44ad4sODeRPVw/mvZ+fRV294YbXl7FxrzW7ijQRWNhz32yjorqOP1x2uqdDUV4uJiyIW0Z258sN+73qy2p/cQWzl+5i4qDOHp3ifPGgLvxzaiq78su4/h8/sbeogtW7C7l/zlrOSIxhxo3DWlyba2TveD6/71yCA/z510+72jXutqKJwKLS9xUzZ8Uebj6nG707ao2AOrG7RvekY2Qw989ZS2nlsV0g7c0Yw2/+vZ66esNvLvL8+j+j+yTw9h0jyDtUxXV/X8odb63itJhQZt2aSljQ8WfaR4cFcuXQrnyatq/Z7iVvp4nAgowxPP35JmJCA3no/L6eDkdZRGx4EC9NHsrugjIem7vR45W17yzfww8Z+Tx26el095LloVO7d2DOnWdTVVtPgJ8f/7p9RKtn4t10djJVtfV8tDrrxCd7GU0EFvTfjQdYvrOQX03oR3TYya3RrnzTWT3j+OWFffl83T7mrPDcF9au/DKe/XIzo/vEc9NZyR6LozmDukYz/+ExfPXgaJemsg48LZozu8Xy7vI9ltt7WhPBcdTU1bN2z0FPh3GUypo6npm3mf6dI5k8POnEL1CqiXvG9mZ0n3ie+jydzfvbfxZRXb3hkY/WEeDvKOzyxsLHuIhgEiJdr8m5+exu7Mwv48ft1lrhVBPBcUz/NoOrXl3Kmz/u9HQoh72+eAfZByt4/LIBBOjGMuok+PkJz/8shajQQO59bw1lVS3Pn3eH13/YwardB3l60kDbLYk+8YzOdAgP4u2fdp/4ZC+i3yQtqKyp493luwny9+PpLzbx3437PR0S+4oqmLEwk0vO6MzI3vGeDkdZWEJkMNMnp7Arv4z/+aT9xgu2HCjh+fnbuHhgZ65M6doun9meggP8uX54Et9uzmFfUYWnw2k1TQQt+CxtHwfLa3jt5mGkJMXw4PtprN7t2W6iZ+dtxhh47BKdLqpO3che8Txwfh/mrt3LR6uyDx83xpB/qIq0rCK+XL+f7IPlbfJ51bX1/PKDdUSFBvDMVYO8skuoLdwwIhkDzFmxx9OhtJquPtoMYwxvLt1Fv06RjOvXkSGJMVzz96X8/K2VfHzPqBarH5dm5vPitxmc3SuOB8b3btOum+U7Cvhi/X4ePL8PibHuX4tF+Yb7x/dhxc5CHv9sI1+nHyDrYDnZBysob7QQ26CuUXx+37mn/MX98vcZbNpfwsybz7T1mlhJHcIY368jc1Zkcf/4Ppao+Pf+CD1gxc5CNu8vYeqo7ogIcRHBzL5tBCLC1DdXHLO7UVZhOdPeXs0NbyxnW24pL32XwfUzl7XZL6naunqe+CydrjGhTDvPvSszKt/i7ye8ODmFAV2i2FtUQbe4cCYPT+bxywbw+i2p/ObifmzcW8L3W3JP6XM27i3m1YXbuWZYIhMGdm6j6L3XTed0I/9QFV+nH/B0KK2iLYJmzF66i5iwwKP6MLvHh/PGralMmbmMO95axZw7z6beGF5dmMnrP+zEX4RHJvTl56N78nX6AX4/dyOXTP+B/7tmMBPP6HJK8cxZmcWWA6W8euMwQoP8T/XylDpKx8gQPr5nVLPPje2XwJwVe5j+XQbj+3c8qVZBfb3h8U83EhsWyOOXDzjVcC3hvD4JJHUI5W3nGkreTlsETewtquDr9ANMHp58zJfusORYpk8eyrrsIm6dtYLxzy1kxoLtXDKoMwseGct94/sQEujPpJSufPnAufSID+fud9fwu483tHrN86YOllXz3PytnNMzjomD7P9LSnmXQH8/7h3bm/XZxSzclndS7zF37V7W7CniNxf3JzrUN+pe/PyEm87qxoqdhWw9UOrpcE5IE0ETDdO+bj6nW7PPXzyoM49fNoAVuwrpFBXCf+4+hxcnD6Vz9NH7q3aLC+ejaSOZdl4v5qzYw6QZS05qzvbz32yjpKKGJ64YYNvBNeXdrh6WSNeYUKZ/m+Hy7KKSyhr+9NUWUpJiuHZYopsi9E7XpSYRFODHO8u8fyqpJoJGKqrreH/lHi4a2JmuMS3Pb75tVA8WPjKWT+4ZxZndOrR4XlCAH49O7M/bd4ygsKyGy15ewtOfb6Kkleu8bN5fwrvLd3Pz2d3o37n12+kp1ZaCAvy4d1xv0rKKWJzhWqHU9G8zKCir4ulJA/Hz860fMh3Cg7hscBc+XpPNoXau1XCVJoJGPk3bS1F5DVNHdj/hud3jw1v9F3t0nwS+eXgMk4cn8ebSnYz/20I+XJV13DJ0YwxPfpZOdGggD1+o6wkpz7r2zIZWwbZWtwq25ZQye+kuJg9PZnBijJsj9E43n92Nsuo6pn+7zdOhHJcmAidjDLOX7uL0LlGM6NHyr/yTFRsexDNXncHn951LcocwfvPv9Vz996Wszy4CHDOD0vcV886y3Tzy0TrOf34Ry3cW8shF/Wy3UbaynqAAP+4e24s1e4pYknniVoExhic+TSciOIBfX9SvHSL0TkOTY7nxrGRe/2En73txXYHOGnJatqOQLQdK+cs17l37ZFDXaP49bSRz1+7lT19tYdKMHxncNZptOYeocG53FxcexNDkGG4+uxuTh3vXglzKd12XmsiMBZlM/zaDc3vHH/ffybwNB/hpRwF/vHIQHcJ9+4fMk1cMZE9hOf/zyUaSO4R55aoAmgicZi/dSWxYIFekuH+ql5+fcM2ZiVw4sBMvf5fB2j1FXD88iaHJMQxLjiUxNlQHhpXXCQ7w5+6xvXj803R+2l7Q4hdaeXUt//vlJgZ0ieKGEfpDJtDfjxk3DuOaV5cy7Z3VfHzPKHp3jPB0WEfRRICjIOybTTlMO68XIYHtN08/KiSQ31/qG/OqlT38LDWJGQsyefG7jBYTwYwFmewvruTlKUPx97EB4pZEhQQya+pwrnr1R26fvZJP7h3VbEvpYFk1763Yw+6CMrrGhJEYG0rX2FASY0PpHBXitoUmfTYR1NTVsy6riB8y8vnvxgOICDed3fyUUaWUQ0igP3ef14snP9/ET9sLGNGjA3mlVRwoqeRAcQV7iyp5ffFOrh7aldTubT/WZmVJHcKYeUsqk2cu4xdvr+Kdn59FcIDjh2dWYTn/XLKTD1ZmUVFTR3xEMPlNVjDw9xO6RIcw78HRRIW0bT2GTyWCHXmH+CEjnx8y8lm2o4BDVbWIwOCu0fzp6jM47ThTRpVSDpNHJDNj4XamvrmCmrp6mk5+6xoTyqMTPb/1pDcalhzLc9cN4f45a3n0PxuYOrI7M3/YwVcb9uPvJ0xK6cqdo3vSr3MkVbV17CuqJNu5/lP2wXL2F1cSGdz2X9vi6e3qTkZqaqpZtWqVy68b85cF7CksJ6lDKOf2TmB0n3hG9orTWTlKuWjBlly+2rifTlEhdIoKoXNUCJ2jHbcOYUE+VzPgqle+z+Bv8x1TSiNDArjxrG5MHdn9mMLUtiYiq40xqU2P+1SL4K/XDqZLdCjJcbp6p1KnYlz/jozr39HTYVjWveN6IyIEB/hx/fAkItu4q8dVPpUIzuoZ5+kQlFIKEeHecb09HcZhWlCmlFI+zq2JQEQuFpGtIpIpIo8287yIyEvO59eLyDB3xqOUUupYbksEIuIPzAAmAgOAKSLSdNL8RKCP83YX8Hd3xaOUUqp57mwRjAAyjTE7jDHVwPvApCbnTAL+ZRyWATEicmq7uCillHKJOxNBVyCr0eNs5zFXzwFARO4SkVUisiov7+Q2yFBKKXUsdyaC5iYSNy1aaM05joPGzDTGpBpjUhMSEk45OKWUUg7uTATZQFKjx4nAvpM4RymllBu5MxGsBPqISA8RCQImA581Oecz4Bbn7KGzgWJjzH43xqSUUqoJtxWUGWNqReQ+4GvAH5hljEkXkWnO518D5gGXAJlAOXBba9579erV+SJyoo1A4wHX9tWzB71u36LX7VtO9bqbXVnTkmsNtYaIrGpuTQ270+v2LXrdvsVd162VxUop5eM0ESillI+zcyKY6ekAPESv27fodfsWt1y3bccIlFJKtY6dWwRKKaVaQROBUkr5ONslghMtfW0nIjJLRHJFZGOjYx1E5BsRyXD+GevJGNuaiCSJyAIR2Swi6SLyoPO43a87RERWiMg653U/5Txu6+tuICL+IrJWRL5wPvaV694lIhtEJE1EVjmPtfm12yoRtHLpazuZDVzc5NijwHfGmD7Ad87HdlIL/MoYczpwNnCv8/+x3a+7ChhvjBkCpAAXO6vx7X7dDR4ENjd67CvXDTDOGJPSqH6gza/dVomA1i19bRvGmMVAYZPDk4C3nPffAq5s16DczBiz3xizxnm/FMeXQ1fsf93GGHPI+TDQeTPY/LoBRCQRuBR4o9Fh21/3cbT5tdstEbR6WWsb69SwXpPzT9vuMC4i3YGhwHJ84Lqd3SNpQC7wjTHGJ64beBH4DVDf6JgvXDc4kv18EVktInc5j7X5tdtt8/pWL2utrE1EIoD/AA8ZY0pEmvtfby/GmDogRURigLkiMsjTMbmbiFwG5BpjVovIWE/H4wGjjDH7RKQj8I2IbHHHh9itRaDLWkNOwy5vzj9zPRxPmxORQBxJ4F1jzMfOw7a/7gbGmCJgIY7xIbtf9yjgChHZhaOrd7yIvIP9rxsAY8w+55+5wFwc3d9tfu12SwStWfra7j4DbnXevxX41IOxtDlx/PT/J7DZGPN8o6fsft0JzpYAIhIKXABswebXbYz5nTEm0RjTHce/5++NMTdh8+sGEJFwEYlsuA9MADbihmu3XWWxiFyCo0+xYenrZzwcktuIyBxgLI6laXOAJ4BPgA+BZGAPcJ0xpumAsmWJyLnAD8AGjvQZP4ZjnMDO1z0Yx8CgP44fcB8aY54WkThsfN2NObuGHjHGXOYL1y0iPXG0AsDRjf+eMeYZd1y77RKBUkop19ita0gppZSLNBEopZSP00SglFI+ThOBUkr5OE0ESinl4zQRKMsQkQecq46+KyJXuLK6rIh0F5EbWnhubMOqlt5ARFKc06CVahd2W2JC2ds9wERjzE7n42OKBUUkwBhT28xruwM3AO+5L7w2kwKkAvM8HYjyDZoIlCWIyGtAT+AzEZkFHARSjTH3ichsHKuwDgXWiMhnwHTnSw0wBvgzcLpz0ba3jDEvtPA5w3HsC3uNMWZHo+MDgTeBIBwt6WuAW4B8Y8x05znP4CjsWw885byfAnyMowDuQSAUuNIYs90ZdyUwEOgE/BKYDzwNhDqL5/4EfAPMcl5/OXCXMWa9iDwJ9AC6AH2drz8bxzLse4HLjTE1IvJn4AocS3jPN8Y80vr/8sonGGP0pjdL3IBdQLzz/lTgFef92cAXgL/z8ec4FusCiMDxg2cs8EUL7zvW+fqRwGoguZlzXgZudN4PwvGF3h1Y4zzmB2wH4pzvV4TjCzoYx5fyU87zHgRebBT3f52v7YNjrayQxtfW6LOfcN4fD6Q57z8JLMGxJPUQHEliovO5uTiWJ+4AbOVI8WiMp/8/6s37bjpGoOziI+NYnRPgR+B5EXkAxxdfc11FTZ2OoyVwuTFmTzPP/wQ8JiK/BboZYyqMMbuAAhEZimMdmLXGmALn+SuNY++EKhwJYr7z+AYcCaTBh8aYemNMBrAD6N/MZ58LvA1gjPkeiBORaOdzXxljapzv648jsTT+nBIcrY43RORqHMlCqaNoIlB2UdZwxxjzZ+DnOH61LxOR5r5cm9qP4wtzaHNPGmPew9G9UgF8LSLjnU+9geMX/G04um8aVDW6X9/ocT1Hd8k2XeOluTVfjre8epUzvnqgxhjTcLweaBgvGYFjtdYrOZIolDpME4GyHRHpZYzZYIz5P2AVjl/ZpUDkcV5WhGMXrGebW/feuQDYDmPMSzgGqQc7n5qLYzno4cDXJxHudSLiJyK9cIwBbG0m1sXAjc44xuIYlyhpzZs7922INsbMAx7CMWah1FF0sFjZ0UMiMg6oAzYBX+H4hVwrIuuA2aaZwWJjTI6IXA58JSK3G8cOYA2uB24SkRrgAI4BXYwx1SKyAChq1DXliq3AIhyDxdOMMZXO93vUObD9JxxjAW+KyHocXTu3tvRmzYgEPhWREBwti4dPIkZlc7r6qFKnQET8gDU4lgLOcPG1s3EMYP/bHbEp1VraNaTUSRKRAUAm8J2rSUApb6ItAqWU8nHaIlBKKR+niUAppXycJgKllPJxmgiUUsrHaSJQSikf9//lTtj9sdNXZAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "prevalence = 0.072\n",
    "\n",
    "num_symptoms, probs_positive = [], []\n",
    "all_symptoms = nigam_df[\"Clinical observation\"].values.tolist()\n",
    "for i in range(len(all_symptoms)):\n",
    "    symptoms = all_symptoms[0 : i+1]\n",
    "    prob_positive = compute_prob_positive(prevalence, symptoms, nigam_df)\n",
    "    num_symptoms.append(i + 1)\n",
    "    probs_positive.append(prob_positive)\n",
    "    \n",
    "plt.plot(num_symptoms, probs_positive)\n",
    "plt.xlabel(\"first k symptoms\")\n",
    "plt.ylabel(\"P(+|symptoms)\")\n",
    "_ = plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### How does the predictions change with change in prevalence?\n",
    "\n",
    "As the prevalence increases, the shape of the chart does not really change, except the P(+|.) value goes higher. This is expected, since all that is changing is the P(D) multiplier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOydd3hUVfrHP2dKJm0SEtITIAk19BJ6k64UXRUL2FGxwE/ddXFdXdey66prWxUFdVXEVdlVsSNFekc6EggklJBCei+TKef3xyQhfVomBe/neeaB3HvOuWeSmfvec973/b5CSomCgoKCwm8XVVtPQEFBQUGhbVEMgYKCgsJvHMUQKCgoKPzGUQyBgoKCwm8cxRAoKCgo/MbRtPUEnCEoKEhGR0e39TQUFBQUOhQHDhzIkVIG1z/eIQ1BdHQ0+/fvb+tpKCgoKHQohBDnGzuubA0pKCgo/MZRDIGCgoLCbxzFECgoKCj8xumQPgIFBYWOg9FoJDU1lYqKiraeym8GT09PoqKi0Gq1drV3qyEQQnwIzAaypJT9GzkvgDeAmUAZcKeU8qA756SgoNC6pKamotfriY6OxvqVV3AnUkpyc3NJTU0lJibGrj7u3hpaAVzZzPmrgJ5Vr4XAMjfPR0FBoZWpqKigc+fOihFoJYQQdO7c2aEVmFsNgZRyG5DXTJNrgJXSyh6gkxAi3J1zUlBQaH0UI9C6OPr7bmsfQSRwodbPqVXHMuo3FEIsxLpqoGvXrq0yOYX2h7RYMBcUYMrOxpSVbf03Jweh0aCNCEcbHo4mLBxNcBBCpcRCKCjYQ1sbgsbMVqMFEqSU7wHvAcTHxytFFC5jjJlZGFPOU5mahjE11fpKS6MyPQ1TVjaYTLYH0WrRhoSgCQpC5eON8PJG5eVlfXl7o9Lr0YaHo42MRBsViTY0FFHLsSalxJyXhzHjIsaMdEzZ2fhOmIBHVJQb37lCWzB37lz++c9/EhsbW+f4ihUr2L9/P0uXLmXp0qX4+Phw11132T3uCy+8wAcffIBarebNN99kxowZDdrk5eVx0003ce7cOaKjo/nf//5HQEAA+/btY+HChYD1s/jMM89w7bXXAnDgwAHuvPNOysvLmTlzJm+88YbLK662NgSpQJdaP0cB6W00F4U2xHDmLMXr1lK0dh2GxMRLJ4RAExqKNioSn+HDUYeGUqRXcV5bxEmRyUHzOY5aUvDGg7Hq3gwX0fQ2BdG5wIL5Yiam3BwspWVYcnKxlJdjKS9HlpVhKS+H2kWZVCo0YaFog0MwFxRgvHgRaTDUmWNOYCBdP/wAzz59WueXouB2jh8/jtlsbmAE6rNgwQLGjh1rtyFISEhg1apVHD9+nPT0dKZOncqpU6dQq9V12r344otMmTKFxx9/nBdffJEXX3yRl156if79+7N//340Gg0ZGRkMGjSIOXPmoNFoeOCBB3jvvfcYNWoUM2fOZO3atVx11VVO/w6g7Q3Bd8BiIcQqYCRQKKVssC2kcHnS2M3fa+hQQv70J3S9euIRFYU2PJw0Qxbb07azI20Hh7N2UFRZBIDeQ8/A4IEsCJpFibGEPRl7+CH/e9CCX6QfI+NHMij4CkK8QwjyCiLEO4Rgr2C8td7IykqMmZkY09LqvjKz0EWE4ztlinXFEB6GJjwcLBZS/+8hzt9+B13fexevwYPb8len4CDnzp3jyiuvZOTIkRw6dIhevXqxcuVKPv30U6655pqadh999BEvvPAC4eHh9OrVC51OB4C3tzfR0dHs27ePESNG2Lzet99+y80334xOpyMmJoYePXqwb98+Ro8e3aDdli1bALjjjju44ooreOmll/D29q5pU1FRUfPEn5GRQVFRUc04t99+O9988037NgRCiM+BK4AgIUQq8DSgBZBSLgfWYA0dTcIaPmr/ukuhQ5O99G1yli4FrDf/0Cf+jH76dLRhYRjNRg5lHWJb6hdsP7SdM4VnAOiq78q0btMYFDyIQcGDiPaPRiXq+gFyynPYm7GXPRl72J2+mw3nNzS4to/Wh3CfcEaFj2Jcl3EMi59FJ41ng3ZSSs4WnWVvxl4O5Rzi+tcfpfNjb3J+wd10eecdfEaNdMNv5vLm2e+Pk5Be1KJj9o3w4+k5/Wy2S0xM5IMPPmDs2LEsWLCAd955h507dzJv3jzAepN9+umnOXDgAP7+/kyaNIkhQ4bU9I+Pj2f79u2MGDGCl19+mU8//bTBNSZMmMCbb75JWloao0aNqjkeFRVFWlpag/aZmZmEh1vjY8LDw8nKyqo5t3fvXhYsWMD58+f55JNP0Gg0pKWlEVVre7KpcR3FrYZASjnPxnkJLHLnHBTaH4Xff0/O0qX4zZ5NyJI/og0NrTn3bdK3vPTLSxRXFqNVaYkPjeeGXjcwPmo83fy62Rw7yCuIWbGzmBU7CyklRZVFZJdlk11e9ar6/9nCs3xx6gv+c+I/eKo9iQ+LZ1zkOIaEDCGpIMlqTNL3kFVu/WKqhZpEv0S++M9KLtx9Dxfuu4+oN9/Ad+JEt/2eFFqWLl26MHbsWABuvfVW3nzzTTIyMggOtopx7t27lyuuuKLm55tuuolTp07V9A8JCeHkyZMALFmyhCVLljR5rcZqwTu6jz9y5EiOHz/OiRMnuOOOO7jqqqtaZNzGaOutIYXfGGUHD5HxxJN4Dx9OxD+eR3h4WI8by3h+7/N8l/wdw0KHcXvf2xkVPgpvrbeNEZtGCIG/zh9/nT89Ano0OF9hqmB/5n52pO1gR9oOXtz3Ys25TrpOjAwfycjwkYwKG8Wvub/y2LbH2FpxjEmffGI1Bov/j8iX/4nflc2lyijUxp4nd3dR/4YphMDLy6tOvH1zN9WKigq8vLwAbK4IoqKiuHDhUkBkamoqERERDdqHhoaSkZFBeHg4GRkZhISENGgTFxeHj48Pv/76K1FRUaSmptoc12GklB3uNWzYMKnQ8TBcSJWJo8fI09OnS2NeXs3xk7kn5ezVs+WAFQPk24feliazqU3ml1KYIr9P/l6eyD0hzRZznXMms0nOWj1L3vDdDdJisUhTUZE8O2++TIjrKwu+/6FN5ttRSEhIaOspyLNnz0pA7tq1S0op5T333CNfeeUVedNNN8kNGzZIKaVMT0+XXbt2lTk5ObKyslKOGzdOLlq0qGaMxYsXy88//9yu6/36669y4MCBsqKiQp45c0bGxMRIk6nh5/qPf/yjfOGFF6SUUr7wwgtyyZIlUkopz5w5I41Go5RSynPnzsnw8HCZnZ0tpZQyPj5e7t69W1osFnnllVfKH3/8sdE5NPZ7B/bLRu6pbX5Td+alGIKOh6m4WCbPniNPxg+XFcnJUkopLRaL/O/J/8qhK4fKSf+dJPem73Vu8KSNUqYeaMHZNs7qU6tl/xX95dYLW6WUUppLS+WZG26UpyZNkhaLxe3X76i0F0MQFxcn77vvPjlgwAB53XXXydLSUrly5Ur55JNP1rT78MMPZc+ePeWECRPkQw89VMcQDBkypOZmbA9///vfZWxsrOzVq5dcs2ZNzfG7775b/vLLL1JKKXNycuTkyZNljx495OTJk2Vubq6UUsqVK1fKvn37ykGDBskhQ4bIr7/+uqb/L7/8Ivv16ydjY2PlokWLmvzsKYZAoV1hMZnk+YULZULffrJk504ppZSllaXy0S2Pyv4r+sv71t8nc8pyHB/YZJRy7RNSPu0n5XPBUp5cY7uPC1SaK+W0L6bJW368pebLl//VapnQu48sO3rMrdfuyLQXQ9CvX78Gx8vKyuTIkSMbfVqvzcGDB+Wtt97qrum5BUcMgZJ6qdAipD78CGeuvob0J54k77PPKD96FEtlJQBZ//wnpVu3EfbUX/AZMwaA5UeXs+H8Bh4Z+gjvTH2Hzl6dHbtgSRasvAZ2L4X4BRDaD1bdAkf/19JvrQatSsuC/gs4kn2EXy7+AoB+8iTQaChev95t11VwH15eXjz77LM2I29ycnL429/+1kqzan2E1Uh0LOLj46VSqrL9YExPJ2nyFDy6d8ecl4c5P996QqtFF90Nw+kkAm6/jbAnngCg0lzJ1C+mMix0GK9Pet3xC174Bf53O5TnwZw3YNDNYCiGz+fBuR0w82UYcW8LvsNLGMwGrvzqSrp36s6/p/8bgJQFd2NMSyN27U+Kpk4jnDhxgri4uLaexm+Oxn7vQogDUsr4+m2VFYGCyxRv3gxA1Ftv0XPXTnps2kjkm2/Q+a670AQHEzB/PqGPPVbTfmPKRvIN+dzQ6wbHLiQl/PJv+OgqUGvh7g1WIwCg08MtX0CvK2HNH2H7qy319uqgU+u4s9+d7M3Yy5HsIwDop0+j8vx5DKdPu+WaCgruRjEECi5TsmkzHtHR6GJjEEKgjYjAb/p0Qv7we7p++CFhf30KobkUqfzlqS+J9I1kVMSoZkath9kI3y2GHx+F2Ctg4RYIH1i3jdYLbvoEBtwIG5+DDX+tKyPRQtzQ6wb8df68f/R9APRTpoAQFK9vmLymoNARUAyBgkuYS0oo3bcP38mT7Wp/rvAc+y7uY26vuQ2ygpu+iBG+XACH/gMTlsD8/4F3YONt1Vq49l2Ivxt2vgFf3gXnd4HFYuc7so231pvb4m5ja+pWTuadRBMcjNfQoRRvUAyBQsdEMQQKLlG6YwcYjVanqR18eepLNELD73r8zr4LmCrhizvhxHcw4wWY/BewJS+tUsGsV2Hi45D4k3Ur6fW+8NPjcGFfixiFeXHz8NX61qwK/KZPw5CYSOX58y6PraDQ2iiGQMElijdtQt2pk10ibJXmSr5N/pZJXScR5BVke3CTweoUPvkDXPUyjH7Q/okJAZP+DEuS4foPIGIo7P8APpgG/xoAPz9rHd9ezKY6BsTPw495feax4fwGzhSeQT91KoCyKuiAzJ07lzNnzjQ4vmLFChYvXgzA0qVL+eijjxwa94UXXqBHjx707t2bdevWNdrmiy++oF+/fqhUKuoHwDTV/8CBAwwYMIAePXrw0EMP0RIBP4ohUHAaaTJRsnUbvhMn1vEBNMXP53+mwFDA3F5zbQ9urID/3gqnfrI+3Y9c6Nwkdb4wYC7M+wyWJFm3jUL7wY7X4LMbwVBie4ykjfDPWNi7vM7hW/veiqfGk4+Pf4w2MhLP/v0pUvwEHQpHZKjffPNNu8etLUO9du1aHnzwQcxmc4N2/fv3Z/Xq1UyYMMHu/tUy1KdPn+b06dOsXbvW7nk1hWIIFJym7OBBLIWFdvsHvjz9JVG+UYwKt+EkNpbDf2+B0+th9r9g+D0tMFvA098aZXTL/+B3y+HsNvjkd1Ce33Sfgyvh0xvAUAjnttc5FegZyLjIcey/aH2S00+bRsXRoxgzFCX19sa5c+fo06cPd9xxBwMHDmTu3LmUlZU1KkPdq1cvJk6cyM6dO2uO15ahtoemZKjrExcXR+/eve3uX1uGWghRI0PtKoronILTlGzajNBq8alSdGyOs4Vn+eXiLzw89OHmncQmgzUf4MwWuHopDL2t5SZcm8HzrKuFLxfAitlw62rQX1JBRUrY9HfY/gp0nwxCBdknGwzTs1NPfj7/M+WmcvTTp5H9+usUb/iZwNvdNO+Ozk+Pw8VjLTtm2AC46kWbzdqjDHVTNNVfq9W6RYZaWREoOIWUkuJNm/AeNQq1r4/N9nY7iQ9/Cmc2w9Vvus8IVBM3xxqBlHcGProSClKsx00GWH2v1QgMvd3aJjIe8s5aVyu16BHQA4nkTOEZdDEx6Hr2VLKM2yn1Zah37NjRpAy1h4cHN910U53+ISEhpKdbCyguWbKEw4cPN3hVbx81tm/vSLJhU/1dHbcplBWBglNUnjmDMSWFznfdabOtwWzgu+TvbDuJpYS970LYQBjSSk/U3SfB7d/Cp3PhwyvhhhXw8zNwfidM+SuM+4PV8RzSB5CQcwrCB9V079HJKm+dlJ9Ev8790E+bRs6yZZhyctAE2eEQ/61hx5O7u2iPMtRN0VR/d8lQKysCBaco3rQJAN9J1rDRxLxETuY13DqBS05im5nEZ7ZYt19GPWC9+bYWXUbAnT+CudIaVZT6izXSaPyjl+YRXJWqn1X3PXbRd8FD5UFSQRIA+hnTQUqKN25qvfkr2EVKSgq7d+8G4PPPP2fcuHHExcWRlGT9240cOZItW7aQm5uL0Wjkiy++qNP/1KlT9O/fH7C9Irj66qtZtWoVBoOBs2fPcvr0abtKXFbTVP/w8HD0ej179uxBSsnKlSvr+DicRTEECk5Rsmkznn37og0LA+DBjQ9yw/c3MP/H+XyT9A0VpktPWV+c+oIu+i6MDLdR2nHvcvAJhv7Xu3PqjRM2ABasgz6z4bZvrJFGtQmMBZWmgZ9Ao9IQ2ymW0wVWeQldr15ou3ZVwkjbIXFxcXz88ccMHDiQvLw8HnjgAWbNmlVTMzg8PJxnnnmG0aNHM3XqVIYOHVqn/86dO5laFSZsi379+nHjjTfSt29frrzySt5+++2awvX33HNPTajo119/TVRUFLt372bWrFnMmDHDZv9ly5Zxzz330KNHD7p37+5yvWJAkaFWcBxjbq5M6BMns95aKqWUMq88zyonveE+OefrObL/iv5yzGdj5Ev7XpJbUrbI/iv6y38f/Xfzg+YkSfm0v5Qb/94K78BJlo6Q8rObGxx+fNvjcsr/ptT8nPnyyzKhX39pKihozdm1WxQZ6rZBkaFWcCslW7aClDXZxNXbIrfF3ca313zLhzM+ZFT4KD4/8TmLNy1Go7LDSbzvfesT9/C73T195wnuA1knGhzu0akHmWWZFFVai7Lrp00Dk4mSqidNhfaLIkNtRXEWKzTAlJtL9ltvETBvPp69ezU4X7J5E5qwMHRVErfJBckAdO/UHSEEw8OGMzxsONll2XyT9A3+Ov/m6w1UFFl1hPpdC/owt7ynFiEkDhK+hcoy8LhUS7lnQE/A+nsYEjIEzwED0ISFUbRuPf4tsH+r4DrR0dH8+uuvjZ6r3o5pjmnTprX0lNoVyopAoQEFq1dTsOq/nJ07l5zly5EmU805i8FAyY6d6CdPqomwSCpIwlfrS6h3aJ1xgr2DuXfgvdzY+8bmL3j4M6gshlH3t/h7aVGCa0UO1aI6cuh0vtVPIFQq/GbOpGTTJlIf+T3GFojzVlBwJ4ohUGhAycZN6Hr2xG/aVLL/9Qbnbp5Xo7VftmcPsrwc30mXsomTC5JrVgMOY7HAvnchajhEDnNqvudzS8ksqrDd0FWC+1j/zU6sczjcJxxvjXfNFhlA8MMPEfTQ/1GyZQvJM2eRvfRtLBWtMEcFBSdQDIFCHUw5OZQfOYL+yhlEvvYakf96HWNaGmevu56cd9+jaMMGVN7eeI+8FAqXXJBc81TsMEkbrAldI51bDVQYzdywfDcPfnrQues7QufuoNJCdl0/gRCCHgE96hgClU5H8IMP0n3Nj/hOnkTO0qWcmTmLorXrWkQkTEGhJVEMgUIdSrZWO4KtT/x+V15J7A/f4ztpEtmvv07hl1/hM24cKg8PAHLLc8k35DtvCPYsA30E9HVuL331wTSyig0cOJ9P4sVi5+ZgL2otdO7RIJcArFITp/NPN7jJayMiiHr9dbp+/DEqX1/SHnmEC3ffg6W8vMEYCgpthWIIFOpQvGkzmvBwdH361BzTdO5M5Bv/IuLVV9B260qnGy4lhtV2FDtM1kmrnMTwu603WQcxWyTvbUumd6geD7WKz/elOD4HRwnp02BFAFY/QYGhgNyK3Ea7+YwcQczqrwhZ8kdKd+1S8gzaEW0pQ/3UU08xcOBABg8ezPTp02skLJrrr8hQK7gVS0UFpTt3op80qdF0fP9Zs+ixbh2+48fVHK9OpHJqRbB3Oah1MOxOp+a79teLnMst45GpPZneL5SvD6VRYWwo9duiBMdB/nlr5FAtegRUSU3U2h6qj9BoCLzrLrQRERR+/4Nbp6lgH20tQ71kyRKOHj3K4cOHmT17Ns8995zN/ooMtYJbKd29G1lRge8U+2Slwboi8PPws6/QTG3K8+HIKhh4A/g4rskjpWTZ1iRig3yY3i+M+SO6Ulhu5Kdf3SwBHdyb5iKHkvKbNgRQFVE0axalu3Zhym189aDQ8rRXGWo/P7+a/5eWltY8gCky1AptRsmmTah8fPAZPtzuPtWOYocjhn79CkzlTjuJdyTl8GtaES9dPwC1SjAqtjPdOnvz+b4LXDskyvYAzhJSpTmUfRIiLlVl6+zZmQBdQLMrgmr85swm9/33KfppLYG33uKumbZLXtr3UpOaVM7SJ7APfxrxJ5vt2qsM9ZNPPsnKlSvx9/dn8+bNgCJDrdBGSIuF4s1b8JkwHlHlCLbZR0qSCpKc8w9kHAWvQKvGjxMs35pMqJ+O3w2JBEClEtw0vAv7zuaRnG1H1TFnCYy1Rg5lNR45VL1V1hyevXqh69WLoh+U7aHWpL3KUD///PNcuHCBW265haVLlzbb35FxHUFZESgAUHHsGOacnJpoIXvIKc+hqLLIOUOQc+pSXL6DHLlQwM6kXJ6Y2QedRl1zfO6wKF5bf4pV+1J4clZfp8a2iVoLQT0bLVLTo1MPvk36FimlzS+n35zZZL/6GpUXLuDRpYt75toOsefJ3V20dxnq+fPnM2vWLJ599tnLS4ZaCHGlECJRCJEkhHi8kfP+QojvhRBHhBDHhRB3uXM+Ck1TvGkzqNX4jh9vd5/qbRCHHcVSWm+kwQ3lK+xh+dZk/Dw1zBvRtc7xEL0nU+NC+epgGgaTG53Gwb2b1BwqM5WRUWrbT+E/cyYART/+2OLTU2ic9ihDffr0pRXkd999R5+qaL3LRoZaCKEG3gauAvoC84QQ9R/TFgEJUspBwBXAq0II+/YlFFqUkk2b8I6PR92pk919qkNHHTYEpTlWZ7ETK4Iz2SWsPX6R20Z3Q+/ZMOT05hFdyCutZENCpsNj201wnLWaWWVpncPVmkP2+Am0kZF4xQ+j8PsflASzVqI9ylA//vjj9O/fn4EDB7J+/XreeOMNm/07lAw1MBpYV+vnPwN/rtfmz8A7gABigCRAZWtsRYa6ZTGkpMiE3n1k7ooVDvV7eufTcvzn4x2/4JltUj7tJ2XSRoe7/unLI7LXk2tkVlFFo+dNZosc88JGOf/93Y7Py16Of2Odf9rBOocLDYX2SW5Xkff55zKhdx9Z3g5kmt2JIkPdNrQXGepI4EKtn1OrjtVmKRAHpAPHgIellJbGBhNCLBRC7BdC7M/OznbHfH+zlFRFKlRXG7MXpx3F1fvrQb0d6pZZVMHqg2ncGN+FYL2u0TbqKqfxzqRczueWNtrGZZqoVubn4Ueod6hdKwIA/YwZoNEoOQUthJQSi9HoUB9FhtqKOw1BY16X+mvgGcBhIAIYDCwVQvg16AVIKd+TUsZLKeOrvfwKLUPxxk3oevbAo2tX242rkFLWiM05TM4p8NCDn2NOrg93nMVksXDv+OaTf26Ij0IlYNUvF5pt5zTVkUONZRjX0xxqDk1AAL7jxlH0449IS6PPPwp2Ik0mKs+fx5CYiLm04QOALRnqrjY++9OmTSM6OrolptoucachSAVqh0NEYX3yr81dwOqqVUsScBZwLpREwSnMhYWU7d9fR03UHjLLMikxljiXUVztKHYg7K24wsine1OYPTCCrp29m20b7u/F5D4hfLE/FaPZDTdYtcYaOdSE5tCZgjOYLKZGOjbEb85sTJmZlP2yv6Vn+ZvBUl6OITkZS2kpQq3GlJGh+F0cxJ2G4BegpxAipsoBfDPwXb02KcAUACFEKNAbaCj6oeA2SrZtB7MZvQPZxOCixlB2osOO4tUH0ygxmLh7XIxd7W8e3pWcEgMbT2Q5Pj97CO7TZAhppaWSC8X2rUb0kyYhvL2VnAInMeXnY6jSCdLFxKCNiMBSUYE5N6+NZ9axcJshkFKagMXAOuAE8D8p5XEhxP1CiOp00r8BY4QQx4CNwJ+klDnumpNCQ0o2b0IdFITnAMcSu5wOHS3Ph5LMKqkG+5BSsnL3OQZ16cSgLvZFNV3RO5hQPx0f7TyLyR2rgpA4KDjfIHLIHs2h2qi8vdFPnULRunVYKitbfJqXK9JiwZiejjEtDZW3N7ru3VF5e6Py80Pl64spKxPpoL/gt4xb8wiklGuklL2klN2llM9XHVsupVxe9f90KeV0KeUAKWV/KeV/3DkfhbrIykpKtm1HP+kKhMqxj0JyQTKBnoEEeAY4dtHsKo0eBxzFu5JzSc4u5fZR3ezuo1GreGBid/aezeOuFb9QWN7CN4UmitTE+sciEDY1h2rjP2cOlqIiSrdvb8kZXrZYjEYqz57DlJeHJigIj+hohMaaGyuEQBsejpQSY6YbQ4gvMxSJid8wZfv3Yykpcdg/ANYn3p6dejp+0ZyqG6cDK4KPd50j0MeDWQPDHbrUnWNjeOn6AexOzuW6d3ZyLqcFo4hqDEHd7SEvjRdd9F3skpqoxmf0aNSBgUr0kB1YysupTD6DxVCBR5cuaMPCGmQDq3Q6NEFBmAsKGnUc16ctZaifeeYZIiMjGTx4MIMHD2bNmjU2+ysy1AothpSS3I9WWEXmxox2uK/TEUPZiaDxgk72RSilFZTz84lMbhreBU+t2naHetw0vCv/uWckuaWV/O6dnexKbqGdx8BYUHs06Sewd2sIrPLUflddRcnmzZhL3KiT1MExFxVhOHsWBOhiY1H7+zfZVhMUhNBqMaU37zhuaxlqgN///vc1mckzqzLOFRlqhVah6PvvKd2+neBHHkHl6elQ34zSDMpMZc7nEAT1AJV9N/XP9p4H4JaR9oe21mdUbGe+XTSWIF8dt3+wr2UK2Kg10LnxyKEeAT1IKUqh0mz/nr//nNlIg4EL999P2cFDrs/vMkJKiSknh8qUFFQ6HR6xsTY/s0KtRhsejsVQgTk3t93KUDvaX5GhVmgxTHl5ZP7jBbwGDyZg/jyH+zvtKAarj6DrSLuaGkxmVu27wJS4UKICmg8ZtUW3zj6sfnAMiz87xJ9XHyMpq4S/zIpzTbkxpA+k/tLgcM9OPTFLM2cLz9I70L4tMK/Bgwl75hmyly7l/Pz5+EycQMjDD+PZ103ieW3ExX/8A8MJ+2WoJVZflmJ8ugMAACAASURBVDQaERoNQqdrkKCki+tD2BNPNOir0uurHMdZWLTaditDvXTpUlauXEl8fDyvvvoqAQEBigy1gvvJ/McLmEtLCf/bcwi149stToeOGkqgMMVuR/GaYxnkllZy+2j7ncTN4eep5cM74rl1VFc+2HGWbaft3ybaeyaX7GJD3YPVmkOGuts5NUVqHNgeAgi4+SZ6rF9H8KN/oPzwEc5edz2pDz+CITnZoXEuFyQSWVFhNQJabaNGoDlqO45NOTntUob6gQceIDk5mcOHDxMeHs6jjz7abH97x3UUZUXwG6Nk61aKfviBoEWL0PV0wtmL9QYX7BWMv67pPdpGqa7qZaej+ONd54kN9mFsd8crmDWFRq3ir7P78XNCFsu2JDGxl+0s9aSsYm5+fw8LJ8Ty56viLp2ofh85iRA5rOZwN79uaFQahw0BWMNJg+69l4CbbybvoxXkrVhB8YYNaMJCEQ7cBrVdutD1g3/XRNO0Fxp7cm8Kw7lzWEpL0UZEoAlwMDqtimrHsSU1tcFvrz3IUIeGhtb8/95772X27NkAl5cMtUL7wlxSSsYzz+LRozud71vYZDuTxcSejD1YGpd9cl5jqMYQ2E4mO5ZayOELBdw2qhsqletPPLXx0Ki4Z3wMe87kcTAl32b7f/18Gikho6Ci7omaamV1Q0i1ai3RftEOhZDWR63XE/zQ/9H95w0E3X8fPiNG4j1ihF0vjx7dKdu7l7IDB52+fltjMRqxlJSgCQ522ghUowkORmg0pFy4wK5du4D2I0OdkXFJsvzrr7+uGb+1Zajb1+OCglvJ/te/MF28SLfPPkXVTBWyLRe28Pstv2d+n/k8PuLxOk9JFmnhbOFZru95vRMTOAkqDQTazg5eufsc3h5qrh/mnrKT80Z05a1NSSzfksx7t8c32e7kxSJ+PGb9smYV1zMEATHWyKFGahP07NSTA1kHqDRX4qF2XlldExhI8EMPOdTHUlbGqTFjKV6/Hp+RDW8+HQFLUREAar9GpcccQqhUqAMD6RMby8cffMD9999Pz549eeCBBwgJCWHLli1MnTq1jgx1eHg4Q4cOrRPps3PnTp5++mm7rllbRlqj0TSQob7//vuJj4/nscce4/DhwwghiI6O5t1337XZf9myZdx5552Ul5dz1VVXtW8Zane+FBlqxyk9eFAm9ImTGX/7u822r+5/VfZf0V/2X9FfvnfkvTrnUopSZP8V/eWXiV86PonPbpZy6QibzfJKDLLXk2vkE6uPOn4NB3h13UnZ7U8/yNOZRU22uf+T/bLfX9fKee/tlpNe2dywwTtjpPzPDQ0Orz+3XvZf0V/et/4+WWYsa8FZ28eFxYvlqQkTpcVsbvVr18cZGeqKM2dk+alTLTaHM8nJsm/PnrLi9GlpsVhqjisy1O6XoVZoJ1gqK8l46ik0YWEEP/KIzfYJOQnEBcYxO3Y2bx56k69OfVVzzmWNoSDbVcn+t/8CBpOF20dHO34NB7hjTDSeWhXLtzYub3U8vZCffr3IgnEx9ArVN3QWgzUfoqhh1Ma0btN4bsxz7ErfxYM/P0ip0U2S2E2gnz4dU2YmFUePtup1WwJpMmEpLW2R1UA1QqVCqNVYKipqVhugyFBXoxiCyxxzQQEXn32WyqRkwp95GrWvT7PtpZQk5CYwIGgAz419jrGRY3luz3NsStkEXIqEcdgQGCsg/6xN/4DZIvlkz3lGxATSO0zv2DUcpLOvjpuHd+WbQ2mkF5Q3OP+vn0+j99Rw97gYgvU6iitMVBjrJQX5hkJx46Upr+15LS+Of5FDWYdYuGEhRZVFjbZzB74TJ4JWS9GGDa12zZbC3ILbQtVER0dzLCEBodNhzMqqE32jyFArhuCyxZSfT9Zrr5M0eQqFX60m8I7brTcHG1wovkCxsZi+nfuiVWl5beJr9Ovcj8e2PcaBzAMkFyQT6h2K3sPBm3ReMkiLzYih/efySM0v5zYHdIVc4Z7xMUjg39vP1jl+LLWQDQmZ3Ds+Fn8vLSFVhXCyiuqtCvRhUJYLpsaTx2bGzuTVK14lITeBe9bdQ36Fbed0S6D288Nn9CiK129oF5LMjszBXFSE8PBAOJjoaAshBNrQUKTBgLmgoEXHbm84+jdXDMFlhikvj6xXXyVpylRy338fn4kTiPnuW0L//Ge7+h/PPQ5Av6B+AHhrvXl7ytuE+4Tzfxv/j18u/uJ8DQKwaQgOVEXxjO/ZciGjzREV4M01gyL4fF8K+aWXbuav/3wKfy8td42NBiDEz3pTauAw9q0K/yttWu56StcpvDX5Lc4UnuGutXeRXdY6Ffb006ZhvHABw0n7E7jcgaenJ7m5uXbdnGpvC7VEfHx9VHo9Ki8vTFlZl20xICklubm5eDpgSJWoocuI3A8+IHvp28iKCvxmziTo/vsczhU4nnMcD5VHna2fAM8A3p32Lrf9dBuZZZnMiJ7h+OSyE0GooHPzRuTg+QJig33o5O18pI2j3DexO6sPpfHx7nM8MrUXh1Ly2XQyiyUzeqP31AJcWhHU9xPoq4TwijPBv+kIp3GR41g2dRmLNi7izrV38sWcL/DWupYtbQv9lClcfPoZijdswDMuznYHN1Ed+25PiVlLWRnmggI0JhMi3z2rJ0uFAXNeLuriYlQ+zW+VdlQ8PT3rZCDbQjEElwmmnByyXnkVnzFjCH3yCXQ2RLSa4njucfoE9kGr0tY5HuEbwfKpy3lo00OMCh/VRO9myE6EgGjQejXZRErJ4Qv5TOwV4vj4LtA7TM/UuBA+3nWOhRNief3n0wT6eHDHmOiaNpe2huqtCPRVK4Im/AS1GR42nJfGv8RDmx9ib8ZeJnV1rEa0o2gCA/GOj6do/XqHQ1BbEq1WS0yMfQWFLjy4iIoTJ+ixaaNbVgRg/Zyl3H4HhrNn6bF+HSpv9xrkjoCyNXSZULJlC0hJyB8fddoIWKSFhNwE+nZuXN+mZ0BPfrr+J8ZHjXd88OxEm9ISqfnl5JRUMqSrfcVnWpIHruhOfpmRx786xrZT2dw3IRZf3aXnpABvDzQq0XBF4Btm/bfkol3XGRM5Bq1Ky4HMAy019WbRT59OZVJyTRWv9oy5pJTSHTvQT5vqNiMAVl9B8O8fwZyTQ14j2cHNUbZ/Pxf/8Y/LroiQYgguE4o3bkIbEYGuj/Mln88VnaPMVFbjH2gxzCbITbLpH6jO8m0LQzCsWyAjogP57kg6Qb4e3FZP30ilEgT56hoaAp9gQFi3huxAp9YxIGgAB7NaJ+tXP20qAMXr23/0UOm2rcjKSvymT3f7tbyHDsV34kRy//1BTZSSLaTZTMZTfyV/5SekL3kM2YSsdEdEMQSXAZayMkp37cJ38mSXnqSO51gdxU2tCJwm/yxYjDYNwaGUAry0anqHujdstCkemGT1i9w/sTveHg13TUP8dA1zCdQaqzGwc0UAMCx0GAm5CZQZy1yarz1oQ0PxGjSI4g4QRlq0fgPqoCC8ail+upPgRx7GUlhI7r8/sKt90ZqfqDx7Ft/Jkylet46Lf/97u4jIagkUQ3AZULp7N9JgcLgAfX0SchPwVHsS6+/c1lKT2BkxdOhCAQOj/NGo2+ZjOal3CN8sGsuCsY3vZ4foG1kRgDWE1M4VAUB8aDxmaeZw9mFnp+oQ+unTqDh+nMpU1+WK3YWlooKSbdvQT53ilCKuM3jGxeF39RzyVqygspaQW2NIs5mcZcvQ9exJ1NK36HzP3RR8voqct99plbm6G7u/cUKIECHEtUKIRUKIBUKIEUIIxZC0A4o3bkKl1+Md37Rmjj0k5CbQJ7APGlULxxBUi7I1k1VcYTSTkF7IkK6uCYy5yuAunZoUuQvWe5JdP3wUqgyBbWdxNYNCBqEW6tbzE0ybBkDxz+13VVC6YweyrKxVtoVqE/Loo6BWk/XPl5ttV7R2LZVnzhC0aBFCpSL40Ufxv+46cpYuJf/zz1tptu7D5o1cCDFJCLEO+BG4CggH+gJ/AY4JIZ4VQrRcCqCCQ0izmZItW/CdMAGh1dru0ARmi5kTeSec9w9knWz6qTg7EfyiQNf0ls/x9CKMZtkm/gF7CdHryC2txGSuF3/uGwol9q8IfLQ+9Ans02qGwKNrV3R9+rRrP0HR+vWo/f3xHj68Va+rDQ0laOG9FK9fT+mevY22qb0a0E+3GlUhBOHPPYvvpElcfO5vFLVAuci2xJ4n+pnAvVLK4VLKhVLKv0gp/yilvBoYBBwCprl1lgpNUn74MOa8PJe3hc4WnqXcVE6/zk4ags9vguVjIa0RJ2j2STv8A1WO4i7t2BD46ZASckrqRYzow6A0Gyz2Ow+HhQ7jWPYxh8pZuoJ++jTKDx3CmNV04ltbYamspGTTZnynTHHpYcZZAu+6C21EBJkvvIA0mRqcL163jsqkZIIefAChunTLFBoNka+9itfQoaQteYzS3btbc9otik1DIKVcIqVstMirlNIkpfxGSvlVY+cV3E/xxk2g1eIzYYJL49RkFDtjCCxmKLhgvRmumAWJP9U6Z4Gc03b5ByI7edVk8LZHQvRNZBfrw6zyGaX2ZwwPCx1GpaWSX3N+bckpNonf9OkgJSUbN7bK9RyhbPduLCUlNU/brY3K05OQxx7DkJhIwZdf1jknLRZyli3Do0d39DMaJlKqvLzo8s7b6KKjSV20GFOO/VXv2hOO+AgeFkL4CSsfCCEOCiFad0NPoQ6y6ovtM2IEal9fl8Y6nnscb4033fyc0PgpyQJphomPW0XlVs2Hfe9bzxWmgKncpiE4nFLQrreFAIKb0huqziVwwE8wNGQoQOttD3XvjkdMTLuMHipavx6Vry8+Y8a02Rz0M6bjPXw42f96A3NhYc3x4vXrMZxOIvjBB+usBmqj9vcn7JmnsZSVUX7sWGtNuUVxxNm7QEpZBEwHgoG7gBfdMisFu6g8c4bK8+fxdXFbCC5lFKtVTkRsFFnruBIxGO78AXrOgDV/hPV/uVS0pRnV0cyiCtIKytvcUWyL6uzi7JJGhOfAocihTp6d6NGpR6sZAiEE+unTKd27D5ObpBucQZrNlGzchO+kSc0WS3I3QghCn/gz5qIict6xRgJJi4Wct9/Go3vjq4Ha6LpbQ48rz55z91TdgiOGoDqUYibwkZTySK1jCm1A8UarNLR+smuGwGQxkZiX6LyjuLjKEOjDwcMHbv4Uht8Lu96C7x+2nmsmYuhQilUJsr2vCIJ8m1oRVMlMOJBLANbtocPZhzFZGu5LuwP9tGlgNlOyZWurXM8eyo8cxVxQgH6ye+U27MEzLo5Oc+eS9+lnGJKTKV6/AcPpJIIeeMBmSKu6UyfUAQFUnj3bbLv2iiOG4IAQYj1WQ7BOCKEHLk/5vg5CycaNePbrhzYszKVxkguSMZgNzjuKq1cEfpHWf1VqmPkyTP+7NZrGJxi8A5vsfiglHw+1in4R7Tv4zEOjItDHo2kFUgdWBGA1BKXGUhLzE203bgE8+/VFEx5O8c8/t8r17KFk+zZQqdp0W6g2wY88jMrLi8wXXiTnnXfwiI3F76or7errERPzmzAEdwOPA8OllGWAB9btIYU2wJSdTfnRoy2yLZSQmwA46SgGqyFQacG786VjQsCY/4Nbv4I5bzbb/VBKAf0i/dBpWieRyBUaTSrTeFjfu4Mrgho/wcVW3B6aMoXSnTuxlDcsxNMWlG7bjtfgwaj9/dt6KoBVqC9o0YOU7tiB4dQpu1YD1XjERGM4d86t83MXdhsCKaUFMAEThBDXARMBJ4TpFVqC4iqROf2UKS6PdTz3OL5aX7r6NV+lqUmK0sEvHBpzpvWYCn1mNtnVaLZwNK2AIV3at3+gmuCmsot9w6DYMUMQ6hNKF32XVvMTAOinTkVWVFCyY0erXbMpTLm5VBw/ju8EJ0QM3Ujg/PlW53qP7vjNtL8wvC42FnNODubiYjfOzj04EjX0IfAhcD0wp+o1203zUrBBycZNaCMj0fWyXQPYFsdzjtO3c19UziaKF6WDPsKprokXi6kwWtq9f6CaEL0n2fWlqMEqR+2gIQDrquBg1sFW06zxjh+G2t+fknawPVRaZYx8xrWMIdiZtpNFGxdRaCi03bgZhIcH0f9dRfSnnzokd+FRJbXdEbeHHPnmj5JSxksp75BS3lX1WuC2mSk0iaWsjNLdu/Gd4prIHIDRbCQxP9E1obnidPBzzhAcakPFUWcI8dORXWJoeOP2DXMou7iaYaHDKDAUcKawdWSihUaD76RJFG/ZijQaW+WaTVGybTvqzp3x7Ot60Zwj2Ud4ZPMjbEvdxpenvrTdwQZqX1+Ht6s8oq2GoCNIftfHEUOwWwjh0N1CCHGlECJRCJEkhHi8iTZXCCEOCyGOCyHaTzhDO6Zk506ryNxk17eFThecxmgxOu8fkLJqa8hZQ1BAsF5HZKemC9a0J4J9dRjNkvyyejdRfZUhcLD8YXyoVR+qdbeHpmApLKRs//5Wu2Z9pNlM6Y4d+I4b12R8vr2cKTjDoo2LCPEOYWDQQD47+RlGS+sbOY8uUaDRdMgQUkf+Ah9jNQaJQoijQohjQoijTTUWQqiBt7HqE/UF5tU3JEKITsA7wNVSyn7ADQ6/g98gJRs3ofL3xzt+mMtjuZRRDFCeD6YK5w3BhQKGdOnk1kIkLUmIX1UuQYOSlWFgMUF5nkPjRemjCPEKYX9m692UfcaORXh6Uvxz22UZVxw7hrmwEB8X/QMXSy+ycMNCtCot7057l3sH3ktWWRY/n2/9rS+h1eIRFXXZbw19CNwGXMkl/8CcZtqPAJKklGeklJXAKuCaem3mA6urJSyklO1PCKWdIU0mq8jcxAkIjesqoQm5Ceg99ETp7a9vWoea0FHHDUF+aSVnc0rbfSJZbZqUmfC1v2RlbYQQDA0dyoHMA63mJ1B5eeEzbizFGze2mZ5+yfYdLoeNFhoKuX/D/ZQaS1k2dRlR+igmRE2gm183Pkn4pE3eW0cNIXXEEKRIKb+TUp6VUp6vfjXTPhK4UOvn1KpjtekFBAghtgghDgghbndgPm6nIjGRtD88isXQSJRIG1GRmIi5oADf8a5pC1VzPOc4/Tr3c/6JvPrG54Sz+PCFjpFIVpuQpmQmnMgurmZY6DCyyrJIK2m9egH6qVMxXbxIxa+to3VUn5Lt2/EaMABNgHMPAeWmchZvXExKcQpvTn6TPoHWzHWVUHFL3C0cyznGkewjLTllu/CIiaHy/PkOV73MEUNwUgjxmRBinhDiuupXM+0bu7PUN9EaYBgwC5gBPCWEaDQMRgixUAixXwixPzvbfnEvVyjZvIWiNWso27OnVa5nD4aT1uQjz/6ul5M0mA2cLjjt/LYQQFHVzcuJFcGhlHxUAgZGtY8Ycnuo3hpqWLvYuexisBoCaGU/wRVXgFrdJttDprw8Ko4dc3pbyGQxsWTrEo5kH+GlCS8xPKyudPU13a9B76FnZcLKlpiuQ+hiY5CVlRgzHFsZtjWOGAIvwIBVa8ie8NFUoEutn6OA9EbarJVSlkopc4BtWKWtGyClfK8qaik+ODjYgWk7T/Ufs3jT5la5nj1UJJ5EeHnh0dXJmP9anM4/jclicq1GcVEGIC49ETvAoQsF9Anza7QsZHvF20ODr07TuAIpOBVC2r1Td/x1/q1qCNSdOuE9fHibZBmX7twJUuLrhGKulJLndj/H1tSt/GXUX5jWraFiqbfWm7m95rIxZWOrrrKg44aQOpJQdlcjr+bCR38BegohYoQQHsDNwHf12nwLjBdCaIQQ3sBI4ISjb8JdGDOsdqtk8+Z2U5vUcDIRXa+eLVLOr1oC2aXQ0aI08A0BtWM68haL7BCKo43RaHax1gs8/Z0yBCqhYmjI0FY1BAD6KVOoTE7GcKZ1b1ol27ajDgjAs5/jDyDLjyzn66SvWThwITf2vrHJdvP7zEcg+PxE61YPqzEEHSyE1JGEsighxNdCiCwhRKYQ4ishRJMeRimlCVgMrMN6c/+flPK4EOJ+IcT9VW1OAGuBo8A+4N9SyrbZtGwEU0YGQqvFlJVFxfGEtp4OUkoqEhPx7N20kqcjY311+iu6+XUjwse5iB+gKpks3OFuydklFBtMHcpRXE2wXkd2fR8BVOUSOG4IwLo9lFKcQnZZ62x7gjWMFKB4Y+utCqTFQumOHfg4ETa6+vRq3jnyDtd0v4bFgxc32zbMJ4xp3abx1emvKDWWujJlh1AHBKDy98dwua4IgI+wPtFHYHX6fl91rEmklGuklL2klN2llM9XHVsupVxeq83LUsq+Usr+Usp/Of4W3Icx46JVflYISja3/faQ6eJFLIWF6Po0r+1vD5subOJk3knuHXCva6GbxRmXxOYcoKMojjaGVWaiqexix53FcMlPsPTwUipMjYztBrTh4Xj260dJK/oJKo4fx5yf77CsxPbU7Ty3+znGRIzh6TFP2/WZva3vbZQYS/gm6Rtnp+swQgh00dEdLpfAEUMQLKX8qKoqmUlKuQJrXYLLEnNxMZaSEjz79sVryBCKN29q6ylRcfIkAJ59XFsRSClZfmQ5XfRdmBU7y7VJFaVZdYYcJDmnBA+1ipjOPq5dvw0I0Xs2zCMA68rIyRVBv879uDXuVlafXs0N39/AsezWKXCinzaV8iNHMGa2TuR2ybZtIAQ+48bZ3ed47nEe3foovQJ68doVr6FV2bcNOTB4IIOCB/GfhP9gdqCMqKt0xBBSRwxBjhDiViGEuup1K5Drrom1NcZ0q6NYGxGO76QrMCScwHjRuS95S2GoMgS6Xq6tCDZf2MzJvJMsHLgQjcoFR21lKVQUOhUxlJpfTmSAFypVx0gkq02In47SSjOlhnp1BHyrVgRO+JOEEPxpxJ94f/r7VJgruO2n23jr0FsYze7NkNVPnQpAyabGVwWWiooWDYUs3bYdTwfCRlOLU1n08yICdAG8PeVtfLSOPTjc1vc2UktS2ZK6xYnZOodHTAymrCzMJa23JeUqDlUoA24ELgIZwFwuYxnqakexNiyspvBLW28PVZxMRNulC2pf55+ia68GZse6qBlYVBUi58TWUGpeGVEBHUNWoj41uQSNZRebDdZsaycZFT6K1VevZlbsLN47+h7z18zndP5pV6bbLB7du+PRrVudMFJZWUnxpk2kPvQwp4aPIPOll1rkWqb8fKt0+nj7toUKKgp44OcHMFqMLJu6jGBvxzcgpnSdQoRPBP9J+I/DfZ3FI7bKYdyBJKkdeRzsIqW8uvYBIcRYoNHC9h0dU9XTvyY8Ak1IMNpuXSnevJmAefOa7COlJOett/CdOBGvQY1GwbqE4eRJPF30D2xN3cqJvBM8N+Y511YDULcymYOk5pczPaLj5A/Upia7uKiCmKBaRrkmlyCz2UI8ttB76Hl+3PNM6TqFZ3c/y00/3MSAoAEO+XKifKP429i/2ewjhEA/bSq5Kz6mdPduijf8TNGaNZgLClB37oxKr8eQ0DKBfKW7dlnDRsfbty301K6nSC9J5/3p7xPbKdapa2pUGubHzeeV/a+QkJvgWoScnehqQkjP4NUC+T6tgSMrgrfsPHZZYEzPAK0WTXCQ9ctyxSTKdu/BUtr0cq94/QZy3llG5ssvt/h8LGVlVKakoHMhYkhKyTuH3yHKN4rZ3VtAQbx+ZTI7KTWYyC2t7LgrgqaSylzIJWiMyV0n8/U1X3N196tRq9SohMquV3FlMd8mf8up/FN2Xcd3yhQwmUi5awEFX32Fz5jRRC1fRs8tm/EZNxZjev30H+co3bYddadOeA4YYLutsZQdaTuYHzefoaFDXbrudT2vw0vjxWcnPnNpHHvRdu0KKlWH8hPYfCQUQowGxgDBQog/1DrlB7T/klJOYszIQBsaWhPi5jt5Mnkff0zJrl34TWuYxCKNRrJfew3Uasr3H6Di5EmXnbq1MZw6BVK6tCLYlrqtZjVgr8OtWWqyih1bEaQVWKtjdVhD0OTWUNXvwQk56qYI9AzkmTHPONTnYulFpn05jd3pu+kdaPvz4jVoEEGLF6MNC0U/YwZqvb7mnDYigqLMNUiz2aXcFWmxULJjh1Xwzo5x9l/cj8liYlyk/U7lptB76Lm6+9V8ffpr/hD/BwI9nV+t2YPKwwNtVFSHCiG1Z0XgAfhiNRr6Wq8irH6CyxJjRnqdWsDeQ4eg8vOjpIks4/wvvqDy/Hki/vE8Qqcj/9OWffqoqJKW0DlpXKSULDuyjEjfyJZZDYDVR+Dpby1Y7wAX8soA6BLo3TLzaGX8vbR4qFUtJjzX0oT5hNHdvzu7M3bb1V6oVAQvXkSnuXPrGAEAbXgEmM2YslyLKjKcOoU5NxcfO7eFdqbvxFPtyZCQIS5dt5r5feZTaankq1Nftch4tvCI6VghpDYNgZRyq5TyWayFaZ4FXgdek1K+JqV0nxerjTFlXEQTcelJV2i1+E6YQMnWrQ2iKMwlpeQsfRvvESPwu/pq/ObMpvCHHzAXulYpqTYViSdR+fqijXTcMQuwPW07x3OPs3DgwpZZDYDTlclS8zv2ikAIYU0qq78i0PmCh6/TuQQtyeiI0RzIPOByToI2wvr3dVU7x5CcDIBnnH179LvTdxMfFo+H2sOl61YT2ymW0eGjWZW4qlVqFeiiY6g8dw7pYH2KtsKhPAIhxDGsWcDHhBBHhBCuC+K3Q6TZjDEz0/o0VAv95EmY8/IoP1K3DEPehx9gzssjZMkfEUIQOH8+srycgq+/brE5GU4mouvd26nkLyklyw5bVwNzujenHO4gTlYmS80vQ6dREeyra7m5tDKNGgKwrgqczCVoScZEjMFgNnAw66BL42irHoaMaa75CYwp1pgSj65dbLSEtJI0zhWdY0yE8xLVjTE/bj5ZZVlsSnF/TpBHTAyyoqIm6KS942g9ggellNFSymhgETYyizsqppwcMJnQhtcVUvMZPx40mjphpMasLHI/WoHfxXmX3QAAIABJREFUzKvwqnKCefbti9fQoeR/9nmLPBFIiwVDYiKevZ3zD2xP286vub9y74B7W241AE5XJruQV05UgFeHKUbTGCF6XUMparD6CdrBimBY6DC0Ki270+3bHmoKbXiVIXBxRVB5PgVNaCgqL9urwOo5t7QhGB85nijfqFZxGleHkHYUP4EjhqBYSrm9+gcp5Q6guOWn1PaYqj701V+CatR6Pd7D4+tkGecsfRtpMhH8yCN12gbcMh9jSkpNgW5XMKamYikrc1hawmwxs+rkKh7f/jiRvpFc3f1q253sHtwIJVnOrQgKyjqsf6CaEL9mZCbawYrAW+vNkJAh7Erf5dI4Km9v1J06YUx3TcWzMiXFbsXcXem7CPUOJdbfuZDRplCr1Nzc52YOZh3kRK57tS1rQkhbWdDPWRwxBPuEEO9W1RieKIR4B9gihBgqhHAtvqudUf30owlvGA2jnzSJyqRkKlNSMCQnU/DllwTcfHODD7nftGmog4PI+/RTl+fjjLTE4azDzPtxHs/vfZ6+nfuyfOpytA4qhDZL8UVAurQi6MiE6D3JLzNSaaq34vMNs/5u2oFa7eiI0ZzKP+WykJ0mItz1FUFKCtputg2B2WJmT8YexkSMccuK8dqe11pDSU+6d1WgDgpC5evbYUJIHTEEg7FWFHsaeAaIwxpW+irwSovPrA25JC/R8CbnO2kSYM0yznr1NVTe3gQ9cH+DdsLDg4Abb6J023YqzzdXyM02hpOJoFKh69nTZtvc8lye2vkUt/10G7kVubwy8RXen/Y+0f7RLs2hAU5WJiuqMFJYbiQqoIOvCKpCSLNL6oeQhoKxDAxtv1iu3lrZk+FaYSVtRASmdOcNgbmkFHNODh5du9ls+2vurxRXFrf4tlA1fh5+zImdw5oza8ivcD4D3BZCCKvm0LnLzBBIKSc185rszkm2NsaMDFS+vqh9fRuc8+jSBV3PHuR+tIKSTZvofO+9aAIbj0vudOONoFaT//kql+ZTkZiIR3S0zf3V75O/Z843c/gh+QcW9F/A97/7nhnRM9yzF+9kZbK0qoihLh3dEFQnlRXVDyGt8iu1YC6Bs/QJ7EOALsDl7SFteATG9HSna3IYL1Q7im2vCHal70IgGBk+0qlr2cP8uKpQ0tPuDSX1iInG0EFCSB2pR9BJCPGQEOI1IcSb1S93Tq6tMF7MaOAfqI3vpMmYLl5EExpK4O23NdlOGxqC3/RpFKxejaWszOn52CMtUWos5a+7/kqMfwxfXf0Vvx/2e7y1brzZOlm0vjqH4HLYGgL3Zxe7gkqoGBUxit3pu10qrKSNiMBSWoql2LlVTuX5KkNgx9bQ7vTd9O3clwBP99Wp6N6pOyPDR7Lq5CpMFpPtDk6ii4nBlJHh0ne/tXBka2gNEA0cAw7Uel12mNIz6uQQ1Ec/fToAwQ8/bPMpPeCWW7AUFVH4ww9OzcVcXIwxLc2mtMS+jH2YLCYeGfqI07osDlGUDhpP8HLsC1udQ9DRncXB1VtD7dgQgHV7KLci1265icaoiRxyUmqisip0VNuleUNQXFnM0eyjbtsWqs0tfW4hsyzT4VBSo8VIZql9q72aamUubg23Bo4YAk8p5R+qahJ8XP1y28zaEGNGBtqwpg2BV/9+9NiymU7XXWtzLK+hQ9H17k3+p5859VRmSKwqVm9jRbAjbQfeGm8GBw92+BpOUV2ZzMFtpwv5ZXh7qAnwbkHHdRvQ2ccDIVq2iL07GB0+GsClMFJtZFVSmZN+gsqU86iDgmyq5u67uA+zNLeKIZgQNYFI30g+PeFYMMeLe1/kqtVXcST7iM22Hal+sSOG4BMhxL1CiHAhRGD1y20zayMs5eWY8/Ob3RoC6shPNIcQgv9n77zD26rOP/45kizvEccjtpN4Ze89SUJIgIQRNjRsCmUVKC1tWWW2rAIlrDLKTCCMEAJh/cIme09nONs7HrEdb2ud3x9HcjxkW5Il27Hv53n0SL6641xZuu897/r2uOpKatPTqd7i/gTKldYSUkrW5K5hUtwk72YGtYSHymTZJad+DQGAQa+jZ7A/hY1TSAPC1Uypk8wIYoNjSQ1PbVOcoK0zAnOGa6mj63LXEWQIYmS09zv3Nkav0zN/0Hy2FmxlX/E+l7Ypqi7ii4NfYLaZ+csvf2k1G8uYmAhCtLsmtCe4YwhMwLPAOk66hTb7YlAdiUN8xq8F15C7hJ93HrqwMEoWu5+yVrNvL/qICAwxMc2uc6TsCDkVOUxNmNqWYbpHWY6HqaNVp3yg2IHTojIhlHuoEwSLHbS13YS+Z0+E0Vin0eEurtYQrM1dy4ReE9rtZubCfhcSaAjk3TTX6mI/Tf8Uk83EczOeo9xczl9+/UuLwkG6gAD84uO73IzgL0A/e2Vxsv3RDs7o9qW5YrK2oAsKIuKiCyn74UdVtewGtfvS8R80qMU76DU5awC80qnRJWw2dcfrZtdRKSU5Jad+DYEDVVTWjIh9J5kRgIoTmGwmtuZ71m5CCIFfXJxHMwJbdTWW/PxWA8VZZVlklWcxJcH3biEH4f7hXDnoSr478l2rMZQaSw0f7/uYGb1ncHbS2Tw+9XG2F27nqY1PtbjdqSJb6Y4h2A10/vB3G2mpmKwtRFxxBZjNlH7uev8habFQe+BAq60lVuesJiU8hfgQ9+/QPaLqOFhNbruGyqotlNdaTvlAsYOYFkXsO48hcLSbaJN7KMGzWgJTVpbavpUZgWNs7REfqM8Nw24gxC+El7e1LK3y1eGvKKkt4bqh1wEwJ2kONw67kSX7l7Bk/5Jmt1O1BEfblLXVHrhjCKzAdnt1cZdNHzXn5oEQ+LXgivEE/5QUgsaPp/TTT13uP2TKyEDW1rYYH6i2VLP52Ob2dQt5qEyWVdI1UkcdxIQGUFRhwmpr9CMP6VyuoSC/IMbEjHG5LbUzDB7OCE42m2u5mGxt7loSQhLoG+paGwpvEe4fzvXDrufXrF/ZXrDd6To2aWPh7oUMjhzMuNhxdcvvHH0nUxOm8uSGJ5vd1pichK2qqs1tvH2NO4bgC+AJYC1dOH3UfCwPQ3Q0wuid9rf1ifjdFZizs6lc49qd2cnWEs3PCDYf24zJZmo/txB4rEyWXWcIusiMIMwfq01SXGlq+EZoLNSWganzTKDb2m7CLy4eS2EhNpOp9ZXrUVdD0ELXUbPNzMZjG5kcP7lDkgiuHnw1kQGRvLTtJad37quyV3G07CjXDb2uwfj0Oj3PTHuGuOA4/vzrnymoanqx9z9FMofcMQTHgUX1U0e7YvqoJS8PQ5xrGUHuEnrmmegjIyn99BOX1q/dlw4GA8bU1GbXWZ2zmgB9AGNj27EjuIfKZFnFXaOq2IGjjXbTWgKHUlnncQ9NjldppJ62m3C0W3G3rbIpMxN9RAT68Ob1qdOK0qgwV7S7W8hBkF8QN4+4mU3HNjmdNS3cs5DYoFjOSjqryXvh/uG8OPNFKs2V/PnXPzfROjCmqDCqQ4+hs+KOIfgdcEAI8W8hxGBfDaijMefmNdEh8BY6o5GIiy+i/OdfMOe3PlWsSd+Hf0oKuhZmJ2ty1zC+13j89e3Y278sD4T+ZM68i2SXVBHqbyAssFWF1FOCk9rFzSmVdR73UFvbTdTpErgZJzBlZrTabG5t7lp0QufTthKtcdmAy4gLjuOlrQ1nBXuP72XjsY1cPfjqZlu49+/Rn/sn3M/Owp3sLGyoVWKIiVHN5w4d9un424o7vYauBkYDh4B3hRDrhBA3CyFCW9n0lEFKqYrJvBwork/E5ZeD1Urp0s9aXVdlDDXvFsoqyyKjLKN93UKgXEMhsaBzT8M2q6Sa3pFBp3wNgYPW20x0rGRlfeq3m7BJ9zUy6pTK3IwTmDIyWo8P5KxleNRwwoxhbo/LWxj1Rm4beRu7j+/mp8yf6pa/v+d9ggxBXDLgkha3dxixQ6UN7/yFEPinplJ7uIsYAgApZRmwFPgYiAMuArYKIe70wdjaHWtpKbK21qeGwNi3L8FTplC65LMmkpf1sZSUYCkoIKCF1hKrc5XWQbsbgjYok3WVQDG00GaiEzWeq4+j3UR6cbrb2xrsBZTu1BLYamux5B1rsYagtKaUtONpHeYWqs/5qeeTHJ7My9texmqzcqzyGCuOrODi/hcTamz5frdXcC8CDYEcPtH0gm9MTaX20EFfDdsruNN07nwhxDLgZ8APmCClnAuMBP7qo/G1K467HV/FCBxE/O4KLHl5VKxc2ew61dtVFkJLM4I1OWvoE9qHvmHtm2nhiTKZlJLskuouEx8ACPDTExZgaNqBNCgSdH6dKoUU1A2DQPBL1i+tr9wIndGIITrarRmBOTsbpGyxhmBt7lps0sa0hGluj8nbGHQG7hh1B4dPHOabI9+weO9ibNi4esjVrW6rEzpSwlOazAgA/FNTsRYWYS0t9cWwvYI7M4LLgBeklCOklM9KKQsApJRVwO99Mrp25mQxmW/z8UNnzkQfHUXpx86DxjV795J3/wMYYmMJHOm8d1CttZaNxza2/2wAVIzATUNQXGmiymTtUjMCgJiwgKauoU5YXQwQFRjF6JjRDVwf7mCIj6v7jbjCyYyh5g3BqpxVRAZEMjRqqEdj8jazE2czOHIwr257lc/2f8aZiWeSEOJadlxqRCqHS5vOCPz7qWSPzuwecscQ3AasBhBCDBBCzBNC+AFIKT37ZnUyzHneby/hDOHnR8Sll1KxciXmnIYSgNVpu8m4/gZEUCCJixY226hra/5Wqi3V7W8IasrAVO62IegqXUcbExvmz7HGMwJQMZROFCNwcEbfM9hfsp+ssiy3t/WLj3crWGzKVF03/RKdxwisNitrctYwNX4qOuGWl9pn6ISOP435E7mVuZSby7luyHUub5sakUpBdQFlprIGy42p/QCoPdh53UPufPorgQAhRALwE3AD8J4vBtVRmPPyEP7+6Hv4rhe6gx6XXgpAyWcng8bVO3aQecMN6ENCSFy4qMU7qdU5q/HT+TUocPEa5cdg2W3Os148VCbrasVkDuLDA8krdVZd3KtTZQ05mNV3FoBHswK/uHjMeXkuV8maMzPRhYWhj4hw+v7u47spqS3pmFltC0yJn8KU+ClMjpvM8OjhLm+XGq7u/BvPCvzi4xCBgZg6cQqpO4ZA2N1AFwMvSykvAob4ZlgdgzkvF79evdolq8UvIYGQ6dM58dlSpNlM1dZtZP7+RvQ9epC4aCHG3i1PR9fkrGFc7DjfiM/s+wZ2LIbPfg/WRsIdHiqTOWYEXc4QRASSX16D2dooEye0F5zIhmrfySF6Qu/Q3gyKHOSZIYiPR9bWYi0udml9k73raHO/p1U5q9AJXacIFNdHCMF/Z/2X12a/5tZ2Dh2QJplDOh3+KSnUHuwihkAIMRm4CvjGvqzFhHAhxBwhRLoQ4qAQ4r4W1hsvhLAKIS51YzxepzVBGm8TccUVWAoLKXhhAVk33YQhKorERQudaiXXJ68ij0MnDvmurUTuNtAZIGM1/PKvhu/VVRW79zlll1QREeRHaMCprUPQmISIQKSEYycazQpGXAHWWvj4KjB71vXTV5zR9wx2FO5wu8r4ZC2BawHj1rqOrspexYioEUQEOJ8xdCR6nR69m+nR8cHxBOgDOHTCScC4X2qnLipzxxDcDdwPLJNS7hZCpADNph8IIfTAq8Bc1MxhvhCiyQzCvt4zwAp3Bu4LzMeO+TxQXJ+QGdMxxMVR/M47GHr1ou+ihfjFtl6k5fO00bztkDwDxl4Pq1+A9O9OvlfmoWuouOt0Ha1PfIQ6p9zS6oZv9JkAF74GGWvgi1tVx9ZOwuy+s5FIt7OHTuoStB4nkCYT5pycZovJiqqL2H18N9N6d3y2kLfQ6/Qkhyc7DRgbU/thOXYMa0VFB4ysddwpKPtNSjlPSvmM/e/DUsq7WthkAnDQvp4JVXtwgZP17kTVJnRoVyZpNmMpKHBZcMYbCL2e6DvuIGjSJBIXvu9yo7s1OWuIC44jJdwHXcDN1VCwF+JHw5xnIG4kLLsFSo6q98tyIDAS/ALc2m12SdfRIahPfIT6HHIaGwKA4ZfCmY/D7mXww0PtPLLm6RfRj76hfd12D7lTVGbOzQWbrdliMkeFc2dIG/UmKREpzc4IgE4bJ2i11l8Ica2L+9oupaxfX50A1E9NyAYa1JDbA88XAWcA41sZx83AzQB9XRC5cBdzfgFI6fOMocZEXHIxEZdc7PL6ZquZ9XnrmZs81zexjPzdYLNA/Ch1sb98IbwxHT69Fn7/vUfKZI4agjMGebeja2eg2RmBgyl3qVjBulcgvDdMuq0dR+ccIQSz+s5i0Z5FlJnKXK7o1YWFoQsKcqmozKFT3FwNwarsVUQHRjMosmUt7lON1PBUvjn8DRWmCkKMIXXL/e39wmoPHiJwpO8V2NzFlRlBsouPxo4+Z1epxukGC4B7pZTNl9g6NpTyTSnlOCnluOjoaBeG7R6WPEcx2UlDsKtwF7f+eCvVlmZ+5B3A3uK9VJor67RovU7uNvUcP1o990iCi96AvB3wf/d5pExWWFFLrcXWZbqO1ifAT09UiJEcZ5lDoGoK5jwNg86D/7sf9nzZvgNshlmJs7BIC79l/ebyNkII/BLiXZoRtFRDYLFZWJO7RhW4dZF2Iw5SI9QF/8iJht1G/Xr3RhiNnTZO0OqMQEr5mIf7zgbq957tDTT+Bo0DPrZ/GaKAc4QQFinlFx4e02PqJCrrGYLfsn9jTc4aVmav5Oyks9t7SE7ZVbQLgBHRI3xzgNxtEBzd8K5/4Fw47c8qXiB0kOBep9O6rqORXS9GAGpW0OyMAFRPpkvegvfnwdI/QHAMJPrIkLvI8KjhRAdG83Pmz5yfer7L2xni4lwSqDFlZqILCkLfs2eT93YW7qTcVN6l4gMOHIbg0IlDDVJPhV6PMSWl07aacLuKQwjhagngJqC/ECJZCGFEdS9dXn8Fu9xlkpQyCfgMuL0jjACcDIDVjxFklqm7mu+Pft8RQ3JKWlEaMYEx9Ar2USwjd7uaDTS+U5v5D0g8DaSt2+sQNCY+vBVDAOAXCPM/hog+8NHvoPJ4+wyuGXRCxxl9z2BN7hq3ZryqqMwV11AGfomJTu/4V+WswiAMTIqb5NaYTwUSQhIw6ozOK4xTUzF10hRST8r5FrmykpTSAtyBygbaC3xqzza6VQhxqwfH9SnmvFz0ERHogk5erDLKVWXkyuyVVJk7h8hIWlEaw6KG+Wbnpioo3AtxTtpa6A1w6TuQOBWS3LuT66o1BA4cM4JWC62Ce6pMoppSOPRz+wyuBWb1nUW1pdqt1tR+cfFYS0uxVbX8ezBnNJ86uip7FaNjR7fayO1UxKAzkBSe5DRgbExNwZyb2+pn1xF4YghcdupJKb+VUg6QUqZKKZ+wL3tdSvm6k3Wvl1K23pvZR5jzGtYQSCnJLMukf4/+1FhrWJWzqqOGVseJ2hMcLTvqVrWjWxzbpe74HfGBxoTGwg3fuu3WyC6pomewkSBj19AhaEx8RACVJitl1ZbWV04YC/7hcLT5hoPtxbhe4wgzhvFzputGqS5zqIWeQ9JiwZST49QQ5Ffmk16S3uWyheqTGp7aTPO5fiAltZ1QrcwlQyCEeEQI8bAQ4hEg1v76YSHEwz4eX7thyTuGX6+ThqC4ppgKcwUXpl5Iz4CencI9tLtoN4DvZgSNA8VeItuuQ9BVSbBnDjlNIW2MTg9Jp8GRjjcEfjo/ZvSewa9ZvzZR1mp2GxcEaszHjoHZ7DRjaHWOqoHpyoYgJSKFnIqcJl6EzpxC6uqM4CiQYX822187Hl2CxoI0meUqPpAcnszsxNmdwj20q2gXAsHQnj7q1Ji3XfXSd7NqGOBElZlai/Pkr6zirqVD0Jh4dwwBQPI0VZdRmum7QbnIrMRZlJnK2JLvmvy4K7UEpgx7szknM4JVOauIC46rC6p2Reoyh8oa3vkb+/YFg6FTtppwyRA00ig+3tU0i63l5djKyxvUEGSUqS9zYlgiZyed3SncQ2lFaSSHJ/vOt5q7TdUPeMAVb65jzoJVZByvbLDcZpPklHbNqmIHrdYSNCZ5uno+0vHuxinxUwjQB/Bjxo9N3iuoKmDJ/iUNhGwM0dGg17dYS2CuqyFoWExmtppZl7uOaQnTulzaaH2aaz4n/PwwJiW2mkJas3cv5T//4nJzP2/g0xjBqYLD32lolDFkEAbiQ+IZEzOGngE9WXG09S4YleZKn/wDpZTsLNrpO7dQbQUUpnvkFqoxW0nPL+dIUSUX/XctWzJONlorKK/FbJVdsqrYQVSIEaNB57ohiB4MQVGdwj0UaAhkasJUfsn8BZu0UVhVyOK9i7nuu+uYvWQ2j697nJe3vVy3vjAYMMTGtDIjyEQEBCijUY+tBVupslR1ybTR+vQJ64NBZ2g2TmBqpR117r33kX377eT86W4sJe3TtNATQzDL66PoYCx1NQQnC6UyyjJICE3AoDOg1+mZnTibVdmrWnQP5VbkcuZnZ7Jg6wKvjzGvMo/immKGR/kqULwTkB4ZgsziKqSEP88eQGiAgSv/t57vdinj2lXbT9dHCEFCRKDrriGd7mScoB3v+ppjVt9ZFFQXMP+b+cxaMounNj5FmamM20bdxvhe45sWR8XHt1hLYMrMxNinD0LX8PKyKnsVfjo/JvSa4JPz6Cz46fxICnOeOeSfmoopKwtbba2TLaFm3z5q9+8neMoUyn/5hcPz5rWoZOgtWjUEotEcTkrZpAdt43VONepqCOIbxgj6hJ6sh2vNPSSl5MkNT1JuKufDvR9SVF3k1TE6Csl8ZghylTSm09TRVjhcqNxBMwdF8/ltUxgaH8bti7fyv5WHySpWhqCrCdI0Jj4iwPUZASj3UHkuFHe8atWMPjOI8I/AZDVx26jb+PKCL1l2wTJuG3kbo2NGk12RjclqqlvfoUvQHKqGwHl8YHyv8b5pnd7JSAlPaV6tzGbDdPSo0+1OfLkc/PyIf/45kpd8iiGiB1k330LeY4/5NO3UlRnBL0KIO4UQDf6zQgijEOIMIcT7gOsyPp0Qc14e6PV1U1kpJRllGSSGnfRxtuYe+inzJ37L/o0rBl6B2WZm4e6FXh1jWlEaRp2RAT0GeHW/deRuUx1FQ1vvftqYo/a4QFJUMD1D/Fn8h0nMHdaLJ77dy/Pf7wdOZtZ0VVRRmRvtppNnqOcjrrd48BVhxjB+u+K3uou/o68+qAuaTdrqiivBXlSWn4+0Nk0OkDYb5sysJs3mcipyOHzicJfOFqpPakQq2RXZ1FgafidaUiuTVitlX39NyPTpGHr0IGDQIJI+W0LkDTdQ+vEnHLnoYqp37PDJeF0xBHMAK/CRECJXCLFHCHEYOADMR+kYv+eT0bUD0majct06jL17I/Sq/3hRdRHVlmr6hp60fS25hypMFTy14SkGRQ7ivgn3MSdpDh+nf0xJjff8e7uKdjGo5yD89D7q55+7zeO00aNFlUSFGAmzaw0E+Ol5Zf4YbpmeQk5pNdGh/gT4udfb/VSjWYGa5uiZqgxvJ4gTAM1KRSaHJwMNM2D84uLAYsFS2FTPwJKfjzSZmtQQ7CpUM9qxse61JzlVSYlQBvRo2dEGy43JSaDTOU0hrVy3HkthIeHz5tUt0/n7E3vv3+n77rvYTCaOXn1Ni7MxT2nVEEgpa6SU/5VSTgUSUTGCMVLKRCnlH6SU270+qnbkxOefU7NzJz1vPVnsXD9jqD4O99DKnIY/3pe3vUxhdSGPTH4Eg87AzSNupsZSw6I9LhVht4rFZmHP8T2+cwvVlMHxAx4bgsNFlST1bKitrNMJ7j9nMAuuGMV9c7pWh0lnNCtQ0xxCqDTSI6s6RZygOZLCkoCGTdT8EppPIa1rNtfINbSveB8GYejSaaP1cWQONQ4Y64xGjH37Ok0hPbH8S3RhYYTMPL3Je8GTJpKy/EsSnnm6QZq7t3AlRhAghLhbCPEKSqe4UEpZ6vWRdACWkhIKnnuewLFjCb/wpFSCo4agb1jDL7PDPVS/uCytKI2P9n3E/EHz6zJ6UiNSmZ04m4/2fdREyNoTDpUeotpS7buMoWP27uEepo4eLaokOSrY6XsXjk7gkrG9PR3ZKYPbtQSg4gRVRUr/oZMS5BdEr+BeDQ1BMwI1tqoqSpcuBZqmjqaXpJMSkYJRb/TxiDsHiWGJ6IXeaeaQ0Ylama2ykvIffiRszhx0RuefkT40lLBzzvHJeF1xDb2P6hK6CzgHeN4nI+kACl9YgLW8nF4PP9wgrzmjLAODzkBccEPL29g9ZLFZeGzdY0QHRnPn6DsbrHvLiFuoMFfw4d4P2zzOtKI0AEZE+bDjKHgUKK6otVBQXktSM4agu+AQqHE7YAydxj3UHMlhyc0YgpMzgorVazh8/jzKvvqKyOuvb9DOHSC9OJ2BPQa2z4A7AUa9kT6hfTh8wlnzuX6YMjKQppMB+PIff0RWVxN+wbwm67cHrhiCIVLKq6WUbwCXAl0i2lO9cyelS5YQefXVBAxsGIDNLMukd0hvDLqmvXHqu4cW713MvuJ93DfxvgYiFAADIwdyep/TWbRnERWmtsnT7SraRZgxrEEWk1fJ3QbhfSDEfZ2Ho0UqUJzS7Q2Bm0VlABF9ISIRjnZ8YVlLJIcrQ+Coj9EFB6MPD8ecl4ulpITc++4n66abEEYjiR9+QOx99za4sTpefZzC6kIGRnYfQwDKM+C8liAFLJY68R6AE8u/wi8hgcDR3m3v4iquGIK6JiT2jqKnPNJq5dijj2GIiiLqzjuavJ9RntEkPuDA4R76aO9HvLL9FWb0nsHsvrOdrnvriFspN5XzcfrHbRrvrqJdDI8a7rtqzDZUFNfPGOrOtCpQ0xzJ05UhsLWqzdRhJIcnU2WpoqDqpJqsISGeyjVrOXze+Zz4+mt63noLyV8sI2hs02BweomqTO6OhiCgfFm0AAAgAElEQVSrPKtB6i2A0aFWdkjNFswFBVSuW0fYvPOb1F60F64cdaQQosz+KAdGOF4LIdruAO8ASj75hJo9e4i57170IQ3v5G3SRlZZVpP4gAOHe2hrwVYAHpj4QLMX6KFRQzkt4TQW7l7ocZ+iKnMVB0sP+i4+UF2qctk9cAsBHLHXEDQOFndHWhWocUbyDKg5cTJO0wlxmjkUH485MxO/uDiSP1tCzN13o/P3d7r9/mKVQuxz19D2xfC/WVCe79vjuEhqeCpWaa1LPnHgn5ICQtSJ1JR9/Q3YbA2yhdobV7KG9FLKMPsjVEppqPfaNbHTToTl+HEKF7xI0KRJTgMvBVUF1FhrSAx1PiMAmJs8F4A/jvoj8SEtyzbeMuIWSmpLWLJ/iUfj3Vu8F5u0+S5jKM+el+xhxtCR45XEhQcQaOza6aGu4JJATWOS7Z7WTtB3qDnqDEG9OEH07bcT9+STJH38EQGDWs4K21eyj5igGHoE9PDdINe/Bl/cBjmb4benfXccN6ivVlYfXWAgfgkJdSI1J5YvJ2DECPyTk9t9jHVj6rAjdxAFzz6HrbqaXg8/5PRO3lE409yMAFQu9NJ5S7l2yLWtHm9UzCgmxk3k3bR3mxSXuIIjUNxZW08fcZI62l1xWaCmPqG9IGpApw4YRwdGE+wX3MAQBAwZQsTFFyEMrWtMpBen+06kXkr47d9KT3vQeTDmOtjyPhR1vCRkYlgiOqFrVq2s9tAhatL3U7tvX4fOBqCbGYKqzZs58cUX9Lz+ejU9c4JDlay5GIGDAT0GuOyzv2XELRyvOc7SA0vdGzAqPpAQkkDPwKbar14hd5sKWAZFerT50aJKkqM1QwAnBWpOVLvW27+O5OmQsRasbm7XTgghmmQOuUqttZYjJ474xi0kJfzwEPzyBIz4HVz2PpzxkJIF/flx7x/PTQIMAfQO6c3B0qZGydgvFdORI5xYtgwMBsLO9U1aqKt0G0MgLRaOPfY4hvg4om5rXikzsywTo87oVU3g8b3GMzZ2LG/vetstfVhQFZk+mw2A0iDwMFBcWmWipMpMsjYjANwUqKlP0jQwV0LOVh+Myjs4Mofc5VDpIazS6v1Asc0KX98Na1+G8X9QEqB6g8p8m3In7PkSsjd795gekBLRTM+h1H5Ik4mSTz4hZNo0DD186DZzgW5jCNDrifrj7cQ9+mgDXeLGZJRl0Ce0T7Ml955y1+i7KKxWLX5dpai6iNzKXN/FB6qKlUBKG9xCoGUMOTiZQuqmC9ChAd0J5CubIzk8mfyqfCrNla2vXA+HloFXZwRWM3z+B9jyHky7B855VnV0dTD5jxAcDT883OFV26nhqWSUZTRRgHOolXVk7UB9uo0hEEIQNmcOIdOnt7heZllmi/EBTxkTO4ZpCdN4J+0dl6uNHdKUvgsU27uDeNpjyJ462lxVcXfDo1oCUKL2scM7dZzAETBu3DunNdJL0gk0BHqvBqa2Aj6aD2lLYfZjMOth1a6jPv6hMONeyFgDB37wznE9JDUiFYu0kFWW1WC5MUUZAl1ICCEzZ3bE0BrQbQyBK9ikjazyrFbjA55y15i7KDOV8V7aey6tv6toF3qh912g7dAv6jlupEebHymqQiegbxdvMe0qPYPdFKipT/J0yNwAZvcTCtoDZ5lDrrCveB/9e/RHr/NCVll5Prx3Dhz6Gc5/EU67u/l1x14PkSnw4yMdWqPh6OTaOE6gDwnGf+BAwi+6qNm02/ZEMwT1OFZ5DJPN5JMZAcCgyEHMTZrLB3s/cEmvIK0ojX4R/XzTvz1vJ6z/Lwy/DAI9808eKaqkd48gjAbtawSq0Z5bAjX1SZ4G1lrY8VGHuzOc0Se0D3qhd8sQSCnZX7yfQT28cCNTdADenq2e53+kLvQtofdTgeOCPbDzk7Yf30NSw1Mx6Ax12X/1SV7yKbH33dsBo2qK9guuR13X0RZqCNrKH0f/EZPVxJs732xxPSklu4p8FCi2mOCL2yEwEub+2+PdHC2q1OIDjXBboMZB0jSIGaICoG+f2encREa9kd6hvd0yBLmVuZSby9seKM5crz4TczVc/w0MONu17YZeBPFj4OcnOmymFWAIYFjPYWwp2NLkPWE01rW+72g0Q1APV2oI2kpiWCIX9b+IJfuXkF2e3fxYyjMpM5UxItoHjeZWPQ/5u+D8BR6njUopOVJUSXJPzS1Un/hwD2cE/iFwy0o4/yUoy4X3z4f353WKzBcH7qaQ1gWK22II9nypPofASLjxB0gY4/q2QsCZj0FZNmxs+cbLl4yNHcueoj0edxdoD1qvBulGZJRnEKAPICYoxqfHuXXErXx16Cte2/EaT5z2RJP3i6qLeHLDkwCMjPbMf98sudth1XMq73rQuR7vpqjCREWtRQsUNyI+IpCC8lpMFpv7LjO9H4y9DkZcAVvehZXPwVuzYOA5kHqGe/uK6Ov6nbOLJIcnsyZ3DRabxWlDxsakF6cjEPSP6O/ZATe/C1//GXqPh/kfq6C6uyRPh36z1c3PmGs8doO2hTGxY3g77W12Fe1iYtzEdj++K2iGoB6ZZZn0CfN+6mhjYoNjmT9oPgv3LOT3w37fQKxjbc5aHlj9AOWmch6a9JB3hTwstcolFBQFc9tWhq81m3OOQ6Amv6zGc51mvwCYdBuMvgY2vAZrXob0b93fzx1bIKqfZ2NwQnJ4MmabmdyKXJdmzekl6SSGJXoW46oth+8fUhfy+R+DsQ0zz9mPwRvTlIvo3Oc834+HjIoZhUCwNX+rZghOBTLKMugX4b0fTkvcOOxGPtv/GS9ve5kFMxdgtpp5adtLvLf7PfpF9OPNs970vj7xb/+Ggt0w/5M23xk5ms1pM4KG1Beo8dgQOPAPgel/g8l3qgujq9SUwn8nw6a32mzw61M/c8gVQ7CveB9Dew717GA7PgZTOcx6pG1GAKDXMJhwM2x4Q822+oxv2/7cJMwYxsDIgU7jBJ0FLUZgx2KzkF2R7dP4QH0iAiK4buh1/JT5E98c/oarv7ua93a/x+UDLmfxuYu9bwRytsLqF2DUVTBwTpt3d+R4JX560eVF6d3FI4Ga1vALUBWzrj6i+sPQC1U3TpN7BWAt4Uy2sjnKTeXkVOR4lvosJWz8n6pv6e0ljeMz/gGhcSoY3wGtPMbEjGFn4c4mhWWdBc0Q2MmrzMNis/ishsAZ1wy5hsiASO5bdR/Z5dksOH0BD01+iECDly+ullrVmTEkFs5+0iu7PFJYSZ/IIAx67StUH4+LyrzN+D9A7QnY+anXdhkREEFkQGSDdtTNsb/E3nrak0Dx0VVQlK7u4r2Ff6iqQM5Pg3Wvem+/LjImdgzVlmr2Hu+csqTar9hOXcZQaPvMCACC/YJ5YOIDnJV4FkvnLWVW4izvH+RYGiy9CQr3wbyXIDDCK7s9eryy26uSOSPAT0/PYA8EarxNnwnQa7i6s/ZiXYKrPYfa1Fpi45sqS2joxe5v2xKDz1MdSn99WrVWaUfGxqqZzdb8ztlPSjMEdupqCNpxRgBK+vL505/3apM7zNXKLfDWmfD6VNi/QpXc9z/TK7u32SRHj2vtp5vDI4EabyOEmhUU7IbMdV7brcuGoCSdCP8I9zPwTmTDvm9Vho9fgIejbIG5z4BOD9/c066Fe1GBUSSGJbIlv3PGCXxqCIQQc4QQ6UKIg0KI+5y8f5UQYqf9sVYI4eVcSdfJLM8kyBBEVGBURw2h7Rw/BN/dC88PVK6g6hLlCrpnH8x8wGuHOVZWQ43ZpmUMNYPH1cXeZvhlEBCuZgVeIjksmdLaUkpqSlpcL704nYGRA92XV938LkgbjLuxDaNsgfDeKl5w8EfY/blvjtEMY2LGsLVgKzZpa9fjuoLPDIEQQg+8CswFhgDzhRBDGq12BJghpRwB/BPosKqPjLIM+ob19Z0usK+pLoX/zYTN70C/M1UF5h2bVCdGD4vGmkMTrG8ZjwRqfIExSKWg7l0O5ce8sktXeg5ZbBYOlBxw3y1kqYWt78OAOdDDhzPzCTcradbv7lO/m3ZiTOwYykxlTgXtOxpfzggmAAellIellCbgY+CC+itIKddKKR23FuuB3j4cT4tklmW2a3zA62x8U2nf3vgDXPo2JJ3WtCujlzistZ9ukfiIAKo8EajxBeN+DzaLatnsBVwxBBllGZhsJvczhvYsh8pCmPCHtgyxdXR61bSuqgh+fNS9bfd+BZ/9Hsry3D7s2JjOGyfwpSFIAOr3Xs22L2uOG4HvmntTCHGzEGKzEGJzYWGhl4aoMNvM5FTktHt8wGvUlqsGcgPP8Vhkxh2OFlUS4KejV5gPfLhdAI8FanxBz1RVWbv5Xa+kTcYFx+Gv92/REOwr3gfgfgr0xjchMhVS2qEtc/womHibquDOcDGGkr8Hlv5BtcB+Yxoc/s2tQ/YO7U1MYEynjBP40hA4ux11OlcWQsxEGYJmW/FJKd+UUo6TUo6Ljo720hAVuRW5WKW13WoIvM6mt1U8YPpf2+VwjkCxTneKutF8jMcCNb5i/B+g4hjs+7rNu9Lr9CSGJbaYQppeko6fzo+UcOdysE7J3Q7ZG9VsQNdOOSwzH1AyrZ9eA8WtBMBrK2DJdSoN9drlKqtp0YWw8lmwuebzF0IwJnYMWwq2tOg2/CnzJ17d/ioFVQXunE2b8OUnng3UV6PoDeQ2XkkIMQJ4C7hASnnch+Nplo7KGPIKpipY9wqkzoIELxXftMJhTbC+RTpNLYGD/meq3kNeChq3ljmUXpxOakQqfno/13e66X/gFwQj53thhC7iHwJXfaZcZx9cDBXNeBqkVD2Pjh9UbteUGfCHn2HYJfDzv2Dx5UrtzwXGxI6hoKqAnIocp+8XVRfx4OoHeX3H65y99GweXvOwU6lLb+NLQ7AJ6C+ESBZCGIHfAcvrryCE6At8Dlwjpdzvw7G0SEfUEHiNre8rv+r0v7XL4SxWG1nFVZpgfQu0SaDGF+j0MP4mpdiVv7vNu0sOTyanIodaa63T99OL090LFFcVw67PVPsHL9W5uEz0ALjyU+XzX3y5uvNvzNb3YdencPoDqvcRKCNy8f/g3P/Akd/gjekudYqtqycocB4neGXbK9Raanl99utc2v9SvjvyHRd8eQF3/nwn2wq2eXyareEzQyCltAB3ACuAvcCnUsrdQohbhRAO9fiHgZ7Af4UQ24UQ7d5zV0rJtoJthPiFEBng3ewan2OugTUvQuJpkDi5XQ6ZU1qN2So1wfoW0OkE8eEBZHcWQwAqe8gQoPoPtZHksGRs0lZ3A1Wfouoijtccd6+iePuHYKnxfZC4OfpMgEvfUdKtS65rGEvJ2wnf/l3NuKfd03A7IWD8jfD7Fer1u3NVCncL9IvoR5gxzGmcYF/xPj4/8DnzB89nasJUHpz0ICsuXcFtI29je8F2rv3uWq7+9mqXRK3cxafOOCnlt1LKAVLKVCnlE/Zlr0spX7e/vklK2UNKOcr+GOfL8Tjj1e2v8n3G91w28LJTL3V0+4dQngczvD8bqKy1OF3uEKzXZgQt0ymKyuoTFKlcGTs+aXPKZEuZQw6dbZczhqwWZZwSp0Kshw3qvMGgc+C8F1R9wfK7lDuopkwZhqCecPGbzccuEsbADd8pA7LrsxYPoxM6RseMbpI5JKXk35v+Tbh/OLeMuKVueWRAJLePup3vL/2eByY+QIhfCD38vd9Ku1tXFr+16y3e2PkGl/S/hLvHtKB/2hmxmmH1AtWrPXmGV3f9S3oBwx5dwSNfplFjbqj36qgh0GIELZMQEcjRokqOVzh3n3QIE28Fc5VqsdAGHLE0hyGQUrIhbwP3/HoPd/9yNyF+Ia7PCLa8q9o9TL6jTWPyCmOvV+6fHYvhp8dg+Z1QkgGXvQvBrRSahvdWxsyFIrUxsWM4Wna0wZ39z5k/s+nYJu4YdQfh/uFNtgk0BDJ/0HxeP/N17+g/N6LbGoIP9nzAi1tf5NyUc3lo0kM+1yDwOjs/gROZMP3vXq8XeGf1EQIMet5fl8EFr6xh37GyuveOFFUS4m8gKsTo1WN2Nc4e2ouKWgtnPP8bH2/MxGbrBDrEcSNUXcHGN1SWjocE+QURFxxHWlEaC3cvZN4X87jp+5vYcGwDVw2+ik/O+4QwY1jrO6ouhV+fUjKdA+d6PB6vMuPvyiCsfgH2fAGzH4G+k1zbdthFqqdX/p4WV3PECRw+f5PVxHObn6NfRD8uGXBJW0bvMafY1c87LNm/hGc2PcOZiWfyr6n/8omF9Sk2q1Jc6jXCa/2DHBwtqmTVgSJuOz2V924Yz/FKE/NeWcO7a44oecrjVSRHBZ96brR2ZvaQWL69axoDe4Vy3+e7uPyNdaQfc0NTwFfMeli5Or7+s/oeeUhyeDK/Zv/Ks5ufJdw/nCdPe5IfL/2Rv47/q+tp2KueU4His5/wWfGj2wgB5zwPo65Wj8l3ur7t4AtA6FqdFQyJHEKAPqAuTvDB3g/Irsjmb+P/5pLymy/odsI0Xx36in+u+yfTEqbxzLRnOuyDbxNpn0PxYbh8kdd/QIs3ZmLQCX43vg8xYQH8393T+NuSHTz21R5+21/I/mPljE8+xYLqHUT/2FA+uXkSS7Zk89S3ezn3pVXcOC2ZP83qT5Cxg753gRFw9lPw+U3KLTP+Jo92c82Qa0gJT+HCfhd61mq6+LASihl1FcR1WIsx5+gNcKEHrapDolVWUdrnMPPBZn+bfno/RkSPYGv+Voqqi3hz55uc3vt0psRPaePAPadbzQhWHF3BP9b8gwlxE3hh5gut5zlLqTp5diZsNnUnFT1YtdT1IjVmK0s2Z3HW0Fhi7FXDUSH+vHP9eB6bN5S1h45zrKxGUyVzAyEEl4/rw0/3nM5FoxN447fDnL1gJYXlHRg7GH6piiv9+DiU53u0i9MSTuPeCfd6Lkz/46OgM6gGcF2JoRdB8SE4trPF1cbGjiW9JJ2nNz5NrbWWv45vn2LQ5ug2hsBsNfPKtlcYFT2Kl2a+hL/ev/WNdnwM/06BvB2+H6CrpH+j/JDT7vF6Bea3u/IoqTJz9cSGhXVCCK6bksRXd5zGnKG9mDvMiy2zuwmRwUaevWwki2+aSF5pDS/9dKDjBiOEyn+3VMP3D7b/8TPWwp4vYerdEBbX/sf3JYPnKQOX1rJ7aEzsGGzSxoqjK7hq0FUdXszabQyBn96P/531P16d9arrYtrbPlBZFktvUhW8HY2UsOo/0CNZ3Xl4mQ/WZ5ASFczk1J5O3x/YK5TXrxnL4DgXAoEaTpnSL4r5E/qyeGMmhwudFC+1F1H94LQ/w64lcOiX9juuzQYrHoDQeJjihv/9VCEoElJOh93LWtQ7GBE1AoMw0MO/BzeP9KISm4d0G0MA0Cu4FyHGENdWLstVlZips6Bof8fcOTXmyG+QuxWm/kn5Mb3IntwytmaWcuXEU7gV9ynCXbP6E2DQ8eyK9I4dyGl/UTcV39yjihPbg12fQu42lY3TVlH6zsrQi6E0Q/1WmyHIL4i7xtzFv077l2sZVj6mWxkCt9i9DJAw998qx3nzO0o5qSNZ9R8I6QWjrvT6rj/ckIG/QcelYzusE3i3ITrUn5unp/Jd2jG2ZLQs8OJT/ALg3OeVT3vNiyeX22xQmqmKqza84ZW2FICaVf/4mBKlH365d/bZGRl0LuiNrbqHbhh2A9N7T2+nQbWMZgiaI22pymaI6qdS7noNh+V3tCzwUbAPFl4I61/z/niyt6gZweQ/gsGF+IYbVNRa+GJbDuePjCciSKsPaA9umpZMVIg/T3+3t2MFbPrNUnewq55XffZfnwZPJcCC4fDBJfDd3+GL271zrHWvQHmuUs1rrw6jHUFghPIk7P7C5c6kHU0X/m+0geIjkLNFleSDuvBe8o66o1l2a9N/rpSqTP7NGXB0NfzffbDiQe9+CVb/BwIiYNwN3tunnWXbcqg0Wbl60inYffUUJdjfwN2z+7PpaAk/7m2/dsNOmfOUqpzN2gQhMaqg6rwFcP23cMZDqgdPWxMmyvNVkdaQCyCx49Ik242hF0FZNmRv6uiRuIRmCJyRtlQ9D7345LLoATDnSTj8ixKBcVB5HD6+UvlZE6fC3buUFN66V2DZLWAxtX08BftUL/kJN6t+6F5ESsmH6zMYlhDGyN5NS9s1fMcV4/uQEhXM09/txWLtwDvH0F7wlz3w511w9VJlGMbdAElTVVM1QwBseb9tx1j3imosN+sR74y5szNwLuj9210X2VM0Q+CMtM+hzySI6NNw+dgbYOC5qg9J3k6VbfHaFOVLPfsp1ds8LE7FFc54SAXGPrrCeWtbd1izQPVqn3hr6+u6yZaMEvYdK+eqiYlakLid8dPr+PucQRwqrGTJluyOHo5zAnuou/hdS8BU6dk+qktUjG3oxUoxzctYbZIf9uRzzdsbmPDEjyzf0UT2pP0JCFNV/7u/aFMFd3uhGYLGFOyFgt0n3UL1EQLmvXxSnWjRhRAQrkQqJt9+0u8phFILm/eKkrN7/7zmRS9aoyQDdn6qpuvBztM628IH6zMI9Tdwwah4r+9bo3XOHhrL2MQe/OeH/VSZnHd87XDGXAe1Zeqi5gmb3gJTBZzm3caORRW1vPrLQab/+xf+sHAzB/IriA71566PtvHktx08ywIYdrFShstc37HjcAHNEDQmbanqFzLkAufvB/eEi15Xd/njboSbf1WBZGeMuQZ+96EyLu+c1bocnjPWvqzG44PujMWVJr7ddYyLxyR0XMuDbo4QgvvnDqKwvJa3V3nw/WgPEqdAz/6wdaH725qqVPJE/7Oa/524gc0mWX/4OH/6eBuTn/qJZ1ekk9gziNevHsPqe2ey7PapXDMpkTdXHub6dzdRUukF16ynDJgDhsBTwj2k/frrI6UyBEnTIDS2+fVSZ8L92WBwIcNm4Fylcbr4cvjfTLjwNdc7LVYUwLZFMPIKCE9wbRsXyS+r4cFluzBZbVylBYk7lHFJkZw1JJY3Vh5mbGIPRvWN6FyGWQgYcy388JCKV8W4qDUA6vtbdVzVLLSBjOOVLN2aw+dbs8kuqSbU38BVExO5elJf+sU0jJv988JhDE8I5x9fpHH+K6t545qxDI3vgPiXMRgGnK2qqOc84/XaH28iOjR1zUPGjRsnN2/2gZhZ7jZ483Tl/hlzrXf3ffwQLLle9SCZfIcKmrVmSH58TGVa3LEJovp7ZRgWq43312Xwwg/7MVlt3HPmAG6Z4X2/rYZ7HCyo4IJXVlNpsqLXCQbEhjK6bwSj+kQwpm8E0SEBbu0v0KjHaPDihL+yCJ4fpBIW5jzp2jZWM7w0GsIS4MYVbh+yotbCNztz+WxLNpuOliAETE2N4pKxCZw9tFerxnJ7Vim3LtpCabWJZy4ZwQWjvHsz5RJ7voRPr4WrlkL/2e1//EYIIbY4EwDTDEF9VjyoCmj+ul+Vinsbcw18/w8l1J0wTsnj9WjmbrwkA14/Tc0+LvdgSu6ELRnFPLgsjX3HypkxIJrHLxhKoiYw02korTKxLbOUbZklbMsqZXtWKeU1nsUNIoL8ePOacUzwZqfYT6+DIyvhnn2u1bJs/wi+uFVpAg84261DHSqs4Nq3N5JTWk1qdDCXjO3NRaMTiAsPdGs/heW1/HHxVjYeKeai0Qk8eO5gokK8W4fTIuZqeHmcqqK+ZZUq4utANEPQGjYbLBim/JhXfuLdfTdm9xdK/UgIuOC/MPg85UvNWAuHfoJDP6vGcjoD3PQTxI/y+FA2myS7pJpXfznIJ5uziAsP4JHzh3D20F5allAnx2aTHC6qZHtWKWXV5tY3sCNRleI5JdW8euUYZg9pwc3pDod+hkUXwSVvqw6mLWGzwX8nqe/wbWvcape+M7uU69/dhE7AK1eOYWJyZJu+q2arjZd/OsBrvx0i2N/AA3MHc9m43u33/T/4E3xwMUy5C876Z/scsxk0Q9AaGevg3Tlw8Vsw4jLv7tsZxYdhyQ2qWCdhLBxLA2utyj1OnKIqPgfMccsllFNazbbMEg4VVHKosIJDhRUcLqyk2mzFoBPceFoyd83qT7B/5/VVaniH4koTN7y7kbTcMp6+eDiXjevT+katYbPBS6PULPa6r1ped983qr7Gzd/T6gNF3LJoMz2CjSy6caJXW54fyC/ngWW72HS0hAnJkTx50bAm8QWfsfwuFS/5/ffQZ3z7HNMJmiFojW/ugW0fwt8Ogr+LjenaiqUWfnpcVSMnToV+Z0DfKW4345JS8sH6DP759V5MVhtCKM3c1OgQ9YgJZnJKT1Ki2+m8NDoFFbUWbl20hdUHi3jgnEHcPN0LsaCVz8LP/4K7tkFkivN1pIS3ZkNVEdyxxeUg6be78rj74+0kRwWz8MYJxIZ5341is0mWbMniyW/3UWWycNuMVG6f2Y8Av9ZVCmstVpZuyWFP3gkmpfRkWr9owoNa0TRxUFOmao4MAXDrKvBzz8XlLTRD0BJWCzw/EJJOg8vbWEHZzpTVmLl/6S6+2ZXHzIHR3HPWQFKjQwg0nmLymxo+odZi5S+f7OCbXXncMiOF++YMaptLpCwPXhgKU++C2Y86X+fISnj/fKV5MP5Gl3b74YYM/vFFGmP79uDt68a7foH1kKKKWv719R6+2J5LRJAfl47pzZUT+zq9Wao2Wfl4UyZv/HaYY2U1+Bt01Fps6HWC0X0iOH1gNKcPjGFIXBg6XQufrcO11pqLSEqwWaA14SwP0AxBSzj+QZcvgiHzvLdfH5OWc4I/Lt5Kdkk1fzt7IDdPS2n5i6jRLbHaJA9/mcaHGzK5fFxvHr9gmEt3wM3y0XzI3qzaUji7WC28UHUsvXtXq8FRq03y8s8HWPDjAc4YFMOrV45p15uYTUeLeW/tUVakHcNik0zt15OrJiZy5pBYai02Plyfwf9WHaaowsSEpEjunNWPySk92Z5Vyq/phfy6v4C0nDIA4sMDeHZBVM4AAA9iSURBVO3qsYzsE9H8Ab/6k6rHaM5FdOhn+O4+OH4QeiRB9EDlHo4aaH89QFUte4hmCBpjs0HWetj1mWo5bbPAXw90eFTfFaSUfLAhk39+tYfIYCMvXzma8UmajrBG80gpeeHHA7z00wESIgK5b+4gzhsR59nsIP3/VOuUKz5QcqnVJVByRBVMFu5T7qPZjyrhmxbYmV3KP75IY2f2CS4encAzl47AT98xNa4F5TUs2ZzN4g2Z5JRWExPqj8lqo7TKzLT+Udwxsx8TU5xX9heU17BqfxELftpPaZWZD2+ayIjezRiD5lxExUdURuG+r+3CUxeqOGLhfmUUbPWSBe49qlp/eIBmCEBNufK2q6KxtGWqO6AhEAbOgUm3Q58JDVYvrjSx6kAhswfHdpoA64lqMw8u28XXO/OYMSCaF64YRWSw1jpawzXWHizin9/sZW9eGaP7RvDQeUMY09fNi4rVojLszNXqN1V7ouH7cSNVMDnAeRHXiWozz3+fzqL1GUSF+PPQeUM431Oj5GWsNslv+wv4aGMWBp3g5ukpjHbx88kuqWL+/9a3bgwO/aLa00y5E06/X9UKrXlJZVhN/2vTVvNWC5QchaJ09Tz5jx6fn2YIpFTFYnnb1QfebzYMu1RV+ToJDueUVnPNWxs4XFRJRJAf101O4ropSR160V17qIi/frqD/PJa/nLmAG6bkaq5gjTcxmqTLN2azbMr0iksr+X8kfHcO2cgvXuoJAWTxUZlrYWKWguVJguJkcFN3TU7P1WPHonqDjYyWT33SGo22UFKyfIdufzz670UV9Zy7eQk/nLWAMICfBsPaE8cxuBElZkPWjIGX/1JdXQN7QXleUqo58zHIMy3Pb80QwCweoESjRg8r8WCscOFFVz91gbKayw8dN4Qftybz/d78gn003PF+D78YXoKCRFqSldjtrIjq5TNGSVsOlrM9qxSEiICmTUohpmDYhjZO6LNF+sas5VnV6Tz9uojpEQF88IVo1r2Q2pouEBlrYU3Vh7mzZWHsNkgJMBARY0FU6NmbX0jg3jz2rEM6uW5bzrzeBX3L9vJmoPHGdk7nH9dOJzhXbTteXZJFb97cz1l1WY+vGmS8/OsLYc3Zqg2FHP/DYmT22VsmiFwkd25J7junY1ICe//fgLDEtQ/8UB+OW+sPMwX23IAmDU4hsLyWnblnMBsVZ/hgNgQRvWJ4EhRJVsySrBJiAoxMmNADLMGxzC1XxThge7d/ezOPcGfP9nO/vwKrp2cyP1zB2sZQRpeJe9ENe+sPkK12Uqwv4FQfwPB/gZC7O7QZ1ekU1Fr4bnLRnLO8Di39m2zSRatz+Dp7/Zh0An+PmcgV05MRN/FZ7IuGQObVTWUbEeXmGYIXGDz0WJueG8Tof4GFt00kVQnqWS5pdW8teoIy3fkkNgzmHFJPZiQFMnYxB4NZB5LKk2sPFDIT3sL+G1/ISeqzQgBA2NDGZ8UybikHoxPiiQ+omE+sc0mKasxU1RhYsXuYyz4cT8RQUaevXQEpw+M8fo5a2i0Rn5ZDbd9sIWtmaXcfnoq95w10KULeebxKv6+dAfrDxczfUA0T188vMn3vStT3xg8dsFQzhsR32HBcAeaIWiFX9MLuPWDLcSHB7Lopol1rh9vYLHa2JZVyvpDx9mUUcLWjBIqalUPmYSIQBJ7BlFSZeZ4RS3FlSYstpP/k7nDevHkRcPpoQWENTqQWouVR5fv5qONWcwYEM1LvxvdbK6/zSb5YIOaBeiF4B/nDebycX06RTC4vckuqeKm9zez71g58eEB/P60ZH43oW/dbKu90QyBnfIaM/llteSX1dgfteSUVvHJpiz6x4Sy8MYJPm9KZbVJ9uaVsfloMZsySsgrrSYy2J+oECM9Q4z0DPanZ4iRPpFBjO4T0S1/QBqdkw83ZPDo8t3ERwTyyPlD0AlBjdlKlclKtdlKtcnKj3vzWX+4mGn9o3jmkhHdahbgDJtN8uv+At747TAbjhQTGqBaaN8wNYnoEH+ySqrYd6yc9GPlpOer5xPVZoKNeoL9DQQbDQT569WzUc8/L/S8DqRDDIEQYg7wIqAH3pJSPt3ofWF//xygCrheSrm1tf16agjOfmEl6fnlTZaHBhiYmBzJ85ePctuHr6HR3dh8tJjbPtxKYXmt0/dDAww8eM5grhjfPWcBLbEjq5Q3Vx7mu7Q89DqBn15HlemklGXfyCAGxIYSHWqkstZKlclCZa2VSpOFylr1evW9MzF46GJqd0MghNAD+4EzgWxgEzBfSrmn3jrnAHeiDMFE4EUp5cTW9u2pIXh3zRFMFhu9wgOICQ2wP/t3mhoBDY1ThZJKE7tzywg06gj0MxBo1BNk1BPgpyfYqPf4QtVdyDxexQcbMjBZbAzqFcrAXqEMiA31+bWoOUPgy6NOAA5KKQ/bB/AxcAGwp946FwALpbJG64UQEUKIOCllni8GdMPUZF/sVkOj29Ej2Mhp/aM6ehinLH17BvHAOYM7ehh1+NJsJwBZ9f7Oti9zdx0AhBA3CyE2CyE2FxZ6KASvoaGhodEEXxoCZ87Bxn4oV9ZRC6V8U0o5Tko5Ljo6us2D09DQ0NBQ+NIQZAP11TB6A7kerKOhoaGh4UN8aQg2Af2FEMlCCCPwO2B5o3WWA9cKxSTghK/iAxoaGhoazvFZsFhKaRFC3AGsQKWPviOl3C2EuNX+/uvAt6iMoYOo9NEbfDUeDQ0NDQ3n+DRXSUr5LepiX3/Z6/VeS8DznqoaGhoaGm1GS/bV0NDQ6OZohkBDQ0Ojm3NK9hoSQhQCGa2sFgUUtcNwOhvaeXcvtPPuXrT1vBOllE3y709JQ+AKQojNzkqpuzraeXcvtPPuXvjqvDXXkIaGhkY3RzMEGhoaGt2crmwI3uzoAXQQ2nl3L7Tz7l745Ly7bIxAQ0NDQ8M1uvKMQENDQ0PDBTRDoKGhodHN6XKGQAgxRwiRLoQ4KIS4r6PH40uEEO8IIQqEEGn1lkUKIX4QQhywP/foyDF6GyFEHyHEL0KIvUKI3UKIP9mXd/XzDhBCbBRC7LCf92P25V36vB0IIfRCiG1CiK/tf3eX8z4qhNglhNguhNhsX+b1c+9ShsAuj/kqMBcYAswXQgzp2FH5lPeAOY2W3Qf8JKXsD/xk/7srYQHukVIOBiYBf7T/j7v6edcCZ0gpRwKjgDn2jr1d/bwd/AnYW+/v7nLeADOllKPq1Q94/dy7lCGgnjymlNIEOOQxuyRSypVAcaPFFwDv21+/D1zYroPyMVLKPCnlVvvrctTFIYGuf95SSllh/9PP/pB08fMGEEL0Bs4F3qq3uMufdwt4/dy7miFwWfqyCxPr0HSwP8d08Hh8hhAiCRgNbKAbnLfdPbIdKAB+kFJ2i/MGFgB/B2z1lnWH8wZl7L8XQmwRQtxsX+b1c/dpG+oOwGXpS41TGyFECLAUuFtKWSaEs39910JKaQVGCSEigGVCiGEdPSZfI4Q4DyiQUm4RQpze0ePpAKZKKXOFEDHAD0KIfb44SFebEWjSl5AvhIgDsD8XdPB4vI4Qwg9lBD6UUn5uX9zlz9uBlLIU+BUVH+rq5z0VmCeEOIpy9Z4hhPiArn/eAEgpc+3PBcAylPvb6+fe1QyBK/KYXZ3lwHX219cBX3bgWLyOULf+bwN7pZT/qfdWVz/vaPtMACFEIDAb2EcXP28p5f1Syt5SyiTU7/lnKeXVdPHzBhBCBAshQh2vgbOANHxw7l2uslgIcQ7Kp+iQx3yig4fkM4QQHwGno1rT5gOPAF8AnwJ9gUzgMill44DyKYsQ4jRgFbCLkz7jB1Bxgq583iNQgUE96gbuUynl40KInnTh866P3TX0Vynled3hvIUQKahZACg3/mIp5RO+OPcuZwg0NDQ0NNyjq7mGNDQ0NDTcRDMEGhoaGt0czRBoaGhodHM0Q6ChoaHRzdEMgYaGhkY3RzMEGqcMQoi77F1HPxRCzHOnu6wQIkkIcWUz753u6GrZGRBCjLKnQWtotAtdrcWERtfmdmCulPKI/e8mxYJCCIOU0uJk2yTgSmCx74bnNUYB44BvO3ogGt0DzRBonBIIIV4HUoDlQoh3gBJgnJTyDiHEe6gurKOBrUKI5cCL9k0lMB14Ghhsb9r2vpTyhWaOMx6lC3uJlPJwveVDgXcBI2omfQlwLVAkpXzRvs4TqMK+ncBj9tejgM9RBXB/AgKBC6WUh+zjrgGGArHAX4DvgceBQHvx3FPAD8A79vOvAm6WUu4UQjwKJANxwAD79pNQbdhzgPOllGYhxNPAPFQL7++llH91/ZPX6BZIKbWH9jglHsBRIMr++nrgFfvr94CvAb39769QzboAQlA3PKcDXzez39Pt208BtgB9nazzMnCV/bURdUFPArbal+mAQ0BP+/5KURdof9RF+TH7en8CFtQb9//Zt+2P6pUVUP/c6h37EfvrM4Dt9tePAqtRLalHoozEXPt7y1DtiSOBdE4Wj0Z09P9Re3S+hxYj0OgqLJGqOyfAGuA/Qoi7UBc+Z66ixgxGzQTOl1L+f3t3z1pFEIVx/P/Exiak8ANYpBELIYVWFmoviCAWCqJYWMZK8QP40loHDAgWIoiNIYKKNgqKogEhRILYqEUgRCTG6H0szl5NwhqSiCC5z6/au3N3dqo9M3PgzPuW9ifABUnngO22522/A2YkDVF1YF7anmn+/8x1dsICFSDuNfcnqADSddN2x/YUMA3saHn3XuA6gO0HwDZJA03bmO3Fpt8tVGBZ+p45atUxIukwFSwilkkgiM3iS/fC9mXgNDVrfyqp7eO60gfqgznU1mj7BrW9Mg+MSzrQNI1QM/iT1PZN18KS686S3x2Wb8murPHSVvNltfLqC834OsCi7e79DtDNl+yhqrUe4negiPglgSA2HUmDtidsXwGeU7Psz0D/Ko/NUqdgXWyre98UAJu2fZVKUu9qmm5T5aB3A+MbGO4RSX2SBqkcwGTLWB8Dx5px7KPyEnNr6bw5t2HA9l1gmMpZRCyTZHFsRsOS9gM/gDfAGDVD/i7pFTDqlmSx7U+SDgJjkk65TgDrOgocl7QIfKQSutj+JukhMLtka2o9JoFHVLL4jO2vTX/nm8T2JSoXcE3Sa2pr58SfOmvRD9yRtJVaWZzdwBhjk0v10Yi/IKkPeEGVAp5a57OjVAL71r8YW8RaZWsoYoMk7QTeAvfXGwQi/idZEURE9LisCCIielwCQUREj0sgiIjocQkEERE9LoEgIqLH/QTIVHrsB7dwlwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def generate_prob_positives(prevalence, nigam_df):\n",
    "    num_symptoms, probs_positive = [], []\n",
    "    all_symptoms = nigam_df[\"Clinical observation\"].values.tolist()\n",
    "    for i in range(len(all_symptoms)):\n",
    "        symptoms = all_symptoms[0 : i+1]\n",
    "        prob_positive = compute_prob_positive(prevalence, symptoms, nigam_df)\n",
    "        num_symptoms.append(i + 1)\n",
    "        probs_positive.append(prob_positive)\n",
    "    return num_symptoms, probs_positive\n",
    "\n",
    "prevalences = [0.03, 0.1, 0.3, 0.5]\n",
    "for p_d in prevalences:\n",
    "    xs, ys = generate_prob_positives(p_d, nigam_df)\n",
    "    plt.plot(xs, ys, label=\"p(d)={:.3f}\".format(p_d))\n",
    "    \n",
    "plt.xlabel(\"first k symptoms\")\n",
    "plt.ylabel(\"P(+|symptoms)\")\n",
    "plt.legend(loc=\"best\")\n",
    "_ = plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What is the symptom profile for patient positive for COVID-19?\n",
    "\n",
    "We generate some large number of \"patients\" as bags-of-symptoms by throwing a \"dice\" 50 times (once for each observation) and including a symptom if the dice has a value higher than the probability. Then using the bag of symptoms, we compute the P(+|.) value. We (somewhat arbitrarily) choose the median of these scores as the cutoff, and consider as positive any bag-of-symptom that scores P(+|.) above the cutoff, and then tabulate their symptoms. Results seem believable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['cough', 'sore throat', 'hypertension']"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def get_random_bag_of_symptoms(nigam_df):\n",
    "    all_symptoms = nigam_df[\"Clinical observation\"].values.tolist()\n",
    "    ps_obs = nigam_df[\"P(observation)\"].values.tolist()\n",
    "    dice = np.random.uniform(size=len(ps_obs))\n",
    "    selected_syms = dice < ps_obs\n",
    "    symptoms = [s for b, s in zip(selected_syms, all_symptoms) if b]\n",
    "    return symptoms\n",
    "\n",
    "get_random_bag_of_symptoms(nigam_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEHCAYAAABfkmooAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAXTElEQVR4nO3df7AdZ33f8fcHCdsiwcXG164qWUi0gmB7ikGqq4SGAVxqAWnktHgikmBB3SpxHUr6Y4pMOi2djDL2P0yqNjZRIbGcpngEAaxADKjiZ4rBXINAyEa1gh1bY9USEMABRlTi2z/Oo3CQrrQr6Z5zr3zfr5kzZ/d79sejHcsf7bO7z6aqkCTpZJ420w2QJM1+hoUkqZNhIUnqZFhIkjoZFpKkToaFJKnT/JluwKhcdNFFtXTp0vHv+P77B98rVox/35J0hu6///6vV9XEsfWnbFgsXbqUycnJ8e84GXzPxL4l6Qwl+Yup6nZDSZI6GRaSpE6GhSSpk2EhSepkWEiSOhkWkqROhoUkqZNhIUnq9JR9KO9MLN3wodNe95Ez2MYjt7zmtPcrSaPkmYUkqZNhIUnqZFhIkjqNLCySPD/JzqHPd5L8RpILk2xP8lD7vmBonZuT7E2yJ8k1Q/UVSXa13zYlR0frkySNw8jCoqr2VNWVVXUlsAL4HvB+YAOwo6qWAzvaPEkuA9YClwOrgduSzGubux1YDyxvn9Wjarck6Xjj6oa6GvjzqvoLYA2wpdW3ANe26TXAXVV1qKoeBvYCVyVZCJxfVfdWVQF3Dq0jSRqDcYXFWuDdbfqSqtoP0L4vbvVFwGND6+xrtUVt+ti6JGlMRh4WSc4Bfh54T9eiU9TqJPWp9rU+yWSSyYMHD55aQyVJJzSOM4tXAV+oqifa/BOta4n2faDV9wGXDq23GHi81RdPUT9OVW2uqpVVtXJi4ri3AkqSTtM4wuJ1/KgLCmAbsK5NrwPuHqqvTXJukmUMLmTf17qqnkyyqt0Fdf3QOpKkMRjpcB9JngG8EvjVofItwNYkNwCPAtcBVNXuJFuBB4DDwE1VdaStcyNwB7AAuKd9JEljMtKwqKrvAc8+pvYNBndHTbX8RmDjFPVJ4IpRtFGS1M0nuCVJnQwLSVInw0KS1MmwkCR1MiwkSZ0MC0lSJ8NCktTJsJAkdTIsJEmdDAtJUifDQpLUybCQJHUyLCRJnQwLSVInw0KS1MmwkCR1MiwkSZ0MC0lSJ8NCktTJsJAkdRppWCR5VpL3JvlqkgeT/HSSC5NsT/JQ+75gaPmbk+xNsifJNUP1FUl2td82Jcko2y1J+nGjPrP4L8CHq+qngBcCDwIbgB1VtRzY0eZJchmwFrgcWA3clmRe287twHpgefusHnG7JUlDRhYWSc4HXgq8C6CqflBV3wLWAFvaYluAa9v0GuCuqjpUVQ8De4GrkiwEzq+qe6uqgDuH1pEkjcEozyyeCxwE/iDJF5O8M8lPAJdU1X6A9n1xW34R8NjQ+vtabVGbPrZ+nCTrk0wmmTx48OD0/mkkaQ4bZVjMB14M3F5VLwK+S+tyOoGprkPUSerHF6s2V9XKqlo5MTFxqu2VJJ3AKMNiH7Cvqj7X5t/LIDyeaF1LtO8DQ8tfOrT+YuDxVl88RV2SNCYjC4uq+r/AY0me30pXAw8A24B1rbYOuLtNbwPWJjk3yTIGF7Lva11VTyZZ1e6Cun5oHUnSGMwf8fbfBPxRknOArwFvZBBQW5PcADwKXAdQVbuTbGUQKIeBm6rqSNvOjcAdwALgnvaRJI3JSMOiqnYCK6f46eoTLL8R2DhFfRK4YnpbJ0nqyye4JUmdDAtJUifDQpLUybCQJHUyLCRJnQwLSVInw0KS1MmwkCR1MiwkSZ0MC0lSJ8NCktTJsJAkdTIsJEmdDAtJUifDQpLUybCQJHUyLCRJnQwLSVInw0KS1GmkYZHkkSS7kuxMMtlqFybZnuSh9n3B0PI3J9mbZE+Sa4bqK9p29ibZlCSjbLck6ceN48zi5VV1ZVWtbPMbgB1VtRzY0eZJchmwFrgcWA3clmReW+d2YD2wvH1Wj6HdkqRmJrqh1gBb2vQW4Nqh+l1VdaiqHgb2AlclWQicX1X3VlUBdw6tI0kag1GHRQEfTXJ/kvWtdklV7Qdo3xe3+iLgsaF197XaojZ9bF2SNCbzR7z9l1TV40kuBrYn+epJlp3qOkSdpH78BgaBtB5gyZIlp9pWSdIJjPTMoqoeb98HgPcDVwFPtK4l2veBtvg+4NKh1RcDj7f64inqU+1vc1WtrKqVExMT0/lHkaQ5bWRhkeQnkjzz6DTwj4CvANuAdW2xdcDdbXobsDbJuUmWMbiQfV/rqnoyyap2F9T1Q+tIksZglN1QlwDvb3e5zgf+Z1V9OMnnga1JbgAeBa4DqKrdSbYCDwCHgZuq6kjb1o3AHcAC4J72kSSNycjCoqq+Brxwivo3gKtPsM5GYOMU9UngiuluoySpH5/gliR1MiwkSZ0MC0lSJ8NCktTJsJAkdeoMiyQXjqMhkqTZq8+ZxeeSvCfJqx0aXJLmpj5h8TxgM/B6YG+S307yvNE2S5I0m3SGRQ1sr6rXAf+cwRAd9yX5ZJKfHnkLJUkzrvMJ7iTPBn6FwZnFE8CbGIzjdCXwHmDZKBsoSZp5fYb7uBf4Q+Daqhp+r8RkkneMplmSpNmkT1g8v72h7jhVdes0t0eSNAv1ucD90STPOjqT5IIkHxlhmyRJs0yfsJioqm8dnamqv+RHr0KVJM0BfcLiSJK/fkdpkudwgteaSpKemvpcs/hN4M+SfLLNv5T2nmtJ0tzQGRbt7XYvBlYBAf51VX195C2TJM0afd+Udy7wzbb8ZUmoqk+NrlmSpNmkz0N5twK/COwGftjKBRgWkjRH9DmzuJbBsxaHRt0YSdLs1OduqK8BTz/dHSSZl+SLST7Y5i9Msj3JQ+37gqFlb06yN8meJNcM1Vck2dV+2+Tot5I0Xn3C4nvAziS/1/5HvSnJplPYx5uBB4fmNwA7qmo5sKPNk+QyYC1wObAauC3JvLbO7QzuwFrePqtPYf+SpDPUJyy2Ab8FfAa4f+jTKcli4DXAO4fKa4AtbXoLg26uo/W7qupQVT0M7AWuSrIQOL+q7m3Djtw5tI4kaQz63Dq7JckCYElV7TnF7f8O8O+BZw7VLqmq/W3b+5McfRp8EfDZoeX2tdr/a9PH1o+TZD3tGZAlS5ZMtYgk6TT0ea3qPwZ2Ah9u81cm2dZjvZ8DDlRVr7MQBs9wHKtOUj++WLW5qlZW1cqJiYmeu5UkdelzN9TbgKuATwBU1c4kfd5h8RLg55O8GjgPOD/J/wCeSLKwnVUsBA605fcBlw6tvxh4vNUXT1GXJI1Jn2sWh6vq28fUOseGqqqbq2pxVS1lcOH6Y1X1Kwyugaxri60D7m7T24C1Sc5tYbQcuK91WT2ZZFW7C+r6oXUkSWPQ58ziK0l+CZiXZDnwrxhc7D5dtwBbk9wAPApcB1BVu5NsBR4ADgM3VdWRts6NwB3AAuCe9pEkjUmfsHgTg8EEDwHvBj7C4O6o3qrqE/yoG+sbwNUnWG4jsHGK+iRwxansU5I0ffrcDfU9BmHxm6NvjiRpNuozNtTHmeIaRVW9YiQtkiTNOn26of7d0PR5wD9lcE1BkjRH9OmGOvY5if899CIkSdIc0Kcb6sKh2acBK4C/ObIWSZJmnT7dUPfzoyepDwMPAzeMslGSpNmlTzdUn6e1JUlPYX26of7JyX6vqvdNX3MkSbNRn26oG4CfAT7W5l/O4AG7bzPonjIsJOkprk9YFHDZ0WHF2+B/v1tVbxxpyyRJs0afgQSXHg2K5gngeSNqjyRpFupzZvGJJB9hMC5UMRhB9uMjbZUkaVbpczfUryf5BeClrbS5qt4/2mZJkmaTPmcWAF8Anqyq/5XkGUmeWVVPjrJhkqTZo89rVf8F8F7g91ppEfCBUTZKkjS79LnAfRODV6R+B6CqHgIuHmWjJEmzS5+wOFRVPzg6k2Q+PV6rKkl66ugTFp9M8lZgQZJXAu8B/mS0zZIkzSZ9wmIDcBDYBfwq8KfAfxhloyRJs8tJwyLJPODOqvrvVXVdVb22TXd2QyU5L8l9Sb6UZHeS/9zqFybZnuSh9n3B0Do3J9mbZE+Sa4bqK5Lsar9tSpIz+DNLkk7RScOiqo4AE0nOOY1tHwJeUVUvBK4EVidZxeBMZUdVLQd2tHmSXMbggb/LgdXAbS2sAG4H1gPL22f1abRHknSa+jxn8QiDt+NtA757tFhVbz/ZSu3s46/a7NPbp4A1wMtafQuDQQnf0up3VdUh4OEke4GrkjwCnF9V9wIkuRO4FrinR9slSdPghGcWSf6wTf4i8MG27DOHPp2SzEuyEzgAbK+qzwGXHB1rqn0fvQ13EfDY0Or7Wm1Rmz62Lkkak5OdWaxI8hzgUeC/ns7GWzfWlUmeBbw/yRUnWXyq6xB1kvrxG0jWM+iuYsmSJafYWknSiZwsLN4BfBhYBkwO1cPgf9bP7buTqvpWkk8wuNbwRJKFVbW/DXd+oC22D7h0aLXFwOOtvniK+lT72QxsBli5cqXPgkjSNDlhN1RVbaqqFwB/UFXPHfosq6rOoEgy0c4oSLIA+IfAV4FtwLq22Drg7ja9DVib5NwkyxhcyL6vdVU9mWRVuwvq+qF1JElj0GfU2RtPc9sLgS3tjqanAVur6oNJ7gW2JrmBQRfXdW0/u5NsBR4ADgM3tW4sgBuBO4AFDC5se3Fbksao76izp6yqvgy8aIr6N4CrT7DORmDjFPVJ4GTXOyRJI9TnCW5J0hxnWEiSOhkWkqROhoUkqZNhIUnqZFhIkjoZFpKkToaFJKmTYSFJ6mRYSJI6GRaSpE6GhSSpk2EhSepkWEiSOhkWkqROhoUkqZNhIUnqZFhIkjoZFpKkToaFJKnTyMIiyaVJPp7kwSS7k7y51S9Msj3JQ+37gqF1bk6yN8meJNcM1Vck2dV+25Qko2q3JOl4ozyzOAz826p6AbAKuCnJZcAGYEdVLQd2tHnab2uBy4HVwG1J5rVt3Q6sB5a3z+oRtluSdIyRhUVV7a+qL7TpJ4EHgUXAGmBLW2wLcG2bXgPcVVWHquphYC9wVZKFwPlVdW9VFXDn0DqSpDEYyzWLJEuBFwGfAy6pqv0wCBTg4rbYIuCxodX2tdqiNn1sfar9rE8ymWTy4MGD0/lHkKQ5beRhkeQngT8GfqOqvnOyRaeo1UnqxxerNlfVyqpaOTExceqNlSRNaaRhkeTpDILij6rqfa38ROtaon0faPV9wKVDqy8GHm/1xVPUJUljMsq7oQK8C3iwqt4+9NM2YF2bXgfcPVRfm+TcJMsYXMi+r3VVPZlkVdvm9UPrSJLGYP4It/0S4PXAriQ7W+2twC3A1iQ3AI8C1wFU1e4kW4EHGNxJdVNVHWnr3QjcASwA7mkfSdKYjCwsqurPmPp6A8DVJ1hnI7BxivokcMX0tU6SdCp8gluS1GmU3VA6RUs3fGjG9v3ILa+ZsX1Lmv08s5AkdTIsJEmdDAtJUifDQpLUybCQJHUyLCRJnQwLSVInw0KS1MmwkCR1MiwkSZ0MC0lSJ8NCktTJsJAkdTIsJEmdDAtJUifDQpLUybCQJHUaWVgk+f0kB5J8Zah2YZLtSR5q3xcM/XZzkr1J9iS5Zqi+Ismu9tumJCd6r7ckaURGeWZxB7D6mNoGYEdVLQd2tHmSXAasBS5v69yWZF5b53ZgPbC8fY7dpiRpxEYWFlX1KeCbx5TXAFva9Bbg2qH6XVV1qKoeBvYCVyVZCJxfVfdWVQF3Dq0jSRqTcV+zuKSq9gO074tbfRHw2NBy+1ptUZs+ti5JGqPZcoF7qusQdZL61BtJ1ieZTDJ58ODBaWucJM114w6LJ1rXEu37QKvvAy4dWm4x8HirL56iPqWq2lxVK6tq5cTExLQ2XJLmsnGHxTZgXZteB9w9VF+b5NwkyxhcyL6vdVU9mWRVuwvq+qF1JEljMn9UG07ybuBlwEVJ9gH/CbgF2JrkBuBR4DqAqtqdZCvwAHAYuKmqjrRN3cjgzqoFwD3tI0kao5GFRVW97gQ/XX2C5TcCG6eoTwJXTGPTJEmnaLZc4JYkzWKGhSSpk2EhSepkWEiSOhkWkqROhoUkqZNhIUnqZFhIkjqN7KE8nV2WbvjQjOz3kVteMyP7lXRqPLOQJHUyLCRJnQwLSVInw0KS1MmwkCR1MiwkSZ0MC0lSJ8NCktTJh/I0o2bqYUDwgUDpVHhmIUnqZFhIkjqdNWGRZHWSPUn2Jtkw0+2RpLnkrLhmkWQe8LvAK4F9wOeTbKuqB2a2ZTqbOXii1N9ZERbAVcDeqvoaQJK7gDWAYaGzjhf1dTY6W8JiEfDY0Pw+4O8fu1CS9cD6NvtXSfac5v4uAr5+Oivm6MStP3eauz6rnPZxmmNmzXHKrTPdgpOaNcdpFhvHMXrOVMWzJSwyRa2OK1RtBjaf8c6SyapaeabbearzOPXjcerH49RtJo/R2XKBex9w6dD8YuDxGWqLJM05Z0tYfB5YnmRZknOAtcC2GW6TJM0ZZ0U3VFUdTvLrwEeAecDvV9XuEe7yjLuy5giPUz8ep348Tt1m7Bil6riuf0mSfszZ0g0lSZpBhoUkqdOcDYuu4UMysKn9/uUkL56Jds60Hsfpl9vx+XKSzyR54Uy0c6b1HY4myd9LciTJa8fZvtmiz3FK8rIkO5PsTvLJcbdxNujx9+5vJPmTJF9qx+mNI29UVc25D4OL5H8OPBc4B/gScNkxy7wauIfBMx6rgM/NdLtn6XH6GeCCNv0qj9PUx2louY8Bfwq8dqbbPRuPE/AsBiMzLGnzF890u2fpcXorcGubngC+CZwzynbN1TOLvx4+pKp+ABwdPmTYGuDOGvgs8KwkC8fd0BnWeZyq6jNV9Zdt9rMMnoGZa/r89wTwJuCPgQPjbNws0uc4/RLwvqp6FKCq5uKx6nOcCnhmkgA/ySAsDo+yUXM1LKYaPmTRaSzzVHeqx+AGBmdjc03ncUqyCPgF4B1jbNds0+e/p+cBFyT5RJL7k1w/ttbNHn2O038DXsDg4eRdwJur6oejbNRZ8ZzFCPQZPqTXECNPcb2PQZKXMwiLfzDSFs1OfY7T7wBvqaojg38Mzkl9jtN8YAVwNbAAuDfJZ6vq/4y6cbNIn+N0DbATeAXwt4HtST5dVd8ZVaPmalj0GT7EIUZ6HoMkfxd4J/CqqvrGmNo2m/Q5TiuBu1pQXAS8OsnhqvrAeJo4K/T9e/f1qvou8N0knwJeCMylsOhznN4I3FKDixZ7kzwM/BRw36gaNVe7ofoMH7INuL7dFbUK+HZV7R93Q2dY53FKsgR4H/D6Ofavv2Gdx6mqllXV0qpaCrwX+JdzLCig39+7u4GfTTI/yTMYjC794JjbOdP6HKdHGZx9keQS4PnA10bZqDl5ZlEnGD4kya+139/B4I6VVwN7ge8xSPI5pedx+o/As4Hb2r+aD9ccGzm053Ga8/ocp6p6MMmHgS8DPwTeWVVfmblWj1/P/55+C7gjyS4G3VZvqaqRDl3ucB+SpE5ztRtKknQKDAtJUifDQpLUybCQJHUyLCRJnQwL6Qy0EWR3JvlKkve0ZwNIsiDJJ5PM67mdtyV5wxT1iXYrqTSjDAvpzHy/qq6sqiuAHwC/1ur/jMGAeEeGF07yhiRv67vxqjoI7E/ykulqsHQ6DAtp+nwa+Dtt+pcZPI08HT7QtifNGMNCmgZJ5jN4n8euNkTDc6vqkWna/CTws9O0Lem0zMnhPqRptCDJzjb9aeBdDAYK/NbRBZI8G9jRZi8EzklybZt/fVXt6tjHAeBvTV+TpVNnWEhn5vtVdeVwIcn3gfOOzreReK9sv70BWFpVbzuFfZwHfP+MWyqdAbuhpGnW3hw4L8l5nQv38zxgTg2mp9nHsJBG46Oc5ougkqxM8s6h0suBD01Lq6TT5Kiz0ggkeRHwb6rq9T2XfxvwSFXdMcVvnwLWDL3rXBo7zyykEaiqLwIf7/tQ3okkmQDeblBopnlmIc0CSV4GfKuqdnYtK80Ew0KS1MluKElSJ8NCktTJsJAkdTIsJEmdDAtJUqf/D5DpIGxVDZTGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ps_pos = []\n",
    "num_examples = 10000\n",
    "for i in range(num_examples):\n",
    "    symptoms = get_random_bag_of_symptoms(nigam_df)\n",
    "    p_pos = compute_prob_positive(prevalence, symptoms, nigam_df)\n",
    "    ps_pos.append(p_pos)\n",
    "\n",
    "median = np.median(ps_pos)\n",
    "\n",
    "plt.hist(ps_pos)\n",
    "plt.axvline(x=median, color=\"red\", linewidth=2)\n",
    "plt.xlabel(\"P(+|.)\")\n",
    "plt.ylabel(\"frequency\")\n",
    "_ = plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('cough', 3860),\n",
       " ('febrile', 3276),\n",
       " ('dyspnea', 3119),\n",
       " ('sore throat', 1051),\n",
       " ('chest pain', 868),\n",
       " ('fatigue', 822),\n",
       " ('chills', 799),\n",
       " ('myalgia', 733),\n",
       " ('malaise', 686),\n",
       " ('headache', 642)]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "symptom_counts = collections.Counter()\n",
    "num_examples = 10000\n",
    "for i in range(num_examples):\n",
    "    symptoms = get_random_bag_of_symptoms(nigam_df)\n",
    "    p_pos = compute_prob_positive(prevalence, symptoms, nigam_df)\n",
    "    if p_pos > median:\n",
    "        for symptom in symptoms:\n",
    "            symptom_counts[symptom] += 1\n",
    "\n",
    "symptom_counts.most_common(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (bayes)",
   "language": "python",
   "name": "bayes"
  },
  "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.8.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}