{
"cells": [
{
"cell_type": "markdown",
"id": "1aaff07b-c682-4c12-92c7-40b06ff882d9",
"metadata": {},
"source": [
"(lecture_17)=\n",
"# Measurement and Misclassification\n",
":::{post} Jan 7, 2024\n",
":tags: statistical rethinking, bayesian inference, measurement error\n",
":category: intermediate\n",
":author: Dustin Stansbury\n",
":::\n",
"\n",
"This notebook is part of the PyMC port of the [Statistical Rethinking 2023](https://github.com/rmcelreath/stat_rethinking_2023) lecture series by Richard McElreath.\n",
"\n",
"[Video - Lecture 17 - Measurement and Misclassification](https://youtu.be/mt9WKbQJrI4)# [Lecture 17 - Measurement & Misclassification](https://www.youtube.com/watch?v=mt9WKbQJrI4)"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "015e420d-9c89-4c5f-8885-d7d619ba13a4",
"metadata": {},
"outputs": [],
"source": [
"# Ignore warnings\n",
"import warnings\n",
"\n",
"import arviz as az\n",
"import numpy as np\n",
"import pandas as pd\n",
"import pymc as pm\n",
"import pytensor.tensor as pt\n",
"import statsmodels.formula.api as smf\n",
"import utils as utils\n",
"import xarray as xr\n",
"\n",
"from matplotlib import pyplot as plt\n",
"from matplotlib import style\n",
"from scipy import stats as stats\n",
"\n",
"warnings.filterwarnings(\"ignore\")\n",
"\n",
"# Set matplotlib style\n",
"STYLE = \"statistical-rethinking-2023.mplstyle\"\n",
"style.use(STYLE)"
]
},
{
"cell_type": "markdown",
"id": "3829d401-d400-4b28-ad3e-4b87372e378e",
"metadata": {},
"source": [
"# Probability, Intuition, and Pancakes 🥞\n",
"\n",
"McElreath poses a re-dressing of the [classic Monty Hall Problem](https://en.wikipedia.org/wiki/Monty_Hall_problem). The original problem uses an old gameshow called \"Let's Make a Deal\" (hosted by Monty Hall) as the backdrop for a scenario where the correct, optimal strategy for winning a game is given by following the rules of probability theory. What's interesting about the Monty Hall problem is that the optimal strategy doesn't align with our intuitions.\n",
"\n",
"In lecure, instead of opening doors to find donkeys or prizes, as was in the case of the Monty Hall problem, we have pancakes that are either burnt or perfectly cooked on either side. The thought experiment goes like this:\n",
"\n",
"- You have 3 pancakes:\n",
" - One (likely the first one you cooked) is burnt on both sides\n",
" - One (likely the next one cooked, after you're improving your cooking skills) is burnt on only one side\n",
" - One is cooked perfectly on both sides\n",
"- Say you are given a pancake at random and it is burned on one side\n",
"- What is the probability that the other side is burned?\n",
"\n",
"Most folks would say 1/2, which intuitively _feels_ correct. However, the correct answer is given by Bayes rule:\n",
"\n",
"If we define $U, D$ as observing upside $U$ or downsid $D$ hot being burnt and $U', D'$ as up or down sides being burnt\n",
"\n",
"$$\n",
"\\begin{align*}\n",
"p(D' | U') &= \\frac{p(U', D')}{p(U')} &\\text{Bayes Rule} \\\\\n",
"p(U', D') &= 1/3 &\\text{only one out of three pancakes has both sides burnt} \\\\\n",
"p(U') &= (1/3 \\times 2/2) + (1/3 \\times 1/2) + (1/3 \\times 0) = 3/6 = 1/2 &\\text{probability of initially observing upside burnt} \\\\\n",
"p(D' | U') &= \\frac{1/3}{1/2} = 2/3\n",
"\\end{align*}\n",
"$$\n",
"\n",
"## Avoid being clever\n",
"\n",
"- Being clever is unreliable\n",
"- Being clever is opaque\n",
"- You'll rarely beat the axioms of probability theory\n",
"- Probability allows us to solve complex problems, if we stick to the rules\n"
]
},
{
"cell_type": "markdown",
"id": "d3f780fa-4a1c-4477-a4d5-77f4d5894ee2",
"metadata": {},
"source": [
"# Measurement Error\n",
"\n",
"- Many variables are proxies for what we want to measure (e.g. Descendant elemental confound)\n",
"- It's common practice to ignore measurement error\n",
" - assumes error is harmless, or will \"average out\"\n",
"- Better to draw out the measurement error, think causally, and fall back on probability theory (don't be clever)"
]
},
{
"cell_type": "markdown",
"id": "38a21cc1-4874-45d2-b1d4-8b7f0e3e7055",
"metadata": {},
"source": [
"## Myth: Measurement error can only decrease an effect, not increase\n",
"\n",
"- This is incorrect, measurement error has no predefined direction of causal effect\n",
"- It can, in some cases, _increase_ an effect, as we'll demonstrate below"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "6fd3cd90-b700-4326-8cc1-3bf8b533b202",
"metadata": {
"jupyter": {
"source_hidden": true
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"P \n",
"\n",
"P \n",
" \n",
"\n",
"\n",
"C \n",
"\n",
"C \n",
" \n",
"\n",
"\n",
"P->C \n",
" \n",
"\n",
"\n",
"P* \n",
"\n",
"P* \n",
" \n",
"\n",
"\n",
"P->P* \n",
" \n",
"\n",
"\n",
"e \n",
"\n",
"e \n",
" \n",
"\n",
"\n",
"C->e \n",
"recall bias \n",
" \n",
"\n",
"\n",
"e->P* \n",
" \n",
"\n",
"\n",
"unobserved \n",
"\n",
"unobserved \n",
" \n",
" \n",
" \n"
],
"text/plain": [
""
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"utils.draw_causal_graph(\n",
" edge_list=[(\"P\", \"C\"), (\"P\", \"P*\"), (\"e\", \"P*\"), (\"C\", \"e\")],\n",
" node_props={\n",
" \"P\": {\"style\": \"dashed\"},\n",
" \"unobserved\": {\"style\": \"dashed\"},\n",
" },\n",
" edge_props={(\"C\", \"e\"): {\"label\": \"recall bias\", \"color\": \"blue\", \"fontcolor\": \"blue\"}},\n",
" graph_direction=\"LR\",\n",
")"
]
},
{
"cell_type": "markdown",
"id": "b88467a6-b18d-450d-bb0f-46c322bd601c",
"metadata": {},
"source": [
"- Child Income $C$\n",
"- Actual Parental Income $P$, unobserved\n",
"- Measured parental income $P^*$\n",
"- Error in Parental income reporting $P^*$ (e.g. due to recall bias)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "d4a8038c-eaa5-4d74-a995-b71993281788",
"metadata": {},
"outputs": [],
"source": [
"def simulate_child_parent_income(beta_P=0, n_samples=500, title=None, random_seed=123):\n",
" np.random.seed(random_seed)\n",
" P = stats.norm().rvs(size=n_samples)\n",
" C = stats.norm.rvs(beta_P * P)\n",
" mu_P_star = 0.8 * P + 0.2 * C\n",
" P_star = stats.norm.rvs(mu_P_star)\n",
"\n",
" with pm.Model() as model:\n",
" sigma = pm.Exponential(\"sigma\", 1)\n",
" beta = pm.Normal(\"beta\", 0, 1)\n",
" alpha = pm.Normal(\"alpha\", 0, 1)\n",
" mu = alpha + beta * P_star\n",
" pm.Normal(\"C\", mu, sigma, observed=C)\n",
" inference = pm.sample()\n",
"\n",
" az.plot_dist(inference.posterior[\"beta\"])\n",
" plt.axvline(beta_P, color=\"k\", linestyle=\"--\", label=\"actual\")\n",
" plt.title(title)\n",
" plt.legend()\n",
"\n",
" return az.summary(inference, var_names=[\"alpha\", \"beta\", \"sigma\"])"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "a87b5077-cb1c-4647-81a5-001c1c7e7d2f",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (4 chains in 4 jobs)\n",
"NUTS: [sigma, beta, alpha]\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "c016896552b34b9fa2a0448127553a8a",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"Output()"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" mean \n",
" sd \n",
" hdi_3% \n",
" hdi_97% \n",
" mcse_mean \n",
" mcse_sd \n",
" ess_bulk \n",
" ess_tail \n",
" r_hat \n",
" \n",
" \n",
" \n",
" \n",
" alpha \n",
" -0.043 \n",
" 0.044 \n",
" -0.123 \n",
" 0.043 \n",
" 0.001 \n",
" 0.001 \n",
" 6736.0 \n",
" 2944.0 \n",
" 1.0 \n",
" \n",
" \n",
" beta \n",
" 0.165 \n",
" 0.038 \n",
" 0.091 \n",
" 0.235 \n",
" 0.000 \n",
" 0.000 \n",
" 6331.0 \n",
" 3019.0 \n",
" 1.0 \n",
" \n",
" \n",
" sigma \n",
" 0.982 \n",
" 0.032 \n",
" 0.925 \n",
" 1.046 \n",
" 0.000 \n",
" 0.000 \n",
" 6682.0 \n",
" 2723.0 \n",
" 1.0 \n",
" \n",
" \n",
"
\n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"simulate_child_parent_income(title=\"Measurement error makes\\nus over-estimate effect\")"
]
},
{
"cell_type": "markdown",
"id": "b64021fe-9c57-41eb-a11e-ca7591efbf5e",
"metadata": {},
"source": [
"In the scenario above, measurement error can cause us over-estimate the actual effect"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "4a57fd6a-cc8f-4a61-b8a2-0588b1a54883",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (4 chains in 4 jobs)\n",
"NUTS: [sigma, beta, alpha]\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "4ea027e87cf44c0b934e262413e50740",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"Output()"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" mean \n",
" sd \n",
" hdi_3% \n",
" hdi_97% \n",
" mcse_mean \n",
" mcse_sd \n",
" ess_bulk \n",
" ess_tail \n",
" r_hat \n",
" \n",
" \n",
" \n",
" \n",
" alpha \n",
" -0.077 \n",
" 0.047 \n",
" -0.163 \n",
" 0.014 \n",
" 0.001 \n",
" 0.001 \n",
" 5738.0 \n",
" 3032.0 \n",
" 1.0 \n",
" \n",
" \n",
" beta \n",
" 0.541 \n",
" 0.036 \n",
" 0.474 \n",
" 0.608 \n",
" 0.000 \n",
" 0.000 \n",
" 6257.0 \n",
" 3416.0 \n",
" 1.0 \n",
" \n",
" \n",
" sigma \n",
" 1.052 \n",
" 0.034 \n",
" 0.994 \n",
" 1.118 \n",
" 0.000 \n",
" 0.000 \n",
" 5247.0 \n",
" 2873.0 \n",
" 1.0 \n",
" \n",
" \n",
"
\n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"simulate_child_parent_income(beta_P=0.75, title=\"Measurement error makes\\nus under-estimate effect\")"
]
},
{
"cell_type": "markdown",
"id": "5d7ce87a-2123-4df0-992d-b28226c8981c",
"metadata": {},
"source": [
"In the scenario above, measurement error can cause us under-estimate the actual effect\n",
"\n",
"- What happens depends on the details of the context\n",
"- **There is no general rule to tell you what measurement error will do to your estimate**"
]
},
{
"cell_type": "markdown",
"id": "11b1a13a-2ffa-4d55-95d7-a093f0f864bf",
"metadata": {},
"source": [
"# Modeling Measurment\n",
"\n",
"## Revisiting the Wafflehouse Divorce Dataset"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "6730d60c-a2d7-4ef0-bb03-17f05d6441f8",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" Location \n",
" Loc \n",
" Population \n",
" MedianAgeMarriage \n",
" Marriage \n",
" Marriage SE \n",
" Divorce \n",
" Divorce SE \n",
" WaffleHouses \n",
" South \n",
" Slaves1860 \n",
" Population1860 \n",
" PropSlaves1860 \n",
" \n",
" \n",
" \n",
" \n",
" 0 \n",
" Alabama \n",
" AL \n",
" 4.78 \n",
" 25.3 \n",
" 20.2 \n",
" 1.27 \n",
" 12.7 \n",
" 0.79 \n",
" 128 \n",
" 1 \n",
" 435080 \n",
" 964201 \n",
" 0.45 \n",
" \n",
" \n",
" 1 \n",
" Alaska \n",
" AK \n",
" 0.71 \n",
" 25.2 \n",
" 26.0 \n",
" 2.93 \n",
" 12.5 \n",
" 2.05 \n",
" 0 \n",
" 0 \n",
" 0 \n",
" 0 \n",
" 0.00 \n",
" \n",
" \n",
" 2 \n",
" Arizona \n",
" AZ \n",
" 6.33 \n",
" 25.8 \n",
" 20.3 \n",
" 0.98 \n",
" 10.8 \n",
" 0.74 \n",
" 18 \n",
" 0 \n",
" 0 \n",
" 0 \n",
" 0.00 \n",
" \n",
" \n",
" 3 \n",
" Arkansas \n",
" AR \n",
" 2.92 \n",
" 24.3 \n",
" 26.4 \n",
" 1.70 \n",
" 13.5 \n",
" 1.22 \n",
" 41 \n",
" 1 \n",
" 111115 \n",
" 435450 \n",
" 0.26 \n",
" \n",
" \n",
" 4 \n",
" California \n",
" CA \n",
" 37.25 \n",
" 26.8 \n",
" 19.1 \n",
" 0.39 \n",
" 8.0 \n",
" 0.24 \n",
" 0 \n",
" 0 \n",
" 0 \n",
" 379994 \n",
" 0.00 \n",
" \n",
" \n",
"
\n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"utils.plot_scatter(WAFFLE_DIVORCE.Marriage, WAFFLE_DIVORCE.Divorce, label=None)\n",
"\n",
"utils.plot_errorbar(\n",
" xs=WAFFLE_DIVORCE.Marriage,\n",
" ys=WAFFLE_DIVORCE.Divorce,\n",
" error_lower=WAFFLE_DIVORCE[\"Divorce SE\"],\n",
" error_upper=WAFFLE_DIVORCE[\"Divorce SE\"],\n",
" error_width=7,\n",
")\n",
"\n",
"utils.plot_x_errorbar(\n",
" xs=WAFFLE_DIVORCE.Marriage,\n",
" ys=WAFFLE_DIVORCE.Divorce,\n",
" error_lower=WAFFLE_DIVORCE[\"Marriage SE\"],\n",
" error_upper=WAFFLE_DIVORCE[\"Marriage SE\"],\n",
" error_width=7,\n",
")\n",
"\n",
"plt.xlabel(\"marriage rate\")\n",
"plt.ylabel(\"divorce rate\")\n",
"plt.title(\"Errors not consistent across states\");"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "ec2c18fc-ea3e-485a-b220-726cb9ae16e4",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAtsAAAHrCAYAAAAe4lGYAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAz3hJREFUeJzs3Xd4HNW5P/DvlO2rXm11F8lyxxUMptm0EEIJSQg1heAk1Jtwk5B+0yC5wT+4oQQnBDDgxKHFCRibYhtswL3Isi3LltWtVV9JW2d25/z+WO9Ko92VVmVXK+n9PA8P1syZmbO7s7PvnDnnPRxjjIEQQgghhBAy6vixrgAhhBBCCCETFQXbhBBCCCGERAkF24QQQgghhEQJBduEEEIIIYRECQXbhBBCCCGERAkF24QQQgghhEQJBduEEEIIIYRECQXbhBBCCCGERAkF24QQQgghhEQJBduETFJ79uxBSUkJ7rjjjrGuChlDDQ0NKCkpweWXXz7kbekcIkP1ox/9CCUlJXjzzTdjdsySkhKUlJTE7HiE9CeOdQUIiYXLL78cjY2Ng5Z79NFHcdNNN8WgRoTEvxdffBE9PT246667kJiYONbVCenEiRP44IMPUFpaitWrV4/qvv/0pz8BAO6///5R3S8ZPW+++SYaGxtx4403Ijc3d6yrQ0hIFGyTSaWwsBCpqalh16elpcWwNoSMPY1Gg6KiImRlZQWtW79+fSCQiedg+6mnnsKNN9446sH2U089BYCC7Xj21ltvYe/evVi2bFnYYLuoqCjGtSJEjYJtMqmsWbOGWq4J6SMrKwtbtmwZ62oQEjV0fpOxRn22CSGEEEIIiRJq2SZkAP5BNSdPnsTWrVuxfv16nDx5Ej09Pfjwww8BAKtWrUJOTg62bduGf/7zn9i4cSPOnDkDh8OBkydPBvZ18OBB/O1vf8PBgwfR3d2N1NRUnH/++VizZg2mT58edOw77rgDe/fuxfr165GYmIhnn30W+/fvR0dHB373u99F3EL//vvv4/nnn0dFRQV0Oh0WLlyI++67b9DtrFYrXnjhBXz44YdoaGgAx3GYOXMmbr75Ztx8883g+dD36mVlZVi/fj3279+PtrY2JCYmIj8/H6tWrcItt9yChIQEAL6+lo888ghuvPFG/PznP8ef//xnbN26FWfPnsXChQvx8ssvq/b5wgsvYP/+/ejs7ERiYiKWLFmCb3/725g9e3ZQHSorK7FlyxZ88sknaGxshNVqRXJyMhYuXIhvfOMbWLRoUci6b9u2Da+++irKy8ths9mQkJCAjIwMLFu2DLfeemvIz2nnzp149dVXUVZWFvhcL7zwQnznO99Bfn5+yLo999xz2L9/P9rb26HX65Gamop58+bh+uuvx8UXXzzoZ/P+++/jvvvuwxVXXBHo6uD385//HBs3boRGo8H+/fuh1+sD6/bs2YM777wTy5YtC7y/DQ0NqnO472fjt2rVKtUx1q9fj+XLl6uWKYqCl19+Ga+99hpqa2uRkJCASy+9FA8//PCAXbdCifQ96jsW46233sJbb70V2Eff19ja2ootW7Zgx44dqK6uRmtrK/R6PUpKSnDzzTfjhhtuUB3/T3/6k+p97T+47sMPP1R1WbBYLPjrX/+KnTt3oqmpCVqtFqWlpbjttttw9dVXD+m1/+hHP8Jbb72FRx99FOeddx6efPJJ7NmzB3a7HdOmTcPtt9+Om2++Oez2O3bswMsvv4zy8nLY7XZkZmbi4osvxpo1azBlypSg8v738MMPP0R1dTXWrVuH48ePgzGGBQsW4N5778WSJUsG3C5U942+16/+50ooLpcLH3zwAT788EMcP34cLS0tAID8/HxcddVV+NrXvgaj0Rgo7z+X/fr+G1CPv+l7He/P4XDgxRdfxLvvvov6+nrwPI/p06fjhhtuwFe+8hWIojpM6vsdeumll0btnCcTGwXbhERg3bp1ePzxx5Geno7CwsKQgy1/8Ytf4B//+AemTJmCadOmoa6uLrBuw4YN+NWvfgXGGNLS0lBSUoK6ujps2rQJW7Zswf/93//h0ksvDXnsffv24bnnnoMoiigqKlL94AzmL3/5C/74xz8CADIyMpCZmYm9e/fi1ltvxXe+852w2506dQrf/OY30dzcDI1Gg4KCAkiShLKyMhw5cgS7du3Ck08+CY7jgo73+OOPgzEGs9mM4uJi9PT0oLy8HIcOHcL8+fODfnhdLhduv/12HD9+HNOmTcOMGTOg1WoD61988UU89thjYIwhOTkZM2fOxNmzZ7F161Zs27YNa9euxZVXXqna5+9+9zt89tlnSExMDLzus2fP4v3338e2bdvw+9//Htddd51qm1deeQW//vWvA+/VrFmzYLPZUFtbi8rKSuTl5QUF27/97W+xfv16AL7+/jNnzkRdXR3efPNNvPfee/jLX/6iCuzLyspwxx13wOVyISEhAdOnT4eiKGhqasLbb78Np9MZUbC9ZMkScByHffv2gTGm+hz2798PAJBlGYcOHcIFF1wQWLdv3z4AwNKlSwfcf1paGhYtWoTy8nJIkoS5c+eqPhP/DVNf//3f/423334bhYWFKCgoQHV1Nd544w2UlZXhzTffVG0/kKG8R3PnzoVGo0FNTQ3S0tJQUFAQ2E9xcXHg36+99hqefPJJ6PV6ZGZmori4GB0dHdi3bx/27duHQ4cO4X/+538C5adMmYJFixbh4MGDABB0c6bT6QL/3rt3L7773e+ip6cHer0eBQUF6Onpwd69e7F371584xvfwA9/+MOIXntftbW1+N3vfge3242ZM2fCarXixIkT+MlPfoKKigr89Kc/Ddrm8ccfx7p16wAA2dnZyMnJwZkzZ/D3v/8dmzdvxvPPP4958+aFPN7mzZuxdu1aJCUlobCwEA0NDfj000+xe/durF27Ftdcc82QX8NQlJeX4/vf/z5EUUR6ejqmTZsGm82G06dPo6KiAh988AE2bNgQuHlMSEjAokWLUFlZCZvNhuLiYpjN5sD+Ihl/09HRgbvuuguVlZXgeR4zZ86Ex+NBWVkZysrK8OGHH+LZZ59Vfd59jdY5TyYBRsgkcNlll7Hi4mL2xhtvDGm74uJiVlxczObMmcM2btzIFEVhjDEmyzKTZZnV19ez4uJiVlpayhYuXMg++OCDwLZOp5Mxxtjx48fZ7NmzWXFxMfvLX/7CvF4vY4wxt9vNfvnLX7Li4mK2ePFi1tzcrDr27bffHtj3z372M+ZwOIL2PZBjx46x0tJSVlJSwl555ZVA3W02G3vooYfYnDlzWHFxMbv99ttV29ntdrZ69WpWXFzMfv3rX7Oenp7AulOnTrFrr72WFRcXs1deeUW13fvvvx+o79/+9jcmSVJgncPhYBs3bmSnT58OLHvjjTcC5a+88krVOpfLxRhj7KOPPmIlJSVs+fLlbOvWrarj/fOf/2SzZ89m5513XtB79+6777KKigrVMkVR2Pvvv88WLlzIFi1apHpdsiyzpUuXstmzZ7P3339ftZ0sy2zbtm1s7969quV///vfWXFxMbv88svZ7t27A8s9Hg975plnWHFxMbv44osDr4UxxtasWcOKi4vZ2rVrmdvtVu2vrKyM/fvf/2aR8n8OJ0+eDCxrb29nxcXFbOXKlay4uJg9+eSTqm3uvPNOVlxczD799NPAMv85fNlllwUdw/+9qa+vD1mH3bt3B74fF110ETty5Ehg3ZkzZ9jFF1/MiouL2YYNGyJ+XUN9j/zn0Q9/+MOw+9y3bx/77LPPmMfjUS0/ceIEu+aaa1hxcTHbs2dP0Hb+7384FouFLVu2jJWUlLA///nPqvoeOHAg8Dls27ZtwNfc1w9/+ENWXFzMZs+eze644w7W3t4eWPfuu+8Gvrfbt29Xbbdt27bAdps2bQos7+npYffee2/gM+5/7fB/xrNnz2aPPvpo4HsryzL7wx/+wIqLi9miRYuCvmODnRv+61ff70bf19f/WtzQ0MA2b97MbDabanlLSwu7//77WXFxMfvTn/4U8XH6Cvc5+vd77bXXstra2sDysrIytmLFClZcXMz+8Ic/qLaJxjlPJj7qs00mlUceeSSQczXUf93d3SG3u+WWW/DlL3850IIoiqLq8aLX68UDDzygetzub4F5/vnn4fF4sGrVKtx9992B7hdarRY///nPMXPmTPT09ODvf/97yGPPnDkTv/zlL2EwGIL2PZAXXngBXq8XV111FW677bZA3U0mEx577LGw2SXeeOMN1NXV4YorrsBPf/pTVWvRjBkz8Mc//hEcx+GFF15Qbff//t//AwDce++9+PrXvw6NRhNYZzAY8OUvfzlkNwyv14u1a9eq1vlbkp544gkwxvDb3/42qPX6S1/6Eu68807Y7Xa89tprqnVXX3110KN/juOwevVq3HXXXbDZbNi+fXtgXWdnJ7q6ulBcXByU0UIURVx22WWq1mBJkvDUU09BEAT86U9/UrXWC4KA73znO7jqqqtgsVjw7rvvBtbV1NQAAL71rW8FtXrNmzcvqLV9IP5H+/7WasDXygr4HuFrtVrVOlmWceTIEWg0Gpx33nkRHycSsizjZz/7GebPnx9YVlRUhLvvvhsA8PHHH0e8r9F8j/yWLFmC888/H4IgqJbPmjULP/vZzwAA//nPf4a83xdeeAFWqxV33XUX1qxZo6rvokWLAq3lL7744pD3zfM8Hn/8cVV3hKuvvhq33347AOCvf/2rqry/RfvWW2/FF77whcBys9mM//3f/0VKSgoaGxvx9ttvhzzejBkz8KMf/SjwvRVFEf/93/+NOXPmwGaz4R//+MeQX8NQ5OTk4JprroHJZFItz8jIwB/+8AdoNJphfUbh1NTU4L333gMA/OEPf1B1+Zo3b17gycGGDRtgs9mCth/Nc55MfBRsk0mlsLAQixYtCvtf/x9jv+uvv37QfYcr88knnwBAyIk/OI4LLPeX6+8LX/hC2P7RA/Hv76tf/WrQOp1Ohy9+8Ysht/P/AH3pS18KuX7WrFnIyclBfX09LBYLAN8j79OnT0Oj0eCuu+4aUj1nzpyJOXPmBC1vbGzEsWPHkJaWFtRn2M8/EUvfoNLv7NmzWLduHR588EHceeed+OpXv4qvfvWrgeC3oqIiUDY1NRVarRY1NTWq5eEcPnwYra2tmD17dsg+4+Hq5u8z2zcAH65ly5YF7d/fhWTlypVYsGABDh8+DEmSAABHjx6F0+nE3LlzI7pZG4qkpKSgmyEAgS4LDQ0NEe9rNN+jvmw2G/75z3/ihz/8Ib7xjW/g1ltvxVe/+tVAN6tIPvf+BvuurFy5EhqNBocOHYLH4xnSvq+88kpkZGQELb/11lsB+MaAOBwOAIDdbsfhw4cBhL7O+G92gfDXGf9++/NfP3bt2jWk+g+Hoij44IMP8D//8z+4++67A5/R17/+dXAch5qaGjidzlE51ieffALGGBYvXhzyO3zVVVchOzsbDocj0J2or9E858nER322yaQy3NR/oVpk+0pJSQk5IKa7uxsdHR0D7mPGjBkAelv0hnrsULq7u9He3j7g9uGWV1ZWAgCefPJJ/PnPfw5ZprOzEwDQ3NyM7OxsVFVVAfC9lr4t4ZGYNm3agPVwu90hbxj86/z16Outt97CL37xi8D6ULq6ugL/FgQBd9xxB55//nnceOONWLRoEZYvX44lS5Zg8eLFQX02/XVrbGwMW7eenp6gut1111349NNP8dOf/hR/+9vfcNFFF2Hx4sVYvnw5UlJSwtY1lHAt20lJSSguLsbSpUuxb98+HDlyJPBvoDdIH015eXkhl/v7zdrt9oj3NZrvkd/x48exZs2awKC7UPqeD5Gw2+2BsRv+1vFw3G43rFYr0tPTI95/uO9FXl4etFotJElCXV0dZs2ahbq6OiiKAq1WG/azGO51xr883Hajpbu7G/fccw8OHTo0aLm+T/mGy/96/O9LfzzPY9q0abBYLKipqQkaSzGa5zyZ+CjYJiQCgw1KDLfe3/IEhB+w4/8BDndxHs4PS9/jhhsVH64+/kemx44dG/Q4LpdLtU2ogXODCffe+YNVm80WsmWpr75BdV1dHX72s59BlmV84xvfwBe+8AXk5eXBZDKB4zi89tpr+OlPfxrU0vjwww8jKysLGzZswP79+wOtxGazGbfeeivuv//+QDcBf906OjoCN1OR1O3SSy/FunXr8Oyzz+LIkSM4c+YM1q9fD1EUsXr1avz4xz8OOblMKJmZmSgsLERNTQ3OnDmDtLQ0VFZW4rLLLgPP84FuL3v37sXSpUsDXUwGGxw5HOE+w/4DaCMxmu8R4Oum9NBDD6GlpQWXXHIJvvWtb2HGjBlITEyEIAiora3FlVdeOeSW575dCwY7P4He70qkwn0/OY5DamoqLBZL4Jrh/39qamrY93yw60y468Rg242WRx99FIcOHUJRURG+973vYcGCBUhJSQl85y6++GI0NzdDluVROZ7/GhnJJGehXvtonvNk4qNgm5Ao6ntBbm9vR2ZmZlCZtrY2AAjqqzhax+3o6Aj5ONrf8h1q2+7ubrz33nuq7A4D8dfdH4SOBv9rWLRoUdj+7KG8++67kGUZ1157bcgsEE1NTSG343ked911F+666y40NDRg//79+Pjjj/Hee+9h3bp1sNvt+PnPf66q23XXXRfohhCpSy65BJdccgmsViv279+P3bt34+2338aWLVtQV1eHf/7zn6r+7gNZsmQJampqsG/fPqSlpYExFmi5Pu+886DRaLBv3z54vV4cOnQIoiiGTXsYT0bzPSorK0NtbS1ycnLw1FNPBfUDD3c+DKbvd6y8vDzi+kQq3E0cYyywzv+98/+/o6MjKDuN32DXmY6OjpCt6f7rRLjtGGMhl/e94R+Mx+MJTDzzzDPPBNXD4/EE6j9a/J/fQDfLg712QiJFfbYJiaLExMRAy4m/q0V/p0+fBuDrTz6ax/W3ypw5cyZkmXDL/Y+NT506FfHxZs6cCcD3WkINJhoO/+PdM2fOQFGUiLfzP9oPNwgwkr65ubm5uOGGG7B27Vo8++yzAHwDR/318NdtKO9Rf8nJyVi9ejV++tOf4u2330ZCQgKOHz+O8vLyiPfhb6Xet29foCXev8xgMGDu3Lk4fPgwDh8+DLvdjtLS0nEVOETyHg3Wkug/H+bMmRMyFdtw+moDvqc4/ptn/3d4NIW7XjQ0NECSJPA8HxjUl5+fD57nIUkS6uvrQ2432HUm3PXAX4/+2/mDVX+Xsv7C1SOUjo4OOBwOJCcnhwz4T506Ba/XG/H+IuF/PeE+O0VRAu/JaF6byeREwTYhUXbRRRcBgGqSFj/GGF555RVVudGyYsUKAAiZRUCSJLzxxhsht/MP+lm/fn3YVqv+8vPzUVxcDFmWQ77O4SgsLERxcTGsViv+9a9/Rbydv391qJawqqoqVRaSSCxcuBCArxuAv1/vkiVLkJKSgoqKCuzZs2dI+wslPT09MDHIQP2K++s7SHLv3r0wm80oLS0NrF+yZAmcTieef/55AEPvQuIfSDlQ3/dYCfce+T/vcN00/K8h1Pkgy3IgT/pA24bbt/+78tJLLw1W/SF77733QtZ5w4YNAHxPfPwBr8lkCtxchvr+uVyuQMaecNcZ/37DLb/wwgtVy/2fxdGjR4O22bp165D6wPvfZ5vNFvK97p95JdS2Q+2mc9FFF4HjOBw4cADHjx8PWv/ee+/BYrHAaDSOi6dBJL5RsE1IlH3jG9+AKIr48MMP8be//S3QOipJEn7729+isrISCQkJYQfaDdfXvvY18DyPd999F3//+98DgbPD4cAjjzwS9sfwK1/5CvLy8rBnzx48/PDDQcGf3W7H5s2b8eijj6qWP/jggwCAp556CuvXr1f1rXQ6nXjttdfCttaF8/DDD4PjOPzqV7/Ca6+9FtSvtr6+Hs8++2wgKwQALF68GADw97//HSdOnAgsr66uxkMPPRTycf/p06fx85//HGVlZaobDEmSAi3bOTk5gQF6Op0ODzzwQOB1v//++0E3JpWVlfjf//1fHDhwILDsv/7rv7Bjx45AhhC/LVu2oLKyEhzHhc1uEsrUqVORk5MDi8WC48ePY/HixaqMOv7g2j8z5FAHR/oHgfn7e8fCUN8jfx392Vb6W7BgAURRxMGDB1U3bT09PXj44YcH7J4w2Ov/1re+heTk5MCMj/1Th1qtVrz++ut45plnBn7RISiKgocffljVzeH9998PBNPf/OY3g+oC+M77vinybDYbfvCDH6CjowM5OTm49tprQx7v1KlT+MMf/hD43no8HqxduxbHjh2DyWTCLbfcoirvHzD417/+VTV4sqysDL/5zW+G1K0mMTExMKHM7373u8Bn7/V6sW7dOmzevDns/vyfUaiMRAMpKCgI3Cz98Ic/VLXEHzt2DL/5zW8AALfddtuQB30T0h/12SaTynPPPReUk7mva665Jmja35EqLS3FT37yE/zqV7/C73//ezz//POYMmUKamtr0d3dDa1Wiz/+8Y8h+1WPxNy5c/HQQw9h7dq1+OUvf4lnnnkGmZmZqK6uhizL+O53v4snn3wyaDuTyYTnnnsO99xzD95++21s3rwZRUVFMJvN6OrqQn19PbxeLxYsWKDabvXq1fj+97+PtWvX4re//S3+7//+DwUFBbDZbGhsbAy0Ig4lu8oll1yCn/3sZ/jtb3+Ln/70p3j00UdRWFgIjuNgsVgCgdIvf/lLVT0WLlyIw4cP44tf/CIKCwshCAJOnTqF9PR0fOc738ETTzyhOo4sy9i4cSM2btyIxMRE5OXlgTGG+vp69PT0QKPR4Be/+IVqm1tvvRVNTU1Yt24d7rvvPiQnJyMvLw+KogSmiAegysG9c+dObN68GVqtFoWFhdDpdLBYLGhtbQUAfPe73w2b5SCcJUuWoLGxEYyxoJZrf/Dt9XrB83zgRiRS11xzDXbs2IFf/vKX2LBhA5KTkwEAP/7xj1Ut6KNpqO/RnDlzAgNFL730UhQVFUGj0WDWrFn4yU9+goyMDNx5553429/+hh/+8Id48sknkZKSgqqqKni9XvzkJz9RnT/9X/+pU6fw7W9/GyUlJYGga+3atcjIyEB2djaeeeYZ3HvvvXjxxRfx6quvoqioCAaDAR0dHWhoaABjDJ/73OeG/D5885vfxIYNG3DppZdixowZsFqtgS4xt956ayC1pN9ll12Ge+65B+vWrcPDDz+Mxx9/HGlpaThz5gwcDgeSkpLwxBNPhE37+OCDD2Lt2rV44403kJeXh4aGBnR2doLnefz6178OGpT6xS9+ERs2bMCpU6dw7bXXYtq0aZAkCTU1Nbj22mvR2to6pJu0733ve/jud7+LjRs3YuvWrcjNzUVjYyM6Ozvx3e9+F5s2bQo5c+/nPvc5vPrqq/jLX/6C999/HxkZGeA4Dt/61rcGnY31l7/8Jaqrq1FZWYmrrroqEPD7u5asWLEC999/f8SvgZBwKNgmk0pNTc2AKazmzp0blePeeuutKCkpwd/+9jccPHgQFRUVSElJwaWXXoo1a9aETT81UmvWrEFhYSGef/55nDx5Em63G0uWLMH9998/4ACm6dOnY9OmTdiwYQM++OADVFVVob6+HhkZGVi6dCkuueSSkDlm77nnHixbtgwvvvgi9u/fj5MnTyIpKQlz587FFVdcETKf9mBuu+02LF26FOvXr8fu3btx+vRpaLVaZGdn4/zzz8cVV1yBSy65JFBeFEU8//zzeOKJJ7B161bU1dUhLS0NN998Mx544IGQ+YILCgrwm9/8Brt27UJFRQWqq6sB+FqOr732Wnzzm99UTXrh9/3vfx+XXXYZXn31VRw4cAAVFRUwmUzIysrCqlWrcNVVV6mmS3/sscfw8ccf49ChQ2hpaYHD4UB2djauuOIK3HXXXcPKFLJ06VJs2rQJQHDLtb9bSXl5OUpKSsJOZBTODTfcgO7ubrz++uuBaesBhJ38aTQM9T3ieR7PPfcc1q5di/3796OsrCyof+8PfvADZGdn4x//+Afq6+vhdDpxwQUX4Dvf+c6A03rfc889UBQF77zzDk6fPh1oce3brWbx4sXYvHkzXnrpJezYsSOQhi8rKwsrV67EZZddFvK7MpjCwkK89tpreOKJJ7B3717YbDbMmjULt912W9i83t///vexaNEivPzyyygvL0dbWxsyMzPxhS98Ad/+9rcDOcxD+dznPofS0lI899xzOH78OBhjOP/883HfffeFPC91Oh1eeuklrF27Ftu3b0dNTQ1yc3Pxwx/+EF/72teGnG//8ssvx1/+8hc8/fTTOHHiBKqrqzFz5kz8+Mc/xhe+8IXAOd7fkiVL8Pjjj+Oll17C6dOnA9f3G2+8cdBjpqamYuPGjXjhhRewZcsW1NTUgOd5zJs3DzfccAO+8pWvjPrAVzI5cSzSTpmEEEIIiaof/ehHgW4pw5kTYKguv/xyNDY24sMPPwz0wyaEjC7qs00IIYQQQkiUULBNCCGEEEJIlFCwTQghhBBCSJRQsE0IIYQQQkiU0ABJQgghhBBCooRatgkhhBBCCIkSCrYJIYQQQgiJEgq2ybjGGIPVag2aKpuQgdB5Q4aDzhsyXHTuTG4UbJNxjy5eZDjovCHDQecNGS46dyYvCrYJIYQQQgiJEgq2CSGEEEIIiRIKtgkhhBBCCIkSCrYJIYQQQgiJEgq2CSGEEEIIiRIKtgkhhBBCCIkSCrYJIYQQQgiJEgq2CSGEEEIIiRIKtgkhhBBCCIkSCrYJIYQQQgiJEgq2CSGEEEIIiRJxrCtACCHxRnG6IdVY4O2yAQCEJDO0hdngDboxrhkhhJDxhoJtQgjpQ7Z0wHm0CvAqgWVeqw1SQwsM86ZDk506hrUjhBAy3lA3EkIIOUdxuoMC7QCvAufRKihOd+wrRgghZNyiYJsQQs6RaiyhA20/r+IrQwghhESIgm1CCDnH30d7pGUIIYQQPwq2CSHkHCbJo1KGEEII8aMBkpMQZVoghBBCCIkNCrYnGcq0QAghhBASOxRsTyKRZFoQkkzUwj0K6OkBIYQQQgAKtieVSDMt6EsLYlepCYieHhBCCCHEjwZITiKUaSH6KE8zIYQQQvqiYHsSoUwL0Ud5mgkhhBDSFwXbhIwienpACCGEkL4o2CZkFNHTA0IIIYT0RcE2IYQQQgghUULBNiGEEEIIIVFCwTYhhBBCCCFRQsE2IYQQQgghUULBNiGEEEIIIVFCwTYhhBBCCCFRQsE2IYQQQgghUULBNiGEEEIIIVFCwTYhhBBCCCFRQsE2IYQQQgghURKXwfamTZvw85//HDfddBPmzp2LkpISvPnmmxFtW19fj/POOw8lJSX4+c9/HuWaEkIIIYQQEp441hUI5cknn0RjYyNSUlKQmZmJxsbGiLZjjOHHP/5xlGtHCCGEEEJIZOKyZfs3v/kNtm3bht27d+OWW26JeLuXX34ZBw8exIMPPhjF2hFCCCGEEBKZuGzZXrFixZC3qa2txdq1a3H33XejtLQ0CrUihBBCCCFkaOKyZXuoFEXBI488gqlTp+Lee+8d6+oQQgghhBACIE5btofqpZdewqFDh7BhwwZotdohbcsYi1Kt4g9jg7/eSMrEE39d46XOE/E9nojCnTf0+ZGBxNv1howfdO5MXBzHDVpm3Afb1dXVeOKJJ3DnnXfivPPOG/L2XV1dk+bkd9lsYE73gGU4xQOP1RqbCo2irq6usa4CgIn9Hk9E/c8b+vxIJOLlekPGHzp3Jp6UlJRBy4zrYNvffSQzMxMPPfTQsPaRlJQ0upWKYzazGQo/8EfOG/UwJyfHpkKjgDGGrq4uJCUlRXR3GW0T8T2eiMKdN/T5kYHE2/WGjB907kxu4zrYXr9+PQ4fPoyXXnoJBoNhWPuYTCc9xw3+eiMpE484jouLek/k93gi6n/exOvnpzjdkGos8HbZAABCkhnawmzwBl1M60F84uV6Q8YfOncmp3EdbFdUVIAxhjvvvDPk+o0bN2Ljxo1YtWoVnnnmmRjXjhBCRk62dMB5tArwKoFlXqsNUkMLDPOmQ5OdOoa1I4QQMphxHWwvXboUgiAELW9tbcVHH32EadOmYdGiRZg9e/YY1I4QQkZGcbqDAu0ArwLn0SoISSZq4SaEkDg2roPtL37xi/jiF78YtHzPnj346KOPsHTpUvzqV78ag5oRQsjISTWW0IG2n1eBVGOBvrQgdpUihBAyJHEZbL/22ms4cOAAAKCysjKwbO/evQCA1atXY/Xq1WNWP0IIiQV/H+2RliGEEDJ24jLYPnDgAN566y3VsoMHD+LgwYMAgJycHAq2CSETHpPkUSlDCCFk7MRlsP3YY4/hscceG/b2y5cvx8mTJ0exRoQQQgghhAxdXAbbhBAyVpjHC3d1E6SqBnh7nGBeBbxeC82UNOhLC8AbaTAiIYSQyFGwTQgh8E064Tx6BrbPysFckmqdYnPA02aFs7wK+hl5MCycCSHRNEY1JYQQMp7wY10BQggZa4wxWN/8CLadZUGBtrog4K5uQssT/4SnjaZdJoQQMjgKtgkhk17PB/th33U04vLeThvantsEZaDAnBBCCAF1IyGETHLeHgd63tsbtFwzNR266TngdBrITe1wV9aDyZ7Aek9bF2wfH0bilctiWd0xRdPGE0LI0FGwTQiZ1Oy7j4F51BPHJFy2CLrpOYG/dYVTYFwwE93v7oa3x9G77aflSFi9BBw/8R8S0rTxhBAyPBP/F4IQQgbg2HtC9bdueq4q0PbjjTqYLl6gWubtssNdURfV+sWDSKaNV5zu2FeMEELGAWrZJiNCj5XJeMY83qCBjvpZeWHLazKSIaYmqpbJlg7oZxdGo3pxg6aNJ4SQ4aNgmwwbPVYm412o2Rf5QVL6CUkm1cBIJk/8GRxp2nhCCBk+CrbJsETyWFlIMlELN4lrnF4LcABY7zJPWxcEkyFkecYY5DYrBLOxzz7G5zk+lKdSNG08IYQMH/XZJsMS6WNlQuIZx/PQ5maqlrmO1YAxFrK8VNcMpcepWqYtyIpa/aJFtnTAtqsMUq0FXqvN90Sq1gLbrjLIlo6xrh4hhEwoFGyTYaHHymSiMF0wR/W3fLYVjr0nwPrdTMrNvgC1L012KrQF2VGv42iiwY6EEBJb1I2EDAs9ViYThWFRCbr+86kqwHQerYLrdAN0RVPA67SQm9ohW9oBAJzYe9k0rZwPjuNiXueRoMGOhBASW9SyTQiZ1HidBslfvixoOXO64TpeA8ehykCg3ZduRg5M588JWh7v6KkUIYTEFgXbhJBJz7hwJpK/dKlvsGQEdNOmIu0b147LyWzoqRQhhMTW+PulIISQKDCvmIekL1wEMSt8ykpOr4Nx6Sykf/t6yrRDCCEkItRnmxBCztHmpCP5ugvh6eiGu6oR3m4HoCjg9VpopqZDWzgFQoIBnIYunYQQQiJDvxiERAFjDJ5WK5QeB5hXAafXQpOZAl6vHeuqkQiIqYlBM0USQgghw0HBNiGjSJG9vkwWJ2rh7bar1nECD23RVBiXlIxR7QghhBASaxRsEzJK5OYOWDd+CE9Hd8j1zKvAfboBUo0FQqIJCVcsHXdp4wghhBAyNDRAkpBR4GnrQutTb/r6+Eag+9096Nm6N8q1IoQQQshYo2CbkBFijKFj/RYoNmfQOt6o9/X9DdGC3b11L1yV9bGoIiGEEELGCHUjIWSE3KcbIdW3qJaJaYkwnT8XYnYqOI7zTZFdfgbOstOqcrYdh6AvzotldQkhhBASQ9SyTcgI2T8pU/3NmwxIvHYFNFPSAn2yeYMOpqWlMC6epSrrqqiFp70rZnUlhBBCSGxRsE3ICDDG4KqoUy0zzJsGXqsJWV4/twicts8DJYag7QnpT3HLkOqa4TrVAPfpBsjNHWCKMtbVIoQQEgHqRkLICDDJA+ZWT22tyc0IW54XRWimpoG5PYFlSk9kgyrJ5ONp74ZjfwXcpxvBvF7VOt5sgH5WAYyLKZUkIYTEMwq2CRkBjg8e+Mi8bMBtmKdfi2SIfRBi23kEnf/cBiZ7Qq5XbE5fIF5ZD/2sfGjzsmJcQ0IIIZGgbiSEjACnEcGb9Kplck1T2PKK0w1PU5tqmZBkjkrdyPjV89FhWN/8GBj4vg0AoDjcaH36Lchn2wYvTAghJOaoZZuQETIsmAH7p+WBv53lZ6CbmQsh0aQqxxiDY1+FquWbE3no5xTFrK4k/kn1zejatDN4BcdBTE6AInug2NRdj5hbRvsLm5H1o9vACUJU6sUUBXJDK1yn6uHt7AHzKOC1IsTMFOhK8qEroJZ1QggJhYJtQkbItGKeKthmsgfWTTthPK8Yuhm54HQaeCztcB6pgtTQAk7s/doZFs6EYDaMRbVJnLLtOBzUoq2fXQTjwhngjb6nKHKrFfbd5fA0dwbKeNq64DpWA8P86aNeJ3eNBbbtB+HtN75AAeDp6IarohaanAxop+dCm5M+6scnhJDxjLqREDJC2px0GOZNUy1jbhn23cfQ8cpWtD//Nrre+QxSgzoXNyfySLh8cSyrSuKct8cB5xF1LnbDvOkwr5gbCLQBQJORjKRrLoCQrO6CZNulTkM5GmyfHEX3u7uDAu3+PM2daP3T63Cfbhj1OhBCyHhGwTYhoyDlq6uhmZoW+QYckPLVK6CZMoRtyITnPtUA5u0dQMsJHAwLZoQsy4kCDAtmqrc/3RB2QGV/it0Jd00TXCfr4K5qhNzeDcbUTerOY9WwvrEjor7jgO8ms+2vb0Nu6Ry8MCGETBLUjYSQUcAbdMi49yY0/eJ5uKvDD5AEAE4rIu3rn4Nh3ug/7ifjW//WYyEjGbxeG7a8NjcD9r4LGOC1OyEmJ4TdRmpohePASUh1lqAgWkxLhL60EIYFM8AUBV1vBQ/SFFMSoSvJA59ogrfVCtfJOigOV28V3DK63/kUaV+/drCXSwghkwIF24SMEt6oR+LnL4R0phHO4zWQaixAn5ZCIdEEfWkBDAtnUqBNQgpKJTlYGklv8MQ2HB/6gSVTFFjf+Ahd//4EzBO69dvT3g3brjK4qxrBmfTwtHer1hsWzIBxyazAzKjIz4JhwQz0bD8IubE3G4qzvBreLhtl2iGEEFCwTcio4jhAMyUdminpvqwRdieYl0HQa8AZ9eA4Drwu9OyShAgp6hZpT5sV3m57UGYbP6lfmklO5INSUQK+TDid//gQjn0VEdXD09qF1sf/ATEjBZzG9zMhpiaqA+3AMQWYVy7wdTfxU3yZdxJWL4noeIQQMpFRsD3JeG1OuCpqITe2QnFK4HgOvNkA3fQcaKdNBQ/dWFdxwuA1IvgBHucT0p9uZh44naZ3VlIG2PceR8KqJUFBruJwwdFvMKV+7vSQqf8ce0+EDLQ5QYCQbILX7gZzuVXr5OZOeHuc0Jfk++pWnBdUBz9er4W2aKpqNlTZ0jH4CyaEkEmAgu1JQrG70P3ePrhO1gb1wfR22yGfbQO3+zhMy0thWrkg7I8qISR6eJ0GpmWlsO3szSoi1VjQ/c5nMCyYAU1OOpjshftMI5yHToE53apUkuaL5gXtkzGGnu0HVcs4QYBxcQl0s/LBazVgiuLry/1ZeW+/cYXB29EN5VwQzodpXfcTkkyqYJvJ8pBfPyGETEQUbE8C3i4bWp9+C+7TjQNmFWCyDPvu4xDT30fKravD9v0khESP+dLzYN97ord1G4BsaYdsaR9wO//Tqf6kqkZVPm4AMK+cD92M3MDfHM9Dl58FMS0R1td3AAzgBB4MvpR+nFYDb5sVyA8/cY2ntRPM44W3xwEmeyA0tsF1ohbawmzwBnpiRgiZvCiamuCY7EHbun/D02qNeBvHgZPofuez6FWKEBKWmJqItDuvBvoPlhxom/QkpN51dcgnUq4Ttaq/hWQTtNNzQu5HMBmgK84DAPAJRgCAt7PHt5+T9WAeb8jtvN12uCrrIVs6oNgcYG4JvEkPqdYC264y6lJCCJnUKNie4Oy7j0E+q24R4wQB+tICJKxaDPPFC6HNzQzarmf7QXjau2JVTUJIH/rZhUhfcz24AdL++YnZKci4/2YI54Lj/vqnE9TmZg3YTcx/PdBkpQBAIG+3YneiZ/tBKG5Jvf9uO7q37oHS5QCYLzsKpxGhm3YuoPcqcB6tguJU9wknhJDJgrqRTGCMMVXfTwAQEo1IvHYFBFPvFOH64jxIdc3ofn9fn40B+6flSLruwlhVlxDSh744D6m3Xwln2Wm4T9TCc66FGQDAAdq8LOhLC6ErzoOQGDrQBoJTATIlOF2gar3X13rNmw3gzQZVdxap1oLOf7RBW5ANIdEIT3s3pLpmX0DfJ82lbkaOOuuOV4FUY4G+tCCSl04IIRMKBdsTmFRrCeo+YloxXxVo+2nzfT/c7sr6wDL73hNI/PwKGixJyBjhtSIMs4tgmF0ExeGC159ByKQHr/UFs4N9PftP6e6uaYJxcQk4IfSDTanGcu5fHHTTpsJda1GtZ7IneEr2Pt1LhAQjjItKgvbr7bINXFFCCJmgqBvJBOZpUncf4RMM0OSkhy2vO5fiy0+xOaH0OKNSN0LI0PBGPTRpiRBTEgKBdiQM89UTKDGHG86y0yHLyq1WuKt6A2neZEDqbVeC0wSnE1TtU/G1agsJRiRec37IAZFMouwkhJDJiVq2JzBFUs8SJySYBmylDvUomtJ3ETK+aaam+1qoz5wNLHMcOAmv1Q793CKI6UlgDhdcJ+vgLKtSZSziRB7JN14MIckE28dHQk7xDvgyl2hyMpD0hYso8wghhPRDwfYExvcbXOVp7wJTlLAp/bxt1hD7oB9OQsa7hCuXwv3nTapl7qoGVSt2KKYV8yAkGCGmJCDxiqXw2p1wn2qAt7MHzOMFpxEhZiaDuT3gNAIF2oQQEgIF2xOYtiBb9Tdzy3BXnYV+Zm7I8s7jNaq/xfREcEb68SRkvNOX5CPpugvR9Z9PIt5GV5wXNEBaMBlgXDgzqKzzWLVqgCQhhJBe1Gd7AtNkp0JbpA647Z+VQ2psVS1jigLHoUpI1U2q5aYL5tHgSEImiITLFyH5S5eCEwb/TuuK85B+9+fBiQP31SaEEDI4atme4MwrF6CjujebAJNkdL+7G2JmCjRT08EkD6RaCxS7eiAkpxVhXFYa6+oSQqLIvGIevF12OA9VwlVRB8XhCqzjBAHa6VN96QQLssBp6OeBEEJGA11NJzjDwpnQH6yEq7xatdzT0glPS2eYreAbFGUOThFIJj7F6YZUYwmkahOSzDTl9gQimPQwLiqB4bxiX8YhtwxOFCCYDdSSTQghUUDBdgyNRRDDcRxS77gKHS9shtwvFWA4SdddCNP5c6JWJxK/ZEsHnEerAG/vxCdeqw1SQwsM86ZDk506hrUjo4njOAgJRggJY10TQgiZ2CjYjpGxDGJ4rQZpd38eilOC4+BJ1aPjvsSMZJhXzEPC5YuiVhcSvxSnO+gcDTg35baQZKIWbkIIIWQIKNiOgXgIYjhBgHHRTOiKcyHVWiA3tkFxuQGOh2A2QDttKjQZyeAp+8ikJdVYQp+jfpNgym1OqwEc7sHLkLhF3aAIIfGGgu0YiKcghuN56IqmQlc0NerHIuNLJNNpT/Qpt4UkM7xW9WtUJBmelk4odt8TId2MXChONwVvcYi6QRFC4hEF2zFAQQwZDyKZTnuiT7mtLcyG1NASCNY8nT2QavvcLAs8mMcL264yCt7iTDw8QRwL1JJPSPyjPNsxQEEMIeMDb9DBMG86IPBQJDko0NYWZIPXaQLBm+IcuMsJiZ1InyBOJLKlA7ZdZZBqLfBabb5W/FoLbLvKIFs6xrp6hJBz4rJle9OmTThw4ADKy8tRWVkJWZbx6KOP4qabblKVk2UZ27Ztw/bt23HkyBE0NTWB4zjMmDEDN9xwA2655RYIAqWyIoRETpOdCiHJhJ5tBwOtg7xRDzEzxRdo+02CPuzjyWR7gjhZW/IJGY/iMth+8skn0djYiJSUFGRmZqKxsTFkubq6OjzwwAMwmUw4//zzcfnll6Onpwfbt2/Hr371K+zcuRPPPvsszYJICBkS3qCDmJYITsgfsNxECt7Gu8n2BDGexgIRQgYWl8H2b37zGxQUFCAnJwfr1q3D448/HrKc2WzGL37xC9x4440wGHonYHE4HLjjjjuwfft2bNmyBddcc02sqk4ImSAmW/BGxpfJ1pJPyHgWl322V6xYgZycnEHLZWVl4dZbb1UF2gBgNBrx9a9/HQCwb9++qNSREEIIGSt0M0jI+BGXwfZoEEVfoz312SaEEEIIIWMlLruRjIY33ngDAHDRRRcNWI4xFvW6MDb4cSIpE0/1iJ/XxFT/H2vx8r4Mx3iu+1BFct7Ew/sRjToMVl6RZMitVmD3MQCAkGSC1yX7Ngy7TwYOI38/RyIWn1c8XW/i4fwkkYunc4eMrkjGBU7IYHvjxo34+OOPcf755+OSSy4ZsGxXV1fUT36XzQY2SIowTvHAY7WOm3rEy2vy6+rqislxBhNv78tQjOe6D9dA5008vB/RqMNA+1SsNnjqWwFRgNtw7uehAZAqGyBmJoNPNofcTnK5AAZ47fZRq+dQxfLziofrTTycn2To4uHcIaMrJSVl0DITLtjesWMHfv3rXyMnJwf/+7//O2j5pKSkqNfJZjZD4Qd+q3mjHubk5HFTj3h5TYwxdHV1ISkpKS6yzsTL+zIc47nuQxXJeRMP70c06hBun4okw93aDUGrAa/TQm8yBdYJWi1Yazd0GangQ0xXz+v14ADVNiOt51DF4vOKp+tNPJyfJHLxdO6Q2JtQwfbOnTtx//33Iy0tDS+99BIyMzMH3SYWJz3HDX6cSMrEUz3i5TX1HouLiwtYvL0vQzGe6z5cA5038fB+RKMO4cp7W62AwsCBO1euXxmFwdtqhZAbfF3lOA5gbMB6jMf3Kvx+xv56Ew/nJxm6eDh3SOxNmAGSH3/8Me69916kpKRg/fr1yMvLG+sqEULIuKHYXWHXcRox4jKEEELUJkSw7Q+0ExMTsX79ehQUUBJ/QggZCubxhl3HG/URlyGEEKI27psi/IF2UlIS1q9fj8LCwrGuEpnkFEmGp6Uz0ArIm85N9R2iryshscZpNYBjkIF1/VqpxcwUeNoHGNjFcxAzU6DYnZBbrWBuGZwogE8y+2bipMfmhJBJLC6D7ddeew0HDhwAAFRWVgaW7d27FwCwevVqrF69GlVVVbj33nshSRKWLVuGd955J2hfOTk5uOmmm2JXeTKpeTq64aqoU02jrNic8LR1QVuQDTElYQxrRwaiON2QaiyBWfekumbwCSbwuol1kyQkmeG1DjyzYP9Wal6ngbYgG7KlPcQOefAGHWwfH4ZUZwH6JXcSUxOhn10Iw4IZI606IYSMS3EZbB84cABvvfWWatnBgwdx8OBBAL4AevXq1Whra4MkSQAQMtAGgGXLllGwTWJCcbohVTepAu0ArwKp1gLepAdv1MW+cmRAsqUDzqNVqs9Obu4Ea2gd05ukiFqhh/jERFuYDamhJfR5CgACDzEzOJWVmJIAMT0J2tzMwA0Jn2iE4/Ap2HYcBvN4Qu7O09EN264yuE83wLioGEJS6PSBhBAyUcVlsP3YY4/hscceG7Tc8uXLcfLkyRjUiJDBSTUWsHABDAB4FXhaOiGGyVVMxobidAcF2gH+mySjfkxauCNphR5q8MobdDDMmx76NQs8tAXZYV8rr9NAX+obE8MYg3XjNjj3RXYN9rR1o/WpN5HxwM0QEoxDqvNwKe5zXbocvi5dYloiFKcbvIFueAkhsROXwTYh45G/tW8gA2VzIGNDqrGEb+UFAjdJ2rzBU4mOtkhaobWF2UPeryY7FUKSSdVtRpOVMqRuM459FbDvOR60nBMECMkmeO1uMJe6Vd7T1oXOf25D+jc/P+Q6D5WnswdSrfqzlT1e2HaVwTBvOjTZqVGvQ7QpLgnuU/Vw11jAzt1QcEYddHlZ0JXk01M0QuIEBduEjBImyYOXGSCbAxkbEd0kOcbmJmmwVmjDvOnDbqXlDbpAKzXgO3+VQbqsBMoyBtv2g6plnMDDuHgWdLPywWs1YIxBrm+B/bNyeHscgXKu8mp4Wq0QM5KHVe9IKG45KNAO8CpwHq2CkGQaty3czOOF7ZOjcB4+Debtd03pccDT3AnHoUro502D8fw5NDibkDE2IVL/EULIcEV0kySH7o8cC5rsVJgvmg9tQTaEZDOEZDO0BdkwXzR/zFpnpTNnIVs6VMvMFy2AYf70QGDHcRy0+VlI+vyKoOwmtk+ORqVeiluGVN8Cx97jkOqa4WnrghLqBter+J5ojEOKJKNt3b/hPFIVHGj3wbwKXOU1aHv2X1BcUgxrSAjpj1q2CSEkzvVvhR5rrhO1qr+FJBO0M3JCluVNBuhmFcDdZxvXiVrghpWjWidPRw9cFbVQ3BLc1U3nnkZw4PQaaIumQDSr+4lH8kQj3jDG0LnhA7hPNUS8jVRjQcdL7yLtni9QCkZCxggF2+NM//RkQpIZ2sLscfs4lBAy/vTtFgIAmryMAQM5bU6GKthW+m0/UorTDWdZFVwVtfBYOuDpsgGstwuJVNME3axCGEp6ZxaO5IlGvJGqGuE8clq9kAO0hVOhO9d3X6pvgbuqEWC9ORhdFXVwHa+BYU5RLKtLCDmHgu1xJFR6Mq/VBqmhZcIM+CGExD+O7xdYe1nogucwpV/f6f7bj5D1jR3o+egQEKZbBZNkuE/UwNPcAU1uJjRR7C8eTbZd6u43nEZE4ucuUL0e3fQcGOZPR9c7n6o+F/uuMgq2CRkj1Gd7nBgsPZnzaBUUZ2SDmwghZCSEfukr3TVNA6a9lGqaBtx+JHo+Ooyud/cETaYTxMvA3DK6N38Gub171I4fK4rD5fsN6MO4qDjkjYOYmgjT0lLVMldF3bjsOkPIREDB9jgRSXqy8TrghxAyvhgWzFT9zZwSnIdPhSwrt1rhPlXfb/vRmU1SamhB16adwddGDuB1OnCaPlk4znUrYbIHtg/3B7e2xzlPqxVQ1HcUuuK80IUBaGfkgBPUTxD6D2olhMQGdSMZJyJpkaBWi7HHaUQw98B9QYc64x8h8UaTnQrd9Bxf3+BzHIcq4emywTBvOsT0JDCHC67KejjLqlStzpzAw3T+nFGph237oaAWbU1eFnR9sp94HS7Ija1gfdIaervtcFdbkHDpqFQjJhRJnRGH02nA67Rhy/OiCM6oVy0by6w6hExmFGyPExGlJxuHA34mGt6oh2JzDliGpqsmE0HClUvhfrZRtUw6cxbSmbMDbme6YO6ozCDp7XHAWaYeLKgtmAJ9cR6kuuZAYCkY9eCn5wZl8HCVnxlxHWKJN6gDa+aW4bU7IZgMIcsrTjcUmxNCnywsvD58cE4IiR7qRkLIKBIzUwBhgK8Vzw1rxr/xRHG64TpRC/vuY7DvPgbXiVoaTzAB6YvzkHT9RUPaRjczd8jbhOM+3QDm6e0KwgkctEVTQpbleA7iFPUAcvlsG5Rx1EChyUoF1292T2d5ddjyrhM16icKIg/N1PQo1Y4QMhAKtgkZRbxOA21BduiAW+ChK5o6odM0ypYO2HaVQaq1wGu1+bLl1Fpg21VG/UUnoIRLz0PKVy4HJw7+U6IrzkP6t64DJwqjcuz+6QP5RBN4TfiHtUKiSb2Ajd3MoMPBaUQYl5SolrnKq+A6WQfWJ80fYwzuqkY4D6tb/Q0LZoLv162EEBIb1I2EkFEmpiSAN+rhaekM/JjzRj3EzBSIKRO3C0kkGXPG8xTZJDTT+XPgsfbAeegUXBV1UOy93ag4QYB2+lToSwuhK8gKmklyRPh+Ab4ySDqSEOvH2yQv5gvnw/5peW+LNQNsO4/AWVYFbX4WwAFyfQs8nT0AAK7Pe2ReOX8MakwIASjYJiQqeJ0G2rzMsa5GTEWaMSeeZkIko0Mw6mE8rxiGhTOh2JxQ3DI4UYBgNoxaS3bQMfulD/R226E43WFv5rxW9QByTuDAh+nvHK80U9KQsHoJet7fr1ru7bLBeTT8AHnzyvm+J26EkDFBwTaZcGiWzbExETPmeLrtkBtbodiccJafAScIEFLM0M3MG7ObqXg+vzmOg5BghJAQ/WPpZuaB12uhuCTfAga4KutCphVUZA9kS7uqG4W2aGrUbgSiKfGa88FcEjr/uT2i8sZlpUi6YWWUa0UIGQgF22RCoVk2x85Eypjj6eiG7dOj8LZ3AQA4XggEZt4uG6QaC3izAQmXngfzxQtjVi86v3vxOg2My0ph+/hIYJmnpROOfSfB6zXgDDowhUHp7IHc3A541E9d9HPH52yKHMch6caL4Wnvhn3PsaAWez8h0QTjstlIuWXVuOsuQ8hEQ8E2mTCozzCJlOJ0Q65sgJ2dhetELSDwEDNTwGs1kJs70L1lD7ydA7fCKzYnurfshbZwChIuPS8mdabzW818yULY9xxXLfNauyH1OAOT2PhxXG//ZU12GjRTxm9mDo7joC/JgyY3Ax5LO9zVFijOc+NDDHpo87OgyUmHYNJToE1IHKBgm0wY1GeYRMLfOuzp6obXZIK3xwHmluFp64KYloie7YfCTv7BvIovjeG5VlJOp4H1jY8gpCTAOEqzIoYTq/Nbcboht3SCuSRwogA+0QQxPSkugzYxNRFpd10Dy69fDLmeKQBkD6AwMIEHL0ngE01IWLUYcfhyhozjOGimpI/rGwdCJgMKtsmEMRH7DJPR5W8dZh5v8Eqvgp7th6A43eD6pG7ktBpoi6ZAPyMX7jNNkBua4XX58oYzRYFs6UDHK1thmFsEToheH+Bon99ySyfse09ArrEETWUuJJmhLy2A4bziYe8/WvSlBUj8/Aq0r98K5vJ9rszjgWJ3gfWZdZHjAA8DNB4F7qpG8EmmcLskhJBRRXm2yYQxkfoMk+gYqHVYkWR4mtpVE/CIaUkwrZwPXdFUMK8CbW4GjOfPgW5Gbu+GTIFU0wzHwcqo1j1a5zdjDF2bP4P19Y8gnTkbFGgDviDevvsYrP/8EHJL55CPEW3a3AwkXLoQ+lkFUNweeK02daCtFX0t9FkpYA7nudeyDZ5W69hVmhAyaYwo2G5tbcWrr76K3/zmN/jJT34SWN7R0YGysjK4XONnwgBCyMQ3UMuvx9IOxhTgXKs3J/DQz50Gnhfg7ewJBOkcfDMVCol9Us8xBd1b9ka17tHS/c5nQankwvF2OdD61JvwdHRHuVZDx2lEeJ0SeJ0IMTURQpIJQrIZQmoihEQjeK169kWv1Y7Wp9+Ex9ozRjUmhEwWw+5G8uqrr+L3v/89JMmXdonjOPz2t78FALS3t+MrX/kK/ud//gdf/vKXR6emJO7IrVa4T9bB09wJRZZ706IV50ObmzHW1SPDwDxeuKubIFU1wNvtAFMYeJ0Gmqnp0M/KB28c34PvBmr59fb4JmNh5yY/8Q+YBADFJalmIOTAQZOTDu/J+sAyubE1GlWOKteJWvR8eCBoOcfzEJLNUFxS0CyLSo8DHS9vRcYDN8dVP265uRNybZPvD54Hdy625jgenE4T+Fz78nbZ0fnye0i/76a4ei2EkIllWMH2tm3b8Otf/xpz587Fvffei48//hj/+Mc/AutnzpyJkpISfPDBBxRsT0BSYxs63/wIckNwcNGbFs2IhMtimxaNDB9jDM6jZ2D79CiYWx2QKjbA094F59Eq6GbmwrBgBoSkiTcTJuvXvUR1YxGi6wlvUE99rbilqNQrmnq2H1Qv4ADDwpkwzJkGXq8FYwweSzvsnx6Dp7O3NVuqsUCqsUBXNCXGNQ5POnO23xIOmuxUCOnJ4EUBnEaAkJoIb1tXYIZFAHCfOQu5vsU3AyMhhETBsLqRPP/885g6dSrWr1+Pyy67DGlpaUFliouLcfr06RFXkMQX9+kGtP7pdXgsA/fbVGwOdL+7B7adRwYsR8YeYwzWNz+CbWdZUKDdn1RtQcsTr03Ivq6cRj240dtlH7C80qNe37+bQryTWzrhPtWgWmZaOhumxbPA67UAerNdJH7+AtWEMADi6rvtae8O6tqizUmHJjsNfJ+Ja8T0ZCReuwJcv/SI8fRaCCETz7CC7RMnTuCSSy6B0WgMWyYrKwvt7e3DrhiJP55WK9qff2fQgKwv65sf+3IDk7hl+/AA7LuORlzea7Whbd0m1UDCiaBvNxHA15rvtTtDlmUKg7u+RbVMzBpfE8q4TtSq/uZ0GuhnF4Ysy+u0MMydplrmrqiLVtWGTKprDvybyV7fkwiOh8dqg9IvjSOv18IwRz2hjSuOXgshZOIZVrDNGIMoDtwDpaOjA1qtdliVIvGpa/Pu3qmRz+GNOhjPK0bClctgOn8OxNTEoO2s/9oZMsMBGXveHge6t+4JWq6Zmg7zygVIWL0E+jlF4DTq77unrRu2jw7HqJaxIWanghP6vE4GOI+chtfer8+yxwtn+Rkwl/pmQz97fOVvV2wO1d/anPQBpy8Xc9TjMBSnO3QKxTGg2J1Q7E54mq2QmzvgtbsgNTRDqjkL17EzcJ9phLfLDsZ8/bY1/V+L3RlYRwgho21YfbaLiopw4EDwoBo/j8eDffv2obg4/nKykuHxdtngLFN3C9LmZyHh8sWqH2j9nCI4Dp6E89Cp3m07euA6URvUmkTGnmPPcbB+01gnXLYIuuk5gb91hVNgXDgT3e/uhre7N0Cz7z6GhCuWRDW3dCzxggAxJ101FkGxO2H/tAycRgPF6YLilOBp7gDzqoNMMTVx/E0s0m9AYP/zIEioG2Z+7AcVenscsO+r8HUjCQTM6sDZ222H0uMEpxFhKC0Mfi00OJIQEkXDatm+7rrrcPz4cTzzzDNB67xeL37/+9+jvr4eN9xww0jrR+KE40Al0Gc0PycKMF+8MKgljOM4GBeVQExTt3A7+k2pTOKDfe8J1d+66bmqQNuPN+hgWrlAtczbZYfr5MR6/K6bPhWCWR+03Gu1wV3VCPlsa1CgzQkCTCsXjLt4TUxJUP0tn20dcJCnVK0egCgkGsHxYztVg+Jwoe2Zt6DY1E8fmOTpP1s7AN/g7p4dB+Hu/1qSTJSNhBASNcNq2b799tuxbds2/OlPf8KmTZsC3UUefPBBlJeXo7GxERdeeCFuvvnmUa0sGTtyU5vqb23RlMAgqv44joOuOB+OfRV9th+d/vuK0w2pxhLIl8wnmaAkG0Zl35MN83qDBjrqSvLCltdkJAfdRHmaOoDZE+eJBS+KMC4thVRjgadj8PzLnCggYfUSaDKSo1+5UaafMw3gdwRuoplHgWNfBUwXzgsKPL1dtqA+3oYoT08fCeubH0O2dIBPVI8f4jgO8HoBQYP+Ubf7ZD0UpxtiWlJgmWHBzJjUlxAyOQ2rWUKj0eD555/HPffcA6vVilOnToExhq1bt6Krqwvf+ta38Oyzz1JLwQTC+g8ySgw/OBYAhAT1+v7bD4ds6YBtVxmkWgu8Vhu8Vl+aQffu45AtHSPe/2TTd4Y9v/6fW9D6flNcj8bnGm94ox5JX7gQ2sLsAbvIiBnJMF0wB9q8zBjWbvQIiUYY5qsDZldFLXre2we5qR2MscD09tZNu4I+a9NF82NZ3SAeaw8cB08CAHiNCN6gB8cLENOSIKQlgtMI4M0GX1DN8b7X022Hu8YCT5u1T5cTwHzhvDF6FcPHRZD9JpIyhJDoG/akNlqtFv/1X/+Fhx56CGfOnEFXVxfMZjOmT58OYYL04SS9eL06VZa3tWvA8v3TcHFhWsEj5f/RDznVtleBq/zMiPY/GXE6DcBB1b3V094VNuBmjMHT1gXe1PskYaSfa7zidb6MFdqCbHgsHXDXWiCY9OAEHrzZ4EspZ9CN+9efsGoxXEerVDnGpfpmSPXNA2zly8WtyUyJdvUGZP/smOrcFZPNYEwBJ/CBCYqUHjsUAMzt8Q0I1WrA67RgCuC1OaABYFxcAjE9KeQx4pmQZIbXGn5GVH8ZQsjYG1awffbsWSQmJsJsNoPjOEyfPj2ojM1mQ3d3N6ZOnTriSpKxpy3Igr1Pv2upvhkeqw1icvDFnHmVoEfO2oLsER1fqrGEDrT9x/R44W3vgpiePKLjTCYcz0OblwmprjeFnetYNbQF2SGfSkn1LfB2O1TB9kg/13iiSPK53Nq+WTO9nT3gDFpoczPACTw0U4LnExjvtLkZSLllFTo2vB/xNprcDKR85fIo1ioy7op+15hpU8EbdXAdrw4qy+lEwKluBFJsTmjzM5H8pcuiWs9o0RZmQ2poCX9dFHhoCyfO95OQ8WxY3UhWrVqFl156acAyGzZswKpVq4ZVKRJ/DOcV+1pC/Rhg+3B/oO+0n+KW0LP9AJR++YlNF8wZ0fH7HydkmX6DpMjgTBfMVf0tN7XDvud40GyKcnNH0MQfYlbKhPgx93bZ4Dx6BvadR+A8XAln2Rn07DgE26flsO86CqnWEvR+TCTGJbOQetc1QZP6hKItyELGvTeFHa8RS94edepCMSMJ2pwMGOZNBycE/7T171IhJJuR/p0bwevGZ1cL3qCDYd50IMRrhcDDMG86+H6T9xBCxsawWrYZY5STdJLh9VqYlpbCtqsssMzT2YPO13dAm5cJMT0J3m6HLzDp17dTm5c54hZQJkUwkY5n4vUfjjbDomJ0/ecTKI7enNGu8jNwn26EbtoUcDoNPE0dkC2+Aa5cn/z65pULxvW4DMYYHAdPwv5JOZgSOl+04nTBVVkH5lUgpJghmCbmYFzjghlIveNqOA+fgquiRpXikeN5aAqyYZhd6Gs9joNAG0BwJpRzN0Sa7DQY5k2Dx9IJud3ae+1gDOAAzqSHYNBBSE+C8/ApCElmaAuzx2VgqslOhZBkUg0aH8+vh5CJath9tgdjsVhgMpkGL0jGjcSrl8NVUQtPW5/+2oxBqmtWzeDWF6cVkfzly8Z1UDaR8VoNkr98OTpefFe1nLnccB2vCbudbnoOTOfPjnLtoks6fBrsZMPgBQEwl+zL1LF89oQNYni9Bob506GfNw3M4YLXJZ/ro64Hf+4mK56+xkKSGZ723rEhcnNHIMMIJ4oQs1IgZCWDyV4w2QOp6iw4ncbXgs8YOI7rHWjd0ALDvOnQZI+vWUABXwu3vnR8TahEyGQTcbD91FNPqf7eu3dv0DIAUBQFFosFmzdvxvz5YztanYwu3qRH+revR9ufN0WUyo/TCEj7xrXQ5o7PbA2ThXHBDLCvXI6WJ1+LqLy2aArSvvG5cT2ZjdTYCqm8Bpp+M2NyHA8+0QjeqAf6pf5jsgxnWRWMy0vBIY6izlHGcRw4k0HVNz8eGRbOgPtMb75suakDminpqvzhHDhwGhGejh4wjxeMMQiiL4+6qg++V4HzaBWEJNOEvZkihIydYQXbHMdh79692Lt3b9jymZmZePjhh0dWOxJ3xLQkZD70ZTT/v41wHasOnfqN8w2cM1+0APqS/NhXkgyZ6fw5SKppgn3X0UCXkf44vQ7G82Yi4zs3BE3fPt64yqqClmmLpkJbkA0hwQj97ELYd5bBWXYa3j7Tmnu7bfB22oImhCGxZ1wyC11vf6pKYek8WAlNQZbvDw5QHBK8rZ3wdHaDebzgeF//bCHRHJypw6tAqrFQKzEhZNRF/Iu5fv16AL5+jnfddRduvPFG3HjjjUHleJ5HcnIypk2bBn6MZxcj0cGb9Ei4eAEM86fDXdUIT6sVzC2DE0UIKWboZuZCMBnAG6mFaDzRTk2H+PkV8HR0w13V6BuApijg9TpopqRBWzgFQoJh3Afa3i4b5LPqSZr0M/OgLZwS+JvjOAhpiTAsLYX903LA3tunXapvoWA7DvAGHcwXL0DPBwcCy5jihVR91pf6r/8Ukn1mwNVOm4JQIhmITQghQxXxr+ayZcsC/77vvvuwfPlyLF26NCqVIuMDr9XAUFoIlI51TchoElMTIaYmDl5wnJIaWtULNCI0+Vkhy/KiAF1BFhx98sZ726xRrN3gFLcET3MnFJcEjuPAJxggZqaM6dTpnFYD9BlkG7bMKEu85nzITe2Qm9rB8VzIKdoDzvXT1k3PgSYjdI7wiAZiE0LIEA2rieq+++4b7XoQQkhMKE51UCimJgwYqAr9pqhnXi+YosQ8uPW0d8OxrwLuqkYwrzp7Cm8yQF9aAOOikpjWyW+sJljheB5pX/sc5LPtsH16dOCMRBwH3Ywc6KbljHo9CCFkICN+HtzU1ISWlhZIkhRyPbV+E0LiSVBmHGXgNKYs5PrYDpC07TyCzn9uCz1GAoBid8KxvwKuynroS/OhzQvdUh8tYznBCicKSLh0ITR5WbBt2x/URYgTxXPTt4vQlRZGpQ6EEDKQYQfb27Ztwx/+8AfU1tYOWO7EiRPDPQQhhIw63qzOsuFp74Iie8CH6Yvuae5U/c1pNOD42AXbPR8dRte/dqqmJg+HOdxoffotZNx/M7Q56dGv3Dn+CVacR6uCA+4YTbCinZKKxM9dAKm6Ca6TdVAkX+pCTiuC43kwxsCL4zeDDiFk/BpWsL1nzx7cd999SE9Px2233YZXXnkFS5cuxbRp03Dw4EGcOnUKl156KebOnTv4ziYJxS3BfaoBUo0lMLsiZ9BBm58FfXE+DSYkJEa0+VnqhNEKg/tUPfSzC4NS+ikOF+R6dQ55MYa5mKWGFnRt2hm8guMgJidAkT1QbOqZFJlbRseL7yDrR7fHND1jPEywIqYkgDfp4e22w9ttBwDwOi2EJBM8Y9zXnhAyeQ0r2F63bh2MRiPefPNNpKen45VXXsHy5ctx3333gTGGdevW4dlnn8WDDz442vUdd5jXC9un5XAeOhXUzxI9DnhaOuE8WAn93GkwLp8zbqcOJmS84A06aKdNgVxRF1gmN7aCuWRoi7IhZqdBcbohnTkL18m6oO+tNi92eeNtOw4HtWjrZxfBuHCGLxc4ALnVCvvuclULvKetG87yahgXzIhZXYH4mGCF12ogZiRDSDCOaT2iTXG6aeZIQsaJYY3wKS8vx+rVq5Ge3vuY0j99O8dxWLNmDUpLS/Hkk0+OTi3HKSZ70P6X/8B5+HRwoN23nKLAdbwGbU+/CcXhimENCZmcDPNnAP26gnjarXDsr0D35t3oePU9OE/UgvUbcCdmpcZsynZvjwPOI6dUy/TzpsO8Ym4g0AYATUYykq65ICgdoX1nWUzqSWJPtnTAtqsMUq2ldxbMWgtsu8ogWzrGunokjihON1wnamHffQz23cfgOlEbNEicRN+wgm2n04msrN4BOFqtFjabeiT6woULcejQoZHVbpzr3LgNrpP1EZeX6lvQ/uK7gRuX0RZJ6q1opOcik4u32w53ta9V2F3VCM8gWSrGgpCSAE1pga+7gd0JxSmBeQf+3vkmuymC1+aAbGmH1NgKuanNl488CtynG8A8vf2fOYGDcf70kGU5UfDdQPTd/kwjFEplN+EoTnfovvFAYCZMCqYIQDdl8WRY3UjS09PR0dH7QWVlZeH06dOqMlarFd4BWnMnOnd1ExwHTqoXcoC2cCp0BVkAz0NqaIF0uhFM6b1ouk81+KaEjsLj37FKzxUOPQadOBhjcJ85C8fByqBsEAAgZqXAMLsI+jlFY1C7XopbhvPIaThP1ELu7AbPcfB29gS+g7xRDyHJFLSdkGyGZmoGnAdOwtvd5zvEC5BqLNDmZsCwuASmc7mcR6Wu/YJ4ISN5wO+GJjdDvYABit0Fnm6gJxSpxhI+6wtAM2ESAJHdlAlJJvq9jZFhBduzZs3CqVO9jzeXL1+Of/3rX3jnnXdw2WWX4cCBA9i8efOkHiBp/+So6m9OIyLx6uXQZPUOrtJNmwrPvOnoeuczwNN7Y2LfVRaVYHss03P1J1s6gi4EXqsNUkMLDPOmQxPDQWhkZBRJRucr76H7vX1B3S78PM2d6GnuhPvMWZiWzx6TC7zU2IbOjdvgtfYElgkJRmimpkPpcfhauB0uKE43eFEEeA5iVho02alwHq2C62T4zEtSQytkSyfAgNSvrh6dWTb75/H2DNx40fem3S+WWVNIbEQyyyXNhEnopiy+DKsbyeWXX44TJ06gsbERALBmzRoYjUY8/PDDWLx4Me655x54vV489NBDo1nXccPXeqbua2lcMEMVaPuJKQkwLZ+tWuY+3QhPnxnrRos/PReEEB97jNJzAfQYdCJhXi/an38HzqNnIiovN7ahbd2/Y969QbZ0+MZE2JxB6zieB59ggCY71TctfXYa+AQjdDNyoZ9dCKmqEUpP8HahOA+dQvv6LSED36ESktVPmTzt3YEMG6FI1U2qvzmRBx+j/uXxhosgxd+o3BCNgUhmuaSZMAndlMWXYV1tbr75Ztx8882Bv/Py8vD666/jhRdeQH19PaZOnYqvfvWrKC2dnPN4e1qtqr6WAKAryQ9bXjdtKhx7j6uWyZb2qEyZHQ/puSbyHbcieyBVNcJ95iwUmwNMAXi9FpqcdOhnFUy4FI89H+yHuzJ4XAKnEX3dlrrtQT/8Uo0F3e98huQbL45JHZmioP3FzSFv4HiTAczjPde9ioETeIAXwGlFuE81wNPcCa8tuE82J4rgjXooruDJvFzl1bB9fAQJl543onrrivPA67W9x2CAfc9xJKxaHDR7peJwwXFE3ZVPP2daREHnRMSb9CFvrFRl+gwyJWSioZuy+DJqt/b5+fn4xS9+MVq7G9f6z/LGaYQBA1lO4MGb1Bd+5o7el2Cs03NNxDtuxhjsB07CsfdE0Oev2BzwtFnhLDsNfUkBDItKIJjHf4sjkz2wfXxEtYzTamBaPhu66TngRAFMUSCdOQvb7uNgrt5g1777GBKvWhaTgMd1ojZoYhoxIxma+UVIKJgK9/EaeLsdkGotkGr7tA57FbirGiD2ma6d02mhL8mHmJni65utEcBpRDj2nlCl6LN9dBjmixeMaEp3XquBcVmp6j2Wai3o3rwbhgUzoMlJB5O9cJ9phPNQJZjTDU7svaSbL5o37GPHs1BjPZR+10sxMwWetq4Bu8yJmSnRriohhAAYZjeS0tJSfP/73x/tukwYvF6r+pvJ3gEf/ypuCd4u9XpeP7FaQPuaaHfcTFHQ+er7cOwJDrTVBX1dhFqffA2ePv2GxyvH4dNQHOrW4sQrlkJfkh9oUeV4HroZuUi85nzVDOdM8sC+ryIm9ew/fkJINCHxc+dDyEgODGbkdRroi/OgLZoaKOe1O6A43WDnAjaO52FaMguarNTAdpwgwFBaiITLFqmO4bXa4CqvHnHdzZeeB65f7n3Z0o7urXvQ/rd30PHyFtg/ORr0OeimTYV2es6Ijx9vwmVXcJWfgaez9zvFazXQFmSH7TKnLcimOQ0IITEzrGDbbDZjypQpo12XCUMMkTXAOcAPr/tErXriCp4LzixA4lb3258GZ54ZgKetC+1/+c/Agfk44DpRo/pbk5MOzZS0kGU1aYnQFqqvGe6K8AMORwvzeoO6uRgWzFC1APelLZzi60oCgDmlc//33fiJ2WlhW+K1+VkQM5NVy1yj8PrElASk3XVNoE4RbZOehNSvXTNqWVHixUBjPZhXgVRrUY0FEFMSoJ9dCDErBbzZAN6oh5iZAt2MXCg9Drgq6+E61QC5qW1U+tgTEo8UtwR3XbMvFeupBshN7XS+j4FhdSOZP38+Kipi0yo1HnGiAOPSWarHv64T1RBSzNDPKgj8CDLGIFU3wXGoEuB6f0wN82dM+NnPJgpPWxd6dgTnk9fmZkI7bSo4rQj5bBvclfVgfbJJyGfbYf+sHOaLF8awtqOrf59YbV5WmJLn1udkQK5vDfwdrfzUfSmO3pZpv4FuZHlRgJBiBhQW+EFiigJO4FXdSULR5GfB09ib9nC0Xp++tABpa66H5XfrwXoGvkETs1OQcf/NE/L6EclYD09LJ7S5vTN88lpN4G/F6YZU3wz3tgPBXf0MOuhn5cO4eFZU6k5IrHnau+HYVwF3VWPQpHq8yQB9aQGMi0rGqHaTz7CC7fvvvx+333473nrrLdx4442jXacJwXThPNh2HultsWa+x9muo2egyc8Cx3OQ61vh6fRlHeHE3mB7ova1nIjsnx5VP5XgOCResRTa/N7AU1c4BYYFM9D9zmeqx/22XUdhWrlg/LZA9qv3QLOk+tarA6WR9GeOVKhjMHnwevb9TPz/7F//IP1S841m2j39zFyk3nYlnEer4D5eG7hu+A7ku9HRlxZCV5wHIXHiBdpAZOM4FHvoGXhdJ+tg234QfJibEOZ0w3noFNyVDdDPzIVuAnbBIZOHbecRdP5zW9inp4rdCcf+Crgq66EvzR+0oYSM3LCC7V27dmH58uX48Y9/jJdffhnz589XTd3ux3Ec7r333iHvf9OmTThw4ADKy8tRWVkJWZbx6KOP4qabbgpZ3maz4U9/+hPee+89tLa2IiMjA1deeSXuv/9+mM2xm6SlL01mChKvOR/dm3erlnu77fCWh0+TZrpwLl3oxwnGWFC/Y31poSrQ9hNMBphWzEPPB/sDyzytVkg1FuiKxmeXLDElAX17CktnzkI/uzDszUP/1HT9U9tFA2fQgtNpVAOO5ZomCGHy2CtuCV6rDWJKgm+Ap+zx9ftlgKe1E9qc0K3ivqdUZ8Ebewe+jvYEUbxWhKG0EIbSQigONxSXG+A58EZ9YOKa8XrfFomIxnqEyEXuOlkH284jYBE8OWduGW3PbUL6d2+ErnB8fi/J5Nbz0WF0/WunuhEoDOZwo/Xpt5Bx/83Q5gTHcGT0DCvYfuqppwL/Pn78OI4fPx6y3HCD7SeffBKNjY1ISUlBZmZmIJ93KA6HA7fffjtOnDiBCy+8ENdeey0qKirw4osvYs+ePdiwYQOMxrFp6UlYvQTM6UbHhg8iKm9cXILkmy6Jcq3IaFFszqCuFAOleNTkpINPUGch8Vjax22wbZg/HfY9vd99T3s33JX10Id4D9w1TZAt7aq+0oYoTNzUH8fzMMybBsf+3j71jqNV0M4IvqFlYHCdrA/8SPFGHRSnBE6nBXNJ8LRa4WnvDtmdxHWsGt4uhyrYNiycOfovCL11m2hpJKPBY7UFDZAFAHCAmJwA5lWCBq8z2YuOF99F1iN3xO0gSm+3Ha6KWsiNbVDcki9XvNkI3fSp1FgziUkNLejatDN4BcdBTEmAInmg9EtlytwyOl7cjKwf3QZOmJypQmNhWMH2+vXrR7seKr/5zW9QUFCAnJwcrFu3Do8//njYsn/9619x4sQJ3H333fjv//7vwPL/+7//w9NPP42//vWveOCBB6Ja33A4jkPSFy6C3GaF47NjqtHyfQmJRhiXliLltivGb5eCSSjUIzohIXxKP47jICSYANbb5KBI43eQpK60AGJaIjztvV0abLuOwNNmhX52EYRkM5RuO5wnauE6ph4gzJsNMQm2AcB80XxVsM3cMqz/2gluVh6Mc6cDjPluFGrOwtvRDfC+HxzeaIAgecAJXKCRyHm4Epr8bGhzM8AZtFC67LB9fASuyjrVjYRmahq04+wmitNqAMfAk0lx42zqd9fx6qDBYPqSfBgWFUM4N+GP3N4Nx57jkM/2GU/QZYfzUCVM58+JaX0H4+12oHvLHrhO1Qe1XHq77ZDPtsK+5zhMF8wZ313UyLDYdhwOOi/0s4tgXDgjMLhbbrXCvrtclQ7V09YFZ3l1VGauJj7DCraXLVs22vVQWbFiRUTlGGN47bXXYDQag1rQ16xZg1deeQWvv/467r///jG96Ohn5EIzJR0eSwfcNU3n+u0y8AYdtHlZ0ORmQDDp6cI4zvA6bdAyb3sX+CmhH8cxRYGnvUs1WVH/NJHjCcdxSLhiKTr/8WHvQubLa+06MXAmjoRVi2M24YomPwv6WflwVdQFljGXBGnfSXQcroK3owdMCe5+wGlF6OcUQW5o6d1OUSDVnIVUc9ZXRhRDdodJuHLZuPs+C0nmc5P7DFxmvFBkD9ynGlTL9CX5MK9coFqmSUtE4tXL0PX2p/B29DaI2HaVxVWw7WnvQuszb8F9pmnALgJMkmH/pByazBQkf/nyUTkPQ+U2j+VEaGRw3h5H0MzVhnnTg2ao1mQkI+maC9C1aSe8fWbGte8qo2A7iqI/QimKampq0NLSgkWLFgV1FdHpdFiyZAmam5tRWxv9FGOD4TgOmilpMF8wF4mrFiNx1RKYV8yDNi9z3P0oEx/epIeYnqRa5jw2QIrHqrNBkxVpC7KjUrdYMS4rhfniBYMX7LvNkhKYL1kYnQqFwHEcUm+/CmLWECYx4TiYLpwH/awCiFmp4cuFuGFIuGLJuPzR0haGyUvtJ/C+MnEo1NTrHkuH+umTwMFwXnHo7Xk+qNuP3NgGb3f0M+ZEQnHLaHvu36qbgcHYdx9Hz3t7R3zscLnNbbvKIFs6Rrx/MjrcpxtUM1dzAhf26SEnCjDMV69zVzWqUmeS0TVqM0iOBX8QXVhYGHJ9QUFBoFy4MoxFMIpghBgb/DiDlRmNfYymkdRnONuG3xcL/H8s3iPjBXPQ/Z9P/HuHVGOB/cBJGM6bqbqJkhpbYf+sHL4mKd/xtUVTIWalxOwzG8xw37/E6y8CZ9Sh48V3Bz4A5+unnXzLqnP7it3r5ow6pN97E5p+8TzkuuZAT55QVeC0AhKuXAoxMwXu4zXQz5sOxR7cPx8cp27Z44DE6y6E+dKFUXlt0T6/Ob0W+rnT4Co/EzTQkBMF6OdOA6fXjtn5OtBr44w6MJsDYCxQRnGqu8SIqQngTfqw+9DkpKPv9xPwdc3oO86i7/UmlmwfH4an1f/Y33dsThSgK86HZkoqFMkD6cxZyI2tfcowdG/dB+Py2cN+IuHPbR5q4Ck8XjiPVoFPNFILdwSife74bgx79y1kJPsGh4c5npibAdX5zgCvzQkuZVyHhWMikgbTcf2u9vT47vLDZRzxL/eXC6WrqyvqF06XzQbmHKQvpOKBx2qN6j5G00jqM5xtB9vG4XBAdjrhsYefqXOgOg2XUjwVkqKAyR54ZRnMo0Deexy28iqIeZmARoDX0gFvW2+KR8XtmyxFs6AI1hh9XpEY0Tm2rBjazi5IZWcgn24Ac/e2KHJaEZppU6ApzgPLTkVXd3fw9jEirFoIdqYJcmU9UN8Kj8cDxesB8zLwJh2EqekQ87Mgp5ggyxKUjER46luBKWkQEwzwNLbBa+kAUxh4swFMUcBrBYgzc6CdNx3e84rQ1dUVlbrH5Bqg56HMzYe3rgXKuYGDfKIJQn4mvHoeGMPzdaDXr5h0kCUJHFOgnLsGyC4XZH/LNseB02thH+D6wGQPZFkOfD8BoNvWA8Ea/DMZrc84ZL0UBV3b9wfq5ZVlcCYDDFcsBgw6+NsihZxUsLpmOD8u870P58q3fLgXhssXhdn7wOTKBni6Bv6+uo6egqY4d1j7n4yidn1w2OHuc+7C4x3wfFfsruDzvacbPDdwalQSLCVl8Kem4zrYHg1JSUmDFxohm9kMhR/4reaNepiTk6O6j9E0kvoMZ9tw2zDG4HA4YDQaIRkM0JlMw6rTsCUDuq9eAevGDyFrNGD+C5XsBTvjS3XHA+DPPebmRAEanRb6OUVIXXleXHUhGuk5Jk7JhJKUCHbRAig2B5jkBacRwCcYe6dFj+E5GrKO5gQoMzRg0/Ng7+qGnvFwV9QBYIBWA44DOJ0WBv95ZDJByUiFQ4EvC0n+FHB6HQSjDuB5cFoRvNkAjuOi/tpieg3Izhy8TIwN+PpNgKekEJ7mDujPfXZyRgq8GtGXiSE1Ad4uO3QuD8S00Nd8V2U9NBoNNH3GYqTkTVHNGsoYQ1dXF5KSkmL23XWdrIPDKQPn6iVrNEhctRhieojuTaXTILR2Q6puCrwOrrwOSTdeNqz62tlZeAe5pgqMh2kMv9PjRbTPHWduNrx9zl25swN6Bb5B+SG4qpvV57sgIGVqdszG0kw24zrYTkhIAODLsx2Kf7m/XCixuGBy3ODHGazMaOxjNI2kPsPZdvD3hxuz98h8/hwwl4S2P2+KoDQHfWkhUu+8GnwMJnUZipG+f/51HMeBTwz9tCmW5+hgx+c0IkSTCZ5Eo6ovvf81+Ak6LbRZvmm+I913NMTbNSDWBnttmtREaDKSoc3NhLfLBiHJBOfBSgAcOFGAt8sBx94TSLx6edBkR4pLgvPQKQDcuf8AfUleIGNJcF24mL3PvqwRfc7HlARostLCltfPKoBUbQls47XaANkLbjhpDCV58NcZSRkSEK1zR1+cD16vg+I611LNAMeeE77B6P3Pd4cLjrIq9D3fDfOmBxqFyOgb1++sv092TU1NyPX+Pt3+coRES8Kl50FuaIHtk/I+fSvVeKMexkUlSLv78zGZPZGQyYbXaaAv7b3eKz0OdG/pHSQon21D138+gWHhTGhzM8G8XkjVTXAcOgXF5lClbzRdND+mdQ+n/2Q+QuIgLc0hWjKZWwbiNGc4GR28TgPjslLYPj4SWCbVWtC9eTcMC2ZAk5MOJnvhPtMI56FKMKdbdb6bL5w7FtWeNMZ1sF1YWIjMzEwcPHgw0JXAz+12Y//+/cjMzKRgm8SENj8LyenJkNuskKrOwmtzAkwBr9dBk5sBbX4WBLMhrgNtxS3D09IJxeGb9po36iFmpsTt5B6EDMS0Yh5sHx85l27Vx9NqRc/7+wbcTpObAf3swijXLjL9Bx962rrAFCXkdURxy3AeOwOPtQcQeHAaEUKCEbxh/KYYJZEzX3oe7HuOq2fMtbRDtrQPuJ1ueg60NBlSVI3oV//999/Hgw8+iOuuuw5XXHFFYHlVVRX+8pe/oLm5ecQVHAjHcfjSl74Eh8OBp59+WrXuueeeQ1dXF770pS9NqEdciluCu64ZrpN1cJ9qgGxpD5q0gYwtTXoyTMtn+1I8rl4K80XzoSucEtdBNgB4OnrgOlHjC7bPzY7paen0LQszIRMh8UxIMCLt7uvAaSLvhyokmZD2zWvj5vuqzc9S/a3YXZBqLUHlPJ2+76/z8CkwyQPmlnxjJ2RZNfEUmbjElASk3XUNuIFSePbfJj0JqXddPaHipHg0rJZtRVHwve99D1u3bgUA6PV6uFyuwPqkpCQ88cQTUBQFa9asGfL+X3vtNRw4cAAAUFlZGVi2d6/vceDq1auxevVqAMDdd9+Nbdu2BWaSnDNnDioqKvDxxx+jtLQUd99993BeYtyR27rh2H8CUtVZMK96tDBvNkI/Kx/GxSVjVDsy3ilON9zVZwFviBs3rwKp1gLeqKcpwsm4oyuagox7b0LTr14MTMoSjpiRhIyHvgQxOfw4n1jT5GVCk5MOubEtsMz2yVHfE7Mpvr7biluGVH0WrsoGeNq7ArOgAoCuOA/Oo1UQkkyUom8S0JcWIG3N9bD8bj1Yz8AzFIvZKci4/2YICcYBy5GRG1aw/eKLL2LLli245ZZb8PDDD+OFF17AM888E1ifnp6OxYsXY8eOHcMKtg8cOIC33npLtezgwYM4ePAgACAnJycQbBuNRrz88st46qmnsHXrVuzduxfp6en42te+hvvuuy9ospvxqGfHIVhf2xZyenAAUGwOOPZXwH2qAfrSQmhzM2JcQzLeSTUWQBkgBaZXgaelE2LK+JlBkBA/bUE2Um67Aq6jZ+A6UQtPm7V3JQdoczKgKy2EviQ/rgJtwPcE17xygWqmVuaS0PXOpxCzUqCZkg65qR3uU/VB/bs5jcY3sNerQKqxqPqzk4lLPzMXqbddCefRKriP18LT2efJBgdo87KgLy2ErjgPQuL4j5HGg2EF22+99Rbmzp2LX/7ylwBCjxAvKCjAjh07hlWpxx57DI899ljE5RMSEvDII4/gkUceGdbx4lnPtoPo+s8nA07P66fYXWh7+k1kPHBzoMWDkEgM1uIHINCPm5DxiBcF6EvyoS/Jh+J0+ya9OZey0T8mgePj81G6ceksOA+fgquiTrXc09wJT3MnvFY7mCd49j/zRfPAa32vLZLvOJk4eK0IQ2khDKWFUBxuKC43wJ8738+dE9RzJHaGFWzX1tbitttuG7BMcnJyXE3YMR5Jdc2+QLs/joOYkgBF8kCxqacTVlwS2l/YjKwf3RY3fQ5J/OvfIhayTJgnK4TEAmOA3NQO18k6eNqsYLIvj7uQmgh9cf65GSAjwxt046pLBcfzSL3rGrQ//zbkphCD3fqP2+EA04q50PUZ9BbJdzwUxe6E62QdpPoWKE4JHAfwJgO006ZANyOXupaNA7xRR5/TGBtWsK3X68PmtvY7e/YsEhMTh1Up4tOz/VDQMv3caTDOnwHeqANjDJ5WK+y7j8HT0ptuztNqhet4DQxzp8WyuiSMvhk+OJ0WYkYKtIXZ4+rHnpCx5K5pgvWf2yA3dwSt81ptkM6chZBkRsKqxWNQu9jg9Vqkr7keiiTDcbAy7GyaQnIC9KUFMMwuGtHxFKcb3R8egOt4je9Opw9vjwOypR2OfRUwLCqG6aL51LhDyACGFWyXlpZi165dkCQJWm1wSiGr1YqdO3diyZIlI67gZOXttsNZdlq1zLBgBkxLSwN/cxwHTWYKkj53AaxvfQzF3vuY376rjILtOODp7PFlDjg38JDzeCHVWiA1tMAwbzo02SFmgSOEBDjLz6DjpXcHzajh7bKh6z+fQF+cB+OSWTGqXWxxogDT4hLoi/Mh1TRBbmqH4nJD1mvBiQLE7DSICUZw+pGl+vPanGh79i24T9YHBdp9MdkD54FKtL+wGWlfuwacQLMPEhLKsG5F77jjDjQ1NeGBBx4ISu9XV1eH++67Dz09PbjjjjtGpZKTkbuyXjVgjRN5GOZPD1mWE4Wgda5TDWAeb8jyE523xwF39VlfesSqxjFLW6e4ZVWgreJV4Dxa5es3SggJSapvRsf6LWCeCNObMqDj7x/AfbohuhUbY5zAQzc9B+aL5iNx9VIYF86EfmYexFHIKsEUBe1/exvy2YFzM/flKq+G9Y2PRnxsMno47eBzI0RShoyOYbVsr169Gvfccw/WrVuHSy+9FAaDb0rbCy64AFarFYwxfPe738UFF1wwqpWdTLw96r7YYmYqeF341gptXibsnx3rXaAwKA7XoLONTRSMMciWDrhPNUBqbA1aL2akQD+7EIZ5sWvt97R0hg60/ShDACFhKU432v7yNqRzKe+8PQ5wWhFishn6knzwqYlQrDa4TtapB/8pDNZ/7UTm92+h3MHD4Nh/8tx07704gYduRi40uRmAR4G7tsmXwagP+2fHYFoxj7JhxQkhyQyvdeDuvkISZZeKlWHPIPm9730Py5cvxyuvvIKysjJIkgRFUbBy5UrccccdWLly5WjWc/Lp1/+tf27t/kK2Yk+SHxrm8cK2/RCch09DSArdsuNp7YTto05IZ87CtLQUvEkf9XpFkr2DMgSQSHFaDeAY+EnIRGmpki0dsH18GK7j1YFlitMNTZIJ5tVLoElL8i3Mz4J+3jTYPyv39S32b9/YBqnWAl3hlBjXfPyzf1Km+ps36pH0+RWqhhvdzFzITW3o3rJXve2nR6H98uUxqScZmLYwG1JDS/gGH4GHtjA7tpWaxEY0XfuFF16ICy+8cLTqQvoQktV3nJ6WTnhtTghmQ8jyfVsZmMcLxeWGs/wMOI6DkGSesAPyGGNw7SwD12yNqLzc1I7WP/8LGfd9MepTkEeSvWO4GQKigROFQTNMTpRgbjyaLC1VitMN59EquE7UqpZzWg30pQWQ61sgmI296fo4Dqbz50CqawFzSYHyjj3HxyzYVpxuSDWWwM30eLkGy5YOSHUtqmWmC+aEfEKqmZIO/bxpcJX33hA5DpxE8k2XgBOp7/ZY4w06GOZNh/NoVXDALfAwzJse9+fjRELDh+OUviQfXN9gkAH2PcfBQgxWUexOOM4NplQcLsiWDoipSVC67L6R+rUW2HaVQbYEj+Qf71xHq+BpaAtazmlEiOnJ6vfwHLmhFV2bdsaieuNKJK39EyGYG6+0hdnAQNMwT5CWKqnGN87B22+shWZqmi/jxbkJlvrieB76mbmqZWN1vZMtHbDtKoNUa4HXahtX12C5SX0t5XQaaAvCn1P9u8AxyQNPB00NHy802akwXzQf2oJsCMlmCMlmaAuyYb5oPg3Oj7FhBdvbt2/HfffdFzQ40q+5uRn33XcfPvqIBkwMF6/XwrhUPaJeqj6L7nc+g1TfAqYoUNwSnMer0fmvnWAuCczrhaejB2AK9LP79QOegAPymKLA2adVBfAF2aYL5yH1tiuRfMNKpN52JRIuXxyUY9Sx70RQv/jJTsxMGTCY4yZIMNcfpxn8AV8kZaLN31IV8jOaQC1V/tbg/l3jhD7f4VBdtPh+ra9j8dTI3yo/XgdFM0n9NE5IMAyY0o836oF+vRX774OMLd6gg760AKbz58B0/hzoSwsmxHVivBnWL8irr76KlpYWZGVlhVyflZWFhoYGbNiwAZdccsmIKjiZJVx6Hhx7T6guXrKlHbIl9ChxxeEGb9BCMzXdFzj1N8EG5Ek1lqBcswmrlqgG6HA8D920qRBSE1Sj5ZlHgWPviQmdl3eoeK2vFStkBhWBh7ZoyoS8SPNGPRSbc+AyMejjHwlNdiqEJNO47KIQKX+Q3P+plLfbDs0U38Q1obpoeVqtqr95fezfD3+rfFhxfg3ufw55rTYosgd8mJtNX+NO/32MLO0gIRPRsFq2KyoqsGDBggHLzJ8/HydOnBhWpYiPmJaE1LuuBiKcQph5vb6JHS5bFHYU/kQakCc1qPsWCqnmsCPhxeQE6KblqJb17xNKADElAfrZhRCzUsCbDeDNBohZvkwuYurEnKRqsBZ9CDzErPh55DpZWqrEjGTV33JTBxR36NZqxS0FpfvTFoRuDIqmSK6v0bwGK7IHcqsVUn0z3KcaIDd3gPWfXXIAmvxMVUs18yhwn6wLW959TP1kkTcbIKQkDLnehEx0w2rZ7urqQmrqwD8+KSkp6OzsHLAMGZxhdhHS7/kCLI++AmYb+PGcmJyApM+vGPCHN54G5I0U65eZQROqNb8PbW6Gr9X2HOpGEhqv1UCbmznW1YgZXjdwi75h3nS4K+ugDJIJhIwu/cxcOA6cDEyqwrwKnEdOwbBgBoQ+rd6Kw42eHQfBJBmc2PuTZjx/TszrHMn1daTX4FBZabxdNrhPN8JzLrjmRBFygy8FqpBghH5WAQyLigfdt5icAH1poSqzi2N/BfgEI7T5WYFGHKYocB2rhquyTvWemy6YQzNJEhLCsILt1NRUVFdXD1imuroaSUlJw6oUUdOX5CP19ivhLDsN94la9SQtHAdtfhb0pYVQnO4J2cIVVr8WfzbQ41sEp0/kInxiQCY+MSUBvFEPT0snmFfxDSbq0z3DXRm+dY9EB28yQFs0BdKZs4Fl3i4b7DvLIE5NA3PL8FhtkOuag77b+ln5g958j1d9s9IwxuA8chr2T8vV70GfbCDeHgfs+07AdaoBhjlF0ExNH3D/5ovmq4Jt5vGi5/19ENMSocnNBPN4IdU2Q7GpGys4gYfpgrkjf4GETEDDCraXLl2KrVu3oqKiArNmBU+LW1FRgW3btuHKK68ccQWJD68VYZhdBMPsIigONxSXG+A58EY9+HPp2Cb6rGn98f3SIMpN7WCMhe1CI9U0qf7un16RTG68TgNtXiZ4ow6mMWgVJcFMy2fDY2lXPVVgTIGnqQOOMF1KeKMOyTdN3LFCffMnOw+fOtf636cAx4VsdFFsTrQ+/SYyHvzSgDciuln5MC6dBce+CtVyT3s3PO3hM40kfn4FROpCQkhIw3re861vfQscx+HWW2/FU089hUOHDuHs2bM4dOgQnnrqKdx6663geR733HPPaNeXwPdjIqYmQkxOCATak4Hi8cB1sg5dm3ej858fwnmsBlJNE7xtVjDJA8XmgutETchtpfoWSA3qmSUNC2bEoNaEkOESTAYkXn2+L+tFBDi9Fulrrg/q7w1MnOmr/Vlp5JZOX6DdF8eBTzBBk5Ec1BgB+LrcdLz0bsgUsr274JDy5cthWDA94jolXrUM5ksWRlyekMlmWC3bJSUlePzxx/GDH/wATz/9NJ5++unAOsYYjEYj1q5dG7LVm5ChYozBcfgU7LuPB/V3ZAqD0mmDu8sBTuRh+7gM3vZu6OcUQUhJgNLjgKuiVjXxAnDuByuCPoyEkLElpiYi+fqLYP33Lsj1rWBK8Gy5HO/LlmO+eAG0+aEHRk6kSYE02alQHC7wZqMvM4tGBMeLMC0thfG8YvAGHRhj8LR0wv7ZMXjarIFt5bPtcJ+sh35Wftj9c6KA1LuugWJ3wb6vIqjLiJ+YlgjTBXORePXy0X6JhEwow04ee8UVV+DDDz/Em2++iaNHj6KnpwcJCQmYP38+brzxxkEHUBISCcYYrBu3wf7pMTCPeoAoxwFiehK8Da0AGBSHG1JtE8AUuAYYQQ8A5svOm1RPBYjPcKc8n0xTpccj3mSAYf4MaIumwtPUBq/DDU1mCjhRhJiWCN2MXPBGXVA+/b4m0vTVno5uuCvre7tteBUY5k+DcVFJoAzHcdBkpSLx2gvQ9eZHUJy9s2vaPikbMNj2b2+YN83Xb76hFXJDy7kc4Rx4kx66aVMhZiRDiJO0mITEs2EF2//617+QlpaGlStX4u677x7tOhES0LN1L+x7joddLySawKc4gS47AIDJXkj1rdAWZocdAGlYOBMJq5dEpb4kvg23dXMitYqOZ7xGhDY/G5xOA8PcaUPbdgJNX+2urFf10+Y0AvTzQnf74DUi9POnw7GnNxWvu6J2wPEtfXE8D11+FnRhnhgQQgY3rD7bP/nJT7Br167RrgshKh5rD7rf3xe0XJubAfPKBUhYvQT62UUQslKhyegd8MMkGUq/qZ79DPOmIfWOKyP6kSETz3CnPJ8sU6VPdBNl+ur+aUs1U9LCTjwDANoc9fwDzKMETQhGCImeYbVsZ2RkQJYnTr5mEp/snx0DlL7NN74ZInWFUwKLtAXZUGZOgfejMl9Dj6LAa7XDY+2BkJoIcL6ptnUzc6GfVQBtbgblgZ3EVK2b/Q3QujmRWkUnO/+kQONZ/8YC5hkk7aknuJ876DpISMwMK9hevXo1du7cCUmSoNXS1Kxk9DEGOPaqZyDVF+erAm0/3qCD/sL5kGpbICQZIaQnA7IHpgvnQpubCd5soACbBGiyU8EnGuE6egoC48EhsinPJ8NU6WR86J+2VG5q8w2YDJO1pX/aU16vBaejMQaExMqwgu0HH3wQBw8exL333osf/OAHmDlz5mjXi0xwcqsV7opaeFo6oUgecKIAIcUMXXG+b8p12RPUR1ZXnBd2f77pxX0/NBwHQCsCCoOQaIrmyyDjFG/QQVOcC1Ny8pC6FE2EVlESXZ7OHrgqan35wd0yOJ6HkGiCbkYOtEVTR+UY+tmF4DQCmHyuxZoB9j0nYL50YdD57O1xwNkvG5NhwQzqSkdIDA0r2L7hhhsgSRJOnDiBXbt2QafTITU1NejLy3EcPvjgg1GpKJkYpMY2dL7xEeTG1qB13i4bpBoLeLMRpvNnB60fKHDmOA6CyeCb7OeckI9OCSEkCjwd3eh6+1O4+8x46efttkNqaAG3+zjMF82D+eKFIzoWb9TDuKhENXjcXdUAxemGYcF0aKakgUkeuE83wHH4dNBU9qaL5o3o+ISQoRlWsM0Yg0ajwZQpU4KWD/Q3mdzcVY1o+8t/4GnuHLCcYnOgZ9tBcAIPTXZaYLmnrQvavMyQ2zDG4O2yqR6NctphZ7Ykw0Dp8chkJTe1o/XZf0GqaxmwHHO5YdtxGNqp6Ui87sIRtS6bL18Ex8GTva3bAOSzrZDPBjdk9KWfUwhtbujrKCEkOoYVjWzbtm2060EmOE+rFe1/fRsszBTL/XE8D4+lA5xWAzE1EQDgOlYdNtiWqpuguGQIfYJtcYApiUloit0JudUK5pbBiQL4JBPEtKSIggJKj0dGE1MUeFqtULodYIoCubnTl0t7gKwbY8Frc6Jt3b+h9ISe+CWUnu2HIKQkwLxywbCPq8lMQertV6H9xc2Rb5OditRbr4i4PN1AEzI64uuqRSasrnc+g+KSVMt4ox76WfkQ0pOh9DjgrqyDp707sF7MSoVU3RSYuEFqaIF93wkYF5eoBjx6mjvg/lSdi1tMT4aYlhTFVzSxSI1tcByogFRrUeXvBQAxJQG62YUwDjK9/USaNISMHU+PA/bth+CqrIficAGiAN6gg7fLDg4MYlYqdDNzx7qaAbZtB4JuMjmtBvqSfIjZqWBON9ynGyFb2lVlut75DMYls0Y0uNYwfzrS774Olt+/GjTpV3+a3Axk3PfFsIMoQ6EbaEJGx4iDbY/Hg5qaGvT09MBsNqOoqAiiSDE86eXtsgWlWtPmZyHh8sXgRCGwTD+7EM6DlXAcqgQAiBnJkOqb4emTM9t55DRclfXQFU0Bp9VAbmyD82wrNP1au/SlBTQAKAJMUWB962N0bdoV9sfa09kDzydHIZ1qgOG84t5Z6/qh9HhkpBxHTqPt2X+ppwf3eOB1S770dgIHuakNnpZO8Eb9mA/0Y7InaNItMSMZiVcvB6/rzdSln1UA18k62HYe6d3WLcOxv2JErduA77qZevuVcJZVwXW8JpApBwDAcdAWZkM/q8A3y+YQZ3ukG2hCRsewo+Kuri788Y9/xNtvvw2XyxVYrtfr8fnPfx7f+973kJJCj/EJ4DhQqcqXzYk8zBcvVAXawLnpgRcVQ6pr9rViiQJ0RVPhaelUtZ4wpxuu4zVhj6fNzYirlq94xRiD9bXtsO8OP0NnX56OHrQ9/SYyHvgShERjyDKUHo8Ml7uuGS3/+3ewUN0WGIPidIM3an1PtRjg2F8BADAuHLtsWM7yM1D61dd8yUJVoO2nL8mHVNcMubEtsMy+5/iIg20A4HUaGOYUQT+7EMzhgtclgxN48CZ9oNvNcO5J6AY6PMXppuscidiwgu2uri585StfQU1NDZKTk7F48WJkZGSgra0N5eXleO2117B3715s3LgRycnJo1xlMt70f3yqLZoKXh86PzvHcdCV5MGx1/dDKmYkgzfr4T03HftgNFPTYb58MeXVjoDzYGXIQJsTBAjJJngd7qBZ5jzt3ejc+CHSv3Vd2P1SeryJI1YBBWMMbX96Q5VNyI/T68CLArw2J5jkAdfn2uHYXwFNVio0U9KCtosF2dKh+lvMToGYHPrJDwDoSvJVwbbH0hHxtOmR4DgOnMkA3mQYlf0BdAMdimzpCLoB8VptkBpaYJg3fVzNRkpiY1jB9jPPPIOamhqsWbMG3/72t2Ew9H6xXS4X/vznP+PPf/4znn32WTzyyCOjVlkyPjFJ3T2BD9Mq6ieY1es1WakwX7wQtk+OwtPcEXIbzqCDfmo6Eq9aDm6gabUJAF9w07P9oGoZJ/AwLiqBrrQAvFYDxhjkhlbYPyuHt7v3Zsd1vAaypYN+UCa4WAYU7oo6SGfbVMuERDN0JfkQkk3gwEE80wi5xQqlR33j7Sw/M2bBdv9rm5A4cP9lIUF9bWNeBVAUQBDCbBEf6Aa6l+J0h27pBwCvAufRKghJpkl7I0JCG1aw/cEHH2D58uX4r//6r6B1er0eDz30EA4fPowPPviAgm0C3qBuxfa2dg1Y3tPRrfqb0+ugzc1A8nUXQm7vhlTVAG+PE2AKOJ0WYooZiaWFkKubJk2gPdIWR6nWomphAwDTinnQl+QH/uY4Dtq8TIhpF6Lzte2+aT3PsX9yFMlfvGQUXgmJR7EOKGyflKmOxRt0MCwpAd8nCOVEEdqp6WByMqQaS2C5VGeB1+aEYB691txIBV3b2qwDtlR72tTXPk4jgIvzQJuoSTWW8H3YAcCrQKqx0M0JURlWZNLS0oKFCxcOWGbBggVoaRk45yiZHLQF6gE0Un0zPGFGuDOvAteJ2n7bZwX+rUlLhGnZbCSuWozE1Uthvmg+xMKsSdVtRLZ0wLarDFKtBV6rzdfaWGuBbVdZ0GPtcPq/x0KiMewMnbxRB/1s9Q+Hq6JmWHUn40OkAcVoYIoC98k61TJt4RRVoN2XpiAbEPt83xlCTpIVC/2vbZ6OHshN7SHLMsbgOq6eybH/9iT+qQagjqAMmVyGFaEkJCSgsbFxwDKNjY1ISAjfd41MHoaFM1WTzYABtg/3B12QFLeEnu0HoNidquWmFXNjUc1xIZIWR8U5cF5cAEE5gTW5mQP2G9XmZqgP1eMMU5JMBLEMKJjzXKaRPkKm7Tz31IoXhaCMOP5zPtY5n3XFeRDTElXLbB8dhtxqVS1jHi/snx6Fp99y0wV0bRtvmDT4XBGRlCGTy7C6kSxbtgxbtmzBTTfdhBUrVgSt/+yzz7BlyxasXr16xBUk4x+v18K0tBS2XWWBZZ7OHnS+vh3avCyIaUnw9jgg1VrAZHUfSE1uBrQF2ZDrmmNd7bg0ao8weXVgPdjU9n1nqQMwqZ4kTEYxDSj44Js8pgSf47xeC8Xuy3zF+n0H/DeKsc75zHEcTBfOQ9e/PwksU+xOdG3aCc2UtHN5tiW4q88GTeglJBphmD89pvUlhIyNYQXb9913Hz766CN885vfxCWXXIKlS5ciLS0N7e3t2Lt3Lz7++GPo9Xrce++9o11fMk4lXr0cropadZ9FBkh1zZDCBNKcVkTKVy6nfNl9jFaLY/+gRKq1gC0rDUrHGFjfr8uAkGQa9BiERILTa9VPvgBfus8idR9sPskMxemG4pbg7eyB0Kd1mzfpxyzns3nlAjjLqoK+I3JTe9guJeCAlK+sCvt9I4RMLMMKtmfMmIHnn38ejzzyCHbs2IEdO3aA4ziwcwOo8vPz8eijj2LmzLHLf0riC2/SI/3b16Ptz5vC/wD1wWkEpH3jWmhzQ0/PPlmNVoujYcEMdG/e3buNW4bjUCVMS0uDysptVrhO1asGchnOo+82GR0cx8G4cKZ60GONBZrsNNUATF4UIGamwH2qQTXLKScK0BZOGbOcz5woIP3u69C2LsJrG88h9bYroZ9dGP3KkahgigKpoRXuU/XwdvSAebzgtSLEjBToZuVTX3wSZNiT2ixatAhbtmzBgQMHcOLECdhsNpjNZpSWlmLx4sXUGkmCiGlJyHzoy2j+fxvhOlYd1GUEAMAB2vxsmC9eoMqMQUaXJjMFuuI8uCvrA8ucR07D2+2AYd40iBnJYA4XXKca4DxyWpWJBDwH0/lzxqDWZKIyXTgX1n/tDPzNPL6ZGXXTpkKckgZOFODt6IHU0Awme8EZdeDOTRxjWlKChFWLxzTVGm/SI/3em6B4XoPz8KmgLiN+mqnpMF00H8bFJTGuIRkt7hoLbNsPwttv3Iti93WPdFXWQTM1HbrpOUFjXcjkNaJ51TmOw5IlS7BkyZLRqg+Z4HiTHgkXL4Bh/nS4qxrhabWCSTI4QYSQYoZuZi4EkwG8kXKURlvilUvReqpe1UooVZ+FVH12wO1MK+ZCSKRuJGT0aPOyoC3Igruqd+A9k2W4TtYCJ/tkzuEFiKkJ4HVGaDJTwGkEpNx2ZVzkNOa1GpiXz4ahtBBSdRNkSzsUtwROEMAnGKGbkQsx2RyVaxun1QChZt7sX4aMiO2To+h+d3fohqI+PC1WtD71BtK+cS30YbI8kcllWMH2T3/6U9xwww0UZJNh47UaGEoLgeBeCyRGdNNzkHzjxbC++fHQtrl+ZRRrRSarhNVL4O3ohrejJ7INeA6pt18FTWZKdCs2RJwoQDczF7qZuTE7ppBkhjdMOtW+ZcjwOY9Vw/rGDlXjxECYW0b7395B5n99GZosmgBsshtWSoHXX38dd9xxBy6//HL8v//3/1BVVTXa9SKExIB55QKk3LIKnDj4pUA3Iwfpa75Ag7pIVPA6DRKvXQGhXyq9UDidBunf/Dxl8zhHW5gdSI0Y0hgNHp0omKKg662PgwJtMTURpvPnIOHKZTAuKgZvUg/qZW4Z3e98FsOakng1rJbtDRs2YNOmTdi6dSuee+45rFu3DnPmzMH111+Pa6+9FqmpdBdHyHhhWj4bns6r4Dx8Cq6KWii23hzanCBAO30q9KWF0BVkgdOMqOcZIQPi9VqYls+G3NQOqaEFnuZO9I1wBLMepgvmwjB/Og0w7IM36GCYNz10Dn6BH7PBo8DIZ7uNB+6KOnja1TMbGxbMgHHJrN7xaflZMMyfgZ7tB1Wz8zrLz8Bj7aEnC5PcsH45Fy1ahEWLFuFnP/sZduzYgX//+9/46KOP8Nvf/ha///3vceGFF+L666/HqlWroNONny8UIZOVYNTBuHAmDAtmQHG4wNwyOIEHbzJQSzaJOTE1EWJqIhSvF8wlA0wBp9GATzTCMKcIvI76H/enyU6FkGSKq8BWtnQE3QB4rTZIDS0wzJsOTfb4aJiz7z2u+ltMSVQH2udwogDzxQtgfX1H70IGOPZVIGE1dbudzEbUTCWKIlavXo3Vq1fDZvv/7d15fFTlvT/wzzmzL8lkhUAgCQESlrDvi6hIQYVSVKqtC7YVW1vtr9bbqrV7+6r0ttfbKqVXaet+XVFuKYpbRVEw7Pu+JZANQvbZz8w5vz+GhJzMTPbZks/79fKFc+Y5Z74zczLznec8z/exY9OmTdiwYQO2bNmCLVu2wGKxYNeuXb0VK1HCURQFvksNkCovQW5yAQLgOVYatyWiBEGAxmIC2lwOJYoGxeeHt6wa0rkLkF0eKIoCUaeFmJYMw5BMVrnqgGgydLygVZR0ZrVbjc2SED3cUlWt6rahcGjYc1E06KHPHwy58Uq1El+b/an/6bVrwlarFV/96lexZMkSvPjii1i9ejUcDkdvHZ4o4UgX6+Buu3y6qIGvrokloohaURQFzv2n4Sg+BF9lLRT5yoqlfrcH/iYHpNJK6AZlwJCXxWpFCaDXVruNA22rj4gdVGPSJFtUybbM5dv7vV5JthVFwdatW7FhwwZ89NFHcLkCYz5nzJjRG4cnSgity2+5jpbAtfs4FH/bZc7VvSEsEUX9naIoqF/3CRxbD0LxtV9STbpQh/oNW5F6yzXRCY66rbdWu40HolGP1p/k/up6IGdg2Pa+S/VB+1P/1qNk+9ChQ9iwYQPeffdd1NTUQFEUjBgxAkuXLsXSpUuRlRVfl8iJIqm5/Jbn3AU4th0MXSIqxPhnloii/qzpo11wbDvU6fayw4WGd7YhaeG0TiUxfWGCXiLqrdVu44E+NwtSxZXVQd3Hz8E0YUTI+Sz+RkfLCqit96f+rVvJ9v/8z//gX//6F86ePQtFUZCRkYG7774bX/nKVzB6NAsnU/+kz8uC9/wFOL84FJRoa6xm6IYOgLEoH3KDHe7j5yE7rlT9aC4Rlf6txVGOmih2/I1ONL2/I2i7Js0GXVYaBL0Wvtom+MqrVVeJ/A0O2D/di+RF7V897SsT9Ci2LLPGwvHF4ZbbstONpo93wzpvouoHn7/BjqZ/71Z9/gt6LVcMpe4l208++SRMJhOWLFmCpUuXYs6cORDFbpXs7hf64upeffE59ZRoMkA0GuC3u1Xb9XmDYSwcCn3eIGhTkwCg3RJR2pSkqMZNFCuO7YehtBnXa5o0Etq0K7W2dZmpkIcNgmv3cfhbfeY4vjiMpAVTIWhCV8vpSxP0KLb0QwcGOlNKqlq2ec9dQN1rH0GfOwiaZDN8NQ3wnr8IKAoE7ZXUyjx1FESTAYrSydVwqEsS5cpVtzLk3//+99i6dSv++Mc/4qqrrmKi3YHO1NdMtBqcffE59QbPmXLostIgWs2ATgttqg2WOUUwjh3WkmgDV0pECbpWicLlElFE/YWzTUk1w4gh0A1KD2on6nUwjB2m2uZvcMBz/HzYY3d2gl5v6kwHQ3/rhOgrUpZfC6FNyUnF54fndBmce0/Ae+4C0Cah1qYnI/mGmdEMs1+Rqmph//wAvKVV8NfbA1etSqtg//xAUAWZWOtWz/ayZct6OYy+TZ+XBW/ZxfAf/G1W9wr1S032xNfYtq4+p/5CqqqFoNVAm5oExe2FeWohDENDT6RhiSiKZ/4mJ9zHSiGVX4Ls9kIQAdFigj5/MAwjhkBEz3qOFJ8fvkvqhUKMhUODFg9ppk22QJNsVm2TqmrDLm4Tiwl6XDa979JnZyBj5RJU/vq5DifyAoFEO+M7X4HGyjKqkZBoV67YJR0Fzat7hVxOt83qXuF+qbkOnYGvrinKkYfXlefUn3SnRFRrLBFFsSY73Gj8YCfq3vg3XPtPwXepHrLdCX+jE1JlDRxbD6LulQ/h2H28R5fGQ02O6+jvpe1niiKF/3uJxQQ9LpvetxlGDEHKV6+FftggIEydbUGnhbEoD5kP3gptZkp0A+xHYnHlqic61bM9atQoiKKId955B8OGDcOoUcErJ4UiCAKOHDnSYbv+oDOre7X7S01W4C2tgmg2xs3qafG4YlmssUQUJTJ/gx3Va9bDc6o8dDWdyxTJB+f2o6h96X2k3bkQQjeGEgpGPSBA9Ti+Sw3hHxMK/I1OiKYrfyOCMb4+Z+J52fRYkL0SfBfrIDsC81hEixHaAakQE3gojTbFiuTrpkJ2uOA5WQZfnR2K3wdBr4M2MwWG/GxoU63s0Y6wRCst2alke9q0aQAAk8mkuk1d09HqXp35pea7WAf90AHtPo6/yQnfpXooXh8ErQaa1CTVhKPeFE8rlsUDloiiRKVIPlz627/gq67v9D6uvSfRkGxByrKruvx4gihCP2RAYFLZZe7DZ6EP8+PUV10P2eVRJdv63PA/ZGOFnRABvtpGuI+dgyL54K+zt1y1E416mCaOgC4jJbYB9pBoMcE0cWSsw+i3mq9KyS4PpIt1UNxeCKIIMdkM7YBUCIIQV6UlO5Vsv/TSS+3ept7RmV9hstMdcruiKJAu1MJzqgzesuqg+7UZKTCOzoVpwvAex0nhsUQUJSpH8WFVdRwgMJHXUDAUusEZUCQ/vGcrAhPBWrFv2QfrnHHdumRumTVWlWxLFZcgO9zQ5w1SLQDlq7fDfbhEta9uUHrc/jjt750QsssD95ESuI6dg6/iUtAYZ9f+UzDPHAvrzLExipASnXSxDs6dx+EtqQiqaKRJMsM4Og+myQUxii5Yry3XTj3XqTGGUvDEDMUvw75lH5x7TkJjM4fYKzBcwf5ZPTxnK2CeMgqapNDtqGd6VCJq2qh+1fNF8UNRFNi3HlRt0ySbkbx4NjSWK5fDjSOHwHv+Iho/aFUbWwHsWw92q3fbNLkQDRu3QW5V0s97thLecxegHZgKQaeFv64J/vrL81XEK1eIrFeN79RwRoq+xo92oWnz3rBj6hVJgnPnUUhnK6HPHQhTUX6UI6REpSgKmj7cifp1n4adqOpvcsKx4wjcJ8/DNHZYyApH0dYnJkgqioIPPvgAd911F+bOnYsJEyZg0aJF+MUvfoHz58OXhuoLFFlG08e74Tld3qn2vqo6XHr6/yC7vRGOrP/qdomo61kiimJDOncBvgt1qm2W2eNViXYz/dABQRVAnDuPdmuypGjQIeXW+UHbFUmCVHYR3rMVVxLtVnTZGTDPGNPlx6PIcx0+i7rX/t3+kEgA8Pmh+PyoefYduA6fjU5wvYDlHWOr6f0daNy0vVNt5SYXqte83aWhcZHSqZ7tn/zkJ906uCAIePzxx7u1b1f853/+J5577jlkZmbiuuuug9VqxbFjx/DGG29g48aNeO2111BQED+XE3qT+/BZeEuDZ9wKOm2gDFWTA0qbsoFSRQ3q129B2tcXRCvMfoUloijRSJU1qtui1QRddkbY9sbCHFWNa9npgdzkDKqu01q4xSfME0ZAufVaXHxqXadi1WWlI/n6Gd2alEmR5W90oval90Mm2uLlH27NK+cq8uUfZwpQ+9L7yPrpioS44sryjrHjPlmGxhArzgqiCE1qoESybHep7pMdbtS8sAkD/uNrMb0S1qlke/369SG3C4IQsjejeXs0ku3q6mq88MILyM7OxoYNG2C1XjnJn3/+eaxatQrPPfccVq1aFdE4YkGRZbgOqXsEBJ0W5umjYRg5BKJWC0WW4T13AY4vDquWB3fuOgbb4lntfjlS9zWXiLJ/tj8wpCTU34lOC+PYQIkoJtoUS22Hp2mSzO1+MYnJwUlR2x/1rXW0bLplVhFspVVwfH4Q/noHIPuDjiHodTDkD0byDTMTIinrjxzFh4LOA+3AdBhHZEM0GwEE5rG4T5XBV32l8ozikeAoPozkL8V/8QWuMRE79s171BuEwGrMpqL8llU6fRfr4PjisKrSl1R+CZ6TZTAWDI1uwK10Ktn+97//rbotyzJ+97vfYf/+/VixYgWmTp2K9PR01NTUYOfOnXjppZcwceJEPPbYYxEJurXy8nLIsozJkyerEm0AuOaaa7Bq1SrU1vbNhUK85y+qEmgASLp2smo2vyCKMOQNgjYtGfVvfXKloawEPtwWTo9StP0PS0RRomhbQs9X2wjFL0MIUzPaXx1cok8IU7ays4tP6AdnQLtkNvTDBsF18AxklwdQ5MDfS2oStAPSIJoNEDQiew7jkCLLcGw7pNqmzUyFaXw+BFz54SaajTCNHw73obOqCeKObYeQdN2UuL9iwfKOseGraYD7WKlqm3nKKJhbVYQRBAG6gWlIvnEm6t/eAqXVcFnH1gPxn2xnZ2erbq9duxYHDhzAP//5TwwYcKUMXX5+PqZNm4ZbbrkFy5Ytw3vvvYd77723dyNuIzc3FzqdDnv27IHdblcl3J9++ikAYObMvjkWVmo1ix8ANKmWsGWzNMkWGIYPgedMRcs299FSJttRwBJRFO/altBTPBI8p8pgLMwJ2d51pER1W5ueDDHMj8auLj5hGDEk0NPeSz2HiqIEygY2OKDIMgSDDrqBaUyIeplUUQN/g0O1TZ8/WJVoNxMgwDBiCDwny1q2+evtkCproW9n+FK8YHnH6HMfPxdUvctYNCxkW1Gvg2lcPpw7j13Z/2hpy4iLWOhWNZJ169bhhhtuUCXarQ0cOBA33HAD3nzzzYgn26mpqfjhD3+IP/zhD7jxxhsxf/58WCwWnDhxAl988QVuu+023HnnnWH378kKaL1NUcLHoygKFCiBNSAut5FdHlUbbUZKu89Hm50Bz5lyNJ+x/kZHTJ5/e88zXJtw+zRvUxSlW8eNpEjEE6nnGG+vXaS1Pm8I0GamQD98MLwtE60DV74EswH6IVc+5xVFgWvfSXjPVlyuHR94/cyzxrbc35avwd7h6xxoE9hf0GuhyxkI6dwFKH71cBJBFGAsyodg1Hd4TL/XB/eB03AdK4XcqE4CIQjQ5w2CeWphF//+eN6E429y4Eo2pAA6LTTJZihhVkcSbZbLE8mv3O9vdEAZHPvKEZ0hGPUwjAr+Mdred3h791P7/I1OtD5XtFnpEDSasK+nLjsT2Hm0ZR9F8kF2eyOycFxnEvhuJdtVVVXQ69sP2GAwoKoqOktl3nPPPRgwYAB+8Ytf4NVXX23ZPmnSJCxduhQ6XfiZwQ0NDXFz8rvtdihtEuhmXrcbikeCIMvwOwJfHB6vF75WYy1ltwcOhyPk/gAg2R2QJAmyJ3BpReP1or6+vveeQCe19zybCbIPvlaxdbSP0+mE5HLB187zD3XcSOrO84zFMSN53HjX0BB+xcL+RpkwDJ4jgTkgfkmC4pNRu3EbNOnJ0AxKAyQ/fOcvtpTpExQZsscLQaeBVDA47GeJu6auU+cWgCvt9CLknEzI1fVQXIG1BQSTEdqcAXAYRaCD89B/qQH1//s+/LXBlUyaSSfOwXmmDC7IMH1pSpd6vHjeBJPsdnguf7f4JQl+yQu3yx1+WXPZD6/LDb/nSk9wk8MOdx/7jGmL5073uJzOlvMLAODzQWzn+97f2KjKdwCgobERgrv3K16npqZ22KZbj5qVlYWPPvoIDz74IAyG4EsmLpcLH374IbKyojNJ4K9//Sv++te/4oEHHsCyZcuQnJyMo0eP4ve//z1WrFiBP//5z1i4cGHIfW02W1Ri7Ay71QpZDP2WiEYjFEGEaNDDaLk8qTHNBnfllfHoSk0jzObwE5uaLtRDp9NBZwj8UDIMykBKSkqvPofOaO95NhPNRlhbxRZuH0VR4HQ6YTab4TWZYLC0P+Gz7XEjqTvPMxbHjORx45WiKGhoaIDNZmOt5suUOROhOVUF94FTkHQ6KMLlXuVGJ5RGJwBAA0CjC5wnglYDnUGPlFvnwzJkUNjjdvbcAqBuZwGQagtq19E56KttRPVLH0Hj8kLUtf+4glYDpfgoRKMJtqVz2m0L8Lxpj5SjwHv5u0XS6SBAhNbhgTYj9Hesv8EJRRRbvo8AIGXoIOj6yGdMWzx3ekY3ZCDkVueKVFULk6gNO2zHebxMle8IZiNSMtMTaxjJ8uXL8d///d/4+te/ju9973uYMmUKUlNTUVdXh927d2PNmjUoLy/HQw891NvxBvniiy/w5JNP4hvf+Abuu+++lu1TpkzBM888g+uuuw6rVq0Km2zH00kvCOHjCWwXWv0/YMgfDHer+qSy0wvP4bMwjQteJdJbXg3vuQuXF1EJ7G+aODImz7+95xmuTUf7CILQreNGUiTiEQRAtrvgPlYKqfwSZLcHgiBAtJqgzx8Mw4gh3XqO8fbaRUvgvOlbz6m7BEFA+p0LUfO8H1JlZyaVC7B9eS6ss4o6OG7nzq3mGHpyLEVRUPfyB5CbXECbscKixQTRqIO/3t5qxbnA56r9k70wjhwSVD88fBw8b9rSDUiFbnAGpIoaKAoge7xwbD0ITboNkBUIGhFikhn67EyINgs8J85BNBrQ/D7pBqdDd3mZ7b6M5073mMfmo0GnhSJd7gRQAvX9rfMmBr2e/kZHYAIuruRN5okjIcZw8m23ku2VK1eipKQEb7/9Nr7//e8DAERRhCwHPsAURcHNN9+MlStX9l6kYTRPgpwxY0bQfWlpaSgsLMTevXtRW1uLtLS0iMcTTdoBqdCmJ8NX09iyzbH9CHx1TTCOGQZtejLkJifcx8/BfeiMal/BoIN5yqhoh0w9IDvcaPxgJ9zHS9F2GKS/yQmpsgbOHUdhnj4alqsm8AOdukzQaZF+z2LILjecu08EVTtqps1MhXVOEZLmT45yhO3zllSpJlsCgDYtGZZZYwNjPAUBstsL95GzcO49qWrXtHlPp5NtCiYIAixzxqHm2XfhPRNYjEhRFPgbndAkmyEYDfDbXfCerYDs8kCTpJ7Qb53LFUEpPNFihGlSAZw7jrZs85wsg+L0wDh+OHSD0gMTu0+Xw7n3BBTJp1qh2TJ7XCzCbtGtZFsURTz++ONYtmwZ1q9fj+PHj7dUAhk1ahSWLVuG6dOjU+VCurwcbLjyfs3bOxpjnogEQYBpUgGaPtql2u45cR6eE+2vnJl07SSIBq5ylSj8DXZUr1kPz6nyoES7NUXywbn9KGpfeh9pdy6M+zJaFH8EUYR54kgYRgyB99wFSOXVgRVnBQEaS+AKii4zBaI5/qouOD4/oLotmo1IXjwLYqvLz6JRD/PkQgg6LZy7T7Rs95wqh1RVC11W3+qUiSbRYoT3dDmUVuNkFZ8PvtrGoLay2wsx2QzdoHRo05NhmlIYzVApASVdOxmuPceh+K5UKvKWV8NbXt3ufsaiYTGvctOjkeLTp0+PWlIdzuTJk/Hyyy/j+eefx6JFi5CUlNRy3/r161FaWoqxY8cG1eDuKwx5g+CfVADn3hMdN77MND4fSQmweAAFKJIPl9Zu6NKSs669J9GQbEHKsqs61V7Q6wBnB5PYuARxv9Jco9+QF348drxpW4fXVJSvSrRbM44ZBtfhkqD9mWx3j7esGnX/+yEMBUPh2x5+YmprvopLkIcNRvrKL0Pk5wt1QJeVhtQ7F6H2hU2d32dwOtJu/1IEo+qc3p+WGWXXX389XnvtNezYsQMLFy7E/PnzkZycjOPHj2Pr1q3Q6/VRWVwnlkyTCyDotUE93KEYi4YhbcX17PFMII7iw5Aq1MtpC1oNDAVDoRucAUXyw3u2At5zF1Rt7Fv2wTpnHLSZKR0+BpcgpkSnyHJLpZRm2uzMsO0FjQh9doZqH7nJGbH4+rqGf22FIvkhmo3QDxsE97FSCLIC0WKCoNMAPj/8DjcUX6vVSvVaaDNs0A7suJoDEQCYJ4yAcM8SXPjDK+pzKQT90Exk3n9zXNQ+T/hkW6PR4B//+AdeeOEFbNq0Ce+88w4kSUJ6ejqWLFmC73znOygoKIh1mBElCAJM44ZDAaB4JbiPnVOV2hJ0WhhGZMM4Og/6IZkQNJrYBUtdoigK7J8fVG3TJFuQvHgWNJYri4gYRw6B9/xFNH6wo9XOgH3bIaR8ZW6Hj8MliKlP6uDLWPG1WRZe5Jjh7pAu1qmGLooGPcxTCqBJSYZUVg1/fWAYiWA1BRa+8UoQk80w5GRBtrvgPVMBw/DscIcnUjGNHYa0OxfCtf8U3EdLWxYWAtBSQ980Jhf6/OyWSkexlvDJNhAYj33vvfdGfAGdeKcxGWAYNxzmyYWQ7S7IXh9EnQai1cSe7ATlLa2C72Kdaptl9jhVot1MP3QAjKPzVF96zh1HYFs6p8OJR1yCmBKdIIrQ2CyqVQy9JVXQDQw9LET2eCGVV0ObeaVXVZOaFLItta/1pDUgMAHfOCYP8PmhH5QOWfJB8QbmVwk6LZw7jkD2+IDL30uOHUeZbFOXaJLMMBXlwzh2GBSXJzCvRBQhWowtJT/jaV4aM7AEEVitrYM2zfVvBQGaJDN06cnQJFuYaCcwX5V64q+YZIKunYkexjYrmslOz+UyaB3TZaXBOnc89LlZ0KRYoUmxQp+bBevc8RzHShHXmTkBHbUxTRihuu0+WgJfiOFRiqLAufOYaqIVRAGmouCyqdQxqUo9zM0wcggEzZXvHVGnhcZigsZigqjXBVb3a8XXZn+ijjQPaxQEAaLZCG1aMrQpVlVt/Xga+tgnerb7A9FihGxvP2kSLfFxuYR6T3NvUDNNkqXdXmoxydzhMdojmgwwjs7tfIBEvaQ35g1Y5oyDfcv+ltuKz4+GDZ8HqquMHALBoIPvQi1cB87Ae/6CqjSYafwIaJKD/36oY4pXPVxHY7NA8YT/3AnU175CbqctUSiJNvSRyXaC0A5Ihe9SO8u8akRoB6aFP/EoIQltvpR8NQ1Q/LKq16g1f4hzRDD2vbKX1Pf0xpenbkAqTJNGwtWqhrbileDYcQSOHUfC7idoRCRdN6VbcRMgmtSfMb7qhnZ/uPga1ctsc4gadVWiDX3k+IIEIep10OdmhU6ymk+sOBqfRL1DnztQdVvxSPCcKgvb3nXkrOq2Nj2ZVzwoLiiyDG/FJbhPlsF9/By85y5AbnXVpfnLE+19xnXiyzP1tuugzxnQ+cAEIPVr10E/JHzlEmqfPlddHtJ7ugxyq1rbrck+P3wVl9rsHz89kJQ4EmnoI3u2E4g2NQnaDBv0Qwa0zL7V2KzQ52VBNBngOXEuxhFSb9MNTIMhfzA8ZypatjmKj0C0mFTJgSLLcO0/Be/ZSvWqWbOKuCpbgpBdHnhLqkL+bScyv9MD1+7jcB0rVVVJAi5XShqeDfO0wGq2uqw0aGyWHr0OokGHjO/eBOlCbWARqHYIOg3S7r4B5jZjvalrzNNGofHdL6Bc7mFU/Aqcu4/DNGaY6n2TvRJcB88ESraJV+YhWWaNjXrM1DckytBHJtsJRjToEuLEigVB1/HpnIgLs1jmTVAl24okofG9YmgzU67U2S6phOx0q/YT9FqYZ4yJdrjUDVJVbdDlUH+9Hd6yizCNGx53vTSd5TlTgbpXPwq77Lsi+eA+VgrP6XLoBqbBMnNsr3x5ikY9km+YBW9JJVyHz0IqqYIiX3ltNckWGEfnwjRhJBPtXqBJMsM0cSScu4+3bJPrHbB/fgDazBRokszwO9zwV9dDkdXlFg0FQ8NWjCHqK5hsU5/RmeES8TQ7ubNM44fDND4frgNnVNt91fXtriqZctM8aKzBJQIpvsguT+hxhwDgl+E6eBoamyXherg9JZW49PT/tTtRroUC1L3+MQDAMrN3ejkFIXBlSDcwDbLPB8XhhuKXIRp0EMzGQBUDY+L9+I5XtqVz4Dld3maSqwJfdR181XUh9xEMOqTccnV0AiSKIY7Zpj5DOzAt9HjPZnE2O7mzBEFA2p2LAnVrO8m2dE6vJS0UWd6SqvYnNvvlQJsEInsl1D6/CYrkD7pPk2yBNjUZCDG6qe7NzZAu1Abf0UOiVguNzQptWnJgRUMOrep1mmQLMr7zFWhSLndodFByVtBrkHLzPOgGcPVI6vuYbFPU9UY93VBEg65XJljFI0GnRfo9i2GZNRZiiAVtmmkHpMK2ZBaSrp0cxeioJ1Srn/WgTTxx7TmhWlwGAHSDM2G76Wqk3jofKbdcjdSvLYBxVJvhIrKiKt1HiUWXlYYBP7wNxtG5EAyhL5wLggjtwHRYZ4+HsSAnZBuivobDSCjqeqOebji9McEqXgmiCPOkQL1g77kLkMqrA6tmCQI0VhP0+YOhy0iBaE7s59nfdKYOeldqpccDx7ZDqtvazBQkXz9dtcCWxmKCde54AFBNZHTuPg7bl+dAZMnKhKRJNiPp2kkwjs5F0+a98Nc3QfH7IWg00FjN0A5KvzyUx5CQVxqJuoPJNkVdpIvRJ8rs5O4SRBGGvEEw5A3quDFRlPmbnPCev6jaZppUEHYlW9OkkfCcvpJsKx4J3jMVXRo2RfFHk2yBdd4EeEtDDJPSiDAMG5zwHSBEncVkm6Iu0YrRE1HnhVrpVp8dvoa1xmIKjOFuxd/k7PW4KPq0qUkQzUb4Lta1VEsSzUZoB6RCm5p4k9WJuovJNsVEV4d7CHod4PQEbQ9qQ0Q91qO/txCTDwPDCMJPEVL8fghodbwOJtdR4hANOuiHdmGRIaI+iMk2xUxXhntEcpw3Ean15O9Nk2INVBpRrmzznq2EsTD0ZDipphH+BgdE85XSnS0VLYiI+gB2H1BC0OdltVvWT9BqONmGqJd09PfW3rwK0agPqjLi3HsCsiu4p1yRZTi3H1Hvn2SGIZ/zEYio72CyTQmheZx3uLJ+xqJ8jvMm6iUd/b11NK/CMmec6rZsd6H+n5/BffwcZK8ExS/DU1qFhn9thVRRrWprnV0EQaMBEVFfwWEklDBCjfMWbRb4UkwJu5w1UbzqSRlN4+hc6IcOUFUlke0u2D/bD3wWvo62aDYEJerdxXkescPXnkiNyTYllLbjvBVFgbu+PnYBEfVh3S2jKYgi0r51I6qfXAepsqZz++g0SL9nCTRJ5i4/Xiic5xE7fO2J1DiMhIiIep02JQmZP1gObWZKh21FixGZ998MQ/7gXnv8now7p57ha0+kxmSbiIgiQpuShJRbrkbSwunQDxkQqFLS+v6MFFivmoDUO74EfW7vJl89HXdO3cfXnkiNw0iIiChiBFGAIWcgDDkDIXulwOImsgLRZGhJuERtZCZE9mTcOfUMX3uiK5hsExFRVIh6HcQoT4zr7rhz6jm+9kQBHEZCRERERBQhTLaJiIiIiCKEw0iIiIiIekh2eThGnUJisk0UIVzYgYiof5CqauE6eBrwyy3b/PV2eMsuwjRuOLQDU2MYHcUah5EQRUhnFm3gwg5ERIlNdnmCEu0Wfhmug6chu9rveKG+jck2UYRwYQcior7PW1IVOtFu5pfhLa2KXkAUdziMhChCmhd2CNnjwYUdqBV/owPuY6WQyi9BdnshiAJEqwmG4dnQD8+GCJ4nRPGqeYx2+20cQJYtCtFQPGKyTRRBXNiB2uNvcqLx/R1wnzgHKMH3SZU1ELYfhWXGaFiumgBBEEIfiIhiRvFKHbfxdNyG+i4m20QRxoUdKBRffRMurXkbntMVQYl2a4okwVF8BNrMD5H69QUQRI7+IyJKJEy2iYiiTPZKqHlmA3yXGju9j3PXcWiSLbB9eU4EIyOi7pKdHkgXa6F4JAiiCNFmgTYzhVekiMk2EVG0ObYdglRVq9om6LQwjBwK3eB0KJIf3jMV8J6/oGrTtHkPLHPGQZuWHM1wiagd3spaOHcdg7ekElDUl6k0NiuMo3NhnDgyRtFRPGCyTUQURYqiwLH1oGqbxmaF7caZEC2mlm3GkUPgOXcBTR/ubLVzIFG3LZkdrXB7jPXmqa9SFAWNG7ehYf0WKD5fyDb+BjscxYfhPlkGS24GkJIS3SApLnDwHxFRFHnPVMB3qUG1zTp3nCrRbmbIGQjjqBzVNseOo1CUdgZ5xxnWm6e+qmHDVjR9vKdTbeUGO5r+/g58dU0RjoriEZNtIqIoki6oh49obGZos9LDtjcUqpNtuckJpYOe4njCevPUF7kOn4X9k71B2wWNCG16MkSzMeg+ucmFupfeT6gfy9Q7OIyEiCiKFK/6crOYZGl3ApUmyRy0TfZKEC3BX+bxiPXmO0d2eVgiNIHY2/ZoCwLMkwtgHJMH0aCHoijwVdXAse0wfHVXJkJ7SyrhLa2CIW9QlCOmWGKyTUQURaJRr7rtq2mAIsthS/q1HXIS6hjxjvXm2ydV1Qb9GPHX2+EtuwjTuOHQZaXFMDpqS6qsgedMhWqbZcYYmIryW24LggDdoAwkL5mFunWfANKVH9mOrQeZbPczTLaJEggnmyU+fa56yITi8sJ7ugKGkUNCtncfKVHd1g5ITcgElfXmQ5NdntC9/gDgl+E6eBoamyUh3/O+yn20VHVbMOlhHJ0Xsq1o0MM8dhic+06G3Z/6Po7ZJkognGyW+HSD0oPGKNu/OARvebVqmyLLcOw+Bm9plWq7ZdbYiMdI0eMtqQqdaDfzy4E2FDdku1N1W5+dCaGdeQnaIQPU+zvcUOR23nPqc9izTZRA9HlZ8JZdDP/lzMlmCcF61QTUtkqgFK+Exk3F0GamQDc4A4rkg7ekEnKbqxiCQQfL9NHRDpciqHlYTU/bUBS1mWOh+P3tt297vxB8DOrb2LPdh3Rm+ACHGCS25slmIas7cLJZwjBNGglj0bCg7b7qerj2n4L7SElQog0AKbdcE7LKASUuxSv1ShuKHk1qkuq2dL4asif8e+Q5W6ne32blqpL9DHu244zf7oL7WCmksmrIbg8EQYBoMUGfPwiGEUMgmsMnUhqbFf769ntAOMQg8XGyWeITBAFpdy5C7fPvQqqs6cQOgG3ZVbBMGxX54IioXaZx+ahfvwWQAyX8FJ8fzt3HYJlVFJRE++qa4DlaGujNbt6fq0n2O0y244TscKPxg51wHy8F2pTg9Dc5IVXVwLnzGMxTR8Eyd3zIygUcYtB/cLJZ4hMNOqSvXALZ5YVzz3HITnfIdtqBqbDOHoekqydGN0AiCkljs8I0Lh+u/adbtrmPlEC2u2Aalw9tVjoUlxeek+fh3HcKit8PQatpaWudXRSLsCmGmGzHAX+jE9Vr3obnVHlQot2aIvng3HkMtS++h7QV1wcl3KxnS5RYBI0G5skjYSgYAu+5C4ErWh4vIAjQWE3Q52dDl2Fr94oWEUVf0oKpcB0809K7DQDecxfgPXeh3f1MUwqhzUyJcHQUb5hsx5ji86Pm7/+C72Jdp/dx7T+NhvWfIeWWq4Pu4xADosQjiCIMeYNYe5coQeiHDEDq1xag7pUPO72PdmgmUpZfG8GoKF4x2Y4xx44j8J6/qNomaDQwFAyBLjsTkPzwXF5xqjX75wdgmTMu5GIHHGJAREQUWZZpoyDotbj4X69C8fnabavPy4LpmzdANLBIQX/EZDuGFEWBY+tB1TbRaoZt8SzVEs2GkUPgLa9G0/s7VG0dWw+G7N0mIiKiyDNPGIG0u66Ha99JuI+Vwt/oaLlP0IjQ5w2GcUwudLlZ8DPR7reYbMeQVHYRUoW6EoF1dpEq0W6mz86Eceww1cpTzl3HYLvpqrDLPBMREVFkiUYdTOOHwzguH4rTDdkjQdCIEC2mlomRLPXXvzHZjiGpSj1OWzQboBs6IExrwFCYo0q2ZbcX/gYHtG1qfhIRUeKQJR88p8rgPV0B2eGCoigQjQbosjNgHJXLCbIJQhAECBYTRIsp1qFQnGGyHUOKpC6CLyab2/31q0kO7vFW2imkT0RE8UtRFDh2H4dzx1EoknrMr2x3wXepHq4Dp2AszIVpciE0ViZxRImI4w9iSDSqeyv8NY1QfOGXffVV1wcfw6Tv7bCIiCjCFFlG7Uvvw7k9ONFWNwQ8p8pR/eSbXLY9TnVq9WaO1+7XmGzHkD53oOq2IvnhPnk+bHvX4bOq25oUK8RkS0RiIyKiyGnYsBWuvSc73d53qQGX1m5oPzGnmOjMyswaG7+r+zMm2zGkTbfBUDBUtc254yg8bYriK7IM576T8J6pUG23zA5eGpaIiOKbr7oe9i37grbrcwYi6epJSFowFcYxeapVBwFAqqiBo/hwlKKkztLnZQGadtIpjQh9Lldv7s84ZjvGrFeNh+fEld5sRfKh6YMdcKXboBuSCUXyw1taBdnhUu0naEVYZoyNdrhERNRD9q0H1asFCwKSF06HvtUEeUPeIJgmjEDjO9sgO72qfS1zx7OjJY50evVmjyv0AajP61PJ9ocffohXXnkFR44cgcvlQkZGBiZOnIgf//jHGDQoPldmM44dBtOkkUGXE301DfDVNITdz7ZsXsgJk0REFL8URYFz51HVNuOYPFWi3UxjMcEyezyaPtrVss13oQ7e0iquNhpnOlq9WVGUDo5AfVmfSLYVRcEvf/lLvP7668jJycGNN94Ii8WCixcvYufOnSgvL4/bZFsQBKR9fQFqJB+kypqOdwBgWzwL1jnjIhwZERH1NrnJCdnpUW0zFuaEba/LzoDYpgqJr6qWyXYc4urNFE6fSLZfeuklvP7667jjjjvw05/+FBqNepybr4NlVGNN0GmR/s0bITs9cO4+Btke+lKTNjMV1tlFSFowNcoREhFRb1C8wd9HYoiFzJoJghC00JniZclXokSS8Mm22+3GmjVrMHToUDz22GNBiTYAaLXx/zQFUYR5wnAYhg+Gt6waUtlFyC4PIAjQWEzQ5w+GLjOFixsQESUwwRhcrtV/qQHioPSQ7RW/DF9tI7Rpya2Owe8BokQS/1loB7Zu3Yr6+nrcdNNNkGUZH3zwAUpKSpCUlITZs2cjNzexLukIoghDzkAYcgZ23JiIiBKKaDFCm5EM36XGlm2uI2ehC5Nse06VBS1eps9jZQuiRJLwyfahQ4cAABqNBkuXLsXZs1dqUYuiiG984xt45JFHwu4fT5MWFKXjeDrTpj9pfi34mlBXxNN5w7/7xNFb54155lg0btzWfFR4z1bCsfs4TJNGqqqMeMsuXi71p6C5fIl+eDa0mSk8HxJMPH3mUO/qTGWghE+2a2oCkwqfe+45jBkzBm+++SaGDx+Oo0eP4uc//zmeffZZDB06FLfffnvI/RsaGuLm5Hfb7VBcnnbbCLIPvvr66ASUQBoawlduIQonHs4b/t0nnp6eN3JhNrwb/FAkP/ySBMUnQ9pxBPZDp6EdOgDQaeCvrIW/JtD7LWhFyJ5A+T/dhGGo57mQsOLhM4d6V2pqaodtEj7Zbk6UdTod1qxZg4EDA8Mvpk6diqeeegpLly7Fc889FzbZttlsUYu1I3arFbLY/lsimo2wpqREJ6AEoCgKGhoaYLPZWHeWOi2ezhv+3SeOXjtvUgDD1xei/vV/Q9LpoAj+wHbJD+VMJYDAinOiLnBeCFoNdAY9jOOGI23OxJifs9R18fSZQ9GX8Mm21RpYJrWoqKgl0W42cuRIDB06FKWlpWhsbERycnLQ/vF00gtCx/F0pk1/JAgCXxfqsng4b0SDrsOebdGgi3mcdEVvnDfWmWOhuDy49MyGzjwijGOGIf2uRRBELvycyOLhM4eiL+H/avPz8wEASUlJIe9v3u52u6MWExFRZ2ls1l5pQ4kn6drJsC2ZBe2A8JehRYsJltljkX7PYgi6hO8fI+qXEv4vd8aMGQCAM2fOBN0nSRLOnTsHs9mMtLS0aIdGRNQhfV4WvGUXg5d5bqYRWX2iD9PnDERKRgqkS/XwnqmA3+4CFAWiyQDd4AzocwZCYzWxR5sogSV8sp2Tk4O5c+fi888/x5tvvomvfvWrLfetXbsWjY2NWLp0aULU2iai/kc0GWAaNxyug6eDE26NCNO44RBNrKvc1+kyUqDLSIl1GEQUAX0iA/3lL3+Jr33ta/jZz36Gjz76CPn5+Thy5AiKi4uRnZ2Nhx9+ONYhEhGFpctKg8ZmgbekCv4GO4DA0BF9XhYTbSKiBNcnku2cnBy89dZbeOqpp/DZZ59h69atyMjIwB133IH7778f6emhFwsgIooXoskA4+jEWoSLiIg61ieSbQAYNGgQVq1aFeswiIiIiIhacMYFEREREVGEMNkmIiIiIooQJttERERERBHCZJuIiChGBL2uV9oQUfxisk1ERBQjXEGUqO9jsk1ERBQj+rwsQNPOVzFXECVKeEy2iYiIYqR5BdGQCTdXECXqE/pMnW0iIqJExBVEifo2JttEREQxxhVEifouDiMhIiIiIooQJttERERERBHCZJuIiIiIKEKYbBMRERERRQiTbSIiIiKiCGGyTUREREQUIUy2iYiIiIgihMk2EREREVGEMNkmIiIiIooQJttERERERBHCZJuIiIiIKEKYbBMRERERRQiTbSIiIiKiCGGyTUREREQUIUy2iYiIiIgihMk2EREREVGEMNkmIiIiIooQJttERERERBHCZJuIiIiIKEKYbBMRERERRQiTbSIiIiKiCGGyTUREREQUIUy2iYiIiIgihMk2EREREVGEMNkmIiIiIooQJttERERERBHCZJuIiIiIKEKYbBMRERERRQiTbSIiIiKiCGGyTUREREQUIUy2iYiIiIgihMk2EREREVGEMNkmIiIiIooQJttERERERBHCZJuIiIiIKEKYbBMRERERRQiTbSIiIiKiCGGyTUREREQUIUy244ig1/VKGyIiIiKKD0y244jGZu2VNkREREQUH5hsxxF9Xhagaect0YiBNkRERESUEJhsxxHRZIBp3PDQCbdGhGnccIgmQ/QDIyIiIqJu0cY6AFLTZaVBY7PAW1IFf4MdQGDoiD4vi4k2ERERUYJhsh2HRJMBxtG5sQ6DiIiIiHqIw0iIiIiIiCKEyTYRERERUYT0yWT7b3/7GwoLC1FYWIh9+/bFOhwiIiIi6qf6XLJ9+vRpPPXUUzCbzbEOhYiIiIj6uT6VbPv9fjzyyCMYNWoUFixYEOtwiIiIiKif61PJ9t/+9jccO3YMjz/+ODQaTazDISIiIqJ+rs8k2ydOnMBf/vIXfPe738XIkSNjHQ4RERERUd+os+3z+fDoo49i+PDh+Pa3v92lfRVFiVBUFA3N7x/fR+oKnjfUHTxvqLt47vRdgiB02KZPJNtPP/00jh8/jjfeeAM6na5L+zY0NPDk7wMaGhpiHQIlIJ431B08b6i7eO70PampqR22Sfhk+9ixY3j66afxrW99C2PHju3y/jabLQJRUbQoioKGhgbYbLZO/bokAnjeUPfwvKHu4rnTvyV8sv3II49g6NCh+P73v9+t/XnS9w2CIPC9pC7jeUPdwfOGuovnTv+U8Mn2sWPHAADjxo0Lef9tt90GAFizZg3LARIRERFRVCV8sr18+fKQ23ft2oWSkhLMnz8faWlpyM7OjnJkFA2CIHRqvBRRazxvqDt43lB38dzp3wSlj84OfPTRR7F+/Xq8/vrrmDhxYqzDISIiIqJ+qM/U2SYiIiIiijdMtomIiIiIIqTPDiMhIiIiIoo19mwTEREREUVIwlcjof7pwIEDWL16Nfbt2wdJkjBixAjcfffd+PKXvxzr0ChO/fOf/8Tu3btx6NAhnDhxApIkYdWqVbj55ptjHRrFsQsXLmDTpk3YsmULzpw5g0uXLsFms2Hy5MlYuXIlJkyYEOsQKQ41NjbiqaeewsGDB1FWVoaGhgakpqZi2LBhuOOOO7Bw4ULW2+5HOIyEEs727dtxzz33QKfTYfHixUhKSsIHH3yAsrIy/PCHP8R9990X6xApDs2fPx/l5eVITU2F2WxGeXk5k23q0H/913/hb3/7G3JycjBt2jSkp6ejtLQUH330ERRFwRNPPIEbb7wx1mFSnCktLcWyZcswYcIE5OTkICUlBTU1Ndi8eTNqampw66234re//W2sw6QoYbJNCcXn8+GGG25AVVUVXn/9dYwZMwYAYLfb8bWvfQ1nz57FO++8g7y8vNgGSnFn27ZtyM3NRXZ2NtauXYsnnniCyTZ16IMPPkBaWhqmTp2q2r5r1y584xvfgMViwWeffQa9Xh+jCCke+f1+KIoCrVY9gMBut+O2227DqVOnsHHjRowcOTJGEVI0ccw2JZTi4mKcO3cOS5YsaUm0AcBqteJ73/sefD4f3n777RhGSPFq9uzZXNyKumzhwoVBiTYATJ06FTNmzEB9fT2OHz8eg8gonmk0mqBEGwh8V82dOxdAoPeb+gcm25RQduzYAQAtH1atzZkzR9WGiCiSmpOpUEkVUSgejwfFxcUQBAEjRoyIdTgUJfyEoIRSUlICAMjNzQ26z2azITU1lb0FRBRxFRUV2LZtGzIzM1FQUBDrcChONTY24oUXXoAsy6ipqcGWLVtQWVmJBx54gMMd+xEm25RQ7HY7ACApKSnk/VarFVVVVdEMiYj6GUmS8PDDD8Pr9eJHP/oRNBpNrEOiONXY2Ii//OUvLbd1Oh0efvhhfOtb34phVBRtTLaJiIg6SZZlPPbYY9i5cyduvfVWLFu2LNYhURwbMmQIjh8/Dr/fj8rKSrz77rv405/+hL179+LPf/4zhyD1ExyzTQnFarUCAJqamkLeb7fbw/Z6ExH1hKIo+NnPfoYNGzZg6dKl+PWvfx3rkChBaDQaDBkyBN/+9rfx4IMP4sMPP8Qbb7wR67AoSphsU0JpHuMWalx2Q0MD6urqQo7nJiLqieYe7bfeegtLlizB73//e4giv0Kp65on+HMyf//BTwpKKNOmTQMAfP7550H3bd26FQAwffr0qMZERH2bLMv46U9/irfffhs33ngj/vCHP3CcNnXbhQsXAIDnUD/CZJsSyqxZszB06FBs3LgRR48ebdlut9vx17/+FVqtFjfddFMMIySivqR1on399dfjj3/8I5Mk6tDRo0dDDnesr6/Hn/70JwDAvHnzoh0WxQhXkKSEU1xcjJUrV0Kn02HJkiWwWq0ty7U/+OCD+O53vxvrECkOvfnmm9i9ezcA4MSJEzh8+DAmT57cMuxowYIFWLBgQSxDpDi0evVq/OUvf4HZbMaKFStCTmhbsGABRo8eHYPoKF797ne/w7p16zBjxgwMHjwYJpMJFRUV+OSTT+B0OrFo0SL8+c9/5lCkfoLTYCnhzJw5E6+88gqeeuopbNq0CZIkYcSIEfjBD36ApUuXxjo8ilO7d+/G+vXrVdv27NmDPXv2AACys7OZbFOQ8vJyAIDT6cTTTz8dsk12djaTbVJZtGgR7HY79u3bh507d8LtdsNms2HKlClYtmwZFi9eDEEQYh0mRQl7tomIiIiIIoTXL4iIiIiIIoTJNhERERFRhDDZJiIiIiKKECbbREREREQRwmSbiIiIiChCmGwTEREREUUIk20iIiIioghhsk1EREREFCFMtomIiIiIIoTJNhHRZWVlZSgsLMSjjz4a61D6jGi8pnzfiCieMdkmIqK4N3/+fMyfPz/WYRARdZk21gEQERH1xMCBA/Huu+8iKSkp1qEQEQVhsk1ERAlNp9Nh+PDhsQ6DiCgkDiMhIuqEiooKPPbYY7jqqqtQVFSEefPm4bHHHkNlZWXI9seOHcO9996LSZMmYcqUKbj33ntx4sQJPProoygsLERZWVmHj7l9+3YUFhZi9erV2LlzJ+68805MmjQJ06dPx3/8x3+gqqoq5H4nT57Egw8+iFmzZqGoqAjz58/H448/jvr6+qC2zcMzGhoa8LOf/QyzZ8/G+PHjsXz5cvz73/8Oat9e/KtXr0ZhYSG2b9/e4XM7dOgQfvOb32DJkiWYMmUKxo8fjy9/+ctYu3YtJElqadc8Hru8vBzl5eUoLCxs+W/16tWqNqHGbHflfbvrrrtQWFgIn8+HNWvWYP78+SgqKsKiRYvwv//7vx0+JyKiUNizTUTUgZKSEtx+++2oqanBtddei5EjR+LkyZN466238Mknn+DVV19Fbm5uS/tjx47h9ttvh9vtxpe+9CXk5ubi8OHDuP322zFq1KguP/6+ffvwzDPP4JprrsFdd92Fw4cPY+PGjdi9ezfWrVuHjIyMlrZ79uzBPffcA6/Xi0WLFiE7Oxv79u3DCy+8gE8//RSvvfYaUlNTVcf3er345je/CbfbjWXLlqGpqQnvvPMO7r//fvzhD3/A0qVLu//ihfHGG29g8+bNmDZtGubNmwe3240dO3bgiSeewMGDB1sS6eTkZDzwwAN44YUXAAB33313yzGmT5/e7mN09X1r9tBDD+HAgQOYN28eRFHEpk2b8Jvf/AY6nQ633nprL74KRNQvKEREpCiKopw/f14pKChQHnnkEdX2FStWKAUFBcprr72m2v7aa68pBQUFyt13363a/vWvf10pKChQNm3apNr+1FNPKQUFBUpBQYFy/vz5DuMpLi5uaf/GG2+o7lu9erVSUFCg/OQnP2nZ5vf7lYULFyoFBQXKli1bVO2feOIJpaCgQHnsscdU26+99lqloKBAWbFiheL1elu2nzp1Shk/frwydepUpampqWX7I488Ejb+5udXXFzcsi3ca1pWVqb4fD7VNlmWlZ/85CdKQUGBsmvXrqA4r7322pCvU2+9b3feeadSUFCgfPWrX1U959OnTytjxoxRFi1aFPLxiYjaw2EkRETtqKysRHFxMUaMGBHUq3nrrbdi+PDh+OKLL1qGJZSXl2P37t0YM2YMrr/+elX7lStXIiUlpcsxDBs2DMuXLw86VlpaGt555x14vV4AgV7tkpISzJs3D1dddZWq/X333YeUlBRs3LixpX1rP/jBD6DT6VpuDx8+HLfccgsaGxtDDifpqezsbGg0GtU2QRBwxx13AAC++OKLHh2/q+9baw899BCsVmvL7fz8fEyePBlnz56F3W7vUVxE1P8w2SYiaseRI0cAANOmTYMgCKr7BEHA1KlTAQSGjrT+d9KkSUHHMplMKCws7HIMkydPDnpso9GIsWPHwu12o6SkRBVrqOEVZrMZRUVFqvbNdDodJk6cGLRP2+fWm7xeL5577jksX74ckydPxqhRo1BYWIibb74ZAHDx4sUeHb+r71trY8eODdo2cOBAAEBTU1OP4iKi/odjtomI2tHck9l6XHRrmZmZAK4kYc3t09LSQrYPd5z2pKent3usto/d2VibpaSkQBSD+16aHzcSCeb/+3//D5s3b0ZeXh5uvPFGpKenQ6vVorGxES+++GLI3veu6O5rASBkCUGtNvB16ff7exQXEfU/TLaJiNrRPJzg0qVLIe9v3t7crvnf2tradtt3RU1NTbvHak4Ouxprs/r6esiyHJRwNz9u6+SzuZc4VNLZ2aT8wIED2Lx5M+bOnYu1a9eqhpPs27cPL774YqeO057uvhZERL2Nw0iIiNoxevRoAMCuXbugKIrqPkVRsHv3blW75moje/fuDTqWy+XC8ePHuxzDnj17gh7b7Xbj8OHDMBqNyMvLAwCMGTMGALBjx46Qj33o0CEYjUYMGzZMdZ8kSdi3b1/QPrt27VI9JwCw2WwAgAsXLgS1P3r0aKeez/nz5wEA11xzTdC47ebHbEsUxS71Knf1fSMiihQm20RE7Rg8eDBmzJiBkydPYt26dar71q1bh5MnT2LmzJkYNGgQgMDEv8mTJ+PIkSN47733VO3/8Y9/hKx13ZGzZ88GPfbf//531NbWYvHixdDr9QACY7tzcnKwZcsWbNu2TdX+mWeeQV1dnap9a08++aSqvvXp06fx1ltvISkpCdddd13L9qKiIgDA+vXrVfu/9957IZP8UAYPHgwALQlvs5MnT2Lt2rUh97HZbKirq4PH4+n0Y3TlfSMiihQOIyEi6sCvfvUr3H777fj5z3+OzZs3Y8SIETh16hQ+/vhjpKWl4Ve/+pWq/c9//nPccccdeOihh/Dee+8hJycHhw8fxv79+zFt2jTs3Lkz5BjpcObOnYtf//rX+PTTT5Gfn4/Dhw/j888/x6BBg/DQQw+1tBNFEatWrcLKlSvx7W9/u6XO9v79+1FcXIycnBz86Ec/Cjp+ZmYmmpqasGzZMlx99dWw2+3YuHEjPB4Pfvvb36qGWixYsABDhgzB22+/jcrKSowePRpnzpxBcXExrr76anz66acdPp/x48dj/Pjx2LRpE6qrqzFhwgRUVlbi448/xtVXX433338/aJ+ZM2fi0KFDuO+++zBlyhTo9XpMnjy5ZaJjKF1934iIIoE920REHcjPz8dbb72Fm266CQcOHMA//vEPHDhwADfddBPWrVsXNCxjzJgxeOWVVzBr1ix8+umnePnllyGKIl555ZWgsd2dMXHiRDz33HOoq6vDiy++iAMHDmDx4sV49dVXgyYATp06Fa+//jrmz5+PrVu34tlnn8X58+dx11134fXXXw85cVOv1+PZZ5/FpEmT8H//939Yv349hg0bhjVr1gQtaGM0GvH888/juuuuw4EDB/Daa6/B4/Hg5Zdfxrhx4zr1fDQaDZ555hnccsstOHfuHF5++WWcOnUKDz/8MH784x+H3Od73/sebr31Vpw8eRJr1qzBE0880WF5wK6+b0REkSAobQezERFRRPj9fnzpS1+C2+0OGuYRyvbt27FixQo88MAD+P73vx+RmObPnw8A+PjjjyNyfCKi/o4920REvczn84WsRrJ27VqUl5djwYIFMYiKiIhigWO2iYh6mdPpxLx58zBnzhzk5eXB5/Nh//79OHjwIDIzM/HAAw/EOkQiIooSJttERL3MaDRi+fLlKC4uxs6dO+H1epGZmYnbbrsN999/PwYMGBDrEImIKEo4ZpuIiIiIKEI4ZpuIiIiIKEKYbBMRERERRQiTbSIiIiKiCGGyTUREREQUIUy2iYiIiIgihMk2EREREVGEMNkmIiIiIooQJttERERERBHy/wFUX+8cT2M+fwAAAABJRU5ErkJggg==",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"utils.plot_scatter(np.log(WAFFLE_DIVORCE.Population), WAFFLE_DIVORCE.Divorce, label=None)\n",
"\n",
"utils.plot_errorbar(\n",
" xs=np.log(WAFFLE_DIVORCE.Population),\n",
" ys=WAFFLE_DIVORCE.Divorce,\n",
" error_lower=WAFFLE_DIVORCE[\"Divorce SE\"],\n",
" error_upper=WAFFLE_DIVORCE[\"Divorce SE\"],\n",
" error_width=7,\n",
")\n",
"plt.xlabel(\"log population\")\n",
"plt.ylabel(\"divorce rate\")\n",
"plt.title(\"Error decreases with state population\");"
]
},
{
"cell_type": "markdown",
"id": "6fcbae6d-070c-4fad-af16-783c59d33f87",
"metadata": {},
"source": [
"## Measured divorce model"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "8028569f-a1cb-48e2-8bcb-0597d65f5beb",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"M \n",
"\n",
"M \n",
" \n",
"\n",
"\n",
"D \n",
"\n",
"D \n",
" \n",
"\n",
"\n",
"M->D \n",
" \n",
"\n",
"\n",
"M* \n",
"\n",
"M* \n",
" \n",
"\n",
"\n",
"M->M* \n",
" \n",
"\n",
"\n",
"D* \n",
"\n",
"D* \n",
" \n",
"\n",
"\n",
"D->D* \n",
" \n",
"\n",
"\n",
"A \n",
"\n",
"A \n",
" \n",
"\n",
"\n",
"A->M \n",
" \n",
"\n",
"\n",
"A->D \n",
" \n",
"\n",
"\n",
"A* \n",
"\n",
"A* \n",
" \n",
"\n",
"\n",
"A->A* \n",
" \n",
"\n",
"\n",
"eD \n",
"\n",
"eD \n",
" \n",
"\n",
"\n",
"eD->D* \n",
" \n",
"\n",
"\n",
"eM \n",
"\n",
"eM \n",
" \n",
"\n",
"\n",
"eM->M* \n",
" \n",
"\n",
"\n",
"eA \n",
"\n",
"eA \n",
" \n",
"\n",
"\n",
"eA->A* \n",
" \n",
"\n",
"\n",
"P \n",
"\n",
"P \n",
" \n",
"\n",
"\n",
"P->D \n",
" potential backdoor \n",
" \n",
"\n",
"\n",
"P->eD \n",
"population confounds \n",
" \n",
"\n",
"\n",
"P->eM \n",
" \n",
"\n",
"\n",
"P->eA \n",
" \n",
"\n",
"\n",
"unobserved \n",
"\n",
"unobserved \n",
" \n",
" \n",
" \n"
],
"text/plain": [
""
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"utils.draw_causal_graph(\n",
" edge_list=[\n",
" (\"M\", \"D\"),\n",
" (\"A\", \"M\"),\n",
" (\"A\", \"D\"),\n",
" (\"M\", \"M*\"),\n",
" (\"D\", \"D*\"),\n",
" (\"eD\", \"D*\"),\n",
" (\"eM\", \"M*\"),\n",
" (\"A\", \"A*\"),\n",
" (\"eA\", \"A*\"),\n",
" (\"P\", \"eA\"),\n",
" (\"P\", \"eM\"),\n",
" (\"P\", \"eD\"),\n",
" (\"P\", \"D\"),\n",
" ],\n",
" node_props={\n",
" \"unobserved\": {\"style\": \"dashed\"},\n",
" \"A\": {\"style\": \"dashed\"},\n",
" \"M\": {\"style\": \"dashed\"},\n",
" \"D\": {\"style\": \"dashed\"},\n",
" },\n",
" edge_props={\n",
" (\"P\", \"eA\"): {\"color\": \"blue\"},\n",
" (\"P\", \"eM\"): {\"color\": \"blue\"},\n",
" (\"P\", \"eD\"): {\"color\": \"blue\", \"label\": \"population confounds\", \"fontcolor\": \"blue\"},\n",
" (\"P\", \"D\"): {\n",
" \"style\": \"dashed\",\n",
" \"label\": \" potential backdoor\",\n",
" \"color\": \"red\",\n",
" \"fontcolor\": \"red\",\n",
" },\n",
" (\"A\", \"A*\"): {\"color\": \"red\"},\n",
" (\"P\", \"eA\"): {\"color\": \"purple\"},\n",
" (\"eA\", \"A*\"): {\"color\": \"red\"},\n",
" (\"A\", \"D\"): {\"color\": \"red\"},\n",
" },\n",
" graph_direction=\"LR\",\n",
")"
]
},
{
"cell_type": "markdown",
"id": "d3f9dd66-c976-44e0-a746-07e62a3cabc0",
"metadata": {},
"source": [
"## Thinking like a graph\n",
"\n",
"- **Thinking like regression**: which _predictors_ do I include in the model?\n",
" - stop doing this\n",
" - To be fair, GLMs enable this approach\n",
" - The scientific world is SO MUCH LARGER than GLMs (e.g. think about the difference equation in tool use example)\n",
"- **Thinking like a graph**:how do I model the network of _causes_?\n",
" - e.g. Full Luxury Bayes\n",
" \n",
"## Let's start simpler"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "5b79ccbe-a2a8-40c5-8fda-e56907fc8c7e",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"M \n",
"\n",
"M \n",
" \n",
"\n",
"\n",
"D \n",
"\n",
"D \n",
" \n",
"\n",
"\n",
"M->D \n",
" \n",
"\n",
"\n",
"M* \n",
"\n",
"M* \n",
" \n",
"\n",
"\n",
"M->M* \n",
" \n",
"\n",
"\n",
"D* \n",
"\n",
"D* \n",
" \n",
"\n",
"\n",
"D->D* \n",
" \n",
"\n",
"\n",
"A \n",
"\n",
"A \n",
" \n",
"\n",
"\n",
"A->M \n",
" \n",
"\n",
"\n",
"A->D \n",
"Divorce Model \n",
" \n",
"\n",
"\n",
"eD \n",
"\n",
"eD \n",
" \n",
"\n",
"\n",
"eD->D* \n",
" Divorce Measurement \n",
"Error Model \n",
" \n",
"\n",
"\n",
"unobserved \n",
"\n",
"unobserved \n",
" \n",
" \n",
" \n"
],
"text/plain": [
""
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"utils.draw_causal_graph(\n",
" edge_list=[(\"M\", \"D\"), (\"A\", \"M\"), (\"A\", \"D\"), (\"M\", \"M*\"), (\"D\", \"D*\"), (\"eD\", \"D*\")],\n",
" node_props={\n",
" \"unobserved\": {\"style\": \"dashed\"},\n",
" \"D\": {\"style\": \"dashed\", \"color\": \"red\"},\n",
" \"M\": {\"color\": \"red\"},\n",
" \"A\": {\"color\": \"red\"},\n",
" \"eD\": {\"color\": \"blue\"},\n",
" \"D*\": {\"color\": \"blue\"},\n",
" },\n",
" edge_props={\n",
" (\"A\", \"M\"): {\"color\": \"red\"},\n",
" (\"A\", \"D\"): {\"color\": \"red\", \"label\": \"Divorce Model\", \"fontcolor\": \"red\"},\n",
" (\"M\", \"D\"): {\"color\": \"red\"},\n",
" (\"eD\", \"D*\"): {\n",
" \"color\": \"blue\",\n",
" \"label\": \" Divorce Measurement\\nError Model\",\n",
" \"fontcolor\": \"blue\",\n",
" },\n",
" (\"D\", \"D*\"): {\"color\": \"blue\"},\n",
" },\n",
" graph_direction=\"LR\",\n",
")"
]
},
{
"cell_type": "markdown",
"id": "fe6342d4-36f2-4701-b949-d809972a0a2c",
"metadata": {},
"source": [
"### 2 submodels\n",
"\n",
"#### Divorce Model \n",
"\n",
"$$\n",
"\\begin{align*}\n",
"D_i &\\sim \\text{Normal}(\\mu_i) \\\\\n",
"\\mu_i &= \\alpha + \\beta_A A + \\beta_M M\n",
"\\end{align*}\n",
"$$\n",
"\n",
"\n",
"#### Divorce Measurement Error Model \n",
"$$\n",
"\\begin{align*}\n",
"D^*_i &= D_i + e_D \\\\\n",
"e_{D,i} &\\sim \\text{Normal}(0, S_i) &S_i \\text{ is the standard deviation estimated from the sample}\n",
"\\end{align*}\n",
"$$\n"
]
},
{
"cell_type": "markdown",
"id": "d004052d-31ed-41a6-8b4d-8dcbafaf310c",
"metadata": {},
"source": [
"### Fit the Divorce Measurement Error Model\n",
"\n",
"- two simultaneous regressions\n",
" - One for $D$\n",
" - Note: these are just parameters in the model, they are not observed 🤯\n",
" - One for $D*$"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "36a1cbe4-ac2b-41cf-83de-d477a4c95c85",
"metadata": {},
"outputs": [],
"source": [
"DIVORCE_STD = WAFFLE_DIVORCE[\"Divorce\"].std()\n",
"DIVORCE = utils.standardize(WAFFLE_DIVORCE[\"Divorce\"]).values\n",
"MARRIAGE = utils.standardize(WAFFLE_DIVORCE[\"Marriage\"]).values\n",
"AGE = utils.standardize(WAFFLE_DIVORCE[\"MedianAgeMarriage\"]).values\n",
"DIVORCE_SE = WAFFLE_DIVORCE[\"Divorce SE\"].values / DIVORCE_STD\n",
"STATE_ID, STATE = pd.factorize(WAFFLE_DIVORCE[\"Loc\"])\n",
"N = len(WAFFLE_DIVORCE)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "f8064b5f-f720-4ec9-9aff-a3ab27907a10",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (4 chains in 4 jobs)\n",
"NUTS: [sigma, alpha, beta_A, beta_M, D]\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "e9c32c493ca84e1d834ed757e039e1e1",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"Output()"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" mean \n",
" sd \n",
" hdi_3% \n",
" hdi_97% \n",
" mcse_mean \n",
" mcse_sd \n",
" ess_bulk \n",
" ess_tail \n",
" r_hat \n",
" \n",
" \n",
" \n",
" \n",
" D[AL] \n",
" 1.202 \n",
" 0.365 \n",
" 0.546 \n",
" 1.897 \n",
" 0.006 \n",
" 0.004 \n",
" 3821.0 \n",
" 2742.0 \n",
" 1.0 \n",
" \n",
" \n",
" D[AK] \n",
" 0.705 \n",
" 0.558 \n",
" -0.325 \n",
" 1.781 \n",
" 0.008 \n",
" 0.007 \n",
" 4484.0 \n",
" 2868.0 \n",
" 1.0 \n",
" \n",
" \n",
" D[AZ] \n",
" 0.439 \n",
" 0.343 \n",
" -0.197 \n",
" 1.076 \n",
" 0.005 \n",
" 0.004 \n",
" 5020.0 \n",
" 2553.0 \n",
" 1.0 \n",
" \n",
" \n",
" D[AR] \n",
" 1.443 \n",
" 0.478 \n",
" 0.587 \n",
" 2.366 \n",
" 0.007 \n",
" 0.005 \n",
" 4344.0 \n",
" 2760.0 \n",
" 1.0 \n",
" \n",
" \n",
" D[CA] \n",
" -0.913 \n",
" 0.130 \n",
" -1.147 \n",
" -0.670 \n",
" 0.002 \n",
" 0.001 \n",
" 5309.0 \n",
" 2611.0 \n",
" 1.0 \n",
" \n",
" \n",
" D[CO] \n",
" 0.684 \n",
" 0.398 \n",
" -0.054 \n",
" 1.416 \n",
" 0.006 \n",
" 0.005 \n",
" 5191.0 \n",
" 2965.0 \n",
" 1.0 \n",
" \n",
" \n",
" D[CT] \n",
" -1.394 \n",
" 0.346 \n",
" -2.057 \n",
" -0.781 \n",
" 0.005 \n",
" 0.004 \n",
" 4951.0 \n",
" 3122.0 \n",
" 1.0 \n",
" \n",
" \n",
" D[DE] \n",
" -0.329 \n",
" 0.504 \n",
" -1.253 \n",
" 0.657 \n",
" 0.008 \n",
" 0.008 \n",
" 4483.0 \n",
" 3030.0 \n",
" 1.0 \n",
" \n",
" \n",
" D[DC] \n",
" -1.916 \n",
" 0.598 \n",
" -3.033 \n",
" -0.761 \n",
" 0.010 \n",
" 0.007 \n",
" 3756.0 \n",
" 3303.0 \n",
" 1.0 \n",
" \n",
" \n",
" D[FL] \n",
" -0.628 \n",
" 0.173 \n",
" -0.978 \n",
" -0.321 \n",
" 0.002 \n",
" 0.002 \n",
" 5325.0 \n",
" 2963.0 \n",
" 1.0 \n",
" \n",
" \n",
"
\n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"def plot_change_xy(x_orig, y_orig, x_new, y_new, color=\"C0\"):\n",
" for ii, (xo, yo, xn, yn) in enumerate(zip(x_orig, y_orig, x_new, y_new)):\n",
" label = \"change\" if not ii else None\n",
" plt.plot((xo, xn), (yo, yn), linewidth=2, alpha=0.25, color=color, label=label)\n",
"\n",
"\n",
"D_true = measured_divorce_inference.posterior.mean(dim=(\"chain\", \"draw\"))[\"D\"].values\n",
"\n",
"# Plot raw observations\n",
"utils.plot_scatter(AGE, DIVORCE, color=\"k\", label=\"observed (std)\", zorder=100)\n",
"\n",
"# Plot divorce posterior\n",
"utils.plot_scatter(AGE, D_true, color=\"C0\", label=\"posterior mean\", zorder=100)\n",
"\n",
"plot_change_xy(AGE, DIVORCE, AGE, D_true)\n",
"\n",
"# Add trendline\n",
"# AGE and DIVORCE are standardized, so no need for offset\n",
"trend_slope = np.linalg.lstsq(AGE[:, None], D_true)[0]\n",
"xs = np.linspace(-3, 3, 10)\n",
"ys = xs * trend_slope\n",
"utils.plot_line(xs=xs, ys=ys, linestyle=\"--\", color=\"C0\", label=\"posterior trend\", linewidth=1)\n",
"\n",
"plt.xlabel(\"marriage age (std)\")\n",
"plt.ylabel(\"divorce rate (std)\")\n",
"plt.legend();"
]
},
{
"cell_type": "markdown",
"id": "2a9ea3ef-601c-4871-8cb8-e60682b9ccab",
"metadata": {},
"source": [
"- black points ignore measurement error\n",
"- red points are posterior means for model that captures measurement error in Divorcd rates\n",
"- thin pink lines demonstrate the movement of estimates by modeling measurment error\n",
"- Modeling measurement error in divorce rates shrinks estimates for states more uncertainty (larger standard errors) toward the main trend line (dashed red line)\n",
" - partial pooling information across states\n",
"- **You get this all for free by drawing the graph, and following the rules of probablity**"
]
},
{
"cell_type": "markdown",
"id": "c74d31b3-3ea9-4215-93dc-86b9e71aecca",
"metadata": {},
"source": [
"### Divorce and Marriage Measurement Error Model"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "cc68c649-279d-46b8-ae9a-f39aa056bb4b",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"M \n",
"\n",
"M \n",
" \n",
"\n",
"\n",
"D \n",
"\n",
"D \n",
" \n",
"\n",
"\n",
"M->D \n",
" \n",
"\n",
"\n",
"M* \n",
"\n",
"M* \n",
" \n",
"\n",
"\n",
"M->M* \n",
"Marriage Measurement \n",
"Error Model \n",
" \n",
"\n",
"\n",
"D* \n",
"\n",
"D* \n",
" \n",
"\n",
"\n",
"D->D* \n",
" \n",
"\n",
"\n",
"A \n",
"\n",
"A \n",
" \n",
"\n",
"\n",
"A->M \n",
"Marriage Model \n",
" \n",
"\n",
"\n",
"A->D \n",
"Divorce Model \n",
" \n",
"\n",
"\n",
"eD \n",
"\n",
"eD \n",
" \n",
"\n",
"\n",
"eD->D* \n",
" Divorce Measurement \n",
"Error Model \n",
" \n",
"\n",
"\n",
"unobserved \n",
"\n",
"unobserved \n",
" \n",
" \n",
" \n"
],
"text/plain": [
""
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"utils.draw_causal_graph(\n",
" edge_list=[(\"M\", \"D\"), (\"A\", \"M\"), (\"A\", \"D\"), (\"M\", \"M*\"), (\"D\", \"D*\"), (\"eD\", \"D*\")],\n",
" node_props={\n",
" \"unobserved\": {\"style\": \"dashed\"},\n",
" \"D\": {\"style\": \"dashed\", \"color\": \"red\"},\n",
" \"M\": {\"color\": \"darkcyan\"},\n",
" \"A\": {\"color\": \"darkcyan\"},\n",
" \"eD\": {\"color\": \"blue\"},\n",
" \"D*\": {\"color\": \"blue\"},\n",
" },\n",
" edge_props={\n",
" (\"A\", \"M\"): {\"color\": \"darkcyan\", \"label\": \"Marriage Model\", \"fontcolor\": \"darkcyan\"},\n",
" (\"A\", \"D\"): {\"color\": \"red\", \"label\": \"Divorce Model\", \"fontcolor\": \"red\"},\n",
" (\"M\", \"D\"): {\"color\": \"red\"},\n",
" (\"M\", \"M*\"): {\"label\": \"Marriage Measurement\\nError Model\"},\n",
" (\"eD\", \"D*\"): {\n",
" \"color\": \"blue\",\n",
" \"label\": \" Divorce Measurement\\nError Model\",\n",
" \"fontcolor\": \"blue\",\n",
" },\n",
" (\"D\", \"D*\"): {\"color\": \"blue\"},\n",
" },\n",
" graph_direction=\"LR\",\n",
")"
]
},
{
"cell_type": "markdown",
"id": "8c2b6ba1-a76d-433e-af64-12a4237d8b0b",
"metadata": {},
"source": [
"### 4 Submodels\n",
"\n",
"#### Divorce Model \n",
"\n",
"$$\n",
"\\begin{align*}\n",
"D_i &\\sim \\text{Normal}(\\mu_{Di}, \\sigma_D) \\\\\n",
"\\mu_{Di} &= \\alpha_D + \\beta_{AD} A_i + \\beta_{MD} M_i\n",
"\\end{align*}\n",
"$$\n",
"\n",
"\n",
"#### Divorce Measurement Error Model \n",
"$$\n",
"\\begin{align*}\n",
"D^*_i &= D_i + e_D \\\\\n",
"e_{D,i} &\\sim \\text{Normal}(0, S_{Di}) &S_{Di} \\text{ is the divorce standard deviation estimated from the sample}\n",
"\\end{align*}\n",
"$$\n",
"\n",
"#### Marriage Rate Model \n",
"$$\n",
"\\begin{align*}\n",
"M_i &\\sim \\text{Normal}(\\mu_{Mi}, \\sigma_M) \\\\\n",
"\\mu_{Mi} &= \\alpha_M + \\beta_{AM} A_i\n",
"\\end{align*}\n",
"$$\n",
"\n",
"#### Marriage Rate Measurement Error Model\n",
"$$\n",
"\\begin{align*}\n",
"M^*_i &= M_i + e_M \\\\\n",
"e_{M,i} &\\sim \\text{Normal}(0, S_{Mi}) &S_{Mi} \\text{ is the marriage standard deviation estimated from the sample}\n",
"\\end{align*}\n",
"$$\n"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "508712b5-1c28-4956-9fa6-4957b2eca90a",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (4 chains in 4 jobs)\n",
"NUTS: [sigma_M, alpha_M, beta_AM, M, sigma_D, alpha_D, beta_AD, beta_MD, D]\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "91c78e9aa55143afaa700989344bbb3b",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"Output()"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"D_true = measured_divorce_marriage_inference.posterior.mean(dim=(\"chain\", \"draw\"))[\"D\"].values\n",
"M_true = measured_divorce_marriage_inference.posterior.mean(dim=(\"chain\", \"draw\"))[\"M\"].values\n",
"\n",
"# Observations\n",
"utils.plot_scatter(MARRIAGE, DIVORCE, color=\"k\", label=\"observed (std)\", zorder=100)\n",
"\n",
"# Posterior Means\n",
"utils.plot_scatter(M_true, D_true, color=\"C0\", label=\"posterior mean\", zorder=100)\n",
"\n",
"# Change from modeling measurement error\n",
"plot_change_xy(MARRIAGE, DIVORCE, M_true, D_true)\n",
"\n",
"# Add trendline\n",
"trend_slope = np.linalg.lstsq(M_true[:, None], D_true)[0]\n",
"xs = np.linspace(-3, 3, 10)\n",
"ys = xs * trend_slope\n",
"utils.plot_line(xs=xs, ys=ys, color=\"C0\", linestyle=\"--\", linewidth=1, label=\"posterior trend\")\n",
"\n",
"plt.xlabel(\"marriage rate (std)\")\n",
"plt.ylabel(\"divorce rate (std)\")\n",
"plt.legend();"
]
},
{
"cell_type": "markdown",
"id": "60a2e6ac-7c2f-433f-a2a0-37974c95ac52",
"metadata": {},
"source": [
"- When modeling measurement error for both Marriage and Divorce, we get estimates that move in multiple dimensions\n",
"- Shrinkage (thin pink line) is toward the posterior trend (dashed red line)\n",
" - the direction of movement is proporitonal to the uncertainty along the respective dimension\n",
" - the degree of the movement is inversely proportional to the \"quality\" (in terms of std dev) of the data points\n",
"- Again, we get this all for free by writing down joint distribution and leaning on the axioms of probability"
]
},
{
"cell_type": "markdown",
"id": "287f9897-f5e6-49f8-9348-1b2f656992e1",
"metadata": {},
"source": [
"## Compare causal effects for models that do and do not model measurement error"
]
},
{
"cell_type": "markdown",
"id": "e4e87b5c-654a-4d72-bd3a-6cc3f2d2adc2",
"metadata": {},
"source": [
"### Fit model that considers no measurement error"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "3c1921ca-68ae-46c7-965f-78da14ab4817",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (4 chains in 4 jobs)\n",
"NUTS: [sigma, alpha, beta_A, beta_M]\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "723101abaf79481cb7619afa505a2c97",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"Output()"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"_, axs = plt.subplots(1, 2, figsize=(10, 5))\n",
"plt.sca(axs[0])\n",
"az.plot_dist(no_measurement_error_model.posterior[\"beta_A\"], color=\"gray\", label=\"ignoring error\")\n",
"az.plot_dist(\n",
" measured_divorce_marriage_inference.posterior[\"beta_AD\"], color=\"C0\", label=\"with error\"\n",
")\n",
"plt.axvline(0, color=\"k\", linestyle=\"--\", label=\"no effect\")\n",
"plt.xlabel(\"$\\\\beta_A$\")\n",
"plt.ylabel(\"density\")\n",
"plt.legend()\n",
"axs[0].set_title(\"Effect of Age on Divorce Rate\")\n",
"\n",
"plt.sca(axs[1])\n",
"az.plot_dist(no_measurement_error_model.posterior[\"beta_M\"], color=\"gray\", label=\"ignoring error\")\n",
"az.plot_dist(\n",
" measured_divorce_marriage_inference.posterior[\"beta_MD\"], color=\"C0\", label=\"with error\"\n",
")\n",
"plt.axvline(0, color=\"k\", linestyle=\"--\", label=\"no effect\")\n",
"plt.xlabel(\"$\\\\beta_M$\")\n",
"plt.ylabel(\"density\")\n",
"plt.legend()\n",
"axs[1].set_title(\"Effect of Marriage on Divorce Rate\");"
]
},
{
"cell_type": "markdown",
"id": "687588ba-353e-4108-af16-a7dc333f7ccf",
"metadata": {},
"source": [
"#### Again, no clear rule for how measurment error will effect causal estimates\n",
"\n",
"- Estimated effect of age is attenuated when accounting for measurement error (left)\n",
"- Estimated effect of marriage on divorce rate increases when modeling measurement error (right)\n",
"- Effects of age and marriage found in previous examples could be due in part to measurement error and/or population confounds"
]
},
{
"cell_type": "markdown",
"id": "769eac77-2a8b-457c-a0b1-8991b9b31014",
"metadata": {},
"source": [
"## Unpredictable Errors\n",
"\n",
"- Modeling measurment error of Marriage increases the estiamted effect of Marriage on Divorce\n",
" - This is not intuitive\n",
" - your intuitions are the devil 😈\n",
" - Often results for nonlinear interactions are not intuitive -- rely on probability theory\n",
" - Likely due to down-weighting the effect of unreliable, high-uncertainty datapoints, improving the estimate\n",
"- Errors can \"hurt\" or \"help\", depending on goals\n",
" - only honest option is to attempt to model them\n",
" - do the best you can, analyze your data like you wish your collegues would analyze their own.\n"
]
},
{
"cell_type": "markdown",
"id": "4d32ea95-daa5-480f-ac93-d029fe06611e",
"metadata": {},
"source": [
"# Misclassification\n",
"\n",
"## [Paternity in Himba pastoralist culture](https://www.science.org/doi/10.1126/sciadv.aay6195)\n",
"\n",
"- Unusual kinships systems (in a Western sense, anyway)\n",
"- \"Open\" marriages\n",
"- **Estimand**: proportion of children fathered by men in extra-marital relationships $p$\n",
"- **Misclassification**: categorical version of measurement error\n",
" - Paternity Test has some false positive rate $f$ (FPR, say 5%)\n",
" - If the rate of extra-marital paternity is small, it may be on the same order as the FPR\n",
" - thus often can't ignore the misclassification error\n",
" - How do we include misclassification rate?\n"
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "681a04cf-32f6-4f77-8fc2-466ab1e1766c",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"F \n",
"\n",
"F \n",
" \n",
"\n",
"\n",
"X \n",
"\n",
"X \n",
" \n",
"\n",
"\n",
"F->X \n",
" \n",
"\n",
"\n",
"T \n",
"\n",
"T \n",
" \n",
"\n",
"\n",
"F->T \n",
" \n",
"\n",
"\n",
"Xstar \n",
"\n",
"Xstar \n",
" \n",
"\n",
"\n",
"X->Xstar \n",
" \n",
"\n",
"\n",
"T->X \n",
" \n",
"\n",
"\n",
"M \n",
"\n",
"M \n",
" \n",
"\n",
"\n",
"M->X \n",
" \n",
"\n",
"\n",
"M->T \n",
" \n",
"\n",
"\n",
"eX \n",
"\n",
"eX \n",
" \n",
"\n",
"\n",
"eX->Xstar \n",
" \n",
"\n",
"\n",
"unobserved \n",
"\n",
"unobserved \n",
" \n",
" \n",
" \n"
],
"text/plain": [
""
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"utils.draw_causal_graph(\n",
" edge_list=[\n",
" (\"F\", \"X\"),\n",
" (\"F\", \"T\"),\n",
" (\"T\", \"X\"),\n",
" (\"M\", \"T\"),\n",
" (\"M\", \"X\"),\n",
" (\"X\", \"Xstar\"),\n",
" (\"eX\", \"Xstar\"),\n",
" ],\n",
" node_props={\"X\": {\"style\": \"dashed\"}, \"unobserved\": {\"style\": \"dashed\"}},\n",
")"
]
},
{
"cell_type": "markdown",
"id": "d7dd3d63-fa31-47a7-98e9-8cf70864a8fe",
"metadata": {},
"source": [
"- social father $F$\n",
"- mother $M$\n",
"- social ties (dyads) $T$\n",
"- _actual_ extra-marital paternity $X$, unobserved\n",
"- _measured_ extra-marital paternity $X^*$\n",
"- misclassification errors $e_X$\n",
"\n",
"### Generative Model\n",
"$$\n",
"\\begin{align*}\n",
"X_i &\\sim \\text{Bernoulli}(p_i) \\\\\n",
"\\text{logit}(p_i) &= \\alpha + \\mu_{M[i]} + \\delta_{T[i]}\n",
"\\end{align*}\n",
"$$\n",
"\n",
"### Measurement Model\n",
"\n",
"$$\n",
"\\begin{align*}\n",
"P(X^*=1 | p_i) = p_i + (1 - p_i)f \\\\\n",
"P(X^*=0 | p_i) = (1 - p_i)(1 - f)\n",
"\\end{align*}\n",
"$$\n",
"\n",
"With\n",
"- error rate $f$\n",
"- probability of extra-marital paternity $p$"
]
},
{
"cell_type": "markdown",
"id": "f4871ae6-f5f6-470b-a7ff-6bc8d681502e",
"metadata": {},
"source": [
"#### Deriving the Misclassification Error Model graphically:\n",
"\n",
"Not being clever, right out all possible measurement outcomes"
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "14afee65-a0f4-4b81-ba78-621d6a493663",
"metadata": {
"jupyter": {
"source_hidden": true
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
" \n",
"\n",
"\n",
"X=1 \n",
"\n",
"X=1 \n",
" \n",
"\n",
"\n",
"->X=1 \n",
" p \n",
" \n",
"\n",
"\n",
"X=0 \n",
"\n",
"X=0 \n",
" \n",
"\n",
"\n",
"->X=0 \n",
" 1-p \n",
" \n",
"\n",
"\n",
"X*=1 \n",
"\n",
"X*=1 \n",
" \n",
"\n",
"\n",
"X=0->X*=1 \n",
" f \n",
" \n",
"\n",
"\n",
"X*=0 \n",
"\n",
"X*=0 \n",
" \n",
"\n",
"\n",
"X=0->X*=0 \n",
" 1-f \n",
" \n",
" \n",
" \n"
],
"text/plain": [
""
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"utils.draw_causal_graph(\n",
" edge_list=[\n",
" (\"\", \"X=1\"),\n",
" (\"\", \"X=0\"),\n",
" (\"X=0\", \"X*=1\"),\n",
" (\"X=0\", \"X*=0\"),\n",
" ],\n",
" node_props={\"\": {\"shape\": \"point\"}},\n",
" edge_props={\n",
" (\"\", \"X=1\"): {\"label\": \" p\"},\n",
" (\"\", \"X=0\"): {\n",
" \"label\": \" 1-p\",\n",
" },\n",
" (\"X=0\", \"X*=1\"): {\n",
" \"label\": \" f\",\n",
" },\n",
" (\"X=0\", \"X*=0\"): {\"label\": \" 1-f\"},\n",
" },\n",
")"
]
},
{
"cell_type": "markdown",
"id": "a23f4ada-7bdf-463a-a4e7-722388f314f3",
"metadata": {},
"source": [
"#### $p(X^*=1) = p + (1-p)f$ \n",
"\n",
"The probability of observing $X=1$ given $p$ is the red path below"
]
},
{
"cell_type": "code",
"execution_count": 22,
"id": "ed3c1c64-7cb3-4bb8-aefd-f3cadc9f5533",
"metadata": {
"jupyter": {
"source_hidden": true
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
" \n",
"\n",
"\n",
"X=1 \n",
"\n",
"X=1 \n",
" \n",
"\n",
"\n",
"->X=1 \n",
" p \n",
" \n",
"\n",
"\n",
"X=0 \n",
"\n",
"X=0 \n",
" \n",
"\n",
"\n",
"->X=0 \n",
" 1-p \n",
" \n",
"\n",
"\n",
"X*=1 \n",
"\n",
"X*=1 \n",
" \n",
"\n",
"\n",
"X=0->X*=1 \n",
" f \n",
" \n",
"\n",
"\n",
"X*=0 \n",
"\n",
"X*=0 \n",
" \n",
"\n",
"\n",
"X=0->X*=0 \n",
" 1-f \n",
" \n",
" \n",
" \n"
],
"text/plain": [
""
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"utils.draw_causal_graph(\n",
" edge_list=[\n",
" (\"\", \"X=1\"),\n",
" (\"\", \"X=0\"),\n",
" (\"X=0\", \"X*=1\"),\n",
" (\"X=0\", \"X*=0\"),\n",
" ],\n",
" node_props={\n",
" \"\": {\"shape\": \"point\"},\n",
" \"X=1\": {\"color\": \"red\"},\n",
" \"X=0\": {\"color\": \"red\"},\n",
" \"X*=1\": {\"color\": \"red\"},\n",
" },\n",
" edge_props={\n",
" (\"\", \"X=1\"): {\"label\": \" p\", \"color\": \"red\"},\n",
" (\"\", \"X=0\"): {\"label\": \" 1-p\", \"color\": \"red\"},\n",
" (\"X=0\", \"X*=1\"): {\"label\": \" f\", \"color\": \"red\"},\n",
" (\"X=0\", \"X*=0\"): {\"label\": \" 1-f\"},\n",
" },\n",
")"
]
},
{
"cell_type": "markdown",
"id": "6b5671df-c027-46c9-b24e-c177e82e9050",
"metadata": {},
"source": [
"#### $p(X^*=0) = (1-p)(1-f)$ \n",
"\n",
"The probability of observing $X=0$ given $p$ is the red path below"
]
},
{
"cell_type": "code",
"execution_count": 23,
"id": "8afd1a64-0070-4483-a1ac-3c2d0ae7d066",
"metadata": {
"jupyter": {
"source_hidden": true
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
" \n",
"\n",
"\n",
"X=1 \n",
"\n",
"X=1 \n",
" \n",
"\n",
"\n",
"->X=1 \n",
" p \n",
" \n",
"\n",
"\n",
"X=0 \n",
"\n",
"X=0 \n",
" \n",
"\n",
"\n",
"->X=0 \n",
" 1-p \n",
" \n",
"\n",
"\n",
"X*=1 \n",
"\n",
"X*=1 \n",
" \n",
"\n",
"\n",
"X=0->X*=1 \n",
" f \n",
" \n",
"\n",
"\n",
"X*=0 \n",
"\n",
"X*=0 \n",
" \n",
"\n",
"\n",
"X=0->X*=0 \n",
" 1-f \n",
" \n",
" \n",
" \n"
],
"text/plain": [
""
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"utils.draw_causal_graph(\n",
" edge_list=[\n",
" (\"\", \"X=1\"),\n",
" (\"\", \"X=0\"),\n",
" (\"X=0\", \"X*=1\"),\n",
" (\"X=0\", \"X*=0\"),\n",
" ],\n",
" node_props={\n",
" \"\": {\"shape\": \"point\"},\n",
" \"X=0\": {\"color\": \"red\"},\n",
" \"X*=0\": {\"color\": \"red\"},\n",
" },\n",
" edge_props={\n",
" (\"\", \"X=1\"): {\"label\": \" p\"},\n",
" (\"\", \"X=0\"): {\"label\": \" 1-p\", \"color\": \"red\"},\n",
" (\"X=0\", \"X*=0\"): {\"label\": \" 1-f\", \"color\": \"red\"},\n",
" (\"X=0\", \"X*=1\"): {\"label\": \" f\"},\n",
" },\n",
")"
]
},
{
"cell_type": "markdown",
"id": "3fdd91b7-4823-4ab5-9b6f-a4d28f988898",
"metadata": {},
"source": [
"### Generating a Proxy Dataset\n",
"\n",
"AFAICT, the dataset in the Scelza et al. paper isn't publicly available, so let's simulate one from the process defined by the generative model."
]
},
{
"cell_type": "code",
"execution_count": 24,
"id": "1477e678-e4fc-4aaf-a64a-48667ad79078",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Actual average paternity rate: 0.48\n",
"Measured average paternity rate: 0.53\n"
]
}
],
"source": [
"from itertools import product as iproduct\n",
"\n",
"np.random.seed(123)\n",
"\n",
"# Generate social network\n",
"N_MOTHERS = 30\n",
"N_FATHERS = 15\n",
"\n",
"MOTHER_IDS = np.arange(N_MOTHERS).astype(int)\n",
"FATHER_IDS = np.arange(N_FATHERS).astype(int)\n",
"\n",
"MOTHER_TRAITS = stats.norm.rvs(size=N_MOTHERS)\n",
"DYADS = np.array(list(iproduct(MOTHER_IDS, FATHER_IDS)))\n",
"\n",
"N_DYADS = len(DYADS)\n",
"DYAD_ID = np.arange(N_DYADS)\n",
"\n",
"RELATIONSHIPS = stats.norm.rvs(size=N_DYADS)\n",
"\n",
"PATERNITY = pd.DataFrame(\n",
" {\n",
" \"mother_id\": DYADS[:, 0],\n",
" \"social_father_id\": DYADS[:, 1],\n",
" \"dyad_id\": DYAD_ID,\n",
" \"relationship\": RELATIONSHIPS,\n",
" }\n",
")\n",
"\n",
"# Generative model\n",
"ALPHA = 0\n",
"BETA_MOTHER_TRAIT = 1\n",
"BETA_RELATIONSHIP = 3\n",
"\n",
"PATERNITY.loc[:, \"mother_trait\"] = PATERNITY.mother_id.apply(lambda x: MOTHER_TRAITS[x])\n",
"p_father = utils.invlogit(\n",
" ALPHA + BETA_MOTHER_TRAIT * PATERNITY.mother_trait + BETA_RELATIONSHIP * PATERNITY.relationship\n",
")\n",
"\n",
"PATERNITY.loc[:, \"p_father\"] = p_father\n",
"PATERNITY.loc[:, \"is_father\"] = stats.bernoulli.rvs(p_father) # unobserved actual paternity\n",
"\n",
"# Measurment model\n",
"FALSE_POSITIVE_RATE = 0.05\n",
"\n",
"\n",
"def p_father_star(row):\n",
" if row.is_father == 1:\n",
" return row.p_father + (1 - row.p_father) * FALSE_POSITIVE_RATE\n",
" return 1 - (1 - row.p_father) * (1 - FALSE_POSITIVE_RATE)\n",
"\n",
"\n",
"PATERNITY.loc[:, \"p_father*\"] = PATERNITY.apply(p_father_star, axis=1)\n",
"PATERNITY.loc[:, \"is_father*\"] = stats.bernoulli.rvs(p=PATERNITY[\"p_father*\"].values)\n",
"\n",
"print(\"Actual average paternity rate: \", PATERNITY[\"is_father\"].mean().round(2))\n",
"print(\"Measured average paternity rate: \", PATERNITY[\"is_father*\"].mean().round(2))"
]
},
{
"cell_type": "markdown",
"id": "d6e41c0d-8e42-418c-b10c-40ccb5a1ce6b",
"metadata": {},
"source": [
"We print out the actual, and observed paternity rates above. In principle, we should be able to recover these with our model"
]
},
{
"cell_type": "code",
"execution_count": 25,
"id": "6d7b2ac1-60e0-4e1c-a5a7-c3ffa303ca67",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" mother_id \n",
" social_father_id \n",
" dyad_id \n",
" relationship \n",
" mother_trait \n",
" p_father \n",
" is_father \n",
" p_father* \n",
" is_father* \n",
" \n",
" \n",
" \n",
" \n",
" 0 \n",
" 0 \n",
" 0 \n",
" 0 \n",
" -0.255619 \n",
" -1.085631 \n",
" 0.135581 \n",
" 0 \n",
" 0.178802 \n",
" 0 \n",
" \n",
" \n",
" 1 \n",
" 0 \n",
" 1 \n",
" 1 \n",
" -2.798589 \n",
" -1.085631 \n",
" 0.000076 \n",
" 0 \n",
" 0.050072 \n",
" 0 \n",
" \n",
" \n",
" 2 \n",
" 0 \n",
" 2 \n",
" 2 \n",
" -1.771533 \n",
" -1.085631 \n",
" 0.001658 \n",
" 0 \n",
" 0.051575 \n",
" 0 \n",
" \n",
" \n",
" 3 \n",
" 0 \n",
" 3 \n",
" 3 \n",
" -0.699877 \n",
" -1.085631 \n",
" 0.039724 \n",
" 0 \n",
" 0.087738 \n",
" 0 \n",
" \n",
" \n",
" 4 \n",
" 0 \n",
" 4 \n",
" 4 \n",
" 0.927462 \n",
" -1.085631 \n",
" 0.845111 \n",
" 1 \n",
" 0.852855 \n",
" 0 \n",
" \n",
" \n",
" ... \n",
" ... \n",
" ... \n",
" ... \n",
" ... \n",
" ... \n",
" ... \n",
" ... \n",
" ... \n",
" ... \n",
" \n",
" \n",
" 445 \n",
" 29 \n",
" 10 \n",
" 445 \n",
" -0.581850 \n",
" -0.861755 \n",
" 0.068670 \n",
" 0 \n",
" 0.115236 \n",
" 0 \n",
" \n",
" \n",
" 446 \n",
" 29 \n",
" 11 \n",
" 446 \n",
" -0.659560 \n",
" -0.861755 \n",
" 0.055178 \n",
" 0 \n",
" 0.102419 \n",
" 1 \n",
" \n",
" \n",
" 447 \n",
" 29 \n",
" 12 \n",
" 447 \n",
" 0.750945 \n",
" -0.861755 \n",
" 0.800764 \n",
" 1 \n",
" 0.810726 \n",
" 1 \n",
" \n",
" \n",
" 448 \n",
" 29 \n",
" 13 \n",
" 448 \n",
" -2.438461 \n",
" -0.861755 \n",
" 0.000281 \n",
" 0 \n",
" 0.050267 \n",
" 0 \n",
" \n",
" \n",
" 449 \n",
" 29 \n",
" 14 \n",
" 449 \n",
" -1.307178 \n",
" -0.861755 \n",
" 0.008299 \n",
" 0 \n",
" 0.057884 \n",
" 0 \n",
" \n",
" \n",
"
\n",
"
450 rows × 9 columns
\n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.subplots(figsize=(10, 4))\n",
"az.plot_dist(paternity_inference.posterior[\"p_father\"], color=\"k\")\n",
"p_father_mean = paternity_inference.posterior[\"p_father\"].mean()\n",
"plt.axvline(\n",
" p_father_mean,\n",
" color=\"k\",\n",
" linestyle=\"--\",\n",
" label=f\"Without Measurement Error: {p_father_mean:1.2f}\",\n",
")\n",
"\n",
"az.plot_dist(paternity_measurement_inference.posterior[\"p_father\"])\n",
"p_father_mean = paternity_measurement_inference.posterior[\"p_father\"].mean()\n",
"plt.axvline(p_father_mean, linestyle=\"--\", label=f\"With Measurement Error: {p_father_mean:1.2f}\")\n",
"plt.xlim([0, 1])\n",
"plt.legend()\n",
"\n",
"plt.legend()\n",
"\n",
"print(\"Actual average paternity rate: \", PATERNITY[\"is_father\"].mean().round(2))\n",
"print(\"Measured average paternity rate: \", PATERNITY[\"is_father*\"].mean().round(2))"
]
},
{
"cell_type": "markdown",
"id": "2b3d595c-ca66-47d0-b4ef-4a8262b87258",
"metadata": {},
"source": [
"- Including measurment error, we're able to accurately recover the _actual_ average paternity rate.\n",
"- Without including the measurement error in the model, we recover the paternity rate observed in the presence of misclassification error.\n",
"- At least for with this simulated dataset, including measurment error gives us a much tigher posterior on the estimated paternity rate, when compared to not modeling the misclassification error in the model."
]
},
{
"cell_type": "markdown",
"id": "8ddee423-2126-4386-90de-b8a26613d37f",
"metadata": {},
"source": [
"## Measurement & Misclassification Horizons\n",
"\n",
"There are a number of problems and solutions related to modeling measurment error and misclassification\n",
"\n",
"- Item response theory (IRT)\n",
" - e.g. NJ wine judging example earlier\n",
"- Factor Analysis\n",
"- Hurdle models\n",
" - measurments need to cross a threshold before they are identifiable\n",
"- Occupancy models\n",
" - considering the existence of phenomena without detection\n",
" - big in applied ecology: just because a species isn't detected doesn't mean it isn't there"
]
},
{
"cell_type": "markdown",
"id": "9a1d7494-5445-4c65-b391-870f235a1805",
"metadata": {},
"source": [
"# BONUS: Floating Point Monsters\n",
"\n",
"- Working on log scale avoids (er, minimizes) issues associated with floating point arithmetic overflow (rounding to one) or underflow (rounding to zero)\n",
"- \"ancient weapons\": special functions for working log scale:\n",
" - `pm.logsumexp`: efficiently computes the log of the sum of exponentials of input elements.\n",
" - `pm.math.log1mexp`: calculates $\\log(1 - \\exp(-x))$\n",
" - `pm.log(1 - p)`\n",
"\n",
"### Logs make sense\n",
"\n",
"$$\n",
"\\begin{align*}\n",
"P(X^*=0) &= (1-p)(1-f) \\\\\n",
"\\log P(X^*=0) & = \\log (1-p) + \\log (1-f) &\\text{put on log scale (used by HMC)}\n",
"\\end{align*}\n",
"$$\n",
"\n",
"### Devil's in the detail\n",
"\n",
"Some terms can become numerically unstable\n",
"\n",
"- e.g. if $p \\approx 0, \\log(1-p) \\approx \\log(1) = 0$\n",
"- e.g. if $p \\approx 1, \\log(1-p) \\approx \\log(0) = -\\infty$"
]
},
{
"cell_type": "code",
"execution_count": 29,
"id": "0a794d91-6cae-4b85-a05a-a8d24454a5a1",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"p: 0.01 , log(1-p): -0.01005033585350145\n",
"p: 1e-10 , log(1-p): -1.000000082790371e-10\n",
"p: 1e-90 , log(1-p): 0.0\n"
]
}
],
"source": [
"for p in [1e-2, 1e-10, 1e-90]:\n",
" print(\"p:\", p, \", log(1-p):\", np.log(1 - p))"
]
},
{
"cell_type": "code",
"execution_count": 30,
"id": "6ff68523-8921-459b-9dff-348458492769",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"p: 0.01 , log1p(1-p): -0.010050335853501442\n",
"p: 1e-10 , log1p(1-p): -1.00000000005e-10\n",
"p: 1e-90 , log1p(1-p): -1e-90\n"
]
}
],
"source": [
"def log1m(x):\n",
" return np.log1p(-x)\n",
"\n",
"\n",
"for p in [1e-2, 1e-10, 1e-90]:\n",
" print(\"p:\", p, \", log1p(1-p):\", log1m(p))"
]
},
{
"cell_type": "markdown",
"id": "86118a97-d54f-4ca8-99c1-cea4abe193b4",
"metadata": {},
"source": [
"### `log1p` and Taylor series approximation for small $p$\n",
"\n",
"- when $p>e{-4}$, just calculate $\\log(1 + p)$\n",
"- when $p 1e-4:\n",
" return np.log(1 + x)\n",
" else:\n",
" # second order Taylor expansion\n",
" return x - x**2 / 2\n",
"\n",
"\n",
"# Compare results to numpy log1p\n",
"for p in [1e-3, 1e-10, 1e-90]:\n",
" print(\"p:\", p, \", numpy.log1p:\", np.log1p(p), \", my_log1p\", my_log1p(p))"
]
},
{
"cell_type": "markdown",
"id": "4eb0ae38-44b2-4469-a080-6f52afa25b60",
"metadata": {},
"source": [
"## `logsumexp`\n",
"\n",
"- Use when you need to take the sum of log of multiple terms.\n",
"- Want the log because it will be numerically stable\n",
"- But logs don't apply to sum, only products, e.g.\n",
"\n",
"$$\n",
"\\begin{align*}\n",
"P(X^*=1) &= p + (1-p)f \\\\\n",
"\\log P(X^*=1) &= \\log [p + (1-p)f] \\\\\n",
"\\log P(X^*=1) &= \\text{pm.math.logsumexp}([\\text{pm.math.log}(p), \\text{pm.math.log}(1-p) + \\text{pm.math.log}(f))]) \\\\\n",
"\\log P(X^*=1) &= \\text{pm.math.logsumexp}([\\text{pm.math.log}(p), \\text{log1m}(p)) + \\text{pm.math.log}(f))])\n",
"\\end{align*}\n",
"$$\n"
]
},
{
"cell_type": "markdown",
"id": "a6917cfe-6350-4cb4-9889-31f04a68f5ae",
"metadata": {},
"source": [
"## Previous Paternity measurement model with log-scaling tricks"
]
},
{
"cell_type": "code",
"execution_count": 32,
"id": "09f90af0-43d3-454d-8fad-3cac68d6093f",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array(-2.9957323, dtype=float32)"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p_father = 0.5\n",
"# pm.math.logsumexp((pm.math.log(p_father), log1m(p_father) + )).eval()\n",
"pm.math.log(FALSE_POSITIVE_RATE).astype(\"float32\").eval()"
]
},
{
"cell_type": "code",
"execution_count": 33,
"id": "2f50eb61-e351-490f-bbfa-5265256be901",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (4 chains in 4 jobs)\n",
"NUTS: [beta_M, beta_T, alpha, T, M]\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "60a00c90ff1542a7a79928cb023a0396",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"Output()"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Original vanilla parameterization\n",
"az.plot_dist(paternity_measurement_inference.posterior[\"p_father\"])\n",
"p_father_mean = paternity_measurement_inference.posterior[\"p_father\"].mean()\n",
"plt.axvline(p_father_mean, linestyle=\"--\", label=f\"Vanilla parameterization: {p_father_mean:1.2f}\")\n",
"\n",
"# Re-parameterization\n",
"az.plot_dist(paternity_measurement_lse_inference.posterior[\"p\"], color=\"C1\")\n",
"p_father_mean = paternity_measurement_lse_inference.posterior[\"p\"].mean()\n",
"plt.axvline(\n",
" p_father_mean, linestyle=\"--\", color=\"C1\", label=f\"Re-parameterization: {p_father_mean:1.2f}\"\n",
")\n",
"plt.xlim([0, 1])\n",
"\n",
"plt.legend()\n",
"plt.title(\"Different Parameterizations,\\nEquivalent results\");"
]
},
{
"cell_type": "markdown",
"id": "9ba9b678",
"metadata": {},
"source": [
"## Authors\n",
"* Ported to PyMC by Dustin Stansbury (2024)\n",
"* Based on Statistical Rethinking (2023) lectures by Richard McElreath"
]
},
{
"cell_type": "code",
"execution_count": 35,
"id": "cc0bd5ef",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Last updated: Thu Feb 20 2025\n",
"\n",
"Python implementation: CPython\n",
"Python version : 3.12.9\n",
"IPython version : 8.32.0\n",
"\n",
"pytensor: 2.27.1\n",
"aeppl : not installed\n",
"xarray : 2025.1.2\n",
"\n",
"arviz : 0.19.0\n",
"xarray : 2025.1.2\n",
"scipy : 1.15.2\n",
"pytensor : 2.27.1\n",
"numpy : 1.26.4\n",
"pymc : 5.20.1\n",
"statsmodels: 0.14.4\n",
"matplotlib : 3.10.0\n",
"pandas : 2.2.3\n",
"\n",
"Watermark: 2.5.0\n",
"\n"
]
}
],
"source": [
"%load_ext watermark\n",
"%watermark -n -u -v -iv -w -p pytensor,aeppl,xarray"
]
},
{
"cell_type": "markdown",
"id": "e8e4f581",
"metadata": {},
"source": [
":::{include} ../page_footer.md\n",
":::"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "default",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.9"
}
},
"nbformat": 4,
"nbformat_minor": 5
}