{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "53b76130",
   "metadata": {},
   "source": [
    "# Differential Splicing Analysis with pylimma (pasilla dataset)\n",
    "\n",
    "This tutorial walks through a differential splicing workflow -\n",
    "testing which **exons** are differentially used between conditions\n",
    "while accounting for overall gene-level expression. The underlying\n",
    "test is a per-exon interaction evaluated via `diff_splice`.\n",
    "\n",
    "\n",
    "**R-parity validation**: the sibling notebook [`pasilla_R_vs_Python.ipynb`](pasilla_R_vs_Python.ipynb) runs R limma 3.66.0 on the same input and compares the outputs numerically.\n",
    "\n",
    "## Dataset\n",
    "\n",
    "**Brooks et al. 2011, *Genome Research* 21:193-202**: *Drosophila melanogaster* S2 cells with RNAi knockdown of the splicing regulator *pasilla* (Ps) vs untreated controls. The shipped subset is the exon-level count matrix from the `pasilla` Bioconductor package.\n",
    "\n",
    "## Pipeline\n",
    "\n",
    "Differential splicing differs from the standard DE workflow:\n",
    "\n",
    "- **Input**: exon-level counts (not gene-level).\n",
    "- **Output**: `diff_splice` + `top_splice` replace `eBayes` +\n",
    "  `top_table`.\n",
    "\n",
    "Steps:\n",
    "\n",
    "1. Load exon counts + metadata\n",
    "2. Library-size QC\n",
    "3. log-count density (KDE, per sample)\n",
    "4. log2 transform\n",
    "5. Design matrix\n",
    "6. Fit per-exon linear models\n",
    "7. `diff_splice` + `top_splice`\n",
    "8. Top-ranked gene's exon-level logFC\n",
    "9. Summary\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "751d58d7",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-04-20T19:13:13.190327Z",
     "iopub.status.busy": "2026-04-20T19:13:13.190228Z",
     "iopub.status.idle": "2026-04-20T19:13:14.507655Z",
     "shell.execute_reply": "2026-04-20T19:13:14.506924Z"
    }
   },
   "outputs": [],
   "source": [
    "import sys\n",
    "from pathlib import Path\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "REPO = Path.cwd().resolve().parents[1]\n",
    "sys.path.insert(0, str(REPO))\n",
    "sys.path.insert(0, str(REPO / 'data'))\n",
    "\n",
    "import generate_data as gd\n",
    "\n",
    "import pylimma\n",
    "\n",
    "pd.set_option('display.width', 120)\n",
    "pd.set_option('display.max_columns', 10)\n",
    "\n",
    "# Shared plotting helpers used by every tutorial.\n",
    "\n",
    "def _log_cpm(mat, lib_size):\n",
    "    cpm = mat.div(lib_size, axis=1) * 1e6 if hasattr(mat, 'div') \\\n",
    "          else (mat / lib_size) * 1e6\n",
    "    return np.log2(cpm + 1e-2)\n",
    "\n",
    "\n",
    "def plot_log_cpm_density(mat, ax, title=\"\"):\n",
    "    '''Kernel-smoothed log-CPM density, one line per sample.\n",
    "\n",
    "    Matches edgeR/limma's standard density plot (limma User's Guide\n",
    "    Figure 15.1): a Gaussian KDE per sample drawn on a shared grid.\n",
    "    '''\n",
    "    from scipy.stats import gaussian_kde\n",
    "    arr = mat.values if hasattr(mat, 'values') else np.asarray(mat)\n",
    "    # Shared evaluation grid spanning the combined range.\n",
    "    lo, hi = np.nanmin(arr), np.nanmax(arr)\n",
    "    grid = np.linspace(lo, hi, 200)\n",
    "    for j in range(arr.shape[1]):\n",
    "        col = arr[:, j]\n",
    "        col = col[np.isfinite(col)]\n",
    "        if col.size < 2:\n",
    "            continue\n",
    "        kde = gaussian_kde(col)\n",
    "        ax.plot(grid, kde(grid), linewidth=0.8, alpha=0.7)\n",
    "    ax.set_xlabel('log2 CPM')\n",
    "    ax.set_ylabel('Density')\n",
    "    ax.set_title(title)\n",
    "\n",
    "\n",
    "def mds_coords(mat, top=500):\n",
    "    '''Pairwise leading-logFC MDS, matching R limma's plotMDS.default.\n",
    "\n",
    "    Returns the (n_samples, ndim) coordinate matrix and the percent of\n",
    "    variance explained per dimension. We reuse pylimma's private\n",
    "    _mds_coordinates for consistency with plot_mds().\n",
    "    '''\n",
    "    from pylimma.plotting import _mds_coordinates\n",
    "    r = _mds_coordinates(np.asarray(mat, dtype=np.float64), top=top,\n",
    "                         gene_selection='pairwise', ndim=2)\n",
    "    lam = np.maximum(r['eigen_values'], 0.0)\n",
    "    coords = r['eigen_vectors'] * np.sqrt(lam)\n",
    "    return coords, r['var_explained']\n",
    "\n",
    "\n",
    "def plot_mds_coloured(mat, groups, ax, top=500, title='MDS'):\n",
    "    '''MDS scatter plot coloured by sample group.'''\n",
    "    coords, var_exp = mds_coords(mat, top=top)\n",
    "    pct = np.round(var_exp * 100).astype(int)\n",
    "    groups = np.asarray(groups)\n",
    "    levels = sorted(pd.unique(groups))\n",
    "    palette = plt.get_cmap('tab10').colors\n",
    "    for k, level in enumerate(levels):\n",
    "        m = groups == level\n",
    "        ax.scatter(coords[m, 0], coords[m, 1],\n",
    "                   s=55, color=palette[k % len(palette)],\n",
    "                   label=str(level), edgecolor='k', linewidth=0.3)\n",
    "    ax.axhline(0, color='grey', linewidth=0.4, linestyle=':')\n",
    "    ax.axvline(0, color='grey', linewidth=0.4, linestyle=':')\n",
    "    ax.set_xlabel(f'Leading logFC dim 1 ({pct[0]}%)')\n",
    "    ax.set_ylabel(f'Leading logFC dim 2 ({pct[1]}%)')\n",
    "    ax.set_title(title)\n",
    "    ax.legend(fontsize=8, frameon=False)\n",
    "\n",
    "\n",
    "def plot_heatmap(E, groups, fit, n_top=50, ax=None):\n",
    "    '''Heatmap of top-n DE rows showing BOTH directions.\n",
    "\n",
    "    Picks the top n/2 genes with the most positive t and the top n/2\n",
    "    with the most negative t (by smallest p within each side), then\n",
    "    stacks them - up-regulated in the contrast at top, down-regulated\n",
    "    at bottom. Row z-scored.'''\n",
    "    t = np.asarray(fit['t']).ravel()\n",
    "    p = np.asarray(fit['p_value']).ravel()\n",
    "    half = n_top // 2\n",
    "    up_pool   = np.where(t > 0)[0]\n",
    "    down_pool = np.where(t < 0)[0]\n",
    "    top_up    = up_pool[np.argsort(p[up_pool])[:half]]\n",
    "    top_down  = down_pool[np.argsort(p[down_pool])[:n_top - half]]\n",
    "    # Sort each block by signed t descending so most-up is top and\n",
    "    # most-down is bottom.\n",
    "    top_up   = top_up[np.argsort(-t[top_up])]\n",
    "    top_down = top_down[np.argsort(-t[top_down])]\n",
    "    ordered  = np.concatenate([top_up, top_down])\n",
    "\n",
    "    mat = E[ordered]\n",
    "    col_order = np.argsort(groups)\n",
    "    mat_sorted = mat[:, col_order]\n",
    "    z = (mat_sorted - mat_sorted.mean(axis=1, keepdims=True)) / \\\n",
    "        (mat_sorted.std(axis=1, keepdims=True) + 1e-8)\n",
    "\n",
    "    if ax is None:\n",
    "        fig, ax = plt.subplots(figsize=(9, 8))\n",
    "    im = ax.imshow(z, aspect='auto', cmap='RdBu_r', vmin=-2, vmax=2)\n",
    "    ax.set_xticks(range(len(col_order)))\n",
    "    ax.set_xticklabels(np.asarray(groups)[col_order], rotation=90, fontsize=6)\n",
    "    ax.set_yticks([])\n",
    "    ax.set_title(f'Top {n_top} DE genes (top {half} up + top {n_top - half} down; row z-scored)')\n",
    "    return ax, im, ordered\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb040651",
   "metadata": {},
   "source": [
    "## 1. Load exon counts and sample metadata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "66bf6f94",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-04-20T19:13:14.512700Z",
     "iopub.status.busy": "2026-04-20T19:13:14.512264Z",
     "iopub.status.idle": "2026-04-20T19:13:14.554108Z",
     "shell.execute_reply": "2026-04-20T19:13:14.553216Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "exon counts: (14599, 7)  (exons x samples)\n",
      "            condition         type\n",
      "untreated1  untreated  single-read\n",
      "untreated2  untreated  single-read\n",
      "untreated3  untreated   paired-end\n",
      "untreated4  untreated   paired-end\n",
      "treated1      treated  single-read\n",
      "treated2      treated   paired-end\n",
      "treated3      treated   paired-end\n"
     ]
    }
   ],
   "source": [
    "data = gd.load_pasilla()\n",
    "counts  = data['counts']\n",
    "targets = data['targets']\n",
    "print(f\"exon counts: {counts.shape}  (exons x samples)\")\n",
    "print(targets[['condition', 'type']])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0411582",
   "metadata": {},
   "source": [
    "## 2. Library sizes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "aa7bb10e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-04-20T19:13:14.559692Z",
     "iopub.status.busy": "2026-04-20T19:13:14.558995Z",
     "iopub.status.idle": "2026-04-20T19:13:14.683921Z",
     "shell.execute_reply": "2026-04-20T19:13:14.682108Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk4AAAEiCAYAAAAPh11JAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAM4tJREFUeJzt3Qd0VFX3+P2NdKSGHnp/6F0w0pQOghQfEES6ghQpUvWRIk16LyJFkKoCAcS/SEeqICAkIF0IJaBSApHOfdc+v3dmJZCECUzKzP1+1pqVmTuTmZOdm5mdc/Y5J4FlWZYAAADgmV569kMAAABA4gQAABAN9DgBAAC4iMQJAADARSROAAAALiJxAgAAcBGJEwAAgItInAAAAFxE4gQAAOAiEicghn399deSIEEC2b9/f6SP+fPPP81j9LEOQ4YMMcf+/vtvr/8d5c6dW9q2bRvXzYCbRHQ+A94iUVw3AIBI1qxZZffu3ZIvXz5bhmPVqlWSOnXquG4GADwTiRMQDyRNmlQqVqzo1ue8c+eOJE+e/IWe499//5UUKVJITCtdurTY1YMHD0zvTKJEvB0DnoChOiCeD20EBQVJkyZNTI9MmjRppFWrVvLXX389NdT15ptvysqVK00SkixZMhk6dKi5b/r06VKlShXJlCmTvPzyy1K8eHEZM2aM+cAOq1q1alKsWDHZvn27+Pn5mYSpffv20qFDB/Hx8TFJ1JPeeOMNKVq0aJQ/28GDB03b9PU1QfT19ZX69evLhQsXIh2q07ZoPCK6hI1RcHCwdOrUSbJnzy5JkiSRPHnymJ/74cOH4dowc+ZMKVmypKRMmVJSpUol//nPf+STTz5x6XeisRoxYoTkzJnTxLVcuXKyadOmpx5/8uRJadmypfPnLFy4sIl9WFu3bjXP+c0338jHH38s2bJlM489depUpO14Vtv1XOjSpYsUKVLEPEZfX38vv/zyS4Q/z9ixY2X06NEm5ppYa6xPnDhhzocBAwaY34+eZ40bN5arV69GeJ5pD2GJEiVMPPLmzStTpkyJMpbRiREQ3/EvDhDP6QdYs2bNpHPnzhIYGCifffaZHD16VPbu3SuJEyd2Pu7AgQNy7Ngx+d///mcSCE2S1OnTp82HlR7T5OL33383icAff/wh8+bNC/daly9fNolZv379ZOTIkfLSSy9J2rRpzeOWLFkiHTt2dD5W27Bly5YoP/hCQ0OlZs2a5rX1cZkzZzbJjn7frVu3Iv2+GTNmSEhISLhj+nPr9xUqVMjc1ud55ZVXTBsHDRpkhjl1uHP48OEmSZg/f7553LJly0xi0b17dxk3bpx5vCYq2n5XTJs2TXLlyiWTJk2Sx48fm0Sqbt26sm3bNnn11VedsdBkU5Or8ePHS5YsWWT9+vXy0UcfmRq1wYMHh3vOgQMHmu+dNWuWaY8mEhFxpe3Xrl0zX/U19HVv375tEhtNiDTB069h6e9Bkx79euPGDZPANWjQQCpUqGDOJ/1dnzt3Tvr06WN+32vWrAn3/YcOHZKePXuaGjx9vcWLF0uPHj3k/v375nsiE90YAfGWBSBGzZ8/39I/tX379kX6mLNnz5rH6GMdBg8ebI716tUr3GMXL15sji9atMh5LFeuXFbChAmt48ePR9mWR48eWQ8ePLAWLlxoHn/t2jXnfVWrVjXPu2nTpqe+T+8rVapUuGMffvihlTp1auvWrVuRvt7+/fvNc/r7+0fZLm1/mzZtIr1/7Nix5nlmz57tPNapUycrZcqU1rlz58I9dty4ceaxgYGB5na3bt2stGnTWtHl+J34+vpad+7ccR4PCQmxfHx8rBo1ajiP1a5d28qePbt18+bNcM+hr50sWTJnnLds2WKes0qVKi614Xna/vDhQ/M7rl69utW4ceOnfp6SJUua88Bh0qRJ5njDhg3DPU/Pnj3N8bA/k/6eEiRIYB06dCjcY2vWrGnOhdDQ0EjPZ1djBMR3DNUB8dy7774b7rb2Pmk9jPa+hKW9CAULFoxwqKxhw4aSPn16SZgwoelVaN26tTx69MgM0YSVLl06M8zzJO1R0J6GnTt3mtvaG6TDTW3atDHDQ5HJnz+/ec7+/fub3hVXe3nCWrp0qekB0560999/33n8hx9+kNdff90MLenQnOOivUFKe4SU9kppz0qLFi1k9erV0Z6lqMOkOiTloMNl2kOjQ5oaw7t375qeHe0Z1OHNsG2pV6+euX/Pnj3hnrNp06YuvbarbdfYlilTxrRTzw39HWubtAfySdom7bly0OEypcOnYTmOnz9/PtxxHZrVocOwtEdTzwnt9YzI88QIiK9InIB4Toc0wtIPRk2C/vnnn6dm5j1JP/QqV64sFy9elMmTJ5u6l3379jmH17SA/FnPod566y1T3+L4Pq0z0mG4rl27Rtl2rZXRBKZUqVKmLkc/dDXR0WGZJ2usIqLJodY+aaI3bNiwcPdduXJF1q5da5KEsBdHzZUjyXjvvfecw0+asOiwmA5LbdiwQZ4n/o5jOjSlw2L6e9AEYOrUqU+1RZOCsG15Vpyf5ErbJ0yYIB9++KE5vmLFCpOA6O+4Tp06T/1+ldarhaXDt1Ed16TGlXioJ89Jh+eJERBfUeMExHNay6NFxA76AaQfRJo8haWFv0/y9/c3CY4WjWudjoP2HkUkoudQ2kOhSZImP1qfojVI1atXd9YbRUWL0bVWx7IsOXz4sEm6Pv/8c1OYrMXIkdHHNmrUSKpWrSpfffXVU/dnyJDB9LJpvVZENEFzaNeunbloLLSnSBM3LXLWHrewcYks/hEd08RCe9v0w1978jTJiSyR1BovV+IckWe1fdGiRaaOSYvIw4qqhuxFRBYP9eQ56aC9jtGNERBfkTgB8ZwW35YtW9Z5+9tvvzXJ05NFvxFxfEDrDCYHTWAiSkSeRQuFtSBYhw6PHz9uZmZFh7ZFh3gmTpxokqfIhnUcPWU65KYztrQXJWwRvIMmDz/++KMpCtcPZldowbw+r/YWaVKmxfbPSpw06dSZaI7hOk1ItKdLe/I0GdChJx0y1CFRTeQcPTXuFlnbNa5hf7+OpFML5XPkyOH2dujr6gSDsMN1OnFAhzB1uDAisRUjIDaQOAGxZPPmzWa215McQxVRfXDr8JzOTnPMqtMPLa11ehb9Hv2Q0hoZrRPSYRftmbh+/Xq026+z63TITL9fP7C1zudZtA5Je6f0g16TIE3a9OfRuh1tW2Q0QdDH6Iw2/ZnD0kQpY8aMptdKh6x0ppbOzNLeL/35NMaaUGndjy5ToHVR2rv12muvmSEy7R0ZNWqUGUYsX778M38GTY60rb179zaz6jRh1Hoex3IPSodBK1WqZJIpHTbTYU1NsHQGnCZZ+rt/Hq60XRNIHcbUnijtndOkVmOjPThPLsvgDtqTpzVzmkRrm7THS38PGpeo1vyKqRgBsY3ECYglWiAdkbNnz0b5fZpo6IeUJizau6AJi06Nd+W/dl3zR3tstLBai5x1KEULeTUJcBRRR0fz5s1NO/SDL2yBcWQKFChgEi6dwn/p0iXTZk1wtMdJC8sj4ygi1zY/SZcZ0Lon/dDWbWw0adAeIV0XSns9NGHQ+h5HL5R+UOvraU+dJow6xKcf4AsXLjQJ2LN069bNJGSanOm6RlpDtW7dOpPMOOgaStqDpm3RWOvj9OfWn/9ZiXFUXGn7p59+atbYmjt3romztkWTRl2SQNeNcjetV9OhQ03UdF0mTaS0zqpXr15Rfl9MxQiIbQl0al2svyoAj6Rr/mjipItyRlbP4i2050qTME3KolqfyE60l0gXSdWeRMCu6HEC8Ew6U0uLkXXYTVfq9vakCQAiQ+IE4Jl0lWutX9F6Gl2ZGwDsiqE6AAAAF7EAJgAAgItInAAAAFxE4gQAAOAiry8O1wXrdP0YXd8lOtscAAAAe7AsyyzIquuSPWuNOq9PnDRpioltBwAAgHfRNep0xwFbJ07a0+QIRurUqeO6OQAAIJ7RbZS0k8WRM9g6cXIMz2nSROIEAAAi40pJD8XhAAAALiJxAgAAcBGJEwAAgItInAAAAFxE4gQAAOAiEicAAAAXef1yBIg7tYet8+rwr/+sflw3AQAQy+hxAgAAcBGJEwAAgItInAAAAFxE4gQAAOAiEicAAAAXkTgBAAC4iMQJAADARSROAAAALiJxAgAAcBGJEwAAgItInAAAAGIycXrw4IEEBQXJ8ePH5dq1a/K8Ro0aJeXLl5dUqVJJpkyZpFGjRuY5w7IsS4YMGSK+vr6SPHlyqVatmgQGBj73awIAAMR44nT79m358ssvTeKSJk0ayZ07txQpUkQyZswouXLlkvfff1/27dsXrRfftm2bdO3aVfbs2SMbNmyQhw8fSq1atSQ0NNT5mDFjxsiECRNk2rRp5vmzZMkiNWvWlFu3bkXvJwUAAHhBCSzt0nmGiRMnyogRI0yy1LBhQ3nllVckW7ZspgdIe5wCAgLkl19+kVWrVknFihVl6tSpUqBAgWg35q+//jI9T5pQValSxfQ2aU9Tz549pX///uYx9+7dk8yZM8vo0aOlU6dOz3zOkJAQk+jdvHlTUqdOHe024fnVHrbOq8O3/rP6cd0EAIAbRCdXSOTKE+7atUu2bNkixYsXj/B+TaTat28vs2bNkrlz55rE53kSJ22w8vHxMV/Pnj0rwcHBphfKIWnSpFK1alXTpogSJ02s9BI2GAAAAO7gUuL03XffufRkmtR06dLluRqivUu9e/eWSpUqSbFixcwxTZqU9jCFpbfPnTsXad3U0KFDn6sNAAAAMTqrTnt0/P395dixYy/0PN26dZPDhw/L0qVLn7ovQYIETyVZTx5zGDhwoOm5cly0iB0AACBOEqdmzZqZQm11584dKVeunDlWokQJWbFixXM1onv37rJmzRozHJg9e3bncS0ED9vz5HD16tWneqHC9nrp+GTYCwAAQJwkTtu3b5fKlSub61oMrr0/N27ckClTpsjw4cOj9Vz6vdrTtHLlStm8ebPkyZMn3P16W5MnnXHncP/+fVND5efnF92mAwAAxG7ipMNfjuLtn376SZo2bSopUqSQ+vXry8mTJ6P1XLoUwaJFi2TJkiVmLSftWdKL9mQpHY7TGXUjR440SZrO3mvbtq15vZYtW0a36QAAADFfHB5Wjhw5ZPfu3SZ50sRp2bJl5vj169clWbJk0XqumTNnmq+6NlRY8+fPNwmS6tevn0mktOhcX6NChQry888/m0QLAAAgXidO2gP07rvvSsqUKc3Cl46kR4fwIluuIDIuLCFlep105XC9AAAAeFTipD0/um6TzlbTFbxfeun/Rvvy5s0b7RonAAAAr06clM6k00tYWuMEAAAgdk+cdGFKV+m+cgAAALZNnA4ePBju9m+//SaPHj2SQoUKmdsnTpyQhAkTStmyZWOmlQAAAJ6SOOnClGF7lHRG24IFCyRdunTmmM52a9eunXN9JwDAi2OjbMAL1nEaP3682Q/OkTQpva6F4XofAACAt3rpefamu3LlylPHdRuUW7duuatdAAAAnp84NW7c2AzLff/993LhwgVz0esdOnSQJk2axEwrAQAAPHE5glmzZkmfPn2kVatW8uDBg/97kkSJTOI0duzYmGgjAACAZyZOuk/cjBkzTJJ0+vRps/p3/vz55eWXX46ZFgIAAHjyAphKE6USJUq4tzUAAADeljjt27dPvvvuOzl//rzcv38/3H0rV650V9sAAAA8uzh82bJl8tprr8nRo0dl1apVps5Jr2/evFnSpEkTM60EAADwxMRp5MiRMnHiRPnhhx8kSZIkMnnyZDl27Jg0a9ZMcubMGTOtBAAA8MTESQvCHRv6Jk2aVEJDQyVBggTSq1cvmT17dky0EQAAwDMTJx8fH+dCl9myZZOAgABz/caNG/Lvv/+6v4UAAACeWhyu+9Ft2LBBihcvbobnevToYeqb9Fj16tVjppUAAACemDhNmzZN7t69a64PHDhQEidOLDt27DCrhn/22Wcx0UYAAADPTJx0qM7hpZdekn79+pkLAACAt4t2jZOjQPx///uftGjRwmzuq3766ScJDAx0d/sAAAA8N3Hatm2bqW/au3evWezy9u3b5vjhw4dl8ODBMdFGAAAAz0ycBgwYIMOHDzfF4LqOk8Prr78uu3fvdnf7AAAAPDdxOnLkiDRu3Pip4xkzZpR//vnHXe0CAADw/MQpbdq0cvny5aeOHzx40KzrBAAA4K2inTi1bNlS+vfvL8HBwWbF8MePH8vOnTulT58+0rp165hpJQAAgCcmTiNGjDB70mnvkhaGFylSRKpUqSJ+fn5mph0AAIC3itY6TpZlyaVLl+Srr76SYcOGyYEDB0yPU+nSpaVAgQIx10oAAABPTJw0QdL1mvRr3rx5Y65lAAAAnjxUpyuFa8LE7DkAAGBH0a5xGjNmjPTt21cCAgJipkUAAADekji1atVKfv31VylZsqQkT57c7F0X9hId27dvlwYNGoivr6+Zoefv7x/u/rZt25rjYS8VK1aMbpMBAADiZpPfSZMmueeVRSQ0NNQkYO3atZOmTZtG+Jg6derI/PnznbfDrlYOAAAQrxOnNm3auO3F69atay5RSZo0qWTJksVtrwkAABBrQ3WxbevWrZIpUyYpWLCgvP/++3L16tUoH3/v3j0JCQkJdwEAAPD6xEl7oxYvXiybN2+W8ePHy759++SNN94wyVFkRo0aJWnSpHFecuTIEattBgAA3ivaQ3WxqXnz5s7rxYoVk3LlykmuXLlk3bp10qRJkwi/Z+DAgdK7d2/nbe1xInkCAABenzg9KWvWrCZxOnnyZJQ1UXoBAACw1VDdk3ThzaCgIJNAAQAAxNsep/bt27v0uHnz5rn84rpJ8KlTp5y3z549K4cOHXKuCTVkyBCzTIEmSn/++ad88sknkiFDBmncuLHLrwEAABDridPXX39thsl0Q1/ds84d9u/fL6+//rrztqM2SZc8mDlzphw5ckQWLlwoN27cMMmTPnb58uWSKlUqt7w+AABAjCROnTt3lmXLlsmZM2dM75OuIB7dlcKfVK1atSiTsPXr17/Q8wMAAMRJjdOMGTPk8uXL0r9/f1m7dq2ZqdasWTOT3LirBwoAAMBrisN1tlqLFi1kw4YNcvToUSlatKh06dLFDOFpvRIAAIA3e+7lCByb7mpv0+PHj8XOag9bJ95s/Wf147oJAAB4Xo+Trti9dOlSqVmzphQqVMgUb0+bNk3Onz8vKVOmjLlWAgAAeFKPkw7JaXF4zpw5pV27duZ6+vTpY7Z1AAAAnpg4zZo1yyRNefLkkW3btplLRFauXOnO9gEAAHhe4tS6dWtT0wQAAGBX0VoAEwAAwM48aq86AACAeJ846arhurmuK3RLlMWLF79ouwAAADxzqC5jxoxSrFgx8fPzk4YNG0q5cuXE19dXkiVLJtevXzeLYe7YscPMtMuWLZvMnj075lsOAAAQHxOnYcOGSffu3WXu3Llmdl1AQEC4+3XT3Ro1asicOXOkVq1aMdVWAAAAzygOz5QpkwwcONBcbty4IefOnZM7d+5IhgwZJF++fMy4AwAAXu+5tlxJmzatuQAAANgJs+oAAABcROIEAADgIhInAAAAF5E4AQAAxGTi9PDhQ9m4caN8+eWXcuvWLXPs0qVLcvv27ed5OgAAAO+cVafLENSpU0fOnz8v9+7dk5o1a5p1nMaMGSN379416zwBAAB4o2j3OPXo0cOsHK4rhidPntx5vHHjxrJp0yZ3tw8AAMBze5x0a5WdO3dKkiRJwh3PlSuXXLx40Z1tAwAA8Owep8ePH8ujR4+eOn7hwgUzZAcAAOCtop04aU3TpEmTnLcTJEhgisIHDx4s9erVc3f7AAAAPHeobuLEifL6669LkSJFTDF4y5Yt5eTJk2bPuqVLl8ZMKwEAADwxcfL19ZVDhw6ZJOnAgQNm6K5Dhw7y7rvvhisWBwAA8DbPtcmvJkjt27c3FwAAALuIdo1TwoQJzVDdtWvXwh2/cuWKuQ8AAMBbRTtxsizLLHypazkFBAQ8dR8AAIC3inbipLPoVqxYIQ0aNBA/Pz9ZvXp1uPsAAAC81XP1OOmQ3OTJk2XcuHHSvHlzGT58+HP1Nm3fvt0kYFpwrkmXv7//U681ZMgQc7/WVVWrVk0CAwOj/ToAAABxVhzu8MEHH0jBggXl7bfflm3btkX7+0NDQ6VkyZLSrl07adq06VP36/53EyZMkK+//tq8jiZouo7U8ePHWWwTAIAwag9b59XxWP9ZffHIxEm3VglbBK69QHv27DE9R9FVt25dc4mI9jbpQpuffvqpNGnSxBxbsGCBZM6cWZYsWSKdOnWK9usBAADE6lDd2bNnJX369OGO5c+fXw4ePChnzpx5ocY8+TrBwcFSq1Yt57GkSZNK1apVZdeuXZF+nxauh4SEhLsAAADESeIUmWTJkpneKHfRpElpD1NYettxX0RGjRoladKkcV5y5MjhtjYBAAB7cylx8vHxkb///ttcT5cunbkd2cXdnpypp0N4Uc3eGzhwoNy8edN5CQoKcnubAACAPSVydX+6VKlSOa/HxrIDWbJkMV+1dylr1qzO41evXn2qFyosHc7TCwDP4s2FrfGlqBVALCVObdq0cV5v27atxIY8efKY5GnDhg1SunRpc+z+/ftm9t7o0aNjpQ0AAAAvNKtON/ZNnDixFC9e3NzWBTDnz58vRYoUMWsuJUmSxOXnun37tpw6dSpcQbhuIKxDfjlz5pSePXvKyJEjpUCBAuai11OkSCEtW7aMbrMBAF7Cm3snFT2UXlYcrssAnDhxwlzXWXS6AKYmM999953069cvWs+1f/9+05vk6FHq3bu3uT5o0CBzW59Pk6cuXbqYLV4uXrwoP//8M2s4AQAAz+hx0qSpVKlS5romS7o8gK6rtHPnTnnnnXfM2kuu0jWgolpxXGuptBdLLwAAAB655crjx4/N9Y0bN0q9evXMdZ3275h5BwAA4I2inTjpkJluffLNN9+YQu369es765Oimu0GAABgu8RJh+K0QLxbt25mOxRdNVx9//334ufnFxNtBAAA8MwapxIlSsiRI0eeOj527Nhwe9gBAACI3ROnqLZcAQAA8GZu26sOAADA25E4AQAAxPZQHQDXsOoxANiox2nr1q0x0xIAAABvS5zq1Kkj+fLlM2s5BQUFxUyrAAAAvCFxunTpkvTo0UNWrlwpefLkkdq1a8u3334r9+/fj5kWAgAAeGri5OPjIx999JFZBFM36S1UqJB07dpVsmbNao7//vvvMdNSAAAAT55Vp5v9DhgwwCROoaGhMm/ePClbtqxUrlxZAgMD3ddKAAAAT02cHjx4YLZY0Q1+c+XKJevXr5dp06bJlStXzJ51uuHvf//7X/e3FgAAwJOWI+jevbssXbrUXG/VqpWMGTNGihUr5rz/5Zdfli+++EJy587t3pYCAAB4WuJ09OhRmTp1qjRt2lSSJEkS4WN8fX1ly5Yt7mgfAACAZw7V6RBdzpw5pUKFCpEmTSpRokRStWpVd7QPAADAMxOnxIkTy6pVq2KuNQAAAN5UHN64cWPx9/ePmdYAAAB4U41T/vz5ZdiwYbJr1y6z9IAWg4elazkBAAB4o2gnTnPmzJG0adPKb7/9Zi5hJUiQgMQJAAB4rWgnTrpOEwAAgB290MrhAAAAdhLtHid14cIFWbNmjZw/f/6pzX0nTJjgrrYBAAB4duK0adMmadiwoeTJk0eOHz9uVg3/888/xbIsKVOmTMy0EgAAwBOH6gYOHCgff/yxBAQESLJkyWTFihUSFBRkFrxkfzoAAODNop04HTt2TNq0aeNcIfzOnTuSMmVK+fzzz2X06NEx0UYAAADPTJx03aZ79+4596Q7ffq0876///7bva0DAADw5BqnihUrys6dO6VIkSJSv359M2x35MgRWblypbkPAADAW0U7cdJZc7dv3zbXhwwZYq4vX77crCg+ceLEmGgjAACA5w3VPXr0yBSC58iRw9xOkSKFzJgxQw4fPmx6nHLlyuXWxmlipquRh71kyZLFra8BAAAQIz1OCRMmlNq1a5sC8XTp0klsKFq0qGzcuDFcGwAAADxiqK548eJy5swZs45TbNCZe/QyAQAAj5xVN2LECOnTp4/88MMPcvnyZQkJCQl3cbeTJ0+a2XuaqL3zzjsmaYuKzviL6TYBAAB7inaPU506dcxXXT1ca44cdOVwva11UO5SoUIFWbhwoRQsWFCuXLkiw4cPFz8/PwkMDJT06dNH+D2jRo2SoUOHuq0NAAAAz504bdmyRWJL3bp1ww0Rvvrqq5IvXz5ZsGCB9O7dO9KVzcPepz1OjmJ2AACAWE2cdGuVuKKLb2oCpcN3kUmaNKm5AAAAxHnipK5fvy5z5841s+t0eK5w4cLSrl078fHxkZik9Uv6mpUrV47R1wEAAHBLcfi2bdskd+7cMmXKFJNAXbt2zVzX4m29z520CF2f8+zZs7J37155++23zdCbY688AACAeN3j1LVrV2nevLnMnDnTuaaSFoR36dLF3BcQEOC2xl24cEFatGhh9sDLmDGj2dJlz549bl9oEwAAIEYSJ93Ud8WKFeEWotTrWpCtM+DcadmyZW59PgAAgFgdqitTpoypM3qSHitVqtQLNQYAAMDje5x0LzqHjz76SHr06CGnTp0yQ2dKh8+mT58uX3zxRcy1FAAAwBMSJ+1J0tlzusilQ79+/Z56XMuWLU39EwAAgG0TJ53VBgAAYHcuJU7MYgMAAHAxcVqzZo3Z/iRx4sTmelR0DzsAAADbJk6NGjWS4OBgyZQpk7keGXdv8gsAAOBxidPjx48jvA4AAGAn0V7HKTJBQUHSvn17dz0dAACA9yZOumfdggUL3PV0AAAA3ps4AQAAeDsSJwAAABeROAEAALhzVp1q0qRJlPffuHHD1acCAADw7sQpTZo0z7y/devW7mgTAACAZydO8+fPj9mWAAAAxHPUOAEAALiIxAkAAMBFJE4AAAAuInECAABwEYkTAACAi0icAAAAXETiBAAA4CISJwAAABeROAEAALiIxAkAAMBFJE4AAAAuInECAABwEYkTAACAi0icAAAAvClxmjFjhuTJk0eSJUsmZcuWlV9++SWumwQAAGwo3idOy5cvl549e8qnn34qBw8elMqVK0vdunXl/Pnzcd00AABgM/E+cZowYYJ06NBBOnbsKIULF5ZJkyZJjhw5ZObMmXHdNAAAYDOJJB67f/++/PbbbzJgwIBwx2vVqiW7du2K8Hvu3btnLg43b940X0NCQmKsnQ/v/ive7HljR1yIC+cLf0O8t/Ce6y4x+TnueG7Lsp79YCseu3jxov4E1s6dO8MdHzFihFWwYMEIv2fw4MHme7gQA84BzgHOAc4BzgHOAYlGDIKCgp6Zm8TrHieHBAkShLutGeGTxxwGDhwovXv3dt5+/PixXLt2TdKnTx/p93gSzYp1qDIoKEhSp04d182JN4gLceF84W+I9xbec5+X5hW3bt0SX1/fZz42XidOGTJkkIQJE0pwcHC441evXpXMmTNH+D1JkyY1l7DSpk0r3kaTJhIn4sL5wt8R7y2858al1F70WZQmTRrPLw5PkiSJWX5gw4YN4Y7rbT8/vzhrFwAAsKd43eOkdNjtvffek3Llysmrr74qs2fPNksRdO7cOa6bBgAAbCbeJ07NmzeXf/75Rz7//HO5fPmyFCtWTH788UfJlSuX2JEOQw4ePPip4Ui7Iy7EhfOFvyHeW3jPjQ0JtEI8Vl4JAADAw8XrGicAAID4hMQJAADARSROAAAALiJxAgAAcBGJEwDAlpgbRVyeB4lTPKTbxACcJ8/Hscn3o0ePOImeQEzC84ZtuBD7WI4gHjl16pTkz5//mfvx2cmWLVtk9erVki5dOnnllVekbt26cd2keEFXz1+yZImkSpVK3njjDWnUqFFcNyle2LRpk3Ts2FG2bt1q1nrTf0Jeeon/D0+fPi358uVzJk+6lZXdbd++3by/ZM2aVUqVKmXeXyCyefNmWbVqleTOnVsqVKgglSpVIixP4B0lnlizZo0ULFjQrJKuNGmye8/TunXrpEePHuYDUN/sJ06cKMePHxe7W79+vQwYMMC8oWnipOcO/s/Bgwfl3LlzUrVqVZMsaNJk9+EYPT9066oPPvjA3Nakye49T//v//0/6dSpkyRLlsz8Pek/IaGhoWJ3P/30k/Ts2VMKFSokGzduNElUWHb/W3LSBTARt06fPm2VK1fOmjdvnpU/f37rvffec9736NEjy44CAwOtYsWKWTt27DC3z507Z9WoUcPatWuXZWcHDx608uTJY23dutXcXrlypdWgQQPr+++/t9atW2fZ1ePHj83XAwcOWN988431+eefW2nTprXu3LljXbt2zQoJCbHs6NSpU1bhwoWtYcOGWf/973+tzp07O+97+PChZUeHDx+2ihYtam3evNn5N5U9e3Zr//79lp1pHAoUKGBt2rTJ3F60aJHVpk0ba+PGjbZ/330SiVM8oR986vr161a2bNmsVq1aWXamH3aLFy82b+6OD8V33nnHmjFjhmVnGotff/3VXL9w4YKVL18+q3Xr1tbo0aMtHx8fa+7cuZadY6MfflWrVjW3BwwYYKVMmdLKkCGDFRQU5DyP7OTevXvWjz/+aBLHvXv3WnXr1g2XPNnRyZMnrWXLloVLHt966y1rw4YNlp3p+8m+ffvM9fPnz1vp06c3iVPHjh2tN9980xkzWBZDdXFo27ZtMmjQIHO9adOm8uDBA0mbNq0EBgaasfdWrVqZ+7799ltTu2GXmAwdOtTUNDVp0sQMKzi6h3U/OkdtxsqVK2XHjh1it3NFh3DLly9vjulm1/369ZMFCxaYr3p/QECA2InjfFEamxIlSkiOHDnM7ffff1+SJEni/LuyU83goUOHzNC2/vxaF6hDunreaKx0KPPDDz80j9NasF9++UXsEpMpU6aYOlJHraSj/k3fY/79919zXYeo9u7dK3ahcZk8ebJky5ZNypUrZ4Zxz549KyNGjJCvv/5axowZY+rjLl26FNdNjTdInOKQvqlNmzZNRo0aZW4nTpxY7t69K2nSpDHJ0549e8ymxv379zcFjHaJif4Rjx492tQfhJ0lpbfz5s1rahJGjhwpmTJlErudK1988YXz2KuvvuqsW1EhISG2K4R2nC+OuOjfkP79dO3aVWrVqmUS7DZt2piE6v79+2IHmigeOXLE/J1MnTrVeVwTRy2C1uRJN04vU6aMKaTXD0y7xEQ3iJ8+fbqkTp3aHNf3W6VJU5YsWcxElE8//dQ277eOuGjNl+Nc0X9OtX5Sa8CU/hObMmVK8/6iLOqcqHGKC2GHDLT7U2sxJk6c6DymdRlq0qRJZpghICDAsntMHEMv5cuXtypWrGgdOXLEsoNnxcVRA6fDmmXKlLGOHTtm2TUu48aNM7enT59uagX9/f2dj7l8+bJlB2FrIqdOnWrVqlXLmjNnzlOP69u3rykJsMPfkSsx+eCDD6y3337bqlSpkqmBsoNnxcVxv9YMlihRwjp+/HictDM+osYpDukHoNYbaHGv1qdoQauDFkGXKlXK+v333y07iSomWjSfIkUKU/BqN0/GZejQoeb4v//+a5KmQoUK2SLBjiwuDRs2tNKlS+dMnnRygXrw4IH5arf6ppkzZ1pNmzY1/2i89tpr4ZJtrXVKlSqV7d5booqJ1vDoP6l2TA4ii4v+A79gwQKrSJEitnxviQqJUxzZs2ePmdkRGhpqLrt37zb/AYb9Yw4ODrbsJLKYOD4MDx06ZJ04ccKym8jioj2S6o8//rDOnj1r2T0u+s9G1qxZzX/PDnZLmNT69etND4H2GFy8eNH0GGjCHXbiwKVLlyw7iSwms2fPdk7O0aJxu4ksLvPnzzf36+zdM2fOxHUz451EcTrAamM3b94048YpUqQwtytWrCjNmjWT3r17m+K8jz/+WDJnzix2EllM+vbta2pZunfvLnYUWVx69epl6lY++ugjsaMn46I1X82bNzfx0L8hXQPMTgXhDlqvowXyWu/m6+srtWvXlu+//17Gjh0rd+7cMfVfWs9jJ5HFZNy4cea9Revg7CiyuGjdrXastGvXLq6bGC/Zq5I0joQtpnMsaqkffhkyZJAvv/zSeV+RIkVM0XODBg3E20U3JvoHbQfRjUudOnXEDlyNS9GiRU1c7LLCfESFujoDShNGnVyiEysyZsxoin1btGjhXGHemxPK6MTk3XfflRo1aogdRDcuOrkCEWPLlRgWdusUfYO/ePGimcnRuXNnM+tFT9hbt25JlSpVZPbs2WYrDW+f5UJMiAvni3v/jubOnWtmPV29etX0FmgvbXBwsGTPnt18GH7zzTfObTS8GTEhLrGBoboY5nhj06nk2gWqyZP2Fug2Iq1btzb/NS9evNi84el6Td6eNCliQlw4X9z3dzRp0iTx9/c3CVPlypXNkIsuz6Dre+kWRTrdfNGiRV6fNCliQlxiRVwXWXk7LU7VorsmTZqYVXx1qnTt2rWt+/fvh3ucnbZWISbEhfPFPXSyhK56rXQShc4uvHv3brjH6PuOnRAT4hLTqHGKAWEX2tPiO11sTbuQP/nkEzM8p13mulCfLkQ3b9488zhvX7iQmBAXzhf3LFgY9m9Ki+P1ohMFdHcB7bXWFfbHjx9vrit9r/FmxIS4xDbv/rSOo5k+mhhpLdOyZcvMCtj6xpY8eXJZvny5qTXQ60uXLjWrGvv5+Ym3IybEhfPlxV27dk127txpapmWLFli3ks0KXr48KHs3r3b/BOmSZMO0emlbNmyXl8ITkyIS1ygxsnN9L9A3c6gZs2a5k3NsZ9at27dnEvZ6ywO3SNK3/z+85//iLcjJsSF88U9/4CsWbPGzBrUf8x0jzkfHx8z+2n79u1m1pzuNaa92vqPmc6Y8nbEhLjEBRInN9MZLC+//LIEBQWZQs2//vrLucGm7g2lvUz6ZqdrEuXJk0fsgJgQF86XF6fvF7rmjs6g69Chg1mrSrVv315eeeUVOXjwoFmTSGfs2iFpUsSEuMQFliNw8xRYxzozv//+u9k48bfffpMBAwaYxGnXrl2mh0kTJ29HTIgL54v7/470HzJNkLQcIGfOnGbBT52lGxAQYG47Nq/1ZsSEuMQ1epzc+EesNU3Xr18304F1Ab706dObeoRhw4aZhfm0PkG70L0dMSEunC/u/Ttat26d3Lhxw9QtNWzY0Ly36G72usSJTkDRf9QWLlwo3o6YEJd4Icbn7dnElClTzAaJI0aMsMqWLWv16tXLCgkJMbuy695Z9erVc248ahfEhLhwvry4yZMnWxUrVrR69OhhFS5c2Jo1a5b18OFDszffZ599ZtWtW9c6fPiwZSfEhLjEJRInN1i2bJlVqVIlc/2LL74wCVTLli2t3r17W9evX7flWirEhLhwvry4b7/91ryfPHjwwGwAromT/hOmyZPuXq+eXLfJ2xET4hLXSJyew5M7ruvu0adOnbK++uorq0aNGuaNbNiwYVaOHDnMf4m6uKW379JOTIgL54t7/470+v79+80CujNnzrRq1apljnfv3t3y9fW15syZY3qevB0xIS7xDTVOL1B3oDUHyjE7Thed0122dS0VrXNq1aqV2and2xe3JCbEhfPFvX9HuuGq3taaJl2nSSeZdOnSxdyns3O1lvLNN980S5x4M2JCXOIjEqdocryxTZgwQX799Vc5c+aMdOzY0SRJOlVYCzR11ouu06TF4lmyZBFvR0yIC+eL+/6OtOh77969Zv9KXZupXbt2ZqaurtW0f/9+M8lElyTInDmzeDtiQlziI+/uCokh06dPN0sN6Oa8uiq4LnKZIkUKqV+/vjRq1Mi84emqvgUKFBC7ICbEhfPlxelm4Lphr27Sq3TdN6U92YkSJZJjx47JlClTbLMGnCImxCXeieuxQk8bY9d6pU8++cS6ffu22VRTN+x1FH5funTJ+RhvR0yIC+eLe/+O9H2kc+fOpkZSZ429+eab5pjevnr1qnmM3WqaiAlxiY/ocXp2YunsLtbtDv744w/zX1/16tXN/lC6voqu1jt27Fjp06eP2V7Em/eGUsSEuHC+uIfjvUIXtNRaSL2tWzL99NNPsnbtWvPeMmfOHLPnpb63eHu9pCImxCXei+vMzVOsXr3aTAvW//5+/vlnq2jRota8efPMffPnz7dKly5tBQQEWHZCTIgL58uL96p89913VqFChaygoCAze87Pz89au3atuW/BggVWiRIlrKNHj1rejpgQF09B4uSCHTt2mDczXXBO6fRgTZZy585tNWnSxCxOZ7fFLYkJceF8eXHbt2+3Wrdube3evdvcPnHihNW3b1/zfqMLW5YpU8Z2/5ARE+IS37FXXQTu3r1rusi1W/zcuXNmjzndNqVkyZLhtkzRDXx1tosWbeoWCN6MmBAXzpcXd+nSJTl79qy89tprMnPmTFmyZIn8/fffZhuVQYMGmQ3CdZsmvdy6dUsyZMhgNsn2ZsSEuHgaEqcn6BvWDz/8IKVKlTIz5zQ5GjhwoJk5p0sNFCxY0CRREW026a2ICXHhfHEPTRIqVqwoJUqUkNu3b8u8efNM7eSBAwfMrNy33npLkiVLJnZCTIiLp/H+SsNo0t3FL168KA0aNDBvam3btpU0adJItWrVzJoq+t/ixx9/bB5rh6RJERPiwvniHrowbqdOneTnn3+WChUqSN68eeWDDz4wywvoP2rLly83i1/aCTEhLp6GxOn/p71HDpo0+fj4mMUrdedx7XHRLnQ9rqv1hoaGmp4ob0dMiAvni/u99957snr1apk8ebIMHz7crAE3dOhQefDggVk8126JkyImxMWjxHWRVXybzeHYOFM31Rw0aJDZUHPjxo3mmO4bpUXRISEhlrcjJsSF8yVmaUF4kiRJrPHjx5uNa8uXL2/2vbQzYkJcPAE1TmHMmDFDtm3bJtmzZ5cxY8aYfaC0vikgIECyZcsm69evN9se6NYqdkFMiAvnS8zRbZt69OhhygH0PUdrn+yOmBCXeC+uM7f4wt/f36yXsm7dOqtkyZJWhw4drODgYHPfokWLrNGjR9tuyQFiQlw4X2LezZs3bdGLHR3EhLjEZ/Q4iZglBnRWy9tvv20KNnV6cOvWrU0vk04RtlMPkwMxIS6cLwDwNFsWh4ctelY6/DZ79mw5efKkWZdJ105ZtGiRBAYGmq1UHj58KN6OmBAXzhcAeDZb9zjp9N+6deua63379jX1TboEQbFixcyx69evm7VW7NTjREyIC+cLAETOtomT9iwVLVrUrKOiG/Wqrl27miE7LYguXbq02A0xIS6cLwAQNdskThGt8q27jevClpkyZRJ/f39zrE2bNhIUFGR2J9dtV7wZMSEunC8AED22SZwcdCiuevXqkitXLmcvS5kyZcxtXZROBQcHm8Uv7YKYEBfOFwBwje2Kw9OlS2fWSjl//ry5rRv56maba9eulc6dO5tjdkqaFDEhLpwvAOCaRGIzjRs3lqRJk0r58uVl7969kjt3bvnnn39kyJAh0rJlS7EjYkJcOF8AwDW2G6pz0IJw7WGqV6+ebNmyxWy6qUmUnRET4sL5AgBRs23ipI4ePWpm0WnvU6FCheK6OfECMSEunC8AEDlbJ04AAADRYbvicAAAgOdF4gQAAOAiEicAAAAXkTgBAAC4iMQJAADARSROAAAALiJxAgAAcBGJEwAAgItInAAAAFxE4gQAACCu+f8AVFNNjilSH3UAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 600x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "lib = counts.sum(axis=0)\n",
    "fig, ax = plt.subplots(figsize=(6, 3))\n",
    "ax.bar(range(len(lib)), lib.values / 1e6, color='steelblue')\n",
    "ax.set_xticks(range(len(lib)))\n",
    "ax.set_xticklabels(counts.columns, rotation=45, ha='right', fontsize=8)\n",
    "ax.set_ylabel('Library size (M reads)')\n",
    "ax.set_title('Library sizes per sample')\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7d7aad2",
   "metadata": {},
   "source": [
    "## 3. log-count density (KDE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "bd3ddb6d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-04-20T19:13:14.689510Z",
     "iopub.status.busy": "2026-04-20T19:13:14.689044Z",
     "iopub.status.idle": "2026-04-20T19:13:14.963151Z",
     "shell.execute_reply": "2026-04-20T19:13:14.962183Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArEAAAFUCAYAAAAzu2SBAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAlWRJREFUeJzs3Qd4VFXaB/D/9N4ySWbSOyEkhBJ6R5pgAcUVG3bXtp91d9XVLeq6rmVdK9iwK2KvqHSkl9CTkN7rpM1kev2ec9nEBAKEOinvz+ca5s6dO2funWTeOfc97+EFAoEACCGEEEII6UP4wW4AIYQQQgghp4qCWEIIIYQQ0udQEEsIIYQQQvocCmIJIYQQQkifQ0EsIYQQQgjpcyiIJYQQQgghfQ4FsYQQQgghpM+hIJYQQgghhPQ5FMQSQgghhJA+h4JYQoLkvffeA4/H61iEQiGio6Nx0003obq6ms7LccTHx+PGG288a8eH7YvtM5jq6+vx8MMPY+jQoVAqlZBKpUhJScG9996LwsLCju3+8Y9/dHnPiMViJCQkcNu1trZ2+97asGHDMc/HJmpMTk7m7p82bVqvOReni72Gzq/Dbrdzx6q71342fPDBBwgLC0NbWxv6i4KCAu79tGfPnmA3hZAeE/Z8U0LIufDuu+9i8ODBcDgc+PXXX/H0009j48aNOHjwIBQKBR30fm7nzp24+OKLucDyD3/4A8aPH88FE/n5+fjoo48wZswYtLS0dHnMzz//DI1GwwVRK1euxEsvvcTtZ+vWrVxg2k6lUmHZsmXHBKrs/VVcXMzd3x8sWbKky20WxD7++OPcv88kSO8O2/df/vIXPPTQQ/3m+DGDBg3Ctddei/vvv597fxDSF1AQS0iQZWRkYNSoUdy/p0+fDp/PhyeffBLffPMN96Fyph+4crn8LLWUnG0WiwXz58/nel5ZAMp64tux4Ov222/HF198cczjsrKyEBoayv171qxZaGpqwocffsjtY+LEiR3bLVq0CB9//DFee+01qNXqjvUssGXBMnv+/mDIkCHn7bnef/997njfeuut6G/Ylyj2t4i9jyZMmBDs5hByUpROQEgvM27cOO5neXk595P10LGepuHDh0Mmk0Gn0+GKK65ASUlJl8exoIcFxKw3l30AseD15ptvPu7zrFu3jnuMXq/n9hsbG4uFCxdygW871ps1duxYhISEcEHQyJEjuQCItenoy8qsN/GHH37AiBEjuP2lpaVxt9svb7PbrGeZ9Szu3r27y+PZJWl2GT0nJwczZszgtmOXa9mHauf2HA8Lxv74xz9yl9ZZL2ZUVBTuu+8+2Gw2nA6n04lHHnmky/7uvvvuLpfsGZfLhQcffBBGo5E73lOmTEF2dnaPL7O/9dZbqKurw7PPPtslgO2MnetTfc+0u/rqq7mfy5cv71hnNpvx5ZdfnvC9cTSPx4M///nPHa9z0qRJXM9vd9jrYcE3ez3t6Q7sfeT1eju2KSsr43qMn3/+ebzwwgvcNuz8s8B6+/btXfbH3udXXXUVIiMjIZFIYDAYuPfIvn37uk0nYPtm7x2GPW97WgU7H5s2beL+3fl4dE4RYPft2rXrhMdi6dKluOSSS6DVarus78nv6aeffso9x6uvvtrlsX//+98hEAiwevXqjnXfffcddzzY8WY9vuzLyrZt27o8rj29hP3esHPNeufZ8WHnlp3nzj7//HPud5ltw/aZmJh4zHuAfTliv6evv/76CY8BIb1GgBASFO+++y6LBAO7du3qsv6ll17i1r/55pvc7dtuuy0gEokCDz74YODnn38OfPLJJ4HBgwcHDAZDoK6uruNxU6dODYSEhARiYmICr7zySmD9+vWBjRs3dvvcpaWlAalUGpg1a1bgm2++CWzYsCHw8ccfBxYvXhxoaWnp2O7GG28MLFu2LLB69WpuefLJJwMymSzw+OOPd9lfXFxcIDo6OpCRkRFYvnx5YOXKlYGxY8dy7f7b3/4WmDhxYuCrr74KfP3114FBgwZxbbfb7R2Pv+GGGwJisTgQGxsbeOqppwKrVq0K/OMf/wgIhcLAxRdffMxzse3b2Wy2wPDhwwOhoaGBF154IbBmzRruGGo0msAFF1wQ8Pv9JzwPbF9sn+3Y9nPmzOGe+69//SvXlueffz6gUCgCI0aMCDidzo5tr7766gCfzw88/PDD3HYvvvgid/zZc3du4/HMnj07IBAIAlarNdATf//737n3hslk6rL+/vvv59azNhz93mLndMyYMR3bLl26lHstFoslkJ6ezr1vToa9Fh6PF/jTn/7EPQc7zlFRUQG1Wt3lddbW1nKvnx3PN954gzsX7D0jkUi491Ln9x9rX3x8fODCCy/k3oNsGTp0aECn0wVaW1s7tk1NTQ0kJycHPvzwQ+79/OWXX3K/C+z93Y69hvbXwc4P+z1h+7/lllsC27Zt45aioiLufnYO2fvxaKNHj+aWE6msrOT2u2TJkmPu6+nv6R133MG919t/79euXcu9hx577LGObdjvInse9v5gx2XFihWBrKws7nGbNm065v3AjhH7PWO/o+zcsON90003dWy3detW7vxdddVV3O/munXruPcIe28c7c477+R+l072e0NIb0BBLCFB0h5obN++PeDxeAJtbW2BH374IRAWFhZQqVTcBx/78GXb/Oc//znmw5QFk3/+85871rEPcbYt+1A8mS+++ILbdt++fT1ur8/n49r5xBNPBPR6fZcPORa0sPZUVVV1rGP7Zs8RERHBBZrt2IcyW//dd991rGOBEFvHgs/OWEDL1m/evLnLc3UOnJ5++mkuCDj6y0D7a2Qf2qcSxLYHQM8++2yX7Vgg0fnLRU5ODnf7oYce6rIdC+LZ+p4EsSzIMRqNgZ5qD1rYe4OdC/aF46OPPuKOPQseHQ7HMUEsC/bYvw8dOsTdxwK19oCyJ0FsXl4e93gWKHfWHmh1fp233357QKlUBsrLy7tsy74EsG3ZMescxLKg1ev1dmy3c+dObj07hkxjYyN3m305OJHOQSzDgnz2OHa8jtZ+bPbu3XvM877//vsnfJ729wD7ne3sVH5PWZDNAumEhIRAbm4uF+SytrcfB/Z7FhkZyR0b9u927O9DeHh4YMKECce8H45+r951113cl9T239H249/5y8HxvPXWW9y27LwT0ttROgEhQcYuBYtEIu6SIbskzy7Z/vTTT9xlQXY5nl0uvO6667jLse0L22bYsGHHjL5mlzAvuOCCjtt+v7/L41i+LcMuebJLvb///e+5HL+jUxM6pxzMnDmTuwTJLneydv7tb3/jcgIbGhq6bMv2yS67t2OXJRl2mbdzXm77+qMvfTNH5wBfc8013M/169cf9/ixY8TSKNjzd36tc+bMOe7o/BNhr5k5Oh3gd7/7HZfmsHbtWu52++CXK6+8sst27BIyqzTRWed2seXodIxTxc4/OxfsfLP3BkvzYIO9WG7t0aZOnYqkpCS888473GBBdrn8VFIJ2o/90eeGve6jXyc7Fyyvm1367/x6586dy91/9IChiy66iHtftcvMzOzy3mBpLKztzz33HJd2sHfvXu49fSbYZffw8HAuT7jdK6+8wqUgsBziE6mpqeF+sscf/bp7+nvKUiI+++wz7neInTf2XmDpDe3HgQ3oY8+zePFi8Pm/fUSzdAuW7sPSLY5Osbn00ku73GbHkaXEtP+Ojh49uuOcsec+UfWT9tdGFVJIX0BBLCFBxnLxWGDBPqDZh9eBAwc6Buew0kvsQ44FtCxo6bywD7PGxsYu+4qIiOhymwUrnR/DcgkZFhisWbOG+8BiuZ7sNlvYKPd2LOdx9uzZHbmbW7Zs4dr56KOPcutYNYXOWMDRGQuST7Sefch2xgIilp/bGQsCGPaBfzzsGLFjdvTxYV8K2LE7+hidDHsu1pb2vMp2LEhh7WlvS/tPdm5O9jqObhv74sCwPGSTyXTKubvs3LFzwfJC2evbvHnzcQc3sXazsm2s0gHLdWSj0CdPntzj52p/ne3n4kSvk52L77///pjXm56ezt1/9Lk4+vEswOv83mJtZ18a2BcSljfMgj52Xu65557TLm/FnoPl7H7yySdcjjM7/iywYwO12p//eNrbdfSXhVP9PWXlzdg5YL8D7MtB59/b9uN99O8yw74csCD+6GoVJzuOLFebDRRlgfX111/P5SuzL37d5Qa3v7ajf78J6Y2oOgEhQcZ6JturExyNjUBnH+RsQEp3H7BHr+tcXql94AcbHNWuc0kg9iHKFtY7ywZasd4oNhiKfRCzgTRsEAr7EGa9TJ0/tNmH4bnAPmDZB3jnD2Q2SKi7D+mjjxEbSMN6Go93/6lgz8XawoKbzoEsC1JYe9p7tdrbxAKYzj3Q7a+js6MHC7GBTAwLzlatWsUFfuyY9xTr3TuV18V6lVkPOgtin3rqKZyK9tfJXvvJXidrE+sFPN5zsCDsVMXFxXGDCdtrmbKAk72v3W73aQ9AuvPOO/Hvf/+be8+wQJK9ljvuuOOkj2s/5s3NzV2CzFP9PX377bfx448/coMc2SAv1gPMBl11Pt61tbXH7Id9yWW9s6wH/lSxKhhsYYMRWWDNSvmxKx1sECIbQNaOvbbOr5WQ3ox6Ygnpxdrrh7JLeyzQPXphxfFPhH1Add4+NTX1mG3YZUz2Adp+ebW92Hn7BAydL/ey3hlWyulcYeWgOmO9ZSer9cmOEat5yj78uztGpzqRQXtvNeu57IyN6Gc9pu33s94tZsWKFV22YyWxOo/EZ45uU3ugcsstt3A9nGzk//Eu33711Vc4Uyz4/NOf/sSNqr/hhhtO6bHtx/7oc8OCyaNfJzsXhw4d4nr1uzsXpxPEdsZ6kR977DHufX+iovxH90QejQWgLD2EVRNggTA7LqxX/GRYPWeGvd9O9/eUpXSwnmTWI8qCXhb0syC2vXeV/Y6y88Xe+53TTth7j70H2ysWnC52bFiKyTPPPMPdZleAOmOpRSxQ7u5vBSG9DfXEEtKLsbQClrfKLgez3lIWOLG8TNZLwy4hsw9H1qt0qtgHN8v9ZDmJ7MOb9Ua192SyHFiG3cfyEFlvDWsD63VjJZFOdsn1dLE0g//85z+wWq1cbyerVfnPf/6Ty6dkJZ2Oh/Uesw93dmxYoXYWFLBLrhUVFVwvJyuB1d7L1ROslBHrIWXF7FnpLnYOWLoCK4PEyoexXEWGXSJn+ZWszSzQZ7nIrNQRu81yiDvnMx4P2+7bb7/lgiC2786THbCZulggvX//flx++eU4U6zn8XSvFLBczxdffJHrmWfvDxaosvdC59qzzBNPPMGViWIl3ligxgIh9t5iZa/YpAzsfXe8UmLdYcedHRMWcLIZzNhxYe9btp7NcHY87IoD68Flx5Z96WApLaxnsfMXGjbLWfv7gk040hNse9brz3oyO+eh9vT3lAWiLC+V9cSzAJq9HvZlgKVJsMeyqxzsfcNSJ1iaAXtfsNQH1nvK8oJZ+sPpnEfWC19VVcUdC3b82X5Y6hA7nyyg7Yy9NpZffjq9vYScd8EeWUbIQHW8Elvdeeedd7iSVaw0EhvtnJSUFLj++usDu3fv7tiGjXBmo817go2mvuyyy7hR+awcD6s2wB7fuWJA+/Oy8j1sm8TERK4SACu5xdrNRpi3Y/u56KKLjnkett3dd9/dZV37yPTnnnuuYx0b4c5e24EDBwLTpk3jXiMrF8bK/Rxdfuro6gQM24aVKGJtZWWIWIkrNrqbjajvXN6oJ9UJGDbKn1UdYOtZ2SRWYYG1pXP5sfaR5g888AA3apyNBh83bhx3bNnzHz2a/0RYG9nzsfMnl8u5483KSrHR/gcPHjxpia3TfW/1tMSWy+XiSkcd/Tq7Oxesbffccw83+p4dO3YeWXmoRx99tONcdvceaNe5qkB9fT1XSYFVcWDvD1b5IDMzM/Df//63S1WDo6sTMKy8F6sCwI7l8apFsBJfaWlpgVPBylINGTLktH5Pr7vuOu78tldpaPf5559zbWSvq3MVD7YvdrzZ/mbMmBHYsmVLl8cd7/3Qfv7bf0dZ1ZO5c+dyZdHY7wc7j/PmzetSrqu9AgJr39FVFgjprXjsf+c/dCaEkK45m+wyPOuF7evaZ81il9/bqyuQ3of15rLcYpZGc9ddd/X4caynlV0pYD2Wp9LD3xew3GPWQ11ZWUk9saRPoCCWEBJ0fTWIZZfO2SxKbKYjdpmZXfpnl3tZmgALkroreUWCi+WzshJef/nLX7iUk6KiolPOMWU5rCw1oH1Guv6A5TezChcsZ7q9AgkhvR0N7CKEkNPEckJZ3i3Lk20vA8VyeFk9VApge6cnn3ySy3tmX5jYVKynM0iK5T2z3tjTLfPVG7HeV5b7zHLICekrqCeWEEIIIYT0OdQTSwghhBBC+hwKYgkhhBBCSJ9DQSwhhBBCCOlzaLKDbrBC6Wx6P1Yw++hpPAkhhBBCyJlhFV7Z4Eg2k19PJofpDgWx3WABbExMzBmeHkIIIYQQcrLKGKcyk19nFMR2g/XAth/Yo6dVJIQQQgghZ4ZN6806DNtjrtNBQWw32lMIWABLQSwhhBBCyLlxJmmbNLCLEEIIIYT0ORTEEkIIIYSQPoeCWEIIIYQQ0udQEEsIIYQQQvocCmIJIYQQQkifQ0EsIYQQQgjpcyiIJYQQQgghfQ4FsYQQQgghpM+hILYXzB38Q8kPqLPVBbsphBBCCCF9BgWxQRZAAF6/F8/uehbrK9ZzQS0hhBBCCDkxCmKDjM/jY0HyAtw57E6sr1yP1/a9Bp/fF+xmEUIIIYT0ahTE9gKehgYka5PxyNhH0GBvwJ6GPcFuEiGEEEJIr0ZBbJAF/H40Ll2KxiVLILa6MDt+Nn4u/Rn+gD/YTSOEEEII6bUoiA0yHp+P8Af/CL5Uhron/4nMWjFcPhf2NewLdtMIIYQQQnotCmKDjA3kyt20Cby5s6G79hpY3v8IF8pG4afSn2iQFyGEEELIcVAQG2QBfwDKcgUKlq1FaUsLpGPHYNCeBti9duw37Q928wghhBBCeiUKYoN9AgR8xF0/DknDRoO3w478Jicc+w9glmwk1lWuC3bzCCGEEEJ6paAHsUuWLEFCQgKkUimysrKwadOm425bW1uLa665BqmpqeDz+bjvvvtOuO9PP/0UPB4PCxYsQG/Gl4sQekkaEq6YgJA2IyyxyUjJrkepuRRt7rZgN48QQgghpNcJahC7YsUKLhB99NFHsXfvXkyePBlz585FRUVFt9u7XC6EhYVx2w8bNuyE+y4vL8cf//hHbp99AQu2ZYP0UKVHwunWw3sgF4M8ehxqPBTsphFCCCGE9DpBDWJfeOEF3HLLLbj11luRlpaGF198ETExMVi6dGm328fHx+Oll17C9ddfD41Gc9z9+nw+XHvttXj88ceRmJiIvsQ4Mx1KhMASnYSRuS4KYgkhhBBCelMQ63a7kZ2djdmzZ3dZz25v3br1jPb9xBNPcD22LEDuCdbDa7FYuizBIlRJoMqIgMOlQ2hREw6bcuDxeYLWHkIIIYSQ3ihoQWxjYyPXY2owGLqsZ7fr6upOe79btmzBsmXL8NZbb/X4MU8//TTXs9u+sN7gYDLOyIBKEAarUIPoeh8KWgqC2h5CCCGEkN6G3xtyQY+um3r0up5qa2vDddddxwWwoaGhPX7cI488ArPZ3LFUVlYimIQqMVTpRnil8chskOJA44GgtocQQgghpLcRBuuJWZApEAiO6XVtaGg4pne2p4qLi1FWVoZLLrmkY53ff2T6VqFQiPz8fCQlJR3zOIlEwi29SeiEQWjeVQxR6SGsMx1EIPWq0w7uCSGEEEL6m6D1xIrFYq6k1urVq7usZ7cnTJhwWvscPHgwDh48iH379nUsl156KaZPn879O9hpAsfDep+PJglXQqxTwu9RQ1zfgoq27is2EEIIIYQMREHriWUeeOABLF68GKNGjcL48ePx5ptvcuW17rjjjo7L/NXV1fjggw86HsOCUcZqtcJkMnG3WUA8ZMgQrtZsRkZGl+fQarXcz6PX9xZ+fwDP/pKPsQkhmJYa1tHbypXcSg6FxRSN4Q3lyGnMQZw6LtjNJYQQQgjpFYIaxC5atAhNTU1cNQE2kQELNFeuXIm4uCPBGlt3dM3YESNGdPybVTf45JNPuO1ZGkFfxOfzcHFmBD7cVo5tJU24YUI8orQy7j798AS07CiCoSwHu8zFwW4qIYQQQkivwQt0dy17gGMltliVAjbIS61Wn5fndHl9+G5fDdbnN+CfC4YiRCFGwB/A4ad/hKXsZ3x7aQD/vPgV8HlBH4tHCCGEEBL0WIsiol5CIhTgd6NiMCxai2/3VXPreHweJAla8DSJCCszo8ZaE+xmEkIIIYT0ChTE9jKXjYjC9pIm1Jod3O2QYfEQ8PVIaBCixFwS7OYRQgghhPQKFMT2MuFqKSYlh+LrvUd6YzWpkZBK1dCYRChtpSCWEEIIIYShILYXumRYJA5WmVHaaANPyIc4WgVJIBzVlXnBbhohhBBCSK9AQWwvpJWLMSPNwA30YhRJ4RDKjeCVVcHsMge7eYQQQgghQUdBbC/Fasbm1JhhdXmhHRwNsVCHWJMApebSYDeNEEIIISToKIjtpUKVEsSEyLG/shWyCB2EUikMJhmKW6leLCGEEEIIBbG9WFacDrvLWrhSW8JwOaQOFSrqC4LdLEIIIYSQoKMgthcbFafjUgrsbi/k8WFcXqy1uAAenyfYTSOEEEIICSoKYnt5ua1IrQz7KluhGRQFsSgUhjo3Ktq6TsVLCCGEEDLQUBDby42K1yG7rAWqBAPEYgUiTVJUtlUGu1mEEEIIIUFFQWwvNyouBIdqzHCxk6UTQ2VRoKqVemIJIYQQMrBRENvLGTVSGNRS7K9qhTQ+FEJhGFrKDge7WYQQQgghQUVBbB8wLFqLQ9VmqJONEEvD4Sotg9fvDXazCCGEEEKChoLYPiDVqEJBfRu0qdGQijQIqfejzlYX7GYRQgghhAQNBbF9QHK4Es02D1oDAE/Oh9EcgiprVbCbRQghhBASNBTE9gFSkQAJoXIU1LVBGCqHwq1ElZkGdxFCCCFk4KIgto8YZFAhv74NsgQDRIIQNFbQzF2EEEIIGbgoiO1jebGquHCIJSGwlhYhEAgEu1mEEEIIIUFBQWwfkRKugqnNBZ9BD7FIBWWVA03OpmA3ixBCCCFkYAaxS5YsQUJCAqRSKbKysrBp06bjbltbW4trrrkGqamp4PP5uO+++47Z5q233sLkyZOh0+m4ZebMmdi5cyf6OplYgNgQBYosLvCkfESatahqo8FdhBBCCBmYghrErlixggtEH330Uezdu5cLPufOnYuKiu4HLblcLoSFhXHbDxs2rNttNmzYgKuvvhrr16/Htm3bEBsbi9mzZ6O6uhp93eD/pRQI9VKonEpUWmhwFyGEEEIGJl4giImVY8eOxciRI7F06dKOdWlpaViwYAGefvrpEz522rRpGD58OF588cUTbufz+bge2VdffRXXX399j9plsVig0WhgNpuhVqvRW+yrbMVnuyvxe7Sh5pf1yL1GilunPBjsZhFCCCGEnJKzEWsFrSfW7XYjOzub6yXtjN3eunXrWXseu90Oj8eDkJAQ9HWDDEo0WJzgRYZCJNHDUkIVCgghhBAyMAmD9cSNjY1cL6nBYOiynt2uqzt7s1E9/PDDiIqK4nJjj4elKbCl87eD3kguFiJKK4NJqYJEpIaoogU2jw0KkSLYTSOEEEIIGVgDu3g8XpfbLLvh6HWn69lnn8Xy5cvx1VdfcQPHjoelLrAu7fYlJiYG51XxesDW2KNN40MVKHP6wJPwEdGsRI215pw3jxBCCCGktwlaEBsaGgqBQHBMr2tDQ8MxvbOn4/nnn8e//vUvrFq1CpmZmSfc9pFHHuFyMtqXyspKnDcsJbk+B/jhPmD3u4Cj5YSbJ7AgttEGgU4MrUODWlvteWsqIYQQQggGehArFou5klqrV6/usp7dnjBhwhnt+7nnnsOTTz6Jn3/+GaNGjTrp9hKJhEsq7rycN6zXecIfgNlPHemN/eEBwHb8+q/xegXKmmyQROshgRYNpvLz11ZCCCGEkF4iqOkEDzzwAN5++2288847yMvLw/3338+V17rjjjs6ekiPriiwb98+brFarTCZTNy/c3Nzu6QQPPbYY9w+4+PjuZ5etrDte6vynCbYBRHA1D8BUVnAnvePu220Tga31w9/RCjE4hCYK4rOa1sJIYQQQgb0wC5m0aJFaGpqwhNPPMFNZJCRkYGVK1ciLi6Ou5+tO7pm7IgRIzr+zaobfPLJJ9z2ZWVlHZMnsMoHV1xxRZfH/f3vf8c//vEP9DYsB7i1zo7czTWIy9AjecjVEK/5I1CzF4j87bW2Ewr4iAmRo0mthFSkhb0sOyjtJoQQQggZsHVie6tg1Im1NDpweFstzCYHJo+qhrTsB2De84BQfMy2H24vh4gHDPlqC7LxM6548jWoxKrz0k5CCCGEkAFdJ5Z0pQ6VYcwliQiJVOBwfSogVgK533Z7mBJDFShttoOv4ENvU9HgLkIIIYQMOBTE9jJDJkWirtSCpshFQMHPgM9zzDZxejkqmuwQ6hVQeVSos1KFAkIIIYQMLBTE9jIypRgpowzIyZHCL1QC1cfmvEZqZNxPj0EHCU+DhvrSILSUEEIIISR4KIjthRKGhcIfAMr4M4CSjcfcz+fzEKuXw6xRQyTWoaW8MCjtJIQQQggJFgpieyG+gI+0CREoaYiBv/YgYG8+ZpsEvQLlIhnEIiXcpdVclQNCCCGEkIGCgtheKjxWBQglaBRlAWWbu525q9TqAU/Eg7I+AIvbEpR2EkIIIYQEAwWxvRSPz0N0qg5V3hFAyYYj09N2Eh+qQEWzHTylAKF2DepsXafvJYQQQgjpzyiI7cWiUnWob9HC3WYBmoq73BeukkAo4MGnlUPhVaPWWhO0dhJCCCGEnG8UxPZiqhAp1GEK1IqmAOVbutzH4/EQrZPDqtdBzFejoa4kaO0khBBCCDnfKIjtA72xVW0JR6ahPUqMTgaTUgWxWAdzGVUoIIQQQsjAQUFsLxeZooXFoYS1yQpYG7rcFx0iRzlfDJFQDl9JfdDaSAghhBByvlEQ28uJpUIYkkJQ7R8J1Ow7pie23OoGxAEoTTxY3dagtZMQQggh5HyiILYPCI9Xw+ROAGq7BrEsJ9bi8AAqEUIcajTYu/bUEkIIIYT0VxTE9gFhsSqYXTq4qvMBn6djvVQkQLhaApdGxlUoqLPVBrWdhBBCCCHnCwWxfYBEJoQmIgSNjkjAdPjY3lidDmKBBqa6sqC1kRBCCCHkfKIgtg/1xpqQfmxebIgctTI5xGINzOVda8kSQgghhPRXFMT2EWGxapjsEQhUdy21Fa2ToSQghEAghbeMZu0ihBBCyMBAQWwfoTXK4ZdoYDbZAFtjx/oYnRxVVg8g9kNU74LX7w1qOwkhhBBCzgcKYvsIPp+H0FgtTIEhQP2hjvWhSjHEQj78UgHUVgUaHb8FuIQQQggh/RUFsX0tL5aV2mo4fNT0szI4lGKofCrU22nSA0IIIYT0fxTE9rEgtsWhhbu28JgKBa1qFaRQo6G5MmjtI4QQQggZMEHskiVLkJCQAKlUiqysLGzatOm429bW1uKaa65Bamoq+Hw+7rvvvm63+/LLLzFkyBBIJBLu59dff43+QKYUQxmmQ3O9B3C0dpm5q1qmgkikQXNFUVDbSAghhBDS74PYFStWcIHoo48+ir1792Ly5MmYO3cuKioqut3e5XIhLCyM237YsGHdbrNt2zYsWrQIixcvxv79+7mfV155JXbs2IH+ICRGh+ZAImDK79ITW8QTQyiUw1HS/bEjhBBCCOlPghrEvvDCC7jllltw6623Ii0tDS+++CJiYmKwdOnSbrePj4/HSy+9hOuvvx4ajabbbdg+Zs2ahUceeQSDBw/mfs6YMYNb3x/oIhRo9sUCpryOdVFaGZq8QEAYAK/KikAgENQ2EkIIIYT02yDW7XYjOzsbs2fP7rKe3d66detp75f1xB69zzlz5pzRPnuTkAgFLK4QeGsLOtbJxALolWL4xIDcIkKbpy2obSSEEEIIOdeECJLGxkb4fD4YDIYu69nturrTL9rPHnuq+2RpCmxpZ7FY0FvJVCJItBq01rQi1G0HxHJufZRWDpdcDE2zEvW2eqjF6mA3lRBCCCGk/w7sYiWiOmOXwo9ed673+fTTT3PpCe0LS2nordjr0MWEoMUXAzT+1hsbpZOhVSGHLKBGfVttUNtICCGEENJvg9jQ0FAIBIJjekgbGhqO6Uk9FUaj8ZT3yfJmzWZzx1JZWdnrUwqa/Wxw12/1YiO1UlQrVRALNGiqLg5q+wghhBBC+m0QKxaLuZJaq1ev7rKe3Z4wYcJp73f8+PHH7HPVqlUn3CcrxaVWq7ssvT2IbXGFwt9p0oNorRzFPCnEIhXaSkuC2j5CCCGEkH6bE8s88MADXAmsUaNGccHnm2++yZXXuuOOOzp6SKurq/HBBx90PGbfvn3cT6vVCpPJxN1mATGrB8vce++9mDJlCp555hnMnz8f3377LdasWYPNmzejv1CFSMGTqmGprIPW5wUEQhg1UjQLhICAB39Zc7CbSAghhBDSf4NYVs+1qakJTzzxBDeRQUZGBlauXIm4uDjufrbu6JqxI0aM6Pg3q27wySefcNuXlZVx61iP66efforHHnsMf/3rX5GUlMTVox07diz6Cx6f5cXq0VwZBm1rOaBPgljIR7hGBo/ID0mTHx6fByKBKNhNJYQQQgg5J3gBKip6DFadgA3wYvmxvTW1oCi7AebtPyJrViQwaA637rX1RcjcsA0W0z6MePZBRCojg91MQgghhJBzEmsFvToBOT26CDmanQYETIUd66J1MpjlUij9atTbTr9MGSGEEEJIb0dBbB+lDZfDw1PCUVvZZeauGrkKYr4aDQ3lQW0fIYQQQsi5REFsHyUQ8qGOCEGryQm4rB21YkuEcohFGrSW/9ZDSwghhBDS3wR1YBc5M9pILVqbYxDZVAREDke4SooWmRQ8vhCuYkonIIT0TT6PHw6rB06bB26HF3wBj/viLpIKuOos7N+EEEJBbB+mMchRuT8G+F8QK+DzEBGigFfghaDWcVZmPyOEkHMt4A+gsdoKU0UbmmtsMJsc4PN5kCqEEMuE8PsD8HsDcNo98Lr9UGol0EcpYEzSIiRSwW1LCBl4KIjt43mxh9wh8JuyO/JCWF6sS8KH0iKBxW2BRqIJcisJIaR7jjY3ynOaUH24BYEAYEhQI2FYKDehi1Qp6nYKcdY7azE5YaqwYO/qciAAxKbrkZAZygW8hJCBg37j+zClTgKeVIm26lpo2CcAj8flxVqkYqiaVaiz1VEQSwjpdewWN4r3NKDqcAvC41UYOj0aoTGqk/aosqBWphRzCwt40ydHcT24JXtNWPdhHhfMpowyQCQRnLfXQggJHgpi+zD2B10bqUdrrQQaawOgMnA9sbtkCgyBGvWtVUgNSQ12MwkhhOP1+FC0uwGlBxphTNBg0pUpXI7rmUz8Ehaj4pbWejvyd9RhwyeHkT4pEhHJWkqnIqSfoyC2j9NGqNBam4A4lhfLglidDCskSgwTqlFfUQgkzgh2EwkhBPVlFuT8Wg2JQoSJC5OhDpX9dlTYlSRbIxwtVWisr0NLswlOjxcOHx/ugBABhQEBTSRU2nDEhSkRoZYe02urNcgx5pIE1Ba1ImdzDaryWzF8RgylGBDSj1EQ28exP9wF3gigqRCInwi9Qow2tRICgQz2osPAtGC3kBAy0HtfczfXora4FUMmRCI6TXekh9RpQaByB5qLdqKtIheORh8srjDwBDqIRTrIxSKoJQII+R7wfYch8NbAJnBiJz8SFZIUiGJHY1hyLEbEaiEXH/koY/uNTNFxqQkHN1Rh02cFGDE7jsuxJYT0PxTE9oMgts2jgbc+mzuZ7I+4MVQFL98LVFmC3TxCyADGqgzsXVUOiVyIKYsGQaYUAdV74Mr+BqZdh9DcHAe/Pxpy4UxIBAIYBW4EeB74Ai4EPHbwwYNAqoBYNwJ82XTIeEIY5Q4EfMVwF72OQxV6PMcbjujBozF3aAQitUd6d8VSIUbOiUPZgUbs/K4EgydEIH5oaLAPByGkNwSxpaWlSEhIONttIadBqhBBqlHBXNsKvc/LZkHg/pC7RH5IW/lw+9wQC8R0bAkh51VNYQsOrK9C0ohwJI/QI5C/Gq0/fIymQiGc/sGAZhFgEEAQ5kcT34YWVz1cTjv3WJFEyqUYeJ1O+ByN4JkLIG3yQdbmRKgyBvqEERBLRmIUWjDauxNVBdvxdv446JOycEVWNAxqKfeFPmFYGDThcuxeWQan1YPUcUbKkyVkoAexycnJmDJlCm655RZcccUVkEpPPzGfnDlNlB6tpTrozZVASAI3uKtNIoLKokC9vR4xqhg6zISQ81bzlQ2wYqWzRs6ORYh9L1r+sRStZRrYdFPRZtDDa/DA6qkATypAeFwSwo0JGGSIgEKng1gmA5//W3UBn9cDW2sLrM1NaKmuQu3e3cjL/xZymwuJ8aMQnjQPEW1m3IetKCvfixfKpmPc8HTMGxoBiVDApRJMuDwZO78v4cpzZU6PBl9AkyUQ0h/wAqzw3ik6dOgQ3nnnHXz88cdwuVxYtGgRF9COGTMG/YHFYoFGo4HZbIZarUZvV5TdAPO2H5B1YSyQPBM5NWbsfPcXpFQWQ/q3GRgV2T/OCyGkd/N5/di3pgKWJidGjhbB98mTaCvwwmq8ADVyLZyKeggUNsQNG4HI1DRoDRGn1TPqsttRsXsHin5ZCUtVJZJiMhE3ZBb4LSYI+OuxRhGPbPkE3DY1BYlhSu4xLIDd+UMpFBoxRsyKpUCWkH4Qa51WENvO6/Xi+++/x3vvvYeffvoJKSkpXDC7ePFihIWFoa/qa0FsY5UV+z/fgBkTaoCxv4fZ7sELr6/BRSWlaL4pBBdmXRXsJhJC+jmPy8ddtvd5fEh3/gTbT6vhNM5BhSwaJkEFQpMUSB03HpEpg8EXdK3jyj6G2OQsTq8TATZ7AQCZUAaVWAU+7/i9puxxjUUFOPDBu2gqL0VqxnREhQ6DyHEQLcpCvOafg+mj0jE340gagdvpxfZvSrga28NZIEszfREycIPYdqw3dsmSJXjkkUfgdrshEom43tlnnnkGERER6Gv6WhDLPjxWvbIJMwethuTSp7g/7H9+fweu3puPggmVuGrRY8FuIiGkH3PZPdjxfSnEXgtitj6PgNsIk2EKDnvqIY8Hxs2ZhfCEJC6Q9Pq9KDWXckuZpQw11hq0OFvgC/gg4v82SxfL5+eBB7VEjWhVNOJUcUjUJiJZmwwh/9hMuLrdO7Hn/WXwun0YOWkxZGYrBOL1+Fg+Gh7jCNwxNQlSkQAuhxc7vimGSi/D8JkxXK1ZQkjfjLXOqDrB7t27ubSCTz/9FAqFAn/84x+5ntiamhr87W9/w/z587Fz584zeQrSA2x2GnmoFuY6C8K9bvCEYoQaNPDxffCXNNExJIScM+wy/fZviiGtz0XkjnfgTrwYhdCiLFCAaVfNRsrwEfAH/DjUeAh7G/biYONBLghN0iYhSZOEqdFTESINgVaihUgg6tivx++BxWVBs7MZVW1VqGirwJaaLXD5XEjXp2O0cTSG6Id09NQaR43BhcNGIPedt7B19RIkpE9FvOgSLDZvwibU4ekfL8A9MwdBr5Rg7IIkbPuqCHlbazFkUiS9Owjpo06rJ/aFF17Au+++i/z8fMybNw+33nor95PP/+2yT1FREQYPHsylHPQ1fa0nltnzSzlURR8g5cpFQNggfLS9HHFfrEJtYD9+9/wrNCKXEHJOAthtX+RDnLsecfW7YI27FHutdeANDcXli+bDKfBga/VWbK7ezG0/yjgKw8KGIV4df1p/k9jHFQtm95v2Y0ftDgh4AkyOnoyJkRMhF8k7tjMfOoRtr7wAKPUYPnQBpJaDKAppxHLhhbh31hDE6uWwmV3Y+mURkkaGI3F4301/I6SvClo6Act9vfnmm3HTTTfBaDR2uw1LK1i+fDluuOEG9DV9MYgt2WdC85afMGpmGJA6F+vzG+B87wcEzHkY85+HECqjGomEkLMbwG75cD/E+35CvMiGBu0obHOUYOSVc5AxPAGry1ZzwWuyLhlToqYgPTT9hPmtp4qlJbBe3fUV61Fnr8PsuNmYEj2lo6Sg12xG9jNPoaLRhNGjroXGWYMmXQGW8ObiD3OGcQO+2FS1278tRuYFMYhM1p61thFCenEQW1ZWhthYlhTf9Q8S21VlZSV3X1/WF4PYpmor9n32K2aMLQEm/AGF9W3YuvRbJNSVQf/PBUgPzQh2Ewkh/QQbILX5zW0Q7FmNhAgtymDANpkZV9x6OXIdu7C+cj1SdamYlzjvnJf4Y587ec15+L74e25w2MKUhRgRPoLr6Q14vSh+83Xs3b0dmaMWIczrRJs6Gy8JL8GdszKRYlBx0+GyCRnGLUiCNvy33lxCSO+PtU7ra3FSUhIaGxuPWd/c3EyTIASJOkwGp18JZ205d5tNeFAhV0PMV6PBdGQdIYScKTZwasurayHYuwHRMXHI9auxOVKIsdePwLslr6KgpQD3jrwXtw+7/bzUqGbBKsuN/fPoP2N+0nx8VvAZXj/wOpocTeAJhUi6825MvPAS5O7+FBV+L5Tm0XjA+x2WrNqHogYrDPFqJGcZkP1TGTfoixDSd5xWEHu8zlur1XrKEx+wqgZs9i/2uKysLGzatOmE22/cuJHbjm2fmJiI119//ZhtXnzxRaSmpkImkyEmJgb3338/nE4n+jORWABFmBbmRifgtkEhEcJjCINIpEJLUUGwm0cI6Qd8Xh82P/cdkLsXEdHx2O/xYGeGEYoJJqyq/B6XJl2KB7IeQJw67ry3jQWzYyLG4K/j/gq1WI2ndz6NXXW7uPWRVy7ClCsXo2r/jyhztUFqHoMHfN/h9dX7UdViR9LIMOgiFNjzcxn8Pv95bzshBOe+OsEDDzzA/WR/FFj1Abn8t0svPp8PO3bswPDhw3u8vxUrVuC+++7jAtmJEyfijTfewNy5c5Gbm9ttSgKb7pYNILvtttvw0UcfYcuWLbjrrru4mrQLFy7ktmETMDz88MNc1YQJEyagoKAAN954I3fff//7X/Rn2kgtzAVRMDSXAsYMaI06BHiAu7A22E0jhPRxfvY3/pkV8JfVIyIqDrv9ZuwbqYAvfAMSZMNwW+bNXQZXBYtCpMC1adciMzQTH+Z9iMPNh3Fl6pXQX3ghJivk2PLumwgMnol48zjcof4BL/8ixB8vGoZh02Ow9X8VC9InRwX7ZRBCznZO7PTp0zt6Q8ePHw+x+EgCPcP+HR8fz5XZYgO/emLs2LEYOXIkli5d2rEuLS0NCxYswNNPP33M9g899BC+++475OXlday74447sH//fmzbto27/Yc//IG7f+3atR3bPPjgg1ypr5P18vblnFimdL8JjVtWYfRUBZC+AJ/tqkTE8p9QIz6ERf9+OdjNI4T0USy3dO+/30djhRNRYSHYKbJg05BmhBnMWJx+HTJ6ac59q7MVH+R+wOXK3jHsDm6Aq3nDBmz98B2EpcxAvFSMet0hfKJciEcuzgTP6cPmzwuROT0GEUmaYDefkH7Ncr7rxK5fv577yaoSvPTSS2cU4LHqBdnZ2VyvaWezZ8/G1q1bu30MC1TZ/Z3NmTMHy5Ytg8fj4SZZmDRpEtdLy4JWNg1uSUkJVq5cecIqCWyyBrZ0PrB9kSZcjiKXHoGmHLDiNVE6GdqEPMitYtg8Nq6HghBCTkXA50PuM2/CVMlHZGgINktbsWrQYUxNTsFNGX/gZtXqrbRSLe4efje+KfoGz+56Frdk3ILUadMw3ufDtuXvQ5gyF7Gtg3Bx4Du8tk6GB2encgHsgfWVUIdKodBIgv0SCCFnOyeW1Yg90x5KNjCMpSAYDIYu69nturq6bh/D1ne3PatF2z7Q7KqrrsKTTz7JBbMsqGWD0FgP8tHBcmes15d9G2hfWB5tX8T+6LoDSrjqqrjbUVoZTGIJlD4V6m31wW4eIaQPBrAl/30dFRV8ROg02Kiqx3fJu3Hb6IvwfyPu7tUBbDsBX4CFgxZiQfICvHHgDWyr2QbtjBkYc+lCVBetQq1HhVSzASl1P+HD7eUwJqoRNUiHvasq4KP8WEJ6tR73xF5++eV47733uOCV/ftEvvrqqx434OiC1yy74URFsLvbvvP6DRs24KmnnuLybFm6Apt04d577+Wmv/3rX//a7T7ZdLnt+b7tPbF9MZAVigRQhmvR2uyF0WlGhFaJarkSyU1q1LVWclM2EkJITwT8flS++gby84HIEDXW6CqwLrYE/57xRwwNT+tzB3FC5ATopXq8dfAtWD1WzLzkUoxqbsGuHT9BmHQJLrAdwicF6/GLZh5mTYjgJkIo2FGHtAk0oxchfT6IZT2U7YEi+/eZCg0NhUAgOKbXtaGh4Zje1nZsYoXuthcKhdDr9dxtFqguXryYm0WMGTp0KGw2G37/+9/j0UcfPaa2LSORSLilP9BGaGAuiIOxqQiSqCx4I8MhqmtBbclhIH5qsJtHCOkDWOdA/bJ3kHPIjQhtGNbrivBrjAVL5j2JWE04+qrUkFSu/Ndr+16D1W3F/Ouuw/DWFuwrXg1x/Dxcw9+Il3ZqEKefheGzYrHl80KEx6mhj1IGu+mEkDMJYlkKQXf/Pl1sIBgrlbV69WpcdtllHevZ7fnz53f7GDaY7Pvvv++ybtWqVRg1ahSXOsDY7fZjAlUWLLM/yqcxr0OfowmToSE/BmgqAaKyoIkOA/YVwpFfCVwQ7NYRQvqC1i8/x6Ht9QjTxWFrSDF+jfHjjUv+DoOq96cPnAyrXftg1oN4ee/L8Aa8uPyOO+H811PIq96IzKipuJ3/C15ep8efFkzA4AkR2LemAlOuSoVIIgh20wkhZyMn1uFwcMFiu/Lycq42KwsoTwW7hP/2229z5bBYRQFWz7WiooKrONB+mf/666/v2J6tZ8/FHse2Z49jg7pYRYR2l1xyCVft4NNPP+VKcrGgmPXOXnrppVww299pwmVodYcg0FTM3Y7WK+CEG6ixBbtphJA+oG31Khz6cQ9Umhjs1ZZiXbQSL174x34RwLYLk4dxPbL7G/bj64ofEHfPvYj2WJHfuB8B2xRc6/geb6zPR1SaDqoQKQ5trA52kwkhZ1qdoB3rKWV5sSyobG1t5aoAsJ5VNrjqhRdewJ133tmj/SxatAhNTU144oknUFtbi4yMDK6SQFzckULZbB0LatuxSRHY/SzYfe211xAZGYmXX365o0Ys89hjj3FpD+xndXU1V0OWBbYsT3YgUIfK4IECjroqyAMBRGnlaBYGILHw4fF7IOIf6bEmhJCj2XfvRv6H30IQMhoFmhqsjjHiyRmLEKfvf5fTWbmt+7Luw4vZL4IHHi76wz1w/OdZlItDkSTMwNC6b/D1Xi0uvSAGv35agLoSM4yJVHaLkD5bJ7ZzPiurFZuens71pL7yyivYu3cvvvzyS24ShM51XPuivlontt2vy/OQYn4DEdc9imqvEuv//QkizOVIf/4uGBXGYDePENILuYqLUfTEP9GsnYxajRnfJutx14R5mJgciv7MZDfhv9n/xZToKZhYLsWuTz9BeOxsRKoK8IE0DDPmXgmdzY/DW2sx5epBEEtPq++HEHIOYq3TSidgqQSq/11aYikErFeW5aGOGzeOu9xPgktjUMLMSwCai2FQSWBSqCCFGnVmuiRGCDmWp6EB1f9+Ag3qMTDLnVgzOAZXZs7o9wFse2rB3SPuxtqKtdiXxEdG1ihU129Biz0T17hz8MX67VBGK6AOkyFvC81+SEhvclpBbHJyMr755htUVlbil19+6ZiAgFUK6Is9l/1xcJfZHwU0FUMo4MMbEwGJUIOGyvxgN40Q0sv4rDbU/+tvKBIPhV8qwcbMGIyNzcKlwwZOaakoZRTuHHYnvi7+BtUXDkOKUo7S5r1w2afhSvsPeHdTATKmRqK22AxTRVuwm0sIOZMglqUMsMFUbJpZVouVVQ1o75UdMWLE6eySnEXacDnMLl3H4C51jAHgi9CWV0LHmRDSIeDxoPG/T+OAIxQyiQEHx8dAr03BDePjTlivuz9idbRvzrgZHxavABZfgrC2alSa6xDqGIH4ss+xpbIVaRMicHBDFbweX7CbSwg53SD2iiuu4AZc7d69Gz///HPH+hkzZuC///0vHdggU4VK4eUrYK+rYQUfERWmhBMu+KvMwW4aIaSXYMMhmt99C/vKmqCVp6F+UhwqvGG4c2oSdwVnIMoIzcClSZfiTdM3CLviUvib98PUpsYkN7Bn62oIo2SQqcTI3979rJKEkPPrtP9SsYkHWK9r55qsrErB4MGDz1bbyGkSCPhQGXQw2xSApYabftbG90HU5BsQtXIJISdn+f477Nu0ESrVBNhHRGKtXYG7pydDIx/YFUymxUxDliEL78iykZCehrrW3bBYs3CVaws+3ngIQ6ZEojK3Gc21VLaQkD4ZxLIZsFjt1QkTJnD5sYmJiV0WEnxagwJmXiI3uIsFsU0CIeQulmZAvbGEDHSOffuwf/lbEOlnwR0fhh9EGiwaHYOksP5XSut0LExZCL0sFD+M9COB70VVWxEEjskYUfcFNlQ0YdAYIw6sq4TP6w92UwkZ0E6rVgib0pWV2GLTu0ZERAy43Km+MrirJi+SG9wVFj8ZzSoVYm1q1NnqoJVqg908QkiQeOrqsPf5v8ATNgcibTg2xUdjeJgC01L77nSyZxufx8cN6Tfg+d3Po2jOIIR+kY0GcTiGidX4aNfPGHH51RAWtaJodz1Sx0UEu7mEDFinFcT+9NNP+PHHHzFx4sSz3yJyVmjC5chzaRFoyua+ZPhZhQKTBQ3VRRisp5QPQgYiv92OPU/8H5o1o6GTRyN/TCp4fh+uHXdkghnyG7lIjt9n/p4LZK8bmYzWvAMwiybhCt6P+HzTENwydTi2fV2MqNQQKHUSOnSE9JV0Ap1Oh5CQkLPfGnLWqEIk8AkUsNfXAz4vlPGR4AskaD1MZbYIGYhYPvzuZ/+Earce4YoMtI5LxwGLHXdNT4ZogA7kOhk2Ocz1Q67HxwlVMIo8qLEXAq4LMKJmBXY2mBE7JAQ5m6pprAEhQXJaf7mefPJJrswWm/SA9E58AR9qoxZmpw4wVyIyXA1nwAVPaVOwm0YICYI97z+P8uI6xOmmwpo2CD+0WXHntGSEKMR0Pk4gMywTkxIuwKoxAqhbCmCyeDDYG4GiLV9CP0QHS5MDtUU01oCQPpNO8J///AfFxcUwGAxcrViRqOto1j179pyt9pEzrBfbak5CZHMxonWjUcP3gN/opmNKyACzf+0KlK9cjfiYq9BmjMWPEj8uHx6DVOORmRfJiV2ceDFeai1GQ2oVFNW5aGsbi7n4Ht9t24N544cgd0sNwuNUEIoFdCgJ6e1B7IIFC85+S8hZpwmTo+pwBNBUgqiMKchhvbN2MZxeJ6RCKR1xQgaAnJwNKHljCYyRl8OlicGOGBUGG5WYmUYDuXpKwBfgpoyb8Kz5n7iyrBH1zhLESGYjpeRjlKc+DoVGgoJd9RgyceDMckZInw1i//73v5/9lpCzThMuQ45Dg0DTFmjlIliUKhjtKtTb6xGnpoEchPR3BdUHkPvc4wgLmwaxJhmHUiLgkPJw/YSBNyPXmdJJdbh22E34tvYlXLjpIJr4oUgSxOLHDV9hwcXXYt93ZYhO1UEdKgt2UwkZME47m7+1tRVvv/02HnnkETQ3N3ekEVRXV5/N9pEzoAyRIiBSwFrfBJ7PA1+0EVK+BvWmMjquhPRzZa2l2PncI9CIMxCiG47K2BhkB9y4a3oSJEK67H26M3oNGTsPBXE+uN35cFjTMNO2A6tyCxCXoadBXoT0hSD2wIEDGDRoEJ555hk8//zzXEDLfP3111xQS3oHPp8HtUEDs9cAtJRBlRQNoUCGpsOHg900Qsg5xOpB/7L0YYS0hiBSPxFVYbFYK3Tj9qlJCFdRKtGZuCTpEtROTYPXUoRGTyWE7gugy/kYnigZ7GY3qvNbztp5JIScgyD2gQcewI033ojCwkJIpb/9QZw7dy5+/fXX09klOYf1Ys38IzN3RUSGwA0PnEW1dLwJ6aeaHE349LO/I/aQBzHh02EKicUWNTB3RBQyojTBbl6fJ+QLccPI27BtvBqo3wWzVYjhXjG2bvoJKeONyNtaC4/LF+xmEjIgnFYQu2vXLtx+++3HrI+KikJdXd3ZaBc5i3mxZo+Rm7krSiuHFS7WTUPHl5B+iE0r/c6qp5G2uhrh4VPhCEvCPp0MkfEazBtqDHbz+g29TI85c+5Csb4Nbd5COO1jMbFxFQ60NUGll6Jwd32wm0jIgHBaQSzrfbVYLMesz8/PR1hY2NloFzmLZbYsThX8jUWI0slgFvAhbhPA4/fQMSakH7F77Hhj24vI/L4IqpCxkBhTkStToUkvwC2TEmgg11k2PHw4hAvmwlW9E83uRih949C240MYh4ei/FATrC2us/2UhJCzEcTOnz8fTzzxBDyeI4EQG+VaUVGBhx9+GAsXLjydXZJzRKmVABIlrI1WKHkuWJVKKHwqNNgb6JgT0k+4fC4s3fsa0n8sgFiQhFDjEOSKw7BH6ccfLhgEqYgGcp0Ll2UtRtEYI+xN22FvC8UYZyM279+OmME65G2tOSfPSQg5wyCWDeYymUwIDw+Hw+HA1KlTkZycDJVKhaeeeup0dknOER43uEsNM2KBxgL4WYUCnga1plI65oT0A+yqylsH3kLsxiJImhWICh+JQmUcdki9uOmCJBg1NJDrXJEJZbjgqj/B5K9Bc6AUPud0xBd/CmeEEC21NjSUH3vFkhAS5CBWrVZj8+bN+Oqrr/Dvf/8bf/jDH7By5Ups3LgRCoXilPa1ZMkSJCQkcCkKWVlZ2LRp0wm3Z8/BtmPbJyYm4vXXXz9mG1Yt4e6770ZERAS3XVpaGte+gZxSYOYnA6Z8KFLjIBQqYMrLDXazCCFnyB/w4/2c9yHdX4SQQ3ZEhY1Fc1Qq9vL9mDwuGiNidXSMz7FB+lQorrocjuK1MLscMHjTcHjTCiRkhSN3Sy38Pj+dA0J6y2QHfr8f7733HhfAlpWVcakELAg1Go0IBAKnlHe1YsUK3HfffVwgO3HiRLzxxhtchYPc3FzExsYes31paSnmzZuH2267DR999BG2bNmCu+66i8vDbU9jcLvdmDVrFtdL/MUXXyA6OhqVlZVcL/FApQmToeywATBtQmTcDLjhhuNwFTAv2C0jhJwu9vd2+eHlaCspxKiNJqhCx4KXlIpdZj+UmRrMH06zR50vcybdiM+3bEagcRcU/OkYyfsGBY7pUELC5ccmDKOxIoQEvSeW/dG89NJLceutt3KTGgwdOhTp6ekoLy/nSm5ddtllp/TkL7zwAm655RZuf6y39MUXX0RMTAyWLl3a7fas15UFt2w7tj173M0338ylN7R75513uMkXvvnmGy4wjouLw6RJkzBs2DAM5AoFFocKflMxYjRiWOAGv54GHRDSV7G/xd8UfYOSyoOYtLoZQlU6dAnJ2NosRl2kCLdNTaKBXOeRSCDCpFsfhbMxF42BavA9F8C9630Yh+u56WjdDu/5bA4hA8YpBbGsB5bVgV27di327t2L5cuX49NPP8X+/fuxZs0arFu3Dh988EGP9sV6TLOzszF79uwu69ntrVu3dvuYbdu2HbP9nDlzsHv37o5BZt999x3Gjx/PpRMYDAZkZGTgX//6F3y+gVu3T6GVgC+Vw+LRIwr1aBXxIbGKucEghJC+Z1X5Kuyq2o4FmwKw+kIQFZWEXV4jcpUB3Dk3FQrJac0oTs5AbPggaH+3EK1lv8BplyLDHcDu/C3QGeUo2EmlJwkJehDLgta//OUvmD59+jH3XXDBBVx1go8//rhH+2psbOQCSxZodsZuH6/WLFvf3fZer5fbH1NSUsKlEbB9szzYxx57DP/5z39OOODM5XJxJcM6L/0JS/HQGuVoEQyGqLkQLn0IV6Gg1kqTHhDS1/xa9SvWlK3G4pxQVNdZkWBMR6lhMPZ7XFh0UQqidfJgN3HAmnzR7eAZlKjxHITPOR4JpV+BFydB5eEWWBodwW4eIQM7iGXTzV544YXHvZ/ls7Je2VNxdA7tyfJqu9u+83qWs8vyYd98801uANhVV12FRx999LgpCszTTz8NjUbTsbCUhv6G9Qa0+GMB02EIk+IgE+pQV10Y7GYRQk7BrrpdXBrBLc1DUbYnD8nhI2EfkoEd1TZkTIrE6PgQOp5BJOALkPX7R+Ar3YIWdxtCPOko3PUZooeEIHdzTcfnFSEkCEEsyzU9uie0M3ZfS0vP5o0ODQ2FQCA4pte1oaHhuM/BBo91t71QKIRer+dus4oEgwYN4vbdjuXPssexFIbuPPLIIzCbzR0LGwjW3+gMCrQ6dFyFAsPgGIAvRPNBqlBASF+xr2EfN5DrVtF0VP6wDrFhWZANTcK6fBcEQzS4Yuyxg2HJ+ReZPAy6aVNR1boBPnsShjVlo1bQjLZmJ+pL+9dVPkL6VBDLLtGzgPF4WODILu33hFgs5npKV69e3WU9uz1hwoRuH8NyXY/eftWqVRg1ahREIhF3mw3mKioq4npk2xUUFHDBLXvO7kgkEq5sWOelv9Ea5HC4ZXDaPEhSu2ALOOApNgW7WYSQHjhgOoAPcj/AzYb5aHz/a+h1wxCeEoFfKuVoNIpw+7xUCPg9rwxDzq0xix+E2tqEal8ZeJ4ZcO5ehqjhocjbUgOfl0puEXK2nFL2P7sUwqoQsKDveLmlp+KBBx7A4sWLuSCUBagsBYDN/HXHHXd09JCyKgjtg8XY+ldffZV7HCuzxQZ6LVu2jMvVbXfnnXfilVdewb333ov/+7//Q2FhITew65577sFAJpIIoNLL0BoYjChvJQ7yvOA10ohZQnq7nMYcvJfzHm5IuBKulz+FUJaM6Gg9NiMRJXDgziuGQkkDuXoVsVKNQdfehuyP3oAjejFSBAHkVe9AmHgwSvc3IjkrPNhNJGTgBbE33HDDSbe5/vrre7y/RYsWoampiZvCtra2lqskwAZjsbJYDFvHgtp2rB4tu//+++/Ha6+9hsjISLz88stdprpl+aysd5Ztk5mZiaioKC6gfeihhzDQaY0KNNcmw9hSBIdShpBWKTffulxEA0EI6Y0ONx/GskPLcF3yVRC+/T2anEoMjjSgOH4EDh6w4PIbhiBKKwt2M0k34mfMR+3alSjx7sUQ50QY8j+DfOYzKNpZj+jBOkgVR64eEkJOHy9AmebHYNUJ2AAvlh/bn1ILKvOaUbkzFxP0X+Kz8umIOlQM41OXICkkOdhNI4QcJb85H28ceANXp14F/WdbUJBXiYyIofBMGILPN3kw6MIYLJhw5As/6Z1sJYXY/Ne7IU24AnpFK/IjRIg3XgahSIBhM/rfAGJCznesdVrTzpK+W6HAbFPAb66DLiUMEqEatcWHgt0sQshRCloKuAD2ytQrEbmuELmHCpFmyIRinAHfbvdBlanD/PE0kKu3UySmIHHqPFS0sUFeCRhUvwttGgtqi1rRWm8PdvMI6fMoiB1gkx4IxCJYJOlI0bbCCz/MB6nMFiG9SU5TDl7f/zoXwKYcakP22nVIixgD/VA+PssLg0snws0LBtOMXH1E4lU3I8XqQXmgEPDMgH3XMkSkU8ktQs4GCmIHEG7SA4McLcI0RLiK0cYqFJT1rCQaIeT8lNFadnAZrku7DkNqxdi2/BMMipqE8Fg7vmsbika7B7dcnwmJ8LcSgqR3E6jVSL3qVvjrd8LuECPGFkCVeQ8cbW6uR5YQcvooiB1gdEYFWrxRkDTlwiYGhC2/TRhBCAmenbU7uTJaN6XfhCFWDba+9iriIyYj0mDGlrAJKC2z4eobMhCq7r46DOm9tNNnYHBYIgr92fA7x0GbuxwhaSrkba2Fz0Mltwg5XRTEDjC6CAVa2pQIOMzgh8ig9KjQ5GwKdrMIGdA2V2/Gp/mf4veZv8cghwa//vspGMLHI0pvRlHaeOzZacGc3w1CYlT/GWg6kPCEQsTfcDsS6ith8rZA4RqKyoIvIVGIULy3IdjNI6TPoiB2AA7ucjn9sKsyEBIJyAU6VNYcDnazCBmw1las5aaSvXv43UhwabHpqcehDx2LuFA3LKPS8PMGN7JmxmJUBtUW7cukqalIGjMDDZ5s+BzxiK3aAb7RhZK9Jjis3c8mSQg5MQpiBxiBkM+lFDTx05GobwD4Ipj2Hgh2swgZcPwBPz4v+Byry1fj/0b8H2K9Gmx+4q9Q60YiLlQI/kglPtmoQnyGHnOmxwe7ueQsCL3yKgwz+1GKQng909G2bxn08UfSCgghp46C2AEoJFKBJncUDO5iWAN22A/3wz+gfh8CHic8dhdsrU7YWl2wW9xw2T3w+ykHmASXy+fCWwfewuGmw/jjqD8iyq/B5n88BrkqHQnhGigzbXgzezB0ITJcdcVgOl39hDAkBNGX/g5qWzFcDhHCzR60+vLQUGZBc60t2M0jpH/P2EX6B32UEpW5UvAVfjiFLghMPvRZTjNQsw9oKYOzoRb1VW60WsRotSpgc8sRCPDB4/nB4/Ph50m4nmcIxZAoJJBq5FCFqaE0hkEVroJSJ4VMJaLSReScMrvMXA1YMV+MB0c9CLHTjy1/fwhiaQqSIqKgTi7Aq4WzIPUFcMNNQ8EXUF9Df6KeNQuDt27GJv4hpDjHQJbzEWRjnkLuphpM/F0y/f0h5BRQEDtA82LdTh/sxuGQqpxAowI2jw0KkQJ9gscJlKwHyrci0FiE2sAIVFmT0dgWDZ1RiZAkJQZHqKHUyyGWAEKBHzyPjRvM5re1wN1sgrOpBo7mNrSVOdGaI0KlLxQ2XwgEUinUoQqoI0KgjjZCHa6ESi/l0jAIOVN1tjos2bcESdokXJN2DXg2F7b+/S8QChOQEpsEbdQOvN5wGYTNPiy+azhkcpqatL/hicUIW3Q1hr7xImrUrdA5BqOh6gdIfDNQld+CmMEhwW4iIX0GBbEDOS9WmAFDyCpYm+NRaSrG4MhM9Ppe1/yfgIJfAHUkTMrpyGv4Hbx+EWJH6TE0VQuZUnycB4eBpwNYdU3Z/xYdW83Ki7ksgLkavuZK2GpqYKlrgiXfgZq9Chz2hcPDU0Kpk0Bt0EATbYAqIgzKECmkSuq1Jac2icF7h97D1JipuCjhIngtFmz6218gFiYiJTENOv0qLHMsgr/Ki6tvHcqlEpD+STp0KCLSRqK0oRgBx0iEV6yEe9w05G+rQ0SiBkIx1QEmpCcoiB3IebGN0UjSlCOfn4qaPXt6bxDLAs3idcC+jwF9CtxjHsSBQwo05dmQkmVAXKYeglO45Orx+VHUYEWxyQpTmwuNVhdsLpYnG80tcr0AmighjEIrYviNCHM1QtBSB5vJipZiL8pdWtj9R3ptFTo5lOFaKCPCoQhVQ6ERQ66RQEQfQqTTAK6VpSuxrmIdrh58NUYbR8Pd2IiN/3gMatkQJKSkQ6f8Eh9hEexFfiy8Ng2GKBUdv/4+8czvrkTmE39Htr4cRudE2PPehcJwB4r2NGDwuIhgN5GQPoGC2IGcF5vXjGGRyXBkt6H1gBm4GL1PWx2wfQlgawLG340mXir2raqA1sjD9OsGQyzt2VuYTehwqNqCDfkNyKuzQCoUINWogkEtRXK4EkrJb/uxu31otXvQbBdhTasUVS16OFxuJGq9iAuzIjLQiGifCe7mFngsNrTVu9G8WwC/KAxeQRT8Qi3ESjnkISoowvVQaGVQaMWQqyWQa8Q9bjPp+9rcbXgv5z20OFu4AVyRykjYKyux8al/IEw7GvGpg6EVf4wvhAvRnCvEpZcPQlwKXU4eCESGcOhmzUHszg1w+iOgabLDHl+O0v0GxKSFQKGhSS0IORn6NB3IebEOLxy6MfDx94Hf4EKvU50NbH0VSJgMTPsLSvOsOLytFGnjIxA3VN+jARAseN1W0oSVB2thdXoxLTUcC0ZEIVonO+njPU4nagoPo7JhD6rzc+F0uOD0AHkuP7w+QCkVQamQQK1WQqgEXK46uB0lEAsDEHoAQWsA/goV2hRxaEA47D4d3AE5RCoV14PLemzbe27ZT4VWApFEQAM7+olScyk3hWy8Jh63jr4VMqEM5twc/Pqf5xEXMQ3RaUnQ8N7BN+KFqDykwryLk5CSGRbsZpPzSD13LqJ37MB2ZQXCHaMgyn0PsvTHcXhbLbIupLJqhJwMBbEDPC+2kZ8GuXg1vHYDPD4PRAJR70gfOPQlkPcdMOZ2BOImoHBXPcoONmHc/ESu3T1RZ3bi/W1lXMrA/OGRGJugh7gHA7TcDjvy1v6CovVrIXE4oXd5MTo6Dqr4VIhCQ8FTKmBxelHdYkdVTTNKaxuh9ToQ4/VC43GDr1PBodfBJhWh0WlFoK0SYSENSNJ6EcJrBdw+2G3RsLnjYDcZ0RjQw+aSwu3wQSgRcD0wqhApNOEyaMJkUIfKaGBZH8K+OG2o3IDvS77HxYkXY3rMdO6LSf3mTdj+zrsYFH8RIjOioHS+ga8EC1CWq8WsmfFIH2UMdtPJecaXSKC94gpkfvwBCqV6SGxJcJjXwmeZgKZqK3fFjBByfBTEDmCh0Uo0NjmgNQgQKNOi1lyF2JCE4DbK7wd2vgnUHwRmPYmANha5m2tQW2zG+MuSuOCuJ0HEqtx6fL2nGlNTw3DvjBRIRScfKBHw+5G38nvkrVoJeZsVIzKzYJw0GZK0NPCEUngbHfA2OeAzu6EKeBAtt0AbZYE11A6Lw4YGrwMlbgf4bjukTjtULW0IsbTB53ejuVGESjkfHoEMIeFGRBq1iNC6kMTfA35TASDywRs9CDbZYNjECWjzSbnakQU76+F1+/4X1Mqhj1RAH62EVNELvmyQY7Q6W/FR3keot9fjD8P/gERtIre+4NOPkffLr8gYdAXCh2khNy/F55iP6vxwzLogDsMmRtLRHKBkI0dC/esmyN2tELYkgF/8MzB6InI212Dy71LA45/8ihMhAxUvwD7xSRcWiwUajQZmsxlqdf+dq9xssmPb1yUYnbITpV/5YblUhQtmLQpeg3weYOsrgLkKuOAxQB6C3C01qCs2Y+z8xB7liLm9fnywrQy5tRbcPT0ZSWE968mwmhqw9ZUXYC8vR+bEaYi9fCF4YjVc5Wa4y9u44FWgFsOp9qHIV4bDzkI0BZph0BgRrjAgXBEOJU8BUUCM5mYPqmva0GSyge90I8rnQajTC4nDDV9bI2zeVpilXtiEbnjkfIRmJCM1IQ4RcjdE5mKg4TDgbAVCUxAwZMKhzoDZHYaWBgeaa2wwmxxc+gHrpdFHKRAarYJYRt9Hg4n9Gd1dvxuf53+OoWFDsTBlIeQiOfweD7Jf+g+aCxoxJO0ShA4XQFj7Oj7xXo7mciNmzojH0AkUwA50npoa1D39NPaGxUPq4UEYVQ+e8kYkjghHXLo+2M0jpNfGWvTJN4Cxy9RCMR+B0LFw+1ehcXcNMCuIAeyvzwOuNmDmPwCpGiX7TKjOb8GEhck9CmAtTg9eXlMIlur6t4uHQCs/XrmtrsrXrcHOj9+DMSQMk598BkJZKOwHGuGuaoA4UgFpqg6tehl+rF2FvQ17MTR8KEYbJmOIfgjEAvEJA5uqFgeyy1vwVWkTrBY3sjTDMNjVhISGWvhrzQiw7IKtfhRvPYgDAjsQpUB4xjwMHpkEtasSvLoDkB/+HnIeDxHGoUDWSHj0mWhuFnCXG4v3mLB3dSWX42yIV8OQoOZya3uSL0zODpPdhE/zP+VqwLLar8PDh3PrbfV12PzMU1B6I5E56kroM9uAkrfwtuNyeOsiMXtWHNJpFDphg7wiI6GcMgWDCvJQb4uAuyEfkshq5G/3IzJZy+XKE0KORT2xA7gnljmwvhJCIR+NK15HM9+FBS8/e/4b4fcBW148UoGA9cCK5agpbMGB9dVcDqzWIO9RAPvcz/mI1Mpwy6SEHuW+Moc//Rj7V36HEXMuRuLFV8K2sw7eBjukg3SQDtHDLwV+LP2Ry3FkpZFmx81GmPz0Bt9Utzqwu6yZC2rrLU5kRGqQZZRhUEsd2g7mwlJUA49dgIBECx+fD4fMBV6MGuFDEpEcr4HYnHdksFtLGddLi6hRQFQWnIIw1Je3cekHjZVWrn4tC2bZhx/Lq6WA9tzw+r1YU7EGq8pWYaxxLC5JuoTrfWVqtm3BrmXvITFsIiJHjoE2JgfOkh/weusVULSGYfbFyUjMDD1HLSN9kd9uR+0//oE8YyQEJhlEIRuAxL9DF65C+uSoYDePkF4Za1EQe44ObF9RV2LmRsKKCz+Hqw6Y8NqfIRGcx9IuLJtl+1KgufhID6xEhZY6G3Z8W4KsufEIiz15vcw2pwfP/5IPo0aG309JhKAHOWQs/3X/K/9Fwd7dmHDD7QgJTYdjvwmSZC3kI8LBlwhRbinHh7kfcr2t1w+5HkbF2Rt4U2s+0kO7u6wFNa0OpEdqMCpehwxlAPyifNTtykZjQRXcLIqW6CEVa+HV8CFM0CEqNRpGEatdmw0eyx1WhgMx44DYcfAqotBYbeNSMOpLLRBLBYhM0SIiRQu1nornny1FLUVc7yufx+dqvyZojuSSez0e7HtrKeqzi5CWMBth04ZByfsBzWXZeM30O0S7Q3DhZcmITOGm2iCkC9u2bWj65hscFsfCwa+HbpABNusUTLg8mRvkSUh/YqEgtvce2L7C6/Fh1ds5SIhpgPWnQkjvzcDQQZPPXwOy3z/Su8gCWHkInDYPNn9eiKQRYUgYdvIeT4fbh2d+PszVe+1xABsIYN9z/0ZxQR6m3PUgxHUqBNx+KCdEQhR+pCdta81WfFHwBWbHz8as2FkQ8M/d5bwGixO7y1u4oLay2c5NmTsqTofhMRpIG+thO3QQJds2o7W6FR6hBmJJKOSKUPBCFJAlhiAixI8Qfx4EjTsBuZ4LZtniU8bAVGlFTWEr6ssskKvFiBqkQ3SqjuutJaeu3lbPVR3IbcrFvIR5mBYzDUL+kays5tIS7Hj5Rei8MYgdPBH6mXEQl7+OsvpWfFB3GVIlWsxaOIhGnJMT/m1qeO55NKgUaC2TwqFaDfXwP8PnUmHC5Ul0VYX0K5b+EMQuWbIEzz33HGpra5Geno4XX3wRkycfP4jauHEjHnjgAeTk5CAyMhJ//vOfcccdd3S77aeffoqrr74a8+fPxzfffNPjNg2kIJbZ8V0J9JFSNL/3PWqTizD/vufOzxMf/hHI+QaY/U9AZYDf5+faIlWKMXxmzEn/YPv8Aby8thDsDXzPBckQ9mDWLvZ2z3v1JeTs242pv38YwjIRRFFKKMdFgCfkc/f/UPIDfq36FbcOvRWpIak4n1g5MBbMZpc3o6zJjlSDiuuhHRGrg0oQgKuwEHU7tqNk305YzR7wWECrDIdKFQlehBIqgx+RwmqorTvBU2qB2PFAzFh4lTGoL2tD1eEWLpeWVaaITgvh0g5OZbazgVx1gM26tatuF8ZGjMXchLnQSDTcfT6vBwc++RDVv+5BimESwseNhHo4ENj1AraZI7CpfgaGhmtxwcIUKmBPTspdWYmGZ59DXvxg+Jo8UEQWwK26G4NGGxEzhCbCIP2Hpa8HsStWrMDixYu5QHbixIl444038PbbbyM3NxexsbHHbF9aWoqMjAzcdtttuP3227FlyxbcddddWL58ORYuXNhl2/Lycm6fiYmJCAkJoSD2BEoPNKK+1Aznmq9hFuRj/n/fZPMi4pyq2H5kJq4Zfwf0SdyqnE3VaKqxYeLlyRCITh5YfbKjArm1ZvxlHqt127MxisXL3kL2lg2YeNW9kDYoIR9pgDQthAuY2fSgLH2gxFyCO4fdeVbTB05Hk7U9oG1BSaMNKeFKjI4PwchYHTRyEXwWC1qys1G6YzOqCwrg9yugUBshVRkh0erBD/VBJ21AtDcHMq3sSEAbNx4OoZELZtnicXoRyXpnB+u4y5WUP9uV3WPH6vLV2Fi1Een6dC7vNVwe3nF/bc5B7Fn2DkLdEYhOHI2QuUMh5e2Cc/eH+NwyFS2tmRg1NBxj5sTTVMSkx1q//BJtBYUobAtBszgHxozpsJtTMO3aVJrxj/Qblr4exI4dOxYjR47E0qVLO9alpaVhwYIFePrpp4/Z/qGHHsJ3332HvLy8jnWsF3b//v3Ytm1bxzqfz4epU6fipptuwqZNm9Da2kpB7AnYzC5s/CQf4pqN8JkqMeeJm4DwNJwzrITUhn8BE+/lBiYx7HL3vtUVmHRlz3qr1uc34Ju91XjsoiEIU/Ush7d+5Y/Y+PnHGH3p7VDbw6GaGg1x9JGcW/Zr8MnhT1DcWoz7su6DWty7euBbbO4jObTlLSg2WY8JaFn7bcVFKN+wHhUHsmFpsXG9sxKNEWJFGPyaACTKBhhwGDF6FYQJExGIGYdmuxZVeS2oLW7lpsVlwWxUqg6SAV6yq8nRxA3m21KzBfHqeMxPno84dVzH/bbWFux55y04DjchMXw0QidkQjU+DPwDb6G29DA+NC1AqN+I6XMSkDAslL4ckFPid7lQ9/gTqEtIhrXID5vuZyiT/gFtmAZDp0XT0ST9gqUvl9hyu93Izs7Gww8/3GX97NmzsXXr1m4fwwJVdn9nc+bMwbJly+DxeCASHcnze+KJJxAWFoZbbrmFC2JPxuVycUvnAzuQsKCRLWJ+HPyNbphyvkfYuQpizdXAr88CI67vCGBZHuz+tZXImBLVowC2tNGGz3ZV4v5Zg3ocwNpyc7H9s48xZMLlRwLY6TEQRx6pIcsCwC8Kv0BBSwHuH3l/rwtgGZ1CjJlDDNzSandzA8J2lDbjk50VnQLaeKTf+nuks2NqakDFujWo2puNqvydUMnCEAiJQ7V8MgqbvfDXFkIj+Bnx4TJkJkxD+ojxqG2QcgEtG+gXHqdGdJoO4bEq8AdQukFlWyXWlK/BvoZ9yAzLxD0j7uGmje08m1vOd9+gZmM24qRpSMmaiJB5mRDz8uFZ8xw2N0Vim+kGDDaGYNqlST2qrEFIdzN56a6+CoF330WrIQ0OywhIzF+junk+YtJC6H1FSLCD2MbGRq7H1GAwdFnPbtfV1XX7GLa+u+29Xi+3v4iICC7FgAW1+/bt63FbWK/v448/joEscpAWDUUxkAhNKNz7PcLGtQIy7dl9EkfLkR7YlNlAysyOAJIFsOFxKq4H8GTsbi+WbijippEdZDh55QLG09iIHS89j9DkcTDKh0A9MxYiw29T1/5c9jMXtNyfdT+00rP8ms8BVv/25AGtDoMWXcMtHpsNVb9uQOXO7agt3QuZQANtaCL4qmnIbfFgc0UeRKKVMIYLMDhhKsRjp6GqUoCcTTU44PEjapCWOzf9Nd2ATbd8oPEAtlRvQZmlDOMixuGxcY91KaXmcTpRsHEtylduRIQ/Gpkxs6C/aARkSRLw9r2LhoLd+Mx0EXjuWFw4NYabwIBN7UzI6ZINHQrJ4DSkuL0oazGirHYj4tJacXBjFSZdQTN5EcIE/Zrh0R+KLKg50Qdld9u3r29ra8N1112Ht956C6GhPa/B+Mgjj3CDxTr3xMbExGAgYaPWC3bWQQ0nTCYDcOhLYPQtZ+8JPE5gwzNA+BAg87dZwUr3NXLpDCMv/O1S7fGwc/3uljKuFuyc9J7lqwa8Xhz6z7NwKyOQET0N6mkxXQJYFryynMcHRz2IUFnfq9vZ04A2Ye5F3OLzeFCzZRMqtm5CZfHPkAQUiDEMgkQ9Ey4L8FNFEeziNdCG+pE8eAxC1XPQVO3H9m+KIVOJucoGLKDt69PechNRtFVhW+027K7bDYVYwQWvN2fcDKX4t1neXHY7CjeuQ+WarQh3hSE9bBJCZgyFIisS/PK1cH23ApuahiDbdBOSYvSYcXESdMbf3l+EnAndlb+D8/HHoUwfg8jCcWiteA3ykEdRntOE+KF97+8VIf0miGVBpkAgOKbXtaGh4Zje1nZGo7Hb7YVCIfR6PVexoKysDJdccknH/X6/n/vJtsnPz0dS0pFBRJ1JJBJuGchY+aWQCCWshQH4rXqgZD0w+CJAdRYGN/m8wOb/AhIlMOb2jkFjrQ12LnBmU8qKxIIe5cGyVIK/XzKkxz2CNZ9+iooWK0aNvBHqCTEdObBMtbUaH+R+wNWAjVL2/WLiPQ1oY6ZdwC0+nxf1O3egfPOvKC34CQKPCCkR6ZDr5sFr46OqpgTrhU+Br3cgIWUINIKpqK/gI39HHfTRSu6LD6tu0JNz15tKZB1sPMhVGTA5TMgyZOH2YbcjUZPY5T1lMTWgYP0aNG07jHBvONJCxiDkojQoxyeC37QfvjWvIKdKjF/qfweJJAzzL0tA6rCwftlTTYJHoNVCM38+eKvXwKpKRl1LElTajcjfMQ3GRE2f/zJJSJ8NYsViMbKysrB69WpcdtllHevZbVYSqzvjx4/H999/32XdqlWrMGrUKC4fdvDgwTh48GCX+x977DGuh/all14acL2rp4r1sh3cpYK2RQVb5Hgo9i8HJt1/ZjtlPeW73gYczUdqwQqEHfVp966uQFJWeI96rtgMV19kV+GeGSlQSXv2h9uRk4P969ZiSPpV0I6O52bhamfz2PDG/jcwK25WxzSh/UlPA9rI8RO5hX3ZM+3bg/JfN6AwbyX8Lj+ijRkYFHoJUC1EQ0MlcgTLUaqsQagxGhGu8ajdmgTBWhmMcVpEJGlgSND0uukxWY9rqaUUB00HccB0AE3OJqTqUrn6riMNI7tM7OFxOVGVl4OKjVvgK22DIWBEmn4kdNOGQDkuEbzGQwhs/AeKy1xYb5qAFq8BoyZFY/LUWAj7UCBP+hY2Ha19+w4kRIrBz03B4ZpVSEnKwqGNVdyEMPTFiQxkQU0nYJfwWYktFoSyAPXNN99ERUVFR91Xdpm/uroaH3zwAXebrX/11Ve5x7EyW2ygF8t/ZSW2GKlUypXg6kyrPZLjePR6cixjkgb7FeFQW1xYy0vBpdXvAE3FHSWwTsvBL4C6A0dqwYp/C1ZZvqVULkLyyN/KFR2P3x/AO5tLMSk5DIONPRt05Wtrw/6lryI8agrChg6CfHhYl8Dmo9yPEK2KxoXxF6K/62lAaxg5ilvY8WnMyUHFhrXIy10Jl9WJ6PA0TAufjdkeJWzWBlTwD+Cg9GtUK0WIaBgOXXkapC4NDDFaJA2KRES8juvdP998fh+qrFVclQlWKq2otYgL0NNDj5THStOndQlc2Ws1lZeiYtdumLNLoXOpEcnXQxmXCO2s4ZCmGsCr2gb3T++ipEKInc1ZqPLpkTzciGsuTIRcTj1h5Nzi8fnQXXctGp55Fsqs6YjNm4Tm+iVw2v6E2mIzN700IQNVUIPYRYsWoampiasmwCY7YIHmypUrERd3JD+SrWNBbbuEhATu/vvvvx+vvfYaN9nByy+/fEyNWHJ6WC+acXgyvHW1KN++D5g5D8h+F5jxWw/qKfXAHvwcKFwNzPgbNxtXu5qiVtSXmDH5qkHg9WCGrVW59bA4Pbh8ZM8v+Ve+8w4cwmgMHjQc6mmxXXorNlVv4kahPzL2kQHXi3GygDaLmylMi7CMDG4ZGQigtaQYlevXoejQWpibmhGuS0KscSQy5HMgEDvRyj+MasGPKNB7UdxkwLbVEdA4IyHXSKCNlnI9tbHxBkSojWd1SmPWm87SA+rtR5YKSwU3MEvAEyBRm8ilCEyLnsZNCdt5xjU2NWxjeSlqc/LQml0EmVUMbUCLUHkUNJPioZo5AkKeBd78dah9/xDKWyKQ3zYJ5UIVYjPDcdfsRGgV5z9AJwOXODoaqlkzwd9/EG3qGNQ1xkCl3Y6cXwXcxCViadCHtxASFEGfsas3GmgzdnXWUG7B/ueXo5ZXiSufeQTSDU8AIUnAmNt6PgECe0vt+xgo2wxc8Big+a2uod3ixqbPCpA5PYa7BH0ytWYHnvg+Fw/OHoTk8J5VI7BlZ2P7mx8iYdCliL1pAoQh0i55sP/Z/R9uMoMUXUrPXs8A0B7Q7qloQWGDFVFaGUbEarmANjZE3hHss9JdlevWombvbphqqiGVG2EwpCJEEwsZTwyp3AJPoAD1/BqUiNRocIXDbtbC7RbCqbBCoPNCFsqHNEwAlVoGpUjJDaRSCBXg8/ngtf/HO/LTG/DC6XVyS5unDRaXBWaXGQ2OBljdVm7WLDb5gEFu4PKak7RJiFBEdPlywtXQbWmGqbQUpn2FcJaYoHTLofDKIJR4oU41QDUtE5IwMax5O9GYk4/GOi9q3YloEqtxWKRAamYYLhkX2+OSboScbWyQav3T/4YvYRDKDzmwT7USaXH3QWeMwvCZx04OREhv1+cnO+itBnIQyy7d//TwexBbW6D78w0Yxa7C//IXYMilRwZ69aQKwe53gPpDR3pgOw0MC/gD2PZNMVQh0h4V7GbTyj69Mo8rpXXl6J7lM/usNux5+BFItZOQdN0FkA/+bQSvx+/BMzuf4XJgL068uEf7G4isLi8OVpmxr7IVh6rNkIkFGBajxYgYLVKNKoj+VzfWZ7Wicfs21G3fClNZKdp8gFwVBb02DmpFOBfUiiU2SEU1CPhMaBTrUc3XodGlgc0mQ0DOg1/thkduh0tmQ0DhPrII/Aj87z8hTwipUMr14LL6vdwiUSNMFsYFr3LRsXVY2TSwzdXVaC2qRFtpPRyVjRDZ+VAEFBB4nZDKfZAPioBoaDLcsMNSUgZztQktzTx4RTrwQ9UoEmuRFxBgbHoYVwkjXP3bFyFCgj0lrWXMTNQcrEJVxFYYFX/GyDnxXG1nQvqSPj3ZAemd+HweosakwLWxFNsP5mDUJVOAKX8C1j0BsNJDCVOO3yPbXApseelIfdlZTwIKfZe7i7Ib4Hb6kDYxokdt+SWnDg6PDwtG9DyNoPbjj+AVJSM0KxWy1K7Pz4rYC/lCzI2f2+P9DURKiRDjk/Tc4vH5kV/Xhr2VrXhvaxlXpzfNqEZGlAbpUWoYZs7iloDPB3dFBZr270VTziHUFm+H3RdAQKqHXBoOjTwJcpECsTwhEvh2CASVCNg98HjVcLbpYEUU2txSuCCGWCuHVKOAVCGETCnmRmBLZUJu8JRQJACfDwSsbthNrWhtq4S9vgWuRgtcTW3wtlgBuw+SgBSCAA8CtxUKgQACvQ5enQxWMR91NhfsZTagcBvkch/UBhVkgwajVhKOLc0euP0BTBkUhsVp4VD3cBAhIeeDOCYGqgvngL8rG9aQFPCaYuBUbsPBjSJMvUpBAwzJgENBLDlG2sWjsXtjITw7DsB30WQIQpOPTBG74w2gdCMw8gZA97+6rqwjnw3+YiW5Sn8FhswH0i9n0XCXfTZWWVG8pwHjL0/mApGTqWqx47t9NfjzhakQ97BovDM3FxX7KhCaOBmGeeldLimz3MlVZau4KWU750eSE2O9rixgZct1Y2NR2ezAoRozdpU1Y/nOCuiVYu6+IRFqpETGIDIhAZELLkfA74enpoYLbFsLC9BaVoaW2ho4nG4E5DpAoASfJ4OI3wYR3wmxoBkGiMAHH6jjgRXG8/P48AR4cILHXSFAgMct7D/87/qRL+CHx++F1++H1+eHl8dDQChBQOiHQOqFRMuDRGSDGCZIvDwY1Coo4vRQRMUAxmQcMIuwsbgJebUWpBj4WDg6hutxFg6gWcpI36KeMweOffsRHcWHP3cINpm+R6YxHYe313GzHhIykFA6QTcGcjpBu3X3v4xmtwjpj12DtIj/5a667UcmQSj8BRCrAJEM8HsBpxmInwQkzwRCEo7Zl8PqxubPCpE6zojYIV17R7vDgpGnVuYhI1KDhVk9myc84PEg5y+Pwy/IQNLtM6BI6FqN4JW9r8CoMOLK1CtP5TCQE3B6fDhc14aD1Wbk11lQZ3ZyubSpRjVSjUqkGFRdejLZefC1tMBrMsHX3AxPUxMcJhNcTY1wtbTC43bC4/Ig4PEj4HaD5/WwP1Bczz/LlxXwfBDxPBAL3BDLeBAphBCq5BAoxRAqxRCHqCEICQFPqgFPqgIkbFEDynBAaeDqFDfb3FyKxN6KVuTUmGHUSDEmIQTjE/XQKynflfQN7qpqNDzzDBzj5qJqfwUOhK9DiuJhjJ2fAv3/ptMmpLejdAJyzsjiZNAX8bF7Tz3SLvpfECuWAyMXH8mNtZkAjwMI+IDwdEDUfc6gz+tH9k/lXA3RngSwzMpDdfD6Arh0eGSP29vy8y9weiMROjauSwDLsML2bPT6rUNv7fH+yMlJRQJu4BdbGLPDg4L6Ni794Nt9NahpdUCvkCA+VIGEUDni9ArE6dWQh/xWqeJExYFYigILYFmJodPBguYWuwdFDVYUFDYhv64MtWYnksIUyIzWYtHoGC6IJaSvEUdHQX3RPPC2boMmIhMRDRloEf+M/WulmLJoEKUVkAGD0glIt2KnjUVN8R7U7CiDZ3YSRJ1TAFi5rE4ls04URBz6tRrs6m/6lJ4FpJXNdqw8UIuH5w7uGEB0Mt6mJpT9vAvS8EzEXjL6mDJMXxV9hUWpi7odBETOHo1MxNWcZUv7ALHyJhvKGu0oNtmwNq+B6wllKQhGtZQbLMV+GrhFArVMBImQ35EGwhP0LO2DpRqwEmyNVhfqzC7UWZzc+6ii2Y42pwfROjk3OHDBiEikRaghF9OfPdL3qWbNgmPvPkSFeeBpTsT25g3QiMqQu0XJVX8hZCCgv+akWxGZ6SjHrwix1mL9z6WYfUnyKR0pFsAe3lbL5cJOuDwJgh4EpCyNYNnmUszJMHC9dz1V9/5yCOQpMMwfDr6k61v6m6JvEKeKw/Cw/jcrV18YIJYeqeGWdizYZD20LPWgweJCbq0F6/Ib0Njm4qpRCPg8KKVCLg2BPV4q4kPAUgnYYK4A4PUHuO3YgD+H24c2p5frAWbvN41chAjNkaB4WIwGlwyL4AJY1mNMSH/DvuSF3HwT6v/1NKJnzcfoXydik/ldIPdPMMSzGfQGZiocGVgoiCXd4vMFcMkciPTYULDfhOGZ4adUwoVVIqjKb8GEy5K5EeY98ePBWi4YuSSz52kE9gMH0Fjmg2iwFmHDu84sxmZtyq7PxqNjHx1wkxr0Viw4VRtFx8y8xs47C0xZUMoW1ovLelFdHj8XuPoDAdahzwW5QgGPC0wVYiEUEgE3gYNWJqLBWGTAERkM0F75O5i//RYhGTMw/NAk5Ik+hHjd7zH1msGQyOgjnvRv9A4nxyWJ10GSa8UaVQC7finHjGsHc+WOToQFI0W7G1C634RxC5Kg0PZssExFkx0/HazDX+al9TgYYYN/Kj9aCchjkHjV5C73ef1eLD+8HHMT5kIv61kuLgke9iWDXeZni4E6kAjpMcWECVxlFmFrIWwhkQhpqkatcCf2r5Fj9MUJ9AWe9GtUR4YcV8rcCyDmKZGkaoBZxse2r4thNtmPu73b4cWuH0pRebgZY+cnQR0q69HRPZJGUIK5Q42I1fc8b7Xl+1/g9kVCOTUOspCus3+trVgLPo+P6THT6QwTQvr1F8CQa6+Fr64WMakyDMYw1JjXoaK0BKX7G4PdPELOKQpiyXGFxifAym9BSmUJDsi8iB6sw7avilGyzwSfj1XyPMJl96DsYCM2rSiAQMTH5CsHQRPWswCW+eFALfeH+KKhPZsEgfE0NKB2YxHcOh4SZ4/vcp/JbsIvZb/gqsFXcZMbEEJIf8aXy6G/9Ra416+CYbQRU60XYJ9jGQ5uKkVrw/E7Hgjp6+gTnhwXCyz9ESKoagNwev1wR0oxJioR+1ZXIG9rLeRqMcRSAcwmB3RGOQaNMyI6VXdKl6/KGm34+VAdHr3oFNIIAgHUv/ctvJJQRP9uPPidRrGz+z4r+AyjjaORqEmks0sIGRAkiYnQzL8UbWvWQJcyBZMLxyCbtwKqX5SYfGUqRBIa4Ej6HwpiyQnFzBoHy3uHMEXrxy859Xhg1iBMXzwYTpsH1hYXXDYP9NHKHg/e6oyNLn/j1xJclBmBmJCepxE4svfDUiuEL10KQ2rXqgl7Gvagqq0KN6XfRGeWEDKgKC+4AK6iYuhsBXDp45DQ0oi8wK9Qr5Mj68I4yo8l/Q6lE5ATSh4xBtaAGdGF+1BisnKzHLGeVha0hsWoED045LQCWNZj+v62MoQqxaeURuB3u1H72Ra0SbxIvWJml/vsHju+KPgCC1MWUk1YQsjAzI9dfB0CzU2IjPVjkHAovPZs5Bzaj9J9lB9L+h8KYskJsdmSXCFu+MosXLD5+e4qLgA9UxvyTdzsTrdNSQSf3/P0g5bPf4HDr4FmdioUWl2X+74v/h6RykhkGbLOuH2EENJX82ND77wD7p1bYMzSY7xrKsqdX2HHukNoqrYGu3mEnFUUxJKT0o4fDLFbhSlGGVe/c0dp8xkdtWKTFZ/trsTtU5K4uqE95a6uReNuE1r1NqRM6TqYq9Rcih11O3BV6lV0yYwQMqCJIiKgv/FGBNb/iPCsaMy0TMFhxwps/i4XjjZ3sJtHyFlDQSw5qWEzLoQ90IbKVetx2YgofLWnCp5O1QlOBZup6aU1hVg4MhqpRlWPH8d6fxveXQWz0IvkKy+AQCg6pibs7PjZCJOH0RklhAx4smHDuKlpxfvXIGRwNKaah2Bv09fY9n0BvB7fgD8+pH+gIJaclEQsg1VpgWV/FcYn6rmC9Cyt4FS12t14YXU+pgwKw8whhlN6rHXjbrQ184FMFYxJKV3u21i5Ef6AHzNju+bIEkLIQKa+aB4k8fHQmPYj1BiLURYVdhT9hL2ry85KWhghwUZBLOkR2eQUCB0yeGobcff0ZOwsbcLavPpTCmD/u7oAqUY1Fo6MOqWj7rM70PhjLuqlDci4dG6X+5ocTfix9EcsSl1ENWEJIeTogV7XXw++gIdweT2iNIOQbDPj113rcXh7LR0r0udREEt6ZNysy9Dkr0fxF6sQppLgDxck44vsKhysMvcoB/aJH3IRq1fghvGnXual+cO1aPW7YZw7GnKN9piasGwgV4qua+8sIYQQgCcWI/SOOxCoqUJkfABJkqHQO/OwavUGlOc00SEifVrQg9glS5YgISEBUqkUWVlZ2LRp0wm337hxI7cd2z4xMRGvv/56l/vfeustTJ48GTqdjltmzpyJnTt3nuNX0f9JxTK0RXvgLLXB7/MjOVyFGyfEY+nGIqzKqeOmjj0ay5tlvbXP/5KPuRkRuHlifI8nNGjnzC1Da34rmo1WpIzvOphrv2k/ysxlWJC84IxfHyGE9FcClQqhf/gDvPt2ISJdhqG8LEic2/HD1+tRX2YJdvMI6ZtB7IoVK3Dffffh0Ucfxd69e7ngc+7cuaioqOh2+9LSUsybN4/bjm3/l7/8Bffccw++/PLLjm02bNiAq6++GuvXr8e2bdsQGxuL2bNno7q6+jy+sv4p4fKZcHndaN2wl7s9NlHPpRb8WmjC377LwaZCE/ZWtOBQtRlf763Cnz7fj7WHG/B/F6Rg1hDDKffA+l0emJZvR42wDhkLL+kymMvhdeDzgs9xecrlUIgUZ/21EkJIfyIyGBD2h7vh274eEVl6jPaPhd+xCV8vX4fWepqalvRNvEAQs7vHjh2LkSNHYunSpR3r0tLSsGDBAjz99NPHbP/QQw/hu+++Q15eXse6O+64A/v37+cC1u74fD6uR/bVV1/F9ddf36N2WSwWaDQamM1mqNXq03pt/ZHP78P3DzyKBGUShv3ztk7rA9hY0IDtJc1ws+lpfX4Y1VLMSAvHkAj1aZe8anxvLWoOFMMxXoWxv7u6y30sgK211uL/RvwfldQihJAecubno/G1JcC0i1G3ow6bhXsgV0/HoltmQhPa85kTCTlTZyPWClpPrNvtRnZ2NtdL2hm7vXXr1m4fwwLVo7efM2cOdu/eDY/H0+1j7HY7d19ISMhZbP3AJOAL4BkVCr+FD3d9S6f1PFww2IC/zEvDPy5Nx78uG4p7ZqQgPVJz2gGmPacC5gO1qNaZMGzeJV3uK2ktwbaabdxgrtPdPyGEDETS1FTob7kZ2PADjKONmOQZAbvlV3z67hq0tTiD3TxC+kYQ29jYyPWSGgxdSy2x23V1dd0+hq3vbnuv18vtrzsPP/wwoqKiuNzY43G5XNw3gs4L6d7oeZfD5K1GzYoN5+wQ+RxuNK/YiQpBJTIWXgqpQtlxn8fnwcd5H+PC+AthUJxamS5CCCHgasiG3Hgj8OsPMIyOwCRPJlxNW/DxO7/A2kqBLOk7gj6w6+ieNJbdcKLete6272498+yzz2L58uX46quvuIFgx8NSF1iXdvsSExNzGq9kYIjTJaAkthGW4lZ46lvP+v7Z+Wx6dyMa2+ogHRWL6LSMLvf/XPYzxAIxZsTOOOvPTQghA4V85AiE3HgDeJt/hHFMBCb40+EzbcP7b/4ISzPlyJK+IWhBbGhoKAQCwTG9rg0NDcf0trYzGo3dbi8UCqHX67usf/755/Gvf/0Lq1atQmZm5gnb8sgjj3A5Ge1LZWXlab+u/o59WUi7cgFqPOWo/mj9Wd+/5ZeDaC2sRHV4C4ZfdGmXLyeVbZVYV7EO16Zdy6U2EEIIOX3ykSOhv+lGYNOPMI4Mw7jAEIhN2Vi29Ds0m9ro0JJeL2hBrFgs5kplrV69ust6dnvChAndPmb8+PHHbM+C1FGjRkEk+m3k+nPPPYcnn3wSP//8M3ffyUgkEi6puPNCjm9C3GTkJtbCUtEKZ5nprB0qZ2E9Wtflokhagqwrr4JMqeqSRvBBzgeYETcD0apoOj2EEHIWyIYPR+gdtwPbfoFxeAiyhEMQ2pyLd5Z8i5pKqiNLeregphM88MADePvtt/HOO+9wFQfuv/9+rrwWqzjQ3kPauaIAW19eXs49jm3PHrds2TL88Y9/7JJC8Nhjj3H3xcfHcz23bLFarUF5jf2RSCDC4EsvQaW3GLWfbjor0xd6muxo+ng7yvyFiJ099ZipZb8p/gYSgYTLhSWEEHL2SNPSEHbvPcDuDTCmipEhG4QESwk+eP0b5OWW0aEmvVZQg9hFixbhxRdfxBNPPIHhw4fj119/xcqVKxEXF8fdX1tb26VmLJsUgd3PasGy7Vlv68svv4yFCxd2mTyBVT644oorEBER0bGw9AJy9kxNmI4DybUw1zTBurXwjPblbXOh8c2NqLUcRiDDgCFTLuhyf05TDrbXbMcN6TfQ1LKEEHIOSBITEf7nP4FXfBDGCDuS1CkY4WjGNx+sxObN+85KZwUh/apObG9FdWJ75vPcFfC9uRlZ4gmIu28exJGaUz7WfqcXpjd+RX3pXlTFeTHtljsgV/+2nzZ3G/6141+Ynzwf4yLGnfL+CSGE9JyvrQ2NS5bCF+CjRRCHpto6bBe6EDMyE7+7YgaEQhqPQM6OPl0nlvR9FyTMxO4pPpTac1H71jr4HN3X6j0en8WFxne2oankAMrCrZi0+KYuAazX78Wyg8swSDcIY41jz8ErIIQQcvQUteH33wdpRBg0pt0wDorBNJ8Orbt24dVXP4fZQql5pPegIJacNr1Mj8mZc5CTZUeNqRim9zdzPas94a6xwvTGJjTmb0NBiAkTr78Z6tDwLtt8WfAlnD4nV42AJjUghJDzgycWQ3f99dBceCHkuethzAjBWGE0IstL8dpzy3Ewt4hOBekVKIglZ2Ruwly0DlKgOsGMqkPZqH91Pdy1x/+m7nd5Yd1ehcZlv6KmZiuKjHaMX3wT9NFda/Nurd6KPQ17cNvQ27i6sIQQQs4f1nGgumA6wu75PwgLdsFgtCNVE4NxVjtWLfsBy79aBa/XR6eEBBXlxHaDcmJPDRt49f7edzCrMBqiQ01I0o+GZmoGJEmhEBkVCPgC8Jld8NS0wbq5EO6aApT5imGLD8WEqxZDHda1BzanMQfLDi3DHcPu4FIJCCGEBI/PakPLhx/AWV4Je/RwmMusOORtQH24AYuunY2EuCg6PSQosRYFsefowA40LHfV7/djbG0EyleuQ4IvDJqQRAiVYWweLsBrh6etHvXuUtSofQgdkoFRF18OiVzeZT95TXl46+BbWDxkMUaEjwja6yGEEPIbNgbctmkTWr/8CoGU4Wis5cFka8FBIRA1MQNXXDQdQjEN+iI9R0HsOUJB7Kkzu8z4z+7/YKRhJEZ5U1C4aQOaS4qgCYgAPg9uER8u+BGamoa0SdMQGht/zD7ym/PxxoE3cPXgqzHaOPqsnEtCCCFnj7exEc0ffgR3QyNs4WmwVLlR6G1AuUaH2QunImtYKh1u0iMUxJ4jFMSennpbPf6757+YFTuLm1nLbjGjvriQJVdxVQfkWi1UIaHdfsPfXL0ZXxd9jStTr6RSWoQQ0tt7Zbdshfnrr+E3JqDZKofF4kJuoA32hCj87oqZiI7qmiZGyNEoiD1HKIg9fRWWCry892VuwNf0mOng8048dtDtc+PTw5/icPNh3JRxE1J0XWfqIoQQ0jv5rFaYv/oa9l274EkYjpYGoMlpRT7PC3FmIn43fwZ0Ib9NH05IZxTEniMUxJ6ZktYSvJ/7PlQiFa5OuxpRymOT/j0+D7bUbMHq8tUIlYXi5oyboZGc+mQJhBBCgstdWYnWL7+Eu6Iajog0mBt4aPC0oogXgDwzBZdfMg0hoRTMkq4oiD1HKIg9c6yH9eeyn7G+Yj3i1HHcEqGMQKuzFXW2OuS35HNB64XxF2JY2DCqA0sIIX08xcCZkwvL99/D3dQKW0gCrC0iNLrbUMrzwZcUhblzJyMxKYL+3hMOBbHnCAWxZ0+joxEFLQVcmkGtrRYh0hAYFAbEq+ORqkulP2aEENLfgtlDObCsXAl3TT3soUmwtElgc7tQGXCgUavCiGnjMG5MGqQKUbCbS4KIgthefGAJIYSQgRzMuktK0LZ6DRyHDsETloBmhxhetxw1vhY0CISQpSRh8rRRiEsJh0BIcy8NNBaqE9t7DywhhBBCAG9zM+zbt8O2bRvcNi9aZeFwulSwe92oC9jRKpZCNzQVE8ePRES8FkIR1ZsdCCwUxPbeA0sIIYSQrr2zroJC2LZthSN7L1zycDQHZPD7NLD7XKgP2NEmkkKRGIMxk8YgOimUUg76MQsFsb33wBJCCCGke36nE87cXDj2H4DjUA7sbglaRBoEoIPH70dTwIY2CODWyZA0MhODM1Khj1JBJKFe2v7CQkFs7z2whBBCCDm5gN/P5c86DhzkAlpHTStaBAo4hVoI+VrY/A60wgU7hPCFSBCbmYq0oZnQRakgkQnpEPdRFMT24gNLCCGEkFPnt9ngKimBq7AIjoJiNJc2oC0ghlesg0gQAg98MAfssPEAt4wHeaQag4ZnIC59CJQhUqp600dQENuLDywhhBBCzlzA7Ya7qhqe6iq4KirRnF+J5pomeHhK8MQhEAlVcAe8sPvdsPG8cCt8EEfIET8sGSkZw6DQh1Bg2wtRENuLDywhhBBCzl0KgtdkgrehAa6qWlTlFaO5vBY+hxBCvgYSkRoBngAOvxuOgAtuvgN+hReKKAUSx6QhLnMURDItnZ4goiC2Fx9YQgghhJx/fpcLrro61B8oRFVeIex1ZsApgBhySAQKCHkiuP1uuHw2uGGDV+iAQOWHNlqJqPRUhKeMhCTUAJ5YTKfvHKIgthcfWEIIIYT0Lj6fD+V55Sjbthe2ilqgzQOxVwgxJJDyJRBCAJ/fDY/XBq+/DT6eDTyxE1KtGNqocOgTk6GOTIAoJAQCrRY8mYxSFYIYawV9WN+SJUvw3HPPoba2Funp6XjxxRcxefLk426/ceNGPPDAA8jJyUFkZCT+/Oc/44477uiyzZdffom//vWvKC4uRlJSEp566ilcdtll5+HVEEIIIaS3EggESMxI5Jaja9hazU4U5Vai4dBhOKp94LfxIPZoIAkI4W8VwmoWwZbXAp+vFl6fDQGPHYGAFRC4IJALINWpoIkyQBUdDbUxCkKtFgKdDnyVCjw+zUh2LgQ1iF2xYgXuu+8+LpCdOHEi3njjDcydOxe5ubmIjY09ZvvS0lLMmzcPt912Gz766CNs2bIFd911F8LCwrBw4UJum23btmHRokV48sknucD166+/xpVXXonNmzdj7NixQXiVhBBCCOnNeDweVFoZRkwYxC2dedw+tDU7UF/VjJqCCrTV1MDdwgfPIYXIo4M4AIi9PAQahfA1i9FysBE+fw38Pgf8XicC7CfPBQg94EsEECjEkIaooQoPgzLCCGVYOKS6EPDVaggUCkpjOAW8APv6ESQsqBw5ciSWLl3asS4tLQ0LFizA008/fcz2Dz30EL777jvk5eV1rGO9sPv37+eCV4YFsKyL+qeffurY5sILL4ROp8Py5ct71C5KJyCEEELIyfh9fjhtXjja3DC3ONBY04DWylrYmprhsVoBhxsCtx8CHyAKACIeWwRHFoiOpCL4ffAFvPD7PPAH3PD7PfDzvPDDgwDfB78gAB7rcpTwIJCKIVRKIdGoIFGpIddoIddqoFBrIFGpIFLKIZDLwZNIen3vb59OJ3C73cjOzsbDDz/cZf3s2bOxdevWbh/DAlV2f2dz5szBsmXL4PF4IBKJuG3uv//+Y7ZhaQrH42JJ4C5XlwNLCCGEEHIifAEfcrWYW/RRSiRmhAFIP2Y7n8cPl8MLt8PL/bS3OdHW0gpLgwnWlma4LBZ4bXb4nSxAEoDvDYDv50EQ4EEYEEDg5UPo5YNv4wPNPHjAgw9OuHgmmNH4v2cJcP/5/X4E4IM/4Af778jaAAIB9m+2sC3b1x/ZArz/3eIFoBoVjck3XtMnTnzQgtjGxkYuwdpgMHRZz27X1dV1+xi2vrvtvV4vt7+IiIjjbnO8fTKs1/fxxx8/o9dDCCGEENIdgYgPuehIsHsE63kMB9A1deF42EVzvz8AvzcAn8cHr9MBl90Bl9MFq9UCa1sbnDYbPGyxsmDYDa/TB7/Lg4DbA7/PB3gDgI/tB+D5AixuBbsWz2O3A7z/LXwEAr4+cxKDPrCL60o/6kQdve5k2x+9/lT3+cgjj3CDxTr3xMbExJzCqyCEEEIIOTdYDCMQsAUQSQSAUgwlNAP+cActiA0NDeVGCR7dQ9rQ0HBMT2o7o9HY7fZCoRB6vf6E2xxvn4xEIuEWQgghhBDSNwQt61csFiMrKwurV6/usp7dnjBhQrePGT9+/DHbr1q1CqNGjeLyYU+0zfH2SQghhBBC+p6gphOwS/iLFy/mglAWfL755puoqKjoqPvKLvNXV1fjgw8+4G6z9a+++ir3OFZmiw3iYoO6OlcduPfeezFlyhQ888wzmD9/Pr799lusWbOGK7FFCCGEEEL6h6AGsawcVlNTE5544glusoOMjAysXLkScXFx3P1sHQtq2yUkJHD3s+oDr732GjfZwcsvv9xRI5ZhPa6ffvopHnvsMW7CAzbZAatHSzViCSGEEEL6j6DWie2tqE4sIYQQQkjvjrV6dyVcQgghhBBCukFBLCGEEEII6XMoiCWEEEIIIX0OBbGEEEIIIaTPCfqMXb1R+1g3lnRMCCGEEELOrvYY60zqC1AQ2422tjbuJ009SwghhBBybmMuVqXgdPx/e/cBGsXahXH82GJsN3YTe++9xW4Ue8UaUTRiQbGgRhErUbFgJWIvWLEEK4oNQY29xYKKFVvsxG4EW+bjvLD7JZpiouu47v8He5Od2d3MnR1nnz1z5h2G2IpHTEyMPHnyRLJkyWKuV4y435w03EdGRqZ4SAz8W9gmwPYA9g9I7meGXgdAM5aO+Z86dcq6W6nExkNXZv78+VO0Qj2FBlhCLNgmwD4CfGYgJbT6+qs5ghO7AAAA4HYIsQAAAHA7hFgkS/r06SUkJMT8BNgmwD4CfGbArhzBiV0AAABwO1RiAQAA4HYIsQAAAHA7hFgAAAC4HUIsUqxw4cJmoOLYtzFjxrBGPcjixYulSJEi4u3tLdWqVZNjx47ZvUiwyaRJk37YH/j6+vJ+eIijR49K27ZtzcD1+t7v3Lkzzny9tKhuIzo/Q4YMEhAQINeuXbNteWH/NtG7d+8f9hm1atVK1t8gxOKXTJkyRZ4+feq8TZgwgTXqIcLCwmT48OEyfvx4uXjxotSvX19atmxprsICz1SuXLk4+4MrV67YvUj4Q6Kjo6VSpUqycOHCeOfPmjVL5s2bZ+afO3fOfMFp2rSp8zLv8LxtQrVo0SLOPmPv3r2SHFyxC79EL81LtcUz6QdS3759pV+/fuZ+aGioHDhwQJYsWSIzZsywe/Fgg7Rp07I/8FD6BVZv8dEqrO4f9Atvx44dzbS1a9dKnjx5ZOPGjTJgwIA/vLSwe5tw0GG2fiVDUInFL5k5c6bkyJFDKleuLNOmTZPPnz+zRj2Avs8RERHSrFmzONP1/smTJ21bLtjr9u3b5tChtph069ZN7t69y1sCuXfvnjx79izO/kLDS8OGDdlfeLgjR45I7ty5pWTJktK/f3958eJFsp5PJRYpNmzYMKlatapky5ZNzp49K2PHjjU7q5UrV7JW/3FRUVHy7ds3U0mJTe/rhxU8j7+/v6xbt858GD1//lymTp0qderUMX2P+kUXnsuxT4hvf/HgwQOblgp20yptly5dpFChQiY7TJw4URo3bmwKJD97IQRCLOLQxvvJkycnula0n6l69eoyYsQI57SKFSuaMNu5c2dndRb/Pm3E//6w4ffT4BliHzasUKGC1K5dW4oVK2YOGwcHB9u6bPg7sL9AbIGBgc7fy5cvb3KFBto9e/Y4206SQohFHEOGDDGHAZMalSA+jrMK79y5Q4j9x+XMmVPSpEnzQ9VVDwV9X22BZ8qUKZMJs9piAM/m6HnU/YWfn59zOvsLxKbbhobY5OwzCLH4IZzoLSX0DHXHhoh/m5eXlxlS6+DBg9KhQwfndL3fvn17W5cNf4dPnz7J9evXzagV8GzaI61BVvcPVapUcfbVh4eHmyN3gHr58qVERkYmK0MQYpEip06dktOnT0ujRo3Ex8fHtBhoe0G7du2kYMGCrFUPoIeIe/bsaQ4B6aHj5cuXm+G1Bg4caPeiwQajRo0yY0Lqv3+tsGlP7Lt37yQoKIj3wwN8+PDBHIVz0B7HS5cuSfbs2c02ocPxTZ8+XUqUKGFu+nvGjBmle/futi437Nkm9Kbti506dTKh9f79+zJu3DhTRItdGEmSBaRARESE5e/vb/n4+Fje3t5WqVKlrJCQECs6Opr16UEWLVpkFSpUyPLy8rKqVq1qhYeH271IsElgYKDl5+dnpUuXzsqbN6/VsWNH69q1a7wfHuLw4cOWRorvb0FBQWZ+TEyM+Yzw9fW10qdPbzVo0MC6cuWK3YsNm7aJjx8/Ws2aNbNy5cpl9hkFCxY00x8+fJisv5FK/+OaDA4AAAC4BuPEAgAAwO0QYgEAAOB2CLEAAABwO4RYAAAAuB1CLAAAANwOIRYAAABuhxALAAAAt0OIBQAAgNshxALAbxIQEGAurwkAcD1CLAD8xbZv3y5NmzaVXLlyyX///Se1a9eWAwcO/NRzt23bZoK1j4+PZM6cWSpWrChTpkyRV69emflr1qyRVKlSOW96DfOuXbuaa5w7FC5c2MzbvHnzD69frlw5M09fBwD+NEIsAPzFjh49akLs3r17JSIiQho1aiRt27aVixcvJvq88ePHS2BgoNSoUUP27dsnV69elblz58rly5dl/fr1zsdpMH769Kk8efJENm7cKJcuXZJ27drJt2/fnI8pUKCArF69Os7rnz59Wp49eyaZMmVywf81ACSNEAsALvL69Wvp1auXZMuWTTJmzCgtW7aU27dvx3nMihUrTEjU+R06dJB58+ZJ1qxZnfNDQ0Nl9OjRJoyWKFFCpk+fbn7u3r07wb979uxZ8zgNrbNnz5Y6deqYiqqGYa3OBgUFOR+rlVRfX19ThdWAHBISYgLvnTt3nI/p0aOHhIeHS2RkpHPaqlWrzPS0adP+xjUGAD+PEAsALtK7d285f/687Nq1S06dOiWWZUmrVq3ky5cvZv6JEydk4MCBMmzYMFMB1ZA5bdq0RF8zJiZG3r9/L9mzZ0/wMRs2bDDtA4MGDYp3fuyQ/L0MGTKYn45lVHny5JHmzZvL2rVrzf2PHz9KWFiY9OnTJ4k1AACuQ4gFABfQiquG15UrV0r9+vWlUqVKJlw+fvxYdu7caR6zYMECU50dNWqUlCxZ0oROvZ8Yra5GR0eb3tXE/nbRokUlXbp0yVrmR48emcpt/vz5zfLEpoFVe181iG/dulWKFSsmlStXTtbrA8DvRIgFABe4fv26OdTu7+/vnJYjRw4pVaqUmadu3rwpNWvWjPO87+/HtmnTJpk0aZKpgubOnTvBx2nQ1DaBn/H27VtTtdXeVm1r+Pz5szmZzMvLK87jWrduLR8+fDA9utpKQBUWgN1oZgIAF9AgmVTAjC9sJvQ8Da59+/aVLVu2SJMmTRL921pFPX78uGkJSKoamyVLFrlw4YKkTp3atA0kdKKWBvKePXuantkzZ87Ijh07En1dAHA1KrEA4AJly5aVr1+/msDn8PLlS7l165aUKVPG3C9durQ5CSs27aGNrwKr/bU6eoBWRJPSvXt3UzVdvHhxvPPfvHnj/F3Da/HixU37QVIjDWj1VU/wat++vTlZDQDsRCUWAFxARxDQsNe/f39ZtmyZqXiOGTNG8uXLZ6aroUOHSoMGDcyIBDps1qFDh8xwWLGrsxpgdYSD+fPnS61atcywVo4TsHT81/hoC4OOaDBy5EjTg6ujHuTNm9eMOLB06VKpV6+eOZksuTR8R0VFmZEUAMBuVGIBwEV0bNVq1apJmzZtzEUKtFVAx3t1HOKvW7euCZUaYvXEr/3798uIESPE29vb+RoagLWiO3jwYDMMluOWVAidOXOmqdxqJVhHFtALEwQHB5sLHsQeYiu5tK/XMYIBANgplZVQAxYA4I/Tyu2NGzfk2LFjrH0ASATtBABgozlz5pjxYbUfVVsJdCzWhHpZAQD/RyUWAGyk470eOXLEXMBAT67SPlm9AAIAIHGEWAAAALgdTuwCAACA2yHEAgAAwO0QYgEAAOB2CLEAAABwO4RYAAAAuB1CLAAAANwOIRYAAABuhxALAAAAt0OIBQAAgLib/wFVtylbxkNFHwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 700x350 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "lc = _log_cpm(counts, lib)\n",
    "fig, ax = plt.subplots(figsize=(7, 3.5))\n",
    "plot_log_cpm_density(lc, ax, 'Per-sample log-CPM density (exons)')\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27f4024f",
   "metadata": {},
   "source": [
    "## 4. log2 transform\n",
    "\n",
    "For differential splicing, limma treats log-transformed exon counts\n",
    "as a microarray-style matrix. voom is not used here; the\n",
    "mean-variance trend is captured implicitly by the per-exon linear\n",
    "model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "179c268c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-04-20T19:13:14.968408Z",
     "iopub.status.busy": "2026-04-20T19:13:14.968046Z",
     "iopub.status.idle": "2026-04-20T19:13:14.980118Z",
     "shell.execute_reply": "2026-04-20T19:13:14.977355Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "log2_counts range: [0.00, 18.46]\n"
     ]
    }
   ],
   "source": [
    "log2_counts = np.log2(counts.values.astype(float) + 1)\n",
    "print(f\"log2_counts range: [{log2_counts.min():.2f}, {log2_counts.max():.2f}]\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0803730c",
   "metadata": {},
   "source": [
    "## 5. Design matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "44ba313a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-04-20T19:13:14.983915Z",
     "iopub.status.busy": "2026-04-20T19:13:14.983779Z",
     "iopub.status.idle": "2026-04-20T19:13:14.995315Z",
     "shell.execute_reply": "2026-04-20T19:13:14.993862Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "            treated  untreated\n",
      "untreated1      0.0        1.0\n",
      "untreated2      0.0        1.0\n",
      "untreated3      0.0        1.0\n",
      "untreated4      0.0        1.0\n",
      "treated1        1.0        0.0\n",
      "treated2        1.0        0.0\n",
      "treated3        1.0        0.0\n",
      "\n",
      "Contrast vector: [-1.  1.]  (treated - control)\n"
     ]
    }
   ],
   "source": [
    "design, C = gd.build_two_group_design(targets['condition'])\n",
    "design_df = pd.DataFrame(design, index=counts.columns,\n",
    "                         columns=sorted(targets['condition'].unique()))\n",
    "print(design_df)\n",
    "print(f\"\\nContrast vector: {C.ravel()}  (treated - control)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f90a16ef",
   "metadata": {},
   "source": [
    "## 6. Fit per-exon linear models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "69e3b470",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-04-20T19:13:15.002299Z",
     "iopub.status.busy": "2026-04-20T19:13:14.999874Z",
     "iopub.status.idle": "2026-04-20T19:13:15.327888Z",
     "shell.execute_reply": "2026-04-20T19:13:15.324008Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/John/miniconda3/lib/python3.11/site-packages/h5py/__init__.py:36: UserWarning: h5py is running against HDF5 1.14.3 when it was built against 1.14.2, this may cause problems\n",
      "  _warn((\"h5py is running against HDF5 {0} when it was built against {1}, \"\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "coefficients shape: (14599, 2)\n"
     ]
    }
   ],
   "source": [
    "fit = pylimma.lm_fit(log2_counts, design)\n",
    "print(f\"coefficients shape: {fit['coefficients'].shape}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23f7d460",
   "metadata": {},
   "source": [
    "## 7. `e_bayes` + `diff_splice` + `top_splice`\n",
    "\n",
    "The shipped exon matrix has no explicit gene annotation, so we\n",
    "simulate a simple geneid grouping: five consecutive exons per\n",
    "\"gene\". Real analyses would use a GTF-derived grouping."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "884ef35d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-04-20T19:13:15.344874Z",
     "iopub.status.busy": "2026-04-20T19:13:15.344364Z",
     "iopub.status.idle": "2026-04-20T19:13:15.488202Z",
     "shell.execute_reply": "2026-04-20T19:13:15.487152Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "top_splice output: (2920, 5), columns: ['GeneID', 'NExons', 'F', 'P.Value', 'FDR']\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/John/Documents/Projects/staged/pylimma/pylimma/pylimma/squeeze_var.py:396: UserWarning: Zero sample variances detected, have been offset away from zero\n",
      "  warnings.warn(\"Zero sample variances detected, have been offset away from zero\")\n",
      "/var/folders/0x/4309q_jn5xbf3zzq_4hcp2100000gn/T/ipykernel_99389/895928089.py:4: UserWarning: diff_splice: 14599 exons, 2920 genes, 0 with 1 exon, mean 5 exons/gene, max 5\n",
      "  ds = pylimma.diff_splice(fit, geneid=geneid)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>GeneID</th>\n",
       "      <th>NExons</th>\n",
       "      <th>F</th>\n",
       "      <th>P.Value</th>\n",
       "      <th>FDR</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>85</td>\n",
       "      <td>5</td>\n",
       "      <td>2583.673589</td>\n",
       "      <td>1.883358e-34</td>\n",
       "      <td>5.499406e-31</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1195</td>\n",
       "      <td>5</td>\n",
       "      <td>2013.896455</td>\n",
       "      <td>5.383258e-33</td>\n",
       "      <td>7.859557e-30</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2847</td>\n",
       "      <td>5</td>\n",
       "      <td>1940.547358</td>\n",
       "      <td>8.867099e-33</td>\n",
       "      <td>8.630643e-30</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2338</td>\n",
       "      <td>5</td>\n",
       "      <td>1809.610233</td>\n",
       "      <td>2.268798e-32</td>\n",
       "      <td>1.656223e-29</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2158</td>\n",
       "      <td>5</td>\n",
       "      <td>1370.151358</td>\n",
       "      <td>9.536298e-31</td>\n",
       "      <td>5.569198e-28</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>301</td>\n",
       "      <td>5</td>\n",
       "      <td>1339.796840</td>\n",
       "      <td>1.288355e-30</td>\n",
       "      <td>6.269993e-28</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>2462</td>\n",
       "      <td>5</td>\n",
       "      <td>1265.720229</td>\n",
       "      <td>2.764835e-30</td>\n",
       "      <td>1.153331e-27</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>883</td>\n",
       "      <td>5</td>\n",
       "      <td>992.929347</td>\n",
       "      <td>7.173435e-29</td>\n",
       "      <td>2.618304e-26</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>2794</td>\n",
       "      <td>5</td>\n",
       "      <td>954.707733</td>\n",
       "      <td>1.213902e-28</td>\n",
       "      <td>3.938439e-26</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>178</td>\n",
       "      <td>5</td>\n",
       "      <td>944.913253</td>\n",
       "      <td>1.393761e-28</td>\n",
       "      <td>3.976064e-26</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  GeneID  NExons            F       P.Value           FDR\n",
       "0     85       5  2583.673589  1.883358e-34  5.499406e-31\n",
       "1   1195       5  2013.896455  5.383258e-33  7.859557e-30\n",
       "2   2847       5  1940.547358  8.867099e-33  8.630643e-30\n",
       "3   2338       5  1809.610233  2.268798e-32  1.656223e-29\n",
       "4   2158       5  1370.151358  9.536298e-31  5.569198e-28\n",
       "5    301       5  1339.796840  1.288355e-30  6.269993e-28\n",
       "6   2462       5  1265.720229  2.764835e-30  1.153331e-27\n",
       "7    883       5   992.929347  7.173435e-29  2.618304e-26\n",
       "8   2794       5   954.707733  1.213902e-28  3.938439e-26\n",
       "9    178       5   944.913253  1.393761e-28  3.976064e-26"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fit = pylimma.e_bayes(fit)\n",
    "geneid = (np.arange(counts.shape[0]) // 5).astype(str)\n",
    "\n",
    "ds = pylimma.diff_splice(fit, geneid=geneid)\n",
    "ts = pylimma.top_splice(ds, number=np.inf, sort_by='p')\n",
    "print(f\"top_splice output: {ts.shape}, columns: {list(ts.columns)}\")\n",
    "ts.head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88c6729a",
   "metadata": {},
   "source": [
    "## 8. Top-ranked gene's exon-level fold changes\n",
    "\n",
    "A gene with differential splicing shows some exons strongly\n",
    "up-regulated while others are flat or down - a relative shift\n",
    "rather than a uniform one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "eeec7464",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-04-20T19:13:15.492290Z",
     "iopub.status.busy": "2026-04-20T19:13:15.492123Z",
     "iopub.status.idle": "2026-04-20T19:13:15.541710Z",
     "shell.execute_reply": "2026-04-20T19:13:15.541418Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Top gene: 85  (exons in gene: 5)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk4AAAFUCAYAAADS/LOVAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAQBRJREFUeJzt3QecTNf7+PGHZVfvJQhLiF4jiB6iBNGFIHqC6CWJEL2HEBIRfEVLI3qNGjVCdImIllUiWInoLNb8Xs95/Wf+M9vMrJndnZnP+/W6dufOnXvP3DvjPnvOc85JZLFYLAIAAIAnSvzkTQAAAEDgBAAA4AJqnAAAAJxE4AQAAOAkAicAAAAnETgBAAA4icAJAADASQROAAAATiJwAgAAcBKBE3zS0aNHpVOnTpI3b15Jnjy5WZ5//nnp0qWL7N+/XxK606dPS5s2bSRXrlym7Po++vXrJ//++6/DdsOHD5dEiRJFWpIlSyb+KHfu3NK+fXtJCF5++WWzPMmDBw+ka9euki1bNgkICJCSJUt65Dhnz541n4158+bJ09q2bZvZ15IlS556X4C3SRLfBQDcbebMmdKjRw8pUKCA9O7dW4oUKWL+kz9+/Lh89913UqZMGROYaDCSEF29elVeeuklSZMmjYwaNcoET4cOHZJhw4bJ1q1b5cCBA5I4sePfPOvXr5e0adPaHkd8HgnXF198YT6zn332mZQuXVpSpUoV30UCEAMCJ/iUn376Sbp16yb16tUzfw0HBgbanqtevbp0795dFi9ebGpxEqqVK1eamqVFixbJK6+8YtZVq1ZNwsLCZNCgQXLkyBEpVaqUw2v0hpspU6Z4KjGexm+//WY+jxrsA0j4+LMUPmXs2LGmuUP/grcPmuy9/vrrkj17dod12nzXoEEDyZAhg2nm0sDk+++/d9hGmzi05kprfd555x0TqGTMmFGaNGkif//9d6TjaOBTvnx5SZkypalFqF27tqk5epKkSZOan/Y1SCpdunTmZ1w0w2nz0ejRo6VgwYISFBQkmTNnlg4dOpjaMKvx48ebmq3Vq1c7vFabylKkSCG//vqrbd2cOXOkRIkSpux6jhs3bmxqACO+Ts+T1gbWrVvX/J4zZ07p37+/CRpj6+bNm/Luu+9Knjx5zGciR44c0qdPH7lz545tG73elStXjvTa8PBws71eY1fOjbP08zR79my5d++erZnV2pR2//59GThwoEO5NfC/fv36E/ern8fmzZtL6tSpzeeoRYsWcvnyZafLdfHiRencubM5/3ps/b40a9ZMrly54rDdw4cP5cMPPzTPaw1pjRo15MSJEw7bbNq0SRo2bCjPPvusuf758uUzTeb//PNPlM3Ox44dk5YtW5pyZ82aVTp27Cg3btxw2FbPgTbF62dJPyf6h9Kff/5pXq/7sXfq1Clp1aqVZMmSxVyvQoUKyeeff+6wzePHj8011VpqDWL1u1a8eHGZOnWq0+cMfsQC+IhHjx5ZkidPbilfvrxLr/vxxx8tgYGBlsqVK1sWLVpkWb9+vaV9+/YW/XrMnTvXtp3+ruuee+45S8+ePS0bNmywzJ4925I+fXpLtWrVHPY5ZswYS6JEiSwdO3a0rFmzxrJs2TJTrpQpU1qOHTsWY3muX79uyZUrl6VKlSqW3377zXLr1i3L9u3bzbr69es7bDts2DBTpmeeecaSOHFiS5YsWSxt2rSxnDt3LtJ+27VrZ7YNCQmJ8fjh4eGWV1991ZR1xIgRlk2bNpn3mSNHDkvhwoUtd+/eNds9fvzYUrduXfP+z549a9bNmTPHHEO3txo7dqxZ17JlS8vatWstCxYsMOcwbdq0lpMnTzqUT69DoUKFLB9//LFl8+bNlqFDh5rzqOVwRnBwsNmP1Z07dywlS5a0ZMqUyTJ58mSzz6lTp5pjV69e3bwHpeu0jPblUevWrTPrV61a5dK5UVWrVjVLTH7++WdzDvVzq7/rEhoaaspVu3ZtS5IkSSxDhgyxbNy40ZwTPW6pUqUs9+/fj/Y4WgY9h/oeP/vsM/M57dWrl/n8RPxMR+Wvv/6yZMuWzeGc6fdCP8vHjx8322zdutXsK3fu3JbWrVub6/rdd9+ZYzz//PPmu2j1xRdfWMaNG2fOoX6O58+fbylRooSlQIEClgcPHkT6LOt6ve56bvX4QUFBlg4dOti202tQqVIlS7JkySzjx48350avhR5XX6/7sdLvmp6HYsWKmc+dbtu/f3/zXRk+fLhtOy1fQECAee2WLVvM/wFTpkxx2AawInCCz7h8+bL5j/ONN96I9Jz+R/7w4UPbYr1hqoIFC5qbka6399prr5kbiP5HbR84devWzWG7CRMmmPWXLl0yj8+fP29ueBpc2dMASAOc5s2bP/G9/P333ybQ0v1al9dff93hhqn0ZqBBmt7gNQDUG0mGDBksWbNmNTdAe3rj05uDNciJjt4A9XhLly51WL9v3z6zfvr06bZ1//zzj+XZZ5+1lC1b1nLw4EFLihQpLG+++abt+f/++88EBRoc2NNzpDfEVq1aRQrsvv/+e4dt9bV6M41N4KQ3RL1JatntLVmyxBxLz5v1fWjQNmjQIIft9FrpubR+Nlw5N84ETtb3rQGRPb1x6/70s2VPAxhdP2vWrGiPo4GKbrNy5UqH17799ttOBU76OUmaNKnl999/j3Yba+AU8brqtdP1GgBGRb93ei41sI9YRmvgFPE96/dNgyTrd1aDNN1O36c9vdYRAycNPvXzeePGDYdte/ToYfZ57do123ddA2zAGTTVwS9oDpA2gVmXSZMmmfXaLPTHH39I69atzeNHjx7ZFm0uunTpUqSmB23Ss6dV+urcuXPm54YNG8zr27Zt67A/baaoWrWq6ZEUk//++880bWgT0zfffCM7duyQ6dOny65du8yxdV9W2vNO857q1Klj8qAGDBggP/zwg2k2mjBhgsN+v/zyS/Pa4ODgGI+/Zs0a01RRv359h/Jrb69nnnnGofzaVKlNkgcPHpQKFSqYRPYZM2bYnv/5559NM1TEnm7aBKQ5Z1u2bHFYr00tetyI59d6bq3NZ/bl0maWmN5L0aJFTdntX6PNpnos63vR96HHnT9/vm1/eh0030yvY5IkSVw+N0/jxx9/ND8jnjdtZtam34jnzZ42JWsTXcTPqTZXOUM/P/pZ0iatJ3nSd0GFhoaaXoN6zfU86vfP+hmM2Fwb3T612VL3o7Zv325+alOkPW3es6ev0fOkzcLadBzxu63P79mzx2xbtmxZkzuo+ZH6/dXvHhAdksPhMzTnSPMT7P/Ttvr222/l7t27JhCy/4/ZmrOhOTC6RCViLobeZO1p3oTSAMF+n9p7LypP6vH20UcfyeHDh8370C7qSvNvNKdGgw0Nptq1axft6/UmkD9/fttNwVVafs0hiS5HLOL5KFeunOm5qDcezf3SG7uVdfgE6/uwp3kxmv9iT29wEXO49PzqTc5KE+atN0+l5yK6Lvb6XjQ4tuaNxfReNJdm6dKlpkwaWGkPTM2tsg9eXD03saXnTYMMzZ+yp8GeBmgRh6WI+FrNDYpIX+cMDbo1H8kZT/ouaBBaq1Ytk3M1ZMgQKVasmPl86HrtOWrdzpV9Ws+N5jfZi/iedTsNkrS3oi4xXS/NJdNyff311ybw1zzJKlWqmO/iiy++6NS5gP8gcILP0P/sNLDYuHGjCZDsb9aFCxe2jWVjz9oTTf/jtE8AtqcJo66w7lN79T2pdicqGjRpInDEYMMaiGkvrCfRZvjYDklgTXrXIQ6iorUZ9nSYBE0E11q9oUOHymuvvSbPPfecw01Qr0dEejONTU9ATfy/deuWQ3mfFExrcnp0z1tpsKTB3Ny5c83v+lODQutnJzbnJrb0GHrT1yDGPnjS66pJ3tEF5dbX/vLLL5HWO5scrsf766+/xB30s6oBtQa29sG+BrNPe26uXbvmEDxFfH/p06c3/ydorawm1UdFE++VBmI6TpouGhhv3rzZ1OTq5+DChQsmoAesCJzgUzQA0qYGbRrQwCW6mgb7oEgHxtT/3LVHnjvof7b6H/GZM2ekadOmLr9eb97axKA9mzSAsm/2Uk+qDdCaJu1J1KtXr1iUXkzgs3DhQtMkpoFDTLR2Zty4cTJ48GDTU02brLQHlw4LobUy2qtQAxf9S16bmaz0xqzNUdpTy1WuBLL6XvS66s3WepOMjvUmO2XKFNm5c6fpaalBWmzPzdPQWjVtatXz1rdvX9t6rRHT3oDWYSqios1s2iN01apVDrWrWuvqDG32/eqrr0wTtat/NESkNWT2tUZWEc+rK7S5W8+NNhFrDaeVXhd7GuzoudCerNrcF10tYUTaFKufS/3+6Wda/9iyD54BAif4lIoVK5quxj179pQXXnjBdKnWZiStfdFaD73xKO06bf+fuN4sNODRZhkNVvSvWc2/0NwdHffJ1dGrR44cabppaxfpV1991fz1q808WhOgTQIjRoyI9vX617E2x9WsWVM++OADkxuif7lrd2ltjrDmYynt4v/mm2+afBRt4tL9T5w40TTLvP/++w771e7bmsOjAV1MNWFvvPGGOb7mgegAotr0pwGoBjuaP6P5V5o3oudTj603Mq110nOsNzNt4tBjawCiNyFtotG/3jVXSPNQtAlF37+WV1/nSXrj02uuZdIARG+g2kx0/vx5UzOpQx3YB0DaXKfNM5oPpAGfBoGxOTdPS6+9fh41Z03zbfRzraPh6/nSoRM0wIuOnudPPvnE/BwzZoz5w2DdunUmd8cZ+tnVPz70nOl10+Y1rYXRWjatkdEmY2fptjrQrH6OtbZMa4h0+IqITbSu0O+Tng+9dnputKZT/6hYsGCBed6+plWHE6hUqZJp6tYgS7+bWlupNV5aDmsumeasaS6cNstpjZs2k+vnV78nev4AB06lkANe5vDhw6YLc548eUzvLe1Bky9fPkvbtm1Nd+OIjhw5YnpQaXd+7VGkvd+0u/qMGTNs21h71UXsoWXtYaQ/7a1YscIMU5AmTRpTBu3x1axZM9O9+0m0h1rjxo1NjyB9rXbff+utt0xvNHvag1Dfl/bK0nLrMbp27Wp65cV2OAKlPZ+0+7t2G9dzlypVKtP7sEuXLpZTp06ZXorak0t7nFl7E1pNnDjRHGf58uW2ddplv3jx4qbnmnYPb9iwYaRhGaLqXWbf2yo2verU7du3LYMHDzY986zH1+7pffv2NT0xI6pQoYI5nnazj825cUevOnXv3j3LgAEDzHvSa6s9PN955x3TU9FeVMfRHpVNmzY1ZUudOrX5fffu3U71qlMXLlwwvev0e6DHzp49u/l+XLlyxeEzv3jxYofX6Wcr4jG0d17NmjVNOXToCu0dqp/jiD3grNf56tWrDvu0fu/sP7faG06/3+nSpTM9OXX/e/bsMdvp0BIRy6TvRYeM0PeSOXNmc41Hjx5t22bSpElmnQ7BoJ8RHVahU6dOT+yBCv+USP9xDKUAAPAu2hSptbHaTKw9PAFPIXACAHgV7fGoOUjajKhNc5rXp03U2oxp3+MS8ARynAAAXkV7L2oyuOb9abK89kDV/ER9DHgaNU4AAABOYuRwAAAAJxE4AQAAOInACQAAwBuSw3XyUu0JceDAATOY3vLly6VRo0YO2+gghDoInPaU0IHrdDBDHRVXJxN1hr5Gp3bQZELrKLYAAABWOjKTDo6qMzc8abqqeA2ctDeEjnzcoUOHKKem0BGOddRXHfFYRxpOmzatCaQiTgIaEw2adORlAACAmOjchE+a1irB9KrT2qCINU46vYFOZ6DzJsXWjRs3zLQPejLsp9kAAABQOn2PVrLo9EJaSeOV4zhpE9vatWvNnFc6Z5NO1KiTdOokrhGb8+yFhYWZxco6i7oGTQROAAAgOs6k9CTY5PDQ0FC5ffu2jB8/3kzqqBNy6uSZTZo0iXFkWJ2pXaNF60IzHQAAcJcE21SnuUk6S73Opq5zEFk1aNDAzC6vQ+47U+NkrX7TJjtqnAAAQEQaK2hlizOxQoJtqsuUKZMkSZJEChcu7LC+UKFCsmvXrmhfFxQUZBYAAAB3S7BNdYGBgVKmTBk5ceKEw/qTJ09KcHBwvJULAAD4r3itcdIcptOnT9seh4SEyOHDhyVDhgxmnKb33ntPWrRoIVWqVJFq1arJ+vXrZfXq1bJt27b4LDYAAPBTLuU46aaamL1z5045e/as3L17VzJnziylSpWSGjVquJyIrQGQBkQRtWvXTubNm2d+nzNnjkn4/uuvv6RAgQJmPKeGDRt6pN0SAAD4n5suxApOBU737t2TTz75RKZPny7//vuvGbRSE7eTJ08u165dk99++80kc9eqVUuGDh0qL730kiQUBE4AACBOk8Pz588v5cqVkxkzZpgxlXRQyojOnTtner9p09rgwYPl7bffdmbXAAAAXsOpGietUSpatKhTO3zw4IEJop5//nlJCKhxAgAA7ooVnOpV52zQZO0Nl1CCJgAAAHdyqqnu6NGjTu+wePHi4m9W580b30XwWvXPnInvIgAA4N7AqWTJkmZk7+ha9azP6c/w8HDnjw4AAOBrgZOOrwQAAODvnAqcGKkbAAAgliOHnzlzRqZMmSLHjx83zXM6f1zv3r0lL7k+AADAh7kcOG3YsEEaNGhg8p4qVqxocpt2794tRYoUMdOh1KxZ0zMlBZxAon7skagPAB4InD744APp27evjB8/PtL6AQMGEDgBAACf5dQ4Tva0ea5Tp06R1nfs2FF+//13d5ULAADA+wMnndT38OHDkdbruixZsrirXAAAAN7fVKdz0HXu3Fn+/PNPqVChgkkO37Vrl3z00UfSv39/z5QSAADAGwOnIUOGSOrUqWXSpEkycOBAsy579uwyfPhw6dWrlyfKCAAA4H2B06NHj+Sbb76Rli1bmgTxW7dumfUaSAEAAPg6l3KckiRJIu+8846EhYXZAiaCJgAA4C9cTg4vV66cHDp0yDOlAQAA8KUcp27dupkk8L/++ktKly4tKVOmdHi+ePHi7iwfAACA9wZOLVq0MD/tE8G1Z52OIK4/w8PD3VtCAAAAbw2cQkJCPFMSAAAAX8txOnfunOTIkUOCg4MdFl2nz7lix44dUr9+fTOcgdZWrVixItptu3TpYrbRyYUBAAC8InCqVq2aXLt2LdL6GzdumOdccefOHSlRooRMmzYtxu00oNq7d68JsAAAALymqc6ayxTRv//+GylR/Enq1KljlphcvHhRevToIRs2bJB69eq5WlwAAIC4D5yaNGlifmrQ1L59ewkKCrI9pwnhR48eNVOwuNPjx4+lTZs28t5770mRIkXcum8AAACPBU5p06a11TjpoJfJkye3PRcYGCgvvfSSmcfOnXT+Ox1005WpXHRwTusAnermzZtuLRMAAPBfTgdOc+fONT9z584t7777rsvNcq46cOCATJ06VQ4ePBhl02B0xo0bJyNGjPBo2QAAgH9yOTl82LBhHg+a1M6dOyU0NFRy5cplap100V57OvimBm/R0YmHNVHduly4cMHjZQUAAP7B5eTwK1eumBqnLVu2mMBGm+7suWsATM1tqlGjhsO62rVrm/UdOnSI9nWae2WffwUAABBvgZMmhp8/f16GDBki2bJlc6kZLaLbt2/L6dOnHQbXPHz4sGTIkMHUNGXMmNFh+6RJk8ozzzwjBQoUiPUxAQAA4ixw2rVrl2lGK1mypDyt/fv3O4z91K9fP/OzXbt2Mm/evKfePwAAQLwGTjlz5ozUPBdbL7/8skv7Onv2rFuOCwAAECfJ4TrlyQcffEAQAwAA/I7LNU4tWrSQu3fvSt68eSVFihQm78heVNOxAAAA+GXgxCS7AADAX7kcOGniNgAAgD9yOXCyjtW0YsUKOX78uBmOoHDhwtKgQQMJCAhwfwkBAAC8NXDScZfq1q0rFy9eNOMpaa+4kydPmt52a9euNblPAAAAvsjlXnU64a4GRzqVic4jd+jQITMgZp48eVyajBcAAMDna5y2b98ue/bsMaN7W+kI3+PHj5eKFSu6u3wAAADeW+Ok88DdunUryulTAgMD3VUuAAAA7w+cXnvtNencubPs3bvX5DfpojVQXbt2NQniAAAAvsrlwOnTTz81OU7ly5eXZMmSmUWb6PLlyydTp071TCkBAAC8MccpXbp0snLlStO7Tocj0BonHY5AAycAAABfFqtxnJQGSgRLAADAn7jcVNesWTPTgy6iiRMnyuuvv+6ucgEAAHh/4KTDEdSrVy/S+ldffVV27NjhrnIBAAB4f+AU3bADSZMmlZs3b7qrXAAAAN4fOBUtWlQWLVoUaf3ChQtNkjgAAICvcjk5fMiQIdK0aVM5c+aMVK9e3azbsmWLfPfdd7J48WJPlBEAAMA7Aycd5HLFihUyduxYWbJkiSRPnlyKFy8umzdvlqpVq3qmlAAAAN46HIEmh0eVIA4AAODLXM5xAgAA8FfxGjjp8AX169eX7NmzS6JEiUwToNXDhw9lwIABUqxYMUmZMqXZpm3btvL333/HZ5EBAIAfi9fA6c6dO1KiRAmZNm1apOfu3r0rBw8eNMno+nPZsmVy8uRJJhIGAADeN+WKO9SpU8csUUmbNq1s2rTJYd1nn30mZcuWlfPnz0uuXLniqJQAAABuqHH66aefJCwsTOLKjRs3TJOeTjQMAADgVYGT1hZdvHhR4sL9+/flgw8+kFatWkmaNGmi3U4DOR3B3H4BAACI98DJYrG4pRBPoonib7zxhjx+/FimT58e47bjxo0zzXzWJWfOnHFSRgAA4PsS/HAEGjQ1b95cQkJCTM5TTLVNauDAgaZJz7pcuHAhzsoKAAB821Mlh8+cOVOyZs0qng6aTp06JVu3bpWMGTM+8TVBQUFmAQAASFCBk+YbPY3bt2/L6dOnbY+1Vunw4cOSIUMGM25Ts2bNzFAEa9askfDwcLl8+bLZTp8PDAx8qmMDAAB41XAE+/fvl2rVqtke9+vXz/xs166dDB8+XFatWmUelyxZ0uF1Wvv08ssvx3FpAQCAv4vXwEmDn5gSzOMq+RwAAMAnksMBAAASCgInAAAAdzbVWXONnNGgQQOntwUAAPC5wKlRo0YOj3XaE/v8I31spb3fAAAA/LapTkfsti4bN240vdx++OEHuX79uhlkct26dfLCCy/I+vXrPV9iAAAAb+lV16dPH5kxY4ZUqlTJtq527dqSIkUK6dy5sxw/ftzdZQQAAPDO5PAzZ86YOeAi0nVnz551V7kAAAC8P3AqU6aMqXW6dOmSbZ2O6N2/f38pW7asu8sHAADgvYHTnDlzJDQ0VIKDgyVfvnxmyZUrlwmkvvzyS8+UEgAAwBtznDRQOnr0qGzatEn++OMP07uucOHCUqNGDYfedQAAAL4mVlOuaIBUq1YtqVKligQFBREwAQAAv+ByU50OSTBq1CjJkSOHpEqVSkJCQsz6IUOG0FQHAAB8msuB0+jRo2XevHkyYcIECQwMtK0vVqyYzJ49293lAwAA8N7AacGCBTJr1ixp3bq1BAQE2NYXL17c5DwBAAD4KpcDp4sXL5oE8aia8B4+fOiucgEAAHh/4FSkSBHZuXNnpPWLFy+WUqVKuatcAAAA3t+rbtiwYdKmTRtT86S1TMuWLZMTJ06YJrw1a9Z4ppQAAADeWONUv359WbRokZnYV4clGDp0qJmfbvXq1VKzZk3PlBIAAMBbx3HSSX11AQAA8Ccu1zg999xz8u+//0Zaf/36dfMcAACAr3I5cDp79qyEh4dHWh8WFmbyngAAAMTfm+pWrVpl+33Dhg2SNm1a22MNpLZs2SK5c+d26eA7duyQiRMnyoEDB8wkwcuXL5dGjRrZntd58EaMGGHGjfrvv/+kXLly8vnnn5uefQAAAAk2cLIGNJoQ3q5dO4fnkiZNaoKmSZMmuXTwO3fuSIkSJaRDhw7StGnTSM/r6OSTJ082I5Xnz5/fjFquCejaiy916tQuHQsAACDOAicdekDlyZNH9u3bJ5kyZXrqg9epU8csUdHapilTpsiHH34oTZo0Mevmz58vWbNmlW+//Va6dOny1McHAADwaI6TTurrjqDJmeNcvnxZatWqZVsXFBQkVatWld27d0f7Os21unnzpsMCAAAQb8MRaBPb9u3b5fz58/LgwQOH53r16uWWgmnQpLSGyZ4+PnfuXLSvGzdunMmLAgAAiPfA6dChQ1K3bl25e/euCaAyZMgg//zzj6RIkUKyZMnitsDJSnOqIjbhRVxnb+DAgdKvXz/bY61xypkzp1vLBAAA/JPLTXV9+/Y1o4dfu3ZNkidPLnv27DE1QKVLl5aPP/7YbQV75plnHGqerEJDQyPVQtnT5rw0adI4LAAAAPESOB0+fFj69+8vAQEBZtGcIq3R0R5wgwYNEnfRJHQNnjZt2mRbp82C2kRYoUIFtx0HAADAY011OvSAtalMa340z6lQoUJmXCf93RW3b9+W06dPOySEa2CmzX+5cuWSPn36yNixY+X55583i/6uTYKtWrVytdgAAABxHziVKlVK9u/fb8ZVqlatmpnkV3OcvvrqKylWrJhL+9L96D6srLlJOk6Ujt30/vvvy71796Rbt262ATA3btzIGE4AACBeJLJotrWLwc6tW7dMwHP16lUT5OzatUvy5csnc+fONQNaJiSaHK61YTdu3PBYvtPqvHk9sl9/UP/MGbfuj2uRcK4FAHgLV2IFl2ucXnzxRdvvmTNnlnXr1sWulAAAAL6eHK4ePXokmzdvlpkzZ5raJ/X333+bnCUAAABf5XKNkw498Oqrr5pEcO1Rp3PH6bxx2qvu/v37MmPGDM+UFAAAwNtqnHr37m2a6zRZW8dxsmrcuLFs2bLF3eUDAADw3honTQT/6aefJDAw0GF9cHCwXLx40Z1lAwAA8O4ap8ePH0t4eHik9X/99RfDBAAAAJ/mcuCkOU1TpkyxPdbBMDUpfNiwYWYOOwAAAF/lclPd5MmTpXr16lK4cGGTDK6jeJ86dUoyZcok3333nWdKCQAA4I2BU44cOcy0KAsXLpQDBw6YprtOnTpJ69atHZLFAQAA/DpwevjwoRQoUEDWrFkjHTp0MAsAAIC/SOzqBL86dpN1kl8AAAB/4nJyeM+ePeWjjz4yo4cDAAD4E5dznPbu3WsGuty4caMUK1ZMUqZM6fD8smXL3Fk+AAAA7w2c0qVLJ02bNvVMaQAAAHwpcJo7d65nSgIAAOBrOU46htP169cjrb9586Z5DgAAwFe5HDht27ZNHjx4EGm9Doa5c+dOd5ULAADAe5vqjh49avv9999/l8uXL9se69x169evN4NjAgAAiL8HTiVLljTjN+kSVZOcjhr+2Wefubt8AAAA3hc4hYSEiMVikeeee05++eUXyZw5s+25wMBAyZIliwQEBHiqnAAAAN4TOAUHB5ufOjddXNFBNocPHy7ffPONaRrMli2btG/fXgYPHiyJE7ucngUAAPBUnIo+fv75Z6d3eOfOHTl27Ji4g45QPmPGDJk2bZocP35cJkyYIBMnTqRJEAAAJNzAqW3btlKzZk35/vvv5fbt21FuownjgwYNknz58snBgwfdUjgN2Bo2bCj16tWT3LlzS7NmzaRWrVqyf/9+t+wfAADA7YGTBkUawAwdOlTSp08vRYoUMYFU/fr1pVKlSpIpUyYpXbq0nDt3TjZt2iRt2rQRd9B96/QuJ0+eNI+PHDkiu3btkrp167pl/wAAAG7PcUqaNKn06NHDLFqbpOM1nT17Vu7duyclSpSQvn37SrVq1SRDhgziTgMGDJAbN25IwYIFTeK5DnswZswYadmyZbSvCQsLM4v9wJwAAADxMuXKCy+8YJa4sGjRIvn666/l22+/NbVchw8flj59+kj27NmlXbt2Ub5m3LhxMmLEiDgpHwAA8C+JLDrGQAKVM2dO+eCDD6R79+62daNHjzbB1B9//OF0jZPuR2uu0qRJ45Fyrs6b1yP79Qf1z5xx6/64FgnnWgCAt9BYIW3atE7FCi7XOMWlu3fvRhp2QJvsYhoSISgoyCwAAADulqADJ00+15ymXLlymaa6Q4cOyeTJk6Vjx47xXTQAAOCHEnTgpFO4DBkyRLp16yahoaEmt6lLly6mdx8AAIBXBU7379+XZMmSiaekTp1apkyZYhYAAID45vK8JZpfNGrUKMmRI4ekSpVK/vzzT7Nea4a+/PJLT5QRAADAOwMn7dU2b948M/2JTu5rVaxYMZk9e7a7ywcAAOC9gdOCBQtk1qxZ0rp1a9PDzap48eLRDhEAAADgl4HTxYsXzXx0UTXhPXz40F3lAgAA8P7ASYcF0ClXIlq8eLGUKlXKXeUCAADw/l51w4YNM5P4as2T1jItW7ZMTpw4YZrw1qxZ45lSAgAAeGONkw5KqXPIrVu3ThIlSmTGVDp+/LisXr1aatas6ZlSAgAAeOs4TrVr1zYLAACAP3G5xgkAAMBfuVzjlD59etNEF5Gu01HEtcdd+/btpUOHDu4qIwAAgHcGTprTpBPv1qlTR8qWLSsWi0X27dsn69evl+7du0tISIi888478ujRI3n77bc9U2oAAABvCJx27dplRg/v2rWrw/qZM2fKxo0bZenSpWYwzE8//ZTACQAA+HeO04YNG6RGjRqR1r/yyivmOVW3bl3bHHYAAAB+GzhlyJDBDD0Qka7T59SdO3ckderU7ikhAACAtzbVDRkyxOQwbd261eQ4aVL4L7/8YsZ1mjFjhtlm06ZNUrVqVU+UFwAAwHsCJ034Lly4sEybNs2MGq7J4QULFpTt27dLhQoVzDb9+/f3RFkBAAC8bwDMihUrmgUAAMCfxCpw0jnqTp8+LaGhoeZ3e1WqVHFX2QAAALw7cNqzZ4+0atVKzp07Z5rp7Gm+U3h4uDvLBwAA4L2Bk47f9OKLL8ratWslW7ZsUY4iDgAA4ItcDpxOnTolS5YsMVOrAAAA+BOXx3EqV66cyW+KKxcvXpQ333xTMmbMKClSpJCSJUvKgQMH4uz4AAAAsa5x6tmzpxlu4PLly1KsWDFJmjSpw/M63Yq7/Pfff6b3XrVq1eSHH36QLFmyyJkzZyRdunRuOwYAAIDHAqemTZuanx07drSt0zwnTRR3d3L4Rx99JDlz5pS5c+fa1uXOndtt+wcAAPBo4BQSEiJxZdWqVVK7dm15/fXXzQCbOXLkkG7dujF5MAAA8I7AKTg4WOKKThT8xRdfSL9+/WTQoEFmapdevXpJUFCQtG3bNsrXhIWFmcXq5s2bcVZeAADg22I1AKb6/fff5fz58/LgwQOH9Q0aNBB30cE1deiDsWPHmselSpWSY8eOmWAqusBp3LhxMmLECLeVAQAAINaBk9YCNW7cWH799VdbbpOyjufkzhwnHSdK58WzV6hQIVm6dGm0rxk4cKCpobKvcdI8KQAAgDgfjqB3796SJ08euXLlihkeQGuAduzYYWqGtm3bJu6kPepOnDjhsO7kyZMxNhdqM16aNGkcFgAAgHgJnH7++WcZOXKkZM6cWRInTmyWSpUqmSYyzT9yp759+5opXrSpTseO+vbbb2XWrFnSvXt3tx4HAADAI4GTNsWlSpXK/J4pUyb5+++/ze9aCxSxduhplSlTRpYvXy7fffedFC1aVEaNGiVTpkyR1q1bu/U4AAAAHslx0gDm6NGj8txzz5lRxCdMmCCBgYGmJkjXudtrr71mFgAAAK8LnAYPHix37twxv48ePdoENZUrVzZToixatMgTZQQAAPDOwEkHpLTSGiYdluDatWuSPn16W886AAAAXxTrcZzsZciQwR27AQAA8K3ASZvpxo8fL1u2bJHQ0FAzSGXEcZ4AAAB8kcuB01tvvWXmjWvTpo0ZoJLmOQAA4C9cDpx++OEHWbt2rRmcEgAAwJ+4PI6TJoGT0wQAAPyRy4GTDkI5dOhQuXv3rmdKBAAA4M1NdaVKlXLIZdLpT7JmzSq5c+eWpEmTOmx78OBB95cSAADAWwKnRo0aeb4kAAAAvhA4DRs2zPMlAQAA8LUcp3379snevXsjrdd1+/fvd1e5AAAAvD9w6t69u1y4cCHS+osXL5rnAAAAfJXLgZPOTffCCy9EmUCuzwEAAPgqlwOnoKAguXLlSqT1ly5dkiRJ3DL1HQAAgG8ETjVr1pSBAwfKjRs3bOuuX78ugwYNMs8BAAD4KperiCZNmiRVqlSR4OBg0zynDh8+bMZ1+uqrrzxRRgAAAO8MnHLkyCFHjx6Vb775Ro4cOSLJkyeXDh06SMuWLSMNhgkAAOBLYpWUlDJlSuncubP7SwMAAOBLOU720qRJI3/++af7SgMAAOCrgZPFYnFfSQAAAHw5cIpr48aNM5MN9+nTJ76LAgAA/NBTBU5vvvmmaa6LCzrVy6xZs6R48eJxcjwAAAC3Bk5ffPGFZMqUSTzt9u3b0rp1a/nf//4n6dOn9/jxAAAAYt2r7tNPP3X67PXq1cvtZ1rnwKtXr57UqFFDRo8eHeO2YWFhZrG6efOm28sDAAD8k1OB0yeffOLw+OrVq3L37l1Jly6dbeTwFClSSJYsWdweOC1cuFAOHjxomuqczYMaMWKEW8sAAADgdFNdSEiIbRkzZoyULFlSjh8/LteuXTOL/q4T/44aNcqtZ/XChQvSu3dv+frrryVZsmROvcY6HYx10X0AAAC4QyKLi2MK5M2bV5YsWWKbbsXqwIED0qxZMxNcucuKFSukcePGEhAQYFsXHh5uetYlTpzYNMnZPxcVbapLmzatCaI8lci+Om9ej+zXH9Q/c8at++NaJJxrAQDewpVYweWRwy9duiQPHz6MtF4DmitXrog7vfLKK/Lrr786rNPpXQoWLCgDBgx4YtAEAADgTkliE8y8/fbb8uWXX0rp0qVN7c/+/fulS5cuJnnbnVKnTi1FixaNNN1LxowZI60HAABIcMMRzJkzx0z0W7ZsWZN3FBQUJOXKlZNs2bLJ7NmzPVNKAAAAb6xxypw5s6xbt05Onjwpf/zxh5l2pVChQpI/f36JC9u2bYuT4wAAADx14GSlgVJcBUsAAABeEzj169fP6R1Onjz5acoDAADg3YHToUOHnNqZJooDAAD4deC0detWz5cEAADAlyf5/euvv+TixYvuKw0AAIAvBU6PHz+WkSNHmhE2g4ODJVeuXGbOOp1uRZ8DAADwVS73qvvwww/N4Jfjx4+XihUrmuEIfvrpJxk+fLjcv3/fzGUHAADgi1wOnObPn28GumzQoIFtXYkSJcygmN26dSNwAgAAPsvlprpr166ZueIi0nX6HAAAgK9yOXDS2qVp06ZFWq/r9DkAAABf5XJT3YQJE6RevXqyefNmKV++vBm7affu3XLhwgUzFQsAAICvcrnGqWrVqmaeusaNG8v169dN81yTJk3kxIkTUrlyZc+UEgAAwFtqnDQwmjdvnqRJk0YWLFggLVq0IAkcAAD4HadqnNasWSN37twxv3fo0EFu3Ljh6XIBAAB4Z42T9pgbOHCgVKtWzYzb9P3335vap6i0bdvW3WUEAADwnsBpxowZ0q9fP1m7dq1JBh88eHCUE/rqOgInAADg14FThQoVZM+ePeb3xIkTm+TwLFmyeLpsAAAA3t2rLiQkRDJnzuyZ0gAAAPjSOE46sS8AAIA/crnGCQAAwF8ROAEAAPhC4DRu3DgpU6aMpE6d2iSjN2rUyIxQDgAAEB8SdOC0fft26d69u+nRt2nTJnn06JHUqlXLNhgnAABAgk0OP3LkiKxevVoyZMggzZs3l0yZMtmeu3nzpvTp00fmzJnjtsKtX7/e4fHcuXNNzdOBAwekSpUqbjsOAACAW2ucNm7cKGXLlpWFCxfKRx99JIUKFZKtW7fanr93757Mnz9fPMk61YsGbtEJCwszQZz9AgAAEKeB0/Dhw+Xdd9+V3377Tc6ePSvvv/++NGjQIFKtkKfoVC86enmlSpWkaNGiMeZFpU2b1rbkzJkzTsoHAAB8n9OB07Fjx6Rjx462qVXee+89mTVrljRr1sw033lajx495OjRo/Ldd9/FuJ3Oqac1U9blwoULHi8bAADwD07nOAUFBcn169cd1rVs2dJMwfLGG2/IpEmTxFN69uwpq1atkh07dsizzz77xHLqAgAAEG+BU8mSJU1OU+nSpR3Wt2jRQh4/fizt2rXzSPOcBk3Lly+Xbdu2SZ48edx+DAAAALcHTu+8846p8YmK1jwpbbpzJx2K4Ntvv5WVK1easZwuX75s1mvuUvLkyd16LAAAgCdJZNFqnQRKc6miosMStG/f3ql9aK86DbQ03ylNmjTiCavz5vXIfv1B/TNn3Lo/rkXCuRYA4C1ciRVcnuQ3LiXgmA4AAPghlwOn9OnTR1kTpOuSJUsm+fLlM7VBHTp0cFcZAQAAvDNwGjp0qIwZM0bq1KljBsTUWqF9+/aZ8Zw0JykkJMTkQ+n0KG+//bZnSg0AAOANgdOuXbtk9OjR0rVrV4f1M2fONKOLL126VIoXLy6ffvopgRMAAPDvSX43bNggNWrUiLT+lVdeMc+punXryp9//umeEgIAAHhr4KTzxEU1Urh18l91584dM3wAAACAXzfVDRkyxOQw6WCYmuOkSeG//PKLrFu3TmbMmGG22bRpk1StWtUT5QUAAPCewEkTvgsXLizTpk2TZcuWmeTwggULyvbt26VChQpmm/79+3uirAAAAPEqVuM4VaxY0SwAAAD+JFaBU3h4uKxYsUKOHz9umuq0BqpBgwYSEBDg/hICAAB4a+B0+vRp02vu4sWLUqBAAdNUd/LkScmZM6esXbtW8jL9CAAA8FEu96rr1auXCY4uXLggBw8elEOHDsn58+clT5485jkAAABf5XKNkyaB79mzxzb0gMqYMaOMHz+evCcAAODTXK5xCgoKklu3bkVaf/v2bQkMDHRXuQAAALw/cHrttdekc+fOsnfvXpPfpIvWQOkULJogDgAA4KtcbqrTOejatWsn5cuXl6RJk5p1OqGvBk1Tp071RBkBAE9hNZ12Yq3+mTN89vB0gVO6dOlk5cqVcurUKfnjjz9MjZMOR5AvXz5XdwUAAOD74zip559/3iwAAAD+wqnAqV+/fk7vcPLkyU9THgAAAO8OnHSsJmfoKOIAAAB+HTht3brV8yUBAADwteEI4sP06dPNyOTJkiWT0qVLy86dO+O7SAAAwA8l+MBp0aJF0qdPH/nwww9Nk2HlypWlTp06ZpoXAACAuJTgAydNNu/UqZO89dZbUqhQIZkyZYqZUPiLL76I76IBAAA/k6ADpwcPHsiBAwekVq1aDuv18e7du+OtXAAAwD/FehynuPDPP/9IeHi4ZM2a1WG9Pr58+XKUrwkLCzOL1Y0bN8zPmzdveqycdx8/9ti+fZ27rwvXIvY8+R1B/OJ7EXt8L/zrOlssFu8OnKIb5kDfWHRDH4wbN05GjBgRab027yEBSps2vksAK64FEBnfC79y69YtSfuEa56gA6dMmTJJQEBApNql0NDQSLVQVgMHDnQYsPPx48dy7do1yZgxo9+NM6URtAaMFy5ckDRp0sR3cfwa1yJh4XokHFyLhMOfr4XFYjFBU/bs2Z+4bYIOnAIDA83wA5s2bZLGjRvb1uvjhg0bRvmaoKAgs0ScX8+f6RfA374ECRXXImHheiQcXIuEw1+vRVonaxcTdOCktPaoTZs28uKLL0r58uVl1qxZZiiCrl27xnfRAACAn0nwgVOLFi3k33//lZEjR8qlS5ekaNGism7dOgkODo7vogEAAD+T4AMn1a1bN7PANdpkOWzYsEhNl4h7XIuEheuRcHAtEg6uhXMSWZzpewcAAICEPQAmAABAQkLgBAAA4CQCJwAAACcROPmw6dOnS548eSRZsmRmPKydO3fGd5H80o4dO6R+/fpmYDUdhHXFihXxXSS/pLMKlClTRlKnTi1ZsmSRRo0ayYkTJ+K7WH5JJ2kvXry4bbwgHWrmhx9+iO9i4f99T/T/qT59+nA+okHg5KMWLVpkPvgffvihHDp0SCpXrix16tQxY2Ahbt25c0dKlCgh06ZN49THo+3bt0v37t1lz549ZhDdR48emQnD9fogbj377LMyfvx42b9/v1mqV69uBjU+duwYlyIe7du3z4yVqEEtokevOh9Vrlw5eeGFF8xfdlaFChUyf2XrXxSIH/qX3PLly811QPy6evWqqXnSgKpKlSpcjniWIUMGmThxonTq1Cm+i+KXbt++be4Z2lIxevRoKVmypEyZMiW+i5UgUePkgx48eCAHDhwwf03b08e7d++Ot3IBCcmNGzdsN2zEn/DwcFm4cKGp+dMmO8QPrY2tV6+e1KhRg0vgCwNgwjX//POP+c8o4kTI+jjihMmAP9Lh63Q6p0qVKpnZCBD3fv31VxMo3b9/X1KlSmVqYgsXLsyliAcauB48eNA01eHJCJx8vFko4s0i4jrAH/Xo0UOOHj0qu3btiu+i+K0CBQrI4cOH5fr167J06VJp166daTYleIpbFy5ckN69e8vGjRtNRyI8GYGTD8qUKZMEBAREql0KDQ2NVAsF+JuePXvKqlWrTG9HTVJG/AgMDJR8+fKZ33USd63tmDp1qsycOZNLEoc0rUPvDdrz2kpbLPT7oR1awsLCzP0E/x85Tj76H5J+CbTnkD19XKFChXgrFxCftMZVa5qWLVsmP/74oxmqAwnr+uhNGnHrlVdeMc2mWvtnXTSQbd26tfmdoCkyapx8lOZvtGnTxnwBNI9Au5jqUARdu3aN76L5ZW+V06dP2x6HhISY/5A0KTlXrlzxWjZ/S3799ttvZeXKlWYsJ2uNbNq0aSV58uTxXTy/MmjQIDM8Ss6cOeXWrVsmx2bbtm2yfv36+C6a39HvQsQ8v5QpU0rGjBnJ/4sGgZOPatGihfz7778ycuRIuXTpkvkCrFu3ToKDg+O7aH5Hx6mpVq2aQ1CrNKdj3rx58Vgy/2IdmuPll192WD937lxp3759PJXKP125csX8Yaf/N2ngquMGadBUs2bN+C4a8ESM4wQAAOAkcpwAAACcROAEAADgJAInAAAAJxE4AQAAOInACQAAwEkETgAAAE4icAIAAHASgRMAAICTCJwAuNXw4cOlZMmST72f3Llzy5QpUySh0BHH+/TpE+M2OhJ8unTpYtxGRylv1KiRm0sHIK4QOAFeSm/AiRIlirS8+uqr8Vqud999V7Zs2SK+RicHHjVq1FMHdlOnTmWqHcCLMVcd4MU0SNK51uwFBQVJfEqVKpVZfI1OyuwOOjcbAO9FjRPgxTRIeuaZZxyW9OnTm+d0tvnAwEDZuXOnbftJkyZJpkyZzOSq6tdff5Xq1atL8uTJzWzonTt3ltu3b0dqVvr4448lW7ZsZpvu3bvLw4cPnW6qc2YfoaGhUr9+fVOOPHnyyDfffBNpvzdu3DDly5Ili6RJk8aU+8iRI+a5q1evmvc+duxY2/Z79+4173/jxo1RlrNp06bSs2dP22NthtMau2PHjpnHjx49MjPHb9iwIVJTnf5+7tw56du3r62mz56+plChQiaA1ODWer7tz4eV7qtXr17y/vvvm+BM34eew5ho2fQ12iyo53PAgAFm0mj7/VosFpkwYYI899xz5ryWKFFClixZYntePx9abq0dfPHFFyVFihRSoUIFOXHihMOxVq9eLaVLl5ZkyZKZfY0YMcIcH/BXBE6Aj7Le6HUWeg06NMj48MMP5X//+58JYO7evWtu6hpo7du3TxYvXiybN2+WHj16OOxn69atcubMGfNz/vz5pplJF1c8aR8aTJw9e1Z+/PFHc3OfPn26Cabsg4B69erJ5cuXZd26dXLgwAF54YUX5JVXXpFr165J5syZZc6cOSbg2L9/vwn+3nzzTenWrZvUqlUr2vOjwYPV9u3bTVCpP5Wek/v370vFihWjbLZ79tlnZeTIkSYosg+M9LxqkPjVV1/Jjh075Pz586b5MiZ6TlKmTGmCPQ12dL+bNm2KdvuPPvrIBJda2/jTTz/JzZs3ZcWKFQ7bDB482Dz/xRdfmGBQgzw9J9b3Z6WfCQ2o9bwlSZJEOnbs6BAA6ms0SPv9999l5syZ5rqNGTMmxvcD+DQLAK/Url07S0BAgCVlypQOy8iRI23bhIWFWUqVKmVp3ry5pUiRIpa33nrL9tysWbMs6dOnt9y+fdu2bu3atZbEiRNbLl++bDtGcHCw5dGjR7ZtXn/9dUuLFi2iLdewYcMsJUqUcChnTPs4ceKERf8r2rNnj+3548ePm3WffPKJebxlyxZLmjRpLPfv33c4Vt68eS0zZ860Pe7WrZslf/78ltatW1uKFi1quXfvXrTlPHr0qCVRokSWq1evWq5du2ZJmjSpZfTo0aZsauzYsZZy5crZtq9ataqld+/etsf6nqzls5o7d64p9+nTp23rPv/8c0vWrFkdzkfDhg0d9lupUiWH/ZQpU8YyYMCAaMuu+5s4caLtsZ7bXLly2far1zRZsmSW3bt3O7yuU6dOlpYtW5rft27dasq6efNmh+uv66znrXLlyuY82Pvqq68s2bJli7ZsgK8jxwnwYtWqVTM1CtHl4mhT1ddffy3FixeX4OBgh2Tm48ePm+Ybremw0tqVx48fm+aarFmzmnVFihSRgIAA2zZaW6VNfK6IaR9aDq3p0OYiq4IFCzr0TtMaJq1F0mYpe/fu3TM1WVZa01O0aFH5/vvvTQ2KNi9FR7fT/WkNTNKkSc25aNCggXz66afmea2Nqlq1qrhKm7zy5s3r8F7ta8+iotfHXkyv0drDK1euSNmyZW3r9Nxqc5peO6W1Q1pbVrNmTYfXPnjwQEqVKhXtsfW4So+dK1cuc9615s2+hik8PNzsW2vW9L0C/obACfBiGvTky5cvxm12795tfmqTli7WQEmbvyLm5ljZr9egIuJz1hu0s2Lah5Yj4jEj0m31pm7ftGZlH2D9+eef8vfff5vtNQcpYkASsQxVqlSx5YJp050GUxoYaFCn5+1Jww84+16t79GV1zzpHEc8X/bHsL527dq1kiNHjhg7D9gf27pP6+v1p+Y0NWnSJNLxYwpKAV9G4AT4MK2N0dwWzWvSWpi2bduaZODEiRNL4cKFTW7NnTt3bMGU5svoc/nz54+zMmoStSYbaw2RtRZFa7yuX79u20bzmTS/SWumdBiAqGhtSuvWraVFixamxqpTp04mALLWnEVFg6VZs2aZwEnzijRwqFy5sqm50tqsqPKbrPQ1GmTFNe2Vp+/pl19+MWVVWo5Dhw7ZkvL12mqApPlVsak1sz/vei2eFJwD/oTkcMCLhYWFmYDCfvnnn39sN1NNDNfk6A4dOphE4d9++80kAisNMrTWQHtj6XpN3NZeZvqamIINdytQoIBJUn/77bdNcrQ2D7311lumJ5hVjRo1pHz58qbXmCYsayK51ghpArQGXNYkZ23G0qY27aGmAZkGTzHRwEkTpzXAsgYhuk4TrzVo0N570dEATpO/L168aDvncUWv07hx42TlypUmsOndu7f8999/thoj7Q2oCekaNGtwrAG0Blaff/65eeysoUOHyoIFC0zSvZ4nbVZdtGiROe+AvyJwArzY+vXrTROW/VKpUiXznOalaIChNSpKu7nPnj3b3PQOHz5s8lM0CNHmuzJlykizZs1ML7Vp06bF+fvQoC5nzpymdkSbhazDDlhpQKC96bRpTXt9aY3YG2+8Yd6fBnna3Kb5W9qTTYMdrTXT33ft2hUpB8yeNs1pTzrNb7IGSVoGDTqfVFOjNVR6fM1n0l59cUmHH2jZsqWpQdSAUoc9qF27tkPzmQ7WqYGPBlgaROrzOrSADvfgLH3NmjVrTA8//Yy89NJLMnnyZJMvB/irRJohHt+FAADEnuYiaXDUvHlzh9HNAbgfOU4A4GU08V0H9tRaMW2u1VrCkJAQadWqVXwXDfB5NNUBgJfRpkgdiFKbzzSBXXO0dPBSrXUC4Fk01QEAADiJGicAAAAnETgBAAAQOAEAALgXNU4AAABOInACAABwEoETAACAkwicAAAAnETgBAAA4CQCJwAAAHHO/wE11ehhckb5UAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 600x350 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "top_gene = str(ts.iloc[0]['GeneID'])\n",
    "gene_mask = (geneid == top_gene)\n",
    "exon_fc = fit['coefficients'][gene_mask, 0]\n",
    "print(f\"Top gene: {top_gene}  (exons in gene: {gene_mask.sum()})\")\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(6, 3.5))\n",
    "ax.bar(range(len(exon_fc)), exon_fc,\n",
    "       color=np.where(exon_fc > 0, 'firebrick', 'steelblue'))\n",
    "ax.axhline(0, color='k', linewidth=0.5)\n",
    "ax.set_xlabel('Exon index within gene')\n",
    "ax.set_ylabel('log2 fold-change (treated - control)')\n",
    "ax.set_title(f'Gene {top_gene}: exon-level fold changes')\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  }
 ],
 "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.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}