{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Matching Model\n",
"To find the expected effect of the intervention on the population, we match each treated individual with one or more untreated individuals which are \"almost the same\" as him or her."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"from causallib.estimation import IPW, Matching\n",
"import matplotlib.pyplot as plt\n",
"import seaborn as sb\n",
"import pandas as pd\n",
"import numpy as np\n",
"from causallib.evaluation.metrics import calculate_covariate_balance\n",
"from sklearn.linear_model import LogisticRegression\n",
"from causallib.preprocessing.transformers import PropensityTransformer, MatchingTransformer\n",
"from causallib.datasets import load_nhefs\n",
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Data:\n",
"The effect of quitting to smoke on weight loss. \n",
"Data example is taken from [Hernan and Robins Causal Inference Book](https://www.hsph.harvard.edu/miguel-hernan/causal-inference-book/)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"When we are looking for nearby data points to match against, the one hot encoding may not be the best choice. Augmented features are also not needed and may introduce bias (eg, if we inlude `age^2`, we will call one pair of subjects more distant in the age variable than another pair if the base age is older, even though the difference is the same ). So we do not augment the continuous variables, and instead of one hot encoding the categorical variables, we binarize them to \"high/low\" values."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def binarize(df, column_name):\n",
" df = df.copy()\n",
" m = df[column_name].median()\n",
" def balance(i): return np.abs(0.5 - (df[column_name] < i).sum()/len(df))\n",
" mstar = min([m-1, m, m+1], key=balance)\n",
" df = df.assign(**{column_name: (df[column_name] < mstar).astype(int)})\n",
" df = df.rename(columns={column_name: column_name + f\"<{mstar}\"})\n",
" return df\n",
"\n",
"\n",
"def get_matching_data():\n",
" data = load_nhefs(onehot=False, augment=False)\n",
" data.X = binarize(data.X, \"education\")\n",
" data.X = binarize(data.X, \"exercise\")\n",
" data.X = binarize(data.X, \"active\")\n",
" return data\n",
"\n",
"\n",
"binarized_data = get_matching_data()\n",
"X, a, y = binarized_data.X, binarized_data.a, binarized_data.y"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
active<1.0
\n",
"
age
\n",
"
education<3.0
\n",
"
exercise<2.0
\n",
"
race
\n",
"
sex
\n",
"
smokeintensity
\n",
"
smokeyrs
\n",
"
wt71
\n",
"
qsmk
\n",
"
wt82_71
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
1
\n",
"
42
\n",
"
1
\n",
"
0
\n",
"
1
\n",
"
0
\n",
"
30
\n",
"
29
\n",
"
79.04
\n",
"
0
\n",
"
-10.093960
\n",
"
\n",
"
\n",
"
1
\n",
"
1
\n",
"
36
\n",
"
1
\n",
"
1
\n",
"
0
\n",
"
0
\n",
"
20
\n",
"
24
\n",
"
58.63
\n",
"
0
\n",
"
2.604970
\n",
"
\n",
"
\n",
"
2
\n",
"
1
\n",
"
56
\n",
"
1
\n",
"
0
\n",
"
1
\n",
"
1
\n",
"
20
\n",
"
26
\n",
"
56.81
\n",
"
0
\n",
"
9.414486
\n",
"
\n",
"
\n",
"
3
\n",
"
0
\n",
"
68
\n",
"
1
\n",
"
0
\n",
"
1
\n",
"
0
\n",
"
3
\n",
"
53
\n",
"
59.42
\n",
"
0
\n",
"
4.990117
\n",
"
\n",
"
\n",
"
4
\n",
"
0
\n",
"
40
\n",
"
1
\n",
"
1
\n",
"
0
\n",
"
0
\n",
"
20
\n",
"
19
\n",
"
87.09
\n",
"
0
\n",
"
4.989251
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" active<1.0 age education<3.0 exercise<2.0 race sex smokeintensity \\\n",
"0 1 42 1 0 1 0 30 \n",
"1 1 36 1 1 0 0 20 \n",
"2 1 56 1 0 1 1 20 \n",
"3 0 68 1 0 1 0 3 \n",
"4 0 40 1 1 0 0 20 \n",
"\n",
" smokeyrs wt71 qsmk wt82_71 \n",
"0 29 79.04 0 -10.093960 \n",
"1 24 58.63 0 2.604970 \n",
"2 26 56.81 0 9.414486 \n",
"3 53 59.42 0 4.990117 \n",
"4 19 87.09 0 4.989251 "
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"binarized_data.X.join(binarized_data.a).join(binarized_data.y).head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can run either a Euclidean or a Mahalanobis metric match and predict the individual outcome using `MatchingIndividualOutcomeEstimator`:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"m_euclid = Matching(metric=\"euclidean\").fit(X, a, y)\n",
"m_mahalanobis = Matching(metric=\"mahalanobis\").fit(X, a, y)\n",
"Y_euclid = m_euclid.estimate_individual_outcome(X, a)\n",
"Y_mahalanobis = m_mahalanobis.estimate_individual_outcome(X, a)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can see that the two metrics lead to very similar results on a population level."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0 1.712532\n",
"1 5.562541\n",
"ATE 3.850009\n",
"dtype: float64"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Y_euclid.assign(ATE=Y_euclid[1]-Y_euclid[0]).mean()"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0 1.862848\n",
"1 5.084078\n",
"ATE 3.221230\n",
"dtype: float64"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Y_mahalanobis.assign(ATE=Y_mahalanobis[1]-Y_mahalanobis[0]).mean()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If we inspect the individual counterfactuals, we find, as expected that both metrics return the same value for the observed outcome but differ in the unobserved outcome:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
0_euclidean
\n",
"
1_euclidean
\n",
"
0_mahalanobis
\n",
"
1_mahalanobis
\n",
"
qsmk
\n",
"
\n",
"
\n",
"
sample_id
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
" \n",
" \n",
"
\n",
"
270
\n",
"
-1.585221
\n",
"
-14.515030
\n",
"
-1.358036
\n",
"
-14.515030
\n",
"
1
\n",
"
\n",
"
\n",
"
1231
\n",
"
-1.478036
\n",
"
-2.152889
\n",
"
-1.478036
\n",
"
6.354588
\n",
"
0
\n",
"
\n",
"
\n",
"
312
\n",
"
4.994384
\n",
"
5.214677
\n",
"
4.994384
\n",
"
2.611378
\n",
"
0
\n",
"
\n",
"
\n",
"
747
\n",
"
8.052550
\n",
"
20.983137
\n",
"
8.052550
\n",
"
-20.751731
\n",
"
0
\n",
"
\n",
"
\n",
"
1286
\n",
"
-0.453769
\n",
"
4.653709
\n",
"
-0.453769
\n",
"
7.375747
\n",
"
0
\n",
"
\n",
"
\n",
"
660
\n",
"
2.044486
\n",
"
12.582741
\n",
"
2.044486
\n",
"
21.207697
\n",
"
0
\n",
"
\n",
"
\n",
"
97
\n",
"
14.738576
\n",
"
7.706729
\n",
"
14.738576
\n",
"
6.917506
\n",
"
0
\n",
"
\n",
"
\n",
"
661
\n",
"
-14.289883
\n",
"
-2.152889
\n",
"
3.518181
\n",
"
-2.152889
\n",
"
1
\n",
"
\n",
"
\n",
"
1097
\n",
"
-2.610177
\n",
"
6.348855
\n",
"
-2.610177
\n",
"
-2.722508
\n",
"
0
\n",
"
\n",
"
\n",
"
1462
\n",
"
8.280703
\n",
"
3.058078
\n",
"
8.280703
\n",
"
3.850703
\n",
"
0
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" 0_euclidean 1_euclidean 0_mahalanobis 1_mahalanobis qsmk\n",
"sample_id \n",
"270 -1.585221 -14.515030 -1.358036 -14.515030 1\n",
"1231 -1.478036 -2.152889 -1.478036 6.354588 0\n",
"312 4.994384 5.214677 4.994384 2.611378 0\n",
"747 8.052550 20.983137 8.052550 -20.751731 0\n",
"1286 -0.453769 4.653709 -0.453769 7.375747 0\n",
"660 2.044486 12.582741 2.044486 21.207697 0\n",
"97 14.738576 7.706729 14.738576 6.917506 0\n",
"661 -14.289883 -2.152889 3.518181 -2.152889 1\n",
"1097 -2.610177 6.348855 -2.610177 -2.722508 0\n",
"1462 8.280703 3.058078 8.280703 3.850703 0"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Y_euclid.join(Y_mahalanobis, lsuffix=\"_euclidean\",\n",
" rsuffix=\"_mahalanobis\").join(a).sample(10)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Propensity Matching\n",
"\n",
"To do propensity score matching, we can supply a transformer that replaces the covariates with a learned propensity model, using a given learner."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"propensity_transform = PropensityTransformer(\n",
" learner=LogisticRegression(\n",
" solver=\"liblinear\",\n",
" class_weight=\"balanced\"),\n",
" include_covariates=False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this case we will want to use the augmented data to improve the accuracy of the propensity model. We can calculate the ATE:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0 1.726202\n",
"1 4.682503\n",
"dtype: float64"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"augmented_data = load_nhefs()\n",
"X, a, y = augmented_data.X, augmented_data.a, augmented_data.y\n",
"matcher = Matching(propensity_transform=propensity_transform)\n",
"matcher.fit(X, a, y)\n",
"matcher.estimate_population_outcome(X, a)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We have also provided a convenience subclass `PropensityMatching` which makes this common task straightforward:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0 1.726202\n",
"1 4.682503\n",
"dtype: float64"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from causallib.estimation import PropensityMatching\n",
"\n",
"pm = PropensityMatching(learner=LogisticRegression(\n",
" solver=\"liblinear\",\n",
" class_weight=\"balanced\"))\n",
"pm.fit(X, a, y)\n",
"pm.estimate_population_outcome(X, a)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Multiple neighbor match (with replacement)\n",
"As long as we permit replacement, we can allow multiple neighbors to match. We now check how the number of neighbors impacts the ATE."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Using 1 neighbors, the effect is: 2.956\n",
"Using 2 neighbors, the effect is: 3.189\n",
"Using 3 neighbors, the effect is: 3.242\n",
"Using 4 neighbors, the effect is: 3.176\n",
"Using 5 neighbors, the effect is: 3.211\n",
"Using 6 neighbors, the effect is: 3.188\n",
"Using 7 neighbors, the effect is: 3.141\n",
"Using 8 neighbors, the effect is: 3.196\n",
"Using 9 neighbors, the effect is: 3.219\n"
]
}
],
"source": [
"for n in range(1, 10):\n",
" matcher.n_neighbors = n\n",
" matcher.fit(X, a, y)\n",
" Y = matcher.estimate_population_outcome(X, a)\n",
" print(f\"Using {n} neighbors, the effect is: {(Y[1] - Y[0]):.3f}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Replacement \n",
"\n",
"Until now, we have executed all of the matching with replacement, meaning that we can select the same treated sample as a match for multiple control samples or vice versa. If we want to only allow each sample to be used once, we must disallow replacement. \n",
"\n",
"If we mix in-sample and out-of-sample data, we would end up generating different estimated counterfactuals for a set of samples if they were checked all at once compared to if they were checked in subsets. Because of this, we have restricted no-replacement matching to operate on a single dataset only as a `PopulationOutcomeEstimator` called `Matching`."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"matcher = Matching(with_replacement=True, propensity_transform=propensity_transform)\n",
"matcher.fit(X, a, y)\n",
"match_df_with = matcher.match(X, a)\n",
"ATE_with_replacement = matcher.estimate_population_outcome(X, a).diff()[1]"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"matcher = Matching(with_replacement=False, propensity_transform=propensity_transform)\n",
"matcher.fit(X, a, y)\n",
"match_df_without = matcher.match(X, a)\n",
"ATE_without_replacement = matcher.estimate_population_outcome(X, a).diff()[1]"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"With replacement we find:\n",
"2.956\n",
"Without replacement we find:\n",
"3.426\n"
]
}
],
"source": [
"print(\n",
" f\"With replacement we find:\\n{ATE_with_replacement:.3f}\\nWithout replacement we find:\\n{ATE_without_replacement:.3f}\")"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"With IPW we find:\n",
"3.490\n",
"and the naive estimate is:\n",
"2.541\n"
]
}
],
"source": [
"ipw = IPW(LogisticRegression(solver=\"liblinear\"))\n",
"ipw.fit(augmented_data.X, augmented_data.a)\n",
"Yipw = ipw.estimate_population_outcome(\n",
" augmented_data.X, augmented_data.a, augmented_data.y)\n",
"ATE_ipw = Yipw[1] - Yipw[0]\n",
"ATE_naive = y[a == 1].mean() - y[a == 0].mean()\n",
"\n",
"print(\n",
" f\"With IPW we find:\\n{ATE_ipw:.3f}\\nand the naive estimate is:\\n{ATE_naive:.3f}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Caliper\n",
"\n",
"Often we want to impose a restriction on the proximity of examples so that we do not permit any match if there is more than a distance $\\kappa$ between the samples. We call this a \"caliper\" and we can see how it impacts the predicted effect here:"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"caliper = np.logspace(-3, 0, 20)\n",
"\n",
"\n",
"def check_caliper(c, with_replacement=True):\n",
" matcher = Matching(propensity_transform=propensity_transform,\n",
" caliper=c, with_replacement=with_replacement)\n",
" matcher.fit(augmented_data.X, augmented_data.a, augmented_data.y)\n",
" Y = matcher.estimate_population_outcome(\n",
" augmented_data.X, augmented_data.a,)\n",
" p = matcher.samples_used_.sum() / len(augmented_data.y)\n",
" return p, (Y[1] - Y[0])\n",
"\n",
"\n",
"p_with, ATE_with = zip(\n",
" *[check_caliper(c, with_replacement=True) for c in caliper])\n",
"p_without, ATE_without = zip(\n",
" *[check_caliper(c, with_replacement=False) for c in caliper])"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA7AAAAF3CAYAAACcz6fnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAB9Y0lEQVR4nO3dd3RURRvH8e9mk930AoSE3juhCYIgKk2kihRFEAEpigIqKsWKCogFFF8V6SAIKhZEQAVBQUVFpITeAwSS0BLSk23vHwvBSAmQsim/zzk5u3Pv7L1P9iaZPDtzZwwOh8OBiIiIiIiISD7n5uoARERERERERK6HElgREREREREpEJTAioiIiIiISIGgBFZEREREREQKBCWwIiIiIiIiUiAogRUREREREZECwd3VAdwou92OzaaVf0REJGd4eBhdHUKBp7ZZRERy0rXa5gKXwNpsDuLikl0dhoiIFBLBwX6uDqHAU9ssIiI56Vpts4YQi4iIiIiISIFQ4HpgRURE5OakpaXRt29f0tPTsdlstG/fnpEjR2aqc/LkScaMGUNCQgI2m41nn32WO++800URi4iIZKYEVkREpIgwmUwsWLAAHx8fLBYLffr04Y477qBBgwYZdaZPn06HDh3o06cPBw8eZOjQoaxbt851QYuIiPyLhhCLiIgUEQaDAR8fHwCsVitWqxWDwXBZncTERAASEhIoWbJknscpIiJyNeqBFRERKUJsNhvdu3fn2LFj9OnTh/r162faP3z4cAYNGsSiRYtISUlh3rx5LopURETkcuqBFRERKUKMRiPffvst69evJzw8nP3792fav3LlSu677z42bNjAzJkzGT16NHa73UXRioiIZKYEVkREpAjy9/enadOm/Prrr5m2f/nll3To0AGAhg0bkpaWRmxsrCtCFBERuYwSWBERkSLi3LlzxMfHA5CamsrGjRupXLlypjqlSpXijz/+AODQoUOkpaVRrFixPI9VRETkSnQPrIiISBFx6tQpxo4di81mw+FwcM8999CqVSumTZtG3bp1adOmDWPHjuXFF19k/vz5GAwGJk+efNlETyIiIq5icDgcjtw48Lhx4/jll18oXrw4K1asuGy/w+Fg4sSJrF+/Hk9PTyZPnkydOnWyPK7FYiMuLjk3QhYRkSIoONjP1SEUeGqbRUQkJ12rbc61IcTdu3dn9uzZV92/YcMGIiIiWL16Na+//jrjx4/PrVBERERERESkEMi1IcRNmjQhMjLyqvvXrl1Lt27dMBgMNGjQgPj4eE6dOqX15kQKCLsdrFawWJyPVqsBTVQqecXLy4Gvr6ujkCLFng5uJldHISJS5LnsHtiYmBhCQ0MzyqGhocTExCiBlSLDYoHEREhMNJCQYPjPcwMJCZfvS0w0kJKSczE4HGCzgcViuJCEXkpGLZaL+65cx27XPXHiOt7eDrZsSURzC0muc9jw3fssnicXk1D7fdJKPeDqiEREijRN4iSSg86fh+PH3Th+3I3ISAPHjjkfY2MvT0xTU68vAfT2duDj48DPD3x9HXh5OcjJ+VTc3cFsduDhAe7uDtzdwcMDjEYybbv05ci03/norOOmec0ljwQHOwgMdHUUUujZUvHfOQjzqe+weVXEf+cQkhP3kFT1ZTDoD56IiCu4LIENCQkhOjo6oxwdHU1ISIirwhHJksMBsbGXEtTjxw1ERjofL26Lj8+cWXp7Oyhb1k6JEg5Kl3bg42PHz8859NH56Mgo+/peLJOx3cfHmTSKiEjeMlji8N/WG4+4P0is8SYpZQfhu/c5vCOmYkzcQ0LYbBzumgBMRCSvuexf49atW7No0SI6derE9u3b8fPz0/BhyRdSU+Gnn9yJiDBc6Em9lKQmJWVOUH19HZQrZ6d8eQfNmlkoV85OuXKOjMdixXK2t1RERHKfW+pJArb2wJi0n4SwuaSF9gAgsdZ7WH1r47t/LIGb2nK+wWfYvSu5OFoRkaIl15bRGTVqFJs2bSI2NpbixYszYsQIrFYrAA8++CAOh4PXXnuNX3/9FS8vLyZNmkRYWFiWx9VU/ZJbHA74/nt3XnnFzNGjzqFhgYHOHtT/JqbORzsBAShBFSngtIxO9hWmttmYtJ+ALfdhsMQRX/9TLMXvuqyOx9mf8Q/vDwY34ustxFKsZd4HKiJSiF2rbc61BDa3FKZGUvKPffvcePFFM+vXu1Ojho1XXknj1ltt+Pu7OjIRyW1KYLOvsLTN7nGbCNjWCwzunG/4FVb/Bleta0w6iP/23hiTD5NY8x1Syz6Sd4GKiBRyLlkHVqQgiIuDF14wc9dd3mzbZmTSpFR+/jmZtm2VvIqIFCWm0z8Q+E8XHO6BxDZZc83kFcDmU5W4JmtJL9YKvz1P4bv3GbBb8iZYEZEiTAmsFEk2GyxY4MFtt/kwZ44HDz1k4c8/kxg82KJJk0REihjziUX4b38Qq08NYpuswe5d+bpe5/AIIL7hFyRXGInX8VkEbO2OwXIul6MVESnalMBKkfPnn0batfPmuec8qV7dzpo1ybz9dhrFixeo0fQiIpJdDgdeR6biv/txLEF3cL7xShzmG5xQ0mAkqfoE4utMxyP2D4L+aoUxcV/uxCsiIkpgpeg4ccLA0KGedO3qTWysgZkzU1i2LIWwMLurQxMRkbzmsOOzbwy+B8eTGtqT8w2XZmtZnLTSfYlrvBKDLYnATa0xnf4xB4MVEZGLlMBKoZeSAu+8Y6J5cx9++MGdZ55J4/ffk+jWzaoZhEVEiiJ7Gn47HsH7+Mckl3+ChLqzwc2U7cNaA5sS2/QXbN6V8d92P14R7zunuBcRkRyju/2k0HI4YMUKd8aPN3P8uBtdu1p45ZU0ypXTPxMiIkWVwRqP//a+mM6tJ7Ha66RUGJmj66HZPcsS1+QH/HY9ju+BF3FP3EVCrWlg9Myxc4iIFGVKYKVQ2rXLuSzO77+7U7u2jW++SaZFC5urwxIRERcypMUQsLUn7ok7ia/zMWml++TOiYw+JITNx+ZTC5/DkzAmH+R8/cU4zCG5cz4RkSJEQ4ilUDl3DsaMMdOmjTe7dxt5881UfvpJyauISFFnTDpI0N/tcE86QHyDz3Mveb3IYCC5yljO11uIe8Iugv66C/f4bbl7ThGRIkAJrBQKVivMmeNBs2a+fPKJBwMHWvjzz0QGDtSyOCIiRZ37+S0E/n03BmsCcY1XkF7i7jw7d3rIvcTeuhoMBgL/bo85+us8O7eISGGkBFYKtPR0WLLEnZYtfRg3zpOwMBvr1iXzxhtpBAW5OjoREXE1j7NrCfynEw6jN3FNVmMNaJznMdj86hF76y9Y/erhv2MA3gcngEMz4IuI3Az1TUmBlJwMixd78OGHJk6ccKNuXRvz5qXQsaNmFhYRESdz1Bf47XoMm09Nzjf8CrtnKZfF4jCXJK7xCnz3PI3PkbcwONJJqvaay+IRESmolMBKgRIfD/PmmZgxw4MzZ9xo2tTKO++k0rq1TYmriIg4ORx4HX0f3wMvkR50O/H1l+DwCHB1VOBmJrH2h2BwxytiGukl2mMJauHqqEREChSDw1GwFiizWGzExSW7OgzJY6dPG5g1y4M5c0wkJBho3drKU0+l06yZJmcSkewJDvZzdQgFXr5qm62J+O0ejmfM16SV7EZ83Zn5bwkbayLF/mwBDjuxt23E4a6fQRGRf7tW26weWMnXTpww8NFHJhYt8iA1Fbp0sTJyZDr16uneIRERycyYtB//7X0xJh0gsep4Uio+BYZ8ON2Huy/xdWcQ+Pc9+Ox/nsTa/3N1RCIiBYYSWMmXDh408L//mVm61Pkj2quXlREj0qhatUANGBARkTxiilmG367Hwc2T842WYSl+l6tDuiZrYDNSKj6Fd8RU0oM7kh7cwdUhiYgUCBpCLPnKjh1uTJtm4rvv3DGb4aGHLDz+eDplyxaoH1MRKUA0hDj7XNo22y34HHgF72MfYAloQny9T7B7lnFNLDfKnkbQX61wSz/Fudv+wmEq7uqIRETyBQ0hlnzvzz+NTJtmYu1ad/z8HIwcmc7QoRaCg5W4iojIlbmlReMXPgBT3EZSyg0lsfokcDO5Oqzr52Ymvu5Mgv66E789TxFf7xM0I6GIyLUpgRWXcTjg55+NvPeeiT//dKd4cTvPP5/GwIHpBOSDySJFRCT/8ojdiF94f9ysCcTXnUVaqQdcHdJNsfnVJanKi/gefAVz9Oeklert6pBERPI1JbDiEj//bGTiRDPh4UbKlLEzcWIqffta8PZ2dWQiIpKvORx4HfsQnwMvYfOqSGyjZdj86rg6qmxJqTgS85nv8d37HJag27F7lnV1SCIi+VY+nJpPCrMdO9zo2dOLBx7wJi7OwLRpKfz1VxJDhih5FRGRazNYE/DbMQDf/c+THtyRuFt/KfDJKwAGI/F1PsbgsDononJopn0RkatRAit54vhxA48/7kmbNj7s2GFkwoRUfv89iQcftGIqQLcriYiIaxgT9xK4qRXmmG9JrPY68fUW4fAoPPeb2L0rk1h9EqZzv+B5fJarwxERybc0hFhyVVwcvPeemdmzPXBzg5Ej0xgxQve4iojI9TNHf4Xf7uE4jN6cv+U7LMVaujqkXJFaZgCmUyvwPfAyluKtsflUc3VIIiL5jnpgJVekpsJHH3lw662+TJ/uQffuVv74I4kXX1TyKiIi18luwWffGPx3DMTqW5fYpr8W2uQVAIOBxDof4jB64rdzKNitro5IRCTfUQIrOcpuhy+/dKdFCx/Gj/ekUSMb69Yl8/77qZQpoyVxRETk+rilRhH4Tye8j00nudxjxDVeid2ztKvDynV2cyiJNd/FI/4fvCOmuDocEZF8R0OIJcds2GDk1VfN7NhhJCzMxtSpydx5p83VYYmISAHjce43/HcMwGBLIj5sLmmhPV0dUp5KC+1O6ukVeB9+k/QSd2P1b+jqkERE8g31wEq27drlRu/eXvTs6ZxZePr0FNasUfIqIiI3yOHAK+J9ArZ0we4eQOyt64pc8npRYs13sJuCnUOJbSmuDkdEJN9QAis37cQJAyNHetK6tTdbthh59VXnzMI9elhx00+WiIjcCIcN//D++B54kfTgzsQ1/Rmbby1XR+UyDo9iJNT5CPekffgcfN3V4YiI5BsaQiw37Px5eP99E7NmmXA44PHHLTz5ZBqBga6OTERECiq3tFO4x/1JYrWJpFQYDgaDq0NyOUvxNqSUHYzXsQ9JD+5QuCewEhG5TgaHw1GgZtaxWGzExSW7OowiKS0N5s/3YOpUM3Fx0LOnlbFj0yhXrkD9CImIZBIc7OfqEAo8tc25yJZE0B8tMDisxN62EYe7v6sjEhHJdddqmzXQU7Jks8Hnn7vTvLkPL73kSViYjZ9+SubDD1OVvIqIiOQmow8JdWfglhqJz75xro5GRMTllMDKVTkc8P337rRq5c2IEV4EBTn4/PNkvvwyhbAwu6vDExERKRKsgU1JrjQKr5MLMZ1a6epwRERcSvfAyhX9/ruRCRPM/POPkSpV7MyenULnzpqcSUSkIEtLS6Nv376kp6djs9lo3749I0eOvKzeqlWr+OCDDzAYDNSsWZMpU7QeqaslVx6L+fSP+O0ewbnAW3GYgl0dkoiISyiBlUy2b3dj4kQzv/ziTunSdqZOTaV3bwvu+kkRESnwTCYTCxYswMfHB4vFQp8+fbjjjjto0KBBRp2IiAhmzpzJkiVLCAgI4OzZs64LWC5xMxEfNougP+/Ab89TxNdbpImuRKRIUn+aAHDwoIHBgz1p186H8HA3xo9P5Y8/knjoISWvIiKFhcFgwMfHBwCr1YrVasXwnyToiy++oG/fvgQEBABQvHjxPI9TrszmW5ukqi9jPvUd5qglrg5HRMQllJoUcSdPGnjnHRNLlnhgNsMzz6Tx+OPp+GlSThGRQslms9G9e3eOHTtGnz59qF+/fqb9ERERAPTu3Ru73c7w4cO54447XBCpXElKhScwnV6F777RWIJaYvcq5+qQRETylHpgi6izZw288oqZpk19+OILDwYNsvD330mMGaPkVUSkMDMajXz77besX7+e8PBw9u/fn2m/zWbj6NGjLFy4kClTpvDSSy8RHx/vomjlMgYjCXWmg8OO3+7HwaFJFUWkaFECW8QkJsI775ho0sSHGTM8uO8+K3/8kcSECWkEB2tJHBGRosLf35+mTZvy66+/ZtoeEhJC69at8fDwoFy5clSsWDGjV1byB7t3JZKqv4Hp3Hq8jn/s6nBERPKUEtgiIi0NZs704NZbfXjrLTN33GFl/fpk3n9fa7mKiBQV586dy+hNTU1NZePGjVSuXDlTnbZt27Jp06aM+hEREZQrp2Gq+U1qmYdJK9EenwPjMSbuc3U4IiJ5RvfAFnI2Gyxd6s5bb5mJjHSjZUsrzz+fwi23aMiRiEhRc+rUKcaOHYvNZsPhcHDPPffQqlUrpk2bRt26dWnTpg0tW7bk999/p2PHjhiNRkaPHk1QUJCrQ5f/MhhIqP0Bxf5oin/4QyTU/gBrYFNXRyUikusMDoejQHW/WSw24uKSXR1GgRAZaaBfPy927TLSoIGNF15I4847ba4OS0QkXwkO1o3/2aW22XU8zv6M385HMaZHkxrSnaRqr2L3quDqsEREsuVabbMS2EJqzx43evf2IjHRwNSpqXTtatVycSIiV6AENvvUNruYNRHvo9PwjngfsJNS/gmSK43C4e7v6shERG7Ktdpm3QNbCP3xh5EuXbxxOGD58mTuvVfJq4iISKHl7ktylRc41+If0kLuwztiKsV+a4Bn5FywW10dnYhIjlICW8isWOHO/fd7ERJiZ+XKZOrU0b2uIiIiRYHdsywJdWcSe+svWH2q47fnKYL+bIHHmZ9cHZqISI5RAluIzJvnwaBBnoSF2fnuu2TNLiwiIlIEWQMacb7x95yvtwiDPYXArd3x39IDY+JeV4cmIpJtSmALAYcDJk82MWaMJ+3a2fjyy2SKFXN1VCIiIuIyBgPpIV051/xvEqtPwuP8JoL+vA3fPU9jSD/j6uhERG5ariawGzZsoH379rRr146ZM2detv/EiRP079+fLl260K9fP6Kjo3MznELJaoVRo8xMnWqmb9905s9Pwdvb1VGJiIhIvuBmJqXCcM612EZK2cF4nphPsd8b4BXxHthSXR2diMgNy7UE1maz8dprrzF79mxWrlzJihUrOHjwYKY6b775Jt26deO7777j8ccfZ8qUKbkVTqGUnAwDBnjx6acmRo1KY+rUNNy1sq+IiIj8h8NUnKSabxN7219YApvje+Bliv1xK6aYb5xDuURECohcS2DDw8OpUKEC5cqVw2Qy0alTJ9auXZupzqFDh2jWrBkAzZo1u2y/XN25c9Czpzdr1hh5881Uxo5N10zDIiIick02n+rEN/yCuEbf4jD6EhDen8C/78b9/GZXhyYicl1yLYGNiYkhNDQ0oxwSEkJMTEymOjVr1mT16tUArFmzhqSkJGJjY3MrpELj+HEDXbp4s2OHG3PmpDJwoMXVIYmIiEgBYineithmv5JQ+wOMKUcI2tQavx2DcEs57urQRESuyaWTOI0ePZq///6bbt26sWnTJkJCQjAaja4MKd/bvduNTp28OXXKjaVLU+jcWeu7iYiIyE0wGEkt8zDnWmwlqdKzmE99R7GNjfDf/jDmqKUYrPGujlBE5DK5dsdkSEhIpkmZYmJiCAkJuazOBx98AEBSUhKrV6/G398/t0Iq8DZuNPLww174+DhYvjyZWrW0xquIiIhkj8Pdj+SqL5NaZiDeEe9iPrUc86llOAwm0ovfRXpwF9JKdsJhKuHqUEVEcq8HNiwsjIiICI4fP056ejorV66kdevWmeqcO3cOu92ZhM2cOZMePXrkVjgF3nffufPAA16EhtpZuVLJq4iIiOQsu1c5EmtN5ewd+4htspqU8o/inrQfvz0jKL6+KgGbO+J1bLqGGYuISxkcjtybem79+vVMmjQJm81Gjx49GDZsGNOmTaNu3bq0adOGH374galTp2IwGGjcuDGvvPIKJpPpmse0WGzExSXnVsj50pw5Hjz/vJnGje0sWpRMUJCrIxIRKTyCg/1cHUKBVxTb5iLD4cCYuANzzHLMp1fgnrgbAIt/Q9JKdiW9ZBdsPtVdHKSIFDbXaptzNYHNDUWpkXQ4YPJkE+++a+aeeyzMmJGKl5eroxIRKVyUwGZfUWqbizpj0gFMp1ZgPrUcj/h/ALD61CCtZBfSS3bF6lefbC+LYEvCLS0aY1oMbmnRuKVH45YWAw4rKeWfwO5ZOge+ExHJz5TAFkBWKzz7rJnFi03065fOm29qjVcRkdygBDb7ikrbLJm5pZ64kMx+h0fsbxiwY/MsT1rJzqSX7IolsCkYLkzO6XBgsMbh9p+k1C0t+kL54vMY3GwJl53LYfBwPhp9SKw5hbTQntlPlEUk31ICW8AkJcHQoV6sWePOs8+m8dxzWuNVRCS3KIHNvqLQNsu1GdLPYjr9PeZTyzGdXYfBkY7dFIzNq5IzKU2PxmBPu+x1Djdv7OYQ7OZQbOZQ7Cbn84vb7KZQ7OZQHB5BGJMP47frUTzO/01qSHcSa07BYSrugu9WRHKbEtgCJC4OHnzQm61b3XjzzTT699caryIiuUkJbPYV9rZZbozBmoDpzBpMp77DzXLWmYxeTEwzElRnkuow+t1YT6rditfRafgcmoTdoxiJtf9HevA9uffNiIhLKIEtQF54wcy8eR7Mnp1Kx45a41VEJLcpgc2+wt42S/5jTAjHf+dQ3BN3k1JmAEnVJ+Jw1++ySGFxrbY515bRkRsXE2Ng4UIP7r/fouRVRERE5CpsfvWIbbqe5IpP43liAUF/tsAj9ndXhyUieUAJbD7y0Ucm0tNh5Mh0V4ciIiIikr+5mUmq9ipxTX4EDARs7ojP/hfAlurqyEQkF2kIcT5x5oyBxo196NjRykcf6Q+viEhe0RDi7MvJtrnbso6Xbeta9T4eqTuEZEsyfVb2vGx/75p96V2zL2dTzjLox36X7R9QZxDdqvXgREIkT6wdetn+YQ1G0L5iBw7GHuDZ9U9etv/pW57jznKt2HEmnJd+G3vZ/uebvsKtpZqyKeovJv316mX7X799MmEl6rH++M+8+8/bl+1/585pVA2qxo8R3zN92/8u2/9hm5mU8SvLsgNfMX/XnMv2z2m/kOJexfls76d8tvfTy/Yv7vQl3h7ezN05i+UHv7ls/7Juq5zn2fo+a47+kGmfp7snn3X+GoApm9/k18j1mfYHeRZj3j2LAJjwx3g2x2zKtL+UT2mmt5sNwIu/jWHnmR2Z9lcJrMqUu94H4JlfRnIo7mCm/XVLhDHh9jcBGLZmMFFJJzPtbxxyKy/eNh6sify99jY6uh3liN2b19Kqs9/uS8uyd/JM4zEA9F7RnVRr5v+x2lW4hycajgT0s6efvZv82QMG/vAQsannMu0vij97F69nTtAQ4gJgxgwPUlLg6afV+yoiIiJyQ9x9+chRn1GpdfDFyizP7QzwOIabw+7qyEQkh6kHNh+IjYVbbvGlbVsrM2eq91VEJC+pBzb7CmPbLAWXwXIO373P4hn9JRb/W0ioOxObTzVXhyUiN0A9sPncrFkmEhMNPPWUel9FREREssPhUYyEsLnEh83HmHKYoD9b4HVsOqg3VqRQUALrYvHxMHOmiY4dLdSurT+sIiIiIjkhLbQ7sbf9RXqxO/DdN4aAf7rilnLc1WGJSDYpgXWxOXNMxMcbGDVKva8iIiIiOcluDiW+wVISar2Pe/wWgv68DfPJT6Fg3UEnIv+iBNaFEhOdkze1a2elXj31voqIiIjkOIOB1LIDiG32O1bfuvjvGob/9j4Y0mJcHZmI3AQlsC40f74H5865MWpUmqtDERERESnU7N6VON94JYnVJmA6s4bivzfA++DrGCznXR2aiNwAzULsIsnJ0LixD3Xq2Fm6NMXV4YiIFFmahTj7CkvbLEWHMekg3ocm4BnzNXaPIJIrPkNKuSFg9HJ1aCKCZiHOlxYt8uDMGTeeeUb3voqIiIjkJZtPVRLqzSe26Qas/o3wPfAixX5viGfkArBbXReYw4ExcR9YE10Xg0g+px5YF0hNhSZNfKha1c4336j3VUTEldQDm32FoW2Wos3j3K/4HByPx/m/sXpXJanqS6SXvBcMedPX45YaiWfUZ5hPLsY9+SAON0/Si7cmrWQX0oM74PAolidxiOQX12qb3fMwDrlg8WIPYmLc+OijVFeHIiIiIlLkWYq1JK7JT5hOr8Ln4GsEhPfH4teQpGqvYCnWCgyGnD+pLQnzqe/wPLkEj3O/YMBBemALEsoPw5h8EPOp7zCfXoXDYMQS1JK0kp1JD+6M3bN0zsciUoCoBzaPpadD06Y+lC7tYMWK5Fz5eygiItdPPbDZV9DbZpFMHDbMUZ/jc2gSxtRjpBe7k6Sqr2ANaJwDx3bgEbcR88nFmGO+wc2WiM2rIqmlepNa6kHs3pUy1XWP34r51HeYTi3HPfkAAJaAJqSV7EJayS7YvatkPyaRfOhabbMS2Dy2aJEHo0Z58tlnybRubXN1OCIiRZ4S2Owr6G2zyBXZ0/CMnIfP4bdws5whrWQXkqq8hM235g0fyi0lAs+TS/CMWoIxJQK70Ze0kG6kle6LJfC26xqqbEzch/n0d5hivsMjYSsAVt86pJXsTFrJrth86+ZOT7GICyiBzSesVrjtNh+Cghz8+KN6X0VE8gMlsNlXkNtmkawYrAl4HfsIr4j3MdiSSCv9IEmVx2H3Kp/l60wx3+IZtRhT7G84MGApdieppR4kLaQrGH1uOia3lGOYT69wJrNxGzHgwOZV8ULPbFesAU3y7P5dkdygBDaf+Pxzd0aM8OKTT5K55x71voqI5AdKYLOvILfNItfLkH4W74ipeB2fCQ4HKeUGk1zpWRymEpcqOex4xP6K58lPMccsx2BPxupVmbTSfUkt1Ru7V7lciOs05lOrMJ1ajuncLxgcFmymENJLdiatZBesAU1wKJmVvGDwADdTjhxKCWw+YLPB7bf74OnpYN069b6KiOQXSmCzr6C2zSI3wy01Eu/Db+J5YiEOozcpFYaTVrIL5phleEZ9hjH1OHZ3f9JCepBaug/WgFvzbGivwXIe05nVzvtmz6zGYNfvpeQdu3sg51ruxOHun+1jKYHNB77+2p3HHvNizpwUunRx4fpiIiKSiRLY7CuobbNIdhiT9uNzcALmU8sAcOCGpXgrUkv3JS24Exi9XBugLQXT2Z8xXpj8SSS32U0lSSt1PxiM2T6WElgXs9vhrru8Afjll2TcNIpDRCTfUAKbfQWxbRbJKe7nt+Aev5X04A5a4kYkh1yrbVYqlQdWrnRn714jTz2VruRVRERcJi0tjZ49e9K1a1c6derE+++/f9W6P/74IzVq1GDHjh15GKFIwWMNaERquUFKXkXyiLurAyjsHA54910TVarYufdeDR0WERHXMZlMLFiwAB8fHywWC3369OGOO+6gQYMGmeolJibyySefUL9+fdcEKiIichXqD8xlq1cb2bnTyJNPpmHM/nBwERGRm2YwGPDxcS7dYbVasVqtGK4wucy0adMYMmQIZrM5r0MUERG5JiWwucjhgKlTzZQvb6dHD/W+ioiI69lsNu69916aN29O8+bNL+tl3bVrF9HR0dx1112uCVBEROQalMDmop9/NrJ1q/PeVw8PV0cjIiICRqORb7/9lvXr1xMeHs7+/fsz9tntdiZPnsyYMWNcGKGIiMjVXXUW4nPnrv3CYsVyI5ysFZSZDh0O6NzZm6goA3/+mYQpZ9b0FRGRHFaUZyH+4IMP8PLyYtCgQQAkJCTQtm3bjGHGp0+fJiAggOnTpxMWFnbV4xSUtllERAqGa7XNV53E6ZZbnGsuOxxw7BgEBTmfx8VB+fJw5EhuhFp4/Pabkb//NjJ5cqqSVxERyRfOnTuHu7s7/v7+pKamsnHjRoYMGZKx38/Pj7/++iuj3K9fP0aPHn3N5FVERCQvXTWBvZigDhkC990HHTs6y99/D8uW5UFkBdzUqSZCQuz06WNxdSgiIiIAnDp1irFjx2Kz2XA4HNxzzz20atWKadOmUbduXdq0aePqEEVERK7pqkOILwoLg/8uAXelbXmlIAxT+vNPI127evP666k8+qgSWBGR/KwoDyHOKQWhbRYRkYLjpoYQX1S6NEyYAA895Cx/+qlzm1zd1KkmSpSw06+fklcREREREZGckuUsxEuWwOnTzmHE3bs7ny9ZkhehFUz//OPGL7+4M2yYBW9vV0cjIiIiIiJSeGQ5hPiipCS4MCmhS+X3YUp9+3rxzz9ubN6chK+vq6MREZGsaAhx9uX3tllERAqWa7XNWfbAbtwItWtDrVrO8vbt8PjjORZboRIe7saaNe48+qhFyauIiIiIiEgOyzKBffpp+PFHKF7cWa5fHzZsyO2wCqapU00EBDgYNCjd1aGIiIiIiIgUOlkmsADlymUuG425EUrBtnu3G6tWeTB4cDr+/q6ORkREREREpPDJchbicuWcw4gNBrBYYNq0S8OJ5ZL33jPh4+Ng6FD1voqIiIiIiOSGLHtgP/4YPvwQTpyAMmVg2zZnWS45cMCNb791Z9CgdIKCXB2NiIiIiIhI4ZRlD2yJEs61X+Xq3nvPhJcXPPaY1n0VERERERHJLVkmsKdPw6xZEBEBVuul7XPn5mJUBUh0tIFvvnHnkUcslChxXSsSiYiIiIiIyE3IMoG9915o2RLatr3xyZs2bNjAxIkTsdvt9OrVi6FDh2baf/LkScaMGUNCQgI2m41nn32WO++888ZO4mILFnhgs8Ejj+jeVxERERERkdyUZQKbnAxvvnnjB7bZbLz22mvMmzePkJAQevbsSevWralatWpGnenTp9OhQwf69OnDwYMHGTp0KOvWrbvxk7lIejp88okHbdvaqFxZva8iIiIiIiK5KctJnDp3hlWrbvzA4eHhVKhQgXLlymEymejUqRNr167NVMdgMJCYmAhAQkICJUuWvPETudDy5e6cPu2mdV9FRERERETywFV7YP38nEvnOBwwaRKYzeDh4SwbDBAff+0Dx8TEEBoamlEOCQkhPDw8U53hw4czaNAgFi1aREpKCvPmzcved5PH5swxUaWKnbvusrk6FBERERERkULvqj2wCQnOJDUhAex2SEm5VM4qeb1eK1eu5L777mPDhg3MnDmT0aNHY7fbc+bguWzrVjf++cfIoEHpuGXZjy0iIiIiIiLZlWXq9c03cP78pXJcHCxblvWBQ0JCiI6OzijHxMQQEhKSqc6XX35Jhw4dAGjYsCFpaWnExsZeV+CuNnu2CR8fBw88oKVzRERERERE8kKWCeyrr0JAwKVyYKBzW1bCwsKIiIjg+PHjpKens3LlSlq3bp2pTqlSpfjjjz8AOHToEGlpaRQrVuyGvgFXOHXKwLffutO7twU/P1dHIyIiIiIiUjRkOQvxlUb0/ns92Kse2N2dl19+mcGDB2Oz2ejRowfVqlVj2rRp1K1blzZt2jB27FhefPFF5s+fj8FgYPLkyRgMhpv5PvLUwoUepKcbNHmTiIiIiIhIHjI4HI5rrv/yyCPOXtcnnnCWP/wQzp2D+fNzP7grsVhsxMUlu+bkgMUCt9ziQ61adj7/PMVlcYiISM4IDtZQmuxyddssIiKFy7Xa5iyHEP/vf2AywQMPQO/e4OkJH32Uo/EVKCtXuhMd7cbgwep9FRERERERyUtZDiFetQomT868belS6NUrt0LK32bP9qBiRTtt2mjpHBERERERkbyUZQ/sG29c37aiYMcONzZtcueRR7R0joiIiIiISF67ag/s9987e19PnICRIy9tj48H9yz7bQun2bNNeHs7ePBBLZ0jIiIiIiKS166aipYuDY0bw/LlcMstl7b7+cG77+ZFaPnLmTMGvv7anQcftGRaVkhERERERETyxlUT2Pr1nV99+oCHR16GlD99+qkHaWkGBg1S76uIiIiIiIgrZDkYOCICxo2D3bshNfXS9sOHczGqfMZqhfnzPWjZ0kqNGldYGFdERERERERyXZZTEQ0cCMOGOe97/flnePhheOihvAgt//j+e3dOnHBj8GD1voqIiIiIiLhKlglsSgq0aQMOB1SoAOPHw8qVeRBZPjJnjgfly9u5+26rq0MREREREREpsrIcQmw2g90O1arBBx9AmTKQmJgXoeUPu3a5sXGjOy+/nIrR6OpoREREREREiq4se2CnTYPkZHj/ffjnH1i4EBYsyIvQ8oe5cz3w8nLQt6+GD4uIiIiIiLhSlj2wTZo4H319Yd683A4nf4mNhS+/9KBnTwtBQa6ORkREREREpGjLMoHdvBkmToSjR52z8V4UHp6bYeUPn37qQUqKls4RERERERHJD7JMYPv2hbffhrAwcMtywHHhYbPB/Pkmmje3Uru2ls4RERERERFxtSwT2OBg6No1L0LJX1avdufYMTdeeSXN1aGIiIiIiIgI15HAvvoqDB7sXErHbL60vXv33AzL9WbP9qBMGTsdOmjpHBERERERkfwgywR23jzYuxcslktDiA2Gwp3A7tvnxq+/uvPCC2m4Z/kOiYiIiIiISF7IMj37+2/Yty8vQsk/Zs/2wGx28NBDmrxJREREREQkv8hyWqbmzWH37rwIJX84fx6WLvWge3crxYs7XB2OiIiIiIiIXJBlD+yff0KDBlCpkvMeWIfDOYS4sC6js2SJB8nJBgYPTnd1KCIiIiIiIvIvWSawP/yQF2HkD3Y7zJ1r4tZbrYSFaekcERERERGR/CTLBLZChbwII39Yu9ZIRIQbzz+vpXNERKTwSUtLo2/fvqSnp2Oz2Wjfvj0jR47MVGfevHksXboUo9FIsWLFmDRpEmXKlHFRxCIiIpkZHA5HgbrR02KxEReXnCvHfuABL/bsceOff5Lw8MiVU4iISD4THOzn6hDyjMPhIDk5GR8fHywWC3369OGFF16gQYMGGXX+/PNP6tevj5eXF4sXL2bTpk2899571zxubrbNNysuNRZPdy883T1dHYqIiNyga7XNWU7iVFQcPGjg55/dGTDAouRVREQKJYPBgI+PDwBWqxWr1YrBYMhUp1mzZnh5eQHQoEEDoqOj8zzO7Eq3pdPuyztp/2UrEi2Jrg5HRERykBLYC+bMMWEyOejXT0vniIhI4WWz2bj33ntp3rw5zZs3p379+let++WXX3LHHXfkYXQ546v9X3A0PoI953bx5LrHKWCDzURE5BqyTGD//BOaNAFfXzCZwGgEf/+8CC3vJCTAZ595cO+9VoKD1ciJiEjhZTQa+fbbb1m/fj3h4eHs37//ivW+/fZbdu7cyeDBg/M4wuyx2W38b+u71C1Rj1dum8B3h5bx3j/vuDosERHJIVkmsMOHw5IlUK0apKTA7NnwxBN5EVre+fxzD5KStHSOiIgUHf7+/jRt2pRff/31sn0bN27k448/Zvr06ZhMJhdEd/NWHfmOg3EHeLLRKB5vMIIe1e5n8qYJrI743tWhiYhIDriuIcRVq4LN5ux9HTiwcC2tY7c7hw/fcouNhg21dI6IiBRe586dIz4+HoDU1FQ2btxI5cqVM9XZvXs3L7/8MtOnT6d48eKuCPOmORwO3vtnCpUDqtC58r0YDAamtvofYcH1GfbTEA7EXrm3WURECo4sl9Hx9ob0dGjQAEaPhlKlnElfYfHLL0YOHXLjo49SXB2KiIhIrjp16hRjx47FZrPhcDi45557aNWqFdOmTaNu3bq0adOGt956i+TkZJ588kkASpUqxccff+ziyK/Pz8fXsuPMdt5r9SFGNyMAXu5ezL/nU+7+8k76f/8gP/RYh785wMWRiojIzcpyGZ2jRyEkxJnEvvsunD/vHEJcpUpehZhZTk/V37evF9u3u7FlSxIFbJSUiIjkgIK0jM5TPzzFe/e8B8C0P6fxZLMnM/YNWDaA+d3muySu/LKMzr3LOnD0fASbHtqOyZi5Uf/z5Ea6L+9Mq3Jt+KTDZxkJroiI5D/ZWkZn2TLw9HRO3PTKKzB1KqxYkZPhuc7hwwZ++snIww9blLyKiEi+t+HohoznC7YvyLQvPCY8r8PJV/6M+oM/Tv7O4w1GXJa8AjQr3ZyJt7/FmqM/8uamiS6IUEREckKWCeyCBZdvmz8/FyJxgXnzTBiN0L+/ls4REZH8z4Hjis8F3v9nCsU9i9O3dv+r1hlQZxAP1erPe1veYfnBb/IwOhERySlXvQd2yRJYvBiOHIGuXS9tT0iAYsXyIrTclZgIixd70LWrlZAQ/RMgIiL5n91hJzYlFrvDnvH8YiJrc9hcHJ3r7DgTzk/HVjPu1pfw8fC5aj2DwcAbd7zD3nN7GLluGFUCq1GnRN08jFRERLLrqgls8+bOCZvOnIFnnrm03c8P6tXLi9By19KlHiQkGBg0SEvniIhIwXA+9Ty3zLwlI2ltNLNRxj4DBleF5XL/2zIVXw8/HgkbkmVds9HMvHsW0e7LO+n/Qx9W9/yZYp4Fa7ZlEZGiLMtJnPKbnJgowuGAO+7wxtMTVq9OxlB023wRkSKvIE3idDTuKBUCK7g6jMu4chKnQ3EHaL64McMbPsVLt7163a/bErOZe5d14NbQZnze5Rvc3bJcmEFERPJItiZx+vNPaNIEfH3BZHKuBevvn6Px5bnISAP79hkZMiRdyauIiBQY931+n6tDyHc+2DoNs9HMo/WfuKHXNQppzNt3vsevJ9bz6sYXcyk6ERHJaVl+3Dh8OHz2GfTqBZs3wyefwP4Cvg542bIO1q1Lok6dQrSgrYiIFHqauCmzk4kn+GLfEvrVHkBJ75I3/PreNfuy4/R2ZoR/RN0S9XigZp9ciFJERHLSdY2XqVoVbDZn7+vAgdCwIbzxRm6HlnsMBqhbV8mriIgULCfiTzDy+5FX3f9+h/fzMBrXm77tf9gddh5vcPX3JCvjm09kz7ndPLv+SaoFVadRSOMcjFBERHJalgmstzekp0ODBjB6tHNiJ7tyPxERkTzn5eHFLaVuueI+QxG7J+ZsylkW7p5Pj+r3U97/5u8L9jB6MOvuBbT/8i4G/vAQq3utJ8Q7JAcjFRGRnJTlPbALFzp7Xz/4AHx84Phx+OqrvAhNRERE/q24V3H6N+h/2VfloMr8FfmXq8PLU7N2TCfFmsLIhqOyfaziXsWZ32Ex59PiGPRDP9JtWqFARCS/yrIHtsKFDzW9vOCVV3I7HBEREbkak9GU8Xxr1FYW71jM0t1LqRRUiR61ergwsryVkB7PnB0z6Vi5C9WL1ciRY9YtEca01h8xZPUAxv36HFPumpYjxxURkZx11QQ2LIxrztAbHp4b4YiIiMjVLOi2gFd/eZUlO5dQwrsED9R5AAcOfu7/s6tDy1Pzd83lfFocTzbKfu/rv91btTs7Tofz/taphJWox4C6g3L0+CIikn1XTWBXrHA+fvih87FfP+fjokXXTmxFREQkd9T6sBYtK7RkRZ8VVC1WFYB3/3zXxVHlrRRrCh9v+4A7y7aiQclGOX78cU1fYtfZHTz/23PULFaLZqWb5/g5RETk5l31HtgKFZxfa9bAW285e2TDwuDNN2H16rwMUURERAC+fuBrSvmWotWCVgxZPoS1h9cWuaV1Ptv7KadTTvHULc/myvGNbkY+bjeHCv4VeeTHfpxIiMyV84iIyM3JchInhwN+//1SeeNGzUIsIiLiCt1qduOznp+x94m9tKrUivf+eo9TSacYtmIYqw8V/k+XLTYLH26dRuOQW2le+vZcO0+AOZAF9ywh1ZrKgB/6kmJNybVziYjIjTE4HI5rfnT7zz/wyCNw/ryzHBgIc+dCo+sYtbNhwwYmTpyI3W6nV69eDB06NNP+SZMm8ddfzlkTU1NTOXv2LJs3b77mMS0WG3FxyVmfXERE5DoEB/u5OoRsiU2JZenupXy+63PWPrzWJTHkVdv8xb4lDF/7KAs7fk77ih1y/Xw/RnxPv1UP0LP6A3zYZmaRW6pIRMRVrtU2Z5nAXnQxgQ0IuL6T2mw22rdvz7x58wgJCaFnz55MnTqVqlWrXrH+woUL2b17N2+88cY1j6sEVkREclJBT2Dzg7xom+0OO3d81hSjwZ2fH/gdN0OWg8hyxJTNb/LmpomMbz6RxxuMyJNziogUdddqm6/7r39AwPUnrwDh4eFUqFCBcuXKYTKZ6NSpE2vXXv2T4ZUrV9K5c+frP4GIiIgUGT8cWcX+2H08ecuoPEteAZ6+5Tk6Ve7K+I0v0P3bzqyJ+AG7Q/dSiYi4Sq61ADExMYSGhmaUQ0JCiImJuWLdEydOEBkZSbNmzXIrHBERESmgHA4H07a8QwX/inStcl+entvN4Mb0trN55bYJHI47RN9V99Nyya18smue7o0VEXGBqyawS5c6H48cyf0gVq5cSfv27TEajbl/MhERESlQNkT+wtZTWxjR8Gnc3a66AmCu8XT35ImGI/n7oXCmt52Nl4c3z65/klsW1uGtTZM4nXw6z2MSESmqrprAXrwVtUePmztwSEgI0dHRGeWYmBhCQkKuWHfVqlV06tTp5k4kIiIihdq0LVMI8Q7lgZp9XBqHh9GDHtXvZ03P9Xxz70puCWnCO5sn02hhbUb9PIL95/a5ND4RkaLgqh9jFi8Od9/t7IHt2vXy/cuXX/vAYWFhREREcPz4cUJCQli5ciVTpky5rN6hQ4eIj4+nYcOGNxy8iIiIFG6bozfx24kNvNp8Emaj2dXhAGAwGGhRpiUtyrTkYOwBPt7+IV/sW8yiPQtoW/5uhjUYwe1l7tCsxSIiueCqsxCnp8OWLdCvH8yeffn+O+/M+uDr169n0qRJ2Gw2evTowbBhw5g2bRp169alTZs2APzvf/8jLS2NZ5+9vgXJNQuxiIjkJM1CnH252TY/vKo3f0X9wT8P78LXwzdXzpETzqScYf7O2czdOYszKaepW6Iew+oP596q3TEZTa4OT0SkQMnWMjqnT0NwMCQmOsu+Lm47lMCKiEhOUgKbfbnVNu85u5s7P2/Gc03G8VyTcTl+/NyQak3lq/1f8PH2D9gXu5dSPqUZFPYo/esMJMAc6OrwREQKhGwlsDt3Onthz50Dh8OZzC5YAHXr5nic10UJrIiI5CQlsNmXW23zsDWD+f7ISrY+vIsgz2I5fvzcZHfY+fnYT3y0/QN+jfwFb3cf+tbqx9D6j1PBv6KrwxMRydeytQ7s0KEwdSocPQrHjsGUKc5tIiIiIrkl4vwRvjn4Jf3rPFLgkldwLr/TpsLdfNV1OWvv/41Olbswb9dsmn7agEE/PsyhuAOuDlFEpEDKMoFNSoJWrS6V77rLuU1EREQkt3ywdRruBneGNRju6lCyLaxEPT5sO5N/HtrJ8AZPsf74z9z9ZSt+Ovqjq0MTESlwskxgK1eG11+HiAjn14QJzm0iIiIiuSE6KYrP9i6id82HCPUp5epwckwp39K8eNt4fnlgIxX9K9F35f1M+2cKWdzNJSIi/5JlAjt3rnMip+7dnWvCnjnj3CYiIiKSGz7e/iFWh5XhDZ90dSi5oqxfOb6770e6Ve3OxL9e5dE1A0myaHibiMj1yHISp/xGkziJiEhO0iRO2ZeTbXNs6jkaflKHeyp15ON2c3LkmPmVw+Hgg23TmPDHK9QpEcaCDosp51fe1WGJiLhctiZxEhEREckrs3fMINmaxMhGo1wdSq4zGAyMaPgUizst5Vj8Ue5eeicbT/zm6rBERPI1JbAiIiKSLyRaEpkVPp17KnakdvE6rg4nz7SpcDc/9lxHMc/i9PyuK3N2zNR9sSIiV6EEVkRERPKFhbvmE5cWVyR6X/+rSmA1fui5jjbl2zHu12cZ9csI0mxprg5LRCTfyfIe2NOnYdYs5wzEVuul7a6ayEn3wIqISE7SPbDZlxNtc5otjcYLw6gWVJ2v712RQ5EVPHaHnbc2TWTqP2/TOORW5t2ziBCfUFeHJSKSp67VNrtn9eJ774WWLaFtWzAaczQuEREREQBOJceQaEnk2cZjXR2KS7kZ3Bjb9CXqlAhjxNrHaPflncy/51MahTR2dWgiIvlClj2wDRrAtm15E8z1UA+siIjkJPXAZl9Otc3ptnRMRlMORFQ47DyzgwHf9yEmOZp37pzGAzX7uDokEZE8ka1ZiDt3hlWrcjQeERERkcsoec2sbokwfuz5C01CmzJi3WO89Ps4rHZr1i8UESnEsuyB9fODpCQwmcDD48KLDBAfnxfhXU49sCIikpPUA5t9aptzl8VmYfzGF5i142PuKNuKWXfPI8izmKvDEhHJNddqm7NMYPMbNZIiIpKTlMBmn9rmvLFkzyKeW/8UpXxL80mHz6hVvLarQxIRyRXZTmCXL4cNG5zP77rLOazYVdRIiohITipKCWxaWhp9+/YlPT0dm81G+/btGTlyZKY66enpjB49ml27dhEYGMi7775L2bJlr3lctc15Z3P0Jgb+8BAJ6Ql80GYGnat0dXVIIiI5Llv3wI4dC9OmQe3azq9p02DcuByNT0RERPKAyWRiwYIFLF++nGXLlvHrr7+y7T8zNS5duhR/f3/WrFnDgAEDeOedd1wTrFxR49BbWdNrPbWK1+KRHx/irU2TsDvsrg5LRCTPZLmMzqpVzlmI3S6kuv37Q8OG8MYbuRyZiIiI5CiDwYCPjw8AVqsVq9WKwWDIVGfdunUMHz4cgPbt2/Paa6/hcDguqyeuE+pTim/uXcXoDU/zzubJrDu2hprFahPqE0qoT2lCfUpdeF6KEl7BuLtl+e+eiEiBcV1/0eLioNiFuQLOn8/FaERERCRX2Ww2unfvzrFjx+jTpw/169fPtD8mJoZSpUoB4O7ujp+fH7GxsRQrpkmD8hNPd0+mtfqIBiUb8dmeRaw7/hOnkmMu6411M7gR7FUyI6kN8b6U3JbyKUWITylCfUpRzLMYboYsB+aJiLhclgnsuHHOHtdWrcDhcN4LO3lyXoQmIiIiOc1oNPLtt98SHx/PE088wf79+6levbqrw5KbYDAYeKTuEB6pOwQAm93G6ZRTRCdFEZ0U7XxMjiLmwvPIhEg2R2/ibOrZy47l4eZBiHcoJb1L4u3hg4ebh/PLaMKU8WjCw+hx4fHSdg83Eyajx4VHU8ZrTUYTwV4lqV28Dv7mgLx+e0SkkMoygX3wQefETX//7Sy/+SaEhuZyVCIiIpKr/P39adq0Kb/++mumBDYkJISoqChCQ0OxWq0kJCQQFBTkwkjlehndjBd6Wktds16aLY1TyTEZiW7MxYT3QrKbZksjxZpMus2CxZ5Oui0di93i/LKlk37hMc2WhoPrW8yivH9F6hSvS53idalboh51StSlvF8FDU0XkRt21QR2716oWRO2bHGWL05AePKk86tRo7wIT0RERHLKuXPncHd3x9/fn9TUVDZu3MiQIUMy1WndujXffPMNDRs25Mcff6RZs2ZKMgoZs9FMOb/ylPMrn+1j2ew20u3pGYmt1W65kPCmk2ZL52RiJLvO7GTnmR3sOruDH46szEh6/Uz+1C5eh7olwqhTPIw6xetSs3htvNy9sh2XiBReV11GZ+hQmDnTOXT4shcZYN263A7tyjRVv4iI5KSitIzO3r17GTt2LDabDYfDwT333MPw4cOZNm0adevWpU2bNqSlpfHcc8+xZ88eAgICePfddylXrtw1j6u2Wa5XkiWJved2s+vMTnad3cHOMzvYfXYXSZZEwHnPbtXAas7e2hJhGT22Jb1DrvuDFIfDgcVuId2WRvqFHuQ0WxrptvQLX2mk2y0U9yxOef8KmIym3PyWReQmZGsd2NRU8PTMelteUSMpIiI5qSglsLlFbbNkh91h52h8hLOn9mw4u8/sZNfZnRxPOJZRp4RXCaoGVseB48LwZWcimnax99eWlrEt3Z5+3ec2GoyU969AlYCqVAmsSuVA52OVgKqU8i2tia1EXCRbCWyjRpeGEV9rW15RIykiIjlJCWz2qW2W3BCXGsvus7syemqPnD+M0WDEZDRhNpoxGc14uHlceO7c5mF0TiRldnPuNxudE05drGNyM2dMOHU65RSH4w5yKO4Qh84f5HDcQZKtl36Ovdy9qBRQhcoBVZxJbWBVKl9IdIt5FtPQepFcdK22+ar3wEZHw4kTkJICW7c6ZyAGiI+HZLVRIiIiIpKLAj2DaF7mdpqXuT1PzudwOIhOiuLQ+YMcinN+HY47yJ5zu/ghYiVWu/VSbObATAltad8yOdZba8BAWb9yVAuqQQmvEjlyTJHC5Ko9sAsWwPz5sHkzNGlyKYH194f+/aF79zyM8l/0Ka+IiOQk9cBmn9pmKewsNgvHE49d6LG98HX+EIfjDnIiMTLXzlvcszjVgmpQLagG1YOqX3isQRnfsuoBlkItW0OIv/oKevTI8ZhumhpJERHJSUpgs09tsxRlyZZkYpKjc+x4doeNo/FH2R+7lwOx+9kfu4/95/YSmxabUcfHw5dqgdUyEtqLjxUDKuHuluUqmSL5XrYS2Oefh9GjITDQWY6NhSlTYMKEnAzx+qmRFBGRnKQENvvUNovkLofDwZmUMxyI3cf+2H0Zj/tj9xGVdDKjnsnNROXAKhd6batTPagGpXxKY0C9tZL7SngFUzWoWo4cK1sJbMOGzntg/02TOImISGGhBDb71DaLuE5Cevylntp/JbdH4yOwO+yuDk+KEAMG9g86SoA5MNvHuqlJnC6y2SAtDcxmZzklxVkWERERERHX8jP50yikMY1CGmfanmpN5VDcQc6mnnFRZFLUBHuVzJHkNStZJrB9+0KbNjBwoLM8b55zEicREREREcmfPN09qVOirqvDEMlxWQ4hBvj+e1i71vm8XTto3z63w7o6DVMSEZGcpCHE2ae2WUREclK27oHNb9RIiohITlICm31qm0VEJCddq23OcsXlP/90rgPr6wsmExiNzrVgRURERERERPJSlgns8OGwZAlUq+acwGn2bHjiibwITUREREREROSSLBNYgKpVnbMRG43OyZx++CG3wxIRERERERHJLMtZiL29IT0dGjSA0aOhVCmwa0kpERERERERyWNZ9sAuXOhMWD/4AHx84Phx+OqrvAhNRERERERE5JJrzkJss8HDD8Onn+ZlSNemmQ5FRCQnaRbi7FPbLCIiOemmZyE2GuHoUecQYhERERERERFXyvIe2MqVoUUL6NrVOYT4olGjcjMsERERERERkcyyTGCrVHF+2e2QkJAXIYmIiIiIiIhc7qoJbL9+zgmcAgPhySfzMCIRERERERGRK7jqPbD//AMnT8LcuRAbC+fOZf66Hhs2bKB9+/a0a9eOmTNnXrHOqlWr6NixI506deKZZ565qW9CRERERERECr+r9sA+9hi0aQOHD8Mtt8C/5yo2GJzbr8Vms/Haa68xb948QkJC6NmzJ61bt6Zq1aoZdSIiIpg5cyZLliwhICCAs2fPZvsbEhERERERkcLpqj2wI0fCnj3wyCPOZPXIkUtfWSWvAOHh4VSoUIFy5cphMpno1KkTa9euzVTniy++oG/fvgQEBABQvHjx7H03IiIiIiIiUmhdcxkdgOnTb+7AMTExhIaGZpRDQkKIiYnJVCciIoIjR47Qu3dv7r//fjZs2HBzJxMREREREZFCL8tZiHOTzWbj6NGjLFy4kOjoaB566CG+++47/P39XRmWiIiIiIiI5ENZ9sDerJCQEKKjozPKMTExhISEXFandevWeHh4UK5cOSpWrEhERERuhSQiIiIiIiIFWK4lsGFhYURERHD8+HHS09NZuXIlrVu3zlSnbdu2bNq0CYBz584RERFBuXLlciskERERERERKcBybQixu7s7L7/8MoMHD8Zms9GjRw+qVavGtGnTqFu3Lm3atKFly5b8/vvvdOzYEaPRyOjRowkKCsqtkERERERERKQAMzgc/14gJ/+zWGzExSW7OgwRESkkgoP9XB1Cgae2WUREctK12uZcG0IsIiIiIiIikpOUwIqIiIiIiEiBoARWRERERERECgQlsCIiIiIiIlIgKIEVERERERGRAkEJrIiIiIiIiBQISmBFRERERESkQFACKyIiIiIiIgWCElgREREREREpEJTAioiIiIiISIGgBFZEREREREQKBHdXByAiIiJ5IyoqitGjR3P27FkMBgP3338//fv3z1QnISGB5557jpMnT2Kz2XjkkUfo0aOHiyIWERHJTAmsiIhIEWE0Ghk7dix16tQhMTGRHj160KJFC6pWrZpR59NPP6VKlSp8/PHHnDt3jnvuuYcuXbpgMplcGLmIiIiThhCLiIgUESVLlqROnToA+Pr6UrlyZWJiYjLVMRgMJCUl4XA4SEpKIiAgAHd3fd4tIiL5g1okERGRIigyMpI9e/ZQv379TNv79u3LsGHDaNmyJUlJSbz77ru4uenzbhERyR/UIomIiBQxSUlJjBw5kueffx5fX99M+3777Tdq1arFr7/+yrJly3jttddITEx0UaQiIiKZKYEVEREpQiwWCyNHjqRLly7cfffdl+3/+uuvufvuuzEYDFSoUIGyZcty+PBhF0QqIiJyOSWwIiIiRYTD4eCFF16gcuXKDBw48Ip1SpUqxR9//AHAmTNnOHLkCGXLls3LMEVERK7K4HA4HK4O4kZYLDbi4pJdHYaIiBQSwcF+rg4hz2zevJm+fftSvXr1jPtaR40axcmTJwF48MEHiYmJYdy4cZw+fRqHw8GQIUO49957r3lctc0iIpKTrtU2K4EVEZEirSglsLlFbbOIiOSka7XNGkIsIiIiIiIiBYISWBERERERESkQlMCKiIiIiIhIgaAEVkRERERERAoEd1cHICIiIiIi18dmsxIbexqrNd3VoYhkm7u7iaCgYIzG609LlcCKiIiIiBQQsbGn8fT0xscnFIPB4OpwRG6aw+EgKSme2NjTlChR6rpfpyHEIiIiIiIFhNWajo+Pv5JXKfAMBgM+Pv43PJpACayIiIiISAGi5FUKi5v5WVYCKyIiIiIieeLAgX388cdvWdZr167lFbcvW/Yl33+/IsfiefHF0Zw4EXlddSdPfp0jRw4D8MknczO2R0WdpF+/+3MspqvZsmUzo0c/levnyQkbNvyS8V4BfPDBe/zzz985cmwlsCIiIiIikicOHNjPH3/8ftOv79atJx06dM6RWA4fPoTNZqdMmbLXVX/s2JeoVKkyAAsXzsvWua1Wa7Zen9/9+usvRERcSmB79nyARYvm58ixNYmTiIiIiIhcl6iokzzzzAjq1Aljx45watWqTceOXZg7dwaxsbG8/PLr1K5dl927dzJt2hTS09Mwmz15/vmXKVWqDLNnf0x6ehrh4dvp128At912O++99zZ79+7GYDAwcOAQ7rqrDQAzZnzIxo2/YTabmTx5CsWKFWfOnBl4eXnTp08/hg8fSu3addm6dTMJCYmMG/cS9es3JDU1lYkTx3PkyCHKlavAmTOneeaZMdSsWTvT97JmzQ+0bHknAOvW/cSuXeGMGDGKL75YwtKln7F06becOBHJhAkvM336XIYPH8rw4U/x889rSUtLY8CAPlSqVJmhQx/Hbrfz5psT2LEjnODgYCZPnoLZ7JnpfBMnjsdkMrF//z7q1atP9+73M2XKm8TFxeLp6cmYMS9SoULFjHp79+4hKSmJESOepkWLzD3SV3p/y5eviM1mY/r0//HXXxtxc3OjS5du9OzZm7179/DBB++SnJxMYGAgzz8/nhIlSjB8+FCqV6/B9u3bSE1N4cUXX2XhwvkcPnyQ1q3bMXTo4wD8+OMqvvzyMywWK7Vr1+GZZ8ZiNBpp164lPXv2znSdTpyI5LffNrBt2xYWLJjLxIlvUaZMWc6fP8/Zs2coXrxEtn4G1QMrIiIiIiLX7cSJSHr3fojFi7/k6NEI1qz5gY8+msMTTzyZ0TNZoUJFPvxwFvPmLWbQoEeZMeNDPDw8GDz4MVq3bsf8+Ytp0+Zu5s+fjY+PL5988jkLFnxGo0ZNAEhJSaFOnTAWLFhCgwYNWb78myvGYrPZmDXrE558chRz584C4Ouvl+Ln58eiRUsZMuQx9u/fe8XX7tixnRo1agFQv34Dtm/fBkB4+FYCAgI4ffoU4eHbqF+/UabXDRs2ArPZzPz5i3nllQkAREYep3v3Xixa9AW+vn788su6K57z9OlTfPzxXEaMGMVbb03k6aefY+7cRTzxxFNMmTI5o15UVBSzZi3g7bff45133iAtLS3Tca70/gIsX/4N0dEnmTdvMQsWfMbdd3fAarXy3ntv8/rrbzJ37iI6derKzJkfZhzL3d2DOXMWcu+9PRg79hlGjRrDJ598zvffr+D8+TgiIo6wdu0apk+fy/z5i3FzM7J69fdXvU5hYfW5/fY7ePzxkcyfvzijh7tGjZrs2LH9iu/LjVAPrIiIiIhIAfT53sUs2bsoR4/5YM2HeKBmn2vWKVWqNFWqVAWgUqXKNG58KwaDgcqVqxIVFQVAYmIiEyaMJzLyGAaD4apDZjdv3sSrr07KKPv7+wPg4eGR0etYo0Yt/v77ryu+/s47W2XUiY4+CcCOHdvo1etBACpXrpoR63+dOXOGwMBAAIoXL0FKSjLJyUnExMTQrl17tm3bwvbtWzPOkdV7Uq1ajQux1CQq6uQV67Vq1Raj0UhycjI7doTz0ktjM/ZZLJdm423dui1ubm6UK1ee0qXLcOxYRKbjXO393bz5L7p164G7uzPN8/cP4PDhgxw+fIinn34CALvdlqkX9Pbb7wCgSpWqVKpUmRIlnPtKly7DqVMxhIdvY9++PQwe/DAAaWmpBAUFAdd/nQACA4M4c+b0td7G66IEVkRERERErpuHh0fGczc3t4yym5sbNpszkZo9+2MaNWrMG2+8Q1TUSUaMePSGzuHu7p4xQ63zuLYr1jOZTBfqGK9a52rMZjPp6ZeSxrp167Fy5XeUL1+B+vUbsmLFcnbu3MHw4U9neazM74kRmy3tivU8PZ3Dih0OO35+vsyfv/iK9S6fnTdz+UbeX4fD+UHDjBlXvm/34ntoMBgynl8s22w2HA4HHTp05rHHhl/22uu9TgDp6emXDau+GUpgRUREREQKoAdq9smyt9RVEhMTCQ4OBmDVqu8ytnt7e5OcnJxRbtKkKV9/vZQnn3wGgPj4+Ixe2JsVFlafdevW0KhRY44cOcyhQwevWK9ixYpERh6nVKnSANSv35DZsz9m4MAhVKtWg61bN2M2m/H19b3stUajO1arNaOn80b5+PhSqlQZ1q37idat2+JwODh48ADVqlUH4Oeff6JDh85ERZ3k5MkTlC9fgV27dmS8/mrvb5MmTfn2269p2LAx7u7uxMefp3z5CsTFxbJzZzh169bDarVy7NhRKleucl2x3nLLrYwb9wwPPNCHoKBixMefJzk5mdDQUld9zX+vM8Dx48do1artdb9HV6N7YEVEREREJEf17fswH3/8IQMH9snUK9eoUWMiIo4wYEAf1q5dTf/+g0hIiKdfv/vp3/9Btm7dnO1z33dfL+LiYnnooV7MmjWdSpWq4ONzeRJ62223s3XrPxnl+vUbcupUDPXrN8RoNFKyZAj16jW44jm6dr2P/v178+qrL950nC+//DorVnxL//4P0q/f/fz22/qMfSEhoQwZ0p9nnhnJs8+Ow2w2Z3rt1d7fzp27ERISyoABD9K//4OsWfMDHh4eTJjwJtOn/4/+/R9kwIA+7NwZft1xVqpUmSFDhvH008Pp3783Tz31BGfOnLnma9q0uZslSxYycGAfTpyIxGq1cuLEcWrWrHXd570ag8PhcGT7KHnIYrERF5ecdUUREZHrEBzs5+oQCjy1zSJ5Jzr6KKGhFVwdRr5ms9mwWq2YzWZOnIjkqaceZ/HirzIN8wXnvZwjRjzG9OlzMBqNLor2chMnjqd589tzpLcyv1i//mf279/LkCHDLtt3pZ/pa7XNGkIsIiIiIiKFxsXE1DmxkYNRo8ZclrwCmM2eDBr0KKdPnyY0NDTvAy1CbDYbvXs/lCPHUg+siIgUaeqBzT61zSJ5Rz2wUtjcaA+s7oEVERERERGRAkFDiEXkihwOB5YkC2mxqaTGppAam0paXOplj2mxqaTGOZ9bk6+8xptITvMK9qLLFz0w+ZmzriwiIiKFhhJYkULOYXeQFp92KdH8byL67yT0P3XsVvtVj+vu7Y450BPPQE/MQZ4EVgnC3cvjv8uUieQKrxLeuHnknwk3REREJG8ogRUpQNIT0kg+lUxqbMoVe0MzekUzPU/DYb/6re4eviY8gzwzktHitUo4n1/cVuxSknrx0Rzgibun/nyIiIiISN7Sf6Ai+YDdZifldDKJUYkkXfhKjEpwPo++tM2SZLnyAQxgDjBnSjz9KwT8q4fUK3NSeuHRHGjGqF4sERERySMHDuzjzJnT3Hbb7des165dS9as+fWy7cuWfYnZ7EmHDp1zJJ4XXxzNsGEjKVOmbI4c72ZFRZ1k9OinWLjwC5fGcT22bNmMh4cHYWH1Afjqq88xmz3p3PnePDm/EliRXGZJtpB4MiEjCb2YkCb+63nyqSQctsy9pG7ubniH+uAT6kvxWiUo37oiPqG+eIf44FnM61KvaJAnJn8zbkbNySYiIiL524ED+9m7d3eWCezVdOvWM8diOXz4EDabPceTV6vVirt74U2ztm79By8v74wEtlOnexk27JHCkcBu2LCBiRMnYrfb6dWrF0OHDs20/+uvv+att94iJCQEgIceeohevXrlZkgiuSI9MZ34I3HEHYnj/JE4zh+JJf7C86TopMvqm/zN+JTyxSfUl2I1imc89y3lm/Hcq4Q3BjfdUCoiIiL5R1TUSZ55ZgR16oSxY0c4tWrVpmPHLsydO4PY2Fhefvl1ateuy+7dO5k2bQrp6WmYzZ48//zLlCpVhtmzPyY9PY3w8O306zeA2267nffee5u9e3djMBgYOHAId93VBoAZMz5k48bfMJvNTJ48hWLFijNnzgy8vLzp06cfw4cPpXbtumzdupmEhETGjXuJ+vUbkpqaysSJ4zly5BDlylXgzJnTPPPMGGrWrJ3pe1mz5gdatrwzo9yuXUt69ux92Tmjok7yxhuvcf58HIGBQYwb98pl68bOmTODkycjOXnyBCVLhvLUU8/xzjuTiImJAWDkyFHUq9cgo15kZCTnz8fRp8/DdO1632Xv8euvv0xqagoATz89OiNZXLRoPqtXf4/B4EazZs0ZNmwEJ05EMmXKm8TFxeLp6cmYMS9SoUJFJk4cj9lsZv/+fcTGxjJu3Ev88MNKdu3aQe3adXnhhfEAbNr0J3PmzMBiSad06bI8//wreHt707NnFzp06Mzvv2/AarXy+utvYjKZ+Pbbr3Fzc2P16u95+unnqF+/IaGhpdm9eye1a9fNuR+2q8i1BNZms/Haa68xb948QkJC6NmzJ61bt6Zq1aqZ6nXs2JGXX345t8IQyTFp51MvJKeXf6Wczrz+oXdJHwIqBVLuror4VwzAr5z/heTUD59QXzx8Ll9MW0RERKQgOHEiktdff5Nx4yozePDDrFnzAx99NIffflvPwoXzeOONKVSoUJEPP5yFu7s7f//9FzNmfMjEiW8zePBj7N27m1GjxgDw0Ufv4+PjyyeffA5AfHw8ACkpKdSpE8ajjz7BRx9NY/nybxgwYPBlsdhsNmbN+oQ//viNuXNnMW3aR3z99VL8/PxYtGgphw8fZODAvlf8Pnbs2E7btu0zylc757vvvk2HDp3p0KEzK1Z8y7Rpb/PGG1MuO96RI0eYPn02ZrMn48e/wP3396V+/QZER0fzzDPD+fTTLwE4ePAgM2fOIyUllUce6Uvz5pl7o4OCivHuux9iNps5fvwY48e/wJw5C/njj9/57bcNzJy5AE9PT+LjzwPw1lsTefbZcZQrV55du3YyZcpk3n//YwASEuKZMWMev/22nrFjn2H69DlUquS8bgcO7CM4OIQFC+bw3nsf4eXlxaJF8/n8808ZOHAIAAEBAcyd+ylff72UJUsWMnbsS9x7b/eMDxEuqlmzFtu3byvYCWx4eDgVKlSgXLlyAHTq1Im1a9delsCK5CWHw4E93YY11er8SrFiu/j8X+X0hHTij57/V5IaS+q51EzH8inlS0ClQCreXZmASoEXvoIIqBiAh6/JRd+hiIiIFBXmk4vxPLkoR4+ZWvoh0kr3uWadUqVKU6WK83/6SpUq07jxrRgMBipXrkpUVBQAiYmJTJgwnsjIYxgMBqzWKy+1t3nzJl59dVJG2d/fHwAPDw9atGgJQI0atfj777+u+Po772yVUSc6+iQAO3Zso1evBwGoXLlqRqz/debMGQIDAzPKVzvnrl3hTJr0NgD33NOJ6dPfv+Lxbr/9Dsxmz4zvKyLiSMa+pKQkkpOdHR4tW96J2eyJ2exJw4a3sHv3LqpVq55R12q18u67b3LgwH7c3IwcP34045gdO3bB09PzwnsVQHJyMjt2hPPSS2MzXm+xpGc8b9HijoxrU6xYsUzXLSoqilOnThERcZhhwwZdOLeFOnXC/vX+ts54P9av//mK3zc4k+6jRyOuuj8n5VoCGxMTk6lrPSQkhPDw8MvqrV69mr///ptKlSoxbtw4SpUqlVshSQHlcDiwJKZfccbd1NhLS76kJ6RlJKK2C8lopvKFbVx9Qt7MDOBbxo+ASoFU7lTNmaBWDiKgUiD+FQLw8FYvqoiIiBQ9Hh6X/gdyc3PLKLu5uWGzORPV2bM/plGjxrzxxjtERZ1kxIhHb+gc7u7uGAyGfx3XdsV6JpPpQh3jVetcjdlsJj39UrJ3vee8Gk9Pr4znDoedGTPmYTZfvl75xXNcKmfe//nnnxIUVJz585dgt9tp06bFVc/pcNjx8/Nl/vzFV9z/72vz3+tms1lxc3OjceOmmT5EyPx65/trNF66tlfiHCqeN2uzu/Tu4latWtG5c2dMJhOfffYZY8aM4ZNPPnFlSJLLHA4HaefTMiYzSou98nqkqedSMi0H898Jjv7t4nqkJj8z7p7uuHu64+FrxivYJ6NsvPDo7vWv51cqe7pj9HLHw8cDv7L+WipGRERE8q200n2y7C11lcTERIKDgwFYteq7jO3e3t4ZPZEATZo05euvl/Lkk88AziHEF3thb1ZYWH3WrVtDo0aNOXLkMIcOHbxivYoVKxIZeZxSpUpf83h169bjp59+5J57OrF69ffUq9cwyxiaNGnGV199Tp8+DwPO2ZerVasBwK+/ruehhwaQmprC1q3/MGzYCCyWSytNJCUlEhwcgpubG99/vyIjkW7SpCnz58/m7rs7ZAwh9vcPoFSpMqxb9xOtW7fF4XBw8OCBTD2611KnThhTp75JZORxypYtR0pKCqdPn6J8+QpXfY23tw/JyZnneDl+/FjGfbq5Ldf+Ow8JCSE6OjqjHBMTkzFZ00VBQUEZz3v16sXbb7+dW+FIHrBZbCSfSibpP8u//HtpmKSYRKzJV/70xuRnurTWaKAnvmX8Llt/VOuRioiIiOR/ffs+zIQJ41mwYE6mGYcbNWrMokULGDCgD/36DaB//0FMnfom/frdj5ubkUceGZIxbPVm3XdfLyZOfIWHHupF+fIVqVSpCj4+vpfVu+2229m69R+aNGl6zeM9/fRoJk16lSVLFmZM4pSVp556jqlT36R//97YbDbq12/Ic889D0CVKlUZOfIxzp+PY8CAwZQoEUxU1MlM8b/44mh++GElTZvehpeXs2e3WbPmHDiwn8GD++Hu7sFtt7Xg0Uef4OWXX+eddyazYMEcbDYrbdrcfd0JbFBQEC+8MJ7x41/IGHo8ZMiwayawLVq05KWXxvDrr+szJnHasWM7jzwy9KqvyUkGh8NxvQMqb4jVaqV9+/bMnz8/YxKnKVOmUK1atYw6p06domTJkgCsWbOGWbNm8cUX1177yGKxEReXfM06kjuSohM5t+9s5qT0X2uUJp9Oumx4rpvJiG/ohZl1L8yue+m5D17FvbUeqYi4VHCwn6tDKPDUNovknejoo4SGXj25EOfETlarFbPZzIkTkTz11OMsXvxVpiG0AGlpqYwY8RjTp8/BaMyb/0P/PYtyYbF//14+//xTXnrp9Zt6/ZV+pq/VNuda15W7uzsvv/wygwcPxmaz0aNHD6pVq8a0adOoW7cubdq0YeHChaxbtw6j0UhAQABvvPFGboUjNyg9IY3T22OI2RLNqa3RxGyJJikqMVMdc5BnxtIvJeoGZ0pOfUN98Snth2cxz8vG+YuIiIiI5JaLialz4igHo0aNuSx5BTCbPRk06FFOnz592bI4cv3On49j8OBheXa+XOuBzS36lDfn2Sw2zu0540xWt0QTszWa2P1nM3pT/SsGENKoFCUbhVKiTnDGWqXuXprESEQKPvXAZp/aZpG8ox5YKWzyTQ+s5E8Oh4P4iPPEbInK6Fk9s/MUtlTnzeGexb0IaRRK1XurE9IolJINQvEs5pXFUUUKNpvNSmzsaazW9KwrS4Hl7m4iKCgYo1FNn4iISEGlVryQSzmb4kxWLw4F3hpNWqxzPVN3L3eC64VQd0ADSjYKJaRRKH7l/DXkV4qc2NjTeHp64+MTqp//QsrhcJCUFE9s7GlKlNBybSIiIgWVEthCxm61E7MlmuPrjnB0XQSnt8eAAwxuBorVLE7lTlUJaRhKyYalKFazOG7ubq4OWcTlrNZ0Ja+FnMFgwMfHn8TEOFeHIiIiItmgBLYQSIpO5Ni6CI6tiyBy/VHSzqdhcDMQ0rgUt45uTunmZQkOK4mHr8nVoYrkW0peCz9dYxERkYJPCWwBZEu3Eb3pREbSenb3GQB8Qn2o1Kkq5VtXouwd5fEM9HRxpCKSU559diSvvDIRgDVrfqB7914AbNmymc8+W8Rbb72Xq+dfteo79u7dzahRY3L1PDlh1arvuPXWZpQoEezqUERECqXHHnuEjz+e6+owpIhSAltAxB87z7F1ERxfF0Hkr8ewJFlw83CjVNMy3PZSS8q1rkjx2iXUwyBSSL3zzvsAREWd5JtvlmYksDfDarXi7l54//yvWvUdlStXUQJ7BVFRUYwePZqzZ89iMBi4//776d+//2X1/vrrLyZNmoTVaiUoKIhFixa5IFoRya+UvIorFd7/YAo4a4qFk3+c4NjPzl7WuAPnAPAr50/1nrUo36YSZW4vh0nDgkUKvMWLP8HDw0SvXr15//0pHDx4gPff/5h//vmbFSu+5ZVXJtCzZxdmz17Ixx//jxMnTjBgQB+aNGnKbbe1IDk5mRdfHM3hw4eoUaMWL7/8+mUfZg0fPpRq1WoQHr6Ntm3b07DhLXzwwbskJycTGBjI88+Pp0SJEgwfPpSqVauzbdsWbDYr48a9TO3adTMd67ffNrBgwRysVgv+/oG88srrFCtWnOTkZN5772327t2NwWBg4MAh3HVXGzZt+pM5c2ZgsaRTunRZnn/+Fby9venZswtt27bnzz83YjQaGT36BWbM+IDIyOP06dOPbt16Zrw/69b9hMWSzh13tGLQoEeJijrJs8+OpF69BuzYEU5wcDCTJ09h48bf2LdvD6+++iJmsyczZszFbNZolIuMRiNjx46lTp06JCYm0qNHD1q0aEHVqlUz6sTHx/Pqq68ye/ZsSpcuzdmzZ10YsYjkR+3atWTNml/ZsmUzc+bMwNvbm8jI4zRq1JhnnhnLL7+sY9eucEaMGMUXXyxh6dLPWLr0W06ciGTChJeZPl0JsNw8JbD5SNr5VA6tOMDh7w5wYuNxbKk2jJ5GyjQvR53+9SjfuiKBVYLUyypSyNSr15DPPltEr1692bt3DxZLOlarle3bt1K/fsNMdR97bASHDx9i/vzFgHMI8YED+1i48AtKlAhm2LBBhIdvp379Bpedx2KxMGfOQqxWK8OHD+WNN6YQFBTE2rWrmTnzQ55//hXAuQD8/PmL2bZtC2+88RoLF37xn3gbMHPmfAwGA999t4xPP/2EESOeZv782fj4+PLJJ58DzkQoLi6OBQvm8N57H+Hl5cWiRfP5/PNPGThwCAAhIaHMn7+Y99+fwqRJ45k+fQ5paek8/PADdOvWk02b/uT48ePMmrUAh8PB2LGj2LZtCyEhoURGHmf8+ImMGfMiL73k/IepffuOfPXVFwwf/hQ1a9bO6UtV4JUsWZKSJUsC4OvrS+XKlYmJicmUwH733Xe0a9eO0qVLA1C8eHGXxCoiBcOePbtYuPALQkNL8cwzI1i/fh316zdg8eJPAAgP30pAQACnT58iPHwb9es3cnHEUtApgXUxa4qFoz8dYf9Xezn60xHs6Tb8KwZQ52Fnwlr6trK4e3m4OkyRIuPzz91ZsiRnf+cefNDCAw9Yr7q/Zs1a7Nu3l6SkRDw8TFSvXpO9e3ezffs2nnrq2SyPX6tWHUqWDAGgWrXqREefvGIC26ZNOwCOHYvg8OFDPP30EwDY7TaKFy+RUa9t2/YANGjQiKSkJBISEjId5/TpU7zyyjjOnj2DxWKhVKkyAGzevIlXX52UUc/f35/ff/+ViIjDDBs2CACr1UKdOmEZdW6//Q4AKleuSkpKCt7ePnh7++Dh4UFCQgKbNv3J33//ycCBfQFISUkmMvIYISGhlCpVmmrVagBQo0ZNoqJOZvleySWRkZHs2bOH+vXrZ9oeERGB1WqlX79+JCUl8fDDD9OtWzfXBCkiWeq2rONl27pWvY9H6g4h2ZJMn5U9L9vfu2Zfetfsy9mUswz6sV+mfcu6rbqh89eqVYcyZcoCzvYjPHw7rVq1JSUlmeTkJGJiYmjXrj3btm1h+/at3Hlnqxs6vsh/KYF1AbvVzonfjrP/qz0cXnkQS2I63iE+hA2sT7UeNQmuH6JeVpEixN3dndKlS7Nq1QrCwupRpUpVtmzZzIkTx6lYsVKWrzeZLt1K4Obmhs1mu2I9Ly8vABwOqFSpMjNmzLtivf/+/flv+d1336J3777cfvudbNmymblzZ141NofDQePGTTMltv/m4WHKiNvD49IHBxe/D4fDwUMPDaBbtx6ZXhcVdfI/9Y3YbGlXjUMyS0pKYuTIkTz//PP4+vpm2mez2di1axfz588nNTWV3r17U79+fSpVyvpnUUSKnsvbDOdj3br1WLnyO8qXr0D9+g1ZsWI5O3fuYPjwp10QpRQmSmDziMPh4NSWaPZ/vZeDy/aRcjoZk5+JKl2rUb17LUq3KIubUWuyirjaAw9Yr9lbmlvq1WvAkiULGTfuZapUqcr//vcuNWrUuuwfA29vb5KTk7N1rvLlKxAXF8vOneHUrVsPq9XKsWNHqVy5CgBr166mUaPGbN++DV9f38sSnKSkREqUcA5D/eGHlRnbmzRpytdfL+XJJ58BnEOI69QJY+rUN4mMPE7ZsuVISUnh9OlTlC9f4bpibdr0NmbNms7dd3fA29ub06dPZTkBlbe3T7bfo8LMYrEwcuRIunTpwt13333Z/tDQUAIDA/H29sbb25vGjRuzd+9eJbAi+dS1eky9Pbyvub+4V/Eb7nH9r927d3Hy5AlCQ0uxbt0auna9D4D69Rsye/bHDBw4hGrVarB162bMZvNlbYrIjVICm8tiD5zjwFd72P/1XuIjzmM0G6nQrjLVutekQttKuHvqEoiIs6H/5JO51K1bDy8vL0wm8xWHAQcEBBIWVp9+/e6nWbMW3HZbixs+l4eHBxMmvMl7771DYmIiNpuN++9/MCOBNZnMDBzYB6vVOYnTfz3yyFBeemksfn5+3HJLE06ePAFA//6DmDr1Tfr1ux83NyOPPDKEO+9szQsvjGf8+BewWNIBGDJk2HUnsLfe2oyIiCM89thAALy8vHn55ddxc7v6B34dO3bm7bcnaRKnK3A4HLzwwgtUrlyZgQMHXrFOmzZteO2117BarVgsFsLDwxkwYEDeBioiBUatWrV59923MiZxuuMO5xDh+vUbcupUDPXrN8RoNFKyZAgVKlR0bbBSKBgcDofD1UHcCIvFRlxc/v5kPTEqgYPf7OPA13s5HX4Kg5uBMreXo1qPmlTuVA2zv9nVIYrIv0RHHyU09PoSqsJu+PChhXoCpCtd6+BgPxdFk/c2b95M3759qV69esaHAKNGjeLkSef9ww8++CAAs2fP5uuvv8bNzY2ePXtmmcAWhLZZpLDIT21WXq1FLoXbjbbN6v7LIalxqRxecYADX+/lxO/HwQElG4bQ4vW7qNqtOj4hGi4hIiKu1bhxY/bt25dlvcGDBzN48OA8iEhEROTGKIHNBofDwcmNkeycu40jPx7Gnm4jsEoQTZ69jWo9ahJYOcjVIYqI3JAPPrj6hEwiIiL/1qhRYxo1auzqMKSIUQJ7EyyJ6ez/ai875m7l3J6zmIM8nTMI96xFcL2SmkFYREREREQkFyiBvQFxh2PZOW87e5fsIj0+jRJhJWk17W6qdauhtVpFRERERERymRLYLDjsDo6tO8KOOds4tjYCN3c3qnStRt1HGhLapJR6W0VERERERPKIEtirSDufyt4lu9gxdxvxEefxLulDk+duo/bDYZqQSURERERExAWUwP7H2d2n2TFnG/u/2oM12UroraVp+vztVO5YFaPJ6OrwRERyXEFaBmHDhl8oV648lSpVdnUoIiJynZYt+xKz2ZMOHTq7OhQpBJTAAnarnSPfH2THnG2c3BiJ0dNIte41CRvUkOCwkq4OT0TkulmtVtzdC++f9l9//YXmzW9XAisiUoB069bT1SFIIVJ4/8u5Dsmnk9m9aAe7Fmwn6WQifuX9ue3lltTqUxfPYl6uDk9EioioqJM8++xI6tVrwI4d4QQHBzN58hTMZk8OHNjH22+/QVpaKqVLl2XcuJfx9/fP9PqJE8djMpnYv38f9erVp3v3+5ky5U3i4mLx9PRkzJgXqVChYka9vXv3kJSUxIgRT9OiRctMx9q9eyfTpk0hPT0Ns9mT559/mfLlK2Kz2Zg+/X/89ddG3Nzc6NKlGz179mbv3j188MG7JCcnExgYyPPPj6dEiRIMHz6U6tVrsH37NlJTU3jxxVdZuHA+hw8fpHXrdgwd+jgAP/64ii+//AyLxUrt2nV45pmxGI1G2rVrSc+evdm48TfMZjOTJ0/hxIlIfvttA9u2bWHBgrlMnPgWZcqUzbPrJCIiTldrt3788XuWL/8Gi8VC2bJleeml1/H09GTOnBl4eXnTokVLJkx4mVmzPsk4zpgxT/PJJ59ftT0R+S83VwfgCjaLjZ+fXs0nDWex6Y3fCapWnA6f3Evfvx6h4fAmSl5FJM9FRh6ne/deLFr0Bb6+fvzyyzoAJkx4hWHDRrBgwWdUqVKVefNmXfH1p0+f4uOP5zJixCjeemsiTz/9HHPnLuKJJ55iypTJGfWioqKYNWsBb7/9Hu+88wZpaWmZjlOhQkU+/HAW8+YtZtCgR5kx40MAli//hujok8ybt5gFCz7j7rs7YLVaee+9t3n99TeZO3cRnTp1ZebMDzOO5e7uwZw5C7n33h6MHfsMo0aN4ZNPPuf771dw/nwcERFHWLt2DdOnz2X+/MW4uRlZvfp7AFJSUqhTJ4wFC5bQoEFDli//hrCw+tx++x08/vhI5s9frOS1kOq2rCOf7f0UAIvNQrdlHVm67zMAki3JdFvWkWUHvgIgPu083ZZ1ZMWh5QCcTTlLt2Ud+THC+XMUkxxDt2UdWXdsDQAnEiLptqwj64//DEDE+SN0W9aRjSd+A+Bg7AG6LevIpqi/ANhzdjfdlnVka8w/AOw4E063ZR3ZcSYcgK0x/9BtWUf2nN0NwKaov+i2rCMHYw8AsPHEb3Rb1pGI80cAWH/8Z7ot68iJhEgA1h1bQ7dlHYlJjgHgx4jv6basI2dTzgKw4tByui3rSHzaeQCWHfiKbss6kmxJBmDpvs/otqwjFpsFgM/2fkq3ZR0z3suFu+fTY3nXjPLcnbPovaJ7Rnnm9o/ot+qBjPKHW99n4A8PZZTf3zKVoasHZJSnbH6TYWsGZ5Qnb5rAyHXDMsoT/hjPM7+MzCi/8vsLjNkwKqP84m9jePG3MRnlMRtG8crvL2SUn/llJBP+GJ9RHrluGJM3TcgoD1szmCmb38woD109gPe3TM0oD/zhIT7c+n5Gud+qB5i5/aOMcu8V3Zm789Lf0B7Lu7Jw9/yMckH62Uu2JHEw9gAp1hQAktITORh7gFRrKgCJ6QkcjD1AmtX5Nz4hPZ6DsQdIt6VfiN9Zvvizcz7tPAdjD2C1WwGIS43jYOwBbBfKsamxF8o25/60OI4fP0a3+3qyaNEXmLxMfL5yMXfe2YrZsz9h6sf/wz80gBUrljnjtSZzNvUMFSpUxGKxsvNQOIfPH2Lt2tW0bt2Ok/EnmPzO6xntye1t7mDqh5eudXRSFEfjIzLKUUknORZ/NKN8MvEkxxOOZZRPJJ4gMuH4pXJCZMbvHUBkwnFOJJ7IKB9POMbJxJMZ5WPxR4lKulQ+Gh9BdFJURjni/BFikqIzykfOH+bUhd9jgMPnD3E6+VRG+VDcQc6knM4oH4w7wNmUM5fKsQc4d+H33u6wO8up5wCw2W0cjD1AbGrshbKVg7EHiEuNA8B6oXz+wt8Ji83CwdgDxKfFA5BuS+dg7AES0p3lNGsaB2MPkJieAECqNZWDsQdISk8EIMWawsHYAyRbkgDnz76znHyhfO2fvbxQJHtgLQnpnNoWQ51+YdR9pAFB1Yq5OiQRySf2fr6bvUt25ugxaz5Yl5oP1L5mnVKlSlOtWg0AatSoSVTUSRITE0lISKBhw1sA6NChMy+9NOaKr2/Vqi1Go5Hk5GR27AjnpZfGZuyzWNIznrdu3RY3NzfKlStP6dJlOHYsItNxEhMTmTBhPJGRxzAYDFitzn9eNm/+i27demQMT/b3D+Dw4YMcPnyIp59+AgC73Ubx4pc+Lb/99jsAqFKlKpUqVc74JL106TKcOhVDePg29u3bw+DBDwOQlpZKUFAQAB4eHhm9wzVq1OLvv/+65vsnIlJUdVvWkfuq9qBFmTuw2C30XtadntXup1npFqRYk3lgxX08UKMPTUKbEp8Wz+DVD9O31sM0LHkL51LO8uhPj/BwnUeoV6I+p5JjKO1b5rrOGxxSkmrVqgNQtVo1zpw6zeHDh5g1azrn4+NISknCrdnlfWWtW7flt/UbuLtHB9atW8Orr77BieORRB49ntGepFvT8A8KyLk3SQoVg8PhcLg6iBthsdiIi0t2dRgiUohERx8lNLQC4JoENirqJKNHP8XChV8AsHjxQlJSknnggb48/PADfP31SgBOnIjkpZfGMHfup5leP3HieJo3v51WrdqSlJRInz49+PbbHy87z8SJ42nQoBGdOjl7ZJ54YghPPfUcCQnxGZM4TZw4nurVa9KrV2+iok4yYsSjfPnld7zwwnN069aDJk2aZRzv0KGDvPXWRGbMmHfZuYYPH8rw4U9Rs2btyyaJurhv585wzpw5w2OPDb/s9e3atWTNml8B+Pnnn9i48TdeeGF8pu/1Zvz7Wl8UHOx3U8eSS9Q2i+Sd//4d67asI71r9qV3zb5YbBZ6fXcvfWs9TK8avUm2JNNnZU8G1BlEt2o9iE87z8PfP8jgsMfoXKUrZ1POMujHfgxrMIL2FTsQkxxDiHdIljFcrd364YeVTJr0DtWqVWfVqu/YuvUfXnhhfMYQ4j59+mW0ZePHT2L8+BeYO3fRNdsTKfxutG0ukj2wIiJXU/OB2ln2luYVX19f/Pz82b59K/XrN+SHH1bSoEGja77Gx8eXUqXKsG7dT7Ru3RaHw8HBgwcyPiX/+eef6NChM1FRJzl58gTly1dg164dGa9PTEwkODgYgFWrvsvY3qRJU7799msaNmyMu7s78fHnKV++AnFxsezcGU7duvWwWq0cO3aUypWrXNf3d8sttzJu3DM88EAfgoKKER9/nuTkZEJDS131Nd7e3iQnK1ESEbloWbdVGc89jB6Zyt4e3pnK/uaATOXiXsUzla8neb2W5OQkSpQogdVqZfXq7wkOvnwy1DJlyuLmZmTBgtm0adMOINvtiRQtSmBFRPKxF18c/69JnMowbtwrWb7m5Zdf5513JrNgwRxsNitt2tydkcCGhIQyZEh/kpKSePbZcZjN5kyv7dv3YSZMGM+CBXO47bbbM7Z37tyN48ePMWDAgxiN7nTt2o0ePR5gwoQ3ee+9d0hMTMRms3H//Q9e9z8clSpVZsiQYTz99HAcDjtGozujRo25ZgLbps3dvPXWRL788jMmTNAkTiIi+cngwcMYOnQAgYGB1K5d96ofOLZu3Y6PPprG0qXOe4g9PDyy1Z5I0aIhxCJS5F1p6EphlN3ht4WBhhDnDrXNInmnqLRZUnTcaNtcJGchFhERERERkYJHQ4hFRIqIF14Y7+oQRERERLJFPbAiIiIiIiJSICiBFREBCth0AHITdI1FpLDQ3zMpLG7mZ1kJrIgUee7uJpKS4vUPQSHmcDhISorH3d3k6lBERLJFbZYUFjfbNuseWBEp8oKCgomNPU1iYpyrQ5Fc5O5uIigo2NVhiIhki9osKUxupm1WAisiRZ7R6E6JEldfe1RERCS/UJslRZ2GEIuIiIiIiEiBoARWRERERERECgQlsCIiIiIiIlIgGByawkxEREREREQKAPXAioiIiIiISIGgBFZEREREREQKBCWwIiIiIiIiUiAogRUREREREZECQQmsiIiIiIiIFAhKYEVERERERKRAUAIrIiIiIiIiBYISWBERERERESkQ3F0dQH506NAhFixYQFxcHM2aNaNPnz6uDklyyU8//cQvv/xCYmIiPXv25Pbbb3d1SJJLjh8/zvTp00lMTOT99993dTiSw5KTk3n11Vfx8PDg1ltvpWvXrq4OSXKY2uaiQ21z0aG2uXDLrba50PXAjhs3jttuu43OnTtn2r5hwwbat29Pu3btmDlz5jWPUaVKFV577TXee+89tmzZkpvhSjbkxLVu27YtEyZM4NVXX2XVqlW5Ga5kQ05c63LlyjFp0qTcDFNy2I1c99WrV9O+fXsmTJjAunXrXBGuXIPa5qJDbXPRoba5aMoPbXOh64Ht3r07Dz30EGPGjMnYZrPZeO2115g3bx4hISH07NmT1q1bY7PZmDp1aqbXT5o0ieLFi7N27VqWLFnCvffem9ffglynnLrWANOnT6dv3755Gr9cv5y81lJw3Mh1j4mJoUaNGgAYjUZXhSxXoba56FDbXHSobS6a8kPbXOgS2CZNmhAZGZlpW3h4OBUqVKBcuXIAdOrUibVr1/Loo48yY8aMKx6nTZs2tGnThqFDh9KlS5dcj1tuXE5ca4fDwTvvvMMdd9xBnTp18iRuuXE59XstBcuNXPeQkBCio6OpVasWdrvdFeHKNahtLjrUNhcdapuLpvzQNhe6BPZKYmJiCA0NzSiHhIQQHh5+1fp//fUXa9asIT09nTvvvDMvQpQccqPXeuHChfzxxx8kJCRw9OhRHnzwwbwIU3LAjV7r2NhY3n33XXbv3s2MGTN49NFH8yJMyWFXu+79+vXj9ddf55dffqFVq1YujFCul9rmokNtc9Ghtrloyuu2uUgksDeqadOmNG3a1NVhSB54+OGHefjhh10dhuSBoKAgXnvtNVeHIbnE29ubN954w9VhSC5S21x0qG0uOtQ2F2651TYXukmcruRi9/VFMTExhISEuDAiyS261kWHrnXRpOteeOhaFh261kWHrnXRlNfXvUgksGFhYURERHD8+HHS09NZuXIlrVu3dnVYkgt0rYsOXeuiSde98NC1LDp0rYsOXeuiKa+vu8HhcDhy7eguMGrUKDZt2kRsbCzFixdnxIgR9OrVi/Xr1zNp0iRsNhs9evRg2LBhrg5VsknXuujQtS6adN0LD13LokPXuujQtS6a8sN1L3QJrIiIiIiIiBRORWIIsYiIiIiIiBR8SmBFRERERESkQFACKyIiIiIiIgWCElgREREREREpEJTAioiIiIiISIGgBFZEREREREQKBCWwIoXA119/zWuvvQbAkiVLWLZsmWsDEhERKeLUNovkDndXByAiOevBBx/MkeNYrVbc3fUnQkREJLvUNovkHP0GiORjy5YtY86cORgMBmrUqEGHDh2YPn06FouFwMBA3nnnHUqUKJHpNf/73//w9vZm0KBB9OvXjxo1avD3339js9mYNGkS9erVIzk5mddff50DBw5gtVoZPnw4bdu25euvv2b16tUkJydjt9tZtGiRi75zERGR/Elts4hrKYEVyacOHDjA9OnTWbJkCcWKFSMuLg6DwcAXX3yBwWBg6dKlzJ49m7Fjx17zOKmpqXz77bf8/fffPP/886xYsYKPP/6YZs2a8cYbbxAfH0+vXr1o3rw5ALt372b58uUEBgbmwXcpIiJScKhtFnE9JbAi+dSff/7JPffcQ7FixQAIDAxk3759PP3005w+fZr09HTKli2b5XE6deoEQJMmTUhMTCQ+Pp7ffvuNdevWMXfuXADS0tKIiooCoEWLFmogRURErkBts4jrKYEVKUAmTJjAgAEDaNOmDX/99RcffPBBlq8xGAxXLL///vtUrlw5077t27fj5eWVcwGLiIgUcmqbRfKWZiEWyaeaNWvGDz/8QGxsLABxcXEkJCQQEhICcN2zGa5atQqAzZs34+fnh5+fH7fffjuLFi3C4XAAzqFJIiIicm1qm0VcTz2wIvlUtWrVeOyxx+jXrx9ubm7Url2b4cOH8+STTxIQEEDTpk2JjIzM8jhms5lu3bphtVqZNGkSAI8//jiTJk2ia9eu2O12ypYty4wZM3L7WxIRESnQ1DaLuJ7BcfFjHhEpdPr168fo0aMJCwtzdSgiIiKC2maR7NIQYhERERERESkQ1AMrIiIiIiIiBYJ6YEVERERERKRAUAIrIiIiIiIiBYISWBERERERESkQlMCKiIiIiIhIgaAEVkRERERERAoEJbAiIiIiIiJSIPwf+R6oRLxezmMAAAAASUVORK5CYII=",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"with sb.axes_style(\"dark\") as s:\n",
" f, ax = plt.subplots(1, 2, figsize=(16, 6))\n",
"\n",
" ax[0].semilogx(caliper, p_with, \"blue\", label=\"with replacement\")\n",
" ax[0].semilogx(caliper, p_without, \"purple\", label=\"no replacement\")\n",
" ax[0].set_ylabel(\"fraction of data matched \", color=\"blue\")\n",
" ax[0].legend()\n",
" ax[0].set_xlabel(\"caliper\")\n",
"\n",
" ax[1].semilogx(caliper, ATE_with, \"green\",\n",
" label=\"matching (with replacement)\")\n",
" ax[1].semilogx(caliper, ATE_without, \"orange\",\n",
" label=\"matching (no replacement)\")\n",
" ax[1].set_ylabel(\"ATE\", color=\"green\")\n",
" ax[1].hlines(xmin=caliper.min(), xmax=caliper.max(),\n",
" y=ATE_ipw, ls=\"--\", color=\"green\", label=\"ipw\")\n",
" ax[1].hlines(xmin=caliper.min(), xmax=caliper.max(),\n",
" y=ATE_naive, ls=\":\", color=\"green\", label=\"naive\")\n",
" ax[1].legend(loc=4)\n",
" ax[1].set_xlabel(\"caliper\");"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Intermediate results\n",
"The `Matching` object implements `IndividualOutcomeEstimator`, specifically the `fit` and `estimate_individual_outcome` methods. Because matching has many uses, `Matching` also outputs a DataFrame of matches and a vector of weights. These can be used for filtering data before applying regression, for example, or for assessing the quality of the covariate balancing (as shown above)."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
\n",
"
distances
\n",
"
matches
\n",
"
\n",
"
\n",
"
match_to_treatment
\n",
"
sample_id
\n",
"
\n",
"
\n",
"
\n",
" \n",
" \n",
"
\n",
"
1
\n",
"
0
\n",
"
[]
\n",
"
[]
\n",
"
\n",
"
\n",
"
1
\n",
"
[]
\n",
"
[]
\n",
"
\n",
"
\n",
"
2
\n",
"
[0.00020765596476873815]
\n",
"
[819]
\n",
"
\n",
"
\n",
"
3
\n",
"
[0.0001142785835505089]
\n",
"
[1564]
\n",
"
\n",
"
\n",
"
4
\n",
"
[6.015153572747067e-05]
\n",
"
[731]
\n",
"
\n",
"
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
\n",
"
\n",
"
0
\n",
"
1623
\n",
"
[]
\n",
"
[]
\n",
"
\n",
"
\n",
"
1624
\n",
"
[]
\n",
"
[]
\n",
"
\n",
"
\n",
"
1625
\n",
"
[]
\n",
"
[]
\n",
"
\n",
"
\n",
"
1627
\n",
"
[0]
\n",
"
[1627]
\n",
"
\n",
"
\n",
"
1628
\n",
"
[2.757121594376688e-05]
\n",
"
[28]
\n",
"
\n",
" \n",
"
\n",
"
3132 rows × 2 columns
\n",
"
"
],
"text/plain": [
" distances matches\n",
"match_to_treatment sample_id \n",
"1 0 [] []\n",
" 1 [] []\n",
" 2 [0.00020765596476873815] [819]\n",
" 3 [0.0001142785835505089] [1564]\n",
" 4 [6.015153572747067e-05] [731]\n",
"... ... ...\n",
"0 1623 [] []\n",
" 1624 [] []\n",
" 1625 [] []\n",
" 1627 [0] [1627]\n",
" 1628 [2.757121594376688e-05] [28]\n",
"\n",
"[3132 rows x 2 columns]"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"matcher = Matching(with_replacement=False, propensity_transform=propensity_transform)\n",
"matcher.fit(X, a, y)\n",
"match_df = matcher.match(X, a)\n",
"match_df"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The matching DataFrame can be understood in the following way: when sample `sample_id` searched for a match with treatment value `match_to_treatment` it found the samples indexed in `matches` which were located at distances according to the `distances` entry. If `matches` is empty it means `sample_id` had no match in treatment class `match_to_treatment`. This allows us to handle the case with multiple matches, as well as with uneven numbers of matches due to the caliper constraint:"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
\n",
"
distances
\n",
"
matches
\n",
"
\n",
"
\n",
"
match_to_treatment
\n",
"
sample_id
\n",
"
\n",
"
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
0
\n",
"
[0]
\n",
"
[0]
\n",
"
\n",
"
\n",
"
1
\n",
"
[0]
\n",
"
[1]
\n",
"
\n",
"
\n",
"
2
\n",
"
[0]
\n",
"
[2]
\n",
"
\n",
"
\n",
"
3
\n",
"
[0]
\n",
"
[3]
\n",
"
\n",
"
\n",
"
4
\n",
"
[0]
\n",
"
[4]
\n",
"
\n",
"
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
\n",
"
\n",
"
1
\n",
"
1623
\n",
"
[]
\n",
"
[]
\n",
"
\n",
"
\n",
"
1624
\n",
"
[]
\n",
"
[]
\n",
"
\n",
"
\n",
"
1625
\n",
"
[]
\n",
"
[]
\n",
"
\n",
"
\n",
"
1627
\n",
"
[0.0007702527722480701]
\n",
"
[1202]
\n",
"
\n",
"
\n",
"
1628
\n",
"
[0]
\n",
"
[1628]
\n",
"
\n",
" \n",
"
\n",
"
3132 rows × 2 columns
\n",
"
"
],
"text/plain": [
" distances matches\n",
"match_to_treatment sample_id \n",
"0 0 [0] [0]\n",
" 1 [0] [1]\n",
" 2 [0] [2]\n",
" 3 [0] [3]\n",
" 4 [0] [4]\n",
"... ... ...\n",
"1 1623 [] []\n",
" 1624 [] []\n",
" 1625 [] []\n",
" 1627 [0.0007702527722480701] [1202]\n",
" 1628 [0] [1628]\n",
"\n",
"[3132 rows x 2 columns]"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"matcher = Matching(with_replacement=True, \n",
" propensity_transform=propensity_transform, \n",
" n_neighbors=3,\n",
" caliper=0.001).fit(X, a, y)\n",
"match_df_3 = matcher.match(X, a)\n",
"match_df_3"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can also output weights which can be useful for comparing with other methods or preparing the data for further processing:\n"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
control_to_treatment
\n",
"
treatment_to_control
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
0.0
\n",
"
0.0
\n",
"
\n",
"
\n",
"
1
\n",
"
0.0
\n",
"
0.0
\n",
"
\n",
"
\n",
"
2
\n",
"
1.0
\n",
"
1.0
\n",
"
\n",
"
\n",
"
3
\n",
"
1.0
\n",
"
1.0
\n",
"
\n",
"
\n",
"
4
\n",
"
1.0
\n",
"
1.0
\n",
"
\n",
"
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
\n",
"
\n",
"
1623
\n",
"
0.0
\n",
"
0.0
\n",
"
\n",
"
\n",
"
1624
\n",
"
0.0
\n",
"
0.0
\n",
"
\n",
"
\n",
"
1625
\n",
"
0.0
\n",
"
0.0
\n",
"
\n",
"
\n",
"
1627
\n",
"
1.0
\n",
"
1.0
\n",
"
\n",
"
\n",
"
1628
\n",
"
1.0
\n",
"
1.0
\n",
"
\n",
" \n",
"
\n",
"
1566 rows × 2 columns
\n",
"
"
],
"text/plain": [
" control_to_treatment treatment_to_control\n",
"0 0.0 0.0\n",
"1 0.0 0.0\n",
"2 1.0 1.0\n",
"3 1.0 1.0\n",
"4 1.0 1.0\n",
"... ... ...\n",
"1623 0.0 0.0\n",
"1624 0.0 0.0\n",
"1625 0.0 0.0\n",
"1627 1.0 1.0\n",
"1628 1.0 1.0\n",
"\n",
"[1566 rows x 2 columns]"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"weights_df = matcher.matches_to_weights(match_df)\n",
"weights_df"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We distinguish between matching from 0 to 1 and from 1 to 0 because in general they are distinct processes which can have different weights. For the case with no replacement the two columns are identical and always only 1 or 0. We see this when we examine the weights obtained with replacement, caliper and 3 neighbors."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
control_to_treatment
\n",
"
treatment_to_control
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
0.0
\n",
"
0.000000
\n",
"
\n",
"
\n",
"
1
\n",
"
1.0
\n",
"
0.333333
\n",
"
\n",
"
\n",
"
2
\n",
"
1.0
\n",
"
0.500000
\n",
"
\n",
"
\n",
"
3
\n",
"
1.0
\n",
"
1.000000
\n",
"
\n",
"
\n",
"
4
\n",
"
1.0
\n",
"
1.000000
\n",
"
\n",
"
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
\n",
"
\n",
"
1623
\n",
"
0.0
\n",
"
0.000000
\n",
"
\n",
"
\n",
"
1624
\n",
"
0.0
\n",
"
0.000000
\n",
"
\n",
"
\n",
"
1625
\n",
"
0.0
\n",
"
0.000000
\n",
"
\n",
"
\n",
"
1627
\n",
"
1.0
\n",
"
0.500000
\n",
"
\n",
"
\n",
"
1628
\n",
"
4.0
\n",
"
1.000000
\n",
"
\n",
" \n",
"
\n",
"
1566 rows × 2 columns
\n",
"
"
],
"text/plain": [
" control_to_treatment treatment_to_control\n",
"0 0.0 0.000000\n",
"1 1.0 0.333333\n",
"2 1.0 0.500000\n",
"3 1.0 1.000000\n",
"4 1.0 1.000000\n",
"... ... ...\n",
"1623 0.0 0.000000\n",
"1624 0.0 0.000000\n",
"1625 0.0 0.000000\n",
"1627 1.0 0.500000\n",
"1628 4.0 1.000000\n",
"\n",
"[1566 rows x 2 columns]"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"weights_df_3 = matcher.matches_to_weights(match_df_3)\n",
"weights_df_3"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The columns no longer match because a given treatment sample may be selected multiple times when matched from the set of control samples. Each selection increases its weight. However, if it is one of $n$ chosen samples it is only increased by $1/n$. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Compare Covariate Balancing with IPW\n",
"Using the weights we can compare how well the matching algorithm balances the covariate distributions, compared to IPW. Even though IPW is a population effect estimator and matching is an individual outcome estimator, their results can both be expressed as weights and compared using `calculate_covariate_balance`."
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"ipw = IPW(LogisticRegression(solver=\"liblinear\"))\n",
"ipw.fit(X, a)\n",
"ipw_binarized_weights = ipw.compute_weights(X, a)"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"