inhouse_GSFA_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",
    "\n",
    "from utils import (\n",
    "    auprc,\n",
    "    factor_enrichment_gsea,\n",
    "    get_gprofiler,\n",
    "    inhouse_preprocess,\n",
    "    perturbation_signal_recovery,\n",
    "    read_aws_csv,\n",
    "    read_aws_h5ad,\n",
    "    read_aws_npz,\n",
    "    read_aws_pickle,\n",
    "    retrieve_stringdb_neighbors,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "39762015-fa94-4ae8-9038-405f7c56a03b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read learned parameters\n",
    "Z = read_aws_csv(\n",
    "    \"s3://pert-spectra/gsfa_checkpoints/inhouse_gsfa_outputs/Z.csv\"\n",
    ").to_numpy()\n",
    "W = read_aws_csv(\n",
    "    \"s3://pert-spectra/gsfa_checkpoints/inhouse_gsfa_outputs/W.csv\"\n",
    ").to_numpy()\n",
    "F = read_aws_csv(\n",
    "    \"s3://pert-spectra/gsfa_checkpoints/inhouse_gsfa_outputs/F.csv\"\n",
    ").to_numpy()\n",
    "beta = read_aws_csv(\n",
    "    \"s3://pert-spectra/gsfa_checkpoints/inhouse_gsfa_outputs/beta.csv\"\n",
    ").to_numpy()\n",
    "lsfr = read_aws_csv(\n",
    "    \"s3://pert-spectra/gsfa_checkpoints/inhouse_gsfa_outputs/lsfr.csv\"\n",
    ").to_numpy()\n",
    "# read gene and perturbation information\n",
    "gene_labels = read_aws_csv(\"s3://pert-spectra/gsfa_checkpoints/inhouse_top_genes.csv\")[\n",
    "    \"x\"\n",
    "].to_numpy()\n",
    "pert_labels = read_aws_npz(\"s3://pert-spectra/gsfa_checkpoints/inhouse_G_labels.npz\")[\n",
    "    \"arr_0\"\n",
    "]\n",
    "pert_labels = [x.split(\"_\")[1] for x in pert_labels[:-1]]\n",
    "pert_labels = pert_labels + [\"ctrl\"]\n",
    "# read inputs and preprocessed input\n",
    "preprocessed_Y = read_aws_npz(\n",
    "    \"s3://pert-spectra/gsfa_checkpoints/inhouse_GSFA_preprocessed.npz\"\n",
    ")[\"array1\"]\n",
    "G = read_aws_npz(\"s3://pert-spectra/gsfa_checkpoints/inhouse_GSFA_inputs.npz\")[\"array2\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "601bbfb0-17ed-4694-850c-fcbf1fef961b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in adata of raw data for reference\n",
    "adata = read_aws_h5ad(\"path to raw inhouse h5ad\")\n",
    "adata = adata[:, gene_labels - 1]\n",
    "adata = inhouse_preprocess(adata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a6ea3e4a-7c2f-41ba-b136-9ef9a6ab3029",
   "metadata": {},
   "outputs": [],
   "source": [
    "# filter adata to perturbations with at least 50 samples for each treatment\n",
    "adata.obs[\"condition\"] = adata.obs[\"condition\"].astype(str)\n",
    "adata.obs[\"Treatment\"] = adata.obs[\"Treatment\"].astype(str)\n",
    "adata.obs[\"pert_treat\"] = adata.obs[\"condition\"] + \"+\" + adata.obs[\"Treatment\"]\n",
    "obs_df = pd.DataFrame(adata.obs[\"pert_treat\"])\n",
    "category_counts = obs_df[\"pert_treat\"].value_counts()\n",
    "filtered_categories = category_counts[category_counts >= 50].index\n",
    "adata = adata[adata.obs[\"pert_treat\"].isin(filtered_categories)]"
   ]
  },
  {
   "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 (on DE genes if available)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "e2885bd3-5433-47d6-88d5-47be83df5b5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# for GSFA, measure reconstruction on training\n",
    "from scipy.stats import spearmanr\n",
    "\n",
    "reconstruction = Z @ W.T\n",
    "# rebuild pert labels\n",
    "G_labels = []\n",
    "for row in G:\n",
    "    pert_idx = np.where(row == 1)[0].tolist()\n",
    "    perts = [pert_labels[i] for i in pert_idx]\n",
    "    G_labels.append(\"+\".join(perts))\n",
    "test_corr_singles = []\n",
    "test_corr_combos = []\n",
    "\n",
    "for pert in set(G_labels):\n",
    "    hold_idx = [i for i, x in enumerate(G_labels) if x == pert]\n",
    "    if not hold_idx:\n",
    "        continue\n",
    "    recon = reconstruction[hold_idx]\n",
    "    # correlation\n",
    "    mean_reconstruction = recon.mean(axis=0)\n",
    "    mean_observed = preprocessed_Y[hold_idx].mean(axis=0)\n",
    "    if (\"+\" in pert) and (\"ctrl\" not in pert):\n",
    "        test_corr_combos.append(\n",
    "            [pert, spearmanr(mean_reconstruction, mean_observed)[0]]\n",
    "        )\n",
    "    else:\n",
    "        test_corr_singles.append(\n",
    "            [pert, spearmanr(mean_reconstruction, mean_observed)[0]]\n",
    "        )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "971dead3-70d4-41a5-98c0-ab485b32da4f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAncAAAHFCAYAAACQIN+BAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABJE0lEQVR4nO3deZzNdf//8edhxpnVMMZszJYYu2iZohhkK0pcLUSWconkKqkuV2lGC1GkktTVZamEukIqjZD1soQra0N00SgGYxvLGMa8f3/0m/N1zAxnmDNn5uNxv93O7ea8P9vr/fl85n2ePud8zrEZY4wAAABgCeU8XQAAAACKD+EOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYyBWFu82bN6tPnz6Ki4uTj4+PAgIC1KRJE40ZM0ZHjhwp7hqvytKlS2Wz2bR06dIiL/vzzz8rOTlZe/bsyTetd+/eio2Nver6imrPnj2y2Wx68803i22d+/btU3JysjZu3Fhs60R+U6dOlc1mczqfEhMTlZiYWOR1jRw5UnPnzi222vLknV9Tp04t9nVfqeL+W0tMTFT9+vWLbX1lRWk8tqVVcnKybDabp8u4Jrj79ScvA/z73/8utnWWxmxwsSKHu3/+85+68cYbtW7dOj377LNKSUnRnDlzdP/992vSpEl69NFH3VGnR/z8888aMWJEgQdw+PDhmjNnTskX5Qb79u3TiBEjCHceMHHiRE2cOLHIy7kr3AGQHnvsMa1evdrTZVwTyuLrT1nIBl5FmXn16tUaMGCA2rRpo7lz58putzumtWnTRs8884xSUlKKpbDTp0/Lz88vX/v58+eVk5PjtG1PqFGjhke3j5LjznOubt26xb5OlA2lZSwr7bKysuTr61ui26xevbqqV69eotssSaXh3MuroayuvzClJRsU6crdyJEjZbPZ9OGHHxZ4UlSoUEH33HOP43lubq7GjBmj2rVry263KzQ0VI888oh+//13p+Xy3iZZvny5mjZtKj8/P/Xt29fxNsKYMWP06quvKi4uTna7XUuWLJEkrV+/Xvfcc4+Cg4Pl4+Ojxo0b6/PPP79sP9avX6+HHnpIsbGx8vX1VWxsrLp166bffvvNMc/UqVN1//33S5Jatmwpm83m9JZGQZdez5w5o2HDhikuLk4VKlRQtWrV9MQTT+jYsWNO88XGxqpjx45KSUlRkyZN5Ovrq9q1a2vy5MmXrf3Cffvaa68pOjpaPj4+uummm7R48eJ88+3cuVPdu3dXaGio7Ha76tSpo/fee88xfenSpbr55pslSX369HH0Mzk5Wd9++61sNpvWrVvnmP/LL7+UzWbT3Xff7bSdhg0bqmvXro7nxhhNnDhRN9xwg3x9fVW5cmX95S9/0f/+9798NS5atEitW7dWxYoV5efnp2bNmuXrS97bJNu2bVO3bt0UFBSksLAw9e3bV8ePH7/s/so7x1asWKFbb71Vvr6+qlatmoYPH67z58875iuuc27NmjVq1qyZfHx8FBkZqWHDhuncuXMF1nXx27LZ2dl6+eWXVadOHfn4+KhKlSpq2bKlVq1aJUmy2Ww6deqUpk2b5jheF64jPT1d/fv3V/Xq1VWhQgXFxcVpxIgR+Qa6ffv26YEHHlBgYKCCgoL04IMPKj09/bL7UpIOHTqkgQMHqm7dugoICFBoaKhatWqlFStWOM134ccIxo0bp7i4OAUEBOi2227TmjVr8q136tSpio+Pd5yrH3/8sUv1SK6PN3kudy5I0vvvv69GjRopICBAgYGBql27tv7xj384zePK/i7svPr8889VoUIFDR8+PF9927dvl81m0zvvvFOkbUlXd2xd2Y4xRnfddZeqVKmitLQ0x3KnT59WvXr1VKdOHZ06dUrS//3t/vTTT+rSpYsqVqyooKAg9ejRQ4cOHXLabt7YOHv2bDVu3Fg+Pj4aMWJEkfp+uWN2+vRpDR061PGxouDgYN10002aMWOGY56C3pYt6uvZunXrdMcdd8jPz0/XXXedXn/9deXm5l5239tsNg0aNEgffPCBatWqJbvdrrp162rmzJlFPk7S5cc0T9dQ2OtP3r4s6GMrF7/+utLHM2fOaMiQIQoPD5evr69atGihn376yWm9VsgGDsZFOTk5xs/PzyQkJLi6iPnrX/9qJJlBgwaZlJQUM2nSJFO1alUTFRVlDh065JivRYsWJjg42ERFRZl3333XLFmyxCxbtszs3r3bSDLVqlUzLVu2NP/+97/N999/b3bv3m1++OEHU6FCBXPHHXeYWbNmmZSUFNO7d28jyUyZMsWx7iVLlhhJZsmSJY62L774wrz00ktmzpw5ZtmyZWbmzJmmRYsWpmrVqo66Dh48aEaOHGkkmffee8+sXr3arF692hw8eNAYY0yvXr1MTEyMY525ubmmXbt2xsvLywwfPtx8//335s033zT+/v6mcePG5syZM455Y2JiTPXq1U3dunXNxx9/bBYsWGDuv/9+I8ksW7bskvs0b59ERUWZ22+/3Xz55Zfmiy++MDfffLPx9vY2q1atcsy7bds2ExQUZBo0aGA+/vhj8/3335tnnnnGlCtXziQnJxtjjDl+/LiZMmWKkWRefPFFRz/37t1rTpw4Yby9vc3IkSMd63z88ceNr6+v8ff3N2fPnjXGGHPgwAFjs9nMxIkTHfP169fPeHt7m2eeecakpKSYzz77zNSuXduEhYWZ9PR0x3yffPKJsdlspnPnzmb27Nnm66+/Nh07djTly5c3ixYtcsyXlJRkJJn4+Hjz0ksvmYULF5px48YZu91u+vTpc8l9Zsyf51iVKlVMZGSkeeedd8yCBQvM4MGDjSTzxBNP5Nu/V3PObdu2zfj5+Zm6deuaGTNmmK+++sq0a9fOREdHG0lm9+7dTnW1aNHC8fzcuXOmZcuWxsvLywwdOtTMnz/fzJs3z/zjH/8wM2bMMMYYs3r1auPr62vuuusux/Hatm2bMcaY/fv3m6ioKBMTE2M++OADs2jRIvPKK68Yu91uevfu7djO6dOnTZ06dUxQUJB59913Hfsjr8YL+1OQ7du3mwEDBpiZM2eapUuXmm+++cY8+uijply5ck5/a3n7MzY21rRv397MnTvXzJ071zRo0MBUrlzZHDt2zDFv3nl47733mq+//tp8+umn5vrrr3f053KKMt64ci7MmDHDSDJPPvmk+f77782iRYvMpEmTzODBgx3zuLq/L3Ve3XfffSYqKsqcP3/eqT/PPfecqVChgsnIyCjRY+vqdjIyMkz16tVNQkKCYyzo1auX8fX1NZs3b3bMl/e3GxMTY5599lmzYMECM27cOMfYmLesMX+OjREREea6664zkydPNkuWLDE//vijyzW5csz69+9v/Pz8zLhx48ySJUvMN998Y15//XXz7rvv5qv5QkU9v2rWrGkmTZpkFi5caAYOHGgkmWnTpl1y3xtjHON73vgxb9480759eyPJfPHFF0U+Tpc69zxdw6ZNmwp9/cnblxeOj3kufv29VB/zMkBUVFS+saVixYrm119/daynLGeDfMfQ1RnT09ONJPPQQw+5NH9qaqqRZAYOHOjUvnbtWiPJ/OMf/3C0tWjRwkgyixcvdpo374DVqFHDaQAwxpjatWubxo0bm3Pnzjm1d+zY0URERDgGyoLC3cVycnLMyZMnjb+/v3n77bcd7V988UWhy158AFNSUowkM2bMGKf5Zs2aZSSZDz/80NEWExNjfHx8zG+//eZoy8rKMsHBwaZ///6F1mnM/+2TyMhIk5WV5WjPzMw0wcHB5s4773S0tWvXzlSvXt0cP37caR2DBg0yPj4+5siRI8YYY9atW1fooH/77bebVq1aOZ5ff/315tlnnzXlypVznGzTp083kswvv/xijPkzfEgyY8eOdVrX3r17ja+vr3nuueeMMcacOnXKBAcHm06dOjnNd/78edOoUSNzyy23ONryBtuL9+/AgQONj4+Pyc3NveR+yzvHvvrqK6f2fv36mXLlyjmORXGccw8++KDx9fV1CrE5OTmmdu3alw13H3/8sZFk/vnPf16yP/7+/qZXr1752vv3728CAgKczi1jjHnzzTeNJEcIfP/99wvdH64EgIvl5OSYc+fOmdatW5v77rvP0Z63Pxs0aGBycnIc7T/++KOR5Ais58+fN5GRkaZJkyZOx3LPnj3G29v7suHuSsaby50LgwYNMpUqVbrkdl3d35c6r+bNm2ckme+//97RlpOTYyIjI03Xrl2LvK2rPbaubscYY1auXGm8vLzMU089ZSZPnmwkmY8++shpuby/3aefftqpPW/c+PTTTx1tMTExpnz58mbHjh1XVJMrx6x+/fqmc+fOl5zn4nB3JefX2rVrneatW7euadeu3SW3a8yfwaqw8eP66693tBXHuVcaarjU609Rw11B68/LAIWNLY899lih+6EsZYOLue2rUPIuh/bu3dup/ZZbblGdOnXyve1WuXJltWrVqsB13XPPPfL29nY837Vrl7Zv366HH35YkpSTk+N43HXXXdq/f7927NhRaG0nT57U888/r+uvv15eXl7y8vJSQECATp06pdTU1Cvprn744QdJ+ft7//33y9/fP19/b7jhBkVHRzue+/j4qFatWk6Xfy+lS5cu8vHxcTwPDAxUp06dtHz5cp0/f15nzpzR4sWLdd9998nPzy/fPjpz5kyBb4tdrHXr1vrPf/6jrKws/fbbb9q1a5ceeugh3XDDDVq4cKGkP99WjY6OVs2aNSVJ33zzjWw2m3r06OG03fDwcDVq1Mhx5/KqVat05MgR9erVy2m+3NxctW/fXuvWrXO8tZPnwrf9pT/fDj5z5owOHjx42b4EBgbmW7579+7Kzc3V8uXL823nSs+5JUuWqHXr1goLC3MsX758eT344IOXrfG7776Tj4+P+vbte9l5C/LNN9+oZcuWioyMdKqxQ4cOkqRly5Y5aixsf7hq0qRJatKkiXx8fOTl5SVvb28tXry4wL+hu+++W+XLl3c8b9iwoSQ5zvcdO3Zo37596t69u9PbYTExMWratOllaynqeOPKuXDLLbfo2LFj6tatm7766itlZGTk266r+zvPxeeVJHXo0EHh4eGaMmWKo23BggXat2+f03lQUse2KH1q1qyZXnvtNY0fP14DBgxQjx49Cr2pLu9vJ88DDzwgLy+vfG8PNmzYULVq1bqimlw5Zrfccou+++47/f3vf9fSpUuVlZV12X1S1PMrPDxct9xyS75+uTq+FzZ+7Nq1y/E2cHGce6W9hqK61PoLG1suPP+skA3yuBzuQkJC5Ofnp927d7s0/+HDhyVJERER+aZFRkY6pucpaL7Cph04cECSNHToUHl7ezs9Bg4cKEkF/lHn6d69uyZMmKDHHntMCxYs0I8//qh169apatWqLv2hF+Tw4cPy8vJS1apVndptNpvCw8Pz9bdKlSr51mG3213efnh4eIFtZ8+e1cmTJ3X48GHl5OTo3XffzbeP7rrrLkmX3kd57rzzTmVnZ2vlypVauHChQkJC1LhxY915551atGiRJGnx4sW68847HcscOHBAxhiFhYXl2/aaNWsc2807jn/5y1/yzTd69GgZY/J9tc7F+y3vs5+u7LcLB6oL95mky56PRTnnDh8+XOjxuZxDhw4pMjJS5cpd2f+7Dhw4oK+//jpfjfXq1ctX46X2x+WMGzdOAwYMUEJCgr788kutWbNG69atU/v27Qs8Fpc7bnn7/0r3W1HHG1fOhZ49e2ry5Mn67bff1LVrV4WGhiohIcHxnxrJ9f2dp6D6vLy81LNnT82ZM8fxGZypU6cqIiJC7dq1K/K2rvbYFrVPDz/8sCpUqKDs7Gw9++yzha734u17eXmpSpUqLr0WuFqTK8fsnXfe0fPPP6+5c+eqZcuWCg4OVufOnbVz585Cay/q+eWu8f3CWorj3CvtNRTVpdZfWH8uPHZWyAZ5XL5btnz58mrdurW+++47/f7775e9kyivwP379+ebd9++fQoJCXFqu9R3Cl08LW/ZYcOGqUuXLgUuEx8fX2D78ePH9c033ygpKUl///vfHe3Z2dlX9R19VapUUU5Ojg4dOuR0EI0xSk9Pd3xotLgU9OHo9PR0VahQQQEBAfL29lb58uXVs2dPPfHEEwWuIy4u7rLbSUhIUEBAgBYtWqQ9e/aodevWstlsat26tcaOHat169YpLS3NKdyFhITIZrNpxYoVBd54k9eWdxzfffdd3XrrrQVuv6AXqSuVF9AulLcfL/6DuppzrkqVKoUen8upWrWqVq5cqdzc3CsKeCEhIWrYsKFee+21AqdHRkY6avzxxx+vqEZJ+vTTT5WYmKj333/fqf3EiRNFrFiOegrbvis1FXW8cfVc6NOnj/r06aNTp05p+fLlSkpKUseOHfXLL78oJibG5f2dp7Bxrk+fPnrjjTc0c+ZMPfjgg5o3b56eeuopp6udJXVsi9Kn8+fP6+GHH1blypVlt9v16KOP6j//+Y8qVKhQ4ParVavmeJ6Tk6PDhw9f9m+vqDVd7pj5+/trxIgRGjFihA4cOOC4itepUydt3769wPUX9fy6Wpf6O8irpbjOvdJcg4+PT4E3zBV2YeJS6y+sP3l9sUo2yFOkr0IZNmyY5s+fr379+umrr77K9wd87tw5paSkqFOnTo63WD/99FOn4tetW6fU1FS98MILV1x0fHy8atasqU2bNmnkyJFFWtZms8kYky90fPTRR/nulCvKVaHWrVtrzJgx+vTTT/X000872r/88kudOnVKrVu3LlKdlzN79my98cYbjrdmT5w4oa+//lp33HGHypcvLz8/P7Vs2VI//fSTGjZsWOBgm+dS/fT29lbz5s21cOFC7d27V6+//rok6Y477pCXl5defPFFR9jL07FjR73++uv6448/9MADDxS63WbNmqlSpUr6+eefNWjQoCvaD0Vx4sQJzZs3z+ntqs8++0zlypVT8+bNL7lsUc65li1bat68eTpw4IAjnJ4/f16zZs26bI0dOnTQjBkzNHXq1Eu+NVvY/+Q6duyo+fPnq0aNGqpcufIla/z8888L3B+usNls+f6GNm/erNWrVysqKsqldVwoPj5eERERmjFjhoYMGeIYpH/77TetWrUq34vExYo63hT1XPD391eHDh109uxZde7cWdu2bVNMTIzL+/ty6tSpo4SEBE2ZMkXnz59Xdna2+vTp4zRPSR3bovQpKSlJK1as0Pfffy9/f381b95czz77rN5+++18806fPl033nij4/nnn3+unJwcl77E+0r2c2HH7EJhYWHq3bu3Nm3apPHjxxf6FVzufD0ryOLFiwscP2rUqOEIl8V17nm6hku9/sTGxuqLL75Qdna2Y77Dhw9r1apVqlixYpG2U9jY8sgjj0iyTjbIU6Rwd9ttt+n999/XwIEDdeONN2rAgAGqV6+ezp07p59++kkffvih6tevr06dOik+Pl5//etf9e6776pcuXLq0KGD9uzZo+HDhysqKsqpk1figw8+UIcOHdSuXTv17t1b1apV05EjR5Samqr//ve/+uKLLwpcrmLFimrevLneeOMNhYSEKDY2VsuWLdO//vUvVapUyWnevG+x//DDDxUYGCgfHx/FxcUVeNm0TZs2ateunZ5//nllZmaqWbNm2rx5s5KSktS4cWP17Nnzqvp7sfLly6tNmzYaMmSIcnNzNXr0aGVmZjq+NkCS3n77bd1+++264447NGDAAMXGxurEiRPatWuXvv76a8dnAWrUqCFfX19Nnz5dderUUUBAgCIjIx0vqK1bt9YzzzwjSY4rdL6+vmratKm+//57NWzYUKGhoY7tNmvWTH/961/Vp08frV+/Xs2bN5e/v7/279+vlStXqkGDBhowYIACAgL07rvvqlevXjpy5Ij+8pe/KDQ0VIcOHdKmTZt06NChfFeGrkaVKlU0YMAApaWlqVatWpo/f77++c9/asCAAU6fcSiMq+fciy++qHnz5qlVq1Z66aWX5Ofnp/feey/f5wcL0q1bN02ZMkWPP/64duzYoZYtWyo3N1dr165VnTp19NBDD0mSGjRooKVLl+rrr79WRESEAgMDFR8fr5dfflkLFy5U06ZNNXjwYMXHx+vMmTPas2eP5s+fr0mTJql69ep65JFH9NZbb+mRRx7Ra6+9ppo1a2r+/PlasGCBS/uyY8eOeuWVV5SUlKQWLVpox44devnllxUXF3dF3y1Vrlw5vfLKK3rsscd03333qV+/fjp27JiSk5NdejuxqOONK+dCv3795Ovrq2bNmikiIkLp6ekaNWqUgoKCHC/wru5vV/Tt21f9+/fXvn371LRp03zvPpTUsXV1OwsXLtSoUaM0fPhwxwvUqFGjNHToUCUmJuq+++5zWu/s2bPl5eWlNm3aaNu2bRo+fLgaNWp0yf8AFrUmV45ZQkKCOnbsqIYNG6py5cpKTU3VJ598ottuu63AYCcV/fy6WiEhIWrVqpWGDx8uf39/TZw4Udu3b3f6KpLiPPc8WcOlXn969uypDz74QD169FC/fv10+PBhjRkzpsjBTpIOHjzoGFuOHz+upKQk+fj4aNiwYZKskw0cinT7xf+3ceNG06tXLxMdHW0qVKjguKX3pZdectwObMyfd8CNHj3a1KpVy3h7e5uQkBDTo0cPx23OeVq0aGHq1auXbzt5d8C88cYbBdaxadMm88ADD5jQ0FDj7e1twsPDTatWrcykSZMc8xR0t+zvv/9uunbtaipXrmwCAwNN+/btzdatW01MTEy+OxDHjx9v4uLiTPny5Z3u6Ln4jhhj/ryr5fnnnzcxMTHG29vbREREmAEDBpijR486zRcTE2PuvvvufP0p7M6ggvbJ6NGjzYgRI0z16tVNhQoVTOPGjc2CBQsKnL9v376mWrVqxtvb21StWtU0bdrUvPrqq07zzZgxw9SuXdt4e3sbSSYpKckxbdOmTUaSqVmzptMyr732mpFkhgwZUmCtkydPNgkJCcbf39/4+vqaGjVqmEceecSsX7/eab5ly5aZu+++2wQHBxtvb29TrVo1c/fddzvdcp9399qFXzlgzP99fcalbus35v/OsaVLl5qbbrrJ2O12ExERYf7xj3843f1aHOecMcb85z//Mbfeequx2+0mPDzcPPvss+bDDz+87N2yxvx5Hr300kumZs2apkKFCqZKlSqmVatWTl9zs3HjRtOsWTPj5+dnJDmt49ChQ2bw4MEmLi7OeHt7m+DgYHPjjTeaF154wZw8edIxX97fQUBAgAkMDDRdu3Y1q1atcumOyuzsbDN06FBTrVo14+PjY5o0aWLmzp1b6F1sBe3Pi88zY4z56KOPHP2uVauWmTx5coF/awUp6nhzuXNh2rRppmXLliYsLMxUqFDBREZGmgceeMDpaz6McW1/X+68MubPryXy9fW95N3SJXFsXdnOvn37TGhoqGnVqpXTV7jk5uaaTp06mUqVKjnO87y/3Q0bNphOnTo5aurWrZs5cOCA03YLGxtd7bsrx+zvf/+7uemmm0zlypWN3W431113nXn66acdXzlzYc0XutrXM1fPY/3/r+SZOHGiqVGjhvH29ja1a9c206dPv6J94sq55+kaLvX6M23aNFOnTh3j4+Nj6tata2bNmlWkcSYvA3zyySdm8ODBpmrVqsZut5s77rgj32tRWc4GF7MZY4x7YiNQeiQmJiojI0Nbt271dCnANSU5OVkjRozQoUOHiv2zaVZks9n0xBNPaMKECdd0Dbg6bvsqFAAAAJQ8wh0AAICF8LYsAACAhXDlDgAAwEIIdwAAABZCuAMAALCQIn2JsdXl5uZq3759CgwMLPLPpAAAAM8wxujEiRNX9dvcVkK4u8C+ffuu6KeTAACA5+3du/eqfpnDKgh3FwgMDJT058lxJT9vAgAASl5mZqaioqIcr+PXOsLdBfLeiq1YsSLhDgCAMoaPVP2JN6YBAAAshHAHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABbi5ekCAKCkpKWlKSMjw9NlFElISIiio6M9XQaAMoRwB+CakJaWptq16ygr67SnSykSX18/bd+eSsAD4DLCHYBrQkZGhrKyTiuhb5IqRsR6uhyXZO7fo7WTRygjI4NwB8BlhDsA15SKEbEKjo73dBkA4DbcUAEAAGAhhDsAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsJAyEe5GjRqlm2++WYGBgQoNDVXnzp21Y8cOp3l69+4tm83m9Lj11ls9VDEAAIBnlIlwt2zZMj3xxBNas2aNFi5cqJycHLVt21anTp1ymq99+/bav3+/4zF//nwPVQwAAOAZXp4uwBUpKSlOz6dMmaLQ0FBt2LBBzZs3d7Tb7XaFh4eXdHkAAAClRpkIdxc7fvy4JCk4ONipfenSpQoNDVWlSpXUokULvfbaawoNDS10PdnZ2crOznY8z8zMdE/BgMWkpaUpIyPD02UUSWpqqqdLAIASUebCnTFGQ4YM0e2336769es72jt06KD7779fMTEx2r17t4YPH65WrVppw4YNstvtBa5r1KhRGjFiREmVDlhCWlqaateuo6ys054u5Yqcyz7r6RIAwK3KXLgbNGiQNm/erJUrVzq1P/jgg45/169fXzfddJNiYmL07bffqkuXLgWua9iwYRoyZIjjeWZmpqKiotxTOGARGRkZyso6rYS+SaoYEevpcly2f8tqbZ33oXJycjxdCgC4VZkKd08++aTmzZun5cuXq3r16pecNyIiQjExMdq5c2eh89jt9kKv6gG4tIoRsQqOjvd0GS7L3L/H0yUAQIkoE+HOGKMnn3xSc+bM0dKlSxUXF3fZZQ4fPqy9e/cqIiKiBCoEAAAoHcrEV6E88cQT+vTTT/XZZ58pMDBQ6enpSk9PV1ZWliTp5MmTGjp0qFavXq09e/Zo6dKl6tSpk0JCQnTfffd5uHoAAICSUyau3L3//vuSpMTERKf2KVOmqHfv3ipfvry2bNmijz/+WMeOHVNERIRatmypWbNmKTAw0AMVAwAAeEaZCHfGmEtO9/X11YIFC0qoGgAAgNKrTLwtCwAAANcQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYSJm4WxYArmWpqameLqFIQkJCFB0d7ekygGsW4Q4ASqms44cl2dSjRw9Pl1Ikvr5+2r49lYAHeAjhDgBKqXOnT0gyuqH786oaV9vT5bgkc/8erZ08QhkZGYQ7wEMIdwBQygWERis4Ot7TZQAoI7ihAgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALCQMhHuRo0apZtvvlmBgYEKDQ1V586dtWPHDqd5jDFKTk5WZGSkfH19lZiYqG3btnmoYgAAAM8oE+Fu2bJleuKJJ7RmzRotXLhQOTk5atu2rU6dOuWYZ8yYMRo3bpwmTJigdevWKTw8XG3atNGJEyc8WDkAAEDJ8vJ0Aa5ISUlxej5lyhSFhoZqw4YNat68uYwxGj9+vF544QV16dJFkjRt2jSFhYXps88+U//+/T1RNgAAQIkrE1fuLnb8+HFJUnBwsCRp9+7dSk9PV9u2bR3z2O12tWjRQqtWrSp0PdnZ2crMzHR6AAAAlGVlLtwZYzRkyBDdfvvtql+/viQpPT1dkhQWFuY0b1hYmGNaQUaNGqWgoCDHIyoqyn2FAwAAlIAyF+4GDRqkzZs3a8aMGfmm2Ww2p+fGmHxtFxo2bJiOHz/ueOzdu7fY6wUAAChJZeIzd3mefPJJzZs3T8uXL1f16tUd7eHh4ZL+vIIXERHhaD948GC+q3kXstvtstvt7isYAACghJWJK3fGGA0aNEizZ8/WDz/8oLi4OKfpcXFxCg8P18KFCx1tZ8+e1bJly9S0adOSLhcAAMBjysSVuyeeeEKfffaZvvrqKwUGBjo+RxcUFCRfX1/ZbDY99dRTGjlypGrWrKmaNWtq5MiR8vPzU/fu3T1cPQAAQMkpE+Hu/ffflyQlJiY6tU+ZMkW9e/eWJD333HPKysrSwIEDdfToUSUkJOj7779XYGBgCVcLAADgOWUi3BljLjuPzWZTcnKykpOT3V8QAABAKVUmPnMHAAAA1xDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhXh5ugDgWpeWlqaMjAxPl+Gy1NRUT5cAALgEwh3gQWlpaapdu46ysk57upQiO5d91tMlAAAKQLgDPCgjI0NZWaeV0DdJFSNiPV2OS/ZvWa2t8z5UTk6Op0sBABSAcAeUAhUjYhUcHe/pMlySuX+Pp0sAAFwCN1QAAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIW4Ndxdd911Onz4cL72Y8eO6brrrnPnpgEAAK5Jbg13e/bs0fnz5/O1Z2dn648//nDnpgEAAK5JXu5Y6bx58xz/XrBggYKCghzPz58/r8WLFys2NtYdmwYAALimuSXcde7cWZJks9nUq1cvp2ne3t6KjY3V2LFj3bFpAACAa5pbwl1ubq4kKS4uTuvWrVNISIg7NgMAAICLuCXc5dm9e7c7Vw8AAICLuDXcSdLixYu1ePFiHTx40HFFL8/kyZPdvXkAAIBrilvD3YgRI/Tyyy/rpptuUkREhGw2mzs3BwAAcM1za7ibNGmSpk6dqp49e7pzMwAAAPj/3Po9d2fPnlXTpk3duQkAAABcwK3h7rHHHtNnn31WLOtavny5OnXqpMjISNlsNs2dO9dpeu/evWWz2Zwet956a7FsGwAAoKxw69uyZ86c0YcffqhFixapYcOG8vb2dpo+btw4l9d16tQpNWrUSH369FHXrl0LnKd9+/aaMmWK43mFChWurHAAAIAyyq3hbvPmzbrhhhskSVu3bnWaVtSbKzp06KAOHTpcch673a7w8PAirRcAAMBK3BrulixZ4s7V57N06VKFhoaqUqVKatGihV577TWFhoYWOn92drays7MdzzMzM0uiTAAAALdx62fuSlKHDh00ffp0/fDDDxo7dqzWrVunVq1aOYW3i40aNUpBQUGOR1RUVAlWDAAAUPzceuWuZcuWl3z79Ycffii2bT344IOOf9evX1833XSTYmJi9O2336pLly4FLjNs2DANGTLE8TwzM5OABwAAyjS3hru8z9vlOXfunDZu3KitW7eqV69e7ty0IiIiFBMTo507dxY6j91ul91ud2sdAAAAJcmt4e6tt94qsD05OVknT55056Z1+PBh7d27VxEREW7dDgAAQGnikc/c9ejRo8i/K3vy5Elt3LhRGzdulCTt3r1bGzduVFpamk6ePKmhQ4dq9erV2rNnj5YuXapOnTopJCRE9913nxt6AAAAUDq59cpdYVavXi0fH58iLbN+/Xq1bNnS8Tzvs3K9evXS+++/ry1btujjjz/WsWPHFBERoZYtW2rWrFkKDAws1toBAABKM7eGu4tvZDDGaP/+/Vq/fr2GDx9epHUlJibKGFPo9AULFlxRjQAAAFbi1nAXFBTk9LxcuXKKj4/Xyy+/rLZt27pz0wAAANckt4a7C38KDAAAAO5XIp+527Bhg1JTU2Wz2VS3bl01bty4JDYLAABwzXFruDt48KAeeughLV26VJUqVZIxRsePH1fLli01c+ZMVa1a1Z2bBwAAuOa49atQnnzySWVmZmrbtm06cuSIjh49qq1btyozM1ODBw9256YBAACuSW69cpeSkqJFixapTp06jra6devqvffe44YKALCw1NRUT5dQJCEhIYqOjvZ0GUCxcGu4y83Nlbe3d752b29v5ebmunPTAAAPyDp+WJJNPXr08HQpReLr66ft21MJeLAEt4a7Vq1a6W9/+5tmzJihyMhISdIff/yhp59+Wq1bt3bnpgEAHnDu9AlJRjd0f15V42p7uhyXZO7fo7WTRygjI4NwB0twa7ibMGGC7r33XsXGxioqKko2m01paWlq0KCBPv30U3duGgDgQQGh0QqOjvd0GcA1ya3hLioqSv/973+1cOFCbd++XcYY1a1bV3feeac7NwsAAHDNcsvdsj/88IPq1q2rzMxMSVKbNm305JNPavDgwbr55ptVr149rVixwh2bBgAAuKa5JdyNHz9e/fr1U8WKFfNNCwoKUv/+/TVu3Dh3bBoAAOCa5pZwt2nTJrVv377Q6W3bttWGDRvcsWkAAIBrmlvC3YEDBwr8CpQ8Xl5eOnTokDs2DQAAcE1zS7irVq2atmzZUuj0zZs3KyIiwh2bBgAAuKa5Jdzdddddeumll3TmzJl807KyspSUlKSOHTu6Y9MAAADXNLd8FcqLL76o2bNnq1atWho0aJDi4+Nls9mUmpqq9957T+fPn9cLL7zgjk0DAABc09wS7sLCwrRq1SoNGDBAw4YNkzFGkmSz2dSuXTtNnDhRYWFh7tg0AADANc1tX2IcExOj+fPn6+jRo9q1a5eMMapZs6YqV67srk0CAABc89z6CxWSVLlyZd18883u3gwAAADkphsqAAAA4BmEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALKTMhLvly5erU6dOioyMlM1m09y5c52mG2OUnJysyMhI+fr6KjExUdu2bfNMsQAAAB5SZsLdqVOn1KhRI02YMKHA6WPGjNG4ceM0YcIErVu3TuHh4WrTpo1OnDhRwpUCAAB4jpenC3BVhw4d1KFDhwKnGWM0fvx4vfDCC+rSpYskadq0aQoLC9Nnn32m/v37l2SpAAAAHlNmrtxdyu7du5Wenq62bds62ux2u1q0aKFVq1YVulx2drYyMzOdHgAAAGWZJcJdenq6JCksLMypPSwszDGtIKNGjVJQUJDjERUV5dY6AQAA3M0S4S6PzWZzem6Mydd2oWHDhun48eOOx969e91dIgAAgFuVmc/cXUp4eLikP6/gRUREONoPHjyY72rehex2u+x2u9vrAwAAKCmWuHIXFxen8PBwLVy40NF29uxZLVu2TE2bNvVgZQAAACWrzFy5O3nypHbt2uV4vnv3bm3cuFHBwcGKjo7WU089pZEjR6pmzZqqWbOmRo4cKT8/P3Xv3t2DVQMAAJSsMhPu1q9fr5YtWzqeDxkyRJLUq1cvTZ06Vc8995yysrI0cOBAHT16VAkJCfr+++8VGBjoqZIBAABKXJkJd4mJiTLGFDrdZrMpOTlZycnJJVcUAABAKWOJz9wBAADgT4Q7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZSZu6WBQDAnVJTUz1dQpGEhIQoOjra02WgFCLcAQCuaVnHD0uyqUePHp4upUh8ff20fXsqAQ/5EO4AANe0c6dPSDK6ofvzqhpX29PluCRz/x6tnTxCGRkZhDvkQ7gDAEBSQGi0gqPjPV0GcNW4oQIAAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCvDxdAFCc0tLSlJGR4ekyXJaamurpEgAAFkO4g2WkpaWpdu06yso67elSiuxc9llPlwAAsAjCHSwjIyNDWVmnldA3SRUjYj1djkv2b1mtrfM+VE5OjqdLAQBYBOEOllMxIlbB0fGeLsMlmfv3eLoEAIDFcEMFAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCLBPukpOTZbPZnB7h4eGeLgsAAKBEWep77urVq6dFixY5npcvX96D1QAAAJQ8S4U7Ly8vrtYBAIBrmqXC3c6dOxUZGSm73a6EhASNHDlS1113XaHzZ2dnKzs72/E8MzPTbbWVtR+0l/7cP3a73dNluCw1NdXTJQAA4HGWCXcJCQn6+OOPVatWLR04cECvvvqqmjZtqm3btqlKlSoFLjNq1CiNGDHC7bWV2R+0t9kkYzxdRZGdyz7r6RIAAPAYy4S7Dh06OP7doEED3XbbbapRo4amTZumIUOGFLjMsGHDnKZlZmYqKiqq2Gsryz9of0P351U1rrany3FJXs05OTmeLgUAAI+xTLi7mL+/vxo0aKCdO3cWOo/dbi/Rtx3L4g/aB4RGl7maAQC4llnmq1Aulp2drdTUVEVERHi6FAAAgBJjmXA3dOhQLVu2TLt379batWv1l7/8RZmZmerVq5enSwMAACgxlnlb9vfff1e3bt2UkZGhqlWr6tZbb9WaNWsUExPj6dIAAABKjGXC3cyZMz1dAgAAgMdZ5m1ZAAAAEO4AAAAshXAHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFeHm6AAAAcGVSU1M9XUKRhISEKDo62tNlWB7hDgCAMibr+GFJNvXo0cPTpRSJr6+ftm9PJeC5GeEOAIAy5tzpE5KMbuj+vKrG1fZ0OS7J3L9HayePUEZGBuHOzQh3AACUUQGh0QqOjvd0GShluKECAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhlgt3EydOVFxcnHx8fHTjjTdqxYoVni4JAACgxFgq3M2aNUtPPfWUXnjhBf3000+644471KFDB6WlpXm6NAAAgBJhqXA3btw4Pfroo3rsscdUp04djR8/XlFRUXr//fc9XRoAAECJsEy4O3v2rDZs2KC2bds6tbdt21arVq3yUFUAAAAly8vTBRSXjIwMnT9/XmFhYU7tYWFhSk9PL3CZ7OxsZWdnO54fP35ckpSZmVmstZ08eVKSdOS3HcrJzirWdbtL5v7fJEnH/9gpby+bh6txDTWXjLJYs1Q266bmkkHNJSMz/c+PSJ08ebLYX2fz1meMKdb1llnGIv744w8jyaxatcqp/dVXXzXx8fEFLpOUlGQk8eDBgwcPHjws8Ni7d29JRI5SzzJX7kJCQlS+fPl8V+kOHjyY72penmHDhmnIkCGO57m5uTpy5IiqVKkim+3S/xPKzMxUVFSU9u7dq4oVK159B0op+mkd10IfJfppNfTTWtzVT2OMTpw4ocjIyGJbZ1lmmXBXoUIF3XjjjVq4cKHuu+8+R/vChQt17733FriM3W6X3W53aqtUqVKRtluxYkVL/yHmoZ/WcS30UaKfVkM/rcUd/QwKCirW9ZVllgl3kjRkyBD17NlTN910k2677TZ9+OGHSktL0+OPP+7p0gAAAEqEpcLdgw8+qMOHD+vll1/W/v37Vb9+fc2fP18xMTGeLg0AAKBEWCrcSdLAgQM1cOBAt2/HbrcrKSkp39u6VkM/reNa6KNEP62GflrLtdJPT7MZw33DAAAAVmGZLzEGAAAA4Q4AAMBSCHcAAAAWQrgDAACwEMJdIY4ePaqePXsqKChIQUFB6tmzp44dO3bJZWbPnq127dopJCRENptNGzduzDdPenq6evbsqfDwcPn7+6tJkyb697//7Z5OuMBd/ZSk1atXq1WrVvL391elSpWUmJiorCzP/LauO/sp/fnt6B06dJDNZtPcuXOLtfaicEc/jxw5oieffFLx8fHy8/NTdHS0Bg8e7PgtZk9w1/HMzs7Wk08+qZCQEPn7++uee+7R77//7p5OuOBK+mmMUXJysiIjI+Xr66vExERt27bNaR4rjEOu9FMqPeOQO/uYN29ZHYMu18/SOAaVdoS7QnTv3l0bN25USkqKUlJStHHjRvXs2fOSy5w6dUrNmjXT66+/Xug8PXv21I4dOzRv3jxt2bJFXbp00YMPPqiffvqpuLvgEnf1c/Xq1Wrfvr3atm2rH3/8UevWrdOgQYNUrpxnTjl39TPP+PHjL/uTdSXBHf3ct2+f9u3bpzfffFNbtmzR1KlTlZKSokcffdQdXXCJu47nU089pTlz5mjmzJlauXKlTp48qY4dO+r8+fPF3QWXXEk/x4wZo3HjxmnChAlat26dwsPD1aZNG504ccIxjxXGIVf6WZrGIXf1MU9ZHoMu18/SOAaVeh77VdtS7OeffzaSzJo1axxtq1evNpLM9u3bL7v87t27jSTz008/5Zvm7+9vPv74Y6e24OBg89FHH1113UXlzn4mJCSYF198sTjLvWLu7KcxxmzcuNFUr17d7N+/30gyc+bMKabKi8bd/bzQ559/bipUqGDOnTt3NSVfEXf189ixY8bb29vMnDnT0fbHH3+YcuXKmZSUlGKr31VX0s/c3FwTHh5uXn/9dUfbmTNnTFBQkJk0aZKjrayPQ672s7SMQ+7sozFlewwqSj8v5MkxqCzgyl0BVq9eraCgICUkJDjabr31VgUFBWnVqlVXte7bb79ds2bN0pEjR5Sbm6uZM2cqOztbiYmJV1l10bmrnwcPHtTatWsVGhqqpk2bKiwsTC1atNDKlSuLo+wic+fxPH36tLp166YJEyYoPDz8aku9Ku7s58WOHz+uihUrysur5L8H3V393LBhg86dO6e2bds62iIjI1W/fv1i33+uuJJ+7t69W+np6U59sNvtatGihdMyZX0ccqWfpWkccuexLOtjkKv9vJgnx6CygHBXgPT0dIWGhuZrDw0NVXp6+lWte9asWcrJyVGVKlVkt9vVv39/zZkzRzVq1Liq9V4Jd/Xzf//7nyQpOTlZ/fr1U0pKipo0aaLWrVtr586dV7zeK+XO4/n000+radOmuvfee69qPcXBnf280OHDh/XKK6+of//+xbbOonBXP9PT01WhQgVVrlzZqT0sLKxY919R6ilqP/Paw8LCnNov7kNZH4dc6WdpGofceSzL+hjkaj8v5OkxqCy4psJdcnKybDbbJR/r16+XpAI/u2CMuerPNLz44os6evSoFi1apPXr12vIkCG6//77tWXLlqta74U83c/c3FxJUv/+/dWnTx81btxYb731luLj4zV58uQrXu/FPN3PefPm6YcfftD48eOveB2u8HQ/L5SZmam7775bdevWVVJSUrGsM09p6qc711sS/bx4+sXLWGUculQ/S2Ic8nQfrTQGXe6czePOMchKrqnrmYMGDdJDDz10yXliY2O1efNmHThwIN+0Q4cO5fvfRVH8+uuvmjBhgrZu3ap69epJkho1aqQVK1bovffe06RJk6543RfydD8jIiIkSXXr1nVqr1OnjtLS0q54vRfzdD9/+OEH/frrr6pUqZJTe9euXXXHHXdo6dKlV7zuC3m6n3lOnDih9u3bKyAgQHPmzJG3t/dVr/NCnu5neHi4zp49q6NHjzpdvTt48KCaNm16xeu9mDv7mfe2XHp6uuPvUPqzD3nLWGEccqWfJTEOebqPVhiDXOlnHnePQZbiiQ/6lXZ5Hwpdu3ato23NmjVX/YHtzZs3G0nm559/dmpv27at6devX7HUXhTu6mdubq6JjIzM90HmG264wQwbNqxYai8Kd/Vz//79ZsuWLU4PSebtt982//vf/4q7G5flrn4aY8zx48fNrbfealq0aGFOnTpVnGUXmbv6mXdDxaxZsxxt+/bt8/gNFUXpZ96H00ePHu1oy87OdvpwuhXGIVf6WZrGIXf10QpjkCv9NKZ0jUFlAeGuEO3btzcNGzY0q1evNqtXrzYNGjQwHTt2dJonPj7ezJ492/H88OHD5qeffjLffvutkWRmzpxpfvrpJ7N//35jjDFnz541119/vbnjjjvM2rVrza5du8ybb75pbDab+fbbb0u0f3nc0U9jjHnrrbdMxYoVzRdffGF27txpXnzxRePj42N27dpVYn27kLv6eTF58E41Y9zTz8zMTJOQkGAaNGhgdu3aZfbv3+945OTklGj/8rjreD7++OOmevXqZtGiRea///2vadWqlWnUqFGZ6ufrr79ugoKCzOzZs82WLVtMt27dTEREhMnMzDTGWGcculw/jSld45C7+nixsjgGXa6fpXEMKu0Id4U4fPiwefjhh01gYKAJDAw0Dz/8sDl69KjTPJLMlClTHM+nTJliJOV7JCUlOeb55ZdfTJcuXUxoaKjx8/MzDRs2zPeVBCXJXf00xphRo0aZ6tWrGz8/P3PbbbeZFStWuL9DhXBnPy9ehycHVnf0c8mSJQVOl2R2795dYn27kLuOZ1ZWlhk0aJAJDg42vr6+pmPHjiYtLa1kOlWAK+lnbm6uSUpKMuHh4cZut5vmzZubLVu2OC1jhXHIlX4aU3rGIXf28eJ1lLUx6HL9LI1jUGlnM8aYor6VCwAAgNLpmrpbFgAAwOoIdwAAABZCuAMAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAACyHcAYAHJScn64Ybbig16wFQ9hHugDLs4MGD6t+/v6Kjo2W32xUeHq527dpp9erVni6t1Pryyy+VmJiooKAgBQQEqGHDhnr55Zd15MgRT5fmMpvNprlz5zq1DR06VIsXL/ZMQQBKFcIdUIZ17dpVmzZt0rRp0/TLL79o3rx5SkxMLBVB5dy5c54uIZ8XXnhBDz74oG6++WZ999132rp1q8aOHatNmzbpk08+ueL1FtTXku5/QECAqlSpUqLbBFBKefr3zwBcmaNHjxpJZunSpZecT5KZOHGiad++vfHx8TGxsbHm888/d5rn999/Nw888ICpVKmSCQ4ONvfcc4/Tbzb++OOP5s477zRVqlQxFStWNM2bNzcbNmzIt53333/f3HPPPcbPz8+89NJLJikpyTRq1Mj861//MlFRUcbf3988/vjjJicnx4wePdqEhYWZqlWrmldffdVpXWPHjjX169c3fn5+pnr16mbAgAHmxIkTjulTpkwxQUFBJiUlxdSuXdv4+/ubdu3amX379hW6H9auXWskmfHjxxe6P/NMnDjRXHfddcbb29vUqlUr3++uXq6vcXFxxmazmdzcXHPs2DHTr18/U7VqVRMYGGhatmxpNm7c6FhX3nKu7uuYmBin39aMiYkpcD3nz583I0aMMNWqVTMVKlQwjRo1Mt99951j+u7du40k8+WXX5rExETj6+trGjZsaFatWlXoPgRQNnDlDiijAgICFBAQoLlz5yo7O/uS8w4fPtxxla9Hjx7q1q2bUlNTJUmnT59Wy5YtFRAQoOXLl2vlypUKCAhQ+/btdfbsWUnSiRMn1KtXL61YsUJr1qxRzZo1ddddd+nEiRNO20lKStK9996rLVu2qG/fvpKkX3/9Vd99951SUlI0Y8YMTZ48WXfffbd+//13LVu2TKNHj9aLL76oNWvWONZTrlw5vfPOO9q6daumTZumH374Qc8995zTtk6fPq0333xTn3zyiZYvX660tDQNHTq00H0wffp0BQQEaODAgQVOr1SpkiRpzpw5+tvf/qZnnnlGW7duVf/+/dWnTx8tWbLksn3dtWuXPv/8c3355ZfauHGjJOnuu+9Wenq65s+frw0bNqhJkyZq3bp1oVdXL7ev161bJ0maMmWK9u/f73h+sbfffltjx47Vm2++qc2bN6tdu3a65557tHPnTqf5XnjhBQ0dOlQbN25UrVq11K1bN+Xk5BS6HwGUAZ5OlwCu3L///W9TuXJl4+PjY5o2bWqGDRtmNm3a5DSPJPP44487tSUkJJgBAwYYY4z517/+ZeLj401ubq5jenZ2tvH19TULFiwocLs5OTkmMDDQfP31107beeqpp5zmS0pKMn5+fiYzM9PR1q5dOxMbG2vOnz/vaIuPjzejRo0qtJ+ff/65qVKliuP5lClTjCSza9cuR9t7771nwsLCCl1Hhw4dTMOGDQudnqdp06amX79+Tm3333+/ueuuuxzPC+urt7e3OXjwoKNt8eLFpmLFiubMmTNO89aoUcN88MEHjuUuvOJ2scL29Zw5c/Jt/8L1REZGmtdee81pnptvvtkMHDjQGPN/V+4++ugjx/Rt27YZSSY1NbXQegCUfly5A8qwrl27at++fZo3b57atWunpUuXqkmTJpo6darTfLfddlu+53lX7jZs2KBdu3YpMDDQcTUwODhYZ86c0a+//irpzxs3Hn/8cdWqVUtBQUEKCgrSyZMnlZaW5rTem266KV+NsbGxCgwMdDwPCwtT3bp1Va5cOae2gwcPOp4vWbJEbdq0UbVq1RQYGKhHHnlEhw8f1qlTpxzz+Pn5qUaNGo7nERERTuu4mDFGNput0Ol5UlNT1axZM6e2Zs2aOfbXpfoaExOjqlWrOp5v2LBBJ0+eVJUqVRz7NiAgQLt373bs24u5uq8vJTMzU/v27XOpHw0bNnT8OyIiwlEDgLLLy9MFALg6Pj4+atOmjdq0aaOXXnpJjz32mJKSktS7d+9LLpcXdHJzc3XjjTdq+vTp+ebJCyq9e/fWoUOHNH78eMXExMhut+u2225zvG2bx9/fP986vL298223oLbc3FxJ0m+//aa77rpLjz/+uF555RUFBwdr5cqVevTRR51uUihoHcaYQvtbq1YtrVy5UufOncu37MUuDoEFBcOC+npxW25uriIiIrR06dJ88+a9DXwxV/e1K1zpx4X74sJzAkDZxZU7wGLq1q3rdIVLktPn2fKe165dW5LUpEkT7dy5U6Ghobr++uudHkFBQZKkFStWaPDgwbrrrrtUr1492e12ZWRkuKX+9evXKycnR2PHjtWtt96qWrVqad++fVe93u7du+vkyZOaOHFigdOPHTsmSapTp45WrlzpNG3VqlWqU6dOkbfZpEkTpaeny8vLK9++DQkJKXAZV/a1t7e3zp8/X+h2K1asqMjIyGLrB4CyhXAHlFGHDx9Wq1at9Omnn2rz5s3avXu3vvjiC40ZM0b33nuv07xffPGFJk+erF9++UVJSUn68ccfNWjQIEnSww8/rJCQEN17771asWKFdu/erWXLlulvf/ubfv/9d0nS9ddfr08++USpqalau3atHn74Yfn6+rqlXzVq1FBOTo7effdd/e9//9Mnn3yiSZMmXfV6ExIS9Nxzz+mZZ57Rc889p9WrV+u3337T4sWLdf/992vatGmSpGeffVZTp07VpEmTtHPnTo0bN06zZ8++5M0ahbnzzjt12223qXPnzlqwYIH27NmjVatW6cUXX9T69esLXMaVfR0bG6vFixcrPT1dR48eLXA9zz77rEaPHq1Zs2Zpx44d+vvf/66NGzfqb3/7W5H7AaBsIdwBZVRAQIASEhL01ltvqXnz5qpfv76GDx+ufv36acKECU7zjhgxQjNnzlTDhg01bdo0TZ8+XXXr1pX052fXli9frujoaHXp0kV16tRR3759lZWVpYoVK0qSJk+erKNHj6px48bq2bOnBg8erNDQULf064YbbtC4ceM0evRo1a9fX9OnT9eoUaOKZd2jR4/WZ599prVr16pdu3aqV6+ehgwZooYNG6pXr16SpM6dO+vtt9/WG2+8oXr16umDDz7QlClTlJiYWOTt2Ww2zZ8/X82bN1ffvn1Vq1YtPfTQQ9qzZ4/CwsIKXMaVfT127FgtXLhQUVFRaty4cYHrGTx4sJ555hk988wzatCggVJSUjRv3jzVrFmzyP0AULbYzKU+pAKgzLPZbJozZ446d+7s6VIAACWAK3cAAAAWQrgDAACwEL4KBbA4PnkBANcWrtwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYyP8DQX5EsTV/CoYAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot train correlation distribution vs test correlation distribution\n",
    "import seaborn as sns\n",
    "\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": 95,
   "id": "cfa5cc5e-a7c8-44b2-92de-2c5d083721a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(\n",
    "    \"figures/reconstruction_spearmans/GSFA_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": [
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "from scipy.spatial.distance import pdist, squareform\n",
    "\n",
    "# Put into df\n",
    "pert_emb_df = pd.DataFrame(beta[:-1], index=pert_labels)\n",
    "pert_emb_df = pert_emb_df.drop(columns=[0]).astype(\"float\")\n",
    "\n",
    "# Compute the pairwise distances\n",
    "df = pert_emb_df.drop(index=[\"ctrl\", \"intergenic\"])\n",
    "\n",
    "distance = \"euclidean\"\n",
    "\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": 10,
   "id": "5c7d0199-28ac-4e7c-b4b5-180844815fe6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiQAAAHMCAYAAADlDLMCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABil0lEQVR4nO3dd1xV9eM/8PdhgwKCA0URBwoIpqioCCgqouTeppm5Z+6F2zR3qeXKrFzlykz9pmZLLdNyISooTlyomYgLLnB5/f7gd8/nXpb3wr0ckNfz8eBRXC/nvjmc8TrvKQGAICIiIlKQmdIFICIiImIgISIiIsUxkBAREZHiGEiIiIhIcQwkREREpDgGEiIiIlIcAwkREREpjoGEiIiIFGdR0B+Ynp4u7t+/L+zt7YUkSQX98URERJQHAMTz58+Fq6urMDMzfn1GgQeS+/fvCzc3t4L+WCIiIjKCO3fuiEqVKhl9uwUeSOzt7YUQGb+Qg4NDQX88ERER5cGzZ8+Em5ubfB83NoMDyb1798SUKVPEwYMHRVJSkqhZs6b48ssvRf369fX6eU0zjYODAwMJERFREWOq7hYGBZKEhAQRGBgomjdvLg4ePCjKlSsnrl+/LkqVKmWSwhEREVHxYFAgWbx4sXBzcxNff/21/FqVKlWMXSYiIiIqZgzqJrtv3z7RoEED0b17d1GuXDnh5+cnvvjii1x/RqVSiWfPnul8EREREWkzKJDcuHFDrF27VtSoUUP89NNPYtiwYWL06NFi8+bNOf7MwoULhaOjo/zFETZERESUmQQA+r7ZyspKNGjQQPz111/ya6NHjxanTp0SJ06cyPZnVCqVUKlU8veaXrqJiYns1EpERFREPHv2TDg6Oprs/m1QDUmFChVErVq1dF7z9vYWt2/fzvFnrK2t5RE1HFlDRERE2TEokAQGBoorV67ovBYbGyvc3d2NWigiIiIqXgwKJOPGjRMnT54UCxYsENeuXRPffvutWL9+vRg5cqSpykdERETFgEGBxN/fX+zZs0ds27ZN+Pr6innz5okVK1aIPn36mKp8REREVAwY1KnVGEzdKYaIiIiMr1B1aiUiIiIyBQYSIiIiUhwDCRERESnO4NV+izoAIilVrXQxqIDZWpqbbIVKIiLKv2IVSACIbutOiDNxCUoXhQpYA3cnsWtYAEMJEVEhVayabJJS1QwjxdTpuATWjBERFWLFqoZE2+kZocLOylzpYpCJvUpRiwbzf1G6GERE9BrFNpDYWZkLO6ti++sTEREVKsWqyYaIiIgKJwYSIiIiUhwDCRERESmOgYSIiIgUx0BCREREimMgISIiIsUxkBAREZHiGEiIiIhIcQwkREREpDgGEiIiIlIcAwkREREpjoGEiIiIFMdAQkRERIpjICEiIiLFMZAQERGR4hhIiIiISHEMJERERKQ4BhIiIiJSHAMJERERKY6BhIiIiBTHQEJERESKYyAhIiIixTGQEBERkeIYSIiIiEhxDCRERESkOAYSIiIiUhwDCRERESmOgYSIiIgUx0BCREREimMgISIiIsUxkBAREZHiGEiIiIhIcQwkREREpDgGEiIiIlIcAwkREREpjoGEiIiIFMdAQkRERIpjICEiIiLFMZAQERGR4hhIiIiISHEMJERERKQ4gwLJnDlzhCRJOl/ly5c3VdmIiIiomLAw9Ad8fHzEL7/8In9vbm5u1AIRERFR8WNwILGwsGCtCBERERmVwX1Irl69KlxdXUXVqlVFr169xI0bN0xRLiIiIipGDKohadSokdi8ebOoWbOmePjwoZg/f75o0qSJuHTpkihdunS2P6NSqYRKpZK/f/bsWf5KTERERG8cg2pIwsPDRdeuXUXt2rVFaGio+PHHH4UQQmzatCnHn1m4cKFwdHSUv9zc3PJXYiIiInrj5GvYb4kSJUTt2rXF1atXc3xPRESESExMlL/u3LmTn48kIiKiN5DBnVq1qVQqERMTI4KDg3N8j7W1tbC2ts7PxxAREdEbzqAakokTJ4qjR4+Kmzdvir///lt069ZNPHv2TPTr189U5SMiIqJiwKAakrt374p33nlHPH78WJQtW1Y0btxYnDx5Uri7u5uqfERERFQMGBRItm/fbqpyEBERUTHGtWyIiIhIcQwkREREpDgGEiIiIlIcAwkREREpjoGEiIiIFMdAQkRERIpjICEiIiLFMZAQERGR4hhIiIiISHEMJERERKQ4BhIiIiJSHAMJERERKY6BhIiIiBTHQEJERESKYyAhIiIixTGQEBERkeIYSIiIiEhxDCRERESkOAYSIiIiUhwDCRERESmOgYSIiIgUx0BCREREimMgISIiIsUxkBAREZHiGEiIiIhIcQwkREREpDgGEiIiIlIcAwkREREpjoGEiIiIFMdAQkRERIpjICEiIiLFMZAQERGR4hhIiIiISHEMJERERKQ4BhIiIiJSHAMJERERKY6BhIiIiBRnoXQBqHACIJJS1UoXI99epaRl+/9Fma2luZAkSeliEBEZFQMJZQFAdFt3QpyJS1C6KEbVYP6vShfBKBq4O4ldwwIYSojojcImG8oiKVX9xoWRN8npuIQ3ovaKiEgba0goV6dnhAo7K3Oli0FCiFcpatFg/i9KF4OIyCQYSChXdlbmws6KhwkREZkWm2yIiIhIcQwkREREpDgGEiIiIlIcAwkREREpjoGEiIiIFMdAQkRERIpjICEiIiLFMZAQERGR4vIVSBYuXCgkSRJjx441UnGIiIioOMpzIDl16pRYv369eOutt4xZHiIiIiqG8hRIXrx4Ifr06SO++OIL4eTkZOwyERERUTGTp0AycuRI0bZtWxEaGvra96pUKvHs2TOdLyIiIiJtBq+atn37dnH27Flx6tQpvd6/cOFCMXfuXIMLRkRERMWHQTUkd+7cEWPGjBFbt24VNjY2ev1MRESESExMlL/u3LmTp4ISERHRm8ugGpIzZ86IR48eifr168uvqdVqcezYMbFq1SqhUqmEubm5zs9YW1sLa2tr45SWiIiI3kgGBZKWLVuKCxcu6LzWv39/4eXlJaZMmZIljBARERHpw6BAYm9vL3x9fXVeK1GihChdunSW14mIiIj0xZlaiYiISHEGj7LJ7MiRI0YoBhERERVnrCEhIiIixTGQEBERkeIYSIiIiEhxDCRERESkOAYSIiIiUhwDCRERESmOgYSIiIgUx0BCREREimMgISIiIsUxkBAREZHiGEiIiIhIcQwkREREpDgGEiIiIlIcAwkREREpjoGEiIiIFMdAQkRERIpjICEiIiLFMZAQERGR4hhIiIiISHEMJERERKQ4BhIiIiJSHAMJERERKY6BhIiIiBTHQEJERESKYyAhIiIixTGQEBERkeIYSIiIiEhxDCRERESkOAYSIiIiUhwDCRERESnOQukCUPEEQCSlqpUuRpHyKiUt2/8n/dlamgtJkpQuBhFlg4GEChwA0W3dCXEmLkHpohRZDeb/qnQRiqQG7k5i17AAhhKiQohNNlTgklLVDCOkiNNxCayZIyqkWENCijo9I1TYWZkrXQx6w71KUYsG839RuhhElAsGElKUnZW5sLPiYUhEVNyxyYaIiIgUx0BCREREimMgISIiIsUxkBAREZHiGEiIiIhIcQwkREREpDgGEiIiIlIcAwkREREpjoGEiIiIFMdAQkRERIpjICEiIiLFMZAQERGR4hhIiIiISHEMJERERKQ4BhIiIiJSnEGBZO3ateKtt94SDg4OwsHBQQQEBIiDBw+aqmxERERUTBgUSCpVqiQWLVokTp8+LU6fPi1atGghOnbsKC5dumSq8hEREVExYGHIm9u3b6/z/UcffSTWrl0rTp48KXx8fIxaMCIiIio+DAok2tRqtdi1a5d4+fKlCAgIyPF9KpVKqFQq+ftnz57l9SOJiIjoDWVwp9YLFy6IkiVLCmtrazFs2DCxZ88eUatWrRzfv3DhQuHo6Ch/ubm55avARERE9OYxOJB4enqKyMhIcfLkSTF8+HDRr18/ER0dneP7IyIiRGJiovx1586dfBWYiIiI3jwGN9lYWVkJDw8PIYQQDRo0EKdOnRIrV64Un3/+ebbvt7a2FtbW1vkrJREREb3R8j0PCQCdPiJEREREhjKohmTatGkiPDxcuLm5iefPn4vt27eLI0eOiEOHDpmqfERERFQMGBRIHj58KPr27Svi4+OFo6OjeOutt8ShQ4dEq1atTFU+IiIiKgYMCiRffvmlqcpBRERExRjXsiEiIiLFMZAQERGR4hhIiIiISHEMJERERKQ4BhIiIiJSHAMJERERKY6BhIiIiBTHQEJERESKYyAhIiIixTGQEBERkeIYSIiIiEhxDCRERESkOAYSIiIiUhwDCRERESmOgYSIiIgUx0BCREREimMgISIiIsVZKF0AojcdAJGUqla6GMXaq5S0bP+flGFraS4kSVK6GFTIMJAQmRAA0W3dCXEmLkHpotD/12D+r0oXodhr4O4kdg0LYCghHWyyITKhpFQ1wwhRJqfjElhrSFmwhoSogJyeESrsrMyVLgaRYl6lqEWD+b8oXQwqpBhIiAqInZW5sLPiKUdElB022RAREZHiGEiIiIhIcQwkREREpDgGEiIiIlIcAwkREREpjoGEiIiIFMdAQkRERIpjICEiIiLFMZAQERGR4hhIiIiISHEMJERERKQ4BhIiIiJSHAMJERERKY6BhIiIiBTHQEJERESKYyAhIiIixTGQEBERkeIYSIiIiEhxDCRERESkOAYSIiIiUhwDCRERESmOgYSIiIgUZ6F0AXIDQCSlqo22vVcpadn+vzHYWpoLSZKMuk0iIqLiotAGEgCi27oT4kxcgkm232D+r8bdnruT2DUsgKGEiIgoDwptk01SqtpkYcQUTsclGLU2h4iIqDgptDUk2k7PCBV2VuZKFyNbr1LUosH8X5QuBhERUZFWJAKJnZW5sLMqEkUloiLO2H3X6H9M2Y+PMhTl/owG3eUXLlwovv/+e3H58mVha2srmjRpIhYvXiw8PT1NVT4iogJj6r5r9D/G7sdHGYpyf0aD+pAcPXpUjBw5Upw8eVL8/PPPIi0tTYSFhYmXL1+aqnxERAWmqPVdI8qsKPdnNKiG5NChQzrff/3116JcuXLizJkzomnTpkYtGBGRkgpz3zWizN6E/oz56piRmJgohBDC2dk5x/eoVCqhUqnk7589e5afjyQiKhDsu0ZUsPI87BeAGD9+vAgKChK+vr45vm/hwoXC0dFR/nJzc8vrRxIREdEbKs+BZNSoUSIqKkps27Yt1/dFRESIxMRE+evOnTt5/UgiIiJ6Q+WpPvKDDz4Q+/btE8eOHROVKlXK9b3W1tbC2to6T4UjIiKi4sGgQAJAfPDBB2LPnj3iyJEjomrVqqYqFxERERUjBgWSkSNHim+//Vbs3btX2NvbiwcPHgghhHB0dBS2trYmKSARERG9+QzqQ7J27VqRmJgoQkJCRIUKFeSvHTt2mKp8REREVAwY3GRDREREZGyFdrVfIiIiKj4YSIiIiEhxDCRERESkOAYSIiIiUhwDCRERESmOgYSIiIgUx0BCREREimMgISIiIsUxkBAREZHiGEiIiIhIcQwkREREpDgGEiIiIlIcAwkREREpzqDVfoneRABEUlqSSbb9KlWt9f9JQkjmRv8MWwtbIUmS0bdLRFSQGEioWAMg3jv4noj8N9I020+3FELME0IIEbKzmZDMUo3+GX7l/MSmNpsYSoioSGMgoWItKS3JZGFECCEks1Rh7z3VZNsXQohzj86JpLQkYWdpZ9LPISIyJQYSov/vSI8jwtbCVuli6C0pLUmE7AxRuhhEREbBQEL0/9la2LKWgYhIIRxlQ0RERIpjICEiIiLFMZAQERGR4hhIiIiISHEMJERERKQ4BhIiIiJSHAMJERERKY7zkBARERVCAERaSrpe701N+d+6WakqtUiFfktJWFiZFZplJxhIiIiIChkA4vulZ8WDG4l6vT9FQIhSGf//1aQ/hZXQL2RUqO4oOk+sVyhCCQMJERFRIZOWkq53GBFCCCshiUlPDV/6Iv56okhLSReW1sZfidxQDCRERESFWP8lQUYPDKkqtfh68p9G3WZ+MZAQEREVYpbW5oWiBsPUOMqGiIiIFMdAQkRERIpjICEiIiLFsQ8JESkCgEBSktLF0JGuNZdD+qskkZ5WuNrtJVvbQjE8k8gUGEiIqMABEHG9+4ikc+eULoqOZHMrIdovEEIIcTUwSNioUxQukS7bevWE+zdbGUrojcRAQkQFDklJhS6MCCGEjTpFHPxhotLFyFHS2bMCSUlCsrNTuihERsdAQkSKqnH8T2Fma/iETsVJelKSuBoYpHQxiEyKgYSIFGVmayvM+MRPVOxxlA0REREpjoGEiIiIFMdAQkRERIpjICEiIiLFMZAQERGR4hhIiIiISHEMJERERKQ4BhIiIiJSHAMJERERKY6BhIiIiBTHQEJERESKYyAhIiIixRkcSI4dOybat28vXF1dhSRJ4ocffjBBsYiIiKg4MTiQvHz5UtSpU0esWrXKFOUhIiKiYsjC0B8IDw8X4eHhpigLERERFVMGBxJDqVQqoVKp5O+fPXtm6o8kIiKiIsbknVoXLlwoHB0d5S83NzdTfyQREREVMSYPJBERESIxMVH+unPnjqk/koiIiIoYkzfZWFtbC2tra1N/DBERERVhJg8kRFR8ABBISnrt+9K13pOux/s1JFtbIUlSnspmbPr+rsaQ1/1lDIVpn9ObzeBA8uLFC3Ht2jX5+5s3b4rIyEjh7OwsKleubNTCEVHRAUDE9e4jks6dM+jnrgYG6f1e23r1hPs3WxW/Qeb1dzUGQ/aXMRSWfU5vPoMDyenTp0Xz5s3l78ePHy+EEKJfv35i48aNRisYERUtSEoy+Q066exZgaQkIdnZmfRzXqcgftfCorDsc3rzGRxIQkJCBABTlIWo2AMgktL0q5LXfp++PyOEELYWpq+Cr3H8T2Fma2u07aUnJRV4zYC+jP27FhaFeZ/Tm4l9SIgKCQDivYPvich/Iw3+2ZCdIXq/16+cn9jUZpNJQ4mZra0wKyZP1MXpdyUyJS6uR1RIJKUl5SmMGOrco3MG1agQERUE1pAQFUJHehwRthbGbQZISksyqCaFiKggMZAQFUK2FrbCzpLNAERUfLDJhoiIiBTHQEJERESKK95NNoAQqa/yt40Utdb/vxJCmOd9W5Z2QnDyISIiKoaKbyABhPiqtRB3/s7ndqyFEF9n/P9SDyEkVd635dZYiAGHGEqIiKjYKb6BJPVV/sOIEMJOUolbNr2NUCAhxJ2TGeWyKmGc7RERERURxTeQaJt4TQgrBUc0pLwSYpmHcp9PREWaKRb6M/WCfly0jzJjIBEiI4ywVuKN8aZMv06kj4JY6M8UU8hz0T7KjIGE3ihv0vTrRPooqgv9cdE+yoyBhN4oBT39Oicvo8KkKCz0x0X7KCcMJKak77DilFfZ///rcJhwrjj9OhU3XOiPijIGElPJ67BiQzq3cphwrjj9OhFR0cGZWk3FSMOKc6UZJkxERFTEsYakIBh7WDGHCRMR0RuGgaQgcFgxERFRrhhIiIiEYZOL5WXSME4ERpQ7BhIiKvbyM7mYvkNYOREYUe7YqZWIir2CmFxMMxEYEWWPNSRERFqMPbnYmzARmDHXyjHlGjlsFivaGEiIiLRwcjFdplwrx9hBjc1iRRubbIiIKEdFaa0cNosVbawhISIivRTWtXLehGYxYiAhIiI9maI5y5j9U4QwTr8U9kVRBgMJEREpwhT9U4xRU8K+KMpgHxIiIlJEYe2fwr4oymANCRERKa4w9E9hXxRlMZAQEZHiONya2GRDREREimMNCRERvbFMvWiiEByVYywMJERE9EYqiEUTheCoHGNhkw0REb2RCmoUD0flGAdrSIiI6I1nilE8HJVjXAwkRET0xuMonsKPTTZERESkOAYSIiIiUhwDCRERESmOgYSIiIgUx0BCREREimMgISIiIsUxkBAREZHiGEiIiIhIcQwkREREpDgGEiIiIlIcAwkREREpjoGEiIiIFMdAQkRERIrLUyBZs2aNqFq1qrCxsRH169cXf/zxh7HLRURERMWIwYFkx44dYuzYsWL69Oni3LlzIjg4WISHh4vbt2+bonxERERUDBgcSD755BMxcOBAMWjQIOHt7S1WrFgh3NzcxNq1a01RPiIiIioGDAokKSkp4syZMyIsLEzn9bCwMPHXX38ZtWBERERUfFgY8ubHjx8LtVotXFxcdF53cXERDx48yPZnVCqVUKlU8veJiYlCCCGePXuW62e9SkkT6apX8nvTrAwq6uulvBRCBfH/P0AIK3XR2b6Jy27qfW/K7b9KfSXUSer/bdsyzWjbNvX2i3LZhRAi/dUr8UL9v+2bpRlv+6bctqm3X5TLburtF+Wym3r7qSq1SEp5KW/b0trcaNvO6/Y1920ARi2LDAa4d+8ehBD466+/dF6fP38+PD09s/2Z2bNnQwjBL37xi1/84he/3oCvO3fuGBId9GbQo2mZMmWEubl5ltqQR48eZak10YiIiBDjx4+Xv09PTxdPnjwRpUuXFpIkGfLxREREpBAA4vnz58LV1dUk2zcokFhZWYn69euLn3/+WXTu3Fl+/eeffxYdO3bM9mesra2FtbW1zmulSpUyvKRERESkKEdHR5Nt2+DG+/Hjx4u+ffuKBg0aiICAALF+/Xpx+/ZtMWzYMFOUj4iIiIoBgwNJz549xX///Sc+/PBDER8fL3x9fcWBAweEu7u7KcpHRERExYAEmKq7LBEREZF+uJYNERERKY6BhIiIiBTHQEJERESKYyAhIiIixTGQEBFRoZdm5GnfqfBhIHmDXb16VRw9etSkn6EZpAXAdOsbFFHcHwVHrTbyWlQKSU9PV7oIhdLVq1fF4sWLxa1btxT5/PyeyzExMeKHH35gqHoNBhIFmfKG9ezZM9GmTRtx7do1k32GEEKe/j8tLc1oSwEUxYtyXFycWL16tViwYIH45ptvhBCiyC6NUJSClGYZC3Nz80IZSgw5lmNjY8W0adNEhw4dxKJFi0RkZKRJypSQkCASEhJMsm1TuHDhgmjcuLF4+vSpSE5OFkIUzDXi5s2bYv369QJAvs7l8+fPCx8fHxEbGyssLIy8SGwmRenczY5p944C0tPThZmZcXNWSkqKsLKyMsq27t69Kx4+fCjq169v0huWg4ODEEKI8uXLm+wz4uLixPbt28Xhw4fF3bt3RXBwsOjSpYt4++23hRDC4BP53r17wtnZWdja2prk76hx//59cebMGaFWq4Wnp6fw9vbO1/YuXLgg2rRpI3x9fUV8fLxITEwUMTExYv78+UYqsS5T7ptbt26J/fv3iwcPHojhw4eLSpUq5Xub165dEzt37hRXrlwRQUFBolWrVqJKlSr53u7ly5fFW2+9JUJCQsThw4flUGJubtxVUQ1x8eJFsWLFCtGnTx/h7e2tc/7ldj5cvHhRNG/eXLRp00ZIkiT27dsnLly4ID799FNRunRpo5Xv33//FZ07dxahoaFi9OjRwtnZOd/bfPLkiYiPjxfm5ubCxcVFODk5GaGkGR49eiR69uwp3n//fbF06VL59cTERKN+TmZRUVEiODhYVK5cWQwaNEhIkpSnYBIVFSUCAgLE5MmTxeTJk41ezps3b4rdu3eLxMRE0aJFC9G8efN8Byh9meRzTLJknwLi4uJw/fp1ndfUanW+txsdHY3+/fvjzJkz+d7W5cuX4eLiAn9/f/zxxx/53l5OUlNT8ezZM7i6uprsc6KiouDl5YXu3bujf//+GDNmDDw8PGBvb481a9YYvL2YmBhIkoQ6derg1atXAIzz98vs/Pnz8PDwgJ+fH9zd3REYGIiYmJg8b+/WrVuoVq0apkyZgvT0dDx8+BDLly9HUFAQ4uPjjVhy4OHDh0hMTARgmn1z4cIF+Pj4YNiwYZgyZYr8elpaGgAgPT09T9ssW7YsevTogYCAAPj7+2PYsGF48eJFvsp6//59BAYGokmTJqhRowbatm2bpbwFLTU1FU2bNoUkSRg3bhy8vb2xY8cOxMbG6rwv8368e/cu3nrrLZ19vmvXLlSoUAEXL140ejmHDBmCOnXqYNGiRXj8+HG+tnXhwgXUrVsXtWrVQpkyZdCvXz9cvnzZSCUFjh8/jgYNGiApKQlqtRqDBw9G8+bN4enpiQULFuDp06dG+yyNyMhI2NnZoWfPnihdujTWrVuXp+1cuXIFkiRh8uTJADL+7hs2bMDYsWPx0Ucf4bfffst3OStWrIjg4GDUq1cPdnZ2+d6mPv777z/5//NyTcjNGxFI7t69C0mSUKpUKXz00UfYvn27zr/n9eJ9/fp1uLu7Q5IktGjRIl8Xh/j4eISEhCAwMBDh4eEICwvDsWPH8ry97MTFxSEpKQkAcPv2bTg6OuLvv/826mcAwLlz51CyZElMnjxZ5+A8efIkevXqBTMzM2zcuFHv7T18+BAhISFo06YNvLy84O/vb5JQEhMTg3LlymHq1Kl4/PgxDh06hJo1a+Lo0aPyeww5wdRqNZYsWYLw8HCdC+OpU6dgb2+PCxcuGK3s0dHRqF69Ovr27YsnT57In2/M7ZcqVQqTJk3S+Zvu3LkTX3/9NV6+fGnwNm/fvo1atWph6tSp8murV69GtWrVcO/evXyVd+fOnejSpQuOHDmCH374AR4eHjqhJDU1NV/bz6s9e/agUaNGOHLkCFasWIHatWsjNDQUkyZNwt27d7MNdzt37kSnTp1w5coVnb+pn58fvvvuuyzvzyuVSiX//7hx4+Dj45OvUHLlyhW4uLhg4sSJuHTpEtavX4/GjRtj/fr1Rivz/v374e/vDwBo2bIlwsPDsWTJEixYsACWlpYYNGhQvj9DW2RkJGxtbREREQEAeOedd9C2bVs8f/7c4N/nu+++gyRJWLNmDVJSUhAcHCyH8sDAQNjb22Pz5s35Kue0adOQnJyMW7duoVmzZvj888/ztD19Xbp0CRYWFhgzZoz8mjFDyRsRSF68eIHOnTtj2LBhmDRpEqpXr4727dtj27ZtWRK0vhfx5ORkREREoGfPnvjnn39QqVIlBAYG6oQSQ/4Qp06dQsuWLXH8+HEcPHjQ6KFEpVKhcePGqF69OpKTk6FSqVC2bFmcO3fOKNvXiI2NhYWFBRYsWAAg69NoVFQUwsPD4ePjgytXrui1zZ9//hnvvvsufvvtN0RGRsLb21snlBjjiff58+cIDw/HsGHDdF5v06YNVqxYgc2bN+Ovv/4CoP/fNT09Hfv378enn34qv6ZWq/HixQu4u7tnGwbzcvLeuXMHjRo1gpeXF4KCgjB06FCjhpJnz55lu28++ugjSJKEmjVrYvPmzXIo0ed3UKvV2LRpEzp16oRbt27J5VSpVKhevTp+/vnnfJU5OTkZ+/fvlz9rz549WUKJEjUlV69eRefOneXf79q1azh69CgkSUJwcDB69eqF69ev48mTJ3j8+DGSk5Nx8uRJnRtJeno6UlJSUKtWrTw/nWvTBI7M+2P06NGoVatWnkLJixcv8M4772DAgAE6rw8cOBCNGzfO13H5+PFj+dw/c+YMLC0tsXz5cnTt2hW3b9+W33f06FGYm5tj165def4sbbGxsZAkCdOmTZNf27ZtG8zNzeVzWZ9j//Hjx/K58tVXX8HMzAwuLi7o1q0bbt68CSDjIXrixIlwdXXFpUuXDCrn1atX4ejoiOHDh+u83qpVK3Tu3BnNmzfH9OnTcfbsWYO2+zr37t1Dw4YNUa9ePZQoUQJjx46V/81YoaTIB5L09HQkJSVh2LBhmDlzJoCMavQhQ4agU6dO8PHxwXfffYfz58/rvU3NybR3715s27YNQEY1VaVKlRAUFISLFy9m+wd43R9Fu9nnxx9/lEOJ9hO65rPzckKfPn0avr6+8Pf3x5UrV1C3bl3s3bsXT58+xa1bt3Dt2jU8ePAA8fHxOHbsGP7991+Dtp+amorp06fD1tYWmzZtylJmjZ07d8Lc3BxHjhzRa7spKSn4/fff5e/PnTsnhxLNia1Wq5Genp7nC11ycjK+++47nZAwf/58SJKEJk2aICgoCJIkyU+k+tIEA0D37+/p6anzdz148GCeyg0AW7ZsQUhICE6cOIElS5YgICBAJ5Tk98b74MED1KhRA7t27ZJ/hx9//BGSJOGXX35Bnz594Ovri02bNhnU1HLw4EGsXr1a/j49PR0vXrxApUqVjHYT0UhOTsYPP/yA6tWr64SSNWvW5KtJLi9GjhwJLy8v+fv+/fujYsWKWLp0Kdq0aQNJktCtWzfUqVMHp06d0vlZ7WMoODgYX375pfz9559/juPHjxtUFk2zQWBgIAYMGIDff/9dp0llxowZ8PT0xKJFi/Do0aMsZcjJgwcPMGzYMOzevRvA/2qkduzYgSZNmuT5PD179iz8/Pxw8uRJ+bWRI0eiatWqcHd3l495zec1bdoUs2fPztNnaUtPT8eDBw+wdu1a+XuN0NBQdO7cWa9awuzKv3HjRvj7++u8BmTUKDs6OuKnn34yqKzr1q1D6dKlMWPGDLmmccGCBbC2tsawYcMwfPhwWFpaokOHDnj27JlB285Jeno6tm7diu7du+P48eP49ttvYW1tjXHjxum8J7+KdCBJT0+Xd8LZs2fh5OSEffv2yf/eokULWFtbo3HjxqhVqxaGDx/+2qf2ixcv4pNPPgGALAfgo0ePstSUpKeny0/WmSUkJODatWuIjY3NNmgcOHAAbdq0QevWreWakjFjxmQ5cHNz5coV+UkRyKjK8/T0hLOzM8zMzODs7Ixy5crB3t4e1tbWKFu2LJycnODs7JynavOYmBhMmDABnp6eOn1F1Gq1/LulpqbCwcFB52Kak5wO4uxqSj7++GODT16VSoWUlJQsVfg//vgjypcvj7179yIpKQmvXr3C8OHD4efnh//++y/Hct26dQvffPMNli5dKp/s2u9Vq9V4/vw5KlasiD///BMAMHPmTEiShDt37hh00iYkJMj//3//93/y/y9evFgOJZrmFe3jytCA8ueff0KSJERHR+tsQzvE9+rVCxUqVMjxWNcu89WrVxEXF6fzuvbv3aBBA+zdu1f+ftOmTa89LzX7fdmyZfJ+z/x7JiUlyc037dq1w+jRoyFJUpa+Zcakvd81v+O///6L0NBQHD58GO+++y7Kly+vU7M6e/ZsWFtbY9KkSTrbynxshIWFyU2fERERsLOzMzhc/fDDD5AkCQ4ODujUqRPKlCkDDw8PdOnSBRs3bsSTJ08wbNgwNG3aFEuWLHntQ4pKpZLPpcjIyCz74fDhw2jYsCFSU1Pl30ffvlSRkZGwsbHJsl9++uknuW+O9nkAZOyj5cuX67X9nFy4cEHu55GdxYsXo1q1arhx4waAnB8Wcyo/kNF/UNOcrvn56Oho+Pr64vTp03qV8+bNm3IN0UcffQQ/Pz/MnTsXU6dORdmyZXHo0CH5vbt374YkSXo/FOrj9u3bOuftN998A2tra6PWlBTJQHLr1i2cOHECQMYO0JwgY8eOxaxZswAA77//PlxdXXH9+nVERUVhzZo1qFKlCm7dupXjdtVqNdq3b49Ro0Zl+TfNZzx8+FCuKTl37hxGjBgBf39/nXZ3ICPYBAUFwc3NDZUrV8bQoUPlbWjfHDXNN23atEGnTp0gSZJBVW0zZsyAJEn4/vvv5f1x9uxZtGzZElZWVjh06BBOnz6NY8eO4fjx4zhx4gROnjypU/VpqOjoaIwdOxY1a9aUnyiA/51ov/zyC2rVqpVrVaQ+T1CaUNK4cWMMGTIEkiRl6SCYm6tXr2LkyJHo168fvv766yyfn7l8M2bMQFBQUI7bO3/+PKpUqYIGDRqgdOnSqFatmtzJVCM1NRVPnz6Fq6srzp8/j8WLF6NEiRJ6X3Q0rly5goCAgCw3dk3Zs6sp+eSTT/R+Mr116xauXr0KIKOZzcbGBp999lmWm7zmWD1x4gTq16+fa78Y7WO+UqVKOse89oWqYcOG+OGHHwAA06ZNg729fa5/V332u2b7KpVKvhg7OzsbpTN6du7fvy/XFmUXjN555x3Y2dmhZs2acrNpeno6IiMjUaJECbmPQk7UajXq1auHr776CgsXLoSNjY3Bx5DG9u3bYW5ujhUrViAqKgo//vgj2rdvj9q1a6NSpUpo27YtJEmCo6MjVq5cmeMxpDmf+vbtq9P3IXN/mKpVq8o33zlz5qBFixZ49epVrjer1+2X//u//0OTJk1gbm6Ozz77DLt378aUKVNQvnx5XLt2LS+7Rfb++++jX79+WV7XlDcpKQnu7u5ZmkgMKX92Jk+ejLp168o1U7l5/PgxypUrhx07dsivzZs3DzVr1oStrS2++eYbABnHYnp6Ok6cOAEPDw/8888/epfHUGlpaVlqSlJTU7F169Y8958rcoEkKSkJ/fv3h5ubmzyCRHPgbNq0CVWqVEHr1q1RqVKlLH8M7U5dOWnatKlOG6K2lJQUABk1JVWqVEGJEiVgY2OT5aIXGRmJkiVLYsSIEdixYwcGDhyISpUqYfLkyXJZtU/O/fv3w8nJCaVKldJ56tCHSqXCxIkTYWlpKVefqtVqnD17Fj4+PmjWrJl8cciL58+f49GjRzh8+DDOnTsnlzs2NhZjxozJEkqAjA5zbdu21WnO0NC+kehz8zx16pRc02NIUIuKioKrqytGjx6NLVu26Pztc6pBGDlyJAYPHgyVSpXl4qnpRDZ9+nQ8evQIN27cQJUqVfDtt99m2U56ejrq16+PgIAAWFtbZ6mW18fGjRvh6+sLQHc/af4/LS1NJ5QMGjRI78Cmubhpjhcgo1raw8MjxyfwiRMnonXr1tn+TQH9jnkg4/ytVq0a9uzZg0WLFr32RmvofgeAwYMHo0SJEga3zesrJiYGZcuWRWhoqHw8Z+6sGh0dDVdXV3z88cfyz8XHx8PS0hKDBw/Wee+MGTN0njI1WrVqhUqVKsHOzi5Px5C2L774ApIkYd68efJraWlp2Lp1K5YvX46aNWuiSpUqOdZU5XY+adu/fz9q1aoFAJg1axYsLS1fGwpz2y8ffPCB/L7IyEhMnToVZcqUwVtvvYWGDRsapY9chw4dMGTIkGz/TXO+zZ8/H35+fnItib7lz+7vGhUVhVGjRsHJyUnv6/29e/fg4uKCixcv6lwPli1bBh8fH0RERODOnTvy69OmTUPt2rXx8OFDvbafHc39Ljepqak6oWT06NGwtLTM9kFKH0UukADAr7/+il69eqFOnTpZOoW2bdsWzs7O2Z7A+lQntWzZUqfdOzPNE9+AAQNQunTpLCNvLl++jBIlSmD69OnyaykpKQgNDUXjxo11boaaPhFjx441eFSG9kGpUqkwbtw4WFpa6vTKj4yMhJeXFxo1amRQh0SNK1euoEePHqhVqxYkSULJkiXRsGFD+aZ17do1jBkzBp6ennIomT17NpydnbMdkRQfH482bdpk6QSa2+84ZMgQ2NraGnRzuXnzJipXrvzaKnGNFy9eYMaMGShbtqxOs4XG1atX5ZuituDgYEyaNAm9e/fGxo0b5YtVQkICSpUqBUtLS4P6LmlbunQp6tatm+2/afZZeno6Fi9eDFtbW5QqVUrvi3N2F7d//vkH1apVQ+3atfHPP//IF6MHDx5g0qRJcHBwQFRUVLbbM+SYV6vVCAoKgo+Pz2tvtIbudyCjxtHLyyvPtQmvc+/ePQQGBqJBgwZo0qQJ2rdvnyWUqNVqvHr1CgMHDsTAgQORkpKC9PR03Lx5E82bN0eVKlXkMi9cuBBOTk46za4amiaWnPZ7dm7duoVVq1bho48+wpYtW3T+bcOGDTAzM8OcOXOyNGH+999/OYZNQ86nI0eOoGXLlpgwYQKsra31+jsYul8ePnyIZ8+e5WvIr3bZO3bs+Np+KBcvXoQkSfjqq6/yVf7Y2FhMnToV/v7+Bj18XrlyBRUrVpSb2bXPKU1YmjRpEhISErBgwQLY2toa/HCr7eLFi+jQoYNe1920tDRs2bIFkiTByckpX+G5SAYSIOPA7969u04oUavVWLt2LRo1aiQ3zbzuKTwuLg6ff/45nj9/DgByNSmge9Bq/79m9EHmJ3a1Wo0BAwbA2dkZO3fu1Pm5efPmoVGjRnj69KnOtqKiolCxYkW9L6DaNQzaB6UmlFhZWcknQVpaGiIjI+Hi4oIWLVrotX2N8+fPo3z58hg1ahR27dqFuLg4fPrpp6hVqxbc3Nzk/hExMTEYO3YsateujeDgYNja2ub4RBQTE4POnTsjODgYX3zxhfx6Tn+jM2fOoHHjxgbfXD755BO0atUK8fHxrw1gv/zyCwYPHowKFSpkWwOjVqsxa9YslCtXDkuWLJFfX7hwIczNzdGnTx8EBwfD0tISI0aMQGJiItLS0rB27Vq9RxlpaN8kVqxYgaZNmwLIvkZHs89GjBgBBwcHg4akZ764ARkB4sCBA/Dx8YG9vT2aNm2K1q1bIzg4GO7u7jnWThl6zKelpSEwMBBOTk65hjVD97vm/FWpVHjw4IHe+8JQO3bsQKtWrXDs2DFs2bIFAQEBOqFE+2+4Y8cOSJKkExTj4uIQGhqKypUrY8KECShXrpxO27+2H374waD+L5pajLCwMNSuXRuVK1fOEuY0oWTevHl6zy9jyPmkGepaokQJg5rL9N0vxug4eevWLSxZsgT3798HkNHXUFOTpV2Drfl/zbk2Y8aMHG/Qhvxdr1y5olczzaNHj+T72OXLl+Hg4CCP0tEuF5ARSho1aoQ6derkq3kPyAhY1atXhyRJqFu37muvY2q1GgMHDoSDg0O2D3SGKBKBRBMaNmzYIPcdAYA//vhDDiWaEQ0JCQmoXLkyRo4cqde2R40aBU9PT6xatQppaWmoW7dutilY219//ZWlaltzgN29exddu3ZF06ZNsWHDBgAZib5EiRJYvHhxttvL3B6ek+xqGLRvVikpKRg1ahRsbGx0LvYXLlwwqJ1V0x46depUnRMzJSUF58+fR8OGDVG1alX5CeXKlSsYNGhQjjeu+/fvy++9dOkS+vXrh4CAAJ1Qov17pKamyhdx7Y6d+urQoQPCw8Oz/TfNSfzixQukpaXh1KlTWLRoUa4n3aNHjzBx4kQ0atQIK1aswKJFi1C2bFkcPHhQ3jeTJk2Cra2t3GRi6DwYd+7cQb169eRyfPTRR+jQoUOu29q9ezesrKz0uvDndnHTvsj/+++/mDJlCnr27ImOHTti+fLl2VZTa7YJGH7Mf/3113L/ldeV2dD9biqawANkdIjW0A4lmmNVu6q7c+fOuHLlis7xHRcXh44dO0KSJLnTanadYw2R2wR99+/f19nmhg0bYG1tjalTp+rVbKrv+aRSqXD+/Hk0bdpU7863huwXY5k9ezaqVauGDz/8EK9evUJISIg8kCE32Z2Hpvq7JiYmwtfXF0OHDsW9e/dw48YNODk5yc0y2e2XuXPnokaNGnmulQUyRqvNmTMHXbp0walTp9CwYUN4e3vnen08cOAAqlatmu9mRaAIBJLz58+jYsWKCAoKQunSpREeHq7TN+To0aPyEDrNLHXz5s1D7dq1kZiY+NqD4Pnz5xg4cCD8/f2xdOlSeVz+d999h127duH777/Hzp07sX37dqxZswa//PJLlm3cvn0bFSpUkIeU3r17F506dULLli2xdOlSVKpUSaejbOZ+JPoeqDnVMGifFI8ePUKrVq0watSoPE0Olbk9NLvRQUePHoWTk5PO5DjXr1/P9sn06dOnaNOmDdq1aydXCWuHEs0kSprPUKlUGDlyJBo2bGhwGElPT4dKpUKnTp3QvXt3ADn3G5o5c6bcvJVbW6lm3/77778YN24cvLy8YG5ujsOHDwOAPALo8OHDqFatWp4nzzt79iwaNmyIGjVq4NatW/joo4/Qrl27HN+v+XvoM1LqdRe3nPrU5HZc5uWYN+SCbOh+N1V/ESAjoIeHh2cbejTV1ZpQogneH3/8Mc6ePZvjRfratWto27YtKlasKNeC5PXmm5aWhsWLFxs0Qd+qVatQunTpXEfVGHo+aYZy59T0oxEfH5/jTdOY+yU3U6ZMQd26dfHhhx/Cx8cHo0ePxvr167Fu3TqsX78ea9aswZo1azB//ny5w2hBl3/VqlVwd3fHuHHj8PXXX8PHxwcPHjzA06dP8d9//yEhIQEJCQn4999/5YeG/DRjaaZW2L17t/y3TEhIeG0ouXfvntFmpS7UgUTToW3q1KlQqVQ4efIkypUrJ/fQ1zhy5Ai6deuG+vXr488//8TDhw9z7VQTHx+Pw4cPy086L1++xHvvvYc6depAkiRYWFjA29sbZcqUgZOTE1xdXVG5cmW4urpmm/yPHDmCypUrQ6VSyQef5gLt6OiIkJAQ+WKcl5Bw9+7dbG/mOYWS7t27o1evXgZ/DqDbHprTCZWeni6PzX+dtLQ0LFu2DMHBwejdu/drQ8nIkSP1furPyfTp02FpaSnfQDLfcO/fv49u3bplGy6BrBddzc8/fvwYEyZMwFtvvSVPDKcxbty4bEdbGSIqKgrNmzdH9erV0bVrVzRu3BhdunRB3759MWTIEPTu3Rs9evRA9+7dMWLECINu8Ppe3B4/foxLly69tiO0KY55pfZ7biIjI2FhYYEZM2Zk+TftYe5bt25FQEAAOnTogMGDB0MIAXt7ezg7O2P8+PH4448/svy9bt68idDQUHk0oPY29aXZZ/v27cNnn32mU7bcJugD9L956Xs+acJibu7fvw9HR0eT7xd9jBs3Dm+99RYkSYKVlRUaNWoENzc3uLq6ombNmvDy8soSdgu6/F988QUqV66M5s2by300ypYtC1dXVzg7O6NUqVJwcXGBu7t7vjqwXrlyRV66IPN5+PjxYzmUaNcAHz9+PF8DJrJTaANJTEwMLCws5MnONFq3bo2hQ4di2LBhmD9/vvz60aNHERYWhqCgoFxH01y8eBH16tVDu3bt0LlzZ/lJ6+XLlxg6dCjq1KmDiRMnyifro0ePkJycDLVanWPTyv79+1GlShWdTm1ARvDp0qULmjZtqjORmCE3Es0smk2aNJFnU7x48SLef//9LDdzzZP+gAEDMH36dDnxGur27duvPaE6deqEjh075rod7REhq1evRkBAQI6hZN26dRg/fjxsbW0NGk1z+/ZtbNmyBWvWrJFDTGxsLGrUqIGaNWtmG0xnz54Nf39/uQ1ZW0xMDOrVq4dFixZle4HVPLE3bNgQc+fOBZDRfluiRAmDq0rj4uKwYcMGLFu2TO4HdenSJXTo0AGSJMHHxwejRo1CmzZt0LFjR3Tp0gXdunVDz54981Qtq+/FrXr16q/th2HsY74g97u+zp8/n22n2szzzgAZx/imTZtga2sLJycn7N69G61atcLmzZvRrVs3hIWFoUmTJvjnn390htzfvn0brVq1grW1tU7/AH1ERUVh6NChePr0qc4EWLlN0HfgwIFct2ns8ymzmJgYk++XzGJjYzFr1iz06NEDGzdu1BmNMnPmTLz11lv44IMP5FCrvS8zz0VlyvLHx8fjjz/+yDKL8caNG1GqVCkEBATI58eRI0ewb98+7N+/H8eOHcv30Oddu3ZBkiSdCRwzN+VqQsnFixcxbNgwNGzYMN9rIWVWaAPJ6tWrIUmSzvC+BQsWQJIk9OvXT57xsHfv3vK/nzhxQudgy+zSpUtwcnLC9OnTdZ6oNBfVly9f4v3330ejRo2wcuVKJCcnA9Ad1aDx7NkzOR1+9913qFChAl68eJGliUO7Kjsv00Cr1Wp89dVXaNq0KcLDw3VCSb9+/dC4cWO509/Dhw8xZ84clCtXLk8LXGVuD80cSjT/fu/ePYSFhckTo2W+2Wh/r/kZ7VDyzjvvyM0x0dHRGDhwIJycnAyuGTl//jzc3d0REBAASZJQq1YtuXp1w4YNqFKlCtzc3LB7925cuHABhw4dwtChQ+Ho6JilB7qmzBEREZAkCcuXL0epUqUwZsyYLGsjPXjwAGPHjkVQUBDq16+v92iCzGV3c3ND/fr1YWFhAXd3d7m/xenTp9GjRw+4ubnlqQ8NYJqLmymO+YLe7/q6du0aJEnSWfAOAJYvX55l5ITmdxg+fDhKlCghN9u1b98eQ4cOBZBxIxs4cCBCQkIQEBCAbdu2yU+08fHx6Nixo0H9YCIjI2FmZiaHs8xlycsEfcY8n3Jjyv2SWVRUFFxcXNC5c2eEh4fD0dFRJygDwPjx41GvXj3MmTNHvr5md803ZfmjoqJQs2ZN1K5dG5IkZVmjZ9OmTahYsSKmTJmS61xaeZGeno6XL1+id+/eCAwMzLEj9ePHjxEQEAAzMzPY2toapc9IZoUukNy+fVs+KObOnQsLCwv88MMPWLx4MZydneVmluTkZCxatAjW1tZ6TaecmJiI1q1bZ1mvQ7v3P5ARSvr374/AwEAsXLgw29qW+/fvo0mTJvL6E5p1NLSrqNPS0uSDOi4uDi1atMjSxvs6mjKpVCp8++23CAkJQfv27eUahujoaEyYMAEVKlSAq6srgoKC4O3tbdDY/NzaQ2/duiWHEu2bVEREBHx9fbN9Wrp8+TKWL1+ucyPVDiVr1qxBkyZNMGbMGPnpIyYmBiNGjDBoFsqoqCjY2tpi5syZSExMxLVr11CuXDm8/fbb8nu2b9+OFi1ayFWyXl5eaNasWbbDKDV/57i4OHh5eeHgwYM4f/483n33XQQFBSEgIAD79u2Tf+fHjx9j+PDh8PDwMHh4nabss2bNwqNHj3D79m20bdsWHh4eci/1v//+Gy1btkSVKlXkC5qmBux1tQ2muLiZ6pgvyP1uCM3staNHj5ZD2OLFi2FmZpbtiqoHDx6Em5sbzpw5I++D69evo2HDhjqhsGHDhnBxcUGpUqXQpk0b9O/fH69evTJodt1z587JC6tp064+T0tLM2iCPmOfT9kx9X7J7Nq1a3Bzc8P06dPl7fTt2xfjxo3LUuM7YcIENGzYUB46W5Dl1/57Xr16FYcPH4YkSVnm0Vq3bh1cXV0xceJEo6yonLnv3K5du1C9enX5QTO7Zqb+/fujdOnSJuuzVagCyb1791CpUiWd6W41qd7c3Fw+ADQXwT179sDd3V2vBHr37l1UrVoVe/bsyfbftXf+ixcv0LVrV4SGhubYQatt27bw9fXFN998gyVLlqBly5a5fn5CQkKutTfatNvctQNRz549YW1tjbZt28qh7dGjR4iOjsbHH3+M//u//zNoQhpD20MfPnwoX9hyuhls2LABkiRh0aJFOjcizYmZnJyMWbNmoXbt2jojLfSZhEdDswhW5hl1Q0JC4Orqirt37+q8fvz4cfz000+4cuVKthebqKgoDBw4EA8fPkRqairGjx8vV9Nr9r+FhYU8edT69esRHR2NtLQ0vYbvabt58yZKlCihU7MHAMeOHYONjQ1+/fVX+bXTp0+jbdu2WYb75caUFzdjH/MFud/z4qeffoKtrS0mTJiA+fPnw9nZOddFAbWbudRqNZ48eYJ27dph4cKFAID33nsPLi4uuHv3LqKjozFv3jx4e3vnOIopO5cuXYKtrS0+/PBDnddXr16Nr776Kkt/HX0m6DP2+ZQbU+2XzFJSUrBgwQKMHDlSZ3DDgAED0KFDBwQHB2PmzJk6S1EMHz4cISEhuXb0NXb5L1++DEtLSyxatAhAxr1NMzT+q6++wrJly3TK+OWXX8LGxgbTp0/P14rWMTEx8Pf3x08//aTTFaF///4oX768/HfVrin69NNPDZ5J3FCFKpDExcXBxcVFZ3VQIKPHuiRJWSb6mTJlCho1aqRXO9bhw4dhZmYmH2zZJVfNXAxARk1J5jbRmzdvYs2aNXKVVvfu3dGoUSOEh4dDkiQ0btwYvr6+aNSoEerXr4/atWvDx8cHLVq00Hv59mvXriEiIiJLbcGyZcvg5OSEKVOmICgoSKf5Jq8MaQ8NDQ2Vn4yye8p69OgRTp06hf/++w/bt2+HJElYsGBBtqFEpVLB1tZWr7VusnP+/HlIkoQBAwbIN7xFixZBkiRUrVoVHTt2RI8ePTB//nwkJCTk+oQeGRkJc3Nznarv77//Hvb29nK4GzhwIFxdXbF//3589NFHcHBwQHBwsN7DtbVdvnwZLi4u6NGjB44fPy4f57/99hscHR11hrUDGQtwde3aVa9hsqa4uJnqmC/o/a6PmzdvYtWqVTojBg4dOgRHR8ds11HRuHv3Lg4ePIgDBw5kWZJh3759cHFxQUhICMqXL6/TJJmammpQp8CEhAT4+vrCx8dHJ5AtWrQIlpaWWSaJfPr0qV4T9BnzfNJ2//79AtkvOYmMjNQJ4rNnz4aNjQ0iIiIwceJENGvWDO3atdMpm3bHUFOXPyUlBUOGDIGVlZXO4qKargmaST7d3Nx0OnNv3bo138PcNQ/6Hh4emDNnjtw8qlKpUKtWLXTp0iXLPFy///67Xteh/Cg0gUQzJ0Tp0qWzvdHOmjUL5ubm8hwhs2bNgp2dnd7Vtjdu3IC9vb3OvAiZawM2b96Mxo0bZ3vCaarBO3furLOA3zvvvANJkhAQEICpU6di3rx5WLp0KZYsWYI5c+Zg2bJlBg0F1cwtMWbMGPmPv2jRIjg5Ock1R9988w2aNGmCDh065DuU6Nseeu/ePQwePDjb6tlLly4hMDBQZ+TNihUr5FCifQNJSUlBfHw8goODdU5CfSQlJck1KX/99Resra0xZMgQTJs2Dc7Ozti/fz/Onz+PU6dOYdSoUfDz84OdnR26dOmCpKSkLH/vS5cuwcbGJttZGt977z25I5yLi4vOcXb27NksT42v8+zZM3kui3PnzsHT0xPt27fH5cuXcefOHZQvXx4TJkyQ369dVk1fptyY4uJmqmO+IPe7vqKiouDp6YnQ0FB5XSiN33//XZ4WX9MJXvvnatSoIY/WaN26tbz/NW3zbdu2RbVq1YzSxDRv3jw0bNgQEyZMQEpKCpYvXw5nZ+ccR7isXr06x+Gaxj6ftBX0ftF49eoVHj9+nKW54cmTJ+jRo4dO/599+/bBxsYm2/VeTF3+//77D2q1GidPnkS/fv3g4+ODU6dOYdWqVfIMr+np6Xj+/DmaNWuGBg0a5GsUTWYqlQrvvvsuQkNDERERgZCQEAwZMgRJSUn47LPPEBAQIE90WJAUDyTanUv/+usvlClTJscqs9mzZ8PW1hbNmzdHyZIlc+3Q9vLlS6Slpck9ph8/fowGDRqgbt268ho4mUVERGDo0KFZmg9iYmLg5OSEqVOnZjvvQ9++feHt7Y3t27frdfPIjqaG4cGDB/jpp59QsWJFTJ8+HRMnTkTp0qV1LjhqtRrbtm1DrVq10KNHjzwNKTO0PfTly5fZPkVfvHgRpUqVwrRp0xAXF6ez71auXCmvn6G932bPng0vLy+9m7CAjNqzBg0a4Ndff5Wr848fPw47O7ssnZ81v196ejq++uqrbDtpXbhwAWXKlIG3t7f8mvbvt3r1atjY2MDLyyvf7aX3799HYGAg1q9fLx+PZ8+elW+ALi4uOgt3Gfr3NMXFzVTHfEHud33FxMSgdOnSmDx5co775cCBA7CxscGQIUPkUHL+/HnY2dkhIiIC9+/fx549e2BlZYXx48fr/OySJUtQoUIF+W+f3yGsixYtQv369REUFAQHBwe506p2QPj6669z3X/GPp+0KbVfoqOj8fbbb6NWrVrw8/PDtm3b5AUQgf/1sdHU1J47dw5169bN0mRp6vJfv34dnp6e8uf+/fff6NOnDypUqAArKyu5aU3zd1m5ciU8PDyMGkhSU1OxceNGDB8+HEePHsXp06fh5eWF7t27Y8mSJahRowb69+9vkuHWuVE0kDx9+hRly5bF1KlTAWRMFe7u7o4nT57k2Blo5syZsLGxybUdKzo6Gu3bt0dQUBDatGkjJ9izZ8/C3t4eQUFBOsPf/vvvP0yaNAmurq5Zpr599eoVunXrlmXm15SUFNy8eVNu1hk8eDBq1KiBLVu26MzqqA9NDUOrVq3kGoaNGzeiXLlysLS01FmpVruD6K5du/I1JC6/7aH//fcfgoKCdBbAAnRvMKtWrYKjoyO6du2KAQMGYODAgVmefPXl4eEBLy8v/PHHH/LJeubMGdjY2OC9997T6T+TW2eyyMhI2NnZyW3ko0ePzrbsgYGB6NOnj8HlzE6rVq1Qp04dbNq0Sb6AnTt3Dj4+PjoLRQKGDQs3xcXNVMe8Evv9dZKSktCpU6csnd3T0tIQHx+P+Ph4ucb0wIEDcHBwQO/evXHu3DlYWFjo1GoBGcNs69Wrp3MjBIA6derINZGGePDgAf7880/8888/OrUzH3/8Mdzc3NC1a1e574r2om4lS5Z87ZTfxjqftEVHRxfIfsksMjISpUqVQr9+/fDhhx8iICAgS5+fzOfVlClT0KRJE52H4oIo/++//w5vb2+d/oF///033n33XVStWlXuR6Yp79ChQxEeHm7wfUXbtWvX0LFjR1y4cEH+HW7fvo2QkBBMnjwZQMbfeMqUKfjggw9Qrlw5SJKUZUSSqSkWSDTJ6+OPP4aNjQ3mz5+Pv/76S+fJWfNUp3kCu3HjBlJTU3PtSBUZGQkHBweMHDkS06dPx9tvv423335b/iP8+uuvKF26NMqWLYuWLVuiS5cuaNWqFSpWrJhtyElJSUFQUJDOpEOHDh3C2LFj4eDggEqVKqFLly4AMi7QLi4uOktEv07mGgbtC/N3332H8uXLY/z48TrV6nnpeW6K9tBLly6hevXqOHLkSJYkrT0Hyt69ezFq1CiEhYVh0qRJBq13oJktUqNhw4aoXr26zkX0+PHjsLa2Rp8+fV5b63Lq1ClYWlpizpw5SEtLw+eff44yZcro3Bw1x9vnn3+OoKAgg0b/5KZLly7w8fHRCSWapoKOHTtm6T+iD1Nc3ExxzCu533Pz/PlzNGzYUKdG4KeffsIHH3wABwcHuLm5oWfPnnKtwL59+1C+fHls3boVkiThk08+kZskNWtc1atXD71798a8efPkYbPvv/8+OnXqpDPHxetojg1fX19IkoTJkyfr1EAuXboUfn5+GDdunBweZs6cmeuQTGOfT9rUajV+/PFHk++XzKKiomBvb59lvphKlSqhW7duWd5/9+5dTJ8+HY6OjjpN0AVV/m3btqFq1apZahY1DxM+Pj7yfCAzZ87MdWFLfajVahw7dgwVKlRAtWrVMHHiRLk7wOnTp2Fvb4+tW7cCyDg+zp07h6lTp8LS0jLPM0/nlSKB5NKlS1iyZAmSk5OhUqmwdu1amJmZoVGjRnB2doaLiwu8vb3h7e2NihUronr16vDw8ICvr2+uHdouXLgAOzs7nbbplStXomPHjnjy5Il8I9YckOHh4Wjfvj0WLVqUYzVkYmIivLy8MHjwYMTExGDBggXw9PRE165dsXLlSnz55Zdwd3eXO+f169dP7wWx9Klh2LJlCypWrIjRo0fnuUORqdpDv/nmG1hYWGRZhErby5cvdQKIIVWAV65cwahRo9C5c2edfg/+/v7ZXkRLliyJ9u3b59rP4OjRozo3wadPn2Z7cwQynirMzc11lpDX182bN7F8+XLMnz9fnqIeyD6UnD17Fj4+PggJCclxVs2cmOLiZopjvqD2u6FSU1Ph6+uLLl264OnTp5g/fz48PT3RuXNnfPHFF1i2bBn8/PywbNkyqNVqREVFoV+/fgAgjzpYvXo15syZAycnJ2zduhW///471q9fj27duqFs2bLw8/NDhw4dDFps8dy5c7Czs5OHZmtWU80c0hYsWAA/Pz9MmzYNI0eOzHVhNVOcTxqaEVOm3i+ZpaamonXr1pAkSb43aB6gevTogd69e+sEsLNnz6Jdu3bw8fHRmR7B1OXXDixffPEFPDw8sn3Q05y39erVQ6dOnWBra5uvuXZiYmIwbtw4+fvJkyejcePGcHV1xc6dO5GcnIw1a9agV69eWR4U8zMNfV4VeCCJjIyEJEk6nUtVKhU2bNgAJycnuLu749tvv8Wnn36K1atXY/Xq1fLCerkNT3zw4AEcHBwQFhamUzU3adIklC5dGjVq1EDJkiUxdepUg4aYAhm1KprJq+zt7bFu3To5HKSkpCAsLCzLME595FbDoL32x9atW1G5cmX079/foNU/AdO2hx4/fhw2NjY6N9zMPv30U7Rq1crgvjWRkZEoW7YsOnXqhF69esHS0lLnIqpZ4E/7Inr06FG4uLjotb4L8L9ag8TExBxvjsuXLzf4KUEz6VlQUBCqVasGOzs7rF69Wv737t27Zwkl//zzD/z9/fV6Ii2Ii5upjnnAdPs9r/bv349y5cqhYsWKcHBwwNq1a3VqJJs0aYJ33nlHHhk0a9Ys+d80nbclScqypEVycjLu3LmDiIgIg5pWo6OjYWVllWVhwqCgIGzatAnr16/X6VO2ePFilC1bFvb29jlOLGjK86mg9ktOrl27Bk9PTwQFBcnz6ty/fx82NjY6553Gb7/9ptMcZery379/H23atJHLsnnzZgQEBMj/rlnDS+PEiRNo3749ypQpk68lNDSfVbNmTZ0Rb//88w+GDh0KMzMzDBgwAFOnTkX//v3lWk5NWYyxsrKhCjSQaNqQM0/oA2RcZL/88ktYWFhg1apVedp+u3bt4OfnJw8PXrZsGUqWLIlNmzbh119/lUccZO6wpc+Ov337Nk6fPp2lw61arUb37t0xY8YMgxcQ07eG4e7du9izZw+8vb0NWlrd1O2hd+/eRbly5dChQwedCba098GECRPkVYP1pZmyW3OcpKWlYdSoURg7dqxODVlISAiqVKmCP//8Uw48eR0uqH1z1A5rhp6UmgCoWX8pMjISPj4+8PX11XnS7Nq1K2rVqoUtW7bIv5M+oa0gL26mOOYzM9Z+19f9+/dx6NAhfP/99/LvlZycjPv37+PIkSNZ5hJJTk5G9+7d8cEHH+iMDNI+V7/88ktIkoSlS5dmOyGgIVJTU9G/f39IkqRzTs2bNw+SJKFly5YoW7YsKlSogBUrVsj/vn79+hyf1E15PmUeMWWq/fI6N27cgIeHB0JDQ/H333+jcuXKGDFihPzvOR2npi5/SkoKHjx4gPbt2yMwMBA7duzAmjVr0L59+xy3k5KSgpiYGKOMKDt8+DDs7e2zndZ/165dCAsLk2vO3d3dTTakXl8FFkiio6NhaWkpd57U2Llzp9xZS1N9ZG5urjPxjyEXvS5duqBOnTro0aMHnJ2ddSZZU6lU8Pb2xvvvv2+E3yhjezNmzICrq2uexoXrU8OwcuVKtGrVCoBhVWgF1R66e/duWFtbo2/fvjq9+l++fImIiAi4u7sbVCV7+/ZtlClTRl5dVKNnz56oU6cOvLy80LJlS3kYarNmzeDk5CT3v8jvzfGLL77Idsrw/JS9RYsWqFixIuLj43VC4HvvvYdy5cph27Zteh3jSl/cgPwf89nJ737XV1RUFGrVqgVPT084OzujVq1aue4XtVqNmTNnwsXFBU5OTllGBmnfvLSHuee1qlvz1H79+nWEh4fLSwesWbMGTk5O8pN6bGwsQkND0bhx49euH2PK8ymnEVPG3i+ZZV6qQPPfmzdvwtPTU15eJPP7C7r8d+7cQbdu3RAZGYn4+Hj07t0bYWFh8PPzg5OTE0JDQ9GsWTM0b94cwcHB8gi4Dh065Loem77UajX+++8/VKlSRR7anHlfXLlyBV9++SUqVKiQJQQrocACyZQpUyBJks5TmmbyHe3OpCqVCp9//rmcTHOTUzt9r1695E5g2hfqV69eITg4OEsoyostW7Zg9OjRcHFxyfPMdfrWMEyaNMmgUFaQ7blqtRrr1q2DhYUFPD090b9/fwwfPhwdOnRAuXLlDN43N2/ehL+/Pzp06CAPZ1y4cCHs7Ozw4YcfYsOGDfD29kaVKlXkC3jLli2NNmHP06dPsXHjxjy1a2dXdk2tXJ06dRAWFobQ0FBMnToVly9fxosXLzBs2DC9FsZS+uIGGOeYz0l+9rs+NLWzU6ZMwY0bN7B582ZYWlqib9++8k1I+/z65ptvMGrUKJQqVQo2NjbZjgzSniofyHh4sLS0lKdfN0RycjIaNmwIDw8PAP9bS6pkyZKws7OTr5uaMs6dOxc1a9Z87cODqc6n3EZMGXO/ZHbp0iUMHjwYf/31l06NonYoqVu3Lpo0aZJrWCuI8l+/fh2NGjXC22+/jYsXL+Lu3bvo06ePXMP14YcfYsSIERg5ciQmT56M0aNH44MPPsCFCxcM2SU6Ll++jF27duHYsWPyw0/VqlWxfPlynfdlnsLhyZMnRmk+yy+TBxLNjTYlJQW9e/dGiRIlEBsbi88++wxlypTJdkIflUqFr776KtfRGK9rp+/Vqxd8fX2xadMmuf1s5syZqFSpUr5vXpcvX0ZISAg6d+5s0IiR7Bi7hkGp9ty///4b3bp1g5+fH4KCgjBlypQ8P0HHxsaiTZs26NChAwYNGoRy5crpzDAaFxcHSZJ0RoEYU35qWTKXvWzZsti9ezcePnyIY8eOYf369fD09ETZsmURFhamd38mJS5u2ox5zOfEVM00169fh52dnU7TpVqtRpUqVdC6dess7//777/Ru3dvNGvWDBYWFrmODMp889JMYmjohIXp6en4448/UKtWLfj7+wPIuHZ27doVpUuXlofda46X4cOHo3379lmaXbNj7PNJnxFTxtov2tLS0hAeHg4rKytUqlQJgwcPlm+02p918+ZNVK9eHc2aNcsyorCgyx8bG4vWrVsjLCwMFy5cwMOHD9GnTx8EBgZmmYAvv9RqNYYMGYLq1aujZMmS8PDwQM2aNVGzZk306tULP/zwAx4/fqz3rOFKMGkgSU5ORqNGjVCtWjWkp6fLbc9mZmawsbExeESBhqHt9N999x2mTp0Ka2vrfHcS0nj48KFRqiCNWcOgdHuuMSfRuXLlClq1agVbW1ssW7YMQMZFOyUlBXfv3kWdOnWwa9cu+fXCRLvs2dXyPX/+HH/99ZfBS4YX5MUtO8Y65gvaJ598Ii8mqGke1tTOvvXWWxg+fDhGjBiBU6dOyefGvXv38OOPP+o1MijzzSun9a9yojl+09PTceLECdSsWVMnlLRq1Qpubm5yp/4ZM2YYPBTUmOeTviOm8rtfsrN27VrMnTsXp06dwsqVK1GxYkWEhYVh3rx5On2dbt68CScnJ7Ru3TrLNa6gy5/5vL137x569+6Nxo0bY+PGjfL78tsfSyMpKQm3b9/GH3/8gYULF6J58+aQJAnVq1eHm5sbqlevjiFDhmDChAl5vgebikkDiSb1+/j4oH79+khPT0dqaiqGDx+uM2mTIX+EvLTTS5KEkiVLGi2MmEJ+axiUas/Vlnntg/y6du0awsLCEB4errNOx8yZM1G1atVsn34KC+2ya096lp8FsYCCv7gVZZqn2aSkJMyfPx+NGjXC9OnTMXfuXJQpUwafffYZ/vnnHyxYsADdunVDmTJl4OrqipkzZ2bZ1utGBqWlpekEC31o145pL2KWXSgJCwuDh4cHBg0alOehoKY4n0yxX3ITExODUqVKYe/evQAy9tvKlStha2uLqlWrYv78+fLvFhcX99rQX1Dl1z5vL168iPv37+O9996Dt7c3tm3bluftamiXLXMAW7JkCTw8PJCQkICff/4ZixcvxqBBg+Dj42PwQ5GpmbzJRq1WyyeYJpSkpaWhR48eKFGiBI4fPy6/Tx95aaefOHGi0aquTSmvNQxKtecWBE11c+vWrXH27FksXrz4tTP1FhbaZdccq8barikvbm+CM2fOwMLCQr6+JCcnY+7cuahduzbMzc3lG5q248eP47PPPnvtcGPtm5f2HA+GuHXrFsaOHauz6J12KDl58iS8vLzQpk0bABmDAlq2bAlbW9t8PViZ8nwyxn7J7ObNm1n+VsuWLUOnTp3kh6o+ffrA29sbU6ZMkRd9HDNmjMEBwhTl16Z93l66dAl3797F4MGD89VUHhsbi0mTJqF3795ZZlXVHE9//PEHqlWrlqW/kSlGO+WX0QNJfHx8lhknU1JS8Pfff6N69eqoV6+e3HzTs2dPlCpVSmckjD4MaacPDQ01eN4RpeSlhkGp9tyCFBsbi3bt2slT6ednoqCCpil748aN8zQTa27bNfbF7U0RGRkJe3t7TJw4EcD/LsxJSUmYN28e/Pz8MGHCBLktPS+dfrVHBmmWvjBEVFQUqlatimHDhukEIE1ZU1JSsHfvXnh5ecmrDJ85c0bvOXZyY8rzKb/7Rdu9e/fkWl/NTKJAxky6fn5+iI+Px5AhQ1C+fHl5Qsf79+/j+++/z3M/J2OWPzuxsbF4++230ahRI8TExOQrFERGRqJcuXLybOSSJOnUkmrcvHkTtra2+PHHH3VeL4w1qEYNJLdv30bp0qUhSRJCQkIQERGBX3/9VWfyp7p166JOnTpyTUl4eDgqVqyYZRXN1zFVO31Ro2R7bkG6fPkyOnToUOBTGRtDTEwMunXrpjMZkzEY8+L2poiKioKtrW2WZhfNiAtNTUmjRo0wZswYOZTkZd/ld2TQuXPnUK9ePQwaNCjbUPL8+XNUrlw5ywRpxmDK88lYI6Z+++03SJIEf39/dOrUSWdNr+7du0OSJFSoUEGnlskYTD3iKyYmBl27ds3X9UB7bhnNcgB9+/ZFREREloD98uVLNGjQQGf/FVZGDSS3bt1C3bp14enpiQYNGqBfv36wsbFB3bp18e6772LHjh3YuXMnatSogZYtWwLIuBDkdY4EU7XTF1UF3Z5b0IpKTVd2jDX0NjNjXNzeFPHx8ahSpQqaNm2q8/qSJUvQtWtXuYpfE0qCgoIwaNAggx+GtOX3PDp79qwcSrRH2aWmpuLVq1do27at3Axn7HPWlOeTsco6YMAA1KlTB127dkWLFi3kGoBTp07B09NTnuTS2KvSmvr6mJ/rQVxcHJycnNC3b1+d17t06YLAwEDUrl0b77zzDnbv3i3/W+vWrfHOO+8AKNzXfqM32Vy9ehWdO3dGx44dcfLkScTFxWHbtm0IDAxEw4YNYWtrKy8WpVmgKz9M1U5f1Jm6PZQKD1OFnaJGc+1p1qwZvvjiCwAZC9A5ODjgl19+AfC/G1dycjKmTJmCVq1aGTT7sSloQkn//v3lvhwqlQqzZs2Cu7u74pNVKUEzx8iPP/6I999/Hz/99BO6dOmC4OBgOaA1adIE7777rpLFLHCvXr3CwYMHUbt2bfTo0UMeabVw4ULY2tpi+vTpWLp0KSpWrIi6devKIXfp0qUmq/ExJpN0ar18+TJat26NVq1ayTPEAUBCQgI2b96M6dOno169ekbrmGiqdvqiztTtoUSFwePHj+WRdZcuXUK/fv3QtGlTdOjQAWXKlMFvv/2m837tUJJ5WnylREZGIjg4GF5eXujYsSO6dOkCNzc3nQXg3nS3b9/Gnj17dF579OgRvLy8sGrVKjx69AhdunRBUFAQDh48iOPHj8PR0RH79+9XpsAFTFMr9O+//2Lnzp1o0aIFevbsiVGjRmWZW+bcuXOQJEmn701RYLJRNppOd61bt86206qxm1VM1U5f1Jm6PZRISf/++y/CwsIwduxYua+aJpQ4OzvrNFlqV+sbu4rfGOLi4vDpp5+ie/fuWLBggdGm5i8KtPsfvv3229ixY4d8zdq3bx+Cg4Px6NEjREdHo0uXLmjVqhUmTJiAzp07F+rh/8ai6ag9atQo+bVvv/0WQUFBsLKywpo1awBkNMunpqbixo0b8PHxkTtEF+ZmGm0mHfar3ZyiGX5nSqy6zl5RORiJDJWSkoIRI0YgMDAQ06dPl0PJ5cuX0a9fPwQGBmL9+vXy+wtjEKGM/ocNGjRAQEAA6tevj0GDBsHd3R3r1q3Djh070K5dOxw4cAAAcPHiRYSGhqJfv35GW6OpMNNMBKpZHFH7en748GE0a9YMnTt3lkcaARmT53l4eBS5/SMBgDChq1evivHjx4vHjx+L5cuXi8aNG5vy44iomFCr1cLc3FykpqaK6dOni2PHjonQ0FAxZcoUYW9vL6Kjo8WSJUtEbGysGDhwoBg4cKDSRaZcXL16VUydOlWkp6eL9957T5iZmYkVK1aIUqVKib179wp/f3/xxx9/CCsrKxETEyPs7e1FpUqVlC62Sd25c0fUq1dPtGjRQuzYsUN+/eOPPxaPHz8WCxcuFNu3bxdffPGFsLe3Fx9//LHYvXu3mD17tjhx4oSoW7eucoXPi4JIPWxOISJj0V76QCM5ORkTJkyAv7+/Tk3JpUuXMHDgQHh7e2eZOIoKn8uXLyM8PBxhYWG4cuUKXrx4gRMnTqBdu3bYvHkzgOJV45vT4ojaHbUBYPv27QgLC0OFChVgZWVVpOZq0mbyGhKNlJQUYWVlVRAfRURvqOjoaOHr6ytCQkJExYoVxfDhw0Xp0qWFp6enUKvVYt68eeLAgQMiNDRUTJ06VTg4OIgLFy6IdevWiUmTJokqVaoo/SvQa1y9elWMGjVKCCHErFmzRGBgoMIlUtbVq1fF6NGjhZWVlXBxcRF79+4VW7ZsEWFhYSI9PV2YmZkJIYTYunWr2LJli/j444+Fr6+vwqXOmwILJERE+ZGeni4OHz4s3n77bVGpUiXh6uoqnjx5Ip4/fy5atWolmjZtKtq0aSMWLFgg7t69K/z8/MTEiROFvb09H4iKGM1NGICYMWOGCAoKUrpIioqNjRWjRo0Sf/75p5g3b56YMGGC/G/aoeTFixeiZMmSShUz3xhIiKjQu3Dhgli5cqXYsGGD+Oqrr8SgQYPEihUrhL+/v3j8+LH45ptvxD///CNKliwpJEkSd+7cEUlJSWLu3Lli0qRJQgghJElS+LcgQ7D/oa7r16+LESNGCHNzczFt2jQ5pGlu4W/C8W2mdAGIiHJz/vx54efnJypWrCiEEGLAgAHi448/FuPGjRO//fabaN++vdi+fbu4fv26mDZtmujdu7dwcXERrq6uomPHjkKSpDfiYl3c1KhRQyxdulSuDSvuqlevLlatWiUAiPnz54vjx48LIcQbdXyzhoSICq3o6GhRv359MWXKFDFnzhyd6umVK1eKcePGifnz54thw4YJZ2dn+ecePnwobGxshKOjo1JFJyNhc5uuN7nmyELpAhARZefixYuiefPmomrVqmLOnDlCiIz2ciGEMDMzE2PGjBGSJImxY8cKSZLEiBEj5ADi4uKiVLHJyBhGdGlqjmbOnPnG1RyxhoSICp3z58+LJk2aiIYNG4rY2FjRrVs3sXLlSiFExvwjkiTJNSWffvqpmDhxopg6daqYOHGicHBwULLoRAXiTaw5Yh8SIipUTp8+Lfz9/cXkyZPFL7/8ImbPni2+/fZbMWbMGCGEEObm5gKAXFsyevRoMW/ePLFq1SqRmpqqZNGJCsybFkaEYA0JERUyx44dE7t375ZrRBITE8WOHTvE9OnTRe/evXOsKUlISBBOTk6KlZuI8oeBhIgKLQBCkiTx7NkzsX379mxDiZmZmZAkSX4vERVN7NRKRIWWJmA4ODiIXr16CSGEmD59ujA3NxeffPKJMDc3z/JeIiqaGEiIqEjQhBIzMzMxZMgQYW1tLRYuXKh0sYjISBhIiKjIcHBwEN27dxeWlpYiICBA6eIQkRGxDwkRFTnsL0L05uGwXyIqchhGiN48DCRERESkOAYSIiIiUhwDCRERESmOgYSIiIgUx0BCREREimMgISIiIsUxkBAREZHiGEiIiIhIcQwkREREpDgGEiIiIlIcAwkREREpjoGEiIiIFPf/AOEiGcXgcQc/AAAAAElFTkSuQmCC",
      "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, labels=distance_matrix.index, color_threshold=2\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "0a73d648-3099-47e8-b7d2-6c9e54b4b8c5",
   "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": 266,
   "id": "fa1a6087-7597-4ac7-8c6b-ebc945044ef8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save grpofiler results\n",
    "with open(\n",
    "    \"figures/pert_embedding_cluster_gprofiler/GSFA_inhouse_pert_emb_gprofiler.pickle\",\n",
    "    \"wb\",\n",
    ") as handle:\n",
    "    pickle.dump(cluster_process, handle, protocol=pickle.HIGHEST_PROTOCOL)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b11d8f8-00ca-47cf-acde-53d5d1c88bef",
   "metadata": {},
   "source": [
    "## PR Curve+AUC\n",
    "- AUCPR using prior graph as binary label"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 268,
   "id": "16fceeb8-b9a3-495f-a799-f733d21cba7b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "               AUC\n",
      "StringDB  0.685749\n",
      "CORUM     0.536772\n",
      "                F1  Threshold\n",
      "StringDB  0.689548   0.011192\n",
      "CORUM     0.632953   0.016780\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/GSFA_inhouse_aucpr.csv\")\n",
    "f1_df.to_csv(\"figures/pert_embedding_recall/GSFA_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": {},
   "outputs": [],
   "source": [
    "factor_to_go = factor_enrichment_gsea(adata, W.T, 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": 98,
   "id": "a9e4b2f0-41d1-434d-b67e-94c6bef8ddd4",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\n",
    "    \"figures/factor_enrichments/GSFA_inhouse_factor_enrichment.pickle\", \"wb\"\n",
    ") as handle:\n",
    "    pickle.dump(filtered_factor_to_go, handle, protocol=pickle.HIGHEST_PROTOCOL)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c48ac5d8-e953-4b82-b6ca-9e59bc8e49bf",
   "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": 14,
   "id": "6867de68-b31b-452d-8f3f-167e5af8d2c9",
   "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/GSFA_inhouse_factor_enrichment.pickle\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "bf3b1040-d325-4154-9f9e-92d4ec8d9ab7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get neighbors for each perturbation\n",
    "pert_neighbors = retrieve_stringdb_neighbors(pert_labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "781275db-dc25-4733-84e6-010ebbc17cdd",
   "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",
    "# subset to biological processes\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 df.index}\n",
    "for goterm in filtered_go_terms:\n",
    "    for pert in df.index:\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": 19,
   "id": "acdb19ec-92d3-4d85-9a0b-23c3711b879b",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_32012/994851414.py:11: RuntimeWarning: divide by zero encountered in log\n",
      "  delta_loading = np.abs(np.log(np.abs(pert_loading))-np.log(np.abs(ctrl_loading)))\n"
     ]
    }
   ],
   "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 pert_emb_df.index:\n",
    "    if pert in [\"ctrl\", \"intergenic\", \"basal\"]:\n",
    "        continue\n",
    "    # get top factors\n",
    "    pert_emb_df = pd.DataFrame(beta[:-1], index=pert_labels)\n",
    "    pert_loading = pert_emb_df.loc[pert].to_numpy()[1:].astype(float)\n",
    "    ctrl_loading = pert_emb_df.loc[\"intergenic\"].to_numpy()[1:].astype(float)\n",
    "    delta_loading = np.abs(np.log(np.abs(pert_loading)) - np.log(np.abs(ctrl_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": 20,
   "id": "26781b0d-61c6-45ca-9b42-ee77c73ceed2",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overlap for RIPK1: 9 out of 38 in researchDB\n",
      "P-value for RIPK1: 1.1308619075120921e-07\n",
      "Overlap for MAP3K7: 0 out of 23 in researchDB\n",
      "P-value for MAP3K7: 1.0\n",
      "Overlap for IKBKB: 7 out of 24 in researchDB\n",
      "P-value for IKBKB: 3.830013304044622e-07\n",
      "Overlap for SKP1: 0 out of 7 in researchDB\n",
      "P-value for SKP1: 1.0\n",
      "Overlap for BIRC3: 3 out of 20 in researchDB\n",
      "P-value for BIRC3: 0.007686716976516679\n",
      "Overlap for TAB2: 0 out of 5 in researchDB\n",
      "P-value for TAB2: 1.0\n",
      "Overlap for NFKBIA: 6 out of 28 in researchDB\n",
      "P-value for NFKBIA: 1.8192895607157648e-05\n",
      "Overlap for SKP2: 1 out of 6 in researchDB\n",
      "P-value for SKP2: 0.11702982854059442\n",
      "Overlap for TRAF6: 4 out of 46 in researchDB\n",
      "P-value for TRAF6: 0.018440536255468803\n",
      "Overlap for FBXW11: 5 out of 13 in researchDB\n",
      "P-value for FBXW11: 3.4478778936666224e-06\n",
      "Overlap for TRAF5: 4 out of 14 in researchDB\n",
      "P-value for TRAF5: 0.00020482544648502187\n",
      "Overlap for CHUK: 7 out of 21 in researchDB\n",
      "P-value for CHUK: 1.84359744108617e-07\n",
      "Overlap for RBCK1: 0 out of 11 in researchDB\n",
      "P-value for RBCK1: 1.0\n",
      "Overlap for TAB1: 0 out of 9 in researchDB\n",
      "P-value for TAB1: 1.0\n",
      "Overlap for IKBKG: 5 out of 17 in researchDB\n",
      "P-value for IKBKG: 2.6711340009329823e-05\n",
      "Overlap for CUL1: 2 out of 7 in researchDB\n",
      "P-value for CUL1: 0.008357243957233893\n",
      "Overlap for TNFRSF1A: 7 out of 21 in researchDB\n",
      "P-value for TNFRSF1A: 2.298789369209975e-07\n",
      "Overlap for TRAF2: 1 out of 29 in researchDB\n",
      "P-value for TRAF2: 0.2850611728304155\n",
      "Overlap for TRADD: 6 out of 13 in researchDB\n",
      "P-value for TRADD: 1.0736863763836888e-07\n",
      "Overlap for BIRC2: 4 out of 30 in researchDB\n",
      "P-value for BIRC2: 0.0024805517666413985\n",
      "Overlap for UBE2N: 0 out of 16 in researchDB\n",
      "P-value for UBE2N: 1.0\n"
     ]
    }
   ],
   "source": [
    "pvals = perturbation_signal_recovery(\n",
    "    pert_to_go,\n",
    "    model_pert_to_go,\n",
    "    list(filtered_go_terms.keys()),\n",
    "    list(pert_emb_df.index),\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 262,
   "id": "c3e93509-ea30-44d2-b978-5627a449c255",
   "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/GSFA_inhouse_hypergeo_neighbors_recovery_pvalues.csv\",\n",
    "    header=False,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c7cbc9b-d955-4778-82ea-9e31b1934f59",
   "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
}