Chapter05StratandPostStrat.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 5: Stratification and Post-Stratification in Randomized Experiments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "matplotlib       : 3.8.0\n",
      "pandas           : 2.1.1\n",
      "scipy            : 1.11.3\n",
      "statsmodels      : 0.13.5\n",
      "matplotlib_inline: 0.1.6\n",
      "numpy            : 1.23.5\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy as sp\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "\n",
    "# viz\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "font = {\"family\": \"IBM Plex Sans Condensed\", \"weight\": \"normal\", \"size\": 10}\n",
    "plt.rc(\"font\", **font)\n",
    "plt.rcParams[\"figure.figsize\"] = (6, 6)\n",
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 1\n",
    "\n",
    "%load_ext watermark\n",
    "%watermark --iversions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>duration</th>\n",
       "      <th>treatment</th>\n",
       "      <th>female</th>\n",
       "      <th>black</th>\n",
       "      <th>hispanic</th>\n",
       "      <th>ndependents</th>\n",
       "      <th>recall</th>\n",
       "      <th>young</th>\n",
       "      <th>old</th>\n",
       "      <th>quarter</th>\n",
       "      <th>durable</th>\n",
       "      <th>lusd</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>18.011343</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>5</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.003399</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>5</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>26.960396</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>7.009044</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>9.022409</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     duration  treatment  female  black  hispanic  ndependents  recall  young  \\\n",
       "1   18.011343          0       0      0         0            2       0      0   \n",
       "4    1.003399          0       0      0         0            0       0      0   \n",
       "5   26.960396          0       0      0         0            0       0      0   \n",
       "6    7.009044          1       0      0         0            0       0      0   \n",
       "12   9.022409          1       0      0         0            0       0      1   \n",
       "\n",
       "    old  quarter  durable  lusd  \n",
       "1     0        5        0     0  \n",
       "4     0        5        0     1  \n",
       "5     0        4        0     1  \n",
       "6     0        2        0     0  \n",
       "12    0        3        0     0  "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "penndata = pd.read_table(\"Penn46_ascii.txt\", sep=\"\\s+\")\n",
    "y, z, block = (\n",
    "    np.log(penndata.duration).values,\n",
    "    penndata.treatment.values,\n",
    "    penndata.quarter.values,\n",
    ")\n",
    "penndata.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>col_0</th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "      <th>5</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>row_0</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>234</td>\n",
       "      <td>41</td>\n",
       "      <td>687</td>\n",
       "      <td>794</td>\n",
       "      <td>738</td>\n",
       "      <td>860</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>87</td>\n",
       "      <td>48</td>\n",
       "      <td>757</td>\n",
       "      <td>866</td>\n",
       "      <td>811</td>\n",
       "      <td>461</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "col_0    0   1    2    3    4    5\n",
       "row_0                             \n",
       "0      234  41  687  794  738  860\n",
       "1       87  48  757  866  811  461"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.crosstab(z, block)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## FRT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def stat_SRE(z, y, x):\n",
    "    xlevels = np.unique(x)\n",
    "    K = len(xlevels)\n",
    "    PiK, TauK, Wk = np.zeros(K), np.zeros(K), np.zeros(K)\n",
    "    for k in range(K):\n",
    "        id = np.where(x == xlevels[k])\n",
    "        zk = z[id]\n",
    "        yk = y[id]\n",
    "        PiK[k] = zk.shape[0] / z.shape[0]\n",
    "        TauK[k] = np.mean(yk[zk == 1]) - np.mean(yk[zk == 0])\n",
    "        Wk[k] = sp.stats.mannwhitneyu(yk[zk == 1], yk[zk == 0])[0]\n",
    "    return np.sum(PiK * TauK), sum(Wk / PiK)\n",
    "\n",
    "\n",
    "def zRandomSRE(z, x):\n",
    "    xlevels = np.unique(x)\n",
    "    K = len(xlevels)\n",
    "    zrandom = z.copy()\n",
    "    for k in range(K):\n",
    "        xk = xlevels[k]\n",
    "        zrandom[x == xk] = np.random.permutation(z[x == xk])\n",
    "    return zrandom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.001\n",
      "0.0\n"
     ]
    }
   ],
   "source": [
    "# observed test statistics\n",
    "stat_obs = stat_SRE(z, y, block)\n",
    "# null distribution\n",
    "MC = int(1e3)\n",
    "statSREMC = np.zeros((MC, 2))\n",
    "for k in range(MC):\n",
    "    zrandom = zRandomSRE(z, block)\n",
    "    statSREMC[k] = stat_SRE(zrandom, y, block)\n",
    "\n",
    "print(np.mean(statSREMC[:, 0] <= stat_obs[0]))\n",
    "print(np.mean(statSREMC[:, 1] <= stat_obs[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA/sAAAQbCAYAAADXm7s2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAB7CAAAewgFu0HU+AABm60lEQVR4nO3de7hWZZ0//veGzRlBQAXkIDn4pSgbatIxUzxUlpqig02NVpozkYcxvzTpUKagaEN9S3GSxEN5yizPh8iyJvOUkzo5pXkYwQOQm1FAQOQMz+8Pfnsnujls2M8+3Lxe17Wva+9n3eten2d7u3je+77XWjWVSqUSAAAAoBgdWrsAAAAAoHkJ+wAAAFAYYR8AAAAKI+wDAABAYYR9AAAAKIywDwAAAIUR9gEAAKAwwj4AAAAURtgHAACAwgj7AAAAUBhhHwAAAAoj7AMAAEBhhH0AAAAojLAPAAAAhRH2AQAAoDDCPgAAABRG2AcAAIDC1LZ2AW3NihUr8sQTTyRJdt5559TW+hUBAABQPWvWrMmrr76aJNlzzz3TtWvXbe5Tkn2LJ554InvvvXdrlwEAAMB26JFHHslee+21zf1Yxg8AAACFMbP/FjvvvHPD94888kgGDhzYitUAAABQurq6uoYV5m/OpNtC2H+LN1+jP3DgwAwePLgVqwEAAGB70lz3jbOMHwAAAAoj7AMAAEBhhH0AAAAojLAPAAAAhRH2AQAAoDDCPgAAABRG2AcAAIDCCPsAAABQGGEfAAAACiPsAwAAQGGEfQAAACiMsA8AAACFEfYBAACgMMI+AAAAFEbYBwAAgMII+wAAAFAYYR8AAAAKU9vaBQDA1ho2YcY29/HilMOboRIAgLbFzD4AAAAURtgHAACAwgj7AAAAUBhhHwAAAAoj7AMAAEBhhH0AAAAojLAPAAAAhRH2AQAAoDDCPgAAABRG2AcAAIDCCPsAAABQGGEfAAAACiPsAwAAQGGEfQAAACiMsA8AAACFEfYBAACgMMI+AAAAFEbYBwAAgMII+wAAAFAYYR8AAAAKI+wDAABAYYR9AAAAKExtaxcAwPZp2IQZrV1Ckuap48UphzdDJQAAzcfMPgAAABRG2AcAAIDCCPsAAABQGGEfAAAACiPsAwAAQGGEfQAAACiMsA8AAACFEfYBAACgMMI+AAAAFEbYBwAAgMII+wAAAFAYYR8AAAAKI+wDAABAYYR9AAAAKEyLhv2lS5fmqKOOSk1NTQYMGJAHH3xwk+0feuihHHPMMdl1113TpUuXDBw4MGPHjt3sfgAAALA9a7Gw//rrr+eQQw7JHXfckXe/+9159NFHs99++220/RVXXJHRo0fnlltuSV1dXTp27Jh58+bl1ltvzQEHHJDp06e3VOkAAADQrrRI2F+3bl0++clP5uGHH8773//+PPDAAxkyZMhG2z/77LM59dRTs27dupx99tmZP39+li1blgULFmTixIlZt25dTjvttDz11FMtUT4AAAC0Ky0S9i+44IL84he/yO6775577rknffr02WT7iy66KKtXr864ceNy3nnnpV+/fkmSvn37ZtKkSTn55JOzZs2afPvb326J8gEAAKBdqXrYf/bZZ3P++eenY8eOueWWWxqC+6bceeedSZLTTz+90e3jx49Pktx1112pVCrNVywAAAAUoOph/5xzzsmqVasybty4jBo1arPtZ82albq6uvTp0ycjR45stM0ee+yRAQMGZP78+ZbyAwAAwFtUNew///zzufnmm9O7d++cd955W7TPc889lySbvKb/zdtnzpy5bUUCAABAYWqr2fl1112XdevW5bOf/WyWLl2aL3/5y7nvvvuyYMGCDB06NJ/85Cfzla98JTvssEPDPrNnz06SdOvWbZN9d+/efYP2W2ru3Lmb3F5XV9ek/gAAAKCtqWrYv/XWW5Mk++yzT/bdd98NgvTTTz+d8847L7fddlvuu+++hpv2LV26tEnHaGr7za0YAAAAgPauasv4X3/99Tz55JPp0qVLbrrpppxwwgl54YUXsmrVqrz00ks5++yz06FDhzzxxBM588wzG/Zbvnx5k46zbNmy5i4dAAAA2rWqzew/8cQTWbduXVauXJkRI0bkG9/4RsO2oUOHNlzDP3ny5Pzwhz/MhRdeuMFy/mqZM2fOJrfX1dVl7733rnodAAAAUC1VC/vz589v+P6UU05ptM0pp5ySyZMnZ8WKFfn973+fAw44YLPX6r9V/bX7W2rw4MFNag8AAADtTdWW8S9ZsiRJ0rNnz+y2226NthkwYED69euX5C83zuvZs2eTjtPU9gAAAFC6qoX9+hn6Hj16bLJdfVhfsWJFkvVL/N/888bUX6tf3x4AAABYr2phv3///kk2f8O9N954I0nSu3fvJMnw4cOTbP7a+vrt9e0BAACA9aoW9t/73vemQ4cOWbJkSf785z832mb+/PkN1/aPHDkyyfrwPnDgwCxcuDDPPPNMo/vNmjUr8+bNS79+/Rr2AwAAANarWtjv1atXDjrooCTJFVdc0Wib+teHDRuWd73rXQ2vH3nkkUmSiy++uNH9pk6dmiQ54ogjUlNT01wlAwAAQBGqFvaT5Gtf+1qS5IILLsi3vvWtLFq0KEmyePHiXHTRRZk4cWKSZMKECRuE9vHjx6dTp06ZPn16Jk2alIULFyZJXnvttUyePDnTpk1LbW1tzjjjjGqWDwAAAO1SVcP+wQcfnPPPPz9r1qzJv/7rv6Zv377p2bNn+vTpky9/+ctZvXp1TjzxxIwbN26D/UaMGJFLLrkkHTp0yLnnnpt+/fqlR48e6du3b84555wkyXe/+11L+AEAAKARVQ37SXLWWWfl17/+dcaMGZOdd945q1atSr9+/fKxj30sN998c77//e83uhR/3LhxuffeezNmzJgMGDAgq1evTv/+/TNmzJjce++9Oemkk6pdOgAAALRLNZVKpdLaRbQlc+fOzZAhQ5Ksv+P/4MGDW7kigDINmzCjtUtoU16ccnhrlwAAtJJq5NCqz+wDAAAALUvYBwAAgMII+wAAAFAYYR8AAAAKI+wDAABAYYR9AAAAKIywDwAAAIUR9gEAAKAwwj4AAAAURtgHAACAwgj7AAAAUBhhHwAAAAoj7AMAAEBhhH0AAAAojLAPAAAAhRH2AQAAoDDCPgAAABRG2AcAAIDCCPsAAABQGGEfAAAACiPsAwAAQGGEfQAAACiMsA8AAACFEfYBAACgMMI+AAAAFEbYBwAAgMII+wAAAFAYYR8AAAAKI+wDAABAYWpbuwAA2p9hE2a0dgkAAGyCmX0AAAAojLAPAAAAhRH2AQAAoDDCPgAAABRG2AcAAIDCuBs/wHbGnfTZlOYYHy9OObwZKgEAtoWZfQAAACiMsA8AAACFEfYBAACgMMI+AAAAFEbYBwAAgMII+wAAAFAYYR8AAAAKI+wDAABAYYR9AAAAKIywDwAAAIUR9gEAAKAwta1dAC3k1VeTXXbZ8LVXXkl23rl16gFgA8MmzNjmPl6ccngzVAIAlMDMPgAAABRG2AcAAIDCCPsAAABQGGEfAAAACiPsAwAAQGGEfQAAACiMsA8AAACFEfYBAACgMMI+AAAAFEbYBwAAgMII+wAAAFAYYR8AAAAKI+wDAABAYYR9AAAAKIywDwAAAIUR9gEAAKAwwj4AAAAURtgHAACAwgj7AAAAUBhhHwAAAAoj7AMAAEBhhH0AAAAojLAPAAAAhRH2AQAAoDDCPgAAABRG2AcAAIDCCPsAAABQGGEfAAAAClPb2gUAsOWGTZjR2iUAANAOmNkHAACAwgj7AAAAUBhhHwAAAAoj7AMAAEBhhH0AAAAojLvxA0AhPK0BAKhnZh8AAAAKI+wDAABAYYR9AAAAKIywDwAAAIUR9gEAAKAwwj4AAAAURtgHAACAwgj7AAAAUBhhHwAAAAoj7AMAAEBhhH0AAAAoTG1rFwAAlGXYhBnbtP+LUw5vpkoAYPtlZh8AAAAK06Jhf8GCBenVq1dqampSU1OTSZMmbbL9k08+meOPPz5Dhw5Nly5dsssuu+TQQw/NXXfd1TIFAwAAQDvUomF/ypQpef3111Nbu/mrB372s5/lAx/4QK699trMmTMnHTt2zKuvvpqf//znOfLII3PWWWe1QMUAAADQ/rRY2K+rq8u0adMyatSoHHvssZts+9prr+W4447LypUr84UvfCFz587NsmXL8vrrr+eSSy5J586d841vfCP33HNPC1UPAAAA7UeL3aDv/PPPz/LlyzNx4sTcfvvtm2x72WWXZdGiRTnkkENy+eWXN7zes2fPnHrqqXnjjTfyr//6r5kyZUoOOeSQKlcOALSkbb3BX+ImfwDQIjP7L730Uq688sqMGjUqY8aM2Wz7O+64I0ly+umnN7r9lFNOSZcuXXL//fdn0aJFzVkqAAAAtHstEvbPPffcrFq1KhMnTkxNTc0m265YsSKPPfZYkmS//fZrtE3Pnj3zvve9L2vXrs1DDz3U7PUCAABAe1b1sP/ss8/m2muv3eJZ/RdffDFr1qxJr1690qtXr422GzJkSJJk5syZzVYrAAAAlKDq1+xPnDgxa9eu3aJZ/SSZPXt2kqRbt26bbNe9e/cN2m+puXPnbnJ7XV1dk/oDAACAtqaqYf+Pf/xjbrzxxi2e1U+SpUuXNukYTW1fvyIAAAAASlXVsH/22WenUqls8ax+kixfvrxJx1i2bNnWlAbQKprjLuMAALA5VQv7jzzySO68884mzeq3hDlz5mxye11dXfbee+8WqgYAAACaX9XC/llnnZUkTZrVTzZ/rf5b1V+7v6UGDx7cpPYAAADQ3lTlbvz33XdffvWrX23VrH7Pnj2r2h4AAABKV5Ww//Wvfz1J02f1k2To0KFJkhUrVmyyXf21+vXtAQAAgPWqsoz/wQcfTJKccMIJjW6vD+pTpkzJ1KlTkySLFi1KkgwbNiwdO3bM4sWLs2TJkvTq1avRPuqvvR8+fHjzFQ4AAAAFqMrMfr3Fixc3+rV69eokycqVKxteq9e1a9fstddeSZLf/va3jfb7xhtv5PHHH0+HDh2y7777VvMtAAAAQLtTlbBfqVQ2+XX88ccnWb/Mv/61N6u/zv/iiy9utP/p06dn5cqVGT16dPr06VONtwAAAADtVlVn9rfWuHHj0rt37/z85z/PSSedlJdffjnJ+hn9Sy+9NF/72teSJF/96ldbs0wAAABok9pk2O/bt2+uv/76dO7cOZdddlkGDRqUHj16ZIcddsgpp5ySVatWZcKECTnkkENau1QAAABoc9pk2E+Sww8/PI8++miOPfbYDBo0KGvWrEnfvn3z0Y9+NLfddlv+7d/+rbVLBAAAgDapKnfj35yrr746V1999Wbbvfe97831119f/YIAAACgIG12Zh8AAADYOsI+AAAAFEbYBwAAgMII+wAAAFAYYR8AAAAKI+wDAABAYYR9AAAAKIywDwAAAIUR9gEAAKAwwj4AAAAURtgHAACAwgj7AAAAUJja1i4AAKC5DZswY5v7eHHK4c1QCQC0DjP7AAAAUBhhHwAAAAoj7AMAAEBhhH0AAAAojLAPAAAAhRH2AQAAoDDCPgAAABRG2AcAAIDCCPsAAABQGGEfAAAACiPsAwAAQGGEfQAAACiMsA8AAACFEfYBAACgMMI+AAAAFEbYBwAAgMII+wAAAFAYYR8AAAAKI+wDAABAYYR9AAAAKIywDwAAAIUR9gEAAKAwwj4AAAAURtgHAACAwgj7AAAAUBhhHwAAAAoj7AMAAEBhhH0AAAAojLAPAAAAhRH2AQAAoDDCPgAAABRG2AcAAIDCCPsAAABQGGEfAAAAClPb2gUAALBxwybM2OY+XpxyeDNUAkB7YmYfAAAACiPsAwAAQGGEfQAAACiMsA8AAACFEfYBAACgMO7GD7R57kQNAABNY2YfAAAACiPsAwAAQGGEfQAAACiMsA8AAACFEfYBAACgMO7GDwDQCE8CAaA9M7MPAAAAhRH2AQAAoDDCPgAAABRG2AcAAIDCCPsAAABQGGEfAAAACiPsAwAAQGGEfQAAACiMsA8AAACFEfYBAACgMMI+AAAAFEbYBwAAgMII+wAAAFAYYR8AAAAKI+wDAABAYYR9AAAAKIywDwAAAIUR9gEAAKAwwj4AAAAURtgHAACAwgj7AAAAUBhhHwAAAAoj7AMAAEBhhH0AAAAojLAPAAAAhRH2AQAAoDDCPgAAABRG2AcAAIDC1LZ2AQDtxbAJM1q7BKCdcd4AoLWY2QcAAIDCCPsAAABQGGEfAAAACiPsAwAAQGGEfQAAACiMsA8AAACFEfYBAACgMMI+AAAAFKZFwv5jjz2WsWPHZpdddkmnTp2y884757DDDsvdd9+9yf0eeuihHHPMMdl1113TpUuXDBw4MGPHjs2DDz7YEmUDAABAu1T1sH/OOedkr732yq233ppXX301HTp0yPz583P33XfnsMMOy7e//e1G97viiisyevTo3HLLLamrq0vHjh0zb9683HrrrTnggAMyffr0apcOAAAA7VJVw/5VV12VyZMnp3PnzvnmN7+Z+fPnZ+XKlXn++eczZsyYJMlXv/rVzJo1a4P9nn322Zx66qlZt25dzj777MyfPz/Lli3LggULMnHixKxbty6nnXZannrqqWqWDwAAAO1SbbU6XrNmTSZNmpQkueiii3LKKac0bHvHO96RG264IbvvvnvmzZuXG264IV//+tcbtl900UVZvXp1xo0bl/POO6/h9b59+2bSpEl55ZVXcumll+bb3/52fvCDH1TrLQAFGTZhRmuXAAAALaZqM/urV6/Oqaeemn333Tfjxo172/Zu3brlwAMPTJL84Q9/2GDbnXfemSQ5/fTTG+17/PjxSZK77rorlUqlGasGAACA9q9qYb9bt24588wz89BDD6W2tvEFBIMHD06SLF68uOG1WbNmpa6uLn369MnIkSMb3W+PPfbIgAEDMn/+fEv5AQAA4C1a9dF7q1atSpLsuOOODa8999xzSZIhQ4Zsct/67TNnzqxOcQAAANBOVe2a/S3x9NNPJ1l/DX+92bNnJ1m/MmBTunfvvkH7LTV37txNbq+rq2tSfwAAANDWtFrYf/nll3PvvfcmST7xiU80vL506dIm9dPU9ptbMQAAAADtXauF/a985StZs2ZN9t9//+y///4Nry9fvrxJ/Sxbtqy5SwP+f81xB/sXpxzeDJUAAABN0Sph//vf/35uuOGG9OjRI1dccUWLHnvOnDmb3F5XV5e99967haoBAACA5tfiYf+ee+7JySefnJqamlxzzTUZMWLEBts3d63+W9Vfu7+l6p8AAAAAAKVq0bD/wAMP5Oijj87q1atz4YUXZuzYsW9r07Nnzyb12dT2AAAAULoWe/TeY489lk984hNZtmxZvvrVr2b8+PGNths6dGiSZMWKFZvsr/5a/fr2AAAAwHotEvaffPLJfPzjH8+SJUty0kkn5Rvf+MZG2w4fPjzJ5q+tr99e3x4AAABYr+phf+bMmfnoRz+aBQsW5NOf/nSmTZu2yfbDhw/PwIEDs3DhwjzzzDONtpk1a1bmzZuXfv36ZeTIkdUoGwAAANqtqob92bNn58Mf/nDmzZuXQw89NNdee206dNj8IY888sgkycUXX9zo9qlTpyZJjjjiiNTU1DRbvQAAAFCCqoX9pUuX5iMf+Uhmz56d/fbbL7fccks6deq0RfuOHz8+nTp1yvTp0zNp0qQsXLgwSfLaa69l8uTJmTZtWmpra3PGGWdUq3wAAABot6oW9ufPn5/nnnsuSfL4449n4MCB2XHHHTf69aMf/ahh3xEjRuSSSy5Jhw4dcu6556Zfv37p0aNH+vbtm3POOSdJ8t3vftcSfgAAAGhEizx674033thsm1WrVm3w87hx4/LOd74zF154YX73u99lwYIF6d+/f/bZZ5+MHz8+BxxwQLXKBQAAgHatamF/2LBhqVQq29TH6NGjM3r06GaqCAAAALYPLfLoPQAAAKDlCPsAAABQGGEfAAAACiPsAwAAQGGEfQAAACiMsA8AAACFEfYBAACgMMI+AAAAFEbYBwAAgMII+wAAAFAYYR8AAAAKI+wDAABAYWpbuwCgbMMmzGjtEgC2e81xLn5xyuHNUAkALcXMPgAAABRG2AcAAIDCCPsAAABQGGEfAAAACiPsAwAAQGGEfQAAACiMsA8AAACFEfYBAACgMMI+AAAAFEbYBwAAgMII+wAAAFAYYR8AAAAKI+wDAABAYYR9AAAAKIywDwAAAIUR9gEAAKAwta1dAFAdwybMaO0SAACAVmJmHwAAAAoj7AMAAEBhhH0AAAAojLAPAAAAhRH2AQAAoDDuxg8AwGY1x1NeXpxyeDNUAsCWMLMPAAAAhRH2AQAAoDDCPgAAABRG2AcAAIDCCPsAAABQGGEfAAAACiPsAwAAQGGEfQAAACiMsA8AAACFEfYBAACgMMI+AAAAFEbYBwAAgMII+wAAAFAYYR8AAAAKI+wDAABAYYR9AAAAKIywDwAAAIUR9gEAAKAwwj4AAAAURtgHAACAwgj7AAAAUBhhHwAAAAoj7AMAAEBhalu7ACjNsAkzWrsEAGiTmuPfyBenHF5MHQDVZGYfAAAACiPsAwAAQGGEfQAAACiMsA8AAACFEfYBAACgMMI+AAAAFEbYBwAAgMII+wAAAFAYYR8AAAAKI+wDAABAYWpbuwBoS4ZNmNHaJQAAm+DfaoAtY2YfAAAACiPsAwAAQGGEfQAAACiMsA8AAACFEfYBAACgMMI+AAAAFEbYBwAAgMII+wAAAFAYYR8AAAAKI+wDAABAYYR9AAAAKExtaxcAzWXYhBmtXQIAwBZpjs8tL045vBkqAUplZh8AAAAKI+wDAABAYYR9AAAAKIywDwAAAIUR9gEAAKAwwj4AAAAURtgHAACAwgj7AAAAUBhhHwAAAAoj7AMAAEBhhH0AAAAojLAPAAAAhRH2AQAAoDBtPuw/9NBDOeaYY7LrrrumS5cuGThwYMaOHZsHH3ywtUsDAACANqlNh/0rrrgio0ePzi233JK6urp07Ngx8+bNy6233poDDjgg06dPb+0SAQAAoM1ps2H/2Wefzamnnpp169bl7LPPzvz587Ns2bIsWLAgEydOzLp163Laaaflqaeeau1SAQAAoE1ps2H/oosuyurVqzNu3Licd9556devX5Kkb9++mTRpUk4++eSsWbMm3/72t1u5UgAAAGhbalu7gI258847kySnn356o9vHjx+fSy+9NHfddVcqlUpqampasrw2Y9iEGVvUru+yxfn9W157/+RfZmH33s1Sx4tTDt/mPrb0vQAA0DyfnXyG21Bz/D5oXm1lnLdHbXJmf9asWamrq0ufPn0ycuTIRtvsscceGTBgQObPn28pPwAAALxJmwz7zz33XJJkyJAhm2xXv33mzJlVrwkAAADaizYZ9mfPnp0k6dat2ybbde/efYP2AAAAQBu9Zn/p0qVVaz937txNbp8zZ07D93V1dU2qozWsWTJ/i9qtWrEkb33nq5YuyJo1q5uljs39XrfElr4XAIDWtq2ffdrK5x6f4TbUHL8PmldzjK/28N/1zdlzzZo1zdJnmwz7y5cvb1L7ZcuWbXHbzV0a8GZ77713k+poy/6c5G3v/KovNVv/Qy5ttq4AANq8Uj77lPI+movfR5na23/XV199NcOGDdvmftrkMn4AAABg67XJmf3NXav/VvXX7m+JNy/Tb8yKFSvyzDPPpH///tl5551TW9v0X1FdXV3DqoBHHnkkAwcObHIf0J4Y82xPjHe2J8Y72xPjnda0Zs2avPrqq0mSPffcs1n6bJNhv2fPnlVrP3jw4M22GT58eJOOvykDBw7comNCKYx5tifGO9sT453tifFOa2iOpftv1iaX8Q8dOjTJ+ln2Tam/Vr++PQAAANBGw379zPrmltzXb2/OmXgAAABo79ps2B84cGAWLlyYZ555ptE2s2bNyrx589KvX7+MHDmyhSsEAACAtqtNhv0kOfLII5MkF198caPbp06dmiQ54ogjUlNT01JlAQAAQJvXZsP++PHj06lTp0yfPj2TJk3KwoULkySvvfZaJk+enGnTpqW2tjZnnHFGK1cKAAAAbUubDfsjRozIJZdckg4dOuTcc89Nv3790qNHj/Tt2zfnnHNOkuS73/2uJfwAAADwFm027CfJuHHjcu+992bMmDEZMGBAVq9enf79+2fMmDG59957c9JJJ7V2iQAAANDm1FQqlUprFwEAAAA0nzY9sw8AAAA0nbAPAAAAhRH2AQAAoDDCPgAAABRG2AcAAIDCCPsAAABQGGEfAAAACiPsAwAAQGGEfQAAACiMsL8FHnrooRxzzDHZdddd06VLlwwcODBjx47Ngw8+uE39Ll26NAcffHBqamoavp555pkt2vfaa6/NgQcemJ122ildu3bN7rvvnpNPPjmzZ8/eppqgWuN9a8bspEmTNvj/Y2Nf//mf/7lNtVG2aozpbTkHO39TTW1lvDt/01KqMeZ9RqcYFTbp8ssvr3To0KGSpJKk0q1bt4bvO3ToULn00ku3qt9FixZVPvjBD1aSVAYNGtTQ59NPP73ZfT/72c82tK+pqal07dq14edevXpVHn300a2qCao13rd2zE6cOLGhjt69e2/0y5hnY6oxprflHOz8TTW1pfHu/E1LqMaY9xmdkgj7m/DMM89UOnXqVElSOfvssyvz58+vVCqVyoIFCxr+Eautra386U9/alK/r776auV973tfJUlln332qTz++ONbfCL5wQ9+0HDCuO666yrLly+vVCqVysyZMytHH310JUll2LBhlVWrVm3dm2a7Va3xvi1jtv64V1111Ta/P7Y/1RjT2zKenb+pprY23p2/qbZqjHmf0SmNsL8JX/ziFytJKuPGjWt0+8knn1xJUvn85z+/xX2+/PLLlZEjR1aSVA488MDK66+/XnnhhRe2+EQyYsSISpLKj370o7dtW7VqVeW9731vJUnlmmuu2eKaoFKpznivVLZtzPqwyLaoxpjelvHs/E01tbXx7vxNtTX3mPcZnRIJ+5swcODASpKN/kXwf/7nfypJKjvttFNl3bp1W9TnFVdcUUlSOeywwxr+4relJ5JnnnmmkqTSv3//ytq1axttc/nll1eSVP7u7/5ui+qBetUY79s6Zn1YZFs095jelvHs/E21taXxXqk4f1N9zT3mfUanRG7QtxGzZs1KXV1d+vTpk5EjRzbaZo899siAAQMyf/78PPXUU1vU7z/90z/lhz/8YW6//fZ07dq1STXV32hk3333TYcOjf+nGz16dJLk/vvvb1LfbN+qNd6NWVpLNcb0toxn/y9QTW1tvEO1VWPM+4xOiYT9jXjuueeSJEOGDNlku/rtM2fO3OK+jzvuuHTq1KkqNdVvmz9/fhYvXtzkY7B9qtZ4N2ZpLdUY09synv2/QDW1tfEO1Vatzy0+o1MaYX8j6h+P0a1bt0226969+wbtW7um+npaqibKUK3x3pxj9k9/+lPGjh2b3XbbLUOHDs2RRx6Z//iP/9iiOtj+VGNMb8t4dv6mmtraeH8r52+aW1v7nO4cT1sl7G/E0qVLq9p+a7TFmihDtcZWc/X74osvZv/998+tt96a2bNnZ86cObnrrrvykY98JBdeeGGTjsH2oRpjelv6dP6mmtraeH8z52+qoa2dU9taPVBP2N+I5cuXN6n9smXLqlTJX7TFmihDtcZWc/U7derUnHDCCZkzZ05WrlyZJ598MkcffXSS5Mwzz8wf/vCHJh2H8lVjTG9Ln87fVFNbG+9v5vxNNbS1c2pbqwfq1bZ2AdV222235atf/eoWt7/yyiuz3377VbEiqJ7Sxvvo0aPz2c9+Nv/n//yffP3rX294/d3vfnduuummjBgxIrNmzcq0adNy+eWXt2KlALyZ8zdA6ys+7C9evDjPPvvsFrevX1azuWuA3urN1+FUS1usibalrY33be334IMPzsEHH9xo244dO+af/umf8tWvfjX33ntvk45D+aoxprelT+dvqqmtjffE+Zvqamvn1LZWD9Qrfhn/CSeckEqlssVfH//4x5MkPXv2bNJxmtp+a7TFmmhb2tp4r/aYHT58eBI3uuHtqjH2tqVP52+qqa2N9y3h/M22aGvn1LZWD9QrPuxvraFDhyZJVqxYscl29dfc1Ldv7ZrefA1QS9REGao13qs9ZnfaaackyapVq7J27dom7UvZqjGmt2U8O39TTW1tvG8J52+2RVv7nO4cT1sl7G9E/V+c58yZs8l29dvr27d2TfXbdtppp/Tu3bvqNVGGao33ao/ZBQsWJFn/F/KOHTs2aV/KVo0xvS3j2fmbampr431LOH+zLdra53TneNoqYX8jhg8fnoEDB2bhwoV55plnGm0za9aszJs3L/369cvIkSOrXtP++++fJHn44YdTqVQabfPggw9u0Ba2RLXG+7aO2Zdffjk//elPN9r/rFmzkiTvfOc7t6geth/VGNPbMp6dv6mmtjbeE+dvqqutfU53jqetEvY34cgjj0ySXHzxxY1unzp1apLkiCOOSE1NTdXrGTFiREaMGJG6urrceOONb9u+Zs2aXHLJJUmSMWPGVL0eylKN8b4tY/aFF17I+9///nzmM5/Jc88997Z9165dmyuvvDJJcuihh25RPWxfmntMb8t4dv6m2trSeHf+piW0pc/pzvG0WRU26plnnql06tSpkqQyceLEyoIFCyqVSqWycOHCynnnnVepqamp1NbWVv70pz9tsN+iRYsq73vf+ypdu3atfPazn93scV544YVKkkqSytNPP73Jtt///vcrSSq9e/euXH/99ZUVK1ZUKpVKZdasWZWxY8dWklR22223ysqVK7fyXbO9qtZ439oxu2bNmsoHP/jBSpLKrrvuWvnJT35SWbZsWaVSqVSef/75yqc+9alKksqOO+5Yefnll6vxK6Gd25oxXa3xvK37wua0pfHu/E1LqMaYfyuf0WnvhP3NuOyyyyodOnRo+B+9e/fuDd/X1NRULr300rft88ADDzS0qa2tfdv2jh07vu2rvv1bXz/33HPftv9xxx3X0L5Dhw6Vrl27Nvzcq1evyiOPPFKV3wXlq8Z4r1S2fsy++uqrlf3222+DGrp167ZBfT//+c+b9XdAWZo6pqs5nrd1X9ictjTenb9pCc095n1GpzTC/ha47777KmPGjKkMGDCg0qlTp0r//v0rY8aMqfzmN79ptP2iRYsqo0aNqnTp0qXRvxrW/0+/JV8TJ05s9Bjf//73K/vtt1+lT58+lS5dulR22223yhe+8IXKCy+80IzvnO1Rc4/3els7ZletWlW55pprKqNHj67ssssulS5dulSGDRtWGTduXGXmzJnb8lbZTjRlTFd7PG/rvrA5bWm8O3/TEppzzPuMTmlqKpWN3EUCAAAAaJfcoA8AAAAKI+wDAABAYYR9AAAAKIywDwAAAIUR9gEAAKAwwj4AAAAURtgHAACAwgj7AAAAUBhhHwAAAAoj7AMAAEBhhH0AAAAojLAPAAAAhRH2AQAAoDDCPgAAABRG2AcAAIDCCPsAAABQGGEfAAAAClPb2gUAAO3T0qVL853vfCeVSiVf+MIXMmjQoNYuCQD4/wn7AMBWOffcc/Ptb387SfLUU0/lxhtvbOWKAIB6NZVKpdLaRQAA7cuf/vSnjBo1KjvuuGP69OmT5557Lvfcc08++tGPtnZpAECEfQBgKxx44IG57777cvPNN2fQoEHZb7/98ld/9Vd54okn0rlz59YuDwC2e27QBwA0yQ9/+MPcd999OfbYYzN27Njss88+OfPMM/M///M/+X//7/+1dnkAQMzsAwBNsHjx4rzzne9Mhw4d8uSTT6ZPnz5JklWrVmWvvfbKc889l6eeeirDhg1r3UIBYDsn7AMAAEBhLOMHAACAwgj7AAAAUBhhHwAAAAoj7AMAAEBhhH0AYJPuvvvu1NTUNPnr0EMPbe3SAWC7JewDAJv08MMPb9V+++67bzNXAgBsKY/eAwA2ae7cuVm0aNEGr61cuTJ77bVXKpVKvvCFL+RLX/rS2/YbMmRIevfu3UJVAgBvVtvaBQAAbdvgwYMzePDgDV57/PHHUz9fcMABB+Q973lPa5QGAGyEZfwAQJM98cQTDd+/973vbcVKAIDGCPsAQJP98Y9/TJJ06tQp73znO1u5GgDgrYR9AKDJ6mf23/nOd6ZTp06tXA0A8FbCPgDQZPUz+3vuuWcrVwIANEbYBwCaZP78+Zk3b14S1+sDQFsl7AMATeLmfADQ9gn7AECTPPPMMw3fv/vd796ifX7xi1/k0EMPzdChQ9OlS5f0798/BxxwQK655ppqlQkA2zVhHwBokpdeeilJ0qFDhwwePHiz7c8444x8/OMfz6xZs3LEEUdk/Pjx+djHPpYXXnghjz/+eLXLBYDtUm1rFwAAtC+rV69OktTU1KRSqWyy7fPPP5/vfOc7Oeqoo3LLLbekQ4e/zDOsW7cuixYtqmapALDdEvbfYsWKFQ3XIu68886prfUrAoA369WrV5Jk7dq1+eIXv5ijjjoqXbt2zV/91V+97TF8DzzwQCqVSgYPHpyXX3650f6WLVtW9ZoBoC1bs2ZNXn311STrn3TTtWvXbe6zprK5P8lvZx599NHsvfferV0GAAAA26FHHnkke+211zb345p9AAAAKIw16m+x8847N3z/yCOPZODAga1YDQAAAKWrq6trWGH+5ky6LYT9t3jzNfoDBw7corsMAwAAQHNorvvGWcYPAAAAhRH2AQAAoDDCPgAAABRG2AcAAIDCCPsAAABQGGEfAAAACiPsAwAAQGGEfQAAACiMsA8AAACFEfYBAACgMMI+AAAAFEbYBwAAgMII+wAAAFAYYR8AAAAKI+wDAABAYYR9AAAAKIywDwAAAIWpbe0CAICtM2zCjBY5zotTDm+R4wAAzcfMPgAAABRG2AcAAIDCCPsAAABQGGEfAAAACiPsAwAAQGGEfQAAACiMsA8AAACFEfYBAACgMMI+AAAAFEbYBwAAgMII+wAAAFAYYR8AAAAKU9vaBQBAiYZNmNHaJdCIlvrv8uKUw1vkOACwMWb2AQAAoDDCPgAAABRG2AcAAIDCCPsAAABQGGEfAAAACiPsAwAAQGGEfQAAAChMbWsXAAAtqaWesw4A0JrM7AMAAEBhhH0AAAAojLAPAAAAhRH2AQAAoDDCPgAAABRG2AcAAIDCCPsAAABQGGEfAAAACiPsAwAAQGGEfQAAACiMsA8AAACFEfYBAACgMMI+AAAAFEbYBwAAgMII+wAAAFAYYR8AAAAKI+wDAABAYYR9AAAAKIywDwAAAIUR9gEAAKAwwj4AAAAUpuphf/HixZk8eXL23nvv7LjjjunUqVP69euX/fffP//+7/+elStXbnTfn/70pzn00EPTv3//dOnSJUOGDMnnPve5PPnkk9UuGwAAANqt2mp2/sILL+SAAw7InDlzGl7r0qVLFi5cmAcffDAPPvhgrr766vzHf/xH+vTps8G+Z599ds4///yGn7t165a5c+fmuuuuy4033pibbropRxxxRDXLBwAAgHapqjP7p59+eubMmZNevXrlmmuuyZIlS7JixYosWLAgU6dOTW1tbR5//PFMmjRpg/1+/etf5/zzz0+nTp1y8cUXZ8mSJVm2bFn+/Oc/56STTsrKlSvzmc98JvPnz69m+QAAANAuVS3sr1y5Mj//+c+TJOeff34+97nPZYcddkiS9O3bN6effnpOO+20JMntt9++wb7f+ta3kiTnnntuvvSlLzXst+uuu+bSSy/NYYcdliVLluR73/tetcoHAACAdqtqYX/hwoVZvXp1kuRDH/pQo23qX583b17Da0uXLs2vf/3r1NbW5pRTTml0v//7f/9vkuSOO+5oxooBAACgDFUL+zvvvHNqa7fslgC77rprw/e/+93vsnr16uy5557p3bt3o+3322+/1NTU5L//+7/z+uuvN0u9AAAAUIqqhf3a2tqGmfuHHnqo0Tb1rx944IENrz333HNJkiFDhmy0727dumWnnXbKunXr8vzzzzdTxQAAAFCGqt6N/5vf/GYOPvjgfP3rX8+OO+6Yo48+Oj179sxrr72W6667Lt/97nez884755xzzmnYZ/bs2UnWB/pN6d69e0P7v/7rv97imubOnbvJ7XV1dVvcFwAAALRFVQ37f/u3f5sHH3wwkyZNyvHHH5/Pfe5z6dq1a1asWJFOnTrl7//+7zN58uS84x3vaNhn6dKlTTpGU9tvasUAAAAAlKCqj95LkjfeeCO9e/dOTU1NkmTFihVJ1i/z79WrV954440N2i9fvrxJ/S9btqx5CgUAAIBCVHVm/7LLLsupp56anj175jvf+U6OOOKI7LrrrnnppZdy/fXX5//9v/+Xa6+9Nj/+8Y9zxBFHVLOUBnPmzNnk9rq6uuy9994tUgsAAABUQ9XC/osvvpjTTjsta9euzY033phDDjmkYds73/nOTJ48ObvuumtOOeWUHH/88ZkzZ0569Oix2Wv136r+2v0tNXjw4Ca1BwAAgPamasv4b7vttqxevTojR47cIOi/2UknnZTu3bvntddeyy9/+cskSc+ePZt0nKa2BwAAgNJVLezXP0Jvjz322Gibmpqahu317YcOHZrkL9f2b0z9tfr17QEAAID1qhb2165du/4AHTZ9iPob99W3Hz58eJJNX1u/fPnyzJ8/PzU1Ndl9992bo1wAAAAoRtXC/sCBA5P8Zca+MZVKJTNnzkyS7LrrrkmSffbZJ506dcoTTzyRJUuWNLrfb3/721QqlYwaNSo77LBDM1cOAAAA7VvVwv5BBx2UJPnTn/6UBx54oNE211xzTZYuXZoOHTrkgAMOSLL+GvyDDz44q1evzqWXXtroflOnTk2SjBkzpvkLBwAAgHauamH/gAMOyD777JNKpZJjjjkmV199dRYuXJgk+d///d9ceOGF+ed//uckyT/8wz9kt912a9j3zDPPTJKcc845ueSSS7J06dIk6x+Ld+qpp+anP/1pevXqlVNPPbVa5QMAAEC7VbVH7yXJT37yk3z0ox/N//zP/+Tzn/98kqRz585ZtWpVQ5t999033/ve9zbY7+CDD85ZZ52VCy64IKeddlq+9KUvpWvXrlm+fHlDHz/84Q+z0047VbN8AAAAaJeqNrOfrL9T/n/913/lW9/6Vj74wQ9mxx13zLp169KvX78cdNBBueyyy/Kb3/wmvXr1etu+559/fm677bYccsgh6devX9auXZtBgwbl2GOPzaOPPpojjjiimqUDAABAu1XVmf1k/TX4Z5xxRs4444wm73vUUUflqKOOav6iAAAAoGBVndkHAAAAWp6wDwAAAIUR9gEAAKAwwj4AAAAURtgHAACAwgj7AAAAUBhhHwAAAAoj7AMAAEBhhH0AAAAojLAPAAAAhalt7QIAAJJk2IQZrV0CABTDzD4AAAAURtgHAACAwgj7AAAAUBhhHwAAAAoj7AMAAEBhhH0AAAAojLAPAAAAhRH2AQAAoDDCPgAAABRG2AcAAIDCCPsAAABQGGEfAAAAClPb2gUAQL1hE2a0dgkAAEUwsw8AAACFEfYBAACgMMI+AAAAFEbYBwAAgMII+wAAAFAYYR8AAAAKI+wDAABAYYR9AAAAKIywDwAAAIUR9gEAAKAwwj4AAAAURtgHAACAwtS2dgEAQNs2bMKM1i4BAGgiM/sAAABQGGEfAAAACiPsAwAAQGGEfQAAACiMsA8AAACFEfYBAACgMMI+AAAAFEbYBwAAgMII+wAAAFAYYR8AAAAKI+wDAABAYYR9AAAAKIywDwAAAIUR9gEAAKAwwj4AAAAURtgHAACAwgj7AAAAUBhhHwAAAAoj7AMAAEBhhH0AAAAojLAPAAAAhRH2AQAAoDDCPgAAABRG2AcAAIDCCPsAAABQGGEfAAAACiPsAwAAQGGEfQAAACiMsA8AAACFadGwv3Tp0hx11FGpqanJgAED8uCDD26y/UMPPZRjjjkmu+66a7p06ZKBAwdm7Nixm90PAAAAtmctFvZff/31HHLIIbnjjjvy7ne/O48++mj222+/jba/4oorMnr06Nxyyy2pq6tLx44dM2/evNx666054IADMn369JYqHQAAANqVFgn769atyyc/+ck8/PDDef/7358HHnggQ4YM2Wj7Z599NqeeemrWrVuXs88+O/Pnz8+yZcuyYMGCTJw4MevWrctpp52Wp556qiXKBwAAgHalRcL+BRdckF/84hfZfffdc88996RPnz6bbH/RRRdl9erVGTduXM4777z069cvSdK3b99MmjQpJ598ctasWZNvf/vbLVE+AAAAtCtVD/vPPvtszj///HTs2DG33HJLQ3DflDvvvDNJcvrppze6ffz48UmSu+66K5VKpfmKBQAAgAJUPeyfc845WbVqVcaNG5dRo0Zttv2sWbNSV1eXPn36ZOTIkY222WOPPTJgwIDMnz/fUn4AAAB4i6qG/eeffz4333xzevfunfPOO2+L9nnuueeSZJPX9L95+8yZM7etSAAAAChMbTU7v+6667Ju3bp89rOfzdKlS/PlL3859913XxYsWJChQ4fmk5/8ZL7yla9khx12aNhn9uzZSZJu3bptsu/u3btv0H5LzZ07d5Pb6+rqmtQfAAAAtDVVDfu33nprkmSfffbJvvvuu0GQfvrpp3Peeefltttuy3333ddw076lS5c26RhNbb+5FQMAAADQ3lUt7L/++ut58skn06VLl9x000054YQTMm7cuAwaNCh1dXW58sorc8EFF+SJJ57ImWeemSuuuCJJsnz58iYdZ9myZdUoH4A3GTZhRmuXAABAE1Qt7D/xxBNZt25dVq5cmREjRuQb3/hGw7ahQ4c2XMM/efLk/PCHP8yFF164wXL+apkzZ84mt9fV1WXvvfeueh0AAABQLVUL+/Pnz2/4/pRTTmm0zSmnnJLJkydnxYoV+f3vf58DDjhgs9fqv1X9tftbavDgwU1qDwAAAO1N1e7Gv2TJkiRJz549s9tuuzXaZsCAAenXr1+Sv9w4r2fPnk06TlPbAwAAQOmqFvbrZ+h79OixyXb1YX3FihVJ1i/xf/PPG1N/rX59ewAAAGC9qoX9/v37J9n8DffeeOONJEnv3r2TJMOHD0+y+Wvr67fXtwcAAADWq1rYf+9735sOHTpkyZIl+fOf/9xom/nz5zdc2z9y5Mgk68P7wIEDs3DhwjzzzDON7jdr1qzMmzcv/fr1a9gPAAAAWK9qYb9Xr1456KCDkqThsXpvVf/6sGHD8q53vavh9SOPPDJJcvHFFze639SpU5MkRxxxRGpqapqrZAAAAChC1cJ+knzta19LklxwwQX51re+lUWLFiVJFi9enIsuuigTJ05MkkyYMGGD0D5+/Ph06tQp06dPz6RJk7Jw4cIkyWuvvZbJkydn2rRpqa2tzRlnnFHN8gEAAKBdqmrYP/jgg3P++ednzZo1+dd//df07ds3PXv2TJ8+ffLlL385q1evzoknnphx48ZtsN+IESNyySWXpEOHDjn33HPTr1+/9OjRI3379s0555yTJPnud79rCT8AAAA0oqphP0nOOuus/PrXv86YMWOy8847Z9WqVenXr18+9rGP5eabb873v//9Rpfijxs3Lvfee2/GjBmTAQMGZPXq1enfv3/GjBmTe++9NyeddFK1SwcAAIB2qbYlDnLQQQc1XL/fFKNHj87o0aOrUBEAAACUq+oz+wAAAEDLEvYBAACgMMI+AAAAFEbYBwAAgMII+wAAAFAYYR8AAAAKI+wDAABAYYR9AAAAKIywDwAAAIUR9gEAAKAwwj4AAAAURtgHAACAwgj7AAAAUBhhHwAAAApT29oFAACUZtiEGVU/xotTDq/6MQBov8zsAwAAQGGEfQAAACiMsA8AAACFEfYBAACgMMI+AAAAFEbYBwAAgMII+wAAAFAYYR8AAAAKI+wDAABAYYR9AAAAKIywDwAAAIUR9gEAAKAwwj4AAAAURtgHAACAwgj7AAAAUBhhHwAAAApT29oFsBVefTXZZZcNX3vllWTnnVunHgCgxQ2bMKNFjvPilMNb5DgANC8z+wAAAFAYYR8AAAAKI+wDAABAYYR9AAAAKIywDwAAAIUR9gEAAKAwwj4AAAAURtgHAACAwgj7AAAAUBhhHwAAAAoj7AMAAEBhhH0AAAAojLAPAAAAhRH2AQAAoDDCPgAAABRG2AcAAIDCCPsAAABQGGEfAAAACiPsAwAAQGGEfQAAACiMsA8AAACFEfYBAACgMMI+AAAAFEbYBwAAgMII+wAAAFAYYR8AAAAKI+wDAABAYYR9AAAAKIywDwAAAIUR9gEAAKAwwj4AAAAURtgHAACAwgj7AAAAUBhhHwAAAArTomF/wYIF6dWrV2pqalJTU5NJkyZtsv2TTz6Z448/PkOHDk2XLl2yyy675NBDD81dd93VMgUDAABAO9SiYX/KlCl5/fXXU1tbu9m2P/vZz/KBD3wg1157bebMmZOOHTvm1Vdfzc9//vMceeSROeuss1qgYgAAAGh/Wizs19XVZdq0aRk1alSOPfbYTbZ97bXXctxxx2XlypX5whe+kLlz52bZsmV5/fXXc8kll6Rz5875xje+kXvuuaeFqgcAAID2o8XC/vnnn5/ly5dn4sSJqamp2WTbyy67LIsWLcohhxySyy+/PIMGDUqS9OzZM6eeemomT56cZP1KAQAAAGBDm19P3wxeeumlXHnllRk1alTGjBmT22+/fZPt77jjjiTJ6aef3uj2U045Jeecc07uv//+LFq0KDvuuGMzVwzQfgybMKO1SwAAoI1pkZn9c889N6tWrdqiWf0VK1bkscceS5Lst99+jbbp2bNn3ve+92Xt2rV56KGHmr1eAAAAaM+qHvafffbZXHvttQ2z+pvz4osvZs2aNenVq1d69eq10XZDhgxJksycObPZagUAAIASVH0Z/8SJE7N27dotmtVPktmzZydJunXrtsl23bt336D9lpo7d+4mt9fV1TWpPwAAAGhrqhr2//jHP+bGG2/c4ln9JFm6dGmTjtHU9vUrAgAAAKBUVV3Gf/bZZ6dSqWzxrH6SLF++vEnHWLZs2daUBgAAAMWq2sz+I488kjvvvLNJs/otYc6cOZvcXldXl7333ruFqgEAAIDmV7Wwf9ZZZyVJk2b1k81fq/9W9dfub6nBgwc3qT0AAAC0N1VZxn/fffflV7/61VbN6vfs2bOq7QEAAKB0VQn7X//615M0fVY/SYYOHZokWbFixSbb1V+rX98eAAAAWK8qy/gffPDBJMkJJ5zQ6Pb6oD5lypRMnTo1SbJo0aIkybBhw9KxY8csXrw4S5YsSa9evRrto/7a++HDhzdf4QAAAFCAqt6Nf/HixY1+rV69OkmycuXKhtfqde3aNXvttVeS5Le//W2j/b7xxht5/PHH06FDh+y7777VfAsAAADQ7lQl7FcqlU1+HX/88UnWL/Ovf+3N6q/zv/jiixvtf/r06Vm5cmVGjx6dPn36VOMtAAAAQLtV1Zn9rTVu3Lj07t07P//5z3PSSSfl5ZdfTrJ+Rv/SSy/N1772tSTJV7/61dYsEwAAANqkNhn2+/btm+uvvz6dO3fOZZddlkGDBqVHjx7ZYYcdcsopp2TVqlWZMGFCDjnkkNYuFQAAANqcNhn2k+Twww/Po48+mmOPPTaDBg3KmjVr0rdv33z0ox/Nbbfdln/7t39r7RIBAACgTarK3fg35+qrr87VV1+92Xbvfe97c/3111e/IAAAAChIm53ZBwAAALaOsA8AAACFEfYBAACgMMI+AAAAFEbYBwAAgMII+wAAAFAYYR8AAAAKI+wDAABAYYR9AAAAKIywDwAAAIUR9gEAAKAwwj4AAAAURtgHAACAwgj7AAAAUBhhHwAAAAoj7AMAAEBhhH0AAAAojLAPAAAAhalt7QIAWsOwCTOqfowXpxxe9WMAVJvzJUD7ZGYfAAAACiPsAwAAQGGEfQAAACiMsA8AAACFcYM+gCppiZtaAQBAY8zsAwAAQGGEfQAAACiMsA8AAACFEfYBAACgMMI+AAAAFEbYBwAAgMII+wAAAFAYYR8AAAAKI+wDAABAYYR9AAAAKIywDwAAAIUR9gEAAKAwwj4AAAAURtgHAACAwgj7AAAAUBhhHwAAAAoj7AMAAEBhhH0AAAAojLAPAAAAhRH2AQAAoDC1rV0AAADbt2ETZrTIcV6ccniLHAegLTCzDwAAAIUR9gEAAKAwwj4AAAAURtgHAACAwgj7AAAAUBhhHwAAAAoj7AMAAEBhhH0AAAAojLAPAAAAhRH2AQAAoDDCPgAAABRG2AcAAIDCCPsAAABQGGEfAAAACiPsAwAAQGGEfQAAACiMsA8AAACFEfYBAACgMMI+AAAAFEbYBwAAgMII+wAAAFCYFgn7jz32WMaOHZtddtklnTp1ys4775zDDjssd9999yb3e+ihh3LMMcdk1113TZcuXTJw4MCMHTs2Dz74YEuUDQAAAO1S1cP+Oeeck7322iu33nprXn311XTo0CHz58/P3XffncMOOyzf/va3G93viiuuyOjRo3PLLbekrq4uHTt2zLx583LrrbfmgAMOyPTp06tdOgAAALRLVQ37V111VSZPnpzOnTvnm9/8ZubPn5+VK1fm+eefz5gxY5IkX/3qVzNr1qwN9nv22Wdz6qmnZt26dTn77LMzf/78LFu2LAsWLMjEiROzbt26nHbaaXnqqaeqWT4AAAC0S1UL+2vWrMmkSZOSJBdddFHOPPPM9OvXL0nyjne8IzfccEMGDBiQNWvW5IYbbthg34suuiirV6/OuHHjct555zXs17dv30yaNCknn3xy1qxZs9FVAQAAALA9q61Wx6tXr86pp56aO+64I+PGjXvb9m7duuXAAw/Mj3/84/zhD3/YYNudd96ZJDn99NMb7Xv8+PG59NJLc9ddd6VSqaSmpqb53wDQKoZNmNHaJQAAQLtXtZn9bt265cwzz8xDDz2U2trG/6YwePDgJMnixYsbXps1a1bq6urSp0+fjBw5stH99thjjwwYMCDz58+3lB8AAADeolUfvbdq1aokyY477tjw2nPPPZckGTJkyCb3rd8+c+bM6hQHAAAA7VTVlvFviaeffjrJ+mv4682ePTvJ+pUBm9K9e/cN2m+puXPnbnJ7XV1dk/oDAACAtqbVwv7LL7+ce++9N0nyiU98ouH1pUuXNqmfprbf3IoBAAAAaO9aLex/5StfyZo1a7L//vtn//33b3h9+fLlTepn2bJlzV0aAAAFaombwL445fCqHwNgS7RK2P/+97+fG264IT169MgVV1zRoseeM2fOJrfX1dVl7733bqFqAAAAoPm1eNi/5557cvLJJ6empibXXHNNRowYscH2zV2r/1b11+5vqfonAAAAAECpWjTsP/DAAzn66KOzevXqXHjhhRk7duzb2vTs2bNJfTa1PQAAAJSuxR6999hjj+UTn/hEli1blq9+9asZP358o+2GDh2aJFmxYsUm+6u/Vr++PQAAALBei4T9J598Mh//+MezZMmSnHTSSfnGN76x0bbDhw9Psvlr6+u317cHAAAA1qt62J85c2Y++tGPZsGCBfn0pz+dadOmbbL98OHDM3DgwCxcuDDPPPNMo21mzZqVefPmpV+/fhk5cmQ1ygYAAIB2q6phf/bs2fnwhz+cefPm5dBDD821116bDh02f8gjjzwySXLxxRc3un3q1KlJkiOOOCI1NTXNVi8AAACUoGphf+nSpfnIRz6S2bNnZ7/99sstt9ySTp06bdG+48ePT6dOnTJ9+vRMmjQpCxcuTJK89tprmTx5cqZNm5ba2tqcccYZ1SofAAAA2q2qhf358+fnueeeS5I8/vjjGThwYHbccceNfv3oRz9q2HfEiBG55JJL0qFDh5x77rnp169fevTokb59++acc85Jknz3u9+1hB8AAAAa0SKP3nvjjTc222bVqlUb/Dxu3Li8853vzIUXXpjf/e53WbBgQfr375999tkn48ePzwEHHFCtcgEAAKBdq1rYHzZsWCqVyjb1MXr06IwePbqZKgIAAIDtQ4s8eg8AAABoOcI+AAAAFEbYBwAAgMII+wAAAFAYYR8AAAAKI+wDAABAYYR9AAAAKIywDwAAAIUR9gEAAKAwwj4AAAAURtgHAACAwgj7AAAAUBhhHwAAAAoj7AMAAEBhalu7AAAAKMWwCTNa5DgvTjm8RY4DtF9m9gEAAKAwwj4AAAAURtgHAACAwgj7AAAAUBhhHwAAAAoj7AMAAEBhhH0AAAAoTG1rFwBsu5Z4pq/n+QIAQPthZh8AAAAKI+wDAABAYYR9AAAAKIywDwAAAIUR9gEAAKAwwj4AAAAURtgHAACAwgj7AAAAUJja1i4AaB+GTZjR2iUAAABbyMw+AAAAFEbYBwAAgMII+wAAAFAYYR8AAAAK4wZ9AADQzrTEjXNfnHJ41Y8BVI+ZfQAAACiMsA8AAACFEfYBAACgMMI+AAAAFEbYBwAAgMII+wAAAFAYYR8AAAAKI+wDAABAYYR9AAAAKIywDwAAAIUR9gEAAKAwwj4AAAAURtgHAACAwgj7AAAAUBhhHwAAAAoj7AMAAEBhhH0AAAAojLAPAAAAhRH2AQAAoDDCPgAAABSmtrULAAAA2p5hE2a0yHFenHJ4ixwHtjdm9gEAAKAwwj4AAAAURtgHAACAwgj7AAAAUBhhHwAAAAoj7AMAAEBhhH0AAAAoTJsP+w899FCOOeaY7LrrrunSpUsGDhyYsWPH5sEHH2zt0gAAAKBNqm3tAjbliiuuyEknnZR169YlSbp165Z58+bl1ltvze23355p06blpJNOauUqaa+GTZhR9WO8OOXwqh8DAKA9a4nPZInPZWx/2uzM/rPPPptTTz0169aty9lnn5358+dn2bJlWbBgQSZOnJh169bltNNOy1NPPdXapQIAAECb0mbD/kUXXZTVq1dn3LhxOe+889KvX78kSd++fTNp0qScfPLJWbNmTb797W+3cqUAAADQtrTZsH/nnXcmSU4//fRGt48fPz5Jctddd6VSqbRYXQAAANDWtclr9mfNmpW6urr06dMnI0eObLTNHnvskQEDBmTevHl56qmn8u53v7uFq4TNa6lr0AAA2H6UdO+pkt5LW9MmZ/afe+65JMmQIUM22a5++8yZM6teEwAAALQXbTLsz549O8n6u+9vSvfu3TdoDwAAALTRZfxLly6tWvu5c+ducvucOXMavq+rq2tSHS1mwYK3v/byy8nKlS1fSzu2Zsn81i4BAIAWsrkc0J60xOfYlvp9lfRetsWbs+eaNWuapc82GfaXL1/epPbLli3b4rabuzTgzfbee+8m1dGqRo1q7QoAAKDNGnJpa1fQvpT0+2pv7+XVV1/NsGHDtrmfNrmMHwAAANh6bXJmf3PX6r9V/bX7W+LNy/Qbs2LFijzzzDPp379/dt5559TWtslfUbOqq6trWMXwyCOPZODAga1cEdsrY5G2wDikLTAOaQuMQ9qC7WUcrlmzJq+++mqSZM8992yWPttkku3Zs2fV2g8ePHizbYYPH96k45dk4MCBW/Q7gmozFmkLjEPaAuOQtsA4pC0ofRw2x9L9N2uTy/iHDh2aZP0s+6bUX6tf3x4AAABoo2G/fmZ9c0vu67dvzzPxAAAA8FZtNuwPHDgwCxcuzDPPPNNom1mzZmXevHnp169fRo4c2cIVAgAAQNvVJsN+khx55JFJkosvvrjR7VOnTk2SHHHEEampqWmpsgAAAKDNa7Nhf/z48enUqVOmT5+eSZMmZeHChUmS1157LZMnT860adNSW1ubM844o5UrBQAAgLalzYb9ESNG5JJLLkmHDh1y7rnnpl+/funRo0f69u2bc845J0ny3e9+1xJ+AAAAeIs2G/aTZNy4cbn33nszZsyYDBgwIKtXr07//v0zZsyY3HvvvTnppJNau0QAAABoc2oqlUqltYsAAAAAmk+bntkHAAAAmk7YBwAAgMII+wAAAFAYYR8AAAAKI+wDAABAYYR9AAAAKIywDwAAAIUR9gEAAKAwwj4AAAAURtgv3IIFC9KrV6/U1NSkpqYmkyZNatL+V199dcO+W/L14osvVuV90L5t6zisV1dXlzPPPDOjRo3KDjvskE6dOmWXXXbJRz/60Vx99dVZt25d8xZOUZprHK5bty7XXHNNPvzhD2ennXZKly5dMmjQoPzDP/xD/uu//qt5i6bd+81vfrNF/35OmTKlyX1fe+21OfDAA7PTTjula9eu2X333XPyySdn9uzZVXgntHfVHItLly7NwQcfvEE/zzzzTBXeBe1dtcZhpVLJTTfdlNGjR6dXr17p0qVLdttttxx//PF56qmnqvRu2r7a1i6A6poyZUpef/311NbWZs2aNVvdT6dOndK9e/dGt1UqlSxZsiRJUlNTs9XHoFzNMQ4fe+yxfPSjH82iRYsaXuvSpUteffXV/OpXv8qvfvWr3HDDDbnrrrvSuXPnZqqckjTHOFy9enWOOeaY3HnnnUnWn/O6deuWl19+OT/+8Y9z44035sorr8znP//55iydAmzq39Ek6dq1a5P6+9znPpfrrrsuyfpx2KVLl7zwwguZPn16fvSjH+U//uM/8oEPfGCbaqZMzT0WFy9enEMPPTQPP/xwBg0alD//+c/bWiLbgeYchytXrsyRRx6Ze+65J8n6c2JtbW1mz56da6+9NjfffHPuvvvujB49epvrbm/M7Besrq4u06ZNy6hRo3LsscduU1/HHntsFi1a1OjX1VdfnSTZeeedM2jQoGaonJI01zg88cQTs2jRogwaNCh33HFH3njjjaxYsSJ1dXU566yzkiT33HNPpk2b1lylU5DmGocXXHBB7rzzzvTo0SNXXnllXn/99bzxxht55ZVXMmHChKxbty5f/OIX86c//akZq6cEm/p3dNGiRfm///f/bnFfV111Va677rr06tUr1113XZYtW5bly5dn5syZOfroo7NkyZJ88pOfzOrVq6v3hmi3mnMszp8/PwcddFAefvjh7LPPPvnpT39avcIpSnOOw3HjxuWee+7JoEGDcvPNN2fp0qVZtmxZHn744YwcOTLLli3LiSeeuF2uABX2C3b++edn+fLlmThxYlVn3OtnFj71qU+lttZiETbUHOPw+eefzxNPPJEkufTSS3PkkUc2/DV4wIABOf/883P00UcnSW6//fZmqZuyNMc4XLZsWaZOnZokueSSS/KP//iP6dGjR5L1f+z8t3/7t4wbNy6rV6/ON7/5zeYqHd6mfnxNnz49n/nMZxpmwP7qr/4qP/nJT/Le9743L774Ym644YbWLJPC1dXV5YADDsjjjz+eAw88ML/85S+z4447tnZZbGeeffbZXHfddenUqVNmzJiRsWPHpnv37qmtrc0+++zTkFNmzZqV3/72t61cbcsT9gv10ksv5corr8yoUaMyZsyYqh1n4cKFmTFjRpLkM5/5TNWOQ/vUXOPwf//3fxu+/9CHPtRom/rX582bt9XHoUzNNQ5/+9vfZvHixenevXuOO+64RtuMGzcuSXLHHXekUqls9bFgY5599tk8++yz6d+/fz71qU+9bXunTp3yz//8z0nWj0OolhkzZuSpp57KYYcdlrvvvjs9e/Zs7ZLYDnXv3j1f+tKX8tnPfjZ//dd//bbt73//+zNgwIAkyR/+8IeWLq/VmYYt1LnnnptVq1Zt86z+sGHDMnbs2Oy1116Nbr/xxhuzatWq7LHHHvnbv/3brT4OZWqucThw4MAtbrvrrrtu9XEoU3ONw6effjpJMnz48HTq1KnRNu9617uSJEuWLMkLL7yQ3XfffauPB4158MEHkyT77rtvOnRofM6m/rrU+++/v8XqYvvzT//0T+nWrVv+/u//fqPnRKi2IUOGNKy625jBgwdn3rx5Wbx4ccsU1YaY2S/Qs88+m2uvvbZZZvUPPPDA3HzzzTn11FMb3V6/NGZjs1xsv5pzHA4dOrQhND300EONtql//cADD9ymY1GW5hyH9Tci7dat20bbvPlmQ1aZUA3PPfdckvUfcDemftv8+fO3yw+3tJzjjjtO0KfNW7VqVZJsl5eZCPsFmjhxYtauXVv1a/XffO2LsM9bNec47NChQy688MJ07NgxJ598cu66664sX748yfol/l//+tdz2223Zfjw4U26oQvla85x2KVLlyTJG2+8sdE2S5cubfi+foxCvYULF+a0007LO9/5zgwcODD7779/rrjiiqxdu3aL+6h/rN6W/tHJY/hoTHOMRdhWLTEOV61aleeffz5J8o53vKPZ+m0vLOMvzB//+MfceOONVb9WP/nLrP7f/u3fZvjw4VU9Fu1LNcbhmDFjcs8992Ty5Mk58sgjk6wPXytXrkyPHj3yxS9+MZMnT07v3r2b5Xi0f809DutnS2fOnJkVK1Y0+ligJ598suH7Xr16bfMxKcfatWvzsY99LI899ljDa/PmzcuDDz6Yn/3sZ7nllls2uiz/zd78B6Ut0dT2lK+5xiJsi5Yah3fccUeWLl2aXr16efQe7d/ZZ5+dSqVS9Vn9JPnhD3+YxI35eLtqjcPly5dvEOZXrlyZJOncuXN69OixyRlXtj/NPQ733Xff1NTUZMWKFbnqqqsabfO9732v4fudd955m49JOe6666707ds3f/zjH7NixYrMmTMn55xzTmpqanL77bfn0ksv3aJ+mrpiZNmyZVtTLgVrrrEI26IlxuHy5cszYcKEJMnpp5/e8ASd7YmwX5BHHnkkd955Z4vM6v/2t7/NrFmzUltbm09/+tNVPRbtS7XG4dlnn51PfOITeeSRR/KDH/wgL774Yt5444388Y9/zKc//elMnTo1H/jAB/LII4802zFpv6oxDocMGdKwqmT8+PG58MILM2/evKxZsyYzZ87MKaec0vAHqD59+mS33XZrluPSvg0bNixf/OIXc/TRR+f222/PnnvumS5dumTw4ME599xzc/zxxydJLr744laulNIZi7QFLTkOv/jFL+b555/Pe97znnzta1/b5v7aI2G/IGeddVaStMis/rXXXpsk+djHPpaddtqpqseifanGOHzooYdy/vnnp3Pnzrnnnnvy+c9/Prvttlu6d++ePffcM9/73vfyla98JQsWLMjnPvc5jzyjaufD733vexk+fHhWrlyZf/mXf8nAgQPTqVOn7LHHHvnNb36Tgw46KEnywQ9+sOrnYdqHYcOGZfr06bnqqqsavc7+i1/8YpL1N97785//vNn+NnWtfmPefP0+27fmHouwNVpqHE6cODHXXXdd+vbtm9tuu63RS++2B8J+Ie6777786le/apFZ/VWrVuXGG29MYgk/G6rWOKwfbx/5yEfynve8p9E2X/7yl5Osv/v69vgcVf6imufDXXfdNY899ljOPvvs/PVf/3V22WWXvOc978mXv/zl3H///fnFL36RJDnmmGOa9biU6833vHnppZc2276pzzL37HO2VFPHIlRDc4zDCy+8MOedd166deuWO++8c7u+t5gb9BXi61//epKWmdX/6U9/mtdeey077LBD1f+wQPtSrXFY/6ipPfbYY6Nt+vfvn969e2fx4sV57rnnMmrUqGY7Pu1Ltc+HvXv3znnnnZfzzjtvg9dfeumlzJgxIz179szYsWOb/biU6c2r41asWLHZ9kOHDt1s2zdfp1/fHjanqWMRqmFbx+Hll1+ef/mXf0ltbW1uuummfOhDH2rO8todYb8QDz74YJLkhBNOaHR7/T/8U6ZMydSpU5MkixYt2qpj1d+F/+ijj27yckLKVq1xWP8Ils3dlbU+2Hl00PatJc+HbzZp0qSsXr0648ePdyd+ttiCBQsavt+Sp4nUz1DNmTNno23qt+20006eUMIWa+pYhGrYlnF4/fXX5+STT05NTU2uueaaHH744c1dXrsj7Bdm8eLFm9y+cuXKhhtIbY0FCxbkZz/7WRJL+Nm45h6HAwcOTPKXGf7GvPrqqw2Bbdddd93ivilXtc+Hb/bLX/4yV199dfr165czzjijWfqkDCtXrswNN9yw0T8+zZo1K8n6P2b+n//zfzbb3/77758kefjhh1OpVBpdvVL/B6/6tpA0/1iErVGtcXj77bfnhBNOyLp16zJt2rQce+yxzVFuu+ea/UJUKpVNftXf2XLixIkNr22Nn/zkJ1m1alUGDhyYD3/4w835FihAtcZh/U3PfvWrX2008E+bNi3J+utT99prr2Z4N7RXLXU+rPfCCy80fKj4zne+46alNFi7dm0OOuigfP7zn89NN93UaJvLLrssSfKhD30oO+yww2b7HDFiREaMGJG6urqG+5m82Zo1a3LJJZckiUvtaFCNsQhNVa1xeM899+TTn/501qxZk8mTJ+eUU05ptprbO2GfrF27Nh/5yEfStWvXfOQjH9nkEuj6Jfz/8A//sNkl1dAUmxqHn/70p7PbbrtlxYoV+fjHP57bb789S5YsSaVSydy5c3PWWWflggsuSLL+OaouL2FrNeV8mCSzZ8/OwQcfnPnz5+fEE09s+EMCJEnHjh0bZtePO+64TJkypWGJ6muvvZZvfOMbufrqq5Nkg8dCLV68OO9///vTrVu3fO5zn3tbv2eeeWaS9Xet/tGPftSwQuX555/Ppz/96fz3f/93dtttt/zDP/xDNd8e7Ui1xiI0RTXG4X//93/n6KOPzsqVK/PlL3+54b49/P8qbBeOP/74SpLKxIkT37Ztzpw5lSQNX3Pnzm20j+eee66hze9///sqV0yJtmUcPv7445UBAwZs0KZz584b/HzUUUdVVq9e3ULvhvaqOc6HlUql8sc//rEyaNCgSpLKJz7xicrKlSurWDXt1Zo1ayqnnnrqBuOqW7duG/x83nnnbbDPAw880LCttra20X6PO+64hjYdOnSodO3ateHnXr16VR555JGWeHu0I9UYix07dnzbV337t75+7rnnttRbpQ1r7nF41VVXNWzr3bv3Zr9eeumllny7rc7ULBk4cGAOPvjgdO7cOQcffHAGDBjQaLv6Wf2RI0fmfe97X0uWyHZgc+Nw1KhR+eMf/5izzz4773vf+7LDDjtk3bp12WWXXXLooYfmxz/+cW699dbU1roVCVtvS8+Hd911Vz70oQ/lz3/+cz7/+c/nlltuSefOnVu4WtqDjh075pJLLsm9996bo446KkOHDs3atWvTv3//HHXUUbn33ntz9tlnb7DPnnvumVGjRqVLly4bnZ3/4Q9/mO9///vZb7/90rt371Qqley22275whe+kD/84Q8uZ+JtqjEW165d+7avjW1bt25d1d8jbV+1zonJ+hUAm/va3sZhTaWyjRcrAgAAAG2KmX0AAAAojLAPAAAAhRH2AQAAoDDCPgAAABRG2AcAAIDCCPsAAABQGGEfAAAACiPsAwAAQGGEfQAAACiMsA8AAACFEfYBAACgMMI+AAAAFEbYBwAAgMII+wAAAFAYYR8AAAAKI+wDAACw3fjXf/3X1NTUNHxNnz696sf8n//5n7zrXe9KTU1N/vZv/zZ//vOfq37M2qofAQAAAFpZpVLJP//zP+d73/te+vTpkxUrVmT58uVVP+5TTz2Vgw46KK+88ko++clP5rrrrkuXLl2qflwz+wAAABRt7dq1OfHEE/O9730vu+yyS+69997ssssuVT/u//7v/+ZjH/tYXnnllfzjP/5jfvKTn7RI0E+EfQDgTX760582LGv80Y9+tMm2v/jFLzZYBrl27dpNtv+7v/u71NTUpFOnTpkzZ05zlg0AG7V69eoce+yxufrqqzN48ODcf//9+eu//usWOfZxxx2XuXPn5uijj87ll1+empqaFjluIuwDAG/Su3fvhu+XLFmyybb//u//vsHPm2o/e/bs3HnnnUmSo446KkOGDNmGKgFgyz3zzDO56667svvuu+eBBx7IiBEjWuS4V111Vf7jP/4ju+22W6699tp06NCy8VvYBwAabGnYnzlzZu6+++4NXttU+0svvbRh5v9LX/rSNlYJAFtuzz33zO23354HHnggw4YNa/L+r7/+es4///yMGjUqPXr0yA477JC99947l1xySVavXt3oPqtWrco555yTJPnWt76Vnj17bstb2Cpu0AcANNjSsH/JJZekUqlkxx13zKJFi5IkixcvbrTtypUrc+WVVyZJRo0alf3337/5CgaALXDIIYds1X7z5s3LwQcfnKeffjpJ0rFjx1QqlTz66KN59NFHM2PGjNx1112prd0wWv/oRz/K3Llz86EPfSh///d/v831bw0z+wBAgy0J+0uXLs1VV12VJDnttNM22/7HP/5x5s+f/7b2ANDWnXjiiXn66aez11575be//W2WLVuWpUuX5mc/+1mGDBmSn//857nooovett8111yTJBk/fnweeOCBHHLIIenfv3923HHH7LvvvvnBD36QSqVS1dqFfQCgQa9evRpuHrSx8H7NNddkyZIl6dy5c0499dR07959k+0vueSSJMlOO+2UY489tgpVA0Dz+6//+q/cfffd2XnnnXPPPffkgx/8YDp37pxu3brl0EMPbQj0U6dO3eAmtfPnz8/999+fHj16pEOHDjn44IPzy1/+Mq+88koWL16chx9+OP/4j/+Yz33uc1WtX9gHABp06NAhO+ywQ5LGw3ulUmkI75/+9KfTv3//9OrVK0njy/j/8z//M4899liS5Atf+EK6du1ardIBoFndcsstSZK///u/z4477vi27QcddFAGDBiQl19+OY8//njD67/73e+ybt26vOc978nEiRNz2WWX5dVXX83KlSvzhz/8IWPHjk2S/PCHP8z1119ftfqFfQBgA/VL+RsL+7/85S/zzDPPJElOP/30JGkI+421r//DQG1tbU4++eSq1AsA1VD/x+q/+qu/2mib3XbbLUny3//93w2v/eEPf0iyPvR/85vfzIknnpiddtopnTt3znvf+97cdNNNGT16dJLke9/7XpWqF/YBgLeoD/uvv/7627Z997vfTZJ86EMfyvvf//4k2ehKgFdeeSU33XRTEo/bA6D9+d///d8kyVlnnZUdd9yx0a/6Pwi88sorDfvV36fmHe94Rw499NC39VtTU5NTTjklyfo/KGzsjv7byt34AYAN1C9VfGt4nzVrVn72s58l+cusfpKNLuO//PLLs2rVqiQetwdA+/PGG28kSZYvX57ly5dvsu3SpUsbvq//9/Pd7373RtvXb1u1alVeeeWVDBo0aFvLfRthHwDYwMaW8U+bNi3r1q3LkCFDcvTRRze83tgy/jVr1uSyyy5L4nF7ALRPPXr0SJJceumlOemkk7Z4v27dum2wf2N69uzZ8P2KFSu2ssJNs4wfANhAY2H/jTfeyA9+8IMkyamnnrrB84QbW8Z/++23Z+7cuUk8bg+A9ql///5JstlZ/a3Zr37VQLLhY2+bk7APAGyg/kPHG2+8kXXr1iVJrr322ixevDjdunXLF77whQ3aN7aMv/7afo/bA6C9+sAHPpAkDTembcwLL7zwtpn5973vfUmSp556aqP7Pf3000mSXXbZJTvttNO2ltooYR8A2ED9NfuVSqXhJn314f0zn/lM+vbtu0H7ty7jf+KJJ3L//fcn8bg9ANqv+kfk3XLLLVm4cOHbtq9atSqf+tSnstdee2XZsmUNrx988MHZYYcdMnPmzPzmN79ptO8rrrgiSfLxj3+8+Qv//wn7AMAG3ryccMmSJfnVr37VMAPR2I323hr2m/K4vV/84hc59NBDM3To0HTp0iX9+/fPAQcckGuuuaZZ3gsAbK2/+Zu/yWGHHZYFCxbkYx/7WH73u9+lUqlk7dq1eeyxx3LooYfm0Ucfzbve9a507969Yb9u3bpl/PjxSZJPfepTufXWWxtuWDtnzpz84z/+Y+6555506tQpX/nKV6pWv7APAGzgrWH/3//935MkH/7wh/Oe97znbe3rr9lfvHhxFi1alOuvvz7J5h+3d8YZZ+TjH/94Zs2alSOOOCLjx4/Pxz72sbzwwgt5/PHHm/MtAbCd+/CHP5za2toNvl566aUkf7kXTf3Xhz/84Yb9rrzyyrzrXe/KY489ln322Sddu3ZNly5dstdee+XXv/51PvjBD+byyy9/2/HOOuusHH744XnllVcyduzYdO/ePT169MjQoUPzgx/8IB06dMj06dOz5557Vu09uxs/ALCBN4f9P/zhD5kxY0aSjT8+780z+1dddVXDTYc29bi9559/Pt/5zndy1FFH5ZZbbkmHDn+Zf1i3bl0WLVq0rW8DABqsXbs2a9eubXRb/f1p3ty23sCBA/PII49k6tSpuemmmzJz5sx07949w4cPz+c+97mNXq7WuXPn3Hnnnbnqqqty9dVX54knnsjy5cszePDgjB49Ov/yL/+S97///c37Jt+iplKpVKp6BACgXbn77rtz2GGHJUn23HPPPPHEE9l9993z3HPPbRDK691+++05+uij07Vr1wwePDgzZ87MqFGjNjk7P2PGjHziE5/I1772tVxwwQVVey8AsL2yjB8A2MCbZ/afeOKJJOsfn9dY0E/+sox/xYoVmTlzZkP7TRk1alS6d++eKVOm5Ljjjsudd965wc2NAIBtI+wDABt46/N+d9hhh5x44okbbV+/jL/eljxub9CgQbn//vtz2GGH5eabb86YMWOy00475fjjj09dXd3WFw8AJBH2AYC3eGvYP/74498W6N/srdu29HF7f/M3f5O77rorCxcuzO23354PfOADufbaa/PZz3526woHABq4Zh8AaBPWrl2bnXbaKd26dcvLL7/c2uUAQLtmZh8AaFG///3v8+KLL77t9fvvvz+LFi3KPvvs0/JFAUBhzOwDAC3qhBNOyLXXXpv99tsvf/M3f5MuXbrkqaeeys9+9rP0798/DzzwQHbffffWLhMA2jVhHwBoUbfffntuuummPPLII6mrq8vq1aszbNiwHHHEEZkwYUJ22mmn1i4RANo9YR8AAAAK45p9AAAAKIywDwAAAIUR9gEAAKAwwj4AAAAURtgHAACAwgj7AAAAUBhhHwAAAAoj7AMAAEBhhH0AAAAojLAPAAAAhRH2AQAAoDDCPgAAABRG2AcAAIDCCPsAAABQGGEfAAAACiPsAwAAQGGEfQAAACiMsA8AAACFEfYBAACgMMI+AAAAFEbYBwAAgMII+wAAAFCY/w+C4ojT6toTtwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 600x600 with 2 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 525,
       "width": 509
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f, ax = plt.subplots(2, 1)\n",
    "ax[0].hist(statSREMC[:, 0], bins=50)\n",
    "ax[0].vlines(stat_obs[0], 0, 50, color=\"red\")\n",
    "ax[0].set_xlabel(r\"$\\hat{\\tau}_s$\")\n",
    "ax[1].hist(statSREMC[:, 1], bins=30)\n",
    "ax[1].vlines(stat_obs[1], 0, 50, color=\"red\")\n",
    "ax[1].set_xlabel(r\"$W_s$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Neyman"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def neyman_SRE(z, y, x):\n",
    "    xlevels = np.unique(x)\n",
    "    K = len(xlevels)\n",
    "    PiK, TauK, varK = np.zeros(K), np.zeros(K), np.zeros(K)\n",
    "    for k in range(K):\n",
    "        id = np.where(x == xlevels[k])\n",
    "        zk, yk = z[id], y[id]\n",
    "        PiK[k] = zk.shape[0] / z.shape[0]\n",
    "        TauK[k] = np.mean(yk[zk == 1]) - np.mean(yk[zk == 0])\n",
    "        varK[k] = np.var(yk[zk == 1], ddof=1) / sum(zk) + np.var(\n",
    "            yk[zk == 0], ddof=1\n",
    "        ) / sum(1 - zk)\n",
    "    return np.sum(PiK * TauK), np.sum(PiK**2 * varK)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sim_cluster(K, n, n1, n0):\n",
    "    x = np.repeat(range(K), n)\n",
    "    y0 = np.random.exponential(1, n * K)\n",
    "    y1 = y0 + 1\n",
    "    # block level assignment vector\n",
    "    zb = np.repeat([0, 1], [n0, n1])\n",
    "    MC = int(1e4)\n",
    "    TauHat, VarHat = np.zeros(MC), np.zeros(MC)\n",
    "    for k in range(MC):\n",
    "        z = np.concatenate([np.random.permutation(zb) for i in range(K)])\n",
    "        y = z * y1 + (1 - z) * y0\n",
    "        TauHat[k], VarHat[k] = neyman_SRE(z, y, x)\n",
    "    plt.hist(TauHat, bins=50)\n",
    "    plt.vlines(1, 0, 500, color=\"red\")\n",
    "    return np.var(TauHat), np.mean(VarHat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.009649300175537943, 0.00978208282388499)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABAoAAAPyCAYAAAD1wj3WAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAB7CAAAewgFu0HU+AABKb0lEQVR4nO39e5yXdYH//z8HBjmIEAc5iUouxopa1i1dP2ZqtllpSi3qlm5mtiHKmtFmUYagdGC7uambB5RsLcs+n81TmunWbmRq/jL7ZGqIAYqIDp8EPBFn5vr9wXdeCzmcnPd7Bpj7/Xab2w3mer2u9+vt5Qzzfsx1Xe+GqqqqAAAAACTp0tELAAAAAHYcQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFI0dvYBdyapVq/LYY48lSfbcc880NvrPCwAAQP2sW7cuL7zwQpLk4IMPTo8ePdq8T69ka+ixxx7LYYcd1tHLAAAAoBN66KGHcuihh7Z5Py49AAAAAApnFNTQnnvuWf780EMPZejQoR24GgAAAHZ1TU1N5cz2jV+TtoVQUEMb35Ng6NChGT58eAeuBgAAgM6kVvfJc+kBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQtGsoWL58eT74wQ+moaEhQ4YMyf3337/F8Q888EBOPvnkDBs2LN27d8/QoUMzduzYrc5r61wAAADorBqqqqra44FeffXVvPe9782DDz6YAw88MHfffXf23nvvzY6fOXNmxo8fn+bm5iRJz549s3LlyiRJly5dctVVV2X8+PE1n9sWixYtKs/p2WefzfDhw2v+GAAAANCiHq9D2+WMgubm5pxyyil58MEH87a3vS333XffFiPBk08+mQkTJqS5uTmTJ0/OkiVLsmLFiixdujRTpkxJc3NzzjvvvMyePbumcwEAAKCza5dQ8JWvfCX/+Z//mf322y8//elP069fvy2Ov+yyy7J27dqMGzcul1xySQYMGJAk6d+/f6ZOnZpzzjkn69aty6WXXlrTuQAAANDZ1f3SgyeffDJvfvObs379+jz88MM55JBDtjpn2LBhaWpqyh/+8IeMHj36Ndvnzp2bN73pTRk4cGD+9Kc/paGhoSZz28qlBwAAALSnnfLSg4suuihr1qzJuHHjtikSzJ8/P01NTenXr1+rL/STZP/998+QIUOyZMmSTS4haMtcAAAAoM6h4KmnnsrNN9+cvn375pJLLtmmOXPnzk2SLd7DYOPt8+bNq8lcAAAAIGms585vvPHGNDc356Mf/WiWL1+ez3zmM7n33nuzdOnS7LPPPjnllFPy2c9+NnvssUeZs3DhwiQb3qlgS3r16rXJ+LbO3RaLFi3a4vampqbt2h8AAADsaOoaCm699dYkyeGHH54jjjhikxfSTzzxRC655JLcdtttuffee8sNDpcvX75dj7Hx+LbM3RZbO1MBAAAAdnZ1u/Tg1VdfzeOPP57u3bvnhz/8Yc4888w8/fTTWbNmTZ555plMnjw5Xbp0yWOPPZbPfe5zZd7KlSu363FWrFhRk7kAAABAHc8oeOyxx9Lc3JzVq1dn1KhR+epXv1q27bPPPuWeBdOmTcv3vve9fOMb39jkEoQd0bPPPrvF7U1NTTnssMPaaTUAAABQe3ULBUuWLCl/Pvfcc1sdc+6552batGlZtWpV/u///b85+uijt3p/gb/Ucr+BZOv3JtjS3G3h7Q4BAADY1dXt0oNXXnklSdK7d+/su+++rY4ZMmRIBgwYkOR/bhTYu3fv7Xqcjce3ZS4AAABQx1DQ8tv93XfffYvjWl6sr1q1KsmGyxI2/vvmtNxfoGV8W+cCAAAAdQwFgwcPTrL1Gwz++c9/TpL07ds3STJy5MgkW78fQMv2lvFtnQsAAADUMRS8+c1vTpcuXfLKK6/kueeea3XMkiVLyr0MRo8enWTDi/ehQ4dm2bJlmTNnTqvz5s+fn8WLF2fAgAFlXlvnAgAAAHUMBX369Mm73vWuJMnMmTNbHdPy+REjRuSAAw4onz/ppJOSJFdccUWr8y6//PIkyYknnpiGhoZNtrVlLgAAAHR2dQsFSfLFL34xSfKVr3wlX//61/PSSy8lSV5++eVcdtllmTJlSpJk0qRJm7xonzhxYrp165YZM2Zk6tSpWbZsWZLkxRdfzLRp03LVVVelsbExF1xwwWsesy1zAQAAoLNrqKqqqucDfOUrX8mXvvSlDQ/W0JBevXplxYoVaXnYs846K9/61rde89v96667Luecc06am5uTpMxr2c/VV1+d8ePHt/qYbZnbFosWLcree++dZMN9ELydIgAAAPVUj9ehdT2jIEkuvPDC/PznP8+YMWOy5557Zs2aNRkwYEDe+9735uabb87111/f6iUA48aNy6xZszJmzJgMGTIka9euzeDBgzNmzJjMmjVriy/02zIXAAAAOrO6n1HQmTijAAAAgPa0U55RAAAAAOw8hAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAICisaMXAABQKyMm3dXmfSyYfkINVgIAOy+hAABgByR6ANBRXHoAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAABFY0cvAABgVzNi0l0dvQQAeN2cUQAAAAAUQgEAAABQuPQAAGAjLhsAoLNzRgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFI0dvQAAgBYjJt3V0UsAgE7PGQUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUjR29AAAA6mPEpLvavI8F00+owUoA2Jk4owAAAAAohAIAAACgEAoAAACAwj0KAICaqMX18ABAx3NGAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAEVdQ8EvfvGLNDQ0bPVj+vTprc5/4IEHcvLJJ2fYsGHp3r17hg4dmrFjx+b+++/f6mO3ZS4AAAB0Vu1yRkG3bt3St2/fzX706NHjNXNmzpyZo446KrfcckuamprStWvXLF68OLfeemuOPvrozJgxY7OP15a5AAAA0Jm1Syg47bTT8tJLL23249Of/vQm45988slMmDAhzc3NmTx5cpYsWZIVK1Zk6dKlmTJlSpqbm3Peeedl9uzZr3mstswFAACAzm6HvEfBZZddlrVr12bcuHG55JJLMmDAgCRJ//79M3Xq1JxzzjlZt25dLr300prOBQAAgM5uhwwFd9xxR5Lk/PPPb3X7xIkTkyR33nlnqqqq2VwAAADo7Ha4UDB//vw0NTWlX79+GT16dKtj9t9//wwZMiRLlizZ5BKCtswFAAAAdsBQMHfu3CTJ3nvvvcVxLdvnzZtXk7kAAABA0theD7Rs2bJMmTIlP/vZz/Lyyy9n5MiROeOMM3LWWWela9euZdzChQuTJD179tzi/nr16rXJ+LbO3RaLFi3a4vampqbt2h8AAADsaNolFKxfvz7vfe978/DDD5fPLV68OPfff39+8pOf5JZbbkmXLhtObli+fPl27Xvj8W2Zuy22dqYCAAAA7OzaJRTceeed+Zu/+Zs8+uijedOb3pQXXnghM2fOzLRp03L77bfnmmuuyYQJE5IkK1eu3K59r1ixovy5LXMBALZX/xUv5/9+8/RNPve2876fZb36dtCKam/EpLvavI8F00+owUoAaC91DQUjRozI2WefndWrV+fqq68ulwQMHz48F198cRYuXJgbbrghV1xxRQkFO7Jnn312i9ubmppy2GGHtdNqAAAAoPbqHgpmzJix2e1nn312brjhhsydOzfPPfdc9tprr63eX+AvtdxvINn6vQm2NHdbDB8+fLvGAwAAwM6mQ9/1YOTIkeXPzzzzTJKkd+/e27WPjce3ZS4AAADQwaFg4MCB5c+rVq1Kkuyzzz6b/H1zWu4v0DK+rXMBAACADg4FS5cuLX/u23fDTX9azjLY2v0AWrZvfFZCW+YCAAAAdQ4Fq1evzg033LDZ7fPnz9+wiC5d8qY3vSnJhhfvQ4cOzbJlyzJnzpzNzlu8eHEGDBiQ0aNHl8+3ZS4AAABQx1Cwfv36vOtd78rHP/7x/PCHP2x1zLXXXpskecc73pE99tijfP6kk05KklxxxRWtzrv88suTJCeeeGIaGho22daWuQAAANDZ1e1dD7p27Zp3vvOdefDBB3P66adn/vz5+eQnP5kBAwbkxRdfzDXXXFPONvjiF7+4ydyJEyfm29/+dmbMmJHBgwfnU5/6VPr3758XX3wxV155Za666qo0NjbmggsueM3jtmUuAAC1N2LSXW3ex4LpJ9RgJQBsi7peevDVr341EyZMyNq1a/OFL3whAwcOTK9evdK/f/9ceOGFaW5uziWXXJL3ve99m8wbNWpUrrzyynTp0iUXX3xxBgwYkN133z39+/fPRRddlCT55je/2eqlA22ZCwAAAJ1dXUNB165dc+WVV2bWrFn54Ac/mH322Sfr16/P4MGD88EPfjCzZs3K5MmTW507bty4zJo1K2PGjMmQIUOydu3aDB48OGPGjMmsWbMyfvz4zT5uW+YCAABAZ1a3Sw82dswxx+SYY47Z7nlHHXVUjjrqqNf1mG2ZCwAAAJ1Vh749IgAAALBjEQoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAorGjFwAAdLwRk+7q6CUAADsIZxQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUjR29AAAA2JoRk+5q8z4WTD+hBisB2PU5owAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAIp2DQVLly5Nnz590tDQkIaGhkydOnWL4x9//PF87GMfyz777JPu3btn0KBBef/7358777xzq4/VlrkAAADQWbVrKJg+fXpeffXVNDY2bnXsT37yk7z97W/Pd7/73Tz77LPp2rVrXnjhhdxzzz056aSTcuGFF9ZlLgAAAHRm7RYKmpqactVVV+WQQw7JaaedtsWxL774Yk4//fSsXr06n/zkJ7No0aKsWLEir776aq688srstttu+epXv5qf/vSnNZ0LAAAAnV27hYIvf/nLWblyZaZMmZKGhoYtjr322mvz0ksv5bjjjst1112XvfbaK0nSu3fvTJgwIdOmTUuy4QyFWs4FAACAzq5dQsEzzzyTb33rWznkkEMyZsyYrY7/0Y9+lCQ5//zzW91+7rnnpnv37vnlL3+Zl156qWZzAQAAoLNrl1Bw8cUXZ82aNdt0NsGqVavy8MMPJ0mOPPLIVsf07t07b33rW7N+/fo88MADNZkLAAAAtEMoePLJJ/Pd7353m88mWLBgQdatW5c+ffqkT58+mx239957J0nmzZtXk7kAAABAsvW3H2ijKVOmZP369dt0NkGSLFy4MEnSs2fPLY7r1avXJuPbOndbLFq0aIvbm5qatmt/AFALIybd1dFLAAB2IXUNBY8++mj+4z/+Y5vPJkiS5cuXb9djbDy+LXO3RcuZCAAAALCrquulB5MnT05VVdt8NkGSrFy5crseY8WKFTWZCwAAANTxjIKHHnood9xxx3adTbCje/bZZ7e4vampKYcddlg7rQYAAABqr26h4MILL0yS7TqbINn6/QX+Usv9Bto6d1sMHz58u8YDAADAzqYulx7ce++9+a//+q/XdTZB7969X/f4tswFAAAA6hQKvvSlLyXZ/rMJkmSfffZJkqxatWqL41ruL9Ayvq1zAQAAgDpdenD//fcnSc4888xWt7e8UJ8+fXouv/zyJMlLL72UJBkxYkS6du2al19+Oa+88kr69OnT6j5a7hcwcuTI8rm2zAUAAADq/K4HL7/8cqsfa9euTZKsXr26fK5Fjx49cuihhyZJfvWrX7W63z//+c/53e9+ly5duuSII46oyVwAAACgTqGgqqotfnzsYx9LsuHShJbPbazlvgZXXHFFq/ufMWNGVq9enaOOOir9+vWr2VwAAADo7Op6RsHrNW7cuPTt2zf33HNPxo8fn+effz7JhrMBrrnmmnzxi19MknzhC1+o6VwAAADo7HbIUNC/f/98//vfz2677ZZrr702e+21V3bffffsscceOffcc7NmzZpMmjQpxx13XE3nAgAAQGe3Q4aCJDnhhBPym9/8Jqeddlr22muvrFu3Lv3798973vOe3Hbbbfna175Wl7kAAADQmdXlXQ+25oYbbsgNN9yw1XFvfvOb8/3vf/91PUZb5gIAAEBntcOeUQAAAAC0P6EAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgaOzoBQAAQHsYMemuNu9jwfQTarASgB2bMwoAAACAQigAAAAACpceAEAHqsWp0AAAteSMAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAACKxo5eAAAA7CxGTLqrzftYMP2EGqwEoH6cUQAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAEVjRy8AAHZWIybd1dFLAACoOWcUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAABFXUPBq6++mn/913/NW97ylvTs2TM9evTI/vvvn/Hjx2fu3LlbnPvjH/8473//+zN48OB07949e++9d84444w8/vjjW33ctswFAACAzqxuoeCpp57K29/+9nz2s5/No48+mvXr12ft2rWZN29err322hxyyCG5//77W507efLknHjiibnnnnvypz/9KV27ds2iRYty44035u1vf3vuvPPOzT5uW+YCAABAZ1eXULBy5cp84AMfyB//+Me8613vym9/+9usWrUqq1atys9//vOMHDkyK1asyCc+8YnXzP35z3+eL3/5y+nWrVuuuOKKvPLKK1mxYkWee+65jB8/PqtXr84//MM/ZMmSJTWdCwAAANQpFDQ2Nuakk07KYYcdlnvuuSdve9vb0qVLl3Tr1i3vete7MnPmzCTJH//4xzzxxBObzP3617+eJLn44ovzqU99KnvssUeSZNiwYbnmmmty/PHH55VXXsnVV1/9msdty1wAAACgTqGgW7dumT59eu67777stttur9l+6KGHlj83NTWVPy9fvjw///nP09jYmHPPPbfVfX/6059OkvzoRz/a5PNtmQsAAABsUNebGbYWCZJkzZo15c8DBw4sf/71r3+dtWvX5uCDD07fvn1bnXvkkUemoaEhjzzySF599dWazAUAAAA26JC3R/zJT36SJBk9enQOOuig8vmWd0LYe++9Nzu3Z8+eGThwYJqbm/PUU0/VZC4AAACwQWN7PthLL72U2267Lf/8z/+c3r1759prr02XLv/TKhYuXJhkwwv6LenVq1cZ/5a3vKXNc7fVokWLtrh948soAAAAYGdU91Bw00035dxzz8369euzfPny7Lbbbjn99NNzwQUX5IADDthk7PLly7dr3xuPb8vcbbWlsxUAAABgV1D3Sw/WrFmTl19+ubwwb25uzoIFC17zbgfJhrdV3B4rVqyoyVwAAABgg7qHgjPPPDNVVZVAcPXVV+fRRx/N2LFjM3369Ho/fE09++yzW/x46KGHOnqJAAAA0Cbtdo+ChoaG7LvvvvnkJz+Z/fbbL3/7t3+bL33pS/nIRz6SfffdN8nW7y/wl1ruN9DWudtq+PDh2z0HAAAAdiYd8q4H7373u/PXf/3XWb9+fW677bby+d69e2/XfjYe35a5AAAAwAYdEgqSZOTIkUn+520Nk2SfffZJkqxatWqLc1vuL9Ayvq1zAQAAgA06LBS0vC1ic3Nz+VxLPHj22Wc3O2/lypVZsmRJGhoast9++9VkLgAAALBB3ULB//7f/ztVVW12+/z585Mkw4YNK587/PDD061btzz22GN55ZVXWp33q1/9KlVV5ZBDDskee+xRk7kAAADABnUJBR//+MfzkY98JJdddlmr23/1q19l9uzZSZJjjz22fL5379459thjs3bt2lxzzTWtzr388suTJGPGjNnk822ZCwAAAGxQl1Bw+OGHJ0k++9nP5otf/GKef/75JMkrr7ySm2++Oaeeemqqqso73/nOvOMd79hk7uc+97kkyUUXXZQrr7wyy5cvT5I0NTVlwoQJ+fGPf5w+ffpkwoQJr3nctswFAAAA6hQKzj777Pzrv/5rGhoa8rWvfS177bVXevTokb59++aUU07Jc889l1GjRuWmm256zdxjjz02F154YdasWZPzzjsvffr0Sa9evTJs2LBcffXV2W233fK9730vAwcOrOlcAAAAoI73KPjMZz6TRx99NB//+MczfPjwNDc3Z4899sihhx6ar33ta3n44YczfPjwVud++ctfzm233ZbjjjsuAwYMyPr167PXXnvltNNOy29+85uceOKJm33ctswFAACAzq6xnjs/8MAD8+1vf/t1zf3gBz+YD37wg+0+FwAAADqzDnt7RAAAAGDHIxQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUjR29AAAA6ExGTLqrTfMXTD+hRisBaJ0zCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAACKxo5eAAB0hBGT7uroJQAA7JCcUQAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAEVjRy8AAADYdiMm3dXmfSyYfkINVgLsqpxRAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABA0djRCwCA7VWL9xAHAKB1zigAAAAACqEAAAAAKNolFDz88MMZO3ZsBg0alG7dumXPPffM8ccfn7vvvnuL8x544IGcfPLJGTZsWLp3756hQ4dm7Nixuf/++7f6mG2ZCwAAAJ1V3UPBRRddlEMPPTS33nprXnjhhXTp0iVLlizJ3XffneOPPz6XXnppq/NmzpyZo446KrfcckuamprStWvXLF68OLfeemuOPvrozJgxY7OP2Za5AAAA0JnVNRT8+7//e6ZNm5bddtst//Iv/5IlS5Zk9erVeeqppzJmzJgkyRe+8IXMnz9/k3lPPvlkJkyYkObm5kyePDlLlizJihUrsnTp0kyZMiXNzc0577zzMnv27Nc8ZlvmAgAAQGdXt1Cwbt26TJ06NUly2WWX5XOf+1wGDBiQJHnjG9+YH/zgBxkyZEjWrVuXH/zgB5vMveyyy7J27dqMGzcul1xySZnXv3//TJ06Neecc07WrVvX6tkIbZkLAAAAnV3dQsHatWszYcKEHHHEERk3btxrtvfs2TPHHHNMkuT3v//9JtvuuOOOJMn555/f6r4nTpyYJLnzzjtTVVXN5gIAAEBnV7dQ0LNnz3zuc5/LAw88kMbGxlbHDB8+PEny8ssvl8/Nnz8/TU1N6devX0aPHt3qvP333z9DhgzJkiVLNrmEoC1zAQAAgA5+e8Q1a9YkSd7whjeUz82dOzdJsvfee29xbsv2efPm1WQuAAAAkLT+q/528sQTTyTZcM+CFgsXLkyy4YyELenVq9cm49s6d1ssWrRoi9ubmpq2a38AAACwo+mwUPD8889n1qxZSZIPfOAD5fPLly/frv1sPL4tc7fF1s5UAAAAgJ1dh1168NnPfjbr1q3LO9/5zrzzne8sn1+5cuV27WfFihU1mQsAAAB00BkF119/fX7wgx9k9913z8yZMztiCa/Ls88+u8XtTU1NOeyww9ppNQAAAFB77R4KfvrTn+acc85JQ0NDvvOd72TUqFGbbN/a/QX+Usv9Bto6d1u0vEsDAAAA7KraNRTcd999+dCHPpS1a9fmG9/4RsaOHfuaMb17996ufW48vi1zAQAAgHa8R8HDDz+cD3zgA1mxYkW+8IUvZOLEia2O22effZIkq1at2uL+Wu4v0DK+rXMBAACAdgoFjz/+eN73vvfllVdeyfjx4/PVr351s2NHjhyZZOv3A2jZ3jK+rXMBAACAdggF8+bNy3ve854sXbo0H/7wh3PVVVdtcfzIkSMzdOjQLFu2LHPmzGl1zPz587N48eIMGDAgo0ePrslcAAAAoM6hYOHChXn3u9+dxYsX5/3vf3+++93vpkuXrT/kSSedlCS54oorWt1++eWXJ0lOPPHENDQ01GwuAAAAdHZ1CwXLly/P3/7t32bhwoU58sgjc8stt6Rbt27bNHfixInp1q1bZsyYkalTp2bZsmVJkhdffDHTpk3LVVddlcbGxlxwwQU1nQsAAACdXd1CwZIlSzJ37twkye9+97sMHTo0b3jDGzb7cdNNN5W5o0aNypVXXpkuXbrk4osvzoABA7L77runf//+ueiii5Ik3/zmN1u9dKAtcwEAAKCza5e3R/zzn/+81TFr1qzZ5O/jxo3LX//1X+cb3/hGfv3rX2fp0qUZPHhwDj/88EycODFHH330ZvfVlrkAAADQmdUtFIwYMSJVVbVpH0cddVSOOuqodp8LAAAAnVW7vD0iAAAAsHMQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoGjs6AUA0LmMmHRXRy8BAIAtcEYBAAAAUAgFAAAAQCEUAAAAAIV7FAAAQCdTi/vFLJh+Qg1WAuyInFEAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQNHb0AgDYeYyYdFdHLwEAgDpzRgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQNHY0QsAAAB2PiMm3dXmfSyYfkINVgLUmjMKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAACKxo5eAADtY8Skuzp6CQAA7AScUQAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABSNHb0AAACgcxox6a4272PB9BNqsBJgY84oAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAot1Cwec///k0NDSUjxkzZmx1zgMPPJCTTz45w4YNS/fu3TN06NCMHTs2999/f13nAgAAQGdV91BQVVUmTJiQr3/96+nXr1969uy5TfNmzpyZo446KrfcckuamprStWvXLF68OLfeemuOPvroLYaGtswFAACAzqyuoWD9+vU566yzcvXVV2fQoEGZNWtWBg0atNV5Tz75ZCZMmJDm5uZMnjw5S5YsyYoVK7J06dJMmTIlzc3NOe+88zJ79uyazgUAAIDOrm6hYO3atTnttNNyww03ZPjw4fnlL3+Zt7zlLds097LLLsvatWszbty4XHLJJRkwYECSpH///pk6dWrOOeecrFu3LpdeemlN5wIAAEBnV7dQMGfOnNx5553Zb7/9ct9992XUqFHbPPeOO+5Ikpx//vmtbp84cWKS5M4770xVVTWbCwAAAJ1d3ULBwQcfnNtvvz333XdfRowYsc3z5s+fn6ampvTr1y+jR49udcz++++fIUOGZMmSJZtcQtCWuQAAAECd71Fw3HHHZdiwYds1Z+7cuUmSvffee4vjWrbPmzevJnMBAACApLGjF/CXFi5cmCRbfXeEXr16bTK+rXO3xaJFi7a4vampabv2BwAAADuaHS4ULF++/HWPb8vcbbG1MxUAAABgZ1fXSw9ej5UrV27X+BUrVtRkLgAAALADnlGwI3v22We3uL2pqSmHHXZYO60GAAAAam+HCwVbu7/AX2q530Bb526L4cOHb9d4AAAA2NnscJce9O7d+3WPb8tcAAAAYAcMBfvss0+SZNWqVVsc13J/gZbxbZ0LAAAA7IChYOTIkUm2fj+Alu0t49s6FwAAANhBQ8HQoUOzbNmyzJkzp9Ux8+fPz+LFizNgwICMHj26JnMBAACAHTAUJMlJJ52UJLniiita3X755ZcnSU488cQ0NDTUbC4AAAB0djtkKJg4cWK6deuWGTNmZOrUqVm2bFmS5MUXX8y0adNy1VVXpbGxMRdccEFN5wIAAEBnV9dQ8O53vzuNjY2bfDzzzDNJkgkTJmzy+Xe/+91l3qhRo3LllVemS5cuufjiizNgwIDsvvvu6d+/fy666KIkyTe/+c1WLx1oy1wAAADo7OoaCtavX/+ajxbNzc2b3ZYk48aNy6xZszJmzJgMGTIka9euzeDBgzNmzJjMmjUr48eP3+zjtmUuAAAAdGaN9dz5L37xizbNP+qoo3LUUUe1+1wAAADorHbIexQAAAAAHUMoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKBo7egEAbJsRk+7q6CUAANAJCAUAAMBOqxYhfcH0E2qwEth1uPQAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACgaO3oBAAAAHWnEpLvavI8F00+owUpgx+CMAgAAAKBwRgFAO6jFbyoAAKA9OKMAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAACKxo5eAMCObsSkuzp6CQAA0G6cUQAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQNHb0AgDqacSkuzp6CQAAsFNxRgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAABFY0cvAAAAYGc3YtJdbd7Hgukn1GAl0HbOKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgaOzoBQAAAJCMmHRXm/exYPoJNVgJnZ0zCgAAAIDCGQXADqsWVR0AANg+zigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAArvegAAALCLqMW7Ri2YfkINVsLOzBkFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUDR29AKAHc+ISXd19BIAAIAO4owCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACjczBB2MW5ECAAAtIUzCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgcDNDAAAAilrcHHvB9BNqsBI6ilAANeIbKgAAbNDWn439XNyxdulLDx544IGcfPLJGTZsWLp3756hQ4dm7Nixuf/++zt6aQAAALBD2mXPKJg5c2bGjx+f5ubmJEnPnj2zePHi3Hrrrbn99ttz1VVXZfz48R28SthULc5KAAAAaItd8oyCJ598MhMmTEhzc3MmT56cJUuWZMWKFVm6dGmmTJmS5ubmnHfeeZk9e3ZHLxUAAAB2KLtkKLjsssuydu3ajBs3LpdcckkGDBiQJOnfv3+mTp2ac845J+vWrcull17awSsFAACAHcsueenBHXfckSQ5//zzW90+ceLEXHPNNbnzzjtTVVUaGhrac3kAAABsgRuFd6xdLhTMnz8/TU1N6devX0aPHt3qmP333z9DhgzJ4sWLM3v27Bx44IHtvEoAAADqSWx4/Xa5Sw/mzp2bJNl77723OK5l+7x58+q+JgAAANhZ7HJnFCxcuDDJhnc52JJevXptMp7OzbsNAAAAbLDLhYLly5fXbfyiRYu2uP3ZZ58tf25qatqudXSEw7/6323ex//vi+/eIdYBAB1hzapX8pc/HaxZvjTr1q3tkPUAUFtbew24I9j4tee6detqss9dLhSsXLlyu8avWLFim8du7XKGjR122GHbtY6d1d7XdPQKAKDjPJfkNT8d/PunOmAlANTDzvZ654UXXsiIESPavJ9d7h4FAAAAwOu3y51RsLV7E/yllnsVbIuNLy1ozapVqzJnzpwMHjw4e+65Zxobd4z/vE1NTeUMh4ceeihDhw7t4BVRK47trslx3XU5trsmx3XX5djuuhzbXVNnPa7r1q3LCy+8kCQ5+OCDa7LPHeOVbA317t27buOHDx++1TEjR47crsdvb0OHDt2m58HOx7HdNTmuuy7HdtfkuO66HNtdl2O7a+psx7UWlxtsbJe79GCfffZJsuG3+1vScm+ClvEAAADALhgKWn6jv7XLBFq27+hnAAAAAEB72iVDwdChQ7Ns2bLMmTOn1THz58/P4sWLM2DAgIwePbqdVwgAAAA7rl0uFCTJSSedlCS54oorWt1++eWXJ0lOPPHENDQ0tNeyAAAAYIe3S4aCiRMnplu3bpkxY0amTp2aZcuWJUlefPHFTJs2LVdddVUaGxtzwQUXdPBKAQAAYMeyS4aCUaNG5corr0yXLl1y8cUXZ8CAAdl9993Tv3//XHTRRUmSb37zmy47AAAAgL+wS4aCJBk3blxmzZqVMWPGZMiQIVm7dm0GDx6cMWPGZNasWRk/fnxHLxEAAAB2OA1VVVUdvQgAAABgx7DLnlEAAAAAbD+hAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqFgJ/TAAw/k5JNPzrBhw9K9e/cMHTo0Y8eOzf3339/mff/xj3/MAQcckIaGhvzN3/xNnnvuuRqsmG1Vq2M7derUNDQ0bPMH9VXrr9m5c+dmwoQJGT16dHbffffstttuGTZsWE488cTcfvvttV08W1TrY7tmzZr827/9W4444oj069cvPXr0yIgRI/KP//iPmTdvXo1Xz9Z8/vOf3+R75YwZM9q0v+9+97s55phjMnDgwPTo0SP77bdfzjnnnCxcuLBGK2Zb1Pq4Ll++PMcee+wm+5wzZ06NVsv2qOWxraoqP/zhD3PUUUelT58+6d69e/bdd9987GMfy+zZs2u4arZFLY/tM888k09/+tN54xvfmN122y29e/fO2972tkydOjWvvPJKDVe9k6vYqVx33XVVly5dqiRVkqpnz57lz126dKmuueaa173vP/zhD9WgQYOqJNUpp5xSrVq1qoYrZ2tqeWynTJlSJam6d+9e9e3bt9WP3r17l/1TP7X+mv3xj39cde/eveyjoaGh2m233crfk1RnnnlmnZ4NG6v1sX3ppZeqww8/fJN9bLzPHj16VPfcc0+dng0ba25urs4999wqSdWvX79yHNryb+xHP/rRTb5ue/ToUf7ep0+f6je/+U0NnwGtqcdxfemll6r/9b/+V5Wk2muvvcoxfeKJJ2q4cram1sd21apV1XHHHbfJ12y3bt3K33v16lXde++9NX4WtKbWx/buu+/e5GfgjX+mSlKNHDmy+n//7//V+FnsnLxC2InMmTOnfJOaPHlytWTJkqqqqmrp0qXlhWFjY2P1hz/8Ybv3vXjx4mr48OFVkuoTn/hE1dzcXOvlswW1PrYtc6ZMmbLZMd/4xjeqJNUhhxxSi6dAK2p9XFeuXFntueeeVZJq9OjR1axZs6pVq1ZVzc3N1dNPP1194hOfKP/Q3XrrrfV8ap1ePb4ff+xjH6uSVHvuuWd1yy23VKtXr66qqqqeeeaZ6qyzzio/JPkBpr7WrVtXnXnmmVWSatCgQdUjjzxS7bvvvm36wfTb3/52CQI33nhjtXLlyqqqqmrevHnVhz70oSpJNWLEiGrNmjW1fCpspB7H9YUXXqje+ta3Vkmqww8/vPrd734nFHSAehzbM844o8Sfm2++ufrzn/9crV27tnrwwQer0aNHV0mqv/qrv6rWr19f42fDxmp9bGfPnl1Cw7nnnlstXLiwqqqqWr58eXXDDTeUbePGjav1U9kpCQU7kbPPPnuL//Oec845VZLq4x//+Hbv+93vfneVpPrQhz7km14HqPWx3ZZQ0PLDzaWXXvp6lsw2qPVx/e///u/yQ+ijjz76mu3Nzc3V2972tipJdcYZZ7Rp7WxZrY/twoULq65du1ZJqlmzZrU6puW3W5MnT369y2Yr1qxZU5166qlVkmr48OHVnDlzqqqq2vyiY9SoUVWS6qabbmr1Md/85jdXSarvfOc7bVo/ravHcX3++efLC8ZjjjmmevXVV6unn35aKGhn9Ti2c+bMKWcQPPLII6/Z/tvf/rYc5/vuu6/Nz4HW1ePY/ulPf6pOOOGE6lOf+lSr21t+fh40aFCb1r6rEAp2IkOHDq2SbPY3VH/84x+rJNXAgQO364yAlt907LvvvtWrr75aq+WyHWp9bLcWCh5//PFyavPzzz/flqWzBbU+rjfddFOVpNpjjz02O+af/umfqiTVcccd97rXzdbV+tjOnDmzSlLtv//+mx1z8803V0mqgw8++HWvmy179NFHq549e1b77bdf9fTTT5fPt/VFR5Jq8ODBmw3x1113XZWk+ru/+7vXu3S2oB7HteVr9vjjjy9niAgF7a8ex3bhwoXV+eefX5111lmbHTNkyJAqSXXllVe+nmWzDepxbFts7uytH//4x+Vr2C9Oq6ox7BTmz5+fpqam9OvXL6NHj251zP77758hQ4Zk8eLFmT17dg488MCt7nfNmjW56KKLkiRf//rX07t375qum62rx7EdPXp0xo4du9n93XjjjUmSY489NkOHDm3bE6BV9Tiu23Oshg0btl3rZdvV49g+8cQTSZIDDjhgs2Natj3xxBNZvXp1unfv/jqfAZtz8MEH5/bbb89BBx1Us6+hlhtbHnHEEenSpfV7SB911FFJkl/+8pc1eUw2VY/j+o//+I/p2bNnTj311HTr1q0m+2T71ePY7r333rn88su3OGb48OFZvHhxXn755Zo8Jq9Vj2PbYnNfs2vWrEmS9OvXb7PfrzsT/wV2EnPnzk2y4ZvXlrRs39a7Y990001ZtGhR3vGOd+TUU09t2yJ5XepxbE899dTcfPPNrR7T5ubmfP/730+S/MM//MP2LpdtVI/j+ta3vjV9+/bNq6++mscee+w126uqyoMPPpgkOeaYY7ZzxWyrehzblrss9+zZc7NjevXqlSRZt25dli5duk1rZfsdd9xxNf2hdFv+f2nZtmTJEi886qTWxzVJTj/9dJFgB1CPY7s1LS8o3/CGN7Tr43Y27X1sf/KTnyRJ/v7v/77dHnNHJhTsJFreOmlLP0Qm//OD5La+1dJ3vvOdJMnEiRNz33335bjjjsvgwYPzhje8IUcccUS+/e1vp6qqNqycranXsd2cWbNmZdGiRenZs2f+7u/+rk37YvPqcVz79u2bL3/5y0mSD3/4w7n33nuzZs2aVFWVZ555JuPGjctvf/vbHH744SJQHdXj2LacHfDnP/95s2OWL19e/rxy5cqt7pMdw7b8/9Ly/8rG44Ed05o1a/LUU08lSd74xjd28Gpoq6qq8vzzz+fCCy/Mt771rey///6ZNm1aRy9rh+DSg53Exj8g1mr8kiVL8stf/jK77757unTpkmOPPTbr1q0r2x988ME8+OCDmTVrVjlVndqrx7HdkpZjedJJJ2WPPfZo077YvHod13/6p3/K4MGDM3369BxzzDFpaGhIt27dsmbNmvTv3z8XXHBBpk6dmq5du76eZbMN6nFsW36j/Pjjj292zMbb+vTps11roOO09/d4oL5+9KMfZfny5enTp0+5bIidz7nnnpubbropq1evzqpVqzJo0KBMmjQpn//8550p8v9xRsFOYnt/e7RixYqtjvn1r3+d5ubmHHTQQZkyZUquvfbavPDCC1m9enV+//vfZ+zYsUmS733ve+VUdWqvHsd2S3NvueWWJC47qLd6HtfVq1eXF4pVVZVTIHv27JkePXq06f8Rtq4ex/bII49MkixYsCB33333a7ZXVZVrrrkmSdK1a9f069dvu9ZAx2nP7/FAfa1cuTKTJk1Kkpx//vnZfffdO3hFvF4rVqzIyy+/nFWrViVJVq1alXnz5mXBggUdu7AdiFDQif3+979PsiEY/Mu//EvOOuusDBw4MLvttlve/OY354c//GEppVdffXVHLpUauf3227N8+fIMHDgw73vf+zp6OWynqqpyxhln5KMf/Wiee+653HzzzXnuuefyyiuv5De/+U3e8Y53ZNq0aTn00EMzf/78jl4u2+Ed73hHDjnkkCTJRz7ykXz729/O0qVLs3bt2jz22GM55ZRTMmTIkCTJQQcdlMZGJwQCtLezzz47Tz31VA466KB88Ytf7Ojl0AY33HBDqqrK6tWrM3v27Jxzzjm5/fbbc/jhh+cXv/hFRy9vhyAU7CS2di3sX9r4esfNWbJkSZIN11e9//3vf832hoaGnHvuuUmShx9+OGvXrt2uNbBt6nFsN+e73/1ukg03afFCo77qcVxvuumm3HjjjenXr1/uu+++jB07NsOGDcsee+yRt7/97fk//+f/5NRTT82CBQvK1y61V6+v2e9973sZNGhQXn755XziE5/YJNy+8MILOeigg5JsiArsPNrzezxQP1OmTMmNN96Y/v3757bbbkuPHj06eknUwG677ZYDDjgg06dPz6RJk7J69eqMGzfOPdoiFOw0tvdtC7dlfMtdtrf0tl0t29asWZM//elP27UGtk09jm1rFi9enP/6r/9K4rKD9lCP4/of//EfSTb8xnnw4MGtjvnMZz6TJPnZz36WF198cbvWwLap19fsgQcemEceeSSf/vSnc8ABB2TQoEE55JBDcvHFF+cnP/lJbr311iTJySefvN1rpuO01/d4oH6+8Y1v5JJLLknPnj1zxx13ZOTIkR29JOrgn/7pn5JseLealjOvOzO/UtxJ7LPPPklSrqPZnJZrG1vGb0nLbzm2dH3Vxj+wbO2xeX3qcWxbc9NNN2X9+vX5q7/6qxx++OGvax9su3oc15a3Wdt///03O2bUqFFJNlymMG/evBx66KHbtF62XT2/ZocOHZrLLrvsNZ//1a9+lUceeST77rtvjj766O1YLR1tW/5/2fi+BK/3ezxQH9ddd13++Z//OY2NjfnhD3/orK5d2ODBg9O7d+8sX748c+fOLZcEdlbOKNhJtJTLZ599dovjWrZvS+ls+Y3klm60tPFbdfXt23er+2T71ePYtqbl3Q5OP/301zWf7VOP47p+/fokSZcum//W3dDQ8Jrx1FZ7fc1u7Etf+lKS5POf//wWjz87nm35/6Vl28CBA/1bCzuQ73//+znnnHPS0NCQ73znOznhhBM6eknUWcu/sc3NzR28ko7np42dxMiRIzN06NAsW7Ysc+bMaXXM/Pnzs3jx4gwYMCCjR4/e6j7f+ta3Jklmz5692TFPPPFEkmTQoEEZOHDg61g5W1OPY/uXHn/88TzyyCNJXHbQXupxXIcOHZrkf84saM0f//jH8udhw4Zt56rZFu3xNbuxb33rW5k1a1be9KY35ayzzmrTvmh/73znO5NseMvhzV3zev/9928yFuh4t99+e84888w0NzfnyiuvzGmnndbRS6KNVqxYkR/96Eeb3b5kyZJyabafoYSCncpJJ52UJLniiita3X755ZcnSU488cRNfqu4Occee2z22GOPzJs3b7N395w5c2aSuEN+ndX62P6llpsYHnbYYVs8bZ3aqvVxfde73pVkw70Kli1b1uqYq666Kkmy3377OYW5jur9NdvioYceynnnnZeGhobMmDEj3bt3f937omOMGjUqo0aNSlNTU7nPyMbWrVuXK6+8MkkyZsyY9l4e0Iqf/vSn+fCHP5x169Zl2rRpbhC8C3jllVdyxBFH5JRTTsmvf/3rVse0vO7ZY489XLqZJBU7jTlz5lTdunWrklRTpkypli5dWlVVVS1btqy65JJLqoaGhqqxsbH6wx/+UOa89NJL1Vvf+taqR48e1Uc/+tHX7POiiy6qklSDBg2qbrnllmr16tVVVVXVwoULq7POOqtKUnXr1q169NFH2+dJdlL1OLYt1q9fX+21115Vkurf/u3f6v5c+B+1Pq6LFy+u+vbtWyWpDjnkkOq///u/qz//+c/V+vXrq3nz5lVnn312laRKUs2cObNdn2tnU8+v2Ra/+93vqv79+1dJqksuuaRuz4Wt23fffask1TXXXNPq9q0d2+uvv75KUvXt27f6/ve/X61ataqqqqqaP39+NXbs2CpJte+++5Z/g2kfbT2uf+npp58u34OfeOKJWi+X7dCWY/u73/2u6tWrV5Wk+sxnPtMey2U7tOXYfvjDH66SVP369auuv/766uWXX66qqqqampqqr33ta1WPHj2qJNXkyZPr/jx2BkLBTubaa6+tunTpUv4havlGlqRqaGh4zRfNfffdV7Y3Nja+Zn+rV6+uTjjhhDKma9eum+yzS5cu1fXXX99eT69Tq/WxbfGzn/2sjPnTn/5U76fBX6j1cf3P//zPao899thkHy0vWFs+JkyY0F5Pr1Or19dsVVXVz3/+8xKFzj777Kq5ubmeT4WNHHvssVXXrl03+dj438SNP3/sscdWVbVtx/b000/fZD8tP5Amqfr06VM99NBD7fk0O516HNe/3N/G+/zLz1988cXt/ZQ7jVof23//938v2/r27bvVj2eeeaYjnnanUOtju2rVqurv//7vN/mZqXv37pv8fcyYMdXatWs74unucISCndC9995bjRkzphoyZEjVrVu3avDgwdWYMWOqX/ziF68Z+9JLL1WHHHJI1b17983W8PXr11ff+ta3qiOPPLLq27dvtdtuu1XDhw+vTjvttOq3v/1tvZ8OG6n1sa2qqjrjjDOqJNXxxx9fz6WzBbU+rgsWLKg+/elPVwceeGC1++67V926dauGDRtWfehDH6ruueeeej8dNlKPr9nrrruuxJ8LL7xQJGhnRx999CY/NG7p4+ijj66qatuP7fXXX18deeSRVb9+/aru3btX++67b/XJT36yevrpp9vnyXVi9Tiu27q/ljOPqI9aH9uNQ8G2fPj6rZ96fT++++67q+OPP77q379/1djYWA0YMKB6z3veU33ve9/zb+5GGqpqM3fWAQAAADodNzMEAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAAiv8/o4nt/uFU3DgAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 600x600 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 505,
       "width": 517
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sim_cluster(5, 80, 50, 30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.010911041476145415, 0.010650749123757628)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABAoAAAPyCAYAAAD1wj3WAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAB7CAAAewgFu0HU+AABLuElEQVR4nO3de7yVdYHv8e+GzVUuAcpNVHJwSNSyTno8Znhp8pSm1KBO6WRmE14Yx2iyKFMQujC9HC+TF7zUWJadM3lLs5yaiUzNV2YnU0MMKER0cxJQELmzn/MHZ/8GdHNzr8UG9vv9eu3Xa7ue3+9Zv+XDvqzPftazGqqqqgIAAACQpFN7LwAAAADYeQgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFA0tvcCdierVq3Kk08+mSTZa6+90tjofy8AAAD1s27durz44otJkkMOOSTdu3dv8z49k62hJ598Mocffnh7LwMAAIAO6NFHH81hhx3W5v146QEAAABQOKOghvbaa6/y+aOPPpohQ4a042oAAADY3TU1NZUz2zd+TtoWQkENbXxNgiFDhmTYsGHtuBoAAAA6klpdJ89LDwAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAorG9FwAA0GL4xPvaNH/etBNrtBIA6Lh26BkFy5cvzwc/+ME0NDRk8ODBeeihh7Y4/uGHH84pp5ySoUOHplu3bhkyZEjGjh271XltnQsAAAAd1Q4LBa+88kqOP/74/OAHP8hBBx2UX//61znqqKM2O/6mm27K6NGjc8cdd6SpqSmdO3fOwoULc+edd+boo4/O9OnT6zIXAAAAOrIdEgqam5tz6qmn5pFHHsk73vGOPPjgg9lnn302O/6ZZ57J+PHj09zcnEsuuSSLFi3KihUrsnjx4kyaNCnNzc254IILMnPmzJrOBQAAgI5uh4SCL3/5y/n3f//37L///vnJT36Sfv36bXH8lVdembVr12bcuHGZMmVKBgwYkCTp379/Jk+enPPOOy/r1q3L5ZdfXtO5AAAA0NHVPRQ888wz+dKXvpTOnTvnjjvuKE/ct+See+5Jklx44YWtbp8wYUKS5N57701VVTWbCwAAAB1d3UPBpZdemjVr1mTcuHE59NBDtzp+7ty5aWpqSr9+/TJq1KhWxxxwwAEZPHhwFi1atMlLCNoyFwAAAKhzKPjjH/+Y22+/PX379s2UKVO2ac7s2bOTZIvXMNh4+5w5c2oyFwAAAEga67nzW2+9Nc3NzfnoRz+a5cuX59Of/nQeeOCBLF68OPvuu29OPfXUfOYzn0nv3r3LnPnz5ydJevToscV99+zZc5PxbZ27LRYsWLDF7U1NTdu1PwAAANjZ1DUU3HnnnUmSI444IkceeeQmT6SffvrpTJkyJXfddVceeOCBcoHD5cuXb9d9bDy+LXO3xdbOVAAAAIBdXd1CwSuvvJKnnnoq3bp1y/e///2cddZZGTduXPbee+80NTXl5ptvzpe//OU8+eST+exnP5ubbropSbJy5crtup8VK1aUz9syFwA6suET72vzPuZNO7EGKwEA2lvdQsGTTz6Z5ubmrF69OiNHjsxXvvKVsm3fffct1yyYOnVqvvOd7+SKK67Y5CUIO6Pnnntui9ubmppy+OGH76DVAAAAQO3VLRQsWrSofH7++ee3Oub888/P1KlTs2rVqvyf//N/cvTRR2/1+gKv1XK9gWTr1ybY0txtMWzYsO0aDwAAALuaur3rwbJly5IkvXr1yn777dfqmMGDB2fAgAFJ/utCgb169dqu+9l4fFvmAgAAAHUMBS1/3d9jjz22OK7lyfqqVauSbHhZwsb/vTkt1xdoGd/WuQAAAEAdQ8GgQYOSbP0Cg6+++mqSpG/fvkmSESNGJNn69QBatreMb+tcAAAAoI6h4K1vfWs6deqUZcuW5fnnn291zKJFi8q1DEaNGpVkw5P3IUOGZMmSJZk1a1ar8+bOnZuFCxdmwIABZV5b5wIAAAB1DAV9+vTJsccemyTlrQ9fq+X24cOH58ADDyy3n3zyyUmSq6++utV5V111VZLkpJNOSkNDwybb2jIXAAAAOrq6hYIk+cIXvpAk+fKXv5yvfe1refnll5MkS5cuzZVXXplJkyYlSSZOnLjJk/YJEyakS5cumT59eiZPnpwlS5YkSV566aVMnTo11157bRobG3PRRRe97j7bMhcAAAA6uoaqqqp63sGXv/zlfPGLX9xwZw0N6dmzZ1asWJGWuz377LNz8803v+6v+zfeeGPOO++8NDc3J0mZ17Kf6667Lueee26r99mWuW2xYMGC7LPPPkk2XAfB2ykCsKsYPvG+Nu9j3rQT230dtVgDAOxK6vE8tLHNe9iKiy++OEceeWSuvvrqPPLII3nppZcyYMCA/Lf/9t/yyU9+MmPHjm113rhx4/KWt7wlV1xxRX71q19l8eLFGTRoUI444ohMmDAhRx999Gbvsy1zAYBd184SPABgV1b3UJAkxx57bLlewfYYPXp0Ro8e/Ybusy1zAQAAoKOq6zUKAAAAgF2LUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFA0tvcCAIDdw/CJ97X3EgCAGnBGAQAAAFAIBQAAAEDhpQcAsItzyj8AUEvOKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAorG9FwAAsDMZPvG+Nu9j3rQTa7ASAGgfzigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAICirqHg5z//eRoaGrb6MW3atFbnP/zwwznllFMydOjQdOvWLUOGDMnYsWPz0EMPbfW+2zIXAAAAOqodckZBly5d0rdv381+dO/e/XVzbrrppowePTp33HFHmpqa0rlz5yxcuDB33nlnjj766EyfPn2z99eWuQAAANCR7ZBQcPrpp+fll1/e7MenPvWpTcY/88wzGT9+fJqbm3PJJZdk0aJFWbFiRRYvXpxJkyalubk5F1xwQWbOnPm6+2rLXAAAAOjodsprFFx55ZVZu3Ztxo0blylTpmTAgAFJkv79+2fy5Mk577zzsm7dulx++eU1nQsAAAAdXWN7L6A199xzT5LkwgsvbHX7hAkTcv311+fee+9NVVVpaGioyVwAgFoYPvG+Nu9j3rQTa7ASANh+O90ZBXPnzk1TU1P69euXUaNGtTrmgAMOyODBg7No0aJNXkLQlrkAAADAThgKZs+enSTZZ599tjiuZfucOXNqMhcAAADYgS89WLJkSSZNmpSf/vSnWbp0aUaMGJEzzzwzZ599djp37lzGzZ8/P0nSo0ePLe6vZ8+em4xv69xtsWDBgi1ub2pq2q79AQAAwM5mh4SC9evX53/+z/+Zxx57rNy2cOHCPPTQQ/nRj36UO+64I506bTi5Yfny5du1743Ht2XuttjamQoAAACwq9shLz249957079//zzxxBNZtWpVnnvuuVx66aVpaGjI3Xffneuvv76MXbly5Xbte8WKFTWZCwC0r/4rlmbeP31gk4/+K5a297IAoMOp6xkFw4cPzznnnJPVq1fnuuuuKy8JGDZsWC677LLMnz8/t9xyS66++uqMHz++nkupieeee26L25uamnL44YfvoNUAAABA7dU9FEyfPn2z288555zccsstmT17dp5//vnsvffeW72+wGu1XG8g2fq1CbY0d1sMGzZsu8YDAADArqZd3/VgxIgR5fNnn302SdKrV6/t2sfG49syFwAAAGjnULDnnnuWz1etWpUk2XfffTf5781pub5Ay/i2zgUAAADaORQsXry4fN63b98k/3WWwdauB9CyfeOzEtoyFwAAAKhzKFi9enVuueWWzW6fO3fuhkV06pS//Mu/TLLhyfuQIUOyZMmSzJo1a7PzFi5cmAEDBmTUqFHl9rbMBQAAAOoYCtavX59jjz02H//4x/P973+/1TE33HBDkuRd73pXevfuXW4/+eSTkyRXX311q/OuuuqqJMlJJ52UhoaGTba1ZS4AAAB0dHULBZ07d8673/3uJMkZZ5yRadOmlZcavPTSS/nKV75Szjb4whe+sMncCRMmpEuXLpk+fXomT56cJUuWlHlTp07Ntddem8bGxlx00UWvu9+2zAUAAICOrq4vPfjKV76S8ePHZ+3atfn85z+fPffcMz179kz//v1z8cUXp7m5OVOmTMn73ve+TeaNHDky11xzTTp16pTLLrssAwYMyB577JH+/fvn0ksvTZJ8/etfb/WlA22ZCwAAAB1dXUNB586dc80112TGjBn54Ac/mH333Tfr16/PoEGD8sEPfjAzZszIJZdc0urccePGZcaMGRkzZkwGDx6ctWvXZtCgQRkzZkxmzJiRc889d7P325a5AAAA0JE17og7OeaYY3LMMcds97zRo0dn9OjRb+g+2zIXAAAAOqp2fXtEAAAAYOciFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUDS29wIAoCMbPvG+9l4CAMAmnFEAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQNHY3gsAgF3V8In3tfcSAABqzhkFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFC4mCEAwE6oFhfLnDftxBqsBICOxhkFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAARWN7LwAAgPoYPvG+Nu9j3rQTa7ASAHYlzigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoGhs7wUAQHsYPvG+9l4CAMBOyRkFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUOzQULF68OH369ElDQ0MaGhoyefLkLY5/6qmn8rGPfSz77rtvunXrloEDB+b9739/7r333q3eV1vmAgAAQEe1Q0PBtGnT8sorr6SxsXGrY3/0ox/lne98Z7797W/nueeeS+fOnfPiiy/m/vvvz8knn5yLL764LnMBAACgI9thoaCpqSnXXnttDj300Jx++ulbHPvSSy/ljDPOyOrVq/PJT34yCxYsyIoVK/LKK6/kmmuuSdeuXfOVr3wlP/nJT2o6FwAAADq6HRYKvvSlL2XlypWZNGlSGhoatjj2hhtuyMsvv5zjjz8+N954Y/bee+8kSa9evTJ+/PhMnTo1yYYzFGo5FwAAADq6HRIKnn322dx888059NBDM2bMmK2O/8EPfpAkufDCC1vdfv7556dbt275xS9+kZdffrlmcwEAAKCj2yGh4LLLLsuaNWu26WyCVatW5bHHHkuSHHXUUa2O6dWrV97+9rdn/fr1efjhh2syFwAAANgBoeCZZ57Jt7/97W0+m2DevHlZt25d+vTpkz59+mx23D777JMkmTNnTk3mAgAAAMnW336gjSZNmpT169dv09kESTJ//vwkSY8ePbY4rmfPnpuMb+vcbbFgwYItbm9qatqu/QEAAMDOpq6h4Iknnsi//du/bfPZBEmyfPny7bqPjce3Ze62aDkTAQAAAHZXdX3pwSWXXJKqqrb5bIIkWbly5Xbdx4oVK2oyFwAAAKjjGQWPPvpo7rnnnu06m2Bn99xzz21xe1NTUw4//PAdtBoAAACovbqFgosvvjhJtutsgmTr1xd4rZbrDbR17rYYNmzYdo0HAACAXU1dXnrwwAMP5D/+4z/e0NkEvXr1esPj2zIXAAAAqFMo+OIXv5hk+88mSJJ99903SbJq1aotjmu5vkDL+LbOBQAAAOr00oOHHnooSXLWWWe1ur3lifq0adNy1VVXJUlefvnlJMnw4cPTuXPnLF26NMuWLUufPn1a3UfL9QJGjBhRbmvLXAAAAKDO73qwdOnSVj/Wrl2bJFm9enW5rUX37t1z2GGHJUl++ctftrrfV199Nb/97W/TqVOnHHnkkTWZCwAAANQpFFRVtcWPj33sY0k2vDSh5baNtVzX4Oqrr251/9OnT8/q1aszevTo9OvXr2ZzAQAAoKOr6xkFb9S4cePSt2/f3H///Tn33HPzwgsvJNlwNsD111+fL3zhC0mSz3/+8zWdCwAAAB3dThkK+vfvn+9+97vp2rVrbrjhhuy9997ZY4890rt375x//vlZs2ZNJk6cmOOPP76mcwEAAKCj2ylDQZKceOKJ+fWvf53TTz89e++9d9atW5f+/fvnve99b+6666589atfrctcAAAA6Mjq8q4HW3PLLbfklltu2eq4t771rfnud7/7hu6jLXMBAACgo9ppzygAAAAAdjyhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAIrG9l4AAAA7r+ET72vzPuZNO7EGKwFgR3FGAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFI3tvQAA2F7DJ97X3ksAANhtOaMAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoGhs7wUAAMDWDJ94X5v3MW/aiTVYCcDuzxkFAAAAQOGMAgAA6qoWZwMAsOM4owAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACga23sBAHQswyfe195LAABgC5xRAAAAABRCAQAAAFDUNRS88sor+ed//ue87W1vS48ePdK9e/cccMABOffcczN79uwtzv3hD3+Y97///Rk0aFC6deuWffbZJ2eeeWaeeuqprd5vW+YCAABAR1a3UPDHP/4x73znO/OZz3wmTzzxRNavX5+1a9dmzpw5ueGGG3LooYfmoYceanXuJZdckpNOOin3339//vznP6dz585ZsGBBbr311rzzne/Mvffeu9n7bctcAAAA6OjqEgpWrlyZD3zgA/nDH/6QY489Nr/5zW+yatWqrFq1Kj/72c8yYsSIrFixIp/4xCdeN/dnP/tZvvSlL6VLly65+uqrs2zZsqxYsSLPP/98zj333KxevTp/+7d/m0WLFtV0LgAAAFCnUNDY2JiTTz45hx9+eO6///684x3vSKdOndKlS5cce+yxuemmm5Ikf/jDH/L0009vMvdrX/takuSyyy7LP/zDP6R3795JkqFDh+b666/PCSeckGXLluW666573f22ZS4AAABQp1DQpUuXTJs2LQ8++GC6du36uu2HHXZY+bypqal8vnz58vzsZz9LY2Njzj///Fb3/alPfSpJ8oMf/GCT29syFwAAANigrhczbC0SJMmaNWvK53vuuWf5/Fe/+lXWrl2bQw45JH379m117lFHHZWGhoY8/vjjeeWVV2oyFwAAANigXd4e8Uc/+lGSZNSoUTn44IPL7S3vhLDPPvtsdm6PHj2y5557prm5OX/84x9rMhcAAADYoHFH3tnLL7+cu+66K//4j/+YXr165YYbbkinTv/VKubPn59kwxP6LenZs2cZ/7a3va3Nc7fVggULtrh945dRAAAAwK6o7qHgtttuy/nnn5/169dn+fLl6dq1a84444xcdNFFOfDAAzcZu3z58u3a98bj2zJ3W23pbAUAAADYHdT9pQdr1qzJ0qVLyxPz5ubmzJs373XvdpBseFvF7bFixYqazAUAAAA2qHsoOOuss1JVVQkE1113XZ544omMHTs206ZNq/fd19Rzzz23xY9HH320vZcIAAAAbbLDrlHQ0NCQ/fbbL5/85Cez//7756/+6q/yxS9+MR/5yEey3377Jdn69QVeq+V6A22du62GDRu23XMAAABgV9Iu73rwnve8J295y1uyfv363HXXXeX2Xr16bdd+Nh7flrkAAADABu0SCpJkxIgRSf7rbQ2TZN99902SrFq1aotzW64v0DK+rXMBAACADdotFLS8LWJzc3O5rSUePPfcc5udt3LlyixatCgNDQ3Zf//9azIXAAAA2KBuoeB//a//laqqNrt97ty5SZKhQ4eW24444oh06dIlTz75ZJYtW9bqvF/+8pepqiqHHnpoevfuXZO5AAAAwAZ1CQUf//jH85GPfCRXXnllq9t/+ctfZubMmUmS4447rtzeq1evHHfccVm7dm2uv/76VudeddVVSZIxY8Zscntb5gIAAAAb1CUUHHHEEUmSz3zmM/nCF76QF154IUmybNmy3H777TnttNNSVVXe/e53513vetcmcz/72c8mSS699NJcc801Wb58eZKkqakp48ePzw9/+MP06dMn48ePf939tmUuAAAAUKdQcM455+Sf//mf09DQkK9+9avZe++907179/Tt2zennnpqnn/++YwcOTK33Xbb6+Yed9xxufjii7NmzZpccMEF6dOnT3r27JmhQ4fmuuuuS9euXfOd73wne+65Z03nAgAAAHW8RsGnP/3pPPHEE/n4xz+eYcOGpbm5Ob17985hhx2Wr371q3nssccybNiwVud+6Utfyl133ZXjjz8+AwYMyPr167P33nvn9NNPz69//eucdNJJm73ftswFAACAjq6xnjs/6KCD8s1vfvMNzf3gBz+YD37wgzt8LgAAAHRk7fb2iAAAAMDORygAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgaGzvBQCw6xg+8b72XgIAAHXmjAIAAACgEAoAAACAwksPAADoEGrx8ql5006swUoAdm7OKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAIrG9l4AADvG8In3tfcSAADYBTijAAAAACiEAgAAAKDYIaHgsccey9ixYzNw4MB06dIle+21V0444YT8+Mc/3uK8hx9+OKecckqGDh2abt26ZciQIRk7dmweeuihrd5nW+YCAABAR1X3UHDppZfmsMMOy5133pkXX3wxnTp1yqJFi/LjH/84J5xwQi6//PJW5910000ZPXp07rjjjjQ1NaVz585ZuHBh7rzzzhx99NGZPn36Zu+zLXMBAACgI6trKPjXf/3XTJ06NV27ds0//dM/ZdGiRVm9enX++Mc/ZsyYMUmSz3/+85k7d+4m85555pmMHz8+zc3NueSSS7Jo0aKsWLEiixcvzqRJk9Lc3JwLLrggM2fOfN19tmUuAAAAdHR1CwXr1q3L5MmTkyRXXnllPvvZz2bAgAFJkje/+c353ve+l8GDB2fdunX53ve+t8ncK6+8MmvXrs24ceMyZcqUMq9///6ZPHlyzjvvvKxbt67VsxHaMhcAAAA6urqFgrVr12b8+PE58sgjM27cuNdt79GjR4455pgkye9+97tNtt1zzz1JkgsvvLDVfU+YMCFJcu+996aqqprNBQAAgI6ubqGgR48e+exnP5uHH344jY2NrY4ZNmxYkmTp0qXltrlz56apqSn9+vXLqFGjWp13wAEHZPDgwVm0aNEmLyFoy1wAAACgnd8ecc2aNUmSN73pTeW22bNnJ0n22WefLc5t2T5nzpyazAUAAACS1v/Uv4M8/fTTSTZcs6DF/Pnzk2w4I2FLevbsucn4ts7dFgsWLNji9qampu3aHwAAAOxs2i0UvPDCC5kxY0aS5AMf+EC5ffny5du1n43Ht2XuttjamQoAAACwq2u3lx585jOfybp16/Lud7877373u8vtK1eu3K79rFixoiZzAQAAgHY6o+Ab3/hGvve972WPPfbITTfd1B5LeEOee+65LW5vamrK4YcfvoNWAwAAALW3w0PBT37yk5x33nlpaGjIt771rYwcOXKT7Vu7vsBrtVxvoK1zt0XLuzQAAADA7mqHhoIHH3wwH/rQh7J27dpcccUVGTt27OvG9OrVa7v2ufH4tswF2NkNn3hfey8BAIAOYIddo+Cxxx7LBz7wgaxYsSKf//znM2HChFbH7bvvvkmSVatWbXF/LdcXaBnf1rkAAADADgoFTz31VN73vvdl2bJlOffcc/OVr3xls2NHjBiRZOvXA2jZ3jK+rXMBAACAHRAK5syZk/e+971ZvHhxPvzhD+faa6/d4vgRI0ZkyJAhWbJkSWbNmtXqmLlz52bhwoUZMGBARo0aVZO5AAAAQJ1Dwfz58/Oe97wnCxcuzPvf//58+9vfTqdOW7/Lk08+OUly9dVXt7r9qquuSpKcdNJJaWhoqNlcAAAA6OjqFgqWL1+ev/qrv8r8+fNz1FFH5Y477kiXLl22ae6ECRPSpUuXTJ8+PZMnT86SJUuSJC+99FKmTp2aa6+9No2NjbnoootqOhcAAAA6urqFgkWLFmX27NlJkt/+9rcZMmRI3vSmN23247bbbitzR44cmWuuuSadOnXKZZddlgEDBmSPPfZI//79c+mllyZJvv71r7f60oG2zAUAAICOboe8PeKrr7661TFr1qzZ5L/HjRuXt7zlLbniiivyq1/9KosXL86gQYNyxBFHZMKECTn66KM3u6+2zAUAAICOrG6hYPjw4amqqk37GD16dEaPHr3D5wIAAEBHtUPeHhEAAADYNQgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEDR2N4LAACAXcXwife1eR/zpp1Yg5UA1I8zCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKBrbewEAANCRDJ94X5vmz5t2Yo1WAtA6ZxQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQNHY3gsA6AiGT7yvvZcAAADbxBkFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQNHY3gsA2NkNn3hfey8BAAB2GGcUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABSN7b0AAABg2w2feF+b9zFv2ok1WAmwu3JGAQAAAFAIBQAAAEAhFAAAAADFDgsFn/vc59LQ0FA+pk+fvtU5Dz/8cE455ZQMHTo03bp1y5AhQzJ27Ng89NBDdZ0LAAAAHVXdQ0FVVRk/fny+9rWvpV+/funRo8c2zbvpppsyevTo3HHHHWlqakrnzp2zcOHC3HnnnTn66KO3GBraMhcAAAA6srqGgvXr1+fss8/Oddddl4EDB2bGjBkZOHDgVuc988wzGT9+fJqbm3PJJZdk0aJFWbFiRRYvXpxJkyalubk5F1xwQWbOnFnTuQAAANDR1S0UrF27NqeffnpuueWWDBs2LL/4xS/ytre9bZvmXnnllVm7dm3GjRuXKVOmZMCAAUmS/v37Z/LkyTnvvPOybt26XH755TWdCwAAAB1d3ULBrFmzcu+992b//ffPgw8+mJEjR27z3HvuuSdJcuGFF7a6fcKECUmSe++9N1VV1WwuAAAAdHR1CwWHHHJI7r777jz44IMZPnz4Ns+bO3dumpqa0q9fv4waNarVMQcccEAGDx6cRYsWbfISgrbMBQAAAOp8jYLjjz8+Q4cO3a45s2fPTpLss88+WxzXsn3OnDk1mQsAAAAkje29gNeaP39+kmz13RF69uy5yfi2zt0WCxYs2OL2pqam7dofAAAA7Gx2ulCwfPnyNzy+LXO3xdbOVAAAAIBd3U4XClauXLld41esWFGTuQAA0FEMn3hfm/cxb9qJNVgJsDPa6ULBzuy5557b4vampqYcfvjhO2g1AAAAUHs7XSjY2vUFXqvlegNtnbsthg0btl3jAQAAYFdT13c9eCN69er1hse3ZS4AAACwE4aCfffdN0myatWqLY5rub5Ay/i2zgUAAAB2wlAwYsSIJFu/HkDL9pbxbZ0LAAAA7KShYMiQIVmyZElmzZrV6pi5c+dm4cKFGTBgQEaNGlWTuQAAAMBOGAqS5OSTT06SXH311a1uv+qqq5IkJ510UhoaGmo2FwAAADq6nTIUTJgwIV26dMn06dMzefLkLFmyJEny0ksvZerUqbn22mvT2NiYiy66qKZzAQAAoKOrayh4z3vek8bGxk0+nn322STJ+PHjN7n9Pe95T5k3cuTIXHPNNenUqVMuu+yyDBgwIHvssUf69++fSy+9NEny9a9/vdWXDrRlLgAAAHR0dQ0F69evf91Hi+bm5s1uS5Jx48ZlxowZGTNmTAYPHpy1a9dm0KBBGTNmTGbMmJFzzz13s/fblrkAAADQkTXWc+c///nP2zR/9OjRGT169A6fCwAAAB3VTnmNAgAAAKB9CAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFA0tvcCAOpp+MT72nsJAACwS3FGAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAAhVAAAAAAFEIBAAAAUAgFAAAAQCEUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUje29AAAAYNczfOJ9bd7HvGkn1mAlQK0JBcBOqxa/gAAAANvHSw8AAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKBobO8FAAAAHdPwife1eR/zpp1Yg5UAGxMKgLqoxQ9+AABgx/PSAwAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoGhs7wUAAAC8UcMn3tfmfcybdmINVgK7D2cUAAAAAIVQAAAAABRCAQAAAFAIBQAAAEAhFAAAAACFUAAAAAAUQgEAAABQCAUAAABAIRQAAAAARWN7LwCoreET72vzPuZNO7EGKwEAAHZFzigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgMLFDIHXqcUFEQEAgF2TMwoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACga23sBAAAA7Wn4xPvavI95006swUpg5+CMAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAoXMwQdiK1uJAOAABAWzijAAAAACiEAgAAAKAQCgAAAIDCNQqgRlxfAAAA2B04owAAAAAonFEAAADQRrU4u3TetBNrsBJoO6EAAABgJyA2sLPw0gMAAACgcEYBAADAbsJZCdSCMwoAAACAQigAAAAAit06FDz88MM55ZRTMnTo0HTr1i1DhgzJ2LFj89BDD7X30gAAAGCntNuGgptuuimjR4/OHXfckaampnTu3DkLFy7MnXfemaOPPjrTp09v7yUCAADATme3DAXPPPNMxo8fn+bm5lxyySVZtGhRVqxYkcWLF2fSpElpbm7OBRdckJkzZ7b3UgEAAGCnsluGgiuvvDJr167NuHHjMmXKlAwYMCBJ0r9//0yePDnnnXde1q1bl8svv7ydVwoAAAA7l93y7RHvueeeJMmFF17Y6vYJEybk+uuvz7333puqqtLQ0LAjl8dOqBZvIwMAALA72O1Cwdy5c9PU1JR+/fpl1KhRrY454IADMnjw4CxcuDAzZ87MQQcdtINXufvwBBsAAHittj5PmDftxBqthDdit3vpwezZs5Mk++yzzxbHtWyfM2dO3dcEAAAAu4rd7oyC+fPnJ0l69OixxXE9e/bcZDwAAAA7x1nDO8Mako57ZsNuFwqWL19et/ELFizY4vbnnnuufN7U1LRd62gPR3zlP9t7CQBQrFm1LK/9Sbtm+eKsW7e2XdYDAFt7Drgz2Pi557p162qyz90uFKxcuXK7xq9YsWKbx27t5QwbO/zww7drHQDQ0T2f5HU/af/1H9phJQCwwT7Xt/cKts+LL76Y4cOHt3k/u901CgAAAIA3brc7o2Br1yZ4rZZrFWyLjV9a0JpVq1Zl1qxZGTRoUPbaa680Ntb+f29TU1M5W+HRRx/NkCFDan4f1Idjt2ty3HZdjt2uyXHbNTluuy7HbtfkuO266nHs1q1blxdffDFJcsghh7R5f8luGAp69epVt/HDhg3b6pgRI0Zs1/23xZAhQ7ZpTex8HLtdk+O263Lsdk2O267Jcdt1OXa7Jsdt11XLY1eLlxtsbLd76cG+++6bZMNf97ek5doELeMBAACA3TAUtPxFf2svE2jZviPPAAAAAICd3W4ZCoYMGZIlS5Zk1qxZrY6ZO3duFi5cmAEDBmTUqFE7eIUAAACw89rtQkGSnHzyyUmSq6++utXtV111VZLkpJNOSkNDw45aFgAAAOz0dstQMGHChHTp0iXTp0/P5MmTs2TJkiTJSy+9lKlTp+baa69NY2NjLrroonZeKQAAAOxcdstQMHLkyFxzzTXp1KlTLrvssgwYMCB77LFH+vfvn0svvTRJ8vWvf93LDgAAAOA1dstQkCTjxo3LjBkzMmbMmAwePDhr167NoEGDMmbMmMyYMSPnnntuey8RAAAAdjoNVVVV7b0IAAAAYOew255RAAAAAGw/oQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhoJ09/PDDOeWUUzJ06NB069YtQ4YMydixY/PQQw+1ed9/+MMfcuCBB6ahoSH//b//9zz//PM1WDEtanXsJk+enIaGhm3+oG1q/TU3e/bsjB8/PqNGjcoee+yRrl27ZujQoTnppJNy991313bxHVytj92aNWvyL//yLznyyCPTr1+/dO/ePcOHD8/f/d3fZc6cOTVefcf2uc99bpPvY9OnT2/T/r797W/nmGOOyZ577pnu3btn//33z3nnnZf58+fXaMW0qPWxW758eY477rhN9jlr1qwarZYWtTxuVVXl+9//fkaPHp0+ffqkW7du2W+//fKxj30sM2fOrOGqqeVxe/bZZ/OpT30qb37zm9O1a9f06tUr73jHOzJ58uQsW7ashqsmqf33yo395Cc/2WTfP//5z2u27y2qaDc33nhj1alTpypJlaTq0aNH+bxTp07V9ddf/4b3/fvf/74aOHBglaQ69dRTq1WrVtVw5dTy2E2aNKlKUnXr1q3q27dvqx+9evUq++eNq/XX3A9/+MOqW7duZR8NDQ1V165dy38nqc4666w6PZqOpdbH7uWXX66OOOKITfax8T67d+9e3X///XV6NB1Hc3Nzdf7551dJqn79+pX/x235+fbRj350k6+57t27l//u06dP9etf/7qGj6Djqsexe/nll6v/8T/+R5Wk2nvvvctxe/rpp2u48o6t1sdt1apV1fHHH7/J11yXLl3Kf/fs2bN64IEHavwoOp5aH7cf//jHm/zuuPHvKkmqESNGVP/3//7fGj+Kjqke3ytf67DDDqs6depUNTQ0VEmqGTNm1GzfW+JZRzuZNWtW+UZ7ySWXVIsWLaqqqqoWL15cnjg2NjZWv//977d73wsXLqyGDRtWJak+8YlPVM3NzbVefodW62PXMmfSpEmbHXPFFVdUSapDDz20Fg+hQ6r1cVu5cmW11157VUmqUaNGVTNmzKhWrVpVNTc3V3/605+qT3ziE+UH8p133lnPh7bbq8f3y4997GNVkmqvvfaq7rjjjmr16tVVVVXVs88+W5199tnlB75fpN64devWVWeddVaVpBo4cGD1+OOPV/vtt1+bfoH65je/WYLArbfeWq1cubKqqqqaM2dO9aEPfahKUg0fPrxas2ZNLR9Kh1OPY/fiiy9Wb3/726sk1RFHHFH99re/FQpqrB7H7cwzzyxh5/bbb69effXVau3atdUjjzxSjRo1qkpS/cVf/EW1fv36Gj+ajqPWx23mzJnlyer5559fzZ8/v6qqqlq+fHl1yy23lG3jxo2r9UPpcOrxNfdad911V5WkuvDCC8u+hYLd3DnnnLPFL9LzzjuvSlJ9/OMf3+59v+c976mSVB/60Id8466DWh+7bQkFLb9cXX755W9kyVS1P27/+Z//WX7JfeKJJ163vbm5uXrHO95RJanOPPPMNq29o6v1sZs/f37VuXPnLf6wbfkL2iWXXPJGl92hrVmzpjrttNOqJNWwYcOqWbNmVVVVtfkXqJEjR1ZJqttuu63V+3zrW99aJam+9a1vtWn9HVk9jt0LL7xQnlQec8wx1SuvvFL96U9/EgpqqB7HbdasWeUMgscff/x123/zm9+UY/jggw+2+TF0RPU4bn/+85+rE088sfqHf/iHVre3/N45cODANq29o6vXz7mNNTc3V4ccckjVvXv36oUXXhAKOoohQ4ZUSTb7F7A//OEPVZJqzz333K4zAlr+2rLffvtVr7zySq2Wy0Zqfey2Fgqeeuqpcmr0Cy+80Jald2i1Pm633XZblaTq3bv3Zsf8/d//fZWkOv7449/wuqn9sbvpppuqJNUBBxyw2TG33357laQ65JBD3vC6O7Innnii6tGjR7X//vtXf/rTn8rtbX3SkqQaNGjQZiP4jTfeWCWp/vqv//qNLr3Dq8exa/maO+GEE8pZIEJBbdXjuM2fP7+68MILq7PPPnuzYwYPHlwlqa655po3suwOrx7HrcXmzqz64Q9/WL72/EHxjavnsWvR8rvmhRdeuMm+d1QoaAw73Ny5c9PU1JR+/fpl1KhRrY454IADMnjw4CxcuDAzZ87MQQcdtNX9rlmzJpdeemmS5Gtf+1p69epV03VTn2M3atSojB07drP7u/XWW5Mkxx13XIYMGdK2B9BB1eO4bc+xGDp06Hatl/9Sj2P39NNPJ0kOPPDAzY5p2fb0009n9erV6dat2xt8BB3TIYcckrvvvjsHH3xwzf79t1y08sgjj0ynTq1fi3n06NFJkl/84hc1uc+OqB7H7u/+7u/So0ePnHbaaenSpUtN9smm6nHc9tlnn1x11VVbHDNs2LAsXLgwS5curcl9djT1OG4tNve1tmbNmiRJv379Nvu9lK2r57FLknXr1mXSpEnp3r17Pve5z9V8/9vCv452MHv27CQbvgFvScv2bb369m233ZYFCxbkXe96V0477bS2LZJW1ePYnXbaabn99ttbPWbNzc357ne/myT527/92+1dLv9fPY7b29/+9vTt2zevvPJKnnzyyddtr6oqjzzySJLkmGOO2c4V06Iex67las89evTY7JiePXsm2fCDevHixdu0VjZ1/PHH1/SXp235t9CybdGiRZ64tEGtj12SnHHGGSJBndXjuG1Ny5PON73pTTv0fncnO/q4/ehHP0qS/M3f/M0Ou8/dVT2P3be+9a3Mnj0755xzTrv9oVAoaActb9+0pV9Sk//6RXVb3+7pW9/6VpJkwoQJefDBB3P88cdn0KBBedOb3pQjjzwy3/zmN1NVVRtWTr2O3ebMmDEjCxYsSI8ePfLXf/3XbdpXR1aP49a3b9986UtfSpJ8+MMfzgMPPJA1a9akqqo8++yzGTduXH7zm9/kiCOOEHnaoB7HruXsgFdffXWzY5YvX14+X7ly5Vb3Sf1ty7+Fln8HG48H6mPNmjX54x//mCR585vf3M6rYUuqqsoLL7yQiy++ODfffHMOOOCATJ06tb2XxWasWbMmU6ZMadezCZLESw/awca/gNZq/KJFi/KLX/wie+yxRzp16pTjjjsu69atK9sfeeSRPPLII5kxY0Y5lZ3tV49jtyUtx+rkk09O796927Svjqxex+3v//7vM2jQoEybNi3HHHNMGhoa0qVLl6xZsyb9+/fPRRddlMmTJ6dz585vZNmkPseu5a/OTz311GbHbLytT58+27UG6mNHf/8FtuwHP/hBli9fnj59+pSX/bBzOf/883Pbbbdl9erVWbVqVQYOHJiJEyfmc5/7nLNAdmI33nhj5s+fnwsvvLBdX3bsjIJ2sL1/nVqxYsVWx/zqV79Kc3NzDj744EyaNCk33HBDXnzxxaxevTq/+93vMnbs2CTJd77znXIqO9uvHsduS3PvuOOOJF520Fb1PG6rV68uTySrqiqnYfbo0SPdu3dv078B6nPsjjrqqCTJvHnz8uMf//h126uqyvXXX58k6dy5c/r167dda6A+duT3X2DLVq5cmYkTJyZJLrzwwuyxxx7tvCJas2LFiixdujSrVq1KkqxatSpz5szJvHnz2ndhbNaKFSvy5S9/ud3PJkiEgt3G7373uyQbgsE//dM/5eyzz86ee+6Zrl275q1vfWu+//3vl9p73XXXtedS2UZ33313li9fnj333DPve9/72ns5vEZVVTnzzDPz0Y9+NM8//3xuv/32PP/881m2bFl+/etf513velemTp2aww47LHPnzm3v5bKRd73rXTn00EOTJB/5yEfyzW9+M4sXL87atWvz5JNP5tRTT83gwYOTJAcffHAaG518B7Cxc845J3/84x9z8MEH5wtf+EJ7L4fNuOWWW1JVVVavXp2ZM2fmvPPOy913350jjjgiP//5z9t7ebTimmuuycKFC9v12gQthIJ2sLXX2r7Wxq+53JxFixYl2fAasfe///2v297Q0JDzzz8/SfLYY49l7dq127UGNqjHsducb3/720k2XGzGE5W2qcdxu+2223LrrbemX79+efDBBzN27NgMHTo0vXv3zjvf+c787//9v3Paaadl3rx55WuP7Vevr7nvfOc7GThwYJYuXZpPfOITm4TVF198MQcffHCSDVGBncOO/P4LbN6kSZNy6623pn///rnrrrvSvXv39l4SW9G1a9cceOCBmTZtWiZOnJjVq1dn3Lhxrl22k1m2bFm+9rWv7RRnEyRCQbvY3rct3JbxLVfx3tLbgrVsW7NmTf785z9v1xrYoB7HrjULFy7Mf/zHfyTxsoNaqMdx+7d/+7ckG/4iPWjQoFbHfPrTn06S/PSnP81LL720XWtgg3p9zR100EF5/PHH86lPfSoHHnhgBg4cmEMPPTSXXXZZfvSjH+XOO+9MkpxyyinbvWbqY0d9/wU274orrsiUKVPSo0eP3HPPPRkxYkR7L4nt9Pd///dJNryTTMsZyewcrrzyyixevHinOJsgcTHDdrHvvvsmSXm90Oa0vL6yZfyWtPylZUuvEdv4l6at3Tetq8exa81tt92W9evX5y/+4i9yxBFHvKF98F/qcdxa3qrtgAMO2OyYkSNHJtnwMoU5c+bksMMO26b18l/q+TU3ZMiQXHnlla+7/Ze//GUef/zx7Lfffjn66KO3Y7XU07b8W9j4ugRv9Psv0Lobb7wx//iP/5jGxsZ8//vfd8bVLmrQoEHp1atXli9fntmzZ5eX4tG+lixZkiuuuGKnOZsgEQraRUt9fe6557Y4rmX7ttTalr9obuliTxu/FVjfvn23uk9erx7HrjUt73ZwxhlnvKH5bKoex239+vVJkk6dNn9iVkNDw+vGs3121Nfcxr74xS8mST73uc9t8fiyY23Lv4WWbXvuuaefc1BD3/3ud3PeeeeloaEh3/rWt3LiiSe295Jog5afbc3Nze28Elo88cQTWbZsWTp37pwDDzyw1TEtZ5B/4AMfSGNjY4466qj88Ic/rNua/AbUDkaMGJEhQ4ZkyZIlmTVrVqtj5s6dm4ULF2bAgAEZNWrUVvf59re/PUkyc+bMzY55+umnkyQDBw7Mnnvu+QZWTj2O3Ws99dRTefzxx5N42UGt1OO4tZwS1nJmQWv+8Ic/lM+HDh26nasm2TFfcxu7+eabM2PGjPzlX/5lzj777Dbti9p697vfnWTD2/1u7nW1Dz300CZjgba7++67c9ZZZ6W5uTnXXHNNTj/99PZeEluwYsWK/OAHP9js9kWLFpUnnH432fmsX78+S5cubfWj5Wffq6++mqVLl9b9bYCFgnZy8sknJ0muvvrqVrdfddVVSZKTTjppk79Kbs5xxx2X3r17Z86cOZu9iulNN92UJK6g30a1Pnav1XIRw8MPP3yLp7WzfWp93I499tgkG65VsGTJklbHXHvttUmS/fff32nQbVDvr7kWjz76aC644II0NDRk+vTp6dat2xveF7U3cuTIjBw5Mk1NTeUaIRtbt25drrnmmiTJmDFjdvTyYLf0k5/8JB/+8Iezbt26TJ061cV5d3LLli3LkUcemVNPPTW/+tWvWh3T8nygd+/eXhK5EznmmGNSVdUWP/bbb78kyYwZM1JVVf3fuaKiXcyaNavq0qVLlaSaNGlStXjx4qqqqmrJkiXVlClTqoaGhqqxsbH6/e9/X+a8/PLL1dvf/vaqe/fu1Uc/+tHX7fPSSy+tklQDBw6s7rjjjmr16tVVVVXV/Pnzq7PPPrtKUnXp0qV64okndsyD3E3V49i1WL9+fbX33ntXSap/+Zd/qftj6UhqfdwWLlxY9e3bt0pSHXroodV//ud/Vq+++mq1fv36as6cOdU555xTJamSVDfddNMOfay7m3p+zbX47W9/W/Xv379KUk2ZMqVuj6Wj22+//aok1fXXX9/q9q0dt2984xtVkqpv377Vd7/73WrVqlVVVVXV3Llzq7Fjx1ZJqv3226/8/KN22nrsXutPf/pT+R759NNP13q5/H9tOW6//e1vq549e1ZJqk9/+tM7Yrn8f205bh/+8IerJFW/fv2qb3zjG9XSpUurqqqqpqam6qtf/WrVvXv3Kkl1ySWX1P1xdES1/l7Z2r5nzJhRg5VunVDQjm644YaqU6dO5QdlyzfjJFVDQ8Pr/oE9+OCDZXtjY+Pr9rd69erqxBNPLGM6d+68yT47depUfeMb39hRD2+3Vutj1+KnP/1pGfPnP/+53g+jw6n1cfv3f//3qnfv3pvso+UJbcvH+PHjd9TD263V62uuqqrqZz/7WYk+55xzTtXc3FzPh9JhHHfccVXnzp03+dj459HGtx933HFVVW3bcTvjjDM22U/LL71Jqj59+lSPPvrojnyYu6V6HLvX7m/jfb729ssuu2xHP+TdQq2P27/+67+WbX379t3qx7PPPtseD3uXV+vjtmrVqupv/uZvNvldpFu3bpv895gxY6q1a9e2x8PdrdTr59zmCAUdzAMPPFCNGTOmGjx4cNWlS5dq0KBB1ZgxY6qf//znrxv78ssvV4ceemjVrVu3zRao9evXVzfffHN11FFHVX379q26du1aDRs2rDr99NOr3/zmN/V+OB1KrY9dVVXVmWeeWSWpTjjhhHouvUOr9XGbN29e9alPfao66KCDqj322KPq0qVLNXTo0OpDH/pQdf/999f74XQo9fiau/HGG0vcufjii0WCGjr66KM3+cV0Sx9HH310VVXbfty+8Y1vVEcddVTVr1+/qlu3btV+++1XffKTn6z+9Kc/7ZgHt5urx7Hb1v21nDnE9qv1cds4FGzLh6+/N6Ze3yt//OMfVyeccELVv3//qrGxsRowYED13ve+t/rOd77jZ12N1PPnXGt2dChoqKrNXBEIAAAA6HBczBAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAohAIAAACgEAoAAACAQigAAAAACqEAAAAAKIQCAAAAoBAKAAAAgEIoAAAAAAqhAAAAACiEAgAAAKAQCgAAAIBCKAAAAAAKoQAAAAAo/h+Wc3yB0VFxfQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 600x600 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 505,
       "width": 517
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sim_cluster(50, 8, 5, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### application to Penn bonus experiment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.0899064590011073, 0.030797749686469634)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "est, var = neyman_SRE(z, y, block)\n",
    "est, np.sqrt(var)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## chong et al (2016)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>class_level</th>\n",
       "      <th>1.0</th>\n",
       "      <th>2.0</th>\n",
       "      <th>3.0</th>\n",
       "      <th>4.0</th>\n",
       "      <th>5.0</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>treatment</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Soccer Player</th>\n",
       "      <td>16</td>\n",
       "      <td>19</td>\n",
       "      <td>15</td>\n",
       "      <td>10</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Physician</th>\n",
       "      <td>17</td>\n",
       "      <td>20</td>\n",
       "      <td>15</td>\n",
       "      <td>11</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Placebo</th>\n",
       "      <td>15</td>\n",
       "      <td>19</td>\n",
       "      <td>16</td>\n",
       "      <td>12</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "class_level    1.0  2.0  3.0  4.0  5.0\n",
       "treatment                             \n",
       "Soccer Player   16   19   15   10   10\n",
       "Physician       17   20   15   11   10\n",
       "Placebo         15   19   16   12   10"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dat_chong = pd.read_stata(\"chong.dta\")\n",
    "pd.crosstab(dat_chong.treatment, dat_chong.class_level)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>class_level</th>\n",
       "      <th>1.0</th>\n",
       "      <th>2.0</th>\n",
       "      <th>3.0</th>\n",
       "      <th>4.0</th>\n",
       "      <th>5.0</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>z</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>15</td>\n",
       "      <td>19</td>\n",
       "      <td>16</td>\n",
       "      <td>12</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>17</td>\n",
       "      <td>20</td>\n",
       "      <td>15</td>\n",
       "      <td>11</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "class_level  1.0  2.0  3.0  4.0  5.0\n",
       "z                                   \n",
       "0             15   19   16   12   10\n",
       "1             17   20   15   11   10"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "use_vars = [\"treatment\", \"gradesq34\", \"class_level\", \"anemic_base_re\"]\n",
    "dat_physician = dat_chong.loc[dat_chong.treatment != \"Soccer Player\", use_vars]\n",
    "dat_physician[\"z\"] = np.where(dat_physician.treatment == \"Physician\", 1, 0)\n",
    "dat_physician[\"y\"] = dat_physician.gradesq34\n",
    "pd.crosstab(dat_physician.z, dat_physician.class_level)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.40589046478271484, 0.04096197815071462)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(\n",
    "    tauS := neyman_SRE(\n",
    "        dat_physician.z.values, dat_physician.y.values, dat_physician.class_level.values\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.4633431335975384, 0.03624964630443229)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dat_physician[\"sps\"] = pd.Categorical(\n",
    "    dat_physician[\"class_level\"].astype(str)\n",
    "    + \"_\"\n",
    "    + dat_physician[\"anemic_base_re\"].astype(str)\n",
    ")\n",
    "(\n",
    "    tauSPS := neyman_SRE(\n",
    "        dat_physician.z.values, dat_physician.y.values, dat_physician.sps.values\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Estimate</th>\n",
       "      <th>Std. Error</th>\n",
       "      <th>z value</th>\n",
       "      <th>Pr(&gt;|z|)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Stratify</th>\n",
       "      <td>0.405890</td>\n",
       "      <td>0.202391</td>\n",
       "      <td>2.005480</td>\n",
       "      <td>0.044912</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Stratify and post-stratify</th>\n",
       "      <td>0.463343</td>\n",
       "      <td>0.190393</td>\n",
       "      <td>2.433609</td>\n",
       "      <td>0.014949</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                            Estimate  Std. Error   z value  Pr(>|z|)\n",
       "Stratify                    0.405890    0.202391  2.005480  0.044912\n",
       "Stratify and post-stratify  0.463343    0.190393  2.433609  0.014949"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "seS = np.sqrt(tauS[1])\n",
    "seSPS = np.sqrt(tauSPS[1])\n",
    "\n",
    "pvalS = 2 * (1 - sp.stats.norm.cdf(abs(tauS[0] / seS)))\n",
    "pvalSPS = 2 * (1 - sp.stats.norm.cdf(abs(tauSPS[0] / seSPS)))\n",
    "\n",
    "pd.DataFrame(\n",
    "    np.r_[\n",
    "        np.c_[tauS[0], seS, tauS[0] / seS, pvalS],\n",
    "        np.c_[tauSPS[0], seSPS, tauSPS[0] / seSPS, pvalSPS],\n",
    "    ],\n",
    "    columns=[\"Estimate\", \"Std. Error\", \"z value\", \"Pr(>|z|)\"],\n",
    "    index=[\"Stratify\", \"Stratify and post-stratify\"],\n",
    ")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "econometrics",
   "language": "python",
   "name": "econometrics"
  },
  "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}