{
"cells": [
{
"cell_type": "markdown",
"id": "0",
"metadata": {},
"source": [
"# Reed Frogs: A Discrete Likelihood Validated Against an Exact Binomial\n",
"\n",
"This notebook reproduces the multilevel tadpole survival analysis from\n",
"McElreath's *Statistical Rethinking* (Ch. 13), using real experimental data\n",
"from Vonesh & Bolker (2005). The setup: 48 tanks of reed frog tadpoles\n",
"(*Hyperolius spinigularis*) with varying initial densities, predator\n",
"treatments, and tadpole sizes. The outcome is the number of survivors.\n",
"\n",
"The survivor count is a **pure discrete** observation, so the estimator here\n",
"is `DiscreteNLE` -- a neural likelihood for discrete data (a\n",
"CategoricalMADE over the count, conditioned on the initial density). Because\n",
"the Binomial likelihood is tractable, we can fit the exact model in PyMC and\n",
"compare it against the learned one. If `DiscreteNLE` has learned the\n",
"likelihood well, both approaches recover matching posteriors. This validates\n",
"the neural likelihood bridge before moving to intractable simulators in the\n",
"next notebook.\n",
"\n",
"**Outline:**\n",
"1. Load and explore the data\n",
"2. Exact PyMC model: hierarchical Binomial with varying intercepts\n",
"3. Add predator effect and observe variance reduction\n",
"4. Train a `DiscreteNLE` on simulated count data and validate it\n",
"5. Replace `pm.Binomial` with the `DiscreteNLE` and compare posteriors"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "1",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-17T11:44:53.469833Z",
"iopub.status.busy": "2026-07-17T11:44:53.469650Z",
"iopub.status.idle": "2026-07-17T11:44:55.737589Z",
"shell.execute_reply": "2026-07-17T11:44:55.737134Z"
}
},
"outputs": [],
"source": [
"import arviz as az\n",
"import jax\n",
"import jax.numpy as jnp\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import pandas as pd\n",
"import pymc as pm\n",
"\n",
"import setu"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "2",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-17T11:44:55.738947Z",
"iopub.status.busy": "2026-07-17T11:44:55.738767Z",
"iopub.status.idle": "2026-07-17T11:44:55.741175Z",
"shell.execute_reply": "2026-07-17T11:44:55.740678Z"
}
},
"outputs": [],
"source": [
"# Shared setu-docs figure style: applied when this notebook is executed as a\n",
"# docs tutorial (setu_docs_style is on the build's PYTHONPATH); skipped\n",
"# harmlessly when the notebook is run standalone from examples/.\n",
"try:\n",
" from setu_docs_style import set_style\n",
"\n",
" set_style()\n",
"except ImportError:\n",
" pass"
]
},
{
"cell_type": "markdown",
"id": "3",
"metadata": {},
"source": [
"## The data\n",
"\n",
"The dataset comes from Vonesh & Bolker (2005), studying predation on\n",
"African reed frog tadpoles. The experiment used a 3 x 2 x 2 factorial\n",
"design:\n",
"\n",
"| Factor | Levels |\n",
"|--------|--------|\n",
"| Initial density | 10, 25, 35 tadpoles per tank |\n",
"| Predator | present / absent |\n",
"| Tadpole size | big / small |\n",
"\n",
"Each of the 12 treatment combinations has 4 replicate tanks, giving\n",
"48 tanks total. Each tank records the number of survivors."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "4",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-17T11:44:55.742187Z",
"iopub.status.busy": "2026-07-17T11:44:55.742114Z",
"iopub.status.idle": "2026-07-17T11:44:55.754026Z",
"shell.execute_reply": "2026-07-17T11:44:55.753629Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
density
\n",
"
pred
\n",
"
size
\n",
"
surv
\n",
"
propsurv
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
10
\n",
"
no
\n",
"
big
\n",
"
9
\n",
"
0.9
\n",
"
\n",
"
\n",
"
1
\n",
"
10
\n",
"
no
\n",
"
big
\n",
"
10
\n",
"
1.0
\n",
"
\n",
"
\n",
"
2
\n",
"
10
\n",
"
no
\n",
"
big
\n",
"
7
\n",
"
0.7
\n",
"
\n",
"
\n",
"
3
\n",
"
10
\n",
"
no
\n",
"
big
\n",
"
10
\n",
"
1.0
\n",
"
\n",
"
\n",
"
4
\n",
"
10
\n",
"
no
\n",
"
small
\n",
"
9
\n",
"
0.9
\n",
"
\n",
"
\n",
"
5
\n",
"
10
\n",
"
no
\n",
"
small
\n",
"
9
\n",
"
0.9
\n",
"
\n",
"
\n",
"
6
\n",
"
10
\n",
"
no
\n",
"
small
\n",
"
10
\n",
"
1.0
\n",
"
\n",
"
\n",
"
7
\n",
"
10
\n",
"
no
\n",
"
small
\n",
"
9
\n",
"
0.9
\n",
"
\n",
"
\n",
"
8
\n",
"
10
\n",
"
pred
\n",
"
big
\n",
"
4
\n",
"
0.4
\n",
"
\n",
"
\n",
"
9
\n",
"
10
\n",
"
pred
\n",
"
big
\n",
"
9
\n",
"
0.9
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" density pred size surv propsurv\n",
"0 10 no big 9 0.9\n",
"1 10 no big 10 1.0\n",
"2 10 no big 7 0.7\n",
"3 10 no big 10 1.0\n",
"4 10 no small 9 0.9\n",
"5 10 no small 9 0.9\n",
"6 10 no small 10 1.0\n",
"7 10 no small 9 0.9\n",
"8 10 pred big 4 0.4\n",
"9 10 pred big 9 0.9"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"frogs = pd.read_csv(\"data/reedfrogs.csv\")\n",
"n_tanks = len(frogs)\n",
"density = frogs[\"density\"].values\n",
"survivors = frogs[\"surv\"].values\n",
"predator = (frogs[\"pred\"] == \"pred\").astype(int).values\n",
"size = (frogs[\"size\"] == \"small\").astype(int).values\n",
"frogs.head(10)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "5",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-17T11:44:55.755106Z",
"iopub.status.busy": "2026-07-17T11:44:55.755036Z",
"iopub.status.idle": "2026-07-17T11:44:55.886376Z",
"shell.execute_reply": "2026-07-17T11:44:55.885881Z"
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA90AAAFUCAYAAAA57l+/AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAhyhJREFUeJzt3Qd4VFX6BvA3vZKEhBJ6770LgqigKAKWVRFdwa5rl9VVbIiuoqIoKiuWtf2FRcVVV1EsKAiCKE16lQ4JJYE00uf/vGeYYZJMQkJmMu39+YzM3Llz58y9kzn3u+c75wRZLBYLRERERERERMTlgl2/SRERERERERFR0C0iIiIiIiLiRmrpFhEREREREXETBd0iIiIiIiIibqKgW0RERERERMRNFHSLiIiIiIiIuImCbhERERERERE3UdAtIiIiIiIi4iYKukVERERERETcREG3iJ9q3rw5rrvuOrdtf8GCBQgKCjL/VuSJJ54w6x0+fNhtZREREQl0rPdHjBgBf3D22Webm7vs3LnTnJu89957bnsPEUcKusWj+GPHHz3eFi9eXOZ5i8WCJk2amOf9pSIR9/jXv/7l8cpzyZIl5iLD0aNH4QueeeYZfP75554uhohIjZxn8BYZGYm2bdvizjvvRGpqqtv3vH5nRYQUdItXYCU4a9asMssXLlyIvXv3IiIiwiPl8mWbN2/GW2+9hUDhLUH3pEmTFHSLiHiZJ598Ev/3f/+H1157DQMGDMDrr7+O/v37Iycnx63vq6D79Hz33XfmJuIvFHSLVxg+fDg++eQTFBYWlljOQLxXr15ITk5GoKvqiQEvVISFhbmtPL4sOzvb00UQEZEadOGFF+Kvf/0rbrrpJnOB9t5778WOHTvwxRdfVHvb7g7cA7HeDQ8PNzcRf6GgW7zCmDFjcOTIEXz//ff2Zfn5+ZgzZw6uvvpqp68pLi7Gyy+/jE6dOpmW8vr16+PWW29Fenp6ifVYoV500UVo2LChCURbtWqFp556CkVFRSXWY9+hzp07Y8OGDTjnnHMQHR2NRo0a4fnnn6/UZ2DZBw4ciISEBMTGxqJdu3Z4+OGHy6S4sR/RqfpG28qyYsUKnHXWWaYs3BZT7Fu2bOn0/XnFvnfv3k77dC9fvty8x/vvv1/mdd9++6157quvvjKPd+3ahdtvv92UPyoqCklJSbjiiivKlLuq2Kf7yiuvRFxcnNnmPffcg9zcXPvzgwcPRrdu3Zy+lmUZNmxYudvmZ12/fr3JjLClENr6gtn2O5/j56pXrx4aN25sf+0333yDQYMGISYmBrVq1TLfFW7L0Zo1a8y+5L7nd40XgW644QbznbVhWvkDDzxg7rdo0cJeDtt+432mM/LiUseOHc2+5TFbu3atef6NN95A69atzfZZdmf7e9myZbjgggsQHx9vvhPcZ7/88ovTPvTbtm0zZeb3ketff/31JU4MuQ5PgvidsJXVnWMAiIh4k3PPPdf8y8Db5sMPPzQX+vn7nJiYiKuuugp79uwp8bry6mdnKvqdrWxda6vD+Fs/fvx41K1b19RXl156KQ4dOnTKz8n3Dg0NtddP5eF5AuvZOnXqmPKwHmM9d6pxXJz1jeZn5HnQ9u3bTaMK69ZrrrnG1IFc7uwiBc8DWbfazs0c+3SzGwA/AzPJnGX18f2ZwUBpaWm4//770aVLF/NePOfgBZc//vjjlPtKxJ1C3bp1kUpi0MQA5D//+Y/5cbQFQ8eOHTOV3iuvvFLmNQyw+SPPYOLuu+82FSd/dFetWmUqJ1srL9fhDy8rK/77448/4vHHH0dGRgamTJlSYpsM2BnUXHbZZSZAZND/4IMPmh9vW7mcYZDGgLhr164mhY3BPYOe0gFRVTCg43vy8/PqPC8q8GRg7Nix+P3339GnTx/7uqy8f/311zKfx4bBOAPGjz/+GOPGjSvx3EcffYTatWvbg1pum2nSfF8Gp6xQmYbHyo8XJHiCcTq4P3mcJ0+ebMrKY8r9/cEHH5jnr732Wtx8881Yt26dOaGxYXm2bNmCRx99tNxt8+LLXXfdZY7vI488YpZxfzniyQ1PVnjsbVfcmWrI/cHP/txzz5kTAX5WXjzh94jltV1Q+fPPP813jScFPN5vvvmm+ZefhRU+vzMsJ7/DL730kjlxIb6nzaJFi/C///0Pd9xxh3nMfcHvzT/+8Q+THs8ycp/wQg9PdvhdteF9fh/4HZg4cSKCg4Px7rvvmhNHbrdv375l9jdPmvgeK1euxNtvv20uOPBz2j47W3z4ultuucUs4wUpEZFAwICQGOzS008/jccee8z8dvK3kQHtq6++agJr1ge8gFlR/exMRb+zVa1rWcexrubvP9dlvccglnV4eVhP3XbbbeaiwD//+c9y1zt48CDOP/98U1899NBD5rPyPf773//idDFzkXUr69MXXnjBfB7WqdOnT8fcuXPNBQYb1r1ffvmlCdZDQkLKbIv7lxeZeQ7Dz++In5+vsW2PdTXHKuFj1oEM2HlRm6/nfmUDjIhHWEQ86N1337Xwa/j7779bXnvtNUutWrUsOTk55rkrrrjCcs4555j7zZo1s1x00UX21y1atMi8bubMmSW2N2/evDLLbdtzdOutt1qio6Mtubm59mWDBw82r/3ggw/sy/Ly8izJycmWv/zlLxV+jpdeesm89tChQ6f8rDt27Cix/KeffjLL+W/pssyYMaPEuseOHbNERERY/v73v5dY/vzzz1uCgoIsu3btsi/jPhs3bpz98YQJEyxhYWGWtLS0Ep8vISHBcsMNN1S4v5YuXVpm3zgrtzMTJ040640aNarE8ttvv90s/+OPP8zjo0ePWiIjIy0PPvhgifXuvvtuS0xMjCUrK6vC9+nUqZPZb+Xt94EDB1oKCwvtyzMzM81nv/nmm0usn5KSYomPjy+x3Nk++c9//mO2+/PPP9uXTZkyxekxJi7nsXN87o033jDL+R3LyMgocawct1NcXGxp06aNZdiwYea+Y7latGhhOe+888rsb8djSpdeeqklKSmpxDLuV8fviIiIv7HVAT/88IOpo/fs2WOZPXu2+T2Mioqy7N2717Jz505LSEiI5emnny7x2rVr11pCQ0NLLC+vfi5Peb+zla1rbeUfOnRoid//++67z5SZdaeN47nStGnTzHnBU089dcoyfvbZZ/ZzsfKUV+eznuJyltOGn5fLHnrooRLrsvyNGjUqc0718ccfl6lPuZ8d63Rbfclj4qhjx46Wc8891/6Y53VFRUVlysj698knn6yw3CLupPRy8Rq8unz8+HGT5pyZmWn+LS+1nCm6TJk977zzTNqy7cZWQLZ2/vTTT/Z1mSZlw+1yPaYT88rqpk2bSmyXr+VVaxv2J+IVal45rYjtCjhT2Zn27gpsLWfLqiNbmhSv9lrjuJNXes844ww0bdq03O2NHj0aBQUFJa5cc5ASjrTN55ztL67PK/pMe+ZnZIvp6bK17jpetaevv/7a/MvjefHFF5uWYttnY5oZP9sll1xi0umqg63ojlfQ2XrNz86UNsfvENfp169fud8hpsRzPe5vqso+GTJkiL31nPg+9Je//MWk35VebvverV69Glu3bjV/DzwetrKyxZ7b/Pnnn8t879i64Yjfeb6WGR4iIoFm6NChpiWXM6KwdZn1/WeffWa6kbFe5G8oz0Mc6wNmNrVp06ZEfVBe/VxVVa1r2VLOrCrH33TWkcx0K43ZUuzCxcymirLESp/D8LyLZXGVv/3tbyUes/xsgWa9n5WVZV/Oep7Hga3i5WE2GVPMHVv2mRnH1mvHcxgeG2aCEfcP96uty191zmFEqktBt3gNVoasFDl4GitA/lhefvnlTtdlAMLUc6bL8nWON/6QM1XKhinA7PvEoI5BK9exBdbchiOmeDlWasR0rtL9xEvjD/6ZZ55p0siYBsUKnYFxdQJwVkDOBhHhe7GP2dKlS+0pcuxb5ljpOMP+0u3bty9RYfE+06BtfduIFz6Ygs0TE1ZefJ77jAFq6f1VFTxxccQUO1aMjv3XmDq/e/duky5NP/zwg0kNY+p5dTHNrPR3iPjZS3+HeDHC8TvEPmI8geGx5YkS17Ftryr7pPRFEX4nifva2XLb985WVqbCly4r08bz8vLKlKP0e/F77LhNEZFAwrRmXmxlAM1AjRc1bd2q+BvLi72sp0r/xm7cuLFEfeCsfubvb0pKiv3GOuNUqlrXVvY3neOXsFscb6fqx23D1Gte/GWfaZaDF8DZfYl1y+ligOw4fooNz1X42dnVinjOxiCcwXjp8y9HLBcvMvPcyvEchu/DgNyG513s4sVj6bhfOTZLdc5hRKpLfbrFq7Aljy2SrLTYouvYh8oRf1QZcM+cOdPp87Z+tKy8WJkw2GZfawZ6HKiKVztZIZUOip31JSLHVmVnGIixtZGVOfsqzZs3z1QGDOgYwHG75VUmpQd0c9ymMyNHjjR9o1jxcNoT/svg1bF/VHlY2bHfGq/gs2WVlR5bellpObZAs7LlyK7sZ88AkGXnhQRXteKTs/3BEyAGthzMhv3o+C9bGngxprpK70/bZ2GfO2ej4zvuE7Z+sO8dT2C6d+9urprz9ez/X5V9Ut7361TfO9t7sM8+398Zlqkq2xQRCSTMWnMcbNQRf2NZJ3EsGWe/naV/X0vXJ7wo6zhQKc87Sg84VlpV69rK/qZzcFme+7Bu49g3pS84O8P35Rg2HKOEfas5wCrHFXnxxRfNMn7+qp7DOLY4O2KWGDO+eO7Ccz6+H4PwUzUcEPcNMwyY/cW6kNtgIG4bQ8U2RRv75rP8HDSXA+KxHNzPrjyHEakqBd3iVdgizUqCP/IVDQ7C4JmtoGxdLi84JVZ6TC1iyzmDOBvH0UpdhT/q/PHnberUqeaHn4N6MRBn0Gi7Ks3K0JGz1LCKMM2ag28xxZ7vw/3ENLPKDA7CSo1Xsj/99FMT3DLVmJWYI1a8bFFlZeuYUl263FXFlgTHyp8DzbECdEy35kkFK2EOfse0OA6GUjotvDwVXSF3xjaYDS/eVBTUsxVh/vz5Zr+xVcLx81S3DFUtKy8eueIChLvLKyLiS/gby+CVdVTbtm2r/HoOhunYNc1W31f0O+uuupYBKLfNVG2ejyxevLjSg4cxIOaNF+eZdcgRx2fPnm2y+Fx1DmO7kD1t2jRzDsJzGJ4H2LpsVYRdzXiOaDs/5OClEyZMKLEOPztnoPn3v/9dYjnL7Rici9Q0pZeLV+HVVI7eyWmP2KJb0Q82r67yKqazETNtlYItWHO8EsypyDhStCs5SyWztUja0rNsgRNbxG34GTi6aFUxeN6/f79JLeY0GJW5QkwdOnQwI7GzwuKtQYMGJS5G2PZZ6SvnHMG1vKvZVUntK71NKj0qPFPJGeiyYmXameOJzKkuRlTlZIWt6gxieXHEWR8221Qszr5DxJFjnZWBqnvSVBrHKuD3hyPAOvaDK13WqqrqPhMR8UdMT+ZvPS+ulv6t52PH6SGd4TSQvCBqu/E3+1S/s+6qa4lp3WyYYAsyx745VflZ55YuS+lzmGbNmpkyO57D0OmcT/GchdtldgAzA3lOVxnMfmTdzRZuXgxgij8D8VPtVzZS7Nu3r8rlFHEltXSL1yk9pZUzTN1iUMbpkJhmxKkuOEUYWx/548orqOwPzvRrXp3lNjmtGK84M+XK1Sm2TF1nRcQ5nlkxsf8XKyJWfLaBQZjyxSu5vCrLIJ0pT6w0eJGgqmzzXnIuSlYw7ItVlcqOLbZMs7/xxhvLpH+xFZ37iKluPJFg33FW3rZpVU4XswtGjRplUrK5TaaOs1W79NzcPXr0MFOG8TjyIkHPnj0rtX2e5PCCDadF4WA0bMF27KteGgNurs8gn+/BFn92S2CfcnYRYBYFp6DjerwwwYFpGJyzLx+7DDjLlrCdaDHDgdvjd5IXj6o7CByPES+w8AIFv0dMr2M5eBLBTAqWkSl6VcXy8tgyY4ItIWzlsQ3iJiISKHhRk3UH62eOM8JAjnUsf+c52BoHMWN9ezrK+511V11rw3qQdRWnIGOgymknWVc4w+CX5yzMNuS+4KCzb731llmf5xvEcrIbGy8M8FyK63HgtdL93SuDdS7Lx7qSwXdlGw6I6/JiPMvLz1W6GyL3K8/JWE/yHHDt2rWmKyKnTRXxKLeOjS5ShSnDKlJ6yjCbN99809KrVy8z7QenG+vSpYvlH//4h2X//v32dX755RfLGWecYdZp2LChef7bb791Ok0Xp50qjVNf8P0rMn/+fMvFF19sth8eHm7+HTNmjGXLli0l1tu+fbuZ9oNTV9SvX9/y8MMPW77//vtKl8XRNddcY59GpLx95myakq1bt5rX8bZ48eIyz6enp1uuv/56S506dSyxsbFmmqpNmzaV2V5VpwzbsGGD5fLLLzfHqXbt2pY777zTcvz4caev4RRofM0zzzxjqSxO9cXvCLfP19qmGjnVd4zl52fkNGGcsqxVq1aW6667zrJ8+XL7OpxShlNucYoxrsfp7Pgd43b5+RxxehZOiRIcHFxi2i/ev+OOO0qsa5uyhFONlS4Tl3/yyScllq9atcpy2WWXmalu+B3iMbnyyivN96/0/i49fZ2zKet4XM866yzzt8HnNH2YiATqeQZ9+umnZnpJTvPFW/v27c3v9ubNm6tUPzsq73e2snVteeV3Vgc7O1datmyZqRdZBmfTlNHKlSvNOUvTpk1N3VKvXj3LiBEjStSDxHqF031xylXW45x+dd26dU6nDOP+q8gjjzxiXte6dWunz5eeMsyG02va9uWHH35Y5nlOGcZpVRs0aGDWO/PMM81UbKW3pynDpKYF8X+eDftFREpipsJ9991nWhwqmgZNRERERMTbKegWEa/C64BMOWeKXem5UUVEREREfI36dIuIV8jOzjZTmDHQZh+sL774wtNFEhERERGpNrV0i4hXYCo5B5jhoCi33367mbJERERERMTXKegWERERERERcRPN0y0iIiIiIiLiJgq6RURERERERNwk4AZSKy4uxv79+1GrVi0EBQV5ujgiIiJOR/HPzMxEw4YNERwcuNfHVWeLiIg/1NcBF3Qz4G7SpImniyEiInJKe/bsQePGjQN2T6nOFhERf6ivAy7oZgu3bcfExcV5ujgiIiJlZGRkmAvEtjorUKnOFhERf6ivAy7otqWUM+BW0C0iIt4s0LtBqc4WERF/qK8Dt6OYiIiIiIiIiJsp6BYRERERERFxEwXdIiIiIiIiIm6ioFtERERERETETRR0i4iIyCn9/PPPGDlypJmLlAPGfP7556d8zYIFC9CzZ09ERESgdevWeO+997SnRUQk4ATc6OUiIiJSddnZ2ejWrRtuuOEGXHbZZadcf8eOHbjoootw2223YebMmZg/fz5uuukmNGjQAMOGDfOaQ2CxWLBk+xF8tyEVqcdy0bpeLC7q2gAdGlR+WtFDmXmYu2Y/lu9KR3hoMM5qUxcXdE5GZFiIy8q5KSUDc9ccwNbULNSPj8R5HerjzNZJAT/CvTfad/S4+T6s2n0UcZFhGNyuLs7rWB9hIWrrkqrLLyzGdxtSsHDLIWTlFqJn0wSM6NYQDeKjtDt9SJCFtY0Hr5pPmTIFK1aswIEDB/DZZ5/hkksuOeVV8/Hjx2P9+vVmTrRHH30U1113XZXmUouPj8exY8c0ZZiIiHglb6+r2NJ9qjr7wQcfxNy5c7Fu3Tr7squuugpHjx7FvHnzvGY/zFmxF/9evAP5hUWICA1BbkEREqLDMGF4B/RpnnjK1x/MzMUjn63D9oNZiAgNRlGxBUUWYGCbJDwyvKMJwqtrxa40PD13E47m5JtAPq+wyGz3+gEtcGWfJtXevrjO7iM5ePTztdiTnoOosBAU8MsAYFinZPz9/La6SCJVUlxswZRvN+P7janghFShwUHILSxGs8Ro/PPSzmhcO1p71MMqW08Fe8NV8+nTp1dqfdtV83POOQerV6/Gvffea66af/vtt24vq4iIiFTe0qVLMXTo0BLL2MLN5d7iSFYe/vPbboQEAc2TYtAgPhLNk6Jx7HgBPli605zwnspXfxzAtoNZaFI7Cg0TotAkMRr1aoVjybYj+H1nWrXLyDK8t2SnCbhZNmsZYxASFISPft+Dw1l51X4PcZ05K/ZgT1oOmiXy+xSFponRSIgKw4+bDmLdvgztaqmS1XuPYsHmg0iMDjPfJf7GMODemZaNz1bt0970IR5NL7/wwgvNrbJmzJiBFi1a4MUXXzSPO3TogMWLF+Oll17yqlQ1T2G2AG+VxRQ/3kRERFwtJSUF9evXL7GMj9kqcPz4cURFlU2NzMvLMzcbrutOm1IyTTDbxKG1iK34STHh2Hk4GwcyctEooeIUTgbWkaHBCHVIHY4OD0VRcR7W7juGM1vXqVYZUzNzsfNwDurEhpdoJa0TG2FaUzceyMCgNnWr9R7iGkwe/W1HGmpFhiEk+OSxqhUZiiM5+Vi3/xi6NI7X7pZKW78/A3mFxWgYGWZfxu9WbHgolv2ZhrvO1c70FaH+cNWcLd4CvPHGG5g0aVKld8XEiRPxxBNPaNeJiIhXmDx5cpXqsepiqmZwUBCKLRaEmORNKzZwc3mYQ+BUHqaU8/WObD33XNGHl9sIDgaKiksu53uyjI7BnXheGLsY5BU6f059uqWq36cTf9/8TXG86Ma/f1d0XZGa41NBty9cNfekW2+9FaNGjbI/5j4ZOHCguc+MgNL7R63cIiLiLsnJyUhNTS2xjI/Z581ZfU0TJkww47Y41tkcv8Vd2OqYHB+JA8dy0bh2lAliC4uLcSQ7H2e2SkLdWhGn3MbANnVNi/bx/CJEhVsHTkvPKTD3+zSvXe0yskW7W+ME/LzlEGIjQkyLOk+4UzJyUS8u0jwn3oFB0dlt6+LDZbtNv3uOEcBg6WBmHmpFhLrk+yCBpXfz2pi1LBSHs/Lt2S4cdyK3oBjntFOGiy/xqaDbF66ae1LpdHH2mbfp3r07YmJiPFQyEREJNP3798fXX39dYtn3339vlpeHU4vxVlOYBv63s1th6ndbsOtIDtg+zcYk9p2+cVDLSg16NbxLMlbvSTepnoUn+oBHh4fg0h6N0KWRa1KJbxzYAnvTj2PnkWywEZ2lqh0TjjvOboWYCL8/lfMpf+nVGOv2Z5gLMRxUj2IjQnFt/+ZolqTzMKma1vVqYUy/ppi5bBd2HMlhmzdCgoPRq1ltXNyjkXanD/GpX2pfuGouIiLij7KysrBt27YSg5tyUNPExEQ0bdrU1Lf79u3DBx98YJ7nVGGvvfYa/vGPf5hpxn788Ud8/PHHZkRzbzKgVR00GR2NRVsPIz0n3/ThHty2rglqKxu4Pz6iE37984jpsxseEmxGPe/aON5lI1UzWJt6ZTczZRCD79rR4RjUpo4ZtE28S0J0OJ65tAsWbzuMLamZZgTz/q2SqjQFndQsDkbIiyT8a+3eJMEcQ29yVZ8m5vfk1z/TTCt3++RaZqwIV05JKO7nU0G3L1w1FxER8UfLly83s4fY2C5ojxs3Du+9954ZyHP37t325znwKQPs++67D9OmTUPjxo3x9ttve+XApwxer+7X9LRfb+bmblvX3NyFgcDF3dWy5QvYtYDzcvMm3oup/5+v3ocPl+7G0eP5ZllSbARuOLMFLuicDG/Bi3edGsabm/gujwbd/nrVXERExN+cffbZ9gHCnGHg7ew1q1atcnPJRESqbuXudLy9aIe5b5vBIDUjF68v2IZmSdHKThCXCvb0VfMePXqYm+2qOe8//vjj5nF5V83Zus35vTl1mLdeNRcRERER75OenY+fNh3EDxtSTZAlgWn+xoMmXTs5LtLMAsBbg/hIZOQWmq4cIn7T0q2r5iIiIiJSU75dn4K3F/2JtGxrOnFcZBjG9G2KK3o3dlkffPENh7LyzLSBjvgdCAkCjmSdnPlIxBU0wZuIiIiI+L3NKZn410/bkJVbaNKJmyZGo7CoGO8u2WkGqZLA0q5+LeQXFZsp+Gw44jwHnW9VN9ajZRP/41MDqYn4InaT4O10p34TERGR6vt5y0GTOtwiKdreqs25zncezsaPm1LNKOMSOIZ1SjYp5juP5CApJtyMWZGWU4AmiVE4t0M9TxdP/IyCbhE3e+ONN6o0V/zEiRPxxBNPuLVMIiIigSY9p8CkeJZOIw8LDcbhLGu6uQQOzlrw+MiOeH/JTmxMyUAQgsxUfOMGNEe9WpGeLp74GQXdIm526623YtSoUfbHx48fx8CBA839xYsXl5ljXq3cIiIirteybiwsJ1KIOWgWMbU4r7DYzH0sgYfzp0++rAuOZOcjOCgIiTEl5+jmQGv/W70fP2xMRVZeIXo0rY1LezRC63rW9HO2jv+0+SDmrjmA/Udz0bxODEZ2a4ABrerU2GdgGX7YeBDfrDuAA8dy0aJODC7u1hD9Wipzw5so6BZxs9Lp4tnZ2fb73bt3R0xMjI6BiIiIm53Tri7mrtlv0olrR4WZFu+0nHzUj4vAhZ3VrStQ8XtQJzaizHJenHl+3iYzknlYSLC5fb32AFbsSsM/L+mM1vVq4ZMVe/HO4h0oLC5GdFgoft9xBGv3HsVdQ9qY9PWa8J/fduP9pbtQXGxBVFgIlv15BGv2HsV9Q9tiSAfNFe8tNJCaiIiIiPi9pNgITBrVGee0q2dt8bZY0L9lEiaO7ISmSdZ5mkVsVu9Jx5LtR0xA3ighCvVqRaB5UrSZZm7Oir04mpOPT5bvMSOgN0uMQd1aEWiWFIOCIosJhNlK7m4chf/TFfsQHhJkBgZkGZonxSC/sBizlu02/4p3UEu3iIiIiAQEBtfsx3vseIFpyawdbW3xFiltc0oWCoqKERtxMlxiCnqtiDCs2n0UW1MzcTSnwMzt7YiDsh3MyMPe9BzTGu5OW1IzcSw3H43io8qUISUjF/uOHjfp5uJ5CrpFREREJKDER4V5ugji5SLDgu39/hls25hU8vBwRIWHmlbuwiILwkJOvq6guNiMGRAR6rDQTZhOHhIcjMJiCxzfjq3tLIPtM4jn6UiIiIiIiDgZoMrd2z/Ve1T3eTl9fVskIi4yDCnHcu1zeefkFyK3oNj0lebgexw4jS3KzJogtoxzJPxODePQuHbJ1md3DQTXLDG6TBmOZOWja6N4JMdpFHZvoZZuERERERHA9IH93x/78e26FKQfz0fHBnG4rGdjdG+S4LL9cygzD5+u2Iuftx4yQXP/VnVwea/GaJgQVWLE7O82pJg0+E4N4/GXno3RpXG8eZ6vWbD5kCnnnrQc87oRXRvgvI71lSrvQo1rR+PWwS3xxsI/sSstB2zrZusx53PnCOahIcG4e0gbTP56I/akW59n2NuqbgxuP7t1jRyL8NATZfhmk70M1KZ+LG4d3ErfBy8SZAmwS2QZGRmIj4/HsWPHEBcXB3/GUbJjY61TGmRlZWmUbC+h4yIipxJIdVVFtB+kJvGU+MXvtpiplzhSdURoMDJzCxEXFYbHRnRAr2aJ1X4PBtETPl1j5oWOCQ8F47LM3CK0rheDZ//S1fTFfW7eJny/IRXhIcEmqMrILTR9zx8b0dFMWcVg+/UF26z9jcNDkZ1fhODgINw4sAWu7N3EFbtCHPDCxtI/j5iLIW3r10LvZrVNwG3DAdV+2XYEh7PykBwfiTNb1ynRD7wmpGfn45fth00Ld4MTZYhxKANbwTmi+a4jOYiLCkXfFkk1XsZAr6e0t0VEREQk4G09mIUFmw+idnS4vc83g+CdaTmY/fse9Gxau9othz9tOohNqZloUjvaBPaUGFOMPw9l49v1KejRpDYWbT2MxOhwE+zby3AkBx8v34O29WPx0W+7TYsmR8w2zwMmvZit55ymSv3VXatJYrS5lSchOhwXdfXslHO1Y8IxomtDp89l5BbguW82YcWudHOhhl/hhvFReGBYe3v2hLifgm4XOHY831wJ9TY52Tn2+/vScxCd772jc9aKDEV8VLiniyEiIiIBiiNB5xQUmamhbBhkc07vbalZyMwrNH18q2PDgQwTMNsCbgoNDkZoSBDW7j1mWr/Zolq/VBkYSG9OycSW1Cwcyc43gbgjBumpmbnYcTjbpanw4vs+/HWXmfqM3ym2fhcWFWPv0eN44bvNmPHXXogKd/+Ab6Kg2yUYcLPvDdN/vEl+7smgm1dowyO9cw7KuMhQjOre0OVBNyst26AS3iQ7r7Dk/VDv+t44so58qR9jERHxfxwJmqEwzx0YBDuOBM3AhOne1RUTHgpnpyZ8TzZA2AIgruNQBNNCyed544jZfOxYP/Mxg3d+BhGb4/lFpv9/rYhQe7o5U+M57zinE1u5O92koov7qaXbRRhws5+ON8nPPVkeli3c4l3lczdWYE9+td70xfEmhXnH7ffvnLUCoRHuH93ydDCV6vERnTxdDBGRatcFS7cfwZLth83F2K6NE3Buh3rVbrGsCrYssaWJ/ULzCorQrUkCzm1fD7UcyrBu3zH8tPkgDmfmoVW9WAztUN8+sFZlZOYWYP7Gg6bfZkRYCAa0SsKAVnXMxdNK278K2Po9kH0IqNsOaDcciHOesuox6TuBzd8AaX8CtRoAbS8A6nc8+fzxdGDLt9bPEh4LtBwMNBsIBJ86YO7dPBH14iKx/1iuCUq47xi0ZOYV4IIuyS65CD2wTZJJI2f/X1trNed6ZiA0qE1ddGkUj7q1IrD/2HEz93LwiTLk5Bfh4u4N0apuLDo0jMPynemmPGwxZ8DNVu4ujRLQpp51LB8R8+dQUGR+czg+gSNeuOGI7PzdkJqhoFv8GgPu7Yey4U2K8k8G3X8ezkZIeLFHyyMi4s8DY/1rwTZ89ccBE/iyL+OCLYcwf2Mqnry4s+kH6W7FxRa8+uM2zFuXYub3ZQi8cMsh/LjpIJ4c1Rnx0WH4Zu0BvL5gu0lf5skw+/Ry9OzHR3ZCu+Rap3wPDuT0+BfrsX7/MbN9tpJy+xd1TsZdQ9pUrh/yuk+BJa8CeZlAcCiw/Udg41fAhc9ZA3BvsG8l8N2jQFaqtYzFhcCmr4CzJwBtzgOyDgFf/x1I3QAEBQOWImuA3m00MOBu5mlXuHmmcN8zpA2mfr/FjARNDLwZjF/Tt5lLPgL7hV/Vpwk+WbHXjIjNRu/osBBc0r0RBrauY4Lsu89tg5d/2IrdDmXg9FWj+zY1x/Kuc9rgqa82YMeRbPMd56dqUScG9wxtY14vYpMQFYbGidHYkpJpHyPA1ljI7x0v4kjNUNAtIiIifmnN3mP4Zm2KGaU3ITrMPiXUuv0ZZgTocQOau70Mq/ak47v1KSfGLjlZBpbtyzX7zVRP7y7ZaVorWyRFm6CKLVAcOOvdX3Zg8mVdThk0f/nHftNSztZZjnZty3D7Zn0KBrWta0a8rlD2YWDZm9YgNqm1dVlxkbU1+bc3gYtehMcVFwNLpwNZB4HElieCagtwbLd1ebMBwJrZQMo6oHZzIOREgJFzBFj7KdDqXCC5yynfpl/LJEy/pqfJjmD3QQazvZvXLtEHuzp4LK/t3wwDWtcxqb2Mupn5wAHSbMeZz7WpX8tkRmTlFpopqHo2O1mGpknRmDamuyljakaeaRnnNFYajVpK40WY0b2bYMq3m7DzSLb5DcorKMbxwiKc37E+WiszosYo6BYRERG/tGp3ukkp5xQ6NgxKo8KCsXjr4ZoJuncfRV5hMRomhJUoA9M9mfLeLCkaR7LyzGjWtqArOCjIpB5vPJCBw1n5JqiqCKcr4vZsATfx5DotO9+8/ymDbrYgH08DEhxac4NDgJg6wP7VQE4aEF396bKq5egu4MhWILaeNeAm7i+mmGccAFLWAtt/AsJjTgbcFJVovXiwd3mlgm6qExuBkd3cl1bP48xgp6KAh8d8VAVliA4PxZAO9d1UQnE1ZiTwd4AXTsrr8sGLcfxKu+oCj81Zbeua95yzYq8JvBNjw80o95wbvibmEvfIBbqiPCA08pTZLTVJQbeIiIj4pYqH0vTsQJvmVNBibaytiKUS5axonep/ykoU0muwnOWdZPvKZxB/s+zPI5izci+2H8wy2Qjnnwh4bWMEcMT52b/txvJd6WA8zm4GV/VtivpxJy8WVhcHS+M4D+zjzQEBHecZ9xuF+cDaT4ANnwPHj1ozXti1pNUQrwi+/XCPi4iIiMC08HJQMceBTpnGfbygGAPb1K2RXdSjaYJpgc5wKANbtNjqxTTiro3jkRgTjoOZufbnmV7OaaE6NIhD3diKW7ltJ9TcHrdrw8/M1u+eTSsxfVTDHkBUbWtfaRtLsTXtnM95upWb2ArP1Heml7NsxIsBmSnW1m+2Yrc8B8jPAooKSg6sFhYNNO7jsqIwe4LzeTNQ+mFDKrIcZiURccQuAE9/vRErd6Wbr+uRrHzTbYTjBrD1e//R43j8i3X4bkMqioqKTer356v3Y+IX610+QDNbtZkh4ZcBNy15BfjlZSBjv3mIA38AP0yyjuvgBdTSLSIiIn6pW+N4XNCpPr5el4L0nHwEIQjFsKBTwzi3pg876tGktuk7OW99igmkmTrOMnRuFG/6cydEh+O6AS0wY+F20+JlpoMqtpiU+OvObF6p9M+RXRvi9x1pWH8gA8H8jBbrlFcXdm6Abo0rEXTH1gX63mIdSO3ItpODlMU3Bvrc7BWtRGb08f53WAdSS9sOBIcDxQVAZLx1OdPKu10F7P0dOLjRuj7TTJlq3uVKoH5nlxSD0yw9/dUGbDmYZc8AYB/rh4d3QOt6px70TgIHg+qPl+8x08M2SzzZfSQjtwC/bDuMTSmZWLr9MPamHzfP29LOOf7E1oOZ5sLOxd0befhT+Ij0XcCmudbfg+gk6zJeLOTyVf9nHWjRsduJByjoFhEREb/Ek9w7z21jWrx5kmubMmxox/r2Qc3czYxGPaSNGQiL04axRbp7k3jTH9c2bdlFXRugaWI0ftyUikOZ+WhVL8akoHJgtMrgKOxPX9bFtLr+sfcYIkODzcBattGwK6XL5UBiC2ALpww7CNTrALS70Bp4e4tGvYCL/2U9ubZNGcYyJp8IqNniPXKatWXLccqwFoNdcuGAQdT0n7ZhY0oGGidEmwwGZk7sPJxjRhufdlWPqk3RJn6NsxHwQlpCVHiJi2ecM5vjLWw/lIW1+zJMRorj94Z9unk9Z2tqpodK7oMOb7XOvMDfMEcMvNnyzSweD/+WKegWcbP8jCMoyEyzPy4uzLPfzzmwHcGhJVMHw2olIjzuxFU6ERGpFp7MciAh3jyF6Zxnt6tnbuXp0jje3E4XA/jLejY2t2oFtbx5M55UD7iz/Od5kt3jGuvNxdjKvXbvMdSJibAPWscAqX5chAmgNqdkomPDOJe/r/gmBtORYcHIzisqsbyo2GKC8JgTsyrwwk3pizvMoYg9cVFOKiEiFggJ5by81gHUbHjOzRZudjHxMAXdIm528Pe52P/j/zl9buOb95VZ1vDca9F4yFgdFxERES9yPL8IhcWWMqNLc2CqwiILcvLVt1tOiggNwbnt6mH273uQHR5iguzComLsO5Zruo/0aZ5oupswA4ZdTxKjw0ywfSgzD7ERISZTRSqpYU+gdgtri3dCU2ugXZBjHdOh02VeMS6Fgm4RN6vX5yLUbt+/0uuzpVtERES8S5PEaNSJDcfhzDxEhZ9M/U/LyUd8dBha1i1/CjAJTGP6NcWutBys2JWOQ1l5ZlyJ+vGRGH9eWzOS+aDWdXBFr8b4YvV+7EzLMWPv14oMw9j+zcy4D1JJoeHAOY8APzwBHN1tna2A0x426Qec8Td4AwXdIm7GVHGli4uIiPg2TvE0uk8TTP9pO3alZSM2Iszeus0pnjgKvYgjBtBPXtwZa/YeNf27+bhfy0T7eA4cc+HGgS1wbvt6WLP3mOkO06tZbTSs5HgO4oBjO1zxHrBrCZBzBKjdDGjc15p27gW8oxQiIiIiIl7uoi4NTJrw/1bvx+60HLSpV8uMQj+sU7KniyZeioE0B3PkzRn272aWhDIlXNS3u+358EYKukVEREREKoEB0jnt6pkbB8TSaOUiUhl+Oju6iIiIiIj7KOAWkcpS0C0iIiIi4mc49RTng9ao6iKep/RyP5KVfsjcbArycu33U3dsQliEw7x1nP+vdl1zExERERH/sWT7YXz02x7sTMtGWHAwBreri2vPaIaEaA32JuIJCrr9yKrvPsEvH89w+tzMR68rs+zMK2/DoNG310DJRERERKQm/L4zDc99swlZeYWoHR2OgqJi/HflPjPw2zOXdikzz7iIuJ+Cbj/S4/wr0KbP2ZVeX63cIiIiIv7l0xV7kZlbiOZJ0WbgN4oJD8Ufe45i5a509GuZ5OkiigQcBd1+ROniIiIiIoGrsKgYWw9mIS4qzB5wU1R4CIoyLdh5JFtBt4gHKL9ERERERMRPRlRPiA5DXkFRieWc3ozio8I8VDKRwKaWbhER8RoHDhwwt8pq0KCBuYmIiHUe8Qs6JWPGwu04mpNvguzCYgv2H8tF/bhInKHUchGPUNAtIiJe44033sCkSZMqvf7EiRPxxBNPuLVMIiK+5OLujbAnLQc/bj6InWk5CEYQGiRE4r7z2mr0cpFADbqnT5+OKVOmICUlBd26dcOrr76Kvn37lrv+yy+/jNdffx27d+9GnTp1cPnll2Py5MmIjCw5HZaIiPieW2+9FaNGjbI/Pn78OAYOHGjuL168GFFRUSXWVyu3iEhJ4aHBJsC+uEcjbE3NQkxECHo2rY2YCI+f9osELI/+9X300UcYP348ZsyYgX79+pmAetiwYdi8eTPq1atXZv1Zs2bhoYcewjvvvIMBAwZgy5YtuO6660wqzdSpUz3yGURExHVKp4tnZ2fb73fv3h0xMTHa3SIip8Bz41Z1Y81NRAJ8IDUGyjfffDOuv/56dOzY0QTf0dHRJqh2ZsmSJTjzzDNx9dVXo3nz5jj//PMxZswY/PbbbzVedhERERERERGvDbrz8/OxYsUKDB069GRhgoPN46VLlzp9DVu3+RpbkP3nn3/i66+/xvDhw2us3CIiIoGKXcJ40ZtdupihdqqL3sxga9eunekW0KRJE9x3333Izc2tsfKKiI87uAlY9iaw+CVg8zdAfo6nSyTiW+nlhw8fRlFREerXr19iOR9v2rTJ6WvYws3XsX+fxWJBYWEhbrvtNjz88MPlvk9eXp652WRkZLjwU4iIiAQGdQkTkRq1dg6wdDqQd+LcPSgEaNgDuGAyEJ2ogyE+xafm6V6wYAGeeeYZ/Otf/8LKlSvx3//+F3PnzsVTTz1V7ms4yFp8fLz9xivtIiIiUjXqEiYiNSZ9F/Dr60BxIZDYCkhqDcQ1APYtB1bP0oEQ/w+6U1NTce2116Jhw4YIDQ1FSEhIiVtlceRxrs/tld5+cnKy09c89thj5r1vuukmdOnSBZdeeqkJwhlYFxcXO33NhAkTcOzYMfttz549VfzEIiIivslVdXZNdQljZhoz0hxvIhKA9iwDco8BtZI5Kpx1WWgkEB4DbP0esFg8XUIR96aXc7RwTtfFAJgjzHJ0xNMRHh6OXr16Yf78+bjkkkvMMgbOfHznnXc6fU1OTo6p5B3ZThqYbu5MRESEuYmIiAQaV9XZNdUljBfRqzJPu4j4qUJb19BSv1lMMS/KAyzF1vsi/hp0c57URYsWmalbqovThY0bNw69e/c2c3NzwBVOD8PRzGns2LFo1KiRqYRp5MiRJr2tR48eZgCXbdu2mRMJLq/KFXsREZFA4Mo6uzpdwmx19j333GO6hLHuLi87jecGNmzpVrcwkQDUoBsQFmVt7Y5KsC5joM3HHS8GgnXeL34edLPyK69VuapGjx6NQ4cO4fHHH0dKSoo5KZg3b579Sjqvzju2bD/66KPmKj3/3bdvH+rWrWsC7qefftol5REREfEnrqqzq9sljNgtjBfWb7nlFjzyyCNlMtdI2WkiYiR3AdpeCGz8Asg9CoSEAfnHgfjGQPdrtJPE//t0szX6oYcews6dO11SAKaS79q1y/TjWrZsmbka7niV/L333rM/Zn+0iRMnmqvlx48fN0E5py9JSDhxBUxERERcXmc7dgmzsXUJ69+/v8u6hImIGOwKM/gB4OwJQKNeQEIzoMc1wMhpQJ3W2kni/y3dbJ1mRdqqVStER0cjLCysxPNpaWmuLJ+IiIicJlfW2eoSJiI1iq3bnS6x3kQCLejmVXMRERHxfq6ss9UlTERE5PQEWQIsx4uDsnC+bk4fFhcX55Jt7k3PwYe/7sax4wUu2V6giY8Kw1/PaIrGtaNdut3svELc/MHv2H4o26XbDRSt6sbgrbF9EBNR5WtzIi7DPsCxsbHmflZWFmJiYgJi77qjrvJF2g8iIuIP9dRpnU1z2pDPP/8cGzduNI87deqEUaNGaQRxERERL6M6W0RExLOqHHRzELPhw4eb0cPbtWtnlnFKL46QOnfuXNNvTERERDxPdbaIiIgPjl5+9913m8B6z549WLlypblxFPEWLVqY50RERMQ7qM4WERHxwZbuhQsX4tdff0ViYqJ9WVJSEp599lmceeaZri6fiIiInCbV2SIiIj7Y0h0REYHMzMwyyznADefxFBEREe+gOltERMQHg+4RI0bglltuwbJly8CBz3ljy/dtt91mBlMTERER76A6W0RExAeD7ldeecX06e7fvz8iIyPNjWnlrVu3xrRp09xTShEREaky1dkiIiI+2Kc7ISEBX3zxBbZu3YpNmzaZZR06dDBBt4iIiHgP1dkiIiKed1rzdFObNm3MTURERLyb6mwREREvD7rHjx+Pp556CjExMeZ+RaZOneqqsomIiEgVqc4WERHxwaB71apVKCgosN8XERER76Q6W0RExAeD7p9++snpfREREfEuqrNFRER8fPTyG264wek83dnZ2eY5ERER8Q6qs0VERHww6H7//fdx/PjxMsu57IMPPnBVuURERKSaVGeLiIj40OjlGRkZsFgs5saWbs7PbVNUVISvv/4a9erVc1c5RUREpJJUZ4uIiPhg0M25PoOCgsytbdu2ZZ7n8kmTJrm6fCIiIlJFqrNdK7+wGL/vTMOetBwkRIfhjJZJSIgOD7zvZUEusHsJcHQPEJ0EtBgERMZXbRtpO4A9y4CiAqBhD6B+J55EuqvEIiK+FXRzYBa2cp977rn49NNPkZiYaH8uPDwczZo1Q8OGDd1VThEREakk1dmucyQrD/+cuxHr9h2zZvwBSI6Pwj+GtUO3JgmB853MOAB8+zBwcANgKbYui28CnPckkNz51K+3WIDVs4Dl/wZyM6zLwqKBDiOAgeOBkEqfkoqI+JxK/8INHjwYhYWFGDduHHr37o0mTZq4t2QiIiJyWlRnu867v+zE6j3paBgfhciwEBQVW7AnPQcv/bAFM/7ayywLCEteBQ78ASQ0AUIjgeJCIH0X8NPTwBXvA6GnaPlPWQP89qb1fmJL67+5R4F1/wXqdwbaD3f/ZxAR8YWB1EJDQzFnzhzTh1tERES8l+rs6svILcCS7YcRHxVuD65DgoNMAL4v/ThW7k5HQMg6BOz51ZpSzoCbgkOBuEZA+k5rMH4qOxYBeVlAbH1rOjlvUbWtrebbvnf7RxAR8anRy5levnDhQveURkRERFxGdXb15OYXobDIgvCQkn2OQ0OCUGyx4Hh+gDRC5GdZW7ZDSrVmh4RZlxdkV24bVLr/NtPKbenmIiJ+qsodaC688EI89NBDWLt2LXr16oWYmJgSz48aNcqV5RMREZHTpDq7epJiI9A4MRpbUzMRGxFqBo2lozkFiIkIRdv6tQLju8kW7VoNgKO7gIjYk8tz0oCIOKBOu1Nvo257IDgYKMw/mYpeXAQU5gKNerqv7CIivhh033777ebfqVOnlnmOlZFSz0VERLyD6uzqYSr5mL5N8Py8zdiVloNakaHILShGYbEFl/ZoiCaJ0QgIDJJ7jgMWPmsdfTyiFlCQY00N5/K4BqfeRushwMYvgf2rrK8PCgZyjwG1mwEdL66JTyEi4jtBd3HxiRErRURExKupzq6+QW3qIjwkGJ+t2odtB7PQJDESw7skY2TXAJuxhQOdhUcDaz4GjmwHklpbg+VOl1bu9Qy0L3weWD3T2oe7qBDoPBTofjUQ39jdpRcR8SjNzyAiIiJSgX4tk8yNI5cHmzHAAnRe6ZZnW29MCw8+jVHbY5KAM+8GBtxlnUKM6eYiIgHgtH7tOJDayJEj0bp1a3NjP+5Fixa5vnQiIiJSLaqzXZtuHrABt6PTCbgdcR8q4BaRAFLloPvDDz/E0KFDER0djbvvvtvcoqKiMGTIEMyaNcs9pRQREZEqU50tIiLig+nlTz/9NJ5//nncd9999mUMvDmw2lNPPYWrr77a1WUUERGR06A6W0RExAdbuv/880+TWl4aU8x37NjhqnKJiIhINanOFhER8cGW7iZNmmD+/PmmL7ejH374wTwnIuILDhw4YG6V1aBBA3MT8SWqs0VERHww6P773/9u0slXr16NAQMGmGW//PIL3nvvPUybNs0dZRQRcbk33ngDkyZNqvT6EydOxBNPPKEjIT5FdbaIiIgPBt1/+9vfkJycjBdffBEff/yxWdahQwd89NFHuPjii91RRhERl7v11ltNtxib48ePY+DAgeb+4sWLzQCRjtTKLb5IdbaIiIiPztN96aWXmpuIiK8qnS6enZ1tv9+9e3fExMR4qGQirqU6W0RExAeDblq+fDk2btxo7nfs2BG9evVyZblERETERVRni4iI+FDQvXfvXowZM8b0405ISDDLjh49avp3z549G40bN3ZHOUVERKSKVGeLiIj44JRhN910EwoKCkwrd1pamrnxfnFxsXlOREREvIPqbBER8QmZKcDSfwEf/RWYcxOwehaQf7LrX8C1dC9cuBBLlixBu3bt7Mt4/9VXX8WgQYOqXIDp06djypQpSElJQbdu3cx2+vbtW+76bFV/5JFH8N///tcE/M2aNcPLL7+M4cOHV/m9RUTkhIJcoLjQ+3ZHnkOFm5cFhFrgtYJDgbBIeBNX19kiIiIul5kKfHUfcHgzEBoNWIqAlD+AvcuBC54FQsMDc55utnSXVlRUhIYNG1ZpWxzxfPz48ZgxYwb69etngudhw4Zh8+bNqFevXpn18/Pzcd5555nn5syZg0aNGmHXrl32NHcRETlNDLjnTQCO7vKuXZjncCFgzvVAxGkPReJeCc2ACybD27iyzhYREXGL9Z8BhzYDiS2sF7ApPwfY9Quw82eg9VCf3/FVPnthq/Rdd91lWqh79+5tH6DlnnvuwQsvvFClbU2dOhU333wzrr/+evOYwffcuXPxzjvv4KGHHiqzPpezdZtX7cPCwsyy5s2bV/UjiIiIMwy4D2/xrn2TV3zy/uFtQESVe0UFNFfW2SIiIm6xeykQFnUy4KbwaCCzEDiwJjCD7uuuuw45OTmmZTo01PrywsJCc/+GG24wNxsGyOVhq/WKFSswYcIE+7Lg4GAMHToUS5cudfqa//3vf+jfvz/uuOMOfPHFF6hbty6uvvpqPPjggwgJCXH6mry8PHOzycjIqOpHFhER8UmuqrNFRETcJiwKKC4qucxyojtZiO+nlp9W0M0UcFc4fPiwSW+rX79+ieV8vGnTJqev+fPPP/Hjjz/immuuwddff41t27bh9ttvN6lzEydOdPqayZMnY9KkSS4ps4iIiC9xVZ0tIiLiNq3OBfatAPKzgPBYa8CdcwQIiwaa9Q/MoHvcuHHwFI6Qzv7cb775pmnZ5tzg+/btM+lz5QXdbElnv3HHlm72cRMREfF3nqyzRUREKqXDSGDfcuDPn62jmFNYDNDtaqBhT7/Yiac1Ig1bqD///HMzVRh16tQJo0aNKjfF25k6deqY9VNTU0ss5+Pk5GSnr2nQoIHpy+34Ph06dDAjnzNdPTy8bPpBRESEuYmIiAQiV9TZNppxRE4pfReQug4IjQAa9wUi47TTROTU6eXnP20dOO3AH0BwmLWFu0F3ICgIARl0M6Wb03Oxhdk2BQlTuNl6zEHQWrVqVantMEBmS/X8+fNxySWX2Fuy+fjOO+90+pozzzwTs2bNMuux/zdt2bLFBOPOAm4REZFA5qo6mzTjiFSI/TGXTreOQpyXaT1RrpUMnPUA0OIs7TwRqVhIGNDybOvND1V5GNi7777bVNJ79uzBypUrzW337t1o0aKFea4qmPb91ltv4f333zdX4P/2t78hOzvbPpr52LFjSwy0xuc50AtHXWWwzROGZ555xgysJiIiIu6rsx1nHOnYsaOZcSQ6OtrMLOKMbcYRtrLzojlnGxk8eDC6deumw+SPNn0FrJ4JBAVbp/2JbwpkHQR+egbI2O/p0omI+FZL98KFC/Hrr78iMTHRviwpKQnPPvusqVSrYvTo0Th06BAef/xxkyLevXt3zJs3zz64Gk8MbC3axCvz3377Le677z507drVzNPNAJyjl4uIiIh76mzNOCKntGmu9d+YOtZ/Q4KBhKZA2g7gzwVA96u1E0UkYFU56Gb/6MzMzDLLs7KyTivFm6nk5aWTL1iwoMwyThnGEwgRERGpmTpbM47IKbFVOzSy5DK2etPxo9qBIhLQqpxePmLECNxyyy1YtmwZLBaLuTEIvu2228zALCIiIuIdPFlnO844wjFcmN32yCOPmLT08rBL2bFjx+w3psWLj0jubJ3uxza3LhXmnUg3b+nJkomI+F7Q/corr5j+YWxxjoyMNDemqLVu3RrTpk1zTylFREQEnqqzT3fGkbZt25Y740h5LfNxcXElbuIjulwBRNW2ppOzZTv7MHB0N1C/kwZSE5GAV6X0cl4h5zzXs2fPNiOh2qYfYSXKClxERES8gyvrbM04IqfUoJt1yp/l7wJHtgBBIUCnS4G+NwPh0dqBIhLQqhx0s6Jev3492rRpo0BbRETES7m6zuaMI+PGjUPv3r3Rt29fvPzyy2VmHOEAp5ySzDbjyGuvvWYGPL3rrruwdetWM+NIVUdNFx/StB/QpC+QkwaEhgMRtTxdIhER3wu6OZI4K+4jR46Yf0VERMQ7ubrO1owjUimcnzsmSTtLRKQ6o5dzmpEHHngAr7/+Ojp37lzVl4uIiEgNcXWdrRlHREREaiDoZvpYTk4OunXrZvp4RUVFlXg+LS3tNIohIiIirqY6W0RExAeDbvbhEhEREe+nOltERMQHg24OoiIiIiLeT3W2iIiIDwbdVFRUhM8++8w+/UjHjh1x8cUXIzT0tDYnIiIibqI6W0RExLOqHCVz6pFRo0YhJSUF7dq1M8uee+451K1bF19++aUGVxMREfESqrNFREQ8L7iqL7jpppvQqVMn7N27FytXrjS3PXv2oGvXrrjlllvcU0oRERGpMtXZIiIiPtjSvXr1aixfvhy1a9e2L+P9p59+Gn369HF1+UREROQ0qc4WERHxwZbutm3bIjU1tczygwcPonXr1q4ql4iIiFST6mwREREfDLonT56Mu+++G3PmzDEp5rzx/r333mv6dmdkZNhvIiIi4jmqs0VERHwwvXzEiBHm3yuvvBJBQUHmvsViMf+OHDnS/pjPccRUERER8QzV2SIiIj4YdP/000/uKYmIiIi4lOpsETltWQeBHT8D+dlA3XZAo95AiAunB2aj3eEtwN7lABvyGvcBklpb74v4mSr/5QwePNg9JRERERGXUp0tIqdl+0/AwueB7EPWxyFhQNP+wHmTgIhargm4l/4LWPuxNainiFig2xig7y0KvMXvVLlPt4iIiIiI+Knsw8DPzwPH04HaLYCkVkB0orXVe9VM17wHt/XHTCA4FEhsab0FBQMrPwB2L3XNe4h4EQXdIiIiIiJitXOxNfCObwwEnwgVwmOBsGhg8zdAcbFrWtIL84GYOtZWbd5i6gKFudbnRPyMgm4REREREbEqyLH+y5ZnR0wxLzwOFBdWf0/lZQDBIWWXB4UAeZk6EuJ3FHSLiIiIiHdgX9/UDcCupUBmiqdLE5g4aFpIOJCfVfK45B4FkrsAoeHVf48G3azBu2MAX1QAWIqBBl2rv30RL+PCIQhFRERERE7T0d3AT88AqeusqcccsKv9CKD/Ha4J9KRyGvQAWgwGtn4H5GZY933uMSAqEeh+jWv2YvuLrKnqHL08Io5RvbWFu157oO0FOlISmEF3jx497HNyn8rKlSurWyYRERE5Tf5YZxcXF5tbafycjp/V2TqOgm39U/18XYvFYm7evK7jsTPrMsj+biJw4A+gVgMgJhI4fhRY/R8ERdVGUO/rTm+7Wrfq+4F3znkEqNMG2PS1NfBuNQToehWCGvaA7S+uKtst8x1mAH/hFOCPWWZQtSD+12GkNaiPTqzw+16Vv3v9RvjRb4TFO9c91W9zlYLuSy65pFIbExEREc/yxzp7+/btiI2NLbM8JiYGjRs3LrFeeSdA0dHRaNKkif3xn3/+iaKiIqfrRkZGolmzZvbHO3fuREFBgdN1IyIi0Lx5c/vj3bt3Iy8vz+m6YWFhaNmypf3xnj17kJub63TdkJAQtG7d2v543759yMk50dfWyYlsmzZtSqybnX1iGiYn2rVrZ79/4MABZGaW34eW27WdKKempuLYsWPlrsvystx08OBBHD16tNx1uR+4P+jw4cNI27wE2HcEiOwMFLLvMJ9JAIrT0Hzd/xDR/WrT4nrkyBFzKw+PG48fpaen49ChE1NeOcHvA78XxM/Fz1eeRo0a2b+DGRkZSEkpP/W9YcOGqFXLOq1WVlYW9u/fX+66ycnJiI+PN/d5zHjsylO/fn0kJCSY+8ePHzffn/LUrVsXiYmJ5j6/j7t27Sp33aSkJNSpU8fcz8/PN993o9YAoHd/630GHtlA4uHDZttUWFho/o7Kw7KyzMS/tW3btpVdqcElQPLFZh8kN2hgFvFveOvWreVul/uW+9imonX1G+FHvxFpaeWuy99g/hZTTf9G8G/cZUH3xIkTK7UxERER8SzV2eKTmL5sKbIO1uUoNBLIPWxNPQ5N8lTpAlcls2aq/R418T4ihXnA1h+A7L3WafASe9TYPgmyVNSm7od4dZJX03i1Ii6OfUiqb296Dj78dTeOHXd+FVwqFh8Vhr+e0RSNa1uvJLlKdl4hbv7gd2w/VP6VPClfq7oxeGtsH8REBMbQD7zia2vF4FVLXh0PKHlZwOyrrf3rvEh2XjFiH7G2YmQ93QYxEV46/medtsBVs4CIsq2x3lJX+SLbfmBrhLP9oNRRP0od3bMc+Ooea9pxuMPv77E9CIpvgqCrZgIhoV6dZqp1y98P5ImUcf1GeNfffbAn08sz9gPfPoygg5sQFGQxgwNaYuvDcs5jQNO+p71d1lO1a9c+ZX1d5bNppoe89NJL+Pjjj00KFdNQHFXU9C8iIiI1x1/qbJ6oOZ6sVbReVbbpr+uWDjR8Yt3GPYGGPYDdS61zN7OF+3i6dTTrrleYgNuryqt1q7QfvOVvQ+sG8Hf4txnWQRoTmlkHBywuQtDR3QhaPAW48v+A8OjT2m5lv1NVbjKYNGkSpk6ditGjR5uIfvz48bjsssvMGz7xxBNV3ZyIiIi4ieps8Rmcs/m8SUC7C63TSOWkAZEJQP/bgU6Xebp0IuLLctKs0xAypdw2EwJ/c+IbAcf2AvvdP6holVu6Z86cibfeegsXXXSRCbLHjBmDVq1aoWvXrvj1119x9913u6ekIiIiUiWqs8WnsIV72NPW+bk5cnl8Y5d12RCRAFaYe2LMiKiSy4NDTYs3CpwPaOlKVW7p5miNXbp0MffZ/9E2Qt2IESMwd+5c15dQRERETovqbPFJtZKt8zUr4BYRV4itD9RuDuQcNn257bI5Y0I8UL8jvC7o5tQcHDqe2ML93Xffmfu///67fah2ERER8TzV2SIiEvCCQ4DeNwDhtYC07UBWKpC+EyjIATpfDsSdnILOa4LuSy+9FPPnzzf377rrLjz22GNmfraxY8fihhtucEcZRURE5DSozhYREQHQ4ixg+PNAqyFARBzQoBtw7iNAv1trZPdUuU/3s88+a7/PwdQ4wfiSJUtM4D1y5EhXl09EREROk+psERGRExr3tt48oMpBd25uLiIjI+2PzzjjDHMTERER76I6W0RExPOqHHTXq1fPpKv99a9/xZAhQ6o0352IiBmRNi/D+3ZEds7J+0f3AAUVz9foUUyLikrwdCnEB6jOFhER8cGg+/3338esWbNw8cUXIz4+3qSYMwDv3dszTfUi4mMYcK+d432B9/H8k/dXfgBEnZjH0RsD7i6XK+iWSlGdLSIi4oNBN1u5ecvMzMScOXPwn//8x6SXt2zZ0gTfjz/+eJULMX36dEyZMsVMbdKtWze8+uqr6Nu37ylfN3v2bDNPOC8AfP7551V+XxHxEAbcx9O9a/fnFjjcPwogzJOlEXEJd9TZIiIiUjWnnRteq1YtXH/99WbKsDVr1iAmJgaTJk2q8nY++ugjjB8/HhMnTsTKlStN0D1s2DAcPHiwwtft3LkT999/PwYNGnS6H0FERCQguKrOFhERkRoMujk4y8cff4xLLrkEPXv2RFpaGh544IEqb2fq1Km4+eabzclAx44dMWPGDERHR+Odd94p9zVFRUW45pprzAkDr9aLiIh/OJBRiJV7c+231fty7c/xvuNzvHF9qbk6W0RERGogvfzbb781fbqZzh0aGorLL7/cXDk/66yzqvzm+fn5WLFiBSZMmGBfxoHZhg4diqVLl5b7uieffNIMDnPjjTdi0aJFFb5HXl6eudlkZHhZP1IREbF7Y+lRTPr+iNM9MvBfe8osm3heEp4YVkd7sAbqbBEREanBPt0jRozABx98gOHDhyMs7PT7PR4+fNi0WtevX7/Ecj7etGmT09csXrwY//73v7F69epKvcfkyZOVQici4iNu7Z+AUZ1iK71+g7gqV2MBxZV1toiIiJyeKp+tpKammr5hnsCBYK699lq89dZbqFOnci0bbEVnn3HHlu4mTZq4sZQiInK6GEQrkHYdT9bZIiIiUoWgm4FqXFycuW+xWCpM0batVxkMnENCQsxJgSM+Tk5OLrP+9u3bzQBqI0eOtC8rLi62fpDQUGzevBmtWrUq8ZqIiAhzExERCQTuqrNFRETEjUF37dq1ceDAAdOPOiEhAUFBQWXWYcXO5UwXr6zw8HD06tUL8+fPN4O72IJoPr7zzjvLrN++fXusXbu2xLJHH33UtIBPmzZNLdgiIhLw3FVni4iIiBuD7h9//BGJiYn2+84q8NPF1O9x48ahd+/eZm7ul19+GdnZ2WY0cxo7diwaNWpk+mZHRkaic+fOJV7PEwoqvVxERCQQubPOFhERETcF3YMHD7bfP/vss+FKo0ePxqFDh/D4448jJSUF3bt3x7x58+yDq+3evduMaC4iIiKerbNFRESk6qoczbZp0wZPPPEEtm7dCldhKvmuXbvM1F7Lli1Dv3797M8tWLAA7733Xrmv5XOcCkVERETcW2dPnz4dzZs3N5lnrKt/++23Sr1u9uzZpsXd1pVMPMRiAdJ3Aoe3AkUFOgwSGPJzgEObgYz9ni6J5GVaj0VmyfG8AkGVg+7bb78dc+fONf2r+/TpY/pSs4VaREREvIsr6+yPPvrIdAmbOHEiVq5ciW7dumHYsGE4ePBgha/jAKj3338/Bg0adJqfQlzi4Cbg878BH48D5twAfDwW+HOBdq7490WmNR8Ds68G5twIzL4GmPcQkKm4pcYVFwEr3gf+M+bEsRgDfP8EkJOGQFHloPu+++7D77//jo0bN5o5P3nVm1NwnX/++WYeUBEREfEOrqyzp06diptvvtmMudKxY0fMmDED0dHReOedd8p9DQdqu+aaazBp0iS0bNnSBZ9ITkv2YeDbCcDe34GwaCAyDjjyJzD/SWD/au1U8U+bvgIWv2z9/kfXBkLCgK3fA989pkyPmrb2E+DXfwG5GdZjERQCbPwfMH+S9eJIADjtztJt27Y1leiWLVuwaNEi0y/bNviZiIiIeI/q1tn5+flYsWIFhg4dal/G8Vb4eOnSpeW+7sknnzSjqN94443V/gxSDdt+AI7uAWq3sAbc4bFA7WbWE+D16qInftzKbSkG4htZLzZF1QbiGgGpa60XoKRmFOZbg+7gMCCugfVYRCcCtRoAe1cAKSVnpgrogdTKw75cs2bNMilnnAf0iiuucF3JRERExGWqU2cfPnzYtFrbBjm14eNNmzY5fc3ixYvx73//G6tXV74llWO78GZT0RzjUgXH9lj/DQ45uYyj2odFAWnbtCvF/+RnA5kHrBeZHPE7z1TnY3s9VbLAczwNyDkCRNYquTw8BshKtR6LBl3h76rc0s2r5OzPxavmZ555pklZe+6555CammoGShERERHv4Kk6OzMzE9deey3eeust1KlTp9Kv4/Sg8fHx9htT4cUFYuqx6c/a6ueo4DgQ31i7WPwPg+uoRGvw7agwDwgKBmLqeqpkgScyHoio5eRYHAdCwoFY/j75vyq3dNsGY7njjjtw1VVXlbnqLSIiIt7BVXU2A+eQkBATrDvi4+Tk5DLrb9++3QygNnLkSPuy4mJrwBcaGorNmzejVatWZV43YcIEM1ibY0u3Am8XaD0EWPMRkL7LmtLJqVizDlnTPNufPEYifoNZHZ0uBX55yfpdZz9iBtxs/a7bDmja39MlDKwLIO1HAr+9eaLFOx4oyLUOaNewO9CgOwJBlYJuppa98cYbuPzyy1G7dm33lUpERESqxZV1dnh4OHr16oX58+fbp/1iEM3HnPbTWbC/dm3JfnqPPvqoaQHnCOrlBdIRERHmJi7G1uwhE4HFL1r7drPFO7oO0Pt6oNkA7W7xT12usAZ5G76wfu/ZqsoA75yHgbBIT5cusPQaB+RlAJu/th6L0AigSV/g3EeAkGr1dvYZVfqUvMp91113mYFTFHSLiIh4L1fX2WyBHjduHHr37o2+ffvi5ZdfRnZ2tn1AtrFjx6JRo0YmRZzzeHfu3LnE6xMSEsy/pZdLDWnaD7jyQ+ugRcUFQL0O1hYnEX/FYG7AndbgO227NcW5Xidrpoe4XuoGYOt31myCpNZAuwtPdl8JjQDOuh/oNgZI32kd1I6/QRxbIkBU+dICK8s///wTLVq0cE+JRERExCVcWWePHj3ajHr++OOPm7m+u3fvjnnz5tlT1nfv3m1GNBcvFhoONO7l6VKI1Kxa9a03cZ8t3wILnwNyj1lT+zljAjMMLpgMJHc5uR5HkuctAFU56P7nP/+J+++/H0899ZRJNYuJiSnxfFxcqVECRURExCNcXWczldxZOjktWLCgwte+9957VXovERHxAXmZwNLp1oEZE1tZW6/ZhSXtT+vyS14PqBZtlwXdw4cPN/+OGjUKQQ470GKxmMfsQyYiIiKepzpbRDyGA5ixD2/qOiAyAWh1LtD0DAVg/ubAH9ZB0ZhKbosNbSPEH9oEZOzTLAmnE3T/9NNPrj9YIiIi4nKqs0XEIzhY1tf3A0e2WdONOTf25m+sg/f1uVEHxZ9YLCfulG7NDir1fGCrctA9ePBg95REREREXEp1toh4xMoPgMNbgcQWQPCJcCPrILDqQ2uLN5eLf2jQFYitC2QdAOIan0wvzz4ENOqpVu7TDbp//vnnCp8/66yzqrpJERERcQPV2SJS44oKgJ2LrKPj2wJuYrox+/nuXa6g25/wOPe7DVj0onWU+JAwoCgfiE0GzvibuhOcbtB99tlnl1nm2LdbfbpFRES8g+psEfE+Sjf2SewicHAjUJQH1GlrnYLNpsNIa4s2RzG3TRnWfoTrL64UFVr7iZsytAMiYuG3QXd6enqJxwUFBVi1ahUee+wxPP30064sm4iIiFSD6mwRqXFs6Ww2EFj/XyA68WRrN9ONGag17q2D4mtS1gE/TwGObAcshUB0HaDXdUDnv5xsyW7Yw3pzlwNrgEUvnChDERBTz1qGTpf6RGt6lYPu+Pj4MsvOO+88hIeHY/z48VixYoWryiYiIiLVoDpbRDyi1zggdY01QGLQXVwIhEZZg6TEljooviQnDfj+MeDobqBWA+vxzDkM/PIyEJ0EtDqnZkbC/+5RIGP/iTKEWC/iLH7Z2m2hxSD4XdBdnvr162Pz5s2u2pyIiIi4iepsEXGrhCbAqOnWKcNS1lqnDGvNKcP6a8f7mu0/Wkejr938ZNZCXCMgbQew4YuaCbq3/2gNuE0ZQqzLmM7OMQI2/s8/g+41a9aUeMz5uQ8cOIBnn30W3bt3d2XZREREpBpUZ4uIx3BEa7Z4i2/jqPNM33YcFI/CooH0nTVUhhTrv7aA216GKCB9F3xBlYNuBtYcOI3BtqMzzjgD77zzjivLJiIiItWgOltERKqlVrJ1rm12EXAMvAuyrVOC1YRaDa3/li5Dfo7PjIRf5aB7x44dJR4HBwejbt26iIyMdGW5REREpJpUZ4uIX2MwyBTjrFRrX18fCcB8CtPHV88C0ncAsfWtQW/2YWsrc+fLaqgM5wJ/sAw7HcpwCAiPsQ6k5o9Bd7NmzdxTEhEREXEp1dki4tcDfC18Dtj9K1BwHAiPBpoPBM76BxAZ5+nS+Y+o2sCwZ6zzcB/aaJ06jK3fva+37u+aEJN0ogxTgUObraOXxzUEet8IND0DviC4sisuXboUX331VYllH3zwAVq0aIF69erhlltuQV5enjvKKCIiIlWgOltE/B6nsNr2g7XFlQO3hUQAm78Bfpnm6ZL5n7ptgUtnAJe/C1z6BnDVTOvc3DWpXgfgsjeBy9+x/jt6JtB+OHxFpYPuJ598EuvXr7c/Xrt2LW688UYMHToUDz30EL788ktMnjzZXeUUERGRSlKdLSJ+jWnGu5da54uOjLemG0clWFtl//wJyDwx8Ja4DgdTS2oFJHe2XujwVBnqtAbqdwLCfKtrc6WD7tWrV2PIkCH2x7Nnz0a/fv3w1ltvmfm5X3nlFXz88cfuKqeIiIhUkupsEfFrnLe5INfap9cRHxfmWvv7ivhi0J2enm7m9bRZuHAhLrzwQvvjPn36YM+ePa4voYiIiFSJ6mwR8Wvxjax9uPMySi7Py7QG3rbRrkV8LehmwG0bBTU/Px8rV64004TZZGZmIiwszD2lFBERkUpTnS0ifo2DaLUaChw/am3V5kBqnE86NwNoO9w68JaILwbdw4cPN323Fy1ahAkTJiA6OhqDBg2yP79mzRq0atXKXeUUERGRSlKdLSJ+b+A9QLcxQFCIdSTzkHCg51jgjL95umQipz9l2FNPPYXLLrsMgwcPRmxsLN5//32Eh4fbn3/nnXdw/vnnV3ZzIiIi4iaqs0XE7zGN/Ky/A72us7Z2x9YDohM9XSqR6gXdderUwc8//4xjx46ZoDskJKTE85988olZLiIiIp6lOltEAgZTyZVOLv4SdNvEx8c7XZ6YqCtLIiIi3kR1toiIiA/16RYRERERERGRqlHQLSIiIiIiIuImCrpFRERERERE3ERBt4iIiIiIiIibKOgWERERERERcRMF3SIiIiIiIiL+HHRPnz4dzZs3R2RkJPr164fffvut3HXfeustDBo0CLVr1za3oUOHVri+iIiIiIiIiM/M0+1qH330EcaPH48ZM2aYgPvll1/GsGHDsHnzZtSrV6/M+gsWLMCYMWMwYMAAE6Q/99xzOP/887F+/Xo0atTII59BRHzPgbRsHEjPtj8+nldov7/6z0OIiij589igdgwaJMbUaBlFRERExPd5POieOnUqbr75Zlx//fXmMYPvuXPn4p133sFDDz1UZv2ZM2eWePz222/j008/xfz58zF27NgaK7eI+LY35q3FpNnOs2QGPjSnzLKJV/XFE1efUQMlExERERF/4tGgOz8/HytWrMCECRPsy4KDg03K+NKlSyu1jZycHBQUFCAxMdHp83l5eeZmk5GR4YKSi4ivu/WCLhjVr2Wl12dLt4iIiIiITwXdhw8fRlFREerXr19iOR9v2rSpUtt48MEH0bBhQxOoOzN58mRMmjTJJeUVEf/BVHGli4uISMDKOgTsWwHAAjTqBcSW7dYpIn6SXl4dzz77LGbPnm36ebN/tzNsRWefcceW7iZNmtRgKUVEREREvMjaOcBvbwLH0wGLBYiuDfS+Ceh6JRAU5OnSifgdjwbdderUQUhICFJTU0ss5+Pk5OQKX/vCCy+YoPuHH35A165dy10vIiLC3EREREREAt6+lcCSVwFLMZDQDEAQkJUK/DodSGoFNO4d8LtIxK+mDAsPD0evXr3MIGg2xcXF5nH//v3Lfd3zzz+Pp556CvPmzUPv3vphEBERERGplG0/APnZQFxDIDiEAyoBcQ2A/Bxg63faiSL+mF7O1O9x48aZ4Llv375myrDs7Gz7aOYckZxTgbFvNnGKsMcffxyzZs0yc3unpKSY5bGxseYmIiIiIiLlyDlsDbZL47LsQ9ptIv4YdI8ePRqHDh0ygTQD6O7du5sWbNvgart37zYjmtu8/vrrZtTzyy+/vMR2Jk6ciCeeeKLGyy8iIiIi4jPqtge2zQeKi04G38XFQHEhUK+jp0sn4pc8HnTTnXfeaW7OcJA0Rzt37qyhUomIiIiI+Jl2w4GNXwJpfwLRSdZlOUeAhCbW50TEv/p0i4iIiIhIDaqVDFz4PNCa0+1arAOqtTrXuiy+kQ6FiL+2dIuIiIhUhcViQWFhIYqKirTjvAxnpgkNDUWQpp7yXnXaAMOfd5gyLNHTJRLxawq6RURExKdwbJcDBw4gJyfH00WRckRHR6NBgwZmphrxYlG1PV0CkYCgoFtEREQqZfr06ZgyZYoZ+LRbt2549dVXzcwjzrz11lv44IMPsG7dOvOYU4Q+88wz5a5fWZxadMeOHaY1tWHDhiaoU4uqd2Ug8KIIB8nlcWrTpk2JAXFFRAKRgm4RERE5pY8++shM8zljxgz069fPTPE5bNgwbN68GfXq1XM6EOqYMWMwYMAAREZGmik/zz//fKxfv95MBXq6GNAx8G7SpIlpTRXvExUVhbCwMOzatcscLx5/EZFApkuPIiIickpTp07FzTffjOuvvx4dO3Y0wTeD3nfeecfp+jNnzsTtt99upgJt37493n77bRMsz58/3yV7W62n3k3HR0TkJAXdIiIiUiG2Vq5YsQJDhw49eQIRHGweL126tFJ7j/2vCwoKkJjouQGbiootWLr9CL5Yvc/8y8eB4KuvvsLZZ599yvVWr16N2bNn10iZREQCidLLRUREpEKHDx82o4TXr1+/xHI+3rRpU6X23oMPPmj6YDsG7qXl5eWZm01GRobLjsy8dQcw6csNOHAs176sQXwkJo7siAs6N4Av4ajtHB3c1Rh0f/7557jqqqu8pkwiIv5ALd0iIiLiVs8++6xpQf3ss88q7N87efJkxMfH22/st+2qgPtvH64sEXBTyrFcs5zPVxcHc7MNFNeiRQu8++679ueWL19u+rZ37drVPP/LL7843cZ1112HG264wazbtm1bjBs3DsePHy/x3FlnnYXOnTubZf/3f/9n+tf37NnTLP/jjz/McmYUMLWfg5jx/X766aeTnzklBeecc44Z2K5Tp0648847Tdr/wYMH8fjjj5t12SXgtttuM+t/++23Zvss++DBg7FhwwZ7n32+/sYbbzTr89iKiIhzCrpFRESkQnXq1DGjhaemppZYzsfJyckVvvaFF14wQfd3331nAreKTJgwAceOHbPf9uzZU+0jwxRytnA7SyS3LePzrkg1j4iIwG+//YZvvvkGd999t2n9ZWr+ZZddhokTJ2LNmjWmb/xf/vIXZGVlOd3GsmXLTKC7ceNGpKWl4aWXXrI/xxT/uXPnmuwCBu7/+c9/8PPPP2PlypV4+umncfXVV5v13nzzTTPAHQetW7x4sXneJiEhAV9++aXZFsuzc+dOfPzxx2YwvCeffNIE5GzxZp99BuLc5vvvv2/WveWWW3D55ZebEcqJZRw7dqxZ/4orrqj2/hMR8VcKukVERKRCnJaLLaOOg6DZBkXr379/ua97/vnn8dRTT2HevHno3bt3pYLWuLi4Erfq+m1HWpkWbkcMH/k816uua665xvzLgeOYas1WZQa/7P/Okd5p4MCBJi2fgaozV155JWrVqmUucrAV+YcffrA/x8CWz9EXX3xhWrbZ0s2W5rvuussE6WwZ53FhMMzjxhtbyB2PG1P9OeVbjx49TCt8eWXhBYAuXbqYm+3z7d+/H/v27TOPW7ZsaVq/RUSkYup8IyIiIqfE6cKY7szgmSnLnDIsOzvbjGZODPI4FRhTxIlThDFdedasWWjevLkJQCk2NtbcasrBzFyXrlcRx9R5Bs1s6XamKvOKO67ruN/Y2szjwZT2qmyDLe1swWZAzfLyuObmnt5nr8njKCLiy9TSLSIiIqc0evRokyrOQJotq2wdZQu2bXC13bt348CBk32jX3/9dZNazXTkBg0a2G/cRk2qVyvSpetVVbt27Uzr8vfff28eL1myxFyA4D50Zs6cOSb1nAPXsV94eQPPjRo1Ch9++KHZ78T3YKs18TV8jn27eQwc+5enp6ebLgEMuFmOTz75xP4cMwuY1m9zxhlnYO3atVi3bp15zH75vLBSnXnWRUQCkVq6RUREpFI46BZvznBgLUfsK+wN+rZINKOUc9A0Z7222QacHB9p1nMHpnf/97//NX28//73v5tgl4F1ea3Effr0Manohw4dMqn79957r9P1Bg0aZNL3L730Unvf8YsuushkInA+dQbKnE+9du3aZl324aZ77rnHXAjhIGilR5MfMmSIuSjCvvcczI39ujnfOrMY+B7cFoP0qrTUi4iIgm4RERHxYyHBQWZaMI5SzlDRMfC2hY58nutVh21wMcdp1mwYCLOFuzIY8L7zzjtllr/33ntllnFqL2fTe4WFheFf//qX0+03bdrUDPbmDEeML13OCy64wNxK47zf5fUFlwCRlwXsXAxkpQCxyUDzgUCEuhyIOKOWbhEREfFrnIf79b/2LDNPd7KPztMt4nFHtgPfPgKk/XniUlYQkNgSGPY0kNTK06UT8ToKukVERMTvMbA+r2OyGaWcg6axDzdTyqvbwu1KzlqzRbwOszoWTgGObAUSmgIh4UBRvvXxz1OAS17n6H2eLqWIV1HQLSIiIgGBAXb/VkmeLoaIbzuyDTi4AYitbw24if/yceoG6/N12ni6lCJeRaOXi4iIiIhI5eRnAcUFJwNum5Aw63I+LyIlKOgWEREREZHKSWoNRCcBOUdKLs9Jsy7n8yJSgoJuERERERGpnIhaQPdrAEsRkL4TyD5s/ZePe/zV+ryIlKA+3SIiIiIiUnndrgKiEoC1c4Bje4H6nYEulwNty04vJyIKukVERCRQFBcBu5YAWanWQZ+aDQCCQ+DNLr/8cowYMQLXXXfdKUc+P+OMM9C+ffsaK5sEMI5O3u5C6624GAhW8qxIRdTSLSIiIv5vw/+AeQ8CGftPLotrCFzwHNBxVI0Vo7CwEKGhrj/9YtCdkJBQ5aDbXeWRAKKAW+SUdFlKRERE/D/g/nhsyYCbMg5Yl/P5agoKCsKjjz6KHj16oG3btpg5c2aJ5yZOnIg+ffpgwoQJyMzMxM0334y+ffuia9euuOWWW5Cfn2/W3bRpEwYMGIBOnTrhkksuQUZGhn07s2bNQr9+/cx7dOvWDV9++aVZ/vbbb2P58uW477770L17d3z99dcoKirCAw88gM6dO5vbXXfdZX8PtprfcMMNOOuss8xzIiLiXgq6RURExL9TytnCDYuTJ08sm/eQdb1qYnC9atUqzJs3zwS5O3futD8XEhKC33//HVOmTMHf//53DBo0CL/99hv++OMPFBcXY9q0aWa9a6+9FjfeeCPWr1+Pp556CgsXLrRvY9iwYfj111/Ne3zxxRcmcM/Ly8NNN92E3r1746WXXsLq1asxfPhwvPnmm+b9VqxYYZZt377dPG/D5XPnzjVBvoiIuJeCbhEREfFf7MNduoW7BAuQsc+6XjUx+KWWLVuaVuSff/7Z/hxblm0+//xzE3yzVZqt1osWLcK2bdtMqzYDZFv/7S5dumDgwIH21+3YsQMXXnihaZ1mK3haWppZ5swPP/xgthMREWHSxxmgf//99/bnr7jiCtSqpVGmReTExcndy4BfZwC/vw2krAMszi5UyulSJx4RERHxXxw0zZXrVbHl2yY2NtZ+32Kx4NNPPzVp6I4cU8mdbeOqq67Cs88+awZXo8TEROTm5la5LKXLIyIBrDAf+PGfwPYfrPd5IXLlB9Zp4freYh00T6pNLd0iIiLivzhKuSvXq8C7775r/mVaOVuvmULuDFupn3vuOTOIGaWnp5uW7ri4ONPy/cEHH5jlTDFfvHix/XVcr0WLFub+hx9+aB7b8LXHjh2zPx46dKjZDvtx833Y7/v888+v9mcUET+z+Wtg6zwgMh5IagUktgKCQ4FVHwL7Vnq6dH5DQbeIiIj4L04LxlHKUV5rTRAQ18i6XjVx8DIGzQxuX3nlFTRv3tzpeuxbHRUVZdLLOZDakCFD7P2/GSizPzZTyDkwG9PUbdjvm63cfA/2627atKn9OQ7G9swzz9gHUuPjnj17mhuXsSz33ntvtT+jiPiZ7fMBSxAQEWd9zJbtmLpAQQ6w6+RFP6kepZeLiIiI/+I83JwWjKOUm8DbsZ/iiUD8gmddMl83B0jj4GelMZ28dGr3a6+95nQbnPJr6dKlTp/761//am42L774ov0+5/LmzdELL7xgbs6mFxMRMQqOl//7V1C57ityamrpFhEREf/Gebiv/ACIa1ByOVvAubwG5+kWEfEqjfsCRXlAsbW7iz3YDgoBkrt4smR+RS3dIiIi4v8YWLe/yDpKOQdNYx9uppS7oIXbWWu2iIhP6HixNcX88FYgPAawFAOFeUCTvkCrczxdOr+hoFtEREQCAwPsFs4HNxMRCUi16gMjXgLWzgF2LgJCwoDWQ4HOfwHCojxdOr+hoFtERER8TnFxsaeLIBXQ8RHxIexqc+bd1pu4hYJuERER8Rnh4eEIDg7G/v37UbduXfO49BzU4jlMs+c0ZYcOHTLHicdHRCTQKegWERERn8FAjnNVHzhwwATe4p2io6PNlGY8XiIigc4rgu7p06djypQpSElJQbdu3fDqq6+ib9++5a7/ySef4LHHHjNzWrZp0wbPPfcchg8fXqNlFhEREc9g6ykDusLCQjM3tniXkJAQhIaGKgNBRMRbgu6PPvoI48ePx4wZM9CvXz+8/PLLGDZsGDZv3ox69eqVWX/JkiUYM2YMJk+ebOajnDVrFi655BKsXLkSnTt39shnEBERkZrFlPKwsDBzExER8WYez/mZOnUqbr75Zlx//fXo2LGjCb6ZkvTOO+84XX/atGm44IIL8MADD6BDhw546qmn0LNnT7z22ms1XnYRERERERERrw26OdDGihUrMHTo0JMFCg42j5cuXer0NVzuuD6xZby89UVEREREREQCMr388OHDpi9W/fr1Syzn402bNjl9Dft9O1ufy53Jy8szN5tjx46ZfzMyMuAqmRk5CCs6jgiL+pWdjrCiQmRmZCAjpBCulJ1XiHqRxSiIdelmAwb3Hf9OiiJc/DORkQkUhgPF0a7dbqDgvuM+DHbdb5iRlwWENwBiCly73UDBfcd6JcI101jZ6iiOBB3IbJ/flXW2iIiIq1S2vvZ4n253Y9/vSZMmlVnepEkTj5RHnLtfO8YrTfN0AaQcj2jPeKU3Xb7FzMxMxMfHI1Dx85PqbBER8Wanqq89GnTXqVPHjHCZmppaYjkfJycnO30Nl1dl/QkTJpiB2myKi4uRlpaGpKSkgBhVk1dfeLKyZ88exMXFebo4coKOi3fScfFOgXhceMWcFXjDhg0RyPj5edxr1apV7To7EL9H7qJ9qf3oTfR91L70hfo61NNTfvTq1Qvz5883I5DbgmI+vvPOO52+pn///ub5e++9177s+++/N8udiYiIMDdHCQkJCDQ8wdBJhvfRcfFOOi7eKdCOSyC3cDuO89K4cWOXbjPQvkfupH2p/ehN9H3UvvTm+trj6eVshR43bhx69+5t5ubmlGHZ2dlmNHMaO3YsGjVqZNLE6Z577sHgwYPx4osv4qKLLsLs2bOxfPlyvPmm69P6RERERERERKrD40H36NGjcejQITz++ONmMLTu3btj3rx59sHSdu/eba502wwYMMDMzf3oo4/i4YcfRps2bfD5559rjm4RERERERHxOh4Puomp5OWlky9YsKDMsiuuuMLc5NSYWj9x4sQyKfbiWTou3knHxTvpuIi+R95Ff5Paj95E30ftS18QZAn0+UhERERERERE3ORk3raIiIiIiIiIuJSCbhERERERERE3UdAt1cZ+95w/9ejRo9qbXkTHpeY0b97czLwg3kXHRURERLyBgm4RERHxe9OnTzcXYiIjI9GvXz/89ttvni6SV/v5558xcuRINGzY0FxY50wxjjgkEGeeadCgAaKiojB06FBs3brVY+X1Vpzytk+fPqhVqxbq1auHSy65BJs3by6xTm5uLu644w4kJSUhNjYWf/nLX5CamuqxMnur119/HV27drXPx92/f39888039ue1H0/Ps88+a/7G7733Xu1LN1LQLfbKs7CwUHvDy+i41Jz8/PwafDepLB0XcYWPPvoI48ePN7N5rFy5Et26dcOwYcNw8OBB7eByZGdnm/3EixXOPP/883jllVcwY8YMLFu2DDExMWafMvCRkxYuXGgC6l9//RXff/89CgoKcP7555v9a3Pffffhyy+/xCeffGLW379/Py677DLtxlIaN25sAsQVK1Zg+fLlOPfcc3HxxRdj/fr12o+n6ffff8cbb7xhLmY40nfSDTh6uXjO4MGDLXfddZflgQcesNSuXdtSv359y8SJE0uss2vXLsuoUaMsMTExllq1almuuOIKS0pKSrnb3LFjB0ekt/znP/+x9O/f3xIREWHp1KmTZcGCBfZ1fvrpJ7PO119/benZs6clLCzMLCsqKrI888wzlubNm1siIyMtXbt2tXzyyScltj937lxLmzZtzPNnn3225d133zXbSk9Pt/gLHRf/OIZ33HGHucXFxVmSkpIsjz76qKW4uNg836xZM8uTTz5pufbaa83f1bhx48zyRYsWWQYOHGi+340bNzZ/n1lZWfbtpqamWkaMGGGe59/Jhx9+aLb10ksveeyz+hIdF/GEvn37mt8CG9Z1DRs2tEyePFkHpBJYx3/22Wf2x/wdTU5OtkyZMsW+7OjRo+Z8g+ceUr6DBw+a/blw4UL7fuM5mOO51saNG806S5cu1a48BZ47v/3229qPpyEzM9Ocz3///fembr7nnnv0nXQjBd0exi85A4InnnjCsmXLFsv7779vCQoKsnz33Xf2E4Pu3bubIGD58uWWX3/91dKrVy/zulMF3QwY5syZY9mwYYPlpptuMoHF4cOHSwTdDKr5Xtu2bbMcOXLE8s9//tPSvn17y7x58yzbt283ATUrUVvAvnv3bvN4/Pjxlk2bNpmAgxcK/DHo1nHx/WMYGxtrKhHbdzU6Otry5ptvmucZKPMYv/DCC+b7b7vx4hYDaP49/vLLL5YePXpYrrvuOvt2L7zwQku3bt3MyRD/JgcMGGCJiopS0K3jIl4qLy/PEhISUiJopLFjx5oL2lL1oJvnB1y2atWqEuudddZZlrvvvlu7tAJbt241+27t2rXm8fz5852eQzVt2tQydepU7ctyFBYWmgs84eHhlvXr12s/ngb+Bt57773mvmPQre+keyjo9jB+yRlQO+rTp4/lwQcfNPcZEPNkgcGuDX9c+AP922+/VRh0P/vss/ZlBQUFJgh/7rnnSgTdn3/+uX2d3NxcE5QsWbKkxPZuvPFGy5gxY8z9CRMmWDp27FjieZbVH4NuHRffP4YdOnSwt2zbvqtcZgu6L7nkkjLf9VtuuaXEMrZ8BwcHW44fP27ZvHlzmb89W4uEWrp1XMQ77du3z/yNlq7bmGHGFnCpetDNC5Jctn///hLrMRPvyiuv1C4tBxtSLrroIsuZZ55pXzZz5kwTOJbGc8F//OMf2pelrFmzxlwc57lxfHy8yb7Ufqw6XrDo3LmzObcpHXTrO+keoe5IWZeqKd2PgoOS2PqZbdy4EU2aNDE3m44dOyIhIcE8x8E5ysMBJmxCQ0PRu3dv8xpHXGazbds25OTk4LzzzivTp7JHjx728nAAmvLex5/ouPi+M844wwwO4vhdffHFF1FUVFTm+09//PEH1qxZg5kzZ9qX8XyzuLgYO3bswJYtW8zfUq9evezPt2/f3vw9io6LiEhF2Ld73bp1WLx4sXbUaWrXrh1Wr16NY8eOYc6cORg3bpzpBy+Vt2fPHtxzzz1mjAEOLCk1Q0G3FwgLCyvxmEECT/JrAgc+scnKyjL/zp07F40aNSqxXkREBAKNjov/c/z+2/4Gbr31Vtx9991l1m3atKkJukXHRXxLnTp1EBISUmY0aD5OTk72WLl8mW2/cR+yocCGj7t37+7BknmvO++8E1999ZUZFZ4DgjnuSzZucNpVxwu4+n46Fx4ejtatW5v7vADOgcCmTZuG0aNHaz9WEgeiY+Nez5497cvYGMHv5muvvYZvv/1W+9INNHq5l+vQoYO5IsWbzYYNG8yPM1u8K8KRMm04Mjn/yLi98nB7DK53795tftAcb7aWdr6+9DQrju8TKHRcfANH1C39XW3Tpo05AXeGFRD/vkp//3ljRc9Wbdvfkg2nftEc9Tou4r34t8uT8/nz59uX8cI2H/trppa7tWjRwgSLjvs0IyPD/OZqn5bEbCkG3J999hl+/PFHs+8c8bvJi/yO+5L1Cs/FtC9PjX/LeXl52o9VMGTIEKxdu9ZkDNhuzPy75ppr7Pf1nXQ9tXR7Oc572aVLF/OH8PLLL5sT/ttvvx2DBw8ukxpbGqf5YIDBAPGll15Ceno6brjhhnLX5xyS999/v5kmgD9iAwcONOk7v/zyi5kPkSk8t912m0nPfeCBB3DTTTeZ4OO9995DoNFx8Q08aeE0QWy95jRBr776qvn+lufBBx80Kek8QeL3my3hDMKZgsWrv0xru+CCC8z2OF8oU805ryXnqBUdF/Fe/B1gHcZ6s2/fvqY+5ZRN119/vaeL5rWY+cNuZzbsYsMT8sTERJP5w9++f/7zn+Y8g4HkY489Zub05jzUUjKlfNasWfjiiy/MeVZKSopZHh8fb+oO/nvjjTea7yj3Lc+37rrrLhNwsz6SkyZMmIALL7zQfP8yMzPNfl2wYIFpmdV+rDx+Dzt37lxiGc93OE+8bbm+k27gpr7iUkmOAxfYXHzxxfbpi6ozZdisWbPMIDEcoIODn/3444/2dWwDqZUe/IyDTr388suWdu3amSks6tataxk2bJh9agv68ssvLa1btzajmA8aNMjyzjvv+OVAajouvn8Mb7/9dsttt91mRinntCIPP/xwiSnDnA1+xkHSzjvvPDPyOf/mOML/008/bX/+wIEDZiAcfv85uuwHH3ygKcN0XMQHvPrqq+ZvlnUi60bOBiLls50nlL7Zzk/4W/rYY4+ZGUz4ezhkyBAz2KSU5Gwf8sbZYWw4mBXrK9ZTHND20ksvNXWNlHTDDTeY+pZ/wzw/5XfONtuP9qNrz3v1nXS9IP7PHcG8eM7OnTvNVedVq1apb5UX0XGpWWeffbb5/rNFS7yHjouIiIgEGvXpFhEREREREXETBd0iIiIiIiIibqL0chERERERERE3UUu3iIiIiIiIiJso6BYRERERERFxEwXdfjASMOfKLE/z5s01erOfeOKJJ0qMRn/ddddpPtQq0t9L4NDfi4iIiHiLUE8XQNzr999/NxPei4j+XkRERESk5ino9nN169b1dBFEfIb+XkRERETE1ZRe7gcKCwtx5513Ij4+HnXq1MFjjz0Gi8XiNL1806ZNGDhwICIjI9GxY0f88MMPCAoKwueff+7BT+Cb5syZgy5duiAqKgpJSUkYOnQosrOz7WnfzzzzDOrXr4+EhAQ8+eST5jg98MADSExMROPGjfHuu++W2N6DDz6Itm3bIjo6Gi1btjTHsaCgwGOfz1/p78Uz9PciIu60c+dOcz6zevXqKr0uPz8frVu3xpIlS7ymTDXVpap0Nxxxn8OHD6NevXrYu3evdnOAUtDtB95//32Ehobit99+w7Rp0zB16lS8/fbbZdYrKioywSCDumXLluHNN9/EI4884pEy+7oDBw5gzJgxuOGGG7Bx40YsWLAAl112mf1ix48//oj9+/fj559/Nsdj4sSJGDFiBGrXrm32/W233YZbb721xI9vrVq18N5772HDhg3mOL711lt46aWXPPgp/ZP+Xmqe/l5EhMFnRTcGgJ4wY8YMtGjRAgMGDNBBcpGabszxtvGLnI25w0axsWPHmvNBCUxKL/cDTZo0McEZf+TatWuHtWvXmsc333xzifW+//57bN++3QSIycnJZtnTTz+N8847z0Ml9+0ggi2mDLSbNWtmlrHV24at2a+88gqCg4PNMXn++eeRk5ODhx9+2Dw/YcIEPPvss1i8eDGuuuoqs+zRRx8tUYHcf//9mD17Nv7xj3/U+OfzZ/p7qXn6exER/g7YfPTRR3j88cexefNm+7LY2Nga30m8UP7aa6+ZbLTqYGt5eHi4y8pVU9v2JH/9XBW5/vrr0atXL0yZMsWcJ0pgUUu3HzjjjDNMwG3Tv39/bN261bRsO2LlxoDDFnBT3759a7Ss/qJbt24YMmSICbSvuOIK0yqdnp5uf75Tp04m4LZhmrljUB4SEmJS0g8ePFjiJOTMM880x4cnHwzCd+/eXYOfKjDo76Xm6e9FRFi32W7sDsfzFttjds265pprTF3J+q9Pnz6m+5sjXoxmty1mmDEzrGnTpiZjrzw8B+K67du3L7cuXbFihWmMuOiii0osZ+PFueeea+8+dssttyArK6tMSyYbLho2bGgurhMzDnv06GG68PXu3RurVq0q857r1q3DhRdeaD4nP++1115rUo8dU8LZZZBp4WwdHTZsWKVex33IllQ+36BBA7z44ouV/tK98cYb5vyQmZBXXnkljh07ZpYzWy8sLAwpKSkl1mfZBg0a5HRbPE506aWXmmNse2xLZWcmJjMLuI/o6NGjuOmmm8yYKnFxcWa///HHH/bt8fhcfPHF5X43uL927dqF++67z541QcwcZPe+r776yhwffrbLL7/cNIAw443lYvbh3XffXeJ8OS8vzzR6NGrUyAxE3K9fP9NYZWPb7rfffosOHTqYMl1wwQX2i0r8nNz+F198YS+P7fU8N+T35bPPPqv0sRH/oaBb5DQwaGbmwDfffGP6xr/66qvmR33Hjh3meVZSjvij62xZcXGxub906VJzwjF8+HBTQbCiZuo/rwSL+Dr9vYhIRRjQsv6bP3++qf8YxIwcObJMsMxA0hbM3n777fjb3/5WorXcMXDiBXH2pV60aJEJ0J3hcxxLhUG8Y/DKQJcBGWeA+eSTT0yQx0DYEcvK9+a5AOttfgZ2I+M5AYN5Bl8M3hwxwGRQycB8+fLlmDdvHlJTU02g64hBG1uBf/nlF5P+XpnXccyYhQsXmmDvu+++M4HeypUrT/nF27ZtGz7++GN8+eWXZru2fUtnnXWWGWPm//7v/+zrc6yZmTNnmgsaznCfEcetYSBqe2x7r08//RT//e9/7f3ceZzYAMHzKe63nj17mkaNtLS0Sn03uC2Ok8NsBb6fY0YFA2xmHTJrkJ+N+4QXA77++mtz4+fiBQeOOWLD48xzMr5mzZo1pnx8TzZmOW73hRdeMK/nhQmWxXas+S+Piy0Q582x6wIbu/i9k8Cj9HI/wD7Cjn799Ve0adPGnOg6YlC4Z88e80PNK4bk+GMoVcOgmS3TvDFNjmnmp3v1kgO48PWOfex55VZcT38vnqG/FxGpKBuGN5unnnrK1Kf/+9//SgS7DL5sASEHH2VXup9++sne0mwL0thyzcCbz7FVvTysZ9ny6GjWrFnIzc3FBx98YJ9ylSnoDPSee+45+/kTn2OrrS1Fmq3uvJD+73//27TislWT47bwwoANt8PAmS32Nu+8845pZd6yZYu5AEA8h2O3NJt//vOfFb6On4Hv++GHH5qA1Ra4Mxg9FdtnZcsusRGB+48XOJiFcOONN5oAmkE9MTjna0pfKCg9Cwhbgx0zK4kNCXwv2zrsYsfsAAbdERERZhmDWfYHZyDMDINTfTeYps3zXV44Kf1+vEDw+uuvo1WrVuYxW7oZKPM8mC3UvEByzjnnmO/J6NGjTfDMz8p/bd8LBtEM2Lnctv+5XV4MsW2X5bB1UeB2mSHB71/p8hC36ywDQvyfgm4/wB+H8ePHm4G5eFWTP5jO0orYd5s/EOPGjTM/5pmZmfZ+xI7p6VK5wI1XXc8//3wzGiUfHzp0yKQa8cpoVbGC5XHklVWmTs2dO1fpR26iv5eap78XEakIA2W2DLPus40Bcfz48TIt3V27drXft6WnO3bTIg5yymCTA5oy+KkI38OW5mzDwVEZ5NkCbuLFdQbUbNm2Bd3sMubYJ5mvY/kct8fufo6YNs0Az1n/daZR24Ju9vutyuv4ORjQMhXahsGo48WI8jALwBZw28ps+6zcv0yl57kiG3TYPYvp1Qy4HfdPZbFxwXFqTn4uHnum8Dvi5+Hnqsp3wxmmlNsCY+KxY1q5437kMtt3iN0KmGpuOw42DKAdy1h6u0znL/09LA+/k2wpl8CjoNsPsA8Pf4CYssKrfffcc4+5Olgan+PVQ/adYWDHlCEO5sCrt6UrHakY+x0xpYijZWZkZJiKhBc62N+KfbOratSoUaY/Eq+W8sedV5k5ZZinRnP1Z/p7qXn6exGRirA1kWnabOXk9F0MTNgqWbqLVUXdtBxbw9niyxRhpmRXhH2mGWidjtMJOhlA2lrMS2PgVt62T/U6pm27CxsW+N5s6WVfbKaBO/Zxrgpnn4vld7Y9tpRX5bvhTFW7+rE8PFdmmnvpbFHHQN3ZNmyz15wK0+YdLzxI4FDQ7eMcf6iYQuNsjkhHHFCE6Tw27C9E/CGTymOLNtONnOFV4IqOU3nHhtkHjulk5Di/JgNwxyDc2ftIxfT34hn6exGRivBchC2q7G9rC35K15GVxXTuzp07m4vZbB0dPHhwuesyZZvnTgyYbBl//L1i/cq+3bYgkeWzzUZSHr6OqctMvbY1ZLB12BH7K7NPM1tbOdVrZZ3qdWx1ZSDIrCJb/3UO7srU84o+P7HFmFOc2tKpWebSn5WNNbYMAr4XW/4rwrKUHsy3vM/FQdr4mWwDrp3Od4MZB5V5v1Ph94HbYat1eQPFVUZF5eGAeBz8TQKPBlILMOwHwyuG/MHiwCBsEeePp2OajIjo70VEpKawi5VtcC2mHF999dVlWrCr4q677jL9oDmwmWNDQ2nsz8sgbv369fZlHNSUQTO74jFAYlo3t8fRwm2p5c6wzAzcOV3rhg0bzEBdbJ11dMcdd5iWTgawHFOHKdQcBZtTSVUUNJ7qdWyFZd9r9rtmWj3LzUDVcRaV8tg+K/c7B/jiaN5MH3fsj8yB5ZixxH3K9zwVBtDsgseA2nFml9KGDh1q0tk5EjwHf+O5Kce44fg2HDCust8Nvh+zD/ft21diRPeqYlo5jz8z4vieHByXfc4nT55sLuBUFsvDroZM0Wd52AecmFbOVnR2TZTAo6A7wLAfN3+82eLNH2SmmXOkSxHR34uIiCdMnTrVjBbOUZ6Zyswgj62g1cFMsUmTJpl0cwZyzrCfLltQORq3Y39dBrQMcnmOxFRmDk7GQdAqwsCXg4wxXZ0tpgwcS6eDszWZLbcMlBl4sV84y8lU6ooC5Mq8jt0F2TrL/cdgduDAgWX6hjvDTMfLLrvM7Cdum/3S//Wvf5VYh+/Bc0a+PwPSU2F3OzbwcKA37ovy8CIFL05wlHQG8wx6r7rqKjPAne0CR2W+GxzEjAE7G5Cqm7rNNHp+xr///e+mtZ8XBHiho7wR8J3hhRe+liPtszy2rFKeb3M71WlFF98VZKlsJwQRERERET/CFkkONMvWY2cDlYkVW9I5YCxHDZfTw4HomEnA1noJPOrTLSIiIiIBiS27bJFmKjFbkKWkY8eOmdZ7TqWmgPv0Mc2cGQXsIiCBSS3dIiIiIiJSBgf9Yr9mTkvLedFF5PQo6BYRERERERFxEw2kJiIiIiIiIuImCrpFRERERERE3ERBt4iIiIiIiIibKOgWERERERERcRMF3SIiIiIiIiJuoqBbRERERERExE0UdIuIiIiIiIi4iYJuERERERERETdR0C0iIiIiIiIC9/h/j3s4yrWbKAQAAAAASUVORK5CYII=",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, axes = plt.subplots(1, 2, figsize=(10, 3.5))\n",
"\n",
"# Left: survival proportion by treatment\n",
"ax = axes[0]\n",
"for i, (pred_val, pred_label) in enumerate([(0, \"no predator\"), (1, \"predator\")]):\n",
" for j, (sz_val, sz_label) in enumerate([(0, \"big\"), (1, \"small\")]):\n",
" mask = (predator == pred_val) & (size == sz_val)\n",
" props = survivors[mask] / density[mask]\n",
" x_pos = i * 2 + j\n",
" ax.bar(\n",
" x_pos,\n",
" props.mean(),\n",
" yerr=props.std(),\n",
" capsize=4,\n",
" color=f\"C{i}\",\n",
" alpha=0.6 + 0.3 * j,\n",
" edgecolor=\"w\",\n",
" )\n",
"ax.set_xticks(range(4))\n",
"ax.set_xticklabels([\"no pred\\nbig\", \"no pred\\nsmall\", \"pred\\nbig\", \"pred\\nsmall\"])\n",
"ax.set(ylabel=\"Survival proportion\", title=\"Mean survival by treatment\")\n",
"\n",
"# Right: raw survival by tank, ordered by treatment\n",
"ax = axes[1]\n",
"colors = [\"C0\" if p == 0 else \"C1\" for p in predator]\n",
"ax.scatter(range(n_tanks), survivors / density, c=colors, s=20, alpha=0.7)\n",
"ax.axhline(0.5, color=\"gray\", ls=\"--\", alpha=0.3)\n",
"ax.set(\n",
" xlabel=\"Tank (ordered by treatment)\",\n",
" ylabel=\"Survival proportion\",\n",
" title=\"Per-tank survival\",\n",
")\n",
"# Manual legend\n",
"ax.scatter([], [], c=\"C0\", label=\"no predator\")\n",
"ax.scatter([], [], c=\"C1\", label=\"predator\")\n",
"ax.legend(fontsize=8)\n",
"\n",
"fig.tight_layout()"
]
},
{
"cell_type": "markdown",
"id": "6",
"metadata": {},
"source": [
"## Model 1: varying intercepts\n",
"\n",
"Following McElreath (Ch. 13), the first model uses a hierarchical Binomial\n",
"with partial pooling over tanks. Each tank $j$ has its own survival\n",
"log-odds $\\alpha_j$, drawn from a shared population distribution:\n",
"\n",
"$$\n",
"S_j \\sim \\mathrm{Binomial}(N_j,\\, p_j), \\qquad\n",
"\\mathrm{logit}(p_j) = \\alpha_j\n",
"$$\n",
"$$\n",
"\\alpha_j \\sim \\mathrm{Normal}(\\bar\\alpha,\\, \\sigma), \\qquad\n",
"\\bar\\alpha \\sim \\mathrm{Normal}(0,\\, 1.5), \\qquad\n",
"\\sigma \\sim \\mathrm{Exponential}(1)\n",
"$$\n",
"\n",
"A non-centered parameterization ($\\alpha_j = \\bar\\alpha + \\sigma z_j$,\n",
"$z_j \\sim \\mathrm{Normal}(0,1)$) helps the sampler."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "7",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-17T11:44:55.887359Z",
"iopub.status.busy": "2026-07-17T11:44:55.887285Z",
"iopub.status.idle": "2026-07-17T11:45:00.840814Z",
"shell.execute_reply": "2026-07-17T11:45:00.840393Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"