{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Heterogenous Effect Mixture Model (HEMM) Demo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float:right;\" src=\"https://i0.wp.com/makerfaire.com/wp-content/uploads/gravity_forms/145-c8d5744ce54ae5e6b8c251641327d92d/2017/04/ibmpos_black.jpg?fit=750%2C500&strip=all\" width=\"20%\"/>\n",
    "\n",
    "# Contents\n",
    "\n",
    "### 1. [Introduction](#Introduction) \n",
    "####    &nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;  1.1 [Subgroup Discovery and Heterogenous Treatment Effect Problem](#introsgdisc)\n",
    "####    &nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;  1.2 [HEMM Description and Plate Notation](#introhemm)\n",
    "\n",
    "       \n",
    "### 2. [Synthetic Data Example](#syndata)\n",
    "\n",
    "\n",
    "####   &nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;   2.1 [Data Description and Generative Process](#syndatagen)\n",
    "\n",
    "####   &nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;   2.2 [Estimation of Counterfactual Outcomes, PEHE Estimation](#syndatapehe)\n",
    "\n",
    "####   &nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;   2.3 [Subgroup Discovery](#syndatasg)\n",
    "\n",
    "####   &nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;   2.4 [Bootstrapping PEHE Estimates](#syndatabs)\n",
    "\n",
    "\n",
    "### 3. [IHDP Data Example](#IHDPdata)\n",
    "\n",
    "####   &nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;   3.1 [Data Description](#IHDPdatadesc)\n",
    "\n",
    "#### &nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;   3.1 Estimation of Counterfactual Outcomes, PEHE Estimation\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='Introduction'></a>\n",
    "## 1. Introduction\n",
    "\n",
    "In a large number of decision problems, estimating Heterogenous Outcomes to a treatment is not sufficient from a decision making perspective. A Neural Network might be able to approximate the Outcome and Corresponding Counterfactuals well, but deployment for real world decision making would be limited by lack of transparency and exaplinability of the Neural Model. \n",
    "The Heterogenous Effect Mixture Model approach was originally proposed in the paper\n",
    "[\"***Interpretable Subgroup Discovery in Treatment Effect Estimation with Application to Opioid Prescribing Guidelines***\"](https://arxiv.org/abs/1905.03297) in order to mitigate this challenge and allow decision makers more insight. \n",
    "\n",
    "<a id='introsgdisc'></a>\n",
    "### 1.1 Subgroup Discovery and Heterogenous Treatment Effect Problem\n",
    "\n",
    "\n",
    "The idea behind HEMM is involves assuming a low dimensional clustering or a latent $\\mathcal{Z}$ for each individual in the dataset. An intuitive visual example of the following phenomenon is below in **Figure A**. Notice that, almost all instances receiving treatment in $\\mathcal{Z}_1$ have a positive outcome, while very few in $\\mathcal{Z}_3$ do. We are interested in recovering such latent subgroups. \n",
    "    \n",
    "\n",
    "<table border=\"1\">\n",
    "\n",
    "<tr><td><center><img src=\"https://ndownloader.figshare.com/files/22805153\" title=\"HEMM\" width=\"90%\"/></center><td><center><img src=\"https://ndownloader.figshare.com/files/22805156\" title=\"HEMM\"  width=\"60%\"/></center></tr>\n",
    "<tr><td><center><h3>Fig A. Example of the Heterogenous Effect Subgroup Problem</h3><center> <td><h3>Fig B. The HEMM Model in Plate Notation </h3></tr>\n",
    "</table>\n",
    "\n",
    "<a id='introhemm'></a>\n",
    "### 1.2 HEMM Description and Plate Notation\n",
    "    \n",
    "In the HEMM model, this is carried out using a finite (${K}$) mixture of Normals and Bernoullis. The contribution of performing the Treatment is then mediated by belonging to one of the $K$ subgroups. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='syndata'></a>\n",
    "## 2. Synthetic Data Example\n",
    "This section will demonstrate **HEMM**'s ability on a synthetic dataset."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='syndatagen'></a>\n",
    "\n",
    "### 2.1 Data Description and Generative Process\n",
    "\n",
    "We will use a Synthetic Dataset included with the **`HEMM`** module in **`causallib`** to demonstrate the ability of **HEMM** to discover subgroups, as well as hetereogenous effects. The synthetic dataset is generated according to the following scheme\n",
    "\n",
    "<!---<img src=\"hemm_demo_resources/hemm_truedist.svg\" title=\"HEMM\" width=\"50%\"/> --->\n",
    "\n",
    "<img src=\"https://ndownloader.figshare.com/files/22805159\" width=\"65%\">\n",
    "\n",
    "\n",
    "<center><h3> Fig C. The Generative Process for Synthetic Data </h3></center>\n",
    "\n",
    "\n",
    "1. We take $\n",
    "\\mathbf{X} = (X_0, X_1) \\in \\mathbb{R}^2$ and sample it from a uniform distribution over $\\mathcal{X} = [0,1]^2$. \n",
    "\n",
    "2. In order to simulate the selection bias inherent in observational studies, the treatment variable depends on $\\mathbf{X}$ as $T \\sim \\text{Bernoulli}(0.4)$ for $x_0 < 0.5$ and $\\text{Bernoulli}(0.6)$ for $x_0 > 0.5$.  \n",
    "\n",
    "3. The potential outcomes $Y(0)$ and $Y(1)$ are also Bernoulli with means given by the functions of $\\mathbf{X}$, $T$ shown in the figure.  The figure shows that $p(Y(1)=1 | \\mathbf{X}) > p(Y(0)=1 | \\mathbf{X})$, i.e. treatment increases the the treated population has a higher probability of positive outcome.  Note that under the conditional exchangeability assumption we have $p(Y(t)=1 | \\mathbf{X}) = p(Y=1 | T=t, \\mathbf{X})$.  \n",
    "\n",
    "4. We model the effect of the confounders $\\mathbf{X}$ by assigning higher probability to the upper triangular region of $\\mathcal{X}$. This together with the distribution of $T$ imply that individuals who are more likely to have positive outcome regardless of treatment (upper triangle) are also more likely to receive treatment (right half-square). \n",
    "\n",
    "5. Lastly, we model the enhanced treatment effect group as a circular region $\\mathcal{S} = \\{x: \\lVert x-c \\rVert_{2} < r \\}$, where $p(Y(1)=1 | \\mathcal{S}) > p(Y(1)=1 | \\mathcal{X}\\backslash\\mathcal{S})$. We set $c = ( \\frac{1}{2}, \\frac{1}{2})$ and $r = \\frac{1}{4}$.  A total of $2,500$ samples $(\\mathbf{x}_i, t_i, y_i)$ are generated as described above.\n",
    "\n",
    "\n",
    "\n",
    "The `gen_montecarlo` from `causallib.contrib.hemm.gen_synthetic_data` allows us to randomly sample multiple instantiations of the synthetic dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from causallib.contrib.hemm.gen_synthetic_data import gen_montecarlo\n",
    "\n",
    "syn_data = gen_montecarlo(5000, 2, 100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='syndatapehe'></a>\n",
    "\n",
    "### 2.2 Estimation of Counterfactual Outcomes and PEHE Estimation\n",
    "\n",
    "We will now fit an **HEMM** model to each realization of the Synthetic Dataset and evaluate the capability of **HEMM** to estimate the **ITE** ie. **Individual Treatment Effect**, also known as the **Conditional Average Treatment Effect** (**CATE**) using the **Precision of Estimating Heterogenous Effects Metric (PEHE)**. The **PEHE** is the Mean Squared Error around the ***ITE*** more formally defined as:\n",
    "\n",
    "$$ \\text{PEHE} = \\frac{1}{n} \\sum^{n}_{i=1}  \\left(f_1(\\mathbf{x_i}) - f_0(\\mathbf{x}_i) - \\mathbb{E}[Y(1) - Y(0)] | \\mathbf{X} = \\mathbf{x}_i] \\right)^2.\n",
    "$$\n",
    "Here $f_0$ and $f_1$ represent the functions used to estimate the counterfactual outcome for an observation with features $\\mathbf{x_i}$. The Second term inside the parenthesis is the **True CATE**. Note that the **True CATE** can only be observed in synthetic settings and not real world observational data and hence using **PEHE** as a metric ofestimating goodness of fit for **CATE** is limited to synthetic datasets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from causallib.contrib.hemm import HEMM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Leave 30% of the Training Data as Validation\n",
    "vsize = int(0.3*syn_data['TRAIN']['x'].shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#setting the Training Dataset\n",
    "Xtr  = syn_data['TRAIN']['x'] [:-vsize,:, 0]\n",
    "Ytr  = syn_data['TRAIN']['yf'][:-vsize  , 0]\n",
    "Ttr  = syn_data['TRAIN']['t'] [:-vsize  , 0]\n",
    "\n",
    "#setting the Dev/Val Dataset\n",
    "Xdev = syn_data['TRAIN']['x'] [-vsize:,:, 0]\n",
    "Ydev = syn_data['TRAIN']['yf'][-vsize:  , 0]\n",
    "Tdev = syn_data['TRAIN']['t'] [-vsize:  , 0]\n",
    "\n",
    "#setting the Test Dataset\n",
    "Xte  = syn_data['TEST']['x'] [:,:, 0]\n",
    "Yte  = syn_data['TEST']['yf'][:  , 0]\n",
    "Tte  = syn_data['TEST']['t'] [:  , 0]\n",
    "\n",
    "#Feature size\n",
    "Xdim = Xtr.shape[1]\n",
    "\n",
    "#Number of Components to Discover.\n",
    "K = 3\n",
    "\n",
    "#Initialize the model with Population means and Std. Devs.\n",
    "#Empirically results in faster, and better convergence.\n",
    "mu  = Xtr.mean(axis=0).reshape(1,-1)\n",
    "std = Xtr.std(axis=0).reshape(1,-1)\n",
    "\n",
    "#Set the Learning Rate for Adam Optimizer and Batch Size\n",
    "learning_rate = 1e-4\n",
    "batch_size    = 100\n",
    "\n",
    "#Indicate the Outcome Distribution (Y), could be 'bin' for \n",
    "#Binary Outcomes or 'cont' for Continuous Outcomes.\n",
    "response = 'bin'\n",
    "\n",
    "#Indicate what kind of a model to be used to adjust for \n",
    "#confounding. You can also pass your own PyTorch Model!!!!\n",
    "outcome_model='linear'\n",
    "\n",
    "#Instantiate an HEMM model\n",
    "model = HEMM(Xdim, K, homo=True, mu=mu, std=std, bc=2, lamb=0.0000,\\\n",
    "                spread=.1,outcome_model=outcome_model,sep_heads=True,epochs=10,\\\n",
    "                 learning_rate=learning_rate,weight_decay=0.0001,metric='LL', use_p_correction=False,\\\n",
    "                 response=response,imb_fun=None,batch_size=batch_size )\n",
    "\n",
    "\n",
    "#Fit the HEMM model on the Training Data, with Early Stopping\n",
    "#on the Dev Set.\n",
    "cd = model.fit(Xtr, Ttr,Ytr, validation_data=(Xdev, Tdev, Ydev))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Extract the True Potential outcomes for the Training \n",
    "#and the Test Dataset.\n",
    "mu1tr = syn_data['TRAIN']['mu1'][:-vsize, 0] \n",
    "mu0tr = syn_data['TRAIN']['mu0'][:-vsize, 0]\n",
    "\n",
    "mu1te = syn_data['TEST']['mu1'][:, 0] \n",
    "mu0te = syn_data['TEST']['mu0'][:, 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Estimate the Individual Potential Outcomes using HEMM\n",
    "inSampleCFoutcomes  = model.estimate_individual_outcome(Xtr, Ttr)\n",
    "outSampleCFoutcomes = model.estimate_individual_outcome(Xte, Tte)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "In Sample PEHE: 0.011010778786300154\n"
     ]
    }
   ],
   "source": [
    "#Compute the True In-sample Individual Treatment Effect (ITE)\n",
    "trueCATE = mu1tr - mu0tr\n",
    "\n",
    "#Compute the In-Sample  ITEs\n",
    "inSampleITE = inSampleCFoutcomes[1]-inSampleCFoutcomes[0]\n",
    "\n",
    "#Compute the In-Sample PEHE\n",
    "inSamplePEHE = ((inSampleITE - trueCATE)**2).mean()\n",
    "print(\"In Sample PEHE:\", inSamplePEHE)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "In Sample PEHE: 0.010915930183215048\n"
     ]
    }
   ],
   "source": [
    "#Compute the True Out of Sample Individual Treatment Effect (ITE)\n",
    "trueCATE = mu1te - mu0te\n",
    "\n",
    "#Compute the Out of Sample  ITEs\n",
    "outSampleITE = outSampleCFoutcomes[1]-outSampleCFoutcomes[0]\n",
    "\n",
    "#Compute the Out of Sample PEHE\n",
    "outSamplePEHE = ((outSampleITE - trueCATE)**2).mean()\n",
    "\n",
    "print(\"In Sample PEHE:\", outSamplePEHE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, **HEMM** returns a **low PEHE** score, justifying its ability to estimate Heterogenous Effects to Treatment."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='syndatasg'></a>\n",
    "### 2.3  Subgroup Discovery"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now demonstrate **HEMM**'s ability to extract the Subgroup with the **Enhanced Treatment Effect**. We will first use the `get_groups_effect` member function of the **`HEMM`** class to extract the absolute values of the Treatment Effect Terms ( **$\\gamma_k$** in Figure B )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.69842033 0.71704322 0.50003658]\n",
      "Subgroup 2 demonstrates a Higher Effect to Treatment\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "gamma = model.get_groups_effect()\n",
    "print(gamma)\n",
    "\n",
    "K = np.argmax(gamma) \n",
    "\n",
    "print (\"Subgroup\", K+1, \"demonstrates a Higher Effect to Treatment\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have identified the Subgroup with Higher Effect to Treatment, we will compute the ability of **HEMM** to extract the subgroup by comparing with the True Subgroup Assignment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "ZinSample  = model.get_groups_proba(Xtr)\n",
    "ZoutSample = model.get_groups_proba(Xte)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have the group assignment scores, we will compute the ROC Curve and Precision-Recall Statistic in Recovering the Subgroups"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "In Sample ROC AUC: 0.9738558368695354\n",
      "Out Sample ROC AUC: 0.9772955843564394\n"
     ]
    }
   ],
   "source": [
    "from sklearn.metrics import roc_auc_score, roc_curve\n",
    "\n",
    "trueInSampleZ  = syn_data['TRAIN']['z'][:-vsize, 0] \n",
    "trueOutSampleZ = syn_data['TEST']['z'][:,0] \n",
    "\n",
    "inSample_roc_auc  = roc_auc_score( trueInSampleZ , ZinSample[K] )\n",
    "outSample_roc_auc = roc_auc_score( trueOutSampleZ, ZoutSample[K])\n",
    "\n",
    "print (\"In Sample ROC AUC:\", inSample_roc_auc)\n",
    "print (\"Out Sample ROC AUC:\", outSample_roc_auc)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us now Plot the ROC Curves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'roc_curve' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-1-308359c939f0>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0minfpr\u001b[0m \u001b[0;34m,\u001b[0m \u001b[0mintpr\u001b[0m \u001b[0;34m,\u001b[0m \u001b[0m_\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mroc_curve\u001b[0m\u001b[0;34m(\u001b[0m \u001b[0mtrueInSampleZ\u001b[0m  \u001b[0;34m,\u001b[0m \u001b[0mZinSample\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mK\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0moutfpr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mouttpr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0m_\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mroc_curve\u001b[0m\u001b[0;34m(\u001b[0m \u001b[0mtrueOutSampleZ\u001b[0m \u001b[0;34m,\u001b[0m \u001b[0mZoutSample\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mK\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mmatplotlib\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mpyplot\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtitle\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"ROC Curve\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'roc_curve' is not defined"
     ]
    }
   ],
   "source": [
    "infpr , intpr , _ = roc_curve( trueInSampleZ  , ZinSample[K] )\n",
    "outfpr, outtpr, _ = roc_curve( trueOutSampleZ , ZoutSample[K] )\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "plt.title(\"ROC Curve\")\n",
    "plt.plot(infpr , intpr, label = \"In-Sample\")\n",
    "plt.plot(outfpr, outtpr, label= \"Out-Sample\")\n",
    "plt.plot(np.linspace(0,1,100),np.linspace(0,1,100),c='k', label='Random', ls=':')\n",
    "plt.xscale('linear')\n",
    "plt.legend()\n",
    "plt.xlabel('FPR')\n",
    "plt.ylabel('TPR')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "In Sample  PR AUC: 0.9486824676876224\n",
      "Out Sample PR AUC: 0.9545712084631546\n"
     ]
    }
   ],
   "source": [
    "from sklearn.metrics import precision_recall_curve, average_precision_score\n",
    "\n",
    "\n",
    "inSample_pr_auc  = average_precision_score( trueInSampleZ , ZinSample[K] )\n",
    "outSample_pr_auc = average_precision_score( trueOutSampleZ, ZoutSample[K])\n",
    "\n",
    "print (\"In Sample  PR AUC:\", inSample_pr_auc)\n",
    "print (\"Out Sample PR AUC:\", outSample_pr_auc)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3gVVfrA8e+bTgkdFAglIEqvAQsWEAVEBbvIuiuuig11dy2LZRVR1PWHZV3RFZW1C4gNVxBRQCwgCQKhQ6gJIIROaCHJ+/tjJuESbpKbkLk35f08z30yc+bMnXcuet8758ycI6qKMcYYk19YqAMwxhhTNlmCMMYY45clCGOMMX5ZgjDGGOOXJQhjjDF+WYIwxhjjlyUIUyGIyDIR6VVEnaYikiEi4UEKy3MiskFELnKXR4rIB6GOyVQcliCMp9wvsEPuF/M2EfmviFQv7eOoajtVnV1EnU2qWl1Vs0v7+O6X81H3PPeIyC8icnZpH+dkiEgNEXlZRDa5caa46/VCHZspmyxBmGC4XFWrA12B7sBj+SuIo7z/9zjRPc96wCzgkxDHk0dEooDvgXZAf6AGcA6wE+hRgveLKNUATZlU3v+HNOWIqm4GpgHtAURktoiMFpGfgYNACxGpKSJvi8hWEdksIk/7NgmJyG0iskJE9ovIchHp6pb7NrX0EJEkEdnnXrW86JY3FxHN/XITkUYiMkVEdrm/pm/zOc5IEZkkIu+5x1omIgkBnmcW8CHQWETq+7znZSKyyOcKo6PPtiYi8pmIpIvIThF51S1vKSIz3bIdIvKhiNQqwcf/J6ApcKWqLlfVHFXdrqpPqepU91gqIqf5xPSOiDztLvcSkTQR+buI/A781/13uMynfoQbY+6/yVnuee4RkcVFNQGasscShAkaEWkCDAAW+hT/ERgGxAIbgXeBLOA0oAvQF7jV3f9aYCTOl10NYCDOL+D8/gX8S1VrAC2BSQWE9DGQBjQCrgGeEZE+PtsHAhOAWsAU4NUAzzPKjXEnsNst6wqMB24H6gJvAFNEJNpNgP9zz7850Ng9LoAAz7oxtgGauJ9BcV0EfKOqGSXYN9epQB2gGc6/2cfADT7b+wE7VPU3EWkMfA087e7zAPCpb8I0ZZ8lCBMMX4jIHuAn4AfgGZ9t76jqMvdXdx3gEuAvqnpAVbcDLwGD3bq3As+raqI6UlR1o5/jHQVOE5F6qpqhqvPyV3CT1bnA31X1sKouAt7CSVi5flLVqW6fxftApyLO8zr3PA8BtwHXuOeFu/6Gqv6qqtmq+i5wBDgLp4mnEfCge96HVfUnAPccZ6jqEVVNB14ELigiDn/qAltLsJ+vHOAJN5ZDwEfAQBGp6m4f4pYB3AhMdT+/HFWdASTh/EAw5YQlCBMMV6hqLVVtpqp3uV8uuVJ9lpsBkcBWt1liD84v7Qbu9ibA2gCOdwtwOrBSRBJ9m0F8NAJ2qep+n7KNOL/ec/3us3wQiHGbUf7gdvJmiMg0nzqTVLUWcAqwFOiW79zuzz0v99yauHE0ATb6JJM8ItJARCa4zW37gA9w+jiKayfQsAT7+UpX1cO5K6qaAqwALneTxECOJYhmwLX5zvfcUojBBJF1NJlQ8x1OOBXnV3U9f1+W7vaWRb6h6hrgBrfT+ypgsojUzVdtC1BHRGJ9kkRTYHMA7/8hTh9DQdt3iMjtQKKIfKSqW93YR6vq6Pz13budmopIhJ/zfhbnM+qoqjtF5AoCbOrK5zvgaRGppqoHCqhzEKjqs34qThNc3qn52Se3mSkMWO4mDXDO931Vvc3PPqacsCsIU2a4X6TfAi+4t2SGuZ20uU0qbwEPiEg3966n00SkWf73EZEbRaS+quYAe9zi425tVdVU4BfgWRGJcTuMb6GQL/5instKYDrwkFv0JnCHiJzpxl5NRC4VkVhgPk7zz3NueYyI9HT3iwUygD1uu/6DJQzpfZwv7U9FpLX72dYVkUdEJLfZZxEwRETCRaQ/gTVlTcDpJ7qTY1cP4FzpXC4i/dz3i3E7uuNKGL8JAUsQpqz5ExAFLMfp4J2M2yyhqp8Ao3G+iPYDX+D0W+TXH1gmIhk4HdaDfZtGfNyA0ym8Bfgcp319Rimey/8Bw0Skgaom4fRDvOqeVwowFMDt47gcp2N+E86v9uvd93gS5/bgvTidvp+VJBBVPYLTUb0SmAHsw0lM9YBf3Wr3uXHsAf6A8/kW9b5bgbk4t8xO9ClPBQYBjwDpOMnpQew7p1wRmzDIGGOMP5bNjTHG+GUJwhhjjF+WIIwxxvhlCcIYY4xfFeY5iHr16mnz5s1DHYYxxpQrCxYs2KGqfodAqTAJonnz5iQlJYU6DGOMKVdExN9wNYA1MRljjCmAJQhjjDF+WYIwxhjjlyUIY4wxflmCMMYY45dnCUJExovIdhFZWsB2EZFXxJnqMTl3mkJ3200issZ93eRVjMYYYwrm5RXEOzijahbkEqCV+xoGvA4gInWAJ4AzcWbaekJEansYpzHGGD88ew5CVeeISPNCqgwC3lNnONl5IlJLRBoCvYAZqroLQERm4CSaj72I82DGXpInPOnFWxtTbDWrRNKmYY2S7Rx/PsSfV7oBmUotlA/KNeb46SbT3LKCyk8gIsNwrj5o2rRpiYI4fDCDHqnjS7SvMV7QNSeWSdF7wdrv4baZHkRkKqtQJgh//81rIeUnFqqOA8YBJCQklGhiizoNGsOTe4quaIzHPklK5cHJySeUX9TmFC5u24BmdatxVov8M6e6PrwOMrZ5HKGpbEKZINJwJmvPFYczs1caTjOTb/nsoEVlTIhcm9CEC1s3IEchTEBEuOyVH/luxTa+W+F8+Y+5thPZOTkczVaysnOoHhPJ1V0bO7+q9qbClHsgJxtystxX9vF/Ndt/+Qnbs6BhJ7j+g5B+Jia0QpkgpgDDRWQCTof0XlXdKiLTgWd8Oqb7Ag+HKkhjgqlu9ejj1r/92wXsPpDJ2z+t551fNvDAJ4tP2Gfl1n080jiBsC2/wZoZEBYBEub8zXvlW5dwiIgpuM7WZEix5qrKzrMEISIf41wJ1BORNJw7kyIBVPU/wFRgAM7cvAeBm91tu0TkKSDRfatRuR3WxlQ21aMjqB4dweOXtWXoOc0JEyEiXIgIE5Zu2cuf30nirZ/Wk9ruQlp3voLM7Bwys5xXtiq3nBtPy/rVi3/g6Y9C0n9L/4RMuVJh5qROSEhQG83VVDa/bdrNVa/9krceGS5EhYcRFibsP5wFwEVtGnDreS0K7r/wJzdBPLqltEM2ZYyILFDVBH/bKsxw38ZURl2b1mblU87jRrmJAUBVufXdJDbtOsh3K7bz3YrtDO7ehOeu7hjKcE05Y0NtGFPOxUSGExMZnpccwOngfntod2b87QLu7NWSxrWqMGWxXQ2Y4rEEYUwF9/f+rbm0Y0MOZmYzZvoqKkqzsvGeJQhjKoFru8UB8OqsFP749nxLEiYgliCMqQRanRLLx7edBcBPKTu4b8Ii1u84EOKoTFlnCcKYSuLslnX5/K5zAJiyeAu9x8xm+77DxX+jrEw4sAN2roWM7aUcpSlL7C4mYyqRLk1r88ODvbh3wiIWp+5h+McLmXT72f4rZx2CD66Gw3vh8D73716nPFdULIzYCGHhwTkBE1SWIIypZJrVrcbnd55Di0emMn/9Lvq/PIcdGUdoXLsqb9zYjVNrxkCTM2HNt3BwF8TUgBqNIKYmRNeAmFpO2YYfYcVXzjAdliAqJEsQxlRCYWHCC9d2YmJSKrWqRLLy9/3syMjkrGe/56azm/HwgEuJaTuw8Dc5st9JEKbCsgRhTCV1dbc4rnbvbtq69xAfzNvI2FlreXfuRv6XvJWep9Vj3+Gj3HNhK7o1szm7KiPrpDbG0LBmFR7s15qv7z0XgP1HskhO28PsVelc/fovbCjsjqdZo+Gr++D7p8Bun61QbCwmY0yB/vxOIjNXbqdmlUgWP9H3+I2LPoYv7nBGhg2PcjqvH94M0SUYHNCETGFjMVmCMMYUKDMrh+6jv+NgZhYDOzXm8NFsHrusDfWrRxMRHgZHMiCyKsx9FWb8Ay78BxzZ59wGG10D+j8LUvR8eCZ0bLA+Y0yJREWEcc+Fp/F/01fx6W9pAHy9ZCsAy57sR7Xcq4WqdZy/M59yribCoyAzA85/EKoVYxRZU6bYFYQxJiBZ2Tm8OGM1783dSMaRLNo0rMFUt88CVWRf2rFbYee/CdMehIRboFp96DXCriTKKGtiMsaUmv2Hj9Jh5LeAM6HRgcwsqkVF8Nld53D6KbFOpZVfw4QhOFPMKzywBqo3OP6NsjLdebQVajUN5ikYH5YgjDGlavW2/Tz99QpiYyL4OnlrXrnT7OS2XB89DIs+gK/vhx63O01O+7fC/m2Q8Tsc3OnuJXD/Sog9NfgnYgpNEJ7e5ioi/UVklYikiMgIP9ubicj3IpIsIrNFJM5nW7aILHJfU7yM0xhTPKefEst7f+7B2CFdSR7Zl05NagHQ7onprNi6z6kUGQOxjZzlpLdh3Ww4tAdqN4O2g6DXI9DlRkAh9VfYsjAk52IK5tkVhIiEA6uBi4E0nDmmb1DV5T51PgH+p6rvisiFwM2q+kd3W4aqBny/nF1BGBM6B45k8Ye3fmVR6h4ANjx36bGNh/dBVHUI8/N7dPmXMOlPx9bvS3YSiAmaUF1B9ABSVHWdqmYCE4BB+eq0Bb53l2f52W6MKQeqRUfkjRQLcNGLP/DAJ4vJyVFn3CZ/yQHgtIth0Gtw1t3OeqYNQV6WeJkgGgOpPutpbpmvxcDV7vKVQKyI5N4TFyMiSSIyT0Su8HcAERnm1klKT08vzdiNMcUkInxwy5m0aViDlO0ZTF6QxmT31tgCRVWFLn+AJj2CE6QpFi8ThL972vK3Zz0AXCAiC4ELgM1AlrutqXvZMwR4WURanvBmquNUNUFVE+rXr1+KoRtjSuLcVvWYdt95fHl3TwD2HTpavDdIS4SU723IjjLCywSRBjTxWY8Djps1XVW3qOpVqtoFeNQt25u7zf27DpgNdPEwVmNMKWpRvxrgzF53MDOriNpARLTz96t74YOrYGfKiXUO74Nty2DVN7BksjPMuPGUl09SJwKtRCQe58pgMM7VQB4RqQfsUtUc4GFgvFteGzioqkfcOj2B5z2M1RhTiiLDnd+es1el0/bx6SQ0q81bNyVQq2qU/x1aXgjXveckgB/+CYlvg+bA3lTYkwp7NzmTFfmq0RiaFTDZkSkVnl1BqGoWMByYDqwAJqnqMhEZJSK5A833AlaJyGrgFGC0W94GSBKRxTid18/53v1kjCnbYiLDmTK8JzVinN+gSRt389y0lQXvEBHt3Poaf4Gz/uvrsPhj2L0RajaGDtfBRU/CNePhkv9z6mRnenwWxh6UM8Z4KjMrh9MfmwbA9L+czxmnxhZcWRX2bHKG7KhSy3+dDT/DOwPgT1OgxQUeRFy5hOxBOWOMiYoIo23DGgD0e3kOl/37R+asLuCuQxHnOYiCkoOvH/4Jb/aBV7vDrvWlGLHJZQnCGOO5KcN7cs+FpwGwdPM+np++krd+XMf89buK/2Y146BqXdi32Wlm2rEadq4t5YgNWBOTMSaIsrJzGPLmr8zfcCwxvHNzd3qd0aCQvQqRmghvXwR/+BRaXVRKUVYu1sRkjCkTIsLD+NcNnfn0znMYeXlbAJZu3lvEXiZULEEYY4KqYc0qdGtWmz+c5Yy5tHHnwRBHZApiCcIYExK5Qy18siCNf323JqSxGP8sQRhjQiIiPIwn3Gaml75bzc6MIyGOyORnCcIYEzI394znprOdpqZ/frOSLXsOleyNksbDR4Ph393guabw+5JSjLLysgRhjAmpv1/SGoBJSWmc89xMBvzrR1b9vj+wnavVAwmDlO9gz0aIbegMybF7o4cRVx6WIIwxIVU1KoLxQxM497R6ACzfuo9+L8/h13U72X+4iNFg68TDiE3w6Fa4ay70eyYIEVceXg7WZ4wxAbmw9Slc2PoUADqOnM6+w1lcP24eAK1PjeWNP3ajWd1q/neO9jN0x08vOsNw+NtmAmZXEMaYMmXxE3156or2NKoZA8DK3/dz3RtzA9u5qjvf2OYFsO4HyEiHrck2NHgJ2ZPUxpgyKzMrh/Yjp5OZlUO/dqdwSo0Y7uzVkoY1qxS808a58N/+EBEDWYedsiGfwOl9gxN0OWNPUhtjyqWoiDCeGtQOgOnLtvHe3I30f/nHwndq1BnaXQWdBkPP+5yyyTfDvq0eR1vxWB+EMaZMu757UwZ0aMi+w1lc+sqP7Dl4lEWpe+jcpIARXyOrwLX/dZYP7Yaf/wWZGbDhJ+h4bfACrwDsCsIYU+bFxkTSuFYVXryuEwBXjP2ZC8fMLnrHKrVh+AJvg6vALEEYY8qNnqfV44rOjQBYt+MApz0ylR2BPoG9dibMHA0HdnoYYcXiaYIQkf4iskpEUkRkhJ/tzUTkexFJFpHZIhLns+0mEVnjvm7yMk5jTPkQHRHOy4O7MPN+Zya5rBzlirE/F75TVFXn7+KPYM7zsGa6x1FWHJ4lCBEJB8YClwBtgRtEpG2+amOA91S1IzAKeNbdtw7wBHAm0AN4QkRqexWrMaZ8aVG/OmtGXwJA2u5D9H95DgXekVmjEdz+ozNFKcAXdzq3wJoieXkF0QNIUdV1qpoJTAAG5avTFvjeXZ7ls70fMENVd6nqbmAG0N/DWI0x5UxkeBjPX9MRcJ6ViH94Ko99sYRDmX6eeWjYEZqfBx2uc9bn/B+8ezm83AE2zQti1OWLlwmiMZDqs57mlvlaDFztLl8JxIpI3QD3RUSGiUiSiCSlpxcwx60xpsK6LqEJPz7Um9gY54bMD+Zt4rOFaf4rh4XBgOchPBrSkpy+iD2bYPuKIEZcvniZIMRPWf5rwAeAC0RkIXABsBnICnBfVHWcqiaoakL9+vVPNl5jTDnUpE5Vlozsxxd39wTg0c+X0nzE17w3dwM5Ofm+NqrUhofWwSOb4cZPgx9sOeNlgkgDmvisxwFbfCuo6hZVvUpVuwCPumV7A9nXGGN8dW5Si4f6n5G3/viXy/h6iZ+H46KrQ1h4ECMrv7xMEIlAKxGJF5EoYDAwxbeCiNQTkdwYHgbGu8vTgb4iUtvtnO7rlhljTIHu6nUa658dwH9v7g7AI58v4eeUHXy1eEvBndi71sFv78HX98N7g2DLoiBGXLZ59iS1qmaJyHCcL/ZwYLyqLhORUUCSqk4BegHPiogCc4C73X13ichTOEkGYJSq7vIqVmNMxSEinBXvDNq3/3AWf3jrVwAaxEZTLzaacBGa16t27Cril1ecv5HV4OgB+GYEDJkEMTVCEX6ZYoP1GWMqpN0HMpmyeAvp+4/w6qyUvPI61aL47R8XOyuLJ0BENJzaEao3gGfdR7EufQG63xqCqIOvsMH6bCwmY0yFVLtaFDed0xxVZevew1SNCmdtega/rN3JWz+u49bzWjgD+vka+jW8cylkZYYm6DLGEoQxpkITEV5wx3D6ZulWflm7k6e/XkFkeBg3ndP8+MqntHf+LvrQSR5V6wQ32DLGxmIyxlQa/ds3ZMy1TrJ4YsoyFqfuOb5CeJTzd9tSmP8mrJkBv7xaaYcKtwRhjKlUrukWxx0XtARg6Za9x2+MqgrD3GE4Zj8DH14D3z4KP75QKWelswRhjKl0/tyzOeA8VPfhrxuPvwW2UWe46Em4eBTc+JlTlvgm/PZu8AMNMeuDMMZUOrWrReUtP/r5UtZuP8AjA1oTEe7+Zj73L8cqD/w3TLkHDla+O+3tCsIYU+lEhoeR+OhFvH9LDwDG/7yev01a7L9yx8H+yysBSxDGmEqpfmw057Wqzwtup/WUxVv4zw9r2X3AbnHNZQnCGFOpXd0tjvFDnefEnpu2ki5PzWDBxsrXnOSPJQhjTKV3YetTuK9PK/q0bgDA8q37T6y0/ItK1w9hCcIYY4C/Xnw6D/VvDcA/vlhKdu5Q4bnjif6+BH56EXJyQhRh8FmCMMYY1+mnVKdJnSoAfLrAnXgoPOLYsxG//BsWjC9g74rHEoQxxrhEhHdvdu5s+vDXjcc2NOoMl77oLB/YGYLIQsMShDHG+GhRvzoAi9P2snDT7mMbut0coohCxxKEMcbkc05LZz6JK1/7hQ07DuTbWjGmSAiEJQhjjMln7JCudG9eG4BeY2bz4rerjm2c/Sws+hiyjoQouuCxBGGMMfnUrhbFhGFn0/uM+gC8MjOFOSk74fwHnQrTHoKnT4HEt0MYpfc8TRAi0l9EVolIioiM8LO9qYjMEpGFIpIsIgPc8uYickhEFrmv/3gZpzHG5BceJvz35h7c2csZ+fVP4+fzee2hcPolEJcAKOzbEtIYveZZghCRcGAscAnQFrhBRNrmq/YYMElVuwCDgdd8tq1V1c7u6w6v4jTGmML8vX9rLml/KgB/nbgYhkyAP34OEh7iyLzn5RVEDyBFVdepaiYwARiUr44CuTOD1wQqdjo2xpRLr9/YjQEdTqWOzyiwlYGXCaIxkOqznuaW+RoJ3CgiacBU4B6fbfFu09MPInKevwOIyDARSRKRpPT09FIM3Rhjjle7ahRhEuoogivgBCEijUXkHBE5P/dV1C5+yvLfH3YD8I6qxgEDgPdFJAzYCjR1m57+BnwkIjXy7YuqjlPVBFVNqF+/fqCnYowxxabAjoxMpi/7PdShBE1ACUJE/gn8jNNn8KD7eqCI3dKAJj7rcZzYhHQLMAlAVecCMUA9VT2iqjvd8gXAWuD0QGI1xhgvnBlfB4CZK7aHOJLgCfQK4grgDFUdoKqXu6+BReyTCLQSkXgRicLphJ6Sr84moA+AiLTBSRDpIlLf7eRGRFoArYB1AcZqjDGlblDnxlSNCmdiUioZR7JCHU5QBJog1gGRxXljVc0ChgPTgRU4dystE5FRIpKbXO4HbhORxcDHwFB1Joc9H0h2yycDd6hq5Rpn1xhT5oS7nRC/rqsc4zEFOif1QWCRiHwP5D0+qKr3FraTqk7F6Xz2LXvcZ3k50NPPfp8CnwYYmzHGBMX7t5zJFWN/RivJaBuBJogpnNg8ZIwxlUq4VK7bmAJKEKr6rtuPkNtRvEpVj3oXljHGmFAL9C6mXsAanCejXwNWB3CbqzHGVEiTklKLrlQBBNpJ/QLQV1UvUNXzgX7AS96FZYwxZc8Zp8YCEBURBpoNP46B35eGOCrvBJogIlU1b7xbVV1NMe9qMsaY8i4qIox61aP5X/JWsqo3dAqXVtz7aQJNEEki8raI9HJfbwILvAzMGGPKos5NagHQbsezqN8BIyqOQBPEncAy4F7gPmA5YCOsGmMqnTHXdiQiTDhCFEc1jPnr0mHbMsg8GOrQSp1oBbmhNyEhQZOSkkIdhjGmkvjstzQu+7IjUZLtFFz42LEJhcoREVmgqgn+thV6BSEik9y/S9wJfY57eRGsMcaUB1d1jSM1/lo+yTqfoxrOkYw9oQ6p1BX1HMR97t/LvA7EGGPKm5ZD3+Dtz5dw7eJzYf6r0Kg9dL4h1GGVmkKvIFR1q7u4A0hV1Y1ANNAJm9zHGGP4x6Vt+SCrj7OyrWLd8hpoJ/UcIEZEGgPfAzcD73gVlDHGlBdVosJ5LOsWDmh0qEMpdYEmCFHVg8BVwL9V9UqceaaNMabSO/90Z8KyRal7QxxJ6Qo4QYjI2cAfgK/dskAH+jPGmApt5OVtEWDLxlUs3VhxJhQKNEH8BXgY+Nyd06EFMMu7sIwxpvxoUb86VeUIA8Ln8824R5i1smIkiYAShKr+oKoDVfWf7vq6ouaCMMaYyiRr4OsA1JYMPvx1Y4ijKR1FPQfxsvv3KxGZkv8VnBCNMabsi+g6BKJiCQsTvluxneS08v9cRFH9CO+7f8eU5M1FpD/wLyAceEtVn8u3vSnwLlDLrTPCnYUOEXkYuAXIBu5V1eklicEYY4Kp9amxsAlSdx2iY1ytUIdzUgpNEKqaOyBfEnBIVXMARCQc53mIArl1xgIXA2lAoohMcacZzfUYzlzVr4tIW5zpSZu7y4OBdkAj4DsROV1Vs4t9hsYYE0StT60Bm+DX9Tu5tGPDUIdzUgLtpP4eqOqzXgX4roh9egApbn9FJjABGJSvjgI13OWaHHv4bhAwQVWPqOp6IMV9P2OMKbsyM6iROpO67GX9jgOhjuakBZogYlQ1I3fFXa5aSH2AxoDvtEtpbpmvkcCNIpKGc/VwTzH2RUSGiUiSiCSlp6cHch7GGOOdiBjCd69lZsxDrE9ZzoZyniQCTRAHRKRr7oqIdAMOFbGPv4HS8w8dewPwjqrGAQOA90UkLMB9UdVxqpqgqgn169cvIhxjjPHYn74EoCb7GRo+ncQNu0Ic0MkJ9GG3vwCfiEhuE1BD4Poi9kkDmvisx3Hi+E23AP0BVHWuiMQA9QLc1xhjypamZ8LwBfBqN8LJCXU0Jy3Q5yASgdY4EwfdBbTx6cAuSCLQSkTiRSQKp9M5/62xm4A+ACLSBogB0t16g0UkWkTigVbA/MBOyRhjQqjeaeREO3cvPTg5mf2Hj4Y4oJILKEGISFXg78B9qroE506jQocAV9UsYDgwHViBc7fSMhEZJSID3Wr3A7eJyGLgY2CoOpYBk3BmrvsGuNvuYDLGlBci5F1B9B7zQ4ijKblAm5j+izMH9dnuehrwCfC/wnZyn2mYmq/scZ/l5UDPAvYdDYwOMD5jjCkz5PAe/hQxgzStx7iMyzmYmUXVqPI3fF2gndQtVfV54CiAqh7Cf0eyMcaY3o8CcF59516eRZv2cCiz/DWCBJogMkWkCu6dRCLSEjjiWVTGGFOeXfAQVKlD83rVABjy1q+0efwbjmSVryQRaIJ4AqcvoImIfIjz4NxDnkVljDHl3aFdNEn5kJda/kbuXfojpywLbUzFVGSCEBEBVuJMFjQUpzM5QVVnexqZMcaUZ52GAHDl5jGs6zyBWA7y8fxUdmSUn8aXIhOEqirwharuVNWvVfV/qrojCLEZY0z5deXrcPEoAMJWfsWoWs5ca+vSy8/T1YE2Mc0Tke6eRmKMMRVNz/tgyCQA+sRu4umIt4nd+kuIgwpcoPdd9QbuEJENwAGcO5hUVbVQz/UAABkDSURBVDt6FZgxxlQIp/eDmk2okb6AGyPg57mRkNAbIqpAWKC/0UMj0ARxiadRGGNMRXb1Wxw4sI/9E4bRYd8ceLYJnH039H0q1JEVqtAE4Y6NdAdwGrAEeNt9QtoYY0ygmp5FNWBllTOoeXATNcJ3w7rZcGgPVCm7kwoVdX3zLpCAkxwuAV7wPCJjjKmgvmr7IleFvQxRVeH3ZPh+VKhDKlRRCaKtqt6oqm8A1wDnBSEmY4yp2K74j/P38F7nr54wm0GZUFSCyBuG0JqWjDHm5B06mg2n94Xa8bBrLXw2DP7ZDL57MtShnaCoTupOIrLPXRagirueexdTjYJ3NcYY4+vAkSyOZitv/LCW28OjYMtC2LUesjJh17pQh3eCQq8gVDVcVWu4r1hVjfBZtuRgjDHFcFXXOACenbaSteeNgT9+AQ+mQK2mIY7Mv7J9E64xxlQgZ7esy929WwKwOqwVtOwN4ZEhjqpgliCMMSaILu/UKNQhBMzTBCEi/UVklYikiMgIP9tfEpFF7mu1iOzx2Zbtsy3/VKXGGFOurdq2P9QhFMmzKY5EJBwYC1yMMwNdoohMcWeRA0BV/+pT/x6gi89bHFLVzl7FZ4wxoVCnahQAv6zdyV8uCnEwRfDyCqIHkKKq61Q1E5gADCqk/g04Q4kbY0yF1aBGDHWqRTF//S627DkU6nAK5WWCaAyk+qynuWUnEJFmQDww06c4RkSSRGSeiFzhXZjGGBNcsTFO483FL/4Q4kgK52WC8DdndUGPCw4GJquq73x8TVU1ARgCvOxOc3r8AUSGuUkkKT09/eQjNsaYIPjqnnMBOJCZzda9ZfcqwssEkQY08VmPA7YUUHcw+ZqXVHWL+3cdMJvj+ydy64xT1QRVTahfv35pxGyMMZ6rERPJg/3OAODsZ2eScaRsDlThZYJIBFqJSLyIROEkgRPuRhKRM4DawFyfstoiEu0u1wN6Asvz72uMMeXVDT2aEua2s2zde5gtO3bCss9hx5rQBubDswThjt00HJgOrAAmqeoyERklIgN9qt4ATHCnNs3VBkgSkcXALOA537ufjDGmvKtTLYq1zwzg+gSnoaXR9h/hk6Ewa3RoA/MhWkZHESyuhIQETUpKCnUYxhhTLKrKa8//nar7N3BjzWQim3SDwR8G7fgissDt7z2BPUltjDEhJCKcfvn9PJl1E9syjsLK/0HiW6EOC7AEYYwxIXdOy7oA/JLdzilYPT2E0RxjCcIYY0KsWnQE9/ZpxUNZt7Ml5jQIKxsD+FmCMMaYMuCuXs6jXocys9ECHxkLLksQxhhTBsREhgNwJCuHGcu3sTPjSIgjsgRhjDFlxuQ7zs5bTtywO4SROCxBGGNMGZHQvA4t6lcLdRh5LEEYY0wZEpaZQd/wBdROTwx1KJYgjDGmLDlavSEAVea/EuJILEEYY0yZEjn4A9K1BjsOZBdd2WOWIIwxpgyJqtmAbVqH2uzlzSdvQcecDq+fC+mrgx6LJQhjjClj4utXp0tYCn/O+RTJ2AbblsD0h4Meh2dzUhtjjCmZan3/weLkhdy6oAm1JYNvo/8OWcF/LsKuIIwxpqw5oz8drh7BtMeuJfzUdizIaUWmBv/r2hKEMcaUQWFhQr3q0VzUpgEA2TnBH37DEoQxxpRh8fVC9+CcJQhjjDF+eZogRKS/iKwSkRQRGeFn+0sissh9rRaRPT7bbhKRNe7rJi/jNMYYcyLP7mISkXBgLHAxkAYkisgU37mlVfWvPvXvAbq4y3WAJ4AEQIEF7r6hH73KGGMqCS+vIHoAKaq6TlUzgQnAoELq3wB87C73A2ao6i43KcwA+nsYqzHGmHy8TBCNgVSf9TS37AQi0gyIB2YWd19jjKkMIvanwcFdQT2mlwlC/JQVdJ/WYGCyquYOPhLQviIyTESSRCQpPT29hGEaY0zZVoVMIvesg6/uC+pxvUwQaUATn/U4YEsBdQdzrHkp4H1VdZyqJqhqQv369U8yXGOMKZuez7rOWcjYHtTjepkgEoFWIhIvIlE4SWBK/koicgZQG5jrUzwd6CsitUWkNtDXLTPGmEql52n1mJ3ThS01OkFEdFCP7dldTKqaJSLDcb7Yw4HxqrpMREYBSaqamyxuACaoqvrsu0tEnsJJMgCjVDW4jW/GGFMG1K0WFbJjezpYn6pOBabmK3s83/rIAvYdD4z3LDhjjDGFsiepjTGmHNi851DQj2kJwhhjyrCI8GNf04eOBneWOUsQxhhTxj1/dUcAgj2gqyUIY4wp42pUcbqLDxzJCupxLUEYY0wZF1+vOgCrt+0P6nEtQRhjTBl3xqmxecvj5qwN2nEtQRhjTDnQtmENANJ2B+9uJksQxhhTDlSLiiAiPIww8TdUnTcsQRhjTHlw9ABnsYT2u74N2iE9fZI61I4ePUpaWhqHDx8OdSiVSkxMDHFxcURGRoY6FGMqjlrN4PcltFz3AfDXIquXhgqdINLS0oiNjaV58+ZIEC/LKjNVZefOnaSlpREfHx/qcIypOK4Zz9JRPTikwRubqUI3MR0+fJi6detacggiEaFu3bp21WZMaYuIplpsTcD5IRYMFTpBAJYcQsA+c2O8ke0+Sv3bpj1BOV6FTxDGGFNR1KvuzAfx+uzgPAthCcJj1atXL1b9nJwc7r33Xtq3b0+HDh3o3r0769ev9yg6R/PmzdmxY4enxzDGnLya7pAbEWHBuUqv0J3U5dHEiRPZsmULycnJhIWFkZaWRrVq1UIdljGmDBCcxPDNst/5fe9hTq0Z4+nxKk2CePKrZSzfsq9U37Ntoxo8cXm7gOrOnj2bkSNHUq9ePZYuXUq3bt344IMPTmiv37p1Kw0bNiQszLm4i4uLy9t25513kpiYyKFDh7jmmmt48sknAecKYMiQIcyaNYujR48ybtw4Hn74YVJSUnjwwQe54447mD17No8//jh169Zl1apVnH/++bz22mt5x8n1wQcf8Morr5CZmcmZZ57Ja6+9Rnh4+Ml8TMaYUlS7ahTshb9/msw7N3f3tM/P0yYmEekvIqtEJEVERhRQ5zoRWS4iy0TkI5/ybBFZ5L5OmMu6PFq4cCEvv/wyy5cvZ926dfz8888n1Lnuuuv46quv6Ny5M/fffz8LFy7M2zZ69GiSkpJITk7mhx9+IDk5OW9bkyZNmDt3Lueddx5Dhw5l8uTJzJs3j8cfPzaB3/z583nhhRdYsmQJa9eu5bPPPjvu2CtWrGDixIn8/PPPLFq0iPDwcD788EMPPgljTIkcPUCbI4v4a8Rkflidzrx13s7E7NkVhIiEA2OBi4E0IFFEpqjqcp86rYCHgZ6qultEGvi8xSFV7Vxa8QT6S99LPXr0yLsi6Ny5Mxs2bODcc889rk5cXByrVq1i5syZzJw5kz59+vDJJ5/Qp08fJk2axLhx48jKymLr1q0sX76cjh2dceIHDhwIQIcOHcjIyCA2NpbY2FhiYmLYs2dP3vFbtGgBwA033MBPP/3ENddck3fs77//ngULFtC9e3cADh06RIMGvv8kxpiQan4ebF3M9bVW8NIOeGbqCr6659yi9yshL5uYegApqroOQEQmAIOA5T51bgPGqupuAFXd7mE8IRcdHZ23HB4eTlZWFr/++iu33347AKNGjWLgwIFER0dzySWXcMkll3DKKafwxRdf0KJFC8aMGUNiYiK1a9dm6NChxz1rkPveYWFhxx0nLCyMrCxnDPn8l6L511WVm266iWeffbZ0T9wYUzr6jYYti6ibnQM7YMnmvRzMzKJqlDdf5V42MTUGUn3W09wyX6cDp4vIzyIyT0T6+2yLEZEkt/wKfwcQkWFunaT09PTSjT5IzjzzTBYtWsSiRYsYOHAgv/32G1u2bAGcO5qSk5Np1qwZ+/bto1q1atSsWZNt27Yxbdq0Yh9r/vz5rF+/npycHCZOnHjC1UufPn2YPHky27c7eXrXrl1s3Ljx5E/SGFOqIsPDuPfC0zw/jpdXEP56TvI//hcBtAJ6AXHAjyLSXlX3AE1VdYuItABmisgSVT3u5l9VHQeMA0hISAjyZHze2L59O7fddhtHjhwBnGah4cOHExMTQ5cuXWjXrh0tWrSgZ8+exX7vs88+mxEjRrBkyRLOP/98rrzyyuO2t23blqeffpq+ffuSk5NDZGQkY8eOpVmzZqVybsaY0lMt2vt7jLw8QhrQxGc9Dtjip848VT0KrBeRVTgJI1FVtwCo6joRmQ10AYI3U0YpycjIAKBXr1706tUrr/zVV1/1W79///7079/f77Z33nnHb/mGDRvylocOHcrQoUP9bqtatSoTJ04sdP/rr7+e66+/3u9xjDGVi5dNTIlAKxGJF5EoYDCQ/26kL4DeACJSD6fJaZ2I1BaRaJ/ynhzfd2GMMcZjnl1BqGqWiAwHpgPhwHhVXSYio4AkVZ3ibusrIsuBbOBBVd0pIucAb4hIDk4Se8737idTfPmvYIwxpiieNmKp6lRgar6yx32WFfib+/Kt8wvQwcvYjDHGFM7GYjLGGOOXJQhjjClPso/Axp+ov2+p54eyBGGMMeVJZFUAeq553vNDWYIIgrS0NAYNGkSrVq1o2bIl9913H5mZmYXu88wzzxS4bdu2bVx22WV06tSJtm3bMmDAgNIO+TgbNmygffv2nh7DGBOgy16CsAhyxOlC/tvExXkTCZU2SxAeU1WuuuoqrrjiCtasWcPq1avJyMjg0UcfLXS/whLE448/zsUXX8zixYtZvnw5zz33XGmHbYwpq+q2hKZnExnhjLL8zbLf2X2w8B+cJVVphvtm2gj4fUnpvuepHeCSwr+cZ86cSUxMDDfffDPgjMH00ksvER8fT3x8PMuXL897aO6yyy7jgQce4JtvvuHQoUN07tyZdu3anTCi6tatW+nbt2/eeu6AfRkZGQwaNIjdu3dz9OhRnn76aQYNGsSGDRvo378/5557LvPmzaNTp07cfPPNPPHEE2zfvp0PP/yQHj16MHLkSNauXcvmzZtJTU3loYce4rbbbjvu2NnZ2YwYMYLZs2dz5MgR7r777ryxpIwxwVO3WhSPDmjD6Kkr2LbvcN5sc6XJriA8tmzZMrp163ZcWY0aNWjatGneIHr5Pffcc1SpUoVFixb5HW777rvv5pZbbqF3796MHj06b+ymmJgYPv/8c3777TdmzZrF/fffnze5eUpKCvfddx/JycmsXLmSjz76iJ9++okxY8Ycd7WSnJzM119/zdy5cxk1alTee+d6++23qVmzJomJiSQmJvLmm296PuOdMeZEgtCoOkSSxa3vJnlyjMpzBVHEL32vqKrfCT0KKg9Ev379WLduHd988w3Tpk2jS5cuLF26lFq1avHII48wZ84cwsLC2Lx5M9u2bQMgPj6eDh2cR0vatWtHnz59EBE6dOhw3FAbgwYNokqVKlSpUoXevXszf/58Onc+Nur6t99+S3JyMpMnTwZg7969rFmzhvj4+BKdizGmhLYtZcC0c6kS047nq/zjpL5TClJ5EkSItGvXjk8//fS4sn379pGamkrNmjXJycnJK/cdvtvX2LFjefPNNwGYOnUqjRo1ok6dOgwZMoQhQ4Zw2WWXMWfOHPbv3096ejoLFiwgMjKS5s2b571n/iHAfYcH972SCWRI8H//+9/069evuB+FMaa0VD8Fso8i4VH0joPet5znycxy1sTksT59+nDw4EHee+89wGnDv//++xk6dCgtWrRg0aJF5OTkkJqayvz58/P2i4yM5OjRo4DTpJQ7JHijRo2YOXMmBw8eBGD//v2sXbuWpk2bsnfvXho0aEBkZCSzZs0q0VDdX375JYcPH2bnzp3Mnj07b/KgXP369eP111/Pi2316tUcOHCgRJ+NMaaErngdHloHjbsi+39HNi/w5DB2BeExEeHzzz/nrrvu4qmnniInJ4cBAwbwzDPPEBUVldf00759e7p27Zq337Bhw+jYsSNdu3Y9oR9iwYIFDB8+nIiICHJycrj11lvp3r078fHxXH755SQkJNC5c2dat25d7Hh79OjBpZdeyqZNm/jHP/5Bo0aNjmuCuvXWW9mwYQNdu3ZFValfvz5ffPFFiT8fY0wJREQBUZB1GPZugi+Hw11zoZSvIiS3E7O8S0hI0KSk4ztqVqxYQZs2bUIUUfkzcuRIqlevzgMPPHDS72WfvTFBsGoabF4AnYdAnRYlegsRWaCqCf622RWEMcaUV2dc4rw8YgnC5Bk5cmSoQzDGlCEVvpO6ojShlSf2mRtTMVToBBETE8POnTvtCyuIVJWdO3cSExMT6lCMMSepQjcxxcXFkZaWRnp6eqhDqVRiYmKIi4sLdRjGmJPkaYIQkf7Av3CmHH1LVU94nFlErgNGAgosVtUhbvlNwGNutadV9d3iHj8yMtKe8DXGmBLyLEGISDgwFrgYSAMSRWSK79zSItIKeBjoqaq7RaSBW14HeAJIwEkcC9x9d3sVrzHGmON52QfRA0hR1XWqmglMAAblq3MbMDb3i19Vt7vl/YAZqrrL3TYD6O9hrMYYY/LxMkE0BlJ91tPcMl+nA6eLyM8iMs9tkgp0X0RkmIgkiUiS9TMYY0zp8rIPwt8z3/lvJ4oAWgG9gDjgRxFpH+C+qOo4YByAiKSLSPEHHzqmHrDjJPYvjyrbOVe28wU758riZM65WUEbvEwQaUATn/U4YIufOvNU9SiwXkRW4SSMNJyk4bvv7MIOpqr1TyZYEUkq6HHziqqynXNlO1+wc64svDpnL5uYEoFWIhIvIlHAYGBKvjpfAL0BRKQeTpPTOmA60FdEaotIbaCvW2aMMSZIPLuCUNUsERmO88UeDoxX1WUiMgpIUtUpHEsEy4Fs4EFV3QkgIk/hJBmAUaq6y6tYjTHGnKjCjOZ6skRkmNunUWlUtnOubOcLds6VhVfnbAnCGGOMXxV6LCZjjDElZwnCGGOMX5UqQYhIfxFZJSIpIjLCz/ZoEZnobv9VRJoHP8rSFcA5/01ElotIsoh8LyIF3hNdXhR1zj71rhERFZFyf0tkIOcsIte5/9bLROSjYMdY2gL4b7upiMwSkYXuf98DQhFnaRGR8SKyXUSWFrBdROQV9/NIFpGu/uoVi6pWihfOnVRrgRZAFLAYaJuvzl3Af9zlwcDEUMcdhHPuDVR1l++sDOfs1osF5gDzgIRQxx2Ef+dWwEKgtrveINRxB+GcxwF3usttgQ2hjvskz/l8oCuwtIDtA4BpOA8anwX8erLHrExXEIGMDTUIyB01djLQR6SUZwEPriLPWVVnqepBd3UezkOJ5Vkg/84ATwHPA4eDGZxHTmbcs/IqkHNWoIa7XJMTH9QtV1R1DlDY7f6DgPfUMQ+oJSINT+aYlSlBBDK+U14dVc0C9gJ1gxKdNwIa08rHLTi/QMqzIs9ZRLoATVT1f8EMzEMnM+5ZeRXIOY8EbhSRNGAqcE9wQguZ4v7/XqQKPWFQPoGM7xTQGFDlSMDnIyI34gyvfoGnEXmv0HMWkTDgJWBosAIKghKPe6aqezyOzSuBnPMNwDuq+oKInA28755zjvfhhUSpf39VpiuIQMeGagIgIhE4l6Xl+QnuQM4ZEbkIeBQYqKpHghSbV4o651igPTBbRDbgtNVOKecd1YH+t/2lqh5V1fVA7rhn5VUg53wLMAlAVecCMTiD2lVUAf3/XhyVKUEEMjbUFOAmd/kaYKa6vT/lVJHn7Da3vIGTHMp7uzQUcc6quldV66lqc1VtjtPvMlBVk0ITbqk4mXHPyqtAznkT0AdARNrgJIiKPC/AFOBP7t1MZwF7VXXrybxhpWli0sDGhnob5zI0BefKYXDoIj55AZ7z/wHVgU/c/vhNqjowZEGfpADPuUIJ8JwLHPesPArwnO8H3hSRv+I0tQwtzz/4RORjnCbCem6/yhNAJICq/genn2UAkAIcBG4+6WOW48/LGGOMhypTE5MxxphisARhjDHGL0sQxhhj/LIEYYwxxi9LEMYYY/yyBGFMMYhItogsEpGlIvKViNQq5fcfKiKvussjReSB0nx/Y4rDEoQxxXNIVTuranucZ2XuDnVAxnjFEoQxJTcXn8HQRORBEUl0x+J/0qf8T27ZYhF53y273J1zZKGIfCcip4QgfmMKVWmepDamNIlIOM4wDm+7631xxjbqgTNo2hQROR/YiTPOVU9V3SEiddy3+Ak4S1VVRG4FHsJ58teYMsMShDHFU0VEFgHNgQXADLe8r/ta6K5Xx0kYnYDJqroDQFVzB3+MAya64/VHAeuDEr0xxWBNTMYUzyFV7Qw0w/liz+2DEOBZt3+is6qepqpvu+X+xrP5N/CqqnYAbscZSM6YMsUShDEloKp7gXuBB0QkEmfQuD+LSHUAEWksIg2A74HrRKSuW57bxFQT2Owu34QxZZA1MRlTQqq6UEQWA4NV9X13SOm57qi4GcCN7gijo4EfRCQbpwlqKM5sZ5+IyGacIcfjQ3EOxhTGRnM1xhjjlzUxGWOM8csShDHGGL8sQRhjjPHLEoQxxhi/LEEYY4zxyxKEMcYYvyxBGGOM8ev/Adu4u1ygpUJ1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "inpr , inre , _ = precision_recall_curve( trueInSampleZ  , ZinSample[K] )\n",
    "outpr, outre, _ = precision_recall_curve( trueOutSampleZ , ZoutSample[K] )\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "plt.title(\"Precision-Recall Curve\")\n",
    "plt.plot(inre , inpr, label = \"In-Sample\")\n",
    "plt.plot(outre, outpr, label= \"Out-Sample\")\n",
    "plt.legend()\n",
    "plt.xlabel('Recall')\n",
    "plt.ylabel('Precision')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.4 Qualitative Assesment of the Discovered Subgroups\n",
    "In the previous section, we attempted to perform some Quantitative assesment of the capability of **HEMM** to discover the subgroups with **Enhanced Treatment Effects**. Since we had access to the **True Subgroup Assignment Label**, we proceeded to compute the **ROC** and **Precision-Recall** characteristic. In this section, we will attempt to perform a more qualitative assessement by plotting the Subgroups with a heatmap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Let's first create a uniform grid over the feature space. \n",
    "\n",
    "d = 100\n",
    "\n",
    "x1 = np.linspace(0, 1, d)\n",
    "x2 = np.linspace(0, 1, d)\n",
    "\n",
    "Xgrid = np.meshgrid(x1, x2)\n",
    "\n",
    "Xgrid[0] = Xgrid[0].reshape(-1,1)\n",
    "Xgrid[1] = Xgrid[1].reshape(-1,1)\n",
    "\n",
    "Xgrid = np.hstack([Xgrid[0], Xgrid[1], np.ones_like(Xgrid[1])])\n",
    "\n",
    "Zgrid = model.get_groups_proba(Xgrid)[K].values\n",
    "\n",
    "Zgrid = Zgrid.reshape(d, d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlYAAAIMCAYAAAA3sOKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOydeZwcRfn/38/s7Oa+L0I4whGOcN9y35dcIqgcKghfDkVFDkFUMPpTPAC5BBUUAQXkhiD3YcJNDgh3kACBBEgICbmP3dmp3x/Vs5mpqpnqncxmN9nnnVe/Nl1dVV3d091TU/XpzyPGGBRFURRFUZQVJ9PeDVAURVEURVld0I6VoiiKoihKjdCOlaIoiqIoSo3QjpWiKIqiKEqN0I6VoiiKoihKjdCOlaIoiqIoSo3QjlWViMgYEVklvSpEZKqITG3vdnQGRGSUiBgR2au921JARIYnbbqxvduiKCsTEbk0ufa3b++2KKsvnb5jldxkrVlO7ABtHh5o1xIR+UxEXhSRP4nI7u3dTqU6RGRjEbleRKYkn+siEflARB4TkYtEZEh7t1FpH8rc+7Flr/ZudxrEcryIPCEin4tIU/L3TRG5SUSObe82Kkoasu3dgA7ALwNpPwL6AFcCc51tk9q8RemZB1yR/D8L9Ae2Ar4LnCEijwHfNsbMdMrtu/KaqLQGEdkHeBDoCrwAPAIsBoYDWwP7A88D7meqdA7mEn5m/SL5G9o2tc1aU1v+BRwHLMTeAx8CvYH1ga8COwG3tVvrFCUlnb5jZYwZ5aYlo1J9gCuMMVNXcpNaw9wy7V8f+DtwAPCIiOxsjFla2G6MeW/lNVFpJX/FdqpONMbc5G4UkS2BL1Z6q5QOgTFmLjDKTReRXyTbvW2rAiJyELZT9R6wq/tjUES6ADoKr6wSdPqpwBVFRLIi8lMReVdElonINBH5vYg0lMm/iYjcmORbJiIzReRWEdm4Vm0yxrwPHAJMxo5ynO60wdNYiUiDiPxQRF4WkS9EZHGS734R2a/McdyQ5FmWTEM+IyLfDeTdV0QeEZE5IrJURP4nIr8TkT5Ovski0igiA0PHJSI/SaY2znDS10qmP99P2jJbREaLyA6BOlo0TyJynIi8JCILA+djJxG5S0RmJG2aJiJ/FZE1y7Rtu+QYF4jI/GQ6Y+dQ3nKIyGBgQ2BeqFMFYIx5zRgzzSlnRGRMmTpvTLYPL7N9ExG5L/lsFonIsyJyQJm8fUTkChGZnnyOk0XkbBFZXwKaraJ9ry8iPxCR18RObY4pypMRkdNFZHzyOSxK/v9dEck49VXUhklA95h8zib53HdOPpd5yef0qKTU2iRljYjcUyHP28n11z9ZFxE5QUSeF5FZyTmbluz3G2n2W0vETqu9ISL9ReRqEflIRHIicm6y/a7kGL37T0QOTbadG9g2SKx26Z3kGL9IjnGvVjRvl+TvvwMj7BhjlhljnnD2+/2kTUcH2tQz2fafcjsUkVOLrskZyf1d7tmzm4g8lVyjc5N7fVsJaLaK9y0ia4udxvxURJqL25psuy75HBrFfhfcIfbHk7v/stowEdk82fYnJ73weQ4RkQtk+XfUR2K/o3qUOzfKitHpR6xqwK3YX1IPA/OBLwPnAYOB7xRnFPur7B6gHngAmAKshR3mPkRE9jbGvFyLRhljFovIpcDfgONZPmVYjhuBY4E3gJuBJcCawG7AQUDLQ01EDgHuBLpgp6puA/pipyHPA/5clPe0ZH1RUuYzYC/gfOAwEdk1+RUOcBNwcdKOqwNt/DbQCPy7qP5tgcew06CPYs/vQOArwLMicqQx5qFAXedgp9UeAP6LHaEs1Pkd4HpgGTAamAaMAP4vafOXjDEfFeXfJTk/Dcn+p2A7tGOApwL7Lsc8IAf0FJGhxphPW1G2GtbDTje+gR0pGwp8A3hYRI4zxtxeyCgiXbHHsi3wCnAL9pz9jPhIwpVJngeBh4Dmom3/xI5UTMNeqwY4ErgWe+0dv0JHuJydgAuwn9M12A7sV4E9ROQAY8wzlQobY14QkXeAQ0VkgDFmdvF2EdkR2AS42xgzJ0n+TbLPD4A7sJ/vUGAH4GvA7ax8ugNPY6/Vh7DTzNMqlqiAiGyEvS6GYe+jB7HTd4cDT4rIt4wxt6aoqnA+N6q2La3kQqwk4nZsm/cGTgX2TO7vFgmIiByIfU6AfYZ9iL0PngUqXTdrAC9hn3l3YgcyZid1boL9HAZhn1v/wt6PR2OvscOMMU/W5EjhOmBX7DW4EPuj+zxgFxHZxxjTVKP9KAWMMbo4C1aTYIDhFfKMSfJMBPoXpffAfrE2A2sUpffDTuF8Dox06toMe8G/nLJ9w5N9T43k2yDJlwOyzvFNLVrvA+SBCUBdoJ4BRf8fiP2CaAT2DORdq+j/62I7J/OBTZx81yZtu64obVhy3iYE6t0hyX93UVo2OddL3bZgO4UfA58CXYrSRyX1LAK2Cexno+TYpgDDnG37JO27tyhNsCODBjjCyX9mkm6AvVJ+tncl+d8DzsV2CLpHyhhgTJltN7rXctH1Y4BLnPzbA03Jtdq7KP3CJP9tgBSlrw3MSrbdWGbfHwPrBdp2bLL9ZaCncw9NSLYdF2j3jWWOdQxgnLS9io71+862I5L0d4FMis/mglA9ybZrkm2HFaXNBqaHPj9gYJrrIe1SOMZIns+TfKOBrhWuPa9twKHJtnOd9PHY58vhTvqA5L6YD/RN0f4NsPekSdrxjSRNKpT5fpL/6MC2nsm2/zjpl7L8/nefw9cl264sSqtPPsM8/jPm3KJra/vAvg3wl9C1he2UGeBMJ32/ZF+fUPrcutTdT9G2zZNtfyrzeX4CDC1Kz2IHAgxwVi2vQ12Sc9zeDeiIC63rWO0X2PbLZNuhRWmFL9kzytR3ebJ9ZIr2DSddx6pr0Q0+2Dm+qUXrvZM8z1V6kCV5z3EfPhXy/izJe3FgW7/kobvEeYA8lpTZzMn/pyT98KK0whfjJWX2XzjnXy5KG5WkXR75HA4ps/1e7BdJr2R91yT/2EDeOmwHrTUdq37A3cnDtfDZNQOvAr8GhgTKVNuxmls4jjJlTihKK/xY8O6Jos/5xjL1nFmmbY8n2w8IbNs32fZUoN03lqlvDOU7VsHOE8vv4z1TfDZrJedgvJPegO1EzaT0B8xs7GhVl1jdK7oUrpVInkLHaoMy21vVsSq69v9Rpr7jk+3fTnkMB2NHg0zRMhc7ovR1nGcTK9ax8p5f2NGjxck+M0naQUn+0YH89cBHlO9YLQT6BMptnGyfXOaavDfZ/tVAu6vpWHmdp6Iyr7f1tdkZF50KXHEmBNIKQ+v9itIKeputRGRUoExhCHxT4K3aNA0p+r8pl8kYM19EHgAOAyaJyN3YIe6XjDGLnexfSv4+nGL/2yZ/vekwY8wXIvIKsAd2CuXVZNON2Cm6E7DD1YjVqx2DHRkpntYrnNN1y5zTEcnfTZ1yAOPKtLlQ554S0Ghhp3jrsJ/XRJYf41g3ozGmWUSexf7yToUx5gvgKLGaqAOxI0g7AFsmy3dF5CBjzPi0dVbgZWPMgkD6GOz53wa4SUR6Y49hmgm/zPFsZD/lzvW22A7kmMC2sdhOzDaRutPyjDEmH0gfA+yZ7Mf7DIsxxkwXkSeB/UVkpDGmcJ8ehp2KvtwYkysqcgvwA+BNEbkzqf8FY8y8FTuUFWKOqd3LK4V7ZVCZ+29Y8nfTNJUZYx4W++LNHtip422wnbcvJ8vxInKUc46rJXS/zhKRt4DtsG8iTmH59edd48aYJhF5CTtqG+KdMp914Zkxpsw1+RRWyrANVlqwooSO9Q0RmQ1sJiINxpjGGuxHSdCO1Qpiiubiiyjc+HVFaQOSv6dEquy5wo1aTkFo3Uz8TbJvYHVPx7H8le2lInIX9ldqQVDaN/n7cYr9F3RL5bRChfS+RWn3YkeyvikiFxhjmrG/lgdg39IsfqgWzunXIu0IndMZZfIW6vxxyjoLx1jO/qDcfiqSdGD+miyIyFrY6dPDsPqvraup1yHW5sKx9Y7kj1k/lDsHfbBf9N5D3RiTE5HPsR3ZWpD2WGPcyPKO//lJ2gnJX/eFg7OwU7onAT9JlpyIPAScY4yZknKftaSq67EMhXvlkGQpR+pnWnK//zdZEBFJ6r4Rq9s6EavFW1HSXg+x+7vStV/puofWPRdXhErHOgDoxXKNm1ID9K3AlUfhl8tWxhipsATfBquSvZO/E2O/8owxS4wxo4wxGwHrAN/E/kr7JnZIuUChIzmMOIVjXqPM9qFOPowxS7Aiy6HYLzAo/8VVKHdE5JyGvH3KjeAV6uwTqXOsk7+caWe5Y28Vxpjp2FG7Ruyo54DizZT/kVTp4Rxrc+HY5kfyxwxLK53r/iJS724QkSxWzze/KLnw674tjzVGcce/TkQGYaewXjXGvFqc0RjTbIy50hizVbL/o5Lyh2NtULqk3GctKTtyTeXzGzq3hXN2cuRe+UHVjbX8B/h/SdI+K9DeYlbGtR97xqR+LtJ2x2qA0Ki1sgJox2rl8WLyd6V4sYhId6weCuyURGqMMdOMMbdgp6LeBXYr+iIvHMfBKap6Jfm7V6B9fbGjLkuBt53NNyZ/T0hefz4YeM0Y45qztsU5bW2dhbc493Q3iEgd9s22WrEM27Fy+YLAdESy/0ojW9uKSK9A+l7J31fAThUD7wPDJGzbUO0xvoJ9Bu0R2LYHdsS3+C3Zwqhr6Fh7U/mNst3EsW9I2KuoLVGKOv5rYoXGx2O/7Cr+IDLGfGaMuccY83XsVM8GWJ1LR6Ls+cVOSbuszGda4cu/WN7Q2vYWE7pfBwEjsR2a95PkwnXhXePJD4KdIvsJUahzz2REzqXwgzjVtU91x7o5drTqTZ0GbAPaW+TVERdaIV4vs+3EpPyJRWkDsDfHZ8COgTIZ0guch1NBvI59bfe/LH/jqouzfSql4vVBwE6Benphh6WbWC7WLn4rcI9AmbWcdjZiR7k2dPJdnbTv+jLH8D+skLTwNlpIgFmP1UEspkig7uTZmaK3slguXg+ea6zeqzHZ/0aB7Q3A7kXrNXsrEPs23IUEBOrJ9sKLA2866YU3fA5w0n9RtP/hRenDi9LLvRU4l9K3Agt1VfNWYPA+wk47G6wGq/gz6p6kGeB4p8zb2Kn2kUVpdSy3ajBO/r2KjnWF3gosKlcQbd+Cvb+aKHo5JMnTBSvAdwXX9dgvVgNsWpQ+NLn2PLFzyjZ5xx7I8znwRoXtJyf1XOekb499ycRQKl4XrM6wCTi2TJ3bAP1StP8I7DR3NrCtb3KeDfDdovTCW8+vUfoCzGCsTrWSeL3SW4FXOZ/Xx9gRoz2c/LG3Av9T4XifT/Kc7qTvnezrU4re3GT5yxyPFF+rWC3YTCqL1yu9FXh2NdebLpUX1VitJIwxsxNzuHuBFxMR7JvYm2gdbAdgAPZNvrT0LRKNZrFi+a2SujLYm/AEY8yySD3Dkja9jX2ATcPqag7FDhdfZRKRszHmcxE5DnvT/ldEHsY+2HpjxdVrYzt2GGOmisiPsK+ivywid2C/hPdM2jiZ5ToVl5uxw/8XYr9IPS8cY8WjX8X6wDwoIs9jQw4tTtqxA/bBMzRJi2KMmSwiJwE3YEXHj2A7WfXYz2n35Bg2SfIbETkZ+4bb3WINJKdgP4f9sJ/BQWn2nezjV8AvRGRccixfYIXRuwJbYL8QTnfKXYodXbxfRG4H5mANF9fD/gDYq8z+ngb+T0R2wr4RWvCxygCnGTtSVeAPWEHtMcDGYsMl9cG+rfV0si0kxC2LMeZWETkiqeNNEbkP+7D/StL2O4wdOS3mEmxUgecSQfhS7JdRPfYFiK3K7O4R4DIROTjJV/CxWoqdykrddmPMcyIyBavtqwceMMZ85mTrhvXMmpoInD/E3tv7Y8Xco40xxSO1v8VOeX+H5SO2K5vbgYuAU0RkA+yzYD3s1OW92M+pheTa/xrwJHCriJyDtV+Yj73/tsHeJ1sQ13hugb3fZ4vIM9h7KJ/Ucyj2R95YivRVxpj3RORerO/ZK8m92i/J/xSVRfNPAOOS++Uz7DW0E/Zev6hoH00icgpwP/B4cs19lBzbnti3mA+gldc+1hPvaeDPInI49l4fjvWxasQ+t5cW5f8v9kWpA7HP6qex9+sRLH9rshzjgdeSthd8rEZi7/mQX6CyorR3z64jLrTBiFXRtuFY64B3sQ/1+dgOxj+Br6Rs33BKX0k2SV2fYYfnrwZ2ixzf1KL1vtiHyVPYX2fLsL+YxmC9hjwLBqz31s1J/kbsr6axwKmBvAdgH0BfJHVPwX5Rl/W3wXZimpNjeyByPgYDv8MaXS7GPjzexXb+vknpK/CjSDGChH3Q34j9QlyG7awUjDT3CeTfDvvlvSBZnsB2HlPtL6kjg+2E/RFrLPgJdjRgAbbzekW5axL75TchuQ5mY01U16Wy3cKN2C+f+5PPZjH2YXtgmX30Ba5K2rUsuW7PAXZM6rvCye/tu8wxfy9p++JkmQicQZlRJOzIyptJG2Ykn8kAKtstjEo+jyew99yC5JrcocpnxM9Zfu8dFdhej32r9WHsF/FSbIf8RWzHuKHMufKeGSnbs8IjVkme9bFvos1NPosXsF/EQR+rouviF9jOwaKk3HtYv6yTCHhmBeoYgn2x5w7sqOTc5NqfmXxmpxAezeqONaAtXJPvAGez3EKm3IjV9sBpwOvJZzMzuY4GlWnfHtjOzSLsiP0j2Lf7Cp/bhkV5oyNWSb51sR3F6dhn6CzsM2vrMvkHJfv7PGnzJKxpcsxuYQjwU+wzcRn2h/MfgB7VXGu6xBdJPoAOiYhcgL14t8P+cvrQGDO8inq+jX1DZxPsQ/UB4AJjzKzatVZROifJL/rrsNMaf23v9hQjNqzKf4FfmlU0jp7ScRGRSVhfql6mNjYQNSN5o/sobGfx8/ZuT2eio4vXL8a+BfIeVQaeFZGzsMLSeVjNy1+x0xljNFaSoqRHAnESRWRtlk/Xlo3LpiirKmJj/3kveYjI97FTzw90tE6V0r50dI3VBsYGFEZE3qCVHk/JG2W/xs4x72usRwoiMh47TH0mtvOmKEqcu5M3oSZip2qGY6eIumNHgNN4mynKqsYmwNhEV/g+9sWE7bGarM8prxNV2gGxMXmvJHmpxRjzO2d7F6yMZTusbOIbJjE+TmbJTsbKUH5ojHk0Se+LnbYtTLueZIx5oVwbOvSIVaFTtQJ8BfvQv7rQqUrqfQB7g3xzBetXlM7EP7G6l6OwOpYDsFqwo9yHl6KsRkzDar+2xOqyTsW+1HM9sJ0x5oN2bJtSRGIxcw3WomckcKyIjHSynQx8YYzZEBvC7PdJ2ZHY2azNsFrXa5P6wHbUHjHGbIIdpXQtgkrb0ZE1VsUURqxao7ESkb9ib4IRxnE5FpFbsMLs3saYhbVsq6IoiqIoKxcR2RkYZYw5MFm/AMAY89uiPI8meV5IjIhnYF8M+Elx3kI+7IsyrwLrm5Qdpg49YlUDCpqQ0BTFx1gfFk83oiiKoijKKscwlsfqBfvGpRslpCVPoo2bh32ruFzZ9bFvbP5DRF4Rkb/F9NkdXWO1onRP/oZ8nJY6eVoQkVOxI11IfY/tugxcbuhsAl3RvJNmMqWdWhPw1s27eZw6Qn684pTJSOXtoXoyEs8j4u4nRRkqty2UJ1aHTavctpBtsbefQJ7Wl4m3LQ2hY4yXWTnUYuzaVNHa0H6rq0ec9TYqY+JlovUGCrll8l6ewFXp5HGLuG1Nncet19tPFWVS7MdtXD708IyUkUAZ9xEmeXd7ijKR9XC9fh73C6EW9Yb3U7o6f+7LnxtjBgVytgkHidT8VcSJdvSo2N/rOmPMdUXroQeHe3bK5SmXnsW6E/zAGPOSiFyJHd26sFw7V/eOVcEUsgvWObiYrk6eFpIP6jqAbmtuazY45bmWbUu7+z5wjV1LPzc3T2NXv8zSbk6eLqV11HdrxqVr19K0hoa8s93fT0N9aZmuDX69DfWl5RqyTpn6QBknT0OdU0ddoEzGLVO6nnWfHkCDVM6Tqozj3ZcJfKtlaX0ZNy1dHqfeFCPLoXrbgnyKzkw+GIGjuA7/l4dbb2wdIOfUE8rjpuVKYp6nqzdnKq8H95OijJcnn6KMk6exua7idpvmtK3ZqSNXWkeaMgCNTZmKeXLN/rltbHLamyvN09gYaL+TJ5dz2h8q01Raps6po2GZXybbVDmPuz2Ulm0q3d6wNLCfxsplwnmk4naAhiWV6w2VcR63PHx3/Yd+a9qOz7HmdLVEYKkxplIIn+mUhv1ZC+txFsozPZkK7IP1KSxXdjow3RjzUpJ+F8m0YTlW96nAwgkNBQwexnK7f0VRFEVRakldprZLnPHACBFZT0QasGL00U6e0dgoB2Cd7p9KtFOjgWNEpIuIrAeMAMYZY2YA00Rk46TMvtiQSWVZ3UesxmOn9HbGun0XsxPwjgrXFUVRFGXVxxiTS/zFHsXaLdxgjHlTRH4FTDDGjMaGxPpnEpZqDrbzRZLvDmynKQecUeQm8APglqSz9j429FRZVpuOlYisg9VLvWeMKQyU3o8NwfF9Ebm1yMfqMGwAz7JzpAWMQK5++VRM3h9ZL9meZj2UVtfFmeLKBqarnDR3KjBbF5oKdPP49XrTetnKU4MQn/pzp/2CeSRexp3qi03zhdJi03yhNHf6za0D/Gm8VFOBNZj6y6QPa1eRvFT+BRicGnQ1MOJO6/lt86b+vDKBKRWnHncKL1Sve95C7c8404WufjB07mP7SfW5u/rH0EcY+UEe0jqC+0BKcW1404PVXE+BzyPv1uPmSXPQpXm8KgNlmt2pzcDz1qsh2lZ/Ki2fqWY/gSls5/S7TckHntG5BjfFrTfwfRGYHlypCFBX4zb4j2EPY8xDwENOWnH8x6XYOJ+hsr8BfhNIn4T1LktFh+5Yici3sPGUwL4O2SAiP0/WPzTG/LMo+83YoJjrYWPhYYyZJSIXYuNDPSEit2GnAM/Bxjm7os0PQlEURVE6HZJ2+q4VpOhZdQA6dMcKa+S1p5P2/5K/Y7GGhRUxxlwmIrOxsQKvwsYKvAP4iU4DKoqiKIpSSzp0x8oYs1ct8hpjbsRGBVcURVEUpa0RINvO05HtxOr+VqCiKIqiKMpKo0OPWHUIpFRsmEaI7ooPcwEherOT1tUVjDcEBNZOGVeI7grVISBED3pSVRaru0J1m1ZZeB70sYqI1d3tEBCve0L0QBlPvF55PVRP1lQWs4MvRM+a+Px/GjH7yvKtwmlvNT5WQQNHr4zjgeR5VMUF79mAf5brW5VG8O4Lz506Ar5obj1eHQFRufe5xrXSqV6GiJGJvJAQJOB1VZ2gPUZI8F75GMPbKwvcg5+7c2vm6t3toed6xaaRCZyifGQ/oba4ZfIBwbcneHee/fmgSL6dw9UJbaCxWjXQjpWiKDWjoTHLkC/6029Rb+py4Ydquse9Z8vf6nrC2yt3BKtxXq+mnjTO8mkc0uN1+N/+frgz/wdB7OXV0PbYMQXLuN9AztSR6VK5HcF6Q87lbtuqcC73ji60n1gZiDqvh9tfun3JsjrenNqdZ9/sxxcLvFcGOw61fitwFUE7Voqi1ISGxiybTF+XIf2G0G1gN+rq6xCRqsLGuARDqHgpaULNVFOm9Z0x75jThIRpZR3p9hMo44bkSdNJ8kLNBCqOheQJdAz9PO721ofOCe3HCwkT2Q5+mJtahbSppt7i/rExBvI51h42n5HDP+a6B4exYE6KHqiy0uic43SKotScIV/0Z0i/IfQa2ItsQxaJhL9RFKX1iAiZTD31XQewzrD+7LbZF+3dpDAi7eG83iHQEasIrkFoGo1VGoPQrJPm6aeyAV1Tg6uXihuEummuniqYx1nPBsQEMU1VSC8V01QFy0Q0Ve52mydiEBow2IxpqkL6qVQGoV6e1scKjNVZLdG4f6HtnnRo+cOu/8JedB/eLfCLuzQhTX/LH9GJjTWlGwlzp4RMilKtb0kgl3fQoWGiSJ7g5xHJEzo87wOK12uceoNH7JTxznUomHtkxDActDw2guiXMW5UePcREJ6zc7KkGHlN44kaIzTKVeaQ6xp6s9l6c3jwOX9byMxaWTlox0pRlJpQ15yhLqtPc0VZWUgmS7eGDmqa2YnF653zqBVFaRN0+k9RVh56v3VMdMRKURRFUZTao28FKiGq0lg5einXswp83ypfYxXQZbnaJ8+jKu5jFdRLxQIqhzypqgioHNNUhfRSDeRK1l39VFBjZVxdlnMOghqriI9V0G/KqTcQMTamh3LrSFOmrYhprsAPmJwveg1KCGtiPOWQlyWuUQrLvSprklIonzy81/IJaXbSqLtieQKtczVIMc0VBDRVaXRZkf2G2pem3hjBk19NvZFzGXzLccXbn+q0uYTFaBUzhVw1Kt5DEvasCnlbrVQK4vVOSOc8akVR2pUZn8xgUN2gZBnoLYOdxeXViZMYkhlQsqxRZtlq2Eiv/GMPPMLQTH/W9JZ+K+PwFYfb/vJnturRlcfvvae9m6IoK4yOWCmKotSAtTN9W5X/khv+zNdOPL6NWlNbxj70MLdf/zdenziRebPn0L1XL/oPHMimW2/F9rvvzjGnndreTVQ6IjoVqCiK0jZ8+smMlv+vseYa7diStuOsi85v+X9hmvKGK//M/HnzOOmH36V33z5F22Hk1lus7CZWxeU/v4jr/nAp9Q0N7HHQgayzwQYsW7qUj6d+yPNPPsWTox/QjpWiFKEdK0VR2pwt1l7eiZjVPIs11lyDWc2zkpQ0ruqlbLXd1szMz25VmWIOOOwgPs3PaSmzZqZ/K0qHOXvUBUVtscd010232o7Vj77L2sPX9bZ3dN57ezLXX3IZfQcM4NYxTzF8oxHA8vY3Nzfz/BNPtmcTlY6KANnOqTbqnEfdGhLxemHJ1xFYTMlSnD9Xb8iGlqy75J0lkKcutuSjS0O22V/qIksmsIizxLZLM10lV7rQVLI0kAsszSWLV8bkvKWru+SboktDvtlZciVL1+am6OKWCS1euZy/NDTnSpZQnlossf00NDd7S6XjF2PIBBZX0J7BlCwSWLw8hsBiSha/XgJLeD+lt7zbnsp1VOPxDxQAACAASURBVNpPobvnbi9u52E77MHIXkNZtmQJl/38V+w1Yis26jKAUT84m4wYfvvjn7J+XS9en/CyV8//3niTDet68ssfnG33KbQsixct5E//72IO3XoHtuw1kG36DuG4PffjsfvuL8kXbL8YRAyvvvQixhh2P2B/1tt4RFEZe3zZugx7HLh/yb7HPvQQm3ftxo1X/LGlHrvY7XusPYyvbr+t14bCuX/yvns4fo9d+dKgfuy1zlAuOPFbfPLh1GA733/7Lc76xpHsufYgdhvaj5MP3JuXxjzJ7ddfy/Z9G3ji/ruhaD/7bTCUb+yyNXPnzuEP55/JIVusz06DunLzn/7Y8mnN/WIOl190HofvsAk7rdmDvUcM4QfHHs7EF57BFP8TuO2Ga9h6SJbH/nM3RihZFi1eyJZDs5zxrcNL0i/51Y/ZYliWN16bwJ23XM+R+2/Ndhv2YM9t1mTU+afzxRefewfq1l2yjdB3EuQzpctKR8ROBdZyWUXQEStFUZR2Jp/Pc9JhX+f9d95ljwP2pd+A/gxbd+2q65s9axbH73MwU96azJY7bM/X/+9EmhobefqRx/nB147jnN/8ktPOPzdaT98BAwD48L33MMa0uW/Sf267lWcfe5T9v3IkX9p7b96YOJFH7rqD8U+P5eb/Ps2wdYe35H3n9dc4+cB9WLxwIXsdchjDN96Ej6a8yw+PPoxdDzio7D6WLV7CqV/eh1xTI7vufxBdu3dnyLC1AJgz6zNOOngPpn/wHlvtuDP7HX4Usz+bweP338ULTz7KqKv/zqFf/2ZNjvWvl/+Gl559ioMO/zp77Ptlxj83hrtuuZ6JLz7NLaOfp3ef1mn2lI6DdqwURVHamaVLlrBowQIeff1F+vRd8S/Ui773I6a8NZlfXHUZ3/zeaS3pSxYv5uRDvsoVF/2KA448gvWSqb1ybL/7rgxaYwivjRvPSQcfyhHHHcsWO2zPuiNGUFdXe5f9sQ89yHUPPsyOe+4F2IDR1//hd1zzq19wyXnncMXtd7fk/fUPv8eiBQv41V9v4NBjl78E8Ni99/CTE44tu4+PP/yA3Q86lN/+4za6dO0KLA+e/Mefn8P0D97juO+eyTm/vrSlzDfPOIsT9t+Fi8/9HjvvtT8DBg9Z4WN98Zknue3hF9lwxGYtaaPOO427b/kb11w2igt+dcUK76PdUbsFRVGUzsmW227FlttuxRbbbtVubTjv4lE16VR9Ov1jHrt3NDvttXtJpwqgW/funP3rUTQ3N/Pg7XdF6+rTrx/X3nMnG44cyYtP/ZcL/u9UDt1qW3YcOITvHHgwd/79Bpqamla4zQX2/PIhLZ2qAif+6GwGrzmMpx9+iDmzrC7vg3cm88aE8YzYfIuSThXAfl85is2337Hifs6++NKWTlWBxQsX8vh9d9GnX39O/8mokm0bbro5X/32KSxdsoSH7761uoNz+OpxJ7PBRqVWID88/9d07dqN0XfdTD7gi6esGuiIVYTC/HWBXMCEM2YI6pqBgm8A6gdYDphlOmkNWXc9ZBDqmn2GDEIdc0xx2hIy+6wqoHLOWU8TUNnJ4xh5NuRL6wTf3NNfD+wnXzlP0CDUScvm44Ga02zPBAxMVwZ5cX9n+cfjmogWr5c1CHW9Gd31UNDcFDNOvmli3FDTj1FsEx6f8FRRibiNaGw/Lq6WKx8Qr2+z47YV82TEkBFXD+b+3/Dqi+MwxpBrauKqX/3G28/ihYsBeG/yZHvuIyaiW+6wHQ9MGsdr4ybw0tixvD3pNV554UVeGjOWl8aM5e5/3MgNjzxI9569knYU6cqkTKVYHVfLESR/tt99d685DV0a2HLHHXnivnt5943X2HHvfZn82qsAbP2lnb38Amyzy268MWFci7aqeM99+vVnrfXXL2mPyQjvTn6d5lyOTbfdnm69e+Lehjvsvje3/vUqJr8+CSP+VZD6fYTCse68R6KRWt6OfgMHsv5GI3nrtYlM//B91llvQ9s+twpnX6GAyyHT0JWK0GlHrLRjpSiK0s50696dnr161aSuL2bPAWDicy8y8bkXy+ZbvHBR6jpFhK122oGtdtoBsJ3Yic89zwUnn8Jr4ydw/SWXceYvR61QuwEGDB4cTh9ip94Wzp+f/J2X5A9PyfUvU09xXS6FugcOCduBFNIXzJtXtu7WMGBQuI0DB9n2LVhQm/20GwXxeiekc3YnFUVROhCVROGZjH1M53L+6Oz8uXO9tF59egPw/Z+fz5Tmhd7ybs4u191/5wq1ebtdd+H8P/wOgJfGjF1+LIX2NvvtbW5uZvHChWXrnP3ZZ+H0mTMB6Nm7d8nf2Z/NDOafU6Ye28DwuW6pc2a4zs9nzijJB8s/m+bAsS6Y7382xcyeFW7j57Ps/nv16hPcrnR8tGOlKIrSgendz+quPp32sbft9YmveGnbfMmOKo1/9vm2bRjQIxllM0VT2n362bBAMwLtffeNN2hctqxsfROeecZLa2ps5LVx48hkMmy0xZYAbLyl1cJNevGFYD2TXngu5REsZ8ORW1CXzfLWKxNYssgfzZvw7BgANtlym5a0Xn3tsc78eJqX/81JEyvub8ILT3tpc2bP4v3/vUWv3n1Ya931W9P8jkldprbLKoJOBUYwYlodhDnrrgcCKruaqqwzH+6uh/PEAyxnMyaeJ6KpCuqlqgqoXDlPg/F/9bmaKlf75G4P5Ynpp4J58pXrAF8fFdJG+Xmc9RQC1bYKyhwLupzP+A+ySmUEE2yra4YZ0lT5heJlYjqs0GZPh5Volm6+7ia7HeHbp57glqq8o+C+3faXHkPG2Qb+MRbybLPjtgDcccPNfPmow1tGST56/wP+8rvLWspmxJA3wtrD1+WArxzGY/c9wA2XX8V3zjyjpYzNLHzwv3fp0q0ba669VigqdsuozvhnnuPTadM44MivlIq9BRqXLeOGP14OwA6779YyELTRFpvRpWtXHr/3Xs785UUtHa1FCxfx+x+fV3Lcy/VF9u/Yhx5k/NNjlr8ViOHGKy7js08+Zs8vH0L/wQMxBjbYdFNGbrsdb708kf/8+5YSAfuT993N6+NfajmMEi0Xyy2gijEYuvfswX5HHMWjd9/O3y75NT/8xW9btr8/+S3uvvl6unTrxkFHH9uijdps2+0B+M+d/+KYU79Pl3p7jubM+oyrL/758voDF+M9t/6do799ChuOWC5gv/p3F7J06RK+etzJSN1yxV1FVZ+UCcLcSafhOgLasVIUpdNz7unntPzf71i1L7vsvQdbbr8NYx99kq/svC877r4Ln82YyeP3P8TehxzIg3f4gYt/e93VTJ/6IRef+1Pu+PtNbLvzTvQbOICZn8zg3bfe5s2XJ3Hd/XfajlUFPvnoI358wv/Ro9eP2H63XVh/442ob2hg5iczeObRx5gzaxYjNhvJKectP389evbkmNNO4aYrr+arO36JfQ49lMalS3n28SdYb6ON6F3hzcc9DzmE0w47jP2PPJJh667LGy+/zItPPsmAwUM4/7LLSvJeeNU1nHzQflx06sk8df99rL/Jpnw45V3GPPgAexx8CE8//CDivZBRmbMuvpTXJ4zj5isv5bWXXmDrL+3K7M9m8sT9d9G4dCkXXnldiQZr7fU2YO9DvsJ/H7yPY/fegV32PoAF8+by9OMPsuNue/PB/94uu68v7bEvxx+8Mwcd/nX6DxzEuOfH8PrL4xi+/kacce4vW9XuDomgGitFUZS2Ym7TTM6/0BpS9ssOiS4/KuroFPjR6efQPzu4ZRlQV7oMrBvUsvz+l3/wyh93+PEMqhtYsgzO2KUjk8lkuOnBuzj6hOP46P2p3HzN9bz75mR+fc1l/PDn5wXL9BswgLufe5KfXnoxPXr14qE77+UfV17LuKefpW///lx01aVsv9vO0X3vfcjBXHLT39jv8EP59KNp3HvzLdzwxysZ89DDrLfRCM77w2+54/mnW0alCvz4dxdzxoU/Q0S4/fq/8dyTT3L48cdzzT13k6ngf3XYscfxh5tuZtp77/OvP/2Jt15+mQOP/hr/HDu2xBwUYNNttuGmJ8ey+0EHM/7pMdz65z8xe+ZMrr77fjbZ2k7X9ezduhcCBgwewk1PPsdx3zuTWZ9+wr+uuZwxD97PNl/ajT/f+xiHH+d3un/955s49tTvs2DeXO74x595dfwLnPTD87nwj3+tuK9Tz/oZ54y6hNcnjeef11/Jxx9N5ehvnsJN9z29mpiDik4FKoqiKLVl3AevAZUnEx8aPyaaZ8CggVz2jz976QZhan5+sEyXrl05+awfcPJZP/DKpKV3374ccfyxHH5cqeFmrI66ujrO+PnPOOPnP1teJpkPe+7j6V7+4777XY777ncpTHrtf+SRAOQjs7AjNt+cK++8t7RtBv79l2sBGL7RxiXbnpr6qc1TYQa+74CBnPWbSzi7yCAUfJuQAt26d+fHF1/Ojy++HMmXnpdJM315QzFf+/apfP2bpQGsJY3XiNKh0Y6VoiiKssrR1NTEgrlz6TdwUEn62If+w9MPP8hm223fEqpGaQc68VSgdqxiSKkpqGsGCnFD0JB4PevlccXsAZG5k+YagrpmoOCbf7rrwTyueD1gEOoJ3j1hekCIHhGrh4ToXh5XVB4q45iGxoTpoTzVmH9mm+N52lOsHiMmZgdf0F5SJgmM7JFk+elF5/LTi/z4dOGjlRJRtRG48i+XcuVflo8ixEZNjFPzbaP/5ZVx922c/ULAjNFvaWDfLmlMRWN54nvKuG0Njn64Dq0hg9BIngqC93KEjtid3Mm7577CMS+aN5cDRqzHTnvtzbojNkJEePeN1xk3dgxdu3fn/EsvX94kp16T4lr32pvmtvQMXAP7kUIb7OKb1oZMditft2GD0FhjlbZCO1aKoijKKkf3nj356onfYdyYMUx68QWWLVlCv0GDOfjrx3DyueezwciR8UqUtkOd1xVFURRl1aFL16787IqrPAuNdg7kEuXsUZdw9qhL2rsZShuiHStFURRFUWqPaqyUEEYgXzSaGZq3jhmCuvopiBuChgxCXU2VbxAaKBPRT0EKs89gQOXKmqqgQWhEUxU0CHW0TZ5+KqSxanbNPeNmn14ep46Q7ilNnpjGKo2eKhN7NapK8q4gx90e0KFU0lgJJqixiqqLgnoXU3E1WG80gcCrXSmMSJ163EDIYa1XawM5h8qsuC4r9BH7uqsUeqmY5iqQJ40xrKt1qsXXcBr5V5r9GPd4nJMpcblk+BNzL+2wZ23lPM62kHF14FG/chHptFOBnfOoFUVRFEVR2gAdsVIUpWYYYyoGFFYUpXYYYzq2qExHrBRFUVYAydGcSzE/oihKbcjnWJzTr/GOho5YKYpSE5rrZ7JwwRr07d+jvZuiKJ2C3LL5TPqke3s3I4wahCrlcIWBIYPQekc07ovXUxiEuuaf9QETS0+s7pQJmH+6Rp7uepo8Wc8SMS5WDxqERsTqrlDdprl5nPVmfz8xQ9CgeN0VvHsCeL9MGiF6yDS0tI644H1l4YrVY+J2W2b5r2VT9wGzZw9GzBr07NWVbDaDSEiyXEpI8J5G1BzVYKfQrnsGoalUzSnUxqmE560t03rBe9CU02lvVSaioc/MFatHXhQI1uNNI8fPrWuoGTLl9FuS5mJx6oi/TxH4CAMvc0jlzyj4XkBRLmMMmeZmcsvm89Gnc3hy6hod0yC0E4vXtWOlKEpNkLolNHcfx2fz1uXzOUMgnyX4LliKzkuqeHaeY3V8R7HuTTWO3Gna6nfg0pSpot4UPcOatL+a/aQpk6LzEqsn9LvET5PI9kCa2wcMdhSdPMEDqNSSdGWWNNbxyqfdeeqDNfh8SYN+kXcw9PNQFKVmSN0S6D6ZPJNb0nJ1pT+dvZGxQCcjViaUJ5dx1/1fy7E8OfF/5rt5GjOlj82c+PtplMp53O0AOUfy2khdxe02T2k9jaau4jpAzmSieRqd4Q5vvTlQptnN4+wnF2iLlydeprEp06p1gMZG53N299MYL9O8rHS9YZl/TWZzpWkNS/16G5x6sk1umUC9KfJ0SDrpVGDnHKdTFEVRFEVpA3TEKoohX6ShkvqAPidiCBrUWMUMQgNmn54OK6K5CqUFDULdoMuuQWiagMqu5qqKgMqufiqU5mqqXDNQ8DVVDTnHVDRYJl8xT0gLlcog1NVhOfWkMwhd8TftXGPPVGWCBqGuDqtCUOYE9zylGY1yz0uaPJk699z6x+yVcfJkAvddJvL7MxMY5XKLZIgLXjKBka8YeVf/WMUAQT5QKDaqGNLfhbVarWyLcdcD+8m7607bAvdq3nkGu3lCn7v7HHf3m8sHri+3LQGjZ9fM030kB/VSbh73PgzsJ9/eo0UaK1BRFEVRFKVGiOhUoKIoiqIoirJi6IiVoiiKoii1R6cClRCuj1VIquLOxbt5XG1UKC22Dr5eylsPeFS5+qmgj5Wj0/D1U/EgzFnj1BHUS1XWVAXLRDRVQR8rN08ajVVzZY1VqEwavVRcYxXyx1o5PlbxIMz+xV6Nxsot457rXODhmyaPr5dydTOBMp4Oq85ZD5179zFZej25vlAhvDxpvm/Ev7bjlLY1+Hk4mqqQxqot5jKCeinnOefmyQfuD09TZdw64vt2NVZBXZazH/dyag6cI08/FZBHZpqdfTvTZbn6UBmnbXVu+0P3XUeOdbN6ox0rRVEURVFqSyd2Xu+c43SKoiiKoihtgI5YKYqiKIpSYzSkjaIoiqIoSm3oxFOB2rGKIVDsA9clGxB/xwxCA+ZtvkFoXLweMwRNE4TZNf8MpWVdYXogCHPWC6hceR2qC6gcE6u7wvRQWhohuidez7VevJ5K4J7C7DMmXq/GMDSNQWiqoMuuWN0Ts6cIT+PUkQkYLXp5Upl9phCvV2Eq6tfhODiGvD/dj8gJg0PgRY2YMCMfePGkKqoxEa3CINS7NkLiddcQNMX1FBW8B4XoThn34wicWt+INE2ZyvcH+MLzmGGoLVO67grcM82h/fj1KCsH7VgpiqIoilJbOrHzeuc8akVRFEVRlDZAR6wURVEURakxKl5XyiGlgZdDAZXdAJ6efiqgy4qWSRFQ2dVPBQ1CXf1UVQahAR2Tq6lKYfbp6bJSBFSOaaqCGqum1mus3DJp9FOuDiukjYoahKYwA61FEOYQMd1VUB8SMwgNlHGv21y2tEw257cjl3WMO4MBbyufy0yNTEVjdQTzxMwZXc0V+LqrVCaijiGoZ/7pX7denjQmoq7GKhS42b023DqCGis3mLBzrYQ+d+fcuuarwUDabhk3SHbgXMfyhMq4pqGufgp809CYYSiU6nxtHnc9EOy5wW/fSkWAFLrN1ZHO2Z1UFEVRFEVpA3TESlEURVGU2tNJpwI751EriqIoiqK0ATpilYLiufTQ/H3MxyoUpNXzrfICLAfKRHyrgvopV4cV8KRyNVSe5irkSWVyldeDAZVdjVU8oHJMU+Vqo0J5XC1UqIynw0qhn6qNxqoDB2EOiEiq0Vi5aZl8ZV8r8INTu2VsWkQvVSPvK6+OgJ7Fp3Lg5nDFjl7Kue+CwZLdemskaalKY+Xmqcb7ylkPaU29azDiawWhoMuVfa1sntath/OkuIfcgMqBZ7/nW+XuJ6TLCvhhrVTUIFRRFEVRFKVWdN63AjvnUSuKoiiKorQBOmKlKIqiKEpt6cRTgTpipSiKoiiKUiN0xCqCiCkRo6cxCHU1msEgzJ5BaGUxO8QNQdOI10Nmn75Y3RF/m0C9Tlo2X3kdQoagrQ+o7ArPUxmEphGvu0GXqxKvB85TrrKCNCh4bw6oYiOIU4+pwpgvX5dGiB4LwhwX67rHHNqPK1bPBO4794WDRlOq8K2VqahXRwqDUJ8qxOzejuP1uoagqcw/AyaiDYFyleoI7StVEGYnj2sIGgpG7wveK28HX1Sec4IW18og1E1zDUPBF6d75yBoKuqe28p1QgcJwpziZZDVkc551IqiKIqiKG2AjlgpiqIoilJbRFRjpSiKoiiKUhMEa7dQyyXNbkUOEpF3RGSKiPwksL2LiNyebH9JRIYXbbsgSX9HRA4sSp8qIq+LyCQRmRBrg45YRRBSGIRGgi6HdAK+pmrFgzC7AZfB10+FDUKbK+YJmn1Ggi4HzT6bKwdHTqOXSmUQ2tRUup+cczyNpdtDeWKaK/A1VdkmP4+rffIEISFiQZeb43V4vxNT/HKsc/QQQYmGowdxtVy5er+Uq7tyz6UbcBnSBauOGZxWYyoafiRWoYeK4u8n796/aaR2rtbGM/JMEYQ5pJdydu7W6z4zwvtOYRDq5HGfcUG9VMQQNBe4P1ypj/vMdg1DwTcNrc4g1M/j6sj8az1Qr6vLcu7nkJ6qjWK3d2hEpA64BtgfmA6MF5HRxpi3irKdDHxhjNlQRI4Bfg98Q0RGAscAmwFrAk+IyEbGtDj17m2M+TxNO3TESlEURVGU2lMntV3i7AhMMca8b4xpBP4NHOHkOQK4Kfn/XcC+IiJJ+r+NMcuMMR8AU5L6Wo12rBRFURRFWRUYKCITipZTne3DgGlF69OTtGAeY0wOmAcMiJQ1wGMiMjGwTw+dClQURVEUpbaItIXdwufGmO0r7TWQ5s4Nl8tTqeyuxphPRGQw8LiITDbGPF2uEdqxiiGlmqmwd4nrY+UGdq0mCHONfKxSaKyyTrBXz6Mq5GMV8a0K+VjFNFXu9lAeT3PV5OulGhojPlaNrfexSqWfCnlWuXoPT4SRQnNVha9VKmJi0JCGybkBxBmer8/5bTVZx/uqrrLGJJRWlcYq5DkX1VgF6gkEc14pOLqZoN7ISfP0kKGmO9WENVaVA0CH/LHc50+aMm5aKo2Ve9tJFd5XEV8rqI2PVSiP620V87UCX0Pla678/bR7EGZoj7cCpwNrF62vBXxSJs90EckCfYA5lcoaYwp/PxORe7FThGU7VjoVqCiKoijK6sB4YISIrCciDVgx+mgnz2jghOT/RwNPGWNMkn5M8tbgesAIYJyI9BCRXgAi0gM4AHijUiN0xEpRFEVRlNpSsFtYiRhjciLyfeBR7JjvDcaYN0XkV8AEY8xo4O/AP0VkCnak6pik7JsicgfwFvZV4DOMMc0iMgS41+rbyQK3GmMeqdSODt2xEpEMcCZwGjAcmAXcAVxkjFmUonxP4IfAsUn5ZcD/gOuAm5JeqqIoiqIoqwHGmIeAh5y0i4r+vxT4WpmyvwF+46S9D2zVmjZ06I4VcDm2Y3QvcBmwabK+jYjsZ0xA/JOQdMoeBnbBvlp5NdAd28n6R1LX+W3aekVRFEXprFQRs3R1oMN2rERkM+AHwD3GmKOK0j8ArsIO391aoYqdgN2AK4wxZxWVvxaYjB0Fi3asXINQ1/wTAgah7npAvO4J3N1170WGFOL1oPlnLponZgjqittDeWIBliEgRE9l9llZrB4SotdCvF7nitXTCNODeSJi9TTC9DQC92qIPfRCw/huGTdPoE5xlLWuEWko6LRrvhoKcB0yFi0tE7jvAoL2OCseQDkk3G51HQFJrGsq6gnGA7893TwhY+FYoOaQ4L0a8bonPE8hXo8Fag4L3kvz5JrdOgIvRzjXnGsYmsuFBO/xF5dyUbPP0IsakTKhYM/tHYRZZKVPBXYUOvJRH4vt11zhpF8PLAa+GSnfO/lb8kZAYhr2ORCdSlQURVEURWkNHXbECtgBG9RhXHGiMWapiExKtldiHDAXOE9EpgIvAd2AE4HtgNNr3F5FURRFUQroVGCHY02sGdiywLaPgV1EpCEZgfIwxnwhIocDf8MK3gssAI4yxtxX8xYriqIoitKp6cgdq+7Yt/hCLC3KE+xYJSzE+k2MBp4H+gNnALeKyBHGmMdDhRLL+lMBMmsMK5knT2UQ6gZUDgVhjgRdDhmExnRYafRT2UDw1JghaCgIc9Zxn/PWA2afnomoaxga0Fj5xp1uQOV4mZjmCqBumVNPKo2Vc74DJqJRTdWqZhAa01SF6nSP0SnjarAA6lMYhMZMRMNlKp/LTPn3YYrqaBsFhRfE2A18nEanlcZUNKKfAsg62ie/TNzAuCpTUfccBK5BX4dVut3VbYH/7PR0WaFgyY5pqBfIOaDXyzuNyeUC16AbADqVQWikTMCI023LSqcd7BY6Ch25Y7UYGFxmW9eiPEFEZAtsZ+osY8xfitJvw3a2rheRDYoiV7dgjLkOa8lAdtOt1JJBURRFUZRUdOTu5CfYgItdAtuGYacJK41WnYXtgN1ZnGiMWQw8CKyL9bZSFEVRFKWmiB2druWyitCRO1bjse3bsThRRLoCWwMTIuULUalDL51mnb+KoiiKotSKwlRgLZdVhI7csbgd+CnwI+CZovRTsNqqWwoJIrIBUG+MmVyU7y1sTJ8TgT8U5e0LHAF8AbwXa4T1sSrSWIk/M+gFVHbWQ2W8OX9XYxUKqBzxscqk0ViFAip7QZhd/VQ8CHPMowpCnlQpgjC7ZRpLfaxCeilXU9WwrLSMBAIF4+Tx9FJp9FOhPOpj5R+jWyZ0fI5HWF1IL+XUm8b7Cuqd9dLPPZPK/Ke0jnxAD1mN15VL3tU5hfymXE1S3tFlBfRTroYnpH1y/e9S6bKI6LLSaKyctuUCAbDdZ2XM1wp8rVnOub5C2lnX28rVS4U8qjx9bQodlhtH3tVP2TTX6ypexou2raw0OmzHyhjzuohcA3xfRO7BWtQXnNfHUmoO+iR2aq/4SroC+Dbwu0Rv9RxWvH4KMBQbB2jFn36KoiiKovisQtN3taTDdqwSfgRMxb6hdwjW2PNqbKzAij/1jTEfisiOwEXAvlin9iXAJOAcY8w9bdhuRVEURVE6IR26Y5W8sXdZslTKN7xM+nvACbVvmaIoiqIoZenEIW06dMdKURRFUZRVFJ0KVIJIPAizK3z0jDwDwkJXgOmafQZF8pGgyyGDUM9ENChed+pxjTzzIbPPCb+S1gAAIABJREFUygL3oEFos2si6gjgUxiEesL0UBDmmFjdFapDXKweEqZXI3B3hdxuIOcQEVPLqgkpdosJmA62hUFoUMBf74r+/TyusWgaU1GXdGL1aih9tOY9kXZAYO2+EBIIZO6ViQjcQwahDU5aPiAQd4XnMfNP8J8/rlg9F3gJPVZv6CUeV8gdMwyFQLBnR+DuBmWGkBDdLVNdEGY3TZxrPZ/zr8mqDEI752BRh0A7VoqiKIqi1JZO7LzeOY9aURRFURSlDdARK0VRFEVRasyq5ZZeS7RjFcU4BqF+Dk9T5eRJE1A5jUFoLOhyMAhzxPwzTZ6QQahrABrTTwXzeAGW42afacpENVVuwGWojcYqpBXy8qziBqGulsPVaTUENEvuMdZ7kYL9Mm5amjyuqahfgoZAWmvJe9q0gGavmnojQZiDRp7NkTKBZ4Krw8qm0lg1V9wOfoB3r45gWyobgoaenbE8rmEoxDVVIVNRLwiz+1yv86/JnCtBDHnsRvKkMwitrLmyedq5U2Pdtdu3De1E5zxqRVEURVGUNkBHrBRFURRFqT3tPWrWTuiIlaIoiqIoSo3QEasIXhDmNIE3vSDMfr2ej5Wr08LfTyzoclBjlSYIs+dblcbHqrJvVaqAyl6A5db7WNWF9FIxTVVIL7W0iiDMMf0UVBeEuRodVjW4Gir312XomGM+ViHPLU9EEvG1goDfV+g8OWld/SwubjDnWmiu/MDOAT+pNHqpXOVf96EyblpMcwV+oOagVxSVdViuLxcEfKzcOqT1flnBIMxuHtejKqTLco/Z1VzlA5oxNwizF7g58Ix2va4Cn2nM66oxFLi5ydVUVQ7KbPO0kTYzLSKqsVIURVEURVFWDB2xUhRFURSl9qjdgqIoiqIoSg0QVLyuKIqiKIqirBg6YhXD0d8FxeueQWg8oHJMrF6NQWhQ8O4GWA4GYXYNQeMGoVUFVK7CINRNq6tFcGRXqJ6mTGg/jSnyxIToIePLmFg9TeBml9AvR7e9MWE6gGuK6NYbOp6YwN01DK0VbiDnAO6e04jZXYPQfMb/3D1xcWQ9mMcVpoeE6G5b3EDtgXu3KoNQU1mYDr44PY1BqBuY2Qs0H3gOxsTqQYNQ57md8wyaQ2aflV9KCgdhdtfjLzt56wHjUVecvkoYhIKK1xVFURRFUZQVQ0esFEVRFEWpLSLBkdnOgHasFEVRFEWpKYZQXM3OgXasIgileqigXio2Zx4qE9FUhcrENFWhAMsZ45apwiA0YPaZMaZinlCZrBMcOVUQZtfs080TKBM1BA0ahFZRxtVYtVUQ5mo0VdXgajLSBGFOo8tyAzNHgicH06oN1NxK6gLtr42JaJyYXiofML507zNXl5UN6I28PCGDUDcovKm8H/CfT55+KoUuK6a5CqV5ZYIGoZWDLod0WZ4ONsWl7pbJBsw+nfj10aDM4GuoYkGZQ2WUlYd2rBRFURRFqTmddSqwc47TKYqiKIqitAE6YqUoiqIoSk0xIjSH5ASdAO1YpaA0CHNgey18rFIEYa6Jj1XA18bVS3m6LBOoN6KpygT0LjFNVUhjJY4uqypPqliA5VCau+7qqdLm8YIJpwnCnMLrqhbEhumDGqtI4GZXTwV++13fqlTnoI0CUafAPaKs50kV0DFV4WPler95gZsDZbKeLsvRT2X8z8PzsQq137g+Vo4WKvRMiARddrVQtkzlQPKuZ5Wt1ynj+liFAjd7XlelxxfUtHq+VY5OKxvQcjneVlX5WAXKNLlxzCOaKwgHZl7Z6FSgoiiKoiiKskLoiJWiKIqiKLVFwHRSu4XOedSKoiiKoihtgI5YKYqiKIpSU6xBaOfUWGnHKoaExYTFuNdO2wVhjhmExgM3B4XokaDLQbPPSBDmsEFoRLxejai8qiDMfoBoL8+SFIJ3zyA0cJ00OvvyzDFDwu0qxOpuPdW8jZMmoHJsP6G2u4J2t45aBWFOY7YaI4Upan2KauLi9ZBxp3M/OG0JCt49gbsjsM77161vIhowCHVF5a5gPGBG7NUbEaYD5CLPvWAZR9DumXKmCNzsnkpXqA6+WD2V8bMnRPeypBKruxhXrF4XNwhV2g/tWCmKoiiKUls0VqCiKIqiKEpt6MyxAjvnUSuKoiiKorQBOmIVwQZhXr6eal49orkKpVVjEBpbh0Cw5GCgZqeeVBorR+/llgmYfbqmoW4ezwwU/CDLtdBYpdFLpTH/bKsgzDGNVSotUaC9LjEdlqu5CpVx84T0Um57Pc1VigDLbRWIugZTFa5hKEBDxjH7dI08MwGNklNPQ1MaXVblezVo/pkiUHPO1Wo5zw3XMBR809BMCr2UayrqmogGTY/det2gzAFTUU+HFTEMhbimKmwWXXm/oXL+eqCMo6HyAywHDEI7QBDmzjoVqCNWiqIoiqIoNUJHrBRFaXNyecPMJU18uriJBcl6U96QyxuaRchmhPqMsO2gHgzunua9O0VROjJGhObA6GFnQDtWiqKsMItyeSbNWcJmA7vTt0vpY2XLu9/ijS+WBCZ1fB44dGMOXa9fSdrBd7/JvKY8Q3s0sGbPBob2bGCd3l3YZp0+bDKwB3WddLpBUTo6nXUqUDtWiqK0ioVNzUyau5SJc5Yw8YslTJyzhMnzl5E3cN8+63HEun1L8tdlSNWpgrBeafyMhcwOeY8B3eszbL1GL7Zfsxfbrd2HQzYZyIDuDa09JEVRlJqhHatWEhQWBgSKse1VCdG9PI6APGAQ6opOqzEIDZWJCdxD4nU3zRWzB0XlrvC5GiF6bD1UTy3KgN/+asTrtTC+DOK0t4wwfXEuz5OzFvHAjIU8N2cpby9YVrajNHHmIo4Y0rMkbWjXLJOAQV2zrNktS98uWeozQr0I2Qxk6jLkkmnBQdlMyXlc1pwv26kCWNyU5/lp83h+2jx4aTqTTt+BAWuUiuONMUhBsF2VcaqrLi7t/ElAxZytKz23DY2lhrOhX/IxU9HQPeUJ0R1Ty2yd/zJB3r13AwahblrOFX+neNa4Ane3DgiYHEeMSQFyprLAPfS89cXqlU1GQ2muiahrIArpzD9jecKmoqXrbq2B9w86hGmojlgpiqIEOHbCJ4yesbBinozAJr260LfBf8Lfsvu69MhmaCh0UNyHbUP5x1B9Rnj3W1vzSWOOTxY28umiRj5Z2Mg7c5Yw8bNFfLKwsSVvl7oMIwf1KCk/b2kTW1/9IgeOGMDhmw5in5FD6Forp3dFUZQA2rFSFAWAtxcs46MlOQ4cXNo5OXhIj5KOVQbYtHcXtuvble36d2O7/t3Yum9XemQzwRGhfl2qf8xkRNiwb1c2dC0aAOrr+HThMibOWMjEGQtZ0Jyn3tn/I/+bzdQvlvLXcR/z13Ef072hjgM2Gsjhmw/miM2G0L9vt6rbpihKeYyA6aQGodqxUpROTGPecM+MhVz74Vye+WIpa3fN8uF+6y+fOgMOW6MnV73/BYev0ZMvD+3Fdv262U4UhL2uViJDe3bh0A27cOiGA6Cr/zh7fMrskvXFjc3c98ZM7ntjJl2yb3LMtmvyvd2Hs+Pwfl5ZRVFWBA1po5TFlMyBx/RUoTxtZhDqGvKlMAgN6qUiBqGeFiqQFlsHXyNS52qUAlqPqKYqFCi4FgahMcPQNGXA10elCdwcq6OaPM5IzkdLc1w3fT7Xf7KAz4raNG1pjldnL2Hr3l1a0oYhvLXLOkk9Ys9F4Xy4I0khXUfMILQN+cvBG3HCFmsw+p3PuX/yLN6ds6Rl27JcnpvGTeemcdPZbq3e/OLAERy27TC/kojGKiSKcd9UdEX5Qb1UJFCzq9sCyGVLz2WaYOie+WcgULNrGuoajbpGnhAPCp+VwH5w660caB7iz86gLss1HvXq8IpE81Sjn0qTJ02ZJuf0h8xAO2unpiOgHStF6STkjeHxOUu4dtp8/vP5Ys8HOytw+KAetTAj7zBk6zLsvm5fdl+3L5ccsCHvzFvKA5NncdurM3j5kwUt+SZOn8/MBY0ValIUpVVI540VqB0rRekE3PvZIn46ZQ6TFzd524Z1qeO0Yb35v3X6MHQF9FCrAhsP6sHGg3pw7u7DGT99HteO/5jbXvmUrtkMx24z1Mv/ydylrLmGGpYqipKe1fspqigKAI/MXux1qvbt343vrdWbwwd2t1NV1VgRrMLssFYf/rHhQC49fFNe/WQ+PZxO5csfzWXHS57h2zuvw6gjRrLOgO7t1FJFWfUw+HYgnQXtWNWAWuilUu0n5n2VRj8V8J+J6bDCQZgr+1aFdVnOvmMBlkNpabyjqvGkqkUQ5qo0Vm0VhLmUXwztxT8/XUhW4KTBPTh9jZ5s0ivRUbUEv3Z9rVIEYXbbGtJPVSNwD2nnWkua/WYyDKivY591+9nPpmge9IL73qI5b/jHcx9y60vTOGOv9fjpQRszwO1ghTqkdRGNVTaksXL0RtkUwdCdeyjmawW+t5UbPBl8HytvPaBj8oMhp3g+RZ6V4cDNKx6EOY0O1vOtqsL7Kqjdivx+CWmsXIyTJ9/OL5EopXSun6iKspozZUkTJ/xvdokgHWDNLnXctfEAPthuKFes349NNB5fRRY3NpMp6qQsy+X54xPvsf7PH+PXo99mYQXTUkVRLPmM1HRZVdCOlaKsBizNG34ydS6bvjKDm2ct5jfT53t5vty/GwPUHDMV3RvqePj0HfnvGV8qsWKYvzTHhfe+yYbnP8zdE6a3YwsVpWNjRMhnMjVdVhVWnZYqihLkpQXL2HbSDH7/8QJyyQzBX2Ys9EatlNaz14gBvHj+Htxz2o5sssbyUD0z5y/j6Gte5Ni/vMTn85e1YwsVReloaMdKUVZRluYN50+dyy6vfcbbS5ZPTe3euwtjtxjM4JXoF7U6IyIcuc2avH7hPvz9W9swtG/Xlm3/fmka37jquXZsnaJ0XJpFarqsKqh4PYJIOjHhilJdEOY0ZVofUDmVqWjEEDRkgOgJ2l3hc0iw7KZVE7jZ229A/N3oaGaqEbyvRPH6Swsb+c4Hc3m7SOvTMyNcsnZvThvU3Tqnx0asYkFag+J1p4wj8A0eT6yDFzreWnQK02gyygSeDuXJAiftuBZHbrsmZ9/5Ojc+/xF1GeH3X9vCvw4jgZrDBqFunkzF7eCbhsYMQ22a8+JJoN5cprSeNM+EWBDmbEBU7gnRnfVcGoPQNEbJVYjXo2WCRp7ueop6vTJeEd9ENEWA5Y4QhLmzoh0rRVmFWJo3/GLafC6dsbDkK2ifXg38fb2+DF/Nfag6Av16NPCPE7fj6G2HMXnGArZfv397N0lROhwGNQhVFGUV4MpPF/CHooDI3iiVstI4ZMs1OGTLNbz0h1/9lIffnMll39qW+mzn/GJRFBBMJ30macdKUVYhzh7ai0fnLuW/Cxp1lKoD8s6n8zn2zy8yb3ETb0ybx51n7caAXl3iBRVFWW3QJ3IrSRVUs5r5e08vFTDyjAVhDugevDqqCsIcMhWtnCdUJhp0OY1BaBqNUkz7lMZUNKaNqrYtaerNlf8c64E71+7DffOXcVK/rggCy6p8+y+mwcim0VhJ5e21aEc5YhqqNKaJrdBYld1v0fa/PfUe8xKH+/++OZMdLniE0WfvxuYbDy4pkg1c677GyjHyzKbQSznarVzAvNQzAA7cq25g5lhQZoCco6lKpdV0zD2reg7WQi8VLNO6OsqleXlSBF32y1Te3iGDMEsHaEM7oePUitJByRnDP2YvJu98IQ3IZji5fzed+uuA/OFrW/Crr4xsWf9g1iJ2/uWTjB43rR1bpSidBxE5SETeEZEpIvKTwPYuInJ7sv0lERletO2CJP0dETnQKVcnIq+IyH9ibdCOlaJ0QObk8hz83hxO+mge/69IU6V0bESECw/blHvO3LUl9uDCpTm+8vsxXHzX65gUo8qKsjpgYwVmarrEEJE64BrgYGAkcKyIjHSynQx8YYzZELgc+H1SdiRwDLAZcBBwbVJfgTOBt9Mcu3asFKWDMWVZjp3+9zlPLGgEYNSMhTy3sLGdW6W0hiN3WIsXRu3L8EE9ADAGfnbrJI67/FkaQ1PeirIa0g4hbXYEphhj3jfGNAL/Bo5w8hwB3JT8/y5gX7HD/0cA/zbGLDPGfABMSepDRNYCDgH+lqYRqrFazUmlc4jkCQdUrpwnVCbqSRX0dHI9qVLosmL1pvGOSuU3VQMdlqOnentpE/u+O4dPc8vzjRrUnZ3rZbnXVsznKi2uB5VLIICvp4dydVgrU2MV01BVo7FKEzW3gsaqOM8Wa/Ri/IX78LVrX2TM5FkA/PvZqSxY3MhdP9mHro5XV6be8XTKufqpgK7J861yyjT7fmBunlxAwBN9JqR5jqTx2XPLVBHAviYegNVoo4KXSjX6qRXXXIVu1dWUgSIyoWj9OmPMdUXrw4DieffpwE5OHS15jDE5EZkHDEjSX3TKDkv+fwVwHtArTSM7z8ehKB2c15Y0sWdRp6qbwF3r9uEXg3uUBARWVh0G9urCY+fszun7b9iS9uDLn3DPCx+2Y6sUZSUgQr7GC/C5MWb7ouU6d6+Blrg91XJ5gukicijwmTFmYtpD146VonQAXl7cxN7vzmZW0qnqmREeXq8fR/XpGimpdHTqsxmuPXkHfnbkZgBceNTmHLfn+u3cKkVZLZkOrF20vhbwSbk8IpIF+gBzKpTdFThcRKZipxb3EZF/VWqETgUqSjvz+pIm9p8ymznJFF/vjPDIen3ZuUdDO7dMqRUiwq+P2Yq9NxvCPpsPae/mKEqbY4Dmle+8Ph4YISLrAR9jxejHOXlGAycALwBHA08ZY4yIjAZuFZE/AmsCI4BxxpgXgAsARGQv4FxjzDcrNUI7VorSjkxemmPfKXNaOlX96oTH1+vHdt3r27llSluw7xa+U3tTLk/eGDLd9DNXlBUh0Ux9H3gUqANuMMa8KSK/AiYYY0YDfwf+KSJTsCNVxyRl3xSRO4C3gBxwhjGmqjdNtGPVQaiNaDNk5OmKTiOCZcoIzyN5/PXAfty0WLDkUJor3A4JuWOi8pDg3ROZR8TsobQ0QZgdM88ZSxpZlOTpnREeW7cP22UzpXWFhPPefqoIqBx7ZIRE5V7QZecXachkPPZCY0xED+H2R4NMp/i1XI1BqGeKGthPvSMa94I02zKNTc0cc8WzNOfN/2fvzMPkqsr8/32rqqs7G9khC0uAJKwqSEDHBXFhYERkHHVExRlXVETH7ec6oqK4jIIroriNGyqjOEYHFxQRRtkCAwjIEkiELJAVsnSqq6vq/P6oSqfqPW/XOX371tKp7+d58qTvvefce7q6uvv0uZ/7ffGTDz0XubpzeWJ6yQjy9AJBlRBvhX9qwT0bFtx1IKgODDXb6CLMzuijHqHXRZjTEtHTKMLsHY8psGydNyCnR4VQq23rO8EKDW03lQ64oc65KwFcqfadV/dxAcBLR+l7AYALmpz7GgDXhMZAx4qQDnLSlDx+c9AMLMhl8KuDpmMZVy16huFSBf/02Wvxs5sexvIVa/Cub9zU6SERkhrVHKvU5fUJASdWhHSYZ0zpwwNLZuFpvP3XU+SygiP3nz6y/aVf/BXf/O19HRwRISQNOLEipM0MGbc8B3q0plYvIyL41CuPxYufcuDIvnO+cj3+dPejHRwVISkhApfJpPpvokDHKoKY8Lix9rfcgU7hFVDWBZZjijt7RZgt9ymFgNCYsM9QKGcSLyvJ2ADPqbps8yDO3ziI5Qfug6X9uVG8rMDrHeNceX0i2mhXqGx0SlJ0OdQm7zs+UeGrwetaAadqnz6vNcHVfbQ/ZX09PKdKnbd2jgyA75zzVKx8ZBtu/9tjGC5V8OILrsbNnz0N8w+c2dBFO1eAESKaCweEltTXwy6ynkJAaJI+qiizdq6AiCLM1nVaULg5Bmsu4BV3jpgvhLwsZx7nH2udYuJMAQmZ4Nw8OIzXrduOe4tlnPDgY7hxcLjTQyJdwJSBHH7+/56FubXMsg2PF3DGJ/6AnQW+P8jEho4VIaRlrB8u4x9XP4ZC7Q/L+bkMDu83VmpIT3LQ3Cn46XtOQq62Injbqi14/WevZdFmMmGhvE4IaRkl5/Divz2OdbXH5WdkBMsP3AfTY+IASM/wzKP2w8VnP3Vk+6fXrcJnL7+jgyMihCSBjtUYadW9+ERjiXChYjyH4HUsB8PKqRorMTlWIS8rNZcrVLjZuI72c4bsHKvPbBzE9bXbfhkAl8+biiUi1f4xvldM0eUk3lXIqTKLGOvJoHahjC46QH4oPLREhPwpwCigrLbN91Og8HdM9lUo16rW5uyTF+OO1Vtw8a/uBQBc+F934JzTDse0yXnPpwJ87yomTy7KffK8yxQcq5hMqkBR5mqf5llXrcqxSiOjysI7bwrOFQBUumBBfCKtMqUJJ1aEtJA7CyV8ZOPOke2PzZ6Ek6cwVoGMzudeswwrVm7G5Ml9+Po7nolpk1naiJCJBCdWhLSIknN4zdrtKNb+uDxhUg7vmcmiyqQ5fbkMrvz352CfuVMhPfoXP5n4OBFUpDd1h67+rEUkIyLvEJF7RKQgIg+LyIUiMmUM55glIp8VkZW1c2wUkT+IyDNbOXZCPrNxECsKJQBAXoBvL5iGHH9RkghmTevnpIpMeHpVXu/2FavPAXgbgJ8BuBDAEbXtY0Xkec41L3wnIgehWtdnKqqFF+8DMB3AEwEsbN2wSa9zd6GEjzy6Y2T7o3On4MiBUXKrCIng7tVbsGT/6egz8qwIId1D106sROQoAG8FcIVz7sV1+1cB+CKqFakvC5zm+6h+jk90zq1PY1wV58+arX2twCqynMp5EwntgSLMlkytRdqIYs/Bws1RBYoTBITGXKc0uvB+cFbwtlmTcNHmXVg2kMO7Zw5UrxuS5K3x6TZGMd5E6PPktMxu9NFhn1rctnQgLbQnkdnNItIBWd3q44nnusByRKJjSEwHwsK72afxtc3W+gyXKvjEf92Bj19+Oz74quPw/lc+ec/QArK6LqYMGKGixuesCzOHZHZrX5zwrs6bSjH6lPp4snrz43aflAT3rr63ZOMAlHu0okQ3f7lejmp07OfV/q8DGARwVrPOInIigGcA+A/n3HoR6RORyS0ZKSGKSRnBZ+ZNxf8eNB3/OX8qbwGSxHz36pX4yA9vQ6ns8Inv34o7V23p9JAIIU3o5onV8ag+y91Q8t05VwBwW+14M55f+/8hEfkFgF0AdorIfSLSdFJGSFr83aQ+HNHftQvDZALw6ucuxlMPmwugunr1ga/f2OERERJBTV5P899EoZtHugDAJuecdaNgLYA5ItLsOeTDav9/HcAsAP8K4HWo3pj4noi8Js3BEkJIK8hmM/j6uU9DpnZb5bc3P4w/3r6uw6MiJIwTSfXfRKGb/5SejNHti0JdGyuSEACm1f7fDuDZzrkiAIjIzwA8COATIvIdS4AXkbMBnA0AuQULxuxQtcu5ahVmAeV2EOVLJQjLTOQ1BbaN8zy4q4QLt+zCv8+ZjPm7faUkBaGHSv6+UJ/A2EwSBYRqIkI59Z8/IecKAMrq2jHSvx6v1cf7nCNculDArBX2qQtLRzl7o3tYRx8wA2c9dzG+e9X9AIAPfeNGXPf5F/q+VKnx87MKN3vuk/H97rlagULt1j4/yDOioHJEcfp0fKnOlQrSjlVoO+qcMcXQSdvo5hWrQQD9oxwbqGszGrtq//9w96QKAJxzWwEsBzAPe1a1GnDOXeqcW+acW5abNWtsoyY9y4c3DeIrjxVw6ANb8K3HCuEOhIyBD73qycj3VX9k33TPRvz8z3/r8IgIGR0HoAJJ9d9EoZsnVutQvd1nTa4WonqbcLTVKgBYU/v/EePY7icEZ45jfISMcMeuYfxgW3WBdZcDDtMrFoSMk4P2m4Y3n37kyPaHvn0zSklKGRFCWko3T6xuRnV8J9TvFJEBAMcAWBHov1t63984tnvfhvEMkJDdvP+RHSM3IE6fmsfTJ7NsDUmf9555DKbV3lv3Pvw4LvvNvR0eESGjw4DQ7uPHAD4A4O0Arqvb/wZU3aof7N4hIocC6HPO3VPX7r8BfAHAWSLycefcjlrb+QD+EcD9zrmVrf0UqqTlXFU6NA+Oca60kyFWn5hiyJqQTxSTSZWkT4wTU9t37c4irtxeXTwVABfMGtjTPuRUWT5VjBMW6hNF6LzW+y1BH+06aadK+1QAUNL5UhHFtr0cK2MswYLKRp9QbpX19dFjiSzC3ID2nIbLmDe5D+9+0dH48A/+DwBwwX/ejJc9+1BMqj15GsqXq+5LkEnVshyr5u/bRJlUlj8V+PZIVITZcKG8rKsEvpSF72GF+1RSujYZO127YuWc+wuAiwH8k4hcISKvF5ELAVwE4I9oDAf9PYC/qv5bAbwb1duGN4jIO0XkfQBuQPVH+7lt+DRID/DhDXtUv1dNy+MJjFcgLeQd/3gU9p1R1UwrFcdcK9Kl9G7cQrf/Bng7gNWoPqF3GoBNAL4E4LxQORugKqGLyCYA7wHwMVT/3L4ewCucc39q1aBJ73BnoYRrBocBAFkAH509qbMDIns9Uyf14ZJznoaKc/j7kw5liRvSlThgQt2+S5Ounlg558qo1gi8MNBuUZNjVwC4It2REVLlki27Rj5+0dQ+LLJKlRCSMi962kEAgAInVYR0HV09sSKkm9leruC7j++JWjtn+kCT1oQQ0kMIUOaKFWknEymTo20kCSZNUoQ5CbrgMoDvPT6EHbUxH9GfxUmTIr6d0hhvEmHfPI/aDgWGAn4Ipx6LmTynzxuQzi1yVhFmPRa9HSGVhwJDrX1aPLeiNULFw6Me7mjcth4IqZfTnXPBosy6T3SbBCJ6VB8tnus+EhMqGn7vJxHeOxUimpbwTjrHxLHBCOkyZmVlJK/qnJmTID361xnpHKVyBb/482q84H1X4oLv39rp4RAywm7HinELhJBozpw+gJft04+VkKQwAAAgAElEQVSrdw5j2aRcwsgDQpLzyz//DWeefxUA4K5VW/BvrzmBMjvpEqRjEUGdpjc/a0JSQkTw3Kl5TLfyjwhpMc9/6oGYN6v6JOr6LYO48rpVHR4RIYQrVgGcAyqVPUuQ9R+niRciGnGZiZTrkRoxLlHI1YrxaNIo3GztK0UUYQ5ep1OBoYD/t5jqU4p443qBoFaQpw5SjSmonODrEeM+hfrEvDdiwnH1viZFmXczkBW87h8OxwW1wNDv/vIuvOhZh4wcNwNCdbHkiDbeccvLSuBYeecIOFcAgj8bY4owJyGmkLMX5BkTPJqCU2UFhnbD+rmbQLfv0qQHfzMTMj5cxC8IQtrFq09ZOvLxtbeuxbadzUqoEkJaDSdWhIyRVz28Dc9+YAs+t2UXNugVKELazEH7TcMxi2cDAIZLFfz+poc6PCJCelte58SKkDFQrDj8YtsQrtk5jHdu2IlHObEiXcALnnrgyMe/+vPqzg2EkDoqkFT/TRToWKWA9q60L5WkCHOSN5E1o09jll/RlUUnGik+rXftziK21ZyURX0ZHN1f573EeDMxx4NOT4yXFTHhCzlKOrPKPG/gHAA8D8sT/a3XIJBRZY0lJsfKa5PAfYpyuVLIsfKOj543dfoJB+Lj3696Vr+94W+oDJeRy2YSZVRZxOVYjf2PDP+86hzGlz2USRV13Qhfqpth1lV3wxUrQsbAL7btSVo/fWqe2VWkK3jy4tlYMHsyAGDLtiHccOcjHR4R6XVcDxdhnjgjJaTDOOewvG5i9cKp+Q6OhpA9iAhOf8qe24FX8nYgIR2DEytCIrmzUMLq4eqtin0yghMn93V4RITsYbdnle/LYMeu4Q6PhhA6VoSQAL/cvucx9n+YlkeetwFJF/HsJ83HDz56Cp57/AGYOomTfkI6BSdWQaQxINQQ0RPJ6Vpwj5iN6zatEtyTyOoVK6FuL+OGwT2rAKdM609Y+DiBeJ6WrD7mPsbX1BLag+cMFTqOCAi1XgNdGDsbI6IHxpJERLeuEyioHCfJRxQXrmszOZvBi56xqHb96n4r/DOGUKHmOOFd9YkKoFXnsMI+A7J6kj7meQKCezuF94koqztJ5+GpiQgnVoREckvd7ZXjJ/dVY/kJIYSYlCfQ7bs02fuXGQhJgUeHy1hb86smCXB4PwvdEkII8eGKFSER3F8sIy9A0QHHTOpDrkeXuEn38+iWQdx630bcvnIz3vqaE5BlgXDSAapxC735c5ITqxYQ40+F/KiKsZiYqE+L3tghD8sZx4Mj6aYgUuW3PGNKHjuO3hd3FUrYldBbsa+T4FytKsKs/akYX8rzZqxf4oFCzUl/8QcDQlMqcB0M7ozwsoygzuB5Qp7WKPzdm36KRzYPAgBOfe5SHH7wrIbjIX8KSFhAOVHRZe1hJXCh2tQn6rwRBZa1jqp/7PWArrrXwy8hIZH0ieCYSX34uynMryLdy7FL5ox8fPu9Gzo4EtLrOEiq/yYKnFgRQshexLFL5458/H/3cGJFOgeT1wkhhEx4uGJFSGehYxXAoVF10AWXgXDR5ZicqyQZVV7+VIRPZc36dQaVbhOTURWVfaXbJHFrEvVJsIRcd51Hi2Xcs6OI+X0ZLMhlMLUbZeA0Ck2nUYQ55mHJqEynmCLMep+6eFRR7ARPd8aMP9QnxtOL8LTE86Ucjj109sj2fau3RDlVYyXJOWKcK6/PBHeuOkmlwz+mHJJlLe4NdOFvCEK6i6u3DeGkB7fisHs349VrtnV6OIQ0Zf6sycjU/ojZum0IQ8Vyh0dESG/BiRUhAdYP7/nFtCDHbxnS3WSzGew3Y9LI9qNbdnZwNKR3SbdO4ERa/Rr3bwkROVVEVojIoIg8JCL/ISJTjXavFBH+6UQmHOvr/uKf38dgUNL9zJ89eeTjR2vRC4S0G06sEiAixwH4JYD5AK4AsBLAuwDcJiJHjH94hHSedUWuWJGJxby6idX6TVyxIqSdjFde/wiABwA8xTn3GACIyNMBXAbgOhE5zTl34ziv0VkcwkWYK1pWD582KLgbk/NWCe6hIsxmn4DQXjEE72yS5LuQeG6J5CkHja4f3iPqzu/roYmVJcSnUYS5VYt+uiizFTcWLKicoAizRehhgiixPjmNK1bNJ1ZWoea2yeltEtqTsDcK7e3EASj3aPL6eH9LPBnApbsnVQDgnPsTgGUAVgP4nYicPM5rENJR6les5nPFikwAnnDwbJz4pPl4yfOWYNGC6Z0eDiE9xXhXrPYB8Jje6ZzbKCInoXqb8Bci8qpxXoeQjtEgr9OxIhOAN51xFN50xlEoDLBKAOkcE8mLSpPxTqxWAXiidcA5t0NETgXwU1RvDf5+nNcipO0MlivYVrutkxdgVpJMLEIIIT3DeCdW1wB4qYi8yzlX0gedcwUROQPVidVLgIl507rBsYqo0RoTEBrypWIKN8f10eGfEY5VwLmy9nkho0Yfb60nSWBoGhObMQR8DtZ9cadlBFL/2njn4UOviSgZ31T6tTXdp5DH1KYfN1HFnhP4U6GizMY+7UalEQZaPU+b/KkEfaLOOzF/9YyLjFYD23x9B/F+//QKwc9aRJ4lIj8RkV+IyGvU4e8C2O1UmdQmXC8D8HkA145nsIS0m3ofOtejIiYhhCShV4swN12xEpH9APwa1edrfg/gEyLyPOfcKwHAObcCwEtDF3HOOQDvHP9wCWk/i/qzGK44zOnGUjaEGKzdtBPfuPIeDIlg31mT8aYXm8YGIaQFhG4FzgDQj+otvLcC2Ajgf0Tk35xzX2jWUUQyzrn0nh8mpAPMy2ex6tj5wOBwp4dCSDQbtu7CBT/4PwDAk5bM4cSKdIReldeb/gnunLsXwAcA/CeADc65LQD+EcC7RWThaP1EZF8AV4nIghTHSgghJIa632eVFnlLhBCboLzunPuU2n5URD4L4FMAzBgF59wGEbkNwNcAnJ7GQLsFU0QPyOqJ+sTI60nCPyPCPkPb1X2BUFErDFSL5/rWmhXsqc+j+1gyu9cmRooPnNc8R3nsbbzrGuOnA9850grpTENWH8/l64T6nHpfahE9k+DGQjeJ6aQ7cejdFaukTwV+BcBfReRI59zd+qCIzAfwPAAHj2dwhBBCxs5w3ZOWemJFSLvgxGoMOOeGReQzAD4iIh8DsLj2bwmAo1FNZO8DcE9aAyWkExQrDr95rIBSoYwMgDNmDHR6SIQEKdXFLPSxWgAhbSVqYiUiOQB/D+BpAJYC2B/AAQAWAHixbl77fweqwjshE5ZdFYcX/nUTAGBKRrBjxrwOj4iQMIW6Mkx5VgsgHcBBUOaKlY2IPAXA5ahOpgCzPDCKAB5GNYn9fgA3A/i5c25rSuPsGM4BpXJ9QGhEEWavKHOEY6Ve1pIzvCbtS6lnD0ri/wD1HauI80a5W80DQa2AUM+XigkI1Q5SVJ+QLxXhZdXYJyvoF2DIATsrDjvKFUzd3TYmrNQbS0RqnxeOGSFd6bG0Kxyzm/C8pogJRVS4ZwteyxZ/fR7Zsmvk4zkzJrX0WmR8eDmwfJZ+whOzYvUdVFenAGANgFsB3Ivqbb7NAL4OYLFzbntLRkhIBxERzMtn8beh6uRm/XAFS+iskC5n/ZbBkY/nzZ7cwZGQXmYihXqmSczEaiGqq1BnOuf+Tx8UkdNRveX3iZTHRkhXsKB+YlUqY8m4K0ER0loe4cSKdAG9Kq/H/Ol9LYD/sSZVNT4G4K0iMju9YRHSPczP77mltH6Y6/Sk+1m3uX5iNaWDIyGkvYjIqSJyr4isFJH3Gcf7ReTHteM3isiiumPvr+2/V0ROqe0bEJGbROR2EblLRD4aGkPMn95vBHCZiOzvnFujDzrnHhKR7wI4H8BbIs43oXDQRZjHnklVqhheU4Icq5JyRmKKMGvvqlVFmEs5dR0jx8qpp5MkJscqlEll9gl4WDEuV12bBf17vk3WlStArta2GDiHtU8/oWUWF9b5WPq8xvjTyl9qet0uh7doR3jJiQfj4HnTsPbxAg5fNKvTw5mQpLHaYju54z+Pdd5uwwEoG35xKxGRLICLAZyMqrp0s4gsV7FQrwOw1Tm3WETOBPBpAC8TkSMBnAngKFQfzPudiCwFMATgOc65HSLSB+B/ReRXzrkbRhtHTEDoGhF5bqDZRwHcLiLHO+durn2CUwCUnHNDoWsQ0s00rlgxvZN0Py858RC85MRDUBjId3oohLSTEwCsdM49CAAi8iMAZwCon1idAeAjtY9/AuDLIiK1/T+qzVlWichKACc4565HNeUAqMZI9aE6bxyVqD/xnHPDzrlRi6U55wZRXa36gYjMEpF+AOsA3C8iM2KuQUi3Uj+xWlfkrUBCCImhAkn1XwQLUU0o2M2a2j6zjXOuBOBxALOb9RWRbK2azAYAVznnbmw2iNTWzp1zvwXwI1SjFt4MYKA2qOPSugYhnWBB/cSKK1aEENIp5ojIirp/Z6vj1uxLry6N1mbUvs65snPuGFRjp04QkaObDTLVx5ucc+eJyCMAPojqctkNtX+ETFgOHtjzbXLXrhKccxDDVSOEEFLFQVoRt7DJObesyfE12BMPBVQnQutGabOmFn4+HcCWmL7OucdE5BoApwK4c7RBpP7cuHPuKwC+IiJ555zWeycg0iivG3dWQ4GgVkCoDgDV2zFFmEtqwTFGTC8ZUnkp0yiel7LNt619fuFmYyxKLs4q4d2Uj3VqdExx5JDgHlW4eU+bxVP7cGA+i8UDOSyb0ofhTAb5jBh9rNvuCYowB2V163ZkTJsUSCKIhyT4TpZcmWiCfgRf/cXd+MaV9+DJS+bgn045HCcdt3+4Uxuwfj617dpteuzf+lnflusa3+7d8M7WIdZt4GYAS0TkYABrUZXRX6HaLAfwrwCuB/ASAFc755yILEf1Qb2LUJXXlwC4SUTmAhiuTaomoVoH+dPNBtGyQJ69Y1JFCJARwepj53GVikwIbvzrBtzx4Bbc8eAWHHnYvl0zsSKk1TjnSiJyLoDfoFp64VvOubtE5HwAK5xzywF8E8D3anL6FlQnX6i1uxxV0b0E4C3OubKIzAfwndoThxkAlzvnftlsHEw6JCQCTqrIROGW+zeNfHzsYXM7OBLS63RiBc85dyWAK9W+8+o+LgB46Sh9LwBwgdp3B4BjxzIGBr8QQshewo5dw7j34ccAAJmM4AmL53R4RIT0HlyxCuEa719HBYR6zpV/2nSKMDd3rgCgpIslG20SFWEOFF3WgaFWm6z2W7RPBfh+lHbE8kafQroBoeZxs40hOoQ8rFJEn5gizMFzJHCuYnyqGOctCVGFs1vw13DMeyPpecZ7zghuf2AzXO0tdvgBMzB5oC/YxyrMHu7TmvHHnLcVvlSSc3bKp4olo34Ptft5Zgeg3BWmV/vhxIqQSH68cRD/u30It+wo4vLDZoPmCuk2bq27DfjkJVytIp1E4Lp88tkqOLEiJJKvPLID126rFhK4Zccw9p8UXg0gpJ00TqxYvpWQTkDHipBIlk3dM5H6w+OFDo6EEB/nHK69Y/3I9nFLKK6TzuHQkeT1roArVgF0EeZS2fCYyrrocrgIs96nt60cKO1QReVYaQ+rRTlWfq6VMX7lXfVpD8tyrPQ+7VRpnyqmz7BxnaIyEHSfcgWnzJ2Ci9ZVS0Yt31rA5w6c3vi0oOUxGedpoN/4FhwqNW7HFGEO5la1yn1K4qK1yMvKxXhZLfhbMi1fyvMJx3be21dtxcMbdwIAZkzNY9lhczFsfL9b35utIMqXaoGrFfMLuJNOVdDJNQs36zapDIW0CE6sCInkWTMGMC0r2F52WFUo4e5dJRw1mbcDSXew/OY9Zc5OPf4A9OUyGLXAKyGtxgHlHnWseCuQkEj6M4JTZ00a2V6+dVcHR0NIIw9t3IndC0AveOqBnR0M6Xl6+VYgJ1aEjIEXzpk88vHyrfSsSPfwjXOfhod++Apc+s5n4pTjDwh3IIS0BN4KJGQMPH/2JGQFKDvgxh1FPFosYz8rS4uQDrDfzMl49SmHdXoYhAAA4xaIjXNAqVRXhNmQBkPyYVSoqNq2AkK9Qs26wDIMyVyyTbcBX2gPyezWPi2mWwGhep9TEq1YUm2ooLI1qdFyeiFwDus8WmavjWNWNoOnTx/AtY8V4AD8z7YCXrvfVHusgB8aGroOEC5KbPXRX6OyTqVNyXZNQ/6OKbocFeoaCuFMcp0Y+T5JIeoUXjfrfesF9aZzEyKmqLrXJyCiJxHV0xLRQ20sMT0kq5t9vHDoNj0oUPav09fuRFAyAm8FEjJGXjh3z+3An2+hZ0UIIT7p+lV0rAjZi3nh3CkjHwsEzhk1iwhpEz++bhXe/JXr8ZfVWzs9FEIIeCuQkDGzZHIfLloyC6dPy2Mx09dJh7no53fhpvs24au/vhffeM9JOOvkpZ0eEiHVpwLpWJHRaAgILUUEhKoQ0aiAUNd822yj/amIIsx62zpPTBHmknJGvO0Ix6qkgjz7YgJCQ9uA77NorylvvO11pexAsOc7lswGdqmUIM9rMtBtdFFmAMiHTxM8r3aqDE8uaryaJOGeiUJFE4R9puVUjZWoot4Rr1uCz2fF/Ztw033VMjb9fRmcekI4ZiGJPxU6h9mmQ05V1DkS/MJPLSA0iZNb0dsTY8LCHCtCCCETikt+de/Ix//8jIMxZ/pAB0dDCAG4YkVIKlSc418ppK1s2TGEy659cGT7nH9gzALpLtwEEs7ThL8LCBkHawolfHjVVhx4/RqsKpTCHQhJif+85kEUatEbxx4yC085jEWXCekGuGIVIEmOle9cGY5SoFCz5Vh5BZUDRZmBSMdK+RLFbOPbopjz3yb5UuMkIqpwcyDryivKDCRzrEJtrOyrYql5G6vActnhTXduwP9sHAQAfHXTID59yMzm/WLCRItqO8q50uOL8Y1izhs6h3pvWxlVSYowx/hGoaLLMe5W6LqAnx8V42WF/KgkhZDrvk8rFYdLrlo5sv3m044AshnPl4rKn7IKNUd4lqE+3nHjPan3VdTPJ8uX0vvScKoS5ViZ+YThNq1AjOtkOuxh9bK8zhUrQsbBm/afPvLxNx/ZgV3WBIyQlPnt7eux8pHtAIDpU/J4+bMO7vCICFE4QTnlfxMFTqwIGQf/MHcyDhqoruhtLlXwpXXbOzwisrfjnMNH/uuOke1/ec6hmDLA2A9CuoWunliJSEZE3iEi94hIQUQeFpELRWRKuLd3rskiskpEnIh8uRXjJb1HVgTvPXjP7b9PPrwNW4dZS4K0jp/fvAY33l+NWMjnMnjni47u8IgI8aneCkz330ShqydWAD4H4CIAdwN4K4D/AvA2AL8QMWSh5pwPYE66wyMEeP3+++DQWlDoY6UK/mPNtg6PiOzNLJo7Bc85eh4A4JzTjsBB+07t8IgIIfV0rbwuIkehOpm6wjn34rr9qwB8EcCZAC6LPNeTAbwdwHsAXDimgTiJCAhVwqUuqFwO9/HkdStUNBAiaorpylA2izDroNFAUWYgLKvHBIQW8423L/JaIAeQ1edJIq9rYVwXaQb80FD955HlTtWeyOrLCj5+2Gy8/LZHAABfWLsdb100Awv6c3a/sRIjs2tZWl+31CL3S8vqUWGZSQJCEwjuifokCCKNKI4MLYjHSPKjSPPHLJ6N359/Mq66bR2eeOQ8VOrOFRP+GSe0Ny/ubInqnnge8bdvGoWbPQG+RcJ7jIgd0ydUqDnmASm9nelStdNNIC8qTbp5xerlAATA59X+rwMYBHBWzElEJFvr82sAV6Q5QEJ288/zp+LYadUZz66Kw/kPsm4baS0nH7OAgaCka9n9VGCa/yYK3TyxOh7VZ8hvqt/pnCsAuK12PIZ3ADgcwLmpjo6QOjIi+OTiWSPb31i3Hfft1EtNhBBC9na6eWK1AMAm59yQcWwtgDki0jThR0QOBvBRAOc751anP0RC9vD3sybh2TOrKwhlB1xM14qkxOODRVz8q3sxxAcjyASiAkn130Shax0rAJMBWJMqACjUtWm2LHAJgFWoCvDRiMjZAM4GAOy7fzAgtP54dVvd889bAaEBxyqmCLMXEGp4TbqNWYQ5UNzZCBD0PCzljBT7jFDR4caixSXl5xSN4siT+tUvkpJ6rNz6RTOQa97G8rImqT7aUbJcqUmNYxEAnzxiLl5863p8ePEsvGb/fYChFqSxG4GzXkFl3SamQHESkoRl6j5WaKruo8NArfNEeVmBa0d5WUn8r4A/ZfbZs/3Bn9yJi397Pz7787vxhdcfjxfWii2HAkHN8E/dJkH4Z1p9Qu6THSoa6JOWC+W1aX48Fj9EVG+HvawYutW76gW6eWI1CGDfUY4N1LUxEZGzAPw9gBOdc8OjtbNwzl0K4FIAkKXHTKCHPEmnecqMATx40iLkkyRrE2Jwzd2P4uLf3g8AWL1hB3YVuWpFuh8HTKhQzzTp5luB61C93ddvHFuI6m1Cc7Wq1uciAFcCeEREFovIYgAH1ZpMr+2b0YqBk96GkyqSFjsKw3jt124c2T79+P3xz09f1LkBEUKCdPPE6mZUx3dC/U4RGQBwDIAVTfpOAjAXwGkA7q/7d03t+Fm17denOmJCDIYqDvcNjmnRlBAAwPt+eDtWbdgJAJgxJY+vvvnvIBG34QjpOE7gUv43UejmW4E/BvABVPOnrqvb/wZU3aof7N4hIocC6HPO3VPbtRPAS41zzgXwFVSjF74J4A6jTSMOqNT5KmaOlcoI8oswWzlWughz86LMAFCsqBwo7UIZrkcixyoix0oXZtbbukgzABT7VG5VThVyNrKvysqH8nKtrFIe2qnSzpVZUDlQLHmUIsyhNrds2YVX37kB20sV/OXpB2JaVEHlAEkcK33cPK8afxIvKyoHKsJR6o8onK33aUcvxt3SY7H8O32emKLeweyucI7VH+7eMHILEAC+ePZTsGBOY9EJP7cqIsdK503FtNG+VILCzXb2VaBPRCZVzHHPl2pRrpWfNxU8bSL0ebNWEWbr50SbaVcR6m6jaydWzrm/iMjFAM4VkStQva13BKrJ639EYzjo71G9zSe1vsMAfqLPKSKLah8+4JzzjhOSJrvKFTz/1vXYUHNi3nvfJnylLpKBkNHYUSjhdV+7YWT79OMW4qyTDungiAghsXTzrUCgulr1bgBHAbgY1bT1LwF4gXOOzzyQrmZSNoOLDttTRemSh7fh6q27OjgiMlF433/9pe4WYB++9oYTeAuQTCh2y+tp/psodPXEyjlXds5d6Jw7zDnX75xb6Jx7p3Nuh2q3yEXcgHXOrXbOiXOOYaGkLbxi/lScse+e2zevunsT1rUigoHsNfx0xRpcfPUDI9tffPUyzJ85qYMjIoSMha6eWBEy0RERfPXIuZjVV/1WW1cs40V/2YBCGjUEyV7JL25bP/LxC5ctxFnPXNS5wRCSFNe7JW261rHqGhyQrQsA1WGggC/oeQWWDYnQbxNRhFmHiCqBtOgMyVyUZA4/A6eYadxXdGo766+w5MvNizBrmR0A8uo8OhA0V/LHpvf5gaFWQGggRNSSjXVh5oo6hyV/xxRqBjBvUg6Xn7AQp/z5YZQdcNP2Is5euRXfOXrfZFnCxgMUfqBpEnk9QdindzxBEWMr/DOV4M4E501SRNoswqzaaOHd6lNr8603PgWHzJuG7//5b/j2OU+D1L63XERB5UQBoVYAcDYgxRu3JfXPo3QCQlvUJ1FAaPNCyNHnDRRUtgNC1bb6Xu3r0mizifQkX5pwxYqQNvDcuVNw0dF78m6/t247LvrbYx0cEelWMhnBeS86Crd95jTMmmbF+BFCuhlOrAhpE289ZCZet3Cfke333LsZv948avEA0uNM7ucNBTJxcUj3NuBEuhXIiRUhbUJEcPGRc/H0GdWKTBUAl2/gxKrXueSaB/GWy27DcIneHSF7A/yTKIA4IDe8Z6Y8bDhWfhHm5kWZAb8AcagoM2AUYQ5sA0b4pzGXDoWIFjP+26SYVR5WVECocqqUF6QDQwE/NHRY/RXfZzlWXkCo8qW0GwUYAZtqWxdpTkg/gJ8+dX+c8MfVeO1BM/Chg2cA2j0JeT9WrTjtR+n3XJQsb7hnY8VylLRDFXKjrDZRAaEJgjujwj4TFJHW+7QfVbvu7+56FG/78R0olR3u3rATP33HMzBrar85tpIRXhoKBLXCP7UPaQULh5wqO+xT9UHzbbvP2H0p/fMrkS8V00e7UJZ2GeVLhRyr8HlFbWesgNAumKdbr1EvwIkVIW1mv4Ec7nzOwZjWl7UnSaQnuO6+jTjji39GqTaJ376rhEnWJI2QCYhzQLlHk9d5K5CQDjDNWH0oO4cdvB3UE9z4wGac9rk/YbA2sV44cxL++13PxCRdlocQMuHgdzEhXcBwxeFfbn8Eawol/OpJ8zA1x7959lZufegxnPql67G9UL31PW/6AK5+z4nYf/bkDo+MkHTp1bgFTqwCiBPk6pyp4rBVhFk5SoGizNa+UFFmACiq7Khipvk24Gdb6cLN1X0q60oab0/lxci+0rlV2eZFmQEgV248T66v8XXSuVaAn2PlbRt9RGddeRlPEUWY2xjg6ZzDmTetxRXrqwUF/v62R/Cr4xdgulUQGLCdpFBBaCvHqhWfY5IcqHY6ViGnyhqL7uNlUoULKu/uc8ODW3DqF6/H47uGAQBzpvXjd+99FpYunA7oAuMxeVOqj7edtbws7UuNvaCyVZg92Me6jval1M8n27Fq7lS1LccqgZdlt1HHI7ysGH+qG4ow9yr8s5iQDiMieGbdasX1jxXwvJvWYgv9q72Ka+/bhJM/96eRSdXMKXlc9Z4TcdT+0zs8MkJaA+MWCCEd4+2LZ+GLT9gTILri8SE8+8Y1WFtgXcG9gRse3IJTv/Bn7KjViZw7rR9/eP9JOOagmR0eGSEkbTixIqRLeOuhs3DpMfNGblrcsb2I4//0EG58vNDRcRhxiaoAACAASURBVJHxc/i8qThgVrWQ8vzpA7jmAyfhSQfO6PCoCGkdDtWnAtP8N1HgxIqQLuINi2bgO0/ab0T1WT9UxrNuWY/vrd/e2YGRcTFjch7L3/JUPGHhPvjju5+JIxfy9h/Zy0n5NuBEuhVIeT2CegnQkgZDgaCW8J7va16EuVjy5dB8VknxMQGhau5cNL7kujCzFtwtKT5fGb+8rkNEi30qyBNALt/4OefUgwEZI4GuX0cW6DZmsl8gILRVGLL3qxbPxsJ9BvDSG9dgS7GMoYrDv9y1EXcUSvjUEXORtcavfayoIszq65pEZg8VZQYi5HXjHDFSeUhWt+R//bBDVJ+xS/I7Kw5T6sNs+7I47IAZuO38k5HJiH0d9TmWdUColtvhBw2HAkOtNlZAqB8iqqRyMyBUC+7hws2e8B4REOoFggbEdMD/2RjVxwsETVAs2QwR1dtjP282Qma3QkNJe+CKFSFdyHP2nYKbnn0wjqwrwnvVxkEU2jXhI4m5+aHHcNj5V+OHK9Z4xzLGRIeQvREHwFXS/TdR4MSKkC7l0Kl5XP/sRTh9vymYm8/i58cvxBTmW3U1l92yFid++QasfayA137/Nqz422OdHhIhpM3wViAhXcw+fVn89/ELsXpwGAdNbrxVWqo4ZATIGLdYSHvZsrOIf/vvu/H9W9aN7BvIZbGtMNzBURHSWSaSF5UmnFgF0EWY6z/ejS7MHFWEWTkw2qnK5/x1z6L2sAKBoYDvS+WzRtincm1CgaHVa6mgTl2UOWs4VrnGNoV8vvG44fjoQNDMgHaswn2yuk03BYRat4aUa5PJCg6pn1TVfKqP/3Ujrtk0iG89eT4OmaT8tEQBoS0qwmz5UfUkKXxs9fPCP40fb8E+EaGianv5PRvxxh/dgUe2DY3sO3zeVCx/y99hyX5Tqzv052P4UnpfKPwTMHwpXdzdCAjVTpUVPBoKBLUCQv1wT13c3egT8KWsovEhD8tyTUPhnlbR+3AR5oQBoUmKMKvfF33qR7IVBprpdAyesz//XoATK0ImILc9VsAF921GyQFPuHoVPn3EXJyzaAZXr9rIlsEi/u3n9+D7t65r2H/W8Qvx5Vcei+mT/YcxCCF7P5xYETIBuWHrLuxegxosO7z1zg34yfrt+NaT5uGQKfmmfcn4WX7XBrzxp3fhke17Vqnm7dOPr535RLzwCfOAfv5oJb2NA1Du0VuBNGEJmYC86eCZuOHEg3DUtD2TqD9u3oUn/HE1LnpgC4baWOuw1xgqVfCO5X9tmFSddfxC3PWBk6qTKkJIT8M/qwJ4jlXJn4HrnCqvCLPZR7lPuiCx4WXltT+hvIBixXCsdKFm57fJ69wqz7kKnzfnlJdluFzascopmaCQ92+daIcqo5651T5VtU+jTzRZTTLE8o2swJl2EFOAeJTtZZP6cMt+U3H+nRvw6b9uRNlVV6/edfdGfGn1Yzj/yLl4xYHTkRUZxStrV1ZXi4owa/cpps9AIMdqwLh9p67TP5DDt17xJJz05Ruqq1SvPAYvfNL85mMLbSOcW2VlUoWKLpvZV7qNVdw5kFtlZVLp8+iiy2YmlfrZ4jtXxtgC+X1mEeaAU2X6Up5TpcZRNsamfm7bbZo7uZaXJF5u1cTIsXJdMIZOwBUrQiYw/dkMLnjSPNzw3ENw1D57Mq9WDw7jX1asw7G/exCbhlhvMCmDxTIuueEhlNVv1Wctno3vnfUk3PXeE/1JFSEEDtXJaJr/JgpcsSJkL2DZrEm45eRDcMkDW/Hxuzdic+3JwRn5LGbnsxPrp1IXMFyu4Fu3rsNH//Ag1m8vYlIui1cvW9jQ5qxl+3dodISQboYrVoTsJfRnM3j70tl48LQl+NDhczAlK/jU0ftC1G2bhwaH4RwnWha7hsv49q3rcNSXbsCblt+D9duLAIDzrlqJwnCnn18nZALh0i3APJGKMHPFipC9jH36sjj/qH3xziWzMUM5REPlCo67ZjX268/inEUzcNbCfbCPVbeux1i5ZRcuuXUdvn3HI9i6q/HW6fxp/fjgcw5BluVoCCERcGIVQJwEA0K10O4JjLooMHxhUQvwWlQH/BDRopbBrYBQLa9bIroSVXOBwFAAKChbMqek8kLWCO6sqPGqEFRdpBkAcioJT4eIZgb8lRctr2cqjfEDA14PGLqrIq1baZ7IbfhPIXndmgjp1ZRyBTMm1X3daqL6T1ZvxaZiGZuKZbzlLxvw3r9uwqsOmo43HzoLT5hhvTIpEJLVrVBRPYmxwj71ecZYULlccfifB7bgKzevxW8e2OI1nT6Qw/tOOhhve86hmFx/7ggR3dun4hec8f0dCgQtGq9BKBA0SYFlwBDRA9uAESoqzWV2IBwIGhUQqsV0qxi9Fty9wFCviye4axHdDPLU9d6tNoGAUOthJy2n6+3csP85dzog1IEBoYSQHmBDoYwpuQx21ib7O0oVXPLAVlzywFY8c85knH3oTLxg/jRvpWtvpFRxeMMv7sGGnY1lZw6ZNQlvfuoBeN2y/TFzcp/9dCEhJIjr0RwrTqwI6SHecfgcvPaQmfje6q245P4tuLuuDMt1mwZx3aZB5AR41twp+H+Hz8Ep86Z2cLTpUK443Lh+O0oVhxMPnTWyvz+XweufvACfuO5vEACnLZ2Nc55+EE5ZMgcZ3vYjhCSEEytCeozp+SzOXToHbzlkJq7dOIivPLAFV6zZhlLtNkbJAb/fsBNvOGRmZwc6DnYOl3HV6sewfOVm/PLBrdg4OIxn7r8Prq2bWAHAG49biIoD3njcAiyaOcnOsSKEjB1n3wrtBTixCuEVYTaCO4cab6wPKu+haPgUuVxjn3yf6mMFhPYFCjcbXpMODbUKNRcqjW+DnLqBX4ThPulQUbWdM7yNgmv8pZVRT6aZjpXal8lrf8r/nLV35TtXvlDRn8Sh0qsaMQWIC8qpsvroCqt5tV005Al9u0q3MQJCBTk8a2o/nnXwTKwfHMa379uMnz30GFZs3oWcAKcePLPhvMMVh6N+fg8OnZbHstmTcdzsSThu9mTsP7nPe/KwgdDrFPO6WW1qYytVHO7ePIhbNg3ilkd24JZHtuP/Ht2BIRWA+qe127Cp7DCnruTPgdMG8MkXHL6nkQ4QBXxfSk++rPI1/aqN9qcML0s7VJ4/FRP2GQgMBQxfymzTvOiyWVBZh31624b7lKAIczAg1Lj95PlSgW0gXHQ5KlQ0qgizOm4UVO5Xbq8uumz5VN0QENqrcGJFCMH8yX34wBP3wweeuB/W7izils27MF1N1v730R24f9sQ7t82hF+v3T6yf+5ArjrJmjUZT549CYum5rH/5D7sOynd1Z+hcvXXbr5ukrW1UMLCr92MXcYDIvXsO6UPpy+dg8HhHv0TmpA20yl5XUROBfAFAFkA33DOfUod7wfwXQDHAdgM4GXOudW1Y+8H8DoAZQBvc879RkQOqLWfB6AC4FLn3BeajYETK0JIAwun5LHQKOT823XbjdbAxkIJv167vWGy9bS5U/Cn5y9paPfHR7bjmvU7MHOgD7kM0JcR5DIy8jRZqeJQqjgUKw4bCiWsGyxi/c5hrBscxvrBYWwplHDF85fiRXW382YO5DBrIIe1O4reuI7edwpeuHQOTl86Gycs3AcZEfsJPkLIXoGIZAFcDOBkAGsA3Cwiy51zd9c1ex2Arc65xSJyJoBPA3iZiBwJ4EwARwFYAOB3IrIUQAnAu5xzt4rINAC3iMhV6pwNcGJFCInigifPx2sWz8KKzbtwy+ZB3LJ5ELdu2YXtxirQgsn+j5bfr9uBj93+yLjGsH6nP4E6br+pyMgOHDdvGo6bN3Xk/31nThrXtQgh48ChE6GeJwBY6Zx7EABE5EcAzgBQPwk6A8BHah//BMCXpeoynAHgR865IQCrRGQlgBOcc9cDWA8AzrntIvJXAAvVORvgxIoQEkVGBEunD2Dp9AG8oia2V5zDym1DuGXzLqzYPIg7Hytg3eAwDtvHz8RaNzjs7aueNy4mLCvANsMvu/wFh6E/l+FqFCFdhIN04lbgQgAP122vAfCU0do450oi8jiA2bX9N6i+DXWsRGQRgGMB3NhsEJxYBRBPXjcCQtW+UGBodZ8SxJUUn7MCQpWYmss4ddyQv3VAqCGv50SNxTW2KYjvyuSgAkEDgaHVsag2SrYv9PnX0YJ7aBswZPWBsXs1/eEmydASdsGYbHiytzQ/DvgzEy2zl42ZixbarTYBMgCWTspj6X7T8PLRGtXG/49LZmG/aXlsLzmUnMNwufp/NiOoODdyazCXEcyZnMeCqXnMn9JX/X9aP+ZM6mtMP699jiNfKz2xsoR3LafHhH1qWT1BQGhZbRfz/ns9FBBqSebFPiW8ZyPEdC28G2GfSQJCtZzuBYbCGIsXCNrYxhTRA/K6KaIHxHMd/lk9T/M25s91r40x/lJzwV2MyYgW0f3AUKPP3lmBaY6IrKjbvtQ5d2ndtjWT0z/YRmvTtK+ITAXwUwBvd85tazZITqwIIW3hBYtm4gWLZoaT2AG/DXOlCJlwGH9fj5dNzrllTY6vAXBA3fb+ANaN0maNiOQATAewpVlfEelDdVL1A+fcFaFBsggzIYQQQvYGbgawREQOFpE8qjL6ctVmOYB/rX38EgBXu2pV+uUAzhSRfhE5GMASADfV/KtvAvirc+6imEFwxYoQQggh6eKAcptL2tScqXMB/AbVuIVvOefuEpHzAaxwzi1HdZL0vZqcvgXVyRdq7S5HVUovAXiLc64sIs8A8CoAfxGR22qX+oBz7srRxsGJVYA4x0r5UcONt3R3GX1ChZqt+/c6NDSnbpfowFDA8LAM58JzrHTYp/jruQX11tFttHMFABmEfCmjj26j0vSssM9Mf/rrz1HOlXW7KuRUWX10QKgusKyPA75jFREQGjTGrT4hLK9JkyRY1QxSDbhO1i3GUNin5UvpAFCvwLLvS5VVHz/8MxwQGtNH+1LFXK7pttVHh39a+0Lhn9Y+/XPEDPsMFF2OKcIccq4AKxA0SdinOm58+3htIgJCtXNl/OgMFl22fCrrd1U76VSOVW3Cc6Xad17dxwUALx2l7wUALlD7/he2fzUqvBVICCGEEJISXLEihBBCSOr0aq1ArlgRQgghhKQEV6wCSAXIFffcXs0XDF9KFUfW97b7ioYv1acKEKvcKl2k2Wyj/QqjXpqXFVXyv+SeHyXN/anqvsY2BZVrlTFEgZxrFAH0+HPOH1uSHKsQlsuVhLxyhSRJbID2p6x9MY6V10bnWEU4VkmcqhhC3lWMm2a5T/r1jsmx8vyogD9l9ml0qrRPBVi+VF/T44DvUBVVrpvOrLL2xWRUFbO5ptsAUMzkmm8bjpXvVMU4VqpNjC8VyK2ycqy8Njpvyih87LeJyLHy3NmIHCt1bV1wGQgXXe7KHCsHuB4tBM2JFSGEEEJSpVPyejfAW4GEEEIIISnBFStCCCGEpItrnV3Q7XDFihBCCCEkJbhiFaAaELpnO6YIc35ICYx9vmCtQ0OLxRh5vVH0zPepYE9DpvQE8Yx/Xh0aqoslFyqGVK7OkxMlphtvLR0a6gWGZo2wz5C8nuB5XitUtCJj/xujoqTrvCFLZ7U8nIa8HtMnRkzXRZfb9Wy0fk2S1goMCe4xIrreNkI4teDuF1S25PXmsrrZJ1BQOSbsMyogNKagsm4jetsKCB172GcwIDRGXo8owuwVVA7I7NV9ughz820gXGDZ2pf1AkLDv2N0m5xVy73DfpOD0LEihBBCCCHjgytWhBBCCEkXBzhjVa8X4MSKEEIIIani0LvyOidWIZw0BITm/Hqr3v1t37myQkWbe1i6KHN1n/KydPFnq8CyCgjVhZsBPzQ0I9qfssI+daho4wuj/SnrvNrTssI+Le8q1CcJ2pdK0qdiuCo6RLRPXycm7LNVjlWSgNBQ4WaL0GtrBXnGFGrWfpT32qbjWA0rx8oL8sz7PxRCTlVM2GdoGwg7VVaBZS8gNBMOCPUdq4iwT8+fMrysgC9VrBjjL6vrxASEBgJBbceqeSBokgLLAFDSAdKeL2WEfQba6ADR6j5vF2kTnFgRQgghJHUorxNCCCGEkHHBFStCCCGEpItrX4JLt8GJVYBqjlWdYzXsOyb5gvallNc07C8M6vMUVPZVMWf0UdlWOtcqZ/hIXo6V1UblS+WUl6GLMgNGBpX2p8TwpUJ9jIwty9VqOJ6SYxXCcrC0U2W3Uc6FziIrlrw+2Vyg6DIdK7tNTI6Vev11JpX++lj7tC9l9Qk5VVGOledLGQWVA05VTIFlM8dKOVR+EWbjvMqhKqpfL/o4ABSd6qOcqmRFmA2vKVERZuVQugh/qqTb+OOXkC9lfBv6ThWablt9OoH+XMdLe37ajx/eCiSEEEIISQmuWBFCCCEkXRyQTTnHyl/f7064YkUIIYQQkhJcsSKEEEJIqghsX6wX4MQqgDiogFCjQKbKB8wXlMBo9NGhoVpyLvbFyOtq2xB8M0pmNaVyTzxvftzuEyGi6zZaZje+C0PyequICf/U+0rGAwchWd0Sn3WbXKnRTPXkdgDoV6+TJ7NPcHndeP294s3qtXTGdUsBWT1GRA/J7IBRHDmFsM9CPu/30WGf+rpGaLDXxhLRRcvqWkw3JP+2FVQOyOtRBZUj+ijxPEZM94swe02QKzWX162HnbSc7vfpwoBQJx0vBN0peCuQEEIIISQluGJFCCGEkNSRTq+adQiuWBFCCCGEpARXrAJUHas927m84VipsE+vCHPB7+P5OboI85A/5y0GHCvtU1X3aQ9r7IWadWBodZ/qI823gYiAUMunatPU3wv7FGm6Hd9Ge1jhgFDdRjtWetval1GFg7NWqGjIl2pXbLLlTyk/yvKlKtnmr60V2BpyqmICQrVLV+wzijArh6qkxqr9KWuft22FfXpOVbjAckn7UlYRZs+pUp+P8Y3pBYLq8E839oBQqwhzKBA0WRHmcECo71hZAaGN13FmOHTzQFDLlwoHhFqOVWf9JnFAlo4VIYQQQggZD1yxIoQQQkjqMG6BEEIIISQFxHX+dmSn4MQqgFeEuei30dlWOtfKzL7yPCzti/j+y5DyrnzHyiiwnKBQc6Icq4A/VWsUOEf4zxvvvAluZsf4Ut7xmBwrw1/znKphlWNl9NFtYhyrjPKlPOcqbzhvqk9GOVX6eKtIXuC6eZs4X6q5PwX4X6NQRlVMmyjHKpBRZbUpZBt/+Fj+VCGj2ohxXmnuVGmfqrqvuVOVVo5VUflSxZK6rpUvFcitsjOpmrcxHSv1c73PaBPKrYopqOxvh/uQ9sGJFSGEEEJSRyivE0IIIYSQ8dDVEysRyYjIO0TkHhEpiMjDInKhiEyJ6LtURM4XkRtEZKOIbBeR20TkgzH9CSGEEJIMcUC2nO6/iUK33wr8HIC3AfgZgAsBHFHbPlZEnuecayblvBbAWwAsB/ADAMMAng3g4wD+WUSe6pzb1crBE0IIIb1J79YK7NqJlYgcBeCtAK5wzr24bv8qAF8EcCaAy5qc4icAPumce7xu31dF5H4AHwTwOgBfDg7ENcqEprye1yJ64/GKLhZr7NOFmvNGQKjOyiv2xcjrWioPh4hqQdwKFQ0WYbbkdX0OHQgas34a832agtCut0uWvB4TEKqFahUUqUV1ICIgtBwREFrJNT1ebdP8a6Rl9pg+MVgiut8mRl5vLqtbDxxoOT0kplvn1V9Dq08o7NPsEyyobIjoAVndFtMjijBrET2wDYRldSvss1jOjmkbCBddtgoqa6E9KiA0EAhqCe/aKTKLIwcCQWMKKuun7eyAUG8XaRPdfCvw5aj+Gv282v91AIMAzmrW2Tm3Qk2qdvPj2v9Hj3uEhBBCCPGpLUqk+W+i0M0Tq+MBVADcVL/TOVcAcFvteBL2r/3/aPKhEUIIIYT4dO2tQAALAGxyzg0Zx9YCeJqI5J1zxs05GxHJAjgPQAnNbyMSQgghJCGC3o1b6OaJ1WQA1qQKAAp1baInVqjeVnwqgA845+4drZGInA3gbAAYmHwgcsU9bw5DE/C8K93GLKpZbO5UWdepqHDPQV0MNsaxMtuogDrlYFh9fKeq+XELs+iybmNcuwHj+9Zzn9B82+yTpMCy5eeofflSOCBUO1S5sipwbRVh9vo0d66q+5oHhFq0IjQ0iXNl9dOvpVmEWTtV2ssyvs5+QWV1DmNsIafKLqgccKysYskBp0qHgVbbqD6WL+U5Vc0LLANGAeVAgWVrX6jAMhAOBLUcK+1DFYezTY9X9wUcK+Pnug4EtQsqNw8EjSnCrANBY37HtB03sZ7kS5NuvhU4CKB/lGMDdW2iEJGPATgXwKXOuU82a+ucu9Q5t8w5tyzfPyf2EoQQQgjpcbp5xWodgCNFpN+4HbgQ1duEUatVIvIRAP8O4NsA3pTqKAkhhBDSgKB3izB384rVzaiO74T6nSIyAOAYACtiTiIiHwbwYQDfBfB651x7CqARQgghpOfo5hWrHwP4AIC3A7iubv8bUHWrfrB7h4gcCqDPOXdP/QlE5DwAHwHwPQCvCQSKmohrdKgs9ymvXA7tQllFmPM6pyfbOLSYorO6UPOQ4Xp4GVWWL5UJFGE2pt8hhyrKsYpo43caexfdJ6YIs+9pGQ5GRvlSVm6S/pqp7ZzhNWkPy/OnssbXWedWeY6V8XV3zR2rGJ8qM/ZvKc9NM9vo76mYHLGIPqEMKtOTCxTbjsmk0ucwHatMoI/lWAWcqpiMqoJZUFkVd3YRjlXAqYrLsWpeYLm6r7lTVRw2vu46t0q5UMVihJel2mSjCiyPvU1UEeZiuE/H4wmcna/VC3TtxMo59xcRuRjAuSJyBYArsSd5/Y9ofKrv9wAOQp3KLCJvAfBRAA8B+B2AV0jjD+JHnXNXtfSTIIQQQkhP0bUTqxpvB7Aa1Sf0TgOwCcCXAJwXsfq0O+fqQADfMY7/EQAnVoQQQkgLkB51rLp6YuWcK6NaI/DCQLtFxr5XA3h1K8ZFCCGEkNGpFmHuzVuB3SyvE0IIIYRMKLp6xaorqEgwINST1XX4Z0yfBNJjTKhoMavldb+NL7g3Pw4AGUOKDfUJkuDdaEr+6u8FLyBU/MLHnvis+pTEN0G9MENTXtdhkqrAcsU4b0Bwt4owaxE9Tl5vLqvrc1rEhIqGsCRzr01EgWs/sNUKFW3+NYsJefUfSAgL7174pyGie310kGeEvB4T/qlldS2qA4aIniDs09uOKKgcCv8EDFk9UGAZ8IX2UPintc+pc/QZP6PzQ+OX16OKMOtiz0YYaMcDQtEFAn2H4IoVIYQQQkhKcMWKEEIIIakizl9Z6xU4sSKEEEJI6hgGRU/AiVWAmIDQSlbd7+7Tx/0+OjRUe1imO5Rt7lRpbwsACl7IoD8WLyBUO1ZifAIBIurqGn1SCsXXgaBJijDrAr/OcKG002O8TiVpHiJqBYSGPCyrj3addBvLhfICQgOBoVabNLD8Ka9NTECo9zUMB7aGnKvqvrE7VknCPvX7yfOljPdXyKmyCyyHwz69QNBAgWWzjRf+aQWENnesrILKvi/VvMCy2abY3LkCgKJyWL0Cy2ZAqPr+Nh0rvR3hZRW1U6WPe116tpxMN8CJFSGEEELSxQnjFgghhBBCyPjgihUhhBBCUkVc78YtcGIVoOpYjS3HyvOlrMLNysvSbXSR5mqb5oWaTQ9F7SoYa5TFiELNIfwsrLF7WRYVp1+nCF/Ky61qvm318f0p/yeENzZjLN55ohyrQI6VkX3l5VipNpYbFSzCHNGnVcR8zTQ6b8outt3862H1CTpWlvukHSvtT0U5VtqXChdULqmbEAUo4RPhjCog7FQVysZYAk6V9qmAcG6VWYRZOVTaubIKKvtOVbiPLrLs5QgaP0yT+VIRmVReoebm29XzeLvaTq96XrwVSAghhBCSElyxIoQQQki6OEAorxNCCCGEkPHAFStCCCGEpIoAyFJeJxb6yQZLCPSEauVb5vJWn8ZtX3j3JeF8ISByW0Vn1Xms4NEhJeOGijLbbdKR1dPAk8o9yT8cEFpC40+EClSqn9XHeKFyaLQ3ddHZnCSQ142wUl9e1yJ6goBQU15vj42qJXO7zdgfUtAiuvc1tEJeA7K6ls6tfbqPFfbptYkI+ywF2phiekRB5UIlEBBqhn0qwd0T040+AVldi+mAX2Q5VGAZ8OV0va0LLAN+keU0CizbbZofB/yn67yAUP/Hkym0k/bAiRUhhBBC0oVxC4QQQggh6SDo3VUzyuuEEEIIISnBFasQXhFmK4hNeUzan9plOD26ULDnXIUDQnVAXcVIY7NCQ0MUxtyjfXj+lDNe20CIaEyQp+dlWX2gCiwbro3nWGl/CoZjFfCwMvDdJ+1d+Y6V4UuheZtOBoRq7LDPgGNl/N0Y8uKsPr4vFVN8u3kfK+xTh3v64Z9WQeXmTpXlT0W1SVBQOeRUaZ8KiAj7NNynUCCoFfapvauS8pj6h/z3V14VYU6jwLK1TweCWrfPvBBRLzDU6uPvaysOMBTSnoArVoQQQgghKcEVK0IIIYSkSi/HLXDFihBCCCHp4qryepr/YhCRU0XkXhFZKSLvM473i8iPa8dvFJFFdcfeX9t/r4icUrf/WyKyQUTujBkDV6wCiANyQ3u27SLMjdv6friZSeWdJ6KPkb/U0MfKZzJcLb+f7tO43UnnqqJehkpf+E+gUHFk08sKOFVRWUXwx5ZT+3Qfy7EqqmLCOadzrPw+nmMV4WV5TlXAuRptXzuIcqzU34kxxbajcqy0fxfz3gjkS2mfKqqPM3yjQEHlKH/KatOKgsrDxnW8NhGOlVdQWZoeB4Ci8qX6VJuckX0VyqRKUmDZbBNwrqx9fq5VOPuqFxCRLICLAZwMYA2Am0VkuXPu7rpmrwOw1Tm3WETOBPBpAC8TkSMBnAngKAALAPxORJY658oA/hPAjY18mAAAFxxJREFUlwF8N2YcXLEihBBCSLrUcqzS/BfBCQBWOucedM4VAfwIwBmqzRkAvlP7+CcAnisiUtv/I+fckHNuFYCVtfPBOXctgC2xnzonVoQQQgjZG1gI4OG67TW1fWYb51wJwOMAZkf2jYK3AgkhhBCSKtWA0NRPO0dEVtRtX+qcu1RdVqP9hdHaxPSNghMrQgghhKSLa0ny+ibn3LImx9cAOKBue38A60Zps0ZEcgCmo3qbL6ZvFJxYBRD15sgN+RNYPzRUB4ZaUrnebuxjhYqGzpE3wkB1qKh199cS8kMUDbl+vFQqxuuUDwj7loje17xNyQxfVW2UKGwVWA6FQAJAThfsVTK7ltur+5S8rkR0q3Czlte1rG4J77qNLrBsCe/eOdoks8eJ6M0lc+s8MSK6luK1eG6L6LpNWF7X4rl3Dkte1yK6lteNIE99npiCylpWjyqorMM/jYDQwlAgINQS0QOBoFYR5mwpgYieRF4PFFi22oQKLFvn8SV5o48htPcANwNYIiIHA1iLqoz+CtVmOYB/BXA9gJcAuNo550RkOYDLROQiVOX1JQBuSjIITqwIIYQQkiotuhXYFOdcSUTOBfAbAFkA33LO3SUi5wNY4ZxbDuCbAL4nIitRXak6s9b3LhG5HMDdAEoA3lJ7IhAi8kMAJ6F6K3INgA8757452jg4sSKEEELIXoFz7koAV6p959V9XADw0lH6XgDgAmP/y8cyBk6sCCGEEJIurjeztABOrMJ4RZj9Jv797XDYp+9Y6T7+dfxCzX4bv0/MffbmqRvW+L3irxFXCWE5Vl4brwhzKdymr9T0uLVPe1iW36K9K9uxqjRtox0swA8N1R6WFSqqz+M5VoaXFSzCHONYJXtoZsyYvpT5EE/dccPL0q5TVIHugFNlFm4OeVkRYZ+eC2X8UAi10QW9gaRhnxEFlXWbFhVUDrWRIb+PLqicV0WXdUF7s09B9xl7geXqPr0ddrdCHpYliffqpKYb4MSKEEIIIenCFStCCCGEkHQQxNf329tg8johhBBCSEpwxYoQQggh6cJbgWQ0RL057NA1HQjaeNTqEwoItdP1m2MJ75ZgGTxPRvexLPnG75g0ZPYoeV21qfhONir5smqjxPSc/91eyqk2Svq1QkU9ed1oE5bXLRFdB4KqwFBDXvdkdTQPDDX7GGMJ9elqeT2iT4yIHupjPtgQaGP2CbRJEvZpy+vhsM+S1yYiIDQgq5vyekBEj5HXtayupfPqvuayuh32qbcjAkIjRHQ/3HPsYZ9RoaK9GRDaFXBiRQghhJBU0YsSvQQnVoQQQghJnV6dWFFeJ4QQQghJCa5YBRBn37+uR4dwWoWa/T7eldR2RKiocqHyu6wr6blz2KOJCxXVn0Bz5woACoFHb01fSu3zA0KNPjrss1xRxyP6qAtbrkpMiKjnWGWTOFbhIsy+L9U8MNTaF+NP6VDR0DnTwvalAsG2EY6V508ZX8NQGzNUNOBLmY6V9vpiHKuKbpNpehzw/aiS8X0ZDvu0vCxRbcYf9mn1CTlVti/VvI12sADDwwq4UdU+MQGhzQs1JwkVtVaGOr1aVL0V2JueF1esCCGEEEJSgitWhBBCCEmdTq+adQpOrAghhBCSLnwqkDSj/s1hvVFCRZhtX6p5G8tz8r2sJNlX1t1f7c2M3cvSzlWmHO5T0FcxcqxC2Va2l6WzrrSXZTgx5cYvbD6nvCwr+0rnYxk+Sy7TeB5dJDdvvKF0H+1UZcR/P4U8LNOXkgSOlXHtUJ8khDKqzD7av4txrCL6eE6V9vFiMqkC/pTZJrBd3ae8Js+fshyrzNjb6ByrBAWVSyXD5UqhoLLvSxl9dAHlobG7UJ4/FZVR5TUJOlX275jm543xskj74MSKEEIIIanSyzlWlNcJIYQQQlKCK1aEEEIISR2uWBFCCCGEkHHBFasQDsiU9mwmecGs4sgh4d0vyuy38a8z9j5VQrK6P//OlPW1wn38czRuFyPkdX/b/5Mo1McKRCz1NQ8V1ZIwAOQyja+BVdw5l208T15t60LO1fNWmm5bgrjXRnRgaDhUNCSzW4Rk9tjzaFomrwfaWA82eAGhESK6HzgbI6LrNvo9Ge4TI6Z7BZYNEd0P+1SBoUZwp5bTY8I+dZ8kBZV1+KdVeD4kq+sw0Op59XZMseTmfaxrx4jo+melDt40A0JL/r520ssBoZxYEUIIISR1eCuQEEIIIYSMC65YEUIIISRdejhugROrALoIs3XbOqMcKv1msu6ZW6Gh9diFkHWIqD6epI/VL0lAaGOfjFHp2PMEKnrb71NUbQq6KLMxNO1tlPLKp8gbHo32WTxXxfA2+sbvYWkHy+oT8qdi2piOlXaqAtsWSfypmPNarlOwT4wv5RXxHruXpb0mM3A20MZ6r2j3yXOuLDdQ+1IRBZZjwj49x8oL+wwXVPacK8OX6lN9Yoojh5wqs0/AqdI+ldlGuVD5XdbYwh5W6PeD7W4FvKyhcB/SPjixIoQQQkiq9HJAKCdWhBBCCEmdXp1YUV4nhBBCCEkJrlgF0MuZZoFM4/52xJnVdhIXSpM0xypUzDmmcHMMoXm8dc7m2T4FwyEJZ1+FvZN8XvlTpt+i8qYMXyrkYWmfyjqP50tFeFkxvlSMu6VJ4mF559Dv9QSZVRYhfyqmTZI+li/ltdH+lJEvpRXDmD6hvKkYLyumOHIoo6rapnFfWTlVfcZ1/EyqmHypsRdUDjlV5nUC+VJ2RlXzPlabGCc31MfMserwalEv3wrkihUhhBBCSEpwxYoQQggh6aKqlvQSnFgRQgghJHV6taQNbwUSQgghhKQEV6xCKAHPCm8r5Ru3k4iFemafG4qRgtMQ05P2GXuIqP8561BRq/ho43mLA0q4NoReXczZD/805FYtq6s2VqioloDzff4XWrfR4nnOKJztt1Gfc8l/siEkvGcMST4cEOp1SUVebxUxIrrfp/k5AP9rGHMdLY2HZHZrn3cOM1Q0IK8b73UvQNcI+wzJ6tZ5nQ7u1FK59X3nSeVaTA8XR/aDPK3gzlAfr4sRKjr+IM+YNjG/L/QDUzGFm9sN5XVCCCGEEDJuuGJFCCGEkNThihUhhBBCCBkXXLGKIDTrDt3/TsbYCzdbffyg0bGHiMb91aF8qbJRUHlA71NekKFphTysUp9RXFiJM14h51w4VDSXU8WTE3hZ1nmiHKtc8zZm2Kd2rAKBodV9ajuqT9jd6hSWg+S1CQWEmmGfzc9hhn3q92nAuQLChb9NL8vzpcbuT9kelhqLcp908WTAKEAc2AZ8h8oPCI0pdJykjzpuhoqOvU8abUwvS/1O0U6uGVxtnKed9LJjxYkVIYQQQtKlhydWvBVICCGEEJISXLEihBBCSOr06ooVJ1YBYu4T63vZXq6VEeuf7IVvXrg5ro9FKLeqNV5WptK4YFrUMgt8V8vzp/qta+s+jcdLOf86nnfSrzwnw8vSHop2owDfw9LnsRylXLb5ea18Kb+NdqGMPkYx52bnqJ4nhSLMEW9J460Q0SfCsdIZZ1GOlW4Tvq5+b3gFlg33yXO3IlyoUEZbUsdKlPvUrx0lyyf0nCqVaxXjPsV4WZ77NPbiyPld4y+o3LIcK+P3hZdz6I3N6uPvI+2BEytCCCGEpArldUIIIYSQtOjhIsxdLa+LSEZE3iEi94hIQUQeFpELRWRKO/oTQgghhIyFbl+x+hyAtwH4GYALARxR2z5WRJ7nnAsVqBtvf0IIIYSMEQFvBXYdInIUgLcCuMI59+K6/asAfBHAmQAua1X/ehrFwbBV6wWGWoWb1XYnZXb/zR8uwuwHgIYD60LXtfoUB0J9DDm3T8vr4VDRUkmHijb2GbaE977xC+46QNTalxnW5zCkctVGS+ZWEGkm07yPRYwUP9ZzxBAjpnt9jD+bQucx5XVPeG88HiOi67GYfQJtksnr4WLJcSJ6hOwdKKCsCyEDYfHclspDxZEjRPRAYKg5thYFhOrbZVF9In529uqkphvo5luBL0f1t/Xn1f6vAxgEcFaL+xNCCCEkCTV5Pc1/E4WuXbECcDyqNU9uqt/pnCuIyG21463sTwghhJCETKTJUJp084rVAgCbnHNW5b21AOaISN44llZ/QgghhJAx0c0rVpMBjFbOuFDXZrRSk4n7i8jZAM6ubQ59FHLnyEFrBq73dbj45QRiDoBNnR5EUrQpZH3ZO/hWmNCvbZezV722hvpk7msTe9Vr22Uc1s6Lrcctv/kIZE7Kp50Q741unlgNAth3lGMDdW1S7++cuxTApQAgIiucc8uaD5Ukga9t6+Br2zr42rYOvratQ0RWtPN6zrlT23m9bqKbbwWuQ/V2nVG0BAtRvc3XbEFgvP0JIYQQQsZEN0+sbkZ1fCfU7xSRAQDHAAjNvsfbnxBCCCFkTHTzxOrHqGosb1f734CqG/WD3TtE5FAROTxp/wCXxg6YjBm+tq2Dr23r4GvbOvjatg6+tm1CnEtQRr5NiMiXAJyLanL6ldiTnP4nAM/ZnZwuIqsBHORcYypfbH9CCCGEkDTo9olVFtUVp7MBLEL1iYAfAzjPObejrt1q2BOrqP6EEEIIIWnw/9u711g7qjKM4/9HKNciiIC0BVpoLSVCYlAU8BaI4YOIIIZEsCgoEqQQrA1GEbm0BNBAIRY1WlCRAELVCigQMLHRiA3lGhQp10OBlgJysZSWS/L6Ya0D2+OutGfW7H3mzPNLJnPOXHbeec+cvd+9ZmatEV1YmZmZmTXJSL7HqhaS3iFppqQHJK2V9ISkCyVt2Yv9R7MquZE0VdJsSYslPStplaR7JH3HuS173knaQtJjkkLSJXXE2yQlcitpW0kXSHo4v8azkv4k6WN1xj7SFXi/HSvpNEn35feE5yTdJukYSRs+gOQoIunbkhZIejT/Lw8M83W+KOluSWskrZR0qaTtC4fbKq0rrICLgLnA/aRBmheQ7ru6QdL65KPq/qNZldx8GZgJPALMBk4FlgLnALdJ2ryuoBui5Hk3m9QRoyWVcitpInAn8CXg18CJwLnAAKlrlzYbdm7z+puAOaSnvGeR3g82An4OnF9f2I1wLnAg6T3zheG8gKSZwOXAS8ApwE+AzwOL/IW2gohozQS8jzR+4G+GLD+Z9AThUXXuP5qnArn9ILB1l+Xn5P1P6vcxNjW3Q/bZG3gD+Ebe95J+H1/Tcwv8BXgCGNfv4xlJU4H3hP3ydhcNWb4J8CjwYr+Psc/53a3j578DAxu4/3bAatJ4uht1LD8k5/20fh9jU6e2tbAcCQi4eMjy+aRe2KfXvP9oVik3EXFHRLzUZdU1eb5n5Qibq8h5lx/mmA/cDPy2ZIANVim3kj4OfBT4fkSskDRG0ha1RNo8Vc/bd+b58s6FkTp2fo5UFLRWRDxa8SUOI3U9NC8i3hyYLSJuIBWubf48q6RthdU+pG9Qt3cujIi1wD15fZ37j2Z15WanPF85/NAar1RuZwLTSF2QWFI1t5/K82WSbgDWAKslPSip7R9MVXN7O/Ai8E1JR0jaRdLuks4DPgCcVT7kVhnM/9+6rFsMTJM0tofxjBptK6zGk4ay6TY481OkIXA2qXH/0ax4bnILyxmkS1dXVQ+xsSrnVtKuwNnA7IgYKB9iY1XN7eDAtvOBbUn3WX2FNP72FZKOLRlsw1TKbUS8AHwGeB64FngceACYAXwuIuaXD7lVxuf5U13WPUVqbRzfZZ29jZE8CHMdtgC6/ZMDrO3YZl1jCFbdfzSrIzcXA/uSrvUvrRBb05XI7Y+Bx0g3EttbquZ2qzxfBRyQL1MhaSHpcsq5ki6PdnZGXOK8fZl0/9D1wG2k4nUGcJWkQyPi1kKxttHgJetuf6O1Q7axDdC2FqtXgG6DMgNs1rFNXfuPZkVzI2kO6ZLVTyPivIqxNV2l3OZLUgcBJ0TE64Vja7qq5+2aPL86OgZ1z60t1wM78larVttUPW/3IhVTt0bEqRGxMCIuI93T9jQwP7dq2/AM5r7b36jtn2eVtK2wWk5qfu52Ik0gNVv/v29PVfcfzYrlRtJZwOmkR6pPKBZhcw07t3mfuaQhnZ6WNEXSFGBi3mTrvGybOgJvgKrn7ZN5/nSXdSvy/F0V4muyqrmdSfqAX9C5MCJeAf5AOocnlQm1lQYfCujWJcgE0pOBy7uss7fRtsJqCemYP9S5UNJmwPuBO2refzQrkhtJZwJnAr8Ejov8/G/LVcnt5sD2wMHAQx3Torx+ev79uKIRN0fV83bwxuyduqwbXPZMlQAbrGpuBz/wu7VKbTxkbhtuSZ7v12Xdh4Gl4aHfhqVthdU1pCr860OWf5V0LfnKwQWSJkuaNtz9W6hqbpF0BulJnyuAY1t6X0o3VXK7Gjiiy3RiXn9z/v36WiIf+aqet78j3V81vfMJKknjSI+zPxQRD9cReANUze39eX5M58LcunooqVPMRwrGO2rlJyqnSRrTsfg60qXskzovqUo6BJhMuz/PKmndWIGS5pHu3VlIujyyB6kn4L8CBw5+mGvdAzuv1/5tVCW3kmYAlwDLgO+SHtPutLLNN6pWPW+7vN4k0s3sP4yIVne/UOA94XhSj9X/AH5G6sDya8A44NMRcUtvjmTkqfieMBG4i3Qp9cq8z7akwmwSMCMiftSrYxlpJB3NW5f0Tyaddxfm3x+PiCs6tl0EfALYtfOpYEmzgAtILdhXk1oJZ5E6vN3HLVbD1O8eSns9kZqVZ5GGS3mV9FjpXGDskO0GUnqGt38bpyq5BX5B+na7rmlRv4+vqbldx+tNwj2vF8stcDip75/VpBasW4CP9PvY+j0VeL+dTBpy5UngdeDfwJ+Bw/t9bP2eSMXQer1fdmw7qcvrHAPcS3oS8BnSl4Md+n18TZ5a12JlZmZmVpe23WNlZmZmVhsXVmZmZmaFuLAyMzMzK8SFlZmZmVkhLqzMzMzMCnFhZWZmZlaICyszMzOzQlxYmZmZmRXiwsrMRoQ8ltmO/Y7DzKwKF1ZmVitJO0u6UdLLkpZK2rfLNt8C/gkMSDqg91GamZXhwsrMaiNpE+Aq4EFgBTAVWJCXd5qa55uSBpQ1M2skjxVoZrXJBdSOEbFM0t7AnXnVIRHx+47t3g3sAPwAeC0iDu59tGZm1bnFysxqExGvRcSy/PNdwJN51f5DtvsXsBLYC1jY0yDNzApyYWVmvXRXnk/tXCjpPcCNwCLgsh7HZGZWjAsrM+ulpXm+E6RLhZJOIV0ivBY4Mnx/gpk12Mb9DsDMWuWBPN9Z0gnAqcBiYP/BS4ZmZk3mFisz66Xn8nw88EngsIj4QmdRJWmspLNzFw3PSorcHYOZ2YjnwsrMaidpF0mX8t/3T82JiPu6bL4dcAbpRva7exGfmVkpLqzMrDaStpT0PWAJ8DgwBXgtr56yjt1WABMiYmfg+PqjNDMrx4WVmdVC0kHAvcA7gT0jYk5EvAQ8kTfpWlhFxKsRsbxHYZqZFeWb182sOEmzgaOA6RGxeMjqx4DJwHt7HpiZWc3cYmVmRUk6ETga2K9LUQXwUJ7v2buozMx6w4WVmZV2OjAPWLWO9UvyfG9JE+DNDkLNzBrPhZWZFSNpU2AccD7wjKStumx2HfAyMAa4SdJc4FeSfGuCmTWeCyszKyYiXiV1+DmG1GL1Rpdtngc+S+osdDdgd+DoiPifbc3MmsbfEM2stIOAfYE7I2JNtw0i4o/AHj2NysysB1xYmVlREbEKuLXfcZiZ9YMLKzMbcSSdBGyTJ4ADOu7Bmpf7wzIzG3HkgeTNbKSRNABMXMfqXSNioHfRmJmtPxdWZmZmZoX4qUAzMzOzQlxYmZmZmRXiwsrMzMysEBdWZmZmZoW4sDIzMzMrxIWVmZmZWSEurMzMzMwKcWFlZmZmVogLKzMzM7NC/gPMjShqE2ltzQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x576 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Notice how HEMM comes close to recover the True Subgroup\n"
     ]
    }
   ],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "\n",
    "plt.figure(figsize=(10,8))\n",
    "\n",
    "axis = plt.gca()\n",
    "\n",
    "plt.imshow(Zgrid**5, cmap='rainbow', interpolation='nearest', origin='lower', extent=[0, 1, 0, 1])\n",
    "\n",
    "truesg = plt.Circle((0.5, 0.5), 0.2, color='k', lw='3', ls='--', fc='none', label='True Subgroup')\n",
    "plt.gca().add_artist(truesg)\n",
    "plt.legend(handles=[truesg], fontsize=20)\n",
    "plt.xlabel(r'$\\mathcal{X}_1$', size=24)\n",
    "plt.ylabel(r'$\\mathcal{X}_2$', size=24)\n",
    "plt.xticks(size=18)\n",
    "plt.yticks(size=18)\n",
    "plt.title(\"The Discovered Subgroup vs. True Subgroup\", size=20)\n",
    "plt.colorbar(orientation='vertical')\n",
    "plt.show()\n",
    "\n",
    "print  (\"Notice how HEMM comes close to recover the True Subgroup\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='syndatabs'></a>\n",
    "### 2.5 Bootstrapping PEHE Estimates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Till now we have only looked at Point Estimates for the Counterfactual Outcomes and the Associated **PEHE** score. For fair comparison we would also like to estimate the uncertainity around our estimates of **PEHE** by computing a **Bootstrapped Confidence Interval**. We will compute this by fitting multiple models on each bootstrapped realization of the Synthetic Dataset. \n",
    "\n",
    "We will wrap the previous experiment into a single function, and use `joblib` to parallelly estimate the **PEHE** in each bootstrapped realization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "###################################################################################\n",
    "# NOTE: For more complicated decision problems, linear adjustment for confounding #\n",
    "# may not be sufficient. HEMM allows passing any PyTorch model for determining    #\n",
    "# the outcome as a function of the features independent of the Treatment          # \n",
    "# assignment. The parameters of this 'Outcome Model' and the underlying,          #\n",
    "# 'Heterogenous Effect Mixture Model' are optimized jointly. HEMM implementation  #\n",
    "# in Causllib, has helper functions to instantiate Neural Network outcome models. # \n",
    "###################################################################################\n",
    "\n",
    "from causallib.contrib.hemm.outcome_models import genMLPModule, genLinearModule, BalancedNet\n",
    "\n",
    "\n",
    "def experiment(data, i, K=2, response='bin', outcomeModel='linear', lr=1e-3, batch_size=100, vsize=0.3, bc=2, epochs=20):\n",
    "    \n",
    "    vsize = int(vsize*data['TRAIN']['x'].shape[0])\n",
    "    \n",
    "    Xtr  = data['TRAIN']['x'] [:-vsize,:, i]\n",
    "    Ytr  = data['TRAIN']['yf'][:-vsize  , i]\n",
    "    Ttr  = data['TRAIN']['t'] [:-vsize  , i]\n",
    "\n",
    "    Xdev = data['TRAIN']['x'] [-vsize:,:, i]\n",
    "    Ydev = data['TRAIN']['yf'][-vsize:  , i]\n",
    "    Tdev = data['TRAIN']['t'] [-vsize:  , i]\n",
    "\n",
    "    Xte  = data['TEST']['x'] [:,:, i]\n",
    "    Yte  = data['TEST']['yf'][:  , i]\n",
    "    Tte  = data['TEST']['t'] [:  , i]\n",
    "    \n",
    "    mu1tr = data['TRAIN']['mu1'][:-vsize, i] \n",
    "    mu0tr = data['TRAIN']['mu0'][:-vsize, i]\n",
    "\n",
    "    mu1te = data['TEST']['mu1'][:, i] \n",
    "    mu0te = data['TEST']['mu0'][:, i]\n",
    "\n",
    "    mu  = Xtr.mean(axis=0).reshape(1,-1)\n",
    "    std = Xtr.std(axis=0).reshape(1,-1)\n",
    "        \n",
    "\n",
    "    #Set the preferred Outcome Adjustment\n",
    "    if outcomeModel == 'MLP':\n",
    "        outcomeModel = genMLPModule(Xte.shape[1], Xte.shape[1], 2 )\n",
    "    elif outcomeModel == 'linear':\n",
    "        outcomeModel = genLinearModule(Xte.shape[1], 2 )\n",
    "    elif outcomeModel == 'CF':\n",
    "        outcomeModel = BalancedNet(Xte.shape[1], Xte.shape[1], 1 )\n",
    "        \n",
    "    Xdim = Xte.shape[1]\n",
    "        \n",
    "    model = HEMM(Xdim, K, homo=True, mu=mu, std=std, bc=bc, lamb=0.0000,\\\n",
    "                spread=.01,outcome_model=outcome_model,sep_heads=True,epochs=epochs,\\\n",
    "                 learning_rate=learning_rate,weight_decay=0.0001,metric='LL', use_p_correction=False,\\\n",
    "                 response=response,imb_fun=None,batch_size=batch_size )\n",
    "\n",
    "    \n",
    "\n",
    "    cd = model.fit(Xtr, Ttr,Ytr, validation_data=(Xdev, Tdev, Ydev))\n",
    "    \n",
    "    inSampleCFoutcomes  = model.estimate_individual_outcome(Xtr, Ttr)\n",
    "    outSampleCFoutcomes = model.estimate_individual_outcome(Xte, Tte)\n",
    "    \n",
    "    inSampleITE = inSampleCFoutcomes[1]-inSampleCFoutcomes[0]\n",
    "    outSampleITE = outSampleCFoutcomes[1]-outSampleCFoutcomes[0]\n",
    "\n",
    "    #Compute the In-Sample PEHE\n",
    "    trueCATE = mu1tr - mu0tr\n",
    "    inSamplePEHE = np.sqrt(((inSampleITE - trueCATE)**2).mean())\n",
    "\n",
    "    #Compute the Out of Sample PEHE\n",
    "    trueCATE = mu1te - mu0te\n",
    "    outSamplePEHE = np.sqrt(((outSampleITE - trueCATE)**2).mean())\n",
    "\n",
    "\n",
    "    return [inSamplePEHE, outSamplePEHE]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "from joblib import Parallel, delayed\n",
    "\n",
    "PEHEs = Parallel(n_jobs=10)(delayed(experiment)(data=syn_data, i=i, K=3,lr=1e-4) for i in range(100))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "In  Sample Root PEHE 0.10267147617234845 , Std. Error 0.0001256753369609759\n",
      "Out Sample Root PEHE 0.1026503798822282 , Std. Error 0.0001542102418122442\n"
     ]
    }
   ],
   "source": [
    "PEHEmean =  np.mean(PEHEs,axis=0)\n",
    "\n",
    "PEHEstd =  np.std(PEHEs,axis=0)\n",
    "\n",
    "print (\"In  Sample Root PEHE\", PEHEmean[0], ',', \"Std. Error\", PEHEstd[0]/np.sqrt(100) )\n",
    "print (\"Out Sample Root PEHE\", PEHEmean[1],',', \"Std. Error\", PEHEstd[1]/np.sqrt(100) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 708,
   "metadata": {},
   "outputs": [],
   "source": [
    "#syn_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='IHDPdata'></a>\n",
    "## 3. IHDP Example\n",
    "In this section, we will experiment with the IHDP dataset. As opposed to the previous dataset the IHDP dataset comes from a real world study, with simulated outcomes and hence is considered 'Semi-Synthetc'.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='IHDPdatadesc'></a>\n",
    "### 3.1 Data Description"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **UN Infant Health Development Program** (**IHDP**) dataset has gained popularity in the causal inference literature dealing with heterogenous treatment effects. The original data includes 25 real covariates and comes from a randomized experiment to evaluate the benefit of **IHDP** on IQ scores of three-year-old children. A selection bias was introduced by removing some of the treated population, thus resulting in **608 control** patients and **139 treated** (**747 total**). The outcomes were simulated using the standard non-linear **‘Response Surface B’** as described in \n",
    "***Hill, J. L. (2011). Bayesian nonparametric modeling for causal inference. Journal of Computational and Graphical Statistics, 20(1), 217-240.***\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**`causallib`**'s **`HEMM`** module contains helper functions to download the **IHDP** dataset courtesy of _Frederik Johansson's_ [Personal Website](https://www.fredjo.com) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "IHDP Data exists\n"
     ]
    }
   ],
   "source": [
    "from causallib.contrib.hemm.load_ihdp_data import loadIHDPData\n",
    "\n",
    "ihdp_data = loadIHDPData()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "from joblib import Parallel, delayed\n",
    "\n",
    "PEHEs = Parallel(n_jobs=10)(delayed(experiment)(data=ihdp_data, i=i,outcomeModel='CF', K=3,lr=1e-3,vsize=0.25, batch_size=10, bc=6, response='cont', epochs=500) for i in range(10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "In  Sample Root PEHE 4.779168419249578 , Std. Error 0.8265074053983431\n",
      "Out Sample Root PEHE 4.300757764814082 , Std. Error 0.7104830767106435\n"
     ]
    }
   ],
   "source": [
    "PEHEmean =  np.mean(PEHEs,axis=0)\n",
    "\n",
    "PEHEstd =  np.std(PEHEs,axis=0)\n",
    "\n",
    "print (\"In  Sample Root PEHE\", PEHEmean[0], ',', \"Std. Error\", PEHEstd[0]/np.sqrt(100) )\n",
    "print (\"Out Sample Root PEHE\", PEHEmean[1],',', \"Std. Error\", PEHEstd[1]/np.sqrt(100) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.7.5 64-bit",
   "language": "python",
   "name": "python37564bitd905ce786ede48728f7c06fb20d97845"
  },
  "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.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}