{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/smart-stats/ds4bio_book/blob/main/book/linearSeparable.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "rIzXXVGuR9YM"
   },
   "source": [
    "# Linear separable models\n",
    "\n",
    "We've now covered two ways to do prediction with a single variable, classification using logistic regression and prediction using a line and least squares. What if we have several predictiors? \n",
    "\n",
    "In both the logistic and linear regression models, we had a linear predictor, specifically, \n",
    "\n",
    "$$\n",
    "\\eta_i = \\beta_0 + \\beta_1 x_i.\n",
    "$$\n",
    "\n",
    "In the continuous case, we were modeling the expected value of the outcomes as linear. In the binary case, we were assuming that the naturual logarithm of the odds of a 1 outcome was linear. \n",
    "\n",
    "To estimate the unknown parameters, $\\beta_0$ and $\\beta_1$ we minimized\n",
    "\n",
    "$$\n",
    "\\sum_{i=1}^n || y_i - \\eta_i||^2 \n",
    "$$\n",
    "\n",
    "in the linear case and \n",
    "\n",
    "$$\n",
    "-\\sum_{i=1}^n \\left[\n",
    "  Y_i \\eta_i + \\log\\left\\{\\frac{1}{1 + e^{\\eta_i}} \\right\\} \\right].\n",
    "$$\n",
    "\n",
    "in the binary outcome case (where, recall, $\\eta_i$ depends on the parameters).  We can easily extend these models to multiple predictors by assuming that the impact of the multiple predictors is linear and separable. That is,\n",
    "\n",
    "$$\n",
    "\\eta_i = \\beta_0 + \\beta_1 x_{1i} + \\beta_2 x_{2i} + \\ldots \\beta_{p-1} x_{p-1,i}\n",
    "$$\n",
    "\n",
    "If we think about this as vectors and matrices, we obtain\n",
    "\n",
    "$$\n",
    "\\eta = X \\beta\n",
    "$$\n",
    "\n",
    "where $\\eta$ is an $n \\times 1$ vector, $X$ is an $n \\times p$ matrix with $i,j$ entry $x_{ij}$ and $\\beta$ is a $p\\times 1$ vector with entries $\\beta_j$. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ApSoCaMXb4x-"
   },
   "source": [
    "Let's look at the voxel-level data that we've been working with. First let's load the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "GMMLqAkYRxb5"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import sklearn.linear_model as lm\n",
    "import sklearn as skl\n",
    "import statsmodels.formula.api as smf\n",
    "import statsmodels as sm\n",
    "\n",
    "## this sets some style parameters  \n",
    "sns.set()\n",
    "\n",
    "## Read in the data \n",
    "dat = pd.read_csv(\"https://raw.githubusercontent.com/bcaffo/ds4bme_intro/master/data/oasis.csv\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "HUih09BMdi9_"
   },
   "source": [
    "Let's first try to fit the proton density data from the other imaging data. I'm going to use the `statsmodels` version of linear models since it has a nice format for dataframes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "Onw6CyaCdrtz"
   },
   "outputs": [],
   "source": [
    "trainFraction = .75\n",
    "\n",
    "sample = np.random.uniform(size = 100) < trainFraction\n",
    "trainingDat = dat[sample]\n",
    "testingDat = dat[~sample]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 507
    },
    "colab_type": "code",
    "id": "fDJOnSsxe2N5",
    "outputId": "140b3f55-7e65-4bcf-91c3-417127a8769d"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                 Results: Ordinary least squares\n",
      "=================================================================\n",
      "Model:              OLS              Adj. R-squared:     0.780   \n",
      "Dependent Variable: PD               AIC:                65.4823 \n",
      "Date:               2024-01-29 06:47 BIC:                84.0222 \n",
      "No. Observations:   75               Log-Likelihood:     -24.741 \n",
      "Df Model:           7                F-statistic:        38.48   \n",
      "Df Residuals:       67               Prob (F-statistic): 4.34e-21\n",
      "R-squared:          0.801            Scale:              0.12678 \n",
      "------------------------------------------------------------------\n",
      "               Coef.   Std.Err.     t     P>|t|    [0.025   0.975]\n",
      "------------------------------------------------------------------\n",
      "Intercept      0.1001    0.1278   0.7832  0.4362  -0.1549   0.3551\n",
      "FLAIR          0.0537    0.0774   0.6934  0.4905  -0.1008   0.2081\n",
      "T1            -0.3114    0.0842  -3.6993  0.0004  -0.4794  -0.1434\n",
      "T2             0.5366    0.0830   6.4681  0.0000   0.3710   0.7022\n",
      "FLAIR_10       0.0698    0.3321   0.2102  0.8341  -0.5931   0.7327\n",
      "T1_10          0.3136    0.1531   2.0480  0.0445   0.0080   0.6193\n",
      "T2_10          0.2117    0.2849   0.7431  0.4600  -0.3569   0.7803\n",
      "FLAIR_20       1.2552    0.7459   1.6827  0.0971  -0.2337   2.7441\n",
      "-----------------------------------------------------------------\n",
      "Omnibus:              1.233        Durbin-Watson:           1.844\n",
      "Prob(Omnibus):        0.540        Jarque-Bera (JB):        0.763\n",
      "Skew:                 0.230        Prob(JB):                0.683\n",
      "Kurtosis:             3.179        Condition No.:           42   \n",
      "=================================================================\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the\n",
      "errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "results = smf.ols('PD ~ FLAIR + T1 + T2  + FLAIR_10 + T1_10 + T2_10 + FLAIR_20', data = trainingDat).fit()\n",
    "print(results.summary2())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 73
    },
    "colab_type": "code",
    "id": "7UrWqsZAfF1u",
    "outputId": "fce680da-2fbb-476f-a5b9-b013192885b8"
   },
   "outputs": [],
   "source": [
    "x = dat[['FLAIR','T1', 'T2', 'FLAIR_10', 'T1_10', 'T2_10', 'FLAIR_20']]\n",
    "y = dat[['GOLD_Lesions']]\n",
    "## Add the intercept column\n",
    "x = sm.tools.add_constant(x)\n",
    "\n",
    "xtraining = x[sample]\n",
    "xtesting = x[~sample]\n",
    "ytraining = y[sample]\n",
    "ytesting = y[~sample]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 71
    },
    "colab_type": "code",
    "id": "MaU1zUQluyFb",
    "outputId": "0bf756bd-d907-413a-a31e-4c92159e4d53"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.264729\n",
      "         Iterations 8\n"
     ]
    }
   ],
   "source": [
    "fit = sm.discrete.discrete_model.Logit(ytraining, xtraining).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 365
    },
    "colab_type": "code",
    "id": "clrxtIgQvDhJ",
    "outputId": "2fc144fe-1bd0-450b-e180-5788b1ca6769"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                           Logit Regression Results                           \n",
      "==============================================================================\n",
      "Dep. Variable:           GOLD_Lesions   No. Observations:                   75\n",
      "Model:                          Logit   Df Residuals:                       67\n",
      "Method:                           MLE   Df Model:                            7\n",
      "Date:                Mon, 29 Jan 2024   Pseudo R-squ.:                  0.6176\n",
      "Time:                        06:47:16   Log-Likelihood:                -19.855\n",
      "converged:                       True   LL-Null:                       -51.926\n",
      "Covariance Type:            nonrobust   LLR p-value:                 2.236e-11\n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const         -4.8293      1.876     -2.574      0.010      -8.507      -1.152\n",
      "FLAIR          2.2957      1.166      1.969      0.049       0.011       4.580\n",
      "T1             2.3580      1.055      2.236      0.025       0.291       4.425\n",
      "T2             2.1419      1.084      1.976      0.048       0.017       4.267\n",
      "FLAIR_10       7.7567      3.677      2.109      0.035       0.550      14.964\n",
      "T1_10          1.2147      1.572      0.773      0.440      -1.866       4.296\n",
      "T2_10         -4.9682      3.051     -1.628      0.103     -10.949       1.012\n",
      "FLAIR_20     -18.1028      7.771     -2.330      0.020     -33.334      -2.872\n",
      "==============================================================================\n"
     ]
    }
   ],
   "source": [
    "print(fit.summary())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "EkCEdWt0xGzn"
   },
   "source": [
    "Now let's evaluate our prediction. Here, we're not going to classify as 0 or 1, but rather estimate the prediction. Note, we then would need to pick a threshold to have a classifier. We could use .5 as our threshold. However, it's often the case that we don't necessarily want to threshold at specifically that level. A solution for evalution is to plot how the sensitivity and specificity change by the threshold. \n",
    "\n",
    "In other words, consider the triplets\n",
    "$$\n",
    "(t, sens(t), spec(t))\n",
    "$$\n",
    "where $t$ is the threshold, `sens(t)` is the sensitivity at threshold $t$, `spec(t)` is the specificity at threshold `t`. \n",
    "\n",
    "Necessarily, the sensitivity and specificity \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 301
    },
    "colab_type": "code",
    "id": "zxlcuTECxCyQ",
    "outputId": "94d8070e-71a3-42f1-a4f4-2a5354602f82"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkoAAAHKCAYAAAD8e2TcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB31UlEQVR4nO3dd3gUVdsG8Htr+qY3ILRAQu8ttAgoNUgv70tC6KEICqIoKk0EUUTEUEJRqoJUUYqASJNiAaULBKSEJKRv+rb5/uDNfq5JIG0zm839uy4uzdmZ2Xv3JNknc86ckQiCIICIiIiI8pGKHYCIiIjIUrFQIiIiIioECyUiIiKiQrBQIiIiIioECyUiIiKiQrBQIiIiIioECyUiIiKiQrBQIiIiIioECyUiIiKiQrBQogptz549CAwMNP5r0KABOnbsiOnTp+Pvv/8WOx4AoGvXrnjrrbfEjpFPVlYW1q5di/79+6N58+Zo1qwZ+vXrhzVr1iArK0vseEW2Zs0aHDt2LF/7hQsXEBgYiAsXLoiQ6qmHDx9iwYIF6NGjB5o0aYKmTZuiT58++PTTTxEfH2/cLiwsDCEhIaLlLI3vvvsOGzduNNvxS/Lzc/HiRXz++edQq9X5HgsLC0NYWFhZxaNKQC52AKKysHjxYtSuXRu5ubm4ePEi1qxZgwsXLuDQoUNwdnYWNVtkZCQcHR1FzfBviYmJGD16NB48eICwsDC88cYbAIDz589j9erVOHDgAL788kt4eHiInPT5oqKi0KNHD7z44osm7Q0bNsSOHTtQp04dUXL99NNPmDFjBlxdXTFixAg0aNAAAHDr1i3s3r0bJ0+exL59+0TJVpa+//573L59G6NGjTLL8Uvy83Pp0iVERkZiwIABUKlUJo/NnTu3LONRJcBCiaxC3bp10bhxYwBA27Ztodfr8fnnn+PYsWMYNGiQqNnyPiDLk16vh16vh1KpLPDxWbNm4e7du9i0aRNatWplbO/QoQOCg4MRHh6OWbNmYcOGDeUVGcDzcxeHo6MjmjVrVvpQJfDw4UPMmDEDNWvWxObNm+Hk5GR8LCgoCCNHjsTRo0fLNZMgCMjNzYWtrW25Pm9J5eTkwNbWtsx/fsQqnKni4tAbWaW8oikpKcmk/cqVK5g4cSLatGmDxo0bo3///jh48GC+/ePj4/Hee+8hODgYjRo1QseOHTFt2jQkJiYat8nIyMCSJUvQtWtXNGrUCJ06dcIHH3yQb9jqn0MHycnJaNSoEZYvX57vOaOjoxEYGIjNmzcb2xISEjBnzhx07twZjRo1QteuXREZGQmdTmfc5tGjRwgMDMS6deuwatUqdO3aFY0bN8b58+cLfG+uXLmCM2fOYNCgQSZFUp5WrVph0KBBOHPmDK5evWpsDwwMxIIFC7B9+3b06NEDjRo1Qu/evXHgwIF8xyht7tzcXHz44Yfo168fWrZsiTZt2mDYsGH5htgCAwORlZWFvXv3Godf84ZVChp6e+utt9C8eXPcv38f48ePR/PmzREcHIwPP/wQGo3G5NhxcXGYNm0amjdvjlatWuH111/H5cuXERgYiD179hT43ubZuHEjsrKyMHfuXJMiKY9EIkH37t3ztV++fBn//e9/0bRpU3Tr1g1r166FwWAwPl7U9yXvvVmwYAG+/vpr9OrVC40bN8bevXsBPD1LM2TIELRp0wYtWrTAgAEDsHPnThR0j/TvvvsOw4YNQ/PmzdG8eXP069cPO3fuBPB0GOvEiROIiYkxGQLPo9FosGrVKvTs2RONGjVCu3bt8PbbbyM5OdnkObp27YqIiAgcOXIE/fv3R+PGjREZGWl87J9DbwaDAatWrTIOZ7Zq1Qp9+/bFpk2bAACff/45PvroIwBAt27djJnyvg8KGnrTaDSIjIw0vk9t27ZFWFgYLl68mO/9oMqHZ5TIKj169AgAULNmTWPb+fPnMW7cODRt2hTz5s2Dk5MTDh48iOnTpyMnJwcDBw4E8LRIGjRoEHQ6HSZOnIjAwECkpKTgzJkzSEtLg4eHB7KzsxEaGoq4uDjjNrdv38aKFStw69YtbNy4ERKJJF8uNzc3vPDCC9i3bx+mTZsGqfT//1bZs2cPFAoF+vbtC+BpsTFkyBBIpVJMmTIF1atXx6VLl7B69WrExMRg8eLFJsfesmULatasiVmzZsHR0RE1atQo8L05e/YsAOQbqvqnbt26YceOHTh79iwaNWpkbD9+/DguXLiAadOmwc7ODl999RVmzJgBmUyGnj17lllujUaDtLQ0jBkzBt7e3tBqtTh79iymTp2KxYsXo3///gCAHTt2IDw8HG3btsXkyZMB4LnDNFqtFpMmTcLgwYMxZswY/Prrr1i1ahUcHR3xyiuvAHg6f2vkyJFIS0vDzJkzUaNGDZw+fRrTp09/5rHznDlzBh4eHsU6o5WQkIA33ngDo0ePxiuvvIKjR4/ik08+gZeXl/H1FvV9yXPs2DH89ttvmDJlCjw8PODu7g4AiImJwbBhw1ClShUAwB9//IGFCxciPj7e+B4AwGeffYZVq1ahe/fuGD16NJycnHD79m08fvwYwNNhrPfeew8PHz40FjZ5DAYDJk+ejN9//x1jx45FixYtEBMTg88//xyXL1/G7t27Tc5uXbt2DdHR0Zg0aRKqVasGOzu7At+n9evXIzIyEpMmTUKrVq2g0+lw9+5dpKenAwCGDBmCtLQ0bNmyBZGRkfD09ARQ+JkknU6HcePG4ffff8fIkSPRrl076PV6/Pnnn4iNjS1K15G1E4gqsN27dwsBAQHCH3/8IWi1WiEjI0M4deqU0KFDB2HEiBGCVqs1btuzZ0+hf//+Jm2CIAgRERFChw4dBL1eLwiCILz99ttCw4YNhTt37hT6vFFRUUK9evWEy5cvm7QfPnxYCAgIEE6cOGFs69KlizBr1izj1z/++KMQEBAgnDlzxtim0+mEjh07ClOnTjW2vffee0KzZs2EmJgYk+fYsGGDEBAQINy+fVsQBEF4+PChEBAQILz44ouCRqN57ns2Z84cISAgQIiOji50mzt37ggBAQHC3LlzjW0BAQFCkyZNhISEBJPcPXv2FF566SWz5tbpdIJWqxVmz54t9O/f3+SxZs2amby/ec6fPy8EBAQI58+fN7bNmjVLCAgIEA4ePGiy7fjx44UePXoYv966dasQEBAgnDx50mS79957TwgICBB27979zLyNGzcWhg4d+sxt/ik0NFQICAgQ/vzzT5P23r17C2PGjCl0v2e9LwEBAULLli2F1NTUZz63Xq8XtFqtEBkZKbRp00YwGAyCIAjCgwcPhPr16wuvv/76M/efMGGC0KVLl3zt33//vRAQECD88MMPJu2XL18WAgIChG3bthnbunTpItSvX1+4e/duvuP8++cnIiJC6Nev3zMzrV+/XggICBAePnyY77HQ0FAhNDTU+PXevXuFgIAA4ZtvvnnmMany4tAbWYWhQ4eiYcOGaNGiBcaNGweVSoVVq1ZBLn960vT+/fu4e/eu8WyNTqcz/uvcuTMSEhJw7949AMCpU6fQtm1b+Pv7F/p8P/30E+rWrYv69eubHKtjx46QSCT45ZdfCt23c+fO8PT0NBm+OXPmDJ48eWIyn+rEiRNo27YtvLy88uUFkO85unbtCoVCUcx37tn+fVYsKCjIZIK3TCZD7969cf/+fcTFxZVp7kOHDmH48OFo3rw5GjRogIYNG2LXrl2Ijo4u9Wvq2rWrSVtgYKDxLAkA/Prrr3BwcDBmzmPOK9M8PT3RpEmTZ+YCive+tGvXrsCLGc6dO4dRo0ahZcuWqF+/Pho2bIgVK1YgNTXVOFx99uxZ6PV6jBgxokSv56effoJKpUKXLl1Mvg/q168PT0/PfN8HgYGBqFWr1nOP27hxY9y8eRPz5s3D6dOnkZGRUaJ8eU6fPg0bGxvR5zKS5eLQG1mFJUuWwN/fH5mZmTh48CB27NiBGTNmYP369QBgnFu0ZMkSLFmypMBjpKSkGP/r7e39zOdLSkrC/fv30bBhw2ceqyByuRwvv/wytm7dCrVaDZVKhT179sDT0xMdO3Y0eY6ffvqpyM+RN8TwPHnDLY8ePULt2rUL3CZv6NLX19ekvaCr4PLaUlNT4ePjUya5jxw5gtdeew09e/bEuHHj4OHhAZlMhq+//hq7d+9+zit8Njs7O9jY2Ji0KZVK5ObmGr9OTU0t8LXmDV09j6+vr/E9LCoXF5d8bf/OVdz3paD39vLlyxg7dizatGmD999/Hz4+PlAoFDh27BjWrFmDnJwcADDOI/Lx8SnW68iTlJQEtVptMnT7TyX9/o2IiIC9vT3279+P7du3QyaToVWrVpg5c6ZxbmJxJCcnw8vLy2QYnOifWCiRVfD39zf+kmzXrh0MBgN27tyJw4cPo2fPnnB1dQXw9JfsSy+9VOAx8v6adXV1NVnjpiCurq6wsbHBokWLCn38WQYNGoQNGzbgwIED6N27N44fP47w8HDIZDKTYwQGBuK1114r8BheXl4mXxc0J6og7du3x7Jly3Ds2LF8Z0zy/Pjjj8Zt/+mfk9n/3Zb3QV8Wuffv349q1aph+fLlJo/nTdg1NxcXF1y+fDlfe0GvvyCdOnXCli1b8Mcff5TplXfFfV8Kem8PHDgAuVyOqKgok4Lx3xPC3dzcADyd1P7vgrkoXF1d4eLiYvxj5d8cHByem7Ugcrkco0ePxujRo6FWq3H27Fl8+umnGDduHE6cOFHo3KbCuLm54ffff4fBYGCxRAVioURW6Y033sCRI0ewYsUKdO/eHbVr10bNmjVx8+ZNzJgx45n7du7cGfv378fdu3cLPePywgsvICoqCi4uLvDz8yt2Pn9/fzRt2hR79uyBwWCARqMxTib/53OcPHkS1atXL9O1oBo3boyOHTti9+7dxqun/um3337D7t270alTp3xnA86dO4fExETj2Ra9Xo+DBw+ievXqxjMPZZFbIpFAoVCYfHgmJCQYC7h/UiqVxrMgZaV169Y4dOgQTp48ieDgYGN7QVf4FSQ8PBy7d+/G/Pnz8y0PADy9VP/YsWOFFu2FKc778qxjyGQyk6IgJycH+/fvN9muQ4cOxrNVzZs3L/R4hb3/L7zwAg4cOACDwYCmTZsWOV9xqFQq9OzZE/Hx8Vi0aBFiYmJQp04d4/IS/zwbV5hOnTrh+++/x549ezB48GCz5KSKjYUSWSVnZ2dMmDABH3/8Mb777jv069cP8+fPx/jx4zF27FgMGDAA3t7eSEtLQ3R0NK5du4YVK1YAAF599VWcOnUKoaGhiIiIQEBAANLT03H69GmMGjUK/v7+CA8Px5EjRxAaGopRo0YhMDAQBoMBsbGxOHPmDMaMGfPcD4dBgwZhzpw5ePLkCZo3b56vKJs2bRrOnj2L4cOHIywsDLVq1YJGo8GjR49w6tQpzJ8/v8TDIkuWLMHo0aMxduxYhIWFISgoCMDTKwM3b96M2rVr48MPP8y3n6urK8LDwzF58mTjVW93797Fp59+Wqa5X3jhBRw5cgTz5s1Djx49EBcXh1WrVsHLyyvfiusBAQH45ZdfcPz4cXh6esLBwaHQAreoBgwYgE2bNuHNN9/Eq6++iho1auDUqVM4c+YMADz3zIOfnx+WLVuG6dOno1+/fggNDUX9+vUBPF0GYvfu3RAEodiFUnHel8IEBwfjyy+/xOuvv45hw4YhNTUVGzZsyLd2VbVq1RAREYFVq1YhJycHISEhcHJywp07d5CSkoJp06YBePr+HzlyBF999RUaNWoEiUSCxo0bo0+fPvjuu+8wYcIEhIWFoUmTJlAoFIiLi8OFCxfQrVu3Yr9+AJg4cSLq1q2LRo0awc3NDTExMdi0aROqVq1qvNIzICAAwNMzbQMGDIBcLketWrUKvCIyJCQEe/bswbx583Dv3j20bdsWgiDgzz//hL+/P/r06VPsjGRdWCiR1QoLC8O2bduwatUqhISEoF27dti5cyfWrFmDRYsWQa1Ww8XFBf7+/ujVq5dxP29vb+zatQsrVqzAunXrkJqaCldXV7Rs2dI4vGRvb49t27Zh7dq12LFjBx49egRbW1v4+vqiffv2qFq16nPz9enTB4sWLUJcXJzJJdl5vLy8sGvXLqxatQobNmxAfHw8HBwcULVqVXTq1CnfisPF4eHhgR07dmDLli04dOgQtmzZAgCoXr06IiIiEB4eDnt7+3z7de3aFXXq1MHy5csRGxsLPz8/LF26FL179y7T3IMGDUJSUhK2b9+O3bt3w8/PDxMmTEBcXFy+y9DfeecdzJ8/HzNmzEB2djbatGljfD0lZW9vj02bNmHRokX4+OOPIZFI0LFjR8ydOxcTJkwocG2kf+vSpQu+++47fPHFF9i+fTtiY2MhlUpRrVo1dOrUCaGhocXOVZz3pTBBQUFYtGgR1q1bh4kTJ8Lb2xtDhw6Fm5sb3nnnHZNt84rErVu3YubMmZDJZKhZs6bJOkQjR47E7du38emnnyI9PR2CIOCvv/6CTCbD6tWrsXnzZnz77bdYu3YtZDIZfHx80Lp1a2MxU1xt27bFDz/8gJ07dyIjIwOenp5o3749Jk+ebLwooG3btoiIiMDevXuxc+dOGAwGbN68GW3bts13PLlcjnXr1iEqKgoHDhzApk2b4ODggHr16qFTp04lykjWRSIIBawwRkT0L4GBgRgxYgTmzJkjdhTRrFmzBsuXL8eJEydKfDaPiCoWnlEiIirA1q1bAQC1a9eGVqvF+fPnsWXLFrz88ssskogqEdELpfv372PDhg34888/cfv2bdSuXRvff/99kfbdu3cvoqKiEBMTgxo1amDKlCkmQyhERCVla2uLTZs24dGjR9BqtfD19cX48eMxadIksaMRUTkSvVC6ffs2Tp48iaZNm8JgMBR4r6GCHD58GG+99RYmTJiADh064NixY5g+fTqcnJxM1qIhorLx119/iR2hXA0ePJhXQRGR+HOU/rl2xVtvvYWrV68W6YxSr169EBAQgM8++8zYNnbsWKSnp+Obb74xW14iIiKqPERfXaskC3w9fPgQd+/ezXc7gZCQEFy+fDnfnamJiIiISkL0Qqkk7t69CwD51krx9/eHIAjGx4mIiIhKo0IWSmlpaQCQbz2WvFWA8x4nIiIiKg3RJ3OXxr/vDZQ33aqo9wwqiCAIpdqfns0gCHiSnCV2DCIisiIeLnaQy8xz7qdCFkr/PHP0zzt8q9VqAPnPNBWHRCKBWp0Nvd5QupBUoFyNHuM/+gkAEDm9M2wUsgK3k8okUDnZQZ2eDYOea6KKjf1hOdgXloN9IZ6c6NtI2LYFPpOnQeHhAXOe3qiQhVLe3KS7d+/C39/f2B4dHQ2JRFLq+zzp9QbodCyUzOGf76tMIoFMWvC3t1wmha2NHNlZUugE9oXY2B+Wg31hOdgX5U8wGJDyw2Ek7t0FGAxI2bcLvhPMu7ZZhZyj5Ofnh9q1a+PgwYMm7d9//z2aNGkCNzc3kZIRERGROejS1YhZsRyJu78BDAY4tWkH75GjzP68op9Rys7OxsmTJwEAMTExyMjIwOHDhwEAbdq0gZubG2bPno19+/bh+vXrxv2mTZuG6dOno3r16mjfvj1+/PFH/Pzzz1i/fr0or4OIiIjMI+vWX4hduxr61FRIFAp4/mcEnDsFl8ucYtELpaSkJLz66qsmbXlf593t2WAwQK/Xm2zTq1cv5OTkYM2aNdiwYQNq1KiBTz/9lKtyExERWZHMq1cQ89kyQBCg9PGF78TJsKnmV27PL/rK3JYoJSWTc5TMJFejx6RlT88grp4RDBtlwZO55XIpXF0d2BcWgv1hOdgXloN9UT4MWg0eLv4ANlWrwWtEGKS2tvm2cXNzgIxXvREREVFlkHPvLmyq14BEJoNUoYTfm29BamsnSpYKOZmbiIiIrI9gMCDx2714sOh9JH2/39guVpEE8IwSERERWQBdaipi10ch++YNAIA+Lc0iFoFmoURERESiyrx2FXHro6BPT4fExhbeI8OhahskdiwALJSIiIhIJIJej6Rv9yL50AFAEGDj5wffiClQ+viIHc2IhRIRERGJQpuQgJRjRwBBgHNwF3gO/w+kCqXYsUywUCIiIiJRKH184B0WDolcAafWbcSOUyAWSkRERFQuBJ0Oifv2wLFZc9jVqQsAUAV1EDnVs3F5ACIiIjI7bVISHn78IVIOH0Ts2tUwaDRiRyoSnlEiIiIis8r44xLivlgPQ1YmpHZ28Bz2X0iVljUXqTAslIiIiMgsBJ0OCbu+QeqxIwAA21q14TthEhSeniInKzoWSkRERFTm9FmZeLRsKXL/vgcAcH2pBzwGDYFEXrFKj4qVloiIiCoEqZ09FK5u0D55Ap8x4+DYrLnYkUqEhRIRERGVCYNWA+gNkNraQiKRwHvUGBhycqBwdxc7WonxqjciIiIqNU18HB4u/gDxm76AIAgAAJmDQ4UukgCeUSIiIqJSUv9yHk82b4QhJwe65GToUpKhcKvYBVIeFkpERERUIgaNBgnbv0LaqRMAALu6AfCZMAkKV1dxg5UhFkr/YhAE5Gr00OkMYkexSrlavdgRiIioDGjiYvF4zSpoHj0EJBK49QmBe9/+kMhkYkcrUyyU/uVJchbGf/ST2DGIiIgslmAwIObz5dDGx0PmpILPuAlwaNhI7FhmwUKJRFGnmjOUCl5LQERUEUmkUniHjULywe/hM2Y85C4uYkcyGxZKBYic3hkyiUTsGFZNqZBCwveYiKjCyI2JgTbhiXE9JPt69WEXWM/qf5ezUCqAjUIGmdS6O56IiKgoBEGA+ufTePLVVkAiQY335kHp4wsAVl8kASyUiIiIqBCGnBzEb92E9PPnAAD2DRtBau8gcqryxUKJiIiI8sl9+BCPo1ZCGxcHSKXw6D8Qrj17QyKtXPNLWSgRERGRidRTJ5Dw9TYIWi3krq7wnTAJdnUDxI4lChZKREREZEL75AkErRYOjZvAZ8x4yJycxI4kGhZKREREBMFgMA6refQfCJsqVeHULqjSDbX9W+V+9URERJWcIAhIOX4Mj5YugaDTAQAkcjlU7TtU+iIJ4BklIiKiSkuflYn4jV8g4+LvAAD1+bNw7thZ5FSWhYUSERFRJZRz7y5io1ZDm5gAyGTwHDIMqg6dxI5lcVgoERERVSKCICD12BEk7PoG0Ouh8PCEb8Qk2NaqLXY0i8RCiYiIqBJJ3L0TKYcPAgAcW7aCd/hoyCrZIpLFwUKJiIioEnHu1BnqM6fh/nI/OHfpViluQ1IaLJSIiIismGAwICc6GnZ16wIAlN4+qPXhx5Da2oqcrGLgdX9ERERWSp+ejpgVy/Hwo0XIunHd2M4iqeh4RomIiMgKZd36C3Hr1kCXkgKJQgFdWqrYkSokFkpERERWRDAYkHzweyR9uxcQBCh8fFAlYgps/PzEjlYhsVAiIiKyErq0NMRtWIus69cAAE5B7eE9YiSH2kqBhRIREZGVyLx6BVnXr0GiVMJrRBicuYBkqbFQIiIishKq9h2gTYiHU5t2sKlSVew4VoFXvREREVVQutRUxK6Pgj4jAwAgkUjg0X8Qi6QyxDNKREREFVDmtauIWx8FfXo6YBDgO2Gi2JGsEgslIiKiCkTQ65H07V4kHzoACAKU1fzg3vdlsWNZLRZKREREFYQ2ORlx69Yg+/YtAIBz8AvwHPZfSJVKkZNZLxZKREREFUD23buIWbEMhowMSG1t4TVyFFRt2okdy+qxUCIiIqoAlN7ekNrYQOHmDt+IyVB6e4sdqVJgoURERGShdOlqyBydIJFIIHNwQLUZb0Du5gapgkNt5YXLAxAREVmgjD8u4e933ob69Cljm9Lbh0VSOWOhREREZEEEnQ5Ptn+Fx5GfwZCVCfX5sxAEQexYlRaH3oiIiCyENiEBj6NWIffvewAAl5d6wHPQEEgkEpGTVV4slIiIiCxA+u+/In7jFzBkZ0Nq7wCfMePg2Ky52LEqPRZKREREItPExyF2zSpAEGDrXwe+EyZC4e4hdiwCCyUiIiLRKb194P5yfxhyc+HRfyAkcn48Wwr2BBERkQjSf7kAm+o1oPTxAQC49+0nciIqCK96IyIiKkcGjQbxWzYidu1qxEathEGrETsSPQPPKBEREZUTTVwsHq9ZBc2jh4BEAocmzSCRysSORc/AQomIiKgcqM+dRfzWTRBycyFzUsFn3AQ4NGwkdix6DhZKREREZmTQaPDkqy1QnzkNALCrVx++4yIgd3ERNxgVCQslIiIiM5JIpdA8fgxIJHDv2w9uIS9DIuUU4YqChRIREVEZEwQBEARIpFJI5HL4RkyCNiEB9vXqix2NiomFEhERURky5OQgfusmyJ1d4DlkGABA4e7BBSQrKIsolO7du4eFCxfi999/h52dHfr06YOZM2fC1tb2mftlZWVh1apVOHz4MBISEuDt7Y2+ffsiIiICSiXvrkxEROUr9+FDPI5aCW1cHCCTwfmFLlB6eokdi0pB9EJJrVYjPDwcVapUwYoVK5CcnIzFixcjNTUVS5cufea+8+bNw7FjxzB9+nTUrVsXly9fxooVK5CWloZ33323nF4BERFVdoIgIO3USSR8vRWCTge5qyt8xk9kkWQFRC+Utm/fDrVajX379sHNzQ0AIJPJMHPmTEyaNAn+/v4F7qfT6XD48GGMGzcOYWFhAIB27drh8ePHOHjwIAslIiIqF/rsbMR9+QXSf7kAALBv1AS+Y8dD5uQkcjIqC6JPuz916hSCgoKMRRIA9OjRA0qlEidPnix0P0EQoNfr4fSvb0SVSvV0Eh0REZGZCQYD7n+46GmRJJXCY/BQVJ32GoskKyJ6oRQdHZ3vrJFSqUT16tURHR1d6H4KhQIDBw7Eli1b8OeffyIzMxPnz5/HN998gxEjRpg7NhERESRSKdx79oLczQ1+s2bDrWdvXvpvZUQfelOr1VCpVPnaVSoV0tLSnrnvvHnzMHfuXAwdOtTYFhYWhldeeaVUmaQyCeQyfqOLSfa/91/GfrAI7A/Lwb4Qnz4zE9rERNjUrgUAcOvYEU4tWkJqYyNysspLIjHfsUUvlAojCAIkz3nlS5cuxYkTJ/D++++jVq1auHbtGlasWAGVSoVp06aV+LlVTnawtbHYt6ZSUansxI5A/8D+sBzsC3Gk376D6I8/gUGjQbPlywDYPe0L9ofVEr0aUKlUUKvV+drT09MLncgNALdu3cIXX3yBVatWoVu3bgCA1q1bQyKR4KOPPsKIESPg7u5eokzq9GxkZ/GvNTHJZFKoVHZQq7Oh1xvEjlPpsT8sB/tCHIIgIPnID4j/Zgeg10Ph6YnUmHh4uTizLyyAs7MdpGYa8hS9UPL39883F0mj0eDBgwcYNGhQofvduXMHAFC/vukqp/Xr14dOp0NMTEyJCyWDXoBO4De9JdDrDdDp2BeWgv1hOdgX5UefkYG4jRuQ+cclAIBjy1bwDh8NherphG32hfjMeQ2X6KdNOnfujPPnzyMlJcXYdvToUWg0GgQHBxe6X9WqVQEA165dM2m/evUqAKBatWpmSEtERJVJdvQd3F8wF5l/XIJELofXf0PhO3EKZPYOYkejciL6GaXhw4dj69atmDx5MiZPnoykpCR8+OGH6Nu3r8nQ2+zZs7Fv3z5cv34dANCoUSM0adIEc+fORWJiImrVqoUrV65g1apV6N27t8lyA0RERCWR+tOP0CUnQeHpBd+Jk2Fbo6bYkaiciV4oqVQqbNq0CQsXLsTUqVNha2uLkJAQzJw502Q7g8EAvV5v/Fomk2HNmjX47LPPsG7dOiQmJsLX1xehoaGYOHFieb8MIiKyQl4jRkLupILby/0hs+OE7cpIInB1RhNxSZmATg+Z1IzXGtJzyeVSuLo6ICUlk2P/FoD9YTnYF+aVdesvpF84D6/Qkc+98pp9YTnc3BzMtmSG6GeUiIiIxCYYDEg+dABJ3+4FDAbY1qwJ506Fz5OlyoOFEhERVWo6tRpx66OQdf3pxUFO7YLg1LqtyKnIUrBQIiKiSivr5g3ErlsDfVoaJEolvP4bClWHTs8ddqPKg4USERFVSinHjyHh622AIEBZpQp8I6bA5n9LzxDlYaFERESVkm2NmoBUClW79vD6byjv1UYFYqFERESVhi4tDXJnZwCAnX8d1Jj7PmyqVBE5FVky0VfmJiIiMjdBr0fi3t24N/tN5MY8MrazSKLn4RklIiKyatqUFMStXY3s27cAABmXLsKmKm9zRUXDQomIiKxW5pXLiNuwDvqMdEhtbeE1chRUbdqJHYsqEBZKRERkdQSdDon79iDl8EEAgE31GvCNmASlt4/IyaiiYaFERERWJ+3nM8YiyblLN3gOHQapQilyKqqIWCgREZHVce7UGZlXL0PVNghOrVqLHYcqMF71RkREFZ6g0yH5h0MwaDUAAIlUiqpTprFIolLjGSUiIqrQtAkJiF27Gjn37kKbmAjvEWFiRyIrwkKJiIgqrPTff0P8xg0wZGdDam8PhwYNxY5EVoaFEhERVTgGrRaJO7cj9fiPAADb2v7wjZgEhbuHyMnI2rBQIiKiCkWT8ASxq1ci98F9AIBrj17wGDAIEjk/0qjs8buKiIgqFIlEAm1iAqSOjvAZMw6OTZqJHYmsGAslIiKyeILBAIn06YXaCg9PVJk8FQovbyjc3ERORtaOywMQEZFF08TF4sH785B55bKxzb5efRZJVC54RomIiCyW+txZxG/dBCE3Fwk7t8O+YSPjmSWi8sBCiYiILI4hNxdPvt4K9ZnTAAC7wHrwHR/BIonKHQslIiKyKLmPYxC7ZhU0j2MAiQRuIS/DvW8/FkkkChZKRERkMbSJCXiwcD4EjQYyZ2f4jouAff0GYseiSoyFEhERWQyFhyec2rSDLikJPuMmQO7sLHYkquRYKBERkahyHz2ETOUMuUoFAPAaEQaJTMahNrII/C4kIiJRCIKA1FMn8OCDBYjbsBaCwQAAkCoULJLIYvCMEhERlTtDTjbiN29C+i/nnzZIJBA0uZDY2okbjOhfWCgREVG5ynlwH7FrVkH7JB6QSuExYDBce/TkWSSySCyUiIioXAiCgLQTx5Gw42sIOh3kbm7wnTAJdnXqih2NqFAslIiIqFwIGg1Sjh6BoNPBoWkz+IweB5mjo9ixiJ6JhRIREZULqY0NfCdORvbNm3B5qTskEonYkYieq0QDwtHR0ZgxYwY6duyIRo0a4dq1awCAyMhInD9/vkwDEhFRxSQIAlKOHUHqT8eNbbbVa8C1ew8WSVRhFLtQunHjBgYPHoxffvkFbdq0gV6vNz6WmZmJ7du3l2lAIiKqePQZGXi8cgUStn+FhB1fQRMXJ3YkohIp9tDb0qVLERgYiC+//BIKhQIHDx40PtakSRMcOXKkTAMSEVHFkh19B7FRq6FLToJELofH0OFQeHuLHYuoRIpdKF28eBEff/wx7OzsTM4mAYCHhwcSExPLLBwREVUcgsGAlKM/IHHPLkCvh8LTC74Rk2Fbs6bY0YhKrESTuRUKRYHtaWlpUCqVpQpEREQVjyAIeLw6EpmXLgIAHFu1gffIUZDZ24ucjKh0ij1HKTAwEMeOHSvwsdOnT6Nhw4alDkVERBWLRCKBXW1/SORyeIWOhG/EJBZJZBWKfUZp5MiReP3112FnZ4d+/foBAGJjY3H+/Hns3r0bK1asKPOQRERkeQSDAfr0dMidnQEArj16wbF5Syh9fERORlR2il0o9e7dGw8ePEBkZCS2bNkCAJg6dSpkMhmmTZuGrl27lnlIIiKyLDq1GnEb1kKXkozq78yF1MYGEqmURRJZnRLNUZo4cSL69++P06dPIykpCa6urujYsSOqVq1a1vmIiMjCZN28gdh1UdCnpUKiVCLn/t+wDwgUOxaRWRS7UPr111/RoEED+Pj4YMiQISaPZWZm4vr162jdunWZBSQiIssgGAxIPvAdkvbvAwQBSt8q8J04BTb8I5msWLEnc48cORLR0dEFPnbv3j2MHDmy1KGIiMiy6NJS8WjZx0j6di8gCFB16ITq785lkURWr9hnlARBKPQxnU4HqbREd0UhIiIL9uSrrci+eQMSpRLeYeFQBXUQOxJRuShSoZSRkQG1Wm38OiEhAY8fPzbZJicnB3v37oWHh0fZJiQiItF5/WcE9JmZ8B4RBqVvFbHjEJWbIhVKGzduxMqVKwE8XSvjlVdeKXA7QRAQERFRdumIiEgU2pQUZF76HS5dXwQAyF1c4TdzlsipiMpfkQqlDh06wN7eHoIg4OOPP0ZoaCiqVDH9i0KpVCIgIABt2rQxS1AiIiofmVcvI279Ougz0iFzdoFTy1ZiRyISTZEKpebNm6N58+YAgOzsbAwZMgTevMEhEZFVEXQ6JO7bg5TDT292buNXHTZVq4mcikhcxZ7MXdiwGxERVVza5CTERq1GTvQdAIBzl67wHDocUgXv30mVW4kWnNTr9Th16hSio6ORk5Nj8phEIsGUKVPKJBwREZlfxuU/EbdhLQyZmZDa2cE7fAycWnE9PCKgBIVSSkoKRowYgbt370IikRiXC5BIJMZtWCgREVUcglYLQ2YmbGrWgm/EJCg9vcSORGQxil0offrpp7CxscFPP/2ELl264JtvvoGLiwu+/vprnDhxAhs3bjRDTCIiKkuCXg+JTAYATydrT3oFDk2aQqpQiJyMyLIUe3XI8+fPY9SoUfDyevoXh1QqRfXq1TFr1iy0b98eS5YsKfOQRERUdjIu/Y6/58yGNiXF2ObUshWLJKICFLtQiouLQ9WqVSGTySCVSpGdnW18rEuXLvj555/LNCAREZUNg1aLJ19vw+OVn0MbH4+UQwfEjkRk8YpdKLm6uiIjIwMA4OXlhVu3bhkfS0tLg16vL7t0RERUJjRPnuDhhx8g9cejAADXHj3hOXS4yKmILF+x5yg1bNgQt2/fxgsvvIDOnTtj1apVcHR0hEKhwLJly9C0aVNz5CQiohJK//UXxG/6AoacHEgdHOAzdjwcmzQTOxZRhVDsQik0NBQPHjwAALz22mv4888/MWvW02Xtq1evjnfeeadsExIRUYmpz59F3Pq1AADbOnXhO2EiFG7uIqciqjiKXSi1b98e7du3BwC4ublh3759uHXrFiQSCWrXrg25vERLMxERkRk4Nm8JZdVqcGzaDO79BhivdCOioin2HKV/k0gkCAwMREBAAGQyGb799tuyyEVERCWUef0aBIMBACC1sUH1d+fAY+BgFklEJVDqQinPwYMH0adPH7z11ltldUgiIioGQ24u4jZ9gZhlHyPlh0PGdt6GhKjkijxOtnbtWmzfvh1JSUmoVasWZsyYgc6dO+PixYuYP38+bt26BXd3d7z33nvmzEtERAXIffwYsVGroIl5BEgkEHQ6sSMRWYUiFUrbtm3DsmXL4OTkhICAAMTGxmLKlCl47733sGDBAsjlckyZMgVjxoyBvb19sUPcu3cPCxcuxO+//w47Ozv06dMHM2fOhK2t7XP3TU1NxfLly3Hs2DGkpaWhSpUqGD16NIYP52WvRFQ5pP18Bk+2bYag0UCmUsF3/ETY128gdiwiq1CkQmn37t1o2bIloqKi4OjoCL1ej3nz5mHu3LmoWrUqNmzYgBo1apQogFqtRnh4OKpUqYIVK1YgOTkZixcvRmpqKpYuXfrMfTMzMxEWFgYbGxvMnj0b7u7uuH//PrRabYmyEBFVJIbcXDzZthnqs08X+rWv3wA+4yZA7uwibjAiK1KkQunevXv4+OOP4ejoCACQyWSYNGkSdu7ciVdffbXERRIAbN++HWq1Gvv27YObm5vx+DNnzsSkSZPg7+9f6L5RUVHIycnBzp07jWef2rZtW+IsREQViSY+Dum/XAAkErj3GwC33iGQSMts6ikRoYiTubOzs433dsvj7e0NAKUqkgDg1KlTCAoKMhZJANCjRw8olUqcPHnymfvu3r0bgwcPLtIQHRGRtbGtXgNeYeGoNnMW3ENeZpFEZAalXvSotOsmRUdHY9CgQSZtSqUS1atXR3R0dKH7PXz4EImJiVCpVIiIiMDPP/8MBwcH9O7dG7NmzSpV8SSVSSCX8ReOmGT/e/9l7AeLwP6wDPrsbDzeuBWKgS9D5ukLAHAPDhY5VeXFnwvLIZGY79hFrnKWLFkCJyenfO2LFi0yDskBT9dVWr16dZEDqNVqqFSqfO0qlQppaWmF7peYmAgA+Oijj9CzZ0+sW7cOd+7cwbJly6DVarFw4cIiZ8j33E52sLXhwpmWQKWyEzsC/QP7QzwZd+/h7sefIOdxLHL/vovmKz7lukgWgj8X1q1I1UCVKlUQGxuL2NjYfO2PHz82aZOUUVknCMIzj2X432Jq/v7+WLx4MQAgKCgIOp0OH330EV599VV4enqW6LnV6dnIzuJfCGKSyaRQqeygVmdDrzeIHafSY3+IRxAEpPx0HPFffQVBp4XCzQ11XpmM9EwN+0Jk/LmwHM7OdpCaaei5SIXS8ePHzfLkwNMzR2q1Ol97enr6Mydyu7i4AADatWtn0t6uXTsYDAZER0eXuFAy6AXoBH7TWwK93gCdjn1hKdgf5UuflYX4zV8i47dfAQAOTZqi6oQJUPn5ICUlk31hIfhzIT5BMN+xRR9f8vf3zzcXSaPR4MGDB/nmLv2Tn58fFApFvnbhf++WuSpLIqLyoE1JwaOPFkGbkADIZPAcNAQuL/WAXMHhNqLyJHo10blzZ5w/fx4pKSnGtqNHj0Kj0SD4GZMUlUolOnTogHPnzpm0nzt3DnK5HHXq1DFbZiIic5M7O0Ph5Q25uzv8Zs2Ga/eeZTa1gYiKTvRCafjw4XBycsLkyZNx+vRp7Nu3D++//z769u1rMvQ2e/ZsNGhgutLslClT8Ndff+HNN9/EmTNnsHHjRnz++ecYMWKEyXIDREQVgT4zE4bcXACARCqF77gI1JizAHa1C5+GQETmJfrQm0qlwqZNm7Bw4UJMnToVtra2CAkJwcyZM022MxgM0Ov1Jm1NmjRBVFQUPvnkE0ycOBEuLi4IDQ3Fq6++Wp4vgYio1LLvRiM2ahXs6zeEz6gxAABZAVcaE1H5kgiCOadAVTxxSZmATg+ZlKe4xSSXS+Hq6sAJqxaC/WE+giAg5chhJO7ZBej1UHh6ovq78yBzcChwe/aF5WBfWA43NwezrWcl+hklIqLKSp+Rgbgv1iHz8p8AAMdWbeA9chRkJbi5OBGZR4kLpfT0dPzxxx9ISUlBcHAwnJ2dyzIXEZFVy759G7FrV0OXkgyJXA7P4f+Fc3AXTtgmsjAlKpRWrlyJdevWIScnBxKJBLt27YKzszPCw8PRoUMHTJgwoaxzEhFZDYNGg8drIqFPS4PC2xu+EZNhW710980kIvMo9oDetm3bsHLlSgwePBhRUVH45xSnLl264MSJE2WZj4jI6kiVSviMGguntkGo8d48FklEFqzYZ5S2bduGUaNG4c0338x3FVqNGjVw//79MgtHRGQtsv66CUNuDhybNAMAODRuAofGTcQNRUTPVexC6eHDh+jUqVOBjzk4OBR4OxIiospKMBiQfOA7JO3fB6mdHWrMmQ+FR8lur0RE5a/YhZKTkxMSExMLfCwmJgbu7u6lDkVEZA10aamIW78WWTeuAwAcmzWHzEklcioiKo5iF0pBQUFYv349unXrBhsbGwCARCKBTqfD119/jY4dO5Z5SCKiiibrxnXErlsDvVoNiVIJrxEj4dyBvx+JKppiF0rTpk3D4MGD0adPH7z44ouQSCTYunUrbty4gcePH2P58uVmiElEVDEIgoCkb/ci+cB3gCBAWbUafCMmwaZKVbGjEVEJFPuqtxo1auDrr79G7dq18fXXX0MQBHz77bdwdXXFV199hSpVqpgjJxFRhSCRSKDPzAQEAapOnVF99nsskogqsFLdwkSj0SAlJQXOzs6wtbUty1yi4S1MLANvDWBZ2B/PJ+j1kMhkAACDVoOsa9fg2Kx5mT8P+8JysC8shzlvYVLso/70008wGJ5+QyiVSnh7e1tNkUREVFyCXo+E3Tvx6NOlEP73u1GqUJqlSCKi8lfsOUqTJk2Ch4cH+vXrh4EDB8Lf398cuYiILJ42OQmxa9cg585tAEDmlctwbNpM3FBEVKaKfUYpKioKrVq1wpYtWxASEoJhw4bhm2++QUZGhjnyERFZpIzLf+D+/DnIuXMbUjs7+EZMZpFEZIVKPEdJrVbju+++w759+3DlyhXY2tripZdewsCBAxEUFFTWOcsN5yhZBo79Wxb2x/8TdDok7tmFlCOHAQA2NWrCN2IylF5e5fL87AvLwb6wHOaco1Sqydx5oqOjsXv3buzbtw+pqam4fv16WWQTBQsly8BfQJaF/fH/4jasg/rczwAAl24vwWPwUEgVinJ7fvaF5WBfWA5zFkrFnqP0b4IgIDY2FnFxccjIyEAZ1F1ERBbLtWcvZN64Bq//hsGpRUux4xCRmZW4ULp//z727NmDb7/9FvHx8fDy8sLo0aMxcODAssxHRCQqg1aLnDu3YV+/AQDApmo11Fr8cbmeRSIi8RS7UNq9ezf27NmDixcvQqFQoGvXrhg4cCA6duwIqdQ8p72IiMSgSXiC2DWrkPvwAfzefBt2deoCAIskokqk2IXSO++8gwYNGuCdd95B37594ezsbI5cRESiSv/tV8Rv+gKG7GxIHRxgyM0VOxIRiaDYhdK+fftQr149c2QhIhKdQatBwjfbkfbTcQCAbZ268J0wEQo3d5GTEZEYil0osUgiImuliY8zDrUBgGuvPvDoNwASeamveyGiCqpIP/2RkZEYMmQIvL29ERkZ+cxtJRIJpkyZUibhiIjKU+a1q8h9+AAyRyf4jBsPh0ZNxI5ERCIrcqHUuXNnFkpEZNVcunSDISMDzp2DIXdxFTsOEVmAIhVKN2/eLPD/iYgqstzHj5G4dxd8xoyHzM4OEokE7i/3FzsWEVkQXs9PRJVS2s9n8GDhPGReuojE3TvFjkNEFqrYhVL9+vVx+fLlAh+7evUq6tevX+pQRETmYsjNRdwX6xD/5XoIGg3s6zeAe9+XxY5FRBaq2JdyPOsWJQaDARIJ75FGRJYpN+YRYtesgib2MfC/YTa3Pn0h4WK5RFSIMr3m9dq1a3BycirLQxIRlYnMq5fxeFUkBI0GMmcX+E6YCPtALndCRM9WpEJp06ZN2Lx5M4D/v6pNqVSabJObm4ukpCT06NGj7FMSEZWSTfWakNrZw6ZuAHzGToBcpRI7EhFVAEUqlNzd3VG37tN7HMXExMDPzw+qf/2SUSqVCAgIwMiRI8s+JRFRCWiTk4wrastVKlR/6x3I3d051EZERVakQikkJAQhISEAgLCwMMybNw/+/v5mDUZEVFKCICDt5E9I2P4VvEeNgapdewCAwtNT5GREVNEUe47Sli1bzJGDiKhM6LOzEb/pS2T89gsAIPPyn8ZCiYiouIpUKD1+/Bienp5QKBR4/Pjxc7evUqVKqYMRERVXzt9/IzZqFbQJTwCZDB4DB8O1e0+xYxFRBVakQqlbt27YsWMHmjRpgq5duz53CYAbN26USTgioqIQBAGpx48hcecOCDod5O7u8J0wCXb+dcSORkQVXJEKpUWLFsHPz8/4/1wriYgsSe79v5Hw9TYAgEPzFvAZNRYyBweRUxGRNShSoTRgwADj/w8cONBsYYiISsK2Zi249ekLmZMKLt1e5B9zRFRmyuQa2dzcXERHR0Ov15fF4YiInkkQBKQcOwptUqKxzWPAILi++BKLJCIqU8UulLZs2YKVK1cav7569SqCg4MREhKCHj16IDY2tkwDEhH9kz4jA48/X46E7dsQG7UKgk4ndiQismLFLpR27txpstjk0qVL4ezsjLfffhuCIGD16tVlGpCIKE/2ndu4v2AOMi//CYlcDlX7joBMJnYsIrJixV5HKTY2FrVr1wYAZGRk4LfffsOyZcvQvXt3qFQqrFixosxDElHlJhgMSPnhEBL37gYMBii8veEbMRm21WuIHY2IrFyxCyWNRgO5/Oluf/zxBwwGA9q3f7qYW7Vq1ZCYmPis3YmIikWfmYnYdWuQdfUKAMCpTTt4jwyH1NZO5GREVBkUe+jN19cXv/32GwDgxx9/RL169eDo6AgASE5ONv4/EVFZkCgV0KelQaJQwGvkKPiMj2CRRETlpthnlF5++WWsXLkSP/74I27evIk333zT+NjVq1dRs2bNssxHRJWQYDAAACRSKaQKJXwnToag0cLmf+u5ERGVl2IXSpMmTYJcLsfFixfx4osvIiwszPjYrVu30L179zINSESViy4tDXHr18K2Th149Hu6hpvS20fkVERUWRW7UJJIJJgwYUKBj61Zs6bUgYio8sq6cR2x69ZAr1Yj5140XLu+CJmTk9ixiKgSK3ahlCcjIwN//PEHUlNT4erqiqZNm3J+EhGViGAwIOm7b5H8/X5AEKCsWg2+EZNZJBGR6EpUKG3YsAGRkZHIycmBIAiQSCSwtbXFtGnTMHr06LLOSERWTJeagth1Ucj+6yYAQNWpM7yGj4DUxkbkZEREJSiU9u3bh48//hidO3fGgAED4OXlhSdPnmDfvn346KOP4Orqiv79+5shKhFZG4NWiweLFkKXnASJjS28R4ZD1TZI7FhEREYSQRCE4uzQv39/1KlTB0uXLs332MyZMxEdHY29e/eWWcDyFpeUCej0kEl5vygxyeVSuLo6ICUlEzqdQew4lZ45+yP11Amk/fQjfCOmQOnDSdvPw58Ny8G+sBxubg6Qycrk9rX5FPuod+/excsvv1zgYy+//DKio6NLHYqIrJc2ORk5D+4bv3buFIzq78xlkUREFqnYQ2+2trZIS0sr8LG0tDTY2tqWOhQRWaeMy38i7ot1kCptUGPOfMgcHSGRSAB5ia8rISIyq2KfUWrZsiUiIyMRHx9v0p6QkICVK1eiVatWZRaOiKyDoNMhYecOPF7xKQwZGZA5OcGQmyt2LCKi5yr2n3EzZszAsGHD0L17dwQFBcHT0xMJCQk4f/485HI5IiMjzZGTiCoobVIiYqNWI+fu02F5l24vwWPwUEgVCpGTERE9X7ELpbp162LXrl2IjIzEhQsXkJqaChcXF3Tr1g2vvPIKatWqZY6cRFQBZVy6iLgv18OQlQWpvT18Ro+FY/OWYsciIiqyYhVKer0eycnJqFatGpYtW2auTERkBQRBQNrZMzBkZcG2dm34TpgEhYen2LGIiIqlSIWSIAhYtmwZtm7dipycHMhkMnTv3h0LFizgatxEVCCJRAKf8DFIrV4Dbr36QMIJ20RUARXpN9fmzZuxbt06+Pn5oWHDhnjw4AEOHjwIhUKBJUuWmDsjEVUQ6b/9iqwb1+EVOhISiQQyR0e49+0ndiwiohIrUqG0Z88eBAcHY+XKlZD/76/CTz75BBs3bsSCBQtgw1sNEFVqBq0GCd9sR9pPxwEA9g0awqklr4AlooqvSMsD/P333xg+fLixSAKAsLAwaLVaPHr0yGzhiMjyaeLj8HDRQmOR5NqrDxybNhM3FBFRGSnSGaXc3Fy4u7ubtOV9ncu1UIgqLfWF84jfvBFCbg5kjk7wGTceDo2aiB2LiKjMmOfGKMV07949jB07Fs2aNUNQUBAWLlyInJycYh3j6NGjCAwMREhIiJlSEtE/Je7Zhbh1ayDk5sAuIBDV5y5gkUREVqfIl6HMnDmzwLlI06dPh1KpNH4tkUiwf//+IgdQq9UIDw9HlSpVsGLFCiQnJ2Px4sVITU0t8Ma7BcnJycHixYvh4eFR5OclotKxb9AQyT8cgluvPnDv2w8SmUzsSEREZa5IhVLr1q2L1V4c27dvh1qtxr59++Dm5gYAkMlkmDlzJiZNmgR/f//nHiMqKgpVqlRBtWrVcPXq1VJnIqKCaRITIXV5+nNqX68+ai1aAoU7/0AhIutVpEJpy5YtZgtw6tQpBAUFGYskAOjRowdmz56NkydPPrdQevDgAb788kts374dGzduNFtOosrMkJuL2ys2IvHns6j+7jwofXwAgEUSEVk90VeAi46OxqBBg0zalEolqlevjujo6Ofu/8EHH6Bfv36oV69emWWSyiSQyyxi+lalJfvf+y9jP4guJ+YRYlatRG5MDCCRIDf6NuyrVRE7VqXFnw3Lwb6wHBKJ+Y4teqGkVquhUqnytatUKqSlpT1z3+PHj+PSpUs4fPhwmWZSOdnB1kb0t4YAqFR2YkeotARBwJMfj+PvqPUwaDRQuLoi8PXX4Ny4kdjRCPzZsCTsC+tmsdWAIAiQPKNEzM3NxaJFizB16lSTYbuyoE7PRnYW/0IQk0wmhUplB7U6G3q9Qew4lY4hJwexmzYi7dxZAIBj48aoP3M6cqRKpKRkipyucuPPhuVgX1gOZ2c7SKXm+dwWvVBSqVRQq9X52tPT0585P2nTpk2QSqXo06ePcX+tVguDwQC1Wg1bW1uTq/GKw6AXoBP4TW8J9HoDdDr2RXlLPnr0aZEklcKj/0B4hoRA6eKEzJRM9oeF4M+G5WBfiE8QzHds0Qslf3//fHORNBoNHjx4kG/u0j/dvXsX9+/fR1BQUL7HWrdujXnz5uE///lPmeclqgxcX+qBnL/vwfXF7rCrGwCJmf5SIyKydKIXSp07d8bq1auRkpICV1dXAE8Xj9RoNAgODi50v/Hjx2PAgAEmbWvXrsW9e/ewePFi1KxZ05yxiayKPjsbKUcOw71PX0jkckjkclSZ9IrYsYiIRFfiQik6Ohq//vorUlJSMHjwYHh6eiI+Ph7Ozs6wtbUt8nGGDx+OrVu3YvLkyZg8eTKSkpLw4Ycfom/fviZDb7Nnz8a+fftw/fp1AE/PRP17aG7v3r2Ij49H27ZtS/qyiCqdnPt/I3bNKmgTnkDQauE5eKjYkYiILEaxCyW9Xo/33nsPe/fuNU647ty5Mzw9PTF37lzUr18fr776apGPp1KpsGnTJixcuBBTp06Fra0tQkJCMHPmTJPtDAYD9Hp9ceMSUSEEQUDqTz8i8ZvtEHQ6yN3d4di8hdixiIgsikQQijcFKjIyEmvXrsX06dPRqVMnhISEYPfu3WjYsCG2bduGvXv3YteuXebKa3ZxSZmATg+Z1IyLMtBzyeVSuLo6IIWTh81Cn5WJ+I1fIOPi7wAAh+Yt4DNqLGQODgVuz/6wHOwLy8G+sBxubg5mW8+q2GeU9u7di8mTJ2P06NH5zvBUq1YNjx49KrNwRFT2ch7cx+NVn0OXmAjIZPAcMhwu3V585nIcRESVVbELpfj4eDRr1qzAx2xsbJCZyTVWiCyZVGkDfXoGFJ6e8I2YDNuatcSORERksYpdKLm7u+Phw4do165dvsfu3bsHn//dA4qILIeg00Eif/rjrvTxQdVXp8Ommh9k9vYiJyMismzFHtALDg7GmjVrEB8fb2yTSCRIT0/Hli1b0KVLlzINSESlk33nNu69+xaybt4wttkHBLJIIiIqgmKfUZo2bRpOnTqF3r17o23btpBIJFi2bBlu374NuVyOyZMnmyMnERWTYDAg5YfDSNy7CzAYkLR/H+wC63EuEhFRMRT7jJKHhwd27dqFPn364Nq1a5DJZLh58yY6d+6M7du3w8XFxQwxiag4dOlqxKxYjsTd3wAGA5zatEPVaa+xSCIiKqYSLTjp4eGBBQsWlHUWIioDWbf+Quza1dCnpkKiUMDrP6FQderMIomIqAREv4UJEZWdnAf38ejjDwFBgNLHF74TJ8Ommp/YsYiIKqxiF0pvv/32Mx+XSCRYtGhRiQMRUcnZ+FWHU+s2kMjk8BoRBmkxbidERET5FbtQunDhQr621NRUZGVlQaVSwcnJqUyCEVHRZP118+ml/g4OkEgk8Bkz3rgUABERlU6xf5seP368wPZz585h/vz5+Oyzz0odioieTzAYkPTdt0j+fj8cmjVHlclTIZFIWCQREZWhMrsxSlBQEEJDQ/HBBx+U1SGJqBC61FQ8+uQjJH/3LSAIkDk4ArxpNBFRmSvTPz39/f1x5cqVsjwkEf1L5rWriFsfBX16OiQ2NvAOC4eqXXuxYxERWaUyLZR+/fVXuLq6luUhieh/BL0eSd/uRfKhA4AgwMbPD74RU6DkbYOIiMym2IVSZGRkvjatVou//voLp06dwtixY8skGBGZMuTkQH3+LCAIcA7uAs/h/4FUoRQ7FhGRVSuTQkmpVKJq1aqYNm0aCyUiM5E5OMA3YjJ0yclwat1G7DhERJVCsQulmzdvmiMHEf2LoNMhcd8eKL294dwpGABg518H8Bc5GBFRJVKsq95ycnLw+uuv47fffjNXHiICoE1KwsOPP0TK4YN48vU26NJSxY5ERFQpFatQsrW1xY8//ghBEMyVh6jSy/jjEu7Pn4Oc6DuQ2tnBZ+wEyJ1dxI5FRFQpFXvorV69erh16xZat25tjjxElZag0yFh1zdIPXYEAGBbqzZ8J0yCwtNT5GRERJVXsQulmTNn4s0330TdunXRpg0nlBKVBUGnw8OPFiHn7l0AgOtLPeAxaAhX2SYiElmRfgv/+uuvaNCgARwcHDB//nxkZmYiPDwcKpUKXl5eJttKJBLs37/fLGGJrJVELod9w8bQxMXDZ8w4ODZrLnYkIiJCEQulkSNHYseOHWjSpAlcXFzg4uJi5lhE1s+g1cCQmQm5y9NFWt379oNL8AvGr4mISHxFKpT+OXl7y5YtZgtDVFlo4uMRG7UKAOD39juQKpSQSKUskoiILAwnQBCVs/RfLiB+85cw5ORA5ugEbVw8bPz8xI5FREQFYKFEVE4MGg0SdnyFtJMnAAB2dQPgM2ESFLw/IhGRxSpyoRQeHg6JRPLc7SQSCX7//fdShSKyNpq4WDxeswqaRw8BiQRufULg3rc/JDKZ2NGIiOgZilwotWnTBm5ububMQmS1nmzbAs2jh5A5qeAzPgIODRqKHYmIiIqgyIXSlClT0KRJE3NmIbJa3uGjkbBjO7xGhEHOq0aJiCqMYt3ChIiKJjcmBik/HjV+rfDwRJUpU1kkERFVMJzMTVSGBEGA+uczePLVFggaDZTePnBo1FjsWEREVEIslIjKiCEnB/FbNyH9/DkAgH3DRrCpXkPkVEREVBpFKpRu3rxp7hxEFVruw4d4HLUS2rg4QCqFR/+BcO3ZGxIpR7eJiCoynlEiKqW0n0/jydbNELRayF1d4TthEuzqBogdi4iIygALJaJSkkhlELRaODRuAp8x4yFzchI7EhERlREWSkQlYNBqIVUoAACqoPaQOTnCvkEjDrUREVkZ/lYnKgZBEJBy/BjuvzcbunS1sd2hURMWSUREVoi/2YmKSJ+VidjVkUj4aiu0iQnGe7YREZH14tAbURHk3LuL2KjV0CYmADIZPIcMg0u3l8SORUREZsZCiegZBEFA6rEjSNj1DaDXQ+HhCd+ISbCtVVvsaEREVA5YKBE9Q+rRH5DwzXYAgGPLVvAOHw2ZvYPIqYiIqLywUCJ6BlWnYKSdOQ2XLl3h/EJXSCQSsSMREVE5YqFE9A+CwYCMSxfh2KIlJBIJZHZ2qDF3ASQymdjRiIhIBLzqjeh/9OnpiFmxHLGrI5F6/JixnUUSEVHlxTNKRACybv2FuHVroEtJgUShgFSpFDsSERFZABZKVKkJBgOSD36PpG/3AoIAhY8PqkRMgY2fn9jRiIjIArBQokpLp1Yjbn0Usq5fAwA4BbWH94iRkNraipyMiIgsBQslqrS0T+KRdfMGJEolvP4bBlWHjryqjYiITLBQokrLrk5deIePhm3N2rCpWlXsOEREZIF41RtVGrrUVMSs+BS5MTHGNucOnVgkERFRoXhGiSqFzGtXEbd+LfTpaujT0+E3+z0OsxER0XOxUCKrJuj1SNq/D8kHvwcEAcpqfvAZO55FEhERFQkLJbJa2uRkxK1bg+zbtwAAzsEvwHPYf7lGEhERFRkLJbJKmtjHeLBkEQwZGZDa2sJ75Gg4tWkrdiwiIqpgWCiRVVJ4ecOmSlUYcnLgGzEZSm9vsSMREVEFxEKJrIY2OQkyJxWkCgUkMhmqTHoFEltbSBUKsaMREVEFxeUByCpk/HEJ9+fNQeLO7cY2mZMTiyQiIioVnlGiCk3Q6ZCw6xukHjsCAMi+excGjYYTtomIqEywUKIKS5uQgMdRq5D79z0AgMtLPeA5aAgkcn5bExFR2eAnClVI6b//hviNG2DIzobU3gE+Y8bBsVlzsWMREZGVYaFEFY4+MxPxm76AITsbtv514DthEhTu7mLHIiIiK2QRhdK9e/ewcOFC/P7777Czs0OfPn0wc+ZM2NraFrpPRkYGvvzyS5w6dQr37t2DXC5Hw4YNMWPGDDRs2LAc01N5kzk4wDt8DHLu3YVH/4EcaiMiIrMR/ao3tVqN8PBwZGZmYsWKFZg1axa+++47vPvuu8/c7/Hjx9ixYwfat2+PTz/9FIsXL4bBYMDw4cNx7dq1ckpP5SX9lwvIvHbV+LVTy1bwHDyURRIREZmV6J8y27dvh1qtxr59++Dm5gYAkMlkmDlzJiZNmgR/f/8C96tWrRqOHj0KOzs7Y1v79u3RrVs3bN26FYsXLy6X/GReBo0G8Vu3Iu3UCcgcnVBj/vuQO7uIHYuIiCoJ0c8onTp1CkFBQcYiCQB69OgBpVKJkydPFrqfvb29SZEEADY2NvD398eTJ0/MlpfKT9ajR7i3YD7STp0AJBI4B78AmaOT2LGIiKgSEf2MUnR0NAYNGmTSplQqUb16dURHRxfrWFlZWbhx4wb69etXqkxSmQRymeg1ZKWmPncWjzdthCEnBzKVClUnTIRjo0Zix6q0ZP/7eZDx50J07AvLwb6wHBKJ+Y4teqGkVquhUqnytatUKqSlpRXrWMuXL0d2djZCQ0NLlUnlZAdbG9HfmkpJ0OtxZ+UaPPnxOADAuXEjBMx4DUo3V5GTEQCoVHbP34jKBfvCcrAvrJvFVgOCIEBSjBLxu+++w6ZNmzBnzhzUqFGjVM+tTs9Gdhb/QhCLRqMFJBL4DR8K5559kCkAmSmZYseq1GQyKVQqO6jV2dDrDWLHqdTYF5aDfWE5nJ3tIJWa53Nb9EJJpVJBrVbna09PTy90Ive//fzzz3j77bcxduxYjBgxotSZDHoBOoHf9OVFEAQIOi2kiqe3HfEYHgqXTp1RtW0LpKRkQqdjX1gKvd7A/rAQ7AvLwb4QnyCY79iinzbx9/fPNxdJo9HgwYMHRSqULl++jFdeeQU9e/bEG2+8Ya6YZCaGnBzEfbEOj1d+DsHw9BeN1MYG9gGBIicjIiKygEKpc+fOOH/+PFJSUoxtR48ehUajQXBw8DP3jY6Oxvjx49GiRQssXry4WEN1JL7chw9xf+E8pJ87i6xrV5Fz767YkYiIiEyIPvQ2fPhwbN26FZMnT8bkyZORlJSEDz/8EH379jU5ozR79mzs27cP169fBwAkJSVh7NixUCgUGDdunMkik0qlEg0aNCj310JFIwgC0k6dRML2bRC0WshdXeEzfiLs/OuIHY2IiMiE6IWSSqXCpk2bsHDhQkydOhW2trYICQnBzJkzTbYzGAzQ6/XGr+/cuYPY2FgAwKhRo0y2rVq1Ko4fP2727FR8+uxsPNmyEem/XAAA2DdqAt+x4yFz4vpIRERkeSSCYM4pUBVPXFImoNNDJuUwnjk8Wv4Jsq5eAaRSeAwcDNfuPSEp4EoFuVwKV1cHTua2EOwPy8G+sBzsC8vh5uZgtvWsRD+jRJWLR/9BiE14Ap8x4znURkREFk/0ydxk3fRZmci8etn4tW3Nmqi5YBGLJCIiqhBYKJHZ5Ny7iwcL5uFx5Ark3P/b2C6RycSKREREVCwceqMyJwgCUo8dQcKubwC9HnIPD4Az4YiIqAJioURlSp+RgbiNG5D5xyUAgGOLlvAeNQYyeweRkxERERUfCyUqM9nRdxAbtRq65CRI5HJ4Dh0O5y7duBAoERFVWCyUqMxk/3UTuuQkKDy94DtxMmxr1BQ7EhERUamwUKIy49qzNyCRwvmFLpDZ2Ykdh4iIqNR41RuVWPbtW3j06VIYcnMBABKpFG69erNIIiIiq8EzSlRsgsGA5EMHkPTtXsBgQPLB7+ExYJDYsYiIiMocCyUqFp1ajbj1Uci6/vQmxE7tguDWq4/IqYiIiMyDhRIVWdbNG4hdtwb6tDRIlEp4/TcMqg4deVUbERFZLRZKVCTqcz8j7ov1gCBAWaUKfCOmwKZqVbFjERERmRULJSoSu3oNIHVwgGPT5vD6byikNjZiRyIiIjI7FkpUKE1cHJQ+PgAAhasras5bCLmLi7ihiIiIyhGXB6B8BL0eiXt34+/33kb6xd+N7SySiIiosuEZJTKhTUlB3NrVyL59CwCQE30bTi1aipyKiIhIHCyUyCjzymXEbVgHfUY6JDa28A4fBVWbdmLHIiIiEg0LJYKg0yFx3x6kHD4IALCpXgO+EZOg9PYRORkREZG4WCgRsv66aSySnLt0g+fQYZAqlCKnIiIiEh8LJYJDw0Zw7dUHtjVqwqlVa7HjEBERWQxe9VYJPR1q2w1daoqxzXPQEBZJRERE/8IzSpWMNiEBsWtXI+feXWTfvo1qM2fxFiRERESFYKFUiaT//hviN26AITsbUnt7uL7YnUUSERHRM7BQqgQMWi0Sd25H6vEfAQC2tf3hGzEJCncPkZMRERFZNhZKVk6bnIzHkZ8h98F9AIBrj17wGDAIEjm7noiI6Hn4aWnlZA4OELRaSB0d4TNmHBybNBM7EhERUYXBQskKGbQaSGRySKRSSG1sUGXKVEiUNlC4uYkdjYiIqELh8gBWRhMXiwcfvI+UHw4b25Q+viySiIiISoBnlKyI+vxZxG/ZBCE3FykZ6XDp2g1SGxuxYxEREVVYLJSsgCE3F0++3gb1mVMAALvAevAdH8EiiYgsksFggF6vEztGqRkMEuTkyKDR5EKvF8SOY7VkMjmkUvEGwFgoVXC5j2MQu2YVNI9jAIkEbiEvw71vP0hE/KYiIiqIIAhQq5ORnZ0hdpQyk5gohcFgEDuG1bOzc4RK5SbK2n8slCowfXY2Hn64CIasTMicneE7LgL29RuIHYuIqEB5RZKjoyuUShurWPBWJpPwbJIZCYIAjSYXGRlPb7nl7Oxe7hlYKFVgMjs7uL/cH5l//gGfcRMgd3YWOxIRUYEMBr2xSHJ0VIkdp8zI5VLodDyjZE5K5dNpJBkZKXByci33YTgWShVM7qOHgADY+PkBAFy6vQiXrt041EZEFk2v1wP4/w89ouLI+77R63WQSpXl+tz8dK0gBEFA6qkTePDBAjxeHQlDTjYAQCKRsEgiogrDGobbqPyJ+X3DM0oVgD47G0+2bET6LxcAAAovLwg6vcipiIiIrB8LJQuX8+A+YtesgvZJPCCVwmPAYLj26MmzSEREROWAn7YWShAEpP70Ix4ueh/aJ/GQu7nB78234darN4skIiILMXLkf9GxYytcvPhbvsdiYx+jY8dW+OmnY/keS01NRceOrXDw4Hcm7TqdDrt378D48eF46aXO6Nq1PUJDh2Lz5i+Qnp5uttfxT0lJiZgz52107x6Mnj1fwPvvz4Fanfbc/QRBwLZtmzBkyMvo0iUIYWFD8eOPR/Jtl5GRgSVLPkCfPt3QrVsHvPLKBNy+/Zc5XkqZ4BklSyUIyPjjEgSdDg5Nm8Fn9DjIHB3FTkVERP9z//7fuHXrJgDg6NHDaNGiVamOp9Fo8MYbr+Hy5Uvo338wxoyZABsbG9y5cwt79uzCo0cPMXv23LKIXiidTofXX58GnU6L996bD51Oh1WrPsdbb72OlSvXPXOu0FdfbcbatasQHj4WjRs3wenTJzFv3juwsbFFx46djdvNn/8Obty4jkmTpsHNzR07dnyFadMmYePGr+Dt7WPW11cSLJQslEQqhc/YCcj47Rc4d+nGCZBERBbmyJFDkMlkaNasJX766UfMmDELCoWixMfbsCEKFy/+io8//gzt2rU3trdo0QoDBgwp8KxVWTt58jju3LmFzZt3oHZtfwCAh4cnJk0aiwsXzpnk+ietVotNm77A4MHDMWbMBABA69btEBcXi3XrVhsLpatXr+DcuZ/x4YfLjG0tWrTCkCEv4+uvt+K112aa/TUWF8dwLIQgCEg5dgTx2zYb2+QqFVy6vsgiiYjIAh09ehgtW7bG8OH/RUZGOs6d+7nEx8rNzcWePTvRqdMLBRYjCoUCbdsGlSZukZw79zP8/esaiyQAaNy4KXx9q+DcuTOF7hcT8whZWZn5MrZtG4To6NuIi4sDANy+/RckEgnatGln3MbW1hZNmzbDzz+fLuNXUzZ4RskC6DMyELdxAzL/uAQAcGrZGvb16ouciojI/ARBgEYr3oKNSoW0RH+MXr16BY8fx2DMmHFo3bodXFxccOTIIXTu/EKJcty8eQPZ2VkICupQov2Bp8NmzyOVSp+5YOP9+/dQs2bNfO01a9bC33//Xeh+Gk0uAEAuNy0rFAql8bg+Pj7QaHILzKBQKBEX9xi5uTmwsbF97usoTyyURJYdfQexUauhS06CRC6Hx9DhsAusJ3YsIiKzEwQBi7dexJ2Y508UNpc61Zzx9ogWxS6Wjh49BKVSiRde6Aq5XI4uXV7CgQP7kZmZAQeH4s8nTUx8AgDw8vIu9r55Xnih3XO36dUrBO+8M6/Qx9PT0+Ho6JSv3clJhb//vlvoftWq+UEqleLGjWsmc7WuXbsCAMbJ4H5+NaDX63Hr1k00aNAIwNObJN+4cR2CICA9PYOFEj0lGAxIOfoDEvfsAvR6KDy94BsxGbYFVPJERFarAs4s0Ov1OH78GIKCOsDR0Qk6nQHdu/fE3r07cfLkT+jdu2+xjykIT+8XV5qpFuvXb37uNs7OLs/dpqAMT/MVns3e3gE9evTGtm2bUbt2HTRs2Bg//3wKx479AADGM0ht2rRDtWrVsXTpYrzzzny4ublh69aNiI2N+d92lvcNwUJJJPFfboD6f+PZjq3awDt8NGR2diKnIiIqPxKJBG+PaFHhht5+/fUCUlKS0aFDZ6Snp0OnM6Bmzdrw8vLGkSOHjIWSTCYD8PSMyb8ZDE8XDc4bqvL0fHomKT4+rsSvpU6dgOdu87z7pDk5OSE9XZ2vPSMjHU5Oz75H39Sp05GUlIQ33ngVAODi4oJx4yZh5crlcHN7ejNbuVyO999fjDlz3kZ4+HAAgL9/HQwZ8h/s2rUdKpXl3bOUhZJIHFu3Qfrvv8Jz2H/g3PkFTtgmokpJIpHARikTO0axHD16CACwaNF8LFo03+SxxMQEJCUlwt3dA87OzpBKpUhKSsp3jKSkRACAq6srAKBevfqwt3fA+fM/o2/f/iXKVRZDbzVq1MLt27fytf/99z20b9/xmcdWqZyxbNnnSExMgFqdhmrVquPMmVNQKBQICAg0ble3biC++mo3Hj16CEEQ4OdXHcuWfYTAwPr55jhZAstLZKUEgwHaJ/FQ+vgCABybNEWtxR9D7mx51TMRERUsJycHp06dRKdOL2DIkOGQyaTQ65+eMUpNTcWcOW/hxx+PYOjQ/8LGxhb16jXAmTMnMXTof0yOc/r0SSiVNqhXryEAwMbGBgMGDMbXX2/Br7+eR+vWpkWPTqfDxYu/mVwt9m9lMfQWFNQBP/xwEH//fQ81a9YC8HTiemzsYwQFPbtQyuPh4QkPD0/o9Xrs27cLXbu+lG/elkQigZ9fdQBASkoKjh8/gkmTphXp+OWNhVI50KnViNuwFjn37qLG3AVQuHsAAIskIqIK5syZk8jOzsKQIcPRokUryOVS6HT/P7T29dcNcOTIYQwd+l8AwNixEXjjjVcxe/Yb6NGjN2xsbHDx4q/YseMrjBgRDien/584PXZsBG7evI5Zs2ZgwIDBaN26HZRKJe7di8aePTvRsGHjZxZK9eo1KPXrCw7uCn//unj33VmYOHEK9Ho9Vq78DE2aNDO59P/LL9dh48b12LFjH3z+dwLgyJFDyM3NRdWq1ZCYmIj9+/fg8eMYzJmz0OQ5Nm3agGrV/ODq6oYHD+5jy5YvERhYv0Rzu8oDCyUzy7p5A7HroqBPS4VEqUTuw4fGQomIiCqWI0cOw9vbB82btyzw8Z49Q/Dppx/hwYP7qF69Btq2DcInn6zApk1fYOHCOdDpdKhWzQ9Tp07HoEHDTPZVKpX45JPPsW/fLhw+fBD79++FXq9H1arVEBzcFcOG/dfsr08ul2Pp0hX47LOlWLBgDiQSoGPHzpg27XWTKSIGgwF6vd44CR14OuF7+/atiI19DDs7O7Rr1wFz5rwPDw/Tz7z09HSsXPkZUlKS4e7ugR49eiM8fOxz50+JRSL881US4pIyAZ0eslLOvBcMBiQf+A5J+/cBggClbxX4TpwMm6rVyiaolZPLpXB1dUBKSqbJX2skDvaH5aiofaHVapCUFAt3d1/j2jrW4N9nlMg8nvf94+bmAJnMPIUWzyiZgS4tFbHropB98wYAQNW+I7xGhEFqYyNyMiIiIioOFkpmkHL0CLJv3oBEqYR3aDhU7Uu+0ioRERGJh4WSGbj36w9dSgrc+vSFTZUqYschIiKiErLMmVMVjDYlBQk7t0P436JiUoUSvuMjWCQRERFVcDyjVEqZVy8jbv066DPSIbWzh3vIy2JHIiIiojLCQqmEBJ0Oid/uRcqhAwAAG7/qcGrdRuRURESWjRdaU0mI+X3DQqkEtMlJiI1ajZzoOwAA5y5d4Tl0OKRWdMkrEVFZyrvvmUaTC6WSVwBT8Wg0uQAAmaz8yxYWSsWUef0aYqNWwZCZCamdHbzDx8CpVWuxYxERWTSpVAY7O0dkZKQAAJRKG6u4x6XBIIFez7Nk5iIIAjSaXGRkpMDOzlGURSlZKBWTXKWCoNHApmYt+EZMgtLTS+xIREQVgkrlBgDGYskaSKVSGAxccNLc7Owcjd8/5Y2FUhEYcnONi0XaVPNDtdffhE2NmpAqFCInIyKqOCQSCZyd3eHk5Aq9Xid2nFKTySRwdrZHWloWzyqZkUwmF/X2JhZRKN27dw8LFy7E77//Djs7O/Tp0wczZ86Era3tc/fdu3cvoqKiEBMTgxo1amDKlCno1atXmWVLv/g7nmzZiCqvvAo7/zoAALs6dcvs+ERElY1UKoVUWvHndMrlUtja2iI7W8/bmFgx0QsltVqN8PBwVKlSBStWrEBycjIWL16M1NRULF269Jn7Hj58GG+99RYmTJiADh064NixY5g+fTqcnJzQsWPHUuUyaLVI3LkDqcePAQBSjv5gLJSIiIiochC9UNq+fTvUajX27dsHN7en448ymQwzZ87EpEmT4O/vX+i+n332GXr27InXX38dANCuXTvcu3cPK1asKFWhpH3yBPHrViP3/t8AANcePeExYHCJj0dEREQVk+grc586dQpBQUHGIgkAevToAaVSiZMnTxa638OHD3H37l2EhISYtIeEhODy5ctITk4uUR6FQYeHC+ci9/7fkDo6osq01+A5ZDgkctFrSiIiIipnon/6R0dHY9CgQSZtSqUS1atXR3R0dKH73b17FwBQu3Ztk3Z/f38IgoC7d++aFF9F5erqgBafLYNEIYdc5QyIOIGsMsu7atjZ2Q5cn0587A/Lwb6wHOwLyyGVmm+pCdELJbVaDZVKla9dpVIhLS2t0P3yHvv3vs7OziaPF5dULoetj3eJ9qWyJ+aVDpQf+8NysC8sB/vCulls7wqCUKTFyP69Td4y59awkBkRERGJS/RCSaVSQa1W52tPT08v8ExTnsLOHOUd61n7EhERERWF6IWSv79/vrlIGo0GDx48eOYVb3lzk/LmKuWJjo6GRCLJN3eJiIiIqLhEL5Q6d+6M8+fPIyXl/5e0P3r0KDQaDYKDgwvdz8/PD7Vr18bBgwdN2r///ns0adKkRBO5iYiIiP5J9EJp+PDhcHJywuTJk3H69Gns27cP77//Pvr27WtyRmn27Nlo0KCByb7Tpk3DoUOH8Omnn+LChQtYtGgRfv75Z0ybNq28XwYRERFZIdGvelOpVNi0aRMWLlyIqVOnwtbWFiEhIZg5c6bJdgaDAXq93qStV69eyMnJwZo1a7BhwwbUqFEDn376aalX5SYiIiICAIkgcPUHIiIiooKIPvRGREREZKlYKBEREREVgoUSERERUSFYKBEREREVgoUSERERUSFYKBEREREVotIUSvfu3cPYsWPRrFkzBAUFYeHChcjJySnSvnv37kXPnj3RuHFjhISE4NChQ2ZOa91K0hcZGRn4/PPPMWTIELRq1Qrt2rXD2LFjce3atXJKbZ1K83OR5+jRowgMDERISIiZUlYepemP1NRUzJs3Dx07dkTjxo3Ro0cPbN++3cyJrVdJ+yIrKwtLly7Fiy++iKZNm6J79+74/PPPodFoyiG1dbp//z7mzJmDfv36oUGDBsX6XVMWn9+iLzhZHtRqNcLDw1GlShWsWLECycnJWLx4MVJTU7F06dJn7nv48GG89dZbmDBhAjp06IBjx45h+vTpcHJy4sKWJVDSvnj8+DF27NiBQYMGYdq0adDpdNi8eTOGDx+O7du3o2HDhuX4KqxDaX4u8uTk5GDx4sXw8PAwc1rrV5r+yMzMRFhYGGxsbDB79my4u7vj/v370Gq15ZTeupSmL+bNm2f8nKhbty4uX76MFStWIC0tDe+++245vQLrcvv2bZw8eRJNmzaFwWBAUZd/LLPPb6ESiIqKEpo2bSokJSUZ2/bv3y8EBAQId+7ceea+PXv2FKZNm2bSNmbMGGHIkCFmyWrtStoXmZmZQlZWlklbTk6O0KFDB+Gtt94yW15rVpqfizzLly8XRowYIcyaNUvo06ePuaJWCqXpj08++UR48cUXhezsbHPHrBRK2hdarVZo3Lix8Nlnn5m0z507VwgKCjJbXmun1+uN/1+c3zVl9fldKYbeTp06haCgIJMb5fbo0QNKpRInT54sdL+HDx/i7t27+U7zhYSE4PLly0hOTjZbZmtV0r6wt7eHnZ2dSZuNjQ38/f3x5MkTs+W1ZiXtizwPHjzAl19+yb+Sy0hp+mP37t0YPHgwbG1tzR2zUihpXwiCAL1eDycnJ5N2lUpV5LMglJ9UWvxSpSw/vytFoRQdHW1yg10AUCqVqF69OqKjowvd7+7duwCA2rVrm7T7+/tDEATj41R0Je2LgmRlZeHGjRv5+oeKprR98cEHH6Bfv36oV6+euSJWKiXtj4cPHyIxMREqlQoRERFo1KgR2rZti/nz5xd7vhk9VdK+UCgUGDhwILZs2YI///wTmZmZOH/+PL755huMGDHC3LHpH8ry87vSzFFSqVT52lUqFdLS0grdL++xf+/r7Oxs8jgVXUn7oiDLly9HdnY2QkNDyypepVKavjh+/DguXbqEw4cPmytepVPS/khMTAQAfPTRR+jZsyfWrVuHO3fuYNmyZdBqtVi4cKHZMlur0vxszJs3D3PnzsXQoUONbWFhYXjllVfKPCcVriw/vytFoVQYQRAgkUieu92/t8k7hVqUfaloitoXeb777jts2rQJc+bMQY0aNcyYrPJ5Xl/k5uZi0aJFmDp1qsnQBJnH8/rDYDAAePqX8uLFiwEAQUFB0Ol0+Oijj/Dqq6/C09OzXLJau6L8nlq6dClOnDiB999/H7Vq1cK1a9ewYsUKqFQqTJs2rZySUp6y+PyuFENvKpUKarU6X3t6enqBfzXkKazyzDvWs/algpW0L/7p559/xttvv42xY8fydHYplLQvNm3aBKlUij59+kCtVkOtVkOr1cJgMECtVvMy6BIqaX+4uLgAANq1a2fS3q5dOxgMhmIPaVPJ++LWrVv44osvMH/+fAwdOhStW7fGqFGj8OqrryIqKgpJSUnmjE3/UJaf35WiUPL398/3y0Kj0eDBgwf5xqH/KW9s899jmdHR0ZBIJJwbUwIl7Ys8ly9fxiuvvIKePXvijTfeMFfMSqGkfXH37l3cv38fQUFBaN26NVq3bo3vv/8e0dHRaN26NXbv3m3u6FappP3h5+cHhUKRrz3vL+eSTISt7EraF3fu3AEA1K9f36S9fv360Ol0iImJKfuwVKCy/PyuFD9BnTt3xvnz55GSkmJsO3r0KDQaDYKDgwvdz8/PD7Vr18bBgwdN2r///ns0adKEww4lUNK+AJ5+g48fPx4tWrTA4sWLOfRZSiXti/Hjx2Pz5s0m/zp27IiqVati8+bN6Nq1a3nEtzol7Q+lUokOHTrg3LlzJu3nzp2DXC5HnTp1zJbZWpW0L6pWrQoA+RbCvXr1KgCgWrVqZkhLBSnTz+9iLSZQQaWlpQmdOnUShg8fLpw6dUrYu3ev0LZtW+H111832e7tt98W6tevb9J28OBBITAwUFi2bJlw/vx54YMPPhACAwOF06dPl+dLsBol7YvExEQhODhY6NChg3D27Fnh0qVLxn/Xrl0r75dhFUrzc/FvXEep9ErTH3/++afQsGFD4Y033hBOnz4tfPnll0LTpk2FDz74oDxfgtUoaV/odDph8ODBQlBQkPDVV18J586dE9auXSs0a9ZMeO2118r7ZViNrKws4dChQ8KhQ4eE0NBQITg42Ph13lpX5vz8rhSTuVUqFTZt2oSFCxdi6tSpsLW1RUhICGbOnGmyncFggF6vN2nr1asXcnJysGbNGmzYsAE1atTAp59+ylW5S6ikfXHnzh3ExsYCAEaNGmWybdWqVXH8+HGzZ7c2pfm5oLJXmv5o0qQJoqKi8Mknn2DixIlwcXFBaGgoXn311fJ8CVajpH0hk8mwZs0afPbZZ1i3bh0SExPh6+uL0NBQTJw4sbxfhtVISkrK972c9/XmzZvRtm1bs35+SwSBq2ARERERFaRSzFEiIiIiKgkWSkRERESFYKFEREREVAgWSkRERESFYKFEREREVAgWSkRERESFYKFEREREVAgWSkRERESFYKFEVEHt2bMHgYGBBf5bsmRJkY/z6NEjBAYGYs+ePWZMW/Bz5v2rV68e2rZti/Hjx+PSpUtmec6wsDCEhYUZv87Ozsbnn3+OCxcu5Ns277199OiRWbIU5sKFCybvS/369dGuXTtMnDgRV65cKfFxt23bVq79S2RNKsUtTIis2eLFi/PdCdvLy0ukNMUTFhaGkJAQ6PV63LlzB5GRkRg5ciR27NiBBg0alOlzzZ071+Tr7OxsREZG4pVXXkHbtm1NHnvhhRewY8cO0d7HGTNmoG3bttDpdLh+/TpWrlyJsLAw7Nu3DzVr1iz28b7++mu4urpi4MCBZR+WyMqxUCKq4OrWrYvGjRuLHaNEfH190axZMwBAy5YtUb16dYwaNQpfffUVFi5cWKbPVadOnSJv6+bmVry7i5exGjVqGN+XVq1aQaVSYdasWdi/fz+mTZsmWi6iyohDb0RW6v79+3j77bfRvXt3NG3aFJ06dcLEiRPx119/PXff5ORkvPfeewgODkajRo3Qrl07DB8+HGfPnjXZ7uzZswgPD0eLFi3QtGlTDB8+HOfOnStx5rzi4PHjx8a2Xbt24eWXX0bjxo3Rpk0bTJkyBdHR0Sb7PXz4ENOnT0fHjh3RqFEjtG/fHuHh4bhx44Zxm38OvT169AhBQUEAgMjISONQ11tvvQUg/9DbBx98gGbNmiEjIyNf5tdeew3t27eHVqs1th08eBDDhg1Ds2bN0Lx5c4wdOxbXr18v8fvSqFEjAEBiYqJJe2RkJIYMGYI2bdqgRYsWGDBgAHbu3Il/3sKza9euuH37Nn755Rfj6+zatavx8YyMDCxZsgRdu3ZFo0aN0KlTJ3zwwQfIysoqcV4ia8IzSkQVnMFggE6nM2mTy+V48uQJXFxc8Prrr8PNzQ1paWnYu3cvhg4dir179+YbrvunN954A9evX8f06dNRs2ZNqNVqXL9+HampqcZtvv32W8yaNQvdunXDkiVLIJfLsWPHDowdOxYbNmwwFiLFcf/+fQCAq6srACAqKgrLli1DSEgIXn/9daSkpCAyMhLDhg3Drl27jMNQ48ePh8FgwBtvvIEqVaogJSUFly5dglqtLvB5vLy8sH79eowbNw6DBw/GkCFDAKDQs0iDBg3C5s2bcejQIeO2AKBWq/Hjjz9ixIgRUCgUAIA1a9Zg+fLlGDhwICZNmgStVosNGzZgxIgR2LlzZ7HObOXJK9hq1apl0h4TE4Nhw4ahSpUqAIA//vgDCxcuRHx8PF555RUAT4upadOmwcnJyTj8qFQqATwdfgwNDUVcXBwmTpyIwMBA3L59GytWrMCtW7ewceNGSCSSYuclsioCEVVIu3fvFgICAgr8p9Vq822v0+kEjUYjdO/eXVi0aJGx/eHDh0JAQICwe/duY1uzZs2EDz74oNDnzsrKEtq0aSNERESYtOv1euHll18WBg8e/Mzsec+5du1aQavVCrm5ucLVq1eFQYMGCQEBAcKJEyeEtLQ0oUmTJsL48eNN9n38+LHQqFEjYcaMGYIgCEJycrIQEBAgbNy48ZnPGRoaKoSGhhq/TkpKEgICAoQVK1bk2zbvvX348KGxbcCAAcKwYcNMttu2bZsQEBAg/PXXX8ZsDRo0EN5//32T7TIyMoQOHToIr7766jMznj9/XggICBAOHDggaLVaITs7W/j999+FHj16CL179xbS0tIK3Vev1wtarVaIjIwU2rRpIxgMBuNjffr0MXnteaKiooR69eoJly9fNmk/fPiwsR+IKjueUSKq4JYsWQJ/f3+TNrlcDp1Oh/Xr12P//v148OCBydDQv4eu/q1JkybYu3cvXFxc0L59ezRs2NB4xgQALl26hNTUVAwYMCDf2axOnTph/fr1yMrKgr29/TOfZ+nSpVi6dKnxaw8PDyxYsADBwcE4efIkcnJyMGDAAJN9fH190a5dO5w/fx4A4OLigurVq2PDhg0wGAxo27Yt6tWrB6m0bGcWDBw4EO+//z7u3r1rPBu3Z88eNG7cGAEBAQCAM2fOQKfToV+/fibvi42NDVq3bl3gFXYFmT59usnXnp6e2L59O1QqlUn7uXPnEBUVhStXruQbFkxKSoKHh8czn+enn35C3bp1Ub9+fZO8HTt2hEQiwS+//ILg4OAiZSayViyUiCo4f3//Aidzf/jhh9i2bRvGjx+P1q1bw9nZGRKJBO+++y5yc3OfecxPP/0Uq1evxq5du/DZZ5/B3t4eL730Et544w14enoa58o8a2JxWlracwulkSNH4uWXX4ZUKoVKpUK1atWMQz15w3yenp759vPy8jLOl5JIJNi4cSNWrlyJ9evX48MPP4SLiwv69u2L1157DY6Ojs/MUFR9+/bFkiVLsHfvXrz++uu4c+cOrly5YnI1Xd77Mnjw4AKPUdTibebMmWjXrh1ycnJw5swZrF27FlOmTMHOnTuNw2aXL1/G2LFj0aZNG7z//vvw8fGBQqHAsWPHsGbNGuTk5Dz3eZKSknD//n00bNiwwMdTUlKKlJfImrFQIrJS+/fvR//+/TFjxgyT9pSUlHxnJv7Nzc0N77zzDt555x08fvwYx48fxyeffIKkpCRs2LDBOIfovffeQ9OmTQs8hru7+3Mz+vj4FHrFnouLCwAgISEh32NPnjwxZgCAqlWrYtGiRQCAe/fu4dChQ4iMjIRGo8GCBQuem6MonJ2d0a1bN+zbtw+vvfYadu/eDRsbG4SEhBi3ycu0YsUK47yhkvDz8zO+L61bt4atrS2WL1+OLVu2YOzYsQCAAwcOQC6XIyoqCjY2NsZ9jx07VuTncXV1hY2NjfG9K+hxosqOhRKRlZJIJCbDZQBw4sQJxMfHo0aNGkU+TpUqVRAaGopz587h4sWLAIAWLVpApVLhzp07CA0NLdPceZo3bw5bW1vs378fvXr1MrbHxcXh/Pnz6NGjR4H71apVC5MnT8aRI0eeeaVZ3pmZopx5yTNw4EAcOnQIJ0+exHfffYeXXnrJpOjs2LEj5HI5Hjx4UGi+khg3bhz27t2LtWvXYtiwYXB0dIREIoFMJjM5S5WTk4P9+/fn21+pVBb4Ol944QVERUXBxcUFfn5+ZZaXyJqwUCKyUi+88ILx6rbAwEBcu3YNGzZsgI+PzzP3S09Px8iRIxESEoLatWvDwcEBV65cwenTp/HSSy8BABwcHPDuu+/irbfeQlpaGnr06AF3d3ckJyfj5s2bSE5Oxvz580uVX6VSYfLkyVi2bBnefPNN9OnTB6mpqVi5ciVsbGyMV3XdvHkT77//Pnr27IkaNWpAoVDg/Pnz+OuvvzBhwoRCj+/o6IiqVavixx9/RFBQEJydneHq6opq1aoVuk/Hjh3h4+OD+fPnIyEhId8CjtWqVcO0adOwfPlyPHz4EJ07d4ZKpUJiYiKuXLkCOzu7Eq2DpFAoMH36dLz22mvYvHkzJk+ejODgYHz55Zd4/fXXMWzYMKSmpmLDhg3GAvCfAgICcODAARw8eBDVqlWDjY0NAgMDER4ejiNHjiA0NBSjRo1CYGAgDAYDYmNjcebMGYwZM6bQM4ZElQULJSIr9c4770Aul2Pt2rXIyspCgwYN8Pnnn+Ozzz575n42NjZo0qQJvv32W8TExECn08HX1xfjx4/HuHHjjNv169cPVapUwfr16zF37lxkZmbCzc0N9evXzzcBu6QiIiLg5uaGLVu24ODBg7C1tUWbNm0wY8YM49IAnp6eqF69Or766ivExcUBeDp0NWvWLJNblhTkgw8+wEcffYRJkyZBo9FgwIAB+PDDDwvdXiqVon///lizZg18fX0LXAIhIiIC/v7+2Lx5Mw4cOACNRgNPT080atQI//nPf0r8XvTq1QtffvklNm7ciLCwMAQFBWHRokVYt24dJk6cCG9vbwwdOtQ4bPpPU6dORUJCAt59911kZmaiatWqOH78OOzt7bFt2zasXbsWO3bswKNHj2BrawtfX1+0b98eVatWLXFeImshEYR/rExGREREREZcmZuIiIioECyUiIiIiArBQomIiIioECyUiIiIiArBQomIiIioECyUiIiIiArBQomIiIioECyUiIiIiArBQomIiIioECyUiIiIiArBQomIiIioEP8HsV3Npc1/VdMAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "phatTesting = fit.predict(xtesting)\n",
    "\n",
    "## See here for plotting\n",
    "## https://stackoverflow.com/questions/25009284/how-to-plot-roc-curve-in-python\n",
    "fpr, tpr, threshold = skl.metrics.roc_curve(ytesting, phatTesting)\n",
    "roc_auc = skl.metrics.auc(fpr, tpr)\n",
    "\n",
    "# method I: plt\n",
    "import matplotlib.pyplot as plt\n",
    "plt.title('Receiver Operating Characteristic')\n",
    "plt.plot(fpr, tpr, 'b', label = 'AUC = %0.2f' % roc_auc)\n",
    "plt.legend(loc = 'lower right')\n",
    "plt.plot([0, 1], [0, 1],'r--')\n",
    "plt.xlim([0, 1])\n",
    "plt.ylim([0, 1])\n",
    "plt.ylabel('True Positive Rate')\n",
    "plt.xlabel('False Positive Rate')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "iaXw2nFg6UoO"
   },
   "source": [
    "## Aside different python packages\n",
    "\n",
    "So far we've explored several plotting libraries including: default pandas methods, matplotlib, seaborn and plotly. We've also looked at several fitting libraries including to some extent numpy, but especially scikitlearn and statsmodels. What's the difference? Well, these packages are all mantained by different people and have different features and goals. For example, scikitlearn is more expansive than statsmodels, but statsmodels functions more like one is used to with statistical output. Matplotlib is very expansive, but seaborn has nicer default options and is a little easier. So, when doing data science with python, one has to get used to trying out a few packages, weighing the cost and benefits of each, and picking one. \n",
    "\n",
    "'statsmodels', what we're using above, has multiple methods for fitting binary models including: `sm.Logit`, `smf.logit`, `BinaryModel` and `glm`. Here I'm just going to use `Logit` which does not use the formula syntax of `logit`. Note, by default, this does not add an intercept this way. So, I'm adding a column of ones, which adds an intercept.\n",
    "\n",
    "Consider the following which uses the formula API\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 420
    },
    "colab_type": "code",
    "id": "yj_2oMxvWC4S",
    "outputId": "22d63911-ebca-429a-b14b-af07c39a59bb"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.264729\n",
      "         Iterations 8\n",
      "                           Logit Regression Results                           \n",
      "==============================================================================\n",
      "Dep. Variable:           GOLD_Lesions   No. Observations:                   75\n",
      "Model:                          Logit   Df Residuals:                       67\n",
      "Method:                           MLE   Df Model:                            7\n",
      "Date:                Mon, 29 Jan 2024   Pseudo R-squ.:                  0.6176\n",
      "Time:                        06:47:16   Log-Likelihood:                -19.855\n",
      "converged:                       True   LL-Null:                       -51.926\n",
      "Covariance Type:            nonrobust   LLR p-value:                 2.236e-11\n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "Intercept     -4.8293      1.876     -2.574      0.010      -8.507      -1.152\n",
      "FLAIR          2.2957      1.166      1.969      0.049       0.011       4.580\n",
      "T1             2.3580      1.055      2.236      0.025       0.291       4.425\n",
      "T2             2.1419      1.084      1.976      0.048       0.017       4.267\n",
      "FLAIR_10       7.7567      3.677      2.109      0.035       0.550      14.964\n",
      "T1_10          1.2147      1.572      0.773      0.440      -1.866       4.296\n",
      "T2_10         -4.9682      3.051     -1.628      0.103     -10.949       1.012\n",
      "FLAIR_20     -18.1028      7.771     -2.330      0.020     -33.334      -2.872\n",
      "==============================================================================\n"
     ]
    }
   ],
   "source": [
    "results = smf.logit(formula = 'GOLD_Lesions ~ FLAIR + T1 + T2 + FLAIR_10 + T1_10 + T2_10 + FLAIR_20', data = trainingDat).fit()\n",
    "print(results.summary())\n"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "include_colab_link": true,
   "name": "Copy of notebook5.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}