inhouse_scETM_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": 2,
   "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",
    "    get_gprofiler,\n",
    "    read_aws_csv,\n",
    "    read_aws_h5ad,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d15b3b79-f031-4029-b4cc-3085534cf1ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in scETM results\n",
    "adata = read_aws_h5ad(\"s3://pert-spectra/scETM_checkpoints/scetm_inhouse/fold_0.h5ad\")\n",
    "# read in pertspectra results to retrieve gene labels\n",
    "ref_adata = read_aws_h5ad(\n",
    "    \"s3://pert-spectra/PertSpectra_checkpoints/pertspectra_inhouse/fold_0.h5ad\"\n",
    ")\n",
    "adata.var_names = ref_adata.var_names"
   ]
  },
  {
   "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": 5,
   "id": "e2885bd3-5433-47d6-88d5-47be83df5b5a",
   "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/scETM_checkpoints/\"\n",
    "    experiment_name = \"scetm_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": 6,
   "id": "e0c30413-c2c1-4058-b3bf-4dc74437c0dd",
   "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",
    "for n in range(n_folds):\n",
    "    adata_n = model_adatas[n]\n",
    "    loss_weights = np.ones(adata_n.shape[0])\n",
    "\n",
    "    for pert in adata_n.obs[\"condition\"].unique():\n",
    "        hold_idx = [\n",
    "            i\n",
    "            for i, x in enumerate(\n",
    "                adata_n[adata_n.obs[\"Treatment\"] == \"TNFA+\"].obs[\"condition\"]\n",
    "            )\n",
    "            if x == pert\n",
    "        ]\n",
    "        recon = adata_n[adata_n.obs[\"Treatment\"] == \"TNFA+\"][hold_idx].uns[\"recon\"]\n",
    "        # correlation\n",
    "        mean_reconstruction = recon.mean(axis=0)\n",
    "        mean_observed = np.squeeze(\n",
    "            np.array(\n",
    "                adata_n[adata_n.obs[\"Treatment\"] == \"TNFA+\"][hold_idx].X.mean(axis=0)\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": 7,
   "id": "e61e2799-ed83-421e-86d4-565da4876c10",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAncAAAHFCAYAAACQIN+BAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABIPElEQVR4nO3de3zP9f//8fubne2QGTuwE7FNWjmUkPMpItKnk8ihfOTwkXSUtKkPovJRhFQOlUNJB/UTlpAawieKRidMsZhkw4xtz98ffff+eNvGpr333l67XS+X9+Xi/Xy/3q/n4/l6vt6v993r/X69ZzPGGAEAAMASqri6AAAAAJQewh0AAICFEO4AAAAshHAHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALCQywp33377rQYNGqTo6Gh5eXnJ19dXTZo00dSpU/XHH3+Udo1/y/r162Wz2bR+/foSP/f7779XYmKi9u/fX+CxgQMHKioq6m/XV1L79++XzWbTCy+8UGrrPHTokBITE7Vjx45SWycKWrBggWw2m8P+1K5dO7Vr167E65o0aZI+/PDDUqstX/7+tWDBglJf9+Uq7ddau3bt1KhRo1JbX0VRHue2vEpMTJTNZnN1GZWCs99/8jPAe++9V2rrLI/Z4EIlDnevvfaamjZtqq1bt+rRRx/VqlWr9MEHH+j222/XnDlzdN999zmjTpf4/vvvNWHChEIncPz48frggw/KvignOHTokCZMmEC4c4FZs2Zp1qxZJX6es8IdAOn+++/Xpk2bXF1GpVAR338qQjZwK8nCmzZt0rBhw9S5c2d9+OGH8vT0tD/WuXNnPfzww1q1alWpFHb69Gn5+PgUaM/NzVVOTo5D365Qr149l/aPsuPMfa5hw4alvk5UDOXlWFbeZWVlydvbu0z7rFOnjurUqVOmfZal8rDv5ddQUddflPKSDUp05m7SpEmy2WyaO3duoTuFh4eHbrnlFvv9vLw8TZ06VbGxsfL09FStWrV077336tdff3V4Xv7HJF988YVatmwpHx8fDR482P4xwtSpU/Xvf/9b0dHR8vT01Lp16yRJ27Zt0y233KLAwEB5eXmpcePGevfddy85jm3btumuu+5SVFSUvL29FRUVpbvvvlsHDhywL7NgwQLdfvvtkqT27dvLZrM5fKRR2KnXM2fOaOzYsYqOjpaHh4dq166tESNG6M8//3RYLioqSj169NCqVavUpEkTeXt7KzY2VvPmzbtk7edv24kTJyoiIkJeXl5q1qyZ1q5dW2C5H3/8UX379lWtWrXk6empuLg4vfLKK/bH169fr+uuu06SNGjQIPs4ExMT9f/+3/+TzWbT1q1b7csvX75cNptNN998s0M/8fHxuu222+z3jTGaNWuWrr32Wnl7e6t69er6xz/+oV9++aVAjZ999pk6duwof39/+fj4qFWrVgXGkv8xye7du3X33XcrICBAwcHBGjx4sE6cOHHJ7ZW/j23cuFE33HCDvL29Vbt2bY0fP165ubn25Uprn9u8ebNatWolLy8vhYWFaezYsTp37lyhdV34sWx2draeeeYZxcXFycvLSzVq1FD79u2VnJwsSbLZbDp16pQWLlxon6/z15GWlqahQ4eqTp068vDwUHR0tCZMmFDgQHfo0CHdcccd8vPzU0BAgO68806lpaVdcltK0tGjRzV8+HA1bNhQvr6+qlWrljp06KCNGzc6LHf+1wimTZum6Oho+fr6qkWLFtq8eXOB9S5YsEAxMTH2ffXNN98sVj1S8Y83+S61L0jS7Nmzdc0118jX11d+fn6KjY3Vk08+6bBMcbZ3UfvVu+++Kw8PD40fP75AfXv27JHNZtPLL79cor6kvze3xenHGKPu3burRo0aSk1NtT/v9OnTuuqqqxQXF6dTp05J+t9r95tvvlGfPn3k7++vgIAA9evXT0ePHnXoN//Y+P7776tx48by8vLShAkTSjT2S83Z6dOn9cgjj9i/VhQYGKhmzZppyZIl9mUK+1i2pO9nW7duVevWreXj46O6devqueeeU15e3iW3vc1m08iRI/Xqq6+qQYMG8vT0VMOGDbV06dISz5N06WOaq2so6v0nf1sW9rWVC99/izPGM2fOaMyYMQoJCZG3t7fatm2rb775xmG9VsgGdqaYcnJyjI+Pj2nevHlxn2L++c9/Gklm5MiRZtWqVWbOnDmmZs2aJjw83Bw9etS+XNu2bU1gYKAJDw83M2bMMOvWrTMbNmww+/btM5JM7dq1Tfv27c17771n1qxZY/bt22c+//xz4+HhYVq3bm3eeecds2rVKjNw4EAjycyfP9++7nXr1hlJZt26dfa2ZcuWmaefftp88MEHZsOGDWbp0qWmbdu2pmbNmva6jhw5YiZNmmQkmVdeecVs2rTJbNq0yRw5csQYY8yAAQNMZGSkfZ15eXmma9euxs3NzYwfP96sWbPGvPDCC6ZatWqmcePG5syZM/ZlIyMjTZ06dUzDhg3Nm2++aVavXm1uv/12I8ls2LDhots0f5uEh4ebG2+80SxfvtwsW7bMXHfddcbd3d0kJyfbl929e7cJCAgwV199tXnzzTfNmjVrzMMPP2yqVKliEhMTjTHGnDhxwsyfP99IMk899ZR9nAcPHjSZmZnG3d3dTJo0yb7OBx54wHh7e5tq1aqZs2fPGmOM+f33343NZjOzZs2yLzdkyBDj7u5uHn74YbNq1SqzePFiExsba4KDg01aWpp9ubfeesvYbDbTu3dv8/7775uPP/7Y9OjRw1StWtV89tln9uUSEhKMJBMTE2Oefvppk5SUZKZNm2Y8PT3NoEGDLrrNjPlrH6tRo4YJCwszL7/8slm9erUZNWqUkWRGjBhRYPv+nX1u9+7dxsfHxzRs2NAsWbLEfPTRR6Zr164mIiLCSDL79u1zqKtt27b2++fOnTPt27c3bm5u5pFHHjErV640K1asME8++aRZsmSJMcaYTZs2GW9vb9O9e3f7fO3evdsYY8zhw4dNeHi4iYyMNK+++qr57LPPzLPPPms8PT3NwIED7f2cPn3axMXFmYCAADNjxgz79siv8fzxFGbPnj1m2LBhZunSpWb9+vXmk08+Mffdd5+pUqWKw2stf3tGRUWZm266yXz44Yfmww8/NFdffbWpXr26+fPPP+3L5u+HvXr1Mh9//LF5++23zZVXXmkfz6WU5HhTnH1hyZIlRpL517/+ZdasWWM+++wzM2fOHDNq1Cj7MsXd3hfbr2699VYTHh5ucnNzHcbz2GOPGQ8PD5Oenl6mc1vcftLT002dOnVM8+bN7ceCAQMGGG9vb/Ptt9/al8t/7UZGRppHH33UrF692kybNs1+bMx/rjF/HRtDQ0NN3bp1zbx588y6devM119/XeyaijNnQ4cONT4+PmbatGlm3bp15pNPPjHPPfecmTFjRoGaz1fS/at+/fpmzpw5JikpyQwfPtxIMgsXLrzotjfG2I/v+cePFStWmJtuuslIMsuWLSvxPF1s33N1DTt37izy/Sd/W55/fMx34fvvxcaYnwHCw8MLHFv8/f3Nzz//bF9PRc4GBeawuAumpaUZSeauu+4q1vIpKSlGkhk+fLhD+5YtW4wk8+STT9rb2rZtaySZtWvXOiybP2H16tVzOAAYY0xsbKxp3LixOXfunEN7jx49TGhoqP1AWVi4u1BOTo45efKkqVatmnnppZfs7cuWLSvyuRdO4KpVq4wkM3XqVIfl3nnnHSPJzJ07194WGRlpvLy8zIEDB+xtWVlZJjAw0AwdOrTIOo353zYJCwszWVlZ9vaMjAwTGBhoOnXqZG/r2rWrqVOnjjlx4oTDOkaOHGm8vLzMH3/8YYwxZuvWrUUe9G+88UbToUMH+/0rr7zSPProo6ZKlSr2nW3RokVGkvnhhx+MMX+FD0nmxRdfdFjXwYMHjbe3t3nssceMMcacOnXKBAYGmp49ezosl5uba6655hpz/fXX29vyD7YXbt/hw4cbLy8vk5eXd9Htlr+PffTRRw7tQ4YMMVWqVLHPRWnsc3feeafx9vZ2CLE5OTkmNjb2kuHuzTffNJLMa6+9dtHxVKtWzQwYMKBA+9ChQ42vr6/DvmWMMS+88IKRZA+Bs2fPLnJ7FCcAXCgnJ8ecO3fOdOzY0dx666329vztefXVV5ucnBx7+9dff20k2QNrbm6uCQsLM02aNHGYy/379xt3d/dLhrvLOd5cal8YOXKkueKKKy7ab3G398X2qxUrVhhJZs2aNfa2nJwcExYWZm677bYS9/V357a4/RhjzJdffmnc3NzM6NGjzbx584wk8/rrrzs8L/+1+9BDDzm05x833n77bXtbZGSkqVq1qtm7d+9l1VScOWvUqJHp3bv3RZe5MNxdzv61ZcsWh2UbNmxounbtetF+jfkrWBV1/LjyyivtbaWx75WHGi72/lPScFfY+vMzQFHHlvvvv7/I7VCRssGFnPZTKPmnQwcOHOjQfv311ysuLq7Ax27Vq1dXhw4dCl3XLbfcInd3d/v9n376SXv27NE999wjScrJybHfunfvrsOHD2vv3r1F1nby5Ek9/vjjuvLKK+Xm5iY3Nzf5+vrq1KlTSklJuZzh6vPPP5dUcLy33367qlWrVmC81157rSIiIuz3vby81KBBA4fTvxfTp08feXl52e/7+fmpZ8+e+uKLL5Sbm6szZ85o7dq1uvXWW+Xj41NgG505c6bQj8Uu1LFjR3311VfKysrSgQMH9NNPP+muu+7Stddeq6SkJEl/fawaERGh+vXrS5I++eQT2Ww29evXz6HfkJAQXXPNNfYrl5OTk/XHH39owIABDsvl5eXppptu0tatW+0f7eQ7/2N/6a+Pg8+cOaMjR45ccix+fn4Fnt+3b1/l5eXpiy++KNDP5e5z69atU8eOHRUcHGx/ftWqVXXnnXdessZPP/1UXl5eGjx48CWXLcwnn3yi9u3bKywszKHGbt26SZI2bNhgr7Go7VFcc+bMUZMmTeTl5SU3Nze5u7tr7dq1hb6Gbr75ZlWtWtV+Pz4+XpLs+/vevXt16NAh9e3b1+HjsMjISLVs2fKStZT0eFOcfeH666/Xn3/+qbvvvlsfffSR0tPTC/Rb3O2d78L9SpK6deumkJAQzZ8/3962evVqHTp0yGE/KKu5LcmYWrVqpYkTJ2r69OkaNmyY+vXrV+RFdfmvnXx33HGH3NzcCnw8GB8frwYNGlxWTcWZs+uvv16ffvqpnnjiCa1fv15ZWVmX3CYl3b9CQkJ0/fXXFxhXcY/vRR0/fvrpJ/vHwKWx75X3GkrqYusv6thy/v5nhWyQr9jhLigoSD4+Ptq3b1+xlj927JgkKTQ0tMBjYWFh9sfzFbZcUY/9/vvvkqRHHnlE7u7uDrfhw4dLUqEv6nx9+/bVzJkzdf/992v16tX6+uuvtXXrVtWsWbNYL/TCHDt2TG5ubqpZs6ZDu81mU0hISIHx1qhRo8A6PD09i91/SEhIoW1nz57VyZMndezYMeXk5GjGjBkFtlH37t0lXXwb5evUqZOys7P15ZdfKikpSUFBQWrcuLE6deqkzz77TJK0du1aderUyf6c33//XcYYBQcHF+h78+bN9n7z5/Ef//hHgeWmTJkiY0yBn9a5cLvlf/ezONvt/APV+dtM0iX3x5Lsc8eOHStyfi7l6NGjCgsLU5Uql/f/rt9//10ff/xxgRqvuuqqAjVebHtcyrRp0zRs2DA1b95cy5cv1+bNm7V161bddNNNhc7FpeYtf/tf7nYr6fGmOPtC//79NW/ePB04cEC33XabatWqpebNm9v/UyMVf3vnK6w+Nzc39e/fXx988IH9OzgLFixQaGiounbtWuK+/u7clnRM99xzjzw8PJSdna1HH320yPVe2L+bm5tq1KhRrPeC4tZUnDl7+eWX9fjjj+vDDz9U+/btFRgYqN69e+vHH38ssvaS7l/OOr6fX0tp7HvlvYaSutj6ixrP+XNnhWyQr9hXy1atWlUdO3bUp59+ql9//fWSVxLlF3j48OECyx46dEhBQUEObRf7TaELH8t/7tixY9WnT59CnxMTE1No+4kTJ/TJJ58oISFBTzzxhL09Ozv7b/1GX40aNZSTk6OjR486TKIxRmlpafYvjZaWwr4cnZaWJg8PD/n6+srd3V1Vq1ZV//79NWLEiELXER0dfcl+mjdvLl9fX3322Wfav3+/OnbsKJvNpo4dO+rFF1/U1q1blZqa6hDugoKCZLPZtHHjxkIvvMlvy5/HGTNm6IYbbii0/8LepC5XfkA7X/52vPAF9Xf2uRo1ahQ5P5dSs2ZNffnll8rLy7usgBcUFKT4+HhNnDix0MfDwsLsNX799deXVaMkvf3222rXrp1mz57t0J6ZmVnCimWvp6j+i1NTSY83xd0XBg0apEGDBunUqVP64osvlJCQoB49euiHH35QZGRksbd3vqKOc4MGDdLzzz+vpUuX6s4779SKFSs0evRoh7OdZTW3JRlTbm6u7rnnHlWvXl2enp6677779NVXX8nDw6PQ/mvXrm2/n5OTo2PHjl3ytVfSmi41Z9WqVdOECRM0YcIE/f777/azeD179tSePXsKXX9J96+/62Kvg/xaSmvfK881eHl5FXrBXFEnJi62/qLGkz8Wq2SDfCX6KZSxY8dq5cqVGjJkiD766KMCL+Bz585p1apV6tmzp/0j1rffftuh+K1btyolJUXjxo277KJjYmJUv3597dy5U5MmTSrRc202m4wxBULH66+/XuBKuZKcFerYsaOmTp2qt99+Ww899JC9ffny5Tp16pQ6duxYojov5f3339fzzz9v/2g2MzNTH3/8sVq3bq2qVavKx8dH7du31zfffKP4+PhCD7b5LjZOd3d3tWnTRklJSTp48KCee+45SVLr1q3l5uamp556yh728vXo0UPPPfecfvvtN91xxx1F9tuqVStdccUV+v777zVy5MjL2g4lkZmZqRUrVjh8XLV48WJVqVJFbdq0uehzS7LPtW/fXitWrNDvv/9uD6e5ubl65513Llljt27dtGTJEi1YsOCiH80W9T+5Hj16aOXKlapXr56qV69+0RrffffdQrdHcdhstgKvoW+//VabNm1SeHh4sdZxvpiYGIWGhmrJkiUaM2aM/SB94MABJScnF3iTuFBJjzcl3ReqVaumbt266ezZs+rdu7d2796tyMjIYm/vS4mLi1Pz5s01f/585ebmKjs7W4MGDXJYpqzmtiRjSkhI0MaNG7VmzRpVq1ZNbdq00aOPPqqXXnqpwLKLFi1S06ZN7fffffdd5eTkFOtHvC9nOxc1Z+cLDg7WwIEDtXPnTk2fPr3In+By5vtZYdauXVvo8aNevXr2cFla+56ra7jY+09UVJSWLVum7Oxs+3LHjh1TcnKy/P39S9RPUceWe++9V5J1skG+EoW7Fi1aaPbs2Ro+fLiaNm2qYcOG6aqrrtK5c+f0zTffaO7cuWrUqJF69uypmJgY/fOf/9SMGTNUpUoVdevWTfv379f48eMVHh7uMMjL8eqrr6pbt27q2rWrBg4cqNq1a+uPP/5QSkqK/vvf/2rZsmWFPs/f319t2rTR888/r6CgIEVFRWnDhg164403dMUVVzgsm/8r9nPnzpWfn5+8vLwUHR1d6GnTzp07q2vXrnr88ceVkZGhVq1a6dtvv1VCQoIaN26s/v37/63xXqhq1arq3LmzxowZo7y8PE2ZMkUZGRn2nw2QpJdeekk33nijWrdurWHDhikqKkqZmZn66aef9PHHH9u/C1CvXj15e3tr0aJFiouLk6+vr8LCwuxvqB07dtTDDz8sSfYzdN7e3mrZsqXWrFmj+Ph41apVy95vq1at9M9//lODBg3Stm3b1KZNG1WrVk2HDx/Wl19+qauvvlrDhg2Tr6+vZsyYoQEDBuiPP/7QP/7xD9WqVUtHjx7Vzp07dfTo0QJnhv6OGjVqaNiwYUpNTVWDBg20cuVKvfbaaxo2bJjDdxyKUtx97qmnntKKFSvUoUMHPf300/Lx8dErr7xS4PuDhbn77rs1f/58PfDAA9q7d6/at2+vvLw8bdmyRXFxcbrrrrskSVdffbXWr1+vjz/+WKGhofLz81NMTIyeeeYZJSUlqWXLlho1apRiYmJ05swZ7d+/XytXrtScOXNUp04d3XvvvfrPf/6je++9VxMnTlT9+vW1cuVKrV69uljbskePHnr22WeVkJCgtm3bau/evXrmmWcUHR19Wb8tVaVKFT377LO6//77deutt2rIkCH6888/lZiYWKyPE0t6vCnOvjBkyBB5e3urVatWCg0NVVpamiZPnqyAgAD7G3xxt3dxDB48WEOHDtWhQ4fUsmXLAp8+lNXcFrefpKQkTZ48WePHj7e/QU2ePFmPPPKI2rVrp1tvvdVhve+//77c3NzUuXNn7d69W+PHj9c111xz0f8AlrSm4sxZ8+bN1aNHD8XHx6t69epKSUnRW2+9pRYtWhQa7KSS719/V1BQkDp06KDx48erWrVqmjVrlvbs2ePwUySlue+5soaLvf/0799fr776qvr166chQ4bo2LFjmjp1aomDnSQdOXLEfmw5ceKEEhIS5OXlpbFjx0qyTjawK9HlF/9nx44dZsCAASYiIsJ4eHjYL+l9+umn7ZcDG/PXFXBTpkwxDRo0MO7u7iYoKMj069fPfplzvrZt25qrrrqqQD/5V8A8//zzhdaxc+dOc8cdd5hatWoZd3d3ExISYjp06GDmzJljX6awq2V//fVXc9ttt5nq1asbPz8/c9NNN5ldu3aZyMjIAlcgTp8+3URHR5uqVas6XNFz4RUxxvx1Vcvjjz9uIiMjjbu7uwkNDTXDhg0zx48fd1guMjLS3HzzzQXGU9SVQYVtkylTppgJEyaYOnXqGA8PD9O4cWOzevXqQpcfPHiwqV27tnF3dzc1a9Y0LVu2NP/+978dlluyZImJjY017u7uRpJJSEiwP7Zz504jydSvX9/hORMnTjSSzJgxYwqtdd68eaZ58+amWrVqxtvb29SrV8/ce++9Ztu2bQ7Lbdiwwdx8880mMDDQuLu7m9q1a5ubb77Z4ZL7/KvXzv/JAWP+9/MZF7us35j/7WPr1683zZo1M56eniY0NNQ8+eSTDle/lsY+Z4wxX331lbnhhhuMp6enCQkJMY8++qiZO3fuJa+WNeav/ejpp5829evXNx4eHqZGjRqmQ4cODj9zs2PHDtOqVSvj4+NjJDms4+jRo2bUqFEmOjrauLu7m8DAQNO0aVMzbtw4c/LkSfty+a8DX19f4+fnZ2677TaTnJxcrCsqs7OzzSOPPGJq165tvLy8TJMmTcyHH35Y5FVshW3PC/czY4x5/fXX7eNu0KCBmTdvXqGvtcKU9HhzqX1h4cKFpn379iY4ONh4eHiYsLAwc8cddzj8zIcxxdvel9qvjPnrZ4m8vb0verV0Wcxtcfo5dOiQqVWrlunQoYPDT7jk5eWZnj17miuuuMK+n+e/drdv32569uxpr+nuu+82v//+u0O/RR0bizv24szZE088YZo1a2aqV69uPD09Td26dc1DDz1k/8mZ82s+3999Pyvufqz/+0meWbNmmXr16hl3d3cTGxtrFi1adFnbpDj7nqtruNj7z8KFC01cXJzx8vIyDRs2NO+8806JjjP5GeCtt94yo0aNMjVr1jSenp6mdevWBd6LKnI2uJDNGGOcExuB8qNdu3ZKT0/Xrl27XF0KUKkkJiZqwoQJOnr0aKl/N82KbDabRowYoZkzZ1bqGvD3OO2nUAAAAFD2CHcAAAAWwseyAAAAFsKZOwAAAAsh3AEAAFgI4Q4AAMBCSvQjxlaXl5enQ4cOyc/Pr8R/JgUAALiGMUaZmZl/629zWwnh7jyHDh26rD+dBAAAXO/gwYN/6y9zWAXh7jx+fn6S/to5LufPmwAAgLKXkZGh8PBw+/t4ZUe4O0/+R7H+/v6EOwAAKhi+UvUXPpgGAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFiIm6sLAHD5UlNTlZ6e7uoyCggKClJERISrywCASolwB1RQqampio2NU1bWaVeXUoC3t4/27Ekh4AGACxDugAoqPT1dWVmn1XxwgvxDo1xdjl3G4f3aMm+C0tPTCXcA4AKEO6CC8w+NUmBEjKvLAACUE1xQAQAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBC3FxdAACUldTUVKWnp7u6DAdBQUGKiIhwdRkALIRwB6BSSE1NVWxsnLKyTru6FAfe3j7asyeFgAeg1BDuAFQK6enpyso6reaDE+QfGuXqciRJGYf3a8u8CUpPTyfcASg1hDsAlYp/aJQCI2JcXQYAOA0XVAAAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAupkOFu8uTJstlsGj16tL3NGKPExESFhYXJ29tb7dq10+7du11XJAAAgAtUuHC3detWzZ07V/Hx8Q7tU6dO1bRp0zRz5kxt3bpVISEh6ty5szIzM11UKQAAQNmrUOHu5MmTuueee/Taa6+pevXq9nZjjKZPn65x48apT58+atSokRYuXKjTp09r8eLFLqwYAACgbFWocDdixAjdfPPN6tSpk0P7vn37lJaWpi5dutjbPD091bZtWyUnJ5d1mQAAAC7j5uoCimvp0qX673//q61btxZ4LC0tTZIUHBzs0B4cHKwDBw4Uuc7s7GxlZ2fb72dkZJRStQAAAK5RIc7cHTx4UA8++KDefvtteXl5FbmczWZzuG+MKdB2vsmTJysgIMB+Cw8PL7WaAQAAXKFChLvt27fryJEjatq0qdzc3OTm5qYNGzbo5Zdflpubm/2MXf4ZvHxHjhwpcDbvfGPHjtWJEyfst4MHDzp1HAAAAM5WIT6W7dixo7777juHtkGDBik2NlaPP/646tatq5CQECUlJalx48aSpLNnz2rDhg2aMmVKkev19PSUp6enU2sHAAAoSxUi3Pn5+alRo0YObdWqVVONGjXs7aNHj9akSZNUv3591a9fX5MmTZKPj4/69u3ripIBAABcokKEu+J47LHHlJWVpeHDh+v48eNq3ry51qxZIz8/P1eXBgAAUGYqbLhbv369w32bzabExEQlJia6pB4AAIDyoEJcUAEAAIDiIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALcXN1AQCsKSUlxdUlOChv9QCAsxDuAJSqrBPHJNnUr18/V5dSqHPZZ11dAgA4FeEOQKk6dzpTktG1fR9XzehYV5djd/i7Tdq1Yq5ycnJcXQoAOBXhDoBT+NaKUGBEjKvLsMs4vN/VJQBAmeCCCgAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAW4ubqAoCKIjU1Venp6a4uwy4lJcXVJQAAyiHCHVAMqampio2NU1bWaVeXUsC57LOuLgEAUI4Q7oBiSE9PV1bWaTUfnCD/0ChXlyNJOvzdJu1aMVc5OTmuLgUAUI4Q7oAS8A+NUmBEjKvLkCRlHN7v6hIAAOUQF1QAAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBA3VxcAXCg1NVXp6emuLsNBSkqKq0sAAKBYCHcoV1JTUxUbG6esrNOuLqVQ57LPuroEAAAuinCHciU9PV1ZWafVfHCC/EOjXF2O3eHvNmnXirnKyclxdSkAAFwU4Q7lkn9olAIjYlxdhl3G4f2uLgEAgGLhggoAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhFSLczZ49W/Hx8fL395e/v79atGihTz/91P64MUaJiYkKCwuTt7e32rVrp927d7uwYgAAANeoEOGuTp06eu6557Rt2zZt27ZNHTp0UK9evewBburUqZo2bZpmzpyprVu3KiQkRJ07d1ZmZqaLKwcAAChbFSLc9ezZU927d1eDBg3UoEEDTZw4Ub6+vtq8ebOMMZo+fbrGjRunPn36qFGjRlq4cKFOnz6txYsXu7p0AACAMlUhwt35cnNztXTpUp06dUotWrTQvn37lJaWpi5dutiX8fT0VNu2bZWcnOzCSgEAAMqem6sLKK7vvvtOLVq00JkzZ+Tr66sPPvhADRs2tAe44OBgh+WDg4N14MCBi64zOztb2dnZ9vsZGRmlXzgAAEAZqjBn7mJiYrRjxw5t3rxZw4YN04ABA/T999/bH7fZbA7LG2MKtF1o8uTJCggIsN/Cw8OdUjsAAEBZqTDhzsPDQ1deeaWaNWumyZMn65prrtFLL72kkJAQSVJaWprD8keOHClwNu9CY8eO1YkTJ+y3gwcPOq1+AACAslBhwt2FjDHKzs5WdHS0QkJClJSUZH/s7Nmz2rBhg1q2bHnRdXh6etp/XiX/BgAAUJFViO/cPfnkk+rWrZvCw8OVmZmppUuXav369Vq1apVsNptGjx6tSZMmqX79+qpfv74mTZokHx8f9e3b19WlAwAAlKkKEe5+//139e/fX4cPH1ZAQIDi4+O1atUqde7cWZL02GOPKSsrS8OHD9fx48fVvHlzrVmzRn5+fi6uHAAAoGxViHD3xhtvXPRxm82mxMREJSYmlk1BAAAA5VSF/c4dAAAACiLcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAAC6kQf6ECAKwsJSXF1SUUEBQUpIiICFeXAeAyEO4AwEWyThyTZFO/fv1cXUoB3t4+2rMnhYAHVECEOwBwkXOnMyUZXdv3cdWMjnV1OXYZh/dry7wJSk9PJ9wBFRDhDgBczLdWhAIjYlxdBgCL4IIKAAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAACyHcAQAAWIhTw13dunV17NixAu1//vmn6tat68yuAQAAKiWnhrv9+/crNze3QHt2drZ+++03Z3YNAABQKbk5Y6UrVqyw/3v16tUKCAiw38/NzdXatWsVFRXljK4BAAAqNaeEu969e0uSbDabBgwY4PCYu7u7oqKi9OKLLzqjawAAgErNKeEuLy9PkhQdHa2tW7cqKCjIGd0AAADgAk4Jd/n27dvnzNUDAADgAk4Nd5K0du1arV27VkeOHLGf0cs3b948Z3cPAABQqTg13E2YMEHPPPOMmjVrptDQUNlsNmd2BwAAUOk5NdzNmTNHCxYsUP/+/Z3ZDQAAAP6PU3/n7uzZs2rZsqUzuwAAAMB5nBru7r//fi1evNiZXQAAAOA8Tv1Y9syZM5o7d64+++wzxcfHy93d3eHxadOmObN7AACASsep4e7bb7/VtddeK0natWuXw2NcXAEAAFD6nBru1q1b58zVAwAA4AJO/c4dAAAAypZTz9y1b9/+oh+/fv75587sHgAAoNJxarjL/75dvnPnzmnHjh3atWuXBgwY4MyuAQAAKiWnhrv//Oc/hbYnJibq5MmTzuwaAACgUnLJd+769evH35UFAABwApeEu02bNsnLy8sVXQMAAFiaUz+W7dOnj8N9Y4wOHz6sbdu2afz48c7sGgAAoFJyargLCAhwuF+lShXFxMTomWeeUZcuXZzZNQAAQKXk1HA3f/58Z64eAAAAF3BquMu3fft2paSkyGazqWHDhmrcuHFZdAsAAFDpODXcHTlyRHfddZfWr1+vK664QsYYnThxQu3bt9fSpUtVs2ZNZ3YPAABQ6Tj1atl//etfysjI0O7du/XHH3/o+PHj2rVrlzIyMjRq1Chndg0AAFApOfXM3apVq/TZZ58pLi7O3tawYUO98sorXFABAADgBE49c5eXlyd3d/cC7e7u7srLy3Nm1wAAAJWSU8Ndhw4d9OCDD+rQoUP2tt9++00PPfSQOnbs6MyuAQAAKiWnhruZM2cqMzNTUVFRqlevnq688kpFR0crMzNTM2bMcGbXAAAAlZJTv3MXHh6u//73v0pKStKePXtkjFHDhg3VqVMnZ3YLAABQaTnlzN3nn3+uhg0bKiMjQ5LUuXNn/etf/9KoUaN03XXX6aqrrtLGjRud0TUAAECl5pRwN336dA0ZMkT+/v4FHgsICNDQoUM1bdo0Z3QNAABQqTkl3O3cuVM33XRTkY936dJF27dvd0bXAAAAlZpTwt3vv/9e6E+g5HNzc9PRo0ed0TUAAECl5pRwV7t2bX333XdFPv7tt98qNDTUGV0DAABUak4Jd927d9fTTz+tM2fOFHgsKytLCQkJ6tGjhzO6BgAAqNSc8lMoTz31lN5//301aNBAI0eOVExMjGw2m1JSUvTKK68oNzdX48aNc0bXAAAAlZpTztwFBwcrOTlZjRo10tixY3Xrrbeqd+/eevLJJ9WoUSN99dVXCg4OLvb6Jk+erOuuu05+fn6qVauWevfurb179zosY4xRYmKiwsLC5O3trXbt2mn37t2lPTQAAIByzWl/oSIyMlIrV65Uenq6tmzZos2bNys9PV0rV65UVFRUida1YcMGjRgxQps3b1ZSUpJycnLUpUsXnTp1yr7M1KlTNW3aNM2cOVNbt25VSEiIOnfurMzMzFIeGQAAQPnl1L9QIUnVq1fXdddd97fWsWrVKof78+fPV61atbR9+3a1adNGxhhNnz5d48aNU58+fSRJCxcuVHBwsBYvXqyhQ4f+rf4BAAAqCqf+bVlnOXHihCQpMDBQkrRv3z6lpaWpS5cu9mU8PT3Vtm1bJScnu6RGAAAAV3D6mbvSZozRmDFjdOONN6pRo0aSpLS0NEkq8D2+4OBgHThwoMh1ZWdnKzs7234//8+lAQAAVFQV7szdyJEj9e2332rJkiUFHrPZbA73jTEF2s43efJkBQQE2G/h4eGlXi8AAEBZqlDh7l//+pdWrFihdevWqU6dOvb2kJAQSf87g5fvyJEjF70qd+zYsTpx4oT9dvDgQecUDgAAUEYqRLgzxmjkyJF6//339fnnnys6Otrh8ejoaIWEhCgpKcnedvbsWW3YsEEtW7Yscr2enp7y9/d3uAEAAFRkFeI7dyNGjNDixYv10Ucfyc/Pz36GLiAgQN7e3rLZbBo9erQmTZqk+vXrq379+po0aZJ8fHzUt29fF1cPAABQdipEuJs9e7YkqV27dg7t8+fP18CBAyVJjz32mLKysjR8+HAdP35czZs315o1a+Tn51fG1QIAALhOhQh3xphLLmOz2ZSYmKjExETnFwQAAFBOVYjv3AEAAKB4CHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBC3FxdAACgfEpJSXF1CQ6CgoIUERHh6jKAco9wBwBwkHXimCSb+vXr5+pSHHh7+2jPnhQCHnAJhDsAgINzpzMlGV3b93HVjI51dTmSpIzD+7Vl3gSlp6cT7oBLINwBAArlWytCgRExri4DQAlxQQUAAICFEO4AAAAshHAHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAAC3FzdQEAABRXSkqKq0soICgoSBEREa4uA7Aj3AEAyr2sE8ck2dSvXz9Xl1KAt7eP9uxJIeCh3CDcAQDKvXOnMyUZXdv3cdWMjnV1OXYZh/dry7wJSk9PJ9yh3CDcAQAqDN9aEQqMiHF1GUC5xgUVAAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYCOEOAADAQipMuPviiy/Us2dPhYWFyWaz6cMPP3R43BijxMREhYWFydvbW+3atdPu3btdUywAAICLVJhwd+rUKV1zzTWaOXNmoY9PnTpV06ZN08yZM7V161aFhISoc+fOyszMLONKAQAAXMfN1QUUV7du3dStW7dCHzPGaPr06Ro3bpz69OkjSVq4cKGCg4O1ePFiDR06tCxLBQAAcJkKc+buYvbt26e0tDR16dLF3ubp6am2bdsqOTnZhZUBAACUrQpz5u5i0tLSJEnBwcEO7cHBwTpw4ECRz8vOzlZ2drb9fkZGhnMKBAAAKCOWOHOXz2azOdw3xhRoO9/kyZMVEBBgv4WHhzu7RAAAAKeyRLgLCQmR9L8zePmOHDlS4Gze+caOHasTJ07YbwcPHnRqnQAAAM5miXAXHR2tkJAQJSUl2dvOnj2rDRs2qGXLlkU+z9PTU/7+/g43AACAiqzCfOfu5MmT+umnn+z39+3bpx07digwMFAREREaPXq0Jk2apPr166t+/fqaNGmSfHx81LdvXxdWDQAAULYqTLjbtm2b2rdvb78/ZswYSdKAAQO0YMECPfbYY8rKytLw4cN1/PhxNW/eXGvWrJGfn5+rSgYAAChzFSbctWvXTsaYIh+32WxKTExUYmJi2RUFAABQzljiO3cAAAD4S4U5cwfnSE1NVXp6uqvLsEtJSXF1CQAAVGiEu0osNTVVsbFxyso67epSCjiXfdbVJQAAUCER7iqx9PR0ZWWdVvPBCfIPjXJ1OZKkw99t0q4Vc5WTk+PqUgAAqJAId5B/aJQCI2JcXYYkKePwfleXAABAhcYFFQAAABZCuAMAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAsxM3VBVQWqampSk9Pd3UZDlJSUlxdAgBYQnk7nmZnZ8vT09PVZRQQFBSkiIgIV5dheYS7MpCamqrY2DhlZZ12dSmFOpd91tUlAECFlHXimCSb+vXr5+pSHNlskjGurqIAb28f7dmTQsBzMsJdGUhPT1dW1mk1H5wg/9AoV5djd/i7Tdq1Yq5ycnJcXQoAVEjnTmdKMrq27+OqGR3r6nIk/e/YXp5qkqSMw/u1Zd4EpaenE+6cjHBXhvxDoxQYEePqMuwyDu93dQkAYAm+tSLKzfE9/9henmpC2eKCCgAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGWC3ezZs1SdHS0vLy81LRpU23cuNHVJQEAAJQZS4W7d955R6NHj9a4ceP0zTffqHXr1urWrZtSU1NdXRoAAECZsFS4mzZtmu677z7df//9iouL0/Tp0xUeHq7Zs2e7ujQAAIAyYZlwd/bsWW3fvl1dunRxaO/SpYuSk5NdVBUAAEDZcnN1AaUlPT1dubm5Cg4OdmgPDg5WWlpaoc/Jzs5Wdna2/f6JEyckSRkZGaVa28mTJyVJfxzYq5zsrFJd99+RcfiAJOnEbz/K3c3m4mr+Uh5rkspnXeWxJom6SqI81iSVz7rKY01S+ayrPNYkSRlpf31F6uTJk6X+Ppu/PmNMqa63wjIW8dtvvxlJJjk52aH93//+t4mJiSn0OQkJCUYSN27cuHHjxs0Ct4MHD5ZF5Cj3LHPmLigoSFWrVi1wlu7IkSMFzublGzt2rMaMGWO/n5eXpz/++EM1atSQzWZTRkaGwsPDdfDgQfn7+zu1flwe5qhiYJ7KP+aoYmCeCmeMUWZmpsLCwlxdSrlgmXDn4eGhpk2bKikpSbfeequ9PSkpSb169Sr0OZ6envL09HRou+KKKwos5+/vz4uonGOOKgbmqfxjjioG5qmggIAAV5dQblgm3EnSmDFj1L9/fzVr1kwtWrTQ3LlzlZqaqgceeMDVpQEAAJQJS4W7O++8U8eOHdMzzzyjw4cPq1GjRlq5cqUiIyNdXRoAAECZsFS4k6Thw4dr+PDhpbIuT09PJSQkFPjoFuUHc1QxME/lH3NUMTBPKA6bMVw3DAAAYBWW+RFjAAAAEO4AAAAshXAHAABgIYQ7AAAAC6lU4W7WrFmKjo6Wl5eXmjZtqo0bN150+VdeeUVxcXHy9vZWTEyM3nzzzQLLLF++XA0bNpSnp6caNmyoDz74wFnlVxqlPU+vvfaaWrdurerVq6t69erq1KmTvv76a2cOwfKc8VrKt3TpUtlsNvXu3buUq658nDFPf/75p0aMGKHQ0FB5eXkpLi5OK1eudNYQLM8ZczR9+nTFxMTI29tb4eHheuihh3TmzBlnDQHlkYv//FmZWbp0qXF3dzevvfaa+f77782DDz5oqlWrZg4cOFDo8rNmzTJ+fn5m6dKl5ueffzZLliwxvr6+ZsWKFfZlkpOTTdWqVc2kSZNMSkqKmTRpknFzczObN28uq2FZjjPmqW/fvuaVV14x33zzjUlJSTGDBg0yAQEB5tdffy2rYVmKM+Yo3/79+03t2rVN69atTa9evZw8EmtzxjxlZ2ebZs2ame7du5svv/zS7N+/32zcuNHs2LGjrIZlKc6Yo7ffftt4enqaRYsWmX379pnVq1eb0NBQM3r06LIaFsqBShPurr/+evPAAw84tMXGxponnnii0OVbtGhhHnnkEYe2Bx980LRq1cp+/4477jA33XSTwzJdu3Y1d911VylVXfk4Y54ulJOTY/z8/MzChQv/fsGVkLPmKCcnx7Rq1cq8/vrrZsCAAYS7v8kZ8zR79mxTt25dc/bs2dIvuBJyxhyNGDHCdOjQwWGZMWPGmBtvvLGUqkZFUCk+lj179qy2b9+uLl26OLR36dJFycnJhT4nOztbXl5eDm3e3t76+uuvde7cOUnSpk2bCqyza9euRa4TF+esebrQ6dOnde7cOQUGBpZO4ZWIM+fomWeeUc2aNXXfffeVfuGVjLPmacWKFWrRooVGjBih4OBgNWrUSJMmTVJubq5zBmJhzpqjG2+8Udu3b7d/9eSXX37RypUrdfPNNzthFCivKkW4S09PV25uroKDgx3ag4ODlZaWVuhzunbtqtdff13bt2+XMUbbtm3TvHnzdO7cOaWnp0uS0tLSSrROXJyz5ulCTzzxhGrXrq1OnTqV+hiszllz9NVXX+mNN97Qa6+95vQxVAbOmqdffvlF7733nnJzc7Vy5Uo99dRTevHFFzVx4kSnj8lqnDVHd911l5599lndeOONcnd3V7169dS+fXs98cQTTh8Tyo9KEe7y2Ww2h/vGmAJt+caPH69u3brphhtukLu7u3r16qWBAwdKkqpWrXpZ60TxOGOe8k2dOlVLlizR+++/X+B/wCi+0pyjzMxM9evXT6+99pqCgoKcXXqlUtqvpby8PNWqVUtz585V06ZNddddd2ncuHGaPXu2U8dhZaU9R+vXr9fEiRM1a9Ys/fe//9X777+vTz75RM8++6xTx4HypVKEu6CgIFWtWrXA/4aOHDlS4H9N+by9vTVv3jydPn1a+/fvV2pqqqKiouTn52d/AwoJCSnROnFxzpqnfC+88IImTZqkNWvWKD4+3mnjsDJnzNHPP/+s/fv3q2fPnnJzc5Obm5vefPNNrVixQm5ubvr555/LYmiW4qzXUmhoqBo0aODwH6e4uDilpaXp7NmzzhuQBTlrjsaPH6/+/fvr/vvv19VXX61bb71VkyZN0uTJk5WXl+f0caF8qBThzsPDQ02bNlVSUpJDe1JSklq2bHnR57q7u6tOnTqqWrWqli5dqh49eqhKlb82W4sWLQqsc82aNZdcJwrnrHmSpOeff17PPvusVq1apWbNmjml/srAGXMUGxur7777Tjt27LDfbrnlFrVv3147duxQeHi4M4dkSc56LbVq1Uo//fSTQ0j44YcfFBoaKg8Pj9IfiIU5a45Onz7tcOyT/jqrZ/66gLJ0B4HyyyWXcbhA/iXnb7zxhvn+++/N6NGjTbVq1cz+/fuNMcY88cQTpn///vbl9+7da9566y3zww8/mC1btpg777zTBAYGmn379tmX+eqrr0zVqlXNc889Z1JSUsxzzz3HT6H8Tc6YpylTphgPDw/z3nvvmcOHD9tvmZmZZT08S3DGHF2Iq2X/PmfMU2pqqvH19TUjR440e/fuNZ988ompVauW+fe//13Ww7MEZ8xRQkKC8fPzM0uWLDG//PKLWbNmjalXr5654447ynp4cKFKE+6MMeaVV14xkZGRxsPDwzRp0sRs2LDB/tiAAQNM27Zt7fe///57c+211xpvb2/j7+9vevXqZfbs2VNgncuWLTMxMTHG3d3dxMbGmuXLl5fFUCyttOcpMjLSSCpwS0hIKKMRWY8zXkvnI9yVDmfMU3JysmnevLnx9PQ0devWNRMnTjQ5OTllMRxLKu05OnfunElMTDT16tUzXl5eJjw83AwfPtwcP368jEaE8sBmDOdpAQAArKJSfOcOAACgsiDcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAshHAHAC6UmJioa6+9ttysB0DFR7gDKrAjR45o6NChioiIkKenp0JCQtS1a1dt2rTJ1aWVW8uXL1e7du0UEBAgX19fxcfH65lnntEff/zh6tKKzWaz6cMPP3Roe+SRR7R27VrXFASgXCHcARXYbbfdpp07d2rhwoX64YcftGLFCrVr165cBJVz5865uoQCxo0bpzvvvFPXXXedPv30U+3atUsvvviidu7cqbfeeuuy11vYWMt6/L6+vqpRo0aZ9gmgnHL13z8DcHmOHz9uJJn169dfdDlJZtasWeamm24yXl5eJioqyrz77rsOy/z666/mjjvuMFdccYUJDAw0t9xyi8MfI//6669Np06dTI0aNYy/v79p06aN2b59e4F+Zs+ebW655Rbj4+Njnn76aZOQkGCuueYa88Ybb5jw8HBTrVo188ADD5icnBwzZcoUExwcbGrWrFngD8+/+OKLplGjRsbHx8fUqVPHDBs2zGRmZtofnz9/vgkICDCrVq0ysbGxplq1aqZr167m0KFDRW6HLVu2GElm+vTpRW7PfLNmzTJ169Y17u7upkGDBubNN98s0Vijo6ONzWYzeXl55s8//zRDhgwxNWvWNH5+fqZ9+/Zmx44d9nXlP6+42/rCv5UcGRlZ6Hpyc3PNhAkTTO3atY2Hh4e55pprzKeffmp/fN++fUaSWb58uWnXrp3x9vY28fHxJjk5uchtCKBi4MwdUEH5+vrK19dXH374obKzsy+67Pjx4+1n+fr166e7775bKSkpkqTTp0+rffv28vX11RdffKEvv/xSvr6+uummm3T27FlJUmZmpgYMGKCNGzdq8+bNql+/vrp3767MzEyHfhISEtSrVy999913Gjx4sCTp559/1qeffqpVq1ZpyZIlmjdvnm6++Wb9+uuv2rBhg6ZMmaKnnnpKmzdvtq+nSpUqevnll7Vr1y4tXLhQn3/+uR577DGHvk6fPq0XXnhBb731lr744gulpqbqkUceKXIbLFq0SL6+vho+fHihj19xxRWSpA8++EAPPvigHn74Ye3atUtDhw7VoEGDtG7dukuO9aefftK7776r5cuXa8eOHZKkm2++WWlpaVq5cqW2b9+uJk2aqGPHjkWeXb3Utt66daskaf78+Tp8+LD9/oVeeuklvfjii3rhhRf07bffqmvXrrrlllv0448/Oiw3btw4PfLII9qxY4caNGigu+++Wzk5OUVuRwAVgKvTJYDL995775nq1asbLy8v07JlSzN27Fizc+dOh2UkmQceeMChrXnz5mbYsGHGGGPeeOMNExMTY/Ly8uyPZ2dnG29vb7N69epC+83JyTF+fn7m448/duhn9OjRDsslJCQYHx8fk5GRYW/r2rWriYqKMrm5ufa2mJgYM3ny5CLH+e6775oaNWrY78+fP99IMj/99JO97ZVXXjHBwcFFrqNbt24mPj6+yMfztWzZ0gwZMsSh7fbbbzfdu3e33y9qrO7u7ubIkSP2trVr1xp/f39z5swZh2Xr1atnXn31Vfvzzj/jdqGitvUHH3xQoP/z1xMWFmYmTpzosMx1111nhg8fboz535m7119/3f747t27jSSTkpJSZD0Ayj/O3AEV2G233aZDhw5pxYoV6tq1q9avX68mTZpowYIFDsu1aNGiwP38M3fbt2/XTz/9JD8/P/vZwMDAQJ05c0Y///yzpL8u3HjggQfUoEEDBQQEKCAgQCdPnlRqaqrDeps1a1agxqioKPn5+dnvBwcHq2HDhqpSpYpD25EjR+z3161bp86dO6t27dry8/PTvffeq2PHjunUqVP2ZXx8fFSvXj37/dDQUId1XMgYI5vNVuTj+VJSUtSqVSuHtlatWtm318XGGhkZqZo1a9rvb9++XSdPnlSNGjXs29bX11f79u2zb9sLFXdbX0xGRoYOHTpUrHHEx8fb/x0aGmqvAUDF5ebqAgD8PV5eXurcubM6d+6sp59+Wvfff78SEhI0cODAiz4vP+jk5eWpadOmWrRoUYFl8oPKwIEDdfToUU2fPl2RkZHy9PRUixYt7B/b5qtWrVqBdbi7uxfot7C2vLw8SdKBAwfUvXt3PfDAA3r22WcVGBioL7/8Uvfdd5/DRQqFrcMYU+R4GzRooC+//FLnzp0r8NwLXRgCCwuGhY31wra8vDyFhoZq/fr1BZbN/xj4QsXd1sVRnHGcvy3O3ycAVFycuQMspmHDhg5nuCQ5fJ8t/35sbKwkqUmTJvrxxx9Vq1YtXXnllQ63gIAASdLGjRs1atQode/eXVdddZU8PT2Vnp7ulPq3bdumnJwcvfjii7rhhhvUoEEDHTp06G+vt2/fvjp58qRmzZpV6ON//vmnJCkuLk5ffvmlw2PJycmKi4srcZ9NmjRRWlqa3NzcCmzboKCgQp9TnG3t7u6u3NzcIvv19/dXWFhYqY0DQMVCuAMqqGPHjqlDhw56++239e2332rfvn1atmyZpk6dql69ejksu2zZMs2bN08//PCDEhIS9PXXX2vkyJGSpHvuuUdBQUHq1auXNm7cqH379mnDhg168MEH9euvv0qSrrzySr311ltKSUnRli1bdM8998jb29sp46pXr55ycnI0Y8YM/fLLL3rrrbc0Z86cv73e5s2b67HHHtPDDz+sxx57TJs2bdKBAwe0du1a3X777Vq4cKEk6dFHH9WCBQs0Z84c/fjjj5o2bZref//9i16sUZROnTqpRYsW6t27t1avXq39+/crOTlZTz31lLZt21boc4qzraOiorR27VqlpaXp+PHjha7n0Ucf1ZQpU/TOO+9o7969euKJJ7Rjxw49+OCDJR4HgIqFcAdUUL6+vmrevLn+85//qE2bNmrUqJHGjx+vIUOGaObMmQ7LTpgwQUuXLlV8fLwWLlyoRYsWqWHDhpL++u7aF198oYiICPXp00dxcXEaPHiwsrKy5O/vL0maN2+ejh8/rsaNG6t///4aNWqUatWq5ZRxXXvttZo2bZqmTJmiRo0aadGiRZo8eXKprHvKlClavHixtmzZoq5du+qqq67SmDFjFB8frwEDBkiSevfurZdeeknPP/+8rrrqKr366quaP3++2rVrV+L+bDabVq5cqTZt2mjw4MFq0KCB7rrrLu3fv1/BwcGFPqc42/rFF19UUlKSwsPD1bhx40LXM2rUKD388MN6+OGHdfXVV2vVqlVasWKF6tevX+JxAKhYbOZiX1IBUOHZbDZ98MEH6t27t6tLAQCUAc7cAQAAWAjhDgAAwEL4KRTA4vjmBQBULpy5AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsJD/D6CgL1YTvUBvAAAAAElFTkSuQmCC",
      "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": 7,
   "id": "4b36f1b2-1b9c-4686-b9fb-43c1d9f64b59",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(\n",
    "    \"figures/reconstruction_spearmans/scETM_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",
    "# aggregate cell embeddings to perturbation embeddings\n",
    "cell_emb = adata.uns[\"cell_emb\"] @ adata.uns[\"topics\"]\n",
    "perts = []\n",
    "pert_emb = []\n",
    "for i in adata.obs[\"condition\"].unique():\n",
    "    if \"+\" not in i and i != \"nan\":\n",
    "        perts.append(i)\n",
    "        pert_emb.append(cell_emb[adata.obs[\"condition\"] == i].mean(axis=0))\n",
    "pert_emb = np.array(pert_emb)\n",
    "pert_emb_df = pd.DataFrame(pert_emb, index=perts)\n",
    "\n",
    "# Compute the pairwise distances\n",
    "df = pert_emb_df.drop(index=[\"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": 9,
   "id": "271cd67a-9441-471b-b125-819e8b7b738a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi4AAAHMCAYAAADyLiPLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABnZ0lEQVR4nO3dd1QU198G8O/SQQWxEVHEjghGsaNgRRSNvSUmauxYYu/YosYexSSWGGMsSYwa0/xFE00xGjWJBWxgjYIFJdiN1OV5/+DdyS7NXVjKyPM5hyMOs7N3d+/OPHPnzr0aABAiIiIiFbAo6AIQERERGYvBhYiIiFSDwYWIiIhUg8GFiIiIVIPBhYiIiFSDwYWIiIhUg8GFiIiIVIPBhYiIiFTDqqALYIzU1FS5ffu2lChRQjQaTUEXh4iIiIwAQJ48eSKurq5iYWGethJVBJfbt2+Lm5tbQReDiIiIcuDGjRtSsWJFs2xLFcGlRIkSIpL2wh0dHQu4NERERGSMx48fi5ubm3IcNwdVBBfd5SFHR0cGFyIiIpUxZzcPds4lIiIi1WBwISIiItXIUXBZu3atVKlSRezs7KRBgwZy+PDhLNc9ePCgaDSaDD8XLlzIcaGJiIioaDI5uOzYsUPGjx8vISEhEhYWJv7+/hIUFCTR0dHZPu7ixYsSExOj/NSoUSPHhSYiIqKiyeTgsnLlShkyZIgMHTpUPD09JTQ0VNzc3GTdunXZPq5cuXLy0ksvKT+WlpY5LjQREREVTSYFl6SkJDl58qQEBgYaLA8MDJSjR49m+1gfHx8pX768tG3bVn799dds101MTJTHjx8b/BARERGZFFzi4uJEq9WKi4uLwXIXFxe5c+dOpo8pX768bNiwQXbv3i1fffWVeHh4SNu2beXQoUNZPs/ixYvFyclJ+eHgc0RERCSSw3Fc0t+PDSDLe7Q9PDzEw8ND+b+vr6/cuHFDVqxYIS1atMj0MTNmzJCJEycq/9cNYENERERFm0ktLmXKlBFLS8sMrSuxsbEZWmGy07RpU7l8+XKWf7e1tVUGm+Ogc0RERKRjUnCxsbGRBg0ayIEDBwyWHzhwQJo1a2b0dsLCwqR8+fKmPDURERGR6ZeKJk6cKP3795eGDRuKr6+vbNiwQaKjoyU4OFhE0i7z3Lp1S7Zu3SoiIqGhoVK5cmXx8vKSpKQk+fTTT2X37t2ye/du874SIiIieuGZHFz69u0r9+7dk/nz50tMTIx4e3vL3r17xd3dXUREYmJiDMZ0SUpKksmTJ8utW7fE3t5evLy85Pvvv5eOHTua71UQERFRkaABgIIuxPM8fvxYnJyc5NGjR+zvQkREpBJ5cfxWxezQVLAASHyytqCLQURU4OytLc060zGZjsGFsgVAeq0/JiejHhR0UYiIClxDd2fZFezL8FKAODs0ZSs+WcvQQkT0/05EPWALdAFjiwsZ7cSsAHGw4RxTRFT0PEvSSsOFPxV0MUgYXMgEDjaW4mDDKkNERAWHl4qIiIhINRhciIiISDUYXIiIiEg1GFyIiIhINRhciIiISDUYXIiIiEg1GFyIiIhINRhciIiISDUYXIiIiEg1GFyIiIhINRhciIiISDUYXIiIiEg1GFyIiIhINRhciIiISDUYXIiIiEg1GFyIiIhINRhciIiISDUYXIiIiEg1GFyIiIhINRhciIiISDUYXIiIiEg1GFyIiIhINRhciIiISDUYXIiIiEg1GFyIiIhINRhciIiISDUYXIiIiEg1GFyIiIhINRhciIiISDUYXIiIiEg1GFyIiIhINRhciIiISDUYXIiIiEg1GFyIiIhINRhciIiISDUYXIiIiEg1GFyIiIhINRhciIiISDUYXIiIiEg1GFyIiIhINRhciIiISDUYXIiIiEg1GFyIiIhINRhciIiISDUYXIiIiEg1GFyIiIhINRhciIiISDUYXIiIiEg1GFyIiIhINRhciIiISDUYXIiIiEg1GFyIiIhINRhciIiISDUYXIiIiEg1GFyIiIhINRhciIiISDVyFFzWrl0rVapUETs7O2nQoIEcPnzYqMcdOXJErKyspF69ejl5WiIiIiriTA4uO3bskPHjx0tISIiEhYWJv7+/BAUFSXR0dLaPe/TokQwYMEDatm2b48ISERFR0WZycFm5cqUMGTJEhg4dKp6enhIaGipubm6ybt26bB83YsQI6devn/j6+ua4sERERFS0mRRckpKS5OTJkxIYGGiwPDAwUI4ePZrl4z755BO5evWqzJ0716jnSUxMlMePHxv8EBEREZkUXOLi4kSr1YqLi4vBchcXF7lz506mj7l8+bJMnz5dPvvsM7GysjLqeRYvXixOTk7Kj5ubmynFJCIiohdUjjrnajQag/8DyLBMRESr1Uq/fv3k7bfflpo1axq9/RkzZsijR4+Unxs3buSkmERERPSCMa4J5P+VKVNGLC0tM7SuxMbGZmiFERF58uSJnDhxQsLCwmTMmDEiIpKamioAxMrKSvbv3y9t2rTJ8DhbW1uxtbU1pWhERERUBJjU4mJjYyMNGjSQAwcOGCw/cOCANGvWLMP6jo6OcvbsWQkPD1d+goODxcPDQ8LDw6VJkya5Kz0REREVKSa1uIiITJw4Ufr37y8NGzYUX19f2bBhg0RHR0twcLCIpF3muXXrlmzdulUsLCzE29vb4PHlypUTOzu7DMuJiIiInsfk4NK3b1+5d++ezJ8/X2JiYsTb21v27t0r7u7uIiISExPz3DFdiIiIiHJCAwAFXYjnefz4sTg5OcmjR4/E0dGxoItTpDxLSpHac34UEZGI+e3FwcbkrEtEpHrcF+ZMXhy/OVcRERERqQaDCxEREakGgwsRERGpBoMLERERqQaDCxEREakGgwsRERGpBoMLERERqQaDCxEREakGgwsRERGpBoMLERERqQbHLCYiItUBIPHJ2nx7vmdJKZn+nl/srS1Fo9Hk+/MWRgwuRESkKgCk1/pjcjLqQYE8f8OFP+f/c7o7y65gX4YX4aUiIiJSmfhkbYGFloJyIupBvrYwFWZscSEiItU6MStAHGwsC7oYeeZZklYaLvypoItRqDC4EBGRajnYWIqDDQ9lRQkvFREREZFqMLgQERGRajC4EBERkWowuBAREZFqMLgQERGRajC4EBERkWowuBAREZFqMLgQERGRajC4EBERkWowuBAREZFqMLgQERGRajC4EBERkWowuBAREZFqMLgQERGRajC4EBERkWowuBAREZFqMLgQERGRajC4EBERkWowuBAREZFqMLgQERGRajC4EBERkWowuBAREZFqMLgQERGRajC4EBERkWowuBAREZFqMLgQERGRajC4EBERkWowuBAREZFqMLgQERGRajC4EBERkWowuBAREZFqMLgQERGRajC4EBERkWowuBAREZFqMLgQERGRajC4EBERkWowuBAREZFqMLgQERGRajC4EBERkWowuBAREZFqMLgQERGRajC4EBERkWowuBAREZFqWBV0AfICAIlP1hZ0MV4Iz5JSMv2dcsfe2lI0Gk1BF4OISHVeuOACQHqtPyYnox4UdFFeOA0X/lzQRXhhNHR3ll3BvgwvREQmeuEuFcUnaxlaqNA7EfWArYJERDmQoxaXtWvXyvLlyyUmJka8vLwkNDRU/P39M133999/l2nTpsmFCxfk2bNn4u7uLiNGjJAJEybkquDGODErQBxsLPP8eYiM9SxJKw0X/lTQxSAiUi2Tg8uOHTtk/PjxsnbtWmnevLl8+OGHEhQUJBEREVKpUqUM6xcrVkzGjBkjL7/8shQrVkx+//13GTFihBQrVkyGDx9ulheRFQcbS3GweeGuhhERERVZJl8qWrlypQwZMkSGDh0qnp6eEhoaKm5ubrJu3bpM1/fx8ZHXXntNvLy8pHLlyvLGG29I+/bt5fDhw7kuPBERERUtJgWXpKQkOXnypAQGBhosDwwMlKNHjxq1jbCwMDl69Ki0bNkyy3USExPl8ePHBj9EREREJgWXuLg40Wq14uLiYrDcxcVF7ty5k+1jK1asKLa2ttKwYUMZPXq0DB06NMt1Fy9eLE5OTsqPm5ubKcUkIiKiF1SO7ipKfwsngOfe1nn48GE5ceKErF+/XkJDQ2X79u1Zrjtjxgx59OiR8nPjxo2cFJOIiIheMCb1XC1TpoxYWlpmaF2JjY3N0AqTXpUqVUREpE6dOnL37l2ZN2+evPbaa5mua2trK7a2tqYUjYiIiIoAk1pcbGxspEGDBnLgwAGD5QcOHJBmzZoZvR0AkpiYaMpTExEREZl+O/TEiROlf//+0rBhQ/H19ZUNGzZIdHS0BAcHi0jaZZ5bt27J1q1bRURkzZo1UqlSJalVq5aIpI3rsmLFCnnrrbfM+DKIiIioKDA5uPTt21fu3bsn8+fPl5iYGPH29pa9e/eKu7u7iIjExMRIdHS0sn5qaqrMmDFDrl27JlZWVlKtWjVZsmSJjBgxwnyvgoiIiIqEHI3ONmrUKBk1alSmf9u8ebPB/9966y22rhAREZFZvHBzFREREdGLi8GFiIiIVIPBhYiIiFSDwYWIiIhUg8GFiIiIVIPBhYiIiFSDwYWIiIhUg8GFiIiIVIPBhYiIiFSDwYWIiIhUg8GFiIiIVIPBhYiIiFSDwYWIiIhUg8GFiIiIVIPBhYiIiFTDqqALQEREpGYAJD5ZmyfbfpaUkunv5mRvbSkajSZPtp0XGFyIiIhyCID0Wn9MTkY9yPPnarjw57zZrruz7Ar2VU144aUiIiKiHIpP1uZLaMlLJ6Ie5FmLUV5giwsREZEZnJgVIA42lgVdDKM9S9JKw4U/FXQxTMbgQkREZAYONpbiYMPDal7jpSIiIiJSDQYXIiIiUg0GFyIiIlINBhciIiJSDQYXIiIiUg0GFyIiIlIN3rdVhOTlsNRknPwYvptMo7bhzomKOgaXIiI/h6Um4+TV8N1kGrUNd05U1PFSURHxIgxLTZQX1DbcOVFRxxaXIkhtw1IT5QW1DndOVNQxuBRBHJaaiIjUipeKiIiISDUYXIiIiEg1GFyIiIhINRhciIiISDUYXIiIiEg1GFyIiIhINRhciIiISDUYXIiIiEg1GFyIiIhINRhciIiISDUYXIiIiEg1GFyIiIhINRhciIiISDUYXIiIiEg1GFyIiIhINRhciIiISDWsCroARFS0AZD4lPh8f95nyVq93+NFNJb5XgZ7K3vRaDT5/rxEasbgQkQFBoAM2DdAwv8Jz//nTrUWkQUiItJqZ0vRWCTnexl8yvnIlg5bGF6ITMDgQkQFJj4lvkBCi4iIxiJZSnhOL5Dn1gmLDZP4lHhxsHYo0HIQqQmDCxEVCgf7HBR7K/uCLka+iE+Jl1Y7WxV0MYhUicGFiAoFeyt7tjwQ0XPxriIiIiJSDQYXIiIiUg0GFyIiIlINBhciIiJSDXbONRUgkvysoEthuiSt3u/PRCT/B9vKFWsHEY51QURU5DG4mAIQ2dRe5MafBV0S08FWRD5J+315dRFNYoEWx2QVG4n0/+bFCS8MYkREOcLgYorkZ+oMLSLioEmU63b9CroYOXfzuMjiCgVdCvNxayoy+AeGFyIiEzG45NTkKyI2hXTMiaRnIiuqF3QpKDs3/kgLwjbFCrokRESqwuCSUzYO6jjoFOaAVRQxVBIR5UqOgsvatWtl+fLlEhMTI15eXhIaGir+/v6ZrvvVV1/JunXrJDw8XBITE8XLy0vmzZsn7du3z1XByUhqCVhERERGMPl26B07dsj48eMlJCREwsLCxN/fX4KCgiQ6OjrT9Q8dOiTt2rWTvXv3ysmTJ6V169bSuXNnCQsLy3XhiYiIqGgxObisXLlShgwZIkOHDhVPT08JDQ0VNzc3WbduXabrh4aGytSpU6VRo0ZSo0YNWbRokdSoUUP27NmT68ITERFR0WJScElKSpKTJ09KYGCgwfLAwEA5evSoUdtITU2VJ0+eSKlSpbJcJzExUR4/fmzwQ0RERGRSH5e4uDjRarXi4uJisNzFxUXu3Llj1Dbeffdd+ffff6VPnz5ZrrN48WJ5++23TSkakXEKegDBpGeZ/57fOI4MEalUjjrnatLt8ABkWJaZ7du3y7x58+Tbb7+VcuXKZbnejBkzZOLEicr/Hz9+LG5ubjkpKtF/CtsAggV5dxHHkSEilTIpuJQpU0YsLS0ztK7ExsZmaIVJb8eOHTJkyBDZtWuXBAQEZLuura2t2NramlI0oudT8QCCZsdxZIhIpUwKLjY2NtKgQQM5cOCAdO/eXVl+4MAB6dq1a5aP2759uwwePFi2b98unTp1ynlpicylqI5vw3FkiEjlTL5UNHHiROnfv780bNhQfH19ZcOGDRIdHS3BwcEiknaZ59atW7J161YRSQstAwYMkNWrV0vTpk2V1hp7e3txcnIy40shMgHHt6FMAJD4lPg8fx7958iP57O3sjfqcj6RGpgcXPr27Sv37t2T+fPnS0xMjHh7e8vevXvF3d1dRERiYmIMxnT58MMPJSUlRUaPHi2jR49Wlg8cOFA2b96c+1dARGQGAGTAvgES/k94vj5vq52t8vw5fMr5yJYOWxhe6IWQo865o0aNklGjRmX6t/Rh5ODBgzl5CiKifBWfEp/voSW/hMWGSXxKvDhYF8HLo/TC4VxFRETpHOxzUOyt7Au6GLkWnxKfLy06RPmJwaWwMcc4I+YcK4TjfVARZG9lz9YJokKKwaUwyYtxRnJ7BwnH+yAiokLE5LmKKA8VxnFGdON9EBERFQJscSmsCnqcEY73QUREhRCDS2HFcUaIiIgy4KUiIiIiUg0GFyIiIlINBhciIiJSDQYXIiIiUg0GFyIiIlIN3lVU1Bg7Mm9OR9/lSLtERJSHGFyKkpyOzGvKeC4caZeIiPIQLxUVJfkxMi9H2iUiojzEFpeiytwj83KkXSIiygcMLkUVR+YlIiIV4qUiIiIiUg0GFyIiIlINXioiIiqkAEh8SnyOH6//2NxsR8feyl40vGOQChiDCxFRIQRABuwbIOH/hJtle612tsr1NnzK+ciWDlsYXqhA8VIREVEhFJ8Sb7bQYi5hsWFmabkhyg22uBARFXIH+xwUeyv7Anv++JR4s7TYEJkDgwsRUSFnb2UvDtZmHHeJSMV4qYiIiIhUg8GFiIiIVIPBhYiIiFSDfVyIXgSAcZNbJj3L/PfnsXbgjN9EVCgwuBCpHSCyqb3pM3+bMimmW1ORwT8wvBBRgeOlIiK1S35memgx1Y0/jGvRISLKY2xxIXqRTL6SNvO3uSQ9M61lhogojzG4EL1IbBxEbIoVdCmIiPIMgwsRvbBMmaQwpxMScuJBovzF4EJEL6TcTFJoyvD2nHiQKH+xcy4RvZDya5JCTjxIlL/Y4kJEL7y8mKSQEw8SFQwGFyJ64XGSQqIXB4MLERHRCwKApCSlGrVucpL2v98TtZIM4/ppWdlYFGifLgYXIiIiPQAE8cb1W0rVO/inPouX1BRLox6nsTf/3WgA5Kvlp+TO34+MWj9JIFIy7fdNU34XGzGuPOWrOUn3yfULLLwwuJB6GDsfT1ZyOk9PVjh/D9ELB4BE9Xtd4sPCjFo/wdJGpPMiERG53NxP7LRJRj3Ovn59cf/sU7Me/FOSUo0OLSIiNqKRKQ9N7/sVc/WRpCSlirWtcSHN3BhcSB1yOh9PVswxGizn7yF64SA+3ujQIiJip02Sfd9MNvl54k+dEsTHi8Yhb/peDVrmZ/ZgkZyolU+m/m7WbeYEgwupQ37Mx2Mq3fw9HKmW6IVU48jvYmFv3rvRUuPj5XJzP7NuMzPWtpYF1iKS1xhcSH3MPR+PqTh/D1GRYGFvLxZ51CJCOcfgQurD+XiIiIosjpxLREREqsEWFyKiIogTUJJaMbgQERUxnICS1IyXioiIihhOQElqxhYXIqIijBNQktowuBARFWGcgJLUhpeKiIiISDUYXIiIiEg1GFyIiIhINRhciIiISDUYXIiIiEg1GFyIiIhINRhciIiISDU4jguRiAggkvzMuHWTnmX++/NYO4hw6HMiolxhcCECRDa1F7nxp+mPXVHd+HXdmooM/oHhhYgoF3ipiCj5Wc5Ci6lu/GF8qw4REWWKLS5E+iZfEbEx8/DnSc9Ma5khIqIsMbgQ6bNxELEpVtClICKiLOToUtHatWulSpUqYmdnJw0aNJDDhw9nuW5MTIz069dPPDw8xMLCQsaPH5/TshIREVERZ3Jw2bFjh4wfP15CQkIkLCxM/P39JSgoSKKjozNdPzExUcqWLSshISFSt27dXBeYiIiIii6Tg8vKlStlyJAhMnToUPH09JTQ0FBxc3OTdevWZbp+5cqVZfXq1TJgwABxcnLKdYGJiIio6DIpuCQlJcnJkyclMDDQYHlgYKAcPXrUbIVKTEyUx48fG/wQERERmRRc4uLiRKvViouLi8FyFxcXuXPnjtkKtXjxYnFyclJ+3NzczLZtIiIiUq8cdc7VpBtAC0CGZbkxY8YMefTokfJz48YNs22biIiI1Muk26HLlCkjlpaWGVpXYmNjM7TC5Iatra3Y2tqabXtERET0YjCpxcXGxkYaNGggBw4cMFh+4MABadasmVkLRkRERJSeyQPQTZw4Ufr37y8NGzYUX19f2bBhg0RHR0twcLCIpF3muXXrlmzdulV5THh4uIiIPH36VP755x8JDw8XGxsbqV27tnleBRERERUJJgeXvn37yr1792T+/PkSExMj3t7esnfvXnF3dxeRtAHn0o/p4uPjo/x+8uRJ+fzzz8Xd3V2uX7+eu9ITERFRkZKjIf9HjRolo0aNyvRvmzdvzrAMQE6ehoiIiMgAZ4cmIiIi1WBwISIiItVgcCEiIiLVYHAhIiIi1WBwISIiItVgcCEiIiLVYHAhIiIi1WBwISIiItVgcCEiIiLVYHAhIiIi1WBwISIiItVgcCEiIiLVYHAhIiIi1WBwISIiItVgcCEiIiLVYHAhIiIi1WBwISIiItVgcCEiIiLVYHAhIiIi1WBwISIiItVgcCEiIiLVYHAhIiIi1WBwISIiItVgcCEiIiLVYHAhIiIi1WBwISIiItVgcCEiIiLVYHAhIiIi1WBwISIiItVgcCEiIiLVYHAhIiIi1WBwISIiItVgcCEiIiLVYHAhIiIi1WBwISIiItVgcCEiIiLVYHAhIiIi1WBwISIiItVgcCEiIiLVYHAhIiIi1WBwISIiItVgcCEiIiLVYHAhIiIi1WBwISIiItVgcCEiIiLVYHAhIiIi1WBwISIiItVgcCEiIiLVYHAhIiIi1WBwISIiItVgcCEiIiLVYHAhIiIi1WBwISIiItVgcCEiIiLVYHAhIiIi1WBwISIiItVgcCEiIiLVYHAhIiIi1WBwISIiItVgcCEiIiLVYHAhIiIi1chRcFm7dq1UqVJF7OzspEGDBnL48OFs1//tt9+kQYMGYmdnJ1WrVpX169fnqLBERERUtJkcXHbs2CHjx4+XkJAQCQsLE39/fwkKCpLo6OhM17927Zp07NhR/P39JSwsTGbOnCljx46V3bt357rwREREVLSYHFxWrlwpQ4YMkaFDh4qnp6eEhoaKm5ubrFu3LtP1169fL5UqVZLQ0FDx9PSUoUOHyuDBg2XFihW5LjwREREVLSYFl6SkJDl58qQEBgYaLA8MDJSjR49m+phjx45lWL99+/Zy4sQJSU5ONrG4REREVJRZmbJyXFycaLVacXFxMVju4uIid+7cyfQxd+7cyXT9lJQUiYuLk/Lly2d4TGJioiQmJir/f/TokYiIPH78+LllfJaUIqmJz5T1U2xMeonZS/pXJBHy/xsXsdGab9tq3z7L/mJuP4/L/iz5mWjjtf+/+ceSYp2iim2rfftqLrtI3u7nU589k6fa/8pukWLesufl9pMTtRKf9K+ybWtbS7NtO6fb1x23AZivIDDBrVu3ICI4evSowfKFCxfCw8Mj08fUqFEDixYtMlj2+++/Q0QQExOT6WPmzp0LEeEPf/jDH/7whz8vwM+NGzdMiRvZMimmlilTRiwtLTO0rsTGxmZoVdF56aWXMl3fyspKSpcuneljZsyYIRMnTlT+n5qaKvfv35fSpUuLRqMxpchERERUQADIkydPxNXV1WzbNCm42NjYSIMGDeTAgQPSvXt3ZfmBAweka9eumT7G19dX9uzZY7Bs//790rBhQ7G2ts70Mba2tmJra2uwrGTJkqYUlYiIiAoBJycns27P5LuKJk6cKBs3bpRNmzZJZGSkTJgwQaKjoyU4OFhE0lpLBgwYoKwfHBwsUVFRMnHiRImMjJRNmzbJxx9/LJMnTzbfqyAiIqIiweQeTX379pV79+7J/PnzJSYmRry9vWXv3r3i7u4uIiIxMTEGY7pUqVJF9u7dKxMmTJA1a9aIq6urvPfee9KzZ0/zvQoiIiIqEjSAObv6EhEREeUdzlVEREREqsHgQkRERKrB4EJERESqweBCREREqsHgolIpZh6GmoiISA2KZHC5fPmy/PbbbwVdjBy7fPmyLF26VK5fv17QRSEV0GrNPDdSJiIjI+Wbb75hoH5BqfXm09TU1IIugmqZ4zPXbQOAWetQkQsujx8/lg4dOsiVK1fy/LkePHggDx48MOs2z549K02bNpWHDx9KQkKCiKjny6nWnZ9a6abasLS0zNPwcvr0afHy8pJLly6JlZUZJzXNhNrqkFq+m/qioqJkzZo1smjRIvnss89ERFQ51cqlS5dk5syZ0qVLF1myZImEh4cXdJFMqg8FUdevXbsmGzZsEABm+cx120hJSTFrHcrbvUwh5OjoKCJpcyjlpX/++Ue6d+8uAQEBMnbsWClVqlSutxkbGyt9+/aVN998U5YvX64sf/TokTg7O+d4u1euXJGdO3fKxYsXxc/PT9q1ayeVK1fOdXlFRG7evCl3796VBg0a5OnO7/79+xITEyOWlpbi4uKSq/cjM7dv35aTJ0+KVqsVDw8P8fT0NOv200tNTRULi5yfV1y4cEFefvlladWqlezfv18JL5aW5p0t9syZM+Lr6ytTp06VqVOnmnXbImk70t27d8ujR4+kTZs20rp1a7PtVI2R0+e6deuWlCpVSuzt7XP9WT5PUlKS2NjYmGVbZ8+elQ4dOoi3t7fExMTIo0ePJDIyUhYuXGiW7aeXV+/NuXPnpHXr1tKhQwfRaDTy3XffydmzZ+W9997Lco68vHDu3DkJDQ2V119/XTw9PQ2OO9nVrevXr8uePXvkzp07MnLkSKlYsWKel/XMmTPi7+8vlSpVkqFDh4pGo8nVdy0qKkq++OIL2b9/v9y8eVP8/f2lR48e0rFjRxHJ+XdL9+AiIzk5GY8fP4arqysOHz6c5883fPhw1K1bF0uWLEFcXFyut3fkyBE0bNgQ8fHx0Gq1GDZsGFq3bg0PDw8sWrQIDx8+NHmbZ8+eRdmyZdGnTx/4+vqiUaNGCA4OxtOnT3Nd3gsXLsDFxQWNGjXK0/f77NmzqFevHmrXro0yZcpg4MCBuHDhgtm2f/r0aVSvXh0+Pj5wd3dH8+bNERkZabbt64uKisLVq1cNlmm1WpO2cfv2bTRv3hzNmjVDjRo10KlTJ+VvKSkpZiknAFy8eBEajQZTp04FAKSmpmLjxo0YP3483nnnHfzyyy+52n54eDgqVKgAf39/1K9fHw4ODrneprHu3bun/J6ammrSYyMjI6HRaFC3bl08e/YMgOmfobEiIiIwaNAgnDx5Mtfbun79OqpWrYpp06YhNTUVd+/exapVq+Dn54eYmBgzlPY/d+/exaNHjwCY/725efMmXn75ZUybNk1ZtmvXLpQvXx7nzp0z63NlJzk5GS1atIBGo8GECRPg6emJHTt24NKlSwbrpa9fZ8+ehZeXF4KDgw1eg+67a2p9NEZ4eDgcHBzQt29flC5dGuvXr8/V9s6cOYNatWqhd+/eGDRoEMaNG4fq1aujRIkSWLt2ba7LWySCS1RUFOLj4wEA0dHRcHJywp9//plnz5eYmKj8PmHCBHh5eZklvOzZsweNGjUCALRt2xZBQUFYtmwZFi1aBGtrawwdOtSk7UVHR6N27dqYPn26smzNmjWoWrUqbt26lauyxsTEoFWrVmjevDmCgoIQGBiIQ4cO5Wqbmbl48SJcXFwwefJknD9/Hhs2bEDTpk2xYcMGALn/kkdGRqJcuXKYPn064uLi8MMPP6BmzZr47bfflHXMtSO5efMmNBoNSpYsiXfeeQdffPGFwd+N3cHv3LkTPXr0wMGDB/HNN9+gevXqBuElOTnZLOX98ssvodFosHbtWiQlJcHf318Jv82bN0eJEiWwdevWHG07PDwc9vb2mDlzJhISEnD9+nW0bNkSH374oVnKnp3z58/DysoK48aNU5YZ+xnfvXsXrVq1QocOHVCrVi00atQoz8LL1atX4e7uDo1GgzZt2uTqoKzVarFs2TIEBQUZnAAdP34cJUqUwNmzZ81RZABpYatatWro378/7t+/rzy/uezcuRPdunXDxYsXDbbr4+ODL7/8EkDeHPwz8/XXX6NJkyY4ePAgQkNDUadOHQQEBGDKlCm4efNmhjASERGBkiVLYsqUKQbheefOnfjkk0/w77//mr2Muu/ajBkzAACvvfYaOnXqhCdPnuTofQoLC0Px4sUxdepUg9fwxx9/4NVXX4WFhQU2b96cqzK/8MElMTERTZs2RbVq1ZCQkIDExESULVsWYWFhZn8uXTBJf1Y7duxY1K5dO0fhJS4uTtnxnTx5EtbW1li1ahV69uyJ6OhoZb3ffvsNlpaW2LVrl1Hb1Wq12LJlC7p164br168rX/DExERUq1YNBw4cMKmc6R0/fhxt27bFkSNHsG/fvjwJL0+fPsVrr72GwYMHGywfMmQImjZtmuud4ZMnTxAUFITg4GCD5R06dEBoaCi2bt2Ko0ePAjDPjvDp06fo3r07goODMWXKFFSrVg2dO3fG9u3bM7SmZffaEhISsGfPHmW9r7/+OkN4yU3LS1xcnLID3bRpEywsLODi4oJevXrh2rVrANJC2OTJk+Hq6orz58+btP3Lly/DyckJI0eONFjerl07dO/eHa1bt0ZISAhOnTqV49eQlVu3bqFx48aoX78+ihUrhvHjxyt/M+YzPnDgAN544w388ssvCA8Ph6enp0F4MVeLV0JCAmbMmIG+ffvir7/+QsWKFdG8eXOD8GJKnUxNTcWePXvw3nvvKcu0Wi2ePn0Kd3f3TE/0clLnb9y4gSZNmqBWrVrw8/PDiBEjzBZe4uLikJCQgD/++MMg4KampiIpKQm1a9fOdUuCqS5fvozu3bsr+9MrV67gt99+g0ajgb+/P1599VVcvXoV9+/fx+PHjzPd37zzzjvQaDSoWbMmtm7dqnz3zLHPuXTpEjQaDWbOnKks2759OywtLZXP3JTnuXTpEqysrLBo0SIAGev7mTNnEBQUBC8vL1y8eDHH5X7hgwsAnDhxAt7e3mjUqBEuXryIevXq4dtvv8XDhw9x/fp1XLlyBXfu3EFMTAwOHTqEf/75x+Tn0DWbN2/eHIMHD8avv/5qcLli1qxZ8PDwwJIlSxAbGwvg+RXi1KlT8PHxwR9//KEsGz16NKpUqQJ3d3flC687g27RogXmzp1rdJn37duHNWvWKP9PTU3F06dPUbFiRaMDUHb0m6+///57Jbzot1bodlY52WnduXMHwcHB2L17N4D/3ocdO3agWbNmud4RJiQk4MsvvzTYaS9cuBAajQbNmjWDn58fNBqNchaXG6mpqYiPj0dwcDBmz54NIK3pfvjw4ejWrRu8vLzw5Zdf4vTp0zl6Hd988w2qVatmEF7Wrl1r8iWvzOrk5s2b0ahRI4NlQNoZlpOTE3788UeTnmP9+vUoXbo0Zs2apbT8LVq0CLa2tggODsbIkSNhbW2NLl264PHjxyZtOzupqan49NNP0bt3bxw5cgSff/45bG1tMWHCBIN1spOUlIRff/1V+X9YWJgSXnQHHK1Wi9TU1BzXT93jvv32W2zfvh1A2qWtihUrws/PD+fOncu0nM8ru25/kn5dDw8Pg+/svn37clRuANi2bRtatWqFY8eOYdmyZfD19TUILzkNdqdOnUK9evVw/Phxg+X6r8Pf3x8ff/yx8v8PP/wQR44cydHzmWL06NGoVauW8v9BgwahQoUKWL58OTp06ACNRoPg4GBERUWhRo0a2LVrl1Lu77//HhqNBj/99BNef/11eHt7Y8uWLWa5lJ+amoo7d+5g3bp1yv91AgIC0L17d5NaeJKTkxESEgJ7e3ts2bJFWZ6+nu/cuROWlpY4ePBgjsv+wgaXixcvKmedQFpzmIeHB0qVKgULCwuUKlUK5cqVQ4kSJWBra4uyZcvC2dkZpUqVytFlkm+++QYajQaOjo7o1q0bypQpg+rVq6NHjx7YvHkz7t+/j+DgYLRo0QLLli17bjgKDw+HnZ0dpkyZYrD8xx9/VK6b/u9//zP4W2BgIFatWpXtdh88eIDLly8jKirKYLl+pW3YsCG+/fZb5f9btmwxKh0/ePAAV65cwaVLlzINJHv37kWHDh3Qvn17peVl3LhxGQ54z5OYmKiElPDwcGW57rn279+Pxo0bIzk5WXldplyjT0xMRFJSUoZLKt9//z1eeuklfPvtt4iPj8ezZ88wcuRI+Pj44N69ezk+A0pNTVUee+rUKTg7O+O7775T/t6mTRvY2tqiadOmqF27NkaOHGnweVy/fh2fffYZVqxYoRzI0x8A4uPjlctGr7zyCsaOHQuNRpOhP012sqqTQFp/Jt3lWN3nEBERAW9vb5w4ccKo7V+7dk1pRXznnXfg4+ODt99+G9OnT0fZsmXxww8/KOvu3r0bGo0mVzu/zERHRxvU/c8++wy2trZGtbxktTyzlpd3333X5EAHAOfOncPKlSsBIMNBJTY2NkPLS2pqqtIqmJ6u3ixfvlypN/qvQavV4smTJ6hQoQJ+//13AMDs2bOh0Whw48YNk+r7gwcPlN/191tLly5VwovusoL+PsOYIJNVvUxfvsDAQOUSxYwZM+Dg4GD2vmr6Zdc9/z///IOAgADs378fb7zxBl566SWDlrGNGzfiypUr+P3336HRaBAREaH8LSUlxeBk5dVXX0X58uWz/EyNdfbsWaVvWmaWLl2KqlWr4u+//87wurITGRmJSZMmwcPDw6Avi1arVbaRnJwMR0dHgxBpqhc2uMyaNQsajQZfffUVgLRKdOrUKbRt2xY2Njb44YcfcOLECRw6dAhHjhzBsWPH8McffxhcfjHVF198AUtLS4SGhuLMmTP4/vvv0blzZ9SpUwcVK1ZEp06doNFo4OTkhNWrV2dZGcLDw1GsWDHlmmN6//vf/9CsWTNYWlri/fffx+7duzFt2jS89NJLuHLlSpblO3fuHPz8/ODm5oaKFStixIgRysFZ/0veuHFjfPPNNwCAmTNnokSJEhk6lGW37UqVKhlsWz8A6C4bdejQAd26dYNGozGpyf/y5csYPXo0+vfvb9B/Qr/8O3fuRJUqVZQD6bx589CmTRs8e/bsuTtb3fYHDhyITz75xOBvWq02w2WPWbNmwc/Pz+jy67t+/TqOHTumlF/3Po0fPx5z5swBALz55ptwdXXF1atXcebMGaxduxaVK1fG9evXAaR1HK5cuTIaNmyI0qVLo2rVqkqnRx3da05MTFQO+KVKlTKpQ+fz6mRmpk6dinr16iktjNmJi4tDuXLlsGPHDmXZggULULNmTdjb2+Ozzz4DkLYjT01NxbFjx1C9enX89ddfRpcnJ1JSUjK0vCQnJ+PTTz9V+n0Ys1PXhZemTZti+PDh0Gg0z/1OpafVatG5c2eMGTMmw990defu3btKy0tYWBhGjRqFRo0aGfQ1AIyrN8nJyXj48CFcXV1x+vRpLF26FMWKFTM6iOpcvHgRvr6+GU6WdK8ps5aXlStXGv2+GlMvtVot6tevj02bNmHx4sWws7Mz+XVk5/bt20orSGYnDa+99hocHBxQs2ZNpZtCamoqrl+/jsuXLwNIu4xiZ2eH999/P8M2dJ/vsWPH0KBBg1z3OXrzzTcxcODADMt1+4r4+Hi4u7tnuFxrjIiICIwfPx41a9ZUWnOA/74nP/30E2rXrm3yJWR9L2xwSUxMxOTJk2Ftba1cStBqtTh16hS8vLzQsmVL5cBmTh999BE0Gg0WLFigLEtJScGnn36KVatWoWbNmqhcuXKWLRgxMTGwtrbGsGHDAPxXkWbNmoW33npLWS88PBzTp09HmTJl8PLLL6Nx48bZ9tsJDw9H8eLFMWrUKOzYsQNDhgxBxYoVMXXqVIODeXx8PKpWrYqvv/4aS5YsMeoLbsy29Z9jz549cHZ2RsmSJQ1aTJ7nzJkzcHV1xdixY7Ft2zaDTtD69uzZg9q1awMA5syZA2tra6MO0tltP6szv9GjR2PYsGFITEw06Qw0Pj4egwYNgpubm3LHle7xW7ZsQeXKldG+fXtUrFgxw8FZVy5dp7qQkBDExsbi77//RuXKlfH5559neD7dtocNG4ZixYqZtNPIrk7qt0TonDlzBmPGjIGzs7PRn++tW7fg4uKCc+fOGRywVqxYAS8vL8yYMQM3btxQls+cORN16tTB3bt3jX4dmUlKSnruOsnJyQbhZezYsbC2tjY4azbmIHv8+HGltTen/XNatGhh0B9Bn+61xMbGonLlyihWrBjs7Owy1H1T602DBg3g6+sLW1vbDJdijLF582Z4e3sDMHyfdL+npKQYhJehQ4caFexMrZft2rVDxYoV4eDgkKPXkZXIyEiULVsWAQEBSvjLrNOtq6sr3n33XeVxurCuOz4BaZdoqlevnmVL0OTJk9G+fXuDy3o50aVLFwwfPjzTv+k+l4ULF8LHx0dpdcnMkydPEBsbi/379yMsLEx5vZcuXcK4ceMyhBcg7YaVTp065eo1vHDBRf+LkZiYiAkTJsDa2tqgN3l4eDhq1aqFJk2a5Lij0/Xr1/HBBx/gnXfewbZt2wz+tnHjRlhYWGDevHkZLjfcu3cv2w/s2rVraN26NSpXrqxUmMWLF8PZ2dng0pfO3bt38fjx42xvhb5w4QKKFSuGkJAQZVlSUhICAgLQtGlTg4OyVquFn58fvLy8jPqCm7Jt3XX98ePHm3ynwrVr11CpUqXnNgcDwMGDB9G2bVtMmjQJtra2Rp1ZmbJ9IK0j7axZs1C2bFmDpl1T/Pzzz3j11VdRt27dDJ2WO3XqhFKlSmX6/qempuLy5cvKwUefv78/pkyZgn79+mHz5s0GO519+/ahVq1aJp9pmlInL126hOnTp6NRo0YmhdKLFy+iQoUKymVa/Xqj24FOmTIFDx48wKJFi2Bvb2/S9jNz7tw5dOnSxagQl5KSgm3btkGj0cDZ2Rn79u1Dhw4dMnRmzYpWq8Xw4cNhb2+fqzPNtm3bGvRLS0+3vxk8eDBKly6d4U4jU+vNgwcPULJkSVhbW+eofxUALF++HPXq1cv0b7r3LDU1FUuXLoW9vT1Klixp1M0Tpu4rdZfwz5w5k6PXkZlbt26hefPmaNiwIZo1a4bOnTtnCC9arRbPnj3DkCFDMGTIECQlJSE1NTXTsP7XX3+hatWqqFOnDv766y8ljN65cwdTpkyBo6Njjsuvvy/r2rXrc/tDnjt3DhqNBps2bcr07xcvXkSfPn1Qu3ZtaDQaFC9eHI0bN1ZC15UrVzBu3Dh4eHgo4WXu3LkoVapUrm9Lf2GCi34zp/5OTxdebGxslMqckpKC8PBwuLi4oE2bNiY/l+7MPDAwEHXq1EGlSpUy7Ah04WXBggUm338fFRWFgIAAVKpUCZMmTUK5cuUMru/rGLM9rVaLwYMHo1SpUti5c6fB4xYsWIAmTZrg4cOHyrKUlBQ0b94czs7Oz91RmbptIO29q1ChgskHz5UrV6Jdu3aIiYl57uvW3aZbrFgxoy+HmLL9n376CcOGDUP58uVzfWfLwYMH0bt3b4PwotVqsW7dOjRp0kS5JJT+THXOnDkoV64cli1bpixfvHgxLC0t8frrr8Pf3x/W1tYYNWoUnjx5AiDtu3Dnzp0cldPYOgmk7dCMuTwUGxurvL4LFy7A0dFRuStJ9zp1Fi5ciCZNmqBu3bpmaea/du0aqlWrBo1Gg3r16j23D5dWq8WQIUPg6OiIiIgIREZGonv37vD398dHH32UaZn1nTx5Ek2bNjW53FFRUfjwww+Vz1B3uQMw/P7r/667CyV93TS13jx69AgpKSlYt26dyXeA6J+whYaGokWLFgAyb7nUvWejRo2Co6OjSQc1U+rlN998Y1KfLmPs2LED7dq1w6FDh7Bt2zb4+voahBf992HHjh3QaDRKKEsf1oG0k769e/fCy8sLJUqUQIsWLdC+fXv4+/vD3d09x/ub69evY9myZbh9+zaAtH5zutYf/VZx3e+6z2TWrFmZBu3Tp0/jpZdewpgxY7Br1y5ERUXhvffeQ+3ateHm5qb0iYqMjMT48eNRp04d+Pv7w97e3ixjDr0QwSUmJibD2Y/+FyQpKQljxoyBnZ2dwcH47Nmz2fYJyUx2gzTdvn3bYAeyceNG2NraYvr06UY1JeuXOSoqCl27doVGo1E6lJl6F4Lu4HHz5k307NkTLVq0wMaNGwGktdQUK1YMS5cuzfC4Tz75RLnuau5tA8hwLd0YXbp0QVBQUKZ/070vT58+RWJiIk6fPo0WLVqY1PHO2O2npKTg+PHjWLJkick7c91BaOPGjUrfFgA4fPiwEl50d288ePAAlSpVwujRozPdVmxsLCZPnowmTZogNDQUS5YsQdmyZbFv3z6lDk6ZMgX29vYm96XQZ0qdNKXV8tGjR/D29saIESNw69Yt/P3333B2dlYuB2VW199++23UqFEjx2f+OgkJCZg3bx569OiB48ePo3HjxvD09Mz289y7dy+qVKmCffv2Ka2b58+fx8CBA+Hr62sQXvTfs+TkZOVApd9B1VhjxoyBh4cHPvjgA6SkpKBevXpZngHrHD16NMu6n5N6Y+q4Pzdu3ED9+vWV9/Odd95Bly5dst3W7t27YWNjY/RBLa/qpbF0QRJI67Svox9edJ+3/uXIjh07KoMopg/r+uX8559/MG3aNPTt2xddu3bFqlWrsr1k8zxz585F1apVMX/+fDx79gytWrVSOnlnJ7PPS9evaPr06QZhJykpCadPn0bjxo1RpUoV5Xty8eJFDB06NFfBK70XIrhkdfajX7ljY2PRrl07jBkzJscDcKWkpGDp0qUmDdL0wQcfoHTp0lneRRQTE5PljvjKlSvo1KkTKlSooJwpGBteoqOjUb58eeVW3ps3b6Jbt25o27Ytli9fjooVKxp08tOvgHmx7cz+NUZqaioSExPRrVs39O7dGwCy7Nsye/Zs5TZuY6+fmrp93SVHY/pG6Dt9+jQqVKgAPz8/lC5dGkFBQQZ9V3777Tf06tULdevWVXZsCxYsQJ06dfDo0SOD90xXr//55x9MmDABtWrVgqWlJfbv3w8Ayp0r+/fvR9WqVU2+NJFXdTK9Dz74AO7u7pgwYQI++eQTeHl54c6dO3j48CHu3buHBw8e4MGDB/jnn3+UnXZORofWp7sVeffu3UpdefDgwXPDy61bt3Dx4kV06NABr7zyilK/9MOLbtBD3fMkJiZi9OjRaNy4cY5CC5B2gBwyZAgaNWqE5cuXK+NBffnll9i1axe++uor7Ny5E1988QXWrl2Ln376KcttmVpvctqcf+rUKTRu3Bg1atTA9evX8c477+CVV17Jcn1d/Xne3Zz5VS+fJzw8HEFBQZmeEOguKerCi66+vvvuu4iMjISXl1eWYT2rfnTmCl7Tpk1DvXr1MH/+fHh5eWHs2LHYsGED1q9fjw0bNmDt2rVYu3YtFi5cqHSGTy99v6LM7h797bff4OzsbDCA49WrV3Pc2psZVQeXmzdvZroDySq89O7dG6+++mqOnkt3MPvuu+/w/vvvK8ufN0gTkPXO9vbt23ByckKpUqUwceJEHD58OEMlvXbtGgICApQ7S3TP+TwHDx5EpUqVkJiYqKyvCxhOTk5o1aqV8lymBrm83HZWQkJCYG1trews0n/Jb9++jV69eik74bzafnYHhqzoOkNOnz4diYmJ+OOPP1CuXDnlzi2dgwcPolevXmjQoAF+//133L17V7kTI32Y0pUvLi4OkyZNwssvv6wM+qQzYcKETO8oyU5e1snMfPTRR6hUqRJat26t9B8pW7YsXF1dUapUKZQsWRIuLi5wd3fPdUfcixcvKkOop38/4+LilPCi39Jw5MgRpRN/SkoKVqxYAX9/f/Tr1++54WX06NEmtSLoxMTEYP/+/cqZ/L///osBAwagbt260Gg0sLKygqenJ8qUKQNnZ2e4urqiUqVKcHV1zdDSkl/1Jr0zZ86gdevWqFatGnr27ImmTZuiR48e6N+/P4YPH45+/fqhT58+6N27N0aNGvXcg3N+18ushIeHw8rKCrNmzcrwN/3bfT/99FP4+vqiS5cuGDZsGDQaDS5fvmx0WI+Li8P58+fNfgPJhAkT8PLLL0Oj0cDGxgZNmjSBm5sbXF1dUbNmTdSqVSvbkx39fkVZvcepqanKODB5RbXBRTfKYLNmzZTRaM+dO4c333wzww5Ed3Y8ePBghISEKGddxjpz5gxGjBiBhw8fGgx4ld0gTXv37n3udiMjI9GuXTts3boVvXr1QmBgIJo1a4a//vrL4Lbs6OhotGvXDra2tgZ9ALKzZ88eVK5c2aCDGJC2U+zRowdatGhhMEiQKe9HXm4bSHu927Ztw9q1a5Wd/qVLl1CjRg3UrFkz09sq586di0aNGinXcAty+/oiIyNhZWWlDCqn0759e4wYMQLBwcFYuHChsvy3335DYGAg/Pz8lINOZGQk6tevjyVLlmQazHRn0I0bN8bbb78NIK1PSLFixUy+rJKXdTImJgaHDx/OMCrz5s2bUbJkSfj6+iqv8eDBg/juu++wZ88eHDp0yORLupnZtWsXNBqNwQBq6ZvndeHl3LlzCA4ORuPGjREXF2dwB8yaNWvg6+ubZXhZv349Jk6cCHt7e5Obxs+dO4f69evjlVdeQffu3ZVWkH///RcjRoxA3bp1MXnyZOVkKDY2FgkJCdBqtRkuweZnvYmKisLGjRuxYsUKpa/W+fPn0aVLF2g0Gnh5eWHMmDHo0KEDunbtih49eqBXr17o27evUc+Vl/XSWKdPn860Y3P6sW+AtHqyZcsW2Nvbw9nZ2aAeGBvWq1WrluNWikuXLmHOnDno06cPNm/ebHBH3uzZs/Hyyy/jrbfeUsKp/nHteYPORUdHPzcgduvWDV27ds1R2Y2h2uCi1WqxadMmtGjRAkFBQQbhZeDAgWjatKnSAe3u3buYN28eypUrZ/Lke+Hh4bCwsFC+2Dr6nZhyM0hT586dMWLECABpX84hQ4agVatW8PX1xfbt25WzzJiYGHTt2jXb/gqPHz9WEvqXX36J8uXL4+nTpxma8/Qv7Rg7BHZeblvf6dOn4e7uDl9fX2g0GtSuXVtptty4cSMqV64MNzc37N69G2fPnsUPP/yAESNGwMnJyai7TPJ6++mtWbMGGo3G4FbTRYsWQaPRYODAgcrImf369VP+fuzYMYO6M2PGDGg0GqxatQolS5bEuHHjMsxjdOfOHYwfPx5+fn5o0KCB0XdTZcacdVLnzJkzqFmzJurUqQONRpNhXq0tW7agQoUKmDZtmtJh15xSU1Px77//ol+/fmjevHmWnTTj4uLg6+sLCwsL2NvbG1zO0wV1/fDy2muvKZeBIiIiMGTIEDg7O+eopeX8+fNwdnZGSEiIQWuH7nn//fdfvPnmm2jSpAlWr16NhIQEAIZ35uj/m1/15vTp03Bzc0ODBg1gZWUFd3d3pX/biRMn0KdPH7i5ueX4cplOXtRLY125cgUajcZg0kMAWLVqVYY7mHTv/8iRI2Fvb4+tW7fma1g/c+YMXFxc0L17dwQFBcHJycngRBIAJk6ciPr162PevHnKsTN9PUovfb+i9OFF9/dbt24hMDBQGYAuL/oYqTK46N6gxMREfP7552jVqhU6d+6snP1ERERg0qRJKF++PFxdXeHn5wdPT0+T5ycKCwtTJnrTp998l5KSkqNBmnSV5OrVq2jcuLFBxW7cuDFcXFxQsmRJdOjQAYMGDcKzZ8+yHUXy9u3baNasmTJHh25+Gv1LNikpKcrzRkVFoU2bNhn66+T3tvWdOXMG9vb2mD17Nh49eoQrV66gXLly6Nixo7LOF198gTZt2ihNnbVq1ULLli2NukUwr7evLzo6WtkhvP3227CyssI333yDpUuXolSpUsplgISEBCxZsgS2trYZhh/XtbhERUWhVq1a2LdvH06fPo033ngDfn5+8PX1xXfffae0EMXFxWHkyJGoXr16jkKWueukjv736PLly9i/fz80Gk2G8WnWr18PV1dXTJ482Wyze6fvi7Rr1y5Uq1ZN2almdilh0KBBKF26NP73v/9h1apVBgdc/fCydu1aNGvWDOPGjVPOUiMjIzFq1CiTR2R99OgR2rdvn2GeGv27/YC08DJo0CA0b94cixcvzrQ/Vn7WG913as6cOYiNjUV0dDQ6deqE6tWrK8ME/Pnnn2jbti0qV66shAnd52LsnZGA+eulKXSj2o4dO1bZ/y9duhQWFhaZzli+b98+5RJnfob1K1euwM3NDSEhIcp70L9/f0yYMCFDXZ80aRIaN26sDDOQmez6FV2/fl0JL/oha8aMGfD29s605dpcVBVc9PtL6H9h+/btC1tbW3Tq1Ek5WMTGxiIiIgLvvvsu/ve//5n8Jp4/fx729vaYP3++wfI1a9Zg06ZNGfpu5GSQJq1Wi/v37+OVV17B4sWLAQADBgyAi4sLbt68iYiICCxYsACenp5G9Sjv1KkTvL298dlnn2HZsmVo27Zttus/ePDAoAmxoLYN/DfZV/pRQVu1agVXV1fcvHnTYPmRI0fw448/4uLFi0adyeX19vXdunULFStWNBiOXtcKZ2lpqex4dTvtr7/+Gu7u7gZniGfOnMGQIUNw9+5dJCcnY+LEiUoTta7uW1lZKQMabtiwAREREUhJSTHqVuSsmLtOXrhwAdbW1liyZInymnW33G/atAkrVqwwGPr+448/hp2dHUJCQnLdPyoyMhKNGjXCjz/+aHAZZdCgQXjppZeUz1X/TPO9995TbiXeuHEjNBoNlixZYhDAdQeEhIQEzJkzB3Xq1DG4C8/UjttAWktllSpV8PXXX2f6d/2DztOnT9GzZ08EBARk6ISen/Xm2rVrKFasmEFrIQAcOnQIdnZ2+Pnnn5VlJ06cQKdOnTLc8m4sc9fLnPjxxx9hb2+PSZMmYeHChShVqlSWk9EWRFhPSkrCokWLMHr0aIMO/YMHD0aXLl3g7++P2bNnG3zfRo4ciVatWmV684ip/Yru3r2rnLjndoyl51FNcLly5QpmzJiR4UxmxYoVcHZ2xrRp0+Dn52dw2SinHjx4AG9vb3h5eRl8mZcsWQJra+sMg4U9fPjwuYM03b59G/v27cPevXszTCvw3XffwcXFBa1atcJLL71k0MScnJycbQeta9euYe3atUpzXe/evdGkSRMEBQVBo9GgadOm8Pb2RpMmTdCgQQPUqVMHXl5eaNOmzXOvZeblttM7ffo0NBoNBg8erASeJUuWQKPRoEqVKujatSv69OmDhQsX4sGDBybfXZLX29cXFRUFFxcXg1m3gbQ7CzQaTYYBC6dNm4YmTZoo9TY8PByWlpYGlye/+uorlChRQgngQ4YMgaurK/bs2YN33nkHjo6O8Pf3N+lW87yqkzpJSUkYPnw4bGxsDCYe1F0q0w2y5+bmZtBB9NNPPzVLM78uLFavXh3z5s1TLpMkJiaidu3a6NGjR4axUH799Vf88ccfOH78OO7du4cvvvgCGo0GixYtyjS8JCYmwt7ePlfzrgBpd/JYWFgoB5DMWgx0Y3wAaS0v6ftb5Ve90blw4QJcXFzQp08fHDlyRKnrv/zyC5ycnAxu+QfSJt3s2bPnc4dayOt6aaxr167hgw8+MJjn7IcffoCTk1Omc8XpREREFFhYDw8PNwhHc+fOhZ2dHWbMmIHJkyejZcuWeOWVVwze16w6vZvSryggIEBppTbnVApZUU1w0d3nP27cOKXiL1myBM7OzsqZ7WeffYZmzZqhS5cuuQ4vCxYsQOPGjTFp0iQkJSVh1apVKFWqVJZ3raxZsybLWynPnDmDGjVqKL2527dvr+zIddfeO3XqhKpVq5o8BH7NmjXRvXt3g0n5XnvtNWg0Gvj6+mL69OlYsGABli9fjmXLlmHevHlYsWLFc291zMtt64uPj1fOUI8ePQpbW1sMHz4cM2fORKlSpbBnzx6cPn0ax48fx5gxY+Dj4wMHBwf06NED8fHxz21qzuvtp6cb56V06dKZ1sE5c+bA0tJSGYtjzpw5cHBwUD738+fPw87OLtNRLQcMGKB0uHNxcTGoK6dOncrQapSdvKqTOvfu3YNWq8Uff/yBgQMHwsvLC8ePH8cHH3ygjGyampqKJ0+eoGXLlmjYsGGu7xpKLzExEW+88QYCAgIwY8YMtGrVCsOHD0d8fDzef/99+Pr6KgMn6pw/fx7Nmzc3uCsiNDRUCS/6B/ikpCTExMTA39/fIJjlxN9//40SJUoYjH2Uvu5t3boVTZs2zTRU51e9AdL6u+nGMQkLC4OHhwc6d+6MCxcu4MaNG3jppZcwadKkTF+Hrl9OVvK6XhrrzJkz8PDwQEBAgDLfnc6vv/6qTHGi6zytUxBh/dmzZwYdyHXu37+PPn36GPTB+e6772BnZ2f0HF/G9iu6desWhg0bZtZRibNT6INLbGwsjh8/jjt37uDHH39EhQoVEBISgsmTJ6N06dIGQUKr1WL79u2oXbs2+vTpk+tb4ZYsWYIGDRrAz88Pjo6OSudb/S/iJ598ku04GadPn4aDgwNmzJiB27dv4+uvv4aNjQ0mTpxosN6yZctQvnx5pXf388oeGRkJZ2dnTJ8+PdPxD/r37w9PT0988cUXz91Z5Oe29UVFRaFhw4b4+eeflWbsI0eOwMHBIUOHVuC/MTg2bdpk1AiYeb19ffqdKY8ePYoyZcpkOXbP3LlzYW9vj9atW6N48eLKGcrZs2dRpkwZeHp6Kuvqn4GtWbMGdnZ2qFWrVq6Gjc+rOqlz9epVeHh4KE3ff/75J15//XWUL18eNjY2yqVU3WeyevVqVK9e3ezBJTk5GZs3b8bIkSPx22+/4cSJE6hVqxZ69+6NZcuWoUaNGhg0aJDyus6dO4eSJUti5syZiIqKMrjks3r1amUOMv3vxNy5c1GrVi2TLosCaS0mKSkpynsbFxeHhg0bol69esrcVenNmDEDI0aMyHApKr/qDZDWGtK8eXNs2LBBKfupU6eUg7yLi4vBxHym7IPzul4aKzIyEqVLl8bUqVOzrJN79+6FnZ0dhg8froSXggjrERER6NixI2rXrg0fHx9s375dmewRgMGt/EBa0KxXr95zL0uZ2q/o33//NdvQF8Yo1MFFd/bTrl075exn8+bNKFeuHKytrQ1m79XvNLdr1y6Tr6PeuXMHv//+O/766y+DFP3uu+/Czc0NPXv2VG5N05/Mq3jx4lm2tERERMDKysrg7ANIu3W6fv36BhUMAOrWrauk2+w8e/YMvXr1yjCqalJSEq5du6Y0IQ8bNgw1atTAtm3bDEZ6LKhtZ6Z69eqoVasWDh8+rBzITp48CTs7OwwYMMCgb1JOOtzl9faBtEuFZcuWxfTp05Xtu7u74/79+1luc/bs2bCzs1NukwwPD4eDg4PS52bs2LHKuvo7hObNm+P111/PUTmBvKuT+n799Vd4enoa9EP7888/8cYbb6BKlSpK3wfd92jEiBEICgrKVT0C0i4nd+3aFWfPnlVeR3R0NFq1aoWpU6cCSPuMp02bhrfeegvlypWDRqPBli1bcO/ePfj5+RlMZAoYvvcffPABnJyc0LNnTwwePBhDhgzJ0IJhjIiICHTu3Bl+fn7o0KGD8vhTp06hRIkS8PPzMxhO4d69e5gyZQpcXV0zzIuVX/VGX7t27VC3bl1s2bJFCQ9hYWHw8vIymDQUMP6Okvyol8aIj49Ht27dMnSSTklJQUxMDGJiYpQWr71798LR0RH9+vVDREREvof18PBwlCxZEgMHDsT8+fPh6+ubod9N+vd/2rRpaNasmVFj9BSGfkVZKbTBJf3Zj/6X8Msvv8RLL72EiRMnGjSv5fTAo2sW9Pb2hkajwdSpUw3OapYvXw4fHx9MmDBBOdDNnj0b9vb2WXbE1Wq1+P7776HRaLBy5UqliVk3j0j9+vXRr18/LFiwQLkd980330S3bt0M7qnPTFJSEvz8/AwGwvvhhx8wfvx4ODo6omLFiujRoweAtIDh4uKCHTt2GPVe5OW2dXSj1eo0btwY1apVMwgXR44cga2tLV5//XWTz2bzevv6dGcm7777Luzs7LBw4UIcPXrU4Cxcd5ala536+++/kZycrHQOPX78OKytrTFv3jykpKTgww8/RJkyZQwOQrrHfvjhh/Dz8zP5rhVdWfOqTurbvn07qlSpkqE1Trcz9/LyUsZSmT17dq4mjtN/bYcOHUL58uVRtWpVTJ48WbmkfOLECZQoUQKffvopgLTPIywsDNOnT1dmeT5//jyqVauGgwcPZjiD1x/36dtvv8WYMWMQGBiIKVOmmDzBZnh4OBwdHTF69GiEhISgY8eO6Nixo3Jg/vnnn1G6dGmULVsWbdu2RY8ePdCuXTtUqFAhw5gw+VVvMtOjRw94eXkZhBfdfrRr164Z+rdkJ7/qpTGePHmCxo0bG7TG/vjjj3jrrbfg6OgINzc39O3bV2mR/e6771C+fHns3r07X8P6mTNnUKJEiQxjylSsWBG9evXKsP7NmzcREhICJyenTL9rhaVfkbEKZXAx5uxn27ZtqFChAsaOHfvczl7ZCQsLg4ODg3JLmm4G2PRf8EWLFsHHxwczZ87E6NGjs53oTdezH4Byl8KaNWswb948ODs749NPP8Wvv/6KDRs2oFevXihbtix8fHzQpUsXo+a/efToEWrVqoVhw4YhMjISixYtgoeHB3r27InVq1fj448/hru7u9JJb+DAgUZf+sjLbQNpo5eOGTMG3bt3N7jG26hRo0zDRfHixdG5c2ejr8Pn9fb1nT9/HsuWLUNCQgISExOxbt06WFhYoEmTJihVqhRcXFzg6ekJT09PVKhQAdWqVUP16tXh7e1t0Ffit99+MzjYPHz4MNODEJDWqmBpaalMkGasvK6T+geQjz76CNWrV890h6YLL/Xr10e3bt1gb2+f6858kZGRmDBhgvL/qVOnomnTpnB1dcXOnTuRkJCAtWvX4tVXX80QNHRnz5999hmsrKwyTDKn799//zV4vKmXKM6ePQsHBweDfiirV69G165dcf/+feWAoTvIBAUFoXPnzliyZEmm37H8qDdAWifVVatWYeHChcqUF0Dm4eXUqVPw8vJCq1atshxJXF9e10tTJScnw9vbGz169MDDhw+xcOFCeHh4oHv37vjoo4+wYsUK+Pj4YMWKFQZzmOVnWE9OTkb79u2h0WiU/Yjuu9anTx/069fPIECdOnUKr7zyCry8vDIdEqSw9CsyRaEMLtmd/ejPqfPpp5+iUqVKGDRoUI5m/YyIiICNjU2GyQD9/PywZcsWbNiwwaAPzdKlS1G2bFmUKFEiy8GldD3758yZoyzTde7TaDQZhnlPSEjAjRs3MGPGDJMub/3888/KYE8lSpTA+vXrlQCXlJSEwMDADLcpFvS2w8PDUbZsWXTr1g2vvvoqrK2tDcKFbnIu/XDx22+/wcXF5bnzmOTH9tM/l0ajMag7iYmJ2LhxI5ydneHu7o7PP/8c7733HtasWYM1a9YoEyxmd31ZV7cfPXqU5UFo1apVJnWAzus6efv2bXTo0AFr1qwBkNaJ1NfXV/m7bt4enWPHjqFz584oU6aMWWaK3bp1K2rWrGlwJ9tff/2FESNGwMLCAoMHD8b06dMxaNAgpXVQVx7d+33kyBHY2dkZHJjTe++999CuXbsc9eu6c+cOHB0dERgYaNB8P2XKFJQuXRo1atRA8eLFMX369BzdTp0X9Qb4b3A5Pz8/VK1aFQ4ODsrnDKTdaZg+vPz1119o1KjRc1sy82tfaao9e/agXLlyqFChAhwdHbFu3TqDlv1mzZrhtddeK9CwfuXKFXh4eMDPz08ZA+b27duws7Mz+Hx0fvnll0yHBCks/YpMVSiDi7FnPzdv3sTXX38NT09Pk4dGTk5OxqBBg6DRaAwG/1mwYAE0Gg3atm2LsmXLonz58ggNDVX+vmHDhiyTfvqe/frl/vjjj6HRaLB8+fJMB7TKiejoaJw4cSJDJ1CtVovevXtj1qxZJk2emJfb1g2XrRvMLyUlBWPGjMH48eMNWh9atWqFypUr4/fff1cOEMY0Reb19vXp+hWkH5gQSGt5+Pjjj2FlZYUPPvjApO2mp38Q0t+RmPJ55nWdTEpKwp07d9C5c2c0b94cO3bswNq1a9G5c+cst5OUlITIyMgctXJlZv/+/ShRokSm0zHs2rULgYGBytmku7t7prf+3rx5E+XKlUOXLl0M9gf67/WkSZOUGXFz4pVXXoGPj49yO/yKFStQvHhxbNmyBT///LNy50n6TuOmPp856g3w30FNN8dWeHg4vLy84O3tbfDZ9ezZE7Vr18a2bduU9/Z54S6/95VZuX37Nn744Qd89dVXyr4uISEBt2/fxsGDBw2OK1qtFgkJCejduzemTp1aoGEdSLvkXL16dQQEBODPP/9EpUqVMGrUKOXvz9s/F5Z+RTlRKIOLMWc/q1evRrt27QCYPmOsLnlevXoVQUFBynDUa9euhbOzs5L0L126hICAADRt2vS5c9Rk1bNf/wupf1tlbme5zUpiYiJmzZoFV1dXsw55nZttR0dHo0yZMsrsyzp9+/ZF3bp1UatWLbRt21a57bply5ZwdnZWrpM/b4eb19vXpxujQddZTWfnzp1Kp27dpQlLS0uDAQxzEiIfPXqEjz76KNPhxp8nr+vkjRs30KtXL4SHhyMmJgb9+vVDYGAgfHx84OzsjICAALRs2RKtW7eGv7+/cidFly5dspyB21RarRb37t1D5cqVlVs805/oXLx4ER9//DHKly+f4URF3+7du2Fra4v+/fsb3H3z77//YsaMGXB3d8/15YkePXqgbt266NOnD0qVKmUwSGFiYiI8PT3x5ptv5uo5gNzVGyDr71SbNm1QoUIFxMTEGBzYBgwYgHLlymH79u3PreeFZV955swZ1K5dGx4eHihVqhRq166dbZjWarWYPXs2KlasiKNHj+ZrWE8/tYru32vXrsHDw0OZQiT9+tltr7D0K8qJQhlcjD37mTJliskHg4SEBDRu3BjVq1cH8N+cC8WLF4eDg4OShnXbfPvtt1GzZs1sP6zsevbrD4UPpAUua2trZdh5c9q2bRvGjh0LFxcXkyd3y8ttX7t2DY0aNUKXLl2UW8oXL14MBwcHzJ8/Hxs3boSnpycqV66shMq2bdsa3Xcpr7evb9q0adBoNAZnTbqB7PTfl8TERHz44YfKmWNuPHz4EJs3bzbpoJkfdfLq1ato0qQJOnbsiHPnzuHmzZt4/fXXlZbK+fPnY9SoURg9ejSmTp2KsWPH4q233sLZs2eNfo7MXLhwAbt27cKhQ4eUg2eVKlWwatUqg/XS3555//79bC8xaLVarF+/HlZWVvDw8MCgQYMwcuRIdOnSBeXKlctRvc+sb8irr76q3ASgf6B79uwZ/P39M4TinMpJvdEve/rvlK5FqG7duggMDERAQACmT5+OCxcu4OnTpwgODn7u/DqFZV+pK8e0adPw999/Y+vWrbC2tkb//v2VEKV/XPnss88wZswYJZDnZ1g/f/48hg0bhqNHjxq0ZOmHl3r16qFZs2ZGTQJb2PoV5UShDC5A3p39pKam4vDhw6hduzYaNWoEIG3OhZ49e6J06dLKbV2668wjR45E586dMzSb6RjTsz/9F1I3cF5uB8nTd+HCBbRq1Qrdu3c3+U6H/Nj2pUuX0KFDB3Tp0gVDhw5FuXLlDEaPjIqKgkajMbibqTBtXxegk5KS0K9fPxQrVgyXLl3C+++/jzJlymQ6MGFiYiI2bdpkls/DlHCen3Xy0qVLaN++PQIDA3H27FncvXsXr7/+Opo3b55h4C5z0Gq1GD58OKpVq4bixYujevXqqFmzJmrWrIlXX30V33zzDeLi4kweuVnfn3/+iV69esHHxwd+fn6YNm2aya2Xz+sb8uqrr8Lb2xtbtmxRyqo7m8/NzQbp5WaCu/TfqbJly2L37t24e/cuDh06hA0bNsDDwwNly5ZFYGDgc/vmFJZ95dWrV+Hg4GBwiUSr1aJy5cpo3759hvX//PNP9OvXD926dcMPP/yQr2E9JSUFQUFBsLGxQcWKFTFs2DAloOu/T9euXUO1atXQsmXLDHcF6Sus/YpMVWiDi7nPfgDD2VOPHTuGmjVrGoSXdu3awc3NTek8OWvWrOf2/ja2Z3/6L2T6OUbM4e7du3nWrGqObV+8eBHt2rWDvb09VqxYASDts0hKSsLNmzdRt25d7Nq1S1leWLafkJCAJk2aoGrVqkhNTVX6+VhYWMDOzs6ouyfyU37XyfTh5datW+jXrx+aNm2KzZs3K+vltL9VZuLj4xEdHY3Dhw9j8eLFaN26NTQaDapVqwY3NzdUq1YNw4cPx6RJk3L0+eSm86GpfUO+/PJLTJ8+Hba2tmbr/2Au+t+pzFoOnzx5gqNHjxo1k3Fh2VeuXLlSmdxQd3lX12r68ssvY+TIkRg1ahSOHz+u9LG5deuWsv/L77C+bt06vP322zh+/DhWr16NChUqIDAwEAsWLDDog3jt2jU4Ozujffv2mV6yKiz9isyh0AYXHXOc/eifCehPqJZZeAkMDET16tUxdOhQk3t/P69nf0pKSoap54uaK1euIDAwEEFBQQZzPs2ePRtVqlTJ9myhoLava6Xz8vJCgwYNkJqaiuTkZIwcOdJgcKnC+JnmV53U35mfO3cOt2/fxoABA+Dp6Ynt27fn7kX8P/3ypd+hLlu2DNWrV8eDBw9w4MABLF26FEOHDoWXl5dRB9XsnsuU9yUnfUM0Gg2KFy9e6EKLjv53Sn9wudyMlFoQ+0pdq018fDwWLlyIJk2aICQkBG+//TbKlCmD999/H3/99RcWLVqEXr16oUyZMnB1dcXs2bMzbCs/w3pkZCRKliyJb7/9FkDaMWz16tWwt7dHlSpVsHDhQmVfFxUVlWl9Lyz9isyl0AcXIHdnP9evX8f48eMNJj/UDy9//PEHatWqhQ4dOgBI63zZtm1b2Nvb52pHov+F1B9jgv5rgm7fvj1OnTqFpUuXGowiWxi3r9VqlaCrCy8pKSno06cPihUrhiNHjijrFVZ5XSf1d+bnz5/HzZs3MWzYsFw3MV+6dAlTpkxBv379sGXLFoO/6d7vw4cPo2rVqhn6ouX3GWNO+oZMnjw51/1+8pr+d0r3uswlP/aVJ0+ehJWVlfI9TUhIwNtvv406derA0tJSCQX6jhw5gvfffz/LW8jzKqxfu3YtQ3lWrFiBbt26KYHi9ddfh6enJ6ZNm6ZMejtu3LhMQ1Jh6VdkTqoILjk9+wHSOiJVqVIFwcHBBhVQ92ElJSXh22+/Ra1atZTZPk+ePGnyuB6Z0e/ZrxsOntJcunQJr7zyijJ9g7lnFM3t9mNiYjKM/pmUlIQ///wT1apVQ/369ZXLRn379kXJkiUN7hAprPK6Tl66dAkdO3ZEkyZNEBkZmevgEB4ejnLlyikjzGo0GoMzWp1r167B3t4e33//vcHygmgFM6VvSEBAQI7GbSkIuu9U06ZNTRoZ1xh5WS/Dw8NRokQJTJ48GcB/+/74+HgsWLAAPj4+mDRpktLXyJROtOYO67du3VJaRnQjPQNpo/f6+PggJiYGw4cPx0svvaQMBnf79m189dVXmfalKyz9isxNFcElt8LCwlC/fn0MHTo00/Dy5MkTVKpUKcNAdOaQm579L7oLFy6gS5cuJg+Kldfbj46ORunSpaHRaNCqVSvMmDEDP//8s8EAW/Xq1UPdunWVlpegoCBUqFAhw2yxhVFe18nIyEj07Nkz0wGvTKE/No9uGof+/ftjxowZGQ4u//77Lxo2bGgwf1lBMmffkMIkMjISvXr1yvVnm5m8qJdnzpyBvb19hss9urtvdC0vTZo0wbhx45TwYkrgNmdY/+WXX6DRaNCoUSN069bNoD737t0bGo0G5cuXN7iCkJ3C0q/I3IpEcAHShj3WhRf9u5SSk5Px7NkzdOrUSWneM/dZWmHs+1BY5PXZZk62f/36ddSrVw8eHh5o2LAhBg4cCDs7O9SrVw9vvPEGduzYgZ07d6JGjRpo27YtgLQvvrkGVMsPeV0ncztOS1RUFJydndG/f3+D5T169EDz5s1Rp04dvPbaa9i9e7fyt/bt2+O1114DUDi+c3nRN6QwMNcYPJkx5+cWExODypUro0WLFgbLly1bhp49eyqXXXThxc/PD0OHDs3RyYe5wjoADB48GHXr1kXPnj3Rpk0bpYXx+PHj8PDwUAYoNPWy9IvUB7PIBBfgv/AyaNAgpb9DYmIi5syZA3d39ywHpqKi5/Lly+jevTu6du2KP/74A1FRUdi+fTuaN2+Oxo0bw97eXpmUUzfpJJnHs2fPsG/fPtSpUwd9+vRR7upbvHgx7O3tERISguXLl6NChQqoV6+eciKyfPnyQteymZd9Qyh7uu9wy5Yt8dFHHwFIqyOOjo746aefAPx38E9ISMC0adPQrl07k0dh18ltoNON0fL999/jzTffxI8//ogePXrA399fOalu1qwZ3njjjVw9D6D+PphFKrgAadc7/f39UatWLXTt2hU9evSAm5tbppNPUdF24cIFtG/fHu3atVNGZQWABw8eYOvWrQgJCUH9+vXNPthfUaY7q/znn3+wc+dOtGnTBn379sWYMWMyjM0TFhYGjUZj0BegMMrLviGUUVxcnHLn1vnz5zFw4EC0aNECXbp0QZkyZfDLL78YrK8fXtJPcZLXoqOj8fXXXxssi42NRa1atfDBBx8gNjYWPXr0gJ+fH/bt24cjR47AyckJe/bsyfVzq7kPZpELLkBaM/R7772H3r17Y9GiRWYfGp9eHLrOd+3bt8+0863am/0LE10nyjFjxijLPv/8c/j5+cHGxgZr164FkNaknZycjL///hteXl5Kp/rC3Lydl31D6D///PMPAgMDMX78eKVPmi68lCpVyuDSiP6lloK4G1C/L13Hjh2xY8cOpcXwu+++g7+/P2JjYxEREYEePXqgXbt2mDRpErp3757rYSN01NoHs0gGFyJT6Df3626nJPPSDdqmm7hSP4Ts378fLVu2RPfu3ZU7KYC0ASKrV6+umr5Fedk3hNIkJSVh1KhRaN68OUJCQpTwcuHCBQwcOBDNmzfHhg0blPULcviC69evo2HDhvD19UWDBg0wdOhQuLu7Y/369dixYwdeeeUV7N27FwBw7tw5BAQEYODAgWav74U58GdFAwBCRNm6fPmyTJw4UeLi4mTVqlXStGnTgi7SC+PGjRtSv359adOmjezYsUNZ/u6770pcXJwsXrxYvvjiC/noo4+kRIkS8u6778ru3btl7ty5cuzYMalXr17BFZ4KDa1WK5aWlpKcnCwhISFy6NAhCQgIkGnTpkmJEiUkIiJCli1bJpcuXZIhQ4bIkCFDCrrIcvnyZZk+fbqkpqbKgAEDxMLCQkJDQ6VkyZLy7bffSqNGjeTw4cNiY2MjkZGRUqJECalYsWJBF7vgFXRyIlILNvfnjawmydTvRAkAX3zxBQIDA1G+fHnY2NiYfewfUif9Iep1EhISMGnSJDRq1Mig5eX8+fMYMmQIPD09MwxmWFAuXLiAoKAgBAYG4uLFi3j69CmOHTuGV155BVu3bgWgzlaRvMQWFyITJCUliY2NTUEX44Vz+fJlGTt2rNjY2IiLi4t8++23sm3bNgkMDJTU1FSxsLAQEZFPP/1Utm3bJu+++654e3sXcKmpoEVERIi3t7e0atVKKlSoICNHjpTSpUuLh4eHaLVaWbBggezdu1cCAgJk+vTp4ujoKGfPnpX169fLlClTpHLlygX9EkQkrf6PGTNGRETmzJkjzZs3L+ASFW4MLkRUKFy6dEnGjBkjv//+uyxYsEAmTZqk/E0/vDx9+lSKFy9eUMWkQiI1NVX2798vHTt2lIoVK4qrq6vcv39fnjx5Iu3atZMWLVpIhw4dZNGiRXLz5k3x8fGRyZMnS4kSJQrlCYguvAOQWbNmiZ+fX0EXqdBicCGiQuPq1asyatQosbS0lJkzZyo7b91uSqPRFGTxqJA4e/asrF69WjZu3CibNm2SoUOHSmhoqDRq1Eji4uLks88+k7/++kuKFy8uGo1Gbty4IfHx8fL222/LlClTRKRw1iX2pTOORUEXgIhIp1q1avLBBx8IAFm4cKEcOXJERNIOMoXxQEP57/Tp0+Lj4yMVKlQQEZHBgwfLu+++KxMmTJBffvlFOnfuLF988YVcvXpVZs6cKf369RMXFxdxdXWVrl27Fuq6VKNGDVm+fLnSgkSZY4sLERU6PPOkzEREREiDBg1k2rRpMm/ePINLiKtXr5YJEybIwoULJTg4WEqVKqU87u7du2JnZydOTk4FVXSTFMZLWYWJVUEXgIgoPd2Z5+zZs3nmSSIicu7cOWndurVUqVJF5s2bJyJp/VxERCwsLGTcuHGi0Whk/PjxotFoZNSoUUpQcXFxKahi5whDS/bY4kJEhRbPPEkk7fJQs2bNpHHjxnLp0iXp1auXrF69WkTSxm/RaDRKy8t7770nkydPlunTp8vkyZPF0dGxIItOeYB9XIio0GJooRMnTkijRo1k6tSp8tNPP8ncuXPl888/l3HjxomIiKWlpQBQWl/Gjh0rCxYskA8++ECSk5MLsuiUR9jiQkREhdahQ4dk9+7dSgvLo0ePZMeOHRISEiL9+vXLsuXlwYMH4uzsXGDlprzD4EJERKoAQDQajTx+/Fi++OKLTMOLhYWFaDQaZV168bBzLhERqYIuiDg6Osqrr74qIiIhISFiaWkpK1euFEtLywzr0ouHwYWIiFRHF14sLCxk+PDhYmtrK4sXLy7oYlE+YHAhIiJVcnR0lN69e4u1tbX4+voWdHEon7CPCxERqRr7sxQtvB2aiIhUjaGlaGFwISIiItVgcCEiIiLVYHAhIiIi1WBwISIiItVgcCEiIiLVYHAhIiIi1WBwISIiItVgcCEiIiLVYHAhIiIi1WBwISIiItVgcCEiIiLVYHAhIiIi1fg/b4mGHwtKiMUAAAAASUVORK5CYII=",
      "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.25,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "89af6f64-194f-4f1c-a5ed-dafd7243d371",
   "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": 57,
   "id": "018bb814-e2b1-47a1-8331-22ec810dfeb3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save grpofiler results\n",
    "with open(\n",
    "    \"figures/pert_embedding_cluster_gprofiler/scETM_inhouse_pert_emb_gprofiler.pickle\",\n",
    "    \"wb\",\n",
    ") as handle:\n",
    "    pickle.dump(cluster_process, handle, protocol=pickle.HIGHEST_PROTOCOL)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "126853e8-6814-45b0-9420-6b822cb8f42c",
   "metadata": {},
   "source": [
    "## PR Curve+AUC\n",
    "- AUCPR using prior graph as binary label"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "b02b6a47-b384-457b-b9ed-7840807a880e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "               AUC\n",
      "StringDB  0.573643\n",
      "CORUM     0.520186\n",
      "                F1  Threshold\n",
      "StringDB  0.679641        0.0\n",
      "CORUM     0.630435        0.0\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/scETM_inhouse_aucpr.csv\")\n",
    "f1_df.to_csv(\"figures/pert_embedding_recall/scETM_inhouse_f1.csv\")"
   ]
  },
  {
   "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": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2024-12-30 15:49:36,590 [INFO] Parsing data files for GSEA.............................\n",
      "2024-12-30 15:49:36,656 [INFO] 17782 gene_sets have been filtered out when max_size=300 and min_size=10\n",
      "2024-12-30 15:49:36,657 [INFO] 1203 gene_sets used for further statistical testing.....\n"
     ]
    }
   ],
   "source": [
    "factor_to_go = factor_enrichment_gsea(adata, np.abs(adata.uns[\"gene_emb\"]), fdr=5e-2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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": null,
   "id": "4b785ba3-f9c1-4207-b784-1ffe2ba74c30",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save latent enrichment results\n",
    "with open(\n",
    "    \"figures/factor_enrichments/scETM_inhouse_factor_enrichment.pickle\", \"wb\"\n",
    ") as handle:\n",
    "    pickle.dump(filtered_factor_to_go, handle, protocol=pickle.HIGHEST_PROTOCOL)"
   ]
  }
 ],
 "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
}