{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Weighted Least Squares"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true,
    "execution": {
     "iopub.execute_input": "2021-02-02T06:51:10.187151Z",
     "iopub.status.busy": "2021-02-02T06:51:10.186254Z",
     "iopub.status.idle": "2021-02-02T06:51:10.581848Z",
     "shell.execute_reply": "2021-02-02T06:51:10.583212Z"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:51:10.588108Z",
     "iopub.status.busy": "2021-02-02T06:51:10.586746Z",
     "iopub.status.idle": "2021-02-02T06:51:11.485663Z",
     "shell.execute_reply": "2021-02-02T06:51:11.486823Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import stats\n",
    "import statsmodels.api as sm\n",
    "import matplotlib.pyplot as plt\n",
    "from statsmodels.sandbox.regression.predstd import wls_prediction_std\n",
    "from statsmodels.iolib.table import (SimpleTable, default_txt_fmt)\n",
    "np.random.seed(1024)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## WLS Estimation\n",
    "\n",
    "### Artificial data: Heteroscedasticity 2 groups \n",
    "\n",
    "Model assumptions:\n",
    "\n",
    " * Misspecification: true model is quadratic, estimate only linear\n",
    " * Independent noise/error term\n",
    " * Two groups for error variance, low and high variance groups"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:51:11.492315Z",
     "iopub.status.busy": "2021-02-02T06:51:11.490816Z",
     "iopub.status.idle": "2021-02-02T06:51:11.501567Z",
     "shell.execute_reply": "2021-02-02T06:51:11.502645Z"
    }
   },
   "outputs": [],
   "source": [
    "nsample = 50\n",
    "x = np.linspace(0, 20, nsample)\n",
    "X = np.column_stack((x, (x - 5)**2))\n",
    "X = sm.add_constant(X)\n",
    "beta = [5., 0.5, -0.01]\n",
    "sig = 0.5\n",
    "w = np.ones(nsample)\n",
    "w[nsample * 6//10:] = 3\n",
    "y_true = np.dot(X, beta)\n",
    "e = np.random.normal(size=nsample)\n",
    "y = y_true + sig * w * e \n",
    "X = X[:,[0,1]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### WLS knowing the true variance ratio of heteroscedasticity\n",
    "\n",
    "In this example, `w` is the standard deviation of the error.  `WLS` requires that the weights are proportional to the inverse of the error variance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:51:11.507524Z",
     "iopub.status.busy": "2021-02-02T06:51:11.505921Z",
     "iopub.status.idle": "2021-02-02T06:51:11.539939Z",
     "shell.execute_reply": "2021-02-02T06:51:11.539296Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            WLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.927\n",
      "Model:                            WLS   Adj. R-squared:                  0.926\n",
      "Method:                 Least Squares   F-statistic:                     613.2\n",
      "Date:                Tue, 02 Feb 2021   Prob (F-statistic):           5.44e-29\n",
      "Time:                        06:51:11   Log-Likelihood:                -51.136\n",
      "No. Observations:                  50   AIC:                             106.3\n",
      "Df Residuals:                      48   BIC:                             110.1\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          5.2469      0.143     36.790      0.000       4.960       5.534\n",
      "x1             0.4466      0.018     24.764      0.000       0.410       0.483\n",
      "==============================================================================\n",
      "Omnibus:                        0.407   Durbin-Watson:                   2.317\n",
      "Prob(Omnibus):                  0.816   Jarque-Bera (JB):                0.103\n",
      "Skew:                          -0.104   Prob(JB):                        0.950\n",
      "Kurtosis:                       3.075   Cond. No.                         14.6\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "mod_wls = sm.WLS(y, X, weights=1./(w ** 2))\n",
    "res_wls = mod_wls.fit()\n",
    "print(res_wls.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## OLS vs. WLS\n",
    "\n",
    "Estimate an OLS model for comparison: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:51:11.544452Z",
     "iopub.status.busy": "2021-02-02T06:51:11.543582Z",
     "iopub.status.idle": "2021-02-02T06:51:11.553387Z",
     "shell.execute_reply": "2021-02-02T06:51:11.554464Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[5.24256099 0.43486879]\n",
      "[5.24685499 0.44658241]\n"
     ]
    }
   ],
   "source": [
    "res_ols = sm.OLS(y, X).fit()\n",
    "print(res_ols.params)\n",
    "print(res_wls.params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare the WLS standard errors to  heteroscedasticity corrected OLS standard errors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:51:11.560610Z",
     "iopub.status.busy": "2021-02-02T06:51:11.557829Z",
     "iopub.status.idle": "2021-02-02T06:51:11.569966Z",
     "shell.execute_reply": "2021-02-02T06:51:11.570908Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=====================\n",
      "          x1   const \n",
      "---------------------\n",
      "WLS     0.1426  0.018\n",
      "OLS     0.2707 0.0233\n",
      "OLS_HC0  0.194 0.0281\n",
      "OLS_HC1  0.198 0.0287\n",
      "OLS_HC3 0.2003  0.029\n",
      "OLS_HC3  0.207   0.03\n",
      "---------------------\n"
     ]
    }
   ],
   "source": [
    "se = np.vstack([[res_wls.bse], [res_ols.bse], [res_ols.HC0_se], \n",
    "                [res_ols.HC1_se], [res_ols.HC2_se], [res_ols.HC3_se]])\n",
    "se = np.round(se,4)\n",
    "colnames = ['x1', 'const']\n",
    "rownames = ['WLS', 'OLS', 'OLS_HC0', 'OLS_HC1', 'OLS_HC3', 'OLS_HC3']\n",
    "tabl = SimpleTable(se, colnames, rownames, txt_fmt=default_txt_fmt)\n",
    "print(tabl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate OLS prediction interval:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:51:11.575334Z",
     "iopub.status.busy": "2021-02-02T06:51:11.573993Z",
     "iopub.status.idle": "2021-02-02T06:51:11.580862Z",
     "shell.execute_reply": "2021-02-02T06:51:11.581747Z"
    }
   },
   "outputs": [],
   "source": [
    "covb = res_ols.cov_params()\n",
    "prediction_var = res_ols.mse_resid + (X * np.dot(covb,X.T).T).sum(1)\n",
    "prediction_std = np.sqrt(prediction_var)\n",
    "tppf = stats.t.ppf(0.975, res_ols.df_resid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:51:11.585679Z",
     "iopub.status.busy": "2021-02-02T06:51:11.584440Z",
     "iopub.status.idle": "2021-02-02T06:51:11.589980Z",
     "shell.execute_reply": "2021-02-02T06:51:11.590860Z"
    }
   },
   "outputs": [],
   "source": [
    "prstd_ols, iv_l_ols, iv_u_ols = wls_prediction_std(res_ols)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Draw a plot to compare predicted values in WLS and OLS:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:51:11.594776Z",
     "iopub.status.busy": "2021-02-02T06:51:11.593508Z",
     "iopub.status.idle": "2021-02-02T06:51:12.004752Z",
     "shell.execute_reply": "2021-02-02T06:51:12.005693Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fd037f68850>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAFnCAYAAAB+YZr1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAACRmUlEQVR4nOzdd1hURxfA4d9l6aKiYkUQe+8VjYq99941UWOL3Wg+001sIPbejcbYUKPG2LGBBcWOXRTBiiAibcv9/hjRmFgQFnbReZ/HJ7rs3jug2bMzc+YcRVVVJEmSJElKWxamHoAkSZIkfYpkAJYkSZIkE5ABWJIkSZJMQAZgSZIkSTIBGYAlSZIkyQRkAJYkSZIkE3hvAFYUZZmiKA8VRbnwj8fKKYpyTFGUM4qiBCiKUiV1hylJkiRJH5ekzIBXAI3/9dhU4CdVVcsB37/4syRJkiRJSWT5vieoqnpIURS3fz8MZHrx+8xAWFJu5uTkpLq5/ftSkiRJkvRxOnXq1GNVVbO/6WvvDcBvMRzYpSiKF2IWXf1tT1QUpT/QH8DV1ZWAgIBk3lKSJEmS0hdFUW6/7WvJTcIaCIxQVdUFGAEsfdsTVVVdpKpqJVVVK2XP/sYPAZIkSZL0yUluAO4F+Lz4/QZAJmFJkiRJ0gdIbgAOA2q/+H1d4JpxhiNJkiRJn4b37gErirIW8ACcFEW5C/wA9ANmKopiCcTxYo83ObRaLXfv3iUuLi65l0g3bG1tyZs3L1ZWVqYeiiRJkmRiScmC7vKWL1U0xgDu3r1LxowZcXNzQ1EUY1zSLKmqSnh4OHfv3iV//vymHo4kSZJkYiavhBUXF0e2bNk+6uALoCgK2bJl+yRm+pIkSdL7mTwAAx998E30qXyfkiRJ0vuZRQCWJEmSpE9NugvAWwJDqTF5P/nH7aDG5P1sCQxN8TU1Gg3lypWjZMmSlC1bFm9vbwwGwztfExwczO+//57ie0uSJEmfpuRWwjKJLYGhfONznlitHoDQyFi+8TkPQOvyzsm+rp2dHWfOnAHg4cOHdO3aladPn/LTTz+99TWJAbhr167Jvq8kSZKUNrYEhuK56wphkbHkcbRjTKOiKYobxpCuZsCeu668DL6JYrV6PHddMdo9cuTIwaJFi5gzZw6qqhIcHEzNmjWpUKECFSpUwM/PD4Bx48Zx+PBhypUrx/Tp09/6PEmSJMm0EidvoZGxqLyavBljBTUl0tUMOCwy9oMeT64CBQpgMBh4+PAhOXLkYM+ePdja2nLt2jW6dOlCQEAAkydPxsvLi+3btwMQExPzxudJkiRJpvWuyZspZ8HpKgDncbQj9A3BNo+jndHvpaoqIAqFDBkyhDNnzqDRaLh69eobn5/U50mSJElpK60mbx8qXS1Bj2lUFDsrzWuP2VlpGNOoqFHvc/PmTTQaDTly5GD69OnkzJmTs2fPEhAQQEJCwhtfk9TnSZIkSWnrbZO01Ji8fYh0FYBbl3dmUtvSODvaoQDOjnZMalvaqEsIjx49YsCAAQwZMgRFUXj69Cm5c+fGwsKC3377Db1eLGNkzJiRZ8+evXzd254nSZIkmVZaTd4+VLpaggYRhI29Zh8bG0u5cuXQarVYWlrSo0cPRo4cCcCgQYNo164dGzZsoE6dOmTIkAGAMmXKYGlpSdmyZendu/dbnydJkiSZVmLMMLcsaCVxrzMtVKpUSf13YlJQUBDFixdPszGY2qf2/UqSJH3KFEU5papqpTd9LV0tQUuSJEnSx0IGYEmSJEkyARmAJUmSJMkEZACWJEmSJBOQAViSJEmSgGvh19L0fjIAS5IkSZ+0M/fPUHtFbYrNLcb1J9fT7L7p7hywsYWHh1OvXj0A7t+/j0ajIXv27ACcOHECa2trUw5PkiRJSgU6g46I2AiyZ8iOnaUdtyNv493Qm9wOudNsDJ98AM6WLdvLVoQ//vgjDg4OjB49+uXXdTodlpaf/I9JkiTpoxCrjWX5meV4+nlSPld5fDr5UNSpKDeG3kBjoXn/BYzIrCLL8OHwIhYaTblyMGPGh72md+/eZM2alcDAQCpUqEDGjBlfC8ylSpVi+/btuLm5sXr1ambNmkVCQgJVq1Zl3rx5aDRp+5coSZIkvVtEbATzA+Yz49gMHsU8olreavQu1/vl19M6+ILcA36rq1evsnfvXqZNm/bW5wQFBbFu3TqOHj36sgvSmjVr0nCUkiRJUlLMODaD8fvHUzFPRQ72Pojf5360LNrSpGMyqxnwh85UU1OHDh3eO5Pdt28fp06donLlyoCoKZ0jR460GJ4kSZL0DlceX8HTz5PWxVrTvEhzvqr6FW2Lt6VsrrKmHtpLZhWAzck/mylYWlpiMBhe/jkuLg4QPYN79erFpEmT0nx8kiRJ0n+dCD3BlKNT2By0GRtLG0rnKA2Ak70TTvZOJh7d6+QSdBK4ublx+vRpAE6fPs2tW7cAqFevHhs3buThw4cAPHnyhNu3b5tsnJIkSZ+yPlv7UHVJVfbf2s//av6P28NvM6zaMFMP663kDDgJ2rVrx6pVqyhXrhyVK1emSJEiAJQoUYJffvmFhg0bYjAYsLKyYu7cueTLl8/EI5YkSfr46Q16Nl/eTPMizbG1tKVe/nqUyl6K/hX7k9Emo6mH916yHWEa+9S+X0mSJGOL08Wx8sxKPP08uRFxg1WtV9GjbA9TD+uN3tWOUM6AJUmSpHRBZ9Dh5efFjGMzePD8AVWcq+DZwJNWxVqZemjJIgOwJEmSZNZitDHYW9mjUTRsCtpEuVzlGFtjLB5uHiiKYurhJZsMwJIkSZJZuhZ+DS8/LzYGbeTqkKtks8+Gby9fMlhneP+L0wEZgCVJkiSzEhAWwJSjU9h0aRPWGmv6lOuDzqAD+GiCL8gALEmSJJmRWxG3qLy4MplsMjHus3EMrTqUXA65TD2sVCEDsCRJkmQyeoMenyAfLjy8wE91fiJ/lvysb7+ehgUbktk2s6mHl6reW4hDUZRliqI8VBTlwr8e/0pRlCuKolxUFGVq6g0x9d29e5dWrVpRuHBhChYsyLBhw0hISMDX15fmzZv/5/nbt2+nfPnylC1blhIlSrBw4UITjFqSJCn9itPFsejUIorNLUbHjR1Zf2k9cTpRZbBDyQ4fffCFpFXCWgE0/ucDiqLUAVoBZVRVLQl4GX9oaUNVVdq2bUvr1q25du0aV69eJTo6mvHjx7/x+Vqtlv79+7Nt2zbOnj1LYGAgHh4eaTtoSZKkdMw32Jf8M/Pz5fYvcbR1ZEOHDVwYeAFbS1tTDy1NvXcJWlXVQ4qiuP3r4YHAZFVV418856HRRvSmYNaxIwwaBDEx0LTpf7/eu7f49fgxtG//+td8fd95u/3792Nra0ufPn0A0Gg0TJ8+nfz581OnTp3/PP/Zs2fodDqyZcsGgI2NDUWLFn3/9yVJkvQJux99n4jYCIpnL06RbEUon6s8o9xHUTd/3XR9lCglklsLughQU1GU44qiHFQUpfLbnqgoSn9FUQIURQl49OhRMm+Xei5evEjFihVfeyxTpky4urpy/fr1/zw/a9astGzZknz58tGlSxfWrFnzWqMGSZIk6ZXrT64zYPsA3Ga4MfivwQDkyZiHv7r9Rb0C9cwj+KZhRch/Sm4SliWQBagGVAbWK4pSQH1DXUtVVRcBi0CUonzvld81Y7W3f/fXnZzeO+N9w/je+A/gbY8DLFmyhPPnz7N37168vLzYs2cPK1as+KD7SpIkfczO3j/LxCMT2XhpI1YWVvQq24sxNcaYelivi4+H1atFL9w//4T8+dP09smdAd8FfFThBGAAzKvPUxKVLFmSf9enjoqKIiQkhIIFC771daVLl2bEiBHs2bOHTZs2pfYwJUmSzJ6qqhhUsSK479Y+/r7+N19X/5rg4cEsbLGQQlkLmXiEL0RGwpQpIuD27QuWlmILM40lNwBvAeoCKIpSBLAG0n70RlCvXj1iYmJYtWoVAHq9nlGjRtG7d2/s7e3/8/zo6Gh8/zHLPnPmjOx+JEnSJ01v0LPx0kaqLKnCqrPivXRApQHcGX6HSfUnmdc53pgYKFQIxo2DUqVg9244fRoqv3UnNdUk5RjSWsAfKKooyl1FUb4AlgEFXhxN+gPo9abl5/RAURQ2b97Mhg0bKFy4MEWKFMHW1paJEycCsG/fPvLmzfvyV2BgIFOnTqVo0aKUK1eOH374QS4/S5L0SfrnUaIOGzoQGRdJJptMANhb2ZvPUaKLF8WMF8RW5sSJIuju3g0NGoCJ9qFlO8I09ql9v5IkfbwarW7E7hu7qZSnEmNrjKVNsTZoLDSmHpagqnDoEHh6wo4dIvBevQrOzmk6jHe1I0zuErQkSZL0ibkffZ/x+8YTGRcJwNgaY9nbYy8n+p6gfYn25hN8r16FatXEsdYTJ+Dnn+HOnTQPvu8jS1FKkiRJ73T9yXW8/LxYcWYFWoOWinkq0rZ4W+rmr2vqob0SGwshIVCkCOTODQYDzJ8PvXqBnZ2pR/dGMgBLkiRJb5SgT6DH5h4vjxL1Lteb0dVHm082M0B4OMybB7NnQ44ccO4cZMwIJ0+aemTvJQOwJEmS9JKqqlx6dImSOUpirbFGVVW+rv41w6oNM69s5uBgmDYNli0Tmc3NmsGYMSZLqEoOGYAlSZIk9AY9m4I2MeXoFM49OMeNoTdwzezK+g7rTT201xkMYGEBR47AwoXQrRuMHg0lS5p6ZB9MJmFJkiR9wv55lKjTxk5EJ0Qzv9l8cmbIaeqhvaKqsGsX1K8PXi96/3TqBLduwfLl6TL4ggzAjBgxghkzZrz8c6NGjejbt+/LP48aNQpvb29KlSr1n9ceO3aMqlWrUq5cOYoXL86PP/6YBiOWJEkyngfRDxj812AcbR3Z2GEjlwZdom+FvthY2ph6aKDVwm+/Qbly0LgxBAVBlizia1ZWZpfV/KE++QBcvXp1/Pz8ADAYDDx+/JiLFy++/Lqfnx81atR442t79erFokWLOHPmDBcuXKBjx45pMmZJkqTkuvfsHmP3jKXDhg4A5HPMx7kB5zjR9wTtSrQzn6NEILrc9ewJOp2Y6d66Bf36mXpURmN2e8AeKzz+81jHkh0ZVHkQMdoYmq75bzvC3uV607tcbx7HPKb9+tfbEfr29n3n/WrUqMGIESMA0RmpVKlS3Lt3j4iICOzt7QkKCiJL4ieuf3n48CG5c+cGRBvDEiVKJOE7lCRJSntXw6/iedSTVedWoTPo6FiyIwn6BKw11hTPbibFge7fh1mzYOBAcHGBYcOga1do0kTs+35kzC4Ap7U8efJgaWnJnTt38PPzw93dndDQUPz9/cmcOTNlypTB2tr6ja8dMWIERYsWxcPDg8aNG9OrVy9sbT+thtKSJJm/TZc20WFDB6w11nxR/gtGuY+iYNa3N5tJc5cvi4zmVavEbLdoUXF+t0oVU48sVclSlEC3bt1o0aIFO3fuZOTIkYSGhuLn50fmzJkJDw9nwIABNG/enAsXLvzntTdu3GD37t388ccfKIryWqOGNzGH71eSpI+bqqrsubkHjaKhXoF6RMZFMs1vGkOqDCGngxklVxkM0LEjbNoEtrZiyXnUKNEs4SMhS1G+R+I+8Pnz5ylVqhTVqlXD39//nfu/iQoWLMjAgQPZt28fZ8+eJTw8PI1GLUmS9DqdQce6C+uouKgijVY3wtPPEwBHW0cm1J1gHsFXr4fDh8XvLSxE1arvv4fbt0Xlqo8o+L6PDMCIfeDt27eTNWtWNBoNWbNmJTIyEn9/f9zd3d/6uh07dpC4gnDt2jU0Gg2Ojo5pNGpJkqRXNl7aSNE5Rem8qTMx2hiWtlzK1s5bTT2sV2Jjxbnd4sWhVi1RsQpEBauffhJVrD4xn/weMEDp0qV5/PgxXbt2fe2x6OhonJyciI6O5sqVK+TNm/fl16dPn86mTZsYMWIE9vb2WFpasmbNGjQaM8oglCTpoxYRG4G1xpoM1hmITogmu312vBp40apYKywUM5lfRUWJxKrZs+HhQ6hUCdatA5m0KveA09qn9v1KkmR8d6PuMt1/OotOL+Inj58Y6T4Sg2pAQUExl1KMWq04qxsRAfnywWefwddfQ+3a6apcZEq9aw9YzoAlSZLSiaBHQXj6ebL63GoMqoFOpTrRoEADAPOZ8Z4+LXrwXr8uWgFmyQI3bkD27KYemdmRAViSJCmdGPzXYI7dPcaXFb9kpPtI8mfJb+ohCaoKu3eLwLtvn+hG9OWXEB8vsptl8H0jswjAqqqaz7JJKkrL5X5JktI3VVXZeX0n049NZ2XrleTJmIf5zeaT1S4r2TOYWUDbuFEcJ8qTB6ZOhf79IXNmU4/K7Jk8ANva2hIeHk62bNk+6iCsqirh4eGyUIckSe+k1WtZd3EdU49O5fzD8+TNlJcbT26QJ2MeijoVNfXwhKgoWLwYnJxEwYyWLUURjU6d4C2Fi6T/MnkAzps3L3fv3uXRo0emHkqqs7W1fS2TWpIk6Z9itbGUnl+aGxE3KJG9BCtaraBL6S5Ya8wkqIWFwcyZ4jjR06fQo4cIwDY24vfSBzF5ALaysiJ/fjPZx5AkSUpj4THh/H39b7qV6YadlR29y/WmbM6yNCvSzHwSq0CUivzmG1FIo317GDNGHCmSks3kAViSJOlTdDvyNt7+3iwJXEKsNpbPXD8jn2M+vq31ramHJqgqHDwIxYpBrlxQqpTY2x05EgoUMPXoPgpm9PFKkiTp4xcaFUrPzT0pNLsQ8wLm0b5Ee84NPEc+x3ymHpqg08H69aIRQp06ojwkQKNGMGeODL5GJGfAkiRJqUxVVZ7GP8XR1hFbS1t239jNkMpDGOE+AtfMrqYe3ivz54ujRLduQeHCsGCB6McrpQoZgCVJklKJQTWw/ep2Jh+ZjNag5UTfE2Szz8adEXfMJ7EqKgoyZRK/P3AAcuYU+70tW4IsrZuq5BK0JEmSkSXoE1geuJxS80rR6o9W3Iu+R6+yvTCoBgDzCL7XrsGAAWJ/99Il8diKFeDnB23ayOCbBuQMWJIkychWn1vNF39+QdmcZVnTdg0dS3bE0sJM3m79/cUy85Yt4sxuz57g4CC+Zm9v0qF9aszkX4QkSVL69SD6AbOOz6JwtsL0LtebrqW74pzRmYYFG5pXgaHISKhXT5SHHD8ehgwRS86SScgALEmSlEw3I27i5efF8jPLidfFM7TqUABsLW1pVKiRiUcHxMWJClUHD8Lq1eDoCH/9BZUrQ4YMaTaMLYGheO66QlhkLHkc7RjTqCityzun2f3NlQzAkiRJyfDzwZ/56eBPWFpY0rNMT0ZXH20+pSLDw2HePHFs6OFDqFhRtAXMmhU8PNJ0KFsCQ/nG5zyxWj0AoZGxfONzHuCTD8IyCUuSJCkJVFVl3819PHouyuZWzlOZ0e6juTXsFotbLjaf4HvkCLi6wvffi0pVBw7AyZMi+JqA564rL4NvolitHs9dV0wynnfxD/Fn0uFJ+If4p8n95AxYkiTpHfQGPT5BPkw5OoVT904xoc4Evq31LU0KN6FJ4SamHp5w8qSY4TZsKGa7vXvDoEFQsqSpR0ZYZOwHPW4qfnf8qLuqLjqDDmuNNft67sPdxT1V7ykDsCRJ0lssPrWYqX5Tuf7kOoWyFmJh84X0LGsmhSkMBti5U2Q0HzwIFSqIAGxnB3Pnmnp0L+VxtCP0DcE2j6Od0e5hjD3mUbtHEa+PB8QxMt9g31QPwO9dglYUZZmiKA8VRbnwhq+NVhRFVRTFKXWGJ0mSlLZita+CxV/X/8LR1pENHTZwefBl+lfsj62lGbQU3bFD1GZu3hxu3hSFMw4cMPWo3mhMo6LYWb1+ptjOSsOYRsZZsk/cYw6NjEXl1R7zlsDQd77uWfwzpvtP52bETQD6VeyHtYU1GkWDtcYaDzcPo4zvXZIyA14BzAFW/fNBRVFcgAbAHeMPS5IkKW2FRoUy/dh0lpxewol+JyiSrQirWq/CwdrBPI4SRUSAhYVodB8fL87wrl4NHTuClZWpR/dWiTPR1MqCftce87/v4R/iz/ar27kbdZc/r/5JZFwkKioj3UfyefnPKe5UHN9gXzzcPFJ99gtJCMCqqh5SFMXtDV+aDnwNbDX2oCRJktJK0KMgPP08WX1uNXpVT6eSndAoYsaW0SajiUcHBAfDjBmwZAmMGgU//QStW4tqVebwwSAJWpd3TrWM56TuMfvd8aP2ytroDDoAPPJ5MLn+ZKrmrfryOe4u7mkSeBMlaw9YUZSWQKiqqmff98lQUZT+QH8AV1czKjouSdInLyo+ikqLK2FQDfSv2J9R7qPIn8VM+pOfPi32dzdsEIG2SxfRhxfETFgC3r/HfOHhBUrlKMXB2wfRG8RMWaNoaFiw4WvB1xQ++G9RURR7YDzwfVKer6rqIlVVK6mqWil79uwfejtJkiSjUVWVv679xbCdwwDIZJOJP9r9wZ3hd5jTdI7pg6+qvvr9zz+Lvd4RI8Q+76pVULq06cZmpt60x2xrpVCn3F1qr6hN6fmlOX73OB5uHtha2v53jzcuTtTATkhI87EnZwZcEMgPJM5+8wKnFUWpoqrqfWMOTpIkyRi0ei1/XPiDqX5TufDwAi6ZXBhfazw5MuSgRdEWph6e2NP9/XeYPh02bRKtAGfNEvu9mTOn6q3Te5WqxLF+v9OHkJgTZLBVibEP4NfjV3DJ5ML0RtMpkb0EGW0ysq/nvld7vDYFxXL+3Lnw6BFkyQKtWqXp2D84AKuqeh7IkfhnRVGCgUqqqj424rgkSZKM4uz9s7T8oyV3nt6hZPaSrGq9is6lOmOlMYPEpchI0XN31iy4d0/McMPDRQBOgy279FKl6n0fEnJku811dSwJ1glE6vUUsCzAb21+o1PJTq/9Pbu7uOOetYxYVVi1Snzwad5c7K3Xrp3m39d7A7CiKGsBD8BJUZS7wA+qqi5N7YFJkiQl16Pnj7j99DaV8lSiUNZClMlZhrlN59K0cFMsFDPZP42LE4H28WOoX18sgzZokKaJVR+SQWwq7/qQULWQBbOOz+L387+ToE9Ar+rRKBq+KP8F3ct0f3URVYXbt8HNTXR8CgyEXr1EIC5WzATflZCULOgu7/m6m9FGI0mSlAK3Im4xzX8aywKX4ZrZlaDBQWSwzsC2LttMPTTh9GnYtg1++EF0JJo6FcqXh3LlTDKc9FCl6k0fEqJ0t/ly+2wilb3oDDo88nnwKOYRCfoErDXW1HGrI56o1YoktmnT4MYNCAmBjBnh2DGz6HcsK2FJkpTuXXp0iV8O/cL6i+uxUCzoUaYHo6uPNo/zu6oKf/8NXl6wf78IAF98AXnzQp8+Jh1aWlSpSqnEDwPxFkHEWZxHUe2JsFqIordiQOUvGOk+kkJZC+Ef4v9qfzdTCfHznjkT7t4Vs1xPT3F2Gswi+IIMwJIkpVOqqqIz6LDSWHHh4QW2Xd3G8GrDGVFtBM6ZzGP5lKAgUSjjwgVwdhYz3v79Uz2xKqnGNCr62vIuGLdKlTHkzmxDUPTvRFr9BhhQsMRB35hi9l8wr1m7l89zd3HHPW81sYR/6hSMGQN16og99iZNzPLolgzAkiSlK3qDni2XtzDl6BRaFGnBd7W/o13xdjQo0IAsdllMPTyRWBUcLJaVXV0hWzZYuRI6d341AzMTqV2lKiUS9An8fv537tpMIjLhKqiAAqqqw07JyfjG1V49+cQJscycKRMsXiwaUly+DEXN54PEm8gALElSuhCni2PV2VV4+Xlx7ck1CmUtRIEsBQDQWGhMH3xv335VscrZWcx+M2QAX1/Tjus9UrNKVXL5BPkwdOdQQp+FUjZnWarnHsmGq3MxqFosFCtG1m5D6zK5YMsWEXiPHBHB96uvXl3EzIMvyAAsSVI60ffPvqw5v4aKuSuyvv162hZvi8bCDPbyLlyAiRNh/fpXFatGjUo3ZSLNxbYr2zgWeozmhZvjaOtIMadiLGu1jAYFGqAoCsNC2r9ep/nHH8U53nz5xPnpL74Q++vpiKL+s/JKKqtUqZIaEBCQZveTJCn9CnsWxoxjMxhYaSD5s+TnzP0zPIl9Qh23OqZPrjIYRIatjQ2sWwf9+sGXX8LQoeDiYtqxpTOXH19mzJ4xbL+6HQUFW0vbN/fivX8f5syBRo2gZk2x4nD8OLRtC5bmO5dUFOWUqqqV3vQ18x21JEmfrOiEaMouKMuT2CeUyF6C/FnyUy5XOVMPSxRuWLNGLHt27Qrjx0O7dtC4sdkkVqUX/iH+TDk6ha1XtmJpYYmCgor63168Fy6In/fvv4sPPRkyiACcL5/4lY6ZX1qYJEmfvMB7gTyOecy69uvoXa63qYcjWgFOmiQKOXzxhWj/V6KE+JqlpQy+SeAf4s/EwxPxD/EHYH7AfA7fOcwPtX9ga+etb67T3KuXqA62fr1YZbh6Fb75xnTfhJHJGbAkSWYnIExsVdV0rWnikbzwxRewebNY/vztN6hXT+7xfoCDwQdp8FsDtAYtNhobDvQ6gGcDT+Y3m08G6wwAok7zjb143DDgnruyeGHVqiKZasAAyJrVhN9B6pABWJIksxNwLwCXTC7kdMhpmgGcPCmWPSdNgvz5RbLPjz9CmTKmGU869TTuKQtPLeSXQ7+gNWgB0Bq0ry8xAzx5gvtvB3CfM1/UxM5ZUdRoHjTIRCNPGzIAS5JkdgZXHkybYm3S9qYGg2j/5+UFhw6JYy1du4oALNsAfjCdQUep+aW4G3WXSnkqcf7BeXQG3etLzNHRMG4cLF8OMTHQsOGrmtifABmAJUkyO9VdqqftDXU6Ubzh3DmRxeztLZadM2VK23GkgrRoN5hYBjKfYz6CHgXxc52fsbSwZEr9KRTNVpSKeSq+KhWZrzbuyotMcXt78WGnY0cYOfKT+6AjjyFJkmRWbkXc4kr4lZcN1FNNeDhs3y4SfQB+/RUKFID27UWS1Ufg352EQJSanNS2tNGCsH+IP3VW1iFeHw+AtcaaswPOUszpX12GdDqxj+7lBdeuicYIGTKIx834GFFKvesYksyCliTJrGy8tJEma5rwPOF56tzg+nUYPFjMdHv3FsEAxJGiLl0+muAL7243aAx3nt6h86bOL4OvgsKY6mNeD77PnommCIULi5lueDhMmPCqIcJHHHzfRwZgSZLMysmwk+R3zE82+2zGvXBIiDizW6SIKBfZpYs4Y1q4sHHvY0ZSo91gvC6eCw8vAJDLIReOto5YWVihUTTYWtrSrHAz8cTE1dWLF2H4cFGe08cHrlwRH4BsU3F1I534dD96SJJklgLCAqiU540rdh9Or4ewMDHbzZxZ9OP95hsYMgRy5zbOPcyYsdoN+of4s/P6Th7HPGbz5c1oFA03h918udz8WivAx7bwTXfx8547F6pVg7NnZQb5G8gALEmS2QiPCedW5C0GVBqQsgvFxIgORN7eolzk+fMioer6dbPpBZsWjNFu8M8rf9JufTt0Bh0AVZyrMLHuRKwsXi3VuztXxf3cE+g1Hg4cAAcH8SEnkQy+byQDsCRJZuPUvVMAVM5TOXkXePBAzLrmzRN7jVWqiL6wqioKZ3xCwRdS1m7QoBqwUCzYHLT5ZfDVKBpaF21NvQL1Xj5vS2AoT0Z/w+f7f+NhJiceDhtPqR9Hg6NjqnxPHxOZBS1JktnQGXRcenSJQlkLYW9ln/QXJgbY334TWc0tW8Lo0VCjhqxY9QFUVeXInSNM9ZtKyewlmVx/Mn53/Kj3Wz20ei3WGmvRKMG2EMyfz8F85Rhw3ZrsD+9SPuwyO4rVxMrWxqhZ1umdbMYgSVK6YGlhSZmcSVyuVFVxhtTLC2rXFgG3UydRvrBIkdQd6EfmyJ0jLAhYwJn7Z7j46CJO9k545PMAoLprdfb33C/2eC3y4/7LCli1CuLiuNzwc2LLt+VOltzcySL21HUvsqzTUwCOjYXAQDhxQvx5+PC0ua8MwJIkmY3/7fsfjQs1pla+Wm9/kk4HmzaJUpEnT4KTk6igBGBtLYPvB0o8x6sz6FBQGOU+ip/r/PzaCoS7izvu3y+CFf8Te+o9e8KIEUxeefON10xJlnVq0+vh0iURbBN/nT8vHgfx+U0GYEmSPin3o+8z6cgknOyd3h2AP/9cLDUXKQILFohgYPdhWb0fi+RWuYqMi2RBwAJaFGmBb7AvBtUAgIViQTa7bCL4arXw55/QurXYOy9TBn74QdRnzpEDgDyO94ySZZ2a7t2DY8de/Tp1Cp6/OGKeObNIExg3DipXFr/y5Em7sckALEmSWTgV9pYErLt3YfZskVXr4iLOkLZvL4r1W3y6pQz+XeUqNDKWb3zOA7w1CN+NusuMYzNYeGoh0QnRaBQNHm4e2GhsSNAniDrNTpXA0xNmzRI/+x07oGlTGDHiP9czRpa1McXHi6XkxGDr7w937oivWVlBuXLi81uVKuJXoUKv/gltCQylw6rULdn5bzIAS5JkFgLCAlBQKJ+7vHjg7FmxzLx2rWiUUKoU9Ogh1gild1a5+nfg8A/xZ+SukZwMOwlAp1KdGFN9DOVylQNetAK8uguPnUG4V2krmiTUqQPz50Pjxm8dQ0qyrI3h/n3w83v169QpSEgQX3N1FUeQhw8X/y1f/u21P5LzYcYYZACWJMksBNwLoHj24jho7MSMa+dOUSt48GAYNkx0JZJeel+VK1VVORl2Ep1eR/3f6hOni0NjoWFd+3W0Ld721QsePBB7vHmqwOCS0KqVaIxQoUKSxtG6vHOaBFy9XhTVOnr0VcC9+WIL2sYGKlWCoUOhenXxGe1DlpI/5MOMMckALEmS6cXHcz/0CpUKuYv9xmLFRGZz//6QJYupR2eW3lblKndma3yCfJh6dCrHQ4/Tt3xfEvQJqKioqsqVx1dENNu2TawwXLkCt2+LffRz50QimxmIixMJUkeOwOHDIuBGRYmv5cwpTpgNGiQCboUKIggnV2qU7EwKGYAlSTKdJ09g4UKYPZuT9+6RcH69eNzb27TjSgcS918j9ReIsziPjaE4imUYwVbbabf+FgWyFGBe03kUdyrOmvNrXu3xng6Hz4uJqmD58sH//vfqoiYMvhERYnabGHADAl4tJ5csKUp316ghfuXPb9zj3cYq2fmhZACWJCntPXokOuIsXfpaI3brkmVNPbJ0o3V5Zy4/OcX4I99iULVYKFbYWdpSLFMhZjaZQtvibdFYiMpf+3ruE+d4nzri3mKQyEBatw7atjVZN6KHD0WgPXhQHOc+d04c7bayEsvJw4fDZ5+JGW42I/fl+DdTJZPJACxJUtp5+lSc/bC0hNWroUMHGDmSmTEHOB66gjU0QNatSpo7T++wLXgKBuJBAUXRMaTqACbVm4SSOD28cAG8vXHPkgX3adNEhAuoItZs07hCWGioCLSJATcoSDxuby+C7E8/Qa1a4rNBWp8qM1UymQzAkiSlLr1enCedNk3Mdk+dEvu6iQ3ZgZ2rv+Ze9L1XgUN6q/MPzuPp58naC2sxGAxoFDHLtdZY06poK/EBZs8e8fPetUtEs8GDxYsVBSpWTJNx3r8v+jL4+or/JrZdzpRJzGx79RLb/BUqmMe2c1olk/2TDMCSJKWOf3Ykun4d3NzEuqJeL2bAL4KvqqoEhAXQulhrU47WbP2z1Z+KSo1lNchglYEhlYcwvNpwwp6FvWoF6OIO338vlvdz5YJffoEBA1J/DRexpOzr+yrgXr4sHs+UScxsBwwQAbdcuU+uJ8ZbyQAsSVLqWL9epKlWqSJ+36bNG/cbbz+9TXhsuPF6ABtJcqtMGdORO0eot0o0QrC1tGVPjz3MajyLbmW6kdUuKwD5DBlxP+IL9tbgAnTrJrKUunZNWWrwe0RFieXkfftg/35RzhFEJ8KaNUXBizp1RMA10Taz2ZM/FkmSjOPSJTHbrVBBBN4uXUSpofd0JAoIEx3SzCkAm6owQ6JYbSwrz65k/L7xJOhFKnCCPoFDtw/xTc1vxJNu3IAZM2DZMrHaoKpiebloUfHLyOLiRGWpffvEr5MnxWKGra1YUu7SRQTcihVFIpX0fu8NwIqiLAOaAw9VVS314jFPoAWQANwA+qiqGpmK45QkyRypqlhv9PIShTNsbUUJIhCzr88+e+8lNIqGqs5VKZ2jdCoPNulMVZgBYNOlTQzcMZBHMY8o7lScaG00eoNeHCFy8xBP6ttXBF5LSzHjHTkSShv352cwiMzk3bvFlvKRIyIIazSiZvK4cVCvHri7v73ClPRuSZkBrwDmAKv+8dge4BtVVXWKokwBvgHGGn94kiSZtQEDYNEiUZz/559h4EDRnegDtCnehjbF26TSAJMnLQsz+If4s+XyFqq7VKdVsVbkzZSXSnkqMbbGWGrlq8Wxu8fwvbkfjzBr3J1flOEsXBi++UbUx86d22hjuXtXBNs9e2DvXnFaDKBECfjySxFwa9USiezGZA7L/abw3gCsquohRVHc/vXY7n/88RjQ3sjjkiTJHEVGwuLFYtaVJ49Yd6xcGbp3T9Y0SFVVDKrh5XlVc5FWhRl+O/sbfbb2Qa/q0SgaDvc5jLuLO391+0s84dkz3H1O4D5jCQQHQ9Yy0KgRjDXOfOf5c7GPu2uXmOkmJk7lzCmOZjdoAPXrg3MqxkJTL/ebkjFaiXwO7HzbFxVF6a8oSoCiKAGPEj9OSZKUvgQHi244Li7w9dewfbt43MNDLIcmcw3y+pPrZJ6cmW1XthltqMYwplFR7Kxe/1BgzMIMB4MP0mRNE3pu6YlefbXU7RvsK34THS2CrIuLyBx3dgYfHxENU0BVxbKyp6e4VNas0KyZ+EyVL5/YSTh7VrTwW71aHBVKzeAL717u/9ilKAlLUZTxgA5Y87bnqKq6CFgEUKlSJTUl95MkKY0ZDGJ2u26d6NvWqROMGiVayxhBQFgAz7XPcc3sapTrGUtqFGbQG/QvZ/rrLq7j9L3TfFnxS1adXfWqTGSWFz9XW1vYuFFMQ0eNSlEHqCdPxOw2cZYbFiYeL1UKvvpKTKhr1jTdPq6p6jC/RqczSap2su+oKEovRHJWPVVVZWCVpI+FwQDHj4vsGgsLcV531CjRaiZvXqPeKiAsAFtLW0pkL2HU6xqDMQoz+If4s+fmHp7FP2Pz5c2sbL2SGq41+KXuL3g38sbW0pZeZXrgu3sRHtsv4D73c7h1SySwXbiQrJJQBoPoibtzJ/z1l/irNBhE7ZMGDUR3wYYNU39mm1SmqsMMiKn+/Pkij+HwYbG3noaSFYAVRWmMSLqqrapqjHGHJEmSSSQWzpg+XZQtunQJihcX65OpJOBeAOVylcNK8/GdW9l1fRfN1zZHZ9ABUNyp+MuvZbXLKlKKVyzB3dsb96AgERETC5XABwXfiAgxu925E/7+Gx48ECe/KlWCb7+FJk3EVr05FsAwSR3m06fFEa4//hCz3+bNQatNvfu9RVKOIa0FPAAnRVHuAj8gsp5tgD0vSscdU1V1QCqOU5Kk1BIRIYLuvHkQHi7eqf/4I9VnA3qDntP3TtO7bO9Uvc/bpGbmrUE10Hlj55fB10KxoHuZ7tRwrfHqSX5+0K+fqFSRWBc7iTUZVVX0xt2xQ2zH+/mJWW7WrGJJuUkT8d8cOYzy7aSqNK/DHBkpik9bWoos/q++SvOZbyIlLVePK1WqpAYEBKTZ/SRJeoe4OLHx9/ixKBNZv75Yav7sszQp1B+jjWHq0anUdK1JvQL1Uv1+//TvzFsQs65JbUsn+43/zP0zrDyzEq+GXmgsNEw+MpmfDv6EVq/FWmPNPo9luK86IMpCTpwooqi/v1jqT8LPOy5OHLnevl0E3tu3xeMVKkDTpuJXlSrmOcs1qagocWb6xAn4/Xfx2K5dYl/d0THVb68oyilVVd9YZUYGYEn6lKiqqBs4bZrIzvH3F2/+4eFpUi/YXNSYvP+N+47OjnYcHVc3ydfxu+PH0sClXHh0gROhJ3CwdsDvcz9K5xRFMfzv+OHruxyPv4JwX3dU7O0OHix+/kkQFiYC7vbtovpUTIzoHtSggchebtrUfPZyzc6NGzB7tgi+z56JD5Y7doji1GnoXQFYlqKUpE9BQoJYVvb2FudMcuYURRwSGyOYIPhef3Kd7PbZyWxr5KoOSWCMzNutl7fSZl0bVMQkZkDFAUysN5EsdllePsd9yd+4T1giipP88IMo0fmOdWFVFX8927aJBlKJ8xU3N1FbuXlz0dBAVp56j7/+Ej8sjQY6d4Zhw8SGuJmRAViSPgWrV8MXX4iSRkuXikL9Jn4X77m5J5YWlhzqcyjN753czNsYbQyXHl2iUp5KXHh44WXw1SgaXDO7kiXBAuZ6iTTjMmXEsa28eaFHjzcmVW0JDGXKjqvcOmePEuKMLjgXj+9boihihXTiRGjRAkqWTPP2velLfLz4gJkxI7RtKz6lfP899O8vCsaYKRmAJeljdPOmyPIsV05Mnbp0EW9EjRqZxTu5zqAj8H4gAysNNMn9PzTzNjwmnLkn5zL7xGwAQkaEUDd/XewO24kzvBZWePx5Fpq5iOXO+Hi26LPhuesBYZHO5Jnp/1piUWQkTJj/hMWrLYi+XhM1wRLFSkeGAuEM7m/Nd4OykDNnqv8Y0r+HD2HBApFA+OABtG4tAnCGDPDjj6Ye3XvJACxJHxN/f7G/uHmzOMP79dficTs7cQDUTASEBRCni6NynsomuX9SM2+3BG1hqt9UTt87Tbw+nuZFmjO2xlhsNDa4u7izr+c+fGcMw2PjKdxDN74sVLJFyfmf8oqjV1xl5/oM3DrlyIEDoNNlxSJDHBmKh2FX+AG2ro+xsDIQaGdHzpxJ34f+ZHl6wnffidlv06biCFcKKoVp9VoOBB+gYcGGxhvje8gALEkfi8GDxUzA0RHGjBHHK8wwQ2fGsRl8vedr7K3sqZmvpsnG8a5CG3qDnhOhJ+i8qTPx+ng0iobVbVbTrUw3sW++ezc0bIi7izvumZpBx9qiUImLCwCek/cTq9WjfexAzLWcxFzLRcI9R24ARYqIZPOVoUexzhP5nwWJNK0AlZ4YDGJvt1o1sadeqBD06SP2d4sVS9Gld13fRb9t/QiJCuHsgLOUyVnGSIN+NxmAJSm9io6G5cvFrCtHDmjVSrwR9ekjuqKbkZOhJ8mbKS+5M+amdI7S9KvQj1HVR5E3k3Era6WEqqrsu7WPqUenUiRbEZwzOr88xwtw5/F18QFn+nS4fl20C6pXTyRXvbyGqPFw8U9Xnl/Jhe6J+Huwzh2BY+3L2Be+z5XFHgAcnhxPaOR/x5EmFaDSk+hoUSBm5kxRIMbLS3yCadNG/Eqm0KhQtAYtbo5uOGdypmDWgsxvNp9SOUoZcfDvJgOwJKU3oaHieMXChWIz0dpa9Ipr2FD8MhOJAW3ykcnsu7WPsTXGMrn+ZOoVqJfm537fxj/En/239qNX9Wy5vIXA+4HkcshFs8LNqOJcBWuNtdjjNSh4fDUNLj17Vaikdm1ATIj9/ESvBB8fuHMHsCiArcsTMlUMxq7wfSwzxgPimFMik1SASk9UVTQdXrgQnj4VWWlr10K7dim67PkH5/Hy9+L387/TvkR71rZbS6kcpTjQ64CRBp50MgBLUnqh14tM5jVrxHJcmzZiJuDubuqR/YdPkA8TD0/k1L1T5HbIjWcDT/pX7G/qYb3GP8SfeqvqEaeLQ0XFNZMrS1osoXuZ7thY2sDTp2KP98Y+PEbNwr1wdVggCpXo9Aq+vqJfwpYtIv/HxkZ8/vnpJ7DI94BJ+8++M7imeQWo9EBV4epVKFpUJAteuSISB0eMEEvPKXAw+CCTjkxi141dZLDKwKBKgxhebbhxxp1MMgBLkjkzGMSaZqVK4kyjTifOkg4bBgUKmHp0r9EZdFhaiLeUjZc2EhUfxeIWi+lRpocIaGbiccxj5p6Yy6OYRyToE1BRsVAs+LLSl3xR/nOxtDxtGly6hPuNG7i7uMPB4WhtHNi/Hzb2Fzlu4eEi2bZpUzEpa9pUnIIRcpPJ0fDe4GqMhg8fBa0WNm0Sy/snT4ol/gIFxGMpKO2l1WuxtLBEURT+uvYXZ+6f4de6vzKg0gBRj9vUVFVNs18VK1ZUJUlKgpgYVV20SFWLF1dVRVHVq1dNPaK3ehr3VJ16ZKqaZ1oe9dz9c6qqquqTmCeqTq8z8ched/PJTXXIjiGq3S92Kj+iDtw+ULX7xU7V/KRR7X6xU/0WfKuqZcqoKqhqrlyq+ssvatyT5+r27arau7eqZskivpQxo6p27aqqmzeLvyYpBZ4+VdXJk1U1b17xwy1cWFXnzFHV6OiUXTbuqTrNb5rq4u2i7ry28+VjsdpYY4z6gwAB6ltiopwBS5I5iYwU53fnzYNHj8Q53lWrRLd0M/Mg+gEzj89k3sl5PI1/Sv0C9TGoBoDXqkGZin+IP77Bvni4ebDq7CoWnV6ERtHQo0wPRlUfRYnsJehRpod4zgNb3DuMhFKl0C1axt4cXfljsw1b8ovtx8yZRY5b+/aiDKSsRJVCCQkid+H5c3GUqGZN0RawaVNxfC6ZQqNCmXV8FgtPLeRp/FNq56tNJhtRejLxv+ZE1oKWJHMQHy82ER88EHUH69UT+7seHmZROOPf4nXx5J2el/CYcNqXaM/YGmOpmKeiqYf10uS9W/nf0Q6oqh4LxYoG+TpS1jk3Q6sOxTmTs6gTPGPGyxKR2gSVwBkHWRBUm81bFCIjRdBt00Y0KapfP8mNiqS3UVVR0HrGDFHUev9+8XhoqFGOy6mqSuHZhbkVeYv2Jdoz2n00lZ1Nc878n2QtaEkyR4mNEby9RfWkQ4dEjebbt82yj9yZ+2dYd2EdE+tNxMbShnlN51E2V1mKZCti6qG9pDPoGLtjITNPfYeqaEEBg6rl9E0NA8oOxfnCbZg2FDZvRrW0JKTZQCb0Ax8fhSdPPMiUScx0O3USM10ZdI0gLk50IZoxA86fF/+2Bw0SSYUaTbKDr6qqHAg+wLLAZSxrtQxrjTWLWizCzdGNAlnMKz/ibWQAlqS09u/GCDlyiCIaiW9IZhR8VVXl4O2DTD4ymV03dpHROiN9K/SlYNaCdCjZwdTDeyleF8+iU4vwPuZNcGQwGjU74u3NgIIlGl1JIkaPg/2r0WbMwv4K3zAmeDDnt+TBwQFathRBt2FDubxsdEuXisYfZcqIc+tduojVnmTS6rVsuLQBLz8vAu8HkjNDTq48vkLpnKWpmz99VRCTAViS0try5aIReIkSsGQJdOtmlu/6wZHBdN7YmeOhx8mRIQcT605kYOWBONo6mnpowKszvHXz16VSnkpM859G3kx5iXnQA1tDZRT9OQo+2kqUXS3uxlZl6YOi3HCsxvTI3hguZqBFC/ihk9h2fEOfBCm5zpwRs9169UQTip49xb91I2yn3I26S41lNbjz9A7FnIqxpMUSupXphq2l+f3/kxQyAEtSart+XVTxqVxZvBl16yaSqsykMcI/JegTuPHkBsWzFye3Q24sLSyZ32w+vcr2ws7KfKLUhosb6LqpKzpVh91hO/b13MfJfifJniE7rb9ZR2PflXQJ/JvMCc/53q4xE2Jr8sDCgGMTCxZ1ETPeV0eGpBQzGETT4unTwddXNC0u9aKiVMaMUKdOsi8d9iyMM/fP0LRwU5wzOtOgQANaFW1FsyLNsFCSn7BlDmQAlqTUoKpw9KhYZt6y5fWeuw4OZtUYAeBZ/DMWnVrE9GPTsbSw5NpX17CxtOHI50eMfq8tgaHJLj5xKuwUU/2msuHihpetABP0CfgG++Lu4s7zXoPYuHoxisHAJtrhzQjO5ihMrpIXmDomKz08zLc1XbrWsaM4s+viAlOnQt++kCVlmfAXHl5gmv801pxbg4O1A2GjwrC1tGVJyyVGGrTpyQAsSamhf3+xvJw1K/zvf2KPN3duU4/qPx49f8Ss47OYc3IOkXGR1HGrQ81cX1B76iHuPY0zenWmLYGh/+kS9I3PeYD33iMgLIDKiyuTySYT3Up3Y2PQRrR6LdZY8ux8ber9CPX3O2LHYLYWGkhwUQWd2x0q5D334nswbvBNyQeJdO/OHXFUbuxYEWj79xfp4m3bgpVVii597sE5xu0dx87rO7G3sufLil8ywn1Eul1mfhcZgCXJGKKiRLJJz55iptu2LZQvD716iXJJZkZVVRRFwf+uP78e/pU2xdswtsZYwh46vwiQccCHBcik8Nx15bXyjACxWj2eu6785/qHbx9mzok55M6YmxmNZ1Axd0WWtlxK+xLtsYmzpt7BvARdXErrK48Yd1dLWCGI/2EibbvA8JcVH1MnQzslHyTStWPHxDLzpk3iz9WqiR68KaxBrjPoiIqPIqtdVvQGPafunWJCnQkMrDSQbPbZUj5uMyUDsCSlxO3bMGsWLF4sjhJlzSqCbpMmph7ZG529f5YpR6dQOGthfqrzE82LNOfKkCsUzlYYgBq/7U9ygEyOt7Xa++fj0QnRfLv/W2Ydn/WyTGSHEh2o7lKDYo864vuZN9VOzaW3+pDzluU4U8+baevcqeiedlvqH/JB4qMQGyuSqvz9xQHpESNEu0tX1xRdNjohmqWnlzL92HRq5avFqjarKJ+7PCEjQrDWfPxnwGQAlqTk0Omge3dRjR/EGZYRI0TNZjOTeJRoytEp/H39bxysHfi6+tcAWCgWL4MvJC1ApkQeRztC33CtxBZ8my5tot+2fkTERbz8moLCj4v2cuP3GoTdsiSY+YTlqcTtoaMoP7IOpa3SPpEttX9OZiEyEg4fhhYtRJp46dLQtSv07p3idpf3nt1j9onZzA+YT2RcJDVda9KxZMeXX/8Ugi9A+k4hk6S0pNfD8ePi95aWomTeyJFw65boUGSGwRdg/P7x1FlZh9P3TvNr3V+5M/wO39X+7o3PfVsvWmP1qB3TqCh2Vq8X19dYPaDHZ2J/r3C2wtR2q800j4VYYYdisECToNJ/5SKKFNCxeJUtDncvUy50B5XH1sXSBMEXUv/nZFLXr4vZbd68Yivl4UPx+MKF4jyvEXpNe/p5MvnIZOrlr8exL45xqM8hmhdpnuLrpjdyBixJ75PY+H7GDAgOFmUM3dxEdR8zFK+LZ835NXzm+hlFshWhQ4kOuGZ2TdJRotTuUZu4PPv9Th+CY3ahWIcSrZ7j6KOufJ5QnZv+ZdD8vp7AHRtZkNuVB25XqPYgM5U69KXDLwlgbwlkNspYUuKj7OV77RqMHg3btokPmF27wvDhKS4Mo6oqh24fwtPPkxHVRlCvQD2+rvE1gysPpmDWgsYZezolA7Akvc3jx+Dl9arxvbs7TJ4sZgZmKCo+6uVRorBnYXxf63t+qvMT5XOXp3zu8km6Rlr0qA1J8OGibjgGKwOo0NS5B9a+U3DuL37k3Rx385u+K3GWRbHtsVAUczCzShkfTS/f+HjR9CNvXvEzPnkSvv1WlIrMlStFl9YZdPgE+eDl58XJsJM42Tvx8LmYTedySNm1PxayGYMk/VtMjCgkEBYGBQuKPbARI8yy8X2iCQcnMM1/Gk/jn1I3f13G1hhLgwINUMyk0Mc/+7LWXF6TI3denC82aMi8fyQ/+elwKpkDx0njaNTAgOWBPaIYcwo640jv8OgRLFggjhKVKgV79ojHdTox+zWCGstq4BfiR+GshRnlPoqeZXuaVTGXtCKbMUjS+xgMsGOHaMRuYSGaJOTJIzq1ZDWDxt1vEPI0hLyZ8qIoCo9iHtGgYAO+rv61WXSASWwFWMW5Cucfnsfb35v5TZYQf6kh2j0/QIGWYJGAjQF23JlGdUVBqT0QmgFYiCphkvFdviyKw/z2m2iS0Lix+HCZKAXB90H0A5YFLmN09dFYaawYXHkwo91H07JoSzQWmvdf4BMkA7D0aXv+XPTbnT5d7IG5uMCwYSIgW1iYZfA9fe80U45OYeOljfj28qVmvprMbDzTbGa7/iH+1F1Vl3hd/MtqVc7aWnTrkJGnFyB37vrMiuxM9LPleDyyx73jIBg6VPzsJeNTVfHvWaMR+7u//SbOqw8bJmo0p9Dlx5fx9vdm1dlVJOgTqO5SndputelauqsRBv9xkwFY+rQtWSISTSpXFh2KjFDJJzWoqsq+W/uYcnQKe2/uJZNNJsZUH/PyCJG5BF+AA8EHiNOJQh6oQMAAnu2ahlfplRRZlI0afYpgeWYQHColShZmMr9G6R+F2FhYvVokD44bJ/bSBwyAPn1EH+QkelvFr4jYCHpv7c2fV/7E1tKW3uV6M9J9pFm1pzR3MgBLn5azZ8USXIMG4hxvnz5QsSLUqGF2jRH+KU4XR5dNXbC0sGRK/Sl8WfFLMtuaPhs40YnQEyw+tYR2tnM5sLYO5LUBCx0WBisW2cbRJ5MrFoHh8PAXsBwvjmyZ6bGtdO/+fbG3O3++yGorV+5VHfKMGT+oC8W/K37djYxmpM92oDkty+UmIjaC72t9z+Aqg8mRwXzaaKYXMgBLHz+DAf7+WwTefftEacgyZcTXMmWCzz4z7fjeIFYby8qzK9lyeQs7uu7AzsqO3d13Uzx78XfWxE2r+sT+If4cCD6AjcaGjee3cez+QSziHVmyfBBZ4t1p1+sAxRJG0mxnAO7BK0WX+1GjxAcdKXW1aAGnTr1KHqxdO9kfLhMrfhmI47lmH1GWWzAoUUz+25nW5Z052PugWa2+pDcyAEsfv86dYcMGcHYWx4j6909xp5bUEhEbwfyA+cw8PpOHzx9SxbkKD54/IE/GPO89SpRW9YkT93hfLjNHZ4cj3tTM+AX/G3CRyIZ3mXkwlsY+uQjL2ZQ900fToFVNo91f+geDAf76SxyVW7NGfKCcNUssMRcu/P7Xv0dI5H2iLLfzzHIHBiUKa0NRHBN6cz/OAJjX1kd6JAOw9PG5f18csRg+HBwdxTJz69aiW4sZ7u8mOv/gPNWXVSc6IZpGBRsxtsZYPNw8kvwml9r1iZ/FP2PbyTMs2HmEOEutqKNnsKCO7SDWtc9G9t9qwY9n6XrLi9Bcxfilbl8A7AKeM8k1NP2dkTVnz5/DypWiz/TVq+Ic79WrYlnfCMflDKoBC8UCx0wRhGjXYqevSiZdG2wMJVFQcP4YKn6ZgfcGYEVRlgHNgYeqqpZ68VhWYB3gBgQDHVVVjXjbNSQpTZw/L7KZ16wBrVYsM7dta7aNEQCCHgVx7ck1WhZtSYnsJehbvi+9yvWiXK5yH3yt1KpPHPz4PsN+n8VfD+ej0xuw2LAZTTdrVBKwsbDg1z/nkv3MYyhZkl/bjSHA6fXqRh91kwJTePgQihWDiAioUgXWroV27Yzy4dI/xB9PP0+y22dnYYuF/NykNaN8HNBrXyVtpfuKX2YkKafcVwD/7h4+DtinqmphYN+LP0uSacTFiXZoZcrAunXQrx9cuSKCr5nyC/Gj1R+tKDGvBEP+GoLeoEdjoWF64+nJCr5g3PrE/iH+fL72awqOb0X+Wfn488lkbMLqMTDjbpbPzE8Jmylki+vEltVWFLQtIfbYz59nSaHaJFj+NxB8VE0K0tCWwFBqTN5Pi94zmdpmBFsCQ0VpyKFD4ehR0R6wc+cUBV+DamDr5a18tuwzqi+rjm+wL3kziWpvrcs7M61tfZwd7VAAZ0c7JrUtLT9MGcl7Z8Cqqh5SFMXtXw+3Ajxe/H4l4AuMNebAJOmdYmNFa7S6dcHWFnLmhEmTxP6uGZ7dTXT87nFG7xnNkTtHyGqXle9rfc+QKkOMUqjAGPWJY2NhwqpDTAltjEGJBysDzs9a8WsdT3rUf8T9H34izvs8P3wxD3ulAF83bU58pixMylGK1ory3m5HUtJtCbjD3qlL8Pb3oerdizyxy0S9orUBaP3jj0a7z0++P/HzoZ9xc3RjZuOZfF7+cxysXzVcaF3eWQbcVJLcPeCcqqreA1BV9Z6iKG/NP1cUpT/QH8A1hb0jJem1IxaRkXDnDuTOLYoLmKkEfQIx2hgcbR3RGrTceXqHGY1m8EWFL157o0uplNQnvhRk4H/LdrIjcio6i2jIlQAWBjSKhsGFMtDr155w7BgOdhnZXLYxNroE4qxsibTLBP9YYv4omxSYwuHDVG7didZP7nE3Uw4m1PmCdWUbEY1Vipfzw2PCmR8wn3r56+Hu4k7vcr0pnr047Uu0x9JCpgWlpVT/aauqughYBKIWdGrfT/pIhYTADz+82t9t0UK0AkxhwfjUFJ0QzeJTi/E+5k2LIi2Y12wen7l+xo2hN1Ltje5DZivx8bDBJ4Fftq7lSjZPyHERexsX2rh1YPuDIBL0CVijwWPi72BdAGbPxv1mLp5b/3c2m7jE/NE0KTCFO3dEHfJixcDFhTD7LPxSsxe7i7ij/8cKSXKX829F3MLb35tlZ5YRo43BoBpwd3Enf5b85M+S31jfhfQBkvsu8EBRlNwvZr+5gYfGHJQkAeKIRUSEKCKgqrBpk6icNGwYFDHfajsPnz9k1vFZzDs5j4i4CGrnq02roq1eft3Us4yNx/2Z/acv5/70IDLnNqg5iVxKab6t8xv93T7Dau4C/LP1w7daLjxca+JePgKaNgWNBsfJ+3n+niVmuWT5gU6cEGfUN24UCYPbtoGbG8MHzzbacv6Qv4YwP2A+GkVDtzLdGOU+ilI5Shlj9FIKJPed4E+gFzD5xX+3Gm1EkpRYQm/6dHB1FQk+rq5w757oUmTmvtv/HYtPL6ZN8TZ8Xf1rquatauohodeL46I//r6d04XbgKUBi7Y2fFdsDdUq7qRJTC4Ub29Y2wcMBty//BL3Md+IF+d7dR25xGxEO3fCL7+An584vztiBHz11csvp+RnbVAN7L25l3r566Gx0FAgSwFGu49maNWhOGeSH47MRVKOIa1FJFw5KYpyF/gBEXjXK4ryBXAH6JCag5Q+EQ8eiP3defNECb3y5UW5SFUVlXzMNPieCjvFVL+pjKw2kqp5qzK+1nhGuo+kqJPpg9KDB7B0KcxZe4V7+b2g/HJQ9KCAoiRgl/cyTdfEwc8/iwphgweLFYb8b16SlEvMKRQVJf4dW1pCYKD4UDljBnz++X9KRCbnZx2vi+f387/j5e/FpUeX2NxpM62LtWak+8jU/K6kZJL9gCXTSwywU6bAN9+82t+tVcts6zOrqsqem3uYenQq+27tI5NNJuY1nUe3Mt2SfU1jlZFUVThyROSpbdwI2kZfQsXFWFvY0KRwE3bd2IlWn4C1xpp9vfbjfhc4fFhkkDs6Jnv80jsEB4sKVUuWwOLF0KmTOD5nZSW6FKVQvC6e6cemM+v4LO5F36NszrKMrj6aTiU7YaUx3+IznwLZD1gyPwYD7Nol9r769IGuXeHLL6FNG7Pe3wURfBv81oB9t/aR2yG3UZojGKOM5PPnIkdtyu9Huem4GPurXzBwYE0c6pVCk/FbhhTsQo6VG/Hf4Itv5jg8GvXE3cUdXDBK9STpDfz9xVbKpk2ivWWHDq9aANq+vaZ3UsVoY7C3ssfSwpJlgcsolaMUK1uvpH6B+rJMZDogA7CUtmJjxZGhGTMgKEg0vde/2ONydDTbGViMNob1F9fTs2xPLBQL2hZvS5dSXehepjs2ljYpvn5Kykhevy5W7ZeujCeq7M/gMQkUFV3FtXTu7Yu7y1dif3FhRYiNxb1JE9xHjRJnqE0gRhvDyjMrORpylAXNFyTrKFZaNZ1IEYNBfLh88ADGjIEhQ0TJSCM4c/8MXn5e7Lu1jxtDb2BvZU9A/wAy2cjWjumJDMBS2mrSBA4eFPu7v/0GHTuCtbWpR/VW4THhzDkxh9knZhMeG06+zPmok78OgyoPMup9PrSMZGKDpzlzYOfuBCzcZ2IzcAZYhb18jt6gxzfYV8xyExKgWzcRiI3QhD0lWqxtwf5b+wFoV7wdbYq3+aDXp1XTiQ/29KlYYl69WizpOziIPQA3N/H7FErc9vD082Tvzb04WDvQv0J/4nXx2FvZy+CbDiWlFKUkJd/Zs2Jp+dkz8efx48HXV7RL697dbINvVHwUQ3cOxXWGKz8e/JHqLtU53OcwHm4eyb5mYlnB/ON2UGPyflFW8IWklpGMjBSr9kWKQLNWcQQGwvffWZK/3VLcixTDu74ndoo1GgNYa/V4xL6okTNnjth7NEHwvRp+lSF/DSEyLhKAb2t+y94ee7GztMM32PeDr/eu1QKTuHVLNP7ImxdGjxYZzQ9fnMwsVcoowRcgICyARqsbcfHhRSbXm0zIiBCmNZpGFjvz7OwlvZ+cAUvGZzCIIxbe3rB/v8j67NIFPDygQQNTj+6dImIjyGKXBXsre/6+/jcdS3ZktPtoSuYomaLrvm/W9r4jJ0FBMHu2aIATU2AtGVp6kzHrLS4ODyargwMjIw6Qednv8O1MqqkJ+FZywqPu57jX7Coulsb7gaqqcjTkKF5+Xvx55U+sNdY0L9KcxoUaUyd/HQBquNbA97bvB187tZpOJMv161C0qNjf7dRJrDBUrGiUSz+Lf8bi04t5Fv+MHzx+oFKeSmzutJkmhZoYZdtDMj0ZgCXjioyEatVEM4R00H8XRLA4EHyAqUencvbBWW4Nu4WtpS0XBl3AWmOcGfr79njfdORkVIOiWIY603As7NkDVgX8yDp0HDG2h3kOWBosOR6ynybFW5LZKqM4U1q2LO6j5uLerJkICiYQo42h3qp6HLt7jKx2Wfm21rcMrjyYnA45X3ueZwNPMlhl+ODrm7TetE4HPj4iq/nrr6FQIZHd3Lq1+PduBGHPwph1fBYLAhbwNP4pjQs1RlVVFEWhdbHWRrmHZB5kAJZSLixMFBNo314kUdWrJ8pGtm9v1v139QY9m4I2MfXoVE7dO0XODDkZWnUoeoMIlMYKvpC0WVtiIH76FJYvh687wY0b4n39q1/OMVtXgyhLOxSdgoqKatBz5vv+NFnXXCxzBgWZrDTn84Tn+IX40aBgA+yt7CmVvRTdS3end7neZLB+c5BNbtcnkxQDSdzfnTVLlIwsVUoclbO0FGenjWTlmZX029YPvaqnXfF2jK4+mirOVYx2fcm8yAAsJV9goDhi8ccf4ixjgwaQOTPMnWvqkSXJ0ZCjdNrYicJZC7Oo+SJ6lO2BreXbj4akJPM2KbO2a9fgf/P9+fOsLwnXq1Oo6g06jbzHb/3GY2VVhhrnfyf7mes0v/YjCahYG1Q8CtQT50nt7U0SfB9EP2DuybnMPTmXqPgoQkaEkMshF4tbLk7S638//zsG1UD3Mt2TfM80LwayaRP07g3R0VC7ttgLaNbMKOd3VVXl0O1DZLPPRqkcpaiatyr9K/ZnRLURFMxa8P0XkNI1WYhD+nDnzolqSb6+Yub1+eeiP2lB837DeBL7hHkn52FQDXxf+/uXWaWJ5fre5d97uCBmXUntjfq2109sU5rMkc7MmAHbzvhDz3pgGQeK+P+yWt5qHOlzRIzPxwfatWN78SzMrFyQyPwdGd+q6wcFHmMd3wmNCuXngz+z8uxKEvQJtCzaktHVR1PDpcYHnT9t+FtD7kff59zAcx88hlSjquL8roOD6DF97Rr89JOY8VaoYJRb6Aw6fIJ88PTzJCAsgF5le7Gi9QqjXFsyL7IQh5Ryz5/Dkyfg4iJKFgYHg6enaI5gpmd3E92OvM30Y9NZcnoJz7XP6Viy48s9tYYFGybpGik5pwv/nbXlypCBirqy/NA7C+fPQ/bsUHrgXM5bvJol9ynamaWXCqHMXwCDB7M1bwX2tRvPjgJVRHecOD7o+E1Kj++oqkp0QjQZbTKiV/X8fuF3epfrzYhqI5JddtPDzYPx+8fzOOYxTvZOybqG0STu73p7w/HjInHw99+hcGFxtMhIVpxZwc8Hf+ZW5C0KZS3E/Gbz6VW2l9GuL6UfMgBL73b3rjjCsmgR1KghOrUULCg2J02U5PMhFgYsZPBfg1EUha6luzLafTSlc5b+4OsYI/O2dXlnquV2Zv58mD8bjj2CIjWC+GUhjOpZnG03WtFl0x+gqlgbFPr9byNKsF4ksQFT998ktNDrFas+5ENAcj9E6A16tlzegqefJ5ltM7Or+y5cM7tyf9T9t+7vJlXisa6DwQdpV6Jdiq6VIkuWwIQJYn+3UCGxjdLLeEHx0fNHZLPPhoViwfUn18nlkItpDafRsmjL966+SB8v838HlUzj9GlRHjJ/fjHTrVtX1GlOZKbBV1VVDtw6QNCjIEAcdRlWdRg3h95kZeuVyQq+kPRzum9z/rxYqXep7s/PByaRq/lCqs9uxdUGJQh0/A5bW+hQsgOHDb2YsNfAvjUa3Bv3g8uXYcECIOUfAj709THaGOadnEfROUVpv6E9j2Me06poKxK3rVIafAEq56mMvZV9ss4Dp1hw8KsqbHfvin/rW7eKn/mgQWKlJ4WuhV9jwPYBuM5wZfvV7QD86PEjfl/40aZ4Gxl8P3FyBiy9kvhmpNHAjh2wfbtoj/bVV2/tjmMu9AY9PkE+TPWbSkBYAP0q9GNRi0WUylGKaY2mvff179sbTU7mrarC7t0wbZo4RmRTyB9DzzqgxHMeyPQsEz9+9h2D7+QUMy9XV9zr9sLdMj8MGABOry/JpvT4zYe+ft7JeYzZM4aqzlWZUn8KrYu1NnrAsNJY8ZnrZ4REhRj1uu/k7y+WmX18RKWqNm3gu+/gxx+Nd4sQfzz9PNlyeQtWGit6lulJieyiCIqp+0FL5kMmYUmiStWyZTBzpuhI1KGDaJsGoqqPmVt5ZiUTDk3gRsQNCmctzJjqY96b0fxPSU2wSmoCU1yc2Dr09oaLFyGXczxDBlkSV2kqv/qPR0XFAoXvNXX5YX6QOMb166/wv/8ZZZzJff3V8Kt4+3tTL389OpTsQERsBBcfXfzgxKoPFa+LT/3CEnr9q/3dY8dE3sKAAeLDZZ48xr2VQU+h2YV4GveUQZUHMaTKEHI5mOZ4mGR6MglLerPbt1+1SIuKEnu82bOLr5l54H0S+4QstllQFIWLjy7iZO+EZwPPZO2pJXVv9J8FM97k8WPRAnDOHFGJsFSlSDrOXMhh7UxKNJ1LLgcPpp20JUEbh7VOpeGSfVCivmjY26jRe8eZ0uM3b3t99mzBtFk3hK2Xt2KtsSa/o1jtyGKXhc9cP0vStVMiVYOvXv/quNA334iKYHPmiP1dI5WIjNPFsfrcalafW83f3f/G1tKWLZ22UDBrwWQ1mpA+HXIG/KlSVVFM4MoVMeMdMQKqmP+B/+DIYLz9vVkauJRNHTfRuFBjEvQJWFlYJXuWln/cDt70f4EC3Jrc7L2vv3FDHIdetgxis/lToOmfuJW5w8mobTxLeEaDAg2YkK83VWt1xT/EH1/vr/B4ngP3wZOgbNlkjdlYem/pzcqzK8lql5VBlcRs7d8Vq9JCd5/uuDm68UvdX4xzwcQPl1u2wIULYGcn9nxdXIxyfhdE2dL5AfOZdXwWD54/oHyu8mzosEGe35VeI2fAkjhisWmTiBI+PiLBZMkSUUDexcXUo3uvwHuBePp5sv7ieiwUC7qV6UaBLAWAlFesSu7e6vHj4OUlfpyWltDwc3/2Otfjpj6Wm4+hQf76TLFoSPl5m8G/GwSWwL2cO+7eJ9O8NnOiWG0sv537jS6lupDRJiMti7akUp5K9CnXxyhJVcn1OOYxZ+6fSXkAPn5cLDNv2iT+3LGjWN2xsxNdiYzkxpMblF1Qlufa5zQq2Igx1cdQN39d2YNX+iAyAH/sIiNFF5zZsyEkRByxuHVLzH7TSRN2vUFPyz9a8jTuKSOqjWBYtWHkzWScvqrwYQlWBoPITfPyEh3nMmVW6TL2CDElF1Aub1F2HkwAQINCnY2nKL9tLxQoIH7+hQqJi5jgTfpxzGPmnpjLnJNzeBzzGFtLW3qW7Unb4m3TfCxv4uHmwTf7vuHh84fkyJAjeRcJDBR1yDNnFkUzvvrKqB8uA+8FcuHhBXqU7UGBLAUYUW0EHUp2oEzOMka7h/RpkQH4Y3bnjmg/9/y56EQ0d64ooWemR4gS6Qw6NlzcwOrzq9ncaTPWGms2ddxEkWxFcLR1NPr9krK3Gh8vajF4eYlTKi6uBvpM3cqFzFNZc+8YTmFOtChWB2uNNQn6BKwT9Hjo8sKmJdCqldGWPT+UVq9l2N/DWHFmBbG6WJoXac6Y6mOo6VrTJON5m8TzwIduH6J9ifZJe9GzZ6JodlQUfPstlCsneky3agUZMxplXKqqsvvGbjz9PNl3ax+5HXLTuVRnrDRWTKg7wSj3kD5dMgB/TFRVNLu/fFlkeLq6wpgx4g2pXDlTj+69nic8Z1ngMryPeRMcGUwxp2LceXqHQlkLpXpB+rclWD19CgsXwowZcO8eFK7jT+uBOzilW8XyZyHkt87P3DLf0HvrbeJXzGF61ymEJARQ2lCIB79+DiZqEB8cGYyboxtWGiuuP7lOl1JdGFV91MujMOamYu6KZLDKgG+w7/sDcEiIWFFYtEj8BTVoIP7tK4roMW0kR+8cZdBfgzj34Bx5MuZhav2p9K/YHyuN+TYYkdIXGYA/BgkJoiHC9Olw5oxYdvv8c9Hs/ocfTD26JLkZcZPKiyvzJPYJNVxqMLPxTJoXaY6FYprZ+r174lTW/PliglWrUQStJq1k5d3/cTNSLDP/7NyDb9aFYrl3Ejo7ezaWakD8E2cyWxXgDh9WJtIYDKqBHVd34OnnyfHQ4wQPCyZ3xtz83f1vk/0ck8pKY0X/iv1xc3R79xOXLhUfLlVVdNsaMQKqVjXaOJ7FP+NZwjPyZMxDZtvMGFQDy1stp2vprkbtjiVJILOg079du6BPHxExSpSA4cPFLMAuDXqjptD1J9c5/+A8bYq3QVVVRuwaQYcSHajhWsNkY7pyRSwzr1ol8taadg4hc+MZbL27iHhdPAbVgF7Vo8GCCXsNfHMzDwwdSuOYYlyO/+/nWWdHO46Oq5uqY47TxbHm3Bq8/L24/PgyrpldGVFtBP0q9DNpYpVR6PVi0z1/ftEYIShIJA8OHQr58hntNv/swdukcBPWtlsL8LJmuCQll8yC/thcuiT2cYsVe/XGtHw5NGxosuzaD3Ei9ARTj07FJ8gHJ3snmhVphrXGmhmNZ6TaPd9XRCMgACZPFhnN1tbQvt8tYqv8yLbbv6PeUulcuA0NrxoYoGwjQRGZ1x69x0KPb8Damivjdrzxvh9SKzq5wp6F0X97f8rkLMOatmvoUKJDul0mjdPFEZ0QjZNqBytWiLX/69fhyy9FSc7ixUVpMSO59OgSXn5erD63+mUP3pHVRr78ugy+UmqSATi9SKxrOH26mPV27Ajr1kGRIvD336YeXZKcCjvFqN2jOHj7IJltMjPus3F8VeWrVF/ae1sXIFWFzJHOTJoEe/eCfTE/Phu/h7HtG1IgjyNVlvgwuFBXRhzRk2+yD8TGUnhQC3w7u+Ph5oG7y6ss8pSWifwQtyNvM+PYDO4/v8/admspkKUAgV8GUjpH6XQdMAyqAdfprnSIL8hc7ysQESHOpq9bB22Nl62duOqnKAorzqzgjwt/yB68kknIAJwerF0Lv/wiZr65comuLV9+aepRJUmCPoFn8c/IZp8NRVG4GXET74be9K3Ql4w2xslUfZ9/V7pSVQi/6ETPFQ5E34WcufU0/mUqu3TfchgDAdunsK/nPu4lfIVDt8liStyjB4wYgXuJErzp8FZyakV/qH+ehVYUhW6lu6E36NFYaNL/UZgLF7AoUYKKeSriey0A6tQRR4mqVzfaqk5ivXBPP09+rvMzjQs1ZmyNsXxd42vTt0KUPkkyAJure/cgRw5xfOXSJbCxERuTHTuK35u5p3FPWXRqETOOz6B+gfqsbL2SCrkrcGvYrTTvAJO4DKzqFZ5fykPU8YJowzOiyRZON+9FHNd48XfEtZfPT9An4Bvsi/tn9eA7K9EZJ+e7q0OltEzk+6w8s5LeW3uT0Tojw6sNZ1jVYbhkNv8CKu9kMIjVG29v2LcPtm3DI58H467/zcNV85N/HvhfYrQxrDizAm9/b25E3KBgloLoDDoAstlnM8o9JCk5ZBKWuTl9Wiwzr1sHGzaII0QJCWBllS72d0OjQpl5fCYLTy0kKj6KevnrMbbGWBoUbGCyMVWb4Mu1w05E3nmEIfsxNJGVyZI/G7aVZ3LXsJKKOcrRJjIXv8bsIkFRsdZYse/zg68tMac1rV7LuovryOWQi/oF6hMeE87SwKX0r9g/Vc5CpymtVuzvenuLI3N5RCIb/ftzIuYaVZdUZV37dXQs2THFt1JVlYqLKhJ4P5CqzlUZU31MqnR1kqS3eVcSFqqqptmvihUrqtIb6PWqunmzqtaqpaqgqhkyqOpXX6nqzZumHtkHG75zuGrxk4XaeWNn9VTYKZOOJSpKVadOVVXHbDqVvH4q39qq/IDKD5Zqvm+91eV+p9X933RWDRkdVBVUvxbl1ImLeqh+t4+YbsxxUeo0v2mqi7eLyo+oXTd1NdlYjC4+XvxXq1XVfPlUtXx5VV29+tXjqqpq9VrVYaKDOnD7wGTf5nr4dXXsnrFqvE5cd3PQZvVQ8CHVYDCkZPSSlCxAgPqWmChnwKZkMIhsZr1eJFPpdKJ8Xt++ol2amVNVlUO3D+Hp58lI95HUzV+X+9H3idXGkj+L6foHR0SIOg0zZojfV2t5nuAqXbmvu/Bi4NC9xDh+6zhJ7O3q9WK/sdKbP6SmlVnHZ/H9ge95Gv+U2vlqM6b6GJoUbmL2Z3jf69IlMdvdu1ec87KxEVssuXK9cVVn7fm1FHMqRvnc5T/oNidDT+Lp58mmoE1oFA0Heh0w6ZE2SQJ5DMn83L4tIsSff8K5c2BrKzKb3dxEVX8zpzfo2Xx5M1OPTuVk2Emy22fn4fOHAKnS9zSpfXgfPBCr9/PmiSqFLVuCRcsv2XJ3EbaqLZZYoBoMWOthUK5S4kUrV5q0NGfQoyDyOebD3soeeyt7GhZsyOjqo1O98leqU1XYv18cGdq5U/wb791blEW1sYHcud/60i6lu3zQrcJjwmm3vt3L7Pox1ccwtOpQ8mQ0bp9fSTI283+3/1ioKvj7iwjh4yM++bdvL5ol5Mr1qlB/OuCx0oMjd45QKGsh5jebT6+yvbCzSp3CH287QgSvEp9CQ8HTE+Zv8ych737K9NCwtO8IKpW3YfmJClSKaMnAlRe5EnkD33KZ8ajVE/dyLcQNTBB8VVXlyJ0jTPWbyvar25nTZA6Dqwymb4W+9K3QN83Hkyr8/aF+fZG8NmGCqF7llLRMY61ey4+717A9MJ6oqLxv/NCVoE/g3INzVMpTiax2WXGwdmBaw2n0q9AvzbLrJSml5BJ0WgkIgMqVxdJy//4wZEi6aAMIopPOijMrGF5tOJYWlqw5twZbS9s0SWapMXn/G8/XOjva8XuXukyZIqoT6vL6Qo+GGBQtAD/U+o4f6/wMT56In3OJEjBqlPjQY6JVBlVVXx6DOR56HCd7J4ZUHsLgKoPT/zGYJ09EoQy9Hr77Tnzg3LQJWrT44Kz9Tadu03Fbcez1HmTTDgbEka5JbUtTp7gDC08tZObxmTyLf8bdkXfJZJMpNb4jSTIKuQRtCk+eiDaA8fHw/fdQsSKsWSOymjOkj/KAN57cwNvfm+VnlhOri6Vi7orUyV+HbmW6pdkY3lRJShthz7mdBSn0HaDRUv6rqZzNPJEEgwi+igra+UvB4yfImhXOnxcVw0yURW5QDVgoFiiKwjT/aTyKecTcpnPpXa439lb2JhmT0Vy/LlZ1VqyAmBhRMCOxMUL7JHY1+hfvPTewMZQkzuL8y8eitY8ZvGMkz/7eybOEZ9TLX4+va3xNRms525XSrxQFYEVRRgB9ARU4D/RRVTXOGANLt65cEVX8V64Ub0gtW756Q+ra1dSjS5LIuEj6beuHT5APlhaWdC/d3WSddP5ZYUobnoGnxwrx/GIesH3GkAEwZowldTb9QaYH2XhqHYsBFWsDWEe48eexG7R0LyT68ZpAeEw4807OY0ngEk72O0mODDnY2HEjOTPk/DiOwcyeDcOGiSNy3bqJxgilS6f4smGRsdhYlibWKgAdj7HECb3yhHu6jXQu3pEx1cdQIXcFI3wDkmRayd4AUxTFGRgKVFJVtRSgAToba2Dp0qxZoj7zsmXQubNIsNq6NV2c3zWoBq6Fi2IUmWwyERoVytfVvyZ4WDBLWy1NUfDdEhhKjcn7yT9uBzUm72dLYGiSXzumUVEsnmbivp+esMt+PLfai3Xfdtj/Lx8/T43A1VWhQ0BzHk0JYcsfGWh4qxx5Y35kecX/MeXgnWSPOSVuRdziq7++wnWGK9/7fk+pHKWIio8CIE/GPOk3+Op04nz6hRfZ5HXqwPjxIqlw2TKjBF8QH7ps9eJaD2z+B4C1WpBKNutY226tDL7SRyOlS9CWgJ2iKFrAHghL+ZDSkdhYsaxcsSKULw/16sGPP8LAgaKKVToQr4vn9/O/4+Xvxb1n97gz4g4O1g4c/fyoUeoKJyWJ6m0uXYJ1U525dTQYenUETTwogIUNX1IZwx9roc8gNueugrbpCLYXr0m85au60mnRCOHf7kffp8icIigodCvTjdHuoymZo2Saj8OooqJEB6KZM+HOHdFxa/p0KFVK/DKyMY2KMs4nlme6+lipeQGxBzy+sekKo0hSakh2AFZVNVRRFC/gDhAL7FZVdfe/n6coSn+gP4Crq2tyb2de7t0TZ10WLIDHj0XT+/LloWRJ8SsdeBr39GUyS9izMMrkLMOsJrOw0YiEGWMV9f93HWaAWK0ez11X3hqAL14UibPr14vt8srDN3HSMl6MS4WxRwz8vPsIfF4E+oCTU2Y2la73n+ukRiOEf1NVld03dnMi9ATf1f6OXA65WNxiMQ0KNMA5U9r0AU5VP/8sjhJFRUHt2jBnDjRrlqq3fFXW85tUKespSeYi2VnQiqJkATYBnYBIYAOwUVXV1W97zUeRBT18uAi+Oh00by72vTw80sUyM7zqb3r87nGqLa1G/QL1GVN9DA0KNEiVTjr5x+3gTf/CFODW5NffyC9cEO/3G7bEYlNlBRXq3mHbsElcjfHHY1lN9Ho91gbY97AJ7oMnQdmywH9n2fAqaza13rS1ei1/XPgDL38vzj04h0smFy4NvoSDtUOq3C9NnTkjfraKAv/7HwQHm0WhEklKj1IrC7o+cEtV1UcvbuIDVAfeGoDTJb1etAFs1EicGc2WTXQiGjoUChc29eiS7NyDc3j6eZLROiPzms2jat6qBA0OophTsVS9b1La9F28CEOn+rM/5C8ss9zH/n9biVEeYWFfjMwJQ3B3cce3hCe+Z7fi0W4k7uVbvnat1G6E8G9H7xyly6YuhESFUDJ7SVa0WkGX0l1Sva1iqjIYYNs2Mds9fFg0SWjUCH79Nd18uJSk9CYlAfgOUE1RFHvEEnQ9IJ1Pb/8hKko0uZ81C27efPWG9N13ph5Zkqmqyr5b+/D082T3jd1ksMrAkCpDXn49tYMvvLtNX1AQ/PQTrPPzg151IH8COgVK6lyYtTcXNY9dRlHXwNdf495hBO4dRrz1Pq3LO6fqEuX96PuEx4RTMkdJCmUtRDGnYsxvNj/9l4pMSBAHqadPh2vXwNVVlI10f7HfKoOvJKWalOwBH1cUZSNwGtABgcAiYw3MZJ49E+d2ly4Vv69eHaZMEQlW6czEwxP59sC35HLIxcS6ExlQaQBZ7LJ80DWSWgbybd40O+1atCQbvHLy+/6zWNvosa/lQ4xGBwpYGKCTbwglYyqibJwDrVt/0HiN7crjK3j5ebHq3CqqOlflUJ9D5HTIye4e/0l3SF+0WnF8CMQsN08e+OMPaNcuXZRDlaSPgayEBeKcblgYODuLvd2iRaFaNbHfW7myqUeXZM/in7Hk9BKqu1Snat6q3Iy4yYFbB+hepjs2lh/eQ9jYe6vXrsHPE1TW+O1H+Wwqhvy7yRlfEYXOPLAej6ImYGmwoNr9gegKtuXouLoffA9jCQgL4NfDv7L18lZsLG3oU64PI91HUihr+ikZ+kYXLojZ7sGDIs3c2vqdjREkSUoZWQnrbRISxLnGGTNEAL59W7whXbqULpreJ7r37B6zjs9iwakFRMZF8s1n31A1b1UKZClAgSzJL0KRnAzmN7l1C4ZM8WfnwwWoeY9Bj6tk12RhxI1C9P/jFO2794PsvxKvnMNGLcPtbMVRTHCEyKAaMKgGLC0s8Q/x59DtQ3xX6zsGVxlstObwJqGqohPRtGmi6YednWiMEBMj/r2/ozGCJEmp59MMwI8fiyNEc+fC/fuieMaPP4o3KkhXwffrPV8z8/hMdAYdbYu3ZbT7aKrmrWqUa7/tHG1Sz9fevQs//BLDit2nMHRtBLniUBSVsZey8YNPOLY57Jjn0ZfHGRyxMeTFhuIvX5sWR4gSxeviWX1uNV7+XoysNpJ+FfvRt0JfPi//ORms00fZ0Hc6dAgaNhSz3F9+EY0RsmUz9agk6ZP3aQVgvR40Gjh7ViRTNWokEq0aNjRpS7oPoaoqR0OO4p7XHY2FBid7J/qW78tI95EUzFrQqPdKSgbzm9y/D99PDmfZ+bnoK86mYCcPgq0T0KsqFnrIpFphu2wVdOpEnouP0PqchzckaaW2yLhIFgaIs9D3ou9RPlf5l2d3U6u7U5qIiBAfMC0sYOxYqFVLHKpu2TJdfbiUpI/dxx+ADQaRwTxjhqja4+0NdeuKms1Fiph6dEmmM+jwCfLBy8+Lk2En2dRxE22Lt+XrGl+n2j3flcH8JuHhMN4zmCUXvdGXWQq1YmgSnY/2ty8wpKQ1CfoErK2t8PDaBK7VgbQ/QvRPbda1wTfYlwYFGrCqzSrq5a+XKmeh08yNG+Lf+bJlYnm5QwfxuKK8+r0kSWbj4w3A0dGwapUon3f1qkiwavGiB6yipJvgm6BPYGHAQqYfm86tyFsUzlqYBc0W0KRQk1S/d1KD49OnMMLbnzVHfUkovhxNhVt0f5ibbzbFUOLpPejRg+Kd5+Mb5o+HmwfuLu7/uU9aBNyLDy8y49gMJtefTDb7bEysOxFbS1vK5y6f6vdOdTNniqIwlpai6cfIkVCmjKlHJUnSO3y8WdBffgmLFoks5hEjRGu0xGMX6UCCPgFrjTV6g56ic4qSI0MOxlQfQ8uiLc2mmP/z5yrDZu5jxa0J6HOeBMsErBUL1q/V0upxNhg0CAYPFk3ZTURVVQ7fOczUo1PZcW0H9lb2+HT0oVGhRiYbk1Ho9bB5s1jVKVZMVK9av170mc6Tx9SjkyTphXdlQX8cAVhV4dgxsfw2dixUqCBmvY8eiXO86WhZ8fLjy3j7e7Pz+k6uDLmCvZU9j2Mem1XD9pg4HYPnbWT1zanosgdio7VGa6XFgIpG0TDBoQXfDFgD9qbtdRurjaXuqrocu3sMJ3snhlYZyqDKg8hmn44TkKKjRd7C9OkivXzkSJHdLEmSWfp4jyFptbBxowi8J06Ao6MoJFChglhiTifLzKqqcuTOEbz8vfjzyp/YWtrSq2wvYrWx2FvZmzz4+of44xvsS01XD87tK8awSxXRZbxFHqss/LDThqL34mnSW0OCBqw11nh0+NpkwTdOF8eRO0eoX6A+dlZ2lM9Vnp5letKrXC/srUz7gSDFfv0VvLwgMlJUqvLyglatTD0qSZKSKf0GYINBdCC6eFEE2rlzoVcv0T4nnTl97zS1VtQim102vq/1vVmdO/UP8afuqrrE6xJAZ4O6Yh91q+ZiyMXbtLz2FE2Hjvh+1YNCt28QEhOAi30lHjx2BZe0HWdEbATzA+Yz8/hMHsc8JnhYMC6ZXZjXbF7aDsTYgoLEErOiiCy3+vVh1ChRKEaSpHQtfS9BL1ki9rsaN043x4gAYrQxrDizgojYCMbXGo+qqqy/uJ4WRVuY1SztxpObtFjejaBnx0T7IoOGbrkn8FvukiiHDsKwYWwJ16R5J6J/evj8IVOOTGHR6UVEJ0TTqGAjxtYYi4ebR/rNaFZV0QBk2jTYswf27ROZ+6qarrZTJEn6FPaA04mHzx8y58Qc5p2cR3hsOHXc6rC3516zK+Z/M+ImX274H3vDNqAYwEIxAGBtYcW+zw++lsVcY/L+N54Vdna0S9VSkvG6eGwsbQiNCqXInCIvi5CUzVU21e6Z6rRaWL1aHJW7cEFUqBo6VCQUZvmwGt6SJJmHj3cPOB1ZfW41ff/sS4I+gZZFWzKm+hiqu1Q3i1maf4g/B4IPUClPJVwSGjLiJzjitpWvTlsz7lgct6oW41Drsng0G/yfI0QprZb1IVRV5dDtQ0z1m0q8Lp69PffinMmZ0JGhONo6Gv1+aSaxQIzBAN98AzlywIoV0KWLKBUpSdJHSQbgVJJYsSqLbRZK5ihJpTyV6FW2FyPdR1LUKfWrPCXVxL0+fHu0E6qqA9UCZcVhHCKqc9ylPsXdLLDaPoo8NWtS4y0fFJJbLetD6A16tl7ZytSjUzkeepzs9tkZWnUoBtWAhWKRfoPvjRsim/nAAVGdzcYGjh8XLQHN4IOZJEmpy7zWPj8CeoOejZc24r7UnZrLazLVbyogeu8ubLHQbILv84Tn9N30M98d7oXKi1aAGMhZdTlzttyjzNnNWO3YKsoYviMYjGlUFDur188lG7uU5OLTi2m3vh2PYx4zv9l8bg+/zbe1vjW7pfsk8/MT2fqFC4uz6lWqiNaXAPnyyeArSZ8IOQM2ohVnVjDh0ARuRtykYJaCzG06l15le5l6WG/kfWgmSy/8QOGHNtx2Ar0CFmiwd3Fl4YkgetZNWoec1CglGREbwYKABRTOVpj2JdrTrXQ3stllo23xtmZThCTZDh4EDw+xpztunCycIUmfMBmAU+jh84c42TthoVhwLfwaOTLkwLOBJ62KtjKbYOEf4s+moE3cirhFl5LdiPBvy7oJbdnPD+R8kI8JDSqxt4yCJeXQWxX/4P1bY5WSvBt1l+n+019mNA+uPJj2JdqT0SYjHUqm01rGz5+LwhkGg0ioqlkTli6Fjh3BwcHUo5MkyYRkFnQyJVasWnV2FRs6bKBF0RZo9VqsNOZV7nJ54HL6beuHXtWDCi3PuPLn1tu4u0Mu53WcLpDhP0ueqZ3B/CYTDk5gwqEJGFQDnUt1Zkz1Mek7o/nePZg9W3QlioiApk1hxw5Tj0qSpDQms6CNJLGusJefF9uubsPW0pbe5XpTInsJALMLvv239Wfx6cWgAgpoDJApOitb18fTor0NW898RtAHdDsypsQktdI5SpPZNjPFnIoxoNIARrqPxM3RLdXvn6oWLBCzXZ0OWreG0aNFSVRJkqR/SKdZLKahV/X03NwT/7v+/FD7B+4Mv8OC5guM3oc3ubR6LWvPr+V5wnMAKly0ZfBxsNWBhV5Bo7Gl/5J5tOxgg6KIpeNJbUvj7GiHgpj5pnYBDYNqYMvlLdRYVoOay2uy5PQSADqU7MCsJrPSZ/BVVdi7V2Q1gyiF2q+fqEfu4yODryRJbySXoN8hOiGa5YHLWXthLft77cfW0pbzD85TMGtBs6lY5R/iz+4bu4mMi2RzkA+3o+4w23EQF4Lm8sfiZ/S1Wknc2OJkr32ChoX/2wowraiqyrLAZXj6eXIl/Apujm6Mdh9Nn/J9zOZn+cG0Wli3TtRkPnsWhg0TdcklSZJekEvQH+jes3vMPjGb+QHziYyLpLpLde5H38fN0Y3SOUsb7T5bAkNTlD18+M5h6q6si86gA6DsIw1z9sCTm9EsBQYOzsjY74aQPTtAPaON+0MktlVUFIU/Lv6BnZUda9utpX2J9lhapON/frNnw5QpEBoKJUqIxKquXU09KkmS0pF0/A6YOq48vkKZBWXQ6rW0Kd6GUe6jqO5i/CXELYGhr9VQDo2M5Ruf8wDvDcJPYp+Q1S4rR24feRl8LQxQ9mIxfr26mLzt3QmaBIUKvQjyS413RCip7j27x8zjM1kauJTT/U/jktmF9e3X42jraBbVv5IlLEyUh1QU0SShaFFYvBgaNUpXtcglSTIPn3wAVlWVA8EHuBlxk74V+lIkWxF+qP0DHUt2pFDWQql2X89dV15LfgKI1erx3HXlrQHyZOhJPP082XH5T25+fpba+TywUa3RqlpUvQ2B9otZ5O/+slFOSoJ8cl15fAUvPy9WnVuFzqCjfYn2Lz8kZLFLp/WMT50SjRHWrwdfX/jsM5g1Cyw/+f99JElKgU/2HUSr17Lh0ga8/LwIvB9IgSwF6FOuDxoLDf+r+b9Uv39SaihvCQzlu798uBW3FcXyNtFcJ3OCBUOPG3j4dAPfnv6W+Gu+OFXyZVR7D8ZOcH/tRFFygnxKPI55TOn5pdFYaPii/BeMch9lNglqH8xggL/+EoHX1xcyZoQRI6BAAfF1GXwlSUqhT/JdZNf1XfTb1o+QqBCKORVjcYvFdC/TPU0LZ7yvhvKWwFCG+6zntmYsaLSgwpDjMOpWEf52+h8Vp3UisxPM+dGd/v3dsXrDCajUbpSgqio7r+/k8O3DTKo/CSd7J9a0XUNtt9pm08842eLjoU8fsLMTQbhvX8iUydSjkiTpI/LJBOC7UXfRGXS4ObrhnMmZglkLMq/ZPJoWbmqSmsJjGhV9Yx/dIfXyMuPYDGb89TdP1ZyAXtRpNsBe+/bMu74Wq2BLRo4VjXMyZ377PVKrUYJWr2XdxXVMPTqV8w/P45LJhbGfjcXR1jH9Vqx6/Bjmz4ddu0S5SDs72L8fihXjjZ9uJEmSUuijD8Bn7p9hmv80/rjwBx1KdOD3dr9TKkcpDvQ6YNJxJS4Bf7/Th5CYAHLZFaFU/gcM3r2CCN0zaoVAeMaxPM1qiWrQY9Bbc/nESDIUfcCFHc64ub3/Hm8L8ikptBEQFkC79e248/QOJbOXZGXrlXQu1RlrTTptm3ftmuhItGIFxMZCkybw5Alkzw6ljZfxLkmS9G8fbQDed3Mfk45MYt+tfThYOzCk8hCGVRtm6mG9JqfTHa6rY4m3iidSZ+DyVWgbBGMuZOJyzpZMzlgPi/310Wc5iVVMebLV1VOg+BXc3JK2f2usRgmPYx4T9iyMMjnLUChrIYo7FWdu07kmWz0wGj8/kVBlZQXdu8PIkVCypKlHJUnSJ+KjCsDxuviXZ07/vv43QY+DmFJ/Cv0r9je7nrEnQk+w6twqEvQJGDCgqDAiyJFpHhO5NKYXs4ZZcOWALZZZnpO9iCN2hR9gbx3NmEYfNitLSaOE4MhgpvlNY2ngUoo5FeNU/1M42jryd/e/k3U9k9PrYetWiIqC3r2halWYNAl69YJcuUw9OkmSPjEfRSWsJ7FPWBCwgNknZrO81XIaF2pMVHwUtpa2ZrU0mpi0NHX/BA7eP0bV51k4lzlOFKtQLNnYch/bFtRg0SKR79Pm80iCsgRyPzomTc/wXnx4kYlHJrLuwjosFAu6l+nO6OqjX9a8TneePxdLzNOni3KRVauCv7/suytJUqr7aCth3Yy4yYxjM1gauJQYbQwNCzYkq11WADLZmEfGqn+IP77BvlgoFqw+sZgLz26Q9yl4H1fo61aPCwO+Yu+dozw45kGXWu48fw6DB8P334OTkyNQJ03GqaoqBtWAxkLD8dDj/HnlT4ZXG87wasPJmylvmowhVaxdK3ruPnkiAu/kydCmjQy+kiSZXLoNwAbVQN2VdQl7FkbX0l0Z6T6SMjnLmHpYr9l7cy8t17YkQZ+AhQouT/SsOmZL58++xGrtSFQXV8J8YPmYWty6Bc2bg6enSLxNK3qDni2XtzDVbypdSnVheLXhdC/TnTbF2qTfwhlBQeLcbt684OICtWq96kgkA68kSWYi3QZgC8WCVW1WUShrIfJkzGPq4bzmfvR9Zh/1xvvYDOLRoaKCouGLXE3osXstZM5MYCAM7wGHDolk2z17oH79tBtjvC6eVWdX4eXvxdXwqxTIUoBcDmIf1FpjjbWd+SzdJ4mqih+mlxds3y5mvbNniySrzz4z9egkSZL+I0UBWFEUR2AJUArRdfZzVVX9jTCuJKmVr1ayX5vSRghvcjX8KtP2/8LKi7+TgB6PW+DnZoFOY4G1xpo63b/jQVxmxo+CZcsgWzbROrZvX9CkXQ0QALr6dMUnyIcKuSuwrv062hVvl6aFSIxq0yaxtBwQII4P/fQTDBxo6lFJkiS9U0pnwDOBv1VVba8oijWQLvrKGaNG8pbA0JdneF3sK/Fzk7ZM2VafQG0Ifc7ASLu6FB78Pf75rfC9fZDqzh4c+t2dRr9CXJw48fLtt+DomErf5L/ce3aPWcdnMazaMHI55GJM9TEMrDSQevnrpc/mCDExYP/in9vff4vM5oULoUcPUURDkiTJzCU7C1pRlEzAWaCAmsSLmEs/4BqT97+xQpSzox1Hx9V97+u3BIYyzOcP7mjGAToUrHDVT2GS7SPqhoWSc9h4KFIEECujmzfDmDFw8ya0aCEqGxYubOzv6s2uhV/D08+TlWdXojPoWNV6Fd3KdEubm6eGsDDRCGHBAlG1qmpVEXwdHGRHIkmSzE5qZUEXAB4ByxVFKQucAoapqvr8XzfvD/QHcHV1TcHtjCclNZIT9AmM3Tadx4b5YCm6/KBqeaqeYY6mF13mvgrg587B8OFw4ICo77B7NzRoYIzv4P30Bj3dfLqx/uJ6rDXW9CnXh9HVR6dqh6dUdeGC+OSyZo04z9uu3avazLJGsyRJ6VBKArAlUAH4SlXV44qizATGAd/980mqqi4CFoGYAafgfi+ldP82uTWSn8U+peTUfITwlALREJIZdBYKYIWtofTLAB4eLo4RLVgglpjnzIEvv0z9BjqqqnLm/hnK5y6PxkJDJptMjPtsHMOqDiOnQ87UvXlqiosTmczx8TBggPhUk9iVSJIkKZ1KSUi4C9xVVfX4iz9vRATgVGWM/duk1EhO3OO9HXMYJ43CtNbetC7vTO/bWch23YnTObuz2z4DcZoL2BpKY2MoTu6M9syZI4JvVBQMGiTygbJmNeIP4A10Bh2bLm1iytEpBN4P5OKgi5TIXoJFLRal7o1Ti1YLGzbAli3wxx9gaws+PlCmTOr/MCVJktJIsgOwqqr3FUUJURSlqKqqV4B6wCXjDe3NjNHj9n01krcEhjLIZw73NVNRLQ1EAd+uzguM5Oe5l9hy+QkLfM5jo9VjoxPVoQx3sxN6vDxfXYe6dWHmTChVymjf9hvF6eJYHrgcL38vbkbcpGi2oixpsYSCWdJpD96oKFiyBGbMgJAQcSA6LEyc5/XwMPXoJEmSjCqli6JfAWteZEDfBPqkfEjvZqwet2+rkXzz8TXG/t6SexkuiwcUUAzwRHP3RZCv+1oAvx0MsUdK8eRiDtzcxImY1C60pKoqiqKwPuAaQ3aNwNKQj6LWPzGxeh/aVnBJvRunprNnxTJzVJQItvPni85EMrFKkqSPVIoCsKqqZ4A3ZnelltTocWtQDYTHhJM9Q3bsHkXyTL1Mh4u2bC6uRWehgmKJaln5tSDfsKgzpzc7M3WFOMP7yy8wapRYLU0tYc/CmO4/nbMPzjKw5DIm7Qgll24Olmou4uIVxm++iIVikSb1oo3i3Dm4fVukhpcsKToS9ekDldL0n5QkSZJJpLvpxZhGRbGzer1gRHJ73B4MPkj7lc0o8GMW2v8s1otzF69M05gZBORfh5NuMo667uRM+BUbQ3HyONqhqrBxIxQvDhMmQNu2cOUKjB+fesH3yuMr9P2zL/ln5sf7mDdO9k5M2XWOWK0eKzU3CmK6nbgUb9ZUVZT9atQIypYVB6JVVWSozZ0rg68kSZ+MdFeK0hg9bqPio/jf+i+Ze+MPscSsQO9nzqgJCSjW1jTv156jPucxaItjYygOiCDfqVBJ6teH/ftFPtBvv4lV07cxRrWtv6//TdM1TbGxtKFv+b6Mqj6KAlkKkH/cjjc+/0OX4tPU/v0i4J49K9r/TZwo0sPTYyEQSZKkFEp3ARhS1uMWYMWSIcx9/MfLP1tYaLDp0gPF2vrl9eFVkM9h60DOaxUZOcmBTJmSdqwoudnaqqqy5+Ye4nXxtCjagtr5avOTx098WelLcmTI8fJ5qbEUnyqiokRWc7Zs4s9arajD2bUr2NiYdmySJEkmlO6WoJPjctg5+nrVYvXqrwH4vN0vLNW0w87KDo2iwVpjjYebx2uvaV3emcNf1+X7os0Inl+b7Wsd+OILuHpVtAt835ned2Vrv4nOoGPdhXVUXFSRRqsbMdVvKgB2VnZ8V/u714IvGHcpPlXcvSvKf7m4iA1ygDp14Px5sc8rg68kSZ+4dDkDTqpFB72Z4TuZIB5hq4XCt6OhOzjkcuXzbzdS/EWvXg83D9xd3F977alToqHOsWPg7g5//QUVKyb93h+Srb318lZG7R7FjYgbFM1WlKUtl9Kt9LvLRRpjKT5VnDsnOhKtXSv2djt2FPWZQSw1y+VmSZIk4CMOwC0nlGCbIQhUsFQVNpWZQNMO/3vtOe4u7v8JvE+eiISqhQshRw5YsULEjw89DfO+JeLIuEgUFDLbZsagGshmnw3PBp60KtYKCyVpN0vpUrzRqOqrwDptmiiaMWQIDBsGbm4mHZokSZK5+miWoOO1cSxf9w0Rj0MAyGKXBUUFFFA1FpzNbfHO2ZfBIGpAFCkCixfD0KEiu7lXr+QdRX3bEvEXtR0Zu2csrtNd8fb3BqB1sdYc++IYbYq3SXLwNQtaLaxeDeXLw+nT4rFJk0QRjenTZfCVJEl6h3Q7A/Z/sXxcKWd5Avf+xsyQDYTZatHevU3/Ub8zoJMXG1bVI0Gf8MY93n8KCBD7uidOiN7tc+eKLOeU+PcScZZM4Tjl2cXAfRvRGXR0LNmRNsXbAKS/doBRUeJTyowZYq+3eHHxGECePCYdmiRJUnqR7HaEyWGsdoT+If7UW1WPOF0s6otZboN79nxdvC/1+k5EyZDh5fPetscLomnC//4nYknOnODpCd26pc42ZZt1bdh5befLrkQFs6bTcpF6vWiEcOeOqFg1Zgw0biwrVkmSJL1BarUjNBnfYF8S9Am8iL0MdmrG7O/+fC0IiDO4sYRFlmG7YyxjGoW+nJUaDOIkzNix8PSpaK7z44/G62qnqir7bu1j6tGpzG06l8LZCjOt4TTmN5tPLodcxrlJWjp3TjRF+PVXUfZr0iSxVi+LZkiSJCVbugzAHm4eWGusxfKypTVdW47/T/B92xncfDgzaJDIbq5ZUyw3ly5tnHHpDXo2X97M5COTOXXvFLkccnEz4iaFsxWmQJZ01j5PVWHfPpHRvGsX2NuLDfGiRcUZXkmSJClF0mUAdndxZ1/PfW9dXn7TGdzn0QoDBxt4eBycnGDVKlF62FjLzXqDnvILy3P+4XkKZy3MouaL6FG2B7aWqVgcOrXcuiVqbJ45I9bmf/1V9OGVrQAlSZKMJl0GYHjzEaJE/zxrq6rw/JIzEQeKYYixYfAgURfC0THlY4iKj+LPK3/SvUx3NBYaepbtiZujG22KtUFjoXn/BczJs2eiykjFiuDsLCpXLVkiNsVTs8OEJEnSJyrdBuB3STyDm/DIgSd7ShEfkg3r3BEU63OeObMrp/j6D58/ZOaxmcw9OZen8U8pm7MspXOWZnT10UYYfRoLC4NZs2DBAnBwgOBgsLaGvXtNPTJJkqSP2kcZgL+qVYzBY+J4cswNC2sdWRudI3ulMH5ql7LN3vCYcH7w/YGlgUuJ18XTtnhbxtYYS+mcRtpETkvXrolkqtWrRWZz27Yio/l9NTYlSZIko/jo3m23bYNxX+XhyW3IUTEMm+oXccmjYUyj0smuGhWdEI2DtQO2lrb4BPnQtVRXvq7xNUWdzKTuclKpqiieYW0NN26IzOb+/WHECCiYTo9FSZIkpVMfTQC+c0dUPtyyRfR2P3QIatbMAyS/MMSRO0eYfGQyNyJucGHgBTJYZ+DG0BvYWZlZx6H30elg0yaR0VynDkydKvrxhoS86lIkSZIkpal0Xz1BqxVxpXhxcVpm8mRRFbFmzeRdz6Aa2H51O58t+4yay2tyPPQ43Up3Q2vQAqSv4Pv8OcyeLc7sdu4sqlUlnrlSFBl8JUmSTChdz4CPHhWnYy5cgBYtRC5RSssPb7+6nVZ/tMI1syuzm8zm8/KfY29lb5Txprlhw2DpUqheHby9oWVLWbFKkiTJTKTLUpQAEyeKrkUuLmKS16pV8q4To41heeByrDRW9K/YH51Bh0+QD22KtcFKY2WUsaaZq1dFN6IhQ8RM9/Jl0d6penVTj0ySJOmT9K5SlOl2OpRYhvjSpeQF34jYCH499CtuM9wYsnMIf137CwBLC0s6luyYvoLv0aPQujUUKwYrV77qTFSsmAy+kiRJZirdLkFXr5782LLizAq+2vkV0QnRNCnUhG8++4bPXD8z7gDTgqqKZKo9e0SVqm+/FbPfHDlMPTJJkiTpPdJtAP5Q18Kv4WDtQO6MuSmQpQAtirRgbI2xlM1V1tRD+zCxsbB1K3TqJBKp6tUTe7t9+sCLLlCSJEmS+Uu3e8BJdfreaSYfmczGSxsZUmUIs5rMStP7G83jxzBvHsyZA48egZ8fuL+5FKckSZJkHj66doRJcTD4IBOPTGT3jd1kssnE2BpjGVZtmKmH9eEiIuC770T/xNhYaNpUbH5Xq2bqkUmSJEkp8FEFYFVVUV60N1pxdgVn7p9hUr1JDKw0kMy2mU08ug/05InY17WzE+W9OnWC0aNFlRFJkiQp3fsolqC1ei2/n/+dKUensKrNKirlqcTjmMdksMqQvgpnGAywcyd4esLt26Jes6UlxMeDjY2pRydJkiR9oI/yGBKIM7yzjs+i4KyC9N7a+//t3X2MVNUZx/HvwwppsNSXLqUUoWA1mDYibPCVlxpbG9kSaTEaCVgSTYiKiommgia+Ro02rdgGRdtiaYut1hZr8CUarDYmSkTlzWCFrRi2UqC0QrVVkH36x7nUyXjvMuwOc+4Zf59kM8M9Z3afhzP3PnvP3L2Hvi19+e+esBRha//WdIrvhx/CAw+Ev92dPBk2boTZs8MtJEHFV0SkCSU7Bb23ay+j7h1Fx786GD9sPAsnL2TSMZP+PwWdlGeegQsvDAV48eJw28h+/WJHJSIiB1GyBbilTws3fP0GRhwxIr2/4e3shPnzobUV5s4NF1YtXx4WSkjxFwgRETlgyRZggAtOuCB2CAdm7dqwcsSDD4abaMyaFbb36QNnnBE3NhERaaikPwNOym23wahR8MgjcOml4XPee+6JHZWIiESS9Blwqe1bg7etDY49Fs48M1zlfMklWgZQRER6fwZsZi1m9pqZLatHQMmrXoN30aKw/cQTw72aVXxFRIT6TEHPAdbX4fuk7/bbYdgwuOIKGDwYli6FW2+NHZWIiJRQrwqwmR0FfBv4WX3CSdDbb4cLqiBc3TxhArzwwsdLBPbRx+wiIvJJva0O84HvA11FHcxslpmtNLOV27dv7+WPK5EVK+Ccc2DEiFBsIUw9P/oojBsXNTQRESm/HhdgM5sMbHP3V7rr5+73u/tYdx87cODAnv64cujqgmXLYOLEsBjCs8/CvHnhIivQ2a6IiNSsN1dBjwPONrN24DPA58zs1+4+oz6hldAHH4R1d/v3h7vugosuggEDYkclIiIJ6nEBdvd5wDwAMzsduLrpiu/OnXDfffD44+Fst3//8HjccdC3b+zoREQkYZozzdPZGdbcHToUrrkmFNsdO0Lb8cer+IqISK/V5UYc7v4c8Fw9vld0L78Mp50Wrmw+77ywBm9bW+yoRESkyehOWO7w/POwdWtY9L6tDa69NnzWO3x47OhERKRJfXqnoPfuDfdlPvnksArRLbeEYtzSAjfdpOIrIiIH1aezAD/1FIwcCeeeC+++CwsXhqlnLQUoIiIN8umZgt6xIyyQMGgQHHpoWIv3zjthypRw1isiItJAzX8G/NZbcPnl4R7NN98cto0fDy++CFOnqviKiEgUzXsG/Npr4Qz34YdDkZ0+PazDC5pqFhGR6JqrALt/XFwXLAg30LjqKpgzB4YMiRubiIhIheaYgt6zB5YsgTFjwiIJEJYB3Lw5nAWr+IqISMmkXYDfew/uvhuOOQZmzIDdu+H990PboEFw2GFx4xMRESmQ7hR0VxeMHg0dHWEN3gULoL1dKxKJiEgS0i3AffqEm2cMHw6nnho7GhERkQOSbgEGmDYtdgQiIiI9ovlaERGRCFSARUREIlABFhERiUAFWEREJAIVYBERkQhUgEVERCJQARYREYlABVhERCQCFWAREZEIVIBFREQiUAEWERGJQAVYREQkAhVgERGRCMzdG/fDzLYDb9fxW7YC/6jj94tJuZRPs+QByqWsmiWXZskD6p/Ll919YF5DQwtwvZnZSncfGzuOelAu5dMseYByKatmyaVZ8oDG5qIpaBERkQhUgEVERCJIvQDfHzuAOlIu5dMseYByKatmyaVZ8oAG5pL0Z8AiIiKpSv0MWEREJEkqwCIiIhEkUYDN7Cwz+4uZbTSzuTntZmY/ztrXmFlbjDj3x8yGmtmfzGy9mb1uZnNy+pxuZjvNbFX2dX2MWGthZpvMbG0W58qc9tKPi5mNrPi/XmVmu8zsyqo+pR0TM1tkZtvMbF3FtiPN7Bkz25A9HlHw2m73q0YryOUHZvZG9v5ZamaHF7y22/dioxXkcqOZ/a3ifdRe8NrSjEtBHg9V5LDJzFYVvLZsY5J7/I26v7h7qb+AFqADOBroB6wGvlrVpx14EjDgFGBF7LgLchkMtGXPBwBv5uRyOrAsdqw15rMJaO2mPYlxqYi3Bfg74Q/nkxgTYCLQBqyr2HYnMDd7Phe4oyDXbverkuTyLeCQ7Pkdeblkbd2+F0uSy43A1ft5XanGJS+PqvYfAtcnMia5x9+Y+0sKZ8AnARvd/a/uvhv4LTClqs8U4JcevAQcbmaDGx3o/rj7Fnd/NXv+b2A9MCRuVAdVEuNS4RtAh7vX825tB5W7/xn4Z9XmKcDi7Pli4Ds5L61lv2qovFzc/Wl3/yj750vAUQ0PrAcKxqUWpRqX7vIwMwPOA37T0KB6qJvjb7T9JYUCPATYXPHvTj5ZtGrpUypmNhwYA6zIaT7VzFab2ZNm9rXGRnZAHHjazF4xs1k57amNy/kUH0xSGROAQe6+BcJBB/hCTp/UxgbgQsKMSp79vRfL4rJsOn1RwVRnSuMyAdjq7hsK2ks7JlXH32j7SwoF2HK2Vf/tVC19SsPMPgv8HrjS3XdVNb9KmAI9AfgJ8GiDwzsQ49y9DZgEzDaziVXtyYyLmfUDzgZ+l9Oc0pjUKpmxATCz64CPgCUFXfb3XiyDe4GvAKOBLYTp22opjcs0uj/7LeWY7Of4W/iynG29HpcUCnAnMLTi30cB7/SgTymYWV/C4C9x9z9Ut7v7Lnd/L3v+BNDXzFobHGZN3P2d7HEbsJQwTVMpmXEhHCRedfet1Q0pjUlm676p/uxxW06fZMbGzGYCk4Hpnn0gV62G92J07r7V3fe6exfwU/JjTGJczOwQYCrwUFGfMo5JwfE32v6SQgF+GTjWzEZkZynnA49V9XkM+F521e0pwM59Uwplkn1m8nNgvbv/qKDPF7N+mNlJhDHa0bgoa2Nmh5rZgH3PCRfLrKvqlsS4ZAp/m09lTCo8BszMns8E/pjTp5b9KjozOwu4Bjjb3f9T0KeW92J0Vdc/fJf8GJMYF+CbwBvu3pnXWMYx6eb4G29/iX1lWi1fhKtp3yRchXZdtu1i4OLsuQELsva1wNjYMRfkMZ4wbbEGWJV9tVflchnwOuEqu5eA02LHXZDL0VmMq7N4Ux6X/oSCeljFtiTGhPBLwxZgD+G39IuAzwPLgQ3Z45FZ3y8BT1S89hP7VQlz2Uj47G3f/rKwOpei92IJc/lVth+sIRy8B5d9XPLyyLb/Yt/+UdG37GNSdPyNtr/oVpQiIiIRpDAFLSIi0nRUgEVERCJQARYREYlABVhERCQCFWAREZEIVIBFREQiUAEWERGJ4H+nFKVGH0Gy8AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "prstd, iv_l, iv_u = wls_prediction_std(res_wls)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8,6))\n",
    "ax.plot(x, y, 'o', label=\"Data\")\n",
    "ax.plot(x, y_true, 'b-', label=\"True\")\n",
    "# OLS\n",
    "ax.plot(x, res_ols.fittedvalues, 'r--')\n",
    "ax.plot(x, iv_u_ols, 'r--', label=\"OLS\")\n",
    "ax.plot(x, iv_l_ols, 'r--')\n",
    "# WLS\n",
    "ax.plot(x, res_wls.fittedvalues, 'g--.')\n",
    "ax.plot(x, iv_u, 'g--', label=\"WLS\")\n",
    "ax.plot(x, iv_l, 'g--')\n",
    "ax.legend(loc=\"best\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Feasible Weighted Least Squares (2-stage FWLS)\n",
    "\n",
    "Like `w`, `w_est` is proportional to the standard deviation, and so must be squared."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:51:12.014507Z",
     "iopub.status.busy": "2021-02-02T06:51:12.013704Z",
     "iopub.status.idle": "2021-02-02T06:51:12.030025Z",
     "shell.execute_reply": "2021-02-02T06:51:12.031054Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            WLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.931\n",
      "Model:                            WLS   Adj. R-squared:                  0.929\n",
      "Method:                 Least Squares   F-statistic:                     646.7\n",
      "Date:                Tue, 02 Feb 2021   Prob (F-statistic):           1.66e-29\n",
      "Time:                        06:51:12   Log-Likelihood:                -50.716\n",
      "No. Observations:                  50   AIC:                             105.4\n",
      "Df Residuals:                      48   BIC:                             109.3\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          5.2363      0.135     38.720      0.000       4.964       5.508\n",
      "x1             0.4492      0.018     25.431      0.000       0.414       0.485\n",
      "==============================================================================\n",
      "Omnibus:                        0.247   Durbin-Watson:                   2.343\n",
      "Prob(Omnibus):                  0.884   Jarque-Bera (JB):                0.179\n",
      "Skew:                          -0.136   Prob(JB):                        0.915\n",
      "Kurtosis:                       2.893   Cond. No.                         14.3\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "resid1 = res_ols.resid[w==1.]\n",
    "var1 = resid1.var(ddof=int(res_ols.df_model)+1)\n",
    "resid2 = res_ols.resid[w!=1.]\n",
    "var2 = resid2.var(ddof=int(res_ols.df_model)+1)\n",
    "w_est = w.copy()\n",
    "w_est[w!=1.] = np.sqrt(var2) / np.sqrt(var1)\n",
    "res_fwls = sm.WLS(y, X, 1./((w_est ** 2))).fit()\n",
    "print(res_fwls.summary())"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}