Skip to content

Instantly share code, notes, and snippets.

@tomonari-masada

tomonari-masada/zero_inflated_negative_binomial.ipynb

Created December 3, 2025 15:05
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save tomonari-masada/e5ce480092acb275f6ff3f653d111c0d to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save tomonari-masada/e5ce480092acb275f6ff3f653d111c0d to your computer and use it in GitHub Desktop.
Download ZIP
zero_inflated_negative_binomial.ipynb
Display the source blob
Display the rendered blob
Raw
zero_inflated_negative_binomial.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"id": "view-in-github",
"colab_type": "text"
},
"source": [
"<a href=\"https://colab.research.google.com/gist/tomonari-masada/e5ce480092acb275f6ff3f653d111c0d/zero_inflated_negative_binomial.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
]
},
{
"cell_type": "markdown",
"id": "12788303",
"metadata": {
"id": "12788303"
},
"source": [
"# Zero-Inflated negative binomial distribution"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e18fe42a",
"metadata": {
"id": "e18fe42a"
},
"outputs": [],
"source": [
"!pip install numpyro"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "18ee6419",
"metadata": {
"id": "18ee6419"
},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import pandas as pd\n",
"import seaborn as sns\n",
"from scipy.stats import nbinom\n",
"\n",
"import jax.numpy as jnp\n",
"from jax import random\n",
"from jax.scipy.special import expit\n",
"import numpyro\n",
"import numpyro.distributions as dist\n",
"from numpyro.infer import NUTS, MCMC\n",
"\n",
"%config InlineBackend.figure_format = 'retina'\n",
"\n",
"rng = np.random.default_rng()\n",
"rng_key = random.PRNGKey(0)\n",
"\n",
"numpyro.set_platform(\"cpu\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ec802eff",
"metadata": {
"id": "ec802eff"
},
"outputs": [],
"source": [
"size = 1000\n",
"\n",
"true_proportion_of_zero_inflation = 0.3\n",
"number_of_successes = 5\n",
"success_probability = 0.4\n",
"\n",
"negative_binomial_samples = nbinom.rvs(\n",
" n=number_of_successes,\n",
" p=success_probability,\n",
" random_state=rng,\n",
" size=size\n",
")\n",
"negative_binomial_samples = jnp.array(negative_binomial_samples)\n",
"\n",
"rng_key, rng_key_data = random.split(rng_key)\n",
"is_zero_inflated = random.bernoulli(rng_key_data, p=true_proportion_of_zero_inflation, shape=(size,))\n",
"data = jnp.where(is_zero_inflated, 0, negative_binomial_samples)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5d9ca6a9",
"metadata": {
"id": "5d9ca6a9",
"outputId": "e0b0f0b5-84c9-494c-e078-40cf329408f7"
},
"outputs": [
{
"data": {
"text/plain": [
"Array([ 0, 0, 2, 5, 0, 0, 0, 4, 0, 3, 5, 1, 0, 0, 4, 15, 0,\n",
" 9, 0, 0, 0, 4, 10, 9, 0, 0, 3, 1, 3, 3, 0, 0, 6, 10,\n",
" 0, 0, 0, 0, 5, 0, 0, 0, 12, 0, 2, 0, 6, 5, 8, 8, 9,\n",
" 0, 6, 0, 8, 0, 10, 5, 0, 5, 10, 11, 10, 13, 8, 0, 24, 0,\n",
" 0, 0, 0, 8, 6, 3, 1, 14, 4, 12, 10, 5, 0, 7, 0, 5, 0,\n",
" 4, 0, 1, 0, 4, 0, 7, 17, 10, 9, 10, 0, 4, 0, 5, 5, 10,\n",
" 0, 7, 0, 3, 9, 3, 6, 0, 0, 11, 0, 2, 0, 5, 7, 13, 15,\n",
" 0, 0, 2, 13, 8, 3, 0, 2, 7, 12, 6, 0, 5, 0, 2, 7, 9,\n",
" 0, 10, 5, 1, 12, 12, 19, 15, 0, 0, 0, 15, 2, 9, 0, 8, 0,\n",
" 10, 9, 7, 3, 6, 3, 11, 5, 7, 0, 3, 9, 0, 0, 9, 4, 0,\n",
" 5, 1, 0, 6, 3, 5, 5, 0, 0, 9, 0, 5, 7, 2, 8, 6, 8,\n",
" 6, 0, 11, 3, 6, 7, 7, 0, 10, 5, 6, 9, 17, 0, 0, 9, 0,\n",
" 0, 5, 4, 0, 11, 0, 0, 2, 9, 6, 9, 4, 11, 0, 0, 0, 8,\n",
" 5, 7, 22, 8, 6, 0, 5, 0, 6, 5, 3, 4, 6, 0, 0, 2, 1,\n",
" 8, 0, 2, 0, 6, 10, 3, 0, 5, 6, 8, 8, 0, 3, 8, 5, 6,\n",
" 1, 9, 2, 9, 0, 4, 14, 4, 0, 0, 0, 0, 2, 13, 0, 0, 9,\n",
" 0, 4, 0, 0, 0, 4, 0, 0, 5, 8, 8, 0, 3, 0, 6, 3, 0,\n",
" 8, 7, 0, 0, 9, 8, 4, 7, 2, 4, 0, 10, 2, 0, 8, 6, 1,\n",
" 7, 8, 0, 0, 0, 0, 0, 13, 5, 8, 7, 0, 0, 0, 5, 1, 4,\n",
" 6, 9, 11, 6, 0, 4, 9, 4, 12, 10, 0, 0, 4, 0, 13, 12, 0,\n",
" 0, 15, 0, 10, 14, 4, 6, 4, 7, 6, 7, 0, 9, 0, 10, 3, 12,\n",
" 7, 4, 0, 0, 8, 10, 10, 0, 0, 0, 2, 4, 5, 13, 0, 1, 5,\n",
" 5, 5, 0, 10, 6, 13, 3, 15, 10, 11, 11, 0, 0, 2, 10, 4, 8,\n",
" 5, 0, 3, 7, 4, 0, 5, 3, 23, 0, 10, 4, 0, 1, 0, 4, 3,\n",
" 2, 5, 6, 0, 2, 0, 3, 0, 0, 8, 0, 0, 0, 9, 18, 7, 14,\n",
" 17, 3, 6, 6, 8, 3, 10, 16, 0, 8, 0, 0, 5, 0, 7, 6, 0,\n",
" 0, 9, 14, 0, 4, 7, 5, 0, 8, 0, 0, 0, 17, 10, 0, 5, 8,\n",
" 5, 5, 6, 9, 0, 2, 5, 0, 3, 9, 0, 0, 1, 0, 0, 0, 0,\n",
" 16, 14, 0, 0, 6, 5, 2, 15, 0, 0, 0, 0, 0, 0, 7, 3, 4,\n",
" 0, 6, 0, 13, 9, 0, 11, 3, 0, 8, 0, 0, 8, 7, 1, 13, 2,\n",
" 11, 4, 2, 1, 12, 11, 5, 2, 0, 11, 5, 7, 0, 14, 7, 11, 0,\n",
" 6, 0, 0, 0, 0, 6, 2, 5, 3, 8, 0, 0, 0, 4, 14, 5, 0,\n",
" 13, 4, 9, 5, 4, 0, 8, 6, 11, 0, 4, 7, 0, 9, 6, 0, 4,\n",
" 5, 0, 6, 6, 6, 4, 6, 1, 4, 2, 4, 7, 0, 15, 11, 0, 8,\n",
" 0, 4, 3, 18, 2, 6, 4, 3, 0, 4, 4, 4, 0, 0, 13, 0, 11,\n",
" 0, 0, 13, 0, 3, 0, 13, 8, 10, 0, 7, 0, 8, 9, 0, 3, 6,\n",
" 10, 0, 7, 8, 3, 0, 3, 15, 9, 0, 0, 2, 8, 6, 10, 10, 3,\n",
" 20, 0, 11, 13, 9, 2, 12, 0, 0, 13, 0, 0, 3, 13, 4, 5, 13,\n",
" 0, 0, 8, 7, 0, 4, 3, 0, 6, 8, 0, 12, 11, 9, 6, 0, 9,\n",
" 0, 0, 8, 8, 10, 4, 0, 0, 5, 5, 3, 8, 2, 12, 0, 0, 9,\n",
" 6, 15, 8, 0, 0, 13, 9, 4, 0, 8, 0, 2, 13, 12, 0, 2, 0,\n",
" 7, 12, 4, 17, 8, 0, 0, 16, 14, 3, 0, 8, 6, 1, 0, 0, 6,\n",
" 21, 4, 7, 0, 0, 11, 0, 3, 7, 9, 3, 3, 14, 7, 9, 3, 0,\n",
" 11, 2, 0, 0, 6, 7, 0, 8, 11, 14, 0, 13, 0, 0, 16, 7, 0,\n",
" 0, 2, 5, 0, 8, 5, 9, 12, 0, 17, 4, 11, 3, 7, 3, 2, 10,\n",
" 7, 0, 0, 0, 0, 9, 7, 5, 13, 14, 3, 0, 3, 4, 17, 5, 0,\n",
" 15, 0, 3, 0, 0, 7, 5, 9, 7, 0, 18, 4, 4, 0, 3, 13, 9,\n",
" 0, 16, 0, 7, 8, 0, 8, 6, 5, 12, 0, 3, 4, 5, 5, 0, 0,\n",
" 5, 0, 5, 10, 1, 9, 4, 11, 9, 6, 2, 0, 0, 5, 0, 2, 9,\n",
" 0, 2, 0, 3, 11, 6, 4, 0, 0, 4, 3, 6, 14, 0, 4, 1, 1,\n",
" 0, 6, 9, 0, 0, 0, 12, 0, 7, 0, 10, 0, 3, 7, 8, 8, 12,\n",
" 4, 0, 2, 15, 0, 0, 4, 9, 7, 3, 5, 10, 0, 10, 10, 0, 3,\n",
" 9, 12, 6, 20, 8, 0, 8, 0, 3, 0, 0, 8, 6, 6, 6, 8, 7,\n",
" 7, 4, 0, 4, 4, 9, 6, 14, 0, 8, 2, 2, 0, 8, 17, 11, 0,\n",
" 8, 10, 2, 6, 3, 2, 8, 2, 7, 0, 0, 8, 0, 6, 0, 13, 6,\n",
" 2, 0, 0, 0, 8, 14, 5, 0, 0, 17, 0, 10, 0, 0, 0, 14, 4,\n",
" 0, 8, 4, 1, 0, 12, 5, 0, 7, 0, 10, 6, 7, 0, 11, 0, 6,\n",
" 3, 6, 15, 4, 5, 2, 5, 0, 2, 0, 11, 0, 0, 1, 11, 6, 4,\n",
" 3, 4, 3, 8, 12, 8, 8, 3, 0, 4, 7, 5, 0, 6], dtype=int32)"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c82450fb",
"metadata": {
"id": "c82450fb",
"outputId": "e55090a7-4fb0-4ac8-c736-b631696a5770"
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABHcAAAOMCAYAAADHXV9PAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAewgAAHsIBbtB1PgAAihhJREFUeJzs3QeYVNX9P/4DIk3siA07Yo8x2LBhiRqxa2JJDDbU2DUGS6xRE6No7DWiSIyaaCyxRWNiQ7GgxoqiKEbU2BULRWT+z+d+fzP/AbayOzN7l9freebZuzt37pydcu/M+57zOR0KhUIhAQAAAJBLHWvdAAAAAABmn3AHAAAAIMeEOwAAAAA5JtwBAAAAyDHhDgAAAECOCXcAAAAAcky4AwAAAJBjwh0AAACAHBPuAAAAAOSYcAcAAAAgx4Q7AAAAADkm3AEAAADIMeEOAAAAQI4JdwAAAAByTLgDAAAAkGPCHQAAAIAcE+4AAAAA5JhwBwAAACDHhDsAAAAAOSbcAaBNWXbZZVOHDh2yy/jx42vdHGixadOmpT/+8Y9pyy23TIsuumjq3Llz6TW+zz77pLbA+27Oc9ppp5We81iutnjtF+9/+PDhVb9/gPZGuANAk2266aaz/WWg/ItEbAfmBFOmTMlCnQMPPDA98MAD6cMPP0zffvttrZsFALQznWrdAACohYceeihtttlm2fKAAQOy36G1nXvuuTO8tuK11qdPn9S1a9fs9/XXX7+GrYM5S/RKW2655bLlZZZZRi81oF0R7gAAVMif/vSn0vJ1112XBg0aVNP2AADtk3AHgDbFmVTai2+++Sa99tpr2XLU2dlrr71q3SQAoJ1ScwcAoAI+++yz0nIUUu7Y0ccuAKAyfMoAAKiA8sLJgh0AoJJ80gCgTWnqlMxfffVVuuKKK9K2226bll566dS9e/c099xzp/nnnz+tvPLKafvtt0+/+93v0ksvvVTnrF3FYsrh4YcfLt1n+SXa0tAX92uvvTbttNNOWWHObt26pfnmmy+ttNJKaf/990///Oc/m/2/33zzzVm7l1xyydSlS5fUu3fvbKalqNUS02k3dfrgutb5/PPP04UXXpg22WSTbPudOnXKro+/l4vZnOL/2nvvvdNaa62VFlpooexxXWCBBbLHdd9990333XffbE+1PHny5HTllVdmM6Ytvvji2XCl+D+jFs0rr7xS5/N86aWXpo022ihbPwoRr7DCCunQQw9NEyZMSJVQKBSy52LPPffM7qtHjx7ZJZZ/+tOfpltuuSVbpz7F/7lYuDW8/fbbzXp9NVUlXodFr776ajrqqKPSqquumm0zLt/73vfSSSedlP73v/81uX3XX3992mWXXdLyyy+fPY7x2pt33nmzwtJbb711OuWUU9JTTz3VpO2NGTMm/frXv07rrrtuaVr5RRZZJK233nrZdt57771mzfpXLHb9/vvvZ/uL2O5iiy2W5pprruw1H+J/Lq5/4403pqaKGdKKt4vXa6X/r3IPPvhg9lqN10S8Z+K9s/HGG6fLLrssGy5YabfffnvacccdZ9mXRQ2q4r6sqSZNmpRt74gjjsj2A8XHJ15L8R7aeeed07Bhw9LUqVPr3UbsB5vynixe6vLMM8+ks846K2233Xal13K0I9qzwQYbpBNPPDH997//bdb/BtCqCgDQRAMGDIhvtNnl1FNPbdZtY/3ibWM79VlmmWVK67311lt1rvP4448XllxyydJ6jV2+/fbbOtvR2CXaUpcnnniisMIKKzR6+y233LLw0UcfNfrYfP7559m6DW1rww03LLz//vuFvffeu/S3a6+9ts7tzbzOyJEjC0sttVSd2/3ss89Kt7vwwgsLc801V5Mem80337zw8ccfN/h/lT/WsTxu3LjCmmuuWe82u3TpUvjHP/5Ruv1TTz3V4PM833zzFUaNGlVoTWPHji2stdZajf7//fr1y/6furT09dVUrfk6nPl9d9VVV2XPR33bXHDBBQt33HFHg9t87bXXCqusskqTH4/XX3+93m1Nnjy5cNBBBzX6+uzWrVvh4osvbvJ+7MEHHyzcfvvt2f8z87bmn3/+bP2zzz679LeBAwc2uO3y9pZvM/ZZlf6/Quzr9ttvvwa3teqqqxZeffXVWd6freHLL7/MHqOG7n+jjTZq8r4sXuM9evRo0utn2WWXLTz77LN1bie239TXYVxmts466zTpdnPPPXf2egGoBQWVAciVd955Jzvb/+WXX2a/R6+SddZZJ+sFEL13vv7666zHz/PPP58mTpw4y+3jzHicRX/33Xezs8FhiSWWyM7+zmzhhRee5W+PPPJI2mabbUpnv+Msb2wzejfEmeMnnngijRs3Lrsuek1suOGGaeTIkdlZ+LpMmTIl/ehHP8puVxTtibPscWY4thW3f+yxx0q9H5rjjTfeyHpffPHFF1lviei5E9uPejDxv5SL3gHfffddthz3s8oqq2TtjjP/0cPnxRdfTC+//HJ2/b///e/0wx/+MGt3nJlvTDwX8biNHTs26wESU4JHD4noAfKvf/0rezzjsYjnIe4nenzE9uN2PXv2zNodz0ecGY/7juvjuuixEkWLo8dWS0XviWjXRx99VPrbGmuskb7//e9nz/Nzzz2Xta14Fj/O1sdj2Ldv3xm2U+ylEa/RESNGZMvx2M88U1Zdr6+mau3XYbk77rgje82E6HkRvSXitRjPXbwOp0+fnr1+fvzjH6c777wzez/OLP73eP7i/VoclhY9weI1FduKdsd7MN6nH3/8cYPtifd03Efcd1H0ourXr19acMEF06effppdF6/f6OVx+OGHZ6+N6AnTmMcffzzrVRavp3g+4nUWr7fowRbPd4geMCeccEL2f99///3Z66Oxx/Gee+4p1VyKfVP//v2r8n/Fa6y8d1H0PopeisX3TvRUih5yAwcOTDvssENqTfEYRk/K8v1KvMfjMY3Xf+yL4jUYl3ifN2VfFo9h9N4LvXr1SquttlrWC2ieeebJXkOxzej1Fb2BYr8f799nn302e8zLxesu3peNvSfrU+yRE/u6aENsP/Y5keVGr68nn3wyex3HY3Dcccdl6x577LHNePQAWkFNIiUAcqkt9Nw56qijStdvvPHGhXfffbfeM9gPPfRQ4Wc/+1lh2rRps1wfZ+yb0p5yn3766Qw9SVZcccXC6NGjZ1nv+uuvz860F9fbfvvt693mSSedVFqvY8eOhXPPPbfw3XffzbBO9BBZd911s3XKe1M0pedOp06dsp+HHnpodla93NSpU2e4r2HDhmW9AyZMmFBve59//vnC2muvXdr+GWec0aTnvNjuAw44oDBx4sQZ1nvnnXcKK6+8cmndaH/0jOnQoUPhtNNOK0yZMmWG9V966aXCYostVlr/N7/5TaGl4j7KexX16tWr8M9//nOW9e67775Cz549S+v94Ac/yB7HusTrt7V66VT6dVj+vuvcuXP2WjzvvPNmeS2+/PLLhdVWW620bjwP0Z6ZXXDBBbP0FKnL9OnTsx5aBx98cOG///1vnesMGjSotK2+fftm792ZxXv8sssuK73OoidMfb1lyvdj8f6I11m8jmd+HqNXTdFmm21Wuk1TetDssssuje4rW/v/GjFixAy9SA477LDCN998M8M67733Xtbrrvg8t2bPndNPP720vXhMf/vb386y743eXMX3Wfn9N9Rz59e//nXhxRdfrPd+P/jgg8LPf/7z0ra22GKLeted3fdkvD7vvvvuWR7Povg/43+YZ555Sj143nzzzSZvH6A1CHcAaLLyL0XRTT0Cg6Zeyru1tyTciS/9TRnG0ZjZCXdOOeWU0m1iyEV9X0bDrbfeOsMXrYcffniWdeJLcdeuXUvrnHXWWfVuL4ZPlT82TQ134jJ48OBCa4phZMVwZfHFF68zPKtrCNxee+1V7zZj6NjMQxwa+sIZwUVxvRj601LXXHPNDEMr6hveESKMKIZmcbnuuuuqGu609uswzPza+v3vf1/vNmNITXnAdfLJJ8+yzq677lq6vq6QrKkeeeSR0nZiCFpjw8vKh9/86Ec/anQ/FpczzzyzWa+P9ddfv9H3R3kIW9d+qrX/rwjhyode7rPPPvVuKwKK733ve01+rzVF/M/du3cvbS9C2fp8+OGH2X6jKfuy5thmm21K23vllVeq+p4suummm0rbP/bYY1t9+wANEe4A0GQzfyma3UtLwp3opVC8Pr5QVCvciR4G5b1F/vCHPzTry8Yee+wxy/WXXHLJDF80ymsD1SVChOaGOxEe1dWzoqXiTHbxPl544YVGw504S/+///2vwW0uvfTSpfUXXXTRWXrszPwFtXjmP3oJzNwbqLnWW2+90n0fccQRzfr/6/uyX4kvkpV4Hc78vltuueUafS1edNFFpfWXWGKJrF3lymtI/ec//ynMrp122qm0naiN0xTFXmDxuqirLlT5fiza3tj/Gr744osZekG98cYb9a77xz/+sdHXRmv/X/fcc09pe9HOxuph3X///a0a7kTvouK2evfu3eB7N0RNp9YOd/7yl7+Uthevz1qEOxF0F2sERa8+gGpScweAXFlqqaXS66+/ni3HbFnF+gaVFvVYijMExUw6TanVMHjw4HTvvfdmy8VZecqV/2333XfPZhFqSNQ4Oeigg7IZp5pqq622yup3NFfUHIm6LfF/R92LqA9SPkPU6NGjS8v/+c9/sto0DYkaQjGrTENWX331Um2LmDUsZqKpT8wKFfVJon3Rrqi30Vgb6hN1OMr/n/32269Jz+3ll1+eLT/99NPZ4xN1QPL4OpxZ1Jhp7LW41157paOPPjqr0RQ1YaLuUcymVv4+LYr3afGxao6oo1Kc7SvqNMUsRU0RNWZipq94XUTNmoZqy8R7qrH/tXj/8Zr861//mv3+5z//OZvFqi5xXfnjVI3/K2bHKop6Oo3Vc4p6SFFPKeoetYby+499WUPv3bDHHnukww47rMEZrmYWNXZinxR1r6LuUbxvizXCQvn/EvukSnnhhReyekyxz4kaSFErrFxxtq1oZ9RpinpTANUg3AFgtpx66qml6a2bItb9zW9+0+L73W233bKCuuH444/PviT97Gc/y6bZjUKblVIsrhpimummFMONIrZF8YU8vgRHMeO6voDElMeNiYLREYCUBxGNieKszRHFViMwizCg/ItTQxoriBui3Y0pD6GiaGljYpr2orqKZzfny1rxf41ivzH1dWOiyHKEORHqxG2jMHAUWK60SrwOZ1ZX8d+6nqu4/+L09dGu8nAn3qfXXHNNKdyJAtR77713VkB45mK3DT0v8fgWC6cfeeSRTbpdhG1FxYLOrfH+iKCmsXBnwoQJ6eGHHy61OYKOavxf5a+Lpjx/EUDEPufWW29NraG59x/FjGOfEMWPGxOFpeOxjkLIxUL6rbFPaq7rrrsu/e53v8sKizdFFFeOQvazE64DzA7hDgC5Er0Q/vGPf5RmuoqZluISll566ayHSJzh3nHHHbNZb1pL+QxKyyyzTJNuEz1VYqapYk+b+MJR/qW6fJvlPR0aEgFWc8KdpsyOVHTfffdlj9vMZ6Ib05QvXE2Zzaq8B0Vz148vUrNr5ueheOa9IXE2PtaNnhSV+jJZrdfhzOJ91BSxXjHcKW9XiBAnZne6+OKLS8FEMZyI9sQMXJtuumk221l9oWyEUEWffPJJuvTSS1NzFWesao33R8xqF/uUePziC378PzFTX7kbbrih1MOtuH41/q/yx785z19rmd37byzcefvtt7PZtoo9+pqqqSFQU8Tzuf/++6drr7222beNdgh3gGrRTxCAXImhKHG2+eqrr86mfS4XXwDijHoEQPHlNX7GWd/WUJyONzRn+E35ujN/4SjfZvTKaYroWdIcMXypqV/OopdBMdiJ4OCss87Kpi2OL6MxJCKGGPy/en1Zz62i+HtjmhKYtGT9tvbc5rmtTX0tNrbNiy66KHuvxhTt5T744IP0t7/9LQt/4gt+DI2q68t79HpoqRgC1Rrvj7p64lx//fWzrFP+t5///Od1bqcS/9fs7Etacxhhpe4/hggWXxvR2yeGAka4/+abb2b3Gb3mivuk8qFhTdknNdUf//jHGYKdCO2iF08Mu4qQLfaZxTbEpTx0bc12ADRGzx0Acie++MeZ1LjEGfQYBhE1KB599NHsQ3+xJ8ewYcOyGiOjRo1q1hn6xkKV4pCKpihfN76czLzN4he9CE+au73WFF9gim1Zc8010yOPPJLVA6lPtcKMaqjEc5vnts7Oa7G+be68887ZJb6gx3vx8ccfz96nxR4/8WU4gp7idX379q3zy38MlYuhb7UWQ7OKPW3+8pe/pD/84Q9Z4Bziy35cij3PokZPXSrxf5W/LmqxL6nEvixeD3Epbj/q7cwc6Fdjn3TuueeWlmNocX21lirdDoDG6LkDQK7Fl8EDDjggDR8+PI0bNy4r7PrLX/6y9IUr/tYatX7Kw6GmDhGIosTlxY9nHqJR/nvU6miKpq7XXMWhbeGkk05qMNgpDpdoL8qf23h8ywtH1yfOyJfXPWnNIYDVfh3OrKnbbc7/Hz10ovhz1N95+eWXs/uI92Wxl0cMT4r3bbnyAtzFItK1tv7665dqBkUPpGJh5Jl77URvpBgKV5dK/F+z87porB5Rre+/fJ8U9ZoaCnYqtU+KNhYL+C+wwALphBNOaHD9qP3V2FBAgEoR7gDQ7sKe8847b4ZA5+9//3uLh/2stdZapeWos9KU4V7Rm6hoscUWm6XOSRTlLXryyScb3d6kSZPSSy+9lCqhvA5IY7NOxVCI8v8t76L3RDEMjLPuxd4XDYneFsWeB3Hb6O1UDZV4Hc4sekg05vPPPy/VGwo/+MEPUnNEvaLoAXHVVVeV/nb//ffPUO8p3h9dunQpBVRvvPFGaguigPvMM2NFIHjjjTc2OEtWJf+v8tdFU56/aG9T9jmVuv8YUtXYvqw5+6QQvQ0b09z9fnkbomB4DM1rSAxjbUo4DFAJwh0A2qXyaYLjDPvMys+qN6UY7yqrrJJ9MS6GG3XV25hZDAsriiLPM4uCskUxC09j9UFi+EoEPJVQPl1vY8Mqoph1W+lJ0RpiSNHaa69d+j16gTXnuY2aMtWYBr1Sr8OZRUjR2ExpEWoU11l88cWzmbNa+j6N92F5WBX1cDbffPPS75dddllqC8qDm3gvxPslhoYWe6JEcDVgwIB6b1+J/6v8eb3nnnsaDf1ixsHW7AVYfv8xXK2xfWqs01jh9ubskyKEueOOOxptZ3P3+81pQ7j88ssbXQegUoQ7AORKU2clKu/y36tXr1muL59C+t13323SGd8DDzyw9Pvpp5/e4O2it9Ddd99d+v0Xv/hFncVCi1823nrrrXT++efXu72oZ3HyySenSll++eUb7OlUXng5ipq2NwcddFBpOWqqxHTV9Ylpva+88soGn9tKqcTrcGYxlLGh12KEpXG/RVH7auYeEbPzPo0v0jNP7X7ccceVlmPmrQceeCA1VaUCyBiWFcOzij1QIuAp9uAp9uxprIdIa/9fW221VWnGvQghjj322HpvH0P0jjnmmNSaYl9WHGIXz+nZZ59d77oxBK+xujXN2SdFyBjvialTpza6zRhaVQxsYl/WWMCz3HLLlZ7L6GlUrOlWX2B11113NdoGgEoR7gCQK1G7I76Ix5ny+mYiianCYyaeom222abOD+3FLyNRq6E4TXNDjjrqqLTkkkuWvqBsscUW6T//+c8s6910001pzz33LP0ehVVjOt+ZLbTQQjPUGTn++OPTBRdcMMv/NX78+GyGlvhZHM7R2sqLv8YsWXX1CIlpi6NHQnx5q1ZPlWqJL+TFoVXxJTGm8i6ffacovoTH66nYyyqGI5U/19XQ2q/DmXXu3DkLHy688MJZXotjxoxJW265ZTacqFg/pq6wr3///tkX/nvvvbfeL91RDD1qqRTF/xH3XS5eb8V14jHfdttts9dn+exMMwcXEbbsuOOOM/QKqmTvnegZdcstt9R5XX1a+/+KoYFnnHHGDG2K10l5raViMBSvgxhWOPNj3RJRQLo8UIrwJgKemXuARf2aeP1ET5vG7j8ek2KwEgW3f/WrX83SczH+n1133TULMJuyT4r954orrpgtR7ATj2lDopZUMciL90LUUoq6buXi7xEIx+xo8TzUV2sJoOIKANBEAwYMiGIC2eXUU09t1m1j/eJtYzv1WWaZZUrrvfXWW7NcX7wuLvPOO29h4403Luy1116Fgw46qLDrrrsWVltttRnWWWSRRQrvvvtunff105/+tLRe9+7dC7vsskvhqKOOKhxzzDHZ5be//e0st3n44YezdYu369ChQ2H99dcv7Lffflk7+vTpM8P9r7jiioUPP/yw3v938uTJhXXXXXeG2yy55JKFPfbYozB48ODC5ptvXujUqVP29/79+xd+9rOflda77rrr6tzm3nvvXVrn2muvrfe+Z25H3759Z2jHKquskrVjn332Kay99tqlv6+55pqFY489ttHXQvlz3pTXS3PbXf56fPDBBwst9corr2Svl/LHIP7XaFdcYrn8ul69ehVee+21ercXr9/iuvG6bk2t/Tosf99dcMEFpeXevXuXXoubbLJJoWPHjqXr4nV59913N7q9bt26FdZbb73Cnnvumb1Pd9tttxleT8V1nn/++Xpfm1tttdUM68f/vtlmm2XPy4EHHphts1+/foUuXbqU1onfK/W6+eijjwpzzz33DG2Ky1prrdXkbbT2/xVi/fLtLbjggtl+7YADDij86Ec/Km1nueWWy/Z1s7s/r8uUKVMKG2644Qz3v/jii5deP5tuumlhrrnmyv4er4fy/W997/dBgwbNsr0ddtgh294Pf/jDQufOnUvHgiuuuKJJx5hf//rXpfXiORw4cGDhiCOOKO3341LugQcemOF1H7eJ/yXea/F4R5uK18Uxo7FjGEClCHcAyFW406NHj1m+UNV3iS/jY8aMqfe+xo8fX1hsscXqvX19X8hHjRpVWH755Ru9//jy0dAX6qLPPvssC3Ea2tYGG2xQeP/992f4QnTbbbe1WrgTIqho7P+KL28TJkxoUnCTt3Cn+BjEF/TGntsf/OAHhTfeeKPBbVUy3Gnt1+HM77vLL7+89MW5rssCCyxQuPXWW+vd3uqrr97k92kEDY899liD7Zs2bVrh5JNPniHQaugSX8APPfTQir5utttuu1nu97zzzmvWNlrz/wpTp06d4X1U12XllVfO9ovNfX82xRdffJGFSI3ty957770mvd+//vrrWQKwmS8RQI4cOTJ7LptyjPn888+zx6Chbc4s3g/FkL2uS4Q/p5xySmH69OnCHaBmOtW65xAANEcMQ4lZUWJYVgylim7+UQMkhh/EMKvevXunfv36ZV31Y/hCeUHMmS2zzDLZ8IRLLrkkm6knhonEbEmNFTaObvoxPCWGLkW3/hgSE8NUYiaVKHa70UYbZcNhog5GU0QdiJj2N4oqjxgxIqvpEgVRY0hAFNCN7v4xxCW2X14oNW7X2jONPffcc9kQg1tvvTUbfhBDauJ/itlqog277bZbaWap9igegxjWF8NsooD1U089VRqCFLWb1ltvvWxoRry+mjvzTmtr7ddhuajNs/HGG2dTl8dQtGLx3WWXXTYb1hPDHqOQcn2iLTFrUgxti8cwXksxFCfqwcT7NNoXs0bFezReU40NN4zXXNT5ifuN90i06ZVXXslq+8Twmvnmmy97P8frNIr7Dhw4cIbpuSsh3pflNVaijc0dotfa/1c891EQPKadj5nIYqa0eE0suOCCWa2geKz322+/1KNHj1QJ0d4Yihf7j2hH7KPL92Ux/DGGrTU261RRvFZiezfccEO67rrrsv1TTDce24uaPPE+3GeffbL/L4ZuNXUIWbQrClnHcK54D8Xsbw3V34n3w4YbbpjVoorXdLyWozB2DI+M4tjxmJbPGAZQCx0i4anJPQMAzRZfJorT80a9iah5AgDAnE1BZQDIiZEjR5aCnZgZR7ADAEAQ7gBADsTwqPJZiWKIFAAABOEOANTYwQcfnK655pqs3k9dXnrppayuQ9SCCVEv45BDDqlyKwEAaKvU3AGAGtt0002zAtFRVDaKzK644opZgBOFQ1944YX08ssvx/Qt2bpRxHfYsGFp3333rXWzAQBoI8yWBQBtxJQpU9KTTz6ZXeoSs2PFTFaGZAEAUE7PHQCosZj16rbbbst678SU0TENckz5HhZeeOG0+uqrpy233DKbbre1pz8HACD/hDsAAAAAOaagMgAAAECOCXcAAAAAcky4AwAAAJBjwh0AAACAHBPuAAAAAOSYcAcAAAAgxzrVugG0DZMnT04vvvhitrzIIoukTp28NAAAAKC1TZs2LX300UfZ8hprrJG6du3a4m36Bk8mgp1111231s0AAACAOcZTTz2V1llnnRZvx7AsAAAAgBzTc4fSUKzy5HDxxRevaXsAAACgPXr//fdLI2fKv4u3hHCHTHmNnQh2evfuXdP2AAAAQHvXqZXq3RqWBQAAAJBjwh0AAACAHBPuAAAAAOSYcAcAAAAgx4Q7AAAAADkm3AEAAADIMeEOAAAAQI4JdwAAAAByTLgDAAAAkGPCHQAAAIAcE+4AAAAA5JhwBwAAACDHhDsAAAAAOSbcAQAAAMgx4Q4AAABAjgl3AAAAAHJMuAMAAACQY8IdAAAAgBwT7gAAAADkmHAHAAAAIMeEOwAAAAA5JtwBAAAAyDHhDgAAAECOCXcAAAAAcky4AwAAAJBjwh0AAACAHBPuAAAAAOSYcAcAAAAgx4Q7AAAAADkm3AEAAADIMeEOAAAAQI4JdwAAAAByrFOtG0Db1m/IiJre/zNDB9X0/gEAAKCt03MHAAAAIMeEOwAAAAA5JtwBAAAAyDHhDgAAAECOCXcAAAAAcky4AwAAAJBjwh0AAACAHBPuAAAAAOSYcAcAAAAgx4Q7AAAAADkm3AEAAADIMeEOAAAAQI4JdwAAAAByTLgDAAAAkGPCHQAAAIAcE+4AAAAA5JhwBwAAACDHhDsAAAAAOSbcAQAAAMgx4Q4AAABAjgl3AAAAAHJMuAMAAACQY8IdAAAAgBwT7gAAAADkmHAHAAAAIMeEOwAAAAA5JtwBAAAAyDHhDgAAAECOCXcAAAAAcky4AwAAAJBjwh0AAACAHBPuAAAAAOSYcAcAAAAgx4Q7AAAAADkm3AEAAADIMeEOAAAAQI4JdwAAAAByTLgDAAAAkGPCHQAAAIAcE+4AAAAA5JhwBwAAACDHhDsAAAAAOSbcAQAAAMgx4Q4AAABAjgl3AAAAAHJMuAMAAACQY3NkuDNx4sR00003pWOOOSYNGDAg9enTJ80///ypc+fOqVevXmnTTTdN55xzTvrkk0+atL3HH3887bXXXmmZZZZJXbt2TYsttljaeuut04033tisdsX6W221VXb72E5sL7Y7atSo2fxPAQAAgPauQ6FQKKQ5zAMPPJC23HLLRtfr2bNnuv7667Ogpj6nnXZaOuOMM9L06dPrvH7bbbdNt9xySxbW1GfSpEnpxz/+cbrnnnvqvL5jx47plFNOSaeeemqqlAkTJqSllloqW37nnXdS7969s+V+Q0akWnpm6KCa3j8AAABU4/t3S8yRPXdCPJCDBg1KF154Ybr11luz3jGPPfZY+stf/pJ+8pOfpLnmmit9/PHHaYcddkjPP/98ndu48sor029+85ss2FlhhRXSsGHD0lNPPZVuv/32tNlmm2Xr3H333Wm//fZrsC1xfTHYidvF7WM7sb3Ybmw/QqSrrrqqAo8EAAAAkGdzZM+d7777LgtvGhIBy84775wtx88IgMp9+umnafnll09ffPFFWnrppdMzzzyT9fQpv4+43Z133pn9/uCDD2bDvWb273//O22xxRbZ8vbbb59uu+22GdoWAVO/fv3Sf//737TAAgukN998My244IKptem5AwAAAJWn504raSzYCTvttFNaaaWVsuVHH310luuvvvrqLNgJZ5999gzBTvE+LrvsstJ9DR06tM77Offcc7OfnTp1mmH9othubD98/vnn2f0CAAAAzNHhTlPNO++82c/JkyfX2bMnzDfffGmXXXap8/aRvv3whz/Mlv/1r3+lL7/8cobr4/f4e4j16kvrYvtxPyF69gAAAAAUCXfq8dprr6X//Oc/2fLKK688w3VTp07NauKE/v37Z7Ns1Sdm4wpTpkxJo0ePnuG6p59+OttW+Xp1ie2vv/76pdt8++23s/1/AQAAAO2LcKfMN998k15//fX0hz/8IQtbpk2blv39qKOOmmG9sWPHZjV16gp+ZlZ+/ZgxY2a47pVXXqlzvYa2E22KNgIAAACETnP6wzB8+PC077771nv98ccfn37605/OUvyoqLHCR8UiScVCSa21nVVXXTU1R/l91eX9999v1vYAAACAtmGOD3fq8/3vfz+benydddaZ5bry2jk9evRocDvzzDNPafmrr76qyHaaojwcAgAAANqPOX5YVsyK9eKLL2aXqKNz4403ZlOYR72dPffcM911112z3Ka8wHJD9XZCly5dSsuTJk2qyHYAAACAOdcc33NngQUWyC5F0VNnjz32SH/605/S3nvvnXbcccc0bNiwtM8++5TW6dq1a2m5WBC5PlFIuahbt24zXNda22mKmYeE1TUsa9111232dgEAAIDamuPDnfr8/Oc/z3rt/PWvf02HHXZY2mGHHdJCCy00wxTpTRki9fXXX9c79Kq1ttMUjdX0AQAAAPJpjh+W1ZDotVMMVv7xj3/UGZQ0Vqi4vMfMzHVvWms7AAAAwJxLuNOARRZZpLT89ttvl5b79u2b5pprrmz51VdfbXAb5devssoqM1xXPuNVU7fTqVOntOKKKzb5fwAAAADaN+FOA9599906h0JF8eNifZpRo0Y1WC/n4YcfLhVEXnvttWe4Lur7FAspF9erS2z/iSeeKN1m7rnnnu3/CQAAAGhfhDsNuPnmm0vLa6yxxiyzbIWJEyemW2+9tc7bx1CrBx54IFveYostZqixE+L3+HuI9eobmhXbj/sJMZMXAAAAwBwd7gwfPnyGacjrcv7556d77rknW15uueXSxhtvPMP1gwcPTvPPP3+2fPzxx6dPPvlkhuu/++67dMghh2Q/w5AhQ+q8n1/96lfZz2nTpqVDDz20tH7Rxx9/nI477rhsOWb1ivsFAAAAmKPDndNOOy0tueSS6cADD0wjRoxIjz32WHr++efTyJEj0+WXX5422mij9Mtf/jJbN4ZNXXXVVaUaO0Uxc9bZZ59dqsez3nrrpWuvvTaNHj06/f3vf09bbrlluvPOO7Pr99xzz7TpppvW2ZbNN988m3o9FG8XP2M7sb31118//fe//82uj/tbcMEFK/rYAAAAAPnSoVAoFNIcZtlll52hQHJ9Yjara665Jgtc6nPqqaemM844I9X3MA4cODD97W9/S127dq13G5MmTUo//vGPSz2FZtaxY8d08sknZ6FUpcSQsOIsXDEzV3Emr35DRqRaembooJrePwAAAFTj+3dLdEpzoPvuuy/dfffdWY+dN954I33wwQfZsKpu3bqlXr16pe9///tpu+22S7vttlvq3r17g9v6zW9+k7beeut06aWXpkcffTTbVgyfWnPNNdO+++6b9dppTNxvtOeGG27IhoxFL6LPP/88LbrootlwsMMOOyz179+/FR8BAAAAoL2YI3vuMCs9dwAAACCfPXfmyJo7AAAAAO2FcAcAAAAgx4Q7AAAAADkm3AEAAADIMeEOAAAAQI4JdwAAAAByTLgDAAAAkGPCHQAAAIAcE+4AAAAA5JhwBwAAACDHhDsAAAAAOSbcAQAAAMgx4Q4AAABAjgl3AAAAAHJMuAMAAACQY8IdAAAAgBwT7gAAAADkmHAHAAAAIMeEOwAAAAA5JtwBAAAAyDHhDgAAAECOCXcAAAAAcky4AwAAAJBjwh0AAACAHBPuAAAAAOSYcAcAAAAgx4Q7AAAAADkm3AEAAADIMeEOAAAAQI4JdwAAAAByTLgDAAAAkGPCHQAAAIAcE+4AAAAA5JhwBwAAACDHhDsAAAAAOSbcAQAAAMgx4Q4AAABAjgl3AAAAAHJMuAMAAACQY8IdAAAAgBwT7gAAAADkmHAHAAAAIMeEOwAAAAA5JtwBAAAAyDHhDgAAAECOCXcAAAAAcky4AwAAAJBjwh0AAACAHBPuAAAAAOSYcAcAAAAgx4Q7AAAAADkm3AEAAADIMeEOAAAAQI4JdwAAAAByTLgDAAAAkGPCHQAAAIAcE+4AAAAA5JhwBwAAACDHhDsAAAAAOSbcAQAAAMgx4Q4AAABAjgl3AAAAAHJMuAMAAACQY8IdAAAAgBwT7gAAAADkmHAHAAAAIMeEOwAAAAA5JtwBAAAAyDHhDgAAAECOCXcAAAAAcky4AwAAAJBjwh0AAACAHBPuAAAAAOSYcAcAAAAgx4Q7AAAAADkm3AEAAADIMeEOAAAAQI4JdwAAAAByTLgDAAAAkGPCHQAAAIAcE+4AAAAA5JhwBwAAACDHhDsAAAAAOSbcAQAAAMgx4Q4AAABAjgl3AAAAAHJsjg13Ro8enU4//fS01VZbpd69e6cuXbqkHj16pL59+6Z99903jRw5stFtDB8+PHXo0KFJl1i3Md98800655xz0jrrrJMWWmihNM8886SVV145HXPMMentt99upf8cAAAAaE86pTnQJptskh599NFZ/j516tT0+uuvZ5cIYwYNGpT++Mc/ps6dO1e8TW+88UYaOHBgdt/lXnvttexy9dVXpz//+c9pu+22q3hbAAAAgPyYI8Od9957L/u5xBJLpJ/85Cdp4403TksvvXT67rvv0qhRo9J5552X3n333TRixIj07bffphtuuKHRbd53333Z9uoTvYPq8+WXX6Ztt922FOwccMABaY899kjdunVLDz74YDrrrLPSxIkT0+67754ee+yx9P3vf3+2/m8AAACg/Zkjw50Y6vS73/0u7brrrmmuueaa4br1118//fznP08bbrhhGjt2bLrxxhvTL37xi6y3T0NiONeyyy47W+0ZOnRodl8hhmUNGTKkdF3//v3TpptumgYMGJAN2zrqqKPSQw89NFv3AwAAALQ/c2TNnbvuuivttttuswQ7RT179sx67xTdcsstFWtL9Ay66KKLsuVVVlklq68zsw022CDtv//+2fLDDz+cnn766Yq1BwAAAMiXOTLcaYrNNtustDxu3LiK3U8Mu/riiy+y5b333jt17Fj3U7LPPvuUlm+77baKtQcAAADIF+FOPaZMmVJarq+HT2son5Urhl7VZ+21107du3fPlqPuDgAAAEAQ7tQjhj8VxXCpxsT06VFQOWbWimFdUbvnpJNOygozN+SVV16ZoRZQfTp16pT69OmTLY8ZM6aJ/wUAAADQ3s2RBZUbM3369PT73/++9HvU52lMeZHjTz75JLs8+eSTWe2eCy64IB100EF13m7ChAnZz3nmmSctsMACDd7HUkstlV544YX00UcfZT2LunTp0uT/qXg/9Xn//febvC0AAACg7RDu1OH8889PTz31VLa8yy67pH79+tW77vLLL5+tE7NaRfgS3nzzzfS3v/0tK8Q8efLkbLatDh06pAMPPLDOadBDjx49Gm1XBEBFX331VbPCnWLbAAAAgPZFuFPHcKzjjz8+W+7Vq1e6/PLL61135513zoogR3BTbp111km77757NitXBD8xI9bRRx+ddthhh7TYYovNsG6EPyGGczWmPMyZNGlSs/83AAAAoP1Rc6fMyy+/nAU206ZNS127dk0333xzFvDUZ/75558l2Cm33XbbpVNOOSVb/uabb9KwYcNmWSfuJ0ydOrVZRZ67deuWmuOdd95p8FLsqQQAAADki3Dn/3nrrbfSVlttlT777LNsdqybbropbbLJJi3ebgzFKgZA5UWai+add97SMKvGfP3116XlpgzjKte7d+8GL4svvniztgcAAAC0DcKdlNJ7772XfvjDH2Y/I4i55ppr0o477tgq246ePwsvvHC2XNfMWRGsFIObzz//vMFtRQ+bsMgiizSr3g4AAADQfs3x4c7HH3+cttxyy6wIcrj44ovToEGDWvU+Ghq6teqqq5aWX3311XrXi6Fi48aNa/LU7AAAAMCcYY4Od7744ou09dZbp1deeSX7PaY/P/TQQ1v1PmLa8giQwhJLLDHL9RtttFFpua5hW0WjR48uDcvacMMNW7WNAAAAQH7NseFOFDjedttt07PPPpv9fuKJJ6bjjjuu1e/nqquuSoVCIVseMGDALNdvuummWWHmcN1115XWndnw4cNLy1H0GQAAAGCODXdiZqoISB577LHs9yOPPDKdeeaZzdrG+PHj03PPPdfgOjEV+umnn16a3WrfffedZZ2YAv2II47IlseMGZPOPffcWdYZNWpUaaatCIhiqnUAAACA0GlOfBj23HPPdP/992fLm2++edp///3TSy+9VO/6EcD07dt3lnBns802S/3790/bb799WnPNNUvTpkf9nltuuSW7FHviRGiz5JJL1rn9IUOGpL/85S9p7Nix6dhjj01vvPFG2mOPPbJA6MEHH0y/+93vspo78fsFF1zQio8EAAAAkHcdCvWNA2rHGipwXJdlllkmC3PKPfTQQ1m405ju3bun888/P5sSvSER6AwcODC9/vrrdV4/33zzpT//+c9pu+22S5UwYcKEtNRSS5Vm5SrO4tVvyIhUS88Mbd3i1gAAAFBL9X3/bok5sudOa+jXr1+6/vrrsyFTUez4/fffzwonRw+bBRdcMK222mppiy22SIMHDy716GlInz59smFel156abr55puzsCeGj8UTHqFPDB2LkAkAAAAgzek9d5iVnjsAAACQz547c2RBZQAAAID2QrgDAAAAkGPCHQAAAIAcE+4AAAAA5JhwBwAAACDHhDsAAAAAOSbcAQAAAMgx4Q4AAABAjgl3AAAAAHJMuAMAAACQY8IdAAAAgBwT7gAAAADkmHAHAAAAIMeEOwAAAAA5JtwBAAAAyDHhDgAAAECOCXcAAAAAcky4AwAAAJBjwh0AAACAHBPuAAAAAOSYcAcAAAAgx4Q7AAAAADkm3AEAAADIMeEOAAAAQI4JdwAAAAByTLgDAAAAkGPCHQAAAIAcE+4AAAAA5JhwBwAAACDHhDsAAAAAOSbcAQAAAMgx4Q4AAABAjgl3AAAAAHJMuAMAAACQY8IdAAAAgBwT7gAAAADkmHAHAAAAIMeEOwAAAAA5JtwBAAAAyDHhDgAAAECOCXcAAAAAcky4AwAAAJBjwh0AAACAHBPuAAAAAOSYcAcAAAAgx4Q7AAAAADkm3AEAAADIMeEOAAAAQI4JdwAAAAByTLgDAAAAkGPCHQAAAIAcE+4AAAAA5JhwBwAAACDHhDsAAAAAOSbcAQAAAMgx4Q4AAABAjgl3AAAAAHJMuAMAAACQY8IdAAAAgBwT7gAAAADkmHAHAAAAIMeEOwAAAAA5JtwBAAAAyDHhDgAAAECOCXcAAAAAcky4AwAAAJBjwh0AAACAHBPuAAAAAOSYcAcAAAAgx4Q7AAAAADkm3AEAAADIMeEOAAAAQI4JdwAAAAByTLgDAAAAkGPCHQAAAIAcE+4AAAAA5JhwBwAAACDHhDsAAAAAOSbcAQAAAMgx4Q4AAABAjgl3AAAAAHJMuAMAAACQY8IdAAAAgBwT7gAAAADkmHAHAAAAIMeEOwAAAAA5NseGO6NHj06nn3562mqrrVLv3r1Tly5dUo8ePVLfvn3Tvvvum0aOHNms7d17771p5513Lm0rfsbv8femmjZtWrriiivSxhtvnBZZZJHUrVu3tMIKK6SDDjoovfzyy7PxXwIAAADtXYdCoVBIc5hNNtkkPfroo42uN2jQoPTHP/4xde7cud51pk+fng488MA0bNiwetcZPHhwuvLKK1PHjvVnaR9//HEaOHBgevrpp+u8PgKjSy65JNtWJUyYMCEttdRS2fI777yThVOh35ARqZaeGTqopvcPAAAA1fj+3RJzZM+d9957L/u5xBJLpCOPPDLdcsst6amnnkqjRo1Kf/jDH9KSSy6ZXT9ixIi0zz77NLitE088sRTsrLXWWunGG2/MthU/4/dw9dVXp5NOOqnebXz33XdZL59isLPLLrtkPX6efPLJdNFFF6VevXqlKVOmZD14mtMTCAAAAGj/5sieO9ttt13WK2fXXXdNc801V529aDbccMM0duzY7PeHH3446+0zs7h+tdVWy4ZTrb322umRRx7JhlIVffPNN2nAgAHZELBOnTqlMWPGpD59+syynWuuuSbtv//+2fIhhxySLr300hmuf+ONN1K/fv3SxIkTs9vHdmJ7rUnPHQAAAKg8PXdayV133ZV22223OoOd0LNnz3TeeeeVfo+ePXW54IILsmAnXHzxxTMEO6F79+7Z30Osd/7559e5nXPPPTf7udBCC6WhQ4fOcn0EOieccEIp6Lntttua+J8CAAAA7d0cGe40xWabbVZaHjdu3CzXR4enO+64I1teeeWV0/rrr1/nduLvK620UrYc68/cUSp6/0RPnBCBUwRCdSkfHibcAQAAAIqEO/WIGjdFdfXweeutt0q1e2LoVUOK17/77rtp/PjxM1xXPitXQ9tZbLHFspm8wmOPPdbk/wMAAABo34Q79Yg6O0WrrLLKLNe/8sorpeXoudOQ8uuLvXRasp0Yk/f11183uC4AAAAwZ2jdqrztRExv/vvf/770ewyXqqsAUlFjxY+KhZKKwUxLtxNDu+J2xeFeTVF+P3V5//33m7wtAAAAoO0Q7tQhCh/HdObFacljpqqZffnll6XlHj16NLi9eeaZp7T81VdfVWQ7jSkPmAAAAID2w7CsOoZjHX/88dlyr1690uWXX17nepMnTy4td+7cucFtdunSpbQ8adKkimwHAAAAmDPpuVPm5ZdfTjvvvHM2bXnXrl3TzTffnAU8dYnri6ZOndrk4swzT5c+83bKf2/Odhoz83CwuoZlrbvuus3aJgAAAFB7wp2y2a+22mqr9Nlnn2WzY910001pk002qXf9eeedt8lDpMqLH8889Grm7TQU7jS0ncY0Vs8HAAAAyCfDslLKpjT/4Q9/mP3s0KFDuuaaa9KOO+7Y5LCksWLF5b1mZq59MzvbiTYKawAAAIAwx4c7H3/8cdpyyy3Tm2++mf1+8cUXp0GDBjV6u1VXXbW0/Oqrrza4bvn1M0+rPjvbiYCovLgyAAAAMOeao8OdL774Im299dbplVdeyX6P6c8PPfTQJt12ueWWS0sssUSpCHNDHnnkkeznkksumZZddtkZrttoo41Kyw1t53//+18aO3Zstrzhhhs2qY0AAABA+zfHhjvffPNN2nbbbdOzzz6b/X7iiSem4447rsm3j6FRxaFb0aPmiSeeqHO9+Huxx02sH7cr17dv31Jvnr/+9a9Zu+oyfPjw0nIUfQYAAACYY8OdmJUqApLHHnss+/3II49MZ555ZrO3c9RRR2XFl8Phhx8+y/Tk8Xv8PXTq1Clbvy6/+tWvsp+ffvppOvbYY2e5fty4cemss87Klvv06SPcAQAAAObs2bL23HPPdP/992fLm2++edp///3TSy+9VO/6nTt3znrYzCz+NmTIkGw41+jRo7PhUtH7Z4UVVsgCmbPPPjs999xz2bqx3oorrljn9vfee++siHOETZdeemk2BOuAAw5ICy64YHrqqafSGWeckSZOnJg6duyYLrrooiwoAgAAAAgdCoVCYU57KGYeGtWYZZZZJo0fP77O66ZPn54FMRHO1CfCo6uuuioLZxoq7Dxw4MD09NNP13l9ly5d0iWXXJIGDx6cKiFm6irO5BWzchVn4+o3ZESqpWeGNl7cGgAAAPKivu/fLTFHDstqTRHYDBs2LN19991ZTZ0oshw9feJn/H7PPfekq6++usFgJ/Ts2TM9/vjj6bLLLsuKLC+88MKpa9euafnll8/Co2eeeaZiwQ4AAACQX3Pk+J5KdFaKXjdxaYkYbnXwwQdnFwAAAICm0HMHAAAAIMeEOwAAAAA5JtwBAAAAyDHhDgAAAECOCXcAAAAAcky4AwAAAJBjwh0AAACAHBPuAAAAAOSYcAcAAAAgx4Q7AAAAADkm3AEAAADIMeEOAAAAQI5VNdzZfPPN0xZbbJHefvvtJt/mvffeK90OAAAAgBl1SlX00EMPpQ4dOqSvv/66ybeZNGlS6XYAAAAAzMiwLAAAAIAca/PhTrGXT9euXWvdFAAAAIA2p82HO/fee2/2s3fv3rVuCgAAAMCcVXNnv/32q/PvJ510UlpggQUavO2UKVPSuHHj0tNPP53V2xkwYECFWgkAAACQXxUNd4YPHz5LIeRCoZDuuOOOJt0+1g0LLbRQOuGEEyrSRgAAAIA8q2i4s/TSS88Q7sQU6PH74osvnuaee+56bxfrRI2dWG+DDTZIBx98cFpiiSUq2VQAAACAXKpouDN+/PgZfu/Y8f9K/Nx///1p1VVXreRdAwAAAMwRKhruzGyTTTbJeuXMM8881bxbAAAAgHarquHOQw89VM27AwAAAGj32vxU6AAAAAC0kZ47dZk4cWL68ssv03fffdekAs0AAAAA1Djc+ec//5kuu+yyNHLkyPTpp5826TZRq2fatGkVbxsAAABAnlQ93DniiCPSpZdemi0XCoVq3z0AAABAu1LVcOeGG25Il1xySbbctWvXtNNOO6V+/fqlhRZaqDRNOgAAAABtNNy58sors59LLbVU+ve//51WWGGFat49AAAAQLtT1e4yL7zwQlY759RTTxXsAAAAAOQt3Pn222+zn2uttVY17xYAAACg3apquLPssstmP7/66qtq3i0AAABAu1XVcGeXXXbJfv7rX/+q5t0CAAAAtFtVDXeOOeaYtPTSS6cLLrggvfrqq9W8awAAAIB2qarhzvzzz5/uu+++tOiii6YNNtggXXbZZemzzz6rZhMAAAAA2pWqToW+/PLLZz+/+eab9Pnnn6fDDz88HXHEEalnz56pe/fuDd42ZtkaN25clVoKAAAAkA9VDXfGjx8/w++FQiG7fPjhh43eNsIdAAAAAGoY7uy9997VvDsAAACAdq+q4c61115bzbsDAAAAaPeqWlAZAAAAgNYl3AEAAADIMeEOAAAAQI5VtebOiBEjWnT7QYMGtVpbAAAAANqDqoY7++yzz2xPaR63E+4AAAAA1DDcCYVCodp3CQAAANBuVTXceeuttxpd5+uvv05jx45NN9xwQ7rlllvShhtumK666qrUvXv3qrQRAAAAIE+qGu4ss8wyTVpv1VVXTTvttFP661//mn7605+mww8/PP3zn/+sePsAAAAA8qZNz5a12267pb333js9+OCD6corr6x1cwAAAADanDYd7hQDnqjTM3z48Fo3BQAAAKDNafPhzqKLLpr9fO2112rdFAAAAIA2p82HO//973+zn99++22tmwIAAADQ5rTpcCcCnXPOOSdb7tOnT62bAwAAADBnz5ZV7IXTkOnTp6fPPvssjR49Ol1yySXppZdeSh06dEh77LFHVdoIAAAAkCdVDXeWW265Zt8miin3798/HX300RVpEwAAAECeVXVYVgQ1zbksuOCC6YQTTkgPPPBA6tKlSzWbCgAAAJALVe25c+211za6TseOHdO8886b9fJZffXV01xzzVWVtgEAAADkUVXDnb333ruadwcAAADQ7rXp2bIAAAAAaJhwBwAAACDHqjosa2bPPPNMViw5pjv/9NNPs78ttNBCWa2dH/7wh6lfv361bB4AAABAm1eTcOfFF19MBx54YHrqqafqXefXv/51Wm+99dKVV16Z1lhjjaq2DwAAACAvqj4sK3rqrLvuulmwU5zyvFOnTmnRRRfNLrFc/PsTTzyRrfuvf/2r2s0EAAAAyIWqhjsff/xx+slPfpKmTJmSOnTokAYPHpyefPLJ9PXXX6f33nsvu3zzzTdZ8HPAAQdk06DHunGbTz75pJpNBQAAAMiFqoY7F154Yfriiy9S586d0913352uuuqqtM4662S9dYoi0Fl77bWz4Vixztxzz53dJm4LAAAAQA3DnQhrosfOYYcdlrbeeutG199qq63S4Ycfng3RitsCAAAAUMNw56233sp+7rDDDk2+TXHdN998s2LtAgAAAMirqoY7kydPzn7OM888Tb5Ncd2ovQMAAABADcOdxRZbLPv53HPPNfk2xXVjJi0AAAAAahjubLzxxln9nN///vdp4sSJja7/5ZdfprPPPjur0xO3BQAAAKCG4c5BBx1Uqr2zySabpNGjR9e7blw3YMCANG7cuBluCwAAAMD/7/+fg7wKNtxww3TIIYekyy67LL344otpvfXWS6uttlr2s1evXlkPnQ8++CA9+eST6eWXXy7dLm4TtwUAAACghuFOuPjii1P37t3TH/7whzR9+vT00ksvzRDkhBi6FTp27Jh+9atfZcO4AAAAAKjxsKwQvXPOOeec9J///CcdfPDBacUVV8zCnPJL/C2ui3WKNXcAAAAAaAM9d4pWX331dOmll2bLU6dOTZ999lm2vOCCC6bOnTvXqlkAAAAAuVKzcKdchDmmOgcAAABoY8Oy7r333vSDH/wgu9xwww3Num2sX7ztAw88ULE2AgAAAORZxcKdqJ1z9NFHp+effz4tssgi6ac//Wmzbr/nnnumnj17ZnV3jjnmmEo1EwAAACDXKhbu/Pvf/05jx47NZrw6//zzm337KKJ8wQUXpLnmmiubUevhhx+uSDsBAAAA8qxi4c7f/va37OeWW26ZVl111dnaRtxu6623zpZvueWWVm0fAAAAQHtQsXDnqaeeynrfbL/99i3aznbbbZcN8XriiSdarW0AAAAA7UXFwp233347+7nSSiu1aDt9+/bNfo4fP75V2gUAAADQnlQs3Pniiy+ynwsttFCLtlO8/cSJE1ulXQAAAADtScXCnfnmmy/7+fnnn7doO8XbzzvvvK3SLgAAAID2pGLhTkx/Hl555ZUWbWfMmDHZz169erVKuwAAAADak4qFO+uuu25WCPnOO+9s0XbuuOOOrDDzOuus02ptAwAAAGgvKhbubLPNNtnP+++/P40cOXK2tvHII49kty/fHgAAAABVCHd23XXXtOyyy2a9d37yk5+k119/vVm3Hzt2bNptt92yXjuxnR//+Met2r4PP/ww3XXXXemUU07JgqOePXtm9xWXffbZp0nbGD58eOk2jV1i3cZ888036Zxzzsl6KUUh6XnmmSetvPLK6ZhjjinNPgYAAABQrlOqkLnnnjude+65WSgTQUq/fv3SGWeckQYPHpyFFvX56quv0tVXX52FLrEcwch5552XOnVq3aYuuuiiqS1544030sCBA2cJwV577bXsEo/Jn//857TddtvVrI0AAADAHBTuhF122SX95je/Saeeemr6+uuv0y9/+ct08sknp4033jgLe6JIcgQ9cd0HH3yQnn322fToo49mv0ePnxC332mnnSrZzLT00ktnPWSKQ8Bmx3333ZeWWGKJeq/v3bt3vdd9+eWXadttty0FOwcccEDaY489Urdu3dKDDz6YzjrrrGwq+N133z099thj6fvf//5stxMAAABoXyoa7oQIcyLYOPzww7NhR9Eb5x//+Ed2qUsx1OnevXu65JJLmjxEqrmiZ1AMf4pL9OIZP358Wm655WZ7e3379s2Gj82OoUOHZsPQQgzLGjJkSOm6/v37p0033TQNGDAge/yOOuqo9NBDD812OwEAAID2pWI1d8rtu+++WXgRPXeitk0EOPVd4vqoMRPrVyrYKfYIiiFOtR6e9e2336aLLrooW15llVWy/31mG2ywQdp///2z5Ycffjg9/fTTVW8nAAAAMIf23CmKIUtRgycuL7/8cnr++efTJ598kg1JmnfeedPCCy+c1lxzzbTaaqulOUkMu/riiy+y5b333jt17Fh33hZB15VXXpkt33bbbaaGBwAAAKob7pSLAGdOC3HqUz5NfAy9qs/aa6+dDVWLoVlRdwcAAACgasOy5gQx9Cx6J3Xu3DkbWrb++uunk046Kb377rsN3u6VV14pLUdR5/rEbGF9+vTJlseMGdOKLQcAAADyrCY9d9qj8iLHMdwsLk8++WQ2jfsFF1yQDjrooDpvN2HChOxnzBq2wAILNHgfSy21VHrhhRfSRx99lKZMmZK6dOnS5PYV76c+77//fpO3BQAAALQdwp0WWn755bMp32NWqwhfwptvvpn+9re/pVtuuSVNnjw5/eIXv0gdOnRIBx544Cy3j5pDoUePHo3eVwRARTHrWHPCnWLbAAAAgPZFuNMCO++8c1YEOYKbclHsePfdd0933XVXFvzEjFhHH3102mGHHdJiiy02w7oR/oQYztWY8jBn0qRJrfZ/AAAAAPml5k4LzD///LMEO+ViqvVTTjklW45CyMOGDZtlna5du2Y/p06d2uj9xVCsom7dujWrre+8806Dl6eeeqpZ2wMAAADaBuFOhcVQrGIA9PDDD89yfUwDXxxm1Zivv/66tNyUYVzlevfu3eBl8cUXb9b2AAAAgLZBuFNhvXr1SgsvvHC2XNfMWRGsFIObzz//vMFtRQ+bsMgiizSr3g4AAADQfgl3qqChoVurrrpqafnVV1+td71p06alcePGZcurrLJKK7cQAAAAyCvhToXFtOUff/xxtrzEEkvMcv1GG21UWq5r2FbR6NGjS8OyNtxww4q0FQAAAMgf4U6FXXXVValQKGTLAwYMmOX6TTfdNCvMHK677rrSujMbPnz4DLN0AQAAAAThzmwaP358eu655xpcJ6ZCP/3000uzW+27776zrBNToB9xxBHZ8pgxY9K55547yzqjRo0qzbQVAVFMtQ4AAAAQOs2pD8PIkSPTG2+8Ufq9OHQqxN/Le8qEffbZZ5ZwZ7PNNkv9+/dP22+/fVpzzTWz4snhzTffTLfcckt2KfbEidBmySWXrLMtQ4YMSX/5y1/S2LFj07HHHpvd/x577JEFQg8++GD63e9+l9Xcid8vuOCCVn0cAAAAgHzrUKhvHFA7F2FNDINqqpkfpoceeigLdxrTvXv3dP7552dTojckAp2BAwem119/vc7r55tvvvTnP/85bbfddqkSJkyYkJZaaqnSrFzFWbz6DRmRaumZoYNqev8AAABQje/fLTHH9txpqX79+qXrr78+GzIVxY7ff//9rPdP9LBZcMEF02qrrZa22GKLNHjw4FKPnob06dMnG+Z16aWXpptvvjkLe6ZOnZo94RH6HHnkkWmZZZapyv8GAAAA5Mcc23OHGem5AwAAAPnsuaOgMgAAAECOCXcAAAAAcky4AwAAAJBjwh0AAACAHBPuAAAAAOSYcAcAAAAgx4Q7AAAAADkm3AEAAADIMeEOAAAAQI4JdwAAAAByTLgDAAAAkGPCHQAAAIAcE+4AAAAA5JhwBwAAACDHhDsAAAAAOSbcAQAAAMgx4Q4AAABAjgl3AAAAAHJMuAMAAACQY8IdAAAAgBwT7gAAAADkmHAHAAAAIMeEOwAAAAA5JtwBAAAAyDHhDgAAAECOCXcAAAAAcky4AwAAAJBjwh0AAACAHBPuAAAAAOSYcAcAAAAgx4Q7AAAAADkm3AEAAADIMeEOAAAAQI4JdwAAAAByTLgDAAAAkGPCHQAAAIAcE+4AAAAA5JhwBwAAACDHhDsAAAAAOSbcAQAAAMgx4Q4AAABAjgl3AAAAAHJMuAMAAACQY8IdAAAAgBwT7gAAAADkmHAHAAAAIMeEOwAAAAA5JtwBAAAAyDHhDgAAAECOCXcAAAAAcky4AwAAAJBjwh0AAACAHBPuAAAAAOSYcAcAAAAgx4Q7AAAAADkm3AEAAADIMeEOAAAAQI4JdwAAAAByTLgDAAAAkGPCHQAAAIAcE+4AAAAA5JhwBwAAACDHhDsAAAAAOSbcAQAAAMgx4Q4AAABAjgl3AAAAAHJMuAMAAACQY8IdAAAAgBwT7gAAAADkmHAHAAAAIMeEOwAAAAA5JtwBAAAAyDHhDgAAAECOCXcAAAAAcky4AwAAAJBjwh0AAACAHBPuAAAAAOSYcAcAAAAgx4Q7AAAAADkm3AEAAADIMeEOAAAAQI4JdwAAAAByTLgDAAAAkGPCHQAAAIAcE+4AAAAA5JhwBwAAACDH5thw58MPP0x33XVXOuWUU9I222yTevbsmTp06JBd9tlnn2Zv7957700777xz6t27d+rSpUv2M36PvzfVtGnT0hVXXJE23njjtMgii6Ru3bqlFVZYIR100EHp5ZdfbnabAAAAgPavU5pDLbrooq2ynenTp6cDDzwwDRs2bIa/v/vuu9nl9ttvT4MHD05XXnll6tix/izt448/TgMHDkxPP/30DH9/880301VXXZWuu+66dMkll2TbAgAAAEhzes+dcksvvXTaaqutZuu2J554YinYWWuttdKNN96Ynnrqqexn/B6uvvrqdNJJJ9W7je+++y7r5VMMdnbZZZesx8+TTz6ZLrrootSrV680ZcqUrAdPc3oCAQAAAO3fHNtzJ4ZjrbPOOtklevGMHz8+Lbfccs3axtixY9O5556bLa+99trpkUceyYZShdjuDjvskAYMGJBGjx6dhg4dmvbbb7/Up0+fWbYTvXJGjhyZLR9yyCHp0ksvLV237rrrZsPG+vXrlyZOnJiOOOKINGbMmNSp0xz71AEAAABl5tieO7/5zW/Sdttt16LhWRdccEFWJydcfPHFpWCnqHv37tnfQ6x3/vnn17mdYkC00EILZSHQzCIQOuGEE7LlN954I912222z3WYAAACgfZljw52WKhQK6Y477siWV1555bT++uvXuV78faWVVsqWY/243cy9f6InTthtt92yQKgu5UWehTsAAABAkXBnNr311lvpvffey5Zj6FVDitdHgeUY/lWuOByrse0stthiqW/fvtnyY4891qK2AwAAAO2Hwi2z6ZVXXiktR8+dhpRfH710ymv7NHc70dPnnXfeSV9//XWaZ555mtzeCRMmNHj9+++/3+RtAQAAAG2HcGc2lYclvXv3bnDdpZZaqrQcwUxLtxNDu+J2xeFeTVHeBgAAAKD9MCxrNn355Zel5R49ejS4bnkPm6+++qoi2wEAAADmTHruzKbJkyeXljt37tzgul26dCktT5o0qSLbaczMPYbqGpYV064DAAAA+SLcmU1du3YtLU+dOrXBdadMmVJannm69Jm3U/57c7bTmMaGfAEAAAD5ZFjWbJp33nmbPEQqih/XN/SqtbYDAAAAzJmEO7OpvCdMYzNRlQ+Jmrmw8exsp0OHDnriAAAAABnhzmxaddVVS8uvvvpqg+uWX7/KKqu0eDsREDVnGnQAAACg/RLuzKblllsuLbHEEtnyww8/3OC6jzzySPZzySWXTMsuu+wM12200Ual5Ya287///S+NHTs2W95www1b1HYAAACg/RDuzKYYGrXjjjuWetQ88cQTda4Xfy/2uIn143bl+vbtW+rN89e//jV98803dW5n+PDhpeWdd9651f4PAAAAIN+EOy1w1FFHpbnmmitbPvzww2eZnjx+j7+HTp06ZevX5Ve/+lX289NPP03HHnvsLNePGzcunXXWWdlynz59hDsAAABAyRw7FfrIkSPTG2+8Ufr9448/Li3H38t7yoR99tlnlm1Er5shQ4ak3//+92n06NHZcKnjjjsurbDCClkgc/bZZ6fnnnsuWzfWW3HFFetsy957752uueaa9Nhjj6VLL700G4J1wAEHpAUXXDA99dRT6YwzzkgTJ05MHTt2TBdddFEWFAEAAACEDoVCoTAnPhQR1lx33XVNXr++h2n69OlZEBPhTH3233//dNVVV2XhTH0iXBo4cGB6+umn67y+S5cu6ZJLLkmDBw9OlRAzdRVn8opZuYqzcfUbMiLV0jNDB9X0/gEAAKAa379bwrCsForAZtiwYenuu+/OaupEkeXOnTtnP+P3e+65J1199dUNBjuhZ8+e6fHHH0+XXXZZVmR54YUXTl27dk3LL798Fh4988wzFQt2AAAAgPyaY3vuMCM9dwAAAKDy9NwBAAAAYAbCHQAAAIAcE+4AAAAA5JhwBwAAACDHhDsAAAAAOSbcAQAAAMgx4Q4AAABAjgl3AAAAAHJMuAMAAACQY8IdAAAAgBwT7gAAAADkmHAHAAAAIMeEOwAAAAA5JtwBAAAAyDHhDgAAAECOCXcAAAAAcky4AwAAAJBjwh0AAACAHBPuAAAAAOSYcAcAAAAgx4Q7AAAAADkm3AEAAADIMeEOAAAAQI4JdwAAAAByTLgDAAAAkGPCHQAAAIAcE+4AAAAA5JhwBwAAACDHhDsAAAAAOSbcAQAAAMgx4Q4AAABAjgl3AAAAAHJMuAMAAACQY8IdAAAAgBwT7gAAAADkmHAHAAAAIMeEOwAAAAA5JtwBAAAAyDHhDgAAAECOCXcAAAAAcky4AwAAAJBjwh0AAACAHBPuAAAAAOSYcAcAAAAgx4Q7AAAAADkm3AEAAADIMeEOAAAAQI4JdwAAAAByTLgDAAAAkGPCHQAAAIAc61TrBgDMrN+QETW9/2eGDqrp/QMAADSHnjsAAAAAOSbcAQAAAMgx4Q4AAABAjgl3AAAAAHJMuAMAAACQY8IdAAAAgBwzFTpAO5uqva23DwAAaF167gAAAADkmJ47APD/6PUEAEAe6bkDAAAAkGPCHQAAAIAcE+4AAAAA5JhwBwAAACDHhDsAAAAAOSbcAQAAAMgx4Q4AAABAjnWqdQOA2ug3ZERN7/+ZoYNqev8AAADthZ47AAAAADkm3AEAAADIMeEOAAAAQI4JdwAAAAByTLgDAAAAkGPCHQAAAIAcE+4AAAAA5JhwBwAAACDHhDsAAAAAOSbcAQAAAMgx4Q4AAABAjgl3AAAAAHJMuAMAAACQY8IdAAAAgBwT7gAAAADkmHAHAAAAIMeEOwAAAAA5JtwBAAAAyDHhTgt16NChSZdNN9200W3de++9aeedd069e/dOXbp0yX7G7/F3AAAAgLp0qvOvVNX06dPTgQcemIYNGzbD3999993scvvtt6fBgwenK6+8MnXsKI8DAAAA/n/CnVZy8MEHp0MOOaTe6+eZZ556rzvxxBNLwc5aa62Vjj322LTCCiukcePGpXPOOSc999xz6eqrr06LLLJI+t3vfleR9gMAAAD5JNxpJb169Uqrr756s283duzYdO6552bLa6+9dnrkkUdSt27dst/XWWedtMMOO6QBAwak0aNHp6FDh6b99tsv9enTp9XbD1At/YaMqNl9PzN0UM3uGwAAKsUYnxq74IIL0rRp07Lliy++uBTsFHXv3j37e4j1zj///Jq0EwAAAGibhDs1VCgU0h133JEtr7zyymn99devc734+0orrZQtx/pxOwAAAIAg3Kmht956K7333nvZcgy9akjx+iiwPH78+Kq0DwAAAGj71NxpJTfffHP661//mgUvc801V1psscXSBhtskPbZZ5+02Wab1XmbV155pbQcPXcaUn79mDFj0nLLLdes9k2YMKHB699///1mbQ8AAABoG4Q7raQ8qAlvvPFGdhkxYkTaaaed0vDhw9P8889fb+DSu3fvBre/1FJLlZbfeeedZrev/PYAAABA+yHcaaEoeBwzWm2xxRZZ75oePXqkjz76KD388MPpiiuuSJ988km6/fbb04477pj++c9/prnnnrt02y+//LK0HLdrSPlU6l999VWF/hsAAAAgb4Q7LRQ1cBZYYIFZ/r7lllumww8/PG2zzTbpueeey8Keyy+/PB1xxBGldSZPnlxa7ty5c4P306VLl9LypEmTmt3Oxnr7xLCsddddt9nbBQAAAGpLuNNCdQU7RYsuumi65ZZbsh493377bTaleXm407Vr19Ly1KlTG7yfKVOmlJZnni69KRob9gUAAADkk3CnwpZffvmsF88999yT1eCJ2bGWWGKJ7Lp55523yUOtvv766yYP4QKgfeo3ZERN7/+ZoYNqev8AANTNVOhVsOqqq84wjKuu3jSNzWZVPqxKcWQAAACgSLhTBR06dGg09Hn11Vcb3Eb59ausskortg4AAADIM+FOladJLw7JCsstt1zp9yi43JBHHnkk+7nkkkumZZddtmJtBQAAAPJFzZ0Ke+utt7Ip0MMKK6yQhTPlPXpiivSYRSt65jzxxBNp/fXXn2Ub8fdiz51Yv76eQLQtamMAAABQDXrutMCdd96Zpk2bVu/1H3zwQdp1111LM2Edcsghs6xz1FFHpbnmmitbjqnTZ57mPH6Pv4dOnTpl6wMAAAAU6bnTAhG6xBTnEeD0798/Gy4V05R//PHH6aGHHkpXXnllthw22mijdOihh86yjb59+6YhQ4ak3//+92n06NFpww03TMcdd1zWy2fcuHHp7LPPTs8991y2bqy34oorVv3/BAAAANou4U4LxdTmF198cXapT4Q/V199derSpUud1//2t79NH374YbrmmmuyIGePPfaYZZ39998/nXnmma3adgAAACD/hDstcN1112WFkEeNGpXefPPNrJfOxIkTU48ePbLpyjfYYIO09957Z716GtKxY8c0bNiwLAS66qqr0tNPP51tq2fPnmmdddZJBx10UNpmm22q9n8BAAAA+SHcaYEBAwZkl9YycODA7AIAAADQVAoqAwAAAOSYcAcAAAAgx4Q7AAAAADkm3AEAAADIMeEOAAAAQI4JdwAAAAByTLgDAAAAkGPCHQAAAIAcE+4AAAAA5JhwBwAAACDHhDsAAAAAOdap1g0AANqHfkNG1Oy+nxk6qGb3DQBQa3ruAAAAAOSYcAcAAAAgx4Q7AAAAADkm3AEAAADIMeEOAAAAQI4JdwAAAAByTLgDAAAAkGPCHQAAAIAcE+4AAAAA5JhwBwAAACDHhDsAAAAAOSbcAQAAAMgx4Q4AAABAjgl3AAAAAHJMuAMAAACQY8IdAAAAgBwT7gAAAADkmHAHAAAAIMeEOwAAAAA5JtwBAAAAyDHhDgAAAECOdap1AwAAKq3fkBE1vf9nhg6q6f0DAO2bnjsAAAAAOSbcAQAAAMgx4Q4AAABAjgl3AAAAAHJMuAMAAACQY8IdAAAAgBwT7gAAAADkmHAHAAAAIMeEOwAAAAA5JtwBAAAAyLFOtW4AAABtW78hI2p2388MHVSz+waAvNBzBwAAACDHhDsAAAAAOSbcAQAAAMgxNXfIrVqO/w9qAAAAANAW6LkDAAAAkGN67gAA1JjeqABAS+i5AwAAAJBjwh0AAACAHBPuAAAAAOSYcAcAAAAgx4Q7AAAAADkm3AEAAADIMeEOAAAAQI4JdwAAAAByTLgDAAAAkGPCHQAAAIAc61TrBgAAwOzqN2RETe//maGDanr/ABD03AEAAADIMeEOAAAAQI4JdwAAAAByTLgDAAAAkGPCHQAAAIAcM1sWAABUiNm8AKgGPXcAAAAAcky4AwAAAJBjwh0AAACAHBPuAAAAAOSYcAcAAAAgx4Q7AAAAADkm3AEAAADIMeEOAAAAQI51qnUDAAAAZtZvyIia3v8zQwfV9P4BmkPPHQAAAIAcE+4AAAAA5JhwBwAAACDHhDsAAAAAOSbcAQAAAMgx4Q4AAABAjgl3AAAAAHJMuAMAAACQY51q3QAAAKA2+g0ZUbP7fmbooJrdN0B7o+dOG/P222+nY445Jq288sppnnnmSQsttFBaZ5110tChQ9M333xT6+YBAAAAbYyeO23InXfemfbaa680ceLE0t8i0Bk9enR2ufrqq9Pdd9+d+vTpU9N2AgAAAG2HnjttxHPPPZd23333LNjp0aNH+u1vf5sef/zx9K9//SsdcMAB2Tpjx45N2267bfryyy9r3VwAAACgjdBzp4048sgj06RJk1KnTp3S/fffn/r371+6bvPNN08rrrhiOvbYY7OA57zzzkunnXZaTdsLAABzslrWKwpqFgHl9NxpA5566qn06KOPZsv777//DMFOUdThWWWVVbLlCy+8MH377bdVbycAAADQ9ui50wbcfvvtpeV99923znU6duyYBg0alE444YT0+eefpwcffDBttdVWVWwlAABAy+n1BK1Pz502YOTIkdnPmB2rX79+9a43YMCA0vJjjz1WlbYBAAAAbZueO23AmDFjsp8xC1bU3KlPTI8+822aasKECQ1e/84775SW33///dLy1C8/TbXUULvbctuC9jXMczv7tK99ti1oX/tsW9C+9tm2oH3ts23toX3bnHlLqpV7T/pxg9d77Fr2+LXl59Zj17T2lX/nnjZtWmoNHQqFQqFVtsRsmTx5curWrVu2HDNh3XXXXQ2uHzNpff3112n99ddPo0aNavL9dOjQocVtBQAAAFq3Bu8666zT4u0YllVj5dOaR3DTmBi6Fb766quKtgsAAADIB8Oy2kDPnaLOnTs3un6XLl2ynzFtenOUD7uqrx2vvvpqWnTRRdMiiyzS4PCwpohuZuuuu24piVx88cVTW9GW2xa0r322LWhf+2xb0L722bagfe2zbUH72m/72nLbgva1z7YF7WufbatE+2Io1kcffZQtr7HGGqk1CHdqrGvXrqXlqVOnNrr+lClTsp/FoVxN1bt370bXiZo/lRAv/Kbcfy205bYF7WufbQva1z7bFrSvfbYtaF/7bFvQvvbbvrbctqB97bNtQfvaZ9tas33LLrtsak2GZdXYvPPOW1puylCrqLfT1CFcAAAAQPsn3GkDPXcWXnjhJlVt/+yzz0rhzlJLLVWV9gEAAABtm3CnDVh11VWzn2+88UaD06BFTZyiVVZZpSptAwAAANo24U4bsNFGG2U/o1fOM888U+96Dz/8cGl5ww03rErbAAAAgLZNuNMG7LTTTqXla6+9ts51pk+fnkaMGJEtL7DAAmmzzTarWvsAAACAtku40wbElGobb7xxtjxs2LA0atSoWdY577zz0pgxY7LlI488Ms0999xVbycAAADQ9pgKvY248MILs6FWkyZNSltttVX69a9/nfXOid9vuummdNVVV2Xr9e3bNx1zzDG1bi4AAADQRnQoFAqFWjeC/3PnnXemvfbaK02cOLHO6yPYufvuu1OfPn2q3jYAAACgbRLutDFvv/121osnQpyYGr1z585ZmPOTn/wkHXbYYal79+61biIAAADQhgh3AAAAAHJMQWUAAACAHBPuAAAAAOSYcAcAAAAgx4Q7AAAAADkm3AEAAADIMeEOAAAAQI4JdwAAAAByTLgDAAAAkGPCHVrd22+/nY455pi08sorp3nmmScttNBCaZ111klDhw5N33zzTU3a9OGHH6a77rornXLKKWmbbbZJPXv2TB06dMgu++yzT6ql0aNHp9NPPz1ttdVWqXfv3qlLly6pR48eqW/fvmnfffdNI0eOrFnbJk6cmG666abs+RwwYEDq06dPmn/++VPnzp1Tr1690qabbprOOeec9Mknn6S25rjjjis9x3F56KGHqt6G8vtv6BKPY1vw3//+N5166qlp7bXXTossskjq2rVrWmqppdLGG2+cvXdeeumlqrUlHpOmPn61fI6nTp2arr766rT11lunxRdfvPT+XWmllbL37+OPP55qafLkyemyyy5LW2yxRfacxnt3iSWWSAMHDsze23nY3957771p5513Lu0f42f8Hn+vVfumT5+eXnnllTR8+PB0yCGHZMe4aFtLX4ut0bY4zt56663p4IMPztq14IILprnnnjstvPDCqX///um0005L//vf/2rWvjFjxqRLLrkk7b333ukHP/hB9nzGviY+Lyy//PJp9913T3fccUcqFAo1aV9Dj2u0r7i9ZZddtupti9dbU/eHsW612zezBx54ILttfHaI5zc+P8Rnmx//+Mfp8ssvT1999VXV2jd+/PhmH1Oa+hy35mMX7YzPL/369UsLLLBA9t6Nz9EbbLBB9lkx7qu5WrN9b731Vjr66KPT6quvnuadd97seV1xxRWz/eDLL7+c2sLn4NY8ZrRG2yp1vGit9lXymNEa7avUMWN0Bb9/tfR40SIFaEV///vfC/PNN1+8u+q89O3bt/D6669XvV31tScue++9d6FWNt544wbbVrwMGjSoMGXKlKq375///GeT2tezZ8/CP/7xj0Jb8dxzzxU6deo0QxsffPDBqrejKY9dXAYMGFCotYsuuqgwzzzzNNjOI488smrticekqY9fXDp27FiYMGFCoZrGjx9fWG211Rpt2+GHH16YPn16odpeffXVwkorrdRg27baaqvCl19+2Sb3t999911h//33b3B7gwcPztardvuGDx/e4HZmd3/T0rY9//zzhR49ejT6mozj9E033VT19oWf/exnTd4vfvzxx1VvX32OOeaYGba3zDLLVL1t1157bZP3ibFutdtX9OmnnxZ23HHHRtsYx+pqte+tt95q8mNXvn+sRtuKRowYUejWrVuD21tooYUK999/f5O32Zrtu/LKKwudO3eud1tx3cUXX1yzz8GtfcxorbZV6njRGu2r5DGjtR6/ShwzNq7w96+WHi9aolP1YiTau+eeey5LTydNmpQlnyeccELabLPNst/jDPEf//jHNHbs2LTttttmaWkk/rWw9NJLZ72K7r///lRr7733XvYzzqT/5Cc/yXpIRPu+++67NGrUqHTeeeeld999N40YMSJ9++236YYbbqh6G6PnRjyPcRYplqN3QpyFmDBhQrrllluytP/jjz9OO+ywQ3rqqafSmmuumWop2nbggQemadOmZb2LZucsV2uLsyFxtqY+cfahls4888x08sknZ8txxuKAAw7Izt7EWdbolRXv7dtuuy117Fi9zp7XXntt+vrrrxtcJ86ExT4nRM+UJZdcskqtS9n7MfZlxTOV3/ve99Ivf/nLrMfOl19+mZ3xifdv/A8XX3xx9h4//vjjq9a+eN1vueWW6Z133sl+j/1LnPWKdsR+57rrrks333xzth/cY489srO6bW1/e+KJJ6Zhw4Zly2uttVY69thj0worrJDGjRuX9RiM12X0mooeSb/73e+q2r7yM4RxhnONNdbIXhMvvvjibLWjtdoWvS2LvSE23HDDtN1222U98eIM7EcffZTtr+NYHOv97Gc/S/PNN192Jr9a7QudOnVK6623Xta+eNwWW2yx7Dn87LPP0quvvpquvPLKrJfgww8/nLbffvvsvTQ7+57WPNbHa+2CCy7IzhbH8x3v8ZZojbbdd9992fu5PnEmuhbt++KLL7J9zzPPPJP9Hj0moqdOvHfnmmuubJ8Uz+3f/va3qrYvjg9NeX+eddZZpc9asc+sRtvCY489lvWiic8w8XqP+95xxx2z5zh61cY++84770yffvpp9vd4j0TPgGq1Lz7HH3TQQdlyfDaIHt2bb7551tsh3h+xT37jjTfSEUcckX322m233ar+Obi1jxmt1bZKHS9ao32VPGa01uNXiWPGexX8/tXax4tmq1qMRLtXTEGjx8Tjjz8+y/XnnHNOKcE89dRTq9q2U045pXDnnXcW/ve//81yBqeWPXe23Xbbwl/+8pfCtGnT6rz+o48+yno7Fdv68MMPV7V99bWr3G233VZq384771yotfPPPz9ry8orr1w44YQTWnxmpCVq9XpvjgceeGCGMxRTp06td91a9B5ryLHHHltq+5/+9Keq3vfNN99cuu/+/fvX+V4ZPXp0Ye65587WWWCBBQrffvtt1dp36KGHNvr6i/1icZ34f9rS/va1114r9b5be+21C998880M13/99dfZ34vHnOb0CG2N9j355JNZb7dRo0YVJk2alP0tHueW7m9a2rbHHnussNtuuxVefvnlete5/fbbCx06dMi2ucIKKzSrV1lrPHaNvQ/ivbTLLruUtnvHHXdUtX11tadfv37ZNk4//fTsDOzsnIltjbaV99yJ27em1nrsfv7zn2e36dKlS4PPXbzumrNPrMbnuHiul1hiiWyb88477yz7nUq2LT4PFm9z6aWX1rnOL3/5y9I6sY9vqpa2L/a3vXr1ytaPXh4vvvjiLOt88cUXhTXWWCNbZ9FFF21yj9DW+hxciWNGa7WtUseL1mhfJY8ZrfX4VeKYsW2Fvn+11vGiJYQ7tIrYcRXfAAcddFCd60Q3yFVWWaX0RaehL5GV1lbCnaaIA3L58I62qDj0I4Zn1dLbb79d6l760EMPtcrBsz2HO/GeXHHFFbM2rrnmmlUNH1qj7UsuuWTpw2Z8cKumo48+uvT8xnDU+kTgWVzvhRdeqErb4sPF/PPPX/pAUd+Hl/j70ksvna0XH0ba0v724IMPLt0mPhDXJf5eXOeQQw6pavvqUon9TaWOVbvuumtpu88880yba1/5c/urX/2qpu0777zzstvHcS4C7tb6sN7Wwp3WaN+jjz5aus3QoUPbXPsaE8PLi9vcd999q9q2BRdcMFt/4YUXrnedzz//vLTdH/zgB1VrX/nJjBNPPLFJQ/mbOzyrpZ+Dq3nMaG7b6lKtz6et9R2itY4ZlWpfax0zWtq2Sh0vmkNBZVrF7bffXlqOIlR1iS5ygwYNypY///zz9OCDD1atfXkWQ6KKomtpW1QcYhfFW2vp0EMPzbqXRnfmKABNw6Jb9uuvv54tRwHH6PqaF//617+yLrMhuvx379696oWUixrqGh9dwuu6TSXFcxpDI0IMj4ihEHWJv8f1IYZQRKHMtiBy0SiOGGL4wPrrr1/nevH3GAYXZrcA75yqrR9Xyodt1/K4EhNERBHacMUVV2QFyalfFD0tDts57LDDUt7EEIyi2RmS1RLF48Nyyy1X7zrxuEYx5PL1qyFKKRQ1NCQnJkKIoSghhu1Xa39Vy2NGW9+Xtlb7KvV/ttZ2K3HM2KyZbWsrxwvhDq2iWFE8aodEbZb6lH/hjvHFNG7KlCml5fq+pNXSa6+9lv7zn/+UDqq18te//jWrGxKzSpx77rk1a0eeRM2VEJX8Y5x1UYzpj4AgfubhQ3gxNK6m4gfE8Oabb9a7XvEDQTzGMaNINZTPXrfooos2uG759Y8++mhqCyJkKo6HbyykLV4fQV/MMkP7OK6Uz+RWy+NK1EqLulk///nP28yshm1VhA3FL9gRGhe/5EcNi6izE+/PWp8AakjUxSieqIyZbTbZZJOaHFMaCtmj7knUOCxfvy0dU+IEUXwGC1G3JGofVmN/VctjRlvfl7ZW+yr1f7bWditxzJjSzLa1leOFcIdWEdPUhZjusqGz/+VvuOJtaFgUCCtaZZVVUlsQU/zFl/8//OEP2YGyeAA/6qijatKe6Al25JFHZstnn3126cxWWwpRVl111ax3SZxdiC/5cVaw1r3XnnjiidIH2WhXFIyLYnVRSC8KK8fP+AAZYVn5Qa7WondWFHgOyyyzTE0OonvuuWdWWLD4mosvMHUV1bv77ruz5Z/+9Kel9SstCtoXFXvw1Kf8+ihQ3RaUt6OxD2mOKe3nuBJfWuML4f77759++9vfZn+LfXkU8ayF+LJwzz33ZNMCR3HNtiZ6SUcx0Dg7HI9T9Eo46aSTSj0aq+35558vhTdxHIkgIj4TRNuiUGn0SImeJxH8tGTq50qJnibx2SbEl7MI5KvpF7/4RSlIibP+dTnjjDNmWb8tHVOiJ0w878WwLwosV2N/VctjRlvcl1aifZX6P1uy3UofMx5uRtva0vEiP33wabPiYF48k9DY7Azxoo/ePZFsFmdxoX4xa8Lvf//70u9NnX2gEoYPH17vkLsQMwHFF9haiBkR/ve//2WV9GMn39bM/KU5PvDEJXqf7LTTTtljGx96q/3ailkGigfDCMcuuuiiWdaLGe6GDBmShSkRVCywwAKp1mKWleJMWnvttVfVP4QXH7M//elPWcgTvRBjdrH4IhOhWIRP8bc4wMcH3B/84AdVPdhHyB4zNMQMD4888kiD65ZfHzOytAUxE19RY8eUmMGvyDGl6V/Ci6FjfAmv5ReSCGbLP0DP/B6L/U4t9jkxC0vxZEUcg2NmlramPCCJQCAuTz75ZLaviZlaijMb1eI4F8eXmHGnOOy3KPaHDzzwQDasNmaliuHAbUWte4Put99+WS/4aEcMMY+hsjELacxQGvvmON4UexbFrFA//OEPq9a28n1EvF/r66EfJzSKMy+FaHdLe1E05XNwrY4ZbekzeiXbV6ljxuy0r1rHjOnNaFtbO17ouUOLlU/xVp7uNzbtc/kBgLqdf/752fTiYZdddmlwyFutfP/738/aGB/UavElO4aSxNSW0WMsznbVog31iZ46Mc10TCMZ7YwPPlHnJj6YRa+YEB/WYlrT+CJeTXH2LQ5eIabjjGAnPkRef/312XCsOIMZB9Di2PXHH388+/DZFtT6Q3hRfPCOD+CDBw/OhiZGb6z+/ftnZ6ZPO+207PmPL1nx3Dc2PKo1xT42pqgNL7zwQrrxxhvrXC/+Xj4Va9Wn62yFY0rxeBIcUxoXPfDi9VrsaVY829nWxHTKcVZ9o402qsn9R6D9wQcfZO/nAw44ILUlUePrV7/6VRZyx7E3LnHWOKbzjeNfnHCLXh1XXXVVVdtVPow3ejNGsPOjH/0oa1+06cMPP0yXX355diIjenjECaHiMK5aixCi+IVxgw02yALyaothHzHdefT0XXPNNbPPNXGMiRMHu+66a/ZZIWqA/POf/0xnnnlmVdsWdXaKvfKjx3bxhG65+DwRn23KtcYxpSmfg2t1zGjrn9Fbo32VPGa05uPX2seM85vRtjZ3vKha6Wbarf/+97+lauIxBWZjllpqqdJ0erWSh9myYran4rSOMQXlBx98UNP2fPbZZ9n0l3F56qmnCjfeeGNpJqB4LqOqfLVFJfqY8jzaMGTIkFmur/VsWfGY1SemJF1rrbVK7bvwwgur2rZ33nmndN9x6d69e+HVV1+dZb2YTjRm0iqu98QTTxRqKdrdsWPHrC3rr79+TdsSr78TTjihsMgii8zwWJZfYurV5kzl3Fpi5oji/iOmYz/jjDOy2eRilsL4Gb/H3zt37lxq6xZbbNEm9rcxfWhx/X/9618NrhvXF9eN/6ka7cvzbFmDBw9u1e21tH1vvvlmdkyJmeQeeeSRwh/+8IdsBr94j2+33XalqZur2b6Y8jam/Y33z/PPPz/L9bWcLStmS2poGuI4Dsf7urhPf//996vWvnj/le/7ttxyyzpn6osZtYr78JhBtanTKre0fQ357W9/W9rWFVdc0aJttaRtr7zySmH77bcv7btnvnTt2rWwxx57FCZMmFD19h122GGl28QU0TE9dkx/HlN7x/Fm6623zq4rP6b86U9/qsrn4GofM5rTtvpU+vNpa32HaO1jRkvbV+ljRnPbVq3jRXMId2ixDz/8sPTG33333RtdP94ose7qq69eqJW2Hu689NJLpWkx42AeO4+2asSIEdmOLXasMU1rNRUPjjGd81dffdXmwp3GjBs3rvRBvE+fPlW9748++miGD41HHHFEveveddddpfViCvBaOuuss0ptufzyy2vWjni9bbzxxlk75pprrsKxxx5bGDNmTBb4xAfe+++/v7DRRhtl18f7I6bHrLZhw4bV+yUhLt26dStccsklpd932mmnNrG/Peecc0rr33vvvQ2ue88995TWPffcc6vSvryGO7/73e9K21pnnXXq3GfWsn1F8WUxPqTHNuNkUAS61Wrf5MmTsylsY/1jjjmmznVqGe40RXnIcuaZZ1atfTH1efn+5dlnn6133R//+Mel9er6QlSJ9jWkeJKoS5cuDZ6UqWTb4ovq/PPPX3ptRTASX1QjkI/3wKWXXlpYaKGFsuuXWGKJ7HNiNdsX742BAwfWezwpnswon5I8AqBqfA6u9jGjNT6jV/LzaWt9h6jEMaM121eJY8ZLzWhbNY8XzWFYFi1WPv1cU7o4FmtlNGUI15woqv5vtdVW2RjO6KYb3a2rPWtDc0ThwegOHl1yY+rTas2wFPViYihYuPjii2foapsX0b2+OBV11OApzvZQ7fdtiNdcfbbYYotSl+ynn3461VLUHQhdunRJu+++e83aEcOuirNLDRs2LBuGELUForhpFE6O5zUKZkc3+jiREt12Y9x6NcUwuqjBsfPOO8/w/ojnMrr7P/vss1ldjPKaaHk7phSPJ8ExpX5XXnll+vWvf50tx+s0Cj+21X1mzLJ07bXXZsMaoyZG1FSrlhhyEDNARl2O3/zmNymPDjzwwNLw5PpqU1T6fRs1J9Zaa6161916661Ly7U+psTQi2L9udgv1qLGUwx9ifptMVx6scUWyyY7iHpyMZw36qdFHZmYiSdqpMX7Iz4rVHuq9jjm3nnnndkw8xiOXz4EvlevXtmQrDgmlk8vPrvHlOZ+Dq7mMaOtf0ZvrfZV6phRicevtY4ZbzWzbW31eKGgMi0Wb6qoHxLF/MqLmtUl3jDFHWt5UTP+Txywo0he/IwD5zXXXJPVY2nroo0xFXk8t//4xz+qUlg5xsNGccYISKI+TPk0iEUvvfRSafnf//53VnQ5bL/99m3mi03MohUHzRCznMTsJ9X6oBYfwD/66KNG34/xHo9CdfH4FdevhdGjR5eKdsbU7bUKI+LDa7w3QxRQru9DdoQoMbtJjAGP8DMKZ8frtpqimPOtt96azWj3/vvvZ++ZJZdcsjRNcdRYKlpttdVSW1BeELOxY0p5QUzHlFRvbaX4YlicXS5qdrS1GQVnFu2LAvnR1qjLEjXJ4ktupUVIG+I4HF9k61L8DBM/i8ed+HJbrHNVa9GW+EwWdVGqOXNW+fuvOUVta3lMaSs13OJzU/G5Ovzww7OApy6xj47QJ+rxRL23OGEQ9XmqpWPHjln9lbhEnZuoMxJfqKO9cV0oL6Idn2+q8Tm4WseMtv4ZvbXaV6ljRiUfv5YeM96bjba11eOFcIdWETvwSOyj90F8iahvOvTi2ZG2OmVgLcUHsTjb/+abb5Z6o9SyWGxzlFeGf/vtt6tyn8WpuePxijNejSmfQjTS+bYS7tSyAHR8UCzOuFLXVN7litfX996u9ofwap+1LBcfaIs91Bo6Ox3Ki/CV7/+qLZ63uj7IxheEonXXXTe1BeVfCBp7zBxTGvb3v/89O45EuBgF02OWosa+eLe140qE93F8jPZXWoSfIc4Cx6Uh0abisWfAgAFtJtyp1XGlPBxu6vGk1seU+AJY/oUrCkDXQvmU3BHIN3ZMiXCnuP+rZrgzc2+ZmXsAx/MakwuEOPHW3EBgdj8HV+OY0dY/o7dW+yp1zKjG4ze7x4yPZ7NtbfV4YVgWraJYnTySyfIvCzMr7yIcCSv/J7riRjflYq+EmEovpsLMi/Kzg4ZGzP70sdXqtVNU3t20eFCry8SJE0uzY0Svj1p/CI8DeMzeUSvlX0YizG5I+SxotfwSU5f4IB69ekIEPzFLTFuw3HLLld4LjQ0rKU7lHq/LZZddtirty4v4UB7Tt8ZrNHpyxBnNFVZYIeWF48rsiZ4wxf11NY8pcYZ/6aWXzpbHjx8/w/CcmY0bN660XKtjSojpnaPXeYgex7XaR7eXY0oMRS4+ns0dNt2Sz8GVPma09c/ordW+Sh0zqvX4zc4x44s2/tzODuEOrWKnnXYqLdeXXkYKXDzzHmOaoxYF/5cwb7vttln9ixDjlo877riUJzF1Z9Eaa6xRlfuMIS7/ryh8vZdTTz11hg8dxb+3lS+B0YMoDp4hDqDV/pAb06sW3XbbbfWuF9cVP6hvvPHGqRbuvffeUvf9Wn4IDwsttFBWVyeMGjWqwQ/j5R804wNoWxK1gmIK4HDQQQdlY8zbSq+DYnfoOMsa9SfqEn8vnoWN9WvZC66tefzxx7PHJHo4xtTT9913X5sZdtcUMbQi3lvF0GDmHgKV0tgxJS7RnmK7in8r9oBsC2IK9OL+Os4Q1+KYEicE4otifYqhcqjVdPdtqTdo+bGhWMstb8eUeM1FLboQw2GaMyV0Sz8HV/KY0dY/o7dW+yp1zKjW4zc7x4xvWti2Nnu8qFrpZtq94swxMTvL448/3mA1+6gSX0ttZbasmFlnq622KrXlyCOPLLQlMftVVKFvSExDWGz/csstV+fUp7VSy9my/v73vxe+/fbbJk+FXovZlMI222yT3X/MdvbAAw/Mcn1Mpdu7d+/SFKctnYJ1du26666lx+qZZ54p1Nqee+5Zas9pp51W5zqffvppYdVVVy2td99991W1jQ09VzEdbMyWVZzWtrH3ebX3t6+99lo2C1lxBpZvvvlmhuvj9/h78ZgzduzYqravLc+W9dxzzxUWWGCB7DbzzDNPYeTIka3SltZoXzyvjU1VHNN9Fz9PxOXkk0+uWvuaolazZcX6Dc1AVZwKvTgVdby/W7K/np3H7u23385mmInbrLHGGtnMgTOLWaCK2912222r2r5yn3zySemxira2pua2LWboiqnrY/155503m+a5vpmeitPIL7nkkoXvvvuuKu0LH3/8cTY7UF3ic98hhxxS2uYpp5xS9c/BlThmVPIzemscL1qrfZU6ZrRG+yp1zJhSpe9ftZgtq2315yPXLrzwwmyo1aRJk7Jq41FlPXrnxO8xnCLOJhULkB5zzDFVbdvIkSOzekBFxS7LIf4evUDK7bPPPlVpV4y/vP/++7PlGH+5//77z1AEeGYxE088ftUSZ2HiuYqzcXF2LXqXRFfHKKT34osvpj//+c/pscceK7UtnuO2cva/1qIoYnSfjseuf//+WW+hbt26Za+9SO1jJoLi6zAe21p1A73ggguysx2ff/55VqT4qKOOSgMHDszaGrOIxIxkxQKFUbeoFl3ooxD7XXfdlS2vvvrqjdYkqIZTTjklK9oXZ37ifRLDUePMb9QZmDx5cnaGMB7bYs+YmHGsoRnJKiEeqzhzH2em4gxcFNGO9kRPrHjvRm/K6IUUxdCLBZbbyv429nMxw1h0kY5C2nFsiTNqsQ+KIR1RyPC5557L1o31Vlxxxaq2L8y8XrHWRLFAagxNKerTp0+Teii0tG3x2EQX83g/hzPPPDM7C9vQcSVqjcSlKVravihWGe+FqBMSPX6jfkgUY42eeFGwPY4n0aOsWPw+XsPHH398k9rWGu2rpJa2LV5P8ZkqjicxKUA8hsXnLYbV3nLLLdml2Gvn3HPPbdb+ujUeuxiWdfrpp2ez1cRnhKjjFe/b733ve1lvnuixc/nll2frRu/H5hSYb+3nNj6XFmtmtLTXTkvbFr3Z43Uex5X4fBVDZOMzRNQBiYkDos5bHG9ipqrYb4fYNxaLGFe6fcXezzEj6h577JEdV+K5jmPdCy+8kH32K+7/Ysh09ICo9ufgShwzWvMzeiWOF63RvkoeM1qjfZU6ZuzZxr9/tUjVYiTmCNFbYb755isloTNf4gzx66+/XvV2xZmJ+tpU16VamtOmaie/5YlzY5fo2XH//fcX2ppa9txp6mMXPVLirF0tPfroo4VFF1203jZ26NChcNJJJ9WsfZdffnmpLdEDsK345z//WejZs2ejz/Hmm2+e9eKptjgD11C7VltttcJ//vOfNru/jbPS++23X4O33X///Zt99rq12tecbTT17HhL2xa9LZt7XGlOT9qWti/2w029bfTq+PDDD5vcttZoXyXPxFbrsYseIFdeeWWz/6/WfOyOP/747LhR32179epVZw/varUvrLfeetl60dsjeqi2RGu0bfr06YWjjjqqwcctLnPPPXdh6NChVW/fzTff3OBtot2xv66vd099mtOuxt53rX3MaM22VeJ40Rrtq+QxozXaV6ljRmrF57Yheu6Qe3E2KVL86MUTherijH+knZFC/+QnP8lS/5g2kXyIMbfxPEYyHmd34uxRFMuLXh2R2n//+9/PentEATbP64yuu+66bGx89IqJs6pxpizOXEbPp2Lx2jhbGGdhay3OEL388svZDAG33357VgsozmjGTAObbrppdgaxsVmhKulPf/pT9jN6hf3sZz9LbUVMfxnj9+OsUdQEiscwzn7FGaU4s7TOOutk9YF22GGHmtSDiRlV4sxU9MCKadC/+uqrrBh1nEWP/XFMqVuN6aVnV5yVjsc2er/FmeGnn346ex/FDCzx2EadoFoW1qb54mx6HFceeOCB7Ox6fEaI40r0gIueHFFDZP3118/Oqpp0YUZxxvr666/Pjinx2MV7Ot4PUfMrendE77w4wx3TVDe1J1alRI/P2O9FL52oIRNtjd6BceY7/h7HlOgdUCsxXfeTTz6ZLUfvmPqmHq+mOEZET6biVOfR2yZmH433RnxuiM/R0WMm9nu16EEQ9faGDh2a/v3vf2fHvXjfxj46ChlHj7J99903rbfeeqmWHDPaH8eM5usQCc9s3A4AAACANsBsWQAAAAA5JtwBAAAAyDHhDgAAAECOCXcAAAAAcky4AwAAAJBjwh0AAACAHBPuAAAAAOSYcAcAAAAgx4Q7AAAAADkm3AEAAADIMeEOAAAAQI4JdwAAAAByTLgDAAAAkGPCHQAAAIAcE+4AAAAA5JhwBwAAACDHhDsAADl02mmnpQ4dOmQXAGDOJtwBAGiigw46qBSo/Pvf/27Wbe+///7SbY888siKtREAmPMIdwAAmmjQoEGl5euvv75Zt/3Tn/5U53YAAFpKuAMA0EQbbrhhWmGFFbLlv/3tb2nSpElNut3XX3+dbrvttmx5tdVWS/369atoOwGAOYtwBwCgGX7+859nPydOnJjuuOOOJt3m1ltvzQKe8tsDALQW4Q4AQDNEOFMsYtzUoVnFIVkdO3ZMe+21V0XbBwDMeYQ7AADNsPzyy2fDs8J9992XPvzwwwbXf++999K//vWvbHnzzTdPSy65ZLb8xBNPpJNOOiltuummabHFFkudO3dO8803X1p11VXTwQcfnF555ZUWtbNYvDlm1WpI3H+sFz8b8sYbb6Sjjz46rbHGGmn++edP3bp1yx6LffbZJ40ePbpFbQUAWka4AwDQTMWCyNOmTUs33XRTg+vecMMNafr06TPcbvjw4al///7pt7/9bXr44YfTBx98kL799tv05ZdfpjFjxqQrrrgife9730uXXXZZagvOPffcLHS64IIL0ksvvZQNSZs8eXJ666230nXXXZfWXXfddMopp9S6mQAwxxLuAAA002677Za6du06yyxYdSle36NHj7TLLruUQqEFF1ww6/VyzTXXpEcffTQ9++yz6a677kqnn3566tmzZ/ruu+/SYYcd1uwp11vb0KFD05AhQ7LwKQKnyy+/PD3wwANZb50///nPWUhVKBTSGWeckS666KKathUA5lTCHQCAZophSTvssEO2HCHHa6+9Vud6L7zwQnYJEezMM8882fI222yTJkyYkK699tq07777po022iittdZaadttt00nn3xyNgQqgpQITU499dRUKzE07MQTT8yWox3/+c9/0i9+8Yu0xRZbZDN+/fSnP00jR44s1RGKdT/77LOatRcA5lTCHQCA2VAcYtVQ753yv5evH3V3unfv3mB4FD14QoQnn3zySaqF8847L+uxs/baa2fhTrGQdLkoEn3xxRenLl26pK+++irdcsstNWkrAMzJhDsAALNh6623Tosuumi2HMOTopdNuaizE/V2Qu/evdNmm21W77ZimvTx48enl19+OatpE5e55567dP3zzz+fauHOO+/Mfu666651BjtFCyywQFZoOYwaNapq7QMA/o9wBwBgNnTq1CkblhQimIkeNuVihqyYKSv87Gc/y3q4lPv444/Tr3/967TSSiuleeedNy233HJp9dVXz0KSuMQQrfJ1q+3tt99OH330UbZ8wgknlGbfqu9SnDHrf//7X9XbCgBzOuEOAEAFhmbVNyQrPPPMM2nllVdOZ511Vho7duwsvX5mNmnSpFRtjU3xXp9vvvmm1dsCADSsUyPXAwBQj+9///tZL5sXX3wx3XzzzaXaMzHM6tZbb83WicLDMY140dSpU7PZtqKOTgy9Ovzww9OOO+6Y+vbtm82gFbcPb775ZlphhRWy5cbCn0qI2bqKYprzn/zkJ026XbFoNABQPcIdAIAWiF45MVX4559/ntWo+fGPf5xuu+22LOApXl8upjaP4CZcdtllafDgwXVu99NPP21Ru2KoVIRCUfunIcV2zmzhhRcuLUcIFUPGAIC2ybAsAIAWiHo6c801V7Z8/fXXzzAkK0KRPffcc4b1o2hy0e67717vdos1bGZX1PEJDU1NHuFPTLtel+WXXz6btSs89thjLWoLAFBZwh0AgBZYfPHF0w9/+MNs+Z577slmuopiyuFHP/pRWmSRRWZYf9q0aY32moneNn/84x9b1K4o0NxYSHTvvfdmPY7qEoHVwIEDs+X7778/jRkzpkXtAQAqR7gDANBCxaFX3377bdpjjz1K9WpmHpIVVlxxxdLy8OHD69xezE717LPPtqhNAwYMyH4++eSTdfa8iVmtot5PQ6IdEfJE2BTDzSZMmFDvuvE/x5TwDa0DAFSGmjsAAC208847Z8Ogvvzyy9KwqyiOvP3228+y7tZbb5169eqVzUZ10kknZdOox+179uyZDZGKHjvR82fDDTds0XCoAw88MKvpEz2Foh1RFHmjjTbKCjrHdv/whz9kYVSETa+//nqd24hi0eeee246+uij0yuvvJLV3Yntbr755mnRRRdNkydPzto/atSodMstt6T3338/Ky7du3fv2W43ANB8wh0AgBbq1q1b1rPl2muvLf0tZsQqznw182xSI0aMSDvttFMWjlx55ZXZpdymm26aLrnkkhYVMV5ttdXSOeeck375y19mdXcioCm30EILpdtvvz2dfPLJ9YY74aijjsraHD+/+OKLNHTo0OxSl86dO6euXbvOdpsBgNljWBYAQCvYe++9Z/i9riFZ5b13ohbOXnvtlZZYYoms8HLU5omhVFdddVXWc6c1phSPQOcf//hHdn/FadajFs+hhx6annvuubTxxhs3aTsHHHBANsPXb37zm6xHUfQy6tSpU9bGmMJ91113TVdccUV69913U58+fVrcbgCgeToUYpoEAAAAAHJJzx0AAACAHBPuAAAAAOSYcAcAAAAgx4Q7AAAAADkm3AHg/2vHDmgAAAAQBtk/tTm+QQwAAIAwuQMAAAAQJncAAAAAwuQOAAAAQJjcAQAAAAiTOwAAAABhcgcAAAAgTO4AAAAAhMkdAAAAgDC5AwAAABAmdwAAAADC5A4AAABAmNwBAAAACJM7AAAAAGFyBwAAACBM7gAAAACEyR0AAACAMLkDAAAAECZ3AAAAAMLkDgAAAECY3AEAAABY1wGRTnwBQQkDlgAAAABJRU5ErkJggg==",
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {
"image/png": {
"height": 454,
"width": 571
}
},
"output_type": "display_data"
}
],
"source": [
"sns.barplot(x=pd.Series(data).value_counts().sort_index().index,\n",
" y=pd.Series(data).value_counts().sort_index().values,\n",
" color=\"C0\")\n",
"plt.xlabel(\"Value\")\n",
"plt.ylabel(\"Count\")\n",
"plt.title(\"Histogram of observed data\")\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ed039183",
"metadata": {
"id": "ed039183"
},
"outputs": [],
"source": [
"class ZeroInflatedNegativeBinomial(dist.discrete.ZeroInflatedProbs):\n",
" arg_constraints = {\n",
" \"gate\": dist.constraints.unit_interval,\n",
" \"number_of_successes\": dist.constraints.positive,\n",
" \"success_probability\": dist.constraints.unit_interval,\n",
" }\n",
" support = dist.constraints.nonnegative_integer\n",
" pytree_data_fields = (\"number_of_successes\", \"success_probability\")\n",
" def __init__(self, gate, number_of_successes, success_probability, *, validate_args=None):\n",
" _, self.number_of_successes = dist.util.promote_shapes(gate, number_of_successes)\n",
" _, self.success_probability = dist.util.promote_shapes(gate, success_probability)\n",
" negative_binomial_component = dist.NegativeBinomial2(\n",
" concentration=number_of_successes,\n",
" mean=number_of_successes * (1 - success_probability) / success_probability\n",
" )\n",
" super().__init__(negative_binomial_component, gate, validate_args=validate_args)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e35abca8",
"metadata": {
"id": "e35abca8"
},
"outputs": [],
"source": [
"def model(y):\n",
" ap = numpyro.sample(\"ap\", dist.Normal(0, 10))\n",
" p = expit(ap)\n",
" r = numpyro.sample(\"r\", dist.Exponential(1.0))\n",
" aq = numpyro.sample(\"aq\", dist.Normal(0, 10))\n",
" q = expit(aq)\n",
" numpyro.sample(\"y\", ZeroInflatedNegativeBinomial(p, r, q), obs=y)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3c6461a2",
"metadata": {
"id": "3c6461a2",
"outputId": "ecc1967e-ebf9-42dd-c42a-2dfc48279e59"
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"sample: 100%|██████████| 3000/3000 [00:05<00:00, 530.21it/s, 31 steps of size 1.42e-01. acc. prob=0.95]\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
" mean std median 5.0% 95.0% n_eff r_hat\n",
" ap -0.77 0.07 -0.77 -0.88 -0.65 819.51 1.00\n",
" aq -0.40 0.10 -0.40 -0.56 -0.23 654.70 1.00\n",
" r 4.75 0.47 4.74 3.93 5.48 637.92 1.00\n",
"\n",
"Number of divergences: 0\n"
]
}
],
"source": [
"rng_key, rng_key_mcmc = random.split(rng_key)\n",
"nuts_kernel = NUTS(model)\n",
"mcmc = MCMC(nuts_kernel, num_warmup=1000, num_samples=2000)\n",
"mcmc.run(rng_key_mcmc, y=data)\n",
"mcmc.print_summary()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "eb19de18",
"metadata": {
"id": "eb19de18",
"outputId": "8005c5f3-105c-4a8a-88fe-9bea9e3b220a"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"proportion of zero inflation: 0.317\n",
"number of successes: 4.754\n",
"success probability: 0.402\n"
]
}
],
"source": [
"post = mcmc.get_samples()\n",
"print(f\"proportion of zero inflation: {jnp.mean(expit(post['ap'])):0.3f}\") # proportion of zero inflation\n",
"print(f\"number of successes: {jnp.mean(post['r']):0.3f}\") # number of successes\n",
"print(f\"success probability: {jnp.mean(expit(post['aq'])):0.3f}\") # success probability"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9f2c91e9",
"metadata": {
"id": "9f2c91e9"
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "python",
"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.12.10"
},
"colab": {
"provenance": [],
"include_colab_link": true
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment