triSBM_survival.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2f343eb0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Author: Filippo_Valle\n",
      "\n",
      "pandas    : 1.3.3\n",
      "numpy     : 1.19.5\n",
      "seaborn   : 0.11.2\n",
      "matplotlib: 3.5.1\n",
      "lifelines : 0.26.4\n",
      "\n",
      "Compiler    : GCC 7.3.0\n",
      "OS          : Linux\n",
      "Release     : 5.13.0-28-generic\n",
      "Machine     : x86_64\n",
      "Processor   : x86_64\n",
      "CPU cores   : 12\n",
      "Architecture: 64bit\n",
      "\n",
      "Git hash: 576214878cace9abdbc3f729d37bf445b11e12ff\n",
      "\n",
      "Git repo: git@github.com:fvalle1/keywordTCGA.git\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%load_ext watermark\n",
    "%watermark -a Filippo_Valle -ppandas,numpy,seaborn,matplotlib,lifelines -g -m -r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "0d73c33a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import os\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "7a44336e",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_dir = \"brca\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ef6ae6b",
   "metadata": {},
   "source": [
    "## Survival"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d4f31804",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "Index: 1092 entries, TCGA-3C-AAAU to TCGA-A2-A25C\n",
      "Data columns (total 27 columns):\n",
      " #   Column                                           Non-Null Count  Dtype \n",
      "---  ------                                           --------------  ----- \n",
      " 0   file_miRNA                                       1074 non-null   object\n",
      " 1   cases.0.demographic.days_to_birth                1058 non-null   object\n",
      " 2   cases.0.demographic.days_to_death                148 non-null    object\n",
      " 3   cases.0.demographic.gender                       1073 non-null   object\n",
      " 4   cases.0.demographic.vital_status                 1073 non-null   object\n",
      " 5   cases.0.diagnoses.0.age_at_diagnosis             1058 non-null   object\n",
      " 6   cases.0.diagnoses.0.days_to_last_follow_up       970 non-null    object\n",
      " 7   cases.0.diagnoses.0.last_known_disease_status    1073 non-null   object\n",
      " 8   cases.0.diagnoses.0.tumor_stage                  1073 non-null   object\n",
      " 9   cases.0.project.project_id                       1074 non-null   object\n",
      " 10  sample_submitter_id_mirna                        1074 non-null   object\n",
      " 11  file_fpkm                                        1092 non-null   object\n",
      " 12  cases.0.demographic.days_to_birth.1              1076 non-null   object\n",
      " 13  cases.0.demographic.days_to_death.1              151 non-null    object\n",
      " 14  cases.0.demographic.gender.1                     1091 non-null   object\n",
      " 15  cases.0.demographic.vital_status.1               1091 non-null   object\n",
      " 16  cases.0.diagnoses.0.age_at_diagnosis.1           1076 non-null   object\n",
      " 17  cases.0.diagnoses.0.days_to_last_follow_up.1     987 non-null    object\n",
      " 18  cases.0.diagnoses.0.last_known_disease_status.1  1091 non-null   object\n",
      " 19  cases.0.diagnoses.0.tumor_stage.1                1091 non-null   object\n",
      " 20  cases.0.project.project_id.1                     1092 non-null   object\n",
      " 21  sample_submitter_id                              1092 non-null   object\n",
      " 22  sample_submitter_id_type                         1092 non-null   object\n",
      " 23  sample_submitter_id_mirna_type                   1074 non-null   object\n",
      " 24  cancer.type                                      1092 non-null   object\n",
      " 25  Subtype_mRNA                                     1092 non-null   object\n",
      " 26  Subtype_Selected                                 1092 non-null   object\n",
      "dtypes: object(27)\n",
      "memory usage: 238.9+ KB\n"
     ]
    }
   ],
   "source": [
    "df_files = pd.read_csv(\"files.dat\", sep=\",\", index_col=0).replace(\"unknown\", np.nan)\n",
    "df_files.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "2fe8576c",
   "metadata": {},
   "outputs": [],
   "source": [
    "bins = np.linspace(0,100,20)\n",
    "avg = df_files[\"cases.0.diagnoses.0.age_at_diagnosis\"].astype(float).mean(skipna=True)\n",
    "df_files[\"age_at_diagnosis\"]=pd.cut(df_files[\"cases.0.diagnoses.0.age_at_diagnosis\"].astype(float).fillna(avg)/365., bins=bins, labels = (bins[1:]+bins[:-1])/2).__array__()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d5ad1143",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_files[\"gender\"] = (df_files[\"cases.0.demographic.gender\"]==\"male\").astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "3e2494b7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['stage x', 'stage iii', 'stage i', 'stage ii', nan, 'stage iv',\n",
       "       'not reported'], dtype=object)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "for letter in [\"a\", \"b\", \"c\"]:\n",
    "    for (old, new) in zip([\"stage %s%s\"%(i,letter) for i in [\"i\", \"ii\", \"iii\", \"iv\"]],[\"stage %s\"%i for i in [\"i\", \"ii\", \"iii\", \"iv\"]]):\n",
    "        df_files.replace(old, new, inplace=True)\n",
    "df_files[\"cases.0.diagnoses.0.tumor_stage\"].unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "365f2ecb",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_files[\"tumor_stage\"]=df_files[\"cases.0.diagnoses.0.tumor_stage\"]\n",
    "for i,stage in enumerate([\"stage i\", \"stage ii\", \"stage iii\", \"stage iv\", \"stage v\", \"stage x\"]):\n",
    "    df_files[\"tumor_stage\"].replace(stage, i+1, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "8a5a8f9e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#0 = Alive\n",
    "df_files[\"vital_status\"]=(df_files[\"cases.0.demographic.vital_status\"]==\"Dead\").astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "dc0e41d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_survival(case):\n",
    "    if case[\"cases.0.demographic.vital_status\"] == 1:\n",
    "        return case[\"cases.0.demographic.days_to_death\"]\n",
    "    else:\n",
    "        return case[\"cases.0.diagnoses.0.days_to_last_follow_up\"]\n",
    "\n",
    "df_files[\"days_survival\"] = df_files.apply(get_survival,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a1530ba-11c9-447f-bbc1-37df1943a43b",
   "metadata": {},
   "source": [
    "## Clusters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "7614dc79-b2c6-447a-84bb-926faa199c89",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/conda/lib/python3.8/site-packages/lifelines/plotting.py:710: UserWarning: For less visual clutter, you may want to subsample to less than 25 individuals.\n",
      "  warnings.warn(\"For less visual clutter, you may want to subsample to less than 25 individuals.\")\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*******Cluster 6 Z:2.854441119394784*******\n",
      "0.18181818181818177 0.37679591836734694 0.06830679926251397\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAALICAYAAABijlFfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABo3UlEQVR4nO3deXgeV3nw/++NQmSHOBvB7GAHCG1YUozCWtNQNwEal6UUSEsIlCW0oL4/lrAYaKtSaFMwL+3btCVpy5IaGloKpuACCS60YmukBChpSGxANjUhNk5wrDiSYqT798eM4omiXSONHun7ua65nlnOOc89z3iIbs6cM5GZSJIkSZLm7h5NByBJkiRJS4UJliRJkiTVxARLkiRJkmpigiVJkiRJNTHBkiRJkqSamGBJkiRJUk1MsCRJkiSpJiZYkiRJklQTEyxJagER8esR8e8RcSAihiJiR0S8KyJOLo9/OCJ6a/7OF0bEy+psc4rvOyoi3hoRO8tz3BMR71+o76/TbK7HRL/3fFzbuZrJv42IODUiuiLi6HmK5Tci4msRcXNEDEbEDRHxjvn6PkmaylFNByBJmlxEvA94HfAh4P3AQeA04HeARwHPm6evfiFwMvDheWp/rA8BG4A/Aq4HHkxxnsvFRL/3HwMrFzyayc3k38azgNdkZtc8xXJv4EvAe4EDwBOALuB+QOc8fackTcgES5IWsYj4NeANwCsy84OVQ/8REZcCZzcT2cxFRBvQlpl3jHPsmcC5wOmZed2CB7eIZeb3m45hjh4LXDtfjWfmJWN2fSkijgNeGxG/l5k5X98tSePxEUFJWtxeD1wzJrkCIDOHM/Nz41WKiC9HxCfG7DszIjIiHl3Z96iI+HxE3BIRhyLiuxHx2oj4MPB84JfKOhkRXZV6vxgR/xERt5ePZv1tRKwa830fjojeiHhuRPwPMAg8cYLzfDnw7zNNrsZ8x/XlI2JfiYjTxpR7YUR8p3z08H8j4t0RcdRM25ru7zpOnE+OiH+NiBvL3/lbEfHi6nczwe893iOCU51P5VzOioj/Lr/zKxHxqCl+z0njnCrWcdr7CcW1fXql7C9PFkNNbgZ8RFBSI+zBkqRFKiLuCTwFeN88fs2/UjyOdx4wBDwSOI7isbSHACcArynL7injeiqwHdgK/AbFI1oXASeW21VrgPcA7wT2An0TxPFE4F8j4mLgfIr/Pn0e6MzMG6c4h4cC/xf4fWCA4hHDL0TEIzJzMCLOBj4OXAa8iaJH5Y/LuH9nJm1NEcdUMX4V+ABFovlU4EMRMZKZ/8gkv/dYMzifh1A8Nvfu8lw2A/8UEY+epFdnqjiZbqwREcCvUTy+9z7gs+Whb41Trm2CeO6UmT+b7HjZQ9oOrAP+D/A39l5JaoIJliQtXvem+IPxh/PReBQTZJwCPDczv1Pu3l45fgtwj8z8xpiqFwFfy8wXVcr+CNhe/vFefRzs3sCvZOa3pgjnfsDLgG9TPCq4iiIx+1REPGmKP5RPBp6TmV8rY7ka+H7Z3gcokrsvZ+ZLy/KfL/6m508j4l2ZuWcGbc1KZl4+ul4mFP8JPAh4FfCPmfn9SX7vsaZ7PicBT83MneX33gP4FEUSff1s4izLTCvWzMyIuAlYAXx2krIvpRh/N5WY4vghivsFjiSfkrTgfERQkha/+fp/4W8B/hf4QES8KCJWT1UhIo4BnkzRE3LU6AJ8BTgMPH5MlR9NI7mC4o/noEhu/i0zPw68hGLCgqkeKds3mhABZOZu4GrgCWWvxjrgn8fU+TjFfwOfPN22pnEOE4qIEyPi/0XEborf6TBwAXDqDNuZyfnsGk2uSqOPXz5ovuOseCzFv9/JxmB9BjhjGstUngKsB94IPAe4eJYxS9Kc2IMlSYvXzRSP7T1kPhrPzJHycbN3Ax8EVkbEV4H/k5nfnKDaiRSPc/11uYz14DHbe6cZzk+BH2TmzZV9XwHuoJhJcPu4tQr7Jth3f4oeqXuOE8fo9kkzaGsuPgw8ieLxuusoZoL8XYpEYCZmcj4HxpQZnVxkxQLEOeqxFNf1tknK3ALcOsv275SZ15SrX4mI/cBHIuJ9S2CSEEktxgRLkhapzDxcJjzPAN4xw+qD3H2Q/9hkgsy8Hnh+Od5rPfBnwLaImKiX4wBFj0QX8G/jHB87Xmq6vW/f5cjjXVUBjExRd7yet9XA/wD7KXphxpa5b/l5ywzagmn+rlURsQI4h2I82Qcq+2fzFMlMz2faao5z1GMpHvucTF2PCFaNJltrKR7xlKQF4yOCkrS4/TnQEREvHXsgIu4RxfTm49kD/NyYfWdN9CWZeTgz/51igof7U0xgcAdjejsy8xDwDeCRmdk7zjLVhBQT+Szw2HJc2KinUfTWTPUH+uqIeMroRkQ8hOIxuqsyc5jiEb8XjKnzQorE7evTbavcNaPftdRO0es3VGl3FfDsMeXu9nuPNYvzmYnpxjmtWEs/xwTjvSrqekSw6qnl50STqkjSvLEHS5IWscz8TET8X+Dvy9n7Pg3cRvGH6+8Auyhm2xvrU8ArIuL9wDbg6RQ9YXeKiMdSzCz3ceAHFI//vQX4dmbeEhHXA8+JiOdSJBY3lgnUmykmtBgBPgH0UzzGeA7w9szcMYtTvZRi5rfPRMSfUExy8WfAFzPzK1PU3Q/8Q0SMzvz3TorH+j5cHv9DipkAPwRcDjyG4hG4vx0zwcV02prydx0rM2+NiB7gDyLiIEUi9FaKx+KOqxSd6PceaybnM20ziHMmsR4EnhYRT6NIyv5r7IQl5WOhN49Td1oi4vPAFyl6GYcpkqs3Ah8ffTwwIs6kmM3w6Zn55dl+lyRNhz1YkrTIZeYbgRcBjwA+BlxJ8QfkdorxMePV2Qa8jWLa9E9RTL/9ujHFbqIYu/N24HMUY6q+y5Eei78GrqAYn9VDMdkBZcLzNOA+wD9Q9EC8mWLCjOmOuRob70GKySx+SpE0/FV5fi+cRvXdFDPGdZV1DwLPGJ1WPTOvoJiZsKOM9XUU04Z3zqKt6fyu4/ktit6Uy4C/AP6lXK8a9/cea4bnM1PTiXPasVIkaCdRTus/T9Om91DM8vjPwD9RTA2/iWKSlFHHlJ/jjbGTpFqFr4iQJLWq8qW3j87MjsXUlhaXiPgj4GmZ+fSmY5G09NmDJUmSlrqnUIwvlKR55xgsSZK0pGXmVBORSFJtfERQkiRJkmriI4KSJEmSVBMTLEmSJEmqiQmWJEmSJNXEBEuSJEmSamKCJUlqVES8LCIyIo6dh7YzIqZ8AW9EvCoidkbEYERcHREbplGnIyI+HBE3RMRI+R6t8co9IiL+JSL2RsTBiPhaRDxzmvE/PiJ+GhHHldu7ynOaaPnydNqd4Lt+JSI+HhG7I+L2iLg2Ijojom0adac8x4jomiDmZ1bKdETEzRFx/GzPQ5Ka5jTtkqSmbQOeDNzexJdHxLnAB4Au4CvAbwOfjYgzMvPaSao+FfhF4BvAqgnaXgVcCfwU+F3gNuAC4DMR8dTMvGqK8N4FfCAzD5bbzwPaxyn3a8DbgH+bor3JXAAcA7wD+F+Kc3sfsBZ440SVZniOtwJjk8vvjq5kZm9EfAt4PcX1kKSW4zTtkqTaRUQA7Zk52HAcCfxeZl48SZkbgK9m5svL7XsA3wa+nZnnTVLvHpk5Uq73Atdm5svGlHkm8DngsZn5nXLfUcCPgA9n5lsmaf8RwA7g1MzcOUm5+wP/DXwH+JXRmGYqIk7OzP1j9v0JRbJzQmYOTVBvWucYEV1AZ2aePEUc5wGbgQdl5s9mcy6S1CQfEZSkZS4inl0+FneofBztvyLil8pja8rHuDaOqfPhMqkY3e6KiP0R8YsR0QMMAi8q23zNON/ZGxH/UK7f5RHBiOiLiPeMU+cTEdFdrt8rIi4uH8+7vazzV6OP0s3g3E8BTgX+aXRfmaD8M/CsyepOM5G5Z/l5a6Xez4BDQExR96XAf0+RXAXwEYr/nr9ktslVGdf+cXZ/E1gBTPa7zuUcx/OvwEnAM2ZRV5IaZ4IlSctYRDwM+ATw7xSPmb0Y+CzFH7gzdQzFH/t/R/EY2DfKtl405jtPAR4PfHyCdv4JeGGZPIzWORb41UqdY4A24O0UidDvA79MkRjNxM+Vn9eP2f9d4KSIuM8M2xtrO7AL2BwRD46IkyLibcBq4MNT1N0AfG2KMm8AzgJemZk/qh6IiLaIOGqKZaq/A54C7M/Mn0xSZibneEKZiB+OiG9GxK+Pbax8HPJ/gF+ZIjZJWpQcgyVJy9vjgP7MfFNl32zH8awE3pCZnx7dERGXA5+IiAdk5o3l7hdRjNe5YoJ2LgfeDDyRIkmDIvlrp0ygyj/4f7fyPUcBfcBXIuIhmfnDacZ8Yvl5YMz+n1aOT5ZcTCozb4+IMyl+09GYDgLPyczrJqpXJpePA7ZMUuYXgD8BLs3MT41T5PvAQ6cI8Y+YYKxTRJwG/A7wZ5M1MINz/B7Fdf0WcCzwauBfIuL5mfnJMc1+G3jCFLFL0qJkgiVJy9t3gOMj4iPARynGIh2aZVtJMRan6nMUkx68APiLct+LgE9l5h3jNpL5zYjYUZb7RqXOlzNz72i5iHgJRQ/OI4B7VZo4lSN/6M8k9qqYYP+MRMS9KJLCnwLPoZjI48UUicXTM/ObE1Q9kSKhHO+xPSJiJcX1+gHFGKnxjCalk7lxvJ0RcSLwLxRju/5ksgame46ZuWVMvc9Q9ND9ATA2wdoPrJ8idklalEywJGkZy8wbIuI5wFspeiAOR8SngP9visfCxvPTsUlTZg5GxKcpEqS/iIhHAqcDbxqvgYqPAy+PiDdQzND3TOD3Rg9GxPOAy4C/oZg97xbg/sCnKMYMTTvm8vMEKmOIym24e8/WTL0COI1iwobRtr5Y/g5/BDx7gnqj5zDuxBIUs/s9HHhSZk40++J1TD0G6m5jtiJiBfBpiuTs2RMlwhWzOsfMzIj4JPBnEdGWmcOVw0PM7DpK0qLhGCxJWuYyc1tmrgfuTfHH8q8Af1keHp0F8Ogx1cYbozVRb8/HgSdFxEMoEq2fUIz5mszlwAMppgp/LsV4q2ovxwuA/8rM12Tm5zLzvziSLM3E6Nirnxuz/+eAW2aRZI71c8DuSuIx6lvAwyapd3P5ecLYA+WEI78LbJqkBwyKRwQPT7H8wZi224CPAY8CnlXtMZzEbM9x1Hj/bk6gSJolqeXYgyVJAiAzbwU+Vs4g+ORy9z6KP8R/frRcOeHEk4Hd02z6Cork54UUCdYnxvRWjBfLdRFxbVl+LXBlZt5cKbKSu/fuvHia8VS/5wfl44gvAL4Ad07T/gLu/rjjbOwG1kTEiZlZTQAfTzExxERxDUXEDynO/U4RcT/ggxS/6fun+O7ZPCL41xS9hb+SmTdMUXfUrM6xHGf2PIrp8Mf+e1hDMUW9JLUcEyxJWsYi4tUUydLnKf7YfgRFcnEZFFORl4/4vT4idlM8MvdGYGC635GZo48dvoHiMb67Tds+gY8D/x9wPPCqMceuBP4qIt4O/BfFDIMbphvTGF3AlojYBXyVYnr0RwC/NVqgTDq3Axsy8z/KffcBfqksciLw0Ij4DYDM/ES5/2OULwAup56/HTiPYgKHu0x9P46vUiQpVR+i+D3+BnhiZaLFUUOVMU/fmaL9uyhn/rsA+FNgJCKeVDl83ejLjiNie9n+6O89rXOMiP+gGNd1PcWYuVcBT6LooRyrgykm15CkxcoES5KWt/+mGCPzfyke+/sx8Lfc9dGxTuBSit6NnwLvppi++9Ez+J7LKR4/vBHonkGdP6boqdo65tglwCkUCdgKioTrtzgyKca0ZeY/lr1yb6GY7v1/gI2ZeW2lWFA8pljNaB7FXaeFPwU4s1KezPzfiHg6xW92SRnr9cBvZOa2KUL7JPChiFiZmaMJ7TPLz/FmDYSyN2mKdidydvm5qVyqng58uVxvqx6YwTl+D3gdRZI9AlwDnJOZd+kpjIjHAffh7hNfSFJLiMw5TZAkSZLmQUQcDewBXpuZM32/V8uKiD8FzshM34MlqSWZYEmStEhFxJso3if1i03HshDKKd93U/R+fbnhcCRpVnxEUJKkxeti4JiIOL6chGSpewjwTpMrSa3MHixJkiRJqonvwZIkSZKkmviIYM1OPvnkXLNmTdNhSJIkSZpHV1999f7MvM/Y/SZYNVuzZg29vb1NhyFJkiRpHpXvh7wbHxGUJEmSpJqYYEmSJElSTUywJEmSJKkmJliSJEmSVBMTLEmSJEmqiQmWJEmSJNXEBEuSJEmSamKCJUmSJEk1McGSJEmSpJqYYEmSJElSTY5qOoClZvctt/Pqy69pOgxJkiSpFpecu67pEFqKPViSJEmSVBN7sGr20JOOMcuXJEmSlil7sCRJkiSpJvZgSZIkSQJgzVu3TavcrovOmedIWpc9WJIkSZJUExMsSZIkSaqJCZYkSZIk1cQES5IkSZJqYoIlSZIkSTUxwZIkSZKkmphgSZIkSVJNTLAkSZIkqSa+aLhmu2+5nVdffk3TYUiSJElqgD1YkiRJklQTe7Bq9tCTjuGSc9c1HYYkSZI0Y2u+ta3pEFqePViSJEmSVBMTLEmSJEmqiQmWJEmSJNXEBEuSJEmSamKCJUmSJEk1McGSJEmSpJqYYEmSJElSTUywJEmSJKkmJliSJEmSVBMTLEmSJEmqiQmWJEmSJNXkqKYDWGp233I7r778mqbDkCRJktQAe7AkSZIkqSb2YNXsoScdwyXnrms6DEmSJGnG1nxrW9MhtDx7sCRJkiSpJiZYkiRJklQTEyxJkiRJqsmSTLAi4q0RkaPLFGXvGxHvi4gbImIgIm6JiO6IeGVExELFLEmSJKn1LblJLiLikcAfTrPs44EvAPcud90GrAJ+sVxeEBHPzsyh+YhVkiRJ0tKypHqwIuIewN8DK4CvT1H2eOCzFMnV9cAZmbkKuBfQCRwGzgbeP58xS5IkSVo6llSCBfwe8FTgo8AVU5S9ELgfMAD8amb2AmTmHZn5VxzpBbsgIk6dp3glSZIkLSFLJsGKiLXAu4GbgddPo8r55eflmdk3zvG/pHhksA14cS1BSpIkSVrSlkyCBfwtxeN9b8jMn0xWsByn9ZBy83PjlcnM24DucvPsuoKUJEmStHQtiQQrIl4FbAC+mJmXTaPKoyvr105SbvTYabONTZIkSdLy0fKzCEbEA4H3UoylevU0qz2gsv6jScqNHjsuIo4te7UmtfuW23n15ddMMwxJkiRJS0nLJ1jAJcDxwFsy8wfTrLOqsn77JOWqx1ZRjMm6m4i4ALgA4PgHPXyaIUiSJEmLyzN+4f5323fJuesaiKR1tXSCFRHnAecA3wL+b1NxZOalwKUAHR0d6T9CSZIkaXlq2TFYEbEa+HNgGHhVZv5sBtX7K+vHTFKueqx/wlKSJEmSRGv3YP0ZxUuC/wa4PiKOHXP86NGVyrE7MvMO4MZKuQcCByf4jgeWnwenM/5KkiRJ0vLWsj1YwNry83cpepfGLpsqZUf3vafcrs4cWJ1RcKzRY9fNNVhJkiRJS18rJ1izlpk3AD8sN585XpmIuBewvty8YiHikiRJktTaWjbByswzMzMmWoA/qpQd3f+6ShOj78s6NyLWjPMVrwWOpRjj9dF5Og1JkiRJS0jLJlg12AzcRDGRxbaIeDxARBwdEb8L/HFZ7tLM3NFQjJIkSZJaSCtPcjEnmXlrRGwEvgCcBvRGRD+wArhnWewK4PUNhShJkiSpxSznHiwy82rgUcD7gZ0UidUh4CvAq4BnZeZQcxFKkiRJaiWRmU3HsKR0dHRkb29v02FIkiRJmkcRcXVmdozdv6x7sCRJkiSpTst2DNZ82X3L7bz68muaDkOSJEmqxSXnrms6hJZiD5YkSZIk1cQerJo99KRjzPIlSZKkZcoeLEmSJEmqiT1YkiRJkgBY89Zt0yq366Jz5jmS1mUPliRJkiTVxARLkiRJkmpigiVJkiRJNTHBkiRJkqSamGBJkiRJUk1MsCRJkiSpJiZYkiRJklQTEyxJkiRJqokvGq7Z7ltu59WXX9N0GJIkSZIaYA+WJEmSJNXEHqyaPfSkY7jk3HVNhyFJkiTN2JpvbWs6hJZnD5YkSZIk1cQES5IkSZJqYoIlSZIkSTUxwZIkSZKkmphgSZIkSVJNTLAkSZIkqSYmWJIkSZJUExMsSZIkSaqJCZYkSZIk1cQES5IkSZJqYoIlSZIkSTU5qukAlprdt9zOqy+/pukwJEmSJDXAHixJkiRJqok9WDV76EnHcMm565oOQ5IkSZqxNd/a1nQILc8eLEmSJEmqiQmWJEmSJNXEBEuSJEmSamKCJUmSJEk1McGSJEmSpJqYYEmSJElSTUywJEmSJKkmJliSJEmSVBMTLEmSJEmqiQmWJEmSJNXkqKYDWGp233I7r778mqbDkCRJktQAEyxJkiRJADzjF+5/t32XnLuugUhalwlWzR560jH+I5QkSZKWKcdgSZIkSVJNTLAkSZIkqSYmWJIkSZJUk5ZPsCJiXUT8YUT8a0RcHxE3R8Th8vOrEfH2iDhpgrpdEZHTWB6+0OclSZIkqfUshUkuXg68trI9CAwAJwFPKZfXRcSzM/PrE7RxGLhlku/4WR2BSpIkSVralkKCdRWwC/gKcH1mHgCIiGOB5wPvBe4DbI2IUzPz1nHa+Fpmnrkg0UqSJElaslo+wcrMyybYfxvwkYj4MfAFYDWwEfjoAoYnSZIkaRlp+TFY0/CNyvqDGotCkiRJ0pK3HBKs9ZX17zcWhSRJkqQlb0kmWBHRHhFrIqIT+Idy9/eAz0xQ5VERcW1EDETEbRFxQ0T8bUQ8bmEiliRJkmZueHCIns4utp2+kZ7OLoYHh5oOadlr+TFYVRExCLSPc+irwG9l5kT/4k6mmHXwAHAccGq5vCIi/iQz3zHdGPp39LH9rPNnFLckSZI0Gwd37mJw734YSQ7t3gPAGRd3NRvUMrfUerBuAvYChyr7vgS8LjN/OE75ncCbgUcCKzLz3sC9gGcAVwMBvD0i3jjZl0bEBRHRGxG9dxw+XMNpSJIkSVM7fKAfRhKA4YEh9nX3NhyRIjObjmFeRMRq4CXA24ETgHdl5h/MoP4K4D+BM4DbgAdNMMX7XXR0dGRvr/+wJUmSNP96Orvo27KV4YEh2la2s/a859qDtUAi4urM7Bi7f6n1YN0pM/dl5vuAZwIJ/H5EbJxB/UHgbeXmscCG+qOUJEmSZm/d5k2sPe+5HPdzD2Ptec9l3eZNTYe07C2pMVjjycyrIuIrwNOAC4DPzqD61yvrp9QamCRJkjRHbSvaa+2xWvPWbdMqt+uic2r7zqVmyfZgjfGj8vPhjUYhSZIkaUlbLgnWaO9T/wzrPamy3ldTLJIkSZKWqJZOsCKiLSJiijIbgCeUm1+u7J+qXjvw7nLzELB99pFKkiRJWg5aOsECHgx8MyJeHRGnVJOmiHhwRLwV+DTFdOu3AO+v1H1aRHwxIs6LiAdV6t2zTMq6gSeWu9+ZmQfm+2QkSZIktbalMMnF6cAHyvU7IuIgsJLifVaj+oDnZ+ZNlX1BMTPgBoCIGKDoqToeuGdZZgS4KDPfM3/hS5IkSVoqWj3BuhF4IXAmRW/T/YGTgWHgh8C3KXqwPpaZA2Pqfge4EHgy8Jiy3gnA7cB1FD1Yl2bmd+b7JCRJkiQtDS2dYGXmHcA/l8tM694MvK/2oCRJkiQtW60+BkuSJEmLxPDgED2dXWw7fSM9nV0MDw41HZK04Fq6B2sx6t/Rx/azzm86DEmSpAV3cOcuBvfuh5Hk0O49ALW+BFdqBfZgSZIkqRaHD/TDSAIwPDDEvu7ehiOSFp49WDVbdepaNlx5WdNhSJIkLbiezi76tmxleGCItpXtrF7f0XRI0oIzwZIkSVIt1m3eBMC+7l5Wr++4c1taTkywJEmSVIu2Fe2OudKy5xgsSZIkSaqJCZYkSZIk1cQES5IkSZJqYoIlSZIkSTUxwZIkSZKkmphgSZIkSVJNTLAkSZIkqSYmWJIkSZJUExMsSZIkzdnw4BA9nV1sO30jPZ1dDA8ONR2S1Iijmg5gqenf0cf2s85vOgxJkqQFdXDnLgb37oeR5NDuPQCccXFXs0FJDbAHS5IkSXN2+EA/jCQAwwND7OvubTgiqRn2YNVs1alr2XDlZU2HIUmStKB6Orvo27KV4YEh2la2s3p9R9MhSY0wwZIkSdKcrdu8CYB93b2sXt9x57a03JhgSZIkac7aVrQ75krCMViSJEmSVBsTLEmSJEmqiQmWJEmSJNXEBEuSJEmSamKCJUmSJEk1McGSJEmSpJqYYEmSJElSTUywJEmSdKfhwSF6OrvYdvpGejq7GB4cajokqaX4ouGa9e/oY/tZ5zcdhiRJ0qwc3LmLwb37YSQ5tHsPgC8QlmbABEuSJEl3OnygH0YSgOGBIfZ19zYckRbSrovOaTqElmeCVbNVp65lw5WXNR2GJEnSrPR0dtG3ZSvDA0O0rWxn9fqOpkOSWooJliRJku60bvMmAPZ197J6fced25KmxwRLkiRJd2pb0e6YK2kOnEVQkiRJkmpigiVJkiRJNTHBkiRJkqSamGBJkiRJUk1MsCRJkiSpJiZYkiRJklQTEyxJkiRJqokJliRJ0gwNDw7R09nFttM30tPZxfDgUNMhSVokfNFwzfp39LH9rPObDkOSJM2jgzt3Mbh3P4wkh3bvAfDlvJIAe7AkSZJm7PCBfhhJAIYHhtjX3dtwRJIWC3uwarbq1LVsuPKypsOQJEnzqKezi74tWxkeGKJtZTur13c0HZKkRcIES5IkaYbWbd4EwL7uXlav77hzW5JMsCRJkmaobUW7Y660JK1567Zpldt10TnzHEnrcgyWJEmSJNXEBEuSJEmSatLyCVZErIuIP4yIf42I6yPi5og4XH5+NSLeHhEnTdHGfSPifRFxQ0QMRMQtEdEdEa+MiFioc5EkSZLU2pbCGKyXA6+tbA8CA8BJwFPK5XUR8ezM/PrYyhHxeOALwL3LXbcBq4BfLJcXlHV9g6AkSZKkSbV8DxZwFfAm4MnAiZm5MjOPo0iSXgb8BDgZ2BoRx1crltufpUiurgfOyMxVwL2ATuAwcDbw/oU5FUmSJEmtrOUTrMy8LDM3Z+Y3MvNAZf9tmfkR4Lxy12pg45jqFwL3o+jx+tXM7C3r3pGZfwX8YVnugog4dT7PQ5IkSVLra/kEaxq+UVl/0Jhj55efl2dm3zh1/5LikcE24MXzEJskSZKkJWQ5JFjrK+vfH12JiEcCDyk3Pzdexcy8DeguN8+el+gkSdKiMDw4RE9nF9tO30hPZxfDgw6/ljRzS2GSi7uJiHbg/hSPBL6z3P094DOVYo+urF87SXPXAs8CTpvOd/fv6GP7WedPXVCSJC0qB3fuYnDvfhhJDu3eA+DLhCXN2JJKsCJiEGgf59BXgd8aMxPgAyrrP5qk2dFjx0XEsWWv1tjvvQC4AGDNilUzC1qSJC0Khw/0w0gCMDwwxL7u3oYjktSKllSCBdwErACOpZgJEOBLwJsz84djylYzodsnabN6bBXFmKy7yMxLgUsBOjo6csOVl80wbEmS1LSezi76tmxleGCItpXtrF7f0XRIklrQkkqwMnPN6HpErAZeArwduCoi3pWZf9BUbJIkaXFbt3kTAPu6e1m9vuPObUmaiSWVYFVl5j7gfRHRDXwd+P2IuCozP1sW6a8UPwY4OEFTx1TW+ycoI0mSWlzbinbHXEmasyU/i2BmXgV8pdy8oHLoxsr6AydpYvTYwfHGX0mSJEnSqCWfYJVGJ6p4eGVfdebA6oyCY40eu67WiCRJkiQtOcslwTql/LzzEb/MvAEYnfjimeNVioh7ceQ9WlfMW3SSJEmSloSWTrAioi0iYooyG4AnlJtfHnN4dLq/cyNizTjVX0sxI+Ew8NHZRypJkiRpOWjpBAt4MPDNiHh1RJxSTbYi4sER8Vbg00AAtwDvH1N/M8XU7scA2yLi8WXdoyPid4E/Lstdmpk75vlcJEmSJLW4pTCL4OnAB8r1OyLiILCSI+/BAugDnp+ZN1UrZuatEbER+AJwGtAbEf0U79K6Z1nsCuD18xi/JEmSpCWi1XuwbgReCPw1cDWwHziO4rx+CHwGeCXwqMz85ngNZObVwKMoerd2UiRWhyhmHnwV8KzMHJrf05AkSZK0FLR0D1Zm3gH8c7nMpZ29wBvKRZIkSZJmpdV7sCRJ0jI1PDhET2cX207fSE9nF8ODPnAiqXkt3YO1GPXv6GP7Wec3HYYkSUvewZ27GNy7H0aSQ7v3AHDGxV3NBiVp2bMHS5IktaTDB/phJAEYHhhiX3dvwxFJkj1YtVt16lo2XHnZ1AUlSdKc9HR20bdlK8MDQ7StbGf1+o6mQ5IkEyxJktSa1m3eBMC+7l5Wr++4c1uSmmSCJUmSWlLbinbHXEladByDJUmSJEk1McGSJEmSpJqYYEmSJElSTUywJEmSJKkmJliSJEmSVBMTLEmSJEmqiQmWJEmSJNXEBEuSJDVieHCIns4utp2+kZ7OLoYHh5oOSZLmzBcN16x/Rx/bzzq/6TAkSVr0Du7cxeDe/TCSHNq9B8AXB0tqeSZYkiSpEYcP9MNIAjA8MMS+7t6GI5K066Jzmg6h5Zlg1WzVqWvZcOVlTYchSdKi19PZRd+WrQwPDNG2sp3V6zuaDkmS5swES5IkNWLd5k0A7OvuZfX6jju3JamVmWBJkqRGtK1od8yVpCXHWQQlSZIkqSYmWJIkSZJUExMsSZIkSaqJCZYkSZIk1cQES5IkSZJqYoIlSZIkSTUxwZIkSZKkmphgSZKkBTc8OERPZxfbTt9IT2cXw4NDTYckSbXwRcM169/Rx/azzm86DEmSFrWDO3cxuHc/jCSHdu8B8KXDkpYEe7AkSdKCO3ygH0YSgOGBIfZ19zYckSTVwx6smq06dS0brrys6TAkSVrUejq76NuyleGBIdpWtrN6fUfTIUlSLUywJEnSglu3eRMA+7p7Wb2+485tSWp1JliSJGnBta1od8yVtAiteeu2aZXbddE58xxJ63IMliRJkiTVxARLkiRJkmqyoAlWRDwkIi6NiN0RMRgRN0bExyLi9IWMQ5IkSZLmw5wTrIg4KSJui4jbI+L6iIgJyj0GuAZ4BfBg4GjgfsCLgKsi4ty5xiJJkiRJTaqjB+tXgGOAduDDmZljC5RJ18eAk8apH8A9gY9ExGk1xCNJkiRJjagjwTq7sv6pCcr8OvAoYDT5+gjwPOClwHfLfUcB76ohHkmSJElqRB3TtP9C+XlrZt4wQZmXVNYvzczfHd2IiM8D1wH3Bs6JiOMz89Ya4pIkSZKkBVVHD9ZDKHqmvjfewYg4CthQ2fW+6vHM/Anw9+XmUcAZNcQkSZIWqeHBIXo6u9h2+kZ6OrsYHhxqOiRJqk0dPVirys+Jep06gHtRJGHXZeZ4idg3KuuPAL5YQ1yN6N/Rx/azzm86DEmSFq2DO3cxuHc/jCSHdu8B8KXDkpaMOnqwRttom+D4L1bWvzRBmZ9U1o+fc0SSJGnROnygH0aKYdnDA0Ps6+5tOCJJqk8dPVi3UoyfesAEx3+5sv6fE5Q5urJ+t1kIW8mqU9ey4crLmg5DkqRFq6ezi74tWxkeGKJtZTur13c0HZIk1aaOBGsH8BTg4RFxv8y8afRARBwPPL1S9ssTtLG6su4EF5IkLWHrNm8CYF93L6vXd9y5LUlLQR0JVjdFghXAHwGvrhx7G8X7sRK4JjP3T9DG6ZX1XTXEJEmSFqm2Fe2OuZK0ZNWRYH0IuJBiLNYrI+Lnga9QJE3PrJT7u0na+KXK+rU1xCRJkiRJC27OCVZm7oiI9wJvpeipemq5VF3LkanY7yIiHgw8qay7JzP3zDUmSZIkSWpCHbMIkplvA94N/IziUcHq0gtszMyfTVD9NWU5gC/M9Lsj4t4R8dsRsSUirouIQxExFBF7ImJrRDxvkrpdEZHTWB4+07gkSZIkLT91PCIIQGb+fkT8JXA28CBgEOjNzK9MUXUl8JFyfdxerincxF3PYxA4DDywXJ4TEZ8DfiMzb5+gjcPALZN8x0TJoSRJkiTdqbYECyAz9wFbZljndXP82qOAq4APA1/IzB8ARMQa4B3AK4BnAZcAL5mgja9l5plzjEOSJEnSMldrgtWQX87Mu73AODN3UUy68TOKmQ3Pi4i3Zeb/LnSAkiRJkpaHOY/BiogPlsvr59BGZ9nGjB8RHC+5GqPapm8ylCRJkjRv6pjk4mXAS4Gz5tDGhrKdl809nLsZrKy3zUP7kiRJkgTUNIvgIndmZf07E5R5VERcGxEDEXFbRNwQEX8bEY9bgPgkSWppw4ND9HR2se30jfR0djE8ONR0SJLUmMU2BivrbCwiTgA2lZvdmXnDBEVPBk4CDgDHAaeWyysi4k8y8x3T/c7+HX1sP+v8WccsSVKrObhzF4N798NIcmh38TrLMy7uajYoSWrIYunBOq78nGga9RmLiHsA/wDcHxgCfm+cYjuBNwOPBFZk5r2BewHPAK6meD/X2yPijVN81wUR0RsRvXccPlzXKUiS1BIOH+iHkeL/Ix0eGGJfd2/DEUlScxrvwYqIo4FfoOi92ltj038BbCzXX5OZ3x5bIDM/Os6+O4ArIuI/gf8EzgC6IuLvMvPW8b4oMy8FLgXo6OjIDVdeVtMpSJK0+PV0dtG3ZSvDA0O0rWxn9XrnlJK0fM0owYqIh0xyeOUUx+/SFMULhh8O/C5wIkWCdc1M4pmw8YjNQGe5+frM/OBM28jMwYh4G3AlcCzFRByfrCM+SZKWknWbi6fx93X3snp9x53bkrQczbQHaxfjj5MK4GlA3xzj+dgc6xMR7wFGH+l7U2b++Rya+3pl/ZQ5tCNJ0pLVtqLdMVeSVJrtI4IxzX0zsSUzPz2XBiLivcCF5eabM3PzHGOSJEmSpGmbTYI110QKil6wQ8BNFI8Ffiwz/3UuDZaPBY72XL05M987txABeFJlfa69c5IkSZKWuBklWJl5t1kHI2KEImH6Qmb+al2BzcSY5OrCzHzfNOpEZk44LXxEtAPvLjcPAdvnHKgkSZKkJa2uadrr6NWa3RdH/BlHkqs3TCe5Kj0tIr4YEedFxIMq7d0zIjYA3cATy93vzMwDtQUtSZIkaUmqY5r23y4/f1RDWzNSzlr45nJzBHhLRLxlkiqbK+OygmJmwA1lWwMUPVXHA/estHlRZr6n7tglSZIkLT1zTrAy8yN1BDJL9xizft8pyh9bWf8OxYQYTwYeA5wMnEDxsuPrKHqwLs3M79QVrCRJkqSlrfEXDc9FZu5ilo8nZubNwHQfJ5QkSZKkKdU1BkuSJC2g4cEhejq72Hb6Rno6uxgeHGo6JEkS89CDFREPBZ4C/DxwInAM0+9lysx8Rd0xLaT+HX1sP+v8psOQJC1xB3fuYnDvfhhJDu3eA+DLfiVpEagtwYqI04H/C5w5x6ZaOsGSJGkhHD7QDyPF20aGB4bY193bcESSloJdF53TdAgtr5YEKyKeD3ysbG8uU7ZP+F6qVrHq1LVsuPKypsOQJC1xPZ1d9G3ZyvDAEG0r21m9vqPpkCRJ1JBglY8E/gPF1OajCdIPga8Ceyhm5ZMkSTVat3kTAPu6e1m9vuPObUlSs+rowXo9sIIiuboNeEVmfqKGdiVJ0gTaVrQ75kqSFqE6EqyzKusvzszP1tCmJEmSJLWcOqZpfzBF79UekytJkiRJy1kdCdbouKvv1dCWJEmSJLWsOhKs/y0/j6mhLUmSJElqWXUkWF+kmJr9URFxzxrakyRJkqSWVEeCdQnwM+BewEtqaE+SJEmSWtKcE6zM/C7wdoperPdFxOPmHJUkSZIktaA6erDIzPcCvw8cB3wlIt4REavraFuSpPkyPDhET2cX207fSE9nF8ODQ02HJElqcXN+D1ZE/Htl86fAScAfAX8UEX3ATcAd02wuM3PDXGNqUv+OPrafdX7TYUiSpuHgzl0M7t0PI8mh3XsAfHmvJGlO6njR8JkcmaqdynoAa8tlOmJMO5IkzavDB/phpPhPz/DAEPu6exuOSJLU6upIsKBIjmayf8ladepaNlx5WdNhSJKmoaezi74tWxkeGKJtZTur13c0HZIkqcXVkWD9dg1tSJK04NZt3gTAvu5eVq/vuHNbkqTZmnOClZkfqSMQSZIWWtuKdsdcSVLFmrdum1a5XRedM8+RtK5aZhGUJEmSJJlgSZIkSVJtTLAkSZIkqSZ1zSJ4FxFxLPAUYB1wH+B4IDLzFfPxfZIkSZK0GNSaYEXEA4E/AM4DVlQPUbzj6m4JVkR8ETitPH5mZu6sMyZJkiRJWii1PSIYERuB/wZeCaykSKpGl8lcBtyvXF5aVzySJEmStNBqSbAi4leATwAnUCRUdwBfBP4c+P4U1f8ZOFSuP7eOeCRJkiSpCXNOsCJiJfAh4Ohy1zbglMw8OzPfAHxvsvqZOQBcQZGY/XxE3HeuMUmSNJnhwSF6OrvYdvpGejq7GB4cajokSdISUccYrJcDD6QYQ3UF8OzMzBm28V/A88r1xwB7a4irEf07+th+1vlNhyFJmsTBnbsY3LsfRpJDu/cA+MJhSVIt6nhE8Ncq6783i+QK4IbK+ilzjEeSpEkdPtAPI8V/roYHhtjX3dtwRJKkpaKOHqxHlZ87MnPSxwEn8dPK+vFzjKdRq05dy4YrL2s6DEnSJHo6u+jbspXhgSHaVrazen1H0yFJkpaIOhKs+1A8HvijObRxjwnWJUmq3brNmwDY193L6vUdd25LkjRXdSRYtwEnAsfMoY0HVNZvnls4kiRNrm1Fu2OuJEnzoo7eopsoZgA8LSKmeufVRJ5WWd8154gkSZIkqQF1JFhfKT9XARtnWjkiTgB+s9y8A/hqDTFJkiRJ0oKrI8H6VGX9fRGxaob1P0iRnCXwb+V7sSRJkiSp5cw5wcrMLwBfo3hM8GHAlyLi1KnqRcSDImIb8JzRpoA/nms8kiRJktSUOia5AHgVxaOCJwCPA74TEZ8D/h2432ihiDgLeASwATgHuCdFYpbAuzPzWzXFI0mSJEkLrpYEKzO/GxG/BmwFTqZInH6NIy8hTopE6vOVatUJMf42M/+wjlgkSZIkqSm1vXMqM78GPBb4NEcSqmoSVd03un8f8KrM/J264pAkSZKkptT1iCAAmXkT8LxyDNZ5wHrgMRSPDt4D+BlFUvUNit6sjzqphSRJkqSlotYEa1Rm7gD+YHS7fD/Wysy8fT6+T5K0dA0PDnHNhX/Kvu5eVq/vYN3mTbStaG86LEmSxjUvCdZYmZnAskiu+nf0sf2s85sOQ5KWjIM7dzG4dz+MJId27wHgjIu7mg1KkqQJ1DYGS5Kk+XD4QD+MJADDA0Ps6+5tOCJJkia2ID1Yy8mqU9ey4crLmg5DkpaMns4u+rZsZXhgiLaV7axe39F0SJIkTcgES5K0qK3bvAngLmOwJElarKaVYEXEv1c2MzM3THBsru7StiRJbSvaHXMlSWoZ0+3BOpMj77HKCY7N1XhtS5IkSVLLmMkkFzHFsbkusxIR946I346ILRFxXUQcioihiNgTEVsj4nnTaOO+EfG+iLghIgYi4paI6I6IV5ZTzEuSJEnSlKbbg/Xbszy2EG7irucxCBwGHlguz4mIzwG/Md57uCLi8cAXgHuXu24DVgG/WC4viIhnZ+bQ/J2CJEmSpKVgWglWZn5kNscWyFHAVcCHgS9k5g8AImIN8A7gFcCzgEuAl1QrRsTxwGcpkqvrgZdkZm9EHA28Cng/cHb5+ZoFOBdJkiRJLWwpvAfrlzPziZn5N6PJFUBm7srMV1IkVgDnRcSDx9S9ELgfMAD8amb2lnXvyMy/Av6wLHdBRJw6v6chSZIkqdW1fIKVmV+aosjfV9bHvjzl/PLz8szsG6fuX1I8MtgGvHh2EUqSJElaLuacYEXEH5TLb82hjReOtjPXeMYxWFlvq3znI4GHlJufG69iZt4GdJebZ89DbJIkSZKWkDp6sLooHqU7bw5t/GalnbqdWVn/TmX90ZX1ayepP3rstLoCkiTd1fDgED2dXWw7fSM9nV0MDzqvkCSpNU13FsGWFBEnAJvKze7MvKFy+AGV9R9N0szoseMi4tiyV2tC/Tv62H7W+ZMVkSSNcXDnLgb37oeR5NDuPQC+XFiS1JKWbIIVEfcA/gG4PzAE/N6YIqsq63ebvn2CY6soxmSN/a4LgAsA1qxYNfawJGkKhw/0w0jxrvnhgSH2dfc2HJEkLU+7Ljqn6RBa3mJJsO5Vfg7U2OZfABvL9ddk5rdrbPsuMvNS4FKAjo6O3HDlZfP1VZK0JPV0dtG3ZSvDA0O0rWxn9fqxcxJJktQaGk+wIiI4Mh5qf01tbgY6y83XZ+YHxynWX1k/Bjg4QXPHTFBHklSTdZuLp7n3dfeyen3HnduSJLWaxhKsiFgBPBz4PxTvokpgzr1MEfEe4I3l5psy888nKHpjZf2BTJxgPbD8PDjV+CtJ0uy0rWh3zJUkaUmYUYIVEcMTHQKeMcnx6frEXCpHxHspXh4M8ObM3DxJ8erMgY8GvjtBudHetevmEpskSZKkpW+m07RH5TMq22P3zXQB+BKwZcZnMPrlxWOB1eTqvZOVL2cU/GG5+cwJ2rwXsL7cvGK2sUmSJElaHmbzHqyYusi0DQJXAa8DnpmZOZtGyuRq9LHAC6dKripGZ6M4NyLWjHP8tcCxwDDw0dnEJkmSJGn5mOkYrLVjtgP4AcX4qf8EXjbNdkaAQ8CBzByZYQx3DSDizziSXL0hM98/g+qbgVdSjAHbFhHnZ+bVEXE08Argj8tyl2bmjrnEKUmSJGnpm1GClZm7x+4rJgEkgIHxjs+niHgI8OZycwR4S0S8ZZIqm6vjsjLz1ojYCHwBOA3ojYh+YAVwz7LYFcDraw9ekiRJ0pJTxyyCf1R+fq+GtmbqHmPW7ztF+WPH7ih7rB4FvIXivVkPpuhduxb4CPDBufaySZIkSVoe5pxgZeYfTV1qfmTmLmoYE5aZe4E3lIskSZIkzcpsJrmQJGnahgeH6OnsYtvpG+np7GJ4cKjpkCRJmjeNvWh4qerf0cf2s85vOgxJWjQO7tzF4N79MJIc2r0HwJcKS5KWrHlJsKKY+eJRwIOA44H26dbNzMumLiVJahWHD/TDSPEWjuGBIfZ19zYckSRJ86fWBCsinkAxjukc4JhZNJEceTdVS1p16lo2XNnSpyBJterp7KJvy1aGB4ZoW9nO6vUdTYckSdK8qS3Bioh3U8zEF0w98UROo4wkaQlYt3kTAPu6e1m9vuPObUmSlqJaEqyI2ARU/4s5SDHN+cOAEznyIuJVwAM5Mp16lp9XlXUkSUtM24p2x1xJUotY89Zt0yq366Jz5jmS1jXnWQQj4oHAH5abCfw/4D6Z+QSKxKk4kPn0zOzIzPsDDwf+BBgoD7cDr8jMp881HkmSJElqSh3TtL8KOJoiufpoZr4uMw9NViEzf5CZ7wBOB/rKz+0RcUIN8UiSJElSI+pIsM6srL9zJhUz8/vAM4Eh4CHAX9QQjyRJkiQ1oo4E6xHl5w8z83sTFYqIccd7lXU+QjHpxQsi4rgaYpIkSZKkBVdHgnUSxeOBu8c5dkdlfeUkbVxRfrYDv1RDTJIkSZK04OpIsEbKz8PjHOuvrD9gkjZ+Ull/0JwjkiRJkqQG1JFg3Vx+Hj/OsRsr66dN0sa9K+vjtSNJkiRJi14dCdYNFOOnHjbOsW9X1n91kjaqx26pISZJkiRJWnB1JFg95ecJEXHKmGNfAH5Wrp8XEU8eWzkizgJ+u7Lrv2qISZI0heHBIXo6u9h2+kZ6OrsYHhxqOiRJklreuDP7zdAXgbeW6+cAfzl6IDP3R8QngHMpJrD4UkR8jCIpS2A98IIyjgS+lZnVXq+W07+jj+1nnd90GJI0pYM7dzG4dz+MJId27wHgjIu7mg1KkqQWV0eC9R8U47DuDbyCSoJVegPwy8B9KF5I/NJyGRXl5xDw6hrikSRNw+ED/TCSAAwPDLGvu7fhiCRJan1zTrAyczgifoNy9r+IWJmZA5XjN0XEmcA/AY/mSEJV9WPgtzKz5f/rvurUtWy48rKmw5CkKfV0dtG3ZSvDA0O0rWxn9fqOpkOSJKnl1dGDRWb+xxTHr4+IXwA2AmcDDwHuSZFYfRn4p8wcrCMWSdL0rNu8CYB93b2sXt9x57YkSZq9WhKs6cjMEeBfy0WS1LC2Fe2OuZIkqWZ1zCIoSZIkSWKee7Aioo3iccATKGYRPAjsy8z98/m9kiRJktSE2hOsiLgP8DLg14HTKRKrsWV+BHQDH87MK+uOQZIkSZKaUNsjghHRFhGbgF3ARcATgBUUswaOXR5E8W6sz0fENyLisXXFIUmSJElNqSXBiogVwGeAdwErR3dPVa1cngB8IyLOqSMWSZIkSWpKXY8Ifgh4JpDl9gDwSeBzwLUULyIeAlYBp1AkVecCjy3rrAD+JSKekJn/XVNMkiRJkrSg5pxgRcTZwIsoEqUAtgGvzswbxyl+M8UjhP8OXBQR5wKXAMcCRwN/Azx1rjFJkiRJUhPqeETw5ZX1zwLPmSC5upvMvJyi52u43PWkiPi5GmKSpGVteHCIns4utp2+kZ7OLoYHh5oOSZKkZaGORwSfVH4m8NryhcLTlplfj4jLOJKoPQm4voa4GtG/o4/tZ53fdBiSlrmDO3cxuHc/jCSHdu8B8KXCkiQtgDp6sO5LkVz9T2b+7yzb2DamPUnSHBw+0A8jxbDY4YEh9nX3NhyRJEnLQx09WLcA9wPm8vLgat2fzi2cZq06dS0brrys6TAkLXM9nV30bdnK8MAQbSvbWb2+o+mQJElaFupIsH4A3B948BzaqNb9/tzCkSSt27wJgH3dvaxe33HntiRJml91JFgfp5j572ERcXpmfnsWbbyw/PwJ8OUaYpKkZa1tRbtjriRJakAdY7A+AvSV638fEatmUjkifhN4NsU4rj/JzOEpqkiSJEnSojTnBCsz+4HnU4zFehzw1Yh4ylT1IuKYiHgnMDpg6cOZ+f/mGo8kSZIkNaWOFw0/rVx9HfB+4NFAd0R8G/g8cC3FC4bvAFYBa4EnAOeU2wCXAx+ptDWuzPzPucYrSZIkSfOljjFYX6Z4vG9UAgH8AnD6JPWiUv5F5TKZpJ54JUmSJGle1JmwBEcSrazsm0g1KZusnCRJkiS1hDoSrB9y12RJkiRJkpalOSdYmbmmhjgkSZIkqeXVMU27JKk0PDhET2cX207fSE9nF8ODQ02HJEmSFpCTRtSsf0cf2886v+kwJDXk4M5dDO7dDyPJod17AHzhryRJy4gJliTV6PCBfhgphqUODwyxr7u34YgkSZq+XRed03QILc8Eq2arTl3Lhisvm7qgpCWpp7OLvi1bGR4Yom1lO6vXdzQdkiRJWkAmWJJUo3WbNwGwr7uX1es77tyWJEnLgwmWJNWobUW7Y64kSVrGnEVQkiRJkmpigiVJkiRJNWn5BCsijomIZ0XEOyLikxGxOyKyXLqmqNtVKTvZ8vAFOh1JkiRJLWwpjMF6AvBvc2zjMHDLJMd/Nsf2JUmSJC0DSyHBAvgpcE1leT9wvxnU/1pmnjkPcUmSJElaRpZCgtWdmSdVd0TERU0FI0mSJGn5avkxWJk53HQMkiRJkgRLIMGSpPkwPDhET2cX207fSE9nF8ODQ02HJEmSWsBSeESwDo+KiGuBhwHDwI+A/wT+OjO/OZOG+nf0sf2s8+chREkL6eDOXQzu3Q8jyaHdewB8gbAkSZqSPViFk4GfB24H2oFTgVcCV0fEu6aqHBEXRERvRPTecfjw/EYqaUEcPtAPIwnA8MAQ+7p7G45IkiS1guXeg7UTeDPwaaAvMw9HxNHAmcCfAI8H3h4RP83M903USGZeClwK0NHRkRuuvGzeA5c0v3o6u+jbspXhgSHaVrazen1H0yFJkqQWsKwTrMz86Dj77gCuiIj/pHhM8AygKyL+LjNvXegYJTVj3eZNAOzr7mX1+o47tyVJkiazrBOsyWTmYES8DbgSOBbYAHyy2agkLZS2Fe2OuZIkLTtr3rptWuV2XXTOPEfSuhyDNbmvV9ZPaSwKSZIkSS3BBEuSJEmSamKCNbknVdb7GotCkiRJUktYtglWRMQUx9uBd5ebh4Dt8x6UJEmSpJa2JBKsiDgxIk4eXThyXsdU90fEsZVqT4uIL0bEeRHxoEpb94yIDUA38MRy9zsz88CCnIwkSZKklrVUZhH8JvDQcfa/qVxGfQR4WbkeFDMDbgCIiAGKnqrjgXuWZUaAizLzPfWHLEmSJGmpWSoJ1mx8B7gQeDLwGOBk4ATgduA6ih6sSzPzO00FKEmSJKm1LIkEKzPXzKLOzcD76o9GkiRJ0nK1JMZgSdJMDA8O0dPZxbbTN9LT2cXw4FDTIUmSpCViSfRgLSb9O/rYftb5TYchaRIHd+5icO9+GEkO7d4DwBkXdzUblCRJWhLswZK07Bw+0A8jCcDwwBD7unsbjkiSJC0V9mDVbNWpa9lw5WVNhyFpEj2dXfRt2crwwBBtK9tZvb6j6ZAkSdISYYIladlZt3kTAPu6e1m9vuPObUmSpLkywZK07LStaHfMlSRJmheOwZIkSZKkmphgSZIkSVJNTLAkSZIkqSYmWJIkSZJUExMsSZIkSaqJCZYkSZIk1cQES5IkSZJqYoIlqRHDg0P0dHax7fSN9HR2MTw41HRIkiRJc+aLhmvWv6OP7Wed33QY0qJ3cOcuBvfuh5Hk0O49AL78V5IktTx7sCQ14vCBfhhJAIYHhtjX3dtwRJIkSXNnD1bNVp26lg1XXtZ0GNKi19PZRd+WrQwPDNG2sp3V6zuaDkmSJGnOTLAkNWLd5k0A7OvuZfX6jju3JUmSWpkJlqRGtK1od8yVJElachyDJUmSJEk1McGSJEmSpJqYYEmSJElSTUywJEmSJKkmJliSJEmSVBMTLEmSJEmqiQmWJEmSJNXEBEvSvBoeHKKns4ttp2+kp7OL4cGhpkOSJEmaN75ouGb9O/rYftb5TYchLRoHd+5icO9+GEkO7d4D4AuGJUnSkmWCJWleHT7QDyMJwPDAEPu6exuOSJIkTWTXRec0HULLM8Gq2apT17LhysuaDkNaNHo6u+jbspXhgSHaVrazen1H0yFJkiTNGxMsSfNq3eZNAOzr7mX1+o47tyVJkpYiEyxJ86ptRbtjriRJ0rLhLIKSJEmSVBMTLEmSJEmqiQmWJEmSJNXEBEuSJEmSamKCJUmSJEk1McGSJEmSpJqYYEmSJElSTUywJNVieHCIns4utp2+kZ7OLoYHh5oOSZIkacH5ouGa9e/oY/tZ5zcdhrTgDu7cxeDe/TCSHNq9B8AXDEuSpGXHHixJtTh8oB9GEoDhgSH2dfc2HJEkSdLCswerZqtOXcuGKy9rOgxpwfV0dtG3ZSvDA0O0rWxn9fqOpkOSJElacCZYkmqxbvMmAPZ197J6fced25IkScuJCZakWrStaHfMlSRJLW7NW7dNq9yui86Z50hal2OwJEmSJKkmJliSJEmSVJOWT7Ai4piIeFZEvCMiPhkRuyMiy6Vrmm3cNyLeFxE3RMRARNwSEd0R8cqIiHk+BUmSJElLxFIYg/UE4N9mWzkiHg98Abh3ues2YBXwi+Xygoh4dmb61lRJkiRJk2r5HqzST4HtwHuB3wRumk6liDge+CxFcnU9cEZmrgLuBXQCh4GzgffPQ8ySJEmSlpil0IPVnZknVXdExEXTrHshcD9gAPjVzOwDyMw7gL+KiOOAPwEuiIg/z8wdNcYtSZIkaYlp+R6szByeQ/Xzy8/LR5OrMf6S4pHBNuDFc/geSZIkSctAyydYsxURjwQeUm5+brwymXkb0F1unr0QcUmSJElqXcs2wQIeXVm/dpJyo8dOm8dYpMYMDw7R09nFttM30tPZxfCg87lIkiTN1lIYgzVbD6is/2iScqPHjouIY8terQn17+hj+1nnT1ZEWlQO7tzF4N79MJIc2r0HgDMu7mo2KEmSpBa1nHuwVlXWb5+kXPXYqvEKRMQFEdEbEb13HD5cS3DSQjl8oB9GEoDhgSH2dfc2HJEkSVLrWs49WLXJzEuBSwE6Ojpyw5WXNRyRNH09nV30bdnK8MAQbSvbWb2+o+mQJEmSWtZyTrD6K+vHAAcnKHfMBHWkJWHd5k0A7OvuZfX6jju3JUmSNHPLOcG6sbL+QCZOsB5Yfh6cavyV1IraVrQ75kqSJKkmy3kMVnXmwEdPWOrIsevmMRZJkiRJS8CyTbAy8wbgh+XmM8crExH3AtaXm1csRFySJEmSWteyTbBKo7NRnBsRa8Y5/lrgWGAY+OhCBSVJkiSpNS2JBCsiToyIk0cXjpzXMdX9EXHsmKqbgZsoJrLYFhGPL9s7OiJ+F/jjstylmbljIc5FkiRJUutaEgkW8E3gJ5XlweX+N43Zf3G1UmbeCmwEbgZOA3oj4iBwG/DXwNEUjwa+fv5PQZIkSVKrWyoJ1qxl5tXAo4D3AzuBewKHgK8ArwKelZlDzUUoSZIkqVUsiWnaM3PNHOvvBd5QLpIkSZI0K8u+B0t3Nzw4RE9nF9tO30hPZxfDg3bgSZIkSdOxJHqwFpP+HX1sP+v8psOYk4M7dzG4dz+MJId27wHwRbSSJEnSNNiDpbs5fKAfRhKA4YEh9nX3NhyRJEmS1BrswarZqlPXsuHKy6YuuIj1dHbRt2UrwwNDtK1sZ/X6jqZDkiRJklqCCZbuZt3mTQDs6+5l9fqOO7clSZIkTc4ES3fTtqLdMVeSJEnSLDgGS5IkSZJqYoIlSZIkSTUxwZIkSZKkmphgSZIkSVJNTLAkSZIkqSYmWJIkSZJUExMsSZIkSaqJCVYLGR4coqezi22nb6Sns4vhwaGmQ5IkSZJU4YuGa9a/o4/tZ50/L20f3LmLwb37YSQ5tHsPgC8EliRJkhYRE6wWcvhAP4wkAMMDQ+zr7m04IkmSJC0luy46p+kQWp4JVs1WnbqWDVdeNi9t93R20bdlK8MDQ7StbGf1+o55+R5JkiRJs2OC1ULWbd4EwL7uXlav77hzW5IkSdLiYILVQtpWtDvmSpIkSVrEnEVQkiRJkmpigiVJkiRJNTHBkiRJkqSamGBJkiRJUk1MsCRJkiSpJiZYkiRJklQTEyxJkiRJqokJ1iI1PDhET2cX207fSE9nF8ODQ02HJEmSJGkKvmi4Zv07+th+1vlzbufgzl0M7t0PI8mh3XsAfMmwJEmStMjZg7VIHT7QDyMJwPDAEPu6exuOSJIkSdJU7MGq2apT17Lhysvm3E5PZxd9W7YyPDBE28p2Vq/vqCE6SZIkSfPJBGuRWrd5EwD7untZvb7jzm1JkiRJi5cJ1iLVtqLdMVeSJElaUGveum1a5XZddM48R9K6HIMlSZIkSTUxwZIkSZKkmphgSZIkSVJNTLAkSZIkqSYmWJIkSZJUExMsSZIkSaqJCZYkSZIk1cQES5IkSZJqYoK1SA0PDtHT2cW20zfS09nF8OBQ0yFJkiRJmsJRTQew1PTv6GP7WefPuZ2DO3cxuHc/jCSHdu8B4IyLu+bcriRJkqT5Yw/WInX4QD+MJADDA0Ps6+5tOCJJkiRJU7EHq2arTl3Lhisvm3M7PZ1d9G3ZyvDAEG0r21m9vqOG6CRJkiTNJxOsRWrd5k0A7OvuZfX6jju3JUmSJC1eJliLVNuKdsdcSZIkSS3GMViSJEmSVJNln2BFxMsiIqex/ErTsUqSJEla3HxE8IgR4CeTHPdFVJIkSZImZYJ1xP9m5pqmg5AkSZLUupb9I4KSJEmSVBcTLEmSJEmqiQmWJEmSJNXEBOuI+0TE1RFxW0QMRMQPImJLRJy5EF8+PDhET2cX207fSE9nF8ODzqkhSZIktRonuTjiGGAd8FPgXsDacnlxRHwIuCAzfzZVI/07+th+1vkz/vKDO3cxuHc/jCSHdu8B8EXDkiRJUouxBwtuBP4IOB1YkZknUSRbTwW+WJb5beD9EzUQERdERG9E9N5x+PCsgjh8oB9GEoDhgSH2dffOqh1JkiRJzVn2PViZeQVwxZh9w8DXIuIZwCeB5wCviYj/l5k7x2njUuBSgI6Ojtxw5WUzjqOns4u+LVsZHhiibWU7q9d3zOJsJEmSJDVp2SdYk8nMkYi4kCLBugfwa8D/nY/vWrd5EwD7untZvb7jzm1JkiRJrcMEawqZ+b2I2A+cDJwyX9/TtqLdMVeSJElSi3MMliRJkiTVxARrChHxMIreK4C+JmORJEmStLgt6wQrImIax99bbo4An533oCRJkiS1rGWdYAEPjYirIuLVEXHKaMIVEfeIiCcBnwOeV5a9JDNvaCxSSZIkSYuek1zAGeUCMBQR/cAqoL1S5kPA/1nowCRJkiS1luWeYO0Ffg94MvALwH2AE4FBivFWXwM+mJlfbSpASZIkSa1jWSdYmTkAXFwukiRJkjQny30M1rwYHhyip7OLbadvpKezi+HBoaZDkiRJkrQAlnUP1nzo39HHv/7cWQzu3Q8jyaHdewB8ibAkSZK0DJhgzYPDB/phJAEYHhhiX3dvwxFJkiRJU9t10TlNh9DyfESwZqtOXcva855D28piEsK2le2sXt/RcFSSJEmSFoI9WPNg3eZNAOzr7mX1+o47tyVJkiQtbSZY86BtRbtjriRJkqRlyEcEJUmSJKkmJliSJEmSVBMTLEmSJEmqiQmWJEmSJNXEBEuSJEmSamKCJUmSJEk1McGSJEmSpJqYYNXs1uu+R09nF8ODQ02HIkmSJGmBmWDVbHhwiL4tW7nmwj9tOhRJkiRJC8wEax4MDwyxr7u36TAkSZIkLTATrHnQtrKd1es7mg5DkiRJ0gIzwapZ24p21p73XNZt3tR0KJIkSZIW2FFNB7DUHH/awznj4q6mw5AkSZJmbM1bt02r3K6LzpnnSFqXPViSJEmSVBMTLEmSJEmqiQmWJEmSJNXEBEuSJEmSamKCJUmSJEk1McGSJEmSpJqYYEmSJElSTUywJEmSJKkmJlg1u/W679HT2cXw4FDToUiSJElaYCZYNRseHKJvy1auufBPmw5FkiRJ0gIzwZoHwwND7OvubToMSZIkSQvMBGsetK1sZ/X6jqbDkCRJkrTATLBq1rainbXnPZd1mzc1HYokSZKkBXZU0wEsNcef9nDOuLir6TAkSZIkNcAeLEmSJEmqiQmWJEmSJNXEBEuSJEmSamKCJUmSJEk1McGSJEmSpJqYYEmSJElSTUywJEmSJKkmJlg1u/W679HT2cXw4FDToUiSJElaYCZYNRseHKJvy1auufBPmw5FkiRJ0gIzwZoHwwND7OvubToMSZIkSQvMBGsetK1sZ/X6jqbDkCRJkrTATLBq1rainbXnPZd1mzc1HYokSZKkBXZU0wEsNcef9nDOuLir6TAkSZIkNcAeLEmSJEmqiQlWKSJWRURXRHwnIm6LiFsjoici3hgRRzcdnyRJkqTFz0cEgYh4KPBlYE2563agHegolxdHxIbM/GkjAUqSJElqCcu+Bysi2oDPUCRXPwbOysx7AccA5wL9wOOAjzYVoyRJkqTWsOwTLOBlwGPK9edn5hcBMnMkMz8OvLo89qyI2NBAfJIkSZJahAkWvLT8/FJmfn2c45cDfeX6+QsTkiRJkqRWtKwTrIg4Bnhqufm58cpkZgKfLzfPXoi4JEmSJLWmZZ1gAT/Pkd/g2knKjR67X0ScNFmDt173PXo6uxgeHKojPkmSJEktZLknWA+orP9oknLVYw+YsBQwPDhE35atXHPhn84pMEmSJEmtZ7lP076qsn77JOWqx1aNPRgRFwAXAJzMUQwPDLGvu7eeCCVJkqQFsuuic5oOoeUt9x6sWmTmpZnZkZkdq+Io2la2s3p9R9NhSZIkSVpgyz3B6q+sHzNJueqx/glLAW0r2ll73nNZt3nTnAKTJEmS1HqWe4J1Y2X9gZOUqx67ccJSwPGnPZwzLu6ibUX7nAKTJEmS1HqWe4L1XWCkXH/0JOVGj92UmbfMb0iSJEmSWtWyTrAy83bgq+XmM8crExEBPKPcvGIh4pIkSZLUmpZ1glX6SPn59Ih44jjHXwCcUq5ftjAhSZIkSWpFkZlNx9CoiDgKuAZ4DMX7rl6amdsj4h7A84G/A44DPpeZvzqN9vqBG+YxZC2sk4H9TQeh2ng9lxav59Li9VxavJ5Li9dzfA/NzPuM3bnsEyyAiFgDfAlYU+66naJ3b0W5/U1gQ2b+dBpt9Wamc7QvEV7PpcXrubR4PZcWr+fS4vVcWryeM+MjgkBm7gIeC7wTuBZI4DBwNXAh8KTpJFeSJEmSlrejmg5gscjMfuAPy0WSJEmSZswerPpd2nQAqpXXc2nxei4tXs+lxeu5tHg9lxav5ww4BkuSJEmSamIPliRJkiTVxARLkiRJkmpigiVJkiRJNTHBGkdErIqIroj4TkTcFhG3RkRPRLwxIo6eY9v3jYj3RcQNETEQEbdERHdEvDIioq5z0BHzcT3L9nIay8PrPp/lKiKOiYhnRcQ7IuKTEbG78jt31fQd3p8LZD6vp/fnwouIe0fEb0fEloi4LiIORcRQROyJiK0R8bwavsP7c4HM5/X0/lx4EbEuIv4wIv41Iq6PiJsj4nD5+dWIeHtEnDTH7/D+rHCa9jEi4qHAl7nrS4fbgY5yeXFETOulw+O0/XjgC8C9y123AauAXyyXF0TEszNzaC7noCPm83qWDgO3THL8Z7NsV3f3BODf5qtx788FN6/Xs+T9uXBu4q5/UwxS/P4PLJfnRMTngN/IzNtn2rj354Kb1+tZ8v5cOC8HXlvZHgQGgJOAp5TL68p76Oszbdz78+7swaqIiDbgMxR/jP8YOCsz7wUcA5wL9AOPAz46i7aPBz5L8Y/veuCMzFwF3AvopPgfmrOB98/5RATM7/Ws+Fpm3m+SZdfczkJj/BTYDrwX+E2KPwLmzPuzMfNyPSu8PxfOUcBVwGuAh2Xmysw8FlgL/H1Z5lnAJTNt2PuzEfN2PSu8PxfOVcCbgCcDJ5bX8ziKJOhlwE+Ak4Gt5f02bd6fE8hMl3IBXgFkuTx5nOO/WTm+YYZt/3FZ73Zg7TjHN5XHfwac2vRvsRSWeb6eXWW9Lzd9nstlAdrG2bervA5dc2zb+3NpXU/vz4W/nk+f4vgHKv97++AZtu39ubSup/fnIlsoEqDR6/niGdb1/hxnsQfrrl5afn4px+8ivRzoK9fPn2Hbo+Uvz8y+cY7/JUWXahvw4hm2rfHN5/XUAsvM4Xls3vtzgc3z9dQCy8wvTVHk7yvrHTNs3vtzgc3z9dTi843K+oNmWNf7cxwmWKWIOAZ4arn5ufHKZJGKf77cPHsGbT8SeMgUbd8GdM+0bY1vPq+nlhbvT2lBDFbW26Zbyftz0ZrV9dSitb6y/v3pVvL+nJgJ1hE/z5Hf49pJyo0eu98MZlx59Dj1J2v7tGm2q4nN5/WselREXFvOmHNbOXvO30bE42bRlprh/bl0eX8uHmdW1r8zg3ren4vTmZX1mVzPKu/PBkVEe0SsiYhO4B/K3d+jGLs+Xd6fEzDBOuIBlfUfTVKueuwBE5aaW9vHRcSx02xb45vP61l1MkUyNzo74anAK4GrI+Jds2hPC8/7c+ny/lwEIuIEinEYAN2ZecMMqnt/LjJzvJ5V3p8NiIjBiEiKXsg+ikf4TgS+SjEefSYz/Xl/TsAE64hVlfXJphytHls1YamFa1vjm+/ffCfwZuCRwIrMvDfFjDnPAK4GAnh7RLxxBm2qGd6fS4/35yIREfeg+H/H7w8MAb83wya8PxeRGq4neH827SZgL3Cosu9LwOsy84czbMv7cwImWNIsZOZHM/O9mbkjMw+X++7IzCso3vnQUxbtmumUp5LmxvtzUfkLYGO5/prM/HaTwWjO5nw9vT+blZlrspgG/1jgvsCFwC8AV0XEOxsNbgkxwTqiv7J+zCTlqsf6Jyy1cG1rfI395pk5CLyt3DwW2FBHu5o33p/LiPfnwomIzRTvwQF4fWZ+cBbNeH8uEjVdz0l5fy6szNyXme8DnkkxlfrvR8TGKapVeX9OwATriBsr6w+cpFz12I0Tlppb2wfLWVc0e/N5PaejOi38KTW2q/p5fy4/3p/zLCLeA4w+4vWmzPzzWTbl/bkI1Hg9p8P7c4Fl5lXAV8rNC2ZQ1ftzAiZYR3wXGCnXHz1JudFjN2XmLdNsuzqzynTavm6a7Wpi83k9tbR4f0o1ioj3Am8qN9+cmZvn0Jz3Z8Nqvp5avEYnonj4DOp4f07ABKuUmbdTzKACRVfp3UREUAzCBLhiBm3fAIwOHJyo7Xtx5D0E025b45vP6zlNT6qsj/fiPS0S3p/LkvfnPCkfI7uw3HxzZr53Lu15fzar7us5Td6fzRjtLZz2I3zenxMzwbqrj5SfT4+IJ45z/AUc+Qd42QzbHi1/bkSsGef4aymeNx4GPjrDtjW+ebmeZWI22fF24N3l5iFg+3TbVmO8P5cI78/mlH+Mjz5GdmGNf4x7fzZgPq6n9+fCi4i2afzuG4AnlJtfnuFXeH+OJzNdygU4CvhvioF+eyjeBwBFIvoC4Nby2L+NU7erPJbAmnGOHw/8uDz+P8Djy/1HA79LMd1pAn/d9O+wVJb5up7ALwFfBM4DHlTZf0+KAblXVeq+uenfYSktFO/qOLmy/LD8nd8zZv+x072e5XHvzyVyPb0/G7uWf1b5XV8/w7ren4tsma/r6f3ZyLVcA3wLeDXF/6kclWMPBt4K3Fb+5jcD95vu9SyPe3+O97s3HcBiW8p/iH2Vf0yHgIHK9jXAiePUm/QfYFnm8cD+SrmDwB2V7S8A7U3/BktpmY/rSfEG+6wstwM/GXMth4F3N33+S20Bdo357SdaPjzd61kp4/25BK6n92cj1/EhY37bm6ZYLpzu9ayU8f5cAtfT+7OR67lmzG8+VP7mt43Z/wPgcePU9/6cxXIUuovM3BURj6V45vjXgbXAYYqs/B+Bv8zMO2bZ9tUR8SjgLRTvkXgwxR/811I8zvbBzByZpAnN0Dxdz++U7T0ZeAzF/8N+AsV/KK4DuoFLM/M7dZyDFob355Lh/bnw7jFm/b5TlD92pl/g/bmg5vN6en8uvBuBF1Ikt0+keEn0yRSJ7A+BbwOfBj6WmQOz+QLvz7uLMvOUJEmSJM2Rk1xIkiRJUk1MsCRJkiSpJiZYkiRJklQTEyxJkiRJqokJliRJkiTVxARLkiRJkmpigiVJkiRJNTHBkiRJkqSamGBJ0hIWEbsiIsvlzKbjaTUR8eCIOFT+fp9qOh7dXUQcFRE7y2u0MyKObjomScubCZYkSRPbDBwDjADvaDgWjSMzfwb8Qbn5cOB1zUUjSSZYktQSIuLDlZ6orqbjWQ4i4gnAC8vNT2Tm/zQZjyZ1OXBDuf72iDixyWAkLW8mWJIkje+PK+vvaSwKTSkzE3hfuXkccGGD4Uha5kywJGkJy8w1mRnl8uWm42kVEfF44Oxy8+uZeXWT8Wha/gE4UK7/XkSsajAWScuYCZYkSXf3usr63zYVhKYvMweBj5Wbq4DfbjAcScuYCZYkSRXl+J3RsVdDwCcaDEcz89HK+u80FoWkZc0ES5IWsdGJLYCXVnb/YWXCi7ss49Sfcpr2iHhZpcyXK/vPiIi/jYjrI+K2cumJiNeNNxV2RJwUEX8QEb0RcTAiBiPi+xHxgYhYO4tzv19EXBgR2yPih2V7P42I/4mIi8tJKObDC4HR8/tiZvZPEN/xlSncMyJ+bbpfEBEPi4iRSt0zpihf228RESsiYmNEvD8i/jMifly2d3tE7ImIL0TEmyPi3tNsr6tyHh+u7P+liPj7iLguIg6Ux7eOU//kiHhD+b0/ioiBiDhc1vmfiPhERPxORJwyjXC+Dvy4XP/5iFg3nXOQpDod1XQAkqTFJSLuQTHBwyYgxhzuKJcXRsTZmXlbWecs4B+BsX+UnwK8GnhJRDw/Mz8/je8P4C3A24FjxxxuB04ATgNeGxGXAa8uHw+ry/Mr65+bqFBm3hoRH+fIo2ivAj4zze94BUd+229nZs94her+LSLiN4FLKB6hG88Dy+Vs4Pcj4v9k5oemdUZHvuNY4K+Bl0yj7AvKeMab9e/4cjmN8ppExBMz86qJ2svMjIgrgfPLXb8BXDOT+CVprkywJGlx+0L5+RjgAeX694HvzeN3vosiuQK4GbgeSOB0jvxh/mSKqbE3RsTTKBKLdmAY+A7wU+ChFAkWFO+S+peIeFRm7proiyPiKODDwIvHHNoJ/Kj8jkdX4jgfWBMRZ2XmHbM417HfvxJYX9n1H1NUuYQjCdavRsQDMvPGKb6jjbv2SI47xmuefouHcdfk6ifALuAgsKI8fr/y2LHAByPi6My8ZLJzqoZN8Zjes8vtAxT/fg6XbVfP72yKf0PVp2luBPqAgfL7Hwrcv3J8Ok/efIkjCdYzgLdNM3ZJqkdmuri4uLgs8oXiD+0sl64Z1NtVqXfmBGVeVilzC8VLdX8CvAhoq5RbAfxFpWwCGyn+KE7gL4F7j2n7WUB/pfxlU8T77krZEeBi4EFjytyTorfotkrZ99X0O59ZafMQcI9p1PlWpc7bplH+2ZXytwMnLNRvQdET9lWKXsUHTVDmcRSJfTXGh07SZlel7MHycy9wLnBUpVwAp1S2/7tS72vA4yZo/37AK4GrgCdN4/c9bczvdtx83ZcuLi4u4y2OwZIkVZ1I8Qf1mZn58cwcHj2QmYOZ+f9x116df6HoYejKzN/LzJurjWXm57jrjHzPj4h7jffFEfE44K2VXS/PzM7M3DOmzcOZ+bfAM4Gflbv/v4h4yExOdAKPr6xfl5kj06hT7d15eflY32ReWVn/RGYeGFtgHn+LP8/Mp2bmJWPbqrT5TYrE+FPlrpXAa6c4p1GrKBLqX8rMyzNzNCay8IPy/B5M0SsLRXJ4Tvm948VzU2b+XWY+ARj3UcoxdlBMTgJFUvf4ScpKUu1MsCRJY707M/9nkuPVhOJo4LsUjxVOZAvFH91QPCr4uAnKvYkj/136aGZ+eLIgM/MrwKXlZhtFr8xcnVZZ//4063yUorcLisfgnj5RwYi4P0XyMmqiKeDn5bfIzEPj7R+n3EgZw6jnTKde6Z2Zef0UZR5UWf+fzPzpNOMankaZnwE/rOw6baKykjQfTLAkSWP9/RTHx04y8KHJ/vDNzCHg25VdPz+2TDn2qTq5xPunCrL0scr6hmnWmcxDK+s/mk6FzDxIMZZo1CsnKksxXmt0/PMNmdk9tsBi+S0y8/vA/nLzERFxwjSq/QyYzqQY1Yk4HlGec52q1+6hE5aSpHngJBeSpKofZOa+Kcr8eMz2N6bRbrXOCeMcfyJHpka/NTOvnkabANdW1h8XEZGZd5uufgZOrqxPq1eldAnFzIAAvx4RJ2XmLdUC5aODL6/s+rsJ2lqQ36J8jPBXgMcCq4HjuPvfBaMzFwbFJCsHpojhu2MfE52oHMWjqMcAJwGfjIjOMqmrQ/W3v09NbUrStJhgSZKq9k5VIDNvHzPMaMo6FH9MjzpmnOOPqawfFRFTTuc+jqMpkoRbZ1F3VHV82MB0K2VmT0R8k+Lxx3aKKcr/Ykyxp3NkJr3DwEcmaG5ef4uIOBX4c4pxW1ONF6s6fhplppUgZeZgRPwN8MZy1zOBnRHRA3wR6Aa+VvYOzkb12o075k+S5osJliSpajZTnc+0znh/1Fffn3Uvium1Z+N45pZgVc0k+YCiF+sD5foruXuC9YrK+qcz8ycTtDNvv0VE/CLweWaXdLRPo8y4L2WewNsoHt/7jXI7gCeUC8BwRPwXxeOXH84JXvg8gZleO0mqjWOwJEmLQV29DHP971p1EoiZjgv6GMWMeACPjognjh6IiBOBX6+UnWhyC5in3yIiVgGfqLR/G8XU+r8GnErR43XPzIzRBdg9w++czqyLAGTmHZn5AuAcihc6Hx5TpA14CvD/gL6IeNEM4qheu2lN7CFJdbEHS5K0GFR7Wr6RmU9uKI5qr9KJM6mYmf0R8THggnLXK4H/KtfPo3iPGBTvJvviJE3N12/xcuC+5foBindK3TBFnVVTHJ+zzPw34N/K6fufCjwN+CWKl1m3lcXuDfxjRAxl5tZpNFu9dhP1FErSvLAHS5K0GFTHcT28sSju2mPzwFnUr05hf25EjE4SUX088INTvF9rvn6LsyvrfzFVclUmPCfU+P2TysxDmXlFZr4jM9dTJINv4cj4vQA2T7O56rWbaS+cJM2JCZYktYbqH+RLcXxJdSbCkyPi9IbiuK6yPuPkJjOvAUZn/TsWeFFEnAGMns8wU09jPl+/RfXlw9N5Ye+TaPDvhMy8OTPfw11fcvywiHjYRHUAIuIo7jo1+3UTlZWk+WCCJUmtYS5jg1rBtdz13UWvaSiO6pTop0VE24QlJ1btxXoVd30v1ucyc88U9efrt7hnZX06U9n/dk3fO1dbx2zfd7xCFY/kyDT3CVxTd0CSNBkTLElqDTdV1pt8hG5elO9r+vPKrpdHxC81EMpVHHkJ7grg0bNo4x85MpveE4HzK8cmevfVnebxt6i+i2z9ZAUj4inAb9bwnRO1P5Ne2OPGbN8ybqkjzqisfzsz65pVUpKmxQRLklpD9f+FPzsiTmsskvnzN8COcv0o4NMR8dypKkXEaRFxaUS8bK4BZOYA8J+VXWfOoo3bgI9Wdo1ObvFjYNs0m5mP3+LLlfXOiR49jIgnAJ9mfv9GeHFEfKx8fHJC5eN+f1bZ9WOO/C4TObOy/oXZhSdJs+csgpLUGrYD+4DVFNNs/3f5YtubKMb1AJCZz20kuhpk5qGIeB7wNYp3OB0PfKp8F9InKR6dO0hx/vejeKnvBo70MtX1KNi/cGRCiGdy9/dZTcclwO+M2fehzPzZdCrP029xKcWkESspxod9LSI+APx72dYDgY3AuRSz930BOA148HRinqGjKHrIfjMivld+19XAjRTTxx8HPBZ4KfDzlXp/OtkEIWXP2FmVXf9Sc9ySNCUTLElqAZl5R0S8Avhnih6RNqCj2ajql5nXRcSTgH8FHlHufmK5LJR/onj3UjvwyxFxXGYenEkDmfmtiLiKIy/NTeDvZ9hGrb9FZv44Il4J/ANF79QxwBvKZaxrgRdz1zFp8+XhTO+x178CLp6izBOBB5TrOzJzOpN5SFKtfERQklpEZn6WYja6P6fooThApfdqqcjM6ynO841MPcX2bcBnKHpDPlzT9x+gSLKgmCzhBbNs6r8r6/+emT+YRSy1/haZ+TGKXqrvT9LGXwJPyMybZxrvDHwJeA9FIjfVhBs9wHMys7McnzaZ36qsXzJhKUmaRzH1/1ZJktSciHgk8HiKxyOPpZhR8SbgeuDazDw8D9/5OI48Zjfjl/1GxAqKx91GX3h7bmZ+vIa4avktytkRn0jxaOHxFBNH/BD4cmbePlndukXEicAvAKdQvFD4aIpE74fA1Zk5rfdYRUQ7xRitE8v6D3KCC0lNMMGSJGkcEfE5ijFYAB2ZOe3H5SLiJcBl5eZ+4IGZeUfNIaqifIR2dJbGP8vMtzYZj6Tly0cEJUka3x9U1t883UoRcQ/gwsquvzO5ml/l5Bajv3k/8N4Gw5G0zJlgSZI0jnKChNHH+n4jIh41VZ3y0bs/ppgBD2CAYsIMza8XAT9Xrv/JPI8fk6RJ+YigJEkTiIgHU4xvOgb49HjT4EfEuRRTm7dTTGv+kMrhd2Xm7y9AqMtW+a6s71LMRPh94DR7DCU1yWnaJUmaQGb+L8W7pibzc8BzxtnfDbyr9qB0F+W7xR4xZUFJWiA+IihJUn2GKKYefztwVmYONRyPJGmB+YigJEmSJNXEHixJkiRJqokJliRJkiTVxARLkiRJkmpigiVJkiRJNTHBkiRJkqSa/P8ygbcvB6rkiwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/conda/lib/python3.8/site-packages/lifelines/plotting.py:710: UserWarning: For less visual clutter, you may want to subsample to less than 25 individuals.\n",
      "  warnings.warn(\"For less visual clutter, you may want to subsample to less than 25 individuals.\")\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*******Cluster 14 Z:3.6533003349425597*******\n",
      "0.6 0.37818965517241376 0.0607150588485832\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAALICAYAAABijlFfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABhjUlEQVR4nO3deZxdZX348c/XAAmBsASMKAoBFRVQagjusdgUioLa1qpptbhD1difCy641KlbqUDVSluDVQHBxtZaXCIiItS4koBYEFnUJIgIESIQkswQJt/fH+cMORlmuzPnzpl75/N+vc7rnuV5nvs993BgvjzPeU5kJpIkSZKkiXtI0wFIkiRJUrcwwZIkSZKkmphgSZIkSVJNTLAkSZIkqSYmWJIkSZJUExMsSZIkSaqJCZYkSZIk1cQES5IkSZJqYoIlSVNMRPx5RHwnIu6KiL6IuDEiPhQR+5bHz4mI1TV/50si4pV1tjnCdz0mIpZFxE8joj8iLh+l/P4RcW9EZETsPhkx1m0812y4a9KO6z9Rrf7zExGHRERPROzShlj+IiJ+EBF3RkRvRNwQEe9tx3dJ0lBMsCRpComIM4H/An4F/DVwLPAx4PnAp9v41S8BXtnG9qsOA54H3FguozkduLetEU1Nw12TDw6zv0mt/vPzXOANmXlfG2LZB7gMeG35PZ8F3gP8Uxu+S5IeZKemA5AkFSLi+cBbgddk5mcrh/43Is6mSLY6QkTMAGYM8wf01zLzK2W5LwH7jtDOIuA44CMUida0l5m/bDqGGjwJuLYdDWfmskG7LouIPYA3RsSbMjPb8b2SNMAeLEmaOt4CXDUouQIgM/sz86KhKkXE5WWiUt13dDmk7vDKvsMi4psRsSEiNkXEzyPijRFxDvAi4A/LOhkRPZV6z4qI/42IzeWwq09HxJxB33dORKyOiD+NiJ8BvcBTh4o3M7eN5ccok7RPAh8A7hhjnWoc15dDxL4XEYcOUfYlEXFNOQzz1xHx4YjYqdW2xvr7Dzr+9Ij4akTcWl6LqyPiZdXvZphrMtQQwdHOZdD5HBMR/1d+7/ci4rBRftNxxzpMe78DXg08p1L+j0aKoQZ3Ag4RlDQp7MGSpCkgInYGngGc2cav+SpwPfByoA94HLAHxZCzA4C9gDeUZW8p43omcClwIfAXFMOvTgP2Lrer5gMfpUiIbgfWTDDevwFmAf8CvGyUslUHUgwHex+wBfh74OKIeGxm9gJExLHAF4HzgLdT9Kh8kOL8/qaVtsbpQOD7wKcoktFnAp+LiG2Z+R+McE0Ga+FcKNs8HfhweT5nAP8ZEYeP0LNTZ6xBMdz1Mop/1r9eHrp6iHIzhonnAZl5/3DHygR9JrAA+Fvg3+y9kjQZTLAkaWrYh+KPwZvb0XgUE2QcDPxpZl5T7r60cnwD8JDM/NGgqqcBP8jMl1bK/ga4tPyjvDrMax/gjzPz6hri3YfiD/eXZ+bW4u/tMdsXeGFm/qBs60rglxTPCH2qLPMB4PLMfEW5/c3yO/4hIj6Umbe00FbLMnP5wHqZTHwXeCTwOuA/MvOXI1yTwcZ6LgBzgWdm5k3ldz8E+B+KZPv6dseamRkRt1Ekzl8fofwrgM+N1NZASCMc20RxT8H25FOS2s4hgpI0tbTr/7BvAH4NfCoiXhoR80arEBGzgadT9HDsNLAA3wO2AkcOqvKbOpKr0oeBH2fmN8ZRd/1AQgSQmeuAK4GnwAM9GwsoJhOp+iLFfxefPta2xisi9o6If46IdRS/5VbgJOCQFttp5VwA1g4kV6Xrys9HtjvWiidR/HM+0jNYXwOOGsMykmcAi4C3AS8EzhpnvJLUEnuwJGlquJNi2N4B7Wg8M7eVQ8k+TDGr2q4R8X3gbzPzJ8NU25timNa/lstgjxq0fXsdsZbPBL0aeHZE7FXunl1+7hkR/Zm5ZYQm1g+z7+Hl+r7Azjw43oHtuS20NV7nAE+j6KW7DrgHeD1FItCKVs4F4K5B2wOTkMwa4TvOoZ5YBzwJ+FVmjjQz5Abg7nG2D0BmXlWufi8i7gDOjYgzu2SSEElTmAmWJE0B5TC47wN/Ary3xeq9PPgB/sF/WJOZ1wMvKp/3WgT8I7AiIobrvbiLoqehBxiqJ+nWwV8x9pBH9FiKpOGHQxy7BfgMxRTcwxmqd24e8LNy/Q6KXpjB5R5Wfm5ooS0Y4+8/ICJmAccDSzPzU5X94xlV0sq5tKzmWAc8CfjpKGXqGCJYNZBsHUQxxFOS2sYES5Kmjo8DX42IV2TmudUD5R+0x2bmN4eodwvw7EH7jhnuSzJzK/CdiPgn4AsUkxPcx6BejMzcFBE/Ah6XmR9o8Vwm4nvAcwbtOw54J8X7s341Sv15EfGMynNTB1AMo/scFDMyls9SvRj4t0q9lwDb2DGxG7GtUku/P8VzQTMoeiwp250DvIAdk9QHXZPBWjyX8agt1orHUwwBHMnAEMG6PLP8nOjEK5I0KhMsSZoiMvNrZdLzmXL2vq9QvGD38RSzwa0Fhkqw/gd4TUR8DFhBkZz8SbVARDyJYsa4L1IkKHtTJCw/zcwNEXE98MKI+FOKhOHWzLwVeAfFhBbbgC8BGymGMR4PvCczx/Ki4B2Uz3Y9r9zcH9gjIgZmJPxGZt4BXD6ozvxydeUoQ8ug6NX5fEQMzPz3AYphfedUyryfYjbAzwHLgSdSDIH79KBJIcbS1qi/f1Vm3h0Rq4C/i4h7KBKhd1EMidujUnS4azLYWM+lZW2IFYohhs+OiGdTJGY/Hjy7X2beSTFstmUR8U3g2xS9jP0UydXbgC9WhwdGxNEUsxk+JzMvH893SdKQMtPFxcXFZQotFO8Uuozij9j7gBspkqP9yuPnAKsH1TmVYhKLjcD5bO9hOLw8Pg/4PEVy1QvcBvwHcEB5fF+KRGFDWa+n0vZTKRK7eyhmZruOYuryPStlHhTTCOc3v/yOoZb5w9R5ZXl891HaPgdYDfx5+bv1UUwxfvgQZV8KXFP+xrdQPJ+20zjbGu333+H3AR4DfKf8PW+mSGR7gDsqZYa8JsNc/xHPZYR6A9fihBF+03HHOkx7z6JIfrYCt7Xh/vkgxQQa91IMc70KeBOw86ByzytjPbTpe97FxaW7lsj0lRCSpO5QvvT28MxcOJXa0tQTEX8PPDszBw9HlaQJcZp2SZI0HT2DoidWkmrlM1iSJGnaycyRJiKRpHFziKAkSZIk1cQhgpIkSZJUExMsSZIkSaqJCZYkSZIk1cQES5IkSZJqYoIlSZpSIuKVEZERsXsb2s6IWDqGcq+LiJsiojciroyIxWNsf3ZE/GNE3FzW/VVEvGOMdb8WEe8v1wd+g5GWo8fS7jDf9amIuD4i7o2I30fEdyPij+s6x2Hi/dGgMisi4n3jPQdJmqqcpl2SNNWsAJ4ObG7iyyNiCfApoAf4HvAq4OsRcVRmXjtCvRnAN4D9gPcAvwYeDewzhu98KvAc4BXlroHfYLCdgOVAAP83tjMa0q7AWcANwC7Aa4CLImJRZv5ouEotnuOZwJcq2xsHHT8N+GpEfDIz7xrneUjSlOM07ZKktouIAGZmZm/DcSTwpsw8a4QyNwDfz8xXl9sPAX4K/DQzXz5Cvb8B/gF4XGaubzGuCwAy82WjlPsI8E7gmMz8TivfMUq7M4A1wIWZ+bcjlBvTOY7ldy7L/QL4RGZ+cnyRS9LU4xBBSdIOIuIF5bC4TeXwsR9HxB+Wx+aXw71OGFTnnIhYXdnuiYg7IuJZEbEK6AVeWrb5hiG+c3VEfL5c32GIYESsiYiPDlHnSxGxslzfLSLOiogbImJzWedfImKPFs/9YOAQ4D8H9mXmNuC/gOeOUv3VwH+OI7maA/wZO/b2DFXu2RTJ1el1JlcAmdkP3EXRmzWScZ3jCP4bOLGmtiRpSjDBkiQ9ICIeTfGH/neA5wMvA74OzB1Hc7OBc4F/B44DflS29dJB33kwcCTwxWHa+U/gJWUv2ECd3YHnVerMBmZQDFt7LvA+4I8oEqNWPL78vH7Q/p8DcyPioUNViohdgCcDt0TEBRGxJSLujojPjSHJewbFkL0fDFcgIvYCPg9cRXFu1WMRETuNtgzR5kC9fSLiLcBjgc+OEEOr59gTEfeXifZnI2Kof4Z+ABwZEXsP972S1Gl8BkuSVPVkYGNmvr2y7xvjbGtX4K2Z+ZWBHRGxHPhSRDwiM28td78U+D3wrWHaWQ68A3gqRZIGRfI3kzKByszfAa+vfM9OFEPevhcRB2TmzWOMeeAP/bsG7f995fjvhqi3D8V/U98BXAq8ADgAOB3YDXjJCN95JHBHZt4+QplPUSS5f5yZWwcdewXwuRHqDohB2y8F/qNc3wS8NDOvGKF+K+d4LvA1it9qIUVSeEREPKXsLRvw0zKuhcAlYzgHSZryTLAkSVXXAHtGxLnABRTPIm0aZ1sJXDRo30XAvcCLgU+U+14K/E9m3jdkI5k/iYgby3I/qtS5vJqURMRfA2+l6InZrdLEIcBYE6xq7FUxzP4BAyNCfg+8eCAJioitwLkR8ejM/OUwdfcD7hgukIh4JcX5vjozbxqiyNeAo4arP4KLy3r7UvRULo+I52Xm5cOUH/M5ZuYrK/W+GxE/p0jUnw9cWDk2cN77jSN+SZqSHCIoSXpAZt4AvBA4mOIP4jsi4gvDDY0bxe8HJ03lJBdfoRwmGBGPA46g6KUayReBF5fD2vagGHL4QJ2I+DPgPOCHFMnb0yieawKY1UrM5edeg/YPbN81Sr3vD+phGnhW6tARvnMW0DfUgXL45D8D/5WZw/VSbQCuHsOyg8z8fWauzsxvZuZfU/x2Hxghzomc4zcpEusFg/YPnHcr10iSpjQTLEnSDjJzRWYuohgS9hrgj4GBWd4GZgEcPBnCUM/XDNfb80XgaRFxAEWi9Tu2/5E+nOXA/sCzgD+leN7qy5XjLwZ+nJlvyMyLMvPHbE8IWjHw7NXjB+1/PLChHIr4IJm5GVg3xKGBnq9tI3znBh6c0A0Mc7yAIqk7aYT6rwC2jmEZzU8oEushTeQcc/uUxYP/mdir/NwwhvgkqSN05RDBiHgXxTSyAGTm4HHn1bIPoxhPfgLFWPItwM8oxo9/pvIfBUmaVjLzbuAL5QyCA+9kWk/xx/oTBsqVE048naH/+B7KtyiSn5dQJFhfGvRczlCxXBcR15blDwIuycw7K0V25cG9QCNOeT7M9/yqHI74YoohdAPTtL+YBw93HOzrwAsjYpdKz91iisTjmhHq3QA8IiJmZmb1HN4PPAV4zijviRrvEMEHlBOIPJ3iubWRjOscI+I4YHfgykGH5pefN7YasyRNVV2XYJXDTd4/xrJHUvwHdOAFifcCcyj+D+mzKIajvGDQf/AkqWtFxMkUf2h/E7iV4nmmF1MMvyMzt0XEV4C3RMQ6it6Vt1H8z6kxycytEfE/FM9LPRx40LTtw/gi8P+APYHXDTp2CfAvEfEe4McUMwwuHmtMg/QA50fEWuD7FD1EjwX+aqBAmXReCizOzP8td58OvBz474j4V+BRwD8Cnx1lko3vAzsDTwRWl+0/AziVYjjlfRHxtCHq/TIzf1cmmncOcXxIEbGI4pp9meLZtH3Kc3waxTNSEzrHiDiJYtKKb1M8Y7UAeC9wBcULlKsWAndT/I9NSeoKXTVEsPy/jJ+hGMv9w1HK7knxf+L2oRgSclRmzqF4MHopxf+hPRb4WDtjlqQp5v+AhwL/RNHT9F7g0xTvXxqwlCIp+FfgXyhmomv1vUzLKZKrW4GVLdTZl6K35MJBx5YBZ1IkYF8GDqSSELUiM/8D+BvglRSJ5pOAEzLz2kqxoBimGJV66yiGU+5VxvBBitEQbxrl+24ErmXH92wdW7b/ZxT/PRtqOX485wf8Grgf+AjF/2T8BMXfA8/KzOqMkeM9x18ChwH/Vrb//ygS9GOH6Kk8jmKCk5GGUEpSR4luGgEXEf8P+DjFmPVfUPZkDTVEMCI+SPGHwxbgsMxcM+j4qRT/8ekHDi3/AyhJUu3K91C9JjMPbzqWyVL+j87bKaae/17T8UhSXbqmBysiDgI+TDFM4i1jqDLw5vjlg5Or0icphgzOYBzj+CVJasHZwEMj4o+bDmQSvR74kcmVpG7TNQkWxRCW3SheajnkLE8Dyue0Dig3h3xoOTPvZfuwlWPrClKSpMHKd429gh3f39Xt7gb+tukgJKluXTHJRUS8juJh5m9n5nljqFIdgnHtsKW2j4kf6d0ekiRNWGZ+s+kYJlNm/lvTMUhSO3R8D1ZE7E8xq9EW4OQxVntEZf03I5QbOLZHOQ2xJEmSJA2rG3qwllFM2fvOzPzVGOvMqaxvHqFc9dgcimeyHqSckvYkgN122+3Ixz9+8PspJUmSJHWTK6+88o7MfOjg/R2dYEXEyymmqb2aYkrhRmTm2RQPKLNw4cJcvXp1U6FIkiRJmgTl+yAfpGOHCEbEPIop2fuB12Xm/S1U31hZnz1CueqxjcOWkiRJkiQ6uwfrHyleEvxvwPVDPCO1y8BK5dh9mXkfxYstB+wP3DPMd+xfft5TziooSZIkScPq2B4s4KDy8/UUvUuDl1MrZQf2fbTcrs4cONJLHQeOXTfRYCVJkiR1v05OsMYtM28Abi43jxuqTETsBiwqN781GXFJkiRJ6mwdm2Bl5tGZGcMtwN9Xyg7sf3OliYH3ZS2JiPlDfMUbgd0pnvG6oE2nIUmSJKmLdGyCVYMzgNsoJrJYERFHAkTELhHxeuCDZbmzM/PGhmKUJEmS1EE6eZKLCcnMuyPiBOBi4FBgdURsBGYBO5fFvgW8paEQJUmSJHWY6dyDRWZeCRwGfAy4iSKx2gR8D3gd8NzM7GsuQkmSJEmdJDKz6Ri6ii8aliRJkrpfRFyZmQsH75/WPViSJEmSVCcTLEmSJEmqybSd5KJd1m3YzMnLr2o6DEmSJKkWy5YsaDqEjmIPliRJkiTVxB6smh04d7ZZviRJkjRNmWBJkiRJAmD+u1aMqdza045vcySdyyGCkiRJklQTEyxJkiRJqokJliRJkiTVxARLkiRJkmpigiVJkiRJNTHBkiRJkqSamGBJkiRJUk1MsCRJkiSpJiZYkiRJklQTEyxJkiRJqokJliRJkiTVZKemA+g26zZs5uTlVzUdhiRJkqQG2IMlSZIkSTWxB6tmB86dzbIlC5oOQ5IkSWrZ/KtXNB1Cx7MHS5IkSZJqYoIlSZIkSTUxwZIkSZKkmphgSZIkSVJNTLAkSZIkqSYmWJIkSZJUExMsSZIkSaqJCZYkSZIk1cQES5IkSZJqYoIlSZIkSTUxwZIkSZKkmphgSZIkSVJNdmo6gG6zbsNmTl5+VdNhSJIkSWqAPViSJEmSVBN7sGp24NzZLFuyoOkwJEmSpJbNv3pF0yF0PHuwJEmSJKkmJliSJEmSVBMTLEmSJEmqiQmWJEmSJNXEBEuSJEmSamKCJUmSJEk1McGSJEmSpJr4HixJkiRJAKw97fimQ+h49mBJkiRJUk1MsCRJkiSpJiZYkiRJklQTn8Gq2boNmzl5+VVNhyFJkiTVYtmSBU2H0FHswZIkSZKkmtiDVbMD5842y5ckSZKmKXuwJEmSJKkmJliSJEmSVBMTLEmSJEmqiQmWJEmSJNXEBEuSJEmSatLxCVZELIiI90fEVyPi+oi4MyK2lp/fj4j3RMTcYer2RESOYXnMZJ+XJEmSpM7TDdO0vxp4Y2W7F9gCzAWeUS5vjogXZOYPh2ljK7BhhO+4v45AJUmSJHW3bkiwrgDWAt8Drs/MuwAiYnfgRcDpwEOBCyPikMy8e4g2fpCZR09KtJIkSZK6VscnWJl53jD77wXOjYjfAhcD84ATgAsmMTxJkiRJ00jHJ1hj8KPK+iMbi0KSJEma4ua/a8WYyq097fg2R9K5On6SizFYVFn/ZWNRSJIkSep6XZlgRcTMiJgfEUuBz5e7fwF8bZgqh0XEtRGxJSLujYgbIuLTEfHkyYlYkiRJUjfoqiGCEdELzBzi0PeBv8rMvmGq7ksx6+BdwB7AIeXymoj4SGa+d6wxrNuwmZOXX9VS3JIkSZK6Q7f1YN0G3A5squy7DHhzZt48RPmbgHcAjwNmZeY+wG7AnwBXAgG8JyLeNtKXRsRJEbE6IlZv3bq1htOQJEmS1IkiM5uOoS0iYh7w18B7gL2AD2Xm37VQfxbwXeAo4F7gkcNM8b6DhQsX5urVq8cVsyRJktQkJ7kYu4i4MjMXDt7fbT1YD8jM9Zl5JnAckMD7IuKEFur3Au8uN3cHFtcfpSRJkqRu0rUJ1oDMvILiJcQAJ7VY/YeV9YPriUiSJElSt+r6BKv0m/LzMY1GIUmSJKmrTZcEa6D3aWOL9Z5WWV9TUyySJEmSulRHJ1gRMSMiYpQyi4GnlJuXV/aPVm8m8OFycxNw6fgjlSRJkjQddHSCBTwK+ElEnBwRB1eTpoh4VES8C/gKxXTrG4CPVeo+OyK+HREvj4hHVurtXCZlK4Gnlrs/kJl3tftkJEmSJHW2bnjR8BHAp8r1+yLiHmBXivdZDVgDvCgzb6vsC4qZARcDRMQWip6qPYGdyzLbgNMy86PtC1+SJElSt+j0BOtW4CXA0RS9TQ8H9gX6gZuBn1L0YH0hM7cMqnsNcArwdOCJZb29gM3AdRQ9WGdn5jXtPglJkiRJ3aGjE6zMvA/4r3Jpte6dwJm1ByVJkiRp2ur0Z7AkSZIkacowwZIkSZKkmnT0EMGpaN2GzZy8/Kqmw5AkSZLUAHuwJEmSJKkm9mDV7MC5s1m2ZEHTYUiSJEktm3/1iqZD6Hj2YEmSJElSTUywJEmSJKkmJliSJEmSVBMTLEmSJEmqiQmWJEmSJNXEBEuSJEmSamKCJUmSJEk1McGSJEmSpJqYYEmSJElSTUywJEmSJKkmJliSJEmSVBMTLEmSJEmqyU5NB9Bt1m3YzMnLr2o6DEmSJEkNsAdLkiRJkmpiD1bNDpw7m2VLFjQdhiRJktSy+VevaDqEjmcPliRJkiTVxARLkiRJkmriEEFJkiRJAKw97fimQ+h49mBJkiRJUk1MsCRJkiSpJiZYkiRJklQTEyxJkiRJqokJliRJkiTVxARLkiRJkmpigiVJkiRJNTHBkiRJkqSa+KLhmq3bsJmTl1/VdBiSJElSLZYtWdB0CB3FHixJkiRJqok9WDU7cO5ss3xJkiRpmrIHS5IkSZJqYoIlSZIkSTUxwZIkSZKkmphgSZIkSVJNnORCkiRJEgDz37ViTOXWnnZ8myPpXPZgSZIkSVJNTLAkSZIkqSYmWJIkSZJUExMsSZIkSaqJCZYkSZIk1cQES5IkSZJqYoIlSZIkSTXxPVg1W7dhMycvv6rpMCRJkiQ1wB4sSZIkSaqJPVg1O3DubJYtWdB0GJIkSVLL5l+9oukQOp49WJIkSZJUExMsSZIkSaqJCZYkSZIk1cQES5IkSZJqYoIlSZIkSTXp+AQrIhZExPsj4qsRcX1E3BkRW8vP70fEeyJi7ihtPCwizoyIGyJiS0RsiIiVEfHaiIjJOhdJkiRJna0bpml/NfDGynYvsAWYCzyjXN4cES/IzB8OrhwRRwIXA/uUu+4F5gDPKpcXl3X72ncKkiRJkrpBx/dgAVcAbweeDuydmbtm5h4USdIrgd8B+wIXRsSe1Yrl9tcpkqvrgaMycw6wG7AU2AocC3xsck5FkiRJUifr+AQrM8/LzDMy80eZeVdl/72ZeS7w8nLXPOCEQdVPAfaj6PF6XmauLuvel5n/Ary/LHdSRBzSzvOQJEmS1Pk6PsEagx9V1h856NiJ5efyzFwzRN1PUgwZnAG8rA2xSZIkSeoi0yHBWlRZ/+XASkQ8Djig3LxoqIqZeS+wstw8ti3RSZIkSeoa3TDJxYNExEzg4RRDAj9Q7v4F8LVKscMr69eO0Ny1wHOBQ8fy3es2bObk5VeNPVhJkiRJXaOrEqyI6AVmDnHo+8BfDZoJ8BGV9d+M0OzAsT0iYveyV2vw954EnASw5yMf01rQkiRJkrpGVyVYwG3ALGB3ipkAAS4D3pGZNw8qO6eyvnmENqvH5lA8k7WDzDwbOBtg4cKFuWzJghbDliRJkpo3/+oVTYfQ8brqGazMnJ+Z+2Xm7sDDKGYJ/APgioj4wIiVJUmSJGmCuirBqsrM9Zl5JnAckMD7IqI6TfvGyvrsEZqqHts4bClJkiRJ017XJlgDMvMK4Hvl5kmVQ7dW1vcfoYmBY/cM9fyVJEmSJA3o+gSrNDBRRXUGiurMgdUZBQcbOHZdrRFJkiRJ6jrTJcE6uPx8YIhfZt4ADEx8cdxQlSJiN7a/R+tbbYtOkiRJUlfo6AQrImZERIxSZjHwlHLz8kGHzys/l0TE/CGqv5FiRsJ+4ILxRypJkiRpOujoBAt4FPCTiDg5Ig6uJlsR8aiIeBfwFSCADcDHBtU/g2Jq99nAiog4sqy7S0S8HvhgWe7szLyxzeciSZIkqcN1w3uwjgA+Va7fFxH3ALuy/T1YAGuAF2XmbdWKmXl3ObPgxcChwOqI2EjxLq2dy2LfAt7SxvglSZIkdYlO78G6FXgJ8K/AlcAdwB4U53Uz8DXgtcBhmfmToRrIzCuBwyh6t26iSKw2Ucw8+DrguZnZ197TkCRJktQNOroHKzPvA/6rXCbSzu3AW8tFkiRJksal03uwJEmSJGnKMMGSJEmSpJp09BDBqWjdhs2cvPyqpsOQJEmSWvYnf/DwB+1btmRBA5F0LnuwJEmSJKkm9mDV7MC5s83yJUmSpGnKHixJkiRJqokJliRJkiTVxARLkiRJkmpigiVJkiRJNTHBkiRJkqSamGBJkiRJUk1MsCRJkiSpJiZYkiRJklQTEyxJkiRJqokJliRJkiTVxARLkiRJkmpigiVJkiRJNdmp6QC6zboNmzl5+VVNhyFJkiTVYtmSBU2H0FFMsCRJkiQBcPHVv33QvvlXr3jQvrWnHT8Z4XQkE6yaHTh3tlm+JEmSOtJQyZRa4zNYkiRJklQTEyxJkiRJqokJliRJkiTVxARLkiRJkmpigiVJkiRJNTHBkiRJkqSamGBJkiRJUk1MsCRJkiSpJiZYkiRJklQTEyxJkiRJqokJliRJkiTVZKemA+g26zZs5uTlVzUdhiRJkqQG2IMlSZIkSTWxB6tmB86dzbIlC5oOQ5IkSWrZ/KtXNB1Cx7MHS5IkSZJqYoIlSZIkSTUxwZIkSZKkmphgSZIkSVJNTLAkSZIkqSYmWJIkSZJUExMsSZIkSaqJCZYkSZIk1cQES5IkSZJqYoIlSZIkSTUxwZIkSZKkmphgSZIkSVJNdmo6gG6zbsNmTl5+VdNhSJIkSWqAPViSJEmSVBN7sGp24NzZLFuyoOkwJEmSpJbNv3pF0yF0PHuwJEmSJKkmJliSJEmSVBMTLEmSJEmqiQmWJEmSJNXEBEuSJEmSatLxCVZE7BMRr4qI8yPiuojYFBF9EXFLRFwYEX82Qt2eiMgxLI+ZzHOSJEmS1Jm6YZr229jxPHqBrcD+5fLCiLgI+IvM3DxMG1uBDSN8x/11BCpJkiRNZWtPO77pEDpex/dgUSRXVwBvAB6dmbtm5u7AQcBnyjLPBZaN0MYPMnO/EZa1bT0DSZIkSV2hG3qw/igzLxu8s0yKXhsR9wMnAy+PiHdn5q8nO0BJkiRJ00PH92ANlVwN8pnK+sJ2xiJJkiRpeuv4BGsMeivrMxqLQpIkSVLX64YhgqM5urJ+zTBlDouIa4FHA/3Ab4DvAv+amT9p5cvWbdjMycuvGk+ckiRJ0pSzbMmCpkPoKF3dgxURewGnlpsrM/OGYYruCzwB2AzMBA4BXgtcGREfGsP3nBQRqyNi9datWyceuCRJkqSOFJnZdAxtEREPAb4CnAD0AU/NzJ8OKvMy4BFluTWZuTUidqHo9foIcGRZ9JTMPHMs37tw4cJcvXp1PSchSZIkaUqKiCsz80FzPHTzEMFPUCRXAG8YnFwBZOYFQ+y7D/hWRHyXYpjgUUBPRPx7Zt7dzoAlSZIkdbauHCIYEWcAS8vNt2TmZ1ttIzN7gXeXm7sDi2sKT5IkSVKX6roEKyI+Cryt3Hx7Zn58As39sLJ+8ATakSRJkjQNdNUQwYg4HTil3HxHZp7RZDySJEmSppeuSbDKYYEDPVfvyMzTa2j2aZX1NTW0J0mSJKmLdcUQwUHJ1SljSa4iIkY5PhP4cLm5Cbh0QkFKkiRJ6nodn2BFxD+yPbl661inUweeHRHfjoiXR8QjK+3tHBGLgZXAU8vdH8jMu2oLWpIkSVJX6ughghFxAPCOcnMb8M6IeOcIVc6oPJcVFDMDLi7b2kLRU7UnsHOlzdMy86N1xy5JkiSp+3R0gsWOPXAPAR42SvndK+vXUEyI8XTgicC+wF7AZuA6ih6sszPzmrqClSRJkqay+e9aMaZya087vs2RdK6OTrAycy1FT9R46t4JjHU4oSRJkiSNquOfwZIkSZKkqWJSE6yIOCAizo6IdRHRGxG3RsQXIuKIyYxDkiRJktphwkMEI2IucDNFsnYz8ITMzCHKPRG4DNib7cP69gNeCrwoIl6RmcsnGk/T1m3YzMnLr2o6DEmSJEkNqKMH64+B2cBM4JxhkqsAvgDMHaJ+UMzad25EHFpDPJIkSZLUiDomuTi2sv4/w5T5c+AwYCD5Ohe4ENgDeBfwhDKWD5VlO9aBc2ezbMmCpsOQJEmSWjb/6rHNIqjh1ZFg/UH5eXdm3jBMmb+urJ+dma8f2IiIb1JMi74PcHxE7JmZd9cQlyRJkiRNqjqGCB5A0TP1i6EORsROlC/zLe0wNXpm/g74TLm5E3BUDTFJkiRJ0qSrI8GaU34O1+u0ENiNIgm7LjOHSsR+VFl/bA0xSZIkSdKkqyPBGmhjxjDHn1VZv2yYMr+rrO854YgkSZIkqQF1JFgDPVePGOb4H1XWvztMmV0q6w+ahVCSJEmSOkEdCdaNFFOtPyYi9qseiIg9gedUdl0+TBvzKutOcCFJkiSpI9WRYK0sPwP4+0HH3k3xfqwErsrMO4Zp44jK+toaYpIkSZKkSVfHNO2fA06hSNZeGxFPAL5HkTQdVyn37yO08YeV9WtriEmSJEmSJt2EE6zMvDEiTqd4YXACzyyXqmvZPhX7DiLiUcDTyrq3ZOYtE41JkiRJkppQxxBBMvPdwIeB+ymGClaX1cAJmXn/MNXfUJYDuLiOeCRJkiSpCXUMEQQgM98XEZ8EjgUeCfQCqzPze6NU3RU4t1wfspdLkiRJkjpBbQkWQGauB85vsc6b64yhaes2bObk5Vc1HYYkSZKkBtQyRFCSJEmSVEMPVkR8tly9JjM/Ns42lgILgMzM10w0piYdOHc2y5YsaDoMSZIkqWXzr17RdAgdr44erFcCrwCOmUAbi8t2XjnxcCRJkiSpGQ4RlCRJkqSaTLUEK5sOQJIkSZLGa6okWHuUn5sbjUKSJEmSJqDxBCsidgH+gKL36vZmo5EkSZKk8WtpFsGIOGCEw7uOcnyHpiheMPwY4PXA3hQJli+QkiRJktSxWp2mfS1DPycVwLOBNROM5wsTrC9JkiRJjRnve7BijPtacX5mfmWCbUiSJElSY8aTYE00kYKiF2wTcBvFsMAvZOZXa2hXkiRJkhrTUoKVmQ+aFCMitlEkTBdn5vPqCkySJEmSOk1dswjW0aslSZIkSR1tvM9gVb2q/PxNDW1JkiRJUseacIKVmefWEUi3WLdhMycvd7Z5SZIkaTpq/EXDkiRJktQt6hgiqIoD585m2ZIFTYchSZIktWz+1SuaDqHj1Z5gRcSBwDOAJwB7A7MZ+yQYmZmvqTsmSZIkSZoMtSVYEXEE8E/A0RNsygRLkiRJasDa045vOoSOV0uCFREvAr5QtjeRKduzjngkSZIkqQkTTrDKIYGfB3Zme4J0M/B94BZg80S/Q5IkSZI6QR09WG8BZlEkV/cCr8nML9XQriRJkiR1lDoSrGMq6y/LzK/X0KYkSZIkdZw63oP1KIreq1tMriRJkiRNZ3UkWAPPXf2ihrYkSZIkqWPVkWD9uvycXUNbkiRJktSx6kiwvk0xNfthEbFzDe1JkiRJUkeqI8FaBtwP7Ab8dQ3tSZIkSVJHmnCClZk/B95D0Yt1ZkQ8ecJRSZIkSVIHqmOadjLz9IjYBfgA8L2I+Afg7MxcX0f7nWTdhs2cvPyqpsOQJEmSarFsyYKmQ+goE06wIuI7lc3fA3OBvwf+PiLWALcB942xuczMxRONSZIkSZKaUEcP1tFsn6qdynoAB5XLWMSgdjrSgXNnm+VLkiRJ01QtQwQpkqNW9kuSJElS16kjwXpVDW1IkiRJUsebcIKVmefWEYgkSZKkZs1/14oxlVt72vFtjqRz1fEeLEmSJEkSXZBgRcQ+EfGqiDg/Iq6LiE0R0RcRt0TEhRHxZ2No42ERcWZE3BARWyJiQ0SsjIjXRoTPkUmSJEkak7omuWjSbex4Hr3AVmD/cnlhRFwE/EVmbh5cOSKOBC4G9il33QvMAZ5VLi+OiBdkZl/7TkGSJElSN2hLD1ZE7B4Rx0bEu8qeoX+PiM+047sokqsrgDcAj87MXTNzd4rp4Qe+87nAsiHi3BP4OkVydT1wVGbOAXYDllIkascCH2tT7JIkSZK6SK09WBGxP/B3wMuBWdVDFO+4es0Qdb4NHFoePzozb2rxa/8oMy8bvDMz1wKvjYj7gZOBl0fEuzPz15VipwD7AVuA52XmmrLufcC/RMQewEeAkyLi45l5Y4uxSZIkSZpGauvBiogTgP8DXgvsSpFUDSwjOY8iydkPeEWr3ztUcjVIteds4aBjJ5afyweSq0E+STFkcAbwslZjkyRJkjS91JJgRcQfA18C9qJIqO4Dvg18HPjlKNX/C9hUrv9pHfEM0ltZnzGwEhGPAw4oNy8aqmJm3gusLDePbUNskiRJkrrIhBOsiNgV+BywS7lrBXBwZh6bmW8FfjFS/czcAnyLIjF7QkQ8bKIxDXJ0Zf2ayvrhlfVrR6g/cOzQugKSJEmS1J3qeAbr1RSz9SVFovSCzMwW2/gxMDCd+hOB22uIi4jYCzi13FyZmTdUDj+isv6bEZoZOLZHROxe9moNa92GzZy8/KqWY5UkSZLU+eoYIvj8yvqbxpFcAVQTn4MnGA8AEfEQ4PPAw4E+4E2DisyprD9o+vZhjs0ZqkBEnBQRqyNi9datW8cTriRJkqQuUEcP1mHl542ZOeJwwBH8vrK+5wTjGfAJ4IRy/Q2Z+dOa2n2QzDwbOBtg4cKFuWzJgnZ9lSRJktQ2869e0XQIHa+OHqyHUgwPHGmYXStx1PFc2BkU77ECeEtmfnaIYhsr67NHaK56bOOwpSRJkiRNe3UkWAPPJI2UpIym+jzUnRNoh4j4KPC2cvPtmfnxYYreWlnff4QmB47dM9rzV5IkSZKmtzoSrNsoZgA8NCJGe+fVcJ5dWV873kAi4nTg7eXmOzLzjBGKV2cOPHzYUtuPXTfeuCRJkiRND3UkWN8rP+ew/ZmnMStn+vvLcvM+4PvjCaIcFnhKufmOzDx9pPLljII3l5vHDdPmbsCicvNb44lLkiRJ0vRRR4L1P5X1MyNiyJn2RvBZiuQsgW+U78VqSZlcDQwLPGW05KrivPJzSUTMH+L4G4HdgX7gglbjkiRJkjS9TDjBysyLgR9QDBN8NHBZRBwyWr2IeGRErABeONAU8MFWvz8i/pHtydVbM/PMFqqfQTHEcTawIiKOLNvcJSJeX4nn7My8sdXYJEmSJE0vdUzTDvA6iqGCewFPBq6JiIuA7wD7DRSKiGOAxwKLgeOBnSkSswQ+nJlXt/KlEXEA8I5ycxvwzoh45whVzqg+l5WZd0fECcDFwKHA6ojYCMwqY4NiaOBbWolLkiRJ0vRUS4KVmT+PiOcDFwL7UiQnz2f7S4iTIpH6ZqVadUKMT2fm+8fx1YOnd3/YKOV3H7wjM6+MiMOAd1I8Q/YoYBPFJBjnAp/NzG3jiE2SJEnSNFNXDxaZ+YOIeBLwb8AL2J5AZeVz8CyD64H3ZOZnxvmda4doczzt3A68tVwkSZIkaVxqS7AAMvM24M/KZ7BeTjED3xMphg4+BLifIqn6EUVv1gXjmdRCkiRJkqaiWhOsAeWEEH83sF2+H2vXzNzcju+TJEmSpKmgLQnWYJmZgMmVJEmSpK42KQnWdLJuw2ZOXn5V02FIkiRJakAdLxqWJEmSJGEPVu0OnDubZUsWNB2GJEmS1LL5V69oOoSON6YEKyK+U9nMzFw8zLGJ2qFtSZIkSeokY+3BOprt77HKYY5N1FBtS5IkSVLHaGWI4Egv9J3wy34lSZIkqdONNcF61TiPSZIkSdK0MaYEKzPPHc8xSZIkSZpOnKZdkiRJkmpigiVJkiRJNZnwe7Ai4u/K1V9k5hfG2cZLgMcDZOYHJhqTJEmSJDWhjhcN91BMr34xMK4EC/hL4IVlOyZYkiRJkjqSQwQlSZIkqSZ19GBJkiRJ6gJrTzu+6RA63lRJsHYrP7c0GkUN1m3YzMnLr2o6DEmSJKkWy5YsaDqEjtL4EMGICODwcvOOJmORJEmSpIlorAcrImYBjwH+FtiPYoKLnzYVT10OnDvbLF+SJEmaplpKsCKif7hDwJ+McHysvjTB+pIkSZLUmFZ7sIKipymGOTYRlwHnT7ANSZIkSWrMeJ7BmmgiVdULXAG8GTguM7PGtiVJkiRpUrXag3XQoO0AfkXRq/Vd4JVjbGcbsAm4KzO3tRiDJEmSJE1JLSVYmblu8L5iEkAC2DLUcUmSJEmaLuqYRfDvy89f1NCWJEmSJHWsCSdYmfn3o5eSJEmSpO7X+IuGJUmSJKlbmGBJkiRJUk3qeAbrQaKY+eIw4JHAnsDMsdbNzPPaEZMkSZIktVutCVZEPAV4K3A8MHscTSRggiVJkiSpI9WWYEXEh4F3UkzZPtrLiHMMZSRJkiSpo9SSYEXEqcCplV29wLXAo4G92f4i4jnA/sDDynJZfl5R1ul46zZs5uTlVzUdhiRJktSyi6/+7ZjKrT3t+DZH0rkmPMlFROwPvL/cTOCfgYdm5lMoEqfiQOZzMnNhZj4ceAzwEWBLeXgm8JrMfM5E45EkSZKkptTRg/U6YBeK5OqCzHzzaBUy81fAeyPic8DFwBHApRHx5My8q4aYGnPg3NksW7Kg6TAkSZKkls2/ekXTIXS8OqZpP7qy/oFWKmbmL4HjgD7gAOATNcQjSZIkSY2oI8F6bPl5c2b+YrhCETFkb1lZ51yKSS9eHBF71BCTJEmSJE26OhKsuRTDA9cNcey+yvquI7TxrfJzJvCHNcQkSZIkSZOujgRrW/m5dYhjGyvrjxihjd9V1h854YgkSZIkqQF1JFh3lp97DnHs1sr6oSO0sU9lfah2JEmSJGnKqyPBuoHi+alHD3Hsp5X1543QRvXYhhpikiRJkqRJV0eCtar83CsiDh507GLg/nL95RHx9MGVI+IY4FWVXT+uISZJkiRJmnR1JFjfrqzv8ErnzLwD+BJFD9dM4LKI+GxEvD4i/iYiLgBWULyPK4GrM7Pa6yVJkiRJHaOOFw3/L8VzWPsArwE+Oej4W4E/Ah5K8ULiV5TLgCg/+4CTa4hHkiRJkhox4QQrM/sj4i8oZ/+LiF0zc0vl+G0RcTTwn8DhbE+oqn4L/FVmrp5oPJIkSZLUlDp6sMjM/x3l+PUR8QfACcCxwAHAzhSJ1eXAf2Zmbx2xSJIkSVJTakmwxiIztwFfLZeutW7DZk5eflXTYUiSJElqQB2TXEiSJEmSaHMPVkTMoBgOuBfFLIL3AOvL2QW70oFzZ7NsyYKmw5AkSZJaNv/qFU2H0PFqT7Ai4qHAK4E/B46gSKwGl/kNsBI4JzMvqTsGSZIkSWpCbUMEI2JGRJwKrAVOA54CzKKYNXDw8khgCfDNiPhRRDyprjgkSZIkqSm1JFgRMQv4GvAhYNeB3aNVK5enAD+KiONHKS9JkiRJU1pdQwQ/BxwHZLm9BfgycBFwLcWLiPuAOcDBFEnVEuBJZZ1ZwH9HxFMy8/9qikmSJEmSJtWEE6yIOBZ4KUWiFMAK4OTMvHWI4ndSDCH8DnBaRCwBlgG7A7sA/wY8c6IxSZIkSVIT6hgi+OrK+teBFw6TXD1IZi6n6PnqL3c9LSIe38qXR8TsiHhuRLw3Ir4cEesiIsulZ5S6PZWyIy2PaSUmSZIkSdNTHUMEn1Z+JvDG8oXCY5aZP4yI89ieqD0NuL6FJp4CfKOV7xzCVmDDCMfvn2D7kiRJkqaBOhKsh1EkVz/LzF+Ps40VbE+wHjaO+r8HrqosHwP2a6H+DzLz6HF8ryRJkiQ9oI4EawNFMjORlwdX6/6+xborM3NudUdEnDaBWCRJkiRpXOp4ButXFJNbPGoCbVTr/rKVipnZP3opSZIkSWq/OhKsL5afj46II8bZxkvKz98Bl084IkmSJElqQB0J1rnAmnL9MxExp5XKEfGXwAsonuP6SEM9UodFxLURsSUi7o2IGyLi0xHx5AZikSRJkmrV39vHqqU9rDjiBFYt7aG/t6/pkLrWhJ/BysyNEfEi4BLgycD3I+JvMvMHI9WLiNnAu4BTy13nZOY/TzSecdoXmAvcBewBHFIur4mIj2Tme8fa0MYb13DpMSe2JUhJkiRpPO65aS29t98B25JN624B4KizepoNqkvV8aLhZ5erb6aYve9wYGVE/BT4JnAtxQuG7wPmAAdRTK1+fLkNsBw4t9LWkDLzuxONd5CbgHcAXwHWZObWiNgFOBr4CHAk8J6I+H1mnjlcIxFxEnASwPxZLXXgSZIkSW239a6NsC0B6N/Sx/qVqxuOqHvVMYvg5RTD+wYkxaQXfwCM9ExWVMq/tFxGktQT7/YGMy8YYt99wLci4rvAd4GjgJ6I+PfMvHuYds4GzgZYuHBhLr7kvDrDlCRJkiZk1dIe1px/If1b+pix60zmLVrYdEhdq86EJdieaGVl33CqSdlI5RqRmb0R8W6KoY+7A4uBLzcblSRJktS6BWcUT+WsX7maeYsWPrCt+tWRYN3MjslSN/lhZf3gxqKQJEmSJmDGrJk+czVJ6pjkYn4NcUiSJElSx6tjmvZu9rTK+pphS0mSJEkS0zjBiogRn/uKiJnAh8vNTcClbQ9KkiRJUkerdVa+pkTE3sCMyq6BxHF2ROxb2d+bmfeW68+OiPcB5wCXZ+YtZVs7A88G/oFiBkGAD2TmXW0KX5IkSZoS1p52fNMhdLyuSLCAnwAHDrH/7eUy4FzgleV6UMwMuBggIrZQ9FTtCexcltkGnJaZH60/ZEmSJEndplsSrPG4BjgFeDrwRGBfYC9gM3AdsBI4OzOvaSpASZIkSZ2lKxKs8cxkmJl3AmfWH40kSZKk6WraTnIhSZIkSXUzwZIkSZKkmphgSZIkaVz6e/tYtbSHFUecwKqlPfT39jUdktS4rngGayrZeOMaLj3mxKbDkCRJart7blpL7+13wLZk07pbADjqrJ5mg5IaZg+WJEmSxmXrXRthWwLQv6WP9StXNxyR1Dx7sGo255CDWHzJeU2HIUmS1Harlvaw5vwL6d/Sx4xdZzJv0cKmQ5IaZ4IlSZKkcVlwxqkArF+5mnmLFj6wLU1nJliSJEkalxmzZvrMlTSIz2BJkiRJUk1MsCRJkiSpJiZYkiRJklQTEyxJkiRJqomTXEiSJEkCYP67Voyp3NrTjm9zJJ3LHixJkiRJqokJliRJkiTVxARLkiRJkmpigiVJkqad/t4+Vi3tYcURJ7BqaQ/9vX1NhySpSzjJRc023riGS485sekwJEnSCO65aS29t98B25JN624B4KizepoNSlJXsAdLkiRNO1vv2gjbEoD+LX2sX7m64YgkdQt7sGo255CDWHzJeU2HIUmSRrBqaQ9rzr+Q/i19zNh1JvMWLWw6JEldwgRLkiRNOwvOOBWA9StXM2/Rwge2JWmiTLAkSdK0M2PWTJ+5ktQWPoMlSZIkSTUxwZIkSZKkmphgSZIkSVJNTLAkSZIkqSYmWJIkSZJUExMsSZIkSaqJCZYkSZIk1cQES5IkSZJqYoIlSZIkSTUxwZIkSV2nv7ePVUt7WHHECaxa2kN/b1/TIUmaJnZqOoBus/HGNVx6zIlNhyFJ0rR2z01r6b39DtiWbFp3CwBHndXTbFCSpgV7sCRJUtfZetdG2JYA9G/pY/3K1Q1HJGm6sAerZnMOOYjFl5zXdBiSJE1rq5b2sOb8C+nf0seMXWcyb9HCpkOSNE2YYEmSpK6z4IxTAVi/cjXzFi18YFuS2s0ES5IkdZ0Zs2b6zJWkRvgMliRJkiTVxARLkiRJkmpigiVJkiRJNTHBkiRJkqSamGBJkiRJUk1MsCRJkiSpJiZYkiRJklQTEyxJkiRJqokJliRJapv+3j5WLe1hxREnsGppD/29fU2HJElttVPTAXSbjTeu4dJjTmw6DEmSpoR7blpL7+13wLZk07pbADjqrJ5mg5KkNrIHS5Iktc3WuzbCtgSgf0sf61eubjgiSWove7BqNueQg1h8yXlNhyFJ0pSwamkPa86/kP4tfczYdSbzFi1sOiRJaisTLEmS1DYLzjgVgPUrVzNv0cIHtiVNTWtPO77pEDqeCZYkSWqbGbNm+syVpGnFZ7AkSZIkqSYmWJIkSZJUExMsSZIkSapJxydYETE7Ip4bEe+NiC9HxLqIyHLpGWMbD4uIMyPihojYEhEbImJlRLw2IqLNpyBJkiSpS3TDJBdPAb4x3soRcSRwMbBPueteYA7wrHJ5cUS8IDN99bwkSZKkEXV8D1bp98ClwOnAXwK3jaVSROwJfJ0iuboeOCoz5wC7AUuBrcCxwMfaELMkSZKkLtMNPVgrM3NudUdEnDbGuqcA+wFbgOdl5hqAzLwP+JeI2AP4CHBSRHw8M2+sMW5JkiRJXabje7Ays38C1U8sP5cPJFeDfJJiyOAM4GUT+B5JkiRJ00DHJ1jjFRGPAw4oNy8aqkxm3gusLDePnYy4JEmSJHWuaZtgAYdX1q8dodzAsUPbGIskSR2rv7ePVUt7WHHECaxa2kN/r/NCSZq+uuEZrPF6RGX9NyOUGzi2R0TsXvZqDWvjjWu49JgTRyoiSVJXueemtfTefgdsSzatuwWAo87qaTYoSWrIdO7BmlNZ3zxCueqxOUMViIiTImJ1RKy+b+vWWoKTJKlTbL1rI2xLAPq39LF+5eqGI5Kk5kznHqzaZObZwNkACxcuzMWXnNdwRJIkTZ5VS3tYc/6F9G/pY8auM5m3aGHTIUlSY6ZzgrWxsj4buGeYcrOHqSNJkoAFZ5wKwPqVq5m3aOED25I0HU3nBOvWyvr+DJ9g7V9+3jPa81eSJE1HM2bN9JkrqUvMf9eKMZVbe9rxbY6kc03nZ7CqMwcePmyp7ceua2MskiRJkrrAtE2wMvMG4OZy87ihykTEbsCicvNbkxGXJEmSpM41bROs0sBsFEsiYv4Qx98I7A70AxdMVlCSJEmSOlNXJFgRsXdE7DuwsP28Zlf3R8Tug6qeAdxGMZHFiog4smxvl4h4PfDBstzZmXnjZJyLJEmSpM7VFQkW8BPgd5XlUeX+tw/af1a1UmbeDZwA3AkcCqyOiHuAe4F/BXahGBr4lvafgiRJkqRO1y0J1rhl5pXAYcDHgJuAnYFNwPeA1wHPzcy+5iKUJEmS1Cm6Ypr2zJw/wfq3A28tF0mSJEkal2nfgyVJkiRJdTHBkiRJLevv7WPV0h5WHHECq5b20N/raHpJgi4ZIjiVbLxxDZcec2LTYUiS1Fb33LSW3tvvgG3JpnW3AHDUWT3NBiVJU4A9WJIkqWVb79oI2xKA/i19rF+5uuGIJGlqsAerZnMOOYjFl5w3ekFJkjrYqqU9rDn/Qvq39DFj15nMW7Sw6ZAkaUowwZIkSS1bcMapAKxfuZp5ixY+sC1J050JliRJatmMWTN95kqShuAzWJIkSZJUExMsSZIkSaqJCZYkSZIk1cQES5IkSZJqYoIlSZIkSTUxwZIkSZKkmphgSZIkSVJNTLAkSZIkqSYmWJIkSZJUExMsSZK6RH9vH6uW9rDiiBNYtbSH/t6+pkOSpGlnp6YD6DYbb1zDpcec2HQYkqRp6J6b1tJ7+x2wLdm07hYAjjqrp9mgJGmasQdLkqQusfWujbAtAejf0sf6lasbjkiSph97sGo255CDWHzJeU2HIUmahlYt7WHN+RfSv6WPGbvOZN6ihU2HJEnTjgmWJEldYsEZpwKwfuVq5i1a+MC2JGnymGBJktQlZsya6TNXktQwn8GSJEmSpJqYYEmSJElSTUywJEmSJKkmJliSJEmSVBMTLEmSJEmqiQmWJEmSJNXEBEuSJEmSauJ7sCRJkiQBsPa045sOoePZgyVJkiRJNTHBkiSpw/X39rFqaQ8rjjiBVUt76O/tazokSZq2HCJYs403ruHSY05sOgxJ0jRyz01r6b39DtiWbFp3CwBHndXTbFCSNE3ZgyVJUofbetdG2JYA9G/pY/3K1Q1HJEnTlz1YNZtzyEEsvuS8psOQJE0jq5b2sOb8C+nf0seMXWcyb9HCpkOSpGnLBEuSpA634IxTAVi/cjXzFi18YFuSNPlMsCRJ6nAzZs30mStJmiJ8BkuSJEmSamKCJUmSJEk1McGSJEmSpJqYYEmSJElSTUywJEmSJKkmJliSJEmSVBMTLEmSJEmqiQmWJEmSJNXEFw1LkjRJ+nv7uOqUf2D9ytXMW7SQBWecyoxZM5sOS5IeMP9dK8ZUbu1px7c5ks5lglWzjTeu4dJjTmw6DEnSFHTPTWvpvf0O2JZsWncLAEed1dNsUJKkWjlEUJKkSbL1ro2wLQHo39LH+pWrG45IklQ3e7BqNueQg1h8yXlNhyFJmoJWLe1hzfkX0r+ljxm7zmTeooVNhyRJqpkJliRJk2TBGacC7PAMliSpu5hgSZI0SWbMmukzV5LU5XwGS5IkSZJqYoIlSZIkSTUxwZIkSZKkmkz7BCsiXhkROYblj5uOVZIkSdLU5iQX220DfjfC8b7JCkSSJElSZzLB2u7XmTm/6SAkSZIkda5pP0RQkiRJkupigiVJkiRJNTHBkiRJkqSamGBt99CIuDIi7o2ILRHxq4g4PyKObjowSVLn6O/tY9XSHlYccQKrlvbQ3+scSZI0nTjJxXazgQXA74HdgIPK5WUR8TngpMy8f7RGNt64hkuPObGtgUqSpq57blpL7+13wLZk07pbADjqrJ5mg5IkTRp7sOBW4O+BI4BZmTmXItl6JvDtssyrgI8N10BEnBQRqyNi9X1bt7Y7XknSFLb1ro2wLQHo39LH+pWrG45IkjSZpn0PVmZ+C/jWoH39wA8i4k+ALwMvBN4QEf+cmTcN0cbZwNkACxcuzMWXnNf+wCVJU9KqpT2sOf9C+rf0MWPXmcxbtLDpkCRJk2jaJ1gjycxtEXEKRYL1EOD5wD81G5UkaSpbcMapAKxfuZp5ixY+sC1Jmh5MsEaRmb+IiDuAfYGDm45HkjS1zZg102euJGka8xksSZIkSaqJCdYoIuLRFL1XAGuajEWSJEnS1DatE6yIiDEcP73c3AZ8ve1BSZIkSepY0zrBAg6MiCsi4uSIOHgg4YqIh0TE04CLgD8ryy7LzBsai1SSJEnSlOckF3BUuQD0RcRGYA4ws1Lmc8DfTnZgkiRJkjrLdE+wbgfeBDwd+APgocDeQC/F81Y/AD6bmd9vKkBJkiRJnWNaJ1iZuQU4q1wkSZIkaUKm+zNYkiRJklQbEyxJ0rj09/axamkPK444gVVLe+jv7Ws6JEmSGjethwi2w8Yb13DpMSc2HYYktd09N62l9/Y7YFuyad0tABx1Vk+zQUmS1DB7sCRJ47L1ro2wLQHo39LH+pWrG45IkqTm2YNVszmHHMTiS85rOgxJartVS3tYc/6F9G/pY8auM5m3aGHTIUmS1DgTLEnSuCw441QA1q9czbxFCx/YliRpOjPBkiSNy4xZM33mSpKkQXwGS5IkSZJqYoIlSZIkSTVxiKAkSZIkANaednzTIXQ8e7AkSZIkqSYmWJIkSZJUExMsSZIkSaqJCZYkSZIk1cQES5IkSZJqYoIlSZIkSTUxwZKkNujv7WPV0h5WHHECq5b20N/b13RIkiRpEvgerJptvHENlx5zYtNhSGrYPTetpff2O2BbsmndLQAcdVZPs0FJkqS2swdLktpg610bYVsC0L+lj/UrVzcckSRJmgz2YNVsziEHsfiS85oOQ1LDVi3tYc35F9K/pY8Zu85k3qKFTYckSZImgQmWJLXBgjNOBWD9ytXMW7TwgW1JktTdTLAkqQ1mzJrpM1eSJE1DPoMlSZIkSTUxwZIkSZKkmphgSZIkSVJNfAZLkiRJEgDz37ViTOXWnnZ8myPpXPZgSZIkSVJNTLAkSZIkqSYmWJIkSZJUExMsSZIkSaqJCZbUsP7ePlYt7WHFESewamkP/b19TYckSZKkcXIWwZptvHENlx5zYtNhqIPcc9Naem+/A7Ylm9bdAsBRZ/U0G5QkSZLGxR4sqWFb79oI2xKA/i19rF+5uuGIJEmSNF72YNVsziEHsfiS85oOQx1k1dIe1px/If1b+pix60zmLVrYdEiSJEkaJxMsqWELzjgVgPUrVzNv0cIHtiVJktR5TLCkhs2YNdNnriRJkrqEz2BJkiRJUk1MsCRJkiSpJiZYkiRJklQTEyxJkiRJqokJliRJkiTVxARLkiRJkmpigiVJkiRJNTHBkiRJkqSamGBJkiRJUk1MsKQJ6u/tY9XSHlYccQKrlvbQ39vXdEiSJElqyE5NB9BtNt64hkuPObHpMDSJ7rlpLb233wHbkk3rbgHgqLN6mg1KkiRJjbAHS5qgrXdthG0JQP+WPtavXN1wRJIkSWqKPVg1m3PIQSy+5Lymw9AkWrW0hzXnX0j/lj5m7DqTeYsWNh2SJEmSGmKCJU3QgjNOBWD9ytXMW7TwgW1JkiRNPyZY0gTNmDXTZ64kSZIE+AyWJEmSJNXGBEuSJEmSamKCJUmSJEk1McEqRcSciOiJiGsi4t6IuDsiVkXE2yJil6bjkyRJkjT1OckFEBEHApcD88tdm4GZwMJyeVlELM7M3zcSoCRJkqSOMO17sCJiBvA1iuTqt8AxmbkbMBtYAmwEngxc0FSMkiRJkjrDtE+wgFcCTyzXX5SZ3wbIzG2Z+UXg5PLYcyNicQPxSZIkSeoQJljwivLzssz84RDHlwNryvUTJyckSZIkSZ1oWidYETEbeGa5edFQZTIzgW+Wm8dORlySJEmSOtO0TrCAJ7D9N7h2hHIDx/aLiLkjNXj3db9g1dIe+nv76ohPkiRJUgeZ7gnWIyrrvxmhXPXYI4YtBfT39rHm/Au56pR/mFBgkiRJkjrPdJ+mfU5lffMI5arH5gw+GBEnAScB7MtO9G/pY/3K1fVEKEmSJE2Stacd33QIHW+692DVIjPPzsyFmblwTuzEjF1nMm/RwqbDkiRJkjTJpnuCtbGyPnuEctVjG4ctBcyYNZODXv6nLDjj1AkFJkmSJKnzTPcE69bK+v4jlKseu3XYUsCehz6Go87qYcasmRMKTJIkSVLnme4J1s+BbeX64SOUGzh2W2ZuaG9IkiRJkjrVtE6wMnMz8P1y87ihykREAH9Sbn5rMuKSJEmS1JmmdYJVOrf8fE5EPHWI4y8GDi7Xz5uckCRJkiR1osjMpmNoVETsBFwFPJHifVevyMxLI+IhwIuAfwf2AC7KzOeNob2NwA1tDFmTa1/gjqaDUG28nt3F69ldvJ7dxevZXbyeQzswMx86eOe0T7AAImI+cBkwv9y1maJ3b1a5/RNgcWb+fgxtrc5M52jvEl7P7uL17C5ez+7i9ewuXs/u4vVsjUMEgcxcCzwJ+ABwLZDAVuBK4BTgaWNJriRJkiRNbzs1HcBUkZkbgfeXiyRJkiS1zB6s+p3ddACqldezu3g9u4vXs7t4PbuL17O7eD1b4DNYkiRJklQTe7AkSZIkqSYmWJIkSZJUExMsSZIkSaqJCdYQImJORPRExDURcW9E3B0RqyLibRGxywTbflhEnBkRN0TElojYEBErI+K1ERF1nYO2a8f1LNvLMSyPqft8pquImB0Rz42I90bElyNiXeV37qnpO7w/J0k7r6f35+SLiH0i4lURcX5EXBcRmyKiLyJuiYgLI+LPavgO789J0s7r6f05+SJiQUS8PyK+GhHXR8SdEbG1/Px+RLwnIuZO8Du8Pyucpn2QiDgQuJwdXzo8E1hYLi+LiDG9dHiIto8ELgb2KXfdC8wBnlUuL46IF2Rm30TOQdu183qWtgIbRjh+/zjb1YM9BfhGuxr3/px0bb2eJe/PyXMbO/5N0Uvx++9fLi+MiIuAv8jMza027v056dp6PUven5Pn1cAbK9u9wBZgLvCMcnlzeQ/9sNXGvT8fzB6sioiYAXyN4o/x3wLHZOZuwGxgCbAReDJwwTja3hP4OsU/fNcDR2XmHGA3YCnFv2iOBT424RMR0N7rWfGDzNxvhGXtxM5Cg/weuBQ4HfhLij8CJsz7szFtuZ4V3p+TZyfgCuANwKMzc9fM3B04CPhMWea5wLJWG/b+bETbrmeF9+fkuQJ4O/B0YO/yeu5BkQS9EvgdsC9wYXm/jZn35zAy06VcgNcAWS5PH+L4X1aOL26x7Q+W9TYDBw1x/NTy+P3AIU3/Ft2wtPl69pT1Lm/6PKfLAswYYt/a8jr0TLBt78/uup7en5N/PZ8zyvFPVf59+6gW2/b+7K7r6f05xRaKBGjger6sxbren0Ms9mDt6BXl52U5dBfpcmBNuX5ii20PlF+emWuGOP5Jii7VGcDLWmxbQ2vn9dQky8z+Njbv/TnJ2nw9Ncky87JRinymsr6wxea9PydZm6+npp4fVdYf2WJd788hmGCVImI28Mxy86KhymSRin+z3Dy2hbYfBxwwStv3AitbbVtDa+f1VHfx/pQmRW9lfcZYK3l/Tlnjup6ashZV1n851kren8MzwdruCWz/Pa4dodzAsf1amHHl8CHqj9T2oWNsV8Nr5/WsOiwiri1nzLm3nD3n0xHx5HG0pWZ4f3Yv78+p4+jK+jUt1PP+nJqOrqy3cj2rvD8bFBEzI2J+RCwFPl/u/gXFs+tj5f05DBOs7R5RWf/NCOWqxx4xbKmJtb1HROw+xrY1tHZez6p9KZK5gdkJDwFeC1wZER8aR3uafN6f3cv7cwqIiL0onsMAWJmZN7RQ3ftzipng9azy/mxARPRGRFL0Qq6hGMK3N/B9iufRW5npz/tzGCZY282prI805Wj12JxhS01e2xpau3/zm4B3AI8DZmXmPhQz5vwJcCUQwHsi4m0ttKlmeH92H+/PKSIiHkLxf8cfDvQBb2qxCe/PKaSG6wnen027Dbgd2FTZdxnw5sy8ucW2vD+HYYIljUNmXpCZp2fmjZm5tdx3X2Z+i+KdD6vKoj2tTnkqaWK8P6eUTwAnlOtvyMyfNhmMJmzC19P7s1mZOT+LafB3Bx4GnAL8AXBFRHyg0eC6iAnWdhsr67NHKFc9tnHYUpPXtobW2G+emb3Au8vN3YHFdbSrtvH+nEa8PydPRJxB8R4cgLdk5mfH0Yz35xRR0/Uckffn5MrM9Zl5JnAcxVTq74uIE0apVuX9OQwTrO1urazvP0K56rFbhy01sbbvKWdd0fi183qORXVa+INrbFf18/6cfrw/2ywiPgoMDPF6e2Z+fJxNeX9OATVez7Hw/pxkmXkF8L1y86QWqnp/DsMEa7ufA9vK9cNHKDdw7LbM3DDGtqszq4yl7evG2K6G187rqe7i/SnVKCJOB95ebr4jM8+YQHPenw2r+Xpq6hqYiOIxLdTx/hyGCVYpMzdTzKACRVfpg0REUDyECfCtFtq+ARh4cHC4tndj+3sIxty2htbO6zlGT6usD/XiPU0R3p/Tkvdnm5TDyE4pN9+RmadPpD3vz2bVfT3HyPuzGQO9hWMewuf9OTwTrB2dW34+JyKeOsTxF7P9H8DzWmx7oPySiJg/xPE3Uow37gcuaLFtDa0t17NMzEY6PhP4cLm5Cbh0rG2rMd6fXcL7sznlH+MDw8hOqfGPce/PBrTjenp/Tr6ImDGG330x8JRy8/IWv8L7cyiZ6VIuwE7A/1E86HcLxfsAoEhEXwzcXR77xhB1e8pjCcwf4viewG/L4z8Djiz37wK8nmK60wT+tenfoVuWdl1P4A+BbwMvBx5Z2b8zxQO5V1TqvqPp36GbFop3dexbWW4uf+ePDtq/+1ivZ3nc+7NLrqf3Z2PX8h8rv+tbWqzr/TnFlnZdT+/PRq7lfOBq4GSK/6kclWOPAt4F3Fv+5ncC+431epbHvT+H+t2bDmCqLeU/iGsq/zBtArZUtq8C9h6i3oj/AJZljgTuqJS7B7ivsn0xMLPp36CblnZcT4o32Gdl2Qz8btC17Ac+3PT5d9sCrB302w+3nDPW61kp4/3ZBdfT+7OR63jAoN/2tlGWU8Z6PStlvD+74Hp6fzZyPecP+s37yt/83kH7fwU8eYj63p/jWHZCO8jMtRHxJIoxx38OHARspcjK/wP4ZGbeN862r4yIw4B3UrxH4lEUf/BfSzGc7bOZuW2EJtSiNl3Pa8r2ng48keL/sO9F8R+K64CVwNmZeU0d56DJ4f3ZNbw/J99DBq0/bJTyu7f6Bd6fk6qd19P7c/LdCryEIrl9KsVLovelSGRvBn4KfAX4QmZuGc8XeH8+WJSZpyRJkiRpgpzkQpIkSZJqYoIlSZIkSTUxwZIkSZKkmphgSZIkSVJNTLAkSZIkqSYmWJIkSZJUExMsSZIkSaqJCZYkSZIk1cQES5K6WESsjYgsl6ObjqfTRMSjImJT+fv9T9Px6MEiYqeIuKm8RjdFxC5NxyRpejPBkiRpeGcAs4FtwHsbjkVDyMz7gb8rNx8DvLm5aCTJBEuSOkJEnFPpieppOp7pICKeAryk3PxSZv6syXg0ouXADeX6eyJi7yaDkTS9mWBJkjS0D1bWP9pYFBpVZiZwZrm5B3BKg+FImuZMsCSpi2Xm/MyMcrm86Xg6RUQcCRxbbv4wM69sMh6NyeeBu8r1N0XEnAZjkTSNmWBJkvRgb66sf7qpIDR2mdkLfKHcnAO8qsFwJE1jJliSJFWUz+8MPHvVB3ypwXDUmgsq63/TWBSSpjUTLEmawgYmtgBeUdn9/sqEFzssQ9QfdZr2iHhlpczllf1HRcSnI+L6iLi3XFZFxJuHmgo7IuZGxN9FxOqIuCcieiPilxHxqYg4aBznvl9EnBIRl0bEzWV7v4+In0XEWeUkFO3wEmDg/L6dmRuHiW/PyhTuGRHPH+sXRMSjI2Jbpe5Ro5Sv7beIiFkRcUJEfCwivhsRvy3b2xwRt0TExRHxjojYZ4zt9VTO45zK/j+MiM9ExHURcVd5/MIh6u8bEW8tv/c3EbElIraWdX4WEV+KiL+JiIPHEM4Pgd+W60+IiAVjOQdJqtNOTQcgSZpaIuIhFBM8nArEoMMLy+UlEXFsZt5b1jkG+A9g8B/lBwMnA38dES/KzG+O4fsDeCfwHmD3QYdnAnsBhwJvjIjzgJPL4WF1eVFl/aLhCmXm3RHxRbYPRXsd8LUxfsdr2P7b/jQzVw1VqO7fIiL+ElhGMYRuKPuXy7HA+yLibzPzc2M6o+3fsTvwr8Bfj6Hsi8t4hpr1b89yOZTymkTEUzPziuHay8yMiEuAE8tdfwFc1Ur8kjRRJliSNLVdXH4+EXhEuf5L4Bdt/M4PUSRXAHcC1wMJHMH2P8yfTjE19gkR8WyKxGIm0A9cA/weOJAiwYLiXVL/HRGHZeba4b44InYCzgFeNujQTcBvyu84vBLHicD8iDgmM+8bx7kO/v5dgUWVXf87SpVlbE+wnhcRj8jMW0f5jhns2CM55DNebfotHs2OydXvgLXAPcCs8vh+5bHdgc9GxC6ZuWykc6qGTTFM7wXl9l0U//xsLduunt+xFP8MVUfT3AqsAbaU338g8PDK8bGMvLmM7QnWnwDvHmPsklSPzHRxcXFxmeILxR/aWS49LdRbW6l39DBlXlkps4Hipbq/A14KzKiUmwV8olI2gRMo/ihO4JPAPoPafi6wsVL+vFHi/XCl7DbgLOCRg8rsTNFbdG+l7Jk1/c5HV9rcBDxkDHWurtR59xjKv6BSfjOw12T9FhQ9Yd+n6FV85DBlnkyR2FdjPHCENnsqZe8pP28HlgA7VcoFcHBl+/8q9X4APHmY9vcDXgtcATxtDL/voYN+tz3adV+6uLi4DLX4DJYkqWpvij+oj87ML2Zm/8CBzOzNzP/Hjr06/03Rw9CTmW/KzDurjWXmRew4I9+LImK3ob44Ip4MvKuy69WZuTQzbxnU5tbM/DRwHHB/ufv/RcQBrZzoMI6srF+XmdvGUKfau/PqcljfSF5bWf9SZt41uEAbf4uPZ+YzM3PZ4LYqbf6EIjH+n3LXrsAbRzmnAXMoEuo/zMzlmTkQE1n4VXl+j6LolYUiOTy+/N6h4rktM/89M58CDDmUcpAbKSYngSKpO3KEspJUOxMsSdJgH87Mn41wvJpQ7AL8nGJY4XDOp/ijG4qhgk8eptzb2f7fpQsy85yRgszM7wFnl5szKHplJurQyvovx1jnAoreLiiGwT1nuIIR8XCK5GXAcFPAt+W3yMxNQ+0foty2MoYBLxxLvdIHMvP6Uco8srL+s8z8/Rjj6h9DmfuBmyu7Dh2urCS1gwmWJGmwz4xyfPAkA58b6Q/fzOwDflrZ9YTBZcpnn6qTS3xstCBLX6isLx5jnZEcWFn/zVgqZOY9FM8SDXjtcGUpntcaeP75hsxcObjAVPktMvOXwB3l5mMjYq8xVLsfGMukGNWJOB5bnnOdqtfuwGFLSVIbOMmFJKnqV5m5fpQyvx20/aMxtFuts9cQx5/K9qnR787MK8fQJsC1lfUnR0Rk5oOmq2/BvpX1MfWqlJZRzAwI8OcRMTczN1QLlEMHX13Z9e/DtDUpv0U5jPCPgScB84A9ePDfBQMzFwbFJCt3jRLDzwcPEx2uHMVQ1NnAXODLEbG0TOrqUP3tH1pTm5I0JiZYkqSq20crkJmbBz1mNGodij+mB8we4vgTK+s7RcSo07kPYReKJOHucdQdUH0+bMtYK2Xmqoj4CcXwx5kUU5R/YlCx57B9Jr2twLnDNNfW3yIiDgE+TvHc1mjPi1XtOYYyY0qQMrM3Iv4NeFu56zjgpohYBXwbWAn8oOwdHI/qtRvymT9JahcTLElS1XimOm+1zlB/1Fffn7UbxfTa47EnE0uwqlpJPqDoxfpUuf5aHpxgvaay/pXM/N0w7bTtt4iIZwHfZHxJx8wxlBnypczDeDfF8L2/KLcDeEq5APRHxI8phl+ek8O88HkYrV47SaqNz2BJkqaCunoZJvrfteokEK0+F/QFihnxAA6PiKcOHIiIvYE/r5QdbnILaNNvERFzgC9V2r+XYmr95wOHUPR47ZyZMbAA61r8zrHMughAZt6XmS8Gjqd4ofPWQUVmAM8A/hlYExEvbSGO6rUb08QeklQXe7AkSVNBtaflR5n59IbiqPYq7d1KxczcGBFfAE4qd70W+HG5/nKK94hB8W6yb4/QVLt+i1cDDyvX76J4p9QNo9SZM8rxCcvMbwDfKKfvfybwbOAPKV5mPaMstg/wHxHRl5kXjqHZ6rUbrqdQktrCHixJ0lRQfY7rMY1FsWOPzf7jqF+dwn5JRAxMElEdHvjZUd6v1a7f4tjK+idGS67KhGevGr9/RJm5KTO/lZnvzcxFFMngO9n+/F4AZ4yxueq1a7UXTpImxARLkjpD9Q/ybny+pDoT4b4RcURDcVxXWW85ucnMq4CBWf92B14aEUcBA+fTz+jTmLfrt6i+fHgsL+x9Gg3+nZCZd2bmR9nxJcePjohHD1cHICJ2Ysep2a8brqwktYMJliR1hok8G9QJrmXHdxe9oaE4qlOiHxoRM4YtObxqL9br2PG9WBdl5i2j1G/Xb7FzZX0sU9m/qqbvnagLB20/bKhCFY9j+zT3CVxVd0CSNBITLEnqDLdV1pscQtcW5fuaPl7Z9eqI+MMGQrmC7S/BnQUcPo42/oPts+k9FTixcmy4d189oI2/RfVdZItGKhgRzwD+sobvHK79Vnph9xi0vWHIUtsdVVn/aWbWNaukJI2JCZYkdYbq/4U/NiIObSyS9vk34MZyfSfgKxHxp6NViohDI+LsiHjlRAPIzC3Adyu7jh5HG/cCF1R2DUxu8VtgxRibacdvcXllfelwQw8j4inAV2jv3wgvi4gvlMMnh1UO9/vHyq7fsv13Gc7RlfWLxxeeJI2fswhKUme4FFgPzKOYZvv/yhfb3kbxXA8AmfmnjURXg8zcFBF/BvyA4h1OewL/U74L6csUQ+fuoTj//She6ruY7b1MdQ0F+2+2TwhxHA9+n9VYLAP+ZtC+z2Xm/WOp3Kbf4myKSSN2pXg+7AcR8SngO2Vb+wMnAEsoZu+7GDgUeNRYYm7RThQ9ZH8ZEb8ov+tK4FaK6eP3AJ4EvAJ4QqXeP4w0QUjZM3ZMZdd/1xy3JI3KBEuSOkBm3hcRrwH+i6JHZAawsNmo6peZ10XE04CvAo8tdz+1XCbLf1K8e2km8EcRsUdm3tNKA5l5dURcwfaX5ibwmRbbqPW3yMzfRsRrgc9T9E7NBt5aLoNdC7yMHZ9Ja5fHMLZhr/8CnDVKmacCjyjXb8zMsUzmIUm1coigJHWIzPw6xWx0H6foobiLSu9Vt8jM6ynO822MPsX2vcDXKHpDzqnp+++iSLKgmCzhxeNs6v8q69/JzF+NI5Zaf4vM/AJFL9UvR2jjk8BTMvPOVuNtwWXARykSudEm3FgFvDAzl5bPp43kryrry4YtJUltFKP/u0qSpOZExOOAIymGR+5OMaPibcD1wLWZubUN3/lktg+za/llvxExi2K428ALb5dk5hdriKuW36KcHfGpFEML96SYOOJm4PLM3DxS3bpFxN7AHwAHU7xQeBeKRO9m4MrMHNN7rCJiJsUzWnuX9R/pBBeSmmCCJUnSECLiIopnsAAWZuaYh8tFxF8D55WbdwD7Z+Z9NYeoinII7cAsjf+Yme9qMh5J05dDBCVJGtrfVdbfMdZKEfEQ4JTKrn83uWqvcnKLgd98I3B6g+FImuZMsCRJGkI5QcLAsL6/iIjDRqtTDr37IMUMeABbKCbMUHu9FHh8uf6RNj8/JkkjcoigJEnDiIhHUTzfNBv4ylDT4EfEEoqpzWdSTGt+QOXwhzLzfZMQ6rRVvivr5xQzEf4SONQeQ0lNcpp2SZKGkZm/pnjX1EgeD7xwiP0rgQ/VHpR2UL5b7LGjFpSkSeIQQUmS6tNHMfX4e4BjMrOv4XgkSZPMIYKSJEmSVBN7sCRJkiSpJiZYkiRJklQTEyxJkiRJqokJliRJkiTVxARLkiRJkmry/wE7weAdYVU67wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from lifelines.plotting import plot_lifetimes\n",
    "\n",
    "l=0\n",
    "\n",
    "for cluster in [\"Cluster %d\"%(c+1) for c in [5,13]]:\n",
    "#for subtype in df_files[\"Subtype_Selected\"].unique():\n",
    "    subset = df_files[(~df_files[\"days_survival\"].isna()) & (df_files[\"days_survival\"]!=\"unknown\") & df_files.index.isin(pd.read_csv(f\"brca/trisbm/trisbm_level_{l}_clusters.csv\")[cluster].dropna())]\n",
    "    #subset = df_files[(~df_files[\"days_survival\"].isna()) & (df_files[\"days_survival\"]!=\"unknown\") & (df_files[\"Subtype_Selected\"]==subtype)]\n",
    "    \n",
    "    CURRENT_TIME = 3\n",
    "\n",
    "    actual_lifetimes = subset[\"days_survival\"].astype(float).to_numpy()/365\n",
    "    observed_lifetimes = np.minimum(actual_lifetimes, CURRENT_TIME)\n",
    "    death_observed = actual_lifetimes < CURRENT_TIME\n",
    "    \n",
    "    survival = 1 - death_observed.sum()/float(len(actual_lifetimes))\n",
    "\n",
    "    ax = plot_lifetimes(observed_lifetimes, event_observed=death_observed, figsize=(18,15))\n",
    "\n",
    "    survivals = []\n",
    "    for stat in range(1000):\n",
    "        subset = df_files[(~df_files[\"days_survival\"].isna()) & (df_files[\"days_survival\"]!=\"unknown\")]\n",
    "        sample_size=min((subset.shape[0],len(pd.read_csv(f\"brca/trisbm/trisbm_level_{l}_clusters.csv\")[cluster].dropna())))\n",
    "        subset=subset.sample(sample_size)\n",
    "        #subset = df_files[(~df_files[\"days_survival\"].isna()) & (df_files[\"days_survival\"]!=\"unknown\")].sample((df_files[\"Subtype_Selected\"]==subtype).sum())\n",
    "        actual_lifetimes = subset[\"days_survival\"].astype(float).to_numpy()/365\n",
    "        observed_lifetimes = np.minimum(actual_lifetimes, CURRENT_TIME)\n",
    "        death_observed = actual_lifetimes < CURRENT_TIME\n",
    "        survivals.append(1 - death_observed.sum()/len(actual_lifetimes))\n",
    "    Z = np.abs(survival-np.average(survivals))/np.std(survivals)\n",
    "    if Z > 2:\n",
    "        print(f\"*******{cluster} Z:{Z}*******\")\n",
    "    print(survival, np.average(survivals), np.std(survivals))\n",
    "    \n",
    "    ax.get_figure().set_size_inches(12,10)\n",
    "    ax.set_xlim(0, CURRENT_TIME*1.1)\n",
    "    ax.set_ylabel(\"patients\", fontsize=35)\n",
    "    ax.vlines(CURRENT_TIME, 0, len(subset), lw=10, linestyles='--')\n",
    "    ax.set_xlabel(\"time (years)\", fontsize=35)\n",
    "    ax.tick_params(labelsize=25)\n",
    "    ax.set_title(f\"{cluster} population at $t={CURRENT_TIME}$,\\n survival {survival.round(2)} (Z={round(Z,2)})\", fontsize=15)\n",
    "    #ax.set_title(f\"{subtype} population at $t={CURRENT_TIME}$,\\n survival {survival.round(2)} (Z={round(Z,2)})\", fontsize=15)\n",
    "    plt.tight_layout()\n",
    "    ax.get_figure().savefig(f\"survival_trisbm_{cluster}_t{CURRENT_TIME}.pdf\")\n",
    "    #ax.get_figure().savefig(f\"survival_topsbm_{subtype}_t{CURRENT_TIME}.pdf\")\n",
    "    plt.show()\n",
    "#print(\"Observed lifetimes at time %d:\\n\" % (CURRENT_TIME), observed_lifetimes)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e962b765",
   "metadata": {},
   "source": [
    "## Unitopic variate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47bb448d",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "import lifelines\n",
    "import sys\n",
    "from survival import fit_cox, add_group_to_subset, save_plot\n",
    "from topicpy.hsbmpy import get_file\n",
    "from scipy.stats import pearsonr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1228e10d",
   "metadata": {},
   "outputs": [],
   "source": [
    "l=1\n",
    "\n",
    "algorithm = \"trisbm\"\n",
    "kind = \"keyword\", \"metadatum\"\n",
    "#kind = \"word\", \"topic\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b46c110",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_clusters = pd.read_csv(\"%s/trisbm/trisbm_level_1_%s-dist.csv\"%(data_dir,kind[1]),index_col=1).drop(\"i_doc\",1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d020839d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "summaries = pd.DataFrame()\n",
    "mask = (~df_files[\"days_survival\"].isna()) & (df_files[\"cases.0.diagnoses.0.tumor_stage\"]!=\"not reported\")\n",
    "subset = df_files[mask]\n",
    "df_files = df_files.reset_index()\n",
    "subset = subset[[\"days_survival\",\"vital_status\", \"tumor_stage\", \"gender\", \"age_at_diagnosis\"]]\n",
    "for topic in df_clusters.columns:\n",
    "    top_set = add_group_to_subset(topic, subset, df_clusters.divide(df_clusters.sum(0),1), 0.35)\n",
    "    #print(top_set[topic].sum())\n",
    "    summary, cph, ax = fit_cox(top_set, topic)\n",
    "    if summary is not None:\n",
    "        if summary.at[summary.index[-1],\"-log2(p)\"]/(np.log2(10)) > 2:\n",
    "            print(topic,\"\\n\",summary.loc[summary.index[-1],[\"coef\", \"p\"]],\"\\n\")\n",
    "            summaries = summaries.append(summary)\n",
    "            #ax=None\n",
    "            if ax is not None:\n",
    "                ax.set_title(ax.title.get_text(), fontsize=35)\n",
    "                save_plot(ax, \"\", topic)\n",
    "        \n",
    "            figcph, axcph = plt.subplots(figsize=(18,15))\n",
    "            cph.plot(ax=axcph,elinewidth=15)\n",
    "\n",
    "            axcph.tick_params(labelsize=35)\n",
    "            axcph.set_xlabel(axcph.get_xlabel(), fontsize=35)\n",
    "            plt.tight_layout()\n",
    "            figcph.savefig(f\"survival_{topic}_HR.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d0960a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "summaries"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "beebc604",
   "metadata": {},
   "source": [
    "## Multivariate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd44934d",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext watermark\n",
    "%watermark -a Filippo_Valle -p pandas,scanpy,requests -m -v -g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4787c99c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from lifelines import CoxPHFitter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "becab8b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "algorithm = \"trisbm\"\n",
    "l = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "954f753e",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_covariates = pd.DataFrame(index=df_files[\"case_id\"])\n",
    "\n",
    "df_topics = pd.read_csv(\"%s/%s/%s_level_%d_metadatum-dist.csv\"%(data_dir,algorithm,algorithm,l), index_col=1).drop(\"i_doc\",1)\n",
    "df_topics_std = df_topics.divide(df_topics.sum(0),1)\n",
    "\n",
    "q = np.linspace(0, 1, 25)\n",
    "#q = [0, 0.35, 1]\n",
    "\n",
    "for topic in df_topics:\n",
    "    quantiles = np.argmin([np.abs(df_topics_std[topic]-q) for q in df_topics_std[topic].quantile(q).values], 0)\n",
    "    df_covariates = df_covariates.join(pd.Series(name=topic, index = df_topics_std.index, data = quantiles), how=\"outer\")\n",
    "    \n",
    "\n",
    "df_covariates[\"vital_status\"] = df_files[\"vital_status\"].values.astype(int)\n",
    "df_covariates[\"time\"] = df_files[\"days_survival\"].values\n",
    "#df_covariates[\"gender\"] = df_files[\"gender\"].values\n",
    "#df_covariates[\"age\"] = df_files[\"age_at_diagnosis\"].values\n",
    "#df_covariates[\"stage\"] = df_files[\"tumor_stage\"].values\n",
    "\n",
    "#mask = (~df_covariates[\"time\"].isna()) & (df_covariates[\"stage\"]!=\"not reported\")\n",
    "df_covariates = df_covariates[mask]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06548231",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_covariates[mask].head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e538dc7",
   "metadata": {},
   "outputs": [],
   "source": [
    "cph = CoxPHFitter()\n",
    "res = cph.fit(df_covariates.dropna(how=\"any\", axis=0), duration_col='time', event_col='vital_status')\n",
    "\n",
    "cph.print_summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64ad0aa7",
   "metadata": {},
   "outputs": [],
   "source": [
    "figcph, axcph = plt.subplots(figsize=(18,20))\n",
    "#cph.plot(ax=axcph, elinewidth=15, capsize=1)\n",
    "\n",
    "data = cph.summary\n",
    "coefs = data[\"coef\"]\n",
    "\n",
    "axcph.errorbar(x = coefs, \n",
    "               y = [t.replace(\"Metadatum\", \"Topic\") for t in data.index], \n",
    "               xerr=np.array([coefs-data[\"coef lower 95%\"],data[\"coef upper 95%\"]-coefs]).reshape((2,-1)),\n",
    "               lw=0,\n",
    "               marker=\"x\",\n",
    "               ms=25,\n",
    "               c=\"gray\",\n",
    "               elinewidth=30,\n",
    "               alpha=0.8\n",
    "      )\n",
    "\n",
    "axcph.vlines(x=0, ymin=0, ymax=data.shape[0], color=\"black\", ls=\"--\", lw=15, alpha=0.9)\n",
    "\n",
    "\n",
    "#data = cph.summary[cph.summary.index.isin([\"Metadatum 3\", \"Metadatum 7\"])]\n",
    "data = cph.summary[cph.summary[\"-log2(p)\"]/np.log2(10)>-np.log10(0.05)]\n",
    "coefs = data[\"coef\"]\n",
    "\n",
    "axcph.errorbar(x = coefs, \n",
    "               y = [t.replace(\"Metadatum\", \"Topic\") for t in data.index], \n",
    "               xerr=np.array([coefs-data[\"coef lower 95%\"],data[\"coef upper 95%\"]-coefs]).reshape((2,-1)),\n",
    "               lw=0,\n",
    "               c=\"red\",\n",
    "               elinewidth=15\n",
    "      )\n",
    "\n",
    "\n",
    "axcph.tick_params(labelsize=35)\n",
    "axcph.set_xlabel(\"Log(HR) (95% CI)\", fontsize=35)\n",
    "plt.tight_layout()\n",
    "figcph.savefig(f\"survival_HR.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9950a16f",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(cph.summary[\"coef\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e676a4e2",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}