GCIMS_gc-ims-tools .ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3dc8aae8",
   "metadata": {},
   "source": [
    "# GCIMS (R package) and gc-ims-tools(Python Package) Comparison for pre-processing urines sampled with GC-IMS\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79ffca1c",
   "metadata": {},
   "source": [
    "## IMPORT LIBRARIES"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "dd80f442",
   "metadata": {},
   "outputs": [],
   "source": [
    "import ims\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import random as rd\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.cross_decomposition import PLSRegression\n",
    "from sklearn import metrics\n",
    "from sklearn import decomposition\n",
    "import multiprocessing\n",
    "from joblib import Parallel, delayed"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0472a858",
   "metadata": {},
   "source": [
    "## LOAD DATASETS"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f7e5d2c",
   "metadata": {},
   "source": [
    "First the peak table obtained with the GCIMS R package is imported along with the annotations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "caafb0ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "dataR=pd.read_csv(\"peak_table_R.csv\",index_col=0)\n",
    "#dataR.index=range(dataR.shape[0])\n",
    "annotations=pd.read_csv(\"annotations.csv\",index_col=0)\n",
    "annotations[\"Sex_b\"]=[int(b) for b in annotations[\"Sex\"]==\"Male\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9f5dfec",
   "metadata": {},
   "source": [
    "Then the raw data obtained with the GC-IMS is imported and pre-processed with the gc-ims-tools Python package"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "4924aad9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Dataset: Urines, 56 Spectra"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ds=ims.Dataset.read_mea(\"Urines\")\n",
    "ds.binning(2).cut_dt(7, 17).cut_rt(0, 1300)\n",
    "ds.interp_riprel()\n",
    "ds.scaling()\n",
    "ds.tophat()\n",
    "X,y=ds.get_xy()\n",
    "dataP=pd.DataFrame(X)\n",
    "dataP.index=annotations.index"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e78955ec",
   "metadata": {},
   "source": [
    "## Get Predictions with optimized Latent Variables"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2bc791b2",
   "metadata": {},
   "source": [
    "Function to find the latent variables that produce the best AUC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "52719fbf",
   "metadata": {},
   "outputs": [],
   "source": [
    "def optim(Xtrain,annotations,j):\n",
    "    predL=pd.DataFrame()\n",
    "    for k in np.unique(annotations[\"ID\"]):\n",
    "        XtrainL=Xtrain.loc[annotations[\"ID\"]!=k].values.tolist()\n",
    "        XtestL=Xtrain.loc[annotations[\"ID\"]==k].values.tolist()\n",
    "        YtrainL=annotations[annotations[\"ID\"]!=k][\"Sex_b\"].values.tolist()\n",
    "        YtestL=annotations[annotations[\"ID\"]==k][\"Sex_b\"].values.tolist()\n",
    "        plsL=PLSRegression(n_components=j, scale=True)\n",
    "        plsL.fit(XtrainL,YtrainL)\n",
    "        L=pd.DataFrame(plsL.predict(XtestL))\n",
    "        L[\"Label\"]=YtestL\n",
    "        predL=pd.concat([predL,L])\n",
    "    auc=metrics.roc_auc_score(predL['Label'],predL[0])\n",
    "    return auc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "022279f3",
   "metadata": {},
   "source": [
    "Cross Validation with optimal number of latent variables for data preprocesed with GCIMS R package (peak table)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "aa8a463a",
   "metadata": {},
   "outputs": [],
   "source": [
    "predR=pd.DataFrame()\n",
    "lvhR=[]\n",
    "for i in np.unique(annotations[\"ID\"]):\n",
    "    Xtrain=dataR.loc[annotations[\"ID\"]!=i]\n",
    "    annotations_optim=annotations.loc[annotations[\"ID\"]!=i]\n",
    "    #Parallelization for optimizing latent variables\n",
    "    AUCL = Parallel(n_jobs=multiprocessing.cpu_count()-2)(delayed(optim)(Xtrain,annotations_optim,j) for j in range(1,7))\n",
    "    LV=AUCL.index(max(AUCL))+1\n",
    "    lvhR.append(LV)\n",
    "    Xtrain=dataR.loc[annotations[\"ID\"]!=i].values.tolist()\n",
    "    Xtest=dataR.loc[annotations[\"ID\"]==i].values.tolist()\n",
    "    Ytrain=annotations[annotations[\"ID\"]!=i][\"Sex_b\"].values.tolist()\n",
    "    Ytest=annotations[annotations[\"ID\"]==i][\"Sex_b\"].values.tolist()\n",
    "    pls=PLSRegression(n_components=LV, scale=True)\n",
    "    pls.fit(Xtrain,Ytrain)\n",
    "    A=pd.DataFrame(pls.predict(Xtest))\n",
    "    A[\"Label\"]=Ytest\n",
    "    predR=pd.concat([predR,A])   \n",
    "AUCR=metrics.roc_auc_score(predR['Label'],predR[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e3eb72d",
   "metadata": {},
   "source": [
    "Cross Validation with optimal number of latent variables for data preprocesed with gc-ims-tools for Python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "b79514fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "#This may take a while since it's working with the full matrix\n",
    "predP=pd.DataFrame()\n",
    "lvhP=[]\n",
    "for i in np.unique(annotations[\"ID\"]):\n",
    "    Xtrain=dataP.loc[annotations[\"ID\"]!=i]\n",
    "    annotations_optim=annotations.loc[annotations[\"ID\"]!=i]\n",
    "    #Parallelization for optimizing latent variables/ uncomment if your computer can handle it\n",
    "    #AUCLP = Parallel(n_jobs=multiprocessing.cpu_count()-4)(delayed(optim)(Xtrain,annotations_optim,j) for j in range(1,7)) \n",
    "    #Comment next 3 lines to use paralelization\n",
    "    AUCLP=[]\n",
    "    for j in range(1,7):\n",
    "        a=optim(Xtrain,annotations_optim,j)\n",
    "        AUCLP.append(a)\n",
    "    LV=AUCL.index(max(AUCL))+1\n",
    "    lvhP.append(LV)\n",
    "    Xtrain=dataP.loc[annotations[\"ID\"]!=i].values.tolist()\n",
    "    Xtest=dataP.loc[annotations[\"ID\"]==i].values.tolist()\n",
    "    Ytrain=annotations[annotations[\"ID\"]!=i][\"Sex_b\"].values.tolist()\n",
    "    Ytest=annotations[annotations[\"ID\"]==i][\"Sex_b\"].values.tolist()\n",
    "    pls=PLSRegression(n_components=LV, scale=True)\n",
    "    pls.fit(Xtrain,Ytrain)\n",
    "    A=pd.DataFrame(pls.predict(Xtest))\n",
    "    A[\"Label\"]=Ytest\n",
    "    predP=pd.concat([predP,A]) \n",
    "AUCP=metrics.roc_auc_score(predP['Label'],predP[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "845960db",
   "metadata": {},
   "source": [
    "## ROC COMPARISON"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75acf89d",
   "metadata": {},
   "source": [
    "Comparison of ROC curves obtained with both preprocessing tehcniques in the same dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "790ce636",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x24e825449a0>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAHFCAYAAAAOmtghAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABbhElEQVR4nO3deVxUVf8H8M/AsA2bIoKgiFu4Lwi5gOaSu2laKabmkmbkrqlpPYn6VD6aW1ZqpmKaC+WWPamJa7j9VMTcSE1JXEDDBRQEWb6/P4x5HGfQGRgYuHzer9e8Xsy55577vYeR+XruueeqRERAREREpBBWlg6AiIiIyJyY3BAREZGiMLkhIiIiRWFyQ0RERIrC5IaIiIgUhckNERERKQqTGyIiIlIUJjdERESkKExuiIiISFGY3BCVYCtXroRKpdK+1Go1vLy80KdPH1y8eNHgPpmZmVi8eDGaN28OV1dXODg4oHbt2pg8eTJu375tcJ+cnBysXr0a7dq1g7u7O2xsbODh4YFXXnkFP//8M3Jycp4ba0ZGBr766iu0aNECZcuWha2tLSpWrIjevXtj//79BeqHkuCvv/6CSqXCypUrLR0KkeIxuSFSgPDwcBw+fBi7du3CyJEjsXXrVrRo0QJ3797VqZeWlob27dtj1KhR8Pf3x7p167Bt2za89dZbWLp0Kfz9/XH+/HmdfdLT09GlSxcMHDgQHh4eWLx4Mfbs2YMlS5bA29sbvXr1ws8///zM+JKSkhAcHIzx48ejXr16WLlyJXbv3o25c+fC2toaL7/8Mn7//Xez90tx4uXlhcOHD6Nr166WDoVI+YSISqzw8HABIMeOHdMpnz59ugCQFStW6JQPGzZMAMj69ev12jp//ry4urpK3bp1JSsrS1v+3nvvCQD57rvvDMZw4cIF+f33358ZZ+fOnUWtVsvu3bsNbj969KhcuXLlmW2UVFlZWZKenm7pMIhKFY7cEClQYGAgAODmzZvassTERKxYsQIdO3ZESEiI3j5+fn744IMPcPbsWWzZskW7z7Jly9CxY0cMGDDA4LFeeOEFNGjQIM9YoqOjsX37dgwZMgRt27Y1WOfFF19E5cqVte/PnDmDV199FWXLloW9vT0aNWqE7777Tmefffv2QaVSYe3atfjggw/g5eUFJycndOvWDTdv3sT9+/cxbNgwuLu7w93dHYMHD8aDBw902lCpVBg5ciS++eYb+Pn5wc7ODnXq1MH69et16v39998YPnw46tSpAycnJ3h4eKBt27aIiorSqZd76Wn27Nn45JNPULVqVdjZ2WHv3r0GL0v9/fffGDZsGHx8fGBnZ4fy5csjODgYu3bt0ml3xYoVaNiwIezt7eHm5oaePXsiNjZWp86gQYPg5OSEP//8E126dIGTkxN8fHzw/vvvIyMjI8/fD5ESMbkhUqC4uDgAjxOWXHv37kVWVhZ69OiR53652yIjI7X7ZGZmPnOf59m5c6dO289z/vx5BAUF4ezZs1i4cCE2bdqEOnXqYNCgQZg9e7Ze/Q8//BC3bt3CypUrMXfuXOzbtw9vvvkmXn/9dbi6umLdunWYNGkSVq9ejQ8//FBv/61bt2LhwoWYMWMGNmzYAF9fX7z55pvYsGGDts6dO3cAAGFhYfjll18QHh6OatWqoXXr1ti3b59emwsXLsSePXswZ84cbN++HbVq1TJ4rm+99Ra2bNmCqVOnYufOnVi2bBnatWunM/dp5syZGDJkCOrWrYtNmzbhiy++wKlTp9C8eXO9eVWZmZno3r07Xn75Zfz00094++23MX/+fMyaNcuovidSDEsPHRFR/uVeljpy5IhkZmbK/fv3ZceOHVKhQgV56aWXJDMzU1v3P//5jwCQHTt25Nnew4cPBYB07tzZ6H2eJzQ0VADIH3/8YVT9Pn36iJ2dncTHx+uUd+7cWTQajdy7d09ERPbu3SsApFu3bjr1xo4dKwBk9OjROuU9evQQNzc3nTIA4uDgIImJidqyrKwsqVWrltSoUSPPGLOysiQzM1Nefvll6dmzp7Y8Li5OAEj16tXl0aNHOvvkbgsPD9eWOTk5ydixY/M8zt27d8XBwUG6dOmiUx4fHy92dnbSt29fbdnAgQMFgPzwww86dbt06SI1a9bM8xhESsSRGyIFaNasGWxsbODs7IxOnTqhbNmy+Omnn6BWq/PVnkqlMnOExtuzZw9efvll+Pj46JQPGjQIaWlpOHz4sE75K6+8ovO+du3aAKA3cbd27dq4c+eO3qWpl19+GZ6entr31tbWCAkJwZ9//olr165py5csWYLGjRvD3t4earUaNjY22L17t97lIQDo3r07bGxsnnuuTZo0wcqVK/HJJ5/gyJEjyMzM1Nl++PBhPHz4EIMGDdIp9/HxQdu2bbF7926dcpVKhW7duumUNWjQAFeuXHluLERKwuSGSAFWrVqFY8eOYc+ePXj33XcRGxuLN998U6dO7pyW3EtWhuRuy00sjNnneUxt4/bt2/Dy8tIr9/b21m5/kpubm857W1vbZ5anp6frlFeoUEHvWLlluceaN28e3nvvPTRt2hQbN27EkSNHcOzYMXTq1AkPHz7U299Q/IZERERg4MCBWLZsGZo3bw43NzcMGDAAiYmJOsfPqz+e7guNRgN7e3udMjs7O71zJlI6JjdEClC7dm0EBgaiTZs2WLJkCYYOHYodO3bozBtp06YN1Gq1drKwIbnb2rdvr93Hxsbmmfs8T8eOHXXafp5y5cohISFBr/zGjRsAAHd393zHYkhuImGorFy5cgCA77//Hq1bt8bixYvRtWtXNG3aFIGBgbh//77BNo0d+XJ3d8eCBQvw119/4cqVK5g5cyY2bdqkHanJPX5e/WHuviBSCiY3RAo0e/ZslC1bFlOnTtUusFehQgW8/fbb+PXXXxEREaG3z4ULFzBr1izUrVtXO/m3QoUKGDp0KH799VesWrXK4LEuXbqEU6dO5RlL48aN0blzZyxfvhx79uwxWOf48eOIj48H8Pgy0Z49e7TJTK5Vq1ZBo9GgWbNmzz1/U+zevVvnrrLs7GxERESgevXqqFSpEoDHyYqdnZ3OfqdOndK7RFYQlStXxsiRI9G+fXucOHECANC8eXM4ODjg+++/16l77do17eU7ItKXvwvyRFSslS1bFlOmTMGkSZOwdu1a9O/fH8Djyyvnz59H//798dtvv6Fbt26ws7PDkSNHMGfOHDg7O2Pjxo2wtrbWtjVv3jxcvnwZgwYNwq+//oqePXvC09MTSUlJiIyMRHh4ONavX//M28FXrVqFTp06oXPnznj77bfRuXNnlC1bFgkJCfj555+xbt06REdHo3LlyggLC8N///tftGnTBlOnToWbmxvWrFmDX375BbNnz4arq6tZ+8rd3R1t27bFxx9/DEdHRyxatAh//PGHzu3gr7zyCv79738jLCwMrVq1wvnz5zFjxgxUrVoVWVlZ+TpucnIy2rRpg759+6JWrVpwdnbGsWPHsGPHDrz22msAgDJlyuDjjz/Ghx9+iAEDBuDNN9/E7du3MX36dNjb2yMsLMwsfUCkOJae0UxE+ZfXIn4ij+98qly5srzwwgs6i/I9evRIvv76a2natKk4OTmJnZ2d1KxZUyZNmiRJSUkGj5OVlSXfffedtG3bVtzc3EStVkv58uWlc+fOsnbtWsnOzn5urA8fPpSFCxdK8+bNxcXFRdRqtXh7e8trr70mv/zyi07d06dPS7du3cTV1VVsbW2lYcOGOncZifzvbqkff/zRqD4JCwsTAPL3339rywDIiBEjZNGiRVK9enWxsbGRWrVqyZo1a3T2zcjIkAkTJkjFihXF3t5eGjduLFu2bJGBAweKr6+vtl7uHVGff/653vk/fbdUenq6hIaGSoMGDcTFxUUcHBykZs2aEhYWJqmpqTr7Llu2TBo0aCC2trbi6uoqr776qpw9e1anzsCBA8XR0VHvuLnnTVSaqERELJdaERFZjkqlwogRI/DVV19ZOhQiMiPOuSEiIiJFYXJDREREisIJxURUavGqPJEyceSGiIiIFIXJDRERESkKkxsiIiJSlFI35yYnJwc3btyAs7OzRR8OSERERMYTEdy/fx/e3t6wsnr22EypS25u3Lih97RhIiIiKhmuXr2qfTRKXkpdcuPs7Azgcee4uLhYOBoiIiIyRkpKCnx8fLTf489S6pKb3EtRLi4uTG6IiIhKGGOmlHBCMRERESkKkxsiIiJSFCY3REREpCilbs6NsbKzs5GZmWnpMIiohLOxsYG1tbWlwyAqVZjcPEVEkJiYiHv37lk6FCJSiDJlyqBChQpcW4uoiDC5eUpuYuPh4QGNRsM/RkSUbyKCtLQ03Lp1CwDg5eVl4YiISgcmN0/Izs7WJjblypWzdDhEpAAODg4AgFu3bsHDw4OXqIiKACcUPyF3jo1Go7FwJESkJLl/UziPj6hoMLkxgJeiiMic+DeFqGgxuSEiIiJFsWhy89tvv6Fbt27w9vaGSqXCli1bnrvP/v37ERAQAHt7e1SrVg1Lliwp/EBLqZUrV6JMmTKWDsNspk2bhkaNGhWojeXLl6NDhw7mCegZiqrv//vf/8Lf3x85OTmFfiwioqJi0eQmNTUVDRs2xFdffWVU/bi4OHTp0gUtW7ZETEwMPvzwQ4wePRobN24s5EhLp5CQEFy4cKHIj2tsolvUMjIyMHXqVHz88cfasmnTpkGlUkGlUsHa2ho+Pj4YOnQo/v77b6PbrVKlChYsWFAIET/fK6+8ApVKhbVr1+a7jQ4dOsDa2hpHjhzR29a6dWuMHTtWr3zLli16l2oePXqE2bNno2HDhtBoNHB3d0dwcDDCw8PzPVclIyMDo0aNgru7OxwdHdG9e3dcu3btuftdv34d/fv3R7ly5aDRaNCoUSNER0cDeDxv5oMPPkD9+vXh6OgIb29vDBgwADdu3MhXjERkfha9W6pz587o3Lmz0fWXLFmCypUra78IateujePHj2POnDl4/fXXCynK0svBwUF7pwcBGzduhJOTE1q2bKlTXrduXezatQvZ2dmIiYnBkCFDcP36dWzfvt1CkZpm8ODB+PLLL9G/f3+T942Pj8fhw4cxcuRILF++HM2aNctXDI8ePULHjh3x+++/49///jeCg4Ph4uKCI0eOYM6cOfD398/XqNvYsWPx888/Y/369ShXrhzef/99vPLKK4iOjs7zrqW7d+8iODgYbdq0wfbt2+Hh4YFLly5pR9LS0tJw4sQJfPzxx2jYsCHu3r2LsWPHonv37jh+/Hi+zp8oP0SAtLTi2NhjGncNVFYWmm8mxQQA2bx58zPrtGzZUkaPHq1TtmnTJlGr1fLo0SOD+6Snp0tycrL2dfXqVQEgycnJenUfPnwo586dk4cPH+b7PCwlJSVF+vbtKxqNRipUqCDz5s2TVq1ayZgxY7R10tPTZeLEiVKpUiWxtbWVGjVqyLJly/JsMzw8XFxdXbXvw8LCpGHDhrJ8+XLx8fERR0dHCQ0NlaysLJk1a5Z4enpK+fLl5ZNPPtFpJywsTHx8fMTW1la8vLxk1KhReR7T19dXAGhfvr6+2m2LFi2SatWqiY2Njfj5+cmqVat09r1y5Yp0795dHB0dxdnZWXr16iWJiYl68efau3evvPjii6LRaMTV1VWCgoLkr7/+yjO2bt26yYQJE/TO7ck2RUQ++eQTsbKykrS0NGnTpo2MGDFCZ3tSUpLY2trK7t27pVWrVjrnm/tPMrfvd+zYIbVq1RJHR0fp2LGj3LhxQ9tOdna2TJ8+XSpWrCi2trbSsGFD2b59u3Z7XFycAJCNGzdK69atxcHBQRo0aCCHDh3Sieevv/4SAHLp0qU8zz0v06ZNkz59+khsbKw4OzvLgwcPdLY//RnMtXnzZnnyz8+sWbPEyspKTpw4oVf30aNHeu0a4969e2JjYyPr16/Xll2/fl2srKxkx44dee73wQcfSIsWLUw61tGjRwWAXLlyxeD2kvy3hYqnnByRoCCRx1lJ8Xw9uGn6v9tnSU5OzvP7+2klakJxYmIiPD09dco8PT2RlZWFpKQkg/vMnDkTrq6u2pePj49pBxUBUlMt8xIxOszx48fj4MGD2Lp1KyIjIxEVFYUTJ07o1BkwYADWr1+PhQsXIjY2FkuWLIGTk5NJ3XHp0iVs374dO3bswLp167BixQp07doV165dw/79+zFr1iz861//0l6i2LBhA+bPn49vvvkGFy9exJYtW1C/fv082z927BgAIDw8HAkJCdr3mzdvxpgxY/D+++/jzJkzePfddzF48GDs3bv3n1+ToEePHrhz5w7279+PyMhIXLp0CSEhIQaPk5WVhR49eqBVq1Y4deoUDh8+jGHDhj3zrpaoqCgEBgY+t48cHByQk5ODrKwsDB06FGvXrkVGRoZ2+5o1a+Dt7Y02bdpg06ZNqFSpEmbMmIGEhAQkJCRo66WlpWHOnDlYvXo1fvvtN8THx2PChAna7V988QXmzp2LOXPm4NSpU+jYsSO6d++Oixcv6sTz0UcfYcKECTh58iT8/Pzw5ptvIisrS7vd19cXHh4eiIqKeu65PUlEEB4ejv79+6NWrVrw8/PDDz/8YFIbT/ZJu3bt4O/vr7fNxsYGjo6OAIDPPvsMTk5Oz3zlnkd0dDQyMzN15kh5e3ujXr16OHToUJ6xbN26FYGBgejVqxc8PDzg7++Pb7/99pnxJycnQ6VSKWqOGhVvaWnAMz7GZNa0qgBgxMjNCy+8IJ999plO2YEDBwSAJCQkGNynwCM3Dx5YMO01LutNSUkRGxsb+fHHH7Vl9+7dE41Go/1f8/nz5wWAREZGGtWmiOGRG41GIykpKdqyjh07SpUqVSQ7O1tbVrNmTZk5c6aIiMydO1f8/PzyHFkzxNBnISgoSN555x2dsl69ekmXLl1ERGTnzp1ibW0t8fHx2u1nz54VAHL06FFt/LmjLLdv3xYAsm/fPqNiunv3rgCQ3377Taf86ZGb2NhYqVGjhjRp0kREHn/+3NzcJCIiQlunUaNGMm3aNO17X19fmT9/vk674eHhAkD+/PNPbdnXX38tnp6e2vfe3t7y6aef6uz34osvyvDhw0XkfyM3T47O5fZJbGyszn7+/v46MRlj586dUr58ecnMzBQRkfnz50twcLBOHWNHbhwcHPRGZQ25ffu2XLx48ZmvtLQ0ERFZs2aN2Nra6rXRvn17GTZsWJ7HsLOzEzs7O5kyZYqcOHFClixZIvb29vLdd98ZrP/w4UMJCAiQfv365dkmR27I3J78arp58/H7fL9uPpAH0Dx+Xb75+L0ZXjnZOWY9Z1NGbkrUCsUVKlRAYmKiTtmtW7egVqvzXFHYzs4OdnZ2RRGexVy+fBmZmZlo0qSJtszV1RU1a9bUvj958iSsra3RqlUrg23UrVsXV65cAQC0bNkyz/kiVapUgbOzs/a9p6cnrK2tYWVlpVOWu9x8r169sGDBAlSrVg2dOnVCly5d0K1bN6jVanz22Wf47LPPtPudO3cOlStXNnjc2NhYDBs2TKcsODgYX3zxhXa7j4+PzshcnTp1UKZMGcTGxuLFF1/U2dfNzQ2DBg1Cx44d0b59e7Rr1w69e/fOc3n8hw8fAgDs7e31tp0+fRpOTk7Izs5GRkYGWrdujaVLlwJ4/Pnr378/VqxYgd69e+PkyZP4/fffjZowrdFoUL16de17Ly8vbb+mpKTgxo0bCA4O1uuT33//XaesQYMGOm0Aj//d1KpVS1vu4OCANBOvty9fvhwhISFQqx//GXnzzTcxceJEnD9/XuezZwwRMWotGDc3N7i5uZnUtqnHysnJQWBgoPaz6e/vj7Nnz2Lx4sUYMGCATt3MzEz06dMHOTk5WLRoUYHiIsovR8fHr4L559+/h1kas7gSdVmqefPmiIyM1CnbuXMnAgMDYWNjUzgH1WiABw8s8zJypWT55/LV03+wc8sBPHdi8LZt23Dy5EmcPHkSy5Yty7Pe0/2sUqkMluXeWuzj44Pz58/j66+/hoODA4YPH46XXnoJmZmZCA0N1R7z5MmT8Pb2fmaMhs4vtyyvL6xnfZGFh4fj8OHDCAoKQkREBPz8/Aze8QMA5cqVg0qlwt27d/W21axZEydPnsS5c+fw8OFD7NmzBzVq1NBuHzp0KCIjI3Ht2jWsWLECL7/8Mnx9fZ95roDhvn7yd5pb9rzzfbKd3G1P3/p9584dlC9f/rkxPVl/y5YtWLRoEdRqNdRqNSpWrIisrCysWLFCW8/FxQXJycl6+9+7dw8uLi7a935+foiNjX3ucU25LFWhQgU8evRI73d269YtvcvbT/Ly8kKdOnV0ymrXro34+HidsszMTPTu3RtxcXGIjIzUOR8isiyLJjcPHjzQfrEBj2/1PnnypPaPyJQpU3T+pxQaGoorV65g/PjxiI2NxYoVK7B8+XKdeQhmp1L9Ly0u6peRq5pWr14dNjY2OHr0qLYsJSVFZ+5F/fr1kZOTg/379xtsw9fXFzVq1ECNGjVQsWLFgvXZUxwcHNC9e3csXLgQ+/btw+HDh3H69Gm4ublpj1mjRg3tCICNjQ2ys7N12qhduzYOHDigU3bo0CHUrl0bwONRmvj4eFy9elW7/dy5c0hOTtbWMcTf3x9TpkzBoUOHUK9evTxviba1tUWdOnVw7tw5g9tq1KiBqlWrGhwlrF+/PgIDA/Htt99i7dq1ePvtt/X2f/p8n8fFxQXe3t7P7BNjpaen49KlSwbnu+RlzZo1qFSpEn7//XedBHXBggX47rvvtHN6atWqZfAOomPHjumM7vTt2xe7du1CTEyMXt2srCykpqYCgF5CbOiVOy8qICAANjY2Ov8hSkhIwJkzZxAUFJTnuQUHB+P8+fM6ZRcuXNBJSHMTm4sXL2LXrl18Fh1RcWPWC2Im2rt3r96dIgBk4MCBIiIycOBAadWqlc4++/btE39/f7G1tZUqVarI4sWLTTrms67ZleTr4kOHDpWqVavKnj175MyZM/L666+Ls7OzjB07Vltn0KBB4uPjI5s3b5bLly/L3r17deaCPC2vu6WeNHDgQHn11Vd1yp6cZxEeHi7Lli2T06dPy6VLl+Sjjz4SBwcHSUpKyvO4L7zwgrz33nuSkJAgd+7cEZHHczRsbGxk8eLFcuHCBZk7d65YW1vL3r17RUQkJydH/P39pWXLlhIdHS3/93//JwEBATqfnyfjv3z5skyePFkOHTokf/31l/z666/i5uYmixYtyjOu8ePHy+uvv65TZqhPDFm6dKnY2tpKmTJl9D5f7du3l+7du8u1a9fk77//1vbbk32f2wdP/pOdP3++uLi4yPr16+WPP/6QDz74QGxsbOTChQsi8r85NzExMdp9cucO5fabyON/h05OTpKamvrc88jVsGFD+eCDD/TKU1JSxM7OTrZs2aKNwcHBQYYPHy4nT56U8+fPy1dffSV2dnbyww8/aPdLT0+Xli1bStmyZeWrr76SkydPyqVLlyQiIkIaN26scw6mCA0NlUqVKsmuXbvkxIkT0rZtW2nYsKFkZWVp67Rt21a+/PJL7fujR4+KWq2WTz/9VC5evChr1qwRjUYj33//vYiIZGZmSvfu3aVSpUpy8uRJSUhI0L4yMjIMxlGS/7ZQ8fTknJt83ExYiI0VHlPm3BSbCcVFRanJjaFbwZs0aSKTJ0/W1nn48KGMGzdOvLy8tLeCr1ixIs82zZHcbN68WZo2bSouLi7i6OgozZo1k127dj3zXLZu3So1atQQtVpdaLeCJyYmSo8ePbR94evrK1OnTtWZGP202NhYcXBwkHv37j2zTwy5f/++aDQa7WTfJx0+fFgaNGggdnZ2ereCP+np5ObJW8FtbGzyvBX8ecnNsGHD5N1339U5lqH/WOQ6fvy4zkTtp3Xr1k26deumU79jx47i4eEhLi4uEhgYKOvWrdPbLz09XWbOnCn169cXe3t7cXNzk+DgYFm5cqV20rKpHj58KCNHjhQ3NzdxcHCQV155RWfSucjjCd1hYWE6ZT///LPUq1dP7OzspFatWrJ06VLtttx+NfR6sl+fjqOk/m2h4onJzbOpREy431gBUlJS4OrqiuTkZL1r5Onp6YiLi0PVqlUNThwtSVJTU1GxYkXMnTsXQ4YMsXQ4itG7d2/tpSxTXL16FVWqVMGxY8fQuHHjQorOdH///bf20lHVqlW15a1bt0br1q0xbdo0ywWnIEr620LFQ2oqkLuSx4MHBZwDbNbGCs+zvr+fVqLulqK8xcTE4I8//kCTJk2QnJyMGTNmAABeffVVC0emLJ9//jm2bt1qdP3MzEwkJCRg8uTJaNasWbFKbIDH89wWLVqkk9jcv38fly5dwn//+18LRkakTOZaCPifaWgG3hS0MWVgcqMgc+bMwfnz52Fra4uAgABERUXB3d3d0mEpiq+vL0aNGmV0/YMHD6JNmzbw8/PDhg0bCjGy/GnSpInOEgIA4OzsrDMxm4jMQwRo0aIQFt/z9ID2Vm4CwORGMfz9/bUP9qPio3Xr1nq3bxNR6VQYqwoH4wA05kpsgoONXoKkuGNyQ0REVMRu3jTDPBlPD2iQBlWBG/uHRmP0EiTFHZMbIiKiImbWVYXN05iilKgViomIiIieh8kNERERKQqTGyIiIlIUJjdERESkKExuqNhZvnw5OnToYPJ+L774IjZt2lQIERERUUnC5EZBEhMTMWbMGNSoUQP29vbw9PREixYtsGTJEqQ9tSRmTEwMevXqBU9PT9jb28PPzw/vvPMOLly4AAD466+/oFKptE9sz32vVqtx/fp1nbYSEhKgVquhUqnw119/acs3btyIpk2bwtXVFc7Ozqhbty7ef//9Z55DRkYGpk6dio8//lhbNm3aNKhUKqhUKlhZWcHb2xv9+vXTW2ju448/xuTJk5GTk2Nq1wF4/ERta2trdOrUSW/bvn37oFKpcO/ePb1tjRo10ntMwfP6Nz/279+PgIAA2Nvbo1q1aliyZMkz669cuVLbb0+/bt26pVf/zz//hLOzM8qUKZPvGImIigMmNwpx+fJl+Pv7Y+fOnfjss88QExODXbt2Ydy4cfj555+xa9cubd3//ve/aNasGTIyMrBmzRrExsZi9erVcHV11UkqDPH29saqVat0yr777jtUrFhRp2zXrl3o06cP3njjDRw9ehTR0dH49NNP8ejRo2e2v3HjRjg5OaFly5Y65XXr1kVCQgKuXbuGiIgInD59Gr1799ap07VrVyQnJ+PXX3995jHysmLFCowaNQoHDhxAfHx8vtoACta/eYmLi0OXLl3QsmVLxMTE4MMPP8To0aOxcePGPPcJCQlBQkKCzqtjx45o1aoVPDw8dOpmZmbizTff1Ot3IvqHyOO1ZQr6ymXOtkhfYT/Fs7hR6lPBO3bsKJUqVZIHeTzRNScnR0REUlNTxd3dXXr06GGw3t27d0VE/4nSue//9a9/yQsvvKCzT82aNeXjjz8WABIXFyciImPGjJHWrVubfB7dunWTCRMm6JQZevL2woULDf4eBw0aJG+99ZbJx33w4IE4OzvLH3/8ISEhITJ9+nSd7Xv37hUA2v55UsOGDbVPlTa2f001adIkqVWrlk7Zu+++K82aNTO6jVu3bomNjY3e09Rz2+/fv7/Bp5FTwZXkvy0kIjk5IkFB/3tydj5fD6D538O3oSlweyXhSd7mZMpTwTly8xzmStbz8zJ21f7bt29j586dGDFiBBzzWMhJ9c+qk7/++iuSkpIwadIkg/Wed0mie/fuuHv3Lg4cOAAAOHDgAO7cuYNu3brp1KtQoQLOnj2LM2fOGHcS/4iKikJgYOAz6yQmJmLTpk2wtraGtbW1zrYmTZogKirKpGMCQEREBGrWrImaNWuif//+CA8Pz9djE0zpXycnp2e+OnfurK17+PBhvXlIHTt2xPHjx5GZmWlUbKtWrYJGo8Ebb7yhU75nzx78+OOP+Prrr408S6JSpjCem2AuCnpkgjlxheLnSEv735Pgi5qxT57/888/ISKoWbOmTrm7uzvS09MBACNGjMCsWbNw8eJFAECtWrXyFZONjQ369++PFStWoEWLFlixYgX69+8PGxsbnXqjRo1CVFQU6tevD19fXzRr1gwdOnRAv379YGdnZ7Dte/fu4d69e/D29tbbdvr0aTg5OSEnJwcPHz4EAIwePVovmatYsSLi4+ORk5MDKyvjc/fly5ejf//+AIBOnTrhwYMH2L17N9q1a2d0GwBM6t/c+Ux5cXBw0P6cmJgIT09Pne2enp7IyspCUlISvLy8nnu8FStWoG/fvjrt3r59G4MGDcL3338PFxeX57ZBVOoV5FEHqQBy/xnfvAWYY1FhBT0ywZyY3CiI6qkP+NGjR5GTk4N+/fohIyMDAMzyEMchQ4agefPm+Oyzz/Djjz/i8OHDyMrK0qnj6OiIX375BZcuXcLevXtx5MgRvP/++/jiiy9w+PBhaAz8TyM3abG3t9fbVrNmTWzduhUZGRn46aef8OOPP+LTTz/Vq+fg4ICcnBxkZGTofIk/y/nz53H06FHtnVZqtRohISFYsWKFycmNKf1bo0YNk9p++vebe6ynyw05fPgwzp07pzdf6p133kHfvn3x0ksvmRQLUallrkcdODqaJ7khg5jcPIdG83gExVLHNkaNGjWgUqnwxx9/6JRXq1YNgO4IgJ+fHwDgjz/+QPPmzfMVV7169VCrVi28+eabqF27NurVq5fnKET16tVRvXp1DB06FB999BH8/PwQERGBwYMH69UtV64cVCoV7t69q7fN1tZWmwzUrVsXFy9exHvvvYfVq1fr1Ltz5w40Go3RiQ3weNQmKytLZ1K0iMDGxgZ3795F2bJltaMaycnJepfu7t27B1dXVwCm9a/Tc4YEW7Zsie3btwN4fJkvMTFRZ/utW7egVqtRrly5557jsmXL0KhRIwQEBOiU79mzB1u3bsWcOXMAPD7vnJwcqNVqLF26FG+//fZz2yYiKm6Y3DyHSlX8n0dWrlw5tG/fHl999RVGjRqV57wbAOjQoQPc3d0xe/ZsbN68WW/7vXv3jLoV+O2338bw4cOxePFio+OsUqUKNBoNUvOY5W9ra4s6derg3Llzz13n5uOPP4afnx/GjRuHxo0ba8vPnDmj8/55srKysGrVKsydO1fvmK+//jrWrFmDkSNH4oUXXoCVlRWOHTsGX19fbZ2EhARcv35de0nQlP415bJU8+bN8fPPP+ts37lzJwIDA/UuCT7twYMH+OGHHzBz5ky9bYcPH0Z2drb2/U8//YRZs2bh0KFDenfAERGVGIU5s7k4UurdUn/++ad4enpKrVq1ZP369XLu3Dn5448/ZPXq1eLp6Snjx4/X1t2yZYvY2NhIt27dJDIyUuLi4uTYsWMyceJECQkJEZG875bKfZ+ZmSl///23ZGZmiohITEyMzt1SYWFhMnHiRNm7d69cvnxZTpw4IYMGDRIHBwf5448/8jyP8ePHy+uvv65TZuhuKRGR1157Tbp27apT1qpVK5kxY4bR/bZ582axtbWVe/fu6W378MMPpVGjRtr37733nlSuXFk2b94sly9flgMHDkirVq2kfv362n4QMa5/TXX58mXRaDQybtw4OXfunCxfvlxsbGxkw4YN2jqbNm2SmjVr6u27bNkysbe3lzt37jz3OLxbqnCU5L8tJI/vRjLDnUlmaqbUMuVuKSY3Tyjpf4Bu3LghI0eOlKpVq4qNjY04OTlJkyZN5PPPP5fU1FSduseOHZPXXntNypcvL3Z2dlKjRg0ZNmyYXLx4UUSen9w87enkZs+ePfL666+Lj4+P2Nraiqenp3Tq1EmioqKeeQ6xsbHi4OCgk2zkldwcPHhQAMiRI0dEROTatWtiY2MjV69e1dbJjXvv3r0Gj/fKK69Ily5dDG6Ljo4WABIdHS0iIunp6TJjxgypXbu2ODg4iK+vrwwaNEgSEhL09n1e/+bHvn37xN/fX2xtbaVKlSqyePFine3h4eFi6P8rzZs3l759+xp1DCY3haOk/20p9ZjcFAumJDcqETPMMC1BUlJS4OrqiuTkZL27Q9LT0xEXF4eqVasanNRKRaN3797w9/fHlClTTNpv4sSJSE5OxtKlS7Vl+/btQ8+ePXH58mWULVvW3KESGYV/W0q41NT/3TZr7G2shddMqfWs7++ncc4NFTuff/45tm7davJ+Hh4emDBhgk7Zjh078OGHHzKxISpMIo/XzVCq1FQIgDRoHt/Onf9mqIhw5OYJ/N8VERUGRf9tEQFatCi+i9yZgQBogQM4hGCztcmRG9OZMnLDFYqJiCj/ivPqvWaSBo1ZExsuKlz4eFmKiIjMoyCr9xZnT6wsbI5T5KLChY/JjQGl7EodERWyUvM3xVyr9xZjpeAUFYGXpZ6QuxhampInxhFRkcv9m/K8BReJyDw4cvMEa2trlClTBrdu3QIAaDQao57bQ0RkiIggLS0Nt27dQpkyZfSeYk9EhYPJzVMqVKgAANoEh4iooMqUKaP920JEhY/JzVNUKhW8vLzg4eGBzMxMS4dDRCWcjY0NR2yIihiTmzxYW1vzDxIREVEJxAnFREREpChMboiIiEhRmNwQERGRojC5ISIiIkVhckNERESKwuSGiIiIFIXJDRERESkKkxsiIiJSFCY3REREpChMboiIiEhRmNwQERGRojC5ISIiIkVhckNERESKwuSGiIiIFIXJDRERESkKkxsiIiJSFCY3REREpChMboiIiEhRmNwQERGRojC5ISIiIkVhckNERESKwuSGiIiIFEVt6QCIiMhCRIC0tIK1kZoKAZAGDZBqlqiKnVSFnpeSMbkhIiqNRIAWLYBDhwrWDIAWOIBDCAY8zRMaUUFZ/LLUokWLULVqVdjb2yMgIABRUVHPrL9mzRo0bNgQGo0GXl5eGDx4MG7fvl1E0RIRKURaWoETG+DxiM0hBJshoOIvOBjQaCwdBRnDoiM3ERERGDt2LBYtWoTg4GB888036Ny5M86dO4fKlSvr1T9w4AAGDBiA+fPno1u3brh+/TpCQ0MxdOhQbN682QJnQESkADdvAo6O+ds3FdoRm4I0UxJoNIBKZekoyBgqERFLHbxp06Zo3LgxFi9erC2rXbs2evTogZkzZ+rVnzNnDhYvXoxLly5py7788kvMnj0bV69eNeqYKSkpcHV1RXJyMlxcXAp+EkREJVFqKuDk9PjnBw/ynZWYqRmi5zLl+9til6UePXqE6OhodOjQQae8Q4cOOJTHUGlQUBCuXbuGbdu2QURw8+ZNbNiwAV27ds3zOBkZGUhJSdF5ERERkXJZLLlJSkpCdnY2PD11Z6B5enoiMTHR4D5BQUFYs2YNQkJCYGtriwoVKqBMmTL48ssv8zzOzJkz4erqqn35+PiY9TyIiIioeLH4hGLVUxcwRUSvLNe5c+cwevRoTJ06FdHR0dixYwfi4uIQGhqaZ/tTpkxBcnKy9mXs5SsiIiIqmSw2odjd3R3W1tZ6ozS3bt3SG83JNXPmTAQHB2PixIkAgAYNGsDR0REtW7bEJ598Ai8vL7197OzsYGdnZ/4TICIiomLJYiM3tra2CAgIQGRkpE55ZGQkgoKCDO6TlpYGKyvdkK2trQE8HvEhIiIisuhlqfHjx2PZsmVYsWIFYmNjMW7cOMTHx2svM02ZMgUDBgzQ1u/WrRs2bdqExYsX4/Llyzh48CBGjx6NJk2awNvb21KnQURERMWIRde5CQkJwe3btzFjxgwkJCSgXr162LZtG3x9fQEACQkJiI+P19YfNGgQ7t+/j6+++grvv/8+ypQpg7Zt22LWrFmWOgUiIiIqZiy6zo0lcJ0bIiJwnRsqcUrEOjdEREREhYHJDRERESkKkxsiIiJSFCY3REREpChMboiIiEhRmNwQERGRojC5ISIiIkVhckNERESKYtEViomIyDQiQFqaGRpKBQDNEz/ns5kC7EtUWJjcEBGVECJAixbAoUPmaM0R2qzG0xztERUfvCxFRFRCpKWZK7Exv+BgQKOxdBREj3HkhoioBLp5s4DPcUpNBTw9/mnsVoEfCqXRACpVgZogMhsmN0REJZCjozkeUvnP5B3Hf15ECsHLUkRERKQoTG6IiIhIUZjcEBERkaIwuSEiIiJFYXJDREREisK7pYieYLbVX6nkK4YfhserATs++aagjREpEpMbon+Yd/VXKvlUKNb3R3t6QHsrNxHp4GUpon8U59VfiZ4UjAPQmCux4dLCpEAcuSEyoMCrv1LJ9uTqvZfjit2HQaPxh0r1wFyNcWlhUhwmN0QGmGf1VyrZ/hkZ8eCHgaik4WUpIiIiUhQmN0RERKQoTG6IiIhIUZjcEBERkaIwuSEiIiJFYXJDREREisJbwYnI8orbow74aAKiEo3JDRFZFp97QURmxstSRGRZxfm5F3w0AVGJxJEbIio+ittzL/hoAqISKV/JTVZWFvbt24dLly6hb9++cHZ2xo0bN+Di4gInJydzx0hEpQWfe0FEZmBycnPlyhV06tQJ8fHxyMjIQPv27eHs7IzZs2cjPT0dS5YsKYw4iYiIiIxi8pybMWPGIDAwEHfv3oWDg4O2vGfPnti9e7dZgyMiIiIylckjNwcOHMDBgwdha2urU+7r64vr16+bLTAiIiKi/DB55CYnJwfZ2dl65deuXYOzs7NZgiIiIiLKL5OTm/bt22PBggXa9yqVCg8ePEBYWBi6dOliztiIiIiITGbyZan58+ejTZs2qFOnDtLT09G3b19cvHgR7u7uWLduXWHESApU3BakBbgoLRGRUpic3Hh7e+PkyZNYv349oqOjkZOTgyFDhqBfv346E4yJ8sIFaYmIqDCpRERM2eG3335DUFAQ1GrdvCgrKwuHDh3CSy+9ZNYAzS0lJQWurq5ITk6Gi4uLpcMplVJTgeK8HFJwMBAVxbXbisyTH4gHD7jODREZZMr3t8kjN23atEFCQgI8PDx0ypOTk9GmTRuDk42J8lLcFqQFuCgtEVFJZ3JyIyJQGfjLf/v2bTgWt28pKva4IC0REZmb0cnNa6+9BuDx3VGDBg2CnZ2ddlt2djZOnTqFoKAg80dIREREZAKjkxtXV1cAj0dunJ2ddSYP29raolmzZnjnnXfMHyERERGRCYxObsLDwwEAVapUwYQJE3gJioiIiIolk+fchIWFFUYcRERERGZhcnIDABs2bMAPP/yA+Ph4PHr0SGfbiRMnzBIYERERUX6Y/PiFhQsXYvDgwfDw8EBMTAyaNGmCcuXK4fLly+jcuXNhxEhERERkNJOTm0WLFmHp0qX46quvYGtri0mTJiEyMhKjR49GcnJyYcRIREREZDSTk5v4+HjtLd8ODg64f/8+AOCtt97is6WIiIjI4kxObipUqIDbt28DAHx9fXHkyBEAQFxcHEx8kgMRERGR2Zmc3LRt2xY///wzAGDIkCEYN24c2rdvj5CQEPTs2dPsARIRERGZwuS7pZYuXYqcnBwAQGhoKNzc3HDgwAF069YNoaGhZg+QiIiIyBQmJzdWVlawsvrfgE/v3r3Ru3dvAMD169dRsWJF80VHREREZCKTL0sZkpiYiFGjRqFGjRom77to0SJUrVoV9vb2CAgIQFRU1DPrZ2Rk4KOPPoKvry/s7OxQvXp1rFixIr+hExERkcIYndzcu3cP/fr1Q/ny5eHt7Y2FCxciJycHU6dORbVq1XDkyBGTk4yIiAiMHTsWH330EWJiYtCyZUt07twZ8fHxee7Tu3dv7N69G8uXL8f58+exbt061KpVy6TjEhERkXKpxMhbnIYPH46ff/4ZISEh2LFjB2JjY9GxY0ekp6cjLCwMrVq1MvngTZs2RePGjbF48WJtWe3atdGjRw/MnDlTr/6OHTvQp08fXL58GW5ubiYfDwBSUlLg6uqK5ORkuLi45KsNKpjUVMDJ6fHPDx4AfExZKccPBBEZwZTvb6NHbn755ReEh4djzpw52Lp1K0QEfn5+2LNnT74Sm0ePHiE6OhodOnTQKe/QoQMOHTpkcJ+tW7ciMDAQs2fPRsWKFeHn54cJEybg4cOHeR4nIyMDKSkpOi8iIiJSLqMnFN+4cQN16tQBAFSrVg329vYYOnRovg+clJSE7OxseHp66pR7enoiMTHR4D6XL1/GgQMHYG9vj82bNyMpKQnDhw/HnTt38rwkNnPmTEyfPj3fcdL/iABpaQVvJzW14G0QERHlxejkJicnBzY2Ntr31tbWcDTD8LFKpdJ5LyJ6ZU/GoFKpsGbNGri6ugIA5s2bhzfeeANff/01HBwc9PaZMmUKxo8fr32fkpICHx+fAsdd2ogALVoAeQyqERERFRtGJzcigkGDBsHOzg4AkJ6ejtDQUL0EZ9OmTUa15+7uDmtra71Rmlu3bumN5uTy8vJCxYoVtYkN8HiOjojg2rVreOGFF/T2sbOz08ZM+ZeWZv7EJjgY0GjM2yYREZHRyc3AgQN13vfv379AB7a1tUVAQAAiIyN1VjaOjIzEq6++anCf4OBg/Pjjj3jw4AGc/pmAeOHCBVhZWaFSpUoFioeMd/OmeeZ8ajRAHoN0RERE+Wb03VKFISIiAm+99RaWLFmC5s2bY+nSpfj2229x9uxZ+Pr6YsqUKbh+/TpWrVoFAHjw4AFq166NZs2aYfr06UhKSsLQoUPRqlUrfPvtt0Ydk3dL5Q9vaKFCww8XERnBlO9vk1coNqeQkBDcvn0bM2bMQEJCAurVq4dt27bB19cXAJCQkKCz5o2TkxMiIyMxatQoBAYGoly5cujduzc++eQTS50CERERFTMWHbmxBI7c5A//c02Fhh8uIjJCoaxzQ0RERFQSMLkhIiIiRWFyQ0RERIqSr+Rm9erVCA4Ohre3N65cuQIAWLBgAX766SezBkdERERkKpPvllq8eDGmTp2KsWPH4tNPP0V2djYAoEyZMliwYEGea9SQacz1qANz4SMTSA+fx0FExZTJyc2XX36Jb7/9Fj169MB//vMfbXlgYCAmTJhg1uBKKz7qgIo9fkiJqBgz+bJUXFwc/P399crt7OyQyv+BmUVhPOrAXPjIBALA53EQUbFm8shN1apVcfLkSe1Ce7m2b9+ufWo4mY+5HnVgLnxkAunh8ziIqJgxObmZOHEiRowYgfT0dIgIjh49inXr1mHmzJlYtmxZYcRYqjk6Fq/khkgPP6REVMyYnNwMHjwYWVlZmDRpEtLS0tC3b19UrFgRX3zxBfr06VMYMRIREREZrUCPX0hKSkJOTg48PDzMGVOhKgmPX+Bq9FTs8UNKREWsUB+/MH36dFy6dAkA4O7uXqISGyIiIlI+k5ObjRs3ws/PD82aNcNXX32Fv//+uzDiIiIiIsoXk5ObU6dO4dSpU2jbti3mzZuHihUrokuXLli7di3SitOqc0RERFQqFWjODQAcPHgQa9euxY8//oj09HSkpKSYK7ZCwTk3VKqZc1VhT8/HP/NDSkRFwJTvb5Pvlnqao6MjHBwcYGtri/v37xe0OSIqLFxVmIhKiXw9ODMuLg6ffvop6tSpg8DAQJw4cQLTpk1DYmKiueMjInPhqsJEVEqYPHLTvHlzHD16FPXr18fgwYO169wQUQnCVYWJSMFMTm7atGmDZcuWoW7duoURDxEVBa4qTEQKZnJy89lnnxVGHERERERmYVRyM378ePz73/+Go6Mjxo8f/8y68+bNM0tgRERERPlhVHITExODzMxM7c9ERERExZVRyc3evXsN/kxERERU3Jh8K/jbb79tcD2b1NRUvP3222YJioiIiCi/TE5uvvvuOzx8+FCv/OHDh1i1apVZgiIiIiLKL6PvlkpJSYGIQERw//592Nvba7dlZ2dj27ZtfEI4ERERWZzRyU2ZMmWgUqmgUqng5+ent12lUmH69OlmDY6IiIjIVEYnN3v37oWIoG3btti4cSPc3Ny022xtbeHr6wtvb+9CCZKIiIjIWEYnN61atQLw+LlSlStXhopLrhMREVExZFRyc+rUKdSrVw9WVlZITk7G6dOn86zboEEDswVHREREZCqjkptGjRohMTERHh4eaNSoEVQqFUREr55KpUJ2drbZgyQiIiIyllHJTVxcHMqXL6/9mYiIiKi4Miq58fX1NfgzERERUXGTr0X8fvnlF+37SZMmoUyZMggKCsKVK1fMGhwRERGRqUxObj777DM4ODgAAA4fPoyvvvoKs2fPhru7O8aNG2f2AImIiIhMYfSt4LmuXr2KGjVqAAC2bNmCN954A8OGDUNwcDBat25t7viISARISyt4O6mpBW+DiKgEMDm5cXJywu3bt1G5cmXs3LlTO1pjb29v8JlTRFQAIkCLFsChQ5aOhIioxDA5uWnfvj2GDh0Kf39/XLhwAV27dgUAnD17FlWqVDF3fESlW1qa+ROb4GBAozFvm0RExYjJyc3XX3+Nf/3rX7h69So2btyIcuXKAQCio6Px5ptvmj1AIvrHzZuAo2PB29FoAK4wTkQKphJDq/EpWEpKClxdXZGcnAwXFxdLh2NQairg5PT45wcPzPN9RiUUPwxERABM+/42eeQGAO7du4fly5cjNjYWKpUKtWvXxpAhQ+Dq6pqvgImIiIjMxeRbwY8fP47q1atj/vz5uHPnDpKSkjB//nxUr14dJ06cKIwYiYiIiIxm8sjNuHHj0L17d3z77bdQqx/vnpWVhaFDh2Ls2LH47bffzB4kERERkbFMTm6OHz+uk9gAgFqtxqRJkxAYGGjW4IiIiIhMZfJlKRcXF8THx+uVX716Fc7OzmYJioiIiCi/TE5uQkJCMGTIEERERODq1au4du0a1q9fj6FDh/JWcCIiIrI4ky9LzZkzByqVCgMGDEBWVhYAwMbGBu+99x7+85//mD1AIiIiIlPke52btLQ0XLp0CSKCGjVqQFNCVjzlOjdUovDDQEQEwLTvb6MvS6WlpWHEiBGoWLEiPDw8MHToUHh5eaFBgwYlJrEhIiIi5TM6uQkLC8PKlSvRtWtX9OnTB5GRkXjvvfcKMzYiIiIikxk952bTpk1Yvnw5+vTpAwDo378/goODkZ2dDWtr60ILkIiIiMgURo/cXL16FS1bttS+b9KkCdRqNW7cuFEogRERERHlh9HJTXZ2NmxtbXXK1Gq19o4pIiIiouLA6MtSIoJBgwbBzs5OW5aeno7Q0FA4PnEHx6ZNm8wbIREREZEJjE5uBg4cqFfWv39/swZDREREVFBGJzfh4eGFGQcRERGRWZj8+AVzW7RoEapWrQp7e3sEBAQgKirKqP0OHjwItVqNRo0aFW6AVLqIPF44rzi9iIjIJCY/fsGcIiIiMHbsWCxatAjBwcH45ptv0LlzZ5w7dw6VK1fOc7/k5GQMGDAAL7/8Mm7evFmEEZOiiQAtWgCHDlk6EiIiKgCLjtzMmzcPQ4YMwdChQ1G7dm0sWLAAPj4+WLx48TP3e/fdd9G3b180b968iCKlUiEtrfgmNsHBAFcCJyIyisVGbh49eoTo6GhMnjxZp7xDhw449IwvmPDwcFy6dAnff/89Pvnkk8IOk0qrmzeL13OcNBpApbJ0FEREJYLFkpukpCRkZ2fD09NTp9zT0xOJiYkG97l48SImT56MqKgoqNXGhZ6RkYGMjAzt+5SUlPwHTaWHo2PxSm6IiMho+bostXr1agQHB8Pb2xtXrlwBACxYsAA//fSTyW2pnvrfqIjolQGPFxHs27cvpk+fDj8/P6PbnzlzJlxdXbUvHx8fk2MkIiKiksPk5Gbx4sUYP348unTpgnv37iE7OxsAUKZMGSxYsMDodtzd3WFtba03SnPr1i290RwAuH//Po4fP46RI0dCrVZDrVZjxowZ+P3336FWq7Fnzx6Dx5kyZQqSk5O1r6tXrxp/skRERFTimJzcfPnll/j222/x0Ucf6TwwMzAwEKdPnza6HVtbWwQEBCAyMlKnPDIyEkFBQXr1XVxccPr0aZw8eVL7Cg0NRc2aNXHy5Ek0bdrU4HHs7Ozg4uKi8yIiIiLlMnnOTVxcHPz9/fXK7ezskGrimhzjx4/HW2+9hcDAQDRv3hxLly5FfHw8QkNDATwedbl+/TpWrVoFKysr1KtXT2d/Dw8P2Nvb65UTERFR6WVyclO1alWcPHkSvr6+OuXbt29HnTp1TGorJCQEt2/fxowZM5CQkIB69eph27Zt2rYTEhIQHx9vaohERERUiqlEREzZITw8HB9//DHmzp2LIUOGYNmyZbh06RJmzpyJZcuWoU+fPoUVq1mkpKTA1dUVycnJxfYSVWoq4OT0+OcHD3jTTpFhxxMRFVumfH+bPHIzePBgZGVlYdKkSUhLS0Pfvn1RsWJFfPHFF8U+sSGFEnm8AF9B8VEHRESKYPLIzZOSkpKQk5MDDw8Pc8ZUqDhyozCF9cgEdjwRUbFSqCM3T3J3dy/I7kQFVxiPTOCjDoiISrR8TSg2tMhersuXLxcoIKJ8M9cjE/ioAyKiEs3k5Gbs2LE67zMzMxETE4MdO3Zg4sSJ5oqLyHR8ZAIRESEfyc2YMWMMln/99dc4fvx4gQMiIiIiKoh8PVvKkM6dO2Pjxo3mao6IiIgoX8yW3GzYsAFubm7mao6IiIgoX0y+LOXv768zoVhEkJiYiL///huLFi0ya3BEREREpjI5uenRo4fOeysrK5QvXx6tW7dGrVq1zBUXERERUb6YlNxkZWWhSpUq6NixIypUqFBYMZVYil8o11wnaE7FtrOIiMhSTEpu1Go13nvvPcTGxhZWPCVWYS2UW2wo/gSJiEgpTJ5Q3LRpU8TExBRGLCWa4hfKLYwTNKdi1VlERGRJJs+5GT58ON5//31cu3YNAQEBcHxq0bQGDRqYLbiSSvEL5ZrrBM2p2HYWEREVNaMfnPn2229jwYIFKFOmjH4jKhVEBCqVCtnZ2eaO0awK68GZin/YpeJPkIiIijNTvr+NTm6sra2RkJCAhw8fPrOer6+v8ZFaAJObfFL8CRIRUXFWKE8Fz82BinvyQkRERKWbSROKn/U0cCIiIqLiwKQJxX5+fs9NcO7cuVOggIiIiIgKwqTkZvr06XB1dS2sWIiIiIgKzKTkpk+fPvDw8CisWIiIiIgKzOg5N5xvQ0RERCWB0cmNkXeMExEREVmU0ZelcnJyCjMOIiIiIrMw+dlSRERERMUZkxsiIiJSFCY3REREpChMboiIiEhRmNwQERGRojC5ISIiIkVhckNERESKYtLjF6iEEgHS0grWRmqqeWIhIiIqZExulE4EaNECOHTI0pEQEREVCV6WUrq0NPMmNsHBgEZjvvaIiIjMjCM3pcnNm4CjY8Ha0GgAPkSViIiKMSY3pYmjY8GTGyIiomKOl6WIiIhIUZjcEBERkaIwuSEiIiJFYXJDREREisLkhoiIiBSFyQ0REREpCpMbIiIiUhQmN0RERKQoTG6IiIhIUZjcEBERkaIwuSEiIiJFYXJDREREisLkhoiIiBSFyQ0REREpCpMbIiIiUhQmN0RERKQoaksHQHkQAdLSCt5OamrB2yAiIipBmNwURyJAixbAoUOWjoSIiKjE4WWp4igtzfyJTXAwoNGYt00iIqJiiCM3xd3Nm4CjY8Hb0WgAlarg7RARERVzFh+5WbRoEapWrQp7e3sEBAQgKioqz7qbNm1C+/btUb58ebi4uKB58+b49ddfizBaC3B0NM+LiQ0REZUSFk1uIiIiMHbsWHz00UeIiYlBy5Yt0blzZ8THxxus/9tvv6F9+/bYtm0boqOj0aZNG3Tr1g0xMTFFHDkREREVVyoREUsdvGnTpmjcuDEWL16sLatduzZ69OiBmTNnGtVG3bp1ERISgqlTpxpVPyUlBa6urkhOToaLi0u+4jYkNRVwcnr884MHBbySZNbGiIiISj5Tvr8tNnLz6NEjREdHo0OHDjrlHTp0wCEjJ9Pm5OTg/v37cHNzK4wQiYiIqASy2ITipKQkZGdnw9PTU6fc09MTiYmJRrUxd+5cpKamonfv3nnWycjIQEZGhvZ9SkpK/gImIiKiEsHiE4pVT010FRG9MkPWrVuHadOmISIiAh4eHnnWmzlzJlxdXbUvHx+fAsdMRERExZfFkht3d3dYW1vrjdLcunVLbzTnaRERERgyZAh++OEHtGvX7pl1p0yZguTkZO3r6tWrBY6diIiIii+LJTe2trYICAhAZGSkTnlkZCSCgoLy3G/dunUYNGgQ1q5di65duz73OHZ2dnBxcdF5ERERkXJZdBG/8ePH46233kJgYCCaN2+OpUuXIj4+HqGhoQAej7pcv34dq1atAvA4sRkwYAC++OILNGvWTDvq4+DgAFdXV4udBxERERUfFk1uQkJCcPv2bcyYMQMJCQmoV68etm3bBl9fXwBAQkKCzpo333zzDbKysjBixAiMGDFCWz5w4ECsXLmyqMMnIiKiYsii69xYAte5ISIiKnlKxDo3RERERIWByQ0REREpCpMbIiIiUhQmN0RERKQoTG6IiIhIUZjcEBERkaIwuSEiIiJFYXJDREREisLkhoiIiBSFyQ0REREpCpMbIiIiUhQmN0RERKQoTG6IiIhIUZjcEBERkaIwuSEiIiJFYXJDREREisLkhoiIiBSFyQ0REREpCpMbIiIiUhQmN0RERKQoTG6IiIhIUZjcEBERkaIwuSEiIiJFYXJDREREisLkhoiIiBSFyQ0REREpCpMbIiIiUhS1pQNQpNRUy+5PRERUijG5MRcRAKrHP3t6AEizZDRERESlFi9LmUtaISQzwcGARmP+domIiBSMIzeF4XIc4OFY8HY0GkClKng7REREpQiTm8Lg6Pj4RUREREWOl6WIiIhIUZjcEBERkaIwuSEiIiJFYXJDREREisLkhoiIiBSFyQ0REREpCpMbIiIiUhQmN0RERKQoTG6IiIhIUZjcEBERkaIwuSEiIiJFYXJDREREisLkhoiIiBSFyQ0REREpCpMbIiIiUhQmN0RERKQoTG6IiIhIUZjcEBERkaIwuSEiIiJFYXJDREREisLkhoiIiBSFyQ0REREpCpMbIiIiUhQmN0RERKQoTG6IiIhIUSye3CxatAhVq1aFvb09AgICEBUV9cz6+/fvR0BAAOzt7VGtWjUsWbKkiCIlIiKiksCiyU1ERATGjh2Ljz76CDExMWjZsiU6d+6M+Ph4g/Xj4uLQpUsXtGzZEjExMfjwww8xevRobNy4sYgjJyIiouJKJSJiqYM3bdoUjRs3xuLFi7VltWvXRo8ePTBz5ky9+h988AG2bt2K2NhYbVloaCh+//13HD582KhjpqSkwNXVFcnJyXBxcSn4Sfwj9VYqnDwdAQAPbqbC0cPRbG0TERGVdqZ8f1ts5ObRo0eIjo5Ghw4ddMo7dOiAQ4cOGdzn8OHDevU7duyI48ePIzMz0+A+GRkZSElJ0XkRERGRclksuUlKSkJ2djY8PT11yj09PZGYmGhwn8TERIP1s7KykJSUZHCfmTNnwtXVVfvy8fExzwkQERFRsWTxCcUqlUrnvYjolT2vvqHyXFOmTEFycrL2dfXq1QJGbJjGXYMHN1Px4GYqNO6aQjkGERERPZ/aUgd2d3eHtbW13ijNrVu39EZnclWoUMFgfbVajXLlyhncx87ODnZ2duYJ+hlUVirOsyEiIioGLDZyY2tri4CAAERGRuqUR0ZGIigoyOA+zZs316u/c+dOBAYGwsbGptBiJSIiopLDopelxo8fj2XLlmHFihWIjY3FuHHjEB8fj9DQUACPLykNGDBAWz80NBRXrlzB+PHjERsbixUrVmD58uWYMGGCpU6BiIiIihmLXZYCgJCQENy+fRszZsxAQkIC6tWrh23btsHX1xcAkJCQoLPmTdWqVbFt2zaMGzcOX3/9Nby9vbFw4UK8/vrrljoFIiIiKmYsus6NJRTWOjdERERUeErEOjdEREREhYHJDRERESkKkxsiIiJSFCY3REREpChMboiIiEhRmNwQERGRojC5ISIiIkVhckNERESKwuSGiIiIFMWij1+whNwFmVNSUiwcCRERERkr93vbmAcrlLrk5v79+wAAHx8fC0dCREREprp//z5cXV2fWafUPVsqJycHN27cgLOzM1QqlVnbTklJgY+PD65evcrnVhUi9nPRYD8XDfZz0WFfF43C6mcRwf379+Ht7Q0rq2fPqil1IzdWVlaoVKlSoR7DxcWF/3CKAPu5aLCfiwb7ueiwr4tGYfTz80ZscnFCMRERESkKkxsiIiJSFCY3ZmRnZ4ewsDDY2dlZOhRFYz8XDfZz0WA/Fx32ddEoDv1c6iYUExERkbJx5IaIiIgUhckNERERKQqTGyIiIlIUJjdERESkKExuTLRo0SJUrVoV9vb2CAgIQFRU1DPr79+/HwEBAbC3t0e1atWwZMmSIoq0ZDOlnzdt2oT27dujfPnycHFxQfPmzfHrr78WYbQll6mf51wHDx6EWq1Go0aNCjdAhTC1nzMyMvDRRx/B19cXdnZ2qF69OlasWFFE0ZZcpvbzmjVr0LBhQ2g0Gnh5eWHw4MG4fft2EUVbMv3222/o1q0bvL29oVKpsGXLlufuY5HvQSGjrV+/XmxsbOTbb7+Vc+fOyZgxY8TR0VGuXLlisP7ly5dFo9HImDFj5Ny5c/Ltt9+KjY2NbNiwoYgjL1lM7ecxY8bIrFmz5OjRo3LhwgWZMmWK2NjYyIkTJ4o48pLF1H7Ode/ePalWrZp06NBBGjZsWDTBlmD56efu3btL06ZNJTIyUuLi4uT//u//5ODBg0UYdcljaj9HRUWJlZWVfPHFF3L58mWJioqSunXrSo8ePYo48pJl27Zt8tFHH8nGjRsFgGzevPmZ9S31PcjkxgRNmjSR0NBQnbJatWrJ5MmTDdafNGmS1KpVS6fs3XfflWbNmhVajEpgaj8bUqdOHZk+fbq5Q1OU/PZzSEiI/Otf/5KwsDAmN0YwtZ+3b98urq6ucvv27aIITzFM7efPP/9cqlWrplO2cOFCqVSpUqHFqDTGJDeW+h7kZSkjPXr0CNHR0ejQoYNOeYcOHXDo0CGD+xw+fFivfseOHXH8+HFkZmYWWqwlWX76+Wk5OTm4f/8+3NzcCiNERchvP4eHh+PSpUsICwsr7BAVIT/9vHXrVgQGBmL27NmoWLEi/Pz8MGHCBDx8+LAoQi6R8tPPQUFBuHbtGrZt2wYRwc2bN7FhwwZ07dq1KEIuNSz1PVjqHpyZX0lJScjOzoanp6dOuaenJxITEw3uk5iYaLB+VlYWkpKS4OXlVWjxllT56eenzZ07F6mpqejdu3dhhKgI+ennixcvYvLkyYiKioJazT8dxshPP1++fBkHDhyAvb09Nm/ejKSkJAwfPhx37tzhvJs85Kefg4KCsGbNGoSEhCA9PR1ZWVno3r07vvzyy6IIudSw1PcgR25MpFKpdN6LiF7Z8+obKiddpvZzrnXr1mHatGmIiIiAh4dHYYWnGMb2c3Z2Nvr27Yvp06fDz8+vqMJTDFM+zzk5OVCpVFizZg2aNGmCLl26YN68eVi5ciVHb57DlH4+d+4cRo8ejalTpyI6Oho7duxAXFwcQkNDiyLUUsUS34P875eR3N3dYW1trfe/gFu3bullpbkqVKhgsL5arUa5cuUKLdaSLD/9nCsiIgJDhgzBjz/+iHbt2hVmmCWeqf18//59HD9+HDExMRg5ciSAx1/CIgK1Wo2dO3eibdu2RRJ7SZKfz7OXlxcqVqwIV1dXbVnt2rUhIrh27RpeeOGFQo25JMpPP8+cORPBwcGYOHEiAKBBgwZwdHREy5Yt8cknn3Bk3Uws9T3IkRsj2draIiAgAJGRkTrlkZGRCAoKMrhP8+bN9erv3LkTgYGBsLGxKbRYS7L89DPweMRm0KBBWLt2La+ZG8HUfnZxccHp06dx8uRJ7Ss0NBQ1a9bEyZMn0bRp06IKvUTJz+c5ODgYN27cwIMHD7RlFy5cgJWVFSpVqlSo8ZZU+enntLQ0WFnpfgVaW1sD+N/IAhWcxb4HC3W6ssLk3mq4fPlyOXfunIwdO1YcHR3lr7/+EhGRyZMny1tvvaWtn3sL3Lhx4+TcuXOyfPly3gpuBFP7ee3ataJWq+Xrr7+WhIQE7evevXuWOoUSwdR+fhrvljKOqf18//59qVSpkrzxxhty9uxZ2b9/v7zwwgsydOhQS51CiWBqP4eHh4tarZZFixbJpUuX5MCBAxIYGChNmjSx1CmUCPfv35eYmBiJiYkRADJv3jyJiYnR3nJfXL4HmdyY6OuvvxZfX1+xtbWVxo0by/79+7XbBg4cKK1atdKpv2/fPvH39xdbW1upUqWKLF68uIgjLplM6edWrVoJAL3XwIEDiz7wEsbUz/OTmNwYz9R+jo2NlXbt2omDg4NUqlRJxo8fL2lpaUUcdcljaj8vXLhQ6tSpIw4ODuLl5SX9+vWTa9euFXHUJcvevXuf+fe2uHwPqkQ4/kZERETKwTk3REREpChMboiIiEhRmNwQERGRojC5ISIiIkVhckNERESKwuSGiIiIFIXJDRERESkKkxsi0rFy5UqUKVPG0mHkW5UqVbBgwYJn1pk2bRoaNWpUJPEQUdFjckOkQIMGDYJKpdJ7/fnnn5YODStXrtSJycvLC71790ZcXJxZ2j927BiGDRumfa9SqbBlyxadOhMmTMDu3bvNcry8PH2enp6e6NatG86ePWtyOyU52SSyBCY3RArVqVMnJCQk6LyqVq1q6bAAPH4QZ0JCAm7cuIG1a9fi5MmT6N69O7Kzswvcdvny5aHRaJ5Zx8nJqVCfSJzryfP85ZdfkJqaiq5du+LRo0eFfmyi0ozJDZFC2dnZoUKFCjova2trzJs3D/Xr14ejoyN8fHwwfPhwnSdQP+33339HmzZt4OzsDBcXFwQEBOD48ePa7YcOHcJLL70EBwcH+Pj4YPTo0UhNTX1mbCqVChUqVICXlxfatGmDsLAwnDlzRjuytHjxYlSvXh22traoWbMmVq9erbP/tGnTULlyZdjZ2cHb2xujR4/WbnvyslSVKlUAAD179oRKpdK+f/Ky1K+//gp7e3vcu3dP5xijR49Gq1atzHaegYGBGDduHK5cuYLz589r6zzr97Fv3z4MHjwYycnJ2hGgadOmAQAePXqESZMmoWLFinB0dETTpk2xb9++Z8ZDVFowuSEqZaysrLBw4UKcOXMG3333Hfbs2YNJkyblWb9fv36oVKkSjh07hujoaEyePBk2NjYAgNOnT6Njx4547bXXcOrUKURERODAgQMYOXKkSTE5ODgAADIzM7F582aMGTMG77//Ps6cOYN3330XgwcPxt69ewEAGzZswPz58/HNN9/g4sWL2LJlC+rXr2+w3WPHjgEAwsPDkZCQoH3/pHbt2qFMmTLYuHGjtiw7Oxs//PAD+vXrZ7bzvHfvHtauXQsA2v4Dnv37CAoKwoIFC7QjQAkJCZgwYQIAYPDgwTh48CDWr1+PU6dOoVevXujUqRMuXrxodExEilXoj+YkoiI3cOBAsba2FkdHR+3rjTfeMFj3hx9+kHLlymnfh4eHi6urq/a9s7OzrFy50uC+b731lgwbNkynLCoqSqysrOThw4cG93m6/atXr0qzZs2kUqVKkpGRIUFBQfLOO+/o7NOrVy/p0qWLiIjMnTtX/Pz85NGjRwbb9/X1lfnz52vfA5DNmzfr1Hn6ieajR4+Wtm3bat//+uuvYmtrK3fu3CnQeQIQR0dH0Wg02qcnd+/e3WD9XM/7fYiI/Pnnn6JSqeT69es65S+//LJMmTLlme0TlQZqy6ZWRFRY2rRpg8WLF2vfOzo6AgD27t2Lzz77DOfOnUNKSgqysrKQnp6O1NRUbZ0njR8/HkOHDsXq1avRrl079OrVC9WrVwcAREdH488//8SaNWu09UUEOTk5iIuLQ+3atQ3GlpycDCcnJ4gI0tLS0LhxY2zatAm2traIjY3VmRAMAMHBwfjiiy8AAL169cKCBQtQrVo1dOrUCV26dEG3bt2gVuf/z1m/fv3QvHlz3LhxA97e3lizZg26dOmCsmXLFug8nZ2dceLECWRlZWH//v34/PPPsWTJEp06pv4+AODEiRMQEfj5+emUZ2RkFMlcIqLijskNkUI5OjqiRo0aOmVXrlxBly5dEBoain//+99wc3PDgQMHMGTIEGRmZhpsZ9q0aejbty9++eUXbN++HWFhYVi/fj169uyJnJwcvPvuuzpzXnJVrlw5z9hyv/StrKzg6emp9yWuUql03ouItszHxwfnz59HZGQkdu3aheHDh+Pzzz/H/v37dS73mKJJkyaoXr061q9fj/feew+bN29GeHi4dnt+z9PKykr7O6hVqxYSExMREhKC3377DUD+fh+58VhbWyM6OhrW1tY625ycnEw6dyIlYnJDVIocP34cWVlZmDt3LqysHk+5++GHH567n5+fH/z8/DBu3Di8+eabCA8PR8+ePdG4cWOcPXtWL4l6nie/9J9Wu3ZtHDhwAAMGDNCWHTp0SGd0xMHBAd27d0f37t0xYsQI1KpVC6dPn0bjxo312rOxsTHqLqy+fftizZo1qFSpEqysrNC1a1fttvye59PGjRuHefPmYfPmzejZs6dRvw9bW1u9+P39/ZGdnY1bt26hZcuWBYqJSIk4oZioFKlevTqysrLw5Zdf4vLly1i9erXeZZInPXz4ECNHjsS+fftw5coVHDx4EMeOHdMmGh988AEOHz6MESNG4OTJk7h48SK2bt2KUaNG5TvGiRMnYuXKlViyZAkuXryIefPmYdOmTdqJtCtXrsTy5ctx5swZ7Tk4ODjA19fXYHtVqlTB7t27kZiYiLt37+Z53H79+uHEiRP49NNP8cYbb8De3l67zVzn6eLigqFDhyIsLAwiYtTvo0qVKnjw4AF2796NpKQkpKWlwc/PD/369cOAAQOwadMmxMXF4dixY5g1axa2bdtmUkxEimTJCT9EVDgGDhwor776qsFt8+bNEy8vL3FwcJCOHTvKqlWrBIDcvXtXRHQnsGZkZEifPn3Ex8dHbG1txdvbW0aOHKkzifbo0aPSvn17cXJyEkdHR2nQoIF8+umnecZmaILs0xYtWiTVqlUTGxsb8fPzk1WrVmm3bd68WZo2bSouLi7i6OgozZo1k127dmm3Pz2heOvWrVKjRg1Rq9Xi6+srIvoTinO9+OKLAkD27Nmjt81c53nlyhVRq9USEREhIs//fYiIhIaGSrly5QSAhIWFiYjIo0ePZOrUqVKlShWxsbGRChUqSM+ePeXUqVN5xkRUWqhERCybXhERERGZDy9LERERkaIwuSEiIiJFYXJDREREisLkhoiIiBSFyQ0REREpCpMbIiIiUhQmN0RERKQoTG6IiIhIUZjcEBERkaIwuSEiIiJFYXJDREREisLkhoiIiBTl/wGOAay7sFUO1QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fprP, tprP, tP= metrics.roc_curve(predP[\"Label\"],predP[0],drop_intermediate=False)\n",
    "fprR, tprR, tR = metrics.roc_curve(predR[\"Label\"],predR[0],drop_intermediate=False)\n",
    "plt.plot(fprP,tprP,label=\"gc-ims-tools (Python), AUC=\"+str(round(AUCP, 2)),color='r')\n",
    "plt.plot(fprR,tprR,label=\"GCIMS (R), AUC=\"+str(round(AUCR, 2)),color='b')\n",
    "plt.xlabel(\"False Positive Rate\") \n",
    "plt.ylabel(\"True Positive Rate\")\n",
    "plt.title(\"ROC Comparison\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc4ed766",
   "metadata": {},
   "source": [
    "## PERMUTATION TEST"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "44ce22a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to generae each model with a random test label\n",
    "def permutation(data,lvh,i):\n",
    "    Xtrain=data.loc[annotations[\"ID\"]!=i].values.tolist()\n",
    "    Xtest=data.loc[annotations[\"ID\"]==i].values.tolist()\n",
    "    Ytrain=annotations[annotations[\"ID\"]!=i][\"Sex_b\"].values.tolist()\n",
    "    groups = np.unique(Ytrain)\n",
    "    Ytest=rd.sample(groups.tolist(),1)*len(Xtest)\n",
    "    LV=lvh[np.unique(annotations[\"ID\"]).tolist().index(i)]\n",
    "    pls=PLSRegression(n_components=LV, scale=True)\n",
    "    pls.fit(Xtrain,Ytrain)\n",
    "    predPer=pd.DataFrame(pls.predict(Xtest))\n",
    "    predPer[\"Label\"]=Ytest\n",
    "    return predPer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "0fcceecb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAHFCAYAAAAHcXhbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABEwElEQVR4nO3deVxUZf//8fewLwqWC+AGqLjmlguKqZRKi7a5ZLnc2lcrbzWzLNO8VexOSS3TMi3LrXJr0erWu9Lb1FzLtUxNKzdMCU0BFWU9vz/4MTWBAgPDzMHX8/E4j8eZ6yzzmSuEd9c5Zy6LYRiGAAAATMrN2QUAAAAUB2EGAACYGmEGAACYGmEGAACYGmEGAACYGmEGAACYGmEGAACYGmEGAACYGmEGAACYGmEGMKFFixbJYrFo165d+W7v1q2bwsLCbNrCwsI0cODAIr3Ptm3bFBsbq6SkJPsKvQGtWLFCjRo1kq+vrywWi/bt25fvfhs3bpTFYtHHH3+c7/bhw4fLYrHkad+zZ486d+6scuXKqUKFCurevbuOHj1akh8BMB3CDHCDWLVqlcaPH1+kY7Zt26ZJkyYRZgrp7Nmz6t+/v2rXrq0vv/xS27dvV926dUvs/D/99JOio6OVnp6uDz/8UAsWLNCRI0fUvn17nT17tsTeBzAbD2cXAKB0NG/e3NklFFlGRoYsFos8PMzxq+rIkSPKyMhQv3791LFjxxI//4QJE+Tt7a3Vq1crICBAktSiRQtFRETolVde0dSpU0v8PQEzYGQGuEH8/TJTdna2XnrpJdWrV0++vr6qUKGCmjRpolmzZkmSYmNj9dxzz0mSwsPDZbFYZLFYtHHjRuvx06ZNU/369eXt7a0qVaroH//4h06dOmXzvoZhaMqUKQoNDZWPj49atmypdevWKTo6WtHR0db9ci+7vP/++xo1apSqVasmb29v/fLLLzp79qyGDh2qhg0bqly5cqpSpYruuOMObd682ea9jh8/LovFounTp2vq1KkKCwuTr6+voqOjrUFjzJgxqlq1qgIDA/Xggw8qMTGxUP33+eefq23btvLz81P58uXVpUsXbd++3bp94MCBuu222yRJvXv3lsVisfl8xZWZmanVq1erR48e1iAjSaGhobr99tu1atWqEnsvwGzM8b87APKVlZWlzMzMPO2GYRR47LRp0xQbG6t//etf6tChgzIyMvTTTz9ZLykNHjxY58+f1xtvvKGVK1cqJCREktSwYUNJ0j//+U/NmzdPw4cPV7du3XT8+HGNHz9eGzdu1J49e1SpUiVJ0rhx4xQXF6fHH39c3bt3V3x8vAYPHqyMjIx8L8GMHTtWbdu21VtvvSU3NzdVqVLFegll4sSJCg4O1qVLl7Rq1SpFR0dr/fr1eULDm2++qSZNmujNN99UUlKSRo0apXvvvVeRkZHy9PTUggULdOLECT377LMaPHiwPv/88+v21dKlS9W3b1/FxMRo2bJlSktL07Rp06zvf9ttt2n8+PFq3bq1hg0bpilTpuj222+3CR3Xkp2dXaj/hr/++quuXLmiJk2a5Nm3SZMmWrduna5evSofH58C3xMocwwAprNw4UJD0nWX0NBQm2NCQ0ONAQMGWF9369bNaNas2XXfZ/r06YYk49ixYzbthw4dMiQZQ4cOtWn/9ttvDUnGCy+8YBiGYZw/f97w9vY2evfubbPf9u3bDUlGx44drW0bNmwwJBkdOnQo8PNnZmYaGRkZRqdOnYwHH3zQ2n7s2DFDktG0aVMjKyvL2j5z5kxDknHffffZnGfkyJGGJCM5Ofma75WVlWVUrVrVaNy4sc05L168aFSpUsWIiorK8xk++uijAj9D7r4FLbm2bt1qSDKWLVuW51xTpkwxJBmnT58u8H2BsojLTICJvffee9q5c2eeJfdyx/W0bt1a33//vYYOHaqvvvpKKSkphX7fDRs2SFKep6Nat26tBg0aaP369ZKkHTt2KC0tTQ899JDNfm3atMnztFWuHj165Nv+1ltv6dZbb5WPj488PDzk6emp9evX69ChQ3n2veeee+Tm9uevtwYNGkiSunbtarNfbvvJkyev8Umlw4cP6/Tp0+rfv7/NOcuVK6cePXpox44dSk1NvebxBZk6dWq+/w3/3me58nvCqTDbgLKMy0yAiTVo0EAtW7bM0x4YGKj4+PjrHjt27Fj5+/vrgw8+0FtvvSV3d3d16NBBU6dOzfecf/XHH39IkvXS019VrVpVJ06csNkvKCgoz375tV3rnDNmzNCoUaM0ZMgQ/fvf/1alSpXk7u6u8ePH5xtmbr75ZpvXXl5e122/evVqvrX89TNc67NmZ2frwoUL8vPzu+Y5rqdWrVr59nflypVtXlesWNGmnr86f/68LBaLKlSoYFcNgNkxMgPcoDw8PPTMM89oz549On/+vJYtW6b4+HjdeeedBY405P5hPXPmTJ5tp0+ftt4vk7vf77//nme/hISEfM+d3+jCBx98oOjoaM2dO1ddu3ZVZGSkWrZsqYsXL17/Q5aAgj6rm5ubbrrpJofXUbt2bfn6+mr//v15tu3fv1916tThfhncsAgzAFShQgX17NlTw4YN0/nz53X8+HFJkre3tyTpypUrNvvfcccdknJCxl/t3LlThw4dUqdOnSRJkZGR8vb21ooVK2z227Fjh3X0pjAsFou1llw//PCDzdNEjlKvXj1Vq1ZNS5cutbkp9/Lly/rkk0+sTzg5moeHh+69916tXLnSJsSdPHlSGzZsUPfu3R1eA+CquMwE3KDuvfde3XLLLWrZsqUqV66sEydOaObMmQoNDVVERIQkqXHjxpKkWbNmacCAAfL09FS9evVUr149Pf7443rjjTfk5uamu+++2/o0U40aNfT0009Lyrms88wzzyguLk433XSTHnzwQZ06dUqTJk1SSEiIzT0o19OtWzf9+9//1sSJE9WxY0cdPnxYL774osLDw/N9Eqgkubm5adq0aerbt6+6deumJ554QmlpaZo+fbqSkpL08ssvO/T9/2rSpElq1aqVunXrpjFjxujq1auaMGGCKlWqpFGjRpVaHYCrIcwAN6jbb79dn3zyid59912lpKQoODhYXbp00fjx4+Xp6SlJio6O1tixY7V48WK98847ys7O1oYNG6yXfGrXrq358+frzTffVGBgoO666y7FxcVZL81I0uTJk+Xv76+33npLCxcuVP369TV37lyNGzeu0Pd4jBs3TqmpqZo/f76mTZumhg0b6q233tKqVaus33vjSH369JG/v7/i4uLUu3dvubu7q02bNtqwYYOioqIc/v656tevr40bN+r5559Xz5495eHhoTvuuEOvvPJKnntsgBuJxTAK8YUUAFCCjh07pvr162vixIl64YUXnF0OAJMjzABwqO+//17Lli1TVFSUAgICdPjwYU2bNk0pKSn68ccfr/lUEwAUFpeZADiUv7+/du3apfnz5yspKUmBgYGKjo7W5MmTCTIASgQjMwAAwNR4NBsAAJgaYQYAAJgaYQYAAJhamb8BODs7W6dPn1b58uWZhA0AAJMwDEMXL15U1apVC/yCzTIfZk6fPq0aNWo4uwwAAGCH+Ph4Va9e/br7lPkwU758eUk5nREQEODkagAAQGGkpKSoRo0a1r/j11Pmw0zupaWAgADCDAAAJlOYW0S4ARgAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJiah7MLAAAzChuzxu5jj7/ctQQrAcDIDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDWnhpnMzEz961//Unh4uHx9fVWrVi29+OKLys7Otu5jGIZiY2NVtWpV+fr6Kjo6WgcOHHBi1QAAwJU4NcxMnTpVb731lmbPnq1Dhw5p2rRpmj59ut544w3rPtOmTdOMGTM0e/Zs7dy5U8HBwerSpYsuXrzoxMoBAICrcGqY2b59u+6//3517dpVYWFh6tmzp2JiYrRr1y5JOaMyM2fO1Lhx49S9e3fdcsstWrx4sVJTU7V06VJnlg4AAFyEU8PMbbfdpvXr1+vIkSOSpO+//15btmzRPffcI0k6duyYEhISFBMTYz3G29tbHTt21LZt2/I9Z1pamlJSUmwWAABQdnk4882ff/55JScnq379+nJ3d1dWVpYmT56sRx55RJKUkJAgSQoKCrI5LigoSCdOnMj3nHFxcZo0aZJjCwcAAC7DqSMzK1as0AcffKClS5dqz549Wrx4sV555RUtXrzYZj+LxWLz2jCMPG25xo4dq+TkZOsSHx/vsPoBAIDzOXVk5rnnntOYMWP08MMPS5IaN26sEydOKC4uTgMGDFBwcLCknBGakJAQ63GJiYl5RmtyeXt7y9vb2/HFAwAAl+DUkZnU1FS5udmW4O7ubn00Ozw8XMHBwVq3bp11e3p6ujZt2qSoqKhSrRUAALgmp47M3HvvvZo8ebJq1qypRo0aae/evZoxY4b+7//+T1LO5aWRI0dqypQpioiIUEREhKZMmSI/Pz/16dPHmaUDAAAX4dQw88Ybb2j8+PEaOnSoEhMTVbVqVT3xxBOaMGGCdZ/Ro0frypUrGjp0qC5cuKDIyEitXbtW5cuXd2LlAADAVVgMwzCcXYQjpaSkKDAwUMnJyQoICHB2OQDKiLAxa+w+9vjLXUuwEqBsKsrfb+ZmAgDc0LZt2yZ3d3fdddddebZt3LhRFotFSUlJebY1a9ZMsbGxNm179+5Vr169FBQUJB8fH9WtW1ePPfaY9fvUSsKFCxfUv39/BQYGKjAwUP3798+3vr+yWCz5LtOnT8+zr2EYuvvuu2WxWPTpp5+WWN2O5NTLTADgTMUZXUHZsWDBAj355JN69913dfLkSdWsWdOu86xevVo9evTQnXfeqSVLlqh27dpKTEzURx99pPHjx2vFihUlUm+fPn106tQpffnll5Kkxx9/XP3799d//vOfax5z5swZm9dffPGFBg0apB49euTZd+bMmdf8+hNXRZgBANywLl++rA8//FA7d+5UQkKCFi1aZHPfZmGlpqbq0Ucf1T333KNVq1ZZ28PDwxUZGVngyElhHTp0SF9++aV27NihyMhISdI777yjtm3b6vDhw6pXr16+x+V+1Umuzz77TLfffrtq1apl0/79999rxowZ2rlzp81Xorg6LjMBAG5YK1asUL169VSvXj3169dPCxculD23kn711Vc6d+6cRo8ene/2ChUqWNfLlSt33eXuu+++5vts375dgYGB1iAjSW3atFFgYOA1p/n5u99//11r1qzRoEGDbNpTU1P1yCOPaPbs2XnCj6tjZAYAcMOaP3+++vXrJ0m66667dOnSJa1fv16dO3cu0nl+/vlnSVL9+vUL3Hffvn3X3e7r63vNbQkJCapSpUqe9ipVqlinACrI4sWLVb58eXXv3t2m/emnn1ZUVJTuv//+Qp3HlRBmAAA3pMOHD+u7777TypUrJUkeHh7q3bu3FixYUOQwU5TRnDp16hTp3H+X3/0s15vm5+8WLFigvn37ysfHx9r2+eef6+uvv9bevXuLVZuzEGYAADek+fPnKzMzU9WqVbO2GYYhT09PXbhwQTfddJP1keDk5GSbS0WSlJSUpMDAQElS3bp1JUk//fST2rZte933LVeu3HW3t2/fXl988UW+24KDg/X777/naT979uw1p/n5q82bN+vw4cN5bkb++uuv9euvv+b5jD169FD79u21cePGAs/tTIQZAMANJzMzU++9955effVVxcTE2Gzr0aOHlixZouHDhysiIkJubm7auXOnQkNDrfucOXNGv/32m/WG25iYGFWqVEnTpk2zuQE4V1JSkjUoFOcyU9u2bZWcnKzvvvtOrVu3liR9++23Sk5OLtQ0P/Pnz1eLFi3UtGlTm/YxY8Zo8ODBNm2NGzfWa6+9pnvvvbfA8zobYQYAcMNZvXq1Lly4oEGDBllHV3L17NlT8+fP1/Dhw1W+fHk98cQTGjVqlDw8PNS0aVOdPn1a48aNU4MGDaxByN/fX++++6569eql++67TyNGjFCdOnV07tw5ffjhhzp58qSWL18uqXiXmRo0aKC77rpLjz32mN5++21JOY9md+vWzeZJpvr16ysuLk4PPvigtS0lJUUfffSRXn311TznDQ4Ozvem35o1ayo8PNzueksLTzMBAG448+fPV+fOnfMEGSlnZGbfvn3as2ePJOm1117T4MGD9cILL6hRo0bq27evwsPDtXbtWnl4/DkmcP/992vbtm3y9PRUnz59VL9+fT3yyCNKTk7WSy+9VGK1L1myRI0bN1ZMTIxiYmLUpEkTvf/++zb7HD58WMnJyTZty5cvl2EYeuSRR0qsFlfBdAYAbljO+tI8pjMACsZ0BgAA4IZBmAEAAKZGmAEAAKZGmAEAAKZGmAEAAKZGmAEAAKZGmAEAAKbGNwADsCrO967w3SkAnIWRGQAAYGqEGQAAYGqEGQAAYGqEGQAAYGqEGQAAYGqEGQAAYGqEGQAAYGqEGQAASlhW1hXt3t1Ku3e3UlbWFWeXU+bxpXkAAJS4bF28uMu6DsdiZAYAAJgaYQYAAJgal5kAlIjizOtUHMwJBYCRGQAAYGqEGQCA023cuFEWi0VJSUmSpEWLFqlChQpOrQnmQZgBgDIuPj5egwYNUtWqVeXl5aXQ0FA99dRT+uOPP5xST3R0tEaOHGnTFhUVpTNnzigwMLBE3+uTTz5Rw4YN5e3trYYNG2rVqlUFHvPhhx+qWbNm8vPzU2hoqKZPn26zPTd4/X356aefbPbz9KwkT89KJfp5kD/CDACUYUePHlXLli115MgRLVu2TL/88oveeustrV+/Xm3bttX58+edXaIkycvLS8HBwbJYLCV2zu3bt6t3797q37+/vv/+e/Xv318PPfSQvv3222se88UXX6hv374aMmSIfvzxR82ZM0czZszQ7Nmz8+x7+PBhnTlzxrpERERYt7m7+6tdu7Nq1+6s3N39S+wzIX+EGQAow4YNGyYvLy+tXbtWHTt2VM2aNXX33Xfrf//7n3777TeNGzfOuq/FYtGnn35qc3yFChW0aNEi6+vnn39edevWlZ+fn2rVqqXx48crIyPDuj02NlbNmjXT+++/r7CwMAUGBurhhx/WxYsXJUkDBw7Upk2bNGvWLOuIxvHjx/NcZsrPf/7zH7Vo0UI+Pj6qVauWJk2apMzMzGvuP3PmTHXp0kVjx45V/fr1NXbsWHXq1EkzZ8685jHvv/++HnjgAQ0ZMkS1atVS165d9fzzz2vq1KkyDMNm3ypVqig4ONi6uLu7X/O8cCzCDACUUefPn9dXX32loUOHytfX12ZbcHCw+vbtqxUrVuT5I3095cuX16JFi3Tw4EHNmjVL77zzjl577TWbfX799Vd9+umnWr16tVavXq1Nmzbp5ZdfliTNmjVLbdu21WOPPWYd0ahRo0aB7/vVV1+pX79+GjFihA4ePKi3335bixYt0uTJk695zPbt2xUTE2PTduedd2rbtm3XPCYtLU0+Pj42bb6+vjp16pROnDhh0968eXOFhISoU6dO2rBhQ4GfAY5DmAGAMurnn3+WYRhq0KBBvtsbNGigCxcu6OzZs4U+57/+9S9FRUUpLCxM9957r0aNGqUPP/zQZp/s7GwtWrRIt9xyi9q3b6/+/ftr/fr1kqTAwEB5eXnJz8+vSCMakydP1pgxYzRgwADVqlVLXbp00b///W+9/fbb1zwmISFBQUFBNm1BQUFKSEi45jF33nmnVq5cqfXr1ys7O1tHjhyxjuScOXNGkhQSEqJ58+bpk08+0cqVK1WvXj116tRJ33zzjfU8WVlXtHdvtPbujWY6g1LA98wAwA0qd0TGy8ur0Md8/PHHmjlzpn755RddunRJmZmZCggIsNknLCxM5cuXt74OCQlRYmJisWrdvXu3du7caTMSk5WVpatXryo1NVV+fn75Hvf3e3AMw7jufTmPPfaYfv31V3Xr1k0ZGRkKCAjQU089pdjYWGvoqlevnurVq2c9pm3btoqPj9crr7yiDh06/P/WbCUnb7Kuw7EYmQGAMqpOnTqyWCw6ePBgvtt/+uknVa5c2foItMViyXPJ6a/3w+zYsUMPP/yw7r77bq1evVp79+7VuHHjlJ6ebnOMp6enzWuLxaLs7OL9Qc/OztakSZO0b98+67J//379/PPPeS4L5QoODs4zCpOYmJhntObvtU6dOlWXLl3SiRMnlJCQoNatW0vKCWnX0qZNG/38889F/2AoEYzMAEAZVbFiRXXp0kVz5szR008/bXPfTEJCgpYsWaJhw4ZZ2ypXrmy9lCLlXKZKTU21vt66datCQ0Ntbhr++30kheHl5aWsrKwiHXPrrbfq8OHDqlOnTqGPadu2rdatW6enn37a2rZ27VpFRUUVeKy7u7uqVasmSVq2bJnatm2rKlWqXHP/vXv3KiQkpNC1oWQRZgCgDJs9e7aioqJ055136qWXXlJ4eLgOHDig5557TnXr1tWECROs+95xxx2aPXu22rRpo+zsbD3//PM2oyx16tTRyZMntXz5crVq1Upr1qwp1Pe2/F1YWJi+/fZbHT9+XOXKldPNN99c4DETJkxQt27dVKNGDfXq1Utubm764YcftH//fr300kv5HvPUU0+pQ4cOmjp1qu6//3599tln+t///qctW7bY9M+qVaus9/ScO3dOH3/8saKjo3X16lUtXLhQH330kTZt2mQ9ZubMmQoLC1OjRo2Unp6uDz74QJ988ok++eSTIvcFSgaXmQCgDIuIiNDOnTtVq1YtPfTQQwoNDdXdd9+tunXrauvWrSpXrpx131dffVU1atRQhw4d1KdPHz377LM296Lcf//9evrppzV8+HA1a9ZM27Zt0/jx44tc07PPPit3d3c1bNhQlStX1smTJws85s4779Tq1au1bt06tWrVSm3atNGMGTMUGhp6zWOioqK0fPlyLVy4UE2aNNGiRYu0YsUKRUZGWvc5d+6cfv31V5vjFi9erJYtW6pdu3Y6cOCANm7caL3UJEnp6el69tln1aRJE7Vv315btmzRmjVr1L179yL3BUqGxSjKM3kmlJKSosDAQCUnJ+e5SQ2ALWdNFlkcxZlo8kadHHPixImaMWOG1q5dq7Zt2zq1lrIqK+uyNm/OCYrt21/ii/PsUJS/31xmAoAbzKRJk6yXeiIjI+XmxiC9I7i55f+EFUoeYQYAbkCPPvqos0so09zd/dWhw2Vnl3HDII4DAABTI8wAAABTI8wAAFDCsrKu6ocfuuqHH7oqK+uqs8sp87hnBgCAEpel8+f/a12HYzEyAwAATI0wAwAATI0wAwAATI0wAwAATI0wAwAATI0wAwAATI1HswEAKGHu7v6Kji7T8zi7FEZmAACAqRFmAACAqRFmAAAoYVlZV3XgQC8dONCL6QxKAWEGAIASl6WzZz/W2bMfi+kMHI8wAwAATI0wAwBwfZcvSxZLznL5srOrgYshzAAAAFMjzAAAAFMjzAAAAFMjzAAAAFNjOgMAKGVhY9bYfezxl7uWYCUm4ucnXbr057qLc3PzU/v2l6zrcCzCDADA9Vkskr+/s6soNIvFInd389RrdlxmAgAApkaYAQC4vrQ0aeDAnCUtzdnVFCg7O02HDg3UoUMDlZ3t+vWandPDzG+//aZ+/fqpYsWK8vPzU7NmzbR7927rdsMwFBsbq6pVq8rX11fR0dE6cOCAEysGAJS6zExp8eKcJTPT2dUUyDAy9fvvi/X774tlGK5fr9k5NcxcuHBB7dq1k6enp7744gsdPHhQr776qipUqGDdZ9q0aZoxY4Zmz56tnTt3Kjg4WF26dNHFixedVzgAAHAZTr0BeOrUqapRo4YWLlxobQsLC7OuG4ahmTNnaty4cerevbskafHixQoKCtLSpUv1xBNPlHbJAADAxTh1ZObzzz9Xy5Yt1atXL1WpUkXNmzfXO++8Y91+7NgxJSQkKCYmxtrm7e2tjh07atu2bfmeMy0tTSkpKTYLAAAou5waZo4ePaq5c+cqIiJCX331lYYMGaIRI0bovffekyQlJCRIkoKCgmyOCwoKsm77u7i4OAUGBlqXGjVqOPZDAAAAp3JqmMnOztatt96qKVOmqHnz5nriiSf02GOPae7cuTb7WSwWm9eGYeRpyzV27FglJydbl/j4eIfVDwAAnM+pYSYkJEQNGza0aWvQoIFOnjwpSQoODpakPKMwiYmJeUZrcnl7eysgIMBmAQAAZZdTw0y7du10+PBhm7YjR44oNDRUkhQeHq7g4GCtW7fOuj09PV2bNm1SVFRUqdYKAHAiPz8pMTFnMcl0BlFRiYqKSmQ6g1Lg1KeZnn76aUVFRWnKlCl66KGH9N1332nevHmaN2+epJzLSyNHjtSUKVMUERGhiIgITZkyRX5+furTp48zSwcAlCaLRapc2dlVFJrFYpGXl3nqNTunhplWrVpp1apVGjt2rF588UWFh4dr5syZ6tu3r3Wf0aNH68qVKxo6dKguXLigyMhIrV27VuXLl3di5QAAwFU4faLJbt26qVu3btfcbrFYFBsbq9jY2NIrCgDgWtLSpGeeyVmfMUPy9nZuPQXIzk7TL7/k1Funzgy5ubl2vWbn9OkMAAAoUGamNGdOzmKS6QxOn56j06fnMJ1BKSDMAAAAUyPMAAAAUyPMAAAAUyPMAAAAUyPMAAAAUyPMAAAAU3P698wAAFAgX1/p2LE/112cm5uvIiOPWdfhWIQZAIDrc3OTwsKcXUWhWSxu8vUNc3YZNwwuMwEAAFNjZAaAqYWNWePsElAa0tOlceNy1idPlry8nFtPAbKz03XsWE694eGT5ebm2vWaHSMzAADXl5EhvfJKzpKR4exqCmQYGYqPf0Xx8a/IMFy/XrMjzAAAAFMjzAAAAFMjzAAAAFMjzAAAAFMjzAAAAFMjzAAAAFPje2YAAK7P11f68cc/112cm5uvWrX60boOxyLMAABcn5ub1KiRs6soNIvFTf7+5qnX7LjMBAAATI2RGQCA60tPl6ZMyVl/4QVTTGdw4kROvaGhLzCdgYPZNTITFhamF198USdPnizpegAAyCsjQ5o0KWcxyXQGJ05M0okTk5jOoBTYFWZGjRqlzz77TLVq1VKXLl20fPlypaWllXRtAAAABbIrzDz55JPavXu3du/erYYNG2rEiBEKCQnR8OHDtWfPnpKuEQAA4JqKdQNw06ZNNWvWLP3222+aOHGi3n33XbVq1UpNmzbVggULZBhGSdUJAACQr2LdAJyRkaFVq1Zp4cKFWrdundq0aaNBgwbp9OnTGjdunP73v/9p6dKlJVUrAABAHnaFmT179mjhwoVatmyZ3N3d1b9/f7322muqX7++dZ+YmBh16NChxAoFAADIj11hplWrVurSpYvmzp2rBx54QJ6ennn2adiwoR5++OFiFwgAAHA9doWZo0ePKjQ09Lr7+Pv7a+HChXYVBQCADR8f6bvv/lx3cW5uPrr11u+s63Asu8JMYmKiEhISFBkZadP+7bffyt3dXS1btiyR4gAAkCS5u0utWjm7ikKzWNwVEGCees3OrqeZhg0bpvj4+Dztv/32m4YNG1bsogAAAArLrpGZgwcP6tZbb83T3rx5cx08eLDYRQEAYCM9XZo1K2f9qadMMZ3BqVM59Vav/hTTGTiYXSMz3t7e+v333/O0nzlzRh4eTPcEAChhGRnS6NE5i0mmMzh6dLSOHh3NdAalwK4w06VLF40dO1bJycnWtqSkJL3wwgvq0qVLiRUHAABQELuGUV599VV16NBBoaGhat68uSRp3759CgoK0vvvv1+iBQIAAFyPXWGmWrVq+uGHH7RkyRJ9//338vX11aOPPqpHHnkk3++cAQAAcBS7b3Dx9/fX448/XpK1AAAAFJndYebIkSPauHGjEhMTlZ2dbbNtwoQJxS4MAJBX2Jg1dh97/OWuJVgJ4DrsCjPvvPOO/vnPf6pSpUoKDg6WxWKxbrNYLIQZAABQauwKMy+99JImT56s559/vqTrAQAgLx8facOGP9ddnJubj5o23WBdh2PZFWYuXLigXr16lXQtAADkz91dio52dhWFZrG466abop1dxg3Dru+Z6dWrl9auXVvStQAAABSZXSMzderU0fjx47Vjxw41btw4z+PYI0aMKJHiAACQlPOtv/Pm5aw//rjk4l8Dkp2doTNncuoNCXlcbm6uXa/ZWQzDMIp6UHh4+LVPaLHo6NGjxSqqJKWkpCgwMFDJyckKCAhwdjmASyvOkzJwfaZ+munyZalcuZz1S5ckf3/n1lOArKzL2rw5p9727S/J3d2163VFRfn7bdfIzLFjx+wqDAAAoKTZdc9MrvT0dB0+fFiZmZklVQ8AAECR2BVmUlNTNWjQIPn5+alRo0Y6efKkpJx7ZV5++eUSLRAAAOB67AozY8eO1ffff6+NGzfK5y/P+3fu3FkrVqwoseIAAAAKYtc9M59++qlWrFihNm3a2Hz7b8OGDfXrr7+WWHEAAAAFsSvMnD17VlWqVMnTfvnyZZtwA8A+zL8DR+DnCmWVXZeZWrVqpTVr/vxHkRtg3nnnHbVt27ZkKgMAIJe3t7R6dc7i7e3sagpksXircePVatx4tSwW16/X7OwamYmLi9Ndd92lgwcPKjMzU7NmzdKBAwe0fft2bdq0qaRrBADc6Dw8pK7mGR1yc/NQxYrmqdfs7BqZiYqK0tatW5WamqratWtr7dq1CgoK0vbt29WiRYuSrhEAAOCa7BqZkaTGjRtr8eLFJVkLAAD5y8iQlizJWe/b1xTTGfz+e069QUF9mc7AwewKM7nfK3MtNWvWtKsYoCxhagCUJU6/eTg9XXr00Zz1Xr1cPswYRroOH86pt0qVXpJcu16zsyvMhIWFXfeppaysLLsLAgAAKAq7wszevXttXmdkZGjv3r2aMWOGJk+eXCKFAQAAFIZdYaZp06Z52lq2bKmqVatq+vTp6t69e7ELAwAAKAy7bwDOT926dbVz586SPCXgVGa878WMNQNAcdgVZlJSUmxeG4ahM2fOKDY2VhERESVSGAAAQGHYFWYqVKiQ5wZgwzBUo0YNLV++vEQKAwAAKAy7wszXX39tE2bc3NxUuXJl1alTRx4eJXrlCgCAnCkMPvzwz3UXZ7F4q2HDD63rcCy7kkd0dHQJlwEAwHV4eOR8v4xJuLl5/P/vl0FpsGs6g7i4OC1YsCBP+4IFCzR16tRiFwUAAFBYdoWZt99+W/Xr18/T3qhRI7311lvFLgoAABuZmdJHH+UsmZnOrqZA2dmZSkz8SImJHyk72/XrNTu7LjMlJCQoJCQkT3vlypV15syZYhcFAICNtDTpoYdy1i9dyrns5MIMI00HD+bU2779JZXwN6Hgb+wamalRo4a2bt2ap33r1q2qWrVqsYsCAAAoLLui4uDBgzVy5EhlZGTojjvukCStX79eo0eP1qhRo0q0QAAAgOuxK8yMHj1a58+f19ChQ5Weni5J8vHx0fPPP6+xY8eWaIEAAADXY1eYsVgsmjp1qsaPH69Dhw7J19dXERER8jbBs/8AAKBsseuemVwJCQk6f/68ateuLW9vbxmGYfe54uLiZLFYNHLkSGubYRiKjY1V1apV5evrq+joaB04cKA4JQMAgDLGrjDzxx9/qFOnTqpbt67uuece6xNMgwcPtuuemZ07d2revHlq0qSJTfu0adM0Y8YMzZ49Wzt37lRwcLC6dOmiixcv2lM2AAAog+wKM08//bQ8PT118uRJ+fn5Wdt79+6tL7/8skjnunTpkvr27at33nlHN910k7XdMAzNnDlT48aNU/fu3XXLLbdo8eLFSk1N1dKlS+0pGwBgVl5e0sKFOYuXl7OrKZDF4qV69RaqXr2Fslhcv16zsyvMrF27VlOnTlX16tVt2iMiInTixIkinWvYsGHq2rWrOnfubNN+7NgxJSQkKCYmxtrm7e2tjh07atu2bdc8X1pamlJSUmwWAIDJeXpKAwfmLJ6ezq6mQG5ungoJGaiQkIFyc3P9es3OrhuAL1++bDMik+vcuXNFugl4+fLl2r17t3bt2pVnW0JCgiQpKCjIpj0oKOi6gSkuLk6TJk0qdA0AAMDc7BqZ6dChg9577z3ra4vFouzsbE2fPl233357oc4RHx+vp556SkuWLJGPj8819/vr7NxSzuWnv7f91dixY5WcnGxd4uPjC1UPAMCFZWZKa9bkLCaZzuCPP9bojz/WMJ1BKbBrZGb69OmKjo7Wrl27lJ6ertGjR+vAgQM6f/58vt8MnJ/du3crMTFRLVq0sLZlZWXpm2++0ezZs3X48GFJeadOSExMzDNa81fe3t48Ig4AZU1amtStW866SaYz2L8/p16mM3A8u0ZmGjZsqB9++EGtW7dWly5ddPnyZXXv3l179+5V7dq1C3WOTp06af/+/dq3b591admypfr27at9+/apVq1aCg4O1rp166zHpKena9OmTYqKirKnbAAAUAYVOSpmZGQoJiZGb7/9drHuTSlfvrxuueUWmzZ/f39VrFjR2j5y5EhNmTJFERERioiI0JQpU+Tn56c+ffrY/b4AAKBsKXKY8fT01I8//njd+1ZKyujRo3XlyhUNHTpUFy5cUGRkpNauXavy5cs7/L0BAIA52HUR7x//+Ifmz5+vl19+uUSL2bhxo81ri8Wi2NhYxcbGluj7AACAssOuMJOenq53331X69atU8uWLeXv72+zfcaMGSVSHAAAQEGKFGaOHj2qsLAw/fjjj7r11lslSUeOHLHZpzQuPwEAAOQqUpiJiIjQmTNntGHDBkk50xe8/vrr131UGgCAYvPykmbP/nPdxVksXoqImG1dh2MVKcz8fVbsL774QpcvXy7RggAAyMPTUxo2zNlVFJqbm6eqVTNPvWZn1/fM5Pp7uAEAAChtRRqZsVgsee6J4R4ZAIDDZWVJmzfnrLdvL7m7O7eeAhhGlpKScuqtUKG9LBbXrtfsinyZaeDAgdbpAq5evaohQ4bkeZpp5cqVJVchAABXr0q5c/9duiT97e+Oq8nOvqrvv8+pt337S3J3d+16za5IYWbAgAE2r/v161eixQAAABRVkcLMwoULHVUHAACAXYp1AzAAAICzEWYAAICpEWYAAICpEWYAAICp2TXRJAAApcrTU5o27c91F2exeKpWrWnWdTgWYQYA4Pq8vKTnnnN2FYXm5ualmjXNU6/ZcZkJAACYGiMzAADXl5Ul7dmTs37rraaYzuDixZx6y5e/lekMHIwwAwBwfVevSq1b56ybZDqDPXty6mU6A8fjMhMAADA1wgwAADA1wgwAADA1wgwAADA1wgwAADA1wgwAADA1Hs0GALg+T09p4sQ/112cxeKp0NCJ1nU4FmEGAOD6vLyk2FhnV1Fobm5eCg+PdXYZNwwuMwEAAFNjZAYA4Pqys6VDh3LWGzSQ3Fz7/8UNI1upqTn1+vk1kMXi2vWaHWEGAOD6rlyRbrklZ90U0xlc0c6dOfUynYHjERUBAICpMTKDMi9szBpnlwAAcCBGZgAAgKkRZgAAgKkRZgAAgKkRZgAAgKlxAzAAwPV5ekrPPvvnuouzWDxVo8az1nU4FmEGAOD6vLyk6dOdXUWhubl5qXZt89RrdlxmAgAApsbIDADA9WVnSydP5qzXrGmK6QyuXs2p18enJtMZOBhhBgDg+q5ckcLDc9ZNMp3Bt9/m1Mt0Bo5HVAQAAKZGmAEAAKZGmAEAAKZGmAEAAKZGmAEAAKZGmAEAAKbGo9kAANfn4SENHfrnuouzWDxUtepQ6zocix4GALg+b2/pzTedXUWhubl5q25d89RrdlxmAgAApsbIDADA9RmGdO5cznqlSpLF4tx6CmAYhjIycur19Kwki4vXa3aEGQCAQ4WNWWP3scdf7pqzkpoqVamSs26K6QxStW1bTr1MZ+B4XGYCAACmRpgBAACmRpgBAACmRpgBAACmRpgBAACmRpgBAACmxqPZAADX5+EhDRjw57qLs1g8FBQ0wLoOx6KHAQCuz9tbWrTI2VUUmpubtxo0WOTsMm4YXGYCAACmxsgMAMD1GUbOtwBLkp+fKaYzyM7OqdfNzY/pDByMkRkAgOtLTZXKlctZckONC8vOTtXmzeW0eXM5a6iB4xBmAACAqRFmAACAqRFmAACAqRFmAACAqRFmAACAqfFoNkwhbMwaZ5cAAHBRhBkAgOtzd5d69vxz3eW5q3LlntZ1OBZhBgDg+nx8pI8+cnYVhebu7qNGjcxTr9lxzwwAADA1p4aZuLg4tWrVSuXLl1eVKlX0wAMP6PDhwzb7GIah2NhYVa1aVb6+voqOjtaBAwecVDEAAHA1Tg0zmzZt0rBhw7Rjxw6tW7dOmZmZiomJ0eXLl637TJs2TTNmzNDs2bO1c+dOBQcHq0uXLrp48aITKwcAlKrLl3PmY7JYctZdXFbWZW3caNHGjRZlZbl+vWbn1HtmvvzyS5vXCxcuVJUqVbR792516NBBhmFo5syZGjdunLp37y5JWrx4sYKCgrR06VI98cQTzigbAAC4EJe6ZyY5OVmSdPPNN0uSjh07poSEBMXExFj38fb2VseOHbVt27Z8z5GWlqaUlBSbBQAAlF0uE2YMw9Azzzyj2267TbfccoskKSEhQZIUFBRks29QUJB129/FxcUpMDDQutSoUcOxhQMAAKdymTAzfPhw/fDDD1q2bFmebRaLxea1YRh52nKNHTtWycnJ1iU+Pt4h9QIAANfgEt8z8+STT+rzzz/XN998o+rVq1vbg4ODJeWM0ISEhFjbExMT84zW5PL29pa3t7djCwYAAC7DqSMzhmFo+PDhWrlypb7++muFh4fbbA8PD1dwcLDWrVtnbUtPT9emTZsUFRVV2uUCAAAX5NSRmWHDhmnp0qX67LPPVL58eet9MIGBgfL19ZXFYtHIkSM1ZcoURUREKCIiQlOmTJGfn5/69OnjzNIBAKUgd14278x0za3VUpL0z4lrlebhVeCxx1/u6tDars9dN998j3UdjmUxDMNw2ptf476XhQsXauDAgZJyRm8mTZqkt99+WxcuXFBkZKTefPNN603CBUlJSVFgYKCSk5MVEBBQUqWjlDHRJICicm6YQXEV5e+3U0dmCpOjLBaLYmNjFRsb6/iCAACA6bjM00wAAAD2cImnmQAAuB7f9KvaPbuvJKnF8CW64uVT4DHFuTxd3EtUWVmXtXVrFUlSu3aJcnf3L9b5cH2EGQCAKfhlpDm7hCLJzk51dgk3DC4zAQAAUyPMAAAAU+MyE0oNj1cDAByBkRkAAGBqhBkAAGBqXGYCALi8bItFO2rcYl13fW4KDOxoXYdjEWYAAC4vzdNbD/d52dllFJq7u6+aN9/o7DJuGMRFAABgaoQZAABgaoQZAIDL802/qt2v99Hu1/vIN/2qs8spUM50BpW1dWtlZWVddnY5ZR73zAAATKHilRRnl1AkGRnnnF3CDYORGQAAYGqEGQAAYGqEGQAAYGqEGQAAYGqEGQAAYGo8zQQAcHnZFou+D46wrrs+N5Uv39K6DscizAAAXF6ap7fuH/Cas8soNHd3X7VosdPZZdwwiIsAAMDUCDMAAMDUCDMAAJfnk3FVW+b+n7bM/T/5ZJhhOoNUbd8epu3bw5SVlerscso87pkBALg8iyFVT0m0rrs+Q2lpJ6zrcCxGZgAAgKkxMoMiCRuzxtklAABgg5EZAABgaoQZAABgaoQZAABgatwzAwBweYZFOlKxpnXd9Vnk59fQug7HIswAAFzeVU8fxQye4+wyCs3d3U+tWx9wdhk3DC4zAQAAUyPMAAAAU+My0w2I74oBYDY+GVf1+eJnJEn3DZihq54+Tq7o+rKyUrV7dytJUosWO+Xu7ufkiso2wgwAwOVZDKnuHyet667PUGrqQes6HIvLTAAAwNQIMwAAwNQIMwAAwNQIMwAAwNQIMwAAwNR4mgkA4PIMi3QqoIp13fVZ5O0dal2HYxFmAAAu76qnj2775wJnl1Fo7u5+atv2uLPLuGFwmQkAAJgaYQYAAJgaYQYA4PK8M9L02eKn9dnip+WdkebscgqUlXVFu3e30u7drZSVdcXZ5ZR53DNjUsyvBOBG4mYYaprws3Xd0YrzO/b4y10lZevixV3/vyW7RGrCtTEyAwAATI0wAwAATI0wAwAATI0wAwAATI0wAwAATI2nmQAApvCHb4CzSygST89Kzi7hhkGYAQC4vCtePmoxYqmzyyg0d3d/tWt31tll3DC4zAQAAEyNMAMAAEyNMAMAcHneGWlavnSMli8dY5rpDPbujdbevdFMZ1AKuGcGAODy3AxDbeJ/tK67vmwlJ2+yrsOxCDMAAJSgsDFr5OV+VfO65LxuMOFLpWf5FOrYnHmdUFRcZgIAAKZGmAEAAKZGmAEAAKbGPTNOFDZmjbNLAAC4kOL8XbiR77chzAAATCHV09vZJRRJWqa56jUzwgwAwOVd8fJRw2c+cXYZhZae5aMn/meees2Oe2YAAICpEWYAAICpcZkJAODyvDPTNXfVFEnSPx98QWkeXk6u6Po83dI1vFlOvbP3vaCMbNeu1+w3HhNmAAAuzy07W3cc3WVdd3UWS7aaVtllXYdjcZkJAACYminCzJw5cxQeHi4fHx+1aNFCmzdvdnZJAADARbj8ZaYVK1Zo5MiRmjNnjtq1a6e3335bd999tw4ePKiaNWs6uzy++A4A4BJu5L9HLj8yM2PGDA0aNEiDBw9WgwYNNHPmTNWoUUNz5851dmkAAMAFuHSYSU9P1+7duxUTE2PTHhMTo23btjmpKgAA4Epc+jLTuXPnlJWVpaCgIJv2oKAgJSQk5HtMWlqa0tLSrK+Tk5MlSSkpKQ6pMTst1SHnBQD8KSv9qnJ/i2elpSrbcO0nhLLcr+ry5f+/npaq7CzXrrc4HPX3Nfe8hmEUuK9Lh5lcFovF5rVhGHnacsXFxWnSpEl52mvUqOGQ2gAApSMwd2XOP5xZRqF1ezV3zRz12itwpmPPf/HiRQUGBl53H5cOM5UqVZK7u3ueUZjExMQ8ozW5xo4dq2eeecb6Ojs7W+fPn1fFihWvGYAcJSUlRTVq1FB8fLwCAgJK9b3NgP65Pvrn+uif66N/ro2+uT5X6R/DMHTx4kVVrVq1wH1dOsx4eXmpRYsWWrdunR588EFr+7p163T//ffne4y3t7e8vW1nKq1QoYIjyyxQQEAA/2Cug/65Pvrn+uif66N/ro2+uT5X6J+CRmRyuXSYkaRnnnlG/fv3V8uWLdW2bVvNmzdPJ0+e1JAhQ5xdGgAAcAEuH2Z69+6tP/74Qy+++KLOnDmjW265Rf/9738VGhrq7NIAAIALcPkwI0lDhw7V0KFDnV1GkXl7e2vixIl5LnshB/1zffTP9dE/10f/XBt9c31m7B+LUZhnngAAAFyUS39pHgAAQEEIMwAAwNQIMwAAwNQIMwAAwNQIM8U0Z84chYeHy8fHRy1atNDmzZuvue/KlSvVpUsXVa5cWQEBAWrbtq2++uqrUqy29BWlf7Zs2aJ27dqpYsWK8vX1Vf369fXaa6+VYrWlryj981dbt26Vh4eHmjVr5tgCnagofbNx40ZZLJY8y08//VSKFZeuov7spKWlady4cQoNDZW3t7dq166tBQsWlFK1pa8o/TNw4MB8f34aNWpUihWXrqL+/CxZskRNmzaVn5+fQkJC9Oijj+qPP/4opWoLwYDdli9fbnh6ehrvvPOOcfDgQeOpp54y/P39jRMnTuS7/1NPPWVMnTrV+O6774wjR44YY8eONTw9PY09e/aUcuWlo6j9s2fPHmPp0qXGjz/+aBw7dsx4//33DT8/P+Ptt98u5cpLR1H7J1dSUpJRq1YtIyYmxmjatGnpFFvKito3GzZsMCQZhw8fNs6cOWNdMjMzS7ny0mHPz859991nREZGGuvWrTOOHTtmfPvtt8bWrVtLserSU9T+SUpKsvm5iY+PN26++WZj4sSJpVt4KSlq/2zevNlwc3MzZs2aZRw9etTYvHmz0ahRI+OBBx4o5cqvjTBTDK1btzaGDBli01a/fn1jzJgxhT5Hw4YNjUmTJpV0aS6hJPrnwQcfNPr161fSpbkEe/und+/exr/+9S9j4sSJZTbMFLVvcsPMhQsXSqE65ytq/3zxxRdGYGCg8ccff5RGeU5X3N89q1atMiwWi3H8+HFHlOd0Re2f6dOnG7Vq1bJpe/31143q1as7rMai4jKTndLT07V7927FxMTYtMfExGjbtm2FOkd2drYuXryom2++2RElOlVJ9M/evXu1bds2dezY0RElOpW9/bNw4UL9+uuvmjhxoqNLdJri/Ow0b95cISEh6tSpkzZs2ODIMp3Gnv75/PPP1bJlS02bNk3VqlVT3bp19eyzz+rKlSulUXKpKonfPfPnz1fnzp3L5DfN29M/UVFROnXqlP773//KMAz9/vvv+vjjj9W1a9fSKLlQTPENwK7o3LlzysrKyjN7d1BQUJ5Zvq/l1Vdf1eXLl/XQQw85okSnKk7/VK9eXWfPnlVmZqZiY2M1ePBgR5bqFPb0z88//6wxY8Zo8+bN8vAou/907embkJAQzZs3Ty1atFBaWpref/99derUSRs3blSHDh1Ko+xSY0//HD16VFu2bJGPj49WrVqlc+fOaejQoTp//nyZu2+muL+bz5w5oy+++EJLly51VIlOZU//REVFacmSJerdu7euXr2qzMxM3XfffXrjjTdKo+RCKbu/EUuJxWKxeW0YRp62/CxbtkyxsbH67LPPVKVKFUeV53T29M/mzZt16dIl7dixQ2PGjFGdOnX0yCOPOLJMpyls/2RlZalPnz6aNGmS6tatW1rlOVVRfnbq1aunevXqWV+3bdtW8fHxeuWVV8pcmMlVlP7Jzs6WxWLRkiVLrLMQz5gxQz179tSbb74pX19fh9db2uz93bxo0SJVqFBBDzzwgIMqcw1F6Z+DBw9qxIgRmjBhgu68806dOXNGzz33nIYMGaL58+eXRrkFIszYqVKlSnJ3d8+TZBMTE/Mk3r9bsWKFBg0apI8++kidO3d2ZJlOU5z+CQ8PlyQ1btxYv//+u2JjY8tcmClq/1y8eFG7du3S3r17NXz4cEk5f6AMw5CHh4fWrl2rO+64o1Rqd7Ti/Oz8VZs2bfTBBx+UdHlOZ0//hISEqFq1atYgI0kNGjSQYRg6deqUIiIiHFpzaSrOz49hGFqwYIH69+8vLy8vR5bpNPb0T1xcnNq1a6fnnntOktSkSRP5+/urffv2eumllxQSEuLwugvCPTN28vLyUosWLbRu3Tqb9nXr1ikqKuqaxy1btkwDBw7U0qVLXep6Y0mzt3/+zjAMpaWllXR5TlfU/gkICND+/fu1b98+6zJkyBDVq1dP+/btU2RkZGmV7nAl9bOzd+9el/glW9Ls6Z927drp9OnTunTpkrXtyJEjcnNzU/Xq1R1ab2krzs/Ppk2b9Msvv2jQoEGOLNGp7Omf1NRUubnZxgV3d3dJOb+jXYJTbjsuI3Ifb5s/f75x8OBBY+TIkYa/v7/1DvgxY8YY/fv3t+6/dOlSw8PDw3jzzTdtHgNMSkpy1kdwqKL2z+zZs43PP//cOHLkiHHkyBFjwYIFRkBAgDFu3DhnfQSHKmr//F1ZfpqpqH3z2muvGatWrTKOHDli/Pjjj8aYMWMMScYnn3zirI/gUEXtn4sXLxrVq1c3evbsaRw4cMDYtGmTERERYQwePNhZH8Gh7P231a9fPyMyMrK0yy11Re2fhQsXGh4eHsacOXOMX3/91diyZYvRsmVLo3Xr1s76CHkQZorpzTffNEJDQw0vLy/j1ltvNTZt2mTdNmDAAKNjx47W1x07djQk5VkGDBhQ+oWXkqL0z+uvv240atTI8PPzMwICAozmzZsbc+bMMbKyspxQeekoSv/8XVkOM4ZRtL6ZOnWqUbt2bcPHx8e46aabjNtuu81Ys2aNE6ouPUX92Tl06JDRuXNnw9fX16hevbrxzDPPGKmpqaVcdekpav8kJSUZvr6+xrx580q5Uucoav+8/vrrRsOGDQ1fX18jJCTE6Nu3r3Hq1KlSrvraLIbhKmNEAAAARcc9MwAAwNQIMwAAwNQIMwAAwNQIMwAAwNQIMwAAwNQIMwAAwNQIMwAAwNQIMwAAwNQIMwCcZtu2bXJ3d9ddd91l075x40ZZLBYlJSXlOaZZs2aKjY21adu7d6969eqloKAg+fj4qG7dunrsscd05MgRSdLx48cLNWMyAHMizABwmgULFujJJ5/Uli1bdPLkSbvOsXr1arVp00ZpaWlasmSJDh06pPfff1+BgYEaP358CVcMwBV5OLsAADemy5cv68MPP9TOnTuVkJCgRYsWacKECUU6R2pqqh599FHdc889WrVqlbU9PDxckZGR+Y7sSNKJEyc0fPhwbdmyRenp6QoLC9P06dN1zz33FOcjAXASRmYAOMWKFStUr1491atXT/369dPChQtV1KnivvrqK507d06jR4/Od3uFChXybR82bJjS0tL0zTffaP/+/Zo6darKlStX1I8AwEUwMgPAKebPn69+/fpJku666y5dunRJ69evV+fOnQt9jp9//lmSVL9+/evuFxYWZhOUTp48qR49eqhx48aSpFq1ahW1fAAuhJEZAKXu8OHD+u677/Twww9Lkjw8PNS7d28tWLCgSOcp6khOrhEjRuill15Su3btNHHiRP3www92nQeAayDMACh18+fPV2ZmpqpVqyYPDw95eHho7ty5WrlypS5cuKCAgABJUnJycp5jk5KSFBgYKEmqW7euJOmnn34q0vsPHjxYR48eVf/+/bV//361bNlSb7zxRjE/FQBnsRj2/q8NANghMzNT1atX1+jRoxUTE2OzrUePHnryySc1YMAAVahQQStWrFDPnj2t28+cOaOaNWvqs88+0z333KPLly8rLCxMt912m80NwLmSkpKued/MX40dO1Zr1qxhhAYwKe6ZAVCqVq9erQsXLmjQoEHWEZZcPXv21Pz58zV8+HA98cQTGjVqlDw8PNS0aVOdPn1a48aNU4MGDawhyN/fX++++6569eql++67TyNGjFCdOnV07tw5ffjhhzp58qSWL1+ep4aRI0fq7rvvVt26dXXhwgV9/fXXatCgQal8fgAlj5EZAKXq3nvvVXZ2ttasWZNn2549e9SiRQvt3r1bjRo10rRp07Rs2TIdP35cVapU0e233664uDgFBwfbHLdr1y7FxcVp8+bNSklJUY0aNXTHHXfoueeeU506dfK8z5NPPqkvvvhCp06dUkBAgO666y699tprqlixosM+NwDHIcwAAABT4wZgAABgaoQZAABgaoQZAABgaoQZAABgaoQZAABgaoQZAABgaoQZAABgaoQZAABgaoQZAABgaoQZAABgaoQZAABgaoQZAABgav8PMMDigQdFxYYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Confidence interval for the permutation test: [0.315 - 0.674]\n",
      "AUC R: 0.7356770833333333\n",
      "p-value: 0.01\n"
     ]
    }
   ],
   "source": [
    "permutR=[]\n",
    "for j in range(1000):\n",
    "    predPer = Parallel(n_jobs=multiprocessing.cpu_count())(delayed(permutation)(dataR, lvhR, i) for i in np.unique(annotations[\"ID\"]))\n",
    "    pred=pd.DataFrame()\n",
    "    for k in range(len(predPer)):\n",
    "        pred=pd.concat([pred,predPer[k]])\n",
    "    permutR.append(metrics.roc_auc_score(pred[\"Label\"],pred[0]))\n",
    "sorted_permut = np.array(permutR)\n",
    "sorted_permut.sort()\n",
    "confidence_lower = sorted_permut[int(0.05 * len(sorted_permut))]\n",
    "confidence_upper = sorted_permut[int(0.95 * len(sorted_permut))]\n",
    "plt.hist(permutR, bins=30)\n",
    "plt.title('Histogram of H0')\n",
    "plt.axvline(x=confidence_upper, ymax=0.7, color='r',linestyle=\"--\")\n",
    "plt.axvline(x=metrics.roc_auc_score(predR['Label'],predR[0]), ymax=0.8, color=\"y\",linestyle=\"--\")\n",
    "plt.ylabel(\"Frecuency\")\n",
    "plt.xlabel(\"AUC's\")\n",
    "plt.text(0.6,68,\"Quantile 0.95\")\n",
    "plt.text(0.7,76,f\"AUC= {round(metrics.roc_auc_score(predR['Label'],predR[0]),2)}\")\n",
    "plt.show()\n",
    "print(\"Confidence interval for the permutation test: [{:0.3f} - {:0.3}]\".format(confidence_lower, confidence_upper))    \n",
    "print(f\"AUC R: {metrics.roc_auc_score(predR['Label'],predR[0])}\")\n",
    "print(f\"p-value: {len(np.where(np.array(permutR)>=metrics.roc_auc_score(predR['Label'],predR[0]))[0])/1000}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "289e71ce",
   "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}