{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "314c0d27-5252-4343-a323-5a69fcdc7959",
   "metadata": {},
   "source": [
    "## Enrichment or depletion of spatial factors in cortical layers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "023fb272-7524-4791-ad87-a4e61c50b54d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scanpy as sc\n",
    "import anndata as ad\n",
    "import umap\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.cm import get_cmap\n",
    "import scipy.io\n",
    "import matplotlib as mpl\n",
    "mpl.rcParams['pdf.fonttype'] = 42\n",
    "mpl.rcParams['ps.fonttype'] = 42\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ed4b9a6-c2ab-4936-87c2-5feb9f4d9bfa",
   "metadata": {},
   "source": [
    "### Load results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "91170c1e-e340-4599-880d-b0105935bba6",
   "metadata": {},
   "outputs": [],
   "source": [
    "res_dir = \"Results/INSPIRE_DLPFC\"\n",
    "adata = sc.read_h5ad(res_dir + \"/adata_inspire.h5ad\")\n",
    "basis_df = pd.read_csv(res_dir + \"/basis_df_inspire.csv\", index_col=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dfb8833f-f546-4e11-a68d-1823da38a731",
   "metadata": {},
   "source": [
    "### Enrichment or depletion analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "d24da7bc-ed12-4375-ad50-f6e71c8b946c",
   "metadata": {},
   "outputs": [],
   "source": [
    "layer_list = list(sorted(set(adata.obs[\"layer\"])))\n",
    "n_layer = len(layer_list)\n",
    "n_topic = 20\n",
    "layer_prop = np.zeros((n_layer, n_topic))\n",
    "for i, layer_name in enumerate(layer_list):\n",
    "    adata_tmp = adata[adata.obs[\"layer\"] == layer_name, :].copy()\n",
    "    prop_tmp = adata_tmp.obs[[\"proportion-\"+str(i) for i in range(n_topic)]]\n",
    "    layer_prop[i, :] = list(np.mean(prop_tmp, axis=0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "9abcf79c-4457-46d3-848e-308cbe9073ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "random.seed(1234)\n",
    "\n",
    "n_try = 1000\n",
    "permuted = np.zeros((layer_prop.shape[0], layer_prop.shape[1], n_try))\n",
    "prop_df = adata.obs[[\"proportion-\"+str(i) for i in range(n_topic)]]\n",
    "layer = adata.obs[\"layer\"].values.astype(str)\n",
    "\n",
    "for p in range(n_try):\n",
    "    region_props = np.ones((n_layer, n_topic))\n",
    "    for j, layer_name in enumerate(layer_list):\n",
    "        idx = np.arange(prop_df.shape[0])\n",
    "        random.shuffle(idx)\n",
    "        prop_permute = prop_df.iloc[idx, :]\n",
    "        prop_region = prop_permute.loc[layer==layer_name, :]\n",
    "        region_props[j, :] = np.mean(prop_region, axis=0)\n",
    "    permuted[:,:,p] = layer_prop - region_props"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f9fd3f4e-93cb-453d-9e73-e251795aad41",
   "metadata": {},
   "outputs": [],
   "source": [
    "diff = np.mean(permuted, axis=2) / np.std(permuted, axis=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "6adddee3-814d-4997-86e4-e3a5737c9d40",
   "metadata": {},
   "outputs": [],
   "source": [
    "diff = diff[[0,1,4,5], :]\n",
    "diff = diff[:, [16,17,19]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "5743b965-311b-4de3-bfa4-ade33d7cc9fc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPAAAAG7CAYAAAD9vRiJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABRnUlEQVR4nO3deVxUVf8H8M8sMCyyI6JiIJsimqXiRmWouKW5RO6C62NPKg/ZU09aBvrSx5bHcEvTNCQ1tRIz08fUX7ihomaW+qgEiIEmmwoIyDbf3x/XmRhmBmZgBrj6fb9e81LOPfecMwPfueeee+49EiIiMMZESdrUDWCM1R8HMGMixgHMmIhxADMmYhzAjIkYBzBjIsYBzJiIcQAzJmLypm4A00+pVOL27duws7ODRCJp6uawBiIiFBUVoU2bNpBKTXPs5ABuxm7fvo127do1dTOYiWVmZsLDw8MkZXEAN2N2dnYAhF+4vb19E7eGNVRhYSHatWun/r2aAgdwM6bqNtvb23MAP0ZMeTrEg1iMiRgHMGMixgHMmIhxADMmYhzAjIkYBzBjIsYBzJiIcQAzJmI8kYMxfbKzgfR0oKwMUCgAb2+gVaumbpUGDmDGVJRK4MgRYNMm4NgxICdHO4+bG9CvHzBzJjBwIGCimxLqi7vQjAFAQgLg4wMMHiz8X1fwAkJ6QoKQz8cH2LOncdtZAwcwe7LdvQuMHw+88gpw86aQVlVV+z6q7TdvAmPGCPvfvWvedurBAcyeXLdvA717A99+K/xs7BoHqvzffiuUc/u2adtnAA5g9mS6exd48UXgxo26j7h1qaoSynnxxUY/EnMAsyfT668LI8yVlaYpr7JSKG/OHNOUZyAehX6M3b4NbN4MXLwIPHwItG4NTJwIhIQApn5CT34+cOAAkJcH2NoKA7UdOpi2DpNJSAB27TJ9uVVVwM6dwNixwOjRpi9fF2LNVkFBAQGggoICo/YrKiKaNIlIKiWSyYgkEiKASC4X/vXxITp61DRtzMoiioggsrAQypbJhH8Bon79iE6cME09RERUVUV08CDRyy8TtWlD5O5ONGAA0bffEpWXG16Gl9dfH4qpXxKJUH5VlVbV9f191saoAE5MTCQAel+JiYkma1h1e/bsoejoaLOUbYh///vfFBYWRu3btycA5OnpqTdvbZ8PAFq6dKnB9dbnF15URNStm2Yg1XxJpUIwHzhgcLE6paQQtWr11xdDzZdMJmz75puG1UNERHfvEj3/vOY3UfVvjE6diDIz6y7nxx/NE7g1X4cOaVVtjgCuVxd63LhxGD58uFZ6QEBAfYqr03fffYf4+HjExMSYpfy6LFy4EM7OzujWrRvu379fa96tW7fqTI+JiUFaWhpGjBhhhhb+5R//AH79tfZxGaVS+CsLCwPS0gB3d+PrKSsDBg0Susz66qqqErrqEyYA/v7A008bXw8AoKICeOkl4OxZ4efq562qylNShHODn38Ganv80KZNgEzW8IGr2shkwOefA6Gh5qtDxZhoVx2Bly9fbrJvEENERESQkU012MOHD6mioqLWPGlpaer/BwYG1noE1iUzM5OkUin16NHDqP2M/cbOztZ/NNR3JDaiQ6Dhq68Mr0cuJ5o6tX71EBHRzp2Gv6H//Kf2stzcGucI3KqVVtXmOAKbfBT67NmzmDp1Kvz9/WFjYwM7OzsEBwdjj54ZK3fu3EFkZCS8vb2hUCjg5uaG0NBQHD58GADg5eWF+Ph4AMLDwFSvo0ePqstISkrCkCFD4OjoCGtra3Tt2hVr1qwB1biuN3XqVEgkEuTm5mL69Olo1aoVrK2tkZWVVet78vb2bsAnAsTFxUGpVGLmzJkNKqcu8fHC0dVQSiXw6afCX5yxPv3U8FmElZXAV18B9+4ZXw8AYM0awypTKoW8+j6E7Gz9M6xMrZHqqlcXuqSkBHl5eRppCoUCdnZ22LNnD1JSUjBhwgR4eHggPz8f8fHxGDNmDLZv346JEyeq98nIyEBwcDCys7MRERGB7t27o7i4GGfOnMGRI0cQGhqKlStX4pNPPsGJEyc0uqeq7vqBAwcwcuRIuLq6IioqCk5OTti9ezciIyNx6dIlbNy4Uav9oaGhaNOmDRYtWoTi4mK0aNGiPh+DQYgIcXFxsLGxwYQJE2rNW1ZWhrKyMvXPhYWFRtV16ZLxo8t//gkUFACOjsbt98svxn1ZlJcLvdxevYyrB0RAcrLhld28KfTr3dy0t6WnG1l5A6Wl6W6HKRlzuK5tEGvkyJFERPTgwQOt/YqLi8nf358CAgI00ocOHUoA6JCOE/6qaqN4+rrQlZWV5OnpSXZ2dpRZbQCjsrKShgwZQgAoKSlJq5zw8HBj3rYGY7vQR44cIQA01YA+ZHR0tM7P1tAu1/jxQi/S2N5eTo7Bb0fN0tL4eo4fN74eqqoyfsRY32BWYmLjdJ9VrxqDus2mCz1jxgwcPnxY47VkyRIAgK2trTpfSUkJ8vPzUVJSgv79++Pq1avqo8rdu3dx8OBBDB48GKE6TvYNWXriwoULuHnzJqZOnarxpHuZTIaFCxcCABISErT2mz9/vnFvuAE2bdoEQPjM6rJgwQIUFBSoX5mZmUbV1aaN8TfHWFoaf/QFhGvKxqrXYgRSqXE7WlsDLVvq3qZQ1KMBDdAI9dWrC+3r64uBAwfq3JaTk4P33nsPe/fuRY6Oc4D79+/D3t4eqampICJ07dq1Pk0AAKQ/6hIFBgZqbevSpYtGnur8/PzqXacx7t27hz179qBjx4547rnn6syvUCigaMAvfdIk4JNPDM8vlwsTOywsjK9rxgwgJsawnq1UKkwVbt/e+HoAALNnA++/X3dlcjkQEaE/cBo4lmE0Hx+zV2HSQSylUonQ0FDEx8cjPDwcu3btwsGDB3H48GH1ua/y0S+B6jNyUkN9y7CxsWlw3YbYtm0bysrKDDr6mkK3bkBQkHAVwxCVlcKMwvqYNcvwepRK4fJWvc2aBbRoUXv3QiIRXpGR+vO0amX+c9JGrsukAXzp0iX89ttveOedd/Dxxx9j7NixGDx4MAYOHIiqGtfd/Pz8IJFIcPHixTrL1bcUhc+jb7grV65obbt8+bJGnqawefNmWFhYIDw8vNHq/Pxz4QBUV3BJJEBUlBDw9eHuLkzTVJVVWz3h4cCrr9avHgBCIBw4IHSP5To6jXK58Ia//hqoay5Cv36Gf/PUl1wOvPCCeet4xKQBLHv0wdQ8Ml6+fFnrMpKzszOGDh2KQ4cOqS8ZVVe9DNUo8b0a1yG6desGT09PxMfH49atW+p0pVKJ5cuXAwBGN9ac1BrOnz+PX3/9FSNGjIBbY33rA+jaFfjpJ8DB4a+DUnWqv/833gBWrGhYXVOmCFOKVfMmqh8gpVKhrvnzgS++MMHc6+Bg4MIF7S6yTCbcy5ucDIwaVXc5M2eadxIHIHRtZs0ybx2PmPRmhoCAAAQGBuKjjz5CSUkJOnTogJSUFGzYsAGdO3fGhQsXNPKvXbsWffv2xbBhw9SXkUpLS5GcnAwvLy98+OGHAIBevXph7dq1mDNnDoYOHQoLCwv0798fbm5uWLduHUaOHImgoCDMnj0bTk5OSEhIwLFjxzBr1iz07du3we9r69atuPnoZu/c3FyUl5dj6dKlAABHR0fMnTtXa5/Njw5P5r72q0uvXkBGBrBtm3BZ9OpVId3GRjga/v3vDZgVVcPYscCIEcA33wgHwJwcobfbv78QK/WZ5aWXv78wk2rFCuD334W+ubc34OpqeBkDBwJeXsLlJhOcxmmRSABPT2DAANOXrYsxQ9aGzMTKyMigsLAwcnV1JWtrawoKCqKEhAT1JZIbN25o5M/KyqLZs2dTu3btyMLCgtzc3Cg0NJSOHDmizlNZWUlRUVHk7u5OUqmUAM151ydOnKBBgwaRvb09KRQK6ty5M61atYqUSqVGXfWd0dWvXz+9l890XVIqKSkhBwcH8vDw0LgcZixTXXaoqCAqKSGq8XE8uRISzHv5KCFBZ7XmuIwkITLH1xAzhcLCQjg4OKCgoICXFzW18eOFJ2mYsjstkwkn+zt26Nxsjt8n39DPnkzr1gndb12DYvUhlwuXjT791DTlGYgDmD2ZnJ2Bo0eFi9MNHZWWyYQvg8REodxGxAHMnlxt2gBnzvx1jcvYoXJV/ldfBU6fFsprZBzA7Mnm7CycsyYkCKPHQN1HZFW329NT2G/HjkY/8qpwADMGCM+wSksDDh0SnvWsbwmVVq2EvIcOCfmbaJ6BCj/UjjEVqVR4iobq5pqcHCFIVWsj+fg03lRMA3EAM6aPm1uzC9iauAvNmIhxADMmYhzAjIkYBzBjIsYBzJiIcQAzJmIcwIyJGAcwYyLGAcyYiHEAMyZiHMCMiRgHMGMixgHMmIhxADMmYhzAjIkYBzBjIsYBzJiIcQAzJmIcwIyJGAcwYyLGAcyYiPFTKZlpVFYK65jevAmUlwOWlsKDzwMCTLf+ENPCnyyrv/v3gS+/FBYi/u034fnJNSkUwmLEkycLixM7OjZ2Kx9r3IVmxissBObNA1q3BqKigPPndQcvIKSfPy/ka90aiIwU9mcmwQHMjHP4MNCxI7B+PfDw4V/LWtdGlefhQ2FZz44dhXJYg3EAM8PFxgKDBgHZ2fVfGLuqSth/0CBg5UqTNu9JxAHMDBMbC8yfL/xfqWxYWar933iDg7iBOIBZ3Q4f/it4Te2NN5ptd5oIuHcP+OMPYbyurjOFpsCj0Kx2BQVARISwcl9Dj7y6SKVC+deuAfb2dWavqgJ+/FEYFysqAuzsgB49gMGD617W11D37wPx8cDatUBq6l/pHTsKY3CTJwv1NgvEmq2CggICQAUFBU3XiLlziWQy1TCUeV4yGdG8ebU2o6SE6IMPiDw8hF3kciILC+FfQEj/4AMhX0McOEBka0skkQiv6s1UpdnbE/30k/Flm+P3aVQAJyYmEgC9r8TERJM1rLo9e/ZQdHS0Wco2xL///W8KCwuj9u3bEwDy9PSsNf/JkyfppZdeojZt2pCVlRX5+PjQa6+9RhkZGUbV2+QBfO8ekZWVeYNX9bKyIrp/X2czcnOJgoKIpNLai5BKhXy5ufV7uwcOCGUYUo9cTmTsn7s5fp/16kKPGzcOw4cP10oPCAioVy+gLt999x3i4+MRExNjlvLrsnDhQjg7O6Nbt264f/9+rXkPHDiAESNGwNfXF5GRkXBxccGvv/6Kzz//HLt378alS5fQSt/q783Nl1/qv75ramVlQn3z5mkkl5YCw4YBFy7U3YNXKoFffhHyHzsGWFsbXv3du0BYmGFXxVTtGDUKuHULsLU1vB6TMybaVUfg5cuXm+wbxBARERFkZFMN9vDhQ6qoqKg1T1pamvr/gYGBtR6BBw0aRBYWFpRb4zAQGxtLAGj9+vUGt63Jj8BBQdr9SHO9JBKinj21mvDBB3UfEXUdIT/80Li3umKF8W9VIiHasMHwOszx+zT5KPTZs2cxdepU+Pv7w8bGBnZ2dggODsaePXt05r9z5w4iIyPh7e0NhUIBNzc3hIaG4vCjkUkvLy/Ex8cDACQSifp19OhRdRlJSUkYMmQIHB0dYW1tja5du2LNmjWgGl+lU6dOhUQiQW5uLqZPn45WrVrB2toaWVlZtb4nb29vg99/QUEBrKys4OzsrJHepk0bAICNjY3BZZlUeblx+SsrhemRjTX0SgT8+qtQ7yNVVcJAkrFjZ0qlsJ+hl6qJgDVr6vdWV682fh9TqlcXuqSkBHl5eRppCoUCdnZ22LNnD1JSUjBhwgR4eHggPz8f8fHxGDNmDLZv346JEyeq98nIyEBwcDCys7MRERGB7t27o7i4GGfOnMGRI0cQGhqKlStX4pNPPsGJEyewdetW9b6q7vqBAwcwcuRIuLq6IioqCk5OTti9ezciIyNx6dIlbNy4Uav9oaGhaNOmDRYtWoTi4mK0aNGiPh+DTgMHDkRycjIiIiLwz3/+Ey4uLvjtt9+wYMECPP300wgLCzNZXQa5cgUYPRr4/XfAzw/YswcIDKx7v6tXG6/7rFJWJoxGd+4MQBhtruO7Va/MTODQIWDo0Lrz5uQAGRnG10EkfLzFxU3YjTbmcF3bINbIkSOJiOjBgwda+xUXF5O/vz8FBARopA8dOpQA0KFDh7T2qaqqUv9fXxe6srKSPD09yc7OjjIzMzXShwwZQgAoKSlJq5zw8HBj3raGurrQpaWlNHPmTLKwsND4fEaNGkVFRUW1lv3w4UMqKChQvzIzMxve5fLz+2sUWSYTfjbEvn2N03Wu+dq3T92ExYv/GmU29iWXC/sbIjW1YU3+80/D6mk2g1gzZszA+PHjNdLc3NwAALbVvopKSkpQWloKIkL//v3x2WefobCwEPb29rh79y4OHjyIwYMHIzQ0VKsOqbTu3v2FCxdw8+ZNzJs3Dx4eHup0mUyGhQsX4uDBg0hISEDfvn019ptvrkkJACwsLODt7Y3Bgwdj1KhRcHFxwenTp7F69Wq88sor+P7776FQKHTuu3z5cixevNh0jSkvF468KlVVws+q2/3q2rcpVKu3qAiQSOpXjEQi7G+IhnbAmvKacL0C2NfXFwMHDtS5LScnB++99x727t2LnJwcre3379+Hvb09UlNTQUTo2rVrfZoAAEhPTwcABOroEnbp0kUjT3V+fn71rrMu4eHhOH36NK5cuQLrR8Ogo0aNQmBgICIiIrBhwwZERkbq3HfBggUaXy6FhYVo165d/RtjaSl0m9PTheCVyQBv77qDV7VvU6hWr51d/U/BiQwPLDc3wMvL+G60RAJ06tS0o9AmHcRSKpUIDQ1FfHw8wsPDsWvXLhw8eBCHDx9Wn/sqH41IkAkGR+pbhrkGkv744w989dVXGD58uDp4VV599VVIpVIkJibq3V+hUMDe3l7j1WB79ghBCwj/fvedYft5eja87vrw8lL/t0cPjTEto1RWAkFBhuWVSISrV/U52uv5Lm40Jp1KeenSJfz22294//33tbqCmzZt0vjZz88PEokEFy9erLNciZ5P1sfHBwBw5coVrW2XL1/WyNMYbt26BQCoqKjQ2lZZWQmlUonK+v5F1ldgIJCSYli3ubqAAOFm/MYcyFIohPmKjwweDHh41G8gq1074YYnQ02dCixaJFx3NuS4IJUKMz8nTTK+baZk0iOw7NFk1JpHxsuXL2tdRnJ2dsbQoUNx6NAh9SWj6qqXoRolvnfvnkaebt26wdPTE/Hx8ergAYSj/PLlywEAo0ePbsA7Mk6HDh0gk8mwd+9erQkfcXFxAICePXs2Wns0GNsllsuFJ2nU9yTUWBIJ0LWrxuN3ZDJg7lwhWIwhlQr7GTM32tkZ+PZboRl11SeVCq/vvmviSRww8RE4ICAAgYGB+Oijj1BSUoIOHTogJSUFGzZsQOfOnXHhwgWN/GvXrkXfvn0xbNgw9WWk0tJSJCcnw8vLCx9++CEAoFevXli7di3mzJmDoUOHwsLCAv3794ebmxvWrVuHkSNHIigoCLNnz4aTkxMSEhJw7NgxzJo1S2sAqz62bt2KmzdvAgByc3NRXl6OpUuXAgAcHR0xd+5cAMKXUlRUFFasWIFnn30Ws2bNgouLC06dOoVt27bB09MTc+bMaXB7Gs3kycJdA41ZXw2RkcDu3cIMK0M6L3I50K2b1oQugwwdCuzfL8zIKikR0qofi1TfZS1aAHv3Av36GV+HyRkzZG3ITKyMjAwKCwsjV1dXsra2pqCgIEpISKDo6GgCQDdu3NDIn5WVRbNnz6Z27dqRhYUFubm5UWhoKB05ckSdp7KykqKiosjd3Z2kUikBmvOuT5w4QYMGDSJ7e3tSKBTUuXNnWrVqFSmVSo266jujq1+/fnovn9W8pKRUKmnXrl30/PPPU6tWrcjCwoKeeuopmj17Nv1p6PWGR5p8JpZI50Ln5TX8ba9aReTrq1l+x45E69cTFRbWr9wmv5mBNa4mD2CiZnU30ocf6r8bqV07YXtD70aqTqkkunuX6OZNIahrHA+MZo7fp4TIBMPBzCwKCwvh4OCAgoIC04xI168RwsBSdrb57gd2dxdmfhl4P/ChQ8C5c3/dDxwUJAxYmep+YHMxx++Tb+hntbO3F+5uN2ZI1xhKJbBli0HBCwhBOnSoYVMknwT8SB1Wt9BQ4ZlY5hAbK5TP6oUDmBkmKuqvIDb2uk5Nqv1jY4VyWb1xADPDRUUJJ6CtWtX/hFMqFfY/dIiD1wQ4gJlxQkOFW/5efx2wshIujtY12UOVx8oKmDNH2J+7zSbBAcyMZ28v3Ml+5w6wapUwDKznDisoFML2VauE/KtXGzxgxerGl5GasWZxGclQlZXCkTUj4695115ewiUoXp0QAF9GYs2ZXC48SePR0zRY4+AuNGMixgHMmIhxADMmYhzAjIkYBzBjIsYBzJiIcQAzJmIcwIyJGAcwYyLGAcyYiHEAMyZiHMCMiRgHMGMixgHMmIhxADMmYhzAjIkYBzBjIsYBzJiI8SN1mDjk5wvrqZw/D1y6BBQUCKs62NkJj/Hp3h3o2VNYpuUJwgHMmi8i4Ngx4NNPgT17hIWRZDIhXbVOk0QCfP+98FA9iURYAmbuXGHtlea+WJIJcBeaNU8pKUDfvkBIiLCSdlWVkF5VpbnIGtFfCwcTAUeOACNGAIGBwhH7MccBzJoXImHJlS5d/lpc3JCVvVVUgZ6aCvTuDSxcaNz+IsMBzJoPpRJ47TVg/nzh2dINCTzVkfqDD4BXXgHKykzXzmaEA5g1D0TCuevGjaYv94cfgHHjHssjMQcwax7i44H1681TtlIpDHR98IF5ym9CHMCs6d26BcybV/ciaQ1BBCxeDPz2m/nqaAJ8GekJoFQKp4QWFk3dEj3mzQMePhSCzJyIgGnThMExA74siouBixeBoiJhPbZnnwWsrc3bRGPxEfgxVVwMfP450LWrELiWloCjoxAr//ufaesiAo4eFVYOHTtWiJGtW4WYrNONG8JlosY4P62qAi5cAE6dqjVbWpqwdLG7O/Dcc8Il5eBgoHVr4K23gD/+MH9TDUZGSExMJAB6X4mJicYUZ7A9e/ZQdHS0Wcquy/Xr12nRokXUq1cvcnV1pRYtWlDXrl1p6dKl9ODBA5373Llzh6ZNm0Zubm6kUCioS5cutHHjRqPrLigoIABUUFBg1H7nzxO1bEkkkQgvIcSEl1wu/PuvfxFVVRndJC2nTxP5+/9VtlT6Vx0ODkSffVZHAW+/TSSTaTbSnC+5nGjCBL3NOXiQyNr6r/dQ8yWTEdnZER09avxnVd/fZ23qFcDjxo2jrVu3ar3u3LljsoZVFxERQUZ+15jMv/71L7K1taXx48fTqlWraP369TR27FgCQE8//TSVlJRo5L937x75+vqStbU1LViwgDZu3EgvvfQSAaCYmBij6q7PL/zyZaIWLQyLibffNqo5Wo4dI7K0FIK2tnqWLdNTgFJJ5ObWeMFbPYhr/N6IiM6cMez9SKVCkF+8aNzn1WwCePny5SZrgCHMGcAPHz6kiooKvdvPnTtH9+7d00p/9913CQCtXbtWI/2dd94hALR7926N9BEjRpCFhQWlp6cb3Lb6/ML79jXugHbhgsFFaygtJXJxqfuPXfU6c0ZHIVlZjR+8tTSoVy/D349MRjRwoHGfmTkC2OTnwGfPnsXUqVPh7+8PGxsb2NnZITg4GHv27NGZ/86dO4iMjIS3tzcUCgXc3NwQGhqKw4cPAwC8vLwQHx8PAJBIJOrX0aNH1WUkJSVhyJAhcHR0hLW1Nbp27Yo1a9aAiDTqmjp1KiQSCXJzczF9+nS0atUK1tbWyMrK0vt+evToAUdHR630sWPHAgAuXbqkkb59+3a0b98eY8aM0UifP38+KioqsGvXLr11NdSlS8LpnWoyUl3k8vpfufn6a+H+guqzGmurZ80aHRt+/rl+lTeURKJV98WLQHKyYe8HED7jI0eA3383ffOMUa9R6JKSEuTl5WmkKRQK2NnZYc+ePUhJScGECRPg4eGB/Px8xMfHY8yYMdi+fTsmTpyo3icjIwPBwcHIzs5GREQEunfvjuLiYpw5cwZHjhxBaGgoVq5ciU8++QQnTpzA1q1b1fsGBAQAAA4cOICRI0fC1dUVUVFRcHJywu7duxEZGYlLly5ho46JAaGhoWjTpg0WLVqE4uJitGjRwujP4NatWwAANzc3ddqdO3eQmZmp8R5V+vTpA4lEgrNnzxpdl6G2bxeCxdDxoMpKYbDps88AqZFf5Zs3C/sY8gdfWSkE/IYNgK1ttQ3Xrws3HBj6jWMqcjlw9apG0o4dxn12gND0nTuBRYtM3D5jGHO4rm0Qa+TIkUREOgd2iouLyd/fnwICAjTShw4dSgDo0KFDWvtUVRth0deFrqysJE9PT7Kzs6PMzEyN9CFDhhAASkpK0ionPDzcmLets97evXuTXC6na9euqdPPnz9PAOhtPSeXLVu2pKCgIL3lPnz4kAoKCtSvzMxMo7pcU6bUbzzo7l3j3j8Rkaen8fXcuFGjkOhoIguLpjkHnjFDoynh4cZ/dhYWRK+/bvhnZo4udL2OwDNmzMD48eM10lRHIttqX7ElJSUoLS0FEaF///747LPPUFhYCHt7e9y9excHDx7E4MGDERoaqlWH1IBDwoULF3Dz5k3MmzcPHh4e6nSZTIaFCxfi4MGDSEhIQN++fTX2mz9/vlHvt6bIyEicOXMGS5cuRYcOHdTpJSUlAITeiC5WVlbqPLosX74cixcvrne7LCzqNxfC0rKJ9lHdGtgUatxq2JifnSnVK4B9fX0xcOBAndtycnLw3nvvYe/evcjJydHafv/+fdjb2yM1NRVEhK5du9anCQCA9PR0AEBgYKDWti5dumjkqc7Pz6/edb733ntYt24dZs6ciYULF2pss7GxAQCU6Zk4X1paCvdabjhfsGCBxpdLYWEh2rVrZ3DbunQxvjfarh3wqNlGCQoSLuEa2uVs2RKodrYhcHJq/O4zIERqjXGNTp0MP/9VqawEHp3JNRmTDmIplUqEhoYiPj4e4eHh2LVrFw4ePIjDhw+rzwuVjz4lMsE3b33LsKnPXyyAmJgYLFu2DOHh4diwYQMkNb6y27ZtCwA6B8UePnyI/Px8jZ5CTQqFAvb29hovY4SHGzfbSioV7h+oz5Hn9dcND16pVMgvr3m46Nq1aY7AFRXAM89oJIWHG3//v5UVMGGC6ZpVHyYN4EuXLuG3337DO++8g48//hhjx47F4MGDMXDgQFTV+Kb18/ODRCLBxYsX6yy3ZqCo+Pj4AACuXLmite3y5csaeRpq8eLFWLx4MSZPnoy4uDidXXx3d3d4eHjg9OnTWtvOnDkDIkJQUJBJ2qOLs7MwC8qQASmpVJgWOG1a/erq2xd44YW6/+hlMmEa4muv6dj47LPmnf9cm+7dNX50dQUmTjQ8iGUyYMYM4Yk+TcqYE+a6rgNfunSJANCiRYu00i0tLQkA3ag2kjFs2DC9g1hKpVL9/7lz5xIAultjtEU1iGVvb09ZWVnq9KqqKnXZugaxjLV48WICQJMmTaLKyspa87799tsE6L4OLJfLKS0tzeB66zPoUVJC1Lt37dczpVJhAObwYYOL1Skvj6hLF/2DP3I5kb29MFtLr2eeMfziq6lerq46p6Hdu0fUoUPdg1kyGVHXrkRFRcZ9Xs1+IkdlZSUFBgaSQqGgN998kzZu3Ej//Oc/yc7Ojrp166YVwOnp6eTu7k5yuZxmzJhB69atoxUrVtDYsWM1RnK3bt1KAGjChAn05Zdf0o4dOyg7O5uIiPbv309yuZxat25NMTExtGrVKurXrx8BoFmzZmm0rz4BvHbtWgJATz31FG3ZskVr9lnNL5+7d++St7c32djY0MKFC+nzzz+n4cOH6/xiq0t9f+HFxUTTp/81tVEVH6rpgR07Ep08aVSRehUVCYPJLVtq/pFbWhJNm0b0++91FLBxY+MGr0xG9P77epuTm0v0/POan1f1LyRAmMChY25PnZp9ABMRZWRkUFhYGLm6upK1tTUFBQVRQkICRUdHawUwEVFWVhbNnj2b2rVrRxYWFuTm5kahoaF05MgRdZ7KykqKiooid3d3kkqlBGjOuz5x4gQNGjSI7O3tSaFQUOfOnWnVqlUaR3Gi+gWwah99r379+mntc/v2bZo6dSq1bNmSFAoFBQYG0vr1642ql6jhv/A7d4iWLycaN45o1Cii114jOn5cmMFoauXlRImJRN98Q3TgAFF+voE7PnhAZGvbeAEslRJVu+Soi1JJdOoU0cSJwheTjY0w43PqVKKzZ+v/GZkjgCVETTWOz+pSWFgIBwcHFBQUGD2gJSrLlgmzIcz9pyiTATNnCjNXmoA5fp98OyFrem+/LVwDM+djYKVSoFUr4OOPzVdHE+AAZk3PwgLYtk0IYHONShMBX37ZDIaNTYsDmDUPXboACQnCkdIcQfz558CAAaYvt4lxALPm46WXgH37hPmJWrM+6kEmE74QtmwRLto+hjiAWfMydCjw669aEy2MJpEAvr7CPYIREaZpWzPEAcyanw4dgKQkYOVKYRI1YNgAl2oKmr098P77whdBjx5ma2ZzwAHMmieZDPjHP4RHzn7zjbBGUm1z2K2sgD59gLg44M4dICYG0HNX2OOEHyvLmjcLCyAsTHgplcKaR5cuAYWFwsiyra2wvGjHjk/EaoQ1cQAz8ZBKAX9/4cUAcBeaMVHjAGZMxDiAGRMxDmDGRIwDmDER4wBmTMQ4gBkTMQ5gxkSMJ3I0Y6qHpRQWFjZxS5gpqH6PpnwIDgdwM1ZUVAQARj3cnTV/RUVFcHBwMElZ/EysZkypVOL27duws7PT+2xsQ6hWeMjMzDT7s7Uaqy4x1kNEKCoqQps2bQxaOsgQfARuxqRSaa0rORirPqs9NPe6xFaPqY68KjyIxZiIcQAzJmIcwE8AhUKB6OhovcueirGux62e+uJBLMZEjI/AjIkYBzBjIsYBzJiIcQAzJmIcwI+Bo0ePQiKR4IMPPqg13/Lly/Hqq6/C29sbEokEXl5eZqkrJSUF77//Pnr37o2WLVvCzs4OzzzzDJYtW4bi4mKT1QMAEolE7+v+/fsmqwcACgoKsGDBAnTo0AFWVlZwdnZG3759sWfPHoPekznwTKwnyMKFC+Hs7Ixu3boZ9MddX1988QXWrl2LESNGYOLEibC0tERiYiLee+89fP311zhz5gysra1NVt/zzz+Pv/3tb1rptra2JqsjMzMTISEhuHv3LqZNm4ZOnTqhpKQE165dwx9//GGyeozFAfwESUtLg7e3NwCgc+fOePDggVnqCQsLwzvvvANHR0d12muvvQY/Pz8sW7YMX3zxBebMmWOy+ry9vTF58mSTlafLlClTUFxcjF9//bVZ3VzCXegniCp4za1Hjx4awasyduxYAMClS5dMXmd5ebn67i1TO3HiBI4dO4Z//etfaNeuHSorKw0+FTA3DmDWaG7dugUAcHNzM2m53377LWxsbGBvbw8XFxfMnDkTd+7cMVn5Bw4cACB8AY4ZMwbW1tZo0aIFvLy8sHbtWpPVUx/chWaNoqqqCkuWLIFcLsekSZNMVm5QUBDCwsLg5+eHkpISJCYmIi4uDocOHUJycjJat27d4DquXbsGAJg5cybat2+PzZs3QyKRYN26dZg3bx7u3buHRYsWNbie+uAAZo0iMjISZ86cwdKlS9GhQweTlXv27FmNnydNmoR+/fohPDwc0dHR2LhxY4PrUHXNbW1tcfz4cfW86HHjxqFTp05Yvnw55s6dCycnpwbXZSzuQjOze++997Bu3TrMnDkTCxcuNHt9U6ZMgZeXF/bv32+S8lQj5hMnTtS4qcHS0hKTJk1CaWkpkpOTTVKXsTiAmVnFxMRg2bJlCA8Px4YNGxr0ZBFjeHl5ITc31yRlqR6qoKs7rkq7e/euSeoyFgcwM5vFixdj8eLFmDx5MuLi4kz2GJm6EBFSU1Ph7u5ukvJ69+4NQLgWXJPqGnCrVq1MUpexOICZWSxZsgQxMTGYNGkStmzZYpbgzc7O1pm+Zs0aZGVl4eWXXzZJPSNHjoS9vT2+/PJLFBQUqNOLiooQHx8PJycn9OnTxyR1GYsHsR4jiYmJqKys1Ep3dHTE3LlzsXXrVty8eRMAkJubi/LycixdulQjjynqkkgkiI6OxlNPPYXQ0FDs2LFDI0+rVq0QGhra4HpSU1Nx5MgRDB8+HJ6enigtLcXRo0exb98++Pn5ISYmxiTvZ+7cuYiNjcWMGTPQs2dPzJw5ExKJBJs3b8aff/6JLVu2wMbGxuC6TIqY6CUmJhIAvS9PT08iIurXr1+deUxRV0RERK15+vXrZ5J69u7dS4MHD6a2bduSQqEgKysrCgwMpHfffZfu379v0s+OiOj777+n4OBgsrW1JRsbG3r++efpwIEDBtVjLvxEDsZEjM+BGRMxDmDGRIwDmDER4wBmTMQ4gBkTMQ5gxkSMA5gxEeMAZkzEOIDZE0H19EljpleKAQfwYyI8PBwSiQTu7u465/Q2RzExMZBIJDh69KhJypNIJHjxxRdNUpZY8M0Mj4HCwkLs3r0bEokE2dnZ2L9/P0aOHNnUzWpWevbsiatXr8LV1bWpm2JSfAR+DOzYsQMlJSV488031XfJME02Njbo2LHjYxfAfDfSYyAoKIgsLS0pPz+fnnvuOZLJZHT79m2tfHh0J1BOTg5NmzaNWrZsSVZWVtSrVy9KTEzUyq+6e6miooKWLFlCXl5eZGlpSX5+fvTpp5/qbEtxcTFFR0dThw4dSKFQkJOTEw0bNoySkpJ0ll3zVf3un59++ommTZtG/v7+ZGtrS7a2ttS9e3fasGGDRlm13VEUFxenkSc6OlqrzZcvX6axY8dSy5YtydLSkry8vCgqKory8/O18np6epKnpyc9ePCA3njjDWrTpg1ZWlpSly5d6JtvvtH5mZgTd6FF7tKlSzh37hxGjx4NZ2dnhIeH4+TJk4iPj8c777yjlf/+/fsIDg6Gvb09Jk2ahJycHOzatQuDBw/Gzz//jM6dO2vtM2HCBCQnJ2Po0KGQyWT4+uuvMWfOHFhYWGDWrFnqfGVlZRgwYADOnDmDbt26ISoqSl3+oUOHsGvXLowZMwYAMHXqVADAsWPHEBERoV7mpfrzpD/88EOkpqaid+/eGD16NO7fv4+DBw9i9uzZuH79OlasWAFAeHxOdHQ0Fi9eDE9PT3XZAPDMM8/U+vmdOnUKgwYNQllZGcLCwuDl5YUzZ85g5cqV2L9/P06fPg0XFxeNfSoqKjBo0CDcvXsXY8aMQUlJCXbu3ImxY8fi4MGDGDRoUK11mlSjf2Uwk/rHP/5BACghIYGIiO7fv09WVlbk5+enlRePjkqvv/46VVVVqdM3bdpEAGj27Nka+VVHyV69elFBQYE6/dq1aySXy6lDhw4a+ZcsWUIAaNKkSaRUKtXpv/76q/poXFhYqE6Pjo4mADqP/kRE6enpWmkVFRUUGhpKMpmMbt68qfX+9N1rrOsIXFVVRX5+fgSADh48qJF/wYIFBIBmzJihke7p6UkAaOTIkVRWVqZOP3LkCAGgwYMH66zfXDiARaysrIxcXFzIyclJ449p3LhxBICOHTumkR8A2draUlFRkUZ6RUUFyeVy6tatm0a6KoB/+uknrbpV26oHpLe3N1lYWFBmZqZW/tmzZxMA2rp1qzqtrgDWZ/fu3QSAtmzZovX+jAng48ePEwAaOnSoVv4HDx6Qi4sLWVtba3y2qgDW9eXi6elJzs7ORr2XhuJBLBH77rvvkJ+fj3HjxsHS0lKdHh4eDkBYZKwmPz8/tGjRQiNNLpejVatWehc869atm1aa6kmNqn0KCwuRnp4OX19f9bbqVJd3Ll68WNfbUisqKkJ0dDS6du2KFi1aqFcdfOWVVwAAt2/fNrgsXX755ReNtlVna2uLHj16oLS0FCkpKRrbHB0d0b59e619PDw8zLponC58DixiqgCdMmWKRvrgwYPh7u6Ob775BqtXr4a9vb16m4ODg86y5HI5qqqqdG7TtY9cLvzpqPYpLCwEoP/pjKonRFZ/KFxtysvL8eKLL+LChQt49tlnMWXKFLi4uEAulyMjIwPx8fEoKyszqCx96tvm2j5DpVLZoDYZiwNYpDIzM3H48GEAQHBwsN58O3fu1Ln0pqmpviT0PSlSlV79y6Q2e/fuxYULFzBz5kx8/vnnGtt27tyJ+Pj4BrQWGm0xVZubAgewSMXFxUGpVOK5557TuVRJeXk5tm7dis2bNzdaAHt7eyM1NRW3bt1C27ZtNbYfO3YMgOaosEwmAwCdR/60tDQA0Plo2BMnTuhsg1Qq1duL0OXZZ58FIEyzfPvttzW2lZSU4Pz587C2tjbpUjCmxgEsQkSEuLg4SCQSfPnllzrPxwDg8uXLOHv2LC5fvqzz8pCpRUREIDo6GgsWLEB8fLx6FYbLly8jLi4ODg4OGDVqlDq/s7MzACArK0urLE9PTwDAyZMnMWLECHX6sWPHtI7I1cvTVZY+wcHB8PHxwX//+18cOXIEAwcOVG9bvnw58vLyMH36dI3xheaGA1iE/u///g8ZGRkICQnRG7wAMG3aNPzyyy/YvHkzYmNjzd6ut99+G/v378fWrVtx9epVDBgwALm5udi1axcqKirw5Zdfws7OTp0/JCQEEokE7777Lq5duwYHBwc4ODjg73//O0aMGAEvLy989NFH6i+g69ev44cffsCoUaOwe/durfr79++Pr7/+GmFhYXj22Wchk8nw0ksvoUuXLjrbK5VKsWXLFgwePBjDhg3Dq6++Ck9PTyQnJ+Onn36Cj48PPvjgA7N9XibRqGPezCTGjx+vdUlGl7y8PLK0tCRXV1cqKyur9TKLaoZRdapLRbqonv1848YNjfQHDx7QokWLyN/fnywtLcnR0ZGGDh1KJ06c0FnOli1bqEuXLqRQKLRmYqWnp9Mrr7xCLVu2JBsbGwoKCqKdO3fqnVX1559/0tixY8nV1ZWkUqnBM7F+++03CgsLI1dXV7KwsCBPT0+KjIyk3Nxcgz4nQz4vc+HnQjMmYnwdmDER4wBmTMQ4gBkTMQ5gxkSMA5gxEeMAZkzEOIAZEzEOYMZEjAOYMRHjAGZMxDiAGRMxDmDGRIwDmDER4wBmTMQ4gBkTMQ5gxkSMA5gxEeMAZkzEOIAZEzEOYMZEjAOYMRHjAGZMxDiAmV4SiURjsWxDZWRkQCKRICYmxuRtYpo4gEXg6NGj6qU19b2Y+Xz33XfN9suIl1bRRakE0tOBwkLA3h7w9gakTf9dN27cOAwfPrzR6istLVUvQPYk++677xAfH98sg5gDuLqiImDTJmDNGuDGjb/Svb2BefOAGTOAamv7NLZnnnkGkydPNmsdZWVlkMlkkMvlsLKyMmtdrOGa/rDSXGRmAt27A2++CWRkaG67cQOYP1/YnpnZJM0zlOq89eTJk3j++edhY2MDV1dXzJw5Ew8ePNDIO3XqVEgkEuTm5mL69Olo1aoVrK2t1Sv86TsHTkxMxEsvvQQXFxdYWVnB29sbM2bMQF5enlbe7777Dt27d4eVlRVat26Nt956C5WVlRp5XnzxRXh5eSEjIwOjR4+Go6MjnJycMHXqVDx48ABKpRL//ve/0b59eygUCjz77LM6lxglIqxfvx7du3eHjY0N7OzsEBISgsTERI181c/R62qfl5eXei3i6qcsR48eNej3YW58BAaEI++AAUKg6loqSpV244aQ7+efm+RIXFJSojNILC0tNRahvnjxIkaOHInp06dj8uTJOHr0KDZv3gypVIqNGzdq7R8aGoo2bdpg0aJFKC4uRosWLfS2YcOGDfj73/+Odu3a4fXXX8dTTz2FP/74A/v27UNWVhZcXV3VeQ8cOIB169bhtddew8yZM7F371785z//gZOTExYuXKhRbnFxMUJCQhASEoIPPvgAP//8MzZt2oTS0lK4urri7NmzmDdvHioqKvCf//wHL7/8Mm7evKnxvqdMmYIdO3YgLCwM06ZNQ1lZGbZv347Q0FAkJCRorTVsSPtWrlyJTz75BCdOnMDWrVvV+wYEBOj9jBpVoy6l1lzFxhJJJERCqNb+kkiIVq5s1OapVtbT9xowYIA6LwCSSCR0+vRpjTKGDRtGcrmcioqK1GmqFQbDw8N11guAIiIi1D9nZmaSpaUlderUiQoKCrTyV1VVERHRjRs3CADZ2NhorF6oVCopMDCQ3N3dNfZTreq3YsUKjfRXXnmFJBIJBQUFUUVFhTp97969BIDWr1+vTtu9ezcBoM8++0yjjIqKCurevTt5eXmRUqmsV/tUn1NzxF1opRJYvdq4fVavFvZrZDNmzMDhw4e1Xh9//LFGvj59+qB3794aaf3790dlZSUyap4eAJg/f75B9X/zzTcoLy/HokWLNI58KtIaA32jRo2Cl5eX+meJRIKQkBDcuXNHqzsvk8kwZ84cjbTg4GAQEWbPng25/K/O4vPPPw8ASE1NVadt374dtra2GDVqFPLy8tSv+/fvY8SIEcjIyMDvv/9e7/Y1V9yFTk/XHLCqC5GwT3o64Otrvnbp4Ovrq7GKvD7e3t5aaS4uLgCA/Px8rW1+fn4G1a8KgK5duxqUv652VO+qt27dGgqFQiOvk5MTAGgEWfX06u/l6tWrKC4uhru7u972ZGdnw9/fv17ta644gAsLG3e/RlDbpR/ScY5vY2NjULm69jVVO2rLq29b9TKICM7Ozti1a5fecjp37lzv9jVXHMA6uoJm3U/EOnToAEAYJGs2gziP+Pv74/r16wgKCoKDg4NJy27OE2X4HNjbG2jfHjD0lySRCPvo6H497sLCwmBpaYmlS5eiUEcPpCmPWlOmTAERYcGCBTrbkZ2dXe+yVV3pe/fu1bsMc+EjsFQKREYK13kNFRnZJDOzLl68iG3btunc9vLLL+scWDIlDw8PrFy5EnPmzEGXLl0QHh4OT09P3Lp1C3v37sUXX3yBZ555xqxt0Ed16Wj9+vW4ePEiRowYAVdXV2RlZeH06dNITU1Fenp6vcru1asX1q5dizlz5mDo0KGwsLBA//794ebmZuJ3YTwOYECYYbVunTCYVWOSgQa5XDjyTp/eeG2rZteuXXrP8a5evWr2AAaAv//97/Dx8cHHH3+M1atXo6ysDG3atMGAAQPQrl07s9dfmy+++AIhISHYuHEjli9fjvLycri7u6Nbt25Yvnx5vcudMGECfv75Z+zcuRO7du2CUqlEYmJiswhgCYnlbN3cMjOFSRqqSxPVPxZV99rPDzhyBGjiP1TGVPgcWKVdO2GG1SefADUuW6B9eyA2Fjh/noOXNSt8BNalmd6NxFhNHMCMiRgfVhgTMQ5gxkSMA5gxEeMAZkzEOIAZEzEOYMZEjAOYMRHjAGZMxDiAGRMxDmBmNNXjaM1hy5Ytzeqxrc0dB7AI1FxaRSaTwdHREQEBAZg4cSISEhJQVVXV1M002MWLFxETE6PzAXvMOHw/sA7N9V4G1dIqRIQHDx7g999/xw8//IAdO3YgKCgICQkJ8PDwaOpm1unixYtYvHix+oHu1U2ZMgXjx4+HpaVl0zROZDiAq2nmK6voXFplxYoV+Oijj/DOO+/gpZdews8//6zxCFaxkclkvB6TMRr9SdTN1B9/EPn5Cc9tr/mMd1Wan5+Qr7GpHuy+fPlyvXnGjx9PAGjbtm3qtIcPH9KyZcuoU6dOpFAoyMHBgYYPH04XLlzQWX5cXBytXr2a/Pz8SKFQkK+vL63U8RB7fQ86v337Nr322mvUrl07srCwoNatW9OsWbMoOztba9+ar+joaCIiiouLIwCUmJioUXZeXh7NmzdPo+wZM2bQ7du39b6Xzz//nAICAsjS0pKeeuop+vDDD/V+fmIl3q9qExLJyiq1+tvf/oadO3fihx9+wKRJk1BRUYEhQ4bg1KlTmDJlCubOnYuCggJs2rQJwcHBOH78OHr06KFRxpo1a3Dnzh3Mnj0bdnZ22LFjB6KiopCfn48lS5bUWv8ff/yBPn36oLy8HDNmzICPjw/S0tKwbt06JCYm4vz583BwcMDs2bOhUCiwceNGLFy4UP10y6efflpv2YWFhXjuuedw/fp1REREoGfPnrh8+TI2bNiAQ4cO4dy5c2jVqpXGPuvXr0dOTg5mzpwJBwcHbNu2Df/617/g4eGBiRMn1vNTboaa+hukOWjmK6sYdATOz88nANStWzciIlqxYgUBoP/+978a+QoKCqhdu3bUr18/rfJbtGhBmZmZ6vSysjIKCgoimUxGGRkZ6nRdR+ARI0aQq6urxv5EROfOnSOZTKY+whLpP8rq2/buu+8SAK3ewLZt2wgAzZo1S+u9tG7dmu7du6dOLy4uJldXV+rdu7dWnWLWDIZmmpaIVlapleqBdqrHvW7fvh1+fn7o0aOHxlIj5eXlCA0NxcmTJ1FaWqpRxqRJkzQGwSwtLfHGG2+gqqoK+/bt01v3/fv3sX//fgwfPhxWVlYa9Xl5ecHX1xeHDh2q93vbs2cPnJ2d8frrr2ukT5w4Eb6+vtizZ4/WPtOmTYOjo6P6ZxsbG/Tu3VtreRWxe+K70CJaWaVWqsBVBfLVq1dRWlqKli1b6t0nLy9P40mSuh7W3qlTJwBAWlqa3nJSUlKgVCqxZcsWbNmyRWceXcuYGCo9PR3PPPMMLCwsNNIlEgkCAwOxd+9eFBYWajyVU9+yKbqWlhGzJz6AH5eVVS5evAgA6NixIwDhIeudOnXCqlWr9O5TM7h1Tc6gRwMAtU3cUOWZMGECput55K61tbX+xjcA6Xki1JMykv3EB/DjsrKKat3f4cOHAxCWGvnzzz/Rv39/rVUD9fnf//6nlXb16lUAtR9BfX19IZFIUFZWZtDia8bO4vL29kZKSgoqKiq0jsL/+9//4Orq2ijPxG6OnvhzYLGvrEJE+Oijj7Br1y4888wzGDt2LABhQkRubq7W0qMqupYa2b59O7KystQ/l5eXIzY2FjKZDCNGjNDbBhcXFwwbNgx79+5FUlKSzjbm5uaqfzZ2qZLRo0fj7t272LBhg0b6zp07kZqaijFjxhhUzuPoiT8Ci2hlFY2lVR48eIDU1FTs27cPKSkp6NmzJxISEtRdx3/84x84fPgw3nnnHRw9ehQDBgyAvb09/vjjD/zf//0frKyskJiYqFG+v78/evXqhddeew12dnb46quvcO7cOSxatAienp61tm39+vV47rnnEBISgilTpqBbt25QKpVIT0/H3r17ER4ejpiYGABAjx49IJVKsXz5cty7dw82Njbo3Lmz1uqBKm+//Ta+/fZbREZG4pdffkFQUJD6MpKHh0edl7gea005BN5cFBYKkzTk8tovIcnlRP7+Qv7GpLo0onpJpVKyt7enDh060IQJE2j37t1UWVmptV9FRQWtWrWKevToQTY2NmRjY0O+vr40ceJE+vHHH7XKj4uLo1WrVpGvry9ZWloaPZEjNzeX/vnPf6ongjg4OFDnzp0pMjKSrly5opF38+bN5O/vT3K53OCJHHPnziUPDw+ysLAgd3d3mjFjBt26dUvnZxUXF2dwu8WMnwv9yJO8ssrRo0cREhKCuLg4TJ06tambw4zwxJ8Dq/DKKkyMnvhz4Ors7ICoKOEctznejcRYTRzAOkilzWuSBmP68DkwYyLGHUPGRIwDmDER4wBmTMQ4gBkTMQ5gxkSMA5gxEeMAZkzEOIAZEzEOYMZEjAOYMRHjAGZMxDiAGRMxDmDGRIwDmDER4wBmTMQ4gBkTMQ5gxkSMA5gxEeMAZkzEOIAbwYsvvgivms+qZcwEOIAbICcnB2+//TY6d+4MOzs7ODg4wM/PD+PHj0dCQkJTN489AfiplPWUmZmJoKAgFBUVYdKkSejatSsAIDU1Ffv374e/vz9++OEHAMIiYUQEhULRlE1mjyEO4HqKjIzEmjVr8P333+tcuS8rK0tjtXvGzIG70PWUkpICAAgJCdG5vXrw1jwH3rJlCyQSid5X9VXuiQjr169H9+7dYWNjAzs7O4SEhGitLMieTLwyQz2pFrz+/PPPERUVZdSi1S+88AK2bt2qlf7hhx/i8uXLaNWqlTptypQp2LFjB8LCwjBt2jSUlZVh+/btCA0NRUJCAl5++eWGvxkmXk21LKLYpaWlkb29PQGgdu3a0cSJEyk2NpbOnz+vlbdfv37k6elZa3krVqwgAPTmm2+q03bv3k0A6LPPPtPIW1FRQd27dycvLy9SKpUmeT9MnDiAG+DGjRv0+uuvU9u2bTXW7+3SpYtGINcVwLt37yapVEpjxozRCMgxY8aQra0t3blzh3JzczVeMTExBICuX79uzrfImjkOYBO5ffs27d69m15++WUCQO7u7pSfn09EtQfw6dOnydramnr16kUlJSUa2wICAjS+GHS9jh8/bu63xpoxHoU2g4kTJ2LHjh3YunUrJk+ejBdffBEZGRnIyMjQyJeWloY+ffqgRYsWOHPmDNzc3DS2BwQEICcnB7t27dJbV/fu3eHk5GSOt8FEgAexzKBPnz7YsWMHbt26pTdPfn4+hg0bhsrKShw4cEAreAHA398f169fR1BQEBwcHMzZZCZSfBmpnhITE1FaWqqVrlQqsW/fPgBAp06ddO5bVlaGUaNGISMjA3v27EHHjh115psyZQqICAsWLICujlJ2dnYD3gF7HPARuJ5WrFiBpKQkDB8+HN27d4eDgwPu3LmD3bt34+eff0ZISAheeuklnfu+//77OHnyJEaNGoXMzExs27ZNY3vfvn3h7e2tvnS0fv16XLx4ESNGjICrqyuysrJw+vRppKamIj09vTHeLmuumvYUXLxOnz5N8+fPpx49epCbmxvJ5XJycHCg3r1704oVK+jhw4fqvDUHsSIiImodmIqLi9Oo68svv6TnnnuO7OzsSKFQkKenJ40ePZp27tzZSO+WNVc8iMWYiPE5MGMixgHMmIhxADMmYhzAjIkYBzBjIsYBzJiIcQAzJmIcwIyJGAcwYyLGAcyYiHEAMyZiHMCMiRgHMGMixgEsAikpKXj//ffRu3dvtGzZEnZ2dnjmmWewbNkyFBcXa+Wv7ZnT9+/fb/w3wMyGbycUgXfeeQdr167FiBEj0KdPH1haWiIxMRFff/01nn76aZw5cwbW1tbq/BKJBM8//zz+9re/aZU1btw4WFhYNGbzmRlxAOuQk5ODpKQk5Ofnw8XFBcHBwTqfWdVYzp8/D19fXzg6Omqkv/fee1i2bBnWrl2LOXPmqNMlEgkiIiI0Vnhgjyd+pE4Np06dQmxsLP7880912ldffYU33ngDffv2bZI29ejRQ2f62LFjsWzZMly6dEnn9vLycpSVlcHOzs6czWNNiM+Bq8nJyUFsbCzy8/MREBCAwMBABAQEID8/H7GxscjJyWnqJmpQPfVSV+/g22+/hY2NDezt7eHi4oKZM2fizp07jd1EZmZ8BK4mKSkJf/75JwICAiCVCt9tUqkUPj4+uHr1KpKSkjB69OgmbqWgqqoKS5YsgVwux6RJkzS2BQUFISwsDH5+figpKUFiYiLi4uJw6NAhJCcno3Xr1k3UamZqHMDV5OfnA4A6eFWkUikkEol6e3MQGRmJM2fOYOnSpejQoYPGtrNnz2r8PGnSJPTr1w/h4eGIjo7Gxo0bG7OpzIy4C12Ni4sLAOHZztUplUoQkXp7U3vvvfewbt06zJw5EwsXLjRonylTpsDLywv79+83c+tYY+IAriY4OBitW7dGWlqaOoiVSiXS0tLQpk0bBAcHN3ELgZiYGCxbtgzh4eHYsGGDUcuaenl5ITc314ytY42NLyPVUH0UWiKRgIjQpk0bREVFNdkotMrixYsRExODyZMnIz4+XqurXxsiwlNPPQWJRII//vjDjK1kjYnPgWvo27cvfH19m9V1YABYsmQJYmJiMGnSJGzZskVv8GZnZ2ssEK6yZs0aZGVlaVwvZuLHR2AR+PTTTzF37lw89dRTWLJkCWQymcb2Vq1aITQ0FAAQFRWFI0eOYPjw4fD09ERpaSmOHj2Kffv2wc/PD6dOnYKrq2tTvA1mBhzAIjB16lTEx8fr3d6vXz8cPXoUAPD9999j3bp1uHz5MvLy8iCRSODj44NRo0bhrbfe4lUOHzMcwIyJGI9CMyZiHMCMiRgHMGMixgHMmIhxADMmYhzAjIkYBzBjIsYBzJiIcQAzJmIcwIyJGAcwYyLGAcyYiHEAMyZiHMB6VFRUoLCwEBUVFU3dFADGL5eyY8cOdO/eHdbW1nB1dcWECRNw8+bNxm84Myt+IkcNqampOHz4MBITE/Hw4UNYWVkhJCQEgwYNgo+PT5O2Td9yKba2tho/r127FvPmzUNwcDBiY2ORl5eHlStX4vjx4zh37hzatGnTWE1mZsb3A1dz7NgxxMbGIjc3F05OTrCyssLDhw9x7949tGzZEvPnz8cLL7zQJG0zdLmU/Px8eHl5wd/fH8nJyZDLhe/o8+fPo2fPnpg+fTo2bdrUCC1mjYG70I+kpqYiNjYWDx48QGBgINq2bQsXFxe0bdsWgYGBKCoqwieffIK0tLQmbWd5eTmKior0bt+7dy8ePHiAyMhIdfACwvIsL7zwAr7++muUl5c3RlNZI+AAfuTw4cPIzc2Ft7e31qNaVY+lyc3NxaFDh5qohYYtl6J6qLuuJ2j27dsXRUVFuHbtWqO0l5kfBzCEAavExEQ4OTnpfc6yRCKBk5MTEhMTm2RgKygoCO+//z6++eYbbNu2DaNHj0ZcXBx69uypsRCbar0kDw8PrTJUaVlZWY3TaGZ2PIgFoLS0VD1gVRsrKyuUlZWhtLS00dfYNXS5lJKSEgCAQqHQKkP1/lR5mPjxERiAtbW1esCqNg8fPoRCodBYTLsp6VouxcbGBgBQVlamlb+0tFQjDxM/DmAAFhYWCAkJwb1796BvUJ6IcO/ePYSEhDSrFe5rLpfStm1bALq7ybV1r5k4cQA/EhoaipYtWyI9PV0riIkIaWlpcHNzw6BBg5qohdqICKmpqXB3d1enBQUFARCWiKnp1KlTaNGiBTp27NhobWTmxQH8iK+vL+bPn48WLVrgypUruHXrFvLz83Hr1i1cuXIF9vb2eOONN5pkMkd2drbOdNVyKS+//LI6beTIkbCxscHq1atRWVmpTj9//jyOHz+OsWPHwtLS0uxtZo2DJ3LUkJaWhkOHDiExMRFlZWVQKBRNPhPL2OVSVq1ahaioKAQHB2PKlCnIy8tDbGwsLCwscP78eXU3mz0GiOlUXl5OBQUFVF5e3tRNob1799LgwYOpbdu2pFAoyMrKigIDA+ndd9+l+/fv69xn27Zt9Oyzz5KVlRU5OzvTuHHjKD09vZFbzsyNj8CMiRifAzMmYhzAjIkYBzBjIsYBzJiIcQAzJmIcwIyJGAcwYyLGAcyYiHEAMyZiHMCMiRgHMGMixgHMmIjxM7F0ICLk5uaqn5PVsmVLvQ+7Y6wpcQBX8/DhQ5w8eRIHDx7EtWvXUFlZCblcjoCAAAwZMgTBwcF1PvjOHGJiYrB48WK92+VyucaTMmv7srl37x4cHR1N2TzWhDiAH7l37x4++ugjnDp1ChKJBC1btoSFhQUqKipw7tw5JCcnIzg4GG+//TacnJwatW1jxoyBr6+vVvpvv/2Gjz/+GCNGjNDaZugyLEzcOIAhHHk/+ugjHD9+HN7e3lpPbXR2dkZJSQmOHz8OAIiOjm7UI/HTTz+Np59+Wit99uzZAIAZM2ZobfP29sbkyZPN3jbWtHgQC0BSUhJOnTqlM3hVbGxs4O3tjaSkJCQlJTVyC7WVlJRg586daNu2LYYMGaIzT13LsDDxe+IDmIjw3//+FxKJpM7nJdvY2EAikeDHH3/U+/jZxvL111+jsLAQ06ZNg0wm09puyDIsTPye+C50bm4url27hpYtWxqU383NDVevXkVubi7c3NzM3Dr9Nm/eDIlEgunTp2ttCwoKQlhYGPz8/FBSUoLExETExcXh0KFDSE5ORuvWrZugxcwcnvgAfvjwISorKw1+WLuFhQWKiorqXMXBnK5fv46TJ09iwIABaN++vdZ2Q5dhYeL3xHehraystC7D1KaiogIWFhZNcjlJZfPmzQCAmTNnGryPrmVYmPg98QHcsmVLdOzYUWN5ktrk5OQgICDA4C63qVVWVuLLL7+Es7MzRo8ebdS+NZdhYeL3xAewRCLB0KFDoVQq61y1r6SkBESEwYMHN9nMrH379iE7OxtTpkzRuQKhPqRjGRYmfk98AANAcHAwgoODkZ6erjeIS0pKkJ6ers7bVFTdZ13XfgHjlmFh4scPdn9ENRMrKSkJEokEbm5u6plYOTk5IKImm4mlcvv2bTz11FPo3r07kpOTdeYxdhkWJm5P/Ci0ipOTE6Kjo5GUlIT//ve/uH79OoqKimBhYYGePXs26VxolS1btqCqqqrWwav+/fvj2rVr2LZtG/Ly8iCRSODj44N3330Xb731FhwcHBqxxczc+AisA9+NxMSCA5gxEeNBLMZEjAOYMRHjAGZMxDiAGRMxDmDGRIwDmDER4wBmTMQ4gBkTMQ5gxkSMA5gxEeObGXQgIhQXF6OsrAwKhQK2trY8F5o1SxzA1ZSUlOD06dP48ccfkZKSol6Zwd/fH4MHD0afPn3qfHIlY42Ju9CPXLlyBXPmzMHixYuRnJwMqVQKGxsbSKVSJCcnY/HixZgzZw6uXLnSJO1bvnw5Xn31VXh7e0MikcDLy6vW/Dt27ED37t1hbW0NV1dXTJgwATdv3mxwXta88N1IEII3JiYG2dnZ8PHxgaWlpVae8vJypKWloVWrVoiJiUFgYGCjtlEikcDZ2RndunXDzz//DHt7e2RkZOjMu3btWsybNw/BwcGYPHky8vLysHLlSigUCpw7dw5t2rSpV17W/DzxAVxSUoI5c+YgLS0NHTt2rPVcl4hw7do1+Pj44NNPP23U7nR6ejq8vb0BAJ07d8aDBw90BnB+fj68vLzg7++P5ORkyOXCWdL58+fRs2dPTJ8+HZs2bTI6L2uenvgu9OnTp5GWlgYfH586B6pUT7dIS0vDmTNnGqmFAlXw1mXv3r148OABIiMj1QEJAD169MALL7yAr7/+GuXl5UbnZc3TEx3ARIQff/wREolEZ7dZF0tLS0gkEhw8eLDJl1fRRfVQ9759+2pt69u3L4qKinDt2jWj87Lm6YkO4OLiYqSkpMDFxcWo/VxcXJCSklLnY2ibwq1btwAAHh4eWttUaVlZWUbnZc3TEx3AZWVl6ktFxpDL5aiqqmrS5VX0UX2p6HpmtOqBfKo8xuRlzdMTHcAKhQJyuRyVlZVG7VdZWQmZTNakT6jURzWwVlZWprWttLRUI48xeVnz9EQHsK2tLfz9/ZGfn2/Ufvn5+fD392+Wf9xt27YFoLvrW7PLbExe1jw90QEskUgwePBgEJHBo63l5eUgIgwZMqRZTq8MCgoCAJw6dUpr26lTp9CiRQt07NjR6LyseXqiAxgA+vTpo740VNeoMhGpLzn17t27kVponJEjR8LGxgarV6/WODU4f/48jh8/jrFjx6pH3I3Jy5qnJ34iB2D8TKzFixejU6dOjdrGrVu3qqc3rlmzBuXl5XjzzTcBAI6Ojpg7d64676pVqxAVFYXg4GBMmTIFeXl5iI2NhYWFBc6fP6/uOhublzU/HMCPXLlyBf/5z3+QlpYGiUQCFxcX9QBXfn4+iAg+Pj546623Gj14AeDFF1/EsWPHdG7z9PTUmpW1fft2rFixAlevXoWNjQ1CQ0OxfPlynQuCG5OXNS8cwNWUlJTgzJkzOHjwIFJSUlBVVQWZTAZ/f38MGTIEffr0gbW1dVM3kzE1DmAdiAglJSXqtZFsbGya5YAVYxzAjInYEz8KzZiYcQAzJmIcwIyJGAcwYyLGAcyYiHEAMyZiHMCMiRgHMGMixgHMmIhxADMmYhzAjIkYBzBjIsYBzJiIcQAzJmIcwIyJGAcwYyLGAcyYiHEAMyZi/w81XyFLTrxzfAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 150x200 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "f, ax = plt.subplots(figsize=(1.5,2))\n",
    "n_topic = diff.shape[1]\n",
    "n_layer = diff.shape[0]\n",
    "\n",
    "for j in range(n_layer):\n",
    "    y = diff[j, :]\n",
    "    if j==0:\n",
    "        plt.scatter((np.ones(n_topic)*j)[y>0], np.arange(n_topic)[y>0], c=\"red\", s=y[y>0]*5, label=\"Enrichment\")\n",
    "        plt.scatter((np.ones(n_topic)*j)[y<0], np.arange(n_topic)[y<0], c=\"blue\", s=-y[y<0]*5, label=\"Depletion\")\n",
    "    else:\n",
    "        plt.scatter((np.ones(n_topic)*j)[y>0], np.arange(n_topic)[y>0], c=\"red\", s=y[y>0]*5)\n",
    "        plt.scatter((np.ones(n_topic)*j)[y<0], np.arange(n_topic)[y<0], c=\"blue\", s=-y[y<0]*5)\n",
    "    \n",
    "plt.yticks(np.arange(diff.shape[1]), [\"Factor 17\",\"Factor 18\",\"Factor 20\"], rotation=0, fontsize=13)\n",
    "plt.xticks(np.arange(diff.shape[0]), [\"L1\",\"L2\",\"L5\",\"L6\"], rotation=0, fontsize=13)\n",
    "plt.xlim(-0.7, n_layer-0.3)\n",
    "plt.ylim(2.5, -0.5)\n",
    "\n",
    "scatter = plt.scatter([0, 0, 0, 0], [0, 0, 0, 0], c=\"white\", alpha=0, s=[25, 50, 75, 100])\n",
    "handles, labels = scatter.legend_elements(prop=\"sizes\", alpha=0.6)\n",
    "legend2 = ax.legend(handles, labels, title=\"Size\", prop={\"size\": 13}, loc=(-0.5, -1.75), frameon=False)\n",
    "plt.xlabel(\"Annotation\", fontsize=14)\n",
    "plt.setp(legend2.get_title(),fontsize=13)\n",
    "\n",
    "f.add_subplot(111, frameon=False)\n",
    "plt.tick_params(labelcolor='none', which='both', top=False, bottom=False, left=False, right=False)\n",
    "plt.scatter([], [], c=\"red\", s=75, label=\"Enrichment\")\n",
    "plt.scatter([], [], c=\"blue\", s=75, label=\"Depletion\")\n",
    "lgnd = plt.legend(scatterpoints=1, loc=(-0.5, -0.8), frameon=False, prop={\"size\": 13})\n",
    "plt.rcParams['legend.title_fontsize'] = 13\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2bc01e2-bcfb-49c2-be04-974e2b9e8f27",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}