MachineSex/notebooks/01_biological_model.ipynb
Giorgio Gilestro c435cfba6e Reproducibility pass: figure map, one-command reproduce.sh, notebooks, Makefile gaps
An audit of the figure pipeline found real sync gaps, now closed:

- `paper/pnas/make_figs.py` (which draws every manuscript figure) was invoked
  by NO Makefile target or script - a manual step. Added `make paper-figures`.
- `configs/llm/epistasis{,_compat}.yaml` were reachable from nothing at all,
  despite producing Fig. 3C-D. Added `make llm-epistasis` (+ its statistics).
- `make figures` never regenerated the MNIST montage that Fig. 2B embeds;
  it now runs with the `mnist` target (it needs torch - it re-simulates).
- Added `make llm-society`, `env-notebooks`, `notebooks`.

New REPRODUCING.md is the authoritative map: every manuscript panel -> the
artifact it plots -> the config that produced it -> that config's seed, plus
the determinism policy (biological tier bitwise; GPU tiers statistical), the
seed-provenance statement, and an artifact-hash verification snippet. All 44
committed bundles currently hash-match their manifests, and figure
regeneration is pixel-identical (verified by comparison).

reproduce.sh delivers the one-command reproduction the paper's Methods
promises, writing REPRODUCED.md with recomputed hashes per bundle.

Two executed notebooks: 01 builds the Wright-Fisher model from scratch and
checks both closed forms interactively (runs in ~1 min on a laptop); 02
verifies artifact hashes then regenerates and displays all seven manuscript
figures. Both execute end-to-end (`make notebooks`).

Also pins `.python-version` to 3.14: the interpreter was previously
unpinned, and a `uv sync` silently switched it to 3.11 mid-session (see
tasks/lessons.md). README rewritten - it still described a Layer-1-only repo
of E1-E6 and pointed at a figure_manifest.md that does not exist.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01BkRLcc18rwT2Lysu6PbG7v
2026-09-07 15:50:50 +01:00

489 lines
110 KiB
Text
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

{
"cells": [
{
"cell_type": "markdown",
"id": "26f42fa9",
"metadata": {},
"source": [
"# 1 — The biological model, from scratch\n",
"\n",
"This notebook builds the paper's biological model in the open and checks it against the three\n",
"closed forms it is required to reproduce. Everything here is pure NumPy/SciPy: it runs on a\n",
"laptop in under a minute and needs no GPU.\n",
"\n",
"The claim being demonstrated is the one the whole paper rests on: **knowledge transmission\n",
"between model generations is a WrightFisher process**, not merely something that resembles one.\n",
"\n",
"Companion: `REPRODUCING.md` (the figure-by-figure map) and `tests/test_scientific_validation.py`\n",
"(the same checks as assertions with publication tolerances)."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "4e02e4f2",
"metadata": {
"execution": {
"iopub.execute_input": "2026-09-07T14:39:33.085488Z",
"iopub.status.busy": "2026-09-07T14:39:33.085304Z",
"iopub.status.idle": "2026-09-07T14:39:33.637432Z",
"shell.execute_reply": "2026-09-07T14:39:33.636220Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ready\n"
]
}
],
"source": [
"import sys, pathlib\n",
"sys.path.insert(0, str(pathlib.Path.cwd().parent / 'src'))\n",
"\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from knowledge.lineage import run_lineage\n",
"from knowledge.metrics import heterozygosity\n",
"from knowledge.truth import make_true_distribution\n",
"\n",
"K, N_SAMPLES = 200, 100 # K items of knowledge; n samples drawn per generation\n",
"print('ready')"
]
},
{
"cell_type": "markdown",
"id": "eca5c488",
"metadata": {},
"source": [
"## Knowledge as an allele-frequency distribution\n",
"\n",
"A model's knowledge is a distribution `p` over `K` discrete items (capabilities, facts, modes of\n",
"behaviour). Reality is a fixed distribution `p*` with a Zipf tail: a few things are common, most\n",
"are rare. In population-genetics terms `p` is a vector of allele frequencies at one locus with\n",
"`K` alleles, and the rare tail is what drift destroys first.\n",
"\n",
"Diversity is **heterozygosity**, `H = 1 Σ pᵢ²` — the probability that two random draws differ."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "6afb0adf",
"metadata": {
"execution": {
"iopub.execute_input": "2026-09-07T14:39:33.640226Z",
"iopub.status.busy": "2026-09-07T14:39:33.639948Z",
"iopub.status.idle": "2026-09-07T14:39:34.167745Z",
"shell.execute_reply": "2026-09-07T14:39:34.166324Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"H* (diversity of reality) = 0.9326\n",
"tail items (rare knowledge) = 70 of 200\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAE+CAYAAACnauueAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlcelbwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAN8lJREFUeJzt3Xl0jPf+B/D3TCKRRRZJMCFkIQhBEFtQrSIk1lJupUXt0qrbjbi9rVaJLqf106a0qlr7ri2KVBAhqboi1mokhCCWRPZEJJPn90du5mZkm32embxf5+Qc88yzfGY5ydt3eySCIAggIiIiEhGpsQsgIiIiehoDChEREYkOAwoRERGJDgMKERERiQ4DChEREYkOAwoRERGJDgMKERERiQ4DChEREYkOAwqRGUhISMDZs2fr3aYvcXFxSEpKMsi19MlcXoc+Gfu7Rg2HhCvJEpm+nj17wtXVFYcOHapz26lTp2BjY4Pu3bvr9Ppt27ZFt27dsGvXLp2e19DM5XXURhefv6rfNSJtWRq7ACLSj379+sHBwUFp2+zZs+Hp6Yn9+/cbqSoyJl18/jV9r4j0gQGFyEytWrXK2CWQGeL3igyFAYVIT1JSUpCSkoKhQ4cCAJKSkpCdnY3Bgwcr9snNzcWlS5cgl8vh7+8PZ2dnpXMkJibiwYMHAAALCwu4uLigc+fOsLKyqvf6CQkJsLKyQo8ePQAAR44cQUFBAR4+fKhoirezs8OAAQMQFxcHKysr9O7du9p5rl27htTUVDz33HMqXTc7OxsXL15E06ZN0blz52rPq/OacnJycO3aNZSVlaFDhw7V3p+aaPOeqfM6KuniM1T1/T99+rTGn1Ndn7+qdQLVv1dEeiMQkV588MEHAgAhKSlJCAgIEHr16iXIZDJBEAShuLhYmDt3rmBtbS107txZ6NGjh2BjYyMsXLhQkMvlinO8//77wrBhw4Rhw4YJgwYNEtzd3YWmTZsKO3fuVLpWjx49hGHDhtW5bcyYMYK9vb3g5uamOOe0adMEQRCEefPmCY0bNxaysrKqvY6ePXsKHTt2rPO1+vj4CC+88ILw008/Cb6+vkKfPn0Ea2trYeDAgUJBQYHSvqq+pn/961+CjY2N0K1bN6F///5Cs2bNhKlTpwqFhYV11qLq+bV9Hbr8DGfMmCHY2toKubm51Wrq16+f4v3X5nOq6/NX531T5btGpAsMKER6UhlQxo4dK9y+fVsQBEG4dOmSIAiCMGnSJKFJkybCqVOnFPufPHlSsLa2Fj799NNazymXy4WIiAjB1tZWuHnzpmK7qn80OnXqJISEhFQ776VLlwQAwueff660/cyZMzVuf5qPj4/QoUMH4e2331b8cT59+rQglUqFyMjIOo+t6TXFxsYKAIS9e/cq9isvLxd+/PFH4cGDB3WeT5Xz6+J16PIz/OOPPwQAwpo1a5T2vXr1qtL7r+3nVNvnr2qdgsCAQobDacZEevbiiy+iZcuWAIBOnTrh8uXL2LZtGxYuXIh+/fop9gsKCsIrr7yCL774Qul4uVyOS5cuISYmBtHR0fDz80NRURHi4uJ0VmOnTp3wzDPPYM2aNRCqTOz75ptv0KhRI7z88sv1niMrKwtLly6FVFrxa6VXr14ICgrCr7/+Wm3f+l5TWloaAMDHx0dxjEQiwZQpU+Dm5lZvLdq8Z6q8Dl1/hr1790bnzp3xww8/KB23bt06NGrUCK+88goA3XxOdTHEd41IVRyDQqRnffr0UXp88uRJAEDjxo1x5MgRCBUtmQAAS0tL3Lt3D/fu3UOLFi2wc+dOzJ8/H6WlpejQoQPs7e1RVlYGALh9+7ZO6wwPD8eLL76I6OhoDBs2DNnZ2di2bRtGjhyJZs2a1Xu8v78/GjdurLStTZs2OH78uNI2VV7T888/D2dnZzz33HMICwvDc889h2eeeUal2SPavmeqvA59fIbTp0/HP//5T1y6dAmdO3dGWVkZNmzYgFGjRimFMm0/p9oY8rtGpAoGFCI9a968udLjwsJCAMDPP/+M33//vdr+w4YNgyAISEtLw+TJkzF16lSsXr0aFhYWAICbN2/C09NT6X/QujB27FjIZDKsXr0aw4YNw48//oji4mJMnz5dpeNrGsBqbW2N4uJixWNVX5O7uzsuXbqE1atX48iRI/jmm29QXl6Ol156Cd9++221AKHu+bV9Hfr4DF9++WUsWrQI69atw5dffon9+/fj/v371d5/bT+nmhj6u0akCgYUIj2TSCRKj729vQEAixYtQkhISK3Hbd68GaWlpZg9e7biDwYAXLlyRWe1VGVpaYlZs2bh448/xq1bt7BmzRq0bNkSw4YN0/h6Tzt16pTKr8nd3R1Lly7F0qVLUVRUhG+//RZvvvkm/Pz8sHDhQq3Prw19fIYuLi4YPXo0Nm3ahE8++QQ//PADWrVqpZgFVkmbz6m2z99Q7xuROjgGhcjAhg4dCplMhi+++KLG/5lWTvWsbNbPyMhQPFdeXo6vvvpK42s3bdoU+fn5tT4/a9YsSCQShIWFITk5GVOnTlX6g6UtVV/Tw4cPlR7b2tpi3rx5kEqldXY36OM9q4m+PsMZM2YgMzMT3333HQ4ePFjr+6/p51Tb52+o941IHWxBITIwW1tb7N69G6NGjUKvXr0wbdo0yGQy3Lp1C8eOHYNcLse+ffswaNAgdOnSBfPmzcPdu3dhY2ODjRs3YsiQITh48KBG13722WexYsUKrF69Gl5eXkrrYAAVrRZjxozBrl27IJFIMG3aNF29bABQ+TXt3bsXX3/9NSZOnIj27dujtLQUGzduhK2tbZ1dGfp4z2qir8/w+eefh6enJ95++23I5fJa339NP6faPn9DvW9E6mBAIdKTtm3bYtiwYTX+z7Zv3764evUqNmzYgJMnT6K4uBienp6YOXMmRowYAQCwsrLCiRMnsGrVKkRHR6NJkyZ499130bdvX8TExCi6GYCalx+vaVtERATs7e1x9OhRFBQUQCaTKQUUAJg5cyZ27dqFQYMGKc2iqcvAgQNr3Ldz585K/2NX9TXNmjULAwcOxJYtW7B9+3ZYWVkhKCgI33//Pdzd3WutQ533TJvXAej+MwQqumD+9a9/YdeuXfD19a2zXk0+p9o+f31814i0xZsFEpGSb775BuHh4di0aRMmT55s7HKoFvycyNwxoBCRkoEDB+Kvv/5Cenp6rbNlyPj4OZG5YxcPEaGsrAxHjhxBQkIC4uLi8NVXX/GPngjxc6KGhC0oRITCwkK88MILsLOzQ3BwMGbOnGnskqgG/JyoIWFAISIiItHhOihEREQkOgwoREREJDocJKum8vJy3L17F02aNKlz2XAiIiJSJggC8vPz4e7urrhjeG0YUNR09+5deHh4GLsMIiIik5Weno5WrVrVuQ8DipqaNGkCoOLN5cqJREREqsvLy4OHh4fib2ldGFDUVNmt4+DgwIBCRESkAVWGSHCQLBEREYkOAwoRERGJDgMKERERiQ4DChEREYkOAwoRERGJDgOKCGTkFiM+NRMZucXGLoWIiEgUOM3YyLafuYWIPRdRLgBSCRA5zh8TA1sbuywiIiKjYguKEWXkFivCCQCUC8DiPZfYkkJERA0eA4oR3cgsVISTSnJBQFpmkXEKIiIiEgkGFCPycrWD9KnF9CwkEni62hqnICIiIpFgQFFRVFQU/Pz8EBgYqLNzyhxtEDnOHxb/XfLXQiLB8nGdIXO00dk1iIiITJFEEASh/t2oUl5eHhwdHZGbm6uze/Fk5BYjLbMInq62DCdERGS21Pkbylk8IiBztGEwISIiqoJdPERERCQ6DChEREQkOgwoREREJDoMKERERCQ6DChEREQkOgwoREREJDoMKERERCQ6DChEREQkOgwoREREJDoMKERERCQ6DChEREQkOgwoREREJDoMKERERCQ6DChEREQkOgwoREREJDoNMqAcO3YMkyZNwqRJk5CQkGDscoiIiOgpDTKgtGnTBmPGjEFGRgZu3rxp7HKIiIjoKZbGLsAYvL294e3tjUOHDhm7FCIiIqqB6ALKvXv3sG7dOmzYsAGWlpa4fPlytX3u37+PRYsWITY2FnZ2dpg8eTLeffddSKUVDULr16/H4cOHqx3Xp08fLFiwQN8vgYiIiLQkuoASEhKCoUOHIiQkBBs2bKj2vFwuR3BwMJydnbF7925kZGQgLCwMhYWFWLp0KQCga9eusLGxqXash4eH3uun6jJyi3EjsxBernaQOVb/XIiIiJ4muoBy5swZSKVSrFy5ssbn9+3bh/Pnz+PWrVto1aoVAgICsGTJEkRERCAiIgK2trbo3r07unfvbtjCqUbbz9xCxJ6LKBcAqQSIHOePiYGtjV0WERGJnOgGyVZ209TmxIkT6NChA1q1aqXYNmzYMBQVFeHs2bMqXePKlSuYNGkSYmNjsWrVKkydOrXWfUtKSpCXl6f0Q6rJyC1WhBMAKBeAxXsuISO32LiFERGR6ImuBaU+d+7cQfPmzZW2VT6+e/euSudwdXXFmDFjMGbMGABAo0aNat03MjISH374oWbFNnA3MgsV4aSSXBCQllnErh4iIqqTyQUUQRBgaalcduVjuVyu0jmaNWuGSZMmqbRvREQE3nzzTcXjvLw8jmVRkZerHaQSKIUUC4kEnq62xiuKiIhMgui6eOrj5uaGzMxMpW0PHz4EUBE8dM3a2hoODg5KP6QamaMNIsf5w0IiAVARTpaP68zWEyIiqpfJtaD06tULa9euRU5ODpycnABUjEuxtLTkwFgRmhjYGgN93ZCWWQRPV1uGEyIiUonJtaC88MILcHV1xbvvvovHjx/j9u3bWL58Of7xj3+gadOmertuVFQU/Pz8EBgYqLdrmCuZow36+rgwnBARkcokgiAI9e9mOC+//DJ+//13FBUVoaCgQNFtExsbi/bt2wMAzp07h2nTpuHq1asoLy/H2LFjsW7dOtjb2+u9vry8PDg6OiI3N5fdPURERGpQ52+o6AJKTk4OHj9+XG27q6trtcGxeXl5sLKyQuPGjQ1VHgMKERGRhtT5Gyq6MSiV40pUwYBARERknkxuDIqxcAwKERGR4Yiui0fs2MVDRESkGXX+hrIFhYiIiESHAYWIiIhEhwGFiIiIRIcBRUUcJGu6MnKLEZ+aybsoExGZEA6SVRMHyZqW7WduIWLPRZQLgFQCRI7zx8TA1sYui4ioQeIgWSJUtJxUhhOg4q7Ki/dcYksKEZEJYEAhs3Ujs1ARTirJBQFpmUXGKYiIiFTGgEJmy8vVDlKJ8jYLiQSerrbGKYiIiFTGgKIiDpI1PTJHG0SO84eFpCKlWEgkWD6uM++qTERkAjhIVk0cJGt6MnKLkZZZBE9XW4YTIiIjMumbBRLpmszRhsGEiMjEsIuHiIiIRIcBhYiIiESHAYWIiIhEhwGFSAtcRp+ISD84SFZFUVFRiIqKglwuN3YpJBJcRp+ISH84zVhNnGZMQEXLSdCKo0or1VpIJDi56FnOGCIiqgXvxUOkZ1xGn4hIvxhQiDTAZfSJiPSLAYVIA1xGn4hIvzhIlkhDEwNbY6CvG5fRJyLSA5VbUK5cuaKYwXL16lWUlpbqrSgiUyFztEFfHxeGEyIiHVM5oCQnJ+P999/HtWvX8N5770EikdR/EBEREZEGVA4oZWVliImJQXBwMAICApCVlaXPukQnKioKfn5+CAwMNHYpZIa44BsRkTKVx6D06tULTk5OkMlk6NixI5ycnPRYlviEh4cjPDxcMYebSFe44BsRUXUqt6Bs27YNs2bNwqZNm7B+/XqUl5frsy6iBiEjt1gRTgCgXAAW77nElhQiavBUDijvvPMOxo0bBzs7O/z666+wsal5UODYsWNRUFAAANi/fz9u3bqlm0qJzBAXfCMiqpnKAaXqoNi6Bsh+8sknsLa2BgDEx8dj0qRJ8PLywujRo7Fs2TKUlJRoUS6ReeGCb0RENdP5Oii+vr6Kfy9fvhxAxQDbCxcu4PTp07zZHlEVlQu+Ld5zCXJB4IJvRET/pfXNAq9fv459+/YhJycH3t7eGD58OFxdXXVVn+jwZoGkDxm5xVzwjYjMnjp/Q7VqQYmPj8ezzz6Lpk2bQiaT4dq1aygvL0dERAQWL14MqZQr6ROpQuZow2BCRFSFVgli7dq1ePHFF3H79m0kJiYiLy8PmzZtwqZNm7BgwQIdlUhEteH6KURkrrRqQSktLcXo0aNhYWEBoGLw7NixYxEUFARfX1+8/fbbaN2a6zkQ6QPXTyEic6ZVC8ro0aMRExNTbXuzZs3Qo0cP/P3339qcXlS4kiyJCddPISJzp1VAKS0tRUxMDN544w3cv39fsf3Bgwe4fPky2rdvr3WBYhEeHo4rV67gzJkzxi6FiOunEJHZ0yqgJCUl4dGjR1i1ahVat26NLl26YPDgwfD19cXkyZPRokULXdVJRFVw/RQiMndaBZRPP/0UWVlZuHnzJnbu3Inx48fD0dERTZs2xZdffgl7e3v4+/sjLCzMrLp7iIytcv0Ui/8umsj1U4jI3OhkobbWrVujdevWGDVqlGJbXl4ekpKSFD/Z2dm6uBQR/dfEwNYY6Oum9vopGbnFuJFZCC9XOwYaIhItrRdqa2i4UBuZMs78ISJjUudvKFdSI2ogOPOHiEwJAwpRA8GZP0RkShhQiBoIzvwhIlOiVUA5fPgwHjx4oKtaiEiP1J35w2X0iciYtJrFs379eowcORIhISF49dVXMWLECMWy90QkPqrO/OFgWiIyNq1aULZs2YL9+/fD2toaEyZMgIeHBxYtWoTk5GRd1UdEOiZztEFfH5c6W044mJaIjE2rgCKVSjF06FBs27YNGRkZWLx4MaKjo9G+fXv0798fP/zwA4qKOACPyJRwMC0RiYHOBsk6OzvjlVdewdy5c+Hh4YE//vgDb731Flq3bo21a9fq6jJEpGccTEtEYqB1QBEEAbGxsZgyZQpkMhmWL1+OGTNmIC0tDffu3cOnn36KBQsWsNuHyERwMC0RiYHWg2SXLVuG9PR0jBo1Cnv37sWQIUMgkfzvv1+vvvoqduzYgZSUFPj6+mpdsLFERUUhKioKcrnc2KUQ6R0H0xKRsWm11P1bb70FmUyGKVOmwM3Nrdb9rly5AplMBmdnZ00vJRpc6p6oQkZuMYJWHFUar2IhkeDkomd5jx8iqpE6f0O1akEZOXIk2rVrV2M4OX78ONq1a4eWLVvCz89Pm8sQkQjVNZiWAYWItKXVGJTPP/8c586dq/W5pKQkbU5PRCKm6mBajlEhIk3oban7e/fuwcnJSV+nJyIjU2Uw7fYztxC04iheWnsaQSuOYvuZW8Yql4hMjEZjUBYuXIhTp07VOrbk/v37yMnJwY0bN2Bvb6+zYsWAY1CIlGXkFtc4mJZjVIjoaXofg9K1a1fY2NjgwYMH6Nq1q9LsHIlEAjc3N4waNcrswgkRVSdztKkxcHCMChFpQ6OA8tJLLwEAAgMD0alTJ3h6euqyJiIyA5VjVJ5uQalpjMqNzEJ4udoxuBCRglazeEJCQnRVBxGZmcoxKov3XIJcEGodo8J1VIioJmqPQZk8eTIOHz6MLVu24KeffsLhw4dr3XfLli0YOnSo1kWKCcegEKmHY1SIqJJex6DMmTMHI0eOROfOnRX/rk3nzp3VPT0RmRltxqiw+4eo4VI7oAwYMEDxb3d3d50WQ0QNR31jVNj9Q9SwqR1QCgsLUVpaqtK+9vb2sLTUapgLEZmpusaoZOQWK8IJUBFiFu+5hIG+bmxJIWog1E4PEydOxIEDB1Tad9++fQgNDVW7KCJqGGq7KWF93T/s+iEyf2oHlJUrV2LJkiUq7duuXTt1T09EDUxNY1Tq6v5h1w9Rw6B2QGnbtq0+6iAiUqit+wcAu36IGgi1A0pKSgpycnLQrl07PHz4EDk5ObXu265dOzg6OmpTHxE1UDV1/8SnZrLrh6iBUDugLFiwAAcOHMC+ffuwZs2aOsejcAwKEWnj6e4fdv0QNRxqL9RWOYvH3t4eJSUldc7oEfMsnoMHD+L06dPo1KkTxo0bBwsLC5WO40JtRMa1/cytal0/A33duOgbkQnQ60JtdnZ2/ztYpOGjPm+99RauXr2KXr16YcmSJUhISMAXX3xh7LKISAXqdv0AYLcPkQnSOmFcv34dkZGRSEhIQHZ2Njw9PTF58mTMnDkTjRo10kWNOjd79mzFHZiHDBmCN99808gVEZE6VO36uXAnB5O//4PdPkQmSKuAcuXKFfTp0wf+/v54+eWX4ejoiJSUFHzwwQc4dOgQfv31V7XPWVRUhK1bt2LDhg2QSqU4duxYtX3y8/OxfPlyxMbGws7ODpMnT8bUqVMVz//88884efJkteO6dOmCV155RRFOAGDbtm14+eWX1a6TiMSjplk/7wa3xycHr9Y644cDaonETauAsmrVKoSGhmLLli1K2+fPnw9/f3+cP38eXbt2VeucvXv3Rs+ePeHr64u9e/dWe14QBISGhiI/Px/Lli1DRkYGwsPD8ejRI0VLSJMmTdCiRYtqxzo7OyudZ9GiRbC3t0d4eLhaNRKR+Dzd9VPXYm8nkh9yQC2RyGkVUPLz82ucpdO6dWt07doVeXl5ap/zjz/+gJ2dHVauXFljQImOjsaJEyeQnJysWAjuwYMHWLp0KV577TVYWVlh8ODBGDx4cK3XePLkCaZPnw5/f3+8++67atdIROL0dNdPTd0+tlZSrqVCZAKk2hzco0cPREdHV9uekZGBa9euaXQ346qDcGsSExODdu3aKa1SGxoaipycHJw9e1ala8yYMQOnT5/GgwcP8Pbbb2P58uW17ltSUoK8vDylHyISv8puHwuJBAAUM34Kn8hrbFk5m5aN+NRMZOQWG6FaInqa2i0ocXFxuHPnDgDA1dUVe/fuxd27dxEaGgoHBwekpqZi/fr1CAwMRHFxsVK3ii6kp6dDJpMpbau8q3J6ejr69u1b7zlCQkLQpUsXxeO6FpOLjIzEhx9+qGG1RGRMNc34ycgtrtayIgEwf9s5dvkQiYjaAWXNmjU4fPiw4nGjRo2QmJiIxMREpf3i4+Nx6dIlRXjQldLSUlhbWyttq3ys6l2WJ06cqPL1IiIilGb55OXlwcPDQ+Xjici4nu72eXpArRSAANTY5QNwijKRsagdUDZv3qyPOlTm4uKCGzduKG3LyspSPKdr1tbW1QIREZm2qi0rWYUleG3LOaXn5YKA9SfT8P3J62xVITISrcagGEOPHj1w5coVFBQUKLYlJCRAKpUiICDAiJURkSmROdqgr48LerRxhlSi/JwUUIQT4H+tKufTOU6FyFC0XqgtNzcXGzduxLVr11BYWKj03Pz585XGeujC+PHjsWjRIixbtgyRkZEoKCjAp59+ipEjR6J58+Y6vVZVUVFRiIqKglwu19s1iMjwalpDZXp/T3wXp9xSKxcEjPkmHgJbVIgMQu178VT14MEDdO3aFe7u7sjIyICXlxeKi4tx7tw5dOvWDT/88IParRrz58/HiRMn8PDhQ9y/f18xE2j37t3w8fEBABw7dgxhYWGQy+XIz89Hz549sWvXLri5uWn6UlTGe/EQmaeM3GLFYFoA1e7t8zQLiQR75vVF4RM5x6gQqUidv6FaBZR///vfSE1NxZYtWxAaGoo5c+YgNDQUO3bswNtvv43U1FS1l7u/ceMGcnNzq23v0KEDGjdurHgsl8uRmpoKW1tbtGrVStOXoDYGFKKGoepNCaUAymvYRyIBW1SI1KDXmwVWdfXqVYwfPx5AxWye4uKKftkXX3wRn332GRITE9G7d2+1zunl5aXSfhYWFkpL1hMR6VLVgbS2VlKM/Sa+WouK8NQYlQ4tmrBFhUhHtAooxcXFioXVZDIZrl27pnhOLpcrAos54BgUooan6hTlp6cmP92iwjEqRLql9SDZSiNHjkRYWBikUinS0tKQmpqq8wGyxhQeHo7w8HBF8xQRNSxsUSEyLK3GoJw5cwZt2rRBs2bNAABfffUVNm/ejMaNG+ODDz7As88+q7NCxYJjUIgI4BgVIk0YbJBsQ8SAQkSVKmf+1NaiUhVn/RAZIaAcPnwYCQkJyM7OhqenJ8aPH2+2y8EzoBBRTdiiQlQ/gwWUgoIChIaG4tSpU+jWrRscHR2RmpqKjIwMrFy5EnPmzNH01KJTdZBscnIyAwoRVaNOi4oUwKqXAtCjjTNbU6jBMFhA+fDDD7Fjxw4cPHgQrVv/738CGzduxMyZM3H9+nWd3yzQ2NiCQkSqUKVFBahoTVkY3AH+rRzZ9UNmz2ABZcKECQgJCcHUqVOrPde7d28sX74cgwcP1vT0osSAQkSqUqdFBWDXD5k/df6GanWzQG9vbzx48KDa9vLycmRlZcHb21ub0xMRmbTKGxJ29XBG5Dh/WEgkde5fLgARuy9i/4W7vCEhNXhqr4NSXl6O8vKKxsoZM2Zg0KBBcHJyQmhoKBwcHJCamooVK1YgICAAnp6euq6XiMgkVa6jcjYtG/O3nau1NaUcwGtbzrHrhxo8tbt4QkNDceDAAZX23bdvH0JDQzUqTGw4SJaIdKXq+BRVsOuHzIVex6BcvHgRWVlZKu3r7+8PFxcXdU4vehyDQkS6UDk+5cKdHHx68O96wwpn/ZA54EJtesSAQkS6lpFbXG/XTyV2/ZApM9jdjIGKGwauXr0aCQkJePToEby8vBAWFoZBgwZpe2oiogZB5miD0K42KHxSVm/XT7kARB68CoBhhcybVi0ojx49Qt++ffHo0SMMHToUTk5OuHbtGmJiYrB48WIsXbpUl7WKAltQiEif1O36qcRxKmQKDNaC8vXXX8PV1RWJiYmws7NTbD9x4gSCg4Mxb948yGQybS5BRNSgyBxtFNOTR3V1V7nrp3KKcocWTdDVw9kwxRLpkVbroJw/fx5z5sxRCicAMHDgQAQEBODy5ctaFUdE1JBVdP24q7SGClAxRXnMN/H4NjYV8amZXEuFTJpWLShOTk5ITU2ttv3Jkye4desWnJyctDm9qFSdZkxEZEiVa6io0vUjVBmjIgEwc4AXpvX34vgUMjlajUGJiYnBiBEjsGTJEowaNQqOjo5ISUnB8uXLkZ6ejkuXLsHCwkKX9Rodx6AQkbFVHafyyW9Xa73PTyUOpiWxMOg0402bNuGdd97BvXv3FNuGDh2KNWvWwMvLS5tTixIDChGJyfn0bIz5Jh6q/iZnqwoZk8ECSnl5OSQSCeRyOdLS0pCdnQ0vLy+4urpqekrRY0AhIrFRd2VagLN+yDgMFlBGjRqFWbNmmc1y9qpgQCEiMVK32weomCWxN7wfZ/2QwRjsbsbNmjWr8W7GRERkWJVTk2cP9MGpiOcwa4A3pPVM/OGsHxIzrVpQEhISMHfuXOzbtw8eHh66rEu02IJCRKZC3VYVjk8hfTNYF88///lPrF27FiUlJfD29kaTJk2Unl+5ciX69++v6elFiQGFiExRRm4x1p9Mw9q466jvlz6DCumLwVaSDQoKQvPmzWt93pxWkeU6KERkymSONlgc0hEhXVrUO+tHAPBd3A2sjbvBoEJGw7sZq4ktKERk6tSd9cMWFdIVg66DAgCxsbFISEhAdnY2PD09MXbsWLRo0ULb04oSAwoRmQNNZv0wqJC2DBZQioqKMHbsWBw5cgR+fn5wdHTE9evXkZubi6+++gqvvvqqpqcWLQYUIjI36oxPARhUSHMGCyjLli3Dhg0bcODAAbRt2xZAxeJt69atw/z583Hjxg2za0lhQCEic6VJUJk3yAdB7Vy5hD6pxGABZcKECRg+fHiNLSW9evVCZGQkBg8erOnpRYkBhYjMnbpBBWCrCqnGYAu1eXp6Ijs7u9p2QRCQk5NjlvfiISIyd5UzfuL/u+BbPeu9AfjfzJ9+kUex/MAVLvpGWtOqBSUlJQWDBg3C8uXLERoaCgcHB6SmpiIyMhJ5eXnYs2ePLmsVBbagEFFDo2mLyqLhHTD7GR99lkYmxmBdPKNHj8a+fftQ0ymkUikkkv/l7l9//RUjRozQ9FKiwYBCRA1VRm4xvopJwZY/b6l8zEu9PfD6c+3Y7UMADBhQEhMTVb4XT48ePeDm5qbppUSDAYWIGrpvT6Tik4NXUa7iXw8JgH/08kC/tq7o0caZYaUBM/g6KA1B1ZVkk5OTGVCIqEGrXEfF1kqKAxfuqdX9M4uDaRssBhQ9YgsKEVF16nb/cIxKw2Swe/EQEREBFTN/lo/zRxtXW6z47Wq9rSkCgMiDV/FXRi6e92vBrh+qhi0oamILChFR3TSZ9QOw66chYBePHjGgEBGpRtOg8lIvD7w+mDN/zBEDih4xoBARqScjtxhn07IR89d97E26q/JxbFExPwwoesSAQkSkuW9PpKo0RqUqtqiYDwYUPWJAISLSDseoNFwMKHrEgEJEpBuadv2wRcV0MaDoEQMKEZHuadL1E8F1VEwOA4oeMaAQEemHJi0qY7rJuI6KCWFA0SMGFCIi/eNgWvPEgKJHDChERIah6RgVdv2IFwOKHjGgEBEZ3rcnUhH521WV92fXjzgxoOgRAwoRkXFU3pBw65+32PVjohhQ9CAqKgpRUVGQy+VITk5mQCEiMpLKdVS+i7uu1nEMKsbHgKJHbEEhIhKHyhaVLX/eUus4BhXjYUDRIwYUIiJxYdeP6WBA0SMGFCIicWLXj/gxoOgRAwoRkbix60e8GFD0iAGFiMg0MKiIDwOKHjGgEBGZFk2DCtdS0T0GFD1iQCEiMk2aBhWAq9PqCgOKHjGgEBGZNk2Dyuhu7lg0vANbU7TAgKJHDChEROaBXT+Gx4CiRwwoRETmhV0/hsOAokcMKERE5knToDLUrzlGdXNni4oKGFD0iAGFiMi8sUVFfxhQ9IgBhYioYdA0qPTzccHC4Pbo6uGsp8pMFwOKHjGgEBE1LJoGlW6tHDFjoDe7fqpgQNEjBhQiooaJXT/aY0DRIwYUIqKGjV0/mmNA0SMGFCIiAjQPKh2a2+OT8V0aZFBhQFHB8ePHcfjwYfj4+CAsLAyNGzdW6TgGFCIiqopBRXUMKPX48ccfsW3bNjz33HOIjo5Gs2bNsGXLFpWOZUAhIqKaZOQW45ODV/Fz0l21jmtIQYUBpR537txBy5YtAQApKSkYMWIEkpOTVTqWAYWIiOrCFpXamXRAycvLw8aNG7FhwwZIpVIkJCRU2ycnJwcffvghYmNjYWdnh8mTJ2POnDmK53fu3InY2NhqxwUEBGD69OkAgN9//x27d+/GhQsX8Oabb2L8+PEq18eAQkRE9dE0qPi42mJqkBee92tudtOTTTqgdOnSBUFBQQAqgkZmZqbS8+Xl5ejfvz/kcjlWrFiBjIwMzJ49G++99x4WLlwIADhx4gQuXLhQ7dxt27ZFcHAwACAxMREnTpzAqVOnYGlpia1bt6pUHwMKERGpQ9OuHwB4qZcHXh/czmyCikkHlJKSElhbW2PlypX4+OOPqwWU3377DSEhIUhNTYW3tzcA4LPPPsPHH3+M+/fvqzTYNTExEd27dwcAlJWVwcHBAdnZ2bC2tq73WAYUIiLShDZBJXyQD94J7qCHqgxLnb+hUgPVpLL6QsKxY8fg6+urCCcAEBISgry8PCQmJqp0jX379uGZZ57BzJkz0aNHD4wcObLW65aUlCAvL0/ph4iISF0yRxusnBSAhIjnMKabu1rHRh1PxTOfHsXGhDRk5BbrqUJxEV1AqU96ejpkMpnSthYtWgAAbt++rdI5PvjgA3z22Wfo27cvoqKisH379lr3jYyMhKOjo+LHw8ND8+KJiKjB0zSo3HxUjH//chl9I4/is0NX9VihOJhcQJHL5bCyslLaVtn6UVZWpvJ5evXqhVdffRX9+/evc7+IiAjk5uYqftLT09UvmoiI6CnatqgM+fI49l+4a7YtKiYXUFxcXKqNS6l87OLiovPrWVtbw8HBQemHiIhIVzQNKtfuF+K1LefQN/Io3th6zuyCiskFlMDAQFy5cgX5+fmKbfHx8ZBKpQgICNDbdaOiouDn54fAwEC9XYOIiBquqkHlpV6t1Tr2l/N3zS6oiG4WT6XaZvHk5OSgbdu2mDJlCj777DPk5uZi0KBB8PX1xc6dO/VeF2fxEBGRIWgz62d0V3csGtFBdNOTTXqa8dy5c3Hs2DFkZ2cjKysLbdu2BVAx86Zdu3YAgFOnTiEsLAy5ubkoLCzEoEGDsHXrVjRt2lTv9TGgEBGRIZlTUDHpgHL79m0UFBRU2+7t7a00OFYQBKSnp8PW1haurq4Gq48BhYiIjCEjtxhzN51FUnqu2seKJaiYdEAROwYUIiIypvPp2VgVk4KYqw/UPtbYQYUBRQ+ioqIQFRUFuVyO5ORkBhQiIjIqbbp+hnZsjg/HdDJ4UGFA0SO2oBARkZiY0hgVBhQ9YkAhIiIxMoUWFQYUPWJAISIiMRNziwoDih4xoBARkSkQY1BhQNEDDpIlIiJTJKauHwYUPWILChERmaKM3GIs+eUyDl+5r/axn7zgj4mB6i2/XxN1/oaa3L14iIiISH0yRxt8+0pPje6evGj3RYPf44cBhYiIqAGpelPCYX7NVTpGAHA2LVu/hT2FAYWIiKgBqtqi0s3Dsd79JRIDFFUFA4qKoqKi4Ofnh8DAQGOXQkREpDMyRxv8HN4fv4T3Q38flxr3kUiA7m2cDVoXB8mqiYNkiYjInD09mFYqASLHGX6QrKXWVyMiIiKzUdn1k5FbjLTMIni62hrl5oIMKERERFSNzNHGaHc9BjgGhYiIiESIAYWIiIhEhwGFiIiIRIcBRUWcZkxERGQ4nGasptzcXDg5OSE9PZ3TjImIiNSQl5cHDw8P5OTkwNGx7sXhOItHTfn5+QAADw8PI1dCRERkmvLz8+sNKGxBUVN5eTnu3r2LJk2aQPLfdX8DAwNx5syZOo+ra5/KRGnOrTKqvEemen1dnVvT82hynKrH6GI/c/9+G/u7re8aTO37rc7+2v7uBvj9VpcgCMjPz4e7uzuk0rpHmbAFRU1SqRStWrVS2mZhYVHvF1OVfRwcHMzyCw6o9vpN9fq6Orem59HkOFWP0eV+5vr9NvZ3W981mNr3W539dfW7G+D3Wx31tZxU4iBZHQgPD9fJPubM2K9fn9fX1bk1PY8mx6l6jK73M0dieO38fmu2P39318+Yr59dPCLA+/uQOeP3m8wZv9/6wxYUEbC2tsYHH3wAa2trY5dCpHP8fpM54/dbf9iCQkRERKLDFhQiIiISHQYUIiIiEh1OMyYio8nJyUFmZiYAoGXLlrCxMd6t3YlIXNiCYgLef/99uLi4wMPDAzt27DB2OUQ6s337dgQHB6Nr1644deqUscsh0qlvv/0WHh4ecHBwwJQpU/DkyRNjl2RSGFBE7vjx49ixYwcuXryIX3/9Fa+//jqysrKMXRaRTsyePRspKSl49tlnjV0KkU5lZWUhKSkJCQkJSElJwV9//YXNmzcbuyyTwoBiAGlpaTh06JCiKftpgiAgKSkJsbGxyM3NVXru999/xyuvvAJ3d3cEBARgwIABiIuLM0TZRPV6/Pgxdu7cie+++67WfdLT07Fp0ybs2LGD4ZpMSklJCY4ePYrExMRa97l37x6OHz+Ov//+W2m7i4sLVq9ejVatWqFZs2bw9vaGi4uLvks2KwwoevTnn38iJCQE/fv3x/Dhw/Gf//yn2j4ZGRno3r07hg8fjgULFqBly5ZKKTszMxNubm6Kx82aNas16BAZ0scffwwfHx989NFHmDdvXo37bN68Ge3bt8fWrVvx9ddfo23btjh58qSBKyVST3FxMRYuXAhvb29MmDAB77//fo37ffjhh/Dy8sK7776L3r17IzQ0FI8fP66238aNGwEAI0eO1Gvd5oYBRY9u3ryJuXPn1pm+58yZA2tra9y4cQPnzp3DJ598gunTp+PmzZsAADc3N9y/f1+x//3795UCC5GxeHl54fz583jrrbdqfD4zMxOzZ8/GsmXLcODAAZw4cQLjxo3DlClTUF5ebuBqiVSXn58PZ2dnnDt3DkOGDKlxnyNHjuCjjz7C4cOH8eeff+Kvv/5CYmIili1bprTfl19+iSNHjmDTpk2KG8ySahhQ9GjChAkIDQ2t9Y6NWVlZ2L9/PxYsWIDGjRsDqOiTt7Ozw7Zt2wAAQ4cOxU8//YQbN27gjz/+QFxcHAYMGGCw10BUm8mTJ8PV1bXW5/fv34/S0lLMmDFDse3111/H9evXFXdHLSwsREpKCoqKinD37l1FMCcypmbNmmHRokVo1qxZrfv8+OOP6Nu3LwYOHAgAkMlkmDp1Kn788UcAFXe+f+ONN5CcnIzvv/9esY1Ux4BiRBcvXkR5eTkCAgIU2ywtLdG5c2ecP38eADBw4EBMmTIFQUFBCAsLw5o1a9C0aVNjlUykssuXL8PDwwNNmjRRbPPz81M8BwBxcXEIDg7GrVu38NFHH2HKlClGqZVIXUlJSUq/uwEgICAAt2/fxqNHj5Ceno6oqCisXbsWNjY2aNy4MT7++GMjVWuauA6KEeXk5ABAtcDh4uKC7OxsxeP33nsP7733niFLI9JaXl4enJyclLZZWVnB1tYWeXl5AIDg4GCkpKQYoToi7eTk5NT4uxsAsrOz4ePjg7KyMmOUZjbYgmJEVlZWACoGZFVVVFSkeI7IVNnY2CA/P19pm1wuR3FxMWxtbY1UFZFuWFlZ1fi7u/I50h4DihF5enoCqJiGWdXt27cVzxGZqnbt2uHOnTsoLS1VbLtx4wYEQUDbtm2NWBmR9jw9PWv83W1lZQWZTGakqswLA4oRdezYEa1bt8bPP/+s2Pb333/j8uXLCA4ONl5hRDowYsQIPH78GL/88oti26ZNm9C0aVMEBQUZsTIi7QUHByM6OlrRagIAu3fvxvPPPw9LS46e0AW+i3p07949JCUlKfrbz549CwDw9vaGr68vJBIJPvvsM0yePBn29vbw8vLCihUrMHjwYAYUEr2jR48iOTkZCQkJEAQBa9asAQCMGTMGLVq0gJeXFyIiIjB9+nScO3cOBQUFWL16NdatWwdra2sjV09Ut6NHj+LJkye4d+8eiouLcejQIVhaWuL5558HAMydOxdr165FaGgoZs6ciePHj+PUqVNc50eHJIIgCMYuwlzFxsYiMjKy2vaJEydi2rRpiscxMTHYsGED8vLy0LdvX8yfP18x7ZhIrNavX4/Tp09X2/7OO+/Ax8dH8fjw4cOIjo6GlZUVxo0bh8DAQEOWSaSRiRMnVlvZ28bGBnv37lU8fvjwIb788ktcvnwZMpkM4eHh8Pf3N3SpZosBhYiIiESHY1CIiIhIdBhQiIiISHQYUIiIiEh0GFCIiIhIdBhQiIiISHQYUIiIiEh0GFCIiIhIdBhQiEinfvnlF9y8edPYZRjM7t27cefOHWOXQWR2GFCISKdmz56NuLg4xeO9e/fi1q1bRqxIv6ZMmYIzZ84Yuwwis8OAQkQ6NWbMGKW7cc+cORPx8fHGK4iITBJvFkhEOjV8+HB4eHgAAH777Tc8efIECQkJAKC4Hw8AlJSUICEhAQUFBejUqRO8vLyUzrN792706dMHAHDhwgU4ODigX79+kEgkuH79Oi5evIg2bdqgW7duddZTeZ6ysjKcO3cO3t7e6NKlCw4fPozs7GxIJBK4u7sjICAA9vb2NR4rlUpx/vx5ODo6onfv3pBKa/+/XWFhIfbv3w9/f3/4+fmp9d4R0f8woBCRTs2ePRuff/452rRpg5iYGJSUlODMmTO4f/8+7OzsMG7cOPznP//B2LFj0bx5c8hkMsTHxyMsLAz/93//pzjPlClT0L17d9y5cwcdO3bEyZMnERQUhJ49e2LLli1o3749YmNj8dprr9V4U86q5+nfvz/+/vtvBAQEYMKECejSpQuOHTuGtLQ0lJeXIyUlBRkZGfj555/Ru3dvpWMHDBiAa9euoWPHjjhz5gw6d+6M6OjoGkPKo0ePEBISAicnJ+zevVu3byxRQyMQEelQ8+bNhY0bNyoeu7i4CFu3blU8Li4uFlq2bCmsXr1ase3u3buCm5ubsHfvXsU2Ozs7YeDAgUJJSYkgCIIQFxcnABCGDBkiPHnyRBAEQTh48KBgaWkpZGVl1VqPnZ2d0L17d6GgoKDOuv/9738LPXv2rHZs//79heLiYkEQBOH27duCtbW18Ntvvynts3fvXuH27dtCp06dhEmTJinqIyLNcQwKERlUdHQ07t+/D2dnZ+zatQs7d+5EXFwcvL29cezYMaV9p06dCisrKwBAnz59IJFIMG3aNDRq1AgA0LdvX5SVleH69et1XnP69Omws7Ortj0tLQ2HDh3C9u3b0bhxY5w7dw6lpaVK+0ybNg2NGzcGALRs2RLt2rXD33//rbRPcnIygoKC8Mwzz2Dz5s2K+ohIc+ziISKDSktLQ6NGjbB3716l7Z6envDx8VHa5uzsrPi3paUlpFKp0jZra2sAwOPHj+u8pkwmq7Zt9uzZ2Lx5M3r16gUXFxcUFBRALpfj0aNHaN68uWK/pk2bKh1nbW1d7XpLlixBq1atsHLlyjrHpxCR6hhQiMigHBwcUFpaip9++kkRMPRNIpEoPY6Pj8cPP/yA1NRUtG7dGgBw/PhxHDp0CIIgqH3+L774AitXrsSkSZOwfft2WFryVyuRthj1iUiv7O3tlVocBg8eDABYt26d0n5lZWV48OCBQWq6d+8e7O3t0apVK8W2Xbt2aXy+Fi1a4NixY7h8+TImTZqEsrIyXZRJ1KAx5hORXvXs2RPfffcdLCwsFLN4Pv/8cyxYsAAXL15Ez549kZ6ejl27duHLL7/EkCFD9F7TgAEDIJVKMWnSJAQHB+PUqVPVupzUJZPJcOzYMQwaNAj/+Mc/sHXrVrakEGmBLShEpFNPL9S2Zs0aBAcH4/fff8eBAwcAAG+88Qbi4+PRpEkTxMXFwcLCAr/88otSOBk/frxSCwcATJw4UWk8iYWFBSZOnAg3N7da66npPG5ubjh9+jTatGmD2NhYtG3bFseOHcPEiRNhY2NT57HDhg1TWt+k6j6VIcXGxgZ79uyp760iojpIBE06XImIiIj0iC0oREREJDoMKERERCQ6DChEREQkOgwoREREJDoMKERERCQ6DChEREQkOgwoREREJDoMKERERCQ6DChEREQkOgwoREREJDoMKERERCQ6DChEREQkOv8PmawTh7o1sh0AAAAASUVORK5CYII=",
"text/plain": [
"<Figure size 600x300 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"td = make_true_distribution(K, 1, 'zipf', 0.5, 1.1, 0, tail_threshold=1e-3)\n",
"p_star = td.p_star\n",
"H_star = heterozygosity(p_star)\n",
"\n",
"print(f'H* (diversity of reality) = {H_star:.4f}')\n",
"print(f'tail items (rare knowledge) = {int(td.tail_mask.sum())} of {K}')\n",
"\n",
"fig, ax = plt.subplots(figsize=(6, 3))\n",
"ax.loglog(np.arange(1, K + 1), np.sort(p_star)[::-1], '.')\n",
"ax.set(xlabel='item rank', ylabel='probability $p^*_i$', title='reality has a heavy tail')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "045375e5",
"metadata": {},
"source": [
"## One generation is a WrightFisher step\n",
"\n",
"A generation is: *draw `n` samples from the parent's distribution, optionally mix in `m` verified\n",
"real samples, refit the child*. With `m = 0` that resampling step **is** the WrightFisher\n",
"process — the textbook model of neutral evolution in a finite population of size `n`.\n",
"\n",
"The config below is the same schema the experiments use, so anything you learn here transfers\n",
"directly to `configs/layer1/*.yaml`."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "79caa8e6",
"metadata": {
"execution": {
"iopub.execute_input": "2026-09-07T14:39:34.170068Z",
"iopub.status.busy": "2026-09-07T14:39:34.169950Z",
"iopub.status.idle": "2026-09-07T14:39:34.203294Z",
"shell.execute_reply": "2026-09-07T14:39:34.202249Z"
}
},
"outputs": [
{
"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>generation</th>\n",
" <th>heterozygosity</th>\n",
" <th>tail_mass</th>\n",
" <th>support_size</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>0</td>\n",
" <td>0.932557</td>\n",
" <td>0.054961</td>\n",
" <td>200</td>\n",
" </tr>\n",
" <tr>\n",
" <th>50</th>\n",
" <td>50</td>\n",
" <td>0.614800</td>\n",
" <td>0.000000</td>\n",
" <td>4</td>\n",
" </tr>\n",
" <tr>\n",
" <th>100</th>\n",
" <td>100</td>\n",
" <td>0.498200</td>\n",
" <td>0.000000</td>\n",
" <td>2</td>\n",
" </tr>\n",
" <tr>\n",
" <th>150</th>\n",
" <td>150</td>\n",
" <td>0.495000</td>\n",
" <td>0.000000</td>\n",
" <td>2</td>\n",
" </tr>\n",
" <tr>\n",
" <th>200</th>\n",
" <td>200</td>\n",
" <td>0.130200</td>\n",
" <td>0.000000</td>\n",
" <td>2</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" generation heterozygosity tail_mass support_size\n",
"0 0 0.932557 0.054961 200\n",
"50 50 0.614800 0.000000 4\n",
"100 100 0.498200 0.000000 2\n",
"150 150 0.495000 0.000000 2\n",
"200 200 0.130200 0.000000 2"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"def cfg(m=0, generations=200, n=N_SAMPLES):\n",
" \"\"\"A lineage config: m = real samples mixed in per generation (m=0 is the dry null).\"\"\"\n",
" return {\n",
" 'truth': {'K': K, 'R': 1, 'tail': 'zipf', 'zipf_s': 1.1, 'tail_frac': 0.5,\n",
" 'tail_threshold': 1e-3, 'init': 'truth'},\n",
" 'dynamics': {'n': n,\n",
" 'teachers': {'K_T': 1, 'rho': 0.0, 'q': 1.0},\n",
" 'grounding': {'m': m, 'policy': 'proportional'},\n",
" 'selection': {'mode': 'none', 'novelty_alpha': 0.0},\n",
" 'remint': {'enabled': False, 'period': None, 'H_gate': None}},\n",
" 'metrics': {'kl_floor': 1e-9, 'support_eps': 1e-9},\n",
" 'generations': generations,\n",
" }\n",
"\n",
"one = run_lineage(cfg(), seed=1) # a run is a pure function of (config, seed)\n",
"one[['generation', 'heterozygosity', 'tail_mass', 'support_size']].iloc[::50]"
]
},
{
"cell_type": "markdown",
"id": "1135a629",
"metadata": {},
"source": [
"## Closed form 1 — neutral diversity decays geometrically\n",
"\n",
"WrightFisher predicts that diversity decays at a rate set only by the sample size:\n",
"\n",
"$$\\mathbb{E}[H_t] = H_0\\left(1 - \\frac{1}{n}\\right)^t$$\n",
"\n",
"This is an *expectation*, so we average lineages. Drift has high variance, and the agreement\n",
"tightens as $\\sqrt{R}$ in the number of replicates."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "0074c5ea",
"metadata": {
"execution": {
"iopub.execute_input": "2026-09-07T14:39:34.205125Z",
"iopub.status.busy": "2026-09-07T14:39:34.205000Z",
"iopub.status.idle": "2026-09-07T14:39:38.029646Z",
"shell.execute_reply": "2026-09-07T14:39:38.028989Z"
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAhgAAAFFCAYAAACufJSkAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlcelbwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAd1FJREFUeJzt3XdUVFf38PHvDFU6gkoRRRQ1KliwoMbee4mxx66PMTEmMTExxSQmPj6JiSUajbEbSzQae+8dRAUVsTdQEJHe633/8HV+GSkOCAzo/qw1aznnnnvvPjPMzPbeU1SKoigIIYQQQhQitb4DEEIIIcSrRxIMIYQQQhQ6STCEEEIIUegkwRBCCCFEoZMEQwghhBCFThIMIYQQQhQ6STCEEEIIUegkwRBCCCFEoTPUdwD6kJWVRWhoKJaWlqhUKn2HI4QQQpQaiqIQHx+Pk5MTanXu1yleywQjNDQUFxcXfYchhBBClFohISFUrFgx1+2vZYJhaWkJPH1xrKys9ByNEEIIUXrExcXh4uKi+S3NzWuZYDy7LWJlZSUJhhBCCFEAL+piIJ08hRBCCFHoJMEQQgghRKGTBEMIIYQQhe617IMhREmWmZlJenq6vsMQQrymDAwMMDQ0fOlpHCTBEKIESUhI4MGDByiKou9QhBCvMTMzMxwdHTE2Ni7wMSTBEKKEyMzM5MGDB5iZmVGuXDmZBE4IUewURSEtLY2IiAju3r2Lu7t7npNp5UUSDCFKiPT0dBRFoVy5cpQpU0bf4QghXlNlypTByMiI+/fvk5aWhqmpaYGOI508C0FaRha/HbnFoauP9B2KeAXIlQshhL4V9KqF1jEKIY7X3vJTd/l+2RZ6dWjJnfsP9B2OEEIIoXeSYBSCoU1cSDjyB4kPb9K+ex8ZASDE/zdr1iw++eSTUnWu4ow5KyuLPn364O/vXyzn06eNGzfSpUsXGjVqxLFjx/Qdjpbn3/Pi/BvQhzlz5jB79uwiP48kGIXAwtSYuYtXoDI2427gOcZ/8LG+QxKiRLh//z63bt0qseeaOXMmn3322Usfp6CWLVtGREQE9evXL5bz6cvNmzcZMmQIgwcPZtGiRdStWzdbnZSUFBYuXEifPn1o164dH3/8MQ8fPszxWCNHjqR58+YMGjQIPz+/AtX5t+ff8+L8G9CHt956i+nTpxMcHFyk55EEo5CM7tacVmO/BWD57wv466+/9BuQECXAlClT+OWXX/QdRq7u37/P7du3tcqKK2ZFUfjxxx95//33i/xc+nbx4kVsbGx45513aNiwITY2Ntnq9O7dm6CgIIYNG8bUqVO5efMmXl5ehIaGauo8fPiQpk2bkp6ezrfffkv58uVp2bIlAQEB+arzIiX97/ZlVapUiTZt2jB//vwiPY+MIikkKpWKFd9OoIG/H9Gn/2bEqNF4eHhQu3ZtfYcmRJEKDQ3lp59+4uLFi9jZ2TF48GD69u0LwIYNGwgPD+fnn38G4LvvviMtLQ0HBwcOHTpEbGwsb7/9NuPGjeO3335jx44dGBsbM3bsWPr06aM5x0cffYSLiwsff/x/VwfnzZvH7du3+fXXX3OMa/ny5SxcuBAAa2trGjRowOeff46dnR0AixYt4p9//iEjI4OGDRsCsHDhQo4dO6YVM8DevXv5448/CA8Px93dnU8//VTrs/3dd9+RkZFB1apV2bZtGwkJCfTo0YOJEyfm2mn3+PHjBAcH06tXL63jFOT1Abhw4YLmNXFxcWHMmDG0a9dO59ejoO140esze/ZsFixYQExMDA0bNsTExIRTp05lO8Y///yDmZmZ5nnLli1xcHBg/fr1TJ48GYBffvkFW1tbVq9ejVqtpkOHDvj7+/P999+zefNmneu8SE5/t7q8LoXxHgCsW7eOlStXYmhoSJs2bTA2Ns72t/6ic+X1uQTo168fH3/8MT/99FPRdSxX9OzJkyfKwoULlW+++UbZtGmTkpmZ+cJ9AgMDlXnz5infffedsm/fvnyfMzY2VgGU2NjYgoScp6VHbyimlesqgFKlarUiOYd4NSUnJytBQUFKcnKyoiiKkpWVpSSmpuvlkZWVpXPc9erVU/r27ascOnRI2bZtm9K/f39l7969iqIoynvvvaf06tVLU3f48OGKkZGRMnToUOXQoUPK/PnzFbVardSvX18ZOXKkcujQIeXnn39W1Gq1cvHiRc1+3bp1UyZNmqR13smTJyudOnXSPH/+XKGhoYqfn5/i5+en7Nu3T+nZs6dSp04dJSMjQ1EURXn48KHSt29fpW3btpp6cXFx2Y6zfft2xcjISJkxY4Zy8OBBZdSoUYqFhYVy9+5drXaZmJgoAwcOVA4cOKCsWrVKsbS0VBYvXpzr6/btt98qXl5eWmUFfX0OHz6s2NraKj/99JNy7Ngx5ffff1dsbGyUzZs36/x6FLQdL3p9goODlenTpyt2dnaKn5+fcv78+VyP9W+ZmZlK+fLllZkzZ2rKPD09lffff1+r3k8//aRYW1vnq87znn/Pc/q7fdHrUljvwbp16xQTExPll19+UQ4ePKgMHDhQMTEx0fpb1+VceX0uFUVRbt++rQDKpUuXcnxNnv8++jddf0P1egXj7t27NG/enOrVq9OoUSMmT57M8uXL2bFjR65DZObPn89nn33GkCFDsLe3Z+zYsbRr147ly5cXc/Q5G9miGns//B9bpr1DsqULcckZyIrwoiCS0zOpNW2fXs4dNL0TZsYv/npISUkhICCA5cuXa/oR9OzZk9TU1Fz3qVKlCqtWrUKtVtO2bVs2btxIbGwsy5YtQ6VS0bZtW9auXcvu3bvx9PQscBscHR1xdHTUPG/Tpg3lypXj+PHjtGnTBicnJ82EZs+uYORk6tSpvPfee3zxxRcAtGvXjsuXLzNjxgyWLFmiqVe5cmXWrl2r+e46e/Ys27dvZ9y4cTke986dO1SsWDFbeUFen48//pipU6fy6aefAk//95+QkMCMGTM0/2t90etR0Ha86PVxcXGhSpUqGBoa5vk6P++PP/4gKiqKHj16aMqCg4NxcnLSqufk5ERsbCxxcXFYWVnpVKcgXvS6FNZ7MH36dD7++GPN1bo2bdpQvXp1rVhedC5dPpfP/vZu376Nh4dHgV6TF9FrgvHZZ5/h6urKoUOHMDAw4N1336VGjRps2LCBQYMGZasfFxfH5MmTmTt3LhMmTADg/fffp2rVqgwaNIgOHToUdxOyUatVLBjVmmsP5xOntmDsX1f4c3Rj7C1M9B2aEIXO1NSU5s2bM2bMGN5//33atGmDq6srJia5/73Xr19f6z8QTk5OuLm5aV2mdXJyIjw8/KViS0pKYtGiRRw9epTHjx+TmZlJamoqd+7c0fpBzUtKSgpBQUFat0sAOnfuzM6dO7XKnm+Xs7Mzvr6+uR47NTU1x2mY8/v6xMbGEhAQQGpqKps3b0ZRFBRFITo6Wqv/gq6vR37akZ/XJz+OHj3Khx9+yE8//aR1Kyo9PT3b39azSemejd7TpU5B5PW6FNZ7kJyczPXr17VudajValq3bs2DBw90Ppcun0sjIyPUajUpKSkFfk1eRG8JRkZGBjt27OCXX37BwMAAADc3N1q1asU///yTY4Jx584d0tPT8fb21pQ5OztTqVIlNm3aVCISDIAKVqb8PbkrQ5ee5WpYHG8vOs1PnZ1p5FH9xTsL8f+VMTIgaHonvZ1bVwcOHGDFihVs3bqVjz76iJo1a7Ju3Trc3NxyrG9omP1rJ6cy5SXXYxkxYgS3bt3i008/pVKlSpiYmNC7d2+Sk5N1PkZCQgKKomBubq5VbmFhQXx8vFbZ821QqVR5tqFChQpcvnw5W3l+X5/ExEQAPvjgg2xXCP6dlOj6euSnHfl5fXR14sQJunfvztSpU/noo4+0ttnZ2REZGalV9uTJEwwNDbG2tta5TkHk9boU1nvw7PX8d18UQOv11fVcL/pcRkdHk5WVhYODQ/5fDB3pLcEIDg4mJSWFatWqaZVXq1aNM2fO5LhP1apVMTY25tixYzRo0AB42gv8/v37XLt2Lddzpaamal0aiouLK4QW5K2mgxV/j2/KoN+OcnbZ1zT77grHT/vS1FOSDKEblUql020KfStTpgwTJkxgwoQJJCUl0bFjR6ZPn87KlSsL7Rzm5uYkJSVplYWFheVaPzMzk+3bt7Nr1y7N/wYTExOJiIjQqqdWq/NMAuzt7bG2tiYoKIgWLVpoyq9cuULVqlUL0hSNRo0asWbNGhRFealOdhUqVMDKykrTiTInur4e+VXYr8+pU6fo2rUrkydP5ptvvsm2vUGDBtmGnPr6+uLh4aFJAHSpU9gK6z0oV64clpaWXL9+naZNm2rKr127pvkb0eVc8OLPpb+/f75vW+WX3oapPsvCnr8fZm1trdn2PEtLS3755Rc+//xz+vfvz4QJE2jXrh01a9bM9uXzbzNnzsTa2lrzcHFxKbyG5KGKvTlrxjVFHR9ORmIMHbr14lZYdLGcW4ji8OjRI37++WfNZVZjY2PUanWB1y7ITd26dTlw4ACxsbEAnD9/ni1btuRa38DAAGtra80PTWZmJpMnTyYtLU2rnoODg+bSc27GjBnDL7/8wqNHT5cCOHv2LBs3bsy1T4KuunTpQlJSEufPn3+p4xgYGDBhwgRmzZrF2bNnNeVXr17VjFjQ9fUoiMJ6fc6cOUOXLl34+OOP+e6773KsM3bsWI4cOcLRo0cBCAwM5O+//2bMmDH5qlPYCvM9GD58OL/88ovmKsyRI0c4ePBgvs6ly+fy4MGDtG/fHgsLi8J8KbToLcF41qhnXxjPxMTE5Nng999/n8uXL9OpUydq1arF/v37qVy5MmXLls11n6lTpxIbG6t5hISEFE4jdFDNsSyH92zH0MySxAfXaN3nHRJSZKZP8WqwtbUlPDwcBwcHateujYODA5mZmTn+7/NlvPfeezg6OlK5cmVq1KjByJEjad26dZ77zJ8/nx9++IHq1atToUIF7t27l61T5eDBgwkODsbNzY2GDRtqfWE/8+2331KnTh3c3NyoXr06rVq14sMPP9Qa8lcQdnZ2DB06tFA6qH///feMGDGC1q1b4+bmhpOTE/3799fqHKjL61EQhfX6jB49mpSUFHbt2kXDhg01jx9++EFTp2vXrnz33Xd06dKFqlWr4uXlxbBhw3j33XfzVacoFNZ7MH36dGxsbKhUqRLu7u6MGzeO9u3ba119edG5XvS5zMjIYM2aNXzwwQdF+probZhqRkaGYm5ursybN0+rvFWrVsrAgQN1Pk5aWppiZ2enfPPNNzrvU5TDVHOzbvM2BZVaAZQ3h01RMjN1HwYoXg95DQsr6VJTU5WgoCAlLCxMq/z+/fvKrVu3NM/v3r2r3LlzR6vO7du3lXv37mmV3bx5UwkODtYqy8rKUu7du6cpDwkJUW7cuJHruRRFURISEpTLly8roaGhiqI8HeL++PHjbLFfv35dM0w1p+MoiqI8evRIuXjxYo7fGzm1KywsTAkKCspW998ePnyo2NvbKyEhIbkeR9fXR1EUJSkpSbl8+XK29+GZF70eBW2HouT9+kRGRioBAQF57n/58mXNEM5/P55vu6I8/Q6/ePGi8uTJk1yPp0udZ55/z3X5u83tdXnZ9+CZW7duKffv31cURVH69OmjjB07Nt/nyu1z+ccffyht27bNcZ9nCmOYql7nwRg6dKji5eWlpKamKoqiKFeuXFEMDAy0xvJu3rxZmTNnjuZ5QECA1lwZM2fOVCwtLZUHDx7ofF59JBiKoiiTvvhOARRUamXy3DXFem5R8pXmBEMU3M2bN5WIiAh9hyFKiNu3b2vNV3Hw4EHF0NBQ2blzZ6Gd48aNGy/8myuMBEOlKC/ZVfslPHz4kBYtWmBra4uXlxfbt2+ndevWrF+/XtOhZcyYMfj4+BAYGAg8nQntt99+o3nz5ly/fp2zZ8+yfv16OnfurPN54+LisLa2JjY2tsBjogtCURTe7NKX0/u2YmBRlnOXgqhXpUKxnV+UbCkpKdy9e5cqVaoUeh8GIUTpEB8fz8CBA/H398fExITw8HC++eabbGvmFLW8vo90/Q3Va4IBT4flbN++nfDwcDw8PGjfvr3W9j179hAaGsro0aM1ZUFBQRw5cgQbGxu6dOmSZ/+LnOgrwYCn46BrNutIWp1e1GvYhK3vNcc0H0MCxatLEgwhxDORkZFERkZSpUoVjIyMiv38r0SCoQ/6TDAAniSk0nnucZ4kpDGgoQv/e8uj6OaCF6WGJBhCiJKiMBIMWU1VD+wtTPilfz3UKvhzzwmGTv7hxTsJIYQQpUjJn8XnFdWqejnGe1nz+eBPWJeegpOzEz99PEquZAghhHglyBUMPfr0reY0aN8bgF++mMjgWZuJSXr5iW+EEEIIfZMEQ49UKhUnt6zCvV4TlLRkNv3vA/rOOUBkQu4rUQohhBClgSQYemZqYszpAztxquhCRnQYZ5Z8xaDFp4hOlCsZQgghSi9JMEoAe3t7du/cgZmZOSn3AzizZhZjVvmRmfXaDfAR4rUwfPjwl14tVoiSThKMEqJu3bqsX7/uaSfPuMecuxPB0hN39B2WEKKQJSUlcerUKenQLV55kmCUID179mT//v38sWYTKkMjfjlwg5vh8foOS4hikZKSQs2aNTl37py+QykUmZmZ1KlTh+PHj2uV9+vXj6ioKLp378769ev1FJ0QRU8SjBKmffv2DPR2pU2NcqSmZzJ+6VHiZfVVUcKNHTs22yy8AL6+vpiamnL16tUXHmPOnDm4u7vTsGFDTdnNmzf55JNPqFChAq6uroUZspaLFy8yfvx4bG1tad68ea71tm/fTr169XQ6poGBAZ9++imTJk3SlGVlZVGjRg06dOjA+PHjadq06cuGLkSJJQlGCaRSqfi+5xskHVnMsR9HM+TXvaSkZ+o7LCFy5ePjg4eHR7bys2fPYmBgoLVkdU7S09P59ddfGT9+vFb58OHDqVChAqNGjSIlJaVQY34mNTWV4cOH4+npSd++fUlNzX0U1/bt23NMpHLTv39/7ty5w5EjRwBQq9UkJyczbNgwunfvXqRJkxD6JglGCWVhkIlF1HUy4yPYN3cy41eckU6fokRKTk7m6tWrWlcenjl37hz16tXDwCDv9XYOHz5MTEwMHTt21Co/ffo0n376KXZ2doUa87+ZmJgQEBDAhAkTsLS0zLWeoijs3r2b7t27A3D16lVMTU3Zu3cvnTt3pmzZsnh4eODr66vZp0yZMnTo0IF169Zpyo4fP07Lli2LrD1ClBSSYJRQNjY2HNi7GysbW9LCbvL37M+ZtffFl5mFKG4XL14kMzMz1wQjp/LnnThxAk9PT70s6qSr8+fPk5yczJtvvglAQEAA6enpzJs3j5kzZ3L9+nXc3Nz45JNPtPZr3Lgxx44d0zw3NzenS5cu/Prrr8UavxDFTaYKL8GqVavGjm1badeuPUnXTzHzh+nUqTib7p5O+g5NFKPExMRctxkYGGgtRJRXXbVaTZkyZV5Y19zcPF/xnT9/Hng6Eup5qampOi0zHRISgoODQ77OW9x27txJp06dMDR8+rV58eJF7O3t+euvv7C2tgagc+fO2RIHBwcHQkJCNM8PHjyIr68vlSpVKr7ghdADuYJRwrVs2ZKlS5cAEHdmA+O+nsPVsDg9RyWKk4WFRa6Pt956S6tu+fLlc63bpUsXrbqurq451suv8+fP06ZNG2JiYrQeW7ZsAcDLywt4eoth5syZ1K9fn+7du3Pr1i3NMRRFeelhm7///jumpqZ5Pv7+++8CH3/nzp2a2yPwNMHo06ePJrkAuHHjBu7u7lr7qVQqrTkvrK2t6dixIzVr1ixwLEKUBpJglALDhw/ns88/B+DR7vmMWHRAZvoUJcb58+fx8vLK9mMeFBSEmZmZ5od0w4YN/P333yxdupT27dtrJUfOzs6Eh4e/VBxjx47NluQ8/3g+IdNVWFgYFy9e1ErSLl68iLe3t1a9gICAbFdywsPDcXZ2LtB5hSjN5BZJKfHfGTMIDnnIRaOahKeZ8v76C6wa2RhDA8kRX3UJCQm5bnu+8+Tjx49zratWa/+t3Lt376XigqdzVwQFBTFlypRs2wICArQ6eG7dupVPPvkELy8vvLy8WLBgAbdv36Zq1ao0b96cefPmkZ6eXuB+GAYGBi/sTFpQu3btokmTJprOphEREYSFhVG/fn2tehcvXuS9997TKvPz89P02xDidSIJRimhVqtZt2Y11x/F02fhKU7dimT+4Vt81CHv4X+i9MtPn4iiqpubS5cukZGRke2HFp4mGG3bttU8f/ToES4uLprnlSpVIiwsjKpVq9KuXTvMzMw4fPgwnTp1eum4CltOt0eMjY2pVauWpuz+/ftER0drzZORlpbGoUOH+PPPP4szXCFKhBLz39+CjHHPa7z6q6qGgyUz+3qQHvOIrz8YzeFL9/QdkniNnT9/njJlylCjRg2t8tTUVK5du6bpfwFPR0bFxMRonkdHR2Nraws8Hc45YcIE/vjjD63j9OzZE1NTU6ZOnUp4eLjm9suNGzcKtR3PbvEsWLCACxcuaM6TmJhIamoqBw8ezJZg1KlTR+tqS0BAABYWFlStWlVTtnXrVsqXL0/nzp0LNV4hSgVFz6ZNm6ZYW1srBgYGSrVq1ZQ9e/bkWT80NFR56623FDMzM8XExEQpX768Mm3aNCUrK0vnc8bGxiqAEhsb+7Lh60VmZqZiX8ldARTral5KREyivkMShSA5OVkJCgpSkpOT9R2KztLT05WUlJQctyUnJyuZmZma5z/99JPSr18/JSsrSzl//rzi6OiopKamarYnJCQoFStWVPz9/TVlaWlpSnJycrZHfj7vukhNTc3xPIqiKHv27FEqV66sVT89PV1JS0vTKsvIyNBqT2ZmplKrVi1l7969hRqrEMUhr+8jXX9D9XqLZP78+cydO5ddu3bRuHFjfv75Z3r37k1gYCDVqlXLcZ+xY8cSFhbGrVu3cHR05MiRI3Tp0gUXFxfGjBlTzC3QD7Vazeb1q2nTug2xt87TpvcgLh3eKosniWJnaGioGbb5vH8PnwUYP348+/fvx9bWFpVKxYoVKzA2NtZsNzc3586dO1p9RYyMjIplbox/x/G852+PADm2+fk+IGq1mgsXLmBiYlJ4gQpRiuj1FsncuXMZM2YMb775JsbGxnzxxRc4ODjw+++/57rPlStX6NGjB46OjgC0adOGGjVqEBQUVFxhlwgtm3kze/FKUKkJPLqd0R9k72QnREliaWnJgQMHuHfvHk+ePKF3797Z6hgZGRVZR82CmjNnDvPmzSvQvpJciNeZ3hKMJ0+ecOfOHVq0aKFV3rJlS62pdp83duxY1q1bh6+vL48ePWL16tUEBwczZMiQog65xJk0cgBdx38FwIoFP/P7H0v0HJEQL2ZjY1Pikoi8lMSkR4jSQG+3SJ4Np7O3t9cqL1++fJ4JxieffMLly5fx9vbWXNZcuHChVmey56Wmpmp1CI2Le3Umqlo/+yveuHqb0KNrmfDuu1SqVImunUteL3whhBCvF72PIsnKytJ6npGRkWdfgrfffps7d+4QHBxMSkoK+/btY9KkSaxatSrXfWbOnIm1tbXm8e+hcqWdlakRfy2eg6VHO4wquLHmBrLyqhBCCL3TW4Lh5PR0PY3nZ+97/PixZtvzwsPD2b59O1988QUuLi6oVCpat27N22+/zeLFi3M919SpU4mNjdU8/r0uwKugRfVy/LNuFZWG/sjp0Aw+WO9Plqy8KoQQQo/0lmDY2NhQp04dDh06pCnLysri8OHDWrPeJScnEx8fD6BZqOn5OTOSk5O1FnF6nomJCVZWVlqPV037Os6s+k8LjA3U7A8KZ+CUnwgNDdV3WKIAFEWSQyGEfj1/d6Eg9DpM9YsvvmDEiBG0aNGCpk2bMmvWLFJTU3n33Xc1dSZOnIiPjw+BgYFYWVnRtWtXpk2bhrOzM25ubuzfv5/NmzezcOFCPbakZGhW1Z5Zb3sy6vOZ/L3vN3y2ruCi3xnNZEaiZDMyMkKlUhEREUG5cuVk2LEQotgpikJaWhoRERGo1eo8h3C/iF4TjEGDBpGcnMyPP/5IeHg4Hh4eHD58WDMEFcDMzEzrisPatWv5/vvvGT9+PJGRkVSuXJnff/+d0aNH66MJJU6ves68N7Qv/z21npDb12nbsQunjh3GzMxM36GJFzAwMKBixYo8ePCgUNYJEUKIgjIzM6NSpUrZ1jDKD5XyGl6PjYuLw9ramtjY2FfydklWlkK/mRvY9v0YslIT6dipCzt3bCuWCYvEy8vMzCQ9PV3fYQghXlMGBgYYGhrmehVV199QWezsFaRWq1jx0Vu0eBhB4JJP2b9vD8NHjGTNn6tfKhsVxaMoVwUVQojiIr82ryhrMyPWfjEMx7e+AJWa9evWMumjj6QDoRBCiGIhCcYrzKOiNRt/mIBDz48BFYfuJBEQEqPvsIQQQrwGpA/GK9gH43k+dyIZ8uNGMm0rAdCnvjM/vuWJsaHkl0IIIfJH199Q+YV5DXi72XH8xxH086qIWgWbfW4y8Js/ZDIuIYQQRUY6eb4mnG3K8PPbdWnjZslbPbuxJewGQ8lg+bfvYmokHQqFEEIULrmC8ZrpUt+VxnXfACWL9T9+xBujfuKbbYEkpmboOzQhhBCvEEkwXjNqtZrD2zfSuE0XyMzg/l/fsXjjLgYv9SUqMU3f4QkhhHhFSILxGjI0NOTE3q107twFJSOViM3TOXv2LG//fprgyCR9hyeEEOIVIAnGa8rY2Jh//tlM69atyUpNIuLvb7h6JZBu809wICj8xQcQQggh8iAJxmusTJkybN++HW9vb2wtylDT0ZL4lAzGrj7HzD1Xych8+dX0hBBCvJ5kHozXYB6MF4mOjubJkydUrlKV/+25xvJTdwFo7FqWBYPrU97KVM8RCiGEKClkHgyhM1tbW9zd3TE2VDOtRy3er5WBSWoMZ+9F0WfhaW6Ex+s7RCGEEKWMJBhCy5EjR5g2fhBp27/B2SiZhzHJvLXoNH73ovQdmhBCiFJEEgyhxc3NjfLly3Pvzm0eb/iCWjaZxKdk8N7aC8QmyxLiQgghdKNzgnHkyBEWLlyIr68vKSkpRRmT0KPKlStz5MgRKlWqxK2bN7m9YgrOxqk8jk/lx73X9B2eEEKIUkLnBOP+/ft88skneHt7Y2lpSd26dRk9erQkHa+gKlWqcOTIEVxcXLhx4zoP131OZkI063yD5VaJEEIIneRrFElmZiZXrlzh3Llz/Pnnn1y4cIHk5GTS09MxNDSkVq1afPrppwwdOrQoY35pMopEN3fu3KFVq1Y8ePAAO2c3yvT9nkpODmwc3xRnmzL6Dk8IIYQe6Pobmq/FzgwMDPD09MTT05OoqCi8vLyYMWMGV65cISAggICAgHxfyXj48CGrV68mPDwcDw8P3nnnHYyNjXOt/+GHH+Z4jjfffLPEJzaljZubG0ePHqVVq1Z4NapLlFM5QmKSGfSHDxv+442jtSQZQgghcvbSq6mamJjQoEEDGjRokO99r127RrNmzXjzzTdp0qQJs2bNYtWqVRw+fBhDw5xD8/DwID39/zobPnz4kB9++IH69esXuA0id1WrVuXMmTM4ODgQkZjBwD98CI5Kov/iM6wY0Yhq5S31HaIQQogSSK/LtX/22WfUrVuXbdu2oVKpGDlyJFWqVGHNmjWMGDEix31Gjx6t9fyHH37A3NycQYMGFUPErycXFxcAnGyMWDumMS0Gf8C96u3os/A0C4c0oIV7OT1HKIQQoqTRuZNnVFQUqamphXbi9PR09u7dy6BBg1CpVAA4OTnRpk0btm/frtMxFEVh+fLlDBgwQPpSFJNFP8/g/r7lxG3+mpjICEau8OPI9cf6DksIIUQJo3OCsXz5ciwtLWnQoAHjxo3j+PHjxMXFkZGRUaAT379/n7S0NKpUqaJVXqVKFW7duqXTMQ4fPszdu3cZO3ZsnvVSU1OJi4vTeoiCGTFiBI6OjsSF3iFh81ekxEYw/s/z+NyJ1HdoQgghShCdE4xBgwaxePFimjdvTmBgIIcOHWLJkiVYWFjQsGFDxo4dy8KFC7l586ZOx0tOTgbA0lL7Hr6VlRVJSbotGb5s2TLq1KmDt7d3nvVmzpyJtbW15vHskr/Iv+rVq3P8+HFcXFyICbtP3N9fkhj5iDGrznHpQYy+wxNCCFFC6JxgODs7M3LkSObPn8/p06eJi4vj8uXLmqTj6tWrTJkyhcWLF+t0vGe3NGJiYrTKo6OjdbrdER0dzZYtW1549QJg6tSpxMbGah4hISE6xShyVq1aNY4fP46bmxvxEQ+J3vgF0Y9CGLb8LNcfybolQgghXqKTp4GBAXXq1KFOnToMHz4cgKysLJ1vP7i4uGBpaUlQUBCdO3fWlAcFBVG7du0X7r927VoA3nnnnRfWNTExwcTERKe4hG5cXV05fvw4bdu25caNG7B5GoYjFjJ0mS+bxjelsp25vkMUQgihR4W6FolarcbGxkbnuv3792fFihWaWyLnz5/n9OnTDBw4UFNvzZo1fP/999n2X7ZsGW+99Ra2traFErvIP2dnZ44dO0a9evVY9Osc3qhYloj4VIYs9SUsNlnf4QkhhNCjfM3kWdgiIiJo06YN6enp1K1blwMHDtC/f3+t2yxjxozBx8eHwMBATdmFCxfw8vLSTAKVXzKTZ+HKyMjA0NCQx/Ep9P/9DHefJFC1vCXjW1alfiUb3CvIXBlCCPGq0PU3VK8JBkBaWhoHDx4kPDwcT09PvLy8tLYfO3aM8PBw+vfvrykLCAjg3LlzjBkzpkDnlASj6Phevk7bTt2wavcuJs41AehVz4mZfT0wM9brtCtCCCEKQalJMPRBEoyiM2zYMP7880+My5jR7N1ZBJu4kpmlUKOCJUuGNaSSnZm+QxRCCPESdP0NzXcfjH379rFx48ZCnXRLvDoWLVpE27ZtSUtO4sxvHzPBLQZ7CxOuh8czbLkvsUnpLz6IEEKIUi/fCUZ8fDyjRo3CycmJSZMmcfHixaKIS5RS5ubm7Nq1i549e5Kamspn44cxonwIzjZluBeZxPvrL5CRmaXvMIUQQhSxfCcY/fr1IywsjB9//BE/Pz/q1auHl5cXv/32G9HR0UURoyhlTE1N2bx5M8OGDSMzM5OJ/xlFy3Q/yhgZcOLmE77fGcRreGdOCCFeKwUapmppacmYMWM4ffo0V69epX379syYMQMnJycGDx7MsWPHCjtOUcoYGhqyYsUKPvjgAwB2b1rHzF5PO32uOnOf73delSRDCCFeYS89D4aVlRW2trZYWFigUql4+PAh7du3p23btsTGxhZGjKKUUqvVzJ07l19//ZV9+/bRu6ErM/t6ALD81F3eX+9P4EP5GxFCiFdRgRKM9PR0Nm/eTLdu3ahUqRIbN25k0qRJhIWFcezYMe7du0dycjLLli0r7HhFKaNSqZg4cSIODg4ADGpcicFOUZCVwa5LYXSff5IRK84SnyKdP4UQ4lWS7wRjx44dODk5MWrUKFxcXPDx8eHChQu89957WFtbA09neOzTpw9hYWGFHrAo3VavXs3MScOo7P87Xd8oi5GBiqPXIxi8xJeoxDR9hyeEEKKQ5HvmI5VKxU8//cSAAQMwM8t9ToNJkybJPXaRja2tLaamphw/uBfSkli1aDXvb77B5Yex9F98hjWjm+BgbarvMIUQQrykfF/BsLa2pkKFCjkmFydOnGD37t3A0wXGTE3lh0Jo69GjB/v27cPKyorjx48zrn835vWohKO1KbceJ/D24tMERybpO0whhBAvKd8Jhq+vL4cPH85125EjR146KPFqa9myJcePH8fR0ZErV64wuEd7vm9phaudGSFRyfT7/bQs+y6EEKVcoa2mqigKly9fply5coV1SPEKq1u3LmfOnOGNN97gwYMHvN29I/N7u1LTwZLH8akM+OMMASEx+g5TCCFEAem8FsmCBQv46quvSE1NRVGUbLc/kpOTMTY25vz581SvXr1Igi0sshZJyREVFUXPnj1p1aoVM2bMIDYpnRErz+IfHIO5sQEbxzeltpO1vsMUQgjx/xX6Ymc3btzg0qVL7Nixg+joaIYNG6a13crKCi8vL+zs7F4u8mIgCUbJkpqairGxMSqVCoCouETe++syZ+5E4lbOnJ0T35SVWIUQooQostVUQ0JCSEtLo2rVqi8dpL5IglFyJScn06FDBxo0aoKvXWfCE9IY0NCFH/t56js0IYQQ6P4bqtN/C5OSkoiLi8Pc3Bw7Ozvi4uJ49OhRjnXNzc2xtLQsWNTitbdnzx5OnTrFqVOn6NgzmPAaQ9hwLgQHa1M+aOeOgVql7xCFEELoQKdOngsXLsTR0ZHvvvtO8+/cHt99911RxyxeYX379mXlypUYGhqyf/smjPb/j8zkeOYdusnQpb48jk/Rd4hCCCF0oNMtksjISMLCwrC3t8fIyCjPGTrt7e0100KXVHKLpOTbt28fb7/9NvHx8ThVdqNM9y/JsKiAq50Za8Y0oaJt7pO8CSGEKDpF1gejKKSkpBATE0P58uVRq3UfORsTE4ORkRHm5ub5Op8kGKXD5cuX6datGyEhIdiWLYvLwO+ItayCk7UpK0Y2poaD3IoTQojiputvaIHmwdi5cycXL14EIDY2ln79+lGnTh0WLFiQr+NkZWUxadIkbGxsqFGjBk5OTmzatOmF+508eZIGDRpQsWJF3NzcGDFiBPHxMjHTq8bDwwNfX18aNmyIWqVi0ZjWuJUzJzQ2hW6/nmDGriASUjP0HaYQQogc5DvBePDgAVOmTNHMdTFr1ixu3brFkCFD+PLLLwkICND5WLNmzWLt2rWcO3eOmJgYpk2bxqBBg7hy5Uqu+1y8eJEOHTrQu3dvoqOjefToER06dODmzZv5bYooBRwdHTl69CiHDx+mWf06bPxPU9q/UZ6MLIUlJ+4yfPlZMrP0fhFOCCHEc/KdYBw5coSmTZtSpkwZALZt28bs2bOZOnUqo0aN4sCBAzofa+HChYwZM4Y6deqgUqmYMGECrq6uLFmyJNd9vvrqKxo0aMC0adMwMjJCpVIxZMgQGjRokN+miFLC3NwcT8+nw1TtLUzoVyGSN26uwdwgi/P3o1l64o6eIxRCCPG8fCcYKSkpZGQ8vSz96NEjbt++zZtvvgmApaUlKSm69fJ//PgxwcHBNGvWTKu8efPm+Pn55bhPeno6Bw8e5K233iIrK4vw8HAyMzPz2wRRisXFxTF48GD2/vMXyp4ZZCbH88uBG9x6LLfIhBCiJMl3gtGsWTO2bNnC4sWLmThxIh07dsTY2BiAs2fP0qRJE52OExERAZBt5k97e3vNtuc9efKElJQUIiIicHNzw8PDA3Nzc0aNGkViYmKu50pNTSUuLk7rIUonKysr1q9fj6WlJVf9fYn961MSH91n0l8BxKWk6zs8IYQQ/1++E4zatWszY8YM5s6dS0REBLNnzwYgICCApKQkOnbsqNuJ//9okWdXQ55JT0/HwMAgx32eTSW9fPly9u3bx+PHj7ly5Qr79u3jiy++yPVcM2fOxNraWvNwcXHRKUZRMnXq1IlTp05RuXJl4h8/4NGayZw7cYiRK/xIlE6fQghRIuhtmGp8fDxWVlb89ddfDBgwQFM+YMAAIiMjOXjwYLZ9MjMzMTc3Z/z48cydO1dT/vnnn7NlyxauX7+e47lSU1NJTU3VPI+Li8PFxUWGqZZyjx8/pl+/fpw4cQJQYdN6JB0HjmblyCaUMc45SRVCCPFyinSYKjxd/OzXX39lypQpzJkzRzNsVVeWlpY0aNCAffv2acrS09M5dOgQrVq10pRFRkYSGhoKgIGBAS1btiQ2NlbrWDExMXlOT25iYoKVlZXWQ5R+5cuX5+DBg4wZMwZQUMU/wudOFOP+PEdKuvTNEUIIvVIK4L///a9iaGioODs7K97e3oqrq6uiUqmUjz76KF/H2bp1q2JoaKgsWLBAOX/+vDJ48GClfPnySkREhKbO6NGjldq1a2uenzx5UjE3N1f++OMP5cqVK8ry5csVU1NTZcmSJTqfNzY2VgGU2NjYfMUrSqasrCxlw4YNyunrYcobX+9RKn+2UxmyxEcJCI7Wd2hCCPHK0fU3NN+3SK5cuULDhg1Zu3Ytffv21ZQfOXKE3r17s3XrVtq0aaPz8TZv3sy8efMIDw/Hw8ODGTNmUKNGDc32zz//HH9/f60rHUeOHOF///sf9+/fp1KlSowZM4b+/fvrfE6ZyfPVdeZ2JMOXnSF05zwsG3SncUMvPuxQndbVy2n68AghhCi4IpsqfNGiRZw8eZK1a9dm2zZ58mQsLCxK/IJnkmC82j764lvmzvwOlaEJdl0nYf5GSxpXKct/+9ShWnmZXlwIIV5GkfXBMDQ0zHWui+TkZAwNdVoBXogi8+1nH9G1a1eUjFSebP+JuOMr8b0dQbdfT7L85F2yZOZPIYQocvlOMNq3b8+uXbtYsGCBZmRGRkYGq1evZvny5XTu3LnQgxQiP6ytrdm+fTuffvopANFnNpG2awZJcdFM3xnEu2vPy3BWIYQoYgUapvrnn3/ywQcfkJCQgL29PZGRkRgZGfHf//6XSZMmFUWchUpukbw+Nm7cqJmIzd7BGbNun6Oyr0JNB0sWDK4vt0yEECKfiny59qioKE6cOMGDBw9wdHSkWbNmODg4FDjg4iQJxuslMDCQPn36EBoayvIt+/nZL5mI+FQM1CqGNqnExx1qYG1mpO8whRCiVCjyBKM0kwTj9RMTE0NAQACtW7cmLDaZr7de4UDQI1QqFVXszVk1sjGV7Mz0HaYQQpR4RZZgHDp0iC1btuS4Ta1WY21tTdOmTencubNmOvCSRhIMcerUKSZO/gyDth8QkWWOvYUxy4Y3oq6Ljb5DE0KIEk3X39B8D/lISkpi79693Llzh/r16+Pk5ERERAR+fn44Ojri5ubGzz//TNeuXdm8efNLNUKIopCVlcXYsWO5evUq5e/dpsqAr3iEK/1+P83HHWowrqUbBmqZM0MIIV5Gvq9ghIaG0qxZMzZv3oyXl5em/Pr163Tp0oWdO3diZmZGkyZNWLVqVYkcVSJXMMTNmzfp06cPV65cwcDAgIb93iOscgdUKhWWpobYmRvToLItX3erha25sb7DFUKIEqPI5sE4duwYbdq00UouAGrUqMGAAQPYv38/rq6uDBkyhAsXLuQ/ciGKgbu7Oz4+PgwaNIjMzEx8N/xKeZ9fKZOVTHxKBvcik/jnwkO6zz/JheBofYcrhBClTr4TjNjYWMLCwnLcFhoaSkxMDABGRkaYm5u/VHBCFCULCwvWrl3LokWLMDY25tzxAyRv/JR1Q2qwYkQjXO3MeBiTTP/fz7D0xB1ew/7QQghRYPlOMDp37szRo0eZNGkSgYGBREVFce3aNb766ivWrl1Lz549SU9PZ+/evXTt2rUoYhai0KhUKsaPH8+ZM2eoUqUKjRs3ommdqrSpWZ7tE9+kq4cDGVkKP+y6yn/+PC+rtAohhI4KNEz10KFDfPDBBwQFBWnKqlSpws8//0zfvn0JDQ0lODgYb2/vQg22sEgfDJGTmJgY1Gq15m8iLi4OgC2BUczYdZW0zCz61Hdmdv+6snCaEOK1VeTzYCiKwr179wgNDaVChQpUqVIFAwODAgdcnCTBEC+iKAr9+/fn0qVL/P3338SbOTFs+VkysxS+6FqTcS2r6jtEIYTQiyLr5Hn8+HEOHz4MPL1q0bx5c6pVq1ZqkgshdBEWFoaPjw83btygSZMmBB7azNfd3gDgf3uu8cv+6ySnye0SIYTITb4TjPv379OxY0eqVq3K999/T0hISFHEJYReOTk54e/vT+fOnUlJSWH8+PHsnPMpfWvbkKXA/MO3aD/7GOfuRek7VCGEKJHynWC88847hISEMG7cONasWYOrqyudOnViw4YNmtVVhXgV2Nvbs2vXLn7++WeMjIzYvHkzm78azPu1MnCyNuVhTDKDl/iy1f+hvkMVQogS56XXIjl58iQrVqxg48aNGBsbM2TIEMaNG0edOnUKK8ZCJ30wRH6dO3eOgQMHcvv2bVxdXfG/HMRnW66w70o4AAMbuTC5Yw3KWZroOVIhhChaRdYH43l169aladOmuLu7ExMTw6FDh/D09KR///4kJSW97OGFKBEaNmzIhQsXeOedd1i1ahU2FmVYNMSL8a2edvb8yy+E1rOOsMEvWM+RCiFEyZDvtUjgaQ/748ePs3z5cjZt2kS5cuUYOXIk27Ztw8XFhStXrjB48GCWLVvGxIkT8zzWzZs3Wbp0KeHh4Xh4eDB+/Pg8J+hatmwZBw4c0CpzcXFh1qxZBWmKEDqzsrJi9erVmudqtYrKMf584A5Hk1249CCWzzZfJiw2hUnt3GUoqxDitZbvKxiHDh3C3d2dDh06kJiYyObNm7lz5w7ffPMNLi4uANSuXZvBgwcTHJz3/+YCAgKoX78+oaGhNGrUiLVr19KyZUvS0tJy3cfPz4979+7Ru3dvzaNdu3b5bYYQL+3mzZuMGTOGyWMGU+naX/ynmTMAcw/e5KutgWRmycyfQojXV777YPz999/cu3ePESNGUK5cuVzrxcTEkJmZiZ2dXa51OnfujEqlYs+ePQA8efKESpUqMWfOHP7zn//kuM/48eN58uQJmzZtyk/YWqQPhigMqampfP7558ydOxeAmjVrMvizWSy/rkJRoHNtB+YOrIepkQzhFkK8OoqsD0bLli3p3LlzjslFeHg4Dx48AMDGxibP5CI1NZVDhw7x9ttva8rs7e1p27Ytu3btyjOGK1euMGrUKD766CP++eef/DZBiEJhYmLCnDlz2Lt3L46Ojly7do3pY/vQKukURiqFvVceMXiJD8GR0hdJCPH6yXeC8eeff+Lp6cmECRNIT0/Ptu3Z/+ZeJDg4mIyMDCpVqqRVXrlyZe7cuZPrfgYGBjRo0ABvb2+srKwYM2aMVpKSk9TUVOLi4rQeQhSWTp06cfnyZfr27UtGRgarfp2J5ZH/YWGk4kJwDJ3nHWeNz31ZLE0I8Vop0CiSDh06sG3bNjp06MCTJ08KdOKUlBTg6YqW/2ZhYaHZlpPp06ezdu1axo0bx3fffceBAwfYvHkzO3bsyHWfmTNnYm1trXk86ysiRGGxs7Nj06ZNrFixAgsLC7q2b82ej1rTpEpZktIy+WprIMOWnyU0JlnfoQohRLEoUILh6enJuXPnSE1NpVGjRly8eDHfx7C2tgYgKkp7JsTIyEhsbGxy3e/52y5eXl5UrlyZs2fP5rrP1KlTiY2N1Txk9lFRFFQqFSNGjODy5ctPOz2XNWP9WG/ebWCJYVoCJ24+odOc4/x1NliuZgghXnkFngfD0dGRo0eP0rp1a5o3b87mzZvztb+Liws2NjYEBgZqlV++fBkPD498HSs+Pj7PIYEmJiZYWVlpPYQoKq6urhgZGQGQkZHOxh8/In7dhzjEBhGfmsHn/1xmyFJfIhNk5lshxKvrpSbaMjExYcWKFfzwww8MHDiQlStX6ryvSqXSzJXxrE/EyZMnOXv2LEOHDtXU++OPP/j0008BSE9PZ8OGDVrHmTt3LpGRkfTo0eNlmiJEkQgLCyMxMZHH4Y/w/X0KFS+vxCgzmdO3Ixn4hw/hcbnfDhRCiNIs3wmGk5MTlStX1ir78MMP2bNnD7GxsZQpU0bnY82YMQNbW1tq1apFx44d6dy5M1OmTNGa1+Ls2bOaYawGBgbs3LkTd3d3evToQd26dfnuu+9YsmQJjRo1ym9ThChylStXxt/fn48//hiVSsWp3ZtIWPchpuGB3HycQP/FZ/C9E6nvMIUQotC99FokLysrKwsfHx/NTJ7VqlXT2u7n50dkZCSdO3fWlD148ICLFy9ia2uLh4cHlpaW+TqnzIMh9OHkyZOMHDmSW7duAeDQuBtGb45CbWRCs6p2jH6zCq1rlMdALTOACiFKLl1/QwucYOzYsYMdO3bw4MEDHB0dad++PQMHDiwV0yNLgiH0JTExkS+++IJff/0Vr0ZN6PjZ7/x9IZT0zKcfQ0drU6b3qkOHWhX0HKkQQuSsSBOMoUOHsmHDBtq1a4eLiwuPHj3i4MGDtGzZkj179qBWv/QaakVKEgyhb0ePHsXJyYnq1avzMCaZJYev8s/5EOIyDTFUq5g/qD5dPBz1HaYQQmSj629ovhc7O3XqFHv27OHixYvUqlVLUx4SEkLz5s3ZsmULb731VsGiFuI10bp1a82/nW3KkHB6HbGbNuE1+DMuZFVi4np/5oMkGUKIUivflxouXLhA7969tZILeDrsdOjQofj7+xdacEK8DhITE/nnn38Ivn+fLTMnYHZmMamJcUxc78+ey2H6Dk8IIQok3wmGtbU19+7dy3HbvXv35JaDEPlkbm7OpUuX+OCDD1CpVFw9voOoVe8RG3SC99ZdYOelUH2HKIQQ+ZbvPhgRERFUq1aNQYMGMWHCBJydnQkPD2flypXMmzePK1euZBsJUtJIHwxRUp05c4bRo0dz9epVAMq4e2PX8T0+6d2EiW2roZYRJkIIPSuy1VTLlSvH7t27OXXqFHXr1sXe3p7atWuzefNmtm3bVuKTCyFKsqZNm+Lv78/XX3+NoaEhWQ8DUbIymXPwBsNXnOVCcLS+QxRCCJ0UeJhqVlYWd+/e1QxTrVq1KgYGBoUdX5GQKxiiNLh8+TK3b98mzdmLr7YGkpaZRUZ8JC3rVeeLrm9Qx9la3yEKIV5DRTZMdfbs2YSGhvLzzz+/dJD6IgmGKG1uRyTw2dzVbP5xEtZN3sa62dsMbFKV6b1rY2JYOhJ7IcSrociGqZYrV46AgICXiU0IkU9Vy1lgH3kJMjOIPb2exKvHWBkygYTUDH4dVF9m/xRClDj57oPRrVs3fH19OXPmTFHEI4TIxaJFi9i4cSOOjo5kRIfyeMNXrPrvZCatOMq5e1GExSbrO0QhhNDI9y2SdevW8emnnxIaGoqbmxvlypXT2j5kyBAmTpxYqEEWNrlFIkqz2NhYvvrqK3777TcURUFtYk7ZjhMwr9WKzrUd+KxLTarYm+s7TCHEK6rIbpG4ubnx7rvv5rr9+Qm4hBCFy9ramvnz5zNs2DAGDhvFnWuBlLe1IlkFe6884uDVcMa2dGNSO3dMjaR/hhBCP/S+mqo+yBUM8arIyMhg586d9O7dmxvh8czcfZU9B49gVM4V94oV+LGfJ41cy+o7TCHEK6TIV1Pdu3cvp0+fpkGDBvTu3ZuHDx8SHR1NnTp1Chx0cZEEQ7yqnjx5QpWq7qRkKFi1HI5l3fYMb1qF99pUo7yVqb7DE0K8AorsFgnA2LFj2bx5M3Z2diQlJdG7d29MTEzo27cvFy5cwMLCosCBCyEK7tGjR1Sq6ERQUBBRe38l4eIe/gh7l1Vn7lPByoROtR2Y2NadcpYm+g5VCPGKy/coEj8/P3bv3s3169cZP368ptze3p6mTZuycePGQg1QCKG7OnXqEBAQwOzZs7G0tCQt7CaP/pxM5N5fCX30mNVn7tN61hF+O3KL1IxMfYcrhHiFFWg11R49elCuXDlUKu2x925ubty6davQghNC5J+RkREfffQRN27cYNiwYaAoJFzcT9yfE6llZ0BiWiaz9l2n67wTnLkdqe9whRCvqHwnGCYmJkRH57weQmBgYLZhq7qIjIzk2rVrpKSk5Gu/xMREAgICePjwYb7PKcSrzsHBgVWrVnHy5Enq1avHkEED2Dm5E3MG1MXewpjbEYkMWuLDDzuD5GqGEKLQ5TvB6NixI/v27ePkyZOaKxhZWVn89ttvbN26lR49euh8rPT0dIYPH46TkxMdOnSgfPnyLF++XOf9hw8fToMGDZg1a1Z+myHEa6N58+acO3eOWbNmoVar6FO/Ir91LU+ZE/NJj3nE0pN36Tn/FNsvhpKemaXvcIUQr4h8d/J0cnJi0aJFdOjQAWNjY4yNjfnjjz9ISkpiwYIF+VpNdcaMGezfv5/r16/j6urKunXreOedd6hfvz7169fPc99FixYRGRlZKkatCKFvBgYGmJv/3+Rb30/7kmun92HkdwTbxr242qgfH6yPx9HalGnda9HFw1GP0QohXgUFHqYaGhrK9u3befjwIXZ2dnTv3j3fS7U7OzszatQovv/+e01ZrVq1aNOmDb/99luu+12+fJlOnTrh6+tLjx49aN26NXPnztX5vDJMVbzuLl++zOTJkzlw4AAA5ta22LV8B2q2Q6U2oEsdB77s9gYVbc30HKkQoqQpsmGqFy9eJDk5GW9vb61RJM9ve5GwsDBCQ0Np3LixVrm3tzcXLlzIdb+kpCQGDhzInDlzcHFxyW/4QgjAw8ODffv2sWfPHiZPnsy1a9dI3PErFS7uRtVsBHuox8Gr4Qxo5MLHHWpQ1txY3yELIUqZfPfBOHDgAJs2bcp12+bNm3U6TmTk097rdnZ2WuV2dnY8efIk1/0++OADGjZsyIABA3SMGFJTU4mLi9N6CPG6U6lUdO3alUuXLjF//nzKli1LePAt3qkOzavZkZ6psMYnmK7zTnDuXpS+wxVClDL5TjDyEhISQtmyuk1LbGRkBDz98f+31NRUzbbn7dmzh02bNjF69GgCAgIICAggOTmZiIiIPJeQnzlzJtbW1pqHXPkQ4v8YGRnx/vvvc+vWLaZPn873X3zC2jHe/DXOm/LpYTwMe8SAP3xYfOw2WVmv3coCQogC0rkPxurVq5k9ezYRERGkp6fj5OSktT0xMZG7d+/i4+NDw4YNX3i8hIQErKysWLNmDYMHD9aU9+vXj/j4ePbt25dtnw0bNjBz5kytsuvXr2NhYYGzszPnz5/HwCD74k6pqalaiUxcXBwuLi7SB0OIPKSnp1PHw5O794Mx8+qNVaPedKjryqy368otEyFeY4W+Fsm5c+c4evQox44dIz4+nu7du2ttt7KyokWLFrzxxhs6B+nt7U21atVYs2YNACkpKTg6OjJ16lSmTJkCPL0qkpSURI0aNXI8Rr169aSTpxBFICQkhD59+nD+/HkADMyssWo2kAqNu9G/cRWGNXOlajlZFkCI102RLXZ25coVkpOTdbpK8SL79++nW7dufPvttzRt2pS5c+cSEBDA5cuXsba2BmDMmDH4+PgQGBiY4zEkwRCi6GRlZbFp0ya+/PJLzSy9hjYO2LR4B7M3WlC/UlmGNKnEWw0qolarXnA0IcSrQNff0Hz3wahduzYNGzZk7969TJs2ja1btwLw8OHDXJOA3HTs2JHdu3fj6+vL1KlTKV++PCdPntQkFwCVKlWiZs2auR6jRo0aODs757cZQggdqNVq+vfvT1BQEIsWLcLBwYGMmEc82TGL9AeBBITE8OmmS7y9+AzXHknnaSHE/ynQPBj/Xk21V69e/Pzzzzx58oRmzZqVitVU5QqGEAWTmJjI3Llz8fPz44/Vf7HpwgN+O3yLuJgoTCxsGN2iCpPauWNmXKCFmoUQpUCR3SLx8/Ojd+/eBAQEsHr1asLCwvj555+Bp1N3t2rVilGjRr1c9EVMEgwhXo6iKJqlAoLuPaR+7TcwdPHE+s0h2DhVwchAjYWJIaPerMJQ70qYGGbvfC2EKJ2K7BaJrKYqhPj3Z/+SzwnSkxNIun6KsOXvEbz5R548vMfDmGS+3xlEh9nH8b0jq7YK8bopEaupCiFKr4EDB3Lx4kX69OkDikJi0FHCl0/A9coqrNOjCI5KYvBSX5aeuEMBVyYQQpRCel1NVQjxavDw8OCff/7h/PnzdO/enczMTI7t/Jubi/5D52rmZGYp/LDrKq1mHeXHvdd4FJui75CFEEWsQJ08169fz6hRozSrqaampmpWU31+fZKSSPpgCFG0fH19+eabb6hUqRKLFy9m9Zn7/G/PNRIT4lCbmONgZcq6sU1w+9c8GrHJ6WRlKdjKJF5ClGhF1snzmcJYTVVfJMEQonikp6drpv4/e8Gfls2bU75xNxTPXlSo4MC4llWIT8nA904U5+5HkaVApbJmeLuVZah3ZTwr2ui3AUKIbIo8wSjNJMEQovh9+eWX/Pe//wVAbWSMuWcnrBq/haGVfa77NKxsy7c9a1PH2TrXOkKI4lWkCUZkZCSrV6/m7t27pKWlaW3r0KEDb731Vv4jLkaSYAhR/BRF4cCBA3zzzTf4+PgAoDYwxLNtbwaPncig9g2xNDUiICSGrf4P2XkplPRMBQO1ijEtqvBxh+oy3FWIEqDIEoyoqCjq1KmDiYkJnp6e2VY+7dmzJ8OGDStY1MVEEgwh9EdRFA4fPsz333/PsWPHAHB2dubevXsYGv7fBF2P41KYvjOInZfCAPCsaM1vgxvgUtZML3ELIZ7S9Tc039Pt7d69GxcXF06fPp3jyqVCCJEXlUpFu3btaNeuHSdOnGDGjBm0a9dOk1xkZWVx7do1atWqxYLBDehR9xGfbb7EpQexdPv1BMOautK7vjPVypfsGYOFeN3le5hqRkYGXl5eklwIIV5aixYt2Lt3L5988ommbMuWLdSuXZs+ffpw7tw5OtV2YNcHLajnYkNcSgYLjtyi/exj9FxwkuUn7xKTlJbHGYQQ+pLvBKNp06acPHmS9PT0oohHCPEa+vfMoP7+/qhUKrZu3UqjRo1o3749QX4n2fgfb34dVJ+2NctjoFZx6UEs03cG0eLHI8w7eBOfO5EcuhrO43iZY0OIkkCnPhgXLlzg+PHjmud//vkniqLQp08fLC0ttep6eXnRokWLwo+0EEkfDCFKtqtXrzJz5kzWrVtHZmYmAPXq1WPKlCkMHDiQqMQ0dl4KY/3ZYK49itfa18LEkFWjGuNV2VYfoQvxyivUTp4rV67ULGj2IiNHjmTy5Mm6R6oHkmAIUTrcv3+fOXPmsGTJEpKSkvD29ub06dP/mkVYYdflMJaevEtCSjop6Vk8jEnGwsSQpcMb0qRK2WxrJgkhXo7Mg5EHSTCEKF0iIyNZtGgRXl5edOnSBYDo6Gjmz5/PhAkTsLd/OpdGUloGo1ee48z/X1zN0dqUTrUdmNq1pgxxFaKQSIKRB0kwhCj9/vvf//Lll19SpkwZRo8ezUcffYSbmxvJaZl8suki+688Ij3z6dfb214V+amfp1zNEKIQFNly7UIIURJ4enrSoEEDkpOTWbBgAe7u7vTt25dzvqdZMKg+l77pxJwBdVGr4O/zD1h28q6+QxbitSJXMOQKhhCllqIoHDlyhJ9++ol9+/Zpyr29vTl58iQGBgYsP3mX6TuDAPBwtqaLhwPDm7pibpLvaYCEEJSiKxj+/v785z//oXfv3nz99ddER0fnWT8rK4uNGzcyevRoBgwYwPTp03n06FExRSuEKElUKhVt27Zl7969XLlyhXHjxmFqaoq7u7tmrp6RzV0Z1rA8ahVcfhjLT3uv0+3XE1wIzvu7RgjxcvR6BcPHx4fWrVszZswYmjZtyu+//05kZCTnzp3DzCzn6YDfeecdDAwMNDP/LVu2jKCgIC5cuICDg4NO55UrGEK8up48eUJKSgoVK1YEIDAwkEaNGtFvwCA8Og1i6z01obEpGKhVdPNwZKh3ZcqaG5OUloFbOQss5MqGEHkqFZ0827Zti7W1NVu2bAEgNjYWZ2dnZs6cycSJE3PcJyYmBhsbG83z5ORkrKysWLp0KcOHD9fpvJJgCPH6mD59Ot98843meZt27TGv141LajdUau2RJcaGalq6l2NwExfa1CgvnUKFyEGJv0WSnJzM8ePH6d27t6bM2tqa9u3bs3fv3lz3+3dyARAUFERGRgbVq1cvokiFEKXZ119/zcmTJ+nXrx9qtZojhw6y85ePSF/3HpVC9mGuysDK1BA7c2PSMrI4eDWcUSvPMWz5WW6Gx7/4BEKIHOntWmBISAiZmZm4uLholVesWJEjR47kue/58+f58ssviYuL49atW2zYsIGmTZvmWj81NZXU1FTN87i4uJcLXghRaqhUKpo3b07z5s25d+8eCxcuZNmyZYQ9fEDWoY0EB/+MsbExiqJwPTyev8894M8z9zlx8wnd5p/k2x61GdTYRa5mCJFPeruCkZb2dIGiMmXKaJWbmZlptuWmcuXKfPjhh0yYMAEPDw+mTZuWZ0fPmTNnYm1trXk8n9QIIV4Prq6u/PTTTzx48IAVK1bwww8/YGxsDDztQP75+GE4P/Zhx7uNaFW9HGkZWXyx5TIjVvjx55l73H2SqOcWCFF66K0PxoMHD3BxcWHnzp1069ZNUz569GguX77M2bNndTpOeno61atXp3///vz444851snpCoaLi4v0wRBCaOzcuZMePXoAYGtry8iRI7H16saKwBQys/7va7JXPSc+6VgDl7I5d0QX4lVX4vtgVKxYEXt7e/z9/bXK/f39qVu3rs7HMTIywtnZmbCwsFzrmJiYYGVlpfUQQoh/a9q0KT/++COurq5ER0cze/Zsvh7SjvInf6K10S0aV3r6vbEtIJR2vxxj6Yk7ZGW9dtMICaEzvc6DMWzYMJYtW8aTJ08A2LdvH/7+/lqjQebOncv48eMBSElJYdmyZfz7osvx48fx8/Ojbdu2xRu8EOKVYmdnx5QpU7h16xY7duygc+fOqFQqfE4cY9UPHzLZy5idE9+keTU70jKz+GHXVYYs9WX6jiD+u/sqf58L4WZ4PK/h3IVC5EivA76nT5/O5cuXcXd3x93dncuXL/Pf//6XN998U1MnMDAQHx8f4OnVCn9/fypWrEjVqlWJiYnh3r17fPbZZzoPURVCiLwYGBjQvXt3unfvzr1791i+fDl+fn40btwYlUrFmtFNGPrRNxy+l8LpjKaahdWeqVbegiFNKtHC3Z6KtmakZWbxKDaF8pYm2JgZ66lVQhS/EjFVeFBQEOHh4dSqVYsKFSpobbty5QoxMTE0b95cUxYbG0tgYCDm5ua4u7tjbm6er/PJPBhCiIJKSEjA0dGRhIQEzCytqd+uF3Xbv0WEYTkuPYghJT0rx/3sLYzZMqG59N0QpV6pmGhLXyTBEEIUVGxsLHPnzmXZsmWEhIRoyr29vRk0dBgm1Zuz+0Y8t8LjSUzLBMDIQEV6poJ7eQs2T2iGlamRvsIX4qVJgpEHSTCEEC8rMzOTffv2sWTJEnbs2EFm5tNk4rvvvmPatGkoikJMUjomRmrikjPo9dtJwuNSqetiQw9PRzycrantbC1Tk4tSRxKMPEiCIYQoTI8ePWLNmjWsXLmSnTt34urqCsD+/fs5efIkI0aMIMnEjrd/P0NyeqZmP5UK3Mtb0MPTib5eFXG2KZPLGYQoOSTByIMkGEKI4tCjRw927twJQMuWLen61iAM3Ly5EZVB4MNYQmNTNHXVKhjRrAqfdqpBGWOD3A4phN5JgpEHSTCEEMVh06ZNLF26lP3792uGr5qZmdG7d2+GDh1KoxZtOXo9gk3nQ/C5EwVAZTszPmjrTlcPR0k0RIkkCUYeJMEQQhSnBw8esHr1alauXMnNmzcBaNKkiWYIPsCR64/54p/LhP3/qxqWpoa0rVmeN6vZ41nRhsp2ZpgaScIh9E8SjDxIgiGE0AdFUTh79ixr1qyhYcOGmvl7oqOjadmyJT1698GoeksOhap5EJ2sta9KBS3dy/FD7zoy1FXolSQYeZAEQwhRkixZsoRx48Zpnjdp0oTmnXpjUasFlyPh5uME4lMyADA3NuDTTjUY2LiSXNEQeiEJRh4kwRBClCTx8fFs3bqVtWvXcuDAAbKynk7WZWBgQMeOHZk9ezaGZSsy9Z9L+N2LBqC8pQnvtq7KsKauGKhlKXlRfCTByIMkGEKIkurRo0ds2LCBtWvX4ufnh1qtJjQ0lAoVKpCZpTBn8zE2BMYTkfp0Kan6lWz45e26uJWz0HPk4nUhCUYeJMEQQpQGN27c4NSpU4wcOVJT1rp1a3x9fanj3Zpwu/pQqQGmZuaMfrMKE1pXxVJmCRVFTBKMPEiCIYQojVJSUqhfvz7Xrl3TlBkYm2BSpSFmNVvg7NGUUa3foFc9JwJCYrj2KJ5OtR2o52Kjv6DFK0cSjDxIgiGEKK0UReHSpUts3LiRDRs2cPv2bc22MtWaUP6tr7Xqq1QwslkVPurgLlc3RKGQBCMPkmAIIV4FiqLg7+/Pxo0b2bhxI92Hv8/dso25/DCWioYJhOxeSKKTF2WqNsbcypqudRypWt6CtIws3CtY0KFWBUwMZSSKyB9JMPIgCYYQ4lWjKAqZmZkYGBiQnJ7JHwsX8NFHHwGgUqsxcfHArHpTyrh7Y2hpD4CtmRGDm1Ti3dbVZNE1oTNJMPIgCYYQ4lV3/fp11q1bx5YtW7h8+bLWNsdqdSjXfTKxJuUAcLAyZVJ7dxpUssWlbBmyFMjMVEjPyiItI4vY5HQSUjMoZ2GCk00ZjA3V+miSKCEkwciDJBhCiNfJrVu32LJlC1u2bOHMmTOYmJgQ9igc3wdJ/Hf3NW74+6AyLoOxQ1VUqryTBwO1ijer2TOyuSst3cuhljk4XjuSYORBEgwhxOsqLCyMCxcu0K1bNwBS0jNxr+XBg1tXMbQoi6lbQ8pUa4Jp5bqojU0xMlBhXcYIM2NDHsenkJKepTmWZ0Vr/tvHgzrO1vpqjtCDUpVg3L17l/DwcGrUqIGtre0L66enp3Pt2jUMDQ2pWrUqxsbG+TqfJBhCCPFUWloagwcPZt++fSQkJGjKTU1Nadu2LQMHDuSdd94BnvbzuPskkTU+wWw8F0JCagZqFQxoVIkRzVyp4WCpr2aIYqTrb6heb6QlJyfTq1cvPDw8GDt2LE5OTsyZMyfX+oqi8N1331GxYkWGDBlCjx49cHV1Zdu2bcUYtRBCvDqMjY3ZtGkTT548Yd++fUycOBFXV1dSUlLYvXs3e/fu1aofee8qX3WryeFPWtGjrhNZCqw/G0ynuccZvMSHc/eiUBSFG+Hx7L4cxr4rj/C5E0lmlt7/LyuKmV67DX/77bf4+/tz+/ZtKlSowI4dO+jZsydNmzbF29s7W31FUVAUhevXr2NjYwPA999/z6BBg7h9+zaOjo7F3AIhhHg1mJiY0LFjRzp27Mi8efMICgpix44dNGrUSFMnMDCQxo0bU758ebp160bXrl3pNbABmwOj2R8UzunbkZy+fQZbMyOik9K1jl+prBkjm7tSy9EKqzJGPElIJSwmBQUFY0M1igIZmQrlrExoUMkW6zIyZ0dpp9dbJBUqVOC9995j2rRpmrK6devi7e3N4sWLdTpGeHg4Dg4O7Nq1i65du+q0j9wiEUKI/Pvnn38YOXIkcXFxmjIDAwO8vb1p1rod8RW92R+skJmlYGqkppbj0+/XO08SiXku4ciLSgX1XGx4x7sy3T2dZNRKCaPrb6jermA8fPiQx48f4+XlpVXu5eWFv7+/zsc5e/YsANWqVcu1TmpqKqmpqZrn//5wCCGE0E3fvn3p3r07J0+eZOfOnezdu5erV69y6tQpTp06xYEDB/hyQDPCYlOwUyVibmpMhQoVSErLYNP5B+y8FMaT+FRiktOxMzfGyaYMRgYqUjOedhw1VKu4F5nE3SeJ+AfH4B8cw8w91xjSpBJDmlSmnKWJnl8BkR96SzCio58uOVy2bFmtcnt7e822F3n8+DHvv/8+gwYNonr16rnWmzlzJt99913BgxVCCAE87bPRtm1b2rZty+zZs7l//z779u3j4MGDvPnmm5iamlLR1oyPP/6BOXPm0KBBAzp37kznzp1ZO8obI6MX3/oIj0vh73Mh/Olzn/C4VOYevMmvh27i4WxN06r2dK7jQN2K1qhUMkS2JNPbLZIbN25Qo0YNDh06RNu2bTXl77//PkeOHOHKlSt57h8dHU27du0oU6YM+/fvx9zcPNe6OV3BcHFxkVskQghRRAYMGMDGjRu1yiwtLWnTpg3t27dn/PjxL0w20jOz2BP4iBWn7uIfHKO1raJtGbp5OtLdw4k6zlaSbBSjEj9MNSkpCUtLS1auXKkZAgXQp08fUlNT2b17d677xsTE0KFDB4yMjNi3bx+WlvkbGiV9MIQQouiFh4ezf/9+9u7dy759+4iMjASgYsWKBAcHa5KCY8eOUa1aNZydnXM9VlhsMmduR3L42mMOX3tMUlqmZltlOzO6eTgyuEklKtqaFW2jRMlPMABatmxJhQoV+PvvvwFISEjA0dGR6dOna+bQv379OnFxcZqezLGxsXTo0AEDAwP27dtXoARBEgwhhChemZmZ+Pv7c+jQIQwNDZk8eTIAWVlZlCtXjqioKGrWrEn79u1p3749rVu3xto65wm8ktMyOXL9MbsuhXHoWrhm8i8zYwO+6laLquXM2X05DIA33cvRyPXpqBS5ylE4SkWCceLECdq1a8ekSZNo2rQpCxYsICQkBH9/fywsLAAYM2YMPj4+BAYGkp6eTsuWLblz5w5LlizRaliNGjV0HqYqCYYQQpQMjx49olevXpw7d46srP+bJVStVtOoUSPGjBnDmDFjct0/KS2Dw9ces/LUPc7dz73/nrmxATUdrfignTutqpcr1Da8bkpFggHg4+PDwoULCQ8Px8PDgylTplC+fHnN9lmzZhEUFMSKFSuIjY2lV69eOR7n008/1Ux9+yKSYAghRMkSHR3N0aNHOXjwIIcOHeL69esAfPnll/zwww/A09vjP/74I61ataJ58+Zat8ezshSWn7rLrH3XMTZQ07mOA6ZGBhy7EUFwVJLWuVpVL8f/3vLA0bpM8TXwFVJqEgx9kARDCCFKtpCQEA4dOoSXlxceHh4AmskY4en8G15eXrRq1YrWrVvz5ptvYmVlRVxKOiaGakwMDTTHSkrLIDQmhb/OBrPqzD3SMxXsLYyZP6gBTavaaZ03OjGNa4/iqeFgSVnz/C1D8bqQBCMPkmAIIUTp4+fnx8KFCzl27Bh3797V2qZWq9m4cSNvvfVWnse4E5HAe+v8uRr2dD4kSxNDrM2MMDZQk56VRUhU8tNyU0O+7VGbvg2cSUnPIiopjejENKzLGOFS9vXuSCoJRh4kwRBCiNItODiYY8eOcfToUY4dO8bt27e5fv26Zk6k33//ncWLF9O8eXPNo1KlSsDTTqJfbwtk0/kHOR7bxsxIM/OosaGatIwsre2ta5Sju6cTCSnpxKVkYGigwsLEkJbu5XC1z33KhFeFJBh5kARDCCFeLQ8fPsTJyUkzUuTtt99m06ZNWnUqVqyoSTZGjBhBlqEpkQlPZxbNzFJQFKhW3gIrU0MWH7/D3IM3SM98+hNpbKDGxsyIiIRU8vrV9KxozbCmrvSq54SRwas5xbkkGHmQBEMIIV5tYWFhnDx5UjONub+/P5mZT+fOMDAwIDY2VjNB4969e1Gr1Xh7e2v9JsQkpRGfkoGtuTHmxgaoVCruRyay4tQ9gkLjsLMwxsbMiIxMhdDYZHzuRGlWja1oW4Y3q9njZFOGZlXt8Kps+8oMk5UEIw+SYAghxOslMTGRs2fPcurUKR4/fsyvv/6q2daoUSPOnTuHWq3G09OT5s2b07RpU5o0aULVqlV1TgwiE1LZeO4By07e4UlCmta2Wo5WdK7jQGU7M4wN1ITHpWBjZkxXD8dSt5ibJBh5kARDCCEEgKIojBo1iqNHj3Lv3r1s2+vVq6e1AGdSUhJmZnl38kxOy2TflUfci0zkdkQi+6880izo9rxKZc2Y1M6dznUcMDfR2/Jg+SIJRh4kwRBCCPG80NBQTp06xcmTJ/H19cXf35+ePXtqZptWFIXy5ctTtmxZmjRponl4enpibJz7kNaYpDS2+D8kKDSO+1FJZGRmUcHKFL970TxJeLpOlomhmmZV7XCvYEmlsmZ4OFvzhqOV5upGWGwyQaFxOFibUstRv2uvSIKRB0kwhBBCvEhaWhoxMTGayR/v3btHlSpVstUzMTGhQYMGvPPOO7z77rs6Hz8pLYMVp+6xwS8k22RgAEYGKsoYPZ3PIy4lQ1PubFOG7nUdece7sl7WXpEEIw+SYAghhCiIyMhIzp49i6+vr+YRHf10ivLPPvuM//3vf5p6/fv3x8vLi4YNG+Ll5YWbm1uOVx4UReHao3h870QSHJXM7YgELj6I0QyVBTBQq3CzNyckOkmz9opaBW1rVqCflzMVbc04EBTOo9gUmlWzo03N8liZ5r1abUFJgpEHSTCEEEIUBkVRuHXrFr6+vtSpU4d69eoBsH//fjp16qRV18bGBi8vL7y8vDTJR17HDYtNITk9k8wshYq2ZTAzNiQ5LZOj1x+zxvc+p25F5rq/oVpFHWdrmlQpS/taFWjkWrZQ2guSYORJEgwhhBBFKTQ0lN27d3P+/HnOnz/PxYsXSUv7v5ElS5Ys0SzidvXqVVatWoWXlxf16tWjatWqqNUvHllyMzyeTRcesNX/ITFJ6bSsXg5XOzMOX3vM7YhETb3Rb1bh6+61Cq1tkmDkQRIMIYQQxSktLY0rV65w/vx5zp07x4cffkjNmjUBWLhwIe+9956mroWFBZ6entSrV4+6devSo0ePPFcLV5Snk4Sp1f93+yUkKgm/e1GcvRtFj7pONK9mX2htkQQjD5JgCCGEKCmOHTvGunXruHDhAoGBgaSkpGhtP3r0KK1atQLgxIkTnDp1inr16lGvXj0cHByKPV5df0NLx6BbIYQQ4hXVqlUrTQKRkZHBjRs3CAgI4OLFiwQEBODp6ampu3XrVmbPnq15Xr58eU2yUbduXbp3715i/uMsVzBKyBshhBBCvMiGDRvYunUrAQEB3Lhxg6ws7Qm8Hjx4gLOzMwCbN2/G0NCQXr16FWoMcgVDCCGEeMUMGDCAAQMGAE9nFQ0MDNRc6bhz5w5OTk6auqtXr6ZBgwaFnmDoShIMIYQQohQyMzOjcePGNG7cOMftXbp0oUGDBsUc1f+RBEMIIYR4BY0fP16v59f7Em7Hjx9nwIABtG7dmokTJxIWFvbCfXx8fBgxYgTe3t5cuHChGKIUQgghRH7oNcE4cuQI7dq1w93dnU8//ZSbN2/SvHlz4uPjc93nyy+/ZNKkSbzxxhv4+voSFxdXjBELIYQQQhd6HUXSvHlzXFxc+Ouvv4CnHVYcHR2ZNm0akydPznGf2NhYrK2tefDgAS4uLhw5coTWrVvn67wyikQIIYQoGF1/Q/V2BSMpKQkfHx969OihKTMzM6N9+/YcPHgw1/2sra2LIzwhhBBCvAS9JRghISFkZWVpDakBcHJy4v79+4V6rtTUVOLi4rQeQgghhCg6eksw0tOfLkNrYmKiVV6mTBnNtsIyc+ZMrK2tNQ8XF5dCPb4QQgghtOktwShb9unSsVFRUVrlkZGR2NnZFeq5pk6dSmxsrOYREhJSqMcXQgghhDa9JRhOTk44ODjg5+enVe7r60v9+vUL9VwmJiZYWVlpPYQQQghRdPQ60daoUaNYunQp48aNw9nZmc2bNxMUFMSKFSs0dWbOnElgYCBr164ttPM+GzgjfTGEEEKI/Hn22/nCQaiKHqWkpChvv/22Ympqqri7uytlypRRfvvtN606o0ePVmrXrq15vnPnTqVJkyZK/fr1FUB54403lCZNmihLlizR+bwhISEKIA95yEMe8pCHPAr4CAkJyfO3tkSsphoaGkp4eDjVqlXD0tJSa9udO3dISEjQLFcbERHB7du3sx2jYsWKVKxYUafzZWVlERoaiqWlJSqV6uUbwNOMzsXFhZCQkFfmFsyr1qZXrT0gbSotpE2lw6vWpqJqj6IoxMfH4+TkhFqde0+LErEWiZOTU7bhqs+4ublpPS9XrhzlypV7qfOp1Wqdk5H8ehX7eLxqbXrV2gPSptJC2lQ6vGptKor26DInld7XIhFCCCHEq0cSDCGEEEIUOkkwComJiQnffPNNtonDSrNXrU2vWntA2lRaSJtKh1etTfpuT4no5CmEEEKIV4tcwRBCCCFEoZMEQwghhBCFThIMIYQQQhS6EjEPRmkXHx/P9evXsbe3x9XVVd/h5FtSUhI3btygXLlyODs7Z9vu5+dHamqqVpmLiwuVK1curhB1lpWVxenTp7OVV69enfLly2uVpaenExgYiKmpKW+88UZxhZhvd+/e5eHDhzlue/PNNwFISEggICAg2/a6detmm7xOn4KCgoiJiaFZs2Y5blcUhaCgIDIyMqhduzaGhtm/onSpU1wURcHf3x+1Wk29evVyrBMTE8Pdu3epWLFitjl8FEXh1KlT2fZxd3enQoUKRRHyC2VkZHDu3DlsbW2pUaOG1rbExET8/f2z7ePp6ZltnoXExESuXbuGra1ttvmMiltiYiIBAQG4uLhQqVIlrW1BQUHZFt2Epyt7e3l5AfDkyROuXbuWrY63t7de/v4UReHWrVtkZmbi5uaGsbFxjvUePXpESEgIbm5uuS4iqkudlwlUvIQlS5YoZmZmSo0aNRRzc3OlS5cuSkJCgr7D0klYWJgyfPhwxdraWqlXr55ia2urNG/eXLl7965WPWdnZ8Xd3V1p3ry55rFo0SL9BP0C8fHxCqB4enpqxbt7926teocPH1YqVKiguLq6KnZ2doqnp6dy7949PUWdt3nz5mm1pXnz5oqNjY1Srlw5TR0/Pz8FULy9vbXqXblyRY+R/5+1a9cqjRo1UmxtbRUTE5Mc61y/fl2pWbOmUr58ecXFxUVxdnZWTp8+ne86xWX27NlK9erVFRsbG6Vu3brZtl+7dk3p0aOHYmtrq9SrV0+xsLBQevbsqcTExGjqJCcnK4Di4eGh9b7t2LGjGFvyVEJCgvL1118rLi4uiqWlpTJgwIBsdfz9/RVAadKkiVa8ly5d0qq3du1axdLSUqlevbpiaWmptGnTRqvdxSU0NFSZOHGi4ujoqJiYmCjffPNNtjpTpkzJ9vkyNTVVmjZtqqmzfv16xdDQMFs9fbRpwYIFiouLi+Lu7q64u7sr9vb2yurVq7XqZGRkKKNHj1ZMTU2VWrVqKSYmJsq0adPyXedlSYLxEi5duqSo1Wpl3bp1iqIoyuPHjxVXV1dl4sSJeo5MNz4+PsqqVauU9PR0RVGe/ji3bNlSadGihVY9Z2dnZcWKFXqIMP+eJRhnzpzJtU5sbKxiZ2enTJkyRVEURUlLS1PatGmTrd0lVUpKilK2bFnl008/1ZQ9SzCio6P1F1gevv76a8XHx0dZsmRJrglG/fr1lR49eigZGRmKoijKf/7zH8XJyUlJTk7OV53ikJGRoXz44YfKtWvXlMmTJ+eYYOzYsUPZvn275nl4eLji7u6ujBw5UlP2LME4ceJEcYSdp7t37yrffvut8vDhQ6Vbt255JhgRERG5HufmzZuKkZGR8scffyiKoijR0dFKjRo1tNpdXE6dOqXMmzdPiY6OVqpWrZpjgvG8kJAQRa1WK4sXL9aUrV+/XrGzsyvCSHX37D16ZtGiRYqBgYHWfybmzp2r2NjYKNevX1cURVFOnDihGBoaKtu2bctXnZclCcZL+Pjjj5Vq1applf3vf/9TrK2tNV+Apc3SpUsVQ0NDJTMzU1Pm7OyszJo1S/Hz81MePXqkx+he7FmC8ffffyvnzp3L8Qd39erVipGRkda2vXv3KoBy8+bN4gu2gNavX68Ami8GRfm/BMPPz0/x9/cvsVfRckswLly4oACKj4+PpiwkJERRqVTKli1bdK6jD7klGDn5/PPPFXd3d83zZwnGX3/9pZw7d06Jiooqoijz50UJho+Pj+Lv76/Ex8dnqzNt2jTF0dFRycrK0pQtWLBAMTU1VZKSkoo07rzommB8//33irm5uRIXF6cpW79+vVK2bFnlypUryuXLl5WUlJQijDR/MjIyFAMDA2XZsmWaMk9PT2X8+PFa9dq3b6/06tUrX3VelnTyfAn+/v6ae3TPNG7cmNjYWO7cuaOnqF6On58fVapUybaAzfTp0xkzZgxubm60bt2au3fv6ilC3bz77ruMGDGCChUqMGDAAGJiYjTb/P39cXNzw8bGRlPWuHFjzbaSbtmyZbRq1Yrq1atn29a3b18GDRqEra0tH374IRkZGXqIMP+eve4NGjTQlFWsWBFHR0fNNl3qlHTnzp2jWrVq2crff/99Ro4ciYODA2+//XaOfQJKkn79+jFo0CDKli3LxIkTSU9P12zz9/enQYMGWgtJNm7cmJSUlBz7MZQkiqKwYsUKBg4cmK3vUlRUFL169aJXr16ULVuWH3/8UU9RavP39yczM1Pzd5Wens6VK1dy/G169jnRpU5hkATjJURFRWXrFPPseUn/gsjJoUOHWLp0KV9//bVW+fTp04mKiiIgIID79++Tnp7OgAEDyMrK0lOkuTM0NGT16tVERERw+fJlrl69ip+fH++9956mTk7vm42NDWq1usS/b/fu3ePQoUOMHTtWq9zW1paDBw8SHBzM1atXOXnyJMuWLWPmzJl6ijR/oqKisLKywsjISKvczs5O857oUqckW7p0KUeOHOGLL77QlKnValasWMHjx4+5dOkS165dw9/fn3fffVePkebO2tqaffv2ERISwtWrVzlz5gyrV6/mhx9+0NQpzd+LR44c4c6dO9k+X1WqVCEgIICbN29y+/Zt1q1bxxdffMH69ev1FOlTiYmJjBkzhlatWtGiRQsAYmNjyczMzPE9ePb661KnMEiC8RKMjIxISUnRKktOTgbItVdvSeXn50ffvn355JNPeOedd7S2jRo1StNT2t7enhkzZuDn58ft27f1EWqeTE1NteJ3c3NjypQpbNq0SfO/+Zzet7S0NLKyskr8+7ZixQpsbGx46623tMqrVq1Ku3btNM8bN27M6NGj+euvv4o7xALJ6T2Bp5+nZ++JLnVKqq1btzJhwgQWLlyoGfkDT78nRowYofnffpUqVfj888/5559/SEtL01e4uapSpQodO3bUPPfy8mLs2LFaf2el+Xtx2bJleHh40KRJE63yJk2aULduXc3zXr160aVLF71+vlJSUujduzdpaWls3LhR8zf0LAHP6T3492fpRXUKgyQYL6Fy5crZhg8+e/78UKiS7Ny5c3Ts2JGxY8fyv//974X1nw2fy23oZElToUIF0tLSePLkCVB637esrCxWrlzJO++8g6mp6QvrV6hQodS8R5UrV9Z6jwAyMzMJDw/XvCe61CmJtm/fzsCBA5k3bx7jxo17Yf0KFSqQkZHB48ePiyG6l/f831lp/XzFxMTwzz//ZLt6kRt9fr5SU1Pp3bs3ISEhHD58WGsIvrW1NTY2Njm+B89ef13qFAZJMF5Chw4dOH78OLGxsZqybdu2Ub9+/cIfT1xEzp8/T4cOHRg1ahQ///xztu2JiYnZyvbv34+BgUGJnDsit3jLlSun+RB26NCB8PBwzp49q6mzbds2LCwsaNq0abHFml/79+8nODg4xy/AnNp94MAB6tSpUxyhvbSWLVtibGzM9u3bNWWHDx8mPj6eDh066FynpNm5cyf9+/dn9uzZOd72yO3vtWzZsjg6OhZHiPmiy99Zhw4d8PX11UqQtm3bhru7e4mcO+eZNWvWADB06NBs255vd2pqKsePH9fL5+tZcnHv3j2OHDmCg4NDtjrt27dnx44dmucZGRns2rVL63OiS52XJRNtvYThw4fz66+/0rNnTz7++GMuXLjA2rVr2blzp75D08m1a9fo2LEj9evXp0+fPpw8eVKzrUmTJhgZGXHy5El++uknhgwZgrOzM6dPn+ann35iypQpepsIKC/Lli3j9OnT9OzZExsbG3bv3s3SpUtZsmSJpuNqkyZN6NOnD0OGDOGHH34gKiqKr7/+mm+++QYzMzM9tyB3y5Ytw9vbO8cvtalTp5KRkUG7du0wNDRkzZo1nDlzhr179+oh0uyuX79OREQEt27dQlEUzd/as4nAypYty2effcbkyZPJysrCzMyMzz77jGHDhlGrVi0AneoUp4CAABISEnj48CGJiYmaNjVt2hQDAwOOHDlCv3796N+/P56enprtBgYGmkR21apVHD16VNNxcM+ePfz+++8sWrQIAwODYm/T6dOnycrKIjo6mpSUFE6ePImxsbGmE/RXX31FcnIy7du3x9jYmHXr1nH8+HF2796tOcaAAQOYPXs2vXr1YsqUKVy9epXFixezYcOGYm9Pamoqfn5+wNPbAcHBwZw8eRJra2s8PDy06i5btox+/fpha2ub7ThDhw6lZs2aNGvWjJSUFH777Teio6P58ssvi6Ud//b2229z6tQpli9fzu3btzW3qitXroyLiwsA06ZNw9vbm3fffZeuXbuyatUqUlJS+PjjjzXH0aXOy5LVVF/SkydP+PHHH/H398fOzo7x48fTpk0bfYelkz179jBjxowct23fvp2yZcsC4OPjw/Lly7l//z4uLi4MGjRI635/SbNjxw7+/vtvIiIicHNzY9y4cVr3T+HpF8/cuXM5fPgwJiYmDBgwgCFDhugp4hdLSUmhW7duvPvuu/Tr1y/b9szMTP7880/27NlDYmIiNWvWZOLEiSXmf4zffvstBw8ezFb+xx9/aJIDRVFYuXIlmzdvJiMjg06dOvH+++9rderUpU5xGTlyJDdv3sxWvnfvXiwsLFiyZAmrVq3Ktt3MzIz9+/drnu/atYuNGzcSHh6Om5sbY8eOpX79+kUae27atWuXbdZeOzs7tm3bBjy9TbdmzRp27dql+Tt7//33s81gHBMTw48//oifnx+2traMGTOGTp06FVczNB49epTj56V+/frMnz9f8zwkJIRBgwYxa9asHK9ipqSksHjxYo4ePQo8nbn0gw8+0MuV6tatW+c4OmzcuHEMGzZM8/zixYvMmTOH4OBgatSowZQpU6hSpYrWPrrUeRmSYAghhBCi0EkfDCGEEEIUOkkwhBBCCFHoJMEQQgghRKGTBEMIIYQQhU4SDCGEEEIUOkkwhBBCCFHoJMEQQgghRKGTBEMI8cq6cOGC1gy1QojiI1OFCyFeCefOnSM1NZXmzZtrypYvX86DBw+0VjAVQhQPSTCEEK+EpUuX8uTJE60Ew8vLq1CnPhZC6E4SDCFecykpKRw/fhxDQ0Pq169PSEgI0dHRtGrVSqve5cuXuX37Ni4uLtSrV09rMS5fX18URcHT0xN/f38SEhJo0qQJNjY22c6ny3Fq167NmTNnSE9Pp1u3bgQFBXHp0iUAzUJVFStWzHbM+Ph4/vrrLwDatGlD3bp1SUpK0jp/ZmYmvr6+hIeH4+7unm3xuPy0RQiRO0kwhHiNBQcH07p1awBq1KjBlStXqFixIlZWVpoEIz4+nv79+xMYGEj9+vW5ceMGlpaW7NixQ7NU9KJFi7h48SIpKSlUrlyZR48eERYWxrFjx6hZs2a+j5OYmIibmxs1a9akW7duXL9+na1btwIQFRXFqVOnmDZtGp999hkAV69e5e7du6Smpmrq1alTJ9stkvDwcDp37kxkZCS1a9fGx8eHDh06sH79ek2io0tbhBA6UIQQr63BgwcrzZs3V1JSUhRFUZTAwEDF2NhY6dSpk6bOuHHjlC5duiipqamKoihKZmam0rdvX2XIkCGaOsOHD1fKlCmjBAUFKYqiKFlZWUrbtm2VUaNG5fs4hoaGir+/f55xX7hwQTE2NlZu376tKfvPf/6jvPXWW1r13nvvPaVXr16a5yNGjFC8vLyUhIQERVEU5c6dO4qVlZWyePHifLVFCPFicgVDiNeUoij8888/rFy5EhMTEwBq165N165dSU5OBiAjI4M1a9YwYcIEtm/fjqIoKIpCxYoV2bRpk9bxWrduzRtvvAGASqWiZcuWHDhwIN/HadGiBfXq1csWb1xcHP7+/jx+/JjMzEysrKw4f/48bm5uOrd348aN/P7775ibmwNQpUoVBg8ezF9//cW4ceN0aosQQjeSYAjxmoqIiCAlJQVXV1etcldXV65evQrA48ePSUpKIiAggJCQEK16z26tPFO2bFmt5yYmJqSkpOT7OI6Ojtli3bx5M6NHj6Zq1apUqlQJExMTUlNTefz4sa7NJSIigqSkpGwJSdWqVdm/f7/ObRFC6EYSDCFeUzY2NqjVamJiYrTKo6OjNf+2tLREpVIxduxY+vfvX+Bz5ec4KpUqW9nEiRP5/vvvmThxoqbM3t4eRVF0jsHW1hYDAwOioqK0yqOiorC3t9f5OEII3chEW0K8poyNjWnYsCHbtm3TlKWmprJv3z7Nc0tLS5o1a8bixYuz/Zg/fPhQ53O9zHEyMzN58uQJNWrU0JQdPXqUyMhIrXoWFhZ5XmUwMjKicePG/PPPP1rH3rJli8yTIUQRkCsYQrzGZsyYQZcuXVCr1Xh4eLBmzRoyMjK0riIsXLiQdu3a0bZtW95++21SUlI4dOgQbm5uzJ8/X+dzFfQ4BgYGdO/enUmTJvHRRx8RERHBvHnzMDMz06rXsGFDli5dysKFCylbtixt2rTJdqxffvmFtm3bYmBgQOPGjdm4cSPx8fFMnTpV53YIIXQjCYYQr7H27dtz7Ngx/vzzT65cucJHH33Enj17iI2N1dTx9PTkypUrrFy5krNnz2Jvb88HH3xAp06dNHWaNGlCVlaW1rFr1aqlVaegxwFYs2YNixYt4vTp09jZ2bFv3z7WrFlD9erVNXX69+9PcnIyZ86cIS4ujjp16mSbaKtp06acP3+elStXcvLkSVq3bs369eu1bpHo0hYhxIuplPzcxBRCvFISEhIwMTHByMgIgLS0NN544w3GjRunmWNCCCEKQq5gCPEaCw8PZ8iQIfTv3x8TExPWrFmDoiiMGTNG36EJIUo5uYIhxGvu4sWL/PXXX0RGRvLGG28wduxYLCws9B2WEKKUkwRDCCGEEIVOhqkKIYQQotBJgiGEEEKIQicJhhBCCCEKnSQYQgghhCh0kmAIIYQQotBJgiGEEEKIQicJhhBCCCEKnSQYQgghhCh0kmAIIYQQotD9PwaeeIXaOGy8AAAAAElFTkSuQmCC",
"text/plain": [
"<Figure size 600x340 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"max relative deviation: 2.10%\n",
"(the publication-tolerance version of this check is a standing assertion in tests/test_scientific_validation.py)\n"
]
}
],
"source": [
"R = 200\n",
"H = np.mean([run_lineage(cfg(), seed=s).heterozygosity.to_numpy() for s in range(R)], axis=0)\n",
"t = np.arange(len(H))\n",
"law = H[0] * (1 - 1 / N_SAMPLES) ** t\n",
"\n",
"fig, ax = plt.subplots(figsize=(6, 3.4))\n",
"ax.plot(t, H, label=f'simulation (mean of {R} lineages)')\n",
"ax.plot(t, law, 'k--', label=r'$H_0(1-1/n)^t$')\n",
"ax.set(xlabel='generation', ylabel='heterozygosity $H$')\n",
"ax.legend(); plt.show()\n",
"\n",
"print(f'max relative deviation: {np.abs(H - law).max() / law.max():.2%}')\n",
"print('(the publication-tolerance version of this check is a standing assertion in',\n",
" 'tests/test_scientific_validation.py)')"
]
},
{
"cell_type": "markdown",
"id": "874cfaaf",
"metadata": {},
"source": [
"**This is model collapse.** Diversity falls geometrically, the rare tail goes first, and the\n",
"lineage converges on its own most common behaviour. Nothing about language models is needed to\n",
"produce it — only the fact that each generation is a finite sample of the last."
]
},
{
"cell_type": "markdown",
"id": "34f43d83",
"metadata": {},
"source": [
"## Closed form 2 — grounding is immigration, and it has an equilibrium\n",
"\n",
"Now mix `m` verified real samples into each generation. In population genetics that is\n",
"**immigration** from a non-drifting source, and drift no longer wins: the population settles at a\n",
"non-zero equilibrium diversity\n",
"\n",
"$$H_{eq} = H^* \\cdot \\frac{m(2n + m - 1)}{n + 2nm + m^2}$$\n",
"\n",
"(The textbook $\\theta/(1+\\theta)$ is the rare-immigrant limit of this; the paper uses the exact\n",
"form for the implemented model.)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "44195b4d",
"metadata": {
"execution": {
"iopub.execute_input": "2026-09-07T14:39:38.033172Z",
"iopub.status.busy": "2026-09-07T14:39:38.032931Z",
"iopub.status.idle": "2026-09-07T14:39:44.881068Z",
"shell.execute_reply": "2026-09-07T14:39:44.880729Z"
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAFFCAYAAADGjxacAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlcelbwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAaCBJREFUeJzt3XdcFMf/P/DX0UHhEJQqKKAiiA3BhooNxRo1do29R2OJxq4xiUnU+ElsMZZYEXsnVhILWEGKIgIqiIUq4NHb3fz+4Hf7daV4Bwd3wPv5eNwj3uzM7HvvCPtmd2ZWwBhjIIQQQghRIWrKDoAQQggh5FOUoBBCCCFE5VCCQgghhBCVQwkKIYQQQlQOJSiEEEIIUTmUoBBCCCFE5VCCQgghhBCVQwkKIYQQQlSOhrIDUAaJRIK4uDjo6+tDIBAoOxxCCCGk2mCMISMjAxYWFlBTq7zrHLUyQYmLi4OVlZWywyCEEEKqrTdv3qBhw4aV1n+tTFD09fUBFH24BgYGSo6GEEIIqT7S09NhZWXFnUsrS61MUKS3dQwMDChBIYQQQsqhsodI0CBZQgghhKgcSlAIIYQQonIoQSGEEEKIyqEEhRBCCCEqp1YOkpXKz88HY4wb6CMWiyEWi6GmpgYNDQ1ePQDQ1NRUaN2CggIwxqChocHNJZdIJCgsLIRAIICmpqZK1S0sLIREIoG6ujrU1dXlrssYQ0FBAQBAS0urUuqW9LnLU1eZ331N/jmpiu++oj8nVfXd0+8IxX+f9DuieN2CQjHuv0zG+8w8mBnWRXsbI6irCRTy3UuPo7LV6gRl8+bNWLVqFerUqQMAuHPnDm7cuIG2bdti8ODBXL3ffvsNBQUFmD9/PgwNDQEAAQEBuHr1Klq2bIlhw4Zxdbds2YLs7GzMnj0bJiYmAICQkBD4+PjA3t4eo0eP5uru2LEDIpEI06ZNg6WlJQAgLCwMZ8+eha2tLb766iuu7p49e5CcnIyJEyeicePGAICoqCgcP34cVlZWmDJlClf3wIEDiIuLw5gxY9CsWTMAQExMDLy8vGBqaopZs2ZxdY8cOYLY2FgMHz4cLVq0AAC8ffsW+/fvh5GREebNm8fVPXHiBJ4/f44vvvgCbdq0AQAkJSVh165d0NfXx6JFi7i6Z8+eRXh4OPr164f27dsDAFJTU7F9+3Zoa2tj2bJlXF0fHx+Ehoaid+/ecHNzAwBkZGTg999/h5qaGlavXs3VvXr1KgIDA+Hu7o7u3bsDAPLy8rBhwwYAwKpVq7hfNP/++y/u3buHTp06oU+fPgCK/uf65ZdfAABLly6Fjo4OAMDPzw+3bt2Ci4sLBgwYwO1vw4YNkEgkWLhwITfj6/79+/D19UXr1q0xZMgQru7//vc/5OXlYe7cuTA2NgYAPHr0CJcvX4ajoyNGjBjB1d22bRsyMjIwc+ZMmJmZAQCePHmC8+fPo2nTphg7dixX96+//kJqaiomT54Ma2trAMCzZ89w6tQpNGrUCJMmTeLq/v3330hMTMT48eNhZ2cHAHjx4gWOHj0KCwsLTJ8+nat7+PBhvHnzBqNGjULz5s0BAK9fv8bBgwfRoEEDzJkzh6t77NgxREdHY+jQoWjVqhUAID4+Hnv37oVQKMSCBQu4uqdOnUJkZCQGDhyIdu3aAQCSk5Oxc+dO6OnpYcmSJVzdCxcu4MmTJ+jbty86duwIABCJRNiyZQs0NTWxYsUKru6lS5cQHByMHj16oFu3bgCA7Oxs/PbbbwCAtWvXcnV9fX3x4MEDdOnSBb169QJQ9MtW+t0vX76cO6ndvHkT/v7+6NChAzw9Pbk+pHUXL15MvyPod0S1+h1xJSwef52/jbaFEYgX18WV/OYwF+pg7SBHvPI/X+HfEbt27UJVoFs8hBBCSA1xJSwes72C8CGLf5UjQZSL2V5ByMgtVFJk8hMwxpiyg6hq6enpEAqFSE5OhrGxsUpdwqPLt3T5tib+nNAtHvodoSo/J9XhdwRjDHl5eSgsLISamhoYY9ytlcLCQmhpaaFevXrc5/vq1auifiQM40/E4H1WIQRgUIcEDAKI//+1CAEACwNN+C5yh7aWZrm/+9TUVDRo0AAikahS1xKr1QlKZX+4hBD5iCUMD2NSkZSRCxN9He6+OVF90pOc9AQqPdlmZWXh/fv33LaPX2KxGHZ2dtztjsTERDx+/Jg7aX/6cnV1ha2tLYCiWw1Xrlwpta6Hhwd3izE6Ohp79uyBWCzm9vvx68svv0Tfvn0BFN3uWLlyZan9fvXVV5gwYQIA4OXLlxgzZgxv+8f9T58+nbtV9ebNG7i6ukIsFkMikRT774wZM7BlyxYARbfFTE1NS/2sJ02ahP3793Ofb926dQEA2lYtYTb2l89+V0end0QnO2P5vuCPVNU5tFaPQSGEqI4rYfFYdzEc8aJcrkx639zTyVyJkVUcYwxisRgFBQUoKChAfn4+jIyMuL9K4+PjkZSUxNv+8b89PDygp6cHAHjw4AGCgoK4Op+e9OfNm8ed3P755x9cuHChxOSgsLAQmzZt4sageHt7Y8eOHSXWKygowJEjR9CpUycAwN69e7Fo0SJenY//1vXx8eHGaZw8eRKTJ08u9bM5fvw4Ro4cCQC4ffs29++S7Nu3j0tQwsLCMHPmzFLr1q1bl0tQ3rx5g19//bXUunZ2dlyCkpaWhhMnTpRaV/oZAEVjWwICAkqt+/79e+7fjDEkJiaWWvfjgaefewCfRCLh1dXW1oaamhp0jM3KbCeVlJH7+UoqgBIUQojSSe+bf3o5V3rffOd45zKTlI9n4wFFg3LT09ORn5+PvLw83quwsBD9+vXj6l66dAnPnz8vVi8vLw/5+fnYuXMn1/eGDRtw48aNYgmE9L9BQUHcgNq5c+fiwIEDXL1PL1YnJCRwicT69euxY8eOUo/v5cuX3In57Nmz3IDPkowcOZLrNygoCLt37y617scDURMSEnD37t1S62ZlZXH/FovFyMjIKLVuYeH/jXPQ0tKCrq4uNDQ0uJf0loq6ujo3CBUA6tWrh1atWvG2f/z6+KqCmZkZvvjii1LrSgd1AoCVlRUWLFhQat2uXbtydRs1aoStW7eWWrdly5ZcXWtra1y8eLHUutKBzdJ4Hz9+DDU1Nairqxf778fPtTEyMkJKSgrU1NR49T7+t5Suri5yc4sSjnsvUzBmz/1SvxcpE32dz9ZRBXSLh27xEKJQeXl5yM7ORk5OTrGXuro6NwsDAI4ePYo3b9/hYKodspgmiu6S8wkAmAl10DTiMB4FBiA3N5dLHqSJRJ06dZCSksK16dWrF/77778S49PQ0ODGGADAkCFDcP78+TKPRzrGYdy4cfD29i61bmpqKjc2YObMmWUmB69fv+aeqr569Wrs3bsXmpqa0NTUhJaWFu/fx44dQ6NGjQAUzao5e/Yst11TU5N38l+yZAn3hFk/Pz/cvHmTt1360tTUxKBBg7iT/osXLxAWFlZiXQ0NDTg4OEAoFAIommn1/v37Uuvq6OhwYzpI1RFLGLps+A8JotxiyT7wf/8v+S/tWaFbp1V1DqUEhRIU8pHaNAYiOTkZmZmZyMrK4v4r/be+vj4GDRrE1V29ejXi4uJ4daT/btSoES5dusTVtbe3R1RUVIn7tLGxQXR0NPe+Xbt2eJpcINN9c9OnR/HQ50iJ23R1dZGdnc29/+KLL/Dvv/9CW1u7xNf9+/e5v0I3b96MgICAUusuW7aMS1Bu376NV69ecQnEp//t2LEjN8gwMTERmZmZxZIN6b/pBE4qg/RqJABekiL9Lfa5q5GyoASlElGCQkqi6mMgpJfV09PTkZ6eDpFIBF1dXTg7O3N1fvzxR6SmpnJ1Pn61bNkSx44d4+qamJggOTm5xH21bdsWQUFB3Hs7OzteYvGxJk2a4Pnz59z7Nm3aIDQ0FOrq6tDV1YWOjg50dXWhq6sLa2trXL9+nau7evVqPEyUINKo82ePf357A7SpV1hqIlHWoEJCapPK/l1GCUologSFfKq0MRCK/KtDijGGmJgYpKWlITU1FWlpadwrNTUVtra23EJZjDE0bdqUG1PxqR49evBuZdSvX593q+Njzs7OePToEffe1tYWCQkJqFu3LurUqYM6depw/7a3t8eff/7J1d22bRsyMzO5eh/XNzQ0RNu2bbm6WVlZ3JUCWch637yiMw8IqU0q82owzeIhpIqIJQzrLoaXeM+WoShJWXcxHB6OZqX+Dy4Wi3Hr1i0kJyeX+HJxceFWPJUmHR+PxP9Yz549uQRFIBDgw4cPvOREW1sbBgYGMDAwgLk5P2maPXs2CgsLue0fvxo0aMCr++LFi8/OFpD6eLXQz5EOEpVVexsjmAt1PnvfvL2NkVz9ElKbqasJqn1CTwkKqfUexqTyLoV+igGIF+ViytKfIY5/hri4OMTHx6NTp07Yt28fgKJEok+fPhCLxSX28fEMEzU1NVhaWkIikaBevXowMjJCvXr1uH87Ojry2vr6+kJPTw/16tWDgYEBtLW1S431xx9/lPm4ZU1OKpu6mgBrBzlitlcQBCj5vvnaQY41diwQIaRklKCQWic7OxuvX79GbGwsYmNjcSsmC0Czz7Y79c81ZD+7zb3/eMyDmpoaOnfuDDU1NTRo0KDYS/psFKnXr1/LHK/0mSY1maeTOXaOdy5239xMhcYAEUKqFo1BoTEoNU5+fj5iYmLw/PlzREdHw8jICOPHjwdQtD6Djo4O70qHrKsvdisMQiszXVhYWMDc3BzW1tbcw7aIYtSmWVSEVFc0BoWQMny8MBdjDAsXLkRkZCSioqIQGxvLS0Dc3Ny4BEVDQwMWFhb48OEDGjVqhEaNGsG6USPcEhQgi2mgrHU49i9dRSfLSlYT7psTQhSDEhSi0iQSCV6/fo2wsDCEhYXh6dOnCAsLg4WFBf755x8AReM7zp49y7ttoqenh6ZNm8LOzo43DRcAwsPDUadOHd64kM+tHUBjIAghpGpRgkJUxqfLlffv3x+3b9/mLbEtFRcXx3u/atUqCAQCNG3aFE2bNoW5uTmvr49JH6z1MRoDQQghqoUSFKIUjDE8f/4c9+7dQ2BgIAIDA5GRkYGwsDCuTn5+PrKysqCpqYnmzZvDyckJTk5OaNGiBZycnHgJzfTp0ysck6eTOTwczWgMBCGEqABKUIjcKjKQ8fDhwzhz5gzu3LlT4iqmqampMDIqWu9i8+bN0NTURNOmTWVe9KuiaAwEIYSoBkpQiFxkXUKZMYanT5/i33//xddffw0NjaIfNX9/f5w7dw5A0YJjrq6uaN++Pdq1awcXFxcYGhpyfbRu3bpKjokQQojqoQSFyKy05eATRLmY7RWEzcMckB/9EJcvX8b169cRHx8PAOjQoQM6duwIABgzZgyaNGmCLl26wNnZucxFxwghhNRelKAQmXxuOXgwhvkHbuPtX1MBVrSEu46ODtzd3fHxUjvdu3dH9+7dqyJkQggh1RglKEQmn1sOHgIB1A0aoElnTwzt3AJ9+/aFm5sbdHR0qi5IQgghNQYlKEQm8WmZMtX7bftufNHGspKjIYQQUtOpxtPCiMp6/fo1Vq5cia+nfiVTfRN9umJCCCGk4ugKCinVxIkT4eXlBYlEAgjUYN0tFYI69VDWcvDtbYyqPE5CCCE1D11BIRyJRMIb0NqwYUNIJBL06tULp0+dxI6pPSCAoFh6QsvBE0IIUTRKUAgYYzh37hxat24NX19frnzRokV49uwZfH19MWzYMAxs0xA7xzvDTMi/jWMm1MHO8c60HDwhhBCFEbCP/2SuJarqUdHVwZ07d/Ddd9/h7t27AICBAwfi4sWLZbapyEqyhBBCqreqOofSGJRaKjY2FgsWLOBWddXV1cXChQuxZMmSz7al5eAJIYRUNqUnKCkpKThx4gQSExPRsmVLDB06FGpqZd95ioyMxNWrV5GWlgZra2t8+eWXtf5KiDz27NmDBQsWIDs7G2pqapgyZQq+//57WFrS9GBCCCGqQaljUGJiYtCyZUscP34cWVlZ+PbbbzFo0KCiWSOl2L9/P1q2bIlHjx5BIpHg77//RpMmTfDixYsqjLx609HRQXZ2Nrp27YrHjx9jz549lJwQQghRKUodgzJy5Ei8ffsWfn5+UFdXR3R0NOzt7XHo0CGMGTOmxDYODg7o1asXtm/fDgAoLCxEkyZNMHbsWPz8888y7be2jUGRSCR4+/YtrK2tARQNivXx8cGAAQM+e7WKEEII+VhVnUOVdnYqLCzExYsXMX78eKirqwMAbG1t4e7ujjNnzpTazsjICAUFBdx7iUSCwsJCGBnR+hslSU5OxoABA9C5c2e8f/8eACAQCDBo0CBKTgghhKgspZ2hXr9+jdzcXDRp0oRX3qRJE0RFRZXa7u+//0ZERASGDh2KBQsWoGvXrhg0aBDmzp1bapu8vDykp6fzXrWBn58f2rZtiytXriAlJQWBgYHKDokQQgiRidISlKysLAAodnlIKBRy20ry9u1bxMbGQk9PDwYGBlBXV0dERAREIlGpbX755RcIhULuZWVlpZiDUGG7du1Cjx498O7dO9jb2+PBgwfw9PRUdliEEEKITJSWoNStWxcAiiUWHz584LZ9Kj8/H6NHj8aUKVNw5MgR/PDDD/D390d6ejqWLl1a6r6WL18OkUjEvd68eaO4A1ExEokEixcvxqxZsyAWizFmzBgEBgaiVatWyg6NEEIIkZnSEhRra2vUqVMHkZGRvPLIyEg4ODiU2CYxMREpKSno2LEjV6ampob27dsjLCys1H1pa2vDwMCA96qpNmzYgM2bNwMAfvjhBxw5cqTUhI8QQghRVUpLUNTV1TF06FAcOnQI+fn5AIDw8HD4+/tjxIgRXL0zZ87gjz/+AABYWlrCwMAAN27c4LYXFhbC39+/1KSmtpk1axbatWuHQ4cOYfXq1RAIaIVXQggh1Y9Spxm/e/cOXbt2Rb169dCuXTtcuHAB3bt3x9GjR7kT67Rp03D//n3uComXlxemT5+OXr16wdbWFjdu3IBIJMKtW7dgY2Mj035r2jRjxhgvERGLxdzMKEIIIUSRquocqvRn8WRmZuLChQvcSrK9e/fmbb98+TLi4uIwdepUruzt27e4efMmUlJS0KhRI3h6ekJHR+fTrktVkxIUiUSCKVOmwNXVFV9//bWywyGEEFLD1ZoERRlqUoKydOlSbNy4ERoaGoiIiICdnZ2yQyKEEFKD1fiF2kjF7d69Gxs3bgRQ9AgASk4IIYTUFJSgVFM3b97EnDlzAADr1q3D+PHjlRwRIYQQojiUoFRD8fHxGD16NMRiMcaNG4fVq1crOyRCCCFEoTSUHQCRjVjC8DAmFQkfsvHr98uRmJSMli1bYvfu3TSVmBBCSI1DCUo1cCUsHusuhiNelFtU0GIMrBr1wYLBTtDT01NucIQQQkgloFs8Ku5KWDxmewX9X3Ly/6nXNcaP/8XjSli8kiIjhBBCKg8lKCpMLGFYdzEcJc0Dl5atuxgOsaTWzRQnhBBSw1GCosIexqQWu3LyMQYgXpSLhzGpVRcUIYQQUgUoQVFhSRmlJyflqUcIIYRUF5SgqDATfdmW75e1HiGEEFJdUIKiwtrbGKGeNsCYpMTtAgDmQh20tzGq2sAIIYSQSkYJigpTEwDs0UkAAuCTRyZJVz5ZO8gR6mq0DgohhJCahRIUFebj44PQfw4i/fJmmOhr8baZCXWwc7wzPJ3MlRQdIYQQUnlooTYVduLECQDArAEd8dMKDzyMSUVSRi5M9Itu69CVE0IIITUVJSgq7ODBgxg8eDB69eoFdTUBOtkZKzskQgghpEpQgqLC1NTUMGLECGWHQQghhFQ5GoOiguLj45GTk6PsMAghhBCloQRFBc2bNw/W1tY4f/68skMhhBBClIJu8aiY6OhonD17FhKJBE2aNFF2OIQQQohS0BUUFfPXX39BIpGgT58+aNGihbLDIYQQQpSCEhQVUlhYiMOHDwMAZs+ereRoCCGEEOWhBEWFXLlyBQkJCWjQoAEGDBig7HAIIYQQpaEERYXs378fADB+/HhoamoqORpCCCFEeShBUREfPnzAxYsXAQCTJ09WcjSEEEKIctEsHhVhaGiI4OBgXL16FS1btlR2OIQQQohSUYKiQlq0aEEzdwghhBDQLR5CCCGEqCBKUFSAt7c3xo0bh//++0/ZoRBCCCEqgRIUFXDkyBF4e3vj3r17yg6FEEIIUQmUoChZeno6fH19AQBDhw5VcjSEEEKIaqAERckuX76M/Px8NGvWDA4ODsoOhxBCCFEJlKAomXTtkyFDhkAgECg5GkIIIUQ1UIKiRIwxbmCsp6enkqMhhBBCVAclKEoUERGB+Ph46OjooFOnTsoOhxBCCFEZtFCbEqWmpqJNmzZo0KABdHR0lB0OIYQQojIoQVEiNzc3BAcHo6CgQNmhEEIIISqFbvGoAHpyMSGEEMJHCYqSpKenIzs7W9lhEEIIISpJ5gQlJSWFTqgKtGvXLtSrVw8rVqxQdiiEEEKIypE5Qdm/fz8MDAzQsmVLTJo0Cdu2bcO9e/coaSmnmzdvIj8/Hw0aNFB2KIQQQojKkTlBGTduHPbt24cePXogKioKS5cuRefOnYslLVFRUXIFcOjQIbRr1w4NGzZEv3798Pjx48+2effuHWbMmIEmTZqgVatW2L59Oxhjcu1XmRhjePjwIQCgc+fOSo6GEEIIUT0yJyjm5uaYMGECtm7dirt37+KHH37AiBEj8Oeff6Jr166IiorCihUrsHv3bpl3fuLECUybNg3z5s3Dv//+i4YNG6JHjx5ITEwstU18fDw6dOiAlJQUnD17FmfPnsXLly/h7+8v836VLTY2Fu/fv4eGhgZat26t7HAIIYQQlVOhacbW1taYMWMG914ikSAzM1Pm9uvXr8fkyZMxadIkAMBff/2FCxcu4M8//8S6detKbLNy5UrUqVMHx48fh4ZGUfi///57tbqCEhAQAABo3bo1rX9CCCGElEChs3jU1NRgYGAgU12RSITHjx+jd+/eXJm6ujp69uwJPz+/EtswxnDmzBmMHTuWS06kqtNzbKQJiqurq5IjIYQQQlST0qYZx8XFAQBMTU155aampoiPjy+xTVJSEkQiEYyNjTFkyBBYWVmhQ4cO+PPPP8u8gpKXl4f09HTeS5keBgRC26oldO274t7LFIgl1efqDyGEEFIVZL7Fc+zYMVy7dg0uLi5wcXFBXl5ehXYsTSjU1dX5AWloQCKRlNhGLBYDAFasWIG//voLW7duxYMHDzBp0iRkZ2dj8eLFJbb75ZdfSr1lVNWuhMXjfef5MOuggTMJwJk992Eu1MHaQY7wdDJXdniEEEKISpA5QTE0NERkZCROnDiBrKwsAEVjUDIyMtC2bVu0bdsWrVq1gq6urkz9SafXJicn88qTk5NLnXprbGwMdXV1jBo1CmPHjuViuHXrFg4fPlxqgrJ8+XIsWrSIe5+eng4rKyuZ4lSkK2HxmO0VBPbJx54gysVsryDsHO9MSQohhBACOW7xeHp64s6dO0hPT8fTp09x6NAhDB06FOHh4ViyZAk6duwIfX19rF+/Xqb+GjRoABsbm2Kzb27fvo0OHTqU2EZbWxtt27YtlgTp6emV+TwbbW1tGBgY8F5VTSxhWHcxHCXdzJGWrbsYTrd7CCGEEJRjDIqamhocHR3x1Vdf4Y8//oCfnx9EIhGePXuGgwcPwsXFRea+vvnmG+zduxcPHjxAYWEhNm7ciLi4OMycOZOrs3DhQnTt2pV7v2jRIhw5cgTBwcEAgCdPnuDgwYMYNmyYvIdSpR7GpCJelFvqdgYgXpSLhzGpVRcUIYQQoqIU8jRjNTU1NG/eHM2bN5er3fz585GUlITevXsjLy8PlpaWOHPmDJo1a8bVycjIQFpaGvd+zJgxSEhIgIeHB7Kzs6Gjo4Np06Zh7dq1ijiUSpOUUXpyUp56hBBCSE0mYCqwgIhEIkFWVhb09fWLbcvMzERhYSEMDQ155YwxZGVloW7dunLvLz09HUKhECKRqMpu99x7mYIxe+5/tt7R6R3Ryc64CiIihBBC5FdV51CFXEGpKDU1tRKTEwClJiACgaBcyYmytLcxQoO6mkjKyINAUPzOmgCAmVAH7W2Mqj44QgghRMUobR2U2kZdTYDhNgyAAPjkopV0ibm1gxyhrlZ9FpwjhBBCKovcCUpCQgIKCwsrI5YaT+d9BJLP/QzNwixeuZlQh6YYE0IIIR+R+xaPl5cXNm3ahK+++gpTpkyBo6NjZcRVIz179gw5UfcwQqcnBk1ZiKSMXJjoF93WoSsnhBBCyP+RO0GZPn069PT0sG/fPmzevBkdO3bElClTMHr06FLHkZAi4eHhAACnFo40EJYQQggpQ4Vm8Tx58gR///03jhw5guzsbIwYMQJTpkxBt27dFBmjwiljFg8AnD9/Ho8ePcLEiRNhZ2dXZfslhBBCFKWqzqEKmWb87t07jB07Frdv3wYAODg44Pvvv8fIkSMrHGBlUFaCQgghhFR3VXUOLfcsnsLCQly8eBFDhw6FjY0N0tLSsGXLFoSHh2P06NGYMmUKLly4oMhYCSGEEFJLyD0G5c2bN9i+fTsOHTqEzMxMjBo1Cn5+frzn56xZswaFhYW4e/cuBg8erNCAq6vHjx8jPj4erVu3hpmZmbLDIYQQQlSa3FdQTpw4gZs3b+LHH39EfHw89u7dW+LD/QYOHIi+ffsqJMiaYO/evfD09MT//vc/ZYdCCCGEqDy5r6B06NABzs7O6NGjR5n12rdvX+6gaiLpDB4HBwclR0IIIYSoPrmvoAQHB+P8+fOVEUuNFhkZCQByP1CREEIIqY3kTlBat26N+/fvQwWeMVht5Ofn4927dwAAW1tbJUdDCCGEqD65b/HY2NhALBZj2LBhGDt2LBo0aMDbbmVlRWt8fOLt27dgjEFHRwcmJibKDocQQghReXInKMePH0dgYCACAwNx7ty5Ytu//fZb/Pbbb4qIrcaIjY0FAFhbW0MgoCXtCSGEkM+RO0H59ttvsWDBglK3q6nRA5I/JU1QGjdurNxACCGEkGpC7gRFIBBAQ0PuZrWau7s7Dhw4AGNjev4OIYQQIotyZxo3btxASEgIUlNTeQNmu3TpAk9PT4UEV1PY2NjAxsZG2WEQQggh1Ua5EpTRo0fj6tWrEAqFEIvFqFOnDiIjI9GgQQNYWloqOkZCCCGE1DJyDxjx9/eHv78/Xrx4gblz52LUqFGIiIiAv78/JBIJ+vTpUxlxVmve3t64evUqsrOzlR0KIYQQUi3InaCEhYVhwIABMDY2hoaGBvLy8gAAbm5uGDduXIkze2ozsViMiRMnwtPTEykpKcoOhxBCCKkW5E5QMjMzoa+vDwAwNTVFTEwMt01PTw8ikUhx0dUA8fHxKCwshIaGBiwsLJQdDiGEEFItVGg6jru7O6ZPn46NGzeiXr162LVrFw4ePKio2GqEV69eAQAaNmwIdXV15QZDCCGEVBNyJyh9+/ZFTk4OAMDCwgJ///03fvjhB6Snp2POnDkYNGiQwoOszmgNFEIIIUR+cicoLVu25L0fNWoURo0apbCAahppgtKoUSMlR0IIIYRUHxW6xZOamor8/HxeWZ06dbgxKoSuoBBCCCHlUa516X/++WcIhUIYGxvD3Nyc91q3bp2iY6zWpGNQ6AoKIYQQIju5r6DcuXMH69evx9atW+Hs7AxNTU3e9vr16yssuJrgp59+wtixY9G1a1dlh0IIIYRUG3InKI8fP8bo0aMxderUyoinxnF1dYWrq6uywyCEEEKqFblv8VhbWyMzM7MyYiGEEEIIAVCOBKV379548eIFTp8+DbFYXBkx1RipqanYt28frl69quxQCCGEkGpF7gRlz549CA8Px/Dhw6GrqwtDQ0Pea/Xq1ZURZ7UUERGBqVOnYvbs2coOhRBCCKlW5B6D0qdPH5iZmZW63d7evkIB1SRJSUkAih4JQAghhBDZyZ2gNGvWDM2aNauMWGocaYJiYmKi5EgIIYSQ6qXcC7X9999/8PHxwdu3b2Fubo7evXvTMvefSExMBEAJCiGEECKvci3UNmvWLPTu3RsPHz6Euro6Hj9+jGHDhmHYsGFgjCk6xmqLbvEQQggh5SP3FZQHDx7g2LFjCAgIQLt27bjyqKgodO/eHRcvXsTgwYMVGmR1RVdQCCGEkPKR+wpKQEAAhg4dyktOgKKxKRMmTMDDhw8VFlx1R1dQCCGEkPKR+wpK3bp1ER8fX+K2+Ph4ODg4VDiomuLnn39GdHQ0OnXqpOxQCCGEkGpFwOQcNJKQkIAmTZpg1qxZmDNnDiwtLZGYmIgDBw7gp59+QmhoqMonKenp6RAKhRCJRDAwMFB2OIQQQki1UVXnULmvoJiZmeH8+fOYOXMmNm/ezJVbWFjg2LFjcicnjDEEBgYiMTERTk5OaNy4scxtExIS4O/vj6ZNm6J169Zy7ZcQQgghqkvuKyhSBQUFCA8P56YZOzo6QkdHR64+0tPT0b9/f7x8+RLNmzfHw4cPsXDhQvz000+fbSsWi9G9e3c8ePAAc+bMwR9//CHXfis7+0tPT8fJkydhYWGBfv36Vco+CCGEkKqmsldQpDQ1NdG6desKXblYtWoVEhMT8ezZMxgaGuLWrVvo3r07evbsiZ49e5bZ9vvvv0fDhg2RkZFR7v1XppcvX2LatGkwMzMrdcwOIYQQQkomU4Jy5coVHDt2DP369YNQKMSxY8dKrduvXz+MGjXqs30yxuDl5YXvvvsOhoaGAAB3d3e4urrCy8urzATlxo0b8PLyQnBwMLp37y7LIVQ5msFDCCGElJ9MCYpEIkFhYSEkEgn377LqyuLt27dIS0tDq1ateOWtWrVCSEhIqe3ev3+PCRMm4PDhw1xi8zl5eXnIy8vj3qenp8vUriJomXtCCCGk/GRKUPr374/+/fvz3leUSCQCANSrV49XbmxsjA8fPpTYhjGGiRMnYvz48XJdOfnll1+wbt268oZaLrRIGyGEEFJ+ci/Udvz4cezcuVPubZ/S1tYGAGRnZ/PKMzMzSx1se+LECdy5cwetWrXCqVOncOrUKYhEIrx48QKnTp0qdZn95cuXQyQSca83b97IFGNF0C0eQgghpPzkHiT75s0bJCQklLgtNjYWqampMvVjZWUFDQ0NvH79mlf++vVr2NralthGKBSid+/eOH36NFeWlpaGyMhIHDt2DMOGDYNAICjWTltbm0uIqgpdQSGEEELKT+YEJTU1FUlJSUhOTkZqaioiIiJ427OysnD16lUMGTJEpv50dHTQo0cPnDp1CpMnTwYApKSk4N9//+Wtr/Lo0SOkpqbCw8MDnp6e8PT05PXTpk0bdO/eXa5pxlWBrqAQQggh5SdzgrJv3z4sWbKEe79///5iddq0aYOxY8fKvPNff/0VXbt2xeTJk9GpUyfs2bMHzZo1w5QpU7g6O3fuxP379xEWFiZzv6rg+++/x/jx42mZe0IIIaQcZE5QZs2ahdGjR2PXrl1ITk7GqlWreNuFQiH09fXl2rmzszMePXqE3bt349atW/jyyy/x9ddf827HuLi4lLkQjIeHB1q0aCHXfqtChw4d0KFDB2WHQQghhFRLcq8kKxaLwRiDhka513hTOnoWDyGEEFI+KruSrLq6OgAgMjISjx8/RmpqKm/2TJs2bdCxY0fFRVgN5eTk4NixYxAKhRg6dGiJA3cJIYQQUrpyXQZZvHgxtmzZgrp164IxBrFYjMzMTBgYGGD58uW1PkFJTEzElClToKOjg5ycHGWHQwghhFQ7cq+D8ujRIxw4cACRkZFYuXIlpk2bhoyMDBw/fhxaWlqYMGFCZcRZrUgXoRMKhUqOhBBCCKme5E5QgoKCMHjwYNja2kJNTQ35+fkAgJEjR2LYsGE4deqUwoOsbihBIYQQQipG7gRFJBLByMgIQNEiZG/fvuW2mZqacut/1GaUoBBCCCEVI3eC8rHOnTvj+vXr8Pb2xuXLl7Fv3z44OTkpKrZqixIUQgghpGLkHiTbuXNnZGZmAgBsbW3x008/Yc6cOcjOzsbYsWMxcuRIhQdZ3UgTFJrCTAghhJRPuRKUj82fPx/z589HYWFhtV4bRZHoCgohhBBSMXJnFF5eXkhJScH8+fP5HWlolLqtthk2bBhsbGzQqFEjZYdCCCGEVEtyJygJCQmlPs343bt3SEtLq3BQ1V3z5s3RvHlzZYdBCCGEVFsyJyiJiYl48+YN3r59i5SUFAQGBvK2Z2VlwcfHB2PGjFF4kIQQQgipXWROUA4fPsx7mrGXlxdvu5qaGjp37izX04xrKl9fX2RkZKBDhw6wsLBQdjiEEEJItSPzwwLz8vKQk5ODP//8E4mJiVi3bh1ve926davNINnKftCRm5sb7t69i9OnT2PYsGEK758QQghRFpV7WKC2tja0tbWxYsWKSgumpqBZPIQQQkjFlGuhtt27dyMgIAAAEB8fDxcXF+jr6/NuAdVmlKAQQgghFSN3gvL8+XNs27YNzs7OAIBff/0VQqEQu3fvxuHDh+Hn56fwIKsbSlAIIYSQipF70Iifnx86d+4MdXV1AMClS5dw+PBhdOzYEaGhobhz5w66du2q8ECrC7FYjIyMDACUoBBCCCHlJfcVFIFAgA8fPgAAoqOjkZiYCFdXV26bjGNuayxpcgJQgkIIIYSUl9wJSvfu3XHx4kUsX74ckydPxpAhQ7irKX5+fnB3d1d4kNWJ9PaOdFAxIYQQQuQn9y0eGxsbeHt7Y9euXbCxscGmTZsAACEhIbCxsSn2rJ7axsjICEeOHEFeXp6yQyGEEEKqLZnXQalJqmoONyGEEFLTVNU5tFzTjAkhhBBCKlP1WPq1GomOjkZoaChsbGzQpk0bZYdDCCGEVEt0BUXBrly5gmHDhuHHH39UdiiEEEJItUUJioLRIm2EEEJIxVGComCUoBBCCCEVRwmKglGCQgghhFQcJSgKRgkKIYQQUnGUoCgYJSiEEEJIxVGComCUoBBCCCEVR+ugKNiKFSsQGxsLZ2dnZYdCCCGEVFuUoChY//79lR0CIYTIjTGGwsJCiMViZYdClExdXR0aGhoQCARKjYMSFEIIqeXy8/MRHx+P7OxsZYdCVISenh7Mzc2hpaWltBgoQVEgiUSCCxcuwMDAAN26dYOGBn28hBDVJpFIEBMTA3V1dVhYWEBLS0vpfzkT5WGMIT8/H8nJyYiJiUHTpk2hpqac4ap0BlWgjIwMDB06FACQk5NDCQohROXl5+dDIpHAysoKenp6yg6HqABdXV1oamoiNjYW+fn50NHRUUocNItHgaQzeLS0tJT2hRJCSHko669koppU4edB+RHUIOnp6QAAAwMDJUdCCCGEVG+UoCiQdIAZXSYlhBDl8vHxgb29vbLDwI0bN9C4ceMy6/j6+qJz586wtLTEmjVrqiawaoAGSShQTk4OgKL7d4QQQpQnNzcXiYmJyg4DeXl5SEhIKHU7YwxjxozBkiVLMH78eOjr61dhdKqNEhQFogSFEEKIPJKSkvD+/XsMGjQIFhYWyg5HpSj9Fk9gYCCmTp2KgQMHYvny5UhJSSmzfkZGBrZs2YJRo0Zh3Lhx2LlzJ/Lz86so2rJRgkIIIVUnIiICY8eOhZ2dHdq3b4/9+/eXWf/OnTvo06cPGjVqBFdXV+zatYu3XSQSYe7cuWjRogWcnJywcOFCbmyhdH+jRo2Cra0tnJ2dsXr1auTm5vL6uHr1Kjp16oRmzZphxIgRiImJKTWeK1euoGXLlgCAbt26wczMDEFBQTLFevz4cTg7O+Po0aNwdXWFpaUl4uPjMX/+fEybNg1Lly6Fi4sLmjZtihUrViAzMxNLlixB8+bN0apVK+zYsePzH7CSKTVBuXv3Ltzc3KCvr4+vvvqKe5+VlVVifYlEgpYtW+LVq1f48ssv4eHhgc2bN6N///4qsfph69atsX37dixYsEDZoRBCSIVkZWWV+vr0pFxWXekfbp+rK68XL16gY8eO0NbWxqlTp/Dnn3/Cz88Pjx8/LrF+VFQUevbsCVdXV1y7dg0LFy7E4sWLsXv3bq7Ot99+i5CQEBw+fBinT5+GnZ0d/vjjDwBATEwM3Nzc0KpVK1y5cgW7du3CjRs3MG3aNK59eHg4Bg0ahAEDBsDHxwd9+vTBokWLSj2GHj164OrVqwCAs2fPIiQkBC1btpQp1pycHISEhGDnzp34888/8ejRI5iamkIkEmHfvn0AAG9vb2zatAm//fYbWrRoAR0dHVy4cAErV67EN998g3v37sn9uVcppkTu7u5s6NCh3Pv09HRWt25d9scff5RYXyKRsNTUVF5ZYGAgA8Du378v835FIhEDwEQiUfkCJ4SQGiInJ4eFh4eznJwcXjmAUl/9+/fn1dXT0yu1rru7O69u/fr1S6wnrylTprB27doxiUTCK5e+P3nyJBMKhbz6Xbt25dX94YcfmIWFBfe+Q4cObOPGjSX2N2PGDDZy5EjetvDwcAaAvX//njHG2OTJk1nv3r15dWbPns20tbVLPY43b94wAOzZs2dyxbp//34GgL1+/ZpXb+LEiaxVq1a8sl69erEOHTrwyjp06MB++umnUuMq7eeCsao7hyrtCkpOTg78/PwwZMgQrkxfXx+9evXCtWvXSmwjEAhQr149XpmxsTEAlCsDJ4QQUj09ePAAvXv3LrbqbWmr4D5+/Bhdu3bllbm7uyMuLg6pqakAgNGjR+PHH3/E/Pnz4ePjg6ysLK6/e/fu4cqVK2jYsCEsLS1hYWGBHj16ACi6miPdh5ubG28fXbp0kfvYZIkVKDr/WVlZFWvv6OjIe1+/fn04ODgUK3v//r3csVUlpQ2SffPmDSQSCRo2bMgrb9iwIW7cuCFzP5s2bUL9+vXRoUOHUuvk5eUhLy+Pe//xPUVFiomJwatXr2BlZYUmTZpUyj4IIaQqZGZmlrpNXV2d9z4pKanUup8u+PXq1asKxSVVWFgIbW1tmevn5+cXe66M9L10HOOCBQvg5uaGc+fO4aeffsKYMWOwc+dOjB8/Hnl5eZg8eTKWLVtWrG/pH8oFBQXQ1NQscR/ykCVWAKUuCPrp91NaGWNM7tiqktISFOmH/OmAUj09PZkHvf7999/YtWsXzp8/jzp16pRa75dffsG6devKH6yMDh8+jLVr12LGjBnFBjQRQkh1Utbv1KqqW5YWLVrgwYMHMte3t7dHcHAwrywoKAgGBgYwMzPjylxdXeHq6or169fj559/xnfffYfx48fDyckJgYGBvLqfatasGUJDQ3llISEhMscob6w1ndJu8Uhv1Xx8uQoAUlJSit3GKYmXlxdmz56NQ4cOYcCAAWXWXb58OUQiEfd68+ZN+QMvg3QwGC3URgghlWvhwoXw9fXFtm3bUFBQgJycHGzbtg0REREl1v/mm2/g4+OD06dPgzGGZ8+e4ddff8U333zD1fn666/x6NEjSCQSFBQUIDY2lksIFi9ejAcPHuD7779Hbm4uGGMICwvD9OnTufazZ8/GuXPn4OPjA6BolurOnTvlPjZZYq0NlJagWFpaokGDBtyUKqmgoCC0adOmzLbe3t6YOnUq9u3bh7Fjx352X9ra2jAwMOC9KgNNMyaEkKrRpUsXnD59Gtu2bUOdOnVgaWmJly9fwsbGptT6O3fuxJw5c1C3bl04OztjwIABWLVqFVdn4MCBmDVrFoRCIYRCIZ48eYKDBw8CADp16gQfHx+cO3cOdevWhYGBAcaPHw9PT0+ufc+ePbF+/XqMGTMGderUwejRozFlypRyHdvnYq0VKnUI7mcsXryYNWrUiCUmJjLGGPvnn38YAObv78/V2bx5M5s2bRr3/ujRo0xLS4t5eXmVe7+VNQJ5xowZDAD74YcfFNovIYRUlrJma1QXHz58YIWFhbyy3Nxc7tzyMYlEwlJSUlh+fn6p/WVnZ7O8vLxSt2dmZpb5eeXn57OMjAzGGGN5eXksISGh1LpisZjFx8cXi/9zsWZnZ7OkpKRi5R8+fGAfPnwoVvbp+S4tLY2lp6eXGpcqzOJR6kqy69atw9OnT9G0aVPY2dnh2bNn2LhxI28UdHh4OO7fvw+gaBGdr776CkKhEPv37+ctyvPtt9+iX79+VX4MH6MrKIQQUvWEQmGxMm1tbZiYmBQrFwgEMDIyKrO/z/0O/9w4Gk1NTW6wrJaWFkxNTUutq6amVuq4krJi1dXVLTHOkj6LksoMDQ1LjUlVKDVB0dPTw6VLlxAVFYXExEQ4ODigfv36vDrffvstRCIRV//y5csl9vXptCploASFEEIIUQyVeBZPs2bN0KxZsxK3fTx3W1NTE717966qsOQmfZoxJSiEEEJIxahEglJTTJs2Dd26dYOrq6uyQyGEEEKqNUpQFGjo0KHKDoEQQgipEZT+NGNCCCGEkE/RFRQFCgwMBFA0bkZRqyUSQgghtRFdQVGgYcOGwdXVFeHh4coOhRBCCKnWKEFRIJpmTAghhCgGJSgKRAkKIYQoX2BgIObMmVOt9lWVMVcXlKAoCGOMEhRCCFEBr169gre3t8ru68GDB5g3b16F+6npKEFRkMLCQkgkEgCUoBBCah+xhOHeyxScD3mHey9TIJYwpcXi6uqKP//8U2n7/5yYmBgcPXqUV6bqMSsDzeJREOnVE4ASFEJI7XIlLB7rLoYjXpTLlZkLdbB2kCM8ncwrbb8+Pj64ceMGBAIBPDw80LdvXwBAcnIy/P39uafd+/v748yZM5g5cybOnDmDN2/eoH379pg4cSIiIyNx+PBhpKenw9PTEwMGDPi/47pyBXfv3sUPP/zAlYWGhmL79u3Ys2dPiTE9f/4c69evBwDo6OigadOmmDJlCurVqweg6PlyO3fuRGZmJiZNmgQAGDx4MKytrXkxA0XPn/v7778RGRkJExMTjB8/Hvb29tx26XEtWLAA3t7eePv2LVxcXDBhwgSoqVX/6w/V/whUhDRBEQgE0NbWVnI0hBBSNa6ExWO2VxAvOQGABFEuZnsF4UpYfKXsd+PGjZg2bRqMjY1hZWWFbdu24bfffgNQ/HbJixcv8Ndff2HgwIEQCAQwNzfH3LlzMXz4cHzxxRcQCoUwMjLC8OHDcfr0aa5dWFgYLly4wNvvmzdvcPjw4VLjEgqF6N69O7p3747WrVvj1q1bcHJyQlpaGgCgXr16cHR0hKamJlfPxsamWMzp6elwdnbGqVOn4OjoiNjYWLRp0wZ37twpdlx9+/aFWCxGo0aNsHTpUnz33XcV+3BVBF1BURBdXV2sX78eBQUFEAgEyg6HEEIqnVjCsO5iOEq6mcMACACsuxgOD0czqKsp9vfi2bNn8d1332HRokUAgPnz5yMuLq7U+jk5OTh9+jRatWoFAPjw4QN+//13REREcM+CS0hIwMGDB/Hll1+WOy4TExPuyggAzJ49G507d8bu3buxdOlSmJubw93dHSdPnuTVe/nyJa+f3377DRKJBDdv3oSWlhaAoj+AFy1ahAcPHvCO69ixY2jdujUAoH79+li8eDGXrFVnlKAoiIGBAVasWKHsMAghpMo8jEktduXkYwxAvCgXD2NS0cnOWKH7dnJywoEDB+Dg4AB3d3fo6enBwsKi1PrGxsZccgIAdnZ2sLCw4D2o1s7OjltwsyKePn2KU6dO4d27d8jPz0dKSgoiIiLk6uPWrVsYOnQol5wAwJgxY9C/f3/k5ORwQwmMjY255AQA7O3tkZqayqtTXdEtHkIIIeWSlFF6clKeevLYsmULRowYgVWrVsHY2Bienp54+vRpqfV1dHR479XU1EosE4vFFYrr/PnzaNeuHZKSktCqVSt0794d5ubmyMjIkKuf5ORkGBvzk7r69euDMYbk5GSurKRjAFDh41AFlKAoSHp6OoKDgxEdHa3sUAghpEqY6Ot8vpIc9eShp6eH1atX49GjR3j79i2EQiHGjBmj0H1oa2sjPz+fVyYdS1Ka3bt34+uvv8aOHTswd+5cTJo0qVgSIcswAGtra7x69YpXFhMTAw0NjTKvFNUklKAoyIMHD+Ds7ExPNCaE1BrtbYxgLtRBaadbAYpm87S3MVL4vk+cOIHCwkIARbc5OnXqhA8fPih0H02bNkV0dDQ3tqWgoAAHDx4ss41AIOAlMffv34evry+vjpGRETIyMlBQUFBqPyNGjMDJkyfx+vVrAEB+fj62bt2KIUOGQEOjdozOqB1HWQVokTZCSG2jribA2kGOmO0VBAHAGywrTVrWDnJU+ABZAAgKCsJ3332HNm3aoKCgAH5+fti2bZtC99GnTx906tQJrq6ucHNzw+PHj9GwYcMy23z77bcYOHAgXr16BT09Pdy/fx+Ojo68Op06dUL9+vXh7u6OZs2aYfDgwcX6mTx5Mi5duoS2bdvC3d0dYWFhEAgEtWoxN0pQFIQSFEJIbeTpZI6d452LrYNiVsnroPz666+YP38+Hj58CHV1dezfvx8mJiYAii961qVLF/z++++89j169ODqS/Xr1483aFZNTQ2+vr7w8/NDWloaNm7cCA0NDdy4cYOr8+m+evTogefPn+Pu3bvQ1NTEwYMHERUVhby8PK5O3bp18ezZM9y+fRspKSmwsbGBkZERrx81NTWcPn0aQUFBiIiIwJw5c9CtWzfeoNmSjsvOzg779+8vdlupOhIwxpS33J+SpKenQygUQiQSwcDAQCF9HjhwAJMnT0a/fv1w6dIlhfRJCCGVLTc3FzExMbCxsanQSU0sYXgYk4qkjFyY6Bfd1qmMKyekapT1c1EZ59CS0BUUBaErKISQ2kxdTaDwqcSkdqNBsgpCCQohhBCiOJSgKAglKIQQQoji0C0eBencuTOWLVsGFxcXZYdCCCGEVHuUoChIjx490KNHD2WHQQghhNQIdIuHEEIIISqHrqAoyNu3b5GTkwNTU9NKnXZFCCGE1AZ0BUVBFi9ejGbNmuHAgQPKDoUQQgip9ihBURCaxUMIIapNLBZDLBYjPT0dACCRSJQcESkLJSgKQgkKIYSotpEjR8Lf3x87duzAzJkzERISouyQSBkoQVEQSlAIIUS17dmzB1u2bMGePXswYcIEODs7y9Tu1q1bCA0NreTogLi4OOzdu5d7HxERgZMnT1aoT8YYDhw4gMzMzIqGV+UoQVEQSlAIIaRqvXr1Cnv37kV+fn6xbUeOHCmWVERFRSEzMxN6enq4e/cuZHkUXUJCAoYPHw6hUKiwuEsTHh6O6dOnc+9v3ryJ5cuXc+/v37+Pa9euydWnQCCAn58ffvzxR4XFWVUoQVGQ7OxsAJSgEEJIVTl58iSWLVvGe8IvAIhEInz11VeIiorilZ85cwa7du3C/Pnz0aBBAwQFBX12H7/88gsGDBiAxo0bKzJ0mTg4OGDkyJHc+2PHjmHr1q1y97N06VJs2bIFSUlJigyv0tE0YwWhKyiEEFK1Hj16hHbt2hUrDwwMBGOs2LaNGzcCAO8qRVmys7Nx4MABnDt3jivLy8vD4cOHMWzYMERHRyM8PBzt27dH8+bNARRdpQkMDIShoSE6duwIIyMjrq2vry9evXoFgUAAExMTuLi4wNzcvNT9m5qaom3btgCAkJAQPH36FElJSdxtoL59++LWrVtwcXHh9i918OBBdO7cGU2bNkWzZs3QunVr7Nu3D8uWLZPp2FUBXUFRkKlTp2Lu3LmwtLRUdiiEEFIrPHr0qMTHi0gTBFtb2wr17+fnh9zcXHTu3Jkry8jIwPTp0zF8+HBMnToVN27cQGJiIgBg3rx5cHNzw7lz57BlyxY0a9YMvr6+XNvIyEjcv38f9+7dw/bt29GsWbMyx5h8fIvn3bt3iI+PR1paGu7fv4/79+8jLS0N//zzD1avXs1r5+/vjylTpkBbW5sr69GjB/75558KfR5Vja6gKMiqVauUHQIhhCiMdFyHpqYmBAIBgP+bpqumpgYNDQ2F1lVXV5crPpFIhJcvX/KuKEhduHChxCsr8goJCYGNjQ3vRC9lZGQEX19fqKkV/Z1/4MAB+Pj4IDIykrtqsnv3bkycOBGxsbHQ0NDA119/zevD29sbc+bMwdChQ3mfUUkGDBiA69ev48WLF7zjnTZtGvr374+UlBQYGxsDAPbt24devXrB2tqaq9eiRYty3R5SJkpQCCGEFPPLL78AKFqEsk6dOgCAO3fu4MaNG2jbti0GDx7M1f3tt99QUFCA+fPnw9DQEAAQEBCAq1evomXLlhg2bBhXd8uWLcjOzsbs2bNhYmICoCgRkDehCAoKAmMMhYWFuH//Pm/bw4cPsWjRIrmP+VMpKSnc8Xxq2rRpXHICAF5eXmjatCkuXLgAxhgYY8jKykJcXBxevnwJe3t7AEB0dDQeP36MlJQUZGZm4v3794iNjYWdnV25YuzZsyesrKzg5eWF+fPnIysrCydPniyWtBkaGiInJwc5OTnVZigCJSgKwBjDq1evoKurCxMTE94PLSGEEMULCgqCqakp9u/fzytPSkrC33//rZAny+vr65c6PdfMzIz3/vXr1xAKhfD39+eVT506lbs6tHTpUuzcuRNdunThkjMAeP/+fbkTFIFAgClTpmD//v2YP38+Tpw4AS0tLQwZMoRXLzMzE5qamtDR0SnXfpSBEhQFyMrK4u51ZmZmcn9tEEJIdSUd+6CpqcmVubm5oWPHjsX+CFu8eHGxuq6urnB2di5Wd/78+cXqtmnTRu74Hj16xA0g/VhwcDAAKOQWT/PmzREbGyvTLShDQ0O0bdsWu3btKnH7mzdvsHHjRgQGBnKxRUdH4+DBgzJNdy7L5MmTsXbtWgQFBWHfvn0YO3ZssdtS0dHRsLe3526rVQf0p74CSGfwADSLhxBSM2hpaUFLS4t3QlNXV4eWllax8RKKqCuv0hKUkJAQhQyQBYDu3bsjNzdXphVnBwwYgBMnTnADZqVevnwJoOgqCQDeuJCjR4/KFY++vj7vfCNlbm6O/v37Y8WKFdwA2U/5+/ujT58+cu1P2VTiCkpubi4+fPgg1+2R8rSpLNIfGC0tLaXHQgghNV1GRgaeP39eaoKiiKsnANCgQQMMGTIER48e/WyfS5Ysga+vL5ydnTFt2jQYGBggICAAsbGxuHfvHpycnODo6IghQ4Zg1KhRCAsLw4ULF+SKp2PHjti8eTM2bNgAY2Nj9O3bF1ZWVgCKxsQMHjwYbdq0Kfa5JCcn47///sP//vc/+T4AJVPq2VQikXCDquzt7WFhYYFTp04pvE1ly8zKhrZVSxi07Il7L1MgllTsch0hhJDSpaSkYMqUKejUqVOxbY6Ojpg6darC9rV69WocOnQIIpEIAKCjo4OpU6dyM2ak9PT0cPPmTWzZsgUZGRlITEzE8OHDuTEpmpqa8PPzw6BBgxAREYEmTZogICAAU6dO5cajWFpa8mL/dKG2AQMGwMvLC4mJidw0Y6k+ffpAU1OzxGPfuXMnvvjiCzg4OCjsc6kSTIl+/fVXZmxszJ48ecIkEgnbsWMH09DQYGFhYQpt8ymRSMQAMJFIVOFjuPwkjjl/f4k1WurDvTr+7MsuP4mrcN+EEFLZcnJyWHh4OMvJyVF2KCpr27Zt7MaNG8oOo0znzp1jenp6LDU1lVcukUjYokWLWGxsrFz9lfVzochzaFkEjFVwdE4FNGrUCGPGjMGvv/7KlTVt2hQDBgzAH3/8obA2n0pPT4dQKIRIJIKBgUG5478SFo/ZXkH49AOU3lndOd4Znk6lrxJICCHKlpubi5iYGNjY2FSrGR6kSGxsLHx8fLB582aMGDECGzZsUEi/Zf1cKOoc+jlKu8WTlJSE169f81boA4pGiQcEBCisTWURSxjWXQwvlpwA4MrWXQyn2z2EEEIqTVpaGkJDQzFv3jz89NNPyg5HoZQ2SDY5ORkAit3Hq1+/Pu7evauwNkDRsxPy8vK49+np6eWK+WMPY1IRL8otdTsDEC/KxcOYVHSyMy61HiGEEFJebdq0we7du5UdRqVQ2hUU6WyXwsJCXnlBQUGpU87K0wYoWhFRKBRyL+mo54pIyig9OSlPPUIIIYT8H6UlKA0bNgQAJCQk8MoTEhJKfeBeedoARQsOiUQi7vXmzZuKhA4AMNGX7V6trPUIIYQQ8n+UlqDo6+vD2dkZV69e5coKCgrw77//wt3dnStLSUlBXFycXG0+pa2tDQMDA96rotrbGMFcqIPS1uQTADAX6qC9jVEpNQghRHUocb4EUUGq8POg1HVQ1qxZg8OHD2PHjh0ICgrCpEmToK6ujtmzZ3N1li5dylv9TpY2VUFdTYC1gxwBoFiSIn2/dpAj1NWqz7LChJDaR7rkfHZ2tpIjIapE+vPw8SMJqppSV5L94osvcOzYMWzZsgVbt25Fy5Ytcfv2bdSvX5+rU79+fd7tG1naVBVPJ3PsHO+MdRfDeQNmzYQ6WDvIkaYYE0JUnrq6OgwNDZGUlASgaMGx6vS8FqJYjDFkZ2cjKSkJhoaG5XoMgaIodR0UZVH0HG6xhOFhTCqSMnJhol90W4eunBBCqgvGGBISEvDhwwdlh0JUhKGhIczMzEpMVqtqHRSVeBZPdaeuJqCpxISQaksgEMDc3BwmJiYoKChQdjhEyTQ1NZV65USKEhRCCCEAim73qMKJiRBAyYNkCSGEEEJKQgkKIYQQQlROrbzFIx0XrIgl7wkhhJDaRHrurOw5NrUyQcnIyAAAhSx5TwghhNRGKSkpEAqFldZ/rZxmLJFIEBcXB319fYXN909PT4eVlRXevHlTqdOuVFFtPnagdh8/HXvtPHagdh9/bT52ABCJRLC2tkZaWhoMDQ0rbT+18gqKmpoa91wfRVPUUvrVUW0+dqB2Hz8de+08dqB2H39tPnbg/x7gW2n9V2rvhBBCCCHlQAkKIYQQQlQOJSgKoq2tjbVr10JbW1vZoVS52nzsQO0+fjr22nnsQO0+/tp87EDVHX+tHCRLCCGEENVGV1AIIYQQonIoQSGEEEKIyqEEhRBCCCEqp1augyKLjIwMREZGon79+mjcuLHC2pSnX2V49uwZcnNz4eTkBE1NzQq3CQ4ORlZWFq/M3NwcdnZ2CotZUZKTk/Hq1Ss0atQIJiYmMrURiUQICwuDra0tzM3NFdZvVSssLMTTp0+hoaEBR0dHmRcyDA4OBmMMzs7Oxbb5+/sXK7Ozsyv1c1KmN2/eIDExEc2aNZNpfQuxWIyoqChoaGjAxsYGGhol/0qVt19lyMrKQkREBOrVqwdbW1uZ2iQlJSEuLg42NjbFVhTNyspCcHBwsTatWrVSyc8gMjISWVlZaNGihUyDP6XfvUAggK2tLbS0tBTSrzKkpKQgOjoaVlZWMDMzk6lNUlIS3r17B2traxgbGxfbFhUVVaxN586d5Vs7hZFi9uzZw/T09Ji9vT2rU6cO69evH8vMzKxwm/L0W9VevXrFWrVqxerXr88aN27MTE1N2Y0bNyrcpkWLFszGxoa5ublxr02bNlXegZTT4sWLmba2NnN0dGTa2tps3rx5TCKRlFo/JiaGTZ8+nZmZmTF1dXW2bds2hfSrDHfv3mWWlpbMysqKmZiYsObNm7PIyMgy22zdupU5ODgwQ0NDZm9vX2IdAMzJyYn33Z8+fboyDqHccnJy2LBhw5iuri5r3rw509XVZVu3bi2zzU8//cRMTU2Zg4MDa9y4MbO0tGQXL16scL/KcOTIEaavr8+aNWvG9PX1WY8ePdiHDx9Kre/v7886d+7MTExMWJs2bZiuri6bOXMmKyws5OoEBwczAKxDhw687/7x48dVcUgyi4uLY+3atWP16tVjtra2zNjYmF26dKnMNlu2bGEWFhasZcuWzMbGhtWvX58dOnSowv0qw5o1a3i/m6ZOncrEYnGp9YODg5m7uzuzsLBgbdu2Zbq6umzkyJEsOzubq7N//36mpaXF+97d3Nx4dWRBCconHj9+zNTU1Ji3tzdjjLGkpCTWuHFjNm/evAq1KU+/ytClSxfWq1cvlp+fzxhj7Ntvv2XGxsZMJBJVqE2LFi3Y77//XqmxV5SXlxfT1dVljx49YowxFhoayvT09Njff/9dapsrV66wXbt2sYyMDCYUCktMUMrTb1XLzs5mFhYWbM6cOYwxxgoLC1n//v2Zs7Nzme0WLlzInj59ylauXFlmgnL9+nWFx6xIy5YtYw0bNmRxcXGMMcbOnj3LALD79++XWL+wsJCtXLmSpaSkcGVr165lenp6LD4+vtz9KsPz58+ZpqYm2717N2OMsbS0NGZvb88mT55capv9+/ezu3fvcu/Dw8OZoaEh27BhA1cmTVCSk5MrL3gF8PT0ZJ07d2Y5OTmMsaLv0cDAoMy4N23axNLS0rj3mzdvZmpqaiwhIaFC/Va1c+fOMU1NTXbnzh3GGGMRERFMKBSyLVu2lNrm7NmzLCAggHv/+vVrZmxszH788UeubP/+/czS0rLC8VGC8olFixaxJk2a8Mp+/fVXJhQKeX8dyNumPP1WtaioKAaA+fr6cmXv379nGhoa7PDhwxVq06JFC7ZmzRr28OFD7pe1qunZsycbPnw4r2z06NHMzc1NpvalJSgV7bcqnDlzhgkEAt534+/vzwCw4ODgz7b/XIJy6NAhFhgYyDuhqxJTU1P2/fff88qcnJzYzJkzZe4jISGBAWD//POPQvutbGvWrGHm5ua8K3rbt29nOjo6cv3FO2TIEDZo0CDuvTRBuX//PgsODmYZGRkKjVsR3r17xwQCATt37hxXlpmZyXR1ddnOnTtl7icoKIgBYCEhIQrtt7INHjyYeXp68sqmTZvGWrduLVc/7du3Z7NmzeLe79+/n5mbm7MnT56wsLAwlpeXV674aJDsJ4KDg9GuXTteWfv27SESiRAdHV3uNuXpt6pJ7xd/HKexsTFsbW1LvJcsb5vNmzdj+vTpaNasGTp06ICIiAhFH0KFlPYdlXbsyu5XkYKDg2FhYcEbF9K+fXtuW0UtXLgQU6ZMgYWFBYYMGYLk5OQK96kocXFxSExMrPB3FBAQAADcuCpF9VvZgoOD4ezszBtv1L59e+Tm5sr8/2hBQQFCQ0PRpEmTYtuGDx+OMWPGwMjICPPmzUNBQYHCYq+okJAQMMZ431GdOnXg4ODw2e/o7du38Pf3x+nTpzF79myMGjUKrVu3rnC/Vam0301hYWFlfk9isRj+/v7w9fXFsmXL8PbtW8ybN49XJz4+Hl9++SUGDhwIIyMj/P7773LHRwnKJ1JTU4sN+JG+T01NLXeb8vRb1VJTU6Gurl5ssJuxsXGZxy5LmyVLliAlJQUhISF4+/YthEIhhg0bhvz8fMUfSDkwxvDhw4cSv6Ps7Gzk5eWpVL+KVtLPp6amJvT19Sv887l7924kJycjNDQUUVFRiIyMxLRp0yrUpyJJj6+k70jWY3///j3mzZuHkSNHwt7eXmH9VgVF/G5avnw5UlJS8M0333BlQqEQV69exZs3b/Ds2TPcu3cPhw4dwk8//aS44CuoIt/RrVu3sHTpUixatAjJycm8n+nq/t2LxWKkp6eX2i4nJwfLli3D4sWLsW3bNkyYMAFNmzbltjdt2hRhYWGIjIxETEwM9u3bh2+//RZnz56VKz5KUD6hqamJ3NxcXllOTg4AlDpKW5Y25em3qmlqakIsFhfLnHNycso8dlnaTJw4kRvBLhQKsXHjRjx79gyhoaEKPoryEQgE0NDQKPU7knUmU1X1q2gl/XwCQG5uboV/PqdPn879dW5tbY2VK1fi4sWLxWZ1KYv0OyjpO5Ll2EUiETw9PWFmZoa9e/cqrN+qUtHfTZs3b8aOHTtw6tQp3sxEGxsb9OnTh3vfrl07TJ8+HceOHVNM4ApQke9o3LhxuHPnDmJjY7FkyRJ4enri6dOnFe63KpX3u69bty78/f0REhKC8PBweHt747vvvuO2u7m5oUWLFtz7kSNHomfPnnJ/95SgfKJRo0Z49+4dr0z63trautxtytNvVWvUqBGAokvTH4uLiyvz2OVtAwCmpqYAUOwzUSZra+sSv6OGDRtW6LHildWvIjVq1AgJCQmQSCRcWVJSEgoKChT+82lqagrGGOLj4xXab3lZWVlBTU2txO/oc8eenp6OPn36QF1dHVeuXIG+vr5C+q1KFfnd9Mcff2DVqlU4e/YsPDw8PrsvU1NTlfp/Xvr7q6Lf0axZs6Cnp4fr168rtN/KVtp3b2hoyPtZ/lwfI0eOxOXLl8usV57vXjV+O6oQDw8P3L59GyKRiCs7f/482rZty10KE4lE8Pf35/4ClKWNLHWUrVOnTqhTpw4uXLjAld27dw9JSUm8Xz5BQUF4+fKlzG2ys7PBPnnk07Vr1wCAl2Urm4eHB3x8fLhYGWO4cOEC79jj4uJKXNejov0qm4eHB9LT03Hz5k2u7Pz589DS0kK3bt24Mn9/f7kSi5Kukly7dg36+voq84taT08PnTt35v0MZ2VlwdfXl/cdvXjxgjd+QJqcAEXH9OltTln7VTYPDw88ePAASUlJXNn58+fRtGlT7kSbmZkJf39/3mX/rVu3Yvny5Thz5gw8PT2L9VvSd3/9+nU4OTlVwlGUT7t27WBkZMT7jp48eYKYmBjedyS9PQkUXWH49PdZfHw8MjMzud/lsvarbB4eHrh06RLEYjFXdv78eV6MiYmJ8Pf35/54Kel7ffHiBe889mmdnJwc+Pv7y//dl2tobQ2Wk5PDHB0dWbdu3di5c+fYmjVrmLq6Ort8+TJX58aNG7zZDbK0kaWOKtiwYQOrU6cO27lzJzt27Bizs7NjQ4cO5dWxt7fnzUL4XJu7d+8yNzc3tmfPHnb16lW2fv16VrduXTZ79uwqOy5ZxMTEsHr16rFx48axCxcusIkTJzIDAwP2/Plzrs7vv//OPv7fJiMjg/n5+TE/Pz9Wt25dtnDhQubn58eePXsmV7+q4KuvvmLW1tbsyJEjbM+ePUwoFLI1a9Zw2wsKChgA3kyl0NBQ5ufnxyZMmMCsra25z0I65Xzv3r1syJAh7NChQ+zSpUts4cKFTENDg23fvr3Kj68sN2/eZJqammzZsmXs/PnzrHfv3szOzo4382Tq1KmsRYsWjDHG8vPzuXVAfHx8uOP28/PjTTWVpV9lKygoYM7Ozqxjx47szJkzbP369UxdXZ2dOnWKqxMQEMAAMD8/P8ZY0fcKgK1cuZJ37NJZLIwxtmDBAjZz5kx28uRJdv78eTZq1Cimra3N/v333yo/xrLs2LGD6ejosC1btrATJ04wBwcH1qdPH16ddu3asXHjxjHGimbsdOrUie3atYtdv36d7du3jzk5ObGWLVvyZj3J0q+yxcXFMRMTEzZ8+HB24cIFNnPmTKanp8eePHnC1dmzZw8DwP3Menp6sjVr1rCLFy+y8+fPs2nTpjENDQ3eGi/9+/dnK1euZBcvXmTHjx9nbm5uzMTEhEVHR8sVHz3NuATv37/Hhg0bEBwcDGNjY8yaNQs9evTgtgcHB2PevHk4cOAAN2r9c21kraMKjhw5guPHjyMvLw89e/bEggULeCsgjhs3Dm3btsXixYtlbhMcHIw9e/bgxYsXaNiwIYYNG4aBAwdW6XHJIioqCps2bcLLly9hY2ODxYsXw8HBgdt+8uRJbNmyhbuKEhkZialTpxbrp0ePHvjxxx9l7lcVFBQUYPv27bh69So0NDTw5ZdfYtKkSdz4EbFYDHd3dyxatAjDhg0DAMycOZO77/6xCxcuwMjICABw9epVHD16FAkJCbCxscGUKVPg6upadQcmozt37mDHjh1ITExEy5YtsWzZMt6qmhs2bMDTp09x6NAhiEQiDBgwoMR+li1bxvvZ/ly/quDDhw/YsGEDAgICUK9ePUybNg19+/bltkt/znfu3ImWLVti9erVuHHjRrF+mjdvzo3DkUgk8PLywj///IOsrCw0b94cc+fOVckVtE+ePIkjR44gOzsb3bp1w6JFi6Cnp8dtnzZtGho1aoTVq1cDAMLCwrBr1y5uVfCuXbti8uTJ0NHRkatfVRATE4ONGzciMjIS1tbWWLhwITcbCQB8fHzw66+/4vr169DV1UV2djb++usv3L59GxKJBM2aNcPMmTN5g2SldW7dugV1dXW0bt0a33zzDerVqydXbJSgEEIIIUTl0BgUQgghhKgcSlAIIYQQonIoQSGEEEKIyqEEhRBCCCEqhxIUQgghhKgcSlAIIYQQonIoQSGEEEKIyqEEhRAVMWDAAO55FvHx8XBxccHbt28rfb8vXrzAV199hY4dO2Lp0qWVvr+P5efnw8XFRWUeGqlqli5dCm9vb2WHIZOtW7di06ZNyg6D1CCUoBCiIkJDQ5GSkgIAyMvLw6NHj0p8wrCiTZgwAQKBAFu2bMHMmTMrbT/79u3D5MmTeWUSiQSPHj1CRkZGpe23OunTpw/3vJ9Hjx7hwIEDGDx4sJKjks2XX36J9evX49WrV8oOhdQQlKAQooIsLCwQEBAAKyurSt9XUFAQxo0bhw4dOsDW1rbS9hMXF4dnz57xyrS1tREQEIA2bdpU2n6ri9DQUNy/f597gOaGDRvw1VdfoW7dugrp38vLi/d4CkWztLSEh4cHtmzZUmn7ILULJSikRvL29kafPn3Qv39/bNq0CVu2bME333zDbZ8wYQL27NmDn376Cd27d+d+caenp2PVqlXo1q0bevfujf/9738oLCzk2i1duhS//PILb187d+7kXXmYMGEC9u7diw0bNqBfv34YOHAgTp06xWuTlZWFZcuWoWvXrhg7diz+++8/3vaUlBTMmjULycnJAIBbt26hW7duCAwMxMSJE9G9e3fMmzePu+IidfXqVQwcOBAeHh5YvXo1vLy8MHz48BI/o2fPnsHFxQV5eXn45ptv4OLigv/++w8+Pj4YMGAAbt26haFDh6JDhw5ISkrC9u3b4eLiAhcXF3h4eGD58uX48OEDr8/c3Fz89ttv6Nu3L/r27Yvdu3eDMYYzZ85g586dCAsL4/rw8fFBQUEBZs2axT0dGwASEhKwYMECdOnSBf369cO+fft4+/jmm2/w+++/Y9u2bRg0aBD69euHAwcOlHiMH5NIJNi6dSt69OiBwYMH4+7du5g5cyZ279792baykkgkcHNzw7Vr1/D999/D3d0do0ePxqtXryASibB8+XK4u7tjwoQJSEhI4LX18fGBh4cHtLS0IBKJcPbsWYwYMYJX559//kH//v3x5MkTTJkyBV27dsU333xT4hNmP5WQkIAXL17IdBx79uzB1KlTcePGDQwfPhzu7u7466+/AABXrlzB0KFD0bNnT3h5efHaDR8+HIcPH+aefEtIhVT4cYiEqBhvb2+mra3NNm/ezHx9fdno0aOZtrY269u3L1enXbt2TFtbm82fP5/dunWLvXjxgjHGWJcuXZizszO7ePEiO3bsGLO0tGRTp07l2n355Ze8JzkzxtjKlSuZu7s7r28dHR22ZMkSduPGDfbbb78xgUDA/P39uTpDhgxhDg4O7Ny5c+zkyZPM1taWaWhosMOHDzPGip6ADIB74vHZs2eZmpoaa9WqFTt16hS7du0aa9euHevXrx/Xp5+fH9PQ0GCrVq1i//77L5s3bx7T1tZm9vb2JX5OWVlZ3FNqd+zYwQICAlhqairbv38/09DQYK1atWJnzpxhAQEBLD8/n719+5YFBASwgIAAdvnyZda3b1/Wvn17JpFIGGNFT8Xt1q0ba968OfP29mZXr15lM2bMYAcOHGBJSUls9uzZzMnJievj/fv3LCcnh/eU3JycHGZnZ8c8PDzYlStXuKcq//jjj1zcffv2Zdra2mzWrFnsv//+Yzt37mRaWlrs/Pnzpf5MiMViNmjQIObs7MxOnz7Nzpw5w5o3b87q1avHjh8/Xqz+jBkzWLt27cp8FRYWFmsXGRnJALBWrVqxv/76i/n6+rIWLVowNzc31r17d7Zv3z527do11qRJEzZp0iRe244dO7J9+/Yxxhjz8fFhmpqaLC8vj1fnu+++Yw0aNGAeHh7swoUL7Pz588zU1JT98MMPpR671KZNm9gXX3zx2XqMMTZixAhmbm7OxowZw65fv87WrVvHALApU6awSZMmMV9fX7Zs2TKmrq7OYmNjuXaxsbEMAAsKCpJpP4SUhRIUUuM0b96cLV++nHsvFouZnZ1dsQSld+/evHYXL15kmpqa7M2bN1zZv//+y9TU1FhUVBRjTPYE5dMTgbu7O1u6dCljjLHg4GAmEAhYWFgYt/3WrVsMQJkJCgDeY9CvX7/O1NTUuJOYp6cnGzNmDG+/3bt3LzVBkQLArl+/zr3fv38/A8AiIiLKbJednc20tbXZw4cPGWOMeXl5MR0dHd7nxxhjubm5jDHGfvzxR9ahQwfetk8TlG3btjFjY2OWlZXF1fn7779ZnTp1mEgkYowVJShubm68fr788kteIvmp3bt3MyMjI5aWlsaV/fHHHwwAe/bsWbH6ERERXCJV2qskJ06cYADYnTt3uLJffvmFqaurs+DgYK5s5cqVrHPnztz7pKQkpqGhwRISEhhjjG3dupVZWloW679v377MycmJ+0wZY2zChAlswoQJpR67lDwJSrNmzVj//v259xKJhGlqarKRI0dyZdLv7r///uPKCgsLmUAgYCdPnpRpP4SURUMJF20IqTQ5OTmIjIxEr169uDI1NTV079692IyYTp068d4/evQIjo6OaNiwIVfWo0cPaGpqIjg4mPc48c9xcXHhvbe0tERiYiIAIDAwEKamptxYAwDo2rVrsUe1f0pbWxtOTk68PiUSCZKTk2FpaYmgoCD8/PPPvDa9evUqdhleFvXq1YO9vT2vLCMjAzt27IC/vz+SkpIgkUggkUgQHR0NV1dX3L59Gy4uLrzPTxq3rB49eoQuXbrwHknv6emJrKwsREREoH379gBK/nyjo6NL7Xfbtm1YtGgRDA0NubIGDRpAV1e3xO/102OXVUhICFxcXNC5c2euLDExEV27duWNs0lMTISlpSX3/tKlS3B2doapqSmAokHSWlpaxfoPDQ3FihUreJ/p69ev0bVr12J1Hz58iDlz5vD2mZmZyfvsZsyYgRkzZvDaZWVl4cWLF9i6dStXlpqaioKCAt5tUunP88fHoa6uDg0NjSoZ3E1qPkpQSI2SmZkJxhjvBAcAderUKVb30zoZGRnF6gkEAujp6ck9y0RDg/+/lkAgAGOMi/HTfQsEAujq6srdJ4Ay+y3puGXxaT8AMGrUKKSmpmLRokVo2LAhtLS00KdPH+Tk5AAoSg4rOqCzpO9A2ufH30FZn29JfT558gTbtm3jlcfExMDJyQnq6urF2sycOROPHj0qM9YHDx4UaxsaGoouXbrwykJCQoqVPX78GF988QX33sfHBwMHDuTem5qa4v3797w2SUlJSEhIKNZXaGgo5s6dWyw+BwcHbtwIUDRINiQkBL/99htXZmFhUazdkydPwBjjJVkhISHQ0tKCq6sr7xjq1KmDJk2acGXp6ekoKCiAmZlZsX4JkRclKKRGadCgAfT19REZGcm7QhIREcGd0EvTpEkTHDp0CIWFhdwJMD4+HmlpabCzswNQdMLPzs7mtYuPj5crRltbW7x79w6ZmZncyTcxMRFpaWly9VNSv5GRkbyyiIiICvUplZOTg6tXr+Lu3bvo0KEDgKK/qj8eJNusWTPcvHkTYrG4xJO+mppaqUmEVJMmTXD16lVe2dOnTwGA+w7KEzvAT2oYY/D29uadhD+2aNGizyalJR1jSEhIsYGtnyYQEokEYWFhWLt2LQCgoKAA165dw4oVK7g6rq6uyMjIQHR0NDezKjQ0FFpaWrwrb7GxsUhLS0Pr1q2LxaKvr8+7WnLz5k28evWq2NWnko7Bzs4O+vr6vDJHR0feVZ2QkBC0atUKamr/N9ciODgY6urqvESGkPKiWTykxpk4cSI2b97MzXC5ceMGfH19P9tuxIgRyM/Px8aNGwEAYrEYS5cuhYODA/dXa+vWrfHff/8hNTUVQNFfkSdPnpQrvr59+6J+/frcCUoikWDZsmVy9VGSSZMmYffu3dw6FGFhYTh69GiF+wWKbtPo6ekhICAAQNFJdcGCBbyEY8KECUhNTcWaNWu4WRzBwcHcZ29mZoa4uLgyZ3hMmjQJYWFh3OJkOTk5WL16NTw8PNC4ceNyxW5iYgJTU1Pus5BIJFizZg3Cw8NLPLEDRbd4pLONSnt9KiUlBe/evUPbtm25MmkC8XFZVFQUsrOzuVs+t2/fhr6+Pq9O8+bN4eDggGvXrnFloaGhaNGiRbEkoW7duuVO3koSGhrKi0W6n0/LQkNDi00P9/X1RY8ePSAUChUWD6m9KEEhNc4PP/wAQ0NDWFtbo2nTppgxYwZ69+5d7LbAp4yNjeHt7Y2tW7fCysoKZmZmePDgAby9vbm2M2fOhK2tLWxsbNC8eXOMGTMG3bt3lys+bW1tHD58GIcPH4aFhQXMzc2Rk5ODBg0alPeQAQBz585Fp06dYG9vD3t7e/Tt2xceHh6fPW5ZqKmpYfv27fjuu+9gb28PU1NTpKeno379+lwda2trXLhwAd7e3qhfvz7s7Owwffp02NjYAAC++OILaGtro1GjRtw04081b94cu3fvxpw5c2BrawszMzOkp6djz549FYp/x44d2LdvHxo3bowmTZogLy8PGhoa3NUgRQgJCYG2tjYcHBx4ZUKhkPsMgKITu4mJCXd7RTqt+1Pz58/nTbEuLXFo1arVZ68OynscsiYoH5eJxWIcPnyYN06FkApR3vhcQirXixcvuCmQQ4cOZdOnT+e2hYeHs/j4+BLbFRQUsKdPn7KoqChuCu2nYmNjub7fvXvHm/FSUt/R0dEsJiaGV5abm8uePHnCPnz4wBhjLDQ0lKWkpDDGGMvLy2MBAQHcbI20tLRiUzdzcnK4KcAfe/36NXvx4gUTi8Vs4cKFzMPDo8RjkAoICGDp6enc+/fv37PQ0NAS66anp7PHjx9zs01CQ0PZ+/fveXUkEgmLioriTT+VKigoYM+fP2eBgYHs/fv3TCKRsICAAJaRkVHiZ/PpZ8YYY1FRUcVmCr1+/Zqb8VSajIwM9uTJE5aZmcmOHj3KLC0tS/1+y+P9+/fs8ePHvLKkpCTebC3GGEtISGDh4eHc+6ZNm7ILFy4U66+goIA5Ojqyy5cvM8aKpjC/e/eOV+fVq1fs5cuXMsWXkJDATacvy8c/h1LBwcG8GVCMMfbo0SNuZhVjjO3bt483m42QiqIEhdQ4L1++ZFeuXOHe+/r6Mg0NDebj46PEqCpfRkYGN02ZsaITjVAoZNu3b1diVMr38uVLdvbsWe59QEAAMzU1Zfv371daTFIRERFMR0eHN636Y3FxcezVq1dVHFX5PH/+nCUnJys7DFKDCBj7zKg1QqqZjIwMjB49GsHBwdDW1kZiYiLWrl1b5Q/Cq2pisRhTp07FP//8A0NDQ7x+/Rpff/01Nm3aVOKAztoiPj4eAwYMQHx8PHR0dJCeno5169aVOPOlqqWlpSE+Ph6Ojo7KDoUQlUMJCqmxUlJSkJKSAhsbG2hqaio7nCqTnp6OuLg4NG7c+LNrq9Qmb9++RWFhIaytrXkzTwghqokSFEIIIYSoHPozghBCCCEqhxIUQgghhKgcSlAIIYQQonIoQSGEEEKIyqEEhRBCCCEqhxIUQgghhKgcSlAIIYQQonIoQSGEEEKIyqEEhRBCCCEq5/8BlfNyN2hL9ZwAAAAASUVORK5CYII=",
"text/plain": [
"<Figure size 600x340 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"m= 0 g=0.000 simulated 0.0000 closed form 0.0000\n",
"m= 1 g=0.010 simulated 0.5855 closed form 0.6196\n",
"m= 2 g=0.020 simulated 0.7443 closed form 0.7438\n",
"m= 5 g=0.048 simulated 0.8485 closed form 0.8455\n",
"m= 10 g=0.091 simulated 0.8905 closed form 0.8859\n",
"m= 20 g=0.167 simulated 0.9057 closed form 0.9077\n",
"m= 50 g=0.333 simulated 0.9225 closed form 0.9215\n"
]
}
],
"source": [
"def H_eq(m, n=N_SAMPLES):\n",
" return H_star * m * (2 * n + m - 1) / (n + 2 * n * m + m * m) if m > 0 else 0.0\n",
"\n",
"ms = [0, 1, 2, 5, 10, 20, 50]\n",
"sim = []\n",
"for m in ms:\n",
" tail = [run_lineage(cfg(m=m, generations=800), seed=100 + s)\n",
" .heterozygosity.to_numpy()[-200:].mean() for s in range(12)]\n",
" sim.append(np.mean(tail))\n",
"\n",
"g = np.array(ms) / (N_SAMPLES + np.array(ms))\n",
"grid = np.linspace(0, 50, 300)\n",
"\n",
"fig, ax = plt.subplots(figsize=(6, 3.4))\n",
"ax.plot(grid / (N_SAMPLES + grid), [H_eq(m) for m in grid], 'k--', label='closed form')\n",
"ax.plot(g, sim, 'o', label='simulation')\n",
"ax.axhline(H_star, ls=':', color='gray', label='$H^*$ (reality)')\n",
"ax.set(xlabel='grounding fraction $g = m/(n+m)$', ylabel='stationary $H$')\n",
"ax.legend(); plt.show()\n",
"\n",
"for m, s in zip(ms, sim):\n",
" print(f'm={m:3d} g={m/(N_SAMPLES+m):.3f} simulated {s:.4f} closed form {H_eq(m):.4f}')"
]
},
{
"cell_type": "markdown",
"id": "3a30c1b0",
"metadata": {},
"source": [
"## The operational threshold, and the floor underneath it\n",
"\n",
"The equilibrium is **smooth** in `g` — there is no phase transition in aggregate diversity, and\n",
"the paper is careful to say so. What there *is* is an operational threshold: the grounding\n",
"fraction that retains a chosen share of reality's diversity. A few percent buys most of it.\n",
"\n",
"But averages hide the thing that matters. A capability of rarity `p` appears in a real-data batch\n",
"of size `m` with probability $1 - e^{-mp}$, so protecting the rarest knowledge is priced\n",
"*per item* at $m \\sim 1/p$ — the floor no amount of average-case grounding removes."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "180153fd",
"metadata": {
"execution": {
"iopub.execute_input": "2026-09-07T14:39:44.884947Z",
"iopub.status.busy": "2026-09-07T14:39:44.884827Z",
"iopub.status.idle": "2026-09-07T14:39:44.887903Z",
"shell.execute_reply": "2026-09-07T14:39:44.887436Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"first swept point retaining 95% of H*: g = 0.091 (m = 10)\n",
"the paper reports g ~ 0.05 for its tested setting; see results/E2 and Fig. 2A\n",
"\n",
"common item (p=1e-02): seen once per batch at m ~ 100 -> P(seen | m=1000) = 100.0%\n",
"rare item (p=1e-03): seen once per batch at m ~ 1,000 -> P(seen | m=1000) = 63.2%\n",
"very rare item (p=1e-04): seen once per batch at m ~ 10,000 -> P(seen | m=1000) = 9.5%\n"
]
}
],
"source": [
"target = 0.95\n",
"i = int(np.argmax(np.array(sim) >= target * H_star))\n",
"print(f'first swept point retaining {target:.0%} of H*: g = {g[i]:.3f} (m = {ms[i]})')\n",
"print('the paper reports g ~ 0.05 for its tested setting; see results/E2 and Fig. 2A\\n')\n",
"\n",
"for p, label in [(1e-2, 'common'), (1e-3, 'rare'), (1e-4, 'very rare')]:\n",
" need = 1 / p\n",
" print(f'{label:10s} item (p={p:.0e}): seen once per batch at m ~ {need:,.0f}'\n",
" f' -> P(seen | m=1000) = {1 - np.exp(-1000 * p):.1%}')"
]
},
{
"cell_type": "markdown",
"id": "2fbabda7",
"metadata": {},
"source": [
"## Where to go next\n",
"\n",
"- `02_paper_figures.ipynb` — regenerate every manuscript figure from the committed artifacts.\n",
"- `configs/layer1/E2.yaml` — the full grounding sweep this notebook miniaturises (Fig. 2A).\n",
"- `REPRODUCING.md` — the map from each paper panel to its config and seed.\n",
"- `tests/test_scientific_validation.py` — these same identities as assertions; if they fail, the\n",
" science is wrong, not just the code."
]
}
],
"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.15"
}
},
"nbformat": 4,
"nbformat_minor": 5
}