{ "cells": [ { "cell_type": "markdown", "id": "7554fd94", "metadata": {}, "source": [ "# Pseudoreplication\n", "\n", "A well of a few hundred cells is not a few hundred independent observations of a treatment. The cells\n", "share the well's confluency, focus and plate position, so they carry the well's nuisance shifts in common.\n", "Counting each cell as a replicate inflates the sample size by the cells-per-well factor and manufactures\n", "significance from noise.\n", "\n", "[Which measurements moved](../phenotypes/which_features_moved.ipynb) covers the aggregate route:\n", "`tl.differential_features` refuses cell-level input and asks you to aggregate to wells first. This page is\n", "about the other route, testing at cell resolution with `tl.hit_calling` and `tl.edistance`, where the fix is\n", "not to aggregate but to make the null resample whole wells. `mt.metrics.diagnose_testing` shows, on your own\n", "screen, what the difference costs." ] }, { "cell_type": "markdown", "id": "cb5f3777", "metadata": {}, "source": [ "## A screen where nothing is real\n", "\n", "We build a plate where no perturbation has any effect, so every hit is a false positive by construction.\n", "The wells still differ, but only through nuisance: a row and column gradient across the plate, and a\n", "confounder that ties some features to the well's cell count. Any test that calls a hit here is wrong." ] }, { "cell_type": "code", "execution_count": 1, "id": "45c639d6", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T01:08:39.552906Z", "iopub.status.busy": "2026-10-21T01:08:39.552520Z", "iopub.status.idle": "2026-10-21T01:08:44.372170Z", "shell.execute_reply": "2026-10-21T01:08:44.370941Z" } }, "outputs": [ { "data": { "text/plain": [ "AnnData object with n_obs × n_vars = 3873 × 10\n", " obs: 'Metadata_Plate', 'Metadata_Well', 'Metadata_Row', 'Metadata_Col', 'Metadata_Perturbation', 'Metadata_Control', 'Metadata_Control_Type', 'Metadata_Batch', 'Metadata_ImageNumber', 'Metadata_CellCount', 'Metadata_ImageID', 'Metadata_ObjectType', 'Metadata_ObjectNumber'\n", " var: 'object', 'feature_group', 'feature', 'channel', 'scale', 'angle', 'gray_levels', 'radial_bin', 'params', 'is_feature', 'feature_kind'\n", " uns: 'mantispy'\n", " layers: None (.X)" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import warnings\n", "\n", "import mantispy as mt\n", "\n", "# Render warnings as their message alone; the kernel temp-file path in the default format is noise here.\n", "warnings.formatwarning = lambda message, category, *args, **kwargs: f\"{category.__name__}: {message}\\n\"\n", "\n", "cells = mt.ds.synthetic_plate(\n", " n_plates=2,\n", " n_wells=48,\n", " n_cells=40,\n", " n_features=10,\n", " n_perturbations=11,\n", " effect_size=0.0,\n", " row_gradient=3.0,\n", " col_gradient=3.0,\n", " confounder_effect=3.0,\n", " seed=0,\n", ")\n", "cells" ] }, { "cell_type": "markdown", "id": "c31a6aa8", "metadata": {}, "source": [ "It is annotated at cell resolution, and each well holds tens of cells that share that well's position on the\n", "plate." ] }, { "cell_type": "code", "execution_count": 2, "id": "b9fcd532", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T01:08:44.375853Z", "iopub.status.busy": "2026-10-21T01:08:44.375645Z", "iopub.status.idle": "2026-10-21T01:08:44.381697Z", "shell.execute_reply": "2026-10-21T01:08:44.380705Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "resolution: object\n", "wells: 96 | median cells per well: 40\n" ] } ], "source": [ "counts = cells.obs.groupby([\"Metadata_Plate\", \"Metadata_Well\"], observed=True).size()\n", "print(\"resolution:\", mt.get.resolution(cells))\n", "print(\"wells:\", counts.size, \"| median cells per well:\", int(counts.median()))" ] }, { "cell_type": "markdown", "id": "8f7d243f", "metadata": {}, "source": [ "## The nuisance is real even though the biology is not\n", "\n", "Read one feature across the plate. Its value drifts with row and column: neighbouring wells look alike and\n", "distant wells look different, for reasons that have nothing to do with any treatment. Every cell in a well\n", "inherits that well's value, which is exactly why the cells are not independent." ] }, { "cell_type": "code", "execution_count": 3, "id": "1edf30fb", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T01:08:44.384914Z", "iopub.status.busy": "2026-10-21T01:08:44.384777Z", "iopub.status.idle": "2026-10-21T01:08:45.286776Z", "shell.execute_reply": "2026-10-21T01:08:45.285643Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAArcAAAEOCAYAAACTsD9XAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjIsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvgI3uAAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAS6dJREFUeJzt3XlcVFX/B/DPMMCMsgoSgiwuuZtaLoH7imik5m6mue+pZVnkY6ZZ+JSVjz6aG4+aZtpClimmKRouqLjiviuCqIgwiDDAzP39EczPSZB7Z+44M/h5v1739XLunPnOuQpnvp753nMUgiAIICIiIiIqBxys3QEiIiIiIrkwuSUiIiKicoPJLRERERGVG0xuiYiIiKjcYHJLREREROUGk1siIiIiKjeY3BIRERFRucHkloiIiIjKDSa3RERERGQkIyMDCQkJyMjIKPF5QRBw5coVXLp0CYWFhU+9f0/C5JaIiIiIAABnz57FsGHDUL9+fYSGhmLXrl2PtVmyZAmCg4PRuXNndOnSBQEBAfjhhx+s0t+SMLklIiIiIgDAyZMn0a5dO5w9e7bUNqmpqThw4ACuXLmCq1ev4oMPPsAbb7yBCxcuPNW+lkYhCIJg7U4QERERke3Iy8tDhQoV8OOPP6Jv375PbFtQUACVSoU1a9ZgyJAhT62PpXG0dgeIiIiISLy8vDzk5+eLbi8IAhQKhdE5lUoFlUolS3+OHTsGQRBQvXp1WeKZi8ktERERkZ3Iy8tD9WBXpN3RiX6Nq6srHjx4YHRu1qxZ+Pjjj83uT05ODkaPHo0OHTqgdevWZseTA5NbIiIiIjuRn5+PtDs6XD0SDHe3sm+d0mTrUb3pdSQnJ8Pd3d1wXo5Z27y8PPTq1QsFBQXYuHGj2fHkwuSWiIiIyM5UcBVQwbXs26YKim6tcnd3N0puzaXVatGrVy8kJydj9+7d8PHxkS22uZjcEhEREdkZPfTQi2wnt+LE9tq1a4iLi0OVKlVkfw9zMLklIiIisjM6QYBOxIJXYto8KjMzE+fOnTPcsHbx4kUkJCSgSpUqqFatGgCgb9++2L9/P9asWYPr16/j+vXrAICgoCD4+/ubdD1y4lJgRERERHZCo9HAw8MD18/5i665Da6biqysLFFlCfv27cO0adMeO9+vXz/D+ZCQkBJfO3HiRJtYCozJLREREZGdKE5ur57zg5uI5DY7W4/qdW+JTm7LA5YlEBEREdkZPQToUfb8pJg25Q2TWyIiIiI7Y6ma2/KAyS0RERGRndEXHWLaPWuY3BIRERHZGR0E6ESUHIhpU94wuSUiIiKyMzrh70NMu2cNk1siIiIiO8OyhNIxuSUiIiKyM3oooINCVLtnDZNbIiIiIjujF/4+xLR71jC5JSIiIrIzOpEzt2LalDdMbomIiIjsDJPb0jG5JUkEQYBer4dSqbR2V4iIzKLX/32rjYND2VuYEtkavaCAXhBRcyuiTXnD32iS5Oeff0aDBg2s3Q0iIrNNmDABkydPtnY3iExSPHMr5njWcOaWjAiCAJ1OBxTNZpg7o2Gpmd7CwkLDn5VKJRSKZ++Xl4ieTK/XG2Zn5RgnONNLtkQHB+hEzFHqnkpvbAt/Q8nIzp074eTkBLVaDWdnZwQEBGDatGnQarUmxdu4cSMaN24s+XV//fUXGjVqBEdHR1SvXh3fffed4bnr169DrVZDrVbDyckJ//nPf0zqGxGVbx9++KFhPFOr1ahbty5Wr15tcrwxY8bgnXfekfy6RYsWoWrVqnByckJISAiOHz9ueC4pKQn9+/eHl5cXKlWqhF69euHKlSsm95GeHUJRWUJZh8CyBKK/5eXloaCgAD/99BNWrVqFf//73yW2EwQBhYWFRjOpjyqe6XhSm39KS0tDREQEBg8ejPv37+Pf//43RowYgYMHDwIAgoODDfHq1Klj8jUSUfnXrl07FBYWIjs7G9OmTcOIESOwb9++Etvq9XoUFhYavr16lCAIhqO0NiX57bffEBkZidWrVyM9PR2tW7dG9+7dkZOTAwD4+uuvMWDAAFy+fBlnz56FXq9Hz549zbxqehawLKF0TG6pVAqFAiEhIRg8eDC2b99eYpu//vrLMCvi6uqKjh074vz58wCAI0eOYMiQITh9+rShTfGHyu+//46XXnoJFSpUQM2aNbFixQpDzO+//x7e3t54//334ebmhv79+6NLly5GbYiIpHB2dsbo0aNRp04d7Nixo8Q2s2fPhlqthkqlgo+PD0aOHIns7GwAwBdffIFVq1Zh0aJFUKvV8PDwAIr+4x4VFYWaNWuiQoUKaNGiBQ4cOGCIuWzZMgwaNAhdunSBh4cHPv30Uzx8+BCbN28GAPzvf/9Dnz59UKlSJVSpUgXTpk3DqVOnoNFonsrfC9kvneAg+njWPHtXTJLp9fpSa9WKZ0UKCwtx69YthISEYODAgRAEAU2bNsV3332HBg0aGNq0atUKe/fuxeuvv45PPvkEmZmZ2LBhA2bOnIlt27YBAE6cOIGmTZsavU/Tpk1x4sSJp3K9RFR+PWk8mz17tmGsOnr0KJKTkzF79mwAwPTp0zFixAhMnjwZhYWFePDgAQBg7ty5+P7777Fp0yZkZmZi0qRJiIiIQHp6OlDCeKZSqdCgQYNSx7O4uDjUrFkT7u7uFrh6Kk/0UEAPBxEHZ26JAAA6nQ75+fmIj4/H+vXrER4eXuZrXFxc8PHHH+PcuXO4efNmqe0WLVqEkSNHomvXrlAqlXjxxRcxZswYrF+/HgCQnZ1tmBUp5unpaZhBISISq7iM4MGDB1i8eDEuXryIsLCwMl/n7++PyMhIwwxraf773//is88+Q7169aBUKvH666+jXr162Lp1KyBxPNu+fTu++OILLFu2TPJ10rOHZQml42oJVCIXFxcolUpUrVoVY8eOxXvvvVdiu7S0NLz11luIi4tDRkYGHBwcoNPpkJKSgsDAwBJfc+3aNfz0009YtGiR0fm2bdsCANzd3ZGVlWX0XGZmJmcyiEiy4tIplUqFmjVr4ttvv0VoaGiJbfft24f3338fJ06cwMOHD6FQKKBWq0uNnZ2djXv37qFXr16PPdetWzdAwngWGxuLQYMGYf369ejUqZOJV0vPErElBzrh2dt/lzO3VKK8vDxotVpcuXIF8+bNg7Ozc4ntpk+fDr1ej6SkJBQWFiI3Nxcqlcpw85iDgwOEf/xiVatWDe+//77h67/iY9euXQCAJk2a4PDhw0avO3z4MJo0aWLRayai8qe4dConJwcnT57EG2+8UWrb/v37o2fPnrh58yZ0Oh127NhhdCPsP8czNzc3eHl5Yfv27Y+NZzNmzACKxrNDhw4ZXpOXl4ekpCSj8eznn3/GoEGDsGHDBt5MRqL9XZYg7njWMLkls+Tm5kKpVMLZ2Rm3b9/GpEmTjJYN8/PzQ0pKClJSUgwfEpMnT8bixYvx448/4v79+7h27RoWLFiAxYsXAwAGDRqErKwszJkzB+np6Vi7di127tyJMWPGGOLqdDpDvOI7nItXZiAikkoQBOTm5sLJyQmOjo44ffo0IiMjjdr4+fnhzJkzePDggWG1hKlTp2LKlCk4cOAAsrOzkZSUhAkTJuDw4cMAgPHjx2PDhg3YvHkz7t69i/feew8eHh6IiIgAAKxduxYjR47Ezz//jM6dO0taWYaebfqidW7LOvQmpnoZGRlISEhARkZGqW3u3LmDY8eOPfbthLUxuSUjDg4OT1zs3MHBAY6O/1/N8umnn+LWrVsIDg5GixYtUK1aNfj6+hpe36pVK/Ts2RONGjVChQoVsG/fPrRq1QoxMTH473//i+effx4dO3bEzZs3DTMqzz33HGJjY/H7778jKCgIn376Kb777js0a9bM8L6+vr5Qq9W4dOkSpk+fDrVajZkzZ1r874eI7IdSqXziBjKPPq9QKLB69WosXboUPj4+GDhwIHr16mU03o0ePRp5eXnw8/Mz1NF++OGHGDt2LMaNG4eAgAAMGTIEjRs3NtxE1r17dyxYsABTpkxBtWrVcPLkScTGxqJixYpA0Rj64MEDdO3a1bCqjFqtNqw6Q1SaAsFR9CHF2bNnMWzYMNSvXx+hoaGGb1UfpdPpMHr0aAQFBWHgwIHw9fXFZ599JuPVmUch/PM7YyIiIiKySRqNBh4eHlh77AVUdCt798+H2ToMeTEJWVlZou5d2bhxIx4+fIhevXrBy8sLP/74I/r27WvUZsGCBZgzZw4SEhJQu3Zt7Nq1C126dMGWLVtE3YBuabyhjIiIiMjOiN9+V9oc5oABA4Ci+vDSrFy5EoMGDULt2rUBAB07dkSbNm0QHR1tE8ktyxKIiIiI7IxecBB9yCkvLw9nz55F8+bNjc6//PLLOHbsmKzvZSrO3BIRERHZGakzt//c9U6lUkGlUkl+38zMTOj1enh5eRmd9/b2xr179yTHswTO3BIRERHZGT0AnaAo8yheRygwMBAeHh6GIyoqyqT3dXJyAgCjlZFQtHpS8XPWxplbIiIiIjujF7nMV3Gb5ORkoxvKTJm1BQAvLy+4uroiNTXV6HxqaiqCg4NNiik3ztwSERER2ZniHcrEHCjaLe/Rw9TkVqFQoGPHjkZbUxcWFmLr1q02s7veU5251ev1SE1NhZubW6nrqBJR+SQIArKzs+Hv7w8HB/v+fzXHMqJnmy2MZ2J3H5O6Q1lmZibOnTuH/Px8AMDFixeRkJCAKlWqoFq1agCAjz/+GC1btsTUqVPRtWtX/O9//0N+fj6mTp1q4tXI66muc3vz5k0EBgY+rbcjIhuUnJyMgIAAa3fDLBzLiAhWGs+K17n9OrElKriWPUeZ+6AQbzfbL3qd23379mHatGmPne/Xr5/R+cTERHz99de4efMm6tSpgw8++AA1atQw4Yrk91Rnbt3c3AAATbt+CEcntWxxc73LXsRYquf+TJY9JgDc6SL/B6JDgewhAQAVb8u/BWSFZPm36Mv8zDLb7nqOui97zPRXasse0/WWZX4AlFPvyBqv8GE+9g9YaRgH7FnxNTy/ciqUFU37aq80D2+6yhoPAJ7fkCN7TADI9a0oe0zXk6kiWkmTW9tX9pgAcKepvP/2AHBozErZY3abMkz2mACQVV3+FELrKXtIAECVQ/myxisszMOh3fOsOp6JXy1B2sxyq1atkJCQUGa7Zs2a4bvvvpMU+2l5qslt8dd3jk5qWZNbpbP8ya2jg/yDFgAoneW77mIOFvpW1NFJ/uTWUVn6otCmUrpYJrl1VDjLHtMS//6OjvL//AOA0sUyvwPl4Wv84mtQVlTJntw6VLDAz4hSJ3tMFI3lsse0wNjr6Ch/PwFAaWLN4pO4u8n/Fbcl/p0AQKmSP4VQWqarcHS0TOmANcczvaCAXhBRliCiTXnD1RKIiIiI7Ixe5MytmBUVyhsmt0RERER2RuzuY3LvUGYPmNwSERER2RkdFNCJWAlBTJvyhsktERERkZ3hzG3pTEpuMzMzsWfPHkRERECptMzNLERET8OlS5eQkpKCdu3aWbsrRESi6UTOylrmdlJ5PXjwAEuXLsXJkyeRkZFh9NzAgQPxxhtvSIpnUnK7bNkyLFiwAMuWLUOPHj1MCUFEZBMmT56M06dP48iRI6hcubK1u0NEJEp5mbnNz89HaGgo1Go1lEol9Ho9qlevjq1bt8Lb2xteXl6SY0q+YkEQsHHjRsTExGDFihWS35CIyFZcv34dubm5eP/997F69Wprd4eISDSp2+/aqi1btsDZ2RkHDx5EREQEwsLCsHHjRpw5cwb5+fkmbZIh+Yrj4uLQsmVLhIaGQqvV4ubNm6W21Wq10Gg0RgcRka2Ijo7GxIkTMXToUPzwww+ltuNYRkS2RijafresQ7DxG8ouXbqEDh06wMHBASqVCjk5f286ExgYiAEDBmDHjh2SY0pObpcvX47AwEBs2rQJderUQXR0dKlto6Ki4OHhYTi4XSUR2QqdTofvv/8eAPDnn3/CxcUFe/bsKbEtxzIisjXlZeZWq9VCrf57946qVavizJkzhuc0Gg10OulVw5KuOD09Hfv27cOBAwewevVqJCcnY/369dDrS94hKjIyEllZWYYjOdkyW9oSEUm1ZcsWODk5Yd26dVi9ejWcnJywfPnyEttyLCMiW1O8Q5mYw15069YN+/fvx5AhQzBu3DisW7cOnTp1khxH0g1lq1evxrBhw/DJJ58YzvXr1w+xsbF45ZVXHmuvUqmgssD2hERE5lq+fDkWLVpkGDj1ej1q1KiBjIyMx25g4FhGRLZGJ3KHMjFtrGnkyJEQBAEAUKlSJezatQtfffUVUlNTsXHjRjRt2lRyTEnJ7Z07dzB8+HCjc+PGjcPp06clvzERkbXk5ubC3d0d7du3N5xzcHDA22+/jePHj6Njx45W7R8RUVnEzsra+sytn5+f0ePmzZsbSsZMJSm5/fzzzx8716lTJ5OmjImIrKVChQpYv379Y+enTJlilf4QEUmlhwP0ImZlxbSxBQ8fPsTZs2fh6emJmjVrmhXLPq6YiIiIiAwK9A6iD1u3fv16+Pv7o1mzZli1ahUA4NChQyZ/i8btd4mIiIjsjCByEwfBxldLuHz5MiZMmIA1a9bgyJEjhkUKWrRoAQDYs2eP5B0kbfuKiYiIiOgxOihEH7YsLi4Or732Gnr27PnYjbtNmzbFvn37JMdkcktERERkZ/SC2OXArN3TJ9NoNHBzcwMAKBTGiXhmZiacnJwkx2RyS0RERGRn9EVlCWIOW/byyy/j999/R05OjlFye/bsWWzYsAGtW7eWHJM1t0RERER2pnh7XTHtbFmrVq3QokULNGnSBJUrV4ajoyMGDRqEX375Bb169UJoaKjkmOUiufXdclX2mOe+8BPRygR35P9+oOIty/yvzOvEQ9ljXuvrI3tMp+2yh/xbuPx99dl0TvaYD0Oelz0mAGT+FCBrPF1+nqzxbIHPcjUcHdWyxnTSyP97pyiQvn2lGLdaKWWPGfTAV/aYD6tI/1pTDKVW/pg1doyQPaZHoGU+6tX35P888zlhgb9UAHdelPf3VKfVA3/KGlJ6HwQFdCLWsBXTpiQpKSm4e/cuqlatCh8f+T8PH7V+/XpER0fjp59+QmpqKpycnDB//nyMHz/epHjlIrklIiIiepaILTmQWpZw7do19OvXD5cvX0ZQUBAuXbqEtm3bYsOGDXB3dzejxyVLSEiATqfD6NGjMXr0aFli2nYhBhERERE9Rg8xN5OJK1141NSpU6FQKJCSkoLjx4/j8uXLOHbsGKKioixyHfHx8di8ebOsMZncEhEREdkZoajmtqxDkJjc3r59G02bNkWFChUAAL6+vqhbty5u375tketo1KgRDh48KGtMliUQERER2ZnimVkx7aT417/+hZEjR+LFF19EnTp1sG/fPpw/fx4LFy40o7elq1+/PjIyMjBo0CAMGDAAlStXNno+KCgIQUFBkmIyuSUiIiKyM1JrbjUajdF5lUr12KYJANCmTRt069YN06dPR1BQEK5du4YpU6agXr16Mvb+/61ZswYnT57EyZMnsWHDhseenzFjBubOnSspJpNbIiIiIjsjdeY2MDDQ6PysWbPw8ccfP9a+X79+0Gq1uHnzJlxdXXH79m20bNkSOTk5+Oqrr2S8gr9FRkbi3XffLfV5R0fpqaqkmludTgdPT094enqicuXKaN26NY4cOSL5TYmIrG3Pnj2G8czPzw+9e/dGRkaGtbtFRCSKmHrbR9fCTU5ORlZWluGIjIx8LGZeXh527NiBkSNHwtXVFSiquR0wYAA2bdpkketQKpVQq9WlHgUFBZJjSkpuBUFAXl4erl27hnPnzqF37954++23Jb8pEZG1FRQUoFmzZrh27RqOHDkCZ2dnLFiwwNrdIiISRdzWu/8/u+vu7m50lFSSoFKp4OrqipSUFKPzqampj9XCWtKDBw8QHR2Nli1b4rPPPpP8epPKEjw9PQEAVatWNdxNVxKtVgut9v8XZP5nvQcRkTU5OjrC09MT7u7u8Pb2LnU841hGRLbGEjeUKRQKTJw4EVFRUXB1dUWDBg2wd+9erFu3DtHR0Wb2uGwJCQmIjo7Ghg0b4OHhgR49eqB///6S40hObrVaLTw9PVFYWAiFQvHEtcmioqIwe/ZsyZ0iInoadu7cCU9PT+Tl5cHf3x8HDhwosR3HMiKyNZZaLeHTTz9Fw4YN8dtvv+Hnn39GQEAAtmzZgq5du5rR29Klp6dj7dq1iI6OxpUrV+Dv749hw4Zh4cKFUChM211N8jq3KpUK165dw40bN/Drr79i8ODByM3NLbFtZGSkUX1HcnKySZ0kIrKEdu3a4dq1a0hOTsa4ceMwcuTIEttxLCMiWyO1LEEsBwcHDB48GBs3bkRcXBzWrl1rkcT25MmTGDBgAKpWrYro6GiMGTMGqampGDZsGDw8PExObGFuWUJwcDAyMjKQnZ1d4td5pS0zQURkC4rLElBUZnXlypUS23EsIyJbIxTdVCamnS367bffEBMTg6+++gpvvfWWrLFNLkvQ6/VwcHDAzJkz8dxzz8naKSKip6G4LKGgoACVK1fGN998Y+0uERGJYqmyhKclIiICR48exbRp07B69WqMHj0agwcPliW2pOTW0dER9+/fB4qWbnBzc5OlE0RET1v79u1x9+5doGhm9kk3xxIR2Rp7T26bNGmCmJgY3L59G2vWrMGCBQvw7rvvomrVqggLC4MgCE+v5rZ4XUgmtkRkz4pLEjw9PZnYEpHdsVTN7dPm6+uL6dOn49y5c4iNjUVoaCj+97//ISgoCBMnTkRSUpLkmJKTWyIiIiKyrvKS3D6qTZs2WL16NW7duoUZM2bg0KFD+OGHHyTH4fa7RERERHZGEBQQRCSuYtrYGnd3d4wbNw7jxo1DTk4OAGDp0qXQarWYMmVKma9ncktERERkZx7dWresdvbMxcUFKFoPNy8vT9RrmNwSERER2Rl7v6HMkpjcEhEREdkZnd4BCn3Zt07pRLQpb6yS3GoCHaFUyffWrhc9ZItVrOZivewxAcDx0mXZYypUzrLHBAC9t7vsMYM3Z8oes9LiW7LHBIB8vVL2mGeeryd7zHx3yyzR7XZV3v/tC2bsNmOrClyUEJzk/TnRqdWyxgOAh5UtM9T7HpJ/nFTGHZU9Zv74UNljAoDLLfmv3+uM/OPOrdaWGSP0d+T/nfZOzJY9JgAErLoqa7xCIR/nZI0oXXmuuTUXZ26JiIiI7IwgsiyByS0RERER2TwBgCBiUt5Wt9+1JCa3RERERHZGDwUUz8BqCcXCwsKg0+lEtWVyS0RERGRnylvN7ZYtW5CQkICXX34ZERERSE5ORnZ2NurXrw8AaNGihehYTG6JiIiI7IxeUEBRTpYCGzJkCGJjY+Hh4QFBEBAREQFHR0e89tprOHHiBNQSb7R99taHICIiIrJzgiD+sGXx8fGIj4/HhQsXMGLECMN5Pz8/NGnSBDExMZJjmjRze+jQIezevRsFBQXo3r07XnzxRVPCEBFZlSAI2Lx5M44dOwZPT0/069cP/v7+1u4WEVGZyktZwtGjR9GrVy94eXlB8Y8lI2vUqIFLly5Jjil55vbjjz9Gnz59cO/ePWi1WowbNw6//vqr5DcmIrImnU6HsLAwzJ07F4Ig4Nq1awgPD8eNGzes3TUiojIVJ7diDlumUqlw//79Ep87deoUfHx8JMeUNHN7/vx5LFmyBKdPnza82axZs3Dz5k3Jb0xEZE3/+9//UFBQgIMHDxpmC0obYImIbE15qbkNDw/HBx98YDQW63Q6LFy4ENu3b8fixYslx5SU3B46dAht27Y1yqKVSiWCg4NLbK/VaqHVag2PNRqN5A4SEVlCQkIC+vbta/Q1WKVKlUpsy7GMiGyN2HpaW6+5rVatGhYsWIB27drB2dkZarUa//nPf5CXl4fly5cjKChIckxJya2joyMKCgpEt4+KisLs2bMld4qIyNKkjGccy4jI1vyd3IqpuTX9Pa5evQqlUmlSginFsGHD0LlzZ2zevBmpqanw8fHBq6++iurVq5sUT1LNbcuWLbF3715cvHjRcC4nJwdJSUklto+MjERWVpbhSE5ONqmTRERy69ixI9asWYPc3FzDuatXr+L27duPteVYRkS2xpI1t3v27EHt2rURGhqKsLAwdOjQAampqRa5jmIBAQEYP348PvnkE0yePNnkxBZSZ26Dg4Mxb948hIaGolOnTlCpVEhMTCy1HkKlUkGlUpncOSIiSxkwYAC2b9+OOnXqoGPHjtBoNLh69SpiY2Mfa8uxjIhsjSBya12pE7dnzpxBeHg4PvzwQ/zrX/+CQqHAvn37cP36dYutJlNQUIDo6Ghs27YNKSkp8Pf3R9euXTFq1Cg4OztLjid5KbDRo0cjPDwc+/fvR0FBAebNm8elc4jILkVHR+PEiRM4duwYKlWqhM6dO8PFxcXa3SIiKpOllgL75JNPUKdOHcycOdNwrlWrVib1UYz8/Hy0bdsWFy5cwGuvvYYXXngBaWlp+Oijj7Bq1SrEx8dL3sTBpHVuAwMDMWDAAFNeSkRkUxo3bozGjRtbuxtERNJYaOp2x44dmDBhAvLz83Hx4kX4+vqicuXKpvayTD///DPu3LmDCxcuGL3P559/jpdffhkbNmzAsGHDJMXk9rtERERE9kZsPW1Rm3+u8lJSuVV+fj7u3buHW7du4fnnn4erqyuSk5PRokULrFu3Dn5+fvJeA4Dbt28jIiLisQS6UqVK6NmzZ4n3QZSF2+8SERER2Rmp2+8GBgbCw8PDcERFRT0WU6lUAgBiYmIQFxeHM2fO4MaNG8jIyMCECRMsch0vvPACEhMTodfr/3F9AhITE9GwYUPJMTlzS0RERGRnpNbcJicnw93d3XC+pJtklUolqlatik6dOqFmzZpA0QzqkCFDMGfOHFn7X6x27dp48OABwsLCMHz4cPj5+SEtLQ3ffvstbt26BXd3d+zduxcAEBQUJGpZMia3RERERPZGUBhKDspsB8Dd3d0ouS1NWFjYYzvPpqWlwdvb2/S+PsGaNWsMS8ru3Lnzsefbtm1r+POMGTMwd+7cMmMyuSUiIiKyM5baoWzGjBlo1qwZZsyYga5du+LkyZNYvHgxvvjiC5P7+iSRkZF49913RbV1dBSXtloluXVJ08HRSSdbvHNjPWWLVcztilL2mABQNcND9ph3v7LMvtHe78u/Z5/i7j3ZY578rZ7sMQHAPz5H9piVqsr3c1/szmtaEa2kK0ivKGs8nWW6aVWuh2/A0UH6GoxPcuHdmrLGAwCXFMuMERmN5I/r4RMqe0ytp2Wuv8BF/rg6J/ljul2VPeTfFPJ/Rtx/0TKzg455XrLGKyzIA36TNaR0FlotoWbNmti/fz8+//xzfPjhh/D398f333+PHj16mNrTJ1IqlYZaX7lw5paIiIjIzlhqnVsAqFevHlatWmViz8xTWFho9NjBwQEODtLWP+BqCURERET2SBBx2IH4+Hi0bt0a7u7ucHJyMjo++ugjyfE4c0tERERkZyw5c/s0paamokePHpg6dSpq1aqFwsJCdO/eHStXrsTVq1cxfPhwyTE5c0tERERkb4pXSxBz2LAdO3agQ4cOmDVrFmrWrInAwEAMGjQI27dvh4+PD86dOyc5JpNbIiIiInsjpiTBDkoTUlJSULduXQCAm5sbMjMzgaIbzTp27Ihjx45JjsnkloiIiMjelJPkVq/XG24Yq127NuLi4pCXlwedToeEhARUqlRJckxJNbfHjx8HAKjVagQEBMDV1VXyGxIRWduDBw9w6dIlKBQKuLi4IDg4GE5OTtbuFhGReBI3cbBV7u7uht3SwsLCUKFCBQQFBcHZ2RmFhYXo37+/5Jiik9vCwkK89NJLaNSoER4+fIibN2+iU6dOWLVqFSpXriz5jYmIrCUhIQG9e/dGjRo1kJGRgczMTIwaNQpffvklFArb/iAgIoIFN3F42iZPnmz4s1KpxN69e7F582bk5OTg1VdfhY+Pj+SYksoSnJ2dcfz4cVy4cAFpaWlwd3fHuHHjJL8pEZG1tWzZEsePH8eNGzdw/Phx7Ny5E4sXL7Z2t4iIxCknZQn/VLFiRQwYMAAjRowwKbGFOUuBubu7Y+HChQgICDCql3iUVquFVvv/2xJpNBpT346IyGJq1KiBjz76CMuWLcOkSZMee55jGRHZnHJSloCiutv4+Hhcu3YNBQUFRs+9+OKLaNq0qaR4Zq1z6+XlBb1eD41GA0/Px7fAjYqKwuzZs815CyKip6JKlSq4c+dOic9xLCMiW6MQxO2AbIFdkmWVl5eHdu3a4ciRI6hSpQocHY1T07feeuvpJrdHjhyBp6dniYktAERGRuKdd94xPNZoNAgMDDTnLYmILOLAgQOoU6dOic9xLCMimyO25MDGk9uYmBhkZWXh5s2bqFKliiwxJSW3giDg+PHjePjwIU6fPo3PPvsMM2bMKLW9SqUy3AFHRGRLsrOzcfz4cWRmZmLv3r2YP38+du7cWWJbjmVEZHPKSVlCeno6unfvLltiCynJrUKhQL169TBs2DA4OzujWrVq+OKLL9C3b1/ZOkNE9DS4ubkhJycHw4cPR8WKFdGgQQP8+eefkr/6IiKymnIyc9u4cWNs3rxZ1piik1ulUmlY55aIyJ69/PLLHM+IyL7ZcXKbnJyM69evA0X55d27dzFmzBj07t37sT0UgoKCEBQUJCm+WTW3RERERGQFdpzcrlmzBjNnzjQ6d+LECaxYseKxtjNmzMDcuXMlxWdyS0RERGRv7LjmNjIyEu+++66otv9cPUHUa0zoExERERFZkT0vBaZUKqFUKi0WX9IOZURERERkA8rBDmUpKSlIS0szPNbpdJg3bx7atWuH119/HRcvXjQpLpNbIiIiInqqBEHAK6+8gnv37hnOzZ8/HzNmzICbmxuSkpLQvn175OfnS47N5JaIiIjIzigeKU144mHtjpbir7/+QuXKldGgQQPDuWXLluFf//oXfv/9dxw9ehReXl7YsmWL5NhWqbnduXoy3N3drfHW1veVtTsgQTdrd4CsaqC84TQaDTy+/lDeoFYWc/HLZ3csI3qGaTQaeHjMFNHSgp7CDWV5eXk4deoUvL29Ub16dZPjlOTEiRNo1KiR4fGNGzdw9epVDBz494ePk5MTunbtiitXrkiOzZlbIiIiInvzFGpuJ02ahBYtWjy2bJcc1Gq1UUnCgQMH4OHhgbp16xq1UyikJ+dMbomIiIjsjYWT240bN+LYsWNo06aN3D0HijbT+eWXXxAfH4+7d+9i0aJF6Nq1q1Eye+rUKaPZXbG4FBgRERGRnbHkUmBXrlzB22+/jbi4OEycONGk/pWlcePGeOONN9C2bVsAgLu7O5YuXWp4/vr167hy5Qo6dOggOTaTWyIiIiJ7Y6EdygoKCjBo0CDMmjULderUMbV3oixZsgQjRozA9evX0bJlS/j5+Rme0+l02Lhxo0nr4TK5JSIiIrI3EpNbjUZjdFqlUkGlUj3W/MMPP4Svry/Gjh0rV0+fqFmzZmjWrNlj52vUqGFyTNbcEhEREdkZUcuAPVK6EBgYCA8PD8MRFRX1WMyEhAQsWbIE48ePR2JiIhITE5GdnY2MjAwkJiaisLDQ4tdVWFhodOj1eskxJM/c6nQ6NG7c2OjcgAEDLHInHRGRJSUkJGDUqFFG5+bPn4/w8HCr9YmISBSJS4ElJycbLV1Y0qytRqNBvXr1jHK6CxcuwMnJCePGjcMff/wBb29vua7AID4+HpGRkTh58iSys7ONnpsxYwbmzp0rKZ7k5FYQBFy8eBFHjhwxnPPy8pIahojI6h48eABPT0+jmxgCAwOt2iciIjEU+r8PMe1QdMNWWetyh4WFISwszOhc586dUaVKFaxbt86s/pYmNTUVPXr0wNSpU1GrVi0UFhaie/fuWLlyJa5evYrhw4dLjmlSza1CoUDDhg3LbKfVaqHVag2P/1nvQURkba6urmWOZxzLiMjmWOiGsqdtx44d6NChA2bNmoW5c+fi4cOHGDRoEPr374+WLVvi3LlzqFmzpqSYJtXc5ufno2HDhoajtGw+KirKqL6DMyJEZGv27dtnNJ6dP3/+sTYcy4jI5oittzUzua1Tp45ZN3eVJSUlxbBxg5ubGzIzMwEASqUSHTt2xLFjxyTHNGnm1snJCRs2bDA8fnTphkdFRkbinXfeMTzWaDT8UCAim9K4cWOjsoTg4ODH2nAsIyKb85RmbhcvXmxegDLo9Xo4OPw911q7dm0sXboUeXl5cHJyQkJCAvr27Ss5pkXLEkpbZoKIyFaIKUvgWEZENqeclCW4u7sbxtewsDBUqFABQUFBcHZ2RmFhIfr37y85pknJbXFZQrGwsDB89dVXpoQiIrKq4rKEYjNmzMCgQYOs2iciorJYcoeyp2ny5MmGPyuVSuzduxebN29GTk4OXn31Vfj4+EiOKTm5dXR0xMmTJ43OeXh4SH5jIiJrCw0NxYEDB4zO+fv7W60/RETPuooVK2LAgAFmxTBp5lZMSQIRka1zcXHheEZE9qmclCUU27JlCxISEvDyyy8jIiICycnJyM7ORv369SXH4va7RERERHamvJQlAMCQIUMQGxsLDw8PCIKAiIgIODo64rXXXsOJEyegVqslxeP2u0RERET2SBBx2Lj4+HjEx8fjwoULGDFihOG8n58fmjRpgpiYGMkxOXNLREREZG/KSVnC0aNH0atXL3h5eUGhMN5OuEaNGrh06ZLkmJy5JSIiIrIzYjZwEFu6YE0qlQr3798v8blTp06ZtFoCk1siIiIieyOmJMEOShPCw8Px66+/4uDBg4aZW51Oh6+//hrbt2/HK6+8IjkmyxKIiIiI7Ex5uaGsWrVqWLBgAdq1awdnZ2eo1Wr85z//QV5eHpYvX46goCDJMZncEhEREdmbclJzCwDDhg1D586dsXnzZqSmpsLHxwevvvoqqlevblI8JrdERERE9qacJLfffvsttFotRo8ejfHjx8sSk8ktERERkZ0pL2UJmZmZuHLliqwxeUMZERERkb0pJzeUdenSBX/88Qc0Go1sMTlzS0RERGRvyklZQmpqKgCgdu3aCA8PR+XKlY2e79y5M8LDwyXFZHJLREREZGfKS1nCnTt3EBwcjODgYKSlpSEtLc3o+UaNGkmOKTm5FQQBX375pdG54cOHw9vbW/KbExFZ040bN/DDDz8YHlesWBETJkywap+IiEQpJzO3gwYNwqBBg2SNKbnmVqfT4YMPPjBk12lpaSgoKJC1U0RET8OFCxewePFiw1h2+/Zta3eJiEiU8rJD2ZOkpqbi6NGjkl9nUlmCo6Mj5s+fb8pLiYhsSp06dTieEZH9KSczt/9UWFiILVu2YOXKlYiNjUVkZCReeuklSTFMSm51Op3hw6BSpUoYOXJkie20Wi20Wq3hsZx3whERyeHatWuG8axRo0YICwt7rA3HMiKyOeUsub148SKio6OxZs0a6PV69OnTB1u3bkWHDh0kxzJpKTBBEAxf4929e7fUdlFRUfDw8DAcgYGBprwdEZHFaLVaw3iWlZVVYhuOZURkaxQSDlPcuHED165dg06nk7nn/y83Nxdr165Fu3btUK9ePRw/fhzNmzfH6NGjsWTJEoSFhcHJyUlyXIuWJURGRuKdd94xPNZoNPxQICKbIqYsgWMZEdkahf7vQ0w7KZYtW4Z58+ZBp9NBEATodDosWrQIffr0MbmvJfnxxx8xevRoVKxYEcOHD8eaNWtQrVo1zJ07F3l5eWbFtuhSYCqVCiqVypJvQURkcRzLiMjmWKgs4erVq9izZw+CgoIAAPPnz8frr7+OU6dOoVatWqb1tQQXLlxATk4OJkyYgDFjxiA4OFi22JLLEhwcHPD222/L1gEiImsJDg5G3759rd0NIiLTWGB3snnz5hkSWwCYMmUKCgoKkJCQIGvX3377bURHR2Pv3r2oUaMGwsPDERMTI8sKXJJnbh0cHBAVFWX2GxMRWVutWrVknYkgInpantYmDidOnIAgCKhWrZp5gf6hYsWKGDp0KIYOHYoLFy4gOjoaEyZMwJ07d9C5c2ccPnwYzZo1g0IhvWrYpBvKiIiIiMiKxMzaPjJ7q9FojI5HV4ApzcOHDzF69Gi0bdsWrVu3ttil1K5dG//+979x8+ZNxMTEwNnZGaGhoQgKCjLaaEcsJrdEREREdkbqJg6BgYFGq76U9S28VqtF7969kZubix9++MGkGVSpHB0d0atXL/z++++4fv06xo0bh8zMTADAvXv3kJ6eLi6OhftJRERERHKTeENZcnIy3N3dDaefdJNsfn4+evfujStXrmD37t3w9fWVo8eSVK1aFTNmzDA8/uabb5CXl4e5c+eW+Vomt0RERER2RmrNrbu7u1FyW5rixPbixYuIi4uDv7+/DL19upjcEhEREdkbCy0F1q9fP8THx2Pt2rW4desWbt26BQAICAhAlSpVTOvrU8bkloiIiMjeWCi5TUlJQa1atTBnzhyj82+99RbefPNNacGshMktERERkZ2x1FJgiYmJJvfJVjC5JSIiIrI3Fpq5LQ+Y3BIRERHZGYUgQCGUnbmKaVPeMLklIiIisjflaOb2+vXrUKlUqFKlCnQ6HZYuXYpLly5h2LBhaNy4MVBU8yuITNS5iQMRERGRnZG6iYOtysrKQkREhCFxXbJkCWbOnIkTJ06gbdu2uHfvHgDAw8MDnp6eomIyuSUiIiKyNxK337VVsbGxaNKkCfz8/AAAa9euxYoVK7Br1y506dIFmzZtkhzTpORWEARRexITEdm6/Px86HQ6a3eDiEiS8jJzm5qaCh8fH6BoFjcpKQnh4eEAgDp16uD27duSY0pKbs+fP4/w8HC4urrC29sbLVq0wP79+yW/KRGRtW3btg2NGzeGh4cHvLy80Lt3b9y8edPa3SIiEqeczNzWrVsXmzdvxtWrV7Fo0SK0aNECLi4uQFHeWbduXckxRSe3Dx8+RPfu3dGxY0ekp6dDo9Fg6dKl2L59u+Q3JSKypqSkJAwePBjz5s3DgwcPcOfOHbz55pvYtWuXtbtGRCRKeZm5DQ8Px/PPP48aNWpg3rx5mD17NlB0k9mpU6fwyiuvSI4perWE3bt3w9fXF9OnTwcAZGdnIygoCJMnTy71NVqt1qh8QaPRSO4gEZHc1q5di7Fjx6Jbt25A0VdhrVq1gkqlKrE9xzIisjnlZLUEBwcHxMbG4ubNm/Dw8ICbmxsAQK1WY9u2baWOy0+MKbZhSkoKqlevbng8bdo01K1bF7Vr1y71NVFRUfDw8DAcgYGBkjtIRCS3f45nbdq0Qc2aNTF8+PAS23MsIyJbZO+zto8KCAgwJLYA4Ovri2rVqpkUS/TMbY0aNbBw4UIIggCFQoHly5fjq6++euIbR0ZG4p133jE81mg0/FAgIqurUaMGTpw4YXh8+vRp/P7771i9enWJ7TmWEZHNEYS/DzHtbExsbCx+/fVXUW27d++OHj16SIovOrlt37491Go1Ro8ejbfeegs+Pj5ISkqCg0Ppk78qlcqk6WQiIksaO3YsmjRpgpo1a+KVV16BWq3GsWPHSh3POJYRka0ROzNri7O3Dx8+RHp6uui2UolObpVKJf744w/MmTMHQ4cOxYMHD/Diiy9i586dkt+UiMiaAgICEB8fj1mzZuGbb76Bs7Mz2rdvjyVLlli7a0RE4thxzW2fPn3Qp08fi8WXtP2ul5cXFixYYLHOEBE9LfXq1cMPP/xg7W4QEZlEof/7ENPuWSMpuSUiIiIiG2DHM7c2U3NLRERERLZBoReg0JeduYpp87TZTM0tEREREdkGe76hzKZqbomIiIjIBthxWYKlid7EgYiIiIhsQ3nZfhcA0tLSMHr0aLzwwgv44osvAADHjx/H8uXLTYrH5JaIiIjI3hRv4iDmsGG5ublo27Yt7t27h2rVqiErKwsA0KBBA3z99ddITk6WHJPJLREREZGdsfTMbWZmJs6ePYucnBy5u25k69at8PPzQ0xMDEJCQgznnZyc0Lp1a9GrKjyKyS0RERGRvREkHFLCCgLeeust+Pr6Ijw8HD4+Pvjqq68sdRW4evUqmjVrBgBQKBRGz3l4eOD+/fuSYzK5JSIiIrIzlpq5Xbx4MdauXYujR4/i+vXr+Pnnn/Hee+/hzz//tMh1VK1aFUlJSX9f0yPJbUFBAbZv345atWpJjsnkloiIiMje6AXxhwTLly/HoEGD0KBBAwBAt27d0LJlS6xYscIil9GjRw+cO3cOH374IVJSUpCZmYk//vgD3bp1Q2ZmJnr27Ck5JpNbIiIiIntjgbIErVaL06dPo0WLFkbnQ0NDcfToUfmvAYCLiwt27NiBPXv2YPHixVi8eDHCw8MNSW6FChUkx+Q6t0RERER2RiFyg4biL/o1Go3ReZVKBZVKZXTu/v370Ov18Pb2Njrv7e0tekcxU9SpUwf79u3DrVu3kJqaCh8fHwQFBZkcjzO3RERERPZG4lJggYGB8PDwMBxRUVGPhXRycgKKZnAfpdVqDc/JSafTobCw0PDYz88PTZs2RVBQ0GPPSSE6udXr9ejdu7fRuV9++QVr1qwx6Y2JiKzlxIkTmDVrltG5OXPmWOxrNyIiuUm9oSw5ORlZWVmGIzIy8rGYXl5ecHFxwa1bt4zOp6amIjAwUPZr6NKlC/bt21fic4mJiWjXrp1JcSUlt1u3bjU6d/nyZZw9e9akNyYispa7d+/i4MGDRucOHTqEO3fuWK1PRESSSKy5dXd3Nzr+WZKAotUK2rdvb5Tv6XQ6xMbGolOnTrJ2/9ixY8jOzi41gX355Zfh4OCAhIQEybFZlkBERERkZxSCIPqQ4qOPPkJcXBw++OADxMXFYejQocjJycHUqVNl7f/x48fRpEmTJ7Zp3Lgxjh8/Ljm2pBvKCgoKEBERYXh85coV9OjRo9T2Wq3WqG7jn8XMRETWcvToUaPx7PDhw6W25VhGRDZHX3SIaSdBixYtEBcXhy+//BJxcXGoU6cO9u/fD39/f1N7WqL09HR4eno+sY2np6dJN7JJSm6VSiXGjRtnePzLL788sX1UVBRmz54tuVNERJYWGBhoNJ49qSSBYxkR2Rqxs7JSZ24BoGXLlmjZsqWJPRMnICAA27Zte2Kbo0ePYsCAAZJjSypLcHBwQEREhOGoV6/eE9tHRkYaFS8nJydL7iARkSX4+PgYjWfPPfdcqW05lhGRzbHQ9rtPS5cuXZCQkIBff/21xOe3bduG3bt3o2vXrpJjW3Sd25LWUCMisjccy4jI5jyyzFeZ7WxQ5cqV8cknn6BXr17o27cvunTpAn9/f9y6dQu7du3Cxo0bERUVhSpVqkiOLTq5VSqVj5Uh9O7dG/n5+ZLflIjImpo0aYI5c+YYnZs1axaqVatmtT4REUnx6DJfZbWzVe+88w4qV66M2bNn46effjKcDw4OxsqVKzFixAiT4opObhUKBbp162Z0rkaNGia9KRGRNVWuXBmVK1c2Ote8eXOr9YeISDI7n7ktNnToUAwdOhTJycnIyMiAl5eX2WvqcvtdIiIiIjuj0P99iGlnDwIDA2XbKILJLREREZG9KSczt5bA5JaIiIjI3ohdCeHZy22Z3BIRERHZG0uuc2vvnmpyKxT9BXN3H6JnT/HvvVAOBlqOZUTPNpsYz1iWUKqnmtxmZ2cDRUXDRPRsys7OhoeHh7W7YRaOZUQEK49nCr0AhU7EzK2eya1F+fv7Izk5GW5ublAoFKW202g0CAwMRHJyMtzd3WV7f0vEfdb7yuvn9YuNKQgCsrOzZd+f3BrEjmXgzwivn9dfLq/fJsYzQeSs7LOX2z7d5NbBwQEBAQGi27u7u8v6y2DJuM96X3n9vH4xMe19xraY1LEM/Bnh9fP6y931W308Y1lCqXhDGREREZG90QN48hdH/9/uGcPkloiIiMjOcLWE0tlkcqtSqTBr1iyoVCqbj/us95XXz+u3xPWXJ/wZ4fXz+p/d67coliWUSiGUh3V5iIiIiJ4BGo0GHh4e6FT/XTgqy07EC3Va7DwzH1lZWRape7ZFNjlzS0RERERPwJnbUjG5JSIiIrI3vKGsVExuiYiIiOwMbygrnc0lt7/99ht+++03AMAHH3yA559/Xpa4+/fvR0xMDLy9vTFixAj4+vqaHTM3NxfLly/H9evXERYWhvDwcFn6CgA3btzAJ598gqlTp6JBgwZmxSosLMS4ceMMj/38/PDJJ5/I0Evg2LFjiImJwa1bt9C9e3f07t3brHgxMTHYunWr0bnPPvsMzz33nFlxr1y5glWrVuHBgwfo1KkTIiIizIqHop1pVq5cievXr6NNmzbo06ePSXHS0tLwr3/9CwDQsWNHvP7664bnLl68iDVr1kClUmHs2LGS/h6ioqJw+fJlqNVq/Pe//y3zvBhxcXH47rvvAACTJk1CkyZNAAB6vR7r169HYmIiGjVqhCFDhsDJyUlS7PLoiy++wPnz56FUKrFs2TJZYgqCgB9++AEHDhxA/fr18eabb8pyA8yNGzewcuVK6HQ6DBkyBHXr1pWlvwCwY8cObNy4EUuXLoWjo3kfO6dPn8bXX39teNy+fXu88cYbMvQS2Lx5M+Li4qDRaPD222+bPfbOnTsX165dMzyuXLky5s2bZ3Y///rrL/z6669wcnLC4MGD8cILL5gd8+rVq1i9ejXy8vLQt29fNG/e3KQ4e/bswdq1awEA48ePR9OmTQ3PxcXFYcuWLahevTrGjh0r+mdBq9Vi4sSJAIAXX3zR8OfSzou1dOlSJCYmAgAWLVqEChUqAEU1rcuXL8etW7cQFhaGrl27Sor71LAsoVQO1u7APwUEBCAkJAT79+9HWlqaLDF/+uknfPDBB/Dz80NSUhLatGkjy37Q/fv3R3JyMvz9/TFmzBhDUm4unU6Hd955B8ePH0dycrLZ8QoLC7Fp0yaEhIQgJCQEjRs3lqWff/zxBzp06AC9Xo+QkBBZtiINDAw09LNatWqIi4uDj4+PWTG1Wi1atWqFvLw8BAUFYcKECfj999/N7mufPn1w8OBB1KpVC8uWLcPixYtNilOhQgWEhISgsLAQ+/fvN5xPT09H69atgaIEuF27dtDpdKLjNmjQAC1atMDq1atFnRfD19cXISEhOHnypNGH9pgxY7B9+3ZUr14dixYtwvTp0yXHLo/q16+PkJAQk/6uSzNlyhRs2rQJ1atXR3R0NCZNmmR2zPT0dPTp0wdqtRqCIKB169a4ffu2LP1NSUnBypUr8f3336OwsNDseMnJyTh79qxhnKhevbos/Zw1axYmT54MPz8/hISEyLJA/wsvvGDo58OHD3Hv3j2zYx45cgT9+vWDr68vVCoV2rdvb/a/VVZWFlq1agWtVgs/Pz+MGDECx44dMylW8Rhx+vRpXL161XB+x44dGDJkCHx9fbF582aMHz9edEylUomQkBC4ublhx44dZZ4Xq1atWggJCcGPP/4IrVZrON+2bVvcvn0bzz33HIYMGYJNmzZJjv1U6AXxx1OQmpqKmzdvypJfmU2wUZ06dRLi4+NliXXz5k1Br9cbHru5uQnZ2dlmx01LSzP8eebMmcKcOXPMjikIgjBr1ixh165dwiuvvCLExsaaHS83N1eoVKmSMGXKFGHu3LlG/TZH8+bNhW+++UaWWCX5+OOPhXnz5pkdJzU1VQgICDA8nj59uvDFF1+YHdfV1VW4d++eIAiCkJiYKNSuXduseCtWrBAmTpxoePzf//5XGDlypOFx+/bthT/++ENSzIKCAsHFxUX0ebEGDBgg/PLLL4bHycnJhj/v2rVLaN++vcmxyyOVSiVbrEf/rg8dOiQ0a9bM7Jg5OTlGY2JoaKiwf/9+s+PqdDph4MCBQlpamuDi4iLk5uaaHTM2NlZo0qSJMGnSJGHJkiWCVqs1O2ZGRoagVquFM2fOmB2rNK1btxYSExPNjrN+/Xph8ODBhsctWrQQDh8+bFbMxMRE4YUXXjA8nj9/vjBmzBizYg4ePFj48ccfDY/79u0rrF27VhCKft4qVaokPHjwQFLM2NhYoWfPnqLPi+Xr6yvcv3/f8PjR37GPPvpI+Pjjj02ObQlZWVkCAKFzjSlCeK3pZR6da0wRAAhZWVkW6c/KlSuFmjVrCn5+foKfn58QEBBg9PlgDTZXlmAJVatWNfx5/fr16NGjB1xdXc2O6+vri1mzZuHs2bO4e/cuvv/+e7NjxsXFQaVSoUOHDvjyyy/NjgcATk5O+Pzzz6HT6XDo0CE0b94cp06dMmtJEEEQkJSUhOeffx5TpkxBYGAgRo8eLdt2hIWFhVizZg0SEhLMjuXn54cxY8agY8eO8PHxgUajMXy1bo6BAweidevWaNq0KdLT03H9+nUIggCFQkyFf9muXLmChg0bGh43atQIV65ckSW23Iq3ohUEAdHR0Rg7dqy1u1RuPbrt74oVK2T5u65YsSI0Gg1GjRqFGzduGGaczfXZZ59h7NixspSBFWvYsCEmTpyIvLw8/Pjjj4iNjTX7W7Nz586hevXqOH/+PL755hs0a9YMb7zxBhwc5PlyMykpCbm5uUZf0ZuqV69eiImJQe/evfHw4UO0aNECzZo1Mytmw4YNIQgC2rdvDz8/P1y7dg0uLi5m9/VRj45nFStWRGBgIJKTk2Utf5FL8e9YVlYWtm3bhnXr1lm7S6UQWZYAy86knj17Ftu3b0eNGjUAAP/+978xcOBAnD59GjVr1rToe5fG5soSLGnt2rX4+eefER0dLVvMxo0bIzQ0FFlZWdizZ4/Z8UaMGIFLly5h1KhROHnyJBYsWIDjx4+bFVOpVGLUqFEYO3YsoqOj0bBhQ7P7qlAooFQqER0djdq1ayMxMdHkmtOSbNq0Cc2bNze71hYA7t+/j59//hmtW7dGy5YtcffuXRw8eNDsuCtWrMCyZcvQrVs3zJ07F2q1WrbEFkX/KSkoKDA8LigogLOzs2zx5VZYWIiRI0eiRYsWGDhwoLW7U67p9XpMnDgR1apVw6hRo2SJ6eTkhJCQELRq1Qr79u3DpUuXzIp37tw5REdHY926dRg1ahS0Wi3Gjx+P3Nxcs+IGBARg1KhRmDRpErZt24b4+HhoNBqzYjo4OOD27dvYsWMHateujUWLFuGzzz4zK+ajFi9eLNt/+I4fP44LFy6gVatWaNOmDfbs2WN2+ZpKpcLhw4cxbdo09O/fH8OHD5dlAuhR9jae3blzB6+99hq+/PJL1KpVy9rdKVlxza2Yw4Lmz59vSGwB4O2330Z+fr5Rmd3T9kzM3KLoJpqzZ89i48aNZt/QAAA5OTm4evWq4QYqX19fxMTEYMCAAWbFnTFjhuHPhw8fRr169VCpUiWz+1tMr9fj7t27hsJ5c9SvXx/vvfceXnrpJQwcOFC22jcUfRjMnDlTllgHDx6El5cX5syZAwBwdXXFDz/8gG7dupkdu02bNkDRzY/F9bFyadiwITZu3Ij33nsPOp0O8fHxGDZsmKzvIZcHDx5g4MCB6Nu3r832sbzIzc3FkCFD0LFjR0yYMEGWmOfOnYOfn58hUb5+/Tp2795t1oe6h4eH0Xi2bt06vPzyy1AqlbL0GQAyMzNRWFhodpJUu3ZtODk5YeHChVAqlfDz88OaNWtk6aNGo8Fvv/2G+fPnyxLv119/RZ8+fTBt2jSg6N9u165dePPNN82Kq1ar8eqrr6KwsBAdOnRAz549ZelvsYYNG+Kvv/5C8+bNkZKSgjt37hh9C2FLLly4gKFDh+Kbb77Biy++aO3ulE4viJuVfUo1t8VOnToFQRAQHBz8VN/3UTa3Q1lCQgJWrlyJ2NhYvPTSS2jSpInZd/avWLECU6dOxcCBAw0za1FRUWbdqJSXl4euXbvCw8MDzs7O2LNnD1atWiXLXfjFIiIiMGnSJLNXYTh8+DCWLVsGnU6HI0eOwNvbGzt27DA7yd+6dSumTJmCkJAQHDlyBJ07d8bChQvNiomirzh69uyJ8+fPyzITmpGRgcaNG6NevXrw9vbGH3/8gVWrVpk9eBev7HDt2jWcP38eO3bsMOkrNr1ejzFjxuD8+fPIyMhAaGgoZs2aBR8fHzRv3hy+vr7Izs6Gj4+PpBvhoqOjsX//fnz77bd488030aVLFwwYMKDU82IU36keFxeHGjVqoGHDhvj666/Rs2dPnDt3zpDsy7kihz1bs2YN4uPjsXr1agwbNkyWO/sHDx6MAwcOoGPHjgAALy8vfP7552bFTEpKwhtvvIH69esjIyMDZ86cwf79+2W5SbSYq6sr0tPToVarzYqzbt067N69G1qtFjt37sSIESMwd+5cs/s3c+ZMbNu2DXXq1MGff/6JhQsXon///mbHXbRoEc6cOYNvvvnG7FgAsHfvXkRERCAsLAy5ubk4ePAgDh8+bHYiERUVhYsXLyIxMRGenp7Yvn27Sf9WZ8+exZdffondu3ejWrVqqF+/PhYuXIhTp06hc+fOaN26NY4cOYKJEyfi3XffFR132rRpOH36NM6dO4fOnTtj8uTJaNSoUannxSgew9evX4++ffuiZcuWGDVqFKpUqYKGDRsaVmuSc0UOORTvUNY5aAIcHUTsUKbX4s8bS5CcnGxUjqhSqUpcaSU1NRV37tx5YszatWujYsWKj53Pzc1F69atUbFiRfz111+yfpsphc3N3Pr4+BjuLgUAT09Ps2M2b94c//nPf4zOmbt0jlqtxo4dO/Dnn38iNzcXixYtgp+fn5k9NTZ58mTUr1/f7Dje3t4ICQmBo6Mj3nzzTbRr106WH7ju3bujTp06OHz4MCZMmIDQ0FCzY+KRelu5fim8vLxw9uxZ7Ny5Ew8fPsTcuXNlqQMqXtmhV69eaNOmjcn1xgqFwuhnHgBcXFygVquRkJCAuLg4ODs7G5IZsWrWrAlBEAz/LsUffqWdF8PDw8Oor8UffiNHjjQaDOX4vS0PatasiYKCAsPflxzfbrzxxhvo0KGD4bEctZEvvPACdu3ahV27dsHFxQUdOnSQ5dudRy1evFiW5eGqV6+OvLw8qNVqvPfee6ITmbJ88sknePXVV3H58mXMmjVLtq+ia9eujR49esgSCwBat26NpKQk7N+/H05OTvj2229l+XavQYMG8PX1xdChQ9GmTRuTZ9j/OUYUf9Y2bNgQR48exb59+/DBBx9IrhNu2rQp6tWrZ3hcfM2lnRfj0dV5UPSzpVAoHluuTc5vJWUlcSmwf/5nddasWfj4448fa75x48Yyv7lYu3btY0vQ5efno2/fvtBoNPj999+tltjCFmduiYiIiKhkhpnbquPEz9ymLBU9c2uK4sT2zJkz2L17t9VLTmxu5paIiIiIyiBx5tbd3d2sVZJKU1BQgH79+tlMYgsmt0RERER2SBC5+5iFv58fOHAgdu/eje+++w7p6elIT08HAPj7+8uy4pEpmNwSERER2Rsb2X738uXLqF69umEb+WJTpkzB8OHDLfrepWFyS0RERGRv9HoAepHtLMfctfgtgcktERERkb2xkZlbW8TkloiIiMjeMLktFZNbIiIiIjsj6HQQBF3Z7fRltylvmNwSERER2RtBELe1LmduiYiIiMjmCYK4db6Y3BIRERGRzdPrAYWIlRAEy66WYIuY3BIRERHZG87clorJLREREZGdEfR6CCJmbgXO3BIRERGRzePMbamY3BIRERHZG70AKJjcloTJLREREZG9EQRx2+8yuSUiIiIiWyfoBQgiZm4FJrdEREREZPMEvciZW95QRkREREQ2jjO3pWNyS0RERGRnCgWtqFnZQhQ8lf7YEia3RERERHbC2dkZVapUwd60raJfU6VKFTg7O1u0X7ZEITyL89VEREREdiovLw/5+fmi2zs7O0OtVlu0T7aEyS0RERERlRsO1u4AEREREZFcmNwSERERUbnB5JaIiIiIyg0mt0RERERUbjC5JSIiIqJyg8ktEREREZUbTG6JiIiIqNz4P8I+AVHLp5IbAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "mt.pl.plate(cells, color=cells.var_names[0])" ] }, { "cell_type": "markdown", "id": "8f3c999e", "metadata": {}, "source": [ "## What treating cells as replicates costs\n", "\n", "`mt.metrics.diagnose_testing` relabels control cells as pseudo-treatments of whole control wells, where the\n", "truth is known to be null, and reports each hit caller's false positive rate under two nulls side by side:\n", "\n", "- the **cell-shuffle** null, which permutes single cells and so treats them as independent replicates, and\n", "- the **well-block** null, which permutes whole wells, the unit the experiment actually randomized.\n", "\n", "A calibrated test stays at or below the chance cutoff. Read the two rates against each other." ] }, { "cell_type": "code", "execution_count": 4, "id": "c9ca6ff1", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T01:08:45.290332Z", "iopub.status.busy": "2026-10-21T01:08:45.289839Z", "iopub.status.idle": "2026-10-21T01:08:46.414018Z", "shell.execute_reply": "2026-10-21T01:08:46.412978Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
checkvalueexpectedverdictnote
0hit_calling well-block null rate0 of 8<= 2passhit_calling called a control-only pseudo-treat...
1hit_calling cell-shuffle null rate8 of 8<= 0 (well-block)FAILthe same hit_calling null permuting single cel...
2edistance well-block null rate1 of 8<= 2passedistance called a control-only pseudo-treatme...
3edistance cell-shuffle null rate8 of 8<= 1 (well-block)FAILthe same edistance null permuting single cells...
\n", "
" ], "text/plain": [ " check value expected verdict \\\n", "0 hit_calling well-block null rate 0 of 8 <= 2 pass \n", "1 hit_calling cell-shuffle null rate 8 of 8 <= 0 (well-block) FAIL \n", "2 edistance well-block null rate 1 of 8 <= 2 pass \n", "3 edistance cell-shuffle null rate 8 of 8 <= 1 (well-block) FAIL \n", "\n", " note \n", "0 hit_calling called a control-only pseudo-treat... \n", "1 the same hit_calling null permuting single cel... \n", "2 edistance called a control-only pseudo-treatme... \n", "3 the same edistance null permuting single cells... " ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "report = mt.metrics.diagnose_testing(cells, n_draws=8, n_permutations=50)\n", "report" ] }, { "cell_type": "code", "execution_count": 5, "id": "a465bf61", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T01:08:46.417108Z", "iopub.status.busy": "2026-10-21T01:08:46.416959Z", "iopub.status.idle": "2026-10-21T01:08:46.605394Z", "shell.execute_reply": "2026-10-21T01:08:46.604368Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhEAAAFOCAYAAAA4rk+tAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjIsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvgI3uAAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAXDlJREFUeJzt3Xl8TFf/B/DPZJvsqzULQRAqEZUQhFiCpkJtEYRStfxaRWk1qLYopUW1pU/7UA9V1CNKaTWWEmoLsTf2WCOykn1fzu+Pyn1MMonJzUQmfN6vV16vzDn3nvs9Z+ZmvrnLuQohhAARERFRJenVdABERERUOzGJICIiIlmYRBAREZEsTCKIiIhIFiYRREREJAuTCCIiIpKFSQQRERHJwiSCiIiIZGESoSX/+te/4OjoiOzs7JoOpcrU9eV56p+2+fn5Yfz48U8te15oo29bt26Fr68vnJ2d4ejoKP1MmjRJa3HqAi8vL0ydOvWZbGvLli3o1q0bnJ2d0alTp0qt+zx/XkurbF9fpLGRQyeSiMTERHz66afo0aMHmjdvDnd3d/Tv3x/Lly/Hw4cPazo8yZdffglHR0fk5eWVqUtPT0dsbCyKi4trJDZtUteX56l/2hYfH4/k5OSnlmlT165da+wLt6p9O3PmDIYPH46BAwfi8OHDiIiIQEREBAwNDWXv7zU5HhVtOy4uDo8ePar2GE6ePIkRI0Zg0KBBOHz4MHbs2FGp9av786pL1PW1ovewNo/Ns9gvDKq1dQ3s2bMHw4cPR4sWLfDee+/Bw8MDRUVFOHbsGD7//HN89913iI6OrukwgSe+SNXNFD558mSMGjUKZmZmNRIb6ZYDBw7AwKD6dq+4uDjUr1+/2tqvSFX7FhYWBoVCgalTp0JfX18qf/L3yqrJ8ajJbZfYs2cP9PT0yowplaXu81vRe1jd+3J1ehafzRodmatXr2LIkCHw9fXFzp07YWhoKNW1bt0awcHBmDlzZk2GqDELCwtYWFjUdBikI2r6S6U6VbVvcXFxMDEx4ZedFiUkJHBMNVTZz+/zvC9rQ42ezliwYAFyc3Pxr3/9SyWBKGFqaopVq1aplF2/fh3jxo3DSy+9hObNm+O1117DwYMHVZaZOnUqvLy8yrQXGhoKR0dH3Lx5UypbvHgxHB0dUVxcjMWLF6Nt27Zo1aoVpk+frnL+f8aMGfjyyy8BAC4uLtI53MuXLwPlXDOgadslIiMjMWDAALi4uMDX1xdhYWHYsWMHHB0dce3ataeOZ1JSEubMmYMOHTqgefPmePXVV/H777+XWSYkJATt27dH06ZN0aNHD2zcuPGpbauTmpqKuXPnwtvbG82bN0fv3r3x/fffo7CwsML1St4HR0dHODk5oVWrVhg2bBiOHz+ulfZLxr2wsBALFy5Eu3bt8PLLL0v1x48fR1BQEFxdXdGiRQsEBwdL72OJknPZ0dHRCAwMRPPmzeHj44PNmzdrNDblnUc9ceIERo4cidatW8PNzQ0TJkzAnTt3pPoRI0ZIY9O4cWN4enri/fffVznM7+zsjDt37uCPP/6QlvXz81PZjiZ9PHv2LEaOHAl3d3e0adMGw4cPx8mTJ2X1rWS8YmJiEBwcjObNm6Nz5874+eefpWWys7Ph6OiIH3/8Ufrd0dER//rXv8rdlq6PhybbBlDhuFQ2zieVjOOGDRvUjqkm41ceTcejsjGXqOw+tm3bNvTt2xcuLi5o27Ytpk6diri4uErHXPrz+7T3sPTygwcPRteuXdXGOGnSJLRt2xZFRUVaG5/bt28jODgYLVq0wMKFC4FnvF88laghRUVFwsrKSnh4eGi8zrlz54SFhYXo2rWrOHz4sDh79qyYOHGiUCgUYt26ddJywcHBwsHBocz669atEwDElStXpLKQkBABQMycOVN8/fXX4tKlS2LTpk3CzMxMTJgwQVru0aNHYvr06QKAiI6OFjExMSImJkbk5+cLIYRYvHixACAyMjIq3bYQQpw8eVIYGxuLfv36iePHj4tTp06J4OBgMXbsWAFA/P333xWOze3bt4W9vb1wdXUVv/zyi7hy5YoICwsT/fv3F6dOnRJCCHH37l3h6OgoPDw8xO7du8WVK1fEd999J8zMzMTHH38staWuL+rK/Pz8RMuWLcWePXtEdHS0OHjwoJg8ebJYuXJlhbFmZWVJ43f37l0REREhxowZI4yMjMT58+er3H7JuE+ePFmsWLFCXLx4UXz66adCPP4M6Ovri7feekucPHlSnDlzRowZM0aYmZmJs2fPSm04ODiIV199VfTs2VPs2bNHXLx4UcyaNUsAEF999ZXK9l566SXx2muvPbVszZo1Qk9PT7z++uvir7/+ElFRUWLt2rWib9++0jJJSUnS2ERHR4tff/1VtG7dWnTo0EEUFBQIIYS4f/++cHZ2Fq+++qq0bHx8vNSGJn2Mjo4WZmZmYuzYseLMmTPi6tWrIjQ0VHTu3FkkJiZWOL7q+ubg4CACAgJEv379RFhYmPj777/FW2+9JQCIw4cPCyGEKC4uFjExMWLUqFHC1NRUij09PV0IIUSzZs3EkCFDVNrV9fF42rY1GZfKxFna08ZUk/FT955qOh5yYn5ybDTdx+bPny8UCoWYM2eOuHDhgti7d69o3bq1aNiwoYiJialUzKX7+rT3sPTy//73vwUAceLECZUYHz16JIyNjcXUqVO1Pj5//PGHOHXqlPQ996z2C03UWBKRlJQkAIjBgwdrvE737t1FnTp1RGZmpkp5nz59hJWVlfQFJyeJWL58ucqy77//vjA0NJR2RiGE+OSTTwQAkZOTU6btipIITdru1q2bcHJykpIS8fgPhLu7u0ZJRP/+/YWlpaVISkoqU1dUVCSEEGLQoEHCzs5OJCcnq9SvWLFCGBgYiHv37pXbl9JlWVlZQqFQiC+//LLc7VVWq1atxBtvvFHl9kvGfdGiRSrlSUlJwtTUVAQHB6uUFxcXC09PT9GzZ0+pzMHBQRgYGIhbt26pLDty5EhhamoqUlJSpDJNkoj4+HhhbGwshg0bVun+nD59WgAQBw4ckMrUfeFWpo+rV68WANR+QT4tnvKSCKVSKX2GhBCisLBQODg4iKFDh6osO2nSJGFmZlam3fL6VJqujUdFcWs6LpX5bKpT3piqo278Sr+nmoxHVWPWdB+7e/eu0NfXF+PHj1dZLiYmRhgZGYnXX39d45jV9VU85T0svXx6erowMzMrE8/KlSsFAHHhwgWtjY++vr64efNmhcuVqI79QhM1djqj5JC0pufwsrOzceTIEbz22mtlLl4cPXo00tLSyhwOr4zAwECV1+3bt0dBQYHKoebqajsrKwtHjx7FoEGDVE7rKBQKDBky5Knt5+TkYM+ePRg8eDDq1KlTpl5PTw+5ubnYvXs3+vXrBzs7O5X6gIAAFBYWljktVBETExM0aNAAa9aswd69e5Gfn6+yvYoUFxfjP//5D1555RU0b94cTk5OcHR0xK1bt3D9+vUqt19i2LBhKq//+OMPZGdnY8yYMSrlCoUC/fr1w19//aVy542XlxeaNGmisuzIkSORnZ2Nw4cPaxRDibCwMOTm5mLixIll6p7sT2JiIkJCQtCxY0c0adIEjo6O6N+/P/D4VN7TaNrHkn598MEHuHr1arnxVEbHjh3h5OQkvdbX14eHh4dGcZenNo9HCU3GpbKfTU3JHT9NxkMbMWuyj+3btw9FRUUIDg5WWc7R0RE9evRAWFiYxjFrg4WFBQIDA/Hf//4XWVlZUvl//vMfeHp6wt3dHdDS+Hh6eqJp06Zlyp/VfqGJGruw0s7ODkqlEg8ePNBo+aSkJBQVFcHR0bFMXUlZfHy87HgcHBxUXltbWwOAVm7PelrbycnJKC4uhr29fZl11ZWVlpSUhIKCAjRu3LjcZRITE5Gfn4+dO3fC2dlZusNECCGdv0tMTNS4TwqFAjt27MDUqVPxyiuvwNjYGJ07d0ZgYCDGjx9f4dXMM2fOxMqVK7F48WIsWbIEtra20NPTw9ChQ5Gbm1vl9ks8+YcbAO7fvw8AeOONN2BgYIDHR+IAABkZGSgsLMSjR4/QsGFDQM379mRZZT9rJZ/zit6jrKwsdOrUCUqlEosWLULr1q1hZmaGBw8eoGPHjtLYVETTPvr5+WHZsmVYunQp1q9fjwYNGqBXr154++230blz50r1rYS68bK2tsaFCxdktVfbx6OEJuNS2c+mJqoyfpqMhzZi1mQfS0hIAJ74O/8kR0dH7N27F0VFRdX6HpY2fvx4rF+/HqGhoRg7dizOnTuHc+fO4bvvvpOW0cb4lP4bhme8X2jyeauxJMLQ0BDdunXDkSNHkJqaKn2xlsfKygoA1F4QlJSUpLKMhYWF2gsXK/rDX16mqu52zsp6WtuWlpZAOX3T5AKokn6X7GzqlGwjKCgIH330UYXtaKpjx444efIk4uPjcezYMezatQtvvfUWIiMjsXbtWrXrFBUV4bvvvsOYMWPw3nvvqdQlJCSoHCWR0/6TlEqlyuuSMVi7di1eeukltevUq1dP+l3d2JfcL17ZsSr5fCckJMDFxUXtMn/88Qdu3bqFQ4cOwdfXVyq/e/euxtupTB/fe+89zJgxA1FRUThy5Ah++OEH+Pj44M8//0TPnj0r1T9Uwz5U28ejhCbjUtnPpiaqOn5PGw9txKzJPvbk3/7S+05SUhLMzc2lI9rV9R6W1qVLF7i6umLt2rUYO3Ys1q5dCxMTE4wYMUJaRhvjU/pvGGpgv3iaGr3FMyQkBPv378fChQuxbNkytcusWbMGEyZMgLW1NV566SXs378fxcXFKjvm3r17YWBgAG9vb+Dxf3spKSlISkpC3bp1peUOHDhQpXiNjY0BAAUFBdLv2mBjY4NWrVqpPZ0QHh7+1PWtrKzg5eWFPXv2IC8vT+0Hz9raGl5eXjhx4gQaNGig1fueGzRogCFDhmDIkCF4+PAhdu/eXe6yxcXFKCgoKHPa5fTp07hz506ZUy2Vbb8ifn5+UCgUOH36NPr27fvU5U+fPo20tDSVhGH//v3Q09Or9H82vXr1gkKhwK+//oouXbqoXabkP4jSYxMaGlpmWWNjYxQUFJQpr2wfFQoF3Nzc4ObmhhEjRsDW1hZ//PGHVv/gylVbxqO8bVdGZePURGXGrzwVjYc2YtZkH/Px8QEe/53v2LGjtFxmZiaOHTtW5k6JZ/Uevvnmm5g5cyYuXryIzZs3Y+jQoSr9qI73FDW0X1SkRm/x7NWrF5YsWYIvv/wS77zzDu7duyfVRUVFYdCgQVi0aJFUtmDBAly7dg0zZsxAdnY2ioqK8OOPP+LHH3/ElClTpPt5hw8fDkNDQ8yZMwcFBQXIz8/H559/DoVCUaV4W7RoATyeHU7bPv74Y5w5cwYLFixAfn4+CgoKsHLlSo0OTQHA8uXLERsbi+DgYOm2p9TUVCxevBgXL14EAHz11Ve4fv06xo4di9jYWODxtSkXL17ExIkTpTJN3Lp1C2+99Rb+/vtv6XRIdHQ0/v77b3h4eJS7nqGhIbp27YpNmzZJt9pevnwZM2bMkM4lVqX9iri6uuLdd9/FwoUL8cMPP0jn/DIzM/Hrr7/i/fffV1m+VatWeOutt5CamgohBH777TesWrUKb7zxBho1alSpbbds2RKTJ0/GN998g++++056Xy9fvixt18fHB0qlEp9//jny8vJQWFiIDRs2qD3H2aJFC0RFRZU54qZpH1evXo1vvvlGOoUlhMDOnTsBQPb4alttGY/ytl0Zlf1saqIy41eaJuOhjZg12cc8PDwwaNAgLFu2DH/88QeEEEhLS8O4ceOQnp6OTz75ROOYyyPnPXz99ddhaGiIESNGICUlBW+++aZKfXW8p3jG+4VGNL4EsxqFh4cLf39/YWpqKmxsbISFhYWws7MTgYGBIjw8XGXZ0NBQ4erqKgwMDISpqamoV6+eWLhwYZkrqNevXy/q1KkjjI2NRZ06dcSXX35Z4d0ZpYWFhQkAKtsvLCwUQ4YMEYaGhqJ+/frCwcFBXLp0SYin3J2hSdvi8dXF9erVE0ZGRsLW1lbMmzdPrF27VgAQV69efeo4Hjt2THTt2lXo6ekJa2trUadOHfHhhx+K7OxsaZkzZ86IV155RRgZGQlra2thamoq2rVrJ9asWSPdGqTJ3Rn5+fni3//+t/Dw8BCmpqbCzs5OWFhYiHHjxj31FsE7d+6IHj16CH19fWFlZSXc3d1FZGSk6NWrl2jfvn2V2y9v3Et8//33wtXVVejp6QlbW1thZWUlhgwZIiIiIqRlHBwcRHBwsNi1a5dwcnIS5ubmwtjYWEyePFnk5eWptKfpLZ7FxcVi+fLlolGjRlLf3dzcxK5du6Rltm3bJl3Rb2lpKUaNGiVu3LghAIgVK1ZIy124cEG4uLgIY2Nj4eDgIHr16lWpPt67d0+8++67okGDBsLKykqYmZmJZs2aPfX22fL6VjJepam7U6oyd2fUhvGoaNuVGRdN4ixPeWOq6fiVfk8rMx5yY67MPpaTkyPeffddYWVlJczNzYW+vr5o3769yt9QTWNW9/mt6D1Ut3yJwYMHCwDCxcWl3H5WdXzUeVb7hSZ0IokoUVRUJJKSksrcwqlOamqqSEpKEsXFxRW2l5iYKAoLC4V4Yn6CJ++PTktLk+4zflJubq6IiYkRubm5ZepycnJEbGysyjwR6enpIiYmRiUeOW2XxFwS46effioAiIcPHz51TJ6Mr/RtnKXl5+eLhIQEaWyepK4v6spK5OXliUePHmkcX4ns7GyRlpYmvU5KShIJCQlVbr+8cS8tIyOj3HF6cgcuLi4WSUlJZf6wlYiPjy/TjrqyJz169KjCz/mjR4+k7RUVFYmYmBiVpK5ESkqKuH//vsr935r2sURqaqra25bLo65vcXFxat+jR48eibi4uDIxx8bGllk2Li6u3M+5Lo9HRduuzLhUNs7S21Y3pk9ur6Lxq+jzqul4VDbmyuxjJUr+Pqp77zWNuaK+qnsPK1o+IyNDxMTEaPT3ubLjU95n50nPcr8oj0Jo48pBqjY+Pj54+PAhrly5UtOhvFAcHR3RvXt32bN5ElHFuI89H3TiKZ70z10LM2bMkG4FzMrKwocffohjx45h1qxZNR0eERFRGbXz0WTPIX19fTRq1AidO3dGWloa0tLS0KhRI6xevbrMhCBERES6gKczdFBmZiYKCwufOncGVZ/4+HgolUrY2NjUdChEzyXuY88HJhFEREQkC6+JICIiIlleyGsiiouL8eDBA1hYWFR5AioiIiJtE0IgIyMD9vb2Wn2AmLa9kEnEgwcP1D7YhIiISJfExMSoffiYrnghkwgLCwvg8ZtT8iASIiIiXZGeng4nJyfp+0pXvZBJRMkpDEtLSyYRRESks3T9lLvunmghIiIincYkgoiIiGRhEkFERESyvJDXRBARVYeioiIUFBTUdBhUCxgaGkJfX7+mw6gyJhFERFUkhEB8fDxSU1NrOhSqRaytrdGgQQOdv3iyIjWeRBQVFSEsLAxXr15FUFCQRvM3ZGdnY/fu3UhISICbmxt8fX2fSaxEROqUJBD16tWDqalprf5SoOonhEB2djYSExMBAA0bNqzpkGSr0SRi27ZtmDlzJpycnHDkyBF4eno+NYmIi4tDt27dYGZmhnbt2mHBggXo06cPn0lPRDWiqKhISiDs7OxqOhyqJUxMTAAAiYmJqFevXq09tVGjSUTdunVx6NAh6OvrazyDZEhICCwsLHDixAkolUpERUWhbdu2GDp0KAYOHFjtMRMRPankGghTU9OaDoVqmZLPTEFBQa1NImr07gxfX180btxY4+WLioqwfft2jBkzBkqlEgDQpk0b+Pj4YOvWrdUYKRFRxXgKgyrrefjM1Pg1EZURExODrKwstGzZUqXc1dUVp06dKne9vLw85OXlSa/T09OrNU4iIqpYQUEBsrOzYWVlVSvaJfVqVRKRkZEBPL6i9UnW1tZSnTqLFy/G/Pnzqy0u51m7q63tF9GdJf1qOgTSYbq2vzlY6GNej3rIN0mHwiC3psOpNHdHaw2W0r6dO3di/PjxWr+jpTLt5uXlwdDQsFJPySxvneLiYrX/oJqbm8PAoFZ91VZKrZpsquT8Uek3Ki0tDWZmZuWuN3v2bKSlpUk/MTEx1R4rERGVz8jIqMaOFpw/fx5eXl4wNzeHmZkZJkyYgNzcihPAp61z79492NjYoFGjRnB2dpZ+Dhw48Ax6VHNkJRGpqamYOXMmOnfuDFdX1zI/1aVx48ZQKpW4efOmSvnNmzfRvHnzctdTKpXSw7b40C0iIqAgPx+ZGf/7h+xpk2Tl5ORACFGmPD8/X+UfuydPHePx7Yz5+fll1vP398fFixfLbUddPCV3wqSmpiInJ6fCeMuTkZEBf39/tGvXDqmpqbhw4QL27duH9957TyvrnD17VooxNTUVffv2lRVnbSEriZgwYQJ27twJf39/vPPOO2V+tGnfvn1Yv349AMDAwAD9+vXDpk2bUFxcDAC4c+cODh8+jEGDBml1u0REz7M/ft2Ggd074tNPP4W9vT3MzMzg4eGBqKgoleV++uknNGrUCJaWlrCyssL06dNVvuA3b96Mli1b4r333oO9vT1MTU3h4eGBy5cvY/r06bCxsYG5uTl69+6NR48eSevt3LlT5cL6zZs3w9XVtcJ4IiMjpf/wbW1t4eTkhO+//75S/d6yZQtSUlKwfPlymJmZoUWLFpg1axb+85//IDMzs8rrCCHKJFLPM1knavbt24czZ87AxcWlShu/dOkSwsLCkJaWBgD473//i9OnT6Nz587o3LkzAGDr1q2IiIjA2LFjAQBffPEFOnfujD59+qBDhw7YsmULevTogeHDh1cpFiKiF01SYjzu3LmD69evQ19fH8HBwZg0aRKOHTsGANLf3vXr12PkyJG4cOEC+vXrBxsbG3z88cdSO/Hx8cjIyEB0dDQKCgrQq1cveHl5Ydy4cYiLi0N2djZ8fX2xZMkSfPHFF+XGExcXV2E83t7e0rUOQggcPHgQAwcOROvWrdGtWzeN+nzy5El4eHjAwsJCKvP19UVubi4uXLiALl26VGkdd3d3FBUVwd7eHu+++y6mTp1aqWsuahtZPbOyskKdOnWqvPGcnBzEx8cjJycH7733HszMzBAfH6+S2fXt2xdvvPGG9LpZs2aIiorCgAEDoFAo8Nlnn+GPP/6otffYEtHzKTc7C6nJiSo/WRn//MNUVFhYpi41OVFaNyM1pUxdft4/59/zcrLL1GWmpciK0cDAACtXroS5uTlMTEzwf//3f4iMjJSO9K5YsQKvvvoqRo8eDX19fbz88suYPXs2vvzyyzLtfP311zA1NYWVlRWCgoIghMDy5cthYmICOzs7DBkyBCdPnqxSPE/KyclB+/bt0bdvX+zatUvjPicmJpb5/qpXr55UJ3cdAwMDfPrpp7h//z6ysrKwfPlyzJkzB8uXL9c4ttpI1pGIoUOHYuXKlfjoo4+qtHFPT094enpWuExgYGCZsrp162Lq1KlV2jYRUXW6ffVvXD2n+qXp1KwlPLu/gpzsTITv/LnMOoPenAYAOPPXPqQkxavUtffti0Yuroi9fQMXThxSqavn0AhdXqn8KV0bu7oqk2RZWlpKt0iam5vj2rVrZSbx8/LyQlpaGhISElC/fv1/tl+vnjQDIwBYWFigYcOGMDIyUikrOepcnpJpw8uLJysrC9OnT0doaCgyMzNhZmaGnJycSp/OLp2UFBYWAk+Zt+Fp6zg6OmLu3LlSfUnStHLlSsycObNS8dUmspKIS5cuYd++fdiyZQtcXFzKDPyvv/6qrfiIiGqlJq5uaNioqUqZ4eNJ8kxMzdHjtRHlrtu+Wx8UFapeWGhq8c8F4Q5NmsO2nuqzFgwMDWXF+LTJjgwMDKQvyxIlr5+8bVFdO3ImUnraOh9//DFOnTqFU6dOSRfTv/7668jOztZ4Gw4ODmXmFSo5mmBvb6+1dfB4DqOYmBgUFhY+t7d5yupV8+bNK7wbgojoRWdsagZjU/W3nusbGMC6Tr1y17Wwtim3TmliCqXJs5li283NDcePH1cpO3r0KBo0aFAjzwk5c+YMBg0aJH3/FBUVISIiAu7u7hq34ePjgzVr1iA5OVk6RbF//35YWFhI7ZTM+VAyx4Mm66hz+vRp2NvbP7cJBOQmEatWrdJ+JEREpFM++OADtG3bFosWLcLo0aNx5swZfP7551i4cGGNxOPu7o6tW7di0KBBMDU1xZIlS3Djxo1KJRFDhw7F/Pnz8cYbb2DZsmW4e/cuFi1ahBkzZsDY2BgAcOvWLTRv3hy//fYbAgICNFpnwYIFMDMzwyuvvAIzMzP88ssvWL16NVasWFFt46ELqpweFRQUwFDmoTQiIqoZhoaGMC81Z46BgQGsrKyk0wqtWrXCnj17MG/ePHz99deoV68e5s+fr3Irv7pJo0rm5nmSsbGxyt0NpddT107peD799FOkp6fj1VdfhaGhIfr27Yvx48erzEPxtEmslEol/vzzT8yYMQPdu3eHhYUFZsyYgTlz5kjL6Ovrw8rKSvpu02Sdd999F0uXLsWIESOQmpqKFi1a4JdffsFrr71WwbtQ+ymEutlDnqKgoACLFi3C2rVrcf/+fWkCkrfffhszZsyo8q2f1S09PR1WVlZIS0vTysRTujYNb23Haa+pIrq2v5VMe13P3hEKAyMN1tAtNTXtNQG5ubm4ffs2mjRpIh3RKKHt76nqIusWz8WLFyM0NLTMrSvdunXDggULtBUbERER6TBZScT69evx888/Y9iwYSrl3bt3x++//66t2IiIiEiHyUoiYmNjpatjn7wlx8DAQPZ85kRERFS7yEoiWrRoIU1D+mQSsXHjxkpdJUtERES1l6y7M2bPno1Ro0ZJs3P98ssv2LNnD9atW4fQ0FBtx0hEREQ6SFYSMXLkSCgUCixcuBDFxcUYOnQoWrRogZ9++olP0yQiInpByJ4nYsSIERgxYgRycnJQXFwMMzP1M7MRERHR86nKk00ZGBhwsikiIqIXkKwLKwsKCjBv3jw4OTmpPKXt7bffRnR0tDbjIyIiIh3FyaaIiEiW8PBw9O7du9zXlbV7924MGDBAK7FduHABI0eORJcuXfDVV18BADZv3oyAgAB06tQJ3t7e8Pb2LvM9RpUj63TG+vXr8euvv8Ld3R1BQUFSeffu3fH2229rMz4iolppwKpjz3R7u97p8ky3BwAPHz5EZGRkua8rKykpCWfPntVKbK+99hoGDhyId955Bw4ODrh48SLGjBmDH3/8EU2b/vOI9pCQENy+fVsr23tRyUoiONkUERHpquTkZNy9excTJkzASy+9BDz+59fR0REjR46UlrOxKf+R66QZTjZFRPSCepichDlz5qBXr14YNmwYDhw4oFKfmJiIkJAQ9OzZE4MHD8batWsh45mNKvLz87Fs2TL069cPAQEB+Prrr1FYWKiyzIkTJzB27Fj4+flh1qxZyMzMlOo2b96M119/XWX5Y8eOwdfXFwDw119/oU+fPgCA4OBg6bTFggULEBcXB29v7wpPuVRHn59nnGyKiOgFlJyYgOD+vfBSK1e88847KCoqwpIlS+Dg4ABXV1ckJSXBy8sLffr0waxZs5CWloaPP/4YUVFRWLFiheztzp49G3v27MGnn34KMzMz7N+/H59++inmz58PPP4Sf+eddzBnzhwYGRnhgw8+wP3797Fx40YAwIMHD3Dx4kWVNlNSUnDy5EkAgJubG+bPn48BAwZg9uzZaNy4MQBgx44d2Lx5M7766isYGKj/6quuPj/PONkUEdELaO23K6A0NsGePXukL9UhQ4YgPz8fALBkyRK0bt0aa9askdZp0qQJOnbsiPnz58t+PPWxY8cwfvx4DB48GADQt29fldPgBQUF2LFjBxo1agQ8flz2m2++qXH7NjY2aNeuHQCgbdu2cHV1BQCcP38eSqUS3t7e5a5bXX1+nslKIiIiIjjZFBFRLXYm4hh8/V5R+a9coVBAqVQCj++0SE5Oho+PD4QQEEKgsLAQxcXFuH79Ojw9PWVtt2fPnli6dCkAoHfv3mjTpg1MTEyk+nr16kkJBAA4OjoiIyMDWVlZ1f49U119fp7JSiJ8fHykc1hPvvlERFQ75OTkwMzCotz6zMxM9O7dGxMmTChT17JlS4220bNnT2RnZwMAOnTogG+++QaLFi1Cu3bt8Ouvv2Lp0qVQKpVYvXq1dJ1C6ckLS667exbXJWijzy8aWUlEo0aNEB0dDRcXF+1HRERE1a5xk6a4ceVSufUuLi548OBBhYf/n+bzzz9HUVERAMDKygp4nBQEBgYiMDAQxcXFmDFjBiZNmoRbt25p1KaZmZmUmJSIj4+XHeOTtNHnF42suzPmz5+PcePG4fz589L5MyIiqj2GjnoD4Xt3Y9euXVLZb7/9Jn2Zv/POO9i7dy/Wr18v1aempuLTTz/VeBteXl7S3RGtWrUCHk9WmJSUBADQ09ODubl5pR6d4Obmhlu3buHChQvA47kpVq5cqfH6FdFGn180spKI119/HUeOHEG7du2gVCqhUChUfoiISLd17+2P9z/5DKNGjYKzszOcnJywdu1a1K1bFwDw6quvYvXq1Xj//ffh4OCAVq1aoWXLljA3N6/SduvUqSNd8NisWTNs2LAB33//vcbr+/j4YNy4cfD29oa7uzvc3d2lCymrqrr6/DxTCBknmv78888K6/38/KoSU7VLT0+HlZUV0tLStHK1rfOs3VqJi/5xZ0m/mg6BdJiu7W8OFvqY16Me6tk7QmFgpMEauqVFHWNcu3YN9evXR4MGDcrUFxYW4vr16zAwMECzZs2gr68v1T169Ai3bt2SLjgs/bo8xcXFuHHjBvT19dGkSROpzaSkJNy/f18lKcjMzERUVBQ6dOgAPb3//d8bFxeH1NRUuLi4ICsrC9evX0eHDh2Ax3NRnD17Fh4eHjA2NgYe3zr64MEDeHh4SG1cv34dRkZGcHZ21rjP2pSbm4vbt2+jSZMmUpwltP09VV1kJRG1HZMI3cYkgiqia/tbbU8i3B2tazqEF9bzkETIOp2Rl5en9mjEn3/+iby8PG3ERURERDpOVhIxd+5cnD9/vkz5uXPn8PHHH2sjLiIiItJxspKIn376qczc5QAwZswYaWpSIiIier7JSiKysrLUnrbIzc1Fenq6NuIiIiIiHScriejWrRsWLFggTSICAEVFRZg3bx58fHy0GR8RERHpKFkzVn7++efo2rUrDh8+jM6dO0MIgePHjyM5ORlHjhzRfpRERESkc2QdiWjTpg0uXLiAQYMGITY2FnFxcRg8eDAuXryINm3aaD9KIiIi0jmyjkTg8fMzPv/8c+l1Xl6e9PQ3IiIiev7JOhJx/fp1zJ07V3odEhICU1NTODs7IyoqSpvxERFRLXD58mUsWbKk3NfqnD9/Hl9++WWl2tUGTbZLmpGVRMyYMQNdu3YFAERHR2PVqlXYunUrAgICMHPmTG3HSEREOo5JxItJ1umMo0ePYsuWLQCAvXv3YsCAARgyZAh69uyJZs2aaTtGIqJax/2Hxs90exfH332m2yOC3CTC0NAQSUlJMDc3R1hYGF577TXg8UNLDAxkX2ZBRETPgBACq79ehn6DAuHu+M8DqTZu3IjExETMmDEDePzgq2XLlmHKlCmws7MDAISHhyMyMhLW1tbo1auX1v5pTElJQVhYGGJiYtClS5enThVQWFiIXbt24dq1a6hXrx4GDhwoxfikv/76CydPnoSdnR0GDhwIW1vbctvcsWMHbty4gWnTpvH6vkqQdTrD398fw4YNw9tvv41Dhw6hf//+wOMPWM+ePSvdXn5+PpKTk1GZZ4EVFhbi4cOHlVqHiIgAhUKB8H27cSDsN6ls5syZeP/995GUlAQAOHLkCJYvXw4rKysUFRVh0KBBmDRpEuLj43Hy5Em0b98e//3vf6scS1paGjw9PbF//35ER0fDz88PX3zxRbnL5+XlwcfHBx988AESExOxYcMGtGzZEpcuXZKWKSwsxMCBAzFkyBBER0fj+PHj8PX1RWxsrNo2ly9fjjfffBNdunRhAlFJsg4bfPvtt1iwYAHu3r2L7du3S4+P3b9/P+bNm6dxO8XFxZg5cya+++476OnpwcrKCqtWrcKgQYPKXefevXuYOHEiwsPDYW5ujvz8fEycOBFLly5VeUwsERGVz9PbB5ERRwEAV65cQU5ODjp27IjDhw9j6NChOHToEHx8fGBgYIBvvvkG169fx8WLF6WnTfr7+2PixIkYNGgQjIzkP700MzMTK1aswPjx4wEAfn5+GDt2LEaPHo2GDRuWWX7lypW4ffs2rl69ChsbGwgh0L9/f7z77rvYv38/AODrr7/GoUOHEBUVBUdHR+Dxo8MNDQ3LtDd79mysX78ehw8fhpubm+x+vKhkJREWFhZYunRpmfI1a9ZUqp3ly5dj3bp1OH78ONzd3bFq1SoEBQXhwoULaNWqldp1JkyYgPT0dMTHx8PGxgaRkZHw9fWFi4sL3nrrLTndISJ64Xh26oIdWzagqKgIhw4dQrdu3eDm5obw8HApiRg8eDAAYNu2bbCxscGyZcsghIAQAllZWUhJScGNGzfw0ksvyY5DX19f5VlMgYGBmDBhAg4dOoQRI0aUWT4sLAyBgYGwsbEBHh9VKUlm8vPzYWRkhO3bt2P48OFSAgGgTEJSXFyMCRMm4ODBgzh27BiaNm0quw8vshr91/3bb7/F+PHj4eHhAT09PUydOhWNGjXC6tWry13nxo0b6Nu3r/QB8vLyQvPmzREdHf0MIyciqt3ad+yM7KwsnDlzBocOHUL37t3RvXt3hIeHIzMzE2fPnkX37t2Bx//FK5VKFBYWoqioCMXFxTAxMcEnn3wCCwsLjbb32WefYd68eZg3bx42bNggldvY2KgcydDT00P9+vXx4MEDte08ePCgTEJgb2+P4uJixMXFAQASExNVEgh1UlNTsXnzZowePZoJRBXU2FWQiYmJuHv3Lrp06aJS7uPjg1OnTpW73tSpU/HNN9+ga9euaNy4Mfbt24fY2FiMHTv2GURNRPR8sLC0QsuX/jny8Ndff2HWrFlo2bIlbt68iW3btsHY2Bjt27cHANSrVw9OTk6VOl2tqZSUFOkIAh4fIUhISIC9vb3a5e3t7aVkoURcXBz09PSk5KJevXq4f/9+hdu1tbXFxo0bERAQABsbG0ybNk1rfXqR1NiRiJKLd+rUqaNSXqdOHalOnYkTJ6Jjx47o06cPOnTogOnTp2PRokUVnsvKy8tDenq6yg8R0YvO09sHa9asQV5eHtq2bQtTU1N4enpi4cKF0vUQADBkyBBs3boV169fV1n/+PHjGm9rzpw50pGIJ09fFBUVqRyZCA0NRUFBAXx9fdW24+/vj23btiEtLQ14fKfJmjVr0LNnTykRGTRoELZs2aJyNCMuLg7JyckqbXXv3h2//fYbPvzwQ3zzzTca94X+p8aORJRcBFlYWKhSXlBQAH19/XLXGzhwIHJzc5GYmAg7OzucPn0avXr1AgBMmjRJ7TqLFy/G/PnztRo/EVFt59mpC35a8y1ee+016W9y9+7d8dlnn2HChAnSclOnTsWJEyfg6emJwMBAWFpaSrd6/v7771WKwdzcHIsXL8axY8dgYGCAn376CZ988km5RyKmTJmCbdu2wdPTE/3798fZs2cRFRWFw4cPS8tMmzYNhw4dgoeHBwYPHoyioiJERERgz549Zdrr0aMHfvvtNwQEBEChUGDKlClV6s+LRuMkIjMzU+NGzc3Nn7qMg4MDACA+Pl6lPCEhQaorLT4+Hvv378euXbuke4I9PT0xdOhQrF+/vtwkYvbs2dK9zwCQnp4OJycnjftDRPQ88vTugk8++US69gEAgoODYWhoiKFDh0plBgYGCA0NRUREBI4fPw4DAwMEBQXB29tbWqZ169aYNWtWua/V8fDwwPz58zF27Fhpnoj9+/dLMyKra0epVOLo0aPYtWsXrl69ilGjRmHQoEEq80QYGhri999/x8GDBxEZGYn69evjiy++kK6l8/DwUPlO6NGjB/bu3YsDBw4gMTER9erVkzmiLx6F0HCiBYVCoXGjms7d0K5dO7Rv3x4//PAD8PioRMOGDTFlyhR8/PHHwBPny+rXr4+MjAxYWVlh48aNGDlypNROYGAgUlNTpdt7niY9PR1WVlZIS0uDpaWlxv0qj/Os3VVug/7nzpJ+NR0C6TBd298cLPQxr0c91LN3hMJA/q2ONcXd0bqmQ3hh5ebm4vbt22jSpIl062wJbX9PVReNj0ScOHFC6xv/6KOPEBQUBC8vL3Tq1AnLli0DAJVbNWfOnImIiAhERUXBwsICAwcOxNy5c2Fra4umTZti37592LFjB9atW6f1+IiIiKh8GicRTx620pbBgwdj48aN+Prrr/HFF1/Azc0Nhw8fRt26daVlbG1tpcmsAGDDhg1YunQp5s2bh4cPH6Jx48bYtGkTgoKCtB4fERERla/GrokoERQUVGECUHr6U3Nzc8yfP58XShIREdUwjZMITScUQSWuiSAiIqLaq0aviSAiIqLaq0aviSAiqu2KBQAIgEdgqZKKi4trOoQqkzXZVOlZv0orPQslEdHzKimrCCnZhTBPTYaJpQ0U+jU2h58subm5NR3CC0cIgfz8fCQlJUFPT69KT0GtabI+7U/ePaEOr4kgohdFoQCWHH2EEW4FcKufA329Gn2uYaUZ5ZjUdAgvLFNTUzRq1EiaLbQ2kpVEnDt3TuV1cXExbty4gQ8//BDTp0/XVmxERLXCo9xi/CsyDRZG6TAz0oOe5nPz1bgD73XXYCnSNn19fRgYGFRqIkddJCuJ8PDwKFP28ssvo0mTJpg2bRomT56sjdiIiGoNASA9XyA9v6imQ6mU0jMlElWGVo+htG7dGlFRUdpskoiIiHSU1pKIoqIirFq1qtwnrxEREdHzRdbpDHV3X6Snp8PQ0BCbNm3SRlxERESk42QlEatWrSpTZmNjg5dffvmpd24QERHR80FWEjF8+HDtR0JERES1SpVmRUlMTMSjR4/KlLu6ulalWSIiIqoFZCURkZGRCA4Oxo0bN9TWc7IpIiKi55+sJOL//u//0LlzZ/zyyy+wsbHRflRERESk82QlEVevXsXBgwdhZWWl/YiIiIioVpA1T0TTpk3x8OFD7UdDREREtYasJGL27NmYNGkSbty4wesfiIiIXlCyTmcEBwcDAFq0aKG2nokFERHR809WErF//37tR0JERES1iqwkws/PT/uREBERUa2i1ad4EhER0YuDSQQRERHJwiSCiIiIZGESQURERLIwiSAiIiJZZD/Fc8eOHQgNDcW9e/dQWFioUhcREaGN2IiIiEiHyToSsWrVKkyYMAH29vY4duwY/Pz8YGlpiZMnT8Ld3V37URIREZHOkZVErFy5EqGhoVi2bBkAYOHChdi3bx+WLFmChIQEbcdIREREOkhWEnHr1i106dIFAKBUKpGZmQkAmDRpEsLDw7UbIREREekkWUlEYWEhjIyMAABOTk44d+4cACA+Ph56erxWk4iI6EUg+8LKEqNHj8bw4cPh7++PAwcOYMCAAdqJjIiIiHSarCTi9u3b0u9z586FnZ0dTpw4gTfffBMzZszQZnxERESko2QlEc7OztLvenp6mDx5MiZPnqzNuIiIiEjH8QIGIiIikoVJBBEREcnCJIKIiIhkYRJBREREslTpFs/ExEQ8evSoTLmrq2tVmiUiIqJaQFYSERkZieDgYNy4cUNtvRCiqnERERGRjpOVRPzf//0fOnfujF9++QU2NjZVCuCvv/7Ct99+i4SEBLi5uWHOnDlo2LBhhetkZWXh22+/xcGDB2FqaooJEybA39+/SnEQERFR5chKIq5evYqDBw/CysqqShsPDw9Hnz59EBISgtdffx0rV65Ely5dcOHCBVhYWKhdJyMjA127doVSqcSsWbNgamqKb7/9Fra2tujYsWOV4iEiIiLNyUoimjZtiocPH1Y5iZg7dy6GDBmChQsXAgB69OiBhg0bYvXq1XjvvffUrjN//nwkJSXh6tWrUqLRt29f5ObmVikWIiIiqhxZd2fMnj0bkyZNwo0bN2Rf/5CdnY2IiAj0799fKjM1NYWfnx/+/PPPctfbuHEjRo0aVeZIhbGxsaw4iIiISB5ZRyKCg4MBAC1atFBbr0liERMTg+LiYtjb26uU29vb48CBA2rXefjwIRISEuDq6or3338fZ86cgb29PcaMGYM+ffqUu628vDzk5eVJr9PT058aHxEREVVMVhKxf//+Km+4oKAAAKBUKlXKTUxMpLrSSk5ZvP/++5g2bRo++ugjnDx5Ev369cO6deswatQotestXrwY8+fPr3LMRERE9D+ykgg/P78qb9jW1hYAyswz8fDhQ9jZ2aldx8bGBgqFAn369MHHH38MAOjZsyeuX7+Ob775ptwkYvbs2SpPF01PT4eTk1OV+0BERPQiq9KMlRkZGTh37hzOnj2LjIyMSq1rb2+PBg0aIDIyUqX85MmTaNeundp1TE1N0apVqzKnQBo2bIiUlJRyt6VUKmFpaanyQ0RERFUjK4nIy8vDtGnTYGdnh5dffhnt27eHnZ0dpk2bpnLtwdOMGzcOP/zwA2JjYwEAv/zyCy5fvoxx48ZJyyxevFi6BgMAJk2ahG3btiE+Ph4AkJSUhK1bt2rl6AgRERFpTtbpjJCQEOzatQsbNmyAt7c3FAoFTpw4gVmzZkGhUOCrr77SqJ2PP/4YN27cgIuLC5ycnHD//n2sWrUKXl5e0jI3b97EhQsXpNfvvPMOrl27hhYtWsDZ2Rm3bt1Cnz598MUXX8jpChEREcmkEDLu0axfvz52794NT09PlfLIyEgEBAQgISGhUu09ePAACQkJcHFxKXPr5q1bt5CZmQl3d3eV8uTkZNy/fx9OTk7lXkNRnvT0dFhZWSEtLU0rpzacZ+2uchv0P3eW9KvpEEiHcX/TLu5vuknb31PVRdaRiNTUVDRr1qxMuYuLC9LS0irdnr29fZnrHEo0bdpUbXmdOnVQp06dSm+LiIiItEPWNREeHh5YsWJFmfLly5ejbdu22oiLiIiIdJysIxFLly6Fv78/tm/fjg4dOgAATp06hdu3byMsLEzbMRIREZEOknUkolu3brh27RoCAgKQlJSE5ORkBAQE4Nq1a+jWrZv2oyQiIiKdI+tIBAA4OjpiyZIl2o2GiIiIao0qTTZFRERELy4mEURERCQLkwgiIiKShUkEERERySI7iSgqKsKZM2cQGhoqlWVlZWkrLiIiItJxspKI2NhYdOjQAd7e3hg2bJhUPnz4cM4TQURE9IKQlURMnz4dbm5uZaa4DgkJwWeffaat2IiIiEiHyZonIjw8HJcuXYKpqalKedu2bREZGamt2IiIiEiHyToSkZOTA0NDQwCAQqGQyhMSEsokFkRERPR8kpVE+Pj4YO3atcATSURubi5mzZqFHj16aDdCIiIi0kmyTmcsW7YMvr6+2LNnD4QQGD16NA4dOoScnBwcO3ZM+1ESERGRzpF1JKJNmza4ePEiOnbsiN69eyM+Ph4jR47E+fPn0bJlS+1HSURERDpH1pGIsWPHYv369Vi0aFG5dURERPR8k3Uk4scff1RbLoTAhg0bqhoTERER1QKVOhKRm5ur9ncAKC4uxtGjR1G/fn3tRUdEREQ6q1JJhImJidrfSygUCixcuFA7kREREZFOq1QSceTIEQBA165dpd9LGBoaolGjRmjYsKF2IyQiIiKdVKkkwsfHBwBw5coVuLq6VldMREREVAvIurCSCQQRERHJusUTAHbs2IHQ0FDcu3cPhYWFKnURERHaiI2IiIh0mKwjEatWrcKECRNgb2+PY8eOwc/PD5aWljh58iTc3d21HyURERHpHFlJxMqVKxEaGoply5YBABYuXIh9+/ZhyZIlSEhI0HaMREREpINkJRG3bt1Cly5dAABKpRKZmZkAgEmTJiE8PFy7ERIREZFOkpVEFBYWwsjICADg5OSEc+fOAQDi4+OhpyerSSIiIqplZF9YWWL06NEYPnw4/P39ceDAAQwYMEA7kREREZFOk5VE3L59W/p97ty5sLOzw4kTJ/Dmm29ixowZ2oyPiIiIdJSsJCI+Ph7Ozs4AAD09PUyePBmTJ0/WdmxERESkw2RdwFAycyURERG9uGQlEY0aNUJ0dLT2oyEiIqJaQ1YSMX/+fIwbNw7nz59Hfn6+9qMiIiIinSfrmojXX38dANCuXTu19UKIqkVFREREOk9WErF//37tR0JERES1iqwkws/PT/uREBERUa3C6SWJiIhIFiYRREREJEuVp72uquzsbOzevRsJCQlwc3ODr6+vxuteu3YNv/32G9q3b48ePXpUa5xERESkqkaPRMTFxaFt27ZYtGgRzpw5g8DAQIwaNUqjdXNycjBkyBDMmzcPO3furPZYiYiISJXsJKKoqAhnzpxBaGioVJaVlVWpNkJCQmBhYYGTJ09i3bp1OHjwIH7++Wf8+uuvT1136tSp6NWrF1xcXGTFT0RERFUjK4mIjY1Fhw4d4O3tjWHDhknlw4cPR1hYmEZtFBUVYfv27RgzZgyUSiUAoE2bNvDx8cHWrVsrXHfr1q2IiIjA559/Lid8IiIi0gJZScT06dPh5uaGtLQ0lfKQkBB89tlnGrURExODrKwstGzZUqXc1dUVV65cKXe927dvY8qUKdi0aROMjY012lZeXh7S09NVfoiIiKhqZF1YGR4ejkuXLsHU1FSlvG3btoiMjNSojYyMDACAtbW1Srm1tbVUV1pBQQFGjBiBWbNmwd3dXeN4Fy9ejPnz52u8PBERET2drCMROTk5MDQ0BAAoFAqpPCEhoUxiUZ6S5UofFUhLS4OZmZnadTZu3IirV6+isLAQy5Ytw7Jly5CUlISzZ89i2bJl5U63PXv2bKSlpUk/MTExGveViIiI1JP9KPC1a9cCTyQRubm5mDVrlsa3WjZu3BhKpRI3b95UKb958yaaN2+udp0WLVpg/PjxSEhIQHx8POLj41FQUICsrCzEx8eXm0QolUpYWlqq/BAREVHVyDqdsWzZMvj6+mLPnj0QQmD06NE4dOgQcnJycOzYMc02bGCAfv36YdOmTZg0aRL09PRw584dHD58GOvWrZOW27dvHx48eICxY8eiS5cu6NKli0o7f/75J7p27Yply5bJ6QoRERHJJOtIRJs2bXDx4kV07NgRvXv3Rnx8PEaOHInz58+XuVCyIl988QVu3LiBPn36YM6cOejZsyd69OiB4cOHS8ts3bqVCQIREZEOknUk4sqVK2jVqhUWLVpUpY03a9YMUVFR+Pnnn5GQkIDPPvsMgYGB0NfXl5bp27cvWrVqVW4bY8aMQbNmzaoUBxEREVWerCSidevWePnllzFq1CiMGDECDRo0kB1A3bp1MXXq1HLrAwMDK1x/+vTpsrdNRERE8sk6nXHhwgX4+flhxYoVcHR0RN++ffHTTz8hMzNT+xESERGRTpKVRLi7u+Pzzz/H3bt3ceDAATRu3BjTpk1D/fr1ERwcrP0oiYiISOdU6QFcCoUCvr6+WL16Nfbv34/mzZtj8+bN2ouOiIiIdFaVkoiYmBh88cUXaNu2LTw9PWFkZISvv/5ae9ERERGRzpJ1YeWaNWuwadMm/PXXX2jWrBlGjhyJbdu2lTtJFBERET1/ZCURH374IYKCgrBkyRJ4e3trPyoiIiLSebKSiAcPHsDAQNaqRERE9JzQOBMouX3T3Nwcubm5FS5rbm5e9ciIiIhIp2mcRFhYWAAAhBDS7+Up70FYRERE9PzQOIk4ceKE2t+JiIjoxaRxEvHkBZTff/891q9fr3a5sWPH8mJLIiKiF4CseSJ+/PFHteVCCGzYsKGqMREREVEtUKlbLJ68oLL0xZXFxcU4evQo6tevr73oiIiISGdVKokwMTFR+3sJhUKBhQsXaicyIiIi0mmVSiKOHDkCAOjatav0ewlDQ0M0atQIDRs21G6EREREpJMqlUT4+PgAAK5cuQJXV9fqiomIiIhqAVmTTTk6Okqv1eFkU0RERM8/TjZFREREsnCyKSIiIpJF1mRT6iaTysvLg1Kp1F5kREREpNNkTTZ1/fp1zJ07V3odEhICU1NTODs7IyoqSpvxERERkY6SlUTMmDEDXbt2BQBER0dj1apV2Lp1KwICAjBz5kxtx0hEREQ6qFK3eJY4evQotmzZAgDYu3cvBgwYgCFDhqBnz55o1qyZtmMkIiIiHSTrSIShoSGSkpIAAGFhYejZsycAoLCwEAYGsvISIiIiqmVkfeP7+/tj2LBh8PLywqFDh/DDDz8AAMLDw6WEgoiIiJ5vso5EfPvtt+jevTuSk5Oxfft2NGjQAACwf/9+zJs3T9sxEhERkQ6SdSTCwsICS5cuLVO+Zs0abcREREREtUCVLmDIyMhAdHQ0hBBo3rz5U2eyJCIioueHrNMZeXl5mDZtGuzs7PDyyy+jffv2sLOzw7Rp05CXl6f9KImIiEjnyDoSERISgl27dmHDhg3w9vaGQqHAiRMnMGvWLCgUCnz11Vfaj5SIiIh0iqwk4ueff8bu3bvh6ekplTVu3BjNmjVDQEAAkwgiIqIXgKzTGampqWonlXJxcUFaWpo24iIiIiIdJyuJ8PDwwIoVK8qUL1++HG3bttVGXERERKTjZJ3OWLp0Kfz9/bF9+3Z06NABAHDq1Cncvn0bYWFh2o6RiIiIdJCsIxHdunXDtWvXEBAQgKSkJCQnJyMgIADXrl1Dt27dtB8lERER6RzZ80Q4OjpiyZIl2o2GiIiIag3ZSURubi5CQ0Nx5coVAEDr1q0RGBgIpVKpzfiIiIhIR8k6nXH+/Hm4uLhgypQpCA8PR3h4OCZPnozmzZvj4sWL2o+SiIiIdI6sJGLSpEnw8/NDbGwsTpw4gRMnTuDBgwfo2bMnJk6cqP0oiYiISOfIOp1x8eJF/P777zAzM5PKzMzM8MUXX6Bx48aVbi8/Px/p6emws7ODQqHQaJ2srCwYGBjw9AkREVENkXUkokmTJkhOTi5TnpycjCZNmmjcTnFxMd577z1YW1vD2dkZjo6O2LFjR4XrbNq0Ce3atUPDhg1hbW2Njh07IjIyUk43iIiIqApkJREffPABhg8fjoMHDyI9PR3p6ek4ePAgRowYgZCQEI3bWb58OdatW4fjx48jPT0dISEhCAoKki7WLK2oqAhhYWFYv3490tLSkJqairZt28Lf3x+PHj2S0xUiIiKSSSGEEJVeScNTDgBQUfPOzs4YNmwYvvjiC6nMxcUF/fv3VzsjpjoPHjyAg4MD9uzZg759+2q0Tnp6OqysrJCWlgZLS0uN1qmI86zdVW6D/ufOkn41HQLpMO5v2sX9TTdp+3uqusi6JmL//v1V3nBiYiLu3r2LLl26qJT7+Pjg1KlTGrcTHR0NAGjQoEGVYyIiIiLNyUoi/Pz8qrzhpKQkAECdOnVUyuvUqYPjx49r1EZ2djamTp2KHj16VPjMjry8POTl5Umv09PTZcdNRERE/5B1TYRWNqz3z6YLCwtVygsKCqCvr//U9fPz8zF06FBkZ2dj8+bNFS67ePFiWFlZST9OTk5VjJ6IiIhqLIlwcHAAAMTHx6uUJyQkSHXlKUkgrl+/joMHDz71VMbs2bORlpYm/cTExGihB0RERC+2GksiLC0t4eHhoXJ9RWFhIQ4cOKDyEK+UlBQkJCRIrwsKChAYGIjLly/j0KFDcHR0fOq2lEolLC0tVX6IiIioamQ/O0MbPvroIwQFBcHLywudOnXCsmXLAABvvfWWtMzMmTMRERGBqKgoFBcXY/jw4YiIiMDOnTsBAPfv3wcA2NrawtTUtIZ6QkRE9OKRfSSiqKgIZ86cQWhoqFSWlZVVqTYGDx6MjRs34scff8SgQYOQnp6Ow4cPo27dutIytra20umKtLQ0nDx5EoaGhhg6dCi8vb2ln5KkgoiIiJ4NWUciYmNjMWDAAFy8eBGFhYXSXBDDhw/H22+/DX9/f43bCgoKQlBQULn1T84hYWNjIx15ICIiopol60jE9OnT4ebmhrS0NJXykJAQfPbZZ9qKjYiIiHSYrCMR4eHhuHTpUplrENq2bcvnWBAREb0gZB2JyMnJgaGhIVBqCuyEhARe3EhERPSCkJVE+Pj4YO3atcATSURubi5mzZqFHj16aDdCIiIi0kmyTmcsW7YMvr6+2LNnD4QQGD16NA4dOoScnBwcO3ZM+1ESERGRzpF1JKJNmza4ePEiOnbsiN69eyM+Ph4jR47E+fPn0bJlS+1HSURERDpH9mRTDg4OWLRokfQ6Ly8PSqVSW3ERERGRjpN1JOL69euYO3eu9DokJASmpqZwdnZGVFSUNuMjIiIiHSUriZgxYwa6du0KAIiOjsaqVauwdetWBAQEYObMmdqOkYiIiHSQrNMZR48exZYtWwAAe/fuxYABAzBkyBD07NkTzZo103aMREREpINkHYkwNDREUlISACAsLAw9e/YEHj+F08CgRp/pRURERM+IrG98f39/DBs2DF5eXjh06BB++OEH4PFMliUJBRERET3fZB2J+Pbbb9G9e3ckJydj+/bt0lM29+/fj3nz5mk7RiIiItJBso5EWFhYYOnSpWXK16xZo42YiIjoWZlnVdMRPD/mpWmw0PNF4yQiMzNT40bNzc3lxkNERES1hMZJhIWFhcaNCiHkxkNERES1hMZJxIkTJ6o3EiIiIqpVNE4ivL29qzcSIiIiqlVk3Z1BREREJHtmqB07diA0NBT37t1DYWGhSl1ERIQ2YiMiIiIdJutIxKpVqzBhwgTY29vj2LFj8PPzg6WlJU6ePAl3d3ftR0lEREQ6R1YSsXLlSoSGhmLZsmUAgIULF2Lfvn1YsmQJEhIStB0jERER6SBZScStW7fQpUsXAIBSqZTmkJg0aRLCw8O1GyERERHpJFlJRGFhIYyMjAAATk5OOHfuHAAgPj4eenq8VpOIiOhFUOVHbo4ePRrDhw+Hv78/Dhw4gAEDBmgnMiIiItJpspKI27dvS7/PnTsXdnZ2OHHiBN58803MmDFDm/ERERGRjtI4iRg6dCi2bdsGADh9+jScnZ0BAHp6epg8eTImT55cfVESERGRztH4AoZff/0VRUVFAIDAwMDqjImIiIhqAY2PRLi4uGDRokXo1KkTAODPP/8sd1k/Pz/tREdEREQ6S+MkYuXKlZgyZQoWLFgAAOjdu3e5y/IpnkRERM8/jZOI3r174+rVqwAAhULBRIGIiOgFJ+vujMjISO1HUgPi4+ORlZUlvTY2NoaNjQ0KCwuRlJRUZvmGDRsCAJKTk1FQUCCV2ymykCGUyIcBlCiAuSJfZb0C6CNdGEMBAVtFdpl2HwkTCOjBQpELIxSp1GUJI+TCEEYohIUiT6WuEHpIEyYAAFtFNhRQTexShQmKoAczRR6Mofp8kxxhiGwYwQBFsFLkqtQVQ4EUYQoAsFFkQ69Uu2nCGIXQhynyYaIoUKnLhQGyhBL6KIa1IkelTkCBR4/btVLkwADFKvUlY5iZmYmMjAyVOqVSCVtbWxQVFSExMbHMGNavXx96enp4+PAh8vNVx9/S0hJmZmbIyclBamqqSp2hoSHq1KkDAIiLiyvTbt26dWFgYICUlBTk5qqOk7m5OSwsLJCXl4dHjx6p1Onr66NevXoAgISEBBQXq/bVzs4ORkZGSE9PV/kMAoCpqSmsrKxQUFCA5ORklTqFQoEGDRoAAJKSkso8t8bGxgbGxsZqx7Dk813eGDZo0AAKhULtGFpZWcHU1BTZ2dlIS0tTqTMyMoKdnR2EEIiPjy/Tbr169aCvr692DC0sLGBubo7c3FykpKSo1BkYGKBu3brA43215B8XO8U/45UqjFEE/Sp9vq0VOdAv9TlMF8YogD5MkA/TUp/vPBggs5zPNwA8FGYAAEtFLgxL7cu6+jcCAFJhgRyYqNSZIRuWyEQeDPEINip1eihGffzz2UxAHRSXurzOFilQogDpMEcWTFXqTJADa2SgAAZIhm2pngo0xD9/e5Ngi8JSX1HWSIMJ8pAJU2TAXKVOiTzYIg1F0EMi6pQZw/pIgh4EHsIa+TBSqbNEBsyQgxwYIxWWKnWGKEAd/PPZjEO9Mu3WxUMYoAgpsERuqb8hVfkbYWhoWGZbukhWEuHp6an9SGrAunXrYGxsLL12c3PD4MGDkZ6ejtWrV5dZ/pNPPgEA7Ny5E/fv35fKBxgDh/Ob4FaRHZrop6CT0T2V9WKLLLEvvwUMUIwBxlfKtLs5py3yoIcOhjFopK/6B/pUgSMuFTaAvV46eihvqdQ9LDbFrrzWAIAA5RXoK1T/QOzIfQmpwgQeBnFoYaD6ZXSxoAHOFDqijl4W/JXXVeqyhCG25rYFAPRW3oBZqT+kYXktEF9siVYGiXA3VP3SuF5YB8cKnGGhyCvT1yKhwIbc9gAAX6PbsNNT/WMZntcUd4pt8ffff2Pfvn0qdS1atMCIESOQm5ur9r2ZNWsWlEolwsLCcPPmTZU6f39/dOjQATdu3MCOHTtU6hwdHfHmm28CgNp2p0yZAltbW4SHh+Pvv/9WqfP19UX37t0RExODTZs2qdTZ2Nhg6tSpAIANGzYgO1u1r+PGjYOTkxNOnDhR5oF1np6e6NevH5KTk8vEZGRkhNmzZwMAQkNDyyS7w4cPR8uWLXHu3DkcPHhQpa5169YIDAxEVlaW2r5++OGHMDAwwG+//Ya7d++q1PXv3x8vv/wyrl69it9++02lrnHjxhg7diyKiorUtjt9+nRYWlrizz//xOXLl1Xqevbsia5du+Lu3bvYsmWLSl3dunXx9ttvA4/31ZLEZsDjXXZXbis8FGZwM4hHKwPVcYgqrI/IAifY6OUgQHlVpS5XGODnXA8AQC+jaFjqqX7x7s1rjgfFVmhpkIR2hqpfCjcLbfFXQVOYKvLV7svrcv7529jV6Dbq6akmh7r6NwIA/oI3zincVOq6iFPww1HEoT5+VAxTqbMQGZiBNQCATRiEDIWFSv0YsRXOuI9T8MAxRQeVunbibwzAfqTACqsVo1Tq9EUh5uIbAMB2+CNeUV+lfqj4DS/hBv6GK/YpuqvUtRA3MQI7kQtlmXYBYJZYBSXyEYaeuKlwVqnzFwfQARdwA02wQ+GvUucoHuBN/PPZVNfuFPEf2CIV4eiCv0t9/qvyN6K23MCgEC/geYn09HRYWVnh2rVrsLD434df7pGIgG+O6Ox/GbXxSETUXF8eieCRiHKPRAR8cwTgkQhJVf9G3DQexSMR2joSMemUSl1Vj0TUrVsXaWlpsLRUjUmXvNBJhLbeHOdZu7USF/3jzpJ+NR0C6TDub9p1x3hkTYfw/JiXpsFCmtH291R14YMuiIiISBYmEURERCQLkwgiIiKShUkEERERycIkgoiIiGSp8STir7/+QlBQELp3744pU6aovc1OG+sQERGRdtVoEhEeHo5evXqhefPmmDlzJm7cuIEuXbqUub+9qusQERGR9tVoEjF37lwMGTIECxcuRL9+/bB9+3Y8fPhQ7ax3VVmHiIiItK/Gkojs7GxERESgf//+UpmpqSn8/PzKfcy4nHWIiIioesh6doY2xMTEoLi4GPb29irl9vb2OHDggNbWAYC8vDzk5f1vOtiSaXvT09Or2It/FOeVnaaW5NPW+0LPJ+5v2pWueOEmLa4+WvzbVfJ3UNcnla6xJKLk2RNKpVKl3MTEROW5FFVdBwAWL16M+fPnlyl3cnKSFTtVL6uvajoCoheHVU0H8DxZov3RzMjIgJWV7r5LNZZE2Nr+8+CV0g8lefjwIezs7LS2DgDMnj0bM2bMkF4XFxfj0aNHsLOzg0KhqFI/SLvS09Ph5OSEmJgYnZ4vnuh5wP1NdwkhkJGRUebIu66psSTC3t4eDRo0QGRkJAICAqTykydPomvXrlpbB4+PXJQ+emFtba2VflD1sLS05B81omeE+5tu0uUjECVq9O6McePG4YcffkBsbCwA4JdffsHly5cxbtw4aZnFixcjODi4UusQERFR9auxIxEA8PHHH+PGjRtwcXGBk5MT7t+/j1WrVsHLy0ta5ubNm7hw4UKl1iEiIqLqV6NJhFKpxNatW/HgwQMkJCTAxcUFFhYWKsvMmTMHmZmZlVqHai+lUolPPvmkzOknItI+7m9UVQqh6/ePEBERkU6q8WdnEBERUe3EJIKIiIhkYRJByM/Px5YtW5CamlruMsXFxdiyZQuSk5OfaWwVxfC010TPo9jYWGzfvl16fffuXezcubNGY6IXF5MIQnp6OkaMGIE7d+6Uu0x+fj5GjBiBq1evAgAKCwuxZcsWPHz48JnFWTqGp70meh5FRkaq3NJ+5MgRTJo0SaN1b9++jV27dlVjdPSiYRJBGtHX10dQUBDq1q0LAMjNzcWIESNw48aNmg5NUjpGoheBs7MzBg4cqNGyhw8fxttvv13tMdGLo0Zv8STdc+fOHVy6dAkNGzbEyy+/LJXr6+tj4MCB0vTiO3bsAAD8+eefuHPnDmxsbNC3b9+ntl/yJNb8/Hx4e3urzBy6fft25OfnQ19fH05OTmjXrl2lbj0rHWN+fj62b9+OV155BampqWr7VeLmzZuIiopCkyZN4ObmhtDQUPj6+qJ+/foab59Ijnv37uHcuXOwtbXFyy+/DDMzM5V6IQQiIiKQkpICd3f3Mus7OTnB399fpezRo0c4ffo09PX14eXlBUtLSzx48AAnT55ETk4OtmzZAgBwd3eHnZ0dwsPDgcfPIWrRogVatWql0t7Nmzdx9epVvPLKK7h48SJiY2Ph4eEBR0fHMvHcv38fZ8+eRb169eDl5QV9ff1K9ZdqGUEvvKSkJAFA9O/fX7Rs2VIEBAQIKysr8cYbb0jL5OTkCADiyJEjQgghRowYIQAIPz8/ERQUJEJCQp66nd9//13Y2tqKNm3aiL59+4pmzZqJQ4cOSfXjxo0TQUFBYujQocLV1VU0b95c3Lhxo9wYnvZak34JIcTy5cuFkZGR6NGjh/Dw8BB9+/YVenp6Yv/+/VUeW6KKTJ8+XdjY2IhXX31VeHt7CwcHBxERESHV5+fni1deeUXY2dkJf39/4ejoKPr06SOsrKykZX766SdRv3596fWuXbuEpaWl8PX1FX379hVNmzYVe/fuFWfPnhUdOnQQJiYmIigoSAQFBYnt27eLv//+W3odEBAgbG1tRVBQkCguLpba/O6770TdunVFp06dRLdu3UT37t2FsbGx2L59u7RMcXGxePfdd4WJiYnw9fUVnTp1El26dBGpqaka95dqHyYRJH3ZBgUFiaKiIiGEEKdPnxYAxOXLl4VQ8wWdkZEhAIgTJ05otI2YmBhhYmIiFi5cKJU9fPhQaq+04uJi8frrr4uhQ4dKZXKTiIr6defOHWFkZCR+/vlnaTuTJk0SAJhEULX6z3/+I5ydnUViYqJUtmTJEuHi4iK9/vbbb0XdunVFbGysEI/3mSZNmlSYRHh7e4sFCxZIr1NSUsTBgweFEEKsW7dOODg4VBhXcnKycHR0FFu2bJHKvvvuOwFApWzWrFmidevW0ut///vfwsTERJw9e1Yqi4iIEHFxcRr3l2ofns4gycSJE6Gn989lMu3bt4e5uTmuXbtW5tCmHKGhoTAzM8OsWbOkMltbW/j4+Kgsd+3aNURHRyMjIwO2trY4dOhQlbddUb927dqF+vXrY/jw4dLyM2fOxL///e8qb5eoIuvWrYO7uzsOHz6Mx//QwczMDNHR0YiJiYGTkxO2bt2K4OBg6UmOtra2mDBhAj7//PNy2zUxMcGNGzeQk5MDExMTWFtbo0ePHhXGUlhYiNOnT+PBgwfIz8+Hk5MTTp06haCgIGkZGxsbldfdu3fH8uXLIYSAQqHAjz/+iODgYLRr105apmPHjpXqL9U+TCJIUvKo9RJKpRK5ublaafvevXto2rRpmfOjJQoKCjB48GAcPXoUXl5esLa2Rnx8PBITE6u87Yr6FRMTg8aNG6vUl35NVB3u3LmDzMxMbNu2TaU8KCgIBQUFwOP9ZtCgQSr1TZo0qbDdr7/+GhMmTEDdunXRuXNn9O/fHxMnTiz3+qJLly7B398fSqUSrq6uMDMzQ1xcXJl9T91+VFBQgKKiIhgYGODevXsVXuCpSX+p9mESQc+EtbV1hbeDhoaG4tSpU7h9+7Z0seX69etx6tSpao3L1ta2zPwYFc2XQaQtlpaW6N69O7788styl7Gzs0NKSopKWenXpbm5uSEiIgLx8fE4cOAAFixYgCNHjmDr1q1ql//oo4/g4+ODzZs3S2UBAQGo7BMRnraPa9Jfqn14iyfJYmJiAn19fY2PVPTp0wc3b97EsWPHVMqTkpIAAPHx8ahfv77K3Rq//PKLlqMuq0uXLrh8+bLKraqcuIeehVdeeQU///wz0tLSVMpjY2Ol3318fLBr1y6VL/QnJ5pSp2T9Bg0aIDg4GO+++y4iIiIAAObm5mX22fj4eLi6ukqvExMTceTIkUr3p0+fPti6dSvy8vKksuzsbGRlZWncX6p9eCSCZNHX14eHhweWL1+O2NhY1KlTp8JbPDt16oS3334b/fr1w9SpU+Hg4IDdu3cjICAAEydOhL+/P2bPno233noLnp6e+P3333H06NFq70fXrl3h7++PV155Be+++y5SUlLwww8/QKFQQKFQVPv26cU1d+5c7N+/H15eXpg4cSLMzMxw6tQpXLt2DcePHwcAfPDBB9i4cSMCAgLw2muv4cCBAzh//nyF7QYGBqJFixbo1KkTcnNzsWLFCgQGBgIA2rVrh9TUVMyfPx8tW7aEu7s7Bg4ciCVLlsDU1BRKpRKrVq2SriGqjA8//BB//PEHvL298cYbbyAvLw8///wzdu3aBTMzM436S7UPj0QQlEolgoKCYGNjo1I+ePBg6foAdRM5/fLLL2jTpg3CwsKk+8wr8u2332Ljxo1ITk7GhQsXMGnSJEycOBEA0KpVKxw/fhxGRkY4cuQIevTogd9++w3Dhg2T1i8dw9Nea9IvANi2bRveeecdnD17FkII/PHHHxBC8BHzVK2sra1x8uRJhISE4OrVq4iKikLPnj1x+PBhaZmGDRsiMjISrVu3xunTp9G9e3fs2rULQ4YMkZYpPdnU4cOH0bNnT5w7dw7R0dFYvnw5li5dCgBo1qwZwsLCkJiYiJ07d+LatWuYOXMmvvnmG1y9ehVXrlzB8uXL8dlnn6FTp05Smy4uLggICFCJv169eggKCpISDltbW2k2zXPnziElJQX//e9/pbkkNOkv1T58FDi98B49eqRy0diWLVvwxhtvICkpCebm5jUaGxGRLuPpDNKKy5cv4+LFi2rrbG1t0adPn2cek6Y+++wzZGRkwMvLCzdv3sTKlSsxe/ZsJhBERE/BIxGkFTt27MB///tftXXNmjXDokWLnnlMmiosLMRPP/2EiIgIWFpaom/fvvDz86vpsIiIdB6TCCIiIpKFF1YSERGRLEwiiIiISBYmEURERCQLkwgiIiKShUkEERERycIkgoiIiGRhEkFERESyMIkgIiIiWZhEEBERkSz/DxeYHP9mRaa3AAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "rates = report[report[\"check\"].str.endswith(\"null rate\")].copy()\n", "rates[\"caller\"] = rates[\"check\"].str.split().str[0]\n", "rates[\"kind\"] = np.where(rates[\"check\"].str.contains(\"cell-shuffle\"), \"cell-shuffle\", \"well-block\")\n", "rates[\"rate\"] = rates[\"value\"].apply(lambda v: int(v.split()[0]) / int(v.split()[-1]))\n", "\n", "callers = [\"hit_calling\", \"edistance\"]\n", "x = np.arange(len(callers))\n", "width = 0.38\n", "fig, ax = plt.subplots(figsize=(5.5, 3.5))\n", "for offset, kind in ((-0.5, \"cell-shuffle\"), (0.5, \"well-block\")):\n", " heights = [rates[(rates[\"caller\"] == c) & (rates[\"kind\"] == kind)][\"rate\"].iloc[0] for c in callers]\n", " ax.bar(x + offset * width, heights, width, label=kind)\n", "ax.axhline(0.05, linestyle=\"--\", color=\"grey\", linewidth=1, label=\"nominal 0.05\")\n", "ax.set_xticks(x)\n", "ax.set_xticklabels(callers)\n", "ax.set_ylabel(\"false positive rate on a null screen\")\n", "ax.set_title(\"Counting cells as replicates inflates the false positive rate\")\n", "ax.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "7cc34218", "metadata": {}, "source": [ "The tall cell-shuffle bars are the pitfall: counting cells as replicates makes pure noise look significant\n", "in 8 of 8 draws for both callers, and the verdict fails. The well-block bars, the rates to trust, stay within the cutoff: 0 of 8 for `hit_calling` and 1 of 8 for `edistance`.\n", "The gap between the two is the pseudoreplication, measured on this screen rather than asserted. Each row of the\n", "table names the block it draws whole and, for `hit_calling`, the method it measured, so the check describes the\n", "test you mean to run." ] }, { "cell_type": "markdown", "id": "eb4f9fdf", "metadata": {}, "source": [ "## The default already draws whole wells\n", "\n", "At cell resolution `tl.hit_calling` and `tl.edistance` resample whole wells for their permutation null without\n", "being asked, so the calibrated well-block rate above is what a plain call gives you. On this small screen the\n", "caller goes further and warns that its null has few wells and is therefore only conservative, pointing to\n", "`tl.aggregate` for an exact test: the tool says where it stands rather than returning a confident wrong\n", "answer." ] }, { "cell_type": "code", "execution_count": 6, "id": "d239c642", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T01:08:46.608414Z", "iopub.status.busy": "2026-10-21T01:08:46.608228Z", "iopub.status.idle": "2026-10-21T01:08:47.807358Z", "shell.execute_reply": "2026-10-21T01:08:47.806109Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "UserWarning: the well-block null draws from 4 block(s) in the null half, too few to calibrate finely, so the test there is conservative rather than anti-conservative. For an exact well-level test, aggregate to wells with mt.tl.aggregate before testing.\n" ] }, { "data": { "text/html": [ " \n", " \n", " " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "mt.tl.hit_calling(cells)\n", "mt.pl.hits(cells)" ] }, { "cell_type": "markdown", "id": "cd01e3b9", "metadata": {}, "source": [ "## The other route: aggregate first\n", "\n", "When the question is about the treatment rather than the cells, aggregate to the well with `tl.aggregate` and\n", "test the wells with `tl.differential_features`, the pseudobulk approach single-cell settled on. The unit of\n", "the test is then the unit that was randomized. That path is worked through in\n", "[Which measurements moved](../phenotypes/which_features_moved.ipynb)." ] }, { "cell_type": "code", "execution_count": 7, "id": "ed6119bf", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T01:08:47.810441Z", "iopub.status.busy": "2026-10-21T01:08:47.810081Z", "iopub.status.idle": "2026-10-21T01:08:48.156648Z", "shell.execute_reply": "2026-10-21T01:08:48.155492Z" } }, "outputs": [ { "data": { "text/plain": [ "AnnData object with n_obs × n_vars = 96 × 10\n", " obs: 'Metadata_Plate', 'Metadata_Well', 'Metadata_CellCount', 'Metadata_Row', 'Metadata_Col', 'Metadata_Perturbation', 'Metadata_Control', 'Metadata_Control_Type', 'Metadata_Batch', 'Metadata_ImageNumber', 'Metadata_ImageID', 'Metadata_ObjectType', 'hits_distance', 'hits_qvalue', 'hits_reference_held_out'\n", " var: 'object', 'feature_group', 'feature', 'channel', 'scale', 'angle', 'gray_levels', 'radial_bin', 'params', 'is_feature', 'feature_kind'\n", " uns: 'mantispy'\n", " layers: None (.X)" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "wells = mt.tl.aggregate(cells)\n", "wells" ] }, { "cell_type": "markdown", "id": "f5bf1a56", "metadata": {}, "source": [ "## The same pitfall on a real screen\n", "\n", "The synthetic plate makes the point with strong, deliberate gradients. The effect is not an artifact of those\n", "gradients. `mt.ds.jump_cells` is one real JUMP plate at single-cell resolution, with eight DMSO control wells\n", "sitting at different positions on the plate. We take a slice of it for speed, relabel its control cells as\n", "pseudo-treatments, and run the same check. The direction is the same as on the synthetic plate, though subtler,\n", "because the wells of a single plate vary less than the gradients we injected above." ] }, { "cell_type": "code", "execution_count": 8, "id": "bd4e0d2b", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T01:08:48.159595Z", "iopub.status.busy": "2026-10-21T01:08:48.159433Z", "iopub.status.idle": "2026-10-21T01:08:49.542048Z", "shell.execute_reply": "2026-10-21T01:08:49.540943Z" } }, "outputs": [ { "data": { "text/plain": [ "AnnData object with n_obs × n_vars = 3600 × 120\n", " obs: 'Metadata_ImageNumber', 'Metadata_ObjectNumber', 'Metadata_CellsObjectNumber', 'Metadata_NucleiObjectNumber', 'Metadata_AbsPositionZ', 'Metadata_AbsTime', 'Metadata_BinningX', 'Metadata_BinningY', 'Metadata_ChannelID', 'Metadata_ChannelName', 'Metadata_Col', 'Metadata_ExposureTime', 'Metadata_FieldID', 'Metadata_ImageResolutionX', 'Metadata_ImageResolutionY', 'Metadata_ImageSizeX', 'Metadata_ImageSizeY', 'Metadata_MainEmissionWavelength', 'Metadata_MainExcitationWavelength', 'Metadata_MaxIntensity', 'Metadata_ObjectiveMagnification', 'Metadata_ObjectiveNA', 'Metadata_PlaneID', 'Metadata_Plate', 'Metadata_PositionX', 'Metadata_PositionY', 'Metadata_PositionZ', 'Metadata_Row', 'Metadata_Site', 'Metadata_Well', 'Metadata_Center_X', 'Metadata_Center_Y', 'Metadata_ObjectType', 'Metadata_ImageID', 'Metadata_Source', 'Metadata_JCP2022', 'Metadata_InChIKey', 'Metadata_Perturbation', 'Metadata_Perturbation_Type', 'Metadata_Control_Type', 'Metadata_Control'\n", " var: 'object', 'feature_group', 'feature', 'channel', 'scale', 'angle', 'gray_levels', 'radial_bin', 'params', 'is_feature', 'feature_kind', 'selected'\n", " uns: 'mantispy'\n", " layers: None (.X)" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "jump = mt.ds.jump_cells()\n", "# A slice for speed: at most 150 cells per well and the 120 highest-variance features.\n", "jump = mt.pp.downsample(jump, n_per_group=150, by=\"Metadata_Well\", seed=0)\n", "variance = np.asarray(jump.X).var(axis=0)\n", "jump = jump[:, np.sort(np.argsort(variance)[::-1][:120])].copy()\n", "jump" ] }, { "cell_type": "code", "execution_count": 9, "id": "b8013e93", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T01:08:49.545241Z", "iopub.status.busy": "2026-10-21T01:08:49.545076Z", "iopub.status.idle": "2026-10-21T01:09:30.582206Z", "shell.execute_reply": "2026-10-21T01:09:30.581333Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
checkvalueexpectedverdictnote
0hit_calling well-block null rate0 of 8<= 2passhit_calling called a control-only pseudo-treat...
1hit_calling cell-shuffle null rate5 of 8<= 0 (well-block)FAILthe same hit_calling null permuting single cel...
2edistance well-block null rate0 of 8<= 2passedistance called a control-only pseudo-treatme...
3edistance cell-shuffle null rate2 of 8<= 0 (well-block)passthe same edistance null permuting single cells...
\n", "
" ], "text/plain": [ " check value expected verdict \\\n", "0 hit_calling well-block null rate 0 of 8 <= 2 pass \n", "1 hit_calling cell-shuffle null rate 5 of 8 <= 0 (well-block) FAIL \n", "2 edistance well-block null rate 0 of 8 <= 2 pass \n", "3 edistance cell-shuffle null rate 2 of 8 <= 0 (well-block) pass \n", "\n", " note \n", "0 hit_calling called a control-only pseudo-treat... \n", "1 the same hit_calling null permuting single cel... \n", "2 edistance called a control-only pseudo-treatme... \n", "3 the same edistance null permuting single cells... " ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mt.metrics.diagnose_testing(jump, n_draws=8, n_permutations=100)" ] }, { "cell_type": "markdown", "id": "cb1ce25a", "metadata": {}, "source": [ "The same signature appears on real data: `hit_calling` under the cell-shuffle null over-calls the control-only\n", "pseudo-treatments and its verdict fails, while the well-block null stays clean. `edistance` shows the same\n", "direction more mildly, 2 of 8 draws, which still passes. The inflation is smaller than on the synthetic plate, as expected for one plate whose\n", "wells differ only by real position and seeding rather than by a gradient we set, but the lesson holds: counting\n", "cells as replicates manufactures significance, and the well-block null does not." ] }, { "cell_type": "markdown", "id": "205e80b1", "metadata": {}, "source": [ "## Takeaway\n", "\n", "- Cells within a well are not independent replicates; the well, or the plate, is the unit that was\n", " randomized.\n", "- Testing at cell resolution is fine as long as the null resamples whole wells. `tl.hit_calling` and\n", " `tl.edistance` do this by default.\n", "- `mt.metrics.diagnose_testing` measures, on your own screen, whether the null is calibrated and what a\n", " cell-shuffle would cost, and names the exact method and block it checked.\n", "- When the question is about the treatment, aggregate to wells with `tl.aggregate` and test with\n", " `tl.differential_features`, which refuses cell-level input for this reason." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.13.15" } }, "nbformat": 4, "nbformat_minor": 5 }