inhouse_pertspectra_metrics.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b8270f66-7a27-4263-aa97-25abeb2cf34d",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "529f64f3-c83e-481d-872c-df89b9dd063b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "\n",
    "from utils import (\n",
    "    auprc,\n",
    "    factor_enrichment_gsea,\n",
    "    generate_k_fold,\n",
    "    get_gprofiler,\n",
    "    perturbation_signal_recovery,\n",
    "    read_aws_csv,\n",
    "    read_aws_h5ad,\n",
    "    read_aws_pickle,\n",
    "    retrieve_stringdb_neighbors,\n",
    "    set_seed,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8adbe3c2-a83b-4a8a-95ec-5f4c103f316c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in trained model outputs generated from ./PertSpectra_load_checkpoints/pertspectra_inhouse.ipynb\n",
    "adata = read_aws_h5ad(\n",
    "    \"s3://pert-spectra/PertSpectra_checkpoints/pertspectra_inhouse/fold_0.h5ad\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b7bf773-5e0c-43b8-8df2-4ce3376be959",
   "metadata": {},
   "source": [
    "# Reconstruction"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b4b092b-2120-48c6-a492-84656e8b7169",
   "metadata": {},
   "source": [
    "## Spearman Coefficient\n",
    "- Spearman correlation between predicted and observed expression of test set\n",
    "- Compute correlation for the mean expression aggregated on perturbation, across all kfolds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "d62bb784-fd4f-4a86-acd4-fd54e784aa73",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_folds = 5\n",
    "model_adatas = []\n",
    "for n in range(0, n_folds):\n",
    "    # new adata\n",
    "    adata_n = adata.copy()\n",
    "    # load model from checkpoint\n",
    "    s3_dir = \"s3://pert-spectra/PertSpectra_checkpoints/\"\n",
    "    experiment_name = \"pertspectra_inhouse/\"\n",
    "    model_name = f\"fold_{n}.h5ad\"\n",
    "    m_adata = read_aws_h5ad(s3_dir + experiment_name + model_name)\n",
    "    model_adatas.append(m_adata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "ac5cff31-fd20-4c61-8d0b-7eb7d25a213b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# iterate through all models and get the losses and correlations\n",
    "# take the mean loss and correlation for the test set\n",
    "from scipy.stats import spearmanr\n",
    "\n",
    "test_corr_singles = []\n",
    "test_corr_combos = []\n",
    "\n",
    "set_seed(0)\n",
    "\n",
    "for n in range(n_folds):\n",
    "    adata_n = model_adatas[n]\n",
    "    train_idx, val_idx, test_idx = generate_k_fold(\n",
    "        adata_n, adata_n.X, adata_n.obs[\"condition\"], fold_idx=n\n",
    "    )\n",
    "    loss_weights = np.ones(adata_n.shape[0])\n",
    "    adata_test = adata_n[test_idx]\n",
    "\n",
    "    for pert in adata_test.obs[\"condition\"].unique():\n",
    "        hold_idx = [\n",
    "            i\n",
    "            for i, x in enumerate(\n",
    "                adata_test[adata_test.obs[\"Treatment\"] == \"TNFA+\"].obs[\"condition\"]\n",
    "            )\n",
    "            if x == pert\n",
    "        ]\n",
    "        recon = adata_test.uns[\"recon\"]\n",
    "        # correlation\n",
    "        mean_reconstruction = recon.mean(axis=0)\n",
    "        mean_observed = np.squeeze(\n",
    "            np.array(\n",
    "                adata_test[adata_test.obs[\"Treatment\"] == \"TNFA+\"][hold_idx].X.mean(\n",
    "                    axis=0\n",
    "                )\n",
    "            )\n",
    "        )\n",
    "        if (\"+\" in pert) and (\"ctrl\" not in pert):\n",
    "            test_corr_combos.append(\n",
    "                [pert, spearmanr(mean_reconstruction, mean_observed)[0], n]\n",
    "            )\n",
    "        else:\n",
    "            test_corr_singles.append(\n",
    "                [pert, spearmanr(mean_reconstruction, mean_observed)[0], n]\n",
    "            )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "62812bb6-f82a-4c0b-ad92-18ec1bb9084d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnQAAAHHCAYAAAA238WJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABRgklEQVR4nO3deZxO9f//8ec1ZjUrhlkwi7HvUTH2UJNURLZUg7aPDyFKJB9Li4+yli0lJD6kkDaUtcUSIoSQrTAaYYxlhpn3749+c31dZoYZzTXXHPO4325zu7nOeV/nvF7nOteZp3Od64zNGGMEAAAAy3JzdQEAAAD4Zwh0AAAAFkegAwAAsDgCHQAAgMUR6AAAACyOQAcAAGBxBDoAAACLI9ABAABYHIEOAADA4m75QDdr1izZbDYdOnQoz5Z56NAh2Ww2zZo1K8+WmVPNmjVT9erV8329+Gey2g+bNWumZs2auaymaznjvfJP2Ww2DR8+PE+WlfG+HTNmTJ4sz0q6deumqKgoV5dR4Lny2I68lXE827x5s9PXVVD2m5sKdAcOHNAzzzyjcuXKydvbWwEBAWrYsKEmTpyoixcv5nWNLjNv3jxNmDDB1WU43euvv64lS5a4ugzkwIULFzR8+HCtWbPG1aUAwE07duyYhg8frm3btrm6lBwr6JnAPbdP+OKLL9ShQwd5eXnp8ccfV/Xq1ZWamqrvvvtOL7zwgnbt2qXp06c7o9Z8N2/ePO3cuVP9+vVzmB4ZGamLFy/Kw8PDNYXlsddff10PP/yw2rZt6+pSCpUVK1bk+jkXLlzQiBEjJKlAnd0DbgW32rG9IDt27JhGjBihqKgo1a5d29Xl5EhBzwS5CnQHDx5U586dFRkZqVWrViksLMw+r1evXtq/f7+++OKLf1yUMUaXLl2Sj49PpnmXLl2Sp6en3Nxc92mxzWaTt7e3y9aP/HPlyhWlp6fL09Mzz5ftjGXCGs6fPy9fX19Xl1GguWIb3erH9vT0dKWmprq0x4xjqrNkZIT8VFD2m1ylojfeeEPJycmaMWOGQ5jLUL58efXt29f++MqVK3rllVcUExMjLy8vRUVF6aWXXlJKSorD86KionT//fdr+fLluv322+Xj46N33nlHa9askc1m0/z58/Xyyy+rdOnSKlq0qJKSkiRJGzdu1L333qvAwEAVLVpUTZs21ffff3/DPj799FO1bt1a4eHh8vLyUkxMjF555RWlpaXZxzRr1kxffPGFDh8+LJvNJpvNZr8GJbvPy1etWqXGjRvL19dXQUFBatOmjXbv3u0wZvjw4bLZbNq/f7+6deumoKAgBQYGqnv37rpw4cINa8+wZcsWNWjQQD4+PoqOjta0adMyjUlJSdGwYcNUvnx5eXl5qWzZsho4cKDD9rfZbDp//rxmz55t77Nbt276+eefZbPZtHTpUod12mw21alTx2E9rVq1Ur169RymffXVV/Zt4e/vr9atW2vXrl2ZatyzZ48efvhhFS9eXN7e3rr99tsd1in937UQ33//vfr376+SJUvK19dXDz30kP78888bbqtu3brJz89Pv/32m+Li4uTr66vw8HCNHDlSxhj7uKuvsZowYYJ9v/3ll19yXKsk7dq1S82bN5ePj4/KlCmjV199NcsDWFbX0F26dEnDhw9XxYoV5e3trbCwMLVr104HDhzQoUOHVLJkSUnSiBEj7K/X1deY5XWNWfn555/VrVs3+yUXoaGh6tGjh06dOuUwLjf7ekpKip577jmVLFlS/v7+evDBB/X777/nqB5JOnnypJ544gmFhITI29tbtWrV0uzZs7MdP378eEVGRsrHx0dNmzbVzp07HeafOHFC3bt3V5kyZeTl5aWwsDC1adMm0/WFOdnPM/a/AwcO6L777pO/v7+6du2q3r17y8/PL8v3fZcuXRQaGupwTMrpe2rJkiWqXr26vL29Vb16dS1evDgnmzDH61m1apXc3Nz0n//8x+F58+bNk81m09SpU+3TbDabevfurblz56pSpUry9vZW3bp1tW7dOofnZuwrv/zyix555BEVK1ZMjRo1ss//8MMPVbduXfn4+Kh48eLq3Lmzjh496rCMffv2qX379goNDZW3t7fKlCmjzp076+zZs/YxX3/9tRo1aqSgoCD5+fmpUqVKeumll+zzXXlsz7g+Oq+O7ddu/2rVqsnLy0vLli3LtoaM38UrVqxQ7dq15e3trapVq2rRokWZxp45c0b9+vVT2bJl5eXlpfLly2v06NEOx5HsjqlTpkzRHXfcIUnq3r27/ViWsd2joqLUrVu3LLfR1cfMG2UE6e9PNZ555hmVKFFCAQEBevzxx3X69GmH5d4KmUCSZHKhdOnSply5cjkeHx8fbySZhx9+2EyePNk8/vjjRpJp27atw7jIyEhTvnx5U6xYMTNo0CAzbdo0s3r1arN69WojyVStWtXUrl3bjBs3zowaNcqcP3/erFy50nh6eprY2FgzduxYM378eFOzZk3j6elpNm7caF/2zJkzjSRz8OBB+7S2bduajh07mjfffNNMnTrVdOjQwUgyzz//vH3MihUrTO3atU1wcLCZM2eOmTNnjlm8eLExxpiDBw8aSWbmzJn28V9//bVxd3c3FStWNG+88YYZMWKECQ4ONsWKFXNY97Bhw4wkc9ttt5l27dqZKVOmmCeffNJIMgMHDrzhNm3atKkJDw83pUqVMr179zZvvfWWadSokZFkZsyYYR+XlpZm7rnnHlO0aFHTr18/884775jevXsbd3d306ZNG/u4OXPmGC8vL9O4cWN7nz/88INJS0szQUFBZsCAAfax48ePN25ubsbNzc2cPXvWvp6AgACHbffBBx8Ym81m7r33XvP222+b0aNHm6ioKBMUFOSwLXbu3GkCAwNN1apVzejRo82kSZNMkyZNjM1mM4sWLcr0Gt52222mefPm5u233zYDBgwwRYoUMR07drzhNouPjzfe3t6mQoUK5rHHHjOTJk0y999/v5Fkhg4dah+X8bpWrVrVlCtXzvz3v/8148ePN4cPH85xrcePHzclS5Y0xYoVM8OHDzdvvvmmqVChgqlZs2am/bBp06amadOm9sdXrlwxLVq0MJJM586dzaRJk8yoUaNM8+bNzZIlS0xycrKZOnWqkWQeeugh++u1ffv2XG3P3NSYlTFjxpjGjRubkSNHmunTp5u+ffsaHx8fc+edd5r09HT7uNzs648++qiRZB555BEzadIk065dO3s9w4YNu249Fy5cMFWqVDEeHh7mueeeM2+99ZZp3LixkWQmTJiQ6fWtUaOGiYqKMqNHjzYjRowwxYsXNyVLljQnTpywj23QoIEJDAw0L7/8snnvvffM66+/bu666y6zdu1a+5ic7ufx8fHGy8vLxMTEmPj4eDNt2jTzwQcfmHXr1hlJ5qOPPnLo5/z588bX19f06tUr1+tavny5cXNzM9WrVzfjxo0zQ4YMMYGBgaZatWomMjLyutsxN+vp1auXcXd3N1u2bDHGGHPs2DFTvHhx07JlS4d9QJKpXr26CQ4ONiNHjjSjR482kZGRxsfHx+zYscM+LmNfqVq1qmnTpo2ZMmWKmTx5sjHGmFdffdXYbDbTqVMnM2XKFPuxNSoqypw+fdoYY0xKSoqJjo424eHh5tVXXzXvvfeeGTFihLnjjjvMoUOHjDF/vz88PT3N7bffbiZOnGimTZtmnn/+edOkSZNM+8itcGzP2P5VqlQxJUuWNCNGjDCTJ082P/30U7Y1REZGmooVK5qgoCAzaNAgM27cOFOjRg3j5uZmVqxYYR93/vx5U7NmTVOiRAnz0ksvmWnTppnHH3/c2Gw207dv30zb89pj6qFDh8zIkSONJPP000/bj2UHDhyw1xEfH5/lNrr6mHm9jJDxe6NGjRqmcePG5q233jK9evUybm5upkmTJg77qZUzwdVyHOjOnj1rJGXaYbKzbds2I8k8+eSTDtOff/55I8msWrXKPi0yMtJIMsuWLXMYm/FilStXzly4cME+PT093VSoUMHExcU5vCgXLlww0dHR5u6777ZPyyrQXb2sDM8884wpWrSouXTpkn1a69atszwIZvXi1a5d25QqVcqcOnXKPm379u3Gzc3NPP744/ZpGS9ejx49HJb50EMPmRIlSmRa17WaNm1qJJmxY8fap6WkpNjXn5qaaoz5O6i5ubmZb7/91uH506ZNM5LM999/b5/m6+ub5ZundevW5s4777Q/bteunWnXrp0pUqSI+eqrr4wxxmzdutVIMp9++qkxxphz586ZoKAg89RTTzks68SJEyYwMNBheosWLUyNGjUctnl6erpp0KCBqVChgn1axmt47S+L5557zhQpUsScOXPmutss4z8Wzz77rMN6WrdubTw9Pc2ff/5pjPm/1zUgIMCcPHnSYRk5rbVfv35GksN/Kk6ePGkCAwNvGOjef/99I8mMGzcuUw8Zff/555/Zhhxn1JiVrN4///vf/4wks27dOvu0nO7rGceKf//73w7jHnnkkRwFugkTJhhJ5sMPP7RPS01NNbGxscbPz88kJSUZY/7v9fXx8TG///67fezGjRuNJPPcc88ZY4w5ffq0kWTefPPNbNeZm/08Y/8bNGiQw9j09HRTunRp0759e4fpH330kcO2zM26ateubcLCwhzeEytWrDCSbhjocrOe8+fPm/Lly5tq1aqZS5cumdatW5uAgABz+PBhh+dKMpLM5s2b7dMOHz5svL29zUMPPWSflrGvdOnSxeH5hw4dMkWKFDGvvfaaw/QdO3YYd3d3+/SffvrJSDILFy7Mtr/x48cbSfb3e1ZutWO7JOPm5mZ27dp1w/Ub83+/iz/55BP7tLNnz5qwsDBz22232ae98sorxtfX1/z6668Ozx80aJApUqSIOXLkiDHm+sfUH3/8MdO2vrqO3AS6azOCMf/3e6Nu3br2bWeMMW+88YbD7yxjrJ0Jrpbjj1wzTmH6+/vnaPyXX34pSerfv7/D9AEDBkhSpmvtoqOjFRcXl+Wy4uPjHa6n27Ztm/bt26dHHnlEp06dUmJiohITE3X+/Hm1aNFC69atu+7HR1cv69y5c0pMTFTjxo114cIF7dmzJ0f9Xe348ePatm2bunXrpuLFi9un16xZU3fffbd9W1ztX//6l8Pjxo0b69SpUw6nirPj7u6uZ555xv7Y09NTzzzzjE6ePKktW7ZIkhYuXKgqVaqocuXK9u2TmJio5s2bS5JWr159w/U0btxYW7du1fnz5yVJ3333ne677z7Vrl1b3377rSTp22+/lc1ms3888vXXX+vMmTPq0qWLw3qLFCmievXq2df7119/adWqVerYsaP9NUhMTNSpU6cUFxenffv26Y8//nCo5+mnn5bNZnOoLy0tTYcPH75hL5LUu3dv+78zPopITU3VN9984zCuffv29o82c1vrl19+qfr16+vOO++0P79kyZLq2rXrDev75JNPFBwcrGeffTbTvKv7zkp+1Sg5vn8uXbqkxMRE1a9fX5K0devWTONvtK9nvD/69OnjMO7aC4+z8+WXXyo0NFRdunSxT/Pw8FCfPn2UnJystWvXOoxv27atSpcubX985513ql69evY6fHx85OnpqTVr1mT6aCZDTvfzq/Xs2dPhsc1mU4cOHfTll18qOTnZPn3BggUqXbp0rt9TGceh+Ph4BQYG2pd39913q2rVqjfcjrnpqWjRopo1a5Z2796tJk2a6IsvvtD48eMVERGRabmxsbGqW7eu/XFERITatGmj5cuXO3ykJWXeVxYtWqT09HR17NjRoabQ0FBVqFDBXlNGv8uXL8/2Y6qgoCBJf3+8ltPLC26FY3vTpk1z9PpnCA8P10MPPWR/nPEx5U8//aQTJ07Ya2jcuLGKFSvmUEPLli2VlpaW6SP1a4+pee3ajHC1p59+2uHLCj179pS7u7vDa2f1TJAhx1+KCAgIkPR3szlx+PBhubm5qXz58g7TQ0NDFRQUlOmXcHR0dLbLunbevn37JP39Imbn7NmzKlasWJbzdu3apZdfflmrVq3KtLGuvt4ipzJ6qVSpUqZ5VapU0fLlyzNd4HvtgS+j1tOnT9u3dXbCw8MzXSxcsWJFSX9/ll+/fn3t27dPu3fvzvZNdPLkyRt09fcOdeXKFa1fv15ly5bVyZMn1bhxY+3atcsh0FWtWtW+02a8NhkHl2tl9LZ//34ZYzR06FANHTo02xqv/sV7vW12I25ubipXrpzDtKu32dWu3d9yU+vhw4czXU8oZb1vXOvAgQOqVKmS3N1z/eXzfKtR+js8jhgxQvPnz8+0H2X1/rnRvp5xrIiJibmpeg4fPqwKFSpk+qJUlSpV7POvVqFChUzLqFixoj766CNJkpeXl0aPHq0BAwYoJCRE9evX1/3336/HH39coaGhknK+n2dwd3dXmTJlMo3r1KmTJkyYoKVLl+qRRx5RcnKyvvzySz3zzDP2EJ/TdWX0mVV/lSpVyjJsXy23PTVs2FA9e/bU5MmTFRcXpx49emT5vOy294ULF/Tnn3/at6mU9bHeGJPlMiTZf1FHR0erf//+GjdunObOnavGjRvrwQcf1KOPPmoPe506ddJ7772nJ598UoMGDVKLFi3Url07Pfzww9l+ye5WOLZf73drVsqXL5/pP5BX1xAaGqp9+/bp559/dloNuXW95V+77/j5+SksLMzhuG/1TJAhV4EuPDw808XDN3KjMwsZskvXWc3L+N/Vm2++me3Xnf38/LKcfubMGTVt2lQBAQEaOXKkYmJi5O3tra1bt+rFF1906rdvrlakSJEsp5urLtL/J9LT01WjRg2NGzcuy/lly5a94TJuv/12eXt7a926dYqIiFCpUqVUsWJFNW7cWFOmTFFKSoq+/fZbh//NZWy/OXPmOByoM2SElYxxzz//fLZnZq/9z4Czt1mG7Pa33NSa3/Kzxo4dO+qHH37QCy+8oNq1a8vPz0/p6em69957s3z/5Nfrlpf69eunBx54QEuWLNHy5cs1dOhQjRo1SqtWrdJtt92W4/08g5eXV5ahoX79+oqKitJHH32kRx55RJ999pkuXryoTp062cfkdl03K7frSUlJsd8P8cCBA7pw4YKKFi36j2rI6r1ns9n01VdfZbkfXX2cHzt2rLp166ZPP/1UK1asUJ8+fTRq1Cht2LBBZcqUkY+Pj9atW6fVq1friy++0LJly7RgwQI1b95cK1asyHY/za2Cdmy/3u/Wf1LD3XffrYEDB2Y5PyMA3mwN2eWGtLS0LLfvP+nxVsoEuToS3H///Zo+fbrWr1+v2NjY646NjIxUenq69u3bZ/+fsiQlJCTozJkzioyMzM2qHWT8Tz4gIEAtW7bM1XPXrFmjU6dOadGiRWrSpIl9+sGDBzONzWkYzehl7969mebt2bNHwcHBefr1+2PHjmVK97/++qsk2b91ExMTo+3bt6tFixY37CO7+Z6enrrzzjv17bffKiIiQo0bN5b095m7lJQUzZ07VwkJCQ7bMeO1KVWq1HVfm4yzZR4eHrl+DW9Genq6fvvtN4cDzbXbLDu5qTUyMtJ+puNqWe0b14qJidHGjRt1+fLlbO9nlN1rlV81nj59WitXrtSIESMcvuWY1fJyKuNYkXGGMjf1ZDz/559/Vnp6ukNoyvio5NpjTVa1/vrrr5n2g5iYGA0YMEADBgzQvn37VLt2bY0dO1YffvhhjvfznOjYsaMmTpyopKQkLViwQFFRUfaPsDPqyMm6Mvr8J/tfTtaTYdiwYdq9e7fGjBmjF198UYMGDdJbb72VaVx227to0aI3/BguJiZGxhhFR0dnCglZqVGjhmrUqKGXX35ZP/zwgxo2bKhp06bp1VdflfT3mfoWLVqoRYsWGjdunF5//XUNGTJEq1evzrJnqx/bb0bG2f6rl51VDcnJyf9o379e7cWKFdOZM2cyTT98+HCmT1puZN++fbrrrrvsj5OTk3X8+HHdd999km6NTJAhV7ctGThwoHx9ffXkk08qISEh0/wDBw5o4sSJkmTfWNfeVTnjfxWtW7e+mXolSXXr1lVMTIzGjBnjcO1JhuvdyiIjBV+delNTUzVlypRMY319fXN0ujUsLEy1a9fW7NmzHXbCnTt3asWKFfZtkVeuXLmid955x/44NTVV77zzjkqWLGm/VqVjx476448/9O6772Z6/sWLF+3XxUl/95nVm0f6O7xt3LhRq1evtge64OBgValSRaNHj7aPyRAXF6eAgAC9/vrrunz5cqblZbw2pUqVUrNmzfTOO+/o+PHj2Y7LS5MmTbL/2xijSZMmycPDQy1atLju83JT63333acNGzZo06ZNDvPnzp17w/rat2+vxMREhzqvrleS/QzIta9XftWY1ftHyvw+z41WrVpJUqYwkNNl3nfffTpx4oQWLFhgn3blyhW9/fbb8vPzU9OmTR3GL1myxOH6zE2bNmnjxo32Oi5cuKBLly45PCcmJkb+/v7220LkdD/PiU6dOiklJUWzZ8/WsmXL1LFjR4f5OV3X1ceha2/VkXHrnevJTU8bN27UmDFj1K9fPw0YMEAvvPCCJk2alOl6RUlav369w8e9R48e1aeffqp77rnnhmfF2rVrpyJFimjEiBGZ9jljjP1WOUlJSbpy5YrD/Bo1asjNzc3+mv3111+Zlp/xCc+1t/vIYPVj+804duyYw61ukpKS9MEHH6h27dr2M7cdO3bU+vXrtXz58kzPP3PmTKbXIisZgSar3z0xMTHasGGDUlNT7dM+//zzTLeqyYnp06c77M9Tp07VlStX7O/3WyETZMjVGbqYmBjNmzdPnTp1UpUqVRz+UsQPP/yghQsX2u8dU6tWLcXHx2v69On2U5qbNm3S7Nmz1bZtW4fEnFtubm5677331KpVK1WrVk3du3dX6dKl9ccff2j16tUKCAjQZ599luVzGzRooGLFiik+Pl59+vSRzWbTnDlzsjytWbduXS1YsED9+/fXHXfcIT8/Pz3wwANZLvfNN99Uq1atFBsbqyeeeEIXL17U22+/rcDAwDz7W5QZwsPDNXr0aB06dEgVK1bUggULtG3bNk2fPt1+Zuexxx7TRx99pH/9619avXq1GjZsqLS0NO3Zs0cfffSR/Z5/GX1+8803GjdunMLDwxUdHW2/xqpx48Z67bXXdPToUYfg1qRJE73zzjuKiopyuDYoICBAU6dO1WOPPaY6deqoc+fOKlmypI4cOaIvvvhCDRs2tAeWyZMnq1GjRqpRo4aeeuoplStXTgkJCVq/fr1+//13bd++Pc+2mbe3t5YtW6b4+HjVq1dPX331lb744gu99NJLObpYN6e1Dhw4UHPmzNG9996rvn37ytfXV9OnT7efRbqexx9/XB988IH69++vTZs2qXHjxjp//ry++eYb/fvf/1abNm3k4+OjqlWrasGCBapYsaKKFy+u6tWrq3r16vlSY0BAgJo0aaI33nhDly9fVunSpbVixYos/zebU7Vr11aXLl00ZcoUnT17Vg0aNNDKlSu1f//+HD3/6aef1jvvvKNu3bppy5YtioqK0scff6zvv/9eEyZMyPRFrvLly6tRo0bq2bOnUlJSNGHCBJUoUcL+8dGvv/6qFi1aqGPHjqpatarc3d21ePFiJSQkqHPnzvbtkNP9/Ebq1Kmj8uXLa8iQIUpJSXH4uDW36xo1apRat26tRo0aqUePHvrrr7/09ttvq1q1aln+5/dm1nPp0iXFx8erQoUKeu211yT9fV/Ezz77TN27d9eOHTsczj5Ur15dcXFx6tOnj/0eZBnPuZGYmBi9+uqrGjx4sA4dOqS2bdvK399fBw8e1OLFi/X000/r+eef16pVq9S7d2916NBBFStW1JUrVzRnzhwVKVJE7du3lySNHDlS69atU+vWrRUZGamTJ09qypQpKlOmjMM9765l5WP7zahYsaKeeOIJ/fjjjwoJCdH777+vhIQEzZw50z7mhRde0NKlS3X//ferW7duqlu3rs6fP68dO3bo448/1qFDhxQcHHzd9cTExCgoKEjTpk2Tv7+/fH19Va9ePUVHR+vJJ5/Uxx9/rHvvvVcdO3bUgQMHHM6M50Zqaqr9/bx3715NmTJFjRo10oMPPijp1sgEdrn6Tuz/9+uvv5qnnnrKREVFGU9PT+Pv728aNmxo3n77bYev+F6+fNmMGDHCREdHGw8PD1O2bFkzePBghzHG/P0V5datW2daT8ZXkrP7KvpPP/1k2rVrZ0qUKGG8vLxMZGSk6dixo1m5cqV9TFa3Lfn+++9N/fr1jY+PjwkPDzcDBw40y5cvN5LM6tWr7eOSk5PNI488YoKCghy+9p/VV5SNMeabb74xDRs2ND4+PiYgIMA88MAD5pdffnEYk/EV5Wu/Op9VnVlp2rSpqVatmtm8ebOJjY013t7eJjIy0kyaNCnT2NTUVDN69GhTrVo14+XlZYoVK2bq1q1rRowYYb+PnDHG7NmzxzRp0sT4+PgYSQ5fF09KSjJFihQx/v7+5sqVK/bpH374oZFkHnvssSzrXL16tYmLizOBgYHG29vbxMTEmG7dujncvsAYYw4cOGAef/xxExoaajw8PEzp0qXN/fffbz7++ONM2+bHH3/MtI5rX7OsxMfHG19fX3PgwAH7/ZtCQkLMsGHDTFpamn1cxuua3e0qclKrMcb8/PPPpmnTpsbb29uULl3avPLKK2bGjBk3vG2JMX9/fX7IkCH290xoaKh5+OGH7fdnMsaYH374wdStW9d4enpmuq1HXteYld9//9089NBDJigoyAQGBpoOHTqYY8eOZaolN/v6xYsXTZ8+fUyJEiWMr6+veeCBB8zRo0dzdNsSY4xJSEgw3bt3N8HBwcbT09PUqFEj0/vz6td37NixpmzZsvZ7MGbcy88YYxITE02vXr1M5cqVja+vrwkMDDT16tXLdL84Y3K2n2fsf9czZMgQI8mUL18+2zE5fU998sknpkqVKsbLy8tUrVrVLFq0yMTHx+foPnQ5WU/G7YKuvu2NMcZs3rzZuLu7m549e9qnSTK9evUyH374oalQoYLx8vIyt912W6b3bHb7ytU9NWrUyPj6+hpfX19TuXJl06tXL7N3715jjDG//fab6dGjh4mJiTHe3t6mePHi5q677jLffPONfRkrV640bdq0MeHh4cbT09OEh4ebLl26ONx641Y7tmds/5zK+F28fPlyU7NmTePl5WUqV66c5e/gc+fOmcGDB5vy5csbT09PExwcbBo0aGDGjBljv03IjY6pn376qalatapxd3fPtN3Hjh1rSpcubby8vEzDhg3N5s2bs71tSVb1ZWz3tWvXmqefftoUK1bM+Pn5ma5duzrcSsQYa2eCq9mMKcBXJgN5oFu3bvr4449veIYCQN6y2Wzq1atXjs9WFmbNmjVTYmJirr94mJeioqJUvXp1ff755y6rATfPdX8QFQAAAHmCQAcAAGBxBDoAAACL4xo6AAAAi+MMHQAAgMUR6AAAACwub/4IoMWlp6fr2LFj8vf3d8qfUgEAAHnPGKNz584pPDw8y7+XXJgQ6PT3nzrJyR+rBwAABc/Ro0cd/mpRYUSgk+x/Gujo0aMKCAhwcTUAACAnkpKSVLZs2Ux/4q8wItBJ9o9ZAwICCHQAAFgMl0vxpQgAAADLI9ABAABYHIEOAADA4gh0AAAAFkegAwAAsDgCHQAAgMUR6AAAACyOQAcAAGBxBDoAAACLI9ABAABYHIEOAADA4gh0AAAAFkegAwAAsDgCHQAAgMW5u7oAAPnryJEjSkxMdHUZ1xUcHKyIiAhXlwEAlkGgAwqRI0eOqHLlKrp48YKrS7kuH5+i2rNnN6EOAHKIQAcUIomJibp48YLq9RimgLAoV5eTpaTjh7Tx/RFKTEwk0AFADhX4QDd8+HCNGDHCYVqlSpW0Z88eSdKlS5c0YMAAzZ8/XykpKYqLi9OUKVMUEhLiinIBSwgIi1LxiEquLgMAkEcs8aWIatWq6fjx4/af7777zj7vueee02effaaFCxdq7dq1OnbsmNq1a+fCagEAAPJXgT9DJ0nu7u4KDQ3NNP3s2bOaMWOG5s2bp+bNm0uSZs6cqSpVqmjDhg2qX79+fpcKAACQ7yxxhm7fvn0KDw9XuXLl1LVrVx05ckSStGXLFl2+fFktW7a0j61cubIiIiK0fv16V5ULAACQrwr8Gbp69epp1qxZqlSpko4fP64RI0aocePG2rlzp06cOCFPT08FBQU5PCckJEQnTpzIdpkpKSlKSUmxP05KSnJW+QAAAE5X4ANdq1at7P+uWbOm6tWrp8jISH300Ufy8fG5qWWOGjUq0xctAAAArMoSH7leLSgoSBUrVtT+/fsVGhqq1NRUnTlzxmFMQkJCltfcZRg8eLDOnj1r/zl69KiTqwYAAHAeywW65ORkHThwQGFhYapbt648PDy0cuVK+/y9e/fqyJEjio2NzXYZXl5eCggIcPgBAACwqgL/kevzzz+vBx54QJGRkTp27JiGDRumIkWKqEuXLgoMDNQTTzyh/v37q3jx4goICNCzzz6r2NhYvuEKAAAKjQIf6H7//Xd16dJFp06dUsmSJdWoUSNt2LBBJUuWlCSNHz9ebm5uat++vcONhQEAAAqLAh/o5s+ff9353t7emjx5siZPnpxPFQEAABQslruGDgAAAI4IdAAAABZHoAMAALA4Ah0AAIDFEegAAAAsjkAHAABgcQQ6AAAAiyPQAQAAWByBDgAAwOIIdAAAABZHoAMAALA4Ah0AAIDFEegAAAAsjkAHAABgcQQ6AAAAiyPQAQAAWByBDgAAwOIIdAAAABZHoAMAALA4Ah0AAIDFEegAAAAsjkAHAABgcQQ6AAAAiyPQAQAAWByBDgAAwOIIdAAAABZHoAMAALA4Ah0AAIDFEegAAAAsjkAHAABgcQQ6AAAAiyPQAQAAWByBDgAAwOIIdAAAABZHoAMAALA4Ah0AAIDFEegAAAAsjkAHAABgcQQ6AAAAiyPQAQAAWByBDgAAwOIIdAAAABZHoAMAALA4Ah0AAIDFEegAAAAsjkAHAABgcQQ6AAAAiyPQAQAAWByBDgAAwOIIdAAAABZHoAMAALA4Ah0AAIDFEegAAAAsjkAHAABgcQQ6AAAAiyPQAQAAWByBDgAAwOIIdAAAABZnuUD33//+VzabTf369bNPu3Tpknr16qUSJUrIz89P7du3V0JCguuKBAAAyEeWCnQ//vij3nnnHdWsWdNh+nPPPafPPvtMCxcu1Nq1a3Xs2DG1a9fORVUCAADkL8sEuuTkZHXt2lXvvvuuihUrZp9+9uxZzZgxQ+PGjVPz5s1Vt25dzZw5Uz/88IM2bNjgwooBAADyh2UCXa9evdS6dWu1bNnSYfqWLVt0+fJlh+mVK1dWRESE1q9fn99lAgAA5Dt3VxeQE/Pnz9fWrVv1448/Zpp34sQJeXp6KigoyGF6SEiITpw4keXyUlJSlJKSYn+clJSUp/UCAADkpwJ/hu7o0aPq27ev5s6dK29v7zxZ5qhRoxQYGGj/KVu2bJ4sFwAAwBUKfKDbsmWLTp48qTp16sjd3V3u7u5au3at3nrrLbm7uyskJESpqak6c+aMw/MSEhIUGhqa5TIHDx6ss2fP2n+OHj2aD50AAAA4R4H/yLVFixbasWOHw7Tu3burcuXKevHFF1W2bFl5eHho5cqVat++vSRp7969OnLkiGJjY7NcppeXl7y8vJxeOwAAQH4o8IHO399f1atXd5jm6+urEiVK2Kc/8cQT6t+/v4oXL66AgAA9++yzio2NVf369V1RMgAAQL4q8IEuJ8aPHy83Nze1b99eKSkpiouL05QpU1xdFgAAQL6wZKBbs2aNw2Nvb29NnjxZkydPdk1BAAAALmTJQAfg1rd7925Xl3BdwcHBioiIcHUZACCJQAeggLl49pQkmx599FFXl3JdPj5FtWfPbkIdgAKBQAegQLl84Zwko9qPvKiS0ZVdXU6Wko4f0sb3RygxMZFAB6BAINABKJD8SkWoeEQlV5cBAJZQ4G8sDAAAgOsj0AEAAFgcgQ4AAMDiCHQAAAAWR6ADAACwOAIdAACAxRHoAAAALI5ABwAAYHEEOgAAAIsj0AEAAFgcgQ4AAMDiCHQAAAAWR6ADAACwOAIdAACAxRHoAAAALI5ABwAAYHEEOgAAAIsj0AEAAFgcgQ4AAMDiCHQAAAAWR6ADAACwOAIdAACAxRHoAAAALM7d1QUAt5IjR44oMTHR1WVka/fu3a4uAQDgBAQ6II8cOXJElStX0cWLF1xdyg1dTkl1dQkAgDxEoAPySGJioi5evKB6PYYpICzK1eVk6fiO9dq5dLquXLni6lIAAHmIQAfksYCwKBWPqOTqMrKUdPyQq0sAADgBX4oAAACwOAIdAACAxRHoAAAALI5r6GAZ3BIEAICsEehgCdwSBACA7BHoYAncEgQAgOwR6GAp3BIEAIDM+FIEAACAxRHoAAAALI5ABwAAYHEEOgAAAIsj0AEAAFgcgQ4AAMDiCHQAAAAWR6ADAACwOAIdAACAxRHoAAAALI5ABwAAYHEEOgAAAIsj0AEAAFgcgQ4AAMDiCHQAAAAWR6ADAACwOAIdAACAxbm7ugAAsKrdu3e7uoTrCg4OVkREhKvLAJAPCHQAkEsXz56SZNOjjz7q6lKuy8enqPbs2U2oAwoBAh0A5NLlC+ckGdV+5EWVjK7s6nKylHT8kDa+P0KJiYkEOqAQKPCBburUqZo6daoOHTokSapWrZr+85//qFWrVpKkS5cuacCAAZo/f75SUlIUFxenKVOmKCQkxIVVAygM/EpFqHhEJVeXAQAF/0sRZcqU0X//+19t2bJFmzdvVvPmzdWmTRvt2rVLkvTcc8/ps88+08KFC7V27VodO3ZM7dq1c3HVAAAA+afAn6F74IEHHB6/9tprmjp1qjZs2KAyZcpoxowZmjdvnpo3by5JmjlzpqpUqaINGzaofv36rigZAAAgXxX4M3RXS0tL0/z583X+/HnFxsZqy5Ytunz5slq2bGkfU7lyZUVERGj9+vUurBQAACD/FPgzdJK0Y8cOxcbG6tKlS/Lz89PixYtVtWpVbdu2TZ6engoKCnIYHxISohMnTmS7vJSUFKWkpNgfJyUlOat0AAAAp7PEGbpKlSpp27Zt2rhxo3r27Kn4+Hj98ssvN728UaNGKTAw0P5TtmzZPKwWAAAgf1ki0Hl6eqp8+fKqW7euRo0apVq1amnixIkKDQ1Vamqqzpw54zA+ISFBoaGh2S5v8ODBOnv2rP3n6NGjTu4AAADAeSwR6K6Vnp6ulJQU1a1bVx4eHlq5cqV93t69e3XkyBHFxsZm+3wvLy8FBAQ4/AAAAFhVgb+GbvDgwWrVqpUiIiJ07tw5zZs3T2vWrNHy5csVGBioJ554Qv3791fx4sUVEBCgZ599VrGxsXzDFQAAFBpOC3TlypXTjz/+qBIlSjhMP3PmjOrUqaPffvstR8s5efKkHn/8cR0/flyBgYGqWbOmli9frrvvvluSNH78eLm5ual9+/YONxYGAAAoLJwW6A4dOqS0tLRM01NSUvTHH3/keDkzZsy47nxvb29NnjxZkydPznWNAAAAt4I8D3RLly61/zvjY9EMaWlpWrlypaKiovJ6tQCALOzevdvVJVxXcHAwf2sWyAN5Hujatm0rSbLZbIqPj3eY5+HhoaioKI0dOzavVwsAuMrFs6ck2fToo4+6upTr8vEpqj17dhPqgH8ozwNdenq6JCk6Olo//vijgoOD83oVAIAbuHzhnCSj2o+8qJLRlV1dTpaSjh/SxvdHKDExkUAH/ENOu4bu4MGDzlo0ACCH/EpFqHhEJVeXAcDJnHrbkpUrV2rlypU6efKk/cxdhvfff9+ZqwYAACg0nBboRowYoZEjR+r2229XWFiYbDabs1YFAABQqDkt0E2bNk2zZs3SY4895qxVAAAAQE7801+pqalq0KCBsxYPAACA/89pge7JJ5/UvHnznLV4AAAA/H9O+8j10qVLmj59ur755hvVrFlTHh4eDvPHjRvnrFUDAAAUKk4LdD///LNq164tSdq5c6fDPL4gAQAAkHecFuhWr17trEUDAADgKk67hg4AAAD5w2ln6O66667rfrS6atUqZ60aAACgUHFaoMu4fi7D5cuXtW3bNu3cuVPx8fHOWi0AAECh47RAN378+CynDx8+XMnJyc5aLQAAQKGT79fQPfroo/wdVwAAgDyU74Fu/fr18vb2zu/VAgAA3LKc9pFru3btHB4bY3T8+HFt3rxZQ4cOddZqAQAACh2nBbrAwECHx25ubqpUqZJGjhype+65x1mrBQAAKHScFuhmzpzprEUDAADgKk4LdBm2bNmi3bt3S5KqVaum2267zdmrBAAAKFScFuhOnjypzp07a82aNQoKCpIknTlzRnfddZfmz5+vkiVLOmvVAAAAhYrTvuX67LPP6ty5c9q1a5f++usv/fXXX9q5c6eSkpLUp08fZ60WAACg0HHaGbply5bpm2++UZUqVezTqlatqsmTJ/OlCAAAgDzktDN06enp8vDwyDTdw8ND6enpzlotAABAoeO0QNe8eXP17dtXx44ds0/7448/9Nxzz6lFixbOWi0AAECh47RAN2nSJCUlJSkqKkoxMTGKiYlRdHS0kpKS9PbbbztrtQAAAIWO066hK1u2rLZu3apvvvlGe/bskSRVqVJFLVu2dNYqAQAACqU8P0O3atUqVa1aVUlJSbLZbLr77rv17LPP6tlnn9Udd9yhatWq6dtvv83r1QIAABRaeR7oJkyYoKeeekoBAQGZ5gUGBuqZZ57RuHHj8nq1AAAAhVaeB7rt27fr3nvvzXb+Pffcoy1btuT1agEAAAqtPA90CQkJWd6uJIO7u7v+/PPPvF4tAABAoZXnga506dLauXNntvN//vlnhYWF5fVqAQAACq08D3T33Xefhg4dqkuXLmWad/HiRQ0bNkz3339/Xq8WAACg0Mrz25a8/PLLWrRokSpWrKjevXurUqVKkqQ9e/Zo8uTJSktL05AhQ/J6tQAAAIVWnge6kJAQ/fDDD+rZs6cGDx4sY4wkyWazKS4uTpMnT1ZISEherxYAAKDQcsqNhSMjI/Xll1/q9OnT2r9/v4wxqlChgooVK+aM1QEAABRqTvtLEZJUrFgx3XHHHc5cBQAAQKHntL/lCgAAgPxBoAMAALA4Ah0AAIDFEegAAAAsjkAHAABgcQQ6AAAAiyPQAQAAWByBDgAAwOIIdAAAABZHoAMAALA4Ah0AAIDFEegAAAAsjkAHAABgcQQ6AAAAiyPQAQAAWByBDgAAwOIIdAAAABZHoAMAALA4Ah0AAIDFFfhAN2rUKN1xxx3y9/dXqVKl1LZtW+3du9dhzKVLl9SrVy+VKFFCfn5+at++vRISElxUMQAAQP4q8IFu7dq16tWrlzZs2KCvv/5aly9f1j333KPz58/bxzz33HP67LPPtHDhQq1du1bHjh1Tu3btXFg1AABA/nF3dQE3smzZMofHs2bNUqlSpbRlyxY1adJEZ8+e1YwZMzRv3jw1b95ckjRz5kxVqVJFGzZsUP369V1RNgAAQL4p8GfornX27FlJUvHixSVJW7Zs0eXLl9WyZUv7mMqVKysiIkLr1693SY0AAAD5qcCfobtaenq6+vXrp4YNG6p69eqSpBMnTsjT01NBQUEOY0NCQnTixIksl5OSkqKUlBT746SkJKfVDAAA4GyWOkPXq1cv7dy5U/Pnz/9Hyxk1apQCAwPtP2XLls2jCgEAAPKfZQJd79699fnnn2v16tUqU6aMfXpoaKhSU1N15swZh/EJCQkKDQ3NclmDBw/W2bNn7T9Hjx51ZukAAABOVeADnTFGvXv31uLFi7Vq1SpFR0c7zK9bt648PDy0cuVK+7S9e/fqyJEjio2NzXKZXl5eCggIcPgBAACwqgJ/DV2vXr00b948ffrpp/L397dfFxcYGCgfHx8FBgbqiSeeUP/+/VW8eHEFBATo2WefVWxsLN9wBQAAhUKBD3RTp06VJDVr1sxh+syZM9WtWzdJ0vjx4+Xm5qb27dsrJSVFcXFxmjJlSj5XCgAA4BoFPtAZY244xtvbW5MnT9bkyZPzoSIAAICCpcBfQwcAAIDrI9ABAABYHIEOAADA4gh0AAAAFkegAwAAsDgCHQAAgMUR6AAAACyOQAcAAGBxBDoAAACLI9ABAABYHIEOAADA4gh0AAAAFkegAwAAsDgCHQAAgMUR6AAAACyOQAcAAGBxBDoAAACLI9ABAABYHIEOAADA4gh0AAAAFkegAwAAsDgCHQAAgMUR6AAAACyOQAcAAGBxBDoAAACLI9ABAABYHIEOAADA4txdXQAAoHDbvXu3q0u4rpSUFHl5ebm6jOsKDg5WRESEq8uACxHoAAAucfHsKUk2Pfroo64u5fpsNskYV1dxXT4+RbVnz25CXSFGoAMAuMTlC+ckGdV+5EWVjK7s6nKydHzHeu1cOr1A15h0/JA2vj9CiYmJBLpCjEAHAHApv1IRKh5RydVlZCnp+CFJBbtGQOJLEQAAAJZHoAMAALA4Ah0AAIDFEegAAAAsjkAHAABgcQQ6AAAAiyPQAQAAWByBDgAAwOIIdAAAABZHoAMAALA4Ah0AAIDFEegAAAAsjkAHAABgcQQ6AAAAiyPQAQAAWByBDgAAwOIIdAAAABZHoAMAALA4Ah0AAIDFEegAAAAsjkAHAABgcQQ6AAAAiyPQAQAAWByBDgAAwOLcXV0ACoYjR44oMTHR1WVka/fu3a4uAQAKtIJ+nAwODlZERISry7hlEeigI0eOqHLlKrp48YKrS7mhyympri4BAAqUi2dPSbLp0UcfdXUp1+XjU1R79uwm1DkJgQ5KTEzUxYsXVK/HMAWERbm6nCwd37FeO5dO15UrV1xdCgAUKJcvnJNkVPuRF1UyurKry8lS0vFD2vj+CCUmJhLonKTAB7p169bpzTff1JYtW3T8+HEtXrxYbdu2tc83xmjYsGF69913debMGTVs2FBTp05VhQoVXFe0RQWERal4RCVXl5GlpOOHXF0CABRofqUiCuwxHM5X4L8Ucf78edWqVUuTJ0/Ocv4bb7yht956S9OmTdPGjRvl6+uruLg4Xbp0KZ8rBQAAcI0Cf4auVatWatWqVZbzjDGaMGGCXn75ZbVp00aS9MEHHygkJERLlixR586d87NUAAAAlyjwZ+iu5+DBgzpx4oRatmxpnxYYGKh69epp/fr1LqwMAAAg/xT4M3TXc+LECUlSSEiIw/SQkBD7vKykpKQoJSXF/jgpKck5BQIAAOQDS5+hu1mjRo1SYGCg/ads2bKuLgkAAOCmWTrQhYaGSpISEhIcpickJNjnZWXw4ME6e/as/efo0aNOrRMAAMCZLB3ooqOjFRoaqpUrV9qnJSUlaePGjYqNjc32eV5eXgoICHD4AQAAsKoCfw1dcnKy9u/fb3988OBBbdu2TcWLF1dERIT69eunV199VRUqVFB0dLSGDh2q8PBwh3vVAQAA3MoKfKDbvHmz7rrrLvvj/v37S5Li4+M1a9YsDRw4UOfPn9fTTz+tM2fOqFGjRlq2bJm8vb1dVTIAAEC+KvCBrlmzZjLGZDvfZrNp5MiRGjlyZD5WBQAAUHBY+ho6AAAAEOgAAAAsj0AHAABgcQQ6AAAAiyPQAQAAWByBDgAAwOIIdAAAABZHoAMAALA4Ah0AAIDFEegAAAAsjkAHAABgcQQ6AAAAiyPQAQAAWJy7qwsoDI4cOaLExERXl5Gt3bt3u7oEAADwDxDonOzIkSOqXLmKLl684OpSbuhySqqrSwAAADeBQOdkiYmJunjxgur1GKaAsChXl5Ol4zvWa+fS6bpy5YqrSwEAADeBQJdPAsKiVDyikqvLyFLS8UOuLgEAAPwDfCkCAADA4gh0AAAAFkegAwAAsDgCHQAAgMUR6AAAACyOQAcAAGBxBDoAAACLI9ABAABYHIEOAADA4gh0AAAAFkegAwAAsDgCHQAAgMUR6AAAACyOQAcAAGBxBDoAAACLI9ABAABYHIEOAADA4gh0AAAAFkegAwAAsDgCHQAAgMUR6AAAACyOQAcAAGBxBDoAAACLI9ABAABYHIEOAADA4gh0AAAAFkegAwAAsDgCHQAAgMUR6AAAACyOQAcAAGBxBDoAAACLI9ABAABYHIEOAADA4gh0AAAAFkegAwAAsDgCHQAAgMUR6AAAACyOQAcAAGBxBDoAAACLI9ABAABYHIEOAADA4m6pQDd58mRFRUXJ29tb9erV06ZNm1xdEgAAgNPdMoFuwYIF6t+/v4YNG6atW7eqVq1aiouL08mTJ11dGgAAgFPdMoFu3Lhxeuqpp9S9e3dVrVpV06ZNU9GiRfX++++7ujQAAACnuiUCXWpqqrZs2aKWLVvap7m5ually5Zav369CysDAABwPndXF5AXEhMTlZaWppCQEIfpISEh2rNnT6bxKSkpSklJsT8+e/asJCkpKSnPa0tOTpYk/XV4r66kXMzz5eeFpOOHJUln/9gnD3ebi6vJGjXmDWrMG9SYN6gxb1iixhNHJP39OzEvf9dmLMsYk2fLtCqbuQW2wrFjx1S6dGn98MMPio2NtU8fOHCg1q5dq40bNzqMHz58uEaMGJHfZQIAACc4evSoypQp4+oyXOqWOEMXHBysIkWKKCEhwWF6QkKCQkNDM40fPHiw+vfvb3+cnp6uv/76SyVKlJDNlvl/N0lJSSpbtqyOHj2qgICAvG+gACgMPUqFo096vHUUhj7p8dbhij6NMTp37pzCw8PzZX0F2S0R6Dw9PVW3bl2tXLlSbdu2lfR3SFu5cqV69+6dabyXl5e8vLwcpgUFBd1wPQEBAbf0m1EqHD1KhaNPerx1FIY+6fHWkd99BgYG5tu6CrJbItBJUv/+/RUfH6/bb79dd955pyZMmKDz58+re/furi4NAADAqW6ZQNepUyf9+eef+s9//qMTJ06odu3aWrZsWaYvSgAAANxqbplAJ0m9e/fO8iPWf8rLy0vDhg3L9DHtraQw9CgVjj7p8dZRGPqkx1tHYemzoLolvuUKAABQmN0SNxYGAAAozAh0AAAAFkegAwAAsDgCHQAAgMUVykA3efJkRUVFydvbW/Xq1dOmTZuyHXv58mWNHDlSMTEx8vb2Vq1atbRs2TKHMevWrdMDDzyg8PBw2Ww2LVmyxMkd5Exe9zlq1Cjdcccd8vf3V6lSpdS2bVvt3bvX2W1cV173OHXqVNWsWdN+Y8zY2Fh99dVXzm7juvK6x6v997//lc1mU79+/ZxQee7kdZ/Dhw+XzWZz+KlcubKz27guZ7yWf/zxhx599FGVKFFCPj4+qlGjhjZv3uzMNq4rr3uMiorK9DrabDb16tXL2a1cV173mZaWpqFDhyo6Olo+Pj6KiYnRK6+84tK/UZrXPZ47d079+vVTZGSkfHx81KBBA/3444/ObqPwMIXM/Pnzjaenp3n//ffNrl27zFNPPWWCgoJMQkJCluMHDhxowsPDzRdffGEOHDhgpkyZYry9vc3WrVvtY7788kszZMgQs2jRIiPJLF68OJ+6yZ4z+oyLizMzZ840O3fuNNu2bTP33XefiYiIMMnJyfnVlgNn9Lh06VLzxRdfmF9//dXs3bvXvPTSS8bDw8Ps3Lkzv9py4IweM2zatMlERUWZmjVrmr59+zq5k+tzRp/Dhg0z1apVM8ePH7f//Pnnn/nVUibO6PGvv/4ykZGRplu3bmbjxo3mt99+M8uXLzf79+/Pr7YcOKPHkydPOryGX3/9tZFkVq9enU9dZeaMPl977TVTokQJ8/nnn5uDBw+ahQsXGj8/PzNx4sT8asuBM3rs2LGjqVq1qlm7dq3Zt2+fGTZsmAkICDC///57frV1Syt0ge7OO+80vXr1sj9OS0sz4eHhZtSoUVmODwsLM5MmTXKY1q5dO9O1a9csxxeUQOfsPo35+0AryaxduzZvis6l/OjRGGOKFStm3nvvvX9e8E1wVo/nzp0zFSpUMF9//bVp2rSpywOdM/ocNmyYqVWrllPqvRnO6PHFF180jRo1ck7BNyE/3pN9+/Y1MTExJj09PW+KvgnO6LN169amR48e1x2Tn/K6xwsXLpgiRYqYzz//3GFMnTp1zJAhQ/K4+sKpUH3kmpqaqi1btqhly5b2aW5ubmrZsqXWr1+f5XNSUlLk7e3tMM3Hx0ffffedU2v9J/Krz7Nnz0qSihcvngdV505+9JiWlqb58+fr/Pnzio2Nzbvic8iZPfbq1UutW7d2WLarOLPPffv2KTw8XOXKlVPXrl115MiRvG8gB5zV49KlS3X77berQ4cOKlWqlG677Ta9++67zmniBvLjPZmamqoPP/xQPXr0kM1my7vic8FZfTZo0EArV67Ur7/+Kknavn27vvvuO7Vq1coJXVyfM3q8cuWK0tLSLPf71EoKVaBLTExUWlpapj8HFhISohMnTmT5nLi4OI0bN0779u1Tenq6vv76ay1atEjHjx/Pj5JvSn70mZ6ern79+qlhw4aqXr16nvdwI87scceOHfLz85OXl5f+9a9/afHixapatarTesmOs3qcP3++tm7dqlGjRjm1/pxyVp/16tXTrFmztGzZMk2dOlUHDx5U48aNde7cOaf2kxVn9fjbb79p6tSpqlChgpYvX66ePXuqT58+mj17tlP7yUp+HHeWLFmiM2fOqFu3bnldfo45q89Bgwapc+fOqly5sjw8PHTbbbepX79+6tq1q1P7yYozevT391dsbKxeeeUVHTt2TGlpafrwww+1fv36Av371EoKVaC7GRMnTlSFChVUuXJleXp6qnfv3urevbvc3G6tTZfbPnv16qWdO3dq/vz5+Vzpzctpj5UqVdK2bdu0ceNG9ezZU/Hx8frll19cVHXu3KjHo0ePqm/fvpo7d26m/ylbSU5ey1atWqlDhw6qWbOm4uLi9OWXX+rMmTP66KOPXFh5zuWkx/T0dNWpU0evv/66brvtNj399NN66qmnNG3aNBdWnnO5Pe7MmDFDrVq1Unh4eD5X+s/kpM+PPvpIc+fO1bx587R161bNnj1bY8aMcUk4vxk56XHOnDkyxqh06dLy8vLSW2+9pS5dutxyv09dpVBtxeDgYBUpUkQJCQkO0xMSEhQaGprlc0qWLKklS5bo/PnzOnz4sPbs2SM/Pz+VK1cuP0q+Kc7us3fv3vr888+1evVqlSlTxik93Igze/T09FT58uVVt25djRo1SrVq1dLEiROd1kt2nNHjli1bdPLkSdWpU0fu7u5yd3fX2rVr9dZbb8nd3V1paWlO7+ta+fW+DAoKUsWKFbV///48rT8nnNVjWFhYprPHVapUcclHy85+HQ8fPqxvvvlGTz75pFPqzyln9fnCCy/Yz9LVqFFDjz32mJ577jmXnEl3Vo8xMTFau3atkpOTdfToUW3atEmXL18u0L9PraRQBTpPT0/VrVtXK1eutE9LT0/XypUrb3iNlLe3t0qXLq0rV67ok08+UZs2bZxd7k1zVp/GGPXu3VuLFy/WqlWrFB0d7bQebiQ/X8v09HSlpKTkSd254YweW7RooR07dmjbtm32n9tvv11du3bVtm3bVKRIEaf2lJX8ei2Tk5N14MABhYWF5VntOeWsHhs2bJjp1kG//vqrIiMj87aBHHD26zhz5kyVKlVKrVu3zvPac8NZfV64cCHTmaoiRYooPT09bxvIAWe/lr6+vgoLC9Pp06e1fPnyAv371FJc/a2M/DZ//nzj5eVlZs2aZX755Rfz9NNPm6CgIHPixAljjDGPPfaYGTRokH38hg0bzCeffGIOHDhg1q1bZ5o3b26io6PN6dOn7WPOnTtnfvrpJ/PTTz8ZSWbcuHHmp59+MocPH87v9uyc0WfPnj1NYGCgWbNmjcNtBC5cuJDf7RljnNPjoEGDzNq1a83BgwfNzz//bAYNGmRsNptZsWJFfrdnjHFOj9cqCN9ydUafAwYMMGvWrDEHDx4033//vWnZsqUJDg42J0+ezO/2jDHO6XHTpk3G3d3dvPbaa2bfvn1m7ty5pmjRoubDDz/M7/aMMc7bX9PS0kxERIR58cUX87OdbDmjz/j4eFO6dGn7bUsWLVpkgoODzcCBA/O7PWOMc3pctmyZ+eqrr8xvv/1mVqxYYWrVqmXq1atnUlNT87u9W1KhC3TGGPP222+biIgI4+npae68806zYcMG+7ymTZua+Ph4++M1a9aYKlWqGC8vL1OiRAnz2GOPmT/++MNheatXrzaSMv1cvRxXyOs+s+pRkpk5c2Y+dZRZXvfYo0cPExkZaTw9PU3JkiVNixYtXBbmMuR1j9cqCIHOmLzvs1OnTiYsLMx4enqa0qVLm06dOrns/mwZnPFafvbZZ6Z69erGy8vLVK5c2UyfPj0/WsmWM3pcvny5kWT27t2bHy3kSF73mZSUZPr27WsiIiKMt7e3KVeunBkyZIhJSUnJr5YyyeseFyxYYMqVK2c8PT1NaGio6dWrlzlz5kx+tXPLsxnjwttQAwAA4B8rVNfQAQAA3IoIdAAAABZHoAMAALA4Ah0AAIDFEegAAAAsjkAHAABgcQQ6AAAAiyPQAYCLDB8+XLVr1y4wywFgXQQ6wKL+/PNP9ezZUxEREfLy8lJoaKji4uL0/fffu7q0AuuTTz5Rs2bNFBgYKD8/P9WsWVMjR47UX3/95erScsxms2nJkiUO055//nmHv7sJoPAh0AEW1b59e/3000+aPXu2fv31Vy1dulTNmjXTqVOnXF2aLl++7OoSMhkyZIg6deqkO+64Q1999ZV27typsWPHavv27ZozZ85NLzerXlNTU/9Jqbnm5+enEiVK5Os6ARQwrv7bYwBy7/Tp00aSWbNmzXXHSTJTpkwx9957r/H29jbR0dFm4cKFDmOOHDliOnToYAIDA02xYsXMgw8+aA4ePGifv2nTJtOyZUtTokQJExAQYJo0aWK2bNmS5XoeeOABU7RoUTNs2DAzbNgwU6tWLTNjxgxTtmxZ4+vra3r27GmuXLliRo8ebUJCQkzJkiXNq6++6rCssWPHmurVq5uiRYuaMmXKmJ49e5pz587Z58+cOdMEBgaaZcuWmcqVKxtfX18TFxdnjh07lu122Lhxo5FkJkyYkO32zDBlyhRTrlw54+HhYSpWrGg++OCDHPf67rvvmqioKGOz2ezLfeKJJ0xwcLDx9/c3d911l9m2bZt9WRnPy+m2joyMdPg7ypGRkVkuJy0tzYwYMcKULl3aeHp6mlq1apmvvvrKPv/gwYNGkvnkk09Ms2bNjI+Pj6lZs6b54Ycfst2GAAo2ztABFuTn5yc/Pz8tWbJEKSkp1x07dOhQtW/fXtu3b1fXrl3VuXNn7d69W9LfZ5fi4uLk7++vb7/9Vt9//738/Px077332s8ynTt3TvHx8fruu++0YcMGVahQQffdd5/OnTvnsJ7hw4froYce0o4dO9SjRw9J0oEDB/TVV19p2bJl+t///qcZM2aodevW+v3337V27VqNHj1aL7/8sjZu3Ghfjpubm9566y3t2rVLs2fP1qpVqzRw4ECHdV24cEFjxozRnDlztG7dOh05ckTPP/98tttg7ty58vPz07///e8s5wcFBUmSFi9erL59+2rAgAHauXOnnnnmGXXv3l2rV6++Ya/79+/XJ598okWLFmnbtm2SpA4dOujkyZP66quvtGXLFtWpU0ctWrTI9iPeG23rH3/8UZI0c+ZMHT9+3P74WhMnTtTYsWM1ZswY/fzzz4qLi9ODDz6offv2OYwbMmSInn/+eW3btk0VK1ZUly5ddOXKlWy3I4ACzNWJEsDN+fjjj02xYsWMt7e3adCggRk8eLDZvn27wxhJ5l//+pfDtHr16pmePXsaY4yZM2eOqVSpkklPT7fPT0lJMT4+Pmb58uVZrjctLc34+/ubzz77zGE9/fr1cxg3bNgwU7RoUZOUlGSfFhcXZ6KiokxaWpp9WqVKlcyoUaOy7XPhwoWmRIkS9sczZ840ksz+/fvt0yZPnmxCQkKyXUarVq1MzZo1s52foUGDBuapp55ymNahQwdz33332R9n16uHh4c5efKkfdq3335rAgICzKVLlxzGxsTEmHfeecf+vKvPrF0ru229ePHiTOu/ejnh4eHmtddecxhzxx13mH//+9/GmP87Q/fee+/Z5+/atctIMrt37862HgAFF2foAItq3769jh07pqVLl+ree+/VmjVrVKdOHc2aNcthXGxsbKbHGWfotm/frv3798vf399+1q948eK6dOmSDhw4IElKSEjQU089pQoVKigwMFABAQFKTk7WkSNHHJZ7++23Z6oxKipK/v7+9schISGqWrWq3NzcHKadPHnS/vibb75RixYtVLp0afn7++uxxx7TqVOndOHCBfuYokWLKiYmxv44LCzMYRnXMsZkO+9qu3fvVsOGDR2mNWzY0L69MmTVa2RkpEqWLGl/vH37diUnJ6tEiRL2bevn56eDBw/at+21crqtrycpKUnHjh3LUR81a9a0/zssLEySrrsdARRc7q4uAMDN8/b21t133627775bQ4cO1ZNPPqlhw4apW7duOXp+cnKy6tatq7lz52aalxFO4uPjderUKU2cOFGRkZHy8vJSbGxspgv/fX19My3Dw8PD4bHNZstyWnp6uiTp0KFDuv/++9WzZ0+99tprKl68uL777js98cQTSk1NVdGiRbNd7vVCW8WKFfXdd9/p8uXLmZ57M7Lq9dppycnJCgsL05o1azKNzfiI91o53dZ55eptYbPZJMn+WgCwFs7QAbeQqlWr6vz58w7TNmzYkOlxlSpVJEl16tTRvn37VKpUKZUvX97hJzAwUJL0/fffq0+fPrrvvvtUrVo1eXl5KTEx0Sn1b9myRenp6Ro7dqzq16+vihUr6tixY/94uY888oiSk5M1ZcqULOefOXNGklSlSpVMt335/vvvVbVq1Vyvs06dOjpx4oTc3d0zbdvg4OAsn5OTbe3h4aG0tLRs1xsQEKDw8PA86wOANRDoAAs6deqUmjdvrg8//FA///yzDh48qIULF+qNN95QmzZtHMYuXLhQ77//vn799VcNGzZMmzZtUu/evSVJXbt2VXBwsNq0aaNvv/1WBw8e1Jo1a9SnTx/9/vvvkqQKFSpozpw52r17tzZu3KiuXbvKx8fHKX2VL19ely9f1ttvv63ffvtNc+bM0bRp0/7xcuvVq6eBAwdqwIABGjhwoNavX6/Dhw9r5cqV6tChg2bPni1JeuGFFzRr1ixNnTpV+/bt07hx47Ro0aLrfuEiOy1btlRsbKzatm2rFStW6NChQ/rhhx80ZMgQbd68Ocvn5GRbR0VFaeXKlTpx4oROnz6d5XJeeOEFjR49WgsWLNDevXs1aNAgbdu2TX379s11HwCsgUAHWJCfn5/q1aun8ePHq0mTJqpevbqGDh2qp556SpMmTXIYO2LECM2fP181a9bUBx98oP/973/2MzVFixbVunXrFBERoXbt2qlKlSp64okndOnSJQUEBEiSZsyYodOnT6tOnTp67LHH1KdPH5UqVcopfdWqVUvjxo3T6NGjVb16dc2dO1ejRo3Kk2WPHj1a8+bN08aNGxUXF6dq1aqpf//+qlmzpuLj4yVJbdu21cSJEzVmzBhVq1ZN77zzjmbOnKlmzZrlen02m01ffvmlmjRpou7du6tixYrq3LmzDh8+rJCQkCyfk5NtPXbsWH399dcqW7asbrvttiyX06dPH/Xv318DBgxQjRo1tGzZMi1dulQVKlTIdR8ArMFmcnq1MADLsdlsWrx4sdq2bevqUgAATsQZOgAAAIsj0AEAAFgcty0BbmFcUQEAhQNn6AAAACyOQAcAAGBxBDoAAACLI9ABAABYHIEOAADA4gh0AAAAFkegAwAAsDgCHQAAgMUR6AAAACzu/wG05JCpt/cppAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot train correlation distribution vs test correlation distribution\n",
    "test_corr = np.array(test_corr_singles + test_corr_combos)[:, 1].astype(float)\n",
    "sns.histplot(test_corr, label=\"test_corr\")\n",
    "plt.title(\"Correlation between predicted and observed expression per perturbation\")\n",
    "plt.xlabel(\"Spearman Correlation\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "2336d296-2911-4607-a902-552a7e7b92fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(\n",
    "    \"figures/reconstruction_spearmans/pertspectra_inhouse_spearman_correlations.csv\",\n",
    "    test_corr,\n",
    "    delimiter=\",\",\n",
    "    fmt=\"%.2f\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64fc0b35-f87e-4d8a-87f1-e6080a3ee83f",
   "metadata": {},
   "source": [
    "# Perturbation Embedding Heatmap"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c36ebe2b-2835-44ed-9f71-a52b522125b1",
   "metadata": {},
   "source": [
    "## Hierarchical Clustering + Enrichment\n",
    "- Perform enrichment tests (gprofiler) on hierarchical clustering of perturbation embeddings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0aaf8dc3-7e5b-4748-8685-fc8560bf7fc7",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.spatial.distance import pdist, squareform\n",
    "\n",
    "# Compute the pairwise distances\n",
    "tnfa_pos_pert_embeddings_df = pd.DataFrame(\n",
    "    adata.uns[\"SPECTRA_pert_scores\"], index=adata.uns[\"Spectra_pert_labels\"]\n",
    ")\n",
    "df = tnfa_pos_pert_embeddings_df.drop(index=[\"basal\", \"ctrl\", \"intergenic\"])\n",
    "\n",
    "distance = \"euclidean\"\n",
    "distances = pdist(df.values, metric=distance)\n",
    "\n",
    "# Convert the distances into a square distance matrix\n",
    "distance_matrix = pd.DataFrame(squareform(distances), index=df.index, columns=df.index)\n",
    "clustermap = sns.clustermap(distance_matrix, cmap=\"viridis_r\")\n",
    "clustermap.fig.suptitle(f\"Pairwise {distance} distance of perturbation latent vectors\")\n",
    "# clustermap.ax_row_dendrogram.set_visible(False)\n",
    "clustermap.ax_col_dendrogram.set_visible(False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "b2181a72-86cf-47a0-aa04-ab154fe0bf95",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjIAAAHLCAYAAADSuXIVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABpSUlEQVR4nO3deXhMZ/8G8O9kX0SESCLWIIgtqUYilloaiaX2XWtJLUWV1k5DKBXqpbS8VKutqv1VtIoWXah9J/Y1CZoQJSGJrPfvD785nUkmZJJJ4sT9ua5c1ZkzZ545c+ac+zzb0QCAEBEREamQWVEXgIiIiCivGGSIiIhItRhkiIiISLUYZIiIiEi1GGSIiIhItRhkiIiISLUYZIiIiEi1LIq6AKaQmZkpd+7cEQcHB9FoNEVdHCIiIsoFAPLo0SNxd3cXM7O81a0UiyBz584dqVixYlEXg4iIiPIgOjpaKlSokKfXFosg4+DgICJPN0TJkiWLuDRERESUGwkJCVKxYkXlPJ4XxSLIaJuTSpYsySBDRESkMvnpFsLOvkRERKRaDDJERESkWgwyREREpFoMMkRERKRaDDJERESkWgwyREREpFoMMkRERKRaDDJERESkWgwyREREpFoMMkRERKRaDDJERESkWgwyREREpFoMMkRERKRaxeLu10UFgCSnZRR1MYhIh62leb7upEtE6sIgk0cApPuyg3I88kFRF4WIdPhWdpKNwwIYZoheEmxayqPktAyGGKIX0LHIB6wpJXqJsEbGBI6FBoqdlXlRF4PopZaUmiG+s3YXdTGIqJAxyJiAnZW52FlxUxIRERU2Ni0RERGRajHIEBERkWoxyBAREZFqMcgQERGRajHIEBERkWoxyBAREZFqMcgQERGRajHIEBERkWoxyBAREZFqMcgQERGRajHIEBERkWoxyBAREZFqMcgQERGRajHIEBERkWoxyBAREZFqMcgQERGRajHIEBERkWoxyBAREZFqMcgQERGRajHIEBERkWoxyBAREZFqMcgQERGRajHIEBERkWoxyBAREZFqMcgQERGRajHIEBERkWoxyBAREZFqMcgQERGRajHIEBERkWoxyBAREZFqMcgQERGRajHIEBERkWoxyBAREZFqMcgQERGRajHIEBERkWoxyBAREZFqMcgQERGRajHIEBERkWoxyBAREZFqMcgQERGRajHIEBERkWrlKcgsWbJEqlSpIjY2NuLv7y9Hjhx55vIbN26UWrVqiY2NjdSrV0+2b9+e47LDhg0TjUYjCxcuzEvRiIiI6CVidJBZv369jBkzRsLCwuTEiRPi7e0twcHBcvfuXYPLHzhwQPr06SODBg2SkydPSufOnaVz584SERGRbdnNmzfLoUOHxN3d3fhPQkRERC8do4PMggULZMiQIRISEiK1a9eWZcuWiZ2dnXz99dcGl1+0aJG0adNGxo8fL15eXjJz5kxp0KCBLF68WG+527dvy3vvvSerV68WS0vLZ5YhJSVFEhIS9P6IiIjo5WNUkElNTZXjx49LYGDgvyswM5PAwEA5ePCgwdccPHhQb3kRkeDgYL3lMzMzpV+/fjJ+/HipU6fOc8sRHh4ujo6Oyl/FihWN+RhERERUTBgVZOLi4iQjI0NcXV31Hnd1dZWYmBiDr4mJiXnu8nPnzhULCwsZNWpUrsoxefJkiY+PV/6io6ON+RhERERUTFgUdQGOHz8uixYtkhMnTohGo8nVa6ytrcXa2rqAS0ZEREQvOqNqZJydncXc3FxiY2P1Ho+NjRU3NzeDr3Fzc3vm8vv27ZO7d+9KpUqVxMLCQiwsLCQyMlLGjh0rVapUMaZ4RERE9JIxKshYWVnJq6++Knv27FEey8zMlD179khAQIDB1wQEBOgtLyKya9cuZfl+/frJmTNn5NSpU8qfu7u7jB8/Xn755RdjPw8RERG9RIxuWhozZowMGDBAfH19xc/PTxYuXCiJiYkSEhIiIiL9+/eX8uXLS3h4uIiIjB49Wpo3by7z58+X9u3by7p16+TYsWOyfPlyEREpU6aMlClTRu89LC0txc3NTWrWrJnfz0dERETFmNFBplevXnLv3j2ZNm2axMTEiI+Pj+zcuVPp0BsVFSVmZv9W9DRu3FjWrFkjoaGhMmXKFPH09JQtW7ZI3bp1TfcpiIiI6KWkAYCiLkR+JSQkiKOjo8THx0vJkiUL5T2TUtOl9rSnTV/nPwoWO6si7zdN9FLjb5JIfUxx/ua9loiIiEi1GGSIiIhItRhkiIiISLUYZIiIiEi1GGSIiIhItRhkiIiISLUYZIiIiEi1GGSIiIhItRhkiIiISLUYZIiIiEi1GGSIiIhItRhkiIiISLUYZIiIiEi1GGSIiIhItRhkiIiISLUYZIiIiEi1GGSIiIhItRhkiIiISLUYZIiIiEi1GGSIiIhItRhkiIiISLUYZIiIiEi1GGSIiIhItRhkiIiISLUYZIiIiEi1GGSIiIhItRhkiIiISLUYZIiIiEi1GGSIiIhItRhkiIiISLUYZIiIiEi1GGSIiIhItRhkiIiISLUYZIiIiEi1GGSIiIhItRhkiIiISLUYZIiIiEi1GGSIiIhItRhkiIiISLUYZIiIiEi1GGSIiIhItRhkiIiISLUYZIiIiEi1GGSIiIhItRhkiIiISLUYZIiIiEi1GGSIiIhItSyKugBELyoAkpyWUdTFoFxKSk03+G9SB1tLc9FoNEVdDFIhBhkiAwBI92UH5Xjkg6IuCuWB76w9RV0EMpJvZSfZOCyAYYaMxqYlIgOS0zIYYogK0bHIB6wBpTxhjQzRcxwLDRQ7K/OiLgZRsZSUmiG+s3YXdTFIxRhkiJ7Dzspc7Kz4UyEiehGxaYmIiIhUi0GGiIiIVItBhoiIiFSLQYaIiIhUi0GGiIiIVItBhoiIiFSLQYaIiIhUi0GGiIiIVItBhoiIiFSLQYaIiIhUi0GGiIiIVItBhoiIiFSLQYaIiIhUK09BZsmSJVKlShWxsbERf39/OXLkyDOX37hxo9SqVUtsbGykXr16sn37dr3np0+fLrVq1RJ7e3txcnKSwMBAOXz4cF6KRkRERC8Ro4PM+vXrZcyYMRIWFiYnTpwQb29vCQ4Olrt37xpc/sCBA9KnTx8ZNGiQnDx5Ujp37iydO3eWiIgIZZkaNWrI4sWL5ezZs/LXX39JlSpVJCgoSO7du5f3T0ZERETFntFBZsGCBTJkyBAJCQmR2rVry7Jly8TOzk6+/vprg8svWrRI2rRpI+PHjxcvLy+ZOXOmNGjQQBYvXqws07dvXwkMDJSqVatKnTp1ZMGCBZKQkCBnzpzJ+ycjIiKiYs+oIJOamirHjx+XwMDAf1dgZiaBgYFy8OBBg685ePCg3vIiIsHBwTkun5qaKsuXLxdHR0fx9vY2uExKSookJCTo/REREdHLx6ggExcXJxkZGeLq6qr3uKurq8TExBh8TUxMTK6W37Ztm5QoUUJsbGzk008/lV27domzs7PBdYaHh4ujo6PyV7FiRWM+BhERERUTL8yopZYtW8qpU6fkwIED0qZNG+nZs2eO/W4mT54s8fHxyl90dHQhl5aIiIheBEYFGWdnZzE3N5fY2Fi9x2NjY8XNzc3ga9zc3HK1vL29vVSvXl0aNWokK1asEAsLC1mxYoXBdVpbW0vJkiX1/oiIiOjlY1SQsbKykldffVX27NmjPJaZmSl79uyRgIAAg68JCAjQW15EZNeuXTkur7velJQUY4pHRERELxkLY18wZswYGTBggPj6+oqfn58sXLhQEhMTJSQkRERE+vfvL+XLl5fw8HARERk9erQ0b95c5s+fL+3bt5d169bJsWPHZPny5SIikpiYKB9//LF07NhRypUrJ3FxcbJkyRK5ffu29OjRw4QflYiIiIobo4NMr1695N69ezJt2jSJiYkRHx8f2blzp9KhNyoqSszM/q3oady4saxZs0ZCQ0NlypQp4unpKVu2bJG6deuKiIi5ublcvHhRVq5cKXFxcVKmTBlp2LCh7Nu3T+rUqWOij0lERETFkdFBRkRk5MiRMnLkSIPP/fHHH9ke69GjR461KzY2NvLDDz/kpRhERET0knthRi0RERERGYtBhoiIiFSLQYaIiIhUi0GGiIiIVItBhoiIiFSLQYaIiIhUi0GGiIiIVItBhoiIiFSLQYaIiIhUi0GGiIiIVItBhoiIiFSLQYaIiIhUi0GGiIiIVItBhoiIiFSLQYaIiIhUi0GGiIiIVItBhoiIiFSLQYaIiIhUi0GGiIiIVItBhoiIiFSLQYaIiIhUi0GGiIiIVItBhoiIiFSLQYaIiIhUi0GGiIiIVItBhoiIiFSLQYaIiIhUi0GGiIiIVItBhoiIiFSLQYaIiIhUi0GGiIiIVItBhoiIiFSLQYaIiIhUi0GGiIiIVItBhoiIiFSLQYaIiIhUi0GGiIiIVItBhoiIiFSLQYaIiIhUi0GGiIiIVItBhoiIiFSLQYaIiIhUi0GGiIiIVItBhoiIiFSLQYaIiIhUi0GGiIiIVItBhoiIiFSLQYaIiIhUi0GGiIiIVItBhoiIiFSLQYaIiIhUi0GGiIiIVItBhoiIiFSLQYaIiIhUi0GGiIiIVItBhoiIiFSLQYaIiIhUi0GGiIiIVItBhoiIiFSLQYaIiIhUi0GGiIiIVItBhoiIiFSLQYaIiIhUi0GGiIiIVItBhoiIiFQrT0FmyZIlUqVKFbGxsRF/f385cuTIM5ffuHGj1KpVS2xsbKRevXqyfft25bm0tDSZOHGi1KtXT+zt7cXd3V369+8vd+7cyUvRiIiI6CVidJBZv369jBkzRsLCwuTEiRPi7e0twcHBcvfuXYPLHzhwQPr06SODBg2SkydPSufOnaVz584SEREhIiJJSUly4sQJmTp1qpw4cUJ++OEHuXTpknTs2DF/n4yIiIiKPaODzIIFC2TIkCESEhIitWvXlmXLlomdnZ18/fXXBpdftGiRtGnTRsaPHy9eXl4yc+ZMadCggSxevFhERBwdHWXXrl3Ss2dPqVmzpjRq1EgWL14sx48fl6ioqPx9OiIiIirWjAoyqampcvz4cQkMDPx3BWZmEhgYKAcPHjT4moMHD+otLyISHByc4/IiIvHx8aLRaKRUqVIGn09JSZGEhAS9PyIiInr5GBVk4uLiJCMjQ1xdXfUed3V1lZiYGIOviYmJMWr5J0+eyMSJE6VPnz5SsmRJg8uEh4eLo6Oj8lexYkVjPgYREREVEy/UqKW0tDTp2bOnAJClS5fmuNzkyZMlPj5e+YuOji7EUhIREdGLwsKYhZ2dncXc3FxiY2P1Ho+NjRU3NzeDr3Fzc8vV8toQExkZKb/99luOtTEiItbW1mJtbW1M0YmIiKgYMqpGxsrKSl599VXZs2eP8lhmZqbs2bNHAgICDL4mICBAb3kRkV27duktrw0xV65ckd27d0uZMmWMKRYRERG9pIyqkRERGTNmjAwYMEB8fX3Fz89PFi5cKImJiRISEiIiIv3795fy5ctLeHi4iIiMHj1amjdvLvPnz5f27dvLunXr5NixY7J8+XIReRpiunfvLidOnJBt27ZJRkaG0n+mdOnSYmVlZarPSkRERMWM0UGmV69ecu/ePZk2bZrExMSIj4+P7Ny5U+nQGxUVJWZm/1b0NG7cWNasWSOhoaEyZcoU8fT0lC1btkjdunVFROT27dvy448/ioiIj4+P3nv9/vvv0qJFizx+NCIiIirujA4yIiIjR46UkSNHGnzujz/+yPZYjx49pEePHgaXr1KligDISzGIiIjoJfdCjVoiIiIiMgaDDBEREakWgwwRERGpFoMMERERqRaDDBEREakWgwwRERGpFoMMERERqRaDDBEREakWgwwRERGpFoMMERERqRaDDBEREakWgwwRERGpFoMMERERqRaDDBEREakWgwwRERGpFoMMERERqRaDDBEREakWgwwRERGpFoMMERERqZZFUReAiEjtAEhyWkZRF0OVklLTDf6bjGNraS4ajaaoi1EkXsogY4qDjql/fC/zTkikZgCk+7KDcjzyQVEXRfV8Z+0p6iKolm9lJ9k4LOClPI+8dEGmIA46pvjxvcw7IZGaJadlMMRQkTsW+UCS0zLEzuqlO62/fEHmRT3ovMw7IVFxcSw0UOyszIu6GPQSSUrNEN9Zu4u6GEXqpT5rvggHHe6ERMWHnZU5L0aICtlL/YvjQYeIiEjdOPyaiIiIVItBhoiIiFSLQYaIiIhUi0GGiIiIVItBhoiIiFSLQYaIiIhUi0GGiIiIVItBhoiIiFSLQYaIiIhUi0GGiIiIVItBhoiIiFSLQYaIiIhUi0GGiIiIVItBhoiIiFSLQYaIiIhUi0GGiIiIVItBhoiIiFSLQYaIiIhUi0GGiIiIVItBhoiIiFSLQYaIiIhUi0GGiIiIVItBhoiIiFSLQYaIiIhUi0GGiIiIVItBhoiIiFSLQYaIiIhUy6KoC0DFHwBJTsso6mIYJSk13eC/1cDW0lw0Gk1RF4OIqFAwyFCBAiDdlx2U45EPirooeeY7a09RF8EovpWdZOOwAIYZInopsGmJClRyWoaqQ4waHYt8oLoaMCKivGKNDBWaY6GBYmdlXtTFKLaSUjPEd9buoi4GEVGhYpChQmNnZS52VtzliIjIdNi0RERERKrFIENERESqxSBDREREqsUgQ0RERKrFIENERESqxSBDREREqsUgQ0RERKrFST2IiIgKWEHdc66w7gv3It/DLU9BZsmSJTJv3jyJiYkRb29v+fzzz8XPzy/H5Tdu3ChTp06Vmzdviqenp8ydO1fatWunPP/DDz/IsmXL5Pjx4/LPP//IyZMnxcfHJy9FIyIieqEU1j3nCvK+cC/yPdyMblpav369jBkzRsLCwuTEiRPi7e0twcHBcvfuXYPLHzhwQPr06SODBg2SkydPSufOnaVz584SERGhLJOYmChNmzaVuXPn5v2TEBERvYCKwz3nXuR7uBldI7NgwQIZMmSIhISEiIjIsmXL5Oeff5avv/5aJk2alG35RYsWSZs2bWT8+PEiIjJz5kzZtWuXLF68WJYtWyYiIv369RMRkZs3b+b1cxAREb3w1HbPOTXcw82oIJOamirHjx+XyZMnK4+ZmZlJYGCgHDx40OBrDh48KGPGjNF7LDg4WLZs2WJ8af9fSkqKpKSkKP+fkJCQ53UREREVFt5zzvSMalqKi4uTjIwMcXV11Xvc1dVVYmJiDL4mJibGqOVzIzw8XBwdHZW/ihUr5nldREREpF6qHH49efJkiY+PV/6io6OLukhERERUBIyq33J2dhZzc3OJjY3Vezw2Nlbc3NwMvsbNzc2o5XPD2tparK2t8/x6IiIiKh6MqpGxsrKSV199Vfbs+XeIV2ZmpuzZs0cCAgIMviYgIEBveRGRXbt25bg8ERERUW4Z3eNozJgxMmDAAPH19RU/Pz9ZuHChJCYmKqOY+vfvL+XLl5fw8HARERk9erQ0b95c5s+fL+3bt5d169bJsWPHZPny5co6//nnH4mKipI7d+6IiMilS5dE5GltTn5qboiIiKh4MzrI9OrVS+7duyfTpk2TmJgY8fHxkZ07dyodeqOiosTM7N+KnsaNG8uaNWskNDRUpkyZIp6enrJlyxapW7eussyPP/6oBCERkd69e4uISFhYmEyfPj2vn42IiIiKuTyNARs5cqSMHDnS4HN//PFHtsd69OghPXr0yHF9AwcOlIEDB+alKERERPQSU+WoJSIiIiIRBhkiIiJSMQYZIiIiUi0GGSIiIlItBhkiIiJSLQYZIiIiUi0GGSIiIlIt3kuc6AUBQJCcnOfXZ6Zm/PvvpGTJTDfPV3k0trai0WjytQ4iooLGIEP0AgAgkX3flOSTJ/O8jifmViIdZouIyJUmTcUmIzVfZbJt0EAqr/6eYYaIXmgMMkQvACQn5yvEiIjYZKTKji3jTFQikeQTJwTJyaKxszPZOomITI1BhugF47n/LzGztS2y989MTpYrTZoW2fsTERmDQYboBWNmaytmrAUhIsoVjloiIiIi1WKNTEEARNKScreszkgTSU0SkVyONLG0E2EnTCIieskxyJgaIPJ1sEj04Vwuby0i3zz997zqIpqU3L2uYiORt3cyzBAR0UuNQcbU0pJyH2JExE6TIjdt+hr/PtGHnr6Xlb3xryUiIiomGGQK0rirIlYm7rSZmiTyn+qmXScREZFKMcgUJCs71pjQC8GYWYMzdZbLNGKmYc4ETERFgUGGqJjLz6zBxswnw5mAiagoMMiQHgCSnJ73+/1klZSWofPvZBFN/u7/o2Vrwav/3DLFrMG5wZmAiagoMMiQAoD039FfTt07Zbp1ZlqKyEwREWmxoblozNJMst5XXF6RlW1WMswYqSBmDeZMwERUlBhkSJGcnmzSECMiojFLEwevSSZdp4jIybsnJTk9WewsefVvDM4aTETFDYMMGfRHzz/E1qLo7veTk+T0ZGmxoUVRF4OIiF4QDDJkkK2FLWs7iIjohcd7LREREZFqMcgQERGRajHIEBERkWoxyBAREZFqMcgQERGRajHIEBERkWoxyBAREZFqMcgQERGRajHIEBERkWoxyBAREZFq8RYFRERELxEAgrTMXC2bmZqh9+9M0eTqdRpLM9FocrdsfjHIEBERvSQAyL1lZyQ1MiFXyycLlH//Peuw2OYyyFhVLillh9UvlDDDIENERPSSQFpmrkOMiIitaOQvKWn0+6RGJgjSMkVjZW70a43FIENERPQSKhfqb/KggdQM+XvWYZOu83kYZIiIiF5CGitzMTNxkMldzxvTYpChQgVAktOT8/x63dfmZz1atha2hdYhjYiITI9BhgoNAOm/o7+cunfKJOtrsaFFvtfxissrsrLNSoaZYgSApKcW3nVhms6ojrSUDElD4e1LFlaFNzKE6EXFIEOFJjk92WQhxlRO3j0pyenJYmdpV9RFIRMAID/MOyEx1+ML7T1TBSKlnv776/F/iVUuR3WYQrlqjtJlXAOGGXqpMchQkfij5x9ia2FbZO+fnJ5skhoderGkp2YWaogREbESjYx/WDT78t/X4iU9NVMsrQt+ZAjRi4pBhoqErYUta0GoQIV80rTYnuDTUjLkmwl/FXUxiF4IDDJEVCxZWpsX2yBDRP/ivZaIiIhItRhkiIiISLUYZIiIiEi12EeGiKgAFcS8NmkpGQb/bSqcn4bUhEGGig1jZg3O6wzBnAmYjFEY89oUxOglzk9DasIgQ8VCfmYNNmY+Gc4ETMYoinltTIHz05CaMMhQsVBYswZzJmDKq4KY18bUzVbpqRmyKvSgiJi+yYrNVVRQGGSo2CmIWYM5EzDll6nntSnoZitTN1kVt+YqAJKekpKrZfXux/UkRdIy05/7Ggtr62KzrQoagwwVO5w1mF4Gamu2Kk7NVQBk3bQJcufyhVwtn6axEKkyRERElg59Uyzx/CDjXrO29J4xl2EmFxhkiIhU7kW+HUNxvJ1CekpKrkOMiIgl0uW9G0uNeo87l85LekqKWNrYGFu8lw6DDBG90IzpB5LXYclq77/B2zEUneHLvxdLa9OFjbSUJ7J06FsmW9/LgEGGiF5Y+ekHYkwtQHHrv0GFx9LahrUmRYxBhoheWIXVD6Q49d8wJVOMijLl5H1qrzmjgsEgQ0SqUBD9QIpj/w1TKYhRUfnd1qw5I0MYZIhIFdgPpHC9iKOiWHNGhjDIEBHRMxX1qCjWnNGzMMgQEdEzFVRtWF764LxMo9EodxhkiIio0OW1Dw5Ho1FWZkVdACIievkURh8cbZ8aKt5YI0NEREXK1H1w2Kfm5ZKnGpklS5ZIlSpVxMbGRvz9/eXIkSPPXH7jxo1Sq1YtsbGxkXr16sn27dv1ngcg06ZNk3Llyomtra0EBgbKlStX8lI0IiJSGW0fHFP+0cvD6CCzfv16GTNmjISFhcmJEyfE29tbgoOD5e7duwaXP3DggPTp00cGDRokJ0+elM6dO0vnzp0lIiJCWeaTTz6Rzz77TJYtWyaHDx8We3t7CQ4OlidPnuT9kxEREVGxZ3SQWbBggQwZMkRCQkKkdu3asmzZMrGzs5Ovv/7a4PKLFi2SNm3ayPjx48XLy0tmzpwpDRo0kMWLF4vI09qYhQsXSmhoqHTq1Enq168v3333ndy5c0e2bNmSrw9HRERExZtRfWRSU1Pl+PHjMnnyZOUxMzMzCQwMlIMHDxp8zcGDB2XMmDF6jwUHBysh5caNGxITEyOBgYHK846OjuLv7y8HDx6U3r17Z1tnSkqKpKSkKP8fH/+0w1hCQsJzP0NSarpkpiQpy6dbmbibUGqiSArk/99AxCp/U3IX5vqT0pIkIznj/1edIOmWz7/V/IuyfjWXXUQkMylJHmf8u36zdNOtvyDXXdDrT0vJkOTURGXdBTGzr1rXr+ayF/T6C7zsT57Ik7S0f9efmqqKdYuIZKZmyKOUf7eNmZVpt42x69eetwHk/U1hhNu3b0NEcODAAb3Hx48fDz8/P4OvsbS0xJo1a/QeW7JkCVxcXAAA+/fvh4jgzp07esv06NEDPXv2NLjOsLAwiAj/+Mc//vGPf/wrBn/R0dHGxBE9qhy1NHnyZL1anszMTPnnn3+kTJkynC+AiIhIJQDIo0ePxN3dPc/rMCrIODs7i7m5ucTGxuo9HhsbK25ubgZf4+bm9szltf+NjY2VcuXK6S3j4+NjcJ3W1tZibW2t91ipUqWM+ShERET0AnB0dMzX643q7GtlZSWvvvqq7NmzR3ksMzNT9uzZIwEBAQZfExAQoLe8iMiuXbuU5T08PMTNzU1vmYSEBDl8+HCO6yQiIiISycOEeGPGjJEBAwaIr6+v+Pn5ycKFCyUxMVFCQkJERKR///5Svnx5CQ8PFxGR0aNHS/PmzWX+/PnSvn17WbdunRw7dkyWL18uIiIajUbef/99mTVrlnh6eoqHh4dMnTpV3N3dpXPnzqb7pERERFTsGB1kevXqJffu3ZNp06ZJTEyM+Pj4yM6dO8XV1VVERKKiosTM7N+KnsaNG8uaNWskNDRUpkyZIp6enrJlyxapW7eussyECRMkMTFRhg4dKg8fPpSmTZvKzp07xcbGxgQfkYiIiIorDZCfMU9ERERERYc3jSQiIiLVYpAhIiIi1WKQISIiItVikCEiIiLVYpChF5La+6BnZJj4HltEJABUf2wg02OQoRfOjRs3ZNu2bcrNQE0lMzPTpOsz5MqVK5KUlCTm5qa9ERs9282bN+Xs2bNFXQwqYBqNhrehoWwYZHJw48YNWb58uYwfP1727t0rDx8+LOoivXB070BuKmfOnJGGDRvK/v37lSCT3yuwuLg4EXl6p/aCrCk5ffq01KxZU1auXFlg70HZ3b9/X5o1ayZXr14t6qKogtpqNO7cuSNr166VgQMHSu/evWXVqlUF+l3HxsbK3bt3RcT4bXXr1i1Zu3atfPnll/Lbb78VRPFUr0COwXm+3WQxdubMGVSoUAEtW7ZErVq1UKpUKSxbtgwAkJmZmad1RkZGYt26dQgPD8e5c+eQlJRkyiLj/v372e4gXpCioqLQoEED7N+/32TrjIyMRKVKlTB+/Hi9x9PT0/O8zvPnz8Pd3R1vv/22SdaXk5MnT8LOzg4TJ040+boBID4+HtHR0YiJiUFKSkqBvEdWed3XDYmNjcXevXtx4MAB3Lp1y2TrBYCkpCRUrlwZR48eNel6taKjo7F+/Xps3rwZly5dKpD3yIkpvoNbt25hy5YtWLhwoQlKZLwLFy7g008/xePHj/HkyRMAQEZGRq5eGxERgVdeeQVt2rSBn58fmjdvDo1Gg1atWmH37t0mL2tERASsrKwwduxYo1975swZVKxYEc2bN4ebmxteeeUVLF682ORlBIDk5ORcb0NTMMV++Pfffyv/NvUxmEEmi5s3b6JatWqYMmWKcsL48MMP4e7ujsePH+dpnWfOnEHlypXRpEkTuLq6onz58vjf//4HwDQ7yLlz51C9enVMnDgRt2/fzvf6ciM9PR3e3t7w9PTEkSNHTLLOtWvXolWrVgCAtLQ0TJ8+Hb169cKgQYOwdetWo9cXHR2Nhg0bom7duvDx8cGwYcP0ym8qFy9ehLm5OWbPnq2s+7fffsNXX32Fo0eP5vs7iYiIQEBAAOrUqYNy5cph2rRpBRpaY2Njs5U5PwfNM2fOwMvLC3Xr1oWbmxvefPNNk+2naWlpuHXrFlxcXHD+/HmTrFPX6dOn4eHhgfr166NixYpo27YtLl++bPL3yer+/fvKv/NzjIiIiICfnx8GDRqkF+YLS2JiIsqVKweNRoOuXbti9OjROHv2rN4yOX2+U6dOoWTJkhg3bhxu3rypPP7NN9+gVq1aCAgIMOmF1MmTJ1GyZElUr14dPj4+uHDhQq5fe/36dVSuXBkTJ07EkydPcOXKFYwaNQpdu3ZVwpupnDt3Dt27d8e+fftMerFh6H20xzQgf/vhhQsXYGlpidatWyuPmfIYzCCjIz09HZ988gn69OmDu3fvKgfvGzduwMPDA1euXDF6ndevX0eVKlUQGhqKR48eAQDeeust1KtXzyRlvnXrFnx9fVGtWjV4eHhg+vTpBR5mUlNTATzdXs2bN4eHh4dJwsz06dPRqVMnAECTJk0QFBSE/v37o3v37jA3N8cXX3xh1PqWLl2K4OBg/P777/j0009Rr149DB8+XHneFD+ktLQ0fPjhh9BoNDh27BgAoHXr1qhTpw5Kly4NFxcX9OzZM8+1BRcuXICLiwvGjBmDP//8E9OnT0fVqlWxffv2fJfdkJs3b8Le3h5ubm5YsWIF/vrrL73njT2YnTt3Ds7Ozpg4cSKio6PxzTffwN3dHWfOnMlXOXVrdU6fPg07OztcvHgxX+vM6tSpU7C1tcXEiRPxzz//YNOmTahatSoOHjxo0vfJ6ty5c7CwsMDo0aOVx/JyEomIiECpUqUQGhqK2NhY5fEdO3bgxo0bJihp7sybNw8zZszAypUr8c477yjhZPXq1XrL6f4ez58/D1tbW8yYMSPbc8DTi56yZcvinXfeQUpKSr5P6NrvOiwsDKdOnYKjoyPWrVsH4PkhPi0tDdOmTUOXLl2QmJioPL59+3bY29vrhbD8un79OqpVqwaNRoNXX30VR44cKZAwc/XqVTg7O0Oj0eR7P7xz5w6aNGmCxo0bw9PTE+3bt1eeM1WYYZD5f9rU/N133+Hjjz/Wey4uLg6lS5fG3r17jVpnamoqZsyYgbfeegsPHjxQvrRTp06hcuXK+Q4cmZmZ2Lx5M9q1a4crV65g7ty5KF++fIGFGW0Q05Weno6mTZvmOcw8efJE+XFs2LAB1apVw2effYagoCDExMQAABISEjBr1iw4Ozvj5MmTuV53UlKScjB6/PgxFixYgHr16pmsZka7z1y+fBkjRoyAo6MjatSoga5du+LkyZPIzMzE+vXr0bRpUwwePNjoK7P4+Hh06tQJQ4cO1Xv8jTfeUAKfqV2+fBnt2rXDuHHj0LNnT/j6+qJfv344fvx4thrJ5x3UHjx4gNdffx0jR47UezwoKAjr1q3Drl27EBERkat16Xr8+DF8fX3RpEkTAE9rL1xcXJQTRmZmZr4P7mfPnoWDgwMmTZqk93hAQAAmT56MMWPGKLWqpnT79m34+fmhQYMGsLe3x/vvv688Z8xniouLQ6NGjTBixAi9x8PDw6HRaODn54eoqCiTlftZVq1aBV9fX6WW6ddff8WMGTOg0WgwcOBAfPXVV0hLS1OWT0pKQsuWLVGxYkWcOHFCeTzr9/rhhx/CwcEh35/j7Nmz0Gg0+PDDD5XHBgwYAC8vL8TFxT339ampqVi+fHm2C63IyEi4u7vn6QI4p/eZMWMGunfvjvPnz6Nu3bqoW7euXpgxRahJSEjAoEGD0L17dyxbtgy2trZ6+5Gx77FhwwZ07doVf/zxB7Zs2YLq1avrhRnd7z6vGGQAHDt2DJ07d9arzgX+/cLi4+NRtWpVvavqX375JVdJe/78+ZgzZ47eY9evX4eDgwNOnz6d7x3v9u3b2LVrl/L/c+bMUcKM7lVrft/n0qVLqFevHgYMGIA1a9bg2rVryg6YmZmJVq1aoWLFijh8+HCu3+vYsWPo2rWr0nYaERGB1q1bo1mzZmjXrp3espcvX4aHh0eempi0Hj58aLBmZuXKlXpXUrkte+fOnXH37l0AwLVr1zBs2DC89tpr2Zo4FixYgFKlSum1EefG1atX0atXL/zyyy8A/q0JmzNnjt6BwJQePHiAoKAgfPbZZwCAo0ePIigoCB07dkTr1q1x5MiRXDdrxcbGYsWKFTh16pTy2KxZs6DRaPDKK6/glVdegYuLi97+mxvJyclYv349atasiQ4dOuD8+fNo2LAhDh06hMzMTDx8+BBxcXF48uQJUlJScPbsWWXb5dbIkSOh0Wjw22+/Kfv5rFmzYG5ujvbt26Np06bQaDTZftv5kZmZie+//x49evTA/v37sWbNGlhbW+ODDz7QWyY3jh07hvr16+s1vaxevRrW1tZYvXo16tevDz8/P0RGRpqs/M/Srl079O3bV/n/bt26oWLFiujSpQv8/Pzg4OCA1atXK59v27ZtaNq0KXr06IF9+/Ypr8vMzFS+j6NHj8LBwUHveWNpt/lHH32k/D8A/Pjjj/D09MRPP/0E4PkXPIaOtY8fP4anp6fe8eDPP//MUzm1F0G//PIL1q9fD+BpTZFumDFUc5SX4358fDw++ugjbNq0CRkZGfjf//4HW1tbvPvuu3la75MnT5TtmJGRgc2bN2cLM/mtmXnpg8ypU6dgY2Ojd7AA/q1OzMzMxOPHj1GnTh2lOnzixIkoV65cjp0Wk5KSDCZ57ZcfExODKlWq4Nq1a8pzf/75J+7du5erMiclJeH+/fs5Vnkaqpn5/PPP89xRMS0tDZMnT4ZGo4FGo0GXLl1gbW2Njh07IiwsDFeuXEFycjLatGmDOnXq5Pij0pXTdp87dy7Mzc1RpkwZ5Wpd+5n9/f2VH4Qht27dwrZt27B69WrlxKUbtoCnJ2rdMPPee+9Bo9EYVf2bU9kvX76Mffv26TW9AcCmTZtQu3Zt/PPPP7lav+4Vyh9//KH8W7u+b7/9FsHBwXqvyc2V4/Nov7Pdu3ejUqVKeidBX19fWFhYoGrVqmjbti1mzpyJ5OTk55ZfNyCuW7cOFhYW+OGHH5CQkIBr166he/fu6Nq1KxITE597cLx06ZJyIkhKSsLmzZvh6ekJd3d3aDQaVKhQAW5ubihTpgxKlSqFcuXKoUyZMihfvnyuf1vR0dHK5+rcuTPKli2Lw4cP4+OPP0bp0qXx888/IyMjA0lJSfjggw9ga2uL69ev52rduREVFaUX1rXhw9iamRUrVsDOzk753BkZGdi3b59SaxobG4s6derA29sbDx48MFn5s5ZPu09t3LgRXbt2RVpaGvr37w9XV1dcuXIF6enpSn+S8+fP653Qtm3bBn9/f/To0UOviVO7/pUrV8LLy0uv2cwY58+fx6ZNm/Q6z+uW3d/fP9vvTEt7rPn++++V33vWGqM7d+7A2dkZ586dAwBMnToVdnZ2Rl/QREREYPjw4UhJScn2m0tNTTVYM/Pbb78Z9R7A0ybNtWvXAtA/nqSmpmLDhg3ZamZSU1Of23Hf0HngyZMn2LJlC6pVq6YXZv773/8a1S9J10sdZE6dOgU7OztMnjz5mcs9ePAA7u7uOHr0KMLCwmBra5tjM8r58+fRsWNH1KlTB61atcLq1auzfZl3795FtWrVlHbqiRMnonr16rnawXXXHxgYqLd+3YOANsyEhYXhrbfeQsmSJfNVxXn+/HlMnDgRTk5OWLt2LQ4dOoTJkyejYsWK8PT0hLe3N6ZMmQKNRoO6desq/UUMed52Dw8PR5kyZdCiRQv89ttvuHLlCqZMmYLKlSsjOjra4GtOnz4NT09P1K1bF2XLloWvr2+2q3Dtdnr48CHmz58Pc3NzODk54fjx47neDs8ru6Eri9GjRyMoKChXncWvXr2KKVOmYMSIEXo1Fbr70Oeffw5vb2/lvT7++GOEhITkqVPh7du3lc6XmZmZSE9Px8OHD9GzZ0989dVXAICBAwfC3d0dkZGR+OWXXzB+/Hi4uLgYbL7Ulv/dd9/Fzp079Z5LT0/XayoAgLfffhuBgYHPLWdGRgZGjx4NjUaDPXv2AHgakn744QcEBASgRIkS2L59O3bu3In//e9/2LZtG7Zs2YKtW7fqXTA8y99//w13d3e97d6+fXtoNBqULFkyW7+k5cuX5/p3m1fp6enZambS0tLw/fffZ+s0e+vWLaVWee3atbCysjLYb0h7sluxYgW8vb1NVv4HDx4ooSBr2IqLi0OVKlXg7OyMSpUqZdsPUlJSkJ6ejgsXLiAhIUF5/Mcff1TCjG7NS1paGoYNG4Z+/foZPQI0MzMTycnJ8PPzQ3h4eLbntb+1n3/+Ge7u7tixY4fe87k51mRmZiI6OhplypRBZGQkZs+eDVtb22ceF3MyZswYdOnSJdvj2guGlJQUJcwcOHAAw4YNQ/369ZXa4txISUlB48aNlT5Jup8DeLofasOMtmZmxIgRGDJkiMELmueF7eTkZKWZ6Y033sCoUaOg0Why/VvN6qUNMjdv3oStra1ycNDuvPPnz8eSJUv0ln348CHq1auHli1bwtraOsedUdtJrG/fvpg2bRo8PT1RuXJl/Pzzz3rL3bp1C05OTrh06dJzg5Gx69c94c2ZM0c5CGc9cOTF+fPn8c4778DR0VFJ/I8ePcKxY8cwadIk9O7dG9bW1tBoNDl2JnzWdl+0aJGy3FdffYW2bdtCo9Ggfv368PDwyPEzaGtIpkyZgps3b+KPP/5ApUqVlBOeLu37vf3223BwcNCr9XmeZ5Xd0DDL6OhoTJgwAU5OTtlOOoZoh/2PGDECCxcuzHGY9ZdffomAgAAAwLRp06DRaIzqO6T16NEjdOnSBc2bN8fp06f1nvvkk09QrVo1dOvWDe7u7tn2T0P9pZ5V/qwBT3v1OnjwYIwdOxbp6enPPfjdu3cP77zzDqytrZWht48fP8YPP/yAWrVqoWfPnkZ9/qwuXrwIV1dX/P3333onpn79+sHW1ha7d+/W+0xjx45F8+bN8fDhwzy/Z26avNLS0vTCzKhRo2BpaanXLKQNYb///juAp6GifPnyaNu2rbJM1v3pvffew5tvvml0s6oh58+fh7+/P8LDw5VArf0+tb+TdevWoXLlyti8ebPea69evYp3330XdevWhYWFBWrXro2RI0cqr9u6dWu2MBMWFgY3N7c8X8EDQP369bF06dIcn4+MjISnp6dy0ZKZmWnUsSYhIQGvvPIK3njjDVhbW+e5w3///v3Rq1cvg89p95/U1FT4+PjA3Nwctra2Rl2caXl7ez9zuHhGRgY2bNiAkiVLolq1ajAzM8t2LtTdl3L6PWsfT0lJwaZNm6DRaFC6dOk8lVnrpQ0y69evR40aNdCzZ0/loPzxxx/DyckpW7Xc3bt34ezsjDJlyui19+s6d+4cHBwcMGXKFOWx5ORklC5dGoMHD9Zb9u+//0bt2rXRtWvXZwajvK4/IyMDGRkZ+OCDD1C6dGmlatMY0dHRWLNmDZYvX653hXrlyhUMGTIEjo6O2LJli95r0tLScPPmzWe2uxuz3ZOSknDy5EmcP38+x+rjixcvwt7eXq+jHvC0KSQ8PByjR4/G7t279ZoW1qxZg7Jlyxr9wzGm7EePHkX37t3h5eWVq5Bx9epVuLu7Z5uHRvdgoP332rVr0a1bN0ybNi3X+09OVq5cieDgYHTs2FFv387IyIC/vz/c3d1ztZ1yU35dqampmDp1KsqVK/fck5FuCLp//z4GDx4Ma2trZb9MTEzE5s2bUa1aNbzxxhvPfe+cHDp0COXKlVMOxrpNZB06dEDZsmWVWpnp06fDzs4uX6OvIiIi0LFjx1z9PtPT07Fq1SpoNBo4OTllOynqhrDMzEykpqZi7ty5sLe3R7du3fSWjY+Px+TJk1GmTJk8HRuyioyMRP369eHs7IymTZvi008/zRZmgKdBt27duvj666+Vz3T69GlUqlQJgwcPxrx587Bnzx5069YNzs7OaNy4sRJmfvrpJ/j7+6Nfv37o27dvnk/WWpmZmfDy8npuh+2lS5fC0tISZ8+eNfpYEx0dDY1GAxsbmxzPGznRvSAdNGiQ0rRoaJ/W7qfDhw/P83eqnU5j06ZNz1wuMTERrVq1QpkyZbJdnEVFRaF37956QfV5v8EhQ4bA3t4+3/vhSxtkUlNT8f3336NRo0bo3bs3wsLC4OLikuOw1jlz5uR49Z6eno7WrVvD0dERhw8fVtYPAD179kSfPn30hghGRUXB0tISZcqUydVJLrfr17V7925oNJo8XQUYmtjp888/V56/du0ahg4dCkdHR2zbtg2Afie8Z8nNds/tCSgjIwMhISEoU6aMMjoJAGbPng0LCwu0a9cO9erVg4WFBT755BPl4JCUlJSnkQ7G7jO7d+/O1ftkZmZi2rRp6NixY7YO54YsW7YMGo0G9vb2+QoxWmvXrsXrr7+uF2bS0tIwadIk+Pv7Z+tnlN/y79ixAyNGjICLi8szawrj4+OVf+cUZrR9ZhITE7F161aULl3aqJqZf/75Rwmlhw8fRpkyZfT6n+m+b4cOHVCxYkV0794939v+xo0bqFq1KjQaDXx8fJ7bfy0jIwODBg1CyZIlDc6XkzWEAU+bcyZNmgQHBwfUrFkToaGhePfdd9GxY0e4urrmKwhoZWZm4rPPPkObNm1w+PBhvP322/D399cLM7rb8MMPP0TJkiVx9+5dnDp1Cvb29srcK1qPHj3C4sWL4ezsjLZt2yr73Y4dO1C9enU4OTnlqYb5+vXrOHjwIJKTk5GWloZKlSplqy3P6tChQ2jRogVu3rxp1LFG20Q7efJko/smRkZGYvTo0UrTbffu3TF16lQA+n03dX+PCxYsgEajMWq7REZGYseOHUrTkIeHh3I816VbsxYaGgqNRmMwmB09ehTe3t7o2LGj3nbN6bjx66+/wtPT0yTHsJcyyGgPzGlpafjuu+/g5+cHjUajJEndH97zOq1q0/fVq1fRoEEDBAUF4ddffwXwtEOdtbW1werL999/P1fNDXldPwCj2ki1njWxk25btDbMODs744cffsjVuo3Z7rkVHR2Nbt264bXXXsPWrVsRHh4OZ2dn7NixQzk49u/fH2XKlMnT9shL2fMyUqB58+bo16+fwee0++CjR4+Qnp6O7du3o2bNmnmqVr916xZ++OEHbN26Va/5b/369UqY0R4MY2NjUapUqWzTEeS1/No+Qjt27MD48eOfWf6///4bbdq0UUZPAdnDzIABA+Do6KhczSUnJ+Pnn3/OdV+wf/75B87Ozvjoo4+QmpqKo0ePwt3dHQ8fPkRGRka2phEASnNnXprytJ48eYKwsDB06dIFR48ehZ+fH7y8vJ55wtu+fTs8PDz0LkyeFcK0ZX/48CG2bduGNm3aoF69emjYsCHGjRtnkhmKtRdTt27dUmo2UlNT8fbbb8PPzw8LFixQjhna38+pU6fg7++PgwcP6jXT6v6+gKfBdMaMGShbtqwySgcA9u7dm+fO1W3atEGpUqWwf/9+pKenw8PDI1v/F0O0nfTzcqzJyyzcS5YsgZeXF9555x3Ex8ejc+fO2WqBsjp9+rRREzVmZmaiY8eOqFmzJrZt24a0tDRUq1bNYJDRSkpKwpw5c7LVQsbGxir9mg4dOoTmzZujffv2OXZ5SE9PR2RkJDIyMvLcUTurlybIxMfH613haTes9sTk7++Prl27Ksvk5qQaFRWl13/g2rVr8Pb2RocOHbBmzRqUL19ebw6NzMxMo2ZIzcv680Nb3Z/biZ2uX7+Ovn37onLlynj8+LHB9y+I7a6lXfbWrVvo1KkTPD09YWtrqwQ97UFk5cqVqF69eo4dhXNSkGXXXWdCQgL8/f0xbtw4vXJn9dFHH+HAgQMA8hZST58+jSpVqsDX1xfm5uZ466239LbJunXrlDCjPWGGhIQgMDAwx47KxpR/xowZSvlzGvGkdeHCBXTp0gXNmjXDl19+qTyuu42vXbuGli1b4sMPP8zz8M2wsDBYWlpi3rx5+Pbbb1G3bl0kJycjJSVFGSWSlJSE5ORkpbYpP51jtSFj06ZN2LhxI4Cn/VmeF2Zu376t977PC2G6YUbr4cOHePLkiUmmtj9z5gxmzpyp9/vQSk5OxqBBg5Qwo/2utTUZDx48UJppe/XqpZwEtd+httxxcXFwc3N77kncGM2aNUOVKlWwfft2eHp6YunSpTh16hQOHDiAffv24dChQzhy5Ag2btyoV7th7LEmv3PbLFq0CH5+fhg6dCj8/PwQEhKC2bNnY+7cuViwYAHmzJmD2bNnY/LkyVixYkWe3kM7z1PDhg2xatUqVK1aFStWrMD+/fvx559/6v2tX78eFy5cyLZP3b9/H8HBwejXr5/SV0xbi9WuXTu9YJSZmYmUlBSMGDEi2zkmv16KIHP9+nVUqFAB/fv3VzqU6n4hqampWLlyJRo1aoROnTpl+2Hl5M8//0SVKlWQlpamLHv16lV4e3vDzMwMnTp10uv1bayCXr8u7ZXQF1988dyJnXS33ZUrV3I8sBfEds/6Q9KW+/bt2+jevTt8fHyyzRj63nvvoXnz5gYPujkpqH0mJ4MGDULp0qWVSQCznmyuXbuGdu3aKUHAWLqjrR4+fIhff/0VGo0m2ySPq1evxuuvv46uXbvi4sWLiIiIyNVIgtyW/9ChQ89cz61bt5Qr4HPnzmHAgAEICAjQCzO6TZjt27fPsSboWXS/yzlz5sDCwgLBwcGwsLBAuXLlULVqVdSpUweVK1dGhQoVUL16ddSvXz9fHXsvXbqkdBzNGvbi4uKUMKO9sk5LS8P+/ftzDH3PC2FJSUlISkpCSkoKoqKiTDLxGPB0X9JoNJg+fXq257T7f1JSkhJm5s+fj2HDhkGj0SjBWbeZtmPHjjleDPj4+OjNLGsKfn5+KFGiBCwtLaHRaFClShU4ODjAysoKZcuWhYuLCxwdHbPV/BTUsQZ4Go61E8Zpa7oWLlyIgIAAZcBGq1at4OnpCU9PT9SvXx8+Pj6oU6dOrmr2te7evYvDhw8rv/uEhAQ0a9YM1atXV6bX8PDwgKOjI6ytrVGmTBm4urqidOnSuHr1arb1aWc2b9q0KUaMGGEwzGhrZjIzMzFy5EiYm5ubpFlT10sRZNauXQt7e3usWbMGVapUQY8ePTBp0iQ8evRIrwf1t99+iyZNmqBFixYGR2Vk9dNPP8HDwyNbe3pkZCReffVVBAUFKaMIAONrTAp6/VoRERGYPXs2UlNT9ar6njWxk+775sTU2z0iIgItW7bEpk2bDP4Qbt++jU6dOuG1117Dt99+C+BpLYaDg4PRnTILap8Bnp6sf/rpJ6xfv16ZXG737t0oV64c2rRpY3BYc1hYGJo2bZqnqtgLFy7AwsJC774pANCyZUvMnDkTYWFheld169atQ4MGDdCnTx+Dwaygyp+QkIC2bduicePGyjwWERERGDhwIAICArB8+XJlWe3Bvn///vjoo49yPZvv3bt3cfbs2Wx9xz799FNoNBpl5M1XX32F1atXK8F+3bp1+RohAzzdrhqNRrmKB/R/s/fu3VPCTEREBIYNGwY/P79scwTlJYS9+uqruZ7H6FlOnjwJW1tbvUEHWtp9Rftfbc2MtbW13sjJrM20WS8GtM9fvHgRjRs3VubVMfb4du3aNSxYsADDhw/Hrl279GoV27VrBysrKyxbtgzXr19HdHQ0Ll26hOjoaNy+fRt79+4ttGPNmTNnUKNGDXh4eECj0aBp06ZK4P/ss8/w2muvoWfPnsr20e1PZMyUC+fOnUPjxo3Rrl07dO7cWQnICQkJaNeuHTw8PDB37lxER0fj/v37uHr1KmJiYnD//n2D8zBpv+e0tDR8/PHHykzShsLM1q1b8f7778PW1tYkI2izeimCDAC88sorWLp0KR4/foyFCxeiffv2qF69ul51fXp6OlauXImWLVvm2AyRkJCg7AD/+9//4O7ujsTERCVsaP977do11K9fH+3atVNmZs2Ngl5/Vtqrq6z9IEw1sZMptru2LP3791emEq9RowY++uijbCekqKgodOrUCa1bt0ZQUFC+RvWYap/RdebMGVStWhU+Pj7QaDRo1KiRMqQ0NDQU7u7uyiy1t27dwoEDBzBixAiUKlUq2xDp3MjMzERYWBg0Go3eBHcff/wxNBoNunXrhnr16mWbQXbz5s0GJwksyPJnZGTg66+/xmuvvYa2bdvqhZkBAwagUaNG+OSTTwA8DSTTp0+Hi4tLrvt7nDlzBvXq1UPt2rWh0Wgwbdo0vec/++wzaDQaLF682OQ3+svMzERCQgJ69eqF5s2b5zg9QVxcHAICAmBmZgZbW1u9/bsoQxjwdMJHMzMzzJs3T/lMwNPO51mbUrTHqWHDhsHJyQmHDh16ZjOtNszoLjNhwgT4+/vnqSnvzJkzcHd3R1BQEPz9/WFpaZntbtmNGjVCzZo1DU62V1jHmjNnzsDW1hahoaE4fvw4/vzzT5QtWxZt2rRRlvnPf/6DRo0aYdCgQUqNZ9ZmuOfR3ndrypQpiI6OVl6nDY2PHj1C8+bN0bhxY/z000/PbH7UfU9tOVJTU5UwM3z4cOV7PHz4MAIDA+Hs7JzvkWbPUuyDjHZDf/311+jZs6deNW358uVRuXJl2NraYtSoUVizZg0A6E3IpOvOnTto3Lix0vSinVkU+HfUTnp6uvKely9fRqVKldC5c+dctQcW9PqzOnnyJOzs7LK1QWcd8puXiZ1Mud21V9+nTp3Ca6+9hu3bt2PPnj1o1KgR2rRpgw4dOuDs2bNKv5GoqCi0bt0aZcuWzVOnTFOWXZf2oDV16lTcvn0bR48ehb29PQYNGqQsM2/ePNSvX1+pTq5Tpw78/f3zFGKio6Pxzz//4N69exg6dChsbGxw5swZfP7553ByclJmSY6Li8OwYcPg7u7+zGGQBVl+7TZPSUnBmjVr0KJFC3To0EGpRbhw4QLGjh0Ld3d3lC9fHk2bNkWtWrVy/f1qaxImT56MgwcP4ssvv4SFhUW2CeNmz54NjUaD//znPyapwcg6T8y6devg4eGh1C4ZOmFoR8fofhdFGcKApye8+fPnQ6PRKMOngaeTV5qbm2e7uSjwtJlao9Hgp59+MqqZNjMzE3PnzkXJkiXztN9fuXIF7u7uCA0NVX67wcHBmDVrVrZlGzVqBE9PT+zZswcZGRmFeqy5fPkyNBpNtnuRLV26FPb29no14AsWLFBqZoydxfv+/fto1qwZ3nvvPb3Hs3ZLSEhIQIsWLdC4cWNs2LDBYG3sxYsX8emnn+rNBq0bZmbPno2AgACEhYUp2/7IkSN44403jJqzy1jFOsjo/lgiIiLg4uKCH3/8EcDT2Urd3Nxw9OhRbNu2Db6+vqhfv76SeHPSvn171K1bF6tXr8Ynn3yC119//ZnLx8TEGNXLvqDXr3X27FnY2dkhLCxM7/EvvvgiW7PRo0ePjJrYyZTbPSIiAh999BHi4+MRGxuLXr16YdmyZQCedla7efMmNBoNvL290apVK/z8889ITk7Go0ePnjt9dkGXXVdOBy0fHx94eXnpXYk+fvwYmzdvxnfffYcjR47kenp9Xbdv30bFihWVE8z9+/fx9ttvQ6PRwNraWhnGrz2RfvXVV6hSpUqONxstqPLr9tvQ7TfSq1cvWFtbo3379sqBOzY2FufOncMnn3yCrVu35vo+QefPn4e1tbVereOdO3cQEBCAbdu2YeXKlXrNAfPmzYNGo8Hnn3+erw70Fy5cQKNGjbB79269Zo0BAwbA3d1d2Wa6Q2q1gUS3+r2oQlhWN27cQFhYGBwcHLB27VosXboUpUuXzjaDs67r168b3UxbtmxZWFlZ5akm9cmTJ/jggw/w7rvv6l2A9OrVCyEhIejSpQuWLl2qVzPg5eUFHx8fHDt2rFCPNadPn4ZGo8GgQYP0+jB98803KF++PK5fv64XdD/++GMEBwfn+l5nWhcvXkS1atXw22+/Pfd+TPHx8fD29kZgYKDBpvKvvvpKuceYbn8x3ebEUaNGoX79+nrHkryM3jJGsQsycXFxuHDhgsGpuefOnYsOHTrgjTfegKurq97BIiYmJseke+PGDSxevFjp7NqjRw/4+/srQzEDAgJQt25dNGrUCA0bNkT9+vVRt25dBAYG5mr67IJev6FtVKZMGTRr1kzv9dqZgLNeXd26deu5EzsVxHbXNnvp9u1YtmwZnJyclL4WAwcORMWKFfHFF18oN/rLOlQ8N9vD1GXPSnvQevvtt5UmqLlz50Kj0aBatWoYOHAghg8fju+//z7X5X6WyMhIlCtXTq956N69exg/fjwsLCyU5iDtAWjMmDFo2bJljvfdKYjyX716FZMnT87W7PGf//wHTk5OmDhxIpo2barXzGSs5ORktG/fHi4uLnpXzNqbVzZp0gTm5uaoXr26Xtk/++yzfE/SNXXqVGg0GlSvXh0zZ85URik9fvwYtWrVQo8ePbLVfv7+++96w8eLKoRp6Q6tBp52Sg0NDYWDg4Ne4Mp6gszasdiYZtrXXnstTzUxWvv27dM7hs2YMQNWVlYYOnQoQkJCUKVKFQwdOlQvfP/888+FcqwBnu6T2pqfAwcOwNraGr1790ZycjKioqJQtmxZvT5Iuts2N/M0ZbV+/XqYmZkp72kozCQmJipzMiUkJGRrVr579y6OHj2K+/fvK329Zs+erRdmtN95cnKycnPSwlKsgszZs2fxyiuvoFatWtBoNJg6dape58J9+/ahfPnyqFmzpt7B4lk/eG1HrC5duujNZNunTx8lZEycOBEfffQR5syZg/DwcEydOhVz587N9bT0Bbn+nAwbNgw1a9ZUJrr7z3/+gzJlyuh1RNR68uQJJk2aZPBEDxTMdj937hxsbW2VGiPtsqmpqejevTu+++479OnTB66urnoH8t9++82o+3UURNl1GTpoDR06FJMnT0bp0qWxefNm7N+/H7t370bXrl1Rr149ODg44P3338eTJ0/ydDJKT0/HkSNH4OzsnO3KPC4uTumAqZ1WPTQ0FPb29gZPHgVZ/k2bNsHKygqjR49Wtu2cOXPg5OSk3Cxz9erVaNy4MTp27Gh0mNF+9i1btqB9+/Zo27Ytrl27hs8++wxOTk748ccflRuwVqlSBe3atcvThUFOUlJS8NZbbyEwMBCTJk1CYGAgRowYgZSUFCxcuBCNGzdWZlI1tJ2KMoQBT2skGzRogMDAQKVWEnha2zdz5kw4ODjojSYzdIIsqGbarFJSUgxOEXDjxg20bdtWudkn8PQeUzY2Nrh69SoyMzML7VgDPL3A8PX1xZ49e5RaigMHDsDKygodOnRApUqVMGzYMGX5rH0jjaH9HIcPH4aNjQ3WrFmT4+9x2bJlaNWqlcFteO7cOTRp0gSBgYHKPZ8WLlyohBndQJiamoro6GgEBATk667kxio2QUY7S+SECRPw22+/ISwsDObm5srdPLVCQkLg4+Oj/P+zDrQXLlyAk5MTJk2aZLDKvV+/fvDy8sK6devy1C5d0Ot/nvfeew/Vq1dHhw4d4OTkpCRy3W2yefNm3Lp1K8dhxQWx3c+ePQtnZ2d4eXkpj+le4U2YMEG50s1Pu2tBlF2XoYPW/v37YWdnB41Gk+19njx5gn/++QezZ882anIrLd2rtQMHDsDZ2dngyT8uLg6DBw+Gg4MDunXrluMstQVVfu3VXUxMDH755ReUL18eH374IcaNG5ctTGdkZGDt2rWoXbs2evbsmesD+rVr11CrVi2lHFu3bkVwcDBq1aoFW1tbZVSI9nNNnjwZ9evXN3rY7LOkpaXh22+/xfDhw/HHH3/gyJEjyi0uFixYgKpVq+r1L9JV1CEsIiICZcqUwbhx4wyerKOiohAaGooSJUrgv//9r/K4oVtqaNdnqmbarC5cuICuXbvC19cXrVu3xsGDB/VGRmm3pfYYtmvXLtSrVw937twptGONrho1aqBWrVrYt2+fXpgpUaIEypcvb7CTvTG05wrtf+/cuYPy5csjKChIbwi17vczbtw4jBs3LtvvS7eTcGRkpF6fr0WLFkGj0WDmzJl6zV3Tp09H7dq189TcllfFIsicP38eFhYWyjTOwNM5G5ycnNC7d2+9ZY8ePYr69etjw4YNz1xncnIyevToodzpUys1NRU3btxQvrghQ4bA09MTq1atyvXw28JYf1ZxcXE4f/58tllPx40bB3Nzc4wYMSJbGp8yZQqcnZ1z/GEVxHbXznfSokULuLu7Y9SoUcpz2h/Ro0eP4OvrizFjxjz/g+egIMpuSPXq1bMdtI4fPw4bGxv0799fr59HfiYqe/jwIcqWLYtJkyYp71GlShU8ePDA4AiHhw8fIiQkBHZ2ds8cSWDq8muv7lq3bq1c3X377bdwcXGBpaUlvvnmG2VZ3eG8GzduzHGkjyG///47ateurXfg/emnn9CiRQv4+fkpn1m7Td566y306NHjuRP1PcvVq1fRqVMnnD17VvktRUVFoUWLFso9qNLT0zF+/HiMGjUKbm5u0Gg02ZrjijqEPXz4EE2aNNEbxaalO6jgwYMH+PDDD1G6dGnMnz8fQOE00+o6deoUSpUqhb59+2L48OHw8vJCtWrV9Jors154jB07FoGBgfjrr78K5VhjSMOGDVGtWjW939WRI0dgbW2NN9980+jJO7UuXLiAvn37IigoCCNGjFBG9P3888+wtbVFjx499GpeHz9+jMmTJ6NSpUrZvrP79++jadOmetsF0A96ixcvhqOjI3r06IHBgwdj8ODBcHFxMfreUvml+iCTlpamdEDSveHVzJkzodFo8Prrr2POnDnYsmWLcmCoVasWQkJCnjl5WVpaGpo1a6Z3j6GdO3fi/fffR8mSJVGhQgV07doVwNOw4erqqjeVdm7KXZDr13X27Fl4e3vDy8tLaT/X9cEHH6Bq1apYtGiRclU/derUbMM/s5bf1Nv96NGjsLS0xPTp05Geno4vvvgCzs7Oej+klJQUpKWlYfLkyXjjjTfydPAuqH1GSzuDpZafn1+2g9b+/fvzfdDS0r0Lt42NDWbNmoUDBw6gVq1ayrozMzPx+PFj5SotMjISt27dMtgRt6DKn/XqTveA+L///Q9ubm4YM2aMXm1OXicYXLNmDTw8PLI1b/34448ICgpC27ZtcfDgQQBP93UnJ6d8XXFnZGRg7969yjwu48aNUy4ajh07BgcHByWwZGRk4NixY5gwYQIsLCyyvW9RhDBd169fR926dZUO4cDTJtapU6eiSpUqaNCgAb788kukpKQgNjYWY8aMQcWKFbF///4CbabN6syZM3BwcEBoaKjyWExMDEqVKoV33nkn2/KxsbEIDQ1FqVKlsG7dukI51gBPw+xXX32F5cuX691c1lCY2b9/P0qUKIEOHToYXaNx6tQpODo6YsiQIRg4cCBee+01jBo1Cunp6UhNTcWKFStgbW2NqlWrom/fvhg8eLDBUKl17tw5VKtWDX/++We2ixTdmaO3bNmCESNGoHXr1hg7dqzBe4EVNFUHGe2IlsePH2PSpEmwtLTEzp07MXfuXJQqVQpLlizBsmXLMHHiRFSrVg01a9ZEz549MWPGjOcetOLj41GrVi0MGTIEFy9exOzZs1GzZk1069YNixYtwooVK1C5cmXMmDEDwNORCMa0lxb0+rVOnjwJe3t7jB8/HocPH1Y6aGZtynr33Xfh4eGBL7/8EmPHjoWNjU2OowYKarv/+eefegeShw8fGjzAAP92BM7avPE8BbnPAE9rdUaOHIkuXbrodRw09UFLSzuK58mTJ0hJScHSpUthZmYGf39/lC5dGq6urvDy8kLt2rVRoUIFVKtWDdWrV4e3t7fBA3NBlT83V3erVq1C+fLlMWrUqFzfL0mXbt+KL7/8Ep6engbvwrx161YEBQWhS5cuePPNN5+5r+fGhQsX9GovJkyYgEaNGsHd3R0bNmzAkydP8N///he9e/fO1nfF0EzBhR3CsoqKikKlSpUwZ84cAE+bEHx9fdGqVStMnToV3bp1g6urq3KR8/fff+O3334r0GbarNLS0hAcHAyNRpNtVuDAwEC8++67eus+cuQIgoOD4enpiZMnTxbKsQZ42km+cuXK8PPzQ5kyZVCtWjWlLxAA+Pv7Z/td/fnnn3B1dc1xBKEhZ86cgb29vV6oGzt2rHJDYW3z2unTp9G3b1/4+fmhRYsWmDJlSo6/tdWrV8PCwsLgPce0EhMT9Wq/THH7i7xQbZAxNKJl/Pjx0Gg0sLCwUDoyal25cgU7d+5Eq1atcl1FvWfPHlhYWKBy5cpwcHDAsmXLlC89NTUVQUFB6Nu3b54/Q0Gv/9y5c7C2tkZ4eLjy2N27d9GkSRP88ssv2Lp1q14nwtGjR0Oj0aBEiRI5NjcUxnYH/j3AxcfH53iAGTdunFHpv6DLfurUKZQtWxadO3dG7969YWlpqfdefn5+8PDwyPdBK+vnmTt3rvJYSkoKvvrqKzg5OaFy5cpYs2YNFi1ahM8//xyLFy/G0qVL8cUXXxicIK0gy/+sqzvdWXm///57VKpUCSEhIUYF9zt37qBNmzZYsmQJAOC7775D48aNleczMjL0+plt3rwZr776KpycnPI9Sdd3332HGjVq6DW7HDlyBMOGDYOZmRnefvttTJo0CSEhIUqtqnb7aT93UYUwQ+Lj4zFmzBhUqFABlStXho2NDT755BO9wQVly5ZV5rIprGbarK5evYqaNWuiadOmSvP37du3YW1trdcJWWvr1q0Gp6ooiGMN8DQ02NnZYdKkSUhMTMSuXbtQvnx5tG/fXi/Aai8S/vrrL+U7N6Z2TTuqNOu5YvTo0fDw8ICXlxcqVaqkzE+W21tV7N+/HzY2Nnqj1rL67LPP0Lp1a5PVBuaVKoNM1l7murRzKWiHOgL5u5liVFQUjh07lq0KPiMjAz169EBoaGiup0cvzPUnJyfjjTfeQIkSJfQ60M2YMUOZBdTa2hp169bVu9KYNWtWjld3hbnddekeYHTbqo2Zm6Cgy3769Gm9qdszMjIwcuRIvP/++3o1Hy1atECVKlXyfNDS0vYlMjRVfEJCAlasWAELCwssXrz4hSh/bq/ubt26hc2bN8PLyyvXHT9TU1MRExODDh06oEmTJli/fj2WLFmCDh06PLNpav/+/UaF65z8+uuvcHBwMDi/x8aNGxEUFKRMEli5cuVsNWFFGcKAp/1aTp8+jUuXLilhLCoqCnv27MGyZcv0Ooimp6cjNjYWTZs2xZo1awq8mfZ5rl+/jurVqyMwMBCHDx9GpUqVsvU7NIYpjjXA0+3n7OyMHj166D3esGFD1KhRAw8fPtQLFK+99hpKliypDEU39vhTv3591K5dWxl2/sknn8DOzg7Lly/HqlWr8M477xi8WHvW+9y6dQsuLi7o2LGjXj9J3deMHTsWkyZNMtmxPq9UF2QM9TLPOnvmhAkTYGlpqVeFp2WKDZ6SkqJMx56X0SUFvX5tn4U9e/agVatWaN68OZKTk/H555+jVKlS+OGHH5CcnIyjR4+iRo0a6N69+3OHPBb1do+Pj8eXX34JjUaDCRMmGPXagi57TgetXr16wcfHB7Vq1cLrr7+ujNho3rw5nJyclOYBY7fN+fPnYWlpqVfTBgAbNmxQRq5omzPMzc3x0Ucf6X2WrO9XGOXPzdXdokWL0Lp1awCGm1wMiY6ORvfu3XHq1Cn8/fffSkdHHx8fODk5oXXr1mjRooXyO2jZsiUaNmyIrl27ZtsH8iIjI0MZOaS9S33WkHbp0iWsWLEC5cqVg0aj0TspFHUIO3v2LHx9feHh4YEKFSrg3Xfffe5cJdOmTUONGjXw66+/FmgzbVZZhyJr/3vjxg3UrFkTGo0GAwYMyLa8sfJzrNG6ceMGGjZsiI4dOyrhQnvB5Ofnhw4dOiAkJAQLFixQbkHTpUuXPDWpajVs2BC1a9fG8OHD4ezsrHdLBu28Uto71OfWpk2bYG1tjX79+uk1iyYmJmLy5MmoXLlyrm8RUpBUFWSeNaIl649+4sSJsLe315tO2xRWrVqFUaNG5dhBqqjXn5iYiEaNGsHf3x8AsHfvXrz22muoXLkySpQooRxstSefESNG4JVXXnnmVfWLsN2Bpye3b7/91qgfTmGU3dBBKzw8HHZ2dpg5cya++uoreHl5oUqVKsoon9dffz3PB62JEydCo9HoXY1rJzPU3WdSUlKUaeK198cpqvLn9upu/PjxRtVAXrt2Df7+/mjXrh0iIiJw69YtvPnmmyhbtizKlSuHGTNmYPjw4Rg+fDjGjRuH9957D++++26+5mC6ePEiNm7ciL179yqjkzw8PPDpp5/qLZe1Cv+ff/7RCx9FHcJOnToFBwcHfPDBBzh06BAGDRoEe3t75QaiWb+H33//He+//77SWbYwmpi1zp07h2HDhuHw4cN621U3zPj4+KBx48ZGz3xrSF6ONVldvnwZbdq0QceOHTF48GCULVsWGzduRGRkJDZv3oxZs2bB1dUVLi4ueOutt4y6oMmpA3GzZs2g0Wiy7YsPHjxAw4YN9YbK50ZGRgaWLVsGCwsLpSZt+PDh6NixI1xcXArkHJgXqgkyuRnRkvXE9O6778LFxcVkwxIvXryIFi1aoEuXLgXSM9sU609LS8PWrVtRt25dBAcHA3gaZpo1a4a6desqP3LtgbBfv37o27dvjkHmRdjuuoz5sRdm2bMetFxcXPRu5hkZGWlwxJgxtAEgNTUVffv2hb29PS5fvozPP/8czs7OBiczTElJwddff/3c/akwyl9QV3eXL19GcHAwgoKCcPbsWcTGxuLNN99EkyZN8MMPP+S5vIZkZGRgyJAhqFq1KkqUKIHq1aujRo0aqFGjBnr37o0tW7YgLi4uV/c+K4oQpnXp0iWUKFEC48ePVx67c+cOLC0ts92GIiMjA6tWrULnzp3RqlUrbNmypVCbmNPT09G2bVtYW1ujYsWKGDZsmMF7Nt24cQPVqlVD8+bNs93EMi9MUXt/6dIltG7dGjY2NgYvJuLi4rBx40ajat6f14G4adOm8PT0xN69e5Wgpx1xlpfb2QBPJ9Xr3r07fHx80KxZM0ycOLFAWiPySjVBJre9zLOemHSH/5lCbGxsrqu9C3v9uneU3bFjB2rVqoWgoCAA/9bMvPbaa8oVUW5GPLwo2z0vCrvs2oOWra0t/vOf/wB4ejBMTU3FrVu34O3trRzgjT1IPnnyBP7+/qhatSoyMzOVPlRmZmawsbHRGyqbVwVZfqBgr+6yhpnbt2+jb9++aNSoEb799ltlufz0Z9OlnU5+3759CA8PR4sWLZTbNVSoUAHVq1fH0KFDMXbs2Gd+N4UZwoB/v7cxY8bAyckJS5cuVfp/aPu19O3bFx9//DG2bt2q/BYeP36MI0eO4M8//yySJualS5dixowZOHLkCD777DNUrFgRbdq0wezZs/Wawm7cuAEnJycEBwebpA+OKVy9elUZbaY7221eatWM7UB84sQJTJ06FdbW1vnuT/WibE9DVBNkdD2vl3lue2UXF7o/CN0ws337dr0w88cff6B58+YICgrC8OHDc3UXa11q3u6FVXbdg9bevXuVx6dOnQoPD488XylmZmZi3759qFOnDl599VXlbujDhw+HlZWVMhQ2vyfpgiq/roK6utMNBREREbhz5w769+8PLy+vPA2dzUp322Y9qH/yySfw9PTEgwcPsHPnTsyZMweDBg1CnTp19DrLPq/cBR3CtCe7R48eYfDgwfD398c333yDmTNnwsnJCVOnTsXmzZvRvn17tGzZEo6OjmjXrh02bNhQpE3M58+fh6OjI7Zu3Qrg6XFu0aJFsLe3R9WqVTFnzhylWTQyMvK527ywaWs8g4ODDd4pPDeM7UCsbWZycHAwSafwnGZtfhGoMsjo0j0xGZqJsri7efMm3n//fb3ZGrPWzNSuXVsZmrdjxw40aNDguTO6Po+at3tBl133oHXixAnMnTsXNjY2+W5PzsjIwMGDB1GjRg0lzKSnp6Nnz56wt7fH/v37leVexPLrKqirO91QcO7cOdy6dQtDhgzJV8fYy5cvY/z48ejbty9Wrlyp95x2W+/btw/Vq1fP1qSU289Z0CEMeDojs7m5uXK/oMTERAwcOBDVq1eHra0tfvrpJ73ltX2sBgwYgI0bNxZqE/ONGzeU0KI1b948dO7cWQljb775Jry8vDB+/HgEBQXBzMzM5DPwmtLly5fxxhtvoFGjRkpHeWMY24E4PT0dISEhL0w/loKk+iAD6Pcy107R/rI4c+YMPDw8MGzYML0mIu0BNjk5GStXrkSdOnWUnf/nn3/O9/08AHVv94Iuu/agpZ16Py9zffz999/ZDnipqak4fPgwqlWrhgYNGijNTL169UKpUqWUmy2+COV/loK8urt8+TLatWsHf39/XLhwIV+h6dSpU3BxcUG7du3Qrl07aDQavVoSrRs3bsDW1hbbt2/Xe9yYz1YQIUzr1KlTKFmyJMaOHQvg3+NDYmIihg4dirp162LJkiXKMG/dbZaamlqozbS3b99Wmq90b9/wyy+/4JVXXkFMTAyGDh0KNzc3ZSr86OhobNiwwSQ3yyxIFy5cQPfu3fVu7WEMYzoQ9+3b94WrOSkoxSLIAKbpZa5WJ0+eRIMGDTB48GC9MKM9qDx8+BAuLi5YuHChyd9bzdu9oMt+8eJFdOzYMU9DTaOiolCmTBloNBq0aNECkydPxp49e5Sr2yNHjsDHxwfe3t5KzUzbtm1Rvnx5k908MD/lL2oXLlxAt27d8nzCAPTn1tH2FerXrx8mT56cbV6RxMRE+Pr6Ggw5xjBlCNM6c+YMbG1t9SasA6DMXaWtmfH398fChQsNhhldBd1M+/vvv8PMzAwNGzZEp06d9O6/1aNHD2g0GpQrV87g3drVwNg5abIqiA7Ealdsggzw4rXbFaYTJ04oYUb3qiQtLQ2PHj1CcHAwNm/eXCDvrebtXtBlz+sw2Zs3b8LHxwc1a9aEr68vBgwYABsbG/j4+KBfv35Yv349NmzYgBo1auD1118H8PS7NvUdZ00xzLeo5OeEERkZCScnJ/Tr10/v8S5duqBJkyaoV68e+vTpozcJXHBwMPr06QMgf/uVKUKYVlRUFMqWLYtOnTrpPf7pp5/igw8+UOaPSkxMREhICJo0aYLw8PBcb7uCaqZ9++234ePjg27duqFVq1ZKQDx69Chq1qypdCQuqinxi5opOxAXB8UqyLzstGFGt100NTUVYWFh8PDwMMmBkQrPlStX0KVLF3Tq1AmHDh1CZGQk1q5diyZNmsDPzw92dnaoV68eNBqNcoNRyr+kpCTs2LED9erVQ8+ePZU+JeHh4bC1tUVoaCjmzZuH8uXLw8fHR7lwmDdvnslq9/J71a51+vRpeHt7o0uXLvj5558B/Dvrq3b+EW3NS2JiIrp3747AwEDl3jy5YcpmWm1t0M8//4yBAwfil19+QdeuXdGsWTOlr1Djxo3x1ltv5et9igNTdCAuLhhkiplTp06hWbNmqF27Nrp06YLu3bujQoUKevdUIvW4ePEigoOD0bp1a2UyQ+DpBFffffcdpkyZgldeeeWl6NBXGLRX/Pfu3cOGDRvQqlUr9OrVCyNHjsw2t87Jkyeh0Wj0+nG8KOLj45W5oQ4fPozmzZuja9eu6N+/P5ydnbNNWKet2UhKSsrThHL5aaaNiorKNtT87t27qFWrFhYvXoy7d++ia9euaNq0KXbs2IH9+/fD0dExW+fkl1F+OxAXFwwyxVBkZCQWLlyIbt26YdasWarsv0L/0nYCDQ4ONtiZ90Ue9q4m2pludSeEW7NmDZo2bQorKytlVtT09HSkpaXh+vXrqFOnDrZt2wbgxWlivXPnDgIDA/H5558r/aUOHTqE5s2bw87OLtstK7Ty20yTl8+v2xesXbt2WL9+vXK8+vHHH9GsWTPcvXsX58+fR9euXdG6dWuMHTsWXbp0MclUAMVBfjsQFwcMMkQqoFuNrB1mTaajnWhMe9NM3ZPyr7/+iubNm6NLly7KKBkACA0NRfXq1U3eLym/njx5guDgYAQEBODLL79UwsyJEyfQokULdOjQQQlfQNEGsJs3b8LX1xcBAQFKH7/KlSvjiy++wPr16/HGG28oI8EiIiIQGBiIAQMGvHDbvKiZqilSrTQAIET0wrty5YqMGTNG4uLi5NNPP5VGjRoVdZGKhejoaGnQoIG0atVK1q9frzw+f/58iYuLk/DwcFm3bp18+eWX4uDgIPPnz5dNmzZJWFiYHDx4UHx8fIqu8FmkpaWJpaWlpKSkSL9+/eTmzZsydOhQefPNN8XW1laOHDkiEydOFDs7Oxk5cqS0bdu2qIssV65ckUmTJklmZqb0799fNBqNLFq0SEqVKiVbt24VPz8/2bt3r1hZWcmFCxfEwcFBKlSoUNTFphdJUScpIso9ViObXk43zSxZsqTeHYTXrVuHoKAglCtXDlZWViafWyc/DN29Pjk5Gd26dYOvr69ezczhw4cRGBiIpk2b6vX5KUoXL15E27ZtERQUhEuXLuHx48c4ePAg3njjDaxatQrAi9N0Ry8e1sgQqUxqaqpYWVkVdTGKlStXrsioUaPEyspKXF1dZevWrbJq1SoJCgqSzMxMMTMzExGR77//XlatWiXz58+XunXrFnGpn7pw4YI0atRI2rZtK1WrVpXBgweLjY2NuLu7S0pKigwePFjOnTsnw4cPlzfffFPs7Oxk//79MmfOHPnvf/8rFStWLOqPICJPv4ORI0eKiMi0adOkSZMmRVwiUgsGGSIiEbl8+bKMHDlS/vrrL5k5c6aMHTtWeU43zDx+/FhKlChRVMXUk56eLt99950MHjxY3NzcxMPDQ65fvy7Ozs7Spk0badOmjTRs2FAGDhwoycnJ0r17d+nbt6/Y2trKkydPxMbGpqg/gh5toAQgoaGh0rRp06IuEqkAgwwR0f+7du2ajBgxQszNzWXKlCnKiVR7mNRoNEVZPD3nzp2TH374Qd5//31ZsmSJfPjhh7JhwwZxdnaW8+fPy3fffSe3bt0SNzc3cXV1ld9++03c3Nxk5syZ8uabbwqAF+rzaLEvGBnLrKgLQET0oqhWrZosXrxYAMisWbNk//79IvI0wLxIJ/3Tp09LvXr1xMLCQhwcHGTSpEkyZswY6devn9y5c0eGDx8uv//+u5w+fVp69uwpderUESsrK0lKSlKCwYv0eXR5enrKvHnzpEKFCuLu7l7UxSEVYI0MEVEWL3KtwPnz58XX11cmTJgg06dP13tu4sSJ8umnn8rXX38tb731lt5z169fFwcHBylbtmwhljbv2BeMcotBhojIgIsXL8rUqVNl/vz5UqlSpaIujoiIRERESMuWLaVs2bJy/vx5Efl3yLXWpEmTZMGCBbJy5Urp06dPURWVqNAwyBAR5eBFqhU4ffq0NG7cWPz8/OTy5cvSvXt3WbRokYiIZGRkiLm5ubLspEmTZPHixfL5559LSEhIURWZqFCwjwwRUQ5elBBz7NgxadiwoUyYMEF2794tYWFhsmbNGhk9erSIiJibm0tGRoay/Jw5c2TgwIEyadIkSUhIKKpiExUK1sgQEb3g9u7dK5s2bVJqYOLj42X9+vXy4YcfSt++fXOsmbl79664uLgUSZmJCguDDBGRimiHTSckJMi6deuyhZn09HSxsLAo4lISFR7u7UREKqIdNl2yZEnp3bu3iIh8+OGHYm5uLgsWLGCIoZcO93giIpXShhkzMzMZOnSoWFtbS3h4eFEXi6hQMcgQEalYyZIlpUePHmJpaSkBAQFFXRyiQsc+MkRExcCLessBooLG4ddERMUAQwy9rBhkiIiISLUYZIiIiEi1GGSIiIhItRhkiIiISLUYZIiIiEi1GGSIiIhItRhkiIiISLUYZIiIiEi1GGSIiIhItRhkiIiISLUYZIiIiEi1/g+u0ePgGyTIzwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.cluster.hierarchy import dendrogram\n",
    "\n",
    "den = dendrogram(\n",
    "    clustermap.dendrogram_col.linkage,\n",
    "    labels=distance_matrix.index,\n",
    "    color_threshold=0.01,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "00c46f2e-04d9-4b79-bc0e-b1362d6af844",
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract clusters and perform gprofiler\n",
    "from collections import defaultdict\n",
    "\n",
    "\n",
    "def get_cluster_classes(den, label=\"ivl\"):\n",
    "    cluster_idxs = defaultdict(list)\n",
    "    for c, pi in zip(den[\"color_list\"], den[\"icoord\"]):\n",
    "        for leg in pi[1:3]:\n",
    "            i = (leg - 5.0) / 10.0\n",
    "            if abs(i - int(i)) < 1e-5:\n",
    "                cluster_idxs[c].append(int(i))\n",
    "    cluster_classes = {}\n",
    "    for c, l in cluster_idxs.items():  # noqa\n",
    "        i_l = [den[label][i] for i in l]\n",
    "        cluster_classes[c] = i_l\n",
    "    return cluster_classes\n",
    "\n",
    "\n",
    "clusters = get_cluster_classes(den)\n",
    "# extract functions for clusters\n",
    "cluster_process = {}\n",
    "for c in clusters:\n",
    "    cluster_df = pd.DataFrame(clusters[c], columns=[\"gene_symbol\"])\n",
    "    res = get_gprofiler(cluster_df)\n",
    "    cluster_process[c] = res[res[\"p_value\"] <= 0.05]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "ebfece22-89bd-4010-bcb9-e99431c63e94",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save grpofiler results\n",
    "with open(\n",
    "    \"figures/pert_embedding_cluster_gprofiler/pertspectra_inhouse_pert_emb_gprofiler.pickle\",\n",
    "    \"wb\",\n",
    ") as handle:\n",
    "    pickle.dump(cluster_process, handle, protocol=pickle.HIGHEST_PROTOCOL)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55eaf3f0-4b9f-4028-a042-f6f03996350f",
   "metadata": {},
   "source": [
    "## PR Curve+AUC\n",
    "- AUCPR using prior graph as binary label"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "ffd66c94-5f10-407c-b019-17d09a11a0ca",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "               AUC\n",
      "StringDB  0.675809\n",
      "CORUM     0.705489\n",
      "                F1  Threshold\n",
      "StringDB  0.697966   0.007078\n",
      "CORUM     0.639556   0.160069\n"
     ]
    }
   ],
   "source": [
    "aucpr_df, f1_df, pr_dict = auprc(distance_matrix)\n",
    "print(aucpr_df)\n",
    "print(f1_df)\n",
    "aucpr_df.to_csv(\"figures/pert_embedding_recall/pertspectra_inhouse_aucpr.csv\")\n",
    "f1_df.to_csv(\"figures/pert_embedding_recall/pertspectra_inhouse_f1.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "c3ef0a88-3b42-4d4c-b7f4-79c018ee3f45",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArMAAAIjCAYAAAAQgZNYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABqDUlEQVR4nO3deVxU9f7H8fcB2RVcAEWDcMnMJDUlf2ppluVWN7NulC1qaYta3axbapatYputlmWl1bUil7pWLpVm5XLT3FLT3FBccKEUFVCEOb8/aCYGZoDBgWFmXs/Hg0fNmXNmPsNxefvl8/1+DdM0TQEAAABeKMDTBQAAAACVRZgFAACA1yLMAgAAwGsRZgEAAOC1CLMAAADwWoRZAAAAeC3CLAAAALwWYRYAAABeizALAAAAr0WYBeA3Bg8erMTERJeuWbJkiQzD0JIlS6qkJm936aWX6tJLL7U93rVrlwzD0PTp0z1WEwD/QpgFUGWmT58uwzBsX6GhoWrZsqVGjhypgwcPerq8Gs8aDK1fAQEBql+/vvr06aMVK1Z4ujy3OHjwoB566CG1atVK4eHhioiIUIcOHfTMM8/o6NGjni4PgBeo5ekCAPi+p556Sk2bNtXJkye1dOlSvfXWW5o3b542btyo8PDwaqtj6tSpslgsLl3TrVs35eXlKTg4uIqqKt9NN92kvn37qrCwUFu3btWbb76pHj16aNWqVUpKSvJYXWdq1apV6tu3r06cOKFbbrlFHTp0kCT98ssvmjhxon788Ud98803Hq4SQE1HmAVQ5fr06aOOHTtKkoYOHaoGDRpo0qRJ+u9//6ubbrrJ4TU5OTmKiIhwax1BQUEuXxMQEKDQ0FC31uGqCy+8ULfccovt8SWXXKI+ffrorbfe0ptvvunByirv6NGjuvbaaxUYGKi1a9eqVatWds8/++yzmjp1qlveqyp+LQGoOWgzAFDtLrvsMklSenq6pKJe1tq1a2vHjh3q27ev6tSpo5tvvlmSZLFY9Morr+j8889XaGioGjZsqLvuuktHjhwp9brz589X9+7dVadOHUVGRio5OVkff/yx7XlHPbOffvqpOnToYLsmKSlJr776qu15Zz2zM2fOVIcOHRQWFqbo6Gjdcsst2rdvn9051s+1b98+9e/fX7Vr11ZMTIweeughFRYWVvr7d8kll0iSduzYYXf86NGj+te//qX4+HiFhISoRYsWeu6550qNRlssFr366qtKSkpSaGioYmJi1Lt3b/3yyy+2c6ZNm6bLLrtMsbGxCgkJUevWrfXWW29VuuaS3n77be3bt0+TJk0qFWQlqWHDhho3bpztsWEYeuKJJ0qdl5iYqMGDB9seW1tbfvjhBw0fPlyxsbE666yzNGvWLNtxR7UYhqGNGzfajm3ZskXXX3+96tevr9DQUHXs2FFz5849sw8NoEowMgug2llDWIMGDWzHCgoK1KtXL1188cV68cUXbe0Hd911l6ZPn64hQ4bovvvuU3p6ut544w2tXbtWy5Yts422Tp8+XbfffrvOP/98jRkzRnXr1tXatWu1YMECDRw40GEd3377rW666SZdfvnleu655yRJmzdv1rJly3T//fc7rd9aT3JyslJTU3Xw4EG9+uqrWrZsmdauXau6devazi0sLFSvXr3UqVMnvfjii/ruu+/00ksvqXnz5rrnnnsq9f3btWuXJKlevXq2Y7m5uerevbv27dunu+66SwkJCVq+fLnGjBmjzMxMvfLKK7Zz77jjDk2fPl19+vTR0KFDVVBQoJ9++kn/+9//bCPob731ls4//3z94x//UK1atfTll19q+PDhslgsGjFiRKXqLm7u3LkKCwvT9ddff8av5cjw4cMVExOjxx9/XDk5OerXr59q166tzz77TN27d7c7Ny0tTeeff77atGkjSdq0aZO6du2qJk2aaPTo0YqIiNBnn32m/v37a/bs2br22murpGYAlWQCQBWZNm2aKcn87rvvzMOHD5t79uwxP/30U7NBgwZmWFiYuXfvXtM0TXPQoEGmJHP06NF21//000+mJHPGjBl2xxcsWGB3/OjRo2adOnXMTp06mXl5eXbnWiwW2/8PGjTIPPvss22P77//fjMyMtIsKChw+hm+//57U5L5/fffm6Zpmvn5+WZsbKzZpk0bu/f66quvTEnm448/bvd+ksynnnrK7jXbt29vdujQwel7WqWnp5uSzCeffNI8fPiweeDAAfOnn34yk5OTTUnmzJkzbec+/fTTZkREhLl161a71xg9erQZGBhoZmRkmKZpmosXLzYlmffdd1+p9yv+vcrNzS31fK9evcxmzZrZHevevbvZvXv3UjVPmzatzM9Wr149s23btmWeU5wkc/z48aWOn3322eagQYNsj62/5i6++OJS9/Wmm24yY2Nj7Y5nZmaaAQEBdvfo8ssvN5OSksyTJ0/ajlksFrNLly7mOeecU+GaAVQP2gwAVLmePXsqJiZG8fHxuvHGG1W7dm19/vnnatKkid15JUcqZ86cqaioKF1xxRXKysqyfXXo0EG1a9fW999/L6lohPX48eMaPXp0qf5WwzCc1lW3bl3l5OTo22+/rfBn+eWXX3To0CENHz7c7r369eunVq1a6euvvy51zd133233+JJLLtHOnTsr/J7jx49XTEyMGjVqpEsuuUSbN2/WSy+9ZDeqOXPmTF1yySWqV6+e3feqZ8+eKiws1I8//ihJmj17tgzD0Pjx40u9T/HvVVhYmO3/s7OzlZWVpe7du2vnzp3Kzs6ucO3OHDt2THXq1Dnj13Fm2LBhCgwMtDuWkpKiQ4cO2bWMzJo1SxaLRSkpKZKkP//8U4sXL9YNN9yg48eP276Pf/zxh3r16qVt27aVaicB4Fm0GQCocpMnT1bLli1Vq1YtNWzYUOeee64CAuz/LV2rVi2dddZZdse2bdum7OxsxcbGOnzdQ4cOSfq7bcH6Y+KKGj58uD777DP16dNHTZo00ZVXXqkbbrhBvXv3dnrN7t27JUnnnntuqedatWqlpUuX2h2z9qQWV69ePbue38OHD9v10NauXVu1a9e2Pb7zzjv1z3/+UydPntTixYv12muvleq53bZtm3799ddS72VV/HvVuHFj1a9f3+lnlKRly5Zp/PjxWrFihXJzc+2ey87OVlRUVJnXlycyMlLHjx8/o9coS9OmTUsd6927t6KiopSWlqbLL79cUlGLQbt27dSyZUtJ0vbt22Waph577DE99thjDl/70KFDpf4hBsBzCLMAqtxFF11k68V0JiQkpFTAtVgsio2N1YwZMxxe4yy4VVRsbKzWrVunhQsXav78+Zo/f76mTZum2267TR988MEZvbZVydFBR5KTk20hWSoaiS0+2emcc85Rz549JUlXXXWVAgMDNXr0aPXo0cP2fbVYLLriiiv08MMPO3wPa1iriB07dujyyy9Xq1atNGnSJMXHxys4OFjz5s3Tyy+/7PLyZo60atVK69atU35+/hkte+ZsIl3xkWWrkJAQ9e/fX59//rnefPNNHTx4UMuWLdOECRNs51g/20MPPaRevXo5fO0WLVpUul4A7keYBVBjNW/eXN999526du3qMJwUP0+SNm7c6HLQCA4O1tVXX62rr75aFotFw4cP19tvv63HHnvM4WudffbZkqTff//dtiqD1e+//2573hUzZsxQXl6e7XGzZs3KPP/RRx/V1KlTNW7cOC1YsEBS0ffgxIkTttDrTPPmzbVw4UL9+eefTkdnv/zyS506dUpz585VQkKC7bi1rcMdrr76aq1YsUKzZ892ujxbcfXq1Su1iUJ+fr4yMzNdet+UlBR98MEHWrRokTZv3izTNG0tBtLf3/ugoKByv5cAagZ6ZgHUWDfccIMKCwv19NNPl3quoKDAFm6uvPJK1alTR6mpqTp58qTdeaZpOn39P/74w+5xQECALrjgAknSqVOnHF7TsWNHxcbGasqUKXbnzJ8/X5s3b1a/fv0q9NmK69q1q3r27Gn7Ki/M1q1bV3fddZcWLlyodevWSSr6Xq1YsUILFy4sdf7Ro0dVUFAgSbruuutkmqaefPLJUudZv1fW0eTi37vs7GxNmzbN5c/mzN133624uDg9+OCD2rp1a6nnDx06pGeeecb2uHnz5ra+X6t33nnH5SXOevbsqfr16ystLU1paWm66KKL7FoSYmNjdemll+rtt992GJQPHz7s0vsBqHqMzAKosbp376677rpLqampWrduna688koFBQVp27Ztmjlzpl599VVdf/31ioyM1Msvv6yhQ4cqOTlZAwcOVL169bR+/Xrl5uY6bRkYOnSo/vzzT1122WU666yztHv3br3++utq166dzjvvPIfXBAUF6bnnntOQIUPUvXt33XTTTbaluRITE/XAAw9U5bfE5v7779crr7yiiRMn6tNPP9W///1vzZ07V1dddZUGDx6sDh06KCcnRxs2bNCsWbO0a9cuRUdHq0ePHrr11lv12muvadu2berdu7csFot++ukn9ejRQyNHjtSVV15pG7G+6667dOLECU2dOlWxsbEuj4Q6U69ePX3++efq27ev2rVrZ7cD2Jo1a/TJJ5+oc+fOtvOHDh2qu+++W9ddd52uuOIKrV+/XgsXLlR0dLRL7xsUFKQBAwbo008/VU5Ojl588cVS50yePFkXX3yxkpKSNGzYMDVr1kwHDx7UihUrtHfvXq1fv/7MPjwA9/LkUgoAfJt1maRVq1aVed6gQYPMiIgIp8+/8847ZocOHcywsDCzTp06ZlJSkvnwww+b+/fvtztv7ty5ZpcuXcywsDAzMjLSvOiii8xPPvnE7n2KL801a9Ys88orrzRjY2PN4OBgMyEhwbzrrrvMzMxM2zkll+aySktLM9u3b2+GhISY9evXN2+++WbbUmPlfa7x48ebFfnj17rM1QsvvODw+cGDB5uBgYHm9u3bTdM0zePHj5tjxowxW7RoYQYHB5vR0dFmly5dzBdffNHMz8+3XVdQUGC+8MILZqtWrczg4GAzJibG7NOnj7l69Wq77+UFF1xghoaGmomJieZzzz1nvv/++6YkMz093XZeZZfmstq/f7/5wAMPmC1btjRDQ0PN8PBws0OHDuazzz5rZmdn284rLCw0H3nkETM6OtoMDw83e/XqZW7fvt3p0lxl/Zr79ttvTUmmYRjmnj17HJ6zY8cO87bbbjMbNWpkBgUFmU2aNDGvuuoqc9asWRX6XACqj2GaZfwMDgAAAKjB6JkFAACA1yLMAgAAwGsRZgEAAOC1CLMAAADwWoRZAAAAeC3CLAAAALyW322aYLFYtH//ftWpU0eGYXi6HAAAAJRgmqaOHz+uxo0bKyCg7LFXvwuz+/fvV3x8vKfLAAAAQDn27Nmjs846q8xz/C7M1qlTR1LRNycyMtLD1QAAAKCkY8eOKT4+3pbbyuJ3YdbaWhAZGUmYBQAAqMEq0hLKBDAAAAB4LcIsAAAAvBZhFgAAAF7L73pmAQCoDoWFhTp9+rSnywBqrKCgIAUGBp7x6xBmAQBwsxMnTmjv3r0yTdPTpQA1lmEYOuuss1S7du0zeh3CLAAAblRYWKi9e/cqPDxcMTExbNADOGCapg4fPqy9e/fqnHPOOaMRWsIsAABudPr0aZmmqZiYGIWFhXm6HKDGiomJ0a5du3T69OkzCrNMAAMAoAowIguUzV2/RwizAAAA8FqEWQAAAHgtwiwAAPAYwzD0xRdfuP1cb7dkyRIZhqGjR49KkqZPn666det6tKaaijALAAA0ePBgGYYhwzAUHBysFi1a6KmnnlJBQUGVvm9mZqb69Onj9nPPRGJiou17ER4erqSkJL377rtV/r7u8P3336tv375q0KCBwsPD1bp1az344IPat2+fp0urMoRZAABqqMzsPC3fkaXM7Lxqeb/evXsrMzNT27Zt04MPPqgnnnhCL7zwgsNz8/Pz3fKejRo1UkhIiNvPPVNPPfWUMjMztXHjRt1yyy0aNmyY5s+fXy3vXVlvv/22evbsqUaNGmn27Nn67bffNGXKFGVnZ+ull16q9Ou6615XFcIsAABVyDRN5eYXuPz10Ypd6jpxsQZO/VldJy7WRyt2ufwarm7aEBISokaNGunss8/WPffco549e2ru3LmSikZu+/fvr2effVaNGzfWueeeK0nas2ePbrjhBtWtW1f169fXNddco127dtm97vvvv6/zzz9fISEhiouL08iRI23PFW8dyM/P18iRIxUXF6fQ0FCdffbZSk1NdXiuJG3YsEGXXXaZwsLC1KBBA9155506ceKE7XlrzS+++KLi4uLUoEEDjRgxokI7s9WpU0eNGjVSs2bN9Mgjj6h+/fr69ttvbc8fPXpUQ4cOVUxMjCIjI3XZZZdp/fr1dq/x5ZdfKjk5WaGhoYqOjta1115re+6jjz5Sx44dbe8zcOBAHTp0qNy6nNm7d6/uu+8+3XfffXr//fd16aWXKjExUd26ddO7776rxx9/XJL0xBNPqF27dnbXvvLKK0pMTLQ9dnSvx44dq06dOpV637Zt2+qpp56yPX733Xd13nnnKTQ0VK1atdKbb75Z6c9UUR5dZ/bHH3/UCy+8oNWrVyszM1Off/65+vfvX+Y1S5Ys0ahRo7Rp0ybFx8dr3LhxGjx4cLXUCwCAq/JOF6r14wvP6DUspvTYfzfpsf9ucum6357qpfDgyv9VHxYWpj/++MP2eNGiRYqMjLSFutOnT6tXr17q3LmzfvrpJ9WqVUvPPPOMevfurV9//VXBwcF66623NGrUKE2cOFF9+vRRdna2li1b5vD9XnvtNc2dO1efffaZEhIStGfPHu3Zs8fhuTk5Obb3XrVqlQ4dOqShQ4dq5MiRmj59uu2877//XnFxcfr++++1fft2paSkqF27dho2bFiFvgcWi0Wff/65jhw5ouDgYNvxf/7znwoLC9P8+fMVFRWlt99+W5dffrm2bt2q+vXr6+uvv9a1116rRx99VB9++KHy8/M1b9482/WnT5/W008/rXPPPVeHDh3SqFGjNHjwYLtzXDFz5kzl5+fr4Ycfdvi8q/22Je+1JKWmpmrHjh1q3ry5JGnTpk369ddfNXv2bEnSjBkz9Pjjj+uNN95Q+/bttXbtWg0bNkwREREaNGhQpT5XRXg0zObk5Kht27a6/fbbNWDAgHLPT09PV79+/XT33XdrxowZWrRokYYOHaq4uDj16tWrGip2XWZ2ntKzctQ0OkJxUb6xeLb1M0UEByonv1BNoyMkqVLHij/nK98fAPB2pmlq0aJFWrhwoe69917b8YiICL377ru2UPef//xHFotF7777rm3N0GnTpqlu3bpasmSJrrzySj3zzDN68MEHdf/999teJzk52eH7ZmRk6JxzztHFF18swzB09tlnO63x448/1smTJ/Xhhx8qIqLo75c33nhDV199tZ577jk1bNhQklSvXj298cYbCgwMVKtWrdSvXz8tWrSo3DD7yCOPaNy4cTp16pQKCgpUv359DR06VJK0dOlSrVy5UocOHbK1Pbz44ov64osvNGvWLN1555169tlndeONN+rJJ5+0vWbbtm1t/3/77bfb/r9Zs2Z67bXXlJycrBMnTlRqe9dt27YpMjJScXFxLl/rSMl7LRXV//HHH+uxxx6TVBReO3XqpBYtWkiSxo8fr5deesmW6Zo2barffvtNb7/9tu+G2T59+rjUyD1lyhQ1bdrU1vdx3nnnaenSpXr55ZdrZJhNW5Wh0XM2yDQlQ9It/5egri2iPV3WGVm2PUv/+V+G3L3beIAhpQ5IUkpygptfGQA8KywoUL895drfUQeyT6rnpB9kKfaHbYAhfTequxpFhbr03q746quvVLt2bZ0+fVoWi0UDBw7UE088YXs+KSnJLtysX79e27dvV506dexe5+TJk9qxY4cOHTqk/fv36/LLL6/Q+w8ePFhXXHGFzj33XPXu3VtXXXWVrrzySofnbt68WW3btrUFWUnq2rWrLBaLfv/9d1uYPf/88+12l4qLi9OGDRskSRMmTNCECRNsz/32229KSCj6e+jf//63Bg8erMzMTP373//W8OHDbaFt/fr1OnHihBo0aGBXU15ennbs2CFJWrduXZmBefXq1XriiSe0fv16HTlyRBaLRVJRoG/dunWFvl/Fmabp1o06St5rSbr55pv1/vvv67HHHpNpmvrkk080atQoSUUDlDt27NAdd9xh97kLCgoUFRXltroc8artbFesWKGePXvaHevVq5f+9a9/Ob3m1KlTOnXqlO3xsWPHqqo8O5nZeRrzV5CVJFPSR//L0Ef/y6iW9/c2FlMaO2ejurWMYYQWgE8xDMPlH/U3i6mt1AFJGjtnowpNU4GGoQkD2qhZjOsjdq7o0aOH3nrrLQUHB6tx48aqVcu+7uLBUZJOnDihDh06aMaMGaVeKyYmRgEBrk3NufDCC5Wenq758+fru+++0w033KCePXtq1qxZrn+YvwQFBdk9NgzDFhzvvvtu3XDDDbbnGjdubPv/6OhotWjRQi1atNDMmTOVlJSkjh07qnXr1jpx4oTi4uK0ZMmSUu9n/XF+WVsZW1skevXqpRkzZigmJkYZGRnq1atXpSdbtWzZUtnZ2crMzCxzdDYgIKBUL7WjHuKS91qSbrrpJj3yyCNas2aN8vLytGfPHqWkpEiSrVd56tSppXprz2Sr2orwqjB74MAB27+0rBo2bKhjx44pLy/P4S+c1NRUuyH+6pKelWP3L2qrlg1rKzI0qPQTXuDYydPaevBE+SdWUqFpaldWLmEWACSlJCeoW8sY7crKVWJ0eLX82RgREWEbfayICy+8UGlpaYqNjVVkZKTDcxITE7Vo0SL16NGjQq8ZGRmplJQUpaSk6Prrr1fv3r31559/qn79+nbnnXfeeZo+fbpycnJswWvZsmUKCAiwTU4rT/369Uu9riPx8fFKSUnRmDFj9N///lcXXnihDhw4oFq1atlNnCruggsu0KJFizRkyJBSz23ZskV//PGHJk6cqPj4eEnSL7/8UqGanbn++us1evRoPf/883r55ZdLPX/06FHVrVtXMTExOnDggN1I7rp16yr0HmeddZa6d++uGTNmKC8vT1dccYViY2MlFeWxxo0ba+fOnbr55pvP6LO4yqvCbGWMGTPGNgQuFY3MWn/hVKWm0REKMGQXaAMNQx/cfpHXhrXM7Dx1nbjYYUh3h0DDUGJ0eNW8OAB4obiosBr9d8bNN9+sF154Qddcc42eeuopnXXWWdq9e7fmzJmjhx9+WGeddZaeeOIJ3X333YqNjVWfPn10/PhxLVu2zK4X12rSpEmKi4tT+/btFRAQoJkzZ6pRo0YOJy/dfPPNGj9+vAYNGqQnnnhChw8f1r333qtbb7211MCXO9x///1q06aNfvnlF/Xs2VOdO3dW//799fzzz6tly5bav3+/bdJXx44dNX78eF1++eVq3ry5brzxRhUUFGjevHl65JFHlJCQoODgYL3++uu6++67tXHjRj399NNnVF98fLxefvlljRw5UseOHdNtt92mxMRE7d27Vx9++KFq166tl156SZdeeqkOHz6s559/Xtdff70WLFig+fPnO/3HSEnW73t+fn6p0Pzkk0/qvvvuU1RUlHr37q1Tp07pl19+0ZEjR+yymLt51dJcjRo10sGDB+2OHTx4UJGRkU6H80NCQhQZGWn3VR3iosKUOiBJgX/9q8f6I6Ka/IdSeUp+JivDKOoJrswx23OS139/AMDfhIeH68cff1RCQoIGDBig8847T3fccYdOnjxp+/t20KBBeuWVV/Tmm2/q/PPP11VXXaVt27Y5fL06dero+eefV8eOHZWcnKxdu3Zp3rx5DtsVwsPDtXDhQv35559KTk7W9ddfr8svv1xvvPFGlXzW1q1b68orr9Tjjz8uwzA0b948devWTUOGDFHLli114403avfu3bYgfemll2rmzJmaO3eu2rVrp8suu0wrV66UVNSCMX36dM2cOVOtW7fWxIkT9eKLL55xjcOHD9c333yjffv26dprr1WrVq00dOhQRUZG6qGHHpJUNKL95ptvavLkyWrbtq1Wrlxpe64irr/+ev3xxx/Kzc0ttQLV0KFD9e6772ratGlKSkpS9+7dNX36dDVt2vSMP1tZDNPVReiqiGEY5S7N9cgjj2jevHm2xm1JGjhwoP78808tWLCgQu9z7NgxRUVFKTs7u1qCbWZ2XrX+iKg6WD9TeHCAcvMtttHUyhz7aMUuzdt4QEO6JGr8P8732GcCAHc5efKk0tPT1bRpU4WGVnyyFuBvyvq94kpe82ibwYkTJ7R9+3bb4/T0dK1bt07169dXQkKCxowZo3379unDDz+UVNSo/cYbb+jhhx/W7bffrsWLF+uzzz7T119/7amPUK6a/iOiynD2mSpzbMHGTElSnVCf73gBAABVwKNtBr/88ovat2+v9u3bS5JGjRql9u3b23apyMzMVEbG37P/mzZtqq+//lrffvut2rZtq5deeknvvvtujVyWCwAAAFXPo8Nhl156aZlb7RXfwaP4NWvXrq3CqgAAAOAtvGoCGAAAAFAcYRYAgCpQQ+ZXAzWWu36PEGYBAHAj625Hld3JCfAX1t8jZ7pDGFPIAQBwo1q1aik8PFyHDx9WUFCQy1u6Av7AYrHo8OHDCg8PL7VtsqsIswAAuJFhGIqLi1N6erp2797t6XKAGisgIEAJCQm2bXUrizALAICbBQcH65xzzqHVAChDcHCwW35yQZgFAKAKBAQEsAMYUA1o5AEAAIDXIswCAADAaxFmAQAA4LUIswAAAPBahFkAAAB4LcIsAAAAvBZhFgAAAF6LMAsAAACvRZgFAACA1yLMAgAAwGsRZgEAAOC1CLMAAADwWoRZeFTOqQJJ0vGTBR6uBAAAeCPCLDwmbVWGZq/ZJ0mavnyX0lZleLgiAADgbQiz8IjM7DyNmbNB5l+PTUlj52xUZnaeJ8sCAABehjALj0jPypHFtD9WaJralZXrmYIAAIBXIszCI5pGRyjAsD8WaBhKjA73TEEAAMArEWbhEXFRYUodkCRrnjUkTRjQRnFRYZ4sCwAAeBnCLDwmJTlB113YRJI0uEuiUpITPFwRAADwNoRZeFRESC1JUp3QWh6uBAAAeCPCLAAAALwWYRYAAABeizALAAAAr0WYBQAAgNcizAIAAMBrEWYBAADgtQizAAAA8FqEWQAAAHgtwiwAAAC8FmEWAAAAXoswC6+VmZ2n5TuylJmd5+lSAACAh9TydAGAqzKz8zRtabre+SldkhRgSKkDkpSSnODhygAAQHUjzMIrZGbnKT0rRxv2Zmvigi0yzb+fs5jS2Dkb1a1ljOKiwjxXJAAAqHaEWdR4aasyNGbOBllM5+cUmqZ2ZeUSZgEA8DOEWdRo6/cc0ejZG1RGjpUkBRqGEqPDq6UmAABQcxBmUePYtRTM31JukJWkCQPaMCoLAIAfIsyiRqlIS0FJ/+51LpO/AADwUyzNhRrD2lJQVpANkDSmbyt9Muz/1L1ltCQpKiyoegoEAAA1DiOzqBE27T+m1xdvL7OlIMCQPh/eRW3j60mSpi9Pr57iAABAjUWYhUflnCqQJC3acqjM8wINQxMGtLEFWQAAAIkwCw9KW5Wh2Wv2OX0+QNIjfVvpgiZ1lRgdzgQvAABQCmEWHpGZnacxc5wvuVWypQAAAMARJoDBI9KzcpxO9Ao0DKUOSCLIAgCAcjEyC49oGh2hAEN2gTbAkF67sb06JNajpQAAAFQII7PwiLioMKUOSFKgYUj6ezT2qraNCbIAAKDCGJmFx6QkJ6hbyxjtysplghcAAKgUwiw8Ki4qjBALAAAqjTYDAAAAeC3CLAAAALwWYRYAAABeizALAAAAr0WYhdc6ebpQkpSdd9rDlQAAAE8hzMIrpa3K0A9bsyRJLy78XWmrMjxcEQAA8ATCLLxOZnaexszZYHtsSho7Z6Mys/M8VxQAAPAIwiy8TnpWjt02uJJUaJralZXrmYIAAIDHEGbhdZpGRyjAsD8WaBhKjA73TEEAAMBjCLPwOnFRYUodkGR7bBjShAFt2EkMAAA/RJiFV0pJTlD3ltGSpIeuPFcpyQkerggAAHgCYRZeKzQoUJIUFRbk4UoAAICnEGYBAADgtQizAAAA8FqEWQAAAHgtwiwAAAC8FmEWAAAAXoswCwAAAK9FmAUAAIDXIswCAADAaxFm4VMys/O0fEeWMrPzPF0KAACoBh4Ps5MnT1ZiYqJCQ0PVqVMnrVy50um5p0+f1lNPPaXmzZsrNDRUbdu21YIFC6qxWtRkaasy1HXiYg2c+rO6TlystFUZni4JAABUMY+G2bS0NI0aNUrjx4/XmjVr1LZtW/Xq1UuHDh1yeP64ceP09ttv6/XXX9dvv/2mu+++W9dee63Wrl1bzZWjplm/54hGz94gi1n02GJKY+dsZIQWAAAf59EwO2nSJA0bNkxDhgxR69atNWXKFIWHh+v99993eP5HH32ksWPHqm/fvmrWrJnuuece9e3bVy+99FI1V46a5Jddf6r/m8tlljheaJralZXrkZoAAED18FiYzc/P1+rVq9WzZ8+/iwkIUM+ePbVixQqH15w6dUqhoaF2x8LCwrR06VKn73Pq1CkdO3bM7gu+4eTpQknSF+v2yyyZZCUFGoYSo8OruSoAAFCdPBZms7KyVFhYqIYNG9odb9iwoQ4cOODwml69emnSpEnatm2bLBaLvv32W82ZM0eZmZlO3yc1NVVRUVG2r/j4eLd+DnhG2qoM/bA1q8xzJgxoo7iosGqqCAAAeILHJ4C54tVXX9U555yjVq1aKTg4WCNHjtSQIUMUEOD8Y4wZM0bZ2dm2rz179lRjxagKmdl5GjNnQ5nn3HdZC6UkJ1RTRQAAwFM8Fmajo6MVGBiogwcP2h0/ePCgGjVq5PCamJgYffHFF8rJydHu3bu1ZcsW1a5dW82aNXP6PiEhIYqMjLT7gndLz8qxTfRyJqFBRPUUAwAAPMpjYTY4OFgdOnTQokWLbMcsFosWLVqkzp07l3ltaGiomjRpooKCAs2ePVvXXHNNVZeLGqRpdIQCDPtjAYb0xk3t9X/N6numKAAA4BEebTMYNWqUpk6dqg8++ECbN2/WPffco5ycHA0ZMkSSdNttt2nMmDG283/++WfNmTNHO3fu1E8//aTevXvLYrHo4Ycf9tRHgAfERYUpdUCSAo2iRBtoGEodkKSr2jZWaFCgh6sDAADVqZYn3zwlJUWHDx/W448/rgMHDqhdu3ZasGCBbVJYRkaGXT/syZMnNW7cOO3cuVO1a9dW37599dFHH6lu3boe+gTwlJTkBHVrGaNdWblKjA5nohcAAH7KME1Hixr5rmPHjikqKkrZ2dn0z/qgwdNWasnvh/XiP9vq+g5nebocAABQCa7kNa9azQAAAAAojjALAAAAr0WYhc/LzM7T8h1ZyszO83QpAADAzTw6AQyoammrMjRmzgZZzKLlu1IHJLGZAgAAPoSRWfis9XuOaPTsDbYNFiymNHbORkZoAQDwIYRZ+KSf0/9Q/zeXq+RSHYWmqV1ZuR6pCQAAuB9tBvBJs37ZWyrISkUbLCRGh1d7PQAAoGowMgufcvJ0oSQ5DLKGpAkD2rDBAgAAPoQwC5+RtipD/9v5p9PnH+zVkslfAAD4GMIsfEJmdp7GzNlQ5jmJDSKqqRoAAFBdCLPwCelZObZVC4p7rN95ap9Qt9rrAQAA1YMwC5/QNDpCAYb9sUDDUN8L4hRSi1/mAAD4Kv6Wh0+IiwpT6oAkBRpFiTbQMJjsBQCAH2BpLviMlOQEdWsZo11ZuUqMDifIAgDgBwiz8ClxUWGEWAAA/AhtBgAAAPBahFkAAAB4LcIsAAAAvBZhFgAAAF6LMAufd6rAIkk6kpvv4UoAAIC7EWbh09JWZWhtxlFJ0uNfbFLaqgzPFgQAANyKMAuflZmdpzFzNtgem5LGztmozOw8zxUFAADcijALn5WelSOLaX+s0DS1KyvXMwUBAAC3I8zCZzWNjlCAYX8s0DCUGB3umYIAAIDbEWbhs+KiwpQ6IMn2OMCQJgxoww5hAAD4EMIsfFpKcoLaJ9SVJD15zflKSU7wbEEAAMCtCLPweSG1in6Z1wsP9nAlAADA3QizAAAA8FqEWQAAAHgtwiwAAAC8FmEWAAAAXoswCwAAAK9FmAUAAIDXIswCAADAaxFmAQAA4LUIs/B5pwoskqQjufkergQAALgbYRY+LW1VhtZmHJUkPf7FJqWtyvBsQQAAwK0Is/BZmdl5GjNng+2xKWnsnI3KzM7zXFEAAMCtCLPwWelZObKY9scKTVO7snI9UxAAAHA7wix8VtPoCAUY9scCDUOJ0eGeKQgAALgdYRY+Ky4qTKkDkmyPAwxpwoA2iosK82BVAADAnQiz8GkpyQlqn1BXkvTkNecrJTnBswUBAAC3IszC54XUKvplXi882MOVAAAAdyPMAgAAwGsRZgEAAOC1CLMAAADwWoRZAAAAeC3CLAAAALwWYRYAAABeizALAAAAr0WYBQAAgNcizMLnnSqwSJKO5OZ7uBIAAOBuhFn4tLRVGVqbcVSS9PgXm5S2KsOzBQEAALcizMJnZWbnacycDbbHpqSxczYqMzvPc0UBAAC3IszCZ6Vn5chi2h8rNE3tysr1TEEAAMDtCLPwWU2jIxRg2B8LNAwlRod7piAAAOB2hFn4rLioMKUOSLI9DjCkCQPaKC4qzINVAQAAdyLMwqelJCeofUJdSdKT15yvlOQEzxYEAADcijALnxdSq+iXeb3wYA9XAgAA3I0wCwAAAK9FmIVfyszO0/IdWSzTBQCAl6vl6QKA6pa2KkNj5myQxSyaFJY6IIleWgAAvBQjs/Ar1o0UrOvPWkw2UgAAwJsRZuFXft75JxspAADgQ2gzgN9YvuMPffxzRqnjbKQAAID3YmQWfsNRkDXERgoAAHgzwix83qkCi9PnhnVrxuQvAAC8GGEWPi1tVYbWZhx1+nzzmNq2/2e5LgAAvA89s/BZ1pULylI/omhXMJbrAgDAOzEyC5+VnpVTauUCSXqs33lqHRdpe8xyXQAAeC/CLHxW0+gIBRj2xwINQ30viFNwrb9/6f+yi+W6AADwVoRZ+Ky4qDClDkhSoFGUaAMNo9TKBT9tO6z7PllX6lqW6wIAwDvQMwuflpKcoG4tY7QrK1eJ0eGlluD6aMVulexEYLkuAAC8B2EWPi8uKqxUMM3/a7kuBy21GtGjOZO/AADwErQZwO+krcrQb5nHnD5/TsM61VgNAAA4E4RZ+BVXlusCAAA1n8fD7OTJk5WYmKjQ0FB16tRJK1euLPP8V155Reeee67CwsIUHx+vBx54QCdPnqymauHtylquq0Vs7dJPAACAGs2jYTYtLU2jRo3S+PHjtWbNGrVt21a9evXSoUOHHJ7/8ccfa/To0Ro/frw2b96s9957T2lpaRo7dmw1Vw5vVdZyXUGBHv+3HQAAcJFH//aeNGmShg0bpiFDhqh169aaMmWKwsPD9f777zs8f/ny5eratasGDhyoxMREXXnllbrpppvKHc0FrCqyXBcAAPAeHlvNID8/X6tXr9aYMWNsxwICAtSzZ0+tWLHC4TVdunTRf/7zH61cuVIXXXSRdu7cqXnz5unWW291+j6nTp3SqVOnbI+PHXM+8Qf+obzlugAAgPfwWJjNyspSYWGhGjZsaHe8YcOG2rJli8NrBg4cqKysLF188cUyTVMFBQW6++67y2wzSE1N1ZNPPunW2uH9HC3XdbqwaLmuP3PyPVESAACoBK9qElyyZIkmTJigN998U2vWrNGcOXP09ddf6+mnn3Z6zZgxY5SdnW372rNnTzVWDG+RtipD2w+dkCT9K22d0lZleLgiAABQER4bmY2OjlZgYKAOHjxod/zgwYNq1KiRw2see+wx3XrrrRo6dKgkKSkpSTk5Obrzzjv16KOPKiCgdDYPCQlRSEiI+z8AfEbJ5bpMUxo7Z6O6tYyhBQEAgBrOYyOzwcHB6tChgxYtWmQ7ZrFYtGjRInXu3NnhNbm5uaUCa2BgoCTJNB3t5QSUz9FyXYWmqV1ZuZ4pCAAAVJhHt7MdNWqUBg0apI4dO+qiiy7SK6+8opycHA0ZMkSSdNttt6lJkyZKTU2VJF199dWaNGmS2rdvr06dOmn79u167LHHdPXVV9tCLeAq63JdxQNtoGEoMTrcc0UBAIAK8WiYTUlJ0eHDh/X444/rwIEDateunRYsWGCbFJaRkWE3Ejtu3DgZhqFx48Zp3759iomJ0dVXX61nn33WUx8BPsC6XNcjs4taDQIMsVwXAABewjD97Ofzx44dU1RUlLKzsxUZGenpclCD9Jz0g7YfOqFXb2yna9o18XQ5AAD4LVfymletZgBUJesOYPUjgj1cCQAAqCjCLAAAALwWYRYAAABeq1ITwAoLCzV9+nQtWrRIhw4dksVisXt+8eLFbikOAAAAKEulwuz999+v6dOnq1+/fmrTpo0Mw3B3XQAAAEC5KhVmP/30U3322Wfq27evu+sBAAAAKqxSPbPBwcFq0aKFu2sBPOp0YVG7zJ85+R6upGiL3eU7spSZnefpUgAAqNEqFWYffPBBvfrqq2whC5+RtipD2w+dkCT9K22d0lZleLSWrhMXa+DUn9V14mKP1gIAQE1XqTaDpUuX6vvvv9f8+fN1/vnnKygoyO75OXPmuKU4oDpkZudpzJwNtsemKY2ds1HdWsZU+y5g1lqsW+taPFgLAADeoFJhtm7durr22mvdXQvgEelZObbwaFVomtqVlVvhAJmZnaf0rBw1jY44o9D5y64/z7gWAAD8SaXC7LRp09xdB+AxTaMjFGDILkQGGoYSo8MrdH3aqgzbaGqAIaUOSFJKcoLLdaStytDo2RtKHbfW4q7ADACALzmjTRMOHz6spUuXaunSpTp8+LC7agKqVVxUmFIHJNkeBxjShAFtKhQYnbUFuDpxKzM7T6Nnb1DJLnRDRbX8uPVwmX20TBgDAPirSoXZnJwc3X777YqLi1O3bt3UrVs3NW7cWHfccYdyc3PdXSNQ5VKSE9QitrYk6eWUdhUeWS2rLaCiMrPz9NI3W0sFWUmadENbdWsZU2ZgZsIYAMCfVSrMjho1Sj/88IO+/PJLHT16VEePHtV///tf/fDDD3rwwQfdXSNQLYICi3471I8IrtD5aasydN8n60odd7VFocvExZq1eq/D55vUC9OX6/Y7DczuGhkGAMBbVapndvbs2Zo1a5YuvfRS27G+ffsqLCxMN9xwg9566y131QfUSNYQ6awtoKItCqPnbFBZK9zd+M7/SgVZ6e/AzIQxAIC/q9TIbG5urho2bFjqeGxsLG0G8HmZ2XmatXqvw5A5tu95FW5R+Oa3gw6D7MO9zrX9v7Mga+2jve/TdQ6fr+jIMAAA3q5SYbZz584aP368Tp48aTuWl5enJ598Up07d3ZbcUBNY+1PfembrQ6fb904skKTsaYtS9f4/24qdTzQMHTl+Y3KrGHO8M6KCKmlR2Y7HtWt6MgwAAC+oFJtBq+++qp69eqls846S23btpUkrV+/XqGhoVq4cKFbCwRqipL9qY78svtP3frez06X6crMztP/dv6hJ7/8rdS1ASoKosu2Z5VZx7VvLndaw4D2jSu1LBgAAN6qUmG2TZs22rZtm2bMmKEtW7ZIkm666SbdfPPNCgtjRAje6XShRZL0Z06+w+cdba4gSY/1O0+vL96uo3mn9cp322yjpSV37yq+Hq0jrw9srwvPrqeuExeXes6QbP25ZYXps+rRXgAA8C+VCrOSFB4ermHDhrmzFsBj0lZlaPuhE5Kkf6Wt08nThaVGOAsKLKWuCzQM9b0gTpOX7JCkUj/2L75MV1lBNtAwdOHZ9ZwG5gkD2mjMnI1O67eG3ciwIKfnAADgiyocZufOnas+ffooKChIc+fOLfPcf/zjH2dcGFBdrO0DVmaJEVVJ+mRlht050t8TseKiwmRxklLLWnXA0etIcrgbWebRkw6vDTCk125sr1mr92rJ1sM6lne6oh8bAACfUOEw279/fx04cECxsbHq37+/0/MMw1BhYaE7agOqhaPR0OLLW2Vm52lsiSAbYBRNxGobX09pqzJ01EGILL7qgKNtaq1BtENiPVuQte5GNnbORhWapgINQw/3OVfPzd9S+noV9eTm5BdoydaiHfheX7xdTeqF0TcLAPAbFQ6zFovF4f8D3q5pdITD0dDE6HBlZufp7R92llpP1mJKufmWUqO6Vh3Orqs3Bl4oSeo6cXGp661B9Kq2jUtdm5KcoG4tY7QrK1eJ0eFOWw8c9diaKj2qDACAL6t0z2xJR48eVd26dd31ckC1sY6GPvLX6GmAIduIqrM+V2vYdRY0Q2oFSpI+X+N4PdrXB7ZXvwtKB9niNRUPo47CtrMeWzZNAAD4k0qtM/vcc88pLS3N9vif//yn6tevryZNmmj9+vVuKw6oLinJCWoRW1uS9HJKO3VrGVNmkLX2uFpHdUs6kpOvrhMX6/mFpdejtQbRirKG7UDDqND7u2PThIqslQsAQE1QqZHZKVOmaMaMGZKkb7/9Vt99950WLFigzz77TP/+97/1zTffuLVIoDoEBRb9265+RHCZy3D1vSDOaY+rYRRNINt84LjD9yg52auiSrYelHz/0bOLttY1jDPfNKH4EmKO1soFAKAmqVSYPXDggOLj4yVJX331lW644QZdeeWVSkxMVKdOndxaIOAJefkFpY5Zl+EqGRSLB82fd/6hVxZtc/iaJYOwq0q2HhR//++3HNaCTQd0b48WZxQ81+85YgvGUum1cmuyzOw8pWflqGl0RI2vFQDgPpUKs/Xq1dOePXsUHx+vBQsW6JlnnpEkmabJSgbwet//fljvL023O1beiKo1aG7OPObweWdB2F3Cgot6dM9kndm0VRl2QdaqJvfgWgPshr3Zmrhgi0xGkwHA71QqzA4YMEADBw7UOeecoz/++EN9+vSRJK1du1YtWrRwa4FAdSsZZIsvw1WetRlHSh2rbGuBK/Lyi/4RWZl1ZjOz87R69xHbBLiSAgwp68RJZWbn1ahA62xHNW8aTQYAnLlKhdmXX35ZiYmJ2rNnj55//nnVrl00cSYzM1PDhw93a4FAdbFuZ1uSdRmu8mRm5+mrXzPtjrkShCsrbVWGFmw6IMn1dWbL22JXKuoBvveTddU+4llW20DJdoiSavJoMgDAvSoVZoOCgvTQQw+VOv7AAw+ccUGAJxTfzrakiq4OkJ6V43Q92qpSavcyVXxU0nptWUHW+ppS9Y54ljUJLW1VhkbPcR5kJfes6AAA8A5sZwu/52zjA8m1FoGyNl+oKmeyzmxZW+w6Y31t63u7e7KVo5aH4iH60LGTZY7ISn+vE8yoLAD4B7azhd+r6DJc5XG0FW1Vh6rKBmjrZK+SAgzpxovi9fHPexxeF2gY+nXvUd387v/cvnRXWS0PhaapacvSNfXH9DKDrFTUFgEA8B+GafrXH/3Hjh1TVFSUsrOzFRkZ6elyUANkZuep68TFpQLh0tE9KhVEM7PzSq0HW5WKr0JgGNLoPq2U1CTK4aipdeRz5MdrS71OgKRH+rbSc/O3OAyURUE3QR//nGF3/Ey+V1br9xxR/8nLnQbVkoG9vOfcURMAwHNcyWuV2gEM8CVl7bBV2dfr3LxBtQWplOQE9Tq/kSSpx7mxem7+Fg2c+rO6TlystFV/B8+0VRnqOnGxwyArFW2xm9QkymFojAytJdNUqSAr2bceVEbaqowyg6xURpCVdMfFTd1eEwDAe1RqAth9992nFi1a6L777rM7/sYbb2j79u165ZVX3FEbUG2c7bDlLazrzH6/5ZDDCVuSypzsVXyLXUcjncdOlt5EwirAkMKDHf+7uKwVCcpbEqw8AYb0+fAuio0M1XtL06u1VxkAUHNUamR29uzZ6tq1a6njXbp00axZs864KMATqntE1Z2s68w62/CgrMlexUeiS45SGxV4b4spXfvmcrtRYOnvkeDKjBKXJ9AwlDogSW3j6ykuKkzj+rUu9hwTwADAn1RqZPaPP/5QVFRUqeORkZHKyso646IAVFzxdWZLCjQM/brnqCYu2FLquQBDeu3G9uqQWM8u+FlHqdfsPqIRFQybJUeBz2RFgvL6Yx3V3L99Ez311W+SpNn3dFa7hPoVqhsA4P0qNTLbokULLViwoNTx+fPnq1mzZmdcFICKKWtZsQBDSrkoXqkLtpQKjgEqWoXgqraNHY5gxkWFqV5EsMPXdTZaa11xwNmIq/X5svpjAw3DYQ9seTV/sXaf7f8HvLWi1CgxAMB3VWpkdtSoURo5cqQOHz6syy67TJK0aNEivfTSS/TLAtXI2bJitQIMFVpMhxO2pKLJXv0uaFzmazta9ivAkK48r6EW/Haw1PkBhspcOivAkN75Md3pc9YRV0mlemCt/bGOdlLLzM7TM1//ZnvMdrYA4F8qNTJ7++2366WXXtJ7772nHj16qEePHvrPf/6jt956S8OGDXN3jQCcsAbOkgosZpmjn9bJXmVxtMrDI31aaeHm0kHWkJTSMb7SKxIUH3F19L7W/lhHyto4AgDg+854ndnDhw8rLCxMtWvXdldNVYp1ZuFr0lZl2DZqKKvfVPp7spcrmxwUXzc3PStHA6f+XOqcHufG6PvfD7tce3kjrhVZXcLd6wQDADzPlbxWqTYDSSooKNCSJUu0Y8cODRw4UJK0f/9+RUZGek2wBXxB8WXFvt9yUO/8VPpH+c4mTlWEdbS0+GuVDMyVCbLWYO1sxLXk+5ZV37h+rW0TwNjOFgD8S6XC7O7du9W7d29lZGTo1KlTuuKKK1SnTh0999xzOnXqlKZMmeLuOgGUwRrc3l3qIMjq7x/ju+N9im/ZWxlnEqwrwr/2NAQAVKpn9v7771fHjh115MgRhYX9/ZfRtddeq0WLFrmtOAAV52wy2OsD27vUVlCelOQEzRneucw1aB318Urlr6JQGSUngJkqmgCWmZ3nltcHANRslQqzP/30k8aNG6fgYPulexITE7Vv3z4nVwGoSo4mg1V0spercvILy5xg9kifVqVqCTCkz0d0cWuwlpgABgD+rlJtBhaLRYWFhaWO7927V3Xq1DnjogC4rmQLQPGdvdzN2bJdxdsH6oYFlarFWX+su2thO1sA8B+VWs0gJSVFUVFReuedd1SnTh39+uuviomJ0TXXXKOEhARNmzatKmp1C1YzgK+r6CoAZ6r4KgrOVkmorlreX5pumwBWtJ1tkttHgAEA1ceVvFapMLtnzx717t1bpmlq27Zt6tixo7Zt26bo6Gj9+OOPio2NrXTxVY0wC7hPdYXV8vyZk68Ln/5WkvTFcLazBQBvV+VhVipamistLU3r16/XiRMndOGFF+rmm2+2mxBWExFmAd9TfGQ2wCiaZMbILAB4ryoNs6dPn1arVq301Vdf6bzzzjujQj2BMAv4FjZNAADf40pec3k1g6CgIJ08ebLSxQGAO7GaAQD4t0otzTVixAg999xzKigocHc9AOASZ0uSsZoBAPiHSi3NtWrVKi1atEjffPONkpKSFBERYff8nDlz3FIcAJSH7WwBwL9VKszWrVtX1113nbtrAYAzxna2AOBfXAqzFotFL7zwgrZu3ar8/HxddtlleuKJJ2r8CgYAfJez7Wy7tYxhdBYA/IBLPbPPPvusxo4dq9q1a6tJkyZ67bXXNGLEiKqqDQDKxQQwAPBvLoXZDz/8UG+++aYWLlyoL774Ql9++aVmzJghi8VSVfUBQJmYAAYA/s2lMJuRkaG+ffvaHvfs2VOGYWj//v1uLwwAKsI6AcyKCWAA4F9cCrMFBQUKDQ21OxYUFKTTp0+7tSgAqCwmgAGAf3FpAphpmho8eLBCQkJsx06ePKm7777bbnkuluYCUF2YAAYA/s2lMDto0KBSx2655Ra3FQMAriprAhhhFgB8n0thdtq0aVVVBwBUinUCWPFAywQwAPAfldrOFgBqCiaAAYB/I8wC8ClMAAMA/0KYBeDVnE0Ay8zO81xRAIBqQ5gF4NWcTQBbs/uIZwoCAFQrwiwAr+ZoBzBJGvnxWqWtyqj+ggAA1YowC8CrlZwAZkW7AQD4B8IsAJ9lXW8WAOC7CLMAvFrJCWDFsd4sAPg+wiwAr+ZoApgV680CgO8jzALwas4mgAEA/ANhFoBXi4sKU+qAJIeBlglgAOD7CLMAvF5KcoJevbFdqeNMAAMA31cjwuzkyZOVmJio0NBQderUSStXrnR67qWXXirDMEp99evXrxorBlDTdEysX2p0lglgAOD7PB5m09LSNGrUKI0fP15r1qxR27Zt1atXLx06dMjh+XPmzFFmZqbta+PGjQoMDNQ///nPaq4cQE1ibTew5llDTAADAH/g8TA7adIkDRs2TEOGDFHr1q01ZcoUhYeH6/3333d4fv369dWoUSPb17fffqvw8HDCLABJRZslFP8vAMC3eTTM5ufna/Xq1erZs6ftWEBAgHr27KkVK1ZU6DXee+893XjjjYqIiHD4/KlTp3Ts2DG7LwC+JzM7T2PmbLA7xgQwAPB9Hg2zWVlZKiwsVMOGDe2ON2zYUAcOHCj3+pUrV2rjxo0aOnSo03NSU1MVFRVl+4qPjz/jugHUPI7Wm2UCGAD4Po+3GZyJ9957T0lJSbroooucnjNmzBhlZ2fbvvbs2VONFQKoLo7Wm2UCGAD4Po+G2ejoaAUGBurgwYN2xw8ePKhGjRqVeW1OTo4+/fRT3XHHHWWeFxISosjISLsvAL4nLipM17ZvYnesf/vGTAADAB/n0TAbHBysDh06aNGiRbZjFotFixYtUufOncu8dubMmTp16pRuueWWqi4TgBfIzM7T52v32R37Yu1+emYBwMd5vM1g1KhRmjp1qj744ANt3rxZ99xzj3JycjRkyBBJ0m233aYxY8aUuu69995T//791aBBg+ouGUANRM8sAPinWp4uICUlRYcPH9bjjz+uAwcOqF27dlqwYIFtUlhGRoYCAuwz9++//66lS5fqm2++8UTJAGoga89s8UBLzywA+D7DNE2/Wo7x2LFjioqKUnZ2Nv2zgI9JW5Wh0bM3yFTRpgkTr0tSSnKCp8sCALjIlbzm8TYDAHAnNk0AAP9CmAXgE9g0AQD8E2EWgE9gAhgA+CfCLACfwKYJAOCfCLMAfEJcVJhSByTJmmcNSRMGtGHTBADwcYRZAD6FCWAA4F8IswB8AhPAAMA/EWYB+AQmgAGAfyLMAvAJTAADAP9EmAXgE+KiwnRt+yZ2x/q3b8wEMADwcYRZAD4hMztPn6/dZ3fsi7X76ZkFAB9HmAXgE+iZBQD/RJgF4BPomQUA/0SYBeAT2DQBAPwTYRaAT2HTBADwL4RZAD6BTRMAwD8RZgH4BCaAAYB/IswC8AnumACWmZ2n5TuyGM0FAC9Sy9MFAIA7WDdNmL3m77VmXdk0IW1VhsbM2SCLKQUYUuqAJKUkJ1RVuQAAN2FkFoBPcLRpwudr9+nL9fvKHWldv+eIRs/eYGtTsJj02wKAtyDMAvAJjnpmLaZ07yfr1HXiYqWtynB4XdqqDPV/c3mp1Q/otwUA70CYBeATHPXMWllMacycDVq/54jtWGZ2nr76db9Gz94g08E6Xmy4AADegZ5ZAD7BUc9scRZT6j95uSZelyRJtv5YZ1zptwUAeA5hFoBPcNQzW5Kp8kOs1Rdr9+uhXucSaAGghqPNAIBPcNQz60hFzpHomQUAb0GYBeATyuqZLY+j6+iZBQDvQJgF4BPiosKUOiBJgYZriTZARWvKXndhE7vj9MwCgHcgzALwGSnJCVo6uocmD2yvikTaAEP6fEQXdWsZU6rf9ou1+1lnFgC8AGEWgE+JiwpTvwsaa+J1SWW2HQQahlIHJKltfD2H/bb0zAKAd2A1AwA+KSU5Qa0a1VH/yfYbIgQY0ms3tleHxHq2NgJrv23xQEvPLAB4B0ZmAfistvH1NPG6v/toraOxV7W174e1rlFbHD2zAOAdGJkF4NNSkhPUrWWMdmXlKjE63GFAdbRGLevMAoB3IMwC8HlxUWFlhtKyemYJswBQs9FmAMDvOVqjlp5ZAPAOhFkAfo+eWQDwXoRZAH7PWc8s68wCQM1HmAXg91hnFgC8F2EWgN+jZxYAvBdhFoDfo2cWALwXYRaA36NnFgC8F2EWgN+jZxYAvBdhFoDfo2cWALwXYRaA36NnFgC8F2EWgN+jZxYAvBdhFoDfo2cWALwXYRaA32saHaESLbMyDNEzCwBegDALAI6Y5Z8CAPA8wiwAv5eelVMqu5oSbQYA4AUIswD8HktzAYD3IswC8HsszQUA3oswC8DvsTQXAHgvwiwAv8fSXADgvQizAPwePbMA4L0IswD8nqOe2V5tGio9K4dWAwCo4Wp5ugAA8DRHPbPzNhzQvA0HFGBIqQOSlJKc4KHqAABlYWQWgN9z1DNrZTGlsXM2MkILADUUYRaA33O0nW1xTAYDgJqLMAsA5WAyGADUXIRZAH7P0Xa2xbGBAgDUXIRZAH7P0dJcxbGBAgDUXIRZAH4vLipMqQOSFGg4TrT0zAJAzcXSXAAgKSU5Qd1axmjN7iMa8fFau+cMQ/TMAkANRZgFgL/ERYXpwrMdPFFWQy0AwKNoMwCAYtKzckodMyXaDACghiLMAkAxTaMjVLJ1lqW5AKDmIswCQDFxUWG6tGWM3TGW5gKAmoswCwDFZGbnacnWw3bHWJoLAGouwiwAFJOelSOzxIQvluYCgJqLMAsAxTSNjlDJ1WZZmgsAai7CLACUUGolrr8OZGbnafmOLFoOAKAGYZ1ZACjG2dJc05al692f0mUxpQBDSh2QpJTkhOovEABgh5FZACjGYZuBpKk/FgVZSbKY0tg5GxmhBYAagDALACWUbDMwHRxjUhgA1AyEWQAoxlGbgSNspAAANQNhFgCKcbQDmCNspAAANQNhFgCKcbQDmCNspAAANYPHw+zkyZOVmJio0NBQderUSStXrizz/KNHj2rEiBGKi4tTSEiIWrZsqXnz5lVTtQB8naMdwBwp2TPLsl0A4BkeXZorLS1No0aN0pQpU9SpUye98sor6tWrl37//XfFxsaWOj8/P19XXHGFYmNjNWvWLDVp0kS7d+9W3bp1q794AD7J0Q5gjhTfSCFtVYbGzNlgW7brkT6tlNQkSk2jI2hFAIAqZphmRf7YrhqdOnVScnKy3njjDUmSxWJRfHy87r33Xo0ePbrU+VOmTNELL7ygLVu2KCgoqFLveezYMUVFRSk7O1uRkZFnVD8A35OZnaeuExfbluGSigKqpcSflIak5WMu06FjJ9V/8vLSGy38dd0dFzfV7Rc3JdQCgAtcyWseazPIz8/X6tWr1bNnz7+LCQhQz549tWLFCofXzJ07V507d9aIESPUsGFDtWnTRhMmTFBhYaHT9zl16pSOHTtm9wUAzsRFhSl1QJIC/5oFFmgYuuPipqXOs26k4CzISkUBeOpP6eqSulhpqzKqrmgA8GMeazPIyspSYWGhGjZsaHe8YcOG2rJli8Nrdu7cqcWLF+vmm2/WvHnztH37dg0fPlynT5/W+PHjHV6TmpqqJ5980u31A/BdKckJ6tYyRruycm2tBO/+lG4XWg1J7/yYXqHXM1W0yUK3ljGM0AKAm3l8ApgrLBaLYmNj9c4776hDhw5KSUnRo48+qilTpji9ZsyYMcrOzrZ97dmzpxorBuCt4qLC1Ll5A1v4dLSRgiscbbJwJpPGmHAGAEU8NjIbHR2twMBAHTx40O74wYMH1ahRI4fXxMXFKSgoSIGBgbZj5513ng4cOKD8/HwFBweXuiYkJEQhISHuLR6AX6noRgplKT5hTCo9aSx1QJJSkhMcXpuZnaf0rBxFBAcqJ79Qv+w6ope/3SpT5V8LAL7OY2E2ODhYHTp00KJFi9S/f39JRSOvixYt0siRIx1e07VrV3388ceyWCwKCCgaVN66davi4uIcBlkAcIem0REy5PporJ2/Ls7MztPPO//QI7M32J6ymM7bED5asUuP/3dTmX25tDAA8GcebTMYNWqUpk6dqg8++ECbN2/WPffco5ycHA0ZMkSSdNttt2nMmDG28++55x79+eefuv/++7V161Z9/fXXmjBhgkaMGOGpjwDATzgLkwGGNLBTfLl/mFonjHWduFj/Sltf6nlH69a+t3SnHisjyDq7FgD8iUfXmU1JSdHhw4f1+OOP68CBA2rXrp0WLFhgmxSWkZFhG4GVpPj4eC1cuFAPPPCALrjgAjVp0kT333+/HnnkEU99BAB+wFmbQYD+/hH/vZedo11Zufp171E9t2CLw6W8ypowZm1DyMzO07Sl6Xrnp4pNLrMKD/aqKRAA4DYeXWfWE1hnFoCrnK09+/nwLmobX6/U+ev3HNE1k5e79B6GpNF9WjkMwhXxybD/U+fmDVy/EABqIK9YZxYAvIWjtWdTByQ5DLKSlJPvfO1rZ0xJqfMrF2RLTi4DAH/i0TYDAPAWJdeeLWuylVsmjDnhaDeyKnkjAPASjMwCQAWVXHu2sgIM1863nl7WbmRrdh85o5oAwFsxMgsAbpaeleN0sDTQMHT7xYmaWs4ErwBDGnpxMw25OFGS7HYjc3TtyI/X6sSpAtabBeB3CLMA4GZNoyNKtQMEGNJrN7ZXh8R6OnTspNMwW/y84iPA1v93tuMXW+YC8Fe0GQCAmzmbMHZV28aKiwpzOkHMutSX9TxHytqNjPVmAfgjRmYBoAqUNWHM2cits6W+iosIDnT6HKsaAPBHjMwCQBVxNmHM1aW+iitz2S9WNQDghxiZBQAPcGWpr+IcjepamSqaKEbPLAB/wsgsAHhIZZb6so7qOlreizYDAP6IMAsAXiYlOUGfD+9S+gnaDAD4IcIsAHghR72zpqSvf810unwXAPgiwiwAeCFnqxo88/VmdZ24WGmrMqq5IgDwDMIsAHihslY1sJhFGyg4GqHNzM7T8h1ZjN4C8BmsZgAAXqhpdIQMOW+TtW6gUHxyWdqqDI2Zs0EWs2hd29QBSWx/C8DrMTILAD4qPPjvP+LX7zmi0bM32Jb0Kmv0FgC8CWEWALxQelZOuYsX7D1SFFTTVmWo/5vLS53P9rcAfAFtBgDghcraPMFq5MdrtWzHH/rk5wyHwde6Lm1mdp7Ss3LUNDqCDRcAeB3CLAB4IevmCWPnbFSh6TjRmpI+/rmMVQ1Mae66/XpuwRb6aAF4LcM0nfwp6KOOHTumqKgoZWdnKzIy0tPlAMAZyczO09e/ZuqZrze75fUCDUNLR/dghBaAR7mS1+iZBQAvFhcVpuTEem57PfpoAXgbwiwAeLmy1py1MoyKv17xVRAAoKbjTywA8HLWyWDOBEga3aeVKppnc/Mt7igLAKoFYRYAvJx1MpijQBtgSJ+P6KJ/tG1codeyrnBQGdWxu5ir78GOZ4DvYzUDAPABKckJ6tYyRtOWpevdH9NlUdFkrgkD2qhtfD0t35FV7rq0kpxvKVYOZ7uLncmyXyWvdXUHs5LnP9KnlZKaRCkiOFA5+YUsRQb4CFYzAAAfk5mdp11ZuUqMDreFtczsPHWduNhuXVrDkBz9DfDJsP9TYnR4hUJoZnaeft75h/6Vtt7ueIAh3dO9ud76YUellv0qGURvvCih1DJjZa28sH7PEfWfXHqjiJI1shQZUDO5ktcIswDgJ9JWZdjWpQ00DN1zaTO98f2OUueNuLT8EJqZnadpS9M19af0Cg/mVmTZr8zsPK3efUT3fry2Qq87eWB7XXh2PaVn5dhGXDfszdbE+VsqdD1LkQE1E2G2DIRZAP6s+KhtelaOBk79udxrSga+4qOmrvpk2P+pc/MGDp+r7OsaqnR3hKSiQNzvgor1FAOoHq7kNXpmAcCPxEWF2ULpoWMnK3RN8bVnV+8+otGzN1Q6PObmn3Z4fP2eyr/umY7IjPh4rU6cKqDdAPBShFkA8FMVWZ/Waum2w7r53R2VGo0tbugHqzW6r/1ErHUZR/X8wt/P7IXP0Og5G9StZQztBoAXIswCgJ+KCA6s8LmTl5Tura0MU1LqvC1ueS13Mk1pze4j6ncBYbYmOZPVMOA/CLMA4KdcGZl1JsCQbvprpYHig7bOVkpwhbUX1tl7OK1J0iN9W+mCJnUVHhyg3HyLft17VM8t2FLmyLInZ5C4O7T5Qgh0dSk2+C/CLAD4KevOYZVpHQgwpNdubK8OifUUFxWmC86KqtBKCRV97c+Hd1FsZKjdEmMXnBVV7gQx67Vt4+vZHe/cvIH+r1n9Mpfriq/vmdDnbD3cyq7NO21put75KV3Smb9eVXMUuq0rWjwye4PtPIspjZ2zkVYQOMRqBgDgx0ou13Vr5wRNX767zGsCJKVe53i5LldXSijJutGDsxG4zOw8u40hjL+Gb80KXCuVvWJCWSstVJWy1sOt6GikNRBu2JutiQu2OB1hrmmjm45CfNbxU5r6VxB3xBP3CJ7BagYAgAqx7hxWPISWFWadjXxK9islWM+tyKhv8baA4hs9OBIXFaaxfVtrSNemtpolldokwpmU5ARF1w7WHR+sLvVceHD17PBuFz7LWA/XYpaemFbWrmjlsZjSmDkb1KpRHYf3rzqVXL3CYlasl7q67hG8C2EWAPxcRUNo8e1xK/KaqQOSbKO+zpQVjl2p2ZUfPYcFO/6rLzff4lINFVE8fEpyeaOJ4hPTSo5kDr+0uSZ/v8OlpcksptR/8nJNdDCyXtWs34t1GUf1wsLfK7Wk2t4jedUSxH2h59ifEGYBADYlQ2igYejhPudWaNS0pOKjvsUnYj2/4Hfba1c0HLuTs1Uc3D3qVzx8Gn8dq9Q6uqbjkczK9iSb+rv/VFKVh7bK7BbnTFU2RtqF7W9+l8nEM69BmAUA2CnZenAmIafkCGrn5g30j3aN3fLaleVsFQd3jsyWDJ9nksGWbT+sT1btdUdZNoWmWdR7/FN6lawWULyVorxVJFzhrkl6FR0xd9TqgZqHMAsAKKVkCPWW166IqhqZtQak9XuKRp/dNYhY0SAbIGl4j4q3Hrzz498Traz9tOHBgeqYWL/S98cdo7ABhtSjVawWbT5U6rkN+7IdjuRXtC2g5GoPFdkKmTWIaz7CLADAr+w5kufw+KLNBxUbGVrp5bBcDXAl18PdsC9b477Y5NJ7216rWO9xfP3wcnuVHbGY0r2frKv0KK0rk9FKKtnO8t3mgw7D7LgvNiknv9BuB7nio7+OareNEu/5a7WHYq9X0VJ3/ZHjlj5aT/fievr9qwpLcwEA/MqX6/fp3k/WOXzO1SBXmQAXYEhDL26mIRcn2gWKr37dr5Efr634C/3F0ZJkmdl5WrP7iEZ+vLZSI6SBhqGlo3tUKPBY14W9t4Lv5WhTi5ItJ2Xdo3Jfv9g6xe7q1ZXsN/Fw5deIo5aLquzFdbZ2b/HvhTf0ArM0FwAATiTUD3f6XEUX57cGuOJ9seUpudFESa6MLZU3MS8uKkz9LgjTiVMFlRotLTRN7crKdelH9uVxFuIdKeselcdiStdMXl7p650pPvmuIn20ZY3YW0xp9Ozyl0krbyTV+nzxUWrrWsOGIQ29uKmia4eU6lu2vn94cKAS6ocrJ7/Qq0drCbMAAL9S3ja+1iAnlZ7pX9mWAuvo6VVtGzs9p2NifYc9nBUZyXQmJTlB4cGBlRrlzM0/reU7skp9/ops0OCodlcm/Lljq2VXldWrW5K1j/bCs+1/jbgy8c1UUeh+7q9l0sqalFZ8JLWivdmmqTI3oDAlu18X3jBa6wxhFgDgV5xNACtuwYZMffjzbrvlmSRVaJSzsuEzLipME69L0pjZG4p2N5M07JKKjWSWpWNi/UptW2zdWMK6O9cff+3OVaFWAhdGYR05k62WXVEybP+y688KhVlJ+va3A7r3k7VFS68Z0iUtovXTtiyXWxpGz96g9XuO6pOVe2TK8TJu1gl6xc9zN29euYGeWQCAX1m+I8vlrXYrMuv9TAOcVfFtgd0VKkpuW/xwn3P1x4lTmvqje/pJrcprpXDFmUwoc1ibyv9HRmZ2nrqkLq6SsOgtJg9sr34XOP8JQnWhZxYAACcqM+pX1qnuDHBS1Sxd5mzt4CFdm2rR5oOVXkWhuAAVjWCX1UrhCmvN05alOw3dFbmPrvwjo+TouD86kpvv6RJcRpgFAPiVim61WxHuDnBVyVFIjosKU93w4DN6XXeNSDsSFxWmsX1ba0jXptqVletwBzlJpcJnZXt1pdI71/133X69v2xXpT9DgCHddFGCPv45wyMjvq6uP+yNCLMAAL9jDSyVHZWsygBX3SrTbXgmYbEyrEHc2Q5yJbdNdufOdRl/5krLXLve0ffngrOiNHrOBrduyWttf7H2VzeoHazn5m9x2HNdfP1hw6jarYGrG2EWAOCX4qLC1CymtkvXuLuloCZwtoqCIzUhxDsbYa6qelz5/pQ1aS8lOUGtGtUpd9kw60jqG9/vcPx8sXsgqVSwd7ZddMlWk9lr9urFhVsr8KlqPsIsAMBvOeqfddaH6U0tBa5wtIrCJef8PTO/ukdhaxpnq0yUNQrqTNv4enqu2GuVNSktvn64w/NK3gNXgn3x56LCgs7gu1KzsJoBAMCvlZzpb+3DtB6rCaOR1aHkKgpVsaqCN3P0/ajs96ii11XlPXC2y9obN7WvEf9gYzUDAAAqyNlMf0fHfFnJEb2q/NG9N3Jne0NFr6vKe+Bsl7X4+t53zwmzAAC/V919mICn7TmS5/D43iN5ZW6xWxMFeLoAAAAAVC9nXabe2HxKmAUAAPAz1lUaSmqfULe6SzljhFkAAAA/Y12loWQQ/Oa3gx6p50zQMwsAAOCHik9+XLXrT036dqtS529Rl+bROrdRHU+XV2GMzAIAAPgp685q917WQpeeG6P8Aovu/3StThUUerq0CiPMAgAA+DnDMPT89ReoQUSwthw4rhcW/O7pkiqMMAsAAADF1gnV89dfIEl6d2m6ftp22MMVVQxhFgAAAJKky89rqJs7JUiSHpq5Xkdy8j1cUfkIswAAALAZ16+1msVE6OCxUxozZ4PTNWlrCsIsAAAAbMKCA/Xaje0VFGhowaYDmvnLXk+XVCbCLAAAAOy0aRKlUVecK0l64stN2pWV4+GKnCPMAgAAoJQ7uzVTp6b1lZtfqHtmrNFP2w4rMzvP02WVQpgFAABAKYEBhialtFNIrQBtzjymW99bqa4TFyttVYanS7NDmAUAAIBDAYaUX2CxPbaY0tg5G2vUCC1hFgAAAA6lZ+Wo5FoGhaapXVm5HqnHEcIsAAAAHGoaHaEAw/5YoGEoMTrcMwU5QJgFAACAQ3FRYUodkKRAoyjRBhqGJgxoo7ioMA9X9rdani4AAAAANVdKcoK6tYzRrqxcJUaH16ggKxFmAQAAUI64qLAaF2KtaDMAAACA1yLMAgAAwGsRZgEAAOC1CLMAAADwWjUizE6ePFmJiYkKDQ1Vp06dtHLlSqfnTp8+XYZh2H2FhoZWY7UAAACoKTweZtPS0jRq1CiNHz9ea9asUdu2bdWrVy8dOnTI6TWRkZHKzMy0fe3evbsaKwYAAEBN4fEwO2nSJA0bNkxDhgxR69atNWXKFIWHh+v99993eo1hGGrUqJHtq2HDhk7PPXXqlI4dO2b3BQAAAN/g0TCbn5+v1atXq2fPnrZjAQEB6tmzp1asWOH0uhMnTujss89WfHy8rrnmGm3atMnpuampqYqKirJ9xcfHu/UzAAAAwHM8GmazsrJUWFhYamS1YcOGOnDggMNrzj33XL3//vv673//q//85z+yWCzq0qWL9u7d6/D8MWPGKDs72/a1Z88et38OAAAAeIbX7QDWuXNnde7c2fa4S5cuOu+88/T222/r6aefLnV+SEiIQkJCqrNEAAAAVBOPjsxGR0crMDBQBw8etDt+8OBBNWrUqEKvERQUpPbt22v79u1VUSIAAABqMI+G2eDgYHXo0EGLFi2yHbNYLFq0aJHd6GtZCgsLtWHDBsXFxVVVmQAAAKihPN5mMGrUKA0aNEgdO3bURRddpFdeeUU5OTkaMmSIJOm2225TkyZNlJqaKkl66qmn9H//939q0aKFjh49qhdeeEG7d+/W0KFDPfkxAAAA4AEeD7MpKSk6fPiwHn/8cR04cEDt2rXTggULbJPCMjIyFBDw9wDykSNHNGzYMB04cED16tVThw4dtHz5crVu3dpTHwEAAAAeYpimaXq6iOp07NgxRUVFKTs7W5GRkZ4uBwAAACW4ktc8PjJb3azZnc0TAAAAaiZrTqvImKvfhdnjx49LEpsnAAAA1HDHjx9XVFRUmef4XZuBxWLR/v37VadOHRmGUeXvd+zYMcXHx2vPnj20NXgp7qH34x56P+6hd+P+eb/qvoemaer48eNq3Lix3dwpR/xuZDYgIEBnnXVWtb9vZGQkv4G9HPfQ+3EPvR/30Ltx/7xfdd7D8kZkrTy6ziwAAABwJgizAAAA8FqE2SoWEhKi8ePHKyQkxNOloJK4h96Pe+j9uIfejfvn/WryPfS7CWAAAADwHYzMAgAAwGsRZgEAAOC1CLMAAADwWoRZAAAAeC3CrBtMnjxZiYmJCg0NVadOnbRy5coyz585c6ZatWql0NBQJSUlad68edVUKZxx5R5OnTpVl1xyierVq6d69eqpZ8+e5d5zVD1Xfx9affrppzIMQ/3796/aAlEuV+/h0aNHNWLECMXFxSkkJEQtW7bkz1MPcvX+vfLKKzr33HMVFham+Ph4PfDAAzp58mQ1VYuSfvzxR1199dVq3LixDMPQF198Ue41S5Ys0YUXXqiQkBC1aNFC06dPr/I6HTJxRj799FMzODjYfP/9981NmzaZw4YNM+vWrWsePHjQ4fnLli0zAwMDzeeff9787bffzHHjxplBQUHmhg0bqrlyWLl6DwcOHGhOnjzZXLt2rbl582Zz8ODBZlRUlLl3795qrhxWrt5Dq/T0dLNJkybmJZdcYl5zzTXVUywccvUenjp1yuzYsaPZt29fc+nSpWZ6erq5ZMkSc926ddVcOUzT9fs3Y8YMMyQkxJwxY4aZnp5uLly40IyLizMfeOCBaq4cVvPmzTMfffRRc86cOaYk8/PPPy/z/J07d5rh4eHmqFGjzN9++818/fXXzcDAQHPBggXVU3AxhNkzdNFFF5kjRoywPS4sLDQbN25spqamOjz/hhtuMPv162d3rFOnTuZdd91VpXXCOVfvYUkFBQVmnTp1zA8++KCqSkQ5KnMPCwoKzC5dupjvvvuuOWjQIMKsh7l6D9966y2zWbNmZn5+fnWViDK4ev9GjBhhXnbZZXbHRo0aZXbt2rVK60TFVCTMPvzww+b5559vdywlJcXs1atXFVbmGG0GZyA/P1+rV69Wz549bccCAgLUs2dPrVixwuE1K1assDtfknr16uX0fFStytzDknJzc3X69GnVr1+/qspEGSp7D5966inFxsbqjjvuqI4yUYbK3MO5c+eqc+fOGjFihBo2bKg2bdpowoQJKiwsrK6y8ZfK3L8uXbpo9erVtlaEnTt3at68eerbt2+11IwzV5PyTK1qf0cfkpWVpcLCQjVs2NDueMOGDbVlyxaH1xw4cMDh+QcOHKiyOuFcZe5hSY888ogaN25c6jc1qkdl7uHSpUv13nvvad26ddVQIcpTmXu4c+dOLV68WDfffLPmzZun7du3a/jw4Tp9+rTGjx9fHWXjL5W5fwMHDlRWVpYuvvhimaapgoIC3X333Ro7dmx1lAw3cJZnjh07pry8PIWFhVVbLYzMAmdg4sSJ+vTTT/X5558rNDTU0+WgAo4fP65bb71VU6dOVXR0tKfLQSVZLBbFxsbqnXfeUYcOHZSSkqJHH31UU6ZM8XRpqIAlS5ZowoQJevPNN7VmzRrNmTNHX3/9tZ5++mlPlwYvxMjsGYiOjlZgYKAOHjxod/zgwYNq1KiRw2saNWrk0vmoWpW5h1YvvviiJk6cqO+++04XXHBBVZaJMrh6D3fs2KFdu3bp6quvth2zWCySpFq1aun3339X8+bNq7Zo2KnM78O4uDgFBQUpMDDQduy8887TgQMHlJ+fr+Dg4CqtGX+rzP177LHHdOutt2ro0KGSpKSkJOXk5OjOO+/Uo48+qoAAxtpqOmd5JjIyslpHZSVGZs9IcHCwOnTooEWLFtmOWSwWLVq0SJ07d3Z4TefOne3Ol6Rvv/3W6fmoWpW5h5L0/PPP6+mnn9aCBQvUsWPH6igVTrh6D1u1aqUNGzZo3bp1tq9//OMf6tGjh9atW6f4+PjqLB+q3O/Drl27avv27bZ/iEjS1q1bFRcXR5CtZpW5f7m5uaUCq/UfJqZpVl2xcJsalWeqfcqZj/n000/NkJAQc/r06eZvv/1m3nnnnWbdunXNAwcOmKZpmrfeeqs5evRo2/nLli0za9WqZb744ovm5s2bzfHjx7M0l4e5eg8nTpxoBgcHm7NmzTIzMzNtX8ePH/fUR/B7rt7DkljNwPNcvYcZGRlmnTp1zJEjR5q///67+dVXX5mxsbHmM88846mP4NdcvX/jx48369SpY37yySfmzp07zW+++cZs3ry5ecMNN3jqI/i948ePm2vXrjXXrl1rSjInTZpkrl271ty9e7dpmqY5evRo89Zbb7Wdb12a69///re5efNmc/LkySzN5c1ef/11MyEhwQwODjYvuugi83//+5/tue7du5uDBg2yO/+zzz4zW7ZsaQYHB5vnn3+++fXXX1dzxSjJlXt49tlnm5JKfY0fP776C4eNq78PiyPM1gyu3sPly5ebnTp1MkNCQsxmzZqZzz77rFlQUFDNVcPKlft3+vRp84knnjCbN29uhoaGmvHx8ebw4cPNI0eOVH/hME3TNL///nuHf7dZ79ugQYPM7t27l7qmXbt2ZnBwsNmsWTNz2rRp1V63aZqmYZqM5wMAAMA70TMLAAAAr0WYBQAAgNcizAIAAMBrEWYBAADgtQizAAAA8FqEWQAAAHgtwiwAAAC8FmEWAAAAXoswCwB+zDAMffHFF5KkXbt2yTAMrVu3zqM1AYArCLMA4CGDBw+WYRgyDENBQUFq2rSpHn74YZ08edLTpQGA16jl6QIAwJ/17t1b06ZN0+nTp7V69WoNGjRIhmHoueee83RpAOAVGJkFAA8KCQlRo0aNFB8fr/79+6tnz5769ttvJUkWi0Wpqalq2rSpwsLC1LZtW82aNcvu+k2bNumqq65SZGSk6tSpo0suuUQ7duyQJK1atUpXXHGFoqOjFRUVpe7du2vNmjXV/hkBoCoRZgGghti4caOWL1+u4OBgSVJqaqo+/PBDTZkyRZs2bdIDDzygW265RT/88IMkad++ferWrZtCQkK0ePFirV69WrfffrsKCgokScePH9egQYO0dOlS/e9//9M555yjvn376vjx4x77jADgbrQZAIAHffXVV6pdu7YKCgp06tQpBQQE6I033tCpU6c0YcIEfffdd+rcubMkqVmzZlq6dKnefvttde/eXZMnT1ZUVJQ+/fRTBQUFSZJatmxpe+3LLrvM7r3eeecd1a1bVz/88IOuuuqq6vuQAFCFCLMA4EE9evTQW2+9pZycHL388suqVauWrrvuOm3atEm5ubm64oor7M7Pz89X+/btJUnr1q3TJZdcYguyJR08eFDjxo3TkiVLdOjQIRUWFio3N1cZGRlV/rkAoLoQZgHAgyIiItSiRQtJ0vvvv6+2bdvqvffeU5s2bSRJX3/9tZo0aWJ3TUhIiCQpLCyszNceNGiQ/vjjD7366qs6++yzFRISos6dOys/P78KPgkAeAZhFgBqiICAAI0dO1ajRo3S1q1bFRISooyMDHXv3t3h+RdccIE++OADnT592uHo7LJly/Tmm2+qb9++kqQ9e/YoKyurSj8DAFQ3JoABQA3yz3/+U4GBgXr77bf10EMP6YEHHtAHH3ygHTt2aM2aNXr99df1wQcfSJJGjhypY8eO6cYbb9Qvv/yibdu26aOPPtLvv/8uSTrnnHP00UcfafPmzfr555918803lzuaCwDehpFZAKhBatWqpZEjR+r5559Xenq6YmJilJqaqp07d6pu3bq68MILNXbsWElSgwYNtHjxYv373/9W9+7dFRgYqHbt2qlr166SpPfee0933nmnLrzwQsXHx2vChAl66KGHPPnxAMDtDNM0TU8XAQAAAFQGbQYAAADwWoRZAAAAeC3CLAAAALwWYRYAAABeizALAAAAr0WYBQAAgNcizAIAAMBrEWYBAADgtQizAAAA8FqEWQAAAHgtwiwAAAC81v8Dg729m7pOOH4AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 800x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the precision-recall curve\n",
    "prior = \"StringDB\"\n",
    "plt.figure(figsize=(8, 6))\n",
    "plt.plot(\n",
    "    pr_dict[prior][\"recall\"],\n",
    "    pr_dict[prior][\"precision\"],\n",
    "    marker=\".\",\n",
    "    label=\"Precision-Recall Curve\",\n",
    ")\n",
    "# Adding labels and title\n",
    "plt.xlabel(\"Recall\")\n",
    "plt.ylabel(\"Precision\")\n",
    "plt.title(\"Precision-Recall Curve\")\n",
    "plt.legend()\n",
    "# Save the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fdff66a-8f20-4268-8520-d33a0fa4969a",
   "metadata": {},
   "source": [
    "# Interpretability of our learned latent space"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ec54f80-8a0f-41f0-8c55-edd2b1fbfb3c",
   "metadata": {},
   "source": [
    "## Factor Enrichment\n",
    "- GSEA on log transformed latent-by-gene factors\n",
    "- Then associate each perturbation to its top latent factors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35eb8169-15b1-4b5b-b926-d5b7238f372a",
   "metadata": {},
   "outputs": [],
   "source": [
    "factor_to_go = factor_enrichment_gsea(adata, adata.uns[\"SPECTRA_factors\"], fdr=1e-2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "814189ce-53f3-4716-834a-a4f0ea6efcf9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# filter and add description to processes\n",
    "go_df = read_aws_csv(\"s3://pert-spectra/references/GO_to_Description.txt\")\n",
    "go_df.set_index(\"Term\", inplace=True)\n",
    "go_dict = go_df.to_dict()[\"Description\"]\n",
    "\n",
    "filtered_factor_to_go = {}\n",
    "for i in factor_to_go:\n",
    "    proc = factor_to_go[i]\n",
    "    proc[\"descr\"] = [go_dict[x] for x in proc[\"GO_ID\"]]\n",
    "    filtered_factor_to_go[i] = proc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "999ea1a6-b9f9-4cdd-a8cc-15e4db23f7f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save latent enrichment results\n",
    "with open(\n",
    "    \"figures/factor_enrichments/pertspectra_inhouse_factor_enrichment.pickle\", \"wb\"\n",
    ") as handle:\n",
    "    pickle.dump(filtered_factor_to_go, handle, protocol=pickle.HIGHEST_PROTOCOL)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78d01afb-6501-4424-99f5-05839747f118",
   "metadata": {},
   "source": [
    "## Overlap with prior knowledge/ground truth (stringdb)\n",
    "- Group A: a set of GO terms associated with a perturbation (either drivers from msigdb, or from literature) and its neighbors in stringdb\n",
    "- Group B: a set of GO terms from the interpretability analysis\n",
    "- Hypergeometric test on the overlap of the two groups"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "d0274191-ea1a-4abf-8704-53f99a08fa3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load precomputed factor_to_go dict if available\n",
    "filtered_factor_to_go = read_aws_pickle(\n",
    "    \"s3://pert-spectra/figures/factor_enrichments/pertspectra_inhouse_factor_enrichment.pickle\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "782c4388-f53b-4963-8ee1-578de81d67de",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get neighbors for each perturbation\n",
    "pert_neighbors = retrieve_stringdb_neighbors(adata.uns[\"Spectra_pert_labels\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "abfc3fd4-806e-459d-af34-0b9f365a8e91",
   "metadata": {},
   "outputs": [],
   "source": [
    "# construct group A - known processes for each perturbation\n",
    "gene_sets = read_aws_pickle(\"s3://pert-spectra/references/GO_to_Gene.pickle\")\n",
    "# BP only go terms\n",
    "go_reference = read_aws_csv(\"s3://pert-spectra/references/GO_terms.txt.gz\", zipped=True)\n",
    "go_bp = go_reference[go_reference[\"go_category\"] == \"biological_process\"]\n",
    "go_bp_ids = set(go_bp[\"go_id\"].values)\n",
    "filtered_go_terms = {key: gene_sets[key] for key in go_bp_ids if key in gene_sets}\n",
    "\n",
    "# GO terms per perturbation AND its neighbors in stringdb\n",
    "pert_to_go = {key: set() for key in adata.uns[\"Spectra_pert_labels\"]}\n",
    "for goterm in filtered_go_terms:\n",
    "    for pert in adata.uns[\"Spectra_pert_labels\"]:\n",
    "        if pert in filtered_go_terms[goterm] and set(\n",
    "            filtered_go_terms[goterm]\n",
    "        ).intersection(pert_neighbors[pert]):\n",
    "            pert_to_go[pert].add(goterm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "3245b161-3e6c-4966-99d9-b7bb1d2ba6e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# construct set B - model identified processes for each perturbation\n",
    "n = 2  # number of top factors to get processes from\n",
    "model_pert_to_go = {}\n",
    "for pert in adata.uns[\"Spectra_pert_labels\"]:\n",
    "    if pert in [\"ctrl\", \"intergenic\", \"basal\"]:\n",
    "        continue\n",
    "    # get top factors\n",
    "    tnfa_pos_pert_embeddings_df = pd.DataFrame(\n",
    "        adata.uns[\"SPECTRA_pert_scores\"], index=adata.uns[\"Spectra_pert_labels\"]\n",
    "    )\n",
    "    pert_loading = tnfa_pos_pert_embeddings_df.loc[pert]\n",
    "    intergenic_loading = tnfa_pos_pert_embeddings_df.loc[\"intergenic\"]\n",
    "    delta_loading = np.abs(np.log(pert_loading) - np.log(intergenic_loading))\n",
    "    top_n_factors = np.argpartition(np.array(delta_loading), -n)[-n:]\n",
    "    # get processes\n",
    "    model_processes = set()\n",
    "    for f in top_n_factors:\n",
    "        proc = filtered_factor_to_go[f]\n",
    "        model_processes = model_processes.union(set(proc[\"Term\"]))\n",
    "    model_pert_to_go[pert] = model_processes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "e6a99270-618f-467b-8c0c-67610171ddbd",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overlap for SKP2: 2 out of 6 in researchDB\n",
      "P-value for SKP2: 0.005521122949289679\n",
      "Overlap for NFKBIA: 4 out of 28 in researchDB\n",
      "P-value for NFKBIA: 0.0018530509193777707\n",
      "Overlap for TRAF2: 12 out of 29 in researchDB\n",
      "P-value for TRAF2: 9.372890092145309e-14\n",
      "Overlap for BIRC2: 12 out of 30 in researchDB\n",
      "P-value for BIRC2: 1.5352819800190863e-13\n",
      "Overlap for TAB1: 6 out of 9 in researchDB\n",
      "P-value for TAB1: 9.880142853042106e-07\n",
      "Overlap for UBE2N: 6 out of 16 in researchDB\n",
      "P-value for UBE2N: 3.5981263241025083e-07\n",
      "Overlap for CHUK: 12 out of 21 in researchDB\n",
      "P-value for CHUK: 8.970401417521286e-16\n",
      "Overlap for RBCK1: 8 out of 11 in researchDB\n",
      "P-value for RBCK1: 4.469976857613242e-09\n",
      "Overlap for IKBKB: 15 out of 24 in researchDB\n",
      "P-value for IKBKB: 1.262837991043831e-14\n",
      "Overlap for TNFRSF1A: 12 out of 21 in researchDB\n",
      "P-value for TNFRSF1A: 3.557360199375414e-11\n",
      "Overlap for TRAF6: 14 out of 46 in researchDB\n",
      "P-value for TRAF6: 8.701845914835026e-14\n",
      "Overlap for BIRC3: 8 out of 20 in researchDB\n",
      "P-value for BIRC3: 2.0442561238209195e-09\n",
      "Overlap for TRADD: 11 out of 13 in researchDB\n",
      "P-value for TRADD: 2.20628781790212e-13\n",
      "Overlap for TRAF5: 9 out of 14 in researchDB\n",
      "P-value for TRAF5: 8.644484310231755e-13\n",
      "Overlap for IKBKG: 12 out of 17 in researchDB\n",
      "P-value for IKBKG: 7.268024700298374e-13\n",
      "Overlap for TAB2: 5 out of 5 in researchDB\n",
      "P-value for TAB2: 2.730409006194343e-07\n",
      "Overlap for RIPK1: 28 out of 38 in researchDB\n",
      "P-value for RIPK1: 3.1419351093251465e-29\n",
      "Overlap for SKP1: 3 out of 7 in researchDB\n",
      "P-value for SKP1: 0.00023036479244599009\n",
      "Overlap for MAP3K7: 10 out of 23 in researchDB\n",
      "P-value for MAP3K7: 4.7379439102193765e-08\n",
      "Overlap for FBXW11: 11 out of 13 in researchDB\n",
      "P-value for FBXW11: 2.0898341048664293e-13\n",
      "Overlap for CUL1: 2 out of 7 in researchDB\n",
      "P-value for CUL1: 0.007218985999647311\n"
     ]
    }
   ],
   "source": [
    "pvals = perturbation_signal_recovery(\n",
    "    pert_to_go,\n",
    "    model_pert_to_go,\n",
    "    list(filtered_go_terms.keys()),\n",
    "    list(adata.uns[\"Spectra_pert_labels\"]),\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "b057ceb0-49b0-4c4d-870d-d9a3fe7b1428",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save as csv for visualization\n",
    "pd.DataFrame.from_dict(data=pvals, orient=\"index\").to_csv(\n",
    "    \"figures/process_recovery_hypergeo_pvals/MODEL_inhouse_hypergeo_neighbors_recovery_pvalues.csv\",\n",
    "    header=False,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99859d27-3a23-42f0-9d01-1087d698d118",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pertspectra",
   "language": "python",
   "name": "pertspectra"
  },
  "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.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}