Skip to content

Instantly share code, notes, and snippets.

@tomonari-masada

tomonari-masada/zero_inflated_poisson.ipynb

Created December 3, 2025 04:10
Show Gist options

  • Select an option

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

Select an option

Learn more about clone URLs
Save tomonari-masada/470541edbf9826ce3aec061c20df791c to your computer and use it in GitHub Desktop.
Download ZIP
zero_inflated_poisson.ipynb
Display the source blob
Display the rendered blob
Raw
zero_inflated_poisson.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"id": "view-in-github",
"colab_type": "text"
},
"source": [
"<a href=\"https://colab.research.google.com/gist/tomonari-masada/470541edbf9826ce3aec061c20df791c/zero_inflated_poisson.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 Poisson distribution"
]
},
{
"cell_type": "code",
"source": [
"!pip install numpyro"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "HOVkLclpZelv",
"outputId": "214db7bd-94a4-4592-c4fe-85b732dd883b"
},
"id": "HOVkLclpZelv",
"execution_count": 1,
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"Collecting numpyro\n",
" Downloading numpyro-0.19.0-py3-none-any.whl.metadata (37 kB)\n",
"Requirement already satisfied: jax>=0.4.25 in /usr/local/lib/python3.12/dist-packages (from numpyro) (0.7.2)\n",
"Requirement already satisfied: jaxlib>=0.4.25 in /usr/local/lib/python3.12/dist-packages (from numpyro) (0.7.2)\n",
"Requirement already satisfied: multipledispatch in /usr/local/lib/python3.12/dist-packages (from numpyro) (1.0.0)\n",
"Requirement already satisfied: numpy in /usr/local/lib/python3.12/dist-packages (from numpyro) (2.0.2)\n",
"Requirement already satisfied: tqdm in /usr/local/lib/python3.12/dist-packages (from numpyro) (4.67.1)\n",
"Requirement already satisfied: ml_dtypes>=0.5.0 in /usr/local/lib/python3.12/dist-packages (from jax>=0.4.25->numpyro) (0.5.4)\n",
"Requirement already satisfied: opt_einsum in /usr/local/lib/python3.12/dist-packages (from jax>=0.4.25->numpyro) (3.4.0)\n",
"Requirement already satisfied: scipy>=1.13 in /usr/local/lib/python3.12/dist-packages (from jax>=0.4.25->numpyro) (1.16.3)\n",
"Downloading numpyro-0.19.0-py3-none-any.whl (370 kB)\n",
"\u001b[2K \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m370.9/370.9 kB\u001b[0m \u001b[31m11.8 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
"\u001b[?25hInstalling collected packages: numpyro\n",
"Successfully installed numpyro-0.19.0\n"
]
}
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "18ee6419",
"metadata": {
"id": "18ee6419"
},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import pandas as pd\n",
"import seaborn as sns\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_key = random.PRNGKey(0)\n",
"\n",
"numpyro.set_platform(\"cpu\")"
]
},
{
"cell_type": "markdown",
"id": "f4134cfc",
"metadata": {
"id": "f4134cfc"
},
"source": [
"## 人工的にデータを合成\n",
"* 確率0.3で頻度が無条件にゼロになる。\n",
"* それ以外の場合、頻度はrateが3.0のポワソン分布に従う。"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "ec802eff",
"metadata": {
"id": "ec802eff"
},
"outputs": [],
"source": [
"size = 200\n",
"\n",
"true_proportion_of_zero_inflation = 0.3\n",
"true_poisson_rate = 3.0\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",
"rng_key, rng_key_poisson = random.split(rng_key)\n",
"poisson_samples = random.poisson(rng_key_poisson, lam=true_poisson_rate, shape=(size,))\n",
"data = jnp.where(is_zero_inflated, 0, poisson_samples)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "5d9ca6a9",
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "5d9ca6a9",
"outputId": "7968605b-eeb2-43ad-c061-cfda56816130"
},
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"Array([0, 0, 4, 1, 0, 0, 0, 2, 0, 2, 1, 4, 0, 0, 5, 3, 0, 1, 0, 0, 0, 2,\n",
" 5, 6, 0, 0, 4, 4, 0, 1, 0, 0, 3, 6, 0, 0, 0, 0, 3, 0, 0, 0, 8, 0,\n",
" 3, 0, 3, 2, 3, 5, 1, 0, 4, 0, 4, 0, 2, 3, 0, 3, 1, 4, 2, 5, 2, 6,\n",
" 6, 0, 3, 7, 0, 5, 5, 2, 8, 1, 1, 5, 4, 2, 0, 2, 0, 3, 0, 4, 0, 1,\n",
" 0, 0, 0, 1, 4, 3, 4, 3, 0, 6, 0, 3, 3, 2, 0, 1, 0, 2, 3, 5, 2, 0,\n",
" 0, 2, 0, 6, 0, 4, 3, 3, 1, 0, 0, 4, 3, 3, 5, 0, 6, 2, 2, 6, 0, 2,\n",
" 0, 0, 1, 0, 0, 5, 1, 1, 3, 2, 1, 1, 0, 0, 0, 6, 2, 3, 0, 1, 0, 3,\n",
" 4, 4, 3, 1, 6, 1, 1, 5, 0, 6, 5, 0, 0, 2, 6, 0, 2, 2, 0, 6, 0, 5,\n",
" 3, 0, 0, 5, 0, 5, 2, 1, 4, 2, 3, 0, 0, 3, 4, 3, 3, 4, 0, 3, 1, 5,\n",
" 4, 4], dtype=int32)"
]
},
"metadata": {},
"execution_count": 4
}
],
"source": [
"data"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "c82450fb",
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 471
},
"id": "c82450fb",
"outputId": "e6cc97b4-327a-4e18-96cd-9ff9e0e7ea6c"
},
"outputs": [
{
"output_type": "display_data",
"data": {
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
],
"image/png": "iVBORw0KGgoAAAANSUhEUgAABGYAAAOMCAYAAAAPMbWZAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAewgAAHsIBbtB1PgAAgsJJREFUeJzs3XmYlWUdP/73ALLjCi4IbiDilhmCGCpqLrkiaqZm4JZlRlJ+1czSrCzX1ChTA7cWs0zxW2iZhqAoIYamiaIkKovizr7J+f3hl/ObkWFmwJl5WF6v65rruuc893M/n2fOmQPnPfdzPxWlUqkUAAAAABpdk6ILAAAAAFhXCWYAAAAACiKYAQAAACiIYAYAAACgIIIZAAAAgIIIZgAAAAAKIpgBAAAAKIhgBgAAAKAgghkAAACAgghmAAAAAAoimAEAAAAoiGAGAAAAoCCCGQAAAICCCGYAAAAACiKYAQAAACiIYAYAAACgIIIZAAAAgIIIZgAAAAAKIpgBAAAAKIhgBoB6sc0226SioiIVFRWZMmVK0eXAJ7ZkyZL8+te/zkEHHZTNNtsszZs3L7/GTznllKLLS+L3bl30gx/8oPyc/+AHP2j0459yyinl4992222NfnyAtZFgBmAdsN9++63yf+QrfwjYb7/9GqQ+WN0sXLgwBx10UM4888w89NBDmTlzZhYvXlx0WQDAWqhZ0QUAwMp45JFHsv/++ydJ+vbtm0ceeaTYglgrXX311VVeW3379k3Xrl3TsmXLJEnv3r0LqgzWPVOmTMm2226bJNl6663NDgPWOoIZAICP+c1vflNu33777RkwYECB1QAAazPBDAD1wl8wWVvMmzcvL774YpKkefPmOfnkkwuuCABYm1ljBgCgkvfee6/c3myzzdKkif8uAQANx/80AAAqqbzIr1AGAGho/rcBQL2o621758yZkxtvvDGHH354ttpqq7Ru3TrrrbdeNthgg3Tv3j1HHnlkfvKTn+S5556rst+yu0MtW/g3SUaNGlU+ZuWvbbbZZoXHX7x4cW699dYcffTR2XrrrdOqVausv/762WGHHXL66afnH//4x0qf+5/+9KcceeSR2XLLLdOiRYt06tQpBx10UG6//fYsWbIkSd1uMVtdn/fffz/XX3999t1332y55ZZp1qxZKioq8v7771fZd+bMmbn11lszcODA7L777tl4442z3nrrZcMNN0z37t1z6qmn5u9//3udzqe62/EuWLAgN910U/bbb79sscUWad68eTp16pQBAwbk+eefX26MOXPm5Je//GX23nvvbLHFFmnZsmW6dOmSs88+O1OnTq1THSurVCrlT3/6U0488cR06dIlbdu2Tdu2bdOlS5ecdNJJufvuu1MqlVa4/7JzXrbIaJK8+uqrK/X6qquGeB0u88ILL2Tw4MHZaaedsv7662f99dfPpz71qXzve9/LG2+8Uef6fvvb3+aYY47Jdtttl7Zt26ZZs2Zp165dunbtmkMOOSQXX3xxxo0bV6fxJk6cmO9+97vp1atX+dbjHTp0yJ577pmLL74406dPr3WMyneXW7Yw84wZM/KTn/wkvXr1yuabb56mTZtmww03TJJ86lOfKve/884761Rnkpx55pnl/c4+++wGP6/KRo4cmZNOOilbb711WrZsmS222CL77LNPbrjhhsybN2+lxloVw4cPT79+/ZZ7L/vNb35Tfi+rq/nz52f48OH55je/mb333rv882nbtm222Wab9O/fP8OGDcuiRYtWOMZtt91Wp9/JZV/Veeqpp/LTn/40RxxxRPm13Lx582y22Wb57Gc/m4suuiivvfbaSp0bQL0rAbDW69u3bylJKUnpkksuWal9L7nkkvK+ffv2XWG/rbfeutzvlVdeqbbP448/Xtpyyy3L/Wr7Wrx4cbV11Pa19dZbV3v8sWPHlrp06VLr/gcddFDprbfeqvVn8/7775cOOuigGsfq06dPacaMGaWBAweWH7v11lurHe/jfR577LFS586dqx33vffeK+93/fXXl5o2bVqnn80BBxxQevvtt2s8r8o/60suuaQ0efLk0m677bbCMVu0aFH629/+Vt5/3LhxNT7P66+/fumJJ56o9ee7MiZNmlTafffdaz3/Hj16lCZPnlztGJ/09VVX9fk6/Pjv3c0331xq0aLFCsfcaKONSvfdd1+NY7744oulHXfcsc4/j5deemmFYy1YsKD01a9+tdbXZ6tWrUpDhgypsa7K72MjR44sDR8+vLTRRhstN9YGG2xQKpVKpSuuuKL82GGHHVbj2JXrrTzm448/3uDnVSqVSosXLy6ddtppNY610047lV544YXlfj/rw+zZs0uHHXZYjcffe++96/xeNnbs2FLbtm3r9PrZZpttSv/+97+rHefWW2+t8+swWf5jTc+ePeu033rrrVe64oor6uVnCbAqLP4LQKN4/fXXc8ghh2T27NlJkvXWWy89e/ZM165d07p168ydOzdTpkzJM888k1mzZi23f69evXL22Wdn2rRpGT58eJKkY8eO6d+//3J9N9lkk+UeGz16dA499NDyX50rKirSq1ev7LTTTlm0aFHGjh2byZMnJ0n+8Y9/pE+fPnnsscfSoUOHas9n4cKF+fznP5+xY8eWH+vYsWP22WeftG3bNpMnT85jjz2WMWPGlGcdrIyXX345gwcPzgcffJB27dpl3333TceOHfPee+9l9OjRVfpOnz49H374YZJku+22y4477pgOHTqkZcuWef/99/Pss8/mv//9b5Lkn//8Zw488MCMHTs2LVq0qLWOWbNm5dBDD82kSZOy/vrrp2/fvtl8883zxhtv5OGHH868efOycOHC9O/fP88++2wWL16cAw88MLNmzUr79u2z7777ZpNNNslrr72Wf/7zn1m8eHFmzZqVo48+Oi+++GI22GCDlfq5VGfixInp27dv3nrrrfJju+66az796U+noqIiEyZMyLPPPpvko7+ef/azn83o0aPTrVu3KuMsmx0xe/bs3HHHHUmSdu3aLXdHpupeX3VV36/Dyu67774MHjw4SbLllltm7733Ttu2bTNp0qSMGTMmS5cuzXvvvZfjjjsuf/nLX3LIIYcsN8bs2bNz4IEH5vXXX0/y0aVcu+++e3bccce0bds28+bNy7Rp0/LMM8/k7bffrrGeuXPn5pBDDsmYMWPKj3Xp0iU9evTIRhttlHfffTdjxozJ9OnTM3/+/AwaNCizZs3Kd7/73VrP9fHHH88PfvCDLF68OJtsskn23XfftG/fPjNnzsyECROSJCeddFIuvPDCLF26NA8++GDeeuutWn+O999/f3mNoa5du2avvfZqlPMaMGBAlVk9G264Yfbff//y784jjzyS559/PocddliOOuqoWn8+K2Px4sU5/PDDq7yvbL755tl3333Trl27vPzyy3nsscfy2GOPpX///nV6L3vvvfcyZ86cJMmmm26anXfeOZ06dUqbNm0yb968vPzyyxk3blyWLFmSKVOmpG/fvvn3v/+drl27Vhlnxx13zNlnn13r7+SKLJsJ06JFi+y8887p2rVrNthgg5RKpcyYMSP/+te/8vbbb2fx4sW54IILkiTnn39+ncYGqFdFJ0MANLzVYcbM4MGDy9v32Wef0rRp06odZ/HixaVHHnmk9KUvfam0ZMmS5baPHDmyTvVU9u6771aZwbH99tuXxo8fv1y/3/72t6VWrVqV+x155JErHPN73/teuV+TJk1KV199denDDz+s0mfy5MmlXr16lZJUmcVQlxkzzZo1KyUpnX322aXZs2dX6bdo0aIqxxo2bFhpyJAhpalTp66w3meeeaa0xx57lMf/0Y9+tMK+lZ/zZXV/5StfKc2aNatKv9dff73UvXv3ct+BAweWevToUaqoqCj94Ac/KC1cuLBK/+eee660+eabl/tfeumlK6yhrhYuXFhlNs+mm25a+sc//rFcv7///e+l9u3bl/t95jOfKS1atKjaMV955ZVyv086O6ayhngdVv69a968ealJkyala665ZrnX4n//+9/SzjvvXO67+eabl959993lxrvuuuuWm6FRnaVLl5bGjRtXOuuss0qvvfZatX0GDBhQHqtbt26lkSNHLtdnyZIlpRtuuKH8OmvatOkKZ6lUfh9r1qxZqaKiovSjH/1ouedxwYIF5fb+++9f3qcuM1eOOeaYWt8r6/u87rjjjiqzN77xjW+U5s2bV6XP9OnTSwcccED5eV7V9/Pq/PCHPyyPV1FRUbrsssuWe+998cUXy79nlY9f04yZ7373u6Vnn312hcd98803S1/+8pfLY33uc59bYd9V/Z0866yzSiNGjFju57nMkiVLSrfeemupTZs25Zkz//vf/+o8PkB9EcwArAMqf6Dp2bNn6eyzz67zV+Wp4J8kmOnRo0d5e02XPtRmVYKZiy++uLzPRhtttMIPkqVSqXTPPfdU+ZA0atSo5fq8++67pZYtW5b7/PSnP13heO+9916Vn01dg5kkpTPOOKNO51dX77//fjkY2WKLLaoNvkql5S8bO/nkk1c45mOPPbbcZQE1fVj87W9/W+634447ftJTKt1yyy1VLkdY0SURpdJHl1gtC7ySlG6//fZq+zVUMFPfr8NSqbTca+vyyy9f4ZgzZsyoEk59//vfX67PscceW95eXcBVV6NHjy6P06VLl1ovyap8ycrnP//5avtUfh9LUvrxj39cax2VXx+9e/euse/7779fJUCt7n2qvs/rww8/rHK54imnnLLCsebNm1f61Kc+Vefftbp4//33S61bty6P94Mf/GCFfWfOnFnaYost6vRetjIOPfTQ8njPP/98tX0a6ndymT/84Q/l8c8///x6Hx+gNoIZgHXAxz/QrOrXJwlmtt9++/L2999/f5XPZWWDmaVLl1aZpfGzn/2s1n0qf1A44YQTltv+i1/8osqHhMpr4VTn9ttvX+lgpmXLltXOaPikzjrrrPIx/vOf/1Tbp3Iw07x589Ibb7xR45hbbbVVuf9mm2223EyZyubNm1f+i3tFRcVys3BW1p577lk+9je/+c1a+1c+/xV9UG+ID4EN8Toslar+3m277ba1vhZ//vOfl/t37NixtHTp0irbK6+Z9PTTT9f9BD/m6KOPLo8zfPjwOu2zbPZVRUVFtesgVX4f69ixY63nWiqVSh988EGV2Ucvv/zyCvv++te/rvW1Ud/ndf/995fHa9WqVa3rPz344IP1GszccMMN5bE6depU4+9uqVQq3XzzzfUezNx1113l8X7+859X26ehg5klS5aU18T5zGc+U+/jA9TGGjMANIrOnTvnpZdeSpLceOON5ev5G9rEiRPLd6Jp2rRpndYmOOOMM/LAAw8kSfnuL5VVfuyLX/ximjWr+Z/T4447Ll/96lezYMGCOtd98MEHZ6ONNqpz/2VmzpyZsWPHZuLEiXnvvfcyd+7cKnciGj9+fLn99NNPZ9ddd61xvH322SebbbZZjX122WWX8loORx55ZJo3b77Cvq1atUqXLl0yceLElEqlTJkypdYaVmT27NlVzue0006rdZ8zzjgjv/rVr5IkTz75ZObOnZs2bdqs0vFXRkO8Dj/upJNOqvW1ePLJJ+db3/pWPvzww0yfPj0vvvhiunfvXt7euXPncvvGG28s/6xWxpIlS8p3lVp//fVzxBFH1Gm//fffPy+88EJKpVLGjBlT41oqxx13XK3nuuz4Rx55ZP74xz8mSX73u9/l4osvrrbv7373u3L75JNPXm57Q5zXyJEjy+3DDjus1vWLDjzwwGy55ZaZNm1anY5dm8rH/+IXv1jj726SnHDCCfnGN75R452UPm7evHkZO3Zsnn322bz11luZPXt2eU2sJFXO5emnn6578SvpP//5TyZMmJApU6Zk1qxZWbhwYZXty+7q9Oyzz2bp0qVp0sTNa4HGI5gBWMdccskl5Vsg18UPfvCDXHrppZ/4uMcff3z++c9/Jkm+853v5B//+Ee+9KUv5aCDDkqnTp0+8fgrsmwh0CTZYYcd6rRwa58+fcrtN954I9OnT0/Hjh3Lj1X+8LDnnnvWOl7r1q2zyy67VAkRatOjR486902S559/PhdccEEeeOCBKh96alLb4q3JR6FLbSoHSDvvvHOt/TfeeONyu7qFnuvqP//5T/lc27Ztm0996lO17vPpT386bdq0ydy5c/Phhx/mmWeeyWc/+9lVrqGuGuJ1+HHVLVT7cRtttFF22GGH8i3OJ0yYUCWYOf7443PLLbck+SiYeeqppzJw4MAccsghyy3MuiL/+c9/Mnfu3CQfLfJ9zjnn1Gm/J598stxetvjwiqzM78fJJ59cazAzderUjBo1KslHNX/xi19crk9DnFfl10Vdnr+Kiorsueeeueeee+p07Nqs7PHbtWuXXXbZJf/+979r7fvuu+/m4osvzh133FFe9L02dXlPWlm33357fvKTn2TSpEl16r948eJ88MEHqxSMA6wqwQwAjeKMM87I3/72t/IdlR5++OE8/PDDSZKtttoq++yzT/bff//069cv7du3r7fjVr5Tz9Zbb12nfTbbbLO0bNmyPMPl7bffrvKBuPKYlWcY1KRTp04rFczU5S48y/z9739Pv379lvsLcG3q8mGpLndNqjxzYWX7L168uNb+K/Lx52HZX7xr0qRJk3Tu3DkvvPBCkob5IFidhngdftxWW21Vp3G32mqrcjBTua4kOeSQQzJo0KAMGTIkyUehwrJgYbPNNsvee++d/fbbL0cfffQKA9Xp06eX2++8805++ctf1qmuypbdGWlFVub34/Of/3zat2+ft99+O5MmTcqTTz6Znj17Vunz+9//vjyzbFn/j2uI86r881+Z56++rOrxawtmXn311ey7777lmXR1VdcApy5KpVJOP/303HrrrSu97+zZswUzQKMyRw+ARtG0adPcc889GTp0aHbaaacq21577bX87ne/yxlnnJGOHTvmjDPOyLvvvlsvx112y9YkK3XJSuW+H/+wUHnM1q1b12m8tm3b1vnYyUeX/NTFW2+9lS9+8YvlUGbrrbfOT3/60zz22GOZPn165s2bl6VLl6b00bpyueSSS8r7Ll26tNbx6xJ2fJL+n0RDPLcNpTFqretrsbYxf/7zn+eee+5Jr169qjz+5ptv5s9//nMGDRqUrbbaKscdd1y1H7w/+OCDOtVRkyVLltS4va6/H8nyM2B++9vfLten8mNf/vKXqx2nIc5rVd5L6vPSu4Y6/kknnVR+bbRr1y7f+ta38re//S3/+9//MmfOnHz44Yfl96TKl1PV5T2prn79619XCWU+//nP5/bbb8+zzz6b9957LwsXLizXUCqVqgSm9VkHQF2YMQNAo6moqMjpp5+e008/PZMmTcqoUaMyZsyYPProo/nf//6X5KMZFMOGDcsjjzySJ554YqX+Ml6dyoHIsssQ6qJy33bt2i035rIPafPmzVvp8erTr3/963Itu+22W0aPHp31119/hf0bK4hoDA3x3DaUxqh1VV6LKxqzf//+6d+/f1577bU88sgjefzxx/Poo4+WZ9qUSqX8+c9/Lm/r1q1bed/KH9w/9alP5ZlnnqlTXQ3p5JNPLs9wueuuu/Kzn/0sTZs2TfLRmiLPPvtsko9mfB155JHVjtEQ51X5dVHEe0lDvJc9/vjjefzxx8vjjx07drkwvrKGek+6+uqry+1LL710hWsLNXQdAHVhxgwAhejWrVu+8pWv5LbbbsvkyZPz4osv5tvf/nb5w9LkyZPrZW2bysFOXafVz5w5s8pCvR+/rKHy91OnTq3TmHXtt7KWXQ6WJN/73vdqDGWSjy4xWFtUfm6nTp1aZZHjFVm6dGmVdT7q87K5mjTE6/Dj6jruypz/VlttlQEDBuTGG2/Mf//737z22mu59NJLy7Mr3nnnnXz729+usk/lxaKXLXhctN69e5fXyHnzzTfLi/gmVWfLHHfccWnZsmW1YzTEea3K66K29XeKPn7l96SBAwfWGMokDfOe9Prrr5cXm99www1z4YUX1th/1qxZtV4+B9CQBDMArBa6deuWa665pkoY83//7/9drt/KXiqz++67l9svvPBCnS6RGjNmTLm9+eabL7eux6c//ely+1//+let482fPz/PPfdcHapdeZXXvajt7kYffvhhlXNb033qU58qB3mzZ88uz3qoyTPPPFP+i3/Tpk2z2267NWiNyzTE6/Djxo4dW+uY77//fnl9nST5zGc+U+s+lXXu3DkXX3xxbr755vJjDz74YJX1jT796U+nRYsWST4Kl15++eWVOkZD+dKXvlRuL7sDU6lUyp133ll+vLq7MS3TEOdV+XVRl+evVCrV6T2noY4/Z86cWt/LVuY9KUlGjx5da5+Vfd+vXEP37t2z3nrr1dj/scceq1OwC9BQBDMArFYq30r2zTffXG575b9m12Xh2B133DGbb755ko+CierWl/i4YcOGldv777//ctv322+/cvuPf/xjreth/PnPf878+fNrPe6qqHxL19ouRRg+fPhqM4OhPrRr1y577LFH+fvbbrut1n0qP7e9evVqlFtlJw3zOvy4O++8s9Y7cv3ud78r99liiy2yww471DpudSr/ni5evLhK0NSqVasccMAB5e9vuOGGVTpGfascugwfPjzz5s3LqFGjyjNAOnfunL59+65w/4Y4r8rP6/33319rYPfPf/6zXmffVT7+XXfdVet76l133VXrIuMr8540ffr03HfffbXWubLv+ytTQ5JVui08QH0SzADQKOp695vK0+Q33XTT5bZXvs3wtGnTah2voqIiZ555Zvn7H/7whzXu93//7//NiBEjyt9/7WtfW67PSSedVP6g8Morr+Taa69d4XgffPBBvv/979da56rabrvtyu3qZhgt89Zbb+Vb3/pWg9VRlK9+9avl9i9/+cv85z//WWHfp556KjfddFP5++qe24bSEK/Dj5s8eXKNr8U333wzP/zhD8vfn3766cvNRFiV39MmTZosd/vvCy64oNweMmRIHnrooTqNmzTc5U9du3ZN7969k3w082P48OHlmTPJRzNqapuZUd/ndfDBB5fv7DZv3rycf/75K9x/wYIFOffcc+t8vLo46aSTypelvf7667niiitW2Pedd96pdZ2WpO7vSR9++GHOPPPMLFq0qNYxN9xww3LY8tZbb9Uazmy77bbl5/K5554rr2FWnbvuuit//etfa60BoCEJZgBoFFtttVW++tWvZtSoUSu848X48eMzaNCg8veHHnrocn223Xbb8geJV199tXwr35oMHjw4W265ZZKPPlx87nOfy9NPP71cvz/84Q858cQTy98feeSR2XfffZfrt/HGG1dZV+M73/lOrrvuuuXOa8qUKfn85z+fKVOmlC+BqG+VFyr96U9/Wu1MjH//+9/p27dvXn/99UabIdJYvvSlL5UvR1q0aFEOOeSQKnd5Weahhx7KoYceWp7d9JnPfKbKc90Y6vt1+HHNmzfPBRdckOuvv3651+LEiRNz0EEHZebMmUk+Wi+luqBur732ykknnZQHHnhghR+YJ02alIEDB5a//9znPpfmzZtX6dO3b99ynyVLluTwww/PT3/60yp3AapswYIFGT58ePr161dlNk59qzxrZtiwYbn77rur3bYi9X1eTZs2zY9+9KMqNQ0ePLjK2kLJR6HOkUcemWeeeWa5n/UnscEGG1QJgy6++OJcccUVy828eumll3LQQQdl+vTptR7/8MMPL4cijzzySP7P//k/y80YfOONN3LsscdmxIgRdXpPatGiRbbffvskH82YGT58eI3927dvXw7hli5dmuOOOy4vvvhilT5Lly7NL3/5y3z5y19O06ZNV7i2EECjKAGw1uvbt28pSSlJ6ZJLLlmpfS+55JLyvn379l1hv6233rrc75VXXllu+7JtSUrt2rUr7bPPPqWTTz659NWvfrV07LHHlnbeeecqfTp06FCaNm1atcc66aSTyv1at25dOuaYY0qDBw8unXvuuaVzzz23dNllly23z6hRo0qtW7cu71dRUVHq3bt36bTTTiudfPLJpa5du1Y5/vbbb1+aOXPmCs93wYIFpV69elXZZ8sttyydcMIJpTPOOKN0wAEHlJo1a1ZKUtprr71KX/rSl8r9br/99mrHHDhwYLnPrbfeusJjf7yObt26Valjxx13LJ1wwgmlU045pbTHHnuUH99tt91K559/fq2vhcrPeV1eLytbd+XX48iRI+t0njV5/vnnSx06dKjyM9htt91KAwcOLA0cOLC02267Vdm26aabll588cUVjvfKK6+U+2699dafuL7K6vt1WPn37rrrriu3O3XqVH4t7rvvvqUmTZqUtzVr1qw0YsSIWsdr1apVac899yydeOKJpa9+9aul448/vsrraVmfZ555ptqxFixYUDr44IOr9G/dunVp//33Lw0cOLB05plnlo4//vhSjx49Si1atCj36dGjR7Xj1cfr5q233iqtt956VWpKUtp9993rPEZ9n1epVCodf/zxVcbbaKONSsccc0zpK1/5Sunzn/98eZxtt922NHjw4FV+P6/OwoULS3369Kly/C222KL8+tlvv/1KTZs2LSUp7bnnnlXef1f0+z5gwIDlxjvqqKNKZ5xxRunAAw8sNW/evPxvwY033linf2O++93vlvutt956pcMOO6z0zW9+s/y+f+6551bp/9BDD1V53a+33nql/fbbr3TaaaeVjj/++NIWW2xR3nbZZZfV+m8YQEMSzACsA1aHYKZt27bLfRha0dduu+1Wmjhx4gqPNWXKlNLmm2++wv1X9GH6iSeeKG233Xa1Hv/AAw+s8cPwMu+9917pgAMOqHGsz372s6UZM2ZU+TBz7733VjveqgQzpVKp9OKLL9Z6Xn369ClNnTq1TqHLmhbMlEof/Qx23333Wp/bz3zmM6WXX365xrEaMpgpler3dfjx37tf/epX5Q+91X1tuOGGpXvuuWeF4+2yyy51/j3ddtttS2PGjKmxviVLlpS+//3vVwmjavpab731SmeffXa1Y9XX6+aII45Y7rjXXHPNSo1Rn+dVKpVKixYtqvJ7VN1X9+7dSxMnTlzp38+6+OCDD0qf//zna30vmz59ep1+3+fOnbtcePXxr06dOpUee+yx0siRI8uP1fRvzPvvv1/q3r17jWN+3K9+9atyQF7dV5MmTUoXX3xxaenSpYIZoFDNAgCN4J133sno0aMzatSoPPnkk3nppZfy5ptvZsGCBWndunU6deqUHj165Nhjj81RRx1VZfHGj9t6663zzDPP5Be/+EUefPDBTJo0KbNnz651Ed7evXtn4sSJ+e1vf5vhw4fn6aefzsyZM7Peeutl8803z957750TTzwxBx98cJ3OacMNN8zDDz+cP/7xj7njjjvy1FNP5d1330379u2z44475stf/nJOOumkrLfeelUW9dxwww3rNH5ddevWLRMmTMgvf/nL3HPPPXnxxRezaNGibL755tl1111z0kkn5fjjjy/fwWht1K1bt4wfPz533313/vznP2fcuHHly3Y23XTT7LnnnjnuuONy7LHHrvQdXupbfb8OK/va176WffbZJzfeeGMeeuih8kKx22yzTY488sgMGjQoW2yxxQr3f/rppzN27NiMHDky48aNy4svvpjp06dn3rx5ad26dTbffPN8+tOfzlFHHZXjjz++1kv0mjZtmh/+8IcZNGhQ7rjjjjz00EN5/vnn8/bbb2fx4sVZf/31s/XWW2fXXXfN/vvvn8MOO6zKLZwbwpe//OUqa4o0bdp0pS9rq+/zWm+99XLbbbdlwIABufnmmzNmzJjMnDkzG220Ubp27Zrjjz8+p512Wtq2bbvK512T9ddfPw888EDuueee3HbbbXnyyServJd96Utfysknn1zr3Y2Wad26dR544IH8/ve/z+23354JEyZk1qxZad++fbbbbrsce+yxOeWUU7LRRhvlkUceqdOYG2ywQZ588snccMMNGTFiRCZOnJj333+/xvVmvva1r6VPnz659tprM3LkyEyfPj2tWrXKlltumQMOOCCnnXZalTtTARSlolRybzgAaGhbbrll+Raub7zxRjbbbLOCKwIAYHVg8V8AaGCPPfZYOZTp3LmzUAYAgDLBDAA0oEWLFlW5+81JJ51UYDUAAKxuBDMAsIrOOuus3HLLLZk9e3a125977rkccMABGT9+fJKkbdu2+frXv96YJQIAsJqzxgwArKL99tsvo0aNSosWLfLpT38622+/fdq2bZtZs2blP//5T/773/9m2T+zFRUVGTZsWE499dSCqwYAYHXirkwA8AktXLgw//rXv/Kvf/2r2u0bbrhhfvnLX7qMCQCA5ZgxAwCr6I033si9996bUaNG5cUXX8zbb7+dd955J0myySabZJdddslBBx2U0047rd5vkQ0AwNpBMAMAAABQEIv/AgAAABREMAMAAABQEMEMAAAAQEEEMwAAAAAFEcwAAAAAFEQwAwAAAFCQZkUXwCezYMGCPPvss0mSDh06pFkzTykAAADUtyVLluStt95Kkuy6665p2bJlvYzrU/wa7tlnn02vXr2KLgMAAADWGePGjUvPnj3rZSyXMgEAAAAUxIyZNVyHDh3K7XHjxmWLLbYosBoAAABYO82YMaN8xUrlz+KflGBmDVd5TZktttginTp1KrAaAAAAWPvV5/quLmUCAAAAKIhgBgAAAKAgghkAAACAgghmAAAAAAoimAEAAAAoiGAGAAAAoCCCGQAAAICCCGYAAAAACiKYAQAAACiIYAYAAACgIIIZAAAAgIIIZgAAAAAKIpgBAAAAKIhgBgAAAKAgghkAAACAgghmAAAAAAoimAEAAAAoiGAGAAAAoCCCGQAAAICCCGYAAAAACiKYAQAAACiIYAYAAACgIIIZAAAAgIIIZgAAAAAKIpgBAAAAKIhgBgAAAKAgghkAAACAgghmAAAAAAoimAEAAAAoiGAGAAAAoCCCGQAAAICCCGYAAAAACtKs6AJofD3Ou6PoEtZqT101oOgSAAAAWEOYMQMAAABQEMEMAAAAQEEEMwAAAAAFEcwAAAAAFEQwAwAAAFAQwQwAAABAQQQzAAAAAAURzAAAAAAURDADAAAAUJB1MpjZb7/9UlFRsVJfjzzyyArHe+CBB9K/f/906tQpLVq0SKdOndK/f/888MADjXdSAAAAwBqnWdEFrAmaNGmS7bfffrnHly5dmjPPPDPDhg2r8vi0adMybdq0DB8+PGeccUZuuummNGmyTmZgAAAAQA3WyWDm1ltvzdy5c2vs8/zzz+eLX/xikuRzn/tcttxyy+X6XHTRReVQZvfdd8/555+fLl26ZPLkybnyyiszYcKEDB06NB06dMhPfvKT+j8RAAAAYI22TgYz2267ba19fvOb35TbAwYMWG77pEmTcvXVVydJ9thjj4wePTqtWrVKkvTs2TNHHXVU+vbtm/Hjx+eqq67Kaaedlq5du9bTGQAAAABrA9fXVGPp0qX53e9+lyRp27ZtjjnmmOX6XHfddVmyZEmSZMiQIeVQZpnWrVtnyJAhSZIlS5bk2muvbeCqAQAAgDWNYKYaDz/8cKZNm5YkOe6449K6desq20ulUu67774kSffu3dO7d+9qx+ndu3d22GGHJMl9992XUqnUgFUDAAAAaxrBTDXuuOOOcru6y5heeeWVTJ8+PUnSt2/fGsdatn3atGmZMmVK/RUJAAAArPEEMx8zZ86c3HvvvUmSrbfeOvvtt99yfZ5//vlyu3v37jWOV3n7xIkT66dIAAAAYK2wTi7+W5M///nP5Ts2nXzyyamoqFiuz9SpU8vtTp061The586dy+3XX399peupfKzqzJgxY6XHBAAAAFYPgpmPqe0ypiSZPXt2ud22bdsax2vTpk25PWfOnJWup3KwAwAAAKxdXMpUydSpU/PII48k+Wjh3m7dulXbb8GCBeV28+bNaxyzRYsW5fb8+fM/eZEAAADAWsOMmUp++9vfZunSpUmSgQMHrrBfy5Yty+1FixbVOObChQvL7Y/fUrsuarv8acaMGenVq9dKjwsAAAAUTzBTyW9+85skH81y+eIXv7jCfu3atSu3a7s8adl6NUntlz1Vp7Y1bAAAAIA1l0uZ/p/x48eX77Z0xBFHZKONNlph38phSW2L81ae8WK9GAAAAKAywcz/U3nR35ouY0qSnXbaqdx+4YUXauxbefuOO+64itUBAAAAayPBTJLFixfnD3/4Q5KkQ4cOOfTQQ2vsv+2226Zjx45JklGjRtXYd/To0UmSLbfcMttss80nLxYAAABYawhmkjzwwAN56623kiQnnXRSmjWreemdioqK9OvXL8lHM2LGjh1bbb+xY8eWZ8z069cvFRUV9Vg1AAAAsKYTzKTqZUwDBgyo0z6DBw9O06ZNkySDBg1a7lbY8+fPz6BBg5IkzZo1y+DBg+unWAAAAGCtsc4HM++9917++te/Jkl22WWXfOYzn6nTft26dct5552X5KOFg/v06ZO77ror48ePz1133ZU+ffpk/PjxSZLzzjsv22+/fcOcAAAAALDGWudvl33XXXdl4cKFSeo+W2aZyy67LDNnzswtt9ySCRMm5IQTTliuz+mnn54f//jH9VIrAAAAsHZZ52fM/OY3v0mSNG3aNF/60pdWat8mTZpk2LBhGTFiRPr165eOHTumefPm6dixY/r165f7778/Q4cOTZMm6/yPGQAAAKjGOj9jZsyYMZ94jMMOOyyHHXZYPVQDAAAArEtM5QAAAAAoiGAGAAAAoCCCGQAAAICCCGYAAAAACiKYAQAAACiIYAYAAACgIIIZAAAAgIIIZgAAAAAKIpgBAAAAKIhgBgAAAKAgghkAAACAgghmAAAAAAoimAEAAAAoiGAGAAAAoCCCGQAAAICCCGYAAAAACiKYAQAAACiIYAYAAACgIIIZAAAAgIIIZgAAAAAKIpgBAAAAKIhgBgAAAKAgghkAAACAgghmAAAAAAoimAEAAAAoiGAGAAAAoCCCGQAAAICCCGYAAAAACiKYAQAAACiIYAYAAACgIIIZAAAAgIIIZgAAAAAKIpgBAAAAKIhgBgAAAKAgghkAAACAgghmAAAAAAoimAEAAAAoiGAGAAAAoCCCGQAAAICCCGYAAAAACiKYAQAAACiIYAYAAACgIIIZAAAAgIIIZgAAAAAKIpgBAAAAKIhgBgAAAKAgghkAAACAgghmAAAAAAoimAEAAAAoiGAGAAAAoCCCGQAAAICCCGYAAAAACiKYAQAAACiIYAYAAACgIIIZAAAAgIIIZgAAAAAKIpgBAAAAKIhgBgAAAKAgghkAAACAgghmAAAAAAoimAEAAAAoiGAGAAAAoCCCGQAAAICCCGYAAAAACiKYAQAAACiIYAYAAACgIIIZAAAAgIIIZgAAAAAKIpgBAAAAKIhgBgAAAKAgghkAAACAgghmAAAAAAoimAEAAAAoiGAGAAAAoCCCGQAAAICCCGYAAAAACiKY+X9ee+21XHLJJdljjz3SoUOHtGzZMp07d84+++yTiy++OM8991yN+z/wwAPp379/OnXqlBYtWqRTp07p379/HnjggUY6AwAAAGBN06zoAlYHQ4YMyYUXXpi5c+dWeXzq1KmZOnVqHnvsscyaNSvXXXfdcvsuXbo0Z555ZoYNG1bl8WnTpmXatGkZPnx4zjjjjNx0001p0kQOBgAAAPz/1vlg5sc//nG+//3vJ0m6deuWr3zlK+nZs2c22GCDvPPOO5kwYULuvffeFYYqF110UTmU2X333XP++eenS5cumTx5cq688spMmDAhQ4cOTYcOHfKTn/yk0c4LAAAAWP1VlEqlUtFFFOXhhx/OgQcemCQZMGBAhg4dmvXWW6/avosWLUrz5s2rPDZp0qTsvPPOWbJkSfbYY4+MHj06rVq1Km+fN29e+vbtm/Hjx6dZs2aZOHFiunbtWq/nMHXq1HTu3DlJ8vrrr6dTp0617tPjvDvqtQaqeuqqAUWXAAAAQD1blc/fdbHOXluzdOnSnHXWWUmS3XbbLcOGDVthKJNkuVAmSa677rosWbIkyUeXQ1UOZZKkdevWGTJkSJJkyZIlufbaa+urfAAAAGAtsM4GMw8++GBeeumlJMkFF1yQZs1W7qquUqmU++67L0nSvXv39O7du9p+vXv3zg477JAkue+++7IOT1ACAAAAPmadDWb+9Kc/JUkqKipyxBFHlB9/991389JLL+Xdd9+tcf9XXnkl06dPT5L07du3xr7Ltk+bNi1Tpkz5BFUDAAAAa5N1NpgZO3ZskmSbbbZJu3bt8vvf/z677rprNtlkk3Tr1i2bbLJJdthhh1x99dVZuHDhcvs///zz5Xb37t1rPFbl7RMnTqynMwAAAADWdOvkXZmWLl2aF154IUnSvn37nHPOOfn5z3++XL9JkyblvPPOy7333psRI0Zkww03LG+bOnVquV3bgj/LFgdKPlogaGVUPk51ZsyYsVLjAQAAAKuPdTKY+eCDD7J06dIkybPPPpsnn3wyW2yxRa666qocdthhadmyZZ588slccMEFGTt2bB5//PGcdtppueeee8pjzJ49u9xu27Ztjcdr06ZNuT1nzpyVqrVyqAMAAACsXdbJS5nmzp1bbi9YsCCtW7fOyJEj86UvfSkbbbRRWrVqlX333Tf//Oc/s9tuuyVJ7r333vzrX/+qst8y1d2xqbIWLVqU2/Pnz6+v0wAAAADWcOvkjJmWLVtW+f6MM84o3zmpslatWuWyyy4rLw581113Zc8991xujEWLFtV4vMpr1Hz8ltq1qe3SpxkzZqRXr14rNSYAAACwelgng5l27dpV+f7ggw9eYd/Pfe5zadasWZYsWZInn3yy2jFquzyp8gyd2i57+rja1q8BAAAA1lzr5KVMLVq0SIcOHcrf17SOS8uWLdO+ffskyVtvvVV+vHJgUtsCvZVnvVgzBgAAAFhmnQxmkmTnnXcutz/88MMa+y7b3qzZ/z/BaKeddiq3l93haUUqb99xxx1Xqk4AAABg7bXOBjP77rtvuf2///1vhf1mzZqVt99+O0my5ZZblh/fdttt07FjxyTJqFGjajzW6NGjy/tvs802q1oyAAAAsJZZZ4OZY489tty+9957V9jv3nvvTalUSpLss88+5ccrKirSr1+/JB/NiBk7dmy1+48dO7Y8Y6Zfv36pqKj4xLUDAAAAa4d1Npj51Kc+lUMPPTRJcuedd+bhhx9ers8bb7yR733ve0k+uiX2qaeeWmX74MGD07Rp0yTJoEGDlrsV9vz58zNo0KAkH10GNXjw4Po+DQAAAGANts4GM0ly3XXXZcMNN8zSpUtzxBFH5MILL8yjjz6a8ePH54YbbkjPnj3LC/v+6Ec/qnIpU5J069Yt5513XpJk/Pjx6dOnT+66666MHz8+d911V/r06ZPx48cnSc4777xsv/32jXuCAAAAwGqtorTsOp111GOPPZbjjjsub775ZrXbKyoqctFFF+VHP/pRtduXLl2ar3zlK7nllltWeIzTTz89N998c5o0qf8cbOrUqeU7Pb3++ut1ur12j/PuqPc6+P89ddWAoksAAACgnq3K5++6WKdnzCTJ3nvvnf/+97+55JJLsttuu2X99ddPy5Yts+222+bUU0/NU089tcJQJkmaNGmSYcOGZcSIEenXr186duyY5s2bp2PHjunXr1/uv//+DB06tEFCGQAAAGDNts7PmFnTmTGz+jFjBgAAYO1jxgwAAADAWkYwAwAAAFAQwQwAAABAQQQzAAAAAAURzAAAAAAURDADAAAAUBDBDAAAAEBBBDMAAAAABRHMAAAAABREMAMAAABQEMEMAAAAQEEEMwAAAAAFEcwAAAAAFEQwAwAAAFAQwQwAAABAQQQzAAAAAAURzAAAAAAURDADAAAAUBDBDAAAAEBBBDMAAAAABRHMAAAAABREMAMAAABQEMEMAAAAQEEEMwAAAAAFEcwAAAAAFEQwAwAAAFAQwQwAAABAQQQzAAAAAAURzAAAAAAURDADAAAAUBDBDAAAAEBBBDMAAAAABRHMAAAAABREMAMAAABQEMEMAAAAQEEEMwAAAAAFEcwAAAAAFEQwAwAAAFAQwQwAAABAQQQzAAAAAAURzAAAAAAURDADAAAAUBDBDAAAAEBBBDMAAAAABRHMAAAAABREMAMAAABQEMEMAAAAQEEEMwAAAAAFEcwAAAAAFEQwAwAAAFAQwQwAAABAQQQzAAAAAAURzAAAAAAURDADAAAAUBDBDAAAAEBBBDMAAAAABRHMAAAAABREMAMAAABQEMEMAAAAQEEEMwAAAAAFEcwAAAAAFEQwAwAAAFAQwQwAAABAQQQzAAAAAAURzAAAAAAURDADAAAAUBDBDAAAAEBBBDMAAAAABRHMAAAAABREMAMAAABQEMEMAAAAQEEEMwAAAAAFEcwAAAAAFEQwAwAAAFAQwQwAAABAQQQzAAAAAAVZZ4OZioqKOn3tt99+tY71wAMPpH///unUqVNatGiRTp06pX///nnggQca/kQAAACANVazogtYky1dujRnnnlmhg0bVuXxadOmZdq0aRk+fHjOOOOM3HTTTWnSZJ3NwAAAAIAVWOeDmbPOOitf//rXV7i9TZs2K9x20UUXlUOZ3XffPeeff366dOmSyZMn58orr8yECRMydOjQdOjQIT/5yU/qvXYAAABgzbbOBzObbrppdtlll5Xeb9KkSbn66quTJHvssUdGjx6dVq1aJUl69uyZo446Kn379s348eNz1VVX5bTTTkvXrl3rtXYAAABgzeb6mlV03XXXZcmSJUmSIUOGlEOZZVq3bp0hQ4YkSZYsWZJrr7220WsEAAAAVm+CmVVQKpVy3333JUm6d++e3r17V9uvd+/e2WGHHZIk9913X0qlUqPVCAAAAKz+BDOr4JVXXsn06dOTJH379q2x77Lt06ZNy5QpUxq6NAAAAGANss4HM3/605+y0047pXXr1mnXrl223377DBw4MCNHjlzhPs8//3y53b179xrHr7x94sSJn7xgAAAAYK2xzi/+WzlkSZKXX345L7/8cu64444cffTRue2227LBBhtU6TN16tRyu1OnTjWO37lz53L79ddfX+n6Kh+rOjNmzFjpMQEAAIDVwzobzLRu3TpHHXVUPve5z6V79+5p27Zt3nrrrYwaNSo33nhj3nnnnQwfPjz9+vXLP/7xj6y33nrlfWfPnl1ut23btsbjVL7d9pw5c1a6zsrBDgAAALB2WWeDmWnTpmXDDTdc7vGDDjoogwYNyqGHHpoJEyZk1KhR+dWvfpVvfvOb5T4LFiwot5s3b17jcVq0aFFuz58//5MXDgAAAKw11tlgprpQZpnNNtssd999d7p3757FixdnyJAhVYKZli1bltuLFi2q8TgLFy4stz9+S+26qO3ypxkzZqRXr14rPS4AAABQvHU2mKnNdtttl4MOOij3339/Xn755UyfPj0dO3ZMkrRr167cr7bLk+bOnVtu13bZU3VqW8MGAAAAWHOt83dlqslOO+1Ubk+bNq3crhyW1LY4b+UZL9aLAQAAACoTzNSgoqKi2scrBzYvvPBCjWNU3r7jjjvWT2EAAADAWkEwU4PKt9JedhlTkmy77bbl70eNGlXjGKNHj06SbLnlltlmm23qv0gAAABgjSWYWYFXXnkl//jHP5IkXbp0yZZbblneVlFRkX79+iX5aEbM2LFjqx1j7Nix5Rkz/fr1W+EMHAAAAGDdtE4GM3/5y1+yZMmSFW5/8803c+yxx5bvuPT1r399uT6DBw9O06ZNkySDBg1a7lbY8+fPz6BBg5IkzZo1y+DBg+upegAAAGBtsU7elWnQoEFZvHhxjj322Oy1117ZZptt0qpVq7z99tt55JFHctNNN+Xtt99Okuy99945++yzlxujW7duOe+883L55Zdn/Pjx6dOnTy644IJ06dIlkydPzhVXXJEJEyYkSc4777xsv/32jXqOAAAAwOpvnQxmkmT69OkZMmRIhgwZssI+xx57bIYOHZoWLVpUu/2yyy7LzJkzc8stt2TChAk54YQTlutz+umn58c//nG91Q0AAACsPdbJYOb222/PqFGj8sQTT+R///tf3n777cyaNStt27ZN586d89nPfjYDBw7MXnvtVeM4TZo0ybBhw3Lsscfm5ptvzpNPPpm333477du3T8+ePfPVr341hx56aCOdFQAAALCmWSeDmb59+6Zv3771Nt5hhx2Www47rN7GAwAAANYN6+TivwAAAACrA8EMAAAAQEEEMwAAAAAFEcwAAAAAFEQwAwAAAFAQwQwAAABAQQQzAAAAAAURzAAAAAAURDADAAAAUBDBDAAAAEBBBDMAAAAABRHMAAAAABREMAMAAABQEMEMAAAAQEEEMwAAAAAFEcwAAAAAFEQwAwAAAFAQwQwAAABAQQQzAAAAAAURzAAAAAAURDADAAAAUBDBDAAAAEBBBDMAAAAABRHMAAAAABREMAMAAABQEMEMAAAAQEEEMwAAAAAFEcwAAAAAFEQwAwAAAFAQwQwAAABAQQQzAAAAAAURzAAAAAAURDADAAAAUBDBDAAAAEBBBDMAAAAABRHMAAAAABREMAMAAABQEMEMAAAAQEEEMwAAAAAFEcwAAAAAFEQwAwAAAFAQwQwAAABAQQQzAAAAAAURzAAAAAAURDADAAAAUBDBDAAAAEBBBDMAAAAABRHMAAAAABREMAMAAABQEMEMAAAAQEEEMwAAAAAFEcwAAAAAFEQwAwAAAFAQwQwAAABAQQQzAAAAAAURzAAAAAAURDADAAAAUBDBDAAAAEBBBDMAAAAABRHMAAAAABREMAMAAABQEMEMAAAAQEEEMwAAAAAFEcwAAAAAFEQwAwAAAFAQwQwAAABAQQQzAAAAAAURzAAAAAAURDADAAAAUJBGDWYOOOCAfO5zn8urr75a532mT59e3g8AAABgbdKsMQ/2yCOPpKKiInPnzq3zPvPnzy/vBwAAALA2cSkTAAAAQEFW+2Bm2eyali1bFlwJAAAAQP1a7YOZBx54IEnSqVOngisBAAAAqF8NusbMaaedVu3j3/ve97LhhhvWuO/ChQszefLkPPnkk6moqEjfvn0boEIAAACA4jRoMHPbbbctt2hvqVTKfffdV6f9S6VSkmTjjTfOhRdeWO/1AQAAABSpQYOZrbbaqkow8+qrr6aioiJbbLFF1ltvvRXuV1FRkZYtW2aLLbbIZz/72Zx11lnp2LFjQ5ZaxQUXXJArr7yy/P3IkSOz33771bjPAw88kJtvvjlPPvlk3nrrrXTo0CE9e/bMmWeemUMPPbSBKwYAAADWRA0azEyZMqXK902afLSkzYMPPpiddtqpIQ+9yp5++un87Gc/q3P/pUuX5swzz8ywYcOqPD5t2rRMmzYtw4cPzxlnnJGbbrqpfP4AAAAASSMv/rvvvvtm3333TZs2bRrzsHW2LGRZsmRJNt100zrtc9FFF5VDmd133z133nlnxo0blzvvvDO77757kmTo0KH53ve+12B1AwAAAGumRg1mHnnkkYwcOTJbb711Yx62zn7+85/nySefTPfu3XP66afX2n/SpEm5+uqrkyR77LFHxowZkxNOOCE9e/bMCSeckMceeyx77LFHkuSqq67Kyy+/3KD1AwAAAGsW19b8P6+99lq+//3vJ0luvPHGNG/evNZ9rrvuuixZsiRJMmTIkLRq1arK9tatW2fIkCFJkiVLluTaa6+t56oBAACANVmDrjFTF7Nmzcrs2bPz4Ycf1tp3q622arA6zj777MyZMycDBw5M3759M3LkyBr7V767VPfu3dO7d+9q+/Xu3Ts77LBDXnzxxdx33335xS9+sdydqgAAAIB1UyHBzD/+8Y/ccMMNeeyxx/Luu+/WaZ+Kiory7JT69sc//jF//etfs/HGG5cvTarNK6+8kunTpydJ+vbtW2Pfvn375sUXX8y0adMyZcqUbLvttp+4ZgAAAGDN1+jBzDe/+c388pe/TPLRrJOivf/++znnnHOSJFdccUXat29fp/2ef/75crt79+419q28feLEiSsVzEydOrXG7TNmzKjzWAAAAMDqpVGDmd///vf5xS9+kSRp2bJljj766PTo0SMbb7xxYbeSPv/88/PGG2+kT58+dVrwd5nKgUmnTp1q7Nu5c+dy+/XXX1+p+irvCwAAAKxdGjWYuemmm5J8FDb885//TJcuXRrz8Mt59NFHM3To0DRr1iw33njjSq39Mnv27HK7bdu2NfatfHvwOXPmrHyhAAAAwFqpUYOZ//znP6moqMgll1xSeCizaNGinHnmmSmVSvnWt76VXXbZZaX2X7BgQbld2x2cWrRoUW7Pnz9/pY5T2wybGTNmpFevXis1JgAAALB6aNRgZvHixUmS3XffvTEPW62f/OQneeGFF7LVVlvlkksuWen9W7ZsWW4vWrSoxr4LFy4stz9+S+3a1HaZFAAAALDmatSFXbbZZpskxV/O88ILL+SnP/1pkmTIkCFVLjWqq3bt2pXbtZ3P3Llzy+3aLnsCAAAA1h2NOmPmmGOOyWWXXZaHH344++yzT2Meuoprr702ixYtynbbbZd58+blD3/4w3J9nnvuuXL7n//8Z954440kyZFHHpk2bdpUmclS252TKl+OZDFfAAAAYJlGDWbOPffc/OY3v8l1112XE044odbbTDeUZZcW/e9//8uJJ55Ya/8f/ehH5fYrr7ySNm3aZKeddio/9sILL9S4f+XtO+6448qWCwAAAKylGvVSpg022CB///vfs9lmm+Wzn/1sbrjhhrz33nuNWUK92XbbbdOxY8ckyahRo2rsO3r06CTJlltuWb6cCwAAAKBRg5ntttsun//85/PBBx/k/fffz6BBg9KhQ4dsvvnm2W677Wr8qs+7ON12220plUo1flVeEHjkyJHlx5cFKxUVFenXr1+Sj2bEjB07ttpjjR07tjxjpl+/fit1S24AAABg7daolzJNmTKlyvfLwo6ZM2fWuu/qGGgMHjw4N998cz788MMMGjQoo0ePrnLXpfnz52fQoEFJkmbNmmXw4MEFVQoAAACsjho1mBk4cGBjHq7BdevWLeedd14uv/zyjB8/Pn369MkFF1yQLl26ZPLkybniiisyYcKEJMl5552X7bffvuCKAQAAgNVJowYzt956a2MerlFcdtllmTlzZm655ZZMmDAhJ5xwwnJ9Tj/99Pz4xz8uoDoAAABgddaoa8ysjZo0aZJhw4ZlxIgR6devXzp27JjmzZunY8eO6devX+6///4MHTo0TZr4UQMAAABVVZRKpVLRRbDqpk6dms6dOydJXn/99XTq1KnWfXqcd0dDl7VOe+qqAUWXAAAAQD1blc/fdWEaBwAAAEBBGnWNmTvu+GQzNQYMMBMBAAAAWHs0ajBzyimnrPJtrysqKgQzAAAAwFqlUYOZJLGkDQAAAMBHGjWYeeWVV2rtM3fu3EyaNCm///3vc/fdd6dPnz65+eab07p160aoEAAAAKDxNGows/XWW9ep30477ZSjjz46f/zjH3PSSSdl0KBB+cc//tHA1QEAAAA0rtX6rkzHH398Bg4cmJEjR+amm24quhwAAACAerVaBzPJR+FMqVTKbbfdVnQpAAAAAPVqtQ9mNttssyTJiy++WHAlAAAAAPVrtQ9mXnvttSTJ4sWLC64EAAAAoH6t1sHM4sWLc+WVVyZJunbtWnA1AAAAAPWrUe/KtGz2S02WLl2a9957L+PHj88vfvGLPPfcc6moqMgJJ5zQCBUCAAAANJ5GDWa23Xbbld6nVCplr732yre+9a0GqAgAAACgOI16KVOpVFqpr4022igXXnhhHnroobRo0aIxSwUAAABocI06Y+bWW2+ttU+TJk3Srl27bLvtttlll13StGnTRqgMAAAAoPE1ajAzcODAxjwcAAAAwGpttb4rEwAAAMDaTDADAAAAUJBGvZTp45566qk89NBDee655/Luu+8mSTbeeOPssssuOfDAA9OjR48iywMAAABoUIUEM88++2zOPPPMjBs3boV9vvvd72bPPffMTTfdlF133bURqwMAAABoHI1+KdNDDz2UXr16Zdy4ceXbYjdr1iybbbZZNttsszRr1qz8+NixY9OrV688/PDDjV0mAAAAQINr1GDm7bffzhe+8IUsXLgwFRUVOeOMM/Kvf/0rc+fOzfTp0zN9+vTMmzcv48aNy1e+8pU0bdo0CxcuzBe+8IW88847jVkqAAAAQINr1GDm+uuvzwcffJDmzZtnxIgRufnmm9OzZ880a/b/X1HVtGnT7LHHHrnpppsyYsSIrLfeevnggw9y/fXXN2apAAAAAA2uUYOZESNGpKKiIt/4xjdyyCGH1Nr/4IMPzqBBg1IqlTJixIhGqBAAAACg8TRqMPPKK68kSY466qg677Os7//+978GqQkAAACgKI0azCxYsCBJ0qZNmzrvs6zvwoULG6QmAAAAgKI0ajCz+eabJ0kmTJhQ532W9d1ss80apCYAAACAojRqMLPPPvukVCrl8ssvz6xZs2rtP3v27FxxxRWpqKjIPvvs0wgVAgAAADSeRg1mvvrVryb5aK2ZfffdN+PHj19h3/Hjx6dv376ZPHlylX0BAAAA1hbNau9Sf/r06ZOvf/3rueGGG/Lss89mzz33zM4775w999wzm266aSoqKvLmm2/mX//6V/773/+W9/v617+ePn36NGapAAAAAA2uUYOZJBkyZEhat26dn/3sZ1m6dGmee+65KiFMkpRKpSRJkyZN8n/+z//J5Zdf3thlAgAAADS4Rr2UKUkqKipy5ZVX5umnn85ZZ52V7bffPqVSqcrX9ttvn7POOitPP/10eY0ZAAAAgLVNo8+YWWaXXXbJL3/5yyTJokWL8t577yVJNtpoozRv3ryosgAAAAAaTWHBTGXNmzd3O2wAAABgndOglzI98MAD+cxnPpPPfOYz+f3vf79S+/7+978v7/vQQw81UIUAAAAAxWmwYKZUKuVb3/pWnnnmmXTo0CEnnXTSSu1/4oknpn379nn66adz7rnnNlCVAAAAAMVpsGDmn//8ZyZNmpQmTZrk2muvXen9Kyoqct1116Vp06Z57rnnMmrUqAaoEgAAAKA4DRbM/PnPf06SHHTQQdlpp51WaYyddtophxxySJLk7rvvrrfaAAAAAFYHDRbMjBs3LhUVFTnyyCM/0ThHHHFESqVSxo4dW0+VAQAAAKweGiyYefXVV5MkO+ywwycap1u3bkmSKVOmfNKSAAAAAFYrDRbMfPDBB0mSjTfe+BONs2z/WbNmfeKaAAAAAFYnDRbMrL/++kmS999//xONs2z/du3afcKKAAAAAFYvDRbMdOjQIUny/PPPf6JxJk6cmCTZdNNNP3FNAAAAAKuTBgtmevXqlVKplL/85S+faJz77rsvFRUV6dmzZz1VBgAAALB6aLBg5tBDD02SPPjgg3nsscdWaYzRo0fnwQcfrDIeAAAAwNqiwYKZY489Nttss01KpVK+8IUv5KWXXlqp/SdNmpTjjz8+FRUV2WabbXLcccc1UKUAAAAAxWiwYGa99dbL1VdfnSSZOXNmevTokeuvvz5z586tcb85c+bkuuuuyx577JGZM2cmSa655po0a9asoUoFAAAAKESDph3HHHNMLr300lxyySWZO3duvv3tb+f73/9+9tlnn/To0SObbrpp2rRpk7lz5+bNN9/Mv//97zz66KOZO3duSqVSkuTSSy/N0Ucf3ZBlAgAAABSiwaehfP/730+nTp0yaNCgzJs3L3PmzMnf/va3/O1vf6u2/7JApnXr1vnFL36RU045paFLBAAAAChEg13KVNmpp56aSZMm5dvf/nbat2+fUqm0wq/27dvn3HPPzaRJk4QyAAAAwFqt0RZu6dixY66++upcffXV+e9//5tnnnkm77zzTmbPnp127dplk002yW677Zadd965sUoCAAAAKFQhK+ruvPPOAhgAAABgndcolzIBAAAAsDzBDAAAAEBBBDMAAAAABRHMAAAAABREMAMAAABQEMEMAAAAQEEEMwAAAAAFEcwAAAAAFEQwAwAAAFAQwQwAAABAQQQzAAAAAAURzAAAAAAURDADAAAAUBDBDAAAAEBBBDMAAAAABRHMAAAAABREMAMAAABQEMEMAAAAQEEEMwAAAAAFEcwAAAAAFEQwAwAAAFAQwQwAAABAQQQzAAAAAAURzAAAAAAUZJ0MZmbNmpU//OEPOffcc9O3b9907do1G2ywQZo3b55NN900++23X6688sq88847dRrv8ccfz8knn5ytt946LVu2zOabb55DDjkkd955ZwOfCQAAALAmqyiVSqWii2hsDz30UA466KBa+7Vv3z6//e1vc8ghh6ywzw9+8IP86Ec/ytKlS6vdfvjhh+fuu+9Oy5YtV7nemkydOjWdO3dOkrz++uvp1KlTrfv0OO+OBqmFjzx11YCiSwAAAKCercrn77pYJ2fMJEnnzp0zYMCAXH/99bnnnnvyxBNPZMyYMbnrrrvyhS98IU2bNs3bb7+do446Ks8880y1Y9x000259NJLs3Tp0nTp0iXDhg3LuHHjMnz48Oy///5JkhEjRuS0005rzFMDAAAA1hDr5IyZDz/8ME2bNq2xz/Dhw9O/f/8kSf/+/XPPPfdU2f7uu+9mu+22ywcffJCtttoqTz31VNq3b1/lGP37989f/vKXJMnIkSOz33771e+JxIyZ1ZEZMwAAAGsfM2bqUW2hTJIcffTR2WGHHZIkjz766HLbhw4dmg8++CBJcsUVV1QJZZYd44Ybbigf66qrrvqkZQMAAABrmXUymKmrdu3aJUkWLFiw3Lbhw4cnSdZff/0cc8wx1e7fqVOnHHjggUmShx9+OLNnz26YQgEAAIA1kmBmBV588cU8/fTTSZLu3btX2bZo0aKMGzcuSbLXXnulefPmKxynb9++SZKFCxdm/PjxDVMsAAAAsEZqVnQBq5N58+Zl2rRp+ctf/pIrr7wyS5YsSZIMHjy4Sr9Jkyblww8/TLJ8aPNxlbdPnDixvChwXU2dOrXG7TNmzFip8QAAAIDVxzofzNx222059dRTV7j9O9/5Tk466aQqj1UOS2pb7GfZwkDJR4sDrazK+wMAAABrl3U+mFmRT3/607n55pvTs2fP5bZVXiumbdu2NY7Tpk2bcnvOnDn1VyAAAACwxlvng5mjjz46e+yxR5Jk/vz5mTx5cv74xz/m3nvvzYknnpjrrrsuRxxxRJV9Ki8GXNP6MknSokWLcnv+/PkrXV9ts2xmzJiRXr16rfS4AAAAQPHW+WBmww03zIYbblj+vmfPnjnhhBPym9/8JgMHDky/fv0ybNiwnHLKKeU+LVu2LLcXLVpU4/gLFy4st1u1arXS9dXXfdEBAACA1Y+7Mq3Al7/85XzhC1/I0qVL841vfCPvvvtueduy22gntV+eNHfu3HK7tsueAAAAgHWLYKYG/fr1S/JRuPK3v/2t/HjlWSy13TWp8qVIFvIFAAAAKhPM1KBDhw7l9quvvlpud+vWLU2bNk2SvPDCCzWOUXn7jjvuWM8VAgAAAGsywUwNpk2bVm5XvgypefPm5QV3n3jiiRrXmRk1alSSjxYBXrbIMAAAAEAimKnRn/70p3J71113rbLt6KOPTpLMmjUr99xzT7X7T506NQ899FCS5HOf+1yVtWkAAAAA1slg5rbbbqtyy+vqXHvttbn//vuTJNtuu2322WefKtvPOOOMbLDBBkmS73znO3nnnXeqbP/www/z9a9/PR9++GGS5Lzzzquv8gEAAIC1xDp5u+wf/OAHOffcc3Psscdm7733TpcuXdK2bdvMnj07zz77bH73u99lzJgxST66bOnmm28urymzzMYbb5wrrrgiX/va1/Lqq69mzz33zEUXXZRdd90106dPz3XXXZeRI0cmSU488cTst99+jX2aAAAAwGpunQxmkuTdd9/Nr3/96/z6179eYZ9OnTrllltuyYEHHljt9q9+9auZPn16fvSjH2Xy5Mk57bTTlutz2GGH5ZZbbqm3ugFoHD3Ou6PoEtZqT101oOgSAABWC+tkMPP3v/89I0aMyJgxY/Lyyy/nzTffzDvvvJNWrVpl0003zac//ekcccQROf7449O6desax7r00ktzyCGH5Je//GUeffTRvPnmm9lwww2z22675dRTT82JJ57YSGcFAAAArGnWyWBmhx12yA477JBvf/vb9TLeZz/72Xz2s5+tl7EAAACAdcc6ufgvAAAAwOpAMAMAAABQEMEMAAAAQEEEMwAAAAAFEcwAAAAAFEQwAwAAAFAQwQwAAABAQQQzAAAAAAURzAAAAAAURDADAAAAUBDBDAAAAEBBBDMAAAAABRHMAAAAABREMAMAAABQEMEMAAAAQEEEMwAAAAAFEcwAAAAAFEQwAwAAAFAQwQwAAABAQQQzAAAAAAURzAAAAAAURDADAAAAUBDBDAAAAEBBBDMAAAAABRHMAAAAABREMAMAAABQEMEMAAAAQEEEMwAAAAAFEcwAAAAAFEQwAwAAAFAQwQwAAABAQQQzAAAAAAURzAAAAAAURDADAAAAUBDBDAAAAEBBBDMAAAAABRHMAAAAABREMAMAAABQEMEMAAAAQEEEMwAAAAAFEcwAAAAAFEQwAwAAAFAQwQwAAABAQQQzAAAAAAURzAAAAAAURDADAAAAUJBmRRcAsLbrcd4dRZewVnvqqgFFlwAAAKvMjBkAAACAgghmAAAAAAoimAEAAAAoiGAGAAAAoCCCGQAAAICCCGYAAAAACiKYAQAAACiIYAYAAACgIIIZAAAAgII0K7oAoG56nHdH0SWs1Z66akDRJQAAAOsgM2YAAAAACiKYAQAAACiIYAYAAACgIIIZAAAAgIIIZgAAAAAKIpgBAAAAKIhgBgAAAKAgghkAAACAgghmAAAAAAoimAEAAAAoiGAGAAAAoCCCGQAAAICCCGYAAAAACiKYAQAAACiIYAYAAACgIIIZAAAAgIIIZgAAAAAKIpgBAAAAKMg6G8yMHz8+P/zhD3PwwQenU6dOadGiRdq2bZtu3brl1FNPzWOPPbZS4z3wwAPp379/eaxOnTqlf//+eeCBBxroDAAAAIA1XbOiCyjCvvvum0cffXS5xxctWpSXXnopL730Um677bYMGDAgv/71r9O8efMVjrV06dKceeaZGTZsWJXHp02blmnTpmX48OE544wzctNNN6VJk3U2BwMAAACqsU4mBdOnT0+SdOzYMeecc07uvvvujBs3Lk888UR+9rOfZcstt0yS3HHHHTnllFNqHOuiiy4qhzK777577rzzzowbNy533nlndt999yTJ0KFD873vfa/hTggAAABYI62TM2a6d++en/zkJzn22GPTtGnTKtt69+6dL3/5y+nTp08mTZqUO++8M1/72tey7777LjfOpEmTcvXVVydJ9thjj4wePTqtWrVKkvTs2TNHHXVU+vbtm/Hjx+eqq67Kaaedlq5duzb8CQIAAABrhHVyxsxf//rXHH/88cuFMsu0b98+11xzTfn7u+++u9p+1113XZYsWZIkGTJkSDmUWaZ169YZMmRIkmTJkiW59tpr66N8AAAAYC2xTgYzdbH//vuX25MnT15ue6lUyn333Zfkoxk4vXv3rnac3r17Z4cddkiS3HfffSmVSg1QLQAAALAmEsyswMKFC8vt6mbWvPLKK+W1avr27VvjWMu2T5s2LVOmTKm/IgEAAIA1mmBmBUaNGlVu77jjjsttf/7558vt7t271zhW5e0TJ06sh+oAAACAtcE6ufhvbZYuXZrLL7+8/P3xxx+/XJ+pU6eW2506dapxvM6dO5fbr7/++krVUvk41ZkxY8ZKjQcAAACsPgQz1bj22mszbty4JMkxxxyTHj16LNdn9uzZ5Xbbtm1rHK9Nmzbl9pw5c1aqlsqhDgAAALB2Ecx8zKhRo/Kd73wnSbLpppvmV7/6VbX9FixYUG43b968xjFbtGhRbs+fP78eqgQAqtPjvDuKLmGt9tRVA4ouAQDWOoKZSv773/+mf//+WbJkSVq2bJk//elP2XTTTavt27Jly3J70aJFNY5beSHhj99Suza1Xfo0Y8aM9OrVa6XGBAAAAFYPgpn/55VXXsnBBx+c9957L02bNs0f/vCH7Lvvvivs365du3K7tsuT5s6dW27XdtnTx9W2fg0AAACw5nJXpiTTp0/PgQcemOnTp6eioiK33HJL+vXrV+M+lQOT2hborTzrxZoxAAAAwDLrfDDz9ttv56CDDsr//ve/JMmQIUMyYEDt10/vtNNO5fYLL7xQY9/K26u79TYAAACwblqng5kPPvgghxxySJ5//vkkyeWXX56zzz67Tvtuu+226dixY5KPFgyuyejRo5MkW265ZbbZZptVLxgAAABYq6yzwcy8efNy+OGH59///neS5KKLLsoFF1xQ5/0rKirKlzu98MILGTt2bLX9xo4dW54x069fv1RUVHzCygEAAIC1xToZzCxatCj9+/fPmDFjkiTnnHNOfvzjH6/0OIMHD07Tpk2TJIMGDVruVtjz58/PoEGDkiTNmjXL4MGDP1nhAAAAwFplnbwr04knnpgHH3wwSXLAAQfk9NNPz3PPPbfC/s2bN0+3bt2We7xbt24577zzcvnll2f8+PHp06dPLrjggnTp0iWTJ0/OFVdckQkTJiRJzjvvvGy//fYNc0IAAADAGmmdDGbuueeecvuf//xnPvWpT9XYf+utt86UKVOq3XbZZZdl5syZueWWWzJhwoSccMIJy/U5/fTTV2lGDgAAALB2WycvZapPTZo0ybBhwzJixIj069cvHTt2TPPmzdOxY8f069cv999/f4YOHZomTfyoAQAAgKrWyRkzpVKp3sc87LDDcthhh9X7uAAAAMDayzQOAAAAgIIIZgAAAAAKIpgBAAAAKIhgBgAAAKAgghkAAACAgghmAAAAAAoimAEAAAAoiGAGAAAAoCCCGQAAAICCCGYAAAAACiKYAQAAACiIYAYAAACgIIIZAAAAgIIIZgAAAAAKIpgBAAAAKIhgBgAAAKAgghkAAACAgghmAAAAAAoimAEAAAAoiGAGAAAAoCCCGQAAAICCCGYAAAAACiKYAQAAACiIYAYAAACgIIIZAAAAgIIIZgAAAAAKIpgBAAAAKEizogsAAGDd1uO8O4ouYa321FUDii4BgBqYMQMAAABQEMEMAAAAQEEEMwAAAAAFEcwAAAAAFEQwAwAAAFAQwQwAAABAQQQzAAAAAAURzAAAAAAURDADAAAAUBDBDAAAAEBBBDMAAAAABRHMAAAAABREMAMAAABQEMEMAAAAQEEEMwAAAAAFEcwAAAAAFEQwAwAAAFAQwQwAAABAQQQzAAAAAAURzAAAAAAURDADAAAAUBDBDAAAAEBBBDMAAAAABRHMAAAAABREMAMAAABQEMEMAAAAQEEEMwAAAAAFEcwAAAAAFKRZ0QUAAABrnh7n3VF0CWu9p64aUHQJQCMwYwYAAACgIIIZAAAAgIIIZgAAAAAKIpgBAAAAKIhgBgAAAKAgghkAAACAgghmAAAAAAoimAEAAAAoiGAGAAAAoCCCGQAAAICCCGYAAAAACiKYAQAAACiIYAYAAACgIIIZAAAAgIIIZgAAAAAKIpgBAAAAKIhgBgAAAKAgghkAAACAgghmAAAAAAqyzgYzM2fOzF//+tdcfPHFOfTQQ9O+fftUVFSkoqIip5xyykqP98ADD6R///7p1KlTWrRokU6dOqV///554IEH6r94AAAAYK3QrOgCirLZZpvVyzhLly7NmWeemWHDhlV5fNq0aZk2bVqGDx+eM844IzfddFOaNFlnczAAAACgGpKCJFtttVUOPvjgVdr3oosuKocyu+++e+68886MGzcud955Z3bfffckydChQ/O9732v3uoFAAAA1g7r7IyZiy++OD179kzPnj2z2WabZcqUKdl2221XaoxJkybl6quvTpLsscceGT16dFq1apUk6dmzZ4466qj07ds348ePz1VXXZXTTjstXbt2rfdzAQAAANZM6+yMmUsvvTRHHHHEJ7qk6brrrsuSJUuSJEOGDCmHMsu0bt06Q4YMSZIsWbIk11577aoXDAAAAKx11tlg5pMqlUq57777kiTdu3dP7969q+3Xu3fv7LDDDkmS++67L6VSqdFqBAAAAFZvgplV9Morr2T69OlJkr59+9bYd9n2adOmZcqUKQ1dGgAAALCGEMysoueff77c7t69e419K2+fOHFig9UEAAAArFnW2cV/P6mpU6eW2506daqxb+fOncvt119/fZWPU50ZM2as1HgAAADA6kMws4pmz55dbrdt27bGvm3atCm358yZs1LHqRzqAAAAAGsXlzKtogULFpTbzZs3r7FvixYtyu358+c3WE0AAADAmsWMmVXUsmXLcnvRokU19l24cGG5/fFbatemtkufZsyYkV69eq3UmAAAAMDqQTCzitq1a1du13Z50ty5c8vt2i57+rja1q8BAAAA1lwuZVpFlQOT2hborTzrxZoxAAAAwDKCmVW00047ldsvvPBCjX0rb99xxx0brCYAAABgzSKYWUXbbrttOnbsmCQZNWpUjX1Hjx6dJNlyyy2zzTbbNHRpAAAAwBpCMLOKKioq0q9fvyQfzYgZO3Zstf3Gjh1bnjHTr1+/VFRUNFqNAAAAwOpNMPMJDB48OE2bNk2SDBo0aLlbYc+fPz+DBg1KkjRr1iyDBw9u7BIBAACA1dg6e1emxx57LC+//HL5+7fffrvcfvnll3PbbbdV6X/KKacsN0a3bt1y3nnn5fLLL8/48ePTp0+fXHDBBenSpUsmT56cK664IhMmTEiSnHfeedl+++0b5FwAAACANdM6G8wMHTo0t99+e7XbxowZkzFjxlR5rLpgJkkuu+yyzJw5M7fccksmTJiQE044Ybk+p59+en784x9/4poBAACAtYtLmT6hJk2aZNiwYRkxYkT69euXjh07pnnz5unYsWP69euX+++/P0OHDk2TJn7UAAAAQFXr7IyZ2267bbnLlT6Jww47LIcddli9jQcAAACs/UzjAAAAACiIYAYAAACgIIIZAAAAgIIIZgAAAAAKIpgBAAAAKIhgBgAAAKAgghkAAACAgghmAAAAAAoimAEAAAAoiGAGAAAAoCCCGQAAAICCCGYAAAAACiKYAQAAACiIYAYAAACgIIIZAAAAgIIIZgAAAAAKIpgBAAAAKIhgBgAAAKAgghkAAACAgghmAAAAAAoimAEAAAAoiGAGAAAAoCCCGQAAAICCCGYAAAAACiKYAQAAACiIYAYAAACgIIIZAAAAgIIIZgAAAAAKIpgBAAAAKIhgBgAAAKAgghkAAACAgghmAAAAAAoimAEAAAAoiGAGAAAAoCCCGQAAAICCCGYAAAAACiKYAQAAACiIYAYAAACgIIIZAAAAgIIIZgAAAAAKIpgBAAAAKIhgBgAAAKAgghkAAACAgghmAAAAAAoimAEAAAAoiGAGAAAAoCCCGQAAAICCCGYAAAAACiKYAQAAACiIYAYAAACgIIIZAAAAgIIIZgAAAAAKIpgBAAAAKIhgBgAAAKAgghkAAACAgghmAAAAAAoimAEAAAAoiGAGAAAAoCCCGQAAAICCCGYAAAAACiKYAQAAACiIYAYAAACgIIIZAAAAgIIIZgAAAAAKIpgBAAAAKIhgBgAAAKAgzYouAAAAAKhZj/PuKLqEtdpTVw0o7NhmzAAAAAAURDADAAAAUBDBDAAAAEBBBDMAAAAABRHMAAAAABREMAMAAABQEMEMAAAAQEEEMwAAAAAFEcwAAAAAFKRZ0QWsTV599dX8/Oc/z4gRI/L666+nRYsW6dKlS44//vicffbZad26ddElAgAA67ge591RdAlrtaeuGlB0CaxhBDP15C9/+UtOPvnkzJo1q/zYvHnzMn78+IwfPz5Dhw7NiBEj0rVr1wKrBAAAAFYnLmWqBxMmTMgXv/jFzJo1K23bts1ll12Wxx9/PA8//HC+8pWvJEkmTZqUww8/PLNnzy64WgAAAGB1YcZMPTjnnHMyf/78NGvWLA8++GD22muv8rYDDjgg22+/fc4///xMmjQp11xzTX7wgx8UVywAAACw2jBj5hMaN25cHn300STJ6aefXiWUWebcc8/NjjvumCS5/vrrs3jx4katEQAAAFg9CWY+oeHDh5fbp556arV9mjRpkgEDPloA6v3338/IkSMbozQAAABgNSeY+YQee+yxJEmbNm3So0ePFfbr27dvuT1mzJgGrwsAAABY/Vlj5hOaOHFikqRr165p1mzFP87u3bsvt09dTJ06tcbtr7/+erk9Y8aMOo25aPa7dT4+K6+252xVed4aVkM9b4nnrqH5nVszed7WTJ63NZPnbc3luVszed7WTHV53ip/5l6yZEm9HbuiVCqV6m20dcyCBQvSqlWrJMnhhx+ev/71rzX2b9u2bebOnZvevXvniSeeqNMxKioqPnGdAAAAQP0ZN25cevbsWS9juZTpE6h86+u2bdvW2r9NmzZJkjlz5jRYTQAAAMCaw6VMn8CCBQvK7ebNm9fav0WLFkmS+fPn1/kYlS9VWlENL7zwQjbbbLN06NChxsup1jQzZsxIr169knyURm6xxRYFV0RdeN7WXJ67NZPnbc3keVszed7WXJ67NZPnbc20Nj9vS5YsyVtvvZUk2XXXXett3LXnU3wBWrZsWW4vWrSo1v4LFy5MkvLlT3XRqVOnWvt07dq1zuOtqbbYYos6/SxYvXje1lyeuzWT523N5HlbM3ne1lyeuzWT523NtDY+b9tss029j+lSpk+gXbt25XZdLk+aO3dukrpd9gQAAACs/QQzn0DLli2zySabJKl9Bef33nuvHMx07ty5wWsDAAAAVn+CmU9op512SpK8/PLLNd4u64UXXii3d9xxxwavCwAAAFj9CWY+ob333jvJR5cpPfXUUyvsN2rUqHK7T58+DV4XAAAAsPoTzHxCRx99dLl96623Vttn6dKlueOOO5IkG264Yfbff//GKA0AAABYzQlmPqFevXpln332SZIMGzYsTzzxxHJ9rrnmmkycODFJcs4552S99dZr1BoBAACA1ZPbZdeD66+/Pn369Mn8+fNz8MEH57vf/W7233//zJ8/P3/4wx9y8803J0m6deuWc889t+BqAQAAgNVFRalUKhVdxNrgL3/5S04++eTMmjWr2u3dunXLiBEj0rVr10auDAAAAFhdCWbq0auvvprrr78+I0aMyNSpU9O8efN07do1X/jCF/KNb3wjrVu3LrpEAAAAYDUimAEAAAAoiMV/AQAAAAoimAEAAAAoiGAGAAAAoCCCGQAAAICCCGYAAAAACiKYAQAAACiIYAYAAACgIIIZAAAAgIIIZlgtvfrqqzn33HPTvXv3tGnTJhtvvHF69uyZq666KvPmzSu6PD5m5syZ+etf/5qLL744hx56aNq3b5+KiopUVFTklFNOKbo8VmD8+PH54Q9/mIMPPjidOnVKixYt0rZt23Tr1i2nnnpqHnvssaJL5GNmzZqVP/zhDzn33HPTt2/fdO3aNRtssEGaN2+eTTfdNPvtt1+uvPLKvPPOO0WXykq44IILyu+ZFRUVeeSRR4ouiUoqPzc1fe23335Fl0oNXnvttVxyySXZY4890qFDh7Rs2TKdO3fOPvvsk4svvjjPPfdc0SWu8/bbb786/755v1w9LVq0KEOHDs0hhxySLbbYovx/yx122CGnnnpqHn/88aJLXG1VlEqlUtFFQGV/+ctfcvLJJ2fWrFnVbu/WrVtGjBiRrl27NnJlrEhFRcUKtw0cODC33XZb4xVDney777559NFHa+03YMCA/PrXv07z5s0boSpq89BDD+Wggw6qtV/79u3z29/+NoccckgjVMUn8fTTT6dnz55ZsmRJ+bGRI0f6kL8aqenfuMr69u3rQ+JqasiQIbnwwgszd+7cFfY555xzct111zVeUSxnv/32y6hRo+rcv0mTJnnttdey5ZZbNmBV1NWrr76aww8/PP/9739r7Ddo0KBcf/31dX5vXVc0K7oAqGzChAn54he/mPnz56dt27a58MILs//++2f+/Pn5wx/+kF//+teZNGlSDj/88IwfPz7t2rUrumQ+Zquttkr37t3z4IMPFl0KNZg+fXqSpGPHjvnCF76QffbZJ1tttVU+/PDDPPHEE7nmmmsybdq03HHHHVm8eHF+//vfF1wxy3Tu3Dn7779/evTokc6dO2eLLbbI0qVLM3Xq1Nx9992555578vbbb+eoo47KuHHjsttuuxVdMiuwdOnSnHnmmVmyZEk23XTTzJw5s+iSqMFZZ52Vr3/96yvc3qZNm0ashrr68Y9/nO9///tJPvrj3le+8pX07NkzG2ywQd55551MmDAh9957b5o0cSFB0W699dYaw7Mkef755/PFL34xSfK5z31OKLOaWLx4cZVQ5lOf+lS+/e1vZ4cddsjs2bPz2GOP5ZprrsncuXMzZMiQdOzYMd/5zncKrno1U4LVyD777FNKUmrWrFnp8ccfX277lVdeWUpSSlK65JJLGr9AqnXxxReX/vKXv5TeeOONUqlUKr3yyivl52ngwIHFFke1Dj/88NJdd91VWrJkSbXb33rrrVK3bt3Kz+OoUaMauUKqs6Lnq7J77723/Lz179+/EapiVV177bWlJKXu3buXLrzwwvLzNnLkyKJLoxL/71hzPfTQQ+Xnb8CAAaVFixatsO/ChQsbsTJW1fnnn19+Tn/zm98UXQ7/z5/+9Kfy87LXXntV+/+V8ePHl9Zbb71SktKGG25YWrx4cQGVrr5Ew6w2xo0bV7604vTTT89ee+21XJ9zzz03O+64Y5Lk+uuvz+LFixu1Rqp36aWX5ogjjshmm21WdCnU0V//+tccf/zxadq0abXb27dvn2uuuab8/d13391YpVGDFT1flR199NHZYYcdkqROl6tRjNdee638V/wbb7zR5YJQz5YuXZqzzjorSbLbbrtl2LBhWW+99VbY3+/g6m/p0qX53e9+lyRp27ZtjjnmmIIrYpnKa8dceOGF1f5/pUePHjniiCOSJO+//34mTpzYaPWtCQQzrDaGDx9ebp966qnV9mnSpEkGDBiQ5KNf6JEjRzZGabBO2n///cvtyZMnF1gJK2vZZZ4LFiwouBJW5Oyzz86cOXMycODA9O3bt+hyYK3z4IMP5qWXXkry0QLbzZpZwWFN9/DDD2fatGlJkuOOOy6tW7cuuCKWWbRoUbm93XbbrbBfly5dqt0HwQyrkWV3gGnTpk169Oixwn6V/wM7ZsyYBq8L1lULFy4st+syU4PVw4svvpinn346SdK9e/dii6Faf/zjH/PXv/41G2+8ca6++uqiy4G10p/+9KckHy3evOyv9Eny7rvv5qWXXsq7775bVGmsojvuuKPcXvaHWlYPy2bqJsn//ve/FfZb9oe+ioqKbL/99g1e15pEMMNqY9l0tq5du9b4V43KHzRMgYOGU/nOCMsuIWT1NG/evLz00kv52c9+lr59+5bv8DN48OBiC2M577//fs4555wkyRVXXJH27dsXXBF19ac//Sk77bRTWrdunXbt2mX77bfPwIEDzd5dTY0dOzZJss0226Rdu3b5/e9/n1133TWbbLJJunXrlk022SQ77LBDrr766ip/iGD1NGfOnNx7771Jkq233tqd61YzJ554YtZff/0kH/3b9uGHHy7XZ8KECRkxYkSS5KSTTir35yPm9LFaWLBgQd5+++0kSadOnWrsu9FGG6VNmzaZO3duXn/99cYoD9Y5S5cuzeWXX17+/vjjjy+wGqpz2223rfCyzyT5zne+k5NOOqkRK6Iuzj///Lzxxhvp06dPTj/99KLLYSU8//zzVb5/+eWX8/LLL+eOO+7I0Ucfndtuuy0bbLBBQdVR2dKlS/PCCy8k+WjNtHPOOSc///nPl+s3adKknHfeebn33nszYsSIbLjhho1cKXX15z//uXzHppNPPtmtllcz7du3z29+85uceOKJGTNmTHr27JnBgwenW7dumTNnTsaMGZNrrrkmixYtymc+85kq6xjyETNmWC3Mnj273G7btm2t/ZfdknLOnDkNVhOsy6699tqMGzcuSXLMMcfUeHkhq5dPf/rTGTduXH7605/6j+tq5tFHH83QoUPTrFmz3HjjjZ6fNUTr1q1zwgkn5Ne//v/au/OgqOvHj+Mv5BC8LY9MK0GlBsGpsYvAQoPIvMP7CGvMLi1pxsowpdPKI8rGMxUVi9Ikx9Qy0Rhw6EBNAbUipAZJCkUxDgHZ3x/89vOFXBBF+Gzu8zHDzAc+7915MQ22+9r3sVJJSUk6cOCAdu7cqcjISF177bWSqvbJGzZsGIcS2IkzZ86osrJSkpSWlqYPPvhAXbp0UWxsrE6dOqXi4mIlJibq7rvvllS1celjjz1mZmRcBMuY7N/QoUO1b98+TZkyRT/99JPCw8Pl7++vkJAQRUVFqUWLFoqOjlZSUhIHhtjAjBnYheobVNZnV/zmzZtLkkpKShotE+CoEhMT9dJLL0mSOnXqpKVLl5qcCLYMHz5ct99+u6Sqfwt/++03ffbZZ4qPj9e4ceMUHR1dY18FmKusrExTp06VxWJRRESEfH19zY6Eejp+/LjNmRQhISGaPn26Bg4cqAMHDigxMVFLly7Vs88+2/QhUYN1ZoVU9RqzRYsW2rNnT419MO69917t3r1b/v7+OnjwoOLj4/X999/rrrvuMiMy6pCTk6Nvv/1WknT33XfL29vb3ECwqaysTOvWrdOWLVtksVguuJ+Xl6fY2Fh5enpq6NChJiS0b8yYgV1wd3c3ruuzQ7d1LbCHh0ejZQIcUUZGhkaMGKGKigq5u7tr48aN6tSpk9mxYEO7du3k6+srX19f3XHHHRo7dqw2b96sdevWKSsrS8OGDVNMTIzZMfH/3nrrLR09elQ33nij5s6da3YcXIK6lrd07txZmzZtMo5hXrx4cROlQl2qv66UpClTptQoZaw8PDz05ptvGt9/+umnjZ4Nly42NtaYARUeHm5yGthSVFSk4OBgzZs3T6dOndILL7ygI0eO6Ny5czpz5ox27typwMBApaamavjw4Vq0aJHZke0OxQzsgvVoV6l+y5Osn4TUZ9kTgPo5duyYHnjgARUUFMjZ2VlxcXG69957zY6FSzRp0iSNGjVKlZWVmjZtGieP2IGjR49q3rx5kqreuFuX4+Lq4OXlpZCQEElV+87k5uaanAjVX1dK0gMPPFDr2Pvvv984dOLHH39s1Fy4POvXr5dUNWN+zJgxJqeBLVFRUUpKSpIkrVq1Su+8845uueUWubm5qU2bNgoJCdGePXvUv39/WSwWzZw5UwcPHjQ5tX2hmIFdcHd3N9Zp5+Tk1Dm2oKDAKGZuuOGGRs8GOILc3FwFBwcrNzdXTk5OWr16tYYNG2Z2LFwm63+7oqIiffXVVyanwXvvvaeysjJ5eXmpuLhYcXFxF3ylp6cb43fv3m38vPqSDNgvHx8f4/r48eMmJoFU9Qa+Y8eOxvd1vV50d3c3Tkf7+++/Gz0bLk1qaqqx8fbgwYPVvn17kxPh3ywWi1avXi1J8vb2rnVWk4uLi15//XVJVRt0M6u3JvaYgd3w8fFRUlKSMjMzVVFRUeuR2dZd9iWO8AWuhPz8fIWEhCgrK0tS1Sf6bKz331b9Dcnvv/9uYhJI/1t+m5WVpXHjxl10vPWFq1Q1k40ZNvaPjZztT+/evY19SWwd3Vud9X5trz1hnuqb/rKMyT7l5eUZs3Nvu+22OsdWP0yi+ns6MGMGdiQwMFBS1Se8+/btq3VcYmKicR0QENDouYCr2ZkzZxQaGmp8GvX222/rmWeeMTkVGqr6J/Ys+QQaX/WjtK+//noTk8Cq+lJc6wcPthQWFio/P1+S1LVr10bPhforLy9XXFycpKoPHAYOHGhyIthSvdCsqKioc2z1k+soQmuimIHdGD58uHG9Zs0am2MqKyuN5rxdu3bq379/U0QDrkrFxcUaNGiQ9u/fL0mKjIzUiy++aHIqXAkbN240rv38/ExMAkmKiYmRxWKp86v6hsB79uwxft69e3fzgqNejh07pm+++UaS1KNHD97c24mwsDDjOj4+vtZx8fHxxgky/fr1a/RcqL8dO3YYy8vGjx/PG3k7dc0116hNmzaSpJSUlDrLmeofsHt6ejZ6tv8SihnYjTvvvNP4H+KqVauUkpJywZiFCxfqyJEjkqTnnnvOOAUBwKUpKyvTiBEjtHfvXklVf09vvPGGyalwMTExMSotLa1zzHvvvaft27dLqnrRwxsN4PJt3bq1zjcZeXl5CgsLM06UfPrpp5sqGi6iT58+xgyLTz75RAkJCReMOXHihGbPni1JcnNz06OPPtqkGVG36suYWGJtv5o1a6ZBgwZJqtqzsPpJZ9UVFBTU+ABw8ODBTZLvv8LJYuuQccAkBw4cUEBAgEpKStSqVSu9/PLL6t+/v0pKShQXF6cVK1ZIqtpYKjU19YJd92GO5ORkZWZmGt/n5+dr5syZkqqWm02ZMqXG+MmTJzdlPNgQFhamzZs3S5IGDBig6OjoOvdIcHNzk7e3d1PFQy26d++us2fPKiwsTIGBgerRo4datWqls2fPKi0tTRs2bDDKNjc3N23btk3BwcEmp0Z9REVF6dVXX5VUNWMmKCjI3ECQVPU3V15errCwMPn7+6t79+7y8PBQfn6+vv32Wy1fvtxYBhMYGKhdu3apefPmJqeG1S+//KK77rpLp0+flru7u2bMmKGHHnpIHh4e+uGHHzRv3jzj0Il33nlHL7zwgsmJYVVQUKAuXbro3Llz8vX1VVpamtmRUIejR4+qb9++Ki4uliQNGTJE4eHh8vLyUmlpqb777jtFR0frjz/+kFR1GtquXbvMjGx3KGZgd7Zu3aqJEyeqsLDQ5n1vb29t27ZNPXv2bOJkqM3kyZO1du3aeo/nnx3zXepGlTfddJOys7MbJwzqrXv37vXazLdbt25avXq1cYQv7B/FjH2q799cWFiYPvroI7Vr167xQ+GSJCcna+TIkcrLy7N538nJSZGRkTU23Yb5li1bpqeeekqS9O677xof+MF+7dq1S+PGjTPK6toMGDBAmzZt4oStf2GhHuzOkCFDdOjQIb3//vvatm2bcnJy5Obmpp49e2rUqFGaNm2aWrRoYXZMAGhyX3/9tbZt26a9e/cqMzNTeXl5OnnypDw8PNSpUyfdeuutGjx4sEaPHs2/k8AVsHbtWiUmJiolJUVZWVnKz89XYWGhWrVqpRtuuEH33HOPwsPD5e/vb3ZU1CIwMFAZGRlavHixvvjiCx07dkxlZWXq0qWLgoKCNH369IueJIOmt379ekmSs7OzJkyYYHIa1EdwcLCOHj2qVatWaceOHcrIyNDp06fl4uKi6667TnfccYfGjx+voUOHcpKdDcyYAQAAAAAAMAmb/wIAAAAAAJiEYgYAAAAAAMAkFDMAAAAAAAAmoZgBAAAAAAAwCcUMAAAAAACASShmAAAAAAAATEIxAwAAAAAAYBKKGQAAAAAAAJNQzAAAAAAAAJiEYgYAAAAAAMAkFDMAAAAAAAAmoZgBAAAAAAAwCcUMAAAAAACASShmAAAAAAAATEIxAwAAAAAAYBKKGQAAAAAAAJNQzAAAADShqKgoOTk5ycnJyewoAADADlDMAACAq94TTzxhlCG7d+++pMfu3LnTeOxzzz3XSAkBAICjopgBAABXvUceecS4jo2NvaTHrl+/3ubzAAAAXAkUMwAA4KoXEBCgHj16SJI+//xzlZSU1OtxRUVFio+PlyT17t1bffv2bbSMAADAMVHMAAAAhzBp0iRJUmFhobZs2VKvx2zevFlFRUU1Hg8AAHAlUcwAAACHMGnSJGPD3fouZ7IuY2rWrJkmTpzYaNkAAIDjopgBAAAOwcvLSwEBAZKkr7/+Wn/99Ved43Nzc5WQkCBJGjBggLp27SpJ+u677zR79mwFBQXpuuuuk5ubm9q0aSMfHx899dRTOnz4cINyWjcajoqKqnNcUFCQnJycFBQUVOe4zMxMRUREyM/PT23btpWHh4e8vLw0efJkpaamNigrAABoOIoZAADgMKyb91ZUVCguLq7OsR9//LEqKytrPC4mJkb+/v568803lZiYqLy8PJWXl+vs2bM6cuSIli1bpj59+mjJkiWN+4vU04IFC+Tj46Po6Gilp6ersLBQpaWlOnbsmNauXas777xTc+bMMTsmAAAOjWIGAAA4jNGjR8vd3V1SzdOWbLHeb9WqlR5++GFJVYVO+/btNXnyZK1evVpJSUnav3+/vvzyS7322mvq0KGDzp8/r2nTpl3ysdxX2vz58zVz5kyVl5erT58+Wrp0qXbt2qXU1FRt2LBB/v7+slgsev311/XBBx+YmhUAAEdGMQMAABxG27ZtNXToUElSamqqfv75Z5vjDh06pEOHDkmSHn74YbVs2VKSNHDgQOXk5GjNmjV69NFHFRgYqNtuu02DBg3SK6+8oszMTPXp00cWi0Vz585tml/KhsOHDysyMlKSNHfuXP3000968skndf/996tv374aP368kpOTjX1zIiMjVVBQYFpeAAAcGcUMAABwKNZlSVLts2aq/7z6+K5du6pFixa1Pnfbtm312muvSZKSk5N18uTJhsa9LAsXLlR5ebluv/12zZ0719j0uLpmzZpp8eLFat68uf755x9t2rTJhKQAAIBiBgAAOJTQ0FB17txZkrRhwwZZLJYa9ysrK/Xxxx9Lkrp166b+/fvX+lxFRUXKzs5WRkaG0tPTlZ6eLldXV+P+wYMHG+E3uLitW7dKksLCwmyWMlbt2rWTn5+fJCklJaVJsgEAgJooZgAAgENxcXHR+PHjJUnZ2dlKTk6ucT8hIUG5ubmSpAkTJqhZs5ovl/Lz8/Xyyy/r5ptvVuvWreXp6SlfX1/5+fnJz89PgwYNqjG2qf3+++/6+++/JUmzZs0yTnmq7ct6MtOJEyeaPCsAAKCYAQAADqiu5Uy1LWOSpH379umWW27RvHnz9Msvv1ww2+bfSkpKrkDaS3OxY8BrU1xcfIWTAACA+nAxOwAAAEBTu/XWW+Xn56e0tDRt3LjR2GulqKhImzdvliT17dtXPj4+xmPKyso0evRonTx5Uq6urpo+fbqGDRsmb29vtW/fXs2bN5ckZWVlqUePHpJ00eKmMZw/f964njNnjkaNGlWvx1k3OAYAAE2LYgYAADikRx55RDNnztTp06e1detWjRw5UvHx8SoqKjLuV7d7925lZWVJkpYsWaIpU6bYfN5Tp041KJeTk5MsFosqKyvrHGfN+W/XXnutce3q6ipfX98G5QEAAI2LpUwAAMAhTZgwQc7OzpKk2NhYSf9bxuTq6qpx48bVGJ+RkWFcjxkzptbnte7Zcrlat24tSXUeX22xWJSZmWnznpeXl9q2bStJ2rt3b4OyAACAxkcxAwAAHFKXLl0UHBwsSdq+fbvS09OVkJAgSXrwwQfVsWPHGuMrKiqM69pmq1RWVmrlypUNyuXp6Smp7oJnx44dOn36tM17zs7OeuihhyRJO3fu1JEjRxqUBwAANC6KGQAA4LCsy5XKy8s1duxYY3+Wfy9jkqRevXoZ1zExMTafb9asWdq/f3+DMt13332SpO+//97mjJcTJ05o+vTpdT7HrFmz5OzsrMrKSo0cOVI5OTm1jj1//rw2bNhQ5xgAANB42GMGAAA4rBEjRqh169Y6e/assVSpffv2GjJkyAVjQ0ND1alTJ/3111+aPXu2srOzNWLECHXo0EGZmZlauXKlEhISFBAQ0KAlRFOnTtWSJUtUUVGhIUOGaM6cOQoMDFRZWZn27t2rRYsWqby8XL169dKvv/5q8zn8/Py0YMECRURE6PDhw/L19dXUqVM1YMAAde7cWaWlpcrOzlZKSoo2bdqkP//8U2lpaerWrdtl5wYAAJeHYgYAADgsDw8PjRw5UmvWrDF+Nnr0aOOEpepatmypdevWafjw4SotLdXy5cu1fPnyGmOCgoL04YcfNmjD3d69e+vdd9/V888/r4KCAkVERNS4f8011+iLL77QK6+8UmsxI0kzZsxQy5YtNWPGDJ05c0bz58/X/PnzbY51c3OTu7v7ZWcGAACXj6VMAADAoYWHh9f43tYyJqvQ0FClpqZq4sSJuv766+Xq6qqOHTvqvvvu04oVK5SQkHBFjp2OiIjQV199pdDQUOMobk9PTz3zzDM6cOCA+vXrV6/nefzxx5WVlaVXX31VAQEB6tChg1xcXNSyZUt5e3srLCxMy5Yt0/Hjx9WzZ88G5wYAAJfOyWKxWMwOAQAAAAAA4IiYMQMAAAAAAGASihkAAAAAAACTUMwAAAAAAACYhGIGAAAAAADAJBQzAAAAAAAAJqGYAQAAAAAAMAnFDAAAAAAAgEkoZgAAAAAAAExCMQMAAAAAAGASihkAAAAAAACTUMwAAAAAAACYhGIGAAAAAADAJBQzAAAAAAAAJqGYAQAAAAAAMAnFDAAAAAAAgEkoZgAAAAAAAExCMQMAAAAAAGASihkAAAAAAACTUMwAAAAAAACYhGIGAAAAAADAJBQzAAAAAAAAJqGYAQAAAAAAMAnFDAAAAAAAgEkoZgAAAAAAAExCMQMAAAAAAGCS/wM4CA5KM+9gOAAAAABJRU5ErkJggg==\n"
},
"metadata": {
"image/png": {
"width": 563,
"height": 454
}
}
}
],
"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": "markdown",
"id": "03048865",
"metadata": {
"id": "03048865"
},
"source": [
"## NumPyroを利用したベイズ推測"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "e35abca8",
"metadata": {
"id": "e35abca8"
},
"outputs": [],
"source": [
"def model(y):\n",
" ap = numpyro.sample(\"ap\", dist.Normal(0, 10))\n",
" al = numpyro.sample(\"al\", dist.Normal(1, 0.5))\n",
" p = expit(ap)\n",
" lambda_ = jnp.exp(al)\n",
" numpyro.sample(\"y\", dist.ZeroInflatedPoisson(p, lambda_), obs=y)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "3c6461a2",
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "3c6461a2",
"outputId": "b6c628f0-a805-4832-efe9-74b51871b806"
},
"outputs": [
{
"output_type": "stream",
"name": "stderr",
"text": [
"sample: 100%|██████████| 3000/3000 [00:05<00:00, 545.58it/s, 3 steps of size 7.62e-01. acc. prob=0.93] \n"
]
},
{
"output_type": "stream",
"name": "stdout",
"text": [
"\n",
" mean std median 5.0% 95.0% n_eff r_hat\n",
" al 1.14 0.05 1.14 1.04 1.22 1566.80 1.00\n",
" ap -0.72 0.16 -0.71 -0.98 -0.45 1158.72 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": "markdown",
"id": "822f8a05",
"metadata": {
"id": "822f8a05"
},
"source": [
"* 頻度がゼロになる確率、それ以外の場合のポワソン分布のrateパラメータ、ともにおおよそ真の値に近い値が得られている。"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "eb19de18",
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "eb19de18",
"outputId": "1d8662ce-67f8-43dc-b704-e4a4ab97bb58"
},
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"0.3292086\n",
"3.119743\n"
]
}
],
"source": [
"post = mcmc.get_samples()\n",
"print(jnp.mean(expit(post[\"ap\"]))) # proportion of zero inflation\n",
"print(jnp.mean(jnp.exp(post[\"al\"]))) # poisson rate"
]
},
{
"cell_type": "code",
"execution_count": 8,
"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