diff --git a/pyfixest/Chapter01CorrAssocSimpsons.ipynb b/pyfixest/Chapter01CorrAssocSimpsons.ipynb
new file mode 100644
index 0000000..e58860e
--- /dev/null
+++ b/pyfixest/Chapter01CorrAssocSimpsons.ipynb
@@ -0,0 +1,491 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Chapter 1: Correlation, Association, and the Yule-Simpson Paradox"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "import pandas as pd\n",
+ "import scipy as sp\n",
+ "import pyfixest as pf\n",
+ "\n",
+ "# viz\n",
+ "from IPython.core.interactiveshell import InteractiveShell\n",
+ "\n",
+ "InteractiveShell.ast_node_interactivity = \"all\"\n",
+ "\n",
+ "\n",
+ "%load_ext autoreload\n",
+ "%autoreload 1\n",
+ "\n",
+ "%load_ext watermark\n",
+ "%watermark --iversions\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Unadjusted and Adjusted Regression\n",
+ "Lalonde Observational Data"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# read CPS data\n",
+ "dat = pd.read_table(\"cps1re74.csv\", delimiter=\" \")\n",
+ "dat[\"u74\"] = np.where(dat[\"re74\"] == 0, 1, 0)\n",
+ "dat[\"u75\"] = np.where(dat[\"re75\"] == 0, 1, 0)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Unadjusted regression"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "pf.feols(\"re78 ~ treat\", data=dat, vcov=\"HC2\").tidy().loc[\"treat\"]\n",
+ "# %%\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Adjusted Regression"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "rhs = list(set(dat.columns) - {\"re78\", \"treat\"})\n",
+ "pf.feols(f're78 ~ treat + {\"+\".join(rhs)}', data=dat, vcov=\"hetero\").tidy().loc[\n",
+ " \"treat\"\n",
+ "]\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Fisher's exact test for contingency tables\n",
+ "Bertrand and Mullainathan (2004) experiment "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "call 0 1\n",
+ "race \n",
+ "black 2278 157\n",
+ "white 2200 235\n"
+ ]
+ }
+ ],
+ "source": [
+ "resume = pd.read_csv(\"resume.csv\")\n",
+ "print(xtab := pd.crosstab(resume[\"race\"], resume[\"call\"]))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "SignificanceResult(statistic=1.5498841922408801, pvalue=4.758747107909523e-05)"
+ ]
+ },
+ "execution_count": 6,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "sp.stats.fisher_exact(xtab)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## simpson's paradox \n",
+ "UCB admissions data"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " Admit | \n",
+ " Gender | \n",
+ " Dept | \n",
+ " Freq | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | 0 | \n",
+ " Admitted | \n",
+ " Male | \n",
+ " A | \n",
+ " 512 | \n",
+ "
\n",
+ " \n",
+ " | 1 | \n",
+ " Rejected | \n",
+ " Male | \n",
+ " A | \n",
+ " 313 | \n",
+ "
\n",
+ " \n",
+ " | 2 | \n",
+ " Admitted | \n",
+ " Female | \n",
+ " A | \n",
+ " 89 | \n",
+ "
\n",
+ " \n",
+ " | 3 | \n",
+ " Rejected | \n",
+ " Female | \n",
+ " A | \n",
+ " 19 | \n",
+ "
\n",
+ " \n",
+ " | 4 | \n",
+ " Admitted | \n",
+ " Male | \n",
+ " B | \n",
+ " 353 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | Admit | \n",
+ " Admitted | \n",
+ " Rejected | \n",
+ "
\n",
+ " \n",
+ " | Gender | \n",
+ " | \n",
+ " | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | Female | \n",
+ " 557 | \n",
+ " 1278 | \n",
+ "
\n",
+ " \n",
+ " | Male | \n",
+ " 1198 | \n",
+ " 1493 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " Admit | \n",
+ " Gender | \n",
+ " Dept | \n",
+ " Freq | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | 0 | \n",
+ " Admitted | \n",
+ " Male | \n",
+ " A | \n",
+ " 512 | \n",
+ "
\n",
+ " \n",
+ " | 1 | \n",
+ " Rejected | \n",
+ " Male | \n",
+ " A | \n",
+ " 313 | \n",
+ "
\n",
+ " \n",
+ " | 2 | \n",
+ " Admitted | \n",
+ " Female | \n",
+ " A | \n",
+ " 89 | \n",
+ "
\n",
+ " \n",
+ " | 3 | \n",
+ " Rejected | \n",
+ " Female | \n",
+ " A | \n",
+ " 19 | \n",
+ "
\n",
+ " \n",
+ " | 4 | \n",
+ " Admitted | \n",
+ " Male | \n",
+ " B | \n",
+ " 353 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
"
+ ]
+ },
+ "metadata": {
+ "image/png": {
+ "height": 1241,
+ "width": 691
+ }
+ },
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "f, ax = plt.subplots(3, 2, figsize=(8, 15))\n",
+ "for i, (y0, y1, res, v) in enumerate(\n",
+ " zip([y0_1, y0_2, y0_3], [y1_1, y1_2, y1_3], [res1, res2, res3], [v1, v2, v3])\n",
+ "):\n",
+ " ax[i][0].plot(y0, y1, \"o\")\n",
+ " ax[i][0].set_xlabel(r\"$Y^0$\", fontsize=12)\n",
+ " ax[i][0].set_ylabel(r\"$Y^1$\", fontsize=12)\n",
+ " ax[i][1].hist(res[:, 0] - tautrue, bins=50, density=True)\n",
+ " ax[i][1].set_xlabel(r\"$\\hat{\\tau} - \\tau$\", fontsize=12)\n",
+ " x = np.linspace(-1, 1, 100)\n",
+ " y = sp.stats.norm.pdf(x, loc=0, scale=np.sqrt(v))\n",
+ " ax[i][1].plot(x, y, label=\"Neyman variance\")\n",
+ "ax[0][1].set_title(r\"$\\tau_i = \\tau$\")\n",
+ "ax[1][1].set_title(r\"Negative Corr between $Y^1$ and $Y^0$\")\n",
+ "ax[2][1].set_title(r\"Uncorrelated $Y^1$ and $Y^0$\")\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " Constant | \n",
+ " Negative | \n",
+ " Independent | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | Variance | \n",
+ " 0.031073 | \n",
+ " 0.006261 | \n",
+ " 0.019115 | \n",
+ "
\n",
+ " \n",
+ " | Estimated Variance | \n",
+ " 0.030088 | \n",
+ " 0.030094 | \n",
+ " 0.030080 | \n",
+ "
\n",
+ " \n",
+ " | Coverage | \n",
+ " 0.941100 | \n",
+ " 0.999900 | \n",
+ " 0.982800 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
"
+ ]
+ },
+ "metadata": {
+ "image/png": {
+ "height": 836,
+ "width": 446
+ }
+ },
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "n, n1, n0 = 1000, 600, 400\n",
+ "tautrue = 1\n",
+ "\n",
+ "# cauchy contamination: prob combination = 0.1\n",
+ "eps = np.random.binomial(1, 0.1, n)\n",
+ "y0_1 = (1 - eps) * np.random.exponential(1, n) + eps * np.random.standard_cauchy(n)\n",
+ "y1_1 = y0_1 + tautrue\n",
+ "tautrue_1 = np.mean(y1_1) - np.mean(y0_1)\n",
+ "z_1 = z_2 = z_3 = np.repeat([0, 1], [n0, n1])\n",
+ "MC = int(1e4)\n",
+ "res1 = np.zeros((MC, 3))\n",
+ "for i in range(MC):\n",
+ " res1[i] = sampling_dist(z_1, y1_1, y0_1, n=n, n1=n1, n0=n0, tautrue=1)\n",
+ "v1 = np.var(y1_1) / n1 + np.var(y0_1) / n0 - np.var(y1_1 - y0_1) / n\n",
+ "v1\n",
+ "\n",
+ "# cauchy contamination: prob combination = 0.3\n",
+ "eps = np.random.binomial(1, 0.3, n)\n",
+ "y0_2 = (1 - eps) * np.random.exponential(1, n) + eps * np.random.standard_cauchy(n)\n",
+ "y1_2 = y0_2 + tautrue\n",
+ "tautrue_2 = np.mean(y1_2) - np.mean(y0_2)\n",
+ "res2 = np.zeros((MC, 3))\n",
+ "for i in range(MC):\n",
+ " res2[i] = sampling_dist(z_2, y1_2, y0_2, n=n, n1=n1, n0=n0, tautrue=1)\n",
+ "v2 = np.var(y1_2) / n1 + np.var(y0_2) / n0 - np.var(y1_2 - y0_2) / n\n",
+ "v2\n",
+ "\n",
+ "\n",
+ "# cauchy contamination: prob combination = 0.5\n",
+ "eps = np.random.binomial(1, 0.5, n)\n",
+ "y0_3 = (1 - eps) * np.random.exponential(1, n) + eps * np.random.standard_cauchy(n)\n",
+ "y1_3 = y0_3 + tautrue\n",
+ "tautrue_3 = np.mean(y1_3) - np.mean(y0_3)\n",
+ "res3 = np.zeros((MC, 3))\n",
+ "for i in range(MC):\n",
+ " res3[i] = sampling_dist(z_3, y1_3, y0_3, n=n, n1=n1, n0=n0, tautrue=1)\n",
+ "v3 = np.var(y1_3) / n1 + np.var(y0_3) / n0 - np.var(y1_3 - y0_3) / n\n",
+ "\n",
+ "v1, v2, v3\n",
+ "\n",
+ "f, ax = plt.subplots(3, 1, figsize=(5, 10))\n",
+ "for i, (res, v) in enumerate(zip([res1, res2, res3], [v1, v2, v3])):\n",
+ " ax[i].hist(res[:, 0] - tautrue, bins=100, density=True)\n",
+ " ax[i].set_xlabel(r\"$\\hat{\\tau} - \\tau$\", fontsize=12)\n",
+ " x = np.linspace(-10, 10, 1000)\n",
+ " y = sp.stats.norm.pdf(x, loc=0, scale=np.sqrt(v))\n",
+ " ax[i].plot(x, y)\n",
+ " ax[i].set_xlim(-4, 4)\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Application"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from causalinference.utils import lalonde_data\n",
+ "\n",
+ "y, z, _ = lalonde_data()\n",
+ "y *= 1000"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "(1794.3430782016621, 670.9967300240978)"
+ ]
+ },
+ "execution_count": 9,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "n0, n1 = np.sum(1 - z), np.sum(z)\n",
+ "tauhat = np.mean(y[z == 1]) - np.mean(y[z == 0])\n",
+ "vhat = np.var(y[z == 1], ddof=1) / n1 + np.var(y[z == 0], ddof=1) / n0\n",
+ "sehat = np.sqrt(vhat)\n",
+ "tauhat, sehat"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Homoskedastic"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "reg_data = pd.DataFrame({\"y\": y, \"z\": z})\n",
+ "pf.feols(\"y ~ z\", data=reg_data).tidy()\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Heteroskedastic"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# PyFixest's heteroskedastic covariance is HC1.\n",
+ "reg_data = pd.DataFrame({\"y\": y, \"z\": z})\n",
+ "pf.feols(\"y ~ z\", data=reg_data, vcov=\"hetero\").tidy()\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "reg_data = pd.DataFrame({\"y\": y, \"z\": z})\n",
+ "pf.feols(\"y ~ z\", data=reg_data, vcov=\"HC2\").tidy()\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
+}
diff --git a/pyfixest/Chapter05StratandPostStrat.ipynb b/pyfixest/Chapter05StratandPostStrat.ipynb
new file mode 100644
index 0000000..2294b9e
--- /dev/null
+++ b/pyfixest/Chapter05StratandPostStrat.ipynb
@@ -0,0 +1,847 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Chapter 5: Stratification and Post-Stratification in Randomized Experiments"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "import pandas as pd\n",
+ "import scipy as sp\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\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " duration | \n",
+ " treatment | \n",
+ " female | \n",
+ " black | \n",
+ " hispanic | \n",
+ " ndependents | \n",
+ " recall | \n",
+ " young | \n",
+ " old | \n",
+ " quarter | \n",
+ " durable | \n",
+ " lusd | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | 1 | \n",
+ " 18.011343 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 2 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 5 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ "
\n",
+ " \n",
+ " | 4 | \n",
+ " 1.003399 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 5 | \n",
+ " 0 | \n",
+ " 1 | \n",
+ "
\n",
+ " \n",
+ " | 5 | \n",
+ " 26.960396 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 4 | \n",
+ " 0 | \n",
+ " 1 | \n",
+ "
\n",
+ " \n",
+ " | 6 | \n",
+ " 7.009044 | \n",
+ " 1 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 2 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ "
\n",
+ " \n",
+ " | 12 | \n",
+ " 9.022409 | \n",
+ " 1 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 1 | \n",
+ " 0 | \n",
+ " 3 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | col_0 | \n",
+ " 0 | \n",
+ " 1 | \n",
+ " 2 | \n",
+ " 3 | \n",
+ " 4 | \n",
+ " 5 | \n",
+ "
\n",
+ " \n",
+ " | row_0 | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | 0 | \n",
+ " 234 | \n",
+ " 41 | \n",
+ " 687 | \n",
+ " 794 | \n",
+ " 738 | \n",
+ " 860 | \n",
+ "
\n",
+ " \n",
+ " | 1 | \n",
+ " 87 | \n",
+ " 48 | \n",
+ " 757 | \n",
+ " 866 | \n",
+ " 811 | \n",
+ " 461 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
"
+ ]
+ },
+ "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": [
+ ""
+ ]
+ },
+ "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": [
+ ""
+ ]
+ },
+ "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": [
+ "\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | class_level | \n",
+ " 1.0 | \n",
+ " 2.0 | \n",
+ " 3.0 | \n",
+ " 4.0 | \n",
+ " 5.0 | \n",
+ "
\n",
+ " \n",
+ " | treatment | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | Soccer Player | \n",
+ " 16 | \n",
+ " 19 | \n",
+ " 15 | \n",
+ " 10 | \n",
+ " 10 | \n",
+ "
\n",
+ " \n",
+ " | Physician | \n",
+ " 17 | \n",
+ " 20 | \n",
+ " 15 | \n",
+ " 11 | \n",
+ " 10 | \n",
+ "
\n",
+ " \n",
+ " | Placebo | \n",
+ " 15 | \n",
+ " 19 | \n",
+ " 16 | \n",
+ " 12 | \n",
+ " 10 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | class_level | \n",
+ " 1.0 | \n",
+ " 2.0 | \n",
+ " 3.0 | \n",
+ " 4.0 | \n",
+ " 5.0 | \n",
+ "
\n",
+ " \n",
+ " | z | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | 0 | \n",
+ " 15 | \n",
+ " 19 | \n",
+ " 16 | \n",
+ " 12 | \n",
+ " 10 | \n",
+ "
\n",
+ " \n",
+ " | 1 | \n",
+ " 17 | \n",
+ " 20 | \n",
+ " 15 | \n",
+ " 11 | \n",
+ " 10 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " Estimate | \n",
+ " Std. Error | \n",
+ " z value | \n",
+ " Pr(>|z|) | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | Stratify | \n",
+ " 0.405890 | \n",
+ " 0.202391 | \n",
+ " 2.005480 | \n",
+ " 0.044912 | \n",
+ "
\n",
+ " \n",
+ " | Stratify and post-stratify | \n",
+ " 0.463343 | \n",
+ " 0.190393 | \n",
+ " 2.433609 | \n",
+ " 0.014949 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " coef | \n",
+ " se | \n",
+ " t | \n",
+ " p | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | unadjusted | \n",
+ " 0.0518 | \n",
+ " 0.078 | \n",
+ " 0.669 | \n",
+ " 0.504 | \n",
+ "
\n",
+ " \n",
+ " | adjusted | \n",
+ " 0.0682 | \n",
+ " 0.074 | \n",
+ " 0.925 | \n",
+ " 0.355 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "