LMM_model_shortChain_fattyAcid.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "e826903a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from paths import *\n",
    "import pandas as pd\n",
    "import os\n",
    "import statsmodels.api as sm\n",
    "from sklearn.preprocessing import LabelEncoder\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.patches as mpatches\n",
    "from tqdm import tqdm\n",
    "import joblib\n",
    "import multiprocessing as mp\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "c27230c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "NCORES = 8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "06ad9031",
   "metadata": {},
   "outputs": [],
   "source": [
    "def lmm_model(data_analyte, maxiter=1000):\n",
    "    le = LabelEncoder()\n",
    "    data_analyte['Disease_label'] = le.fit_transform(data_analyte['Group Name'])\n",
    "    data_analyte['Gender_label'] = le.fit_transform(data_analyte['Gender'])\n",
    "    data_analyte['site'] = data_analyte['Client Sample ID'].str[:3]\n",
    "    data_analyte['site'] = data_analyte['site'].astype(\"int\")\n",
    "    data_analyte = data_analyte[[\"Result\", \"Disease_label\", \"Age\", \"BMI\", \"Gender_label\", \"site\", \"house\"]]\n",
    "    '''\n",
    "    In the LMM model here, we consider site and house as crossed random effects (because each site holds household pairs)\n",
    "    This crossed random effect is implemented in the model (as shown below) based on these references:\n",
    "        https://www.statsmodels.org/stable/mixed_linear.html\n",
    "        https://stackoverflow.com/questions/59359568/from-2d-to-1d-how-to-pass-a-second-random-effect-in-mixed-model-python-statsm)\n",
    "        https://www.statsmodels.org/devel/generated/statsmodels.regression.mixed_linear_model.MixedLM.from_formula.html#statsmodels.regression.mixed_linear_model.MixedLM.from_formula\n",
    "    '''\n",
    "    data_analyte[\"group\"] = 1    # all in the case group \n",
    "    re_formula = \"~(1|site) + (1|house)\" \n",
    "    vcf = {\"site\": \"0 + C(site)\", \"house\": \"0 + C(house)\"}  # formula selected for the variance-covariance structure of the random effects in LMM\n",
    "    formula = \"Result ~ Disease_label + Age + BMI + Gender_label\"\n",
    "    model = sm.MixedLM.from_formula(formula, \n",
    "                                 data_analyte, \n",
    "                                 groups=data_analyte['group'], \n",
    "                                 re_formula=re_formula,\n",
    "                                 vc_formula=vcf)\n",
    "    print(\"Fitting the model...\")\n",
    "    mdf = model.fit(maxiter=maxiter)\n",
    "    print(\"Completed model fitting!\")\n",
    "    return mdf\n",
    "\n",
    "\n",
    "def lmm_model_2(data_, analyte_, maxiter=1000):\n",
    "    data_analyte = data_[data_[\"Analyte\"] == analyte_]\n",
    "    le = LabelEncoder()\n",
    "    data_analyte['Disease_label'] = le.fit_transform(data_analyte['Group Name'])\n",
    "    data_analyte['Gender_label'] = le.fit_transform(data_analyte['Gender'])\n",
    "    data_analyte['site'] = data_analyte['Client Sample ID'].str[:3]\n",
    "    data_analyte['site'] = data_analyte['site'].astype(\"int\")\n",
    "    data_analyte = data_analyte[[\"Result\", \"Disease_label\", \"Age\", \"BMI\", \"Gender_label\", \"site\", \"house\"]]\n",
    "    '''\n",
    "    In the LMM model here, we consider site and house as crossed random effects (because each site holds household pairs)\n",
    "    This crossed random effect is implemented in the model (as shown below) based on these references:\n",
    "        https://www.statsmodels.org/stable/mixed_linear.html\n",
    "        https://stackoverflow.com/questions/59359568/from-2d-to-1d-how-to-pass-a-second-random-effect-in-mixed-model-python-statsm)\n",
    "        https://www.statsmodels.org/devel/generated/statsmodels.regression.mixed_linear_model.MixedLM.from_formula.html#statsmodels.regression.mixed_linear_model.MixedLM.from_formula\n",
    "    '''\n",
    "    data_analyte[\"group\"] = 1    # all in the case group \n",
    "    re_formula = \"~(1|site) + (1|house)\" \n",
    "    vcf = {\"site\": \"0 + C(site)\", \"house\": \"0 + C(house)\"}  # formula selected for the variance-covariance structure of the random effects in LMM\n",
    "    formula = \"Result ~ Disease_label + Age + BMI + Gender_label\"\n",
    "    try:\n",
    "        model = sm.MixedLM.from_formula(formula, \n",
    "                                     data_analyte, \n",
    "                                     groups=data_analyte['group'], \n",
    "                                     re_formula=re_formula,\n",
    "                                     vc_formula=vcf)\n",
    "        mdf = model.fit(maxiter=maxiter)\n",
    "    except:\n",
    "        mdf = None\n",
    "    out_dict = {}\n",
    "    out_dict[\"analyte\"] = analyte_\n",
    "    out_dict[\"patient_sample_count\"] = data_analyte.shape[0]\n",
    "    out_dict[\"model\"] = mdf\n",
    "    return out_dict\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e57547f0",
   "metadata": {},
   "source": [
    "## Loading data of Targeted compounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "d763a3de",
   "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>Unique Tube Label ID</th>\n",
       "      <th>Client Sample ID</th>\n",
       "      <th>Client Matrix</th>\n",
       "      <th>Sample Description</th>\n",
       "      <th>Group Name</th>\n",
       "      <th>Treatment</th>\n",
       "      <th>Gender</th>\n",
       "      <th>Age</th>\n",
       "      <th>BMI</th>\n",
       "      <th>Race / Ethnicity</th>\n",
       "      <th>Analyte</th>\n",
       "      <th>Result</th>\n",
       "      <th>Unit</th>\n",
       "      <th>Dilution Factor</th>\n",
       "      <th>LLOQ</th>\n",
       "      <th>ULOQ</th>\n",
       "      <th>Analysis Comment</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1000</th>\n",
       "      <td>FS10431265</td>\n",
       "      <td>71701-0004</td>\n",
       "      <td>serum</td>\n",
       "      <td>serum</td>\n",
       "      <td>MS</td>\n",
       "      <td>Treated</td>\n",
       "      <td>F</td>\n",
       "      <td>52.0</td>\n",
       "      <td>19.422970</td>\n",
       "      <td>Caucasian/Hispanic</td>\n",
       "      <td>2-Methylbutyric acid</td>\n",
       "      <td>73.66340</td>\n",
       "      <td>ng/mL</td>\n",
       "      <td>1</td>\n",
       "      <td>25</td>\n",
       "      <td>2500</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1001</th>\n",
       "      <td>FS10431273</td>\n",
       "      <td>71601-0150</td>\n",
       "      <td>serum</td>\n",
       "      <td>serum</td>\n",
       "      <td>MS</td>\n",
       "      <td>Off</td>\n",
       "      <td>F</td>\n",
       "      <td>51.0</td>\n",
       "      <td>25.501750</td>\n",
       "      <td>Caucasian/Hispanic</td>\n",
       "      <td>2-Methylbutyric acid</td>\n",
       "      <td>48.47522</td>\n",
       "      <td>ng/mL</td>\n",
       "      <td>1</td>\n",
       "      <td>25</td>\n",
       "      <td>2500</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1002</th>\n",
       "      <td>FS10431281</td>\n",
       "      <td>71602-0150</td>\n",
       "      <td>serum</td>\n",
       "      <td>serum</td>\n",
       "      <td>Control</td>\n",
       "      <td>Control</td>\n",
       "      <td>M</td>\n",
       "      <td>52.0</td>\n",
       "      <td>34.970939</td>\n",
       "      <td>Caucasian/Hispanic</td>\n",
       "      <td>2-Methylbutyric acid</td>\n",
       "      <td>38.87322</td>\n",
       "      <td>ng/mL</td>\n",
       "      <td>1</td>\n",
       "      <td>25</td>\n",
       "      <td>2500</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1003</th>\n",
       "      <td>FS10431289</td>\n",
       "      <td>71601-0164</td>\n",
       "      <td>serum</td>\n",
       "      <td>serum</td>\n",
       "      <td>MS</td>\n",
       "      <td>Treated</td>\n",
       "      <td>M</td>\n",
       "      <td>56.0</td>\n",
       "      <td>26.544467</td>\n",
       "      <td>Caucasian/Hispanic</td>\n",
       "      <td>2-Methylbutyric acid</td>\n",
       "      <td>129.41031</td>\n",
       "      <td>ng/mL</td>\n",
       "      <td>1</td>\n",
       "      <td>25</td>\n",
       "      <td>2500</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1004</th>\n",
       "      <td>FS10431297</td>\n",
       "      <td>71602-0164</td>\n",
       "      <td>serum</td>\n",
       "      <td>serum</td>\n",
       "      <td>Control</td>\n",
       "      <td>Control</td>\n",
       "      <td>F</td>\n",
       "      <td>53.0</td>\n",
       "      <td>31.176140</td>\n",
       "      <td>Caucasian/Hispanic</td>\n",
       "      <td>2-Methylbutyric acid</td>\n",
       "      <td>85.74194</td>\n",
       "      <td>ng/mL</td>\n",
       "      <td>1</td>\n",
       "      <td>25</td>\n",
       "      <td>2500</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     Unique Tube Label ID Client Sample ID Client Matrix Sample Description  \\\n",
       "1000           FS10431265       71701-0004         serum              serum   \n",
       "1001           FS10431273       71601-0150         serum              serum   \n",
       "1002           FS10431281       71602-0150         serum              serum   \n",
       "1003           FS10431289       71601-0164         serum              serum   \n",
       "1004           FS10431297       71602-0164         serum              serum   \n",
       "\n",
       "     Group Name Treatment Gender   Age        BMI    Race / Ethnicity  \\\n",
       "1000         MS   Treated      F  52.0  19.422970  Caucasian/Hispanic   \n",
       "1001         MS       Off      F  51.0  25.501750  Caucasian/Hispanic   \n",
       "1002    Control   Control      M  52.0  34.970939  Caucasian/Hispanic   \n",
       "1003         MS   Treated      M  56.0  26.544467  Caucasian/Hispanic   \n",
       "1004    Control   Control      F  53.0  31.176140  Caucasian/Hispanic   \n",
       "\n",
       "                   Analyte     Result   Unit  Dilution Factor  LLOQ  ULOQ  \\\n",
       "1000  2-Methylbutyric acid   73.66340  ng/mL                1    25  2500   \n",
       "1001  2-Methylbutyric acid   48.47522  ng/mL                1    25  2500   \n",
       "1002  2-Methylbutyric acid   38.87322  ng/mL                1    25  2500   \n",
       "1003  2-Methylbutyric acid  129.41031  ng/mL                1    25  2500   \n",
       "1004  2-Methylbutyric acid   85.74194  ng/mL                1    25  2500   \n",
       "\n",
       "     Analysis Comment  \n",
       "1000              NaN  \n",
       "1001              NaN  \n",
       "1002              NaN  \n",
       "1003              NaN  \n",
       "1004              NaN  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "'''\n",
    "Note that, here we are taking all samples without applying the threshold of \"Analysis Comment\". This is because,\n",
    "if we apply the threshold, each analyte will have different number of patient samples.\n",
    "Hence, there could be a sample effect on the parameter estimation, for e.g. one analyte's parameter may be estimated\n",
    "from N1 samples and the other analyte's parameter estimated from N2 samples. If N2<<N1, then the estimated model parameters\n",
    "maynot be comparable. Hence, to ensure that all analyte's model parameters are estimated from the same set of patient samples,\n",
    "we are not applying this \"Analysis Comment\" threshold.\n",
    "'''\n",
    "\n",
    "sample = \"serum\"\n",
    "\n",
    "filename = SHORT_CHAIN_FATTY_ACID_DATA_FILENAME\n",
    "mapping_filename = \"short_chain_fatty_acid_spoke_map.csv\"\n",
    "file_path = os.path.join(DATA_ROOT_PATH, filename)\n",
    "mapping_filepath = os.path.join(OUTPUT_PATH, mapping_filename)\n",
    "data = pd.read_excel(file_path, engine='openpyxl')\n",
    "data = data[data[\"Client Matrix\"]==sample]\n",
    "mapping_data = pd.read_csv(mapping_filepath)\n",
    "analytes = mapping_data[\"name\"].unique()\n",
    "data.head()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2873ced",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Data will be pre-processed as follows:\n",
    "    1. Remove patients with missing values (remove their paired partner too even if they have all values. This is to balance the data in terms of the diet)\n",
    "    2. Remove patients with outlier values in the specific analyte concentration (also their paired partner, for the same reason mentioned above)\n",
    "\n",
    "Build 2 models (for each sample serum/feces):\n",
    "    Model 1 - Stringent model \n",
    "        This model uses the fully pre-processed data as mentioned above. Hence, it does not use outlier and missing values.\n",
    "    Model 2 - Less stringent model\n",
    "        This model uses partially pre-processed data. i.e. data obtained after removing missing values (point 1 of pre-processing).\n",
    "        Hence, it uses outlier values and no missing values.\n",
    "\n",
    "Even though we compute Model 1 for all analytes, for the downstream analysis, we may only consider those Model 1 models of analytes for which there exists at least 800 patient samples after pre-processing.\n",
    "\"\"\"\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ac3beaa",
   "metadata": {},
   "source": [
    "## Check count of patient samples for each analyte after pre-processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "171442d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "data['house'] = data['Client Sample ID'].str[:3] + data['Client Sample ID'].str[-4:]\n",
    "data['house'] = data['house'].astype(\"int\")\n",
    "data_exclude_outlier_threshold_column = data.drop(\"Analysis Comment\", axis=1)\n",
    "house_to_exclude = data_exclude_outlier_threshold_column[data_exclude_outlier_threshold_column.isna().any(axis=1)].house.values\n",
    "data_nan_removed = data[~data[\"house\"].isin(house_to_exclude)]\n",
    "\n",
    "data_outlier_removed = data_nan_removed[data_nan_removed[\"Analysis Comment\"].isna()]\n",
    "data_outlier_removed.drop(\"Analysis Comment\", axis=1, inplace=True)\n",
    "\n",
    "\n",
    "#Removing unpaired house for each analyte after removing the outliers, so that we get paired MS and control patients\n",
    "\n",
    "#Step 1: Group by \"Analyte\" and \"house\" and count the number of occurrences of each group\n",
    "group_counts = data_outlier_removed.groupby(['Analyte', 'house']).size()\n",
    "\n",
    "#Step 2: Get the list of index values where the count is 2\n",
    "valid_index = group_counts[group_counts == 2].reset_index().set_index(['Analyte', 'house']).index\n",
    "\n",
    "#Step 3: Filter the DataFrame to keep only the rows with paired \"house\" values for each \"Analyte\"\n",
    "data_outlier_removed = data_outlier_removed.set_index(['Analyte', 'house']).loc[valid_index].reset_index()\n",
    "\n",
    "analyte_sample_count_after_pre_processing = data_outlier_removed.groupby('Analyte')['Unique Tube Label ID'].nunique().to_dict()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "98f959bf",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'2-Methylbutyric acid': 996,\n",
       " 'Acetic acid': 846,\n",
       " 'Butyric acid': 856,\n",
       " 'Hexanoic acid': 340,\n",
       " 'Isobutyric acid': 996,\n",
       " 'Isovaleric acid': 994,\n",
       " 'Propionic acid': 398,\n",
       " 'Valeric acid': 8}"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "analyte_sample_count_after_pre_processing"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cba9adc4",
   "metadata": {},
   "source": [
    "## Plotting the final count of patient samples after pre-processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ac8554fb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoQAAAGJCAYAAAATwSIMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAACf2ElEQVR4nOydd5gURdPAf8Vx5AySQRAVELOYE6gYXsWEATNmzKCYE0YMICoGzAq+ioqiiOk1gIo58YkBBRVUQJAsEg64+v6o3ru5vQUOdvf2dqnf88xzNz09u929PT3V1VXVoqo4juM4juM4Gy6VMl0Ax3Ecx3EcJ7O4QOg4juM4jrOB4wKh4ziO4zjOBo4LhI7jOI7jOBs4LhA6juM4juNs4LhA6DiO4ziOs4HjAmGOICJdRWS8iCwVkXkiMlxEmiTIt62IvCkii0VkkYiMFpFNV/OZHUXkBRGZEz73JxG5KP21cRwnlfj44DjO2nCBMAcQkT2B/wELgB7ARcBewLsiUjWSbzPgQ6AucAJwKtAG+EBEGsd9ZmfgM6AqcAbwH2AQkJfe2jiOk0p8fDDSIRRH7rlCRFRExqevBo6TXsQDU2c/IvIONnB3UNWVIa0z8AVwnqo+ENIeBY4C2qjqgpDWEpgC3Kuql4W0SsB3wE+qekT51sZxnFTi40ORUPwe8BbwANAQuBn4B9hBVZeHfJsBX2P1uw2oDFwPNAK2VdXZCT57E+Bb4F9gsqrukfYKOU4acIEwBxCRxcBwVT0nLn0O8I2qdgvn3wF/qOpBcfm+BOqq6mbhfB/gXWAvVf2wPOrgOE568PEh9UJx3Ge/BUwF2gOVXSB0shVfMs4NVgEFCdKXA1uWMV87EakWzmMDWjUR+VREVojIbBG5V0Sqp6zUjuOUBz4+wC7A2zFhEEBVvwTmAkfE5fskJgyGfH9iGsNS2lAROR7YHrgyPcV2nPLDBcLc4CdsICtCRDYGmgEN4vLtICL5kXy1gU6AAPVDcvPw9znM9qgbcAdmK/RMGsrvOE768PEh9UIxIlIfGAxcpqrzUlhWx8kILhDmBvcAO4nIzSLSWEQ6AMOBwnDEuBdoAQwVkRbhpfAEUCtcj+WN9YunVfU6VR2nqgOBG4DDRaRjuivkOE7K8PEh9UIxwJ3Az8CT6Smy45QvLhDmAKr6X8xA+hJgFvADMB14HZgZyTceOA+zkfkTs3upCzyFzYpjs9y54e/bcV/1v/B3u1TXwXGc9ODjA5BioTg4qZwMnKNuiO/kCC4Q5giqei3mCbc10ExVjwM2A8bH5XsAaIwtk7QOBuXNgc9UdUXI9v1avq5wLdcdx6lAbOjjQxqE4oeAx4A/RaSeiNTDPJLzwnlROB/HyRZcIMwhVPVfVZ2oqrNE5ECgAzA0Qb7lqvq9qv4hIlsB+wEPRrK8gdnMHBB364Hh75dpKL7jOGlkQx8fUiwUdwR6A/Mjx+7YsvR8oIRHt+NkA5UzXQAneURkO+AgLH4WmBfgpcAdqvpxJF9LbKD6GBvQO2PecS+p6rOxfKo6V0QGANeKyCIsfldn4DrgKVWdkv5aOY6TCnx8KEZV/wUmAkSE4tMT5FtO0IRGhOKTI1m6Jvj4u7HA3BdgYWocJ7tQVT/W8zjggAO0IvDdd9/p7rvvrnXr1tVq1arpdtttp48//nipfH/99Zfuu+++2rBhQ61SpYp27NhRBw4cqCtWrCiVt7CwUAcNGqTt2rXT/Px8bd26tV577bVaUFBQHlVyHCdF+Pig+vXXX+stt9yib7zxhr7xxht69dVXa5UqVfSyyy4rke+PP/7Qq666SseMGaNvv/22DhgwQOvUqaPHHHPMWr9j77331t133z1dVXCcVLFamcYDUydB586d9csvK+TqiOM4jhP4/vvvOfvss/nuu+9Yvnw5HTt25IILLuDUU08tkW/WrFmccMIJTJgwgX/++Yd27dpx+umnc9FFF1G58poX1Lp06cLKlSsZP953r3MqNLLaCy4Qrj8uEDqO4ziOk0WsViB0pxLHcRzHcZwNHBcIHcdxHMdxNnAqlEAoIpuKyEMi8q2IrBKRcQnyiIhcJSJ/iMhSEflARLZNkG8LEXlXRJaIyAwRuVFE8tbnsxzHcRzHcXKZCiUQYtsD/QfbPujn1eS5ArgWuB3oDiwG3hGRprEMYY/JdwAFDgNuxAKS3rCun+U4juM4jpPrVDSB8FVVbaWqR5MgGn7YWPwKYICq3qeq7wBHY4Lf+ZGsvYHqwJGq+raqDsWEwYtFpM46fpbjOI7jOE5OU6ECU6vq2rY82g2oAzwfuedfEXkVC7x6TUg+CHhLVRdF7h2BaQL3Bl5dh89yHMdxnA2WRYsWMXv2bFasWLH2zE5GyM/Pp3HjxtSpU2e9P6NCCYRloAOwCpgcl/4jcGxcvveiGVT1dxFZEq69ug6f5TiO4zgbJIsWLWLWrFm0aNGC6tWrI7LaqCVOhlBVli5dyvTp0wHWWyisaEvGa6M+sFhVV8WlzwdqiEiVSL4FCe6fH66ty2eVQETOEpEvReTLv//+e33q4DiO4zhZwezZs2nRogU1atRwYbCCIiLUqFGDFi1aMHv27PX+nGzTEGYcVX0YeBgsMHW6v++ln2am+yvSwpHtm2W6CI6T8/j4YHg7pI8VK1ZQvXr1TBfDKQPVq1dPalk/2zSE84Fa8eFjMG3fElUtiOSrm+D++uHaunyW4ziO42ywuGYwO0j2d8o2gXASkAdsGpfeIVyL5usQzSAirYAakXxl/SzHcRzHcZycJtsEwo+BRVh4GABEpAYWQ/CNSL43gANEpHYk7VhgKfD+On6W4ziO4zhOTlOhbAiDQPafcNoCqCMiR4Xz11V1iYjcBlwrIvMxTd7FmGA7JPJRQ4ELgZdE5HZgE6A/cFcsFI2qLivjZzmO4ziOE+G1KbNYvmptkeJST9W8Shy8aZP1vv/FF1/k/vvv5+uvv2bp0qVsvPHGHHLIIfTr14/mzZunrJw///wzzzzzDH369KFevXop+cxx48bRtWtXJk6cyJZbbpmSz4xSoQRCoDHwQlxa7LwtMBW4DRPargQaAl8C3VR1VuwGVZ0vIvsC92EhZhYAgzGhMMpaP8txHMdxnJJkQhhM9nsvueQS7r77bk499VT69u1LnTp1+OGHHxg6dCi//fYbo0aNSlk5f/75Z2644QZ69eqVMoEw3VQogVBVpwJrtIpUVQVuCcea8v0A7JOKz3Icx3EcJ3t59dVXueuuu3jsscc47bTTitL33ntvzjrrLP73v/9lrGzLli2jWrVqGfv+GNlmQ+g4Thn46KOP2H///WncuDG1a9dm++235/HHHy+RR0QSHhMmTCj1edOnT+e0006jadOmVK1albZt23LllVeWU20cx3GSY/DgwWy//fYlhMEYeXl5HHTQQQDMmTOHU045hYYNG1KjRg26dOnCl19+WSJ/mzZt6NevH4MHD6Zly5bUr1+fnj17smDBAsCWdrt37w5A27ZtERHatGkDwJNPPomI8Pnnn9OlSxeqV6/OnXfeCcB7773HzjvvTLVq1WjSpAnnnnsuixcvTlOLlKZCaQgdx0meb7/9lv32249ddtmFRx55hBo1ajBy5EhOP/10li9fzjnnnFOUt1evXpx99tkl7t98881LnE+dOpXdd9+dtm3bcu+999KkSROmTp3KlClTyqU+juM4ybBixQo+/vhjLrnkkrXmPfzww5kyZQoDBw6kUaNG3HnnnXTt2pVvvvmGTTctDkry/PPPs/XWW/Pwww/z559/cvHFF3PVVVfxwAMPsP322zNw4ED69evHSy+9RLNmzahatWqJ7znuuOM499xzuf7666lXrx7ff/89Bx54IN26dePFF1/kjz/+4IorruDXX3/lzTffTHmbJMIFQien+Oijj7jhhhuYMGECS5cuZbPNNuP8889POCsEuO2227jyyivZfffdGT9+fKnr06dP59prr+X1119n/vz5NG/enJ49ezJgwIB0V2W9GTFiBKtWreLVV1+lVq1aAHTr1o1vv/2WYcOGlRAIW7RowS677LLGz+vduzctWrRg7Nix5OfnA7bM4jiOkw3MnTuX5cuX07p16zXme/PNN/noo48YN25c0Ri3zz770KZNG+68804eeuihorz5+fm8/PLLVK5sYtQPP/zAiBEjeOCBB6hTpw7t27cHYLvttivSDka58MILueiii4rOe/bsycYbb8zo0aPJy7PwyA0aNODYY4/lk08+Ydddd02qDcqCC4ROzrAumjGAX3/9lZtvvpnGjRsn/Lxs1YwVFBSQn59faneBunXrMn/+/NXclZhffvmFt956i2HDhhUJg47jONnI2gI3f/755zRu3LjEhLdmzZoccsghpRQGXbt2LRIGAbbYYgtmz57NihUryjRWHnzwwaW++6ijjioSBgF69OhB5cqVGT9+fLkIhG5D6OQMUc3YYYcdRrdu3XjooYfYZZddGDZsWKn855xzDieccAIdO3ZM+HlRzdgxxxzD3nvvzSmnnMJNN92U7qokRa9evQCbgc6YMYMFCxbwyCOP8O6779K3b98SeR988EGqVq1KjRo12Gefffjwww9LXP/oo48A2xKpW7duVK1alfr163PyySczd+7ccqmP4zhOMjRs2JCqVavy+++/rzHfzJkzEyoImjRpwrx580qkxXsOV6lSBVVl+fLlZSpTkyYlQ+fMnDmzVFpeXh4NGzYs9d3pwgVCJ2dYk2assLBkqIJnnnmGr7/+erVLvzHN2AUXXJB1mrEtt9yScePG8corr9CiRQvq16/Peeedx9ChQ+nZs2dRvhNPPJEHHniAd955h4cffpi5c+eyzz77MG7cuKI8M2bMAOC0005j880354033uD222/ntdde44ADDijVro7jOBWN/Px8dt99d95666015mvWrBmzZ88ulT5r1iwaNGiQ0jLFaysTffeqVauYO3duyr97dbhA6OQMZdWMzZ8/n759+3LHHXes9kHLZs3Y5MmT6dGjB506deLVV1/lnXfeoXfv3vTu3Zv//ve/RfmGDx/Osccey5577smJJ57I+PHjad68Oddcc01RnpjA16VLF+6//3722WcfzjrrLB544AG++uqrtQ6wjuM4FYE+ffrw5Zdf8tRTT5W6VlhYyJtvvsnOO+/M7Nmz+eCDD4quLVmyhNdee4099thjnb6vSpUqgIWUKQs777wzo0aNYtWqVUVpL730EitXrlzn715f3IbQyRlimrEjjjiCBx54ALCZYbxm7NJLL2XzzTcvEiATEdWMnXTSSVx55ZVMmTKFK6+8kh9++IHPP/+cSpUq5nzqqquuIj8/nzFjxhRpN/fdd1/mzp3LRRddxHHHHZew7LVr1+bggw/mscceK0pr2LAhYE4pUfbff38Avvnmm6JwDY7jOBWV7t27c/HFF3P66afz0Ucfcdhhh1GrVi0mTZrE0KFDadOmDaNGjWK33Xbj2GOP5bbbbqNhw4YMHDiQpUuXcumll67T98WcSh566CF69uxJjRo12GqrrVab/5prrmG77bbj8MMP55xzzuHPP//k8ssv54ADDigX+0FwgdDJIaKasaFDh1K9enVeeeUVevfuTbVq1TjhhBP48MMPGTZsGF9//fUaDYzjNWNg3mZ169alZ8+evPXWWxVWEJo4cSLbbLNNqaXunXbaiWeeeYbZs2fTtGnT1d4fbZdOnTqt8bsqqlDsOE56qZpXKWNb160vgwYNYrfdduO+++7j+OOPZ+nSpbRp04ZDDz2Ufv36AfDyyy9zySWX0KdPH5YtW8ZOO+3Ee++9VyLkTFnYeOONGThwIPfeey9DhgyhZcuWTJ06dbX5O3XqxBtvvMFVV13FkUceSZ06dTjuuOO444471ru+64oLhE7OUBbN2Nlnn83pp59Oy5Yti4KIrly5klWrVrFgwQKqV69O1apVs1oz1rRpUyZMmEBBQUHRsgXAZ599RrVq1Va7TL5o0SLGjBnDTjvtVJS2yy670LRp0yJ7yhixuFg77rhjmmrhOE5FJpn9hDNJjx496NGjx2qvb7TRRgmdEKMkEux69epVatXpkksuKRX7MFG+GPvuuy+fffbZar+3S5cu2AZr6cEFQidnKItm7Mcff+THH39k6NChpe6vX78+gwcPpk+fPlmtGTv//PM5+uij6d69O+eeey7Vq1dn9OjRPPvss/Tt25cqVaowcOBAfvrpJ7p27Urz5s2ZNm0aAwcO5K+//iphZ1i5cmVuu+02evXqRe/evTnyyCOZMmUKV199NV26dGGffda4O6TjOI6TJbhA6OQMZdGMjR07ttR9ffr0YdWqVQwZMqRoWSCbNWNHHXUUr7/+OrfffjtnnHEGy5Yto127dtx///1Fu5K0b9+eUaNGMWrUKBYuXEidOnXYfffdeeyxx0poCAFOOeUUKlWqxO23384TTzxBgwYNOPHEExkwYMBa43o5juM42YELhE7OUBbNWJcuXUrdV69ePVauXFniWrZrxg466KA1Lml37969aK/NsnDSSSdx0kknpaJojuM4TgXEBUInZyiLZmxdcM2Y4ziOs6Gw3gKhiOQBVVV1SVz6PsBhwBLgYVX9LbkiOk7ZWZtmLBHRQMzxuGbMcRzH2RBIxjJ+IDBPROrGEkSkJ/A2cAFwOfC5iLRKroiO4ziO4zhOOklGINwLGKuqCyNp1wMLgJOBy4B6wMVJfIfjOI7jOI6TZpIRCFsBU2InIrIJ0B4YoqpPq+pA4A3gwOSK6DiO4ziO46STZATCOsCiyPnugAJvRtK+B1om8R2O4ziO4zhOmklGIJwJtI2c7wcsBb6KpNUCVibxHY7jOI7jOOtN//79ERE222yzhNc322wzRIT+/fsXpT355JPssMMO1K5dm/r167Pddttx8cW5bQGXjED4KXCoiBwiIvsBRwHvqeqKSJ62wPRkCug4juM4TgWjaVMQKf9jDfuwr4lq1arx22+/8eWXX5ZI/+KLL5g6dSrVqlUrShswYABnnHEGBxxwAC+99BLDhg3jsMMOY/To0Uk1WUUnmTiEt2LhZV4J54XALbGLIlIN2BMYmcR3OI7jOI5T0Zg1K6u+t2bNmmy//faMGDGCzp07F6WPGDGCffbZh6++Kl7cvO+++zj77LO59dZbi9K6d+/O9ddfv/7lzgLWWyBU1YkisjNwSkh6TlW/iGTZDngPeDaJ8jmOE+Gln2ZmugjrxZHtm2W6CI7jbOD07NmT/v37c+eddyIiqCrPP/88N954YwmBcMGCBTRNoInM9Q0JklkyRlUnqmq/cHwRd+0TVT1CVd9JroiO4ziO4zjJceSRRzJr1izGjx8PwIcffsjff//NkUceWSLf9ttvz5AhQ3jqqaeYO3duJoqaEVK2dZ2I1AdqqeofqfpMx4nhmjHHcRwnGerVq8eBBx7IiBEj2HPPPRkxYgQHHnggdevWLZHv/vvv5/DDD6dXr16ICB07dqRHjx7069ePOnXqZKj06ScpDaGI1BKRQSLyFzAH+C1ybWcReV1Etk+2kI7jOI7jOMnSs2dPRo4cyfLlyxk5ciQ9e/YslWfrrbfmxx9/ZPTo0Zx77rmoKjfddBOdO3dm8eLFGSh1+bDeAmHYsu4ToC8wA/gRiC6wT8ScSo5LpoCO4ziO4zip4NBDD2Xx4sVcffXV/Pvvv3Tv3j1hvqpVq9K9e3fuu+8+fvjhBx599FEmT57MY489Vs4lLj+S0RBeDXQCeqnq9sAL0YuqugR4H9g3ie9wHMdxHMdJCTVr1uSQQw5h8ODBdO/enZo1a5bpvtNPP50GDRowadKkNJcwcyQjEB4JvKWqw9aQZxrQIonvSIiI9BSRr0VksYhMF5FhItI8Lo+IyFUi8oeILBWRD0Rk2wSftYWIvCsiS0RkhojcKCJ5qS6z4ziO4ziZ55xzzqF79+707t074fXZs2eXSvv7779ZuHAhTZo0SXfxMkYyTiUtgRfXkmcxUHctedYJETkUC2VzP3Ap0Ay4GXhNRHZQ1cKQ9Qrg2pBnEnAx8I6IbKmqf4XPqg+8A/yAxVRsBwzCBOVrUllux3Ecx3EyT5cuXejSpctqr2+11VYcdthh7L///jRu3Jhp06YxcOBAatSowSmnnLLa+7KdZATCf4DGa8nTFnM2SSXHA1+r6vmxBBFZhAXIbg/8GIJiXwEMUNX7Qp5PgKnA+RQLe72B6sCRqroIeFtE6gD9ReSOkOY4juM4TpQmTTITnLocNHTXXXcdr7zyChdeeCHz5s2jadOm7Lbbbjz33HO0bdt27R+QpSQjEH4BHCIitVX1n/iLItIM+A8wJonvSEQ+sDAubUHsa8Pf3YA6wPOxDKr6r4i8ChxEsUB4ELbsHRX8RgC3A3sDr6a05I7jOI6TC/z1V6ZLUGb69+9fYp/iRMyZU6y7Ou+88zjvvPPSXKqKRzI2hPcADYHXRaRj9EI4fwGoBtybxHck4nFgTxE5WUTqiMjm2JLxe6r6Q8jTAVgFTI6798dwjUi+Ehaiqvo7sCQun+M4juM4Ts6y3gKhqr4F3ADsDnwHXAkgInPC+W7Alar6cQrKGf3e14BewMOYpvAnIA/oEclWH1isqqvibp8P1BCRKpF8CxJ8zfxwrRQicpaIfCkiX/7999/rW420c9MZx9OjQ3Oeufv2Eum//fgdN51xPCdsvykn7rA5A845hZnTfit1/98z/mTI5RdxdtfOHLfNJpx/wB48c/ftLFuypLyq4DhOGvCxwXGcRCS1U4mq3iAiHwAXArtgGkMFXgcGq+p7yRexJCLSFRiKaSjfAJoA/YFRIrJfAiEwpajqw5gwSufOnTWd37W+fDhmFFN/+qFU+oypv3LNCUfQerMOXHTnfaxauYrn77+La088kkEvv03dho0AWLZkCTec2pNVK1fQ88LLaNSsBb98N4Hnhgxi5rRfuWTwQ+VdJcdxUoCPDcXcdMbxTBg/jh69L+L4PpcXpf/243c8PehWJn39OSKV6LTTrvS6oj/NNi62HXtuyECev/+uhJ+bX6UqI74tLUg7TkUn6a3rVHUsMDYFZSkrg4DRqlr0BIvIBGzp9zDgJUzDV0tE8uIExPrAElUtCOfzSewFXT9cyzoWL1zAk7f1p9cV/bm7X0kbiJcfvZ9KeXlc88jT1Kxj1d58m+04b//deeXxBzn50msBmPT158yc9ivXPvoM2+7RBYCtdtmdfxYuYPTjQ1m+dAlVq9co13o5jpMcPjYUk6xgvN/Rx7Pdnl1L3Lts6RJuPvMEOu+zf7nUwXFSTVJb12WIDsCEaIKq/gQsxcLGgAmHecCmCe6N2gxOIs5WUERaATXi8mUNwwfeQqvN2rPnIUeUuvbz/31N+213KBrwARo2bU7rzdrz2TtvFqWtXLECgBq1ape4v2btOmhhIVoh9aKO46wJHxuMIsH48utLXYsKxjvteyC7HnAw1zw8nMULF/DK4w8W5WvYtDmbb7tDiWPe7FmsWrmSrocfXZ7VKRc0G35YJ+nfKRsFwmlAif2RgxNLdSysDMDHwCLg6EieGkB3bJk5xhvAASISHd2OxYTL91Nd8HTz41ef8f4rIznz2lsTXq9UKY/K+fml0itXqcKs36dSsHwZAFvvtifNNt6E4QNv4Y8pP7P033+Z+Ol4Xhv+GPv3PJlqNSq+BsBxnGJ8bCgmFYJxIsa9/Dz1Gm1UpDnNFfLz81m6dGmmi+GUgaVLl5Kf4DkuK2VeMhaRQsw+cF1RVU16aTrCUGCwiMyg2IbwOkwYfD184TIRuQ24VkTmUxyYuhIwJO6zLgReEpHbgU0we8S7si0G4YqCAoZefzmHntabFpvEK0aN5m3b8dM3X7JyxYqiwX/p4sX8MeVnVJV/Fy6kSuNqVKlajVueeZk7LzyDPod0Kbp/v6OP54xrbymP6jiOkyJ8bCgmJhgPevnthNfLIhhXqVqt1PU5M6fz/Wcfc/DJZ5JXOZWvu8zTuHFjpk+fTosWLahevToisvabnHJFVVm6dCnTp09PaieVdem5H7B+AmGquRcoAM7BAksvAMZjHs3/RvLdhgmAV2LOLl8C3VS1KJKmqs4XkX2B+7CYgwuAwZhQmFW88tgDFCxbRo/eF642z8EnncYnb77KQ/0vp+eFl7Jq5Sqeuv0Gli2xZpNKpjAuWL6MQX17s3DuHC68YwiNmrVgyrff8MIDg6mUV5mz+99WLnVyHCd5fGwwUikYx/PB6JcoLCykyxG5t1xcp04dAGbMmMGKYDLgVDzy8/Np0qRJ0e+1PpRZIFTVLuv9LSlEbZH8wXCsLd8t4VhTvh+AfVJWwAzw94w/eXHovZxz80BWFhSwsqCg6NqKguX8u2gh1WrWouMOO3Pmdbfy37sG8N6LIwBbAupy+NF8MPolatWtB8C7I5/l+88/5v7/fUzT1m0A6LTjLtSoXYeh113KAT1Pok2HTuVdTcdx1hEfG4pJpWAcz7hXXqDtFlvSpv0WaSl7pqlTp05SgoaTHeSWbnsDZdYfv1OwfBn3XHp+qWujHx/K6MeHMnDU/2jbcUsOPL4X+x51HDOnTaVGrVo0ataCm888gc222a5oRvz7z5OoVbde0YAfY7OttwXgz18mV9hB33GcYnxsMFItGEeZ/O03TP91CqdedWN5Vcdx0kJKBEIRyQc6YiFcFgI/qqrrlsuJth07ccNTI0ulX3/KUex1aA/27XEcTVsXx9DKr1KV1pu1B2DaTz/y7ScfcsFt9xRdr9doIxYvXMDMab+ViL318/99A0CDJs3SVRXHcVKIjw1GqgXjKONefp7K+fkJnVQcJ5tISiAUkYaYrd7x2DZ1MZaJyDOYXd+chDc7KaNmnbpsufNuCa9t1Lxl0bW5f83gzWeH0WG7zlSuUoVfvvs/Rj18Hzt3O6jEYNb1iGN59cmHueWsk+jR+0ILPvv9/zHygbtp12lrOmy/Y7nUy3Gc5PCxwUi1YBxjRUEB418bzXZ7dqVug4bpq4DjlAPrLRCKSBPgI8wzdyHwOfAX0BTYFjgd6Coiu0cdOZzMkVc5n8nffs3bzw1n6b//0rT1xhx9bl8OPvmMEvkat2zFgOfG8Nx9g3j2ntv5Z/58GjZtTrdjTqBH74uotBo7GsdxspNcHxtSLRjH+GrcOyxeOJ8uhx+T1vI7TnmQjIbwVkwYvBvoHw3TIiJ1sH2OL8KcOs5I9AFOenlx0owS5/UabUT/J54v072tNt2cfndnzzZUjuOUHR8bElNWwTjGuJefp1bd+uzQZb9yLqnjpJ5kBMJDgA9V9eL4C0E47CsinbFg0I7jOI5ToUhGMAa44oEnU1wix8kcyej3a2Px/9bEh0CtJL7DcRzHcRzHSTPJCISTgLW5lDUDfkriOxzHcRzHcZw0k4xAeA9wrIhsneiiiGwLHIPZGDqO4ziO4zgVlGRsCH8D3gY+F5Fh2NZ2s7C9hfcGTsL2Gp4qIntFb1TVD5L4XsdxHMdxHCeFJCMQjsP2NhbMi/j0yLXY7teHhiOevCS+13Ecx3Ecx0khyQiEN2ICoeM4juM4jpPFrLdAqKr9U1gOx3Ecx3EcJ0NUzLDyjuM4juM4TrmR1F7GjuNUfG4643gmjB9Hj94XcXyfywH49pMPGfvSc/w04Svmz55F/cZN2Gb3vel5QT/qNmxU4v6/Z/zJiHvu5LvPP2LRvHk0bNqc3Q7qzpFnXUC1GjUyUSXHcRwnxSQlEIpIS6AvtndxSyA/QTZV1XbJfI/jOOvHh2NGMfWnH0ql/2/EcJYt+Zejel9Ek1YbM3Parzw3ZBATxo/jrlfepXrNmgAsW7KEG07tyaqVK+h54WU0ataCX76bwHNDBjFz2q9cMjg3tzBzHMfZ0FhvgVBEugCvA9WAlVjImZWJsq7vdzhOsiSjHZs9/U8ev+Uafpv0PYvmzqVq9eq02rQ9h595HjvsvW+mqlRmFi9cwJO39afXFf25u995Ja6def0A6jZoWHTeaaddadamHdeddCQfvzmafXscB8Ckrz9n5rRfufbRZ9h2jy4AbLXL7vyzcAGjHx/K8qVLqFrdtYSO4zjZTjIawjuw8DEnA8+oamFqiuQ4qSF57di/1K7fgOMuupyGTZqx9N9/eOf5Z7j17JO49N5H2WX//5R3ldaJ4QNvodVm7dnzkCNKCYRRYTDGplttA8C8WX8Vpa1csQKAGrVql8hbs3YdtLAQ9TgDjuM4OUEyAuFWwLOq+nSqCuM4qSIV2rHWm7XnvFvuKnHvDnvvxzn77cx7L42o0ALhj199xvuvjGTQy2+X+Z4fvvgUgJabbFaUtvVue9Js400YPvAWzup/G42atWDKxG94bfhj7N/zZLchdBzHyRGS8TKeD8xLVUEcJ5VEtWPxlFU7loi8ypWpUbsOeZUrrj/WioIChl5/OYee1psWm2xapnuWLl7M47deR8t2m7HTfgcWpVepWo1bnnkZ1UL6HNKFE3fYjP69jqFzl/0449pb0lUFx3Ecp5xJ5q02BtuiznEqFKnSjsUoLCxECwtZNH8ebz//NDOn/sppV92YsvKmmlcee4CCZcvo0fvCMuVftXIlgy85l3mz/uKWZ18pIewWLF/GoL69WTh3DhfeMcQ0hN9+wwsPDKZSXmXO7n9buqrhOI7jlCPJCIRXAZ+KyP3AZar6b4rK5DjrTSq1YzGG33kTo58wb9pqNWrSd9CDbL3rniktd6r4e8afvDj0Xs65eSArCwpYWVBQdG1FwXL+XbSQajVrkZdnu0cWFhYy5IqL+PaTD7nqoWG0ab9Fic97d+SzfP/5x9z/v49p2roNAJ123IUatesw9LpLOaDnSbTp0Knc6uc4juOkh2R2KpkjIgcCnwEni8jPwMLEWbXiu2Q6OUEqtWMxDjnlTHb/z+EsmDObca+M5O5+59Hvnofp3LVbqoufNLP++J2C5cu459LzS10b/fhQRj8+lIGj/kfbjlsC8ND1l/PRG6Ppd8/DCYXc33+eRK269YqEwRibbb0tAH/+MtkFQsdxnBwgmbAznYCxQP2QtN1qsrofolMupFo7FqNh0+Y0bNocgM5du3HdST146o4bK6RA2LZjJ254amSp9OtPOYq9Du3Bvj2Oo2nrtgA8edsNvDvyGc6/7R523u+ghJ9Xr9FGLF64gJnTfqPZxm2L0n/+v28AaNCkWRpq4TiO45Q3ySwZ3wU0BK4DngJmqOqqlJTKcdaDVGvHVke7LbdmzLBHU1buVFKzTl223Hm3hNc2at6y6NqoR+7j1ScfYp8ePWm+cVt+nvBVUb46DRoWaQS7HnEsrz75MLecdRI9el9ogam//z9GPnA37TptTYftd0x7nRzHcZz0k4xAuCvwkqrenKrCOE4ypFo7lojCwkJ+/PoLmrbaOGXlzgTffDAWgPdeHMF7L44oca3L4cdwwW13A9C4ZSsGPDeG5+4bxLP33M4/8+fTsGlzuh1zAj16X0SlSr4duuM4Ti6QjEBYAExNUTkcJ2lSrR17bshAFi9cQIftd6Reo8YsmDObd0c+y5Rvv6HPwPvTXp9U8uKkGSXObxz+YpnvbbXp5vS727eocxzHyWWSEQjHATulqByOU26UVTu2SaetGPPUo4x//RWW/PMP9TbaiDbtt+Dm/46iw/be9R3HcZzcIRmB8DLgMxG5Arhd1Texciom66sd23GfA9hxnwPSUSTHcRzHqVAkYwB0DfAdcAswRUReFJHHExyPpaaoxYhIZRG5QkQmi8hyEflTRAbH5RERuUpE/hCRpSLygYhsm+CzthCRd0VkiYjMEJEbRSQv1WV2HMdxHMepqCSjIewV+b9tOBKhwOlJfE8ingT2AW4AJgGtgPiYIVcA1wKXhjwXA++IyJaq+heAiNQH3gF+AA4D2gGDMEH5mhSX2XEcx3Ecp0KSjEC4OgEwrYRg2McC26jqD6vJUw0TCAeo6n0h7RPMCeZ8ioW93kB14EhVXQS8LSJ1gP4ickdIcxzHcRzHyWmS2alkWioLsg6cBry3OmEwsBtQB3g+lqCq/4rIq8BBFAuEBwFvxQl+I4DbsX2aX01lwR3HcRzHcSoi2RhEbGfgZxG5T0QWBdu/l0SkeSRPB2AVMDnu3h/DtWi+SdEMqvo7sCQun+M4juM4Ts6SzJJxEcEJoxFQNdH1IGSliqaY/eL/AT2B2sAdwCgR2SV4O9cHFifYOWU+UENEqqhqQci3IMF3zKd4S74SiMhZwFkArVu3TroyjuM4juM4mSYpgVBEtgJuA7qyGmEQcypJieAZ+9pwHKaqc0M5ZgLvY44m76bwu0qhqg8DDwN07tzZQ+04juM4jpP1rPeSsYh0BD4G9gLexoS0b8P/c8P5OGB40qUsyXxgYkwYDIzHdk7ZIpKnVoLwMfWBJUE7GMtXN8F31A/XHMdxHMdxcp5k4xDmA7up6mEhbZSqHoh5ID+BCWjXJVfEUvyICZvxCFAY/p8E5AGbxuWJtxmcRJytoIi0AmrE5XMcx3Ecx8lZkhEIuwBjVHViJE3APHqBszEt201JfEcixgBbiUijSNpemHD6f+H8Y2ARcHRRwURqAN2BNyL3vQEcICK1I2nHAkuxJWjHcRzHcZycJxmBsBElvXhXYpo1AFR1JTAW2D+J70jEw9iS9Ksi0l1EjseWpd9R1fHhu5dhto1Xich5IrIv8AJW3yGRzxoKLAdeEpH9gsNIf+Auj0HoOI7jOM6GQjLOHvOAWpHzOUC8220BiW301htVXSQi+wD3YjEDC4BXgL5xWW/DBMArgYbAl0A3VZ0V+az5QVi8D4s5uAAYjAmFjuM4juM4GwTJCIS/AG0i518B3USksarOFpGa2HZwvyXxHQlR1SnAf9aSR7F9lm9ZS74fMO9kx3Ecx3GcDZJkloz/B3QNgh/Y8msD4BsReQGYCGwMPJpcER3HcRzHcZx0koxA+AhwOrYXMKr6GrZsWx3oATTGtoC7N8kyOo7jOI7jOGkkmb2MZwLPxaXdIyL3YQ4ns8OyreM4juM4jlOBSeUOIgCE7eJmrTWj4ziO4ziOUyFYb4Ew7AJSVVWXxKXvgzmTLAEeVtWUO5U4juM4juM4qSMZG8KBwDwRKQorIyI9sa3rLgAuBz4PO384juM4juM4FZRkBMK9gLGqujCSdj0Wy+9k4DKgHnBxEt/hOI7jOI7jpJlkBMJWwJTYiYhsArQHhqjq06o6ENsa7sDkiug4juM4juOkk2QEwjrYfsExdgcUeDOS9j3QMonvcBzHcRzHcdJMMgLhTKBt5Hw/YCm2Y0mMWtgex47jOI7jOE4FJZmwM58Ch4rIIcAy4CjgXVVdEcnTFpiexHc4juM4juM4aSYZDeGt4f5XgLeAKkT2DRaRasCewGfJFNBxHMdxHMdJL8nsVDJRRHYGTglJz6nqF5Es2wHvAc8mUT7HcRzHcRwnzSS1U4mqTgT6rebaJ8ARyXy+A0d2aJ7pIqwfvmuh46QdHx8Mb4ey89Zbb3H77bfzww8/MH/+fDbaaCN22203+vfvzxZbbFGUb+zYsVx77bV89dVXVK9enYMPPpiBAwfSpEmTEp/3/fffc+211/Lpp5+ycOFC2rRpw6mnnkqfPn2oXDnlm6E5acR/LcdxHMfZQJg3bx477LAD5557LhtttBG///47t912G7vssgsTJ05k44035sMPP2T//ffngAMO4MUXX2Tu3Llcc8017Lvvvnz11VdUrVoVgBkzZtClSxdatGjB3XffTaNGjXj33Xe57LLL+Pvvv7n99tszXFtnXXCB0HEcx3E2EI477jiOO+64Emk77bQTHTp0YOTIkVxyySXccMMNbLzxxrz88stFWr6OHTuy44478thjj3HuuecCMGbMGObMmcNHH33E5ptvDsA+++zDL7/8wrBhw1wgzDKScSpxHMdxHCfLadiwIUCR8Pfpp5/SrVu3Eku+nTt3pmHDhowaNaooraCgAIA6deqU+Lx69epRWFiY7mI7KcYFQsdxHMfZwFi1ahUFBQVMnjyZs88+m6ZNmxZpDvPy8qhSpUqpe6pWrcp3331XdH700UfTqFEjzj//fH777TcWLVrEqFGjGD58OJdcckm51cVJDb5k7DiO4zgbGDvvvDNffWX7SGy66aa89957NG7cGID27dvz6aeflsg/bdo0Zs6cSX5+flFakyZN+OSTTzjssMPYZJNNABAR+vfvz2WXXVZONXFShWsIHcdxHGcDY/jw4Xz66ac888wz1KlTh27dujF16lQALrroIj7//HOuueYaZs+ezaRJkzjppJOoVKkSlSoViw1///03Rx55JDVr1mTkyJGMHTuWa665hptvvtntB7OQ9dYQisjJwARV/XYNebYEtlfVYev7PY7jOI7jpJaOHTsCpik86KCDaNOmDbfddhtDhw7lhBNOYNKkSQwcOJBbbrkFEeHYY4/lP//5T4kl4zvuuIOpU6cybdo06tevD0CXLl1YtWoV1157LaeffjqNGjXKSP2cdScZDeGTwOFryXMY8EQS3+E4juM4ThqpV68em266KVOmTClKu+mmm5gzZw7ffvstM2fO5Nlnn2Xy5MnsscceRXkmTpzIpptuWiQMxthpp51YsWJFic9zKj7pXjLOAzxCseM4juNUUGbNmsWkSZNo165difSaNWuy1VZb0aRJE958800mTZpE7969i643bdqUKVOmMH/+/BL3ffaZ7VjbokWL9BfeSRnpdirZHJi/1lyO4ziO46SdI444gu23356tt96aOnXq8PPPPzN48GAqV65c5Bn8zTff8MYbb7D99tsDMH78eO68804uu+wydtttt6LP6t27N//973/Zf//9ufTSS2nYsCHjxo1j4MCBHHHEEbRq1SojdXTWj3USCEXk8bikw0WkTYKseUBrYE/gtfUrmuM4juM4qWSXXXbh+eefZ9CgQRQUFNCqVSu6dOnClVdeSZs2bQCoUqUKr7/+OnfccQfLly+nY8eODB06lFNPPbXUZ3344YfceOONXHTRRSxatIg2bdpw3XXXediZLER0HfZSFJFopEkFZA3ZFfgMOFFVf12/4lVsOnfurF9++WV6v0TW1MQVmBTv0fnSTzNT+nnlxZHtm6X087wdyk5Z92yNcuCBB/LWW29x9dVXc/PNN5e4Nn/+fC699FJefvllli5dyq677srgwYPZaqutyqM6ifHxwfB2cJyystqHZV2XjNtGPvBX4G7gngT5VgHzVfXfdfx8x3GclFCWPVujPPvss/zf//1fws9SVbp3787UqVMZMmQI9evXZ8CAAXTt2pUJEybQsmXL8qiS4zhO2lgngVBVp8X+F5EbgLHRNMdxnIpCWfZsjTF//nz69u3L4MGDOf7440t91ujRo/noo49477336Nq1KwC77rorbdu25Y477uDee+9Nb2Ucx3HSzHo7lajqDaksiOM4a+fIDs0zXYT1o4IsjcXv2Rrj8ssvZ8stt+S4445brUDYvHnzImEQoG7dunTv3p1XXnnFBULHcbKepL2MRSQPaA/Ux5xJSqGqHyT7PWv4/hbAT0BNoLaqLg7pAlwJnAM0Ar4ALlTVCXH3bwEMAXYFFgCPAjeo6qp0ldlZd1wQctaXVatWsWrVKqZNm8YVV1xRYs9WMA/KYcOGrXa5GOD7779nyy23LJXeqVMnhg0bxuLFi6lVq1Zayu84jlMeJCUQisi1QF+g7lqyJhQUU8SdwGJMIIxyBXAtcCkwCbgYeEdEtlTVvwBEpD7wDvADFkS7HTAIi894TRrL7DhOObGmPVsLCgo4++yz6devH+3bt1/tZ8ybN6/IAzNKgwYNAFtydoHQcZxsJpmt6y4DbgAWAsOBP4CVKSpXWcuwF3AgcCsmGMbSq2EC4QBVvS+kfQJMBc6nWNjrDVQHjlTVRcDbIlIH6C8id4Q0x3GymOHDh7No0SJ+/fVXBg4cSLdu3Rg/fjxt2rThjjvuYOnSpVx99dWZLqbjOE5GSUZDeCYwHdur+O8UlafMhKXqIcCN2FJvlN2AOsDzsQRV/VdEXgUOolggPAh4K07wGwHcDuwNvJqWwjuOU26sbs/Wq666iltuuYVHH32U5cuXs3z58qJ7li9fzoIFC6hduzZ5eXnUr1+/1G4MYJpDoNTWXY7jONlGMlvXtQJezoQwGOgNVAXuT3CtAxb6ZnJc+o/hWjTfpGgGVf0dWBKXz3GcHCC6Z+uvv/7KsmXLOPHEE6lfv37RATBw4EDq16/PxIkTAbMV/P7770t93g8//EDr1q19udhxnKwnGQ3hrCTvX29EpCFwExb0eoWUDkpaH1icwDFkPlBDRKqoakHItyDBV8wP1xJ991nAWQCtW7de7zo4jlP+xPZsPeGEE9h2220ZO3ZsqTxdu3blxBNP5PTTT2fTTTcF4NBDD+WJJ57g/fffZ++99wZg0aJFvPrqqwm9kh2nvMnWoPWQmcD1TmmSEeieB44QkaqqunytuVPLLcCnqvp6OX8vqvow8DDYTiXl/f2O45SNte3ZWq9ePbp06ZLw3o033rjEtUMPPZRdd92VE088kTvvvLMoMLWqctlll5VPhRzHcdJIMgLh9cAuwEgRuVBVf0tRmdaIiHQCTgP2EpF6IblG+FtXRFZhGr5aIpIXpyWsDywJ2kFCvkQe0vXDNcdxspSy7NlaVipVqsSYMWPo168f5557LsuWLWPXXXdl7NixtGrVKj0VcBzHKUfWaS/jEjeK/ArkA7EAcQtJvPyqqtpuvb4k8fceDoxaQ5bHgGeAd4EOqvpT5N7HgG1VdYdw/gEwXVWPi+RpBfwOHKqqa3Qq8b2M14DvVWp4OxgejzE9eH8wvB18ydgpKynbyzhKJSzMzO9r+aJUP6njga5xaQcClwP/wfZYngYsAo4GbgYQkRpAd8Jyb+AN4FIRqa2q/4S0Y4GlwPspLrfjOI7jOE6FJJmt69qksBzr8r1zgHHRNBGJleXDyE4ltwHXish8igNTV8JC1cQYClwIvCQitwObAP2BuzwGoeM4juM4GwoZ8RIuJ27DBMArgYbAl0A3VZ0Vy6Cq80VkX+A+LObgAmAwJhQ6juM4juNsEKRMIAzbwNVS1T9S9ZllRVWfBJ6MS1PMG/mWtdz7A7BPusrmOI7jOI5T0UkmMDUiUktEBonIX8Ac4LfItZ1F5HUR2T7ZQjqO4ziO4zjpY70FQhGpC3wC9AVmYLuARB1IJgJ7AseVvttxHMdxHMepKCSjIbwa6AT0UtXtgReiF1V1Ceapu28S3+E4juM4ToqY+9cMHr3paq48tjvHbbsJPTo0Z/afpS29Zv35O3deeCYn7diB47drx3UnH8WUif9XKt+i+XO5/6q+nLrrlhy3zSZccczBfPPhuPRXxEk5yQiERwJvqeqwNeSZBrRI4jscx3Ecx0kRM6dN5eM3X6Vm3bp03GHnhHn+mT+Pa44/nN8nT6L3DbfTd9CDAFx/ylH8+cvkonwrCpZz/SnH8M2H4zip3zVcOuRRGjZtzoBzTua7zz4uj+o4KSQZp5KWwItrybOYxDuBOI7jrDfZGoTXA/A6mWaLHXfh8Y++BeCdF/7L/31UOuTuWyOGsWDu3wx5+iWatm4DwFa77MG53XZhxJCB9Lv7IQA+fnMMv//8Izc8NZItd94NgO327MrFh+3H8IE3c/sL5b67rJMEyWgI/wEaryVPW8zZxHEcx3GcDFOp0tpf+z9P+JpmG7ctEgYBqtWoQccdduarcW+zauXKkO8rqlSrViQMAogI2+6+F1MmTmDurOycuG2oJCMQfgEcIiK1E10UkWbYziHjk/gOx3Ecx3HKkUp5laicX6VUen6VKhQsW8Zfv08N+fKoXDm/VL7KVaoC8PvPP5W65lRckhEI78ECPr8uIh2jF8L5C0A14N4kvsNxHMdxnHKkedt2zJz2K//Mn1eUVlhYyOSJEwBYvHABAC3atmPJ4n9K2BWCaQ4t3/xyKa+TGtZbIFTVt4AbgN2B77AdQRCROeF8N+BKVXXLUsdxHMfJEvbveTJaqNx7xUX89ftU5s+exWM3X8PsP38HQMKy856HHEGd+g0YcsVFTPvpRxbNn8uLQ+/lhy8/Bcq2PO1UHJL6tVT1BiyszGhgPrAKUOB1YD9VvTPpEjqO4ziOU240bbUxF915H79+/y3n7b8bZ+y1HT9P+IpDTjkTgPobmftAzTp1uXTIYyyaP4+LD9uXU3fdivdeGsGx518CQL2NmmSsDs66k/TWdao6FhibgrI4juM4jlMB2PWAg9lpvwOZOfVXKufn07R1Gx7qfwWNmjVno+Yti/Jt0XlnHnj7E2ZO+43CwlU0b9OOVx57gCrVqtGu09YZrIGzrqRsL2PHcZyKxNy/ZjDqkfv55btvmfrT9xQsW8aD73xG45ativI8N2Qgz99/V8L786tUZcS3RbtxUlhYyMuP3M//nhvOgjl/07xtO44+ty+7HnBw2uviOJkgLy+Plu02A2DerL/4+I3RHHbaOaXyiQjN22wCwNJ//+WdF55h70OPolqNGuVaXic5fIE/RzkQ20fwmrj0CeFaLaAOcCgwJcH9vwOnAK2B6sDm4bP+TU9xHSfllCUA735HH8+AEa+WOK5/4jnyKlem8z77l8j77D138Nx9gzjohFO5+uGn2Xyb7RnU5yy+ev/d8qhOyvCxwUimHfqHexMd1dJV4BTyyZtj+OTNMfzy/UQAvv7wPT55cwzff/4JACtXrOCJAdfz2TtvMPHT8bw+/DEuO+ogWm26Od1PPbvEZz096FY+eXMM3332Me+88F8u63EAeZUrc8LFV5Z7vZzkKLOGUEQKgUJgC1X9OZxrGW5VVXVNZDnyLFB6gyGYjG0uvSXwX2Al5hW0FzYIxoJK/gvsB6wAbsIG/i+A68NnPJe+ojtOyihLAN6GTZvTsGnzEmnjXhnJqpUr6Xr40UVpC+fOYfTjQznizPM47HTTkGy1y+789ftUnh50KzvsnR07dPrYYCTbDmdgQmOUf0PaoakvbsoZ2OesEueP3GDCW6cdd+XG4S8iIsyc+hsfjhnFv4sW0bBpM/bp0ZMeZ19AfpWS4WgWzv2bxwdcz6J5c6jToBE773cgx17Qj9r16pdbfZzUsC6C2geYALgk7typQMwH+gKDgePjrt0O5AFvAPVC2s7ApsBA4I6Q9hE2ML4FxHQkXYF5Id8SwBcCnIrO+no4jnv5eeo12oht9+hSlPbN+HGsXFHAXof2KJF3r+5Hcv/VFzPrz99p0rJ1MsVNOz42GKloh5bhiDIcEyBPSXmJU8+Lk2as8Xpe5cpc9dCadqUt5rxbB6eiSE4FoMwCoap2WdO5UzG4HJvdHkfpwe5TYFeKBzqwQW1LYBTFg11B+Fsn7v56mIrYZwFOrjJn5nS+/+xjDj75TPIqFw+Pf0z5ifwqVWm2cdsS+Vtt1h6AP6f8XOEFQh8bjFS0QyKeApoAB6SqoI5TzrgNYQ4xHhgG3L+a63lA6djzUBX4BVgWzvcDNsMGzh+wDanfwyKR9wZqpq7IjlOh+GD0SxQWFtLliKNLpC9esICadeogIiXSa9WtB8A/IVBvRcXHBiNV7RDPH1iojRNwT00ne1lvgVBEfhWRC9eS5zwR+XV9v8MpOwXA2UA/oP1q8rQHvsLsf2L8A3yPzexjMeWrYQNnIdAJqI0FmzwEuC/VBXecCsS4V16g7RZb0qb9FpkuSsrwscFIZTvE8zTWJtmwXOw4qyMZDWEbSmrWE1EP2DiJ73DKyB3AUuDqNeS5EJiOzeSnA9OAU7FZPhR3hmXAscBszC7mfeBOzGD8vFQX3HEqCJO//Ybpv06hy+HHlLpWq25d/l20CNWSi6KxLbxqB01hRcTHBiOV7RDPMGA7wKPuOdlMupeMa1NsduKkid+BWzCvv+XAgnAQOV8F7IEtlYzE7GLaAAuxWW0VoEG45zFgHLbdzImYh10/YBAwlMTeeU7FJVF4jXex37YdFjqkHXAO9qKPJ1fCjKyNcS8/T+X8fPY85IhS11pt2p4VBcv56/epJdL/mPIzAC033bw8irjO+NhgpLodonwOTMK1g072s04CoYi0jh0hqV40LXK0FZG9gR6ALxmnmV+xmfuJQP3IAeYZVx+YGM7PxV7632GD5NvADMyTLj/kmRjuaRf3PTuFvz+mvAbpIxlhaBpwGKbirg40AvbGXobZwurCawwF5mLt8ia2EfloYBeKtSFQHGbkA+xl+joWcmMQcFraSl3+rCgoYPxro9luz67UbdCw1PXt9uxK5fx8Pnz1pRLpH7z6Iq0361BhHUp8bDBS3Q5Rngrp8Q4qjpNtrKv961RKOpJdFI7VIcDF6/gdzjqyLYn3DuyKDYCnY2ETYlTF7H/ABsF3sCWPGE0xW5kpcfd9Fv62SLrE5cOahKHFmDC0CRZG43oslMa3WEBaQp5GwM2YtmAR8AhwMPAicGQay54K1hRe4wFgo8j53pjmb2/geYqFvWwPM/LJm2MASgTgrVu/IXUaNKTTTrsW5ftq3DssXjg/4XIxQN2GjTik11m89PB9VKtZi0222IqP3hjNd59+xBUPPJn2eqwv2+JjA6S+HWIUACOAgyj5PDlONrKuAuEwTCAU4GTs/TkhQb5VmALiXVX9XzIFdNZOPaDLaq5tHLn2J/AgsBs24H0JDMAEm+Mi9/QC7gL+g9nbtA55bwJ2AHZPXdHTRiqEoU7YElmUg4G2wBNUfIFwTeE1Er28dgx/p0fSsj3MyNoC8MYY9/Lz1Kpbnx267Lfazzq+zxVUr1GT14Y/yoK/beu6SwY/ROeu3dJT+BRQDx8bIPXtEGMMNjny5WInF1gngVBVe8X+F5GTgVGqemOqC+Wkh3xsJv8Q5jnXDriO0ireNlg8rv6YFm0O0Ao4C3sJZEOsolQIQ4moDNSl4oeWiIXXWBebrtg+Hh0jadEwIw9iAsDnZE+YkbUF4I1RFi1fXl4eR53Th6PO6ZNcoSogG9LYsCbK2g4xnsLsCg8pl9I5TnpZ7/eaqmb7s5/zxGtvmmBLH2VhC0xblo2kShiKEduzcQ7wMPAzJhBVVMoSXiOef4A+WP0Pj6THwoz0oHgJDcyOsKKHGXFWz4Y6NsSTTDsAvJLCsjhOpqnoig7HWSdSKQzFuAxzogCzLxyBxV6rqJQlvEaUlZgmdTpmMxgdFOLDjMQ0hDeGfA+mpsiO4zhOhklaIBSRHbHdelpgZhfxqKqenuz3OE5ZSKUwFKMP0BP4C9M8Ho+FpaiIy0Sx8BqPYuE0lkeuxcJr1MZ2ZIDiYLrvAK9ROo5aLMzIFIo9S/fCls3PwpaNt0ltFRzHcZwMsN4CodgeTk9iTlpCsbNJDI2ku0DopJ1UC0MxohvZH4IZoPejYgqE0fAa8QwMxzeY1yWYQPccJuAm0nqWJcyIC4SO4zjZTzJ2gOcDJ2ErSZ0x4e9uzEHrKmwlbgQW2SNliMjRIjJaRKaLyGIR+UpESjmAiciZIjJZRJaFPKXedyLSQkRGicg/IjJHRO4TkYocRcNZA+sSawyKhaF1XQLujGnMKiLbYuE14g+wdhlLcXiNSzDh+QkSL5VDyTAjUSp6mBHHcRxn3UhmyfgU4KeY53HY9H2Bqn4KfCoib2EOaW9j75xUcTHwGxZVZA4WAeEZEWmkqkNCWY7Dws31x2ziTwXGiMiOqvpdyJOPhVcrwFYE62ERFeqRWMHiVHC2peyxxmLC0FOsXhhKRCHWoeI1ZhWFepQtvMbtWGc/DfMi/jSSbyOK69eL7Awz4jiO46wbyQiEHbD3acLPU9VvRGQMFvg9lQJhd1WdEzl/T0SaY4LikJDWH3hKVW8CEJH3sa0mr6BY2DsK8yPYVFV/C/lWACNE5AZVnZzCMjvlQD1SKwz1x2KM7Y5pyv7CbOo+B55JVaEzxBvh7+PhiHIKZgsCuR9mxHEcxzGSdSpZGPn/X0pv9TiZ4g0OUkKcMBjjGywyBiKyCRZn+KLIPYUi8gIlw0kdBHwREwYDL2MawwND2Z0cpKzC0PaYDcQIrKM3xezlPiT7NGPx4TXGrcO9uRRmxHEcx0lMMgLhdEqaEP2KrSJF2QwTFNPNrlh4ODDNJdh+41F+BBqIyEaq+nfI90M0g6oWiMgvkc9wcoD1FYYODYfjOI7j5DrJCISfU1IAfAO4VESuBV7CVugOw3b3SRvBWeRwincci/kRLIjLOj9y/e/wNz5PLF/9BOmx7zsLWzGjdeuKuaG94ziOs2FxZIfmmS7C+qMVfRPMDYNkTIBeBPJEpG04vwOYBtyA7XE8BBO4rkimgGtCRNpg5lyvqOqT6fqeKKr6sKp2VtXOG23k25k7juM4jpP9JLN13cuYzV3sfJ6IbAecidnlTwWGqerM5IqYGBFpgGklpwEnRC7FNIF1KakBrB93fX7IE0991m3XM8dxHMdxnKxmvQRCEWkN7IiZZ32hqn8AqOpCLORbWgmxAscAVYBDVHVJ5HLMdrADJiwSOZ8X7Adj+UrYCopIFSxu4tB0lNtxHMdxHKciss5LxiIyEHMgeR54AfhNRO5MdcHW8P2Vw/duBhyoqrOj11X1V8zB5OjIPZXC+RuRrG8AO4rIxpG0Q7Ht995MT+kdx3Ecx3EqHuukIQwBny/GNIOTsN1J2gMXi8jXqvps6otYigewOLkXAQ1FpGHk2jequhwLm/a0iEzFtqg9BRMgj4/kHYmFUnspOMLUBQYDz3gMQsdxHMdxNiRE18G7R0TeBfYEDlDVsSFtP0zb9r6q7peWUpYsw1QsznAi2qrq1JDvTOByLI7u98Clqvpu3Ge1BO4D9sO2ux0R8kWXoFdL586d9csvv1yPWqwDImvPUxFJtdeYt4Ph7WB4OxjeDoa3Q/a2AbiXcfmy2o6yrgLh38A4VT06Ln0k0EVVG613EbMQFwjXgA/4hreD4e1geDsY3g6GC4SGC4TlyWo7yrraENandMBnQlq9dfwsx3Ecx3EcpwKwrgJhJWBFgvQVrEHqdBzHcRzHcSou6xOY2nW7juM4jpOF/AlcgO33WgPT5ExNkO834Chs6a8m0BVIZCA1B9smbCOgOrAz8FaKy+yUD+sjEPYXkVXRA7gOID49HCtTW2THcRzHcdaHKVjMuPqYh2gi5gJ7AN8BD2HelmBC4Y+RfMuBfbA4bXdge9a2Ag6h7HvGOxWH9QlMva5Lw76U7DiO4zgVgL2AWeH/R4H/JcjzYMjzAbbtGJjgtwlwPSZQggUEngiMBbqEtAOBbYDLgM9TW3QnzayThlBVK63Pka7CO47jrI6yLI31D+mJjmpxeQuBAUCbcG0bbEN3x8kmyvJC/hQL3NsuklYT0yiOAVZG8lWnWBgEe3b2B74ApidZVqd8cWHNcZycpCxLY2cAn8Qd72BLJ4fG5b0WEyDPxwKv7oJtf/R6isvtOJkmD9sXNp6qwFLgl0i+/NXkA1tydrKH9drL2HEcp6JTlqWxluGIMhzTgJwSSZuNbdJ+BdAvpHXFhM4rsK2THCdXaA+8jdkSxrYCK6R4CXheJN8izK6wY+T+T+LyOdmBawgdx8lJ1ndwewpoAhwQSXsLKABOjMt7ImZD9dt6fpfjVER6YwLgyZg2cCZwIcX9PPZsHQ80wiZPEzGP41sx28NoPic78N/LcRwn8AdmIH8CJZdPvseWwTaNy98p/P0h/UVznHJjE+C/wFdYn2+Oaf36huvNwt96mGfxHGBrLPTM45hpRTSfkx24QOg4jhN4GtOMnBKXPg97+cWHTGgQue44uUQPzCnkB8w04itgMRZWpnUk356YFvFnbOn4Z8yusDqwQzmW10keFwgdx3ECw4DtMG2H42zo5GG2ge2AGcBzwDkJ8gnmldwBWAI8ApyEeSZnOx999BH7778/jRs3pnbt2my//fY8/vjjmS5WWnCnEsdxHMxgfhJwd4Jr9YEF2DZNUS1hTDPYIP4Gx6nAjAx/vwp/38CWezcC9sb2or0s/F8HM5kYgJlIXBL3WVdimsBGmCbxTkxDOCB9xS83vv32W/bbbz922WUXHnnkEWrUqMHIkSM5/fTTWb58Oeeck0g8zl5cIHQcx8GcSfIxQ/l4OmG7MvxCSTvCmO3gFuktmuOklKPjzs8Nf/fGdhgRYDLwDDYRaoltT3cVpcPRzAL6YJ74jYEjgBvIjUnSiBEjWLVqFa+++iq1atUCoFu3bnz77bcMGzbMBULHcZxcowDbnusgTEsSz4GYsPhfbKeGGE8DWwJt011Ax0khupbrlbEA1GUhNxdPjYKCAvLz86levXqJ9Lp16zJ//vwMlSp9uA2h4zg5y8hwRJfGRgLvx+Ubgy3/xjuTxGgMXIwtg92FaVHOAd4jN5bGHMcpTa9evQC48MILmTFjBgsWLOCRRx7h3XffpW/fvmu+OQsR1bXNFZzV0blzZ/3yyy/T+yWSpVtBp7pfeTsY3g5GGdthdbliS2MxDgPGY/HWEu3QALAKE/4eAf7CgvJeBxxVppIEvD8Y3g5GKtshW9sAUt8fUsgXX3zBEUccwfTpthFffn4+Dz74IKeffnqGS7berLaj+JKx4zg5S1lfM6+UIU8ecE04HMfJfSZPnkyPHj3o1KkTQ4cOpXr16rzyyiv07t2batWqccIJJ2S6iCnFNYRJ4BrCNeAaAMPbwfB2MLwdDG8HwzWERgWVQ44++mi+/vprJk2aRH5+8a7NJ5xwAm+99RazZ8+mUqWss7xbbUfJupo4juM4juOkm4kTJ7LNNtuUEAYBdtppJ+bOncvs2bMzVLL04AKh4ziO4zhOHE2bNmXChAkUFBSUSP/ss8+oVq0aDRrkQnCdYlwgdBzHcRzHieP888/nt99+o3v37rzyyiv873//Y+ONN+bZZ59l6623pkqV1bmgZScuEDqO4ziO48Rx1FFH8frrr7N8+XLOOOMMDj30UP744w8AKlfOPZ9cFwgdx3Ecx3EScNBBBzFu3Dh+/vlnqlWrRp06dQCQbHbiWQ0uEDqO4ziO46yBPn36sHTpUh588MFMFyVtuEDoOI7jOI6zGsaPH89///tftt9+e4477rhMFydtuEDoOI7jOI6TgIKCAk466SQAnnzyycwWJs24QOg4juM4jpOAAQMGMH36dC655BLat2+f6eKkldxzk3Ecx3Ecx0mS33//nZtvvpn69evTp08fFixYUHRt1apVLFiwgNq1a5OXl5e5QqaQDV5DKCJbiMi7IrJERGaIyI0ikhu/ruM4juM468Unn3zCypUr+fvvv2nevDn169enfv36AHz66afUr1+fCRMmZLaQKWSD1hCKSH3gHeAH4DCgHTAIE5R9D3vHcRzH2UCpWbPmWvMsX768HEpSPmzQAiHQG6gOHKmqi4C3RaQO0F9E7ghpjuM4juNsYOyxxx6MHTu2VHrXrl3ZcsstGTJkCFtvvXUGSpYeNnSB8CDgrTjBbwRwO7A38GpGSuU4juM4zjrz0k8zU/uBzRI7kqysUp15zdrzv+n/AP8k/TVHtm+W9Gcky4YuEHYA3osmqOrvIrIkXHOB0HEcx3GyhCM7NC+X7+nx9Rep/S7V1H3WerKhC4T1gQUJ0ueHa6UQkbOAs8LpYhH5KT1FKxcaAXPS8snZta2Pt4Ph7WB4OxjeDoa3g+HtEGmDW8KRMsqvDd5U1QMTXdjQBcJ1RlUfBh7OdDlSgYh8qaqdM12OTOPtYHg7GN4OhreD4e1geDvkfhts6GFn5gN1E6TXD9ccx3Ecx3Fyng1dIJyE2QoWISKtgBrhmuM4juM4Ts6zoQuEbwAHiEjtSNqxwFLg/cwUqVzJiaXvFODtYHg7GN4OhreD4e1geDvkeBuIVgDPlkwRAlP/AHyHhZrZBLgLuFtVPTC14ziO4zgbBBu0QAi2dR1wH7Ar5nH8KNBfVVdlslyO4ziO4zjlxQYvEDqO4ziO42zobOg2hDmFiGzwv6eIVBGRRJ7jzgaIiFQWkSqZLkemCc5yjuM4q2WDFyByBRHZHjhjQxaGwov/ReBOEWmQ6fJkEhHZVkRyNl5WWRCRasArQO/w/waJiLQAPhCRnDaIXxsiUl1Ejs90OTKNiGwjItdmuhyZRkRqiMi+G/LYEI8LhLnDmcBQ4DgRqZPpwmQCVS0AFmL7UF8pIg0zXKSMICJVsSD6n4rIjpkuTwZZAbQELgBO2IAH/lXYNpwHicjgTBcmg9wBDBaRCzJdkEwhIvnAScC1InJ9psuTYe4G3gb+E8bMDR7fqSRHUNVzRCQPuAeoJCLPqOqCDBer3BCRSqpaqKonish9wGFYO9yqqnMzXb7yRFWXi8hl2PP9uogcoqqfZbpc5YmI5KnqKhHZAXgLuAbrD0+r6tIMF69cUdW/RORG4F9swqiqenGmy5UBbgfaAGeF/nF3ZotT/qjqChEZCuQB54hIFVW9OtPlyhC9gbbAA8B5IvKaqi7LcJkyimsIc4DY7EZVzwIeB24CeopIvUyWq5yJbgR5GzANOBi4PIQX2qBQ1e+BC4HxwKhgUrDBoaorgaOx2KJnY5rCDUYbEDRCqOoc4DPgW+AiEbkpowUrZ0Skqqr+CZwOTAVOFJE+GS1UBhCRfFWdAgwBngTOFJErM1uq8if0h0JV7QZ8g4Wb6y4i1TNctIziGsIsJ8x0lwfh7xpgMbb13o12Wf6rqosyWcZ0I6byWBX+fxGohmkCqgJ9Q/ptqjovY4UsJyKasWrAEcAsoCnwoogcrapfZraE6SeuPzyLLZnWxJaPbwYKgwY9p7UBoR1WBBOSD4G/sF2YFmImFTVU9ZKMFrIciI2R4XQH4Cdgd6CviBSo6gOZK135Edoh1h+OBzoCVYBbQl+5NbMlLB+i/UFEjsbMKQ7A3pkqIq9G+ssGhWsIswwRiWrCCC//GsDH2GA3GTgfmIDNek7MVZvCsESOhthJQeuxJ2Y/1xUb8F4EemIvwJzXFIb+UBOYCHQD5mC2pZWA10Rkp0yWL53EvOwj/eFuoAvwBKYt3gOYgQ38x+e6plBVNTwjTwDLgQtUdU9snHgEOElEBmWyjOVBZHLwPHAvJgQ9gi2bXiYiF2aweOVG5F3xOWZnPR44D/gSuEREbshk+cqLuP5wG7Ax1h8qYZrTQ3J9bFgtqupHFh1A1QRpR2Oz/50iaXnYksAS4BygXqbLnsI2qAHUTJD+LPBmfBuFdliK2RA1yHT5U9wWeQnSLgN+BdpH0nYExmEaw86ZLneK26BafDtgL/13MROKypH0fGzpdDa2fFgt0+VPYTt0BGrHpQn2wn8kLr0FNlFYBdyY6bKnqT0qRf4/EPgH+A/F8Xebhz4yBTg30+VNQ/0lQdopYQzYOpLWPvSFucBVmS53ebRH6A8LgEMiaU2wLWv/Anrk0thQ1sM1hFlECCPycIKQKvWBRsD3IV9ltVnQGZjm8HrgmFwIxSIilYHXgEGRtLzwbz2gjhYvB8RsK3sBP2OC86250A4AIrIHiUOq1MWWy/+MJajqF5igWAA8nys2hSHU0JvApZE0ARQb4Guq2RHG7IZWYEvpKzGj8rNyQRsQPOo/wvp4LK0SUCt2GtLygwPWdGyCtAjTkj1WzkVOCyJSW0TuBFDVwsil2pgGaKKqaugLM4ATgULgOhE5r/xLnB5EZHfgoTBeRsnDJkWLQ75KqvoTMBizrbxZRC4vz7KmExGpJSLXBfOI6C4cG2ETyU9DvjxVnQUcAszEno2Dg0Z1g8EFwuyiG/CLlraF+xJbGjwlPOArg1C4EhgDNMZmgF3KtbTpYwDQD4pe8rFtBh8HOojIOVDkbRsbEP/GhKTtyB3b2TOAzbW0Ldx07CXXMWpioKqfA29g9pVfim3bmO3kA//Ftp+MveA0CH7PAjuLyIFQ1B8E0wzMBrYFTgVywZB8AdBVVR8XkarhBVioqv8AzwCnicheoV1ipha/AZOwpcM6khuB7U/CxoD4SdJk7HfeBYr6QlVVnQlcha069I2NHTlAG2BJbDIUYTpQB+gQzitFhML7gHnAABE5u9xKml4uw+oavxXtd8AyzHYQteX0/PC8PApsAjwE7FeOZc08mVZR+rHuBzaw3QxsGs4bAu8Bn2C2c9GlknMxr+OLiSydZeNB3BIIcCvwC2GZDGgHPI8JyKdG8tXHXopdgcaZrkca2qUGcCxhORyb/c4AXsK0ptGlktuxMAs35WB/GAC8HKsX9vL/CjMa3yeSrxkmRG4NtMp0PVLcJpUwB5J3gLohrUlog7nAXpG8W2Bx2A6keBm1UnmWNw31rwNUCf+fFknfCLMn/hzYPe6ec8LYORLYJNN1SFU/CH9rYCtENSPXnsUEv60iaQJcjQVyP5MEpijZeGCTn2rh/54E0ymgVfjNXwe2i7vndGBEGEval2d5M31kvAB+lPGHKinknYlpBJ4C2oW0TTFV9yeYU0k9zIj+U2BA5N6sFQLi2qAScBa2FPwxtlQMsBumFf0jtM/1wGhgPtAy03VIU1vcgmkEzwUahrR9MZupUcDO4UW5JaYN6peD/aEaZhA+GdMU54f0IzFTim9CXzgNmzTMIkcmB/Evb8yDdF6oZ62Qtmt4DgrDc/EUFoLm84jwUMrmLJsOSk58/hPq+ngk7TBssvg+cHBIaw08BvQnCJK5dETa4RGgekjbCfggjImHYjaEncPYcPXq+lW2HdHyY6YBC7Fg1PVCWhcsNucrBFvC0B+GA3dmuvwZabNMF8CPMvxIoWPHvQCvxpZ7nqZYU9gO0wzMCoPA38DX2fzSj2+D8H+nWBpwAhZG4jOKhcItMJuyyZjB+KdEjKiz/Yj0h+gLcBjmOHM+xZqhfTBN4Z/YZOH3IATkWn9oHf7WxjyIJweBJyYUdsGWgRaEdvge2CbTdUhRO8S0oVXDiz42ITg8vABHRgSBxsBF4Vn5BNOSxtoo2zWD8UJxQ+DyMBYOi6Qfg9kgr8BWF37DBKOcGB+iY0Pk/xOC4PNE5PfeEVs1WRmeizlYZIpYf8r2yUEiJ7P7gR+xzRvqh/T9wngwK4yPP4X+sFV5l7kiHLFlAqeCErMFDKFEbgW+VdXHwrWrgZOBL4D+qjolOEy0BrbHOvZoNfuImE1h1hGLrRf+fxjYBrhbVZ8NgXePBa7D6ttNI3EXRaQ2UKiq/2ag6Ckn0h9qYBrSaao6KlwbhrXFJcBwVV0oIo0wj7lmmAPBveH+ojbNNuL6w33AVpin7Lvh974ca4ePgTPUYq/lY8tndYB/NAd28YnEnKyNOdbMw4SfF4Lt7CGYh/27wCmqGnMkqBl9HrJ5bIBS/eFMzMP+PWyV5GzseXhTVU8KedpjdmXdsAnCC6r6cwaKnlLi3hXXYMLNy9hKQU9MS/g8cKaaLSkisi82NhQCz2X7uwIS9offVPWd8EwMBvbHnpf+qjpfRDpiY8jemGD4nJpN5YZHpiVSP1Z/ULyUUxubxXwADCTM+MO1q7EH/2nC8nGCz8lq1X+kHs9gg/0pRGy/MMeCE0M7fEqxTWFO1Hs1/WEiFkbm8rj+MByLOXc+QVuUw/3hOUwDfD7QNpJeG7OxnYwJRDmh9VhNG9QA/g+zBdyXSMglTIN+OKYBeoFizXF0pSGr2ySuLiMwjV9ReCmgAXAFtloyLBNlLM92iLwrPsWczmL2lJUxTeESTFNYfTWfk9VjQ4L+MDmMEzHNeWXMeeYnTFNYL9NlrkhHxgvgx1p+IBN2XsXsO9pEHvxobLXY8vFThOXjXDuwWd0f4QUXXS6MtUdMKPwBWxaolekyp6kdqgBjMQ1IB4qXhfIjeWLLx+cRlkZy7cA0Hn9gXoJFNnBxL8abQ18YSQ4sk6+mHc7Glj63iaRFX4qVMdu5eZimsEamy5ymdrgXM43Yg+LlwNgKWFQofCzTZU1jG1TG7KfHYeZD8cumMaHwX0xbWCqWa64coT/8gdmUx/pDbKyMCoWDCBMlPzwOYTbQAgs4O0xVp6pqYdhmaGUspIqq3oJphg7FDMpzkXbY9mMfqi1rxMKpKNim7Zj33CBswGuUkVKmny2AzYD7VHVSrC20eFkUVT0Z0wgNwewIc5GOWNiIjzTEm9Mw2of//8F2IXgT28B+o0wUshzYAihQ1f+Dou3qCiO7tqxU1VewycEqrM1yChHZCHOcehj4VFXnQ9FOLZXUwnQ9DNwBnCoi92eutGmlIeYg8gK2TBodJ1FbBh6BmZqcju11nnOISEtsx6p7VPXjSH+ILof3Ad7CTK4uj7bThkyuxGPLZephe9HOgZL2PkEorKWqi1X1FhGZjgmGuchy7IXWFJgbBnuJCQEisjdmLP8UZhOUq/s3t8Z2WPgGStiRxYTCWH84WUR+wTzocobIb14N6w+Vo/0gJhyKyLaqOkFErsPCTvyduVKnlb+BhiLSXFVnRJ6LwmBnehjwlqo+i02YYrEaC9f0oVlGLcyu+HGN2L7F2gFAVeeJyKOYM8nrmSlm2mmCTZx/jk0Kor9z6A9VVfW/IjIH0xjnItWwSfNcKB4zIoqU2N9LsPfKY9HJ5IaMawgrGAlmKtOxqPL7Q5EQmBfyVsY2Jj81XHsyCAd5ZClrCI77Lbb0c4yIVI4+5CKyLRZcds+gEclVYRBsSXw55igSm/XmhbbIBy4UkcPCtRuimuRsJPo8xIJOh9OPgc2Bg0Ld8yL5dgfuDELhP7kgDK7hufgRG8dPEZF6UEJTugdwHLZXaxE5JgyC7b4zA+gkkV1nIpPF40XkP0FTdK/mgAPJavgT6w8niEiDqKY4jAEnYfv0VlLVt7J9bFgDyzGTmU3DucSNHReLyEGqukJVL1XVXzJTzIqHC4QVhNgLLX6mEl5mtwNnikifkBbzDt0YCx/QIe6ebPYejWl4GojIpuH/Sqr6JeZJfA22S0mLcFtbbAuyTTCbkJxgDUL9Qizg9CkiciSU+L3bY0GGm0dv0Cz1GIwJuuH/KtiSGABhGfRp4AkROTRyT0vMlrQpthtJ1hNWBQrFtp1rIyLbinmPo6ovYmFUrgXOE5Hm4Z4dQlrM6STrWd0zobYN30tYjMl94yYHDTEzmu4iUiUXhOE1tMM8zLa4O3B0RCjMwzSHZ2CBtwsj92Tl2ABrbIc/sF1GrhCRQ9R27Cl6rwB7AUeISBVfKi6Jh52pAEjJUCJ9sDAAlYDrVXVOeMldiwWkfhyLNdgQGwABds7mBxtKhQoYgmk3OmLecg9g4RNWADdgQuEPmGYALOr8fjE7qmwn0h+qA0dh/WEe8ETQCHbFQhDVxhxIXsa2YeuLeZbumq2Tghhx/eF2zDh8c6zv36Gq/ycibbDJ0tHYMmAhto/zVkAXVf02E2VPJVIytMxoTEveEbN/eklVnwj5ngQOxsaNudiy2Sxgt2BKkNXLxHH94XAspuK/wFhVnSG2Vd3zWOiQmzHHq/pYRIL9sR1aJmWi7KlESoaWOQsLNv8d8I2qjgt53sBiUr6OmQlsjgnFlYBdsv1dAaX6w2GYaVUVVX0kpLUH7sR++36YoNwEe2cegK0m5YwCIVW4QJhhYgN1GPDHY+ruaZgNRCXgYrUYSi2wZcLLMZuZv7DwAseGAT9r48pFEZFnsMHsDuzl/yFmJ/UQ8Gio6wFYDLHm2GD4gqpOzlCRU0qcAPABZidXFRN2FgEHquq/ItIFG9yOCtdnYMtFB2d7f4jaBIrICGAHbCL0GaYNex+4XVXHhiXCXljMvZpYf3ggF17+McLLfzymHb4C86gfgfWNQap6T8h3OPbyb4Dt4POU5kZcuSJhVkSGYxqeBdiS4DjgQVUdIyJ1sLBcJ2DPxEzM3KZnLkwW494VH2ET5JnY710TGKiqw0Pee7CA7FthfeFH4JhsHxsgYX/YCQuw3QgLrn1ZmDBuhe3edBY2eVhEiMmYC/0hLWgFcHXe0A9sNv8uNotpHNJewYSAP4D9I3kbYhqClhQL9DkRUgPbU3QCYa/RcL4Ke8nPwJaGY/tS5kSdV9MO1TEB4D1sabwyJhwWYvEHY2EUGmChiPbHNAWlQhJl84FNfiZiGk+wwX0VFoB8IuZBXbRna/ib1XHUErRBJUzT8R7QNKQ9jWn/3sYmj73XcH/OtAc2KfyVsB8zcFd4JsYTtqIL6Z2xCeNOwEaZLneK2yD2rngHaBHSxmKCzhRK7t9cL7wrGuXauyLU5eHQ//cM54NCf/gM2DaSbztMmbI30CTT5a7IR8YL4IeCLfV8DHQM5yOxbXT2xmaCv4WXfqlgomT5llOReuQFAeDGcH4eFiLjUGy2/xMmHJ8bEwpz9cBs4N4HNgvnz4f+cG74+xWriZ2VQ/2hKhY254ZwfiGmPT8Isxf9OzwbB5BDgZYTtENtTCt4dDh/EpscbYppTv/EzAnOynRZ09wOXTBt4P7h/CrMZOQKTEv2f8DhmS5nObTD/gQTkXD+UhgXTwpjxnRMAxbLH93eMifGhlCXfbF9qfcN55eF8aEPFpfzC8y+PmcE4HJp10wXwI+iQf9UTBtwU3jp7xyuHYypw78KwlHOPNSraYdWmFbsl/CQ1wrXemFR9n/FtiPLeHnT2A6dgZPC/3eFl/52WFDqmyjWiuRk0OlIO2yEaUDbhxfdhRTvy3tHaIfvgL0zXdY0t0O7ICDvDUzFnAZimtHrgoC4EDgq02VNYxu0wTSl1TB7uAXYdnxgmuIVWMzJnpkqYzm1Q0eCFhDoH/rDbuH80PBMTAHOzXRZ09wOTbB9y2sGYXgRcGK4diS2kvAaZnucU5PEdB7uZVwBUAui+yQmEO6BbdH2ebg8CXvAN8W0AFlrGL4mgt3YP2oeYm2x5dCxGvZfDecfYEvK4zJSyPLjW+DZ4CG5D/YinKCqBdg2TL9iA12uBtiN9Ye/VXUq0AnTIL+mqktDluXAKOwF+GdmSlk+qOovqrociypQFfhVi71HW2NLZHdg7ZGThH5wnaouw5b/Xqc4xub3mMa4C+Z9XysTZSwPVPVHLNYqQFesDb4M559iz8JiYI9c9qBV1VnALWp7cvfAbGpfDpe/xNrhIMzprEomypiNuEBYQVCb2lTDPEpbhHMw4eh7TEtw6Gpuz3oi9QUb3GthwjEi0hQLPPueqh6pqlMyUMRyQ1UL1JwAamPhU6pF2mc7zMN6F2xmnJPE9Yd/MDuoLQBEpAkWammkqm6tG04csYWYZmTzED5jY2yi+Lyq3qJZHoN0bajq0uBQsTXwj6ouCJc2wiaLBwLnRSaRuYqKSGPsnbA0TBTBtIdfY5r0E1UtTmumCpluVHV5iMzRAVs5iP3uTTCFylbAyWEy5ZSBXAxKWSGJek6ugQJs3+KDReQ2bJn4YsxGKLY7R9Z6iMWXfQ1tMhVbKr1TRE7Blsw3wTRlGxILMW3ggSIyDXMkOAuYpKqfQ+k2zSbiQ6GsoT98j2lCHhGRz7Blom2wUExZT1l/Q1V9RUQewmyMJwF1sBAzL4TPkWztC1C2MVJV/xGRkcBFIvIqFmvyTGyy8EXQGGU1a2uH8MzMFpHnsNiTf2CBmE/DHK4+DO+KbA81VJZ35kosBFMXETkZW0E6GzMx+FNVF6a1kDmGh50pByKxoypjXm8z15C3HWYX9B+ss/8EdFMLF1CWB6RCEhc36jjgS11DqJgQZmcfzB7kL2CIqv5QLoVNM5HQMnmUnNlG88R2YumIxZ9rgtlJ/Yx5WeZSf+gGTFTVv9aQvzNwBGZMPgNbPvyuXAqbRuJikJ6MhVn6SVcTJkZsN5rTMa/yucBN4f6snRhAqf4Qez4SCjQisjU2Rh5JcfDxAzQHQolE6l4ZCyr+T/wzHhkbagODMUFwFjAZc7LI6rEBSvWHmmsS9EVkZ8yecA9sNUGx8FxZ3x/KGxcI00xcXLnHMMeIu1V1QoK8sQe9HmYz1wD4OtgLZW0ssbi4Uf/FNqL/BJvZL1/TwBWWPPKyte7xRASAWthgPgVb+iy17BnpO80wDUg+8LZmeWy5uP4wDFsCHAdcqqor4vKWeLGJ7ViSF7ElzFri4sq9DwjwIBZvs8yanWzuC1Dq5X8jNvnps6bfWETqYl7WNbGgzFlvRxp53mthe9L/jk2E12giIyLbYDa1P2f7uwJK9YersfA5D62pHURkEyx2byNMQ/p7eZQ11/Al4zQSW8IJD/jnmKbrAVazxVrsxRdsYxZgy4WxF0fWPuCRl/9DWFDZ4zEtyLIy3KuYpjTrCf1hZaQ/zMe8hRO+zCJakplYaI3Y52S1gBzpD09ijgBnAN/GC4Mhb1QYlIi9VNYTXt5VsSWvBcAFwOR1XebL5r4QXeYOS8FbYu3RFAu3lSi/hqXA98q1sGkk7l3xBfa8f4+FlFnTPRrVhGX7uyKuP7yACf0Ps5Z3gKr+SnhfOuuPC4RpIDbDCdo+wTydFmFLPb/F0suq0s9mO5AYYlsJ7YFtKzU+m5cz1pWYJij87nnAI5gAcDzwRxAMEvaHRL99Ni8NxgjLPHsCl2Jaz7X2h2zvM2J76cYLtNtjGrELgR/X1BdykVg9ReQybOXgeODTRJODaP5cI/KuGIBp+86MrRqsYWwo03iRDUTHyHB+FRZY/ATg81yaCFZk3Ms4hYjIZlDyhR06+FaYEfjUiBYwJz3ARKS2iNwSbJ2i1MQ8ImfGa33i7s+ZPikizaHUIF0F8wZ8B/g9di2H+0MdETkjwaX6WLzJXzeE/iAiewE3haXOKC2wSAI/RYXBWDuISOVgT5YTJBgXYuyMacY+igqD0d8/l54PEWkt5ikcJR97V3wVNSHJ1bEBikxAEgmyO2HhY76MCoO52g4VhZwYbCsCItIK+EjMYSKWli9mLN4aWBgG/KI2Dw96HRHZLgNFThcHALticePiWYXZeMRedNE9a/uKSMdsneHGIyKdgNdEZO9ImmCeoZsBC+Jn+KE/1IxNLHKE84BesYE/gmB9pDKYVh1KaIzOFJEWudIfMIGntpb2elyA2RV3DlqSotUDEamO2dluXs5lTQthefwnETkikiYiUg0TjBfHbOBi1yPmBVVCm2T1OyvUtw22vLlD3OV8bNvKUkJzZGzIib4ARf3hGxE5IZJWObwzt8SUB8skEkopMj40L/cCbwBk9cNVwViKxTx6NjLrWaGqS7AwEceIyB5RoTAICAcB5yeYLWYrozCPv2UicnJ4qaGqX2OBQweJSAdVXRl5uBthXtW9JHfiqC0BHlbV92MvuFDfRZgnaXcR2TR6Q6j7McBhYVDMBR7GPB8LRGT/SPrHmB3lgyLSVEuGI2qJeZAekSsaAVW9U1XPDS/1viKycUh/B4sreTnQNiIUVsI06r2xmJO5QD7F+zID9kyo2RJ/CxwkIrXU7GyLJhBBYHgkCIVZPUEI9Z0KHKmqbwSlQYNw7V9MUNw7flIYnoOTgXPFnJBygSbY+DA6lhDeC0uA/2HPf3st9roGQEQOBK4SC9zvpBAXCFOEqs5R1TdDx/2fiDwfufwWFiqjv4jsGZkFbwZchGlJ/i7/UqcWKbadXCEie2K7rwwOGgCA+7DQCO+LyIEi0lIsnMggbEb4aI7Yx1VS1d9U9cEwC35ORK4EC64LvIhpUXvHBv4gALQFzgHahEEx61HVuWoBZI8G3hSRviF9IWZLWQV4UUQ2E5FqItIBCyGxFbYzSdbbjMWWScNvfCDW388NmiKwsCF1sX7RU0S2wILQPwz8S/HOFFmNqi5W1QdVdaGIPCoiV0QE/nuwSfUbIpIfWyYMk8UuWBDumhkpeAqJKQNUdXToF58DV4tFEgDoh72XHxazu46xOWZP1wTbiSTrUfMEvlcttuQDInJr5PKbmPb8ZhFpp8FRJgiBx2FxBrPWeabCohVg/7xsPwgbaGNxo+phTiSzMA1RLM/pWKDp2Zin8bNYEM1vIvdn3Z6L2LJfS6BqJG3z8PdU7IX2EMV7r+6DbTtViMVR+y0c22a6LmnoD9VDn/go/M59Inn6A8uwHRYuB24Leb4mBzdkD33kdsxs4NKQJpgGbAImDPyKeeD/CWyT6TKnqN6xfl8PuDn83zf0/zuB5iGtE7YF3cJw7RcsFE9+uJ6X6bok2Q6VI/83Ae7FXujnhrR8bL/yP0M/OA+4BltxWABsmek6pKldhoR26A80C2l7YyFn/gjviWHAj9n+rlhDf6gX6rgQuCKSfnWo90RMO3op8DwWnWGrTNchF4+MFyDbj9iDiW0zNhvoiW0/d3UQeB6N5N0X8yL7CXgbuDvygGelEIB5Do8FDgnn74UBPWYL0ysIPo9EXo6VgIOxXTcOx7bqy3hdUtwf6oTBbE+gFTAGWxbrG8l7BrYzzUJMUzAs0h+yVgBYXdkxO7E7gsBzWSR9E0wwvDX0iY0zXYcUtUOsv1fFHCbeANqFtMtDOwwkCIUhfU+gO7Bj5P6sHBtW0yaPY3aRLcJYuAq4INJOe4dnZTo2UXyDHBEGY79n+D8v8v+toS/cgG1cQHiHPIXFp/wftnNTVr8rVtMmo4FDMPOIuzGTmmsi10/GdilaBkzDnPFcGEzX75HpAmTzEXuog4BzQXh4twxpjbAZbgmhMFyrFneetQ94EHyexyLEfxYG8R0j16NC4cNENIm5dkSFOeB+4F1gu5C2CYmFwqqhr1SN/5xsPOJedAdgNpHHR9qmIaYZKyRoCnP5wJbEN8cE/+0IGr9wLSYU3slqhOCoEJGNR1x/OD2MhzuHZ6Q2xVrjC+PuaxvGllqZrkOK2iHW//OxGIst4q4PiAiFzSLpVSmpTcvasSFBf+iFKVH2j/zm98YLheHaxtj+9jnRHyrqkTPhDDKBmrFrdeBKzE3+TcxAHFWdIxaIGaCviBSq6lnhPH4/36y1hVDVRSJyLKbl2hYTgr+KXF8htjsJwFBghYj00xzYaSIeLd6C7EBshv8k8H/h2q8iciE24J0a+sM9avZ1BRpGvWzuD6HssaCyIzBniI3C5atF5HJM4zMA217qdhFZqaqD4z4j6+0GochJaCwWZmgqFoN0hYSdJFT1dhEpxIQiFZH7VXVa9DM0+50oYv3hNCzU0CBV/Sxc/kdEbgn/DxYRBe4Lv//UHOoHsYD0tYHXsKgTtUXkduARVZ2vqlcGc8prgUIReVhVZ6rq8vjPyUglUkSkP/TE7KgfxEwjUNXfROTukPWyMEbG7AqnZ3vds4JMS6TZeFB6ljMXmAMcH9IqU7x0uBG2fDwbeC7TZU9Te2yM7bgxFhMMDyduJovNjE/GZsF3ZbrMKa5/dCkotvwzD9g7dp1ibfImmLboG+DKTJc9Te0xCFvy2x9zFto19I+5wFEhT/NIW52X6TKnsO5Fqwbh79GYiUhMMybRfOH/frnWDnFtchAWbLkQ6BfSqkSu18GE4mXkmNY40h/yMMHnXeAKbHl0FTZJbhXJfytmU3gXUD/T5U9Tm5yIOVHOBE4NaVHN+SbYxHkucGOmy7shHRkvQLYemI1czB7o0iAITaTYoaJSZPBvFAa8N8jyJaBQn1I2YkA1zIHihdAWR8YLhSHfSUDHTNchhW0R+42rRdLuDS+/p4DG8e2GLY18gtkMZrVxeIL2qI1NDAbG1y2k/wrUC+etMGP6LTJd7hS3QTXMNmqHcH4wNiEcE6t7tD+E/4uW1XPtwCbFl2GThHGR9OhSaG0sCsGcXBOEQn/YDXgU2DqSfgKwArOvjgqF92OOaDk1NkTql4ft474cCz/VMJYeydMWszf9HTMzycm2qGhHxguQrQe2PPxe5PwSzDPwGWCTkBYVCutF/s9aoTDuoT0csws6lGATQwgjEoTCwzDNYGXMm27HTJS5PNokDOCPRdIewmb6F8UJATHNUfPI/zkz2GGTgsnAg5G0mKfsJpiH4HWJ+lOuHNhS+TzgpHBeFXMUmQ+8BNRdXf2zXShc3e+JTYovxbSA/01U3yAUNk5n+TLRHpjpyN/huagZbSfMCXEFZl/dMnKfRP9m67GG/pCHrST8gU0E6sXnx1aemmS6DhvS4TaE68/TwG0iso+qvqeqg0JcqV5Y7KRr1OzGKokIqroAiuxAstIuKM5GbCTmYbwMCykyVkSeU9VHQyDZp7A2egILM3E08FhmSp528jFnkc1EpJWq/qGqZ4fguncCiMhToQ9oiFM4I6RXyvL+oJHzSpgQ/AWwnYhspqqTtXgrslmYpqxB7B7NgbiT8ajqpyLyGnCdiLygFqT9DcxkYhjwuIicpqoL4+uvWWwnFYtDGv4/A/MkLgBeVdXvRGRoyHqpiAxX1ZPUbOvy1YL4/4M5p+UManbm72OToV2wvas/xOwERVVHBNvBJ4AGInKuqs5WXbf97isicf3hBEzr9zu2XeV4sf2r8zHns0IRuU5VF0RsbKet/tOdtJBpiTRbD2ALzKP2yrj0K7BwI8MJy8e5dgA3hbrvgRmKN8I8rOcA+4Q81TCD4Z+wATCnQgVQejm0CzbTPzou/QlMaL4AaJDpcqew/tGZfJ24a11DW9wPtI6kN8P2J706URtm40FpDV/V8Hc3YApwWuRaJUxT+DdmT1Yz0+VPU5u8iIUI+Q6bKC0DLsRsBetgmsLZwJOZLmsa6p5w9Qc4CnMw+5awUoLF4YxpAk8NY2jWrh7F1Vci/4/EzAV+DO+Dv4GLwrU8LCj5T9gyck6ZC2TbkfECVPSDNSzrAddhS0Nt49Ivw5aH+me6/GlqkzFYcO0a4bwRFj3/kVhaJG9zcihUQAIBIOpQ8jgWeqdx3ID4GGZTeHR5lLEc2iBa59sxQ/lXgGOx/XpjL7hVmFfl2ZhN6bPYpKFdpuuQ4vaoDhwRl1YTCzr+enzbYdryN3Pl5R9Xv0uxZcDdCROgMC4UAgeG8/qYic0qYGimy5zCusdCy1TBNIF7UdJm8BjgUyzmaCmhMNpHMl2XFLbJ5ZjyoAtmTlIbWx4vBLrF6huEwb+xAP1ZP1HM1iPjBaioR3hQo4FlP8LswXaK5OmAzYJjOy9EPedOjBcesvGIE2wE22LrT+CGkNYeE4qfjwiIx5AjO00kaosgADyDeU82jFw/GVsW3TnWbyLXriLL7cNCPaLC4GOhvk9hM/zfsMDTdcL1AzEBeT6mIfgm2/tFZEwomiiGF9yCUL+jgU3Dtb2wmGrHJupH8e2ZTQe21FcvQfqTmKlITDjaLPSRZ4HqkXwbYVrDnFhFifSH2pjQ90sQev4CHorkOxZzKPsM6BzfH7L1iB/bImPl86E/xDTnbbFJ4X8xAbHI6x6bXObUZDHbjowXoKId2JJG1GO0GrZ7wAPhIZ+GLYXFAg7fD0yM5M+P+7ysFQqJWxaMDHpPYssb+4aH+7mIELAztnfzAZkuf4raoDolBbt8zBD8bWwp7DVCKI1w/S3g7ch5/ECZ9UJhqEdDbFuxbpG0hzDB8G6KjcQbYcbhHcjy5SDM9vFIgpc8Fij3csxbehdsqfTncJyH7brxGhZCpCj0ULYfmAbsayxESqNIeqVQ31fCeTtssvgcxc4UVwGbxfJnui4pbpeqmPbvfWyy2AUzr/kHGB3JdxwWhmkq0CHT5U5BvfOx2LNXRtIk9JNxhC1cMTOr+P5wCbB7puvgR/jdMl2AinSEDnw68EIk7RuKY2dtjWmBfsNCzLwInIK5z5+b6fKnuC2iwuCd4aW2Wzg/BROMl2IG47F89bHQCl8Q2Y4rWw9s2e8s4KpwXgkLnXJQOD8ImygsDC/Im8IL72tg35An62f/CdrlIcwe6iNK2ghWwSZIMaGwTqbKmKZ6twp1G40Fop8afutakTx7YbtN/I1NGlaF/3NCExap5+OYw8iVlBQK78NWTfbD4shFJ4ubYvbEZ2f7c0GCHZewfdp/B7pE0moAR2BC4dBI+inhOcr6SUKo433hPRi/48z9oU22Cc/BCxSblWwZnqUTybHJQbYeGS9ARTqwWc3x2Ax/fBD83iVuW6kg+ByP7TH5TxgYX87FTh0e4F8xTUiLSHp/bLY3EtOg9sSWARaQQw4kmE3L75iNy69Y3Kzo3rN5mLPEfZiwWBiOmzNd9jS1h4QB/I/wwt8+1g7hb354CXyPLaXmjP1oqN822ERoASYQN43WP5JvC2wy8X7oDw8SMSnJ1oOSJgODQ92uJISLwbZl+yOkv0RYbcGWiB/DJtKty7vcKW6D3TDteMO49P2AJRRPnGPLoTUxbep0EsTcjO872XhgS+W3YBOgiyJ13yy8TwuJbMyArRw8hk0sW2WizH4k+B0zXYCKeGDLPYXYcugWIa3EDgSRvIdhQXgLgP9kuuwpbofzsWjye5BgY3WgL7ZEugzzIHub3NmIPmrn9TjmNfsbYYmH0rZkeZh95XmYM8EcciDuIomdqapjMShjE4JacW2Rj9kVfkEOxRGLPAMrw/jwVvQFTwIHNIq9KCeTI8IxqxcKm4S0fTDzmi8x56KLsNWUeUScLLL1wJaCByZI74w5150fa6eIYNQttNOumS5/GtulNib4rgL6hLTKmDZ0YhAMDwj9YQRmW5z1/SGXjowXoKIclFwivRKzBZmMzfDrJ8gTHRTrYKrvYdiyWU5oCjGtxnuU1n5EtxnKwxxL6hCWAnLhCIN57AX/HKYlnIoZPteN9oF4oQnzMJwEnJ3oerYccf29GiWdAqpjoTT+wRxsEgmFTTNdh1S2Q6Rux4ZjKeZd3Wk198UEyDoh7xmZrkuS7RAVdKMOdHcHYecqij2L24ex4xtsf/cRq2unbD2wpdLBQPu4tlgM7BWX91AsDFHOCEBQOng2NimOCYV9Y30FW0Uag00UfsBWnnKqP+TCkfECVKQjvOS2DP83AM4MnfeDyECXF/0buXcwpv7OevV/pE6jgM8i5/GCTydyxEkirl6xga5q6AexWIsPhUH9DoonCQl3FMAmCO+SpZODOGHwBuB1zP7rkcgzUHVNQmEuHWFsGI2FU4kJhl0SCYWhv1SJtQUWemkycHqm65GK/hDO453n7okIhbHlYwltUZscWC5P0CadKV4aj3mWN8f2Kv8X86LeC1tK/hxzsMiJZyNBf4gqSOpHhMKL4/K1xgTpGukuox/r8btmugAV5QgD92PYEmnMIaAyZgcUEwpjL7ymmBaxRRj0Kod7vyaLDOnDg9mRuAC5ESHnRixkxJ4JBJ5u2FJ5s0zXI01tUxnTcLwc97J/mGKhMGYc3RTTANQJ/aE6Zl/6PBGP9Ww5KDnjfw7TjN6AxZj7BxMMY3WPCYXzwoswV4MtNw5jw3eYDVlsO76umN3YKMzDfmNsYnhjuF4pjCGFBO/abDsoOTm4FAu6/y62FBiNyBATCouWj3P1iIyRXTCnstFAm5DWDFMQLAp9Y1p4f8T6TNYKhcR5ywN9MPOQ4UQcSrCtWouWj1nNxNmPinVkvAAV6QB6hIFuAsVBM/Mp1hR+gzmTfIbN+GKagi0we5ltM12HdahrlfCiLwwPdI8EeTYCZmDL5jvEpT+BxdtqWF5lzkAbXRHq/ziR+HmYUPgzpjHcP/SF9yKD3rahXbfLdB2SrP91mG1oLK7i5Zi96KIg9MS8R6sCJ2DOBC0yVd4U172EHWD42yLU+0dKCoVdwov/D8zO9FtKmlV0ImiQsu2Ia4fnMS/rJzHTiUJsohDdg/cezNv0plwaG1jDSghmM7kImxC1iaRvh2kI94y8K7JyRYXEXtXPY0vAwzBHskLgzsj1ephQuAy4ItN18KMMv3OmC5Cxiq9+0+2DMdX+/1FSKDwRM5Kfju0yUGK2Rxbaz2GCbSEmAP+LGX4fE5dnB2ybqSnYcuGtmPPIfHLLmzg66406zpwf6v84JXcduDe8/GdhAnOsP8SEwnqZrlOS7VEHCzd0STjvhzlOHYc5lCzDPPGjQmHWPQOrqXuplzelhcJJlBQKt8Y87y+j2HYwv7zKXA5tMgATBncJ533D2LEKGELJCASPYR7oOSEQRn77WkHwGYFNoneN9P99MaFwNKsR/lf3zqnoR3j/fQ7cFUm7PfSHXcP5tRRHWHgkkq8eNkmYS5bHId0QjowXIKOVtyXTPeMHbkoKhV1DWiUsGO+W0RcGWagKjwxw3YNA8yi2/DM9DGrjMc1Xi5CvFWYj9i22ZJaTBsGhPzQM/0eN5i8I7fQkES9qTFjeM9KeWdkfQnnjvefzgsCzUajnn8B5sbphNoWFWCienBAE4+pfFdP4XB9tk/C3RXgOJmI2hbFdGCrF583WI64uLTE7yd7h/BLM6/5I4FxMKLyDkjEpc2LJOPI8V8cmxT9iqwHfY5EEbiKETcE0hQsw84GciTuJ2U9/FN6JjUN/+C/FTnOXhv7QK/xfSEnhsS6wUabr4UcZfutMFyCjlbfYcYWYABgvFB6BbTv0HbB3gnuz1g4kUoc22DL460EAaI7Z/3yBaYA+JbKUTPHWdbloIJ6HhVBZRHH4jKhQGNOIPEGC7deyXQCI1OMpoHv4PzbxOQPzsm4XyXcfZl85nhzcbio8G2MwLcjFkfSY9q8LplUfi9kRZv14EKljdJk4FkdwH0wQ3hvT9pyBTYAaYFuxFWBRCZrFf0a2HpSc+B+I7cLSkmKt8D2Y2c2NFGsKu4Rx4vZMlz/FbXEgJvifGM4PwRQkXbGA0zHhsAU2SSwEhme63H6s4++c6QJktPLm/TYGWxI8JIFQeHcY/GYS9p3MtQPbeWUVcHhc+u9Y+IRCbIm8f6bLmsI6N8e0XztS7BGZh8XImoDZByYSCkeFvvAi0DbT9UhDuzTBDODPCOexF995mIYwtjzUIAiD55KFTjOrqXspgQ6zDR6O2UldEnetDWZXXCLgbi4d2DLx45TUFl6N2Us3iKSNwCaRc8iRUEORulXDVgY+wjTGEtce94Z6R/e4344stRVcTRvEHOVexrTiTSnWnN4Q+kN0t5oXMC3qH+TAjlUb0lGJHEZEJPJ/Xvx1Vf0H24z+G2zgO0BEqkburYppQIaHPLnI69gM/1oRaQAgIs9g0fVPCEcLoLeItMxYKVOEiOyK2UA+jw1aQ0VkW1VdBbyDBU0tAMaLSBNVLQj3VcEE518wm5ppmSh/KhGRSpH/RVVnYVrxvQBUdUW4/BqgwF0i8hhmS7oH8IaqLivfUqceEamiqoUiki8irUSkYUj7Adup5hPgPBG5LHJbY8wGtwXmaJZTiEhloC1mLlAlcqky5kmdH/K1xMaKi7Gl07/KuajpphlmJrQ5NkHS0FeqAKjqhZhAeDYUPUffqOrK0IZZT6jzUkwg3hSbGGoYP2pipjbVAUSkBWZ//BC25/eMDBXbWQ9yWiCkeNCqqqqrRKSGiJwjIteIyH9EZOPQ0Q/FNENPAEeKSGtMA7AJ8JiqXhbuLyVUZjuqOgd70DsA24nIc5im7Fhsn+JnMXvCLVX1z8yVNHlEZE9MCPwQcxK6GmgHnB7rI9gE4ALMU/IjEdlERGpgWsUawGmqemh4KWTt8xNeXIXh/2qqNrXHQidtFMsjInmqOhVbMlqCBd3eCLOt/a38S54aRGQzEWknIpVUtUBE6mKa8PGYAPiqiLRT1e+xJcH3gYtFZLiI9MX2sG4G/JVrY0Nok5WYI9EmmKNMjE+wScMwEbkSCz21HfBHGEuzmqgSASD08RMwLdg+ItIvpBeE56MaJhDGBESN3Luy3AqeIkSkuojsHJeWD6Cqj2GToKvDeSFmV9gauDn0h1sxzfpHqrq4HIvupIJMqyjTdWCd8nWKA4bWwzwDZ2OhRGJLoceG63mYSnw5pv2ZijmVxGyGst4mJkEbRW2FJmAasJmYHUxefJ5sPjCPwAJsGaxuJP1R4OvwfyyYcB7mLPA5FmPsQ8yg/JscbJeHsBfaKCzMzlOYYXyiPVdrYS++rN6CDdNqzAi/6yahTmMxYacP5ijwdWiX/4R7NsPC7vyJmRREIw1kdV8ggf0rNpnOAwZhQkDHyLVzMAH5L2ypeJtM1yFF7RAb6ythS8VVI79xR2yrwh+AS2PthmlRfwEGZ7r8Kai/YEG2CzH7yK7Ra+Hv8ZjN4BmRa6eFZ+LP8M7Mif6wIR6xHzlnCJoPFZEzsaj5kzF1fhesM/fBDMW7ANdgguIAVX0u3N8bs5EoBG5VU/3nqWmPco5Ie/XB4s4NBm7WHOoYItIYiy/ZDLNxWikiNVX136DtORGbBOQBL6jq0+G+Zlh/aYkJB5fkWn8QkbMxQa8j9kxUwrTjfwBfYROocZizzaeqOi8T5Uw1ItIFC7r9BTZJuB64VlU/C9dbYvZhe2H2Yb8GjXAtzIZyWnhuKmsWaoISISK9gN9V9b1I2j6Yh/FFqvp4JL0htky4JBf6ROyZFpFamPZ3E6x+nwF3qOpUEemECcj7YQLxSmzMaIrFoM36fiAibTGHymuwFYHPsGfjD1VdKCKNsPHgd1X9T+S+Vlh7FKjq3HIvuJMSclEgrK5h6SK87M7BXub/AhNV9ZpI3h0wAehf4By1pbH4z8uZl/+aCA/058AHqnpspsuTSsJy4LlYX/hMVY8O6TtiWqJ3MW1HbWzXjatU9bbVfFbWCgBr68vBfnYrzGb2O+y52AHbiqom0EHNzjAnEJE9KHYqWwbsoaqLIpOkVpg5xQzgMC22qYzdX0nDsnu2IyJ7YyskBZjT1OOY5rxQRO4DDsPaJ+ttZ1eHiNTEJkH/YDuLtMIcz2oC+6nqtyLSAXtn7IQJS4NU9d1wf9aODfGISBtsonwSZiIyFrhbVT8UkQMxTeLZqjo8c6V0Uk6mVZSpPIBtMDvAqMfXuVj4lH+APiGtCsXCcDdMG7h/psuf6QMLJVGIDX4ZL0+K61YH0/bNxCLrbxX6xMOEINKYtngQ9lLcPtNlTnH9o4G3j8PsJE8lbjs1bKlsBiHmXEhrRI7GEcO8zX8P/f5ISu/RehtmalIv02VNcb0TeVVvhcUj/ROLszcK05QdBHyMaQkT3pvNB8XhZa7HYq1uErnWBfMwnkXx1nTtseXjT+Oek1xrl7zwrrw51HUFFnT8EmzyMJKIt7kf2X9krVH8atgOMwA+T0S2B1DVBzC7qLnAZcFQvIDgcIJpiOZg241t6LyJzXqnZrgcKUdVF2Faj9uxXQX+D3haVc/CNGGo6gLMo7YQWybOCYLGa1X4/0VMw3EutjT2hIhcEsur5jX8GaYZjKXNUdW/y7fUqSeRE5Cqfgwcg9lFXQh0iHMQmUuwmyyXQpYDQVMccyhqJiJbBO3WRFV9iuKJdUtMEDwIs5XrCUXOBFlP7HeO1KcVxXbUhGvjsF04lmJORfmq+hMWjWAhcGYwO8mZdolQqKoFaqtqR2GTyD2w5+RQbAK1SQbL56SYnBIIVfVJzF7wQOCSiFD4IDbL+Qd4KiIUgnnX5gFZbweTLGpexF1UdUqmy5IOIkLhQEwr1DSkr4h50mEv/2mYIJDVxAQgVZvui8h1mJfwEcDOmNf0YuA6ETk4cutEoEtYQssJgsBTKCJVggC0j4jUDd7ln2L7mG+NCck9xELP7IHZHU/CdvHJesIyd2xy8BBmDzYB+FlEYhPmuVgM1l2wAOQbYzEqdxaRppkod6qJtYOINBaR+0PyYkxDvizkiXnXvoc5lO2MOV6gqpMwAWkFFpmiXvnWIP2oqsa8rlX1T1UdisXrvYLiMGwLMlQ8Jw3kjA2hWNywWMy40zBD8XcwG4+vQ3pvbMeJOtiA1wDzJq2FBZ7OCfsPZ80Em8JTMa/R8VpsU9gOs59bCBycrTP+oMVYEf4vsnMTkZHYpOh8NYeajTAPyReAC1R1Sch3NrYF1Q6qujAjlUghEYeB2ljkgU0IIWMwk4Fhak4je2D2cxthNpR/Y4byh4RJQy7ZDD6FLYcOxDRiR2C7kbyLBaGfEsnbEtMa/hIEoawm9juKxQl8Cwubsif2u38OPKSqfeLuGQJ0Ag5S1eURO9PNgOWq+nv51iKzhMlmE1WdudbMTvaQ6TXrVB6YU0AszMzZ2ID/XyL2YJid3K/AfMww9lSKww3kxPZjfpSpr0RtCp/HBIT3MUEgFmoi62yCMA3nGMwzMpaWj2k+fgbuDWmbY1rx54AaIe14LAhvfXJsO7pQ/y+xSeIRQGcs5NCfwNMEuzFMK/Yr5miyHZHtyzJdhxS2RXtMC94rLv3KkD4A87DNufEwMtZXwfZyfy30hbxwXIqZkNwTzmuGZ2Ui8EjcZ2V1uKEk2jDn+oUfduREJPUID2Mz2S1U9SERWYUtFSMig1T1a1V9NNiOXI65zj8RrueMh5izdtS8SR/HduC4BFsS/BELH7Eii/vDRlgszcNF5B9VvUlNW7hCRN4AdhOR7piN2NvAmaq6JHjc9wKGqup32IQpl9gBCy5+uQavUOAMEbkCW/o7UURuU9VPReRULGTV/2kIQJ6lfWF1NMDs5b6BosD9y1V1QNAGnoFNHHJC+yMWPBpVXRb5HcdgNpLLVPXLSN6nsTHhRsz0qAATDAuwKAVFobpUNTeW19YR3QCibmyo5JQNIRYiIk9E9gNQ1UexqOr7UdKm8CEsCn/UmD6XBnynDKjZFD6BbcX2PBZQNZuFQVR1Ota3PwBODXaDMT7EPIZHAp+r6rFBMK4P9A7XPi/vMpcTzTGb0ZgQFBMSbsOcqc6meLeJ91X1gIgwmLXLxJJ4B5WJ2HL4cQBqS6BVw7XbMA3xXuVTwvQSxvzngK9F5CsRiW0z+Dg2edpWRA6K5Q9C8D1YuJnXsCX0YQSTojA2bJCCoJP7ZK1AKIHY/yH5PWxpqHssn9p2O1di9jF9RWSXkP6S5tiWU866E4TCwcBxkQE/64TBmANJ0F78BtyBCTq9xLaUQlVHAkOwINOIyGFiwYgfwPb07qVZvj0hrFYI+gTTep4Hpi2KCYXYhKAZtrpQgmwXBrXYgeSUYCML5kX/KrCfiBwDJhSGa20wh6qsd6IRkb2w+HkFmPPMAsyp8GJVHQGcRegTYjFJY6xS1e9V9WJV7aOqt2lxQPqsGxscp6xkpUAYZu1KsceXhr9/YdtOHR8T/EL649gS0AnYPr1Errn6ewNHVRermkddtg34Mc1O0GblhXrkq+rP2E4b9bCJUP+Q7y7sWViGOdD0w+wp91DVbzNQhZQSBPrYvuXXBqN/MGeaN4BDIkLQsnCtNRZ7MetD68SI8yYeji2B9hUL3L8EC78kIe3SkG8bLBDxYsyOMmsRkV0xe9H7MbOIc1V1X2AEcK2ItFHVUcD5mDnBpSLSGYqfpYiigZDu7wonp8laL+PwIhyLLY29oqqfhPStMK/Jx1T1zjjv44OBN/3BdnKB4CX5JoCq7hfSqkW0X59h2qBJWMiM4ap6fchXFVsinoeNA0syUIWUEudN/D5QGTMFGBDS2wPPYqFCXgeGYo4jsThz+2ezRhBK754ShMG9sIDTP6vqjIiHbAfgBmB/LNLCTMxe7mBVnVD+pU8NkniryvxgDrIF5jh2Zlg9QkSOA+7CTCruiNoUOs6GRFZqCAOtsGWAE4CHRf6/vXuPtnO+8zj+/kjCMCWk1dT9Wgy6tA11DTLRDkqqaRrXDhFClOjI0JlFF9qa0lWTkVawKNK6yzKtW5cuipm6M0NdWiND1C1WiURGCOE7f3x/O3mynZyTk5xz9tl7f15rZe1znufZz3n25uzzfb6/3+/71b9L2j4iniSrqZ8maa2IeF9LCpDe5mFiayG1m6JtJN0Ai4dCay24FgB7kVMmfgd8S9KZ5biFwOyIeLfZg8HK7/eHklYng8G3yILTP6ndAEYWFD6YDJBPJgOg2sKi/WpzBhvwElZaXaZ4YNk2nGyxdlRE3FOCwXWA/ZS9nGeSVRZ2I8txTQJ2beZgsFgIXEP+/38tZK3Rsm9R2T+/lgGMiGuBU8j34bwSKJu1nabJEGoZfViVzbh3J3+hh5B11e4ghwKmknUIm+NFmnWTpLWAo8lV83eRQ36Pkxmvb0bES+W4zckFVnsAN0XEPzfkgnuQpM0j4vnydS0DNJZsQTYqIv637Psi2YloLplRXUQunNiFLE31SC2QarYpA9Bxprhs35lsPzeWzBbvRxaaHkBmzy4DJkfE/L6+5t5W93txf0R8o2z/DbB+ROxQvq/W6RxP9mw+qNkzxWYroikCQpWyCOXu/wCy3uCtwNuVeUBIOoIsMDqezH7eXeaNLC4V0PdXb9a76v74DSHr7e0PzK+ulC03T+eSddW+HBFvNOyiV1IZFnwIeDYi9q1sP4gc/htDzgs8hpwz+QZZZuRi4JTq50Z5Xoc3nM2gZIS/Q5ZFuT8ixpbtnyffiyCzyduTK25/St4YTAN2iojH+v6qe1/l9+KfyEz6X5PvwVcj4o+VKQYf+9tQP/Ru1g76bUBY7uonRMTx5fvBZJPxdcu/F8mJ0VfX3+EqOw58lVxRODkiLu3Lazfra+WP33gyM/5sROxftg8sc6hq88Y2Ad6PJq8xV17vCeXfQ7Gk28xuwDlkYLwKmQk7gxxCPJJcUb5NtFh7xrqbgvsiYkzZPoqcQ7iILDV0U9m+AznX+rBWnjNXeV8mkSuot4+IZ+pvAJwwMKN/FqYu83iGARPKL+6xwGlkKYTx5J3/ReSd31qSLoqI+bVf8oj4vaSnyT7FwyVdBktWI5u1mli60PZ3Jd0YEd+MJeUyavPoXmzslfaM8nqnkaulvytpRkSMiYj7JZ1LZoIWAvdEFtpG0kvkgoKmzAR2pvLfH/L9uCkiRkfEzZJurVtosh5ZjHsheWPdssr7cgW5ovofgbOAsSUzuHiKgP82mPXvDOFg4BCySOjPydpYMyPiF5VjbiBXT14ILBUUlv1nkENGO0QL9GQ160pdpqjap7lph0Q7U/d6H4iI0R0csyqwBVlv8C1yfmH//OBbScuaO1fZPxYYDYwERkYLlBpaHnXvy39WhtU9NGxW9MsVdeXObR7Zh/gUssfqGeSqsWqXgdpk6YnAcZIGV4LB1cgM4ZuUeoVmrS6y0Pbl5HSKXSTdUba3XDAIH3u9u0qaUd0v6ZPk1JEryHl0Xy9D5/3ys29l1b0fu0m6sbavrCw+CRgM7NUuwSAs9b6cS74vd5XtDgbNin7zoShpSFkJSRnmWhP4XERMI4eG3wKOKPvfK3f9taDwQbJn8QGVU64DbEfWm5rbZy/ErMEqf/ymAZso+9O2rLogaPdqEESWEtkXeIIsqVLrRtOygUDd+7FH7f2IiHvIFn2HR8QzjbvCxoglrSovBt5p1ZsCsxXVL4aMSz2oA8l2WxPINlMvkJOgR5fh48PIFXPXRMT48rxq0ekfAGdVV43VVic34jWZNVq5qRoUEXMafS19oZNhwU2BF8tnQksOnXek8n5MJhca7dPFU9pCWZW9oJYpbuWbA7Pu6BcBIYCkjcl0/teAeWQD9qMj4pWyfy2yCPUUMig8umxfKuirm0PolWNmbaR8TowDTgWei4gRlX1t93lQ3o9vkyus94kW6FXdU9rx/wezzvSbgBBAWSH+GbJEwviI+GXd/lpQ+K/AVWX1sZnZYuVzYhLZpaPtiwy3W6bYzFZMvyg7U0nbfwL4JVlD7BJJb0bE7bU7uVJC4GqytMY0Sc9HxI8aee1m1r+Uz4kpeFgQgGjBTiRm1vMamiHsaD6Psg3TUHJo+ACy/dZtlf21+YGjgZujCVtNmVnf8LCgmdnyaVhAWOmgsAZZK3AI2Vf01xHxWmlNdSHZceSgiPitpM+Q/Td/EBEPVc/TkBdhZmZm1gIaEhBWsnxrAo+Q5W/eBbYs318ZEVdKWh84n+xLehWwA9mPcjsHgWZmZmY9o08DwupcHkkDgBvIvsTjyBZKATwPzAcOjoinJa1LthzaE3iOXHm8qKPhZjMzMzPrvj5ZVCJpSETMiYiPKkHhamSz8auBWaV+4L7ARsDkEgwqIv5C9uZcJyLeKufzMLGZmZlZD+n1Su2SdgRukzQeslVQqRA/GNgWmF+CwUOB24HTI2JKKR1xnKQNyvNqwaAcDJqZmZn1nF7NEEpaG5gO/A3wf2Xq4OUlKHwLuBcYXnoTXwCcTrZbAhhB1hx8DHildk6vGDQzMzPrWb2aISw9hG8p3w4G/l7SkWXfe2RG8DAyGDyvVlNQ0mfJeYOvkwGhmZmZmfWSXssQVuYKTiWHhl8HtgdOkkRETI+IqSWLeCawoaTDgfWAscAgYETdvEMzMzMz62G9liGsBHBzgL+Ur8cB7wGTJI0rx30f+Acyg3ghcCDwFLBTWU080MGgmZmZWe/p0bIzZQHJVsBtETGvsn1L4A/AEcDDwPXAqsC0iLiiHLMasDbwRq2cjFcTm5mZmfW+HgsIJW1B1gmE7Ef8MvC9St3Bc8jC0geSQ8cXk0HhzyJiejmmWqfQLafMzMzM+kBPDhkPBX5Xvh4E7AzMlHSipM3ITiPDgJER8SQwEVgITJR0PCw1zOzVxGZmZmZ9pCczhAJ2A84mM4Ejgf2BrwBbk6uGjycLUu8bEXMlbQfcBPw+Isb3yIWYmZmZWbf0SEBYHd6VtBtZRmYIsAewABgPnEhmJNcARgEPlxXEmwF/LsWpPUxsZmZm1sd6NENYCQp3IVcMDwX2joiZkrYBdgS2AKZExNt1z3FvYjMzM7MG6OlVxtUAb1eyBuGGwJ4R8VxpWRcREa4taGZmZtY/dGtRSZknuEy1QK98/QBwAvAScK+kLUoAWNvvYNDMzMysH+hWQFg/v0/SgLrvVeYF7ixp24h4BJgE/Bm4W9LWHhY2MzMz61+WOyCUtJWkIyWdUykT86GkgWW/SobwUOABYLNyzIPAycCHwI97/BWYmZmZ2UpZrjmEZeXwDDLTtynZA/nZiNi97K8Fg2OAG4DTgfOqw8KlxMyfnCE0MzMz61+6zBBK2ha4BrgOGANsTmb8Npf0k8px65Ct6U6lEgzW5h1GxNMlozgAMzMzM+s3Os0Qlv7C3ycLTU8EZpVM4BrAdGBoROxZOf4zwOuuJWhmZmbWPLrKEK4CfBb4Q0S8UILBARGxALge2F7ShrV5hBEx28GgmZmZWXMZ2NnOiHhX0onAXFg8V7A2B3Au2XXkw4hY1JsXaWZmZma9p8s5hBHxaskI1pedmU2uHF6rtkHSUEnr9fhVmpmZmVmv6VYdwjpvAIOATwFI2hi4F5jSA9dlZmZmZn1kZQLCgcAiYEBZTHId8AHwrZ64MDMzMzPrGysUEJbSMR8A7wDbAjcCawNfjIgPaotMzMzMzKz/W67C1B0+MUvSzAQ2AJ4BvlALBr3IxMzMzKx5rEwm76+AecDLwPCIWORg0MzMzKz5rHCGEEDSRsCrtZ7GDgbNzMzMms9KBYSLT5LFqt2j2MzMzKwJ9UhAaGZmZmbNa2XKzpiZmZlZC3BAaGZmZtbmHBCamZmZtTkHhGZmZmZtzgGhmVk/JOkoSSHpqEZfi5m1PgeEZmYVkk4vgVhI2rrR17Oyyuu4p9HXYWb9mwNCM7NCkoBjgFo9rmMbeDlmZn3GAaGZ2RJfATYFpgOzgSMlrdrQKzIz6wMOCM3MlqhlBC8FrgY+BXy9/iBJZ5Wh2L0ljZH0sKQFkuZIuk7SBh08Z5ikCyQ9UY57T9Jzks6XtE5XFyZpgKSXJL0t6RPLOOan5brG1OYgll17VYbBQ9JZdc/bWdIMSbMlvV9+ziWS1u/qusysNTggNDMDJA0FRgH/ExH3A1eWXRM6edoJwFXALOBC4CngYOBOSavVHXsscAjwLHAFcBHwGnAKcJ+kNTu7vtIe9FJgTeDQDq5/deAIMrP5a+Bx4Oyy+8Xyde3fPZXnHQ3cB+wH3A38G/AoOXT+qKSNO7suM2sNAxt9AWZm/cQ4YBAlEIyIpyQ9BoyQtGVEzOzgOfsCO0XEk7UNkq4hA7avATdUjv0R8O36vu+SxgOXkcHleV1c46XA94DjytdVBwNrA/8SER+QAeHjks4EZkXEWfUnk7QVcDEZ0O4VEa9U9o0EfgtcQAdZUjNrLc4Qmlnbqywm+Qj4RWXXlYBY9uKSqdVgsKgFal+qboyIF+uDweJy4G3g77q6zoh4DfgVMEzSsLrdx5Xrrw8UOzORDIJPrgaD5WfdBdwMHNhV9tLMmp8zhGZm8LfAFsAddYHRNcD5wFGSziiZt6pHOzjXS+VxqXmBkgaRQdshwLbAYJa+Kf/YvMNlmAaMKeeaUM79OWAX4DcRMWs5zwOwa3ncS9JOHez/NDAA2Ap4rBvnNbMm44DQzGzJPMErqxsjYo6kW4BvkEPAM+qeN7eDcy0qjwPqtl9PDr0+T87xmw0sLPu+A9TPOexQRNwt6Y/AoZImR8T8yvVfsjznqPhkeTy1i+M6XMRiZq3DAaGZtTVJ6wIHlW+vlXTtMg6dwMcDwuX9GTuSweCdwH4RsaiybxXgtG6e8mJybt/hkqaTi0leAW7t5nnmlcfBEfF2N59rZi3EAaGZtbsjgVXJIdHHl3HMKGAfSZtFxAsr8DO2LI83V4PB4kvA6t0833RykcoE4D1yMcnUZcxR/IiPZytrHgSGAcOB27p5DWbWQryoxMzaXW3ByAkRcUxH/8ih2NrCkxUxqzzuXd0o6dNkuZpuiYh55PzGLwA/BGolaTryJrDRMvb9DPgAmFJWHC9F0qqShnf3+sys+TggNLO2JWlvcsHEkxHxcCeH/pxsZzdO0oqMrDxC1vobLel+ST8uQ71PAQuAV1fgnNPK4wbA7RHx8jKOuwvYRNItks6WdIakPQEi4k/A0WR3lqfLMedLmirpV2SdxO6sWjazJuUhYzNrZ7Xs4GWdHRQRsyTdCXwZOLC7PyQiPpQ0iszm7Q9MIuf8XVa2PbMC5/xvSY8Dn6fzxSQnk8HsyPKzVyGLU/9HOc9Vkp4AJgMjyPZ975BB6gxyMYyZtThFRNdHmZlZv1JqA74KzAE2i4iPGnxJZtbEPGRsZtacJpLlYKY5GDSzleUMoZlZk5A0mAwENyCHu+cAW5dahGZmK8wBoZlZk5C0KfACWdD6MeCkiPivhl6UmbUEB4RmZmZmbc5zCM3MzMzanANCMzMzszbngNDMzMyszTkgNDMzM2tzDgjNzMzM2tz/A9Pakc40TA2YAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saving the figure ...\n"
     ]
    }
   ],
   "source": [
    "fontsize = 20\n",
    "figsize = (10,4)\n",
    "\n",
    "fig = plt.figure(figsize=figsize)\n",
    "\n",
    "ax = plt.gca()\n",
    "\n",
    "# group the data by analyte and group name, and get the unique count of unique tube label ID\n",
    "grouped = data_outlier_removed.groupby(['Analyte', 'Group Name'])['Unique Tube Label ID'].nunique().reset_index()\n",
    "\n",
    "# pivot the table to get the count for each group name as separate columns\n",
    "pivoted = grouped.pivot(index='Analyte', columns='Group Name', values='Unique Tube Label ID')\n",
    "\n",
    "# create the bar plot for each analyte\n",
    "for i, analyte in enumerate(pivoted.index):\n",
    "    ms_count = pivoted.loc[analyte, 'MS']\n",
    "    control_count = pivoted.loc[analyte, 'Control']\n",
    "    total_count = ms_count + control_count    \n",
    "    \n",
    "    # plot red bar for MS count\n",
    "    ax.bar(i, ms_count, color='red', width=0.7, align='center')\n",
    "    # plot light blue bar for control count\n",
    "    ax.bar(i, control_count, color='lightblue', width=0.7, align='center', bottom=control_count)\n",
    "    # add the total count above the bars\n",
    "    ax.text(i, total_count+50, f'{total_count}', ha='center', fontsize=fontsize-4)\n",
    "    # add the count for each group name on top of the bars\n",
    "    ax.text(i, control_count/2, f'{control_count}', ha='center', fontsize=fontsize-4, color='black')\n",
    "    ax.text(i, control_count+ms_count/2, f'{ms_count}', ha='center', fontsize=fontsize-4, color='black')\n",
    "\n",
    "# set the x-axis labels to be the analyte names\n",
    "ax.set_xticks(range(len(pivoted.index)))\n",
    "ax.set_xticklabels(pivoted.index, fontsize=12, rotation=45, ha='right')\n",
    "\n",
    "# set the y-axis label and tick labels\n",
    "# ax.set_ylabel('Unique Tube Label ID Count', fontsize=12)\n",
    "ax.tick_params(axis='y', labelsize=12)\n",
    "\n",
    "# remove the top and right spines\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "\n",
    "# set the x label and title\n",
    "plt.xlabel('Analyte',fontsize=fontsize)\n",
    "plt.ylabel('Patient samples',fontsize=fontsize)\n",
    "\n",
    "plt.xticks(rotation=45, fontsize=fontsize-5)\n",
    "plt.yticks(fontsize=fontsize-5)\n",
    "plt.gca().spines['top'].set_visible(False)\n",
    "plt.gca().spines['right'].set_visible(False)\n",
    "\n",
    "# add legend box\n",
    "blue_patch = mpatches.Patch(color='lightblue', label='Control')\n",
    "red_patch = mpatches.Patch(color='red', label='MS')\n",
    "plt.legend(handles=[blue_patch, red_patch], fontsize=fontsize-5)\n",
    "plt.show()\n",
    "\n",
    "print('Saving the figure ...')\n",
    "fig_filename = os.path.join(OUTPUT_PATH, \"figures\", \"shortChain_fattyAcid_patient_counts_{}_sample.svg\".format(sample))\n",
    "fig.savefig(fig_filename, format='svg', bbox_inches='tight')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b3a42fb",
   "metadata": {},
   "source": [
    "## Linear Mixed Effects Models (LMM) with fully pre-processed data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb1d616e",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "data_list = [data_outlier_removed]*len(analytes)\n",
    "analytes_list = list(analytes)\n",
    "maxiter_list = [10000]*len(analytes)\n",
    "arg_list = list(zip(data_list, analytes_list, maxiter_list))\n",
    "p = mp.Pool(NCORES)\n",
    "analyte_stringent_model_list = p.starmap(lmm_model_2, arg_list)\n",
    "p.close()\n",
    "p.join()\n",
    "    \n",
    "joblib.dump(analyte_stringent_model_list, os.path.join(OUTPUT_PATH, \"shortChain_fattyAcid_stringentModel.joblib\"))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4494aada",
   "metadata": {},
   "source": [
    "## LMM using data with outlier values, but excluding missing data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca596e5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "data_list = [data_nan_removed]*len(analytes)\n",
    "analytes_list = list(analytes)\n",
    "maxiter_list = [10000]*len(analytes)\n",
    "arg_list = list(zip(data_list, analytes_list, maxiter_list))\n",
    "p = mp.Pool(NCORES)\n",
    "analyte_less_stringent_model_list = p.starmap(lmm_model_2, arg_list)\n",
    "p.close()\n",
    "p.join()\n",
    "\n",
    "joblib.dump(analyte_less_stringent_model_list, os.path.join(OUTPUT_PATH, \"shortChain_fattyAcid_less_stringentModel.joblib\"))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bcf7fa8b",
   "metadata": {},
   "source": [
    "## LMM on MS patients who undergo treatment (along with their control pairs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "ed20fdf2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There 267 MS patients that undergo treatment\n",
      "Hence, total samples =  534\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "house_to_include = data_outlier_removed[data_outlier_removed.Treatment==\"Treated\"][\"house\"].unique()\n",
    "data_outlier_removed_ms_treated = data_outlier_removed[data_outlier_removed.house.isin(house_to_include)]\n",
    "print(\"There are {} MS patients that undergo treatment\".format(house_to_include.shape[0]))\n",
    "print(\"Hence, total samples = \", data_outlier_removed_ms_treated[\"Unique Tube Label ID\"].unique().shape[0])\n",
    "\n",
    "data_list = [data_outlier_removed_ms_treated]*len(analytes)\n",
    "analytes_list = list(analytes)\n",
    "maxiter_list = [10000]*len(analytes)\n",
    "arg_list = list(zip(data_list, analytes_list, maxiter_list))\n",
    "p = mp.Pool(NCORES)\n",
    "analyte_stringent_model_ms_treated_list = p.starmap(lmm_model_2, arg_list)\n",
    "p.close()\n",
    "p.join()\n",
    "\n",
    "\n",
    "joblib.dump(analyte_stringent_model_ms_treated_list, os.path.join(OUTPUT_PATH, \"shortChain_fattyAcid_stringentModel_MS_treated.joblib\"))\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77572d2a",
   "metadata": {},
   "source": [
    "## LMM on MS patients that are not undergoing treatment (along with their control pairs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "95e6b48c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are 231 MS patients that do not undergo treatment\n",
      "Hence, total samples =  462\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "house_to_include = data_outlier_removed[data_outlier_removed.Treatment==\"Off\"][\"house\"].unique()\n",
    "data_outlier_removed_ms_not_treated = data_outlier_removed[data_outlier_removed.house.isin(house_to_include)]\n",
    "print(\"There are {} MS patients that do not undergo treatment\".format(house_to_include.shape[0]))\n",
    "print(\"Hence, total samples = \", data_outlier_removed_ms_not_treated[\"Unique Tube Label ID\"].unique().shape[0])\n",
    "\n",
    "data_list = [data_outlier_removed_ms_not_treated]*len(analytes)\n",
    "analytes_list = list(analytes)\n",
    "maxiter_list = [10000]*len(analytes)\n",
    "arg_list = list(zip(data_list, analytes_list, maxiter_list))\n",
    "p = mp.Pool(NCORES)\n",
    "analyte_stringent_model_ms_not_treated_list = p.starmap(lmm_model_2, arg_list)\n",
    "p.close()\n",
    "p.join()\n",
    "    \n",
    "joblib.dump(analyte_stringent_model_ms_not_treated_list, os.path.join(OUTPUT_PATH, \"shortChain_fattyAcid_stringentModel_MS_not_treated.joblib\"))\n",
    "        \n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de4a5a66",
   "metadata": {},
   "source": [
    "## LMM on MS patients that undergo treatment (with outliers included)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca3e8223",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "house_to_include = data_nan_removed[data_nan_removed.Treatment==\"Treated\"][\"house\"].unique()\n",
    "data_nan_removed_ms_treated = data_nan_removed[data_nan_removed.house.isin(house_to_include)]\n",
    "print(\"There are {} MS patients that undergo treatment\".format(house_to_include.shape[0]))\n",
    "print(\"Hence, total samples = \", data_nan_removed_ms_treated[\"Unique Tube Label ID\"].unique().shape[0])\n",
    "\n",
    "data_list = [data_nan_removed_ms_treated]*len(analytes)\n",
    "analytes_list = list(analytes)\n",
    "maxiter_list = [10000]*len(analytes)\n",
    "arg_list = list(zip(data_list, analytes_list, maxiter_list))\n",
    "p = mp.Pool(NCORES)\n",
    "analyte_less_stringent_model_ms_treated_list = p.starmap(lmm_model_2, arg_list)\n",
    "p.close()\n",
    "p.join()\n",
    "\n",
    "joblib.dump(analyte_less_stringent_model_ms_treated_list, os.path.join(OUTPUT_PATH, \"shortChain_fattyAcid_less_stringentModel_MS_treated.joblib\"))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c6d8d51",
   "metadata": {},
   "source": [
    "## LMM on MS patients that undergo treatment (with outliers included)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90d5c5f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "house_to_include = data_nan_removed[data_nan_removed.Treatment==\"Off\"][\"house\"].unique()\n",
    "data_nan_removed_ms_not_treated = data_nan_removed[data_nan_removed.house.isin(house_to_include)]\n",
    "print(\"There are {} MS patients that do not undergo treatment\".format(house_to_include.shape[0]))\n",
    "print(\"Hence, total samples = \", data_nan_removed_ms_not_treated[\"Unique Tube Label ID\"].unique().shape[0])\n",
    "\n",
    "data_list = [data_nan_removed_ms_not_treated]*len(analytes)\n",
    "analytes_list = list(analytes)\n",
    "maxiter_list = [10000]*len(analytes)\n",
    "arg_list = list(zip(data_list, analytes_list, maxiter_list))\n",
    "p = mp.Pool(NCORES)\n",
    "analyte_less_stringent_model_ms_not_treated_list = p.starmap(lmm_model_2, arg_list)\n",
    "p.close()\n",
    "p.join()\n",
    "\n",
    "joblib.dump(analyte_less_stringent_model_ms_not_treated_list, os.path.join(OUTPUT_PATH, \"shortChain_fattyAcid_less_stringentModel_MS_not_treated.joblib\"))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7219d971",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb1d12bc",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c22960d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51673874",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "id": "b072a883",
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Get predicted values and residuals from the mixed effects model\n",
    "# predicted_values = mdf.predict()\n",
    "# residuals = mdf.resid\n",
    "\n",
    "# # Create the plot\n",
    "# plt.scatter(predicted_values, residuals, alpha=0.5)\n",
    "# plt.xlabel(\"Predicted Values\")\n",
    "# plt.ylabel(\"Residuals\")\n",
    "# plt.title(\"Residual vs Predicted Value Plot\")\n",
    "# # plt.axhline(y=0, color='r', linestyle='--') # Add a horizontal line at y=0 for reference\n",
    "# plt.show()\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# def get_data(data_type, sample):\n",
    "#     data_chemical_annotation = []\n",
    "#     data_sample_metadata = []\n",
    "#     if data_type == \"targeted\":\n",
    "#         filename = SHORT_CHAIN_FATTY_ACID_DATA_FILENAME\n",
    "#         mapping_filename = \"short_chain_fatty_acid_spoke_map.csv\"\n",
    "#         file_path = os.path.join(DATA_ROOT_PATH, filename)\n",
    "#         mapping_filepath = os.path.join(OUTPUT_PATH, mapping_filename)\n",
    "#         data = pd.read_excel(file_path, engine='openpyxl')\n",
    "#         data = data[data[\"Client Matrix\"]==sample]\n",
    "# #         data = data[data[\"Analysis Comment\"].isna()]\n",
    "#         mapping_data = pd.read_csv(mapping_filepath)\n",
    "#         analytes = mapping_data[\"name\"].unique()\n",
    "#     elif data_type == \"global\":\n",
    "#         sheet_name = [\"Chemical Annotation\", \"Sample Meta Data\", \"Log Transformed Data\"]\n",
    "#         mapping_filename = \"global_metabolomics_compound_spoke_map.csv\"\n",
    "#         if sample == \"serum\":\n",
    "#             filename = GLOBAL_SERUM_DATA_FILENAME\n",
    "#         else:\n",
    "#             filename = GLOBAL_STOOL_DATA_FILENAME\n",
    "#         file_path = os.path.join(DATA_ROOT_PATH, filename)\n",
    "#         mapping_filepath = os.path.join(OUTPUT_PATH, mapping_filename)\n",
    "#         data_chemical_annotation = pd.read_excel(file_path, engine='openpyxl', sheet_name=sheet_name[0])\n",
    "#         data_sample_metadata = pd.read_excel(file_path, engine='openpyxl', sheet_name=sheet_name[1])\n",
    "#         data = pd.read_excel(file_path, engine='openpyxl', sheet_name=sheet_name[2])\n",
    "#         mapping_data = pd.read_csv(mapping_filepath)\n",
    "#         analytes = mapping_data[\"CHEM_ID\"].unique()\n",
    "#     return data, mapping_data, analytes, data_chemical_annotation, data_sample_metadata\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0b83b59",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}