Torch Exploration.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import torch\n",
    "# torch.cuda.is_available()\n",
    "# We first import all the packages we will need for the project\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pylab as plt\n",
    "import seaborn as sns\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "sns.set_style('white')\n",
    "plt.set_cmap('gist_gray')\n",
    "\n",
    "\n",
    "# We also install a package to read NiFTI files\n",
    "import nibabel as nib\n",
    "\n",
    "import torch\n",
    "import torchvision.models as models\n",
    "import torch.optim as optim\n",
    "import torch.nn as nn\n",
    "import torch.nn.functional as F"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We open the demographics table using the pandas package (pd)\n",
    "df = pd.read_csv('demographics.csv')\n",
    "# We display the 10 first lines\n",
    "subjects_id = []\n",
    "missing_subjects = []\n",
    "#figure out the subject ids of the invalid scans\n",
    "for i in df.path_process_data[0:563]:\n",
    "    subject_num = int(i[3:6])\n",
    "    subjects_id.append(subject_num)\n",
    "for j in range(663):\n",
    "    if j not in subjects_id:\n",
    "        missing_subjects.append(j)\n",
    "# print(missing_subjects)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are 1597 patients\n",
      "There are 25 hospitals\n",
      "There are 2 datasets\n",
      "There are 55.48% of women\n",
      "Age statistics : \n",
      "count    1597.000000\n",
      "mean       35.599132\n",
      "std        17.491271\n",
      "min        18.000000\n",
      "25%        22.000000\n",
      "50%        27.000000\n",
      "75%        48.052019\n",
      "max        86.318960\n",
      "Name: age, dtype: float64\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABOAAAAJZCAYAAAAArkr3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3XmcXWV9+PHPAIGAQJgEZAsQHeTrgoAsblWLxBA33DcIClq3KGgro7XWhdYFa2MVl8atFZRR0VYFFJhMEzdcCqhsVr6290dYBEQyEwQJ0Unm98c5A5fJ7DN37vZ5v155zb1neZ7vOfdmznee8zzP6RgaGkKSJEmSJElSbWxX7wAkSZIkSZKkVmYDnCRJkiRJklRDNsBJkiRJkiRJNWQDnCRJkiRJklRDNsBJkiRJkiRJNWQDnCRJkiRJklRDNsBJ2kZEHBgR90TE9vWOpZFExGci4j11qvvUiLis6v09EfHwWSr7XRHxhfL1kogYiogdZqnspvwuzfZ5mGSdT42InMH+50TEB2YzJklSY2rW62ujqb52zvQ6PErZl0TEKeXrB+Vxs1D2iohYM1vlzZXZPg+TrPP+PHea+6+PiGfMZkxqX3P2h4WkxhMR64G9gS1Viw/JzJuAXesR00gRsSPwFeBo4CDg6Zn5/SnsPwTcCwwBdwHnA2/PzC0T7Hcq8NrMfMrwssx841TjH6PsY4HzMnPxdMvIzAk/n8nWk5kfmm4co9S5nuK8/VdZdsN8lxpdZv4IiOH3I8+lJKn9NEmu9kTg/cBRFHF+H3hLZt42hTI6gApwX2Y+uhZxTmTkdXgsEXEmcHBmnjxBec+ajbgiYglwAzAvMwfLsnuAntkov9VV57mjnUtpLtkDTtIJmblr1b9ba1nZNHsUXQacDNw+zWoPLxuslgInAa+bZjktZy57eDUqz4EkqcE1eq7WCXwOWEJxs/Ru4ItTLONpwEOBh0fEMVPct6FEREdEtP3f2eZX0rb8TyFpGyPvDkXEw4BzgccB/w0ksCAzTx6tl1V1z53yLuGhwH3A84C3RcS/A++gaAjbA1gLvDEz+0fGkpl/Aj5eljtur7WJZOb1EfGjMh4i4p1lDA8Fbgb+PjO/FRGPAj4DzIuIe4DBzNwjIs4BbsnMd5f7Pxf4AEXC+T/lMVxTdQ4+BbyKIhm9FDgF2B64BNipLBuKO9kPSqYjYhFF8noscD3QO2L9EPCIzPy/iHg2sAo4APgD8DFg9Wj1AK9n289jMdveyX1N+dl1AKsy86NlvSPPwbGUn39EfBk4ELio/Kz+Efg6D/4u7Vee26cA/cA/Zebny7LOBB5dxvZC4CbglMy8klGU5+CtwF8Du5fn628zc2u5/jXA24F9gMuB12fmjVX7nlbuuwPwsNHqAFZExPuBXYCPZeYHy/13Av4JeFm53dfLujdHxJ7AOeUxbgV+BfxlZm4tvxefBV4J7At8G1iZmfdNdC4z8yMR8Q3gqcDOwNXlvr8aI3ZJUotqsFztkhGxfQr4wRQP6RTgAorr2ynAFVXljXls5fonAv9CkUPcCLx1rNESEfE44N+ARwAXU4yQGF53LFXnKSL+FngLRY5xK/AmYB7wLqAjIl4AVDLz8Ij4PvBjirztSOCx5bDH8zJzePhjR0R8kiI3vA14c2auLetaT1Wv9xG97H5Y7r8xIgCWUfTUu3+kRkQ8GTibItf7TXkOflKu+z7wI+A44DDgp8BJmXnnKOfnWOA84F+BtwH3UOTHPeX6nYAPUuQ/OwHfAv4mMzdV7ftJ4G+APop8Z7TPYRXwV8BG4E3D36EJ8sTHl3EdAmwCejLzbVX/F94AnMm2uetE5/IO4PPA4RTfh97ys9k4WuzSTLR9y7ykSfkKRQPGIooL26gX03E8H/gPigSuhyKZeQHwl8B+wADw6VmKdUwR8WiKxotflosq5fsFwD8A50XEvpn5a+CNwE/LO817jFLWkcC/U1zsF1E0qlxYJibDXgY8k6Jx5zDg1Mz8I/As4NYJ7mR/miIR3hd4TflvLP8GvCEzd6NIoNdNUM/Iz2M0T6dITo8H3jmZuS8y85UUjWbDd+o/MspmXwVuofjcXwJ8KCKWVq1/HvC1MrYLKRoxx/NCiuHJR5bH9RqAMil+F/AiYC+KxPOrI/Z9AfAEioR9LE+hSHKXAu8tG2cB/h54InAERcL2eODd5bozymPci2LY0LuoSvCBFcByoIsiiXw3I4xzLi+h+FweCvwCh59IkgqNlKs9jeLm06RExC4UOcHwsMpXlFOQDBvz2CJif+C7FDdEFwLdwH9GxF6j1LMjxY2vL5fbfgN48RgxBcWNumPK/Go5sD4zLwU+BJxfXp8Pr9rtlRQ3OnejaAgc6QnA/wP2BN4HfDMiFo55Yh7wtPLnHmWdPx0R60KKc/AJinP0L8B3y5u5w04CXk2RP+xIcZ7Gsk8Z4/4UjaGfK88HFDcfD6HIfw4ut3nviH0XUtx8fv0Y5T+BohF1T+AjwL+VQ5Bh/DzxbODszNydIof6+ohyJ5O7jnYuO4CzyjofRXFD+8wxYpdmxB5wkr4dEcNzIHw/M19QvTIiDgSOAZaWvdEui4gLp1jHTzPz2+XrTRHxBuC0zLylrONM4KaIeGWN5mP4RdmLqB/4AuWwiMz8RtU250fE31E0pFwwiTJfB3w2M/+7fH9uRLyLolFm+K7vJ4YbvSLiIopkZUJRTKj8YuCxZUPadRFxLg8kDSP9GXh0RFydmQMUSfJ4Rn4eo23zD2Xd10bEF4ETgRnNRRYRB1A0aD03M+8DrirvDr+S4s46wGWZeXG5/ZcpeqiN55/Ku/H9EfHxMs4vUDSMnlU2phIRHwLeFREHDfeCK9dvcyd/hH/IzE3A1RFxNUVj268pGtFOz8w7yvL/gaIR9j0Un8e+wEGZ+X8UjX/VPpWZN5f7fZDiTvE2jXCjycx/H35d/r8ZiIgFmXnXZPaXJDWlpsnVIuIwigaZ50+h7hcBm4E1FCMFdgCeA3xrEsd2MnDxcO4A9EXElcCzKXrNVXsiRQ+2j2fmEPAfEfG2MWLaQtHD69ER8fvMXD+J4zinulf6KPnVHVV1nx8RZ5TH+eVJlD2e5wD/m5nD5Xw1It4CnEDRIx/gi5n5mzKur1Pc8BzPezJzM/CDiPgu8LIoHlbxOuCw4fypzK++Avxdud9W4H3lvmO5sapX27kUvdr2joh5jJ8n/hk4OCL2LHvv/WxEudPKXctc7f/Kt7+PiH+haCCVZp0NcJJekONP8r4f0J+Z91Ytu5ni7tBk3Tzi/UEUSdXWqmVbKHoL/XYK5U7WkeXF9UEi4lUU3euXlIt2pbgbNxkHAadExOlVy3akOF/Dquesu3fEuvHsRfH7ufq8jXYnddiLKRpwPhwR1wDvHHl3dISRn8dE29wIPHYS+0xk+Lt094iyj656P/KczY+IHcZJ9kfGOXyODwLOjoiPVq3voLhTe+Mo+45lZDzDE17vx4M/k+q6/5nizumaMvn+XGZ+eBIxj6tsmP0g8FKK78jw/589KR4wIklqTU2Rq0XEwRQ9td+axQMNJusU4OvltX4wIr5ZLvsWEx/bQcBLI+KEqvXzgO+NUs9+wG/LBrBho+ZXWUzx8dcU1/PHREQv8LYxRi1UxzWe0eqebG44npE5yXDZ+1e9HyufGc1A2ZBVXdZ+FLnHLsDPqxoXOygaTYf9vmw8G8/9sWTmvWVZu1L03hsvT/wriulNro+IGyga3L5Tte20cteIeChF78GnUvRe3I6Jb2ZL02IDnKSJ3AYsjIhdqpKf6oTujxQXY+D+RoKR3f6HRry/GXhNZv54toOdrIg4iGK+h6UUd323RMRVFIkEbBvzSDcDH8xyTrApmqjs3wODFOf5+nLZgWNtnJlXAM8v7xyeRtEl/4Bx6pmofkapezjhfNDnTTHUYLJl30rxXdqtKrk6kJk1uh7AA8NcquMc/nzGG6I5mfMwllspkv5t6i6P7QzgjIh4DPC9iLhieJ4XHvz/pzrmieI7iaJHwTOA9RRDpwd44DsrSWpPdc/Vyrzqv4D3V/XEmsx+iynmJnt8RAwPB92F4gbcnkx8bDcDX87MyTxg6zZg/4joqGoIO5BiSpJtZOZXgK9ExO4Uvdz/iaI31nTzq9HqHu7NN15+NVG5wzlJtQMp5h+ejs6IeEhVI9yBwHXAnRRzrz0mM8fK3WaaW42ZJ2bm/wInRvGAixdR9GCsHmY7Vu46UXxnlcsPy8wN5TQmE02BIk2Lc8BJGlc5XO9K4MyI2DEinkTRpX3YbyiSpOeUDUDvpuiyP57PAB8skzUiYq+IGHOoQkTsFBHzy7c7RsT8qrkipushFBfb35d1vJry4Qyl3wGLR8xBUu3zwBsj4glRPO3qIeU52G0Sdf8OWBQRC0ZbmZlbgG9SnPNdopi77pTRti0/kxXlMMQ/UzyEYfhhFePWM4H3lHU/hmLOkPPL5VcBz46IhRGxD9sOEf0d8PAxjutm4CfAWeVneBjF3cyZzGP29ojoLIe3vrUqzs8Af1fGT0QsiIiXzqCekb4KvLv87u5JMdzmvLKu50bEweV3dPjzqH6AyJsjYnEUc7a8qyrmkUaey90ohuhsoEjSPzSLxyNJalL1ztWimIdtHfDpzPzMFMN/ZRlfUEzVcQTFHGO3ACdO4tjOA06IiOURsX2ZXxxbNuyN9FOKG5xviYgdIuJFFFOPjHZMERHHRTG3730UDU/V+dWSmPqTTh9a1j2vzEkeRfEgCCjyq1eU646mmP9s2O8per2Pml+VZRwSESeVx/VyivltvzPG9pPxD+X5firwXOAbWTzk6vPAx8peY0TE/hGxfAb13G+iPDEiTo6Ivco4hh+QUJ1fjZW7VhvtXO5G8bCJjeV3+e2zcTzSaGyAkzQZK4AnUfzh/wGKC9pmgHLuqTdRzLv1W4o7eLdMUN7ZFHf81kTE3RRzODxhnO2TIvHZn+LJRJso7/RFxLsi4pJx9h29wMz/AT5KkYz9jqKbevVd3nUUvZtuj4htnhKVxZM5X0dxh2yAYu6IUydZ9/UUDTj/LyI2RvHEp5FOo+iOfzvF/B1fHKfIVwLrI+IPFA+POHkK9YzlBxTHtJbiSVJryuVfpnj65nqKuVpGJjdnUTRMbYyI0Sb4PZFiyO+tFENL3peZfVOIa6QLgJ9TJK7fpXggBZn5LYo71V8rz8t1FA+lmC0foPiD4BrgWooHInygXPcIil4A91B8v/41H/w0tq9QnLv/V/77AKMbeS6/RDGk4rcUT90dOfeJJKl91TNXey1Fg8b7IuKe4X/DKyfI1U6huE7eXv2PogFw+ObjeMd2M0Xv8HdRNK7cTNGAss3fueUcci+iyNcGgJdT3PAczU7Ahyl6fd1O0Xj2rnLd8BzCGyLiF2PsP5r/psgR7qSYUuIlmbmhXPceigcLDFA8GOwrVXHfW27/4zIneOKI49pA0Uh2BsU5egfFPGrb5K+TdHsZx60UjV9vLHNKgL+lyA9/VuZX/0XReDpbxssTnwn8qvxunQ28YsRw17Fy1/uNcS7/geJhXndR5JJjfSekGesYGppJL1FJ7Sgizgeuz0wnKFXdRMQQ8IjR5vdrVBGxHnjtBHP5SJI0I62cq7XysdVbRBwLnJeZo/UgbEgRsQS4AZhXo4e5SbPGOeAkTSgijqF4gugNFI/2fj7FXUFJkiTVWSvnaq18bJLaiw1wkiZjH4ru2IsohiyszMxf1jckSZIklVo5V2vlY5PURhyCKkmSJEmSJNWQD2GQJEmSJEmSasghqG2ofJz2McBtPPjRzZIkqTVsD+wLXJGZm+sdjMy/JElqE2PmYDbAtadjgB/VOwhJklRzTwUuq3cQAsy/JElqJ9vkYDbAtafbAHp6ethnn33qHYskSZplt99+OytWrIDymq+GYP6lmjnvvPO46aab6h2Gpuiuu+4CYMGCBXWOpP4OPPBATj755HqHIc3YeDmYDXDtaQvAjv/1M3ZasEe9Y5EkqS3ttXJO/tBwqGPj2AKwzz77sHjx4nrHohZz5513cuuNN3LQgoX1DkVT8MeNGwHYc6i9p2a/8a5+dtllF383qtVsk4PZACdJkiRJTe6gBQt591OX1zsMTcEHftQL0Paf2/B5kFpdeze1S5IkSZIkSTVmA5wkSZLUhvr6+ujr66t3GJIkzZl6XvscgipJkiS1od7eYtjXsmXL6hyJJElzo57XPhvgGkRE7A18DHgiMAD8CfhIZn6rroFJkiS1MHMwSZI0FxyC2gAiogP4NvDDzHx4Zh4FvALwMTCSJEk1Yg4mSZLmij3gGsNxwJ8y8zPDCzLzRuCTEXEqcHRmngYQEd8BVgFdwKGZ+Tfl8tcBjwLeA3ydInHcHnh/Zp4/h8ciSZLULMzBJEnSnLABrjE8BvjFFPf5GnBNRLwjM/8MvBp4A/BM4NbMfA5ARCyY1UglSZJaR1vnYAMDA/T399Pd3V3vUDRDlUqF3Tu2r3cY0rRsvG8Tf6hU/F2kOVGpVFi4cGFd6nYIagOKiE9HxNURccVY22TmH4F1wHMj4pHAvMy8FrgWeEZE/FNEPDUz75qjsCVJkpqaOZgkSaoVe8A1hl8BLx5+k5lvjog9gSuBQR7cUDq/6vUXgHcB1wNfLPf9TUQcBTwbOCsi1mTmP9Y4fkmSpGbU1jlYZ2cnnZ2drFq1qt6haIa6u7v5862/q3cY0rTsMX9n9tpvb38XaU7Us6elPeAawzpgfkSsrFq2S/lzPXBERGwXEQcAjx/eIDP/GzgAOAn4KkBE7Afcm5nnUcxTcmTtw5ckSWpK5mCSJGlO2AOuAWTmUES8APhYRLwD+D3wR+BvgR8DN1AMa7iObecp+TpwRGYOlO8fC/xzRGwF/gysRJIkSdswB5MkSXPFBrgGkZm3UTz2fjQrxtn1KcDHqsrpBXpnMTRJkqSW1c452PLly+sdgiRJc6qe1z4b4JpUROwBXA5cnZlr6x2PJElSO2ilHGzZsmX1DkGSpDlVz2ufDXBNKjM3AofUOw5JkqR2Yg4mSZKmwwY4SZIkSWpyN97Vzwd+1FSjoNvejRv7Adr+c7vxrn4O3m/veoch1ZwNcG1s4ckvYK/Fi+sdhiRJkqQZ6OrqqncImoZFO+8IwLzOzjpHUl8H77e332G1BRvgJEmSJKmJrVzpQ3clqdFtV+8AJEmSJEmSpFZmA5wkSZIkSZJUQw5BlSRJkqQaWL16NZVKpd5hqDQwMABAZ5vPuTbXurq6HCYtYQNcW7vzvC8wb8Fu9Q5DUovYe+UZ9Q5BkqSGUqlU+L9f/w8HmnM3hDs33g3AQzbdU+dI2sdNd91d7xCkhmEDnCRJkiTVyIELduPvn/r4eoch4IM/uhzAz2MODZ9zSc4BJ0mSJEmSJNWUDXCSJElSm+rr66Ovr6/eYUiS2lC7XYMcgipJkiS1qd7eXgCWLVtW50gkSe2m3a5B9oBrMhFxz4j3p0bEp8rXZ0ZE94j16yNiz7mMUZIkqdWYg0mSpJmwAU6SJEmSJEmqIYegSpIkSW1qYGCA/v5+uru7J95YU1apVNi9Y6jeYUh1c9d9m/lDpeLvGI2qUqmwcOHCeocxZ2yAaz47R8RVVe8XAhdWvf+biDi56v1+cxOWJElSSzMHkyRJ02YDXPPZlJlHDL+JiFOBo6vWfywzV1WtXz9nkUmSJLWulszBOjs76ezsZNWqVRNvrCnr7u7mT7feXO8wpLpZMH8n9trvAH/HaFTt1jPSOeAkSZIkSZKkGrIBTpIkSZIkSaohG+AkSZIkSZKkGnIOuCaTmbuOeH8OcE75+sxRtl8yB2FJkiS1tFbNwZYvX17vECRJbardrkE2wEmSJEltatmyZfUOQZLUptrtGuQQVEmSJEmSJKmG7AEnSZIkSTVy011388EfXV7vMATcuPFuAD+POXTTXXdz8H71jkJqDDbAtbE9T34tey9eXO8wJEmSpJbU1dVV7xBUZc+dBwDYsbOzzpG0j4P38/+BNMwGOEmSJEmqgZUrV9Y7BElSg3AOOEmSJEmSJKmG7AEnSZIkqe2sXr2aSqVS7zCmbGCgGEbZ6TDKOdHV1WVPRkmzwga4NnbbuX/N0O471TsMjeKA03vqHYIkSVJLq1Qq/O//XMuBC5rrT6I77xoEYJd7b69zJK3vpvJcS9JsaK6rjSRJkiTNkgMX7MA7nrxHvcOYko/8ZCNA08XdjIbPtSTNBueAkyRJkiRJkmrIBjhJkiRJkiSphmyAkyRJktpQX18ffX199Q5DktSgvE7MLhvgJEmSpDbU29tLb29vvcOQJDWoVr1ObNiwgTPOOIP+/v45rdeHMDSJiLgnM3ctXz8G+CSwGOgAvgR8IDOHIuJU4J+B3wI7Ah/LzM/XJ2pJkqTmZg4mSVJr6enp4brrrqOnp4fTTz99zuq1B1yTiYidgQuBD2fmIcDhwJOBN1Vtdn5mHgEcC3woIvae80AlSZJaiDmYJEnNb8OGDaxZs4ahoSF6e3vntBecPeCaz0nAjzNzDUBm3hsRpwHfBz5dvWFm3hERFeAg4HdzHagkSVILabkcbGBggP7+frq7u+sdSl1UKhV269ha7zDUwO7avJW7K5W2/T8iVSoVFi5cWO8wZlVPTw9btxa/+7du3TqnveDsAdd8HgP8vHpBZlaAXSNi9+rlEfFw4OHA/81deJIkSS3JHEySpCa3bt06BgcHARgcHGTt2rVzVrc94JpPBzA0xrrh5S+PiKcAm4E3ZObcziwoSZLUelouB+vs7KSzs5NVq1bVO5S66O7uZvNvf13vMNTAFuy0HQ/dv6tt/49Irdj787jjjuPSSy9lcHCQHXbYgaVLl85Z3TbANZ9fAU+rXlDeZb0nM++OCCjmHzmtHsFJkiS1KHMwSZKa3IoVK1izZg0A2223HStWrJizuh2C2nx6gKdExDPg/gmBPwF8pK5RSZIktTZzMEmSmtyiRYs4/vjj6ejoYPny5XM6x5094JpMZm6KiOcDn4yITwPbA18GPlXfyCRJklpXK+Zgy5cvr3cIkqQG1qrXiRUrVnDjjTfOae83sAGuaWTmrlWvr6V4vP1o250DnDMnQUmSJLW4Vs7Bli1bVu8QJEkNrFWvE4sWLeKjH/3onNfrEFRJkiRJkiSphmyAkyRJkiRJkmrIIaiSJEmS2tJNdw3ykZ9srHcYU3LTXYMATRd3M7rprkEesX+9o5DUKmyAa2P7nvJxFi9eXO8wJEmSpDnX1dVV7xCmZc9dBgDYqbOzzpG0vkfs37zfE0mNxwY4SZIkSW1n5cqV9Q5BktRGnANOkiRJkiRJqiF7wEmSJEkCYPXq1VQqlXqH8SADA8WQy85ZHHLZ1dVlDzhJ0pyyAa6N/arntdyxYN6slHXkGy+alXIkSZJUP5VKhd/8+hr2XdBR71Du9/uNQwDssOmWWSnvtruGZqUcSZKmwgY4SZIkSffbd0EHb3jajvUO436f/eGfAGYtpuHyJEmaS84BJ0mSJEmSJNWQDXCSJElSE+vr66Ovr6/eYagN+F2TpOlzCKokSZLUxHp7ewFYtmxZnSNRq/O7JknTZw84SZIkSZIkqYZsgCtFxAsjYigiHjkLZV0QET+dYJszI6J7pnVJkiQ1M3MwSZLUDmyAe8CJwGXAK2ZSSETsARwJ7BERDxtjmxkN/Y2I7WeyvyRJUgMxB5MkSS3POeCAiNgV+Avg6cCFwJkRcSxwJnAncCjwc+DkzByKiPXAucAJwDzgpZl5fVnci4GLgN9RJJJnlXWcA/QDjwN+AdwNHB4R64ADgI9k5ucjogP4CPAsYAj4QGaeX8bzPuA24Ajg0WMcy7fL8uYDZ2fm52Z6fiRJkmrBHGx2DAwM0N/fT3f3zDv2VSoVdu4YmoWoGtfd9w1xR6UyK+er3VQqFRYuXFjvMCSpKdkDrvAC4NLM/A3QHxFHlssfB/w1RaL1cIoEcdidmXkksBqovnqfCHy1/HfiiHoOAZ6RmWeU7w8DngM8CXhvROwHvIgiuTsceAbwzxGxb7n944G/z8xRE7/SazLzKOBo4C0RsWgyJ0CSJKkOzMEkSVJbsAdc4UTg4+Xrr5Xvvwtcnpm3AETEVcASiiESAN8sf/6cImEjIvYGDgYuK+/SDkbEoZl5XbntNzJzS1W9F2TmJmBTRHyPIrl7CvDVcrvfRcQPgGOAP5Tx3DDBsbwlIl5Yvj4AeASwYQrnQpIkaa6Yg82Czs5OOjs7WbVq1YzL6u7u5u5br52FqBrXbvM72G+/rlk5X+3GXoOSNH1t3wBX3p08Djg0IoaA7SmGHVwMbK7adAsPPl+bR1n+cqATuCEiAHanGALx7nL9H0dUP7J//xDQMU64I/cfeSzHUtyxfVJm3hsR36cYBiFJktRQzMEkSVI7cQgqvAT4UmYelJlLMvMA4AaKu6BTdSLwzLKcJcBRjD+h8PMjYn6ZgB4LXAH8EHh5RGwfEXsBTwMun2T9C4CBMvF7JPDEaRyDJEnSXDAHkyRJbcMGuCJh+9aIZf8JnDSVQiJiCXAg8LPhZeVQhT9ExBPG2O1yimEWPwPen5m3lrFcA1wNrAPekZm3TzKMS4EdIuIa4P3VsUiSJDUYc7BZsnz5cpYvXz6XVapN+V2TpOnrGBpq7accaVtlonrDJ17zMPZaMG9WyjzyjRfNSjmSJGnmbrnlFpYuXQrwsMxcX+dwxAP519q1a1m8eHG9wxnT8Bxwb3jajvUO5X6f/eGfAGYtps/+8E/stt9jnQNOkjTrxsvB7AEnSZIkSZIk1VDbP4ShGZXzlawdZdXSzPSJp5IkSTVgDiZJkqbLBrgmVCZ4R9Q7DkmSpHbSLjnYbXcN3T/ssxHctrGYMme2YrrtriF2229WipIkadJsgGtjj1nxhYaeg0SSJElzq6urq94hbGNw5wEAduvsnJXydtuvMY9TktTabICTJEmSBMDKlSvrHYIkSS3JhzBIkiRJkiRJNWQPOEmSJKnNrV69mkqlUu8wJLWZgYFiiHnnLA0xr7euri57EmtMNsC1sR+f/2oWLZhX7zBmzTNee3G9Q5AkSWpKlUqF/PU1PHSPekciqZ3csbF8cd/NdY1jNtx/LNIYbICTJEmSxEP3gBOf7p+x9iXoAAAgAElEQVQHkubOV783CLTG757hY5HG4hxwkiRJkiRJUg3ZACdJkiS1qL6+Pvr6+uodhiSpxXm9mVjz9/OUJEmSNKre3l4Ali1bVudIJEmtzOvNxOwBJ0mSJEmSJNWQDXCliHhhRAxFxCNnoawLIuKnE2xzZkR0z7SuEWUuiYjrZrNMSZKkWjIHkyRJ7cAGuAecCFwGvGImhUTEHsCRwB4R8bAxtpnR0N+I2H4m+0uSJDUQczBJktTynAMOiIhdgb8Ang5cCJwZEccCZwJ3AocCPwdOzsyhiFgPnAucAMwDXpqZ15fFvRi4CPgdRSJ5VlnHOUA/8DjgF8DdwOERsQ44APhIZn4+IjqAjwDPAoaAD2Tm+WU87wNuA44AHj3G4ewQEeeW9fwGeFVm3jujEyRJklQD5mC1NzAwQH9/P93d43f6q1Qq7OSteUmatj/eB/2VyoS/b1tVpVJh4cKF9Q6joXmZLbwAuDQzfwP0R8SR5fLHAX9NkWg9nCJBHHZnZh4JrAaq/4edCHy1/HfiiHoOAZ6RmWeU7w8DngM8CXhvROwHvIgiuTsceAbwzxGxb7n944G/z8yxEj+AAD6XmYcBfwDeNInjlyRJqgdzMEmS1BbsAVc4Efh4+fpr5fvvApdn5i0AEXEVsIRiiATAN8ufP6dI2IiIvYGDgcvKu7SDEXFoZg7PCfKNzNxSVe8FmbkJ2BQR36NI7p4CfLXc7ncR8QPgGIpE7vLMvGGCY7k5M39cvj4PeAuwagrnQpIkaa6Yg9VYZ2cnnZ2drFo1fijd3d0M3HbNHEUlSa3nIfNh8b5dE/6+bVXt2vNvKtq+AS4iFgHHAYdGxBCwPcWwg4uBzVWbbuHB52vzKMtfDnQCN0QEwO4UQyDeXa7/44jqh0Z53zFOuCP3H81oZUqSJDUUczBJktROHIIKLwG+lJkHZeaSzDwAuIHiLuhUnQg8syxnCXAU408o/PyImF8moMcCVwA/BF4eEdtHxF7A04DLpxDDgRHxpKp4LhtvY0mSpDoxB5MkSW3DBrgiQfrWiGX/CZw0lUIiYglwIPCz4WXlUIU/RMQTxtjtcophFj8D3p+Zt5axXANcDawD3pGZt08hlF8Dp0TENcBCivlRJEmSGo052BxYvnw5y5cvr3cYkqQW5/VmYh1DQ/aObzdlonrDP77+QBYtmFfvcGbNM157cb1DkCSpIdxyyy0sXboU4GGZub7O4YgH8q+1a9eyePHieoezjeE54E58etvPUCNpDn31e4MALfG756vfG6Rz38Padg44FcbLwewBJ0mSJEmSJNVQ8zczt6FyvpK1o6xampkb5joeSZKkdtDqOdgdGx/ojSJJc+GOjcXPVvjdc8dG6Ny33lGokdkA14TKBO+ImZbzFy//YkMOgZAkSWpEs5WDNaKurq56hyCpHc0fAKCzs7POgcxc577+LtX4bICTJEmS2tzKlSvrHYIkSS3NOeAkSZIkSZKkGrIBTpIkSZIkSaohh6BKkiRJqonVq1dTqVTmtM6BgdrNKdXV1eVwXUnStNgA18Yu+cYpdO4xr95hTOglr7603iFIkiRpGiqVCtf/+hoWzeH86huK9jcGN99ck3IlSZoOG+AkSZIk1cyiTjhh6dzNfHPR2q3A7Nc5XK4kSdPhHHCSJEmSJElSDdkAJ0mSJGlS+vr66Ovrq3cYTc1zKEntyQY4SZIkSZPS29tLb29vvcNoap7D+tmwYQNnnHEG/f39DRlHLeJrlGOWZANczUXECyNiKCIeWb5fEhEnVa0/NiK+M4Py3zUbcUqSJLUSczBJI/X09HDdddfR09PTkHHUIr5GOWZJNsDNhROBy4BXlO+XACeNufXUmfxJkiRtyxxM0v02bNjAmjVrGBoaore3t249wsaKoxbxNcoxSyr4FNQaiohdgb8Ang5cCJwJfBh4VERcBZwL/LJq+8cDHwd2BjYBr87MjIhTgecBuwBdwLcy8x0R8WFg57KsX2Xmioj4NnAAMB84OzM/NycHK0mS1CDMwWpnYGCA/v5+uru7J7V9pVJhhxa55X/vpuJ4JnvsY6lUKixcuHCWotJk9fT0sHVr8STbrVu30tPTw+mnn94wcdQivkY5ZkmFFrkcNqwXAJdm5m+A/og4Engn8KPMPCIzPzZi++uBp2Xm44D3Ah+qWncE8HLgscDLI+KAzHwnsKksa0W53Wsy8yjgaOAtEbGodocnSZLUkMzBJD3IunXrGBwcBGBwcJC1a9c2VBy1iK9RjllSwR5wtXUixd1UgK+V7787zvYLgHMj4hHAEDCvat3azLwLICL+BzgIuHmUMt4SES8sXx8APALYMO0jkCRJaj7mYDXS2dlJZ2cnq1atmtT23d3d/P72a2oc1dzYZWfYa5+uSR/7WGbag07Tc9xxx3HppZcyODjIDjvswNKlSxsqjlrE1yjHLKlgD7gaKe96Hgd8ISLWA2+nuHvaMc5u7we+l5mHAidQDGEYtrnq9RZGaTyNiGOBZwBPyszDKYZWzB+5nSRJUqsyB5M0mhUrVrDddsWfv9tttx0rVqyYYI+5jaMW8TXKMUsq2ABXOy8BvpSZB2Xmksw8ALgB2ArsNsY+C4Dflq9PnWQ9f46I4bu0C4CBzLy3fOLXE6cXuiRJUtMyB5O0jUWLFnH88cfT0dHB8uXL6zYP31hx1CK+RjlmSQWHoNbOiRST/Vb7T4oncQ1GxNXAOVRNAAx8hGL4w9uAdZOs53PANRHxC+A1wBsj4hoggZ9NP3xJkqSmZA4maVQrVqzgxhtvrHtPsLHiqEV8jXLMkqBjaGio3jFojkXEEuCG7jcupnOPeRNtXncvefWl9Q5BkqSmcssttwzP9fOwzFxf53DEA/nX2rVrWbx4cb3Dmba+vj4Ali1bNqnth+eAO2Hp3A28uWht8dTH2a7zorVb2Wufw2Y8B9xUz6EkqXmMl4PZA06SJEnSpNhoNHOeQ0lqT84BJ0mSJEmSJNWQPeAkSZIk1cyGgQeGhc5VfTD7dW4YgL32mdUiJUltxAa4Nvasl57b1HOQSJIkqbF1dXXNeZ077FS0wHV2ds5quXvtU5/jkSS1BhvgJEmSJNXEypUr6x2CJEkNwTngJEmSJEmSpBqyB5wkSZKkhrN69Woqlcqc1DUwMDvDVru6uuz1J0kalQ1wbezr33wVC/aYV+8w+KtX9dY7BEmSJDWYSqXCr399DXvM7lRuo9pYPrjhvs03z7gMSZJGYwOcJEmSpIa0Ryccd3xHzetZt2YImFldw2VIkjQa54CTJEmSJEmSasgGOEmSJEmSJKmGbICTJEmS2lBfXx99fX31DkMNyO+GJM0+54CTJEmS2lBvb/EgrGXLltU5EjUavxuSNPtsgJuCiNgCXAt0AFuA0zLzJxPsc09m7hoRS4AnZ+ZXyuVHA6/KzLfUOGxJkqSmZf4lSZJagUNQp2ZTZh6RmYcDfwecNYV9lwAnDb/JzCtrkfxFxPazXaYkSVIdmX9JkqSmZw+46dsdGACIiF2BC4BOYB7w7sy8YMT2HwYeFRFXAecCvwS6M/O5EXEm8DBgX+AQ4G3AE4FnAb8FTsjMP0fEUmAVxed2BbAyMzdHxHrg34HjgU9FxBXAp4G9gHuB12Xm9bU4CZIkSXPI/GsWDQwM0N/fT3d3d71DGVWlUmG7JuoucN+mIuZGPZ9TUalUWLhwYb3DkKSW0kSXtIawc0RcFRHXA18A3l8uvw94YWYeCTwd+GhEdIzY953Aj8o7uB8bpewu4DnA84HzgO9l5mOBTcBzImI+cA7w8nL5DsDKqv3vy8ynZObXgM8Bp2fmUUA38K8zPnJJkqT6MP+SJElNzx5wU7MpM48AiIgnAV+KiEMp5iT5UEQ8DdgK7A/sDdw+hbIvKe+yXgtsD1xaLr+WYvhEADdk5m/K5ecCbwY+Xr4/v4xrV+DJwDciYrjsnaZ4nJIkSY3C/KtGOjs76ezsZNWqVfUOZVTd3d3cdvs19Q5j0ubvDPvu09Ww53MqWqEXnyQ1GhvgpikzfxoRe1IMM3h2+fOoMolbD8yfYpGby3K3RsSfM3OoXL6V4nMaeUd3pD+WP7cDNg4nqpIkSa3C/EuSJDUrh6BOU0Q8kuJO6QZgAXBHmfw9HTholF3uBnabQZXXA0si4uDy/SuBH4zcKDP/ANwQES8t4+yIiMNnUK8kSVJDMP+SJEnNyh5wU7NzOYkvFHdET8nMLRHRA1wUEVcCV1EkayNdAwxGxNUUc4n8cioVZ+Z9EfFqiqENw5MAf2aMzVcAqyPi3RSTEn8NuHoq9UmSJDUI868aWb58eb1DUIPyuyFJs69jaGho4q3UUiJiCXDDG968Pwv2mFfvcPirV/XWOwRJklrKLbfcwtKlSwEelpnr6xyOeCD/Wrt2LYsXL653OE1heA64446faCTwzK1bU/xNNJO61q0ZYt99DmuJOeAkSdMzXg7mEFRJkiRJkiSphmyAkyRJkiRJkmrIOeAkSZIkNaSNAw8MD611PTCzujYOwL77zFJAkqSWYwNcG3vZi77kHCSSJElqSF1dXXNW1/ydiha4zs7OaZex7z5zG7MkqbnYACdJkiSp4axcubLeIUiSNGucA06SJEmSJEmqIXvASZIkSXqQ1atXU6lUABgYmPzwzK6uLnuuSZI0Chvg2tjnL3wVu3XOm9M6u0/sndP6JEmSNHWVSoX/uf4adlsId28olt3z55vH3efu/jkITJKkJmUDnCRJkqRt7LYQjnlWB1dcUjwZ9JhndYy7/fB2kiRpW84BJ0mSJEmSJNWQDXCSJElSi+rr66Ovr29O6tq86YH54iRJ0oM5BFWSJElqUb29xfy7y5Ytq3ldf7oX+gedCE6SpNHYA66GImJLRFwVEVdHxC8i4sn1jkmSJKnVmYO1tw0bNnDGGWfQ379tY+B46xrddGJv5uOVpFZjA1xtbcrMIzLzcODvgLPqEURE2NNRkiS1E3OwNtbT08N1111HT0/PlNY1uunE3szHK0mtxqRg7uwODABExK7ABUAnMA94d2ZeEBFLgEuAy4AnA78Fnp+ZmyLiGODfgD+W65+VmYdGxHxgNXA0MAi8LTO/FxGnAs8B5gMPAY6bqwOVJElqIOZgbWTDhg2sWbOGoaEhent7WbFiBQsXLpxwXaObTuzNfLyS1IpsgKutnSPiKooEbF8eSMDuA16YmX+IiD2Bn0XEheW6RwAnZubrIuLrwIuB84AvAq/PzJ9ExIer6ngzQGY+NiIeCayJiEPKdU8CDstM+5xLkqR2Yg5WGhgYoL+/n+7u7intV6lUGGrCvxR6enrYunUrAFu3bqWnp4fTTz99wnWNbjqxN/PxSlIrcghqbQ0Pf3gk8EzgSxHRAXQAH4qIa4D/AvYH9i73uSEzrypf/xxYEhF7ALtl5k/K5V+pquMpwJcBMvN64EZgOPnra4TET5IkaY6Zg7WpdevWMTg4CMDg4CBr166d1LpGN53Ym/l4JakVNeF9reaUmT8t77TuBTy7/HlUZv45ItZT3KEF2Fy12xZgZ4pkcSzjrfvjtAOWJElqAe2eg3V2dtLZ2cmqVaumtF93dzc333FNjaKqneOOO45LL72UwcFBdthhB5YuXTqpdY1uOrE38/FKUiuyB9wcKYcmbA9sABYAd5SJ39OBg8bbNzMHgLsj4onloldUrf4hsKKs4xDgQCBnOXxJkqSmZA7WXlasWMF22xV/4my33XasWLFiUusa3XRib+bjlaRWZANcbe0cEVeVc5CcD5ySmVuAHuDoiLiSInG7fhJl/RXwuYj4KcUd17vK5f8KbB8R15Z1nJqZm8coQ5IkqR2Yg7WpRYsWcfzxx9PR0cHy5csf9NCB8dY1uunE3szHK0mtyCGoNZSZ24+x/E6KyXlHc2jVdtVjBX6VmYcBRMQ7gSvLbe4DTh2ljnOAc6YRtiRJUlMzB3vA8uXL56yuHXeBhbvXv5FnxYoV3HjjjaP2+BpvXaObTuzNfLyS1GpsgGsez4mIv6P4zG5klIRPkiRJs66pc7Bly5bNWV077VzMOVdvixYt4qMf/eiU1zW66cTezMcrSa3GBrgmkZnnUwxvkCRJ0hwxB5MkSbPBBjhJkiRJ27i7H664ZIi7NxTvr7hkaMLteWjt45IkqRnZANfGXve8L7F48eJ6hyFJkqQG09XVdf/rgXkDwCSGlz70wftJkqQH2AAnSZIk6UFWrlxZ7xAkSWop29U7AEmSJEmSJKmV2QAnSZIkSZIk1ZBDUCVJkqQ2tHr1aiqVyqjrBgbGnvetq6vLIaqSJE2RDXBt7B8vPYVdOufVOwxJmpKPv/jSeocgSS2hUqlw3fXXMH/RtuvuK598umHw5lGXS5KkqbEBTpIkSWpT8xfBkudtOyvN+gu3AtuuG14uSZKmxjngJEmSJEmSpBqyAU6SJElqYn19ffT19dU7jGlr9vglSZoMh6BKkiRJTay3txeAZcuW1TmS6Wn2+CVJmgx7wEmSJElqOBdddBHHH388F198cb1DkSRpxtqmB1xEbAGuBTqALcBpmfmTCfb5SWY+OSL2Az6RmS+JiKOBV2XmW2Ypru2Bn49YvBhYm5kvn406JEmS6sUcTNP1qU99CoCzzz6bZz/72XWORpKkmWmbBjhgU2YeARARy4GzgL8cb4fMfHL581bgJeXrK4ErZyuozNwCHDH8PiL2BS4H3j9bdUiSJNWROZim7KKLLmJoaAiAoaEhLr74YhvhJElNrZ0a4KrtDgwMv4mItwMvA3YCvpWZ7yuX35OZu0bEEuA7mXloRBwLdGfmcyPiTOBA4OHlz49n5ifKfd8DrABuBu4Efp6Zq8YLKiI6gHOBf87M68pl7wVOAHYGfgK8ITOHIuL7wH8DTwf2AP4qM38UEbsA5wCPBH4NLAHeXCatkiRJ9WQOVgMDAwP09/fT3d09pf0qlQqDU/xrYPDeYr+p1jVRHAsXLnzQsuHeb8PsBSdJanbtNAfczhFxVURcD3yB8u5mRBwPPAJ4PMVd0KMi4mlTKPeRwPJy//dFxLxyiMSLgccBLwKOnmRZfwMMAp+sWvapzDwmMw+lSACfW7Vuh8x8PPDXwPvKZW8CBjLzsPIYj5rCsUiSJM02czBN2XDvt7HeS5LUbNqpB1z18IcnAV+KiEOB48t/vyy325UiGfzhJMv9bmZuBjZHxB3A3sBTgAsyc1NZ30UTFRIRh1MkccdkZnWG8fSIeAewC7AQ+BUwXN43y58/p7jLSln32QCZeV1EXDPJ45AkSaoFc7Aa6+zspLOzk1Wrxu3ot43u7m7+7/dTC3OHXaBrr64p1zVRHCN1dHQ8qNGto6Nj1uqTJKke2qkB7n6Z+dOI2BPYi2JC4LMy87PTLG5z1estFOd0ShlCROwM9ABvyszfVS2fD/wrcHRm3lwOt5g/St3D9TLVuiVJkuaKOZgm67TTTuOTn3ygQ+Jb3/rWOkYjSdLMtdMQ1PtFxCOB7YENQC/wmojYtVy3f0Q8dIZVXAacEBHzy3KfM8H2q4AfZOZ3RiwfTvTuLMt5ySTrfhlARDwaeOzkw5YkSaodczBN1gknnHB/r7eOjg7nf5MkNb126gG3c0RcVb7uAE4pn361JiIeBfw0IgDuAU4G7hix/6QnnsjMKyLiQuBq4EaKJ3bdNdq2EbEfxZwh11fFB/CrzFwREZ8HrgXWA1dMovp/Bc4thz38ErhmrLolSZLmgDlYjS1fvnwuqqmZseIf7gVn7zdJUivocELTiUXEUcC/ZOZfTmGfXTPznvKJWD8EXp+Zv6hZkA/Uuz0wLzPvi4guYC1wSGb+qWqbJcANf/n2xezSOa/WIUnSrPr4iy+tdwhSw7vllltYunQpwMMyc32dw5m2VsrBhvOvtWvXsnjx4lqHMynDc8Ated62g2LWX7gVYJt16y/cysF7HTarc8BJktQqxsvB2qkH3LSUT9P6CvDOKe76uXL4wXzg3LlI/Eq7AN+LiHkUd5lXVje+SZIkNQNzMEmS1EpsgJtAZl4JHDKN/U4auSwiPg38xYjFZ2fmF6cZ3mj13g0cPVvlSZIk1YM52Ny4b8MDvd1GLodt1923geIRGpIkaUpsgJtDmfnmesdQ7b3PPLdhhkBIkiTVSqPlYI2iq6trzHUDOwwA0NnZ+eAVe42/nyRJGp0NcJIkSVIbWrlyZb1DkCSpbWw746okSZIkSZKkWWMDnCRJkiRJklRDDkGVJEmS1HBWr15NpVKZ8n4DA2PMXzeHurq6HOIrSXoQG+Da2Kl9f8u8hfPrHcaoLnn+v9U7BEmSJNVRpVLhmuuvhT13nNqOd/4JgJu33FGDqCZfvyRJ1WyAkyRJktSY9tyR7V/40CntsuVbRcPbVPebLcP1S5JUzTngJEmSJEmSpBqyAU6SJEnSjPX19dHX11fvMDRJfl6SNLdsgJMkSZI0Y729vfT29tY7DE1So3xeGzZs4IwzzqC/v7/eoUhSTdkAJ0mSJEmqi56eHq677jp6enrqHYok1VTDN8BFxJaIuCoiro6IX0TEk6dRxvqI2LN8/ZNpxvH9iLgpIjqqln07Iu6ZTnlTrHtJRFxXvj46Ij5R6zolSVJ7MwczB5NqbcOGDaxZs4ahoSF6e3vtBSeppTXDU1A3ZeYRABGxHDgL+MvJ7Fgmah3VyzJzysljlY3AXwCXRcQewL5T2Xk4nszcOt0AMvNK4Mrp7i9JkjRJ5mBVzMEmNjAwQH9/P93d3bNSXqVSgXlbZqWsOXXvFiqVyqydh1qpVCosXLiwrjH09PSwdWvx33Lr1q309PRw+umn1zUmSaqVZmiAq7Y7MDD8JiLeDrwM2An4Vma+LyKWAJcA3wOeBLyguoCIuCczd42IY4EzgTuBQ4GfAydn5tA49X8NeAVwGfAi4JvAY8pydwUuADqBecC7M/OC0eKJiGcAfwvcCvwvsDkzT4uIc4DvZOZ/VMc6Iv5jge7MfG5EPAT4JPBYis/yzLLOxwBfBHak6OX44sz833HPrCRJ0tjMwczBpFm3bt06BgcHARgcHGTt2rU2wElqWc3QALdzRFwFzKe423kcQEQcDzwCeDzFHdYLI+JpwE1AAK/OzDeV245V9uMokrdbgR9T3lkdJ5a1wOcjYnuKJPD1wHvKdfcBL8zMP5RDLX4WEReW6+6PJyL2K/c5ErgbWAdcPfnT8SB/D6zLzNeUd4Mvj4j/At4InJ2ZPRHx/9m79/i6qjrh/5+UlmsBT0sBpx0oBvwiFARFHxjFQUoJKlfBCxMFRPSxCqg/44OI84ijI6BFZcQp3hXtOI83hiJCCGWUi4JXLKB8hQMFiwo0iUiBgYbm98fZqadpkiZtTk6S83m/Xnmx99p7rfU9pyl8WXuttbcEttjE9iVJUuMyBxucOdgASqUSpVKJRYsWjUp7bW1tLF+Vo9LWmNp2C5p3ah6176FWxsMMvcMPP5xrr72Wnp4epk6dyvz58+sdkiTVzLjfA45i+UNm7g0cBVxeLCM4svj5NfArYG8qySDAA5l56zDa/llmriyWI9wOzN3I/c9SSQ7fAGyTmSuqrjUBH4+I5cD1wGxglwHieSnw48zsysw1wHeGEedgjgQ+UCTHP6KSIO8G/BT4YEScA+yemU9tRh+SJKkxmYMNzhxMGgWtra1MmVL5X9IpU6bQ2tpa54gkqXYmwgy4dTLzp8WTzVlUkq0LMvPz1fcUyw2eGGaTT1cdP8vwvo//BK6gsnSiWmsR14szc01ErKCSjNEvniYG10MxKFokuFtuJJYmKksb+j8a/F1E3Aa8BmiPiDMy84aNtCVJkjQgc7ANmINJo2DmzJkceeSRXH311bS0tNR9TzpJqqWJMANunYjYm8pU/k6gHTi92PeDiJgdETuPQRg3UdmE+Fv9yncEHikSv1cCuw9S/2fAP0ZEKSKmAidWXVsBvLg4Po7KPiZDaQfO6nsrWEQcWPzzecB9mflvwFJg/+F8MEmSpIGYg23AHEwaJa2trcybN8/Zb5ImvYkwA65v/xGoPG08NTOfBa6LiBcAPy32F1kNvInKU9SaKTYIHmhDhyXAVRHxCypLKe4epP5DEfFx4DYq+578FnisuPxF4MqI+BmVvU429hT5o8BngOVFArgCOJrK8ow3RcQa4M/Avwz7A0qSJFWYgw3OHGwALS0t9Q5BIzBe/rxmzpzJxRdfXO8wJKnmmnp7h3rhlGohIqZn5uri6esVwFcy84ox7H8ucP/ff+Agps3YemO318U1x3253iFIkjRhrVy5sm8z8z367ZfW0OqZg/XlX8uWLWPOnDlj0eWE1/cShi1OGNkEy2eveARgxPVGy7NXPML+O8W4fwmDJGn0DZWDTaglqJPI+cUT5TuB+4H/qnM8kiRJjcAcTJIk1cVEWII6piLiCmCPfsXnZGb7aPWRmfV/57ckSdI4Yg6mAa16Zt2MtpHUAUZeb7SsegZ2qk/XkqTxywG4fjLzhHrHMFa+tuAil0BIkqRxoZFyMA1Pc3PzJtXr3qIbgFKpNJrhDN9Omx67JGnycgBOkiRJ0rizcOHCeocgSdKocQ84SZIkSZIkqYacASdJkiRNUosXL6ZcLo+4Xnf34Ms4m5ubnZ0mSdIIOQDXwN7SfinTZkyvdxiSNCI/POFD9Q5BkiaMcrnM8rvvgp1GmPOtWg3AH57964DlkiRpZByAkyRJkiaznaYz9bgXjqhKz5W/AdigXl+5JEkaGfeAkyRJkiRJkmrIAThJkiRJkiSphhyAkyRJkiawjo4OOjo66h3GqJusn0uS1JgcgJMkSZImsPb2dtrb2+sdxqgb6nN1dnbyvve9j66urhG1ee+993L88cdz3333DatckqTRUtMBuIjYNSL+MyLKEfHbiPhhRDx/M9s8PyLahnnvTwYp/1pEnLQ5cfRr77CIeCwibo+I5RFxfUTsvIltjWpskiSpsTz66KN89KMf7TvtMP8aVkfDKEcAACAASURBVFvmXxPMkiVLuPPOO1myZMmI6l100UU8+eSTXHDBBcMqlyRptNRsAC4imoArgB9lZnNm7gN8ENilVn1W9b0FQGb+Q637qnJTZh6QmfsDPwfeNYZ9DykifNutJEkNoLe3lzPPPJMDDjigr2gB5l91Yf5VO52dnVx33XX09vbS3t4+7Flw9957Lw888AAADzzwwLrZboOVS5I0mmqZGLwSWJOZl/UVZObtETE9IpYBJWAa8KHMvDIi5gLXAjcDBwO/Ab4KfATYGWjNzJ8VTb0wIm4A/h74RGZ+MSIOAz4M/Ak4ANgnIlZn5vRiMPCzwOHA/UBTX0wR8WrgU8Aq4FfA8zLz6IjYrqizH5Xv6fzMvHJjH7roa3vg3uL8fGB1Zi4qzu8Ejs7MFRFxCtAG9ALLM/PN/dr6aPEZTwfuA74OHFN8b6/LzLsHizMiTgNeA2wNbFd8dkmSNIndeuutTJ06lWOOOYbPfOYzgPlXcT6p86/u7m66urpoa9twkmK5XIZpa0evsyefoVwuD9jXaCuXy8yYMWOD8iVLlrB2beUzrV27liVLlnDWWWdttL2LLrpovfMLLriAL37xi4OWS5I0mmq5BHUe8MsByv8HOCEzX0RlkO7iImkC2BO4BNgf2Bv4J+DlVJKkD1a1sT+V5OYQ4P9GxN8V5S8Fzitm21U7AQgqSdLbgH8AiIitgc8Dr8rMlwOzquqcB9yQmS8p4vxkkWwN5tCIuB14EDgC+MoQ9xIR+xZ9HJ6ZLwTe3e/6J6gkvm/JzL6saVXxvS2m8p1sLM5DgFMz08E3SZIawD333MO+++470CXzL8y/JosbbriBnp4eAHp6eli2bNmw6vXNcut/Pli5JEmjqR5T45uAj0fEK4C1wGz+tizi/sy8AyAi7gKWZWZvRNwBzK1q48rMfAp4KiL+m0ri9xfgZ5l5/wB9vgL4VmY+C/yxeHoLlSTzvqo63wLeXhwfCRxbtd/J1sBuwO8G+Vw3ZebRReznAJ8A3jHE93A48N3MXAWQmdVz5/8ZuC0z396vzveLf/4SeO1G4gTo6NeuJElqTOZfFZMy/yqVSpRKJRYtWrTBtba2NpavGsUBpW23pHmn3Qfsa7QNNsvu8MMP59prr6Wnp4epU6cyf/78YbW3++67rze4tvvuuw9ZLknSaKrlANxdwECb2bZSedL54sxcExErqCQtAE9X3be26nwt68fa26/NvvMnhoinfx2oWgoxyLUTMzOHuGcwS4HvFcc9rD/TsO+zNg0SE1T2MHlxRMzol8D1fR/P8rfvY8A4I+J/MfT3IUmSJpm99tprsLdGmn/9rX3zrwmutbWV6667DoApU6bQ2to6rHrnnHMO73znO9edn3vuuUOWS5I0mmq5BPUGYKuIeFtfQUS8BNgdeKRI/l5ZnI/UcRGxdUTMBA6jkjAN5UbgjRGxRUQ8l8pSAYC7gecV+58AvKGqTjtwVt/yjIg4cATxvRwoF8crgBcVbbwI2KMoXwa8vvgMRET1BhfXAhcCV0fE9hvpa3PilCRJk8jBBx/MM888w9VXX72uzPzL/GuymTlzJkceeSRNTU20tLQMuE/cQPbcc8/1Zr0973nPG7JckqTRVLMBuMzspbL3x4KIKBdLGs4HfggcFBG/oPI09u5NaP5nwNXArcBHM/OPG7n/CuAe4A4q+3f8uIjxKeCdwLURcTPwMPBYUeejVDbbXV5s3PvRjfRxaETcHhG/Ad4MvK8o/x4wo9ifZCHw+6Lvu4B/BX5c1PlUdWOZ+R3gi8DSiNhmiH5HGqckSZqkmpqauPTSS/nlL9dtw3sd5l+TPv9qaWmhpaWlniHUxFCfq7W1lXnz5g179lufc845h2233XaDWW6DlUuSNFqaensHm4XfGCJiemauLp5gfg64JzM/Xe+4aql44nz/bv/nBKbNmF7vcCRpRH54wofqHYI07q1cubJvX6w9MnNFncPZQCPnX8uWLWPOnDlj1m/fHnBTj3vhiOr1XPkbgA3q9Vz5G/Yfoz3gJEmaaIbKwWq5BHWieFvxdPQuYEcqb+WSJElS7Zh/SZKkhlKPt6COK8XT1mE9cY2IFuCifsX3Z+YJox6YJEnSJGX+JUmSGk3DD8CNRGa2U9l0V5IkSWPA/GsUrFq9bknpSOoAG9ZbtRp2GqW4JElqIA7ANbCvtpw5pnuQSJIkaWw1NzdvUr3uLboBKJVK61/YadPblCSpkTkAJ0mSJE1SCxcurHcIkiQJX8IgSZIkSZIk1ZQz4CRJkiSNW4sXL6ZcLtc7DI2B7u5Blj5r1DU3NztDVhpjDsA1sNPbL2daaYd6h6Fx4OrXnlnvECRJkgZULpdZfvdvaZq5Y71DUY31dj4GwMqeJ+scyeTW9z1LGlsOwEmSJEka15pm7sjU415W7zBUYz1X3gLgn3WN9X3PksaWe8BJkiRJkiRJNeQAnCRJkjRJdXR00NHRUe8wJEkaE+P5v3suQZUkSZImqfb2dgAWLFhQ50gkSaq98fzfPQfgxomIWJ2Z04vj5wOfAZ4PrAHuAM4CXgBcCdxPZfbiI8A/ZeYjEXEa8FXgiMxcVrRzAvB94HWZ+d2x/USSJEnjnzmYJEkaCy5BHWciYmvgamBxZu6ZmS8AFgOziltuyswDMnN/4OfAu6qq3wGcXHX+RuA3YxC2JEnShGYOJkmSaskZcOPPPwE/zcyr+goy878BIuKwvrKIaAK2B+6tqnsTcGhETAO2AvYEbh+DmCVJkiY6czBJklQzDsCNP/OAXw5x/dCIuB2YCTwBfLDqWi9wPdAC7AgsBfaoUZySJEmTyaTMwbq7u+nq6qKtra3eoWyycrlM79SmeochTRq9Tz5NuVye0P9ekAZTLpeZMWNGvcMYkEtQJ56+5Q9/T2W/kU/0u/6fVJY9vBH41lgHJ0mSNEmZg0mSpE3mDLjx5y7gH4d571Lge9UFmfmziJgHPJWZv4+I0Y5PkiRpMpqUOVipVKJUKrFo0aJ6h7LJ2trauOPRh+odhjRpNG27Fc2zZk/ofy9IgxnPMzudATf+/AfwDxHxmr6CiDgqIvYb4N6XA+UBys9l/WURkiRJGpo5mCRJqhlnwI0zmflURBwNfCYiPgOsAZYD76ay50jf/iNNwGPAGQO0cc0YhixJkjThmYNJkqRacgBunMjM6VXHdwNHDXDbw1Q29h2o/teArw1QftqoBChJkjQJTfYcrKWlpd4hSJI0Zsbzf/eGHICLiKuovNVpQJl57KhHJEmS1ODe8Y53DHn9sssuG6NINNEtWLCg3iFIkjRmxvN/9zY2A85dGSVJksbY6aefXu8QJEmSNIqGHIDLzB+PVSCSJEmqeOlLX1rvEKRxpbfzMXquvKXeYajGejsfA/DPusZ6Ox+DWbPrHYbUcIa1B1xE7AVcAOwDbN1XnpnPq1FcGgNfaTmFOXPm1DsMSZI0iBUrVvCpT32Ke++9l6effnpd+bJly+oYlTS2mpub6x2Cxkj31G0BKJVKdY5kkps1279XUh0M9yUMXwU+DHwaeCXwFipvgJIkSVKNnHvuuZx99tl8/OMf5/LLL+f73/8+vb2Dbs8rTUoLFy6sdwiSJG22KcO8b5vMXAY0ZeYDmXk+cHjtwpIkSdLTTz/NIYccAsDs2bM566yzuPXWW+sclSRJkkZquDPg/icipgD3RMSZwEPAzrULS5IkSVtuuSVr165l991355vf/Ca77LILnZ2d9Q5LkiRJIzTcAbj3ANsCZwMfpbIM9ZRaBSVJkiT44Ac/yFNPPcWHPvQhLrnkEm677TY+8YlP1DssNYDFixdTLpdHVKe7uxvY+P5dzc3NLiuVJDWc4Q7Azc3MnwOrqez/RkS8DritVoGp9t567XeZVnrOZrXxgxNPG51gJEnSBh566CH2339/tttuOy644AIArrnmGl74whfWOTJNduVymeV3/46mmcPfDL+3szIAt7Ln6Y3eI0lSoxnuANy5wHeGUSZJkqRR8oUvfIFXvepVGy2TaqFpZolpxx4x7PvXLL0eYMg6ffdIktRohhyAi4hXAa8GZkfEv1Vd2gHoqWVgkiRJjerHP/4xN954Iw8//DAf+9jH1pWvXr2aLbbYoo6RSZIkaVNsbAbcH4FfAMcCv6wqfxx4b62CkiRJamS77LIL8+bN44YbbmDfffddV77ddttx7rnn1jEyTSYdHR0ALFiwoM6RjNxEjl2S1JiGHIDLzN8Av4mIJZnpjDdJkqQxsPfee7P33ntzzDHHMHXqcHcMkUamvb0dmJiDWBM5dklSY9rYEtRvZ+brgV9HRG//65m5f80ikyRJalDvfve7ueSSSzjhhBMGvH7VVVeNcUSSJEnaHBt7pPru4p9H1zqQsRARqzNz+jDvPRbYJzMv3Jx2htnXJ4FjgGeAMvCWzPxLRMwEvgu8BPhaZp5ZVedk4INAL5Wlwm/KzFURMQP4f8BcYAXw+sz0dVOSJE0g5513HgCXXXZZnSMZHeZgkiSp0W1sCeqfin8+EBG7Ai+lkmz8PDP/PAbx1U1mLgWWjlF3HcC5mdkTERdRecPsOcD/AP8MzCt+AIiIqcAlVJLTVRHxCeBM4HzgA8CyzLwwIj5QnJ8zRp9DkiSNgp133hmA2bNn8+ijj7J8+XKamprYb7/9mDVrVp2jqy1zsLHT3d1NV1cXbW1tG1wrl8v0Tp0y6n32PvkU5XJ5wD5HolwuM2PGjFGKSpKk2hvWpiIRcQbwf4EbgCbgsxHxL5n5lVoGNxYiYhZwGbBbUfSezLwlIk4DDsrMMyNiD+A/qHxf1/ar/37g9cBWwBWZ+eGImAtcA9wM/APwEHBcZj41UAyZeV3V6a3ASUX5E8DNEbFnvypNxc92EdFJ5a209xbXjgMOK46/DvyIcZ78SZKkgX3nO9/hc5/7HAcffDC9vb187GMf453vfCcnnXRSvUPbbOZgkiSpkQx3V9/3AwdmZidAMS3/J8CEH4Cj8hTz05l5c0TsBrQDLxjgnsWZeXlEvKuvMCKOBPaiMjOwCVgaEa8AHizKT87Mt0XEt4ETgW8OI57TqSxfGFRmromIhcAdwBPAPUBfXLtUzVz8U0TsPIw+JUnSOPSlL32JK664glKpBFRmLL3xjW+cFANwmIPVXalUolQqsWjRog2utbW1ccejo7/gpWnbbWieteuAfY7E5s6gkyRprA13XvlK4PGq88eBP4x+OHVxBHBpRNxOZbnDDhGxfb97XgZ8qzj+RlX5kcXPr4FfAXtTSfoA7s/M24vjX1LZD2RIEXEe0AMs2ch904CFwIHA3wHLqSyZkCRJk8iuu+7Kdtttt+58u+2247nPfW4dIxpV5mCSJKlhDHcG3EPAbRFxJZU94I4DfhYR/x9AZn6qRvGNhSnAIf2XJkRE//s2eAsslSeuF2Tm5/vVnQs8XVX0LLDNUEFExKlUXnYxPzMH6qvaAQCZWS7qfpvKPiMAD0fEc4snr88FHtlIW5IkaZzaZZddeP3rX8/8+fNpampi2bJl7Lfffnz1q18F4C1veUudI9ws5mCSJKlhDHcGXBn4L/6WAF0J/AnYvviZyK6jsnkuABFxwAD33AK8sThurSpvB06PiOlF3dmbstwgIo6iskfIsZn55DCqPATsU+ydArAA+F1xvBQ4tTg+lcqflSRJmoB22203jjjiCJqamgCYP38+O++8M0888QRPPPFEnaPbbOZgkiSpYQxrBlxmfqTvOCKmANMz8681i6p2to2IlVXnnwLOBj4XEcupfB83Au/oV+/dwH9ExLuB7/UVZuZ1EfEC4KfF09rVwJuoPG0diUupbCDcUbRza2a+AyAiVlDZ4HfLiDgeODIzfxsRHwFujIg1wAPAaUVbFwLfjoi3UtkH5XUjjEWSJI0TZ565bnyKtWvX8uSTTzJ9+vQ6RrTJzMHGoZaWlnqHsMkmcuySpMbU1Nu7sZn2EBH/QSUhepbKXho7Ap/KzE/WNjzVQrE84/6571/ItNJzNqutH5x42miEJEmSBvC+972Pj3zkI0yZMoXXvva1rF69mtNOO40zzjhjyHorV65k/vz5AHtk5oqxiFVD68u/li1bxpw5c+odzkb1vYRh2rFHDLvOmqXXAwxZZ83S69lvFF7CIEnSeDRUDjbcJaj7FDPejgd+SOV18W8ezSAlSZK0vnvvvZfp06dz/fXX84//+I/893//N1de6cpGSZKkiWa4L2GYVrz16Xjg0uIV7BufOqf1RMTnqLzNq9olmfnVesQjSZLGt56eHtasWcP111/Pm970JqZNm7ZuPzgNnznYpunt7F43q2249wND1unt7IZZu252bJIkTTTDHYD7PLAC+A2VPS92BybiHnB1lZnvqncM1b581EkTYgmEJEmN6g1veAOHH344e++9Ny95yUt46KGHJuoecHU13nKwiaC5uXnEdbqnbgVAqVQa/KZZu25S25IkTXTDfQnDvwH/VlX0QES8sjYhSZIkCeCUU07hlFNOWXc+e/ZsLr/88jpGpEaxcOHCeocgSdKkMqwBuIjYEfgw8Iqi6MfAvwCP1SguSZKkhvf4449z6aWX8vOf/xyAl770pbzrXe9i++23r3NkkiRJGonhvoThK8DjwOuLn78C7pkhSZJUQx/84AfZbrvtuOSSS7jkkkuYPn065557br3DkiRJ0ggNdw+45sw8ser8IxFxey0CkiRJUsWDDz7IZz/72XXnZ555Jscdd1wdI9JksnjxYsrl8oDXursrL1QYaD+35uZml6hKkjRCwx2AeyoiXp6ZNwNExMuAp2oXlsbCW6+5mmlDbZI7yn5w0uvHrC9JkiaDrbfeml/84hccdNBBAPzyl79k6623rnNUmizK5TLL776bppkzNrjW29kJwMqeNf3Ku8YkNkmSJpvhDsAtBL5e7AUH0A2cWpuQJEmSBHD++edzzjnnsHr1agB22GEHLrzwwjpHpcmkaeYMph1z9Abla676AcAG1/rKJUnSyAx3AO53wCeAZuA5VF6+cDywvEZxSZIkNbzm5mbOOOMMHnzwQR5//HG23357rr/+evbee+96hyZJkqQRGO4A3JXAX4BfAQ/VLhxJkiT1WbhwITvssAP77LMP2267bb3D0TjV0dEBwIIFC+ocyaaZ6PFLkjQcwx2Am5OZR9U0EkmSJK3n4Ycf5stf/nK9w9A4197eDkzcAazB4u/s7OTjH/845513HjNmbLhP3WA2tZ4kSbU0ZZj3/SQi9qtpJJIkSVrPgQceSGbWOwypLpYsWcKdd97JkiVLxqSeJEm1NOQMuIi4A+gt7ntLRNwHPA00Ab2ZuX/tQxw9EbE6M6cP895jgX0yc4OdjkfSzjD7+iRwDPAMUAbekpl/iYgFwIXAlsW192fmDUWdFwNfA7YBfgi8OzN7I2IG8P+AucAK4PWZ2T1asUqSpNo75phjAHj22Wf5/ve/z5w5c9hyyy3XXb/qqqvqFdomMQfTSHV2dnLdddfR29tLe3s7ra2tw5rNtqn1JEmqtY0tQd3wlUgNIjOXAkvHqLsO4NzM7ImIi4BzgXOAVcAxmfnHiJgHtAOzizqLgbcDt1JJ/o4CrgE+ACzLzAsj4gPF+Tlj9DkkSdIouOyyy+odQt2Yg41cd3c3XV1dtLW1jaheuVymd+pwd6Sp6H3yKcrl8oj72lgc/QfJlixZwtq1awFYu3YtS5Ys4ayzztpoW5taT5KkWhvyv7iZ+cBYBVIvETELuAzYrSh6T2beEhGnAQdl5pkRsQfwH1S+r2v71X8/8HpgK+CKzPxwRMylkojdDPwDlRdXHJeZTw0UQ2ZeV3V6K3BSUf7rqvK7gK0jYitgBrBDZv60iOFyKm+lvQY4DjisqPN14Ec4ACdJ0oQye/bsjd80wZmDaSg33HADPT09APT09LBs2bJhDaRtaj1JkmptZI+8JqdLgE9n5s0RsRuVJ5wvGOCexZl5eUS8q68wIo4E9gJeSmVZ7tKIeAXwYFF+cma+LSK+DZwIfHMY8ZxOZflCfycCv87MpyNiNrCy6tpK/vZUdpfM/BNAZv4pInYeRp+SJEljzRxslJRKJUqlEosWLRpRvba2Nu549JER1WnadhuaZ+084r42Fkd/hx9+ONdeey09PT1MnTqV+fPnD6utTa0nSVKtDfclDJPZEcClEXE7leUOO0TE9v3ueRnwreL4G1XlRxY/vwZ+BexNJekDuD8zby+Of0llP5AhRcR5QA+wpF/5vsBFwP8uipoGqN67sfYlSZLGEXMwDaq1tZUpUyr/qzJlyhRaW1trWk+SpFpzBlxlEPKQ/ksTIqL/fQMlV03ABZn5+X5151J5WUWfZ6ls1DuoiDiVyp578zOzt6p8DnAFcEpmlovilcCcqupzgD8Wxw9HxHOLJ6/PBUb2WFOSJGlsmINpUDNnzuTII4/k6quvpqWlZdgvUtjUepIk1Zoz4OA64My+k4g4YIB7bgHeWBxXP0ZrB06PiOlF3dmbstwgIo6iskfIsZn5ZFX5c4CrqWwOfEtfebG84fGIODgimoBTgCuLy0uBU4vjU6vKJUmSxhNzMA2ptbWVefPmjXgW26bWkySplhptBty2EVG9b8engLOBz0XEcirfx43AO/rVezfwHxHxbuB7fYWZeV1EvAD4afG0djXwJipPW0fiUiobCHcU7dyame+gkpTuCfxzRPxzce+RmfkIsBD4GpWnutcUPwAXAt+OiLdS2QfldSOMRZIkabSZg9VQS0vLWHVVE4PFP3PmTC6++OIRt7ep9SRJqqWm3l63rWg0xfKM++e2/X9MK5XGrN8fnPT6MetLkqRGtnLlyr7N5/fIzBV1Dkf8Lf9atmwZc+bM2djtY6LvJQzTjjl6g2trrvoBwAbX1lz1A/Yb5ZcwSJI0WQyVg7kEVZIkSZIkSaqhRluCWlcR8Tkqb/OqdklmfrUe8UiSJDUCc7DB9XZ2rZvttn55J8AG13o7u2DWiLfbkySp4TkAN4Yy8131jqHal1/1mnGzBEKSJKlWxlsONl40NzcPeq176jQASv23K5m185D1JEnSwByAkyRJkhrQwoUL6x2CJEkNwz3gJEmSJEmSpBpyBpwkSZIkqaEtXryYcrk8Jn11d3cDAyzxrpHm5mZnvErjgANwDext1/yILUsz6h3GuHLlSa+udwiSJEmSxli5XGb53XczZeasmve1tvNRAB7qWTtmfUmqPwfgJEmSJEkNb8rMWWx1zOtq3s/TV30HYEz7klR/7gEnSZIkSZIk1ZADcJIkSZIkSVINOQAnSZIkSaqLjo4OOjo66h2GJhh/bzQRuQecJEmSJKku2tvbAViwYEGdI9FE4u+NJqKGmgEXEb0R8Y2q86kR8WhE/KA4PzYiPlAcnx8RbTWKY25E3DnR2pYkSdoU5mCSJKnRNdoMuCeAeRGxTWY+BSwAHuq7mJlLgaX1Ck6SJGmSMgeTJEkNrdEG4ACuAV4DfBc4GfgWcChARJwGHJSZZ/bdHBE7A9dk5osj4oXA7cDumflgRJSB/YD5wIeALYFOoDUzH46I84HVmbmoaOtO4OjqYCLiecD3gLcDvwIuBA4DtgI+l5mfL+57P/D6ovyKzPzwEJ9xakR8HTgQ+D1wSmY+OfKvSpIkadSYg0naQHd3N11dXbS11WTi67CVy2V6p25Z1xhqoffJJymX/1L373e0lctlZsyYUe8wpBFpqCWohf8E3hgRWwP7A7cNdXNmPgJsHRE7UEkSfwEcGhG7A48USdXNwMGZeWDR/v8ZTiAREVQSv7dk5s+BtwKPZeZLgJcAb4uIPSLiSGAv4KXAAcCLI+IVQzUNfCEz9wf+CrxzOPFIkiTVkDmYJElqWA03Ay4zl0fEXCpPXn84zGo/AV4GvAL4OHAU0ATcVFyfA/y/iHgulSew9w+jzVnAlcCJmXlXUXYksH9EnFSc70gl6Tuy+Pl1UT69KL9xkLb/kJm3FMffBM4GFg0jJkmSpJowB5M0kFKpRKlUYtGi+v5VaWtr485HO+saQy00bbstzbP+vu7f72ibbDP61BgabgCusJRKMnQYMHMY999E5cnr7lQStnOAXuAHxfXPAp/KzKURcRhwflHew/qzDLeuOn4M+AOVpLIv+WsCzsrM9urOI6IFuKBvKcQw9G7kXJIkqR7MwSRJUkNqxCWoAF8B/iUz7xjm/TcCbwLuycy1QBfwaqDvCeeO/G0j4VOr6q0AXgQQES8C9qi69gxwPHBKRPxTUdYOLIyIaUWd50fEdkX56RExvSifXeyLMpjdIuKQ4vhkKsszJEmS6s0cTJIkNaSGnAGXmSuBS0Zw/4rKViHrlhvcDMzJzO7i/HzgOxHxEHArf0vyvkclubsd+DmVzXir230iIo4GOiLiCeBLwFzgVxHRBDwKHJ+Z10XEC4CfFnGsppKMPjJIyL8DTo2IzwP3AIuH+1klSZJqxRxMUn8tLS31DkETkL83moiaenudGd9oiv1X7n9e2z+zZck3x1S78qRX1zsESZI228qVK5k/fz7AHpm5os7hiL/lX8uWLWPOnDn1DkdSP317wG11zOtq3tfTV30HYMz6mjdr5qTbA04ar4bKwRp1CaokSZIkSZI0JhpyCepkEBEzgWUDXJqfmZPv9T2SJEnjgDmYJEnaFA7ATVBFgndAveOQJElqJOZg0uS1tvPRdctDa90PMHZ9zRrOS6cl1ZoDcA3si686zD1IJEmSJDW85ubmMeure2plJ6hSqVT7zmbNHNPPJmlwDsBJkiRJkhrawoUL6x2CpEnOlzBIkiRJkiRJNeQMOEmSJEmSNGEsXryYcrlck7a7u7uBsVki3Nzc7OzLBuIAXAP739f+gi1LK2rS9hUnvrwm7UqSJEmSGlu5XGb53cmUmbuOettrOx8B4KGe2i4YXNv555q2r/HHAThJkiRJkjShTJm5K9sce+qot/vU0q8D1KTtgfpR43APOEmSJEmSxGA22wAAIABJREFUJKmGHICTJEmSJEljrqOjg46OjnqHMe75PU0OLkGVJEmSJEljrr29HYAFCxbUOZLxze9pcnAG3CaIiN6I+EbV+dSIeDQifjDK/Rw22m1KkiRNVOZgkiRponIAbtM8AcyLiG2K8wXAQ3WMR5IkqRGYg0mSpAnJJaib7hrgNcB3gZOBbwGHAkTE+cDqzFxUnN8JHF1V72bgH6gkjMdl5lMRsSdwGTALeBZ4XXH/9Ij4LjAP+CXwpszsjYgXA58CpgOrgNMy808R8SPgNuCVwHOAt2bmTbX6EiRJksaYOZgkSZpwHIDbdP8J/N9iecL+wFcokr+N2As4OTPfFhHfBk4EvgksAS7MzCsiYmsqsxP/HjgQ2Bf4I3AL8LKIuA34LJXE8dGIeAPwr8DpRR9TM/OlEfFq4MPAEaPzkSVJkurOHEySJonu7m66urpoa2sbUb1yuUzv1K1rFNXY6H1yNeXyqmF99nK5zIwZM8YgKtWSA3CbKDOXR8RcKk9efziCqvdn5u3F8S+BuRGxPTA7M68o2v4fgIgA+FlmrizObwfmAn+h8jS2o7hnC+BPVX18v7r9EX40SZKkccscTJIkTUQOwG2epcAi4DBgZlV5D+vvr1c9NP901fGzwDZA0xB99L9/anH/XZl5yEbq9N0vSZI0mZiDSdIkUCqVKJVKLFq0aET12trauPPRx2oU1dho2nY6zbNmD+uzj3SGoMYnX8Kweb4C/Etm3tGvfAXwIoCIeBGwx1CNZOZfgZURcXxRZ6uI2HaoKsCsiDikuH9aROy7aR9BkiRpwjEHkyRJE4oDcJshM1dm5iUDXPoeMKNYrrAQ+P0wmnszcHZELAd+Auw6RL/PACcBF0XEb4DbqWwoLEmSNOmZg0mSpInGqfGbIDOnD1D2I+BHxfFTwJGDVJ9XVWdR1fE9wOH97r2vr83injOrjm8HXjFAHIdVHa/C/UckSdIkYQ4mSZNLS0tLvUOYEPyeJgcH4CRJkiRJ0phbsGBBvUOYEPyeJgeXoEqSJEmSJEk15Aw4SZIkSZI0oazt/DNPLf16TdoFatL2Bv3M2rGmfWh8cQCugX3+qIOYM2dOvcOQJEmSJGnYmpuba9Z299S1AJRKNR4cm7VjTT+Hxh8H4CRJkiRJ0oSxcOHCeocgjZh7wEmSJEmSJEk15ACcJEmSJEmSVEMuQZUkSZI0LIsXL6ZcLg94rbu7G4BSqTTg9ebmZpeNSZIalgNwDezMa+9ly9Jfa97Pt0/cp+Z9SJIkqfbK5TLL7/49W8ycvcG1ZzsfBeCPPVsOcO2hmscmSdJ45gCcJEmSpGHbYuZstj/2rA3KH1/6WYAhr0mS1KjcA06SJEmSJEmqIQfgJEmSJK2no6ODjo6Oho9BkqTR4hJUSZIkSetpb28HYMGCBQ0dgyRJo8UZcDUQEb0RcXHVeVtEnD9Gfc+NiDvHoi9JkqTxxBxMkiSNVw7A1cbTwGsjYqex7DQithjL/iRJksYZczBJkjQuuQS1NnqALwDvBc6rvhARs4DLgN2Kovdk5i0RcQdwKPAYsAp4b2ZeHhHfAL4O3At8A9iuqHdmZv4kIg4DPgz8CTgAeHVVX88Dvge8PTN/XosPKkmSNI6Yg42S7u5uurq6aGtrW6+8XC6zdup2g9Qa3NonH6dc/vMG7Q2lXC4zY8aMEfclSdJ45Ay42vkc0BoRO/YrvwT4dGa+BDgR+FJRfgvwMmBf4D4qiSDAwcCtwCPAgsx8EfAG4N+q2nwpcF5m7tNXEBFBJfF7y0RN/CRJkjaBOZgkSRp3nAFXI5n514i4HDgbeKrq0hHAPpXcDIAdImJ74CbgFcADwGLg7RExG+jKzNVFEnlpRBwAPAs8v6rNn2Xm/VXns4ArgRMz864afDxJkqRxyRxsdJRKJUqlEosWLVqvvK2tjbsefWLE7U3ZdnuaZ+26QXtDGclsOUmSxjsH4GrrM8CvgK9WlU0BDsnM6oSQiLgReBeVZRHnAScAJ1FJCqGylOJh4IVFG/9TVb1/FvQY8AcqT3MndPInSZK0CczBJEnSuOIS1BrKzC7g28Bbq4qvA87sOymeppKZfwB2AvbKzPuAm4E2/pb87Qj8KTPXAm8Ghtrs9xngeOCUiPin0fk0kiRJE4M5mCRJGm8cgKu9i6kkdX3OBg6KiOUR8VvgHVXXbgN+XxzfBMymkgQC/DtwakTcSmXpw5Bz/zPzCeBo4L0RcdxmfwpJkqSJxRxMkiSNGy5BrYHMnF51/DCwbdX5Kiob+A5U781Vxz+haoA0M+8B9q+6/dyi/EfAj6ruWwHMK47/ArxkMz6KJEnShGEONnpaWlrqHcK4iEGSpNHiAJwkSZKk9SxYsKDeIYyLGCRJGi0uQZUkSZIkSZJqyBlwkiRJkobt2c6HeHzpZwcsBwa/Nuv5NY9NkqTxygG4BnbpUXsyZ86ceochSZKkCaK5uXnQa91TZwFQKm234cVZzx+yriRJk50DcJIkSZKGZeHChfUOQZKkCck94CRJkiRJkqQacgBOkiRJkiRJqiGXoEqSJEkadxYvXky5XK53GBNKd3c3AKVSaUz6a25udlmyJA2TA3AN7DPtjzB9xviaBHn+CX9X7xAkSZI0DpTLZe66+x6222m3eocyYTyxahUAjz277Rj09WDN+5CkycQBOEmSJEnj0nY77ca8486tdxgTxp1XXgAwJt9ZX1+SpOEZX9OfJEmSJEmSpEnGAThJkiRJAHR0dNDR0VHvMKRxz78rkkbKJaiSJEmSAGhvbwdgwYIFdY5EGt/8uyJppJwBJ0mSJEmSJNXQuJgBFxG9wKcy833FeRswPTPP38T2dgZuAw7JzD8XZf8OPJiZFw6zjSOAMzPz+E2JYRjtfw44Btg9M3uLsjOAC4CVwHSgDJyfmbcW178JfDcz/6s4nwqsysznRMSewB3A3cDWwF+BSzPzG7WIX5IkTXzmYOZgkiRpbIyLATjgaeC1EXFBZq7a3MYy85GIuAhYBLwpIl4EvBx48XDqF0lVzUTEFsCxwB+BlwE3V11ekpnvKe47ArgyIg7NzN8Po+nMzAOLunsCV0QEJoCSJGkQ5mB/Yw4GdHd309XVRVtbW71DoVwu8+y07esdhgbxzJOPUS4/Pi5+V+qhXC4zY8aMeochaQIZLwNwPcAXgPcC51VfiIhZwGXAbkXRezLzloi4AzgUeAxYBbw3My+PiG8AXy/aOzUiXgn8K5UnqWsiYpuivRcBa4r2biyefB5B5annVsBFVTH8L2Ax8NrMXNE/+Ig4GPg0laeeTwKnZeY9Q3zeI4BfA1cCJ7N+8rdOZl4fEV8G3ga8f4j2Bqp7b0S8r/jsEy75kyRJY8IcbADmYJIkabSNlwE4gM8ByyPiE/3KLwE+nZk3R8RuQDvwAuAWKk8uHwDuo5IIXg4cDCzMzLURsRC4AViamTcW7Z0NPJOZ+0XEvsAPI2Kv4tohwAGZ2V08+SQiDqWS2B2bmSsHif13wMsz89mIOAr4GPCGIT7rycC3gGuAj0TEuzOzZ5B7fwWcWnX+6Yg4f4i2+9fde5j3SpKkxmQONrCGzMFKpRKlUolFixbVOxTa2tpYserpeoehQWy57Y7M3WnncfG7Ug+NOvNP0qYbNy9hyMy/Uknezu536Qjg0oi4HVgK7BAR2wM3Aa8ofhYD+0XEbKArM1cXbd4O3An8e1V7L6d4GpmZd1FZgrBnce26zOyuundeUffoIRI/gOcA34+IO6ksudh3sBsjYivgSCoJ6V+oJGjzh2i7qd/5ezPzgMw8ADhoiHoD1ZUkSVqPOdigzMEkSdKoGTcDcIXPAG8Ftqsqm0JlI98Dip/Zmfk4cCOVJ66HAj8CHgVOopIUVltb/PQZKiF6ot/5H4FngAM2Eve/Au2ZOQ84nsoyiMG8BtgRuCsiVlB54nvyEPcfSOXp7qbYnLqSJKlxmINtyBxMkiSNmnE1AJeZXcC3qSSAfa4Dzuw7iYgDinv/AOwE7JWZ91HZw6ONDZO//m4EWou2XgA8F7h3kHu7gKOBTxTLIAazI/BQcXzaRvo/mcr+JHMzcy7wPOBVEbFBwljsnXI68OWNtLmBiHge8EngsyOtK0mSGos52PrMwSRJ0mgbVwNwhYupJHV9zgYOiojlEfFb4B1V124D+t5MdRMwm0E2063yWWCbYgPhJcApmfnMYDdn5p+ovC3r8xEx2HKDi4BPRsQtQ3UcEdOpLHW4pqr9x4vP8ZqiqDUibo+I3wPnAMdnZm7kM1V1Eb+OiLuB/wQunohv35IkSXVhDmYORktLCy0tLfUOQxr3/LsiaaSaent76x2DxlhEzAXuP+b/fIvpM3atdzjrOf+Ev6t3CJIkTXgrV65k/vz5AHsM9PZQjb2+/GvZsmXMmTOn3uFMCH0vYZh33Ln1DmXCuPPKCwDG5Du788oLmLvTVg37EgZJGshQOdh4nAEnSZIkSZIkTRpT6x3ARBIRZ1C1F0rhxszs/9YwIuLVwMf7Fd+bmSfVKj5JkqTJyByscT2x6sF1s7q0cU+sehBgTL6zJ1Y9CDvtVfN+JGmycABuBDLzS8CXhnnvD4Ef1jaizfOelp2ZM8cln5IkaXybbDmYhqe5ubneIUw43VtUtnEslbaqfWc77eWfkSSNgANwkiRJksadhQsX1jsESZJGjXvASZIkSZIkSTXkDDhJkiRJo2bx4sWUy+XNbqe7uxuAUqm02W3VS3NzszP5JEmAA3AN7cpru3jOWOwP0cBaT5xV7xAkSZLGVLlc5nd338uMmbtvVjtdnZ0APNMzfTTCGnNdnQ/UOwRJ0jjiAJwkSZKkUTVj5u4cdeyHNquNa5d+DGCz26mXvvglSQL3gJMkSZIkSZJqygE4SZIkSZIkqYYcgJMkSZIaREdHBx0dHfUOQxpV/l5LmggcgJMkSZIaRHt7O+3t7fUOQxpVw/297uzs5H3vex9dXV0DXjv77LM5++yzB7wuSZvLAbh+IqI3Ii6uOm+LiPM3o72dI+L+iNi1quzfI+IDI2jjiIj4r02NQZIkabwzB5NUa0uWLOHOO+9kyZIlA167++67ufvuuwe8LkmbywG4DT0NvDYidhqNxjLzEeAiYBFARLwIeDlw8VD1+kSEb6qVJEmNwBxMUs10dnZy3XXX0dvbS3t7+3qz3Do7O9ebQdf/uiSNBhOLDfUAXwDeC5xXfSEiZgGXAbsVRe/JzFsi4g7gUOAxYBXw3sy8PCK+AXy9aO/UiHgl8K/AmZm5JiK2Kdp7EbCmaO/GiDgDOAKYDmxFJXnsi+F/AYuB12bmiv7BR8TBwKeBrYEngdMy857N/1okSZJqyhxsDHR3d9PV1UVbW1vN+iiXy2wxdYeatT9RPPXkXyiX/1rT71oV5XKZGTNmDHnPkiVLWLt2LQBr165lyZIlnHXWWeuu9fT0rLt3zZo1612XpNHgDLiBfQ5ojYgd+5VfAnw6M18CnAh8qSi/BXgZsC9wH5VEEOBg4NbMXAssBL4H/D4zbyyunw08k5n7AW8GvhERWxbXDgHenJkL+jqPiEOL2I4dKPEr/A54eWYeCHwU+NhIP7wkSVKdmINJqokbbrhh3SBbT08Py5YtW+9ab2/vuvPe3t71rkvSaHAG3AAy868RcTmV5OypqktHAPtERN/5DhGxPXAT8ArgASpPRt8eEbOBrsxcXbR5e0TcCfx7VXsvBz5ZXL8rIv4I7Flcuy4zu6vunVfUXZCZfx4i/OcAl0dE80g/tyRJUj2Zg9VeqVSiVCqxaNGimvXR1tbGw4+uqVn7E8U22z6HXWbNqul3rYrhzDI8/PDDufbaa+np6WHq1KnMnz9/vWtXX331ukG4pqam9a5L0mhwBtzgPgO8FdiuqmwKcEhmHlD8zM7Mx4EbqTxxPRT4EfAocBKVpLDa2uKnT9MQ/T/R7/yPwDPAARuJ+1+B9sycBxxPZRmEJEnSRGEOJmnUtba2MmVK5X9/p0yZQmtr63rXpk7929yUadOmrXddkkaDA3CDyMwu4NtUEsA+1wFn9p1ExAHFvX8AdgL2ysz7gJuBNjZM/vq7EWgt2noB8Fzg3kHu7QKOBj5RLIMYzI7AQ8XxaRvpX5IkaVwxB5NUCzNnzuTII4+kqamJlpaW9faMmzlzJi0tLevO+1+XpNHgANzQLqaS1PU5GzgoIpZHxG+Bd1Rduw34fXF8EzCbShI4lM8C2xQbCC8BTsnMZwa7OTP/BBwLfD4iDhrktouAT0bELRvpW5IkabwyB6uRlpaW9QYapMlguL/Xra2tzJs3b8DZba2trey9997svffezn6TVBNN1ZtNqjFExFzg/rPf/x2eU3puvcOZ1FpPnFXvECRJDWjlypV9+xftMcRLAzSG+vKvZcuWMWfOnHqHU1N9e8AddeyHNquda5dW3mOxue3Uy7VLP8Yus6a5B5wkNZChcjBnwEmSJEmSJEk15FtQJ6iIOIOqvVAKN2bm2fWIR5IkqRGYg0mSpE3hANwElZlfAr5U7zgkSZIaiTnY8HR1PrBuCenmtAFsdjv10tX5ALvM2rPeYUiSxgkH4BrYcUfNYM4c9yiTJEnS6Glubh6VdracOhOAUmnaqLQ31naZteeofReSpInPAThJkiRJo2bhwoX1DkGSpHHHlzBIkiRJkiRJNeQMOEmSJEnrWbx4MeVyecz77e7uBqBUKo1531pfc3OzsxklaRQ5ANfAblnaxcznbDWmfR7xT+45J0mSNN6Vy2Xyd/ey84zdx7TfR7o6KwfPTB/TfrW+R7oeqHcIkjTpOAAnSZIkaQM7z9idk1/1oTHt81vXVN54Otb9an19fw6SpNHjHnCSJEmSJElSDTkAJ0mSJE1gHR0ddHR01DsMacz4Oy9pInIJqiRJkjSBtbe3A7BgwYI6RyKNDX/nJU1EzoCTJEmSJEmSasgBuH4iojciLq46b4uI8zejvZ0j4v6I2LWq7N8j4gMjaOOIiPivYd77sYh4z6bE+v+3d+dhclVl4se/IZBAgIROCIsEiHTiq4AaBBUGdNhiBkcFFEeZuIG45BFx69FR+CnjgIgTUNziKALiMKIiosOMhJ4ogoioKIOAvpgCEoNsoZskbMGQ/v1Rt7HT6e50h6q61V3fz/Pk6brbuW+d3Kp++9xzzpUkSSqLOZgkSRrLbIDb2FrgtRGxYy0Ky8wHgLOBhQAR8SLgEOCcoY7rFREOE5YkSa3AHEySJI1ZJhYbWwd8FfgAcGrfDRExHfgKsEex6v2ZeX1E/A54GbAKWAl8IDMvjohvAt8oyntrRBwGnAmcnJl/iYhtivJeBPylKO/aiDgJOBLYDphINXnsjeGlwCLgtZl59yDvYb+I+AkwAzgrMy94RjUiSZJUf+Zgm6m7u5uuri46OjpqVmalUmHi+Mk1K0+jy6OPP0xXZXVNr6laqlQqTJ06tewwJGlE7AE3sC8B8yNiSr/15wGfzcwXA68Dzi/WXw8cDOwD3Ek1EQQ4EPhFZq4HFgDfA+7IzGuL7acAT2bm84E3A9+MiAnFtoOAN2fm0zOLRsTLitheM0TiB/B84Kgipk9GxM4jefOSJEklMQeTJEljkj3gBpCZqyPiYqrJ2eN9Nh0J7B0RvcuTI2J74Drg5cAyqndG3xkRuwFdmflIUebNEXEr8OU+5R0C/Fux/baI+DMwq9h2dWZ299l33+LYuZl53ybewhWZ+QTwRERcC7wYuHL4NSBJktR45mCbp62tjba2NhYuXFizMjs6Oui+/y81K0+jy7bb7MCMnafX9JqqpWbtmSdJQ7EBbnCfA34DXNhn3RbAQZnZNyGkSLDeQ3VYxKnAscBxVJPCvtYX/3qNG+L8j/Zb/jOwLTAHuGoTsfdsYlmSJKlZmYNJkqQxxyGog8jMLuA7wNv7rL4aOLl3ISLmFPv+CdgRmJ2ZdwI/AzrYOPnr71pgflHW84BdgaWD7NsFvAr4TDEMYijHRMTEYhLjlwG/3sT+kiRJTcEcTJIkjUU2wA3tHKpJXa9TgAMi4paIuB14d59tNwJ3FK+vA3ajmgQO5QvANsUEwpcAb8nMJwfbOTPvBV4D/HtEHDBEub8CfgTcAHwiM+/fRBySJEnNxBxsBObNm8e8efMacSqpKXjNSxqNxvX02DO+1UTETOCuT77vu0zbYdeGnvvIf5ze0PNJktSKVqxYwRFHHAHw7E08NEAN0pt/LVmyhBkzZpQdzib1zgF3/FGnNfS83/rRGQANP6829K0fnUHbzls17RxwktSshsrB7AEnSZIkSZIk1ZEPYRilIuIk+syFUrg2M08pIx5JkqRW0Eo52ANdy57ukdbIcwINP6829EDXMtp2nrXpHSVJw2YD3CiVmecD5z+TMg5+zVRmzHBIqCRJ0nDVIgcbDdrb28s58YRpALS1bVXO+QVA286zyrsGJGmMsgFOkiRJ0gYWLFhQdgiSJI0pzgEnSZIkSZIk1ZENcJIkSZIkSVIdOQRVkiRJ0gYWLVpEpVLZaH13dzc77bQTZ511VglRSZI0etkA18Ju+/ZDPDhlQtlhqIXsd9JOZYcgSZKGoVKpcMfvl/KsHfbYYP2KlffQ1dVVUlSSJI1eNsBJkiRJ2sizdtiDdx122gbrPvH9d5QUjSRJo5tzwEmSJEmSJEl1ZAOcJEmS1II6Ozvp7Oxs2vIkSRpLHIIqSZIktaDFixcDMHfu3KYsT5KkscQGuCYTETOALwF7U+2heCXwT5n55BDH3A0ckJkrGxGjJEnSWGMOJkmS6skhqE0kIsYBlwNXZOZs4DnAdsCZpQYmSZI0hpmDSZKkerMHXHM5HHgiMy8EyMynIuIDwF0RcRewd2aeDBARVwILM/Oa3oMjYiZwZWbuWyx3ANtl5ukNfReSJEmjS0vmYN3d3XR1ddHR0bHRtkqlwqQtJm+0/qn1T/Hk42sHPWbq1Kl1iVWSpNHOHnDNZR/gpr4rMnM1sBwbSyVJkurFHEySJNWVCUVzGQf0jGC9JEmSnrmWzMHa2tpoa2tj4cKFG23r6OjgkXs3nv5u/BbjmTBhq0GPkSRJA7MHXHO5DTig74qImAzsDqxiw/+vrQc4ft0w9pEkSdKGzMEkSVJd2QDXXJYAkyLiLQARMR44B7gIuBOYExFbRMTuwEsGOP5+YKeImBYRE4FXNSZsSZKkUc0cTJIk1ZUNcE0kM3uAY4HXR8QfgTuAJ4CPAdcDdwG/AxYCvxng+L8AnwRuBK4E/tCYyCVJkkYvczBJklRvzgHXZDLzT8CrB9k8f5BjZvZ5/Xng87WPTJIkaexqxRxs3rx5TV2eJEljiQ1wkiRJUguaO3duU5cnSdJY4hBUSZIkSZIkqY7sASdJkiRpI39+eDn//pMzNli3dt1aJk3YpqSIJEkavWyAa2H7vGEaM2bsVHYYkiRJajLt7e0Drp+x9W7stJP5oyRJI2UDnCRJkqQNLFiwoOwQJEkaU5wDTpIkSZIkSaoje8BJkiRJGlUWLVpEpVIZ8XHd3d0AtLW1jei49vZ2ewVKkp4RG+Ba2L3nP0jP9uVdArt/aJfSzi1JkqTRq1Kp8Mfbl7LH9nuM6LiVa1YCMOmRbYd9zPI1y0d0DkmSBmIDnCRJkqRRZ4/t9+DDL/nYiI75zC8/BTCi43qPkSTpmXAOOEmSJEmSJKmObICTJEmSJEmS6sgGOEmSJEkb6OzspLOzs+wwWo71Lkljlw1wkiRJUgt66KGH+NCHPkRXV9dG2xYvXszixYtLiKq1We+SNHb5EIYGiYgZwJeAvak2fF4J/FNmPjnI/hOB/wZ2BM7KzG9vxjlnAldm5r6bG7ckSdJoZg42uEsuuYRbb72VSy65hPe+971lhyNJ0phmD7gGiIhxwOXAFZk5G3gOsB1w5hCH7QdslZlzhpv4RcT4ZxysJEnSGGEONriHHnqIq6++mp6eHhYvXjxgLzhJklQ79oBrjMOBJzLzQoDMfCoiPgDcFRHdwKzMPDEing98C3gl8B/A9Ii4GXgdMBNYSPX/7FfAgsxcGxF3AxcArwC+GBF/LJYfA37WuLcoSZLUdMzBBnHJJZewfv16ANavX79RL7ju7m66urro6OgoK8QhVSoVtu+Z3JBzrVq7ijWV1Q2pi0qlwtSpU+t+HklS49kDrjH2AW7quyIzVwPLqQ6DmBURxwIXAu/KzOXAScB1mTkHuAe4CHhDZj6fagK4oE9xT2TmIZl5aVHGKZl5UJ3fkyRJUrMzBxvEj3/8Y9atWwfAunXrWLJkSckRSZI0ttkDrjHGAT1DrH8bcAvw75l5/QD7BXBXZt5RLH8DeA/wuWL52wARMQXYITN/Wqz/JnBULd6AJEnSKGQONojDDz+cq666inXr1rHllltyxBFHbLC9ra2NtrY2Fi5cWFKEQ+vo6GDtnwacxq/mpkycwk67T29IXTRrj0NJ0jNnD7jGuA04oO+KiJgM7A5UgNnAI8CzBjl+3CbKf7TPfgMlmZIkSa3IHGwQ8+fPZ4stqn8KbLHFFsyfP7/kiCRJGttsgGuMJcCkiHgLPD1R7zlUhzRsBZwHvByYFhHHDXD8H4CZETGrWH4z8NP+O2Xmw8CqiDikWGUmJUmSWpk52CCmTZvGK17xCsaNG8e8efOcd0ySpDqzAa4BMrMHOBZ4fTFB7x3AE8DHgM8CXy6GNrwd+HRE7NTv+CeAE4DvRsTvgPXAVwY53QnAlyLiBuDxerwfSZKk0cAcbGjz589n3333tfebJEkN4BxwDZKZfwJePcCmE/vt03uH9QHgmj7blgD7DVDuzH7LNwEv7LPq9M0MWZIkadQzBxvctGnTOOeccwbcNm/evAZHI7DeJWksswFOkiRJ0gbmzp1bdggtyXqXpLHLIaiSJEmSJElSHdkDTpIkSdKos3zNcj7zy0+N8JhlACM6bvma5cx+eoSyJEmbxwa4FrbrSdOZMWOXssOQJEmSRqS9vX2zjtuxe0cBo+8DAAAYn0lEQVQAJrZNGPYxs5m12eeTJKmXDXCSJEmSRpUFCxaUHYIkSSPiHHCSJEmSJElSHdkDTpIkSVLTWbRoEZVKpeww6O7uBqCtra3kSGqvvb3d3oSS1CA2wLWwB7++lC23X13Xc+zywb3rWr4kSZLGpkqlwtLb72CP7XcrNY6Vax4EYNtHhj9v3GiwfM09ZYcgSS3FBjhJkiRJTWmP7XfjYy99b6kxfOrGLwCUHket9b4vSVJjOAecJEmSJEmSVEc2wEmSJEmSJEl1ZAOcJEmSpGHp7Oyks7Oz7DDU4rwOJY1GzgFXIxHxFPA7qnX6e+CtmfnYEPvfDRyQmSs3Ue6xwOXA8zLzD8W6mcCVmblvRMwBnpWZ/1OL9yFJkjSamIM11uLFiwGYO3duyZGolXkdShqN7AFXO49n5pzM3Bd4Enh3jco9HvgZ8MZBts8BXlmjc0mSJI025mCSJKnp2QOuPq4DXgAQEVcAuwNbA+dl5lf77xwRHwROLBbPz8zPFeu3Aw4GDgN+CJze77gJwCeBbSLiEOAsoBO4ANgLeAx4Z2beUuP3J0mS1IzMwSRJUlOyAa7GImJL4CjgqmLViZnZFRHbAL+KiO9l5kN99t8fOAF4KTAOuDEifpqZvwWOAa7KzDsioisiXpSZv+k9NjOfjIiPUx1GcXJR3heA32bmMRFxOHAx1Tu0kiRJY5Y5WGN0d3fT1dVFR0dH3c9VqVSY3LNt3c/TqlatXcPqyn0N+b+stUqlwtSpU8sOQ5JGxCGotbNNRNwM/BpYDny9WH9KRPwf8Auqd2Fn9zvuEOD7mfloZj5Cda6RlxXbjgcuLV5fWixvyiHANwEy88fAtIiYsnlvSZIkqemZg0mSpKZnD7jaeTwzN7jLGRGHAkcCB2XmYxFxDdVhEH2NG6iwiJgGHA7sGxE9wHigJyI+vIk4BiqvZ9PhS5IkjUrmYA3U1tZGW1sbCxcurPu5Ojo6ePJPj9b9PK1qysTtmb77Lg35v6y10dhrT5LsAVdfU4DuIvF7LnDgAPtcCxwTEZMiYlvgWKrzlxwHXJyZe2bmzMzcHbiL6t3VvtYA2/crbz48nXyuzMzVtXxTkiRJTc4cTJIkNRUb4OrrKmDLiLgF+FeqQyA2UMwnchHwS+BGqhMA/5bqUIfv99v9e8A/9lv3E2DviLg5It5AdZLgA4pzfhp4a83ejSRJ0uhgDiZJkpqKQ1BrJDO3G2DdWqqTAQ+0/8w+r88Fzu23/dABjvl8n8V9i3VdwIv77Xr0MMOWJEka1czBGmvevHllhyB5HUoalWyAkyRJkjQsc+fOLTsEyetQ0qjkEFRJkiRJkiSpjmyAkyRJkiRJkurIIaiSJEmSmtLyNffwqRu/UGoMy9bcA1B6HLW2fM09zOI5ZYchSS3DBrgWNv3ts9hlxoyyw5AkSZI20t7eXnYIAOzYPR2ACW3blhxJbc3iOU1Tx5LUCmyAkyRJktR0FixYUHYIkiTVjHPASZIkSZIkSXVkDzhJkiRJNbVo0SIqlcozKqO7uxuAtra2Abe3t7fbS06SNGrYANfCVl54E1tNXtaw8+38voMbdi5JkiSVp1KpsPT2ZI/Ju2x2GStXPwDAto9uPGhn+er7NrtcSZLKYAOcJEmSpJrbY/IunHrQ2zb7+DNvuAhgwDJ6t0mSNFo4B5wkSZIkSZJURzbASZIkSS2os7OTzs7OssMY1axDSdJwOQRVkiRJakGLFy8GYO7cuSVHMnpZh5Kk4Wr5HnARMSMifhARf4yIOyPiixExseSYHqlROTMj4tZalCVJklRL5mBqVUuXLuWYY47hzjvvLDsUSVIDtXQDXESMAy4HrsjM2cBsYBvgM6UGJkmSNIaZg6mVnX322Tz22GOcddZZZYciSWqgVh+CejjwRGZeCJCZT0XEB4BlEfFH4LmZeTJARFwJLMzMayLiFcC/ABOBCnBCZj4SEfsD5wLbASuBt2XmvRFxDXAjcBiwA/D2zLwuIvYBLgQmUG0MfV1m/rE3uIjYDvgB0AZsBZyWmT+IiJnAj4CfAX8D3AMcnZmPFzFcADxWbJckSWo25mBqSUuXLmXZsmUALFu2jDvvvJO99tqr5KgkSY3Q6g1w+wA39V2Rmasj4m4GqZuI2BE4DTgyMx+NiI8AH4yIs4AvUE3CHoyINwBnAicWh26ZmS+JiFcCnwCOBN4NnJeZl0TEBGB8v9M9ARxbxLQj8IuI+GGxbTZwfGa+IyK+A7wO+A+qyeR7M/OnEfFvm1sxkiRJdWQO1gS6u7vp6uqio6Oj5mVXKhUms3XNy+21au0jrK6srEvsI1GpVJg6deqw9z/77LM3WD7rrLP42te+VuuwJElNqNUb4MYBPYOsH8yBwN7A9REB1TunNwAB7At0FuvHA/f2Oe7y4udNwMzi9Q3AqRExA7i8753XPnF8KiJeDqwHdgN2LrbdlZk39y0zIqYAO2TmT4v13wSOGuK9SJIklcEcTC2pt/fbYMuSpLGr1RvgbqN61/JpETGZaoL1EPCcPpt6b+GNAzoz8/h+xz0fuC0zDxrkXGuLn09R1Htm/mdE3Aj8PbA4Ik7KzB/3OWY+MB3YPzP/UtwV3rpfeb1lbsPgyawkSVIzMQdrAm1tbbS1tbFw4cKal93R0cGTK1bVvNxeUyZux/QZu9Ul9pEYaQ+8Pffcc4NGtz333LPWIUmSmlRLP4QBWAJMioi3AETEeOAc4IvAXcCciNgiInYHXlIc8wvg4IiYVRwzKSKeAyQwPSIOKtZvVcwvMqiI2Au4MzM/D/wQeEG/XaYADxSJ32HAkL+hM/NhYFVEHFKsmr/pKpAkSWo4czC1pI985CMbLH/0ox8tKRJJUqO1dANcZvYAxwLHFRP+PgSsz8wzgeupJoC/AxYCvymOeRB4G/CtiLiFajL43Mx8EjgOODsi/g+4merkvEN5A3BrRNwMPBe4uN/2S4ADIuLXVBO5PwzjbZ0AfCkibgAeH8b+kiRJDWUOplY1a9asp3u97bnnnj6AQZJayLienlHXW75uIuJvgG8Br83Mmza1/2hVPMHrrm+/6TPsOnnHhp135/cd3LBzSZLUylasWMERRxwB8OzMvLvkcDapFXKw3vxryZIlzJgxo+xwAOjs7ARg7ty5NS+7dwjqqQe9bbPLOPOGiwAGLOPMGy5iwowppQ9B3Zw6XLp0KR0dHZx77rk2wEnSGDNUDtbqc8BtIDN/ziaGGEiSJKm2zMHKUY+Gt1azOXU4a9YsrrjiijpEI0lqZi09BFWSJEmSJEmqN3vASZIkSaq55avve3oY6eZYtvo+gAHLWL76PmYxZbPLliSp0WyAa2E7nrA/OzfJHCSSJEkaO9rb259xGTt2rwdgQtvGDW2zmFKTc0iS1Cg2wEmSJEmqqQULFpQdgiRJTcU54CRJkiRJkqQ6sgFOkiRJkiRJqiOHoEqSJEl6xhYtWkSlUqlZed3d3QC0tbUNuL29vd2hrpKkUcMGuBb20MXXMGHy1LLDGLGdTn5l2SFIkiSpn0qlwtLb/8AeU6bXpLyVqx4EYNvH1m+0bXmxTZKk0cIGOEmSJEk1sceU6Zx28OtrUtYZ138XYMDyerdJkjRaOAecJEmSJEmSVEc2wEmSJEnaQGdnJ52dnWWHUXNj9X1JkpqfQ1AlSZIkbWDx4sUAzJ07t+RIamusvi9JUvOzAW4TIuJU4B+Bp4D1wLsy88ZB9v0kcG1m/u8Q5R0KPJmZPy+WTwceycyFzzDOTZ5bkiRptDAHkyRJY4kNcEOIiIOAVwEvysy1EbEjMGGw/TPz48Mo9lDgEeDnNQlyZOeWJElqeuZgkiRprLEBbmi7Aiszcy1AZq4EiIiPA68GtqGaxL0rM3si4iLgysy8LCLuBr5R7LcV8HrgCeDdwFMR8SbgvYOdOCI+CJxYLJ6fmZ+LiJlF+fsW+3QA22Xm6f3O/WngNcA64OrM7KhRfUiSJDWCOVjJuru76erqoqNj+G+hUqkwedyg7aQ1tWrtY6yuPDyi+KAa49SpU+sUlSRJg/MhDEO7Gtg9Iu6IiC9HxN8W67+YmS8ukrBtqN6hHcjKzHwRsAjoyMy7ga8An83MOZl53UAHRcT+wAnAS4EDgXdExH7DCTgipgLHAvtk5guAM4b1TiVJkpqHOZgkSRpT7AE3hMx8pEjEXgYcBnw7Iv4ZWBMRHwYmAVOB24D/GqCIy4ufNwGvHcGpDwG+n5mPAkTE5UUMPxzGsaup3uU9PyL+G7hyBOeVJEkqnTlY+dra2mhra2PhwuFPkdfR0cGT9zxUx6j+asrESUzfbfcRxQeMuMecJEm1Yg+4TcjMpzLzmsz8BHAyMB/4MnBcZj4f+Bqw9SCHry1+PsXIGjvHDbJ+HRv+n2103sxcB7wE+B5wDHDVCM4rSZLUFMzBJEnSWGID3BCianafVXOALF6vjIjtgONGWOwaYPtN7HMtcExETIqIbakOZ7gOuB/YKSKmRcREBhh2UcQ0JTP/B3h/EbMkSdKoYQ4mSZLGGoegDm074AsRsQPVO59LgXcCDwO/A+4GfjXCMv8LuCwijuavEwCfFhHv790hM2cUE/r+slh1fmb+Fp5+1P2NwF3AHwYof3vgBxGxNdW7uB8YYXySJEllMweTJEljyrienp6yY1CDFU/yuus7b/1/7Dp59D0FaqeTX1l2CJIkNbUVK1ZwxBFHADy7eACBStabfy1ZsoQZM2aUHc4mdXZ2AjB37txhH9M7B9xpB7++JjGccf13AQYs74zrv8uE3aaNeA64zXlfkiQN11A5mD3gJEmSJG1grDZQjdX3JUlqfs4BJ0mSJEmSJNWRPeAkSZIk1cTyVQ8+PXT0mVq26kGAActbvupBZu02rSbnkSSpEWyAa2HT3nIoO42COUgkSZLU/Nrb22ta3o6TqoN1JrS1bbRt1m7Tan4+SZLqyQa41jQe4L777is7DkmSVAd9fsePLzMObWDM51+vfvWrG37OFStWNPyckiQNZqgczAa41rQrwPz588uOQ5Ik1deuQKXsIASYf0mS1Eo2ysFsgGtNvwJeBtwLPFVyLJIkqfbGU038flV2IHqa+ZckSWPfoDnYuJ6ensaHI0mSJEmSJLWILcoOQJIkSZIkSRrLHILaYiLi74DzqHaLPD8zP11ySE0pIu4G1lAdIrIuMw+IiKnAt4GZwN3AP2Rmd0khli4iLgBeBTyQmfsW6waso4gYR/W6eyXwGPC2zPxNGXGXaZA6Ox14B/BgsdvHMvN/im0fBd5O9To8JTMXNzzokkXE7sDFwC7AeuCrmXme19rAhqiv0/E6G1REbA1cC0ykmhtdlpmfiIhnA5cCU4HfAG/OzCcjYiLVet4feAh4Q2beXUrw0jCN9Pu0rDgbZaSf+/IibZyIGA/8GrgnM1/V4nVxN/4tAEBE7ACcD+wL9AAnAklr1kVQfd+99gI+TvW7tRXr4wPASVSvi98BJ1AdetmS3xvDYQ+4FlL8Uv0ScBSwN3B8ROxdblRN7bDMnJOZBxTL/wwsyczZwJJiuZVdBPxdv3WD1dFRwOzi3zuBRQ2KsdlcxMZ1BvDZ4lqb06dRZG/gjcA+xTFfLj7DrWYd8KHMfB5wIPCeom681gY2WH2B19lQ1gKHZ+YLgTnA30XEgcDZVOttNtBNtaGS4md3Zs4CPlvsJzW7kX6fjnUj/dy3gvcBv++z3Mp1Af4t0Os84KrMfC7wQqrXSEvWRVbNycw5VG/CPQZ8nxasj4jYDTgFOKDoWDCeak7Z6t8bQ7IBrrW8BFiamXcWrdCXAkeXHNNocjTwjeL1N4BjSoyldJl5LdDVb/VgdXQ0cHFm9mTmL4AdImLXxkTaPAaps8EcDVyamWsz8y5gKdXPcEvJzHt7e7Bl5hqqSd9ueK0NaIj6GozXGVBcL48Ui1sV/3qAw4HLivX9r7Pe6+8y4Iii96XUtDbj+3RM24zP/ZgWETOAv6fa04niO60l62IILfdZiYjJwMuBrwNk5pOZ+TAtWBcDOAKoZOYyWrc+tgS2iYgtgUlUHzLk98YQbIBrLbsBf+qzvIKh/zBrZT3A1RFxU0S8s1i3c2beC9UkFtiptOia12B15LU3tJMj4paIuCAi2op11lk/ETET2A+4Ea+1TepXX+B1NqSIGB8RNwMPAJ1UHxv/cGauK3bpWzdP11uxfRUwrbERS5tvmN+nY94IP/dj3eeAD1MdngzV77RWrQvwb4Fee1GdvuLCiPhtRJwfEdvSmnXR3xuBbxWvW64+MvMeYCGwnGrD2yrgJlr7e2OTbIBrLQPdnfcxuAM7ODNfRHU423si4uVlBzTKee0NbhHQTnX4y73AOcV666yPiNgO+B7w/sxcPcSu1hsD1pfX2SZk5lPFkJIZVHsBPm+A3XrrxnrTqDWC79Mxb4Sf+zErInrnp72pz+pW/57zb4GqLYEXAYsycz/gUVpgeOWmRMQE4DXAd8uOpSzFzdyjgWcDzwK2pfp56a+Vvjc2yQa41rIC2L3P8gzgzyXF0tQy88/Fzweojut/CXB/71C24ucD5UXYtAarI6+9QWTm/cUfAOuBr/HX4X/WWSEitqL6x+IlmXl5sdprbRAD1ZfX2fAVQ2uuoTpH1g7FsArYsG6errdi+xSGP7xcKs0Iv09bxjA/92PZwcBrigcPXEp1CNnnaM26APxboI8VwIrM7O1NfxnVBrlWrIu+jgJ+k5n3F8utWB9HAndl5oOZ+RfgcuBvaOHvjeGwAa61/AqYHRHPLlrt3wj8sOSYmk5EbBsR2/e+Bl4B3Eq1rt5a7PZW4AflRNjUBqujHwJviYhxxQTHq3q7abe6fvOTHUv1WoNqnb0xIiYWTyGbDfyy0fGVrZiD5uvA7zPz3D6bvNYGMFh9eZ0NLSKmF095IyK2oZpU/h74CXBcsVv/66z3+jsO+HFmeodXTW0zvk/HtM343I9ZmfnRzJyRmTOp/n3w48ycTwvWBfi3QF+ZeR/wp+Lpn1Cd9+x2WrAu+jmevw4/hdasj+XAgRExqfj90ntttOT3xnCN6+kxX2wlEfFKqne0xgMXZOaZJYfUdCJiL6p3uqDa7fo/M/PMiJgGfAfYg+oXzuszs2V7PETEt4BDgR2B+4FPAFcwQB0VX8pfpPqUxceAEzLz12XEXaZB6uxQqsMCe6g+tvxdvQ1GEXEq1Ue9r6M6VOhHDQ+6ZBFxCHAd1Ueb985L8zGq8xZ5rfUzRH0dj9fZoCLiBVQnCh5P9ebkdzLzk8Xvg0uBqcBvgTdl5tqI2Br4JtU5tLqAN2bmneVELw3PSL9PSwmygUb6uS8v0saKiEOBjsx8VavWhX8LbCgi5lB9OMcE4E7gBIrPDC1WFwARMYnqPLB7ZeaqYl2rXhv/AryBag75W+AkqnO+tdz3xnDZACdJkiRJkiTVkUNQJUmSJEmSpDqyAU6SJEmSJEmqIxvgJEmSJEmSpDqyAU6SJEmSJEmqIxvgJEmSJEmSpDqyAU6StIGIeH/xiHVJkiRJUg3YACdJ6u/9gA1wkiRJklQj43p6esqOQZI0QhHxFqAD6AFuAU4DLgCmAw8CJ2Tm8oi4CLgyMy8rjnskM7eLiEOB04GVwL7ATcCbgPcCC4EEVmbmYQ18W5IkSaNKRFwB7A5sDZyXmV+NiLcDHwH+DPwRWJuZJ0fEdOArwB7F4e/PzOvLiFtS49kDTpJGmYjYBzgVODwzXwi8D/gicHFmvgC4BPj8MIraj2pvt72BvYCDM/PzVJPFw2x8kyRJ2qQTM3N/4ADglIjYDfh/wIHAXOC5ffY9D/hsZr4YeB1wfqODlVQeG+AkafQ5HLgsM1cCZGYXcBDwn8X2bwKHDKOcX2bmisxcD9wMzKxDrJIkSWPZKRHxf8AvqPaEezPw08zsysy/AN/ts++RwBcj4mbgh8DkiNi+4RFLKsWWZQcgSRqxcVSHng6ld/s6ipstETEOmNBnn7V9Xj+FvxMkSZKGrZjS40jgoMx8LCKuoTqNx/MGOWSLYt/HGxOhpGZiDzhJGn2WAP8QEdMAImIq8HPgjcX2+cDPitd3A/sXr48GthpG+WsA78ZKkiQNbQrQXTS+PZfqsNNJwN9GRFtEbEl1qGmvq4GTexciYk5Do5VUKhvgJGmUyczbgDOBnxZDHs4FTgFOiIhbqA59eF+x+9eoJoG/BF4KPDqMU3wV+FFE/KTmwUuSJI0dVwFbFvnXv1IdhnoP8CngRuB/gduBVcX+pwAHRMQtEXE78O7GhyypLD4FVZIkSZKkGomI7TLzkaIH3PeBCzLz+2XHJalc9oCTJEmSJKl2Ti8etHArcBdwRcnxSGoC9oCTJEmSJEmS6sgecJIkSZIkSVId2QAnSZIkSZIk1ZENcJIkSZIkSVId2QAnSZIkSZIk1ZENcJIkSZIkSVId2QAnSZIkSZIk1dH/B+93PsYHPcCmAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1440x720 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# We can use the packages matplotlib (plt) and seaborn (sns) to plot figures\n",
    "# Figure 1 shows that patients are not equally distributed across hospitals\n",
    "# Figure 2 shows that ages     are not equally distributed across hospitals\n",
    "\n",
    "# We first print some statistics\n",
    "print('There are {} patients'.format(len(df)))\n",
    "print('There are {} hospitals'.format(len(df.hospital.unique())))\n",
    "print('There are {} datasets'.format(len(df.dataset.unique())))\n",
    "print('There are {:.2f}% of women'.format(100*np.mean(df.sex == 'f')))\n",
    "print('Age statistics : ')\n",
    "print(df.age.describe())\n",
    "\n",
    "# We plot the patient count over hospitals\n",
    "plt.figure(figsize=(20,10)) \n",
    "plt.subplot(121) \n",
    "sns.countplot(y='hospital', data=df)\n",
    "plt.title('Figure 1. Patient distribution per hospital')\n",
    "\n",
    "# We plot the age distribution over hospitals\n",
    "plt.subplot(122)\n",
    "sns.boxplot(x='age', y='hospital', data=df)\n",
    "plt.title('Figure 2. Age distribution per hospital');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# In this cell, we plot one slice of the MRI of a subject\n",
    "\n",
    "# We first load the 2 files associated to raw data and segmentation\n",
    "import os\n",
    "from IPython.display import clear_output\n",
    "\n",
    "# df = pd.read_csv('demographics.csv')\n",
    "\n",
    "#Collect all pre-processed images from .anat directories\n",
    "skullstripped_img_list = []\n",
    "ori_img_list = []\n",
    "seg_img_list = []\n",
    "file_id = []\n",
    "sample_size = 0\n",
    "root_path = \"/shared/ixi-dataset/\"\n",
    "for i in range(19):\n",
    "    dir_name = root_path + \"IXI-T1-\" + str(i)\n",
    "    anat_directories = [x for x in os.listdir(dir_name) if x.endswith('.anat')]\n",
    "    for anat_directory in anat_directories:\n",
    "        ori_file_name = dir_name + '/' + anat_directory + '/T1_biascorr.nii.gz'\n",
    "        seg_file_name = dir_name + '/' + anat_directory + \"/T1_fast_seg.nii.gz\"\n",
    "        if anat_directory[0:3] == 'IXI' and os.path.exists(ori_file_name) and os.path.exists(seg_file_name):\n",
    "            ori_img = nib.load(ori_file_name).get_data()\n",
    "            seg_img = nib.load(seg_file_name).get_data()\n",
    "            skullstripped_img = ori_img * (seg_img > 0)\n",
    "            ori_img_list.append(ori_img)\n",
    "            seg_img_list.append(seg_img)\n",
    "            skullstripped_img = (skullstripped_img - np.mean(skullstripped_img)) / np.std(skullstripped_img)\n",
    "            skullstripped_img_list.append(skullstripped_img)\n",
    "            file_id.append(anat_directory[3:6])\n",
    "            sample_size += 1\n",
    "            clear_output(wait=True)\n",
    "            #print(sample_size)\n",
    "\n",
    "print(len(skullstripped_img_list))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Validate the sample\n",
    "dimension = 1\n",
    "s = 90\n",
    "sample = 5\n",
    "\n",
    "plt.figure(figsize=(10, 5))\n",
    "print('Shape of the MRI : {}'.format(ori_img_list[sample].shape))\n",
    "plt.subplot(131)\n",
    "plt.imshow(np.take(ori_img_list[sample], s, dimension))\n",
    "plt.title('Original');\n",
    "\n",
    "plt.subplot(132)\n",
    "plt.imshow(np.take(skullstripped_img_list[sample], s, dimension))\n",
    "plt.title('Skull Stripped');\n",
    "\n",
    "plt.subplot(133)\n",
    "plt.imshow(np.take(seg_img_list[sample], s, dimension), cmap='Paired')\n",
    "plt.title('Segmented');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape of the MRI : torch.Size([536, 3, 224, 224])\n",
      "torch.Size([536, 3, 224, 224])\n"
     ]
    }
   ],
   "source": [
    "#Slice the most central/important parts of the images (so that they have the same shape) \n",
    "#and convert the array into a tensor\n",
    "new_img_list = []\n",
    "\n",
    "for i in range(len(skullstripped_img_list)):    \n",
    "    new_img = skullstripped_img_list[i][16:240,90,:]\n",
    "    #modify the dimensions by cropping and zero-padding so that all the images' sizes become 224*224\n",
    "    d_2 = new_img.shape[1]\n",
    "    if d_2 % 2 == 0:\n",
    "        pad_width = int((224 - d_2) / 2)\n",
    "        new_img = np.pad(np.asarray(new_img),((0,0),(pad_width, pad_width)),\"constant\") \n",
    "    else: \n",
    "        pad_width_0 = int((224 - d_2) / 2)\n",
    "        pad_width_1 = pad_width_0 + 1\n",
    "        new_img = np.pad(np.asarray(new_img),((0,0),(pad_width_0, pad_width_1)),\"constant\") \n",
    "    final_img = np.tile(new_img,(3,1)).reshape((3,224,224))\n",
    "    new_img_list.append(final_img)    \n",
    "new_img_list = torch.from_numpy(np.asarray(new_img_list))\n",
    "print('Shape of the MRI : {}'.format(new_img_list.shape))\n",
    "\n",
    "#Saves the tensor\n",
    "torch.save(new_img_list,\"data.pt\")\n",
    "data = torch.load(\"data.pt\")\n",
    "print(data.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "torch.Size([536, 3, 224, 224])\n"
     ]
    }
   ],
   "source": [
    "data = torch.load(\"../../../data/data.pt\")\n",
    "print(data.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([array([0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         2.53105114e-03, 2.18895901e-02, 3.86578514e-02, 5.35672665e-02,\n",
       "         6.41067102e-02, 7.25896771e-02, 7.69794641e-02, 8.07760367e-02,\n",
       "         8.14483464e-02, 7.93720958e-02, 7.76320842e-02, 7.56941251e-02,\n",
       "         7.33411218e-02, 7.71772023e-02, 7.29851534e-02, 7.58523566e-02,\n",
       "         7.70981279e-02, 8.06178680e-02, 8.58776811e-02, 9.42420179e-02,\n",
       "         1.04860558e-01, 1.14431086e-01, 1.30645617e-01, 1.45594604e-01,\n",
       "         1.60207497e-01, 1.74128267e-01, 1.78458686e-01, 1.86012283e-01,\n",
       "         1.88780618e-01, 1.89314613e-01, 1.86308891e-01, 1.88088534e-01,\n",
       "         1.81147925e-01, 1.83916259e-01, 1.79170543e-01, 1.79487119e-01,\n",
       "         1.80455841e-01, 1.82452997e-01, 1.76995423e-01, 1.78617070e-01,\n",
       "         1.76955684e-01, 1.76817651e-01, 1.79407634e-01, 1.79348702e-01,\n",
       "         1.81622299e-01, 1.79526666e-01, 1.86091177e-01, 1.88128285e-01,\n",
       "         1.90125443e-01, 1.96630609e-01, 1.93803377e-01, 1.93229517e-01,\n",
       "         2.02345674e-01, 2.05291531e-01, 2.03729843e-01, 2.14506119e-01,\n",
       "         2.19924615e-01, 2.22178833e-01, 2.27181121e-01, 2.38610374e-01,\n",
       "         2.40153771e-01, 2.47310844e-01, 1.24179427e-01, 2.47172427e-03,\n",
       "         7.90955195e-05, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 1.97738799e-05, 1.18643279e-04,\n",
       "         3.55926752e-04, 2.05648351e-03, 1.45139020e-02, 1.78162658e-02,\n",
       "         1.71044061e-02, 1.68669733e-02, 1.44349323e-02, 1.54631741e-02,\n",
       "         1.46918654e-02, 1.27343786e-02, 1.27343786e-02, 1.31099687e-02,\n",
       "         1.18445541e-02, 1.17060354e-02, 1.06976690e-02, 1.10733727e-02,\n",
       "         1.07371237e-02, 1.00253571e-02, 9.84739218e-03, 9.49138006e-03,\n",
       "         9.01688923e-03, 9.98580934e-03, 9.19477443e-03, 7.94909971e-03,\n",
       "         8.64118551e-03, 8.34450497e-03, 8.32480343e-03, 7.61287775e-03,\n",
       "         8.54231611e-03, 7.19769228e-03, 7.71174630e-03, 7.94909971e-03]),\n",
       "  array([0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         2.53105114e-03, 2.18895901e-02, 3.86578514e-02, 5.35672665e-02,\n",
       "         6.41067102e-02, 7.25896771e-02, 7.69794641e-02, 8.07760367e-02,\n",
       "         8.14483464e-02, 7.93720958e-02, 7.76320842e-02, 7.56941251e-02,\n",
       "         7.33411218e-02, 7.71772023e-02, 7.29851534e-02, 7.58523566e-02,\n",
       "         7.70981279e-02, 8.06178680e-02, 8.58776811e-02, 9.42420179e-02,\n",
       "         1.04860558e-01, 1.14431086e-01, 1.30645617e-01, 1.45594604e-01,\n",
       "         1.60207497e-01, 1.74128267e-01, 1.78458686e-01, 1.86012283e-01,\n",
       "         1.88780618e-01, 1.89314613e-01, 1.86308891e-01, 1.88088534e-01,\n",
       "         1.81147925e-01, 1.83916259e-01, 1.79170543e-01, 1.79487119e-01,\n",
       "         1.80455841e-01, 1.82452997e-01, 1.76995423e-01, 1.78617070e-01,\n",
       "         1.76955684e-01, 1.76817651e-01, 1.79407634e-01, 1.79348702e-01,\n",
       "         1.81622299e-01, 1.79526666e-01, 1.86091177e-01, 1.88128285e-01,\n",
       "         1.90125443e-01, 1.96630609e-01, 1.93803377e-01, 1.93229517e-01,\n",
       "         2.02345674e-01, 2.05291531e-01, 2.03729843e-01, 2.14506119e-01,\n",
       "         2.19924615e-01, 2.22178833e-01, 2.27181121e-01, 2.38610374e-01,\n",
       "         2.40213093e-01, 2.52629994e-01, 2.59392632e-01, 2.74539358e-01,\n",
       "         2.85495278e-01, 3.04793263e-01, 3.26702627e-01, 3.55672778e-01,\n",
       "         3.79162503e-01, 4.20212900e-01, 4.57526049e-01, 4.95513656e-01,\n",
       "         5.40753949e-01, 5.83445571e-01, 6.19594752e-01, 6.59198976e-01,\n",
       "         6.80851281e-01, 6.96353935e-01, 7.03277812e-01, 6.97757875e-01,\n",
       "         6.79862591e-01, 6.56235752e-01, 6.28371631e-01, 6.02645925e-01,\n",
       "         5.74942839e-01, 5.52304239e-01, 5.29898136e-01, 5.08067867e-01,\n",
       "         4.95551056e-01, 4.80267995e-01, 4.68856434e-01, 1.43320460e-01,\n",
       "         5.93216397e-04, 1.77963376e-04, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "         0.00000000e+00, 0.00000000e+00, 1.97738799e-05, 1.18643279e-04,\n",
       "         3.55926752e-04, 2.05648351e-03, 1.45139020e-02, 1.78162658e-02,\n",
       "         1.71044061e-02, 1.68669733e-02, 1.44349323e-02, 1.54631741e-02,\n",
       "         1.46918654e-02, 1.27343786e-02, 1.27343786e-02, 1.31099687e-02,\n",
       "         1.18445541e-02, 1.17060354e-02, 1.06976690e-02, 1.10733727e-02,\n",
       "         1.07371237e-02, 1.00253571e-02, 9.84739218e-03, 9.49138006e-03,\n",
       "         9.01688923e-03, 9.98580934e-03, 9.19477443e-03, 7.94909971e-03,\n",
       "         8.64118551e-03, 8.34450497e-03, 8.32480343e-03, 7.61287775e-03,\n",
       "         8.54231611e-03, 7.19769228e-03, 7.71174630e-03, 7.94909971e-03]),\n",
       "  array([0.        , 0.        , 0.        , 0.        , 0.        ,\n",
       "         0.        , 0.        , 0.        , 0.        , 0.        ,\n",
       "         0.        , 0.        , 0.        , 0.        , 0.        ,\n",
       "         0.        , 0.        , 0.        , 0.        , 0.        ,\n",
       "         0.        , 0.        , 0.        , 0.        , 0.        ,\n",
       "         0.        , 0.        , 0.        , 0.        , 0.        ,\n",
       "         0.        , 0.        , 0.        , 0.        , 0.        ,\n",
       "         0.        , 0.        , 0.        , 0.        , 0.        ,\n",
       "         0.        , 0.        , 0.        , 0.        , 0.        ,\n",
       "         0.        , 0.        , 0.        , 0.        , 0.        ,\n",
       "         0.        , 0.        , 0.00253105, 0.02188959, 0.03865785,\n",
       "         0.05356727, 0.06410671, 0.07258968, 0.07697946, 0.08077604,\n",
       "         0.08144835, 0.0793721 , 0.07763208, 0.07569413, 0.07334112,\n",
       "         0.0771772 , 0.07298515, 0.07585236, 0.07709813, 0.08061787,\n",
       "         0.08587768, 0.09424202, 0.10486056, 0.11443109, 0.13064562,\n",
       "         0.1455946 , 0.1602075 , 0.17412827, 0.17845869, 0.18601228,\n",
       "         0.18878062, 0.18931461, 0.18630889, 0.18808853, 0.18114792,\n",
       "         0.18391626, 0.17917054, 0.17948712, 0.18045584, 0.182453  ,\n",
       "         0.17699542, 0.17861707, 0.17695568, 0.17681765, 0.17940763,\n",
       "         0.1793487 , 0.1816223 , 0.17952667, 0.18609118, 0.18812829,\n",
       "         0.19012544, 0.19663061, 0.19380338, 0.19322952, 0.20234567,\n",
       "         0.20529153, 0.20372984, 0.21450612, 0.21992462, 0.22217883,\n",
       "         0.22718112, 0.23861037, 0.24021309, 0.25262999, 0.25939263,\n",
       "         0.27453936, 0.28549528, 0.30479326, 0.32670263, 0.35567278,\n",
       "         0.3791625 , 0.4202129 , 0.45752605, 0.49551366, 0.54075395,\n",
       "         0.58344557, 0.61959475, 0.65919898, 0.68085128, 0.69635394,\n",
       "         0.70327781, 0.69775787, 0.67986259, 0.65623575, 0.62837163,\n",
       "         0.60264592, 0.57494284, 0.55230424, 0.52989814, 0.50806787,\n",
       "         0.4955906 , 0.48066347, 0.47389875, 0.464961  , 0.45379077,\n",
       "         0.45408344, 0.45242637, 0.4585365 , 0.47312552, 0.50253338,\n",
       "         0.54562067, 0.60021115, 0.69060276, 0.78874052, 0.89276316,\n",
       "         0.94351068, 0.92940384, 0.84278253, 0.68668753, 0.52817553,\n",
       "         0.37888731, 0.26188527, 0.1752346 , 0.11749639, 0.08162658,\n",
       "         0.05576186, 0.04085284, 0.03132183, 0.02669451, 0.02183036,\n",
       "         0.01971439, 0.01848858, 0.01736147, 0.01692629, 0.01443493,\n",
       "         0.01546317, 0.01469187, 0.01273438, 0.01273438, 0.01310997,\n",
       "         0.01184455, 0.01170604, 0.01069767, 0.01107337, 0.01073712,\n",
       "         0.01002536, 0.00984739, 0.00949138, 0.00901689, 0.00998581,\n",
       "         0.00919477, 0.0079491 , 0.00864119, 0.0083445 , 0.0083248 ,\n",
       "         0.00761288, 0.00854232, 0.00719769, 0.00771175, 0.0079491 ])],\n",
       " array([-2.    , -1.9725, -1.945 , -1.9175, -1.89  , -1.8625, -1.835 ,\n",
       "        -1.8075, -1.78  , -1.7525, -1.725 , -1.6975, -1.67  , -1.6425,\n",
       "        -1.615 , -1.5875, -1.56  , -1.5325, -1.505 , -1.4775, -1.45  ,\n",
       "        -1.4225, -1.395 , -1.3675, -1.34  , -1.3125, -1.285 , -1.2575,\n",
       "        -1.23  , -1.2025, -1.175 , -1.1475, -1.12  , -1.0925, -1.065 ,\n",
       "        -1.0375, -1.01  , -0.9825, -0.955 , -0.9275, -0.9   , -0.8725,\n",
       "        -0.845 , -0.8175, -0.79  , -0.7625, -0.735 , -0.7075, -0.68  ,\n",
       "        -0.6525, -0.625 , -0.5975, -0.57  , -0.5425, -0.515 , -0.4875,\n",
       "        -0.46  , -0.4325, -0.405 , -0.3775, -0.35  , -0.3225, -0.295 ,\n",
       "        -0.2675, -0.24  , -0.2125, -0.185 , -0.1575, -0.13  , -0.1025,\n",
       "        -0.075 , -0.0475, -0.02  ,  0.0075,  0.035 ,  0.0625,  0.09  ,\n",
       "         0.1175,  0.145 ,  0.1725,  0.2   ,  0.2275,  0.255 ,  0.2825,\n",
       "         0.31  ,  0.3375,  0.365 ,  0.3925,  0.42  ,  0.4475,  0.475 ,\n",
       "         0.5025,  0.53  ,  0.5575,  0.585 ,  0.6125,  0.64  ,  0.6675,\n",
       "         0.695 ,  0.7225,  0.75  ,  0.7775,  0.805 ,  0.8325,  0.86  ,\n",
       "         0.8875,  0.915 ,  0.9425,  0.97  ,  0.9975,  1.025 ,  1.0525,\n",
       "         1.08  ,  1.1075,  1.135 ,  1.1625,  1.19  ,  1.2175,  1.245 ,\n",
       "         1.2725,  1.3   ,  1.3275,  1.355 ,  1.3825,  1.41  ,  1.4375,\n",
       "         1.465 ,  1.4925,  1.52  ,  1.5475,  1.575 ,  1.6025,  1.63  ,\n",
       "         1.6575,  1.685 ,  1.7125,  1.74  ,  1.7675,  1.795 ,  1.8225,\n",
       "         1.85  ,  1.8775,  1.905 ,  1.9325,  1.96  ,  1.9875,  2.015 ,\n",
       "         2.0425,  2.07  ,  2.0975,  2.125 ,  2.1525,  2.18  ,  2.2075,\n",
       "         2.235 ,  2.2625,  2.29  ,  2.3175,  2.345 ,  2.3725,  2.4   ,\n",
       "         2.4275,  2.455 ,  2.4825,  2.51  ,  2.5375,  2.565 ,  2.5925,\n",
       "         2.62  ,  2.6475,  2.675 ,  2.7025,  2.73  ,  2.7575,  2.785 ,\n",
       "         2.8125,  2.84  ,  2.8675,  2.895 ,  2.9225,  2.95  ,  2.9775,\n",
       "         3.005 ,  3.0325,  3.06  ,  3.0875,  3.115 ,  3.1425,  3.17  ,\n",
       "         3.1975,  3.225 ,  3.2525,  3.28  ,  3.3075,  3.335 ,  3.3625,\n",
       "         3.39  ,  3.4175,  3.445 ,  3.4725,  3.5   ], dtype=float32),\n",
       " <a list of 3 Lists of Patches objects>)"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD1CAYAAABJE67gAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAADg9JREFUeJzt3XGsnfVdx/F3uWzMbAiZLCnSGklWv70VFomsmMwJBpxlmVSTyVpqwhwu+gcqmZiJGNbUmLAtURutyQwQ2UJWGVPXmSqLUbLNjKVMNytcvqaWCZfSMXFUySJQvP5xTuH0cG/vw+1zznO+97xfSZPzPOfHPd+bls/93u/ze85Zs7CwgCSpnjO6LkCStDIGuCQVZYBLUlEGuCQVZYBLUlFnjuuFIuIs4O3AU8BL43pdSSpsBjgfOJCZzw8/ObYApxfeXxrj60nSavFO4MvDJ8cZ4E8B3HPPPaxdu3aMLytJNR09epQdO3ZAPz+HjTPAXwJYu3Yt69atG+PLSlJ5i46dvYgpSUUZ4JJUlAEuSUUZ4JJUlAEuSUUZ4JJUlAEuSUUZ4JJaN7dxtusSpoIBLklFGeCSVJQBLmkk5jbOOkoZMQNckooywCWpKANckooywCWpKANckooywCWpKANckooa50eqSVrl3Pc9XnbgklSUAS5JRRngklSUAS5JRRngklSUAS5JRRngklSUAS5JRRngklSUAS5JRRngkkbK2+tHxwCXpKIMcEkqygCXpKIMcEkqygCXpKIafaBDRGwBdgMzwB2ZefvQ8z8A3A2c21/zW5m5v+VaJUkDlu3AI2IG2ANcDWwCtkfEpqFlvwPcm5mXANuAP2m7UEnSyZqMUDYDhzLzcGa+AOwFtg6tWQC+t//4HOBIeyVKkhbTZIRyAfDEwPE8cNnQmp3AFyLiV4E3Ale1Up0kaUlNOvA1i5xbGDreDvxZZq4D3g18KiK8QCpJI9QkZOeB9QPH63j1iOQG4F6AzPwK8AbgvDYKlCQtrkmAHwA2RMSFEfF6ehcp9w2teRy4EiAiZukF+LfbLFTSZPM9T8Zv2QDPzOPAjcD9wBy93SYPR8SuiLimv+w3gA9GxDeATwPvz8zhMYskqUWN9oH393TvHzp328DjR4B3tFuaJOlUvNAoSUUZ4JJUlAEuSUUZ4JJUlAEuSUUZ4JJUlAEuSUUZ4JJUlAEuSUUZ4JJUlAEuSUUZ4JJUlAEuSUUZ4JJUlAEuSUUZ4JJUlAEuFXbx3Rd3XYI61OgTeSRNlsHgPvH44PUHuypHHbEDl6Si7MAlnRY/jb47duCSRm5u46xBPwIGuFTMUhcuvaA5fQxwSSrKAJekogxwSSrKAJekotxGKBXR5CKlN/VMFztwSSrKAJekogxwSSrKAJekogxwSSrKAJdWIW+rnw4GuCQVZYBLBdhRazGNbuSJiC3AbmAGuCMzb19kzbXATmAB+EZmXtdinZKkIct24BExA+wBrgY2AdsjYtPQmg3ALcA7MvOHgZtGUKskaUCTEcpm4FBmHs7MF4C9wNahNR8E9mTmdwAy8+l2y5QkDWsyQrkAeGLgeB64bGjNDwFExD/SG7PszMy/baVCSSvi+6Ksfk068DWLnFsYOj4T2ABcAWwH7oiIc0+vNEmTzo9J61aTAJ8H1g8crwOOLLLmc5n5YmY+BiS9QJckjUiTAD8AbIiICyPi9cA2YN/Qmr8CfhIgIs6jN1I53GahkqSTLTsDz8zjEXEjcD+9+fZdmflwROwCHsrMff3n3hURjwAvAb+Zmc+MsnBpGrj/W6fSaB94Zu4H9g+du23g8QLwof4fSdIYeCemtMrZxa9eBrgkFWWAS1JRBrgkFWWAS1JRBrg0BS6++2IvZq5CBrgkFdVoH7gkDfI9UCaDAS5NIMcdasIRiiQVZYBLU6Trzt7RS7sMcEkqygCXpKIMcEkqygCXpKIMcGnKeFfm6mGAS3pN3EkyOQxwSSrKAJcmzLjGG45R6jPAJakoA1ySijLApSnmjpTaDHBJjXaWzG2cdQfKhDHAJXHtLb6zdEUGuCTADrsiA1ySivL3JmlCTPLFRDvzyWSASzqJYV2HIxRJKsoAl6SiDHBJKsoAl6SiDHBJKsoAl6SiDHBpAkzyHnBNrkb7wCNiC7AbmAHuyMzbl1j3XuAzwNsz86HWqpQkvcqyHXhEzAB7gKuBTcD2iNi0yLqzgV8Dvtp2kZKkV2syQtkMHMrMw5n5ArAX2LrIut8FPgb8b4v1SZKW0CTALwCeGDie7597WURcAqzPzL9usTZJ0ik0mYGvWeTcwokHEXEG8AfA+1uqSZLUQJMOfB5YP3C8DjgycHw2cBHwQER8E/gxYF9EXNpSjZKkRTTpwA8AGyLiQuBJYBtw3YknM/MYcN6J44h4ALjZXSiSNFrLduCZeRy4EbgfmAPuzcyHI2JXRFwz6gIlSYtrtA88M/cD+4fO3bbE2itOvyxpOkzjDTwn3m989tG5jiupzzsxJakoA1ySijLAJakoA1ySijLApY5M4wVMtcsAl6SiDHBJKsoAl6SiDHBJKsoAl6SiDHBJKqrRe6FIao/bB9UWO3BJKsoAl6SiDHBJKsoAl8bI+bfaZIBLUlEGuCQVZYBLUlEGuCQV5Y080hh48VKjYAcuSUUZ4NKI2X1rVAxwSSrKAJekoryIKY2Io5NTm9s4y+yjc12XUZoduCQVZYBLI2D3rXEwwKUROPjY412XoClggEtSUV7ElNqy85yuK9CUsQOXpKLswKXTZeetjtiBS1JRduDSSth1awI0CvCI2ALsBmaAOzLz9qHnPwT8EnAc+Dbwgcz8j5ZrlSaD4a0JsewIJSJmgD3A1cAmYHtEbBpa9s/ApZn5NuA+4GNtFypJOlmTGfhm4FBmHs7MF4C9wNbBBZn5D5n53f7hg8C6dsuUOmbXrQnUZIRyAfDEwPE8cNkp1t8A/M3pFCVNjMHgNsQ1YZoE+JpFzi0stjAifgG4FLj8dIqSOnEioHce67YOqaEmAT4PrB84XgccGV4UEVcBtwKXZ+bz7ZQndcBOW0U0CfADwIaIuBB4EtgGXDe4ICIuAT4BbMnMp1uvUqvfcqE5iq7YoFZxy17EzMzjwI3A/cAccG9mPhwRuyLimv6yjwNvAj4TEV+PiH0jq1jTa+c5i4fuUueX+hrSKtFoH3hm7gf2D527beDxVS3XpWmxkuBd6r8xxDVlvJVe3TFIp97cxlnmNs52XUZZBrgkFWWAS1JRvpmVxs/RidQKO3BJKsoAl6SiDHCNl+MTqTUGuCQArr3lTK69xctilRjgklSUAS5JRfn7ksbD2bfUOjtwaUQOPvY4Bx97vOsytIoZ4JJUlAGu0XN8omX4hlYrY4BLUlFexNTo2HlLI2UHLklFGeCSVJQBrtFwfCKNnDNwtcvglsbGDlyv9lo+5V1qiZ+P+drZgatnscA+6ZPgj518fvhY0tgZ4GoWwMNrDG2pc45QJKkoO/BpZhctlWYHLklFGeCSJoo7UZpzhDKNHJ1Iq4IduCQVZYBLmjje1NOMAT5tHJ9Iq4YBLmli2YWfmgEuaaI5Tlmau1CmhaMTadUxwCWVMNiFz2470ju39/uZfXSuq5I61yjAI2ILsBuYAe7IzNuHnj8L+CTwo8AzwPsy85vtlipJr4T3y4+X+u1ysXfMHDy3Ciwb4BExA+wBfgqYBw5ExL7MfGRg2Q3AdzLzrRGxDfgo8L5RFCxJjSz3Fskr+prHTvG1j439B0WTDnwzcCgzDwNExF5gKzAY4FuBnf3H9wF/HBFrMnNhYM0MwNGjR0+3Zq3EczNdVzC1Fp5dWH7RBPnWwmTXe3aX/5ZvfnP/wSI13PzmV86/vK7vpn9Z0csN5OWi33STAL8AeGLgeB64bKk1mXk8Io4B3wf858Ca8wF27NjR4CXVvrd0XcAUe7HrAl6TX+66gOV8vuC/5c9febpf4Xzg34dPNgnwNYucG/4R3WTNAeCdwFPASw1eV5Km3Qy98D6w2JNNAnweWD9wvA44ssSa+Yg4EzgH+K/BBZn5PPDlZjVLkvpe1Xmf0CTADwAbIuJC4ElgG3Dd0Jp9wPXAV4D3An8/NP+WJLVs2QDvz7RvBO6n187flZkPR8Qu4KHM3AfcCXwqIg7R67y3jaLYiPg48DPAC/R+Kv1iZj47itcat4j4eXoXgmeBzZn5ULcVnb7ltp9WFhF3Ae8Bns7Mi7qup00RsZ7etuC1wP8Bf5qZu7utqh0R8Qbgi8BZ9PLvvsz8SLdVrdyahQm/4jwoIt5Fr7s/HhEfBcjMD3dcVisiYpbe/yyfAG6uHuD97af/xsD2U2D70PbTsiLiJ4DngE+uwgA/Hzg/M/8pIs4Gvgb87Gr4u4uINcAbM/O5iHgdvbHur2fmgx2XtiKl7sTMzC8MHD5Ib1yzKmTmHEBEdF1KW5psPy0rM78YET/YdR2jkJlP0dtsQGb+T0TM0dtpVv7vrj/afa5/+Lr+nzpd7JBSAT7kA8Cfd12EltRk+6kmXP+H1CXAVzsupTX93w6/BrwV2JOZZb+3iQvwiPg7erO3Ybdm5uf6a24FjgP3jLO209Xke1tFmmwt1QSLiDcBnwVuysz/7rqetmTmS8CPRMS5wF9GxEWZ+a9d17USExfgmXnVqZ6PiOvpXTy6stpOl+W+t1WmyfZTTaj+fPizwD2Z+Rdd1zMKmflsRDwAbAEM8FHr72r4MHB5Zn6363p0Sk22n2oC9S/03QnMZebvd11PmyLiLcCL/fD+HuAqeu/dVFK1XSiH6G3/eaZ/6sHM/JUOS2pNRPwc8Ef07nl/Fvh6Zv50t1Wdnoh4N/CHvLL99Pc6Lqk1EfFp4ArgPOBbwEcy885Oi2pJRPw48CXgIL2dUQC/nZn7u6uqHRHxNuBuev8mzwDuzcxd3Va1cqUCXJL0Cj9STZKKMsAlqSgDXJKKMsAlqSgDXJKKMsAlqSgDXJKKMsAlqaj/B/SERwKZzoH2AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(skullstripped_img_list[sample][seg_img_list[sample] > 0], range=(-2,3.5),bins=200, density=True)\n",
    "plt.hist([skullstripped_img_list[sample][seg_img_list[sample] == i] for i in [1, 2, 3]], range=(-2,3.5),stacked=True,bins=200, density=True)\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'df' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-2-8525a3d7bddd>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;31m#Filter the csv file with only IXI filenames\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0msheet\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvalues\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      3\u001b[0m \u001b[0mfiltered\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msheet\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0msheet\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m5\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'IXI'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'df' is not defined"
     ]
    }
   ],
   "source": [
    "#Filter the csv file with only IXI filenames\n",
    "sheet = df.values\n",
    "filtered = list()\n",
    "for i in range(0, len(sheet)):\n",
    "    if sheet[i, 5][0:3] == 'IXI':\n",
    "        filtered.append((sheet[i, 4], sheet[i, 3]))\n",
    "sheet = filtered\n",
    "\n",
    "#Create the labels for the skullstripped images\n",
    "age = list()\n",
    "for i in range(0, len(file_id)):\n",
    "    for patient in sheet:\n",
    "        if patient[0][3:6] == file_id[i]:\n",
    "            age.append(patient[1])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Next we compute the histogram\n",
    "histograms = np.empty((len(skullstripped_img_list), 200))\n",
    "        \n",
    "for i in range(0, len(skullstripped_img_list)):\n",
    "    histogram = skullstripped_img_list[i][seg_img_list[i] > 0]\n",
    "    histograms[i], _ = np.histogram(np.ravel(histogram), bins=200, density=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Next we normalize the vector for better prediction\n",
    "# ready_made_histogram = np.load('../BrainAge/histograms.npy')\n",
    "# print(np.mean(ready_made_histogram[512]))\n",
    "# print(np.std(ready_made_histogram[512]))\n",
    "# means = np.mean(histograms, axis=1, keepdims=True)\n",
    "# std = np.std(histograms, axis=1, keepdims=True)\n",
    "\n",
    "# histograms = histograms - means + 5/9\n",
    "# histograms = histograms / std * 0.684\n",
    "# print(np.mean(histograms[512]))\n",
    "# print(np.std(histograms[512]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Some other sanity-check to make sure that input and label and properly aligned\n",
    "X = histograms\n",
    "Y = np.array(age)\n",
    "# print(file_id)\n",
    "# print(age)\n",
    "# df.head(590)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'X' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-46-5c25b3cccc88>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     34\u001b[0m \u001b[0;31m# print('Size of the test set : {} patients'.format(len(test)))\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     35\u001b[0m \u001b[0;31m# We use our function on this split\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 36\u001b[0;31m \u001b[0mtrain_linear_regression\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtrain\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtest\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     37\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     38\u001b[0m \u001b[0;31m# 2. Hospital split\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-46-5c25b3cccc88>\u001b[0m in \u001b[0;36mtrain_linear_regression\u001b[0;34m(train, test)\u001b[0m\n\u001b[1;32m      9\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     10\u001b[0m     \u001b[0;31m# Train on the training set : fit the models with pairs (X, Y)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 11\u001b[0;31m     \u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mtrain\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mY\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mtrain\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     12\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     13\u001b[0m     \u001b[0;31m# Predict on the test set\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'X' is not defined"
     ]
    }
   ],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.linear_model import RidgeCV\n",
    "\n",
    "# We create a function to train and test a linear model\n",
    "\n",
    "def train_linear_regression(train, test):\n",
    "    # Create model\n",
    "    model = RidgeCV()\n",
    "\n",
    "    # Train on the training set : fit the models with pairs (X, Y)\n",
    "    model.fit(X[train], Y[train])\n",
    "\n",
    "    # Predict on the test set\n",
    "    predictions = model.predict(X[test])\n",
    "    true_values = Y[test]\n",
    "    \n",
    "    # Plot the results and compute the MAE\n",
    "    plt.scatter(predictions, true_values, s=10)\n",
    "    plt.grid()\n",
    "    plt.xlabel('Predictions')\n",
    "    plt.ylabel('True values')\n",
    "    mae = np.mean(np.abs(predictions - true_values))\n",
    "    plt.title('Mean Absolute Error : {:.2f} years'.format(mae))\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "# We don't perform 5 fold cross validation here, but a single split\n",
    "\n",
    "# 1. Random split\n",
    "# We randomly split the data between train and test with 80/20 proportions\n",
    "train, test = train_test_split(np.arange(sample_size), test_size=0.25,\n",
    "                              random_state=0)\n",
    "# print('Random Split')\n",
    "# print('Size of the test set : {} patients'.format(len(test)))\n",
    "# We use our function on this split\n",
    "train_linear_regression(train, test)\n",
    "\n",
    "# 2. Hospital split\n",
    "# We split the hospitals using Guys hospital as the test set\n",
    "# train = np.where(df.hospital != 'Guys')[0]\n",
    "# test = np.where(df.hospital == 'Guys')[0]\n",
    "# print('Hospital Split')\n",
    "# print('Size of the test set : {} patients'.format(len(test)))\n",
    "# train_linear_regression(train, test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor(281.9251, grad_fn=<MseLossBackward>)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f7eaed3d668>]"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD1CAYAAAC1BoUqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAG2BJREFUeJzt3XtwXOWZ5/Fvt+5X62bLsiWwAfGAIcFkGGOGQMhwCTCZMTOVIrAMmCwFyQYqm2xqdxIqVUzNhBomlcB6axMmIXhjMoTLQBK8wQMBh0DY4eJwv9gP2NhgGcmyLRnLFtate/84p+22fJFsST7SOb9PlUqn33Pp58jyr1+95+3TqWw2i4iIxFM66gJERGTiKORFRGJMIS8iEmMKeRGRGFPIi4jEWGHUBeQzsxLgT4F2YCjickREpooCoAlY5e59+SsmVcgTBPwfoi5CRGSKOgd4Nr9hsoV8O8C9997LzJkzo65FRGRK6Ojo4KqrroIwQ/NNtpAfApg5cybNzc1R1yIiMtXsN8ytC68iIjGmkBcRiTGFvIhIjCnkRURiTCEvIhJjCnkRkRiLTcj/4oUPuPruF6IuQ0RkUolNyG/s7uW5ddvQh6CIiOwVm5CvLS9iMJNlZ99g1KWIiEwasQn5mvJiALb3DkRciYjI5BGbkK8LQ75rV3/ElYiITB6xCfnaiiIAunsV8iIiObEJeQ3XiIjsLzYhr+EaEZH9xSbkq8uKSKVgu4ZrRET2GPF+8mbWAtwDNAJZ4CfuvsTM/h64HtgSbnqzu68I9/k2cB3BvY2/5u6Ph+0XA0sIPqrqp+5+23idSEE6xbSyIro1XCMissdoPjRkEPimu79sZlXAS2b2RLjuDnf/fv7GZjYPuAI4BZgFPGlmJ4arfwhcCLQBq8xsubu/PR4nAsGQjS68iojsNWLIu3s74UdKuXuPma0GZh9il0XA/eGHya43s7XAgnDdWnd/D8DM7g+3HbeQrykvUsiLiOQ5rDF5M5sDnA7kbhJzk5m9bmZLzaw2bJsNbMzbrS1sO1j7uKktL6Z7l4ZrRERyRh3yZlYJPAx83d13AHcCxwPzCXr6P5iQCg9DbUWxLryKiOQZ1Qd5m1kRQcDf6+6/BHD3zXnr7wJ+Ez7cBLTk7d4ctnGI9nFRW15El0JeRGSPEXvyZpYC7gZWu/vtee1NeZv9NfBmuLwcuMLMSsxsLtAKvAisAlrNbK6ZFRNcnF0+PqcRqCkvZvdAht0D+31guYhIIo2mJ382cDXwhpm9GrbdDFxpZvMJplVuAL4M4O5vmdmDBBdUB4Eb3X0IwMxuAh4nmEK51N3fGsdzoTZ8Q1R3bz9N08rG89AiIlPSaGbXPAukDrBqxSH2uRW49QDtKw6131jVhfev6dqlkBcRgRi94xV0/xoRkeFiFfL5wzUiIhK3kM/dblg3KRMRAWIW8jVluZ68hmtERCBmIV9cmKaypFDDNSIioViFPARDNhquEREJxC/ky4s1XCMiEopdyNdVFOvToUREQrEL+fqKErbt7Iu6DBGRSSF2Id9QWczWXf1ks9moSxERiVzsQr6+spj+wQw7+wajLkVEJHLxC/mKEgC27dS4vIhI/EK+MnhD1LZdGpcXEYldyDdUBj35rerJi4jEN+Q1XCMiEsOQr6sIh2s0jVJEJH4hX1yYprq0kG16Q5SISPxCHoIhm63qyYuIxDPk6yuLNSYvIkJcQ76iRFMoRUSIa8irJy8iAsQ25Evo6u1nKKP714hIssUy5Bsqi8lm9YHeIiKxDHndv0ZEJBDPkA/vX6NplCKSdLEM+QaFvIgIENOQ13CNiEggliE/rayIgnRKc+VFJPFiGfLpdIr6imK29qgnLyLJFsuQB5hRXUJnz+6oyxARiVR8Q76qlM07NFwjIskW45AvobNHIS8iyRbrkN+2q4/BoUzUpYiIRKZwpA3MrAW4B2gEssBP3H2JmdUBDwBzgA3A5e7ebWYpYAlwKdALXOvuL4fHWgx8Jzz0d9192fiezl7Tq0vJZmHbrn4aq0sn6mlERCa10fTkB4Fvuvs8YCFwo5nNA74FrHT3VmBl+BjgEqA1/LoBuBMgfFG4BTgTWADcYma143gu+2isCubKd2pcXkQSbMSQd/f2XE/c3XuA1cBsYBGQ64kvAy4LlxcB97h71t2fB2rMrAn4HPCEu3e5ezfwBHDxuJ5Nnhlh710zbEQkyQ5rTN7M5gCnAy8Aje7eHq7qIBjOgeAFYGPebm1h28HaJ8SMXE9eF19FJMFGHfJmVgk8DHzd3Xfkr3P3LMF4/aTRUBmE/OYd6smLSHKNKuTNrIgg4O9191+GzZvDYRjC751h+yagJW/35rDtYO0TorgwTV1FsXryIpJoI4Z8OFvmbmC1u9+et2o5sDhcXgw8ktd+jZmlzGwh8FE4rPM4cJGZ1YYXXC8K2ybMjKoSXXgVkUQbcQolcDZwNfCGmb0att0M3AY8aGbXAe8Dl4frVhBMn1xLMIXySwDu3mVm/wisCrf7B3fvGpezOIjpVSVs0YVXEUmwEUPe3Z8FUgdZff4Bts8CNx7kWEuBpYdT4FjMqCplbefOo/V0IiKTTmzf8QrBTcq29PSR0Qd6i0hCxTrkG6tKGMxk9YHeIpJYsQ753BuidDdKEUmqeIf8njdE6eKriCRTzEM+d2sD9eRFJJniHfLVQU++4yP15EUkmWId8qVFBdRXFNOukBeRhIp1yAM01ZTS/tHHUZchIhKJ2If8zOoy2rerJy8iyRT7kJ+lnryIJFjsQ75pWhk7dg+yq28w6lJERI662If8rJpgGqV68yKSRLEP+ZnVuZDXuLyIJE/sQ35WTRmALr6KSCLFPuQbw578hxquEZEEin3IFxemaags0bteRSSRYh/yEFx8/VAhLyIJlIiQb5pWSvt2DdeISPIkJOTLNFwjIomUkJAvpadvkJ7dA1GXIiJyVCUj5HPTKNWbF5GESUTIzw7f9bpJ4/IikjCJCPnm2nIA2rp6I65EROToSkTIT68sobgwzcZu9eRFJFkSEfLpdIrm2jLautWTF5FkSUTIQzBks7FLPXkRSZbEhHxLbRkb1ZMXkYRJTsjXlbO9d0Bz5UUkUZIT8rkZNrr4KiIJkpiQb64N3hC1UdMoRSRBEhPyLXVBT17TKEUkSRIT8rXlRVQUF2gapYgkSmJCPpVK0VKnaZQikiyFI21gZkuBzwOd7n5q2Pb3wPXAlnCzm919Rbju28B1wBDwNXd/PGy/GFgCFAA/dffbxvdURqY3RIlI0owY8sDPgP8N3DOs/Q53/35+g5nNA64ATgFmAU+a2Ynh6h8CFwJtwCozW+7ub4+h9sPWXFvO8+91kc1mSaVSR/OpRUQiMeJwjbs/A3SN8niLgPvdvc/d1wNrgQXh11p3f8/d+4H7w22Pqpa6cnb2DdLdq7nyIpIMYxmTv8nMXjezpWZWG7bNBjbmbdMWth2s/aia2xDMsFm/ddfRfmoRkUgcacjfCRwPzAfagR+MW0UTaE59BQAbFPIikhCjGZPfj7tvzi2b2V3Ab8KHm4CWvE2bwzYO0X7UtNSVU5BOqScvIolxRCFvZk3u3h4+/GvgzXB5OfALM7ud4MJrK/AikAJazWwuQbhfAfynsRR+JIoK0rTUlrF+m0JeRJJhNFMo7wPOAxrMrA24BTjPzOYDWWAD8GUAd3/LzB4E3gYGgRvdfSg8zk3A4wRTKJe6+1vjfjajMKehQsM1IpIYI4a8u195gOa7D7H9rcCtB2hfAaw4rOomwJz6Cl5cr2mUIpIMiXnHa85x0yvo7R9iS09f1KWIiEy4xIV8boaNLr6KSBIkLuTnNoTTKHXxVUQSIHEhP6umjOKCNO+pJy8iCZC4kC9IpzimvlwzbEQkERIX8hCMy2tMXkSSIJEhf/z0CjZs7WVwKBN1KSIiEyqRId/aWEX/UIYP9HmvIhJzyQz5GZUAvLN5Z8SViIhMrESG/AlhyK/t7Im4EhGRiZXIkK8oKWR2TZl68iISe4kMeYATGyt5t1MhLyLxltiQb22sYt2WnQxlslGXIiIyYZIb8jMq6R/M8L5ubyAiMZbckG+sAtCQjYjEWnJDPpxh8+5mzbARkfhKbMhrho2IJEFiQx7AZlbhHerJi0h8JTrk5zVVs3bLTnYPDEVdiojIhEh0yJ8yq5qhTJZ3NC4vIjGV6JCfN6sagLc+3BFxJSIiEyPRId9SW05VSSFvK+RFJKYSHfLpdIqTm6p568OPoi5FRGRCJDrkIRiyWd3eo9sbiEgsKeRnVfPxwBAbdHsDEYmhxIf8Kbr4KiIxlviQb51RRVFBSuPyIhJLiQ/54sI0JzdV89rG7VGXIiIy7hIf8gDzW2p4ve0jXXwVkdhRyAOnH1NDb/+Q3vkqIrGjkAfmt9QC8KqGbEQkZhTywJz6cmrKi3j1A4W8iMSLQh5IpVLMb6nhlY3dUZciIjKuCkfawMyWAp8HOt391LCtDngAmANsAC53924zSwFLgEuBXuBad3853Gcx8J3wsN9192XjeypjM7+lhqff2ULP7gGqSouiLkdEZFyMpif/M+DiYW3fAla6eyuwMnwMcAnQGn7dANwJe14UbgHOBBYAt5hZ7ViLH0+nH1NLNgtvtGm+vIjEx4gh7+7PAF3DmhcBuZ74MuCyvPZ73D3r7s8DNWbWBHwOeMLdu9y9G3iC/V84IjW/pYZUCv74voZsRCQ+jnRMvtHd28PlDqAxXJ4NbMzbri1sO1j7pDGtrIiTZlbzwvptUZciIjJuxnzh1d2zQCzeRXTm3Dpeer+b/sFM1KWIiIyLIw35zeEwDOH3zrB9E9CSt11z2Haw9kll4XF17B7I8HqbplKKSDwcacgvBxaHy4uBR/LarzGzlJktBD4Kh3UeBy4ys9rwgutFYduksmBuPQAvrB9+CUJEZGoaMeTN7D7guWDR2szsOuA24EIzexe4IHwMsAJ4D1gL3AV8FcDdu4B/BFaFX/8Qtk0qdRXFWGMVz7+ncXkRiYcR58m7+5UHWXX+AbbNAjce5DhLgaWHVV0EzjyujodeamNgKENRgd4rJiJTm1JsmDPn1tPbP8Trmi8vIjGgkB/mz46vJ52CP7y7JepSRETGTCE/TG1FMZ9sDm5xICIy1SnkD+AzJ07ntY3b6d7VH3UpIiJjopA/gM/YdDJZeHbt1qhLEREZE4X8AZzWXMO0siIN2YjIlKeQP4CCdIpzWht4+p0tZLOxuGODiCSUQv4gPmsz2NLTp6mUIjKlKeQP4vyTZ1CYTvHYWx1RlyIicsQU8gdRU17MWcfX89ibHRqyEZEpSyF/CJ87ZSbrt+7i3c6dUZciInJEFPKHcNG8RlIpeOxNDdmIyNSkkD+EGdWl/Mkxtax4o33kjUVEJiGF/Aj+4pNNrOnowTt6oi5FROSwKeRH8JenzaIwneKXr7RFXYqIyGFTyI+gobKEz5w4nV+/somhjGbZiMjUopAfhb/5VDObd/Tx3Dp9YpSITC0K+VE4/+QZVJUW8vDLGrIRkalFIT8KpUUFXDZ/No++0U6Xbj8sIlOIQn6Urj7rWPoHMzz4x41RlyIiMmoK+VE6sbGKM+fWce8L7+sCrIhMGQr5w3DNWXPY2PUxT7/TGXUpIiKjopA/DBed0khjdQl3PbM+6lJEREZFIX8YigrSXH/OcTz33jZe/qA76nJEREakkD9MVy44hpryIn701LqoSxERGZFC/jBVlBSy+Kw5PLl6M2s6dkRdjojIISnkj8C1fzaHqpJCvv/4O1GXIiJySAr5I1BbUcxXzjueJ1dvZtWGrqjLERE5KIX8EfrS2XOYUVXCP61YrY8HFJFJSyF/hMqLC/nGhSfy8gfbWf7ah1GXIyJyQAr5Mbj8jBY+2TyN7z66mo8+Hoi6HBGR/Sjkx6AgneLWyz7Btp19/OC3HnU5IiL7KRzLzma2AegBhoBBdz/DzOqAB4A5wAbgcnfvNrMUsAS4FOgFrnX3l8fy/JPBJ5qncc1Zc1j23AYu/UQTC4+rj7okEZE9xqMn/1l3n+/uZ4SPvwWsdPdWYGX4GOASoDX8ugG4cxyee1L4758zjq0r55sPvqZhGxGZVCZiuGYRsCxcXgZcltd+j7tn3f15oMbMmibg+Y+6ipJC7vjifDp27OY7v35Ts21EZNIYa8hngd+a2UtmdkPY1uju7eFyB9AYLs8G8m/G3ha2xcLpx9TyjQta+b+vfcj/+X8boi5HRAQYe8h/2t0/RTAUc6OZnZu/0t2zBC8EifDV807gonmN3LpiNf+xbmvU5YiIjC3k3X1T+L0T+BWwANicG4YJv+duvr4JaMnbvTlsi410OsUPLj+NuQ0VfPnnL7G6Xfe2EZFoHXHIm1mFmVXlloGLgDeB5cDicLPFwCPh8nLgGjNLmdlC4KO8YZ3YqCot4mdf+lPKiwtYvPRFPtjWG3VJIpJgY+nJNwLPmtlrwIvAo+7+GHAbcKGZvQtcED4GWAG8B6wF7gK+OobnntSaa8v5+XVn0jeY4fIfP8fazp1RlyQiCZWaTDNBzGwOsH7lypU0NzdHXc6YrW7fwdV3v0gmm+We/7yAU2dPi7okEYmhtrY2zj//fIC57r4hf53e8TqBTm6q5t++chZlRQV88cfP8dibHVGXJCIJo5CfYHMbKnjov5zFCY1VfOVfX+J7j61hKDN5/noSkXhTyB8FTdPKePDLC7lyQQs/+v06vvAv/6FxehE5KhTyR0lJYQH/9DefZMkV81m/dReX/q8/8KPfr2X3wFDUpYlIjCnkj7JF82fz22+cy2dtOt97zLng9qdZ/tqHZDSEIyITQCEfgRlVpfz46jP41+vOpLKkkK/d9woX3vE0D6z6gL5B9exFZPwo5CP06dYGHv3aOSy5Yj4lhQX83cNv8Ol/forvPbZGY/YiMi7GdD95GbuCdIpF82fzV6fN4tm1W1n67Hr+5el1/Oj36zitpYZLTp3JBSfP4PjplaRSqajLFZEpRiE/SaRSKc5pnc45rdPp7NnNI698yK9f3cRt/76G2/59DcfUlXP2CfWccWwdC+bW0VxbptAXkREp5CehGVWlXH/ucVx/7nF8uP1jfremk6fWdPLo6+3c92Jwt+bpVSXMa6rmpKYqTp5Zjc2s4tj6csqL9U8qInspESa5WTVl/O3CY/nbhceSyWR5p7OHVeu7eOWD7azp6OG5ddvoH8rs2b6hspiWunJaastpri1jRlUJDVUlNFQGX9MrS6guK9RfASIJoZCfQtLpFCfNrOakmdVcfVbQNjCUYcPWXazp6OGDrl42dvWysbuXVzZ28+gb7Qd8d21RQYrq0iIqSwupKi2kqmTvcnVpEWXFBZQWFlBSlKa0ME1pUW45/3sBRQUpCtIpCtNpCgtSFKaDx0UF6bA9RWFBek97YTqlFxeRo0whP8UVFaRpbayitbFqv3WZTJbu3n627uxn684+tu7sY0tPH9t29bPj4wF6dg/Ss3uAnX2DbOzqpWf3IDt2D7B7YIiBoYmZt18QBn46BSmC7+lUilQquC6x93G4Tfg4t80Bv7Pv43R677FTuWOHz58Kt0+Fz59bkd+WCp83aE/tXQ73Zc+2+x4vd7C964bvv/c5c/vu+9z7Pgf59eTvn1sOn3Pf8xu+z94X1eHr9h4nPOrw4+dts9/Pgn1XHmjbfR8PWz+srlHtM2z93v33/Tnuu++Bjzmaug+2joM9x37nePC6D7RPdWkR57Q2jHtHSCEfY+l0ivrKEuorSzD2fxE4lKFMlr7BIXYPZNg9METf4L7fdw8MMZTJMjCUZSiTZTCTYXDPcpahTCZvXZbBoUzYnmUgk4EsZLJZslnI7FnOkiVYzmQhm80SbJrN2+bQ3/fsFx4r/y+ZbDY4Vja7d3lvO2QzkCVD7sasWdhTU26b3A7ZAxwvf5/859vnOcJ9Oehz5G2f2yZvfbBPNu94ew82/DmG7092WH3Dz2HP+e1bx37nlffzlPGTSsGT/+0zHD+9clyPq5CXAypIpygvLqS8OOpKZKrIf3GD/BeD/V8chr+IkLcu2Hb0x8o1HO4+ue3Zb/tRHOsI6t7/xXHfbStKCplVU8Z4U8iLyLgYPjSSt+ao1yJ76R2vIiIxppAXEYkxhbyISIwp5EVEYkwhLyISYwp5EZEYm2xTKAsAOjo6oq5DRGTKyMvMguHrJlvINwFcddVVUdchIjIVNQHr8hsmW8ivAs4B2gF9Dp6IyOgUEAT8quErUlndgEJEJLZ04VVEJMYm23DNETOzi4ElBH+2/NTdb4u4pHFhZkuBzwOd7n5q2FYHPADMATYAl7t7t5mlCH4GlwK9wLXu/nIUdR8pM2sB7gEaCe7g9BN3XxLncwYws1LgGaCE4P/lQ+5+i5nNBe4H6oGXgKvdvd/MSgh+Tn8CbAO+6O4bIil+DMysAPgjsMndPx/38wUwsw1AD8GQ9KC7nzGRv9+x6MmHvyg/BC4B5gFXmtm8aKsaNz8DLh7W9i1gpbu3AivDxxCcf2v4dQNw51GqcTwNAt9093nAQuDG8N8yzucM0Af8ubufBswHLjazhcA/A3e4+wlAN3BduP11QHfYfke43VT0X4HVeY/jfr45n3X3+e5+Rvh4wn6/YxHywAJgrbu/5+79BD2BRRHXNC7c/Rmga1jzImBZuLwMuCyv/R53z7r780CNmTUdnUrHh7u353oq7t5DEACzifE5A4T17wwfFoVfWeDPgYfC9uHnnft5PAScH/b6pgwzawb+Avhp+DhFjM93BBP2+x2XkJ8NbMx73Ba2xVWju7eHyx0EQxsQs5+Dmc0BTgdeIAHnbGYFZvYq0Ak8QTAVbru7D4ab5J/bnvMO139EMMQxlfxP4H8AuQ8prife55uTBX5rZi+Z2Q1h24T9fscl5BPL3bPkfahPXJhZJfAw8HV335G/Lq7n7O5D7j4faCb46/SkiEuaMGaWu870UtS1RODT7v4pgqGYG83s3PyV4/37HZeQ3wS05D1uDtvianPuT7bwe2fYHoufg5kVEQT8ve7+y7A51uecz923A08BZxH8eZ6bIJF/bnvOO1w/jeCC5FRxNvBX4UXI+wmGaZYQ3/Pdw903hd87gV8RvKBP2O93XEJ+FdBqZnPNrBi4AlgecU0TaTmwOFxeDDyS136NmaXCi3Yf5f0JOCWE46x3A6vd/fa8VbE9ZwAzm25mNeFyGXAhwfWIp4AvhJsNP+/cz+MLwO/CHuCU4O7fdvdmd59D8P/1d+5+FTE93xwzqzCzqtwycBHwJhP4+x2LKZTuPmhmNwGPE0yhXOrub0Vc1rgws/uA84AGM2sDbgFuAx40s+uA94HLw81XEEy1Wksw3epLR73gsTsbuBp4IxyfBriZeJ8zBO9WXBbOFEsDD7r7b8zsbeB+M/su8ArBCyDh95+b2VqCC/NXRFH0BPg74n2+jcCvzAyC/P2Fuz9mZquYoN9vveNVRCTG4jJcIyIiB6CQFxGJMYW8iEiMKeRFRGJMIS8iEmMKeRGRGFPIi4jEmEJeRCTG/j8nlaDh8sRrDgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import torch \n",
    "from torch.autograd import Variable \n",
    "  \n",
    "class LinearRegressionModel(torch.nn.Module): \n",
    "  \n",
    "    def __init__(self): \n",
    "        super(LinearRegressionModel, self).__init__() \n",
    "        self.linear = torch.nn.Linear(200, 1)\n",
    "  \n",
    "    def forward(self, x): \n",
    "        y_pred = self.linear(x) \n",
    "        return y_pred \n",
    "  \n",
    "model = LinearRegressionModel() \n",
    "  \n",
    "criterion = torch.nn.MSELoss() \n",
    "optimizer = torch.optim.SGD(model.parameters(), lr = 0.001) \n",
    "# optimizer = torch.optim.Adam(model.parameters(), lr=0.1, betas=(0.9, 0.999), eps=1e-08, weight_decay=0.1, amsgrad=False)\n",
    "\n",
    "X = X.astype('float32')\n",
    "Y = Y.astype('float32')\n",
    "\n",
    "loss_hist = []\n",
    "for epoch in range(500): \n",
    "  \n",
    "    # Forward pass: Compute predicted y by passing  \n",
    "    # x to the model\n",
    "    inputs = Variable(torch.from_numpy(X))\n",
    "    labels = Variable(torch.from_numpy(Y))\n",
    "    outputs = model(inputs)\n",
    "    loss = criterion(outputs, labels)\n",
    "    loss_hist.append(loss)\n",
    "    optimizer.zero_grad() \n",
    "    \n",
    "    loss.backward() # back props\n",
    "    optimizer.step() # update the parameters\n",
    "print(loss_hist[499])\n",
    "plt.plot(range(0,500),loss_hist)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAETCAYAAADXmaY8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztvX2cHFWV8P/t7rxMSIbEwCQQZhICkxxmxDXEl2SXBIiQjeu6QZ+f65qgi8HVxdXFXVcQd/dBl9VnUdhVUYRdFUQhBHXXh/zUNZuwvCSsmQgBJc5wkyEhM3kbEobA5GUySXc/f1RVp7q6qrqqu2u6q/t+P598Mt1dfeuequpzzz333HMS2WwWjUaj0dQ3yWp3QKPRaDTRo5W9RqPRNABa2Ws0Gk0DoJW9RqPRNABa2Ws0Gk0DoJW9RqPRNABa2WsqiohkRaS9wm0+LiJ/Vsk2NZpGY0y1O6DxRkReAmYAM5RSh2zvPwvMA2YrpV6qQr9mAy8C/6qU+vhon98PEXkceEAp9Z0yvr8QOGV7+zGl1B+V37vKICLjgNXAW4FZwBKl1OMex/0aaFZKtXq0dQXw38Ax29ufUErdb37+AHAlMBE4AHyl1GurqS7asq99dgErrBci8ibgjOp1B4A/BV4F/kRExle5L1HwSaXUJNs/V0UvIgXGktt7foQ93sYm4IMYCtiLG4GDAdra55D3fttn/wScr5Q6E1gOfFFE3lJinytCGdesodEXrfb5AYZy/Yb5+lrg+8AXrQNMhfsl4P3AeOAnwF8rpY6LyBvMNhZg3O+ngOuVUnvM7z4ObATeAfwO8EtgpX0mYUdEEmZ//h74AvBHwI8dh71LRP4KOBO4D/isUipjune+izErOQk8qpT6E7Pd3wO+DswFtgOfUkr9j8v5vwC0K6U+aL4+H2NAHAv8A7AYWCgiXwO+p5T6pIhcZF6/t2Aov/+tlPqhm3x+mFbwA2Zbfw2sF5HvOt8DPiQiHwU+C0zFUMzXK6X2me1kgU8Cf4VxT2aH6YdSagT4mtlW2qOvszEGg08D3w7TvuNcv7W9zJr/LgSecZxvHMbAc7lS6nnzvWnAS8AspdRBEXk3xnN7PtCNcU1+Yx57M/BRYBrQD/ydUuon5mcfNj/bgvHs3S0i38PjWdK4oy372mczcKaIdIhICvgAhnKxcxuGkpwHtAPnAbeYnyUxFO4sYCZwHPim4/srgVUYP7RxwGd8+rMIaAXWAD/EGHycvBfDxTAfuBq4znz/H4H/At5gtvENABGZCvwMuBM4C/gX4GcicpZPPwpQSv0dxsBlWeafFJGJGAp4tSnfB4BviUhnmLZtnIOhwGcBH3N7T0TegWERvx84F9iNcb3svAdjAHbth4j8RkRWlthHMK7t32Lc72JME5EBEdklIl81r5m9L98SkWPAC8B+4OfOBswBaA3GAGOxAkMJHxSRS4B7gT/HuMf/Cqy1zQxfxBioJ2MM2g+IyLm2thYAO4HpGIaN67Ok8UYr+3hgWfdLgR5gr/WBaWl/DMOSH1RKDQH/B0OpoZR6RSn170qpY+ZnXwIud7R/n1Jqu1LqOIYCn+fTl2uB/1RKvYqhQN9pWnB2vmz2pQ/DArXcUCcxFOIMpdSwUmqT+f4fAjuUUj9QSp1SSj2EoVgq4Sd/N/CSUuo+s+1ngX8H/tjnO3eKyGHbv3+0fZYBPq+UOmFeL7f3rgHuVUptVUqdAD4H/K45C7H4J/MauSpjpdTvKKVWlyKwiLwXSFmWcRFewLjf52LM7t6CMdja+/IXQDOGMv4P4IRHW/cDK8xnEuBDGM8uGM/ovyqlupRSadNVdAJjfQSl1I+UUvuUUhml1MPADuDttrb3KaW+Yd7D43g/SxoPtBsnHvwAeBJjuv99x2ctGD78Z0TEei8BpABE5Azgq8A7MawggGYRSSmlLBeA3e97DJjk1gkRmYChJP8MQCn1SxHpw5gZfM12aL/t790Yi8wAN2FYZFtE5FXgn5VS95qf73acbjfGDKVcZgELROSw7b0xnFZCbtzgswh5UCk1XOS9GcBW64VS6oiIvIIhz0vm2/ZrVDFMq/wrwLuCHK+UOsDp+79LRG4CfophgduPSwObROSDwMcxZmHOtrrMGcAVIrIfY5a51vx4FnCtiPyl7SvjMJ8NEflTDJfT+eZnk4Czbcc6r5fXs6TxQCv7GKCU2i0iuzB+wB9xfHwIY6r+RqXU3oIvw98AAixQSh0QkXnAsxgDQljei+GH/5aIWNPmKRjWvl3ZtwGWr3cmsM+U4wCG7xURWQRsEJEnzc9nOc41E/iFSx+Okr9AfY7jc2ca137gCaXUUl/JguOWJtb5Xp48pgI+C9uMzKOdSjAHQ2FuNAf/ccBkETkALAwQvZXFf8Y/BsNn78X9nF44/rFtEOwHvqSU+pLzCyIyC2Nd4Urgl0qptIg8R/4zmne9vJ4lpVSvn3CNjFb28eEjwBuUUkft0Qjmwue3ga+KyCeVUi+LyHnAxUqpdRjT7+PAYdM3/vky+nAtht/172zvnQf8SkTeZC3MATeKSBeGdfYpTLeAiPwxxo95D0Y0TxbDBfJz4Bumj/qHwP+H4cv+qUsfngM+KyIzgdcwXCR2BoALbK9/CtwmIh/itN98HnBEKdUTUv6gPAQ8JCKrMdxu/wfoqmSYrOnrtpThOBFpwnCLbMMYbC1+D2ONZj4ukTkisgTDF96H4fu+DXjE/GwahmvnpxjP0FUYLrkVznZsPIAR7jmE4cax+DbwExHZgLHQegZwBcaMdSLGs3DQPO8q4OIi8ns9SxoPtM8+JiilXlRKPe3x8WeBXmCziLwObMCw5sGwuCdgzAA2424tF8UcQK4EvqaUOmD794zZpn2h9hGMaI3nMBZev2u+/zagS0SOYEzvP6WU2qmUegXDt/43wCsYU/R3u0UEKaXWAw8DvzHP4RwQvg68T0ReFZE7zXWK38dYw9iHYXF+GSNqyYtvisgR279nfI4tQCm1AfjfGGsD+zEs4Q+EaUNEfisi1/idBkMBnwesM/+eZfq0c/cHGAQy5uu02fYREVlstnMJ8D8YM6b/AZ4HbjA/y2K4bCyFegfwV0opyzXjJns/hgsri7FYbr3/NIYl/k2zrV7gw+Zn3cA/Y0SCDQBvwoga88P1WSrynYYmoYuXaDSaSiIi92IsqP59tfuiOY1242g0mophRhz9L4wZg6aG0G4cjUZTEcwQ1W3A7UqpXdXujyYf7cbRaDSaBkBb9hqNRtMA1JzP3gwpextGFINr3g+NRqPRFJDC2An9K3Pndh41p+wxFP3GokdpNBqNxo3FGMn38qhFZb8f4MEHH+TIkSO0t1e0DkbV6e3trTuZoD7lqkeZoD7l0jLBgQMHuOaaa8DUoU5qUdmnAc455xyGhoZobXWtuRBb6lEmqE+56lEmqE+5tEx5uLq/9QKtRqPRNABa2Ws0Gk0DoJW9RqPRNABa2Ws0Gk0DoJW9RqPRNAC1GI2j8WB99wAbdxxk8ZwWlnZOr3Z3NA70/dHUMtqyjwnruwe44aFn+f4vd3PDQ8+yvnug2l3S2ND3R1PraGUfEzbuOMjxk0b47PGTaTbuKCg6pKki+v5oah2t7GPC4jktTBibAmDC2BSL57RUuUcaO/r+aGod7bOPCUs7p3Pniku0T7hG0fdHU+toZR8jlnZO10qkhtH3R1PLaDeORqPRNADastc0DDo0UtPIaGWvaQis0MjjJ9P86Ok93LnikoZS+PaBrjVR7d5oqoF242gagkYOjXTuAdjcd7TaXdJUAa3sNQ1BI4dGOge6rfuOV7lHmmqg3TiahqCRQyMXz2nhR0/v4fjJNBPGppg/Y0K1u6SpAlrZaxqGRguNtPvp7QNda2Kw2l3TVAGt7DWaMqnFKB+3Belbr74YgJ4erewbEe2z12jKoFYToDXygrTGHa3sNZoyqFWl2sgL0hp3tBtHoykD5+JnrSjV0V6QrkVXliafSJW9iPw18GdAFngeWAWcC6wBzgKeAT6klBqJsh8aTVTUcpTPaC1IN/qGtbgQmRtHRM4DbgDeqpS6GEgBHwC+DHxVKdUOvAp8JKo+aDSjwdLO6dx69cUNq+Bq1ZWlySdqn/0YYIKIjAHOAPYD7wB+bH5+P/CeiPug0WgiRK8PxIPI3DhKqb0icgfQBxwH/gvDbXNYKXXKPGwPcF5UfdBoKoH2R/tTy64szWkiU/Yi8gbgamA2cBj4EfDOoN/v7e1l8uTJ9PT0RNTD6jA8PFx3MkF9yjU8PMx9657mtidf5kQ6y8Nb+rj5smksnDmx2l0riyjuVWsCVsxNAYNVieOv1+cvjEwDA/5hv1Eu0F4F7FJKHQQQkf8ALgWmiMgY07pvBfa6fbm9vZ2hoSE6Ojoi7OLo09PTU3cyQX3K1dPTw67hNCfSWQBOpLPsGm5iVczlrNd71egyNTc3+34epc++D1goImeISAK4EugGHgPeZx5zLfBIhH3Q2FjfPcAtj2yrmY0/caCa/mh9vzSVJEqffZeI/BjYCpwCngX+DfgZsEZEvmi+992o+qA5jQ6PK42g/uhK+/XjdL/0mkY8iDTOXin1eeDzjrd3Am+P8ryaQtzC4/QPMxjF4tWjUMxxuV9xGpQaHZ0uoUHQ4XHREUWceVzul46xjw86XUKDoMPjoiOKlAlxuV+1mi5CU4hW9g1Eo+VzHy2KKeZSfdpxuF9xGZQ0WtlrNBXBSzE3gk87DoOSRvvsNZpI0T5tTa2glb1GU4Ry4t3LWWhtpDj7RpK1Wmg3jkbjQ7lumFJ92l7nrceY9kZwddUCWtlrND6UEu/uVMil+LS93D/1qBTjsqcg7mg3jkbjQ1g3TKVq0i6e08K4lPHzHJdKsnhOS936/+OypyDuaMteo/EhrBumEjMBL+o1pl2Hb44OWtlr6oao/Nlh3DDNTWNJJROkM9lQMwGna2bjjoOMpDMAjKQzbNxxkFuvvjiUUoyTf1+Hb0aPduNo6oJKuU+CnsstcmR99wD3btpFOpMllYDrFs0uaSYA3q6NoCUQR/N6aOKBtuw1dcFoLfLZLfEHu/pY9sZzOHvSuAKfejoLQ8Mni7bVP3iMcakkI+lMgVIvx7WhFz01TrSy19QFUfuz13cPsLprN9sHhk4r9EyWnz+/H4AfPb2H6xbNZsLYVNE+WG091fsKI+kM41JJlkgLKxfMylPI5bg26tW/rykdrew1dUGUi3zruwf4xINbcz50N46fTDM0fLJoH+wzA4uRdIa2qWdUfJ2h1Pj+uPj5NeHQyl5TN0S1yGdfLPViwtgUzU1jiypKu3vF/t0glrdTERdTzGGvR9DNTXpAiCda2Ws0HmzuO8pD27fR3DSWMckEpzJGLdpxqSQfvewChoZP0tw0Nvf/vZt2FVWUdvfKuFSSS9vPKnDfuGGfXazZ0s9HL7sg73zXLZrN0PBJ3xmFpaD39BtyOY8N4ufXu13ji1b2Go0L67sHuO3JlzmRzuY2NwEkE/DRyy7gxmWSd/wtj2wLtCBaqntlddfuvFDM7z21K+989zzeSzqLqwK2K+g1W/rJZDOcyhQeG8TPrxd+44sOvdRoXNi44yAn0oYlP5LO5Kz6jEeUTZhdoEHDJ/04OnLaFZRKJjC76rqz1q6gDVncj7UGoj/93Vm+MxO92zWeaMteo3Fh8ZwWHt7Sl2fZj6QzpBLGxiknUe8CXblgVi56x45Mn8RVnefkXDpuCtjpOrIse7dji/n59W7X+KKVvabhCLLAuLRzOjdfNo1dw000N43lye0v89t9r5POwr2bdjGvbUrBd+2KstKLmEs7p3PXNfPzQjYnjE3xmWUXsbRzOvPapniez6mg9/T3s2u4qeS+6d2u8SQyZS8iAjxse+sC4Bbg++b75wMvAe9XSr0aVT80GjthFhgXzpxIa3ZqQahkMV91KYuYQQcgrygc+2e3PFK4+GpX0D2JQVZ1dPj2p9Q+amqXyJS9UkoB8wBEJAXsBX4C3Aw8qpS6TURuNl9/tlLn1Q+kxo+wC4xuoZKpZMLVlRP2HNazGjSSx6IaJRDd2rZk1b+1eDBaC7RXAi8qpXYDVwP3m+/fD7ynUifR+UA0xQi7wGg/fkwyQRJj5+y9m3Z5Pl9BzmF/Vu95vDdU6mKv3DxRpkB2tr26a7f+rcWM0fLZfwB4yPx7ulJqv/n3AcDVJOjt7WXy5Mn09PQEPsnarkN5D+Taru20JgZL73UEDA8Ph5IpLsRFrj39R7l42jgA3jX3TFoTg/T0uD8jw8PDtDYNctOis9m67zgHjpzkV3uPA/7PV2uC3Hfmz5jgeg77s5rOQiph/D8+lWB2k3EtN/cdzbWxcOZENvcd5efbX+e5/cc5mYGHt/Rx82XTWDhzIgCzm4YZn0pwIp3Na8dNrrD3yt52KgH9B1+rqd9aXJ6/MISVaWDAf8CNXNmLyDhgOfA552dKqayIZN2+197eztDQEB0hfIvLswNs2PlsLiph+YK5dHTU1vSyp6cnlExxIQ5yre8e4Cubdueej49d2eb7fNy37ml2DadZPKeNVcsMf/i2h4I9Xx0dsMqnL85n1bkpyt7XDTuPct2iZu7ddCjPpXQinWXXcFPO/97RAa1thW5Mp2uzp6eHPdmprO7aDVB0U9f67gF2Dae5svMc1m3bTzoLfa+dykvgVu3fWhyev7CElam5udn389Gw7P8A2KqUsoadARE5Vym1X0TOBV6u1Il0WJjGjzD+evumqjVb+nM7XSv1fBV7Vp193dB9IFCaBac/383Xvqf/KP/05Eu5MM6nel/hrmvmF10HsGYfYIShLpEW2qaeoX9rMWE0lP0KTrtwANYC1wK3mf8/UsmT6bAwjRdhMkE6N1U9pg6yeecgd664hFuvvjj0udd3D/D1R3fw6pETvGd+KzcuE99n1dnXqzrPoW9wV6A0C3ZL3m2Ae3XweF68vlUcxa0tZ9pme2GWIGkeNLVDpMpeRCYCS4E/t719G/BDEfkIsBt4f5R90Ggswsz87JuqLEpND7C+e4CPP/BMbhfuXY/1AhSkXCjWV79Yevu57Ja8W9rlPf3DrOs9klP441JJmpvGFoRsuuXbL5aDx+8a6Bl3dYlU2SuljgJnOd57BSM6R1OH1PqPOujMz9pU9eS+TN4mplLSA2zccTCn6C02dB/IKXt7CKZdkTr7GqTvq7t251nybmmXexKD3HVNW85nP2HcGO554kXSmWxeWKU1aDjz7Vv9tfpUDJ08rTbQO2g1FSNuP+piA9PCmRNZtayD29cpNnQf4KrOc0qSZ/GcFlZ39eUp/Ks6z8n1wblpq9Rrt757gKd6X8m9tix2Nxntm7Cu/8HTrrl17Pl0rHz7pdxjnTytNtCJ0DQVI8o470oTdE+GVVdWDRzh3k27uH2dco1x92Np53Tu/uBbuPi8yZw3uYlPLGnPWfV2S9zCfu1uX6dY9tUnuH2dKnoeZ959OWcS927a5Svj6q7d2DxVpJIJFs9p8dwrUMo91snTagNt2WsqRpxK4QW1Np3H+aUS9sMtSmZ112427jhUcKyVbO32dSrn31cD+X5+t1mJ8/qfPWk8z+993VNG50wgCVx/+YW5Y9zWN0q5xzpKrjbQyl7jSim+9zj9qIMqLftxViQKeCvPILJ7lTmcOfUM9h4+ntuhO3VifkoGy8/v5UpxXn+AzTsHXWVc3z3AHeteyOvD5dKSt2jstkbgdg63XDxOdJRc9dHKXlNAOb73uPyogw5M9uPsOWzclGexa2YNBs/2vVqg6CeMTXFhy0T6Bo8BxmAydeJE9h4ezh1j+fn9ZiXO6+8m4+a+o7kNW/bzr1wwK/C1K9V/r6keWtlrCmiUBbUwkTnWcV7hj8WumV0xJhztnze5iS+YsfuWJT4ulSRLgllTJ3Aqk+U9l7TmrO4wrhQ3GbfuO56n6GX6pFyq5DA0ynNSL+gFWk0BekHNm6Wd7lWmil0zu2K0B2GOSyX5gq29hRdM5U3nnUkmm2Xb3tfYPXicg0MjzGubkkuABhStKOXH/BkT8vpaiqIPIrOmttCWvaaAOPnea4WlndO5btFszxDNxXNaeLCrL+fzh3yL2istARihj19Yu42DQyOMpDM5l0mpO3m37jte8uYop8z6OYkPWtlrXBlt33utb8YqhhWiefxkmr7BwkpWSzunc/3lF+aieZwWtTMtQYL8GYDdd1/OTt7cRqneI1zaflbxLxUhLms0Gq3s64J6UJRxX+gL4r++cZl4+vztfngwFH0SmHLGWAaP5Rc493OZOJ8Frzw5znw/cbvemvBoZR9z6lFR3rHuBSDYVnwvRnsADLpo6mUJW26gBze/xOHjpwDIAG9um5K3aFssAZozL469ApY9T46FXlhtHLSyjzn1EBHhtGrVwBFueOjZkgeuUgdAtwHC2vwE3nnfre+V4we3u4EsrHDIlQtmBRq4/NIi2/Pk/Nujv+XXB07k8v04k6DFfaaocUcr+xqh1B9YXHat+slnLfTdse4F1MARoLyBK0gYpJtSLyXvu/17E8amSh6gnLVuneGQlhL228C0eE4La7b0M5LOMC6VzEuLbD0bSzun05oYZE92qmv9W+dsoJg8emCIDzr0sgYop3aupShLDcMbDYLIt7RzOp9ZdlHZoXz2tLxghDb2Dx7LndOrL24DxNZ97nnf7VQqH5AzjNEZDhn0Gclks7n/57VN8Xw2rBDS7n2v+c4G/OTRNZ/jhVb2NUC5CsMr9rtWCCpfuQOXpXweU0b7F583GYDH1MGcMvLqi1vM+PwZExiTzN8C1dyUn8KgUrHmluxLpIWFF0wt+DzINVzdtTuXWfNUJsudj273tbrdsmRe1XlOYHnilPhOo5V9TVDvm1PCyGcfuCy3RVCL0RltcvJUOmeZW8rIqy/OgQaMnabnTWnKO0f3vtcK+lvJmdVTva/wmDrIJx7cmid3Kc/Ib/e+7mt1O7NkXtp+Fjcuk8DyOPtk+f61hV+baJ99DVDvm1OCyOcWMhh2kdWvlJ/dZ+3VF7ecL0GsoVJize3ygmGV/3rP4ZzyHUlnWN21O89nX+warlwwK1doJZmAjCNHvdtGL/v1WrlgVigfvFfeoLhGhdU7RZW9iHwS+L5S6nUR+VfgEuBzSqlHI+9dA1Hvm1P85HNT7KVEGbkpRCuu3SriUawvkD9DyDg+65wxOaDE3tjlXbOln0w2W1DJyo1i/V7aOZ27rplfNGnb5r6jPLTdWOh1ZrAMO8BafbrlkW2xjwqrd4JY9h9TSn1TRH4fmA58FLgXeEukPdM0DG6K3S3KKIjV6VSI1t9hlJhXWmOAoeGTnt8rRV5n9kuLcalk4CyUdoolbbt9neLuxwbIQEHahXIUdq1FhekooUKCzFKtJ/1dwA+UUr8O+D2NJhB2368VPQMU+NBLjfwIu5BozRD+SM7k+ssvrPh6il3eMclEXhbMMckES6SlIMQzKPZ1DufC/fruAe5+vDc3W3Fei3LWjmopKkxHCbkTxLL/tYj8HJgL/K2ITCI/bYdGUxaWoljdtTu3QGlt4w9rdQap4BS0ulJrYpCODu8UB2Gx982Sd+OOQ2QxcuG88bzJfOrKOYFi6r3at2YwD27ezfVXtOcVI1ndtRu7tygJedei3LUj+6yimpZ1PWw0jIIgyn4VhsumVyl1TETOBj4SpHERmQJ8B7gYY4C4DlDAw8D5wEvA+5VSr4bueQPQSFPRpZ3TzeiQ08Wu7T/SIAo7aAWncpRYqbj17eCRkZyvPgu0TBpX8uI0FCZTu+eJFwsSstl543lnum5wi0LW0Xx+a82lVCsUdccopdLABcDHzbcmBPmeydeBXyilLgLeDPQANwOPKqXmAI+arzUOGnEq6udG8HIT2N0Wfu6aYnsRwoZ5hsXZt9Vdu/nt3vwwzkNHTrge63Q72ftq/3vxnJa8H2Y6k+Wv1mzl3d/YxPruAVYumJW32eyGK+cW9LMS16Ha8fe15FKqJYJE43wTGAtcBnwJOArcA7ytyPcmm9/5MIBSagQYEZGrgSvMw+4HHgc+W0rn65lGnIoWs8CdVqdb4i8r0VcYiy4qS9Q+M3Nam1DoCz170vi8HcBW7hqv8odrtvQD5HLcL7loGllHbuSjIxm27X2Njz/wDHd/8C3cdc18/u3R3zJp0iTf6+DmBgoqb3PT2JLuQyWp9+i2Ugjixvk9pdR8EXkWQCk1KCLjAnxvNnAQuE9E3gw8A3wKmK6U2m8ecwAjwqeA3t5eJk+eTE9PT4BTxYfh4eFAMs1uGmZ8KsGJdJbxqQSzm4J9r1oElasYrQlYMTcFDNLTM+h77NquQ3kDYt++l7lp0dls3Xec+TMm0Joo3oZbO2u7ttOaGCwq0+a+o7lzLZw5seCz2558mRPpLA9v6ePmy6bl9Q1g4w44Za6WphIwfdwIn3zwGU6ks4xNwtvOm8C75p6ZJ4e9r/ZInuMn0/z8+f14cSqTZW3XdubPmMBvDhznRPo4/9N7iJsvm4Y6NMzm/mOMSSby3UCP93I2QwWyeV0LS97xqQRvb51A/2snWdh2RuD7UA6Vev5qibAyDQz4z8aCKPuTIpLEtBdE5CwKw4+92p4P/KVSqktEvo7DZaOUyoqI62Jve3s7Q0NDdHR0BDhVfOjp6QkkU0cHtLbFx2cfVK5Ksb57gCMM5VnAyxfMZWnndFaFbGt5doANO08nM1u+YC4dHdN9ZVrfPZAr2r1h51HuXNGWd48e2r6NE2a5qRPpLLuGm7j16ovz+tbaNsCdj25n8OgI77mklaHhk5xIG66dkxnomDmNVcvyq1HZ+2q5ZJybqNwYk0ywfMFcNu44yAkz39qJdJb/2H6c5/e+7vqddBZ2DTexKsB9dcr7P33HSGdh4OgRrrxkTuTP72g/f6NBWJmam5t9Pw+i7O8C/h1oEZF/AN4P/EOA7+0B9iiluszXP8ZQ9gMicq5Sar+InAu8HKCthkRPRd3Jq7iUSrJEWjzTDxdrx0pfHDY9sZsP3v69IIuEz/Ufpnvf66SzcO+mXa5uKLtrxOqfcyOUFdVjJUGzGJNMcOaEMcyYcgafunIOAP2DxxibNAaTCWNTDB4d8ZQxlUwEdsF47U1oFBdkHCiq7JVS3xeRZ4BWWWd7AAAgAElEQVSrMCLE/lgptS3A9w6ISL+IiFJKAVcC3ea/a4HbzP8fKUcATePh3JTUNvWMkhT9Jx7c6pu+2A97OmHr+1Zse5D89uu7B7jniRdztWbt+eYt5W6FotrdNdaawuI5LbmBCnDdgXsqk+XdvzODW6++OG+AHJMkN0A+13+Yux7rLfhuErj+8gsDXw+v1Ak6GqZ2CLJAOwN4FfiR/T2l1L4A7f8l8KDp49+JEcaZBH4oIh8BdmPMFBqORgqrrDSVCK1zJgGz0heHUW6Xtp+Vy7BpT38cJL/9xh0H83bmphLkPQtWG06sWYR9EBiTTOTcWXbXjv3a2AfIUxlyA6R1PnuFLIDLpSXU4qx1Tfx272qqSxA3zqOcXt+fALQBLwJFnwSl1HPAW10+ujJoB+uRaschx51KJI5zWubjUsnQg8bKBbNyJQOtrI93rHshUBSV0+1ht6KdhUzsWJE89oHqVCbLEjmbtqln5Cl3r41l41P57hmrNq59kColVYMde1I5XQWrNgjixslbIRCRtwN/FlmPGoAowipr9UcUVb/KXc+wkoYVKzlYrA0314VFsVq0XgNWc9NYUgljgdSqO9s5Y3LOJQTkWfZWHh23nEBu55vdNOz7eaXuVbGauNrIGV1CpzhWSm0Rke9E0ZlGodI7/Gp1plCr/bKoxAK4W9ZHKCwrGPT867sH+PaTO0lnIZmAj152gas7JehA5Rxsl3YaUUZug3DY6+FM0+xsr1hNXL1wO7oE8dnfYHuZxEidUP/bOSOk0lZUrW7AqtV+RYFzAC+m6L1Y3bU7Z7FnsoXFUiyKKWYr0siaAdgH2819R3Nho6UOwn6bu6z2gtQX0IweQSx7+x05BWzAtlirKY1KhlXWai6QWupX1G6uKNwgpWJXxBb2wXbrvuNlD8JeaZrt7bldE71wWz2C+Oz/92h0RFM6taRo7NRKv0bLnVSJAdxebarUnPZuC7z2wXb+jAls2Hm0rEHYPpB7RQCBe32BWnk+Gw1PZS8iP8EnlbFS6n9F0iNNSYz2jyiopVwLP+5quJP8ro9bCUb7a6vaVJDvumGPNBqTTLB4ztl5fv2FMydy54o21wpeQXEO5FDos68EtRp4EEf8LPtvjlovNLHCzVJuTRT/XrUYbXeS30wiaISKm2ILOkN5rv8wp0zXSjKRKFjA3dx3lF3D6bLrxrpZ7ZWk1hf444anstc1ZjVeuFnKRvKy2mS03Ul+M4lyIlTcUjQ4ZbJ25lpedOdmsfXdA7mEZVZ4p729Sl2bSljkjbTAPxoUzUsvIheKyBoR+Y2IbLf+jUbnNLVJOeXrnLjlT48it/zSTv989pXE7/o4P7twWjOpZML1WL92x6WSPNX7SkG9A6+duRZGIjTjcyu808JK+VAO67sHWHXfFj7x4NayazFU8jnTBIvG+R7wReAO4A8wUh7osoQNjJulXEoKW7dpOoQrDl6rLLxgKlAYA++2ESudMazs6xbN9nTf2MsZbtxxkP7BY7lUDXar1+6vTybg+ivafTdsTWpKMXjUKKI+ks5wx7oXcv0MS7EoIDd5iq331MICf70QpOLUGUqpdQBKqReVUn+PofQ1DUwlLGW3aXq1qxyVi6XwrDq6bljXbmj4ZF7++KHhk57tWVYywK1XX8zKBbOKWr1jkknmtU3Ja+veTbsMix7IZLM5RW+hBo6UbI0XiwJyk6fYeUZzRlbvBFH2J8x89i+KyPUi8keAf+JkjSYAbtP0uE/dgwxWlpvKqugEhpXdP3isQPl5tWdZvc7Se/YEb/bkbAB3Pro911aG/EyZUyacnuSXOsg63UxLpKVgZhb3wTzOBHHj/DUwEbgBoyzhmRiFwzWasvCapsd16l6spKB1jD3h2HWLZtO97zWe6n0lNxuwK8jmprF537e/tkfDFCsJePs6lVekJJkwLH+rn9csPL/ktMRON5NfGoda2mjXaARR9seVUkPAEPChiPujaTDcwgyrEZtfTvSIMzWBX0EVp2U7NHyStqlnMJIu9L9DYboEt/QJbgOIM4/+hu4Ded85d3ITX1h+cdm7W91CSa1MoM6BC7QfvpoEqlQlIm/ASJHwsFLqhYj7pNGMGn45ZNyOg3yL1W1R0q+gipdlG9TaPXTkRMF7bgPIrVfnlzO8qvMc1MDpIiXvuaS1IrtbSwklrYWNdo1IkHQJi0XkPOBPgPvNQiQPK6Vui7x3mpqkXnY1Bo0euX2d4u7HenOx6/aqVkEWJe2EdV2tXDCLjTsO5fzr6sCRXEUsiyCuEStz5obuA1wyLRW6MIkXOtlZfAiU4lgptRf4FxH5T+BzwD9ilBXUNBj1tKsxaPSIfZMS5G9UcuaIubT9rKK58cO4roxznF1QEcuZZiGIa+TGZcKNy4Senp5il8YVr7TIOtlZPAiS4ngOhlX/PuAI8DDw2Yj7palRvHZxzm4apqOjyJdrjCCK2rlJCfKrWo2GD9pZEWvxnBZuX6eMGraZbG7QdbpuKonfIB9FsrN6mT3WEkEs+9XAGmC5Uqov4v5oahyngjR83QcZn0rQ2jYQqx9mEEVtlzeZgDfOOJMbrpw7qpkc3ZKO3fN4b16qg6hTCZSTuiCs4q6n2WMtEcRn/7bR6IgmHtgVj30X54l0ti5zl9RK9Ih9QLnlkW05RQ+QSiYi942XGjJZiuLWOXGiIXRZwjCIyEsYIZtp4JRS6q0iMhXDFXQ+8BLwfqXUq1H2Q1NZLMWzvnsg515wFrGOA0EVUTmWexTuiLxi5QnyipWXi1d/3Qa9ILKVorh1LH40RKrsTZYopQ7ZXt8MPKqUuk1EbjZf6zWAGGJXAG5FrGudUhRRGOUdlTsiqtlGsf46N3EFka0UxV0rs6l6I7CyF5HxSqnCIN/wXA1cYf59P/A4WtnniNvClKUA7BEecZEhrCIKq7yjdEdEsU4Qpr9Bjy1VcetY/MoTJMXx20XkeWCH+frNIvKNgO1ngf8SkWdE5GPme9OVUvvNvw8A+o6ahE0SVYvESQav/DJehM3rErc8P4vntORKDNojjryODSqbTmZWGwSx7O8E3g38XwCl1K9FZEnA9hcppfaKyDRgvYjk7b5VSmVFxDVdcm9vL5MnTy45JrhWGR4e9pRpbdehPGWytms7rYnwqYOrgSVX3GRoTWAWXhksSNPsvFezm4YZn0pwIp1lfCrB7Cbj8819R9m67zjzZ0xg4cyJeW3ftOjs3GeticJzeOHVZiXwegb39B8lmzV2FGSzGfb099Pjce/8ZIuy7174/a7iSliZBgb8Dasgyj6plNotkrfjLu11sB1zMxZKqZfNmrZvBwZE5Fyl1H4RORd42e277e3tDA0N0RG34O0i9PT0eMq0PDvAhp2nc5wsXzCXjo54WEOWXHGWwYnzXnV0QGtbYf3Yr2zabaQK2HmUO1e05VmwHR1GAYgwuLUJlavx6vUMPrR9GyfN3WMnM7BruIlVPr8/N9mKXY+o8PtdlUu13JJhZWpu9k9GHETZ94vI24GsiKSAvwSKVqoSkYkYA8WQ+ffvA7cCa4FrMXbgXgs8EqAPDUE9LEzVgwx+OH3Jfr7rUpWE28Y1K+opyrjzIGsYxWSqt7DJeor5D5LP/uPAp4GZwACw0HyvGNOBTSLya2AL8DOl1C8wlPxSEdkBXIVOu5BHPfg360GGoHj5rstZu3C2CXiuFVSyhOPSzulct2g2Mn2Sa9WsIDLFbZ2iGPWUfz/IpqqXgQ+EbVgptRN4s8v7rwBXhm1Po6kEdssUyneNeM1knEoiTLk/tx2zznQJliylWJ2b+47y0PZtBXJblayOn0zTN7iLeW1TAs9iil2PuFJPMf9BcuN8G5eas0qpj7kcrtHULHbluGZLP4BvWuOguIUJ2pUEnC735zyP3yYm++sgA0rQfQK3PfkyJ9LZArmLtRdU8dVT2GQ9DV5B3DgbgEfNf08B04BKxNtrNKOKXZmNpDO58n1RTM8tJSHTJ+Xec3PBBHX1uLnGSnGZbNxxkBNmrgVnf4q1FzZU1U4l3U2jTb24JYO4cR62vxaRHwCbIuuRRhMRziRugGf5wEpgKQd7FSn7edwWYsMolFKszsVzWnh4Sx8n0tmC/gRprxSr3c/dFFWkS1w29o0mpaRLmI3eCKWJIW6+8KgVgp8CXTynhTVb+nMzjKd6XykoTBKkfSCvEHmx42++bBq7hptc5XZT5uUqTi/3UFSRLvUUQVNJgvjsX+W0zz4JDGLks9FoYodTmY2GEvCyhpd2TufS9rMKCpM813+YDd0HuKrznKIVpcIoNktpz24yBpogA0Sx9r0GAvv7Xr7+qMI06y38s1L4KnsRSWBE1Ow138oopVx3vGo0jUYlXAXOwiSHjozw/V8atW6tmrF+Cj+oYrMr7bFJSCQOBlqcLraPwG0gcHvfbXZTSm6iINe7niJoKomvsjfTGfxcKRVdCRyNJma4FSm/btFshoZP5imiIMrJ6eaxQjQtNnQf8FX2QRWbXWkbu2TzF6e9+ufXvtdA4Pa+2wJnkDUC6xo2N43NhYUWG6DqKYKmkgTx2T8nIpcopZ6NvDeamkUveBl4FSl3lggEQrk/rM+e6z+cs+gBruo8x7c/forNy5UyNglZEpzKZIsmPCu25uA2EISxrP0WfO3XOpUgVGWuSoZ/1suz76nsRWSMUuoUcAnwKxF5ETgKJICsUmr+KPVR48JoPoB6wes0dqvVIpUgV6fWHs4Y1v0Bp102QX32kF9M5pZHtuWUq5crZde+gzzVdyywzH5rDm4DQaUsa/u1TmeNilzpTGEUUZTU07PvZ9lvAeYDy0epL5qAjPYDqBe8TuNWpLxzxuSci8GuiMK6PyxuXCaBlLwd5zOx8IKprq4UgAd+udt04pxeFA4bM+82K7FTCcvaOUNwc5VFTT09+37KPgGglHpxlPqiCchoP4B6wes0XlbrvLYpBe8FdX80N43NWeSlhj06nwkwNka5RcBkbN8LW792NA2NWvC919Oz76fsW0Tk014fKqX+JYL+aAIw2g9gLfzooiaMW8wttt3Nkg3i/vBbeHQqVr9F4OamsXnKfeWCWaxcMMs1AsbaVFVK/drRNjQq6Xsv9fz18uz7KfsUMAnTwtfUDtV4AKv9o4uSsNZqJaxb63re8sg2T+XpVKz3PN5LOovrIrCbm8MtNUGxTVXFqCdLNyj18uz7Kfv9SqlbR60nmlDUywPoxWguQIe1Vitp3fopT/tn1uKk/ZzW39b/Q8Mnc355vwFp4cyJvkVJ/KiGoVEv0TDVxi8RmrboNVVhtOvYFksA5kziVcmc7ZbydEsuZv/s+ssvzOXzSQLNTWN9++GXh31z39GykpIt7Ry9xGBxqmlc6/hZ9jrnvKZsSrHKquEX9otVd7OQK2nd+s3S7GGVmaxh2WeAbz+5k7uumR86Bt4vxXEtUk/RMNXGU9krpWq3SrSm5nAriFGqbzuqbfR+eClcL2Uz2m60jTsOcipzOlOJFTLpZWF7DUhuKY5rWXk24hpBVJSS9VKjycPLWizVKnMqKsA3NDHKUEA/C7nYrtVK98OeIbPYzldwH8Cam8b6vh4tgl6neoqGqTZBipdoGpggRSe8CmKU49u2/MKAr8/WzzddCZZ2FtZl9fIjB/Uvl1LIY2nndO66Zj5LpIUl0sJd18wvqvjczjM0fDLvGOfrcvsZhLB++NFcI6hntGWvAdwtraBWs1dBjEpYZcVmB84drf2Dx0LnhPfDrS5r0ARg9qRm9llKKTMR6/6sXDAr1x54pyf2unde98p5/0sJRw16n6vlh2/0qB6t7DWeP+ygP0q/2O1yfdvFfLbWgGJloXxMHWTzzsGKuXPcrkGQBGBg1J39xINbgdO1bt3SGASx0MPWzvWa8WzccZD3dp7JuElTfBW73/eDLmJ7UQ0/fD3luCmVyN04IpISkWdF5Kfm69ki0iUivSLysIiMi7oPGn+8fthh3DALZ06MZKrtF5poP6Zt6hmR1JR1uwZefbLet9eddda6tdqxt1cM+/0JWjvX2e/mprE518lPul8vWLS13/87H91O/+CxXKin8/tO10tYV1qQe1pponb3xYHR8Nl/Cuixvf4y8FWlVDvwKvCRUeiDxgcvpV6NH6UbQXy2lYx9d57bS7F75Wj/zLKLcn0Zl0rmKc2VC2aFvqZ22Zztecnp7Hf3vtdyyu5EOpun7BbPacm1CfD83tdz1bOWSAt3rriEoeGTnsqylGs/2n74qJ6POBGpG0dEWoE/BL4EfNqsfPUOYKV5yP3AF4C7o+yHxh8/33pcdupGGbUR9hq4RROt7tpd0faCRrJYbpqnel/JvT82SZ6yW9qZXx7RYiSdoW3qGblzeLleiu1TqAU/uY7qid5n/zXgJqDZfH0WcNjMkw+wBzgv4j5oAhAXpe5HLclg78v67oFc6cFS1xOcsoX5/sYdB3OuH4B5504o+L69PKKF32I75IfDul37WvOT19LzUQ0iU/Yi8m7gZaXUMyJyRdjv9/b2MnnyZHp6eoofHCOGh4frTiaoT7mGh4e5b93TbN13nPkzJrBw5sSS2lnbdSjPBbK2azutidHbszi7aZjxqQQn0lnGpxLMOjPFX963MU+m1gTctOhstu47zsRxCY6OZJk/YwKtiUF6egZzx6yYm2JzX39uX8XDW/q4+bJpuXY29x3NXa+t+46Pmtz1+vyFkWlgwD+ENUrL/lJguYi8C2gCzgS+DkyxVcFq5XQx8zza29sZGhqio8SETbVKT09P3ckE9SnXfeue5iubDEW9YedR7lzRVpJluDw7wIadp7NTLl8wl46O8O3Y0xlb2S3htIto5YJZebMJywpftayD1rbT3/3Oky9yIp1lw86jXLeomaHhk2abTSxf0AYYs4HWthbXfj60fVtuX8WJdJZdw02s6uhgffcAX9m0O3e9rls0mwk7j5YtdxCieP6q7YIKK1Nzc7Pv55Epe6XU54DPAZiW/WeUUteIyI+A9wFrgGuBR6Lqg0ZTDk7LtNR48LD+4mJ7HizWbOknk83m0ig81fsKd11jVAt1c58s7TRSKts3wFlpk+1tgn9op1fopDPiZWj4ZGz95LXmgqoE1dhB+1mMxdpeDB/+d6vQB42mKPNnTKhodssg0Sdeu0vtitRiJJ0pyJezums3d6x7wTdyZnzKSGibSibyFL3VRrHQTmvwWiItLLxgal7bbmGqcdz9Wo+hmqOyqUop9TjwuPn3TuDto3FejaYcFs6cyJ0r2kbVMvXayObcsAUwJmkobUvhj0kmeKr3lbzFWLfIGWsDnL1KloUVgjmSzhQd4KxzWTOKeop4qccEbHoHraZuiTIbZlT4KZmFF0zl4JERXtj/OqcyWZKJBB+97AK6972WO8YZPmnl87FjL15ipX9wrgMUu26ru3bnBhVrRuEVlRNH6mngstDKXlOXxNXn6qZk7LKkEuRcLyPpDEPDJ7lvlTFRXt89wJPbD+a5ZoaGT/oOevZYfGdN3WpT7QXSWrkOlUJnvdTUJXH2uTr93HZZ0lnD1w7uLprrr2jPfZ5KwKEjI0UzTJZSDWrlglk5N9KYZCKXoK1S6ApVlUcre01dYl8sTCWql7fdSSlpg53pEjrOPTOXxsBped64TLj+8gtzi6/rtu33HfTWdw/4Luj6kUwk8v4vFbdrEmawXt89wLc2H9IDQhG0stfUJUs7jTz0ltK7d9OuspRBJXK7l2qt2qNfALbtfY3NO703Jw0Nn8wVJ/ebCVj9UQNHcu+FSc5m99kXU8Ze187rmgTNZWN9//9Xr3PDQ89y+zoVSQ7+ekAre02OqIpVVAu70ivHlWNXSNc/8Ay3r1MltVOOaylMZs/mprF5Cv76yy90TbzmDOeU6ZNKSs4WRBmHLT4TNAmf8/v3PN4ba9dPlL9Brew1QH36SMMoJL8fWJ7PPJPlnsd7S7o+fv0J8iMPIo9VbCWdyZJKGNE4Ny4T11h3Z3ufWXZRqAXJhRdM9XQnWRQb4PxkCpvt1L5vIG7rNFD4G9zcd7Si7etoHA1QvepBURIkfC5I1M7iOS082NWX5xop5fp49Sdo5FAQeZyLuW5lB+3F4UsJL7T310rb7EXQ4jOlRt1Y31/btZ2ZM6bl9g3EMTbe+Rvcuu84qyrYvlb2GqA+N5FA8fC5IIPc0s7pXH/5hbnUAuVcH7f+rO7aHXigLSZPsfvoVhzeqvUblDCGQRBlXm6I49LO6bQmBunokNy+gTjGxjvv3fwZEyravlb2GqA+N5EEIeggd+OyaBSJM9f8uFSy7LQMfvfRrTh8WFm8rplXXPxoxqvHOTbeee8qnSFUK3tNjjj8UCq90cZrE9PGHQeZ3TSMPelgJa6Ps//OXPOXtp9V9jm8+rm+e4D+wWOMScKpDIEUtVf7fhu/4rSJLWrCPq/2e2ellq4UWtlrYkNUCsX+A7OfY3wqQWvbQEWteGf/nVZyKZuTnAqlWNbMsUmj3KCVErmU6+ocUEpd8yl38PYamKM4V9i2a20A1MpeU5Rqb1u3CKtQSum3/RxWrdZylJBfYe+NOw5y69UXl+U+cyqUJRdNY922/aSz5CkY+7lPZsgrN1iJxflS1nzKVYZhBuYoFa9X27UW9KBDLzW+hAnJjDpOP0zR6FJDSe3nGJ9KlOQ/D7tRKEiIods5bnlkW8Hi7s+f3+8afugnVyWKcTvj4oGiz0K5KS3cBmavZzDK9BlebddakXNt2Wt8CWqdbO47mqtSFNWUNcwicqlWlf0cs5uGS5LB69yVWgS3W5LjUknGpZKMpDNGnHkmP0G9lSbCT65K9cuSMagV7TYbCDMbWzynhTVb+hlJZ0ja8gA5z2utVVjXqdKK12tWU2tBD1rZa3wJOj2vVFWnYgRdJC0nlNQ6R6k1Tf3OXYlFXvtgMpLOsERaaJt6Bs1NY7nniRfzFL49zt5PrjD9ciuPGMZ/b1foziLmYV0tmWzW/B9+sW0/GZdZjX1gtK9VVAo/pV5LQQ9a2Wt8CWqdzJ8xgQ22eqPVnrJW06oKc+5S1hXcFnXt37XvB2huGsstj2yr2DVwK4/oVMx+g52b1W/F+d/yyLZQBsPGHQfzKnVlsuRSQFvndQ6M9rWKSlJLSt0Lrew1RQnyIEdd1akUpWhfgLS/rjRufQtyzUpZNLTOdd2i2a5WtX0/gL0SVaVca858OlComP0GOz+rP+xszO7GAWOPwkcvu6DgutTjZsFS0MpeUzGism5KjaRY3z3AJx7cykg6w5ot/bnSeVH1bc2Wfi5tPyuwm6CU6CJ7mgL7dXAOOEs7pxe1lMudVdhxppD2ehaKubi8BgmvAfWua+azums3R44c4WNXvrHgnLXmN68mOhpHU/OUGknhVjovyr6NpDM8pg4Gjv4JG63hdR1KSRNcbrplmT4p7323HDzWeezRMc6oHTfl7IxM8uvr0s7p3Lfq7dx61bmeiryUaKd6RCt7Tc1TayFsdux9swg6IBVTfJCvLL2ug3MQsAY1v/bLCUVc2jmdzyy7KFAGTjclHVb5Rhk22UhE5sYRkSbgSWC8eZ4fK6U+LyKzgTXAWcAzwIeUUiNR9UMTf0qdiq9cMIunel9hJJ1hXCpZ8dJ59r6t7tqdO1eYAcnP9eXmvnK7Dk4XysYdRtUmq+2w7pQwcgfNwOnlpgriSiq24Bt0B22jE6XP/gTwDqXUEREZC2wSkf8EPg18VSm1RkTuAT4C3B1hPzR1QCnrAZZPN2p/rT2+vJLn8tpx62zb6UI5lfHe+bu+e4C1XYdYvmBqyb5su5xWJM367oHcjMJasyg2oNy+TuUih0pJ7Rxlaot6JDJlr5TKAla9s7HmvyzwDmCl+f79wBfQyl4TEaMZEleJc9kVaVDre/GcFlZ39eXCEMckE/QPHstZ9/a2LeW4YeezJaU3trfx4ObdXH9FO/PapuQWwgGe6n0ltxjut+B6zxMv5u32vWPdC4B71JTbtQ2S2qJWUn3UApFG44hICsNV0w7cBbwIHFZKnTIP2QOcF2Uf6hX9EMcXr01JQd02bhhFv7NYtb8fUwfZvHMwL0SzErla7G2ks3DPEy9y2Zyz8zJ3WjVp/dxIG3ccLNjtqwaOcMNDz4YqjWgNhm6pLWotEVm1iVTZK6XSwDwRmQL8BLgo6Hd7e3uZPHlyybsYa5Xh4eGyZdrcdzRXgOLhLX3cfNk0Fs6cWKEelkYl5Kokm/uOsnXfcebPmFDytYlCJvu9s7DuoXMX8tqu7fzFwrNZMTcFDHqmvF3bdSinbLNZOJU9na/ecpM8vKWP93aeyfhUghPpLONTCWY35csX5JrNbhomCViqPZ3JcuTIEcYmjQRrAGOSFLTt1o7VFzuW3EFyubcm4KZFZ7N133EubhlDa2KQ+9b152Rwu56VzhEPlXnW3Aj7/A0M+EdUjUqcvVLqsIg8BvwuMEVExpjWfSuw1+077e3tDA0N0VFnqy49PT2uMoWx1B/avi33IzmRzrJruIlVVb5OXnJVg/XdA7k8PRt2HuXOFW0lWXRRyGS/dxbWPVy+oI0NO0/H0S9fMJeOjuK7b2fOmMYEc/fyuJQRYDeSzuR2k1rnePblNH922YX07XuZ5QvmFrhVglyzjg741aGt/Pz5/YARifOxK98IUOCz96OjA1rbTs9w7OUEi8ntbGcVxr3ak52aJ8N1i2bnrovV7h6zpGQldxRX4llzI+zz19zc7Pt5lNE4LcBJU9FPAJYCXwYeA96HEZFzLfBIVH2IC2Gnm+VGUtQ7TneFny84CvwGbrdNSdY9tPu4m5vG+u78dW6wsrtrgAIlCoabpG9wFzctOtvX/+3n4lnfPcBjL7wMQDIB7dMm8lz/YYaGT4bOOWN38VSiCphThqHhk2Xn3ilmhNVaGmM/orTszwXuN/32SeCHSqmfikg3sEZEvgg8C3w3wj7EgrAPjN4V6I9ToYb1BQellIIVToXu3Npv/V9MKZYcsGUAAA8KSURBVLkpNvtiq12J3rHuBdTAkdyxboWs7ddsXCrpusDrPG8mC8/vfZ3n974OkFuwvXGZhL6WlVjcdjOC7O2Gzb0TtBh9XAyvKKNxfgNc4vL+TuDtUZ03jpTywIxmlEncsBSqU8lV0uoqp2BFsXsXpI2gz4xz8EglE0wcl3A97rpFs/m/W/cwMHQit8DrVHBe6RLg9ILtvLYpkYSiFqOYERT2dxb0XsbF8NK5cWqAOD0wccGp5CptdXkpgkpYekHaCPPMWIr8nsd7SWey/KT7dc5apwoigewuH6dczvPaN5HZSZsx/hDeZVIJ/AbSsL+zMANqHH6zWtnXCHF5YOKCW3ZIoGLpfqMsWBG0jWLPjN2yHho+mbdY69zM5JbN0m+gsVvuh46MsO63B0hnsq6phWvJlx3md1ZvRphW9ppRJ+rpvVt2SKispemnCPwUSlDZw7QRZO3gukWzmTA2ZbhybFE6liJ2+uyt7J3gPUDa++jWh7j4sv2oJyNMK3vNqDIaG128EmdV2tIspuycVEJ2NyXulrPeLzJl5MhhHnnhSMFCpnPwCtNfp1KMo1Vc7xsVddZLzagyGhkM3bJDRpk5M2i64ErI7mxjQ/cB1zbd5F3aaWSblLObWHjBVJZIS54Ctz63XpfbX2d7tUypKZ/jhLbsNaPKaISqeVmVpViaftae9Vn/4LFAs4YoFm+v6jyHvsFdBW16XYPb1ynufmyADORKGnrJW2p/42gh1+oaQyXRyl4zqozW9N7N1xrW7bK572hud6TTjWF3cYxLJRmXShZNbxxUdr++ubXhtSHJeQ2s5GNW/IxTqZWTm8d+jnJcVV7rD1Gv8fQPHgt0D+OMVvaaUafai15BFZIzt4pdMdotwZF0hiXSQtvUM8paeA3aNzf/eJDr6Uw+lkqQp9SCplQudg57G6u7dgdW1G7rEd37XsuFeEaxxuMctJdIS+idwHFBK/s6Jo7T6SgJ63aZP2MCG2y5VSzF6GYJVkpBROlOsLtlUskE119+YV7bi+ecLuA9LpUs2800LpU0FfXBnKK2ZHTbPeyU/Z4nXswbnKJwrzgH7bapZ9Ttb0Ur+zqlkdK7uhXOcDsmrNtl4cyJ3LmizTNCZVwqyZvOO5OzJ42vmCyl+MnDhHNet2g2P3t2N394yayiaQ2e6z8c2liwu5n6B4/xmDodCbW6azebdw7mxfPbn828wShBQQrkKNwrcUp3UC5a2dcpjbDgBIai8yqcYadSbhdnO937h0hnXndNLeDX5zA++WJtBR3U7btk7920i3ltU3IyWRuh7AXaLcs6rLFg33RlKXcrMsi5ccv+bNpltydxs8f9V/oZjmOIaKloZV+nNIrFYldQkF84w47zepSqOLysz6ADaik+eb/BIYyP3O1YSxm7bryyyRbG926Xw5l10mnZez2b89qmjJoSrvYa0mihlX2d0igWi93PDHj6mit1Pbysz6ADatgZV7HBwT74jEkmeGL7QTIeNV2dA57VB+t/+8Yrp2Xt9L2HtfIt/DJ+usnqLJtYrXWoelj/0sq+jmkEi2Vpp1FUPEjhjEpdD3s7YfOwVzrzojX4rO7azRPqoGdYpf3YtV3bWb5gLkCem8WeEnh99wALL5ia+67d9+43QBVTin73oJiso7EOVUra6riglb0m9lRzUAt77rAzjKAZMDfuOK3oAVLJwpqs1rGticFcJSi3vrgVRrEGhVQyQXPTWNe+lqsUi8ka9TpUOWmr44BOl6DR+LC+e4BbHtlW0e3zSzuDpxGwBoc//d1ZvsrTnh4hCXSe61+izq8vbnl1rls0O+fHv3fTLtfr4bYmEIZiskaZ8gK800NEfd7RQit7jcaDzX1HayJfSpDBwVKUS6SFMakkz+99veQ+uyk3e4pkrzw5i+e05GrgghEZFfb8frIGHfhKxUupR33e0UK7cTQaD/x20NYilsthJB3Mv+7Xjpt7J4g76dL2s3L+fa/IqHKI0mXn52Krh/Uvrew1Gg+8dtBWg6DRIJUKuXUqt6BrDSsXzCpY9A3T/2pTD0rdC63sNQ1HUMXjtoO2GvgtfDpliTLkNogidDt/vUSzxB2t7DUNRVjFU21Lb333AHese8HVneQlS7X77Dx/vUSzxJ3IlL2ItAHfB6YDWeDflFJfF5GpwMPA+cBLwPuVUq9G1Q+Nxk6cFI9dmVvYXSNxkaVRdnPXOlFG45wC/kYp1QksBD4hIp3AzcCjSqk5wKPma00ViSK8sFaJUxidXZkDTJkwhusWzc4p9NGWpdTnpJrRLI30bBcjMsteKbUf2G/+PSQiPcB5wNXAFeZh9wOPA5+Nqh8afxrNnxqnNBJ2ixjg8PFTuQRmUfvnnZT7nFTDtdRoz3YxRsVnLyLnA5cAXcB0cyAAOIDh5imgt7eXyZMn09PTMxpdHDWGh4drSqa1XYfyXAFru7bTmhgM3U6tyeVHawJWzE0Bg/T0eMtabZlaE3DTorO5/9lBXjp8Eii8R0FlsVOKXJV6TqLCTaZa73Mxwt6ngQH/2Uvkyl5EJgH/DvyVUup1kdM5tJVSWRHJun2vvb2doaEhOjo6ou7iqNLT01NTMi3PDrBh5+mt8csXzM1tpQ9DrclVCWpBpo4OaG3LT19Q6j2yKEWuSj0nUeEmU633uRhh71Nzs/+u6UiVvYiMxVD0Dyql/sN8e0BEzlVK7ReRc4GXo+yDxp84uTUalVq4R7XQh7DEsc9REmU0TgL4LtCjlPoX20drgWuB28z/H4mqD5pgVDtUT1OcWrhHtdCHsMSxz1ERpWV/KfAh4HkRec58728xlPwPReQjwG7g/RH2QaPRaDREG42zCUh4fHxlVOfVaDQaTSE666VGo9E0AFrZazQaTQOglb1Go9E0ALWYCC0FcODAAY4cOVI0djRuDAwM1J1MUJ9y1aNMUJ9yaZkMnWmScvu8FpX9uQDXXHNNtfuh0Wg0ceRc4EXnm7Wo7H8FLMbIq5MucqxGo9FoDFIYiv5Xbh8mslnXbAUajUajqSP0Aq1Go9E0ADXhxhGRJuBJYDxGn36slPq8iMwG1gBnAc8AH1JKjVSvp+Hwket7wOXAa+ahH1ZKPefeSm0iIingaWCvUurdcb9X4CrT94j/fXoJGMJwiZ5SSr017gWEPGT6AvBR4KB52N8qpX5elQ6WiIhMAb4DXIxR8Ok6QFGhe1Urlv0J4B1KqTcD84B3ishC4MvAV5VS7cCrwEeq2MdS8JIL4Eal1DzzX6wUiMmnAHv+1bjfKyiUCeJ/nwCWmP1/q/m6HgoIOWUC4/mz7lWsFL3J14FfKKUuAt6M8SxW7F7VhLJXSmWVUkfMl2PNf1ngHcCPzffvB95The6VjI9csUZEWoE/xLBCrKR3sb5XTpnqnKsx7hHE8F7VIyIyGbgMI3kkSqkRpdRhKnivakLZgzGFNhOmvQysxwgdOqyUOmUesgej0lWscMqllOoyP/qSiPxGRL4qIuOr2MVS+BpwE5AxX59F/O+VUyaLON8nMIyL/xKRZ0TkY+Z7gQoI1TBuMgF80rxX94rIG6rVuRKZjeGCuk9EnhWR74jIRCp4r2pG2Sul0kqpeUAr8Hbgoip3qSI45RKRi4HPYcj3NmAqMSrLKCLvBl5WSj1T7b5UCh+ZYnufbCxSSs0H/gCjDvRl9g+VUlniN9t0k+lu4EIMd+l+4J+r2L9SGAPMB+5WSl0CHMXhsin3XtWMsrcwpy6PAb8LTBERaxG5FdhbtY6ViU2udyql9psunhPAfRiDW1y4FFhuLpKtwXDffJ1436sCmUTkgZjfJwCUUnvN/18GfoIhw4BZOIg4FhByk0kpNWAaVhng28TvXu0B9thm/j/GUP4Vu1c1oexFpMVciUZEJgBLMRYnHgPeZx4Wu0InHnK9YLt5CQwf3Lbq9TIcSqnPKaValVLnAx8A/lspdQ0xvlceMn0wzvcJQEQmikiz9Tfw+xgyWAWEIGb3yksm616ZvJeY3Sul1AGgX07Xbb0S6KaC96omQi8xdn3db4a+JYEfKqV+KiLdwBoR+SLwLObiRYzwkuu/RaQFI9//c8D11exkhfgs8b5XbjwY8/s0HfiJqT/GAKuVUr8QkV8R3wJCXjL9QETmYbg5XgL+vHpdLJm/xHjmxgE7gVWYeqMS90rvoNVoNJoGoCbcOBqNRqOJFq3sNRqNpgHQyl6j0WgaAK3sNRqNpgHQyl6j0WgaAK3sNXWBiKRF5DkR2SYiPxKRM8po6woR+an593IR8Uw+JSJTROQvbK9niMiPvY7XaKpFrcTZazTlctxMS4GIPIgRE/8v1ofmxqiEucMyMEqptRgbW7yYAvwF8C3z+H2c3lym0dQMWtlr6pGNwO+IyPnAOqALeAvwLnOH4j9g1Bh4EVillDoiIu/ESIZ2DNhkNSQiHwbeqpT6pIhMB+4BLjA//jhwA3ChmexuPXAX8FOl1MVmPYO7gbcCp4BPK6UeM9tcDpyBkc/lJ0qpm8zNd981j88C9yqlvhrFBdI0HtqNo6krzPw8fwA8b741B/iWUuqNGMml/h64ykyk9TTwaVMpfxv4I4xB4RyP5u8EnjDrE8wHfouRrOpFM4f6jY7jPwFklVJvAlZg7KZuMj+bB/wJ8CbgT0SkzXzvPKXUxeZ37ivnWmg0drSy19QLE0zr+mmgj9PpGnYrpTabfy8EOoGnzGOvBWZhZLbcpZTaYWYWfMDjHO/AsNStbKaveRxnschqSyn1AsZ297nmZ48qpV5TSg1j5ECZhbFF/gIR+YY503g9uPgajT/ajaOpF3I+ewszf8pR21sJjJoCKxzH5X1vlDhh+zsNjFFKvSoibwaWYaw5vB+jNJ1GUzbastc0EpuBS0WkHXIZFOcCLwDni8iF5nErPL7/KIaf3ipKMxmjFmqzx/EbgWvM4+cCMzFqiroiImcDSaXUv2O4m+aHkE2j8UUre03DoJQ6CHwYeEhEfgP8ErjIdKV8DPiZiGzFO2f4p4AlIvI8RlH1TqXUKxhuoW0icrvj+G8BSfP4hzEKlp/Am/OAx00X0wMYxVM0moqgs15qNBpNA6Ate41Go2kAtLLXaDSaBkAre41Go2kAtLLXaDSaBkAre41Go2kAtLLXaDSaBkAre41Go2kAtLLXaDSaBuD/Acp2bPcNYmWXAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "predict = model(Variable(torch.from_numpy(X)))\n",
    "predict = predict.data.numpy()\n",
    "labels = np.asarray(labels)\n",
    "# Plot the results and compute the MAE\n",
    "plt.scatter(predict, labels, s=10)\n",
    "plt.grid()\n",
    "plt.xlabel('Predictions')\n",
    "plt.ylabel('True values')\n",
    "mae = np.mean(np.abs(predict - labels))\n",
    "plt.title('Mean Absolute Error : {:.2f} years'.format(mae))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "def train_mine(model, optimizer, epochs=1): \n",
    "    \"\"\"Train a model on CIFAR-10 using the PyTorch Module API.\n",
    "    Inputs:\n",
    "    - model: A PyTorch Module giving the model to train.\n",
    "    - optimizer: An Optimizer object we will use to train the model\n",
    "    - epochs: (Optional) A Python integer giving the number of epochs to train for\n",
    "    Returns: Nothing, but prints model accuracies during training.\n",
    "    \"\"\"\n",
    "    model = model.to(device=device) # move the model parameters to CPU/GPU for e in range(epochs):\n",
    "    # move to device, e.g. GPU\n",
    "    for t in range(500):\n",
    "        model.train() # put model to training mode\n",
    "        x = new_img_list\n",
    "        x = x.to(device=device, dtype=dtype) # move to y = y.to(device=device, dtype=torch.long)device, e.g. GPU\n",
    "        scores = model(x)\n",
    "        loss = F.cross_entropy(scores, np.asarray(Y[i]))\n",
    "\n",
    "        # Zero out all of the gradients for the variables which the optimizer # will update.\n",
    "        optimizer.zero_grad()\n",
    "        # This is the backwards pass: compute the gradient of the loss with # respect to each parameter of the model.\n",
    "        loss.backward()\n",
    "        # Actually update the parameters of the model using the gradients\n",
    "\n",
    "        # computed by the backwards pass.\n",
    "        optimizer.step()\n",
    "        if t % print_every == 0:\n",
    "            print('Iteration %d, loss = %.4f' % (t, loss.item())) \n",
    "    #         check_accuracy(loader_val, model)\n",
    "    #         print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def flatten(x):\n",
    "    N = x.shape[0] # read in N, C, H, W\n",
    "    return x.view(N, -1) \n",
    "\n",
    "class Flatten(nn.Module):\n",
    "    def forward(self, x):\n",
    "        return flatten(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "channel_1 = 64\n",
    "channel_2 = 32\n",
    "model = nn.Sequential(\n",
    "    nn.Conv2d(in_channels=256, out_channels=channel_1, kernel_size=5, stride=1, \\\n",
    "              padding=2, bias=True),\n",
    "    nn.Conv2d(in_channels=channel_1, out_channels=channel_2, kernel_size=5, stride=1, \\\n",
    "              padding=2, bias=True),\n",
    "    nn.ReLU(),\n",
    "    nn.MaxPool2d(kernel_size = 2, stride=2),\n",
    "    nn.Conv2d(in_channels=channel_2, out_channels=channel_1, kernel_size=3, stride=1, \\\n",
    "              padding=1,bias=True),\n",
    "    nn.Conv2d(in_channels=channel_1, out_channels=channel_2, kernel_size=3, stride=1, \\\n",
    "              padding=1,bias=True),\n",
    "    nn.ReLU(),\n",
    "    nn.MaxPool2d(kernel_size = 2, stride=2),\n",
    "    nn.Conv2d(in_channels=channel_2, out_channels=channel_1, kernel_size=3, stride=1, \\\n",
    "              padding=1,bias=True),\n",
    "    nn.Conv2d(in_channels=channel_1, out_channels=channel_2, kernel_size=3, stride=1, \\\n",
    "              padding=1,bias=True),\n",
    "    nn.ReLU(),\n",
    "#     nn.MaxPool2d(kernel_size = 2, stride=2),\n",
    "    Flatten(),\n",
    "    nn.Linear(channel_2 * 25 * 25, 1),\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "using device: cpu\n"
     ]
    }
   ],
   "source": [
    "USE_GPU = True\n",
    "\n",
    "dtype = torch.float32 # we will be using float throughout this tutorial\n",
    "\n",
    "if USE_GPU and torch.cuda.is_available():\n",
    "    device = torch.device('cuda')\n",
    "else:\n",
    "    device = torch.device('cpu')\n",
    "\n",
    "# Constant to control how frequently we print train loss\n",
    "print_every = 100\n",
    "\n",
    "print('using device:', device)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "resnet50 = models.resnet50()\n",
    "optimizer = optim.SGD(resnet50.parameters(), lr=1e-3,\n",
    "                     momentum=0.9, nesterov=True)\n",
    "\n",
    "train_mine(resnet50,optimizer)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}